遙感方法在寧武煤田煤火識(shí)別中的應(yīng)用

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

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

0 引 言

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

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

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

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

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

1.1 研究區(qū)域概況

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

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

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

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

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

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

2 研究方法

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

圖1 研究流程Fig.1 Research process

2.1 溫度反演

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

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

式中：Ri為第i(i=10，11，···，14)波段的地表輻射亮度；Lgrd為包含大氣下行輻射的地表輻射亮度；εmax為初始最大發(fā)射率，取0.96；Latm↓為大氣下行輻射亮度；λi為波段i(i=10，11，···，14)的波長(zhǎng)；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ā)射率計(jì)算相對(duì)發(fā)射率值。

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

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

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

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

2.2 地表沉降

短基線集時(shí)序分析技術(shù)(SBAS-InSAR)是一種基于多幅SAR 影像的時(shí)間序列方法，克服了傳統(tǒng)D-InSAR 中存在的時(shí)間、空間失相關(guān)問(wèn)題，通過(guò)時(shí)空基線較短的干涉對(duì)提取地表形變信息[18]。

對(duì)于在tA，tB(A，B=0，1，···，N，A≠B)時(shí)刻生成的第K幅干涉圖，其任意像元的干涉相位值為[19-20]:

式中：Φ[tA，x，y]，Φ[tB，x，y]分別為tA，tB時(shí)刻相對(duì)于初始時(shí)刻t0的形變相位；x，y分別為方位向與距離向坐標(biāo);λ為雷達(dá)波長(zhǎng)；d(tA，x，y)和d(tB，x，y)為相對(duì)于初始時(shí)刻t0的視線方向的形變量。

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

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

3 試驗(yàn)結(jié)果分析

3.1 實(shí)測(cè)火區(qū)驗(yàn)證范圍

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

3.2 溫度異常分析

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

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

3.3 地表形變異常分析

深部煤層長(zhǎng)時(shí)間燒空后會(huì)造成地表塌陷，因此利用時(shí)序InSAR 的方法進(jìn)行地表形變分析。采用密度分割法，提取SBAS-InSAR 形變結(jié)果中的持續(xù)變形區(qū)域。如圖3 所示，在實(shí)測(cè)的10 處火區(qū)范圍中，有7 處發(fā)生了明顯的沉降，這在一定程度上證實(shí)了利用時(shí)序InSAR 技術(shù)分析地表形變從而識(shí)別煤田火區(qū)范圍的可行性。但在實(shí)測(cè)范圍中可以看到，有部分持續(xù)變形區(qū)域并不屬于火區(qū)范圍，可能是因?yàn)閱渭兊牡刭|(zhì)運(yùn)動(dòng)或開(kāi)采活動(dòng)，所以僅利用地表形變分析來(lái)識(shí)別火區(qū)也并不完全準(zhǔn)確。

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

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

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

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

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

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

表3 試驗(yàn)結(jié)果對(duì)比Table 3 Comparison of experimental results

式中：Ra為準(zhǔn)確率；R0為重疊率；S0為試驗(yàn)提取的火區(qū)范圍與實(shí)測(cè)火區(qū)范圍重疊的面積；Sr為試驗(yàn)提取的火區(qū)的面積；Sm為實(shí)測(cè)火區(qū)的面積。

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

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

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

4 結(jié) 論

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

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

3)由于缺乏長(zhǎng)期實(shí)地的地表形變監(jiān)測(cè)資料、煤火燃燒時(shí)間未知等，導(dǎo)致利用地表形變信息來(lái)識(shí)別煤火的方法準(zhǔn)確性較低，今后的工作方向?qū)⒓杏诘乇硇巫儽O(jiān)測(cè)煤田火區(qū)的方法研究。