基于时序InSAR技术的新井煤矿地表沉降监测.pdf

第39卷第4期2023年8月山西大同大学学报（自然科学版）Journal of Shanxi Datong University(Natural Science Edition)Vol.39 No.4Aug.2023基于时序InSAR技术的新井煤矿地表沉降监测马飞1，李日印1，刘小鹏2，侯景怡1，崔喆森1，张垚杰1（1.长治学院 计算机系，山西 长治 046011；2.自然资源部 第二地形测量队，陕西 西安 710054）摘要：为探测内蒙古新井煤业有限公司露天矿区边坡失稳前地表形变情况，采用SBAS-InSAR技术对该矿区2022年1月至2022年11月的15景Sentinel-1A数据进行干涉处理，解算出地表2 475个形变点，最大形变量超过160 mm，最大形变速率达209 mm/a。证明SBAS-InSAR技术可以高效地识别研究区变形位置及其形态，为矿区地质灾害预警提供数据支撑。关键词：时序InSAR；新井煤矿；地表形变中图分类号：TP79文献标识码：Adoi：10.3969/j.issn.1674-0874.2023.04.024大规模开发利用矿产资源会直接影响矿区安全1。例如：2023年2月22日，内蒙古新井煤业有限公司露天煤矿发生大面积边坡失稳事故，造成人员财产损失。因此，对矿区地表变形进行监测具有重要的现实意义2-4。差分干涉合成孔径雷达(DInSAR)、小基线集In-SAR(SBAS-InSAR)和永久散射体 InSAR(PS-InSAR)技术都是监测采矿沉陷、冻土沉陷、建筑物沉陷和滑坡的常用方法5-7，具有空间覆盖面大、成本低、精度高、全天候观测能力等优点，而常规测量技术（如水准测量8、三维激光扫描9和全球导航卫星系统10等）无法做到这一点。为了满足开采沉陷长期监测的需要，削弱大气相位的影响，时序InSAR提出了一种简单有效的方法，即通过对一系列展开的干涉图按时间跨度加权和平均来估计平均视线开采沉陷速度11-12。小基线集 InSAR(SBAS-InSAR)可以很好地估计变形结果13-15，在许多开采沉陷区得到了应用。现对内蒙古新井煤矿应用SBAS-InSAR技术，识别该矿区地表沉陷范围，对沉陷区下沉数值进行解算，分析灾前微小形变规律，以预防同类事故的再次发生。1 SBAS-InSAR技术原理SBAS-InSAR通过设置特定的时间基线阈值来合并小的基线干扰对，根据最小参数准则执行最优的形变相位序列，以限制时间和空间消相干。假设已经采集了在时间序列(t0，t1，tn)相同区域中获得的（n+1）个SAR图像。设置特定的时间基线阈值后，可以形成m个干涉对，且m满足：n2 m n（1）SBAS-InSAR在进行参数反演之前，需要将所有图像与超级主影像配准，并通过多视处理、地形差分和相位解缠来生成m幅差分干涉图图像。考虑在时刻ta和tb（假设ta tb）采集的两幅差分干涉图的SAR图像，在不考虑大气延迟等误差的情况下，方位向和距离方向上具有坐标(x，r)的像素的干涉相位组成为：(x)(tb,x)-(ta,x)+n=4 d(tb,x)-d(ta,x)+n=4v(x)Ti+n(2)式中：()tb,x 和()ta,x分别为相位形变值；为雷达波长；n为随机噪声；d()tb,x 和d()ta,x分别为卫星视线方向的形变量；v(x)为像元的形变速率；Ti为干涉对的时间基线长度。这样可以组成m个方程组：A=（3）式中：A为m n的系数矩阵，每行对应一幅干涉图，收稿日期：2023-03-05基金项目：山西省基础研究计划项目202103021223381作者简介：马飞（1988-），回族，山西长治人，博士，讲师，研究方向：遥感技术应用。E-mail：文章编号：1674-0874（2023）04-0123-05山西大同大学学报(自然科学版)2023年每列对应一个时间上的SAR图像，主影像所在列+1，辅图像所在列-1，其余列为0。如果m n，且A的秩是n，则最小二乘解得的估计值：=(ATA)-1AT（4）由式（4）可知，SBAS技术的关键在于求时间序列干涉图间形变速率的最小范数最小二乘解，可以利用时序干涉图中解算出的形变相位值，作为时间域的约束条件解，简化形变速率解算的复杂度，达到反演形变速率的目的。假设V为形变速率，P为形变参数，以下约束条件被满足：V=QP（5）将式（5）代入BV=得：BQP=（6）在干涉图中任意像元（x,r）的解算模型转换为：(ti)=-v(ti-t0)+12-a(ti-t0)2+16-a(ti-t0)3（7）式中：LOS方向形变的平均加速度变化率a、相位平均速率v、平均加速度a 用下式即可解得：PT=v a a T=1t1-t02()t1-t0361t2+t1-2t02()t2-t03-()t1-t036()t2-t11tN+tN-1-2t02()tN-t03-()tN-1-t036()tN-tN-1(8)2 研究区概况及实验结果2.1 研究区概况图1为内蒙古新井煤业露天矿区地理位置，该矿区位于贺兰山南段，矿区面积 1.344 8 km2，地势北高南低，海拔高程1 159.03 1 490.70 m。矿区为低山丘陵区，属高原侵蚀性丘陵地貌，地形切割一般，基岩裸露，植被稀疏，为荒漠地区。气候属大陆性干旱荒漠气候，日温差大，夏季炎热，冬季严寒，气候干燥，雨量稀少，年最高气温 38，最低-28，平均 8 9，年降水量为146 198 mm，年蒸发量最高2 000 mm以上。地表植被属于荒漠化草原向草原化荒漠过渡地带。地表多沙质化、砾石化和有龟裂结皮。图1 新井煤业露天煤矿地理位置通过分析新井煤业露天矿区地表植被和地形地貌特征，该矿区地表属于荒漠地区，植被稀少，适用C波段数据，但地面建筑物、工矿构筑物、大型岩石等永久散射体物较少，适合采用SBAS-InSAR对矿区地表沉降情况进行监测。2.2 SAR数据选取及预处理表1为覆盖该矿区的Sentinel-1A数据列表，包括卫星拍摄时间、影像模式、极化方式和轨道方向等信息。根据欧空局(ESA)网站查询哨兵-1A数据情况发现，该卫星2022年能够覆盖研究区的影像只有15景，观测时间是2022年1月到2022年11月。根据干涉处理实验需要，下载了该区域干涉宽幅模式数据，单幅影像幅宽为250 km，距离向分辨率为4 m，方位向分辨率为13 m，多视比例设置为3 1。同时下载到所有影像的精密轨道数据，用于干涉数据处理过程中的图像配准、去平地相位等流程，提高数据处理精度。另外，为提高数据处理效率，对原始SAR影像进行了裁切。表1 覆盖矿区的Sentinel-1A数据列表ID123456789101112131415时间202201072022011920220131202202122022022420220308202203202022040120220413202205072022051920220612202207302022081120221115模式IWIWIWIWIWIWIWIWIWIWIWIWIWIWIW极化方式VVVVVVVVVVVVVVVVVVVVVVVVVVVVVV轨道方向升轨升轨升轨升轨升轨升轨升轨升轨升轨升轨升轨升轨升轨升轨升轨1242023年(a)20220401-20220107(b)20220401-20220119(c)20220401-20220131(d)20220401-20220212(e)20220401-20220224(f)20220401-20220308(g)20220401-20220320(h)20220401-20220413(i)20220401-20220519(j)20220401-20220612-图3 以20220401为主影像生成的的干涉图2.3 时序InSAR 技术数据处理流程2.3.1短基线集干涉图组合图2为干涉对时空基线图。将选取的15景影像按照如下策略进行自由组合，形成干涉对：时间基线T 100d，空间基线B 200m，以 20220401 为超级主影像。最终得到43个干涉对。图2 SAR干涉对时空基线图2.3.2干涉图生成图3是在新井煤业露天矿区生成的以20220401为主影像的差分干涉图，该图像为SAR斜距坐标系统。图中可以看到，矿区中心位置形变条纹非常清晰，随着时间的推移，矿区中心位置形变量在不断变化；时间基线越长，干涉条纹周期变化越多；在图(i)、(j)等干涉对的中心位置出现干涉条纹叠加和失相干。2.3.3 SBAS时间序列形变和形变速率结果图 4 为新井矿业露天矿区形变速率图，红色代表形变方向与卫星视线方向一致，可解算为下沉值，蓝色代表形变方向与卫星视线方向相反，可解算为抬升值。图中可见，该矿区大部分区域相对稳定，只有在此次滑坡事故中的核心区发生剧烈形变，整体形变速率非常大，沉降中心位置最大形变速率达到-160 -209 mm/a，说明此处滑坡 2022 年开始已经出现不稳定状态。另外，在该区域的东南侧也出现一个小范围的形变区，该区域的形变速率在-60-110 mm/a，需要引起相关部门的高度关注，避免发生新的地质灾害。马飞等：基于时序InSAR技术的新井煤矿地表沉降监测125山西大同大学学报(自然科学版)2023年图4 灾变核心区形变速率图为进一步分析新井煤业露天矿区核心区灾前形变情况，将该形变区域数值提取出来并进行处理。该区域共解算得到2 475个形变点，图4显示该区域自边坡的坡顶向坡底发生形变量由大到小，2022年1月至11月坡顶累计形变量超过170 mm，坡底累计形变量超过150 mm。图 5为边坡在 2022年的形变速率剖面图，横坐标为距离边坡坡顶GCP1号点的距离，纵坐标为形变速率值。为进一步分析该矿区边坡失稳前的形变规律，在边坡选取了四个点命名为GCP1 4进行分析，分别分析其剖面线和形变时序图。图中可知，边坡顶部下沉速率远远高于坡底位置，自坡顶向坡底形变速率逐渐下降，形变速率从-180 -140 mm/a。由于该区域边坡角度较大，坡顶区松散层不断下沉，经过一段时间的累积，需要采取任何进行边坡加固。图5 形变速率剖面图图6为新井煤业露天矿区灾变核心区内边坡形变点的时序下沉变化情况。由图可知，在边坡顶部的GCP1号点累积形变量最大，在底部的GCP4号点累积形变量稍小，这与图5中的剖面图一致；该边坡从2022年2月开始出现下沉，2022年3月至6月累积形变趋势开始增大，这可能与该地区前期转暖冻融解冻有关，也可能与矿区露天开采活动有关。图6 边坡灾变点形变时间序列图3 结论(1)利用SBAS-InSAR方法可以有效识别出露天采矿区地表变形，这是一种快速、高效、低廉的技术，可以为矿区地表边坡失稳、地面塌陷、地面沉降、地裂缝、滑坡、崩塌等地质灾害提供预警。(2)内蒙古新井煤业露天采矿区核心灾变区在2022 年 1-11 月持续发生下沉形变。分析 SBAS-In-SAR方法获取的2 475个形变点发现该区域最大形变量超过 160 mm，沉降速率超过 200 mm/a；在核心灾变区东南侧亦有一个缓慢形变区，累积形变量稍小，但需要有关部门关注，避免发生新的地质灾害。(3)内蒙古新井煤业露天采矿区失稳边坡在灾前一年间持续发生下沉，且边坡坡顶形变速率高于坡底速率，坡顶形变速率超过180 mm/a，顶部持续微小的下沉形变是导致边坡应力失衡的重要原因。(4)由于缺少地表监测点，故 SBAS-InSAR 的形变监测值无法得到验证，但该技术获得的时间序列累积变形能更直观地反映沉陷盆地的时间序列变化趋势，在探测变形区域范围和形态方面具有优势，有助于定性确定形变类型。致谢：本实验的 Sentinel-1A 数据和 SRTM 外部DEM数据均来自欧空局网站，在此表示感谢。1262023年参考文献1 黄乐亭，王金庄.地表动态沉陷变形的 3个阶段与变形速度的研究 J.煤炭学报，2006，31（4）:420-424.2 郝延锦，吴立新，戴华阳.用弹性板理论建立地表沉陷预计模型 J.岩石力学与工程学报，2006（S1）:2958-2962.3 赵晋帅.GPS RTK技术在滑坡动态实时变形监测中的应用 J.山西大同大学学报（自然科学版），2018，34（6）:72-74.4 袁媛.基于 WEB-GIS的大同市矿山地质环境动态遥感监测系统 J.山西大同大学学报（自然科学版），2020，36（6）:91-94.5 BERARDINO P，FORNARO G，LANARI R，et al.A new algorithm for surface deformation monitoring based on small base-line differential SAR interferograms J.IEEE Transactions on geoscience and remote sensing，2002，40（11）:2375-2383.6 DENG Z，KE Y，GONG H，et al.Land subsidence prediction in Beijing based on PS-InSAR technique and improved Grey-Markov model J.GIScience&remote sensing，2017，54（6）:797-818.7 LI Z W，YANG Z F，ZHU J J，et al.Retrieving three-dimensional displacement fields of mining areas from a single InSAR pairJ.Journal of Geodesy，2015，89（1）:17-32.8 HAO Y，WU L，DAI H.Establishing a ground settlement prediction model with elastic slab theory J.Yanshilixue YuGongcheng Xuebao/Chinese Journal of Rock Mechanics and Engineering，2006，25:2958-2962.9 YINGHUA W，CHANG X.Using genetic artificial neural network to model dam monitoring data C/2010 Second Internation-al Conference on Computer Modeling and Simulation.IEEE，2010，2:3-7.10 张小红，胡家欢，任晓东.PPP/PPP-RTK新进展与北斗/GNSS PPP定位性能比较 J.测绘学报，2020，49（9）:1084-1100.11 YIN H，ZHU J，LI Z，et al.Ground subsidence monitoring in mining area using DInSAR SBAS algorithm J.Acta Geodeticaet Cartographica Sinica，2011，40（1）:52-58.12 WASOWSKI J，BOVENGA F.Investigating landslides and unstable slopes with satellite multi temporal interferometry:Cur-rent issues and future perspectives J.Engineering Geology，2014，174:103-138.13 LI Y，YANG K，ZHANG J，et al.Research on time series InSAR monitoring method for multiple types of surface deformationin mining area J.Natural Hazards，2022，114（3）:2479-2508.14 NG A H M，GE L，YAN Y，et al.Mapping accumulated mine subsidence using small stack of SAR differential interferogramsin the Southern coalfield of New South Wales，Australia J.Engineering Geology，2010，115（1-2）:1-15.15 VADIVEL S K P，KIM D，KIM Y C.Time-series InSAR analysis and post-processing using ISCE-StaMPS package for mea-suring bridge displacements J.Korean Journal of Remote Sensing，2020，36（4）:527-534.Monitoring of Surface Subsidence in Xinjing Coal Mine Based on Time Series InSAR TechnologyMA Fei1,LI Ri-yin1,LIU Xiao-peng2,HOU Jing-yi1,CUI Zhe-sen1,ZHANG Yao-jie1(1.Department of Computer,Changzhi University,Changzhi Shanxi,046011;2.The Second Topographic Surveying Brigaade of SBSM,Xian Shaanxi,710054)Abstract:In order to detect surface deformation before slope instability in open-pit mining area of Inner Mongolia XinjingCoal Industry Co.,LTD.,15 Sentinel-1A data from January 2022 to November 2022 were interfered by SBAS-InSAR technology,and 2 475 surface deformation points were calculated,with the maximum shape variable exceeding 160 mm.The maximum defor-mation rate is 209 mm/a.It is proved that SBAS-InSAR technology can efficiently identify deformation location and shape in thestudy area,and provide data support for geological disaster warning in mining area.Key words:time series InSAR;Xinjing coal mine;surface deformation责任编辑 王东马飞等：基于时序InSAR技术的新井煤矿地表沉降监测127