基于SBAS-InSAR技术的矿区沉降监测研究

摘要：本文介绍了SBAS-InSAR 技术的数据处理流程，结合ALOS PALSAR-2影像数据，分析了矿区沉降形变特征。研究结果表明：矿区沉降区域空间模式主要表现为椭圆形，且随着时间推移沉降范围在逐渐扩大，沉降量级也逐渐增加，研究区域最大垂直向形变量达到了-0.20m，垂直向的最大年平均形变速率可达-0.15m/y。从空间属性上来看，最大形变峰值处位于沉降盆地几何中心，符合矿山开采沉陷理论，也较好地反映了监测结果的可靠性。

关键词：SBAS-InSAR；矿区；沉降监测

一、引言

传统的矿区沉降监测方法有GPS与水准测量，这两种传统的测量方法具有毫米级的监测精度，但整体监测成本高，维护工作量大，且监测目标为单点，只能通过单点预估出大概的沉降 [1]。合成孔径雷达干涉测量技术（Interferometric Synthetic Aperture Radar，InSAR）是用于获取地表三维信息的主动式微波遥感新兴技术[2]，它是通过获取同一地物目标的两幅合成孔径雷达（Synthetic Aperture Radar，SAR）复数图像，并将两幅 SAR 影像共轭相乘来得到包含有地表高程信息的干涉相位图。然后通过数据处理，去除干涉相位图中的高程信息来获得干涉相位图中的形变信息[4-6]。该方法能提供区域的大范围连续性形变监测，具有非接触、全天时、全天候、空间连续覆盖等优点。

本文将InSAR技术应用于矿区地表形变监测，结合矿区特点采取SBAS处理方法，并且结合实地观测资料对其检验分析，证明SBAS-InSAR技术在矿区地表沉陷监测的可靠性。

二、SBAS-InSAR的基本原理

SBAS-InSAR的数据处理流程如图1所示。①将首先构建小基线干涉对；②根据二轨差分干涉原理，对所有干涉对进行干涉、去平地效应、去地形、滤波、相位解缠等处理；③在高相干点上建立线性形变速率，高程误差的线性模型，估算线性形变速率和高程残差；④采用合适滤波算法对大气相位和地形残余相位进行去除；⑤将形变反演结果转换到合适地理坐标系下；⑥通过LOS向和垂直向的转换因子，获取垂直向的形变信息。

三、研究区域概况

从地形地貌上来看，矿区地处山脉北段，山脊线呈长条状，该山脊的延伸方向与背斜轴向基本一致。区内山脊较为平稳，两侧为迭瓦状的顺向坡，地形坡度为9°-25°，一般表现为西坡陡，东坡缓。此外，矿区内无大的地表积水体，多为流量较小的溪沟。

另外，根据多年气象观测资料，该区属于亚热带湿润型气候，气候温和潮湿，四季分明，雨量适中。区域年平均气温为17.3℃，最高气温达42.3℃，最低气温为-4.7℃。年平均总降水量为1208.3mm，年均水面蒸发量1085.0mm。

四、数据处理及结果分析

（一）数据来源

本实验数据为2018年10月18日至2020年5月28日的5景ALOS PALSAR-2影像（其中2019年9月5日的影像为SM1成像模式，其余四景为SM3成像模式）。

（二）数据处理

1.根据小基线原则生成SAR数据干涉对和时空基线连接图

本實验选取2018年10月18日的影像为主影像，通过设置时空基线阈值，共生成6个小基线干涉对。干涉对的具体信息见表2。

2.差分干涉处理

干涉对组合确定后，根据二轨差分干涉原理，对所有干涉对进行干涉、去平地效应、去地形、滤波、相位解缠等处理。

由于本实验所选用的ALOS PALSAR-2扫描模式（SM1）幅宽为50KM，且研究区域地形多起伏，故在干涉时采用10：10多视，以牺牲分辨率为代价提高信噪比及处理效率。考虑到时效性，选用TanDEM-X 90m DEM数据模拟地形相位。根据计算，影像覆盖区域最大高差约 480米。

在获得地形模拟相位后，即可从干涉图中扣除地形相位以获取地表形变相位（实际上还混叠有其他误差）。然而由于各种失相干作用导致差分干涉图中存有不同程度的噪声，影响之后相位解缠质量和沉降监测精度，此处采用Goldstein 进行滤波。滤波后的影像能够有效削弱干涉相位中的绝大部分噪声，从而为后续相位解缠提供更高质量的基础数据。在解缠处理中，本实验选用最小费用流解缠分别对各差分干涉图进行解缠。

3.线性形变建模与估计

对所有构建的干涉对实现差分处理后，依据相干系数选择时序稳定的高相干目标。其后在高相干点上建立线性形变速率，高程误差的线性模型，然后在线性模型的基础上构建方程组，利用矩阵的奇异值分解（SVD）法，求解方程组并估算线性形变速率和高程残差。

4.残余相位估算和误差去除

通过形变速率去除后，残余相位中主要包含了噪声和大气延迟湍流部分。为此，本实验采用高通滤波和低通滤波削弱大气延迟，从而提高时序形变精度。其中大气低通滤波算法和大气的空间特性相关，由于大气在空间上具有高相干性，采用空间窗口进行滤波时，窗口越小，其滤波效果越强；大气高通滤波算法和大气的时间特性有关，由于大气在时间上具有低相干性，采用时间窗口进行滤波时，窗口越大，其滤波效果越强。因此，本实验选择空间域的低通滤波和时间域的高通滤波算法对大气相位和地形残余相位进行去除。

5.地理编码

上述过程均是在SAR坐标系下进行的，地理编码可以将形变反演结果转换到地理坐标系下，以便和实测结果进行对比和分析。图2、图3为2018年10月18日至2020年5月28日矿区的LOS向的累计形变图以及形变速率图。

6.垂直向形变分解

InSAR技术获取的是雷达视线方向的形变，该方向的形变往往不能准确反映矿区的形变特征，因此，可以通过LOS向和垂直向的转换因子（入射角的余弦值 cosθ），获取垂直向的形变信息。为顾及InSAR成像几何结构对沉降监测的影响，此处忽略水平移动将LOS向形变直接转化至垂直向，形变结果如图4、图5所示。

（三）结果分析

图2中白色实线标注的区域为矿区所在区域。结果清晰地显示了20余个地表沉降区域（红色区域），其空间模式主要表现为椭圆形，且随着时间推移沉降范围在逐渐扩大，沉降量级也逐渐增加，最大LOS向形变量级达到了-0.16m。可以推测，这些地表沉降很可能为地下开采所致，因为其时空演化模式与地下开采形变演化模式非常相似。

图3为矿区及周边地区2018年10月18日至2020年5月28日期间的LOS向年平均形变速率图。从中可以看出，形变速率的空间分布与形变累积分布情况基本相同，LOS向最大年平均形变速率可达-0.12m/y。

图4为垂直向形变图。由图可知，主要形变区域的空间分布与LOS 向形变区域分布基本相同，垂直向的最大累积形变量可达-0.20m。

图5为矿区及周边地区2018年10月18日至2020年5月28日期间的垂直向年平均形变速率图。由图可知，垂直向的最大年平均形变速率可达-0.15m/y，数值高于LOS向的形变速率。

为直观分析矿区开采引起地表形变时序特征，本实验在矿区开采形成的地表沉陷盆地上方绘制剖面线D-D（图6所示），随后提取各期沉陷量并绘制剖面图，具体见图7。

如图7所示，（a）和（b）分别表示时序LOS向累积形变及时序垂直向累积形变，从时间属性上来看，自2018年10月18日开始，剖面累积形变随时间推移逐渐加大，在LOS 向上形变峰值可达-0.10m，转化至垂直向后最大峰值可达-0.12m，并且对照（a）和（b）可知垂直向累积形变趋势与LOS向累积形变趋势一致。此外，2019年9月5日至2020年5月28日期间，剖面累积形变增加显著（LOS向最大形变由0.06m增至0.10m）。

从空间属性上来看，结合累积形变与距离的关系可知最大形变峰值处大致位于沉降盆地几何中心，较为符合矿山开采沉陷理论。此外，对照时序InSAR获取的全盆地沉降演化图可知剖面线横跨两个地表移动盆地，在圖8中体现出了两个沉降峰值，也较好地反映了监测结果的可靠性。

五、结束语

本文基于ALOS PALSAR-2数据，利用SBAS-InSAR技术，得到某矿区2018-2020年各时间段的地表形变情况。结果显示，研究区域存在20余个地表沉降区域，其空间模式主要表现为椭圆形，且随着时间推移沉降范围在逐渐扩大，沉降量级也逐渐增加，垂直向的最大累积形变量可达-0.20m，垂直向的最大年平均形变速率可达-0.15m/y。监测结果符合矿山开采沉陷理论并与实际采空区位置基本吻合，证明监测结果具有一定的有效性。

作者单位：况长虹    四川 成都

参  考  文  献

[1] 刘冬，刘星年，陈超.基于CR-InSAR的煤矿区地表沉陷监测研究[J].现代测绘，2015，38（4）：26-30.

[2] 杨成生.基于D-InSAR技术的煤矿沉陷监测[D].长安大学， 2008.

[3] 刘永坦.雷达成像技术[M]. 哈尔滨工业大学出版社， 2014.

[4] 王远坚，齐麟，姜岳.基于D-InSAR技术的金矿老采空区下沉监测与位置反演[J]. 金属矿山， 2019（10）：68-73.

[5] 范军林，刘飞鹏，涂梨平.基于D-InSAR技术的南昌市地表形变监测研究[J].矿产与地质，2019，33（3）：546-550.

[6] 陈优良，王之玉，栗正刚，等.基于D-InSAR的矿山地表沉陷监测研究[J].化工矿物与加工，2019，48（8）：10-14.