2020, Vol. 41
2. 中国黄金集团, 北京 100000;
3. 辽宁科技大学 土木工程学院, 辽宁 鞍山 114053
2. China Gold Group, Beijing 100000, China;
3. School of Civil Engineering, Liaoning University of Science and Technology, Anshan 114053, China
近年来, 我国各类地质灾害事件时有发生, 如何为易发生地质灾害区域的实时预警提供准确形变监测数据已成为一个亟待解决的问题.地基合成孔径雷达(GB-SAR)是集成了步进连续波、合成孔径雷达干涉测量、PS-InSAR等技术为一体的新型地面形变监测设备.由于GB-SAR具有全天候、高精度等形变监测的技术优势, 已被广泛应用于露天矿边坡、大坝、桥梁、公路、冰川等形变监测领域[1-2].
雷达数据处理结果的准确性是露天矿边坡滑坡危害预警的重要保障.在雷达数据处理的研究中, Iglésias等[3]针对陡峭山区的雷达数据大气相位误差改正提出了多元回归模型且模型高度相关, 但该模型所用气象数据源是沿山路近地表环境获取, 不适合露天矿坑环境; 曲世勃等[4]分析了GB-SAR形变监测误差来源并验证了GB-SAR形变监测数据精度; 徐亚明等[5]利用永久散射体改正网对气象误差进行改正.虽然上述方法对GB-SAR形变监测误差都有不同程度改正效果, 但均忽略了方位向的轨道误差[6]对GB-SAR形变监测精度的影响.本文以马兰庄露天铁矿GB-SAR原始影像为数据源, 基于小周期处理方法首先对雷达原始影像进行干涉处理生成干涉图, 然后对干涉图进行小波去噪处理, 并采用三重阈值算法提取了雷达影像中的PS点和高质量PS点.其中利用PS点进行形变分析, 而高质量PS点由于其在时间序列上所表现出的高稳定性特征(理论干涉相位为零), 因此可认为其干涉相位是由于大气误差相位和轨道误差相位共同影响所致.本文通过对马兰庄GB-SAR形变监测误差的规律分析, 提出了基于高质量PS点的相位与像元坐标建立多元回归模型进行相位误差改正的方法, 结果表明所建模型可准确改正GB-SAR形变监测误差, 为边坡的稳定性分析提供了准确形变数据源.
1 GB-SAR原始影像采集及处理流程本次形变监测使用型号为LKR-05-KU国产地基合成孔径雷达, 其方位分辨率在1 km处为4 m, 距离分辨率为0.3 m.地基合成孔径雷达被安置在露天采场东坡和西坡之间的稳定区域, 如图 1所示, 并对采场东坡进行时长为60 h的形变监测, 数据采集周期设置为12 min, 共采集了300幅原始影像.本文采用小周期处理算法, 其处理流程如图 2所示.
|
图 1 GB-SAR边坡监测 Fig.1 GB-SAR slope monitoring (a)—被监测边坡; (b)—GB-SAR监测设备. |
|
图 2 数据处理流程 Fig.2 Data processing flow |
由于本实验GB-SAR观测距离为250~1 200 m, 且在观测过程中GB-SAR安置位置及扫描轨道参数均未发生改变, 因此用于干涉计算的两幅原始影像数据有严格的配准精度, 不存在空间基线解算等问题.在干涉计算过程中只需将两幅原始影像矩阵共轭相乘, 即为两幅影像的干涉图, 设回波矩阵中任意信号为I1, 对应相邻影像中相同位置的回波信号为I2, 则干涉信号I12为
|
(1) |
干涉相位Δφ为干涉信号I12对应的相位角, 其包含雷达视线形变相位φlos、大气误差相位φatm、频率偏移误差相位φf、轨道误差相位φtrack及噪声误差相位φnoise.故干涉相位Δφ可表示为
|
(2) |
由于雷达天线工作频率比较稳定, 且在1 km处由于频率偏移引起的误差量级为10-4 mm, 故频率偏移误差相位φf可以不予考虑.
噪声误差相位主要体现在干涉图中的斑点噪声.由小波去噪原理可知, 噪声在频域内多出现在图像的高频处, 图像的纹理特征多出现在图像的低频处.因此本文利用小波去噪算法[7]对干涉相位进行了三层小波分解, 并基于Stein的无偏似然估计自适应阈值算法对干涉相位进行滤波处理.通过处理前后效果对比, 该去噪方法可有效剔除干涉图中的噪声误差相位φnoise, 且能保持干涉图的纹理特征.因此, 本研究所用干涉图都经过小波去噪处理.
经过小波去噪处理后的GB-SAR形变监测误差主要来自于大气误差相位和轨道误差相位.因此文中将对上述两项误差的剔除方法进行论述.
2.2 GB-SAR数据PS点选取由于影像中的部分像元不是边坡稳定散射体的反射信号, 故需要提取稳定PS点来反映真实的边坡形变量[8-9].此外, 为了准确建立误差改正模型, 还需要提取少量高稳定性的强散射点, 称之为高质量PS点.本文利用雷达原始影像计算了平均幅度、幅度离差指数和平均相干系数三个参数, 提出了基于三重阈值算法来筛选特定PS点的方法.
平均幅度M是幅度在整个时间序列的均值, 其值越大, 反射信号越强, 其计算公式为
|
(3) |
|
(4) |
式中:Mt为单幅影像幅度均值; m和n分别为影像方位向和距离向尺寸; A(i, j)为边坡反射信号; k为影像数量.
幅度离差指数D是衡量一个像元内边坡在时间序列内对雷达信号反射强度的稳定性, 其值越小代表越稳定.计算公式为
|
(5) |
式中:σi, j为像元的幅度标准差; Mi, j为像元的幅度均值.
平均相干系数L是雷达数据中每个像元在整个时间序列内相干系数的均值, 其值越接近1, 相干性越好, 其计算公式为
|
(6) |
式中, It, 1和It, 2为像元复信号; I-t, 2为It, 2的复共轭.
平均幅度、幅度离差指数和平均相干系数的计算结果如图 3所示, 其中平均幅度在雷达图像上用雷达反射强度表示.通过以下步骤选取PS点与高质量PS点.
|
图 3 GB-SAR数据参数 Fig.3 GB-SAR data parameters (a)—雷达反射强度; (b)—幅度离差; (c)—相干系数. |
步骤1 设置幅度均值阈值M1=1.5M, 并提取像元平均幅度Mi, j>M1的像元作为初始PS点, 再提取像元平均幅度Mi, j>4M1的像元作为初始高质量PS点.
步骤2 设置提取PS点的幅度离差阈值D1=0.5, 提取高质量PS点的阈值D2=0.1, 并基于初始PS点提取幅度离差Di, j＜D1的像元作为过渡PS点, 基于初始高质量PS点提取幅度离差Di, j＜D2的像元作为过渡高质量PS点.
步骤3 设置相干系数低阈值L1=0.85, 高阈值L2=0.98, 并基于过渡PS点提取平均相干系数Li, j>L1的像元作为最终PS点, 基于过渡高质量PS点提取平均相干系数Li, j>L2的像元作为最终高质量PS点.
2.3 大气误差相位和轨道误差相位改正由于微波中的Ku波段在传播过程中对空气介质的气象数据比较敏感, 且空气介质在两次采集的时段内会发生一定变化, 故干涉相位经过滤波后, 还包含一定的大气延迟误差相位[10-14].由于雷达每次的运行状态和轨道会存在一定的偏差, 故干涉相位中还包含轨道误差相位.其中轨道误差长期存在, 大气延迟误差只在大气环境变化剧烈的情况下影响较大, 上述两种误差是干涉图误差的主要来源.干涉图误差相位分布如图 4所示.
|
图 4 干涉图误差相位 Fig.4 Interferogram error phase |
假设微波在露天矿坑中传播的空气介质均匀, 设主影像的平均大气延迟系数为α1, 副影像的平均大气延迟系数为α2, 雷达天线距离PS点的观测距离为Ld, 则两次观测的大气延迟误差φatm可表示为
|
(7) |
由于雷达距离向分辨率为0.3 m, 故每个像元在距离向对应的长度为0.3 m.设雷达起始采样距离为Ls, 故Ld可以用距离向的像元坐标x表示为
|
(8) |
由雷达影像的方位向分布特征可知, 轨道误差与方位向像元坐标y呈较好的线性分布, 因此轨道误差φtrack可表示为
|
(9) |
故滤波后的GB-SAR形变监测误差相位φGB-SAR与像元坐标x, y可建模为
|
(10) |
式中:α为距离向改正系数; β为方位向改正系数; ζ为回归常数.
利用上述模型对干涉图误差相位进行改正, 得到真实的形变相位, 如图 5所示.
|
图 5 干涉图改正后相位 Fig.5 Corrected interferogram phase |
为了验证模型的准确性, 本文针对具有较大误差的6幅干涉图进行了分析, 分别计算了影像改正前后的平均相位误差, 以及基于多元回归模型的拟合优度, 如表 1所示.
|
表 1 影像误差及拟合优度统计表 Table 1 Statistics of image error and fit goodness |
多元回归模型中自变量和因变量的拟合优度越高, 则误差修正精度越高.由表 1可以看出, 6幅影像中像元坐标和误差相位的拟合优度均大于0.9, 说明两者高度相关.经过改正后平均相位误差明显变小, 平均相位误差降低了83.3%.
2.4 边坡形变计算及分析由于实验数据采集周期为12 min, 因此可认为一个小周期内的边坡形变量小于λ/4, 即4.5 mm, 故可认为处理后的干涉相位即为真实形变相位.根据雷达参数可以推导出相位差和形变量之间的关系[15]为
|
(11) |
式中:ΔL为雷达视线方向的形变量; φlos为雷达视线上的形变相位; c为光速; f为雷达信号发射频率.
本文基于Matlab编写了GB-SAR原始数据处理程序, 实现了GB-SAR原始数据的自动处理.计算出了每个小周期内边坡在雷达视线方向上的形变量, 并对形变图中空缺的像元(距离PS点小于6 m)进行插值, 最后通过在时间序列上叠加计算分析了10, 20, 30, 40, 50, 60 h共6个阶段的形变结果, 如图 6所示.
|
图 6 各时间段形变图 Fig.6 Deformation map of each time period (a)—10 h; (b)—20 h; (c)—30 h; (d)—40 h; (e)—50 h; (f)—60 h. |
由图 6可以看出, 在60 h内边坡逐渐出现了一个明显的形变区域.通过对图中形变区域对应的形变数据进行统计分析, 该区域内最大形变量为-53.243 mm, 形变区域大于10 mm的点共3 099个像元, 形变面积约1 860 m2, 平均形变量为-21.806 mm, 平均形变速率为-0.36 mm/h.
通过对GB-SAR形变监测数据与马兰庄露天矿DEM数据进行三维配准, 确定了形变区域位置, 经过现场勘查发现在图 6f蓝色区域对应的实际位置出现了一个宽度为7 cm左右的裂缝, 如图 7所示.由此可见利用本文提出的方法处理露天矿GB-SAR监测数据可以为边坡的稳定性分析提供准确形变数据源.
|
图 7 边坡裂缝 Fig.7 The slope crack |
1) 基于三重阈值算法可以准确提取用于形变分析和误差改正的PS点和高质量PS点.
2) 露天矿GB-SAR形变监测数据的误差主要来源于大气延迟误差和轨道误差.基于多元回归模型的误差改正方法可以准确改正露天矿GB-SAR形变监测的误差相位, 依据本文实验数据可使GB-SAR形变监测误差降低83.3%.
3) 基于本文对露天矿GB-SAR数据的小周期处理方法可以准确提取形变区域及形变量, 为边坡的稳定性分析提供准确形变数据源.
| [1] |
Noferini L, Pieraccini M, Mecatti D, et al. Using GB-SAR technique to monitor slow moving landslide[J]. Engineering Geology, 2007, 95(3): 88-98. |
| [2] |
Serrano-Juan A, Vázquez-Suè E, Monserrat O, et al. GB-SAR interferometry displacement measurements during dewatering in construction works:case of La Sagrera railway station in Barcelona, Spain[J]. Engineering Geology, 2016, 205: 104-115. DOI:10.1016/j.enggeo.2016.02.014 |
| [3] |
Iglésias R, Rubén X, Fàbregas A, et al. Atmospheric phase screen compensation in ground-based SAR with a multiple-regression model over mountainous regions[J]. IEEE Transactions on Geoscience & Remote Sensing, 2014, 52(5): 2436-2449. |
| [4] |
曲世勃, 王彦平, 谭维贤, 等. 地基SAR变形监测误差分析与实验[J]. 电子与信息学报, 2011, 33(1): 1-7. (Qu Shi-bo, Wang Yan-ping, Tan Wei-xian, et al. Deformation detection error analysis and experiment using ground based SAR[J]. Journal of Electronics and Information Technology, 2011, 33(1): 1-7.) |
| [5] |
徐亚明, 周校, 王鹏, 等. GB-SAR构建永久散射体网改正气象扰动方法[J]. 武汉大学学报(信息科学版), 2016, 41(8): 1007-1012. (Xu Ya-ming, Zhou Xiao, Wang Peng, et al. A method of constructing permanent scatterers network to correct the meteorological disturbance by GB-SAR[J]. Journal of Wuhan University(Information Science), 2016, 41(8): 1007-1012.) |
| [6] |
Peng J, Cai J, Yang H. A rail central displacement method about GB-SAR[J]. International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 2016, 51(7): 865-871. |
| [7] |
杜立志, 殷琨, 张晓培, 等. 基于最优小波包基的降噪方法及其应用[J]. 工程地球物理学报, 2008, 5(1): 25-29. (Du Li-zhi, Yin Kun, Zhang Xiao-pei, et al. The signal denoising algorithm based on optimization wavelet packet and its application[J]. Chinese Journal of Engineering Geophysics, 2008, 5(1): 25-29. DOI:10.3969/j.issn.1672-7940.2008.01.006) |
| [8] |
Monserrat O, Crosetto M, Luzi G. A review of ground-based SAR interferometry for deformation measurement[J]. ISPRS Journal of Photogrammetry & Remote Sensing, 2014, 93(7): 40-48. |
| [9] |
Iglesias R, Mallorqui J J, Lopezdekker P. D-InSAR pixel selection based on sublook spectral correlation along time[J]. IEEE Transactions on Geoscience & Remote Sensing, 2014, 52(7): 3788-3799. |
| [10] |
Iannini L, Guarnieri A M. Atmospheric phase screen in ground-based radar:statistics and compensation[J]. IEEE Geoscience & Remote Sensing Letters, 2011, 8(3): 537-541. |
| [11] |
Noferini L, Pieraccini M, Mecatti D, et al. Permanent scatterers analysis for atmospheric correction in ground-based SAR interferometry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2005, 43(7): 1459-1471. DOI:10.1109/TGRS.2005.848707 |
| [12] |
Huang Z, Sun J, Li Q, et al. Time- and space-varying atmospheric phase correction in discontinuous ground-based synthetic aperture radar deformation monitoring[J]. Sensors, 2018, 18(3): 1-14. |
| [13] |
Rödelsperger S, Becker M, Gerstenecker C, et al. Digital elevation model with the ground-based SAR IBIS-L as basis for volcanic deformation monitoring[J]. Journal of Geodynamics, 2010, 49(3/4): 241-246. |
| [14] |
Xing C, Huang J J, Han X Q. Research on the environmental effects of GB-SAR for dam monitoring[J]. Advanced Materials Research, 2014, 919/920/921: 392-397. |
| [15] |
Hu C, Deng Y, Tian W, et al. A PS processing framework for long-term and real-time GB-SAR monitoring[J]. International Journal of Remote Sensing, 2019, 40(16): 1-17. |