王呈璋
(山東理工大學交通與車輛工程學院,山東淄博 255000)
隨著調水工程輸水隧洞、深埋長距離交通隧道、深部礦產能源等大型地下工程的不斷開發建設[1-5],越來越多的工程面臨較為復雜的地質構造環境,而工程區域初始地應力場的分布特征受區域地質構造的直接影響。如何能夠準確地獲得工程所需的初始地應力場量值與分布規律,一直都是地下工程研究的基本問題之一。在地下工程的設計與施工中,地應力都被視為最基礎的荷載因素,工程設計是否有效可靠,施工過程是否合理安全,均受其直接的影響[6-8]。
可靠的地應力實測數據可以為工程建設提供有效的參考[9],但是由于工程區域復雜環境的限制,經費的要求以及試驗技術的局限性,無法開展大量的地應力實測工作,而有限離散的地應力實測數據往往無法滿足實際工程的需要[10-11]。因此,在地應力實測資料的基礎上,結合地質構造背景、地形地貌特征、巖性分布規律等地質資料,以數學力學方法為數值計算的手段,反演出有效的工程區域初始地應力場,是初始地應力場反演的基本思路[12]。在眾多的初始地應力場反演方法中,比較具有代表性的有應力函數法、側壓系數法、邊界荷載法、邊界位移法、神經網絡法、遺傳算法、變差函數法、徑向基函數法等[13-17],其中應力函數法是較為經典的初始地應力場反演方法。應力函數法在工程實踐過程中,被不斷的改進發展,最初的應力函數法較為單一,僅僅為根據彈性理論設定的多項式,通過對多項式的系數求解,進而求解出應力函數的具體形式。隨著地應力反演計算相關理論的發展,以及計算條件和計算手段的進步,有關學者通過構建更為復雜多樣的應力函數,如三維正交多項式函數,應力趨勢函數等,利用計算機高效的計算效率進行初始地應力場的反演。近些年來,伴隨著大型通用計算軟件的發展,應力函數的形式也變得多樣化,應力函數從最初的反演目標,逐漸轉變成為反演手段,將應力函數與有限元聯合進行初始地應力場反演,成為目前較為主流的反演方法之一。
基于相關學者在初始地應力場反演過程中對于應力函數法的研究成果,本文從基本思想、實現過程、工程應用特點和優缺點等方面,著重介紹了該方法在不同階段的重要特征。
初始地應力場的應力函數法回歸反演分析,即根據有限且離散的地應力實測數據,結合工程區域的地質計算模型,以數學力學方法為主要手段,構建三維應力場函數或應力分量函數的表達式,利用數值計算方法,推導出工程區域的初始地應力場[18]。在初始地應力場反演分析的方法中,應力函數法的基本概念提出的時間較早,經過大量學者的豐富研究,現在已經具有多樣的形式,如多元線性回歸分析法、應力函數趨勢分析法、應力函數數值分析法、邊界荷載分析法等,以上的方法均采用不同的應力函數形式作為數值計算方法的核心。
數值計算方法多種多樣,其中較為典型的應力函數法最早是由郭懷志等[19]提出,他認為巖體在彈性工作狀態下,初始應力場可視為因變量,若干自變量共同作用形成應力場這一復雜的因變量,這些自變量分別為彈性模量、泊松比、容重、地溫等物理因素,以及地質構造特征、地形地貌特征、構造作用下巖體的三維坐標等構造因素,各種因素在地質體內產生的應力視為基本初始應力,而初始應力場即為各基本初始應力的線性疊加,其表達公式為
(1)

基于這種回歸反演分析理論,初始地應力場的反演問題即轉變為對回歸系數bi的求解。具體實現過程為,首先將各自變量因素作用于地質計算模型,得到對應狀態下的基本初始應力計算值,根據線性疊加原理求得初始地應力場的計算值,結合地應力實測數據,將二者的差值利用最小二乘法原理使殘差的平方和最小,即可得到回歸系數bi的最優解。
張有天[20]等提出了一種根據趨勢對地應力場進行分析的應力函數法,即地應力場的趨勢分析法。該方法的核心思想是假設高階的多項式應力函數,結合有限的地應力實測數據以及邊界約束條件,根據趨勢對地應力場進行反演分析。該方法首要的任務是假定合理的四次多項式函數φ1、φ2、φ3,三個應力函數應該滿足Beltrami-Michell方程,應力分量由應力函數表示的公式如下:
σx=?2φ3/?y2+?2φ2/?z2
σy=?2φ1/?z2+?2φ3/?x2
σz=?2φ2/?x2+?2φ1/?y2
τyz=-?2φ1/?z?y
τzx=-?2φ2/?z?x
τxy=-?2φ3/?x?y
(2)
由式(2),六個應力分量σx、σy、σz、τxy、τyz、τzx均可表示為關于空間坐標(x、y、z)的二次多項式。對于地應力實測數據,每一個測點都包含六個應力分量,當實測點為n個時,可以建立6n個獨立的方程。作為邊界約束條件的地表點,由于其法向應力與剪切應力均為零,因此可以選取若干地表點建立方程組。當六個應力分量的多項式系數個數與所有方程的個數相同時,即可求得六個應力分量的具體表達式。
肖明等[21]在對錦屏二級水電站壩址處山體進行初始地應力場反演回歸分析時,提出了三維應力函數擬合方法。該方法主要分為兩個部分,首先根據工程區域范圍內地勘資料所呈現的地質體地表剝蝕風化產生的地貌特征、時空變化下的構造演化特征和巖土物理力學性質,采用有限元反演出初始地應力場。然后根據測點位置的反演計算值和實測值,運用三維有限元對計算后的地質體模型進行插值運算,利用三維正交多項式構建出能夠合理反映三維應力場的應力函數。假定一組多項式序列p0(x,y,z),p1(x,y,z),...,pm(x,y,z)滿足以下的正交條件:
(3)
式中:n為擬合區域離散點的個數;當i=j時,p為常數;當i≠j時,p為零。
由正交多項式的性質構造的正交多項式系列如式(4)所示:
p0=1
p1=(z-a1)p0
p2=(z-a2)p1+(x-b2)p0
?
pi+1=(z-ai+1)pi+(x-bi+1)pi-1+(y-ci+1)pi-2+di+1pi-3
(4)
式中:ai+1、bi+1、ci+1、di+1為待定系數。
由數學歸納法求得空間坐標(x、y、z)表示的多項式系數和多項式pi的具體表達形式,從而擬合出三維應力函數的具體表達式:
Qi=αi0p0+αi1p1+…+αimpm
(5)
式中:αi0,αi1,…,αim為待定的應力系數。
通過建立以地應力實測數據與擬合邊界的應力值為極值條件的目標函數,利用最小二乘法的原理對目標函數求極值,即可得到式(5)中的應力系數,進而得到三維應力場的應力函數表達式。
邊界荷載反分析法的核心思想是地質體模型的邊界荷載與初始地應力狀態是相對應的關系,如果已知地質體模型的真實邊界荷載,則可得到初始地應力場。因此,邊界荷載法的實現過程就是通過不斷地調整地質體模型的邊界荷載,使得地應力測點位置處的應力回歸計算值與實測值之間的誤差達到最小。
劉允芳[22]等對于該方法的具體實現過程為,根據地質勘測資料,選擇合理的工程區域范圍,構建地質力學計算模型。對時空變化下的構造演化過程進行分析,確定可能對工程區域地應力場的形成產生影響的主要地質構造因素,一般選擇自重應力、兩個水平方向的構造應力和水平面內的剪切應力,將以上四種因素作為待定因素。對于每一種待定因素,將其作用于地質計算模型,求得在已知測點位置處的地應力值。將測點位置處所有計算求得的地應力值與現場實測的地應力值之間建立回歸方程:
(6)

利用最小二乘法原理使得地應力測點的實測值與回歸計算值之間的殘差達到最小值,得到對應于四個待定因素的回歸系數,即可得到應力函數的回歸方程。
由于地應力場具有時間性和空間性二重屬性,所以就整體而言,地應力場是非線性且非穩定的應力場[23]。在復雜漫長的地質歷史年代,地應力場受多重因素的作用,如地殼的運動、地溫的不均勻變化、地幔的隨機對流、地表的風化剝蝕和地形地貌的變化、巖體自重的影響、巖體內孔隙水壓力以及節理結構面的影響等[24-25],如果在地應力場反演工作中全面考慮這些因素的影響,理論上是可以以較高的精度反演出與實際情況相符合的地應力場,但是在實際操作中可行性較低。一方面,地應力場的非線性特征很難用單一的非線性函數或者非線性模型進行模擬表征。另一方面,地應力場又是整體上的非穩定場,無法準確的表示出地質歷史年代整個應力場的變化過程。但是,在研究具體的工程問題時,需要考慮的地質變化過程一般較短,而在所考慮的時間尺度和空間尺度的范圍內,地質體是相對穩定的狀態。所以,在對具體的工程問題進行初始地應力場反演時,完全可以對各種因素進行簡化處理,將工程地質體視為相對穩定的非線性模型。
應力函數法就是在這樣的背景下提出來的,不同形式的應力函數均可以在一定程度上模擬地質體的非線性特征。應力函數趨勢分析法是最為直接的方法,該方法用四次的函數來模擬三維的應力場,基于彈性力學的假定,根據平衡方程、物理方程、相容方程來確定應力分量的具體表達式,對于地質體模型中的任意一點,均可由其坐標值得到對應的應力值,從而得到整個過程區域的初始地應力場。該方法所需的計算工作量較小,簡便易行,具有較強的可操作性,表示的是連續空間的地應力場。缺點在于將地質體視為簡單的彈性模型,并未考慮巖體的塑性變形。而且地質體內部存在的斷層、褶皺等地質構造特征均無法充分的考慮,所以該方法的代表性不強,反演精度較低,適合在巖性比較均勻,斷層發育較少、褶皺等構造地質構造作用不夠強烈的地殼深部工程使用,如王拉才[26]等在撫順龍鳳煤礦、孔廣亞[27]等在某黃金礦的的地應力場反演中均采用該方法,取得了良好的效果。
應力函數數值分析法從地應力實測點出發,雖然離散且有限的應力實測數據無法滿足三維有限元計算的要求,但是完全可以利用正交多項式的特性,根據實測數據構建正交多項式序列,并且結合插值運算擬合出一個應力函數。該方法在錦屏二級水電站地應力場反演中首次運用,并且分別考慮了存在斷層影響和不存在斷層影響兩種情況。反演結果表明采用此方法反演回歸的初始地應力場在兩種情況下的結果基本一致,僅僅在靠近斷層附近的主應力值有所差別,說明改方法具有一定的精度,并且便于操作。但是,該方法的局限性在于建立三維地質體模型的過程中對地質特征進行了較多的簡化處理,對于地形地貌特征以及構造演化特征較為復雜的地質體,一個擬合的應力函數無法準確描述整個初始應力場的真實狀態。趙辰[28]等以該方法為基礎,保留其核心的思想,即構建正交多項式系列,通過正交多項式序列擬合出函數的具體形式,從而提出了考慮地層剝蝕作用的側壓力系數法這一地應力場反演方法。該方法擬合的函數不再是應力函數,而是側壓力系數函數。對于復雜的地質構造環境而言,初始地應力場的變化可能比較復雜,但是在一個相對穩定的區域內,側壓力系數這一代表水平應力與自重應力比值的量值保持在穩定的范圍內,因此完全可以通過建立側壓力系數擬合函數來進行地應力場的反演。利用有限元軟件通過逐次開挖來模擬地表風化剝蝕的過程對初始地應力場的影響,以應力實測點的側壓力系數數據為基礎,采用正交多項式序列構造出三維空間的側壓力系數擬合函數,由水平應力、自重應力和側壓力系數之間的數學關系,即可得到整個地質模型的初始地應力場。該方法的顯著優點就是可以充分考慮地形地貌特征、地質構造條件對初始地應力場的影響,反演方法具有一般性,適用于各種復雜地質體,反演結果具有較高的精度和合理性。
多元線性回歸分析初始地應力場的方法,由于考慮多種影響地應力場的因素,并且每一種因素都對應著唯一的應力狀態,將每一種應力狀態視為一個應力函數,因此通過線性疊加得到的對應于三維空間初始地應力場的應力函數,具有唯一性[29]。邊界荷載反分析法不僅具備多元線性回歸分析法的優點,而且其獨特的優點就是僅考慮對初始地應力場產生顯著影響的四個因素,包括自重應力因素、水平方向兩個構造應力因素以及水平面內的一個剪切構造應力因素,反演過程只需求得各因素作用于有限元地質模型的最優系數,即可獲得各因素共同作用下的唯一初始地應力場。該方法最早應用于水利水電行業,由于其采用有限元計算軟件建立逼真的工程地質模型,能夠反映地形地貌特征、斷層和褶皺等地質構造特征,并且采用科學可靠的數理統計理論做計算支撐,因此從提出來之后便被行業廣泛認可并普遍使用,如抽水蓄能電站[30]、引水隧洞工程[31]、水電站高邊坡[32]等大型水利水電工程的初始地應力場反演均采用該方法。雖然在其他工程領域,利用多元線性回歸分析法和邊界荷載反分析法對工程區域初始地應力場反演的研究工程相對較少,但是由于一些交通隧道工程正在面臨超長深埋等較為突出的地質問題,利用上述方法進行地應力場反演的優勢也正在凸顯。汪波[33]、代聰[34]和尤哲敏[35]等分別在在蒼嶺特長公路隧道、藍家巖深埋特長公路隧道和大坪山深埋公路隧道的初始地應力場反演分析中運用多元回歸分析原理,獲得了較為合理可靠的隧道軸線處初始地應力場分布規律。因此,多元線性回歸分析法,以及以此為基礎發展而來的邊界荷載法,由于能夠更為真實的反應工程區域的地質構造特征,更為合理可靠的反演出初始地應力場的分布規律,無論是規模較大的水利水電工程,還是深埋超長的交通隧道工程,都將其視為較為科學的方法并成功應用于工程實踐。
本文以應力函數法為基礎,詳細介紹了該方法的提出背景,系統地總結了該方法在具體的工程實踐中發展而來的多元線性回歸分析法、應力函數趨勢分析法、應力函數數值分析法和邊界荷載反分析法,具體闡述了每一種方法的基本原理和數值計算方法,對比分析了各方法的工程應用特點以及優缺點,具體對比結果如表1所示。主要結論如下:

表1 各方法對比分析
(1)四種方法的力學假定均是彈性力學理論,即假定工程區域的巖體僅發生彈性變形,不考慮巖體的塑性狀態。
(2)應力函數趨勢分析法是反演初始地應力場的直接法,直接假定應力函數進行求解。而另外三種方法均屬于間接法,多元線性回歸分析法和邊界荷載法通過線性疊加各種因素作用下的應力場建立初始地應力場的應力方程,應力函數數值分析法利用正交多項式的性質構建多項式應力函數。
(3)多元線性回歸分析法和邊界荷載法的基本原理一致,邊界荷載法相比前者而言僅考慮四種因素,即自重因素、兩個水平構造因素和一個剪切構造因素。邊界荷載法的可操作性較強,更加科學合理。
(4)邊界荷載法相對于其他三種方法,具有較大的優勢,也是目前最流行的方法。因為其考慮地表風化剝蝕作用、斷層和褶皺等地質構造作用,更能反映工程區域的真實地質狀態,適用于大型高陡邊坡工程、深埋超長隧道工程和水利水電工程區域的初始地應力場反演工作。
基于應力函數法的發展過程和具體工程實踐過程,仍有部分問題有待進一步的發展完善,使其能夠更加滿足工程實際的需要。
(1)目前許多大型工程逐漸面向深地,越來越大的埋深伴隨著較高的地溫,而地溫對于巖體初始地應力場的影響尚不明確。
(2)實際的工程地質模型,邊界荷載不是簡單的線性荷載,如何采用非線性荷載模擬邊界條件值得進一步研究。
(3)由于地殼淺部中的巖體多處于彈塑性狀態,忽略巖體的塑性狀態,僅僅以彈性力學假定進行初始地應力場的反演,反演結果與真實狀態相比仍有較大的誤差。因此,結合巖體的彈塑性狀態進行考慮,才能使反演結果更加合理。