捷联航空重力测量中的卡尔曼滤波

2019-06-25蔡体菁王新宇

压电与声光 2019年3期

蔡体菁，颜颖，王新宇

(1.东南大学仪器科学与工程学院，江苏南京 210096；2.中国电子科技集团公司第二十六研究所，重庆400060)

0 引言

近十几年来航空重力测量技术发展迅速，国际上不断出现新一代航空重力仪，如GT-2A[1]。国内也有多家科研单位研制出二轴、三轴稳定平台海空重力仪、捷联式航空海洋重力仪[2-3]。在航空重力测量数据处理方面，海空重力仪GT系列均采用了卡尔曼滤波技术，但具体算法未公开报道。2015年，东南大学给出了处理平台式重力仪测量数据的平滑正、反卡尔曼滤波法[4]。近年来，一些作者对重力数据处理进行了研究与讨论[5-7]。本文针对国产捷联式航空重力测量系统，根据马尔可夫过程的重力异常模型，给出基于该模型的平滑综合卡尔曼滤波法，它把速度作为观测量，重力异常及其导数作为状态变量，进行组合导航和自由空间重力异常计算，并把计算结果与同机工作的GT-2A航空重力仪重力测量结果进行比较，两者测量误差小于0.85 mGal。

1 平滑综合的正反卡尔曼滤波

卡尔曼滤波是一种最小化方差估计的递推数据处理算法，它通过前一时刻的状态估计值和当前时刻的观测值，递推出当前时刻的状态估计值。平滑是利用当前时刻前、后的状态值计算出当前时刻的状态值。常用的平滑方法有固定区间平滑、固定点平滑和固定滞后平滑3种。本文采用固定区间平滑方法。

随机线性离散系统的状态方程和观测方程可写为

(1)

式中：Xi为ti时刻的状态向量；Φi+1，i为ti～ti+1时刻的状态转移矩阵；Bi为系统输入矩阵；ui为系统输入向量；Gi为系统噪声输入矩阵；wi为系统噪声向量；Zi为ti时刻的观测向量；Hi为观测矩阵；vi为观测噪声向量。wi和vi均是零均值高斯白噪声且相互独立，其噪声强度分别为Qi和Ri。

式(1)的正向卡尔曼滤波为

(2)

初始条件：

(3)

P0=P(t0)

(4)

式(1)的反向卡尔曼滤波为

(5)

初始条件：

(6)

(7)

式(1)的平滑综合滤波为

(8)

2 GPS/捷联惯性组合导航卡尔曼滤波方程

由于捷联式航空重力仪和运载体固连，其3个加速度计输入轴方向会随测量载体运动而不断发生变化，因此，要获得加速度计比力的垂向分量，就必须知道每一时刻运载体的姿态，即载体坐标系相对导航坐标系的航向角、俯仰角和横滚角。对于重力数据后处理而言，利用捷联式航空重力仪和GPS基站的信息，应用卡尔曼滤波进行GPS/捷联惯性组合导航计算是获得高精度运载体姿态最有效的方法。

GPS/捷联惯性组合导航系统卡尔曼滤波状态方程为

(9)

系统噪声阵G=[ωgx,ωgy,ωgz,ωax,ωay,ωaz]T，其中ωgx、ωgy、ωgz分别为3个方向陀螺仪自身的白噪声；ωax、ωay、ωaz分别为3个方向加速度计的一阶马氏过程驱动白噪声。式(9)的具体形式如下：

(11)

(12)

(13)

(14)

ε=εb+ωg

(15)

(16)

(17)

GPS/捷联惯性组合导航系统卡尔曼滤波观测方程为

z(t)=H(t)x(t)+v(t)

(18)

式中z=[LI-LG,λI-λG,hI-hG,VIE-VGE,VIN-VGN,VIU-VGU]T，下标I为捷联惯性导航计算值，G为GPS测量值。

3 重力异常卡尔曼滤波方程

捷联式航空重力仪的输出信息通常为载体坐标系下的3个加速度、3个角速度、温度和GPS移动站的星历、伪距、多普勒频移、载波相位等原始数据。利用以上这些信息及GPS基站的原始信息，能够计算得到测线上的自由空间重力异常。自由空间重力异常Δg可以通过对惯性导航比力方程的变形得到，其方程为

(19)

地球重力场重力异常模型可以用二阶马尔可夫过程来描述：

(20)

式(19)、(20)可写为卡尔曼滤波状态方程形式：

(21)

用差分GPS的速度信息作为卡尔曼滤波的观测信息，观测方程为

Vup=HX+ν

(14)

式中：Vup为垂向速度；H为测量矩阵；ν为测量噪声。

4 计算结果

东南大学开展重力测量数据处理研究的捷联式航空重力仪由高精度惯性测量单元、GPS接收机和数据采集器组成。惯性测量单元包含3个高精度单轴石英挠性摆式加速度计、3个高精度单轴激光陀螺仪、电源板、温控板、计算机板、I/F板等电子线路，数据采集器采集惯性测量单元和GPS接收机输出的信息，并保证信息精确同步。

根据捷联式航空重力仪的测量数据，获得测线上的重力异常，要进行3步计算：

1)差分GPS位置、速度计算。

2)GPS/捷联惯导组合导航解算。

3)重力异常解算。

航空重力测量数据处理流程如图1所示。

图1 航空重力测量数据处理流程图

为了检验平滑综合的正、反卡尔曼滤波及重力异常模型的正确性和有效性，利用了捷联式航空重力仪与俄罗斯的GT-2A航空重力仪同时飞行的测量数据。运用上述计算方法得到自由空间重力异常，并把计算结果与GT-2A航空重力仪软件处理的结果对比，对比结果如图2所示。调整后，两航空重力仪输出的重力异常符合精度约为0.8 mGal。调整前的误差主要是两个航空重力仪的前校基准点数据有差别，通常相差一个常数。由图可见，在测线的细节上，GT-2A的计算结果更丰富。