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

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

0 引言

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

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

據(jù)統(tǒng)計，礦山630 m 中段及以上現(xiàn)存各類礦柱84 個，地質(zhì)儲量共193.30 萬t，其中，間柱地質(zhì)儲量162.50 萬t，頂?shù)字刭|(zhì)儲量30.80 萬t，具體數(shù)據(jù)見表1。

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

1 礦柱回收地表變形及沉陷規(guī)律

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

剖面函數(shù)法、典型曲線法、分布函數(shù)法和概率積分法等是預(yù)測地表變形及沉陷的常用方法。其中，概率積分法將顆粒介質(zhì)移動看作隨機過程，又可稱為隨機介質(zhì)理論法[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為主要影響半徑，其值與覆巖層性質(zhì)和開采深度有關(guān)。

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

1.2 預(yù)測計算參數(shù)選取

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

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

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

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

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

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

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

表2 概率積分法相關(guān)參數(shù)取值Table 2 Values of related parameters of probability integral method

1.3 地表變形沉陷范圍計算

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

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

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

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

圖3 基于安全系數(shù)評價法的礦柱穩(wěn)定性分級結(jié)果Fig.3 Pillar stability grading results based on safety coefficient evaluation method

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

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

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

圖4 為基于概率積分法計算的地表移動盆地范圍圖。由圖4 可知，隨著三貴口鉛鋅礦630 m 及以上中段礦柱自上而下進行回收，所引起地表最大下沉量為24.4 mm，且最大下沉量位于礦柱集中區(qū)域，該位置上方地表沉陷最為嚴重，但未達到引起地表沉陷災(zāi)害的范圍，并確定地表移動盆地邊界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，均小于采礦相關(guān)規(guī)定的礦床開采時地表位移與變形沉陷的允許值，使地下開采對地表造成的影響控制在合理范圍之內(nèi)。

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

2 數(shù)值模擬及地表變形沉陷分析

2.1 模型建立及模擬方案

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

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

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

2.2 地表沉陷范圍及結(jié)果分析

2.2.1 地表下沉結(jié)果分析

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

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

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

2.2.2 地表水平移動和變形結(jié)果分析

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

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

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

由數(shù)值模擬計算結(jié)果可知，三貴口鉛鋅礦礦柱待回收區(qū)域整體水平位移較小，略大于概率積分法計算所得的最大水平移動理論值，數(shù)值模擬結(jié)果與理論計算結(jié)果大致相符。由于礦柱回收與采空區(qū)的存在，地表對應(yīng)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。由地表水平變形模擬結(jié)果可知，通過回收穩(wěn)定性分級結(jié)果為Ⅰ級礦柱方案，地表變形沉陷區(qū)域的水平變形遠小于地表的安全臨界值；同時，模擬結(jié)果表明水平變形在礦柱回收過程中發(fā)展較為穩(wěn)定，回采結(jié)束后地表移動盆地范圍內(nèi)的水平變形逐漸趨于穩(wěn)定。

2.2.3 地表傾斜及曲率分析

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

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

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

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

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

3 結(jié)論

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

1）采用概率積分法計算地表理論位移變形結(jié)果：地表最大下沉值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，均在相關(guān)規(guī)范要求的合理范圍之內(nèi)。

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

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