GNSS/SINS组合的捷联式航空重力测量系统在实测过程中由于惯性器件误差随时间积累及外界因素干扰，严重影响重力扰动的估算精度^[1-4]。如何有效去除惯性器件误差等因素的影响、提高重力扰动的估算精度，一直是国内外研究的热点。20世纪90年代Jekeli等^[5]开始利用GPS与惯性导航系统组合的方式进行航空重力测量，水平和垂向分量的外符合精度可分别达到7~8 mGal、3 mGal。文献[6]采用波估算法对加拿大卡尔加里大学的数据进行处理，重复测线重力扰动水平分量的内符合精度为4~8 mGal，垂向为5~9 mGal，进一步用波数相关滤波方法处理后，水平和垂向分量的精度可达到2~4 mGal。本文对国产捷联式航空重力矢量测量系统的重复性精度进行测试，利用卡尔曼滤波、波数相关滤波、比力线性校正方法对实测数据进行处理，估计重复测线的重力扰动矢量并评定其内符合精度。

1 数学模型 1.1 GNSS/SINS组合重力测量的基本原理

捷联式航空重力测量的基本原理是：将载体的运动加速度从比力观测量中扣除，即可得重力。当地水平坐标系下某点处重力的表达式为^[7]：

$ {g^n} = {\dot v^n} - \mathit{\boldsymbol{C}}_b^n{f^b} + (2\omega _{ie}^n + \omega _{en}^n) \times {v^n} $

(1)

式中，上角标n、b分别代表当地水平坐标系和载体坐标系；gⁿ为所求的重力值；vⁿ、$ {\dot v^n} $分别为载体在某位置的运动速度和加速度；C_bⁿ为载体坐标系b到当地水平坐标系n的转换矩阵，f^b为原始比力观测值；(2ω_ieⁿ+ω_enⁿ)×vⁿ为由地球自转和载体运动引起的哥氏加速度。

由于重力可以表示成正常重力γⁿ与重力扰动δgⁿ之和，故式(1)又可以写成：

$ \delta {g^n} = {\dot v^n} - \mathit{\boldsymbol{C}}_b^n{f^b} + (2\omega _{ie}^n + \omega _{en}^n) \times {v^n} - {\gamma ^n} $

(2)

1.2 GNSS/SINS组合卡尔曼滤波算法

为抑制惯性器件误差随时间的积累，可利用卡尔曼滤波方法动态估计SINS误差并及时修正。15状态参数包括位置、速度和姿态误差以及陀螺漂移和加速度计偏置，记为X =[δr, δv, Ψ, b_g, b_a]。在当地水平坐标系下SINS的误差方程为^[8-9]：

$ \delta \dot r = \delta v $

(3)

$ \begin{array}{l} \delta \dot v = [{f^n}] \times \psi + \mathit{\boldsymbol{C}}_b^n\delta {f^b} - (2\omega _{ie}^n + \omega _{en}^n) \times \\ \;\;\;\;\;\;\delta v - (2\delta \omega _{ie}^n + \delta \omega _{en}^n) \times v - \delta {g^n} \end{array} $

(4)

$ \dot \psi = - \omega _{in}^n \times \psi + \delta \omega _{in}^n - \mathit{\boldsymbol{C}}_b^n\delta \omega _{ib}^b $

(5)

$ {\dot b_g} = 0 $

(6)

$ {\dot b_a} = 0 $

(7)

用式(3)~(7)可构建卡尔曼滤波系统的连续状态方程^{[2, 9-10]}：

$ \left[ {\begin{array}{*{20}{c}} {\delta \dot r}\\ {\delta \dot v}\\ {\dot \psi }\\ {{{\dot b}_g}}\\ {{{\dot b}_a}} \end{array}} \right] = \mathit{\boldsymbol{F}}\left[ {\begin{array}{*{20}{c}} {\delta r}\\ {\delta v}\\ \psi \\ {{b_g}}\\ {{b_a}} \end{array}} \right] + \mathit{\boldsymbol{GW}} $

(8)

式中，F为连续系统的状态转移矩阵，G为噪声系数阵，W为过程白噪声向量。F和G中的矩阵元素可参考文献[9]。

对式(8)进行离散化处理，得任意时刻i的状态参数为^[9]：

$ {\mathit{\boldsymbol{X}}_i} = {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{i, i - 1}}{\mathit{\boldsymbol{X}}_{i - 1}} + {\mathit{\boldsymbol{\omega }}_i} $

(9)

式中，X_i、X_i-1分别为i和i-1时刻的状态参数向量；Φ_{i, i-1}为离散化后的状态转移矩阵；ω_i为动力学噪声向量。

以GNSS和SINS的位置、速度之差作为观测量，则卡尔曼滤波系统的量测方程为^[11]：

$ {\mathit{\boldsymbol{Y}}_i} = {\mathit{\boldsymbol{H}}_i}{\mathit{\boldsymbol{X}}_i} + {\mathit{\boldsymbol{e}}_i} $

(10)

式中，Y_i=[r_GNSS-r_INS, v_GNSS-v_INS]^T，$ {\mathit{\boldsymbol{H}}_i} = \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{I}}_{3 \times 3}}}&{{0_{3 \times 3}}}&{{0_{3 \times 9}}} \end{array}}\\ {\begin{array}{*{20}{c}} {{0_{3 \times 3}}}&{{\mathit{\boldsymbol{I}}_{3 \times 3}}}&{{0_{3 \times 9}}} \end{array}} \end{array}} \right] $，e_i为观测量白噪声。

1.3 重力扰动矢量的误差补偿 1.3.1 波数相关滤波方法

波数相关滤波(WCF)方法可用于消除重复测线重力扰动(尤其水平分量)中存在的低频趋势项和系统偏差，处理过程如下。

1) 先将2条反向测线的重力信号通过傅里叶变换从时域转到频域，利用下式计算不同波数k所对应的频域数据F_x(k)、F_y(k)的相关系数：

$ {\rm{C}}{{\rm{C}}_k} = {\rm{cos}}(\Delta {\theta _k}) = \frac{{{F_x}(k) \cdot {F_y}(k)}}{{\left| {{F_x}(k)} \right|\left| {{F_y}(k)} \right|}} $

(11)

式中，Δθ_k为F_x(k)、F_y(k)之间的相位差，Δθ_k=θ_x-θ_y；“·”为标量积，“||”为矢量求模。

2) 给定一个阈值T并与上面求得的相关系数进行比较。如果CC_k大于T，则认为该成分是重力信号，予以保留，反之予以剔除^[12]：

$ {\bar F_x}(k) = \left\{ {\begin{array}{*{20}{c}} {{F_x}(k), }&{{\rm{C}}{{\rm{C}}_k} \ge T}\\ {0, }&{{\rm{C}}{{\rm{C}}_k} < T} \end{array}} \right. $

(12)

$ {\bar F_y}(k) = \left\{ {\begin{array}{*{20}{c}} {{F_y}(k), }&{{\rm{C}}{{\rm{C}}_k} \ge T}\\ {0, }&{{\rm{C}}{{\rm{C}}_k} < T} \end{array}} \right. $

(13)

3) 将处理后的频域信号经傅里叶逆变换转到时域，得到该条重复测线上的重力信号G：

$ G = \frac{1}{2}{\rm{IFFT}}({\bar F_x} + {\bar F_y}) $

(14)

1.3.2 比力线性校正方法

由于重复测线重力扰动的垂向分量受陀螺姿态角误差的影响较小，可利用下式对任意时刻比力的垂向分量依时间进行线性修正：

$ \Delta {f_i} = \frac{{\Delta \delta \bar g}}{{{t_e} - {t_s}}}({t_i} - {t_s}) $

(15)

$ {\tilde f_i} = {f_i} - \Delta {f_i} $

(16)

式中，Δδg为首尾两条重复测线上对应点处重力扰动垂向分量的差值均值；t_s、t_e分别为有效重复测线的起始与结束时刻；f_i、Δf_i、$ {\tilde f_i} $分别为t_i时刻比力垂向分量的原始观测值、改正数以及修正后的值。

1.4 重力扰动矢量的估算流程

用捷联式重力仪原始观测数据估算重力扰动矢量的基本流程如图 1所示。

1) 利用陀螺和加速度计的观测量对SINS进行导航解算，得到惯导在当地水平坐标系下的位置、速度、姿态角和比力矢量。

2) 对GNSS的原始载波相位观测量进行定位解算，然后经两次差分分别求得载体的运动速度和加速度^[13-15]。

3) 以GNSS和SINS的位置、速度之差作为卡尔曼滤波量测方程的观测量，根据式(8)、(10)进行GNSS/SINS组合卡尔曼滤波，估计SINS的误差并及时修正，并对当地水平坐标系下的比力进行校正：

$ {\tilde f^n} = {f^n} - \mathit{\boldsymbol{C}}_b^n{b_a} $

(17)

4) 用GNSS的位置、速度计算哥氏加速度和正常重力，用式(2)计算全部测线的重力扰动矢量。

5) 由于式(2)求得的重力扰动中含有高频噪声，而重力信息为低频信号，故需用低通滤波器去噪^[16-17]，得到减弱高频噪声后的重力扰动矢量。

6) 由步骤1)~5)求得的重复测线的重力扰动中存在低频趋势项(尤其是水平分量)和系统误差，可用WCF或比力线性校正方法予以减弱。

2 算例分析 2.1 试验概况

2015-06-15用我国自行研制的捷联式航空重力矢量仪在某地进行飞行试验，图 2为测线轨迹，图 3为其中的6条东-西向重复测线(分别记为L₁、L₂、…、L₆)，每条测线长度约117 km，飞行速度约60 m/s，飞行高度约1 550 m。GNSS的采样率为2 Hz，SINS的采样率为200 Hz。

2.2 重力扰动估计结果与精度分析

根据所用捷联式重力仪性能指标，表 1给出卡尔曼滤波各参数的初始中误差。由卡尔曼滤波估算得到的测线重力扰动矢量如图 4(a)、4(b)、4(c)所示。可以看出，重力扰动矢量中含有大量高频噪声，经250 s低通滤波器处理后如图 5(a)、5(b)、5(c)所示。因飞行速度约为60m/s，故滤波后的重力扰动矢量的半波长分辨率约为7.5 km。

从图 5(a)、(b)中分离出6条重复测线的重力扰动水平分量，见图 6，任意两条重复测线上对应点处的重力扰动水平分量的差值均值列于表 2(单位mGal)。从图 6和表 2可以看出，水平分量(尤其是北向分量)存在明显的偏差和低频趋势。在给定阈值为零的条件下，用WCF方法对飞行方向相反的测线进行处理，结果如图 7所示，较差统计于表 3(单位mGal)。对比图 6与图 7、表 2与表 3可知，重力扰动水平分量中残留的误差得到有效剔除，测线间的系统偏差明显减小；重力扰动东向分量的平均内符合精度从9.77 mGal提高到5.95 mGal，北向分量的平均内符合精度从9.18 mGal提高到3.38 mGal。

图 8为图 5(c)中分离出的6条重复测线的重力扰动垂向分量。可以看出，6条重复测线间存在明显的系统偏差。因为从L₁至L₆为按时间顺序飞行，相对L₁，系统偏差逐渐增大，至L₆最大，达到8.19 mGal。将L₁与L₆之间的差值均值依时间等比例分配，对当地水平坐标系下比力的垂向分量进行线性补偿，结果如图 9所示，较差统计于表 4(单位mGal)。可以看出，测线间的系统性偏差明显减小，最大差值由8.19 mGal减小到2.54 mGal，平均内符合精度由1.27 mGal提高到1.10 mGal。图 10给出进一步用WCF方法对重力扰动垂向分量处理的结果，较差统计见表 5(单位mGal)，平均内符合精度达到0.59 mGal。

在没有重复测线的情况下，可以从飞机起飞和降落两时刻起对GNSS/SINS组合观测数据分别利用卡尔曼滤波进行正、逆向解算，通过直接求差法求得两组重力扰动矢量。然后对单条测线的正、逆向重力扰动矢量用WCF方法进行处理，以提高重力扰动的估算精度^[2]。

3 结语

本文利用卡尔曼滤波、波数相关滤波、比力线性校正的方法对国产捷联式航空重力矢量仪的某次试验数据进行处理，通过6条重复测线重力扰动矢量的对比分析，得到如下初步结论：

1) 由于姿态测量误差对重力扰动水平分量的影响要远大于垂直分量，因此，重复测线重力扰动水平分量的内符合精度比垂向分量要低。

2) 重复测线重力扰动的水平分量经WCF处理后，其间的低频趋势项和系统偏差得到极大减弱，且内符合精度有所提高，东向分量从9.77 mGal提高到5.95 mGal，北向分量从9.18 mGal提高到3.83 mGal。

3) 重复测线重力扰动的垂向分量间存在明显的系统偏差，经过对比力矢量的垂向分量依时间进行线性补偿后，系统误差明显减弱，重力扰动垂向分量的平均内符合精度由1.27 mGal提高到1.10 mGal。进一步用WCF方法处理后，平均内符合精度达到0.59 mGal。