CN102589551B - 一种基于小波变换的船用光纤陀螺信号实时滤波方法 - Google Patents

Abstract

本发明提供的是一种基于小波变换的船用光纤陀螺信号实时滤波方法。(1)利用陀螺信号采集回路实时测量光纤陀螺的敏感信息x(i)；(2)对滤波器进行赋初值和对称周期性拓展运算；(3)对a(i)进行Mallat多尺度小波分解，得到第j层小波分解系数；(4)小波阈值求取；(5)对分解系数d_j，n细节系数代入下式的阈值函数，求小波系数(6)利用小波系数进行小波重构，得到第J层重构；(7)计算滤波后输出信号值；(8)由陀螺采集回路采集下一个光纤陀螺输出值返回到步骤1，重新按步骤(1)-步骤(7)进行，完成光纤陀螺信号的实时滤波。本发明具有效率高，灵活性好，适用性强等优点。

Description

一种基于小波变换的船用光纤陀螺信号实时滤波方法

技术领域

本发明涉及的是一种陀螺信号噪声的处理，尤其是能对船用光纤陀螺信号进行实时滤波，滤除其噪声误差信号的方法。

背景技术

以光纤陀螺作为核心惯性元件构建的捷联式惯导系统具有体积小、重量轻、可靠性高、成本低，且能提供运载体角运动信息等优点，正逐步在各个领域替代传统的平台式惯导系统，成为未来惯性技术的发展趋势。由于光纤陀螺与载体的固联，其输出信号受环境影响包含了大量的随机漂移误差，影响了捷联惯导系统的方位、速度和位置精度。光纤陀螺漂移误差主要包括常值漂移和随机漂移两类，对于常值漂移，可通过器件或系统级标定予以补偿；而随机漂移受环境因素影响，表现出时变与弱非线性，很难建立其准确理想的数学模型。

由惯导系统误差分析可知，随机漂移激励的导航信息误差均方根值具有随时间发散的特性。因此，必须采取合适的数据滤波算法对陀螺随机漂移误差进行滤波处理，这对于提高惯导系统长时间导航定位精度，具有重要的作用。

现阶段国内外一些参考文献针对光纤陀螺噪声误差信号处理给出了一些数据器处理方案，例如低通滤波器、小波变换阈值滤波等，但是存在以下的问题：

（1）低通滤波器的方案对信号进行数字滤波，能提取有用的信号并滤去噪声信号。但是由于低通滤波器具有较大的延迟性，不能满足船用光纤陀螺信号处理的实时性要求。

（2）采用小波变换阈值滤波方案，多采用VisuShink阈值或SureShink阈值。VisuShink阈值选取准则在信号样本长度较大时滤波效果实时性差。SureShink阈值选取准则是均方差准则下的无偏估计，所求阈值更加能反映噪声特性；信号长度N较大时，SureShink阈值趋近于理想阈值，去噪效果较好；但在滤波初始阶段由于信号量少，无法采取SureShink阈值，因此影响到滤波效果。同时，随着时间的积累陀螺信号长度N会越来越大，会导致计算量增大，实时性差，滤波效率低。

发明内容

本发明的目的在于提供一种实时性好、去噪性强的基于小波变换的船用光纤陀螺信号实时滤波方法。

本发明的目的是这样实现的：

采用递推滑动窗的小波变换快速算法，处理边界问题时采用边界值重复的对称周期延拓方法；同时确定滑动窗口上、下限：最小窗口宽度L=64；最大窗口宽度为L=128。具体实施步骤如下：

（1）利用陀螺信号采集回路实时测量光纤陀螺的敏感信息x(i)；

（2）对滤波器进行赋初值和对称周期性拓展运算：

$a\left(i\right)\left\{\begin{array}{cc}x\left(i\right)& k<L,i=\mathrm{1,2},...,L\\ x(k-L+i)& k\&GreaterEqual;L,i=\mathrm{1,2},...,L\\ x(k+L-i)& k\&GreaterEqual;L,i=L,L+\mathrm{1,2},...,2L\end{array}\right.$

式中，L为窗口宽度，a(i)为对敏感信息x(i)进行初值和周期延拓处理后的信号数列；

（3）对a(i)进行Mallat多尺度小波分解，得到第j层小波分解系数，

$\left\{\begin{array}{c}{c}_{j,n}=\underset{l}{\mathrm{\&Sigma;}}{h}_{l-2n}{c}_{j+1,l}\\ d_{j,n}=\underset{l}{\overset{l}{\mathrm{\&Sigma;}}}{g}_{l-2n}{c}_{j+1,l}\end{array}\right.$

（4）小波阈值求取

其中，N为一个海浪周期内的信号采样数，f为陀螺信号采样频率，[]为取整；λ_VisuShink为基于VisuShink阈值选取准则求取的阈值，λ_SureShink为基于SureShink阈值选取准则求取的阈值；在输入数据少于N时，求取λ_VisuShink用于数据滤波；当输入数据达到N或者N的整数倍时，求取出λ_SureShink，将λ_SureShink用于下一周期内N组数据的滤波；依此类推由第i个周期N个数据求取一个阈值用于第i+1周期数据的滤波处理；

（5）对分解系数d_j,n细节系数代入下式的阈值函数，求小波系数

${\hat{d}}_{\mathrm{jk}}=\left\{\begin{array}{cc}0& \left|d_{j,k}\right|\&le;{\mathrm{\&lambda;}}_1\\ \mathrm{sgn}\left(d_{j,k}\right)\frac{{\mathrm{\&lambda;}}_2\left(\right|d_{j,k}|-{\mathrm{\&lambda;}}_1)}{{\mathrm{\&lambda;}}_2-{\mathrm{\&lambda;}}_1}& {\mathrm{\&lambda;}}_1\&le;{|d}_{j,k}|\&le;{\mathrm{\&lambda;}}_2\\ d_{j,k}& \left|d_{j,k}\right|>{\mathrm{\&lambda;}}_2\end{array}\right.;$

（6）利用小波系数进行小波重构，得到第J层重构为：

C_J=H*C_J-1+G*D_J-1

其中，H*，G*分别为H，G的共轭转置矩阵，滤波后的信号为：

$\hat{a}\left(1\right),\hat{a}\left(2\right),...\hat{a}(2L-1),\hat{a}\left(2L\right)$

（7）计算滤波后输出信号值：

$\hat{x}\left(k\right)=\frac{\hat{a}\left(L\right)+\hat{a}(L+1)}{2}$

其中，即为小波实时滤波输出值；

（8）由陀螺采集回路采集下一个光纤陀螺输出值x(k+1)，返回到步骤1，重新按步骤（1）-步骤（7）进行，完成光纤陀螺信号的实时滤波。

本发明的特点：

（1）采用小波变换来进行滤波处理，可以对信号进行任意细节的分解，能够很好处理微弱信号，同时实现滤波处理的实时性，克服低通滤波器的延迟性。

（2）采用VisuShink阈值选取准则与SureShink阈值选取准则相结合的阈值选取方法；克服了VisuShink阈值的大小受样本长度N影响较大、滤波效果较差的缺点。同时也克服了在滤波初始阶段，由于船用光纤陀螺信号输出信号较少，无法采用SureShink阈值准则的问题。

（3）由于受海浪扰动的船用光纤陀螺信号具有频率低、周期性的特性，且SureShink阈值是均方差准则下的无偏估计，这样在一个海浪扰动周期内用SureShink阈值准则求出阈值，将此阈值用于下一个即将到来的一个周期内的数据滤波；由于所求阈值更加能反映噪声特性，同时也减小了阈值求取范围，减少了计算量、提高了滤波效率。

本发明提供了分阶段不同阈值准则选取的方法，在各个阶段充分的发挥了各个阈值准则的优点；同时提出了在一个海浪周期内采用基于均方差准则的SureShink阈值选取准则求一次阈值的方案，这样能更好的滤除噪声，使所设计的实时滤波算法针对性更强，去噪效果更好。即该发明具有效率高，灵活性好，适用性强等优点。

附图说明

图1为船用光纤陀螺信号的频谱曲线；

图2为船用光纤陀螺信号特征分析图；

图3为本发明滤波器实时性验证曲线；

图4为阈值准则选取流程图；

图5为周期性SureShink阈值求取及应用流程；

图6为基于小波变换实时滤波算法流程图；

图7（a）-图7（c）为采用本发明中滤波算法进行滤波的曲线仿真图；其中图7（a）为东向陀螺输出信号及滤波后信号曲线仿真图，图7（b）为北向陀螺输出信号及滤波后信号曲线仿真图，图7（c）为方位陀螺输出信号及滤波后信号曲线仿真图。

具体实施方式

下面结合附图举例对本发明做更详细的描述：

1.光纤陀螺信号的分析

与运载体固联的光纤陀螺，输出信号为运载体相对惯性坐标系的旋转角速率，在动基座下，陀螺的输出信号为：

${\mathrm{\&omega;}}_{\mathrm{ib}}^b=C_n^b({\mathrm{\&omega;}}_{\mathrm{ie}}^n+{\mathrm{\&omega;}}_{\mathrm{en}}^n)+{\mathrm{\&omega;}}_{\mathrm{nb}}^b+{\mathrm{\&epsiv;}}_0+\mathrm{\&delta;\&epsiv;}---\left(8\right)$

式中，为地球旋转引起的角速率；为舰船线运动引起的角速率；为海洋环境因素扰动引起的角速率；ε₀是陀螺的随机常值误差项，一般与陀螺的本身设计有关系；δε是陀螺的随机噪声，δ为噪声强度，C_n ^b为地理坐标系与载体坐标系之间的转换矩阵。

由（8）式所示信号的特性可知，地球旋转角速率为恒定值，在固定纬度处的地理坐标系上的投影为恒定值。由于舰船的运动速度低，因此舰船线运动引起的角速率在地理坐标系上的投影也属低频信号，一般分布在10^-6～10^-5Hz频带范围内。而受到海浪等扰动的影响，舰船的运动可由一系列幅值和频率相近的正弦波来描述：

$\underset{i=0}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\sin (2\mathrm{\&pi;}{f}_i+{\mathrm{\&beta;}}_i),0<{\mathrm{\&beta;}}_i<2\mathrm{\&pi;}---\left(9\right)$

式中，频率f_i的范围为0.07～0.22Hz，幅值α_i的范围为5～37mrad。结合实际航行时海况级别，在短时间内，可以近似认为舰船仅仅受某一频段内的海浪扰动的影响。

对实测的船光纤陀螺数据进行处理，分析信号的频谱可知，2Hz以下主要是低频的有用信号，而在10～50Hz范围内分布着强度很高的噪声。

将船用光纤陀螺信号频段分为有用信号段、过渡段和噪声信号段，如图2所示，其有用信号频率低于2Hz，但在设计滤波器时，为留有一定的裕度，因此设置一过渡段，以避免对特殊运动状态下陀螺信号的错误处理。

2.基于船用光纤陀螺信号阈值选取准则的设计

小波变换阈值滤波方案中主要的阈值选取准则有：VisuShink、RiskShink、SureShink等。由于船用光纤陀螺信号受海浪扰动信号频率低且短时间内具有周期性，而SureShink阈值是均方差准则下的无偏估计。因此首先确定选取SureShink阈值作为其阈值选取准则；但是取SureShink阈值选取准则最佳选择范围一般较大，对于实时小波滤波无法实现，因此在实际应用中，可以适当的缩小阈值选择范围，得到SureShink准则下的次优解。选取SureShink阈值选取准则的范围为一个海浪扰动周期内采集到的船用光纤陀螺信号的数据量；由于海浪扰动信号具有周期性，因此在一个海浪周期内基于SureShink阈值选取准则求出的阈值具有代表性，能够充分反映出一个周期内信号的特点，将其用于下一个周期船用光纤陀螺信号的滤波处理时效果好。阈值范围的求取公式：

N=f/f'     （10）

式中，f'为船舶受海浪摇摆运动的频率；f船用光纤陀螺信号的采样频率。

在滤波处理初始阶段，由于船用光纤陀螺输出信号数据量少，采用SureShink阈值会导致阈值求取不准确影响滤波效果。因此在采集信号量较少时采用VisuShink阈值选取准则，弥补在滤波初始阶段SureShink阈值无法有效滤波的缺点，其计算公式如下

$\mathrm{\&lambda;}=\mathrm{\&sigma;}\sqrt{2\log N^{\&prime;}}---\left(11\right)$

式中，σ噪声信号的标准差、N'为信号的长度。

将两种阈值选取准则有效地结合起来，在滤波的各个阶段充分发挥各自的优点，使所设计的实时滤波算法针对性更强，去噪效果更好。

3.实时滤波算法流程

为满足对船用光纤陀螺采集数据的实时处理，同时船用光纤陀螺信号具有周期性，首先确定基于递推滑动窗的小波变换快速算法，处理边界问题时采用边界值重复的对称周期延拓方法；同时确定滑动窗口上、下限：最小窗口宽度L=64；最大窗口宽度为L=128。

其实时滤波算法流程图如附图6所示，现对其步骤进行详细解释：

SureShink阈值选取准则求出的阈值范围的求取公式：

N=f/f'     （12）

式中，f'为船舶受海浪摇摆运动的频率；f船用光纤陀螺信号的采样频率。

已知舰船受海浪扰动信号的频率范围为0.07-0.22Hz，则确定一个海浪周期内信号采集范围为：

$\frac{1}{0.22}f~\frac{1}{0.07}f$

即确定一个海浪周期内信号采集范围为4.5f～14.3f。为同时保证小波分解数据量为2的整数次幂以及信号采集不失真条件，一次海浪周期内陀螺信号采样量为：

$N=2^{\left[{\log }_2(2\&times;14.3f)\right]}$

式中，[]为取整运算。

（1）利用陀螺信号采集回路实时测量光纤陀螺的敏感信息x(i)。

（2）对滤波器进行赋初值和对称周期性拓展运算：

$a\left(i\right)\left\{\begin{array}{cc}x\left(i\right)& k<L,i=\mathrm{1,2},...,L\\ x(k-L+i)& k\&GreaterEqual;L,i=\mathrm{1,2},...,L\\ x(k+L-i)& k\&GreaterEqual;L,i=L,L+\mathrm{1,2},...,2L\end{array}\right.---\left(13\right)$

式中，L为窗口宽度，a(i)为对敏感信息x(i)进行初值和周期延拓处理后的信号数列。

（3）对a(i)进行Mallat多尺度小波分解，得到第j层小波分解系数，

$\left\{\begin{array}{c}{c}_{j,n}=\underset{l}{\mathrm{\&Sigma;}}{h}_{l-2n}{c}_{j+1,l}\\ d_{j,n}=\underset{l}{\overset{l}{\mathrm{\&Sigma;}}}{g}_{l-2n}{c}_{j+1,l}\end{array}\right.---\left(14\right)$

式中，c_j+1,_n为近似系数，d_j+1,_n为细节系数，把（14）式表示为矩阵形式为：

（4）小波阈值求取流程如附图4所示：

式中，λ_VisuShink为基于VisuShink阈值选取准则求取的阈值，λ_SureShink为基于SureShink阈值选取准则求取的阈值，其中λ_VisuShink按（17）式求取：

${\mathrm{\&lambda;}}_{\mathrm{VisuShink}}=\mathrm{\&sigma;}\sqrt{2\log (N-1)}$

式中，是σ噪声信号的标准差、N为一个周期内船用光纤陀螺信号的采集量。

SureShink阈值求取，其阈值具体求取过程如下：

1）对前N个数据进行小波分解，求出小波系数；

2）将第j层的小波按系数平方由大到小排列，得到一个新的向量：

$\begin{array}{ccc}S\&Element;\{s_1{s}_2& ...& s_{2^j}\},\end{array}$ 其中 $s_1\&le;s_2\&le;...\&le;s_{2^j},$ 2^j为小波系数的个数；

3）计算风险向量其中

$r_i=\frac{2^j-2i+(2^j-i)s_i+\underset{m=1}{\overset{i}{\mathrm{\&Sigma;}}}{s}_m}{2^j}---\left(18\right)$

4）以R中的最小元素r_B作为风险值，由r_B的相应位置B求出的s_B，则阈值为λ_SureShink=σ_n(s_B)^1/2；σ_n为噪声的均方差；

依此类推由第i个周期N个数据求取一个阈值用于第i+1周期数据的滤波处理，具体流程图如图5所示。

（5）对分解系数dj,n细节系数代入（19）式所示的阈值函数，求小波系数

${\hat{d}}_{\mathrm{jk}}=\left\{\begin{array}{cc}0& \left|d_{j,k}\right|\&le;{\mathrm{\&lambda;}}_1\\ \mathrm{sgn}\left(d_{j,k}\right)\frac{{\mathrm{\&lambda;}}_2\left(\right|d_{j,k}|-{\mathrm{\&lambda;}}_1)}{{\mathrm{\&lambda;}}_2-{\mathrm{\&lambda;}}_1}& {\mathrm{\&lambda;}}_1\&le;{|d}_{j,k}|\&le;{\mathrm{\&lambda;}}_2\\ d_{j,k}& \left|d_{j,k}\right|>{\mathrm{\&lambda;}}_2\end{array}\right.---\left(19\right)$

（6）利用式（19）所示的小波系数进行小波重构，得到第J层重构公式为：

C_J=H*C_J-1+G*D_J-1     （20）

式中，H*，G*分别为H，G的共轭转置矩阵。则滤波后的信号可以表示为：

$\hat{a}\left(1\right),\hat{a}\left(2\right),...\hat{a}(2L-1),\hat{a}\left(2L\right)---\left(21\right)$

（7）计算滤波后输出信号值：

$\hat{x}\left(k\right)=\frac{\hat{a}\left(L\right)+\hat{a}(L+1)}{2}---\left(22\right)$

式中，即为小波实时滤波输出值。

（8）由陀螺采集回路采集下一个光纤陀螺输出值x(k+1)，返回到步骤1，重新按步骤（1）-步骤（7）进行，完成光纤陀螺信号的实时滤波。

实例分析：

1.实时性验证

近似地模拟一正弦海浪模型，频率为0.2Hz，并加入白噪声信号，采样频率为100Hz。滤波后的信号如附图3所示，其中单一采样点滤波运行时间为3.8ms小于采样间隔10ms，滤波进行30s后，可以得出其延迟时间非常小，采用此算法对光纤陀螺信号处理后基本不发生延迟，可知采用本发明设计的小波实时滤波算法完全满足要求。

2、小波实时滤波算法的验证

对某型号光纤陀螺进行舰船系泊试验，采样频率为100Hz，x轴输出的原始信号以及经过本文设计的小波滤波器滤波后的信号如附图7所示。由滤波前后的数据比对结果可知，滤波效果比较明显，噪声得到了有效抑制，实时性好。

Claims (1)

1.一种基于小波变换的船用光纤陀螺信号实时滤波方法，其特征是具体实施步骤如下：

(1)利用陀螺信号采集回路实时测量光纤陀螺的敏感信息x(i)；

(2)对滤波器进行赋初值和对称周期性拓展运算：

$a\left(i\right)=\left\{\begin{array}{cc}x\left(i\right)& k<L,i=\mathrm{1,2},...,L\\ x(k-L+i)& k\&GreaterEqual;L,i=\mathrm{1,2},...,L\\ x(k+L-i)& k\&GreaterEqual;L,i=L,L+1,...,2L\end{array}\right.$

式中，L为窗口宽度，a(i)为对敏感信息x(i)进行初值和周期延拓处理后的信号数列；

(3)对a(i)进行Mallat多尺度小波分解，得到第j层小波分解系数，

$\left\{\begin{array}{c}{c}_{j,n}=\underset{l}{\mathrm{\&Sigma;}}{h}_{l-2n}{c}_{j+1,l}\\ d_{j,n}=\underset{l}{\mathrm{\&Sigma;}}{g}_{l-2n}{c}_{j+1,l}\end{array}\right.$

式中，j,l,n∈Z，c_j,n为第j层小波尺度分解系数，d_j,n为第j层小波细节分解系数，h_l-2n为低通滤波器冲击响应，g_l-2n为高通滤波器冲击响应，

(4)小波阈值求取

$\mathrm{\&lambda;}=\left\{\begin{array}{cc}{\mathrm{\&lambda;}}_{\mathrm{VisuShink}}& k<N\\ {\mathrm{\&lambda;}}_{\mathrm{SureShink}}& k\&GreaterEqual;N\end{array}\right.$ 且k％N＝＝0

其中，N为一个海浪周期内的信号采样数，，f为陀螺信号采样频率，此处[]为取整；λ_VisuShink为基于VisuShink阈值选取准则求取的阈值，λ_SureShink为基于SureShink阈值选取准则求取的阈值；在输入数据少于N时，求取λ_VisuShink用于数据滤波；当输入数据达到N或者N的整数倍时，求取出λ_SureShink，将λ_SureShink用于下一周期内N组数据的滤波；依此类推由第i个周期N个数据求取一个阈值用于第i+1周期数据的滤波处理；

(5)对分解系数d_j,n细节系数代入下式的阈值函数，求小波系数

${\hat{d}}_{\mathrm{jk}}=\left\{\begin{array}{cc}0& \left|d_{j,k}\right|\&le;{\mathrm{\&lambda;}}_1\\ \mathrm{sgn}\left(d_{j,k}\right)\frac{{\mathrm{\&lambda;}}_2\left(\right|d_{j,k}|-{\mathrm{\&lambda;}}_1)}{{\mathrm{\&lambda;}}_2-{\mathrm{\&lambda;}}_1}& {\mathrm{\&lambda;}}_1\&le;\left|d_{j,k}\right|\&le;{\mathrm{\&lambda;}}_2\\ d_{j,k}& \left|d_{j,k}\right|>{\mathrm{\&lambda;}}_2\end{array}\right.;$

式中，λ₁为设计的小波阈值，λ₂为下阈值；

(6)利用小波系数进行小波重构，得到第J层重构为：

C_J＝H^*C_J-1+G^*D_J-1

其中，H^*，G^*分别为H，G的共轭转置矩阵；H，G均为滤波系数矩阵；滤波后的信号为：

$\hat{a}\left(1\right),\hat{a}\left(2\right),\&CenterDot;\&CenterDot;\&CenterDot;\hat{a}(2L-1),\hat{a}\left(2L\right)$

其中，式中，C_J＝{c_j,k}_k∈Z，D_J＝{d_j,k}_k∈Z；

(7)计算滤波后输出信号值：

$\hat{x}\left(k\right)=\frac{\hat{a}\left(L\right)+\hat{a}(L+1)}{2}$

其中，即为小波实时滤波输出值；

(8)由陀螺采集回路采集下一个光纤陀螺输出值x(k+1)，返回到步骤1，重新按步骤(1)-步骤(7)进行，完成光纤陀螺信号的实时滤波。