CN103557856A - 一种光纤陀螺随机漂移实时滤波方法 - Google Patents

Abstract

本发明公开了一种光纤陀螺随机漂移实时滤波方法，与现有技术相比，本发明在船用捷联航姿系统中，对采集到的光纤陀螺信号进行实时降噪处理：利用第二代小波的提升算法将滑动数据窗中的数据进行指定层数的分解，得到各层的小波系数和最后一层的尺度系数；对分解后的各层小波系数，建立相应的阈值规则，对其进行阈值量化处理；将处理后的小波系数，结合最后一层的尺度系数，逐级重构各层尺度系数，得到降噪处理后的信号；本发明有效提高了光纤陀螺信号实时降噪的处理精度和反应速度，继而抑制了姿态信息误差；在应用过程中，滑动数据窗宽度、小波分解层数等参数的设置较为简便，可根据仿真实验得到的经验值进行设定。

Description

一种光纤陀螺随机漂移实时滤波方法

技术领域

本发明涉及一种信号的消噪方法，尤其涉及一种应用于捷联航姿系统中的光纤陀螺随机漂移实时滤波方法。

背景技术

光纤陀螺其是以Sagnac效应为基础而发展起来的新型全固态陀螺仪,作为主要惯性敏感元件之一，由于其具有潜在的最好的性能、成本比而倍受青睐,已成为新一代中低精度捷联惯性导航、制导系统中的理想惯性器件。但是，由于原理、加工与装配工艺不完善等造成陀螺漂移，成为系统误差的主要来源，严重影响着惯导系统的初始对准精度及导航精度。

光纤陀螺的漂移包含两种成分：一种是具有确定性质的、有规律的、有系统性的漂移；另一种是由非确定性的随机性质干扰引起的随时间变化的随机漂移。对应前者，可以对这种确定性漂移进行建模,通过数学模型的软件补偿进行消除，这是比较容易做到的。而对于后者，目前采用较多的处理技术包括两种：一、对光纤陀螺的随机漂移进行时间序列建模，依据模型设计卡尔曼滤波器，通过卡尔曼滤波对随机漂移进行抑制；二、采用数字滤波器对光纤陀螺的漂移进行抑制，通过滤波器对光纤陀螺的漂移进行抑制。但实际情况中，由于光纤陀螺的输出受包括探测器噪声、光路噪声、电路噪声和环境变化引起的大量噪声影响,且这些噪声是时变的,故事先不能得到准确的误差统计特性，难以得到准确的漂移模型。因此，利用第一种方法无法得到理想的抑制效果。另外，由于随机漂移的频域特性不明确，增加了数字滤波器的设计难度，第二种方法也有其固有的缺陷。

针对光纤陀螺信号中的随机漂移和噪声，从抑制噪声的角度出发，消除噪声,对光纤陀螺的随机漂移也有抑制和补偿作用，是比较可行的一种方法。在实际系统中，为满足降噪的实时性要求，通常需要降噪的方法具有在线处理能力。因此，如何设计一种既优化了降噪效果，又提高了算法运算速度、简化了降噪处理过程的实时消噪方法，对光纤陀螺的不确定性随机漂移和噪声进行实时直接滤除，对捷联航姿系统精度的提高起着至关重要的作用。

发明内容

本发明的目的是要提供一种应用于捷联航姿系统中的光纤陀螺随机漂移实时滤波方法。

为达到上述目的，本发明是按照以下技术方案实施的：

本发明包括以下步骤：

步骤1：确定滑动数据窗宽度及预测器系数P和更新器系数U的维数，并离线求得系数P和系数U；

步骤2：捷联航姿系统预热结束后，采集光纤陀螺的输出数据，当采集数据量少于滑动数据窗宽度时，直接将之作为当前姿态解算的输入角速度信息ω_ib，解算后输出姿态信息（即开始采集数据后的小于第一个滑动数据窗宽度时间内的数据不作处理）；

步骤3：待采集陀螺数据量达到窗口宽度时，将第一段采集到的数据输入至去噪算法程序，进行降噪处理，在下一个采样点到来之前，输出降噪处理后的最后一个数据值，作为当前时刻陀螺信号的更新值；

步骤4：将更新后的当前采样点作为姿态实时解算的输入角速度信息ω_ib，进行姿态的实时结算，输出当前姿态信息；

步骤5：进入下一个采样周期，滑动数据窗实时地平移一个点，总包含最新的采样点作为窗口的最后一个值，重复步骤3、步骤4和步骤5，完成捷联航姿系统中光纤陀螺输出信号的实时降噪；

本发明所述步骤1中的离线求得预测器系数P、更新器系数U满足如下等式：

V^(M,M)P^(M,1)＝δ^(M,1)

W^(N,K)·Q^(K,N)·U^(N,1)＝δ^(N,1)

其中M和N为偶数，M_d＝M/2-1,N_d＝N/2-1,且有：

δ^(M,1)＝[100…0]^T

M₀＝2M_d+1,N₀＝2(M_d+N_d+1)

m＝0,1,…,M-1,n＝-M₀,-M₀+2,…,-M₀+2(M-1):V^(M,M)(m+1,(n+M₀)/2+1)＝n^mi＝1,2,…,2M-1:q(i)＝0;q(2+2M_d)＝1

i＝1,2,…,2M:q(2i-1)＝p(i)

L＝length(q)

K＝2(M-1)+L

Q^(K,N)＝0;(n＝1,2,…,N:Q(2n-1:2n-1+L-1,n)＝q)

m＝0,1,…,N-1,n＝-N₀,-N₀+1,…,-N₀+(K-1):W^(N,K)(m+1,n+N₀+1)＝n^m

本发明所述步骤3中的去噪算法程序，其特征为：

建立以第二代小波变换的阈值去噪法为基础的降噪算法，利用提升算法对原始信号进行小波分解，分为分裂、预测和更新：

$\left\{\begin{array}{cc}{x}_e\left[n\right]=x\left[2n\right]& n=\mathrm{1,2}\·\·\·\\ x_o\left[n\right]=x[2n-1]& n=\mathrm{1,2},\·\·\·\end{array}\right.$

d[n]＝x_o[n]-P(x_e[n])

c[n]＝x_e[n]+U(d[n])

其中原始信号为x[n]，x_o[n]为奇样本，x_e[n]为偶样本，P和U分别为离线求得预测算子和更新算子，d[n]和c[n]分别为得到的小波和尺度系数。进行逐级分解至L层得到各层细节（小波系数d_l,d_l-1,…,d₁）和最后一层概貌部分（尺度系数^c _l）。

x[n]＝{c_l,d_l,d_l-1,…,d₁}

建立阈值降噪法的改进阈值处理策略，对分解后的各层小波系数进行阈值量化处理，得到处理后的各层小波系数及未作处理的最后一层尺度系数。

采用提升方案的第二代小波变换对信号的逐级重构，由最后一层开始的，分为：恢复更新、恢复预测和合并：

$c_{l-1}^e\left[n\right]=c_l\left[n\right]-U\left(d_l\right[n\left]\right)$

$c_{l-1}^o\left[n\right]=d_l\left[n\right]+P\left(c_{l-1}^e\right[n\left]\right)$

$\left\{\begin{array}{cc}{c}_{l-1}\left[2n\right]=c_{l-1}^e\left[n\right]& n=\mathrm{1,2}\·\·\·\\ c_{l-1}[2n-1]=c_{l-1}^o\left[n\right]& n=\mathrm{1,2},\·\·\·\end{array}\right.$

其中c_l[n]、d_l[n]分别为最后一层尺度系数和经阈值量化处理后的小波系数，和

分别为l-1层尺度系数c_l-1[n]的偶序列和奇序列。逐级重构后得到复原的降噪处理后的信号x'[n]。

本发明所述建立阈值降噪法的改进阈值处理策略包括：

确定每层阈值：采用通用阈值法，第j层阈值T_j由下式决定：

$T_j={\mathrm{\σ}}_j\sqrt{2\ln \left(2^j\right)}$

其中，σ_j为第j层信号噪声标准差，该值由下式决定：

${\mathrm{\σ}}_j=\frac{\mathrm{median}\left(d_j\right)}{0.6745}$

其中，median(d_j)为第j层小波系数d_j的中值。

确定阈值函数：采用硬阈值函数：

$d_{j,k}=\left\{\begin{array}{cc}\left|d_{j,k}\right|,& \left|d_{j,k}\right|>T_j\\ 0,& |{\text{d}}_{j,k}|\≤T_j\end{array}\right.$

建立改进的阈值处理策略：将分解后的第一层小波系数强制置零，而后分解的小波系数仍按照硬阈值函数的策略进行处理。

$d_{j,k}\text{=}\left\{\begin{array}{cc}0,& j=1\\ \left|d_{\text{j,k}}\right|,& \left|d_{j,k}\right|>T_j,j>1\\ 0,& \left|d_{j,k}\right|\≤T_j,j>1\end{array}\right.$

本发明的有益效果是：

本发明是一种光纤陀螺随机漂移实时滤波方法，与现有技术相比，本发明在船用捷联航姿系统中，对采集到的光纤陀螺信号进行实时降噪处理：利用第二代小波的提升算法将滑动数据窗中的数据进行指定层数的分解，得到各层的小波系数（细节信号）和最后一层的尺度系数（概貌信号）；对分解后的各层小波系数，建立相应的阈值规则，对其进行阈值量化处理；将处理后的小波系数，结合最后一层的尺度系数，逐级重构各层尺度系数，得到降噪处理后的信号；更新滑动数据窗最前端的陀螺采样点，将之作为姿态解算的输入陀螺信息，进行实时姿态解算输出；进入下一个采样周期，窗口平移一个采样点，继续进行上述降噪处理。本发明有效提高了光纤陀螺信号实时降噪的处理精度和反应速度，继而抑制了姿态信息误差；在应用过程中，滑动数据窗宽度、小波分解层数等参数的设置较为简便，可根据仿真实验得到的经验值进行设定。

附图说明

图1是本发明的第二代小波分解、重构示意图；

图2是本发明的第二代小波阈值法实时降噪流程图；

图3是本发明的正弦仿真信号降噪效果比较图；

图4是本发明的滑动数据窗宽度的影响图；

图5是本发明的陀螺信号降噪时域图；

图6是本发明的陀螺信号降噪频域图；

图7是本发明的姿态解算输出图。

具体实施方式

下面结合附图以及具体实施例对本发明作进一步描述，在此发明的示意性实施例以及说明用来解释本发明，但并不作为对本发明的限定。

如图1所示：

本发明的技术方案不限于上述具体实施例的限制，凡是根据本发明的技术方案做出的技术变形，均落入本发明的保护范围之内。

本发明包括以下步骤：

步骤1：确定滑动数据窗宽度及预测器系数P和更新器系数U的维数，并离线求得系数P和系数U；

P^(M,1)＝[p₁ p₂ p₃…p_M]^T

U^(N,1)＝[u₁ u₂ u₃…u_N]^T

其中P和U满足如下等式：

V^(M,M)P^(M,1)＝δ^(M,1)

W^(N,K)·Q^(K,N)·U^(N,1)＝δ^(N,1)

其中，M和N为偶数，M_d＝M/2-1,N_d＝N/2-1,其中：

δ^(M,1)＝[100…0]^T

M₀＝2M_d+1,N₀＝2(M_d+N_d+1)

m＝0,1,…,M-1,n＝-M₀,-M₀+2,…,-M₀+2(M-1):V^(M,M)(m+1,(n+M₀)2+1)＝n^m

i＝1,2,…,2M-1:q(i)＝0;q(2+2M_d)＝1

i＝1,2,…,2M:q(2i-1)＝p(i)

L＝length(q)

K＝2(M-1)+L

Q^(K,N)＝0;(n＝1,2,…,N:Q(2n-1:2n-1+L-1,n)＝q)

m＝0,1,…,N-1,n＝-N₀,-N₀+1,…,-N₀+(K-1):W^(N,K)(m+1,n+N₀+1)＝n^m

步骤2：捷联航姿系统预热结束后，采集光纤陀螺的输出数据，当采集数据量少于滑动数据窗宽度时，直接将之作为当前姿态解算的输入角速度信息ω_ib，解算后输出姿态信息（即开始采集数据后的小于第一个滑动数据窗宽度时间内的数据不作处理）；

步骤3：待采集陀螺数据量达到窗口宽度时，将第一段采集到的数据输入至去噪算法程序，进行降噪处理，在下一个采样点到来之前，输出降噪处理后的最后一个数据值，作为当前时刻陀螺信号的更新值；

1、采用提升方案的第二代小波变换对信号进行多尺度分解：具体步骤为（1）分裂、（2）预测、（3）更新。（1）中，分解层数越多，进行降噪处理的频带越宽，但同时可能会影响陀螺的动态输出特性，且增加了算法的复杂程度，故分解层数不宜超过5层。（2）中，随着信号的逐级分解，其对小波系数的阈值量化处理也随之进行，如图2所示。

（1）分裂：将原始的离散信号分解为两组相互关联的部分，通常分解为奇样本x_o[n]和偶样本x_e[n]：

$\left\{\begin{array}{cc}{x}_e\left[n\right]=x\left[2n\right]& n=\mathrm{1,2}\·\·\·\\ x_o\left[n\right]=x[2n-1]& n=\mathrm{1,2},\·\·\·\end{array}\right.$

（2）预测：构造预测算子P与x_e[n]作用，得到小波系数d[n]，作为由x_e[n]预测x_o[n]的误差：

d[n]＝x_o[n]-P(x_e[n])

或 $d\left(n\right)=x_o\left(n\right)-\underset{k=1}{\overset{M}{\mathrm{\Σ}}}{p}_K{x}_e(n-M_d+(k-1))$

（3）更新：构造更新算子U与预测的小波系数d[n]作用，并叠加到x_e[n]上，得到尺度系数c(n)，c(n)代表了原信号的粗略近似：

c[n]＝x_e[n]+U(d[n])

或 $c\left(n\right)=x_e\left(n\right)+\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{u}_{K}d(n-N_d-1+(k-1))$

至此完成了第二代小波的第一层分解，重复步骤1）至3）可以继续进行分解，不同的是重复步骤1）时，原信号x[n]由上一层尺度系数c_j(n)代替之。继续分裂、预测和更新L次后，有

x[n]＝{c_l,d_l,d_l-1,…,d₁}

其中，c_l代表了信号的低频概貌部分（尺度系数），其余为信号的高频细节部分（小波系数）。

2、建立阈值降噪法的改进阈值处理策略，对分解后的各层小波系数进行阈值量化处理。具体步骤分为：（1）确定每层阈值、（2）确定阈值函数、（3）建立改进的阈值处理策略：

（1）确定每层阈值：采用通用阈值法，第j层阈值T_j由下式决定：

$T_j={\mathrm{\σ}}_j\sqrt{2\ln \left(2^j\right)}$

其中，σ_j为第j层信号噪声标准差，该值由下式决定：

${\mathrm{\σ}}_j=\frac{\mathrm{median}\left(d_j\right)}{0.6745}$

其中，median(d_j)为第j层小波系数d_j的中值。

（2）确定阈值函数：在选取基于二代小波变换的阈值函数时，应用硬阈值不仅与软阈值降噪效果相差不大，且其拥有更为简洁的处理策略，故选择硬阈值函数：

$d_{j,k}=\left\{\begin{array}{cc}\left|d_{j,k}\right|,& \left|d_{j,k}\right|>T_j\\ 0,& |{\text{d}}_{j,k}|\≤T_j\end{array}\right.$

（3）建立改进的阈值处理策略：改进的阈值处理策略结合了强制降噪对信号处理反应速度的提升，即将步骤3中的（1）处理分解后的第一层小波系数强制置零，而后分解的小波系数仍按照硬阈值函数的策略进行处理。

3、采用提升方案的第二代小波变换对信号的逐级重构，这种重构是由最后一层开始的。具体步骤分为：（1）恢复更新、（2）恢复预测、（3）合并。

（1）恢复更新：利用最后一层尺度系数c_l[n]和经阈值量化处理后小波系数d_l[n]，可恢复出l-1层尺度系数c_l-1[n]的偶序列：

$c_{l-1}^e\left[n\right]=c_l\left[n\right]-U\left(d_l\right[n\left]\right)$

或 $c_{l-1}^e=c_l\left(n\right)-\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{u}_K{d}_l(n-N_d-1+(k-1))$

（2）恢复预测：由1）得出的和进行阈值量化处理后的d_l[n]预测出l-1层尺度系数c_l-1[n]的奇序列：

$c_{l-1}^o\left[n\right]=d_l\left[n\right]+P\left(c_{l-1}^e\right[n\left]\right)$

或 $c_{l-1}^o\left(n\right)=d_l\left(n\right)+\underset{k=1}{\overset{M}{\mathrm{\Σ}}}{p}_K{c}_{l-1}^e(n-M_d+(k-1))$

（3）合并：由1）和2）获得的c_l-1[n]的奇序列和偶序列，合并恢复出j-1层尺度信号c_l-1[n]：

$\left\{\begin{array}{cc}{c}_{l-1}\left[2n\right]=c_{l-1}^e\left[n\right]& n=\mathrm{1,2}\·\·\·\\ c_{l-1}[2n-1]=c_{l-1}^o\left[n\right]& n=\mathrm{1,2},\·\·\·\end{array}\right.$

进行L次的恢复更新、恢复预测、合并之后，可以复原出原始信号x[n]。

4、更新滑动数据窗内的末端数据，即将降噪处理后返回的此段信号的最后一个采样点的值更新当前采样点，完成对当前陀螺数据的实时更新。

步骤4：将更新后的当前采样点作为姿态实时解算的输入角速度信息ω_ib，进行姿态的实时结算，输出当前姿态信息；

步骤5：进入下一个采样周期，滑动数据窗实时地平移一个点，总包含最新的采样点作为窗口的最后一个值，重复步骤3、步骤4和步骤5，完成捷联航姿系统中光纤陀螺输出信号的实时降噪；

本发明所述步骤1中的离线求得预测器系数P、更新器系数U满足如下等式：

V^(M,M)P^(M,1)＝δ^(M,1)

W^(N,K)·Q^(K,N)·U^(N,1)＝δ^(N,1)

其中M和N为偶数，M_d＝M/2-1,N_d＝N/2-1,且有：

δ^(M,1)＝[100…0]^T

M₀＝2M_d+1,N₀＝2(M_d+N_d+1)

m＝0,1,…,M-1,n＝-M₀,-M₀+2,…,-M₀+2(M-1):V^(M,M)(m+1,(n+M₀)/2+1)＝n^mi＝1,2,…,2M-1:q(i)＝0;q(2+2M_d)＝1

i＝1,2,…,2M:q(2i-1)＝p(i)

L＝length(q)

K＝2(M-1)+L

Q^(K,N)＝0;(n＝1,2,…,N:Q(2n-1:2n-1+L-1,n)＝q)

m＝0,1,…,N-1,n＝-N₀,-N₀+1,…,-N₀+(K-1):W^(N,K)(m+1,n+N₀+1)＝n^m

本发明所述步骤3中的去噪算法程序，其特征为：

建立以第二代小波变换的阈值去噪法为基础的降噪算法，利用提升算法对原始信号进行小波分解，分为分裂、预测和更新：

$\left\{\begin{array}{cc}{x}_e\left[n\right]=x\left[2n\right]& n=\mathrm{1,2}\·\·\·\\ x_o\left[n\right]=x[2n-1]& n=\mathrm{1,2},\·\·\·\end{array}\right.$

d[n]＝x_o[n]-P(x_e[n])

c[n]＝x_e[n]+U(d[n])

其中原始信号为x[n]，x_o[n]为奇样本，x_e[n]为偶样本，P和U分别为离线求得预测算子和更新算子，d[n]和c[n]分别为得到的小波和尺度系数。进行逐级分解至L层得到各层细节（小波系数d_l,d_l-1,…,d₁）和最后一层概貌部分（尺度系数^cl）。

x[n]＝{c_l,d_l,d_l-1,…,d₁}

建立阈值降噪法的改进阈值处理策略，对分解后的各层小波系数进行阈值量化处理，得到处理后的各层小波系数及未作处理的最后一层尺度系数。

采用提升方案的第二代小波变换对信号的逐级重构，由最后一层开始的，分为：恢复更新、恢复预测和合并：

$c_{l-1}^e\left[n\right]=c_l\left[n\right]-U\left(d_l\right[n\left]\right)$

$c_{l-1}^o\left[n\right]=d_l\left[n\right]+P\left(c_{l-1}^e\right[n\left]\right)$

$\left\{\begin{array}{cc}{c}_{l-1}\left[2n\right]=c_{l-1}^e\left[n\right]& n=\mathrm{1,2}\·\·\·\\ c_{l-1}[2n-1]=c_{l-1}^o\left[n\right]& n=\mathrm{1,2},\·\·\·\end{array}\right.$

其中c_l[n]、d_l[n]分别为最后一层尺度系数和经阈值量化处理后的小波系数，

和

分别为l-1层尺度系数c_l-1[n]的偶序列和奇序列。逐级重构后得到复原的降噪处理后的信号x'[n]。

本发明所述建立阈值降噪法的改进阈值处理策略包括：

确定每层阈值：采用通用阈值法，第j层阈值T_j由下式决定：

$T_j={\mathrm{\σ}}_j\sqrt{2\ln \left(2^j\right)}$

其中，σ_j为第j层信号噪声标准差，该值由下式决定：

${\mathrm{\σ}}_j=\frac{\mathrm{median}\left(d_j\right)}{0.6745}$

其中，median(d_j)为第j层小波系数d_j的中值。

确定阈值函数：采用硬阈值函数：

$d_{j,k}=\left\{\begin{array}{cc}\left|d_{j,k}\right|,& \left|d_{j,k}\right|>T_j\\ 0,& |{\text{d}}_{j,k}|\≤T_j\end{array}\right.$

建立改进的阈值处理策略：将分解后的第一层小波系数强制置零，而后分解的小波系数仍按照硬阈值函数的策略进行处理。

$d_{j,k}\text{=}\left\{\begin{array}{cc}0,& j=1\\ \left|d_{\text{j,k}}\right|,& \left|d_{j,k}\right|>T_j,j>1\\ 0,& \left|d_{j,k}\right|\≤T_j,j>1\end{array}\right.$

如图1所示：第二代小波变换对信号进行分解

首先确定维数为M的预测器P和维数为N更新器U：

P^(M,1)＝[p₁ p₂ p₃…p_M]^T

U^(N,1)＝[u₁ u₂ u₃…u_N]^T

其可由如下方程求得：

V^(M,M)P^(M,1)＝δ^(M,1)

W^(N,K)·Q^(K,N)·U^(N,1)＝δ^(N,1)

对于偶数M和N，以及M_d＝M/2-1,N_d＝N/2-1,利用上式即可求得P和U，其中

δ^(M,1)＝[100…0]^T

M₀＝2M_d+1,N₀＝2(M_d+N_d+1)

m＝0,1,…,M-1,n＝-M₀,-M₀+2,…,-M₀+2(M-1):V^(M,M)(m+1,(n+M₀)2+1)＝n^m

i＝1,2,…,2M-1:q(i)＝0;q(2+2M_d)＝1

i＝1,2,…,2M:q(2i-1)＝p(i)

L＝length(q)

K＝2(M-1)+L

Q^(K,N)＝0;(n＝1,2,…,N:Q(2n-1:2n-1+L-1,n)＝q)

m＝0,1,…,N-1,n＝-N₀,-N₀+1,…,-N₀+(K-1):W^(N,K)(m+1,n+N₀+1)＝n^m

预测器系数P和更新器系数U确定后，采用提升方案的第二代小波变换实施信号的分解，过程如下：

步骤（1）：分裂

将原始的离散信号分解为两组相互关联的部分，通常分解为奇样本x_o[n]和偶样本x_e[n]：

$\left\{\begin{array}{cc}{x}_e\left[n\right]=x\left[2n\right]& n=\mathrm{1,2}\·\·\·\\ x_o\left[n\right]=x[2n-1]& n=\mathrm{1,2},\·\·\·\end{array}\right.$

步骤（2）：预测

构造预测算子P与x_e[n]作用，得到小波系数d[n]，作为由x_e[n]预测x_o[n]的误差：

d[n]＝x_o[n]-P(x_e[n])

或 $d\left(n\right)=x_o\left(n\right)-\underset{k=1}{\overset{M}{\mathrm{\Σ}}}{p}_K{x}_e(n-M_d+(k-1))$

步骤（3）：更新

构造更新算子U与预测的小波系数d[n]作用，并叠加到x_e[n]上，得到尺度系数c(n)，c(n)代表了原信号的粗略近似：

c[n]＝x_e[n]+U(d[n])

或 $c\left(n\right)=x_e\left(n\right)+\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{u}_{K}d(n-N_d-1+(k-1))$

至此完成了第二代小波的第一层分解，重复步骤（1）至（3）可以继续进行分解，不同的是重复步骤（1）时，原信号的位置由上一层尺度系数c_j(n)代替之。继续分裂、预测和更新L次后，有

x[n]＝{c_l,d_l,d_l-1,…,d₁}

其中，c_l代表了信号的低频概貌部分，其余为信号的高频细节部分。

确定每层阈值

采用通用阈值法确定每层阈值，通用阈值法是Donoho和Johnstone在高斯噪声模型下，应用多维独立正态变量决策理论得出的阈值。即当维数趋于无穷时,噪声系数的幅值大于阈值的概率趋于零，其中σ_j为第j层信号噪声标准差。

如图5和图6所示：对于陀螺信号来说,噪声的标准差是未知的,这就需要对噪声水平进行估计,噪声的标准差可以用各尺度细节小波系数绝对值的中值来进行估计：

${\mathrm{\σ}}_j=\frac{\mathrm{median}\left(d_j\right)}{0.6745}$

式中，d_j为第j层小波系数。

确定阈值函数

确定了各层的阈值之后，需选取一种合适的阈值量化方法，常用的有硬阈值法、软阈值法和兼顾这两种方法优点的半软阈值法。硬阈值法往往具有较大的方差，而软阈值法往往带来较大的偏差。虽然半软阈值法是以上两种方法的折中，不仅保留了较大的系数，而且具有连续性，但需要确定两个阈值，增加了运算复杂度。通过实验分析，传统的一代小波变换使用软阈值函数对光纤陀螺输出信号进行滤波，降噪效果要优于硬阈值函数。但在基于第二代小波变换的前提下，两者的滤波效果相差无几，且硬阈值函数具有较为简洁的处理手段，故采用硬阈值函数：

$d_{j,k}=\left\{\begin{array}{cc}\left|d_{j,k}\right|,& \left|d_{j,k}\right|>T_j\\ 0,& |{\text{d}}_{j,k}|\≤T_j\end{array}\right.$

建立强制降噪与硬阈值函数组合应用机制：

实际情况下，由舰载捷联惯导系统的实际应用环境，决定了光纤陀螺的动态特性频段要远低于装置在如飞机这种高速运动载体上的频段（特别是天向陀螺），根据这种有用信号多集中在中低频段的特征，可以将第二代小波变换后的第一层小波系数，也就是高频段信号，进行强制降噪处理，即将之置零。而继续分解的第二层、第三层等小波系数仍采用硬阈值函数进行降噪处理。这种小波系数处理策略继承了硬阈值函数与强制去噪的计算简便、快速的优点，且基本保留了原始信号中低频段的动态特性。

如图1所示：第二代小波变换重构信号，采用提升方案的第二代小波变换实施信号的重构：

步骤（1）：恢复更新

给定某层尺度系数c[n]和小波系数d[n]，可恢复出偶序列：

x_e[n]＝c[n]-U(d[n])

或 $x_e\left(n\right)=c\left(n\right)-\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{u}_{K}d(n-N_d-1+(k-1))$

步骤（2）：恢复预测

由步骤（1）得出的x_e(n)和给定的d[n]预测出奇序列：

x_o[n]＝d[n]+P(x_e[n])

或 $x_o\left(n\right)=d-\left(n\right)+\underset{k=1}{\overset{M}{\mathrm{\Σ}}}{p}_K{x}_e(n-M_d+(k-1))$

步骤（3）：合并

由步骤（1）和步骤（2）获得奇序列和偶序列，将其合并恢复出原信号（或上一层尺度信号）。

$\left\{\begin{array}{cc}x\left[2n\right]=x_e\left[n\right]& n=\mathrm{1,2}\·\·\·\\ x[2n-1]=x_o\left[n\right]& n=\mathrm{1,2},\·\·\·\end{array}\right.$

至此完成了第二代小波的一层信号的重构，重复步骤（1）至（3）可以继续进行重构，不同的是重复步骤（3）时，原信号的位置由上一层尺度系数c_j(n)代替之。继续进行分解层数次的恢复更新、恢复预测、合并之后，可以复原出原始信号x[n]。作为实时降噪的应用，此时只输出重构信号的最后一个值，即只更新最新的采样点。

至此，实时降噪处理完毕，可将实时输出的最新采样点作为姿态结算的输入陀螺信号，进行实时姿态结算。

如图3所示：本发明较传统基于一代小波的实时降噪方案，去噪后的模拟信号偏离期望的程度，即方差小于传统方案去噪后的结果，且去噪后的信噪比较传统方案提高约3db，较含噪信号提高约6.7db，采用文中提到的改进方案的实时降噪效果优于传统一代小波的实时降噪方案。

方案类型：染噪信号的方差为1.4992、信噪比为-3.0198；传统方案的方差为0.9260、信噪比为0.6642；改进方案的方差为0.7117、信噪比为37120.

如图4所示：在以下条件下，为比较滑动数据窗宽度的影响，利用该方案对一组染噪正弦信号进行模拟在线实时降噪实验：

正弦信号频率：0.1Hz；幅值：1；采样点数：10000；采样频率：100Hz；噪声特征：每个采样点叠加一个标准正态分布的随机数；小波分解层数：3层；滑动数据窗宽度分别为：128、512、1024；

结果如图4所示，可以看出，滑动数据窗对本发明的实时降噪效果影响很小。

在以下条件下，利用该方案对光纤陀螺输出信号进行模拟在线实时降噪实验：

光纤陀螺捷联航姿系统放置于单轴转台上，转台无基准；采集静态数据5小时，采样频率：98Hz；小波分解层数：3层；滑动数据窗宽度：1024；数据处理长度：10000；

利用本发明所述方法得到结果如图5、图6所示，利用Allan方差分析法对本发明处理后的信号误差系数进行分析，比较了基于一代小波的实时降噪方案（传统方案）与本发明采用的基于二代小波变换的实时降噪方案（改进方案）的去噪性能：

降噪方案：原始信号的方差为4.7625、零偏稳定性/(⁰)·h^-1为0.1147392、量化噪声/μrad为1.5247672、速率斜坡/((⁰)·h^-2)为0.2481738；

传统方案的方差为2.2447、零偏稳定性/(⁰)·h^-1为0.1028541、量化噪声/μrad为1.0404798、速率斜坡/((⁰)·h^-2)为0.1753006；

改进方案的方差为0.8212、零偏稳定性/(⁰)·h^-1为0.0848103、量化噪声/μrad为0.6125475、速率斜坡/((⁰)·h^-2)为0.1069661。

在以下条件下，利用该方案对光线陀螺信号进行模拟实时降噪后，进行姿态解算实验：

光纤陀螺捷联航姿系统放置于单轴转台上，转台无基准；采集静态数据5小时，采样频率：98Hz；小波分解层数：3层滑动数据窗宽度：1024；降噪起始位置：3小时处；降噪数据长度：10000（约100s）；

系统初始参数如下：

载体初始位置：北纬45.7796°，东经126.6705°；地球自转角速度：7.2921158×10^-5rads；重力加速度：9.78049ms²；地球半径：6378393.0m；

利用本发明方案从陀螺输出3小时处进行约100s的降噪处理，而后进行姿态解算，输出姿态如图7所示。虽然只进行了约100s的降噪处理，但该方案仍在降噪后减小纵横摇震荡幅值约7角秒。可以预见，实际应用时除开始阶段的窗口长度数据不做处理外，所有的陀螺原始数据均进行降噪处理，等效于增加降噪处理的时间，将进一步提高输出姿态的精度，抑制系统静态输出的姿态信息震荡。

Claims (4)

1.一种光纤陀螺随机漂移实时滤波方法，其特征在于，包括以下步骤：

步骤1：确定滑动数据窗宽度及预测器系数P和更新器系数U的维数，并离线求得系数P和系数U；

步骤2：捷联航姿系统预热结束后，采集光纤陀螺的输出数据，当采集数据量少于滑动数据窗宽度时，直接将之作为当前姿态解算的输入角速度信息ω_ib，解算后输出姿态信息（即开始采集数据后的小于第一个滑动数据窗宽度时间内的数据不作处理）；

步骤3：待采集陀螺数据量达到窗口宽度时，将第一段采集到的数据输入至去噪算法程序，进行降噪处理，在下一个采样点到来之前，输出降噪处理后的最后一个数据值，作为当前时刻陀螺信号的更新值；

步骤4：将更新后的当前采样点作为姿态实时解算的输入角速度信息ω_ib，进行姿态的实时结算，输出当前姿态信息；

步骤5：进入下一个采样周期，滑动数据窗实时地平移一个点，总包含最新的采样点作为窗口的最后一个值，重复步骤3、步骤4和步骤5，完成捷联航姿系统中光纤陀螺输出信号的实时降噪。

2.根据权利权利要求1所述的光纤陀螺随机漂移实时滤波方法，其特征在于：所述步骤1中的离线求得预测器系数P、更新器系数U满足如下等式：

V^(M,M)P^(M,1)＝δ^(M,1)

W^(N,K)·Q^(K,N)·U^(N,1)＝δ^(N,1)

其中M和N为偶数，M_d＝M/2-1,N_d＝N/2-1,且有：

δ^(M,1)＝[100…0]^T

M₀＝2M_d+1,N₀＝2(M_d+N_d+1)

m＝0,1,…,M-1,n＝-M₀,-M₀+2,…,-M₀+2(M-1):V^(M,M)(m+1,(n+M₀)/2+1)＝n^m

i＝1,2,…,2M-1:q(i)＝0;q(2+2M_d)＝1

i＝1,2,…,2M:q(2i-1)＝p(i)

L＝length(q)

K＝2(M-1)+L

Q^(K,N)＝0;(n＝1,2,…,N:Q(2n-1:2n-1+L-1,n)＝q)

m＝0,1,…,N-1,n＝-N₀,-N₀+1,…,-N₀+(K-1):W^(N,K)(m+1,n+N₀+1)＝n^m。

3.根据权利权利要求1所述的光纤陀螺随机漂移实时滤波方法，其特征在于：所述步骤3中的去噪算法程序，其特征为：

建立以第二代小波变换的阈值去噪法为基础的降噪算法，利用提升算法对原始信号进行小波分解，分为分裂、预测和更新：

$\left\{\begin{array}{cc}{x}_e\left[n\right]=x\left[2n\right]& n=\mathrm{1,2}\·\·\·\\ x_o\left[n\right]=x[2n-1]& n=\mathrm{1,2},\·\·\·\end{array}\right.$

d[n]＝x_o[n]-P(x_e[n])

c[n]＝x_e[n]+U(d[n])

其中原始信号为x[n]，x_o[n]为奇样本，x_e[n]为偶样本，P和U分别为离线求得预测算子和更新算子，d[n]和c[n]分别为得到的小波和尺度系数。进行逐级分解至L层得到各层细节（小波系数d_l,d_l-1,…,d₁）和最后一层概貌部分（尺度系数c_l）；

x[n]＝{c_l,d_l,d_l-1,…,d₁}

建立阈值降噪法的改进阈值处理策略，对分解后的各层小波系数进行阈值量化处理，得到处理后的各层小波系数及未作处理的最后一层尺度系数。

采用提升方案的第二代小波变换对信号的逐级重构，由最后一层开始的，分为：恢复更新、恢复预测和合并：

$c_{l-1}^e\left[n\right]=c_l\left[n\right]-U\left(d_l\right[n\left]\right)$

$c_{l-1}^o\left[n\right]=d_l\left[n\right]+P\left(c_{l-1}^e\right[n\left]\right)$

$\left\{\begin{array}{cc}{c}_{l-1}\left[2n\right]=c_{l-1}^e\left[n\right]& n=\mathrm{1,2}\·\·\·\\ c_{l-1}[2n-1]=c_{l-1}^o\left[n\right]& n=\mathrm{1,2},\·\·\·\end{array}\right.$

其中^c _l[n]、d_l[n]分别为最后一层尺度系数和经阈值量化处理后的小波系数，

和

分别为l-1层尺度系数c_l-1[n]的偶序列和奇序列；逐级重构后得到复原的降噪处理后的信号x'[n]。

4.根据权利权利要求3所述的光纤陀螺随机漂移实时滤波方法，其特征在于：所述建立阈值降噪法的改进阈值处理策略包括：

确定每层阈值：采用通用阈值法，第j层阈值T_j由下式决定：

$T_j={\mathrm{\σ}}_j\sqrt{2\ln \left(2^j\right)}$

其中，σ_j为第j层信号噪声标准差，该值由下式决定：

${\mathrm{\σ}}_j=\frac{\mathrm{median}\left(d_j\right)}{0.6745}$

其中，median(d_j)为第j层小波系数d_j的中值；

确定阈值函数：采用硬阈值函数：

$d_{j,k}=\left\{\begin{array}{cc}\left|d_{j,k}\right|,& \left|d_{j,k}\right|>T_j\\ 0,& |{\text{d}}_{j,k}|\≤T_j\end{array}\right.$

建立改进的阈值处理策略：将分解后的第一层小波系数强制置零，而后分解的小波系数仍按照硬阈值函数的策略进行处理；

$d_{j,k}\text{=}\left\{\begin{array}{cc}0,& j=1\\ \left|d_{\text{j,k}}\right|,& \left|d_{j,k}\right|>T_j,j>1\\ 0,& \left|d_{j,k}\right|\≤T_j,j>1\end{array}\right..$

Images