﻿ 抗差自适应卡尔曼滤波模型及其在塌陷区监测中的应用
 文章快速检索 高级检索
 大地测量与地球动力学  2019, Vol. 39 Issue (12): 1265-1269  DOI: 10.14075/j.jgg.2019.12.010

### 引用本文

HE Han, TAO Tingye, FENG Jiaqi, et al. Adaptive Robust Kalman Filtering Model and Its Application in Subsidence Area Monitoring[J]. Journal of Geodesy and Geodynamics, 2019, 39(12): 1265-1269.

### Foundation support

Natural Science Foundation of Anhui Province, No. 808085MD105.

### Corresponding author

TAO Tingye, PhD, associate professor, majors in GNSS technology theory and application, E-mail:czytty@163.com.

### 第一作者简介

HE Han, postgraduate, majors in GNSS deformation monitoring automation, E-mail: hh19931121@sina.com.

### 文章历史

1. 合肥工业大学土木与水利工程学院，合肥市屯溪路193号，230009

1 地表快速沉降动力学模型与测量方程的建立

 $\left[ {\begin{array}{*{20}{c}} {{X_k}}\\ {X_k^\prime } \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 1&{\Delta {t_{k,k - 1}}}\\ 0&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{X_{k - 1}}}\\ {X_{k - 1}^\prime } \end{array}} \right] + {\mathit{\boldsymbol{\omega }}_k}$ (1)

 ${\mathit{\boldsymbol{\omega }}_k} = \left[ {\begin{array}{*{20}{c}} {\frac{1}{2}\Delta t_{k,k - 1}^2}\\ {\Delta {t_{k,k - 1}}} \end{array}} \right]{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_{k - 1}}$ (2)

 $L\left( {{t_k}} \right) = \left[ {\begin{array}{*{20}{c}} 1&0 \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {{X_k}}\\ {X_k^\prime } \end{array}} \right] + {\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_k}$ (3)

 ${\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{k,t - 1}} = \left[ {\begin{array}{*{20}{l}} 1&1\\ 0&1 \end{array}} \right],{\mathit{\boldsymbol{A}}_k} = \left[ {\begin{array}{*{20}{l}} 1&0 \end{array}} \right],{\mathit{\boldsymbol{W}}_k} = \left[ {\begin{array}{*{20}{c}} {0.5}\\ 1 \end{array}} \right]$ (4)
2 滤波模型

2.1 对观测量进行抗差处理 2.1.1 稳定沉降状态采用抗差M估计

M估计由极大似然估计推出，Huber将其定义广义化[8]，最终得到：

 $\sum\limits_{i = 1}^n {\frac{{\partial \rho \left( {{\mathit{\boldsymbol{l}}_i},\mathit{\boldsymbol{\hat X}}} \right)}}{{\partial \mathit{\boldsymbol{\hat X}}}}} = 0$ (5)

 $\mathit{\boldsymbol{V}} = \mathit{\boldsymbol{A\hat X}} - \mathit{\boldsymbol{L}}$ (6)

 ${\mathit{\boldsymbol{v}}_i} = {\mathit{\boldsymbol{a}}_i}\mathit{\boldsymbol{\hat X}} - {\mathit{\boldsymbol{l}}_i}$ (7)

 $\sum\limits_{i = 1}^n {\mathit{\boldsymbol{a}}_i^{\rm{T}}} {\mathit{\boldsymbol{P}}_i}\frac{{\rho \left( {{\mathit{\boldsymbol{v}}_i}} \right)}}{{{\mathit{\boldsymbol{v}}_i}}}{\mathit{\boldsymbol{v}}_i} = 0$ (8)

$\omega_{i}=\frac{\rho\left(\boldsymbol{v}_{i}\right)}{\boldsymbol{v}_{i}}$为权因子, Pi=Piωi为等价权函数，将式(7)改写成矩阵形式：

 ${\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{\bar PV}} = 0$ (9)

 $\mathit{\boldsymbol{\hat X}} = {\left( {{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{\bar PA}}} \right)^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{\bar PL}}$ (10)

 $\omega \left( v \right) = \left\{ {\begin{array}{*{20}{l}} {1,}&{\left| \mathit{\boldsymbol{v}} \right| < 1.5\sigma }\\ {\frac{1}{{\left| \mathit{\boldsymbol{v}} \right| + k}},}&{1.5\sigma < \left| \mathit{\boldsymbol{v}} \right| < 2.5\sigma }\\ {0,}&{\left| \mathit{\boldsymbol{v}} \right| \le 2.5\sigma } \end{array}} \right.$ (11)

 ${{\mathit{\boldsymbol{\bar P}}}_i} = \left\{ \begin{array}{l} 1,\;\;\;\;\;\;\;\;\left| \mathit{\boldsymbol{v}} \right| < 1.5\sigma \\ {\mathit{\boldsymbol{P}}_i} \cdot \frac{1}{{\left| \mathit{\boldsymbol{v}} \right| + k}},\;\;\;\;\;1.5\sigma \le \left| \mathit{\boldsymbol{v}} \right| < 2.5\sigma \\ 0,\;\;\;\left| \mathit{\boldsymbol{v}} \right| \ge 2.5\sigma \end{array} \right.$ (12)

 $\left| {{{\mathit{\boldsymbol{\hat X}}}^k} - {{\mathit{\boldsymbol{\hat X}}}^{k - 1}}} \right| \le \varepsilon$ (13)

2.1.2 快速沉降状态检测与处理

 $U\left( t \right) = \sum\limits_{i = 1}^t {\sum\limits_{j = t + 1}^n {\rm{sgn}} } \left( {{x_j} - {x_i}} \right),1 \le t \le n$ (14)

 $K = \max \left| {U\left( t \right)} \right|$ (15)

t时刻为突变时刻，同时计算统计量：

 $P \approx 2\exp \left( {\frac{{ - 6{K^2}}}{{{n^3} + {n^2}}}} \right)$ (16)

2.2 对状态模型采用自适应因子调节

 $\begin{array}{*{20}{c}} {{a_k} = }\\ {\left\{ {\begin{array}{*{20}{l}} {1,\left| {\Delta {{\mathit{\boldsymbol{\tilde X}}}_k}} \right| \le {c_0}}\\ {\frac{{{c_0}}}{{\left| {\Delta {{\mathit{\boldsymbol{\tilde X}}}_k}} \right|}}{{\left( {\frac{{{c_1} - \left| {\Delta {{\mathit{\boldsymbol{\tilde X}}}_k}} \right|}}{{{c_1} - {c_0}}}} \right)}^2},{c_0} < \left| {\Delta {{\mathit{\boldsymbol{\tilde X}}}_k}} \right| \le {c_1}}\\ {0,\left| {\Delta {{\mathit{\boldsymbol{\tilde X}}}_k}} \right| > {c_1}} \end{array}} \right.} \end{array}$ (17)

 $\Delta {{\mathit{\boldsymbol{\tilde X}}}_k} = \left\| {{{\mathit{\boldsymbol{\tilde X}}}_k} - {{\mathit{\boldsymbol{\bar X}}}_k}} \right\|/\sqrt {\rm{tr}\left( {{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{{\bar X}_k}}}} \right)}$ (18)

2.3 抗差自适应卡尔曼滤波

 ${\mathit{\boldsymbol{V}}_{{{\bar X}_k}}} = {{\mathit{\boldsymbol{\hat X}}}_k} - {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{k,k - 1}}{{\mathit{\boldsymbol{\hat X}}}_{k - 1}}$ (19)

 ${\mathit{\boldsymbol{V}}_k} = {\mathit{\boldsymbol{A}}_k}{{\mathit{\boldsymbol{\hat X}}}_k} - {\mathit{\boldsymbol{L}}_k}$ (20)

 $\mathit{\boldsymbol{V}}_k^{\rm{T}}{{\mathit{\boldsymbol{\bar P}}}_k}{\mathit{\boldsymbol{V}}_k} + {a_k}\mathit{\boldsymbol{V}}_{{{\bar X}_k}}^{\rm{T}}{\mathit{\boldsymbol{P}}_{{{\bar X}_k}}}{\mathit{\boldsymbol{V}}_{{{\bar X}_k}}} = \min$ (21)

 $\mathit{\boldsymbol{\hat X}} = {\left( {\mathit{\boldsymbol{A}}_k^{\rm{T}}{{\mathit{\boldsymbol{\bar P}}}_k}{\mathit{\boldsymbol{A}}_k} + {a_k}{\mathit{\boldsymbol{P}}_{{{\bar X}_k}}}} \right)^{ - 1}}\left( {\mathit{\boldsymbol{A}}_k^{\rm{T}}{{\mathit{\boldsymbol{\bar P}}}_k}{\mathit{\boldsymbol{L}}_k} + {a_k}{\mathit{\boldsymbol{P}}_{{{\bar X}_k}}}{{\mathit{\boldsymbol{\bar X}}}_k}} \right)$ (22)

 图 1 滤波解算流程 Fig. 1 Flow chart of filtering solution
3 实例分析

3.1 实验数据及滤波初始值的确定

 图 2 自动化监测点位置示意图 Fig. 2 Position sketch map of the automatic monitoring points

3.2 数据处理与分析

 图 3 抗差卡尔曼滤波值与测量值及真值的比较 Fig. 3 Comparison of robust Kalman filtering values with measured value and true values

 图 4 抗差自适应卡尔曼滤波值与测量值及真值的比较 Fig. 4 Comparison of adaptive robust Kalman filtering values with measured values and true values
4 结语