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

楊 潔，曾 強

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

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

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

1 D-InSAR技術原理

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

圖1 D-InSAR基本原理

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

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

φint=φtop+φdef+φflat+φatm+φnoi+2kπ

(1)

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

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

(2)

(3)

(4)

(5)

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

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

2 實驗數(shù)據(jù)選擇

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

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

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

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

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

2)模擬地形相位：將SRTM DEM數(shù)據(jù)與主圖像進行配準，轉換到雷達坐標系下，從而模擬得到地形相位。

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

4)濾波處理：選用Goldstein自適應性濾波方法進行濾波處理。

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

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

(a)干涉圖

(b)差分干涉圖

(c)去平地效應圖

(d)相位解纏效果圖

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

4 差分干涉處理和結果分析

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

表2 干涉對參數(shù)

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

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

圖3 地面沉降分布變化

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

表3 2014—2015年干涉對沉降分布范圍

圖4 不同干涉對沉降變化折線圖

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

5 融合遙感綜合研究

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

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

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

表4 2014年沉降變化信息

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

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

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

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

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

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

6 結語

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

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