基于隨機場理論的頂管隧道施工地表變形特性分析

張軒煜，施成華，孫曉賀，彭立敏，鄭可躍，王祖賢，肖國慶

(中南大學 土木工程學院，湖南 長沙，410075)

隨著城市地下交通快速發展，頂管法在地下交通建設中的應用逐漸增多。但頂管頂進對周圍土體存在擾動[1-2]，導致周圍地層應力場變化及土體變形，進而引起地表變形，嚴重影響周邊建筑物及地表交通的安全，因此，分析和預測頂管隧道施工引起的地表變形具有重要意義。目前在頂管施工引發地表變形的研究中，主要方法有經驗公式法[3]、理論解析法[4]和數值分析法[5]。PECK[6]根據工程經驗提出了沉降槽呈近似正態分布的PECK公式，在工程實際中得到了廣泛應用；蔣凱等[7]采用普式壓力拱理論與土柱理論分析了長距離曲線頂管隧道施工中的力學行為；喻軍等[8]通過建立頂管施工數值模型，考慮施工參數的影響，分析了地表變形影響規律。這些方法通常將周圍地層土體視為均質材料，但實際上巖土體并非嚴格均質，在同一地層中，不同位置的土體物理力學性質存在不同程度的偏差，2種土體的物理力學性質偏差呈一定的相關性，即巖土體具有空間變異性和自相關性的特點[9]。為體現巖土體的這一特性，可采用隨機場理論進行研究。隨機場理論能夠將地層參數隨機化，通過土體參數的均值、方差、變異系數和相關距離表征土體空間變異性和自相關性，分析其對頂管施工地表變形的影響，因此，有必要采用隨機場理論研究頂管施工引發的地表變形。

既有研究在使用隨機場理論分析工程問題時，通常建立二維模型進行模擬分析。趙煉恒等[10]通過隨機場理論建立二維模型，分析了土體的空間變異性對分層邊坡的可靠度的影響；易順等[11]針對基坑開挖引起地表變形的問題，采用隨機場理論建立平面數值模型，研究了土體參數的空間變異性對基坑變形規律的影響；程紅戰等[12]采用隨機場理論建立二維模型，研究了地層土體的空間變異性對盾構施工所引發的地表變形的影響及規律。然而，二維隨機場無法描述巖土體在3個方向上的空間變異性和自相關性[9]，同時，頂管頂進過程是三維動態過程，二維隨機場模型不能體現頂管施工頂進過程對地表縱向變形的影響，因此，有必要基于三維隨機場理論充分考慮地層參數在各個方向上的空間變異性和自相關性，建立模擬施工頂進過程的三維頂管模型，研究頂管施工對地表變形的影響。

本文基于隨機場理論，通過有限元分析軟件建立模擬施工頂進過程的三維頂管模型，在此基礎上，采用高斯型自相關函數、KL級數展開法的隨機模擬方法對模型中的地層參數進行隨機性分析，揭示頂管頂進施工引發的地表隨機變形特性，分析地層參數、變異系數以及相關距離對地表變形的影響，提出考慮地層空間變異性的簡化計算方法。

1 隨機場基本原理及三維隨機場的實現

1.1 隨機場概述

隨機場理論最早由VAN MARCKE提出[13]，目前已被廣泛用于表征地層土體的空間變異性和自相關性[14]，這種特性在地層土體中表現為物理力學參數具有波動范圍，服從概率分布，并與鄰近土體存在關聯。隨機場理論可以很好地描述這一特性，并將同一地層中同一物理力學參數的隨機變量定義為一個隨機場[15-16]。

假定存在n個數據點(x1,x2,…,xn)，并定義K為這些數據點的協方差矩陣，則矩陣K對稱且正定，由此可將K進行Cholesky分解。

式中：L和U分別表示上三角矩陣和下三角矩陣；角標T 表示矩陣的轉置。定義隨機向量y=Lv，其中，v為服從標準正態分布N(0,1)的隨機向量，則yyT(N階方陣)的期望值由下式得到：

由于v為服從N(0,1)分布的隨機向量，故

其中，I為單位矩陣。通過協方差矩陣K求得隨機向量y，再根據物理力學參數具體值即可求出點x1，x2，…，xn處的隨機模擬值。

為得到協方差矩陣K，需確定合適的自相關函數，以描述地層中任意兩單元之間的相關性[17]，同時，隨機場理論上可以離散為無限個單元點。為了減小計算量并控制精度，需選擇合適的離散方法將隨機場離散為n個數據點。

1.2 自相關函數

采用隨機場理論建立模型，需使用自相關函數確定隨機模型中任意2點之間的相關性。在既有二維隨機場的研究中[11-12]，計算土體相關性時，只考慮了2個方向的相關性，存在一定的局限性。三維隨機模擬方法考慮了土體3個方向(x、y、z)的相關性，能夠更好地模擬地層土體在空間上的變異性。

由于高斯型自相關函數具有良好的連續性和平順性，因此，選用高斯型自相關函數表征地層中任意兩土體之間的相關性大小，其表達式[18]如下：

式中：ρ(τx,τy,τz)表示地層中2 點之間的相關性；τx、τy、τz分別表示地層中2點在x、y和z3個方向上的距離；sx、sy、sz分別表示地層中2 點在x、y和z3個方向上的相關距離。隨機場中任意2個單元的相關性只取決于它們之間的距離。假設將1個隨機場離散為n個隨機場單元，每個單元中心坐標為 (xi，yi，zi)，則整個隨機場中所有單元之間的自相關系數將組成如式(6)所示協方差矩陣：

式中：τx12表示點1 與點2 在x方向上的距離，因此，τx12=τx21，ρ(τx12,τy12,τz12)=ρ(τx21,τy21,τz21)，協方差矩陣是元素為1的對稱矩陣。

1.3 基于KL級數展開法的隨機場離散

目前，隨機場離散的方法較多，其中，KL 展開法具有精度和效率均高的優點[18-20]，因此，采用KL展開對隨機模型進行離散。

1.3.1 KL級數展開公式

首先將研究對象離散為一組獨立標準正態隨機變量[19]，即

式中：xi為需要進行離散的地層參數；xi(x，y，z)為隨機場xi的均值；ΔHxi(x，y，z;θ)為隨機場xi的波動部分，其本質是具有原隨機場相關特性的一個零均值隨機場；δXi為隨機場xi的均方差；M為隨機場的級數展開項數；λj和fj(x)分別為H(x，y，z；θ)的特征值和特征函數。

1.3.2 特征值和特征向量的求解

為確定相關函數的特征值λj和特征函數fj(x)，采用Fredholm二類積分方程對展開的ΔHxi(x，y，z;θ)的協方差函數C(x1，y1，z1;；x2，y2，z2)求解，如式(8)所示[21-22]：

式中：Ω為積分區域。

1.3.3 級數展開項數的確定

隨機場的離散理論上可以進行無窮多項展開，但為了減小工作量并保證計算精度，通常需確定級數展開過程中的截斷項數n。對此，通過以下方法[23]求解隨機場的離散誤差ε：

式中：N為隨機場模型離散后的單元數目；ε(xj)為離散后單元xj的均方誤差。

1.4 三維隨機場的實現

既有研究在使用二維隨機場分析工程問題時，只考慮了2個方向上的空間變異性。為表征地層土體空間變異性在三維空間上對頂管施工地表變形的影響，考慮3個方向的空間變異性，建立三維隨機場模型分析工程問題。

圖1所示為隨機場理論的實現流程。通過地勘報告獲取各地層物理力學參數的均值、方差、變異系數和相關距離，根據均值建立模擬施工頂進過程的三維頂管模型，對確定參數的模型進行分析，得到確定性分析結果。同時，根據均值、方差、變異系數和相關距離，選擇合適的自相關函數，生成表征同一地層中各土體之間相關性的協方差矩陣，并通過地層土體的均值和變異系數，利用KL 級數展開法對隨機場進行n項展開，最終得到該地層土體的隨機參數。將生成的隨機參數賦予確定性分析模型中的單元，生成地層土體具有空間變異性和自相關性的隨機性分析模型，即可通過計算得到隨機模型不同頂進階段的橫向和縱向變形，進一步分析頂管頂進過程和土體隨機性對地表變形的影響規律。

圖1 隨機場理論實現流程Fig.1 Implementation processes of random field theory

2 基于三維隨機場的頂管案例分析

為驗證隨機模擬方法用于研究頂管施工地表變形的適用性以及相比于確定性分析的優越性，以南寧市地鐵1號線金湖廣場站Ⅲ號出入口頂管工程[24]為例，采用有限元軟件建立頂管隧道施工三維有限元模型。通過隨機場理論各生成500組[11]地層參數，計算隨機模型地表沉降并評估地表變形超標概率，驗證隨機場理論的可靠性。

2.1 工程概況

該工程采用土壓平衡多刀盤頂管機施工，頂進長度為37.2 m，管節覆土厚度約為5 m，管節外輪廓長×寬為6.9 m×4.9 m，厚度為0.5 m，管節長度為1.5 m，南寧地鐵1 號線金湖廣場站Ⅲ號出入口概況如圖2所示。

圖2 南寧地鐵1號線金湖廣場站Ⅲ號出入口概況Fig.2 Entrance III overview of Nanning metro line 1 in Jinhu Square Station

2.2 模型假設和建立

2.2.1 確定性分析模型的建立

為便于模型建立和分析，在建模時將地層簡化為均勻厚度的土體。假設土體是密實均質的彈塑性材料，采用Mohr-Coulomb本構模型，地層損失沿頂進方向均勻分布，在頂管頂進過程中不考慮頂管的緩慢推進，每次頂進距離為管節的長度。在頂進過程中，頂進面設置均布荷載模擬頂推力，且在頂進過程中保持不變。圖3(a)所示為該頂管工程的三維有限元數值模型，由7個地層組成，設置模型側面邊界為法向約束，模型底面為3個方向的固定約束，上表面為自由面不設置約束。管節采用彈性本構模型，各地層土及管節物理力學參數如表1所示。

表1 各地層土及管節物理力學參數Table 1 Physical and mechanical parameters of soil and pipe joints in each stratum

圖3 三維有限元數值模型Fig.3 3D finite element numerical model

2.2.2 隨機性分析模型的建立

為驗證隨機場理論分析頂管施工中地表變形的可靠性，根據上述隨機場理論，擬定對地層單元土體的天然密度、內摩擦角、黏聚力和彈性模量進行隨機化，同時，根據工程實際地層土體屬性，取水平相關距離為10 m，豎直相關距離為2 m[25-27]，地層土體參數中重度、黏聚力、內摩擦角和彈性模量的變異系數分別為0.02、0.10、0.10和0.30[28]，并生成500組隨機場進行分析。為體現隨機性分析中單元土體參數的變化規律特征，以泥巖地層中單元土體的彈性模量(E)為例說明各地層土體所存在的空間變異性，如圖3(b)所示(其余地層單元土體的各物理力學參數也表現出類似規律)。

2.3 三維隨機場分析

頂管施工段地表沉降現場監測點布置如圖4所示。通過數值模擬得到三維沉降模型，并對此進行橫向沉降和縱向沉降分析，其地表沉降示意圖及縱向、橫向分析截面如圖5所示。

圖4 頂管施工段監測點布置圖Fig.4 Monitoring point layout of pipe jacking construction section

圖5 確定性分析地表沉降及分析截面Fig.5 Deterministic analysis surface of settlement and cross sections for analysis

2.3.1 縱向分析

圖6所示為頂管頂進到不同位置時的縱向地表沉降。從圖6可以看出：確定性分析曲線與隨機性分析曲線的變化規律一致，隨機沉降曲線在確定性分析沉降曲線附近波動，確定性分析曲線基本位于隨機曲線中間位置。通過比較不同頂進距離下的地表沉降曲線可知：各橫斷面的地表沉降范圍隨著頂管的頂進逐漸增大，并在頂管頂進完成后沉降不變。以距始發井0 m的橫斷面為例，當頂管頂進5 m 時，其地表沉降的隨機范圍為-19.8～-6.7 mm；隨著頂管的逐步頂進，其確定性分析的沉降有所增大，隨機性分析結果的沉降和包絡范圍相應增大；當頂管頂進完成時，該斷面的隨機沉降范圍達-29.3～-19.8 mm；隨著頂管的逐步頂進，各斷面地表沉降和隨機沉降曲線的主要包絡范圍逐漸增大。由此可見，地層參數的空間變異性對頂管施工引發的地表變形產生較大影響。

圖6 隨機性縱向地表沉降Fig.6 Stochastic longitudinal surface settlement analysis

在不同的頂進距離下，地表變形形態均不同；當頂管開始頂進時，頂進面前方一定范圍內的土體會發生不同程度的地表沉降，且距離頂進面越遠，地表沉降越小。隨著頂管的頂進，地表變形向頂進方向延伸，未頂進段開始產生地表沉降，已頂進段的地表沉降變化較小，并在頂管頂進完成后地表沉降不變。如頂管頂進5.0 m時，距離始發井0～21.0 m 的縱向范圍內均出現地表沉降，最大沉降為-19.8 mm；頂進15.0 m時，縱向0～27.0 m的范圍內出現了地表沉降，最大沉降為-24.7 mm。

2.3.2 橫向分析

圖7所示為頂進完成后各監測斷面的監測點沉降和確定性分析、隨機性分析結果，表2所示為頂管頂進完成后各監測斷面確定性分析、隨機性分析的各監測斷面沉降模擬值和監測點監測值。從圖7和表2可以看出：確定性分析和隨機性分析所得出的地表沉降趨勢與監測的沉降趨勢較符合，但由于模型建立時受地層的簡化等因素的影響，確定性分析所得沉降與斷面監測值在最大變形處存在差異。隨機性分析在確定性分析的基礎上，給出了地表變形的范圍，能夠將監測值較好地包絡，可對地表變形進行較好預測。如對于監測斷面A，其監測點沉降變化趨勢與確定性分析和隨機性分析的沉降曲線變化趨勢一致，確定性分析和隨機性分析所得到沉降的最大沉降基本位于頂管中軸線附近，監測值為23.8 mm，確定性分析最大沉降為-26.4 mm，隨機沉降曲線的包絡范圍為-33.1～-17.3 mm，確定性分析所得到結果與監測結果存在一定的差距；隨機性分析所得到的包絡范圍能夠將監測值包絡，較好地預測頂管施工中的地表變形。

表2 頂管頂進完成后各監測斷面沉降Table 2 Settlements of each monitoring section after completion of pipe jacking mm

圖7 橫斷面沉降隨機性分析結果Fig.7 Cross-sectional settlement randomness analysis results

2.4 三維隨機場超出變形控制值概率分析

相比于二維隨機場，三維隨機場考慮了3個維度的空間變異性和頂管施工的頂進過程，對此分別從縱向和橫向分析地層參數空間變異性對地表變形的影響。為進一步量化施工過程中地表變形超出控制值的概率，通過式(11)[12]計算超出控制值的概率P：

式中：N為每個工況中隨機計算的次數(N=500)；Nf為計算結果大于最大沉降控制值的次數。

對頂管頂進過程中各監測斷面地表變形的超標概率進行計算，計算結果如圖8 所示。從圖8可以看出：根據頂管頂進過程中各斷面地表沉降控制值，可得到其可能超出控制值的概率；以頂管頂進到25 m 為例，當監測斷面A的控制值為-20 mm 時，其失效概率高達90%，當控制值為-25 mm 時，失效概率為13%；對于橫斷面C，沉降控制值7.5 mm 所對應的失效概率為91%，控制值-11 mm 所對應的失效概率為18%。由此可在施工頂進過程中，預測各監測斷面地表沉降超出控制值的概率。

圖8 地表沉降超標概率分析結果Fig.8 Analysis results of probability of surface settlement exceedance

3 頂管隧道地層參數空間變異性的影響分析

為分析不同地層參數隨機化對頂管施工中地表變形的影響，在上述案例確定性分析模型的基礎上，分別對地層土體參數(彈性模量、內摩擦角、黏聚力)、相關距離和變異系數單獨隨機化，對比分析各因素對地表變形的影響。隨機性分析參數如表3所示[11,29-30]。為便于比較各地層參數隨機化對地表變形的影響，取頂進完成后B測點橫斷面地表沉降進行比較分析。

表3 隨機性分析參數[11, 29-30]Table 3 Parameters of randomness analysis

3.1 地層土體參數不確定性的影響

圖9 和圖10 所示分別為地層土體內摩擦角、黏聚力和彈性模量隨機化下頂管施工地表沉降曲線和地表水平位移曲線。從圖9和圖10可以看出：確定性分析曲線基本位于隨機沉降范圍的中間位置，隨機沉降曲線在確定性沉降曲線附近波動，其沉降變化規律與確定性曲線變化規律相符；在相同變異系數下，彈性模量的隨機化對地表變形有較大影響，能夠使隨機曲線更加離散，隨機曲線所包絡的范圍更廣，其隨機沉降包絡范圍達-25.9～-16.1 mm；黏聚力和內摩擦角對地表變形的影響相對較小，其隨機沉降包絡范圍分別為-25.1～-19.6 mm 和-23.8～-17.6 mm；地層土體物理力學參數的隨機化不會改變地表變形的變化規律和變化趨勢。由于工程實際中土體重度的變異性較低，其變異系數相比于彈性模量的變異系數要小很多，地表橫向離散曲線與沉降曲線基本一致。由此可知，土體物理力學參數中彈性模量的不確定性對頂管施工地表沉降的影響最大，其次是內摩擦角和黏聚力，且這兩者影響程度相近。

圖9 土體物理力學參數不確定性分析所得橫向地表沉降曲線Fig.9 Soil transverse surface settlement curves based on uncertainty analysis of physical and mechanical parameters

圖10 土體物理力學參數不確定性分析所得地表水平位移曲線Fig.10 Surface horizontal displacement curves based on uncertainty analysis of physical and mechanical parameters of soil

3.2 相關距離的影響

圖11 所示為地層參數在不同相關距離下的地表沉降曲線，圖12 所示為不同相關距離下的地表水平位移曲線。從圖11和圖12可以看出：相關距離對地表變形的離散程度存在一定的影響，當Sx=Sy=10，Sz=2 時，隨機沉降曲線的隨機范圍最大，為-28.9～-14.1 mm；隨著水平或豎直相關距離減小，隨機沉降曲線的沉降范圍有所減小。這是由于相關距離減小，地層中土體的相關性減弱，物理力學參數的離散性增大，導致隨機計算結果更加離散。不同相關距離也會在一定程度上改變地表變形的形態，這是由于相關距離的改變會使得地層土體出現整體均質、局部異質的特點，進而導致地表變形曲線不對稱。因此，相關距離對地表變形的影響不可忽視，相關距離越大，對地表變形的影響越大。

圖11 不同相關距離下橫向地表沉降曲線Fig.11 Transverse surface settlement curves at different correlation distances

圖12 不同相關距離下地表水平位移曲線Fig.12 Surface horizontal displacement curves at different correlation distances

3.3 變異系數的影響

地層參數的空間變異性通過變異系數予以表征。圖13 所示為不同變異系數對地表沉降曲線的影響，圖14 所示為不同變異系數下地表水平位移曲線。從圖13和圖14可以看出：變異系數不會改變既有的地表變形規律和趨勢，其對地表變形的影響主要表現為隨機曲線離散程度不同；當變異系數為0.1時，隨機性分析的變化趨勢與確定性分析一致，沉降曲線相對集中，最大沉降為-24.4～-18.7 mm，變化范圍較??；隨著變異系數增大，地表沉降趨勢與確定性分析結果無明顯差別，但曲線分布更加離散；當變異系數為0.4時，地表最大沉降為-27.1～-14.0 mm。這是由于變異系數增大，地層參數的不均勻性增加，使得計算得到的地表變形結果更加離散，因此，變異系數越大，地表變形結果更加離散。不同變異參數下的地表最大沉降統計如圖15所示。

圖13 不同變異系數下橫向地表沉降曲線Fig.13 Transverse surface settlement curves under different coefficients of variation

圖14 不同變異系數下地表水平位移曲線Fig.14 Surface horizontal displacement curves under different coefficients of variation

圖15 不同變異參數下地表最大沉降統計Fig.15 Statistics of the maximum surface sedimentation with different coefficients of variation

4 結論

1) 基于高斯型自相關函數、KL級數展開法的三維隨機模擬方法能夠有效表征地層參數變異性和自相關性，得到頂管施工所引起的地表縱、橫向變形特性，反映地表變形規律和變形范圍的影響。

2) 相比于確定性分析方法，隨機性分析所得結果能夠充分表征頂管施工所引發地表變形的趨勢、數值和波動范圍，更加貼合工程實際中較復雜的情況，所得結果較確定性分析結果更具有借鑒作用和參考價值。

3) 與二維隨機分析結果相比，三維隨機分析考慮了地層中3個維度的空間變異性以及頂管施工的動態影響，可從縱向和橫向2個角度分析不同頂進距離下各縱斷面和橫斷面的地表變形，并能得出地表變形超標概率。

4) 彈性模量的空間變異性對地表變形隨機曲線的離散程度和變形范圍有較大影響，但不改變變形的變化趨勢和規律；變異系數和相關距離增大能夠使得地表變形隨機曲線更加離散。