準(zhǔn)南煤田水西溝火區(qū)地表溫度及植被覆蓋度時(shí)空變化特征

楊 潔，曾 強(qiáng)

(1.新疆大學(xué)資源與環(huán)境科學(xué)學(xué)院，新疆 烏魯木齊 830046；2.新疆大學(xué)干旱生態(tài)環(huán)境研究所，新疆 烏魯木齊 830046；3.綠洲生態(tài)教育部重點(diǎn)實(shí)驗(yàn)室，新疆 烏魯木齊 830046)

地下煤火是指在自然環(huán)境或人為因素下，地下破碎煤體與空氣接觸后因氧化聚熱引發(fā)燃燒并不斷發(fā)展形成的大面積煤田火災(zāi)[1-5]。據(jù)不完全統(tǒng)計(jì)，我國(guó)新疆、內(nèi)蒙古、山西、寧夏、陜西等產(chǎn)煤大省(區(qū))的煤田火災(zāi)達(dá)200多處，造成的直接經(jīng)濟(jì)損失數(shù)以億計(jì)，同時(shí)燃燒對(duì)地下及近地表造成的破壞也很大，主要體現(xiàn)為煤火燃燒所引起的地表大面積的斷裂、塌陷，對(duì)土地、土壤資源及地表植被造成直接影響；煤火的熱效應(yīng)會(huì)使火區(qū)及其周圍的地表溫度上升，加速土壤涵養(yǎng)的水分蒸發(fā)，導(dǎo)致地表植被生存條件破壞甚至死亡[6]。遙感數(shù)據(jù)被廣泛應(yīng)用于煤火監(jiān)測(cè)和分析：JIANG等[7]、HUO等[8]基于Landsat TM/ETM影像，反演內(nèi)蒙古烏達(dá)火區(qū)地表溫度，分析了煤火蔓延方向；MARTHA等[9]、CHATTERJEE等[10]利用Landsat遙感數(shù)據(jù)計(jì)算出Jharia火區(qū)地表溫度，并通過(guò)反射紅外數(shù)據(jù)和煤火橫向蔓延識(shí)別出了火區(qū)位置；KUENZER等[11]利用高分辨率多光譜遙感影像和全色Quick bird數(shù)據(jù)反演烏達(dá)火區(qū)2005—2012年在滅火施工期間火區(qū)地表溫度變化以及土地覆蓋度變化；康高峰等[12]通過(guò)遙感技術(shù)對(duì)新疆的奇臺(tái)進(jìn)行了地下煤火的監(jiān)測(cè)，研究和實(shí)驗(yàn)都表明利用遙感和GIS手段可以初步圈定煤田火區(qū)的范圍；陳雅君等[13]、陳浩等[14]、夏安全等[15]基于Landsat 8影像分別反演得到了福州市、新沂市、濟(jì)南市的地表溫度并分別對(duì)研究區(qū)的土地類型進(jìn)行了解譯；楊菊等[16]選取5期遙感影像數(shù)據(jù)，以像元二分模型估算提取植被覆蓋度，得到赤水丹霞植被覆蓋度時(shí)空變化規(guī)律。

目前，針對(duì)不同研究區(qū)域地表溫度或植被覆蓋因素的研究做了許多分析，算法方法相對(duì)成熟，但表征火區(qū)方面的研究還很少，而新疆地區(qū)煤火自燃災(zāi)害嚴(yán)重，掌握火區(qū)地表溫度和植被覆蓋變化情況及趨勢(shì)具有一定的現(xiàn)實(shí)意義。本文嘗試融合地表溫度和植被覆蓋情況表征火區(qū)基本特征，掌握研究區(qū)歷年監(jiān)測(cè)數(shù)據(jù)，并應(yīng)用GM(1,1)灰色預(yù)測(cè)模型進(jìn)一步描述火區(qū)基本特征的變化趨勢(shì)，以期為運(yùn)用遙感方法監(jiān)測(cè)火區(qū)，推算火區(qū)基本特征時(shí)空變化提供技術(shù)支撐。

1 研究區(qū)概況及數(shù)據(jù)來(lái)源

水西溝位于新疆吉木薩爾縣西南20 km，東經(jīng)88°55′～88°58′之間，北緯43°55′～43°57′之間，火區(qū)面積74 975 m2，燃燒年損失儲(chǔ)量21.39萬(wàn)t，隸屬于準(zhǔn)南煤田(圖1)。導(dǎo)致水西溝煤田火災(zāi)的主要原因是該地區(qū)地質(zhì)活動(dòng)較為劇烈，地層中的水平煤層經(jīng)過(guò)多次地質(zhì)運(yùn)動(dòng)，大多變?yōu)閮A斜煤層，煤層露頭受到雨水沖刷和風(fēng)蝕作用，并與空氣中氧氣發(fā)生氧化作用，積熱增溫，使煤火發(fā)生概率大大增加?；饏^(qū)燃燒煤層厚度分別為10.49 m和3.23 m，燃燒深度為40～60 m[3]。

圖1 研究區(qū)地理位置Fig.1 Geographical location of the study area

考慮季節(jié)、云和霧的影響，選取1998—2018年(2012年缺失)覆蓋研究區(qū)夏季、云量少且圖像清晰的TM影響和OLI/TIRS影像，作為反演溫度和植被蓋度的數(shù)據(jù)源。

2 方法基本原理

2.1 地表溫度反演

考慮到精度及參數(shù)以及Landsat數(shù)據(jù)熱紅外通道數(shù)目的限制，采用覃志豪等[17]改進(jìn)的單窗算法，地表真實(shí)溫度可由式(1)計(jì)算求得。

LST=[a×(1-C-D)+(b×(1-C-D)+

C+D)×Tsensor-D×Ta]/C

(1)

式中：LST為最終反演得到的真實(shí)地表溫度，℃；a、b為擬合系數(shù)；C=ε×τ；D=(1-τ)×(1+(1-τ)×τ)；Ta為大氣平均作用溫度；ε為地表比輻射率；τ為大氣透過(guò)率；Tsensor為傳感器在衛(wèi)星高度上觀測(cè)熱輻射強(qiáng)度相對(duì)應(yīng)的地表溫度，包含有大氣和地表對(duì)熱輻射傳導(dǎo)的影響，需將其演算去除才能得到真實(shí)地表溫度，計(jì)算見(jiàn)式(2)。

Lλ=Gain×DN+Offset

(2)

式中：Lλ為輻射亮度值，W/(m2·μm·sr)；Gain為遙感圖像的增益；Offset為遙感圖像的偏移；DN為遙感圖像的灰度值。經(jīng)過(guò)輻射定標(biāo)后得到輻射亮度值Lλ，由三部分組成：大氣向上輻射亮度L↑，地面的真實(shí)輻射亮度經(jīng)過(guò)大氣層之后到達(dá)衛(wèi)星傳感器的能量，大氣向下輻射到達(dá)地面后反射的能量L↓，表達(dá)式見(jiàn)式(3)。通過(guò)NASA官網(wǎng)查詢所需要參數(shù)，可通過(guò)式(4)求得溫度為T的黑體在熱紅外波段的輻射亮度B(TS)。

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

(3)

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

(4)

獲取溫度為T的黑體在熱紅外波段的輻射亮度后，根據(jù)Planck公式反函數(shù)求得地表真實(shí)溫度，見(jiàn)式(5)。

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

(5)

式中，K1、K2為常數(shù)，可通過(guò)遙感影像頭文件查詢。

2.2 溫度異常區(qū)的提取

按照統(tǒng)一的人工閾值法[18]確定煤火提取區(qū)域的溫度閾值，低于分割點(diǎn)溫度閾值的稱為背景區(qū)，高于分割點(diǎn)溫度閾值的稱為煤火異常區(qū)。溫度閾值具體由式(6)～式(8)確定。

(6)

(7)

T閾=Tm+2Tδ

(8)

式中：TN為地表溫度圖像中任一像元溫度值；N為地表溫度圖像像元總數(shù)；Tm為地表溫度平均值；Tδ為地表溫度標(biāo)準(zhǔn)偏差；T閾為分割點(diǎn)閾值。

2.3 植被覆蓋度反演

基于與歸一化植被指數(shù)NDVI之間存在的較高相關(guān)性，采取像元二分模型[19]估算植被覆蓋度，通過(guò)式(9)和式(10)確定。

(9)

(10)

式中：NDVIsoil為完全是裸土或無(wú)植被覆蓋區(qū)域的NDVI值；NDVIveg為完全被植被所覆蓋的像元的NDVI值；NIR為近紅外波段的反射值；R為紅光波段的反射值。

2.4 GM(1,1)灰色預(yù)測(cè)模型

灰色預(yù)測(cè)模型(grey models)是根據(jù)過(guò)去及現(xiàn)在已知或非知的信息，建立從過(guò)去引申到未來(lái)的GM模型，以發(fā)現(xiàn)和掌握系統(tǒng)的發(fā)展規(guī)律，對(duì)系統(tǒng)的未來(lái)做出科學(xué)的定量預(yù)測(cè)。其基本原理是對(duì)原始數(shù)據(jù)累加淡化數(shù)據(jù)序列隨機(jī)性，提高數(shù)據(jù)序列的內(nèi)在規(guī)律，建立動(dòng)態(tài)微分方程[16]。

設(shè)X(0)為非負(fù)序列，X(0)={X(0)(i)，i=1，2，…，n}為某一預(yù)測(cè)對(duì)象的非負(fù)單調(diào)原始數(shù)據(jù)列，n為序列長(zhǎng)度，X(1)為X(0)的一次累加生成一次累加序列，見(jiàn)式(11)。

X(1)={X(1)(k),k=1,2,…,n}

(11)

式中，X(1)(k)由式(12)確定。

(12)

對(duì)X(1)可建立白化形式的微分方程，見(jiàn)式(13)。

(13)

式中，e、f為待識(shí)別參數(shù)。

(14)

(15)

(16)

(17)

3 結(jié)果與分析

3.1 地表溫度時(shí)空變化

將經(jīng)過(guò)預(yù)處理的遙感數(shù)據(jù)代入單窗算法反演模型，統(tǒng)計(jì)研究期間地表溫度信息，得到各期地表溫度的最小值、最大值、變化值、平均值、標(biāo)準(zhǔn)差和溫度閾值，結(jié)果見(jiàn)表1。

由表1可知，溫度隨時(shí)間持續(xù)發(fā)生動(dòng)態(tài)變化，溫度最高值出現(xiàn)在2017年(53.17 ℃)，變化值為16.47 ℃；溫度最小值出現(xiàn)在2018年(7.08 ℃)，最大值僅為2 186 ℃，變化值為14.78 ℃。1998—2018年期間溫度最小值、最大值及其平均值分別減小18.54 ℃、22.28 ℃、3.74 ℃和19.73 ℃。研究期間溫度最小值、最大值、平均值及溫度閾值變化趨勢(shì)基本一致，如圖2所示。

由圖2可知，1998—2006年四個(gè)溫度信息均先增加后降低，峰值均出現(xiàn)在2003年區(qū)間附近；2007—2013年溫度信息呈現(xiàn)為增加后減小再增加；2014—2018年溫度信息呈現(xiàn)為減小后增大再減小?？傮w看來(lái)，火區(qū)地表溫度呈現(xiàn)整體上升趨勢(shì)，除2018年在原先上升趨勢(shì)下表現(xiàn)為驟減。

表1 地表溫度信息統(tǒng)計(jì)結(jié)果Table 1 Statistical results of surface temperature information 單位：℃

圖2 1998—2018年溫度最小值、最大值、平均值及閾值變化折線圖Fig.2 Line chart of minimum,maximum,average andthreshold changes from 1998 to 2018

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

按照上述算法及流程，利用ArcGIS軟件對(duì)地表數(shù)據(jù)進(jìn)行處理，結(jié)果見(jiàn)圖3。由圖3可知，研究區(qū)溫度異常范圍主要不連續(xù)分布在西中東三個(gè)區(qū)域(經(jīng)度等分為西中東三個(gè)區(qū)域)，隨著年份不斷增加，溫度異常區(qū)面積因煤火燃燒狀態(tài)呈波動(dòng)變化。1998年火區(qū)三個(gè)異常區(qū)面積都不大，分布較為零散但均有存在，至2018年三個(gè)異常區(qū)面積不斷擴(kuò)展，西部向東南方向擴(kuò)展，中部向西北方向，并且不斷向東西方向延伸，東部面積卻明顯減少。這說(shuō)明溫度異常區(qū)存在新生火區(qū)，危害程度增加，同時(shí)也有部分火區(qū)消亡。

3.3 植被覆蓋度

根據(jù)水西溝火區(qū)植被實(shí)地考察資料，研究區(qū)以稀疏草地、裸地等低植被覆蓋指數(shù)土地類型為主，因此將研究區(qū)植被覆蓋度等分為5個(gè)等級(jí)。其中，0%～20%為一級(jí)；20%～40%為二級(jí)；40%～60%為三級(jí)；60%～80%為四級(jí)；80%～100%為五級(jí)。一級(jí)為低植被覆蓋區(qū)，即重點(diǎn)研究區(qū)域，主要包括裸地、水域和戈壁區(qū)域。反演得到一級(jí)植被覆蓋面積和溫度異常提取面積見(jiàn)表2。由表2數(shù)據(jù)繪制研究區(qū)一級(jí)植被覆蓋面積及溫度異常區(qū)面積關(guān)系折線圖，見(jiàn)圖4。

由表2和圖4可知，1998—2018年研究區(qū)一級(jí)植被統(tǒng)計(jì)面積從6.89 km2增加至23.14 km2，溫度異常區(qū)統(tǒng)計(jì)面積從0.71 km2增加至0.83 km2，二者變化較為一致。 一級(jí)植被覆蓋面積最大值則出現(xiàn)在2017年，為24.15 km2，一級(jí)植被覆蓋面積的最小值則出現(xiàn)在2001年，為0.65 km2。 2009—2018年溫度異常區(qū)、一級(jí)植被面積較2000—2008年均有增加，截至2018年分別增加了0.12 km2和16.25 km2。

圖3 溫度異常區(qū)分布圖Fig.3 Distribution of temperature anomalies

表2 一級(jí)植被覆蓋面積及溫度異常區(qū)面積統(tǒng)計(jì)表Table 2 Statistics of first-level vegetation coverage area and temperature anomaly area

圖4 溫度異常區(qū)和一級(jí)植被覆蓋面積關(guān)系圖Fig.4 Relationship between temperature anomaly areaand first-level vegetation coverage area

3.4GM(1,1)模型預(yù)測(cè)分析

(18)

(19)

由式(18)和式(19)分別計(jì)算可得2018—2038年水西溝火區(qū)溫度異常區(qū)面積和一級(jí)植被覆蓋度面積，見(jiàn)表3。

表3 水西溝火區(qū)2018—2038年溫度異常區(qū)面積及一級(jí)植被覆蓋度面積預(yù)測(cè)表Table 3 Forecast of temperature anomaly area and first-level vegetation coverage areain Shuixigou fire area from 2018 to 2038

圖5 溫度異常區(qū)面積及一級(jí)植被覆蓋度面積預(yù)測(cè)趨勢(shì)圖Fig.5 Prediction trend of temperature anomaly areaand first-level vegetation coverage area

由表3可知，水西溝火區(qū)溫度異常區(qū)面積和一級(jí)植被覆蓋度面積預(yù)測(cè)值，如圖5所示。由圖5可知，水西溝火區(qū)2018—2038年期間，隨著溫度異常區(qū)面積的持續(xù)增加，一級(jí)植被面積覆蓋度面積減小后持續(xù)增加。

4 結(jié) 論

1) 1998—2018年，準(zhǔn)南煤田水西溝火區(qū)溫度要素不斷變化，其中，溫度最小值、最大值、平均值及溫度閾值變化趨勢(shì)相似，相互之間影響較大?；饏^(qū)東西中三個(gè)區(qū)域有增有減，火區(qū)西部向東南方向擴(kuò)展，中部向西北方向，并且不斷向東西方向延伸，東部面積有明顯減少。

2) 一級(jí)植被覆蓋度面積受溫度異常區(qū)面積的擴(kuò)大減小影響較大，總體波動(dòng)變化，呈現(xiàn)增長(zhǎng)趨勢(shì)。2009—2018年期間溫度異常區(qū)和一級(jí)植被面積較前十年增長(zhǎng)速度均有增加，截至2018年分別增加了0.12 km2和16.25 km2?；疑A(yù)測(cè)模型結(jié)果表明在預(yù)測(cè)期內(nèi)，隨著溫度異常區(qū)面積的增大，一級(jí)植被覆蓋度面積也增加。

3) 衛(wèi)星遙感影像反演地表溫度和植被覆蓋度可融合表征火區(qū)基本特征，動(dòng)態(tài)監(jiān)測(cè)準(zhǔn)南煤田火災(zāi)的變化情況，以期為相關(guān)部門治理火災(zāi)提供科學(xué)依據(jù)。