捷联式航空重力测量系统的数据处理方法比较与分析

修 睿，郭 刚，薛正兵，李东明，李海兵

（1.青岛海洋科学与技术试点国家实验室,青岛266237；2.北京航天控制仪器研究所,北京100039）

0 引言

航空重力测量是利用机载重力测量系统进行的空中重力测量，是重力测量的重要手段。目前，国际上的航空重力测量系统主要为平台式系统。近年来，国内自主研制出了基于捷联惯性技术的新型重力测量系统，其被称之为捷联式重力测量系统。捷联式航空重力测量系统采用了捷联惯性/GPS组合导航算法技术，可计算出重力传感器敏感轴与当地地理水平面的不垂直度并予以数学校正，即可建立与物理平台等价的数学平台，具有体积小、质量小、可靠性高的优点。

在飞行作业过程中，由飞机机体高速运动带来的高频振动不可避免地对重力传感器产生了高频干扰。为尽可能消除或减弱这些高频噪声的影响，需采用重力数据滤波方法对重力传感器输出进行处理，从而获得高精度的重力异常值。

目前，比较常用的滤波方法主要包括有限脉冲响应FIR低通滤波、无限脉冲响应IIR低通滤波、Kalman滤波、RC滤波、Gauss滤波等。本文采用FIR低通滤波、零相移滤波、标准Kalman滤波和正反Kalman滤波4种方法，结合实际飞行试验数据，对4种滤波方法的处理结果进行了比对。同时，捷联式重力测量系统的数学平台精度也直接影响着重力传感器敏感轴指向的校正精度，从而影响着最终的重力测量精度。本文对基于GPS辅助的SINS/GPS与基于DGPS辅助的SINS/DGPS 2种数学平台计算方法获得的重力测量数据结果进行了比对。

1 系统组成及数据处理过程

1.1 捷联式重力测量系统的测量原理

捷联式海空重力仪（SAG）将重力传感器固联集成于SINS/GPS组合导航系统，重力传感器为高精度石英挠性加速度计。捷联式海空重力仪与显控记录装置、UPS一起构成了捷联式海空重力测量系统。

重力传感器敏感轴、捷联惯性系统天向轴、载体天向轴三者平行，重力传感器敏感轴的输入信息包含了重力信息沿载体天向轴方向的分量及沿载体天向轴方向的运动加速度。通过SINS/GPS组合导航计算出载体姿态，即建立数学平台，将重力传感器输出进行指向校正，可获得垂直于当地地理水平面的重力与垂向运动加速度之和；再利用GPS系统测量信息估算出垂向运动加速度，即可获得重力信息；再经滤波处理滤除振动等产生的干扰噪声，以及厄缶改正、高度改正、正常重力改正等，最终可获得自由空间重力异常值。

1.2 数据处理过程

捷联式重力测量数据处理流程包括导航解算部分和重力异常提取部分，重力异常提取过程在导航解算的基础上进行。导航原始数据（陀螺、加速度计的原始脉冲）经误差补偿和纯惯性导航解算后，与GPS数据进行组合导航和滤波处理，获得纯惯性导航的姿态误差、速度误差、位置误差，修正得到东北天向的比力值和修正后的载体速度、位置信息。将载体位置信息进行2次差分处理后，可得到载体的运动加速度。结合差分GPS速度、位置信息和3个比力信息，可以得到重力异常粗值。经过滤波处理和各项改正后，可得到精确的自由空间重力异常值，其计算公式为[1]

式中，δg为自由空间重力异常值，gb为基准点参考重力值，为比力和比力初值，为运动加速度改正，δaF为空间改正，δaE为厄缶改正，γ0为正常椭球面上的重力值。

自由空间重力异常提取由事后处理软件完成，其流程如图1所示。

图1 自由空间重力异常提取处理流程图Fig.1 Flowchart of free-space gravity anomaly extraction and processing

2 重力数据处理方法介绍

2.1 FIR低通滤波

目前，FIR滤波器主要涉及窗函数设计法、频率抽样设计法和Chebyshev逼近法，这里主要对常用的窗函数设计法进行介绍。

有限脉冲响应FIR数字滤波器要求用有限长的单位冲击响应h（n）来逼近无限长的理想滤波器的单位冲击响应hd（n），最常用和有效的方法就是用1个长度为N的窗函数序列ω（n）来截取hd（n）的主要成分

实际上，窗函数设计FIR滤波器的核心是由给定的频率特性通过加窗确定有限长单位脉冲响应序列h（n）。 工程中常见的窗函数包括以下几种：矩形窗、三角窗、Hanning窗、Hamming窗、Blackman窗、Kaiser窗，其性能如表1所示[2-4]。

表1 常用窗函数的性能比较Table 1 Performance comparison of commonly used window functions

本文中用到的FIR滤波器采用Hanning窗，采样频率fs=1Hz，归一化截止频率fc=0.005Hz，阶数N=600。

2.2 零相移滤波

零相移滤波是指信号序列在经过滤波器滤波之后，信号序列相位不发生变化，系统函数的相位响应为0，滤波后效果更平滑且信号无位移。零相移滤波算法流程如图2所示，需进行正向、反向2次滤波。

假设输入信号序列为 [x（n）]，n=1，2，…，n0，[y（n）]为信号[x（n）]时域翻转后的序列，将 2个信号分别进行沿x轴、y轴的延拓后，进行双边Z变换，得到

图2 零相移算法原理图Fig.2 Schematic diagram of zero phase shift algorithm

由算法流程，设传递函数为H（z），则滤波器输出信号序列为

经过二次滤波和翻转输出的信号序列为

本文中零相移滤波器的采样频率fs=1Hz，归一化截止频率fc=0.005Hz，阶数N=600。

2.3 标准Kalman滤波

设tk时刻的被估计状态Xk受系统噪声序列Wk-1驱动，其驱动机理由下述状态方程描述[5-7]

对Xk的量测满足线性关系，量测方程为

式中，Φk，k-1为tk-1时刻至tk时刻的一步转移阵，Γk-1为系统噪声驱动阵，Hk为量测阵，Vk为量测噪声序列，Wk-1为系统激励噪声序列。

Kalman滤波方程为：

（1）状态一步预测

（2）状态估计

（3）滤波增益

（4）一步预测均方误差

（5）估计均方误差

式中，Pk为估计误差方差阵，Kk为加权矩阵，Pk/k-1为一步预测均方误差，Xk/k-1为状态一步预测，Rk为量测噪声方差阵。

2.4 正反Kalman滤波

正反Kalman滤波是连续进行2次方向相反的Kalman滤波，并将正向、反向2次滤波结果进行加权求和，能够充分利用正向Kalman滤波的状态估计量、一步状态估计、误差协方差矩阵等参数信息。该方法克服了标准离散Kalman滤波的精度缺陷问题，优化了滤波处理结果[8-10]。

正反Kalman滤波的正向过程即为标准离散型Kalman滤波。

反向过程的滤波方程为

预测过程方程为

加权过程方程为

正反Kalman状态量是正向状态估计量xk与反向状态估计量xf（k/k-1）的加权求和，权重为Xzf（k）和Pzf（k）。

3 飞行试验结果比对

3.1 飞行试验概况

采用上述FIR低通滤波、零相移滤波、标准Kalman滤波、正反Kalman滤波4种滤波方法，对SAG某架次的飞行试验数据进行了处理，并将处理结果进行了精度比对。

在本次试验过程中，飞机平稳飞行，飞行平均高度为800m，飞行平均速度为60m/s，飞行轨迹为8条东西重复测线、2条南北测线。测线结果如图3所示，图中的经纬度进行了人为偏移。

图3 东西重复测线和南北测线结果Fig.3 Results of East-West repeat survey lines and North-South survey lines

3.2 基于SINS/GPS、SINS/DGPS组合导航的结果比对

分别采用GPS辅助导航和DGPS辅助导航进行导航解算，2种组合导航解算得到的比力值、航姿、速度、位置在测线段的曲线分别如图4～图7所示，图中的经纬度进行了人为偏移。其中，红色曲线为GPS辅助导航的解算结果，蓝色曲线为DGPS辅助导航的解算结果。对2种组合导航东西重复线上的比力、速度、航姿信息进行均值和标准值的比对，其结果如表2所示。

图4 2种组合导航解算测线段比力值的结果比对Fig.4 Comparison of specific force value between two kinds of integrated navigation solution survey lines

图5 2种组合导航解算测线段航姿的结果比对Fig.5 Comparison of attitude angle between two kinds of integrated navigation solution survey lines

图6 2种组合导航解算测线段速度的结果比对Fig.6 Comparison of velocity between two kinds of integrated navigation solution survey lines

图7 2种组合导航解算测线段位置的结果比对Fig.7 Comparison of position between two kinds of integrated navigation solution survey lines

表2 2种组合导航测线段比力、姿态、速度的信息比对Table 2 Comparison of specific force，attitude angle and velocity between two kinds of integrated navigation solution survey lines

由图4～图7所示的2种组合导航的比力值、航姿、速度、位置结果的比对曲线可以看出，红色曲线和蓝色曲线均完全重合，即GPS辅助导航结果与DGPS辅助导航的结果基本完全一致。由表2的比对数据可以看出，各项信息的均值和标准值基本没有差异，GPS辅助导航结果与DGPS辅助导航结果一致，采用单点GPS或差分GPS对于导航信息和比力值的坐标分解没有影响。

3.3 基于SINS/GPS、SINS/DGPS组合导航的重力数据处理结果比对

分别对基于SINS/GPS和SINS/DGPS组合导航的重力测量数据，采用FIR低通滤波、零相移滤波、标准Kalman滤波、正反Kalman滤波方法对测区重力异常粗值进行滤波处理。其中，在重力数据处理阶段的各项改正中代入的速度、位置、高度数据均为差分GPS数据。处理结果分别如图8、图9所示，蓝色曲线为FIR低通滤波结果，红色曲线为标准Kalman滤波结果，绿色曲线为正反Kalman滤波结果，黑色曲线为零相位滤波结果。

图8 基于SINS/GPS的4种滤波方法处理重力异常结果比对Fig.8 Comparison of four filtering methods for gravity anomaly results based on SINS/GPS

图9 基于SINS/DGPS的4种滤波方法处理重力异常结果比对Fig.9 Comparison of four filtering methods for gravity anomaly results based on SINS/DGPS

对8条东西重复测线上的重力异常值进行滤波处理，结果比对如图10、图11所示。蓝色曲线为FIR低通滤波结果，红色曲线为标准Kalman滤波结果，绿色曲线为正反Kalman滤波结果，黑色曲线为零相移滤波结果。

为了更好地比较4种滤波方法及不同组合导航条件下的重力异常数据处理结果，表3列出了这8条东西重复测线在水平调整后的内符合精度。其中，最后1列为所有测线的内符合总精度。

图10 基于SINS/GPS的8条东西重复测线的4种滤波方法比对Fig.10 Comparison of four filtering methods on eight kinds of East-West repeated survey lines based on SINS/GPS

图11 基于SINS/DGPS的8条东西重复测线的4种滤波方法比对Fig.11 Comparison of four filtering methods on eight kinds of East-West repeated survey lines based on SINS/DGPS

表3 重复测线内符合精度评估Table 3 Accuracy evaluation in repeated survey lines

在表3中，以每条重复测线的内符合精度及所有重复测线的总内符合精度进行评判，分别对采用单点GPS进行组合导航后的相关重力测量数据处理结果及采用差分GPS进行组合导航后的相关重力测量数据处理结果进行比对，并列出了采用4种滤波方法所得的重复测线内符合精度，均为水平调整后的结果。

每条测线的均方根误差及标准差计算方法如下

式中，δij为第j条重复公共段各点观测值与该点各重复线观测值的平均值之差，即

m为重复线数目，n为重复线公共段数据点数。

所有重复线测试数据的总内符合精度计算公式为

从表3基于SINS/GPS组合导航的重力数据处理结果的最后一列，即所有重复测线的内符合总精度可以看出：2种Kalman滤波方式相比，正反Kalman滤波的均方根误差最小，为1mGal左右，标准Kalman滤波的均方根误差最大，达到近1.38mGal，两者相差0.38mGal；2种低通滤波相比较为接近，FIR低通滤波的均方根误差为1.19mGal，零相移滤波的精度为1.08mGal，略优于FIR低通滤波；而零相移滤波与正反Kalman滤波的结果相近，后者精度略高。

从表3基于SINS/DGPS组合导航的重力数据处理结果的最后1列数据可以看出：正反Kalman滤波效果最好，精度达到了0.98mGal；零相移滤波其次，达到了1.07mGal；FIR低通滤波的精度相对较差，为1.2mGal；标准Kalman滤波的精度结果为1.39mGal。

将两者比对可以看出：采用差分GPS数据进行组合导航处理之后得到的重力测量数据，相比于采用单点GPS数据进行组合导航得到的重力测量数据，精度基本接近，相差约为0.01mGal。

从以上分析可以看出：正反Kalman滤波方法的处理精度最高，优于2种低通滤波和单向标准Kalman滤波；正反Kalman滤波得到的重力异常曲线与另外3种曲线相比最为平滑，对高频噪声的去噪效果最好；零相移滤波与FIR低通滤波相比，其优点在于可进行正反2次低通滤波，对高频噪声的去噪效果更好，并抵消了正反相位移动的偏差，相对而言结果更为准确。同时，采用SINS/DGPS组合导航得到的重力测量数据，相比于采用SINS/GPS组合导航得到的重力测量数据，处理结果基本一致，这说明采用单点GPS或差分GPS进行组合导航，对导航解算结果和重力信息数据的处理结果并不影响。

4 结论

本文对FIR低通滤波、零相移滤波、标准Kalman滤波、正反Kalman滤波4种重力异常提取方法进行了比对，并采用以上方法分别对基于SINS/GPS和SINS/DGPS组合导航的实际SAG飞行试验数据进行了处理。重复测线的内符合精度统计结果表明，正反Kalman滤波与零相移滤波的结果精度相对较高，分别达到了0.98mGal和1.07mGal，正反Kalman滤波略优于零相移滤波；标准Kalman滤波的精度相对较差，精度与正反Kalman滤波相差近0.4mGal；FIR低通滤波的精度适中，与零相移滤波相比略有差距。本文验证了这4种滤波方法的实用性和优劣性，以及重力异常模型和窗函数设计的合理性。同时，也验证了单点GPS辅助导航与DGPS辅助导航的结果相同，即采用何种组合导航方式不影响后续重力数据的处理精度。本次试验也从侧面反映了SAG重力仪的测量精度，其内符合精度可达0.98mGal。