一种基于多源时空数据的地质灾害灾毁耕地提取方法与流程

本发明属于自然资源管理技术领域，具体是一种基于多源时空数据的地质灾害灾毁耕地提取方法。

背景技术：

目前，地质灾害灾毁耕地的识别提取主要有外业实地调查方法和利用灾后遥感影像进行人工目视交互解译方法，这两种方法具有各自的不足：外业实地调查方法需要投入大量人力物力财力，且因灾后地形地貌改变，可能难以获取灾毁耕地的准确面积；利用灾后遥感影像进行人工目视交互解译的方法效率虽比外业调查高，但仍存在耗费时间长的问题，同时遥感解译结果的正确性依赖于解译人员的经验，存在一定主观性。

技术实现要素：

本发明的目的在于克服现有技术的缺陷，提供一种基于多源时空数据的地质灾害灾毁耕地提取方法。

为实现上述目的，本发明采用了如下技术方案：

一种基于多源时空数据的地质灾害灾毁耕地提取方法，包括具体步骤如下：

s1、筛选ndvi变化目标：以灾害灾毁前后的光学遥感影像为数据源，通过灾害灾毁前后差值指数变化为阈值筛选ndvi变化目标；

s2、确定时序非稳定目标：以灾毁发生时序上的sar影像为基础，通过计算相干系数方差得到时序非稳定目标；

s3、确定地质灾害灾毁耕地信息：将步骤s1中得到的ndvi变化目标、步骤s2中时序非稳定目标与土地利用数据进行叠加分析，得到地质灾害灾毁耕地信息。

优选的，所述步骤s1中筛选ndvi变化目标的具体步骤如下：

以灾害灾毁前后的光学遥感影像为数据源，经过相对辐射校正消除不同时相影像之间的辐射差异；

通过正射纠正得到带有统一地理坐标的影像；

将灾害灾毁前后监测区内的影像分别镶嵌，得到灾害灾毁前后的完整影像；

分别计算灾害灾毁前后的ndvi影像，两期影像相减，以灾害灾毁前后差值指数变化为阈值筛选ndvi变化目标。

优选的，所述步骤s2中确定时序非稳定目标的具体步骤如下：

以灾害灾毁发生时序上的sar影像为基础，对sar影像进行聚焦、滤波、配准等预处理；

将相邻sar影像进行两两干涉，计算相干系数并进行地理编码；

以时序方差为阈值筛选初始相干系数变化目标；

通过小区域去除和数学形态学滤波，得到时序非稳定目标。

优选的，所述ndvi为归一化植被指数。

优选的，所述光学遥感影像是获取可见光到红外线这部分电磁波信息进行成像，得到的影像。

优选的，所述光学遥感影像是通过航天或航空成像方式，拍摄的地表的影像，是获取地面目标的反射、散射或辐射的电磁波信息。

本发明一种基于多源时空数据的地质灾害灾毁耕地提取方法，本发明采用时序sar数据、多时相遥感影像、数字高程模型dem和土地利用现状数据，建立了基于多源时空数据的地质灾害灾毁耕地自动提取方法，该方法在对地质灾害灾毁耕地提取时人工参与少、自动化程度高，提取结果准确性高，且可以快速准确的提取地质灾害灾毁耕地信息，为自然资源管理工作提取技术支撑。

附图说明

图1是本发明一种基于多源时空数据的地质灾害灾毁耕地提取方法的整体流程图。

具体实施方式

以下结合附图1，进一步说明本发明一种基于多源时空数据的地质灾害灾毁耕地提取方法的具体实施方式。本发明一种基于多源时空数据的地质灾害灾毁耕地提取方法不限于以下实施例的描述。

一种基于多源时空数据的地质灾害灾毁耕地提取方法，如图1所示，具体步骤如下：

s1、筛选ndvi变化目标：以灾害灾毁前后的光学遥感影像为数据源，通过灾害灾毁前后差值指数变化为阈值筛选ndvi变化目标；

s2、确定时序非稳定目标：以灾毁发生时序上的sar影像为基础，通过计算得到时序非稳定目标；

s3、确定地质灾害灾毁耕地信息：将步骤s1中得到的ndvi变化目标、步骤s2中时序非稳定目标与土地利用数据进行叠加分析，得到地质灾害灾毁耕地信息。

所述步骤s1中筛选ndvi变化目标的具体步骤如下：

以灾害灾毁前后的光学遥感影像为数据源，经过相对辐射校正消除不同时相影像之间的辐射差异；

通过正射纠正得到带有统一地理坐标的影像；

将灾害灾毁前后监测区内的影像分别镶嵌，得到灾害灾毁前后的完整影像；

分别计算灾害灾毁前后的ndvi影像，两期影像相减，以灾害灾毁前后差值指数变化为阈值筛选ndvi变化目标。

所述步骤s2中确定时序非稳定目标的具体步骤如下：

以灾害灾毁发生时序上的sar影像为基础，对sar影像进行聚焦、滤波、配准等预处理；

将相邻sar影像进行两两干涉，计算相干系数并进行地理编码；

以时序方差为阈值筛选初始相干系数变化目标；

通过小区域去除和数学形态学滤波，得到时序非稳定目标。

所述ndvi为归一化植被指数。

所述光学遥感影像是获取可见光到红外线这部分电磁波信息进行成像，得到的影像；所述光学遥感影像是通过航天或航空成像方式，拍摄的地表的影像，是获取地面目标的反射、散射或辐射的电磁波信息。

实施例1

参照图1，本实施例公开了一种基于多源时空数据的地质灾害灾毁耕地提取方法中灾害前后光学遥感影像的处理方法，包括相对辐射校正、正射纠正、影像镶嵌、ndvi指数计算、ndvi变化目标提取，其中：

相对辐射校正

相对辐射校正是以一幅清晰的影像为参考影像，利用其它时相影像(目标影像)和参考影像对应像素的灰度值之间的关系，建立各波段之间的校正模型，并将目标影像逐波段的归一化到参考影像，使目标影像和参考影像具有相同的辐射尺度，相对辐射校正能纠正因大气状况变化造成的差异，能减小传感器状态产生的噪声，

以不同时相的两景影像为例，波段数为n，向量x＝[x1…xn]，y＝[y1…yn]分别表示两幅影像的像素值，对x、y分别做线性变换u＝a^tx，v＝b^ty，其中a＝[a1…an]^t，b＝[b1…bn]^t，则两景影像之间的变化信息可以集中到差值影像d＝u-v中，不同的系数a和b得到不同的差值影像d，d的像素值表示变化的大小，要最大限度的突出变化，则d的方差应最大，如果对a和b乘上一个常数c，则方差会变c²倍，因此对a和b作出限制，令u和v具有单位方差，即在式(1)和式(2)的约束下，满足式(3)，

var{a^tx}＝var{b^ty}＝1(1)

cov(u,v)>0(2)

var{d}＝var{a^tx}+var{b^ty}-2cov{a^tx,b^ty}＝2(1-corr{a^tx,b^ty})→max(3)

根据典型相关理论，ρ＝corr{u,v}是典型相关系数，假设向量具有零均值，则z的协方差矩阵可表示为其中var(u)＝a^tσ11a，var(u)＝b^tσ22b和cov(u,v)＝a^tσ12b，利用拉格朗日乘数法解算cov(u,v)极值，即可求解满足约束条件的极值a和b，此时利用a和b即可获得典型变量u和v，则mad变量及其方差为

式中i＝1…n，

令z表示任一像素的所有变化信息，对于未变化像素，mad变量应近似为0，则z近似服从n维自由度的卡方分布，由此可通过卡方分布检验获取未变化像素点(设置置信度大于0.95提取未变化像素点)，

根据目标影像和参考影像未变化像素点的值建立线性回归模型，即可实现参考影像的相对辐射处理，本发明采用正交回归分析建立回归方程，假设在二维平面有n个数据点pi(xi,yi),i＝1…,n，线性正交回归通过寻找直线l：

y＝α+βx(6)

使得目标函数最小，其中d(pi,l)是点pi到l的几何距离，线性正交回归较最小二乘法而言，建立的回归方程考虑了x和y方向误差，使得解算结果更加可靠和稳定。

正射纠正

以监测区域内含坐标信息的遥感影像和dem数据为基础，实现相对辐射校正后影像的倾斜改正和投影差改正，首先在基础底图和待纠正影像上选择均有的同名明显特征地物点，如道路交叉处、球场角、围墙角等位置，且全图均匀分布，然后根据纠正过程中控制点残差文件，检查正射纠正控制点点位精度，平地和丘陵地残差中误差1个像素以内，山地和高山地残差中误差2个像素以内，最后选择二次多项式重采样方法对影像进行正射纠正，并对结果进行目视对比检查。

影像镶嵌

以灾毁发生前后为时间节点，分别对影像进行镶嵌，首先，检查影像重叠区域内影像几何精度，点位误差不大于2个像素，然后在相邻影像重叠区勾绘镶嵌线，镶嵌线尽量沿着道路、河流、林地等地物走向勾绘，避免穿越房屋，最后按照镶嵌线裁切多余数据，并按灾毁前、灾毁后分别镶嵌监测区完整影像，灾毁前影像(记为i1)，灾毁后影像(记为i2)。

ndvi指数计算

由于植物叶片组织对蓝波段和红波段辐射有强烈吸收，但对近红外波段辐射反射强烈，而且近红外波段反射会随着植被覆盖的增加而增强，归一化植被指数(normalizeddifferencevegetationindex，ndvi)通过数学变换增强红波段和近红外波段的差异，能够有效反映地面植被覆盖情况，可作为植被覆盖变化检测的指标，其定义如下：

ndvi的值域范围为[－1，1]，负值表示地面覆盖为云、水、雪等，对可见光高反射；0表示地表覆盖有岩石或裸土等；正值越大，表示植被覆盖度越好，ndvi变化目标提取

首先根据上一步骤公式7分别计算灾毁前后监测区域内的ndvi值，然后采用差值指数评价是否存在地表变化，差值指数计算见式(8)，最后以变化值小于－0.05为阈值提取变化目标，

δndvi＝ndvii2-ndvii1(8)

δndvi的值域范围为[－2，2]，ndvii2表示灾毁后影像归一化植被指数，ndvii1表示灾毁前影像归一化植被指数，通常，当ndvi值大于0.1时，其代表的像素才被认为是植被像素，故这里仅选择ndvii1取值大于0.1的区域参与计算。

实施例2

参照图1，本实施例公开了一种基于多源时空数据的地质灾害灾毁耕地提取方法中时序sar影像数据处理，包括数据预处理、相干系数计算、地理编码、初始变化目标提取和小区域和数学形态学滤波，其中：

数据预处理

本方案涉及sar影像的预处理主要包括影像配准、影像多视、影像滤波等，首先将多期相同传感器影像按时间序列进行排列，选择日期最早或日期居中的影像作为主影像，并将其余影像与主影像进行配准，然后通过增强lee滤波、kuan滤波、自适应滤波等方法，消除相干电磁波所形成的斑点噪声影响，最后根据距离向分辨率、方位向分辨率和中心入射角设置sar影像的多视比，使地距方位向和距离向分辨率一致，

相干系数计算

sar影像相干系数是同时利用sar影像对的幅度和相位信息形成的，能够直观反映两次sar成像过程中目标稳定性的估计，其计算公式可见式(9)，

对于两景sar复数影像s1和s2，其相干系数的定义公式如下：

式中，*表示共轭复数；s1和s2为经过预滤波的单视复数sar影像，相干系数γ的值域为[0，1]，0为失相干，1为完全相干，

高相干性表示稳定的地表硬目标，低相干性则说明两次成像期间地表发生了变化，根据相干系数的变化这一重要特征，可进行地物变化情况评定。

地理编码

sar影像的坐标系统为斜距/多普勒坐标系，为了消除sar产品的扭曲和畸变，并方便与光学影像、灾毁地矢量等常规地学产品联合应用及对照分析，就必须将sar相干系数结果转换到地图投影坐标系下，地理编码是借助轨道状态矢量、多普勒方程、椭球面方程和斜距方程来计算sar影像像素对应的大地坐标，再以dem数据投影参数为准，将其投影至地理参考系中，

初始变化目标提取

利用灾毁前后两景sar影像计算相干系数，由于地表地物发生变化，必将导致相干系数与未发生灾毁之前、未发生灾毁地区的相干系数有所不同，根据这一特性，我们提出基于时序sar数据的变化目标提取方法，若在监测时间段内地表地物发生了灾毁，则相干系数在时序上就会突然降低或升高，以时序上相干系数的方差为阈值，可实现自动提取地表突变的目标点，

假设有时序sar影像拍摄时间t1,…,tn，基于小基线集技术干涉理论，设置时间基线和垂直基线阈值，可生成多个干涉对每个干涉对即可生成1个相干系数图为消除不同时间跨度造成的相干系数差异，在时序相干系数变化探测之前，我们基于时间长度对相干系数进行归一化处理，如式(10)，

计算归一化后的时序相干系数方差，并设置方差阈值，即可探测到时序上发生突变的变化目标点。

小区域和数学形态学滤波

依据方差阈值自动提取的变化目标在影像上存在小面积噪声区域和边界不连续的情况，为得到平滑的变化目标范围，先对小面积斑块进行去除，保留下大面积变化区域，然后利用数学形态学方法进行滤波，滤波方法为开运算，即对影像先腐蚀后膨胀，实现平滑边界，消除细小的尖刺，断开窄小的连接，使变化目标更加完整，最后根据变化目标范围情况，重复使用上述滤波方法，最终得到完整平滑的时序非稳定目标范围。

实施例3

参照图1，本实施例公开了一种基于多源时空数据的地质灾害灾毁耕地提取方法中地质灾害灾毁地耕地信息提取，具体方法如下：

将灾害前后光学遥感影像处理得到的ndvi变化目标、时序sar影像数据处理得到的时序非稳定目标和当年土地利用数据中耕地数据均统一到相同坐标系下；

以时序非稳定目标为源数据，通过叠加分析提取与耕地数据有重叠的范围，并存储为耕地非稳定目标；

以耕地非稳定目标为源数据，叠加ndvi变化目标数据，有重叠范围的斑块即为地质灾害灾毁耕地范围。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。