一种基于交叠M带离散小波变换的光纤陀螺信号去噪方法

摘要

一种基于交叠M带离散小波变换的光纤陀螺信号去噪方法，首先建立光纤陀螺输出信号模型，再对光纤陀螺输出信号进行OMDWT小波分解，然后对OMDWT小波系数进行门限阈值处理，最后对光纤陀螺输出信号进行OMDWT小波重构，得到去噪后的陀螺信号。本发明在利用OMDWT对光纤陀螺输出信号进行去噪处理时，OMDWT的小波系数和尺度系数随着光纤陀螺信号循环移动而循环移动相应的节拍；在变换过程中，陀螺信号按多通道进行分解，按不同频带分解的速度较快，对陀螺信号有更细的频带划分，去噪效果优于传统小波去噪效果。本发明对于惯导系统导航性能的提高具有重要意义。

说明书

技术领域

本发明涉及一种基于交叠M带离散小波变换(Overlap M-band Discrete  Wavelet Transform，OMDWT)的光纤陀螺信号去噪方法，属于惯性导航技术 领域。

背景技术

光纤陀螺是以Sagnac效应为基础而发展起来的新型全固态陀螺仪，受 到世界各国的高度重视，已成为新一代中低精度捷联惯性导航、制导系统中 的理想惯性器件。捷联式惯导系统中的一个关键技术是如何完全自主地实现 快速初始对准，即快速自主地确定载体的初始方位角和水平姿态角。影响方 位对准(寻北)精度的主要因素是陀螺漂移，有效地消除陀螺漂移是保证方 位对准精度的关键。

光纤陀螺漂移可分为系统性漂移和随机漂移两种类型。由于随机漂移是 弱非线性的、慢时变的，同时受外部环境等不确定性因素的影响，无法建立 其精确系统模型，在惯导系统中不能用简单的方法加以补偿，因此随机漂移 成为衡量陀螺仪精度的重要指标。为了提高初始对准精度，必须采用有效的 信号处理方法对陀螺信号处理，来抑制这些随机漂移。目前对陀螺信号的处 理有两种思想：(1)建立陀螺随机漂移模型，使用Kalman滤波、Wiener 滤波等方法进行补偿；(2)直接对陀螺输出信号进行消噪处理。常用方法 有数字低通滤波、自适应滤波和小波变换阈值滤波等。鉴于随机漂移的特性， 事先不可能得到精确的统计特性，需要采用直接滤波方法对陀螺信号进行去 噪处理。

小波具有较好的时频分辨特性，因此被应用到光纤陀螺信号的去噪中。 现有的基于小波理论的光纤陀螺信号分析研究中，都是基于2带小波的，去 噪结果具有一定的局限性，主要体现在两方面：一方面，随着分解的进行， 小波系数的长度以2倍关系递减，分解速度慢；另一方面，光纤陀螺信号循 环移动一定节拍后，其离散小波变换(Discrete Wavelet Transform，DWT)小波 系数和尺度系数不能循环移动同样节拍，去噪效果差。

发明内容

本发明的技术解决问题是：本发明的技术是针对已有的光纤陀螺信号分 析方法中的不足，将交叠M带离散小波变换引入陀螺信号去噪。一方面， 克服了基于2带小波的光纤陀螺信号处理分解速度慢的不足，M带小波理论 为小波基的选取提供了更广泛的选择余地；在分解子带数相同条件下，M带 离散小波变换(M-band Discrete Wavelet Transform，MDWT)中信号是按多通 道进行分解的；另一方面，克服了基于2带小波进行光纤陀螺信号去噪时， 光纤陀螺信号循环移动一定节拍后，DWT小波和尺度系数不能循环移动同 样节拍的不足。

本发明的技术解决方案是：本发明提出了一种基于交叠M带离散小波 变换的光纤陀螺信号去噪方法，给出了计算有效的光纤陀螺信号分解、门限 阈值处理和重构算法。具体步骤如下：

(1)建立光纤陀螺输出信号模型

其中，X_t为光纤陀螺输出信号；t＝0，...，N-1，N表示光纤陀螺输出信号 的长度；ω_ie为地球自转角速率，ω_ie＝7.27×10^-5rad/s；为测试点的地理纬度； ε_t为陀螺漂移，陀螺漂移由常值分量、周期分量和白噪声组成，即：

ε_t＝ε_d+Ω_dsin(2πf_d+θ₀)+W_t

其中，ε_d为常值分量，短时间内近似为一个常数；Ω_d为周期分量的幅值； f_d为周期分量的频率；θ₀为初始相位；W_t为零均值高斯白噪声。式中的 Ω_dsin(2πf_d+θ₀)和W_t两项构成了光纤陀螺输出的随机漂移项。

(2)对光纤陀螺输出信号作OMDWT小波分解

对步骤(1)得到的长度N的光纤陀螺输出信号X_t作OMDWT变换，将 OMDWT变换所得到的第j尺度的OMDWT小波系数和尺度系数分别用N维 矢量和表示，则对应的元素和对应的元素分别为：

${\underset{‾}{\overset{~}{W}}}_{j,t}^i\equiv \underset{l=0}{\overset{L_j-1}{\Sigma }}{\tilde{h}}_{j,l}^i{X}_{t-l\mathrm{mod}N},$${\tilde{V}}_{j,t}\equiv \underset{l=0}{\overset{L_j-1}{\Sigma }}{\tilde{g}}_{j,l}{X}_{t-l\mathrm{mod}N}$

式中，t＝0，...，N-1，N表示光纤陀螺信号的长度；l＝1，...，L，L＝2M表示 第J₀尺度上MDWT滤波器的长度，M表示小波带数，即M带小波； i＝1，...，M-1；j＝J₀，J₀+1，...，J，J₀表示小波变换的初始尺度，J表示小波变换 的终止尺度；和分别为第j尺度上OMDWT的第i个小波滤波器系 数和尺度滤波器系数；lmodN表示l对N取余。

${\tilde{h}}_{j,l}^i\equiv h_{j,l}^i/M^{j/2},$${\tilde{g}}_{j,l}\equiv g_{j,l}/M^{j/2}$

式中，和{g_j，l}分别表示第j尺度上MDWT的第i个小波滤波器系数 和尺度滤波器系数。M^j/2表示M的j/2次方；第j尺度上滤波器的长度是 L_j≡(M^j-1)(L-1)+1。

对光纤陀螺信号X_t分别用周期滤波器进行滤波，可利用光纤陀螺信号在 第j-1尺度上的尺度系数矢量计算出光纤陀螺信号在第j尺度上的小波 系数矢量和尺度系数矢量的元素。和对应的元素分别用 和表示，则：

${\underset{‾}{\overset{~}{W}}}_{j,t}^i\underset{l=0}{\overset{L-1}{\Sigma }}{\tilde{h}}_l^i{\tilde{V}}_{j-1,t-M^{j-1}\mathrm{mod}N},t=\mathrm{0,1},...,N-1$

${\tilde{V}}_{j,t}=\underset{l=0}{\overset{L-1}{\Sigma }}{\tilde{g}}_l{\tilde{V}}_{j-1,t-M^{j-1}\mathrm{mod}N},t=\mathrm{0,1},...,N-1$

式中，和分别表示OMDWT的第i个小波滤波器系数和尺度滤 波器系数；表示t-M^j-1lmodN时刻第j-1尺度的尺度系数， j＝J₀，J₀+1，...，J。

在上述变换中，OMDWT小波滤波器系数和尺度滤波器系数的求解 如下：

${\tilde{h}}_l^i\equiv h_l^i/\sqrt{M},$${\tilde{g}}_l\equiv g_l/\sqrt{M}$

式中，l＝1，...，L；和{g_l}分别为MDWT的第i个小波滤波器和尺度滤 波器系数。

步骤(2)对光纤陀螺信号作OMDWT小波分解时，按不同频带分解从 而速度较快，且对光纤陀螺信号有更细的频带划分；光纤陀螺信号循环移动 后，其OMDWT小波系数和尺度系数会循环移动相应节拍，有利于光纤陀 螺的信号去噪。

(3)对OMDWT小波系数作门限阈值处理

对步骤(2)中光纤陀螺信号分解后得到的小波系数作门限阈值处理。 利用门限阈值处理，将绝对值小于门限阈值的小波系数置零，但对绝对值大 于门限阈值的小波系数作了收缩。

${\tilde{W}}_{j,t}^i({\underset{‾}{\overset{~}{W}}}_{j,t}^i,T_s)=\left(\begin{array}{cc}{\underset{‾}{\overset{~}{W}}}_{j,t}^i-T_s,& {\underset{‾}{\overset{~}{W}}}_{j,t}^i\ge T_s\\ 0& \left|{\underset{‾}{\overset{~}{W}}}_{j,t}^i\right|<T_s\\ {\underset{‾}{\overset{~}{W}}}_{j,t}^i+T_s& {\underset{‾}{\overset{~}{W}}}_{j,t}^i\le -T_s\end{array}\right)$

式中，T_s是门限阈值，10^-5＜T_s＜10^-3；表示在尺度j上经过门限阈值处 理后的第i个小波系数。

(4)对门限阈值处理后的光纤陀螺输出信号进行OMDWT小波重构

${\tilde{V}}_{j-1,t}=\underset{l=0}{\overset{L-1}{\Sigma }}\underset{i=1}{\overset{M-1}{\Sigma }}{\tilde{h}}_l^i{\tilde{W}}_{j,t+M^{j-1}l\mathrm{mod}N}^i+\underset{l=0}{\overset{L-1}{\Sigma }}{\tilde{g}}_l{\tilde{V}}_{j,t+M^{j-1}l\mathrm{mod}N}$

式中，表示t时刻第j-1尺度的尺度系数；和分别 表示第t+M^j-1lmodN时刻尺度j的尺度系数和门限阈值处理后的第i个小波系 数。

利用上式经过计算，最后在第J₀尺度上得到重构的去噪后的光纤陀螺信 号$\tilde{X}={\tilde{V}}_{J_0,t}.$

本发明与现有技术相比的优点在于：本发明在进行去噪时，按不同频带 分解从而速度较快，且对光纤陀螺信号有更细的频带划分；(2)光纤陀螺 信号循环移动后，其OMDWT小波系数和尺度系数会循环移动相应节拍， 有利于光纤陀螺的信号去噪。

附图说明

图1为本发明的计算流程图。

具体实施方式

如图1所示，本发明的实现方法如下：

1、建立光纤陀螺输出信号模型

其中，X_t为光纤陀螺输出信号；t＝0，...，N-1，N表示光纤陀螺输出信号的 长度；ω_ie为地球自转角速率，ω_ie＝7.27×10^-5rad/s；为测试点的地理纬度；ε_t为 陀螺漂移，陀螺漂移由常值分量、周期分量和白噪声组成，即

ε_t＝ε_d+Ω_dsin(2πf_d+θ₀)+W_t                     (2)

式中，ε_d为常值分量，短时间内近似为一个常数；Ω_d为周期分量的幅值； f_d为周期分量的频率；θ₀为初始相位；W_t为零均值高斯白噪声。式中的 Ω_dsin(2πf_d+θ₀)和W_t两项构成了光纤陀螺输出的随机漂移项。

2、对光纤陀螺输出信号作OMDWT小波分解

a.对光纤陀螺输出信号作OMDWT一步分解

令

$X_t=\left(\begin{array}{c}{X}_0\\ X_1\\ .\\ .\\ .\\ X_{N-3}\\ X_{N-2}\\ X_{N-1}\end{array}\right)$

定义

OMDWT算法的一步分解过程为：用小波系数分别对X_t，进行MDWT分解，得到和V₁，并按式(5)排 列。也就是说，进行M次MDWT分解后，把M组系数排列起来，即

其中，中的上标i表示此矩阵由M带小波变换的第i个小波系数组成。

例如，当M＝3，L＝6时

把式(6)中的小波滤波器系数替换为g_l，得到

和V₁的元素可表示为：

$W_1^i={[M^{1/M}{\tilde{W}}_{1,M-1}^i,M^{1/M}{\tilde{W}}_{\mathrm{1,2}M-1}^i,M^{1/M}{\tilde{W}}_{\mathrm{1,3}M-1}^i,...,M^{1/M}{\tilde{W}}_{1,N-1}^i]}^T---\left(8\right)$

$V_1={[M^{1/M}{\tilde{V}}_{1,M-1},M^{1/M}{\tilde{V}}_{\mathrm{1,2}M-1},M^{1/M}{\tilde{V}}_{\mathrm{1,3}M-1},...,M^{1/M}{\tilde{V}}_{1,N-1}]}^T---\left(9\right)$

可以看出，用小波滤波器和尺度滤波器{g_l}与时间序列X_t循环卷积， 得到和V₁。和V₁分别由长度为N的序列和的第 M-1，2M-1，...，N-1个元素组成。

把式(5)中的X_t替换为循环移动矢量和则有

定义

.

.

.

则

例如，当M＝3，L＝6时，有

把替换为g_l，可得到和的元素可表 示为：

的元素可表示为：

同样，的元素可表示为：

的元素可表示为：

定义

${\tilde{W}}_1^i\equiv {[{\tilde{W}}_{\mathrm{1,0}}^i,{\tilde{W}}_{\mathrm{1,1}}^i,{\tilde{W}}_{\mathrm{1,2}}^i,...,{\tilde{W}}_{1,N-1}^i]}^T---\left(21\right)$

${\tilde{V}}_1\equiv {[{\tilde{V}}_{\mathrm{1,0}},{\tilde{V}}_{\mathrm{1,1}},{\tilde{V}}_{\mathrm{1,2}},...,{\tilde{V}}_{1,N-1}]}^T---\left(22\right)$

将和的元素除以系数M^1/M，并将它们按式(21) 排列，可得到用类似方法由V₁，和得到。这里， 和的元素实际上是用OMDWT的滤波器和对X_t滤波后，得到的 输出和通过把矩阵和中的每个替换为 并将这些矩阵中的行按式(23)排列，得到N×N矩阵

例如，当M＝3，L＝6时

故

把式(23)中的替换为可得到则最后，把OMDWT 算法的一步分解表示为：

其中

通过这个过程，对光纤陀螺信号作OMDWT小波分解时，按不同频带 分解从而速度较快，且对光纤陀螺信号有更细的频带划分；光纤陀螺信号循 环移动后，其OMDWT小波系数和尺度系数会循环移动相应节拍。

b.求光纤陀螺信号X_t在第j尺度上的OMDWT变换

对步骤1得到的长度N的光纤陀螺信号X_t作OMDWT变换，将第j尺度 的OMDWT小波系数和尺度系数分别用N维矢量和表示，则对应的 元素和对应的元素分别为：

${\underset{‾}{\overset{~}{W}}}_{j,t}^i\equiv \underset{l=0}{\overset{L_j-1}{\Sigma }}{\tilde{h}}_{j,l}^i{X}_{t-l\mathrm{mod}N},$${\tilde{V}}_{j,t}\equiv \underset{l=0}{\overset{L_j-1}{\Sigma }}{\tilde{g}}_{j,l}{X}_{t-l\mathrm{mod}N},i=1,...,M-1---\left(27\right)$

式中，t＝0，...，N-1，N表示光纤陀螺信号的长度；l＝1，...，L， j＝J₀，J₀+1，...，J，J₀表示小波变换的初始尺度，J表示小波变换的终止尺度； 和分别为第j尺度上OMDWT的第i个小波滤波器系数和尺度滤波 器系数；lmodN表示l对N取余。

${\tilde{h}}_{j,l}^i\equiv h_{j,l}^i/M^{j/2},$${\tilde{g}}_{j,l}\equiv g_{j,l}/M^{j/2}---\left(28\right)$

和{g_j，l}分别表示第j尺度上MDWT的第i个小波滤波器系数和尺度 滤波器系数。可以推得，滤波器的长度是L_j≡(M^j-1)(L-1)+1。

利用光纤陀螺信号在第j-1尺度上的尺度系数矢量计算出光纤陀螺 信号在第j尺度上的小波系数矢量和尺度系数矢量的元素。和 对应的元素分别用和表示，则：

${\underset{‾}{\overset{~}{W}}}_{j,t}^i\underset{l=0}{\overset{L-1}{\Sigma }}{\tilde{h}}_l^i{\tilde{V}}_{j-1,t-M^{j-1}\mathrm{mod}N},t=\mathrm{0,1},...,N-1---\left(29\right)$

${\tilde{V}}_{j,t}=\underset{l=0}{\overset{L-1}{\Sigma }}{\tilde{g}}_l{\tilde{V}}_{j-1,t-M^{j-1}\mathrm{mod}N},t=\mathrm{0,1},...,N-1---\left(30\right)$

式中，和分别表示OMDWT的第i个小波滤波器系数和尺度滤 波器系数；表示t-M^j-1lmodN时刻第j-1尺度的尺度系数， j＝J₀，J₀+1，...，J。

c.求OMDWT小波滤波器系数和尺度滤波器系数

${\tilde{h}}_l^i\equiv h_l^i/\sqrt{M},$${\tilde{g}}_l\equiv g_l/\sqrt{M}---\left(31\right)$

式中，l＝1，...，L；i＝1，...，M-1；和{g_l}分别为MDWT的第i个小波滤波 器和尺度滤波器系数。

3、对OMDWT小波系数作门限阈值处理

对步骤2中光纤陀螺信号分解后得到的小波系数作门限阈值处理。 利用门限阈值处理，将绝对值小于门限阈值的小波系数置零，但对绝对值大 于门限阈值的小波系数作了收缩。

${\tilde{W}}_{j,t}^i({\underset{‾}{\overset{~}{W}}}_{j,t}^i,T_s)=\left(\begin{array}{cc}{\underset{‾}{\overset{~}{W}}}_{j,t}^i-T_s,& {\underset{‾}{\overset{~}{W}}}_{j,t}^i\ge T_s\\ 0& \left|{\underset{‾}{\overset{~}{W}}}_{j,t}^i\right|<T_s\\ {\underset{‾}{\overset{~}{W}}}_{j,t}^i+T_s& {\underset{‾}{\overset{~}{W}}}_{j,t}^i\le -T_s\end{array}\right),t=\mathrm{0,1},...,N-1---\left(32\right)$

式中，T_s是门限阈值，10^-5＜T_s＜10^-3；表示在尺度j上经过门限阈值处 理后的第i个小波系数。

4、对门限阈值处理之后的光纤陀螺输出信号进行OMDWT小波重构 光纤陀螺输出信号的OMDWT重构算法用下式描述：

${\tilde{V}}_{j-1,t}^i=\underset{l=0}{\overset{L-1}{\Sigma }}\underset{i=1}{\overset{M-1}{\Sigma }}{\tilde{h}}_l^i{\tilde{W}}_{j,t+M^{j-1}l\mathrm{mod}N}^i+\underset{l=0}{\overset{L-1}{\Sigma }}{\tilde{g}}_l{\tilde{V}}_{j,t+M^{j-1}l\mathrm{mod}N}---\left(33\right)$

式中，t＝0，1，...，N-1；和分别表示第t+M^j-1lmodN时刻 尺度j的尺度系数和门限阈值处理后的第i个小波系数。

利用上式经过计算，最后在第J₀尺度上得到重构的去噪后的光纤陀螺信 号$\tilde{X}={\tilde{V}}_{J_0,t}.$

本发明基于OMDWT对光纤陀螺信号进行去噪，与基于2带小波的光 纤陀螺信号去噪方法相比，计算速度快，噪声压缩比大。

本发明说明书中未作详细描述的内容属于本领域专业技术人员公知的 现有技术。