蘭 洋，袁利偉，李素敏，馮家寧，劉大榮，王 華

(1.昆明理工大學(xué) 公共安全與應(yīng)急管理學(xué)院，云南 昆明 650093；2.昆明理工大學(xué) 國土資源工程學(xué)院，云南 昆明 650093；3.昭通永善金沙鉛鋅礦，云南 昭通 657000)

0 引言

因采空區(qū)導(dǎo)致的地質(zhì)災(zāi)害事故頻發(fā)，采空區(qū)穩(wěn)定性研究越來越受到專家學(xué)者的關(guān)注。AGIOUTANTIS等[1]運用有限元法和邊界元法研究了采動影響下覆巖垮落的開采條件和垮落高度、覆巖產(chǎn)生離層裂縫的力學(xué)條件及離層裂縫的位置和高度；李夕兵等[2]對金屬礦地下采空區(qū)的探測、處理與安全評判進行了研究；JAISWAL等[3]為解決采空區(qū)的頂板控制問題，在開采前對其結(jié)構(gòu)進行了模擬，并對地層參數(shù)進行了校準。萬文[4]綜合運用遺傳算法、強度折減法和FLAC3D軟件研究了采空區(qū)對邊坡穩(wěn)定性的影響；劉曉明等[5]使用CMS激光對采空區(qū)進行了探測，利用獲取的采空區(qū)三維邊界數(shù)據(jù)，通過自編的程序接口導(dǎo)入FLAC3D軟件進行采空區(qū)圍巖穩(wěn)定性分析。

德國、波蘭等國家率先將InSAR技術(shù)應(yīng)用于礦區(qū)開采引起的地表沉陷監(jiān)測。我國在InSAR監(jiān)測地表變形方面的研究起步較晚，且前期主要將其用于城市地表沉陷監(jiān)測。張學(xué)東等[6]利用SBAS-InSAR技術(shù)生成差分干涉圖，以監(jiān)測地面沉降；獨知行等[7]綜合比較了傳統(tǒng)測量與InSAR測量的優(yōu)勢和缺陷，提出了將傳統(tǒng)測量手段與InSAR技術(shù)結(jié)合共同監(jiān)測礦區(qū)地表沉陷的構(gòu)想；陳繼偉等[8]利用Sentinel-1A TOPS 模式影像對黃河三角洲地區(qū)進行了變形監(jiān)測，證明了TOPS 模式 SBAS 時序分析方法應(yīng)用于變形監(jiān)測的可行性和有效性。

目前多采用單一方式對采空區(qū)穩(wěn)定性進行評價，缺乏系統(tǒng)的綜合監(jiān)測。基于此，本文擬采用InSAR與FLAC3D協(xié)同分析采空區(qū)的穩(wěn)定性。

1 工程概況

永善縣金沙廠207鉛鋅礦體(簡稱“207礦體”)處于金沙鉛鋅礦官房礦段西北部，形態(tài)呈似層狀，總體傾向305°～46°，傾角11°～38°；礦體分布標高1 005～1 062 m，礦體厚3.04～6.80 m，平均厚5.94 m。

207礦體采礦方法為房柱法，沿礦體走向布置連續(xù)礦塊，每個中段高20 m，礦塊寬27 m。中段礦柱、礦塊間柱和采場頂(底)柱尺寸均為4 m，礦房寬23 m，礦房高為礦體垂直厚度。隨著開采的逐步推進，PD207、PD207-1坑道之間形成了采空區(qū)。經(jīng)過多年回采，采空區(qū)的形態(tài)不斷變化，空區(qū)規(guī)模逐漸擴大，其中一個采空區(qū)的最大跨度達67 m，采空區(qū)頂板暴露面積超過2 000 m2。采空區(qū)統(tǒng)計見表1，采空區(qū)位置見圖1。

表1 采空區(qū)統(tǒng)計

由于采空區(qū)和上覆巖層在地應(yīng)力作用下發(fā)生了移動[9]，導(dǎo)致207礦體上方出現(xiàn)了一定程度的地表變形，可能會引發(fā)地面塌陷、滑坡等災(zāi)害。因此，對207礦體地表進行了InSAR變形監(jiān)測。以開采境界為依據(jù)，對礦區(qū)范圍內(nèi)受地下開采影響的地表區(qū)域進行變形監(jiān)測。

圖1 采空區(qū)位置

2 InSAR地表變形監(jiān)測

2.1 InSAR技術(shù)原理

星載合成孔徑雷達(InSAR)可實現(xiàn)對礦山采空區(qū)沉陷、邊坡滑坡等重大災(zāi)害的遠程監(jiān)測、監(jiān)控、排查、預(yù)警等[10]。SBAS-InSAR監(jiān)測技術(shù)原理是將SAR影像在一定時間和空間上用很多小的集合來表示，然后采用最小二乘法將每一個小集合的地表相位進行解算，采用SVD聯(lián)合求解的方式得到最小二乘解[11-12]。

本文通過小基線集(SBAS-InSAR)技術(shù)對礦區(qū)地表進行變形監(jiān)測，具有能克服時空基線失相干問題、顯著提高影像利用率等優(yōu)點[13]。其主要通過對衛(wèi)星參數(shù)及研究區(qū)地形參數(shù)進行雷達可視區(qū)劃分，并結(jié)合SBAS-InSAR技術(shù)反演得到研究區(qū)的垂直變形場，進一步獲取有效的觀測區(qū)域，再結(jié)合礦區(qū)與地表的空間關(guān)系，綜合分析識別潛在的風險區(qū)域。

2.2 數(shù)據(jù)選擇與處理

本文采用的影像數(shù)據(jù)是Sentinel-1(哨兵1號)衛(wèi)星獲取的時間跨度為2018年3月22日-2020年9月19日的升軌Sentinel-1A數(shù)據(jù)。為了縮短數(shù)據(jù)處理時間，需對原始數(shù)據(jù)進行裁剪，裁去非研究區(qū)域；同時在數(shù)據(jù)裁剪時，為了提高數(shù)據(jù)處理的可靠性，應(yīng)適當擴大研究區(qū)范圍。為了更好地顯示礦區(qū)每年的地表形態(tài)變化，3年的數(shù)據(jù)分3次處理。

2.3 數(shù)據(jù)分析

通過對3年礦區(qū)變形量圖的分析，獲得207礦體上方地表沉降量分布云圖(見圖2-圖4)。由圖2-圖4可以看出：207礦體上方地表沉降量呈逐年增大的趨勢，且采空區(qū)上方沉降量明顯更大；隨著時間的推移，2#和3#采空區(qū)上方地表出現(xiàn)了較大變形，2020年其沉降量達30 mm。

圖2 2018年3月-2019年1月207礦體地表沉降量分布云圖

圖3 2019年1月-2020年1月207礦體地表沉降量分布云圖

圖4 2020年1月-2020年9月207礦體地表沉降量分布云圖

通過對InSAR監(jiān)測數(shù)據(jù)的分析，得到了整個207礦體上方地表的變形情況，識別其變形較嚴重的區(qū)域為2#、3#采空區(qū)上方地表。因此，在地表選取1-1′、2-2′兩條剖面線(見圖1)，在剖面線上每隔10 m選取一個特征點，通過對每個特征點的沉降時間序列圖計算累計沉降量，擬合各特征點累計沉降量便可得到地表剖面線的沉降量曲線(見圖5、圖6)。

圖5 1-1′剖面沉降量曲線

圖6 2-2′剖面沉降量曲線

由圖5、圖6可知，地表沉降量最大值為38 mm，在剖面線上呈現(xiàn)出不均勻的凹形沉降，局部下沉非常明顯。地表沉降區(qū)域位于2#、3#采空區(qū)上方，推測是地下開挖導(dǎo)致的地表沉降所致，故對2#、3#采空區(qū)開展采空區(qū)穩(wěn)定性分析。

3 采空區(qū)穩(wěn)定性分析

采用Midas/GTX建立地表和地下采空區(qū)三維可視化模型，然后將模型的單元、節(jié)點和組信息導(dǎo)入FLAC3D定義本構(gòu)模型和材料參數(shù)并施加載荷及邊界條件[14-15]。

3.1 采空區(qū)三維模型的建立

首先用Midas/GTX將AutoCAD地表等高線生成地表模型，然后確定2#、3#采空區(qū)與地表的空間關(guān)系，生成對應(yīng)的地下采空區(qū)，建立實體模型。對實體模型進行分組，劃分網(wǎng)格。

為提高計算效率，需對實際條件作適當簡化：視緩傾斜的207礦體為平行于XY的平面，2#、3#采空區(qū)存在于礦體層。采空區(qū)高度取平均值5.5 m，頂板暴露面積按實際面積確定。最終構(gòu)建了X軸長300 m，Y軸長260 m的實體模型。采空區(qū)單元尺寸為3 m×3 m×3 m，采空區(qū)礦柱單元尺寸為2 m×2 m×2 m，圍巖單元尺寸為5 m×5 m×5 m；模型共劃分為91 174個單元、87 538個節(jié)點(見圖7)。

圖7 三維實體模型

3.2 數(shù)值分析模型

將實體模型導(dǎo)入FLAC3D，綜合分析采空區(qū)的應(yīng)力、位移、塑性區(qū)的分布規(guī)律[16-17]。數(shù)值模擬模型如圖8所示。

圖8 數(shù)值模擬模型

1)模型基本假設(shè)

視所有巖體為各向同性的連續(xù)均勻介質(zhì)，且僅受重力作用，忽略節(jié)理、裂隙、斷裂帶等因素對巖體穩(wěn)定性的影響。

2)選取力學(xué)模型

摩爾-庫倫模型適用于在剪應(yīng)力下屈服但剪應(yīng)力僅取決于最大、最小主應(yīng)力，而第二主應(yīng)力對屈服不起作用的材料[18]。本文選用摩爾-庫倫模型模擬采空區(qū)開挖。

3)巖體力學(xué)參數(shù)

根據(jù)巖石力學(xué)試驗結(jié)果確定各項物理力學(xué)參數(shù)，并采用Hoek-Brown強度折減法對參數(shù)進行折減[19-20]。最終計算模型所采用的巖體力學(xué)參數(shù)見表2。

表2 巖體力學(xué)參數(shù)

4)邊界條件和初始條件

計算模型四周采用水平約束，限制模型側(cè)面橫向位移；底部采用固定約束，限制模型豎向位移；頂部為自由面。設(shè)置重力加速度為-9.82 m/s2，指向Z軸負方向。

3.3 數(shù)值模擬結(jié)果分析

分別對2#、3#采空區(qū)頂板位移、最大主應(yīng)力、壓應(yīng)力和剪應(yīng)力分布情況進行分析，以判斷其是否存在失穩(wěn)風險。將數(shù)值模擬得到的地表沉降過程曲線與InSAR監(jiān)測得到的沉降量曲線進行對比，以判斷采空區(qū)是否會繼續(xù)引發(fā)地表沉降。

3.3.1 監(jiān)測與數(shù)值模擬協(xié)同分析

根據(jù)1-1′、2-2′剖面線上的特征點在計算模型中布置相應(yīng)的位移監(jiān)測點。通過擬合數(shù)值模擬計算結(jié)果得到1-1′、2-2′剖面沉降過程曲線(見圖9)，并與InSAR監(jiān)測的剖面沉降過程曲線(見圖10)作對比分析。

圖9 1-1′剖面沉降過程曲線

圖10 2-2′剖面沉降過程曲線

由圖9、圖10可知：2#、3#采空區(qū)在開挖初期就對地表產(chǎn)生了影響，導(dǎo)致地表沉降；在數(shù)值模擬計算結(jié)果逐漸穩(wěn)定的過程中，采空區(qū)上方巖體變形和地表變形也逐漸趨于穩(wěn)定；在采空區(qū)開挖的初期和中期，采空區(qū)對地表產(chǎn)生的沉降趨勢與InSAR監(jiān)測的沉降趨勢具有一致性，但隨著數(shù)值模擬計算結(jié)果穩(wěn)定以后，2#、3#采空區(qū)導(dǎo)致地表變形的沉降量明顯大于InSAR監(jiān)測結(jié)果，表明采空區(qū)沉降尚未穩(wěn)定，可能會導(dǎo)致地表沉降量繼續(xù)加大。

3.3.2 位移分析

相關(guān)研究表明，采空區(qū)頂板位移會隨著采空區(qū)跨度的增大而增大。因此，在完成數(shù)值模擬計算后，對2個采空區(qū)最大跨度方向取剖面，根據(jù)位移分布云圖判斷采空區(qū)的穩(wěn)定性。

根據(jù)采場開挖后在無支護條件或臨時支護條件下的現(xiàn)場經(jīng)驗[23]，得出頂?shù)装鍑鷰r位移和采場穩(wěn)定性的對應(yīng)關(guān)系(見表3)。

表3 圍巖位移和采場穩(wěn)定性的對應(yīng)關(guān)系

數(shù)值模擬計算后得到的采空區(qū)頂板位移分布云圖如圖11所示。由圖11可以看出，2#、3#采空區(qū)直接頂板的位移分別為30～60 mm、30～70 mm，最大位移集中在采空區(qū)跨度中心區(qū)域，與工程經(jīng)驗相符。根據(jù)最大容許位移判斷，當頂板位移大于50 mm時，采空區(qū)由穩(wěn)定狀態(tài)變?yōu)椴环€(wěn)定狀態(tài)，存在失穩(wěn)風險。2#采空區(qū)直接頂板最大位移達到了60 mm，3#采空區(qū)直接頂板最大位移達到了70 mm，均超過了50 mm。因此，判斷2個采空區(qū)均存在因頂板位移過大而失穩(wěn)的風險。

(a)2#采空區(qū)縱剖面位移云圖

(b)3#采空區(qū)縱剖面位移云圖

3.3.3 剪應(yīng)力及壓應(yīng)力分析

Mohr-Coulomb強度準則認為巖體材料會因受剪應(yīng)力作用而破壞，且其破壞所需剪應(yīng)力和作用在破壞面上的正應(yīng)力之間存在函數(shù)關(guān)系。該強度準則表達式為

(1)

式中：σ1、σ3分別表示巖體發(fā)生破壞時的最大主應(yīng)力、最小主應(yīng)力，MPa；C表示黏聚力，MPa；φ表示內(nèi)摩擦角。

當作用在礦柱上的最大剪應(yīng)力和最大壓應(yīng)力大于其抗剪、抗壓強度時，礦柱存在失穩(wěn)風險。最大剪應(yīng)力和最大壓應(yīng)力模擬結(jié)果分別見圖12、圖13。由圖12、圖13可知，剪應(yīng)力和壓應(yīng)力主要分布在礦柱上，且隨著礦柱尺寸的減小而增大。2#、3#采空區(qū)礦柱上的剪應(yīng)力為5.50～8.13 MPa，壓應(yīng)力為12.00～23.63 MPa。根據(jù)表2和式(1)可得礦柱的抗剪強度和抗壓強度分別為7.29、18.65 MPa，采空區(qū)部分尺寸較小的礦柱上的剪應(yīng)力和壓應(yīng)力超過了其抗剪、抗壓強度，存在失穩(wěn)風險。當某個礦柱失穩(wěn)時，原有載荷將由其相鄰礦柱承擔，使得相鄰礦柱上的剪應(yīng)力和壓應(yīng)力極易超過其自身強度，繼而誘發(fā)大面積的礦柱破壞，因此2#、3#采空區(qū)存在礦柱失穩(wěn)風險。

圖12 剪應(yīng)力分布云圖

圖13 壓應(yīng)力分布云圖

4 結(jié)論

a.利用InSAR技術(shù)監(jiān)測金沙鉛鋅礦207礦體上方2018-2020年的地表沉降發(fā)現(xiàn)，2#、3#采空區(qū)上方地表的沉降量較大，且出現(xiàn)了不均勻沉降，存在采空區(qū)失穩(wěn)風險；1#、4#采空區(qū)上方地表沉降量較小且均勻，暫無失穩(wěn)風險。

b.數(shù)值模擬與InSAR監(jiān)測協(xié)同分析結(jié)果顯示，數(shù)值模擬中的2#、3#采空區(qū)導(dǎo)致的地表沉降量明顯大于InSAR監(jiān)測結(jié)果，說明地表沉降量還可能會繼續(xù)增大。

c.數(shù)值模擬結(jié)果顯示，2#、3#采空區(qū)頂板最大位移分別達到了60、70 mm，且采空區(qū)內(nèi)部尺寸較小的礦柱上的剪應(yīng)力和壓應(yīng)力超過了其抗剪、抗壓強度，故2#、3#采空區(qū)存在失穩(wěn)風險。

d.建議在后期采用廢石充填的方式對2#、3#采空區(qū)進行治理，控制周圍爆破作業(yè)的一次起爆藥量，加強地表沉降監(jiān)測，一旦沉降量出現(xiàn)突變，應(yīng)立即停止采空區(qū)周圍采場的回采作業(yè)。