航空重力垂直加速度改正滤波方法对比研究

郑 崴，王建玲，陈学锋

航空重力垂直加速度改正滤波方法对比研究

郑 崴，王建玲，陈学锋

（河南工学院 电子信息工程学院，河南 新乡 453003）

根据航空重力测量的基本原理，阐述了用于航空重力垂直加速度改正的频域滤波器和时域滤波器的设计方法，以此为依据设计了用于垂直加速度改正的低通滤波器和卡尔曼平滑滤波器，并利用滤波器对航空重力实测数据进行了处理试验。试验结果表明：两种滤波方法都能够实现垂直加速度改正，而时域滤波器处理精度较高，与解算参考标准的误差低于1mGal。

航空重力测量；垂直加速度改正；低通滤波；卡尔曼滤波平滑

0 引言

航空重力测量是以飞机为载体，测量近地空中重力场信息的动态测量技术。它能够在人力无法达到的区域完成高精度、大面积的重力测量任务[1-4]。在测量过程中，它综合运用了航空重力仪、全球定位系统GPS（Global Positioning System）、惯性导航系统INS（Inertial Navigation System）等多种技术设备。航空重力仪测得的数据是重力信息与多种扰动加速度的总和，为获得高精度的重力异常信号，需对测量数据进行相应的改正，以消除扰动加速度的影响。其中，一部分的扰动加速度（如厄特缶斯加速度、水平扰动加速度）可以利用较为严密的解析公式通过计算来消除。而飞机垂直运动产生的扰动加速度无法直接观测，也没有解析公式可以直接计算得出，因此，为了实现垂直加速度改正，需要采用滤波技术来消除干扰[5-7]。

基于所采用的滤波技术的差别，垂直加速度改正又可分为基于频域滤波器的处理方法和基于时域卡尔曼滤波器的处理方法。前者是先利用数字差分对GPS定位信息进行处理以获得载体的垂直运动加速度，再将其与航空重力仪测得的数据求差，最后通过低通滤波器获取高精度的重力异常[8-10]。该种处理方法得以实现的依据在于：只有在信号低频端很小的范围内，重力信号的能量会高于噪声干扰，因此通过低通滤波器限制信号带宽，降低数据中的高频噪声干扰，可获取满足精度要求的重力异常。目前，常用到的处理方法包括无限脉冲响应（Infinite Impulse Response, IIR）滤波法、有限冲击响应（Finite Impulse Response, FIR）滤波法、快速傅里叶变换滤波法、小波变换滤波法等[11-13]。

而基于时域滤波器的处理方法则指的是基于卡尔曼滤波的处理方法，它是通过利用GPS/INS及稳定平台的观测数据对重力测量数据进行修正，在建立系统方程的基础上，以最小方差准则结合观测数据对重力异常进行估计，进而实现垂直加速度改正。俄罗斯国立大学的Bolotin等人基于航空重力仪沿飞行测线的动力学方程，建立了基于卡尔曼滤波的垂直加速度改正处理模型，实验结果表明该模型可获得精度在毫伽级别的重力异常[14-16]。由于种种原因，Bolotin等人只在其相关文献资料里介绍了该种方法的设计思路，但对于状态空间模型构建及相应的滤波器设计细节并未提及。张贵宾教授的课题组通过摸索和钻研，构建了基于重力异常模型的信息提取方法，并设计了固定区间平滑器，进一步提高了信息提取的精度[17][18]。

本文在航空重力测量、重力异常提取基本原理的基础上，根据航空重力测量的基本数学模型，依据不同滤波器处理的特点，分别设计了用于垂直加速度改正的频域滤波器和时域卡尔曼平滑滤波器，通过实测数据验证对比两种滤波器所获取的重力异常的精度，并探讨它们的特点。

1 航空重力测量基本原理

航空重力测量基本原理是从航空重力仪测得的比力数据（重力加速度与飞机垂直加速度之和）中减去垂直扰动加速度，最终获取所需的重力异常。依据牛顿第二定律，航空重力仪传感器质点在垂直方向的方程可表示为[19]：

正常场及高度改正项的计算则有[19]：

为方便分析，对式（1）进行调整，得：

2 卡尔曼滤波平滑处理

卡尔曼滤波是指在时域内，依据构造的状态空间方程对所需的状态向量进行估计的处理方法。为利用卡尔曼平滑滤波实现垂直加速度改正，首先需要依据式（1）建立航空重力场测量的状态空间模型，它由系统方程和量测方程组成。

在航空重力测量中，将航空重力仪沿测线的实时高程作为卡尔曼滤波器的量测方程，有：

在测量时，航空重力仪会因为稳定平台的舒拉或阻尼振荡而偏离水平方向，可以利用测得的稳定平台偏离水平方向的失准角对其进行修正，则式（4）可写为：

航空重力测量中，重力异常可以看作是沿测线分布的随机过程，因此用随机统计模型对其进行描述，有[20]：

将式（6）和（7）联立，可得：

将式（8）作为卡尔曼滤波的系统方程，将式（5）作为卡尔曼滤波的量测方程，并用离散系统状态空间形式表示为：

3 FIR低通滤波处理

与卡尔曼滤波处理算法不同，基于低通滤波器的垂直加速度改正处理方法是一种较为传统的处理方法。大量的研究表明，具有严格线性相位特征的FIR滤波器在航空重力数据处理中具有较好的效果。在信号处理技术高度发展的今天，滤波器的设计已经不是复杂的问题，重要的是如何根据信号特性选择合适的滤波器参数，因此本部分着重讨论载体垂直加速度的确定方法和FIR低通滤波器参数的选取方法，有关FIR低通滤波器的设计方法可参见相关文献[8][22]。

为实现基于FIR低通滤波的垂直加速度改正，首先要确定载体的垂直运动速度，惯用方法是对GPS测定的速度序列求一次差分运算，而GPS测速方式主要有以下三种[22]：

（1）位置微分法，利用GPS提供的位置数据进行差分运算，从而获得载体的运动速度；

（2）利用多普勒原始观测值直接计算载体速度；

（3）利用由载波相位中心差分获得的多普勒观测值导出载体速度。

已有研究指出，三种测速方法的精度基本相当，而位置微分法和载波相位中心差分测速精度略优于多普勒原始测速的精度[23]。

本文利用GPS提供的垂直速度进行一次差分来计算载体的垂直加速度。由于重力异常信号处于0.01 Hz的低频段内，因此采用两点中心差分即可满足精度要求[22]。两点中心差分期输入与输出关系有：

其传递函数为：

求得载体垂直加速度后，由于利用式（4）计算所获得的结果中还包含着大量的噪声干扰，还需利用FIR低通滤波器对其进行处理后，才能够获得高精度的重力异常。利用窗函数法设计FIR低通滤波器时，需要确定三个参数：窗函数的类型、窗函数的阶数和截止频率f。具体来说，在选取窗函数类型时，需要根据原始重力异常的频域特点来选择满足要求的窗函数，在航空重力测量数据处理中，海明窗、汉宁窗、布雷克曼窗和凯泽窗函数的性能比较优越[24]；对于来说，为保证滤波器拥有较好的幅频特性，使其能够准确地分离噪声和信号，应选择较宽的滤波窗口，但这会导致处理后有效数据点的减少，反之如降低滤波窗口宽度可减弱边界效应的影响，但滤波器的幅频特性便不能保证；对于截止频率f的选择来说，要综合考虑重力异常精度和分辨率。航空重力测量中，低通滤波器截止频率与分辨率的关系有[22]：

4 实测数据处理试验及分析

利用航空重力仪在我国南部某省进行了近海重力测量，本试验所用数据为第五架次2号测线的测量数据，该测线是一条东西向的测线。测量选取在气候条件较好的夜间进行，以减小气流对测量结果的影响。测量过程中，设计飞行高度为400 m，飞行速度为60 m/s，实际飞行测量高度在设计高度的±5米范围内波动，图1为该条测线的飞行高度变化示意图。

分别利用卡尔曼平滑滤波器和FIR低通滤波器对测量数据进行了处理，其中在FIR低通滤波器设计中采用了海明窗，窗口阶数选取为600，截止频率选取在0.004 Hz。此外，由于暂时无法获得地面重力测量数据，试验过程中以商用软件的处理结果作为参考标准，将本试验的处理结果与之进行对比研究。

根据本文第2部分和第3部分的方法对该测线处理后，所获重力异常结果如图2所示。图中，横坐标为GPS测量时间；纵坐标为处理所得重力异常值；实线为商用软件处理所得结果，在本试验中作为参考标准；虚线是FIR低通滤波的重力异常曲线，为方便对比，受边界效应影响的测量点已经去除；点划线为卡尔曼平滑滤波处理所得重力异常曲线。

图1 飞行高度变化示意图

图2 FIR低通滤波与卡尔曼平滑滤波处理重力异常对比图

表1 FIR低通滤波与卡尔曼平滑滤波处理重力异常误差对比表

由图2可见，经过处理后，测量结果中的噪声已经被抑制，重力信号被保留了下来，获得了高信噪比的重力异常信号。对比两种处理方法的结果，FIR低通滤波的结果与参考标准的偏差较为明显。与此同时，还应注意到无论是FIR处理的结果还是卡尔曼滤波平滑处理的结果，重力异常曲线在测线末端均出现了较为强烈的波动，且与参考标准出现了明显的偏差，分析其原因可能是因为飞行过程中产生了较为明显的垂直方向的机动而引入了部分的噪声。

表1为两种滤波方法处理所得重力异常与参考标准的误差统计。由表可见，利用FIR低通滤波器与参考标准的均方差为1.2329 mGal，而卡尔曼平滑滤波的处理结果的均方差仅为0.7219 mGal，远远优于FIR低通滤波的处理效果，这表明相比较于传统的低通滤波处理方法，卡尔曼平滑滤波处理所获得的重力异常具有更高的精度。

5 结论

本文依据航空重力测量基本原理，研究了两种用于垂直加速度改正处理的滤波方法，并对它们的处理结果进行了对比试验。由试验结果得出如下结论：

（1）本文所提出的卡尔曼平滑滤波状态空间模型，经实测数据验证，可以获取高精度的重力异常信号，误差小于1 mGal；

（2）尽管基于低通滤波的处理方法所获得的重力异常误差高于卡尔曼平滑滤波处理的结果，但是其处理方法和步骤相对简单，也不失为一种在实际应用中可以选择的处理策略；

（3）无论采用何种滤波处理方法，测量过程中，飞行载体都应尽可能保持平稳的飞行状态，这有助于提高航空重力测量的精度。

[1] 夏哲仁,石磬,孙中苗.航空重力测量系统CHAGS[J].测绘学报,2004,33(3):216-220.

[2] 夏哲仁,孙中苗.航空重力测量及其技术应用[J].测绘科学,2006,31(6):43-46.

[3] 熊盛青.我国航空重磁勘探技术现状与发展趋势[J].地球物理学进展,2009,24(1):113-117.

[4] 熊盛青,周锡华,郭志宏，等.航空重力勘探理论方法及应用[M].北京:地质出版社,2010.

[5] 黄谟涛,宁津生,欧阳永忠，等.航空重力测量厄特弗斯改正公式注记[J].测绘学报,2015,44(1):6-12.

[6] 黄谟涛,宁津生,欧阳永忠,等.海空重力测量平台倾斜改正模型等价性证明与验证[J].武汉大学学报（信息科学版）,2016,41(6):738-744.

[7] 孙中苗, 石磬, 夏哲仁，等. 利用GPS和数字滤波技术确定航空重力中的垂直加速度[J]. 测绘学报, 2004, 33(2): 110-115.

[8] HAMMADA Y. A comparison of filtering techniques for airborne gravimetry[D]. Calgary: University of Calgary, 1996.

[9] 徐泽扬. 航空重力数据滤波处理方法研究[D].南京:东南大学,2013.

[10] 蔡劭琨. 航空重力矢量测量及误差分离方法研究[D].长沙:国防科技大学,2014.

[11] 郭志宏, 罗锋, 王明，等. 航空重力数据无限脉冲响应低通数字滤波器设计与实验研究[J]. 地球物理学报, 2011, 54(8): 2148-2153.

[12] 李宜龙, 曾敏, 孙科，等. 航空重力测量中FIR低通滤波与高斯滤波的比较[J]. 海洋测绘, 2014, 34(6):40-42.

[13] 郎骏健，梁星辉，柳林涛，等．航空重力傅里叶基追踪低通滤波方法研究[J]．地球物理学报, 2018, 61(12):4737-4745.

[14] BOLOTIN Y V, POPELENSKY M Y. The analysis of accuracy of airborne gravimetry with stochastic models[J]. Aerokosmichekoe priborostroenie, 2003, 4: 42-48.

[15] BOLOTIN Y , DOROSHIN D . Adaptive filtering in airborne gravimetry with hidden Markov chains.[J]. IFAC proceedings volumes,2011,44(1):9996-10001.

[16] BOLOTIN Y V. GOLVAN A A. Methods of inertial gravimetry[J]. Moscow university mechanics bulletin,2013,68(5):117-125.

[17] 郑崴，张贵宾.自适应卡尔曼滤波在航空重力异常解算的应用研究[J].地球物理学报,2016,59(4):275-1283.

[18] 郑崴,张贵宾,陈涛,等.Sage自适应滤波在航空重力数据处理的应用[J].物探与化探,2015,39(S1):84 -90．

[19] BOLOTIN Y V. Mathematics behind GTGRAV[R]. Moscow: Moscow Lomonosov State University,2009.

[20] 蔡体菁,周薇,鞠玲玲.平台式重力仪测量数据的卡尔曼滤波处理[J].中国惯性技术学报,2015,23(6):718-723.

[21] 秦永元,张洪钺,汪叔华.卡尔曼滤波与组合导航原理（第二版）[M]．西安:西北工业大学出版社,2012．

[22] 孙中苗.航空重力测量理论、方法及应用研究[D].郑州:解放军信息工程大学, 2004.

[23] 何海波,杨元喜,孙中苗.几种GPS测速方法的比较[J].测绘学报,2002,31(3):217-221.

[24] 郭志宏,罗锋,安战锋.航空重力数据窗函数法FIR低通数字滤波试验[J]．物探与化探,2007,31(6):568-571.

Comparing Research for Vertical Acceleration Correction Filtering Methods in Airborne Gravimetry

ZHENG Wei, WANG Jian-ling, CHEN Xue-feng

(School of Electronic Information Engineering, Henan Institute of Technology, Xinxiang 453003, China)

Based on the principle of airborne gravimetry, the vertical acceleration correction methods are demonstrated in this paper. Then, the methods based on low-pass filter and Kalman smoother are designed for the correction respectively. According to the diversities of vertical acceleration correction methods, the filter and the smoother are applied in measurement data processing which is from practical airborne gravimetry. The result of the processing shows that both of the filter and smoother can process the measured data to obtain the gravity anomaly and the errors is about 1 mGal. Specifically, the error based on Kalman smoother is less than 1 mGal.

airborne gravimetry; vertical acceleration correction; low-pass filtering; Kalman smoothing

V241

A

2096–7772(2020)01–0028–05

2019-12-14

郑崴（1984―），男，河南正阳人，讲师，博士，主要从事航空重力测量技术方法、复杂信号处理与提取技术研究。

（责任编辑吕春红）