盧欣奇，李學(xué)峰，張勤斌，黃海斌

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

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

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

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

1 PS-InSAR技術(shù)

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

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

2 實驗區(qū)概況與數(shù)據(jù)處理

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

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

圖1 PS-InSAR技術(shù)數(shù)據(jù)處理流程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矢量點信息導(dǎo)入ArcGIS進行坐標變換、重新投影坐標系轉(zhuǎn)化、克里金面生成、與采空區(qū)分布平面圖進行疊加，最終將PS點位信息可視化，得到研究區(qū)域內(nèi)采空區(qū)的沉降量分布云圖、平均沉降速率云圖。

3 試驗結(jié)果與討論

3.1 試驗結(jié)果

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

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

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

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

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

3.2 剖面線情況討論

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

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

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

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

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

由圖6可知，6號采空區(qū)剖面線西側(cè)較為平緩東側(cè)較陡且漏斗面較為平滑，整體呈現(xiàn)“V”形沉降梯度分布，由采空區(qū)沉降量分布云圖可知，該沉降區(qū)域的最大沉降點F為6號采空區(qū)的東南方向，應(yīng)屬于6號采空區(qū)頂板的薄弱部位。在最大沉降點西側(cè)平緩部位，并未出現(xiàn)明顯漏斗區(qū)域，平均沉降量為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點的沉降速率出現(xiàn)了三次峰值。該點沉降速率一直緩慢增加，達到極值后逐漸下降，持續(xù)一定時間又出現(xiàn)了加速狀態(tài)，沉降速率均勻，在達到第二次峰值后又出現(xiàn)了沉降速率減小的情況。大體上呈現(xiàn)“平衡狀態(tài)-非平衡狀態(tài)-新平衡狀態(tài)”這一開采沉陷規(guī)律，當(dāng)采空區(qū)頂板在沉陷過程中達到平衡狀態(tài)后，又會出現(xiàn)新的失穩(wěn)狀態(tài)，當(dāng)頂板及礦柱經(jīng)歷過失穩(wěn)過程后又一次達到平衡時即進入到新的平衡狀態(tài)。

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

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

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

4 數(shù)值模擬驗證

4.1 巖體力學(xué)參數(shù)及礦柱平均強度

巖體力學(xué)參數(shù)是依據(jù)巖石力學(xué)參數(shù)特性測試結(jié)果，并考慮了巖體的結(jié)構(gòu)效應(yīng)、地下水、節(jié)理裂隙等因素，對巖石力學(xué)參數(shù)按照POPOV的RMR巖體質(zhì)量分類法進行適當(dāng)?shù)男拚齕16]。根據(jù)長沙礦山研究院在“金屬礦山安全技術(shù)國家重點實驗室”對礦體及圍巖采用MTS實驗系統(tǒng)進行實驗得出的巖石力學(xué)參數(shù)，選用飽和巖樣實驗結(jié)果，使用加拿大Rocscience公司開發(fā)的Roclab1.0軟件計算礦體及頂板、底板圍巖的巖體力學(xué)參數(shù)，從而確定基于數(shù)值模擬的礦山采空區(qū)穩(wěn)定性分析所需的巖體力學(xué)參數(shù)，巖體及后期采用的充填體力學(xué)參數(shù)見表2。

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

(1)

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

各采空區(qū)的平均礦柱強度見表3。

表2 礦巖介質(zhì)的力學(xué)參數(shù)Table 2 Mechanical parameters of rock medium

表3 平均礦柱強度統(tǒng)計表Table 3 StatisticalTable of average strength of pillar

4.2 計算分析

4.2.1 最大主應(yīng)力分析

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

圖8所示為各采空區(qū)頂板及礦柱縱剖面、橫剖面的最大主應(yīng)力云圖。

為確保采空區(qū)經(jīng)治理后達到長期穩(wěn)定，礦柱安全系數(shù)采用保守方式計算，計算公式見式(2)。

F=QP/σPmax

(2)

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

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

n=1+t×σ

(3)

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

采空區(qū)頂板、礦柱的應(yīng)力狀態(tài)、不穩(wěn)定礦柱個數(shù)等信息詳見表4。

表4 采空區(qū)應(yīng)力狀態(tài)統(tǒng)計表Table 4 StatisticalTable of stress state in goaf

圖8 最大主應(yīng)力云圖Fig.8 Maximum principal stress nephogram

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

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

4.2.2 最小主應(yīng)力分析

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

4.3 結(jié)果對比

PS-InSAR測量得出的沉降量較大區(qū)域與數(shù)值模擬得出的薄弱區(qū)域?qū)φ招畔⒃斠姳?。

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

圖9 最小主應(yīng)力云圖Fig.9 Minimum principal stress nephogram

表5 數(shù)值模擬及實測沉降區(qū)域?qū)φ毡鞹able 5 Numerical simulation and measured settlement area check list

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

5 結(jié) 論

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

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

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

4) 3號采空區(qū)與6號采空區(qū)交接位置受到開采擾動影響出現(xiàn)沉降速率較快增加的現(xiàn)象，達到峰值后整體沉降趨勢較為穩(wěn)定，但在數(shù)據(jù)截至?xí)r沉降速率依舊在加快。

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

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