InSAR技术在高寒矿区地灾治理监测中的应用

吴启红 谢飞鸿 杨何 董建辉 魏占玺 毋远召 黄程

摘 要：采用InSAR技术对大通煤矿沉陷区进行监测，得到了该矿区2014年10月至2021年11月的地面沉降变形情况和年均沉降速率等信息.结果表明，矿区主要沉降变形区集中在东北部和东南部2个片区，其沉降速率最大超过220 mm/yr；2017年后矿区沉降变形范围和变形量值都呈显著下降趋势且沉降变形逐渐平稳，前期治理工程对控制矿区沉降变形起到了重要作用，可为高寒矿区防灾减灾及治理工程评价提供重要参考.

关键词：InSAR；矿区治理；地面沉降；沉降速率

中图分类号：P694

文献标志码：A

0 引 言

矿区地质灾害主要包括崩塌、滑坡、泥石流、地面塌陷、地裂缝、矿井突水和矸石堆灾害等；矿区环境灾害主要包括水污染、土污染、水土流失、石漠沙化、地貌景观破坏和植被破坏等.据不完全统计，在1970—2005年，各类矿山地质灾害共计达到2 031次，而在2001—2010年仅10年时间就发生了矿山次生灾害12 367起，其中尤以地面塌陷、滑坡和崩塌居多［1］.矿区内剥离大面积土壤植被及形成的煤矸石山，产生大量“黑色”污染和“灰色”污染，形成潜在不稳定斜坡，对矿区的生产设施和施工人员造成严重威胁.而高寒矿区由于高海拔、低氧和低温，其地质环境尤为复杂及脆弱，如玉龙煤矿、木里煤矿和红河沟煤矿等.因此，开展高寒矿区的变形监测对保障矿区安全至关重要.

传统的矿区沉降监测方法主要采用精密水准测量和全球定位系统（GPS）测量等.这些监测方法存在监测工作量大、费时、费财和测点难以保存等缺陷，同时其测量精度并不高，从而导致从其测量数据获得的结果会出现一定的偏差，对后续工作造成影响［2］.合成孔径雷达干涉测量（InSAR）因在地表变形监测上具有全天时、全天候、高分辨率、高覆盖率、直观和形象等其他测量技术不具备的优点，可以在部分领域代替传统的测量方式［3-4］.同时为去除大气误差的影响，加之合成孔径雷达（SAR）影像数据的增加，在原有InSAR方法的基础上，逐渐发展出Stacking-InSAR、D-InSAR、PS-InSAR和SBAS-InSAR等高级时序InSAR处理方法［5-6］.

InSAR技术最先被用于地球表面大范围形变场的观测，如地震位移测量和冰川漂移等.随着理论发展与数据处理技术的进步，InSAR技术逐渐在国内外各行各业中得到广泛应用，如煤田、石油、天然气开产矿区以及地下水开采区地面沉降监测［7］.高海英等［8］通过小基线集InSAR技术识别黔西南州贞丰县和安龙矿山地表形变；丁刘建等［9］基于SBAS-InSAR对滕州市附近矿区地面沉降进行监测；王志勇等［10］基于InSAR对济宁矿区沉降进行精细化监测.这些研究主要针对矿区生产运营阶段进行变形监测，而对矿区治理过程及后期的变形监测及应用研究不足.InSAR技术无论是在矿区生产运营阶段、闭矿阶段还是在治理过程中对灾害的预防及评价方面都具有重要意义，同时对于指导我国在类似地区的防灾减灾工作也具有重要参考价值.因此，为了研究InSAR技术在高寒矿区治理中应用的有效性和实用性，本研究选取青海省大通煤矿区作为研究对象，采用InSAR技术对大通沉陷区进行处理，得到了该矿区治理过程及治理后的地面沉降变形情况和年均沉降速率等信息，分析变形规律及评价治理工程的有效性.

1 研究区概况及数据

1.1 研究区位置

研究区位于青海省西宁市大通回族土族自治县，地处青藏高原和黄土高原的过渡地带，地理位置东经101°37′～101°56′，北纬36°55′～36°58′之间，海拔2 450～2 850 m，距省会西宁市35 km.

1.2 治理工程情况

2012—2015年共分4期对沉陷区内涉及桥头镇及良教乡的7个村（桥头镇的元树尔、小煤洞和大煤洞3个村，以及良教乡的下甘沟、煤洞沟、白崖和上甘沟4个村）范围内因煤炭开采造成地表土地资源压占、挖损、植被破坏、粉尘污染和地下含水层破坏等，以及后期剩余煤炭资源开采可能引发的采空塌陷范围进行治理，治理建设规模为1 385.02 hm2.虽然经过4期地质环境治理，但区内依然存在地质灾害和生态环境破坏等问题，因此，2020—2022年，又对历史遗留和新产生的塌陷区进行了综合整治.

1.3 数据情况

1.3.1 Sentinel-1数据及获取

Sentinel-1卫星数据是欧洲航天局为接替Envisat-ASAR卫星发射的科研SAR卫星，数据免费开放.该卫星成像质量高，轨道稳定，有多极化能力，单颗卫星重访周期24 d，继Sentinel-1A卫星后，Sentinel-1B卫星的发射使常用重访周期变为12 d，遇到突发情况时可以2颗卫星合作将观测重访周期缩短至6 d.查询可知，该地区Sentinel-1数据覆盖良好，数据日期自2014年10月至今.Sentinel-1数据参数见表1.

Sentinel-1卫星执行全球对地观测极化，以干涉宽幅成像为主，该模式有益于InSAR的干涉处理，且250 km×250 km的幅宽保证单景数据覆盖观测区.该卫星工作波段为C波段，虽然对植被的穿透能力较弱，但卫星数据量多、成像質量高，可以一定程度上弥补数据穿透能力弱的问题.本次计算使用2014年10月7日至2021年11月23日间共计173期降轨Sentinel-1数据，相对轨道号33，数据重访周期为12/24 d.

1.3.2 DEM数据

本次计算使用ALOS World 3D DEM，其空间分辨率为30 m，该数据集是高精度全球数字地表模型数据，水平分辨率为 30 m，高程精度5 m，由高级陆地观测卫星ALOS上搭载的全色遥感立体测绘仪PRISM获取.是目前世界上最精确的3D地图之一，覆盖全球所有的土地尺度.

2 InSAR数据处理

2.1 InSAR技术原理

InSAR技术主要基于雷达影像的相位信息获取目标点至雷达传感器之间的距离，通过获取2期甚至多期次影像并对其进行干涉即可丈量目标点在一定时间段内的位移变化量，其属于定量遥感技术范畴［11-12］.目前，常用星载SAR平台获取的SAR影像进行差分干涉处理，其具有覆盖面积大、飞行轨道稳定和获取数据成本较低等优势.星载SAR平台的基本工作模式为重轨飞行，即1颗雷达卫星在相邻轨道上对同一地区2次SAR成像，每幅SAR图像中的每个像素均记录了强度值和绝对相位值.针对同一地区获取的2幅SAR影像进行差分干涉处理，可得到该地区地表形变信息.

2.2 数据处理

2.2.1 SAR数据配准

在进行干涉计算前，首先要对生成的单视复数数据（SLC）在外部DEM數据的参与下进行数据配准.Sentinel-1数据成像模式较为特殊，为多个子带拼接构成的，为保证子带间干涉处理时不发生相位的跃变，SLC数据配准精度要求较高，尤其是针对Sentinel-1数据而言，距离向（Azimuth）配准精度需要达到1/1 000像元级别.

2.2.2 差分干涉处理

如图1所示，选择2018年1月19日SAR影像为主影像进行图像配准计算.其中，Azimuth方向配准精度小于1/1 000像元，满足Sentinel-1数据配准要求.干涉像对最小时间间隔为12 d，最大时间间隔为120 d，最小垂直基线长度为0.3 m，最大垂直基线长度为216.8 m，完全满足干涉计算需求，共有干涉像对513个.

2.2.3 自适应滤波处理

为平滑干涉图噪声，使用窗口为32的自适应滤波进行处理，增强干涉图低质量区的相干性，方便下一步相位解缠处理.

2.2.4 相位解缠

使用最小费用流方法进行相位解缠处理，相干性阈值设置为0.2，解缠相位均过度平滑，未出现相位跳跃的情况，但部分解缠相位图受大气效应干扰严重.

2.2.5 大气相位估计

估算高程相关的大气延迟相位，并利用快速空间滤波的方式初步估算随机大气相位，而后去除大气延迟相位，大气相位干扰情况将被明显改善.

2.2.6 Stacking-InSAR和SBAS-InSAR计算

利用去除大气相位后的计算结果进行2014—2021年度的Stacking-InSAR计算，估算年均变形速率，并分析工作区逐年沉降变形情况.同时通过SBAS-InSAR方法计算典型形变点的累计变形曲线.

2.3 误差源分析

InSAR计算中，会受到多种因素的影响，而且在不同地区，使用不同类型数据，同样的影响因素造成的误差程度不同，在计算过程中要结合实际情况加以分析，针对不同地区和不同数据来调节参数.误差源主要有系统误差、DEM误差、失相干和大气误差.

3 InSAR结果分析

受制于Sentinel-1卫星数据本身分辨率，工作区InSAR观测结果空间分辨率为近15 m.工作区观测到的变形主要是采矿后沉降变形，所以将InSAR计算得到的雷达视线向（LOS）变形转为垂直向变形.

3.1 2014—2021年逐年变化分析

利用Stacking-InSAR方法计算工作区2014—2021年间逐年地面沉降变形情况，2014年的InSAR数据是从2014年10月开始，经过运算后，图中代表的是整个自然年的变形分布，如图2所示.

由图2可知，大通矿区地面沉降情况自2014—2021年呈现在变形范围和沉降量值上均表现为减弱的趋势.矿区的沉降区主要集中在东北部和东南部，可以分为2个大的沉降区，其他零星分布.沉降时间主要发生在2016年前，在2014—2015年间年均沉降速率大，最大可达220 mm/yr以上，自2017年后大通矿区沉降变形呈现为明显减弱趋势，2018年后明显沉降区主要集中在矿区东南角.

3.2 2014—2021年多年期时间序列分析

采用SBAS-InSAR方法针对大通矿区获取其2014—2021年的形变情况，SBAS-InSAR年均形变速率结果如图3所示.

由图3可知，大通矿区主要沉降区可分为1#沉降区和2#沉降区，整个矿区面上为零星分布点状变形，高变形区主要分布在矿区东南角，未出现大面积高速沉降区，矿区2014—2021年间年均变形速率约为58 mm/yr.选取矿区6处典型变形点做时间序列分析，其累计变形曲线如图4所示.SBAS1和SBAS2位于矿区东北部的1#沉降区.SBAS1位置累计变形量约340 mm，SBAS2位置累计变形量约320 mm，2条曲线都表明在2014—2017年初变形速率较大，后变形速率逐渐平稳.SBAS3和SBAS4位于矿区东南部的2#沉降区.SBAS3位置累计变形量约270 mm，曲线表明2014—2017年初变形速率较大，后变形速率逐渐平稳，到2020年变形速率又开始加大；SBAS4位置累计变形量已达400 mm，曲线表明该点2014—2017年初变形速率较大，后逐渐平稳，至2021年4月后短暂出现加速变形情况，后逐渐平稳.SBAS5位置累积变形量约160 mm，曲线表明2014年6月至2017年6月间变形速率较大，后逐渐平稳，并在2021年1月至2021年6月间出现加速变形情况.SBAS6位置累积变形量约40 mm，曲线表明该点2014—2019年变形较为平稳，而2020年开始出现加速变形的情况，但形变量值较小.

在矿区内做5条年均沉降速率剖面，剖面线位置如图3所示，剖面结果如图5所示.剖面A-A′表明沿该剖面线方向，地面沉降速率逐渐变大，在剖面线上1 200 m位置处出现转折，变形速率开始减小；剖面B-B′表明沿该剖面方向变形速率较为平稳，局部出现较大沉降变形区；剖面C-C′表明沿该剖面线方向约1 300 m前变形速率平稳，1 300 m后出现较大速率变形区，在1 700m位置处变形速率最大；剖面D-D′表明沿该剖面线方向约1 000 m前变形速率平稳，1 000 m后出现较大速率变形区；剖面E-E′表明沿该剖面线方向整体变形平稳，局部单点存在较大变形，所以剖面呈离散状.

3.3 2020年11月至2021年11月地面變形分析

2014—2021年间时间跨度较长，大通矿区地面变形在长时间跨度下可能会出现阶段性的大幅度变化，所以使用2020年11月至2021年11月间33期影像对工作区做1年期的SBAS-InSAR沉降变形速率观测，其结果如图6所示.

由图6可知，大通矿区2020年11月至2021年11月间最大年均变形速率不超过52 mm/yr，高变形区主要分布在矿区东南角，整个矿区面上为零星分布点状变形，未出现大面积高速沉降区.同样选取矿区6处典型变形点做时间序列分析，其累计变形曲线如图7所示.SBAS1与SBAS2位置累积变形量较小，数据较离散，变形较为平稳，与图4中趋势一致；SBAS3和SBAS4位于矿区东南部，其变形曲线显示该区域于2021年5月后开始加速变形，加速变形于2021年10月结束，后变形较为平稳；SBAS5点位于矿区中部位置，其变形曲线表明该位置于2021年5月后开始加速变形，并于2021年7月后加速变形停止，后变形逐渐平稳；SBAS6点位于矿区西部，其变形曲线表明该点一直处于匀速变形状态.

矿区内沿坡向做5条年均变形速率剖面，剖面线位置如图6所示，剖面线结果如图8所示.与图5年均速率变形分布图一致，沿剖面线方向沉降变形表现为零星点状分布，并未连接成为整体，所以剖面线方向上速率点较为离散.

3.4 对比分析

根据前期勘察设计资料和现场调查，2个集中沉降区主要位于煤矿采空区及影响区内.2012—2015年对矿山进行4期次的地质环境治理示范工程，主要治理措施包括煤矸石弃渣处置、废弃矿井井口封堵、塌陷区土地整治、泵站及田间灌溉配套、沟道防洪、小煤洞渣堆景观区、引水蓄水灌溉系统、排水系统、不稳定斜坡治理、田间道路、林网及矿山公园等工程.经过治理，沉降区范围及沉降量逐渐减小，特别是2017年沉降基本趋于稳定，可见治理工程对控制沉降起到了重要作用，治理效果明显.2020年矿区内部分沉降变形继续增加，但变形速率较之前减小较多，减小约3/4.根据图6和图7变形特征图，结合2021—2022年度综合治理实施方案（见图9）进行对比，可见实施方案中12个不稳定斜坡、3个滑坡和5个沉陷区在InSAR解译图中都有对应的变形解译特征点.

4 结 论

1）利用SBAS-InSAR观测大通矿区2014—2021年间长时间序列变化，表明矿区主要沉降变形区集中在东北部和东南部2个片区，大通矿区2014—2016年间沉降速率最大超过220 mm/yr，自2017年后矿区沉降变形范围和变形量值都呈显著下降趋势且沉降变形逐渐平稳，但是存在局部或某时间段内发生突然变形加速的现象.

2）2014—2021年的5条沉降速率剖面表明主体沉降位于矿区中前部，后部（西部）沉降速率很小；2020年11月至2021年11月间的InSAR结果表明，高变形区主要分布在矿区东南角，整个矿区面上为零星分布点状变形，未出现大面积高速沉降区.

3）通过对比分析，表明前期治理工程对控制矿区沉降变形起到了重要作用，近期变形有所增加，但变形速率减小较大；可见InSAR技术对于隐患点的识别及变形监测具有大范围和高精度的特点，可以为防灾减灾及治理工程评价提供依据.

参考文献：

［1］陈家彪，沈冰，刘飞燕.西南地区矿山环境地质问题研究［J］.矿产综合利用，2007，146（4）：43-46.

［２］刘国林，张连蓬，成枢，等．合成孔径雷达干涉测量与全球定位系统数据融合监测矿区地表沉降的可行性分析［J］.测绘通报，2005，51（11）：13-16．

［3］蒲川豪，许强，赵宽耀，等.利用小基线集InSAR技术的延安新区地面抬升监测与分析［J］.武汉大学学报（信息科学版），2021，46（7）：983-993.

［4］陆会燕，李为乐，许强，等.光学遥感与InSAR结合的金沙江白格滑坡上下游滑坡隐患早期识别［J］.武汉大学学报（信息科学版），2019，44（9）：1342-1354.

［5］李承航，张文春.基于InSAR技术的监测应用及研究进展［J］.北方建筑，2020，5（4）：15-19.

［6］范洪冬，邓喀中，承达瑜.InSAR拓展和融合技术在矿山开采监测中的应用［J］.金属矿山，2008，382（4）：7-10.

［7］张训虎，章磊，郝树宾，等.合成孔径雷达干涉（InSAR）测量技术应用及展望［J］.北京测绘，2014，115（2）：28-31.

［8］高海英，赵争，章彭.时序InSAR的贵州地质灾害监测［J］.测绘科学，2020，45（7）：91-99.

［9］丁刘建，陶秋香，高腾飞，等.SBAS InSAR技术在矿区地面沉降监测中的应用［J］.中国科技论文，2019，14（3）：320-325.

［10］王志勇，张继贤，黄国满.基于InSAR的济宁矿区沉降精细化监测与分析［J］.中国矿业大学学报，2014，43（1）：169-174.

［11］Bamler R，Hartl P.Synthetic aperature radar interferometry［J］.Inverse Probl，1998，14（4）：1-54.

［12］Rosen P A，Hensley S，Joughin I R，et al.Synthetic aperture radar interferometry［J］.P IEEE，2000，88（3）：333-382.

（实习编辑：姚运秀）

Abstract：

InSAR technology is used to monitor the subsidence area of Datong coal mine，and the information of land subsidence deformation and average annual subsidence rate in the mining area from October 2014 to November 2021 are obtained.The results show that the main subsidence deformation areas in the mining area are concentrated in the northeast and southeast，and the maximum subsidence rate is more than 220 mm / yr.After 2017，the settlement deformation range and deformation value of the mining area has shown a significant downward trend，and the settlement deformation has been gradually stable.The preliminary treatment work has played an important role in controlling the settlement and deformation of the mining area.The study can provide an important reference for disaster prevention and reduction and treatment engineering evaluation in alpine mining areas.

Key words：

InSAR;mining area treatment;land subsidence;sedimentation rate