《1 引言》

1 引言

惯性导航系统通常是参考椭球体来近似描述重力场, 对于地球大地水准面来讲具有良好的近似, 但对于局部地区, 如地形复杂的山地或海洋, 这种近似会与实际重力加速度存在一定误差。为了表述局部地区实际重力加速度与参考模型之间的误差, 引入了重力异常的概念。所谓重力异常是指重力平面某一点上的重力加速度与椭圆体表面大地测量学法线点上的标准重力加速度之间的差值, 目前它是高性能惯性导航系统的最大误差源。为了有效地消除重力异常对高精度惯性导航系统性能的影响, 必须利用安装在运动载体上的高精度重力仪实时测量重力异常, 并计算出重力加速度的垂线偏差, 以修正惯性导航系统的位置。重力仪测量的重力异常数据通常是强噪声背景下的微弱信号数据, 如何对被噪声严重污染的重力异常信号采取有效的方法进行检测处理, 是亟待解决的问题 [1]。笔者将自适应小波阈值去噪方法应用到重力仪信号处理中, 并将其与小波阈值去噪方法以及自适应卡尔曼滤波方法进行了对比。理论分析和仿真实验表明, 自适应小波阈值去噪方法、传统的小波阈值去噪方法和自适应卡尔曼滤波方法, 都能在一定程度上消除噪声信号对重力仪测量信号的影响, 但在相同情况下, 自适应小波阈值去噪方法具有明显的优越性。

《2 小波阈值及自适应小波阈值去噪原理 [2,3,4,5,6,7,8]》

2 小波阈值及自适应小波阈值去噪原理 [2,3,4,5,6,7,8]

《2.1小波阈值去噪原理及方法》

2.1小波阈值去噪原理及方法

小波变换是一种变尺度时频分析方法。设有函数f (x) ∈L2 (R) , L2 (R) 为二次方可积函数组成的Hilbert空间。通过连续小波变换, 可将函数f (x) 从时域变换到小波域, 并以不同尺度下一系列小波系数的形式来表征f (x) 。连续小波变换式为

Wψf(b,a)=|a|-1/2-+f(t)ψ¯(t-ba)dt(1)

其中f (t) ∈L2 (R) ;“-”表示复数共轭;ψa,b=|a|-1/2ψ(t-ba)称为小波函数, a为尺度因子, b为移位因子。根据Parseval恒等式x(t),g(t)12πx^(ω),g^(ω), 可得小波变换的反变换公式:

f(t)=Cψ-1-+-+a-2[Wψf(b,a)]ψa,b(t)dadb(2)

小波基不是唯一的, 只要满足容许性条件即可定义任一特定信号的小波基。为了数学上的方便, 小波变换也可以表示为

Wsf(t)=f\5ψs(t)=1s-f(x)ψ(t-xs)dx(3)

式中ψs(t)=1sψ(ts)s仍然是尺度参数。在实际应用中, 小波变换的尺度参数不必连续取值, 而是按照某种方式把连续小波及其变换做离散化处理。通常对尺度参数 s 进行二进制离散化, 即取 s=2j, jZ, 则 f (t) 在尺度 2j下的小波变换为

W2jf(t)=f(t)\5ψ2j(t)=12j-f(x)ψ(t-x2j)dx(4)

式 (4) 给出了第j个倍频程的局部信息。f (t) 的小波分解和重构可按Mallat塔式算法进行, 即小波变换分解公式为

{cj,n=kΖhk-2n*cj-1,kdj,n=kΖgk-2n*cj-1,k(k=0,1,2,,Ν-1)(5)

相应的重构公式为

cj-1,n=kΖhn-2kcj,k+kΖgn-2kdj,k(6)

式中cj, k为尺度系数;dj, k为小波系数;h, g为一对正交镜像滤波器组 (QMF) ;j 为分解层数;N为离散采样点数。

在利用小波进行去噪过程中, 首先要解决的问题是如何区分信号和白噪声。因此, 需要研究突变信号和白噪声在小波变换下的特性。设 X 为信号 f (x) 的局部突变点 (奇异点) , 则在该点处 f (x) 的小波变换取得模极大值, 在离散二进小波变换中有

|W2jf(x)|Κ(2j)α(7)

此时 j 为尺度参数, X 取离散值。由上式可知:

lb|W2jf(x)|lbΚ+αj(8)

若信号在X处的奇异性大于0 (即α>0) , 由式 (7) 可知随着尺度 j 的增加, 小波变换模极大值的对数也变大。若信号具有负的奇异性, 则情况恰恰相反, 负的奇异性意味着信号具有比不连续 (有界, α=0) 更差的奇异性, 这就是白噪声的情况。

n (x) 是实的方差为σ2 的宽平稳白噪声, W2jn (x) 是其二进小波变换, 小波ψ (x) 是实函数, 则W2jn (x) 也是一随机过程, 其方差:

E[|W2jn(x)|2]=--E[n(u)n(v)]ψ2j(x-u)ψ2j(x-v)dudv=-σ2ψ2j(x-u)du(9)

由于ψ2j(x)=12jψ(x2j), 由式 (9) 可知

E[|W2jn(x)|2]=ψ22jσ2(10)

式 (10) 表明|W2jn (x) |2的平均幅度反比于尺度 j, 即白噪声具有负奇异性。由此, 式 (8) 、式 (10) 成为区分信号和噪声在多尺度空间中模极大值传播行为的重要依据。

小波变换阈值法也称为“小波收缩”, 目前在信号去噪中, 应用最广的是Donoho等人提出的软阈值去噪法:假设1个叠加了高斯白噪声的有限长信号可以表示为

yi=xi+σzi(i=0,1,,n-1)

式中 zi 为一个标准的高斯白噪声;σ为噪声级。若要从被噪声污染的信号 yi 中恢复原始信号 xi, 则Donoho的去噪方法分为3个步骤:

1) 计算含噪声信号的正交小波变换。选择合适小波和小波分解层数 j, 将含噪信号运用式 (5) 进行小波分解至 j 层, 得到相应的小波分解系数 wjk, 1≤jJ

2) 采用软阈值法对分解得到的小波系数进行阈值处理:

w^j,k={sign(wj,k)\5(|wj,k|-λ)0|wj,k|λ|wj,k|<λ

式中λ=σ\5 (2logaN) 1/2;w^j, k为经阈值处理后的小波系数。

3) 进行小波逆变换。将经阈值处理过的小波系数用式 (6) 重构, 得到恢复的原始信号估计值。

《2.2自适应小波阈值去噪原理及方法》

2.2自适应小波阈值去噪原理及方法

虽然用Donoho的软阈值消噪法能获得不错的效果, 由于白噪声具有负的奇异性, 其幅度和稠密度随尺度增加而减少, 而信号则是相反。因此随着尺度级数的增加, 由噪声所控制的模极大值的幅度和稠密度会快速减少, 而信号控制的模极大值的幅度和稠密度会明显增大。因此在每一级尺度上都采用同一阈值显然不合适;在较低尺度上, 会去除有用信号, 而到了最大尺度级上, 则会留下一部分噪声。因此可以考虑选择自适应阈值来克服这种缺点。Donoho选择的单一阈值为

λ=lb(1+2Ν)J+ΖA

式中N 为预设的噪声功率;J 为所取的最大尺度;Z 为一常数, 取2;A 为最大的极值点幅度。可以将每一级尺度上都视为相互独立, 寻找一个与之最匹配的阈值来进行除噪, 此阈值的 N 还为预设功率, J 为本级尺度, Z 不变。A 为在本级尺度上搜索到的极值点幅度。最后, 将各级尺度上除噪后保留的模极大值点来重构信号。

与单一阈值法相比, 自适应阈值法可更有效地去除噪声并保留有用信号, 表现出更好的视觉效果。从L2 范数误差的观点看, 经该方法能够得出比单一法更高的信噪比。

《3 重力异常状态方程》

3 重力异常状态方程

根据文献[9], 地球重力场的变化可以看作零均值的随机过程, 该随机过程可以用白噪声激励的线性微分方程表示。设具有功率谱密度函数ϕu (ω) 的信号u (t) 通过传输函数为H () 的线性系统, 则输出信号的功率谱密度函数为ϕx (ω) =|H () |2ϕu (ω) 。若输入信号u (t) 为一白噪声过程, 则ϕu (ω) =1;同时当非确定随机信号功率谱函数ϕx (ω) 已知时, 可以计算得到H () 。因此, 如果已知随机过程的期望特征和协方差函数, 该随机过程完全可以用白噪声激励的线性系统来表示。假设一维随机过程的均值为零, 协方差为R (t) , 则该过程的状态方程为

{x(t)=F(t)x(t)+u(t)E[u(t)u(t+τ)]=Q(t)δ(t-τ)(11)

设随机过程是时间稳定的, 则F (t) =F (0) =F为常数, FQ可用下式计算:F=AP-1;Q=- (A+A′) ;A=ddt[R(t)]|t=0;P=R (t) |t=0。由于地球重力协方差模型是无法确切知道的, 因此其状态方程模型也是无法精确得到的。为了解决这一问题, 许多学者提出了各自的重力协方差模型的经验公式, 笔者选用形式简单并且比别的模型更接近实际情况的二阶高斯-马尔可夫异常位模型

CΤΤ(r)=σΤ2+(1+βγ)e-βr(12)

式中σΤ2, β为参数。由式 (12) 可计算重力异常协方差函数得:GΤrΤr(r)=2σΤ2β2(1-β2r)e-βr

假定舰艇行驶方向为x轴, 速度为v, 与之垂直并远离地心的方向为y轴, 则有x=vt, y=0。因此, 重力异常协方差函数变为

CΤrΤr(t)=2σΤ2β2(1-β2vt)eβvt

计算得

{A=d[C(t)]dt|t=0=-3σΤ2β3v;Ρ=C(t)|t=0F=AΡ-1=-1.5βv;Q=-(A+A)=6σΤ2β3v

重力异常状态方程为

{Δg(t)=-1.5βvΔg(t)+u(t)E[u(t)u(t-τ)]=6σΤ2β3vδ(t-τ)

其中u (t) 为重力异常畸变对重力异常产生的干扰, 其形式为白噪声序列, 离散化的重力异常状态方程为

{X(k)=-1.5βvX(k-1)+w(k)Y(k)=Ι\5X(Κ)+u(k)(13)

《4 仿真实验比较》

4 仿真实验比较

在仿真对比试验中, 自适应卡尔曼滤波器采用的是Sage-Husa滤波器, 该滤波器是一种改进型的卡尔曼滤波器。设重力异常状态方程为

{X(k+1)=φk,k-1X(k)+w(k)Y(k)=Ι\5X(k)(14)

假定R未知, 即 r=0, q=0, Q为常数, 则仿真试验中使用的该滤波算法的简化递推方程为 [10]:

计算加权系数dk= (1-b) / (1-bk) ;

状态一步预测方程Xk/k-1=φk,k-1X^k-1;

一步预测均方误差方程

Ρk/k-1=φk,k-1Ρk-1φk,k-1Τ+Qk-1

新信息序列方程vk=Yk-HkXk/k-1;

估计量测噪声

Rk=(1-dk)Rk-1+dk[vkvkΤ-ΗkΡk/k-1ΗkΤ]

滤波增益方程

Κk=Ρk/k-1ΗkΤ(ΗkΡk/k-1ΗkΤ+R)-1

状态估计方程X^k=X^k/k-1+Κkvk;

估计均方差方程Pk= (I-KkHk) Pk/k-1

式中 b 为遗忘因子, 一般取0<b<1。采用遗忘因子的目的在于限制滤波的记忆长度, 当X0, P0R0给定后, 可以运用上述迭代公式得到系统状态的估计值。

小波阈值去噪与自适应小波阈值去噪方法分别采用2.1节和2.2节中给出的方法。根据离散化的重力异常状态方程, 在仿真实验中假定受到白噪声污染的正弦波信号为重力测量信号。采用小波阈值去噪、自适应小波阈值去噪及自适应卡尔曼滤波分别对该重力仪测量信号进行处理, 处理前后的信号曲线如图1所示, 虚线表示处理前信号, 实线表示处理后信号。以处理后信号的信噪比作为评价处理方法好坏的依据, 经计算得, 原始信号的信噪比为-3.2028 dB, 经自适应卡尔曼滤波处理后信号的信噪比为-1.2651 dB, 经小波阈值去噪处理后信号的信噪比为0.8685 dB, 经自适应小波阈值去噪处理后信号的信噪比为2.8754 dB。因此, 仿真实验表明:自适应小波阈值去噪方法、小波阈值去噪方法以及自适应卡尔曼滤波都能在一定程度上消除噪声信号对重力仪测量信号的影响, 但在相同情况下自适应小波阈值去噪方法明显优于小波阈值去噪方法和自适应卡尔曼滤波。

《图1》

图1处理前和处理后的重力仪测量信号曲线比较

图1处理前和处理后的重力仪测量信号曲线比较  

Fig.1 The original gravity signal and processed gravity signal

《5 结论》

5 结论

为了有效地消除异常重力场对高性能惯性导航系统性能的影响, 抑制惯导系统给出的运动载体位置和速度误差的发散, 提高惯性导航系统的精度, 必须采用高精度重力仪实时测量重力异常值。由于重力异常信号通常是强噪声背景下的微弱信号, 因此必须采用有效的信号处理方法提取重力异常值。笔者在分析小波阈值及自适应小波阈值去噪算法的基础上, 将其应用到高精度海洋重力仪系统数据处理中, 并与自适应卡尔曼滤波进行对比, 以处理后信号的信噪比作为衡量信号处理方法优劣的依据。理论分析和仿真实验表明:自适应小波阈值去噪方法、小波阈值去噪方法以及自适应卡尔曼滤波都能在一定程度上消除噪声信号对重力仪测量信号的影响, 但在相同情况下自适应小波阈值去噪方法明显优于小波阈值去噪方法和自适应卡尔曼滤波。