基于自适应卡尔曼滤波的地表移动变形预报研究

陈长坤,焦宝文,乔方,石长伟,肖明,苏迪

(1.安徽理工大学 测绘学院,安徽 淮南 232001； 2.广东省重工建筑设计院有限公司,广东 广州 510700)

0 引 言

卡尔曼滤波通过建立状态方程和量测方程描述系统的动态过程,依据滤波增益矩阵的变化,从量测数据中定量识别和提取有效信息,修正状态参量,无须存储各个不同时刻的观测数据,便于实时数据处理以及预报[1-2]。因此,卡尔曼滤波广泛应用于各种工程变形、滑坡监测、大坝监测等动态测量系统中[3-8]。卡尔曼滤波的应用需要动态系统的数学模型和噪声的先验知识,但在许多条件下,它们是未知的或近似已知的,即需研究动态系统的数学模型和系统中相关噪声的统计学特性。动态系统的数字模型和相关噪声统计特性两者的先验信息一般总是存在各种误差,自适应卡尔曼滤波正是为了克服这一缺点而提出的。它可以解决滤波的发散问题,从而减弱模型误差的影响,使滤波结果更接近于真实值[9-10],从而更好地预报各种工程变形例如地表移动变形监测[11]。

本文利用GNSS CORS地表移动自动化监测系统[12]对地表移动变形进行实时监测和分析。GNSS CORS在变形监测中具有实时、高效、高精度、自动化程度高的优点,现在已经广泛应用在地铁、大坝、桥梁以及滑坡监测、矿山开采沉陷等不同领域。对GNSS CORS连续运行实时监测站监测点的位置坐标序列进行实时滤波预报，并对比标准卡尔曼滤波预报值、自适应卡尔曼滤波预报值、实测值，对GNSS CORS地表移动实时监测进行实时预报分析。

1 标准卡尔曼滤波模型

1.1 卡尔曼滤波状态方程

设某一GNSS CORS连续运行监测网由n个监测站构成,可测量得到监测站监测点的3维位置坐标序列(本文直接得到3维位置坐标序列是适用于矿区的BJ-54坐标系下的高斯平面直角坐标和正常高),将t时刻(历元)的3维位置坐标及其速率构成状态向量。GNSS CORS连续运行监测站i在时刻(历元)t的位置坐标为ξi(t),瞬时速率为λi(t),瞬时加速率为Ωi(t),可将瞬时加速率看作一种随机干扰,则ξi(t)、λi(t)、Ωi(t)有以下微分关系:

(1)

记i点的状态向量为

=Xi(t)Yi(t)Zi(t)

则该监测点的状态方程为

(2)

把n个监测点的状态方程组合可得全网的状态方程为

Xk+1=Φk+1,kXk+Γk+1,kΩk.

(3)

1.2 卡尔曼滤波观测方程

测量获得第i个GNSS CORS监测站监测点的3维位置序列Li,k+1为观测值,某一监测站监测点在第k+1历元的观测方程为

Li,k+1=ξi,k+1+Δti,k+1λi,k+1+Δi,k+1,

(4)

式中:Δti,k+1=ti,k+1-tk+1,ti,k+1为Lij的观测时刻;tk+1为本次k+1历元观测的中心时刻,在地表移动变形观测的过程中,短时间内(数个历元,在本文中可以忽略不计)地表的移动变形速率可以忽略不计,可得全网的观测方程为

Lk+1=Bk+1Xk+1+Δk+1.

(5)

状态方程和观测方程两者共同组成GNSS CORS连续运行实时监测网的卡尔曼滤波模型:

Xk=Φk,k-1Xk-1+Γk,k-1Ωk-1,

Lk=BkXk+Δk,

(6)

式中:Φk,k-1为k-1到k历元的系统一步转移矩阵;Γk,k-1为系统噪声矩阵;Ωk-1为k-1历元的系统噪声;Bk为k历元系统的观测矩阵;Δk为k历元系统的观测噪声;Xk为k历元系统待估状态参数;Lk为k历元系统观测向量矩阵。Xk,Lk均为GNSS监测站监测点的3维位置和速度向量。

1.3 卡尔曼滤波算法

卡尔曼滤波的函数模型由式(6)给出,其随机模型为

E(Ωk)=0,

E(Δk)=0,

Cov(Ωk,Ωj)=DΩ(k)δkj,

Cov(Δj,Δj)=DΔ(k)δkj,

其中:DΩ(k)为系统的动态噪声方差阵,是非负定方差阵;DΔ(k)为系统的观测噪声方差阵,Δ为正定方差阵;δkj为Kronecker函数:

动态噪声与观测噪声完全不相关,即:

Cov(Ωk,Δj)=0

用标准卡尔曼滤波方程处理GNSS监测站各历元位置序列,其滤波方程为

X(k/k)=X(k/(k-1))+Jk

[Lk-BkX(k/(k-1))] ,

(7)

DX(k/k)=(I-JkBk)DX(k/(k-1)) .

(8)

X(k/k)为状态滤波方程;DX(k/k)滤波误差协方差阵,其中:

X(k/(k-1))=Φk,k-1X((k-1)/(k-1))

(9)

(10)

(11)

式中:X(k/(k-1))为一步预测值;Jk为滤波增益矩阵;Lk-BkX(k/(k-1))为预测残差;X(k/(k-1))为预测方差阵;DX(k/(k-1))为预报误差协方差阵。

式(7)和式(8)即是GNSS CORS连续运行监测站监测点的三维位置序列的标准卡尔曼滤波的递推计算公式,在确定初始状态后即可求得监测点三维位置的滤波值,X(k/k)即预报值。

2 方差补偿自适应卡尔曼滤波模型

自适应卡尔曼滤波就是在对观测数据进行递推滤波的同时,不断地对未知的或不确定的模型参数和噪声的统计特性进行适当的估计和修正,以减小模型误差,使滤波结果更接近于真实值[11]。方差补偿自适应卡尔曼滤波是利用预测残差对动态噪声的协方差向量进行修正,计算出更接近实际的状态向量。它的公式推导过程如下:

由第1章中标准卡尔曼滤波方程式(7)～(11)

假定{Ωk}和{Δk}为正态序列,X0为正态向量,{Ωk}和{Δk}为互不相关的零均值白噪声序列。

定义i步的预测残差为

Vk+i=Lk+i-L(k+i)/k,

(12)

式中:Lk+i、L(k+i)/k为第k+i历元的观测值和它的最佳预测值;Vk+i为预测残差。

Lk+i=Bk+iΦ(k+i)/kXk+Δk+i,

则Vi+1的方差阵Dvv为

(13)

记

Bk+iΦ(k+i)/kΓγ/(γ-1)=A(k+i,γ)=[ahj(k+i,γ)],

式中,r=1,…,N;k=1,…,n,上标k+i,r表示与k+i,r有关。假定DΩr-1Ωr-1在观测时间段tk+1,tk+2,…,tk+N上为常值对角阵,即

(14)

并记

有

记

(15)

其中ηk+i为零均值随机变量,r=1,…,N.

令

(16)

又记

E=[Ek+1,…,Ek+N]T,

η=[ηk+1,…,ηk+N]T,

(17)

则有

E=AdiagDΩΩ+η.

(18)

式(18)为关于diagDΩΩ的线性方程组。当N≥r时,有唯一解。记diagDΩΩ的LS估计为

diagDΩΩ=(ATA)-1ATE.

(19)

根据上述各式求得任意时间段长度的DΩΩ,并作为动态噪声协方差阵的实时估计。

由此得出的协方差可以对模型进行实时改正。自适应卡尔曼滤波在进行滤波的同时,能够实时地按照相应的自适应方法对模型进行修改,有效地降低滤波过程中出现的发散现象。

3 算例分析

3.1 试验简介

为了验证本文提出的自适应卡尔曼滤波预报的正确性,以淮南矿区朱集东矿1222(1)工作面上的GNSS CORS自动化监测系统进行试验。该系统由1个基准站(CZJDM)和2个连续运行实时监测站(CORS1和CORS2)组成。

本文选用GNSS CORS1、CORS2连续运行监测站在2017年6月的监测数据进行横向分析,选取两个实时监测站每天1:00～2:00时段各历元的标准卡尔曼滤波预报值和自适应卡尔曼滤波预报值(共60组数据,每组数据3600个历元)进行滤波预报精度评价。

选取GNSS CORS1连续运行监测站2017年6月18日凌晨1:00～2:00的数据(观测数据采样率为1 s,共3600个历元)进行纵向分析,分别采用标准卡尔曼滤波算法和自适应卡尔曼滤波算法对监测点的3维位置序列进行滤波预报,并获取该实时监测站的实测3维位置序列,然后分别对比分析标准卡尔曼滤波、自适应卡尔曼滤波预报值与实测值的偏差。

选取GNSS CORS1、CORS2连续运行监测站2017年6月每周第一天1:00～5:00的监测数据进行内符合精度分析(共8组数据,即选取8个时间段,每组数据14400个历元),分别采用标准卡尔曼滤波算法和自适应卡尔曼滤波算法对监测点的3维位置序列进行滤波预报,分别求出相应滤波值的内符合精度,并对比分析标准卡尔曼滤波、自适应卡尔曼滤波预报值的内符合精度。

3.2 数值结果与分析

横向数据分析中,对各个预报时段所有观测历元的平面坐标中误差和高程中误差进行精度评定;通过该时段中每个历元3维位置序列的预报偏差m求得该预报时段i整体的预报中误差。

(20)

式中:mP为滤波预报的平面坐标中误差；mH为滤波预报的高程中误差；k为该预报时段的历元个数；mx、my、mH分别为该预报时段历元k的预报偏差。

表1 滤波预报平面精度信息统计

表2 滤波预报高程精度信息统计

从表1和表2中可以看出:对60个时段预报的精度进行对比分析,自适应卡尔曼滤波预报相较于标准卡尔曼滤波,其平面位置精度和高程位置精度得到大幅度提升(大约提高了2倍左右),同时精度分布得更为均匀,即自适应卡尔曼滤波预报后的成果更为稳定可靠。

纵向数据分析中,表3～表5以及图1～图3给出了本次试验中标准卡尔曼滤波和自适应卡尔曼滤波在GNSS CORS连续运行监测站一个预报时段的平面X、Y坐标序列偏差和高程序列偏差信息。

表3 X坐标滤波预报偏差统计信息

从表3可以看出，对于平面位置X坐标序列，标准卡尔曼滤波预报最大偏差值3.21 mm，平均偏差值0.53 mm，偏差值小于0.5 mm的比率占57.67%；自适应卡尔曼滤波预报最大偏差值为0.80 mm，平均偏差值0.18 mm，预报偏差值相较于标准卡尔曼滤波平均缩减四分之一，偏差值小于0.5 mm的比率为99.94%.

表4 Y坐标滤波预报偏差统计信息

从表4可以看出，对于平面位置Y坐标序列。标准卡尔曼滤波预报最大偏差值为1.83 mm，平均偏差值0.33 mm，偏差值小于0.4 mm的比率占67.81%；自适应卡尔曼滤波预报最大值为0.80 mm，平均偏差值0.18 mm，预报偏差值相比较于标准卡尔曼滤波平均缩减四分之一，偏差值小于0.4 mm的比率为99.97%.

表5 高程滤波预报偏差统计信息

从表5看出，对于高程序列。标准卡尔曼滤波预报最大偏差值为8.27 mm，平均偏差值1.45 mm，偏差值小于1 mm的比率占43.39%；自适应卡尔曼滤波预报最大偏差值为1.29 mm，平均偏差值0.58 mm，预报偏差值相比较于标准卡尔曼滤波平均缩减四分之一，偏差值小于1 mm的比率为99.94%.由图3可知高程序列偏差值分布更为均匀，用自适应卡尔曼滤波预报的成果稳定性更高。

3.3 内符合精度分析

对于某一时段,若观测了k个历元,取各历元观测值的平均值,即可求得该监测站的平均位置,并根据各观测值与平均值之差,对该监测站的内符合精度进行评定,即

(21)

(22)

从表6可以看出2个GNSS CORS 连续运行监测站(共8个时段)的标准卡尔曼滤波预报值的平面内符合精度在4.3 ～6.2 mm之间,平均值5.4 mm,高程内符合精度在7.0 ～9.1 mm之间平均值为8.3 mm;从内符合精度可知,GNSS CORS连续运行监测站的平面位置精度满足开采沉陷精度要求(根据《煤矿测量规程》,为了保证解算的开采沉陷的关键参数的精度,要求相邻两期间平面点位中误差≤±20 mm,最弱的高程中误差≤±10 mm,这就要求一次测量平面点位相对中误差≤±14 mm,最弱点高程中误差≤±7 mm),但是高程内符合精度比平面位置测量精度低一倍,标准卡尔曼滤波预报值大于7 mm,不满足开采沉陷的精度要求。而自适应卡尔曼滤波预报值的平面内符合精度在3.2 ～4.3 mm之间,平均值为3.9 mm,高程内符合精度在5.9 ～7.4 mm之间,平均值6.6 mm;从内符合精度可知,GNSS CORS连续运行监测站的平面位置精度完全满足开采沉陷精度要求,高程方向均基本满足开采沉陷精度要求。且自适应卡尔曼滤波预报值相比较于标准卡尔曼滤波预报值,平面位置精度最低提高14.0%,最高提高37.9%,平均提高26.6%,高程方向精度最低提高10.0%,最高提高33%,平均提高20.1%.

从以上分析中可以看出在GNSS CORS连续运行监测站的位置坐标序列的预报中,其自适应卡尔曼滤波预报值相较标准卡尔曼滤波预报精度得到大幅提高,自适应卡尔曼滤波更能真实体现地表移动变形动态数据的实时性;对于标准卡尔曼滤波预报偏离实测值大,并且分布不均匀,自适应卡尔曼滤波预报偏差大幅度降低,预报偏差分布更为均匀,内符合精度也得到提高,使滤波过后的每个历元的测量成果更为稳定和可靠,这为更进一步获得稳定的移动变形量提供了坚实基础。

表6 滤波内符合精度

4 结束语

本文中的自适应卡尔曼滤波算法,对GNSS CORS连续运行实时监测站各测点的3维位置序列进行预报,根据该研究区域的实测数据表明,相比于标准卡尔曼滤波算法,其预报值相较于实测值的偏差值缩小约四分之一,预报精度提高了约2倍,同时分布更为均匀,表明本文采用的自适应卡尔曼滤波对GNSS变形监测具有很好的剔除噪声作用;平面位置内符合精度和高程内符合精度均提高了20%左右,自适应卡尔曼滤波预报的成果更为稳定和可靠,为获得准确的实时变形量提供基础。

但是本文中的自适应卡尔曼算法存在一定问题,还有待进一步改进。由图3可知虽然自适应卡尔曼滤波预报偏差值相比与标准卡尔曼滤波偏差值明显减小,但是自适应卡尔曼滤波预报偏差值最终趋近于某一值,而不是随着历元增加预报偏差值变更小。在某些历元标准卡尔曼滤波的预报偏差值会小于自适应卡尔曼滤波,这正是标准卡尔曼滤波预报的随机性导致的。本文的自适应卡尔曼滤波算法对于GNSS CORS连续运行实时监测站监测的地表移动变形监测的精度是足够了。

另外,由于该自适应卡尔曼滤波算法通过建立状态方程和量测方程来描述系统的动态过程,依据滤波增益矩阵的变化,从测量数据中定量识别和提取有效信息,修正状态参量,以此对实时监测数据进行处理及预报。但由于GNSS CORS连续运行实时监测站本身存在系统误差,或者是其他因素导致的粗差,该模型并不能探测并剔除粗差,这也是后续需要考虑的问题。

[1] 崔希璋,於宗俦,陶本藻,等.广义测量平差[M].2版.武汉:武汉大学出版社,2009.

[2] 尹晖.时空变形分析与预报的理论和方法[M].北京:测绘出版社，2002.

[3] 余学祥,吕伟才.GPS监测网动态数据处理抗差Kalman滤波模型[J].中国矿业大学学报,2000,29(6):553-557.

[4] 余学祥,吕伟才.抗差卡尔曼滤波模型及其在GPS监测网中的应用[J].测绘学报,2001,30(1);27-31

[5] LU W C,XU S Q.Reasearch on Kahnan filtering algorithm for deformation informantion series of similar single-difference model[J].Journal of China University of Mining and Technology,2004,14(2):189-194.

[6] 余学祥,吕伟才.GPS监测网动态数据处理抗差Kalrnan滤波模型[J].中国矿业大学学报,2000,29(6):553-557.

[7] 杨元喜.动态Kalman滤波模型误差的影响[J].测绘科学,2006,31(1):17-18

[8] 尹航,朱纪洪,周池军,等.基于Kalman预报观测器的增量动态逆控制[J].清华大学学报(自然科学版),2011.

[9] 贾志军,单甘霖,程兴亚,等.GPS动态定位中的自适应扩展卡尔曼滤波算法[J].军械工程学院学报,2001(2):39-43.

[10] 高雅萍,张勤.方差补偿自适应卡尔曼滤波在GPS滑坡监测中的应用研究[J].水土保持研究,2008(12): 150-152.

[11] 张福荣.自适应卡尔曼滤波在变形监测数据处理中的应用研究[D].西安:长安大学，2009.

[12] 余学祥.煤矿开采沉陷自动化监测系统[M].北京:测绘出版社，2014.