準(zhǔn)南煤田大泉湖火區(qū)地表溫度反演及時(shí)空變化特征

董敬宣，曾 強(qiáng)，李根生，趙龍輝

(1.新疆大學(xué)資源與環(huán)境科學(xué)學(xué)院，新疆 烏魯木齊 830046；2.新疆大學(xué)干旱生態(tài)環(huán)境研究所，新疆 烏魯木齊 830046)

地下煤層在自然條件或受人類活動(dòng)的影響，與氧氣接觸，發(fā)生從低溫氧化到劇烈燃燒的自燃稱為地下煤火，是一種特殊的自然災(zāi)害[1-2]。在我國西北地區(qū)，尤其新疆地區(qū)由于煤層較厚，賦存較淺，干旱少雨，加上不科學(xué)開采，地下煤層氧化自燃，使得煤火普遍存在[3]。煤火的發(fā)生對(duì)資源、生態(tài)環(huán)境危害十分嚴(yán)重?；饏^(qū)溫度變化是火區(qū)精確識(shí)別、監(jiān)測(cè)與早期預(yù)警的重要參數(shù)?；饏^(qū)溫度測(cè)量與監(jiān)測(cè)方法眾多，其中基于遙感數(shù)據(jù)的溫度反演是重要手段之一[4]。Landsat熱紅外系列數(shù)據(jù)是地表溫度反演重要的數(shù)據(jù)源。針對(duì)遙感數(shù)據(jù)的地表溫度反演，國內(nèi)外已有學(xué)者開展了大量研究：蔣大林等[5]基于Landsat8熱紅外數(shù)據(jù)反演并分析了滇池流域的地表溫度及其分布特征。Ellyett等[6]采用機(jī)載熱紅外光譜數(shù)據(jù)監(jiān)測(cè)了燃燒的煤火山。Yu 等[7]對(duì)不同算法的反演精度進(jìn)行了定量對(duì)比分析，結(jié)果表明輻射傳輸方程法精度最高；李如仁等[8]根據(jù)Landsat8衛(wèi)星影像數(shù)據(jù)采用輻射傳輸方程、單窗算法和劈窗算法等方法反演了烏達(dá)礦區(qū)的地表溫度；Stefan等[9]通過反演煤田火區(qū)地表溫度，分析地面及地下煤火時(shí)空變化；蔣衛(wèi)國等[10]通過地表溫度反演與實(shí)測(cè)數(shù)據(jù)相結(jié)合得出夜間熱紅外光譜能夠有效反演火區(qū)地表溫度；從季節(jié)角度看，冬季熱紅外光譜是反演火區(qū)地表溫度及提取火區(qū)地表溫度異常區(qū)的最佳時(shí)相，能有效區(qū)分火區(qū)與背景區(qū)。張春森等[11]基于ETM+遙感影像數(shù)據(jù)，采用普適性單通道法反演煤田火區(qū)地表溫度。

本文基于輻射傳輸方程法嘗試通過Landsat8衛(wèi)星影像數(shù)據(jù)，定量分析了2013～2016年連續(xù)4年11月份大泉湖火區(qū)地表溫度時(shí)空變化特征[12-13]，以期為該火區(qū)監(jiān)測(cè)、預(yù)警提供一定指導(dǎo)。

1 研究區(qū)概況

大泉湖火區(qū)位于烏魯木齊西8 km處，中心地理坐標(biāo)為東經(jīng)87°24′00″，北緯43°47′07″，屬低山丘陵區(qū)，火區(qū)所屬區(qū)域年平均風(fēng)速2.8 m/s，最大風(fēng)速達(dá)26.7 m/s?；饏^(qū)燃燒煤層為B7、B8，厚度分別為7 m、5 m，火區(qū)燃燒深度40～134 m。每年燃燒損失煤炭資源儲(chǔ)量約21.7萬t，受威脅的煤炭儲(chǔ)量為1 128萬t。該火區(qū)地處烏魯木齊市城市規(guī)劃區(qū)內(nèi)，每年向大氣排放CO2664 620 t，CO 4 926 t，總烴977 t，NOx786 t，SO21 805 t，煙塵434 t。火區(qū)燃燒產(chǎn)生的有毒有害氣體對(duì)烏魯木齊及周邊大氣環(huán)境造成嚴(yán)重污染，研究區(qū)地理位置見圖1。

圖1 研究區(qū)地理位置

2 數(shù)據(jù)獲取與處理

2.1 數(shù)據(jù)來源

考慮夏季地表吸收太陽熱輻射較多，會(huì)影響火區(qū)地表溫度。冬季太陽熱輻射的影響小，故選取每年11月份的遙感影像進(jìn)行分析研究。因2015年11月份研究區(qū)地面有積雪覆蓋和云量較大，故選取影像質(zhì)量良好且時(shí)間相近的10月30日的遙感影像圖。Landsat8衛(wèi)星于2013年2月11日發(fā)射后，于當(dāng)年3月18日發(fā)回了第一幅影像。選取2013年11月16日、2014年11月19日、2015年10月30日、2016年11月1日4個(gè)時(shí)期的遙感影像數(shù)據(jù)。所選取的遙感影像質(zhì)量良好，圖像清晰，研究區(qū)上空基本無云、霧影響。

2.2 算法流程

算法流程見圖2。

圖2 算法流程圖

3 分析方法

3.1 地表溫度反演

地表溫度反演算法主要有以下三種：大氣校正法(也稱為輻射傳輸方程)、單通道算法和劈窗算法。本文采用的大氣校正法基本原理(圖(2))，首先估計(jì)大氣對(duì)地表熱輻射的影響,然后從衛(wèi)星傳感器所觀測(cè)到的熱輻射總量中把這部分大氣影響減去,從而得到地表熱輻射強(qiáng)度,再將其轉(zhuǎn)化為相應(yīng)的地表溫度。

遙感數(shù)據(jù)通常采用ENVI軟件進(jìn)行處理。ENVI軟件對(duì)遙感影像數(shù)據(jù)自動(dòng)按照波長分為五個(gè)數(shù)據(jù)集，即多光譜數(shù)據(jù)(1～7波段)，全色波段數(shù)據(jù)(8波段)，卷云波段數(shù)據(jù)(9波段)，熱紅外數(shù)據(jù)(10，11波段)和質(zhì)量波段數(shù)據(jù)(12波段)，其中TIRS第11波段定標(biāo)精度的不確定性遠(yuǎn)大于第10波段，故本文對(duì)第10波段計(jì)算[14]。衛(wèi)星傳感器接收到的熱紅外輻射亮度值Lλ由三部分組成，即大氣向上輻射亮度，地面的真實(shí)輻射亮度經(jīng)過大氣層之后到達(dá)衛(wèi)星傳感器的能量和大氣向下輻射到達(dá)地面后反射的能量。Lλ的表達(dá)式可寫為式(1)(輻射傳輸方程)。

Lλ=[εB(Ts)+(1-ε)L↓]τ+L↑

(1)

式中：ε為地表比輻射率；Ts為地表真實(shí)溫度(K)；B(Ts)為黑體熱輻射亮度；L↓為大氣向上輻射亮度；L↑為大氣向下輻射亮度；τ為大氣在熱紅外波段的透過率。

溫度為T的黑體在熱紅外波段的輻射亮度B(Ts)表達(dá)為式(2)。

B(Ts)=[Lλ-L↑-τ(1-ε)L↓]/(τε)

(2)

式中，Ts可以用普朗克公式的函數(shù)獲取，見式(3)。

Ts=K2/ln(K1/B(Ts)+1)

(3)

式中，K1、K2為常數(shù)，可從影像MTL文件中獲得。地表比輻射率ε可以通過式(4)計(jì)算求得。

ε=0.004Pv+0.986

(4)

式中，Pv為植被覆蓋度，可用式(5)計(jì)算：

(NDVIVeg-NDVISoid)]

(5)

式中：NDVI為歸一化植被指數(shù)；NDVISoid為完全裸地或無植被覆蓋區(qū)域的NDVI值；NDVIVeg為完全被植被所覆蓋的像元的NDVI值。取經(jīng)驗(yàn)值NDVIVeg=0.70和NDVISoid=0.05，即當(dāng)某個(gè)像元的NDVI>0.70時(shí)，Pv=1；當(dāng)NDVI<0.05時(shí)，Pv=0。

3.2 溫度異常區(qū)提取

煤層在燃燒過程中，熱量以傳導(dǎo)方式，通過煤火上層巖石向上傳導(dǎo)，在地面和低空中形成溫度相對(duì)高于其周圍環(huán)境溫度的區(qū)域，即為溫度異常區(qū)[9]。溫度異常區(qū)的提取是探測(cè)煤火和圈定煤火范圍的前提。不同時(shí)期遙感影像會(huì)因遙感成像時(shí)間、氣候季節(jié)、地質(zhì)條件、煤火燃燒狀態(tài)等多種因素影響反演結(jié)果。因此有必要根據(jù)影像信息特征，采取統(tǒng)一方法提取溫度閾值。高于閾值區(qū)為溫度異常區(qū)，低于閾值區(qū)為背景區(qū)。相關(guān)溫度按式(6)～(8)計(jì)算。

(6)

(7)

T閾=Ta+2Tσ

(8)

式中：Ta地表溫度平均值，單位為℃；Ti為研究區(qū)地表溫度圖像中任一像元的地表溫度值，單位為℃；N為地表溫度圖像中像元總數(shù)；Tσ地表溫度標(biāo)準(zhǔn)偏差；T閾為溫度閾值，單位為℃[10]。

4 結(jié)果分析

4.1 火區(qū)不同溫度區(qū)間面積變化

通過實(shí)地考察，在遙感影像圖上截取大泉湖火區(qū)范圍，見圖3。對(duì)研究區(qū)進(jìn)行地表溫度反演，使用Arcgis數(shù)據(jù)處理軟件對(duì)地表溫度數(shù)據(jù)圖處理，不同地表溫度區(qū)間面積數(shù)據(jù)見表1，其分布見圖4。其中，Ⅰ、Ⅱ、Ⅲ、Ⅳ、Ⅴ分別代表溫度為-5～0℃、0～5 ℃、5～10 ℃、10～15 ℃、15～20 ℃。結(jié)合當(dāng)?shù)丨h(huán)境氣溫?cái)?shù)據(jù)， 從表1和圖4中可知：2013～2016年Ⅰ、Ⅱ溫度區(qū)間面積有相似變化趨勢(shì)，2013～2015年逐年減少，從2015年起有增長趨勢(shì)，都呈先減少后增加變化特征；Ⅲ溫度區(qū)間面積與Ⅰ、Ⅱ溫度區(qū)間面積變化趨勢(shì)相反，2013～2015年逐年增加，從2015年起有減少趨勢(shì)，呈先增加后減少變化特征；Ⅳ、Ⅴ溫度區(qū)間面積逐年增加，表明火區(qū)燃燒規(guī)模持續(xù)增加。

圖3 研究區(qū)影像圖

表1 火區(qū)不同溫度區(qū)間面積

時(shí)間各溫度區(qū)間面積/(104m2)Ⅰ-5～0℃Ⅱ0～5℃Ⅲ5～10℃Ⅳ10～15℃Ⅴ15～20℃2013?11?16156．15357．12131．942．790．002014?11?19134．64249．03229．6834．650．002015?10?300．00127．35437．3183．340．002016?11?0129．25154．89322．83130．6810．35

圖4 火區(qū)不同溫度區(qū)間面積變化

4.2 溫度異常區(qū)變化

為了減小環(huán)境溫度誤差，對(duì)反演結(jié)果進(jìn)行了溫度校正。統(tǒng)計(jì)研究期間地表溫度數(shù)據(jù)，得到各期地表溫度的最大值(Max)、最小值(Min)、平均值(Mean)、標(biāo)準(zhǔn)偏差(Stdev)，由公式(7)計(jì)算溫度閾值，結(jié)果見表2。由表2數(shù)據(jù)繪制研究區(qū)地表溫度變化折線圖，見圖5。

由表2、圖5可知：研究區(qū)地表溫度反演結(jié)果2013～2016年最高溫度從13.41℃升高到20.06 ℃，平均溫度從4.19 ℃升高到9.01 ℃，二者都是在2013～2014年增長速率最大，2013～2014年次之、2015～2016年火區(qū)平均溫度增長速率趨于緩和。溫度閾值從10.83 ℃逐漸升高到16.68 ℃?？傮w來看，最高溫度、平均溫度和溫度閾值變化趨勢(shì)基本一致，最高值都出現(xiàn)在2016年。

4.3 地表溫度分布空間變化

根據(jù)表2統(tǒng)計(jì)結(jié)果，利用ArcGIS軟件對(duì)地表溫度數(shù)據(jù)進(jìn)行處理,結(jié)果見圖6。圖6(a)、圖6(b)、圖6(c)和圖6(d)分別為2013～2016年研究區(qū)地表溫度反演結(jié)果。

結(jié)合表2、圖5和圖6可以看出，溫度異常區(qū)呈東西走向不連續(xù)分布，主要在研究區(qū)西部、中部和東部，共三處，正是大泉湖火區(qū)。溫度異常區(qū)位置基本不變。相比于2013年，2014和2015年三處溫度異常區(qū)范圍都有所減小，但是2016年西部和東部溫度異常區(qū)面積有擴(kuò)大趨勢(shì)，中部異常區(qū)消失。產(chǎn)生此種現(xiàn)象的原因可能與煤田滅火工作進(jìn)度有關(guān)。該火區(qū)于2015年10月開始治理，2016年開始全面實(shí)施剝離、鉆探、注水、注漿工作，使得地表以下土壤露出，增加地下煤炭的供氧量，使得地表溫度升高。

表2 地表溫度及溫度異常區(qū)統(tǒng)計(jì)結(jié)果

圖5 地表溫度變化

圖6 火區(qū)地表溫度異常區(qū)變化

4.4 溫度異常區(qū)面動(dòng)態(tài)變化值

單位時(shí)間各個(gè)等級(jí)地表溫度面積的面動(dòng)態(tài)變化值，其表達(dá)式見式(9)[15]。

(9)

式中：Sa、Sb分別表示研究期初和研究期末的各等級(jí)地表溫度面積，單位為m2；T代表兩期監(jiān)測(cè)數(shù)據(jù)相隔年數(shù)，單位為a。

根據(jù)式(9)和表2計(jì)算統(tǒng)計(jì)出溫度異常區(qū)面動(dòng)態(tài)變化值，見表3，其面積變化見圖7。由表3、圖7可知，研究區(qū)溫度異常區(qū)面積總體呈U型變化，2013年溫度異常區(qū)面積最大，2014～2016年溫度異常區(qū)面積都有所減小。2013～2014年火區(qū)溫度異常區(qū)面積從14.22×104m2減小到5.31×104m2，呈現(xiàn)大幅下降，2014～2015年面積也呈現(xiàn)減小趨勢(shì)但減少幅度相對(duì)較小，2015～2016火區(qū)溫度異常區(qū)面積呈現(xiàn)增加趨勢(shì)，且增長幅度相對(duì)較大。

圖7 動(dòng)態(tài)變化

5 結(jié) 論

1) 采用衛(wèi)星遙感影像反演地表溫度可為分析火區(qū)動(dòng)態(tài)變化提供一種手段。

2) Ⅳ和Ⅴ溫度區(qū)間面積與溫度異常區(qū)面積變化不同，不能單從高溫區(qū)面積角度分析火區(qū)燃燒狀況。

3) 2013～2016年大泉湖火區(qū)地表溫度呈現(xiàn)動(dòng)態(tài)變化：2016年火區(qū)溫度最高，火區(qū)西部和東部溫度異常區(qū)面積有擴(kuò)大趨勢(shì)，中部異常區(qū)消失。其原因主要來自該火區(qū)滅火施工。中部注水可降低地表溫度，西部、東部地表剝離增加了火區(qū)供氧通道，同時(shí)剝離工程使更多的高溫區(qū)域暴露，從而增加了火區(qū)溫度異常區(qū)面積。

[1] Pone J D N,Hein K A A,Stracher G B,et al.The spontaneous combustion of coal and its by-products in the Witbank and Sasolburg coalfields of South Africa[J].International Journal of Coal Geology,2007,72(2):124-140.

[2] Wu J J,Liu X C,Jiang W G,et al.Index system of coal fire risk assessment[C]∥EOGC Proceedings of the 2nd International Conference on Earth Observation for Global Changes.2009:1406-1412.

[3] 齊德香,蔡忠勇,曹建文,等.新疆維吾爾自治區(qū)第三次煤田火區(qū)普查報(bào)告[R].新疆煤田滅火工程局,2009.

[4] Anderson M C,Allen R G,Morse A,et al.Use of Landsat thermal imagery in monitoring evapotranspiration and managing water resources[J].Remote Sensing of Environment,2012,122:50-65.

[5] 蔣大林,匡鴻海,曹曉峰,等.基于Landsat8的地表溫度反演算法研究——以滇池流域?yàn)槔齕J].遙感技術(shù)與應(yīng)用,2015,30(3):448-454.

[6] Ellyett C D,Fleming A W.Thermal infrared imagery of The Burning Mountain coal fire[J].Remote Sensing of Environment,1974,3(1):79-86.

[7] Yu X L,Guo X L,Wu Z C．Land surface temperature retrieval from Landsat 8 TIRS-comparison between radiative transfer equation- based method,split window algorithm and single channel method[J]．Remote Sensing,2014,6(10):9829-9852.

[8] 李如仁,賁忠奇,李品,等.基于Landsat-8的煤火監(jiān)測(cè)方法研究[J].煤炭學(xué)報(bào),2016,41(7):1735-1740.

[9] Voigt Stefan,Tetzlaff Anke,Zhang Jianzhong,et al．International Journal of Coal Geology,2004,59:121．

[10] 蔣衛(wèi)國,武建軍,顧磊,等.基于夜間熱紅外光譜的地下煤火監(jiān)測(cè)方法研究[J].光譜學(xué)與光譜分析,2011,31(2):357-361.

[11] 張春森,徐肖雷,陳越峰.基于ETM+數(shù)據(jù)的煤田火區(qū)溫度異常信息提取[J].國土資源遙感,2017,29(2):201-206.

[12] 夏楠,塔西甫拉提·特依拜,張飛,等.新疆準(zhǔn)東露天煤礦區(qū)地表溫度反演及時(shí)空變化特征[J].中國礦業(yè),2016,25(1):69-73,96.

[13] 姬洪亮．多尺度熱紅外遙感數(shù)據(jù)在煤田火區(qū)信息提取中的應(yīng)用[D].烏魯木齊:新疆大學(xué),2012．

[14] 徐涵秋.Landsat 8熱紅外數(shù)據(jù)定標(biāo)參數(shù)的變化及其對(duì)地表溫度反演的影響[J].遙感學(xué)報(bào),2016,20(2):229-235.

[15] 武慧智.基于灰色模型的納木錯(cuò)面積動(dòng)態(tài)變化及趨勢(shì)預(yù)測(cè)[J].大慶師范學(xué)院學(xué)報(bào),2014(6):10-12.