基于雙軌D-InSAR技術的礦區沉降監測

朱文德，王長栓

（1.廣西壯族自治區地理信息測繪院，廣西 柳州 545006）

D-InSAR 是以合成孔徑雷達復數據提取的相位信息為信息源，利用雷達圖像的相位差信息來精確測量出圖像上每一點的三維位置和提取地面目標點微小地形變化的技術[1]。該技術不同于水準、GPS 測量等傳統的地表沉降監測方法，其具有覆蓋范圍廣、觀測連續、高效率、高分辨率、低成本、無須布設大量觀測點等優點。在地震位移測量、監測火山運動、冰川漂移、地表沉降及山體滑坡等方面展現了其獨特的優越性。1989 年，Gabriel[2]等首次利用D-InSAR 技術實現cm 級的地表形變監測。1993 年，Massonnet[3-4]等在Nature 上發表了震驚國際地震界的成果（利用D-InSAR 技術成功獲取了1992 年Landers 地區7.3 級大地震的同震形變場，與GPS 測量所獲得的結果吻合度非常高）。二十世紀九十年代以來，國內的眾多學者也在合成孔徑雷達干涉測量領域取得顯著的成果。2001 年單新建[5]等基于InSAR 技術提取了西藏瑪尼地區的DEM，所獲結果優于1∶20 萬DEM；2002 年利用D-InSAR 技術通過三通差分干涉處理獲取西藏瑪尼地震同震形變場[6]。2001 年，劉國祥[7]等利用D-InSAR 技術獲取了香港赤臘角機場近乎1 a 內的非均勻沉降場，所獲結果與離散水準監測結果相比吻合較好。2004 年，葛大慶[8]等將合成孔徑雷達干涉測量技術應用于礦山沉降監測，為煤礦區地表時空演變過程研究和開采沉陷實時動態監測提供了新的技術方法。這些研究在很大程度上推動了我國D-InSAR 技術的發展。

論文闡述了利用雙軌D-InSAR 技術進行礦區沉降監測的基本原理和技術方法，選用覆蓋山東某礦區2017-02-17、2017-03-11、2017-04-02 三景雷達SAR 影像，獲取該時段礦區的沉降信息，并針對沉降情況進行精細化分析。

1 D-InSAR技術原理及方法

1.1 D-InSAR 測量地表形變基本原理

D-InSAR 技術是在單天線模式下，進行重軌干涉測量，即在大致相同的軌道上，于2 個時刻獲取同一地區的數據。技術方法包括二軌法（Two-Pass）、三軌法（Three-Pass）和四軌法（Four-Pass）[9]。本文以二軌法為例，簡要介紹D-InSAR 測量地表形變的基本原理。

二軌法的基本思想是對同一地區變形前、變形后的兩幅SAR 圖像進行干涉處理，再引入已知的外部DEM，模擬地形相位，與干涉結果進行差分以去除干涉圖中的地形因素，獲取地表形變信息。該方法所需的SAR 數據量小，且通過引入外部DEM 消除地形相位，減少數據處理的工作量。但外部DEM 本身存在的誤差會傳遞到最終的形變結果中，此外，利用外部DEM 模擬的SAR 圖像需要與主影像進行配準，配準的精度也會影響最終形變監測精度。二軌法監測地表形變的基本原理可以表示為[10]:

式中，φ d表示研究區域發生形變前后兩幅SAR 影像生成的干涉圖；φsim,t表示采用外部DEM 模擬生成的地形相位；ΔRtow即為去除地形因素后的形變結果。二軌法數據處理流程如圖1 所示。

圖1 二軌法數據處理流程

1.2 D-InSAR 沉降監測數據處理流程

1）復圖像配準。復圖像配準的主要目的是使主、輔影像的像素對應同一地面單元。基本操作包括同名點偏移量計算、偏移量擬合及輔影像重采樣[11]。主要配準方法有:①相干系數法，選擇搜索窗口后，逐一計算每個像元上的相干系數，取其最大值對應的點作為配準點；②平均波動函數法，計算每個像元的干涉相位變化梯度，計算其平均值，值最小的點為配準點。

2）生成干涉圖及干涉質量分析。主、輔影像完成數據配準后，進行共軛相乘，生成干涉圖，此時的相位是纏繞的，數值范圍為[-π,π]。為確保最終結果的精度，需要保證干涉的質量，通常以相干系數作為評定干涉質量的指標:

式中，M×N代表目標窗口大小；u1(m,n)和u2(m,n)表示目標窗口中(m,n)位置上的2 個復型數據；u*2(m,n)表示u2(m,n)的共軛復數。相干系數的值越大，代表干涉的質量越好。

3）去平地效應及相位濾波。基線的存在致使沒有高程變化的平地產生干涉相位，為了降低解纏難度，將每一像素初始相位中斜距方向上的平地相位剔除。此外，干涉相位會受到很多噪聲的影響，噪聲的存在會降低相位解纏的精度，甚至會對最終形變結果產生影響。因此，在相位解纏前必須進行濾波處理。

4）相位解纏。差分干涉圖中的相位被纏繞在[-π,π] 范圍內，必須進行相位解纏來獲取反映地面高程的真實相位。常用的相位解纏方法有3 種:①Goldstein 枝切法；②最小費用流法；③加權解纏法。

5）形變圖生成及地理編碼。通過公式ΔRtow=λ/(4π×φreal)，將解纏后的真實相位φreal轉換成雷達視線向上的形變ΔRtow。為了更為直觀的分析研究區域的形變，通常將雷達視線向上的形變圖進行地理編碼，轉換到地圖投影坐標系下。

2 D-InSAR在礦區沉降監測中的應用實例

2.1 實驗區概況

本文選用礦產資源豐富的山東某礦區作為實驗區域，位置及范圍如圖2 所示。該礦區位于山東省西南部，含煤面積約4 826 km2，是全國重點開發的煤炭基地之一。煤炭開采會導致礦區地表出現沉降、房屋坍塌、道路裂縫等問題，為研究礦區的沉降情況，本文選取覆蓋研究區域的三景SAR 影像，基本參數如表1 所示。

圖2 研究區域示意圖

表1 干涉像對參數情況列表

2.2 數據處理關鍵步驟及成果

引入外部90 m 精度的DEM 數據，采用二軌法對上述3 個相對進行數據配準及差分處理，得到差分干涉圖（如圖3 所示）。

由圖3 可以看出干涉相對1 所得的干涉效果較清晰，沉降區域較明顯；干涉相對2、3 存在明顯的噪聲，影響干涉結果。為過濾噪聲，提取精確的相位信息，采用Goldstein 濾波對上述差分干涉圖進行濾波處理，濾波后的干涉圖如圖4 所示。

圖3 差分干涉圖

濾波后的干涉圖，去除了噪聲的影響，圖像更加清晰，標記位置的干涉條紋發生突變，條紋密集，色度跨越大，可見此區域由于煤礦長期開采，發生了地表沉降，為了獲取真實相位值，進一步做相位解纏處理。利用式ΔRtow=λ/(4π×φreal)，將解纏后的真實相位φreal轉換成雷達視線向上的形變ΔRtow，再通過地理編碼轉換為地圖投影坐標系下，圖5 為所得形變圖。

需要說明的是圖5 中的形變與圖4 表示的區域是一致的，由于所采用的SAR 數據為升軌數據，所以未經地理編碼的干涉圖的經緯度與實際地理經緯度相反，經過地理編碼后的形變圖所表示的經緯度與實際地理經緯度一致。

圖4 濾波后的干涉圖

圖5 形變結果

3 礦區沉降結果分析

由圖5 可以得知，2017-02-17～2017-04-02 期間，研究區域內出現A、B、C、D 四處明顯的沉降盆地，是由于煤礦的長期開采，導致地面出現沉降，且對周邊地區造成一定程度的影響。其中D 處的煤礦沉降尤為嚴重，為更直觀的分析礦區沉降情況，以該煤礦為例，作圖示剖面線，進行精細化分析，如圖6 所示。

圖6 剖面線沉降圖

通過對比2017-02-17 ～2017-03-11、2017-03-11～2017-04-02、2017-02-17～2017-04-02 三期成像的剖面圖可以得到:D 煤礦由于長期開采及工作面間的相互影響，形成2 個沉降漏斗。3 個干涉相對均呈現這一特征，且沉降位置及整體沉降趨勢一致。2017-02-17 ～2017-03-11 期間2 個沉降漏斗最大沉降量分別達到3.8 cm、4.7 cm，2017-03-11 ～2017-04-02 期間2 個沉降漏斗最大沉降量分別達到3.0 cm、4.8 cm，2017-02-17 ～2017-04-02 期間2 個沉降漏斗最大沉降量分別為3.0 cm、4.5 cm。干涉相對3 的沉降結果應為干涉相對1、2 之和，但由于時間基線過長，加上各種誤差的影響，該時段失相干嚴重，沉降漏斗周邊出現大量不連續的點，造成3 個成像期間沉降量有所差異，說明時間基線過長會降低雙軌D-InSAR 的監測精度。

4 結 語

本文論述了基于雙軌道D-InSAR 技術進行礦區沉降監測的基本原理和技術方法，采用2017-02-17、2017-03-11、2017-04-02 三景SAR 影像提取該時段礦區的沉降信息，并對沉降情況作進一步的精細化分析。結果表明，在2017-02-17、2017-03-11、2017-04-02成像期間，出現4 個明顯的沉降區域，其中D 處煤礦沉降最為嚴重，形成2 個連續的沉降漏斗，并呈現沉降區域擴大，沉降情況加重的趨勢。實驗證明D-InSAR 技術可以反映可靠的礦區沉降情況，能夠用于礦區沉降的實時監測，但由于外界條件制約，各種系統誤差的影響，還需通過限制時間基線，結合水準數據或GPS 等測量結果，以提高成果的精度。