自燃煤矸石山的遙感識別
——基于Landsat 8熱紅外波段地表溫度反演數據

李嘉琪 趙艷玲 任 河 肖 武

(1.中國礦業大學(北京)土地復墾與生態重建研究所，北京 100083;2.浙江大學公共管理學院，浙江 杭州 310058)

煤矸石是煤炭開采和加工過程中排放出的廢棄巖石，通常占采煤量的10%～15%，是我國最主要的工業固體廢棄物之一。目前煤矸石積存已超過80億t，并仍以年排放量1.5億t的速度增長[1]。煤矸石應用廣泛，可以用于道路填充、建筑材料生產等工業活動，但我國西部煤矸石的綜合利用程度很低，大部分煤矸石僅進行露天堆放[2-5]。由于煤矸石孔隙度較大，在松散的堆積下其內部殘煤和黃鐵礦極易與空氣中的氧氣和水接觸，發生低溫氧化，釋放大量熱量并誘發自燃。在燃燒過程中，大量有毒性氣體(如CO、SO2、NOX)排放，嚴重污染了大氣環境。此外，伴隨著降雨進行，會進一步對水環境和土壤環境造成污染。因此，對自燃煤矸石山的有效識別和快速監測對于礦區生態環境保護和礦區安全生產十分重要。

自燃煤矸石山通常會表現出地表熱異常，因此，地表溫度是識別和監測自燃煤矸石山的主要依據。基于熱紅外遙感技術能夠快速識別溫度信息，相比于原位測量表面溫度更經濟、更高效。陸地表面溫度(Land Surface Temperature，LST)是地表熱紅外輻射的綜合定量形式[6]，近年來在基于遙感數據反演地表溫度，進而識別礦區煤火和自燃煤矸石山領域也產生了許多研究成果。RAJU等[7]于2013年提出了一種利用波段Landsat TM/ETM+數據進行火災探測的最大太陽反射閾值方法，可以繪制出煤礦表面火災分布圖。HUO等[8]利用Landsat 7 ETM+影像的多時相夜間無云TIR波段數據反演地表溫度，煤表面火災大于150°C的溫度范圍很容易被監測到。2014年，NADUDVARI等[9]基于Landsat4-5 TM和ETM數據，并將熱圖與NDSI相結合提取自燃煤矸石山，證明圖上的熱異常與缺乏積雪覆蓋的區域有很好的相關性，并于2019年[10]利用Landsat ETM+影像獲取地表溫度圖和NDSI，探測了煤矸石山火點的分布、強度和演變等信息。可見，將 NDVI、NDSI、NDSI、熱圖與 Landsat TM、ETM+、ASTER 和 Sentinel 2圖像相結合，可以識別和定位煤矸石山中的熱點。以往研究主要圍繞地表溫度反演算法、閾值提取方法和高溫區域提取方法，大多利用地溫結果，或者將地溫結果和植被參數等遙感參量相結合提取自燃煤矸石山信息，忽略了夜晚同樣高溫的地物干擾，例如建筑物、工廠、道路等地物。

本研究在已有成果的基礎上，將地溫反演結果和監督分類結果相結合，提取出裸地和煤矸石山地類，去除其他高溫發熱地物的影響，利用閾值分析法確定裸地和煤矸石山地類中的高溫區域，最終識別出自燃煤矸石山區域。研究結果可為礦區大尺度范圍內精確識別自燃煤矸石山提供基礎數據和方法參考。

1 研究區和數據來源

1.1 研究區概況

陽泉市地處山西省中部東翼，太行山中段西側，位于東經 112°5′～114°4′，北緯 37°40′～38°31′區域，總面積4 452 km2，占山西省總面積的3%，北接忻州市，東臨河北石家莊市，西南方連太原市和晉中市。地處黃土高原東緣，地貌類型以山地為主，其余為丘陵和平原，地勢西高東低，區內地形復雜，山巒起伏，溝壑縱橫。本區屬暖溫帶半濕潤大陸性季風氣候，四季分明，日照充足，晝夜溫差較大，風向多變，年降水量為450～550 mm，主要集中在每年6—9月，年均氣溫一般為8～12℃。

陽泉市下轄城區、礦區和郊區3個市轄區和盂縣、平定縣2個縣。陽泉地區礦產資源豐富，開發歷史悠久，就煤炭資源來看，陽泉市位于沁水煤田東北隅，境內石炭紀、二疊紀煤系分布于區域西南角，煤炭資源以高變質無煙煤及貧煤為主[11]。

陽泉市境內含煤面積1 051 km2，煤炭地質儲量104億t，原煤年生產量350萬t。巨大的煤炭開采量產生了大量的煤矸石。據統計，截至2014年，陽泉市矸石累計堆積量約2.8億t[12]，約占全國的4%，并以900萬t/a的速度增長，已知自燃矸石山42座，絕大多數矸石山未采取綠化措施。由于陽泉境內有眾多的矸石山，代表性較強，基于陽泉境內42座自燃煤矸石山的基礎數據，選擇陽泉市作為研究區(圖1)。

圖1 研究區地理位置和行政區劃概況Fig.1 Overview of the geographical location and administrative division of the study area

1.2 數據來源和預處理

目前常用的載有熱紅外波段的衛星傳感器包括:NOAA搭載的 AVHRR3，TERRA搭載的 MODIS，CBERS 搭載的 IRMSS，Landsat5/7搭載的TM/ETM+，Landsat 8搭載的 TIRS，FY-2和 FY-3搭載的 VIRR，HJ-1B搭載的IRS等(表1)。其中，Landsat 8因具有更高的分辨率從而能更有效地分辨小地塊中不同地類的地溫差異，在精確反演地表溫度方面優勢明顯，應用較廣泛[13]。

表1 攜帶熱紅外通道傳感器衛星的熱紅外波段比較Table 1 Comparison of thermal infrared bands of satellites carrying thermal infrared channel sensors

Landsat 8衛星有可見光(OLI)和熱紅外(TIRS)兩個傳感器，可見光傳感器的空間分辨率為30 m(波段1～7和波段9)和15 m(全色波段8)，熱紅外傳感器的分辨率為100 m(波段10～11)。Landsat 8衛星于2013年2月發射升空，與之前發射的同系列衛星相比，不僅繼承了熱紅外通道高空間分辨率的特征，而且在其搭載的熱紅外傳感器上首次引入了量子學紅外光子檢測技術，增強了TIRS對熱紅外波段的敏感性，從而更加有利于從大氣溫度中將地表溫度分離出來。因而與ETM+/TM6相比，Landsat 8更加適合于用來對地表溫度、熱空間分布進行精確分析[14]。

陽泉市Landsat 8影像通過美國地質調查局(USGS)網站下載，3景影像覆蓋研究區全境，軌道號分別為 124-34、125-33、125-34，云量均低于 10%，獲取時間為2018年12月18日和2019年7月30日。陽泉市和北京市都位于東八區時區，沒有時差，故衛星過境當地時間即為北京時間，均為凌晨3時左右。Landsat8數據使用UTM/WGS84投影坐標系，數據產品Level1T已進行了基于地形的幾何校正，在ENVI5.3軟件中對TIRS與OLI數據進行輻射定標，然后利用FLAASH模型對OLI數據進行大氣校正，得到的多光譜影像和熱紅外影像利用陽泉市邊界線的矢量數據進行鑲嵌和裁剪。

反演地表溫度數據需要衛星過境時研究區域盡量精確的氣象數據，包括大氣上行、下行輻射亮度以及大氣透射率等。在NASA的官網中(網址為http:∥atmcorr.gsfc.nasa.gov/)根據影像拍攝時間和影像中心的經緯度可以確定相應的大氣透過率、大氣上行和下行輻射亮度。

其他數據包括在山西省測繪地理信息局官網下載得到的山西省市級、縣級和河流的矢量數據，Google Earth提供的研究區正射影像圖，現有的自燃矸石山點位坐標及其相關信息由陽泉煤業集團提供。

2 研究方法

2.1 技術路線

本研究基于Landsat 8熱紅外波段數據反演陽泉市地表溫度，進而識別自燃煤矸石山。為便于進行冬夏季對比，消除季節的影響，選取研究區2018年冬季和2019年夏季兩天夜晚的影像進行分析。將反演的地表溫度圖與監督分類結果相結合，得到裸地和煤矸石山地類的地表溫度分布圖。基于裸地和煤矸石山地類地表溫度分布的灰度值直方圖確定一個高溫閾值，根據閾值提取出的高溫區域來判讀精度，分析已知的42座自燃煤矸石山被識別進入高溫區域的比例。最后分析識別出的自燃煤矸石山周邊地表溫度分布趨勢以及自燃煤矸石山高溫影響范圍與邊界(圖 2)。

圖2 研究技術路線示意Fig.2 Diagram of technical route of the study

2.2 大氣校正法

目前，地表溫度的遙感反演算法主要有輻射傳輸方程法、單窗算法、分裂窗算法、多通道多角度算法等，不同算法適用于不同遙感傳感器的熱紅外數據。

Landsat8熱紅外波段為第10、11波段，因而它既能像MODIS/EOSS一樣采用劈窗算法反演地表溫度，又可采用單波段算法反演得到地表溫度。由于TIRS10波段相對于TIRS11波段而言位于一個較低的大氣吸收區域，具有更高精度的大氣透過率值[15-16]，并且現階段TIRS11波段定標數據仍不理想，故本研究利用TIRS10波段反演地表溫度。當前應用單熱紅外波段獲得的輻射能量反演LST常見的有3種算法，即輻射傳輸方程法、單窗算法和單通道算法。無論是單窗算法還是單通道算法，都是從輻射傳輸方程法出發對Planck函數進行簡化，從熱紅外傳感器獲得地物的亮度、溫度，從而推導出適用于單個波段數據的地表溫度。輻射傳輸方程法即大氣校正法的基本理論依據是普朗克定律和維恩位移定律[17]，是傳統地表溫度反演用途最廣泛、精度較高的方法之一，因此本研究采用該方法反演陽泉市地表溫度。

大氣校正法的基本原理是首先估計大氣對地表熱輻射的影響，然后將這部分大氣影響從衛星傳感器所觀測到的熱輻射總量中減去，從而得到地表熱輻射強度，再將這一熱輻射強度轉化為相應的地表溫度。衛星傳感器接收到的熱紅外輻射亮度值Lλ由3個部分組成，即大氣向上輻射亮度L↑，地面的真實輻射亮度經過大氣層之后到達衛星傳感器的能量Ts，大氣向下輻射到達地面后反射的能量L↓。輻射傳輸方程為

式中，Lλ為傳感器接收到的熱紅外輻射亮度值，W·(m2·sr·μm);ε為地表比輻射率;Ts為地表真實溫度，K;B(Ts)為黑體熱輻射亮度，W·(m2·sr·μm);τ為大氣在熱紅外波段的透過率。

地表比輻射率是同一溫度下地表發射的輻射量與黑體發射輻射量的比值，是利用大氣校正法進行溫度反演時的一個重要參數，與地表組成成分、地表粗糙度、波長等因素有關。本研究利用覃志豪等[18]提出的混合像元分離法(HHXY)計算地表比輻射率。求取該參數之前需要解算歸一化植被指數NDVI和植被覆蓋度Pv。NDVI可進行如下計算:

式中，NIR和R分別為近紅外波段和紅光波段的反射率。

Pv采用混合像元分解法計算，將整景影像的地類大致分為水體、植被和建筑，計算公式為

式中，通過對NDVI影像進行概率統計分析得到，2018年取NDVImax=0.992 1和NDVImin=-0.003 9。2019年取NDVImax=0.930 1和NDVImin=-0.001 1。當某個像元的NDVI＞NDVImax時，Pv取 1;當NDVI＜NDVImin時，Pv取 0。

本研究將遙感影像分為水體、城鎮和自然表面3種類型，采用基于土地覆蓋類型方法計算地表比輻射率。水體像元的比輻射率賦值為0.995，自然表面和城鎮像元的比輻射率可進行如下計算:

根據式(1)，溫度為T的黑體在熱紅外波段的輻射亮度B(Ts)為

Ts可以用普朗克公式的反函數獲取，

對于TIRS波段10，K1=774.89 W/(m2·μm·sr)，K2=1 321.08 K。

3 結果與分析

3.1 陽泉市地表溫度反演結果

通過上述研究方法在ENVI5.3軟件中計算處理得到2018年12月18日和2019年7月30日兩天夜晚的地表溫度反演圖，如圖3所示。兩期衛星遙感影像均為凌晨3時左右獲取，冬季全區地溫范圍為-4.33～22.12℃，平均地溫為8.56℃;夏季全區地溫范圍為16.68～43.92℃，平均地溫為29.39℃。在陽泉市城區、盂縣縣城和平定縣縣城等一些居民活動區地表溫度明顯較高，市區中間一條線狀低溫區為桃河。已知的42座自燃矸石山點位在地表溫度圖上的分布如圖3所示，可以大致判斷已知的自燃煤矸石山點位大多位于圖中的高溫區域。由于缺乏實地觀測地表溫度數據，本研究無法計算反演結果精度，主要對得到的研究區地表溫度圖進行高溫異常區域分析。

圖3 陽泉市地表溫度反演結果Fig.3 Inversion results of land surface temperature in Yangquan City

3.2 陽泉市土地利用類型分類結果

對研究區范圍內進行監督分類以提取煤矸石山區域，首先對Landsat 8影像進行波段融合，將低分辨率多光譜影像(分辨率30 m)與高分辨率全色波段影像(分辨率15 m)融合生成一幅高分辨率多光譜影像(分辨率為15 m)，處理后得到的影像既有較高的空間分辨率，又具有多光譜特征。波段融合方法有IHS變換、Brovey變換、乘積運算(CN)、PCA變換和Gram-schmidt Pan Sharpening(GS)方法，其中 GS法是通過統計分析方法對參與融合的各波段進行最佳匹配，改進了傳統融合方法中某些波段信息過分集中的不足，不受波段限制，可有效保持空間紋理信息，尤其能高保真保持光譜特征。

利用GS法進行波段融合，對獲得的融合影像利用支持向量機法進行監督分類，結果如圖4所示。土地利用類型分為4類，即植被、水體、建設用地及裸地和煤矸石山。裸地和煤矸石山地類中也包括工礦用地中被煤粉覆蓋的區域，該類地類的地表紋理特征與裸地相近。對監督分類得到的高空間分辨率影像進行重采樣，采用最鄰近內插法不會更改輸入像元的值，適用于離散(分類)數據的重采樣。設置重采樣的土地利用類型影像像元大小為30 m，與地表溫度圖的分辨率相同。將重采樣后得到的裸地和煤矸石山一類與地表溫度圖相結合，得到裸地和煤矸石山區域的地表溫度分布，如圖5所示。

圖4 可見光影像監督分類結果Fig.4 Supervised classification results of visible light image

圖5 裸地和煤矸石山地表溫度Fig.5 Land surface temperature of bare land and gangue dumps

3.3 基于地表溫度異常的自燃煤矸石山提取結果

圖5中裸地和煤矸石山區域的地表溫度分布灰度直方圖如圖6所示。灰度值呈正態分布，根據直方圖統計法[19]，確定峰值右側一倍標準差的位置對應的橫坐標溫度為高溫閾值(圖7)。夏季直方圖峰值對應的溫度為30.55℃，標準差σ為1.859 1，確定夏季高溫閾值為32.409 1℃。冬季直方圖峰值對應的溫度為9.55℃，標準差σ為2.203 3，確定冬季高溫閾值為11.753 3℃。根據閾值提取出裸地和煤矸石山地類中的高溫異常區域，將提取出的高溫區域疊加到Google Earth影像中(圖8)，利用已知的42座自燃矸石山區域進行驗證，結果顯示，2018年冬季有32座顯示出了高溫異常，提取精度為76.19%;2019年夏季有27座顯示出了高溫異常，提取精度為64.28%。

圖6 裸地和煤矸石山地表溫度圖的灰度直方圖Fig.6 Grayscale histograms of land surface temperature maps of bare land and gangue dumps

圖7 高溫閾值提取示意Fig.7 Schematic of high temperature threshold extraction

圖8 裸地和煤矸石山地類中的熱區與Google Earth疊加顯示Fig.8 Overlay display of hot areas in bare land and gangue dumps with Google Earth orthophoto

最終提取出的熱區中存在一些錯分現象，部分非自燃矸石山區域被提取出來，這些區域大多是一些發熱地物周圍的裸地，如電廠、焦化廠、蔬菜大棚周邊的裸地等，一些礦區工礦用地、廠房附近的裸地和煤粉覆蓋的地面溫度也較高。這些地表溫度較高的區域影響了識別精度，需通過Google Earth影像進行識別和去除。去除上述干擾因素后，可以明顯識別出自燃煤矸石山，證明通過Landsat 8熱紅外波段反演地表溫度進而識別自燃煤矸石區域的方法是可行的。

識別結果也存在一定的漏分現象，部分已知的自燃煤矸石山未能被識別成高溫區域，漏分原因可能是:①部分矸石山建成時間較長已被復墾，表面綠植覆蓋情況較好，無裸露地表，在監督分類中被識別為植被。②矸石山自燃范圍小于影像單個像素面積，自燃點無法在影像上顯示;或者采取了重新覆土等煤矸石山復墾措施，復墾效果較好，地表已沒有較大的溫度異常。

通過冬夏季高溫異常區域對比，可以發現2019年比2018年自燃矸石山減少了5座，已經自燃的矸石山高溫異常面積有所減少，說明大范圍的熱紅外遙感反演地表溫度有助于實時、準確地分析自燃煤矸石山狀態(初期、中期和晚期)及發展趨勢，可為準確評估煤矸石山造成的經濟損失和周邊生態環境影響提供相對精確的熱力學指標[20-23]。

3.4 自燃煤矸石山影響范圍劃定

根據上述獲取的自燃煤矸石山的點位信息提取出點位附近100～1 000 m范圍內每100 m的平均地表溫度，如圖9和圖10所示。煤矸石山不同距離范圍內的地表溫度均值變化表明，矸石山區域的地表溫度比區域整體環境地表溫度高，隨著與矸石山區域距離的增加，平均地表溫度逐漸降低，降低幅度逐漸減小。從變化曲線來看，2018年冬季矸石山周邊500 m空間范圍平均地表溫度變化較大，500 m以外平均地表溫度變化較小，逐漸趨于平緩，2019年夏季矸石山周邊600 m空間范圍平均地表溫度變化較大，600 m以外趨于平緩。總體來說，煤矸石山對周邊地區地表溫度的影響距離超過了500 m。

圖9 典型自燃煤矸石山多環緩沖區示意Fig.9 Schematic of muliple ring buffer of typical spontaneous combustion gangue dumps

圖10 不同距離分區地表溫度均值變化Fig.10 Variation of mean land surface temperature in different distance zones

統計了冬季提取出的32座自燃矸石山周邊1 000 m緩沖區范圍內的土地利用情況，將植被和水體面積占比劃分為3個等級，大于50%、20%～50%和小于20%，分別是18座、8座和6座，統計出3類矸石山隨著距離的增加周邊地表溫度的變化特征，如圖11所示。由圖11可知:3類自燃煤矸石山對周邊地區地表溫度的影響距離均超過了500 m，但植被和水體占比超過50%的自燃矸石山500 m范圍內降溫最快，緩沖區內整體地溫最低。隨著植被和水體面積占比逐漸減少，矸石山周邊降溫速度逐漸減緩，緩沖區內整體地溫逐漸升高，說明植被和水體對自燃矸石山地表有降溫作用，周邊地物分布會對自燃矸石山地表高溫影響范圍和強度產生影響[24-27]。

圖11 冬季矸石山周邊不同植被和水體面積占比的地表溫度變化Fig.11 Variation of land surface temperature of different vegetation and water proportion around gangue dumps in winter

4 結 論

本研究提出了一種基于熱紅外遙感數據反演地表溫度進而識別自燃煤矸石山的方法。基于Landsat 8遙感影像，利用大氣校正法反演得到山西省陽泉市2018年和2019年兩天夜晚的地表溫度，將地表溫度反演結果與監督分類結果相結合，去除了夜晚高溫的其他城市地物信息，在此基礎上利用閾值分析法提取出裸地和煤矸石山中的高溫區域。得出如下結論:

(1)基于Landsat 8熱紅外數據進行地表溫度反演，結合監督分類提取出矸石山地類區，利用閾值分析法提取矸石山地類區中的高溫區域，進而識別自燃煤矸石山的方法可行且精度較高。

(2)自燃煤矸石山對周邊地區地表溫度的影響距離超過500 m，植被和水體對自燃矸石山地表具有降溫作用。

(3)本研究識別結果仍存在一定的錯分漏分現象，部分自燃煤矸石山表面高溫范圍小于熱紅外影像像元大小，存在混合像元問題出現漏分自燃煤矸石山現象，進而影響了識別精度。后續工作中，可以考慮熱紅外波段數據與高空間分辨率的可見光/熱紅外數據融合，利用地表像元分解方法獲取空間分辨率較高的地表溫度信息。受到高溫裸地和工礦用地影響，監督分類方法精度仍有一定的提高空間，后續工作中可以結合紋理信息與光譜信息進行監督分類，提高煤矸石山地類的分類精度。