基于多閾值目標提取的時序InSAR礦區地表沉降監測研究

劉澤洲，盧才武，章 賽，李 萌，和鄭翔

(1.西安建筑科技大學資源工程學院，陜西 西安 710055；2.西安建筑科技大學建筑學院，陜西 西安 710055)

礦區資源長期地下開采會導致地表沉降和形變，從而引發地表沉陷、裂縫以及滑坡等地質災害，不但帶來嚴重的經濟損失,而且還會危害居民正常生活及人身安全。為了有效控制礦區災害的發生，必須通過監測手段持續準確地對礦區地表變化情況進行監測，全面并及時地獲取地表沉降情況。傳統煤礦形變的監測方法雖然能精確獲取測量點的形變信息[1]，但這種僅依賴有限離散觀測數據的方法無法真實地反映地表形變的動態進程。 近年來，合成孔徑雷達差分干涉測量技術(differential interferometric synthetic aperture radar，D-InSAR)以其突出優勢發展迅猛[2-4]，尤其是多種時間序列InSAR方法的提出[5-7]，有效彌補了D-InSAR方法監測精度和可靠性受時間、空間去相干以及大氣效應等因素影響的問題[8-10]。

受自身特殊條件影響，礦區地表沉降呈現出地形條件復雜、梯度大、沉降速率大、形變范圍小等特點，使得監測難度較大[11-12]。短基線集(small baseline subsets，SBAS)技術能夠有效解決礦區地表緩慢形變監測問題[13]。在SBAS技術處理過程中，一般會通過人工選取目標點進行軌道精煉與沉降反演[14]，這種方式雖然可以有效去除殘留相位與平地效應，但若沒有提前對研究區進行考查，獲取地面目標點信息就盲目地選擇，反而會引入誤差[15]。由于礦區開采極易發生大量級地表沉陷現象，使得地表像元相干性較差，無法輕易準確獲取穩定像元點，因此目標點的選取是當下礦區InSAR監測中的重點研究問題。

針對以上問題，本文提出一種基于多閾值目標提取的SBAS礦區地表監測方法。通過設定離差閾值參數、區域窗口閾值參數與相干性閾值參數提取最為穩定的目標點，并將其應用于SBAS-InSAR處理中。以蘭州窯街礦區為研究區，探究多閾值目標提取的時序InSAR方法在礦區地表沉降監測中的應用效果。

1 方法原理

1.1 基于多閾值目標點的提取方法

地物目標的干涉相位中一般含有大量噪聲相位，通過衡量噪聲貢獻程度即可選取干涉相位較為純凈的地物目標。SAR影像是由多個像元組成，假設任意像元的高斯噪聲標準差為σn，相應的振幅?信息為Rice分布，可得到式(1)。

(1)

式中：A為某一像元；I為Bessel函數；g為地面目標的正實數反射能量。

若該分布近似接近于高斯分布，則需要g/σn>4。此時，振幅離差指數可以近似為式(2)。

(2)

式中：σ?為地物目標所對應的振幅標準差；σnI為虛部標準差；m?為地物目標所對應的時序振幅均值；D?為時序振幅離差指數。

相干系數數學表達式見式(3)。

(3)

式中：S1、S2為一組干涉對影像；M、N為局部窗口大小；i、j為提取像元的像素坐標。

本文提出的方法首先通過設定合理的離差指數閾值參數對地物目標進行首次篩選，提取干涉相位較為純凈的地物目標。其次，由于相干性能直接反應出數據處理生成的干涉相位質量的高低程度，因此設定相干性閾值參數對地物目標進行第二次篩選，消除大部分相干性差的像素點，最終獲得探測結果。

在計算相干系數時，局部移動窗口大小的設定是保證相干系數精度的重要因素之一。一般情況下，窗口越大，可靠性越高，但在一定程度上降低了分辨率，易將附近的不穩定點誤判為穩定目標點，且易遺漏離散分布單獨存在的個別目標點。因此，為了平衡圖像分辨率和相干系數估計精度這兩個互為影響的因素，在進行第二次篩選之前對研究區域進行窗口劃分。本文采用通過設定局部移動窗口閾值參數將整個區域劃分成若干個子區域，再對每個子區域的地物目標進行相干性篩選，提取大于相干性閾值的目標點為參考點，最后對所有子區域進行鑲嵌。

1.2 SBAS-InSAR方法原理

假設SAR影像在tA和tB(tA>tB)兩個時刻干涉得到第j幅干涉圖，在方位向和距離向的坐標為(x，r)，則在該點的相位可以表示為式(4)。

δφj(x)=φ(tB,x)-φ(tA,x)≈

(4)

式中：λ為雷達波長；d(tB，x)和d(tA，x)為雷達視線方向(LOS)的累計形變，d(t0,x)=0。

在(x，r)這一像元點，假設主影像IE、輔影像IS對應生成M期干涉圖可表示為式(5)。

(5)

第j(j=1，2，…，M)幅干涉圖的主影像獲取時間較輔影像晚，即IEj>ISj，則其相位組成的觀測方程見式(6)。

δφj=δφ(tIEj)-δφ(tISj)

(6)

式(6)表示含有N個未知干涉相位的M個觀測方程，其矩陣形式見式(7)。

Aφ=δφ

(7)

A是一個M×N的近似關聯矩陣，數據中生成的干涉圖可以決定矩陣A。數據都在一個小基線子集里，當M接近或者小于N時，采用奇異值分解的方法可以求得該方程最小范數意義上的最小二乘解，最終獲取時序形變量。

2 研究區概況與數據

2.1 研究區概況

窯街煤礦地處甘肅青海兩省交界處，經度范圍為102°50′～102°56′，緯度范圍為36°20′～36°26′，行政區劃屬蘭州市紅古區，該區域地處黃土高原西部，海拔1 800～2 400 m，地形較復雜[15]。

該礦區現有窯街三礦、金河煤礦、海石灣煤礦等三個開采礦區。由于多年來不斷的地下開采，礦區內地表沉陷、滑坡等地質災害頻發，嚴重影響礦區內各類建設設施以及采礦生產安全[16]，研究區信息如圖1所示。

圖1 研究區影像圖Fig.1 Image map of the study area

2.2 實驗數據準備

本研究選用具有較短時間及空間基線的哨兵一號衛星(Sentinel-1 A)干涉寬幅模式(IW)下的單視復數數據，該衛星搭載5.6 cm的C波段合成孔徑雷達，重訪周期為12 d，影像分辨率為5 m×20 m。本文選取的實驗影像時間跨度為2018年1月—2020年4月，影像雷達入射角為39.06°，方位向和距離向分辨率為13.91 m與3.70 m。此外，為提高影像配準的精度，加入精密軌道星歷數據降低軌道誤差。參考DEM數據采用90 m分辨率SRTM-DEM數據。

3 實驗數據處理

3.1 多閾值目標點提取

傳統SBAS-InSAR方法采用人工選取目標點，但在未掌握地面目標點信息的情況下進行選擇很有可能會增大誤差。因此，本文提出一種基于多閾值的目標點提取方法。對原始單視復數數據進行多視配準之后，得到研究區強度信息圖，采用本文提出的方法對所有影像數據進行分析，提取地面目標點。首先，為了最大程度提取最優點作目標點，設定離差閾值參數為3.2，進行第一次篩選，結果如圖2(a)所示，實驗共篩選出目標點29 453個。其次，由于研究區窗口較大，計算像元相干性時，在較高相干性區域易產生簇成團的控制點，在較低相干性區域易將附近的不穩定點誤判為控制點，且易遺漏離散分布單獨存在的個別控制點。因此，為了平衡高圖像分辨率與高相干系數精度，設定局部窗口閾值參數將整個區域分解成若干個子區域，再進行相干性計算，經過反復實驗，設定區域窗口閾值為16 sqkm，將整個區域劃分為132個子區間，在每個子區間內進行篩選。最后，設定相干性閾值參數為0.9，在第一次篩選結果的基礎上進行第二次篩選，篩選后的地面控制點如圖2(b)所示，實驗最終總共篩選地面控制點103個，將最終篩選結果作為后續軌道精煉的目標點用以進行地面沉降反演。

圖2 篩選的目標點Fig.2 Target point for screening

3.2 改進SBAS-InSAR處理

對研究區數據進行SBAS-InSAR處理，空間基線閾值設定臨界極限的5%，時間基線閾值設定為180 d，得到的時間和空間基線連接圖如圖3和圖4所示。

圖3 空間基線連接圖Fig.3 Diagram of spatial baseline connection

圖4 時間基線連接圖Fig.4 Diagram of time baseline connection

在干涉處理結束后，檢查所得到的所有干涉像對，去除相干性差的像對，避免形變提取時產生誤差。而后將基于多閾值方法提取的目標點代替人工選點方法，用作SBAS-InSAR軌道精煉與重去平處理。最后，利用奇異值分解法估算形變速率與殘余地形，并通過空間域低通濾波算法與時間域高通濾波算法去除大氣相位與殘余地形相位，獲取時序形變速率。

4 結果分析

4.1 沉降速率分析

經過上述處理后，得到研究區2018年1月—2020年4月期間的雷達視線方向(LOS)年平均速率圖。由式(8)可以得到研究區年平均垂直沉降速率情況，如圖5所示。

Δh=Δr/cosθ

(8)

式中：θ為雷達入射角；Δr為雷達視線方向沉降量；Δh為垂直方向沉降量。

圖5 研究區年平均垂直沉降速率圖Fig.5 Annual average vertical settlement rate map

由圖5可知，研究區內大部分地區地表沉降較為穩定，年平均沉降速率主要集中在-24～13 mm/a。A、B、C三處區域存在明顯地表沉降，具體表現為由邊緣區域向中心逐漸遞增的沉降漏斗狀，沉降速率為-25～-156 mm/a沉陷和14～87 mm/a抬升，這與文獻[16]得出的結果基本一致。結合光學影像資料對比，圖5中三處沉降漏斗與實際該煤礦開采工作面基本吻合。圖6為沿HD方向剖面側視圖，結合礦區地下開采資料進行對比分析可知，礦區地下開采是導致三個沉降區發生沉陷的主要原因。

圖6 HD方向剖面側視圖Fig.6 Side view of section in HD direction

4.2 時序沉降分析

為分析礦區地表沉降時序演變規律，對整個研究區進行時間序列分析計算后，獲得研究期內各個時間點相對于第一期影像(2018年1月2日)的沉降累計量，取其中較為重要的11個時間節點沉降累計圖(圖7)。由圖7可知，2018年5月之前研究區沒有明顯沉降，沉降量集中在-5～-15 mm之間。從2018年5月開始，研究區出現較為明顯沉降，并且隨著時間推移，沉降值和沉降范圍逐漸增大，到2020年4月沉降值達到-376 mm。其中，2018年9月、2019年1月、2020年4月沉降程度相對較大，沉降值較前分別有58 mm、54 mm和45 mm的增長。

圖7 研究區時序沉降量Fig.7 Time series subsidence in the study area

為了進一步分析開采沉陷區沉降發育過程，在研究區內選取三處(A區、B區、C區)沉降最為明顯的區域，對其做時序沉降趨勢分析，分析結果如圖8所示。與研究區煤礦地下開采資料對比分析，三處區域分別位于窯街煤礦區域中的窯街三礦、金河煤礦、海石灣煤礦三個開采礦區。三處礦區在整個研究時段內大體呈現線性下沉趨勢，其中A區和B區整體下沉速率較為穩定，C區呈階段性變化趨勢，在2018年3月之前下沉趨勢較為嚴重，沉陷值可達-78 mm，2018年4月之后下沉趨勢較為穩定。

圖8 時序沉降趨勢圖Fig.8 Time series subsidence trend

4.3 監測結果可靠性分析

為驗證多閾值目標提取的SBAS方法的有效性，使用相同的原始數據，利用傳統SBAS-InSAR方法獲取A區、B區、C區的時序地表沉降值，將兩方法監測結果進行對比分析，對比結果見圖9。

圖9 兩種方法監測結果對比Fig.9 Comparison of monitoring results of two methods

由圖9可知，A區兩種方法監測結果中差值最大為9.30 mm，最小為0.12 mm，差值的平均值為0.76 mm。對24期數據的差值離散程度進行計算，得到標準差范圍在0.06～4.66 mm之間，標準差平均值為2.06 mm。B區兩種方法監測結果中差值最大為18.80 mm，最小為0.39 mm，差值的平均值為6.80 mm。對24期數據的差值離散程度進行計算，得到標準差范圍在0～9.40 mm之間，標準差平均值為3.96 mm。C區兩種方法監測結果中差值最大為42.20 mm，最小為0.80 mm，差值的平均值為28.70 mm。對24期數據的差值離散程度進行計算，得到標準差范圍在0.40～21.10 mm之間，標準差平均值為12.30 mm。

綜上所述，基于多閾值目標提取的SBAS方法與傳統SBAS-InSAR得到的監測結果大致相吻合，基于以往SBAS-InSAR方法的監測精度，可認為本文提出方法得出的結果是可靠的。

5 結 論

由于礦區地表環境復雜，很難有效提取穩定的目標點用以SBAS處理。因此，本文提出了一種基于多閾值目標提取的時序InSAR方法，利用此方法反演出研究區的年平均沉降速率和時序累計沉降值等時間序列沉降信息，并與傳統SBAS方法監測結果進行實例對比。研究得出結論如下所述。

1) 兩種方法監測結果對比分析，三處礦區沉降差值平均值分別為0.76 mm、6.80 mm、28.70 mm，標準差平均值分別為2.06 mm、3.96 mm、12.30 mm，兩種方法監測結果較為吻合，證明基于多閾值目標提取的SBAS方法與傳統SBAS-InSAR方法在礦區地表沉降提取精度均可達到亞厘米級。

2) 傳統SBAS-InSAR方法為保證提取到更多地表相位較為穩定的目標點，在選擇研究區方面表現較為苛刻，這使得其在實際應用過程中存在許多局限性。而本文提出方法可在保證精度的同時，克服傳統SBAS-InSAR方法選取目標點在實際應用中存在的局限性，尤其在礦區發生大量級地表沉陷現象，導致地表像元較差的情況下，本文提出方法具有更好的適用性，可為后續礦區沉降監測分析提供新的監測手段。

3) 監測結果表明，研究區范圍內有三處開采沉陷盆地，分別位于三礦(A區)、金河煤礦(B區)、海石灣煤礦(C區)。最大沉陷值為-376 mm，最大年平均沉降速率為-156 mm/a，且沉降下沉趨勢在逐步增加，因此需對該礦區進行長時間序列的監測，用以提供災害預警。