基于DE-FEM 隧道反分析可視化平臺(tái)研究及應(yīng)用

姜諳男，李 鵬，唐述林，王軍祥

（1. 大連海事大學(xué) 道路與橋梁研究所，遼寧 大連 116026；2. 中鐵21 局集團(tuán)有限公司，蘭州 730000）

1 引 言

由于地質(zhì)體環(huán)境的復(fù)雜性和不確定性，隧道預(yù)設(shè)計(jì)施工方案不可避免地存在盲目性。利用施工觀測(cè)信息，對(duì)隧道圍巖性質(zhì)和狀態(tài)進(jìn)行反饋分析，是工程界提倡的新奧法施工、信息化施工的精髓。隧道反饋分析包括定性經(jīng)驗(yàn)分析和定量計(jì)算分析，后者主要是參數(shù)反分析，即從力學(xué)模型層次上進(jìn)行定量識(shí)別分析，具有重要的學(xué)術(shù)價(jià)值和工程意義[1]。

巖土反分析的基本思想最先由Kavangh、Gioda和Maier 等提出[1]。日本學(xué)者櫻井春輔進(jìn)行了巖體彈性模量及初始地應(yīng)力的線彈性有限元位移反分析研究[2]。此后國(guó)內(nèi)外許多學(xué)者在此方面進(jìn)行了大量的研究工作，發(fā)展了彈塑性、黏彈性、黏塑性等非線性反分析，以及有關(guān)量測(cè)誤差處理、優(yōu)化、校驗(yàn)等技術(shù)[3－5]。馮夏庭[6]將智能科學(xué)引入巖石力學(xué)與工程反分析取得良好效果。近年來國(guó)內(nèi)也出現(xiàn)了用于隧道工程的反分析軟件，例如同濟(jì)曙光、BMP 反分析軟件等[7]。

盡管反分析已經(jīng)取得較多學(xué)術(shù)理論成果，但也應(yīng)看到，實(shí)際工程中反分析的應(yīng)用仍然不盡如人 意[8]。制約反分析應(yīng)用的主要原因是：（1）巖土工程的非線性特點(diǎn)使傳統(tǒng)優(yōu)化方法容易限于局部最小，而且線彈性簡(jiǎn)化條件與實(shí)際相差太大；（2）傳統(tǒng)反分析數(shù)據(jù)抽象，對(duì)專業(yè)知識(shí)和理論有較高要求，現(xiàn)場(chǎng)人員不易使用。因此，尋求魯棒性強(qiáng)的優(yōu)化算法并開發(fā)界面友好的反分析可視化平臺(tái)對(duì)于促進(jìn)反分析工程應(yīng)用，提高隧道安全保障具有重要的意義。

2 反分析原理及DE-FEM 步驟

2.1 反分析原理

反分析是指利用現(xiàn)場(chǎng)觀測(cè)物理量反算出圍巖初始地應(yīng)力、力學(xué)參數(shù)等。反分析分為解析法和數(shù)值法，數(shù)值法又分正反分析、逆反分析和圖譜法3種。其中正反分析法可以利用現(xiàn)有正演計(jì)算程序，具有很強(qiáng)的適應(yīng)性。其原理步驟如下：選取一組參數(shù)計(jì)算得到圍巖變形，如果計(jì)算變形與監(jiān)測(cè)變形數(shù)據(jù)相差較大，則調(diào)整輸入?yún)?shù)重新計(jì)算，直到計(jì)算與監(jiān)測(cè)數(shù)據(jù)相差很小，這時(shí)對(duì)應(yīng)參數(shù)即為識(shí)別的參數(shù)。由此可見，參數(shù)反分析本質(zhì)上是優(yōu)化問題。根據(jù)隧道圍巖參數(shù)的特定物理意義設(shè)定上下限，如果區(qū)域內(nèi)有m 個(gè)觀測(cè)值，那么有約束的優(yōu)化問題為

2.2 基于DE-FEM 的反分析過程

2.2.1 差異進(jìn)化的優(yōu)化步驟

目前以遺傳算法、粒子群算法和模擬退火算法等為代表的進(jìn)化優(yōu)化算法，以其良好的魯棒性、全局尋優(yōu)能力和不需要對(duì)目標(biāo)函數(shù)求導(dǎo)數(shù)等優(yōu)點(diǎn)顯示出強(qiáng)大生命力。其中差異進(jìn)化算法（Difference Evolutionary，DE）由美國(guó)加州大學(xué)伯克利大學(xué)的Kenneth Price 和Rainer Storm 提出的。該算法雖源于遺傳算法，由于不進(jìn)行編碼和解碼操作，使用上大為簡(jiǎn)化。同時(shí)它對(duì)初始值無要求，收斂速度快，對(duì)各種非線性函數(shù)適應(yīng)性強(qiáng)[9]。

式（1）的Yi是以圍巖參數(shù)xi（x1, x2,…, xn）為自變量的復(fù)雜非線性隱式函數(shù)，引進(jìn)差異進(jìn)化算法實(shí)現(xiàn)對(duì)參數(shù)xi全局優(yōu)化。假設(shè)需要識(shí)別n 個(gè)巖石力學(xué)參數(shù)，則首先產(chǎn)生維數(shù)為n 的Np個(gè)向量，Np稱為種群規(guī)模，每個(gè)向量也就是一組潛在的解，稱為個(gè)體。對(duì)每一個(gè)個(gè)體向量按式（1）進(jìn)行計(jì)算，以計(jì)算結(jié)果作為評(píng)價(jià)值。根據(jù)評(píng)價(jià)值大小按照DE 規(guī)則迭代，逐漸接近最優(yōu)解。將有限元算法嵌入到差異進(jìn)化算法中，形成DE-FEM 耦合的反分析計(jì)算流程見圖1。如圖所示，DE-FEM 的計(jì)算過程如下：

（1）產(chǎn)生初始種群。選擇需要反分析的參數(shù)（自變量）個(gè)數(shù)為n，首先隨機(jī)產(chǎn)生滿足自變量上、下界約束的NP個(gè)n 維向量，公式如下：

圖1 巖土反分析DE-FEM 計(jì)算流程 Fig.1 The DE-FEM calculation process of back analysis

（2）變異操作?？s放種群中任意2 個(gè)向量個(gè)體之間的差值并疊加到種群中的第3 個(gè)向量個(gè)體上，形成新的變異變量。第G +1 代第i 個(gè)變異向量的第j 個(gè)分量為

式中：下標(biāo)r1，r2，r3為1～NP中的隨機(jī)整數(shù)且互不相同；F 為縮放因子，用來調(diào)節(jié)向量差異的步長(zhǎng)幅值，在0～2 內(nèi)取值。式（3）是差異進(jìn)化算法的基本變異模式，被稱作DE/rand/1 模式。差異進(jìn)化還有其他模式，如DE/best/1、DE/best/2、DE/rand/2等[9]。

（3）交叉操作。將目標(biāo)向量xi(G)與變異向量vi(G +1)按照如下規(guī)則雜交，生成新的試樣向量ui(G +1)，其第j 個(gè)分量表示為

式中：rj∈[0，1]，為與向量第j 個(gè)分量對(duì)應(yīng)的隨機(jī)數(shù)；CR∈[0，1]，為雜交概率常數(shù)；rni為在1，2，…，n 中隨機(jī)挑選一個(gè)整數(shù)，以確保變異向量Vi(G +1)中，至少有一個(gè)分量被試樣向量ui(G +1)采用。

（4）選擇。采用貪婪搜索方法進(jìn)行選擇操作。將試樣向量ui(G +1)和目標(biāo)向量xi(G)比較，調(diào)用有限元程序計(jì)算，采用式（1）的目標(biāo)函數(shù)作為評(píng)價(jià)函數(shù)，如果ui(G+1)對(duì)應(yīng)較小的目標(biāo)函數(shù)值，則選擇ui(G +1)；反之，如果xi(G)對(duì)應(yīng)較小的目標(biāo)函數(shù)值，則保留xi(G)。

（5）循環(huán)迭代。重復(fù)（2）～（4）步的計(jì)算，直到i 從1～Np、j 從1～n 循環(huán)完成，即完成種群的一個(gè)迭代。循環(huán)直到達(dá)到最大迭代步數(shù)或者適應(yīng)值小于設(shè)定值，則中止迭代，輸出當(dāng)前對(duì)應(yīng)的向量即為所識(shí)別的參數(shù)。

2.2.2 應(yīng)力回映彈塑性算法

考慮到線彈性不能反映巖土體的材料非線性特點(diǎn)，在式（1）圍巖位移Yi的計(jì)算中采用彈塑性有限元方法求解。塑性有限元應(yīng)力和應(yīng)變之間的非線性關(guān)系導(dǎo)致平衡方程是應(yīng)變（也是節(jié)點(diǎn)位移）的非線性方程，變形歷史取決于材料的非線性本構(gòu)關(guān)系，因此，采用增量分析方法跟蹤位移、應(yīng)變和外部作用引起的應(yīng)力的發(fā)展[10]。采用Newton Raphson 方法求解非線性方程組，采用應(yīng)力回映算法實(shí)現(xiàn)彈性預(yù)測(cè)-塑性修正過程。彈性預(yù)測(cè)-塑性修正的步驟如下（見圖2）：

圖2 應(yīng)力回映算法 Fig.2 The algorithm of stress return mapping

（1）彈性預(yù)測(cè)

當(dāng)塑性因子 0γΔ = 時(shí)，

式中： γ ( t)為塑性因子； α ( t)為內(nèi)部變量；A(t)為強(qiáng)化參數(shù)；Φ (t )為屈服函數(shù)；為密度；ψ ( t)為應(yīng)變能函數(shù)。由式（5）～（7）得出的結(jié)果是彈性的應(yīng)力狀態(tài)，但還不是最終本構(gòu)方程的解，如果滿足式（8），即是在彈性范圍內(nèi)或是屈服面上，則所得到的解是方程的最后解：

（2）塑性修正

如果不滿足式（8），則利用下面的式子進(jìn)行計(jì)算，即返回映射方程。

式中：H 為定義硬化變量變化的綜合硬化模量；N為流動(dòng)矢量函數(shù)；φ 為標(biāo)量屈服函數(shù)。此時(shí)也要滿足 0γΔ ＞ 的要求，最終得出解。

3 基于VTK 可視化平臺(tái)開發(fā)

基于DE-FEM 算法，采用微軟VC++研制了隧道反分析可視化軟件系統(tǒng)（Engineering Emulation Optimization System，EEOS）。該平臺(tái)分為4 個(gè)模塊，即前處理模塊、FEM 計(jì)算模塊、DE 算法模塊和后處理模塊，各模塊按照功能要求采用面向?qū)ο蟮木幊碳夹g(shù)封裝設(shè)計(jì)了相應(yīng)的類，見圖3。前處理模塊包括三角單元剖分類Delaunaymesh，柵格法四邊形單元剖分類Quadrilateralmesh 和Ansys 網(wǎng)格輸入類，實(shí)現(xiàn)基本模型圖形的網(wǎng)格剖分和復(fù)雜模型網(wǎng)格的導(dǎo)入。FEM 計(jì)算模塊包括區(qū)域類（Domain）、單元（Element）類、節(jié)點(diǎn)（Node）類、邊界（Boundary）類、材料（Material）類、時(shí)步（Time step）類和數(shù)學(xué)求解（Math）類等，實(shí)現(xiàn)彈性、摩爾-庫(kù)侖彈塑性的有限元計(jì)算。DE 算法模塊包括種群（Population）類、進(jìn)化（Evolution）類、適應(yīng)函數(shù)（Fitness）類和數(shù)據(jù)輸入輸出（Inoutdata）類。其中適應(yīng)函數(shù)類之中調(diào)用FEM 計(jì)算模塊計(jì)算位移，差異迭代操作產(chǎn)生的新參數(shù)變量輸入FEM 計(jì)算模塊的材料類，以此實(shí)現(xiàn)DE-FEM 交互耦合。

圖3 反分析可視化系統(tǒng)EEOS 的程序模塊 Fig.3 The program module of back analysis visualizing system EEOS

能否將抽象的結(jié)果數(shù)據(jù)以可視化形式呈現(xiàn)給現(xiàn)場(chǎng)用戶，在很大程度上決定了反分析能否在工程中得到推廣?？茖W(xué)計(jì)算可視化已成為學(xué)術(shù)和工程界的熱點(diǎn)問題[11－13]。直接進(jìn)行可視化底層代碼編寫需要消耗大量時(shí)間和精力，質(zhì)量也不容易保證。VTK (Visualization Toolkit)是一個(gè)開放源碼、跨平臺(tái)、支持平行處理的圖形應(yīng)用函數(shù)庫(kù)，以此編寫可視化程序其質(zhì)量和效率有顯著提高，在國(guó)外的大型研究機(jī)構(gòu)如Sandia, Los Alamosn 及Livermore 國(guó)家實(shí)驗(yàn)室等都得到應(yīng)用。

圖4 EEOS 反分析可視化系統(tǒng)界面 Fig.4 The interface of back analysis visualizing system EEOS

本文采用VTK 來進(jìn)行處理模塊編寫。后處理模塊包括等值線（Contour）類、數(shù)據(jù)輸出（Dataout）類和三維建模（3D model）類。利用VTK 的vtlPlane，vtkCutter，vtkProbeFilter 類實(shí)現(xiàn)云圖，利用VTK 的vtkPlane，vtkContourFilter，vtkProbeFilter 類實(shí)現(xiàn)等值線圖。采用輪廓線建模，即在二維數(shù)據(jù)切片(DataSlice)中逐一提取閉合的等值線，然后將相鄰切片的等值線相連接，形成曲面網(wǎng)格逼近等值表面，這樣便可實(shí)現(xiàn)三維模型。該平臺(tái)界面友好易用，采用了分割窗口的界面技術(shù)，界面窗口左側(cè)為樹狀和按鈕式控制欄，右側(cè)為顯示窗體。底部為狀態(tài)欄分析過程狀態(tài)信息，如圖4 所示。

4 工程應(yīng)用

4.1 工程概況

陜西省延安市的某隧道全長(zhǎng)2 165 m，范圍內(nèi)地表之處三面高，中部低，地層主要為第四系上更新統(tǒng)風(fēng)積砂質(zhì)黃土、中更新統(tǒng)黏質(zhì)黃土和下伏的侏羅系砂巖夾頁(yè)巖。地表第四系覆蓋層較厚，地勢(shì)陡峻，地形起伏較大。隧道洞身通過主要地層為風(fēng)化砂巖夾頁(yè)巖，巖質(zhì)較差，產(chǎn)狀近于水平，節(jié)理較發(fā)育，成份以石英、長(zhǎng)石、云母為主。隧道設(shè)計(jì)為城門洞型，高為8 m，寬為7.5 m。隧道施工過程進(jìn)行了圍巖的監(jiān)控量測(cè)，工程地質(zhì)剖面和監(jiān)測(cè)斷面收斂布置見圖5，其中典型監(jiān)測(cè)斷面KD+500、KD+550和KD+600 對(duì)應(yīng)的埋深分別為76、110、146 m，各斷面測(cè)線的監(jiān)測(cè)值見表1。

圖5 隧道的地質(zhì)剖面和監(jiān)測(cè)斷面 Fig.5 Geologic section and monitoring sections of tunnel

表1 主要斷面監(jiān)測(cè)數(shù)據(jù) Table 1 The monitoring data of main sections

4.2 參數(shù)反分析

調(diào)用本文EEOS 平臺(tái)分別對(duì)KD+500、KD+550和KD+600 斷面相關(guān)參數(shù)進(jìn)行反分析，采用摩爾-庫(kù)侖的本構(gòu)模型，對(duì)彈性模量E、黏聚力c 和內(nèi)摩擦角φ 進(jìn)行識(shí)別。圍巖區(qū)域之外的巖層對(duì)隧道位移和穩(wěn)定性影響不大，各斷面采用均質(zhì)模型。初始垂直應(yīng)力按照自重應(yīng)力場(chǎng)考慮，模型區(qū)域左右邊界分別距離隧道中心為22 m，上邊界為自由地表，下邊界距離隧道為8 m。泊松比v 按照先驗(yàn)信息取為0.25。根據(jù)勘測(cè)資料，待反演參數(shù)圍巖彈性模量E、黏聚力c 和內(nèi)摩擦角φ 的取值范圍分別為：100 MPa ≤E≤800 MPa、0.1 MPa≤c≤1.6 MPa、15°≤ φ ≤45°。3 個(gè)斷面反分析參數(shù)的計(jì)算位移結(jié)果與監(jiān)測(cè)位移對(duì)比見表2。計(jì)算位移與監(jiān)測(cè)位移較為接近，最大相對(duì)誤差為4.76%，反分析精度令人滿意。

表2 各斷面監(jiān)測(cè)信息與反演結(jié)果 Table 2 The monitoring and back analysis results of each section

4.3 DE 參數(shù)對(duì)反分析影響

在EEOS 差異進(jìn)化算法進(jìn)行參數(shù)反分析過程中，CR 和F 取值對(duì)差異進(jìn)化算法收斂有不同的影響，一般說來，CR 和F 取值都宜大于0.5。以KD+600斷面為例，變化差異進(jìn)化算法的參數(shù)，縮放因子F為0.7，雜交概率常數(shù)CR 分別為0.9、0.8 和0.7，研究?jī)?yōu)化算法的收斂性能，所對(duì)應(yīng)的收斂曲線見圖6。這3 種情況下CR=0.9 所對(duì)應(yīng)的收斂性能最好，在第9 代適應(yīng)值就可收斂；CR=0.7 收斂性最差，在第15 代才收斂。總體上這些參數(shù)都能保證算法收斂，較快地搜索到力學(xué)參數(shù)值，體現(xiàn)了差異進(jìn)化算法具有的良好的尋優(yōu)能力。

圖6 不同參數(shù)下DE 收斂曲線（KD+600） Fig.6 The DE convergence curves with different parameters(KD+600)

進(jìn)一步將上述問題的差異進(jìn)化與遺傳算法搜索過程(GA)進(jìn)行比較（種群規(guī)模和迭代次數(shù)相同），收斂曲線的對(duì)比見圖7。由圖可見，差異進(jìn)化的收斂速度和參數(shù)識(shí)別精度均優(yōu)于遺傳算法。基于EEOS 反演的參數(shù)計(jì)算KD+600 斷面，水平位移等值線云圖見圖8，可視化效果良好，工程人員可直觀地獲得該斷面水平位移分布規(guī)律。

圖7 GA 與DE 收斂曲線對(duì)比（KD+600） Fig.7 The comparisons between the convergence curves of GA and DE

圖8 KD+600 斷面水平位移云圖 Fig.8 The horizontal displacement contour map of KD+600 section

BMP 采用彈性的假設(shè)條件，同濟(jì)曙光軟件采用遺傳算法進(jìn)行參數(shù)反分析。本文在吸收上述反分析軟件優(yōu)點(diǎn)的基礎(chǔ)上，引入優(yōu)化能力更好的差異進(jìn)化算法、收斂性更好的彈塑性應(yīng)力回映算法，而且采用VTK 類庫(kù)實(shí)現(xiàn)計(jì)算結(jié)果可視化，使得反分析程序具有更好的易用性和可實(shí)現(xiàn)性。該軟件上手簡(jiǎn)單，現(xiàn)場(chǎng)技術(shù)人員容易掌握使用，采用一般PC 機(jī)求解一般規(guī)模隧道工程（1 000 個(gè)單元）反分析問題，需要5～6 h，可以滿足工程現(xiàn)場(chǎng)實(shí)時(shí)反饋的要求。

5 結(jié) 論

（1）差異進(jìn)化算法是具有全局優(yōu)化特性的計(jì)算智能算法，解決傳統(tǒng)優(yōu)化求導(dǎo)數(shù)困難以及容易限于局部最小化的問題。與線彈性模型相比，本文應(yīng)力回映的彈塑性有限元方法能夠反映巖土材料非線性，收斂快速?；贒E-FEM 隧道圍巖反分析能獲得較高的反分析精度。

（2）科學(xué)計(jì)算可視化使結(jié)果數(shù)據(jù)直觀表達(dá)，是計(jì)算方法在工程現(xiàn)場(chǎng)推廣應(yīng)用的重要決定因素。本文采用主流的VTK 可視化開發(fā)技術(shù)，實(shí)現(xiàn)了高質(zhì)量反分析可視化平臺(tái)，解決了一般反分析的輸入、輸出數(shù)據(jù)繁雜、抽象的問題。VTK 可視化開發(fā)方法加快了軟件的設(shè)計(jì)和實(shí)現(xiàn)，降低了系統(tǒng)開發(fā)難度，減少了開發(fā)時(shí)間。

（3）初步工程應(yīng)用表明，本文方法和系統(tǒng)是成功的?？梢暬脚_(tái)界面友好、結(jié)果表達(dá)直觀，方便現(xiàn)場(chǎng)技術(shù)人員使用。本文方法和思路為類似的工程軟件開發(fā)提供了參考，對(duì)促進(jìn)隧道反分析技術(shù)推廣應(yīng)用、推進(jìn)國(guó)產(chǎn)軟件發(fā)展都具有積極的作用。

[1] 楊林德. 巖土工程問題的反演理論與工程實(shí)踐[M]. 北京: 科學(xué)出版社, 1996.

[2] SAKURAI S, TAKEUCHI K. Back analysis of measured displacements of tunnels[J]. Rock Mechanics and Rock Engineering, 1983, 16: 173－180.

[3] PIERPAOLO ORESTE. Back analysis techniques for the improvement of the understanding of rock in underground constructions[J]. Tunneling and Underground Space Technology, 2005, (20): 7－21.

[4] 文建華, 吳代華, 陳軍明, 等. 地下洞室黏彈性位移反分析模式分層運(yùn)算[J]. 巖土力學(xué), 2010, 31(3): 967－970. WEN Jian-hua, WU Dai-hua, CHEN Jun-ming, et al. Parameter inversion of viscoelastic cavern displacements based on hierarchical pattern search[J]. Rock and Soil Mechanics, 2010, 31(3): 967－970.

[5] 張志增, 李仲奎. 橫觀各向同性巖體中圓形巷道反分析的惟一性[J]. 巖土力學(xué), 2011, 32(7): 2066－2072. ZHANG Zhi-zeng, LI Zhong-kui, Uniqueness of displacement back analysis of a circular tunnel in transversely isotropic rock mass[J]. Rock and Soil Mechanics, 2011, 32(7): 2066－2072.

[6] 馮夏庭. 智能巖石力學(xué)導(dǎo)論[M]. 北京: 科學(xué)出版社, 2000.

[7] 李世輝, 隧道支護(hù)設(shè)計(jì)新論[M]. 北京: 科學(xué)出版社, 1999.

[8] 張魯渝, 歐陽(yáng)小秀, 鄭穎人. 國(guó)內(nèi)巖土邊坡穩(wěn)定分析軟件面臨的問題及幾點(diǎn)思考[J]. 巖石力學(xué)與工程學(xué)報(bào), 2003, 22(1): 166－169. ZHANG Lu-yu, OUYANG Xiao-xiu, ZHENG Ying-ren, Problems and thoughts of development of slope stability analysis software in China[J]. Chinese Journal of Rock Mechanics and Engineering, 2003, 22(1): 166－169.

[9] JARMO ILONEN, AMARAINEN JOWI-KRISTIAN, LAMPIEN J. Differential evolution training algorithm for feed-forward neural networks[J]. Neural Processing Letters, 2003, 17(1): 93－105.

[10] NETO EA DE SOUZA, PERI'C D, OWEN DRJ. Computational methods for plasticity: theory and applications[M]. [S. l.]: John Wilen & Sons, Inc., 2008.

[11] PUSCH R. Practical visualization of rock structure[J]. Engineering Geology, 1998, 49: 231－236.

[12] 柴賀軍, 黃地龍, 黃潤(rùn)秋, 等. 巖體結(jié)構(gòu)三維可視化模型研究進(jìn)展[J]. 地球科學(xué)進(jìn)展, 2001, 11(4): 55－59. CHAI He-jun, HUANG Di-long, HUANG Run-qiu, et al. New progress in the study of rock structure 3-D visualization model[J]. Advance in Earth Sciences, 2001, 11(4): 55－59.

[13] WILLIAM J S. The visualization toolkit users guide, version 4.0[M]. New York: Kitware Inc, 2001.