三貴口鉛鋅礦礦柱回收對地表沉陷影響研究

王長軍，盧磊芬，秦仕文，李 輝，溫嘉明，王新龍

（1.低品位難處理黃金資源綜合利用國家重點實驗室，福建 上杭 364204；2.紫金（長沙）工程技術有限公司，湖南 長沙 410012；3.烏拉特后旗紫金礦業有限公司，內蒙古 巴彥淖爾 015500）

0 引言

隨著礦產資源需求量急劇增加，以及礦山開采程度不斷提升，因開采造成的殘留礦柱回采問題愈加突出[1-2]。得益于采礦技術及工藝水平不斷提高，礦柱回收已經成為大多數礦山拓展發展空間、延長生產年限的有效技術手段[3-5]。然而，殘留礦柱回采引起周圍巖體應力變化，以及地表移動沉陷的過程復雜多變，預測和控制因采動造成的地表沉陷和災害至關重要[6]。針對采礦活動引起的地表沉陷問題，相關領域專家學者從不同方面開展大量的研究工作[7-8]。劉寶琛等[9]在隨機介質理論模型的基礎上將概率積分法應用到地下開采引起的地表沉陷預測中，此方法仍是目前開采沉陷預測的主要方法；何國清[10]提出了可預測開采沉陷區非對稱分布位移變化及下沉最大位移的威布爾預計法。同時，學者們還針對典型曲線法、概率積分法、剖面函數法等數學方法、理論模型在開采引起地表沉陷及下沉預測中的應用進行研究[11-14]。近年來，隨著數字化技術水平的快速發展，包括有限元法、離散元法和邊界元法等數值模擬方法逐漸應用于礦山開采地表沉陷等方面，并實現了地表變形沉陷規律及空間分布的精確預測，為前人研究中難以解決的問題提供了新的思路[15-16]。目前，雖然國內外學者針對開采引起的地表沉陷機理及規律進行了一些研究，分析得到了影響地表沉陷的主要因素及方式，建立了預測函數模型，并根據研究成果進行改進，但是由于礦山開采所引起的地表沉陷本身所具有的時間上和空間上多因素影響的復雜性和動態性的特點，為了更好地預測采動引起的地表沉陷及控制沉陷災害，上述地表沉陷預測模型仍需根據不同開采環境和開采條件進一步深入研究。劉玉成[17]針對開采引起的地表動態沉陷及災害預防控制問題，結合彈性薄板的關鍵巖層理論建立地表下沉盆地的力學模型，通過改進時間函數模型分析開采過程中地表沉陷參數的動態變化，分析結果與礦區測量數據相一致；楊曉玉等[18]針對概率積分法計算參數的選取存在的準確性問題，提出用穩健遺傳算法進行參數反演，解決了計算參數抗粗差能力差的問題，保證了計算結果的穩定性和準確度；朱權潔等[19]基于D-InSAR 技術分別構建了SAR 影響數據技術的分析方法和處理流程，以及礦區地表沉降預測模型，根據礦區實測數據對工作面區域對應的地表沉降規律進行分析，驗證了該方法的可靠性，以及對開采可能造成的地表沉陷危害預防提供支撐。

三貴口鉛鋅礦礦區內礦石資源量豐富，礦體數量多，存在分支復合現象，且礦體厚度主要為緩傾斜中厚至厚礦體，空間產出關系復雜，屬于極難采類型的礦體。區域地質構造復雜，裂隙發育，局部含水豐富，圍巖基本穩固，礦體存在厚度、走向、傾向等變化大的特點。目前，三貴口鉛鋅礦采用分段空場法進行開采，采場內遺留大量不規則間柱和頂底柱，嚴重制約礦山規模化生產和穩產目標，采場內遺留礦柱如圖1 所示。

據統計，礦山630 m 中段及以上現存各類礦柱84 個，地質儲量共193.30 萬t，其中，間柱地質儲量162.50 萬t，頂底柱地質儲量30.80 萬t，具體數據見表1。

表1 礦山各中段礦柱分布情況Table 1 Distribution of pillars in each level

1 礦柱回收地表變形及沉陷規律

1.1 地表變形及沉陷基本原理

剖面函數法、典型曲線法、分布函數法和概率積分法等是預測地表變形及沉陷的常用方法。其中，概率積分法將顆粒介質移動看作隨機過程，又可稱為隨機介質理論法[20]。概率積分法在垂直于開采斷面的方向開采斷面為ds×1（寬×高，m×m）的無限長條，開采引起水平表面形成單元下沉盆地，如圖2 所示。

圖2 地下單元開采引起的地表單元與下沉盆地及水平移動曲線Fig.2 Surface units and subsidence basins and horizontal displacement curves due to underground unit mining

試驗證明，地表單元下沉盆地曲線在形狀上與概率密度曲線相似，單元斷面開采所引起的地表單元下沉表達式見式（1）。

式中，r為主要影響半徑，其值與覆巖層性質和開采深度有關。

采用概率積分法計算地表變形及沉陷范圍時，主要對開采所引起的地表下沉量、地表傾斜值、地表曲率值變化、地表水平移動以及地表水平變形等對地表移動和變形沉降有重要影響的參數進行綜合評價[21-23]。

1.2 預測計算參數選取

概率積分法在實際預測地表變形沉陷模型時主要涉及五個參數，體現了上覆巖層結構、巖性、底層傾角、采礦方法、頂板管理辦法等因素的綜合影響[24-25]。

1）下沉系數。地表下沉系數的取值主要由采礦方法、采空區處理方法，以及上覆巖層結構和性質決定。當采用崩落周圍巖體來處理采空區時，因開采引起上部巖層移動與地表變形最為顯著，地表下沉系數一般取0.60～0.90；當采用空場法預留礦柱來處理采空區時，其下沉系數一般取0.05～0.25；當采用充填法處理采空區時，充填的存在降低了巖層下沉時間，下沉系數主要由采空區充實率決定，一般為0.05～0.50。同時，下沉系數和巖層剛度、堅硬程度負相關，和巖層整體破碎程度也符合相同規律。

2）水平移動系數。水平移動系數的取值主要由地層傾角和上覆近地表地層的特性所決定。一般來說，水平移動系數和巖層傾角符合正相關，即巖層傾角越大，水平移動系數越大；當開采近水平礦體時，水平移動系數主要受表土層厚度和性質兩個主要因素決定，即表土層流動性越好，厚度越厚，該系數取值越大，一般為0.20～0.35。

3）主要影響角正切值。根據充分采動條件下的巖層移動及分布規律可知，地表移動變形主要集中發生于開采區域邊界上方，圈定的影響半徑則為變形的主要范圍。該系數隨上覆巖層堅硬程度增加而減小，一般取值范圍在1.5～3.5。

4）開采影響傳播角。開采影響傳播角θ取值主要由上覆巖層性質和巖層傾角所決定，用于表征不同巖層特征開采條件下下沉盆地向下山方向偏移程度，且巖層傾角越大，開采影響傳播角越小。

5）拐點偏移距。拐點偏移距的取值主要取決于上覆巖層的特性和埋深，反映了下沉曲線上拐點在所在礦體平面上偏離開采邊界的距離。巖層傾角為近水平時，最大下沉拐點處于開采區域的上方，且隨開采巖層傾角的變化偏移逐漸偏離開采區域。

三貴口鉛鋅礦礦區內礦巖及主要圍巖的單軸抗壓強度介于35～45 MPa 之間，普氏系數為3.5～4.5，覆巖類型為中硬，多數礦柱深度介于400～550 m 之間（取450 m），礦體傾角取45°。經實測研究，結合巖石力學試驗結果及采空區、礦柱參數，并基于類比法和公式法得出的概率積分法進行預測，所需的主要參數見表2。

表2 概率積分法相關參數取值Table 2 Values of related parameters of probability integral method

1.3 地表變形沉陷范圍計算

針對表1 中礦山各中段礦柱分布情況統計結果，基于面積承載理論計算三貴口鉛鋅礦630 m 及以上中段礦柱安全系數F，見式（2）。

式中：Sp為礦柱強度，MPa；σp為礦柱平均應力，MPa；SL為礦柱的抗壓強度，MPa；W0為礦房寬度，m；Wp為礦柱寬度，m；γ為上覆巖層巖石容重，kN/m3；z為礦柱埋藏深度，m；h為礦柱高度，m；α為常數，其取值與礦柱的寬高比i有關，當i＞5 時，α=1.4，當i≤5，α=1.0。

基于安全系數法評價待回收區域84 個不規則礦柱的穩定性分級結果，如圖3 所示。根據安全系數評價結果，對安全性等級為Ⅰ的39 個礦柱進行開采，各中段礦柱回收結果具體見表3。由表3 可知，由于目前礦山采用空場法開采，每個采場留頂柱、間柱，加之礦山廢石充填進度較慢，使得回采結束的采空區遲遲未充填，造成大量的井下采空區，且隨著時間的延續及下部中段礦體的回采，上部礦柱的回采條件逐漸變差，導致礦柱整體回收率僅為46.43%；同時，隨著回采深度的加深，礦柱所受地應力的影響效應更加明顯，導致部分礦柱風險加大，嚴重影響回采過程中的安全性，在深部礦柱回收的過程中，應加強礦柱周圍穩定性監測及預警措施。

表3 礦山630 m 及以上中段礦柱回采方案Table 3 Pillar recovery scheme for each level above 630 m

圖3 基于安全系數評價法的礦柱穩定性分級結果Fig.3 Pillar stability grading results based on safety coefficient evaluation method

結合表2 概率積分法相關參數取值結果和礦山630 m 及以上中段礦柱回收方案，得出地表變形沉陷影響因素計算結果，見表4。

表4 地表變形沉陷主要影響因素計算結果Table 4 Calculation results of major influential factors of surface subsidence

由表4 可知，地表最大傾斜、地表最大曲率均在地表構建筑物的變形控制要求范圍之內，且其他參數均小于下沉量的影響范圍[26]。

圖4 為基于概率積分法計算的地表移動盆地范圍圖。由圖4 可知，隨著三貴口鉛鋅礦630 m 及以上中段礦柱自上而下進行回收，所引起地表最大下沉量為24.4 mm，且最大下沉量位于礦柱集中區域，該位置上方地表沉陷最為嚴重，但未達到引起地表沉陷災害的范圍，并確定地表移動盆地邊界x坐標：422 125～423 525、y坐標：4 553 322～4 554 450（圖4圈定范圍）；同時，確定地表最大傾斜i0、地表最大曲率值K0、地表最大水平變形值ε0分別為0.114 mm/m、0.000 81 mm/m2、0.052 17 mm/m，均小于采礦相關規定的礦床開采時地表位移與變形沉陷的允許值，使地下開采對地表造成的影響控制在合理范圍之內。

圖4 基于概率積分法計算地表移動盆地范圍Fig.4 The range of surface subsidence basin based on probability integral method

2 數值模擬及地表變形沉陷分析

2.1 模型建立及模擬方案

數值模擬主要以三貴口鉛鋅礦630 m 中段及以上區域為研究對象。由于礦體的幾何形態極其復雜，對不規則礦體應進行簡化處理，并將簡化后的模型利用前處理軟件建立采空區、礦柱及礦體地質模型，并對各中段采空區、礦柱等重點區域進行網格加密處理，最終建立的模型如圖5 所示。

圖5 FLAC3D 數值模擬計算模型Fig.5 Numerical simulation model of FLAC3D

采用Mohr-Coulomb 本構模型，且僅考慮自重作用，礦柱回收范圍之外的其他特殊地質體的影響則不考慮。同時，設置邊界約束條件為：模型四周節點約束水平位移；模型底部所有節點進行垂直位移約束；模型上部地表設置為自由表面。數值模擬過程采用自上而下依次回采，即從880 m 中段向630 m 中段逐步進行回收處理。

2.2 地表沉陷范圍及結果分析

2.2.1 地表下沉結果分析

當地下開采引起上覆巖層變形傳遞到地表時，原有地面標高將沿著某一中心發生方向豎直向下沉降，并逐漸形成大于采空區的下沉區域，該區域范圍即下沉盆地。下沉盆地中地表點的下沉量為地表移動向量的垂直分量，反映地表點在垂直方向的變化量。圖6 為礦柱回收后的地表下沉等值線及下沉盆地范圍圈定圖。由圖6 可知，下沉等值線最大下沉點位于礦柱及采空區密集位置，且地表最大下沉量為36 mm，略大于概率積分法計算所得的最大下沉量24.4 mm，數值模擬與理論計算趨勢變化基本相符。受到礦柱回收的影響，地表移動沉陷趨勢平穩過渡，未出現急劇變化的現象，表明回采過程中地表一致處于基本穩定狀態。

圖6 地表沉陷等值線及下沉盆地圈定Fig.6 Surface subsidence contour line and subsidence basin delineation

圖6（a）等值線圖中下沉量為10 mm 的地表點所形成的邊界曲線為下沉盆地邊界。圖6（b）為下沉盆地圈定范圍。由圖6（b）可知，下沉盆地為沿空區和礦柱分布的不規則橢圓形，區域范圍為x方向：422 125～423 550、y方向：4 553 330～4 554 472，該區域范圍與理論計算結果基本一致，驗證了理論計算與數值模擬結果的一致性。

2.2.2 地表水平移動和變形結果分析

圖7、圖8 分別為地表水平移動和水平變形的數值模擬結果。由圖7（a）可知，地表x方向（傾向）水平位移大致沿礦體走向對稱分布，兩側各形成水平位移中心，東西兩側水平位移最大數值均為8 mm，西側為正、東側為負，整體向遠離礦柱和采空區的區域偏移，符合礦柱開采特征規律。由圖7（b）可知，地表y方向（走向）水平位移大致沿礦體傾向對稱分布，兩側各形成水平位移中心，南北兩側水平位移最大數值分別為12 mm 和8 mm，南側為正、北側為負，整體也向遠離礦柱和采空區的區域偏移。

圖7 地表水平移動數值模擬等值線圖Fig.7 Numerical simulation contour line of surface horizontal movement

圖8 地表水平變形數值模擬等值線圖Fig.8 Numerical simulation contour line of surface horizontal deformation

由數值模擬計算結果可知，三貴口鉛鋅礦礦柱待回收區域整體水平位移較小，略大于概率積分法計算所得的最大水平移動理論值，數值模擬結果與理論計算結果大致相符。由于礦柱回收與采空區的存在，地表對應680 m 中段、730 m 中段、780 m 中段礦柱密集處的水平位移等值線密度較大，但由于水平位移值較小，地表最大下沉點位置變化不大。

由圖8 可知，x方向水平變形εx隨該方向處水平移動值Ux變化，沿走向存在2 條0 變形值線，0 變形值線周圍水平變形方向相反，其最大值為0.000 23 mm/m；而y方向水平變形εy隨該方向處水平移動值Uy變化，沿傾向存在2 條0 變形值線，0 變形值線周圍水平變形方向相反，其最大值為0.000 80 mm/m。由地表水平變形模擬結果可知，通過回收穩定性分級結果為Ⅰ級礦柱方案，地表變形沉陷區域的水平變形遠小于地表的安全臨界值；同時，模擬結果表明水平變形在礦柱回收過程中發展較為穩定，回采結束后地表移動盆地范圍內的水平變形逐漸趨于穩定。

2.2.3 地表傾斜及曲率分析

圖9、圖10 分別為地表傾斜數值模擬結果和地表曲率數值模擬結果。由圖9 可知，地表傾斜ix曲線等值線與iy曲線等值線隨水平移動曲線變化，其中，ix曲線沿礦體走向對稱，在近礦柱、采空區處上方地表出現極值，西側極小值為-0.000 1 mm/m，東側極大值為0.000 1 mm/m；iy曲線沿礦體傾向對稱，在近礦柱、采空區處上方地表出現極值，北側極大值為7×10-5mm/m，南側極小值為-9×10-5mm/m。結合下沉量等值線及傾斜等值線可知，礦柱回采后地表沿10 mm 下沉量等值線形成下沉盆地，且盆地傾斜率較小，較為平緩。模擬結果表明，地下礦柱回收及采空區的處理方案所引起的地表傾斜不會對地表造成影響，圈定的范圍能夠為工程提供支撐作用。

圖9 地表傾斜值數值模擬等值線圖Fig.9 Numerical simulation contour line of surface tilt value

圖10 地表曲率值數值模擬等值線圖Fig.10 Numerical simulation contour line of surface curvature value

由圖10 可知，研究范圍內地表點曲率Kx基本接近于0，僅部分位置處（礦柱、采空區密集處）出現數量級為10-6mm/m2的曲率極值；而研究范圍內地表點曲率Ky基本接近于0，僅部分位置處（礦柱、采空區密集處）出現數量級為10-6mm/m2的曲率極值。結合圖9（a）和圖9（b）可知，曲率Kx和曲率Ky的極值位置和數量級基本重合，均為礦柱、采空區密集處。模擬結果表明，地表曲率變形值遠未達到研究范圍規定的臨界值，充分說明在礦柱回采期間所引起的巖層移動及地表沉陷處于相對穩定狀態，不會因礦柱回收而對地面構筑物產生較大影響。

由數值模擬結果可知，與概率積分法計算結果相比，數值模擬的下沉量和水平移動值基本一致，且二者所確定的地表影響范圍也大致相同；而數值模擬中傾斜、曲率和水平變形相對于理論分析結果則較小，但都處于規定的合理范圍之內，不會因礦柱回收產生較大的地表變形沉陷。

3 結論

研究三貴口鉛鋅礦630 m 中段及以上礦柱回收后覆巖及地表的變形及沉陷規律，利用概率積分法和FLAC3D軟件分析預測地表的變形及沉陷分布情況，綜合地表下沉量、水平移動、傾斜、曲率、水平變形等評價指標系統研究礦柱回收對地表的沉陷影響，主要結論如下所述。

1）采用概率積分法計算地表理論位移變形結果：地表最大下沉值W0=24.4 mm；地表最大傾斜值i0=0.114 mm/m；地表最大曲率值K0=±0.000 81 mm/m2；地表最大水平移動值U0=7.32 mm；地表最大水平變形值ε0=0.052 17 mm/m，均在相關規范要求的合理范圍之內。

2）采用FLAC3D模擬分析礦柱回采后地表的位移變形，數值模擬結果顯示地表最大下沉量為36 mm，位于采空區及礦柱回收密集區域上方地表，且該區域水平移動也符合相同規律；地表傾斜ix曲線與iy曲線隨水平移動曲線變化，分別沿礦體走向和傾向對稱，且均在近礦柱、采空區處上方地表出現極值；地表點曲率基本接近于0，僅礦柱、采空區密集處出現數量級為10-6mm/m2的極值點；水平變形在x方向、y方向隨水平移動值而變化，最大值分別為0.000 23 mm/m 和0.000 80 mm/m。

3）結合概率積分法與數值模擬結果可知，理論和模擬的結果中地表下沉量和水平移動值基本一致。理論計算或數值模擬中地表最大傾斜遠小于1 mm/m、最大曲率遠小于0.05 mm/m2、最大水平變形遠小于1 mm/m，都處于規定的合理范圍之內，不會因礦柱回收產生較大的地表變形沉陷造成地表塌陷等災害的發生。