基于DE-FEM 隧道反分析可視化平臺研究及應用

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

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

1 引 言

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

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

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

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

2.1 反分析原理

反分析是指利用現場觀測物理量反算出圍巖初始地應力、力學參數等。反分析分為解析法和數值法，數值法又分正反分析、逆反分析和圖譜法3種。其中正反分析法可以利用現有正演計算程序，具有很強的適應性。其原理步驟如下：選取一組參數計算得到圍巖變形，如果計算變形與監測變形數據相差較大，則調整輸入參數重新計算，直到計算與監測數據相差很小，這時對應參數即為識別的參數。由此可見，參數反分析本質上是優化問題。根據隧道圍巖參數的特定物理意義設定上下限，如果區域內有m 個觀測值，那么有約束的優化問題為

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

2.2.1 差異進化的優化步驟

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

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

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

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

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

式中：下標r1，r2，r3為1～NP中的隨機整數且互不相同；F 為縮放因子，用來調節向量差異的步長幅值，在0～2 內取值。式（3）是差異進化算法的基本變異模式，被稱作DE/rand/1 模式。差異進化還有其他模式，如DE/best/1、DE/best/2、DE/rand/2等[9]。

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

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

（4）選擇。采用貪婪搜索方法進行選擇操作。將試樣向量ui(G +1)和目標向量xi(G)比較，調用有限元程序計算，采用式（1）的目標函數作為評價函數，如果ui(G+1)對應較小的目標函數值，則選擇ui(G +1)；反之，如果xi(G)對應較小的目標函數值，則保留xi(G)。

（5）循環迭代。重復（2）～（4）步的計算，直到i 從1～Np、j 從1～n 循環完成，即完成種群的一個迭代。循環直到達到最大迭代步數或者適應值小于設定值，則中止迭代，輸出當前對應的向量即為所識別的參數。

2.2.2 應力回映彈塑性算法

考慮到線彈性不能反映巖土體的材料非線性特點，在式（1）圍巖位移Yi的計算中采用彈塑性有限元方法求解。塑性有限元應力和應變之間的非線性關系導致平衡方程是應變（也是節點位移）的非線性方程，變形歷史取決于材料的非線性本構關系，因此，采用增量分析方法跟蹤位移、應變和外部作用引起的應力的發展[10]。采用Newton Raphson 方法求解非線性方程組，采用應力回映算法實現彈性預測-塑性修正過程。彈性預測-塑性修正的步驟如下（見圖2）：

圖2 應力回映算法 Fig.2 The algorithm of stress return mapping

（1）彈性預測

當塑性因子 0γΔ = 時，

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

（2）塑性修正

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

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

3 基于VTK 可視化平臺開發

基于DE-FEM 算法，采用微軟VC++研制了隧道反分析可視化軟件系統（Engineering Emulation Optimization System，EEOS）。該平臺分為4 個模塊，即前處理模塊、FEM 計算模塊、DE 算法模塊和后處理模塊，各模塊按照功能要求采用面向對象的編程技術封裝設計了相應的類，見圖3。前處理模塊包括三角單元剖分類Delaunaymesh，柵格法四邊形單元剖分類Quadrilateralmesh 和Ansys 網格輸入類，實現基本模型圖形的網格剖分和復雜模型網格的導入。FEM 計算模塊包括區域類（Domain）、單元（Element）類、節點（Node）類、邊界（Boundary）類、材料（Material）類、時步（Time step）類和數學求解（Math）類等，實現彈性、摩爾-庫侖彈塑性的有限元計算。DE 算法模塊包括種群（Population）類、進化（Evolution）類、適應函數（Fitness）類和數據輸入輸出（Inoutdata）類。其中適應函數類之中調用FEM 計算模塊計算位移，差異迭代操作產生的新參數變量輸入FEM 計算模塊的材料類，以此實現DE-FEM 交互耦合。

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

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

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

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

4 工程應用

4.1 工程概況

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

圖5 隧道的地質剖面和監測斷面 Fig.5 Geologic section and monitoring sections of tunnel

表1 主要斷面監測數據 Table 1 The monitoring data of main sections

4.2 參數反分析

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

表2 各斷面監測信息與反演結果 Table 2 The monitoring and back analysis results of each section

4.3 DE 參數對反分析影響

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

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

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

圖7 GA 與DE 收斂曲線對比（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 采用彈性的假設條件，同濟曙光軟件采用遺傳算法進行參數反分析。本文在吸收上述反分析軟件優點的基礎上，引入優化能力更好的差異進化算法、收斂性更好的彈塑性應力回映算法，而且采用VTK 類庫實現計算結果可視化，使得反分析程序具有更好的易用性和可實現性。該軟件上手簡單，現場技術人員容易掌握使用，采用一般PC 機求解一般規模隧道工程（1 000 個單元）反分析問題，需要5～6 h，可以滿足工程現場實時反饋的要求。

5 結 論

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

（2）科學計算可視化使結果數據直觀表達，是計算方法在工程現場推廣應用的重要決定因素。本文采用主流的VTK 可視化開發技術，實現了高質量反分析可視化平臺，解決了一般反分析的輸入、輸出數據繁雜、抽象的問題。VTK 可視化開發方法加快了軟件的設計和實現，降低了系統開發難度，減少了開發時間。

（3）初步工程應用表明，本文方法和系統是成功的?？梢暬脚_界面友好、結果表達直觀，方便現場技術人員使用。本文方法和思路為類似的工程軟件開發提供了參考，對促進隧道反分析技術推廣應用、推進國產軟件發展都具有積極的作用。

[1] 楊林德. 巖土工程問題的反演理論與工程實踐[M]. 北京: 科學出版社, 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] 文建華, 吳代華, 陳軍明, 等. 地下洞室黏彈性位移反分析模式分層運算[J]. 巖土力學, 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]. 巖土力學, 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] 馮夏庭. 智能巖石力學導論[M]. 北京: 科學出版社, 2000.

[7] 李世輝, 隧道支護設計新論[M]. 北京: 科學出版社, 1999.

[8] 張魯渝, 歐陽小秀, 鄭穎人. 國內巖土邊坡穩定分析軟件面臨的問題及幾點思考[J]. 巖石力學與工程學報, 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] 柴賀軍, 黃地龍, 黃潤秋, 等. 巖體結構三維可視化模型研究進展[J]. 地球科學進展, 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.