综合地表与深部位移监测数据的滑坡多目标加权位移反分析方法

戴 粤,戴吾蛟,余文坤

1. 中南大学测绘与遥感科学系，湖南 长沙 410083； 2. 湖南省精密工程测量与形变灾害监测重点实验室，湖南 长沙 410083

数值模拟计算结果的准确性与数值分析模型中材料分布、力学参数的定义和赋值息息相关。岩土力学参数统计特征一般通过开展室内试验、原位测试确定，但室内试验、原位测试都存在耗时长、数量少及费用昂贵等局限性，可获得的现场和室内试验数据往往非常有限[1]。自20世纪70年代以来，位移反分析方法作为岩土力学参数求解方法的一种补充，在岩土工程反演理论研究和工程实践中得到了极大的重视与发展[2-4]。传统方法需要通过大量的迭代计算进行求解，且存在容易陷入局部最优的缺陷。20世纪90年代起，提升全局寻优搜索能力的智能算法开始引入岩土参数反演中。文献[5]提出将神经网络用于位移反分析，可以解决传统方法计算量大、耗时长的问题。针对传统方法容易陷入局部最优的缺点，文献[6]将神经网络与遗传算法相结合，提出了一种位移反分析进化神经网络方法。文献[7]将支持向量机与遗传算法相结合，提出了一种位移反分析进化支持向量机方法，并验证了方法的有效性。随着人工智能技术的迅速发展，蚁群算法和粒子群优化算法等随机搜索方法也逐渐用于岩土参数反演[8-9]，其中粒子群优化算法具有收敛速度快、所需参数少、算法实现简单及全局搜索能力强等优点，近年来开始受到研究人员的重视[10-11]。

上述研究虽然大多针对算法容易陷入局部最优的缺点进行改进，但未认识到造成该问题的根本原因是当利用多点位或多类型变形监测数据进行反演，观测信息和待反演参数就会由相互冲突和影响的多个目标函数组成，即存在多目标优化问题(multi-objective optimization problem，MOP)。多目标优化问题的各个子目标之间往往是矛盾的，一个子目标的改善可能会引起其他子目标性能的降低，问题存在一组Pareto最优解集[12]，解并不唯一。因此，岩土参数反演已然需要面对和解决多目标优化问题。此外，对于滑坡位移反分析而言，深部位移信息量是否充足决定反演结果优劣，然而深部位移测点布设由于成本高昂，工程实际中往往点位布设较少。近年来，GNSS、InSAR和无人机视觉等观测技术迅速发展[13-15]，地表位移监测信息逐渐丰富，为该问题的解决提供了大量观测数据。

为解决以上问题，本文将多目标加权求和法用于岩土参数反演，并结合滑坡变形监测数据的特点，提出综合地表与深部位移监测数据构建滑坡多目标加权位移反分析模型，引入Helmert方差分量估计算法来优化加权位移反分析模型的权参数，并结合滑坡实例，采用对比分析与模拟粗差的方法验证了本文方法的有效性。

1 模型与方法

本文方法首先利用正交试验法计算岩土力学参数对位移的影响度，通过统计和分析确定位移反分析模型中的目标未知数；然后采用抗差Helmert方差分量估计方法计算不同类型滑坡数据的抗差验后随机模型，由此构建滑坡多目标加权位移反分析模型；最后通过正反分析迭代求解目标函数算得各等效力学参数。具体算法流程如图1所示。

图1 综合地表与深部位移监测数据的滑坡多目标加权位移反分析算法流程Fig.1 Algorithm flowchat of the proposed method

1.1 基于正交试验的岩土力学参数影响度分析法

利用正交试验法分析岩土力学材料参数对位移的影响度(敏感度)，在许多岩土工程问题中都取得了成功应用[16-17]，本文采用该方法来确定位移反分析模型中的目标未知数，即分析和剔除影响度小的参数，减少反演计算量。影响度的概念和计算步骤如下。

(1) 试验方案设计。选定能用于衡量试验效果的特征值作为试验指标，确定试验因素及其水平，以及根据因素和水平选用合理的正交表进行表头设计。以2水平3因素L4(23)正交表为例(表1)，列号对应为试验因素，列号下的数字对应为该因素的取值水平，表1中的试验是指给定一组材料参数所进行的一次数值模拟计算，计算位移值亦即试验指标。

表1 正交表L4(23)

(2) 执行方案，记录数据。正交表中的每一行代表一种水平组合，对每一组水平组合做一次试验，按照第k行的水平组合做第k次试验，所得到的位移计算值记为Xk，正交表若有t行，则需要做t次试验得到t个位移计算值。

(3) 列方差分析表，计算影响度。如果正交表安排了m个因素的试验，n次试验对应的试验结果记为C1,C2,…,Cn，因素水平数为j，每个水平做k次试验，对应的试验结果记为Cjk，方差分析计算步骤见表2。

表2 方差分析

总离差平方和ST、各因素离差平方和SM和试验误差的离差平方和Se的计算公式为

(1)

(2)

(3)

(4)

(4) 以特征点的计算位移值作为试验指标，岩土力学参数作为试验因素，根据上述步骤计算岩土力学参数对位移的影响度，影响度小于1的参数则不作为未知数参与反演计算。

1.2 基于抗差Helmert方差分量估计的多目标加权位移反分析模型

随着变形监测信息的逐渐丰富，利用多点位或多类型变形监测数据进行反演计算，必然涉及综合多种变形监测信息的多目标反演问题，多目标优化问题在其他工程应用中普遍存在，如投资、生产调度、物资运输等[12]。对于MOP问题，其解并不唯一，多目标加权求和法通过主观定义权重系数，能够将复杂的多目标优化问题转化为单目标问题，求得问题唯一解[19-20]，故本文提出将多目标加权求和法引入位移反分析中。

滑坡位移反分析的基础数据无非包括地表位移和深部位移监测数据两大类，现假设用于位移反分析的基础数据包含m类(n个)相互独立的观测值L=(L1L2…Lm)T，权矩阵P=diag(P1P2…Pm)=diag(p1p2…pn)，则对应的多目标函数模型即为

(5)

式中，Vobj为目标函数误差值，即正演计算位移值Si(X)与观测值Li之间误差的平方和；n和pi分别对应观测值的个数及其权重；X为目标未知数。

一般而言，给出的m类观测值的先验随机模型P不一定是恰当的。此外，当观测数据受到粗差污染时，该数据的权值就会跟原随机模型规定值不同[21]。文献[22—24]研究表明，方差分量估计方法可以解决不同类型观测量权比失调的问题，其有效性在数据处理领域已经得到广泛证实。为此，本文提出采用抗差Helmert方差分量估计方法来计算各类观测量的抗差验后随机模型，控制定权误差(包含粗差)对位移反分析计算结果的影响。现将抗差Helmert方差分量估计的迭代计算步骤归纳如下，为不失一般性，仍然将观测值分为m类，误差方程表示为

(6)

处理步骤如下：

(7)

(8)

由式(9)重新定权

(9)

2 试 验

2.1 滑坡数据

2.1.1 地质勘查数据

本文试验分析对象为位于常吉高速公路的湘西朱雀洞滑坡，该滑坡属于地质构造较简单的岩土层自然滑坡，预处理后的滑坡等高线地形图如图2所示，坡顶高程210 m，坡脚高程151 m，坡度范围为10°～20°。三维网格建模是数值模拟工作的前提，采用基于SURFER平台的FLAC3D三维地质建模方法建立该滑坡三维网格模型，模型如图3所示，计算范围采用的是中国地质科学院建议的范围，坡顶到底部边界距离设计为2H(H为滑坡高度)，坡底左右边界距离设计为3H。定义数值分析模型：该岩土层滑坡主要由填筑土层、粉砂岩层和基岩层组成，分析模型中岩土体的力学参数指标(初值)见表3，计算过程中将各地层定义为连续分布的材料，本构模型均采用摩尔库伦模型，边界条件除坡面设为自由边界外，模型底部和四周设置为固定边界约束，初始地应力场由分阶段弹塑性求解法生成。

表3 岩、土体力学参数指标

图2 滑坡地形Fig.2 Topographic survey data of landslide

图3 滑坡三维网格模型Fig.3 3D mesh model of landslide

2.1.2 监测数据

为掌握该滑坡的地表位移情况，项目布设了13个沉降监测点，于2008年4月2日—12月20日使用全站仪采用三角高程测量方法(对向观测)进行了25次沉降变形观测。需说明的是，由于施工原因，剔除了部分被破坏的点位，最终用于位移反分析监测点位，如图4所示。同时，本项目于2008年4—5月在边坡的关键部位布设了3个深部位移监测孔，点位布设情况见图4。表4、表5为各监测点原始数据和滤波后结果，滤波方法采用抗差Kalman滤波。

图4 监测点位布设示意Fig.4 Schematic diagram of monitoring point layout

表4 4月至12月地表位移观测数据

由表4可知，截至12月底，滑坡大部分点高程数据无变化，地表位移基本稳定。在数值模拟中可以认为在滑坡在天然状况下已趋于稳定，不再发生变形，即可利用12月的位移数据作为本次确定性位移反分析试验的基础数据，来反演滑坡等效力学参数。

2.2 正交试验分析

2.2.1 试验方案设计

首先，通过简单的试算可知，各岩土层材料参数剪胀角和抗拉强度以及基岩层材料参数对滑坡变形的影响可以忽略不计，可不作为试验因素参与分析，故试验方案中选取的试验因素为弹性模量E1、E2，泊松比v1、v2，黏结力c1、c2和内摩擦角φ1、φ2，密度ρ1、ρ2，下标为1的试验因素为填筑土层的岩土材料，下标为2的试验因素为粉砂岩层的岩土材料。其次，由于试验因素较多，不考虑内摩擦角和黏结力之间的交互作用，选用13因素3水平的正交表L27(313)进行表头设计，各因素水平按等间隔原则进行确定，见表6。最终试验设计方案见表7。

表5 测斜孔CX01-CX03累计位移数据(截至12月底)

表6 试验因素及其水平

表7 正交试验方案(表头设计)

2.2.2 试验结果分析

对于此分析模型，由表8可知，填筑土层和粉砂岩层的材料参数E1、v1、E2、v2、ρ1、ρ2对地表位移的影响度均大于1，材料参数E2、v2、c1、φ1对深部位移的影响度也均大于1，说明这些材料对地表位移(沉降)和深部位移均有显著影响，在进行位移反分析研究时可利用这些类型的位移观测量对该类岩土参数进行反演。而粉砂岩层的岩土材料c2和φ2对地表位移和深部位移的影响度均小于1，说明这两种参数对地表位移和深部位移均无显著影响，位移反分析研究时可不作为反演参数。

为避免偶然性，重新选取了10处节点(地表位移和深部位移特征点各5个)的位移值作为试验指标计算影响度，统计各材料参数影响度大于1的点位所占比例，如图5所示。不难发现，统计结果与表8结果基本一致，说明了上述分析结果的准确性，因此，根据正交试验分析结果最终确定位移反分析模型中的目标未知数为X=(E1,v1,E2,v2,c1,φ1,ρ1,ρ2)。

图5 影响度统计结果Fig.5 Statistical results of the influence degree

表8 影响度计算结果

2.3 位移反分析试验

为了验证本文方法的有效性和抗差性，设计两组试验分别进行对比分析。

第1组试验，位移反分析基础数据为10个地表位移监测点的沉降观测数据和3个测孔的深部位移数据。

方案Ⅰ：单类型输入变量加权位移反分析方法(仅深部位移数据作为输入向量)。

方案Ⅱ：多目标加权位移反分析方法，观测量的随机模型为初权模型。

方案Ⅲ，基于Helmert方差分量估计的多目标加权位移反分析方法，观测量的随机模型为Helmert方差分量估计验后模型。

第2组试验在第1组试验数据的基础上，随机选取5个地表位移监测点(DB01、DB03、DB05、DB07和DB09)人为加上6 mm的模拟粗差，验证本文方法的抗差性。

方案Ⅳ：基于Helmert方差分量估计的多目标加权位移反分析方法。

方案Ⅴ：基于抗差Helmert方差分量估计的加权位移反分析方法。

上述方案中，位移反分析模型中的目标未知数为3.2节分析结果：X=(E1,v1,E2,v2,c1,φ1,ρ1,ρ2)，各方案计算出的观测量的随机模型均列于表9，力学参数的搜索区间为

E1=400～600,E2=400～1400(MPa)

v1=0.2～0.4,v2=0.2～0.4

c1=16.0～26.0 kPa,φ1=32.0～42.0°

ρ1=1.90～2.10,ρ2=2.50～2.70(kg·cm-3)

搜索方法采用分层优化求解法[2]，计算程序采用C#语言编程实现，正反分析迭代计算至目标函数误差值最小时，各方案得到的等效力学参数(即反演结果)见表9。

表9 各方案反演结果

各方案正演计算位移值与实测位移值对比如图6—图9所示；由于地表位移与深部位移监测精度不同，同时深部位移监测精度也随孔深变化，故本文试验中选用平均绝对百分比误差(mean absolute percentage error，MAPE)来比较各方案反演结果的精度，误差统计结果见表9。

图6 实测位移值与各方案计算位移值对比(地表位移监测点)Fig.6 Displacement comparison between measured value and the calculated value of each scheme (for the surface displacement monitoring points)

图8 实测位移值与各方案计算位移值对比(地表位移监测点)Fig.8 Displacement comparison between measured value and the calculated value of each scheme (for the surface displacement monitoring points)

图9 实测位移值与各方案计算位移值对比(深部位移监测点)Fig.9 Displacement comparison between measured value and the calculated value of each scheme(for the deep displacement monitoring points)

2.3.1 第1组试验结果分析

(1) 比较方案Ⅰ与方案Ⅱ、方案Ⅲ计算结果，由图6、图7不难发现，仅利用深部位移监测数据来反演整个滑坡体的岩土等效力学参数，尽管由反演结果计算得到的深部位移与实测位移的吻合度很高，但此时地表位移计算结果却偏离于实际，说明当深部位移信息不足时，单输入变量加权位移反分析方法无法得到较好的岩土等效力学参数。而方案Ⅱ和方案Ⅲ相较方案Ⅰ而言，虽然深部位移计算精度有所损失，但方案Ⅱ和方案Ⅲ地表位移计算精度却能得到大幅提高(表8)，说明综合地表与深部位移信息的反演结果更为准确，能够有效解决深部位移信息量不足致使反演误差较大的问题。

(2) 对于方案Ⅱ和方案Ⅲ而言，方案Ⅲ深部位移计算精度与整体计算精度相较方案Ⅱ分别提高了32.6%和27.2%，而由表9观测值的随机模型可以看出，方案Ⅲ计算的验后随机模型中深部位移观测值的权重较其初始权重更大，说明方差分量估计出的观测值验后随机模型能更合理地描述加权位移反分析模型中不同类型观测数据应贡献的权重，即基于Helmert方差分量估计构建的加权位移反分析模型的准确性更优。

2.3.2 第2组试验结果分析

(1) 方案Ⅳ中，由图8、图9可以看出，各地表位移监测点和深部位移测孔的计算位移值与实测位移值均存在明显偏差。其中，地表位移计算精度为10.35%，DB08监测点绝对百分比误差高达27.86%；深部位移计算精度仅为44.64%，CX02测孔最大绝对百分比误差达到122.41%，精度均低于添加粗差前的计算结果(与方案Ⅲ比较，同一方法不同数据)，说明当位移反分析基础数据含有异常粗差时，基于Helmert方差分量估计的多目标加权位移反分析结果易受到粗差的污染。

(2) 对比方案Ⅳ和方案Ⅴ的计算结果，方案Ⅴ的计算精度均有大幅提升，计算位移值与实测位移值基本吻合，且方案Ⅴ反演得到的等效力学参数和反演结果精度与方案Ⅲ计算结果基本一致，说明基于抗差Helmert方差分量估计的多目标加权位移反分析方法能够有效抵制异常粗差对反演结果的影响。由图10可以看出，添加粗差的地表位移测点的权值在抗差Helmert方差分量估计后均迭代至零，正由于这些含有粗差的测点均被零权淘汰，才使得该方法具有良好的抵制粗差的能力。

图10 含有异常粗差监测点的权重迭代计算结果(方案Ⅴ)Fig.10 Weight iterative calculation results of monitoring points with abnormal gross errors (scheme Ⅴ)

值得说明的是，地表位移数据未添加粗差前，基于Helmert方差分量估计的多目标加权位移反分析方法计算位移值与实测位移值基本吻合，添加粗差后，地表位移计算结果没有明显受到异常粗差的影响，这是因为位移反分析方法的本质是一种基于最小二乘的误差寻优方法，部分地表位移计算误差会代入深部位移计算结果中(图11)，因此获取准确可靠的位移观测数据是位移反分析成功应用的关键。由图12可以看出，当观测数据含有异常粗差时，方案Ⅳ计算位移值较方案Ⅲ结果整体增大，且涨幅基本一致，这是由于建立滑坡数值分析模型时，将各岩土层材料定义为均匀且连续分布的理想情况导致的。尽管本文所提出的方法可以合理确定地表位移与深部位移监测数据的权重，且具有较好的抗差性能，仍然有部分点位计算结果存在偏差，除位移观测数据本身可能存在的误差外，岩土力学参数存在的空间变异特性也可能影响计算结果精度的原因。因此，将位移反分析方法扩展到顾及岩土力学参数空间变异性的情况也是后续研究的重要工作之一。

图11 方案Ⅲ、方案Ⅳ测斜孔计算结果偏差Fig.11 The calculation deviation of the scheme Ⅲ and scheme Ⅳ

3 结 论

(1) 对于滑坡岩土力学参数反演问题，深部位移信息量不足会导致位移反分析结果出现严重偏差，综合地表与深部位移监测信息进行反演计算，能够有效弥补滑坡深部位移监测点位稀疏的不足，提高计算结果精度。

(2) 本文解释了岩土参数反演存在的多目标问题并提出了应对方法。基于抗差Helmert方差分量估计的加权位移反分析方法考虑了不同类型观测数据在位移反分析模型中理应贡献的权重，以及观测数据中异常粗差对反演结果的影响。试验结果表明，无论是在求解精度还是在求解的抗差性上，该方法相较传统位移反分析方法均有较大幅度提升，可以更好地应用于滑坡岩土力学参数反演。

图12 实测位移值与方案Ⅲ、方案Ⅳ计算位移值比较(地表位移监测点)Fig.12 Displacement comparison between measured value and the calculated value of scheme -Ⅲ，scheme -Ⅳ (for the surface displacement monitoring points)