海洋重力测量中滤波方法比较研究

刘啸星，蔡体菁，金伟明

(东南大学 仪器科学与工程学院,江苏 南京 210096)

0 引言

海洋重力辅助惯性导航系统是工作在强噪声环境下的动态测量平台，采集的数据中夹杂了不同频谱特性的噪声[1-2]。海洋重力测量的数据中噪声信号常比重力异常信号高出几百倍甚至上千倍。此外，重力异常信号主要集中在低频段，此频段内的噪声信号幅值通常低于重力异常幅值，而噪声信号主要集中在高频段并占取了大部分的带宽[3-4]。为了获得实际的重力异常信号，需要对测量数据进行低通滤波处理。

目前常用的滤波方法有有限脉冲响应(FIR)低通滤波、无限脉冲响应(IIR)低通滤波、Kalman滤波和小波滤波等[5-6]。周薇等[7]针对三轴惯性平台海空重力仪的特点，应用平滑卡尔曼滤波法对海洋重力测量数据进行了滤波处理，根据实测数据分析，其重复线内符合精度为0.2～0.4 mGal(1 Gal=10-2m/s2),空间分辨率达650 m。郎骏健等[8]基于FIR低通滤波器存在边缘效应而浪费部分有效数据的问题，提出了一种可抑制边缘效应的傅里叶基追踪低通滤波方法，该方法采用凸优化中的内点算法，可将低频与高频分量有效分离。根据仿真实验分析和实测数据验证，该方法东西测线和南北测线的均方根误差为0.7 mGal和1.4 mGal。本文研究比较了不同窗函数下FIR滤波器处理海洋重力测量数据的效果。

1 FIR低通滤波器

FIR滤波器的系统函数：

H(z)=b0+b1z-1+…+bNz-N=

(1)

式中：N为FIR滤波器单位脉冲响应h(n)的长度。H(z)在Z平面上有N-1个零点，在原点处有1个N-1重极点，所以H(z)永远稳定。假设理想低通数字滤波器在通带内的幅频特性|Hd(ejω)|=1，相频特性φ(ω)=0，截止频率ωc=2πfc，具有线性相位，群时延为α，则对应的时域滤波函数为hd(n)，频率响应为

(2)

FIR数字滤波器要求用有限长的单位脉冲响应h(n)去逼近无限长的理想滤波器的单位脉冲响应hd(n)。目前常用的方法是使用一个有限长的窗函数序列w(n)来截取hd(n)的主要部分h(n)=hd(n)·w(n),n=0,1,2,…,N-1,通过这种方式得到的频率响应e-jω近似于理想滤波器频率响应Hd(e-jω)。线性相位滤波器要求h(n)必须是偶对称的，对称中心是其长度的1/2，即h(n)=h(N-1-n)，而且α=(N-1)/2，所以要求窗函数w(n)也必须是中心偶对称的，即w(n)=w(N-1-n)。

由式(1)、(2)可知，基于窗函数的FIR低通滤波器的设计中，滤波器的截止频率、滤波长度以及窗函数的选取是关键因素。

2 滤波参数确定和窗函数选取

低通滤波器截止频率fc的选择应兼顾载体航行速度v和重力测量的分辨率λ，其关系为

(3)

fc一般用滤波周期Tc的倒数表示。根据采样频率fT可确定归一化截止频率为

(4)

考虑到虽然偶数阶滤波器的幅频响应比奇数阶滤波器的幅频响应好，但会产生非整数的时间延迟，滤波后需要进一步进行内插处理，通常滤波器的长度选择为奇数。令归一化截止频率等于窗函数主瓣宽度的一半，则可估算得到滤波器的长度，但该值只是一个估计值，需要以该值为基础，根据实际情况进行对比修正。对4种窗函数进行选取。

1) 三角形窗。长度为N的三角窗函数：

(5)

三角形窗的主瓣宽度为8π/N，虽然三角形窗的主瓣宽度大于矩形窗的主瓣宽度(4π/N)，但三角形窗的频谱密度函数永远是正值。

2) 汉宁窗又称为升余弦窗。长度为N的汉宁窗函数：

(6)

以上两部分之和使边瓣互相抵消，能量更集中在主瓣，边瓣衰减速度较快。但汉宁窗的主瓣宽度比矩形窗的主瓣宽度大1倍，即为8π/N。

3) 海明窗又称为改进的升余弦窗。当将汉宁窗改进为海明窗时，旁瓣可更小，长度为N的海明窗函数：

(7)

海明窗可将99.963%的能量集中在窗谱的主瓣内。与汉宁窗相比，主瓣宽度都是8π/N，边瓣的最大峰值变小了，小于主瓣峰值的1%，但是边瓣衰减速度相对较慢。

4) 布莱克曼窗又称为二阶2升余弦窗。长度为N的布莱克曼窗函数：

(8)

布莱克曼窗加上了余弦的二次谐波分量来进一步抑制旁瓣峰值，此时的主瓣宽度为12π/N。

窗函数的性能主要由主瓣宽度、旁瓣最大峰值及旁瓣渐近衰减速度来确定。理想的FIR低通滤波器应具有最小的主瓣宽度、最小的旁瓣峰值和最大的旁瓣衰减速度。在实际工程中不太符合理论要求，因此需要综合截止频率、滤波长度和窗函数的参数指标来确定滤波器的设计。4种窗函数的性能指标如表1所示。

表1 窗函数的基本参数

根据对上述4种窗函数的频率特性和窗谱性能指标的描述，三角形窗、汉宁窗和海明窗具有相同的主瓣宽度，旁瓣峰值依次减小。与三角窗相比，汉宁窗和海明窗的旁瓣衰减速度较大，但三角形窗的旁瓣峰值最大。因此，汉宁窗和海明窗更适合于处理海洋重力测量信号。布莱克曼窗具有最小的旁瓣峰值，旁瓣衰减速度也较快，满足最重要的处理条件，但主瓣宽度最大。因此,虽然汉宁窗、海明窗和布莱克曼窗均初步满足海洋重力测量数据的处理要求，但无法仅通过性能指标判断窗函数的优劣，需要通过实测数据的滤波处理结果进行对比来选取最合适的窗函数。

下面针对相同截止频率下，基于三角窗、汉宁窗、海明窗和布莱克曼窗4种窗函数设计的滤波器，通过对两组实测数据的比较与分析，选取最适用于海洋重力测量低通滤波的窗函数。

3 试验结果比较

选取某两个海区重力测量数据进行滤波方法比较。载体的航行速度均为6 m/s，采样频率fT=1 Hz。根据本次采集数据的空间分辨率的要求及式(3)可得截止频率为0.003 33 Hz，即滤波周期为300 s，根据式(4)可计算出相应的归一化截止频率为0.003 33 Hz。由窗函数的主瓣宽度和归一化截止频率可计算出各个滤波器的滤波长度。由于本次采用奇数长度的滤波器，这里三角窗、汉宁窗和海明窗的长度为601，布莱克曼窗的长度为901。分别使用以上各参数对测线数据进行滤波处理，分析不同的窗函数对于海洋重力测量数据处理的适用性，最佳窗函数的选择根据重复线内符合精度和测线交叉点不符值来进行评估。

图1、2分别为第一、二组数据滤波结果。由图1、2可看出，三角窗函数进行滤波后，测线数据仍有较明显的振荡，数据中高频噪声未能被完全滤除，这验证了前文所述的三角窗的旁瓣峰值大、旁瓣衰减速度小而不适合处理海洋重力测量数据的结论。汉宁窗、海明窗和布莱克曼窗由于旁瓣峰值小及旁瓣衰减速度快，对测线滤波后的数据相对平滑，可有效地去除高频噪声，由图还可看出三者的滤波效果水平相当。

图1 第一组数据滤波结果

图2 第二组数据滤波结果

基于以上4种窗函数滤波后的结果，分别对两组数据测区的航迹进行了截取，如图3、4所示。由图可看出，第一组测区的测线呈均匀的纵横分布，其中东西测线5条，南北测线6条，共计算交叉点30个。第二组测区的测线分布不均匀且形状不规则，此处截取东西测线4条，南北测线2条，共计算交叉点不符值8个。

图3 第一组数据测区截图

图4 第二组数据测区截图

对以上两组测区数据交叉点不符值进行计算，表2、3为4种不同滤波窗函数交叉点不符值的统计结果。根据交叉点不符值来选取合适的滤波窗函数。

表2 第一组数据交叉点不符值统计结果 单位：mGal

表3 第二组数据交叉点不符值统计结果 单位：mGal

由表2、3可看出，虽然三角窗、汉宁窗和海明窗的滤波长度一致，但三者中三角窗的滤波效果最差，其中误差分别为0.73 mGal和1.41 mGal。汉宁窗、海明窗和布莱克曼窗的滤波效果相当，其中汉宁窗和海明窗效果最好，其中误差均为0.47 mGal和0.41 mGal。布莱克曼窗的滤波长度最长，与另外两个窗函数的滤波结果相比，其在滤波过程中舍弃了测线头尾各一个滤波周期的数据，造成了有效数据的浪费。

在第一组数据的测区内选取一组南北方向的重复测线(见图5)，计算各窗函数滤波结果的重复线内符合精度，其统计结果如表4所示，汉宁窗、海明窗和布莱克曼窗仍表现一致，均符合海洋重力测量数据处理的精度需求，其中汉宁窗略优。

图5 第一组测区重复测线示意图

表4 重复线内符合精度统计结果

综上所述，基于汉宁窗、海明窗和布莱克曼窗设计的滤波器对于船测重力数据的滤波效果都适用，综合考虑有效数据的保留和交叉点不符值的结果，汉宁窗和海明窗对于海洋重力测量数据的处理是一个理想的选择。

4 结束语

采用低通滤波器来提取海洋重力测量数据中混合在高频噪声中的低频重力异常信息。本文给出了FIR低通数字滤波器截止频率和滤波器长度的确定方法，结合两组实测的海洋重力数据，通过对不同窗函数的滤波结果计算其重复线内符合精度和交叉点不符值，汉宁窗、海明窗和布莱克曼窗的计算精度均优于0.5 mGal，符合海洋重力测量数据处理的要求，综合考虑有效数据的保留和交叉点不符值的结果，选择汉宁窗和海明窗设计的FIR低通数字滤波器来提取海洋重力场信息。