楊繼華， 齊三紅， 郭衛(wèi)新， 張黨立

(黃河勘測規(guī)劃設(shè)計有限公司， 河南 鄭州 450003)

0 引言

目前，國內(nèi)外水電站地下廠房呈現(xiàn)出大跨度(20～30 m)、高邊墻(40～80 m)的技術(shù)特點，部分斷面面積大于1 000 m2，為超大斷面地下工程。水電站大型地下廠房結(jié)構(gòu)復雜，其受巖性、斷層、節(jié)理、巖脈、地下水及地應力等因素的影響遠大于普通隧洞[1-2]。

圍巖的力學參數(shù)及地應力場特征(如彈性模量、泊松比、側(cè)壓力系數(shù)等)是地下廠房支護設(shè)計和穩(wěn)定性分析的重要依據(jù)。巖體作為一種地質(zhì)體，其成因極為復雜，且一般經(jīng)歷過多次構(gòu)造作用，具有不連續(xù)性、非均質(zhì)性和各向異性的特點，在力學性質(zhì)上表現(xiàn)出明顯的非線性特征[3]。

目前，巖體力學參數(shù)的確定主要有解析法、試驗法及理論與實測相結(jié)合的方法等3種[4]。解析法先引入一定的假定條件，但假設(shè)往往與實際情況會有所差別，因此計算的參數(shù)存在偏差的可能性較大； 試驗可分為室內(nèi)試驗和現(xiàn)場試驗，但試驗結(jié)果只能反映試驗點處的巖體力學參數(shù)，無法代表整個工程區(qū)域的巖體特性，且試驗費用昂貴、時間較長； 理論與實測相結(jié)合的方法一般利用實測的數(shù)據(jù)，采用理論分析的方法確定巖體的力學參數(shù)，又稱為反分析法，目前常用的反分析法主要是位移反分析法。位移是反映巖土體受力變形的重要物理量，是巖土體穩(wěn)定性最直接和最宏觀的表現(xiàn)，位移反分析法融合了最優(yōu)化理論和數(shù)值計算等先進方法和計算技術(shù)，具有理論基礎(chǔ)完善、計算精度高、可以進行非線性分析、時間和經(jīng)濟成本低等優(yōu)點，能彌補傳統(tǒng)研究方法的不足，具有很大的實用價值[5-8]。劉英棨等[9]提出以水平收斂與拱底沉降相結(jié)合的多位移反分析法，建立平面簡化有限元模型，實時分析隧道支護結(jié)構(gòu)受力狀態(tài)； 江宗斌等[10]提出了基于位移-應力的多元信息聯(lián)合知能反分析方法； 張研等[11]將粒子群優(yōu)化算法與高斯過程機器學習方法相融合，結(jié)合FLAC3D數(shù)值計算程序，提出隧洞圍巖位移優(yōu)化反分析的粒子群-高斯過程智能協(xié)同優(yōu)化方法； 文建華等[12]采用分層運算方法對地下洞室進行了黏彈性位移反分析。綜上所述，目前研究多集中在隧洞等小斷面洞室方面，對大斷面地下廠房涉及較少。另外，在位移反分析中，多采用監(jiān)測的位移數(shù)據(jù)，但實際上在洞室開挖后位移監(jiān)測儀器安裝之前，大部分位移已經(jīng)發(fā)生且無法監(jiān)測到。因此，采用位移增量進行反分析更能反映圍巖的物理力學性質(zhì)。

最小二乘支持向量機做為一種新的機器學習方法，在處理小樣本、非線性問題上有獨特的優(yōu)勢。粒子群算法作為一種新型仿生進化算法，有全局優(yōu)化、收斂速度快等優(yōu)點[13]。本文以正交設(shè)計、最小二乘支持向量機和粒子群算法等現(xiàn)代數(shù)學方法為基本手段，建立圍巖力學參數(shù)分析方法，以多點位移計實測位移增量值為依據(jù)，以厄瓜多爾Coca-Codo Sinclair水電站(簡稱CCS水電站)大型地下廠房為工程背景，采用二維彈塑性有限元方法建立地質(zhì)結(jié)構(gòu)分析模型，對CCS水電站地下廠房區(qū)域巖體力學特性及地應力場特征進行反分析。

1 最小二乘支持向量機與粒子群算法位移反分析方法

1.1 目標函數(shù)

位移反分析一般是通過某種方法，使得一組待分析的參數(shù)及其相應的位移值逼近實測位移值。對于工程實際，通過逼近達到總體上的最優(yōu)效果，因此目標函數(shù)可用式(1)表示：

(1)

式中：x=(x1,x2,…,xk,…,xm)，xk為待反分析的巖體參數(shù)，如彈性模量E、黏聚力c、內(nèi)摩擦角φ等，m為參數(shù)的個數(shù)；fi(x)為第i個測點上位移計算值；ui為相應的實測位移值；n為測點數(shù)。

1.2 最小二乘支持向量機

最小二乘支持向量機采用不同的優(yōu)化目標函數(shù)，并且用等式約束取代不等式約束，可用于回歸或模式識別等問題的解決[14-16]。假定訓練集為(xi,yi)，其中：i=1, 2,…,k；xi∈Rn為n維系統(tǒng)輸入向量；yi∈R為相應的輸出值。如果問題為非線性，可將非線性映射輸入向量映射到高維空間，轉(zhuǎn)化為類似線性問題加以解決。最小二乘支持向量機優(yōu)化問題的最小化函數(shù)為：

(2)

式中：w為權(quán)向量；C為懲罰因子，即調(diào)節(jié)常數(shù)。

約束條件為：

yi-wφ(xi)=b+ξi。

(3)

對于最小二乘支持向量機問題，根據(jù)式(2)和式(3)定義Lagrange求解方程：

(4)

式中ai為Lagrange乘子。

參數(shù)a和b的最優(yōu)值可以通過KKT條件獲得：

(5)

由式(5)可得：

(6)

消去式(6)中的w和ξ，優(yōu)化問題轉(zhuǎn)化為求解如下方程：

(7)

(8)

核函數(shù)K(x,xi)為任意對稱函數(shù)，其應滿足Mercer條件，常用的核函數(shù)主要包括線性函數(shù)、徑向基函數(shù)及多項式函數(shù)等。最小二乘支持向量機的建模精度通常受2個因素的影響，即核函數(shù)參數(shù)與懲罰因子C的取值。本文采用粒子群算法確定核函數(shù)與懲罰因子的取值。

1.3 粒子群算法

粒子群算法(particle swarm optimization)是一種迭代優(yōu)化工具[17-18]。在粒子群算法中，把問題的解看作為搜索空間中的粒子。所有粒子都有被優(yōu)化函數(shù)決定的適應值，同時所有粒子通過速度決定其運動的方向和距離，其他粒子追隨當前的最優(yōu)粒子在解空間中搜索。粒子群算法首先產(chǎn)生一組初始化的隨機粒子，隨后采用迭代的方法尋求最優(yōu)解。在迭代過程中，粒子通過對2個極值的跟蹤更新自己。其中一個是粒子在每次搜索中的最優(yōu)解，稱為個體極值Pbest，另外一個是粒子群全部粒子在每次搜索中的最優(yōu)解，稱為全局極值gbest。粒子群中第i個粒子在n維空間的位置可用xi=(xi1,xi2,…,xin)表示，其速度可用vi=(vi1,vi2,…,vin)表示，第i個粒子的個體極值可表示為Pbest=(Pi1,Pi2,…,Pin)，粒子群的全局極值表示為gbest=(g1,g2,…,gn)。在搜索到這2個極值后，用式(9)和(10)來更新粒子的速度和位置。

vi(k+1)=uvi(k)+c1rand1(Pbest-xi(k))+

c2rand2(gbest-xi(k));

(9)

xi(k+1)=xi(k)+vi(k+1)。

(10)

式(9)—(10)中：c1、c2為學習因子，其取值在(0，2)之間； rand1和rand2為隨機數(shù)，取值在(0，1)之間；u表示動量系數(shù)，其值隨迭代改變。

1.4 最小二乘支持向量機和粒子群算法的位移反分析步驟

1)首先根據(jù)工程實際資料，確定待反分析的巖體力學參數(shù)及其取值范圍，采用正交設(shè)計的方法確定計算方案；

2)建立有限元模型對全部計算方案進行計算，得到各監(jiān)測點對應的位移值，將待反分析的巖體力學參數(shù)作為輸入?yún)?shù)xi，計算位移值作為輸出參數(shù)yi，形成學習的樣本；

3)初始化粒子群算法，主要包括粒子群的規(guī)模、各粒子的權(quán)重因子、計算迭代次數(shù)、產(chǎn)生的隨機粒子群向量及各粒子向量對應的最小二乘支持向量機懲罰因子C與核函數(shù)的參數(shù)σ2，把學習樣本集作為訓練樣本和檢驗樣本，將各粒子的個體極值設(shè)為當前位置，代入最小二乘支持向量機進行訓練并獲取對應的預測位移值；

4)計算各粒子相應的預測值與真實值的平均相對誤差，并將其作為各粒子的適應值。然后進行迭代計算，并更新各粒子的速度和位置參數(shù)，同時記憶個體與全局所對應的最優(yōu)適應值，滿足初始設(shè)置的最大迭代次數(shù)時計算終止，最后記憶最優(yōu)的懲罰因子和核函數(shù)參數(shù)(C，σ2)；

5)將上一步粒子群算法計算得到的最優(yōu)參數(shù)(C，σ2)代入最小二乘支持向量機模型，建立待反分析巖體力學參數(shù)與位移值之間的非線性映射關(guān)系；

6)利用待反分析巖體力學參數(shù)與位移值之間的非線性映射關(guān)系代替正分析中的有限元計算，將位移反分析的目標函數(shù)值作為每個粒子的適應值，利用粒子群算法計算與位移實測值最接近的待反分析巖體力學參數(shù)。

2 工程應用

2.1 工程概況[19]

厄瓜多爾CCS水電站位于納波省和蘇昆比奧斯省境內(nèi)的Coca河下游，為引水式電站，發(fā)電廠房為地下式，總裝機容量為1 500 MW。地下廠房洞室群由主廠房、主變室、施工排水洞及引水壓力管道等組成。主變室和主廠房2大洞室平行布置，洞室走向為NW45°，洞室斷面為圓拱直墻型，其中主廠房開挖尺寸為212.8 m(長)×27.5 m(寬)×46.8 m(高)，主變室開挖尺寸為192.0 m(長)×17.0 m (寬)×34.0 m(高)。

2.2 工程地質(zhì)條件

CCS水電站地下廠房埋深為300～400 m，區(qū)內(nèi)未見較大規(guī)模的斷層，但局部可見小斷層及節(jié)理、裂隙密集帶。地層巖性主要為侏羅系-白堊系迷撒華林地層(J-Km)紫紅色、青灰色及淺紅色的火山凝灰?guī)r、火山角礫巖及流紋巖等，巖體結(jié)構(gòu)以整體塊狀為主，局部呈塊狀、次塊狀。廠房區(qū)地下水以基巖裂隙水為主。廠房區(qū)內(nèi)受構(gòu)造應力影響不明顯，為中—低應力值水平。

2.3 地下廠房施工方案

主廠房和主變室施工開挖采用鉆爆法，由于主廠房和主變室的斷面大，無法采用全斷面開挖，故采用分層開挖的方法，其中主廠房分6層開挖，主變室分4層開挖，采用開挖一層支護一層的方法，支護方式為錨桿、掛網(wǎng)及噴混凝土。

2.4 有限元模型建立

地下洞室的數(shù)值模擬分析主要有二維模擬分析和三維模擬分析，二維平面分析的優(yōu)點是可快速建模、計算速度快，對于CCS地下廠房，其縱向尺寸遠大于橫向尺寸，故可采用二維平面模型。地下廠房洞室群中除了主廠房和主變室之外，還有母線洞、尾水洞等洞室，但是在地下廠房基本開挖完成后才施工的，且其尺寸遠小于地下廠房尺寸，所以對地下廠房影響較小，因此在本文的計算中不作考慮。

采用有限元程序Phase 2，根據(jù)地質(zhì)概化模型建立8#機組剖面數(shù)值模型，計算范圍如下：X方向為垂直主廠房和主變室軸線方向，以8#機組中心向上游方向延伸175 m、向下游方向延伸200 m，以指向下游方向為正；Y方向為豎直方向，以主廠房底板向下延伸90 m、向上延伸至地表面，以向上為正。計算范圍大于5倍洞室，根據(jù)彈性力學理論，能滿足計算邊界的要求。模型考慮了對主廠房和主變室圍巖穩(wěn)定性有影響的F10、F11、F13、F14等小規(guī)模斷層及地下水位。在數(shù)值計算過程中，斷層采用接觸單元模擬，斷層參數(shù)未參與反分析，主要是因為斷層開挖揭露后，對斷層帶物質(zhì)進行了現(xiàn)場及室內(nèi)試驗，得出了其基本參數(shù)，數(shù)值計算中采用的斷層基本參數(shù)如下： 切向剛度，3 500 MPa/m； 法向剛度，9 000 MPa/m； 黏聚力，0.10 MPa； 內(nèi)摩擦角，25°。廠房區(qū)8#機組有限元模型共劃分了6 391個節(jié)點、11 588個單元，網(wǎng)格劃分后的有限元計算模型如圖1所示。模型邊界條件如下： 兩側(cè)為X向位移約束，底部為固定約束，頂部為自由邊界。

圖1 8#機組剖面有限元模型

2.5 位移增量

在CCS水電站主廠房和主變室的施工開挖過程中，為監(jiān)測圍巖的變形，選取了多個斷面布置多點位移計，以8#機組斷面為例，在洞室的項拱、拱角和邊墻位置共布置了12條多點位移計測線，如圖2所示。由于多點位移計是在洞室開挖到相應高程后安裝的，因此其所測的位移值是該處多點位移計安裝后圍巖發(fā)生的位移，而測前發(fā)生的位移已被損失掉[20-21]，關(guān)于測前位移量損失的大小目前尚無好的確定方法，為避免反分析的隨意性，本文采用主廠房和主變室第2層開挖完成與第1層開挖完成后的頂拱中心的多點位移計實測位移增量進行反分析計算。多點位移計孔口實測BX1-29增量位移值為0.48 mm，BX1-34增量位移值為0.35 mm。

圖2 主廠房和主變室多點位移計布置(單位： m)

Fig. 2 Layout of multi-point extensometers of main powerhouse and main converter room (unit: m)

2.6 圍巖參數(shù)的取值范圍

根據(jù)CCS水電站地下廠房勘察設(shè)計資料及地下廠房圍巖的實際情況，待反分析的圍巖參數(shù)取值范圍如下： 彈性模量E∈[9，25]，GPa； 泊松比μ∈[0.20，0.28]； 黏聚力c∈[0.60，2.20]，MPa； 內(nèi)摩擦角φ∈[40，52]，(°)； 側(cè)壓力系數(shù)kx∈[0.30，0.70]、kz∈[1.30，1.70]。

2.7 反分析樣本設(shè)計

理論上，反分析參數(shù)的組合有無窮多個，如果對每組參數(shù)都進行計算，實施困難，且沒有必要，因此本文采用正交試驗的設(shè)計方法，選取一定有代表性的組合進行計算。正交試驗方法是根據(jù)數(shù)理統(tǒng)計學與正交性原理，從大量試驗點中選取合適的有代表性的點，再按照“正交表”安排試驗組合方案，由于正交表具有“均衡分散性”和“整齊可比性”的構(gòu)造原則，因此，按照此方法設(shè)計的試驗次數(shù)少，并且能反映客觀事物的變化規(guī)律。正交表是正交試驗設(shè)計的關(guān)鍵，它必須滿足以下2個條件： 每一列(因素)的不同水平在試驗中出現(xiàn)的次數(shù)相同，以保證其均勻性； 任意2列(因素)的不同水平組合組成的數(shù)對在試驗中出現(xiàn)的次數(shù)相同，以保證試驗點分布的均勻性。具體到本文的計算方案組合，圍巖參數(shù)共有6個，在每個參數(shù)的取值范圍內(nèi)均勻選取5個值，各因素的取值范圍并不是隨意選取的，其中巖石(體)的力學參數(shù)是根據(jù)現(xiàn)場和室內(nèi)試驗得到的范圍值，而地應力的側(cè)壓力系數(shù)是根據(jù)前期的地應力測試得到的范圍值，如kx的范圍值為0.3～0.7，這里將其分為5個水平，即0.3、0.4、0.5、0.6、0.7,以0.1為等級，能滿足精度的要求。由于最小二乘支持向量機在處理小樣本數(shù)據(jù)具有優(yōu)勢，25個不同圍巖參數(shù)組合能代表圍巖情況，滿足參數(shù)反分析的要求，不同計算方案下的參數(shù)組合及有限元計算所得的位移增量見表1。

2.8 反分析結(jié)果及檢驗

在有限元計算中，輸入不同的圍巖參數(shù)，即可得出相應的位移值，即圍巖力學參數(shù)與位移增量之間已有明確的映射關(guān)系，但實際的位移監(jiān)測值與位移計算值并不對應，所以實際上的圍巖參數(shù)是未知的，因此需要通過建立圍巖力學參數(shù)與位移增量之間的模糊映射關(guān)系，根據(jù)實際的監(jiān)測位移值采用逼近的方法來反分析圍巖參數(shù)。

本文首先根據(jù)反分析的參數(shù)樣本，進行有限元計算，各參數(shù)及位移增量值計算結(jié)果見表1。再根據(jù)最小二乘支持向量機和粒子群算法位移反分析理論和步驟，采用MATLAB編制相應的計算程序。最后將表1中的數(shù)據(jù)作為學習樣本，建立圍巖力學參數(shù)與洞室頂拱位移增量之間的非線性映射關(guān)系，按照1.3節(jié)的方法，經(jīng)反分析得到圍巖的力學參數(shù)反分析值如下： 彈性模量E=14.15 GPa，泊松比μ=0.23，黏聚力c=1.46 MPa，內(nèi)摩擦角φ=49°，側(cè)壓力系數(shù)kx=0.55、kz=1.45。為檢驗反分析所得圍巖力學參數(shù)的合理性，利用上述參數(shù)進行正分析有限元計算，得到主廠房和主變室分層開挖各多點位移計測點位移增量的計算值，選取主廠房第Ⅵ層與第Ⅰ層開挖、主變室第4層與第1層開挖的計算值與實測位移值增量進行對比，結(jié)果見圖3和圖4?？梢钥闯?，計算值和實測值相對誤差較小，最大誤差不超過10%，如BX1-34測點的誤差為6.5%，BX1-35測點的誤差為2.4%。在實際工程中，受監(jiān)測儀器本身精度、施工方法、鄰近爆破等的影響，監(jiān)測結(jié)果會有一定的誤差，但結(jié)果誤差在10%以內(nèi)是可以接受的，這也說明了反分析圍巖參數(shù)的合理性。

表1 圍巖力學參數(shù)及有限元計算結(jié)果

圖3 主廠房多點位移計位移增量計算值與實測值對比

Fig. 3 Comparison of displacement increment between calculation results and monitoring values of main powerhouse

圖4 主變室多點位移計位移增量計算值與實測值對比

Fig. 4 Comparison of displacement increment between calculation results and monitoring values of main converter room

3 結(jié)論與討論

1)最小二乘支持向量機可有效降低計算的復雜性，加快求解速度，增強抗干擾能力，具有良好的優(yōu)化性能，在處理小樣本、非線性問題上有較大的優(yōu)勢。粒子群算法作為一種新型仿生進化算法，有收斂速度快、參數(shù)較少易確定等優(yōu)點。最小二乘支持向量機和粒子群算法的有機結(jié)合為巖土工程位移反分析提供了一種新的方法。

2)本文以正交設(shè)計和有限元計算為基礎(chǔ)，采用最小二乘支持向量機和粒子群算法相結(jié)合的方法，建立CCS水電站地下廠房洞室頂拱實測位移與圍巖力學參數(shù)之間的非線性映射關(guān)系，成功地反分析出圍巖彈性模量、泊松比、黏聚力、內(nèi)摩擦角和側(cè)壓力系數(shù)。

3)多點位移計一般都是在洞室開挖后進行安裝，因此無法測量到測前洞室發(fā)生的位移，由于確定測前位移損失較為困難，本文采用主廠房和主變室第2層和第1層開挖所產(chǎn)生的位移增量進行反分析，簡化了反分析過程，有效提高了圍巖參數(shù)反分析的精度。

4)對CCS水電站地下廠房反分析所得圍巖參數(shù)進行檢驗，有限元計算結(jié)果表明： 主廠房第Ⅳ層與第Ⅰ層開挖和主變室第4層與第1層開挖所產(chǎn)生的位移增量計算值與多點位移計實測值吻合較好，最大誤差不超過10%，說明了反分析參數(shù)的合理性，同時也說明了采用基于最小二乘支持向量機和粒子群算法的位移反分析方法在工程上是可行的，且效果較為理想。

5)在目前的地下工程巖體參數(shù)反分析中采用的多是位移反分析，實際上地下洞室開挖后，巖體的應力、位移、塑性區(qū)等均有一定程度的改變，在以后的研究中，可考慮采用應力、位移及塑性區(qū)等多元信息的綜合分析，以便相互驗證，提高分析精度。