小樣本數(shù)據(jù)下三維地應力反演分析

陳正林， 何國志， 張劼超， 劉加柱， 王紅娟， 候圣均

(1.中電建路橋集團有限公司， 北京 100048； 2.北京科技大學土木與資源工程學院， 北京 100083)

對于隧道和地下工程，巖體的地應力場為受力狀態(tài)分析和支護加固提供數(shù)據(jù)支持，在了解工程區(qū)域原巖應力狀態(tài)的前提下，可以提前規(guī)劃工程設計和施工方案，有效地規(guī)避變形、坍塌等常見地質災害，減少人員傷亡和經(jīng)濟損失。

為了了解一個區(qū)域的地應力分布狀態(tài)，常采用地應力反演的方法，即通過區(qū)域內某幾個實測應力的數(shù)據(jù)和數(shù)學推演計算來推斷區(qū)域內三維地應力場?，F(xiàn)場測量方法一般采用水壓致裂法、應力解除法、聲發(fā)射法、超聲波譜法和發(fā)射型同位素法等[1-2]，本文研究現(xiàn)場測量采用空心包體應力解除法。關于獲取到實測數(shù)據(jù)后，采用的地應力反演方法，1983年，郭懷志等[3]和余云燕等[4]提出多元回歸分析法，同時考慮由埋深確定的巖體上覆重力和依據(jù)實際工程地質條件模擬的應力，在二者的基礎上構建了不同工況下應力場。方明禮等[5]大幅改進了地應力函數(shù)反分析法，提出用三維有限元反演地應力場方法來解決地應力計算。1999年尚岳全[6]提出直接邊界調整法，通過不斷調整計算區(qū)域邊界上的荷載大小去分析研究，使用數(shù)值方法來運算獲取應力場的分布變化規(guī)律[3-6]。21世紀初，隨著計算機的飛速發(fā)展，人工智能方法[7-8]逐漸被應用在地應力反演領域。

為了克服某些工程區(qū)域地質情況復雜，實測應力數(shù)據(jù)較少問題和減少現(xiàn)場測量的耗費，現(xiàn)減少測點，提出一種邊界荷載和人工神經(jīng)網(wǎng)絡的聯(lián)合反演方法，在小樣本下對區(qū)域內地應力狀態(tài)進行計算，綜合兩者的優(yōu)點，并對建(個)元高速五老峰隧道部分區(qū)域進行地應力反演，在測點較少的情況下，利用人工神經(jīng)網(wǎng)絡的數(shù)據(jù)處理和非線性優(yōu)勢模擬地應力分布狀態(tài)，希望可以給其他相關工程提供一定的理論依據(jù)和指導。

1 反演理論及實現(xiàn)

1.1 邊界荷載法

在研究礦區(qū)、隧道、邊坡、硐室等區(qū)域的構造應力場時，直接邊界調整法具有廣泛的應用性，可以很好地反演待演區(qū)域的初始地應力場和邊界條件。直接邊界調整法的具體步驟如下。

(1)根據(jù)待演區(qū)域的地勘報告、地形地貌特征、物探報告等地質信息，結合三維建模軟件建立包括實際測點在內的區(qū)域三維地質力學模型。

(2)選定一組邊界條件進行有限單元計算，監(jiān)測記錄測點的計算應力結果，與現(xiàn)場實測應力進行對比，分析出入原因，結合反演區(qū)域地質條件和模型本身，進行邊界條件的大小和作用方式的調整。

(3)對三維模型重新施加調整后的邊界條件和巖石力學參數(shù)，再次進行軟件計算，然后再進行對比，循環(huán)往復，直至得到一組最佳邊界條件使得測點的計算結果和實測結果相近[9]。

1.2 BP神經(jīng)網(wǎng)絡

反向傳播神經(jīng)網(wǎng)絡(BP)屬于前向型網(wǎng)絡結構的一種，全稱為“前向多層誤差逆向傳播學習算法”，下文簡稱為BP網(wǎng)絡。這種人工神經(jīng)網(wǎng)絡通過反復改變網(wǎng)絡拓撲結構中的各神經(jīng)元之間的連接權值和閾值，使網(wǎng)絡的輸出不斷逼近預設的期望輸出。相比于其他網(wǎng)絡結構，其最大的特點就是BP網(wǎng)絡是將網(wǎng)絡輸出與期望輸出之間的誤差通過網(wǎng)絡結構一層一層反向傳輸，各層神經(jīng)元結構通過返回的誤差值進行權值何閾值的調整，這種算法因此也被稱為反向學習算法。

BP網(wǎng)絡的結構和連接方式如圖1所示，主要結構由輸入、輸出和隱含層組成。每一層有大量的節(jié)點，相鄰層節(jié)點之間進行全連接，相隔層節(jié)點之間則無連接。網(wǎng)絡的學習方式也就是網(wǎng)絡對數(shù)據(jù)的處理方式，在BP網(wǎng)絡里，外界的數(shù)據(jù)按照從左至右的順序，即從輸入層經(jīng)過隱含層最后在輸出層得到網(wǎng)絡輸出[10-11]。

圖1 BP網(wǎng)絡結構示意圖Fig.1 Schematic diagram of BP network structure

1.3 地應力反演因素分析

一個地區(qū)的初始應力場的形成，往往經(jīng)歷了幾百甚至上千萬年的地質運動和外界因素的影響，且時刻在變化，不過對于某個地質年代來說，地應力場可以認為是穩(wěn)定的。地應力場在形成過程中，影響因素眾多，不同的區(qū)域，決定性因素和主導因素也不盡相同，研究人員也并沒有定論。但是經(jīng)過研究人員多年的現(xiàn)場考證和實驗理論研究，證明初始應力場主要受以下幾個因素影響。

(1)重力因素。巖體的自重應力是引起初始應力場中垂直應力的主要原因，區(qū)域任意一點的豎直應力大約等于單位面積上的上覆巖層的重量。

(2)地殼構造運動。地殼的構造運動是三維應力場中水平應力形成的主要原因，中國的大陸板塊由于受四周的板塊擠壓和約束，內部主要產(chǎn)生水平的擠壓應力。

(3)巖石力學性質。有學者研究證明，應力場的大小和形成與區(qū)域巖石的力學性質有一定關系，如區(qū)域巖層的彈性模量、泊松比、剛度等因素[12]。

(4)地形因素。數(shù)萬年的地球活動造成了現(xiàn)在的山川大河，地表的高低起伏，不同的形地貌，地形對地應力的影響并不確定，不同地形處不同深度的地應力狀況也并無規(guī)律可言。

(5)其他因素。除了以上對地應力影響較為明顯的因素外，溫度、滲流、冰川以及某些暫時不明的因素在某些特定情況下也能對區(qū)域內應力場產(chǎn)生一定的影響。

1.4 邊界荷載法-BP神經(jīng)網(wǎng)絡聯(lián)合反演

對于以上眾多地應力影響因素，可以用式(1)[13]表示：

σ=F(P,R,ΔG,ΔH,T,W,Q)

(1)

式(1)中：σ為區(qū)域某點地應力；F為一種非線性函數(shù)關系；P為測點的坐標位置，可在現(xiàn)場測量獲得；R為區(qū)域的巖石力學參數(shù)；ΔG為重力因素；ΔH為構造運動因素；T、W、Q為溫度、滲流和其他暫時不明因素。

通過式(1)，將地應力表示成重力、構造應力、溫度等因素的某種函數(shù)關系。在實際地應力反演過程中，測點的坐標位置是不可改變的，溫度、滲流和其他的不明因素也不考慮。所以式(1)可以進一步簡化，得到

σ=F(R,ΔG,ΔH)

(2)

在地應力場多元線性回歸反演理論中，初始地應力場可以簡化為各種因素產(chǎn)生的應力場的線性疊加，最終上述地應力函數(shù)可簡化為

σ=C0+C1σG+C2σH+C3σR+…+Cnσn

(3)

式(3)中：C0～Cn為回歸系數(shù)；σG為由在重力因素下產(chǎn)生的應力場；σH為在構造應力下產(chǎn)生的應力場；σR為由于巖石力學性質導致的應力場；σn為其他因素產(chǎn)生的應力場。

但是在實際上，一個區(qū)域內的地應力場是復雜多變的，尤其是在埋深較大或者地質情況較為復雜的地質區(qū)域，地應力場并不是簡單地疊加，而是具有明顯的非線性特征。因此引入BP網(wǎng)絡技術，如何利用BP網(wǎng)絡處理非線性問題的特性來解決地應力場的反演是需要解決的問題，但是BP網(wǎng)絡需要大量樣本訓練網(wǎng)絡和實際地應力測量數(shù)據(jù)稀少又構成了顯著的矛盾。

首先，地應力的數(shù)值反演基本上是在有限元軟件上進行，基本思路就是建立待反演區(qū)域的三維地質模型，賦予合適的參數(shù)和邊界條件，使實測點的計算數(shù)據(jù)和實際測量數(shù)據(jù)相近，這樣就可以得到一個區(qū)域的大致應力場。在多元線性回歸反演方法中，可以通過回歸后的方程推算邊界條件的大小和作用方式，但是在非線性反演中，并沒有一個確定的方程或者函數(shù)關系，因此就無法獲得準確的邊界條件。

在數(shù)據(jù)稀少的情況下，通常采用直接邊界調整法，通過試算的手法來慢慢調整邊界條件，但是這樣非常耗費時間和精力，因此結合BP網(wǎng)絡的非線性特征，提出以下思路。

首先建立待反演區(qū)域的三維地質模型，通過直接邊界調整法進行反演計算，在試算的基礎上得到大概的巖體參數(shù)和邊界條件范圍(實測點的地應力值被包含在內)。將試算出的巖體參數(shù)和邊界條件按照正交試驗構造原則，構造出若干組組合邊界條件。將組合邊界和巖石力學參數(shù)賦予模型進行“正算”，得到實測點的計算數(shù)據(jù)。這樣就可以得到若干組BP網(wǎng)絡的訓練學習樣本，以測點的地應力值作為輸入，巖石力學參數(shù)邊界條件作為輸出。構建BP網(wǎng)絡，帶入數(shù)據(jù)進行訓練學習，得到地應力值和邊界條件的非線性映射關系。將測點的實際測量數(shù)據(jù)帶入訓練完成的網(wǎng)絡，“反算”出最佳邊界條件，代入模型中，進行模擬計算驗證。具體流程如圖2所示。

圖2 地應力非線性反演流程Fig.2 Nonlinear inversion flow of ground stress

2 工程實例分析

2.1 五老峰隧道地應力實測

五老峰隧道位于建水至元陽段，隧道左線全長8 305 m，起點樁號Z3K22+670,終點樁號Z3K30+975;右線全長8 360 m，起點樁號K22+665,終點樁號K31+025。隧道一般埋深50～900 m，最小埋深25 m，最大埋深929 m，屬典型的深埋特長分離式隧道。本次現(xiàn)場測量采用新型數(shù)字化無線式瞬接續(xù)采型原位空心包體應變計[14-15](圖3)。由于實驗環(huán)境復雜，本次只進行一個點的地應力測試，具體計算過程限于篇幅原因，不作具體介紹，測量結果如表1所示。為了便于計算，將主應力進行應力分解，結果如表2所示。

圖3 包裹空心包體的巖芯Fig.3 A core enclosing a hollow inclusion

表1 主應力值

表2 實測三維應力分量

2.2 三維地質建模

通過對五老峰隧道地區(qū)的深入的地質調查，結合本次實驗測點位置和模型建立原則[6]。試驗點位于五老峰隧道出口段加寬道位置，樁號K27+982～932，埋深大約900 m，根據(jù)隧道圖紙，確定模型y軸與隧道軸線平行，x軸垂直于隧道軸線，埋深方向為z軸。隧道區(qū)域實際位置處于模型中間位置，遠離邊界，實際測點位于隧道中段，同樣處于模型中央位置。五老峰隧道布置圖及模型計算邊界范圍如圖4所示。

通過FLAC3D建立五老峰隧道部分區(qū)域的三維地質模型，如圖5所示。模型總計165 000個網(wǎng)格、174 216個節(jié)點、21 100個面。此次建立三維模型既具有計算精度上的優(yōu)勢，又避免了冗長的計算時間。

根據(jù)地勘報告和現(xiàn)場的取樣室內試驗，得知反演區(qū)域的花崗巖部分巖石力學性質。按照室內試驗結果和以往模擬經(jīng)驗，在參數(shù)選擇上，重度取 2 650 N/m3，黏聚力0.5 MPa，內摩擦角35°，本構模型采取彈塑性巖體模型(莫爾庫倫模型)。其中彈性模量和泊松比在本節(jié)不進行選取，這兩個參數(shù)需要利用BP網(wǎng)絡和FLAC進行推算，需要后期試算確定范圍。

病理學主要是研究疾病的病因、發(fā)病機制、病理變化及過程，是以形態(tài)學為基礎，通過病變大體標本和細微結構的觀察認識，掌握疾病的本質及發(fā)生規(guī)模對許多疾病的病理特點的觀察、掌握，要從宏觀和微觀的角度同時的觀察大體標本和顯微境下的特點，因此，病理學是基礎醫(yī)學和臨床醫(yī)學之間的紐帶[2]。

本文在模型邊界條件上主要選擇自重因素和構造運動這兩類主要因素施加，地形地貌在三維地表建模中得以實現(xiàn)，巖石力學性質則在模型賦參數(shù)這一環(huán)節(jié)得以實現(xiàn)。具體到模型上分為重力值、x向位移、y向位移，巖體的彈性模量和泊松比。

圖4 五老峰隧道區(qū)域三維示意圖Fig.4 Three dimensional schematic diagram of Wulaofeng tunnel area

圖5 五老峰隧道三維地質模型Fig.5 Three dimensional geological model of Wulaofeng

2.3 邊界荷載法人工試算

由于數(shù)據(jù)的稀少，本文研究只能采取人工試算的方法使得實測點的計算值逼近測量值，又因為直接邊界調整法本就繁雜，要想僅通過邊界調整法找到理論上最佳的反演邊界條件，使得測點的計算值很好擬合測量值會耗費大量的時間和精力，當反演區(qū)域地質條件復雜，測點數(shù)量較多時，這種僅靠人工試算的方法更加不現(xiàn)實。所以，需要在直接邊界調整法的基礎上，通過人工試算確定邊界條件的大致范圍，使得實測數(shù)據(jù)處在這個計算范圍之內，構建訓練數(shù)據(jù)，然后通過BP網(wǎng)絡反推最佳邊界。

經(jīng)過若干次邊界調整計算，確定了兩種巖體參數(shù)(巖體的彈性模量和泊松比)，3種邊界條件(重力值、x向位移、y向位移)的大致范圍。每個因素在波動范圍內取4個水平。具體值如表3所示。

表3 邊界條件水平

2.4 邊界荷載-BP神經(jīng)網(wǎng)絡聯(lián)合反演

2.4.1 網(wǎng)絡訓練樣本構建

2.3節(jié)確定了邊界條件和部分模型參數(shù)的范圍并且給每個因素劃分了4個水平，BP網(wǎng)絡訓練需要大量樣本，所以利用這5個因素4個水平(表3)人為構造訓練樣本，這樣一來，邊界荷載法和BP神經(jīng)網(wǎng)絡的優(yōu)勢互補，既避免了邊界荷載法大量試算的缺點，又通過前期試算給BP網(wǎng)絡提供的樣本數(shù)據(jù)。樣本構造則按照正交實驗法則[16]，具體如表4所示。

上述16組均勻設計組合充分地考慮了5個因素的4個水平，但是為了增加BP網(wǎng)絡訓練的可靠性和精度，另外的增加了4組隨機水平組合(和上述16組不重復)，如表5所示，從而使訓練樣本數(shù)量達到了20組，基本滿足BP網(wǎng)絡的訓練樣本要求。

表4 16種邊界條件組合

表5 4組隨機邊界組合

通過上述步驟，得到了20組有完整輸入輸出指標的BP神經(jīng)網(wǎng)絡訓練數(shù)據(jù)，接下來就可以構造關于三維應力和邊界條件的非線性網(wǎng)絡結構，利用這些樣本數(shù)據(jù)對網(wǎng)絡進行訓練。

2.4.2 網(wǎng)絡構造及訓練

網(wǎng)絡的構造和訓練則利用MATLAB軟件實現(xiàn)，網(wǎng)絡模型采用最高效的三層模型，具體參數(shù)如表7所示。

分析本文訓練好的網(wǎng)絡(圖6和圖7)，訓練、驗證和測試的擬合優(yōu)度R都在0.9以上(R越接近1說明擬合越好)，網(wǎng)絡訓練的誤差基本集中分布在0.03附近，極少數(shù)樣本誤差到了1或者2，但都在允許范圍內。這樣的結果說明本次訓練的BP網(wǎng)絡達到了期望效果，可以進行下一步的最佳邊界條件反演，具有工程指導意義。

表6 BP輸入、輸出樣本

表7 神經(jīng)網(wǎng)絡參數(shù)

2.4.3 網(wǎng)絡構造及訓練

2.4.2節(jié)中通過BP網(wǎng)絡得到了三維應力分量和邊界條件的非線性關系，雖然這種關系無法利用函數(shù)關系式表達，但是可以利用保存的網(wǎng)絡映射輸入與輸出之間的關系。本文最終目的是利用訓練完成的網(wǎng)絡，將測點的真實三維應力分量作為輸入帶入網(wǎng)絡，就能得到三維模型在有限差分計算軟件中的最佳邊界條件和部分模型巖體參數(shù)。

將表2代入訓練好的網(wǎng)絡，利用MATLAB命令a=sim(net，b)，其中a為輸出(邊界條件、巖石力學參數(shù))，b為輸入(三維應力分量)，net為訓練完成的網(wǎng)絡，可以得到最佳邊界條件。由于BP網(wǎng)絡的非確定性，決定進行了10次邊界條件反推，具體結果如表8所示。

圖6 BP網(wǎng)絡擬合回歸圖Fig.6 Regression diagram of BP network fitting

通過訓練好的BP網(wǎng)絡進行了10次反推，得到了上述10組邊界條件和巖石力學參數(shù)。到這一步，通過直接邊界調整法和BP網(wǎng)絡的結合成功反推出邊界條件，但是這些邊界條件是否真的滿足反演精度要求？與實測值相比又如何？還需要將這些邊界條件和巖石力學參數(shù)賦予三維地表模型，通過有限元軟件計算來檢驗這些條件的有效性，進一步驗證網(wǎng)絡的有效性。具體的計算結果如表10所示。

從表9中可以看出，10組邊界條件下三維應力計算結果與實測值總體上是相一致的。除去個別誤差稍大外，自重應力(σzz)平均模擬誤差1.24%，最大2.29%，最小0.09%；x向水平應力(σxx)平均模擬誤差1.04%，最大3.03%，最小0.04%；y向平應力(σyy)平均模擬誤差0.55%，最大1.22%，最小0.03%。其中第8組邊界條件下模擬結果與實測應力很是接近，三個方向上的應力誤差都控制在1%以下，y向水平應力甚至與實測應力相同，總體模擬結果精確在可接受范圍內。

圖7 BP網(wǎng)絡誤差分布圖Fig.7 Error distribution of BP network

表8 邊界條件BP網(wǎng)絡輸出

表9 計算應力與實際應力對比

2.5 五老峰隧道區(qū)域應力場分析

利用FLAC3D的切片(cutting tool)功能，將沿五老峰隧道軸線方向的三個主應力進行豎向切片分析，如圖8所示。三維模型區(qū)域的隧道總長約1 400 m，走向與坐標y軸平行，最大埋深約900 m。

(1)從圖8(a)可以看出，五老峰隧道區(qū)域的垂直應力總體呈梯度分布，在接近地表區(qū)域的應力等值線起伏較大，這說明工程區(qū)域的地形地貌對淺層區(qū)域的豎直應力分布是有一定影響的，這主要是地形的起伏造成同一水平下上覆巖層的自重力不盡相同，而豎直應力大部分又是由重力造成。沿隧道軸線的豎直應力分布在23～28 MPa范圍內，波動較大且略大于重力，原因是因為深部構造運動強烈，與自重力共同作用下形成了豎直應力，但總體上豎向應力是符合應力分布規(guī)律的。

圖8 五老峰反演區(qū)域云圖Fig.8 Inversion area cloud map of Wulaofeng

(2)從圖8(b)可以看出，五老峰隧道區(qū)域的最小主應力也是呈梯度分布，與豎直應力不同的是，地表的高低起伏對水平的應力的分布影響不大，只在地表約100 m的位置造成應力起伏。隨著深度的增加，應力等值線越來越趨向水平，這說明五老峰區(qū)域的最小主應力主要是受水平構造運動影響，重力只會在山體表層對水平應力造成影響。沿五老峰隧道軸線最小主應力分布在28～29 MPa，應力基本在一條等值線上，方向呈水平，平均應力在數(shù)值上大于豎直應力，因為隧道埋深較大，越到深部，水平應力與豎直應力的差距就越大。

3 結論

為了解決在某些工程區(qū)域地應力測量難，實測數(shù)據(jù)缺乏的問題，采用邊界荷載法和BP神經(jīng)網(wǎng)絡的聯(lián)合反演法，既能解決邊界荷載法人工試算時間過長，又能通過神經(jīng)網(wǎng)絡產(chǎn)生大量相關數(shù)據(jù)樣本，最后以五老峰隧道為工程實例進行驗證，得到如下結論。

(1)巖體的原位三維應力場并非是簡單相關影響因素的線性疊加，在測量數(shù)據(jù)較少和地質條件復雜情況下，多元線性回歸法不適用于反演研究，邊界荷載-BP神經(jīng)網(wǎng)絡聯(lián)合方法可以有效地克服并反演原巖地應力狀態(tài)。

(2)采用人工試算邊界條件范圍可以有效減少單一邊界荷載法的計算時間，BP神經(jīng)網(wǎng)絡利用正交實驗法則提供的足量數(shù)據(jù)可以精確地反算最佳的邊界條件，進而模擬三維應力狀態(tài)。

(3)聯(lián)合反演的三維應力模擬結果與實測相比，自重應力平均誤差1.24%，x向水平應力平均誤差1.04%，y向水平應力平均誤差0.55%，模擬結果較為精確。