CN108537360A - 一种顾及多因子的自适应Kalman滤波滑坡预报方法 - Google Patents

Abstract

本发明公开了一种顾及多因子的自适应Kalman滤波滑坡预报方法，该方法首先构建Kalman滤波状态方程和观测方程，并建立基于多因子的Kalman滤波递推方程组，引入指数型衰减因子进行加权衰减记忆Kalman滤波，最终实现对滑坡预测的预报。本发明方法考虑滑坡外部诱发因子的随机干扰，减少了滑坡变形模型在线性化中产生的误差；同时引入了衰减因子，有效抑制了经典Kalman滤波易发散的问题，极大提高了滤波的准确性和滤波精度，从而进一步保证了采用顾及多因子的自适应Kalman滤波方法进行滑坡灾害预报的可靠性。

Description

一种顾及多因子的自适应Kalman滤波滑坡预报方法

技术领域

本发明涉及一种滑坡预报方法，在顾及了与滑坡密切相关的多种诱发因子(如大气降雨、温度变化等)对滑坡影响的基础上，提出了一种顾及多因子的自适应Kalman滤波滑坡预报方法，可便捷有效地实现对滑坡的精准预报。

背景技术

滑坡是一种破坏性很强的自然灾害，具有灾害影响大、发生次数频繁、分布范围广以及灾害隐蔽性大等特点。近年对我国经济可持续发展与人民生命财产安全造成了巨大的危害，同时也对交通、水利、采矿等造成了严重的影响，特别是一些大型、复杂滑坡的影响更甚。因此，对滑坡进行有效的预报成为减少滑坡灾害损失的重要技术手段，对于切实保障国民经济可持续发展、人民生命财产安全等具有十分重要的理论与现实意义。

滑坡发生的位移量与位移速率是评价滑坡稳定性最直接的判定依据，且随着多种高精度监测手段的迅速发展，利用多类滑坡高精度监测数据并顾及多种滑坡外部诱发因子，通过构建有效的方法模型对滑坡发展进行分析，已成为滑坡精确预报的有效手段与必然趋势。目前，国内外众多学者利用优化灰色模型、多元时间序列分析、粗糙集理论、Kalman滤波方法等对滑坡变形位移进行了预报，取得了较好的成果。然而，上述预报模型主要仅基于位移参数对滑坡进行预报，而滑坡的发生与发展不仅会受到滑坡自身土体参数、外部地貌特征的控制，以及一些人类活动(工程开挖)等影响，特别会受到外部环境(大气降雨、温度变化等)的影响。因此，仅利用单一位移参数的预报模型已不能满足滑坡预报高可靠性、高精度的迫切需要。此外，在已有利用Kalman滤波开展的滑坡预报方法中，存在着当描述系统动力学特性的数学模型和噪声统计模型不准确、或递推过程中舍入误差逐渐积累使滤波增益阵的计算值逐渐失去合理的加权作用时，则依据标准Kalman滤波估计状态变量时必然会导致滤波发散的严重问题。

发明内容

本发明的目的是提供一种顾及多因子的自适应Kalman滤波方法，用于滑坡预报，拟解决仅顾及单一因子且不考虑外部多影响因子预报模型进行滑坡预报精度低，以及采用经典Kalman滤波预报易发散的问题。

为了实现上述任务，本发明采用以下技术方案：

一种顾及多因子的自适应Kalman滤波滑坡预报方法，包括以下步骤：

步骤1，根据诱发滑坡的多种外部因子，建立顾及多因子的自适应Kalman滤波状态方程和观测方程，将所述的状态方程离散化、线性化，从而建立状态方程的近似线性方程；

步骤2，通过递推方程组对所述的近似线性方程求解，得到状态估计量及其误差方差矩阵；

步骤3，在所述递推方程组中引入衰减因子λ_k，得到改进后的递推方程组；所述的衰减因子λ_k的表达式为：

λ_k＝(1-b^k)/(1-b)

上式中，0<b<1，k表示观测时刻；

步骤4，给定滤波初始值和观测数据，通过所述改进后的递推方程组求得k时刻的状态估计量及其误差方差矩阵，根据k时刻到k时刻之后的某一时刻的观测数据确定出状态转移矩阵，从而求得某一时刻的滑坡变形量的预报值。

进一步地，所述的多种外部因子包括月降雨量r、温度T。

进一步地，步骤1中所述的状态方程和观测方程分别为：

x＝x(t,r,T)

上式中，t为观测时间，r为t时间的月降雨量，T为t时间的温度，x为滑坡变形量，为k时刻的状态估计量，L_k为k时刻的观测向量，B_k＝[1 0 0 0 0 0 0]，Δ_k为观测噪声。

进一步地，步骤1中所述的近似线性方程为：

上式中，分别为k时刻、k-1时刻的状态估计量，Φ_k,k-1为k-1时刻到k时刻的状态转移矩阵，Γ_k,k-1为动态噪声系数矩阵，Ω_k-1为动态噪声向量。

进一步地，步骤2中所述的递推方程组为：

上式中，和分别为k时刻的一步预测量和状态估计量，为k-1时刻的状态估计量，I为单位矩阵，J_k为k时刻的滤波增益矩阵，L_k为k时刻的观测向量，为一步预测误差协方差阵，为k时刻估计误差协方差矩阵，为k-1时刻估计误差协方差矩阵，Φ_k,k-1为k-1时刻到k时刻的状态转移矩阵，B_k为观测值系数矩阵，Γ_k,k-1为动态噪声系数矩阵，为k-1时刻动态噪声协方差矩阵，为k时刻观测噪声协方差矩阵。

进一步地，步骤3中所述的改进后的递推方程组为：

上式中各个参数的含义与递推方程组中的含义相同，只是在改进后的递推方程中对参数上标*加以区分。

本发明与现有技术相比具有以下技术特点：

1.本发明首先采用基于多因子及泰勒级数展开的Kalman滤波状态方程，考虑滑坡外部诱发因子(大气降雨、温度变化等)的随机干扰，减少了滑坡变形模型在线性化中产生的误差；同时引入了衰减因子，构建了指数加权衰减记忆自适应滤波，有效抑制了经典Kalman滤波易发散的问题，极大提高了滤波的准确性，从而进一步保证了采用顾及多因子的自适应Kalman滤波方法进行滑坡灾害预报的可靠性。

2.本发明提出的顾及多因子的自适应Kalman滤波滑坡预报模型结构简单，可操作性强，模拟效果好，可实现对滑坡的实时跟踪预报。顾及多因子的自适应Kalman滤波通过引入指数型衰减因子能够较好地抑制状态异常的影响，滤波精度较高，能更好地为滑坡防灾减灾服务。

附图说明

图1为本发明方法的流程图；

图2为实施例中MD09点垂直方向位移监测数据图；

图3是采用不同方法对滑坡监测数据滤波残差图，其中，(a)图为经典Kalman滤波残差图，(b)图为使用本发明顾及多因子的自适应Kalman滤波残差图。

具体实施方式

本发明提供了一种顾及多因子的自适应Kalman滤波滑坡预报方法，如图1所示，具体包括以下步骤：

步骤1，根据诱发滑坡的多种外部因子，建立顾及多因子的自适应Kalman滤波状态方程和观测方程，将所述的状态方程离散化、线性化，从而建立状态方程的近似线性方程；本方案中，考虑到的多种外部因子包括大气降雨、温度变化等。

步骤1.1，建立顾及多因子的自适应Kalman滤波状态方程和观测方程

顾及滑坡外部多因子影响，构建状态方程如下：

x＝x(t,r,T)

其中t为观测时间，r为t时间的月降雨量，T为t时间的温度，x为滑坡变形量。

建立滑坡变形Kalman滤波观测方程：

边坡上某点变形观测为：

其中，L_k分别为k时刻监测的状态估计量和观测向量，B_k＝[1 0 0 0 0 0 0]，Δ_k为观测噪声，上式即为滑坡变形Kalman滤波的观测方程。

步骤1.2，对所述的状态方程离散化、线性化，从而建立状态方程的近似线性方程；

由于观测时间较短且变形量级较小，可将变形量x在t_k时刻进行泰勒级数展开为t、r、T的函数，舍去三次以上的项，可得如下公式：

x_k(t_k,r_k,T_k)＝x_k-1(t_k-1,r_k-1,T_k-1)+v_k-1(t_k-t_k-1)+a_k-1(t_k-t_k-1)²+b_k-1(r_k-r_k-1)+c_k-1(r_k-r_k-1)²+d_k-1(T_k-T_k-1)+e_k-1(T_k-T_k-1)²

上式中，x_k-1为t_k-1时刻的滑坡变形量，为t_k-1时刻滑坡变形速度，为t_k-1时刻滑坡变形加速度，为月大气降雨量变化对滑坡影响，为月大气降雨量变化对滑坡的随机干扰，为温度变化对滑坡影响，为温度变化对滑坡的随机干扰。

考虑滑坡外界诱发因子(大气降雨量、温度变化)的随机干扰，且将位移加速度视为动态噪声，则动态噪声系数阵为：

(Φ(t_k,τ)为τ时刻到t_k时刻的状态转移矩阵，F(τ)＝[0 0 1 0 1 0 1]^T)；

可得以下状态方程：

令其中，Δt_k-1＝t_k-t_k-1，Δr_k-1＝r_k-r_k-1，ΔT_k-1＝T_k-T_k-1，Ω_k-1为动态噪声向量(该值服从均值为零、方差为1的正态分布)，Φ_k,k-1为状态转移矩阵：

则得到Kalman滤波状态方程的近似线性方程为：

步骤2，通过递推方程组对所述的近似线性方程求解，得到状态估计量及其误差方差矩阵；

步骤1中建立起了顾及多因子的自适应Kalman滤波状态方程和观测方程，下一步则需要求解出状态估计量及其误差方差矩阵本方案中通过建立递推方程组对这些参数进行求解。

根据步骤1得到顾及多因子的Kalman滤波的近似线性方程和观测方程分别为：

则状态估计量及其误差方差矩阵通过以下递推方程组进行求解：

状态一步预测方程为：

一步预测误差方差矩阵为：

滤波增益矩阵为：

状态估计量为：

误差方差矩阵为：

上式中，和分别为k时刻的一步预测量和状态估计量，为k-1时刻的状态估计量，I为对角元素为1其余元素为零的单位矩阵，J_k为k时刻的滤波增益矩阵，L_k为k时刻的观测向量，为一步预测误差协方差阵，为k时刻估计误差协方差矩阵，为k-1时刻估计误差协方差矩阵，Φ_k,k-1为k-1时刻到k时刻的状态转移矩阵，B_k为观测值系数矩阵，Γ_k,k-1为动态噪声系数矩阵，为k-1时刻动态噪声协方差矩阵，为k时刻观测噪声协方差矩阵，T表示矩阵的转置。

步骤3，在所述递推方程组中引入衰减因子λ_k，得到改进后的递推方程组；

经典Kalman滤波随滤波次数及观测数据的增加，会导致滤波增益阵J_k趋向无穷大(发散)，使得状态估计量与实测值不符。针对这个问题，本方案中进一步引入衰减因子λ_k，建立衰减自适应Kalman滤波，从而为滑坡灾害预报提供更可靠的形变量预报值。

为抑制经典Kalman滤波发散的问题，引入衰减因子λ_k，可得改进后的衰减记忆自适应Kalman滤波递推方程组：

上式中各个参数的含义与递推方程组中的含义相同，只是在改进后的递推方程中对参数上标*加以区分，即：和分别为的一步预测量和状态估计量，I为单位矩阵，为k-1时刻的状态估计量，为k时刻滤波增益矩阵，L_k为k时刻的观测向量，为一步预测误差协方差阵，为k时刻估计误差协方差矩阵，为k-1时刻估计误差协方差矩阵，Φ_k,k-1为k-1时刻到k时刻的状态转移矩阵，B_k为观测值系数矩阵，Γ_k,k-1为动态噪声系数矩阵，为k-1时刻动态噪声协方差矩阵，为k时刻观测噪声协方差矩阵，λ_k为k时刻衰减因子。

取衰减因子为指数形式：

λ_k＝(1-b^k)/(1-b)

上式中，0<b<1，b为自然实数，b值的选取可根据实际问题中模型预报与实测数据达到的最佳拟合度而确定。

将改进后的递推方程组中的合并后得到：

相比较与标准Kalman滤波，改进后的衰减记忆自适应Kalman滤波由于引入了衰减因子λ_k且λ_k>1，衰减记忆滤波使滤波增益阵变大，因此新近观测值权值增大，而的系数减小，即削弱了陈旧观测值对滤波的影响。

由衰减因子λ_k形式可知其满足：1)λ_k>1，符合衰减因子大于1的基本要求；2)λ_k为指数型增函数，使得新近观测值值权重逐渐增大，体现了自适应性的优点。

步骤4，给定滤波初始值和观测数据，通过所述改进后的递推方程组求得k时刻的状态估计量及其误差方差矩阵，根据k时刻到k时刻之后的某一时刻的观测数据确定出状态转移矩阵，从而求得某一时刻的滑坡变形量的预报值。

步骤4.1，根据实际观测数据及其精度，可以确定出变形初始值X₀及其协方差阵观测噪声方差阵动态噪声向量阵的协方差阵一般取单位阵I。

步骤4.2，根据步骤4.1中确定的初值以及观测数据，包括滑坡变形量、大气降雨量、温度及观测时间；将这些量值代入到步骤3建立的改进后的递推方程组中，即可计算得到k时刻的状态估计量及其协方差矩阵然后根据k时刻到k时刻之后的某一时刻的观测数据确定出状态转移矩阵，从而求得某一时刻的滑坡变形量的预报值。

例如，需要确定k+1、k+2、k+3、k+4时刻的预报值时，通过k时刻到k+1、k+2、k+3、k+4时刻的降雨量、温度、时刻变化，确定出k+1、k+2、k+3、k+4时刻的状态转移矩阵Φ_k+1,k、Φ_k+2,k、Φ_k+3,k、Φ_k+4,k，最终改进后的递推方程组即可计算出k+1、k+2、k+3、k+4时刻的滑坡变形量的预报值，如下：

实施例：

以陕西某地滑坡为例，选用研究区内变形明显的MD09点2015.07-2016.10垂直方向位移监测资料进行计算(图2)，2016.11-2017.02监测资料用于滑坡变形量预报。

2.建立该点的Kalman滤波状态方程和观测方程，并将状态方程线性化，取至二次项，考虑滑坡外部诱发因子的随机干扰，且将位移加速度视为动态噪声，通过积分获得动态噪声系数阵，并引入指数衰减因子，得到改进后的顾及多因子的自适应Kalman滤波状态方程和观测方程。

3.采用MATLAB编程分别进行标准Kalman和顾及多因子的自适应Kalman滤波

依据构建的Kalman滤波观测方程及状态方程，即可进行经典Kalman滤波和顾及多因子的自适应Kalman滤波，获得两种滤波方法的残差图，如图3所示。图3中(a)图所示，采用经典Kalman滤波处理，残差值最大为0.0518mm，而(b)图采用顾及多因子的自适应Kalman滤波，残差值最大0.00083mm，精度有显著提高，且随着滤波次数的增加，滤波效果越来越好，有效抑制了经典Kalman滤波易发散的现象，为滑坡预报提供了可靠的数据。

4.依据已有研究区监测数据结合所建顾及多因子的自适应Kalman滤波模型进行滑坡预报，结果如表1所示。由表1可知，采用本发明的顾及多因子的自适应Kalman滤波能够有效抑制滤波发散，取得较好滑坡变形量预报效果。

表1滑坡预报结果

Claims (6)

1.一种顾及多因子的自适应Kalman滤波滑坡预报方法，其特征在于，包括以下步骤：

步骤1，根据诱发滑坡的多种外部因子，建立顾及多因子的自适应Kalman滤波状态方程和观测方程，将所述的状态方程离散化、线性化，从而建立状态方程的近似线性方程；

步骤2，通过递推方程组对所述的近似线性方程求解，得到状态估计量及其误差方差矩阵；

步骤3，在所述递推方程组中引入衰减因子λ_k，得到改进后的递推方程组；所述的衰减因子λ_k的表达式为：

λ_k＝(1-b^k)/(1-b)

上式中，0<b<1，k表示观测时刻；

步骤4，给定滤波初始值和观测数据，通过所述改进后的递推方程组求得k时刻的状态估计量及其误差方差矩阵，根据k时刻到k时刻之后的某一时刻的观测数据确定出状态转移矩阵，从而求得某一时刻的滑坡变形量的预报值。

2.如权利要求1所述的顾及多因子的自适应Kalman滤波滑坡预报方法，其特征在于，所述的多种外部因子包括月降雨量r、温度T。

3.如权利要求1所述的顾及多因子的自适应Kalman滤波滑坡预报方法，其特征在于，步骤1中所述的状态方程和观测方程分别为：

x＝x(t,r,T)

上式中，t为观测时间，r为t时间的月降雨量，T为t时间的温度，x为滑坡变形量为k时刻的状态估计量，L_k为k时刻的观测向量，B_k＝[1 0 0 0 0 0 0]，Δ_k为观测噪声。

4.如权利要求1所述的顾及多因子的自适应Kalman滤波滑坡预报方法，其特征在于，步骤1中所述的近似线性方程为：

上式中，分别为k时刻、k-1时刻的状态估计量，Φ_k,k-1为k时刻的状态转移矩阵，Γ_k,k-1为动态噪声系数矩阵，Ω_k-1为动态噪声向量。

5.如权利要求1所述的顾及多因子的自适应Kalman滤波滑坡预报方法，其特征在于，步骤2中所述的递推方程组为：

上式中，和分别为k时刻的一步预测量和状态估计量，为k-1时刻的状态估计量，I为单位矩阵，J_k为k时刻的滤波增益矩阵，L_k为k时刻的观测向量，为一步预测误差协方差阵，为k时刻估计误差协方差矩阵，为k-1时刻估计误差协方差矩阵，Φ_k,k-1为k-1时刻到k时刻的状态转移矩阵，B_k为观测值系数矩阵，Γ_k,k-1为动态噪声系数矩阵，为k-1时刻动态噪声协方差矩阵，为k时刻观测噪声协方差矩阵。

6.如权利要求1所述的顾及多因子的自适应Kalman滤波滑坡预报方法，其特征在于，步骤3中所述的改进后的递推方程组为：

上式中各个参数的含义与递推方程组中的含义相同，只是在改进后的递推方程中对参数上标*加以区分。