郭延輝，侯克鵬

(1.云南農(nóng)業(yè)大學(xué)建筑工程學(xué)院，昆明 650201；2.昆明理工大學(xué)國土資源工程學(xué)院，昆明 650093)

地應(yīng)力是地質(zhì)構(gòu)造和自重共同作用在地質(zhì)體內(nèi)形成的原始應(yīng)力，初始地應(yīng)力的確定是工程分析與設(shè)計的基礎(chǔ)，對于工程的穩(wěn)定與安全至關(guān)重要[1]。礦區(qū)巖體地應(yīng)力測量是提供地應(yīng)力場的一種最為直接的方法，然而面對廣域的工程廠址，只能開展有限測點(diǎn)的地應(yīng)力測量，而這些數(shù)據(jù)只能反映測量點(diǎn)處局部地應(yīng)力場的情況，且有時會由于某些原因致使測量數(shù)據(jù)存在一定誤差[2-3]。為了能夠?qū)こ淘O(shè)計和施工提供更可靠的地應(yīng)力數(shù)據(jù)，根據(jù)現(xiàn)場的地質(zhì)構(gòu)造特征，結(jié)合有限的地應(yīng)力實測數(shù)據(jù)，對初始應(yīng)力場進(jìn)行反演計算和分析，獲得更為準(zhǔn)確、適用范圍較大的地應(yīng)力場十分必要[4]。

近年來，研究人員提出了多種地應(yīng)力場反演方法。李飛等[5]提出批數(shù)據(jù)處理(GMDH)神經(jīng)網(wǎng)絡(luò)算法，結(jié)合地應(yīng)力場分布隨埋深的非線性特征及局部地質(zhì)構(gòu)造處地應(yīng)力場的非連續(xù)性特點(diǎn)，構(gòu)建出形成復(fù)雜地質(zhì)體地應(yīng)力場的邊界條件模式，并基于現(xiàn)場的稀少樣本測點(diǎn)數(shù)據(jù)對杏山鐵礦地應(yīng)力場進(jìn)行反演和重構(gòu)；汪吉林等[6]提出以構(gòu)造指標(biāo)優(yōu)化人工神經(jīng)網(wǎng)絡(luò)(ANN)模型的建模方案，并開展了基于構(gòu)造控制的地應(yīng)力人工神經(jīng)網(wǎng)絡(luò)反演方法研究；易達(dá)等[7]研究了巖體初始應(yīng)力場的遺傳算法與有限元聯(lián)合反演法；戚藍(lán)等[8]用系統(tǒng)工程灰色理論進(jìn)行了地應(yīng)力分析，在處理實測資料少的小樣本問題取得了較好的效果；李欣等[9]采用改進(jìn)的粒子群優(yōu)化算法，對地應(yīng)力場進(jìn)行了反演計算，結(jié)果表明改進(jìn)后的算法計算精度和效率更高；上述這些非線性反演方法雖然能夠反演分析復(fù)雜區(qū)域地應(yīng)力場，卻存在收斂速度較慢、參數(shù)選擇過于依賴經(jīng)驗等問題，并可能導(dǎo)致同樣的研究對象而反演結(jié)果卻并不唯一的現(xiàn)象。數(shù)值模擬方法以其三維化、可視化、成本低和能夠計算復(fù)雜問題等優(yōu)點(diǎn)廣泛應(yīng)用于深部地應(yīng)力場的模擬計算[10-12]。然而，目前關(guān)于地應(yīng)力場的反演研究大多是以自重應(yīng)力場為主，針對初始高水平構(gòu)造應(yīng)力場的反演方法方面的研究仍然不足[13-14]。為此，以獅子山礦區(qū)為研究背景，在模擬邊界討論的基礎(chǔ)上，采用FLAC3D反演分析各種邊界條件下計算域內(nèi)水平構(gòu)造應(yīng)力與垂直應(yīng)力的分布特征，最終采用基于初始應(yīng)變能理論的方法，對獅子山礦區(qū)初始高水平構(gòu)造地應(yīng)力場進(jìn)行反演及驗證，期望能對其他類似工程有一定的借鑒作用。

1 關(guān)于地應(yīng)力反演方法的討論

通常數(shù)值計算過程中賦地應(yīng)力的方法主要有：①模型內(nèi)部部分或某些點(diǎn)的地應(yīng)力實測值已知，采用位移邊界條件，以一定的特征和規(guī)律將已知點(diǎn)應(yīng)力差值分布到模型中，最后對應(yīng)力場反演分析并生成最終應(yīng)力場；②在計算域內(nèi)根據(jù)某種回歸方程直接計算每個單元體的初始應(yīng)力，并轉(zhuǎn)化為計算單元的節(jié)點(diǎn)力，然后根據(jù)位移和應(yīng)力邊界條件，平衡計算模型系統(tǒng)內(nèi)的應(yīng)力場；③在模型位移邊界條件和應(yīng)力邊界條件的基礎(chǔ)上生成初始應(yīng)力場，然后根據(jù)幾個點(diǎn)實測的地應(yīng)力對計算域的邊界條件進(jìn)行調(diào)整，并最終擬合到新的初始應(yīng)力場，使之與實測點(diǎn)的地應(yīng)力值較為吻合；④基于初始應(yīng)變能理論的初始地應(yīng)力場的反演方法，其基本思路為：認(rèn)為將有構(gòu)造應(yīng)力存在的地質(zhì)體作為部分具有彈性體屬性的意義賦上去考慮是合理的，因此在數(shù)值模擬反演初始地應(yīng)力場時，可以將構(gòu)造應(yīng)力作為內(nèi)力施加在模型內(nèi)部，即對模型內(nèi)所有單元體采用三角形應(yīng)力分布初值設(shè)置，而并非在邊界上施加，以此來模擬構(gòu)造應(yīng)力場[15-16]。

2 計算模型構(gòu)建及方案

在應(yīng)力場反演方法討論分析的基礎(chǔ)上，擬通過對比分析位移邊界條件、應(yīng)力邊界條件、混合邊界條件等情況下，計算模型內(nèi)垂直應(yīng)力與水平應(yīng)力的分布情況，從而探討模擬高水平構(gòu)造應(yīng)力場的最佳方法。在實際數(shù)值模擬過程中，為了盡可能消除模型尺寸過小所產(chǎn)生的邊界效應(yīng)，通常模型左右邊界向外擴(kuò)展一定距離，計算作為平面應(yīng)變問題考慮。通過反復(fù)試算，當(dāng)計算模型長3 000 m，高1 400 m，寬50 m時，在這一尺寸下，邊界效應(yīng)對模型的影響可以忽略。由于地表地形對模擬結(jié)果影響較小，因此，在模型的概化過程中，將地表概化成水平地表。計算采用Mohr-Coulomb本構(gòu)模型。

為了充分研究邊界條件對初始地應(yīng)力場分布的影響，計算模擬采用以下邊界條件[15-16]：①位移邊界條件，如圖1所示，即在數(shù)值模擬計算中，模型左右邊界的水平方向的位移固定，模型底部的位移固定，模型頂部設(shè)置為自由表面，其中G表示自重體積力；②應(yīng)力邊界條件，圖2所示，模型左右施加應(yīng)力梯形；③混合邊界條件，模型底部固定豎向位移，左右邊界分別為混合的三角形分布的水平分布面力和水平向位移約束；④基于初始應(yīng)變能理論的初始地應(yīng)力場的反演方法，其主要為位移邊界條件的基礎(chǔ)上，將構(gòu)造應(yīng)力作為內(nèi)力施加在模型內(nèi)部，即對模型內(nèi)所有單元體采用三角形應(yīng)力分布初值設(shè)置。計算所采用的力學(xué)參數(shù)為青灰色白云巖的巖體力學(xué)參數(shù)，如表1所示。

圖1 位移邊界條件示意圖

圖2 應(yīng)力邊界條件示意圖

3 應(yīng)力場計算反演結(jié)果對比分析

3.1 基于位移邊界條件的初始地應(yīng)力場的反演

在位移邊界條件下，反演得到垂直應(yīng)力和水平應(yīng)力分布分別如圖3所示。計算結(jié)果表明：應(yīng)力等值線水平展布，走向穩(wěn)定，等值線間間隔幾乎相等，從地表至模型底邊呈線性分布。但在計算模型內(nèi)，同一位置的自重應(yīng)力明顯大于水平應(yīng)力，約為3倍左右。計算的最大主應(yīng)力為豎直方向，垂直應(yīng)力為最大主應(yīng)力，表明此種方法僅能用于計算分析以自重應(yīng)力場為主導(dǎo)的模型，而不能模擬高水平構(gòu)造應(yīng)力場。

圖3 位移邊界條件下的應(yīng)力分布

3.2 基于應(yīng)力邊界條件的初始地應(yīng)力場的反演

一般情況下，水平地應(yīng)力場是在自重體積力存在的條件下由構(gòu)造運(yùn)動產(chǎn)生的一種應(yīng)力場，尤其是在高構(gòu)造應(yīng)力區(qū)，巖體的原巖應(yīng)力狀態(tài)特征總是與顯著的構(gòu)造應(yīng)力狀態(tài)特征相聯(lián)系。大量的實測結(jié)果表明，水平構(gòu)造應(yīng)力場一般最大主壓應(yīng)力σ1與最小主壓應(yīng)力σ3方向均接近水平方向，而中間主應(yīng)力σ2一般為豎直方向，在數(shù)值上等于γh(γ為重度，h為埋深)。在數(shù)值模擬反演中，自重應(yīng)力通常作為一種體積力進(jìn)行考慮并賦值，然而對于構(gòu)造應(yīng)力條件的處理辦法，特別是高水平構(gòu)造應(yīng)力場的反演方法存在一定的分歧。初始應(yīng)力場模擬是指復(fù)原其采礦工程活動前(如開挖)的應(yīng)力狀態(tài)，如果將構(gòu)造應(yīng)力作為邊界應(yīng)力加載在模型邊界上，正如2節(jié)討論(圖2)，計算時以三角形應(yīng)力施加在模型的左右邊界，并且保持左右邊界三角形應(yīng)力不變，通過邊界單元與模型內(nèi)部單元接觸面的節(jié)點(diǎn)力向內(nèi)部傳遞，并且使得模型外表面上的分布力與其內(nèi)力平衡。

上述應(yīng)力邊界條件下，模型的垂直應(yīng)力分布如圖4(a)所示。由圖4(a)可以看出，垂直應(yīng)力等值線圖在模型中部基本上呈水平走向，且線性增加，但模型兩端自重應(yīng)力場明顯大于理論值，且應(yīng)力等值線分布并不理想。圖4(b)為應(yīng)力邊界條件下水平應(yīng)力分布圖，在模型中部，同一位置水平應(yīng)力高于自重應(yīng)力，但模型上部和下部均未呈水平走向分布，并且模型中間偏下部分的最大主壓應(yīng)力遠(yuǎn)低于理論值，而模型中間偏上部分的最大主壓應(yīng)力又遠(yuǎn)高于理論值?？梢?，左右側(cè)邊界上采用應(yīng)力邊界條件不能夠較好地計算模擬水平高構(gòu)造應(yīng)力場。

圖4 應(yīng)力邊界條件下的應(yīng)力分布

3.3 基于混合邊界條件的初始地應(yīng)力場的反演

圖5(a)、圖5(b)分別為分布應(yīng)力施加在模型左側(cè)時混合邊界條件下垂直和水平應(yīng)力分布，圖6(a)、圖6(b)分別為分布應(yīng)力施加在模型右側(cè)時混合邊界條件下垂直和水平應(yīng)力分布。從圖5(a)、圖6(a)可以看出，遠(yuǎn)離模型施加三角形分布面力的端部，自重應(yīng)力等值線呈線性增加，一旦靠近該端部，自重應(yīng)力等值線便發(fā)生變異，尤其是在施加分布面力一側(cè)，豎直應(yīng)力等值線起伏較大，自重應(yīng)力遠(yuǎn)大于理論值。圖5(b)、圖6(b)中模型施加分布面力一側(cè)的較小范圍內(nèi)，模擬結(jié)果與預(yù)期結(jié)果較好，一旦離開該端部，水平應(yīng)力等值線迅速變異。具體表現(xiàn)為在水平位移約束側(cè)淺部等值線過于密集，但在深部等值線又過于稀疏。礦體附近的水平應(yīng)力低于預(yù)期結(jié)果，水平應(yīng)力等值線越往深部越偏斜，等值線間隔在不同位置變化大，表明模擬得到的水平應(yīng)力梯度變化不均勻。綜上分析可知，混合邊界條件既不適于模擬自重應(yīng)力場，也不適于模擬水平構(gòu)造應(yīng)力場。

圖5 混合邊界條件下應(yīng)力場分布(左側(cè)應(yīng)力邊界、右側(cè)位移邊界)

圖6 混合邊界條件下應(yīng)力場分布左位移邊界、右應(yīng)力邊界

3.4 基于初始應(yīng)變能理論的初始地應(yīng)力場的反演分析

在地質(zhì)力學(xué)中，構(gòu)造應(yīng)力場通常是指地殼內(nèi)一定范圍內(nèi)某一瞬時的應(yīng)力狀態(tài)。構(gòu)造應(yīng)力場中的應(yīng)力分布和變化是連續(xù)并且有規(guī)律的。根據(jù)傳統(tǒng)的彈性理論和彈塑性理論，在自重體積力和構(gòu)造運(yùn)動的共同作用下，產(chǎn)生了位移并且做了相應(yīng)功，在這一過程中，除塑性變形功不能轉(zhuǎn)化為應(yīng)變能外，地質(zhì)體主要以彈性變形的方式儲存應(yīng)變能。近年來，諸多學(xué)者從能量的角度和觀點(diǎn)出發(fā)，對巖石變形破壞過程中發(fā)生的能量積聚，耗散與釋放進(jìn)行了分析[15]，使得能量觀點(diǎn)在巖體變形失穩(wěn)和破壞分析中的應(yīng)用愈加廣泛。

假若將巖體視為彈性體，那么巖體中的某一點(diǎn)的微小單元的單位體積具有的應(yīng)變能可表示為[16-17]

(1)

式(1)可改寫為

Gεijεij=WV+WF

(2)

式中:E為彈性模量；ν為泊松比；G為剪切模量；W為應(yīng)變能；WV為體積變形應(yīng)變能；WF為形狀改變應(yīng)變能；θ為應(yīng)變，θ=ε11+ε22+ε33。

如式(1)、式(2)所示，工程開挖后，初始應(yīng)力場發(fā)生變化，二次應(yīng)力場重新分布形成過程中，受到開挖影響的地質(zhì)體中某單元體將會同時出現(xiàn)體積改變和形狀改變。即微小單元體體積改變彈性比能WV和形狀改變彈性比能WF一般都會隨之變化。也就是說，巖體發(fā)生的變形破壞，其本質(zhì)是能量轉(zhuǎn)化的過程[13-15]，在工程開挖過程中，受影響的地質(zhì)體中所儲存的彈性應(yīng)變能必然會隨之變化。

通過前面的討論，將有構(gòu)造應(yīng)力存在的地質(zhì)體作為部分具有彈性體屬性的意義上去考慮，在理論上是比較合理的。因此，在數(shù)值法反演初始地應(yīng)力場時，可以將構(gòu)造應(yīng)力作為內(nèi)力施加在模型內(nèi)部，即對模型內(nèi)所有單元體采用梯形應(yīng)力分布初值設(shè)置，以此來模擬高水平構(gòu)造應(yīng)力場，并采用位移邊界條件的情況下進(jìn)行初始高水平構(gòu)造應(yīng)力場的反演計算。

圖7(a)為初始地應(yīng)力場中垂直應(yīng)力等值線分布。從圖7(a)可以看出，垂直應(yīng)力從地表向下呈線性增加，應(yīng)力跡線水平展布，走向穩(wěn)定，自重應(yīng)力計算值與理論值較為吻合。圖7(b)為初始地應(yīng)力場中水平壓應(yīng)力等值線分布，水平應(yīng)力從地表向下呈線性增加，應(yīng)力跡線較平直，水平展布，等值線間隔也很均勻，走向穩(wěn)定，水平應(yīng)力模擬結(jié)果與預(yù)期吻合較好。

圖7 初始應(yīng)地應(yīng)力場中的應(yīng)力分布

圖8為初始地應(yīng)力場主壓應(yīng)力跡線分布。從局部放大圖上看到最大主應(yīng)力均為水平方向，鉛直方向為中間主應(yīng)力，滿足預(yù)期要實現(xiàn)的水平構(gòu)造應(yīng)力為最大主壓應(yīng)力的初始應(yīng)力場特征量的條件。

圖8 模型中部的主壓應(yīng)力邊跡線分布

4 基于FLAC3D的獅子山礦區(qū)高水平構(gòu)造應(yīng)力場的反演與驗證

4.1 獅子山礦區(qū)地應(yīng)力場測試結(jié)果

分別在礦區(qū)4個不同中段對應(yīng)的4個測點(diǎn)進(jìn)行了實測，通過對實測數(shù)據(jù)綜合分析，得出如下結(jié)果。

(1)4個測點(diǎn)的最大主應(yīng)力均位于近水平方向，各測點(diǎn)的最大主應(yīng)力方向均接近為NNW-SSE向(近垂直于礦體走向方向)，說明礦區(qū)以近水平構(gòu)造應(yīng)力場為主。

(2)各測點(diǎn)最大主應(yīng)力、垂直主應(yīng)力和最小主應(yīng)力的大小均隨深度基本成線性增長，并且垂直主應(yīng)力和最小主應(yīng)力的隨埋深的增加速率要遠(yuǎn)小于最大主應(yīng)力的增加速率，進(jìn)一步證實了礦區(qū)最大主應(yīng)力的影響具有明顯的主導(dǎo)性。通過實測結(jié)果分析，得到最大主應(yīng)力、中間主應(yīng)力和最小主應(yīng)力三個方向地應(yīng)力(單位：MPa)隨深度H(單位：m)的線性回歸方程：

σh,max=-0.016 3+0.051 1H

(3)

σh,min=-0.388 3+0.021H

(4)

σv=-0.601 8+0.030 7H

(5)

式中：σh,max為最大主應(yīng)力；σh,min為最小主應(yīng)力；σv為垂直主應(yīng)力。

(3)結(jié)合室內(nèi)巖石力學(xué)試驗結(jié)果和三維地應(yīng)力實測結(jié)果，按照高地應(yīng)力定量標(biāo)準(zhǔn)[18]，獅子山銅礦屬于高地應(yīng)力區(qū)。

4.2 獅子山銅礦地應(yīng)力場反演數(shù)值模型的構(gòu)建與計算參數(shù)

計算模型以礦體為中心，計算域上邊界取至地表。計算模型X方向(垂直礦體走向方向)長度1 400 m；模型Y方向(沿礦體走向方向)，長度1 500 m；模型Z方向為豎直方向，模型底部標(biāo)高687 m，模型最大高度1 430 m。模型共包括394 713個單元體和411 825個節(jié)點(diǎn)，計算模型和礦體形態(tài)分別如圖9、圖10所示。采用Mohr-Coulomb彈塑性本構(gòu)模型[15]。計算模型X方向兩端約束X方向位移，模型Y方向兩端約束Y方向位移，模型底部固定位移，模型頂部為自由邊界。地應(yīng)力按3.4節(jié)討論的高水平構(gòu)造應(yīng)力場的施加方法進(jìn)行反演，即在位移邊界條件下，對模型內(nèi)所有單元體采用三角形應(yīng)力分布初值設(shè)置，最大主應(yīng)力為X方向，最小主應(yīng)力為Y方向，豎直方向為中間主應(yīng)力。礦區(qū)各巖性宏觀巖體力學(xué)參數(shù)如表1所示。

圖9 礦區(qū)三維有限差分?jǐn)?shù)值模型

圖10 礦體形態(tài)

圖11分別為初始應(yīng)力平衡后模型最大主應(yīng)力和最小主應(yīng)力分布。初始應(yīng)力場反演計算表明，X方向(垂直于礦體走向方向)為最大主應(yīng)力方向，Y方向(沿礦體走向方向)為最小主應(yīng)力，Z方向(豎直方向)為中間主應(yīng)力方向，與現(xiàn)場三維地應(yīng)力實測主應(yīng)力方向一致。根據(jù)三維地應(yīng)力實測的4個測點(diǎn)位置，在模型內(nèi)部相應(yīng)位置布設(shè)4個監(jiān)測點(diǎn)，以分析初始應(yīng)力平衡后各監(jiān)測點(diǎn)的主應(yīng)力數(shù)值與地應(yīng)力實測值之間的關(guān)系。

圖11 初始應(yīng)力平衡后模型主應(yīng)力分布

為檢驗獅子山礦區(qū)地應(yīng)力反演結(jié)果的可靠性，繪制實測地應(yīng)力與反演地應(yīng)力對比曲線，如圖12所示。各測點(diǎn)最大主應(yīng)力反演值與實測地應(yīng)力值之間的最大相對誤差為-8.10%，為測點(diǎn)D3；4個測點(diǎn)中，D1測點(diǎn)最小主應(yīng)力反演值與實測地應(yīng)力值之間的相對誤差最大，為-18.29%。各測點(diǎn)中，中間主應(yīng)力反演值與實測地應(yīng)力之間最大相對誤差為9.25%，為測點(diǎn)D3。綜上可見，該地應(yīng)力場的反演方法具有一定的可靠性。

圖12 獅子山礦區(qū)4個測點(diǎn)實測地應(yīng)力與反演地應(yīng)力對比

5 結(jié)論

(1)在地應(yīng)力邊界條件討論的基礎(chǔ)上，采用FLAC3D對4種地應(yīng)力場反演方法和邊界條件進(jìn)行數(shù)值模擬對比分析，結(jié)果表明位移邊界條件在模擬自重應(yīng)力場中效果較好，但不能較好地模擬高水平構(gòu)造應(yīng)力場?；旌线吔鐥l件既不適于模擬自重應(yīng)力場，也不適于模擬水平構(gòu)造應(yīng)力場。

(2)采用基于初始應(yīng)變能理論的初始地應(yīng)力場的反演方法，即在位移邊界條件下，對計算域內(nèi)所有單元體采用梯形應(yīng)力分布初值設(shè)置，即設(shè)定初始彈性應(yīng)變能狀態(tài)，使?jié)M足或擬合應(yīng)力隨深度變化的實測值，模擬效果與實際吻合最好，能夠較好地模擬初始高水平構(gòu)造應(yīng)力場。

(3)建立了獅子山礦區(qū)復(fù)雜三維力學(xué)模型，以初始應(yīng)變能理論的初始地應(yīng)力場的反演方法進(jìn)行設(shè)置，數(shù)值反演表明，各測點(diǎn)主應(yīng)力的方向與實測地應(yīng)力方向一致，各測點(diǎn)主應(yīng)力的反演值與實測地應(yīng)力回歸值之間的誤差在允許范圍內(nèi)，表明礦區(qū)水平構(gòu)造應(yīng)力場反演具有一定的可靠性。