基于DinSAR技术的宁东煤炭基地区域性地表沉降监测

程霞 ，王辉 ，王国瑞 ，杜灵通 ，魏采用 ，周峰

（1. 宁夏回族自治区国土资源调查监测院，宁夏 银川 750002；2. 宁夏大学西北土地退化与生态恢复省部共建国家重点实验室培育基地，宁夏 银川 750021；3. 宁夏大学生态环境学院，宁夏 银川 750021；4. 宁夏回族自治区遥感调查院，宁夏 银川 750021；5. 宁夏回族自治区矿产地质调查院，宁夏 银川 750021）

近年来，干旱-半干旱区煤炭资源大规模开采所引起的区域地质环境恶化被广泛关注（董英等，2008；杜臻等，2023），宁东煤炭基地开发对区域植被生态的影响也引起了一些学者注意，他们利用不同手段从不同方面开展评估（王景升等，2013；杜灵通等，2018），其中煤炭开采导致的地表沉降塌陷一直备受关注（谢晓深等，2018）。王国瑞等（2018）利用地面GPS站对宁东金凤煤矿的地面变形进行监测，结果表明GPS站可以精准的监测地面点的沉降过程，特别是能捕捉到突然沉降的时间点；高冠杰等（2018）利用无人机光学遥感技术获取的多期地面高程数据对宁东羊场湾煤矿的地表沉降量进行监测，结果表明无人机光学遥感技术不仅能在小区域范围内对重点采矿区的地表沉降进行监测，还能解译提取地表变形引起的裂缝信息；张永庭等（2018a）利用无人机激光雷达对宁东马莲台煤矿重点沉降区开展监测，得出无人机激光雷达技术能很好地获取地表三维沉降数据，但受制于无人机飞行范围的限制。综上所述，目前的研究主要从单点开展宁东地表沉降监测，或在小范围开展沉降监测实验，而对于整个宁东规划区的地表沉降监测工作尚未开展，这导致政府部门无法对区域性地质环境变化情况进行准确评估。

随着合成孔径雷达（Synthetic Aperture Radar，SAR）技术的发展，特别是合成孔径雷达差分干涉测量技术（Differential InSAR，DInSAR）的出现，推动了利用卫星遥感数据开展大范围地表沉降监测的应用（李广宇等，2018；冯旻譞等，2023），为大型煤炭基地地表沉降监测提供了数据和技术手段。为此，笔者利用欧空局哨兵一号卫星（Sentinel-1A）获取的覆盖宁东煤炭基地规划区数据，监测2015～2019年该区域地表的沉降特征，以期为准确评估煤炭开发对该区域的地质环境影响提供数据支撑和理论依据。

1 研究区概况

宁东煤炭基地位于陕、甘、宁、蒙毗邻地区，地理坐标E 106°21′39″～106°56′34″、N 37°04′48″～38°17′41″，东起鸳鸯湖、马家滩、萌城矿区的深部边界，西至白芨滩东界，南起韦州矿区和萌城矿区南端的宁夏与甘肃省界，北至宁夏与内蒙古省界，东西宽约为16～41 km，南北长约为127 km，规划总面积为3 484 km²，是国家14个亿吨级煤炭基地之一。行政区划覆盖灵武市、盐池县、同心县、红寺堡开发区等4个县市（区），包括灵武市临河镇、宁东镇、马家滩镇、白土岗乡，盐池县惠安堡镇、冯记沟乡，同心县韦州镇、下马关镇，红寺堡太阳山镇等9乡镇的部分区域（图1）。基地矿区包括8个部分，位于北部的横城和南部的韦州为两个古生代石炭纪—二叠纪煤矿区，位于中部的自西向东为石沟驿、碎石井、鸳鸯湖、马家滩、积家井、萌城6个中生代侏罗纪煤矿区，矿区井田总面积为1 591.21 km²，占基地总面积的45.67%，其中以马家滩、积家井矿区面积最大。规划核心区处于灵武境内，核定面积为766 km²，约占规划区总面积的22%。

2 数据与方法

2.1 数据来源

2.1.1 Sentinel-1A干涉宽幅模式数据

2014年4月3日，欧空局发射了哨兵系列卫星之一的Sentinel-1A，该卫星使用近极地太阳同步轨道，轨道高度693 km，重访周期为12天，Sentinel-1A搭载的C波段合成孔径雷达，可以提供覆盖全球任何地方的合成孔径雷达（SAR）图像（李广宇等，2018）。Sentinel-1A可实现不同方式的极化和条带模式、超宽幅模式、干涉宽幅模式及波模式等4种工作模式，其中，以干涉宽幅模式获得的IW（Interferometric Wide Swath）影像扫描幅宽为250 km，距离向分辨率为5.0 m，方位向分辨率为20.0 m，可实现HH+HV、VV+VH、HH和VV 4种极化模式，能够实现大范围覆盖的干涉测量测量（欧阳伦曦等，2017）。笔者选取覆盖宁东煤炭基地规划区的5景Sentinel-1A IW斜距单视复数产品（SLC），下载网址为https://scihub.copernicus.eu/dhus，获取时间分别为2015年5月20日、2016年5月14日、2017年5月21日、2018年5月16日和2019年5月23日，成像时间集中在5月中下旬，植被生长等季相一致，便于地表形变干涉测量分析。

2.1.2 精密轨道星历数据

轨道信息是InSAR数据处理中非常重要的信息，从最初的图像配准到最后的形变图像生成都有着重要的作用。Sentinel-1A卫星采用精密轨道控制技术，其卫星质量控制子系统网站提供最精确的轨道数据，即POD Precise Orbit Ephemerides（POD精密轨道星历数据），该数据距离GNSS下行21天之后才可提供（李广宇等，2018），每天产生一个文件，每个文件覆盖26个小时，定位精度优于5 cm。本研究为提高形变测量的精度，使用Sentinel-1A卫星精密轨道星历数据对轨道信息进行修正，有效去除因轨道误差而引起的系统性误差（Sentinel卫星精密轨道星历数据获取网站https://qc.sentinel1.eo.esa.int/）。

2.1.3 数字地形高程模型（DEM）

笔者使用了SRTM的DEM数据，该数据由美国国家航空航天局（NASA）和国防部国家测绘局（NIMA）联合测量，SRTM覆盖地球80%以上的陆地表面。SRTM数据产品2003年开始公开发布，经历多次修订，目前的数据修订版本为V4.1版本。由于覆盖宁东煤炭基地规划区的Sentinel-lA数据跨约经度范围较大，因此下载了跨越2个经度方格的SRTM_57_05和SRTM_58_05数据，并进行了无缝拼接，数据来源于中国科学院计算机网络信息中心地理空间数据云平台（http://www.gscloud.cn）。

2.2 方法

差分干涉测量（DInSAR）是利用同一地区不同时相的SAR影像，通过差分干涉，获取该地区地表微小形变信息的技术手段，常用的有两轨法、三轨法和四轨法。DInSAR测量中要雷达波长、传感器高度、天线基线距和波束视向等的几何关系，开展解算，从而提取研究区域地面目标的微小形变信息，其基本原理如下（郑美楠等，2017）：

式中：λ为雷达波长；R1、R2为传感器到地面目标的距离；ϕf为地球曲面平地效应引起的干涉相位；ϕt为地形引起的干涉相位；ϕd为地表形变引起的干涉相位；ϕn为相位噪声；ϕa为大气对流层及电离层延迟引起的相位。

SARscape是专业雷达图像处理软件，支持Sentinel-1A数据的导入和处理，其中包括DInSAR干涉测量处理。干涉测量之前，需先设置专门针对Sentinel-1A数据做InSAR处理时的系统参数，再导入数据（包括雷达数据、精密轨道数据和DEM数据），在计算干涉相对的时空基线的基础上，开展DInSAR干涉测量处理（李旺，2018），详细步骤流程如图2所示。

图2 基于SARscape的DinSAR测量技术流程图Fig. 2 Flowchart of differential interferometry synthetic aperture radar by SARscape

3 结果与分析

3.1 相位干涉结果

选取2015年5月20日的影像作为主影像，对其余副影像相对于主影像的时空基线进行计算（表1）。其中，空间基线最长为101.11 m，最短为16.95 m，均小于临界基线，DInSAR的一个相位变化周期代表的地形变化为0.028 m。为了监测年际间的地表沉降速率和沉降特征，遥感影像选择按整年间隔，为2015～2019年每年5月中旬的Sentinel-1A数据，做干涉测量时，2幅主副影像的时间基线在360 d或372 d。

表1 Sentinel-1A干涉像对及其时空基线表Tab. 1 Temporal and spatial (perpendicular) baselines of Sentinel-1A interferograms

选择VV极化方式，做差分干涉处理，并生成干涉图。从2015年5月20日与2016年5月14日的干涉结果可以看出（图3），宁东煤炭基地规划区总体干涉信号较弱，未出现较强的干涉条纹。但去平后的干涉图开始出现干涉条纹，特别是经过滤波处理和重去平处理后，干涉条纹明显增强，在宁东北部的马莲台、羊场湾等煤矿和南部的韦州煤矿周围出现较强的干涉条纹，即在这些区域出现了因采煤沉降地表形变所引起的雷达波相位干涉。由此可见，在宁东煤炭基地规划区大范围的沉降监测中，需要对雷达影像进行深度的增强处理，以精确提取相干信息。

图3 不同处理下的干涉条纹图（2015～2016年）Fig. 3 Interferogram with different processing (from 2015 to 2016)

3.2 地表年际沉降特征

利用2015～2019年的5期Sentinel-1A数据做干涉测量，监测每年宁东煤炭基地规划区范围内的地表沉降情况。从监测结果可以看出（图4），宁东煤炭基地规划区内，不同年份间的沉降特征存在较大差异。2015～2016年，北部地区整体沉降量较大，最大沉降量达到了0.13 m；2016～2017年，整体沉降量最小，其中南部地区无沉降，北部仅在羊场湾、任家庄、梅花井、清水营等煤矿一带沉降较为明显，最大沉降量为0.10 m；2017～2018年，沉降较明显的地区也主要围绕煤矿开采活动区及其周边，最大沉降量为0.14 m；2018～2019年，宁东煤炭基地区域性普遍发生地表沉降，但沉降量较大的依然在煤矿开采活动区，最大沉降量为0.11 m。不同年份间的地表沉降差异，可能与煤矿开采过程有关，如王国瑞等（2018）监测的金凤煤矿在2015～2016年有一个快速的沉降过程，而2016年以后基本进入沉降稳定期，地表沉降开始变得缓慢。从沉降区域分布来看，宁东煤炭基地的地表沉降主要与宁东大型矿井开采活动区范围密切相关（图4），特别是煤矿开采集中的北部地区，也是沉降最为明显的区域。

图4 不同时间段宁东地表沉降及主要煤矿开采区Fig. 4 Subsidence map of Ningdong region in different periods and mainly coal mining areas

3.3 地表累积沉降特征

将2015～2019年的沉降量累加，得到宁东煤炭基地区域近4年的总沉降量（图5）。从图中可以看出，总沉降量最大的区域可达到0.37 m沉降中心在北部的清水营、灵新、梅花井、羊场湾、石槽村、麦垛山煤矿一带。为进一步分析沉降中心和沉降状况，选择了横穿宁东煤炭基地规划区东西和南北的两条剖面AB和CD，绘制沉降剖面图。从AB剖面可以看出，横穿羊场湾和梅花井矿区的剖面，在羊场湾和梅花井矿附近出现了多个超过0.1 m深的沉降中心，从剖面和沉降图均可看出，这一区域是宁东煤炭基地地表沉降最严重的区域。从CD剖面来看，宁东煤炭基地中北部沉降较为严重，且沿着CD剖面出现了4个沉降中心，到南部以后，沉降量均小于0.05 m，沉降幅度较小。

图5 2015～2019年累积沉降量及沉降切面图Fig. 5 Cumulative subsidence and its profile of Ningdong region from 2015 to 2019

利用ARCGIS空间统计功能，对宁东煤炭基地规划区的沉降面积和沉降量进行统计计算，结果表明，宁东地区2015～2019年沉降量超过0.1 m的区域面积为578 km²，为主要的沉降区；沉降量在0.05～0.1 m的沉降区有2 484 km²，面积占比最大；沉降量小于0.05 m的微小沉降变形区面积为313 km²，还有109 km²的区域未发生地表沉降或发生地表抬升。结合宁东各矿分布和煤炭开采情况，对比沉降区分布特征和沉降量，可以得出沉降量超过0.1 m的沉降区主要为煤炭开采引发的地表沉降，其他微量地表沉降可能一方面来自雷达影像的相干失真，另一方面与地表其他风沙侵蚀等变形有关。

3.4 监测精度分析

为评价DinSAR监测的宁东煤炭基地区域地表沉降的精度，本研究利用布设在宁东金凤、灵新、梅花井、清水营、任家庄和羊场湾等煤矿的26个野外GPS沉降监测点的数据，对DinSAR监测沉降量的精度进行验证，考虑到与遥感面上数据的匹配性，选择的这些沉降监测点均未在监测期内发生过短期突然沉降事件。首先根据GPS沉降监测点的经纬度坐标，在2015～2019年累积沉降量遥感监测图像中提取对应坐标像元的数值，再将其与野外GPS沉降监测点的对应实测沉降数据进行对比（图6）。结果表明，Din-SAR提取的像元地表沉降与地面实测沉降有极显著的相关性（R2=0.835 7，p＜0.01），即DinSAR监测出的沉降强度与实测一致。但从图6依然看出，实测沉降量要普遍大于DinSAR监测的沉降量，这可能与二者的监测机理不同有关，GPS监测的是地表某一点的沉降，而DinSAR获取的是对应坐标点所在区域几百平方米区域内的整体雷达反射，差分干涉测量出的结果也是一个区域面的沉降，这导致了验证点对应数据之间的差异。

图6 DinSAR监测沉降量的精度验证图Fig. 6 Precision validation of subsidence by DinSAR

笔者利用DinSAR技术和Sentinel-1A数据开展了宁东煤炭基地规划区2015～2019年的沉降监测。但卫星合成孔径雷达技术尚在发展之中，其中还存在一些技术难点和精度较低的问题，在地表形变监测中还需进一步完善，主要表现在以下2个方面：①时间失相干问题。由于两幅影像成像时间间隔较长，地面覆盖物雷达后向散射特性发生较大变化。本研究使用了近1 a的时间基线，一定程度上会造成两幅影像的雷达散射信息不相干，导致干涉条纹不连续，信噪比降低。由于Sentinel-1A卫星使用波长较短的C波段，对地物的穿透能力较弱，在春季和夏季植被茂密时期，如农田、森林等覆盖类型区（白泽朝等，2019），时间失相干会比较明显。但宁东煤炭基地规划区地处西北干旱-半干旱荒漠草原区，地表植被覆盖度低，规划区内农田、森林等土地利用类型的占比非常小（张永庭等，2018b），故在一定程度上会减弱时间失相干问题。②空间失相干。由于两幅影像成像时侧视角不同引起的失相干，即干涉对空间基线较长时，成像时刻雷达侧视角变化大，回波信号的相干性信号较低，从而导致两幅影像的失相干，往往发生在相对空间基线长度大于临界基线情况下。笔者所用的5期影像，计算出的两两相对空间基线均在临界基线以内，故最大程度的降低了空间失相干现象。

4 结论

（1）宁东不同区域的地表沉降存在较大差异，其中北部地区整体沉降量较大，年最大沉降量达到了0.13 m，南部整体沉降量较小。

（2）由于煤炭开采阶段不同，导致不同年份间的沉降量存在差异，2015～2016年沉降量较大，而2016～2017年沉降量较小。

（3）地表沉降区与煤矿开采活动区的空间分布一致，2015～2019年沉降量超过0.1 m的区域面积为578 km²，主要为煤炭开采引发的地表沉降。

（4）微量地表沉降可能一方面来自雷达影像的相干失真，另一方面与地表风沙侵蚀等变形有关。