高欧阳， 牛朝阳， 刘 伟， 张浩波， 邹玮琦， 张雅歌

(中国人民解放军战略支援部队信息工程大学数据与目标工程学院, 河南郑州 450001)

0 引 言

利用光学或合成孔径雷达(Synthetic Aperture Radar, SAR)遥感卫星影像检测灾后滑坡的位置和边界，调查滑坡的空间分布特征，能够为灾害防治与灾区重建提供重要信息[1]。滑坡的形成常常伴随有恶劣的天气情况，其中降雨是滑坡最主要的诱发因素之一[2]，阴雨天气云层的遮挡会导致光学遥感卫星探测手段的失效。SAR作为一种主动式对地观测系统，通常工作在微波波段，波长较长，基本不受太阳光照、气候条件等环境因素的影响，能够进行全天时、全天候的不间断观测[3]。利用极化合成孔径雷达(Polarimetric Synthetic Aperture Radar，PolSAR)图像进行滑坡检测已经成为重要的研究方向。

PolSAR通过发射和接收极化状态正交的电磁波，能够获取地物的全极化信息[4]，基于极化分解技术对地物散射机理进行建模与分析，能够准确地理解地物的物理特性[5]。在极化分解的基础上进行滑坡检测主要有两种思路[6]，一是利用多阈值法进行滑坡检测，即对多项极化特征分别设置阈值，同时满足所有阈值条件的区域被确定为滑坡区域，此类方法具有简单易行的优点，但是在实施仅限于特定环境并需要人工阈值的干预，普适性不强；二是利用监督分类或非监督分类的方法对灾后PolSAR图像进行分类以区分出属于滑坡的区域。相比于多阈值法，后者得益于分类器的复杂结构，可以提供更高的滑坡检测准确率，因此分类法比多阈值法更常用于滑坡检测。Niu等[6]在总结前人工作的基础上，首次将变化检测(Change Detection，CD)与层次分析法(Analytic Hierarchy Process，AHP)相结合实现复杂地貌背景下的滑坡检测。该方法用CD检测发生变化的区域，即预筛选滑坡区域，然后通过AHP剔除其他变化区域对滑坡检测造成的干扰，但是从复杂地貌背景中检测滑坡的准确率为73.25%，田地对滑坡检测造成的虚警率为14.74%，还有较大的改善空间。

综上，极化测量技术可以充分利用极化信息来区分不同散射特性的地貌，但是对于散射特性相似的地貌难以取得理想的效果。因此，本文在前人研究的基础上，首先通过优劣解距离法(Technique for Order Preference by Similarity to an Ideal Solution，TOPSIS)综合多个极化特征参数，充分利用PolSAR图像所包含的极化信息提取出与滑坡散射特征相似的区域，实现对疑似滑坡区域的初步筛选；然后通过相干图对地形地貌变化非常敏感的特性来确定发生变化的区域，从而抑制疑似滑坡区域中未变化区域造成的虚警，实现从复杂地貌背景中有效检测滑坡的目的。

1 研究方法

本文方法的处理流程如图1所示，包含3个阶段。第一阶段是预处理过程，对灾后PolSAR数据进行辐射定标、多视和地形校正。第二阶段为滑坡检测过程，首先通过Yamaguchi分解、H/A/α分解和相关系数计算得到表面散射ps、二次散射pd、体散射pv、螺旋体散射ph[7]、极化熵H、平均散射角α、反熵A[8]和同极化分量相关系数的实部Re(ρhh-vv)8个极化特征参数，然后利用带权重的TOPSIS法充分融合上述极化特征参数信息，经阈值分割和形态学处理得到与滑坡极化散射特性相同的区域(下文均称为“疑似滑坡区域”)。第三阶段为虚警抑制过程。利用灾前单极化SAR数据与灾后PolSAR数据的相同极化组合图像构建双时相数据，经干涉处理生成相干图以提取出变化区域，然后将通过逻辑“与”运算将两个结果进行融合，从而去除与滑坡表面散射行为一致的未变化区域造成的大量检测虚警。

图1 本文方法实现滑坡检测的流程图

本文方法中的形态学运算为先“开”后“闭”运算，即首先通过“开”运算去除二值化图像黑色背景中的白色斑点(如椒盐噪声)，然后通过“闭”运算来消除白色区域内的黑色空洞、连接邻近的白色区域并平滑边界[9]，从而使检测结果图的区域连通性更好，边界特征更加明显。下文重点阐述极化特征参数提取、TOPSIS检测与相干图计算等主要处理步骤。

1.1 极化特征参数提取

1.1.1 Yamaguchi分解

Yamaguchi分解是一种基于散射模型的非相干分解算法，可以很好地匹配地物的基本散射机制。该分解方法从极化相干矩阵T中提取出各个散射分量对应的功率大小，即表面散射功率Ps、二次散射功率Pd、体散射功率Pv和螺旋体散射功率Ph。最后，极化相干矩阵T可以表示为各散射功率的加权组合形式：

T=fsTs+fdTd+fvTv+fhTh

(1)

式中：Ts，Td，Tv，Th分别为表面散射、二次散射、体散射和螺旋体散射分量对应的极化相干矩阵；fs，fd，fv，fh分别为Ps，Pd，Pv和Ph的分解系数[7]。

1.1.2H/A/α分解

H/A/α分解是一种利用协方差矩阵或者极化相干矩阵的特征值和特征向量定义极化特征的非相干分解算法。在单站互易条件下，极化相干矩阵T为3×3的半正定厄密特矩阵[10]，因此可以作如下特征值分解：

(2)

式中，特征值λi(1≤i≤3)满足λ1≥λ2≥λ3>0，λi对应的特征向量ui相互正交，H表示共轭转置。则极化熵H、平均散射角α和反熵A可由特征值λi和特征向量ui计算得到：

(3)

1.1.3 相关系数计算

同极化分量相关系数ρhh-vv为目标的两个同极化散射矢量之间的相关系数，取值越大表明目标的散射类型越单一，例如Re(ρhh-vv)的正值部分常用于反映地物的表面散射特性[11]。通过极化相干矩阵T计算Re(ρhh-vv)的公式如下：

Re(ρhh-vv)=

(4)

式中，Re(·)表示取实部，T11，T12和T22表示极化相干矩阵T中的相应元素。

1.2 TOPSIS提取疑似滑坡区域

1.2.1 TOPSIS用于PolSAR图像滑坡检测的基本原理

TOPSIS属于多属性决策分析中的常用方法，也可以称为逼近理想解法，其对数据分布及样本含量没有严格限制，能够充分利用原始数据信息，并且评价结果比较精确可靠。该方法根据一系列属性条件(即评价指标)对有限个评价对象与理想解的接近程度进行排序[12]。假设采用M个评价指标对N个评价对象分别评价，其中第i个(i=1,2,3,…,N)评价对象在第j个(j=1,2,3,…,M)评价指标上的数值为zNM，可以得到评估矩阵如下：

(5)

进一步定义评估矩阵的正、负理想解分别为

(6)

(7)

第i个评价对象与正、负理想解的距离分别为

(8)

(9)

式中wi表示每个评价指标的权重，可以根据不同的需求通过不同的方法来确定。

最后计算评价对象的最终得分，即其与正理想解的相对贴近度为

(10)

相对贴近度Si的大小反映了评价对象与正理想解的差距，Si越大说明评价对象越接近正理想解，Si越小说明评价对象越远离正理想解。

利用TOPSIS将8个极化特征参数(ps，pd，pv，ph，H，A，α，Re(ρhh-vv))有机融合，可以在充分利用极化信息的基础上精确计算PolSAR图像每个像素与滑坡的相似度，即相对贴近度。如图2所示，将每一个像素点作为一个评价对象，极化特征参数作为评价对象得分的一个评价指标，然后利用TOPSIS计算出每个像素与理想解的相对贴近度Si，然后选择合适的阈值对Si进行二值化处理即可得到疑似滑坡区域。其中，为了确保TOPSIS方法对不同PolSAR图像数据的适用性，确定正理想解为1，负理想解为0。

图2 TOPSIS提取疑似滑坡区域的关键步骤

1.2.2 极化特征参数的归一化和同趋势化处理

因为8个极化特征参数(ps，pd，pv，ph，H，A，α，Re(ρhh-vv))作为评价指标的指向性、计量单位、数值评价方式和评价值范围存在的差异会对评价结果产生不利影响，所以需要在构建评估矩阵之前对这些极化特征参数进行归一化和同趋势化处理，使各项参数数值都在区间[0,1]内，并且与滑坡散射特性成正相关关系，即数值越大，越能体现滑坡的散射特性。

1) 对Yamaguchi四分量进行归一化处理[6]。

px=Px/(Ps+Pd+Pv+Ph)

(11)

式中，下标x∈{s,d,v,h}，s，d，v和h分别表示表面散射、二次散射、体散射和螺旋体散射。由于滑坡发生后地表覆盖的裸露土壤[11]通常以表面散射的形式所呈现，与其他地物有很大的区别[6]，所以ps值越大、pd,pv,ph的值越小，越能体现滑坡地貌的表面散射类型，因此本文对归一化处理后的pd,pv,ph进行正向化处理。

px=1-px

(12)

2) 对于极化熵H、平均散射角α，其中当H属于0.52～0.63且α属于29～37的取值范围时更能体现滑坡的散射特性[13]，所以本文利用区间型归一化方式来对H和α的值进行处理。

(13)

rij=1

(14)

3) Re(ρhh-vv)的正值部分常用于识别表面散射地物类型，其值越大表明地物表面散射类型更为显著，所以将Re(ρhh-vv)的负值部分赋值为0。

1.2.3 AHP方法确定极化特征参数的权重

不同的极化特征参数对于最终评分的影响不同，因此需要对它们赋予不同的权重。AHP是美国运筹学家Thomas提出的一种分析复杂决策的半定量方法，Niu等[6]将AHP运用到PolSAR图像滑坡检测中。在此，本文使用AHP来确定极化特征参数的权重。根据滑坡灾后的地貌特征，ps最能体现滑坡的表面散射特性，H和α可以较好地区分滑坡与其他地貌[13]，Re(ρhh-vv)常用于反应地貌的表面散射特性，pd，pv，ph和A将对识别滑坡作用较小。因此将ps作为最重要的评价指标，H和α作为比较重要的评价指标，Re(ρhh-vv)作为一般重要的评价指标，pd，pv，ph和A作为相对不重要的评价指标，然后根据文献[14]中的成对比较矩阵重要性标度表确定不同评价指标之间的重要性比较值，构建判断矩阵如表1所示。为了尽可能消除判断矩阵中存在的主观性，需要引入一致性比率(CR)作为参考，当CR<0.1时判断矩阵被认为是合理的。经过计算，表1所示的判断矩阵CR=0，说明各个评价指标的权重分配是合理的。

表1 AHP确定权重的判断矩阵

1.3 基于相干图抑制滑坡检测的虚警

在利用TOPSIS方法从复杂地貌背景中检测滑坡的结果中，田地、裸土等其他以表面散射类型为主的地物地貌，可能对滑坡检测结果带来较多的虚警。滑坡作为一种突发性自然灾害，灾害区域具有十分显著的突变性，而田地、裸土等地物地貌一般不具有这种特征，因此可通过对地物地貌的突变性分析抑制滑坡检测的虚警。地形地貌的显著变化通常会在SAR干涉处理中引起严重的失相干，即相干性较低，而未发生变化的区域将保持较高的相干性，因此可以通过相干图来判断某一区域是否发生变化，进而抑制未变化地物造成的虚警，使滑坡检测结果更加准确。

利用双时相SAR图像可以得到相干图，相干图中灰度值的大小表示相干性的高低，相干性|γ|计算如下：

(15)

其中，E[ ]表示数学期望，S1，S2分别为滑坡前后的SAR影像(极化组合相同)，相干性|γ|的取值范围为[0,1]，|γ|越大说明一个区域的变化程度越低，|γ|越小说明该区域的变化越明显[15]。

在不同极化组合的SAR图像中，垂直极化后向散射更能体现以水平结构为主的地貌，因此本文采用能突出裸土、裸岩区域特征[16]的VV极化组合图像生成研究区域的相干图，然后对相干图进行阈值分割和形态学处理得出变化区域。

2 研究区域与数据源

研究区域包括两个。区域A：2015年12月20日由于暴雨引起的山体滑坡发生于深圳恒泰裕工业园区，本次滑坡发生在复杂的地貌背景中(包含滑坡、水系、森林、居民地、田地五类基本地物)，并且存在其他自然或人类活动引起的变化，因此使用该区域作为研究对象具有一定的代表性。区域B：2018年9月6日北海道地震引起的山体滑坡，滑坡主要发生在植被覆盖的山区，其周边也包含滑坡、水系、森林、居民地、田地五类基本地物类型，但是由于滑坡区域和其他地物距离较远，在同一幅图上展示会因比例尺过小导致目标区域难以判读，为了分析不同地貌对滑坡区域检测的影响，本文从同一景数据中截取不同地貌区域合并为同一幅图像。两个研究区域均使用日本ALOS-2卫星的全极化数据，距离向分辨率为5.13 m，方位向分辨率为3.21 m，具体参数如表2所示。

表2 卫星数据基本参数

两个研究区域的Yamaguchi分解假彩色图像和光学图像如图3所示。其中假彩色图像的红色代表二次散射，绿色代表体散射，蓝色代表表面散射。受灾前后的光学图像与PolSAR图像的采集时间存在一定的差别，但是其仍然能够作为分析讨论滑坡引起地形地貌变化的参考。通过目视判读后在假彩色图上利用白色框选出滑坡、水系、森林、居民地、田地五类地物样本，以便于在下文作出检测正确率分析。

(a) 研究区域A的假彩色图像和光学图像

(b) 研究区域B的假彩色图像和光学图像图3 研究区域基本情况图

3 实验结果与分析

3.1 滑坡检测性能的定量分析

为了定量分析本文方法的有效性，通过图3标记区域的不同地物像素样本计算滑坡检测正确率和总虚警率(即滑坡以外地物样本区域内被检测为滑坡的像素总数除以样本区域像素总数)，与文献[6]的CD+AHP方法作出对比，并在图4中描绘出不同方法在两个研究区域的正确率和总虚警率关系曲线。其中4条曲线的横坐标虚警率均小于1，这是由于CD+AHP方法检测的虚警被限制在变化检测区域范围内，本文方法检测的虚警被限制在相干图分割区域范围内。因为地形地貌上的差异导致相同方法在不同研究区域有着不同的性能表现，对于每个研究区域，本文方法和CD+AHP方法相比，在相等虚警率的情况下滑坡检测正确率均有大幅提升，表明了其检测滑坡的有效性。

为了进一步分析不同方法的滑坡检测性能，选择合适的阈值计算本文方法和CD+AHP方法在相同正确率下不同地物造成的虚警率，以及相同总虚警率下滑坡检测正确率和不同地物造成的虚警率，并与文献[11]的MT方法作出对比，如表3所示。其中，CD+AHP方法的滑坡检测正确率和虚警率均由文献[6]建议的阈值得到，MT方法中的多阈值为文献[11]设定的经验值，该阈值无法调整得到与另外两种方法相同的滑坡检测正确率和总虚警率。

图4 不同方法的滑坡检测正确率和总虚警率关系曲线图

表3 不同方法滑坡检测正确率和虚警率对比

对于研究区域A，MT方法的检测结果给出了最高的滑坡检测正确率，但是田地和水系等以表面散射为主的地物造成了极大的虚警率，使得该方法几乎失效；CD+AHP方法则以检测正确率小幅度降低为代价，较好地抑制了各类地物的虚警率，尤其是田地和水系的虚警率大幅度降低；与CD+AHP方法相比，在相同的正确率(73.46%左右)下，本文方法将总虚警率和各类地物的虚警率控制在极低的水平(均低于1%)；在相同的总虚警率(3.45%左右)下，本文方法的滑坡检测正确率提高了约11%，并且较好地抑制了大部分地物的虚警率，只有田地造成的虚警率较高，这是由于田地区域的农作物非常容易受到自然生长或者生产活动的影响而发生变化，因此通过相干图并不能完全剔除田地区域带来的虚警。

对于研究区域B，由于滑坡区域内大量存在以体散射和二次散射为主的残破树木、树干的干扰，导致3个方法的滑坡检测正确率均有明显降低。其中，MT方法检测正确率仅为55.45%，表明该方法的检测效果较差；CD+AHP方法提供了差强人意的68.61%检测正确率，同时对虚警的抑制也有明显提升；在相同的正确率(68.61%左右)下，本文方法与CD+AHP方法相比，总虚警率和各类地物虚警率都降低至理想水平，在与CD+AHP方法相同的总虚警率(8.18%左右)下，本文方法的检测正确率提高了约16%，达到了较为合理的数值(84.48%)。

3.2 滑坡检测结果图的分析讨论

为了更加直观地分析不同方法滑坡检测的效果，将表3对应的结果图在图5中给出，其中本文方法的结果图对应表3中上标②所在行的数值。

(a) 研究区域A不同方法的滑坡检测结果图

(b) 研究区域B不同方法的滑坡检测结果图图5 不同方法滑坡检测结果图对比

通过对比Yamaguchi假彩色图像和光学图像目视解译得到滑坡区域范围(在图5中用黄色标记)。对于研究区域A，三种方法均能将滑坡与二次散射为主的居民地和体散射为主的森林有效地区分开，同时在滑坡区域存在散落分布的小面积黑色空洞，这是由于救灾活动中该区域存在大量的大型机械(比如挖掘机、铲车、运土车等)，机械臂、车体与地面之间形成二次散射结构，从而对滑坡检测造成了一定的干扰，但是并不影响大面积滑坡范围的确定。MT方法检测结果中，田地和水系等以表面散射为主的地物造成了大量的虚警，对滑坡区域的提取造成了很大影响；CD+AHP方法检测结果中，森林、居民地等地物类型造成的虚警均较少，但是明显将少部分田地和水系错误地检测为滑坡，同时滑坡头部(区域1)存在明显空洞，这是由于滑坡前后该区域均为裸土，主要散射类型没有发生变化所致。本文方法得出的滑坡区域更加完整(包括区域1)且连通性更强，而且背景中大部分地物的虚警得到了更好的抑制，但是本文方法在田地(区域2)的虚警稍多于CD+AHP方法，这是由于农作物的生长导致了该区域相干性较低，通过相干图不能完全抑制该区域的虚警；区域3的森林覆盖率较低，中间夹杂了少量裸土区域，因为这些区域在滑坡前后未发生明显变化，因此MT方法仅通过滑坡后单时相PolSAR数据错误地将此类区域错判为滑坡，而CD+AHP方法和本文方法的虚警抑制效果较好，这是因为两种方法采用不同的方式(CD和相干图)来检测研究区域中的变化，进而剔除未变化区域的虚警。

对于研究区域B，MT方法和CD+AHP方法均难以在提取出完整的滑坡区域信息，本文方法则可以在复杂背景中有效地抑制虚警的同时提供较好的滑坡检测结果，两个研究区域检测结果的一致性证明了本文方法可以实现复杂背景下检测滑坡的目的。

4 结束语

针对现有方法在复杂背景下滑坡检测正确率低、虚警率高的问题，本文提出了一种利用TOPSIS融合极化信息的滑坡检测方法，通过综合多个极化特征参数，以不同地貌的散射特性为基础来检测滑坡，进一步使用相干图剔除虚警，可以得到清晰的滑坡区域范围，大幅度提高滑坡检测正确率，而且可以有效抑制大部分地貌造成的虚警。随着PolSAR技术的不断发展，利用PolSAR数据进行灾害检测将会发挥越来越重要的作用，本文方法在不依赖先验信息的条件下可以准确快速地提取出滑坡发生后灾区的信息，为复杂背景下的滑坡检测提供了具体可行的方案，为灾后的救援、二次灾害的预警和灾后重建等活动提供可靠的依据，为今后相关方面的研究提供了一定的参考。