海洋重力数据处理中的正反FIR滤波算法

蔡体菁,邵锦江,胡啸林

(东南大学 微惯性仪表与先进导航技术教育部重点实验室,江苏 南京 210096)

0 引言

在海洋重力场信号测量时,由于环境和仪器影响,大量的测量噪声也会被传感器所采集。目前在重力数据处理上有很多滤波技术,除了基于频域滤波的有限冲击响应(FIR)滤波、无限脉冲响应(IIR)滤波外,基于特征的小波滤波和基于贝叶斯估计的Kalman滤波方法[1]。

对重力测量数据滤波处理时,由于重力异常信号一般频率非常低,所以数字低通滤波器的应用比较广泛。FIR滤波具有线性相位,易于设计,IIR相位非线性,设计较复杂。因此,常采用FIR低通滤波器去剔除重力测量信号中的高频噪声[2-3]。郭圣焕[4]在比较汉宁窗FIR、切比雪夫等纹波FIR和巴特沃斯IIR低通数字滤波器后发现,得出切比雪夫纹波FIR低通滤波器具有更高的滤波精度。但FIR对于与有用信号频段接近的噪声无能为力,并且低通滤波需要一段数据才能进行一次滤波,因此导致了滤波结果的滞后性,这是不利于实际应用的。若能实现全频段线性相位,则滤波器不会出现非线性失真,这对后续的滤波组合算法有十分重大的意义。

1 线性相位FIR滤波器设计

1.1 基本原理与结构

FIR滤波器的卷积型表达式为

(1)

式中:N为FIR滤波器单位冲激响应h(k)的有限长度;阶数M=N-1。其系统函数H(z)在|z|>0处收敛,z平面内只有零点,N-1重极点在z=0处,因此是一个因果稳定的系统;无输出到输入的反馈,一般为非递归结构。

文献[5]指出单位冲激响应序列是偶对称或奇对称(即h(n)=±h(N-1-n))的滤波器具有线性相位,反之亦然。由此设计的线性相位网络结构滤波器,最多只需要进行N/2次乘法,降低了滤波器的计算复杂度。当N为奇数,h(n)为偶对称时,h(n)的相位函数和群延迟为

(2)

即y(n)表示(n-α)时刻数据的滤波结果,经分析可知FIR滤波器的相位延迟等于滤波器阶数的1/2。FIR滤波器的相位特性主要受h(n)的奇偶对称性影响,同时也受H(ω)影响。如果幅度函数H(ω)>0,则系统的相位等于-(N-1)/2·ω;如果H(ω)<0,则系统的相位等于-(N-1)/2·ω+π,从这个意义上来讲,这样的系统并不满足严格意义上的线性相频特性[6]。只有当对任意ω都有H(ω)>0时,系统才是严格意义上的线性相位系统。

1.2 基于窗函数的FIR滤波器设计方法

给定一个理想的零相位低通滤波器,设滤波的截止频率为ωc,其频率特性可表示为

(3)

对应的单位冲激响应[7-8]为

(4)

对应FIR滤波器设计的目标是寻找一个频响函数H(ejω)去逼近Hd(ejω),逼近的方式有频率采样法和时域的窗函数法。时域的窗函数设计思路是设计有限长的h(n)逼近hd(n)。最简单的逼近方法是用一个窗口截取一段hd(n),此时对应的有限长度的FIR冲激响应h(n)可以表示为hd(n)和一个窗函数w(n)的乘积[9],即:

h(n)=w(n)hd(n)

(5)

窗函数w(n)除了取最简单的矩形脉冲函数外,还可以取其他形式对hd(n)作一定的权重变化处理,来改善滤波器的特性。此外,H(ejω)的通带和阻带范围内的波动幅度取决于窗函数的傅里叶变换W(ejω)的旁瓣大小和数量,而过渡带宽则由其主瓣决定。例如汉宁窗表达式如下:

(6)

式中RN(n)为矩形窗函数。经计算可知,汉宁窗阻带最小衰减为-44 dB,旁瓣峰值为-31 dB,过渡带宽6.2 π/N。其他的窗函数还包括布莱克曼窗、三角窗、海明窗等等。实际应用中,根据允许的过渡带宽以及阻带衰减等参数,初步选择窗函数的形式以及FIR滤波器的阶数,再根据仿真结果进一步调整参数直至满足设计要求。

2 实时正反FIR滤波算法设计

2.1 实时正反FIR滤波公式推导

设N为奇数,h1(n)为偶对称序列,用小写z表示n时刻输入数据对应的滤波输出结果。结合式(1)与群延迟计算公式可知:

(m≥0)

(7)

z′1[m+(N-1)/2]表示m+(N-1)/2时刻的数据的滤波输出结果,当滤波输出数据量m≥(N-1)时,对输入的数据用同一滤波器进行反向FIR滤波,滤波过程实际上是一种线性卷积,考虑到相位延迟,设m时刻输入数据的滤波结果为z1(m),由于阶数M=N-1,根据式(1)形式有:

(8)

将式(7)代入式(8)可得:

(9)

式中m不小于j,根据多重求和的运算性质,有:

(10)

其中m≥(N-1),再令:

ai,j=h1(j)h1(i)

(11)

xi,j=x(m+i-j)

(12)

A=(ai,j)∈N×N

(13)

X=(xi,j)∈N×N

(14)

则:

z1(m)=Tr(AXT)

(15)

式中Tr()表示求迹,由于h1(n)是偶对称序列,则ai,j=aM-I,j=ai,M-j=aM-i,M-j=aj,i,为了简化式(10)的运算,由:

xi+δ,j+δ=x(m+i-j)

(16)

可得矩阵X所有行列连续的子矩阵对角线上的元素均相等,令k=i-j,则k∈{-M≤k≤M,k∈Z},则:

(17)

计算可知:

(18)

此时h″1(k)为非因果序列,化简h″1(k),同时用(-k+M)替换k可得:

(19)

取n=m+M,取N1=2N-1,有M1=2M,那么:

(20)

(21)

式中M=α1=(N1-1)/2表示滤波器的延迟。至此得到了正反FIR滤波器y1(n)式(1)形式的表达式,能够实时滤波,运算量从式(10)中的N2次乘法运算缩短为2N次乘法运算,进一步提高了正反FIR滤波的效率。分析冲激响应h′1(M1-k)=h′1(k),可知h′1(n)也为偶对称序列。因此,本节详细推导了一个窗长N1为奇数,冲激响应h′1(n)为偶对称序列的线性相位实时正反FIR滤波器表达式。

2.2 正反FIR滤波频率特性分析

采用单正向FIR滤波器y2(n)做为对比,其重启响应h2(n)也为偶对称序列,且采用与h1(n)相同类型的窗函数设计得到,与h′1(n)长度相同。设整个FIR滤波器的截至频率ωc=0.25π,滤波器长度N1=61,窗函数类型为汉宁窗,设计单正向FIR滤波器y2(n),其幅频与相频特性曲线如图1所示。

图1 汉宁窗60阶ωc=0.25π单正向FIR滤波器特性曲线

设计正反FIR滤波器时,考虑到h′1(n)需要与h2(n)序列长度相同,则h1(n)的序列长度N=(N1+1)/2=31,同样取截止频率ωc=0.25π,再通过式(18)计算得到h′1(n),则h′1(n)为长度N1=61的正反FIR滤波器y1(n)冲激响应。由于在计算过程中截止频率发生偏移,需要微调h1(n)对应的截止频率,使得y1(n)的截止频率与y2(n)的相同,最终得到正反向FIR滤波器y1(n),其幅频与相频特性曲线如图2所示。

图2 汉宁窗60阶ωc=0.25π正反FIR滤波器特性曲线

分析可知,正反FIR滤波最显著的特点是滤波器在所有的采样频谱都满足严格意义上的线性相频特性,不会对信号造成非线性失真,有利于后续对信号的进一步处理。其次正反FIR滤波的阻带衰减大,几乎是单正向FIR滤波的2倍,需要改进的地方是过渡带宽大,实际应用中要考虑选择合适的参数以满足设计要求。

2.3 正反FIR实时滤波方案设计

实时处理算法与后处理算法的区别:实时处理算法输入数据在经历过滤波器的延迟时间后即刻输出,而后处理算法则需要读取全部数据,处理完再输出,显然实时处理算法更适合处理船载试验中的海洋重力异常数据。

1) 实时正反FIR滤波算法首先要根据滤波器的截止频率,通带平稳度,允许的过渡带宽以及阻带衰减,初步选择窗函数与滤波器阶数。

2) 调参得到正反FIR滤波器的冲激响应h′1(n),将h′1(0),h′1(1),…,h′1(N1)存储到嵌入式计算的内存中(设为数组Resp),方便滤波时调用。

3) 开辟一个长度N1的内存空间(设为循环数组DinW),为当数据传输到FIR滤波器,按照滑窗的形式依次进入数组DinW,输入的第(k+1)个数据覆盖存储在DinW[k/N1]处。

4) 将数组Resp中的数据与DinW[(k+1)/N1]～DinW[(k+N1)/N1]组成的新的数组对应位置中的数据求和:

DinW[(k+i)/N1]

(22)

即为延迟了(N1-1)/2输出的滤波结果。综上所述,实时正反FIR滤波算法的示意图如图3所示。

图3 实时正反FIR滤波算法示意图

3 船载海洋重力试验

3.1 海洋重力数据特征

截止频率影响了滤波后重力异常的分辨率及数据的内符合精度。截止频率的选择应该兼顾载体航行速度和重力测量的分辨率,三者之间的关系为

λ=v/(2fc)

(23)

(24)

式中:v为载体的航行速度;λ为重力测量数据的空间分辨率;fc=1/Tc为低通滤波器的截止频率,Tc为滤波周期。

3.2 海洋重力数据处理船载试验

针对相同截止频率下,基于汉宁窗设计的实时正反FIR滤波器对实际数据滤波效果进行研究分析,其中评价窗函数性能通过重复线内符合精度进行评估。对比实验组1为实时单正向FIR滤波器,对比实验组2为离线(事后处理)正反FIR滤波器。

试验选取测线数据中载体航行速度为6 m/s,采样频率为fT=1 Hz,海洋重力测量要求重力异常数据空间分辨率达0.9 km,根据式(24)可以得到截止频率为0.003 33 Hz,即滤波周期为300 s,计算出相应的归一化截止频率为0.003 33 Hz。计算得到FIR滤波器的阶数N1=601,调参得到汉宁窗单正向FIR滤波器和实时在线和事后离线FIR滤波器。

用以上滤波器对某次海上实验采集到重力异常数据进行滤波处理,实时单正向FIR滤波结果如图4所示,实时和事后正反FIR滤波结果去程重复线如图5所示,返程重复线如图6所示。

图4 汉宁窗600阶fc= 0.003 33 Hz 单正向FIR滤波效果

图5 汉宁窗在线、离线正反FIR去程重复线

图6 汉宁窗在线、离线正反FIR返程重复线

分析可知单正向FIR滤波在输入数据噪声方差较大时,效果比较差,而实时正反FIR滤波则能很好的抑制输入噪声,并且其处理重力异常数据的重复线精度为0.09 mGal,与事后处理基本相当,符合海洋重力数据实时处理的要求。

4 结束语

本文提出一种实时正反FIR滤波算法,该算法设计出的滤波器在所有的采样频谱内都满足严格意义上的线性相频特性,不会对信号造成非线性失真,有利于后续对信号的进一步处理。该正反FIR滤波器的阻带衰减大,几乎是单正向FIR滤波的两倍,需要改进的地方是过渡带宽大,实际应用中要考虑选择合适的参数以满足设计要求。船载试验表明,将该实时正反FIR滤波器应用到海洋重力数据处理领域,滤波效果优于传统的单正向FIR滤波器,与事后正反FIR滤波效果相当。