CN108332649B

CN108332649B - 一种滑坡形变综合预警方法及***

Info

Publication number: CN108332649B
Application number: CN201810123208.0A
Authority: CN
Inventors: 王守华; 陆明炽; 孙希延; 纪元法; 李云柯; 吴孙勇; 肖建明; 邓洪高; 符强; 严素清
Original assignee: Guilin University of Electronic Technology
Current assignee: Guilin University of Electronic Technology
Priority date: 2018-02-07
Filing date: 2018-02-07
Publication date: 2020-04-24
Anticipated expiration: 2038-02-07
Also published as: CN108332649A

Abstract

本发明公开一种滑坡预警方法及***，包括滑坡形变监测网和综合预警单元；所述滑坡形变监测网由卫星定位基准站、3个以上的卫星定位监测站；3个以上的卫星定位监测站分布在监测目标的地表的不同监测点处；每个卫星定位监测站分别与卫星定位基准站连接，实现RTK定位解算；所有卫星定位监测站的定位结果均传输到综合预警单元。本发明通过监测目标整体姿态和位移量的精确解算，实现形变位移的预测，并通过融合多源监测信息，实现综合预警。本发明具有预警精确高，算法简单、易实现，有较强的实用性的特点，并能满足山体滑坡形变预警的实际要求。

Description

一种滑坡形变综合预警方法及***

技术领域

本发明涉及地质预警技术领域，具体涉及一种滑坡形变综合预警方法及***。

背景技术

由于水电等重点工程开工建设或环境气候变化等因素的影响，我国山体滑坡、地面变形、地震、土地退化等地质灾害也日益严重。滑坡、崩塌、地面沉降等己成为我国西南山区等地区主要的地质灾害之一，而积极有效地防治滑坡等地质灾害的重要措施之一是对滑坡等形变隐患点进行长期的观测。滑坡隐患点一旦发生形变，若不能及时监测和有效预警，则可能给当地居民带来巨大经济损失，甚至危害居民生命安全。

发明内容

本发明提供一种滑坡形变综合预警方法及***，其具有实现简单和预警精确高的特点。

为解决上述问题，本发明是通过以下技术方案实现的：

一种滑坡形变综合预警方法，具体包括步骤如下：

步骤1、滑坡形变监测网内各个卫星定位监测站解算出其所处监测目标的监测点处三维位置坐标，并根据各卫星定位监测站初始位置求解出位移量；

步骤2、利用卡尔曼滤波器对各个卫星定位监测站所得到的位移量进行平滑滤波后得到形变位移量的实测值；并根据采样时间间隔及前后时刻形变位移量的实测值的位移差，求解出监测目的形变速度和形变加速度的实测值；

步骤3、提取步骤2所得各个监测点的形变位移量的实测值，并利用最小二乘法解算出监测目标的姿态，即监测目标的航向角、俯仰角和横滚角的实测值；

步骤4、提取步骤2所得各个监测点的形变位移量、形变速度和形变加速度的实测值，并利用扩展卡尔曼滤波器解算出监测目标的综合位移的实测值；

步骤5、设定预测滑动窗的起始历元及长度，对监测目标的综合位移进行等时间间隔采样，利用迭代灰色模型实现监测目标的形变位移量的预测值，并根据采样时间间隔及前后时刻形变位移量的预测值的位移差，计算形变速度和形变加速度的预测值；

步骤6、利用步骤3和4所解算出的形变位移量、形变速度和形变加速度的实测值，根据预定的实测预警规则，对监测目标的当前状态进行实测滑坡等级预警；同时，利用步骤5所预测出的形变位移量、形变速度和形变加速度的预测值，根据预定的预测预警规则，对监测目标的后续状态进行预测滑坡等级预警。

上述步骤5的具体过程为：

步骤5.1、设定预测滑动窗起始历元，对监测目标的综合位移的实测值进行采样，得到当前预测滑动窗的观测数列；

步骤5.2、将观测数列进行一次累加，得到数据生成序列，并根据数据生成序列计算其紧邻均值生成序列；

步骤5.3、利用紧邻均值生成序列构建雅克比矩阵，并根据雅克比矩阵构建最小二乘估计参数系数，即发展系数和灰色作用量；

步骤5.4、利用发展系数和灰色作用量对观测数列进行一次累减，得到当前预测滑动窗的预测数列；

步骤5.5、利用预测滑动窗将起始的观测数据不断往前推移，并通过重复步骤5.2-5.4进行迭代计算，得到各个预测滑动窗的预测数列；

步骤5.6、将各个预测滑动窗的预测数列进行加权累加后，得到监测目标的形变位移量的预测值；

步骤5.7、根据采样时间间隔及前后时刻形变位移量的预测值的位移差，计算形变速度和形变加速度的预测值。

上述步骤3的具体过程为：

步骤3.1、提取步骤2所得各个监测点的形变位移量的实测值；

步骤3.2、利用滑坡形变监测网内任意3个以上个不共线的卫星定位监测站的三维位置坐标构建基线向量，并任意选择2条基线向量解算俯仰角、航向角和横滚角的初始值；

步骤3.3、根据构建的姿态系数矩阵计算基线对应权重，并利用最小二乘法估计出姿态改正数；

步骤3.4、根据姿态改正数及俯仰角、航向角和横滚角的初始值解算监测目标的姿态，即监测目标的航向角、俯仰角和横滚角的实测值。

上述步骤4的具体过程为：

步骤4.1、提取步骤2所得各个监测点的形变位移量、形变速度和形变加速度的实测值；

步骤4.1、利用形变位移量和形变速度构建扩展卡尔曼滤波器的状态量矩阵、全局观测矩阵和测量关系矩阵，利用采样时间间隔构建扩展卡尔曼滤波器的状态转移矩阵，根据形变加速度解算扩展卡尔曼滤波器的过程噪声协方差矩阵和过程噪声协方差矩阵；

步骤4.1、利用扩展卡尔曼滤波器对时间进行更新，预测当前***状态，并通过测量更新，校正预测状态先验估计值，得到监测目标整***移的全局最优估计，即监测目标的综合位移的实测值。

作为改进，所述滑坡形变综合预警方法还进一步包括：步骤7、滑坡形变监测网内的雨量计采集降雨量采集降雨量和裂缝传感器采集裂缝位移量，并将降雨量和裂缝位移量与步骤4和步骤5所得到的形变位移量、形变速度和形变加速度的实测值和预测值进行加权运算，并根据预定的综合预警规则，对监测目标的进行综合滑坡等级预警。

实现上述方法的一种滑坡形变综合预警***，包括滑坡形变监测网和综合预警单元；所述滑坡形变监测网由卫星定位基准站、3个以上的卫星定位监测站；3个以上的卫星定位监测站分布在监测目标的地表的不同监测点处；每个卫星定位监测站分别与卫星定位基准站连接，实现RTK定位解算；所有卫星定位监测站的定位结果均传输到综合预警单元。

上述方案中，卫星定位监测站为GPS卫星定位接收机、BDS卫星定位接收机和/或GNSS卫星定位接收机。

作为改进，所述滑坡形变综合预警***的滑坡形变监测网还进一步包括雨量计和2个以上的裂缝传感器；雨量计分布在监测目标的地表上；2个以上的裂缝传感器分布在监测目标的地表的不同监测点处；雨量计和所有裂缝传感器的监测数据均传输到综合预警单元。

上述方案中，卫星定位监测站的数量与裂缝传感器的数量相同；此时，监测目标的地表的每个监测点处均设有1个卫星定位监测站和1个裂缝传感器。

与现有技术相比，本发明通过监测目标整体姿态和位移量的精确解算，实现形变位移的预测，并通过融合多源监测信息，实现综合预警。本发明具有预警精确高，算法简单、易实现，有较强的实用性的特点，并能满足山体滑坡形变预警的实际要求。

附图说明

图1为滑坡形变综合预警***的原理框图。

图2为监测目标姿态算法流程图。

图3为多卫星定位监测站扩展卡尔曼滤波数据融合算法流程图。

图4为迭代灰色模型形变预测算法流程图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实例，并参照附图，对本发明进一步详细说明。

一种滑坡形变综合预警***，如图1所示，主要由滑坡形变监测网和综合预警单元组成。

上述滑坡形变监测网包括卫星定位滑坡形变监测网和滑坡监测辅助网两部分。卫星定位滑坡形变监测网由1个卫星定位基准站和N个卫星定位监测站组成，其作用是完成监测目标地表位移的监测。卫星定位基准站设置在地基稳固、无信号遮蔽、无大功率无线电发射源的位置，N个卫星定位监测站分别设置在监测目标潜在形变位移方向的不同监测点上。卫星定位基准站和各个卫星定位监测站之间的距离一般不能超过5000米。每个卫星定位监测站均与卫星定位基准站连接，并分别实现RTK精密定位解算，获得每个卫星定位监测站的三维位置坐标、形变位移量、速度和加速度等信息。在本发明优选实施例中，卫星定位滑坡形变监测网为BDS、GPS或GNSS滑坡形变监测网。滑坡形变监测辅助网由一个雨量计和若干个裂缝传感器组成，其作用是分别完成当地降雨量信息的收集及裂缝大小的测量。雨量计设置在监测目标内，若干个裂缝传感器设置监测目标的不同监测点上。裂缝传感器的数量和安放位置不一定要与卫星定位监测站的数量和安放位置相一致。但在本发明优选实施例中，为了简化***的安装和维护，裂缝传感器的数量与卫星定位监测站的数量相一致，且监测目标的地表的每个监测点处均设有1个卫星定位监测站和1个裂缝传感器。

上述综合预警单元实现滑坡形变的综合预警，其中包括：利用卡尔曼滤波器对各个监测点的RTK定位信息数据平滑，剔除野值，提高监测精度，提取有效位移量等信息；利用姿态解算算法对监测目标的实时姿态解算，并根据形变预测信息实现姿态预测；利用扩展卡尔曼滤波器实现监测网内多站点的位移形变、速度和加速度数据融合，实现综合位移的最优估计；利用MGM(1，1)灰色模型算法实现形变位移的预测；滑坡形变综合预警模块通过对姿态角信息、实时形变报警等级信息、预测形变报警等级信息、降雨量信息、裂缝大小信息等的融合，做出滑坡形变等级的判断，从而实现滑坡形变的综合预警。

(1)滑坡体多卫星定位监测站数据融合算法

在实际的滑坡监测中，面对同一监测目标，为了预警其整体的形变趋势，滑坡监测网往往采用一个滑坡卫星定位基准站配套多个卫星定位监测站的布网方式。形变预警分析时，为了避免人为误差和其他不确定因素，则需要融合多个卫星定位监测站的位移监测数据进行综合预警，而不是仅仅只根据某单一卫星定位监测站而进行预报。因此，设计扩展卡尔曼滤波器对多卫星定位监测站定位数据进行融合，实现监测目标整***移的全局最优估计。具体设计如下：

监测目标开始产生形变的阶段，可建立状态方程和观测方程如下：

X_E(k+1)＝F_EX_E(k)+G_EW_E(k) (1)

Z_E(k+1)＝H_EX_E(k)+V_E(k) (2)

式中：

k表示k时刻，x(k)为滑坡体累计位移量，v(k)为滑坡体速度；

T为RTK定位***采样时间间隔；

W_E(k)＝a(k)为滑坡体相邻两个时刻加速度；Z_E(k)为监测目标形变位移量；H_E＝[10]；V_E(k)为高斯白噪声序列，其方差为σ²、均值为零，且与W_E(k)不相关。监测目标运动模型可描述为

假设需要融合的卫星定位监测站有N个，则全局观测方程(2)中Z_E(k)、H_E(k)、V_E(k)可描述为：

Z_E(k)＝[Z₁ ^T(k),Z₂ ^T(k),…Z_N ^T(k)]^T (4)

H_E(k)＝[H₁ ^T(k),H₂ ^T(k),…H_N ^T(k)]^T (5)

V_E(k)＝[V₁ ^T(k),V₂ ^T(k),…V_N ^T(k)]^T (6)

扩展Kalman滤波算法应用于该模型的具体算法如下：

P_E(k+1|k)＝F_E(k)P_E(k|k)F_E ^T(k)+G_E(k)Q_E(k)G_E ^T(k) (9)

K_Ei(k+1)＝P_E(k+1)H_Ei ^T(k+1)R_Ei ^-1(k+1) (10)

过程噪声为W_E，可设W_E(k)＝a(k)，解算过程噪声协方差矩阵Q_E，Q_E＝Cov(W_E)＝E(W_EW^T _E)，Q_E为2N×2N的对称矩阵。测量噪声向量为V_E，解算过程噪声协方差矩阵R_E，R_E＝Cov(V_E)＝E(V_EV^T _E)，R_E为2N×2N的对称矩阵。

由此，可实现多监测点位移信息数据融合得到监测目标的综合位移信息。

(2)监测目标姿态的求解

在变形监测时，监测目标的各项姿态角指标可以有效反应形变体的变化情况，故可利用多天线(多监测点)前后历元的位置坐标变化信息求解监测目标的俯仰角、横滚角、航向角等姿态值。求解姿态的涉及到两个坐标系：地理坐标系、载体坐标系。地理坐标系(L系)以主天线的相位中心为原点，以当地北子午线为Y轴，X轴指向东且垂直于Y轴，Z轴与X、Y轴构成右手坐标系。载体坐标系(b系)以主天线的相位中心为原点；Y轴指向M号与M+1号天线的连线方向，X轴与Y轴垂直且指向监测目标的右侧，Z轴与X、Y轴构成右手坐标系。

由于监测网中选定某个监测站作为主监测站，主监测站可以和其它监测站(不共线)分别在地理坐标系和载体坐标系构成的N-1条基线，根据坐标系转换的基本原理可得：

式中：α_i为第i条基线在载体坐标系中的基线矢量、S_i为第i条基线在地理坐标系中的基线矢量，C^b _L为坐标旋转矩阵。

其中

ψ为航向角，θ为俯仰角，γ为横滚角。

在地理坐标系中任意选择两条基线向量：S₁₂＝[x₁₂ y₁₂ z₁₂]^T、S₁₃＝[x₁₃ y₁₃ z₁₃]^T。由姿态角的定义可以求初始姿态角，俯仰角θ°和航向角ψ°初始值分别如下：

ψ°＝-arctan(x₁₂/y₁₂)(13)

横滚角初始值γ°可以通过以下变换得到：

γ°＝-arctan(z₁′₃/x₁′₃)(15)将式(11)写成误差形式，可得：

式中，δψ、δθ、δγ为待估参数，即姿态角改正数。观测值权阵

分别为S_i和α_i的协方差阵。

根据最小二乘法的原理可估算姿态改正数：

所以最小二乘估算态姿角后，监测目标的航向角、俯仰角和横滚角分别为

(3)滑坡形变预测灰色MGM(1，1)模型

灰色预测模型的数据模型最佳的长度是40组数据，组数过多则会造成误差积累，预测精度会受到影响。而在BDS/GPS滑坡监测的实际应用中，一般是一个历元解算一次，即每个历元都会得到一个定位结果，若对每个定位结果都进行预测分析，则会造成计算量过大，数据冗余量增加。现对各个历元的定位数据进行采样，采样间隔设为T_c，T_c可以根据实际需求而设定。现T_c设为4个小时。采样的方法是取一段时间内位移量的平均值，以避免野值等因素的影响。采样后的观测值数列X⁽⁰⁾如下：

X⁽⁰⁾＝{x⁽⁰⁾(1),x⁽⁰⁾(2),…,x⁽⁰⁾(n)} (21)

式中，n为观测数列的长度，即数据模型的长度，现将n设为40；x⁽⁰⁾(k)为形变位移量，在实际的滑坡预测模型中，因为形变位移一般是单向的，所以x⁽⁰⁾(k)为一个非负的数值。根据灰色GM(1，1)模型的基本理论，现将观测数列X⁽⁰⁾一次累加，得到数据生成序列为X⁽¹⁾，具体的公式如下：

X⁽¹⁾＝{x⁽¹⁾(1),x⁽¹⁾(2)，…,x⁽¹⁾(n)} (22)

其中

由数据生成序列X⁽¹⁾可计算其紧邻均值生成序列Z⁽¹⁾，具体公式如下：

Z⁽¹⁾＝{z⁽¹⁾(1),z⁽¹⁾(2),z⁽¹⁾(3),…,z⁽¹⁾(n)} (23)

其中，z⁽¹⁾(k)为白化背景值，取值为前后生成两个生成数据X⁽¹⁾的平均值；z⁽¹⁾(k)＝0.5x⁽¹⁾(k)+0.5x⁽¹⁾(k-1),(k＝2,3,…,n)；

对X⁽¹⁾构建模仿白微分方程：

符号替换并整理得式(24)的微分方程：

x⁽⁰⁾(k)+a_iz⁽¹⁾(k)＝b_i (25)

式(25)即为GM(1，1)模型。式中，a_i、b_i都是方程需要待定参数，a_i能反映形变位移量x⁽⁰⁾的发展趋势，称为发展系数；b_i能反映灰色信息覆盖的***作用大小，称为灰色作用量。

构建式(25)的白色微分方程和灰色微分方程的时间响应函数，具体公式如下：

为还原原始数列可做一次累减，累减后公式变量替换和整理，则有

为求解参数a_i和b_i，可利用最小二乘估计对式(25)的进行解算，具体公式为：

其中，

由此可得发展系数a_i和b_i灰色作用量，将其代入公式(27)，得到各个时刻的预测值。

现以滑动窗口的形式进行预测，即灰色预测模型的数据模型长度固定，但起始的观测数据不断往前推移，从而实现预测值的迭代计算。每一个起始观测数据，可以通过式(22)至式(28)计算出一组参数a_i和b_i，其对应的预测为

具体公式如下。

对于同一时刻的滑坡形变的预测，可以通过迭代预测后有多个预测值。这些预测值

在一定范围内波动，可以通过加权的方式进行预测值拟合。若设当前预测值的个数为M_now，预测模型滑动窗结束时的预测点个数为M_end，则灰色模型预测值个数为M_all：

M_all＝M_now-M_end (30)当M_all大于某个特定数值E_num会后，

精度可能会不断降低，所以M_all在一定范围内

才有效。

式中，w_i为预测值

对应的权重。

具体来说，上述装置所实现的一种滑坡形变综合预警方法，其具体包括步骤如下：

步骤1，采集滑坡形变监测网内的各个卫星定位监测站采用RTK定位算法解算得到的当前三维位置[x_L(k),y_L(k),z_L(k)]_N，并将其保存在特定的数据文件中，以供后面做综合预警分析使用。读取数据文件卫星定位监测站位置信息，先根据各卫星定位监测站初始位置[x(0),y(0),z(0)]_N，求解各个时刻对应的形变位移量x_bin(k)_N；再根据采样时间间隔及前后时刻形变位移量的位移差求解出滑坡体的形变速度和形变加速度；上述各数据信息均为mm级信息。

步骤2，根据采集的多个卫星定位监测站形变位移量x_bin(k)_N，利用卡尔曼滤波器对x_bin(k)_N进行平滑滤波，剔除或修复异常监测值，提高监测数据精度。

利用卡尔曼滤波器对x_bin(k)_N进行平滑。针对其中一个x_bin(k)有：状态量x_p(k)与观测向量x_bin(k)的关系(测量方程)为：x_pin(k)＝C_px_bin(k)+v_bin(k)；状态方程为：x_p(k)＝A_px_p(k-1)+B_pu(k-1)，其中，测量方程的关系矩阵为：C_p＝[10]；状态方程的状态转移矩阵和输入关系矩阵分别为：

利用KF时间更新和测量更新的基本原理实现x_bin(k)的平滑，其中，通过调节过程噪声Q_p和测量噪声R_p实现最优的滤波结果，从而得到形变位移量x_pin(k)，消除野值。

同时，根据前后位移差及采样时间间隔求解滑坡体的形变速度v_pin(k)_N＝[v_pin(k),v_pin(k),v_pin(k)]_N、形变加速度a_pin(k)_N＝[a_pin(k),a_pin(k),a_pin(k)]_N等值。

步骤3，根据监测网内任意N′(N′≥3)个不共线的卫星定位监测站的形变监测位置信息构建基线向量S_i，利用最小二乘姿态解算算法计算监测目标航向角ψ，俯仰角θ，横滚角γ。

提取各卫星定位监测站平滑后的三维位置坐标x_pin(k)＝(x_k,y_k,z_k)_pin等信息，根据需求实现采样，以减少数据的计算量。对各监测采样后的数据进行时间同步匹配，并判断是否实现时间同步，若没有，则继续匹配，若实现同步则进入下一步。选定某一监测站作为主监测站，并与其它不共线的监测站位置构建基线向量S_i，利用下式实现坐标系的转换。

其中，

ψ为航向角，θ为俯仰角，γ为横滚角。

在地理坐标系中任意选择两条基线向量S₁₂、S₁₃，解算俯仰角、航向角和横滚角的初始值θ°、ψ°、γ°。构建姿态系数矩阵G_e、L_e等矩阵，并计算基线对应权重P_e。

其中，

分别为S_i和α_i的协方差阵。

根据最小二乘法的原理可结算姿态改正数：

根据姿态改正数及俯仰角、航向角和横滚角的初始值解算监测目标的姿态。

在本发明中，利用最小二乘姿态解算算法计算监测目标的姿态角，具体算法流程图如图2所示。具体的实现步骤如下：

步骤3.1：提取各卫星定位监测站平滑后的三维位置坐标x_pin(k)＝(x_k,y_k,z_k)_pin等信息，根据需求实现采样，以减少数据的计算量。

步骤3.2：对上一步各监测采样后的数据进行时间同步匹配，并判断是否实现时间同步，若没有，则继续匹配，若实现同步则进入下一步。

步骤3.3：选定某一监测站作为主监测站，并与其它不共线的监测站位置构建基线向量S_i。

步骤3.4：上一步的基线向量是在地理坐标系中的，利用下式实现坐标系的转换。

其中，

ψ为航向角，θ为俯仰角，γ为横滚角。

步骤3.5：在地理坐标系中任意选择两条基线向量S₁₂、S₁₃，解算俯仰角、航向角和横滚角的初始值θ°、ψ°、γ°。

假设基线向量为S₁₂＝[x₁₂ y₁₂ z₁₂]^T、S₁₃＝[x₁₃ y₁₃ z₁₃]^T。由姿态角的定义可以求初始姿态角，俯仰角θ°和航向角ψ°初始值分别如下：

ψ°＝-arctan(x₁₂/y₁₂)

横滚角初始值γ°可以通过以下变换得到：

γ°＝-arctan(z₁′₃/x₁′₃)

步骤3.6：构建姿态系数矩阵G_e、L_e矩阵，并计算基线对应权重P_e。

分别为S_i和α_i的协方差阵。

步骤3.7：根据最小二乘法的原理可解算姿态改正数：

步骤3.8：根据姿态改正数及俯仰角、航向角和横滚角的初始值解算监测目标的姿态。

步骤4，提取各个卫星定位监测站平滑后的位置x_pin(k)、速度v_pin(k)、加速度a_pin(k)，进行卫星定位监测站时间同步对准，构建扩展卡尔曼滤波器状态转移矩阵G_E、测量关系矩阵H_E等，利用EKF融合监测网内多个卫星定位监测站数据，实现监测目标综合位移信息最优估计。

提取各卫星定位监测站平滑后的位置x_pin(k)、速度v_pin(k)、加速度a_pin(k)，进行卫星定位监测站时间同步对准，并判断各卫星定位监测站时间是否对准。构建状态量矩阵X_E(k)：

k表示k时刻，x_pin(k)为滑坡体累计位移量，v_pin(k)为滑坡体速度。X_E(k)需根据实际融合卫星定位监测站的个数N进行扩维；

构建全局观测矩阵Z_E(k)及测量关系矩阵H_E(k)：其中，Z_E(k)为监测目标形变位移量，H_E(k)＝[10]，需根据实际融合卫星定位监测站的个数N进行扩维；

构建状态转移矩阵F_E(k)及矩阵G_E(k)：

T为RTK定位***采样时间间隔；

需根据实际融合卫星定位监测站的个数N进行扩维；

过程噪声为W_E，可设W_E(k)＝a(k)_pin，解算过程噪声协方差矩阵Q_E，Q_E＝Cov(W_E)＝E(W_EW^T _E)，Q_E为2N×2N的对称矩阵。测量噪声向量为V_E，解算过程噪声协方差矩阵R_E，R_E＝Cov(V_E)＝E(V_EV^T _E)，R_E为2N×2N的对称矩阵。

时间更新，预测当前***状态；测量更新，校正预测状态先验估计值，得到监测目标整***移的全局最优估计。

在本发明中，利用滑坡体多卫星定位监测站数据融合算法实现监测目标整***移的全局最优估计，具体算法流程图如图3所示。具体的实现步骤如下：

步骤4.1：提取各卫星定位监测站平滑后的位置x_pin(k)、速度v_pin(k)、加速度a_pin(k)，进行卫星定位监测站时间同步对准，并判断各卫星定位监测站时间是否对准。

步骤4.2：构建状态量矩阵X_E(k)：

步骤4.3：构建全局观测矩阵Z_E(k)及测量关系矩阵H_E(k)：其中，Z_E(k)为监测目标形变位移量，H_E(k)＝[10]，需根据实际融合卫星定位监测站的个数N进行扩维；

步骤4.4：构建状态转移矩阵F_E(k)及矩阵G_E(k)：

需根据实际融合卫星定位监测站的个数N进行扩维；其中，T为RTK定位***采样时间间隔；

步骤4.5：构建过程噪声协方差矩阵Q_E、测量协方差矩阵R_E：过程噪声为W_E，可设W_E(k)＝a(k)_pin，解算过程噪声协方差矩阵Q_E，Q_E＝Cov(W_E)＝E(W_EW^T _E)，Q_E为2N×2N的对称矩阵。测量噪声向量为V_E，解算过程噪声协方差矩阵R_E，R_E＝Cov(V_E)＝E(V_EV^T _E)，R_E为2N×2N的对称矩阵。

步骤4.6：时间更新，预测当前***状态；测量更新，校正预测状态先验估计值，得到监测目标整***移的全局最优估计x_all(k)。

步骤5，对监测目标综合位移信息x_all(k)进行等时间间隔采样，采样时间间隔为T，设定预测滑动窗的起始历元及长度n，利用迭代灰色模型相应原理实现监测目标滑坡形变

预测。

设定预测滑动窗起始历元，对监测目标的综合位移信息x_all(k)采样，采样间隔4小时，共采样40组，采样的方法是取一段时间内位移量的平均值。采样后的数据为：

X⁽⁰⁾＝{x⁽⁰⁾(1),x⁽⁰⁾(2),…,x⁽⁰⁾(n)}

将观测数列X⁽⁰⁾一次累加，得到数据生成序列为X⁽¹⁾：

X⁽¹⁾＝{x⁽¹⁾(1),x⁽¹⁾(2)，…,x⁽¹⁾(n)}

其中

由数据生成序列X⁽¹⁾可计算其紧邻均值生成序列Z⁽¹⁾：

Z⁽¹⁾＝{z⁽¹⁾(1),z⁽¹⁾(2),z⁽¹⁾(3),…,z⁽¹⁾(n)}

其中，z⁽¹⁾(k)为白化背景值，取值为前后生成两个生成数据X⁽¹⁾的平均值。构建雅克比矩阵B：

利用最小二乘估计参数系数：

一次累减还原原始数列：

起始的观测数据不断往前推移，实现预测值的迭代计算，得到N_win组参数a_i和b_i，其对应的预测值为

偏离平均值的比例计算各预测值

的权重w_i，且

其中，E_num为有效预测值的个数。

据权重计算预测值：

在本发明中，利用滑坡形变预测灰色MGM(1，1)模型实现监测目标形变位移量的预测，具体算法流程图如图4所示。具体的实现步骤如下：

步骤5.1：设定预测滑动窗起始历元，对监测目标的综合位移信息x_all(k)采样，采样间隔4小时，共采样40组，采样的方法是取一段时间内位移量的平均值。采样后的数据为：

X⁽⁰⁾＝{x⁽⁰⁾(1),x⁽⁰⁾(2),…,x⁽⁰⁾(n)}

步骤5.2：将观测数列X⁽⁰⁾一次累加，得到数据生成序列为X⁽¹⁾：

X⁽¹⁾＝{x⁽¹⁾(1),x⁽¹⁾(2)，…,x⁽¹⁾(n)}

其中

步骤5.3：由数据生成序列X⁽¹⁾可计算其紧邻均值生成序列Z⁽¹⁾：

Z⁽¹⁾＝{z⁽¹⁾(1),z⁽¹⁾(2),z⁽¹⁾(3),…,z⁽¹⁾(n)}

其中，z⁽¹⁾(k)为白化背景值，取值为前后生成两个生成数据X⁽¹⁾的平均值。

步骤5.4：构建雅克比矩阵B：

步骤5.5：利用最小二乘估计参数系数：

其中：

a_i为发展系数，b_i为灰色作用量。

步骤5.6：一次累减还原原始数列。

现以序列的第一个数据点为例进行说明：

步骤5.7：起始的观测数据不断往前推移，实现预测值的迭代计算，得到多组参数a_i和b_i，其对应的预测值为

步骤5.8：计算各预测值

的权重w_i，根据权重计算预测值：

步骤6，利用步骤3和4所解算出的形变位移量、形变速度和形变加速度的实测值，根据预定的预警规则，对监测目标的当前状态进行实测预警；同时，利用步骤5所预测出的形变位移量、形变速度和形变加速度的预测值，根据预定的预警规则，对监测目标的后续状态进行预测预警。

上述预定的预警规则根据具体情况进行人为设定，如可以采用将形变位移量、形变速度和形变加速度进行加权的方式综合后进行预警判断，也可以采用与和或的逻辑关系将形变位移量、形变速度和形变加速度综合后进行预警判断等。

如利用监测目标综合位移X_all(k)和监测目标姿态进行实时预警等级G_lealdm判断，一般，地表水平位移超过15mm，垂直位移超过20mm，可判断为蓝色等级；地表水平位移超过30mm，垂直位移超过35mm，可判断为黄色等级；地表水平位移超过40mm，垂直位移超过45mm，可判断为橙色等级；地表水平位移超过50mm，垂直位移超过55mm，可判断为红色等级。实际等级划分的阈值须根据地质条件进行设置。

如利用形变预测量

速度、加速度等进行等级判断，预测位移量与预测姿态的等级G_prem判断和实时位移量与姿态的等级判断一样。速度切线角和加速度的预警等级G_preva的判断：加速度小于0，且速度切线角β小于45°时，可判断G_preva为蓝色等级；当加速度等于0或0值附近震荡，且速度切线角为45°<β<80°时，可判断G_preva为黄色等级；当加速度为较大正值，且速度切线角为80°≤β<85°时，可判断G_preva为橙色等级；当加速度为较大正值并增加，且速度切线角为85°≤β，约89°下滑时，可判断G_preva为红色等级；

如以加权方式进行后续状态预警判断时，预测位移量预警等级权重为w_prem，预测速度切线角和加速度的预警等级权重为w_preva，则监测目标的后续状态预测的预警等级G_pregrade为：

G_pregrade＝w_premG_prem+w_prevaG_preva

其中G_prem为预测位移量和预测姿态预警等级，G_preva为速度切线角和加速度的预警等级。

步骤7，滑坡形变监测网内的雨量计采集降雨量和/或裂缝传感器采集裂缝位移量，同时结合形变位移量、形变速度和形变加速度的实测值和预测值，根据预定的综合预警规则，对监测目标的进行综合预测预警。

利用降雨量、裂缝位移量等滑坡因子，进行预警等级判断，其中预定的预警规则根据具体情况进行人为设定。在本发明中，利用雨量计测量降雨量，采集降雨量L_r信息，降雨量滑坡预警因子w_r与等级G_r及降雨的天数d_r有着密切的关系。等级G_r越高，则w_r越大。具体的解算关系式，须根据当地降雨量和滑坡数量统计率而定。统计降雨天数d_r，分析滑坡与降雨天数d_r之间的滑坡预警因子w_rd：

其中：设监测区周围地区历史滑坡数为x，x_i为x的第i个值，

为x的均值；设监测时刻前7天降雨量为y，y_i为y的第i个值，

为y的均值。

降雨量的预警等级G_rain判断如下：当w_rd<0.3且L_r≤20mm时，G_rain可判断为蓝色等级，当0.3<w_rd<0.5且L_r≤30mm时，G_rain可判断为黄色等级；当0.5<w_rd<0.8且L_r≤40mm时，G_rain可判断为橙色等级；当w_rd>0.8且L_r≤90mm时，G_rain可判断为红色等级。

选用可以长期埋设在监测目标地表裂缝下的硅微式倾斜仪。根据硅微式倾斜仪的读数的大小计算监测目标倾角变化量：

△θ＝K₁F_i+K₂F² _i+K₃F_i ³-(K₁F₀+K₂F₀ ²+K₃F₀ ³)

其中：△θ为监测目标倾角变化量；K_i为倾斜仪的三次方拟合参数，i＝1,2,3，i为次方项；F_i为倾斜仪测量读数；F₀为倾斜仪基准读数。计算水平位移：

△x(k)＝Lsin△θ

其中，L为导轮间距，x_l为裂缝水平距离，d为倾斜仪测量次数。计算垂直位移：

△z(k)＝Lcos△θ

其中，z_l为垂直位移。

地表裂缝的内部位移量的预警等级G_fissure判断如下：当地表裂缝的内部位移量x_l≤10mm且Z_l≤15mm时，G_fissure可判断为蓝色等级，当10<x_l≤20mm且15<Z_l≤25mm时，G_fissure可判断为黄色等级；当20<x_l≤40mm且25<Z_l≤45mm，G_fissure可判断为橙色等级；当x_l>40mm且Z_l>45mm时时，G_fissure可判断为红色等级。实际等级划分的阈值须根据地质条件进行调整。

步骤8，根据上述各个预警等级进行实时综合预警及后续状态综合预警的等级判断，并设立的容错机制，进行综合预警。

由步骤6中监测目标综合位移X_all(k)和监测目标姿态实时预警等级G_lealdm、步骤7降雨量滑坡预警等级G_rain和地表裂缝内部位移量的预警等级G_fissure加权解算，得到实时综合预警等级G_lealgrade为：

G_lealgrade＝w_lealdmG_lealdm+w_rain1G_rain+w_fissure1G_fissure

其中，w_lealdm为G_lealdm的权重因子，w_rain1为降雨量滑坡预警等级在实时状态的权重因子，w_fissure1为地表裂缝内部位移量预警等级在实时状态的权重因子。

由步骤6监测目标的后续状态预测的预警等级G_pregrade、步骤7降雨量滑坡预警等级G_rain和地表裂缝内部位移量的预警等级G_fissure加权解算，得到实监测目标的后续状态预测的综合预警等级G_pg为：

G_pg＝w_pregradeG_pregrade+w_rain2G_rain+w_fissure2G_fissure

其中，w_pregrade为G_pregrade的权重因子，w_rain2为降雨量滑坡预警等级在后续状态预测的权重因子，w_fissure2为地表裂缝内部位移量的预警等级在后续状态预测的权重因子。

容错机制设置人为设定。容错机制设置为：当某一时刻出现预警等级变化时，保存预警结果，查询其前后若干时刻的记录是否满足预警等级变化要求，若持续监测后都满足，则进行预警，若不满足则重新判断。

本发明对滑坡的形变预警应用RTK多模高精度定位技术，并依据滑坡体(监测目标)形变的位移量、加速度、速度及当地降雨量等因素进行综合预警。预警级别的阈值因不同的地质、降雨量、监测目标姿态等不同而不同。利用预测算法，根据监测目标的位移形变量解算滑坡监测目标的未来形变趋势，结合裂缝计等数据可实现滑坡发生的提前警报。预警等级可分为蓝色、黄色、橙色、红色四个等级，滑坡形成程度逐级增加。

需要说明的是，尽管以上本发明所述的实施例是说明性的，但这并非是对本发明的限制，因此本发明并不局限于上述具体实施方式中。在不脱离本发明原理的情况下，凡是本领域技术人员在本发明的启示下获得的其它实施方式，均视为在本发明的保护之内。

Claims

1.一种滑坡形变综合预警方法，其特征是，具体包括步骤如下：

步骤6、利用步骤2所解算出的形变位移量、形变速度和形变加速度的实测值，根据预定的实测预警规则，对监测目标的当前状态进行实测滑坡等级预警；同时，利用步骤5所预测出的形变位移量、形变速度和形变加速度的预测值，根据预定的预测预警规则，对监测目标的后续状态进行预测滑坡等级预警。

2.根据权利要求1所述的一种滑坡形变综合预警方法，其特征是，步骤5的具体过程为：

3.根据权利要求1所述的一种滑坡形变综合预警方法，其特征是，步骤3的具体过程为：

步骤3.1、提取步骤2所得各个监测点的形变位移量的实测值；

4.根据权利要求1所述的一种滑坡形变综合预警方法，其特征是，步骤4的具体过程为：

5.根据权利要求1所述的一种滑坡形变综合预警方法，其特征是，还进一步包括如下步骤：

步骤7、滑坡形变监测网内的雨量计采集降雨量采集降雨量和裂缝传感器采集裂缝位移量，并将降雨量和裂缝位移量与步骤4和步骤5所得到的形变位移量、形变速度和形变加速度的实测值和预测值进行加权运算，并根据预定的综合预警规则，对监测目标的进行综合滑坡等级预警。

6.实现权利要求1所述方法的一种滑坡形变综合预警***，其特征是，包括滑坡形变监测网和综合预警单元；所述滑坡形变监测网由卫星定位基准站、3个以上的卫星定位监测站；3个以上的卫星定位监测站分布在监测目标的地表的不同监测点处；每个卫星定位监测站分别与卫星定位基准站连接，实现RTK定位解算；所有卫星定位监测站的定位结果均传输到综合预警单元。

7.根据权利要求6所述的一种滑坡形变综合预警***，其特征是，卫星定位监测站为GPS卫星定位接收机、BDS卫星定位接收机和/或GNSS卫星定位接收机。

8.根据权利要求6所述的一种滑坡形变综合预警***，其特征是，滑坡形变监测网还进一步包括雨量计和2个以上的裂缝传感器；雨量计分布在监测目标的地表上；2个以上的裂缝传感器分布在监测目标的地表的不同监测点处；雨量计和所有裂缝传感器的监测数据均传输到综合预警单元。

9.根据权利要求8所述的一种滑坡形变综合预警***，其特征是，卫星定位监测站的数量与裂缝传感器的数量相同；此时，监测目标的地表的每个监测点处均设有1个卫星定位监测站和1个裂缝传感器。