近年来，随着我国经济建设的快速发展，国内城市化建设快速提高，地下隧道的建设在城市的发展中占比越来越大。由于地形限制，地下隧道施工过程中不可避免会受到既有建筑物、既有管道等影响，盾构掘进路线往往无法按照理想状态进行施工，很多情况下需要采用曲线隧道的方式来完成隧道掘进任务。相比于直线盾构隧道，曲线盾构隧道开挖需要进行纠偏超挖，更容易出现地表土体变形，过大的变形会严重影响地面建筑、交通的安全性，因此合理预测曲线盾构施工引起的地表土体变形尤为重要。
目前国内外对于盾构开挖引起的地表附加沉降的理论研究成果不少。魏纲等［1-2］基于Mindlin解，考虑开挖面附加推力、盾壳摩阻力以及地层损失，推导出盾构施工引起的地表变形计算公式。姜忻良等［3-4］在Sagaseta提出的镜像法的基础上采用数值积分的方法由地层损失推导出隧道周围土体位移解。冯浩等［5］在Peck公式的基础上，通过现场实验、拟合分析、二维数值计算等方法得出了更适合预测水平曲线盾构施工引起地表沉降的计算公式。路林海等［6］结合现场实测数据对Peck公式进行修正，得到了水平曲线盾构地表沉降预测公式，在黏性土地层有较好的适应性。潘泓等［7］通过对比分析直线隧道与水平曲线隧道临近土体分层沉降、水平位移、孔隙水压力的实测数据，揭示了小半径水平曲线盾构施工对邻近土体的扰动规律。赵伯明等［8-9］基于修正源汇法及Mindlin解推导出考虑盾构掘进角度的地表变形公式。胡达等［10］在Loganathan公式与Terzaghi理论的基础上，结合Maxwell-Betti 功互换定理，建立了考虑“土拱效应”的盾构隧道施工地表沉降预测模型。邓皇适等［11-12］建立了盾构施工水平曲线段地层损失模型，基于镜像法与Mindlin解推导小转弯半径盾构施工引起地表沉降计算公式。孙捷城等［13］基于Mindlin解， 考虑到盾构水平曲线掘进时开挖面不均匀附加推力与盾壳不均匀摩擦力，通过引入差异系数，推导出水平曲线盾构引起地表沉降。杨果林等［14］构建了考虑主动铰接机构的曲线隧道开挖力学分析模型，基于Mindlin解、镜像法以及孔口扩张理论推导出曲线线性条件下盾构掘进时的地表位移理论解。
上述研究中大多是分析直线隧道与水平曲线隧道盾构施工对地表土体变形的影响，目前关于竖向曲线盾构引起地表土体变形的研究还比较少。本文以北京东六环路改造工程为背景，考虑竖向曲线盾构转弯方向与转弯角度，对开挖面附加推力、盾壳与周围土体之间的摩擦力、盾尾注浆压力以及地层损失等影响因素进行分析，基于Mindlin解推导出竖向曲线盾构施工引起的地表变形计算公式，结合工程实例研究竖向曲线盾构施工引起的土体变形规律。
1 理论分析
1.1 竖向曲线盾构掘进计算模型
基于 Kelvin 解，Mindlin推导出在竖向集中力和水平集中力作用下弹性半空间内坐标为[(x’,y’,z’)]的任意点处的竖向位移[wv]和[wh]计算公式，其计算模型如图1所示。以[(0,0,c)]为集中力作用点处坐标，竖向集中力[Pv]引起土体内任一点处竖向位移[wv]则为：
[wv=Pv16πG(1-μ)3-4μR1+81-μ2-3-4μR2+]
[z’-c2R13+3-4μz’+c2-2cz’R23+6cz’z’+c2R25] （1）
水平向集中力[Ph]引起任一点土体竖向位移[wh]为：
[wh=Phx’16πG1-μz-cR13+3-4μz’-cR23-]
[6cz’z’+cR25+41-μ1-2μR2R2+z’+c] （2）
其中，[R1=x’2+y’2+z’-c2] ，[R2=x’2+y’2+z’+c2]，G为地层剪切模量（MPa），[μ]为泊松比。
以盾构机向上转弯掘进为例，竖向曲线盾构推进计算模型如图2所示，图中[x]轴正方向为盾构掘进方向，转弯角度为[α]，盾构直径为[D]，盾构机盾壳长度为[L]，开挖面附加推力为[q]，盾尾同步附加注浆压力为[p]，盾壳摩阻力为[f]。
计算模型基本假定如下：
（1）假定土体为各向同性线弹性半无限体，且不考虑土体固结排水影响；
（2）盾构机在竖向曲线开挖过程中，开挖面附加推力呈现不均匀性。假定开挖面附加推力以刀盘中线为分界线，上侧、下侧均满足均匀分布。其中上侧推力[q1]和下侧推力[q2]满足[q2=mq1]（[m]为上下两侧附加推力差异系数，[m]≥1）。
（3）盾构机在曲线掘进过程中，盾构机盾壳内侧受到的挤压作用比外侧更加强烈。假定盾壳摩阻力以隧道中线为分界线，上侧、下侧均满足均匀分布。其中上侧盾壳摩阻力[f1]和下侧盾壳摩阻力[f2]满足[f2=nf1]（[n]为上下两侧摩阻力差异系数，[n]≤1）。
（4）注浆压力均匀分布在管片外围，且仅影响一环管片的宽度。
<G:\武汉工程大学\2024\第5期\王慧-1.tif>[地面][X′][R2][R1][R1][（x′，y′，z′）][O′][c][c][Y′][Z′]
图1 Mindlin解计算模型
Fig. 1 Mindlin solution calculation model
<G:\武汉工程大学\2024\第5期\王慧-2.tif>[x′ x″][x′ x″][正面附加推力q][X][Z][O][α][H][y][y″][y′ ][z″][z′ ][上侧摩阻力f1][上侧推力q1][上侧推力q2][刀盘中线][下侧摩阻力f2][切口][盾尾][L][f][f][盾尾同步附加
注浆压力p]
图2 竖向曲线盾构掘进计算模型
Fig. 2 Calculation model of shield tunneling with
vertical curve
1.2 向上掘进引起地表沉降
（1）开挖面不均匀附加推力[q]
向上掘进时，开挖面附加推力引起土体变形的计算示意图如图3所示。图中[dA]为盾构机开挖面处任意微元体，该微元体面积大小为[dA=rdrdθ]。[r]为微元体距开挖面中心点的距离，[θ]为微元体与开挖面中心水平线之间的夹角，[H]为盾构机转弯前开挖面中心点处埋深。开挖面推力分解为竖向分力[qsinα]以及水平分力[qcosα]，则该微元体所受的集中力可以表示为[dPv+dPh=qsinαrdrdθ+qcosαrdrdθ]，本文利用计算软件Mathematica进行公式推导以及积分计算。
<G:\武汉工程大学\2024\第5期\王慧-3.tif>[x′ x″][X][Z][z″][z′ ][O][α][H][y][y″][y′ ][地表][c][z″][z][y″][o′ ][y ][o ][y′ ][dA][dr][r][θ][dθ][q1][q2][z′ ]
图3 向上掘进开挖面附加推力计算示意图
Fig. 3 Schematic diagram of additional thrust calculation on the upward excavation face
求解该微元体所受到的集中力在坐标系[xyz]某一点[x,y,z]处引起的地表竖向位移，需要基于Mindlin解对坐标系进行转换，在盾构机转过[α]角度后，以[xyz]坐标系中o点为圆心，建立坐标系[x’’y’’z’’]，平面 [y’’oz’’]与开挖面处于同一平面。对[x’’y’’z’’]坐标系进行变换为：
[x’’=x-Htanαy’’=yz’’=z] （3）
在开挖面处建立第二个局部坐标系[x’y’z’]，同理进行坐标变换为：
[x’=x’’y’=y’’-rcosθz’=z’’c=Hcosα-rsinθ] （4）
将式（4）代入式（1）中，由于曲线盾构开挖面上侧与下侧附加推力不相等，将开挖面划分为上下两侧同时积分，可推得盾构机向上掘进附加推力竖向分量引起的地表土体变形计算公式[wq1]。
将式（4）代入式（2）中，开挖面附加推力分为上下两侧同时积分，可推得盾构机向上掘进附加推力水平分量引起的地表土体变形计算公式[wq2]。
盾构向上掘进开挖面不均匀附加推力引起的地层变形计算公式为：
[w1=wq1+wq2] （5）
（2）盾壳摩阻力
盾构机向上掘进盾壳摩阻力引起的土体变形计算示意图如图4所示，图中[dA]为盾壳上任意微元体面积，面积大小为[dA=Rdldθ]，[R]为盾壳半径，[l]为该微元体所在位置沿盾壳至开挖面的轴向距离，盾壳摩阻力分解为竖向分力[fsinα]以及水平分力[fcosα]，则该微元体所受的集中力可以表示为[dPv+dPh=fsinαRdldθ+fcosαRdldθ]。
同理对坐标系进行转换可得：
[x’=x-Htanα+ly’=y-Rcosθz’=zc=Hcosα-Rsinθ] （6）
将式（6）代入式（1），盾壳摩阻力分为上下两侧同时积分，同理得到向上掘进盾壳摩阻力竖向分量引起的地表土体变形[wf1]的公式。
将式（6）代入式（2），盾壳摩阻力分为上下两侧同时积分得到向上掘进盾壳摩阻力水平分量引起的地表土体变形[wf2]的计算公式。
盾构向上掘进盾壳摩阻力引起的地层变形[w3]的计算公式为：
[w3=wf1+wf2] （7）
<G:\武汉工程大学\2024\第5期\王慧-4.tif>[x′ x″][X][Z][z″][z′ ][O][α][H][y][y″][y′ ][切口][dl][f1][f2][dA][l][盾尾][L][dθ]
图4 向上掘进盾壳摩阻力计算示意图
Fig. 4 Schematic diagram of friction resistance calculation of upward driving shield shell
（3）盾尾注浆压力
注浆压力[p]引起的地表土体变形的计算示意图如图5所示。图中[dA]为注浆液表面上任意微元体面积，[S]为注浆压力影响范围。由于管片外径小于盾构机外径，在盾构机盾尾脱出后管片与周围土体之间存在间隙，需要填充浆液来保证地层的稳定性。假定浆液沿管片表面均匀分布，注浆压力方向垂直于管片表面，可分为竖向和水平两个分量。已有的研究表明［15］盾尾注浆压力水平分量力引起的竖向位移较小，本文仅考虑注浆压力竖向分量引起的竖向位移。微元体表面受到的注浆压力如图5所示。在xoz平面注浆压力[p]可以分为竖向分力[p1=pcosα]与水平分力[p2=psinα]，在yoz平面又可把[p1]分为竖向分量[p3=p1sinθ]以及水平分量[p4=p1cosθ]，由于仅注浆压力考虑竖向分量，微元体受到的集中力为[dph= pcosαsinθRdsdθ]。同理对坐标系进行变换：
[x’=x-Htanα+L+sy’=y-Rcosθz’=zc=Hcosα-Rsinθ] （8）
将式（8）代入式（1）积分后得到盾尾注浆压力[w5]的计算公式。
<G:\武汉工程大学\2024\第5期\王慧-5.tif>[x′ x″][X][Z][z″][z′ ][α][H][y][y″][y′ ][dθ][地表][切口][盾尾][z″][z][S][o][y′ y″][y][o o′][p4][p1][p3][θ][z][dθ][p1][p][dA][ds][p2]
图5 向上掘进盾尾注浆压力计算示意图
Fig. 5 Schematic diagram of grouting pressure calculation at the end of the upward driving shield
1.3 向下掘进引起地表沉降
（1）开挖面附加推力
盾构机向下掘进时开挖面附加推力引起土体变形的计算示意图如图6所示。
同理可得盾构机向下掘进坐标系变换为：
[x=x’+Htanαy’=y’-rcosθz’=z’c=Hcosα-rsinθ] （9）
将式（9）代入式（1）后，开挖面附加推力分为上下两侧同时积分得到向下掘进附加推力竖向分量引起的地表土体变形[wq3]的计算公式。
将式（9）代入式（2）后，开挖面附加推力分为上下两侧同时积分得到向下掘进附加推力水平分量引起的地表土体变形[wq4]的计算公式。
盾构向下掘进开挖面不均匀附加推力引起的地层变形计算公式为：
[w2=wq3+wq4] （10）
（2）盾壳摩阻力
盾构机向下掘进盾壳摩阻力引起的土体变形的计算示意图如图7所示。
同理将坐标系变换为：
[x=x’+Htanα+ly’=y’-rcosθz’=z’c=Hcosα-rsinθ] （11）
将式（11）代入式（1），盾壳摩阻力分为上下两侧同时积分，得到盾构向下掘进盾壳摩阻力竖向分量引起的地表土体变形[wf3]的计算公式。
将式（11）代入式（2），盾壳摩阻力分为上下两侧同时积分，得到盾构向下掘进盾壳摩阻力水平分量引起的地表土体变形[wf4]的计算公式。
盾构向下掘进盾壳摩阻力引起的地层变形计算公式为：
[w4=wf3+wf4] （12）
（3）盾尾注浆压力
盾构机向下掘进注浆压力引起土体变形的计算示意图如图8所示。
<G:\武汉工程大学\2024\第5期\王慧-8.tif>[x′ x″][X][z″][α][H][y][y″][y′ ][盾尾][dθ][p1][z′ ][z ][切口][dA][o][地表][z ][z″][z′ ][S][y′ y″][y][o o′][θ][p3][p4][p1][p][p2][dθ][ds]
图8 向下掘进注浆压力计算示意图
Fig. 8 Schematic diagram of grouting pressure calculation in down-turning curve
同理将坐标系变换为：
[x’=x+Htanα+L+sy’=y-Rcosθz’=zc=Hcosα-Rsinθ] （13）
将式（13）代入式（1）后积分得到盾构向下掘进引起地表土体变形[w6]的计算公式。
1.4 地层损失引起的地表沉降
在盾构机曲线掘进的过程中，对内侧土体超挖和盾尾间隙会造成地层损失，C. Sagaseta推导出地层损失引起的地表竖向位移[w7]计算公式为：
[w7=Vloss2πHy’2+H21-x’x’2+y’2+H2] （14）
式（14）中[Vloss]为隧道单位长度的地层损失量。孙捷城等［13］在C. Sagaseta地层损失理论的基础上，推导出曲线盾构施工地层损失引起的地表竖向位移公式：
[w7=Vloss2πHy+0.5δ2+H2]
[1-x+Lx+L2+y+0.5δ2+H2] （15）
式中[δ]为线路上侧额外超挖量，由于本文为竖向曲线盾构，因此用[H-0.5δ]代替式（15）中的[H]，得到竖向曲线盾构地层损失引起的地表沉降公式为：
[w7=Vloss2πH-0.5δy2+H-0.5δ2]
[1-x+Lx+L2+y’2+H-0.5δ2] （16）
式中[Vloss=πR2-πR’2=πRg+δ-g+δ24]，根据赵丹［16］关于曲线盾构施工所需间隙的研究，线路上侧超挖量[δ]为：
[δ=12Rc-R-Rc-R2-L12] （17）
式（17）中：[Rc]为曲线盾构线路的曲线半径；[L1]为单段盾壳长度，考虑盾构机铰接装置的作用，取[L1=L/2]。
综上所述，将式（5）、式（7）、[w5]的计算公式、式（16）相加得到盾构向上掘进引起的总沉降变形计算公式为：
[w=w1+w3+w5+w7] （18）
将式（10）、式（12）、[w6]的计算公式、式（16）相加得到盾构向下掘进引起的总沉降变形计算公式为：
[w=w2+w4+w6+w7] （19）
2 计算结果及影响因素分析
以北京市通州区东六环路改造工程试掘进段为例进行计算分析。结合工程资料与施工现场监测数据，计算所需的各参数取值为：转弯半径[Rc]为1 500 m，隧道中心处埋深[H]为17.5 m，盾构开挖面半径为7.985 m，泊松比[μ]为0.28，摩阻力差异系数[n]为0.9，附加推力差异系数[m]为1.1，盾壳长度L为16 m，剪切模量G为22 MPa，附加推力[q]为100 kPa， 盾构转弯角度为3°，盾壳摩阻力[f]为80 kPa，盾尾注浆压力[p]为200 kPa，土体损失参数g为30 mm，注浆压力影响范围[S]为2 m。
2.1 横向地表变形
开挖面前10 m、开挖面后16 m以及开挖面后26 m处地表横向变形曲线图如图9所示。从图9（a）可以看出在盾构机到达监测断面前各施工荷载引起的地表变形均为隆起，地表隆起最大值在中心线处，向两侧逐渐减小。总地层变形表现为中轴线附近呈现隆起，两侧为沉降。轴线处土体隆起值为0.49 mm。图9（b）与图9（c）可以看出盾尾注浆压力在管片脱出盾尾的位置隆起值最大。整体来看，开挖面前土体隆起主要与盾壳摩阻力有关，开挖面后土体隆起主要是盾尾注浆压力造成的。附加推力与注浆压力对地表变形的影响较小，地表沉降主要是由地层损失导致的，其次是盾壳摩阻力。各断面的总地表变形数据与实测数据之间误差较小，变化曲线基本一致，验证了本文计算公式的可靠性。
<G:\武汉工程大学\2024\第5期\王慧-9-1.tif><G:\武汉工程大学\2024\第5期\王慧-9-2.tif><G:\武汉工程大学\2024\第5期\王慧-9-3.tif>[-20 -15 -10 -5 0 5 10 15 20
距中轴线距离 / m][2.0
1.5
1.0
0.5
0.0
-0.5
-1.0][沉降值 / mm][盾壳摩阻力
附加推力
注浆压力
地层损失
总地表变形
实测数据
][（a）][-20 -15 -10 -5 0 5 10 15 20
距中轴线距离 / m][盾壳摩阻力
附加推力
注浆压力
地层损失
总地表变形
实测数据
][（b）][-20 -10 0 10 20
距中轴线距离 / m][0
-1
-2
-3
-4
-5
-6
-7
][沉降值 / mm][盾壳摩阻力
附加推力
注浆压力
地层损失
总地表变形
实测数据
][（c）][0
-1
-2
-3
-4
-5
-6
][沉降值 / mm]
图9 横向地表变形曲线：（a）开挖面前10 m处，
（b）开挖面后16 m处，（c）开挖面后26 m处
Fig. 9 Transverse surface deformation curves：（a）10 m
before excavation face，（b） 16 m behind the excavation face，
（c） 26 m behind the excavation face
2.2 纵向地表变形
图10为隧道轴线处纵向地表变形曲线，从图10中可以看出开挖面附加推力与盾壳摩擦力引起的地表变形均呈反对称分布。盾尾注浆压力在注浆处隆起值最大，向两边逐渐减小。最大隆起值为0.29 mm。总地表竖向变形在刀盘前方呈现隆起，在刀盘后开始加速沉降，在刀盘后方20 m左右沉降开始放缓，随后曲线逐渐平稳。该理论计算值与现场实测值趋势基本吻合。
<G:\武汉工程大学\2024\第5期\王慧-10.tif>[-20 -10 0 10 20
距中轴线距离 / m][0
-1
-2
-3
-4
-5
-6
-7
][沉降值 / mm][盾壳摩阻力
附加推力
注浆压力
地层损失
总地表变形
实测数据
]
图10 隧道轴线处纵向地表变形曲线
Fig. 10 Vertical surface deformation curves at the tunnel axis
2.3 地表横向沉降影响因素分析
2.3.1 开挖面推力差异系数分析 图11为在x=10 m处，不同推力差异系数下开挖面不均匀附加推力引起的地表横向变形曲线，从图11中可知随着差异系数m的增大，断面处地表隆起值也随之增大，其中隧道中心处增加最为明显，距离中轴线越远，变化越小。差异系数m为1.1、1.3、1.5、2.0、3.0时地表土体最大隆起值分别为0.19、0.21、0.23、0.28、0.41 mm。
<G:\武汉工程大学\2024\第5期\王慧-11.tif>[0.4
0.2
0.0][沉降值 / mm][-20 -10 0 10 20
距中轴线距离 / m][m=1.1
m=1.3
m=1.5
m=2.0
m=3.0]
图11 不同开挖面推力差异系数横向地表沉降曲线
Fig. 11 Lateral surface settlement curve of different
excavation face thrust difference coefficient
2.3.2 盾壳摩阻力差异系数分析 图12为在x=
-16 m处，不同盾壳摩阻力差异系数下盾壳摩阻力引起的地表横向变形曲线，从图12中可以看出不同差异系数地表沉降最大值均在隧道中轴线处，向两侧逐渐减小。地表沉降值整体随着差异系数的减小而减小。差异系数n为0.9、0.8、0.6、0.4、0.2时地表土体最大沉降值分别为-1.80、-1.72、-1.56、
-1.40、0.41 mm。
<G:\武汉工程大学\2024\第5期\王慧-12.tif>[0.0
-0.5
-1.0
-1.5
-2.0][沉降值 / mm][-20 -10 0 10 20
距中轴线距离 / m][n=0.9
n=0.8
n=0.6
n=0.4
n=0.2]
图12 不同盾壳摩阻力差异系数横向地表沉降曲线
Fig. 12 Lateral surface settlement curves of different shield friction coefficient
2.3.3 地层损失分析 图13为在x=-26m处，不同转弯半径下地层损失引起的横向地表变形曲线，从图13中可以看出随着转弯半径的增大，地层损失引起的地表沉降逐渐减小，转弯半径[Rc]为1 000、1 300、1 500、1 700、2 000 m时地表最大沉降值为-4.93、-4.64、-4.50、-4.30、-4.16 mm。
<G:\武汉工程大学\2024\第5期\王慧-13.tif>[0
-1
-2
-3
-4
-5][沉降值 / mm][-20 -10 0 10 20
距中轴线距离 / m][Rc=1 000
Rc=1 300
Rc=1 500
Rc=1 700
Rc=2 000]
图13 不同转弯半径横向地表沉降分析
Fig. 13 Lateral surface settlement analysis with different turning radius
3 结 论
（1）考虑盾构竖向掘进转弯方向与转弯角度，建立竖向曲线盾构掘进计算模型，基于Mindlin解建立局部坐标系并对原有坐标系进行转换，推导出竖向曲线盾构施工引起的地表沉降计算公式，结合工程实例分析各施工参数对土体变形的影响规律，得到的计算结果验证了公式的可行性。
（2）公式计算结果表明，开挖面前土体隆起主要与盾壳摩阻力有关，开挖面后土体隆起主要是盾尾注浆压力造成的。地表沉降主要受盾壳摩阻力与地层损失的影响较大，其中地层损失的影响最大。隧道轴线处纵向地表变形曲线中开挖面附加推力与盾壳摩阻力引起的地表变形呈反对称分布，盾尾注浆压力在注浆处隆起值最大，向两边逐渐减小。总地表竖向变形在开挖面前方呈现隆起，在开挖面后方开始加速沉降，在开挖面后方20 m左右沉降开始放缓，随后曲线逐渐平稳。
（3）不同推力差异系数、摩阻力差异系数以及曲线半径引起的土体变形最大值均在隧道轴线处，其中推力差异系数和摩阻力差异系数大小均与地表沉降值呈正相关，曲线半径大小与地表沉降值呈负相关。