软土地基具有承载力低、沉降大且持续发展等特性. 水泥土搅拌桩(deep mixing column, DMC)作为低成本地基加固的方法，已被广泛应用^[1-2]，设计优化可以合理分配资金，降低工程成本^[3]. 复合地基结构设计参数存在相互影响，地基土体具有变异性. 传统工程设计方法通过试算法不断调整设计参数组合，以验算承载能力极限状态(ULS)和正常使用极限状态(SLS)，在满足地基承载力及工后沉降要求的条件下，较难获得成本较低的优化方案^[4].

地基承载力可以依照规范估算，但利用经验或半经验方法难以考虑土体参数变异且精度有限，依据地基沉降计算结果确定的地基处理方案存在一定的安全风险. 有限元方法(FEM)常被作为复合地基沉降预测的重要手段，但需要设置合理的本构模型与参数. 修正剑桥模型(MCC)作为经典的弹塑性模型，可以很好地描述饱和正常固结或弱超固结黏土的基本特性，模型相对简单且参数较少^[5]. MCC的参数(临界状态应力比M、压缩指数 $\lambda $和回弹指数 $\kappa $)须通过固结不排水三轴剪切试验与各向等压固结试验测定，非常耗时. 塑性指数I_p与M、 $\lambda $和 $\kappa $具有密切联系^[6-7]，但拟合公式精度往往不高，忽略了土体参数的变异性^[8]. 人工神经网络ANN和Monte Carlo dropout (MCD)结合的ANN_MCD模型，可以在每次计算过程中随机切断固定比例的神经元连接，得到不同的模型架构，在随机计算次数后得到预测结果的高斯分布^[9]，实现采用贝叶斯方法揭示预测结果不确定性的目的^[10]. Zhang等^[11]将ANN_MCD模型应用于地基土体蠕变系数和渗透系数的不确定性预测，经FEM计算发现，模拟结果与实测值吻合良好.

通过FEM计算所有DMC复合地基结构参数与土性参数组合下的工后沉降，将耗费大量资源. 代理模型作为高效的数学近似方法，可以极大地提高计算效率^[12]. 常用的代理模型有Kriging模型、支持向量机模型、随机森林模型和神经网络(ANN)模型等^[13]. 其中ANN模型具有高容错性、并行分布和自学习等能力，对高维复杂的非线性关系具有较好的拟合能力. 采用ANN模型学习结构、土性参数与工后沉降量间的高维非线性映射关系，形成代理模型，快速预测不同结构参数组合下因土性参数不确定性导致的地基工后沉降概率分布. 采用地基工后沉降分布上不同置信度对应的界限值进行SLS验算，对于实际工程而言，对应的风险程度不同，应依据规范和工程等级确定.

综上所述，针对DMC复合地基工后沉降可靠性分析及设计优化的问题，采用ANN_MCD模型建立地基软黏土I_p 与MCC参数间的不确定性关系. 经多次ANN_MCD模型随机计算，获得符合高斯分布的统计数据，表征土体参数的变异性. 结合FEM与代理模型，高效预测地基土体与结构参数组合下的工后沉降. 以路基工程正常使用的极限状态设计可靠指标 $\;\beta $=1.65为标准，确定相同的结构参数组合下具有95%单侧置信度的工后沉降. 依据成本效能指标，以工程建造成本最低为目标，DMC复合地基的结构参数为优化变量，ULS和SLS为约束，开展设计方案筛选，实现考虑土体参数变异性的DMC复合地基设计优化.

针对如图1所示的加筋垫层DMC复合地基加固方案，建立分析模型^[14]. 加筋垫层由土、加筋材料、石灰及添加剂组成，具有较好的整体性，抗压和抗剪强度有一定的增强，能够有效地减小地基不均匀沉降和工后沉降量^[15]. 路基在横断面上为对称结构. 路堤高6 m，宽41.6 m，边坡坡率为1∶1.8. 地基表层1 m厚的回填土下为性质较差的2层深软黏土，共11 m深，下部分别为3 m厚硬黏土层及15 m厚黏砂土层. 地下水位(GWT)在地表以下1 m处.

加筋垫层DMC复合地基主要由可提高地基承载力、降低加固区压缩性的桩体和可减少地面不均匀沉降的垫层组成，地基处理成本主要由两者的材料费和相应施工设备、人工费构成. 依据现场试验^[16]，DMC桩按正方形布设，直径d = 1.2 m，桩长l_p = 10 m，桩间距a = 1.9 m；将加筋垫层铺设于路堤底部，厚h = 0.5 m.

为了避免遗漏最低建造成本的结构参数组合方案，结合地基条件和施工要求确定所有可能的设计参数可行域. 其中，桩长受最浅处理深度控制且以一定增量加长，桩间距、桩径和垫层厚度受荷载水平和施工工艺约束，导致设计参数为不连续变量. 基于建造成本的复合地基设计优化本质上为离散变量寻优问题^[17]，即以l_p、d、a、h参数为变量的结构优化. 当DMC长为12 m时，桩端处于下卧硬黏土层面；当DMC长为6 m时，可以保证桩体穿过土性最差的上软土层且进入性质相对较好的下软土层2 m. l_p宜取 6~12 m，增量为1 m；考虑施工设备和工程经验，d 取0.5~1.5 m，a取1.6~2.5 m，h取0.2~0.9 m，增量各取0.1 m^[14].

传统人工神经网络的权重、偏差值和模型结构是确定的，ANN_MCD可以按概率随机切断网络内一定比例神经元间的连接，形成不同的模型结构，输出符合高斯分布的随机值^[9]，如图2所示. MCD通常被作为模型训练中防止过拟合的手段，仅在训练过程中被激活，但ANN_MCD在模型训练与测试过程中均被激活. ANN_MCD网络内的权重和偏差虽然不变，但模型结构具有灵活性. 不同模型结构的输出结果不同，经过N次计算可得输出y = [y₁, y₂, ···, y_N]，最终输出结果的均值E (y)和方差Var (y)可以表示为

(1) $ {E}\;({\boldsymbol{y}})=\frac{1}{N} \sum_{i=1}^{N} {\boldsymbol{y}}_{i}\left({\boldsymbol{x}}, {\boldsymbol{W}}_{i}, {\boldsymbol{b}}_{i}\right) , $

(2) $ \operatorname{Var}\;({\boldsymbol{y}})=\frac{1}{N} \sum_{i=1}^{N}\left[{\boldsymbol{y}}_{i}^{{\rm{T}}}\left({\boldsymbol{x}}, {\boldsymbol{W}}_{i}, {\boldsymbol{b}}_{i}\right) {\boldsymbol{y}}_{i}\left({\boldsymbol{x}}, {\boldsymbol{W}}_{i}, {\boldsymbol{b}}_{i}\right)-{E}({\boldsymbol{y}})^{{\rm{T}}} {E}({\boldsymbol{y}})\right] . $

式中：x为输入变量，W_i和b_i分别为第i次计算中的权重和偏置，y_i为第i次模型的输出结果。

图1所示地基中的2层软黏土性质最差，是地基处理的主要区域，须对土体参数进行精细描述并考虑不确定性的影响. 采用修正剑桥模型(MCC)表达软黏土的力学行为，依据文献[16]，取压缩指数 $ {\lambda} $为回弹指数 $\kappa $的10倍. 仅选取土体临界状态应力比M和回弹指数 $\kappa $进行参数不确定性分析，该方法适用于更多不同参数的考量.

如图3所示，根据文献[8]共收集93组数据，其中 $ \kappa $样本共69组，M样本共24组；2个数据集的影响因素均是I_p. 将2个数据集中80%的样本用于ANN_MCD模型训练，20%的样本用于模型测试. 所有数据采用Min-Max方法进行归一化，消除各指标间的量纲影响. 采用10折交叉验证手段形成10个平行ANN_MCD模型，取平均值作为最终的输出结果.

采用网格搜索法，对ANN_MCD模型的隐藏层数(1, 2, ···， 5)、隐藏层神经元节点数(10，20，···，100)进行搜索，模型在隐藏层数为3、节点数为50时，在训练集和测试集上的预测效果表现较好. Dropout决定每次模型训练中随机切断神经元连接的数量，反映模型结构的灵活程度. 对模型在Dropout率为0.1~0.5的输入进行1 000次随机计算，以平均值作为输出结果，用R²表征训练集和测试集上的训练效果，如图4所示. 随着Dropout率的增大，神经元间更多的连接被切断，导致模型结构虽然更灵活，但映射能力减弱. 当Dropout率为0.1时，训练集和测试集的预测效果较稳定且R²较大. 如图3所示为ANN_MCD模型对应不同I_p时分别预测1 000次的结果，M和 $\kappa $的变异系数C均在8%~26%内波动，依据《地质灾害防治工程勘察规范》(DB 50/143—2003)可以划分为中低变异性等级，预测参数的95%置信区间与试验值吻合良好.

将加筋垫层DMC复合地基等效为平面应变问题^[18]，采用Plaxis 2D软件进行数值仿真，网格按高精度的15节点三角形单元进行划分，如图1所示. 基于等刚度(EA)法，将DMC等效为厚度为0.6 m的沿线路纵向布置的连续墙，桩长和桩间距保持不变，如图5所示. 根据对称性，取结构的一半进行建模，模型地基深度取30 m，长度为路堤底宽的3倍，约为80 m. 模型地基两侧约束法向位移且设置为不排水边界，下边界进行全约束并按排水边界考虑，上边界为自由排水面. 路堤填筑高度包括加筋垫层共6 m，按0.5 m 为1层，共分为12层进行填筑，每层填筑间隔时间平均约为8.3 d，共需100 d完成填筑^[14].

将DMC、路堤填土、加筋垫层和地基黏砂土视为理想弹塑性材料，采用摩尔库仑(MC)模型进行建模，地基表层回填土、硬黏土和2层软黏土采用MCC模型，材料参数如表1所示. 表中，E为弹性模量，v为泊松比，c为黏聚力，φ为内摩擦角，e₀为孔隙率，k_v为竖向渗透率，k_h为水平向渗透率. 考虑土性不确定性的2层软黏土参数M和 $ \kappa $，根据软黏土1和软黏土2土层范围内的I_p平均值(分别为57和30，均为碟式液限仪测定^[16])，采用2.2节训练的ANN_MCD模型进行1 000次蒙特卡洛随机计算，结果如图6、7所示. 图中，F为频数. 可知，软黏土1和软黏土2的M₁、M₂均值 $ \;{\mu} $分别为1.07、1.19，C分别为17.9%、8.2%；κ₁、κ₂的 $ \;{\mu} $分别为0.086、0.037，C分别为8.1%、9.6%. DMC与地基土界面设无厚度接触单元，剪切刚度由邻近土体性质决定，抗剪强度折减系数取0.8.

针对如图1所示的DMC复合地基结构参数，依据表1的材料参数开展FEM计算，得到地基沉降时程曲线，如图8所示. 图中，S、S_a分别为地基沉降的FEM计算值和实测值. 由于现场测试断面未设垫层，FEM计算中的垫层材料参数同路堤填土. 可知，采用ANN_MCD模型所预测的2层软黏土M和 $ \kappa $均值，获得的路堤填筑期(0~100 d)桩间土与桩顶位置沉降计算值与实测值几乎重合，填筑完成后的沉降发展趋势较吻合，尤其桩间土的工后沉降具有高度一致性，验证了建立的等效二维FEM模型适用于地基多层软土条件下的DMC复合地基模拟. 考虑基于ANN_MCD模型获得的土性参数不确定性，计算所得沉降分布的95%置信区间如图8所示. 可知，阴影区域基本覆盖了长期实测数据，吻合较好，表明利用ANN_MCD模型预测土性参数的有效性和准确性. 在计算过程中，以地基沉降小于0.1 mm/a为计算终止标准^[19]，地基固结时间取15 a.

根据1.1节确定结构参数l_p、d、a、h的可行域和增量，采用Monte Carlo方法对地基2层软黏土的M和 $ \kappa $随机参数进行1 000次抽样，得到如表2所示的结构和土性参数统计值. 所有参数可组合7(l_p)×8(d)×10(a)×8(h)×1 000 ( $ M_{1}, $ $ \kappa _{1}, $ $ M_{2}, $ $\kappa _{2} $) = 4.48×10⁶种形式，形成数据集. 若使用有限元逐个计算，则将消耗大量资源. 采用ANN模型建立DMC结构和加筋垫层尺寸、2层软黏土的M和 $ \kappa $与工后沉降的非线性映射关系，形成代理模型，快速预测其余样本结果. 采用正交实验设计方法，抽取297组典型样本进行FEM计算，但部分结构参数组合下的地基变形过大，有限元计算无法收敛. 剔除FEM无法计算的样本，将剩余175组样本用于ANN模型训练，其中随机抽取36组样本作为测试集，剩余的139组作为训练集.

结构参数l_p、d、a和h共同决定了DMC复合地基的加固稳定程度. 随着l_p、d和h的增大和a的减小，地基稳定性增强. 土体与结构间存在复杂的相互作用，可以将土体参数变异性视为正向或负向激励，对复合地基的整体稳定性产生影响. 如图9所示，设计ANN网络结构，结构参数和土性参数分别经过隐藏层进行特征提取，将提取特征交叉融合后，输出最终工后沉降预测值. 该结构较简单输入层的ANN模型更明晰、有意义，降低了模型的冗余性和特征提取的复杂程度，提高了模型的训练和计算效率.

采用网格搜索法，对ANN代理模型的最优结构进行搜索. 当ANN模型的隐藏层层数为2，每层的神经元节点数为32时，模型的学习效果较好. 模型采用Relu激活函数，使用Adma优化器，设置学习率为10⁻⁴，使用早停法(early stopping)监控验证集误差，连续100次都没有降低时停止模型训练. 使用10折交叉验证手段形成10个平行的代理模型，取平均值输出为最终的预测结果. 采用决定系数R²、均方根误差RMSE和平均绝对百分误差MAPE作为模型性能的评价指标，如下所示：

(3) $ R^{2}=1-\frac{\displaystyle\sum\limits_{i=1}^{n}\left(S_{i}-\hat{S}_{i}\right)^{2}}{\displaystyle\sum\limits_{i=1}^{n}\left(\hat{S}_{{i}}-\bar{S}_{i}\right)^{2}} , $

(4) $ \rm{R M S E}=\sqrt{\frac{1}{n} \sum_{i=1}^{n}\left(\hat{S}_{{i}}-S_{i}\right)^{2}} , $

(5) $ \mathrm{MAPE}=\frac{1}{n} \sum_{i=1}^{n}\left|1-\frac{\hat S_{i}}{S_{i}}\right| . $

式中： $\hat S_i$、 $ {\overline S_{i}} $分别为地基工后沉降预测值、均值，n为样本量.

如图10所示为两输入层ANN代理模型预测值与FEM计算值的对比情况. 训练集和测试集的R²均达到0.999，模型在整个结果空间内的预测精度普遍较高，说明该模型很好地学习到了结构、土体参数与地基工后沉降间的复杂非线性映射关系. 训练集的RMSE和MAPE分别为0.010 m和2.617%，测试集的RMSE和MAPE分别为0.014 m和3.868%，测试集的精度略低于训练集，说明模型没有过拟合，预测精度很高.

通过代理模型开展多参数敏感性分析，计算结构和土性参数对工后沉降的影响程度，即敏感度 $ \varepsilon _{i} $. 所有因素的 $ \varepsilon _{i} $总和为100%，因素i的 $ \varepsilon _{i} $越大，说明i对工后沉降的影响越显著. 如图11所示，结构参数对工后沉降的影响显著大于土性参数，l_p对工后沉降的影响最大，占65.3%；d的影响次之，为18.6%；a和h的影响接近，约为4.5%；土性参数的 $ \varepsilon _{i} $≤2.5%.

随着l_p、d和h的增加，地基工后沉降呈减小趋势，a与之相反. d和a对工后沉降的影响主要由面积置换率m = πd²/(4a²)控制，分别计算d和a对m的影响率绝对值，如下所示.

(6) $ \left|\frac{\partial m(d, a))}{\partial d}\right|=\frac{{\text{π}} d^{2}}{2 a^{2}} ; $

(7) $ \left|\frac{\partial m(d, a)}{\partial a}\right|=\frac{{\text{π}} d^{2}}{2 a^{3}} ; $

(8) $ \left|\frac{\partial m(d,a)/\partial d}{\partial m(d, a)/ \partial a}\right|=\frac{a}{d} > 1,\; a > d . $

从式(8)可知，d对m的影响大于a，所以d对工后沉降的影响相对a更显著.

每组结构参数对应的沉降分布计算步骤如下. 1）依据表3，选定复合地基的结构参数l_p、a、d和h的组合形式. 2）将对应土层I_p输入ANN_MCD模型，随机计算1 000次，得到土性参数的组合. 3）在结构参数不变的条件下，将1 000组土性参数循环输入ANN代理模型，得到对应结构参数组合下的复合地基工后沉降分布.

以图1所示的现场结构参数为例，工后沉降预测值的分布如图12所示. 经K-S检验可知，沉降的分布满足高斯分布. 当软黏土1和软黏土2的M₁、M₂变异系数分别为17.9%、8.2%，κ₁、κ₂的变异系数分别为8.1%、9.6%时，现场工况的工后沉降FEM计算值 $ \mu $= 0.230 m、 $ \sigma $= 0.020 8 m，工后沉降在较合理的范围内波动，说明ANN_MCD模型可以依赖大量的可靠试验数据，学习到较合理的参数变异性作为输出. 《铁路工程结构可靠性设计统一标准》(GB 50216—2019)^[20]要求铁路路基在正常使用极限状态下的可靠度指标 $ \;\beta $为1.0~2.5，以具有95%单侧置信度的工后沉降作为控制目标，取 $ \;\beta $ = 1.65. 如图12所示，现场工况对应的地基工后沉降S_p为 $\; \mu $ + 1.65 $ \sigma $= 0.265 m.

复合地基的优化问题通常有3类约束：现场条件、施工工艺和设计规范. 在2.1节，以现场条件和施工工艺作为线性约束条件，确定了结构参数l_p、a、d和h的可行域. 《复合地基技术规范》(GB/T 50783—2012)中规定了ULS (地基承载力、单桩承载力)和SLS (工后沉降)的非线性约束^[21].

路堤底面和桩端以下土体的地基承载力可以根据下式进行验算：

(9) $ p_{z}+p_{{\rm{c}} z} \leqslant f_{{\rm{a}} z} . $

式中：p_z为附加压力，p_cz为地基土自重压力，f_az为地基容许承载力特征值.

在路堤底面，桩与土体所提供的承载力按下式计算：

(10) $ f_{{\rm{a}} z}=\beta_{\mathrm{p}} m \frac{R_{{\rm{a}}}}{A_{\mathrm{p}}}+\beta_{\mathrm{s}}(1-m) f_{{\rm{sk}}} . $

式中：β_p和β_s分别为桩和桩间土承载力折减系数；m为面积置换率；A_p为单桩截面积；f_sk为桩间土承载力特征值； $ R_{{\rm{a}}} $为单桩承载力特征值，按下式计算取较小值：

(11a) $ R_{{\rm{a}}}=U_{{\rm{p}}} \sum_{i=1}^{n} q_{{\rm{s}}i} l_{i}+\alpha q_{{\rm{p}}} A_{{\rm{p}}} ,\tag{11a} $

(11b) $ R_{{\rm{a}}}=\eta f_{\rm{cu}} A_{{\rm{p}}} . \tag{11b}$

式中：f_cu为与DMC水泥土配比相同的室内加固土试块在标准养护条件下90 d龄期的立方体抗压强度平均值， $ \eta $为桩身强度折减系数，U_p为桩的周长，q_si为桩周第i层土侧侧阻力特征值，l_i为桩长范围内第i层土的厚度，q_p为桩端地基土未经修正的承载力特征值， $\alpha $为桩端天然地基土的承载力折减系数.

在复合地基处理范围以下存在软弱下卧层时，可以将桩体与桩间土视为等效基础，按式(9)进行下卧层承载力验算. 依据SLS的约束，应保证地基工后沉降S_p≤容许沉降量[ΔS].

DCM复合地基总建造成本Q由安装费和材料费共同组成，复合地基沿路堤纵向每米的平均建造成本可以表示为

(12) $ Q(x)=\left(c_{1} V_{{\rm{c}}}+nc_{2} l_{{\rm{p}}}\right) / a . $

式中：V_c为一个桩间距范围内垫层的体积；c₁为每立方米垫层的建造成本，包含1%纤维的费用^[15]；c₂为单根DMC的建造成本；n为沿路堤横断面的桩数. 根据《公路工程预算定额》(JTG/TB06-02-2007)^{[19, 22]}的规定可知，c₁= 20.8元/m³；当桩径为0.5 m时，l_p小于10 m的部分c₂ =330元/(10 m)，大于10 m的部分c₂=379元/(10 m)；当桩径 > 0.5 m时， c₂须乘以系数1.05^(d-0.5)/0.05. 考虑地基土体变异性的建造成本优化本质为寻找一组合理的结构参数，使其在一定可靠度下满足约束条件，同时令Q最小.

将表3中所有的结构参数组合按式(9)~(11)验算SLS约束条件，在满足要求后，按式(12)计算复合地基沿线路纵向每延米所需的建造成本. 现场结构参数组合下工后沉降具有95%单侧置信度的界限值为0.265 m，令[ΔS]=0.265 m进行优化. 如图13所示为符合约束条件的不同设计方案每延米建造成本与工后沉降间的关系，以现场结构对应的特征点为原点建立局部坐标系. 第Ⅰ象限内所有方案的建造成本和工后沉降量均大于原方案，位于Ⅰ象限的设计方案应避免. 第Ⅱ象限内所有方案的建造成本虽然低于原方案，但工后沉降较大，当建造成本有限且对工后沉降要求不高时，可以根据需求选择. 第Ⅲ象限内的方案建造成本和工后沉降均小于原方案，为优选区间. 第Ⅳ象限与第Ⅱ象限的情况相反，仅对沉降要求较严时予以考虑.

如图13所示，不同地基工后沉降对应最左侧的标识点为不同S_p对应的最佳结构参数组合. 采用Logistic曲线描述工后沉降随建造成本的变化趋势，成本效能分界点C_v ≈ 8 300万元/m对应Logistic曲线的曲率最大点. 依据C_v可以将曲线分为2个区段：建造成本＜C_v范围内为高效费比区段，工后沉降随着建造成本的增加而显著降低；建造成本＞C_v范围内为低费效比区段，随着工后沉降的进一步降低，需要投入的费用大幅增加，表明采用的设计方案与技术标准不匹配，宜选用其他类型的复合地基加固方案或线路结构形式.

与现场方案相比，第Ⅲ象限内不同工后沉降对应的优化方案共有20组，如表3所示. 从效费比最大化的角度考虑，全局优化方案应在C_v附近，对应表3中的方案15，结构参数l_p = 12 m，a = 2.5 m，d = 1.2 m，h = 0.9 m. 该方案的地基工后沉降为92 mm，仅需要投入7 691元/m；原方案的工后沉降为272 mm，建造成本却高达8 996元/m. 方案15与原方案相比，每延米建造成本约降低14.5%.

从表3可知，不同工后沉降对应的优化方案中l_p大多为12 m，使得DMC桩端处于硬土层. 较桩端处于软土层而言，可以显著地增大桩端承载能力，减小下卧层变形；在建造成本低且施工设备允许的情况下，应尽可能地增大l_p，控制工后沉降. 随着d的增大，地基工后沉降不断减小，在桩端处于良好持力层的条件下，宜优先考虑扩大桩径. h通常选择较大值，以增强加筋垫层的刚度，有助于垫层将路堤荷载更多地传递到DMC上，减少桩间土体所承担的荷载，减小土体变形. 考虑桩体的建造成本较高，应在满足ULS和SLS约束的条件下，尽可能选用较深桩长（端承桩），以控制工后沉降. 优化沉降敏感度相对较低的桩径和桩距可知，较小的桩径和较大的桩间距可以相对降低建造成本.

(1) 基于ANN_MCD模型结构的随机性，建立地基软黏土的塑性指数I_p与修正剑桥模型参数M和 $ \kappa $的不确定性关系. 经过 1 000次模型随机计算所得的M和 $ \kappa $统计参数符合高斯分布，95%置信区间的预测值与试验数据吻合良好.

(2) 通过双输入层ANN代理模型，对DMC复合地基的结构和土体参数独立进行特征提取，建立与地基工后沉降的非线性映射关系，可以有效地降低模型特征提取的复杂度，提升模型的训练和计算效率，实现地基工后沉降的高效、高精度预测，显著地减少FEM计算工作量.

(3) 地基工后沉降与最低建造成本呈Logistic曲线形式，依据成本效能分界值对应的曲率最大点，可以划分为高效费比和低效费比区段. 该算例表明，选用较深的桩长、较小的桩径和较大的桩间距，可以在满足ULS和SLS条件下有效地控制建造成本.