1 引言
在我国常见的各类地质灾害中,滑坡数量占地质灾害总量的70%以上是发生在山区最主要的地质灾害类型[1]。滑坡作为一种极具破性的地质灾害,严重制约着灾害多发地区的经济发展,并对社会公共安全构成严重威胁。开展滑坡体形变监测对了解坡体演化规律及准确预测坡体发展趋势具有重要意义[2]。目前对滑坡表面形变进行测量的方式有很多,然而,传统的精密水准测量、全球卫星定位(Global Positioning System, GPS)系统、伸缩计、倾斜仪等技术仅能获取单个监测点的变形信息,很难满足监测的实际需要,星载平台的合成孔径雷达(Synthetic Aperture Radar,SAR)系统,受雷达侧视角度、重访周期和大气延迟的限制,难以满足实时性需求。作为一种主动式微波成像传感器,地基SAR具有全天时、全天候、连续监测、方便携带等优点,通过在不同时间点对目标区域的重复观测获取时间序列干涉数据用于形变监测,可达到精度为毫米级的形变监测结果,是对局部区域形变监测的一种技术手段[3]。
永久散射体(Permanent Scatterers, PS)干涉测量技术是地基SAR克服时间去相干等非理想因素获取高精度形变信息的有效手段。常用方法有Ferretti等[4-5]提出的利用幅度离差信息选择PS点,Berardino等[6]提出相关系数法,以上两种PS点探测方法均属于单阈值法,但这种基于单一信息的PS点选择方法只强调PS点的某一个特性,故容易造成错判或漏判。国内学者陈强[7],朱茂[8],胡波[9],黄其欢[10]等将幅度离差、相关系数与其他SAR图像信息,如相位误差,振幅指数等相结合,分别提出了双阈值和多阈值的PS点探测方法,这些多阈值的方法有效的减少了错选和误选现象。为了能在不同场景的观测过程中都能高质量地选出 PS 点,需要在形变监测开始前重新调试 PS 选择阈值,虽然多阈值法提升了所选PS点集合的质量,但是人工设置多个阈值易出现阈值不合理现象,若设置的阈值不合理则会导致 PS 选取的质量和数量出现偏差, 进而影响形变反演结果的精确度,且反复调试多个阈值的过程限制了地基SAR 形变测量技术的推广。因此,需要研究多阈值PS点探测的优化算法。
由于多阈值PS点探测方法人工设置阈值不具最优性,本文在基于幅度离差、相干系数和相位误差的多阈值算法基础上提出一种多阈值优化算法,目的是自动确定最优阈值,将严苛阈值条件下筛选出的稳定点形变作为参考形变,与阈值放宽后稳定点邻域平均形变比较,当二者差值超过σ时将当前阈值定为最优阈值,此方法极大程度的省去繁琐的多阈值调试过程。自动确定的阈值筛选效果能达到比手动设置的阈值更优的效果,能满足形变监测需求,此算法对地基SAR形变测量技术的推广具有重大意义。本文其余部分安排如下,第2节介绍PS点选择算法的三种子方法,第3节介绍本文提出的阈值优化算法的流程和框架,最后利用四川茂县滑坡数据验证算法的正确性。
2 幅度离差、相干系数和相位误差的三种PS点选择方法
将同时满足幅度离差、相干系数和相位误差所设门限的点选为PS点,这种用多重信息融合来进行选择的方法充分考虑了PS点的多种特性,选出的PS点集合具有相位质量高且稳定的特点。
2.1 幅度离差PS点选择
幅度离差法是选择PS点的经典方法,求取雷达时序影像中每个像元的幅度均值和标准差,幅度离差指数计算公式如下[4]:
(1)
其中,DA表示目标点时间序列的幅度离差指数,δA表示目标点幅度时间序列的标准差,mA表示目标点幅度时间序列均值。幅度离差法只利用雷达时序影像的幅度信息,效率高,计算量小,幅度离差越小,表示幅度信息越稳定,设置幅度离差阈值TD,将满足条件DA<TD的像素点选为PS点。
2.2 相干系数PS点选择
相干性是衡量干涉相位噪声的最直观的标准,相干系数法是根据目标像素四周临近的像素值来估计其相干系数,具体的相干系数表达式为[11]
(2)
其中,γ为第x个像素相干系数,M(i, j)和S(i, j)分别为构成第k个干涉对的主从影像,*表示复数共轭,m和n为滑动窗口大小,本文窗口大小设置为5×5。在N个地基SAR影像干涉对的情况下,由式(2)可计算各干涉对的相干系数(γ1,γ2,γ3,...,γN),为了选取PS候选点,将逐像元计算平均相干系数,即
(3)
相干系数是依赖移动窗口进行计算的,移动窗口的大小对相干系数的估算有非常大的影响,选择的窗口越大,相干系数越可靠,但是会导致分辨率降低,一些孤立的点会被漏选,而接近PS点的一些非稳定目标会被错选为PS点。设置的阈值越高,所选PS点的可靠性就越高,但漏选的概率也越高,设置合适的相干系数阈值
选择
的点为PS点。使用相干系数法可有效去除植被、水体等失相干严重的像元。
2.3 相位误差阈值
SAR影像的干涉相位Δφ可以建模为[4]:
Δφ=φdef+φatm+φnoise
(4)
其中,φdef为形变相位,φatm为大气相位,φnoise为噪声相位。在小规模中,变形阶段具有一定的空间相关性。 同时,也可以认为大气相在空间范围内缓慢变化。 但对于噪声阶段,不同目标之间的统计属性是完全独立且无关紧要的。 因此,可以通过空间滤波滤除相位噪声。 然后从原始相位中减去滤波后的相位,得到噪声相位,数学描述为[8]:
φnoise=Δφ-Filter(Δφ)
(5)
其中, φnoise为噪声相位,Δφ为干涉相位,Filter表示均值滤波。设置相位噪声阈值TN,选择φnoise<φN的像素点为PS点。
3 三阈值优化算法
为使阈值设置更具合理性,减少误选,本文提出一种多阈值优化算法,算法首先在严苛阈值条件筛选出稳定点,并计算出稳定点形变量作为参考形变量,接着,逐步放宽阈值,使稳定点邻域出现其他PS点,选定合适窗口计算此窗口内的形变均值,形变均值与参考形变值之差不超过σ的情况下PS点数量最多的阈值组合即为最优阈值,此算法可减少阈值设置的不合理性,自动确定最优阈值,算法实现步骤如下:
(1)根据公式(1)计算时间序列幅度离差,根据公式(2)、(3)计算时间序列数据的相干性,根据本公式(4)、(5)计算时间序列相位误差;
(2)设置严苛幅度离差阈值DA=T′D,设置严苛相干系数阈值
设置严苛相位误差阈值φnoise=T′N,选出稳定的参考点,计算其形变量,作为参考形变量D′;
(3)选择0.1为步长放宽幅度离差和相关系数阈值,设置
使参考点邻域出现新的PS点,相位误差阈值不变φnoise=T′N;
(4)选择合适窗口计算阈值放宽后参考点及其窗口中参考点邻域新出现的PS点的平均形变量D;
(5)重复步骤(3)～(4),得到每种阈值组合情况后,判断每种组合情况的平均形变量|D′-D|是否小于σ,选择|D′-D|<σ情况下PS点数量最多的阈值组合为最优阈值组合;组合的个数等于幅度离差阈值放宽的次数乘以相干系数阈值放宽的次数;
(6)此时最优阈值可反演高精度形变。
流程图如图1所示。
图1 三阈值优化算法流程图
Fig.1 Three threshold optimization algorithm flow chart
4 实验分析
4.1 实验数据
2017年6月24日茂县西北部的新磨村,发生山体滑坡。造成河道堵塞2公里,100余人被掩埋。现场滑坡示意图如图2 所示,采用地基雷达对边坡连续监测,监测时间为2017年6月26日0时03分28秒至2017年6月26日4时48分28秒,历时4 h 45 min。共获取地基SAR影像20景。
凤山露天矿测试场位于北京市昌平区,现场图片如图3所示,监测时间为2018年6月9日9时54分59秒至2018年6月9日10时54分20秒,历时59 min 34 s,共获取地基SAR影像20景。
图2 茂县山体滑坡示意图
Fig.2 Schematic diagram of Mao county landslide
图3 凤山露天矿测试场
Fig.3 Feng mountain opencast mine test site
4.2 茂县滑坡数据多阈值优化实验
一般认为相干性在0.8以上则表示具有较高相干性,幅度离差在0.15以下则具有较稳定的幅度信息。故筛选参考点所用严苛阈值设置如下:幅度离差阈值 T′D=0.1,相干系数阈值T′γ=0.9,相位误差阈值 T′N=0.2。用上述阈值筛选出的稳定点分布示意图如图4所示,图中共有453个稳定点。
图4示了严苛阈值下筛选出的稳定点的示意图,从图中可以看出,筛选出的点均为具有强后向散射性的碎石,而水体、植被等均被滤除。选出稳定的参考点后,计算其形变量,作为参考形变量D′。选择以0.1为步长放宽阈值,步长太大不能有效收敛,太小增加计算复杂度。放宽阈值后稳定点邻域出现新的PS点,选择5×5大小的窗口计算窗内形变均值D,这个窗口大小与计算相干系数时的窗口大小一致。令σ=0.5×D′,若|D′-D|<σ,则继续放宽阈值,反之则停止,当前阈值即为最优阈值。σ等于参考点的形变乘以一个常数,这个常数不能太大,太大会影响精度,把非有效点引入,反之,若太小则影响PS点数量。坐标为(1524, 161)的稳定点在阈值放宽过程中PS点增加情况如图5所示。可以看到,随着阈值的放宽,选出的点从2个变为16个。
经过处理,最终确定阈值为:TD=0.45,Tγ=0.55,TN=0.2,自动阈值的最终PS点图如图6(b)所示。将用本文算法自动确定的阈值所筛选的PS点与现场光学图片进行对比,发现裸土,碎石等强散射性的点均可被选中,水体、植被等失相干严重的点可以被有效滤除,PS点图轮廓与光学图片滑坡轮廓一致,由此证明本文算法的正确性。手动调试阈值具有主观性,容易忽略不同数据之间的差异性,根据自身经验对阈值进行调整可能会造成漏选。图6(c)为手动设置阈值的结果,根据经验将阈值设置为TD=0.35,Tγ=0.65,TN=0.2,此时滑坡整体轮廓基本完整出现,植被水体等非相干像元也被滤除,然而与使用本文算法自动确定的阈值相比,滑坡顶部有明显的漏选。两种方式所筛选出的PS点数量如表1所示,手动设置阈值所筛选的点数为64157个,本文算法自动筛选的PS点数量为78205个,比人工手动设置多筛选出了14048个PS点。
图4 稳定点分布示意图
Fig.4 Stable point distribution diagram
图5 坐标为(161,1524)的稳定点及阈值放宽后其邻域的PS点
Fig.5 The stable point of the coordinate (161, 1524) and the PS point of the neighborhood after the threshold is relaxed
图6 PS点分布示意图
Fig.6 PS point distribution diagram
表1 两种方式所筛选的PS点数量(茂县)
Tab.1 Number of PS points screened in two ways(Mao county)
方式PS点数量/个手动64157自动(本文算法)78205
图7显示了茂县滑坡形变图,可以看到滑坡相对稳定,没有明显的下降趋势,大多数点的累积形状变量约在0 mm左右。 此外,滑坡顶点点形状变量最大,最大可达20 mm左右。大多数PS点未发生形变证明了此算法筛选的PS点的正确性。
图7 茂县滑坡形变图
Fig.7 Mao county landslide deformation map
4.3 凤山露天矿多阈值优化实验
使用北京昌平凤山露天矿的数据进行验证性实验,算法中新引入的变量设置情况与茂县实验时一致,自动确定阈值所筛选的PS点示意图如图9所示,自动确定的PS点筛选阈值为:Tγ=0.85,TD=0.3,与光学图片图8对比,可以看到裸土、岩石等强反射体被选为PS点,植被等失相干严重的区域被滤除,观察图10形变发现所筛选出的PS点的累积形变大多在0 mm左右,由此验证了此算法在不同场景下的适应性。
图8 凤山试验场光学图片
Fig.8 Optical picture of Feng mountain test site
图9 本文算法自动确定的阈值所筛选的PS点
Fig.9 PS point screened by the threshold determined automatically by the algorithm in this paper
图10 凤山试验场形变图
Fig.10 Feng mountain test site deformation map
5 结论
本文提出了提取地基合成孔径雷达永久散射体的多阈值优化法,选用茂县滑坡数据对算法进行了验证,得到了稳定点分布示意图、自动确定的幅度离差和相干系数的最优阈值以及根据自动确定的阈值探测到的PS点分布示意图。实验结果表明:使用稳定点及形变与其邻域形变均值之差作为判断标准来确定阈值的方法可以正确自动确定最优阈值,省去了变换场景后繁复的调制过程,避免了人工调试阈值的不合理性,推动了地基SAR技术的推广。
[1] 李媛, 孟晖, 董颖, 等. 中国地质灾害类型及其特征——基于全国县市地质灾害调查成果分析[J]. 中国地质灾害与防治学报, 2004(2): 32-37.
Li Yuan, Meng Hui, Dong Ying, et al. Main Types and characterisitics of geo-hazard in China ——Based on the results of geo-hazard survey in 290 counties[J]. Chinese Journal of Geological Hazard and Control, 2004(2): 32-37.(in Chinese)
[2] 廖明生. 滑坡变形雷达遥感监测方法与实践[M]. 北京: 科学出版社, 2017.09.
Liao Mingsheng. Method and practice of remote sensing monitoring of landslide deformation radar[M]. Beijing: Science Press, 2017.09.(in Chinese)
[3] 曲世勃, 王彦平, 谭维贤, 等. 地基SAR形变监测误差分析与实验[J]. 电子与信息学报, 2011, 33(1): 1-7.
Qu Shibo, Wang Yanping, Tan Weixian, et al. Deformation Detection Error Analysis and Experiment Using Ground Based SAR[J]. Journal of Electronics & Information Technology, 2011, 33(1): 1-7.(in Chinese)
[4] Ferretti A, Prati C, Rocca F. Permanent scatterers in SAR interferometry[J]. IEEE Transactions on Geoscience & Remote Sensing, 2001, 39(1): 8-20.
[5] Ferretti A, Prati C, Rocca F. Nonlinear subsidence rate estimation using permanent scatterers in differential SAR interferometry[J]. IEEE Transactions on Geoscience & Remote Sensing, 2000, 38(5): 2202-2212.
[6] Berardino P, Fornaro G, Lanari R, et al. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms[J]. IEEE Transactions on Geoscience & Remote Sensing, 2003, 40(11): 2375-2383.
[7] 陈强, 刘国祥, 李永树, 等. 干涉雷达永久散射体自动探测——算法与实验结果[J]. 测绘学报, 2006, 35(2): 112-117.
Chen Qiang, Liu Guoxiang, Li Yongshu, et al. Automated Detection of Permanent Scatterers in Radar Interferometry:Algorithm and Testing Results[J]. Acta Geodae Sinica Sinica, 2006, 35(2): 112-117.(in Chinese)
[8] 朱茂. 基于动态PS的地基合成孔径雷达高精度形变测量技术研究[D]. 北京理工大学, 2016.
Zhu Mao. High Precision Deformation Measurement Using Ground Based Synthetic Aperture Radar (GBSAR) Based on Dynamic Persistent Scatterer (PS) Technique[D]. Beijing Institute of Technology, 2016.(in Chinese)
[9] 胡波, 汪汉胜, 贾路路. PSInSAR永久散射体的自动探测[J]. 测绘科学, 2011, 36(2): 50-52.
Hu Bo, Wang Hansheng, Jia Lulu. An improvement method of parameters estimation in P-norm distribution[J]. Journal of Surveying and Mapping, 2011, 36(2): 50-52.(in Chinese)
[10] 黄其欢, 付五洲, 岳建平, 等. 海量时序地基SAR影像相干目标选取[J]. 数据采集与处理, 2016, 31(6): 1199-1204.
Huang Qihuan, Fu Wuzhou, Yue Jianping, et al. Selection of Coherent Targets in Large Number of Time Series Ground Based SAR Images[J]. Data Acquisition and Processing, 2016, 31(6): 1199-1204.(in Chinese)
[11] Touzi R, Lopes A, Bruniquel J, et al. Coherence estimation for SAR imagery[J]. Geoscience & Remote Sensing IEEE Transactions on, 1999, 37(1): 135-149.
