基于改進組合彈簧模型的礦井地應力場計算方法

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

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

地應力場是影響巖層運移和原位巖體力學性質的重要因素, 對采礦、隧道開挖、油氣開采等影響顯著。地應力是采礦工程中覆巖運動的力源, 在構造應力強烈的地區進行采礦活動, 可能會引起沖擊地壓、巖爆等災害, 影響礦山的安全生產。因此準確獲取地應力場對于提高開采設計的科學性及保障礦山生產安全有著重要意義。

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

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

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

1 研究區域與地應力測量

1.1 研究區域與取樣

布爾臺礦是神東煤炭集團主力生產礦井之一,位于內蒙古自治區鄂爾多斯市伊金霍洛旗烏蘭木倫鎮。大地構造分區屬于華北地臺鄂爾多斯臺向斜伊盟隆起中東部, 伊陜斜坡東北部。基本構造形態為一南西傾斜的近水平產狀的單斜構造, 井田構造簡單, 層位穩定。

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

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

1.2 基于Kaiser效應的地應力測量

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

圖2 循環加載模式示意Fig.2 Schematic diagram of cyclic loading mode

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

1.3 古地磁巖芯定向

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

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

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

表1 布爾臺礦區部分鉆孔巖芯Kaiser法地應力測試結果Table 1 Stress measurement results in Buertai mine area based on Kaiser method

1.4 地應力測量結果驗證

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

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

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

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

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

圖5 布爾臺礦區側壓系數隨埋深變化規律Fig.5 Lateral pressure coefficient variation law of the Buertai mine area with the change of depth

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

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

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

2 基于改進組合彈簧模型的地應力計算方法

盡管經典組合彈簧模型考慮了巖層物理力學參數對地應力的影響, 能較為準確反映地應力場的分布規律。但該方法假設地層在水平方向最大和最小主應變為恒定值。因礦區地域范圍廣, 地質構造復雜, 地應力場受巖體成因、巖體結構和地質構造運動等多因素的影響, 表現出顯著的非均質性。故組合彈簧模型無法準確描述區域地應力場分布。

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

2.1 地應力測點處應變計算

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

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

2.2 研究區域任意點應變插值計算

按測定的最大水平主應力方向將應變分解為x方向上的正應變εx,y方向上的正應變εy及剪應變γxy。基于地應力測點處的應變, 使用克里金插值方法求取應變場空間分布。克里金插值法可考慮空間數據的各向異性、非線性和多變量等特征, 其計算式如下:

使用克里金插值法得到研究區域內任一點x方向上的正應變εx,y方向上的正應變εy及剪應變γxy。

2.3 研究區域內任意點地應力計算

根據研究區域內任一點的應變張量, 進而可計算出該點最大水平主應變與最小水平主應變及方向, 計算公式為

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

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

3 結果及分析

3.1 地應力反演結果分析

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

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

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

3.1.1 側壓系數擬合結果

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

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

基于實測地應力計算得到布爾臺礦側壓系數為1.034。側壓系數擬合得到研究區域任意點的水平方向地應力為

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

3.1.2 線性函數擬合結果

線性擬合法是一種常用的地應力反演方法, 它利用最小二乘法對實測數據進行線性擬合, 得到地應力的空間分布, 其擬合結果如圖7所示。

圖7 布爾臺礦地應力線性擬合結果Fig.7 Linear fitting results of in-situ stress in Buertai Mine

使用線性擬合方法得到的地應力擬合公式為

3.1.3 分段函數擬合結果

分段函數擬合法對地應力數據進行分段回歸,每100 m埋深設定為一個分段, 回歸表達式為

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

當H<100 m時,xi1,…,xi5均為0; 當100

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

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

3.1.4 經典組合彈簧模型反演結果

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

3.1.5 改進組合彈簧模型反演結果

使用筆者所提出的改進組合彈簧模型得出的研究區域地應力各分量隨埋深變化曲線如圖9所示。

圖9 布爾臺礦地應力改進組合彈簧模型反演結果Fig.9 Inversion results of the improved composite spring model for the in-situ stress in Buertai Mine

3.2 地應力反演效果比較

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

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

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

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

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

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

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

3.3 數據敏感性分析

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

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

3.4 模型應用效果評價

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

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

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

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

3.4.1 模型建立及巖石力學參數賦值

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

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

3.4.2 地應力分塊精細化加載

為將地應力精確加載到數值模型上, 將工作面的側面應力邊界劃分為4 312個矩形區域, 劃分情況如圖12所示。

圖12 數值模型應力精細加載方法示意Fig.12 Diagram of numerical model based on the stress fine loading method

在FLAC3D中, 應力邊界條件可表示為

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

任意矩形區域的應力邊界如圖13所示。

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

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

式中, (xi,yi)為邊界點Pi的坐標;σi為邊界點Pi處的應力。

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

3.4.3 工作面支架支撐力分布

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

隨后將基于改進組合彈簧模型法、經典組合彈簧模型法與側壓系數法計算得到的3種應力邊界條件下的數值模擬結果與現場實測數據進行對比,結果如圖14所示。

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

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

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

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

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

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

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

4 結 論

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

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

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