一种利用数据同化改进的滑坡位移预测方法与流程

本发明涉及滑坡位移预测领域，尤其涉及一种利用数据同化改进的滑坡位移预测方法。

背景技术：

目前，滑坡位移预测模型主要利用智能算法、动力学方程和数值模拟等方法实现。其中，基于智能算法的方法主要利用历史监测数据对模型进行训练，类似一个“黑箱子”，只关心其输入和输出，难以通过物理力学原理描述和解释滑坡的演化过程和规律，而且其稳定性不高，容易在外界影响因素等输入中包含大量噪声时做出不符合物理力学的激烈响应，导致其适用性也较差。传统的基于动力学方程和数值模拟的方法主要被局限于应用在特定工况条件下的滑坡稳定性分析，因为该类方法是对滑坡演化过程的简化模拟，其精度受限且难以实现动态运行。

技术实现要素：

为了解决上述问题，本发明提供了一种利用数据同化改进的滑坡位移预测方法；一种利用数据同化改进的滑坡位移预测方法，主要包括以下步骤：

s101：根据目标滑坡的历史滑带完整性指标数据确定目标滑坡的滑带完整性指标值范围和滑带完整性指标初始值；

s102：获取目标滑坡第t个月的月降雨量rt和库水位月下降值kt，并根据外界影响因素-滑带完整性指标回归模型计算得到第t个月的滑带完整性指标st；

s103：采用条分法将目标滑坡条分为多个滑块，并根据布西尼斯克方程计算目标滑坡浸润线的位置，进而根据监测点所在的滑块与滑坡浸润线的相对位置确定监测点所在滑块的滑动面类型；所述滑动面类型包括天然滑动面和饱和滑动面；

s104：根据所述滑块的滑动面类型和第t个月的滑带完整性指标st，采用滑带抗剪强度动态变化指标公式计算得到第t个月的粘聚力ct和内摩擦角

s105：采用剩余推力法计算得到监测点所在滑块第t个月的内部推力ft；

s106：根据第t个月的剩余推力ft，采用推力-位移回归模型计算得到第t个月的滑坡累计位移值dt；

s107：根据所述滑带完整性指标值范围和第t个月的滑带完整性指标st，采用粒子滤波算法计算获得第t个月的最优滑带完整性指标和最优累积位移预测值；

s108：判断条件t＝sum是否成立？若是，则到步骤s109；否则，将t更新为t+1，返回步骤s102；其中，sum为预设的预测总月数，且sum>0；

s109：结束滑坡位移预测程序，得到各个月的最优累积位移预测值。

进一步地，步骤s102中，根据外界影响因素-滑带完整性指标回归模型计算得到第t个月的滑带完整性指标st的计算公式如公式(1)所示：

st＝α0+α1rt+α2rtkt+α3kt+α4st-1(1)

上式中，st-1为第t-1个月的最优滑带完整性指标，a0、a1、a2、a3和a4是对应的回归系数，均为预设值；t的初始值为1，当t＝1时，st-1的值为所述滑带完整性指标初始值。

进一步地，步骤s103中，根据监测点所在的滑块与滑坡浸润线的相对位置确定监测点所在滑块的滑动面类型；具体为：当监测点所在的滑块的滑带岩土体位于滑坡浸润线以下时，滑动面属于饱和滑动面；当监测点所在的滑块的滑带岩土体位于滑坡浸润线以上时，滑动面属于天然滑动面。

进一步地，步骤s104中，根据所述滑块的滑动面类型和第t个月的滑带完整性指标st，采用滑带抗剪强度动态变化指标公式计算得到第t个月的粘聚力ct和内摩擦角的计算公式如公式(2)所示：

上式中，cf和是峰值粘聚力和峰值内摩擦角，cr和是残余粘聚力和残余内摩擦角；若为饱和滑动面类型，则cf、cr、和分别为饱和峰值粘聚力、饱和残余粘聚力、饱和峰值内摩擦角及饱和残余内摩擦角；若为天然滑动面类型，则cf、cr、和分别为天然峰值粘聚力、天然残余粘聚力、天然峰值内摩擦角及天然残余内摩擦角；其中，两种类型的cf、cr、和均为从历史数据中获得的已知数据。

进一步地，步骤s106中，根据第t个月的剩余推力ft，采用推力-位移回归模型计算得到第t个月的滑坡累计位移值dt，计算公式如公式(3)所示：

dt＝β0+β1ft-1+β2dt-1+β3ft(3)

上式中，ft-1为第t-1个月的内部推力，dt-1为第t-1个月的累计位移值；β0、β1、β2和β3为对应的回归系数，且均为预设值。

进一步地，步骤s108中，根据所述滑带完整性指标值范围和第t个月的滑带完整性指标st，采用粒子滤波算法计算获得第t个月的最优滑带完整性指标；具体包括：

s201：在初始阶段，即当t＝1时，以所述滑带完整性指标初始值为均值初始化一组满足高斯分布的粒子；当t大于或者等于2时，不再执行步骤s201，直接到步骤s202；

s202：将每个粒子作为一个滑带完整性指标，并依次采用公式(1)、公式(2)和公式(3)相同的方法计算每个粒子所代表的滑带完整性指标在第t个月所对应的累计位移预测值；

s203：根据步骤s202中计算得到的各累计位移预测值和第t个月滑坡监测点的位移检测值，采用似然函数计算得到每个粒子所对应的权重，粒子权重值越大，则其所对应的滑带完整性指标与滑坡的真实情况越相似；

s204：将各粒子对应的累计位移预测值和其对应的权重进行加权求和，得到累计位移的最优估计，即第t个月的最优累计位移预测值；同时对相应的滑带完整性指标进行加权求和，得到滑带完整性指标的最优估计值，即第t个月的最优滑带完整性指标；

s205：利用重采样方法处理所有粒子，减少小权重的粒子，增加大权重的粒子，以避免对小权重的粒子进行无效的计算。

本发明提供的技术方案带来的有益效果是：本发明所提出的技术方案具备以下优点：

1)数值模拟是在滑坡的几何信息和土体参数等信息上进行的，是基于动力学方程进行迭代计算的，可以保证滑坡预测模型不会对外界影响因素的剧烈噪声做出过激响应；

2)针对滑坡位移预测模型存在结构误差和参数误差的问题，基于数据同化理论，利用粒子滤波算法从新的监测数据中提取有效信息，逐步优化模型状态和校正关键参数，可以在很大程度上提高滑坡位移预测模型的预测精度，并且能有效地让滑坡位移预测模型的预测误差收敛。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1是本发明实施例中一种利用数据同化改进的滑坡位移预测方法的流程图；

图2是本发明实施例中浸润线计算模型图；

图3(a)是本发明实施例中滑坡剖面图；

图3(b)是本发明实施例中利用剩余推力法对滑块进行受力分析图；

图4是本发明实施例中白水河滑坡条分效果图；

图5是本发明实施例中白水河动态抗剪强度参数的示意图。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。

本发明的实施例提供了一种利用数据同化改进的滑坡位移方法；应用于如图4所示的滑坡中，所述滑坡上设置有滑坡数据监测点；

请参考图1，图1是本发明实施例中一种利用数据同化改进的滑坡位移预测方法的流程图，具体包括如下步骤：

s101：根据目标滑坡的历史滑带完整性指标数据确定目标滑坡的滑带完整性指标值范围和滑带完整性指标初始值；

s102：获取目标滑坡第t个月的月降雨量rt和库水位月下降值kt，并根据外界影响因素-滑带完整性指标回归模型计算得到第t个月的滑带完整性指标st；

s103：采用条分法将目标滑坡条分为多个滑块，并根据布西尼斯克方程计算目标滑坡浸润线的位置，进而根据监测点所在的滑块与滑坡浸润线的相对位置确定监测点所在滑块的滑动面类型；所述滑动面类型包括天然滑动面和饱和滑动面；

s104：根据所述滑块的滑动面类型和第t个月的滑带完整性指标st，采用滑带抗剪强度动态变化指标公式计算得到第t个月的粘聚力ct和内摩擦角

s105：采用剩余推力法计算得到监测点所在滑块第t个月的内部推力ft；

s106：根据第t个月的剩余推力ft，采用推力-位移回归模型计算得到第t个月的滑坡累计位移值dt；

s107：根据所述滑带完整性指标值范围和第t个月的滑带完整性指标st，采用粒子滤波算法计算获得第t个月的最优滑带完整性指标和最优累积位移预测值；

s108：判断条件t＝sum是否成立？若是，则到步骤s109；否则，将t更新为t+1，返回步骤s102；其中，sum为预设的预测总月数，且sum＞0；

s109：结束滑坡位移预测程序，得到各个月的最优累积位移预测值。

步骤s102中，根据外界影响因素-滑带完整性指标回归模型计算得到第t个月的滑带完整性指标st的计算公式如公式(1)所示：

st＝α0+α1rt+α2rtkt+α3kt+α4st-1(1)

上式中，st-1为第t-1个月的最优滑带完整性指标，a0、a1、a2、a3和a4是对应的回归系数，均为预设值；t的初始值为1，当t＝1时，st-1的值为所述滑带完整性指标初始值。

步骤s103中，根据监测点所在的滑块与滑坡浸润线的相对位置确定监测点所在滑块的滑动面类型；具体为：当监测点所在的滑块的滑带岩土体位于滑坡浸润线以下时，滑动面属于饱和滑动面；当监测点所在的滑块的滑带岩土体位于滑坡浸润线以上时，滑动面属于天然滑动面。

步骤s104中，根据所述滑块的滑动面类型和第t个月的滑带完整性指标st，采用滑带抗剪强度动态变化指标公式计算得到第t个月的粘聚力ct和内摩擦角的计算公式如公式(2)所示：

上式中，cf和是峰值粘聚力和峰值内摩擦角，cr和是残余粘聚力和残余内摩擦角；若为饱和滑动面类型，则cf、cr、和分别为饱和峰值粘聚力、饱和残余粘聚力、饱和峰值内摩擦角及饱和残余内摩擦角；若为天然滑动面类型，则cf、cr、和分别为天然峰值粘聚力、天然残余粘聚力、天然峰值内摩擦角及天然残余内摩擦角；其中，两种类型的cf、cr、和均为从历史数据中获得的已知数据。

步骤s106中，根据第t个月的剩余推力ft，采用推力-位移回归模型计算得到第t个月的滑坡累计位移值dt，计算公式如公式(3)所示：

dt＝β0+β1ft-1+β2dt-1+β3ft(3)

上式中，ft-1为第t-1个月的内部推力，dt-1为第t-1个月的累计位移值；β0、β1、β2和β3为对应的回归系数，且均为预设值。

步骤s108中，根据所述滑带完整性指标值范围和第t个月的滑带完整性指标st，采用粒子滤波算法计算获得第t个月的最优滑带完整性指标；具体包括：

s201：在初始阶段，即当t＝1时，以所述滑带完整性指标初始值为均值初始化一组满足高斯分布的粒子；当t大于或者等于2时，不再执行步骤s201，直接到步骤s202；

s202：将每个粒子作为一个滑带完整性指标，并依次采用公式(1)、公式(2)和公式(3)相同的方法计算每个粒子所代表的滑带完整性指标在第t个月所对应的累计位移预测值；

s203：根据步骤s202中计算得到的各累计位移预测值和第t个月滑坡监测点的位移检测值，采用似然函数计算得到每个粒子所对应的权重，粒子权重值越大，则其所对应的滑带完整性指标与滑坡的真实情况越相似；

s204：将各粒子对应的累计位移预测值和其对应的权重进行加权求和，得到累计位移的最优估计，即第t个月的最优累计位移预测值；同时对相应的滑带完整性指标进行加权求和，得到滑带完整性指标的最优估计值，即第t个月的最优滑带完整性指标；

s205：利用重采样方法处理所有粒子，减少小权重的粒子，增加大权重的粒子，以避免对小权重的粒子进行无效的计算。

布西尼斯克方程计算目标滑坡浸润线的位置，可参考文献：[易贤龙.降雨与库水位作用下白水河滑坡渐进破坏概率研究[d].中国地质大学,2016]；浸润线计算模型如图2所示；可以得到计算滑坡浸润线位置的简化公式，如公式(4)：

上式中，h(x,t)是x位置的浸润线在t时刻的高度，h(0,0)是库水位变化前距离隔水底板的距离，v0是库水位上升或下降的速率，μ是给水度，w是降雨强度，k是含水层的渗透系数，hm是含水层的平均厚度。

剩余推力法(residualthrustmethod，rtm)是极限平衡法中的一种实现，因其计算简单、方便、有效，得到了广泛的应用。剩余推力法基于条分法进行滑坡条块的推力计算，主要应用于推测和分析滑坡的稳定性系数。滑坡的防护和治理工程通常根据稳定性系数来分析滑坡是否处于稳定状态，为工程的设计和建设提供依据。除此之外，剩余推力法还可以用于计算和分析滑坡体的内部受力分布情况。

运用剩余推力法时，需要知道滑坡剖面的几何形状，并且滑坡表面和滑动面可以近似表示为若干平面组成的折线。剩余推力法还做以下基本假定：假定滑坡体划分的每个条块单元均为刚体，自身的互相挤压和摩擦不计；滑坡体的每个条块单元整体沿折线形滑动面下滑；把条块单元的分析近似为平面应变问题；假定每个条块单元对下一条块单元的不平衡推力在分界面处的平行于该条块的滑动面。

假设滑坡的剖面如图3(a)，利用条分法将滑坡体近似等分为n块。利用剩余推力法对第i个条块进行受力分析，如图3(b)。

基于上述的假设，结合第i个条块的受力分析，得到剩余推力的计算公式如式(5)：

上式中，fi和fi-1分别表示第i条块和第i-1条块的剩余下滑力，w1i和w2i分别表示条块中浸润线以上部分的重力和浸润线以下部分的浮重，αi和αi-1表示第i条块和第i-1条块的滑动面倾角，ci和分别表示滑坡体第i条块的粘聚力和内摩擦角，fs表示滑坡的稳定系数，di＝γwaisinσi为渗透压力。其中，γw为水的容重，ai为条块中位于浸润线以下的面积，βi是水力坡降的方向和水平方向的夹角。第二个表达式中，ψi-1是剩余推力的传递系数；公式(5)的第一个表达式右边一共由三项组成，其中第一项为该条块的下滑力，第二项为该条块的抗滑力，第三项为该条块受到的剩余推力。在求第一个表达式的解时，需要进行迭代计算，首先将第一个条块的第三项设置为0，即第一个条块没有受到传递来的推力，然后依次计算各个条块的剩余下滑力。通过第二个表达式在fs的值区间内不断迭代不同的值，以最后的条块的剩余推力为0作为停止迭代的条件，停止迭代时的fs即为滑坡当前状态的稳定系数，各个条块的剩余推力表示滑坡此时的受力分布情况。

在本发明实施例中，以湖北省秭归县白水河滑坡为实验对象，对本发明的技术方案做了进一步解释说明：

以湖北省秭归县白水河滑坡为实验对象，验证本发明提出的滑坡位移预测方法，并和现有的基于bp神经网络、支持向量回归机和长短期记忆神经网络算法的三种滑坡位移预测模型进行对比。

选取2013/5～2015/5的深部位移数据来确定滑带完整性指标，并利用2013/5～2016/6的降雨和库水监测数据对滑带完整性指标进行回归分析。

利用条分法将滑坡体剖面划分为84个条块，效果如图4。

根据浸润线计算公式计算白水河滑坡不同时刻的浸润线位置，根据条块滑动面和浸润线相对位置由滑带完整性指标得计算条块的抗剪强度参数。由于位移监测点位于第39块，所以实验中主要求第39条块的动态抗剪强度参数，如图5。

利用剩余推力法对滑坡的内部受力分布进行迭代分析，结合监测点的累计位移利用对应条块的剩余推力进行滑坡位移预测。对滑带完整性指标和位移预测值进行数据同化分析，得到校正和优化后的滑带完整性指标和位移预测值，并作为预测的初始值，进而得到最后的位移预测值。

为了评价本申请提出的利用数据同化改进的滑坡位移预测方法的效果，引入两种评价标准：均方根误差rmse和平均绝对误差mae，如式(4)所示。

上式中，yi为监测点的位移实际监测值，为位移预测值。rmse可以衡量预测值与监测值的误差，mae可以衡量整体预测误差的平均水平。

表1给出了4种方法的预测误差的mae值和rmse值：

表1不同方法的预测误差对比

可以看出，利用粒子滤波数据同化方法改进的滑坡位移预测模型能够明显地提高对滑坡位移的预测精度。粒子滤波同化算法通过同步估计滑带完整性指标能够降低模型参数的不确定性，从而降低模型的结构误差，更好地提高了预测精度。利用粒子滤波数据同化方法改进的滑坡位移预测模型在预测后期的误差比较小，而且相对比较稳定。

本发明的有益效果是：本发明所提出的技术方案具备以下优点：

1)数值模拟是在滑坡的几何信息和土体参数等信息上进行的，是基于动力学方程进行迭代计算的，可以保证滑坡预测模型不会对外界影响因素的剧烈噪声做出过激响应；

2)针对滑坡位移预测模型存在结构误差和参数误差的问题，基于数据同化理论，利用粒子滤波算法从新的监测数据中提取有效信息，逐步优化模型状态和校正关键参数，可以在很大程度上提高滑坡位移预测模型的预测精度，并且能有效地让滑坡位移预测模型的预测误差收敛。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。