環(huán)北部灣廣東水資源配置工程深埋隧洞段初始地應(yīng)力場(chǎng)反演分析

代志坤,鄭 曉,游思琴

(中水珠江規(guī)劃勘測(cè)設(shè)計(jì)有限公司,廣東 廣州 510630)

地應(yīng)力是天然狀態(tài)下巖體中賦存的應(yīng)力,是地球發(fā)展運(yùn)動(dòng)的結(jié)果,它是地球歷次構(gòu)造運(yùn)動(dòng)產(chǎn)生的。巖漿活動(dòng)、風(fēng)化剝蝕、地殼運(yùn)動(dòng)等都是地應(yīng)力形成的重要因素,地應(yīng)力是以上多種因素的多次綜合作用之后的結(jié)果。目前還不可能做到根據(jù)地應(yīng)力形成的過程來求解可供工程使用的初始地應(yīng)力場(chǎng)[1-2]。但是初始地應(yīng)力是地下洞室設(shè)計(jì)的主要荷載之一,正確選取地應(yīng)力的量值和方向是設(shè)計(jì)工作的關(guān)鍵,且地應(yīng)力也是工程計(jì)算的重要初始條件之一,工程實(shí)踐中,初始地應(yīng)力場(chǎng)的確定是巖石力學(xué)的一個(gè)重要研究?jī)?nèi)容,準(zhǔn)確的查明初始地應(yīng)力分布規(guī)律對(duì)地下工程的建設(shè)有重大意義[3-4]。獲得工程區(qū)地應(yīng)力數(shù)據(jù)的最直接有效的方式是進(jìn)行地應(yīng)力測(cè)試,進(jìn)行大量的地應(yīng)力測(cè)試需要大量的時(shí)間及預(yù)算,且測(cè)量成果受測(cè)量誤差影響會(huì)有一定的離散性,只能反應(yīng)局部地應(yīng)力情況[5],尤其是在工程范圍較大時(shí),例如大型引調(diào)水工程、交通工程及鐵路工程等,線路動(dòng)輒幾十、上百千米,地形地貌及地質(zhì)構(gòu)造差異較大,通過大規(guī)模的地應(yīng)力測(cè)試來確定工程區(qū)內(nèi)的初始地應(yīng)力場(chǎng),所耗費(fèi)的時(shí)間和資金是難以承受的,在工程勘察期往往僅能在局部進(jìn)行有限的地應(yīng)力測(cè)試[6]。隨著計(jì)算機(jī)技術(shù)發(fā)展,利用數(shù)值模擬技術(shù),并結(jié)合有限的實(shí)測(cè)數(shù)據(jù),對(duì)工程區(qū)地應(yīng)力進(jìn)行反演分析,是一種獲得初始地應(yīng)力場(chǎ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)勢(shì)[9],本文采用多元線性回歸分析法,計(jì)算結(jié)果表明,使用該方法能夠得到合理的初始地應(yīng)力場(chǎng)數(shù)據(jù)。

本文利用有限差分軟件FLAC3D對(duì)環(huán)北部灣水資源配置工程輸水隧洞埋深最大的洞段建立三維地質(zhì)模型并進(jìn)行地應(yīng)力反演計(jì)算,結(jié)合勘察期地應(yīng)力測(cè)試成果,同時(shí)利用多元線性回歸分析法對(duì)工程區(qū)地應(yīng)力進(jìn)行回歸分析,并建立初始地應(yīng)力歸回方程,最后根據(jù)回歸分析結(jié)果利用FLAC3D重構(gòu)工程區(qū)初始地應(yīng)力,比對(duì)各監(jiān)測(cè)點(diǎn)的反演值及實(shí)測(cè)值后可知,采用該方法可以有效的反演分析工程區(qū)的初始地應(yīng)力分布,為該工程的深埋隧洞支護(hù)設(shè)計(jì)及施工開挖提供有效的理論及數(shù)據(jù)支撐,對(duì)工程建設(shè)實(shí)施有重要的指導(dǎo)意義。

1 反演分析的方法及原理

1.1 地應(yīng)力反演分析思路

基于地形地質(zhì)資料,結(jié)合對(duì)工程區(qū)地應(yīng)力場(chǎng)條件的認(rèn)識(shí),建立三維模型,分析并確定主要的地應(yīng)力影響因素及邊界條件,然后進(jìn)行數(shù)值模擬分析,選擇相應(yīng)的分析方法,建立地應(yīng)力實(shí)測(cè)值及模擬觀測(cè)值之間的聯(lián)系,進(jìn)而確定工程區(qū)初始地應(yīng)力場(chǎng)的分布規(guī)律,具體流程見圖1。

圖1 地應(yīng)力反演思路流程

1.2 模擬分析因素及邊界條件

本文選擇6種因素作為模擬巖體自重和地質(zhì)構(gòu)造力作用的基本因素[1]:①自重應(yīng)力;②X向均勻擠壓構(gòu)造運(yùn)動(dòng);③Y向均勻擠壓構(gòu)造運(yùn)動(dòng);④水平面內(nèi)的剪切構(gòu)造(XY面);⑤垂直平面內(nèi)的剪切變形構(gòu)造(XZ面);⑥垂直平面內(nèi)的剪切變形構(gòu)造(YZ面)。

對(duì)于邊界條件,有添加初始應(yīng)力和初始位移兩種方法,本文采取設(shè)置初始位移法來模擬地質(zhì)構(gòu)造運(yùn)動(dòng)。FLAC3D軟件無法添加初始位移,可以通過施加節(jié)點(diǎn)速度分量的方式來設(shè)置位移邊界條件。根據(jù)計(jì)算收斂情況設(shè)置在各種工況下的節(jié)點(diǎn)速率及計(jì)算步數(shù)。

工況1:僅考慮自重作用,見圖2a。

工況2:僅考慮x軸向擠壓構(gòu)造運(yùn)動(dòng),見圖2b,計(jì)算時(shí)以邊界位移為控制對(duì)象,控制單位取為1 cm,設(shè)置節(jié)點(diǎn)速率xvel= 2×10-6m/步,計(jì)算步數(shù)為5 000步。

工況3:僅考慮y軸向擠壓構(gòu)造運(yùn)動(dòng),見圖2c,計(jì)算時(shí)以邊界位移為控制對(duì)象,控制單位取為1 cm,設(shè)置節(jié)點(diǎn)速率yvel=2×10-6m/步,計(jì)算步數(shù)為5 000步。

工況4:僅考慮xy平面剪切構(gòu)造運(yùn)動(dòng),見圖2d,根據(jù)文獻(xiàn)[22],長(zhǎng)邊和短邊剪切位移應(yīng)呈相應(yīng)比例,長(zhǎng)邊位移取3.6 mm,短邊位移取6 mm;邊界節(jié)點(diǎn)速率分別為yvel=7.2×10-7m/步、xvel=1.2×10-6m/步,計(jì)算步數(shù)為5 000步。

工況5:僅考慮xz平面剪切構(gòu)造運(yùn)動(dòng),見圖2e,地形面自由,邊界剪切位移取9 mm,邊界節(jié)點(diǎn)速率zvel=1.6×10-6m/步,計(jì)算步數(shù)為5 000步。

工況6:僅考慮yz平面剪切構(gòu)造運(yùn)動(dòng),見圖2f,地形面自由,邊界剪切位移取9 mm,邊界節(jié)點(diǎn)速率zvel=1.6×10-6m/步,計(jì)算步數(shù)為5 000步。

a)自重作用

1.3 多元線性回歸分析原理[9]

(1)

假設(shè)有m個(gè)觀測(cè)點(diǎn),則最小二乘法的殘差平方和為式(2):

(2)

解此法方程式,得n個(gè)待定回歸系數(shù)L=(L1,L2,L3,…,Ln)T,則計(jì)算域內(nèi)任一點(diǎn)P的回歸初始應(yīng)力,可由該點(diǎn)邊界荷載工況有限元計(jì)算值迭加而得式(3):

(3)

2 工程區(qū)地應(yīng)力反演分析

2.1 工程概況

環(huán)北部灣廣東水資源配置工程位于廣東省西南部,該工程區(qū)主體所屬一級(jí)大地構(gòu)造單元為華南褶皺系,二級(jí)大地構(gòu)造單元為粵西隆起區(qū)(Ⅱ1)、雷瓊坳陷(Ⅱ6)、粵北、粵東北-粵中坳陷帶(Ⅱ5),其中主干線位于粵西隆起區(qū),分為3個(gè)三級(jí)大地構(gòu)造單元:大瑤山隆起(Ⅲ1)、羅定坳陷(Ⅲ2)、云開大山隆起區(qū)(Ⅲ3)。該工程自西江干流取水,采用隧洞、暗涵、倒虹吸、渡槽等多種建筑物形式,穿越云霧山山脈、羅定盆地及云開大山,輸水干線長(zhǎng)約200 km,其中輸水隧洞總長(zhǎng)占比約80%。該工程干線輸水隧洞屬于超長(zhǎng)深埋大洞徑水工隧洞,隧洞單段最長(zhǎng)約60 km,埋深一般均超過200 m,最大埋深910 m,開挖洞徑將近10 m。

2.2 計(jì)算范圍及計(jì)算模型

本文選取輸水隧洞埋深最大的區(qū)域進(jìn)行研究。根據(jù)工程區(qū)地質(zhì)條件、水文條件及實(shí)測(cè)點(diǎn)的分布情況,為了消除人工邊界誤差在重要結(jié)構(gòu)部位的影響,計(jì)算范圍確定如下:以實(shí)測(cè)點(diǎn)為中心區(qū)域,垂直隧道走向?yàn)閄軸,平行于隧洞走向?yàn)閅軸,豎直向上為Z軸;X軸及Y軸的計(jì)算范圍為600 m×1 000 m,Z軸方向從高程-90 m至自然地面。計(jì)算區(qū)域X軸及Y軸方向分別為S60°E、N30°E。計(jì)算區(qū)域范圍及測(cè)點(diǎn)平面位置見圖3。共劃分六面體單元120 000個(gè),節(jié)點(diǎn)128 061個(gè),計(jì)算網(wǎng)格劃分情況見圖4。假設(shè)工程區(qū)巖體為連續(xù)彈性體,采用線彈性本構(gòu)模型模擬巖體力學(xué)行為。FLAC3D線彈性本構(gòu)模型涉及的力學(xué)參數(shù)包括彈性模量、泊松比及密度(表1)。

圖3 計(jì)算區(qū)域及測(cè)點(diǎn)平面布置

表1 工程區(qū)區(qū)巖體力學(xué)參數(shù)

圖4 計(jì)算網(wǎng)格

2.3 工程區(qū)地應(yīng)力測(cè)試成果及初步分析

為了解工程沿線深埋段初始地應(yīng)力情況,在工程可行性研究階段,本工程勘察設(shè)計(jì)單位采用水壓致裂法對(duì)線路深埋段進(jìn)行了地應(yīng)力測(cè)量,對(duì)輸水隧洞深埋洞段的地應(yīng)力場(chǎng)有了基本認(rèn)識(shí)。計(jì)算范圍內(nèi)進(jìn)行地應(yīng)力測(cè)試的鉆孔為FSGZ24號(hào)孔,孔口高程963.1 m,孔深888.8 m。根據(jù)測(cè)孔巖芯判別,在巖體相對(duì)完整的深度區(qū)間進(jìn)行水壓致裂測(cè)試,成功獲得18段測(cè)試數(shù)據(jù),對(duì)其中具有代表性的5個(gè)測(cè)段進(jìn)行最大水平主應(yīng)力的方向印模,測(cè)試成果見表2。

表2 FSGZ-24鉆孔水壓致裂法地應(yīng)力測(cè)試結(jié)果

由實(shí)測(cè)地應(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個(gè)主應(yīng)力主要呈σz>σH>σh特征,說明該區(qū)地應(yīng)力場(chǎng)以自重應(yīng)力為主。依據(jù)壓裂縫方向印模結(jié)果,鉆孔最大水平主應(yīng)力方向穩(wěn)定在N44°W～N72°W。測(cè)試方向比較一致,主要呈NW—NWW。

水平主應(yīng)力的大小與孔深關(guān)系見圖5?？梢钥闯?除了局部測(cè)點(diǎn),測(cè)試結(jié)果總體沿孔深增加而增大,應(yīng)力測(cè)試數(shù)據(jù)與深度的擬合關(guān)系見式(4):

圖5 FSGZ-24鉆孔的主應(yīng)力量值與孔深關(guān)系

(4)

由于建立模型劃分網(wǎng)格的精細(xì)程度所限,部分計(jì)算監(jiān)測(cè)點(diǎn)與實(shí)測(cè)地應(yīng)力位置偏差較大,為確保實(shí)測(cè)數(shù)據(jù)的可靠性,剔除與模型計(jì)算監(jiān)測(cè)點(diǎn)位置偏差較大測(cè)點(diǎn),選取以下13個(gè)點(diǎn)作為反演計(jì)算區(qū)域內(nèi)初始地應(yīng)力場(chǎng)的依據(jù),并根據(jù)彈性力學(xué)理論,將實(shí)測(cè)地應(yīng)力值轉(zhuǎn)成計(jì)算坐標(biāo)系下地應(yīng)力實(shí)測(cè)值(表3)。

表3 計(jì)算坐標(biāo)系下地應(yīng)力實(shí)測(cè)值

3 反演計(jì)算及結(jié)果分析

三維模型建立之后,按自重作用、X向擠壓、Y向擠壓、XY平面內(nèi)水平剪切、XZ垂直平面內(nèi)剪切、YZ垂直平面內(nèi)剪切6種因素單獨(dú)作用分別進(jìn)行模擬計(jì)算,在6種工況下沿隧洞軸線方向垂直剖面上的最大最小主應(yīng)力情況見圖6—11。在計(jì)算過程中,在實(shí)測(cè)地應(yīng)力點(diǎn)處設(shè)置計(jì)算監(jiān)測(cè)點(diǎn),監(jiān)測(cè)在各因素作用下的地應(yīng)力計(jì)算值。

a)最大主應(yīng)力

a)最大主應(yīng)力

a)最大主應(yīng)力

a)最大主應(yīng)力

a)最大主應(yīng)力

a)最大主應(yīng)力

3.1 線型回歸分析步驟

步驟一:采用最小二乘法進(jìn)行多元線型回歸確定回歸方程的系數(shù)。

步驟二:通過F檢驗(yàn)是判斷是否存在顯著的線性關(guān)系。

步驟三:通過決定系數(shù)R2值分析模型擬合情況,同時(shí)對(duì)VIF值進(jìn)行分析,確定模型的共線性。

步驟四:分析其顯著性。

步驟五:確定模型公式(圖12)。

圖12 回歸分析流程

3.2 線性回歸分析結(jié)果

對(duì)圖2中確定的 6 種工況利用 FLAC3D分別進(jìn)行模擬計(jì)算,然后以實(shí)測(cè)地應(yīng)力結(jié)果(表2)在計(jì)算坐標(biāo)系下的應(yīng)力分量(表3)為回歸目標(biāo)。利用最小二乘法進(jìn)行多元線性回歸,求得6個(gè)自變量的回歸系數(shù)L1=1.275,L2=8.968,L3=2.247,L4=-7.113,L5=-85.48,L6=185.77,觀測(cè)誤差εk=3.71。曲線回歸的擬合程度R2=0.991,表明回歸公式的相關(guān)性較好。則該工程區(qū)初始地應(yīng)力場(chǎng)回歸方程為:

σd=1.275σz+8.968σx+2.247σy-7.113τxy-85.48τxz+185.77τyz+3.715

從F檢驗(yàn)的結(jié)果分析可以得到,顯著性P<0.001,水平上呈現(xiàn)顯著性,因此模型基本滿足要求,對(duì)于變量共線性表現(xiàn)顯著度檢驗(yàn)觀測(cè)值F=830.123,可以認(rèn)為該 6個(gè)自變量的總體效果顯著。

通過應(yīng)力場(chǎng)平衡計(jì)算,得到計(jì)算坐標(biāo)系下各測(cè)點(diǎn)的回歸應(yīng)力值,測(cè)點(diǎn)的地應(yīng)力實(shí)測(cè)值與回歸主應(yīng)力值對(duì)比見圖13。

圖13 地應(yīng)力實(shí)測(cè)值與回歸值對(duì)比直方

3.3 初始地應(yīng)力重構(gòu)

根據(jù)線性回歸系數(shù)及各工況構(gòu)造相應(yīng)模型重構(gòu)研究區(qū)域初始地應(yīng)力模型,從而得到工程區(qū)任意一點(diǎn)的地應(yīng)力狀態(tài)。各種工況疊加過程中,需要及時(shí)調(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)算,各測(cè)點(diǎn)反演值及實(shí)測(cè)值對(duì)比見表4。

表4 各測(cè)點(diǎn)實(shí)測(cè)地應(yīng)力值與反演值對(duì)比

a)最大主應(yīng)力

3.4 結(jié)果分析

由圖13、14及表4可以看出,回歸結(jié)果、反演結(jié)果與實(shí)測(cè)結(jié)果基本吻合,反演水平最大主應(yīng)力方向?yàn)镹46°W～N54°W,實(shí)測(cè)地應(yīng)力最大水平主應(yīng)力方向?yàn)镹44°W～N72°W,反演結(jié)果與實(shí)測(cè)結(jié)果符合程度較好,說明本次反演較好地模擬了該工程區(qū)初始地應(yīng)力場(chǎng)分布情況。埋深較小處地應(yīng)力大小受地形起伏影響較大,隨著埋深的逐漸增大,初始地應(yīng)力分布受地形影響越來越小。

4 結(jié)論

a)基于最小二乘法的多元線性回歸,采用有限差分軟件FLAC3D進(jìn)行反演分析本工程的初始地應(yīng)力場(chǎng)分布,反演值與實(shí)測(cè)值擬合程度較好,可為工程設(shè)計(jì)及施工提供有意義的初始地應(yīng)力數(shù)據(jù)。

b)實(shí)測(cè)最大水平主應(yīng)力方向的側(cè)壓系數(shù)(σH/σz)為0.7～1.0,反演分析最大水平主應(yīng)力方向的側(cè)壓系數(shù)(σH/σz)為0.81～1.00,且3個(gè)主應(yīng)力均呈σz>σH>σh特征,說明該區(qū)地應(yīng)力場(chǎng)以自重應(yīng)力為主。

c)實(shí)測(cè)最大水平主應(yīng)力方向?yàn)镹44°W～N72°W,反演水平最大主應(yīng)力方向?yàn)镹46°W～N54°W;反演區(qū)域輸水隧洞走向?yàn)镾30°W,根據(jù)實(shí)測(cè)成果,隧洞走向與最大水平主應(yīng)力方向夾角為74～102°,根據(jù)反演結(jié)果,隧洞走向與最大水平主應(yīng)力方向夾角為76～84°,均為大角度相交,這對(duì)隧洞圍巖穩(wěn)定性不利,隧洞開挖時(shí),兩側(cè)洞壁存在較為突出的穩(wěn)定性問題,支護(hù)設(shè)計(jì)時(shí)應(yīng)予以考慮。

d)由反演成果可知,地表凹陷處應(yīng)力集中明顯(圖14),隨著埋深的加大,地形起伏對(duì)地應(yīng)力分布的影響逐漸減小,應(yīng)力集中現(xiàn)象逐步消失,地應(yīng)力分布趨于平穩(wěn)。鑒于此種情況,對(duì)于自重應(yīng)力為主的深埋工程,在建立三維地質(zhì)模型進(jìn)行數(shù)值模擬分析時(shí),可以根據(jù)實(shí)際情況對(duì)地質(zhì)模型進(jìn)行適當(dāng)簡(jiǎn)化,如不考慮地形起伏對(duì)分析成果的影響,可以極大降低模型建立的難度。