一种联合灾前光学与灾后PolSAR图像的滑坡检测方法及系统

一种联合灾前光学与灾后polsar图像的滑坡检测方法及系统
技术领域
1.本发明涉及滑坡检测技术领域，尤其涉及一种联合灾前光学与灾后polsar图像的滑坡检测方法及系统。


背景技术：

2.由于暴雨事件是滑坡最常见的触发因素之一，在滑坡发生后出现云层覆盖的可能性很高。在大多数情况下，灾害发生后，研究区域的polsar(极化合成孔径雷达)数据比无云光学图像更早可用，基于polsar图像的滑坡区域检测方法能够更快更准确的响应灾难。但是仅通过灾后polsar图像容易出现虚警较多的情况，而在滑坡前不久采集的polsar图像往往是不可用的。当前光学遥感卫星发展较为成熟，可以为滑坡检测提供灾前信息，因此基于灾前光学和灾后polsar图像的滑坡区域检测方法可以为滑坡灾后的快速响应提供支持。
3.ainsworth等[ainsworth t,kelly j,lee j s.classification comparisons between dual-pol,compact polarimetric and quad-pol sar imagery[j].isprs journal of photogrammetry and remote sensing,2009,64:464-471.]认为包含至少两种极化方式的polsar数据可以提供更多的地面信息，从而能够更好地进行土地覆盖分类和滑坡区域检测，全极化(即包含水平(h)和垂直(v)极化波的所有四种组合)数据则可以为滑坡区域检测提供最完备的检测条件。在滑坡区域检测任务中，仅利用滑坡后采集到的单时相polsar图像，会因为缺少灾前信息，容易出现较多的虚警；同时利用滑坡前后的双时相polsar图像，理论上基于更多的极化信息可以降低滑坡区域检测的虚警，但基于灾前和灾后polsar图像的滑坡区域检测方法，要求采集到的两幅图像具有相同的几何结构，然而，在大多数情况下，在滑坡前不久采集的polsar图像，尤其与下一个可能事件后所采集图像具有相同几何结构polsar图像是不可用的。由于当前polsar任务，如terrasar-x、cosmo-skymed、radarsa-2和高分三号，并不能系统地覆盖整个世界，每次采集都必须手动编程。此外，由于卫星上的磁盘空间有限，尤其是由于下行传输速率有限，这些卫星无法在短时间内提供全球覆盖，因此，通常在滑坡灾害发生前不久，并没有可用的记录存档图像。为了解决这一问题，plank等[plank s,twele a,martinis s.landslide mapping in vegetated areas using change detection based on optical and polarimetric sar data[j].remote sensing,2016,8(4).]提出了一种基于灾前光学遥感图像和灾后polsar图像的半自动滑坡区域检测方法。在实际应用中，光学图像通过开放平台免费获取，polsar图像可以在滑坡发生后对卫星进行人工编程采集得到，避免因图像几何结构限制导致灾前polsar图像不可用的问题。
[0004]
plank等提供的方法中，在使用遥感图像进行滑坡检测的基础上，还需要通过数字高程模型dem进行滑坡筛选，增加了数据获取的难度。因为在实际应用中，合适的dem数据通常难以获取。同时，该方法作为少数公开的联合灾前光学与灾后polsar图像的滑坡区域检测方法之一，其精度还有待进一步提高。


技术实现要素：

[0005]
本发明针对plank等提供的方法中存在的问题，提出一种联合灾前光学与灾后polsar图像的滑坡检测方法及系统，该方法首先基于灾后polsar图像的极化特征参数，通过自动阈值法确定分割阈值，提取疑似滑坡区域，同时对灾前光学遥感图像进行面向对象分类以得出容易造成检测虚警的干扰区域，进一步通过逻辑运算进行虚警剔除，得出滑坡区域检测结果。该方法在不依赖其他类型数据辅助的条件下，通过灾前光学和灾后polsar图像可以以较高的准确度完成滑坡区域检测。
[0006]
为了实现上述目的，本发明采用以下技术方案：
[0007]
本发明一方面公开一种联合灾前光学与灾后polsar图像的滑坡检测方法，包括：
[0008]
基于polsar提取疑似滑坡区域：对灾后polsar图像进行预处理，然后计算出滑坡区域检测过程中所需的极化特征参数，通过自动阈值确定方法确定分割阈值，进而对其进行形态学处理得到疑似滑坡区域；
[0009]
基于光学遥感图像提取干扰区域：以灾后polsar图像为基准对光学遥感图像进行几何校正，并进行多尺度分割，最后通过设计的分类规则筛选出容易对滑坡检测造成虚警的干扰区域；
[0010]
在疑似滑坡区域的基础上，结合光学遥感图像得到的干扰区域进行虚警剔除，确定最终滑坡区域。
[0011]
进一步地，所述计算出滑坡区域检测过程中所需的极化特征参数包括：
[0012]
对预处理后灾后polsar图像的每一个像素计算香农熵。
[0013]
进一步地，所述基于polsar提取疑似滑坡区域还包括：
[0014]
通过yamaguchi分解得到的二次散射来剔除建筑物带来的虚警。
[0015]
进一步地，按照如下方式确定分割阈值t：
[0016]
t＝αt1+(1-α)t2[0017]
其中α为阈值的权重，t1、t2分别代表滤波与未滤波图像得到的阈值，式中
[0018]
t
x
＝argmax
0＜k＜l-1
(σ
*2
)
[0019]
其中x取值为1或2，σ
*2
为最大类间方差
[0020]
σ
*2
＝p0(μ
0-u)2+p1(μ
1-u)2[0021][0022][0023]
其中p0、p1分别表示前景和背景对应的概率；μ1、μ2分别表示前景和背景的均值；pi表示第i个灰度级出现的概率；ni表示第i个灰度级所包含的像素数量。
[0024]
进一步地，所述基于光学遥感图像提取干扰区域具体包括：
[0025]
以polsar图像为基准，通过地面控制点法对光学遥感图像进行几何校正。
[0026]
进一步地，所述基于光学遥感图像提取干扰区域具体包括：
[0027]
基于ecognition软件对几何校正后的光学遥感图像进行多尺度分割。
[0028]
进一步地，在所述多尺度分割过程中，通过esp2工具进行分割尺度的优选。
[0029]
进一步地，所述分类规则包括：
[0030]
判断不同地物类型的差值植被指数是否小于1500，若是，则将其分类为耕地及裸土区域，若否，则结合绿色波段，判断不同地物类型的差值植被指数是否位于1500至4000之间，若是，则将其分类为稀疏植被覆盖区。
[0031]
本发明另一方面提出一种联合灾前光学与灾后polsar图像的滑坡检测系统，包括：
[0032]
疑似滑坡区域提取单元，用于基于polsar提取疑似滑坡区域：对灾后polsar图像进行预处理，然后计算出滑坡区域检测过程中所需的极化特征参数，通过自动阈值确定方法确定分割阈值，进而对其进行形态学处理得到疑似滑坡区域；
[0033]
干扰区域提取单元，用于基于光学遥感图像提取干扰区域：以灾后polsar图像为基准对光学遥感图像进行几何校正，并进行多尺度分割，最后通过设计的分类规则筛选出容易对滑坡检测造成虚警的干扰区域；
[0034]
最终滑坡区域确定单元，用于在疑似滑坡区域的基础上，结合光学遥感图像得到的干扰区域进行虚警剔除，确定最终滑坡区域。
[0035]
进一步地，所述计算出滑坡区域检测过程中所需的极化特征参数包括：
[0036]
对预处理后灾后polsar图像的每一个像素计算香农熵。
[0037]
进一步地，所述基于polsar提取疑似滑坡区域还包括：
[0038]
通过yamaguchi分解得到的二次散射来剔除建筑物带来的虚警。
[0039]
进一步地，按照如下方式确定分割阈值t：
[0040]
t＝αt1+(1-α)t2[0041]
其中α为阈值的权重，t1、t2分别代表滤波与未滤波图像得到的阈值，式中
[0042]
t
x
＝argmax
0＜k＜l-1
(σ
*2
)
[0043]
其中x取值为1或2，σ
*2
为最大类间方差
[0044]
σ
*2
＝p0(μ
0-u)2+p1(μ
1-u)2[0045][0046][0047]
其中p0、p1分别表示前景和背景对应的概率；μ1、μ2分别表示前景和背景的均值；pi表示第i个灰度级出现的概率；ni表示第i个灰度级所包含的像素数量。
[0048]
进一步地，所述基于光学遥感图像提取干扰区域具体包括：
[0049]
以polsar图像为基准，通过地面控制点法对光学遥感图像进行几何校正。
[0050]
进一步地，所述基于光学遥感图像提取干扰区域具体包括：
[0051]
基于ecognition软件对几何校正后的光学遥感图像进行多尺度分割。
[0052]
进一步地，在所述多尺度分割过程中，通过esp2工具进行分割尺度的优选。
[0053]
进一步地，所述分类规则包括：
[0054]
判断不同地物类型的差值植被指数是否小于1500，若是，则将其分类为耕地及裸土区域，若否，则结合绿色波段，判断不同地物类型的差值植被指数是否位于1500至4000之间，若是，则将其分类为稀疏植被覆盖区。
[0055]
与现有技术相比，本发明具有的有益效果：
[0056]
针对当前基于灾前光学和灾后polsar图像的滑坡检测方法需要地学数据辅助且检测精度有待提升的问题，本发明提出了一种仅利用灾前光学遥感图像和灾后polsar图像的快速滑坡检测方法及系统。该方法通过自动阈值方法确定滑坡分割阈值，进而分割得到疑似滑坡区域，可以在避免polsar图像噪声影响的情况下得到较为完整的滑坡区域；然后通过面向对象分类的方法筛选出容易对滑坡检测造成干扰的区域；最后通过逻辑运算结合光学和polsar图像的阶段性结果，得到滑坡区域。在数据集的实验中，该方法总体精度为87.80％，滑坡生产者精度为71.66％，说明了本发明方法在无需地学数据辅助的情况下，可以有效检测出滑坡区域。
附图说明
[0057]
图1为本发明实施例一种联合灾前光学与灾后polsar图像的滑坡检测方法的流程图；
[0058]
图2为本发明实施例对象层次结构图；
[0059]
图3为本发明实施例灾前光学图像分类规则；
[0060]
图4为本发明实施例数据的地理位置；
[0061]
图5为本发明实施例数据集的基本情况图；
[0062]
图6为本发明实施例不同方法的检测结果柱状图；
[0063]
图7为本发明实施例不同方法的检测结果对比；
[0064]
图8为本发明实施例一种联合灾前光学与灾后polsar图像的滑坡检测系统框图；
[0065]
图9为本发明实施例系统主界面图；
[0066]
图10为本发明实施例系统运行结果示例。
具体实施方式
[0067]
下面结合附图和具体的实施例对本发明做进一步的解释说明：
[0068]
一种联合灾前光学与灾后polsar图像的滑坡检测方法(landslide detection method combining optical image and polsar image，ldop)，包括：
[0069]
1ldop检测方法流程
[0070]
本发明方法的流程框架如图1所示，在polsar数据处理部分，首先对灾后polsar数据进行预处理，然后计算出滑坡区域检测过程中所需的极化特征参数，通过自动阈值确定方法确定分割阈值，进而对其进行形态学处理得到疑似滑坡区域。在光学数据处理部分，首先以灾后polsar图像为基准对光学图像进行几何校正，并进行多尺度分割，最后通过设计的分类规则筛选出容易对滑坡检测造成虚警的干扰区域。最后在疑似滑坡区域的基础上，结合光学遥感图像得到的干扰区域进行虚警剔除，确定最终滑坡区域。
[0071]
本发明采用的形态学处理方法为先开运算后闭运算。由于地形校正后也无法完全消除polsar侧视成像造成的图像固有几何失真，同时由于光学遥感传感器平台俯仰角，传感器摇摆扫描，地形起伏等因素，光学遥感图像也存在一定程度的失真。为了使光学遥感图像和polsar图像更加吻合，本节以polsar图像为基准，通过地面控制点法(ground control-point,gcp)对光学遥感图像进行几何校正。
[0072]
2基于polsar提取疑似滑坡区域
[0073]
2.1极化特征参数计算
[0074]
从信息论的角度来说，确定事件(即发生的概率接近于1)几乎不包含有效信息；而对于不确定事件，其发生概率越接近于0，包含的信息量就越多。因此，克劳德
·
艾尔伍德
·
香农提出信息熵(香农熵)的概念：将离散随机变量所有可能状态的自信息求数学期望，得到的平均信息量被称为香农熵。
[0075]
假设一个离散型随机变量及其概率分布为：
[0076][0077]
则香农熵可表示为：
[0078][0079]
其中pi是概率密度函数，其取值范围为[0,1]，且
[0080]
香农熵的提出补充了信息量在量化度量上的缺失，在遥感的应用上，香农熵被用于衡量遥感图像所包含信息量。在polsar图像的相关研究中，香农熵可以反应地物目标散射的随机性，某一区域的香农熵值越大，表明该区域包含的信息量越多，散射特性越复杂，如森林区域；其值越低，表明所包含的信息量越低，散射特性越单一，如建筑区域。
[0081]
polsar图像的每一个像素都可以用均值为0且服从三维圆高斯分布的散射梯度矢量lj表示：
[0082][0083]
其中，det(
·
)表示对矩阵取行列式，mj表示极化协方差梯度矩阵，则散射特征imj可表示为：
[0084]
imj＝tr(mj)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0085]
极化特征可表示为：
[0086][0087]
在lj的概率密度函数已知的情况下，香农熵可定义为：
[0088]
s[p(lj)]＝∫p(lj)loge[p(lj)]dljꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)
[0089]
其中∫dlj表示对lj求积分，对于圆对称高斯过程，香农熵由散射分量和极化分量两部分构成：
[0090]
[0091][0092]
因此，
[0093][0094]
根据香农熵的定义，滑坡区域常常无植被覆盖，由表面散射为主的裸土构成，散射特性较为单一，常常呈现较低的熵值，而如果滑坡背景存在的建筑物主要由单一的二次散射构成，无植被覆盖，同样呈现较低的熵值，仅通过香农熵并不能很好的避免建筑物的干扰，因此本发明引入yamaguchi分解得到的二次散射来剔除建筑物带来的虚警。
[0095]
yamaguchi分解是一种基于散射模型的非相干分解算法，可以很好的匹配地物的基本散射机制。该分解方法从极化相干矩阵t中提取出各个散射分量对应的功率大小，即表面散射功率ps、二次散射功率pd、体散射功率pv和螺旋体散射功率ph，极化相干矩阵t可以表示为各散射功率的加权组合形式：
[0096]
t＝fsts+fdtd+fvtv+fhthꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(10)
[0097]
其中，ts、td、tv和th分别为表面散射、二次散射、体散射和螺旋体散射分量对应的极化相干矩阵；fs、fd、fv和fh分别为ps、pd、pv和ph的分解系数。
[0098]
在yamaguchi分解产生的结果中，ps主要表征了水面、草地和裸土等表面散射类型为主的地物散射功率大小；pd主要表征了两种不同电介材料(如粗糙墙壁和平整路面)所形成垂直结构的二次散射功率大小；pv主要表征了植被冠层不规则结构产生的体散射功率大小；ph主要表征了散射不对称的复杂目标或结构产生的螺旋体散射功率大小。
[0099]
2.2疑似滑坡区域分割阈值的确定
[0100]
根据式(9)可知香农熵sej仅与极化协方差梯度矩阵mj有关，通过各极化通道的散射矢量计算可以得到mj的行列式，因此，香农熵可以反应丰富的极化信息。根据不同地物的散射特性，植被区域的散射特性比较复杂，滑坡区域的散射特性比较单一，因此森林区域的香农熵要大于滑坡区域，根据植被区域和滑坡区域香农熵的差异，选取合适的阈值，对所有像素的香农熵进行判定，大于该阈值的像素区域被判定为森林区域，小于该阈值的像素区域被判定为滑坡初筛区域。
[0101]
对于滑坡周围存在较多建筑物的区域，仅通过香农熵并不能完全区分建筑物和滑坡，因此，根据滑坡和建筑物二次散射的差异，再次选择阈值进行剔除。具体阈值确定方法如下：
[0102]
假设大小为a
×
b的图像共分为{0，1，
…
，l-1}个灰度级，则第i个灰度级出现的概率为：
[0103][0104]
其中，ni表示第i个灰度级所包含的像素数量，设阈值为t，则前景和背景的灰度分布、均值和对应的概率分别为：
[0105][0106]
则最大类间方差σ
*2
为：
[0107]
σ
*2
＝p0(μ
0-u)2+p1(μ
1-u)2ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(13)
[0108]
其中进一步阈值可定义为：
[0109]
t
x
＝argmax
0＜k＜l-1
(σ
*2
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(14)
[0110]
即t
x
使得类间方差最大化，其中下标x取值为1或2，分别代表滤波与未滤波图像得到的阈值，则最终阈值t由如下公式计算得出：
[0111]
t＝αt1+(1-α)t2ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(15)
[0112]
其中α为阈值的权重，取值范围为[0,1]。其中，α越大，滤波的作用越明显，过度的滤波可能导致滑坡细节无法展现，同时可能导致滑坡的边界信息不明显；α越小，滤波作用越弱，图像可能受到噪声的干扰。因此，α的大小是决定了噪声平滑和滑坡目标信息保留之间的平衡，经过多次实验，当α值为0.6时，滑坡目标的信息不仅可以相对完整的保留，而且可以较好的解决噪声的干扰。
[0113]
3基于光学遥感图像的干扰区域提取
[0114]
3.1基于区域合并算法的多尺度分割
[0115]
本发明方法所使用的多尺度分割方法是ecognition软件中包含的区域合并法，其对象结构层次如图2。该方法基于空间和光谱的同质性，当分割尺度的平方大于所有同质性参数的总和时，两个像素被认为是同一个对象，反之则认为两者不属于同一对象，合并停止，影像对象的同质性定义如式(16)：
[0116][0117]
其中，wc和ws分别为颜色同质性和形状同质性所占的权重，hc为颜色同质性，它与组成对象的像素数目及每个波段的标准差有关，如式(17)：
[0118][0119]
其中，wc表为波段权重，nm为合并区域的像素数，为合并后区域的标准差，n
o1
、n
o2
分别为两个相邻区域的像素个数，分别为两个相邻区域的标准差；hs为形状同质性，它与紧凑度h
cm
与光滑度h
sm
有关，如式(18)：
[0120][0121]
其中w
cm
、w
sm
分别为紧凑度和光滑度对应的权重，h
cm
可表示为：
[0122]
[0123]hsm
可表示为：
[0124][0125]
其中lm为合并区域的周长，bm为合并区域最小外包矩形的周长。
[0126]
不同的地物类型都有其特定的空间尺度，在多尺度分割的过程中，分割参数的调整是非常重要的步骤，本文通过esp2工具进行分割尺度的优选。该方法通过计算在不同分割尺度下同一对象内部的局部方差来求取最优分割尺度，首先计算出分割对象的同质性局部方差(local variance,lv)，然后计算lv的变化率(rates of change,roc)。
[0127][0128]
通常情况下，lv曲线随着分割尺度不断变大而不断升高，roc曲线随着分割尺度不断变大呈现总体下降的趋势，在此过程中，roc曲线在不断下降的过程中出现局部极大值，即突然的峰值点，对应效果较好的分割尺度。
[0129]
3.2分类规则的建立
[0130]
本文光学图像的使用主要用于去除容易对滑坡检测产生干扰的农田、裸地和低植被区域。对于低植被和土壤背景而言，差值植被指数(difference vegetation index,dvi)的变化较为敏感且在低植被生物群落中具有更高敏感度和更广的动态取值范围，但是当植被覆盖度较高时，他对植被的灵敏度降低。dvi计算方式如式(22)：
[0131]
dvi＝ρ
nir-ρ
red
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(22)
[0132]
其中ρ
nir
为近红外波段值，ρ
red
为红波段值。
[0133]
进一步结合绿色波段来设计规则剔除滑坡初筛区域中的虚警，如图3，首先根据不同地物类型的dvi，将易对滑坡检测造成干扰的耕地及裸土区域提取出来。然后结合绿色波段(具体的，选择灰度级大于900且小于1600的绿色波段)，进一步剔除易对滑坡检测造成干扰的稀疏植被覆盖区域。
[0134]
为验证本发明效果，进行如下示例：
[0135]
4实验数据集的构建与介绍
[0136]
2018年9月6日北海道地震引起的山体滑坡(北纬42.7000
°
～42.7486
°
，东经141.9161
°
～141.9847
°
)，位置如图4。这次地震引发了至少9295处山体滑坡，这些滑坡分布广泛，造成了强烈的破坏。据报道本次地震事件造成44人死亡，其中36人死于地震引发的滑坡。此次地震引发的滑坡密度较高，主要发生在植被覆盖的山区，多为连续分布的浅层边坡破坏，大部分受灾地区植被茂密，滑坡其周边也包含滑坡、水系、森林、居民地、田地五类基本地物类型。
[0137]
数据的光学遥感图像为sentinel-2卫星获取的多光谱数据。灾前图像获取时间为2018年6月7日，分辨率为10m；灾后polsar数据为alos-2卫星搭载的palsar-2雷达获取到的全极化数据，获取时间是2018年9月8日，入射角为28.4
°
；具体参数如表1所示。
[0138]
表1多源遥感图像具体参数
[0139][0140]
如图5所示为本实施例研究区域的基本情况，其中灾后光学遥感图像是未经过几何校正的图像，用于参考比较；灾前光学遥感图像是经过几何校正的，用于剔除polsar图像预筛选滑坡区域中的虚警；灾后yamaguchi图像假彩色图像用于目视解译滑坡区域，目视解译结果为通过人工判读得出的滑坡区域结果，用于其他方法得出滑坡检测结果的评价基准。
[0141]
5 实验分析与讨论
[0142]
5.1 性能评价指标
[0143]
为了定量说明本发明方法的优越性能，通过总体精度(overall accuracy,oa)，生产者精度(producer’s accuracy,pa)，用户精度(user’s accuracy,ua)和kappa系数作为定量评价指标。其中oa指所有的正确分类的像素数占总体像素数的百分比，表示分类的总体准确度；pa表示从分类结果中任意取出一个随机样本，其代表的类别和真实样本一致的概率，在本实施例中，pa又分为滑坡区域生产者精度(producer’s accuracy of landslide,pal)和背景区域生产者精度(producer’s accuracy of background,pab)；ua表示真实样本数据被正确分类的概率，在本实施例中，ua又分为滑坡区域用户精度(user’s accuracy of landslide,ual)，背景区域生产者精度(user’s accuracy of background,uab)；kappa系数表示分类结果和样本之间的一致性，取值范围为[0,1]，取值越大说明分类得到的结果与样本真值的一致性越强。各项指标可作如下表示：
[0144][0145][0146][0147][0148]
式中，tp、fp、fn和tn是混淆矩阵中的元素。
[0149]
5.2滑坡检测结果对比
[0150]
为了说明本文方法的有效性，本文选择了文献[2](马战林,薛华柱,刘昌华,等.基于主被动遥感数据和面向对象的大蒜识别[j].农业工程学报,2022,38(02):210-222.)提出的多源遥感分类方法及文献[1](plank s,twele a,martinis s.landslide mapping in vegetated areas using change detection based on optical and polarimetric sar data[j].remote sensing,2016,8(4).)提出的滑坡区域检测方法进行对比试验。相比于文献[1]的方法，本文方法没有使用地学数据作为辅助，进一步降低了对数据的要求，因为合适的dem数据往往难以获取，这导致文献[1]的方法在实际应用中会受到限制。三种方法的
定量指标如表2，并绘制了柱状图如图6。
[0151]
表2不同方法的检测结果
[0152][0153]
由于文献[2]的方法是基于灾后光学遥感和sar图像的监督分类，该方法检测出了较多的滑坡区域且没有将滑坡区域过多的错判为背景，得到了较高的pal和uab，分别为83.81％和90.59％，同时由于难以避免其他类似滑坡区域造成的干扰，导致背景中存在大量区域被错误的检测为滑坡，因此该方法的pab和ual较低，仅为64.52％和49.44％。相较于文献[2]的方法，文献[1]的方法由于使用了滑坡前和滑坡后两个时相的图像，可以较好的剔除与滑坡散射特性类似但没被滑坡影响区域带来的干扰，因此文献[1]的方法得到了更高pab和ual，分别为94.94％和84.86％，同时由于dem的使用，使得部分滑坡区域未被成功检测，导致pal和uab较低，仅为43.89％和72.38％。相对于文献[2]提出的基于监督分类的方法，本文方法检测到的滑坡区域相对不够完整，得到了71.66％的pal值，但是本文方法大幅度降低了滑坡检测虚警，得到了与标签图像一致性更强的滑坡区域，所以本文方法的总体精度(87.80％)和kappa系数(85.81％)均高于文献[2]的方法和文献[1]的方法。
[0154]
图7为文献[1]的方法与本文方法的滑坡区域检测结果图。由图7可以看出，在文献[1]的方法中，由于dem的影响，滑坡周边坡度较小的区域并未被成功检测，而本文方法检测出了更多的滑坡区域。由于一些山坡上区域和滑坡的散射特性非常相似，同时具有较大的坡度，文献[1]的方法无法有效避免该类区域对滑坡检测的干扰，而本文方法则对可能对滑坡产生干扰的进行了更加细致的分类，对滑坡虚警的剔除作用更加明显。
[0155]
在上述实施例的基础上，本发明还公开一种联合灾前光学与灾后polsar图像的滑坡检测系统，作为一种可实施方式，系统需要通过图形化的界面更直观地显示滑坡检测结果，因此系统需包含以下几个功能：
[0156]
(1)输入遥感图像。通过菜单选择需要测试的遥感图像，在gui界面显示遥感图像。
[0157]
(2)显示滑坡检测结果。加载输入的遥感图像，用本文提出的方法进行滑坡检测，在gui界面显示滑坡检测结果。
[0158]
(3)显示评价分析结果。为直观表达模型性能，需对滑坡检测结果进行评价分析，在gui界面显示评价分析结果。
[0159]
系统利用matlab进行程序的开发和编程，gui界面采用matlab app designer进行设计显示。根据系统的功能需求，遵循操作简单便捷的原则，系统框图如图8所示。本系统主要包含三个模块，即数据输入、滑坡检测和评价分析。
[0160]
如图9所示为系统主界面图，根据右侧操作流程依次选择polsar数据输入和光学数据输入，然后通过图像大小的设置可对图像进行裁剪，使得目标区域视觉效果更佳，根据不同的数据类型和任务需求也可以改变高斯滤波的窗口大小和σ值，改变阈值的权重设置，改变形态学处理的窗口大小。如果没有现成标签数据，可以通过系统中的感兴趣区域选择功能进行标签制作，然后进行标签输入；如果有现成标签数据则可以直接输入。最后点击运
行，得到滑坡检测图像和评价指标，如图10，其中figure 1窗口显示为polsar假彩色图像，figure 2窗口显示为标签图像，figure 3窗口显示为滑坡检测结果。
[0161]
进一步地，滑坡检测模块具体包括：
[0162]
疑似滑坡区域提取单元，用于基于polsar提取疑似滑坡区域：对灾后polsar图像进行预处理，然后计算出滑坡区域检测过程中所需的极化特征参数，通过自动阈值确定方法确定分割阈值，进而对其进行形态学处理得到疑似滑坡区域；
[0163]
干扰区域提取单元，用于基于光学遥感图像提取干扰区域：以灾后polsar图像为基准对光学遥感图像进行几何校正，并进行多尺度分割，最后通过设计的分类规则筛选出容易对滑坡检测造成虚警的干扰区域；
[0164]
最终滑坡区域确定单元，用于在疑似滑坡区域的基础上，结合光学遥感图像得到的干扰区域进行虚警剔除，确定最终滑坡区域。
[0165]
进一步地，所述计算出滑坡区域检测过程中所需的极化特征参数包括：
[0166]
对预处理后灾后polsar图像的每一个像素计算香农熵。
[0167]
进一步地，所述基于polsar提取疑似滑坡区域还包括：
[0168]
通过yamaguchi分解得到的二次散射来剔除建筑物带来的虚警。
[0169]
进一步地，按照如下方式确定分割阈值t：
[0170]
t＝αt1+(1-α)t2[0171]
其中α为阈值的权重，t1、t2分别代表滤波与未滤波图像得到的阈值，式中
[0172]
t
x
＝argmax
0＜k＜l-1
(σ
*2
)
[0173]
其中x取值为1或2，σ
*2
为最大类间方差
[0174]
σ
*2
＝p0(μ
0-u)2+p1(μ
1-u)2[0175][0176][0177]
其中p0、p1分别表示前景和背景对应的概率；μ1、μ2分别表示前景和背景的均值；pi表示第i个灰度级出现的概率；ni表示第i个灰度级所包含的像素数量。
[0178]
进一步地，所述基于光学遥感图像提取干扰区域具体包括：
[0179]
以polsar图像为基准，通过地面控制点法对光学遥感图像进行几何校正。
[0180]
进一步地，所述基于光学遥感图像提取干扰区域具体包括：
[0181]
基于ecognition软件对几何校正后的光学遥感图像进行多尺度分割。
[0182]
进一步地，在所述多尺度分割过程中，通过esp2工具进行分割尺度的优选。
[0183]
进一步地，所述分类规则包括：
[0184]
判断不同地物类型的差值植被指数是否小于1500，若是，则将其分类为耕地及裸土区域，若否，则结合绿色波段，判断不同地物类型的差值植被指数是否位于1500至4000之间，若是，则将其分类为稀疏植被覆盖区。
[0185]
综上，针对当前基于灾前光学和灾后polsar图像的滑坡检测方法需要地学数据辅助且检测精度有待提升的问题，本发明提出了一种仅利用灾前光学遥感图像和灾后polsar图像的快速滑坡检测方法及系统。该方法通过自动阈值方法确定滑坡分割阈值，进而分割
得到疑似滑坡区域，可以在避免polsar图像噪声影响的情况下得到较为完整的滑坡区域；然后通过面向对象分类的方法筛选出容易对滑坡检测造成干扰的区域；最后通过逻辑运算结合光学和polsar图像的阶段性结果，得到滑坡区域。在数据集的实验中，该方法总体精度为87.80％，滑坡生产者精度为71.66％，说明了本发明方法在无需地学数据辅助的情况下，可以有效检测出滑坡区域。
[0186]
以上所示仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

技术特征：
1.一种联合灾前光学与灾后polsar图像的滑坡检测方法，其特征在于，包括：基于polsar提取疑似滑坡区域：对灾后polsar图像进行预处理，然后计算出滑坡区域检测过程中所需的极化特征参数，通过自动阈值确定方法确定分割阈值，进而对其进行形态学处理得到疑似滑坡区域；基于光学遥感图像提取干扰区域：以灾后polsar图像为基准对光学遥感图像进行几何校正，并进行多尺度分割，最后通过设计的分类规则筛选出容易对滑坡检测造成虚警的干扰区域；在疑似滑坡区域的基础上，结合光学遥感图像得到的干扰区域进行虚警剔除，确定最终滑坡区域。2.根据权利要求1所述的一种联合灾前光学与灾后polsar图像的滑坡检测方法，其特征在于，所述计算出滑坡区域检测过程中所需的极化特征参数包括：对预处理后灾后polsar图像的每一个像素计算香农熵。3.根据权利要求1所述的一种联合灾前光学与灾后polsar图像的滑坡检测方法，其特征在于，所述基于polsar提取疑似滑坡区域还包括：通过yamaguchi分解得到的二次散射来剔除建筑物带来的虚警。4.根据权利要求1所述的一种联合灾前光学与灾后polsar图像的滑坡检测方法，其特征在于，按照如下方式确定分割阈值t：t＝αt1+(1-α)t2其中α为阈值的权重，t1、t2分别代表滤波与未滤波图像得到的阈值，式中t
x
＝argmax
0＜k＜l-1
(σ
*2
)其中x取值为1或2，σ
*2
为最大类间方差σ
*2
＝p0(μ
0-u)2+p1(μ
1-u)
22
其中p0、p1分别表示前景和背景对应的概率；μ1、μ2分别表示前景和背景的均值；p
i
表示第i个灰度级出现的概率；n
i
表示第i个灰度级所包含的像素数量。5.根据权利要求1所述的一种联合灾前光学与灾后polsar图像的滑坡检测方法，其特征在于，所述基于光学遥感图像提取干扰区域具体包括：以polsar图像为基准，通过地面控制点法对光学遥感图像进行几何校正。6.根据权利要求1所述的一种联合灾前光学与灾后polsar图像的滑坡检测方法，其特征在于，所述基于光学遥感图像提取干扰区域具体包括：基于ecognition软件对几何校正后的光学遥感图像进行多尺度分割。7.根据权利要求6所述的一种联合灾前光学与灾后polsar图像的滑坡检测方法，其特征在于，在所述多尺度分割过程中，通过esp2工具进行分割尺度的优选。8.根据权利要求1所述的一种联合灾前光学与灾后polsar图像的滑坡检测方法，其特征在于，所述分类规则包括：判断不同地物类型的差值植被指数是否小于1500，若是，则将其分类为耕地及裸土区
域，若否，则结合绿色波段，判断不同地物类型的差值植被指数是否位于1500至4000之间，若是，则将其分类为稀疏植被覆盖区。9.一种联合灾前光学与灾后polsar图像的滑坡检测系统，其特征在于，包括：疑似滑坡区域提取单元，用于基于polsar提取疑似滑坡区域：对灾后polsar图像进行预处理，然后计算出滑坡区域检测过程中所需的极化特征参数，通过自动阈值确定方法确定分割阈值，进而对其进行形态学处理得到疑似滑坡区域；干扰区域提取单元，用于基于光学遥感图像提取干扰区域：以灾后polsar图像为基准对光学遥感图像进行几何校正，并进行多尺度分割，最后通过设计的分类规则筛选出容易对滑坡检测造成虚警的干扰区域；最终滑坡区域确定单元，用于在疑似滑坡区域的基础上，结合光学遥感图像得到的干扰区域进行虚警剔除，确定最终滑坡区域。

技术总结
本发明公开一种联合灾前光学与灾后PolSAR图像的滑坡检测方法及系统，在PolSAR数据处理部分，首先对灾后PolSAR数据进行预处理，然后计算出滑坡区域检测过程中所需的极化特征参数，通过自动阈值确定方法确定分割阈值，进而对其进行形态学处理得到疑似滑坡区域。在光学数据处理部分，首先以灾后PolSAR图像为基准对光学图像进行几何校正，并进行多尺度分割，最后通过设计的分类规则筛选出容易对滑坡检测造成虚警的干扰区域。最后在疑似滑坡区域的基础上，结合光学遥感图像得到的干扰区域进行虚警剔除，确定最终滑坡区域。本发明方法在无需地学数据辅助的情况下，可以有效检测出滑坡区域。出滑坡区域。出滑坡区域。


技术研发人员：牛朝阳 高欧阳 王世举 吴亚杰 李卓越 王艳云
受保护的技术使用者：中国人民解放军战略支援部队信息工程大学
技术研发日：2023.06.09
技术公布日：2023/10/15