基于SBAS-InSAR技術的礦區沉降監測研究

摘要：本文介紹了SBAS-InSAR 技術的數據處理流程，結合ALOS PALSAR-2影像數據，分析了礦區沉降形變特征。研究結果表明：礦區沉降區域空間模式主要表現為橢圓形，且隨著時間推移沉降范圍在逐漸擴大，沉降量級也逐漸增加，研究區域最大垂直向形變量達到了-0.20m，垂直向的最大年平均形變速率可達-0.15m/y。從空間屬性上來看，最大形變峰值處位于沉降盆地幾何中心，符合礦山開采沉陷理論，也較好地反映了監測結果的可靠性。

關鍵詞：SBAS-InSAR；礦區；沉降監測

一、引言

傳統的礦區沉降監測方法有GPS與水準測量，這兩種傳統的測量方法具有毫米級的監測精度，但整體監測成本高，維護工作量大，且監測目標為單點，只能通過單點預估出大概的沉降 [1]。合成孔徑雷達干涉測量技術（Interferometric Synthetic Aperture Radar，InSAR）是用于獲取地表三維信息的主動式微波遙感新興技術[2]，它是通過獲取同一地物目標的兩幅合成孔徑雷達（Synthetic Aperture Radar，SAR）復數圖像，并將兩幅 SAR 影像共軛相乘來得到包含有地表高程信息的干涉相位圖。然后通過數據處理，去除干涉相位圖中的高程信息來獲得干涉相位圖中的形變信息[4-6]。該方法能提供區域的大范圍連續性形變監測，具有非接觸、全天時、全天候、空間連續覆蓋等優點。

本文將InSAR技術應用于礦區地表形變監測，結合礦區特點采取SBAS處理方法，并且結合實地觀測資料對其檢驗分析，證明SBAS-InSAR技術在礦區地表沉陷監測的可靠性。

二、SBAS-InSAR的基本原理

SBAS-InSAR的數據處理流程如圖1所示。①將首先構建小基線干涉對；②根據二軌差分干涉原理，對所有干涉對進行干涉、去平地效應、去地形、濾波、相位解纏等處理；③在高相干點上建立線性形變速率，高程誤差的線性模型，估算線性形變速率和高程殘差；④采用合適濾波算法對大氣相位和地形殘余相位進行去除；⑤將形變反演結果轉換到合適地理坐標系下；⑥通過LOS向和垂直向的轉換因子，獲取垂直向的形變信息。

三、研究區域概況

從地形地貌上來看，礦區地處山脈北段，山脊線呈長條狀，該山脊的延伸方向與背斜軸向基本一致。區內山脊較為平穩，兩側為迭瓦狀的順向坡，地形坡度為9°-25°，一般表現為西坡陡，東坡緩。此外，礦區內無大的地表積水體，多為流量較小的溪溝。

另外，根據多年氣象觀測資料，該區屬于亞熱帶濕潤型氣候，氣候溫和潮濕，四季分明，雨量適中。區域年平均氣溫為17.3℃，最高氣溫達42.3℃，最低氣溫為-4.7℃。年平均總降水量為1208.3mm，年均水面蒸發量1085.0mm。

四、數據處理及結果分析

（一）數據來源

本實驗數據為2018年10月18日至2020年5月28日的5景ALOS PALSAR-2影像（其中2019年9月5日的影像為SM1成像模式，其余四景為SM3成像模式）。

（二）數據處理

1.根據小基線原則生成SAR數據干涉對和時空基線連接圖

本實驗選取2018年10月18日的影像為主影像，通過設置時空基線閾值，共生成6個小基線干涉對。干涉對的具體信息見表2。

2.差分干涉處理

干涉對組合確定后，根據二軌差分干涉原理，對所有干涉對進行干涉、去平地效應、去地形、濾波、相位解纏等處理。

由于本實驗所選用的ALOS PALSAR-2掃描模式（SM1）幅寬為50KM，且研究區域地形多起伏，故在干涉時采用10：10多視，以犧牲分辨率為代價提高信噪比及處理效率。考慮到時效性，選用TanDEM-X 90m DEM數據模擬地形相位。根據計算，影像覆蓋區域最大高差約 480米。

在獲得地形模擬相位后，即可從干涉圖中扣除地形相位以獲取地表形變相位（實際上還混疊有其他誤差）。然而由于各種失相干作用導致差分干涉圖中存有不同程度的噪聲，影響之后相位解纏質量和沉降監測精度，此處采用Goldstein 進行濾波。濾波后的影像能夠有效削弱干涉相位中的絕大部分噪聲，從而為后續相位解纏提供更高質量的基礎數據。在解纏處理中，本實驗選用最小費用流解纏分別對各差分干涉圖進行解纏。

3.線性形變建模與估計

對所有構建的干涉對實現差分處理后，依據相干系數選擇時序穩定的高相干目標。其后在高相干點上建立線性形變速率，高程誤差的線性模型，然后在線性模型的基礎上構建方程組，利用矩陣的奇異值分解（SVD）法，求解方程組并估算線性形變速率和高程殘差。

4.殘余相位估算和誤差去除

通過形變速率去除后，殘余相位中主要包含了噪聲和大氣延遲湍流部分。為此，本實驗采用高通濾波和低通濾波削弱大氣延遲，從而提高時序形變精度。其中大氣低通濾波算法和大氣的空間特性相關，由于大氣在空間上具有高相干性，采用空間窗口進行濾波時，窗口越小，其濾波效果越強；大氣高通濾波算法和大氣的時間特性有關，由于大氣在時間上具有低相干性，采用時間窗口進行濾波時，窗口越大，其濾波效果越強。因此，本實驗選擇空間域的低通濾波和時間域的高通濾波算法對大氣相位和地形殘余相位進行去除。

5.地理編碼

上述過程均是在SAR坐標系下進行的，地理編碼可以將形變反演結果轉換到地理坐標系下，以便和實測結果進行對比和分析。圖2、圖3為2018年10月18日至2020年5月28日礦區的LOS向的累計形變圖以及形變速率圖。

6.垂直向形變分解

InSAR技術獲取的是雷達視線方向的形變，該方向的形變往往不能準確反映礦區的形變特征，因此，可以通過LOS向和垂直向的轉換因子（入射角的余弦值 cosθ），獲取垂直向的形變信息。為顧及InSAR成像幾何結構對沉降監測的影響，此處忽略水平移動將LOS向形變直接轉化至垂直向，形變結果如圖4、圖5所示。

（三）結果分析

圖2中白色實線標注的區域為礦區所在區域。結果清晰地顯示了20余個地表沉降區域（紅色區域），其空間模式主要表現為橢圓形，且隨著時間推移沉降范圍在逐漸擴大，沉降量級也逐漸增加，最大LOS向形變量級達到了-0.16m。可以推測，這些地表沉降很可能為地下開采所致，因為其時空演化模式與地下開采形變演化模式非常相似。

圖3為礦區及周邊地區2018年10月18日至2020年5月28日期間的LOS向年平均形變速率圖。從中可以看出，形變速率的空間分布與形變累積分布情況基本相同，LOS向最大年平均形變速率可達-0.12m/y。

圖4為垂直向形變圖。由圖可知，主要形變區域的空間分布與LOS 向形變區域分布基本相同，垂直向的最大累積形變量可達-0.20m。

圖5為礦區及周邊地區2018年10月18日至2020年5月28日期間的垂直向年平均形變速率圖。由圖可知，垂直向的最大年平均形變速率可達-0.15m/y，數值高于LOS向的形變速率。

為直觀分析礦區開采引起地表形變時序特征，本實驗在礦區開采形成的地表沉陷盆地上方繪制剖面線D-D（圖6所示），隨后提取各期沉陷量并繪制剖面圖，具體見圖7。

如圖7所示，（a）和（b）分別表示時序LOS向累積形變及時序垂直向累積形變，從時間屬性上來看，自2018年10月18日開始，剖面累積形變隨時間推移逐漸加大，在LOS 向上形變峰值可達-0.10m，轉化至垂直向后最大峰值可達-0.12m，并且對照（a）和（b）可知垂直向累積形變趨勢與LOS向累積形變趨勢一致。此外，2019年9月5日至2020年5月28日期間，剖面累積形變增加顯著（LOS向最大形變由0.06m增至0.10m）。

從空間屬性上來看，結合累積形變與距離的關系可知最大形變峰值處大致位于沉降盆地幾何中心，較為符合礦山開采沉陷理論。此外，對照時序InSAR獲取的全盆地沉降演化圖可知剖面線橫跨兩個地表移動盆地，在圖8中體現出了兩個沉降峰值，也較好地反映了監測結果的可靠性。

五、結束語

本文基于ALOS PALSAR-2數據，利用SBAS-InSAR技術，得到某礦區2018-2020年各時間段的地表形變情況。結果顯示，研究區域存在20余個地表沉降區域，其空間模式主要表現為橢圓形，且隨著時間推移沉降范圍在逐漸擴大，沉降量級也逐漸增加，垂直向的最大累積形變量可達-0.20m，垂直向的最大年平均形變速率可達-0.15m/y。監測結果符合礦山開采沉陷理論并與實際采空區位置基本吻合，證明監測結果具有一定的有效性。

作者單位：況長虹? ? 四川 成都

參? 考? 文? 獻

[1] 劉冬，劉星年，陳超.基于CR-InSAR的煤礦區地表沉陷監測研究[J].現代測繪，2015，38（4）：26-30.

[2] 楊成生.基于D-InSAR技術的煤礦沉陷監測[D].長安大學， 2008.

[3] 劉永坦.雷達成像技術[M]. 哈爾濱工業大學出版社， 2014.

[4] 王遠堅，齊麟，姜岳.基于D-InSAR技術的金礦老采空區下沉監測與位置反演[J]. 金屬礦山， 2019（10）：68-73.

[5] 范軍林，劉飛鵬，涂梨平.基于D-InSAR技術的南昌市地表形變監測研究[J].礦產與地質，2019，33（3）：546-550.

[6] 陳優良，王之玉，栗正剛，等.基于D-InSAR的礦山地表沉陷監測研究[J].化工礦物與加工，2019，48（8）：10-14.