D-InSAR技術(shù)用于水西溝火區(qū)地表沉降的分析

楊 潔，曾 強(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)

煤炭開采過(guò)程中，煤易發(fā)生自燃，進(jìn)而導(dǎo)致地下煤火的發(fā)生。地下煤火的燃燒會(huì)影響煤礦及火區(qū)周邊地區(qū)的社會(huì)和諧和可持續(xù)發(fā)展，地下煤火已成為煤礦開采區(qū)的主要災(zāi)害之一[1-4]。干旱炎熱氣候條件、易自燃煤炭屬性及落后的開采方法，使新疆地區(qū)煤火不斷發(fā)生[5-6]。燃燒火區(qū)沉降會(huì)導(dǎo)致地表產(chǎn)生裂隙，土壤結(jié)構(gòu)退化，不斷加劇水分蒸發(fā)，并對(duì)植物根系造成嚴(yán)重的機(jī)械損傷，使本已受到破壞的生態(tài)環(huán)境變得更加脆弱[7-10]，煤火造成的裂隙既是煤火燃燒的供氧通道，也是煙氣逸散的通道?；饏^(qū)地表沉降與裂隙變化在一定程度上可反映火區(qū)的動(dòng)態(tài)變化，以及火區(qū)空氣進(jìn)入和煙氣逸散的特征。有學(xué)者研究表明：通過(guò)分析研究區(qū)干涉圖，可以得到地下煤火燃燒導(dǎo)致的新沉降區(qū)及地下煤火的燃燒和地表面的沉降相互影響規(guī)律。其中合成雷達(dá)差分干涉(Differometric Synthetic Aperture Radar，D-InSAR)技術(shù)[11-12]是得到干涉圖的有效方法之一。美國(guó)學(xué)者Gabriel等[13]使用D-InSAR技術(shù)監(jiān)測(cè)地表形變，獲取SEASET L波段影像數(shù)據(jù)經(jīng)差分干涉處理后得到了研究地區(qū)的形變圖，驗(yàn)證了該技術(shù)在灌溉區(qū)沉降監(jiān)測(cè)的可行性；Donald等[14]利用機(jī)載干涉合成孔徑雷達(dá)系統(tǒng)的數(shù)字地形DEM數(shù)據(jù)，對(duì)斷層崖進(jìn)行分析并獲得地形地貌；李旺[15]、吳立新[16]等分別用 D-InSAR 技術(shù)對(duì)礦區(qū)、火山地區(qū)進(jìn)行沉降監(jiān)測(cè)，得到了研究區(qū)形變圖，這些研究成果很好地指導(dǎo)了實(shí)踐。

筆者采用雙軌差分干涉測(cè)量方法對(duì)準(zhǔn)南煤田水西溝火區(qū)的Sentinel數(shù)據(jù)進(jìn)行處理，提取火區(qū)的地面沉降變化信息，嘗試確定火區(qū)空氣進(jìn)入和煙氣逸散區(qū)域，為監(jiān)測(cè)火區(qū)的動(dòng)態(tài)變化提供一種可行的途徑和方法。

1 D-InSAR技術(shù)原理

D-InSAR技術(shù)是以雷達(dá)干涉測(cè)量理論為基礎(chǔ)，通過(guò)復(fù)數(shù)據(jù)提取的相位信息為信息源獲取地表三維信息和變化信息的一項(xiàng)技術(shù)，其基本原理如圖1所示。

圖1 D-InSAR基本原理

圖1中：S1和S2分別為主、輔圖像傳感器；B為基線距離；α為基線與水平方向的傾角；θ為主圖像入射角；H為主圖像傳感器相對(duì)地面高度；R1和R2分別為主、輔圖像斜距；ΔR為主輔圖像斜距變化量；P為地面目標(biāo)點(diǎn)；h為高程；B∥和B⊥分別為基線B在斜距向平行和垂直方向上的投影分量。

雷達(dá)干涉相位φint的組成如下[17]：

φint=φtop+φdef+φf(shuō)lat+φatm+φnoi+2kπ

(1)

式中：φtop為地形起伏引起的地形相位；φdef為地表引起的形變相位；φf(shuō)lat為參考橢球面引起的參考相位，即平地相位；φatm為大氣延遲引起的延遲相位；φnoi為處理誤差及熱噪聲引起的噪聲相位；k為整周模糊度；π為周期。

各相位分量見(jiàn)式(2)～(5)：

(2)

(3)

(4)

(5)

式中：λ為波長(zhǎng)；Δh為目標(biāo)物兩次成像的形變量；Δr為沿雷達(dá)視線向的地面形變量；R為衛(wèi)星到地面的斜距；r為外部DEM高程；δr為兩個(gè)成像時(shí)刻在斜距方向上的大氣延遲。

由式(1)～(5)可知，將其余相位去除可得到地表形變信息。

2 實(shí)驗(yàn)數(shù)據(jù)選擇

研究區(qū)域位于準(zhǔn)南煤田水西溝火區(qū)，位于東經(jīng)88°55′～88°58′，北緯43°55′～43°57′之間。水西溝火區(qū)北部表現(xiàn)為20～30 m深的塌陷狀態(tài)，而火區(qū)南部地表反潮且露頭區(qū)有芒硝狀結(jié)晶；火區(qū)燃燒煤層為6#、7#煤層，每年燃煤損失量為12.96萬(wàn)t[8]。采用NASA發(fā)布的90 m SRTM數(shù)據(jù)作為外部DEM數(shù)據(jù)用以消除干涉相位圖中的地形因素影響，并獲取了14景Sentinel衛(wèi)星Level-1產(chǎn)品中單視復(fù)數(shù)圖像(Single Look Complex，SLC)數(shù)據(jù)(https://earthdata.nasa.gov/)進(jìn)行差分干涉處理。Sentinel-1(哨兵1號(hào))衛(wèi)星系統(tǒng)及數(shù)據(jù)主要參數(shù)見(jiàn)表1。

表1 Sentinel-1衛(wèi)星數(shù)據(jù)主要參數(shù)

3 數(shù)據(jù)處理流程

使用雙軌D-InSAR技術(shù)對(duì)2014年10月7日至2015年10月3日共14景單視復(fù)數(shù)圖像(SLC)數(shù)據(jù)，進(jìn)行差分干涉處理得到水西溝火區(qū)形變場(chǎng)，并提取分析了水西溝火區(qū)地表變化信息。數(shù)據(jù)處理流程如下：

1)獲取干涉圖：經(jīng)配準(zhǔn)、干涉、干涉圖濾波、精密基線估計(jì)等步驟實(shí)施后得到干涉圖。

2)模擬地形相位：將SRTM DEM數(shù)據(jù)與主圖像進(jìn)行配準(zhǔn)，轉(zhuǎn)換到雷達(dá)坐標(biāo)系下，從而模擬得到地形相位。

3)獲取形變量：選擇主、輔圖像，利用參數(shù)文件進(jìn)行估算計(jì)算基線。

4)濾波處理：選用Goldstein自適應(yīng)性濾波方法進(jìn)行濾波處理。

5)相位解纏：采用最小二乘法，將圖像劃分為若干個(gè)小正方形格網(wǎng)，對(duì)圖像中的所有區(qū)域所有像元進(jìn)行處理，對(duì)相干性過(guò)低且小于閾值的部分掩膜進(jìn)行處理，將干涉相位和干涉圖的成像幾何關(guān)系聯(lián)系起來(lái)從而獲得地表高程。

2014年10月7日至2014年10月31日經(jīng)過(guò)差分干涉處理的干涉圖、差分干涉圖、去平地效應(yīng)圖和相位解纏效果圖見(jiàn)圖2。

(a)干涉圖

(b)差分干涉圖

(c)去平地效應(yīng)圖

(d)相位解纏效果圖

2014年10月7日至2014年10月31日解壓后的干涉圖如圖2(a)所示；差分處理后得到的差分干涉圖如圖2(b)所示，圖中可看到上部有明顯的干涉條紋，在其他部分表現(xiàn)不是很明顯；經(jīng)過(guò)去平地效應(yīng)之后去除一部分噪聲得到圖2(c)，圖中條紋更加清晰可見(jiàn)，線條更加完整；相位解纏處理得到解纏效果圖如圖2(d)所示，整幅圖像均出現(xiàn)了條紋，干涉信息得到進(jìn)一步處理。

4 差分干涉處理和結(jié)果分析

差分干涉對(duì)時(shí)間間隔和垂直基線、入射角及中心經(jīng)緯度參數(shù)如表2所示。

表2 干涉對(duì)參數(shù)

由表2可知，除2015年4月18日至2015年 6月 5日干涉對(duì)時(shí)間基線為48 d外，其余均為24 d，而2014年10月31日至2014年11月24日空間基線最長(zhǎng)為213.27 m。基線長(zhǎng)度是相干的一個(gè)重要限制因素，垂直基線過(guò)長(zhǎng)可能導(dǎo)致失相干，過(guò)短效果不明顯。每一幅圖像中心經(jīng)緯度相差不大，均覆蓋研究區(qū)域。

按照上述數(shù)據(jù)處理流程共得到11組干涉對(duì)，通過(guò)ArcGIS軟件對(duì)得到的沉降量圖像進(jìn)一步分析處理，獲取該研究區(qū)域在各時(shí)間段內(nèi)的具體地面沉降分布，其變化如圖3所示。

圖3 地面沉降分布變化

由圖3可知，2014年和2015年Sentinel數(shù)據(jù)經(jīng)過(guò)差分干涉處理生成的地面形變圖中沉降范圍大體一致，說(shuō)明了雙軌法可以從整體上確定礦區(qū)的沉降分布情況。引入外部DEM數(shù)據(jù)去除相關(guān)誤差，得到具體沉降范圍如表3所示。根據(jù)表3數(shù)據(jù)進(jìn)一步繪制干涉對(duì)沉降變化折線圖，如圖4所示。

表3 2014—2015年干涉對(duì)沉降分布范圍

圖4 不同干涉對(duì)沉降變化折線圖

由表3和圖4可知，不同干涉對(duì)的時(shí)間間隔比較均勻(14 d或14×2 d)，并且沉降變化明顯。由圖4 可知，隨著影像干涉對(duì)的時(shí)間推移，研究區(qū)沉降高低值呈現(xiàn)相似變化規(guī)律。2014年10月7日至2015年4月18日期間，沉降高、低值呈波動(dòng)變化；隨后于2015年6月29日分別減少至-21.25 mm(沉降高值)、-84.19 mm(沉降低值)；2015年6月29至2015年10月3日總體呈現(xiàn)增長(zhǎng)趨勢(shì)。

5 融合遙感綜合研究

通過(guò)單窗算法和植被模型反演2014年遙感影像數(shù)據(jù)得到地表溫度和植被覆蓋圖像。利用ArcGIS軟件重分類計(jì)算得到高溫區(qū)域和不同級(jí)別植被覆蓋區(qū)域(將植被蓋度五等分，其中[0,20%)為一級(jí)植被覆蓋區(qū)域，[20%,40%)為二級(jí)植被覆蓋區(qū)域，[40%,60%)為三級(jí)植被覆蓋區(qū)域，[60%,80%)為四級(jí)植被覆蓋區(qū)域，[80%,100%]為五級(jí)植被覆蓋區(qū)域)影像圖，疊加分析高溫區(qū)域和一級(jí)植被覆蓋區(qū)域初步得到火區(qū)裂隙可能存在的位置，火區(qū)裂隙圈定如圖5所示。

圖5 火區(qū)裂隙圈定圖

結(jié)合地表溫度和植被覆蓋場(chǎng)去接受或者拒絕一個(gè)潛在的熱異常和植被異常像素，從而圈定出A、B、C、D、E 5個(gè)火區(qū)裂隙范圍。通過(guò)shp圖像提取A、B、C、D、E 5個(gè)范圍2014年的沉降變化，具體信息見(jiàn)表4。

表4 2014年沉降變化信息

由表4數(shù)據(jù)進(jìn)一步繪制所圈定的5個(gè)火區(qū)位置沉降變化折線圖，地表沉降最大值變化折線如圖6所示，沉降最小值變化折線如圖7所示。

圖6 地表沉降最大值變化折線圖

圖7 地表沉降最小值變化折線圖

由圖6可知，位置A、B、C地表沉降最大值隨著時(shí)間變化一直減??；位置D稍有增加后驟減；位置E先明顯增大后減小。

由圖7可知，位置A地表沉降最小值變化不明顯；位置B、C均隨時(shí)間變化先明顯減小后增大，位置C表現(xiàn)為沉降，無(wú)隆升量；位置D先明顯減小后略減?。晃恢肊先增大后顯著減小。

已有Sentinel數(shù)據(jù)處理后表明2014年研究區(qū)裂隙可疑位置沉降趨勢(shì)增加明顯：2014年10月 7日至2014年10月31日干涉對(duì)僅位置C出現(xiàn)沉降，其他位置地表均為隆升；2014年12月7日至2014年12月31日干涉對(duì)5個(gè)位置沉降最小值及均值均為負(fù)值，表示發(fā)生沉降。C位置沉降幅度變化最小，相對(duì)較穩(wěn)定，但其沉降結(jié)果為負(fù)值表示沉降。

6 結(jié)語(yǔ)

1)采用雙軌法差分干涉處理方法測(cè)量水西溝火區(qū)Sentinel數(shù)據(jù)，驗(yàn)證了差分干涉測(cè)量(D-InSAR)技術(shù)用于煤田火區(qū)地面沉降監(jiān)測(cè)的可行性。

2)遙感反演處理得到水西溝火區(qū)溫度異常區(qū)及一級(jí)植被覆蓋疊加圖，可確定火區(qū)裂隙可能存在A、B、C、D、E共5個(gè)區(qū)域，為遠(yuǎn)程宏觀監(jiān)測(cè)火區(qū)裂隙區(qū)變化提供了實(shí)用方法，為進(jìn)一步分析火區(qū)燃燒系統(tǒng)動(dòng)態(tài)演化提供了可行途徑。