基于DS-InSAR的乌达煤田火区长时序地表形变监测与分析

时间：2024-07-28

李柱，范洪冬，高彦涛，许耀宗

(1.中国矿业大学矿山生态修复教育部工程研究中心，徐州 221116； 2.河南省地质矿产勘查开发局测绘地理信息院，郑州 450006)

0 引言

煤火一般是指地下煤层火灾、近地表煤层火灾、煤矸石火灾和煤堆火灾[1]。大多数煤火是由于不适当的采矿造成的，也有一些是因处理不当的矿山废弃物、闪电和森林火灾引起的[2]。作为一种全球性灾害，煤火对人类健康、煤矿安全和基础设施造成了严重的负面影响。煤火的燃烧不仅消耗了大量宝贵的煤炭资源，同时引发地面沉降、裂缝、裂隙等地质灾害，释放的大量有毒气体(CO，NO2，SO2，H2S)和微量元素(Hg，F，As，Se)还会造成严重的环境问题[3]。因此，全面、详细以及准确地监测煤火，对煤火的演变趋势、自燃规律和燃烧状态进行研究以及开展煤火治理活动具有重要的意义。

遥感方法具有获取信息速度快、监测区域面积大、周期短、全天时、全天候和经济效益高的特点，这些独特的优势可以弥补传统煤火探测与监测技术的不足[4]。利用遥感方法进行地面沉降监测和温度反演，已经成为煤火探测与监测领域的一个重要方向。国内外学者利用热红外遥感技术对乌达煤田煤火区进行了大量研究[5-8]，可以较高精度地探测出煤火区及其演变趋势，但热红外遥感技术很难实现对深部煤火、煤火燃烧阶段以及煤火引起的地表形变的监测，这使得热红外遥感技术监测煤火时存在一定的不足。

差分合成孔径雷达干涉测量(differential interferometric synthetic aperture Radar，D-InSAR)技术作为一种先进的对地观测手段已在各种变形监测领域得到了广泛的应用，但是该技术易受时空失相关及大气延迟等因素影响，难以在长时序上得到高精度的地表形变结果[9]，而近年来发展的分布式目标合成孔径雷达干涉测量(distributed scatterer interfero-metric synthetic aperture Radar，DS-InSAR)方法则有效地解决了长时序上低相干地区监测点位密度不足的问题[10]。目前，国内外一些学者利用合成孔径雷达干涉测量(interferometric synthetic aperture Radar，InSAR)技术对煤火区的探测与监测进行了相关研究，大多数研究侧重于将InSAR监测的地表形变结果作为辅助信息与地表热异常信息相结合进行煤火区域的探测[11-14]，而对煤火引起的地表形变规律研究较少[15-16]；另外，对于乌达煤田煤火的研究多数是利用热红外遥感技术进行探测与监测，只有少数研究采用InSAR技术[17-19]，并且缺乏近些年的地表形变监测结果。

针对上述问题以及煤田火区的地表形变复杂、失相干严重等现象，本文研究了一种基于DS-InSAR技术的煤田火区监测方法，利用该方法对乌达煤田火区2017年3月—2019年4月的63景Sentinel-1A影像进行了处理，通过监测结果对乌达煤田及其煤火区的地表形变进行了分析，验证了该方法的适用性。

1 研究区域概况及数据源

1.1 研究区域概况

乌达煤田位于中国内蒙古自治区西南部的乌海市乌达区，如图1所示。该区东邻黄河与鄂尔多斯黄土高原，南部是乌达工业园区，西邻贺兰山脉北部边缘，北邻乌兰布和沙漠。煤田呈南北向斜状，面积约为55 km2，包括五虎山、黄白茨和苏海图3个煤矿，海拔从1 100～1 300 m不等，平均高程为1 200 m。气候为典型的温带大陆性气候，年降水量为168 mm，蒸发量为3 500 mm左右，常年干旱少雨。

图1 研究区域Fig.1 The study area

1.2 数据源

Sentinel-1A卫星由欧州航天局于2014年4月3日发射，采用近极地太阳同步轨道，轨道高度为693 km，重访周期为12 d，搭载了一个C波段的合成孔径雷达(synthetic aperture Radar，SAR)，能提供超宽幅、干涉宽幅、条带以及波浪4种成像模式。本文选取2017年3月—2019年4月期间的63景干涉宽幅模式成像下的SAR影像为数据源，其距离向和方位向的像元尺寸分别为2.3 m和14.0 m。为去除干涉图中的地形相位，借助空间分辨率为90 m的SRTM DEM数据作为外部数字高程模型(digital elevation model，DEM)进行差分干涉处理。

2 DS-InSAR方法原理

Ferretti等[20]于2011年提出了第二代永久散射体技术SqueeSAR。拉开了DS-InSAR方法研究的序幕。与传统的时序InSAR方法不同的是，DS-InSAR方法在其基础上增加了同质点选取和相位优化2个步骤，有效地解决了传统时序方法相干点选取数目少和空间分布不均匀的问题。

2.1 同质点选取

(1)

在匀质区域的单视复数SAR影像服从瑞利分布，其变异系数CV为常数，即

(2)

假设像元P后向散射性稳定，则式(1)可以转换为：

(3)

2.2 特征值分解相位优化

(4)

式中: 向量y=[y1y2,…,yN]为分布式目标的同质点在N景SAR影像上的复数观测向量x=[x1,x2,…,xN]经过振幅归一化处理得到；NSHPs为分布式目标同质点的数目； H为Hermitian共轭转置。

(5)

DS-InSAR数据处理具体流程如下： ①根据设置的时空基线阈值选择干涉对，并生成干涉图； ②利用FaSHPS方法对强度影像进行同质点识别，并通过特征值分解法对干涉图进行相位优化； ③利用外部DEM去除干涉图中的地形相位生成差分干涉图，并利用优化后的相位计算同质点的时空相干性，设定相干性阈值选取分布式目标(distributed scatterer，DS)点； ④利用差分干涉图对DS点上的相位进行解缠，并建立解缠相位与地表形变速率、DEM误差和残余相位之间的观测方程； ⑤利用奇异值分解法估算地表形变相位和DEM残余相位，并对残余相位进行时空滤波分离出非线性形变相位、大气相位和噪声相位； ⑥将线性形变相位与非线性形变相位相加解算研究区域视线向(line of sight，LOS)的地表时序形变。数据处理流程如图2所示。

图2 数据处理流程Fig.2 Data processing flow chart

3 监测结果与分析

3.1 监测结果对比

分别采用TCP-InSAR与DS-InSAR方法对覆盖乌达煤田的63景Sentinel-1A数据进行处理，获得了2017年3月—2019年4月时间段内的地表时序形变信息，年平均形变速率如图3所示。从图3可以看出，二者虽在空间上形变趋势具有较好的一致性，但DS-InSAR监测结果能够更好地反映出煤田的地表形变情况，并对研究区域形变明显的8个区域进行A—H编号。

(a) TCP-InSAR (b) DS-InSAR

TCP-InSAR方法共选出了156 853个监测点； DS-InSAR方法共选出了194 585个DS点。从选取的监测点数量而言，DS-InSAR能够选出更多的DS点，监测点位密度比TCP-InSAR提高1.24倍，并且在低相干区域也可以获得足够数量的监测点；从监测点位的空间分布来看，DS-InSAR方法选出的DS点在保证点密度的同时分布更均匀，在一定程度上提高了形变监测结果的解算精度，为研究区提供了更为详细准确的形变信息，有效地弥补了在低相干区域TCP-InSAR等传统时序形变监测方法中存在的缺陷。

为验证DS-InSAR监测结果的可靠性，通过对比分析TCP-InSAR和DS-InSAR获得的形变速率值来进行交叉验证。本文以TCP点为基准，选取最邻近的DS点作为同名点，共筛选出95 760对同名点，通过同名点对的形变速率值及其差异的统计分析，绘制了形变速率相关性和差异分布直方图，如图4所示。通过形变速率相关性(图4(a))可知，DS-InSAR与TCP-InSAR方法获得的地表形变速率之间具有较强的相关性，利用Pearson计算的相关系数R为0.84。由形变速率差异分布直方图(图4(b))可知，2种方法获得的形变速率差异较小，两者之间的标准差为5.83 mm/a。因此，在缺少实测水准数据的情况下，与TCP-InSAR方法监测结果的交叉验证，表明了DS-InSAR地表形变监测结果是较为可靠的。

3.2 乌达煤田地表形变分析

为了能够全面分析乌达煤田地表形变特点，本文利用DS-InSAR对乌达煤田进行了长时序的地表形变监测，时序形变累积如图5所示。图5中的14幅图是每隔一段时间相对于起始时间2017年3月17日的累积形变量。从图5中可以看出，在监测时间段内乌达煤田存在严重的地表形变现象。在空间上，沉降区域分布不均匀。在时间序列上，随着时间的推移研究区域内逐渐出现了形变量级与沉降范围大小不同的多个沉降区域。其中，A，B，C这3个区域的形变现象最为严重，沉降范围最广，这3个区域的地表形变在整个监测时间段内持续存在，地表形变量与沉降范围越来越大，导致B和C区域从2018年8月15日之后逐渐连接形成一个大的沉降区域； G区域的地表形变是从2017年5月16日开始产生，前期形变变化较小，到后期形变变化加快； D和F区域从2017年8月20日之后地表开始产生形变，并逐渐形成多个沉降中心； H区域的地表形变从2017年12月18日开始发生地表形变； E区域产生地表形变的时间最晚，但地表形变变化较快，从2018年10月14日—2019年4月12日半年的时间里产生了与D区域相同量级的形变。为更详细具体地分析地表形变情况，利用图3(b)地表形变速率图对其进行定量分析。从图5中可知，地表形变最为严重的A，B，C这3个区域的最大形变速率分别为-178.4 mm/a，-166.3 mm/a和-215 mm/a，在这3个沉降区域的沉降中心，由于形变梯度大，地表形变量在监测时间段内超过时序InSAR的监测能力，失相干现象严重，无法得到足够密度的DS点；另外，F和G这2个区域同样存在较为严重的形变，最大形变速率分别为-96.5 mm/a和-100.9 mm/a；相比之下，D，E和H区域的形变量相对较小，最大形变速率分别为-46.5 mm/a，-53.9 mm/a和-48.5 mm/a。

图5 乌达煤田地表形变时序累积图

3.3 煤田火区地表形变分析

乌达煤田严重的地表形变现象主要与煤火燃烧及人工采挖有关[17,19]，根据文献[7]中2018年1月底实地调查获取的部分区域的明显烟点和煤火点，这些实测的煤火点附近地区伴随着地表形变现象。为了研究煤火燃烧引起的地表形变，对实测煤火区的地表形变进行重点分析。

3.3.1 煤田火区地表形变时序分析

为分析煤田火区的时序地表形变规律，将实测煤火点编号为H1—H6(图3(b))，对应D，A和G这3个沉降区域。选取3个区域中实测煤火点附近且累积形变量最大的5个DS点作为特征点分析地表形变，编号记为P1—P5(图3(b))。为了验证特征点时序形变的可靠性，从TCP-InSAR监测结果中选取距离特征点最近的TCP点作为验证点对其进行验证，由于TCP-InSAR选点的原因，只在P3和P5点旁选了P3′和P5′点。图6(a)—(c)分别为D，A和G这3个沉降区域的特征点及验证点时序形变图，特征点与验证点之间的时序形变趋势基本吻合，具有很好的一致性，2点间的时序形变值整体存在偏大或偏小的问题，但考虑到2点之间具有一定的距离，属于合理范围，所以可以利用选取的特征点进行地表时序形变分析。

由图6可知，3个沉降区域的形变量级与时序形变趋势各不相同，不同时间段地表形变变化快慢也不相同。P1从2018年5月底才开始发生形变并持续到监测时间结束，这段时间内地表形变变化快，煤火燃烧剧烈。处于同一沉降区域的P2，P3和P4开始发生形变的时间各不相同，P2点最晚，从2017年9月开始，并且中间一段时间因煤火处于休眠阶段未发生明显的地表形变现象； P3点地表形变变化缓慢，但在整个监测时间段内形变持续发生；相比之下，P4点从2017年6月开始出现形变现象，且在之后时间段内出现了3个未发生形变的时期。P5点的时序形变趋势更加复杂，多个时间段内出现形变停滞现象，但在发生形变的时期内形变变化较快。煤火区的这些复杂地表时序形变现象归因于煤火燃烧阶段的复杂性、燃烧强度的非均质性，并且可以从这些特征点时序上的地表形变变化快慢推断煤火的燃烧强度及燃烧阶段。从整体上分析这5个特征点时序地表形变规律，可以发现秋冬季节的地表形变变化相对较快，这是由于乌达煤田处于干旱半干旱地区，秋冬季节干旱少雨，煤层易发生火灾[23]，导致地表形变变化较快。

3.3.2 典型煤火区地表形变特征分析

为分析煤田火区的地表形变特征，选取3个沉降区域中形变最严重的A区作为典型煤火区。为清晰直观地表征区域A地表形变的空间形态，将地表累积形变图(图7(a))进行克里金插值，并将插值之后的图转换成三维形变图(图7(b))。从图7(b)可知，区域A的地表形变在空间上变化较快，整体形变规律复杂，边缘形状不规则，且具有多个发育程度不同的沉降中心。从三维形变图的形变等值线中可知，这些沉降中心附近的等值线密集，向四周逐渐变得稀疏，对应的沉降中心凹陷较深，向四周逐渐平缓，边缘近似长椭圆形，并且这些沉降中心的形变延伸方向随着煤火蔓延和燃烧强度而各不相同。

(a) 剖面线示意图 (b) 三维形变图

为定量分析典型煤火区的地表形变特征，在图7(a)区域A的沉降中心附近给出了2条剖面线a1和a2，并沿剖面线提取地表累积形变量，结果如图8所示。由图8可知，剖面线a1穿过了3个形变量级与沉降范围大小不同的沉降中心，最严重的沉降中心累积形变量达到-311.6 mm；剖面线a2经过2个发育程度相似的沉降中心，最大累积形变量从左到右分别为-369.5 mm和-355.2 mm；这些沉降中心到边缘的形变量变化快，形变梯度较大。

图8 典型煤火区剖面累积形变图