一种监管公路临时占地利用与恢复的卫星遥感方法与流程

本发明涉及公路建设环境保护监管技术领域，特别涉及一种监管公路临时占地利用与恢复的卫星遥感方法。

背景技术：

公路临时占地包括取土场、弃土场、施工营地、拌合站、预制场等，按照环保法律法规的要求，公路建设完工后临时占地需要进行恢复，尤其是恢复地表的植被。由于公路临时占地的确切选址和具体利用方式在工程设计阶段很难敲定，往往是在施工阶段确定，因此建设前的环评报告只能在原则上要求临时占地的利用方式要环保且利用后要进行恢复，很难通过“建前环评——建后验收”的传统模式保证公路临时占地的恢复效果。尤其是在我国西部偏远地区甚至是无人区，生态环境十分脆弱，但人工监管成本又极高。卫星遥感技术“站的高、看的远”，具有大范围对地观测的能力，有利于降低监管成本。

技术实现要素：

本发明的目的是提供一种基于卫星遥感影像实现对监管公路临时占地利用与恢复状况进行监测的方法。

为此，本发明技术方案如下：

一种监管公路临时占地利用与恢复的卫星遥感方法，步骤如下：

s1、分别采集监管公路在施工前、施工期和竣工后的三幅多光谱卫星遥感影像，并通过几何校正和统一裁剪使三幅多光谱卫星遥感影像中同一地物在不同多光谱卫星遥感影像中的位置完全相同；

s2、在经过步骤s1处理后的三幅多光谱卫星遥感影像中任选一副作为参考影像，其余两幅作为调整影像，利用线性回归分析方法将两副调整影像的亮度调整为与参考影像一致，使三幅多光谱卫星遥感影像上未发生变化的同一地物的亮度一致；

s3、对经过步骤s2处理后的施工前的多光谱卫星遥感影像与施工期的多光谱卫星遥感影像之间的像素变化检测采用变化强度单峰直方图自动阈值分割的方法，并利用t-point法获得区分变化强度单峰直方图中的变化像素与未变化像素的最优阈值；然后剔除通过最优阈值判定的变化像素，对其余的像素重复步骤s2～s3直至在每个类型波段下两幅影像间线性回归方程的斜率差值△kn小于预设的终止计算阈值；将最后一次计算得到的在每个波段下调整影像与参考影像间线性回归方程的斜率kn和截距bn作为标准亮度值修正参数，重新对调整影像每个像素在每个波段的亮度值进行调整，并将其与参考影像之间的变化像素采用变化强度单峰直方图自动阈值分割方法和t-point法进行区分，识别变化像素；

s4、将经过步骤s3识别到的变化像素根据是否为相邻像素进行地块划分，并计算每一个地块在施工前、施工期和施工后的植被指数，通过不同时期的植被指数计算得到施工结束后该地块的植被的恢复率，进而对该地块的恢复情况是否符合要求进行判断。

进一步地，在步骤s1中，三幅多光谱卫星遥感影像为在相同季节拍摄得到。

进一步地，在步骤s2中，将两幅调整影像的亮度调整为与参考影像一致的具体步骤如下：

s201、提取参考影像中每个像素所包含的n个波段的亮度，使每个像素对应得到一个至少包含蓝光波段亮度b1x、绿光波段亮度b2x、红光波段亮度b3x和近红外波段亮度b4x的亮度值集合x＝{b1x,b2x,b3x,b4x,…,bnx}；

s202、提取每幅调整影像中每个像素所包含的n个波段的亮度，且调整影像每个像素提取的波段类型与步骤s201中对参考影像每个像素提取的波段类型一致，使调整影像每个像素同样对应得到一个至少包含蓝光波段亮度b1y、绿光波段亮度b2y、红光波段亮度b3y和近红外波段亮度b4y的亮度值集合y＝{b1y,b2y,b3y,b4y,…,bny}；

s203、以调整影像的第n个波段的亮度为x轴、参考影像的第n个波段的亮度为y轴建立平面直角坐标系，并将调整影像和参考影像上位于相同位置上的像素的第n个波段的亮度值分别作为x坐标和y坐标在直角坐标系中绘制若干个坐标点；通过对若干个坐标点进行线性回归分析，得到一元线性回归公式：bny＝kn×bnx+bn，求得在第n个波段下调整影像与参考影像间线性回归方程的斜率kn和截距bn，进而将调整影像每个像素的在第n个波段的亮度值bnx代入一元线性回归公式中，求得调整影像每个像素在该光波波段亮度的修正值bny；在该步骤中，n的取值依次为1,2,3,4,…，实现对调整影像的每个像素在每一个波段的亮度值进行修正，进而将调整影像的亮度修改为与参考影像的亮度一致。

进一步地，步骤s3的具体方法为：

s301、将经过步骤s2处理后的施工前的多光谱卫星遥感影像与施工期的多光谱卫星遥感影像的每两个位于相同位置上的像素在不同类型波段的亮度值代入至公式：变化强度图＝sqrt((bnx-b1y)²+(b2x-b2y)²+(b3x-b3y)²+(b4x-b4y)²+…+(bnx-bny)²)中，得到能够表现像素变化强度的数值图像；

s302、将经过步骤s301得到的数值图像变换为x轴表示变化强度、y轴表示像素个数的单峰直方图，进而通过将单峰直方图的波峰陡降部分和拖尾部分分别进行线性拟合，得到两条线性拟合的回归线残差和最小处对应的变化强度值，即用于判定变化像素的最优阈值；

s303、利用步骤s302得到的最优阈值剔除步骤s301中得到的数值图像中的变化像素，并对其余判定为未变化的像素重复上述步骤s203，计算得到的新的kn值；

s304、依次重复上述步骤s203、s301、s302和s303，直至在在每个类型波段下计算得到的新的kn值与前一次计算得到的kn值之间的差值小于预设的计算终止阈值；

s305、将最后一次计算得到的在每个波段下调整影像与参考影像间线性回归方程的斜率kn和截距bn作为标准亮度值修正参数，将调整影像每个像素在第n个波段的亮度值bnx代入一元线性回归公式：bny＝kn×bnx+bn中，求得调整影像每个像素在第n个波段下的标准亮度修正值bny，得到亮度修正影像；

s306、将经过步骤s305得到的亮度修正影像与参考影像重复上述步骤s301和s302，得到的准确的变化像素。

进一步地，步骤s4中恢复率的具体计算方法为：

s401、对每一个地块在施工前、施工期和施工后的多光谱卫星遥感影像中的植被指数进行计算，其中每个像素的植被指数计算公式为：ndvi＝(b4-b3)/(b4+b3)；对应地块的植被指数为ndvi的平均值；

s402、剔除步骤s401的计算结果中ndvi施工前＜ndvi施工期的地块；

s403、通过恢复率计算公式：恢复率％＝(ndvi竣工后-ndvi施工期)/(ndvi施工前-ndvi施工期)计算每个地块的恢复率，并根据恢复率的计算结果判断植被是否恢复到位。

进一步地，在步骤s401之前还包括步骤s400，根据实际临时施工地况，人为设定噪音地块像素个数阈值n，将包含像素的个数≤n的地块剔除。

与现有技术相比，该监管公路临时占地利用与恢复的卫星遥感方法通过对施工前、施工期和施工后三个阶段的三幅多光谱卫星遥感影像通过依次进行的图像校正和剪裁、亮度一致调整、变化像素判断、植被指数计算多个步骤以获得植被的恢复率，进而能够快速判断出公路临时占地利用与恢复，且判断结果准确，为临时施工监管单位的监管工作提供便利，有效降低监管管理成本，提高监管工作的效率。

附图说明

图1为本发明的监管公路临时占地利用与恢复的卫星遥感方法流程示意图；

图2为本发明的相对辐射纠正方法中的回归线受变化像素影响的示意图；

图3为本发明的t-point变化检测方法的原理示意图；

图4为本发明的实施例相对辐射纠正和变化检测不断迭代过程中的最优阈值变化趋势；

图5为本发明的实施例相对辐射纠正和变化检测不断迭代过程中的6个波段的线性回归斜率变化；

图6为本发明的实施例相对辐射纠正和变化检测不断迭代至稳定后的变化强度直方图及最优阈值；

图7(a)为本发明的具体实施例的施工前和施工后植被指数ndvi的变化图；

图7(b)为本发明的具体实施例的施工前的卫星遥感影像图；

图7(c)为本发明的具体实施例的施工期的卫星遥感影像图；

图7(d)为本发明的具体实施例的施工后的卫星遥感影像图。

具体实施方式

下面结合附图及具体实施例对本发明做进一步的说明，但下述实施例绝非对本发明有任何限制。

如图1所示，以依托公路工程为例对该监管公路临时占地利用与恢复的卫星遥感方法进行具体说明。

步骤一、卫星遥感影像挑选与预处理：

依托公路工程建设始于2007年下半年，竣工于2010年上半年。选择美国的landsattm卫星遥感影像对该工程的施工状况进行分析，其影像的分辨率为30m，在后续分析过程中，选取前六个波段的亮度值数据进行分析。

分别采集依托公路施工前、施工期和施工后的多光谱卫星遥感影像，三幅影像的拍摄时间分别是2007年9月3日、2009年8月23日和2010年8月10日，拍摄季节均为夏季，以保证农作物的物候一致。

对依托公路施工前、施工期和施工后的三幅多光谱卫星遥感影像通过几何校正和统一裁剪使三幅多光谱卫星遥感影像中同一地物在不同多光谱卫星遥感影像中的位置完全相同。

步骤二、在经过步骤s1处理后的三幅多光谱卫星遥感影像中任选一副作为参考影像，其余两幅作为调整影像，利用线性回归分析方法将两副调整影像的亮度调整为与参考影像一致，使三幅多光谱卫星遥感影像上未发生变化的同一地物的亮度初步调整为一致；

在该步骤中将两幅调整影像的亮度调整为与参考影像一致的具体步骤如下：

s201、提取参考影像中每个像素所包含的6个波段的亮度，使每个像素对应得到一个包含有六种光波波段的亮度值集合x＝{b1x,b2x,b3x,b4x,b5x,b6x}；

s202、提取每幅调整影像中每个像素所包含的6个波段的亮度，且调整影像每个像素提取的波段类型与步骤s201中对参考影像每个像素提取的波段类型一致，使调整影像每个像素同样对应得到一个包含有六种光波波段的亮度值集合y＝{b1y,b2y,b3y,b4y,b5y,b6y}；

s203、依次对每一幅调整影像进行如下的亮度值调整：

如图2所示，以调整影像的第1个波段，即蓝光波段的亮度为x轴、参考影像的第1个波段的亮度为y轴建立平面直角坐标系，并将调整影像和参考影像上位于相同位置上的像素的蓝光波段的亮度值分别作为x坐标和y坐标在直角坐标系中绘制若干个坐标点；接着，通过对若干个坐标点进行线性回归分析，得到一元线性回归公式：b1y＝k1×b1x+b1，求得在蓝光波段下调整影像与参考影像间线性回归方程的斜率k1和截距b1，进而将调整影像每个像素的在蓝光波段的亮度值b1x代入一元线性回归公式中，求得调整影像每个像素蓝光波段亮度的修正值b1y，实现对该幅调整影像在蓝光波段的亮度值的修正；

以此类推，依次对该幅调整影响的第2个波段、第3个波段、第4个波段、第5个波段和第6个波段进行如上述第1个波段的亮度调整，同样通过线性回归分析，得到第2～6个波段相应的斜率k2、k3、k4、k5和k6和截距b2、b3、b4、b5和b6；进而实现对每个像素在6个波段上的亮度值均进行调整，最终实现将调整影像的亮度修改为与参考影像的亮度一致。

步骤三、识别经过步骤s2处理后的施工前的多光谱卫星遥感影像与施工期的多光谱卫星遥感影像之间的像素变化；具体地，

s301、将经过步骤s2处理后的施工前的多光谱卫星遥感影像与施工期的多光谱卫星遥感影像的每两个位于相同位置上的像素在不同类型波段的亮度值代入至公式：变化强度图＝sqrt((bnx-b1y)²+(b2x-b2y)²+(b3x-b3y)²+(b4x-b4y)²+(b5x-b5y)²+(b6x-b6y)²)中，得到能够表现像素变化强度的数值图像；

s302、如图3所示将经过步骤s301得到的数值图像变换为x轴表示变化强度、y轴表示像素个数的单峰直方图；进而通过将该单峰直方图的波峰陡降部分和拖尾部分分别进行线性拟合，得到两条线性拟合的回归线残差和最小处对应的变化强度值，即用于判定变化像素的最优阈值；

s303、利用步骤s302得到的最优阈值剔除步骤s301中得到的数值图像中的变化像素，并对其余判定为未变化的像素重复上述步骤s203，计算得到的新的kn值；

s304、依次重复上述步骤s203、s301、s302和s303，直至在在每个类型波段下计算得到的新的kn值与前一次计算得到的kn值之间的差值小于预设的计算终止阈值0.01；

具体地，在重复步骤s203的过程中，斜率不断变化直至6个波段的斜率分别稳定在0.985、0.985、0.976、1.226、0.922、0.96；由于landsattm影像的第4波段是近红外波段，对植物长势很敏感，遥感影像中绝大多数像素是植物，因此该波段的斜率稳定较慢，在第10次重复计算后才趋于稳定，而其余的5个波段的斜率在第5次重复计算后趋于稳定，因此，如图5所示，经过11次重复计算后，得到最终稳定的kn值；而在s301和s302的过程中，如图4所示，通过t-point法变化检测到的最优阈值从21.77不断下降直至9.598，说明变化强度图的单峰直方图越来越窄，拖尾部分越来越往左缩进。实际上当重复计算至第五次时，最优阈值已经相对稳定；

s305、如图2所示，线i为亮度调整前作图得到的对若干个坐标点进行线性回归分析绘制的具有斜率k1和截距b1的的一条实际回归线，而线ii则为通过亮度值调整后按照同样方法获得的理想回归线；线i由于受到变化像素的影响，其相对于线i发生的偏移，因此在上述的步骤s304中，通过不断的剔除变化像素重新进行线性回归分析，一条理想回归线，进而获得区别变化像素与未变化像素的准确判断阈值；

因此，将最后一次计算得到的在每个波段下调整影像与参考影像间线性回归方程的斜率和截距作为标准亮度值修正参数，将调整影像每个像素在六个波段的亮度值分别带入对应的一元线性回归公式：bny＝kn×bnx+bn中，求得调整影像每个像素在六个波段下的标准亮度修正值，进而得到最终的亮度修正影像；

s306、如图5所示，将经过步骤s305得到的亮度修正影像与参考影像重复上述步骤s301和s302，得到的准确的变化像素；经过11次重复计算后，变化强度图的最优阈值是23，以变化强度值23为阈值，对全部像素进行分割，变化强度值小于23的像素判定为未变化像素，而变化强度值大于23的像素判定为变化像素。

步骤四、将经过步骤s3识别到的变化像素根据是否为相邻像素进行地块划分，识别到的临时占地30余处，人工剔除掉公路路基这一长条形的变化地块以及周边其它建设造成的变化，并通过现场核实其余未能识别的临时占地有的是租用村庄内现有大空地(未破坏植被)，还有的是直接设置在河滩地上(砾石为主)，最终筛出公路临时占地13处；

对每一个地块在施工前、施工期和施工后的多光谱卫星遥感影像中的植被指数进行计算，其中每个像素的植被指数计算公式为：ndvi＝(b4-b3)/(b4+b3)；对应地块的植被指数为ndvi的平均值；然后剔除步骤s401的计算结果中ndvi施工前＜ndvi施工期的地块；最终通过恢复率计算公式：恢复率％＝(ndvi竣工后-ndvi施工期)/(ndvi施工前-ndvi施工期)计算每个地块的恢复率，并根据恢复率的计算结果判断植被是否恢复到位；

表1：

由表1所示，该13处临时占地的ndvi随着公路施工先降低后增加；其中3处的恢复率不足15％，可以认为基本没有恢复，ndvi的少量恢复来自于杂草自然生长；其中5处的恢复率在20～60％之间，有一定程度的恢复；其中5处的恢复率在70％以上，可以认为恢复到位，甚至有两处超过100％。