遥感方法在宁武煤田煤火识别中的应用

高玉荣 ，隋 刚 ，张新军 ，孔嘉嫄 ，张和生

（太原理工大学 矿业工程学院, 山西 太原 030024）

0 引 言

山西省煤炭资源丰富,2022 年全省规模以上原煤产量为130 714.6 万t，同比增长8.7%。大规模开采带来严重的生态问题，其中煤的自燃危害尤其严重。煤火燃烧不仅造成巨大的能源浪费和经济损失[1-2]，还会产生大量有害气体，危害人民生命健康，同时也会引起地表裂缝和塌陷，对当地的基础设施产生影响。因此，准确识别煤田火区范围对于煤火监测与治理具有重要意义。

近年来，大量学者对新疆、内蒙古、宁夏等地较大范围煤田火区进行了探测与研究。自20 世纪60年代以来，国内外已经发展了多种煤田火区探测技术，主要包括磁探法、物探法、化探法和遥感法等[3]。其中利用煤层燃烧的物理、化学特性进行探测的方法由于技术限制、成本较高，均不适用于大面积火区[4]。而遥感技术具有探测范围大、获取周期短、时效性强、经济效益高等优势，已成为煤田火区识别领域的重要发展方向[5]。

利用煤火燃烧产生热量，传导至地表造成温度异常这一特性，学者们通过地表温度反演提取热异常信息对煤火进行识别和监测。地表温度反演算法总体分为单通道法、分裂窗法、日夜法、TES、分裂窗和温度比辐射率分离结合法[6]。邱程锦等[7]通过大气校正的反演方法对Landsat TM/ETM 进行处理，并提取温度异常区域来判定煤火范围；李峰等[8]利用TES 算法反演内蒙古乌达矿区4 个时期的地表温度，并采用自适应梯度阈值法提取对应火区范围来监测并评估乌达煤火的治理效果。但深部煤火产生的热量可能传导不到地面，易造成煤火区漏判，另外砂岩吸热、城市热岛效应等因素也会导致与煤层燃烧无关的高温异常区[9]，从而造成误判，所以仅利用热红外遥感提取热异常区域的方法识别煤田火区范围存在一定缺陷。针对煤火燃烧会造成地下空洞，引起地表塌陷变形等特性，少部分学者利用In-SAR 技术进行煤火识别和监测。如JIANG L 等[10]通过PS-InSAR、Stacking 和D-InSAR 三种合成孔径雷达方法对乌达煤田进行了监测，证实了利用In-SAR 手段进行地表形变分析，可以有效判断地下煤火的燃烧情况；RIYAS 等[11]为了表征印度Jharia 煤矿火灾的时空动态，利用新小基线子集(N-SBAS)技术计算了该煤田2017-2020 年地表变形时间序列，但煤矿开采[12]、地质灾害等也会引起地表形变，从而对煤火识别产生影响。

综上所述，现有方法虽然在煤火识别方面取得了较多研究成果，但识别特性单一，易受地物吸热、煤矿开采等因素的影响。因此采用融合热异常信息和地表形变信息来识别煤田火区可以有效克服单一方法的不足，提高煤火识别准确性。笔者利用多源遥感手段，通过融合卫星热红外技术与雷达技术，首次对山西省宁武煤火进行识别分析研究。

1 研究区域概况与数据源

1.1 研究区域概况

研究区位于宁武县东寨镇境内，地理坐标为东经112°03′57″～112°12′39″，北纬38°46′41″～38°53′47″。区内地势高峻，山岭纵横，海拔在2 000 m左右，地形总体西高东低，最高处位于西部，海拔2 020 m，最低处位于东部深沟底，海拔1 670 m，最大高差为350 m。属温带大陆性气候，寒冷多大风，昼夜温差较大，年平均气温7 ℃。研究区域内的煤矿主要分布在中部，呈南北走向，包含寺耳沟煤矿，小西沟煤矿，车道沟煤矿，三马营煤矿等多个煤矿区。由于煤矿开采引发了地面塌陷、地裂缝、最终形成火风压，造成了多个煤火自燃区域[13]，加之小煤窑不规范化的开采方式，产生的煤火问题未能得到及时解决。

1.2 数据源

研究采用2019-08-14 在宁武过境的Terra 卫星上的ASTER(Advanced Spaceborne Thermal Emission and Reflection Radiometer)传感器获取的夜间遥感数据。由于晚上太阳辐射消失，所获取的热红外遥感影像能更准确地提取热异常信息[14]。卫星过境时间为14:34:12(UTC)，云覆盖量为2%，数据产品级别为 L1T，轨道号为 222/211，时相和云覆盖量满足试验要求。所用的ASTER 数据参数见表1。

表1 热红外数据参数Table 1 Thermal infrared data parameters

采用的雷达影像数据为Sentinel-1A 数据，用于煤田火区的地表变形监测。试验选用工作模式为IW、极化方式为VV 的SLC 数据。选取2018 年6月至2020 年1 月的26 景影像(SAR 数据成像时间见表2)用于差分干涉处理。外部参考DEM 选用空间分辨率为12.5 m 的ALOS PALSAR 数据。

表2 SAR 影像成像时间统计Table 2 Statistics of SAR imaging time

2 研究方法

以山西宁武部分煤田为研究区域，ASTER 夜间热红外影像、Sentinel-1A 影像为试验数据，采用多源遥感融合的方法进行煤火区域识别。首先采用温度比辐射率分离算法(ASTER-TES)进行地表温度反演，通过一定的阈值提取研究区的热异常范围，结合SBAS-InSAR 技术提取地表持续变形信息，将热异常信息与地表形变信息空间叠加分析后得到研究区内煤火的疑似分布区域，再根据部分实测煤火范围进行验证分析。研究流程如图1 所示。

图1 研究流程Fig.1 Research process

2.1 温度反演

温度反演中2 个至关重要的因素是温度和地表比辐射率。传统的温度反演算法常常假设比辐射率已知，来求解地表温度，而TES 算法则根据一定的先验知识作为约束条件，同时求解比辐射率和温度。由于ASTER TES 算法是针对ASTER 数据反演的官方算法[15]，充分利用了ASTER 数据的5 个热红外波段，且吸收了发射率归一化NEM、光谱比值RATIO、和最大最小发射率差值MMD 3 个模块的优点，精度较高，故采用此算法进行地表温度反演，其包括3 个部分:

1)NEM 模块：初步估算目标表面温度并且从辐射亮度观测中减去反射的大气辐射。

式中：Ri为第i(i=10，11，···，14)波段的地表辐射亮度；Lgrd为包含大气下行辐射的地表辐射亮度；εmax为初始最大发射率，取0.96；Latm↓为大气下行辐射亮度；λi为波段i(i=10，11，···，14)的波长；TNEM为NEM模块的输出温度；c1为第一辐射常数，c2为第二辐射常数，其中h为普朗克常数，6.626 176×10-34J·s；c为真空光速，2.997 924 58×108m/s；k为玻尔兹曼常数，1.380 6×10-23J/k。

2) RATIO 模块：利用NEM 模块估算的发射率计算相对发射率值。

式中：βi为第i(i=10，11，···，14)波段的相对比辐射率。

3) MMD 模块：进一步估算发射率和温度。

其中：εmin为发射率最小值；max 为发射率εi最大值(εmax)所对应的波段。当MMD＜0.03 时，灰体的精度很低，不再使用MMD 方法，此时直接将εmin设为0.983。数值根据水体和浓密植被的性质确定[16]。

对比发现，热红外数据与雷达数据在地理位置上存在偏差，但雷达数据难以选取控制点，故以与雷达数据空间位置较为匹配的Landsat 影像为基准影像，利用Arcgis 对ASTER 数据进行地理配准。随后采用ENVI 软件对地理配准后的ASTER 数据进行辐射定标，大气校正处理，并基于IDL8.7 平台实现该算法。重复式(1)—式(9)，直到迭代计算的相邻2 次的温差＜0.3 K 或迭代次数＞12 次为止，得到相应的地表温度值，通过一定的阈值提取温度异常区域。经统计，温度反演的结果符合正态分布，将(μ+2σ)作为高温阈值进行地表热异常的提取，其中，μ为温度反演结果中统计值的数学期望，σ为统计值的标准差[17]。

2.2 地表沉降

短基线集时序分析技术(SBAS-InSAR)是一种基于多幅SAR 影像的时间序列方法，克服了传统D-InSAR 中存在的时间、空间失相关问题，通过时空基线较短的干涉对提取地表形变信息[18]。

对于在tA，tB(A，B=0，1，···，N，A≠B)时刻生成的第K幅干涉图，其任意像元的干涉相位值为[19-20]:

式中：Φ[tA，x，y]，Φ[tB，x，y]分别为tA，tB时刻相对于初始时刻t0的形变相位；x，y分别为方位向与距离向坐标;λ为雷达波长；d(tA，x，y)和d(tB，x，y)为相对于初始时刻t0的视线方向的形变量。

由于研究区域内植被较茂密，为避免完全空间失相关，将空间基线阈值和时间基线阈值分别设置为45 m 和365 d，对输入的SAR 数据进行干涉像对的配对，并对部分像对进行3D 解缠；采用Delaunay MCF 的方法进行解缠，Goldstein 进行滤波；选择合适的GCP，估算和去除残余的恒定相位和经解缠后还存在的相位坡道;经过两次反演得到形变速率结果，并将形变速率结果从斜距投影转换为地理投影。由于研究区域的煤火初始燃烧时间未知，缺乏地表监测资料，根据相关文献，将形变阈值确定为5 mm/a，沉降速率大于该阈值的作为持续变形区域。

研究区包含城市、裸岩、矿区等，单一的温度反演与地表形变方法易受地物类型的影响，利用Arcgis 软件对实验所提取的温度异常区域与沉降异常区域进行空间叠加，得到疑似火区范围。

3 试验结果分析

3.1 实测火区验证范围

利用测氡的实地勘测方法确定的煤火范围作为验证数据，来分析煤火识别方法的准确性。在地下火区燃烧过程中，燃烧区岩层及其上覆岩层处于高温高压环境中。煤系地层在高温高压作用下，氡的析出量不断增加[21]，此外，煤炭燃烧使煤系地层中孔隙水或裂隙水的温度和矿化度升高，导致氡的溶解度降低，使煤系地层中自由氡的数量进一步增加，这必然在火区上方地表浅层形成一个氡浓度高值区。测氡法虽然能较为准确地圈定煤火范围，但由于其易受气压、降水等因素影响，且监测规模有限，不适用于较大区域的火区识别。故试验选取部分区域进行实地测氡，用来验证融合方法的有效性。测网布置后在测区开展试验工作，经过现场调查，选择在南部测区寺耳沟村附近布置3 条测线，对所测的氡值剖面图进行分析，氡值显示出跳跃式变化，后进行质量检测，主要通过室外重复测量的方法，即在相同点位置重复布设活性炭吸附装置，在相对一致的地质条件和环境条件下埋置5 d，取出测量，将重复试验结果与初次试验结果进行比较，分析显示数据质量良好。沿煤层露头走向，地表岩层含有裂缝，且多处有高温热气流涌出，并伴有异常的刺激性气味。综合分析认定沿煤层露头线方向圈定了10 个区域为地下火区，以此作为试验的验证数据。

3.2 温度异常分析

煤火燃烧时会产生高温，高温产生的热量会以热辐射的形式向地表传导，在火区地表形成高于周围环境温度的温度异常区。通过一定的阈值，提取的温度异常区范围如图2 所示，其中添加了实测煤火范围作为验证。由图2 可以看到，绝大多数确定的火区范围均发生了温度异常，仅有1 处已确认的火区范围(9 号火区)没有提取出温度异常。但有部分产生温度异常的区域在测氡手段下并没有检测出煤火，这是由于地物自身的物理特性，例如砂岩、裸土等地物吸热以及城镇用地的热岛效应等导致出现了与煤火无关的温度异常，这也是单纯利用提取地表热异常信息来识别煤火的缺陷。

图2 温度异常图Fig.2 Temperature anomaly

3.3 地表形变异常分析

深部煤层长时间烧空后会造成地表塌陷，因此利用时序InSAR 的方法进行地表形变分析。采用密度分割法，提取SBAS-InSAR 形变结果中的持续变形区域。如图3 所示，在实测的10 处火区范围中，有7 处发生了明显的沉降，这在一定程度上证实了利用时序InSAR 技术分析地表形变从而识别煤田火区范围的可行性。但在实测范围中可以看到，有部分持续变形区域并不属于火区范围，可能是因为单纯的地质运动或开采活动，所以仅利用地表形变分析来识别火区也并不完全准确。

图3 地表沉降异常图Fig.3 Surface subsidence anomaly

综合分析所提取的温度异常区域与持续变形区域，在实测的10 个火区范围中，2、4、6、7、8、10 号火区均检测出了温度异常与持续变形信息。而1、3、5 号火区虽然显示沉降速率较小，但均出现了热异常，9 号火区发生了较大沉降，但温度仍在正常范围内。

3.4 疑似火区范围分析

将获取的温度异常区域与持续变形区域进行空间叠加分析，根据地表类型进行筛选，得到研究区域的疑似煤火区域。如图4 所示，共提取出11 个煤火区域，且采用笔者的研究方法提取的疑似火区与实测的火区在地理位置上具有较高的一致性。1～7 号疑似火区均位于实测煤火范围中，证明了研究提出的识别煤火方法的可行性。

图4 疑似火区范围Fig.4 Suspected fire area

利用实测的火区范围对试验结果进行对比验证。通过Arcgis 软件进行像元统计，得到各类方法的面积值，经统计，测氡法所得的实测火区面积Sm为183 663.222 672 m2，疑似火区面积为表中的融合法的面积，其值为43 806.258 4 m2，3 种方法的试验结果见表3。其中：

表3 试验结果对比Table 3 Comparison of experimental results

式中：Ra为准确率；R0为重叠率；S0为试验提取的火区范围与实测火区范围重叠的面积；Sr为试验提取的火区的面积；Sm为实测火区的面积。

可以看到，单一的温度异常和沉降异常提取方法的准确率较低，分别为50.49%和31.55%，而融合地表温度信息和沉降信息之后，火区识别的准确率大幅提高，高达93.78%，较单一方法提高了43.29%和62.23%，这大幅增加了煤火治理的有效性。但不足的是，有大量的火区仍未被识别，主要原因在于利用地表变形信息来识别煤火的研究较少，不够深入，缺乏长期的实地监测资料，难以获取较为准确的形变阈值，导致煤火识别准确性较低，进而影响融合结果。

对于实测范围外提取的火区，缺少实地验证数据，故采用影像对比分析的方法进行验证。如图5所示，从谷歌影像中可以看到，8、9 号疑似火区均在矿区内，且存在疑似煤火燃烧后形成的黑色区域，由此判定该区域确实存在煤火。

图5 谷歌地球中截取的部分疑似火区影像Fig.5 Google Earth footage of what appears to be a fire zone

4 结 论

1)协同热红外遥感技术和时序InSAR 技术提取煤田火区，克服了单一遥感技术的缺陷，显著提高了煤火识别的精度，为煤田火区的治理范围提供了有力参考。

2)试验结果表明，仅利用温度信息和沉降信息来识别煤火区域均存在一些缺陷。由于砂岩等地物吸热，城镇区热岛效应等原因，会提取出与煤火无关的温度异常范围。另外，由于矿区开采等造成的地表沉降也会影响利用沉降信息提取煤火的准确性。

3)由于缺乏长期实地的地表形变监测资料、煤火燃烧时间未知等，导致利用地表形变信息来识别煤火的方法准确性较低，今后的工作方向将集中于地表形变监测煤田火区的方法研究。