滑坡遥感识别（滑坡专刊戴粤）

本文内容来源于《测绘学报》2022年第10期（审图号GS京（2022）1079号）

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

戴粤^1,2

,戴吾蛟^1,2

,余文坤^1,2

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

基金项目：国家自然科学基金(41074004;42174053);湖南省自然科学基金(2021JJ30805);中南大学研究生创新项目(206021703)

摘要：针对滑坡参数反演存在的多目标优化问题, 同时为弥补滑坡深部位移监测点位稀疏的不足, 提出了一种综合地表与深部位移监测数据的滑坡多目标加权位移反分析方法。该方法利用地表与深部位移监测数据构建多目标位移反分析模型, 采用抗差Helmert方差分量估计法计算各类位移观测量的抗差验后随机模型, 优化反演模型的权参数。试验结果表明: 深部位移信息量不足会导致位移反分析结果出现严重偏差, 综合地表与深部位移信息的反演结果更为准确; 基于抗差Helmert方差分量估计的多目标加权位移反分析方法不仅可以合理确定不同类型观测数据的权重, 还能够有效抵制异常粗差对反演结果的影响, 提高了反演计算精度。

关键词：滑坡变形监测 多目标优化问题 加权位移反分析模型 抗差Helmert方差分量估计

引文格式：戴粤, 戴吾蛟, 余文坤. 综合地表与深部位移监测数据的滑坡多目标加权位移反分析方法[J]. 测绘学报，2022，51(10)：2149-2159. DOI: 10.11947/j.AGCS.2022.20210163

DAI Yue, DAI Wujiao, YU Wenkun. A landslide multi-objective weighted displacement back analysis method synthesizing ground and underground displacement monitoring data[J]. Acta Geodaetica et Cartographica Sinica, 2022, 51(10): 2149-2159. DOI: 10.11947/j.AGCS.2022.20210163

阅读全文：http://xb.chinasmp.com/article/2022/1001-1595/20221014.htm

引 言

数值模拟计算结果的准确性与数值分析模型中材料分布、力学参数的定义和赋值息息相关。岩土力学参数统计特征一般通过开展室内试验、原位测试确定，但室内试验、原位测试都存在耗时长、数量少及费用昂贵等局限性，可获得的现场和室内试验数据往往非常有限^[1]。自20世纪70年代以来，位移反分析方法作为岩土力学参数求解方法的一种补充，在岩土工程反演理论研究和工程实践中得到了极大的重视与发展^[2-4]。传统方法需要通过大量的迭代计算进行求解，且存在容易陷入局部最优的缺陷。20世纪90年代起，提升全局寻优搜索能力的智能算法开始引入岩土参数反演中。提出将神经网络用于位移反分析，可以解决传统方法计算量大、耗时长的问题。针对传统方法容易陷入局部最优的缺点，将神经网络与遗传算法相结合，提出了一种位移反分析进化神经网络方法。将支持向量机与遗传算法相结合，提出了一种位移反分析进化支持向量机方法，并验证了方法的有效性。随着人工智能技术的迅速发展，蚁群算法和粒子群优化算法等随机搜索方法也逐渐用于岩土参数反演^[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因素L₄(2³)正交表为例(表 1)，列号对应为试验因素，列号下的数字对应为该因素的取值水平，表 1中的试验是指给定一组材料参数所进行的一次数值模拟计算，计算位移值亦即试验指标。

表 1 正交表L₄(2³)Tab. 1 Orthogonal table for three factors with two levels

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

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

表 2 方差分析Tab. 2 The formula of variance analysis

总离差平方和S_T、各因素离差平方和S_M和试验误差的离差平方和S_e的计算公式为

(1)

(2)

(3)

根据方差分析结果计算各因素统计量

，对于某因素，给定显著水平α从F分布表查找对应的分位数F1－α(f_M, f_e)，若对于某因素有F(f_M, f_e)≥F1－α(f_M, f_e)成立，则该因素对试验结果的影响作用显著，因此，影响度计算公式定义为^[18]

(4)

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

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

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

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

(5)

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

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

(6)

处理步骤如下：

(1) 按经验或公式计算各类观测值的初始权矩阵P₁,P₂, …,Pm，进行第1次平差，求得各类改正数的平方和Vi^TPiVi；

(2) 按式(7)进行第1次(第k次)方差分量估计，求得各类观测值单位权方差的第1次估值

(7)

式中，

_ii=n_i－2tr(N^－1N_i) tr(N^－1N_i)²，n_i为该类观测值个数，N_i=Bi^TPiBi，N=N₁ N₂ … N_i，非对角线元素Sij=tr(N^－1NiN^－1Nj)，tr(A)表示为矩阵A的迹;W=[V₁^TP₁V₁V₂^TP₂V₂…Vm^TPmVm]^T为各类观测值改正数的平方和;

(3) 采用IGGⅢ方案计算观测值抗差等价权矩阵Pi(k)

(8)

由式(9)重新定权

(9)

式中，以标准化残差

^[21]，其他误差判别统计量的计算方法可参考；IGGⅢ方案中降权因子k₀一般取1.0~1.5，淘汰因子k₁一般取2.0~3.0；c为任一常数。

(4) 重复步骤(2)—步骤(3)，直至

2 试验

2.1 滑坡数据

2.1.1 地质勘查数据

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

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

图选项

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

图选项

表 3 岩、土体力学参数指标Tab. 3 Mechanical indicators of rock and soil(Original value)

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月地表位移观测数据Tab. 4 Surface displacement data from April to December

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

Tab. 5 Cumulative displacement data of inclinometer CX01-CX03 up to the end of December

mm

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

2.2 正交试验分析

2.2.1 试验方案设计

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

表 6 试验因素及其水平Tab. 6 Experimental factors and levels

表 7 正交试验方案(表头设计)Tab. 7 Header design in orthogonal experiment scheme

2.2.2 试验结果分析

根据正交表需要安排27次FLAC3D数值计算试验，在该滑坡模型中任意选择两处节点10 205和18 253(数字表示为节点ID号)，以节点10 205的竖向位移和节点18 253的水平位移作为试验指标，各组数值试验结果见表 7，单位为mm。按照2.1节中所述的正交试验法对表 7的数据进行计算，各参数对位移的影响度的计算结果见表 8，计算过程中给定显著水平α=0.05。以泊松比v₁对深部位移的影响度分析为例：①根据式(1)—式(3)计算总离差平方和ST=14.856，各因素离差平方和SM和试验误差的离差平方和Se=0.393，其中弹性模量v₁的离差平方和Sv₁=0.256；②由表 2公式计算因素v₁和试验误差的自由度fv₁=2、f_e=6，及其对应的平均离差平方和Sv₁=0.128、SE=0.065，计算得到因素v₁的统计量F(fv₁, f_e)=1.95；③查询分位数表可知F_1－0.05(fv₁, f_e)=5.14，根据式(4)计算得到泊松比v₁对深部位移的影响度I(v₁)=0.38。按照上述步骤同样计算其他参数的影响度，计算结果见表 8。

表 8 影响度计算结果Tab. 8 Results of the influence degree

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

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

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

图选项

2.3 位移反分析试验

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

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

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

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

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

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

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

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

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

表 9 各方案反演结果Tab. 9 Inversion results of each scheme

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

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

注：为将各方案计算位移值绘制在一幅图中进行比对分析，特将方案Ⅱ、Ⅲ计算结果依次增加常数3 mm。图 6 实测位移值与各方案计算位移值对比(地表位移监测点)Fig. 6 Displacement comparison between measured value and the calculated value of each scheme (for the surface displacement monitoring points)

图选项

图 7 实测位移值与各方案计算位移值对比(深部位移监测点)Fig. 7 Displacement comparison between measured value and the calculated value of each scheme (for the deep 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 Ⅳ

图选项

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

图选项

3 结论

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

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

作者简介

第一作者简介：戴粤(1997—), 男, 博士生, 研究方向为滑坡监测预警。E-mail: yuedai@csu.edu.cn通信作者：戴吾蛟, E-mail: daiwujiao@163.com

初审：张艳玲

复审：宋启凡

终审：金 君

资讯

,