利用SBAS-InSAR技术监测菏泽市地面沉降

张亚凤 刘国林 牛 冲 陈 洋 周一鸣

1 山东科技大学测绘与空间信息学院，青岛市前湾港路579号，266590

2 山东省地质测绘院，济南市二环东路11101号，250014

研究地面沉降空间分布形态和产生因素对社会稳定及灾害防控具有重要的实用价值和科学意义。合成孔径雷达差分干涉测量(D-InSAR)技术是一种新型的空间对地观测技术，具有全天候、全天时、分辨率高、监测范围广等特点[1-2]。但这种技术易受时空失相干、大气延迟和植被覆盖等影响，在地形复杂地区监测能力受限[3]。短基线集合成孔径雷达干涉测量(SBAS-InSAR)技术能有效解决时空失相干的影响，实现长时间连续监测[4]。目前已有大量学者采用SBAS-InSAR技术监测地面沉降，结果表明，该方法具有较高的准确性和可靠性[5-7]。

本文利用SBAS-InSAR技术对菏泽市2017-05-20～2021-05-23的Sentinel-1A SAR影像进行处理，获取该时段的年平均沉降速率和累积沉降量，并结合实际煤矿开采情况对研究区域累积沉降的时空演化和时序沉降场进行精细化分析。

1 研究区域概况与数据来源

1.1 研究区域概况

菏泽市地处山东省西南部，共有2个市辖区、7个县、2个开发区，面积约12 239 km2，地理范围为34°39′～35°52′N、114°45′～116°25′E。图1为研究区域的具体位置。

1.2 数据来源

Sentinel-1A SAR影像数据由欧空局提供，搭载C波段的合成孔径雷达，对地观测模式为干涉宽幅模式(IW)，幅宽为250 km，空间分辨率为5 m×20 m，极化方式为VV，重访周期为12 d，具有高重访频率、高覆盖能力以及极好的时效性和可靠性。本文选取覆盖菏泽市2017-05-20～2021-05-23期间65景Sentinel-1A SAR影像，同时引入SRTM3 DEM数据(分辨率为90 m)进行计算，以消除地形相位的影响。

2 数据处理方法

SBAS-InSAR技术是一种用于研究低分辨率、大尺度形变的时间序列分析方法。该方法首先通过设置时空基线阈值将影像自由组合生成一系列干涉图，利用最小二乘法求解地表形变信息；然后通过矩阵的奇异值分解(SVD)法进行时间序列分析，求得最小范数最小二乘法解，去除残余地形，解决时间不连续的问题；最后将多个短基线集数据联合起来进行反演，得到观测时段内的地表形变信息和平均形变速率。具体计算方法见文献[8]，本文不再赘述。

根据短基线原则，SBAS-InSAR技术将收集到的影像数据组合成若干个短基线数据集，能有效提高时空的相关性，提取到更可靠的时序沉降结果。其主要数据处理步骤如下：

1)将原始数据转换成SLC格式数据，并用实验区矢量范围裁剪影像，选取2017-07-19为超级主影像，其他影像依次与主影像配准。

2)设置时空基线阈值和多视比例值(距离向为4，方位向为1)，根据短基线原则生成143对干涉像对。利用外部DEM模拟地形相位去除干涉图中的地形相位，得到差分干涉图。

3)选用Goldstein方法进行滤波，采用基于Delaunay三角网的最小费用流法MCF对滤波增强后的差分干涉图进行相位解缠去除噪声，减少由轨道、大气和植被等因素引起的误差。

4)在地形平坦、没有相位跃变的区域选取多个控制点进行重去平，消除相位偏移量。采用SVD方法估算平均形变速率和残余地形相位。通过大气滤波估算和去除大气延迟相位，从而得到最终地表形变结果。将生成的结果编码至WGS-84坐标系下，得到LOS上的平均形变速率和累积形变量。

数据处理流程如图2所示。

3 SBAS-InSAR沉降结果与时序分析

3.1 SBAS-InSAR沉降结果

采用SBAS-InSAR技术处理65景Sentinel-1A数据，获取研究区域2017-05-20～2021-05-23期间的地面沉降信息(图3)。由图可知：

1)在监测时段内，研究区域存在一处明显的沉降盆地，沉降中心位于郓城县及周边区域，最大沉降速率为-311 mm/a(本文中负号仅表示方向)，最大累积沉降量达-1 269 mm。结合相关资料分析可知，郓城县周围煤矿众多，煤矿工作面的持续开采导致地面沉降加剧，出现沉降盆地。

2)沉降盆地周围的沉降速率为-40～-20 mm/a，沉降量为-150～-80 mm；距沉降中心较远的区域，沉降速率为-20～0 mm/a，沉降量为-80～0 mm。除郓城县这一主要沉降盆地外，巨野县西北部也有较为严重的沉降，最大年平均沉降速率为-60 mm/a，最大累积沉降量为-250 mm；鄄城县、牡丹区和曹县部分地区有轻微沉降，年平均沉降速率约为-20 mm/a，累积沉降量为-80～0 mm。

3.2 累积沉降分析

为验证SBAS-InSAR监测结果的可靠性，选取煤矿分布集中、沉降最严重的郓城区域进行详细分析。以2017-05-20的影像为基准获取郓城沉降的时空分布累积沉降图，时间间隔为0.5 a(图4)。其中，A、B、C、D为郓城区域4个主要煤矿区，监测时段内主要开采工作面如图5所示(D区因数据资料不足，未给出工作面)。

根据煤矿开采进度，将监测过程分为4个阶段进行分析。

1)阶段1：时间为2017-05-20～2018-05-15。由图4(a)、4(b)可见，煤矿A和煤矿C交界处存在明显的沉降盆地，沉降中心最大累积沉降量达-150 mm，周边沉降量普遍为-100～-20 mm。该时段内位于王楼村东南处煤矿A的工作面1开始回采，截至2018-05该工作面已回采1 530 m。SBAS-InSAR监测到该时段内此工作面附近沉降最为严重，最大累积沉降量达-300 mm。郓城其他区域也存在沉降，主要分散在各个煤矿的中心及周边城镇。位于陈河口村附近煤矿B的工作面1和2于2018年初开始开采，监测到附近最大累积沉降量为-150 mm。煤矿C前期已完成多个工作面的开采，部分地面沉降主要位于车楼村附近[9]。2018年初煤矿C的工作面1和2开始开采，导致该地区煤矿沉降最为严重，最大累积沉降量达-600 mm。煤矿D处前期开采导致地面缓慢沉降，累积沉降量和沉降范围较小，沉降量普遍为-60～0 mm。

2)阶段2：时间为2018-05-15～2019-05-10。由图4(c)、4(d)可见，随着时间推移，沉降范围逐渐增大，向北延伸至王楼村附近，向东延伸至煤矿B西侧，向西延伸至后黄岗村附近，向南延伸至煤矿C车楼村附近。沉降中心最大累积沉降量为-300 mm，周边累积沉降量集中在-250～-60 mm。2019-05煤矿A的工作面1开采结束，回采推进至2 730 m，最大累积沉降量为-600 mm。煤矿B的工作面1和2以及煤矿C的工作面1和2仍在开采中，且沉降范围逐步扩大，沉降增加，最大累积沉降量均达-800 mm。2018年底煤矿C的工作面2结束开采，2019年初工作面3接替开采。煤矿D的最大累积沉降量为-150 mm。

3)阶段3：时间为2019-05-10～2020-05-04。由图4(e)、4(f)可知，沉降范围持续扩张，南北向最为明显，王楼村和车楼村处沉降连接。沉降中心最大累积沉降量达-400 mm，周边沉降量基本为-300～80 mm。煤矿A的工作面2开采时段为2019-06～2020-04，共开采1 010 m，其附近最大累积沉降量为-800 mm。由于煤矿B的工作面1、2、3和煤矿C的工作面3、4、5持续开采(其中,2019年底煤矿C的工作面3结束开采，2020年初工作面4、5开始开采)，造成沉降范围继续扩大，部分地区最大累积沉降量超过-800 mm。煤矿D的沉降量同样持续增加，最大累积沉降量为-600 mm。

4)阶段4：时间为2020-05-04～2021-05-23。由图4(g)、4(h)可知，沉降范围往南北方向扩张更加明显，尤其在郭屯镇和王楼村附近。沉降中心最大累积沉降量达-600 mm，周边累积沉降量基本为-400～-100 mm。煤矿A的工作面3开采时段为2020-05～2021-02，共计开采650 m。之后煤矿B的工作面1、3、4以及煤矿C的工作面4、5、6接替开采。煤矿A、B、C部分区域最大累积沉降量达-1 269 mm，煤矿D最大累积沉降量为-800 mm。

综上所述，随着时间的推移，地面沉降愈加明显，无论是沉降量还是沉降范围都呈现逐渐增大的趋势。沉降主要集中在各个煤矿开采的工作面附近，最大累积沉降量可达-1 269 mm。说明郓城区域的沉降主要由煤矿开采引起，随着工作面的持续开采，沉降量逐渐增大。

为更深入地了解沉降区域的时序变化情况，选取不同煤矿工作面附近沉降较为严重的6个特征点(图4(h))，以2017-05-20的影像为基准，分析其时序累积沉降量(图6)。由图可见，随着时间的推移，P1～P6点累积沉降量逐渐增大，沉降趋势近似呈线性。其中，P1～P4位于沉降中心，监测初期沉降较小，2018年初随着煤矿开采量的增加，累积沉降量也随之增大，沉降曲线相对于其他特征点波动较多。截至2021-05-23，P1～P4的累积沉降量分别达到-939 mm、-453 mm、-820 mm和-865 mm。P5、P6位于沉降边缘地区，沉降曲线变化平缓，且沉降量较小，累积沉降量分别为-412 mm和-326 mm。综上可知，P1点的累积沉降量最大，P6点的累积沉降量最小，符合实际地质情况。

3.3 精度验证

收集郓城水准测量数据，验证SBAS-InSAR监测沉降的精度。水准观测时段为2017-05-13～2018-06-23，水准点布设情况如图7所示。选取2017-05-20～2018-06-08的SBAS-InSAR观测结果与水准测量结果进行对比，并利用线性插值法获取与SAR影像日期一致的水准测量结果，解决水准与SAR影像数据获取时间不一致的问题。

考虑到SBAS-InSAR的结果不连续，分别在郓城区域内选取沉降量较小的S1～S4点及接近沉降中心处沉降较严重的S5、S6点，对比分析不同沉降程度的SBAS-InSAR和水准测量的累积沉降结果(图8)。

由图8(a)～(d)可知，水准点S1～S4位于沉降盆地边缘地带，累积沉降量较小。由于2018-04～06期间植被生长茂密，导致SAR影像干涉质量差，影响沉降监测的精度。截至2018-06-08，4个水准点的水准测量结果分别为-35 mm、-273 mm、-200 mm和-351 mm，SBAS-InSAR监测结果分别为-45 mm、-113 mm、-55 mm和-82 mm，二者累积沉降量绝对差值分别为10 mm、160 mm、145 mm和269 mm。可以看出，SBAS-InSAR监测到的沉降趋势和沉降量与水准测量结果基本相符。

由图8(e)、8(f)可知，水准点S5、S6接近沉降盆地中心区域，累积沉降量较大。2017-05-20～2018-02-08期间S5的SBAS-InSAR监测结果与水准测量结果相差不大，且变化趋势基本一致；2017-05-20～2017-09-17期间S6的SBAS-InSAR监测结果与水准测量的沉降曲线保持着较好的一致性。随着时间的推移，水准测量结果大幅度增加，而SBAS-InSAR监测结果增加缓慢，截至2018-06-08，S5、S6的水准测量结果分别为-900 mm和-2 325 mm，SBAS-InSAR监测结果分别为-122 mm和-185 mm。由于受到植被、大气等误差因素的影响，二者结果相差较大，绝对差值分别为778 mm和2 140 mm。

根据6个点的SBAS-InSAR监测结果与水准测量结果可知，SBAS-InSAR可以监测到与水准测量相符的沉降趋势，但在沉降量上与水准测量结果存在一定的差异。在沉降量较小的区域，SBAS-InSAR与水准测量结果较为吻合；在沉降量较大的区域，由于InSAR技术本身的局限性，当单个像元内的形变超过影像的半个波长时，就无法探测出正确形变，因此InSAR在沉降中心的监测能力不足。由于SBAS-InSAR(面单元)和水准监测(点单元)的数据类型不同，且监测时段不完全一致，SBAS-InSAR监测到的沉降量与水准实测数据相差较大，说明SBAS-InSAR技术在沉降较为严重区域的精度还有待提升。

4 降雨量对SBAS-InSAR监测精度的影响

菏泽地属温带季风型大陆性气候，四季分明，春秋多风、夏热多雨，多年平均降雨量约为662.7 mm，年内降水分布不均，主要集中在6～8月。为探究降雨量对SBAS-InSAR监测精度的影响，收集菏泽市2017-05～2021-05的月均降雨量，并计算相应时段SAR干涉对的月均相干系数，分析降雨量和相干系数的关系，结果如图9所示。由图9可见，菏泽市全年降雨主要集中在6～8月，其中，8月月均降雨量最大，为199 mm；1～5月和9～12月降雨量较少，其中，12月月均降雨量最小(9 mm)。在降雨量较大的6～8月，充沛的雨量和茂盛的植被导致SAR影像质量较差，相干性较低，基本为0.37～0.40，其中，8月月均相干系数最小(0.37)；在降雨量较小的1～5月和9～12月，植被稀疏，SAR影像质量较好，相干性较高，基本为0.40～0.63，其中1月月均相干系数最大(0.63)。由此可知，SAR干涉对的相干性与降雨量有直接关系，二者呈负相关，即降雨量越小，SAR干涉对的相干性越高，监测精度就越高。

5 结 语

1)在监测时段内，郓城地区沉降较为严重，主要是由煤矿开采造成的，最大年平均沉降速率为-311 mm/a，最大累积沉降量达-1 269 mm；巨野县西北部地区也有较为严重的沉降，最大年平均沉降速率为-60 mm/a，最大累积沉降量为-250 mm；鄄城县、牡丹区和曹县沉降量较小，年平均沉降速率基本为-20 mm/a，累积沉降量基本为-80～0 mm。

2)煤矿分布集中、沉降最严重的郓城地区累积沉降量最大的区域主要在王楼村、车楼村和陈河口村附近，与实际地质情况相符，最大累积沉降量可达-1 269 mm。

3)SBAS-InSAR技术可以准确探测到沉降位置，反映大范围的沉降情况和沉降变化趋势。在沉降量较小的边缘地区，其监测结果与水准测量结果基本相符，监测精度较高；在靠近沉降中心区域，受到SAR影像特征和技术自身等限制，SBAS-InSAR监测结果与水准测量结果存在较大差异，监测精度仍有待提升。

4)根据月均降雨量和SAR影像的月均相干系数可知，SAR干涉对的相干性与降雨量有直接关系，二者呈负相关，即降雨量越小，SAR干涉对的相干性越高，监测精度就越高。