受海洋重力测量的启发，针对人类难以到达的区域，科学家们将重力仪安装在飞行器上进行重力测量即航空重力测量.20世纪50年代，美国率先开展了航空重力测量实验，所用仪器为L&R型海洋重力仪，仪器搭载于KC-135喷气式飞机，采用摄影测量与多普勒进行定位.利用当时已有的导航定位系统并结合光学、雷达距离跟踪等地面控制技术，航空重力测量的精度仅能达到10 mGal(1 mGal=10^-5 m·s^-2)(Thompson and LaCoste, 1960；Nettleton et al., 1960).直到20世纪90时代，随着GPS技术特别是载波相位测量和差分GPS定位取得突破后，航空重力测量才开始得到广泛应用.目前，GPS的定位、测速精度已经能够满足航空重力测量1 mGal的精度要求(Brozena et al., 1989；Klingelé et al., 1997；孙中苗, 2004).

航空重力测量平台处于不断高速运动状态下，重力传感器输出的总加速度中，除了含有重力信号外，还包含因载体运动、动态环境等因素引起的扰动加速度，其中一部分能用解析式表示，如厄特弗斯效应，另一部分是与载体非匀速运动有关的非规则影响，如垂直和水平扰动加速度等(孙中苗，2004；欧阳永忠，2013).需要特别指出的是，航空重力测量中重力异常信号往往不足100 mGal，而扰动加速度和观测噪声的幅值可能比重力信号高出几百倍甚至几千倍.重力信号在低频段的很窄部分比噪声强而在高频部分则完全被噪声湮没，重力信号比噪声强的低频段范围没法精确确定，需要对重力信号以及各项改正做频谱分析，以此来确定低频段的大致范围，获取航空重力测量的频谱窗口(孙中苗，2004)，为了消除或者减弱观测噪声的影响以获取真实的重力信号，需要将原始重力扰动经过低通滤波处理(孙中苗和翟振和，2009；梁星辉，2012；梁星辉等，2013).

低通滤波是航空重力测量数据处理过程中的重要环节，也是直接影响航空重力测量成果质量的关键技术之一.在航空重力数据处理中，目前已经普遍使用FIR低通滤波器(Schwarz and Li, 2000；Bruton, 2001；欧阳永忠，2013)，第二代L&R海空重力仪也已经使用FIR滤波器取代传统的RC滤波器.除此之外，国内外学者还提出了一些数据滤波方法.采用巴特沃斯、切比雪夫逼近设计的IIR(infinite impulse response)低通滤波器，在一些航空重力测量实验中得到了应用(Forsberg et al., 2001；Meyer et al., 2003；郭志宏等，2011)；Kwon(2000)将波数相关滤波器(wavenumber correlation filter，WCF)应用于航空重力测量数据处理；柳林涛和许厚泽(2004)利用连续小波滤波方法对航空重力测量数据进行滤波处理；孙中苗和翟振和(2009)采用实际的航空重力测量数据，从小波选择、阈值施加方式及分解层次确定等方面讨论了小波阈值滤波方法的适用性；张开东(2007)联合使用固定区间的Kalman滤波平滑算法和迭代算法，对捷联惯导航空重力测量数据进行滤波估算；Alberts(2009)提出一种基于频率数据加权处理的新方法，在能够获得比较可靠的噪声模型条件下，使用分频段加权策略.

采用窗函数法实现的FIR滤波器在航空重力测量数据处理中已得到普遍应用，但由于FIR滤波器的长度一般较长，实际处理时以窗口中点作为输出计算点，则一边损失半个滤波窗口数据，这些数据由于存在边缘效应，因而是不准确的，实际应用时应该去掉(郭志宏等，2007).假设飞机速度100 m·s^-1，数据采样率2 Hz，截止频率取0.01 Hz，FIR窗函数法滤波器的长度根据主瓣宽度的一半等于截止频率来估算(孙中苗，2004)，采用汉宁窗，滤波器的长度需要取400，由此估算测线两端都要舍弃半个滤波窗口10 km的数据，对于成本昂贵的航空重力测量，是一个不小的损失.

针对航空重力数据滤波处理中的边缘效应问题，本文引进了一种新的滤波方法——傅里叶基追踪低通滤波方法，通过该方法得到的低通滤波器能有效减少有限时间序列造成的谱污染和谱泄漏，很好地减弱边缘效应.根据傅里叶基追踪低通滤波方法设计傅里叶基追踪低通滤波器，傅里叶基追踪低通滤波器的英文全称为Fourier basis pursuit low pass filter，选其每个单词的首字母并大写作为英文缩写，则英文缩写为FBPLPF.本文先介绍FIR低通滤波器和傅里叶基追踪低通滤波器的原理，然后基于仿真数据和实测数据验证新方法的有效性，将新方法与现在普遍使用的FIR低通滤波器以及常规快速傅里叶(fast Fourier transformation，FFT)滤波方法作比较，分析了新方法的性能.

IIR系统既有零点又有极点而FIR系统只有零点(胡广书，2012)，因此FIR系统不像IIR系统那样易取得比较好的通带与阻带衰减特性，要想取得好的衰减特性，需要滤波器的长度较长.但FIR系统有突出的优点，其一是系统总是稳定的，其二是易实现线性相位(胡广书，2012).航空重力测量不可避免的含有各种粗差，此时递归运算中的粗差影响比非递归要大，这时候系统的稳定就显得比较重要，此外，还需要所设计的滤波器具有精确的线性相位，因此FIR滤波器在航空重力测量数据处理中得到了广泛的应用(孙中苗，2004).目前，FIR数字滤波器的设计方法主要是建立在对理想滤波器频率特性作某种近似的基础上的.近似方法有窗函数法、频率抽样法及最佳一致逼近法，本文采用航空重力数据滤波中常用的窗函数法.

假设理想低通数字滤波器的截止频率为w_c，且具有线性相位，群延迟为α，该滤波器的单位冲激响应为

(1)

式(1)表示的单位冲激响应是无限长的，而实际中要求滤波器为有限长.窗函数法的设计思想是选取某种效果好的窗函数对h_d(n)予以截短.实际工程应用中常用的窗函数有矩形窗、三角窗、汉宁窗、哈明窗、布莱克曼窗和凯塞窗.根据孙中苗(2004)、郭志宏(2007)、欧阳永忠(2013)等学者的研究，汉宁窗和哈明窗的效果最好，故本文采用汉宁窗.滤波器长度的确定需要顾及两方面的因素，即长度的奇偶性和单位抽样响应的绝对长度.偶数长度的单位冲激响应比同等的奇数长度具有更好的幅频响应，但会产生非整数的时间延迟，滤波后需将非整数历元内插至整数历元，考虑到内插误差可能较大，宜采用奇数长度(孙中苗，2004).FIR窗函数法低通滤波器的长度根据主瓣宽度的一半等于截止频率来估算.

通过傅里叶变换实现信号低通滤波是目前常用的滤波方法，然而由于数据长度有限，离散傅里叶变换谱存在严重的谱泄露和谱污染，使得滤波得到的信号存在严重的边缘效应，常用方法是采用加窗傅里叶变换予以抑制，但还是难以有效克服由于谱泄露产生的边缘效应.信号的稀疏表达在最近这些年得到了广泛的关注和研究，在很多方面有应用，比如信号分离、去噪、编码、图像修复及提取极移信号(Olshausen et al., 2001; Starck et al., 2005; Li et al., 2004; Elad and Aharon, 2006; Wang et al., 2016)，本文将其应用于航空重力数据低通滤波.

信号稀疏表达的原理是借助变换域中原子的特征，用尽量少的基函数准确地描述原始信号，从而抓住信号的本质，以字典中与信号匹配的原子的特征反应信号的特征.稀疏表达主要涉及稀疏表达算法及原子字典两方面，学者们已经提出了许多稀疏表达算法，常用的有匹配追踪算法(Mallat and Zhang, 1993)、基追踪算法(Chen et al., 2001)、框架方法(Daubechies, 1988)、最佳正交基算法(Coifman and Wickerhauser, 1992)、正交匹配追踪算法(Pati et al., 1993)；针对信号的不同特征，学者们提出了许多原子字典，如傅里叶字典，小波包字典，余弦包字典(Coifman and Meyer, 1991)，离散余弦变换字典，Chirplet字典(Mihovilovic and Bracewell, 1991；Mann and Haykin, 1992)、Warplets字典(Baraniuk and Jones, 1993)等.Mallat(1999)和Chen等(2001)证明，基追踪算法得到的结果最为稀疏，而傅里叶字典是最为简单、最为基本的字典，基于这两点本文提出了基于基追踪的傅里叶低通滤波方法.

傅里叶基追踪低通滤波的基本原理是基于傅里叶原子库，通过基追踪准则，将低频信号挤压在低频上，实现低频信号和高频信号的有效分离.通过该方法得到的低通滤波器能有效减少有限时间序列造成的谱污染和谱泄漏，因此该方法能有效抑制边缘效应，具体算法如下.

离散傅里叶变换逆变换的公式：

(2)

以式(2)为基础建立矩阵Φ，令J=[0, 1, 2, …, M-1], 则

(3)

式(3)中M表示傅里叶字典中傅里叶原子数目，N表示时间序列的长度，矩阵Φ的每一列代表傅里叶字典中的傅里叶原子，傅里叶字典中的每个傅里叶原子代表着不同的频率特征.将式(2)写成矩阵形式：Φx=f，其中f为原信号，x为每个傅里叶原子对应的系数.为了使字典超级完备，要求M>N，即傅里叶原子数大于时间序列的长度(Chen et al., 2001)，此时原方程秩亏，求解需要找到一种准则.由于基追踪算法选择全局优化，因此能获得比匹配追踪算法更高的分辨率，且能得到比框架方法更稀疏的结果(Chen et al., 2001)，故采用该准则，准则具体为：在l₁范数下，满足方程组解的模和最小(Mallat, 1999)，即

(4)

式(4)的求解是有关凸优化的问题，要将其转化为线性规划问题求解(Mallat, 1999).在过去的50多年里，研究者们对线性规划问题做了很多的研究工作，其中最重大的突破是发明了内点算法(Chen and Donoho, 1994)，本文在求解时采用内点算法.

在l₁范数下采用内点算法求得式(4)的解x^*，x^*表示傅里叶基追踪谱.若x^*取频率范围[0, w₁]，则可以得到对应的低通滤波结果为

(5)

式(5)中，z即为频率范围在[0, w₁]的滤波值.

选取秦岭地区GT-2A航空重力仪实测的某条测线重力扰动数据，该数据是由仪器自带商用软件GTGRAV经过一系列复杂的数据处理过程得到的，数据如图 1所示.从现有的资料可得知GTGRAV软件同时采用了FIR和Kalman滤波，可选择不同的滤波长度，其他的参数如安装误差，重力仪尺度因子，输出时间延长可以根据Kalman滤波估计，详细的算法过程和其他参数尚未公开发表，本文选取的是滤波长度为100 s，即截止频率为0.01 Hz的结果.航空重力测量中，高频部分的噪声远远大于低频部分的噪声，因此在实测数据中加入调相白噪声，该噪声的幅值随着频率的增大而增大，加入噪声得到模拟的原始测线重力扰动数据如图 2所示.图 1、图 2的横坐标表示测线上测量的点数，点与点之间的距离大约30 m.从图 1和图 2可以看出，重力扰动信号幅值不到100 mGal，加入的高频噪声的幅值达到5000 mGal，重力扰动信号完全淹没于幅值很大的高频噪声中，这与实际航空重力测量的噪声情况类似.重力信号频谱和仿真噪声频谱如图 3所示，横坐标表示频率，纵坐标表示幅值，从图中可以看出，重力信号在低于0.01 Hz频段时比噪声要强，而在高于0.01 Hz的频段噪声远远大于重力信号.模拟噪声的特性与实际情况比较符合(梁星辉，2012)，只有在甚低频部分，重力信号强于噪声.

图 1

Fig. 1

图 1 测线空间重力扰动数据 Fig. 1 The data of survey line gravity disturbance

图 2

Fig. 2

图 2 加入噪声的空间重力扰动数据 Fig. 2 The data of survey line gravity disturbance with noise

图 3

Fig. 3

图 3 重力信号与噪声频谱图 Fig. 3 The spectrum of gravity signal and noise

分别采用窗函数法设计的FIR低通滤波器、常规FFT滤波方法和傅里叶基追踪低通滤波器对加入噪声的重力扰动数据低通滤波处理，滤波的截止频率取0.01 Hz.FIR窗函数法采用汉宁窗，滤波器的长度根据主瓣宽度的一半等于截止频率来估算，估算出滤波器的长度为400，考虑到奇数长度的滤波器产生整数长度的延迟，故滤波器的长度取401.在傅里叶基追踪低通滤波方法中，时间序列的长度为实际的测线点数，傅里叶原子数选取的基本要求是大于时间序列长度，综合考虑算法的时间和精度，傅里叶原子数取为时间序列长度的3倍.在本文的算例中，傅里叶原子数都取为时间序列长度的3倍.然后根据截止频率，选取频率范围0~0.01 Hz的系数，根据公式(5)即可重构出低通滤波后的结果.为了减小滤波的边缘效应，对FIR窗函数法和常规FFT滤波方法采用镜像延拓在测线两端各补长200个数据，而傅里叶基追踪低通滤波方法不需要对原始数据补长.

从图 4可以看出，采用FIR窗函数法和FFT滤波法滤波，边缘效应非常明显，两端的数据与真实的重力扰动数据相差较大，这样的数据在实际应用中是必须要舍弃的.同时从图中可以看出傅里叶基追踪低通滤波方法边缘效应非常小，两端的数据与真实的重力扰动数据相差较小，基本上不需要考虑舍弃边缘部分的数据，而且滤波后得到结果与真实的重力扰动数据更加的接近.表 1详细列出了采用傅里叶基追踪低通滤波方法、FFT滤波法与采用FIR窗函数法滤波得到的结果和真实重力扰动数据的比较.从表中可以看出，采用FIR窗函数法滤波在两端的误差将近20 mGal，即使是去除了边缘数据的影响，在两端各舍弃200个数据，与真实值差值的均方差(RMS)也达到0.9 mGal；采用FFT滤波法得到的结果比FIR窗函数法稍差，去除两端200个数据后与真实值差值的RMS达到1.1 mGal；而采用傅里叶基追踪低通滤波方法在两端的误差最大值小于3 mGal，RMS为0.28 mGal，精度要高于去除了边缘影响的FIR窗函数法和FFT滤波法.图 5展示了这三种方法得到结果的频谱图与真实重力扰动数据频谱图的比较，从图中可以看出傅里叶基追踪低通滤波方法得到的频谱更加接近真实值，证明新方法能更好的压制高频噪声，有效分离低频信号和高频信号，减少有限时间序列的谱污染和谱泄漏.

图 4

Fig. 4

图 4 FIR窗函数法、FFT滤波法与FBPLPF方法比较 Fig. 4 Comparison of the result of FIR filters using window function and FBPLPF method

表 1 Table 1

表 1 FIR窗函数法、FFT滤波法、FBPLPF方法滤波结果与真实值比较(mGal) Table 1 Comparison of the result of FIR filter, FFT filter, FBPLPF method and real data (mGal) 滤波方法 Max Min Mean RMS FIR窗函数法 2.55 -19.32 -0.65 2.85 FIR窗函数法 (去除两端各200个数据) 2.55 -2.29 -0.03 0.94 FFT滤波法 2.87 -20.36 -0.66 3.17 FFT滤波法 (去除两端各200个数据) 2.17 -2.55 -0.05 1.07 FBPLPF方法 2.83 -0.24 0.04 0.28

表 1 FIR窗函数法、FFT滤波法、FBPLPF方法滤波结果与真实值比较(mGal) Table 1 Comparison of the result of FIR filter, FFT filter, FBPLPF method and real data (mGal)

图 5

Fig. 5

图 5 FIR窗函数法、FFT滤波法与FBPLPF方法结果的频谱图 Fig. 5 The spectrums of the result of FIR filter, FFT filter and FBPLPF

为了进一步验证傅里叶基追踪低通滤波方法滤波边缘效应小，结果精度高的特点，不失一般性，选取秦岭地区实测的东西测线T和南北测线L为例来探讨该方法的效果.

对于航空重力测量而言，其原理为式(6)，即

(6)

式(6)中，g_b为停机坪处的重力值，f、f⁰为比力观测值及其初值，为载体加速度，δg_Eötvös为厄特弗斯改正，δg_tilt为平台倾斜改正，δg_F为空间重力异常改正，δg_z为偏心改正，δg_k为零漂改正，γ₀为椭球面上的正常重力.根据原理式(6)做各项改正之后得到测线原始的空间重力扰动数据，其中比力观测值和载体垂向加速度的量级达到几万毫伽，包含大量的高频噪声，需要低通滤波才能提取出重力信号.

图 6中左图表示T测线原始重力扰动数据，从图中可以看出，经过各项改正得到的原始重力扰动数据包含有大量的噪声，噪声幅值最大能达到几万毫伽(mGal)，而真实的重力扰动幅值只有几十毫伽到一百多毫伽.右图比较了三种滤波方法的效果，其中红色双划线是GTGRAV软件经过一系列预处理，最后通过100 s低通滤波得到的测线重力扰动.飞行时已经考虑了拐弯处数据质量不好，因此飞机换测线时为了保证数据质量尽量拐大弯，而且在测区外多飞几公里，后期处理时将测区外的数据舍弃，因此认为处理得到的结果是真实可靠的，可以将其作为真实参考值.与仿真实验一样，滤波前将原始重力扰动数据采用镜像延拓在两端各补长200个数据，而采用傅里叶基追踪低通滤波方法不需要对数据补长.蓝色实线是采用FIR窗函数法经过低通滤波得到的结果，所有的参数设置与仿真实验一致；粉红色实线是采用FFT滤波法得到的结果；黑色点划线是傅里叶基追踪低通滤波方法的滤波结果，三种方法的截止频率都设为0.01 Hz，截止频率与GTGRAV软件一致.从图中可以看出，虽然采用了镜像延拓延长了数据，但是FIR窗函数法和FFT滤波法仍存在非常明显的边缘效应，边缘处与真实结果的误差达到100 mGal，而傅里叶基追踪低通滤波方法未采用数据补长，在边缘处与真实结果的误差小于4 mGal.图 7展示了L测线的结果，效果比T测线结果稍微差些.

图 6

Fig. 6

图 6 T测线原始重力扰动及滤波结果比较 Fig. 6 T survey line raw gravity disturbance and comparison of the result with different filters

图 7

Fig. 7

图 7 L测线原始重力扰动及滤波结果比较 Fig. 7 L survey line raw gravity disturbance and comparison of the result with different filters

表 2列出了采用三种滤波方法与GTGRAV软件比较的结果.由于FIR窗函数法两端的边缘效应太大，统计时把两端半个滤波窗口的数据去除，为了便于比较，FFT滤波法两端也各去除200个数据.从表中可以看出傅里叶基追踪低通滤波方法与FIR窗函数法舍弃了两端边缘数据后的结果精度相当，东西测线RMS为0.7 mGal，南北测线RMS为1.4 mGal，两种方法都能取得与GTGRAV软件相当的精度，而采用FFT滤波法即使舍弃边缘上的数据，精度仍要低于其它两种方法.南北测线精度较差可能是因为南北测线较短.精度明显要差于仿真实验，推测可能由于对仪器的很多硬件不太了解，处理数据时很多改正没有考虑完善.从图 6、图 7和表 2可以看出，傅里叶基追踪低通滤波方法能够有效减少有限时间序列的边缘效应，极大的提高数据利用率，对于花费昂贵的航空重力测量来说，这一点显得非常重要.

表 2 Table 2

表 2 FIR窗函数法、FFT滤波法、FBPLPF方法与GTGRAV软件处理结果比较(mGal) Table 2 Comparison of the result of FIR filter, FFT filter, FBPLPF method and the reference data of the GTGRAV software (mGal) 测线 滤波方法 Max Min Mean RMS 东西侧线T FIR窗函数法 (去除两端各200个数据) 3.37 -2.08 0.06 0.76 FFT滤波法 (去除两端各200个数据) 4.28 -5.48 -0.11 1.80 FBPLPF方法 3.75 -3.02 0.04 0.73 南北测线L FIR窗函数法 (去除两端各200个数据) 2.24 -2.73 -0.21 1.38 FFT滤波法 (去除两端各200个数据) 6.03 -5.05 -0.08 2.06 FBPLPF方法 6.91 -2.07 0.05 1.38

表 2 FIR窗函数法、FFT滤波法、FBPLPF方法与GTGRAV软件处理结果比较(mGal) Table 2 Comparison of the result of FIR filter, FFT filter, FBPLPF method and the reference data of the GTGRAV software (mGal)

针对航空重力测量数据滤波中的边缘效应问题，本文提出了一种新的滤波方法，即傅里叶基追踪低通滤波方法.结合仿真例子和实测的航空重力数据，比较了当前航空重力测量滤波中常用的FIR窗函数法，FFT滤波法和新方法的效果，结果表明新方法能够有效减少有限时间序列的谱污染和谱泄漏，在很大程度上降低滤波边缘效应，而且与FIR窗函数法去除了边缘效应后的结果精度相当，极大的提高航空重力实测数据的利用率.新方法的不足是计算速度比较慢，在处理器主频为2.2 GHz的普通计算机上采用Matlab语言编写程序，处理了51条测线，平均每条测线需要0.5 h.为了达到毫伽级的航空重力测量精度，需对多项误差进行事后改正，新方法计算速度较慢对于事后高精度数据处理是可以接受的.近十年来，我国自主研发的海空重力仪发展迅速，在硬件技术快速发展的同时，进一步完善自主研发的软件迫在眉睫，如何将新方法应用于国产航空重力仪数据处理中有待进一步研究.