代志坤,鄭 曉,游思琴
(中水珠江規(guī)劃勘測設(shè)計有限公司,廣東 廣州 510630)
地應(yīng)力是天然狀態(tài)下巖體中賦存的應(yīng)力,是地球發(fā)展運(yùn)動的結(jié)果,它是地球歷次構(gòu)造運(yùn)動產(chǎn)生的。巖漿活動、風(fēng)化剝蝕、地殼運(yùn)動等都是地應(yīng)力形成的重要因素,地應(yīng)力是以上多種因素的多次綜合作用之后的結(jié)果。目前還不可能做到根據(jù)地應(yīng)力形成的過程來求解可供工程使用的初始地應(yīng)力場[1-2]。但是初始地應(yīng)力是地下洞室設(shè)計的主要荷載之一,正確選取地應(yīng)力的量值和方向是設(shè)計工作的關(guān)鍵,且地應(yīng)力也是工程計算的重要初始條件之一,工程實踐中,初始地應(yīng)力場的確定是巖石力學(xué)的一個重要研究內(nèi)容,準(zhǔn)確的查明初始地應(yīng)力分布規(guī)律對地下工程的建設(shè)有重大意義[3-4]。獲得工程區(qū)地應(yīng)力數(shù)據(jù)的最直接有效的方式是進(jìn)行地應(yīng)力測試,進(jìn)行大量的地應(yīng)力測試需要大量的時間及預(yù)算,且測量成果受測量誤差影響會有一定的離散性,只能反應(yīng)局部地應(yīng)力情況[5],尤其是在工程范圍較大時,例如大型引調(diào)水工程、交通工程及鐵路工程等,線路動輒幾十、上百千米,地形地貌及地質(zhì)構(gòu)造差異較大,通過大規(guī)模的地應(yīng)力測試來確定工程區(qū)內(nèi)的初始地應(yīng)力場,所耗費(fèi)的時間和資金是難以承受的,在工程勘察期往往僅能在局部進(jìn)行有限的地應(yīng)力測試[6]。隨著計算機(jī)技術(shù)發(fā)展,利用數(shù)值模擬技術(shù),并結(jié)合有限的實測數(shù)據(jù),對工程區(qū)地應(yīng)力進(jìn)行反演分析,是一種獲得初始地應(yīng)力場分布規(guī)律的有效方法[7-14]。
目前,常用的地應(yīng)力反演方法有邊界荷載調(diào)整法、應(yīng)力(位移) 函數(shù)法、多元線性回歸分析法、偏最小二乘法、位移反分析法、快速應(yīng)力邊界法(S-B法)、神經(jīng)網(wǎng)絡(luò)、灰色理論、遺傳算法、變差函數(shù)法、模糊理論等方法[15-21]。
鑒于多元線性回歸分析法具有方便快捷、易于掌握、惟一解的優(yōu)勢[9],本文采用多元線性回歸分析法,計算結(jié)果表明,使用該方法能夠得到合理的初始地應(yīng)力場數(shù)據(jù)。
本文利用有限差分軟件FLAC3D對環(huán)北部灣水資源配置工程輸水隧洞埋深最大的洞段建立三維地質(zhì)模型并進(jìn)行地應(yīng)力反演計算,結(jié)合勘察期地應(yīng)力測試成果,同時利用多元線性回歸分析法對工程區(qū)地應(yīng)力進(jìn)行回歸分析,并建立初始地應(yīng)力歸回方程,最后根據(jù)回歸分析結(jié)果利用FLAC3D重構(gòu)工程區(qū)初始地應(yīng)力,比對各監(jiān)測點(diǎn)的反演值及實測值后可知,采用該方法可以有效的反演分析工程區(qū)的初始地應(yīng)力分布,為該工程的深埋隧洞支護(hù)設(shè)計及施工開挖提供有效的理論及數(shù)據(jù)支撐,對工程建設(shè)實施有重要的指導(dǎo)意義。
基于地形地質(zhì)資料,結(jié)合對工程區(qū)地應(yīng)力場條件的認(rèn)識,建立三維模型,分析并確定主要的地應(yīng)力影響因素及邊界條件,然后進(jìn)行數(shù)值模擬分析,選擇相應(yīng)的分析方法,建立地應(yīng)力實測值及模擬觀測值之間的聯(lián)系,進(jìn)而確定工程區(qū)初始地應(yīng)力場的分布規(guī)律,具體流程見圖1。
圖1 地應(yīng)力反演思路流程
本文選擇6種因素作為模擬巖體自重和地質(zhì)構(gòu)造力作用的基本因素[1]:①自重應(yīng)力;②X向均勻擠壓構(gòu)造運(yùn)動;③Y向均勻擠壓構(gòu)造運(yùn)動;④水平面內(nèi)的剪切構(gòu)造(XY面);⑤垂直平面內(nèi)的剪切變形構(gòu)造(XZ面);⑥垂直平面內(nèi)的剪切變形構(gòu)造(YZ面)。
對于邊界條件,有添加初始應(yīng)力和初始位移兩種方法,本文采取設(shè)置初始位移法來模擬地質(zhì)構(gòu)造運(yùn)動。FLAC3D軟件無法添加初始位移,可以通過施加節(jié)點(diǎn)速度分量的方式來設(shè)置位移邊界條件。根據(jù)計算收斂情況設(shè)置在各種工況下的節(jié)點(diǎn)速率及計算步數(shù)。
工況1:僅考慮自重作用,見圖2a。
工況2:僅考慮x軸向擠壓構(gòu)造運(yùn)動,見圖2b,計算時以邊界位移為控制對象,控制單位取為1 cm,設(shè)置節(jié)點(diǎn)速率xvel= 2×10-6m/步,計算步數(shù)為5 000步。
工況3:僅考慮y軸向擠壓構(gòu)造運(yùn)動,見圖2c,計算時以邊界位移為控制對象,控制單位取為1 cm,設(shè)置節(jié)點(diǎn)速率yvel=2×10-6m/步,計算步數(shù)為5 000步。
工況4:僅考慮xy平面剪切構(gòu)造運(yùn)動,見圖2d,根據(jù)文獻(xiàn)[22],長邊和短邊剪切位移應(yīng)呈相應(yīng)比例,長邊位移取3.6 mm,短邊位移取6 mm;邊界節(jié)點(diǎn)速率分別為yvel=7.2×10-7m/步、xvel=1.2×10-6m/步,計算步數(shù)為5 000步。
工況5:僅考慮xz平面剪切構(gòu)造運(yùn)動,見圖2e,地形面自由,邊界剪切位移取9 mm,邊界節(jié)點(diǎn)速率zvel=1.6×10-6m/步,計算步數(shù)為5 000步。
工況6:僅考慮yz平面剪切構(gòu)造運(yùn)動,見圖2f,地形面自由,邊界剪切位移取9 mm,邊界節(jié)點(diǎn)速率zvel=1.6×10-6m/步,計算步數(shù)為5 000步。
a)自重作用
(1)
假設(shè)有m個觀測點(diǎn),則最小二乘法的殘差平方和為式(2):
(2)
解此法方程式,得n個待定回歸系數(shù)L=(L1,L2,L3,…,Ln)T,則計算域內(nèi)任一點(diǎn)P的回歸初始應(yīng)力,可由該點(diǎn)邊界荷載工況有限元計算值迭加而得式(3):
(3)
環(huán)北部灣廣東水資源配置工程位于廣東省西南部,該工程區(qū)主體所屬一級大地構(gòu)造單元為華南褶皺系,二級大地構(gòu)造單元為粵西隆起區(qū)(Ⅱ1)、雷瓊坳陷(Ⅱ6)、粵北、粵東北-粵中坳陷帶(Ⅱ5),其中主干線位于粵西隆起區(qū),分為3個三級大地構(gòu)造單元:大瑤山隆起(Ⅲ1)、羅定坳陷(Ⅲ2)、云開大山隆起區(qū)(Ⅲ3)。該工程自西江干流取水,采用隧洞、暗涵、倒虹吸、渡槽等多種建筑物形式,穿越云霧山山脈、羅定盆地及云開大山,輸水干線長約200 km,其中輸水隧洞總長占比約80%。該工程干線輸水隧洞屬于超長深埋大洞徑水工隧洞,隧洞單段最長約60 km,埋深一般均超過200 m,最大埋深910 m,開挖洞徑將近10 m。
本文選取輸水隧洞埋深最大的區(qū)域進(jìn)行研究。根據(jù)工程區(qū)地質(zhì)條件、水文條件及實測點(diǎn)的分布情況,為了消除人工邊界誤差在重要結(jié)構(gòu)部位的影響,計算范圍確定如下:以實測點(diǎn)為中心區(qū)域,垂直隧道走向為X軸,平行于隧洞走向為Y軸,豎直向上為Z軸;X軸及Y軸的計算范圍為600 m×1 000 m,Z軸方向從高程-90 m至自然地面。計算區(qū)域X軸及Y軸方向分別為S60°E、N30°E。計算區(qū)域范圍及測點(diǎn)平面位置見圖3。共劃分六面體單元120 000個,節(jié)點(diǎn)128 061個,計算網(wǎng)格劃分情況見圖4。假設(shè)工程區(qū)巖體為連續(xù)彈性體,采用線彈性本構(gòu)模型模擬巖體力學(xué)行為。FLAC3D線彈性本構(gòu)模型涉及的力學(xué)參數(shù)包括彈性模量、泊松比及密度(表1)。
圖3 計算區(qū)域及測點(diǎn)平面布置
表1 工程區(qū)區(qū)巖體力學(xué)參數(shù)
圖4 計算網(wǎng)格
為了解工程沿線深埋段初始地應(yīng)力情況,在工程可行性研究階段,本工程勘察設(shè)計單位采用水壓致裂法對線路深埋段進(jìn)行了地應(yīng)力測量,對輸水隧洞深埋洞段的地應(yīng)力場有了基本認(rèn)識。計算范圍內(nèi)進(jìn)行地應(yīng)力測試的鉆孔為FSGZ24號孔,孔口高程963.1 m,孔深888.8 m。根據(jù)測孔巖芯判別,在巖體相對完整的深度區(qū)間進(jìn)行水壓致裂測試,成功獲得18段測試數(shù)據(jù),對其中具有代表性的5個測段進(jìn)行最大水平主應(yīng)力的方向印模,測試成果見表2。
表2 FSGZ-24鉆孔水壓致裂法地應(yīng)力測試結(jié)果
由實測地應(yīng)力成果可知,在孔深342.2~882.2 m(高程80.9~620.9 m)范圍內(nèi),最大水平主應(yīng)力為7.1~22.4MPa,最小水平主應(yīng)力為5.9~16.8 MPa,自重應(yīng)力為9.1~23.4 MPa。最大水平主應(yīng)力方向的側(cè)壓系數(shù)(σH/σz)為0.7~1.0,3個主應(yīng)力主要呈σz>σH>σh特征,說明該區(qū)地應(yīng)力場以自重應(yīng)力為主。依據(jù)壓裂縫方向印模結(jié)果,鉆孔最大水平主應(yīng)力方向穩(wěn)定在N44°W~N72°W。測試方向比較一致,主要呈NW—NWW。
水平主應(yīng)力的大小與孔深關(guān)系見圖5??梢钥闯?除了局部測點(diǎn),測試結(jié)果總體沿孔深增加而增大,應(yīng)力測試數(shù)據(jù)與深度的擬合關(guān)系見式(4):
圖5 FSGZ-24鉆孔的主應(yīng)力量值與孔深關(guān)系
(4)
由于建立模型劃分網(wǎng)格的精細(xì)程度所限,部分計算監(jiān)測點(diǎn)與實測地應(yīng)力位置偏差較大,為確保實測數(shù)據(jù)的可靠性,剔除與模型計算監(jiān)測點(diǎn)位置偏差較大測點(diǎn),選取以下13個點(diǎn)作為反演計算區(qū)域內(nèi)初始地應(yīng)力場的依據(jù),并根據(jù)彈性力學(xué)理論,將實測地應(yīng)力值轉(zhuǎn)成計算坐標(biāo)系下地應(yīng)力實測值(表3)。
表3 計算坐標(biāo)系下地應(yīng)力實測值
三維模型建立之后,按自重作用、X向擠壓、Y向擠壓、XY平面內(nèi)水平剪切、XZ垂直平面內(nèi)剪切、YZ垂直平面內(nèi)剪切6種因素單獨(dú)作用分別進(jìn)行模擬計算,在6種工況下沿隧洞軸線方向垂直剖面上的最大最小主應(yīng)力情況見圖6—11。在計算過程中,在實測地應(yīng)力點(diǎn)處設(shè)置計算監(jiān)測點(diǎn),監(jiān)測在各因素作用下的地應(yīng)力計算值。
a)最大主應(yīng)力
a)最大主應(yīng)力
a)最大主應(yīng)力
a)最大主應(yīng)力
a)最大主應(yīng)力
a)最大主應(yīng)力
步驟一:采用最小二乘法進(jìn)行多元線型回歸確定回歸方程的系數(shù)。
步驟二:通過F檢驗是判斷是否存在顯著的線性關(guān)系。
步驟三:通過決定系數(shù)R2值分析模型擬合情況,同時對VIF值進(jìn)行分析,確定模型的共線性。
步驟四:分析其顯著性。
步驟五:確定模型公式(圖12)。
圖12 回歸分析流程
對圖2中確定的 6 種工況利用 FLAC3D分別進(jìn)行模擬計算,然后以實測地應(yīng)力結(jié)果(表2)在計算坐標(biāo)系下的應(yīng)力分量(表3)為回歸目標(biāo)。利用最小二乘法進(jìn)行多元線性回歸,求得6個自變量的回歸系數(shù)L1=1.275,L2=8.968,L3=2.247,L4=-7.113,L5=-85.48,L6=185.77,觀測誤差εk=3.71。曲線回歸的擬合程度R2=0.991,表明回歸公式的相關(guān)性較好。則該工程區(qū)初始地應(yīng)力場回歸方程為:
σd=1.275σz+8.968σx+2.247σy-7.113τxy-85.48τxz+185.77τyz+3.715
從F檢驗的結(jié)果分析可以得到,顯著性P<0.001,水平上呈現(xiàn)顯著性,因此模型基本滿足要求,對于變量共線性表現(xiàn)顯著度檢驗觀測值F=830.123,可以認(rèn)為該 6個自變量的總體效果顯著。
通過應(yīng)力場平衡計算,得到計算坐標(biāo)系下各測點(diǎn)的回歸應(yīng)力值,測點(diǎn)的地應(yīng)力實測值與回歸主應(yīng)力值對比見圖13。
圖13 地應(yīng)力實測值與回歸值對比直方
根據(jù)線性回歸系數(shù)及各工況構(gòu)造相應(yīng)模型重構(gòu)研究區(qū)域初始地應(yīng)力模型,從而得到工程區(qū)任意一點(diǎn)的地應(yīng)力狀態(tài)。各種工況疊加過程中,需要及時調(diào)整邊界條件。各工況疊加順序及約束條件如下。
a)自重應(yīng)力,約束X向兩側(cè)邊界的X向位移、約束Y向兩側(cè)邊界的Y向位移、約束底部邊界的Z向位移。
b)X向擠壓,釋放X向一側(cè)邊界約束,并添加X向速率。
c)Y向擠壓,清除X向速率,重新約束X向側(cè)邊界,釋放Y向一側(cè)邊界約束,并在此添加Y向速率。
d)XY向剪切,清除Y向速率,釋放X向及Y向兩側(cè)約束,分別在X向兩側(cè)添加大小相同、方向相反的Y向速率,在Y向兩側(cè)添加大小相同、方向相反的X向速率。
e)X向垂直平面內(nèi)剪切,釋放底部Z向約束,清除X向及Y向速率,在X向兩側(cè)添加X向約束,在Y向兩側(cè)添加Y向約束,在X向兩側(cè)分別添加大小相同、方向相反的Z向速率。
f)Y向垂直平面內(nèi)剪切,清除Z向速率,在Y向兩側(cè)分別添加大小相同、方向相反的Z向速率。
按以上逐步疊加6種工況并進(jìn)行運(yùn)算,各測點(diǎn)反演值及實測值對比見表4。
表4 各測點(diǎn)實測地應(yīng)力值與反演值對比
a)最大主應(yīng)力
由圖13、14及表4可以看出,回歸結(jié)果、反演結(jié)果與實測結(jié)果基本吻合,反演水平最大主應(yīng)力方向為N46°W~N54°W,實測地應(yīng)力最大水平主應(yīng)力方向為N44°W~N72°W,反演結(jié)果與實測結(jié)果符合程度較好,說明本次反演較好地模擬了該工程區(qū)初始地應(yīng)力場分布情況。埋深較小處地應(yīng)力大小受地形起伏影響較大,隨著埋深的逐漸增大,初始地應(yīng)力分布受地形影響越來越小。
a)基于最小二乘法的多元線性回歸,采用有限差分軟件FLAC3D進(jìn)行反演分析本工程的初始地應(yīng)力場分布,反演值與實測值擬合程度較好,可為工程設(shè)計及施工提供有意義的初始地應(yīng)力數(shù)據(jù)。
b)實測最大水平主應(yīng)力方向的側(cè)壓系數(shù)(σH/σz)為0.7~1.0,反演分析最大水平主應(yīng)力方向的側(cè)壓系數(shù)(σH/σz)為0.81~1.00,且3個主應(yīng)力均呈σz>σH>σh特征,說明該區(qū)地應(yīng)力場以自重應(yīng)力為主。
c)實測最大水平主應(yīng)力方向為N44°W~N72°W,反演水平最大主應(yīng)力方向為N46°W~N54°W;反演區(qū)域輸水隧洞走向為S30°W,根據(jù)實測成果,隧洞走向與最大水平主應(yīng)力方向夾角為74~102°,根據(jù)反演結(jié)果,隧洞走向與最大水平主應(yīng)力方向夾角為76~84°,均為大角度相交,這對隧洞圍巖穩(wěn)定性不利,隧洞開挖時,兩側(cè)洞壁存在較為突出的穩(wěn)定性問題,支護(hù)設(shè)計時應(yīng)予以考慮。
d)由反演成果可知,地表凹陷處應(yīng)力集中明顯(圖14),隨著埋深的加大,地形起伏對地應(yīng)力分布的影響逐漸減小,應(yīng)力集中現(xiàn)象逐步消失,地應(yīng)力分布趨于平穩(wěn)。鑒于此種情況,對于自重應(yīng)力為主的深埋工程,在建立三維地質(zhì)模型進(jìn)行數(shù)值模擬分析時,可以根據(jù)實際情況對地質(zhì)模型進(jìn)行適當(dāng)簡化,如不考慮地形起伏對分析成果的影響,可以極大降低模型建立的難度。