基于PS-InSAR技術的老采空區地表沉陷監測與分析

盧欣奇，李學峰，張勤斌，黃海斌

(1.北京礦冶科技集團有限公司，北京 100160； 2.廣西大學資源環境與材料學院，廣西 南寧 530004)

傳統的礦區開采沉陷監測方法存在監測成本高、效率低、受天氣影響以及無法實現實時監測等問題，并且各類開采沉陷監測方法與開采沉陷預計算法難以進行有機集成[1-6]。合成孔徑雷達干涉測量(interferometric synthetic aperture radar,InSAR)是通過對相同研究區域內的兩景或多景SAR影像數據進行相干性處理，獲得研究區域內高精度的三維地形信息[7]。對其進行時間序列分析可以得到研究區域內地形微小變形信息的一種區別于傳統測量的新技術。該技術可以同時達到對研究區域的大范圍、低成本、高精度、高效率的變形監測，其理論精度可達到毫米級，這些優勢是傳統監測手段無法達到的。衛星雷達干涉測量技術在地面沉降監測方面已顯示了獨特的優勢和高精度，無疑將成為未來監測地表形變的主要技術之一[8]。

2001年FERRETTI等[9]提出PS-InSAR(permanent scatters InSAR)技術，通過永久散射體具有較高相干性，受空間、時間影響比較小的特性，為降低信號噪聲的影響提供了新的解決方式。2017年白澤朝等[10]利用Sentinel-1A雷達影像采用PS-InSAR技術對天津地區地面沉降進行了監測，綜合水文地質情況對天津地區的地表沉陷情況及原因進行了詳細的分析說明，其監測效果較為突出，為非商業SAR影像進行地表沉降監測提供了新的思路。李曼等[11]等通過PS-InSAR技術對唐山市地面沉降災害較嚴重的地區之一曹妃甸新區的地面沉降發育特征及其影響因素進行了詳細的分析，得出了該地區地面累積沉降量分布及演化狀況。地下水超采、大規模工程擾動是誘發和加劇地面沉降的外在動力的結論，為PS-InSAR研究地表沉降與實際擾動情況相結合提供了新思路。黃佳璇[12]對烏東德區域蠕動型滑坡進行檢測中使用了PS-InSAR技術與ArcGIS軟件結合的方式，對該地區的沉降變化規律進行總結，使PS-InSAR技術測量的結果直觀化、整體化，使沉降規律的總結數字化，對沉降的發展情況更具預測性。鄒昊等[13]應用PS-InSAR技術的對老采空區地表沉降監測進行了研究，將PS-InSAR技術應用到了煤礦的采空區監測當中，綜合工作面的布置情況對采空區地表變形規律進行了總結和預測。PS-InSAR技術監測大范圍內采空區地表沉降具有較大優勢，為采空區地表沉降的預防、規律總結及治理提供重要的數據基礎。

本文針對廣西平南錫基坑鉛鋅礦礦區范圍進行研究，利用2015年6月26日至2018年1月4日期間覆蓋概況區部分的21景Sentinel-1影像開展基于小數據集的PS-InSAR技術礦區監測地面沉降研究，分析其變形規律，為預測和評價老采空區殘余變形提供基礎。

1 PS-InSAR技術

PS-InSAR技術是20世紀末由意大利學者首先提出的，以解決常規干涉中大氣影響、失相干、DEM誤差等問題，極大地拓展了InSAR技術的應用前景，為精確研究地殼形變提供了強有力工具[14]。一般情況來說，PS(永久散射體)點位數量城市區域每平方千米可在數十個點以上，而郊區部分也可達每平方千米內有幾個點，這樣的資料密度，已經遠超過多數地區GPS的站位密度。

其原理為利用永久散體在長時間內保持穩定反射的特性[15]，其反射信號大大高于信號噪聲，將這些永久散體提供的信號在時間序列上的相位進行排列，在一定特征尺度范圍內大氣效應一致的假設前提下可以將InSAR結果中的大氣延遲誤差和DEM誤差從信號中消除，從而達到利用永久散射體對地表變形的精密觀測。

2 實驗區概況與數據處理

試驗區位于貴港市某礦區西北部，礦區范圍內有羅馬村、農昌村等村莊。全礦區內均為第四系所覆蓋，早在1998年進行開采活動，主要為出露地表的礦體及至地表以下30～40 m之間的鉛鋅礦體，礦區范圍內+20 m以上主要礦體已基本采完，目前地下存在大面積采空區。礦石以鉛鋅礦為主呈層狀、似層狀產于中泥盆統東崗嶺組下段白云巖中，頂板、底板均為白云巖。

實驗選取該礦1～6號采空區進行觀測，觀測期間6個工作面均已停采，該實驗是對老采空區地表沉降的變形監測。其中6號采空區面積最大，3號采空區距離地表最近，PS-InSAR技術數據處理流程見圖1。實驗數據為21景高分辨率Sentinel-1A衛星IW模式影像數據，影像選擇及時間基線見表1。采用PS-InSAR技術對數據進行分析得出超級主影像，對21幅SAR單視復數影像，經配準、輻射定標、PS探測和干涉處理，并借助分辨率為90 m、高程精度約為10 m的SRTM-3的DEM文件進行差分干涉處理，得到20幅干涉和差分干涉圖、研究區域內的5 200個PS點以及各PS點的差分干涉相位集。在考慮地表變形、高程誤差、大氣誤差、DEM誤差等其他失相干因素影響的前提下，得到每個PS點在每幅差分干涉圖上的差分干涉相位組成，其中，對形變速率增量和高程誤差增量積分，可以得到每個PS點相對于主參考點的形變速率和高程誤差。對影像進行靈活的時序分析疊加，將相鄰影像進行干涉處理生成影像干涉對以保證影像的相干性。

圖1 PS-InSAR技術數據處理流程Fig.1 PS-InSAR technology data processing flow

表1 影像的選擇及時間基線Table 1 SAR selection and SAR time baseline

序號時間時間基線序號時間時間基線12015-06-25-122017-02-143622015-08-1248132017-04-034832015-09-2948142017-05-214842015-11-1648152017-06-021252016-01-0348162017-07-083662016-05-26144172017-08-254872016-08-0672182017-09-182482016-09-2348192017-10-122492016-11-1048202017-11-2948102016-12-2848212018-01-0436112017-01-0912

最后將雷達坐標系下的PS矢量點信息導入ArcGIS進行坐標變換、重新投影坐標系轉化、克里金面生成、與采空區分布平面圖進行疊加，最終將PS點位信息可視化，得到研究區域內采空區的沉降量分布云圖、平均沉降速率云圖。

3 試驗結果與討論

3.1 試驗結果

通過對覆蓋該礦區的Sentinel-1數據進行處理取得了較好結果，在實驗區域內有5 200多個相位穩定PS點被選取出來。將時序處理得到的沉降信息結果文檔導入到ArcGIS中并繪制成整個區域的沉降信息圖，調節沉降數值的顯示分布，疊加強度影像圖作為底圖，得到各采空區沉降量分布云圖，如圖2所示，并根據累計沉降量得出研究期間平均沉降速率圖，如圖3所示。

由圖2可知，在采空區范圍內，最大沉降量主要集中在6號采空區北部，最大沉降量為205 mm；其次為4號采空區和5號采空區，最大沉降量為195 mm。與其他采空區的沉降量相比，1號采空區沉降量較小，沉降量為127～80 mm；3號采空區和6號采空區交接也出現了部分沉降，沉降量為146 mm。其他部分區域表現出地表抬升現象，抬升區域主要是3號采空區中部(該礦區主井口附近)，為該礦區的廢石堆積處。

由圖3可知，1號采空區沉降速率小，低于4.27 mm/a；3號采空區中部抬升速率高于10 mm/a屬于人為堆積造成；1號采空區、2號采空區、4號采空區、5號采空區相比6號采空區的其他區域沉降速率較小，均低于11.73 mm/a；6號采空區的北部區域主要表現出較為嚴重的沉降情況，沉降速率較大，可達到13 mm/a左右，最大處可達到15.46 mm/a。整體來看，6號采空區沉降速率較大，沉降區域范圍較大且連成一片，形成一個典型的沉降漏斗。

圖2 2015～2018年礦區沉降量分布云圖Fig.2 The settlement distribution nephogram of the mining area from 2015 to 2018

圖3 2015～2018 年礦區平均沉降速率云圖Fig.3 The average settlement rate nephogram of the mining area from 2015 to 2018

3.2 剖面線情況討論

為了詳細研究沉降漏斗區域的分布特征及其時序規律，在圖2中選取礦區內采空區沉降漏斗中心處分別繪制平行于勘探線的剖面線，分析沉降漏斗中心線兩側區域的沉降變化情況，6號采空區北部沉降量最大，1號采空區沉降量最小，3號采空區與6號采空區交接位置也有沉降產生，故對1號采空區、6號采空區及3號采空區與6號采空區交接位置進行詳細分析，其他剖面情況在以下討論中不予列出。

由圖4可知，1號采空區剖面線西側呈現三“漏斗狀”并列分布，最大沉降點東側“漏斗面”較為平滑。由采空區沉降量分布云圖可知，該沉降區域的最大沉降點A為1號采空區的東南方向，應屬于1號采空區頂板的薄弱部位，在最大沉降點西側有兩個明顯的小型漏斗區域存在，沉降量為51 mm、43 mm，為最大沉降量點A的47%、39%。

由圖5可知，3號采空區與6號采空區交接部位剖面線東側呈現兩“漏斗狀”梯度分布，西側呈“V”型階梯狀分布。由采空區沉降量分布云圖可知，該抬升區域的最大抬升點C為3號采空區與6號采空區交接部位，應屬于3號采空區與6號采空區連接部位，在最大沉降點西側有明顯的小型漏斗區域存在，沉降量為46 mm，為最大沉降量點C的32%。

圖4 Ⅰ-Ⅰ等值線剖面沉降情況Fig.4 Ⅰ-Ⅰ contour section settlement

圖5 Ⅲ-Ⅲ等值線剖面沉降情況Fig.5 Ⅲ-Ⅲ contour section settlement

由圖6可知，6號采空區剖面線西側較為平緩東側較陡且漏斗面較為平滑，整體呈現“V”形沉降梯度分布，由采空區沉降量分布云圖可知，該沉降區域的最大沉降點F為6號采空區的東南方向，應屬于6號采空區頂板的薄弱部位。在最大沉降點西側平緩部位，并未出現明顯漏斗區域，平均沉降量為87 mm，為最大沉降量點F的41%。

3.3 最大沉降點沉降速率情況討論

為了詳細研究沉最大沉降點沉降速率情況，對最大沉降點進行沉降速率分析，其各點速率變化曲線如圖7所示。

圖6 Ⅵ-Ⅵ等值線剖面沉降情況Fig.6 Ⅵ-Ⅵ contour profile settlement

圖7 最大沉降點沉降速率曲線Fig.7 Settlement rate curve of the maximum settlement point

由點A的沉降速率變化曲線可知，A點的沉降速率出現了三次峰值。該點沉降速率一直緩慢增加，達到極值后逐漸下降，持續一定時間又出現了加速狀態，沉降速率均勻，在達到第二次峰值后又出現了沉降速率減小的情況。大體上呈現“平衡狀態-非平衡狀態-新平衡狀態”這一開采沉陷規律，當采空區頂板在沉陷過程中達到平衡狀態后，又會出現新的失穩狀態，當頂板及礦柱經歷過失穩過程后又一次達到平衡時即進入到新的平衡狀態。

由點C的沉降速率變化曲線可知，C點的沉降速率出現了三次峰值。該點沉降速率最初緩慢減少，第一次達到極值后沉降速率急速上升，達到第二次極值后沉降速率又減小，持續一定時間到達第三次極值后沉降速率又開始增加。造成這種情況的原因為達到第一次峰值后受到了下部開采擾動的影響，使原本即將到達穩定狀態的地表產生了變化，產生了一段時間沉降速率較大的變形，達到第二次峰值。但這種變化并未持續較長時間，就產生了逐漸趨于穩定的狀態，直到第三次峰值的出現。在第三次峰值之后因受到其他擾動影響沉降速率先快速增加，后緩慢增加，且有即將達到峰值的趨勢。

由點F的沉降速率變化曲線可知，該點沉降速率最初為較小，后沉降速率逐漸增大，達到極值后逐漸下降，持續一定時間又出現了加速狀態，且加速度逐漸增加。

C點與F點均受到擾動影響，其沉降速率不同的原因是下部礦體開采工作面的推進導致距離這兩點的距離逐漸改變導致的。當下部開采擾動距離上部采空區較近時則該區域先出加速狀態，在擾動過后，采空區整體的穩定性下降，新的平衡狀態將會后延，當開采擾動影響過后采空區依舊會達到新的平衡。

4 數值模擬驗證

4.1 巖體力學參數及礦柱平均強度

巖體力學參數是依據巖石力學參數特性測試結果，并考慮了巖體的結構效應、地下水、節理裂隙等因素，對巖石力學參數按照POPOV的RMR巖體質量分類法進行適當的修正[16]。根據長沙礦山研究院在“金屬礦山安全技術國家重點實驗室”對礦體及圍巖采用MTS實驗系統進行實驗得出的巖石力學參數，選用飽和巖樣實驗結果，使用加拿大Rocscience公司開發的Roclab1.0軟件計算礦體及頂板、底板圍巖的巖體力學參數，從而確定基于數值模擬的礦山采空區穩定性分析所需的巖體力學參數，巖體及后期采用的充填體力學參數見表2。

考慮到礦柱尺寸效應的影響，確定了上述的巖體力學參數后，采用式(1)進行礦柱強度估算。

(1)

式中：QP為礦柱強度，MPa；Qr為巖體強度，MPa；B/H為寬高比。

各采空區的平均礦柱強度見表3。

表2 礦巖介質的力學參數Table 2 Mechanical parameters of rock medium

表3 平均礦柱強度統計表Table 3 StatisticalTable of average strength of pillar

4.2 計算分析

4.2.1 最大主應力分析

由于篇幅限制故只對以上沉降較為明顯的位置進行分析，其他采空區模擬結果不予列出。

圖8所示為各采空區頂板及礦柱縱剖面、橫剖面的最大主應力云圖。

為確保采空區經治理后達到長期穩定，礦柱安全系數采用保守方式計算，計算公式見式(2)。

F=QP/σPmax

(2)

式中：F為礦柱安全系數；QP具體數值詳見表2；σPmax為礦柱的最大垂直應力，MPa，由數值模擬給出。

礦柱的最佳強度安全系數由式(3)確定[5]。

n=1+t×σ

(3)

式中：n為礦柱最佳安全系數；t為相當于規定的可靠性系數，根據該礦山實際情況，確定t值為3.48；σ為強度安全系數的均方差，根據該礦山地質條件，計算得出σ=0.25。將數值帶入(3)式得n值約為1.9，即當F

采空區頂板、礦柱的應力狀態、不穩定礦柱個數等信息詳見表4。

表4 采空區應力狀態統計表Table 4 StatisticalTable of stress state in goaf

圖8 最大主應力云圖Fig.8 Maximum principal stress nephogram

由圖8(c)和圖8(f)可知，在礦柱上的最大主應力呈現壓應力狀態，且最大主應力出現于礦柱的中上部位，因各采空區礦柱眾多，礦柱的最大主應力云圖將按上述規律在礦柱的中上部位做橫剖面圖，各采空區橫剖面的最大主應力云圖見圖8(b)、圖8(d)和圖8(e)。

由最大主應力可知，在1號采空區的中部、3號采空區整體及6號采空區北部及南部均有大量不穩定的礦柱集中，因此這些部位為各區域的薄弱位置。

4.2.2 最小主應力分析

由圖9可知，最小主應力呈現出壓應力及拉應力狀態。1號采空區、2號采空區、4號采空區、5號采空區及6號采空區僅頂板局部小范圍內出現拉應力，拉應力范圍是0～0.20 MPa；3號采空區頂板絕大部分為拉應力狀態，拉應力范圍是0～0.27 MPa，而頂板圍巖抗拉強度僅0.19 MPa，說明3號采空區頂板有較大幾率出現拉伸破壞。

4.3 結果對比

PS-InSAR測量得出的沉降量較大區域與數值模擬得出的薄弱區域對照信息詳見表5。

由表5可知，該礦區1號采空區和6號采空區通過PS-InSAR技術測量得到的沉降量較大區域與數值模擬所得的薄弱區域基本一致；3號采空區數值模擬得出的薄弱區域為整個采空區范圍，而通過PS-InSAR技術測量得出的較大沉降量的部位為3號采空區的北部。該情況產生的原因為3號采空區南部的廢渣堆放導致了該部位地表變形出現上升。在廢渣堆放處邊緣有明顯下沉邊界，該邊界外部為3號采空區與6號采空區交接部位(即3號采空區與6號采空區交界位置)，在該區域出現了沉降量較大的區域，綜合PS-InSAR測量得出的沉降量較大區域與數值模擬得出的薄弱區域進行對比可知，通過PS-InSAR技術對采空區沉降監測的數據具有較高的真實性，可以為采空區地表沉陷防治及沉降規律總結、采空區治理提供數據支持。

圖9 最小主應力云圖Fig.9 Minimum principal stress nephogram

表5 數值模擬及實測沉降區域對照表Table 5 Numerical simulation and measured settlement area check list

采空區編號PS-InSAR測量得出的沉降量較大區域數值模擬得出的薄弱區域1號中部中部3號北部整體6號北部及南部北部及南部

5 結 論

1) 利用PS-InSAR技術對地表沉降進行檢測，繪制出地表沉降量分布云圖、等值線剖面、觀測點沉降速率變化曲線。與傳統方法相比大大提高了地表沉降研究的直觀性、整體性及預測性。

2) 通過數值模擬對PS-InSAR沉降監測結果進行驗證，對比結果基本一致，因此PS-InSAR技術獲得的采空區沉降數據可以為地表沉陷防治及沉降規律總結、采空區治理提供支持。

3) 1號采空區沒有受到下部采動影響，其沉降趨勢較為穩定，符合“平衡狀態-非平衡狀態-新平衡狀態”這一開采沉陷規律。

4) 3號采空區與6號采空區交接位置受到開采擾動影響出現沉降速率較快增加的現象，達到峰值后整體沉降趨勢較為穩定，但在數據截至時沉降速率依舊在加快。

5) 基于分析結果，1號采空區正處于新的平衡狀態后的減速下沉狀態。6號采空區在受到下部開采擾動的影響后出現了加速沉降狀態，沉降的極限情況有待研究。

6) 1號采空區整體沉降變化較為穩定，暫時不需要進行處理。但3號采空區與6號采空區交接位置及6號采空區沉降部位在數據截止時沉降速率持續增加，需對這兩個部位進行追蹤研究，為避免發生大規模采空區塌陷等災害，應對這兩個沉降速率持續增加的區域進行加固處理。