SBAS-InSAR和PS-InSAR技术在鲁西南某线性工程沿线地面沉降成因分析中的应用

张凯翔，张占荣，于宪煜

（1.中铁第四勘察设计院集团有限公司, 湖北 武汉 430063；2.湖北工业大学土木建筑与环境学院 , 湖北 武汉 430068）

0 引言

鲁西南某线性工程连接着山东省西南部经济隆起带的临、枣、济、菏发展轴，是鲁西南经济发展的重要支撑。鲁西南尤其是曲阜至菏泽区域采煤历史悠久，随着该地区以煤矿开采为主的矿产资源开发规模逐年增大，采空区上方岩层在重力作用下发生弯曲、离层乃至冒落，从而引发一系列规模不等、形状各异的采煤地面沉陷问题[1]；并且由于受到不同的矿层条件、顶板岩性特征、地质构造和采高、开采条件的影响，地表变形程度和发展过程在空间上呈现不规则分布[2]。

永久散射体干涉测量技术（persistent scatterer interferometric synthetic aperture radar, PS-InSAR）和短基线集干涉测量技术（Small baseline subsets interferometric synthetic aperture radar, SBAS-InSAR）是近年来基于合成孔径雷达干涉测量技术（interferometric synthetic aperture radar, InSAR）开展都市圈和矿区地面沉降监测分析的热门研究方向[3−6]：在都市圈应用方向，李海君等[7]采用PS-InSAR技术对覆盖河北廊坊北部地区的33景Envisat ASAR影像进行处理，并将处理得到的廊坊北部地区地表时序形变速率场，通过耦合分析时序形变速率场与沉降影像数据，分析得到该地区的地面时序沉降情况和影响因素；雷坤超等[8]采用PS-InSAR技术获取天津地面沉降时间序列演变特征，从而实现对不同含水层系统地下水开采、活动断裂等地面沉降诱发因素的分析研究；徐小波等[9]收集11期Envisat/Alos数据，采用PSInSAR技术得到安阳市区水冶镇两年间大面积地表形变动态演化信息，分析得到超采地下水与地面沉降之间的动态变化关系；张红峰等[10]提出了一种基于经验模态分解的改进PS-InSAR技术，并成功应用于非城区形变监测研究；而冉培廉、潘超、高二涛等分别对比分析了PS-InSAR与SBAS-InSAR技术处理西安地区31幅Sentinel-1A SAR影像、成都地区30景Sentinel-1A影像和南京地区23幅Sentinel-1A影像的结果，实验结果表明两种方法得到的时序性沉降变化趋势相关性极高，与实际变化情况相符，能有效应用于都市圈地面沉降监测[11−13]。在矿区沉降监测分析应用方向，陈国浒等[14]通过总结北京地区采空塌陷变形特征，针对PS-InSAR技术的特点和局限性进行可行性分析，详细比对传统形变监测手段和PS-InSAR技术的优劣，肯定了该项技术在矿区时序地表形变监测中的应用前景；卢欣奇等[15]基于PS-InSAR技术对广西平南锡基坑老采空区进行沉降监测分析，并将处理结果与数值模拟计算结论对比，证明PS-InSAR技术在采空区沉降监测分析研究中的有效性和准确性；潘光永等[16]则利用SBAS-InSAR技术对济阳井田矿区40景Sentienl-1A数据进行处理，矿区持续开采会直接影响矿区及周边地区地表形变情况，地面沉降量随开采时间的推移持续增大；李达等[17]对比了SBAS-InSAR技术和传统差分干涉技术在处理陕西某矿区13景TerraSAR-X数据的结果，实验结果表明两种方法差值在5 mm以内，SBAS-InSAR技术在矿区地表沉降监测与分析方面具有良好应用前景。

文中针对鲁西南线性工程走廊带沿线地区地面沉降开展研究，收集覆盖该地区2016年7月至2020年7月28景RadarSAT-2影像和2019年2月至2020年7月32景Sentinel-1A影像，分别采用SBAS-InSAR技术和PS-InSAR技术进行处理，最终得到沉降时间序列与沉降速率信息。然后借助GIS空间数据分析技术，结合地质、气象和矿产资料详细分析该地区地面沉降原因，并与同期地面GPS监测结果进行了可靠性验证。

1 研究区域及数据源

研究区主要包括鲁西南曲阜市、济宁市兖州区、任城区、汶上县和嘉祥县、菏泽市巨野县、郓城县、定陶区和牡丹区。研究区以京杭大运河为界，主要划分为汶泗河冲洪积平原和黄河冲积平原两大地貌单元（图1）。

图1 研究区地貌单元划分图Fig.1 Division map of geomorphic units in the study area

研究区位于华北地层大区（Ⅴ）、晋察冀豫地层区（Ⅴ4）、鲁西地层分区（Ⅴ410）的济宁地层小区，区内地层发育较为齐全，由老至新依次发育古生代寒武系、奥陶系、石炭系、二叠系；中生代侏罗系；新生代古近系以及第四系。上覆第四系地层广泛发育于平原区和东部山间洼地及河谷地带，厚度分布不均一，由东向西依次加深，厚度范围在50～400 m。沿线性工程附近5 km范围内，主要有古城煤矿、星村煤矿等多处煤矿采空区（图2）。

图2 研究区采空区分布图Fig.2 Distribution map of goaf in the study area

研究区位于华北板块（I）鲁西地块（II）、鲁中隆起区和鲁西南潜隆起区（Ⅲ级）、泗水—平邑坳陷和尼山隆起及菏泽-兖州隆起（Ⅳ级）构造单元内，主要区域断裂包括菏泽断裂、曹县断裂、巨野断裂等。其中，菏泽断裂、曹县断裂、巨野断裂属于第四纪活动断层，嘉祥断裂为早中更新世活动断层，孙氏店断裂和峄山断裂属于第四纪不活动断层，研究区整体位于菏泽——济宁断块缓慢倾斜沉降平原区。

根据覆盖研究区的星载雷达影像的存档情况，本文分别选用RadarSAT-2和Sentinel-1A数据进行研究区地面沉降监测，具体数据的基本参数信息见表1。

表1 研究区星载SAR数据源基本参数Table 1 Basic parameters of satellite SAR data sources in the study area

2 InSAR数据处理

综合考虑研究时间跨度和研究区面积，为了维持干涉图的相干性，本文实验基于GAMMA软件综合使用PS-InSAR技术和SBAS-InSAR技术实现对雷达影像的处理。SBAS-InSAR和PS-InSAR技术处理流程如图3和图4所示。

图3 SBAS-InSAR数据处理技术Fig.3 SBAS-InSAR data processing technology

图4 PS-InSAR数据处理技术Fig.4 PS-InSAR data processing technology

传统D-InSAR数据处理主要包括原始信号成像聚焦、成像质量检查、数据预处理、差分干涉计算、形变量计算等几个主要步骤。PS-InSAR数据处理则是在D-InSAR数据处理基础上再进行时间/空间域形变估算处理——PS-InSAR数据处理的核心在于PS点目标选取，以及对于PS点目标的时间/空间域形变估算。和前两种方法不同，SBAS-InSAR数据处理方法需要通过差分干涉相位和相邻点间参数相减得到残余相位，并通过对残余相位进行低通滤波、奇异值分解等处理最终得到PS点目标的时间序列形变相位。

根据数据获取时间和数据获取途径不同，本文采用SBAS-InSAR数据处理技术对Radarsat-2影像，采用PS-InSAR数据处理技术对Sentinel-1A影像进行处理。以Radarsat-2影像为例，综合考虑数据时相上的分布密度，选择2019年3月4日获取的雷达影像作为配准主影像，采用30 m SRTM DEM(Shuttle Radar Topography Mission Digital Elevation Model)数据与主影像进行地理配准，并将所有其他SAR影像配准到主影像上。为进行基线估计以确定差分像对，实验设置时间基线为1250 d，空间基线为480 m（图5）。

图5 RadarSAT-2数据基线图Fig.5 RadarSAT-2 data baseline map

3 研究结果与分析

3.1 区域时序形变分析

应用InSAR数据处理方法提取研究区2017年到2020年地面年平均沉降速率，结果见图6。数据处理结果表明研究区在2017—2020年间出现3处主要形变区域，分别是线路280 km附近古城—星村—红楼煤矿重点沉降区域P1、线路370～380 km附近郭屯—赵楼煤矿重点沉降区域P2和线路400～420 km附近菏泽市城郊工业园重点沉降区域P3。

图6 研究区2017—2020年平均沉降速率分布图Fig.6 Distribution map of the average settlement rate of the study area (2017—2020)

监测结果中P1区域位于曲阜市西北部古城—星村—红楼煤矿采矿区域内，2017—2020年最大年平均沉降速率为−136.5 mm/a，P2区域位于郓城县以南郭屯—赵楼煤矿采矿区域内，2017—2020年最大年平均沉降速率为−130 mm/a，P3区域位于菏泽市东北城郊工业园区，2017—2020年最大年平均沉降速率为−45 mm/a。如图7所示，P1和P2区域沉降中心位于线路附近的煤矿采空区范围内，沉降速率向外逐渐减小，形成沉降漏斗，中心最大沉降量达到480 mm和622 mm；而P3区域则呈现出整体沉降形变的状态，无明显沉降中心。

图7 区域沉降漏斗Fig.7 Regional settlement funnel

如图8所示，研究区主要沉降区域P1沉降中心在2017—2020年间逐渐向西移动，该区域单年累计最大沉降量保持在220 mm左右，同时该区域整体年沉降量在该时间段内呈现出先加剧后减缓的趋势，沉降区域范围缓慢扩大；主要沉降区域P2的沉降中心在2017—2020年间基本保持不变，该区域单年累计最大沉降量保持在210 mm左右，而该区域整体单年沉降量和沉降区域范围在该时间段逐步扩大；主要沉降区域P3在2017—2020年间整体单年沉降量和沉降范围基本保持不变。

图8 研究区2017—2020年逐年累计形变量分布图Fig.8 Distribution of annual cumulative deformation in the study area (2017—2020)

3.2 诱发因素分析

研究区位于鲁西南济宁至菏泽区域，由于工农业生产、民政生活、矿业开采等因素的影响，该研究区存在长期地面沉降历史。该地区地面沉降诱发因素可分为人为因素和自然因素两大类：人为因素主要包括矿业生产活动、地下水开采等，自然因素主要包括区域构造运动影响。

3.2.1 煤矿采空区塌陷

地下煤层开采后，采空区围岩应力会进行重新分布，导致围岩向临空方向出现变形；当采空区的直接顶板变形过大形成崩落后，会继续诱发上部岩层产生新的应力重分布过程；这个过程的重复发生导致掩体变形不断扩展，直至发展到地表形成地表塌陷、地裂缝或地表连续变形。因此，煤矿采空区诱发区域地面沉降的移动、变形、扩展等过程会随着采煤工作面的扩大、开采深度和范围推进而逐渐发生，并最终形成稳定的沉降漏斗；随着开采区域的转移和采空区回填措施，沉降中心会随之转移，原沉降区域沉降速率则逐渐减缓[18−20]。

①古城—星村—红楼煤矿重点沉降区域P1

根据山东省鲁南地质工程勘察院收集矿区资料显示，古城—星村—红楼煤矿重点沉降区域P1主要开采面集中在古城煤矿矿区，红楼煤矿暂属于探矿阶段，星村煤矿开采面不在研究区范围内。古城煤矿开采深度为1 000～1 200 m，煤层厚度约8 m，采用综合采放顶煤工艺条带式开采。

2017年古城煤矿3107和3209等工作面开始作业，导致以3107和3209工作面为中心的沉降漏斗初步形成，中心沉降量达到180～200 mm，边缘沉降量为50 mm左右；随着开采作业的持续，该沉降漏斗加速发展，截至2018年12月，沉降中心沉降量达到410 mm，边缘沉降量达到150 mm左右，同时影响范围逐渐扩大；为避免采空区地面沉降对线路方案施工造成影响，该两处工作面开采工作于2018年年底逐步停止并向矿区西北转移，沉降漏斗的发展也随之发生改变。截至2020年7月，沉降中心沉降量为480 mm，沉降速率自2019年3月起沉降速率大幅减缓；沉降漏斗向西发展，原沉降中心东侧沉降量保持在150 mm左右，沉降范围未扩大，原沉降中心西北侧沉降量由250 mm递增到350 mm，沉降范围也逐渐扩大（图9）。

图9 2017—2020年古城—星村—红楼煤矿重点沉降区域P1累计沉降量分布图Fig.9 Cumulative settlement distribution map of key settlement area P1, Gucheng—Xincun—Honglou (2017—2020)

②郭屯—赵楼煤矿重点沉降区域P2

根据山东省鲁南地质工程勘察院收集矿区资料显示，郭屯—赵楼煤矿重点沉降区域P2主要开采面集中在郭屯煤矿矿区，赵楼煤矿开采面不在研究区范围内。郭屯煤矿以开采山西组含煤地层为主，其中3上煤层赋存于唐官屯村东，煤层厚度1.61m，处于研究区边缘位置；3下煤层赋存于西付庄北，煤层厚度约4.77m。

截至2017年底该区域尚无明显沉降中心，唐官屯庄开采面累计沉降最大值仅150 mm，平均沉降量为100 mm；而西付庄煤矿赋存处最大沉降量为50 mm，表明2017年内该区域沉降受煤矿开采导致的采空区塌陷影响较小。截至2018年底，唐官屯庄开采面累计最大沉降值达到400 mm，该区域整体沉降速率处于加速状态，并逐渐形成以唐官屯庄为中心的沉降漏斗；同时，西付庄煤矿赋存处最大沉降量达到200 mm，整体沉降速率较2017年也有所提高，表明2018年内该区域开始受到煤矿开采导致的采空区塌陷影响。从2018年底至2020年7月，通过采集郭屯—赵楼煤矿重点沉降区域P2每隔4个月的累计沉降数据发现，唐官屯庄沉降中心沉降量以每期50 mm沉降量增长，整体沉降深度持续加深，同时沉降范围向郭屯矿区东南角不断扩大，最终形成中心沉降量达到622 mm的沉降漏斗；此外，西付庄煤矿赋存处和郭屯煤矿东南角区域中心沉降量也由200 mm增长到300 mm，形成另一处小型沉降漏斗，表明2019—2020年间，该区域受到煤矿开采导致的采空区塌陷影响日趋明显。2017—2020年间，唐官屯庄沉降漏斗边缘距线路方案最近仍有1km，而西付庄煤矿赋存处漏斗边缘最近距线路方案仅200 m，有侵入线路方案范围的迹象（图10）。

图10 2017—2020年郭屯—赵楼煤矿重点沉降区域P2累计沉降量分布图Fig.10 Cumulative settlement distribution map of key settlement area P2, Guotun-Zhaolou (2017—2020)

3.2.2 地下水开采

太沙基（Terazghi）固结理论认为：饱和多孔介质可以看作等效的连续介质，其总应力与孔隙压力的差被视为有效应力。土颗粒本身的压缩量是很微小的，在研究中可以忽略，骨架的压缩只有通过颗粒的排列变化（拓扑变化）来实现，即只有通过颗粒接触点传递的有效应力，才能引起土的变形和影响土的强度，抽水造成的饱和土孔隙水压力消散引起粒间有效应力的增加，是土层发生压缩变形的基本机理[21]。

①古城—星村—红楼煤矿重点沉降区域P1

古城—星村—红楼煤矿重点沉降区域P1范围内无集中供水水源地，部分村庄开采中深层孔隙水和浅层孔隙水作为生活饮用和农业灌溉。资料收集得到沉降中心东北约3 km的孔村中深层孔隙水水位长期观测资料，地下水位标高保持在49～50 m，在2017～2018年P1区域下沉量变化较大时，水位未出现剧烈下降的趋势，水位整体属于自然波动状态，因此该区域地面沉降基本不受地下水开采影响。

②郭屯—赵楼煤矿重点沉降区域P2

通过实地调查，郭屯—赵楼煤矿重点沉降区域P2范围内无集中供水水源地，也无深层地下水开采井，生活用水取自黄河水，农田灌溉取用浅层地下水，人类生产生活过量开采地下水引起地面沉降可能性小。

③菏泽市城郊工业园重点沉降区P3

早在20世纪80年代，菏泽城区部分水准点高程出现下降；之后经过1979、1986、2002、2014四次水准测量发现菏泽市整体出现不同程度的沉降。该区域沉降趋势呈现出东大西小，城区东侧沉降较困、城区南北两侧沉降相对较慢，沉降速率相对较快的区域主要集中城东郊区。

菏泽市城郊工业园重点沉降区P3呈现整体缓慢沉降趋势，2017—2020年间年平均沉降基本保持在30 mm/a左右，年最大沉降速率为30.82 mm/a，主要沉降区域集中在菏泽市牡丹化工园区和唐庙镇—皇镇区域（图11）。

图11 2017—2020年菏泽市城郊工业园重点沉降区P3累计沉降量分布图Fig.11 Cumulative settlement distribution map of key settlement area P3, Suburban Industrial Park, Heze City (2017—2020)

菏泽地区属于黄河冲积平原前部，区内地层主要为第四系全新统冲积层，以粉土、粉质黏土、粉砂为主，第四系地层厚度为300～500 m，局部发育有湖沼积淤泥质土，埋深35～50 m。如表2所示，该地区地下水以深层承压水为主，深层含水层与浅层孔隙水无水力联系；加之中间有中层咸水层的阻隔，补给来源较少，限制措施无法快速恢复地下水位，黏性土所失去的有效应力也难以恢复。根据菏泽市水文观测数据显示，2013—2018年间，菏泽城区地下水呈现逐年下降的趋势，深层地下水水位下降超过20 m；同时随着开采量持续增大，第四系深层孔隙水水位也逐渐下降，从而导致菏泽市城郊工业园重点沉降区P3出现整体沉降。

表2 菏泽地区第四系含水岩组水文地质特征[22]Table 2 Hydrogeological characteristics of the Quaternary water-bearing rock group in Heze[22]

3.3 外部检校与分析

为了验证InSAR数据处理得到的研究区地面沉降数据精度，将同期InSAR数据处理结果和CPI水准点测量值进行对比。文中研究选择了线路方案两侧InSAR数据处理结果覆盖范围内共59个CPI水准点与水准点所在位置的InSAR数据采样值，采用中误差、平均误差和相关系数进行验证，最终结果如表3所示。

表3 CPI测量值与InSAR数据采样值精度评估Table 3 Accuracy evaluation of CPI measurement value and InSAR data sampling value

取同期线路方案中线InSAR数据处理成果和高铁水准CPI点测量结果进行对比（图12、图13），同期曲线变化趋势基本一致。CPI点测量为剔除了部分采空区影响的测量结果，InSAR结果为处理成果，部分地区出现差异较大的情况主要为煤矿采空区范围内，差异主要为采矿造成地面变形引起，例如DK275+700处由于现场地形及施工等因素影响，未设置足够的CPI水准测量点，因此导致该段数据量较少。

图12 线路方案InSAR地面形变曲线Fig.12 InSAR deformation curve of Line plan

图13 线路方案附近CPI监测形变曲线Fig.13 CPI monitoring curve near the line plan

4 结论

文中采用时间序列InSAR技术对鲁西南线性工程穿越区域内的RadarSAT-2及Sentinel-1A数据进行处理，得到了研究区的2017—2020年年平均形变速率图和时间序列形变。研究结果表明，研究区存在3个主要沉降区域，在该时间段沉降范围及沉降量级均发生了巨大的变化。其中，沉降区域P1和P2位于煤矿矿区开采范围内，最大沉降速率分别达到136.5 mm/a和130 mm/a，累计最大沉降量为482 mm和634 mm；沉降区域P3位于城郊工业园区范围，最大沉降速率为45 mm/a，累计最大沉降量为158 mm左右。

经过论证分析得到以下结论：

（1）研究区地面沉降主要诱发因素是矿产资源开采引起的采空区塌陷和地下水超量开采，区域构造运动对该地区地表形变影响较小。

（2）沉降区域P1和P2地面沉降形变主要由矿产资源的开采引起，地下水开采对地面沉降的影响不明显；随着开采工作面的推进和转移，沉降中心会随之发生改变或移动，沉降范围逐步扩大，形成明显的沉降漏斗。

（3）沉降区域P3受深层地下水过度开采影响，由于地下水开采使用点平均分散，因此该区域未形成沉降漏斗，而是出现和地下水开采使用范围相近的沉降带。

（4）比较与InSAR数据处理成果同时期线路方案附近CPI水准点观测数据发现，两者平均误差小于1cm，相关系数达到70%以上，表明InSAR数据处理结果有效可靠，能够及时准确反映出线路方案穿行研究区内各处地面沉降变化，为线路方案规划和地质灾害整治提供有效合理参考。