五龍溝礦區時序InSAR地表形變監測

楊 利，薛東劍，王海方，付 林，張 婷

(1.成都理工大學地球科學學院，四川 成都 610059；2.西南交通大學地球科學與環境工程學院，四川 成都 611756)

礦產資源的開發利用是我國經濟的重要組成部分。然而礦產資源開發引發的環境地質問題日益突出，特別是礦山開采易造成崩塌、滑坡、礦震、地面沉陷等地質災害，嚴重威脅到礦區周邊人民的生產生命安全以及國家的財產安全[1]。對礦區地表進行動態監測可以有效獲取地表移動變化規律，最大限度地降低地表沉陷帶來的損失[2]。因此對礦區地表沉降進行有效監測，分析沉降特征，對開采沉陷預計與治理具有重要意義[3]。

合成孔徑雷達(synthetic aperture radar，SAR)具有全天時、全天候、多波段、多極化、高分辨率、不易受氣象條件制約、對地表具有一定的穿透力等優勢[4]。 干涉合成孔徑雷達(interferometric synthetic aperture radar，InSAR)是利用目標與兩天線的位置幾何關系，通過計算兩次過境時SAR影像的相位差來獲得數字高程模型的技術[5]。GABRIEL等[6]提出的差分干涉合成孔徑雷達(differential InSAR，D-InSAR)技術則是引入外部DEM或三軌差分、四軌差分實現對地表厘米級的微小變形的監測，高精度、高時效、低成本、大范圍是該技術的監測優勢。然而時空去相干及大氣延遲等因素使得干涉條紋無法生成或不連續，極大制約著該技術監測的精度和運用，較難完成高精度、長時間間隔地表監測[7]。為克服這些因素的影響，BERARDINO等[8]提出了小基線集方法(small baseline subset，SBAS)。 SBAS-InSAR技術通過影像的自由組合得到較多干涉，選擇空間和時間基線較小、具有高相干性的干涉圖，以減小地形對差分的影響，提高變形監測的準確性，能夠更好地消除D-InSAR技術中的時空失相關、大氣延遲相位以及地形相位誤差等影響，獲取研究區域毫米級的形變信息[9-10]。

本文針對礦區的地表沉降問題，結合地形考慮Sentinel-1A影像的可視情況，采用SBAS-InSAR技術獲取研究區的地表形變信息。利用地學統計方法，分析沉降區在空間上的分布，結合形變速率和累計形變量，對工作區范圍內的典型沉降區進行時序分析。

1 研究區概況及數據準備

1.1 研究區概況

研究區處于螢石溝-紅旗溝脆韌性剪切帶及其所形成的斷裂構造集中帶之中東段，是五龍溝礦區的一個重要構造發育區，隸屬五龍溝礦區三大主要控成礦構造區帶之一。研究區內出露地層主要以下元古界金水口群、上元古界青白口系丘吉東溝組和下古生界奧陶系祁曼塔格群變火山巖組為主，中元古界長城紀小廟組次之，溝谷和山前有大面積的第四系沖洪積物分布。 研究區中心坐標：95°56′4.14″E，36°12′24.2″N，面積約5 144 377 m2，總體構造線呈北西-南東向展布。礦石采礦均采用地下平硐開拓方式，隨著地下礦石開采量的增加，研究區水閘東溝-黃龍溝方向(北西-南東向)形成了大范圍的地下采空區，導致研究區及其附近存在崩塌、泥石流、地面塌陷等地質災害。

圖1 研究區地質構造簡圖

1.2 數據源

本文采用SRTMDEM 90M分辨率原始高程數據作為參考DEM為SBAS-InSAR提供地形信息和地理位置參考，利用從歐洲航天局質量控制中心下載的精密軌道數據對軌道信息進行修正，去除因軌道誤差引起的系統性誤差，以保證結果的準確性。由于觀測幾何的限制，不同方向形變的貢獻值不相同。結合地形考慮影像的可視性，采用的Sentinel-1A衛星以IW模式獲取的30期降軌單視復數產品進行時序反演，該數據輻射精度為1 dB，相位誤差為5°，地面分辨率為5 m×20 m(距離向×方位向)，時間跨度為2017年11月18日—2019年9月21日。

2 礦區形變監測

SBAS技術是將輸入影像自由組合，以多主影像的干涉對為基礎，基于高相干點恢復研究區域的時間序列形變信息[12]。原理是根據每個圖像中基線的時間和位置，將研究區中的SAR圖像組合成小的基線子集，進行差分干涉處理，得到較多的干涉像對，保證了較高的空間密度和時間采樣[13]。然后利用奇異值分解算法聯合多個小基線子集求解地表形變相位最小范數意義下的最小二乘解[14]。該算法有效解決了總體法方程秩虧和單個集合時間采樣率較低的問題[15]。本文利用SBAS-InSAR技術進行數據處理的流程如圖2所示。

圖2 SBAS-InSAR技術流程圖

基于青海省的地形、氣候和研究區的植被茂密程度分析，將時間基線閾值設為550 d，以避免出現無相干數據對；為避免完全空間失相關，將空間基線閾值設為45%，并進行了3D解纏，以便估算相干性低的區域、建立新的連接。利用參考DEM模擬并去除干涉圖中的地形相位，采用自適應濾波方法平滑相干斑噪聲，去除大氣延遲相位，提高干涉條紋的清晰度。選擇Delaunay MCF算法對干涉圖進行相位解纏，并剔除解纏結果不理想的像對。選取相位穩定、未發生形變的68個地面控制點對干涉結果進行軌道精煉和重去平，消除了斜坡效應，對相位偏移進行修正，估算和去除殘余的恒定相位。去除地形殘余相位、大氣延遲相位、噪聲相位是SBAS技術的核心，常分成兩步進行。第一次形變反演是建立線性模型計算出所有像對的殘余地形及地表形變速率，利用奇異值分解法解算相關參數，得出干涉像對平均形變速率和高程誤差；第二次形變反演是通過大氣濾波估算和去除大氣相位和地形殘余相位，獲取更加準確的時間序列最終位移圖。經過SBAS兩次反演得到的地表形變信息是在SAR坐標系中，將反演結果通過地理編碼轉化為地理坐標，最終獲得研究區的地表形變圖。

3 實驗結果分析與討論

在雷達圖像干涉處理的基礎上，利用GIS工具，結合研究區SBAS反演的年平均沉降速率、累計形變等信息，進而分析工作區的沉降分布、沉降量及典型區域沉降中心的形變時間序列，并對沉降區的剖面進行時序形變分析。

3.1 礦區整體形變分析

圖3 研究區形變速率圖

在監測時間范圍內，獲得的五龍溝礦區年平均速率如圖3所示，該范圍內地表移動速率相對來說比較穩定，但沉降區表現明顯。 工作區的平均沉降速率為-1.40 mm/a，最大沉降速率為-11.74 mm/a。如圖4(a)所示，-2～0 mm/a為主要的形變速率區間，占工作區面積的72.35%，而形變速率為-11.8～-10 mm/a的強烈形變區域僅占工作區面積的0.54%。

五龍溝礦區以2017年11月18日為基準，在2018-05-05期開始出現較明顯的沉降區域，隨著時間的推移，采礦活動的繼續開展及礦區自身重力的影響，沉降區域沉降量增加，沉降范圍逐漸擴大，到2019年9月21日，五龍溝礦區范圍內最大沉降量可達-21.68 mm。如圖4(b)所示，工作區范圍內的整體型變量不大，主要集中在-3～0 mm，占該工作區面積的61.57%，其中，沉降量超過-6 mm的面積約為364 344 m2，占工作區面積的7.09%。

圖4 形變面積占比圖

3.2 典型區形變分析

目前工作區有原礦堆存場、選礦廠、尾礦庫、地下采礦區(水閘東溝采區和黃龍溝采區)、礦區道路、硐口工業場地等設施。如圖5所示，工作區范圍內主要存在兩個明顯沉降區，且分布在工作區的西部，這與實際情況相符。在此基礎上圈定兩個典型沉降區域，分別為沉降區A和沉降區B。

圖5 典型沉降區分布圖

沉降區A的總面積約為124 575 m2，占工作區面積的2.42%，該沉降區形變較小，主要集中在-7～-10 mm之間，其所占面積是沉降區A面積的53.67%；沉降區B的總面積為155 311 m2，占工作區面積的3.02%，形變量明顯大于沉降區A，其形變量在-22～-19 mm范圍內的面積大約是沉降區A的7倍。工作區內兩沉降區的形變量面積占比情況見表1。

表1 典型沉降區形變量面積占比

以2017年11月18日為基準，在監測時間范圍內，地表監測點的沉降沒有出現突變現象，整體沉降位移量最大不超過30 mm。沉降區A的最大形變點位于沉降區中部偏左下方，如圖6(a)所示，其最大沉降值為-20.27 mm。該點以近似穩定的速率呈現線性沉降的趨勢。沉降區B的最大形變點位于沉降區的左上方，如圖6(b)所示，其最大沉降值達到-21.68 mm。該點的沉降呈現非線性特征變化。沉降區A主要采用分段空場法、無底柱分段崩落采礦法回采。在2017-11-18—2019-04-06時間段內，該區域空區已形成，隨著礦區地下開采的擴展，采空區的暴露面積、圍巖所承受的應力以及圍巖蠕變量相應變大，誘發頂板出現冒落，造成地表微小形變。2019-04-06期后較平穩的原因是該區通過自然冒落或強制崩落圍巖方式進行處理。沉降區B采礦方式與沉降區A相同。 在監測段內正開采3 505 m中斷，平衡拱承受的應力沒有達到破壞極限，頂板未出現冒落，因此地表呈現的沉降變化不明顯。隨著工作的推進，平衡拱受空區傾向和走向跨度擴展的影響，拱的受壓逐漸接近極限，超過極限后平衡拱被破壞，且由于礦區自身重力的影響，冒落高度增大，頂板巖層移動變形造成地表微小沉降。相較于沉降區A，沉降區B在局部地段未支護，其穩固性較差，地表沉降趨勢較為明顯，其中2019-08-04處形變起伏最為顯著。

地下空區形成后，采礦區會存在對應的沉陷中心。豎直方向上的變形將導致地面路基、土質的松弛，影響其承受能力，繼而引起和加大水平方向的傾斜和拉伸變形。沉降區A以最大形變點位為中心逐漸向四周擴散，其中向東北向擴散的趨勢較明顯；沉降區B以最大形變點位為中心呈帶狀沿西北-東南方向逐漸向四周擴散。在沉降區擴散趨勢范圍內分別繪制沉降區A和沉降區B的剖面線，通過時序分析發現，兩沉降區中心沉降(圖7)均呈現逐年增強趨勢。沉降區A呈漏斗狀逐漸向四周擴散，且向四周擴散趨勢明顯。沉降區B則是呈現盆狀逐漸向外圍擴散，擴散趨勢較為平穩。

圖6 沉降區最大形變點形變序列圖

4 結 論

1) 利用SBAS-InSAR方法得到了青海省五龍溝礦區2017年11月至2019年9月期間各時間段的地表微小緩慢形變，獲取了礦區地面沉降的時間序列圖及其空間展布。監測結果顯示工作區整體形變不大，但存在兩個明顯沉降區。

2) 形變干涉結果表明五龍溝礦區總體比較穩定。在監測時間范圍內，工作區地表整體平均沉降速率為-1.40 mm/a，其中形變速率為-2～2 mm/a的區域占實驗區面積的81.54%；地表形變量為-3～3 mm的區域占工作區面積的79.02%。工作區內最大形變速率為-11.8 mm/a，最大沉降量為-21.68 mm，其中沉降量超過-6 mm的面積為364 344.35 m2，僅占工作區面積的7.09%。

3) 沉降區A呈現線性特征，以最大沉降點為中心向四周擴散，且向東北向擴散的趨勢較明顯。沉降區B呈現非線性特征，以最大沉降點位為中心，沉降區呈帶狀沿西北-東南方向逐漸向四周擴散。