开挖引起的隧道位移动态多目标优化反演预测

何 维， 孙宏磊， 陶袁钦， 蔡袁强

(1. 浙江大学 建筑工程学院，杭州 310058； 2. 浙江工业大学 土木工程学院，杭州 310004)

控制基坑开挖导致的邻近既有隧道变形对设计和施工至关重要，因此有必要对基坑开挖引起的隧道变形进行预测.然而，基坑开挖时存在多种不确定性因素，开挖卸荷的影响范围广，使得传统的理论计算和数值模拟等正向分析方法存在众多局限，主要包括以下3个方面：① 基坑开挖时的土体参数难以确定[1]，采用现场勘察或室内实验所得土体参数计算得到的预测值往往和实际监测数据存在较大偏差；② 基坑开挖存在时间效应[2]，使得邻近隧道的位移发展滞后于开挖进度，这种时间效应不仅仅取决于软土的流变性质，还受开挖速率的影响，难以对其进行量化计算；③ 开挖卸荷影响范围内包含多种结构响应，传统正向分析方法难以对隧道、围护结构和土体变形等多目标响应同时进行准确预测.

文献[2]统计了上海软土基坑的监测资料，提出了基坑时间效应的修正公式.然而，经验公式在应用于新的工程时，参数取值的可靠性较低.而基于数据融合的反演分析方法[3-5]能够有效识别基坑开挖过程中的未知参数，提升预测模型的准确度，其中多目标优化算法[6]能够得到多个目标函数之间的“权衡”解，使得预测结果能够同时满足多目标响应的准确度要求.因此，本文采用多目标优化方法对基坑开挖中的隧道、围护结构等多目标响应进行预测研究.

在实际优化设计问题中，计算成本较低的代理模型常被用于近似替代计算昂贵的数值模型[7].然而，对于输入参数维度高、非线性强的模型近似，静态代理模型不能保证全局的近似精度，且静态代理模型的训练样本数随着输入参数的增加而快速增加[8].近年来，在优化进程中不断更新代理模型的动态优化方法得到了快速发展[9].这种动态代理模型方法只需要获得寻优路径上局部参数空间的高精度近似，减少了在其他参数空间中的数值模型计算.

工程优化中动态代理模型方法的计算效率有待进一步提高.加点准则是动态代理模型方法的核心技术，现有的多目标加点准则中样本更新点随着迭代次数增加而固定新增，不能识别无效的代理模型更新.在面向测试函数的多目标加点准则[10-11]中，全局最优解和Pareto前沿都是已知的，故常采用一些与最优解或解的目标函数值相关的指标作为优化收敛判别准则[12].相比测试函数，工程实例的全局最优解和Pareto前沿都是未知的，通常采用最大迭代次数作为优化收敛判别准则，为保证获得最优解，最大迭代次数往往设有冗余，此时使用不具备收敛判别能力的加点准则[13-14]将导致新增的样本更新点陷入解的多样性寻优空间中.这种多样性寻优空间中的样本更新并不能增强算法的寻优能力，反而很大程度上浪费了计算资源，限制了动态多目标优化方法在工程中的应用.

针对上述问题，本文基于多目标优化理论，融合基坑开挖中的多目标监测数据，反演识别了土体参数，量化修正了隧道位移的时间效应.同时为减少优化迭代时数值模型的计算次数，提出了一种基于自适应加点准则的动态多目标优化(Dynamic Multi-objective Optimization with Adaptive Infill Criterion，DMO-AIC)方法.首先通过两个常用的多目标优化测试函数，说明了方法的计算效益；再通过虚拟数值算例，验证了方法的多目标响应预测能力；最后将该方法应用于上海外滩596基坑工程案例，准确预测了基坑分步开挖引起的既有隧道竖向位移.

1 DMO-AIC方法

1.1 多目标优化的数学描述

一个经典的无约束、多目标优化问题可以表示为

min [f1(x)f2(x) …fk(x)]T

(1)

s.t.l≤x≤u

式中：l和u为参数向量x的寻优空间边界；k为目标函数f(x)的数目.

多目标优化问题所得是一个非劣解集或非支配解集，称为Pareto解集，其在目标空间的映射称为Pareto前沿，如图1所示.Pareto支配定义为假设x和x*是参数搜索空间中的两个可行解，对于∀i∈(1, 2, …,k)都有fi(x)≤fi(x*) ∧F(x)≠F(x*)，则称x支配x*.

图1 Pareto解在二维目标函数空间上的分布及映射Fig.1 Distribution and mapping of Pareto solutions on two-dimensional objective functions space

1.2 Kriging代理模型

代理模型常被用来在工程优化中近似计算昂贵的数值模型.常用的代理模型包括Kriging模型、多项式响应面及径向基函数等[7]，其中Kriging插值模型具有良好的非线性拟合能力，近似精度较高，基于Kriging的加点准则应用也较为广泛.

采用Kriging方法对数值模型进行近似，假设对于设计变量未知点x，其对应的函数响应值y(x) 可表示为随机函数：

(2)

(3)

式中：R(θ,x,x*)为与距离有关的相关函数；θ为函数的宽度控制参数.Kriging插值的公式推导过程和实现见文献[15].

1.3 动态多目标优化

动态优化是在优化进程中不断更新代理模型的优化方法，可分为动态单目标和多目标优化.动态多目标优化的框架如图2所示.

图2 动态多目标代理模型优化框架Fig.2 Framework of dynamic multi-objective optimization using surrogate models

步骤1基于全局抽样方法抽取一定数目的初始样本集，调用数值模型计算模型响应，训练初始代理模型.

步骤2构建多目标优化函数，融合监测数据，调用代理模型计算目标函数适应值.

步骤3基于多目标优化算法寻找目标函数的最优解.

步骤4基于相应的加点准则，得到最优的样本点并更新样本集，详见1.4节.

步骤5调用数值模型计算新样本点的模型响应，更新代理模型.

步骤6判别优化过程的收敛准则，当满足收敛准则时，停止优化并输出Pareto最优解，否则，重复步骤2～5.

1.4 基于距离的自适应加点准则

为减少优化过程中数值模型的计算次数，DMO-AIC方法在基于距离的加点准则[13]中加入了收敛判别策略，使得算法可以在种群陷入多样性寻优空间时停止更新代理模型，提升了多目标优化方法的计算效益及其对于工程优化的适用性.DMO-AIC的自适应加点准则包含3个策略：样本产生、收敛判别及样本选取.

策略1采用多目标粒子群(Multi-Objective Particle Swarm Optimization, MOPSO)算法[16]优化过程中每一代的当前非支配解作为候选样本集.当前非支配解由上一代解以及种群中其他粒子在经过局部飞行和全局变异两种机制后演化而来，兼具较强的局部搜索和全局搜索能力.MOPSO算法详见文献[16]，本文将重点介绍基于距离的收敛判别策略和样本选取策略.

策略2首先，定义某一代种群的非支配解粒子在目标空间的平均相邻距离：

(4)

式中：n为种群中的粒子总数；f(xi)为排序后的第i个非支配解粒子在目标空间的函数值.

(5)

(6)

(7)

根据式(6)，当α=1时，认为当前非支配解与上一代非支配解相比的进化效果能够满足预期，则根据一定的选点策略从当前的候选支配解选取新的代理模型训练样本点.当α<1时，则认为当前的候选非支配解中没有能够满足更新预期的粒子.

图3 最优样本更新点在目标函数空间中的位置Fig.3 Location of best new sample point on objective function space

(8)

为综合考虑两种寻优机制，定义一个相对参数空间最小距离DΩ，则第i个粒子的相对参数空间最小距离为

(9)

再定义一个相对目标空间最小距离D*，则第i个粒子的相对参数空间最小距离为

(10)

基于式(4)和(5)，被选取为代理模型新样本点的粒子将使得一个距离因子WD最大，如下：

WD=c1DΩ+c2D*

(11)

式中：c1和c2为两个权重系数，满足c1+c2=1.本文采用等权分配，即c1=c2=0.5.

2 DMO-AIC计算效益

为验证所提自适应加点准则的计算效益，用两个常用的已知最优解的多目标优化函数对DMO-AIC进行测试，并与无收敛判别策略的加点方法(简称DMO)的测试结果进行对比.

2.1 测试函数

测试函数ZDT2：

(12)

s.t. 0≤xi≤1,i=1, 2, …,n

ZDT2真实最优解集为

{x=[x1x2…xn]T|x1∈[0, 1],

(13)

xi=0,i=2, 3, …,n}

测试函数FON：

(14)

s.t. -4≤xi≤4,i=1, 2, …,n

FON真实最优解集为

{x=[x1x2…xn]T|x1=x2=…=

(15)

基于已知最优解，可以引入两个距离指标[17]对所得的优化解集进行收敛性和多样性的评价.其中，收敛性指的是优化得到的解与真实Pareto前沿的靠近程度，可通过优化解与Pareto前沿的样本点的最小距离DG来评价，计算如下：

(16)

式中：di为优化解中第i个个体与真实Pareto解集样本之间的最小目标空间距离.

多样性是指优化解的分布均匀性，多样性评价指标为

(17)

式中：df和ds为优化解的两端个体与真实Pareto解集样本之间的最小目标空间距离.

2.2 测试结果

对两个测试函数的主要优化参数设置如表1所示.Kriging代理模型的基函数和相关函数分别采用二项式函数和高斯指数函数.

表1 两个测试函数的优化参数设置Tab.1 Parameter settings for optimization of two test functions

为避免优化过程的偶然性，对每个测试函数进行3次优化计算并取其平均值.图4所示为两个测试函数的优化历史.使用DMO-AIC和DMO的收敛性、多样性寻优速度都表现较快.相比收敛性，优化解的多样性寻优需要更多的计算时间，通常在种群收敛性稳定后其多样性才能趋于稳定.

图4 使用DMO-AIC和DMO的优化历史比较Fig.4 Comparison of optimization histories by using DMO-AIC and DMO

表2所示为两个测试函数优化结果的关键评价指标.可以看到，使用DMO-AIC方法和使用DMO方法的最终收敛性、多样性结果接近，且都能满足计算精度要求.

表2 DMO-AIC和DMO的优化质量比较Tab.2 Comparison of optimization quality of DMO-AIC and DMO

结合表2和图4的计算结果可以发现，相比DMO，DMO-AIC寻优性能和收敛速度相近，但是使用DMO-AIC方法所需要的样本更新点平均数是30和97，而使用DMO所需要的样本更新点平均数是93和214，说明DOM-AIC的自适应加点准则能够在优化过程的多样性寻优阶段判断、筛选无效(进化效果差)的新样本点，从而提高计算效率.这种自适应加点准则在面对最优解未知的昂贵数值模型优化问题时，能够保证充分的寻优时间(较大迭代次数)，同时避免了数值模型的无效调用，有利于动态多目标优化方法在工程实践中的应用.

3 多目标响应更新

3.1 虚拟数值模型

利用虚拟数值模型来验证DMO-AIC对基坑开挖-既有隧道中多目标响应的更新和预测能力，在虚拟模型的模型响应值中加入合理的人工误差后，可认为最优解存在于模型参数真实值附近[18].同时该算例仅考虑单层土，输入参数较少，能较快寻优到其最优解.

虚拟模型如图5所示，整体宽30 m，深25 m.基坑采取一步开挖，宽5 m，深6 m，有限元模型采用对称设置.地连墙和顶管隧道结构采用梁单元模拟，均假设为弹性材料，弹性模量分别为20和30 GPa，泊松比均为0.2，墙体厚度为 0.6 m，管壁厚度为0.3 m，隧道内径为3 m. 地连墙以及隧道和土体之间采用绑定接触.支撑结构采用弹簧单元模拟，弹簧刚度上下分别设为12 和 50 MN/m.土体的重度为18 kN/m3，水平土压力系数K0=0.5.土体参数设置采用修正剑桥模型，参数取值如表3所示.数值模型边界约束水平位移，底部同时约束水平和竖向位移.

图5 虚拟数值模型的土层及结构示意图(m)Fig.5 Schematic of soil layer and structure of virtual numerical model (m)

文献[19]指出相对敏感度较高的修正剑桥模型的土体参数包括临界状态应力比M、回弹参数κ和泊松比ν.因此对土体的M，κ，ν这3个参数进行反演分析，其参数寻优范围见表3括号内取值.

表3 土体修正剑桥模型参数Tab.3 Soil parameters of modified Cam clay model

数值模型的参数优化需要对模型计算值和监测值进行融合，采用一个基于最小二乘法的目标函数对每个监测点的计算误差进行等权分配，并采用多个目标函数分别融合地连墙的侧移和隧道位移数据.目标函数值可根据下式计算：

(18)

式中：X为待优化的参数向量；N为该目标函数所利用的监测点个数；Dm、Um分别为第m个监测点的监测值和模型计算值.虚拟监测值通过在基于参数真实值计算得到的模型响应值中加入人工的高斯随机误差[18]来模拟，高斯随机误差e由仪器系统误差ε和高斯观测误差δ两部分组成：

e=ε+δ

(19)

式中：ε服从均值为0，标准差为ξU*的正态分布，ξ为系统误差系数，本文取0.01，U*为模型计算值；δ的均值为0，其标准差σδ取决于仪器精度，取仪器精度为±0.1 mm时，按95%的置信区间设置σδ为0.05.

3.2 多目标响应更新结果

为实现基坑开挖-隧道工程的多目标响应更新，基于式(17)，建立两个目标函数f1和f2，分别融合了地连墙和隧道的位移数据.优化过程的种群规模设为50，初始代理模型训练样本数为16，最大迭代次数为100，折减系数初始值为0.5.

为避免优化的偶然性，共进行3次优化计算，代理模型更新平均增加样本数约为12.图6给出了其中一次优化计算得到的Pareto最优解，解呈现凹的形态，且其在目标函数空间上的最大适应值在5%左右，误差较小，说明多目标函数之间的“权衡”较好[20]，因此将最优解集的每个粒子都输入到数值模型中计算模型响应再取其平均值作为位移更新的结果.

图7(a)和7(b)分别给出了基于3组Pareto最优解计算得到的地连墙侧向位移UL、隧道位移U的更新均值和误差条，图中Dep为地连墙的测点深度.可以看到地连墙和隧道的位移更新平均值与模拟监测值都吻合很好，且其误差条较窄，表明使用图6中的全部最优解能够对多目标响应同时进行准确更新.

图6 虚拟案例的Pareto最优解Fig.6 Pareto optimal solutions of virtual case

图7 地连墙及隧道位移的更新均值和误差条Fig.7 Updated mean value and error bar of wall deflection and tunnel displacement

表4所示为3组Pareto最优解得到的土体参数更新结果.更新后的土体参数均值与真实值基本一致，且更新后的土体参数标准差较小，说明Pareto最优解在参数空间中唯一存在于参数真实值附近.

表4 更新后的虚拟案例土体参数Tab.4 Updated soil parameters of virtual case

4 上海外滩596基坑工程

4.1 工程概况

上海外滩596基坑[21]整体分为S1-A区、S1-B区及S2-A区、S2-B区，工程具体开挖步及时间进度如表5所示.采用Abaqus 2018软件对基坑开挖过程进行模拟，二维数值模型如图8所示，土体采用莫尔-库仑模型，各土层参数如表6所示.地连墙和隧道衬砌结构采用梁单元模拟，假设为弹性材料，弹性模量分别为20和30 GPa，泊松比均为0.2，墙体等效厚度为0.8 m，衬砌结构厚度为0.3 m，隧道内径为6 m.文献[21]中未介绍支撑参数，而本算例目的为预测和更新下卧隧道的竖向位移，不预测基坑挡墙侧移，故支撑结构由等效位移约束替代模拟[22-23]，在挡墙上施加监测所得侧移值[21].模型两端施加水平方向位移约束，底部同时施加水平和竖直方向位移.

表5 上海外滩596基坑的开挖进程Tab.5 Excavation stages of Shanghai Bund 596 project

图8 上海外滩596项目的有限元模型(m)Fig.8 FE model of Shanghai Bund 596 project (m)

表6 土体莫尔-库仑模型参数Tab.6 Soil parameters of Mohr-Coulomb model

4.2 时间效应

文献[2]基于上海软土基坑开挖引起的下方隧道位移实测资料，提出了如下时间效应修正系数：

(20)

(21)

基于式(20)和式(21)，可得到隧道在不同开挖步时的累积位移：

(22)

图9所示为基于土体参数及时间效应参数初始值计算得到的隧道竖向位移(UV)，图中ES为开挖步，初始计算值与监测值有显著偏差，说明本文依据同类软土基坑的文献资料所决定的土体参数可靠性较低，且可推测隧道的实际位移相比计算值较为滞后.但位移计算值和监测值在不同的开挖区域施工时都有着相似的发展趋势，说明了有限元模型的合理性.

图9 使用模型参数初始值计算得到的隧道竖向位移Fig.9 Tunnel vertical displacement calculated by using initial values of model parameters

4.3 优化参数设置

实例计算中，采用DMO-AIC方法预测更新下卧隧道的竖向位移.考虑深厚土层的开挖对下卧隧道的影响较大，隧道位移可能随之快速发展，有必要对近期的隧道位移监测数据重点关注，因此多目标优化函数向量可设为F(X)=[f1(X)f2(X)]T，其中f1(X)、f2(X)分别为融合了全部开挖步数据和近期开挖步数据的目标函数.多目标优化函数所融合的监测数据构成如表7所示.

表7 两次更新计划所融合的监测数据Tab.7 Monitoring data merged from two update schemes

为识别关键的影响参数，降低未知参数维度，对上海外滩596基坑有限元模型进行参数敏感性分析，计算公式[19]如下：

(23)

式中：Spq为第p个开挖步的隧道位移对应于第q个输入参数的相对敏感性系数；Up为第p个开挖步的隧道位移值；ΔUp为第p个开挖步的隧道位移差值；xq为第q个参数的输入值；Δxq为第q个参数的差值.

图10 土体参数及时间效应参数的敏感度Fig.10 Sensitivity of soil parameters and time effect parameters

4.4 位移预测更新

相比单目标优化得到的唯一“最优解粒子”，多目标优化得到的最优解集能体现出多目标函数之间的“权衡”.为了避免工程案例的多目标最优解两端延伸过广，采用一个误差限制对最优解进行二次择优，该误差限制代表了决策者对各个目标函数的容许误差.

图11(a)和11(b)分别给出了两次更新计划I、II得到的Pareto解集.黑色粒子是优化计算得到的Pareto解，红色粒子表示对目标函数设置一定的“误差限制”后得到的二次择优解.解呈现出凹的形态，说明多目标函数间的“权衡”较好.但粒子在Pareto前沿上的分布比较不均匀，解的多样性较差，且由于工程案例最优解的未知性，端部粒子在目标函数f1上的误差值较大，为了提升预测准确度，有必要对目标函数设置一定的“误差限制”来对Pareto解集进行二次择优.本文案例中的“误差限制”可以设置为f1≤10%且f2尽可能小.

图11 两次更新计划下使用DMO-AIC得到的Pareto解集Fig.11 Pareto solutions of DMO-AIC in two updated schemes

图12所示为在两次更新计划下使用DMO-AIC方法的二次择优解计算得到的隧道位移预测均值和误差，可以看出，两次更新中后续开挖步的位移预测值与监测值都很接近.具体来说，更新计划 I 对近期开挖步10～14的监测数据赋予了较高权重，其对后续开挖步15～17的隧道位移预测值比更新计划 II 的计算值更为准确；而更新计划 II 对开挖步17的监测数据赋予了较高权重，其对后续开挖步18～20的隧道位移预测准确度高于更新计划 I，说明近期监测数据更能反映隧道位移的发展趋势.当隧道位移在某些开挖步突然快速发展时，需要及时融合近期的监测数据并更新预测.

图13所示为在更新计划 II 下使用DMO-AIC和DMO方法对隧道位移的预测准确度，图中Ta为所有开挖步完成后的后续施工时间.可以看出，二者对后续开挖步18～20的隧道位移预测误差接近且都较小.更新后的时间效应表明，隧道位移在开挖至坑底后的3个月内将会快速发展，其位移最大值为 25 mm 左右，远远大于开挖至基底时的12 mm监测值，说明如果基坑被暂停施工足够长时间，隧道位移将会有持续的较大发展，应及时进行后序的施工步骤.

图13 更新计划 II 下DMO-AIC、DMO的预测精度比较Fig.13 Comparison of prediction accuracy of DMO-AIC and DMO in Scheme II

表8所示为使用DMO-AIC和DMO方法对上海596基坑案例的计算效率.使用DMO-AIC所需的平均样本更新点个数为11和12个，使用DMO所需的平均样本更新点为72和79个，说明以最大迭代次数作为工程优化问题的收敛判别条件时，使用DMO-AIC方法能很好地避免产生过多的无效代理模型样本更新点，大大提高寻优的计算效率.

表8 使用DMO-AIC和DMO时两次更新计划的计算效率Tab.8 Calculation efficiency of two updated schemes by using DMO-AIC and DMO

5 结论

为提高工程优化中的动态代理模型方法的计算效率，提出一种基于自适应加点准则的动态多目标优化方法(DMO-AIC)，将其应用于上海外滩596基坑案例中，成功预测了基坑开挖引起的既有隧道位移，具体结论包括：

(1) 在两个多目标测试函数以及上海外滩596基坑案例的计算结果中，DMO-AIC和DMO方法都具有较高的计算精度，但DMO-AIC所需的样本更新点平均数分别仅为DMO的31%、45%和15%，有效避免了优化进程中原模型的无效调用.DMO-AIC的计算效率提升对于复杂工程优化问题具有重要意义.

(2) 虚拟基坑案例计算结果表明，最优解在目标函数空间上的误差值较小，表明最优解能较好满足多目标函数之间的“权衡”，此时使用所有最优解粒子能够对基坑开挖中的多目标响应同时进行准确更新.

(3) 以上海外滩596基坑为例，利用DMO-AIC方法分步融合了监测数据，成功准确预测了后续开挖步的邻近隧道竖向位移.其中，近期监测数据更能反映隧道位移的发展趋势，应对新的监测数据及时进行融合并赋予更高的权重.

(4) 上海外滩596基坑的更新后参数表明，隧道位移在开挖至坑底后的3个月内将会快速发展，最大值为25 mm，远大于开挖至基底时的12 mm监测值，说明该隧道位移的时间效应显著，在开挖至基底后应及时进行后序的施工步骤.所提时间效应修正模型及反演方法可为同类软土基坑开挖实践提供参考.

(5) 实际工程案例的Pareto最优解在目标函数空间上可能分布过广，通过对目标函数设置一定的“误差限制”可以对Pareto最优解二次择优，从而提升解的预测准确度.