基于改進(jìn)組合彈簧模型的礦井地應(yīng)力場計算方法

趙毅鑫, 趙良辰, 楊東輝, 邊 華, 楊 哲, 宮智馨

(1.中國礦業(yè)大學(xué)(北京) 共伴生能源精準(zhǔn)開采北京市重點實驗室, 北京 100083; 2.中國礦業(yè)大學(xué)(北京) 煤礦災(zāi)害預(yù)防與處置應(yīng)急管理部重點實驗室, 北京 100083; 3.中國礦業(yè)大學(xué)(北京) 能源與礦業(yè)學(xué)院, 北京 100083; 4.山西大同大學(xué) 煤炭工程學(xué)院, 山西 大同 037003)

地應(yīng)力場是影響巖層運移和原位巖體力學(xué)性質(zhì)的重要因素, 對采礦、隧道開挖、油氣開采等影響顯著。地應(yīng)力是采礦工程中覆巖運動的力源, 在構(gòu)造應(yīng)力強(qiáng)烈的地區(qū)進(jìn)行采礦活動, 可能會引起沖擊地壓、巖爆等災(zāi)害, 影響礦山的安全生產(chǎn)。因此準(zhǔn)確獲取地應(yīng)力場對于提高開采設(shè)計的科學(xué)性及保障礦山生產(chǎn)安全有著重要意義。

要獲得區(qū)域地應(yīng)力場需基于地應(yīng)力測定獲取部分點位的地應(yīng)力值及方向。目前, 常用的地應(yīng)力測量方法主要包括: 水壓致裂法[1]、套孔應(yīng)力解除法[2]、非彈性恢復(fù)法[3]和聲發(fā)射法[4]等。其中, 聲發(fā)射法通過鉆取的巖芯在實驗室開展測試即可獲得巖芯所對應(yīng)位置的三維地應(yīng)力, 較其他測量方法具有測量深度大[5]、成本低和效率高等特點。聲發(fā)射法通過識別巖石Kaiser效應(yīng)點確定原巖應(yīng)力狀態(tài),該方法由KANAGAWA等[6]于1977年首次提出。聲發(fā)射法測定地應(yīng)力有單軸測試與三軸測試兩種方式。單軸測試方法較為簡便, 但試樣亦會在出現(xiàn)Kaiser效應(yīng)點前發(fā)生破壞[7]; 而三軸測試方法能相對準(zhǔn)確地測定地應(yīng)力[8–9]。

然而, 地應(yīng)力測量因測點數(shù)或樣品數(shù)有限, 對礦井區(qū)域地應(yīng)力場的描述仍需結(jié)合應(yīng)力場計算來開展。目前常用的地應(yīng)力場計算方法包括: 估算法、正演算法及反演算法[10]。估算法是基于地應(yīng)力計算模型對地應(yīng)力場進(jìn)行計算。具有代表性的地應(yīng)力計算模型包括金尼克模型[11]、Anderson修正模型[12]、黃榮樽模型[13]、組合彈簧模型[14]和葛洪魁模型[15]等。正演算法主要指通過調(diào)整邊界荷載或位移, 最小化測點處計算值與實測值間的差異, 從而計算出所需區(qū)域地應(yīng)力場的方法[16], 其主要包括邊界載荷調(diào)整法和邊界位移調(diào)整法等[17]。反演算法主要指通過假設(shè)模型邊界荷載類型與分布, 建立邊界荷載與測點應(yīng)力間的回歸關(guān)系, 通過回歸計算, 調(diào)整待定系數(shù)使回歸方程誤差最小以獲得最優(yōu)的模型邊界載荷, 從而獲取所需區(qū)域地應(yīng)力場的方法[18]。反演算法一般包括線性擬合法[19]、非線性擬合法[20]和神經(jīng)網(wǎng)絡(luò)法[21–22]等。

筆者以神東礦區(qū)布爾臺礦42105工作面為工程背景, 結(jié)合對布爾臺礦5個鉆孔中89個地應(yīng)力測點的Kaiser法實測結(jié)果, 提出一種基于改進(jìn)組合彈簧模型的區(qū)域地應(yīng)力場計算方法。該方法考慮了實測數(shù)據(jù)的變化規(guī)律和巖層物性對區(qū)域地應(yīng)力計算中的控制作用, 提高了反演精度和可靠性。同時,將區(qū)域地應(yīng)力反演方法融入FLAC3D軟件建模, 對布爾臺礦42105工作面開采過程中礦壓顯現(xiàn)規(guī)律進(jìn)行了數(shù)值模擬并與實測數(shù)據(jù)進(jìn)行對比, 進(jìn)而證明了所提出方法的優(yōu)越性。

1 研究區(qū)域與地應(yīng)力測量

1.1 研究區(qū)域與取樣

布爾臺礦是神東煤炭集團(tuán)主力生產(chǎn)礦井之一,位于內(nèi)蒙古自治區(qū)鄂爾多斯市伊金霍洛旗烏蘭木倫鎮(zhèn)。大地構(gòu)造分區(qū)屬于華北地臺鄂爾多斯臺向斜伊盟隆起中東部, 伊陜斜坡東北部?；緲?gòu)造形態(tài)為一南西傾斜的近水平產(chǎn)狀的單斜構(gòu)造, 井田構(gòu)造簡單, 層位穩(wěn)定。

試驗巖芯取自布爾臺礦的白堊系及侏羅系地層, 取芯鉆孔包括BK212, BK209, BK213, BK207和BK220, 共計5個鉆孔。圖1為現(xiàn)場取芯鉆孔位置及巖層綜合柱狀圖等。將取芯鉆孔巖樣沿垂直于其軸線方向再按逆時針0°, 45°和90°鉆取尺寸為φ25 mm×50 mm的巖芯。

圖1 布爾臺礦區(qū)域位置及取芯情況Fig.1 Location of Buertai Mine in Shendong mining area

1.2 基于Kaiser效應(yīng)的地應(yīng)力測量

試驗采用WDW–100E型試驗機(jī)加載系統(tǒng)對巖樣施加載荷, 采用美國物理聲學(xué)公司(PAC)生產(chǎn)的PCI–II型AE監(jiān)測系統(tǒng)監(jiān)測聲發(fā)射信號。WDW–100E試驗機(jī)提供的最大軸向載荷量程為100 kN, 試驗力測量范圍為0.4%~100%FN, 變形測量范圍為2%~100%FS。試驗力測量精度、位移測量準(zhǔn)確度及變形示值相對誤差均為±1%以內(nèi)。PCI–II型AE監(jiān)測系統(tǒng)具有40 MHz的采樣率和18位的A/D轉(zhuǎn)換率, 其頻帶范圍為1 kHz~3 MHz。AE數(shù)據(jù)由AE win for PCI–2軟件進(jìn)行采集、存儲和分析。在試驗過程中, AE監(jiān)測系統(tǒng)的前置放大器增益設(shè)置為40 dB,閾值設(shè)定為45 dB, 采樣率為1 MHz, 并采用雙通道傳感器進(jìn)行數(shù)據(jù)采集, 諧振頻率設(shè)定為20~400 kHz。使用雙通道傳感器采集數(shù)據(jù), 傳感器沿試樣兩側(cè)軸向中部布置。試驗時控制室內(nèi)溫度為26 ℃, 濕度為58RH%, 試驗中監(jiān)測應(yīng)力、應(yīng)變、AE計數(shù)等變量。循環(huán)加載模式為載荷控制與位移控制相結(jié)合的多次循環(huán)加載。前4次循環(huán)采用載荷控制, 末循環(huán)采用位移控制。其中前2次循環(huán)主要用于AE法分析, 之后的循環(huán)用于DRA法分析。首循環(huán)峰值載荷設(shè)定為σhmax的80%。之后3次循環(huán)既要減少首循環(huán)加載摩擦型AE的干擾, 又要避免巖石擴(kuò)容損傷產(chǎn)生的影響, 其循環(huán)峰值載荷應(yīng)大于σhmax, 且不超過巖石抗壓強(qiáng)度σc的63%~70%。試驗加載速率設(shè)置為0.05 kN/s, 卸載速率為0.25 kN/s。加載方式如圖2所示。

圖2 循環(huán)加載模式示意Fig.2 Schematic diagram of cyclic loading mode

采用AE–DRA法[23]識別Kaiser點。在巖石試件進(jìn)行單軸壓縮試驗時, 記錄聲發(fā)射事件的數(shù)量、能量和頻率等參數(shù), 然后利用數(shù)據(jù)重構(gòu)算法(DRA)對聲發(fā)射信號進(jìn)行分析, 找出聲發(fā)射響應(yīng)曲線的最大曲率點, 并將其作為Kaiser點。進(jìn)一步根據(jù)Kaiser點對應(yīng)的應(yīng)力值計算地應(yīng)力。

1.3 古地磁巖芯定向

由于鉆取巖芯時采用非定向取芯方法, 故采用古地磁巖芯定向法確定巖芯在地下的原始方位。其基本原理為: 巖芯在形成時會記錄當(dāng)時的地磁場方向, 形成穩(wěn)定的剩磁。通過實驗室的高靈敏度磁力儀系統(tǒng), 分離和測定巖芯的剩磁方向, 進(jìn)而確定巖芯在地下的空間狀態(tài)[24]。

試驗使用美國2G超導(dǎo)磁力儀和英國MMDT80熱退磁儀。在退磁后進(jìn)行剩磁測量, 獲得巖芯黏滯剩磁的平均磁偏角Da與磁傾角Ia, 結(jié)合標(biāo)志線與古地磁方向, 恢復(fù)巖芯在地下的原始方位; 進(jìn)一步計算出最大水平主應(yīng)力方向。

布爾臺礦區(qū)部分鉆孔巖芯Kaiser法測定的地應(yīng)力值及方向詳見表1。

表1 布爾臺礦區(qū)部分鉆孔巖芯Kaiser法地應(yīng)力測試結(jié)果Table 1 Stress measurement results in Buertai mine area based on Kaiser method

1.4 地應(yīng)力測量結(jié)果驗證

由表1可知, 同一鉆孔中的地應(yīng)力測點, 隨著深度的增加, 最大水平主應(yīng)力σH、最小水平主應(yīng)力σh和豎直主應(yīng)力σv均逐漸增加。使用最小二乘法線性擬合各點地應(yīng)力數(shù)據(jù), 結(jié)果如圖3所示?？梢园l(fā)現(xiàn): 豎直主應(yīng)力擬合曲線斜率最大, 即隨埋深增加豎直方向主應(yīng)力增加最快; 埋深127 m以淺的地層,地應(yīng)力呈σv<σh<σH的逆斷型應(yīng)力狀態(tài); 埋深1 725 m以深的地層, 地應(yīng)力呈σh<σH<σv的正斷型應(yīng)力狀態(tài);埋深處于127~1 725 m的地層, 地應(yīng)力呈σh<σv<σH的走滑應(yīng)力狀態(tài)。該結(jié)果與其他針對神東礦區(qū)構(gòu)造應(yīng)力場的研究結(jié)果相一致[25–26]。需要說明: 采用線性函數(shù)對水平方向地應(yīng)力數(shù)據(jù)點進(jìn)行擬合時, 離散性較大。因此, 嘗試將巖石彈性模量E引入對地應(yīng)力的線性回歸擬合中, 進(jìn)行多元線性回歸。多元線性回歸擬合結(jié)果如圖4所示。

圖3 布爾臺礦區(qū)地應(yīng)力隨埋深變化規(guī)律Fig.3 In-situ stress variation law of the Buertai mine area with the change of depth

圖4 布爾臺礦區(qū)地應(yīng)力隨埋深與彈性模量變化規(guī)律Fig.4 In-situ stress variation law of the Buertai mine area with the change of depth and modulus of elasticity

由圖4可知, 引入彈性模量E后, 擬合優(yōu)度R2明顯上升。對最大水平主應(yīng)力的多元線性回歸R2為0.612 01, 較線性回歸的0.501 41提升22.06%; 最小水平主應(yīng)力的多元線性回歸R2為0.631 05, 較線性回歸的0.506 27提升24.65%。結(jié)果表明, 引入彈性模量E對于改善水平方向地應(yīng)力擬合效果具有顯著作用。同時也說明水平方向地應(yīng)力與彈性模量間具有一定關(guān)聯(lián)性。

圖5分析了側(cè)壓系數(shù), 即(σH+σh)/2與σv的比值k,隨埋深變化的規(guī)律。能夠看出, 隨著埋深的不斷增加, 側(cè)壓系數(shù)k的分布范圍逐漸收窄。當(dāng)埋深大于300 m時,k值收斂于1附近。k值隨埋深的變化規(guī)律為k=88/H+0.72, 且35/H+0.29

圖5 布爾臺礦區(qū)側(cè)壓系數(shù)隨埋深變化規(guī)律Fig.5 Lateral pressure coefficient variation law of the Buertai mine area with the change of depth

由于鉆孔位置集中在布爾臺礦三、四盤區(qū), 為驗證地應(yīng)力測量結(jié)果在全礦井的適用性, 選取布爾臺礦一、二盤區(qū)3個鉆孔的空心包體法地應(yīng)力測試結(jié)果[27]與三、四盤區(qū)埋深相近的Kaiser法效應(yīng)測點測量得到的3組地應(yīng)力值進(jìn)行對比。對比結(jié)果見表2, 表2中BK209–21, BK209–22與BK209–23為位于三盤區(qū)的Kaiser法地應(yīng)力測點。BK1號, BK4號與BK5號為位于一、二盤區(qū)的空心包體法地應(yīng)力測點。

表2 一、二盤區(qū)與三、四盤區(qū)地應(yīng)力測試結(jié)果對比Table 2 Comparison of in-situ stress test results between the first and second plates, and those from the third and fourth plates

由表2可知, 一、二盤區(qū)σH測量值與三、四盤區(qū)σH測量值間相對誤差平均值為14.90%,σh相對誤差平均值為19.67%,σv相對誤差平均值為15.47%。因此, 一、二盤區(qū)地應(yīng)力測試結(jié)果與三、四盤區(qū)地應(yīng)力測試結(jié)果誤差較小, 同時考慮到測試方法的不同, 誤差在工程合理范圍內(nèi), 能夠證明Kaiser法地應(yīng)力測量結(jié)果在全礦井的適用性。

2 基于改進(jìn)組合彈簧模型的地應(yīng)力計算方法

盡管經(jīng)典組合彈簧模型考慮了巖層物理力學(xué)參數(shù)對地應(yīng)力的影響, 能較為準(zhǔn)確反映地應(yīng)力場的分布規(guī)律。但該方法假設(shè)地層在水平方向最大和最小主應(yīng)變?yōu)楹愣ㄖ?。因礦區(qū)地域范圍廣, 地質(zhì)構(gòu)造復(fù)雜, 地應(yīng)力場受巖體成因、巖體結(jié)構(gòu)和地質(zhì)構(gòu)造運動等多因素的影響, 表現(xiàn)出顯著的非均質(zhì)性。故組合彈簧模型無法準(zhǔn)確描述區(qū)域地應(yīng)力場分布。

為提高地應(yīng)力場反演計算精度, 將克里金插值方法引入到經(jīng)典組合彈簧模型的應(yīng)變計算過程中,提出一種基于Kaiser法實測地應(yīng)力數(shù)據(jù)與改進(jìn)組合彈簧模型的地應(yīng)力場計算方法。該方法在考慮巖層物性的同時克服了組合彈簧模型理論中應(yīng)變參量無法反映圍巖變形非均質(zhì)性的問題。

2.1 地應(yīng)力測點處應(yīng)變計算

基于Kaiser法測定的地應(yīng)力數(shù)據(jù), 使用組合彈簧模型, 結(jié)合實測地應(yīng)力值與巖石物理力學(xué)參數(shù)求取各地應(yīng)力測點處的應(yīng)變。最大水平主應(yīng)變與最小水平主應(yīng)變的計算公式為

式中,Hε為該點最大水平主應(yīng)變;hε為該點最小水平主應(yīng)變;ν為該點巖石泊松比;Hσ為該點最大水平主應(yīng)力, GPa;hσ為該點最小水平主應(yīng)力, GPa;vσ為該點豎直方向主應(yīng)力, GPa;E為該點巖石楊氏模量, GPa。

2.2 研究區(qū)域任意點應(yīng)變插值計算

按測定的最大水平主應(yīng)力方向?qū)?yīng)變分解為x方向上的正應(yīng)變εx,y方向上的正應(yīng)變εy及剪應(yīng)變γxy?；诘貞?yīng)力測點處的應(yīng)變, 使用克里金插值方法求取應(yīng)變場空間分布?？死锝鸩逯捣煽紤]空間數(shù)據(jù)的各向異性、非線性和多變量等特征, 其計算式如下:

使用克里金插值法得到研究區(qū)域內(nèi)任一點x方向上的正應(yīng)變εx,y方向上的正應(yīng)變εy及剪應(yīng)變γxy。

2.3 研究區(qū)域內(nèi)任意點地應(yīng)力計算

根據(jù)研究區(qū)域內(nèi)任一點的應(yīng)變張量, 進(jìn)而可計算出該點最大水平主應(yīng)變與最小水平主應(yīng)變及方向, 計算公式為

式中,θ為最大主應(yīng)變方向;xε為x方向上的正應(yīng)變;yε為y方向上的正應(yīng)變;γxy為xy方向上的剪應(yīng)變。

在得到任一點的水平方向最大和最小主應(yīng)變大小和方向后, 可結(jié)合該點的楊氏模量E(GPa)、泊松比ν等巖石物理力學(xué)性質(zhì), 根據(jù)組合彈簧模型計算出該點的地應(yīng)力, 具體算公式為

3 結(jié)果及分析

3.1 地應(yīng)力反演結(jié)果分析

為驗證改進(jìn)的組合彈簧模型地應(yīng)力場反演方法的有效性, 筆者將其與側(cè)壓系數(shù)法[28]、線性擬合法[29]、分段擬合法[30]及經(jīng)典組合彈簧模型方法[31]進(jìn)行了對比, 圖6為地應(yīng)力反演方法對比評估流程。

圖6 地應(yīng)力反演方法對比評估流程Fig.6 Comparative evaluation process of in-situ stress inversion methods

基于經(jīng)古地磁巖芯定向的15個地應(yīng)力測點數(shù)據(jù), 使用各方法進(jìn)行地應(yīng)力反演, 隨后測試各方法對71個非定向地應(yīng)力測點的反演計算精度, 從而對比各擬合方法的計算準(zhǔn)確性。

3.1.1 側(cè)壓系數(shù)擬合結(jié)果

側(cè)壓系數(shù)法是一種簡單的地應(yīng)力反演方法, 該方法假設(shè)巖石在加載過程中保持彈性, 水平方向地應(yīng)力大小與豎直方向地應(yīng)力大小成正比, 比例系數(shù)為側(cè)壓系數(shù)λ。側(cè)壓系數(shù)計算公式為

式中,n為地應(yīng)力測點數(shù)量;σHi為第i個測點的最大水平主應(yīng)力, GPa;σhi為第i個測點的最小水平主應(yīng)力, GPa;σvi為第i個測點的豎直方向主應(yīng)力, GPa。

基于實測地應(yīng)力計算得到布爾臺礦側(cè)壓系數(shù)為1.034。側(cè)壓系數(shù)擬合得到研究區(qū)域任意點的水平方向地應(yīng)力為

式中,H為埋深, m;ρ為覆巖密度, kg/m3;g為重力加速度, 取9.8 m/s2。

3.1.2 線性函數(shù)擬合結(jié)果

線性擬合法是一種常用的地應(yīng)力反演方法, 它利用最小二乘法對實測數(shù)據(jù)進(jìn)行線性擬合, 得到地應(yīng)力的空間分布, 其擬合結(jié)果如圖7所示。

圖7 布爾臺礦地應(yīng)力線性擬合結(jié)果Fig.7 Linear fitting results of in-situ stress in Buertai Mine

使用線性擬合方法得到的地應(yīng)力擬合公式為

3.1.3 分段函數(shù)擬合結(jié)果

分段函數(shù)擬合法對地應(yīng)力數(shù)據(jù)進(jìn)行分段回歸,每100 m埋深設(shè)定為一個分段, 回歸表達(dá)式為

式中,β0,…,β6為函數(shù)系數(shù);xi1,…,xi5為虛擬變量。

當(dāng)H<100 m時,xi1,…,xi5均為0; 當(dāng)100

圖8 布爾臺礦地應(yīng)力分段擬合結(jié)果Fig.8 Piecewise fitting results of the in-situ stress in Buertai Mine

使用分段擬合方法得到的擬合式為

3.1.4 經(jīng)典組合彈簧模型反演結(jié)果

經(jīng)典組合彈簧模型地應(yīng)力反演方法通過求取最佳水平方向地應(yīng)力系數(shù)確定地應(yīng)力分布, 基于地應(yīng)力實測數(shù)據(jù), 使用經(jīng)典組合彈簧模型求得的平均最大主應(yīng)變與平均最小主應(yīng)變分別為1.04×10–3與4.56×10–4。使用經(jīng)典組合彈簧模型方法得到的地應(yīng)力擬合式為

3.1.5 改進(jìn)組合彈簧模型反演結(jié)果

使用筆者所提出的改進(jìn)組合彈簧模型得出的研究區(qū)域地應(yīng)力各分量隨埋深變化曲線如圖9所示。

圖9 布爾臺礦地應(yīng)力改進(jìn)組合彈簧模型反演結(jié)果Fig.9 Inversion results of the improved composite spring model for the in-situ stress in Buertai Mine

3.2 地應(yīng)力反演效果比較

使用均方根誤差(RMSE)與平均絕對百分比誤差(SMAPE)指標(biāo)對各類地應(yīng)力反演方法得到的σH與σh的絕對誤差與相對誤差進(jìn)行評估。RMSE與SMAPE的計算方法為

式中,n為樣本個數(shù);yi為第i個樣本的實際值;為第i個樣本的預(yù)測值。

表3給出了各反演計算模型對比分析部分結(jié)果。

表3 各反演方法計算誤差Table 3 Errors of different inversion methods.

由表3可知, 改進(jìn)組合彈簧模型反演得到的Hσ與實測值間的均方根誤差為1.849 8, 較側(cè)壓系數(shù)法的3.153 2降低41.3%; 較線性擬合法的3.272 8降低43.5%; 較分段擬合法的3.066 3降低39.7%, 較經(jīng)典組合彈簧模型反演的2.305 1降低19.8%。實測數(shù)據(jù)反演法得到的hσ與實測值間的均方根誤差為1.811 9, 較側(cè)壓系數(shù)法的3.183 8降低43.1%; 較線性擬合法的2.607 4降低30.5%; 較分段擬合法的2.448 5降低26.0%, 較經(jīng)典組合彈簧模型反演的1.981 5降低8.6%。

此外在評價反演結(jié)果與實測數(shù)據(jù)間的相對誤差時, 改進(jìn)組合彈簧模型反演得到的Hσ與實測值間的SMAPE值為17.27%, 較側(cè)壓系數(shù)法的26.86%降低35.7%, 較線性擬合法的23.52%降低26.6%, 較分段擬合法的25.12%降低31.3%, 較經(jīng)典組合彈簧模型的18.14%降低4.8%。改進(jìn)組合彈簧模型反演得到的hσ與實測值間的SMAPE值為21.62%, 較側(cè)壓系數(shù)法的34.77%降低37.8%, 較線性擬合法的31.14%降低30.6%, 較分段擬合法的28.76%降低24.8%, 較經(jīng)典組合彈簧模型的22.14%降低2.3%。

綜上, 改進(jìn)組合彈簧模型計算的最大水平主應(yīng)力及最小水平主應(yīng)力結(jié)果與實測值間的均方根誤差為1.830 9, 相較其他方法減少14.6%~42.2%, 平均絕對百分比誤差為19.45%，相較其他減少3.4%~36.9%。所提出的改進(jìn)組合彈簧模型方法計算誤差最小, RMSE和SMAPE值均最小, 說明其擬合效果和穩(wěn)定性最好, 能有效反映水平地應(yīng)力與應(yīng)變的非線性關(guān)系。除改進(jìn)組合彈簧模型外, 經(jīng)典組合彈簧模型的RMSE與SMAPE最小。相比之下, 側(cè)壓系數(shù)法和線性擬合法的RMSE與SMAPE較大, 說明其擬合效果與穩(wěn)定性較差。

3.3 數(shù)據(jù)敏感性分析

為分析數(shù)據(jù)量對改進(jìn)組合彈簧模型擬合精度的影響, 進(jìn)行了改進(jìn)組合彈簧模型的數(shù)據(jù)敏感性分析。由于克里金插值需要最少3個數(shù)據(jù)點, 故分析了數(shù)據(jù)點數(shù)量從3增加到15的過程中擬合精度的變化。結(jié)果顯示, 隨著數(shù)據(jù)點的增多, 擬合精度逐漸上升。其中, 最大水平主應(yīng)力的RMSE值從4.88下降至2.56, SMAPE值從26.91%下降至25.83%; 最小水平主應(yīng)力的RMSE值從6.13下降至1.51, SMAPE值從49.75%下降至18.93%。同時, 增加同樣數(shù)量數(shù)據(jù)點對擬合精度的提升逐漸下降。在數(shù)據(jù)量增加至12~15時, 各擬合精度指標(biāo)均不再發(fā)生顯著下降。因此繼續(xù)增加數(shù)據(jù)點對精度提升作用有限, 故15個數(shù)據(jù)點能夠較好地擬合該區(qū)域地應(yīng)力。各精度指標(biāo)隨數(shù)據(jù)量變化的曲線如圖10所示。

圖10 改進(jìn)組合彈簧模型擬合誤差隨數(shù)據(jù)量變化曲線Fig.10 Error curves of the improved combined spring model varies with the amount of fitted data

3.4 模型應(yīng)用效果評價

為進(jìn)一步驗證所提出方法的有效性和優(yōu)越性,將地應(yīng)力反演計算方法融入FLAC3D軟件建模。分別以改進(jìn)組合彈簧模型方法與側(cè)壓系數(shù)方法計算得到的結(jié)果作為數(shù)值模型應(yīng)力邊界條件進(jìn)行模擬,并與實測礦壓數(shù)據(jù)進(jìn)行對比。從而評估不同地應(yīng)力計算方法對模擬準(zhǔn)確性的影響。

使用FLAC3D軟件模擬布爾臺礦42105工作面的開采過程, 在45號及85號支架位置處設(shè)置支架阻力測線。

布爾臺礦42105工作面布置如圖11所示。

圖11 布爾臺礦42105工作面平面布置Fig.11 Layout of No.42105 working face in Buertai Mine

3.4.1 模型建立及巖石力學(xué)參數(shù)賦值

布爾臺礦42105工作面位于4–2煤一盤區(qū), 走向長5 231 m, 傾向長230 m。基于SJ–3–1鉆孔與SJ–3–2鉆孔信息構(gòu)建布爾臺42105工作面數(shù)值模型, 鉆孔位置如圖11所示。數(shù)值模型長900 m, 寬400 m, 高487.62 m, 共包含2 380 800個單元, 2 440 125個網(wǎng)格點。根據(jù)布爾臺礦42105工作面鉆孔巖芯的物理力學(xué)性質(zhì)測試數(shù)據(jù)及以往文獻(xiàn)[32]中的神東礦區(qū)巖石物理力學(xué)性質(zhì), 確定了數(shù)值模型的巖石物理力學(xué)參數(shù), 具體見表4。

表4 布爾臺礦42105工作面巖石物理力學(xué)參數(shù)Table 4 Physical and mechanical parameters of No.42105 working face in Buertai Mine

3.4.2 地應(yīng)力分塊精細(xì)化加載

為將地應(yīng)力精確加載到數(shù)值模型上, 將工作面的側(cè)面應(yīng)力邊界劃分為4 312個矩形區(qū)域, 劃分情況如圖12所示。

圖12 數(shù)值模型應(yīng)力精細(xì)加載方法示意Fig.12 Diagram of numerical model based on the stress fine loading method

在FLAC3D中, 應(yīng)力邊界條件可表示為

式中,Cσ為坐標(biāo)點(0,0)處的應(yīng)力值; ?σ/?x為應(yīng)力沿x方向的梯度; ?σ/?y為應(yīng)力沿y方向的梯度。

任意矩形區(qū)域的應(yīng)力邊界如圖13所示。

圖13 FLAC3D應(yīng)力邊界條件示意Fig.13 Diagram of the stress boundary condition in the model of FLAC3D

已知矩形區(qū)域4個邊界點處的應(yīng)力, 求取該區(qū)域的應(yīng)力邊界條件等價于求解以下矩陣:

式中, (xi,yi)為邊界點Pi的坐標(biāo);σi為邊界點Pi處的應(yīng)力。

對于任意一個矩形區(qū)域, 其應(yīng)力邊界條件可以由3條FLAC3D指令確定, 分別用于指定坐標(biāo)原點處σx,σy與τxy的值及其沿x,y,z方向的梯度。使用Python編寫邊界條件求解程序, 運行后共生成12 936條邊界條件加載指令。最后使用FLAC3D軟件逐條運行邊界加載指令, 以此實現(xiàn)應(yīng)力邊界條件的精細(xì)化加載。

3.4.3 工作面支架支撐力分布

采用FLAC3D模擬布爾臺礦42105工作面開采過程, 從開切眼位置開始每次向工作面推進(jìn)方向開挖3.5 m, 共開采70 m。對鄰近工作面的一列頂板單元施加支撐力從而體現(xiàn)液壓支架對頂板的支撐作用。FLAC3D模型每運算50步后, 根據(jù)頂板下沉量更新支撐力, 直到模型最大不平衡力比率小于10–5。記錄每一步開挖運算至平衡時支撐力的值。為便于與現(xiàn)場實測數(shù)據(jù)比較, 根據(jù)布爾臺礦42105工作面使用的ZFY18000/25/39D型液壓支架參數(shù)[33]將支撐力換算為液壓支架立柱壓強(qiáng)。

隨后將基于改進(jìn)組合彈簧模型法、經(jīng)典組合彈簧模型法與側(cè)壓系數(shù)法計算得到的3種應(yīng)力邊界條件下的數(shù)值模擬結(jié)果與現(xiàn)場實測數(shù)據(jù)進(jìn)行對比,結(jié)果如圖14所示。

圖14 支架支柱壓強(qiáng)實測值與不同應(yīng)力邊界條件模擬結(jié)果Fig.14 Measured values of support pressure and the simulated results of different stress boundary conditions

由圖14可知, 使用側(cè)壓系數(shù)法得到的應(yīng)力場作為數(shù)值模型的水平應(yīng)力邊界條件時, 最大支架阻力在開采至35 m處出現(xiàn), 來壓時45號支架立柱壓強(qiáng)為40.8 MPa, 85號支架立柱壓強(qiáng)為40.6 MPa。使用經(jīng)典組合彈簧模型方法得到的應(yīng)力場作為數(shù)值模型的水平應(yīng)力邊界條件時, 最大支架阻力在開采至35 m處出現(xiàn), 來壓時45號支架立柱壓強(qiáng)為64.7 MPa,85號支架立柱壓強(qiáng)為63.0 MPa。使用改進(jìn)組合彈簧模型法得到的應(yīng)力場作為數(shù)值模型的水平應(yīng)力邊界條件時, 最大支架阻力出現(xiàn)在開采至56 m處時,來壓時45號液壓支架立柱壓強(qiáng)為47.8 MPa, 85號支架立柱壓強(qiáng)為46.6 MPa。根據(jù)布爾臺礦42105工作面液壓支架實測數(shù)據(jù)[34], 首次基本頂破斷引起的礦壓顯現(xiàn)發(fā)生于工作面開采到52.0~61.5 m時, 來壓期間45號支架峰值立柱壓強(qiáng)為49.9 MPa, 85號支架峰值立柱壓強(qiáng)約為48.0 MPa。因此, 相較于使用側(cè)壓系數(shù)計算結(jié)果作為數(shù)值模型的水平應(yīng)力邊界條件,利用改進(jìn)組合彈簧模型的計算結(jié)果能更準(zhǔn)確地模擬現(xiàn)場實際情況。

采用均方根誤差(RMSE)與平均絕對百分比誤差(SMAPE)指標(biāo)對使用3種應(yīng)力場計算結(jié)果作為應(yīng)力邊界條件的模擬結(jié)果與實測礦壓數(shù)據(jù)間的差異進(jìn)行定量評估, 指標(biāo)定義見式(14), 評估結(jié)果見表5。

表5 各計算方法數(shù)值模擬結(jié)果與實測數(shù)據(jù)間誤差Table 5 Error between numerical simulation results and measured data

由表5可以發(fā)現(xiàn), 使用改進(jìn)組合彈簧模型法得到的地應(yīng)力場作為模型水平方向應(yīng)力邊界時, 礦壓模擬值與實測值間的誤差較小。

對于45號支架的監(jiān)測數(shù)據(jù)而言, 改進(jìn)組合彈簧模型法與實測值間RMSE為68.10, 較側(cè)壓系數(shù)法的75.53低9.84%, 較經(jīng)典組合彈簧模型的103.16低33.99%。SMAPE為16.95%, 較側(cè)壓系數(shù)法的18.37%低7.73%, 較經(jīng)典組合彈簧模型的20.90%低18.90%。此外, 對于85號支架的監(jiān)測數(shù)據(jù)而言, 改進(jìn)組合彈簧模型法與實測值間RMSE為65.00, 較側(cè)壓系數(shù)法的68.76低5.47%, 較經(jīng)典組合彈簧模型的93.89低30.77%。SMAPE為14.54%, 較側(cè)壓系數(shù)法的15.35%低5.28%, 較經(jīng)典組合彈簧模型的18.81%低22.70%。

綜上, 模擬得到的工作面礦壓數(shù)據(jù)與實測數(shù)據(jù)間均方根誤差平均值為66.55, 較使用側(cè)壓系數(shù)法降低7.75%, 較使用經(jīng)典組合彈簧模型法降低32.45%。平均絕對百分比誤差平均值為15.75%, 較使用側(cè)壓系數(shù)法降低6.61%, 較使用經(jīng)典組合彈簧模型法降低20.67%。通過對比以改進(jìn)組合彈簧模型與側(cè)壓系數(shù)法的計算結(jié)果作為數(shù)值模擬邊界條件的模擬結(jié)果與實測值間的誤差, 證明了所提出的方法較其他模型能更準(zhǔn)確地模擬現(xiàn)場實際情況。從而為使用數(shù)值模擬方法進(jìn)行工作面礦壓顯現(xiàn)分析, 強(qiáng)礦壓防治以及巷道優(yōu)化布置等提供一種新的數(shù)值模型應(yīng)力邊界條件確定方法, 具有重要的理論指導(dǎo)意義和工程應(yīng)用價值。

4 結(jié) 論

(1)以布爾臺礦為研究對象, 基于Kaiser法測量地應(yīng)力數(shù)據(jù)。測量結(jié)果顯示: 同一鉆孔中的地應(yīng)力測點, 隨深度的增加, 最大水平主應(yīng)力σH、最小水平主應(yīng)力σh和豎直主應(yīng)力σv均逐漸增加。隨著埋深不斷增加, 側(cè)壓系數(shù)k的分布范圍逐漸變小。當(dāng)埋深大于300 m時,k值收斂于1附近。k值隨埋深的變化規(guī)律為k=88/H+0.72, 且35/H+0.29

(2)將克里金插值方法引入經(jīng)典組合彈簧模型,提出了一種同時考慮巖層物性與圍巖變形非均質(zhì)性的改進(jìn)組合彈簧模型。通過比較改進(jìn)組合彈簧模型、側(cè)壓系數(shù)法、線性擬合法、分段函數(shù)擬合法和經(jīng)典組合彈簧模型與實測數(shù)據(jù)間的誤差, 評估了各方法的計算精度。結(jié)果表明改進(jìn)組合彈簧模型計算的最大水平主應(yīng)力及最小水平主應(yīng)力結(jié)果與實測值間的均方根誤差為1.8309, 相較其他方法減少14.6%~42.2%, 平均絕對百分比誤差為19.45%,相較其他減少3.4%~36.9%。

(3)以布爾臺礦42015工作面為工程背景, 使用改進(jìn)組合彈簧模型方法計算得到的結(jié)果作為FLAC3D數(shù)值模型應(yīng)力邊界條件進(jìn)行模擬。模擬得到的工作面礦壓數(shù)據(jù)與實測數(shù)據(jù)間均方根誤差平均值為66.55, 較使用側(cè)壓系數(shù)法降低7.75%, 較使用經(jīng)典組合彈簧模型法降低32.45%。平均絕對百分比誤差平均值為15.75%, 較使用側(cè)壓系數(shù)法降低6.61%, 較使用經(jīng)典組合彈簧模型法降低20.67%。結(jié)果表明采用改進(jìn)組合彈簧模型的計算結(jié)果作為數(shù)值模型邊界條件較其他模型能更準(zhǔn)確地模擬現(xiàn)場實際情況。