劉小姣,曾 強
(1.新疆大學 生態(tài)與環(huán)境學院,新疆 烏魯木齊 830017; 2.新疆大學 干旱生態(tài)環(huán)境研究所,新疆 烏魯木齊 830017;3.綠洲生態(tài)教育部重點實驗室,新疆 烏魯木齊 830017;
煤火災害是指在自然環(huán)境或人為因素下,煤炭氧化聚熱引發(fā)燃燒并不斷發(fā)展形成的大面積煤田火災[1]。新疆地區(qū)煤炭儲量豐富,開采規(guī)模大,加之該地區(qū)氣候干燥少雨,煤的自然發(fā)火期較短,導致煤自燃問題尤為嚴重。地下煤炭的大規(guī)模燃燒不僅浪費了大量的能源,還破壞生態(tài)環(huán)境。因此,對煤自燃狀況的精準監(jiān)測,對于煤火治理具有重要意義。
近30年來,國內(nèi)外學者針對煤火監(jiān)測、探測的研究取得了豐富的成果。煤火探測主要根據(jù)其燃燒所導致的周圍環(huán)境及其氛圍變化(包括電場、熱場及磁場等)來進行監(jiān)測,形成了化探、物探、熱探、鉆探和遙感探測五大類方法[2-6]。其中:物探、化探、鉆探由于技術(shù)限制及成本較高,對于大面積火區(qū)的監(jiān)測識別存在一定的弊端;遙感探測手段具有探測范圍大,獲取周期短和時效性高的特點,成為了煤火探測的有效手段。趙龍輝等[7]對Landsat影像進行處理后獲得礦區(qū)植被分布變化;邱程錦等[8]對4景Landsat TM數(shù)據(jù)采用大氣校正法進行處理,提取了烏達礦區(qū)的溫度異常區(qū)域,且對比實際火區(qū)分布基本一致;李峰等[9]對ASTER熱紅外波段進行TES算法反演并提取煤火燃燒范圍來評估火區(qū)治理效果;PANDEY等[10]采用Landsat TM、ETM+和Landsat 8影像TIR波段進行處理,根據(jù)熱異常的范圍區(qū)分了地表和地下煤火區(qū)域,相對于背景溫度,溫度變化較小的像素被指定為地下煤火,而溫差較大的像素被認為是地表煤火。由于地表沙礫吸熱及熱島效應所導致無關(guān)高溫異常區(qū)域的圈定,使其出現(xiàn)煤火區(qū)域誤判,因此僅通過熱效應對煤火進行監(jiān)測探測存在一定的缺陷。在此基礎(chǔ)之上,有學者針對煤火長時間燃燒會造成地表坍塌變形的特點,通過地表形變對其進行監(jiān)測。唐日斐等[11]使用“二軌法”D-InSAR技術(shù)對ALOS PALSAR數(shù)據(jù)進行處理提取研究區(qū)的地表沉降分布數(shù)據(jù),對礦區(qū)存在的開采沉陷問題進行了分析;王洪明等[12]通過時間序列InSAR技術(shù)對內(nèi)蒙古霍林河露天礦區(qū)進行監(jiān)測分析,證明了時序InSAR方法的實用性;劉曉帥等[13]利用D-InSAR和SBAS-InSAR技術(shù)對研究區(qū)的28景Sentinel-1A 影像進行處理,得到地面沉降信息,分別從定量和定性2個角度對數(shù)據(jù)進行對比分析,驗證了SBAS-InSAR技術(shù)在微小形變監(jiān)測中精度更高;馬頂?shù)萚14]基于SBAS-InSAR技術(shù)對礦區(qū)進行動態(tài)監(jiān)測,結(jié)果表明其在一定程度上滿足實際的監(jiān)測和管理需求;SYED T H等[15]采用N-SBAS技術(shù)對研究區(qū)2017—2020年的地表形變進行了監(jiān)測,表明該方法可以最大限度地減少大氣效應對監(jiān)測結(jié)果的影響,同時闡述了煤礦開采及自然地質(zhì)災害也會導致地表產(chǎn)生形變,影響該方法對煤火識別的精準性;楊潔等[16-18]利用遙感數(shù)據(jù)對水西溝煤火燃燒進行了時空監(jiān)測,同時利用D-InSAR技術(shù)驗證了其監(jiān)測煤火燃燒區(qū)域沉降的可行性,并基于地表溫度、植被覆蓋度與地表沉降,提出了空氣滲入通道與煙氣逸出通道的確定方法。
綜上所述,采用單一方法監(jiān)測識別煤火時存在缺陷,易受其他因素的影響。因此,融合植被覆蓋度、區(qū)域熱異常結(jié)果,以及地表形變等多類特征信息,對新疆準南煤田四棵樹火區(qū)進行煤火監(jiān)測識別研究,以克服單一方法存在的不足,提高地下煤火識別的準確性。
準南煤田四棵樹火區(qū)位于新疆天山北麓低中山區(qū),地勢西南高、東北低,海拔1 050~1 600 m,地形切割劇烈,溝谷縱橫,平均年降雨量245.6 mm,平均年蒸發(fā)量1 857.7 mm。區(qū)內(nèi)出露地層有石炭系、三疊系、侏羅系和第四系,其中侏羅系為該區(qū)唯一含煤地層。燃燒火區(qū)位于硫磺溝西側(cè),呈東西條帶狀展布于煤層露頭區(qū),位于山的南半坡,地面有明火,地表有結(jié)晶硫磺、芒硝,火區(qū)的發(fā)展趨勢由東向西,溝東部煤層露頭已燃燒殆盡,溫度正常,溝底與溝西煤層露頭高差最大達120 m。明火區(qū)著火煤層為A5號煤,煤層厚度為7 m,頂板為礫巖,底板為含礫粗砂巖。古火區(qū)燃燒煤層有A3、A4、A5號煤,煤層總厚7 m,傾角60°,燃燒深度為80 m?;饏^(qū)每年燃燒損失的煤炭資源量為7.63萬t。研究區(qū)地理位置如圖1所示。
圖1 研究區(qū)地理位置圖
研究區(qū)影像數(shù)據(jù)選取Landsat系列衛(wèi)星在2006—2020年內(nèi)(無2012年數(shù)據(jù))14景遙感數(shù)據(jù),用于植被覆蓋度及地表溫度反演。為減少太陽輻射對地表溫度反演結(jié)果的影響,選取過境時間均處于凌晨的影像數(shù)據(jù)。由于數(shù)據(jù)獲取的局限性,2006—2011年選取Landsat 5衛(wèi)星數(shù)據(jù),2013—2020年選取Landsat 8衛(wèi)星數(shù)據(jù),云覆蓋量為0~5%,軌道號為145/29。2012年間Landsat 7影像存在條帶問題,且其他系列不含該時期影像,故不對該年度數(shù)據(jù)進行研究。
地表形變監(jiān)測數(shù)據(jù)采用Sentinel-1A中工作模式為IW、極化方式為VV的單視復數(shù)圖像(Single Look Complex,SLC)數(shù)據(jù)。對2014—2019年共29景雷達影像進行時序差分處理,數(shù)據(jù)處理中所使用DEM數(shù)據(jù)為ASTER GDEM 30M 分辨率數(shù)字高程數(shù)據(jù)。雷達數(shù)據(jù)成像日期見表1。
表1 SAR數(shù)據(jù)成像時間
植被反演是以植被在不同波段具有不同的吸收和反射光譜特征為基礎(chǔ),采用相關(guān)算法演算得到植被覆蓋度。通過像元二分模型計算研究區(qū)的植被覆蓋度,反映該地區(qū)的植被覆蓋情況。
歸一化植被指數(shù)通過遙感數(shù)據(jù)波動中的近紅外波段與紅外波段進行計算,其計算公式如下:
(1)
式中:INDV為歸一化植被指數(shù);RNI為近紅外波段光譜值;R為可見光紅光波段光譜值。
基于歸一化植被指數(shù),使用像元二分模型對植被覆蓋度進行計算:
(2)
式中:Pv為植被覆蓋度;INDV,min為純土壤植被指數(shù);INDV,max為純植被指數(shù);INDV,min、INDV,max分別選取置信區(qū)間內(nèi)5%和95%處的數(shù)值。
獲取Pv值后,結(jié)合該地區(qū)實地植被覆蓋度的情況,使用ArcGIS軟件對各時期植被覆蓋度進行柵格重分類,統(tǒng)計各個區(qū)間內(nèi)像元個數(shù)及面積。
針對Landsat 5與Landsat 8傳感器吸收波段調(diào)整所帶來的差異性,為提高反演結(jié)果的準確性,選擇不同算法對二者地表溫度進行反演。采用經(jīng)典單窗算法[19]對Landsat 5影像進行處理。對于Landsat 8數(shù)據(jù),研究表明在山區(qū)且地勢起伏較大的地區(qū),劈窗算法可以較好地減弱或消除大氣含水量對溫度反演的影響[20],其原理為利用傳感器內(nèi)2個相鄰的熱紅外通道(11 μm附近和12 μm附近)對大氣吸收作用的不同,結(jié)合二者觀測值的各種組合,剔除大氣對溫度反演結(jié)果的影響。基于輻射方程的劈窗算法公式如下:
Ts=A0+A1T10-A2T11
(3)
式中:Ts為地表反演溫度,K;T10、T11為Landsat 8傳感器中第10和11波段的亮度溫度,K;A0、A1、A2分別為系數(shù)。
短基線集時序分析技術(shù)是以多主影像的干涉對為基礎(chǔ),基于高相干點恢復研究區(qū)域的時間序列形變信息。基本原理為假設(shè)雷達衛(wèi)星在不同時刻獲取某一地區(qū)的N+1幅SAR影像,在給定的時間和空間閾值內(nèi),生成M幅差分干涉圖,影像獲取時間為(t0,t1,…,tN),選取其中時間為ta、tb(ta>tb)時刻獲取的影像生成去除地形相位的第K幅干涉圖。在此將全部獲取的差分干涉圖的線性方程表示為M個含有N個未知數(shù)的方程,其矩陣表現(xiàn)形式如下:
Aφ=δφ
(4)
式中:A為M×N的系數(shù)矩陣;φ為各點未知形變相位構(gòu)成的參數(shù)矩陣;δ為M個相對具體的相位值構(gòu)成的矩陣。
將式(4)改寫即可獲取形變速率值:
Bv=δφ
(5)
(6)
式中:B為M×N的系數(shù)矩陣;v為平均速率矩陣,mm/a;vn為第n個干涉圖對應的形變速率;tn為第n個干涉圖獲取所對應的時間。
在此基礎(chǔ)之上,利用最小二乘法及奇異值分解法計算得到地表平均形變速率的相位值,與時間間隔相乘后累加計算即可得到地表累計線性形變量。
在四棵樹火區(qū)實地考察后,結(jié)合該區(qū)域的植被生長情況,將反演所得的植被覆蓋度進行登記劃分,其中: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)過ArcGIS處理后,植被的反演結(jié)果如圖2所示。
圖2 植被反演結(jié)果示意圖
由圖2分析可知,低植被覆蓋度區(qū)域主要分布于研究區(qū)的西北及正北方向,高植被覆蓋度區(qū)域分布于西部的山坳及其兩側(cè)。在時間區(qū)間內(nèi)較高及高植被覆蓋度區(qū)域面積明顯增加。對各個植被等級像元個數(shù)統(tǒng)計可知:2006—2020年,低植被覆蓋度區(qū)域占比由38.6%降低至13.9%,面積減少了1.246 0 km2;較低植被覆蓋度區(qū)域占比由29.2%降低至22.9%,面積減少了0.318 0 km2;較高植被覆蓋度區(qū)域由16.6%增長至34.2%,面積增加0.887 0 km2;高植被覆蓋度區(qū)域占比由15.6%增長至29.0%,面積增加了0.676 0 km2。
地下煤火燃燒出現(xiàn)高溫,產(chǎn)生的熱量以熱輻射的形式傳導至地表,使存在煤火燃燒的區(qū)域地表溫度高于周圍環(huán)境的溫度。對研究區(qū)2006—2020年內(nèi)的遙感影像進行地表溫度反演,結(jié)果如圖3所示。
圖3 溫度反演結(jié)果示意圖
由圖3可見,研究區(qū)中高溫區(qū)域呈東北—西南直線分布,低溫區(qū)主要出現(xiàn)在西面山坳及山坳西側(cè)。時間區(qū)間內(nèi)反演溫度最高為51.85 ℃,最低為-0.23 ℃。
結(jié)合區(qū)域降雨量對反演結(jié)果進行準確性分析。通過網(wǎng)格地理坐標選取研究區(qū)某點,由于研究區(qū)范圍較小,其某點處月降雨量可以代表研究區(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 ℃。本次研究影像獲取時間內(nèi)降雨量總體波動程度較小,反演溫度數(shù)據(jù)較為準確。
高溫異常值是在反演溫度數(shù)據(jù)基礎(chǔ)上得到的。采用人工閾值法作為提取高溫閾值的方法[21],即采取地表溫度平均值與2倍標準偏差之和作為煤火區(qū)與背景區(qū)的最佳閾值分割點。計算得到的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ù)異常高溫閾值計算結(jié)果,在ArcGIS中,對各時相反演溫度進行重分類處理,提取高于溫度閾值區(qū)域,并結(jié)合低植被覆蓋度區(qū)域變化進行分析。高溫異常區(qū)提取結(jié)果如圖4所示,低植被覆蓋度區(qū)域面積及高溫異常區(qū)面積變化曲線如圖5所示。
圖4 高溫異常區(qū)提取結(jié)果示意圖
圖5 低植被覆蓋度區(qū)面積與高溫異常區(qū)面積關(guān)系圖
由圖4可見,高溫異常區(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,二者總體的變化趨勢趨于一致,均呈現(xiàn)波動減小的趨勢。
在植被覆蓋度與地表溫度反演結(jié)果的基礎(chǔ)上,將低植被覆蓋度區(qū)域與高溫區(qū)域疊加區(qū)域作為煤火燃燒的煙氣逸散通道,將高植被覆蓋度區(qū)域與低溫區(qū)域疊加區(qū)域作為新鮮空氣滲入通道,疊加結(jié)果如圖6所示。
圖6 裂隙疊加結(jié)果示意圖
由圖6可以看出,2014—2019年,低植被覆蓋度區(qū)域+高溫異常區(qū)域疊加點,以及高植被覆蓋度區(qū)域+低溫區(qū)域疊加所展現(xiàn)出的裂隙點位置隨著時間推移均處于波動變化中。高溫異常區(qū)域+低植被覆蓋度區(qū)域疊加區(qū)域主要分布在研究區(qū)中部及東部,疊加點位置數(shù)量在2015年最少,僅存在2個點位;2019年其數(shù)量最大,存在73個點位。低溫區(qū)域+高植被覆蓋度區(qū)域疊加區(qū)域主要分布在西部,在北部存在較少點位,其數(shù)量在2015年為最小值,存在165個點位,在2019年為270個點位,為其最大值。2014—2019年,低溫區(qū)域+高植被覆蓋度區(qū)域疊加區(qū)域位置變化較小,說明其空氣滲入通道較為穩(wěn)定;高溫異常區(qū)域+低植被覆蓋度區(qū)域疊加區(qū)域出現(xiàn)了由點轉(zhuǎn)面的變化現(xiàn)象,說明煙氣逸散通道存在擴大的趨勢,煤火燃燒情況加劇。
在植被覆蓋度與地表溫度結(jié)合分析的基礎(chǔ)之上,將SBAS-InSAR地面沉降監(jiān)測數(shù)據(jù)進行融合分析,在可疑裂隙點位上進行沉降信息的提取,進一步提高裂隙確定的準確性。研究區(qū)在2014—2019年內(nèi)的地表平均形變速率如圖7所示。
圖7 地表平均形變速率圖
結(jié)合地表平均形變速率沉降信息圖,剔除疊加點位中無效點位置,點位數(shù)量見表2。在實地考察基礎(chǔ)上,選取出2類典型通道區(qū)域進行沉降信息處理分析,見圖6。
表2 裂隙點數(shù)量
由表2可知,高溫區(qū)域+低植被覆蓋度區(qū)域疊加區(qū)域所得到的裂隙數(shù)量在2014—2017年呈現(xiàn)出上升的趨勢,在2017—2019年先降低后增高,總體數(shù)量增加了8倍左右;低溫區(qū)域+高植被覆蓋度區(qū)域疊加區(qū)域所得裂隙數(shù)量變化趨勢與之相似,總體數(shù)量表現(xiàn)為增加。在2類典型通道范圍內(nèi)選取有效點位,獲取其沉降信息見表3,折線圖如圖8所示。
表3 典型通道代表點位沉降信息
圖8 典型通道代表點位沉降折線圖
由表3和圖8可知,2類典型通道在時間尺度內(nèi)其沉降量變化趨勢基本趨于一致。2014年12月及2017年7月出現(xiàn)2次較大的突變。2014年前后低植被覆蓋度區(qū)域面積占比由6.28%上升至12.80%后降低,高植被覆蓋度區(qū)域面積占比由22.23%降低至18.52%后增高,高溫區(qū)面積占比由20.47%增加至45.14%,并且2014年地表最高溫度達到43.88 ℃,高于2006—2013年的最高溫度;在2017年達到最高溫度49.59 ℃,低植被覆蓋度區(qū)域呈現(xiàn)增加的趨勢,高溫區(qū)域表現(xiàn)為先增大再減小的趨勢,均與理論較為一致。
結(jié)合以上數(shù)據(jù)分析,研究區(qū)內(nèi)煙氣逸散通道及空氣滲入通道數(shù)量均處于增加的趨勢,且隨著煤火不斷燃燒,裂隙通道下沉量也在逐漸增大,裂隙不斷擴大。
1)準南煤田四棵樹火區(qū)在2006—2020年各等級植被覆蓋度區(qū)域分布位置變化較小。低植被及較低植被覆蓋度區(qū)域分布位置基本處于西北及北方向,隨著時間推移二者面積呈現(xiàn)減小的趨勢,低植被覆蓋度區(qū)域面積占比減少了24.7%,較低植被覆蓋度區(qū)域面積占比減少了6.3%;高植被覆蓋度區(qū)域分布于西部的山坳與其兩側(cè),面積占比增加了13.4%。
2)四棵樹火區(qū)地表高溫區(qū)域分布位置呈東北—西南直線分布,低溫區(qū)域主要分布在西向山坳及山坳西側(cè)。研究區(qū)溫度表現(xiàn)為波動變化,總體呈現(xiàn)先上升后下降的趨勢,峰值出現(xiàn)在2017年,溫度為49.59 ℃,最小值出現(xiàn)在2007年,為-0.23 ℃。結(jié)合低植被覆蓋度區(qū)域面積與高溫異常區(qū)面積,二者變化趨勢基本趨于一致,總體表現(xiàn)為減小的趨勢。
3)通過融合植被覆蓋度、地表反演溫度,以及SBAS-InSAR技術(shù)確定了研究區(qū)煙氣逸散通道及空氣滲入通道位置及數(shù)量,分析結(jié)果表明該區(qū)域裂隙數(shù)量呈現(xiàn)出增加的趨勢,地下煤火燃燒加劇。