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

楊 潔，曾 強

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

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

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

1 D-InSAR技術原理

D-InSAR技術是以雷達干涉測量理論為基礎，通過復數據提取的相位信息為信息源獲取地表三維信息和變化信息的一項技術，其基本原理如圖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為衛星到地面的斜距；r為外部DEM高程；δr為兩個成像時刻在斜距方向上的大氣延遲。

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

2 實驗數據選擇

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

表1 Sentinel-1衛星數據主要參數

3 數據處理流程

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

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

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

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

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

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

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

(a)干涉圖

(b)差分干涉圖

(c)去平地效應圖

(d)相位解纏效果圖

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

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

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

表2 干涉對參數

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

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

圖3 地面沉降分布變化

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

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

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

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

5 融合遙感綜合研究

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

圖5 火區裂隙圈定圖

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

表4 2014年沉降變化信息

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

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

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

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

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

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

6 結語

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

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