高玉榮 ，隋 剛 ，張新軍 ，孔嘉嫄 ，張和生

（太原理工大學 礦業(yè)工程學院, 山西 太原 030024）

0 引 言

山西省煤炭資源豐富,2022 年全省規(guī)模以上原煤產(chǎn)量為130 714.6 萬t，同比增長8.7%。大規(guī)模開采帶來嚴重的生態(tài)問題，其中煤的自燃危害尤其嚴重。煤火燃燒不僅造成巨大的能源浪費和經(jīng)濟損失[1-2]，還會產(chǎn)生大量有害氣體，危害人民生命健康，同時也會引起地表裂縫和塌陷，對當?shù)氐幕A設施產(chǎn)生影響。因此，準確識別煤田火區(qū)范圍對于煤火監(jiān)測與治理具有重要意義。

近年來，大量學者對新疆、內蒙古、寧夏等地較大范圍煤田火區(qū)進行了探測與研究。自20 世紀60年代以來，國內外已經(jīng)發(fā)展了多種煤田火區(qū)探測技術，主要包括磁探法、物探法、化探法和遙感法等[3]。其中利用煤層燃燒的物理、化學特性進行探測的方法由于技術限制、成本較高，均不適用于大面積火區(qū)[4]。而遙感技術具有探測范圍大、獲取周期短、時效性強、經(jīng)濟效益高等優(yōu)勢，已成為煤田火區(qū)識別領域的重要發(fā)展方向[5]。

利用煤火燃燒產(chǎn)生熱量，傳導至地表造成溫度異常這一特性，學者們通過地表溫度反演提取熱異常信息對煤火進行識別和監(jiān)測。地表溫度反演算法總體分為單通道法、分裂窗法、日夜法、TES、分裂窗和溫度比輻射率分離結合法[6]。邱程錦等[7]通過大氣校正的反演方法對Landsat TM/ETM 進行處理，并提取溫度異常區(qū)域來判定煤火范圍；李峰等[8]利用TES 算法反演內蒙古烏達礦區(qū)4 個時期的地表溫度，并采用自適應梯度閾值法提取對應火區(qū)范圍來監(jiān)測并評估烏達煤火的治理效果。但深部煤火產(chǎn)生的熱量可能傳導不到地面，易造成煤火區(qū)漏判，另外砂巖吸熱、城市熱島效應等因素也會導致與煤層燃燒無關的高溫異常區(qū)[9]，從而造成誤判，所以僅利用熱紅外遙感提取熱異常區(qū)域的方法識別煤田火區(qū)范圍存在一定缺陷。針對煤火燃燒會造成地下空洞，引起地表塌陷變形等特性，少部分學者利用In-SAR 技術進行煤火識別和監(jiān)測。如JIANG L 等[10]通過PS-InSAR、Stacking 和D-InSAR 三種合成孔徑雷達方法對烏達煤田進行了監(jiān)測，證實了利用In-SAR 手段進行地表形變分析，可以有效判斷地下煤火的燃燒情況；RIYAS 等[11]為了表征印度Jharia 煤礦火災的時空動態(tài)，利用新小基線子集(N-SBAS)技術計算了該煤田2017-2020 年地表變形時間序列，但煤礦開采[12]、地質災害等也會引起地表形變，從而對煤火識別產(chǎn)生影響。

綜上所述，現(xiàn)有方法雖然在煤火識別方面取得了較多研究成果，但識別特性單一，易受地物吸熱、煤礦開采等因素的影響。因此采用融合熱異常信息和地表形變信息來識別煤田火區(qū)可以有效克服單一方法的不足，提高煤火識別準確性。筆者利用多源遙感手段，通過融合衛(wèi)星熱紅外技術與雷達技術，首次對山西省寧武煤火進行識別分析研究。

1 研究區(qū)域概況與數(shù)據(jù)源

1.1 研究區(qū)域概況

研究區(qū)位于寧武縣東寨鎮(zhèn)境內，地理坐標為東經(jīng)112°03′57″～112°12′39″，北緯38°46′41″～38°53′47″。區(qū)內地勢高峻，山嶺縱橫，海拔在2 000 m左右，地形總體西高東低，最高處位于西部，海拔2 020 m，最低處位于東部深溝底，海拔1 670 m，最大高差為350 m。屬溫帶大陸性氣候，寒冷多大風，晝夜溫差較大，年平均氣溫7 ℃。研究區(qū)域內的煤礦主要分布在中部，呈南北走向，包含寺耳溝煤礦，小西溝煤礦，車道溝煤礦，三馬營煤礦等多個煤礦區(qū)。由于煤礦開采引發(fā)了地面塌陷、地裂縫、最終形成火風壓，造成了多個煤火自燃區(qū)域[13]，加之小煤窯不規(guī)范化的開采方式，產(chǎn)生的煤火問題未能得到及時解決。

1.2 數(shù)據(jù)源

研究采用2019-08-14 在寧武過境的Terra 衛(wèi)星上的ASTER(Advanced Spaceborne Thermal Emission and Reflection Radiometer)傳感器獲取的夜間遙感數(shù)據(jù)。由于晚上太陽輻射消失，所獲取的熱紅外遙感影像能更準確地提取熱異常信息[14]。衛(wèi)星過境時間為14:34:12(UTC)，云覆蓋量為2%，數(shù)據(jù)產(chǎn)品級別為 L1T，軌道號為 222/211，時相和云覆蓋量滿足試驗要求。所用的ASTER 數(shù)據(jù)參數(shù)見表1。

表1 熱紅外數(shù)據(jù)參數(shù)Table 1 Thermal infrared data parameters

采用的雷達影像數(shù)據(jù)為Sentinel-1A 數(shù)據(jù)，用于煤田火區(qū)的地表變形監(jiān)測。試驗選用工作模式為IW、極化方式為VV 的SLC 數(shù)據(jù)。選取2018 年6月至2020 年1 月的26 景影像(SAR 數(shù)據(jù)成像時間見表2)用于差分干涉處理。外部參考DEM 選用空間分辨率為12.5 m 的ALOS PALSAR 數(shù)據(jù)。

表2 SAR 影像成像時間統(tǒng)計Table 2 Statistics of SAR imaging time

2 研究方法

以山西寧武部分煤田為研究區(qū)域，ASTER 夜間熱紅外影像、Sentinel-1A 影像為試驗數(shù)據(jù)，采用多源遙感融合的方法進行煤火區(qū)域識別。首先采用溫度比輻射率分離算法(ASTER-TES)進行地表溫度反演，通過一定的閾值提取研究區(qū)的熱異常范圍，結合SBAS-InSAR 技術提取地表持續(xù)變形信息，將熱異常信息與地表形變信息空間疊加分析后得到研究區(qū)內煤火的疑似分布區(qū)域，再根據(jù)部分實測煤火范圍進行驗證分析。研究流程如圖1 所示。

圖1 研究流程Fig.1 Research process

2.1 溫度反演

溫度反演中2 個至關重要的因素是溫度和地表比輻射率。傳統(tǒng)的溫度反演算法常常假設比輻射率已知，來求解地表溫度，而TES 算法則根據(jù)一定的先驗知識作為約束條件，同時求解比輻射率和溫度。由于ASTER TES 算法是針對ASTER 數(shù)據(jù)反演的官方算法[15]，充分利用了ASTER 數(shù)據(jù)的5 個熱紅外波段，且吸收了發(fā)射率歸一化NEM、光譜比值RATIO、和最大最小發(fā)射率差值MMD 3 個模塊的優(yōu)點，精度較高，故采用此算法進行地表溫度反演，其包括3 個部分:

1)NEM 模塊：初步估算目標表面溫度并且從輻射亮度觀測中減去反射的大氣輻射。

式中：Ri為第i(i=10，11，···，14)波段的地表輻射亮度；Lgrd為包含大氣下行輻射的地表輻射亮度；εmax為初始最大發(fā)射率，取0.96；Latm↓為大氣下行輻射亮度；λi為波段i(i=10，11，···，14)的波長；TNEM為NEM模塊的輸出溫度；c1為第一輻射常數(shù)，c2為第二輻射常數(shù)，其中h為普朗克常數(shù)，6.626 176×10-34J·s；c為真空光速，2.997 924 58×108m/s；k為玻爾茲曼常數(shù)，1.380 6×10-23J/k。

2) RATIO 模塊：利用NEM 模塊估算的發(fā)射率計算相對發(fā)射率值。

式中：βi為第i(i=10，11，···，14)波段的相對比輻射率。

3) MMD 模塊：進一步估算發(fā)射率和溫度。

其中：εmin為發(fā)射率最小值；max 為發(fā)射率εi最大值(εmax)所對應的波段。當MMD＜0.03 時，灰體的精度很低，不再使用MMD 方法，此時直接將εmin設為0.983。數(shù)值根據(jù)水體和濃密植被的性質確定[16]。

對比發(fā)現(xiàn)，熱紅外數(shù)據(jù)與雷達數(shù)據(jù)在地理位置上存在偏差，但雷達數(shù)據(jù)難以選取控制點，故以與雷達數(shù)據(jù)空間位置較為匹配的Landsat 影像為基準影像，利用Arcgis 對ASTER 數(shù)據(jù)進行地理配準。隨后采用ENVI 軟件對地理配準后的ASTER 數(shù)據(jù)進行輻射定標，大氣校正處理，并基于IDL8.7 平臺實現(xiàn)該算法。重復式(1)—式(9)，直到迭代計算的相鄰2 次的溫差＜0.3 K 或迭代次數(shù)＞12 次為止，得到相應的地表溫度值，通過一定的閾值提取溫度異常區(qū)域。經(jīng)統(tǒng)計，溫度反演的結果符合正態(tài)分布，將(μ+2σ)作為高溫閾值進行地表熱異常的提取，其中，μ為溫度反演結果中統(tǒng)計值的數(shù)學期望，σ為統(tǒng)計值的標準差[17]。

2.2 地表沉降

短基線集時序分析技術(SBAS-InSAR)是一種基于多幅SAR 影像的時間序列方法，克服了傳統(tǒng)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的視線方向的形變量。

由于研究區(qū)域內植被較茂密，為避免完全空間失相關，將空間基線閾值和時間基線閾值分別設置為45 m 和365 d，對輸入的SAR 數(shù)據(jù)進行干涉像對的配對，并對部分像對進行3D 解纏；采用Delaunay MCF 的方法進行解纏，Goldstein 進行濾波；選擇合適的GCP，估算和去除殘余的恒定相位和經(jīng)解纏后還存在的相位坡道;經(jīng)過兩次反演得到形變速率結果，并將形變速率結果從斜距投影轉換為地理投影。由于研究區(qū)域的煤火初始燃燒時間未知，缺乏地表監(jiān)測資料，根據(jù)相關文獻，將形變閾值確定為5 mm/a，沉降速率大于該閾值的作為持續(xù)變形區(qū)域。

研究區(qū)包含城市、裸巖、礦區(qū)等，單一的溫度反演與地表形變方法易受地物類型的影響，利用Arcgis 軟件對實驗所提取的溫度異常區(qū)域與沉降異常區(qū)域進行空間疊加，得到疑似火區(qū)范圍。

3 試驗結果分析

3.1 實測火區(qū)驗證范圍

利用測氡的實地勘測方法確定的煤火范圍作為驗證數(shù)據(jù)，來分析煤火識別方法的準確性。在地下火區(qū)燃燒過程中，燃燒區(qū)巖層及其上覆巖層處于高溫高壓環(huán)境中。煤系地層在高溫高壓作用下，氡的析出量不斷增加[21]，此外，煤炭燃燒使煤系地層中孔隙水或裂隙水的溫度和礦化度升高，導致氡的溶解度降低，使煤系地層中自由氡的數(shù)量進一步增加，這必然在火區(qū)上方地表淺層形成一個氡濃度高值區(qū)。測氡法雖然能較為準確地圈定煤火范圍，但由于其易受氣壓、降水等因素影響，且監(jiān)測規(guī)模有限，不適用于較大區(qū)域的火區(qū)識別。故試驗選取部分區(qū)域進行實地測氡，用來驗證融合方法的有效性。測網(wǎng)布置后在測區(qū)開展試驗工作，經(jīng)過現(xiàn)場調查，選擇在南部測區(qū)寺耳溝村附近布置3 條測線，對所測的氡值剖面圖進行分析，氡值顯示出跳躍式變化，后進行質量檢測，主要通過室外重復測量的方法，即在相同點位置重復布設活性炭吸附裝置，在相對一致的地質條件和環(huán)境條件下埋置5 d，取出測量，將重復試驗結果與初次試驗結果進行比較，分析顯示數(shù)據(jù)質量良好。沿煤層露頭走向，地表巖層含有裂縫，且多處有高溫熱氣流涌出，并伴有異常的刺激性氣味。綜合分析認定沿煤層露頭線方向圈定了10 個區(qū)域為地下火區(qū)，以此作為試驗的驗證數(shù)據(jù)。

3.2 溫度異常分析

煤火燃燒時會產(chǎn)生高溫，高溫產(chǎn)生的熱量會以熱輻射的形式向地表傳導，在火區(qū)地表形成高于周圍環(huán)境溫度的溫度異常區(qū)。通過一定的閾值，提取的溫度異常區(qū)范圍如圖2 所示，其中添加了實測煤火范圍作為驗證。由圖2 可以看到，絕大多數(shù)確定的火區(qū)范圍均發(fā)生了溫度異常，僅有1 處已確認的火區(qū)范圍(9 號火區(qū))沒有提取出溫度異常。但有部分產(chǎn)生溫度異常的區(qū)域在測氡手段下并沒有檢測出煤火，這是由于地物自身的物理特性，例如砂巖、裸土等地物吸熱以及城鎮(zhèn)用地的熱島效應等導致出現(xiàn)了與煤火無關的溫度異常，這也是單純利用提取地表熱異常信息來識別煤火的缺陷。

圖2 溫度異常圖Fig.2 Temperature anomaly

3.3 地表形變異常分析

深部煤層長時間燒空后會造成地表塌陷，因此利用時序InSAR 的方法進行地表形變分析。采用密度分割法，提取SBAS-InSAR 形變結果中的持續(xù)變形區(qū)域。如圖3 所示，在實測的10 處火區(qū)范圍中，有7 處發(fā)生了明顯的沉降，這在一定程度上證實了利用時序InSAR 技術分析地表形變從而識別煤田火區(qū)范圍的可行性。但在實測范圍中可以看到，有部分持續(xù)變形區(qū)域并不屬于火區(qū)范圍，可能是因為單純的地質運動或開采活動，所以僅利用地表形變分析來識別火區(qū)也并不完全準確。

圖3 地表沉降異常圖Fig.3 Surface subsidence anomaly

綜合分析所提取的溫度異常區(qū)域與持續(xù)變形區(qū)域，在實測的10 個火區(qū)范圍中，2、4、6、7、8、10 號火區(qū)均檢測出了溫度異常與持續(xù)變形信息。而1、3、5 號火區(qū)雖然顯示沉降速率較小，但均出現(xiàn)了熱異常，9 號火區(qū)發(fā)生了較大沉降，但溫度仍在正常范圍內。

3.4 疑似火區(qū)范圍分析

將獲取的溫度異常區(qū)域與持續(xù)變形區(qū)域進行空間疊加分析，根據(jù)地表類型進行篩選，得到研究區(qū)域的疑似煤火區(qū)域。如圖4 所示，共提取出11 個煤火區(qū)域，且采用筆者的研究方法提取的疑似火區(qū)與實測的火區(qū)在地理位置上具有較高的一致性。1～7 號疑似火區(qū)均位于實測煤火范圍中，證明了研究提出的識別煤火方法的可行性。

圖4 疑似火區(qū)范圍Fig.4 Suspected fire area

利用實測的火區(qū)范圍對試驗結果進行對比驗證。通過Arcgis 軟件進行像元統(tǒng)計，得到各類方法的面積值，經(jīng)統(tǒng)計，測氡法所得的實測火區(qū)面積Sm為183 663.222 672 m2，疑似火區(qū)面積為表中的融合法的面積，其值為43 806.258 4 m2，3 種方法的試驗結果見表3。其中：

表3 試驗結果對比Table 3 Comparison of experimental results

式中：Ra為準確率；R0為重疊率；S0為試驗提取的火區(qū)范圍與實測火區(qū)范圍重疊的面積；Sr為試驗提取的火區(qū)的面積；Sm為實測火區(qū)的面積。

可以看到，單一的溫度異常和沉降異常提取方法的準確率較低，分別為50.49%和31.55%，而融合地表溫度信息和沉降信息之后，火區(qū)識別的準確率大幅提高，高達93.78%，較單一方法提高了43.29%和62.23%，這大幅增加了煤火治理的有效性。但不足的是，有大量的火區(qū)仍未被識別，主要原因在于利用地表變形信息來識別煤火的研究較少，不夠深入，缺乏長期的實地監(jiān)測資料，難以獲取較為準確的形變閾值，導致煤火識別準確性較低，進而影響融合結果。

對于實測范圍外提取的火區(qū)，缺少實地驗證數(shù)據(jù)，故采用影像對比分析的方法進行驗證。如圖5所示，從谷歌影像中可以看到，8、9 號疑似火區(qū)均在礦區(qū)內，且存在疑似煤火燃燒后形成的黑色區(qū)域，由此判定該區(qū)域確實存在煤火。

圖5 谷歌地球中截取的部分疑似火區(qū)影像Fig.5 Google Earth footage of what appears to be a fire zone

4 結 論

1)協(xié)同熱紅外遙感技術和時序InSAR 技術提取煤田火區(qū)，克服了單一遙感技術的缺陷，顯著提高了煤火識別的精度，為煤田火區(qū)的治理范圍提供了有力參考。

2)試驗結果表明，僅利用溫度信息和沉降信息來識別煤火區(qū)域均存在一些缺陷。由于砂巖等地物吸熱，城鎮(zhèn)區(qū)熱島效應等原因，會提取出與煤火無關的溫度異常范圍。另外，由于礦區(qū)開采等造成的地表沉降也會影響利用沉降信息提取煤火的準確性。

3)由于缺乏長期實地的地表形變監(jiān)測資料、煤火燃燒時間未知等，導致利用地表形變信息來識別煤火的方法準確性較低，今后的工作方向將集中于地表形變監(jiān)測煤田火區(qū)的方法研究。