融合遙感技術(shù)在四棵樹(shù)火區(qū)煤火識(shí)別中的應(yīng)用

劉小姣,曾 強(qiáng)

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

煤火災(zāi)害是指在自然環(huán)境或人為因素下,煤炭氧化聚熱引發(fā)燃燒并不斷發(fā)展形成的大面積煤田火災(zāi)[1]。新疆地區(qū)煤炭?jī)?chǔ)量豐富,開(kāi)采規(guī)模大,加之該地區(qū)氣候干燥少雨,煤的自然發(fā)火期較短,導(dǎo)致煤自燃問(wèn)題尤為嚴(yán)重。地下煤炭的大規(guī)模燃燒不僅浪費(fèi)了大量的能源,還破壞生態(tài)環(huán)境。因此,對(duì)煤自燃狀況的精準(zhǔn)監(jiān)測(cè),對(duì)于煤火治理具有重要意義。

近30年來(lái),國(guó)內(nèi)外學(xué)者針對(duì)煤火監(jiān)測(cè)、探測(cè)的研究取得了豐富的成果。煤火探測(cè)主要根據(jù)其燃燒所導(dǎo)致的周?chē)h(huán)境及其氛圍變化(包括電場(chǎng)、熱場(chǎng)及磁場(chǎng)等)來(lái)進(jìn)行監(jiān)測(cè),形成了化探、物探、熱探、鉆探和遙感探測(cè)五大類(lèi)方法[2-6]。其中:物探、化探、鉆探由于技術(shù)限制及成本較高,對(duì)于大面積火區(qū)的監(jiān)測(cè)識(shí)別存在一定的弊端;遙感探測(cè)手段具有探測(cè)范圍大,獲取周期短和時(shí)效性高的特點(diǎn),成為了煤火探測(cè)的有效手段。趙龍輝等[7]對(duì)Landsat影像進(jìn)行處理后獲得礦區(qū)植被分布變化;邱程錦等[8]對(duì)4景Landsat TM數(shù)據(jù)采用大氣校正法進(jìn)行處理,提取了烏達(dá)礦區(qū)的溫度異常區(qū)域,且對(duì)比實(shí)際火區(qū)分布基本一致;李峰等[9]對(duì)ASTER熱紅外波段進(jìn)行TES算法反演并提取煤火燃燒范圍來(lái)評(píng)估火區(qū)治理效果;PANDEY等[10]采用Landsat TM、ETM+和Landsat 8影像TIR波段進(jìn)行處理,根據(jù)熱異常的范圍區(qū)分了地表和地下煤火區(qū)域,相對(duì)于背景溫度,溫度變化較小的像素被指定為地下煤火,而溫差較大的像素被認(rèn)為是地表煤火。由于地表沙礫吸熱及熱島效應(yīng)所導(dǎo)致無(wú)關(guān)高溫異常區(qū)域的圈定,使其出現(xiàn)煤火區(qū)域誤判,因此僅通過(guò)熱效應(yīng)對(duì)煤火進(jìn)行監(jiān)測(cè)探測(cè)存在一定的缺陷。在此基礎(chǔ)之上,有學(xué)者針對(duì)煤火長(zhǎng)時(shí)間燃燒會(huì)造成地表坍塌變形的特點(diǎn),通過(guò)地表形變對(duì)其進(jìn)行監(jiān)測(cè)。唐日斐等[11]使用“二軌法”D-InSAR技術(shù)對(duì)ALOS PALSAR數(shù)據(jù)進(jìn)行處理提取研究區(qū)的地表沉降分布數(shù)據(jù),對(duì)礦區(qū)存在的開(kāi)采沉陷問(wèn)題進(jìn)行了分析;王洪明等[12]通過(guò)時(shí)間序列InSAR技術(shù)對(duì)內(nèi)蒙古霍林河露天礦區(qū)進(jìn)行監(jiān)測(cè)分析,證明了時(shí)序InSAR方法的實(shí)用性;劉曉帥等[13]利用D-InSAR和SBAS-InSAR技術(shù)對(duì)研究區(qū)的28景Sentinel-1A 影像進(jìn)行處理,得到地面沉降信息,分別從定量和定性2個(gè)角度對(duì)數(shù)據(jù)進(jìn)行對(duì)比分析,驗(yàn)證了SBAS-InSAR技術(shù)在微小形變監(jiān)測(cè)中精度更高;馬頂?shù)萚14]基于SBAS-InSAR技術(shù)對(duì)礦區(qū)進(jìn)行動(dòng)態(tài)監(jiān)測(cè),結(jié)果表明其在一定程度上滿足實(shí)際的監(jiān)測(cè)和管理需求;SYED T H等[15]采用N-SBAS技術(shù)對(duì)研究區(qū)2017—2020年的地表形變進(jìn)行了監(jiān)測(cè),表明該方法可以最大限度地減少大氣效應(yīng)對(duì)監(jiān)測(cè)結(jié)果的影響,同時(shí)闡述了煤礦開(kāi)采及自然地質(zhì)災(zāi)害也會(huì)導(dǎo)致地表產(chǎn)生形變,影響該方法對(duì)煤火識(shí)別的精準(zhǔn)性;楊潔等[16-18]利用遙感數(shù)據(jù)對(duì)水西溝煤火燃燒進(jìn)行了時(shí)空監(jiān)測(cè),同時(shí)利用D-InSAR技術(shù)驗(yàn)證了其監(jiān)測(cè)煤火燃燒區(qū)域沉降的可行性,并基于地表溫度、植被覆蓋度與地表沉降,提出了空氣滲入通道與煙氣逸出通道的確定方法。

綜上所述,采用單一方法監(jiān)測(cè)識(shí)別煤火時(shí)存在缺陷,易受其他因素的影響。因此,融合植被覆蓋度、區(qū)域熱異常結(jié)果,以及地表形變等多類(lèi)特征信息,對(duì)新疆準(zhǔn)南煤田四棵樹(shù)火區(qū)進(jìn)行煤火監(jiān)測(cè)識(shí)別研究,以克服單一方法存在的不足,提高地下煤火識(shí)別的準(zhǔn)確性。

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

1.1 研究區(qū)概況

準(zhǔn)南煤田四棵樹(shù)火區(qū)位于新疆天山北麓低中山區(qū),地勢(shì)西南高、東北低,海拔1 050～1 600 m,地形切割劇烈,溝谷縱橫,平均年降雨量245.6 mm,平均年蒸發(fā)量1 857.7 mm。區(qū)內(nèi)出露地層有石炭系、三疊系、侏羅系和第四系,其中侏羅系為該區(qū)唯一含煤地層。燃燒火區(qū)位于硫磺溝西側(cè),呈東西條帶狀展布于煤層露頭區(qū),位于山的南半坡,地面有明火,地表有結(jié)晶硫磺、芒硝,火區(qū)的發(fā)展趨勢(shì)由東向西,溝東部煤層露頭已燃燒殆盡,溫度正常,溝底與溝西煤層露頭高差最大達(dá)120 m。明火區(qū)著火煤層為A5號(hào)煤,煤層厚度為7 m,頂板為礫巖,底板為含礫粗砂巖。古火區(qū)燃燒煤層有A3、A4、A5號(hào)煤,煤層總厚7 m,傾角60°,燃燒深度為80 m。火區(qū)每年燃燒損失的煤炭資源量為7.63萬(wàn)t。研究區(qū)地理位置如圖1所示。

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

1.2 數(shù)據(jù)來(lái)源

研究區(qū)影像數(shù)據(jù)選取Landsat系列衛(wèi)星在2006—2020年內(nèi)(無(wú)2012年數(shù)據(jù))14景遙感數(shù)據(jù),用于植被覆蓋度及地表溫度反演。為減少太陽(yáng)輻射對(duì)地表溫度反演結(jié)果的影響,選取過(guò)境時(shí)間均處于凌晨的影像數(shù)據(jù)。由于數(shù)據(jù)獲取的局限性,2006—2011年選取Landsat 5衛(wèi)星數(shù)據(jù),2013—2020年選取Landsat 8衛(wèi)星數(shù)據(jù),云覆蓋量為0～5%,軌道號(hào)為145/29。2012年間Landsat 7影像存在條帶問(wèn)題,且其他系列不含該時(shí)期影像,故不對(duì)該年度數(shù)據(jù)進(jìn)行研究。

地表形變監(jiān)測(cè)數(shù)據(jù)采用Sentinel-1A中工作模式為IW、極化方式為VV的單視復(fù)數(shù)圖像(Single Look Complex,SLC)數(shù)據(jù)。對(duì)2014—2019年共29景雷達(dá)影像進(jìn)行時(shí)序差分處理,數(shù)據(jù)處理中所使用DEM數(shù)據(jù)為ASTER GDEM 30M 分辨率數(shù)字高程數(shù)據(jù)。雷達(dá)數(shù)據(jù)成像日期見(jiàn)表1。

表1 SAR數(shù)據(jù)成像時(shí)間

2 研究方法

2.1 植被反演

植被反演是以植被在不同波段具有不同的吸收和反射光譜特征為基礎(chǔ),采用相關(guān)算法演算得到植被覆蓋度。通過(guò)像元二分模型計(jì)算研究區(qū)的植被覆蓋度,反映該地區(qū)的植被覆蓋情況。

歸一化植被指數(shù)通過(guò)遙感數(shù)據(jù)波動(dòng)中的近紅外波段與紅外波段進(jìn)行計(jì)算,其計(jì)算公式如下:

(1)

式中:INDV為歸一化植被指數(shù);RNI為近紅外波段光譜值;R為可見(jiàn)光紅光波段光譜值。

基于歸一化植被指數(shù),使用像元二分模型對(duì)植被覆蓋度進(jìn)行計(jì)算:

(2)

式中:Pv為植被覆蓋度;INDV,min為純土壤植被指數(shù);INDV,max為純植被指數(shù);INDV,min、INDV,max分別選取置信區(qū)間內(nèi)5%和95%處的數(shù)值。

獲取Pv值后,結(jié)合該地區(qū)實(shí)地植被覆蓋度的情況,使用ArcGIS軟件對(duì)各時(shí)期植被覆蓋度進(jìn)行柵格重分類(lèi),統(tǒng)計(jì)各個(gè)區(qū)間內(nèi)像元個(gè)數(shù)及面積。

2.2 溫度反演

針對(duì)Landsat 5與Landsat 8傳感器吸收波段調(diào)整所帶來(lái)的差異性,為提高反演結(jié)果的準(zhǔn)確性,選擇不同算法對(duì)二者地表溫度進(jìn)行反演。采用經(jīng)典單窗算法[19]對(duì)Landsat 5影像進(jìn)行處理。對(duì)于Landsat 8數(shù)據(jù),研究表明在山區(qū)且地勢(shì)起伏較大的地區(qū),劈窗算法可以較好地減弱或消除大氣含水量對(duì)溫度反演的影響[20],其原理為利用傳感器內(nèi)2個(gè)相鄰的熱紅外通道(11 μm附近和12 μm附近)對(duì)大氣吸收作用的不同,結(jié)合二者觀測(cè)值的各種組合,剔除大氣對(duì)溫度反演結(jié)果的影響。基于輻射方程的劈窗算法公式如下:

Ts=A0+A1T10-A2T11

(3)

式中:Ts為地表反演溫度,K;T10、T11為L(zhǎng)andsat 8傳感器中第10和11波段的亮度溫度,K;A0、A1、A2分別為系數(shù)。

2.3 地表沉降

短基線集時(shí)序分析技術(shù)是以多主影像的干涉對(duì)為基礎(chǔ),基于高相干點(diǎn)恢復(fù)研究區(qū)域的時(shí)間序列形變信息。基本原理為假設(shè)雷達(dá)衛(wèi)星在不同時(shí)刻獲取某一地區(qū)的N+1幅SAR影像,在給定的時(shí)間和空間閾值內(nèi),生成M幅差分干涉圖,影像獲取時(shí)間為(t0,t1,…,tN),選取其中時(shí)間為ta、tb(ta>tb)時(shí)刻獲取的影像生成去除地形相位的第K幅干涉圖。在此將全部獲取的差分干涉圖的線性方程表示為M個(gè)含有N個(gè)未知數(shù)的方程,其矩陣表現(xiàn)形式如下:

Aφ=δφ

(4)

式中:A為M×N的系數(shù)矩陣;φ為各點(diǎn)未知形變相位構(gòu)成的參數(shù)矩陣;δ為M個(gè)相對(duì)具體的相位值構(gòu)成的矩陣。

將式(4)改寫(xiě)即可獲取形變速率值:

Bv=δφ

(5)

(6)

式中:B為M×N的系數(shù)矩陣;v為平均速率矩陣,mm/a;vn為第n個(gè)干涉圖對(duì)應(yīng)的形變速率;tn為第n個(gè)干涉圖獲取所對(duì)應(yīng)的時(shí)間。

在此基礎(chǔ)之上,利用最小二乘法及奇異值分解法計(jì)算得到地表平均形變速率的相位值,與時(shí)間間隔相乘后累加計(jì)算即可得到地表累計(jì)線性形變量。

3 火區(qū)數(shù)據(jù)處理結(jié)果分析

3.1 植被覆蓋度分析

在四棵樹(shù)火區(qū)實(shí)地考察后,結(jié)合該區(qū)域的植被生長(zhǎng)情況,將反演所得的植被覆蓋度進(jìn)行登記劃分,其中:Pv值在[0,0.3]內(nèi)為低植被覆蓋度區(qū)域;Pv值在(0.3,0.5]內(nèi)為較低植被覆蓋度區(qū)域;Pv值在(0.5,0.7]內(nèi)為較高植被覆蓋度區(qū)域;Pv值在(0.7,1]內(nèi)為高植被覆蓋度區(qū)域。經(jīng)過(guò)ArcGIS處理后,植被的反演結(jié)果如圖2所示。

圖2 植被反演結(jié)果示意圖

由圖2分析可知,低植被覆蓋度區(qū)域主要分布于研究區(qū)的西北及正北方向,高植被覆蓋度區(qū)域分布于西部的山坳及其兩側(cè)。在時(shí)間區(qū)間內(nèi)較高及高植被覆蓋度區(qū)域面積明顯增加。對(duì)各個(gè)植被等級(jí)像元個(gè)數(shù)統(tǒng)計(jì)可知:2006—2020年,低植被覆蓋度區(qū)域占比由38.6%降低至13.9%,面積減少了1.246 0 km2;較低植被覆蓋度區(qū)域占比由29.2%降低至22.9%,面積減少了0.318 0 km2;較高植被覆蓋度區(qū)域由16.6%增長(zhǎng)至34.2%,面積增加0.887 0 km2;高植被覆蓋度區(qū)域占比由15.6%增長(zhǎng)至29.0%,面積增加了0.676 0 km2。

3.2 溫度異常分析

地下煤火燃燒出現(xiàn)高溫,產(chǎn)生的熱量以熱輻射的形式傳導(dǎo)至地表,使存在煤火燃燒的區(qū)域地表溫度高于周?chē)h(huán)境的溫度。對(duì)研究區(qū)2006—2020年內(nèi)的遙感影像進(jìn)行地表溫度反演,結(jié)果如圖3所示。

圖3 溫度反演結(jié)果示意圖

由圖3可見(jiàn),研究區(qū)中高溫區(qū)域呈東北—西南直線分布,低溫區(qū)主要出現(xiàn)在西面山坳及山坳西側(cè)。時(shí)間區(qū)間內(nèi)反演溫度最高為51.85 ℃,最低為-0.23 ℃。

結(jié)合區(qū)域降雨量對(duì)反演結(jié)果進(jìn)行準(zhǔn)確性分析。通過(guò)網(wǎng)格地理坐標(biāo)選取研究區(qū)某點(diǎn),由于研究區(qū)范圍較小,其某點(diǎn)處月降雨量可以代表研究區(qū)內(nèi)月降雨量。2006—2019年9月份降雨量(2020年未找到該數(shù)據(jù))分別為7.81、3.84、25.94、15.61、11.34、10.37、13.55、10.51、12.67、14.06、0、8.80、4.24 mm。2008年降雨量最多為25.94 mm,2017年降雨量最少為0。2008年出現(xiàn)反演最高溫度為25.40 ℃,2017年出現(xiàn)反演最高溫度為51.85 ℃。本次研究影像獲取時(shí)間內(nèi)降雨量總體波動(dòng)程度較小,反演溫度數(shù)據(jù)較為準(zhǔn)確。

高溫異常值是在反演溫度數(shù)據(jù)基礎(chǔ)上得到的。采用人工閾值法作為提取高溫閾值的方法[21],即采取地表溫度平均值與2倍標(biāo)準(zhǔn)偏差之和作為煤火區(qū)與背景區(qū)的最佳閾值分割點(diǎn)。計(jì)算得到的2006—2011、2013—2020年溫度閾值分別為28.15、25.24、18.23、21.74、33.53、33.61、38.06、44.40、34.17、32.31、47.74、43.90、34.71、29.79 ℃。

根據(jù)異常高溫閾值計(jì)算結(jié)果,在ArcGIS中,對(duì)各時(shí)相反演溫度進(jìn)行重分類(lèi)處理,提取高于溫度閾值區(qū)域,并結(jié)合低植被覆蓋度區(qū)域變化進(jìn)行分析。高溫異常區(qū)提取結(jié)果如圖4所示,低植被覆蓋度區(qū)域面積及高溫異常區(qū)面積變化曲線如圖5所示。

圖4 高溫異常區(qū)提取結(jié)果示意圖

圖5 低植被覆蓋度區(qū)面積與高溫異常區(qū)面積關(guān)系圖

由圖4可見(jiàn),高溫異常區(qū)域呈現(xiàn)不規(guī)則分布,東部的溫度異常區(qū)數(shù)量較多。由圖5分析可得,2006—2020年,低植被覆蓋度區(qū)面積由1.947 0 km2減少至0.701 0 km2,高溫異常區(qū)面積由0.061 2 km2減少至0.001 8 km2,二者總體的變化趨勢(shì)趨于一致,均呈現(xiàn)波動(dòng)減小的趨勢(shì)。

3.3 煤火煙氣通道分析

在植被覆蓋度與地表溫度反演結(jié)果的基礎(chǔ)上,將低植被覆蓋度區(qū)域與高溫區(qū)域疊加區(qū)域作為煤火燃燒的煙氣逸散通道,將高植被覆蓋度區(qū)域與低溫區(qū)域疊加區(qū)域作為新鮮空氣滲入通道,疊加結(jié)果如圖6所示。

圖6 裂隙疊加結(jié)果示意圖

由圖6可以看出,2014—2019年,低植被覆蓋度區(qū)域+高溫異常區(qū)域疊加點(diǎn),以及高植被覆蓋度區(qū)域+低溫區(qū)域疊加所展現(xiàn)出的裂隙點(diǎn)位置隨著時(shí)間推移均處于波動(dòng)變化中。高溫異常區(qū)域+低植被覆蓋度區(qū)域疊加區(qū)域主要分布在研究區(qū)中部及東部,疊加點(diǎn)位置數(shù)量在2015年最少,僅存在2個(gè)點(diǎn)位;2019年其數(shù)量最大,存在73個(gè)點(diǎn)位。低溫區(qū)域+高植被覆蓋度區(qū)域疊加區(qū)域主要分布在西部,在北部存在較少點(diǎn)位,其數(shù)量在2015年為最小值,存在165個(gè)點(diǎn)位,在2019年為270個(gè)點(diǎn)位,為其最大值。2014—2019年,低溫區(qū)域+高植被覆蓋度區(qū)域疊加區(qū)域位置變化較小,說(shuō)明其空氣滲入通道較為穩(wěn)定;高溫異常區(qū)域+低植被覆蓋度區(qū)域疊加區(qū)域出現(xiàn)了由點(diǎn)轉(zhuǎn)面的變化現(xiàn)象,說(shuō)明煙氣逸散通道存在擴(kuò)大的趨勢(shì),煤火燃燒情況加劇。

在植被覆蓋度與地表溫度結(jié)合分析的基礎(chǔ)之上,將SBAS-InSAR地面沉降監(jiān)測(cè)數(shù)據(jù)進(jìn)行融合分析,在可疑裂隙點(diǎn)位上進(jìn)行沉降信息的提取,進(jìn)一步提高裂隙確定的準(zhǔn)確性。研究區(qū)在2014—2019年內(nèi)的地表平均形變速率如圖7所示。

圖7 地表平均形變速率圖

結(jié)合地表平均形變速率沉降信息圖,剔除疊加點(diǎn)位中無(wú)效點(diǎn)位置,點(diǎn)位數(shù)量見(jiàn)表2。在實(shí)地考察基礎(chǔ)上,選取出2類(lèi)典型通道區(qū)域進(jìn)行沉降信息處理分析,見(jiàn)圖6。

表2 裂隙點(diǎn)數(shù)量

由表2可知,高溫區(qū)域+低植被覆蓋度區(qū)域疊加區(qū)域所得到的裂隙數(shù)量在2014—2017年呈現(xiàn)出上升的趨勢(shì),在2017—2019年先降低后增高,總體數(shù)量增加了8倍左右;低溫區(qū)域+高植被覆蓋度區(qū)域疊加區(qū)域所得裂隙數(shù)量變化趨勢(shì)與之相似,總體數(shù)量表現(xiàn)為增加。在2類(lèi)典型通道范圍內(nèi)選取有效點(diǎn)位,獲取其沉降信息見(jiàn)表3,折線圖如圖8所示。

表3 典型通道代表點(diǎn)位沉降信息

圖8 典型通道代表點(diǎn)位沉降折線圖

由表3和圖8可知,2類(lèi)典型通道在時(shí)間尺度內(nèi)其沉降量變化趨勢(shì)基本趨于一致。2014年12月及2017年7月出現(xiàn)2次較大的突變。2014年前后低植被覆蓋度區(qū)域面積占比由6.28%上升至12.80%后降低,高植被覆蓋度區(qū)域面積占比由22.23%降低至18.52%后增高,高溫區(qū)面積占比由20.47%增加至45.14%,并且2014年地表最高溫度達(dá)到43.88 ℃,高于2006—2013年的最高溫度;在2017年達(dá)到最高溫度49.59 ℃,低植被覆蓋度區(qū)域呈現(xiàn)增加的趨勢(shì),高溫區(qū)域表現(xiàn)為先增大再減小的趨勢(shì),均與理論較為一致。

結(jié)合以上數(shù)據(jù)分析,研究區(qū)內(nèi)煙氣逸散通道及空氣滲入通道數(shù)量均處于增加的趨勢(shì),且隨著煤火不斷燃燒,裂隙通道下沉量也在逐漸增大,裂隙不斷擴(kuò)大。

4 結(jié)論

1)準(zhǔn)南煤田四棵樹(shù)火區(qū)在2006—2020年各等級(jí)植被覆蓋度區(qū)域分布位置變化較小。低植被及較低植被覆蓋度區(qū)域分布位置基本處于西北及北方向,隨著時(shí)間推移二者面積呈現(xiàn)減小的趨勢(shì),低植被覆蓋度區(qū)域面積占比減少了24.7%,較低植被覆蓋度區(qū)域面積占比減少了6.3%;高植被覆蓋度區(qū)域分布于西部的山坳與其兩側(cè),面積占比增加了13.4%。

2)四棵樹(shù)火區(qū)地表高溫區(qū)域分布位置呈東北—西南直線分布,低溫區(qū)域主要分布在西向山坳及山坳西側(cè)。研究區(qū)溫度表現(xiàn)為波動(dòng)變化,總體呈現(xiàn)先上升后下降的趨勢(shì),峰值出現(xiàn)在2017年,溫度為49.59 ℃,最小值出現(xiàn)在2007年,為-0.23 ℃。結(jié)合低植被覆蓋度區(qū)域面積與高溫異常區(qū)面積,二者變化趨勢(shì)基本趨于一致,總體表現(xiàn)為減小的趨勢(shì)。

3)通過(guò)融合植被覆蓋度、地表反演溫度,以及SBAS-InSAR技術(shù)確定了研究區(qū)煙氣逸散通道及空氣滲入通道位置及數(shù)量,分析結(jié)果表明該區(qū)域裂隙數(shù)量呈現(xiàn)出增加的趨勢(shì),地下煤火燃燒加劇。