闫 方，胡平华，赵 明，姜作喜，罗 锋，唐江河，刘东斌

(1.北京自动化控制设备研究所，北京 100074；2.中国国土资源航空物探遥感中心，北京 100083)

0 引言

地球重力是地球引力和地球自转引起的惯性力的合力，是地球的基本物理特征之一。精确的局部地球重力场信息对自然资源勘察、地质科学研究、自然灾害监测及预警、地球极区勘探、高精度战略武器系统研制保障具有重要的意义[1-3]。

航空重力测量是一种以航空器为载体，综合运用重力敏感器、高精度姿态稳定系统、定位系统，获取近地面或近海重力场的方法[4]。

结合重力测量的发展历程，当前还有走行式地面重力测量、海洋船载重力测量和卫星测高重力测量等方法[5]。其中走行式方法借助静态重力仪在地面进行，这种测量方式耗时长、效率低、人工成本高，同时受山川、湖泊、沼泽等自然条件的制约；船载方式更适合远洋应用场景，在近海岸区域由于水流波动较大，容易造成测量精度下降；借助卫星测高资料反演得到的重力场只能反映地球重力场的中低频信息，限制了卫星重力资料的应用领域[6]。而采用航空重力仪系统的机载测量方式，可以克服上述测量手段中的缺点，高效快速地获取较大区域内的重力资料，同时能够保证测量的精度与分辨率[7-8]。

对于航空重力数据处理，目前主要有频域有限脉冲响应(Finite Impulse Response，FIR)低通滤波、时域Kalman滤波以及时频域小波分解去噪方法[9]。其中Kalman滤波方法对于重力异常模型的精度要求较高，小波方法在重力测量领域缺乏明确的理论依据，因此实际应用都受到了限制。而频域FIR低通滤波方法物理含义清晰，实际应用广泛。

除了上述常规手段，国外还有具有特色的滤波方法。Von Frese等提出了波数相关滤波器(Wavenumber Correlation Filter，WCF)，该方法通过比较测量数据与协同数据频谱的相关程度来实现噪声分离[10]。Alberts提出了频域加权的方法，对测量数据的不同频段赋予不同权值，以此进行降噪处理[11]。Bolotin在时域处理方法的基础上，将地球重力场的非均匀一致性建模为多状态隐性马尔可夫模型，对结果的细节有所提升[12]。

当前，国内对采用加速度计式重力敏感器的三轴稳定平台重力仪数据后处理方法研究较少。东南大学蔡体菁带领的研究团队采用Kalman滤波的方法对平台重力仪实测数据进行了处理，得到了较好的结果，此外该团队还提出了巴特沃斯与Kalman平滑滤波级联的处理方法。

本文采用频域信息处理的方法，结合平台重力仪的特点，系统地给出了一套在频域内有严格规范和明确含义的航空重力数据后处理方法。

1 航空重力测量的基本原理与数学模型

航空重力测量围绕2个基本问题展开，一个是如何在运动的环境下稳定重力敏感器的指向，使其严格保持垂向；另一个是如何从重力敏感器输出的惯性加速度中分离得到重力加速度[13]。以下将结合平台式航空重力仪的系统组成和重力测量的基本数学模型进行简述。

1.1 平台式航空重力仪系统组成

平台式重力仪系统硬件部分主要由重力仪主机、差分全球定位系统(Differential Global Positioning System，DGPS)、不间断电源、显控存储装置和减振器组成。

其中重力仪主机的核心装置为惯性稳定平台，主要由1个高精度石英挠性加速度计式重力敏感器、2个动力调谐陀螺、2个导航级石英挠性加速度计和3个环架组成。该惯性稳定平台在水平面内稳定精度优于10″，可以隔离载体角运动对重力敏感器的影响，同时可以跟踪当地地理坐标系，严格保持重力敏感器指向垂向[14]。

差分GPS是航空重力测量系统中另一关键部分，由安装在载体上的移动站和固定在地面的基站组成。原始差分卫星测量数据经过解算可以得到载体高精度的位置与速度信息，用来计算载体垂向方向的运动加速度以及重力相关的改正项[15]。

1.2 航空重力测量的数学原理

航空重力测量的数学模型是由惯性导航原理中比力方程推导而来的。本文研究的是标量重力仪数据后处理，以下给出标量形式重力测量的基本方程

(1)

将式(1)中后三项单独记为δaE，即为厄特弗斯改正项

(2)

由于重力测量一般都是为了获取被测区域的重力异常值，因此从重力测量值中减去标准重力数值，得到以下计算重力异常的公式

(3)

其中，g0为正常重力场，根据不同的需求采用相应的正常重力模型进行计算。常用的正常重力公式主要有赫尔默特公式和1985年国际正常重力公式。需要注意的是，这些公式计算值为测点在地球旋转椭球体表面处的正常重力值，而高度每增加1m，重力值就降低约0.3086mGal，因此在航空重力数据处理中使用的正常重力场公式需包含高度校正部分。

以下给出高度校正部分的公式[16]：

1)当飞行海拔高度h≤700m

δgh=0.3086h

(4)

2)当飞行海拔高度h>700m

δgh=0.3086h-0.073h2

(5)

经过上述流程的计算，即可实现重力敏感器输出中惯性加速度和重力加速度的分离，同时还对重力异常信号进行了必要的改正。在实际航空重力测量中，重力敏感器的输出还包含由飞机发动机以及空气湍流造成的随机振动加速度，这部分噪声的强度是重力异常的数万倍。因此，需要利用滤波的方法进一步去除重力异常中的高频噪声[17]。

2 航空重力数据后处理流程与关键技术

根据上述重力测量的基本原理，结合平台式重力仪的结构特点，给出如图1所示的航空重力数据后处理流程。

图1 航空重力数据处理流程图Fig.1 Flow chart of gravity measurement data processing

结合上述流程图，对于平台式航空重力仪的数据处理，有以下几个问题需要注意：

1)重力敏感器和差分GPS的输出频率分别为100Hz和2Hz，且2个通道独立工作，因此需要实现通道之间的采样率统一；

2)在对各个通道内原始测量数据进行处理的过程中，还需考虑每一步操作对信号相位的影响，从而避免在作差环节因相位不同步而引入误差；

3)在对计算得到的重力异常原始值进行滤波时，同样需要考虑信号的相位变化，保证解算得到的重力异常值可以归算到正确的地理位置。

采用低通FIR滤波方法的前提是认为重力异常信号具有低频特性，因此数据处理的每个环节都需要有明确的频域含义，避免计算过程对原始测量信号的低频信息产生破坏，进而保证数据处理的精度。

针对上述问题，提出了以下几点数据处理的关键技术。

2.1 采用降采样率方法实现通道间采样率统一

多采样率数字信号处理方法有着严格的频域理论基础，在满足一定条件的情况下可以提升或者降低数据的采样率。以下简述该方法的基本思想。

首先介绍降低采样率的方法。现有原始信号序列x(n)，信号中最高频率为Fc，原始采样频率为Fs。降低采样率M倍最直接的方法是在序列中每M个点抽取一个点，组成一个新的序列x′(n)。根据奈奎斯特采样定律，新序列不发生频谱混叠的必要条件是新的采样率大于信号内最高频率的2倍，即

Fs/M≥2Fc

(6)

由于M是可变的，为了保证式(6)始终成立，可以首先对原始信号进行低通滤波以实现抗混叠，然后再进行抽取。其流程图如图2所示。

图2 降采样流程图Fig.2 Flow chart of downsampling

对于提升L倍采样率，首先在每2个相邻点之间插入(L-1)个零值，然后通过一个低通滤波器，将插入的零值点计算得到相应的抽样值。该方法在理论上的依据可以参考相关文献[18]。其实现的流程图如图3所示。

图3 升采样流程图Fig.3 Flow chart of upsampling

综合使用上述方法，可以将通道间采样统一到相同频率。本文采用的方案是将重力仪的采样率直接降低到2Hz，为此设计了截止频率为1Hz的FIR低通滤波器，其幅频特性如图4所示。

图4 FIR低通滤波器幅频响应图(1Hz)Fig.4 Amplitude frequency response curve ofFIR low pass filter (1Hz)

2.2 采用中心差分器求解载体运动加速度

差分GPS可以输出高精度定位信息，对其中的高度序列做二次差分，可以计算得到载体运动加速度。

理想差分器的频率响应为

Hd(ejω)=jω

(7)

由式(7)可知，理想差分器的幅频特性为线性增长。

工程应用中，通常采用在低频范围内线性度较好的牛顿中心差分器[19-20]。将差分GPS输出的高度序列记为x(n)，则牛顿中心差分计算过程为

(8)

中心差分运算可以通过FIR滤波器形式实现，其幅频特性曲线如图5所示，相频特性曲线如图6所示。

图5 差分器幅频特性图Fig.5 Amplitude frequency response curve of differential filter

图6 差分器相频特性图Fig.6 Phase frequency response curve of differential filter

从幅频特性图中可以看出，中心差分器在低频段内线性度较好。

由相频特性图可以看出，其相位具有线性递减特性。在数字信号处理中，群延迟定义为相位变化与频率变化的比值，即：

TGroup=-dθ/dω

(9)

群延迟反映了信号中所有频率成分的时间延迟，对于线性相位的系统，群延迟是常值。借助这一方法，可以精确计算得到中心差分器造成的时延，并在计算过程中予以补偿。

2.3 采用零相位滤波器保证重力异常的相位精度

在数据处理的过程中，采用FIR低通滤波器会产生相位延迟，导致输出信号相对于输入信号在时间上发生平移。体现在重力数据的处理结果上，会使得滤波解算得到的重力异常位置偏离真实地理位置。因此在设计低通滤波器的过程中，需要采用零相位滤波进行处理[21]。

本文采用Forward-Backward零相位滤波器，其实现方法是：首先将输入信号x(n)按正向顺序通过滤波器，即与数字滤波器冲激响应序列h(n)进行卷积运算

u(n)=x(n)*h(n)

(10)

将得到的结果u(n)进行时间反转为v(n)

v(n)=u(N-1-n)

(11)

再次通过滤波器得到序列w(n)

w(n)=v(n)*h(n)

(12)

最后将这一序列再次进行时间反转，即得到精确的零相位结果y(n)

y(n)=w(N-1-n)

(13)

为了进一步说明该方法的零相位特性，以下给出这一过程的频域描述，其中滤波器的频域表示为H(ejw)，输入信号正向通过滤波器后有

U(ejω)=X(ejω)H(ejω)

(14)

然后进行时间反转

V(ejω)=e-jω(N-1)U(e-jω)

(15)

再次通过滤波器

W(ejω)=V(ejω)H(ejω)

(16)

最后再次将所得结果进行时间反转，有

Y(ejω)=e-jω(N-1)W(e-jω)=X(ejω)|H(ejω)|2

(17)

由式(17)可以看出，输入与输出之间没有相位变化，可以实现精确的零相位操作。图7所示为零相位滤波流程图。

图7 零相位滤波流程图Fig.7 Flow chart of zero-phase filter

3 航空重力实测数据处理结果与精度评估

3.1 试验概况

为了验证平台式重力仪的实际工作性能，中国国土资源航空物探遥感中心在海南岛附近某海域组织完成了机载航空重力测量试验。本次试验采用Cessna208b固定翼飞机，该机型配有自动驾驶仪，重力测线海拔高为600m，飞行速度为60m/s。

飞行试验共进行了4个架次，其中东西方向在同一测线安排有3个架次，每架次均取得4条有效重复测线；在南北方向有1个架次共4条重复测线。此外，试验还在东西测线上进行了1次起伏飞行。为更好地评价平台式重力仪测量效果，试验提供了同测线GT-2A的事先测量结果。

3.2 航空重力测量结果精度评估方法

航空重力测量数据精度评价指标主要有内符合精度和外符合精度。其中内符合精度是计算同一台仪器在重复测线测量值的均方根(Root-Mean-Square，RMS)。外符合精度是以其他手段获取的重力异常值作为标准值来统计的标准差(Standard Deviation，STD)。采用文献[22]中给出的计算方法，分别给出滤波周期为140s和100s的统计结果。

3.3 滤波周期140s的重复测线结果

采用滤波周期为140s，即截止频率为0.00714Hz的FIR低通滤波器，可以得到4个架次试验的测量结果。结果中重力异常的地理分辨率为4.2km。

4个架次的重复测线处理结果及其对应的GT-2A测量结果分别如图8～图11所示。

图10 东西方向第3架次测线及GT-2A结果(140s)Fig.10 Third east-west repeated measurementlines and the result of GT-2A(140s)

图11 南北方向架次测线及GT-2A结果(140s)Fig.11 South-north repeated measurementlines and the result of GT-2A(140s)

上述处理曲线的内符合和外符合精度结果统计如表1所示。

将3个架次东西方向12条重复测线的结果综合到一起进行分析，并与其对应的GT-2A测量结果进行对比，其结果如图12所示。

表1 四架次重复测线内、外符合精度(140s)

3.4 滤波周期100s的重复测线结果

采用滤波周期为100s，即截止频率为0.01Hz的FIR低通滤波器，可以得到4个架次试验的测量结果。结果中重力异常的地理分辨率为3km。

4个架次的重复测线处理结果及其对应的GT-2A测量结果分别如图13～图16所示。

图13 东西方向第1架次测线及GT-2A结果(100s)Fig.13 First east-west repeated measurementlines and the result of GT-2A(100s)

图14 东西方向第2架次测线及GT-2A结果(100s)Fig.14 Second east-west repeated measurementlines and the result of GT-2A(100s)

图15 东西方向第3架次测线及GT-2A结果(100s)Fig.15 Third east-west repeated measurementlines and the result of GT-2A(100s)

图16 南北方向架次测线及GT-2A结果(100s)Fig.16 South-north repeated measurementlines and the result of GT-2A(100s)

将3个架次东西方向12条重复测线的结果综合到一起进行分析，并与其对应的GT-2A测量结果进行对比，其结果如图17所示。

图17 东西方向所有测线及GT-2A结果(100s)Fig.17 All east-west repeated measurementlines and the result of GT-2A(100s)

上述处理曲线的内符合和外符合精度结果统计如表2所示。

表2 四架次重复测线内、外符合情况(100s)

3.5 起伏飞行条件下的测试结果

本次试验在东西方向第2架次飞行中进行了2次起伏飞行，最大起伏高度为100m，垂向速度为2m/s左右。起伏飞行的高度与垂向速度变化曲线如图18所示。

图18 起伏飞行高度和垂向速度曲线图Fig.18 Figure of height and vertical velocity underrise-and-fall flight

分别给出140s和100s滤波周期下，起伏飞行测得的重力异常与GT-2A平稳飞行测量结果及同架次其他4条重复线的均值结果的对比图，分别如图19和图20所示。

图19 起伏飞行测量结果与重复线均值及GT-2A结果比较(140s)Fig.19 Comparison of result under rise-and-fall flight conditionwith mean of repeated lines and the result of GT-2A(140s)

图20 起伏飞行测量结果与重复线均值及GT-2A结果比较(100s)Fig.20 Comparison of result under rise-and-fall flight conditionwith mean of repeated lines and the result of GT-2A(100s)

在140s滤波周期下，起伏测量结果与GT-2A的标准差为0.81mGal，与同滤波尺度下重复测线均值的标准差为0.60mGal。

在100s滤波周期下，起伏测量结果与GT-2A的标准差为1.11mGal，与同滤波尺度下重复测线均值的标准差为1.03mGal。

上述指标及具体曲线表明，平台式重力仪具备一定的起伏飞行能力。

4 结论

本文针对平台式重力仪航空测量数据后处理的要求，结合平台式重力仪的系统原理和结构特点，设计了一套数据后处理的流程，并对其中的关键处理环节进行了频域分析。通过对飞行试验实测数据进行处理，得到如下结论：

1)试验中4个架次的重复线内符合精度达到：12条东西测线0.43mGal/140s和0.84mGal/100s，4条南北测线0.39mGal/140s和0.79mGal/100s，表明平台式重力仪性能稳定可靠；

2)以GT-2A单次测量结果(其本身也存在一定误差)为标准值统计外符合精度，12条东西测线为0.72mGal/140s和0.98mGal/100s，4条南北测线为1.41mGal/140s和1.53mGal/100s，东西测线精度较好，南北测线精度较差。经分析，原因可能是GT-2A的此次南北测线测量结果存在较大误差，后续将针对此问题进行重点验证；

3)在起伏高度为100m，垂向速度为2m/s的条件下，重力异常的内符合及外符合精度没有显著变大，表明平台式重力仪具备一定的起伏飞行能力；

4)通过对以上实测数据的处理，验证了本文设计的航空重力数据后处理方法的正确性。