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

摘要
申请专利号：	CN201310508747.3	申请日：	2013.10.25
公开号：	CN103557856A	公开日：	2014.02.05
当前法律状态：	撤回	有效性：	无权
法律详情：	发明专利申请公布后的视为撤回IPC(主分类):G01C 19/64申请公布日:20140205\|\|\|实质审查的生效IPC(主分类):G01C 19/64申请日:20131025\|\|\|公开
IPC分类号：	G01C19/64	主分类号：	G01C19/64
申请人：	哈尔滨工程大学
发明人：	孙枫; 夏健钟; 奔粤阳; 李敬春; 杨晓龙; 阮双双; 赵维珩
地址：	150001 黑龙江省哈尔滨市南岗区南通大街145号
优先权：
专利代理机构：	北京科亿知识产权代理事务所(普通合伙) 11350	代理人：	汤东凤
PDF下载：	PDF下载

内容摘要

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

权利要求书

权利要求书
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为偶数，Md＝M/2-1,Nd＝N/2-1,且有：
δ(M,1)＝[100…0]T
M0＝2Md+1,N0＝2(Md+Nd+1)
m＝0,1,…,M-1,n＝-M0,-M0+2,…,-M0+2(M-1):V(M,M)(m+1,(n+M0)/2+1)＝nm
i＝1,2,…,2M-1:q(i)＝0;q(2+2Md)＝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＝-N0,-N0+1,…,-N0+(K-1):W(N,K)(m+1,n+N0+1)＝nm。

3.  根据权利权利要求1所述的光纤陀螺随机漂移实时滤波方法，其特征在于：所述步骤3中的去噪算法程序，其特征为：
建立以第二代小波变换的阈值去噪法为基础的降噪算法，利用提升算法对原始信号进行小波分解，分为分裂、预测和更新：
xe[n]=x[2n]n=1,2···xo[n]=x[2n-1]n=1,2,···]]>
d[n]＝xo[n]-P(xe[n])
c[n]＝xe[n]+U(d[n])
其中原始信号为x[n]，xo[n]为奇样本，xe[n]为偶样本，P和U分别为离线求得预测算子和更新算子，d[n]和c[n]分别为得到的小波和尺度系数。进行逐级分解至L层得到各层细节（小波系数dl,dl-1,…,d1）和最后一层概貌部分（尺度系数cl）；
x[n]＝{cl,dl,dl-1,…,d1}
建立阈值降噪法的改进阈值处理策略，对分解后的各层小波系数进行阈值量化处理，得到处理后的各层小波系数及未作处理的最后一层尺度系数。
采用提升方案的第二代小波变换对信号的逐级重构，由最后一层开始的，分为：恢复更新、恢复预测和合并：
cl-1e[n]=cl[n]-U(dl[n])]]>
cl-1o[n]=dl[n]+P(cl-1e[n])]]>
cl-1[2n]=cl-1e[n]n=1,2···cl-1[2n-1]=cl-1o[n]n=1,2,···]]>
其中cl[n]、dl[n]分别为最后一层尺度系数和经阈值量化处理后的小波系数，和分别为l-1层尺度系数cl-1[n]的偶序列和奇序列；逐级重构后得到复原的降噪处理后的信号x'[n]。

4.  根据权利权利要求3所述的光纤陀螺随机漂移实时滤波方法，其特征在于：所述建立阈值降噪法的改进阈值处理策略包括：
确定每层阈值：采用通用阈值法，第j层阈值Tj由下式决定：
Tj=σj2ln(2j)]]>
其中，σj为第j层信号噪声标准差，该值由下式决定：
σj=median(dj)0.6745]]>
其中，median(dj)为第j层小波系数dj的中值；
确定阈值函数：采用硬阈值函数：
dj,k=|dj,k|,|dj,k|>Tj0,|dj,k|≤Tj]]>
建立改进的阈值处理策略：将分解后的第一层小波系数强制置零，而后分解的小波系数仍按照硬阈值函数的策略进行处理；
dj,k=0,j=1|dj,k|,|dj,k|>Tj,j>10,|dj,k|≤Tj,j>1.]]>

说明书

说明书一种光纤陀螺随机漂移实时滤波方法
技术领域
本发明涉及一种信号的消噪方法，尤其涉及一种应用于捷联航姿系统中的光纤陀螺随机漂移实时滤波方法。
背景技术
光纤陀螺其是以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为偶数，Md＝M/2-1,Nd＝N/2-1,且有：
δ(M,1)＝[100…0]T
M0＝2Md+1,N0＝2(Md+Nd+1)
m＝0,1,…,M-1,n＝-M0,-M0+2,…,-M0+2(M-1):V(M,M)(m+1,(n+M0)/2+1)＝nmi＝1,2,…,2M-1:q(i)＝0;q(2+2Md)＝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＝-N0,-N0+1,…,-N0+(K-1):W(N,K)(m+1,n+N0+1)＝nm
本发明所述步骤3中的去噪算法程序，其特征为：
建立以第二代小波变换的阈值去噪法为基础的降噪算法，利用提升算法对原始信号进行小波分解，分为分裂、预测和更新：
xe[n]=x[2n]n=1,2···xo[n]=x[2n-1]n=1,2,···]]>
d[n]＝xo[n]-P(xe[n])
c[n]＝xe[n]+U(d[n])
其中原始信号为x[n]，xo[n]为奇样本，xe[n]为偶样本，P和U分别为离线求得预测算子和更新算子，d[n]和c[n]分别为得到的小波和尺度系数。进行逐级分解至L层得到各层细节（小波系数dl,dl-1,…,d1）和最后一层概貌部分（尺度系数cl）。
x[n]＝{cl,dl,dl-1,…,d1}
建立阈值降噪法的改进阈值处理策略，对分解后的各层小波系数进行阈值量化处理，得到处理后的各层小波系数及未作处理的最后一层尺度系数。
采用提升方案的第二代小波变换对信号的逐级重构，由最后一层开始的，分为：恢复更新、恢复预测和合并：
cl-1e[n]=cl[n]-U(dl[n])]]>
cl-1o[n]=dl[n]+P(cl-1e[n])]]>
cl-1[2n]=cl-1e[n]n=1,2···cl-1[2n-1]=cl-1o[n]n=1,2,···]]>
其中cl[n]、dl[n]分别为最后一层尺度系数和经阈值量化处理后的小波系数，和分别为l-1层尺度系数cl-1[n]的偶序列和奇序列。逐级重构后得到复原的降噪处理后的信号x'[n]。
本发明所述建立阈值降噪法的改进阈值处理策略包括：
确定每层阈值：采用通用阈值法，第j层阈值Tj由下式决定：
Tj=σj2ln(2j)]]>
其中，σj为第j层信号噪声标准差，该值由下式决定：
σj=median(dj)0.6745]]>
其中，median(dj)为第j层小波系数dj的中值。
确定阈值函数：采用硬阈值函数：
dj,k=|dj,k|,|dj,k|>Tj0,|dj,k|≤Tj]]>
建立改进的阈值处理策略：将分解后的第一层小波系数强制置零，而后分解的小波系数仍按照硬阈值函数的策略进行处理。
dj,k=0,j=1|dj,k|,|dj,k|>Tj,j>10,|dj,k|≤Tj,j>1]]>
本发明的有益效果是：
本发明是一种光纤陀螺随机漂移实时滤波方法，与现有技术相比，本发明在船用捷联航姿系统中，对采集到的光纤陀螺信号进行实时降噪处理：利用第二代小波的提升算法将滑动数据窗中的数据进行指定层数的分解，得到各层的小波系数（细节信号）和最后一层的尺度系数（概貌信号）；对分解后的各层小波系数，建立相应的阈值规则，对其进行阈值量化处理；将处理后的小波系数，结合最后一层的尺度系数，逐级重构各层尺度系数，得到降噪处理后的信号；更新滑动数据窗最前端的陀螺采样点，将之作为姿态解算的输入陀螺信息，进行实时姿态解算输出；进入下一个采样周期，窗口平移一个采样点，继续进行上述降噪处理。本发明有效提高了光纤陀螺信号实时降噪的处理精度和反应速度，继而抑制了姿态信息误差；在应用过程中，滑动数据窗宽度、小波分解层数等参数的设置较为简便，可根据仿真实验得到的经验值进行设定。
附图说明
图1是本发明的第二代小波分解、重构示意图；
图2是本发明的第二代小波阈值法实时降噪流程图；
图3是本发明的正弦仿真信号降噪效果比较图；
图4是本发明的滑动数据窗宽度的影响图；
图5是本发明的陀螺信号降噪时域图；
图6是本发明的陀螺信号降噪频域图；
图7是本发明的姿态解算输出图。
具体实施方式
下面结合附图以及具体实施例对本发明作进一步描述，在此发明的示意性实施例以及说明用来解释本发明，但并不作为对本发明的限定。
如图1所示：
本发明的技术方案不限于上述具体实施例的限制，凡是根据本发明的技术方案做出的技术变形，均落入本发明的保护范围之内。
本发明包括以下步骤：
步骤1：确定滑动数据窗宽度及预测器系数P和更新器系数U的维数，并离线求得系数P和系数U；
P(M,1)＝[p1 p2 p3…pM]T
U(N,1)＝[u1 u2 u3…uN]T
其中P和U满足如下等式：
V(M,M)P(M,1)＝δ(M,1)
W(N,K)·Q(K,N)·U(N,1)＝δ(N,1)
其中，M和N为偶数，Md＝M/2-1,Nd＝N/2-1,其中：
δ(M,1)＝[100…0]T
M0＝2Md+1,N0＝2(Md+Nd+1)
m＝0,1,…,M-1,n＝-M0,-M0+2,…,-M0+2(M-1):V(M,M)(m+1,(n+M0)2+1)＝nm
i＝1,2,…,2M-1:q(i)＝0;q(2+2Md)＝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＝-N0,-N0+1,…,-N0+(K-1):W(N,K)(m+1,n+N0+1)＝nm
步骤2：捷联航姿系统预热结束后，采集光纤陀螺的输出数据，当采集数据量少于滑动数据窗宽度时，直接将之作为当前姿态解算的输入角速度信息ωib，解算后输出姿态信息（即开始采集数据后的小于第一个滑动数据窗宽度时间内的数据不作处理）；
步骤3：待采集陀螺数据量达到窗口宽度时，将第一段采集到的数据输入至去噪算法程序，进行降噪处理，在下一个采样点到来之前，输出降噪处理后的最后一个数据值，作为当前时刻陀螺信号的更新值；
1、采用提升方案的第二代小波变换对信号进行多尺度分解：具体步骤为（1）分裂、（2）预测、（3）更新。（1）中，分解层数越多，进行降噪处理的频带越宽，但同时可能会影响陀螺的动态输出特性，且增加了算法的复杂程度，故分解层数不宜超过5层。（2）中，随着信号的逐级分解，其对小波系数的阈值量化处理也随之进行，如图2所示。
（1）分裂：将原始的离散信号分解为两组相互关联的部分，通常分解为奇样本xo[n]和偶样本xe[n]：
xe[n]=x[2n]n=1,2···xo[n]=x[2n-1]n=1,2,···]]>
（2）预测：构造预测算子P与xe[n]作用，得到小波系数d[n]，作为由xe[n]预测xo[n]的误差：
d[n]＝xo[n]-P(xe[n])
或d(n)=xo(n)-Σk=1MpKxe(n-Md+(k-1))]]>
（3）更新：构造更新算子U与预测的小波系数d[n]作用，并叠加到xe[n]上，得到尺度系数c(n)，c(n)代表了原信号的粗略近似：
c[n]＝xe[n]+U(d[n])
或c(n)=xe(n)+Σk=1NuKd(n-Nd-1+(k-1))]]>
至此完成了第二代小波的第一层分解，重复步骤1）至3）可以继续进行分解，不同的是重复步骤1）时，原信号x[n]由上一层尺度系数cj(n)代替之。继续分裂、预测和更新L次后，有
x[n]＝{cl,dl,dl-1,…,d1}
其中，cl代表了信号的低频概貌部分（尺度系数），其余为信号的高频细节部分（小波系数）。
2、建立阈值降噪法的改进阈值处理策略，对分解后的各层小波系数进行阈值量化处理。具体步骤分为：（1）确定每层阈值、（2）确定阈值函数、（3）建立改进的阈值处理策略：
（1）确定每层阈值：采用通用阈值法，第j层阈值Tj由下式决定：
Tj=σj2ln(2j)]]>
其中，σj为第j层信号噪声标准差，该值由下式决定：
σj=median(dj)0.6745]]>
其中，median(dj)为第j层小波系数dj的中值。
（2）确定阈值函数：在选取基于二代小波变换的阈值函数时，应用硬阈值不仅与软阈值降噪效果相差不大，且其拥有更为简洁的处理策略，故选择硬阈值函数：
dj,k=|dj,k|,|dj,k|>Tj0,|dj,k|≤Tj]]>
（3）建立改进的阈值处理策略：改进的阈值处理策略结合了强制降噪对信号处理反应速度的提升，即将步骤3中的（1）处理分解后的第一层小波系数强制置零，而后分解的小波系数仍按照硬阈值函数的策略进行处理。
3、采用提升方案的第二代小波变换对信号的逐级重构，这种重构是由最后一层开始的。具体步骤分为：（1）恢复更新、（2）恢复预测、（3）合并。
（1）恢复更新：利用最后一层尺度系数cl[n]和经阈值量化处理后小波系数dl[n]，可恢复出l-1层尺度系数cl-1[n]的偶序列：
cl-1e[n]=cl[n]-U(dl[n])]]>
或cl-1e=cl(n)-Σk=1NuKdl(n-Nd-1+(k-1))]]>
（2）恢复预测：由1）得出的和进行阈值量化处理后的dl[n]预测出l-1层尺度系数cl-1[n]的奇序列：
cl-1o[n]=dl[n]+P(cl-1e[n])]]>
或cl-1o(n)=dl(n)+Σk=1MpKcl-1e(n-Md+(k-1))]]>
（3）合并：由1）和2）获得的cl-1[n]的奇序列和偶序列，合并恢复出j-1层尺度信号cl-1[n]：
cl-1[2n]=cl-1e[n]n=1,2···cl-1[2n-1]=cl-1o[n]n=1,2,···]]>
进行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为偶数，Md＝M/2-1,Nd＝N/2-1,且有：
δ(M,1)＝[100…0]T
M0＝2Md+1,N0＝2(Md+Nd+1)
m＝0,1,…,M-1,n＝-M0,-M0+2,…,-M0+2(M-1):V(M,M)(m+1,(n+M0)/2+1)＝nmi＝1,2,…,2M-1:q(i)＝0;q(2+2Md)＝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＝-N0,-N0+1,…,-N0+(K-1):W(N,K)(m+1,n+N0+1)＝nm
本发明所述步骤3中的去噪算法程序，其特征为：
建立以第二代小波变换的阈值去噪法为基础的降噪算法，利用提升算法对原始信号进行小波分解，分为分裂、预测和更新：
xe[n]=x[2n]n=1,2···xo[n]=x[2n-1]n=1,2,···]]>
d[n]＝xo[n]-P(xe[n])
c[n]＝xe[n]+U(d[n])
其中原始信号为x[n]，xo[n]为奇样本，xe[n]为偶样本，P和U分别为离线求得预测算子和更新算子，d[n]和c[n]分别为得到的小波和尺度系数。进行逐级分解至L层得到各层细节（小波系数dl,dl-1,…,d1）和最后一层概貌部分（尺度系数cl）。
x[n]＝{cl,dl,dl-1,…,d1}
建立阈值降噪法的改进阈值处理策略，对分解后的各层小波系数进行阈值量化处理，得到处理后的各层小波系数及未作处理的最后一层尺度系数。
采用提升方案的第二代小波变换对信号的逐级重构，由最后一层开始的，分为：恢复更新、恢复预测和合并：
cl-1e[n]=cl[n]-U(dl[n])]]>
cl-1o[n]=dl[n]+P(cl-1e[n])]]>
cl-1[2n]=cl-1e[n]n=1,2···cl-1[2n-1]=cl-1o[n]n=1,2,···]]>
其中cl[n]、dl[n]分别为最后一层尺度系数和经阈值量化处理后的小波系数，和分别为l-1层尺度系数cl-1[n]的偶序列和奇序列。逐级重构后得到复原的降噪处理后的信号x'[n]。
本发明所述建立阈值降噪法的改进阈值处理策略包括：
确定每层阈值：采用通用阈值法，第j层阈值Tj由下式决定：
Tj=σj2ln(2j)]]>
其中，σj为第j层信号噪声标准差，该值由下式决定：
σj=median(dj)0.6745]]>
其中，median(dj)为第j层小波系数dj的中值。
确定阈值函数：采用硬阈值函数：
dj,k=|dj,k|,|dj,k|>Tj0,|dj,k|≤Tj]]>
建立改进的阈值处理策略：将分解后的第一层小波系数强制置零，而后分解的小波系数仍按照硬阈值函数的策略进行处理。
dj,k=0,j=1|dj,k|,|dj,k|>Tj,j>10,|dj,k|≤Tj,j>1]]>
如图1所示：第二代小波变换对信号进行分解
首先确定维数为M的预测器P和维数为N更新器U：
P(M,1)＝[p1 p2 p3…pM]T
U(N,1)＝[u1 u2 u3…uN]T
其可由如下方程求得：
V(M,M)P(M,1)＝δ(M,1)
W(N,K)·Q(K,N)·U(N,1)＝δ(N,1)
对于偶数M和N，以及Md＝M/2-1,Nd＝N/2-1,利用上式即可求得P和U，其中
δ(M,1)＝[100…0]T
M0＝2Md+1,N0＝2(Md+Nd+1)
m＝0,1,…,M-1,n＝-M0,-M0+2,…,-M0+2(M-1):V(M,M)(m+1,(n+M0)2+1)＝nm
i＝1,2,…,2M-1:q(i)＝0;q(2+2Md)＝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＝-N0,-N0+1,…,-N0+(K-1):W(N,K)(m+1,n+N0+1)＝nm
预测器系数P和更新器系数U确定后，采用提升方案的第二代小波变换实施信号的分解，过程如下：
步骤（1）：分裂
将原始的离散信号分解为两组相互关联的部分，通常分解为奇样本xo[n]和偶样本xe[n]：
xe[n]=x[2n]n=1,2···xo[n]=x[2n-1]n=1,2,···]]>
步骤（2）：预测
构造预测算子P与xe[n]作用，得到小波系数d[n]，作为由xe[n]预测xo[n]的误差：
d[n]＝xo[n]-P(xe[n])
或d(n)=xo(n)-Σk=1MpKxe(n-Md+(k-1))]]>
步骤（3）：更新
构造更新算子U与预测的小波系数d[n]作用，并叠加到xe[n]上，得到尺度系数c(n)，c(n)代表了原信号的粗略近似：
c[n]＝xe[n]+U(d[n])
或c(n)=xe(n)+Σk=1NuKd(n-Nd-1+(k-1))]]>
至此完成了第二代小波的第一层分解，重复步骤（1）至（3）可以继续进行分解，不同的是重复步骤（1）时，原信号的位置由上一层尺度系数cj(n)代替之。继续分裂、预测和更新L次后，有
x[n]＝{cl,dl,dl-1,…,d1}
其中，cl代表了信号的低频概貌部分，其余为信号的高频细节部分。
确定每层阈值
采用通用阈值法确定每层阈值，通用阈值法是Donoho和Johnstone在高斯噪声模型下，应用多维独立正态变量决策理论得出的阈值。即当维数趋于无穷时,噪声系数的幅值大于阈值的概率趋于零，其中σj为第j层信号噪声标准差。
如图5和图6所示：对于陀螺信号来说,噪声的标准差是未知的,这就需要对噪声水平进行估计,噪声的标准差可以用各尺度细节小波系数绝对值的中值来进行估计：
σj=median(dj)0.6745]]>
式中，dj为第j层小波系数。
确定阈值函数
确定了各层的阈值之后，需选取一种合适的阈值量化方法，常用的有硬阈值法、软阈值法和兼顾这两种方法优点的半软阈值法。硬阈值法往往具有较大的方差，而软阈值法往往带来较大的偏差。虽然半软阈值法是以上两种方法的折中，不仅保留了较大的系数，而且具有连续性，但需要确定两个阈值，增加了运算复杂度。通过实验分析，传统的一代小波变换使用软阈值函数对光纤陀螺输出信号进行滤波，降噪效果要优于硬阈值函数。但在基于第二代小波变换的前提下，两者的滤波效果相差无几，且硬阈值函数具有较为简洁的处理手段，故采用硬阈值函数：
dj,k=|dj,k|,|dj,k|>Tj0,|dj,k|≤Tj]]>
建立强制降噪与硬阈值函数组合应用机制：
实际情况下，由舰载捷联惯导系统的实际应用环境，决定了光纤陀螺的动态特性频段要远低于装置在如飞机这种高速运动载体上的频段（特别是天向陀螺），根据这种有用信号多集中在中低频段的特征，可以将第二代小波变换后的第一层小波系数，也就是高频段信号，进行强制降噪处理，即将之置零。而继续分解的第二层、第三层等小波系数仍采用硬阈值函数进行降噪处理。这种小波系数处理策略继承了硬阈值函数与强制去噪的计算简便、快速的优点，且基本保留了原始信号中低频段的动态特性。
如图1所示：第二代小波变换重构信号，采用提升方案的第二代小波变换实施信号的重构：
步骤（1）：恢复更新
给定某层尺度系数c[n]和小波系数d[n]，可恢复出偶序列：
xe[n]＝c[n]-U(d[n])
或xe(n)=c(n)-Σk=1NuKd(n-Nd-1+(k-1))]]>
步骤（2）：恢复预测
由步骤（1）得出的xe(n)和给定的d[n]预测出奇序列：
xo[n]＝d[n]+P(xe[n])
或xo(n)=d-(n)+Σk=1MpKxe(n-Md+(k-1))]]>
步骤（3）：合并
由步骤（1）和步骤（2）获得奇序列和偶序列，将其合并恢复出原信号（或上一层尺度信号）。
x[2n]=xe[n]n=1,2···x[2n-1]=xo[n]n=1,2,···]]>
至此完成了第二代小波的一层信号的重构，重复步骤（1）至（3）可以继续进行重构，不同的是重复步骤（3）时，原信号的位置由上一层尺度系数cj(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、零偏稳定性/(0)·h-1为0.1147392、量化噪声/μrad为1.5247672、速率斜坡/((0)·h-2)为0.2481738；
传统方案的方差为2.2447、零偏稳定性/(0)·h-1为0.1028541、量化噪声/μrad为1.0404798、速率斜坡/((0)·h-2)为0.1753006；
改进方案的方差为0.8212、零偏稳定性/(0)·h-1为0.0848103、量化噪声/μrad为0.6125475、速率斜坡/((0)·h-2)为0.1069661。
在以下条件下，利用该方案对光线陀螺信号进行模拟实时降噪后，进行姿态解算实验：
光纤陀螺捷联航姿系统放置于单轴转台上，转台无基准；采集静态数据5小时，采样频率：98Hz；小波分解层数：3层滑动数据窗宽度：1024；降噪起始位置：3小时处；降噪数据长度：10000（约100s）；
系统初始参数如下：
载体初始位置：北纬45.7796°，东经126.6705°；地球自转角速度：7.2921158×10-5rads；重力加速度：9.78049ms2；地球半径：6378393.0m；
利用本发明方案从陀螺输出3小时处进行约100s的降噪处理，而后进行姿态解算，输出姿态如图7所示。虽然只进行了约100s的降噪处理，但该方案仍在降噪后减小纵横摇震荡幅值约7角秒。可以预见，实际应用时除开始阶段的窗口长度数据不做处理外，所有的陀螺原始数据均进行降噪处理，等效于增加降噪处理的时间，将进一步提高输出姿态的精度，抑制系统静态输出的姿态信息震荡。

资源描述

《一种光纤陀螺随机漂移实时滤波方法.pdf》由会员分享，可在线阅读，更多相关《一种光纤陀螺随机漂移实时滤波方法.pdf（20页珍藏版）》请在专利查询网上搜索。

1、(10)申请公布号 CN 103557856 A (43)申请公布日 2014.02.05 CN 103557856 A (21)申请号 201310508747.3 (22)申请日 2013.10.25 G01C 19/64(2006.01) (71)申请人哈尔滨工程大学地址 150001 黑龙江省哈尔滨市南岗区南通大街 145 号 (72)发明人孙枫夏健钟奔粤阳李敬春杨晓龙阮双双赵维珩 (74)专利代理机构北京科亿知识产权代理事务所 ( 普通合伙 ) 11350 代理人汤东凤 (54) 发明名称一种光纤陀螺随机漂移实时滤波方法 (57) 摘要本发明公开了一种光。

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

3、得到的经验值进行设定。 (51)Int.Cl. 权利要求书 2 页说明书 11 页附图 6 页 (19)中华人民共和国国家知识产权局 (12)发明专利申请权利要求书2页说明书11页附图6页 (10)申请公布号 CN 103557856 A CN 103557856 A 1/2 页 2 1. 一种光纤陀螺随机漂移实时滤波方法，其特征在于，包括以下步骤：步骤 1 ：确定滑动数据窗宽度及预测器系数 P 和更新器系数 U 的维数，并离线求得系数 P 和系数 U ；步骤 2 ：捷联航姿系统预热结束后，采集光纤陀螺的输出数据，当采集数据量少于滑动数据窗宽度时，直接将之作。

4、为当前姿态解算的输入角速度信息 ib，解算后输出姿态信息（即开始采集数据后的小于第一个滑动数据窗宽度时间内的数据不作处理）；步骤 3 ：待采集陀螺数据量达到窗口宽度时，将第一段采集到的数据输入至去噪算法程序，进行降噪处理，在下一个采样点到来之前，输出降噪处理后的最后一个数据值，作为当前时刻陀螺信号的更新值；步骤 4 ：将更新后的当前采样点作为姿态实时解算的输入角速度信息 ib，进行姿态的实时结算，输出当前姿态信息；步骤 5 ：进入下一个采样周期，滑动数据窗实时地平移一个点，总包含最新的采样点作为窗口的最后一个值，重复步骤3、步骤4和步骤5，。

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 为偶数， Md M/2-1,Nd N/2-1, 且有： (M,1) 1000T M0 2Md+1,N0 2(Md+Nd+1) m 0,1,M-1,n -M0,-M0+2,-M0+2(M-1):V(M,M)(m+1,(n+M0)/2+1) nm i 1,2,2M-1:q(i) 0;。

6、q(2+2Md) 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 -N0,-N0+1,-N0+(K-1):W(N,K)(m+1,n+N0+1) nm。 3. 根据权利权利要求 1 所述的光纤陀螺随机漂移实时滤波方法，其特征在于：所述步骤 3 中的去噪算法程序，其特征为：建立以第二代小波变换的阈值去噪法为基础的降噪算法，利用提升算法对原始信号进行小波分解，分为分裂、预测和更新： dn xon-P(xen) cn xen+U(。

7、dn) 其中原始信号为 xn， xon 为奇样本， xen 为偶样本， P 和 U 分别为离线求得预测算子和更新算子， dn 和 cn 分别为得到的小波和尺度系数。进行逐级分解至 L 层得到各权利要求书 CN 103557856 A 2 2/2 页 3 层细节（小波系数 dl,dl-1,d1）和最后一层概貌部分（尺度系数 cl）； xn cl,dl,dl-1,d1 建立阈值降噪法的改进阈值处理策略，对分解后的各层小波系数进行阈值量化处理，得到处理后的各层小波系数及未作处理的最后一层尺度系数。采用提升方案的第二代小波变换对信号的逐级重构，由最后一层开始的，分为：。

8、恢复更新、恢复预测和合并：其中 c ln、 dln 分别为最后一层尺度系数和经阈值量化处理后的小波系数，和分别为 l-1 层尺度系数 cl-1n 的偶序列和奇序列；逐级重构后得到复原的降噪处理后的信号 xn。 4. 根据权利权利要求 3 所述的光纤陀螺随机漂移实时滤波方法，其特征在于：所述建立阈值降噪法的改进阈值处理策略包括：确定每层阈值：采用通用阈值法，第 j 层阈值 Tj由下式决定：其中， j为第 j 层信号噪声标准差，该值由下式决定：其中， median(dj) 为第 j 层小波系数 dj的中值；确定阈值函数：采用硬阈值函数：建立。

9、改进的阈值处理策略：将分解后的第一层小波系数强制置零，而后分解的小波系数仍按照硬阈值函数的策略进行处理；权利要求书 CN 103557856 A 3 1/11 页 4 一种光纤陀螺随机漂移实时滤波方法技术领域 0001 本发明涉及一种信号的消噪方法，尤其涉及一种应用于捷联航姿系统中的光纤陀螺随机漂移实时滤波方法。背景技术 0002 光纤陀螺其是以 Sagnac 效应为基础而发展起来的新型全固态陀螺仪 , 作为主要惯性敏感元件之一，由于其具有潜在的最好的性能、成本比而倍受青睐 , 已成为新一代中低精度捷联惯性导航、制导系统中的理想惯性器件。但是，由于原理、。

10、加工与装配工艺不完善等造成陀螺漂移，成为系统误差的主要来源，严重影响着惯导系统的初始对准精度及导航精度。 0003 光纤陀螺的漂移包含两种成分：一种是具有确定性质的、有规律的、有系统性的漂移；另一种是由非确定性的随机性质干扰引起的随时间变化的随机漂移。对应前者，可以对这种确定性漂移进行建模 , 通过数学模型的软件补偿进行消除，这是比较容易做到的。而对于后者，目前采用较多的处理技术包括两种：一、对光纤陀螺的随机漂移进行时间序列建模，依据模型设计卡尔曼滤波器，通过卡尔曼滤波对随机漂移进行抑制；二、采用数字滤波器对光纤陀螺的漂移进行抑制，通过。

11、滤波器对光纤陀螺的漂移进行抑制。但实际情况中，由于光纤陀螺的输出受包括探测器噪声、光路噪声、电路噪声和环境变化引起的大量噪声影响 , 且这些噪声是时变的 , 故事先不能得到准确的误差统计特性，难以得到准确的漂移模型。因此，利用第一种方法无法得到理想的抑制效果。另外，由于随机漂移的频域特性不明确，增加了数字滤波器的设计难度，第二种方法也有其固有的缺陷。 0004 针对光纤陀螺信号中的随机漂移和噪声，从抑制噪声的角度出发，消除噪声 , 对光纤陀螺的随机漂移也有抑制和补偿作用，是比较可行的一种方法。在实际系统中，为满足降噪的实时性要求，通常需要降噪的方法具有。

12、在线处理能力。因此，如何设计一种既优化了降噪效果，又提高了算法运算速度、简化了降噪处理过程的实时消噪方法，对光纤陀螺的不确定性随机漂移和噪声进行实时直接滤除，对捷联航姿系统精度的提高起着至关重要的作用。发明内容 0005 本发明的目的是要提供一种应用于捷联航姿系统中的光纤陀螺随机漂移实时滤波方法。 0006 为达到上述目的，本发明是按照以下技术方案实施的： 0007 本发明包括以下步骤： 0008 步骤 1 ：确定滑动数据窗宽度及预测器系数 P 和更新器系数 U 的维数，并离线求得系数 P 和系数 U ； 0009 步骤 2 ：捷联航姿系统预热结束后，采集。

13、光纤陀螺的输出数据，当采集数据量少于滑动数据窗宽度时，直接将之作为当前姿态解算的输入角速度信息ib，解算后输出姿态信说明书 CN 103557856 A 4 2/11 页 5 息（即开始采集数据后的小于第一个滑动数据窗宽度时间内的数据不作处理）； 0010 步骤 3 ：待采集陀螺数据量达到窗口宽度时，将第一段采集到的数据输入至去噪算法程序，进行降噪处理，在下一个采样点到来之前，输出降噪处理后的最后一个数据值，作为当前时刻陀螺信号的更新值； 0011 步骤 4 ：将更新后的当前采样点作为姿态实时解算的输入角速度信息 ib，进行姿态的实时结算，输出当前姿态。

14、信息； 0012 步骤 5 ：进入下一个采样周期，滑动数据窗实时地平移一个点，总包含最新的采样点作为窗口的最后一个值，重复步骤3、步骤4和步骤5，完成捷联航姿系统中光纤陀螺输出信号的实时降噪； 0013 本发明所述步骤 1 中的离线求得预测器系数 P、更新器系数 U 满足如下等式： 0014 V(M,M)P(M,1) (M,1) 0015 W(N,K)Q(K,N)U(N,1) (N,1) 0016 其中 M 和 N 为偶数， Md M/2-1,Nd N/2-1, 且有： 0017 (M,1) 1000T 0018 M0 2Md+1,N0 2(Md+Nd+1) 0019。

15、 m 0,1,M-1,n -M0,-M0+2,-M0+2(M-1):V(M,M)(m+1,(n+M0)/2+1) nmi 1,2,2M-1:q(i) 0;q(2+2Md) 1 0020 i 1,2,2M:q(2i-1) p(i) 0021 L length(q) 0022 K 2(M-1)+L 0023 Q(K,N) 0;(n 1,2,N:Q(2n-1:2n-1+L-1,n) q) 0024 m 0,1,N-1,n -N0,-N0+1,-N0+(K-1):W(N,K)(m+1,n+N0+1) nm 0025 本发明所述步骤 3 中的去噪算法程序，其特征为： 0026 建立以第二代小波变换的。

16、阈值去噪法为基础的降噪算法，利用提升算法对原始信号进行小波分解，分为分裂、预测和更新： 0027 0028 dn xon-P(xen) 0029 cn xen+U(dn) 0030 其中原始信号为 xn， xon 为奇样本， xen 为偶样本， P 和 U 分别为离线求得预测算子和更新算子， dn 和 cn 分别为得到的小波和尺度系数。进行逐级分解至 L 层得到各层细节（小波系数 dl,dl-1,d1）和最后一层概貌部分（尺度系数 c l）。 0031 xn cl,dl,dl-1,d1 0032 建立阈值降噪法的改进阈值处理策略，对分解后的各层小波系数进行阈值量化处理。

17、，得到处理后的各层小波系数及未作处理的最后一层尺度系数。 0033 采用提升方案的第二代小波变换对信号的逐级重构，由最后一层开始的，分为：恢复更新、恢复预测和合并： 0034 说明书 CN 103557856 A 5 3/11 页 6 0035 0036 0037 其中 cln、 dln 分别为最后一层尺度系数和经阈值量化处理后的小波系数，和分别为 l-1 层尺度系数 cl-1n 的偶序列和奇序列。逐级重构后得到复原的降噪处理后的信号 xn。 0038 本发明所述建立阈值降噪法的改进阈值处理策略包括： 0039 确定每层阈值：采用通用阈值法，第 j 层阈值 Tj。

18、由下式决定： 0040 0041 其中， j为第 j 层信号噪声标准差，该值由下式决定： 0042 0043 其中， median(dj) 为第 j 层小波系数 dj的中值。 0044 确定阈值函数：采用硬阈值函数： 0045 0046 建立改进的阈值处理策略：将分解后的第一层小波系数强制置零，而后分解的小波系数仍按照硬阈值函数的策略进行处理。 0047 0048 本发明的有益效果是： 0049 本发明是一种光纤陀螺随机漂移实时滤波方法，与现有技术相比，本发明在船用捷联航姿系统中，对采集到的光纤陀螺信号进行实时降噪处理：利用第二代小波的提升算法将滑动数据窗中。

19、的数据进行指定层数的分解，得到各层的小波系数（细节信号）和最后一层的尺度系数（概貌信号）；对分解后的各层小波系数，建立相应的阈值规则，对其进行阈值量化处理；将处理后的小波系数，结合最后一层的尺度系数，逐级重构各层尺度系数，得到降噪处理后的信号；更新滑动数据窗最前端的陀螺采样点，将之作为姿态解算的输入陀螺信息，进行实时姿态解算输出；进入下一个采样周期，窗口平移一个采样点，继续进行上述降噪处理。本发明有效提高了光纤陀螺信号实时降噪的处理精度和反应速度，继而抑制了姿态信息误差；在应用过程中，滑动数据窗宽度、小波分解层数等参数的设置较为简。

20、便，可根据仿真实验得到的经验值进行设定。附图说明 0050 图 1 是本发明的第二代小波分解、重构示意图；说明书 CN 103557856 A 6 4/11 页 7 0051 图 2 是本发明的第二代小波阈值法实时降噪流程图； 0052 图 3 是本发明的正弦仿真信号降噪效果比较图； 0053 图 4 是本发明的滑动数据窗宽度的影响图； 0054 图 5 是本发明的陀螺信号降噪时域图； 0055 图 6 是本发明的陀螺信号降噪频域图； 0056 图 7 是本发明的姿态解算输出图。具体实施方式 0057 下面结合附图以及具体实施例对本发明作进一步描述，在此发明的示意。

21、性实施例以及说明用来解释本发明，但并不作为对本发明的限定。 0058 如图 1 所示： 0059 本发明的技术方案不限于上述具体实施例的限制，凡是根据本发明的技术方案做出的技术变形，均落入本发明的保护范围之内。 0060 本发明包括以下步骤： 0061 步骤 1 ：确定滑动数据窗宽度及预测器系数 P 和更新器系数 U 的维数，并离线求得系数 P 和系数 U ； 0062 P(M,1) p1 p2 p3pMT 0063 U(N,1) u1 u2 u3uNT 0064 其中 P 和 U 满足如下等式： 0065 V(M,M)P(M,1) (M,1) 0066 W(N,K)Q(。

22、K,N)U(N,1) (N,1) 0067 其中， M 和 N 为偶数， Md M/2-1,Nd N/2-1, 其中： 0068 (M,1) 1000T 0069 M0 2Md+1,N0 2(Md+Nd+1) 0070 m 0,1,M-1,n -M0,-M0+2,-M0+2(M-1):V(M,M)(m+1,(n+M0)2+1) nm 0071 i 1,2,2M-1:q(i) 0;q(2+2Md) 1 0072 i 1,2,2M:q(2i-1) p(i) 0073 L length(q) 0074 K 2(M-1)+L 0075 Q(K,N) 0;(n 1,2,N:Q(2n-1:2n-1+L-。

23、1,n) q) 0076 m 0,1,N-1,n -N0,-N0+1,-N0+(K-1):W(N,K)(m+1,n+N0+1) nm 0077 步骤 2 ：捷联航姿系统预热结束后，采集光纤陀螺的输出数据，当采集数据量少于滑动数据窗宽度时，直接将之作为当前姿态解算的输入角速度信息ib，解算后输出姿态信息（即开始采集数据后的小于第一个滑动数据窗宽度时间内的数据不作处理）； 0078 步骤 3 ：待采集陀螺数据量达到窗口宽度时，将第一段采集到的数据输入至去噪算法程序，进行降噪处理，在下一个采样点到来之前，输出降噪处理后的最后一个数据值，作为当前时刻陀螺信号的更新值；。

24、 0079 1、采用提升方案的第二代小波变换对信号进行多尺度分解：具体步骤为（1）分裂、（2）预测、（3）更新。（1）中，分解层数越多，进行降噪处理的频带越宽，但同时可能会影说明书 CN 103557856 A 7 5/11 页 8 响陀螺的动态输出特性，且增加了算法的复杂程度，故分解层数不宜超过 5 层。（2）中，随着信号的逐级分解，其对小波系数的阈值量化处理也随之进行，如图 2 所示。 0080 （1）分裂：将原始的离散信号分解为两组相互关联的部分，通常分解为奇样本 xon 和偶样本 xen ： 0081 0082 （2）预测：构。

25、造预测算子 P 与 xen 作用，得到小波系数 dn，作为由 xen 预测 xon 的误差： 0083 dn xon-P(xen) 0084 或 0085 （3）更新：构造更新算子 U 与预测的小波系数 dn 作用，并叠加到 xen 上，得到尺度系数 c(n)， c(n) 代表了原信号的粗略近似： 0086 cn xen+U(dn) 0087 或 0088 至此完成了第二代小波的第一层分解，重复步骤 1）至 3）可以继续进行分解，不同的是重复步骤 1）时，原信号 xn 由上一层尺度系数 cj(n) 代替之。继续分裂、预测和更新 L 次后，有 0089 xn。

26、 cl,dl,dl-1,d1 0090 其中， cl代表了信号的低频概貌部分（尺度系数），其余为信号的高频细节部分（小波系数）。 0091 2、建立阈值降噪法的改进阈值处理策略，对分解后的各层小波系数进行阈值量化处理。具体步骤分为：（1）确定每层阈值、（2）确定阈值函数、（3）建立改进的阈值处理策略： 0092 （1）确定每层阈值：采用通用阈值法，第 j 层阈值 Tj由下式决定： 0093 0094 其中， j为第 j 层信号噪声标准差，该值由下式决定： 0095 0096 其中， median(dj) 为第 j 层小波系数 dj的中值。 0097 。

27、（2）确定阈值函数：在选取基于二代小波变换的阈值函数时，应用硬阈值不仅与软阈值降噪效果相差不大，且其拥有更为简洁的处理策略，故选择硬阈值函数： 0098 0099 （3）建立改进的阈值处理策略：改进的阈值处理策略结合了强制降噪对信号处理说明书 CN 103557856 A 8 6/11 页 9 反应速度的提升，即将步骤 3 中的（1）处理分解后的第一层小波系数强制置零，而后分解的小波系数仍按照硬阈值函数的策略进行处理。 0100 3、采用提升方案的第二代小波变换对信号的逐级重构，这种重构是由最后一层开始的。具体步骤分为：（1）恢复更新、（2）恢。

28、复预测、（3）合并。 0101 （1）恢复更新：利用最后一层尺度系数cln和经阈值量化处理后小波系数dln，可恢复出 l-1 层尺度系数 cl-1n 的偶序列： 0102 0103 或 0104 （2）恢复预测：由 1）得出的和进行阈值量化处理后的 dln 预测出 l-1 层尺度系数 cl-1n 的奇序列： 0105 0106 或 0107 （3）合并：由 1）和 2）获得的 cl-1n 的奇序列和偶序列，合并恢复出 j-1 层尺度信号 cl-1n ： 0108 0109 进行 L 次的恢复更新、恢复预测、合并之后，可以复原出原始信号 xn。 011。

29、0 4、更新滑动数据窗内的末端数据，即将降噪处理后返回的此段信号的最后一个采样点的值更新当前采样点，完成对当前陀螺数据的实时更新。 0111 步骤 4 ：将更新后的当前采样点作为姿态实时解算的输入角速度信息 ib，进行姿态的实时结算，输出当前姿态信息； 0112 步骤 5 ：进入下一个采样周期，滑动数据窗实时地平移一个点，总包含最新的采样点作为窗口的最后一个值，重复步骤3、步骤4和步骤5，完成捷联航姿系统中光纤陀螺输出信号的实时降噪； 0113 本发明所述步骤 1 中的离线求得预测器系数 P、更新器系数 U 满足如下等式： 0114 V(M,M)P(M,。

30、1) (M,1) 0115 W(N,K)Q(K,N)U(N,1) (N,1) 0116 其中 M 和 N 为偶数， Md M/2-1,Nd N/2-1, 且有： 0117 (M,1) 1000T 0118 M0 2Md+1,N0 2(Md+Nd+1) 0119 m 0,1,M-1,n -M0,-M0+2,-M0+2(M-1):V(M,M)(m+1,(n+M0)/2+1) nmi 1,2,2M-1:q(i) 0;q(2+2Md) 1 0120 i 1,2,2M:q(2i-1) p(i) 0121 L length(q) 0122 K 2(M-1)+L 说明书 CN 103557856 A 。

31、9 7/11 页 10 0123 Q(K,N) 0;(n 1,2,N:Q(2n-1:2n-1+L-1,n) q) 0124 m 0,1,N-1,n -N0,-N0+1,-N0+(K-1):W(N,K)(m+1,n+N0+1) nm 0125 本发明所述步骤 3 中的去噪算法程序，其特征为： 0126 建立以第二代小波变换的阈值去噪法为基础的降噪算法，利用提升算法对原始信号进行小波分解，分为分裂、预测和更新： 0127 0128 dn xon-P(xen) 0129 cn xen+U(dn) 0130 其中原始信号为 xn， xon 为奇样本， xen 为偶样本， P 和 U 分别。

32、为离线求得预测算子和更新算子， dn 和 cn 分别为得到的小波和尺度系数。进行逐级分解至 L 层得到各层细节（小波系数 dl,dl-1,d1）和最后一层概貌部分（尺度系数 cl）。 0131 xn cl,dl,dl-1,d1 0132 建立阈值降噪法的改进阈值处理策略，对分解后的各层小波系数进行阈值量化处理，得到处理后的各层小波系数及未作处理的最后一层尺度系数。 0133 采用提升方案的第二代小波变换对信号的逐级重构，由最后一层开始的，分为：恢复更新、恢复预测和合并： 0134 0135 0136 0137 其中 cln、 dln 分别为最后一层尺度系数和经阈。

33、值量化处理后的小波系数，和分别为 l-1 层尺度系数 cl-1n 的偶序列和奇序列。逐级重构后得到复原的降噪处理后的信号 xn。 0138 本发明所述建立阈值降噪法的改进阈值处理策略包括： 0139 确定每层阈值：采用通用阈值法，第 j 层阈值 Tj由下式决定： 0140 0141 其中， j为第 j 层信号噪声标准差，该值由下式决定： 0142 0143 其中， median(dj) 为第 j 层小波系数 dj的中值。 0144 确定阈值函数：采用硬阈值函数： 0145 0146 建立改进的阈值处理策略：将分解后的第一层小波系数强制置零，而后分解的小说明书。

34、 CN 103557856 A 10 8/11 页 11 波系数仍按照硬阈值函数的策略进行处理。 0147 0148 如图 1 所示：第二代小波变换对信号进行分解 0149 首先确定维数为 M 的预测器 P 和维数为 N 更新器 U ： 0150 P(M,1) p1 p2 p3pMT 0151 U(N,1) u1 u2 u3uNT 0152 其可由如下方程求得： 0153 V(M,M)P(M,1) (M,1) 0154 W(N,K)Q(K,N)U(N,1) (N,1) 0155 对于偶数 M 和 N，以及 Md M/2-1,Nd N/2-1, 利用上式即可求得 P 和 U，其中 01。

35、56 (M,1) 1000T 0157 M0 2Md+1,N0 2(Md+Nd+1) 0158 m 0,1,M-1,n -M0,-M0+2,-M0+2(M-1):V(M,M)(m+1,(n+M0)2+1) nm 0159 i 1,2,2M-1:q(i) 0;q(2+2Md) 1 0160 i 1,2,2M:q(2i-1) p(i) 0161 L length(q) 0162 K 2(M-1)+L 0163 Q(K,N) 0;(n 1,2,N:Q(2n-1:2n-1+L-1,n) q) 0164 m 0,1,N-1,n -N0,-N0+1,-N0+(K-1):W(N,K)(m+1,n+N0+1)。

36、 nm 0165 预测器系数 P 和更新器系数 U 确定后，采用提升方案的第二代小波变换实施信号的分解，过程如下： 0166 步骤（1）：分裂 0167 将原始的离散信号分解为两组相互关联的部分，通常分解为奇样本 xon 和偶样本 xen ： 0168 0169 步骤（2）：预测 0170 构造预测算子 P 与 xen 作用，得到小波系数 dn，作为由 xen 预测 xon 的误差： 0171 dn xon-P(xen) 0172 或 0173 步骤（3）：更新 0174 构造更新算子 U 与预测的小波系数 dn 作用，并叠加到 xen 上，得到尺度系。

37、数 c(n)， c(n) 代表了原信号的粗略近似：说明书 CN 103557856 A 11 9/11 页 12 0175 cn xen+U(dn) 0176 或 0177 至此完成了第二代小波的第一层分解，重复步骤（1）至（3）可以继续进行分解，不同的是重复步骤（1）时，原信号的位置由上一层尺度系数 cj(n) 代替之。继续分裂、预测和更新 L 次后，有 0178 xn cl,dl,dl-1,d1 0179 其中， cl代表了信号的低频概貌部分，其余为信号的高频细节部分。 0180 确定每层阈值 0181 采用通用阈值法确定每层阈值，通用阈值法是Dono。

38、ho和Johnstone在高斯噪声模型下，应用多维独立正态变量决策理论得出的阈值。即当维数趋于无穷时 , 噪声系数的幅值大于阈值的概率趋于零，其中 j为第 j 层信号噪声标准差。 0182 如图 5 和图 6 所示：对于陀螺信号来说 , 噪声的标准差是未知的 , 这就需要对噪声水平进行估计 , 噪声的标准差可以用各尺度细节小波系数绝对值的中值来进行估计： 0183 0184 式中， dj为第 j 层小波系数。 0185 确定阈值函数 0186 确定了各层的阈值之后，需选取一种合适的阈值量化方法，常用的有硬阈值法、软阈值法和兼顾这两种方法优点的半软阈值法。硬阈值法往往具有。

39、较大的方差，而软阈值法往往带来较大的偏差。虽然半软阈值法是以上两种方法的折中，不仅保留了较大的系数，而且具有连续性，但需要确定两个阈值，增加了运算复杂度。通过实验分析，传统的一代小波变换使用软阈值函数对光纤陀螺输出信号进行滤波，降噪效果要优于硬阈值函数。但在基于第二代小波变换的前提下，两者的滤波效果相差无几，且硬阈值函数具有较为简洁的处理手段，故采用硬阈值函数： 0187 0188 建立强制降噪与硬阈值函数组合应用机制： 0189 实际情况下，由舰载捷联惯导系统的实际应用环境，决定了光纤陀螺的动态特性频段要远低于装置在如飞机这种高速运动载体上的频段（。

40、特别是天向陀螺），根据这种有用信号多集中在中低频段的特征，可以将第二代小波变换后的第一层小波系数，也就是高频段信号，进行强制降噪处理，即将之置零。而继续分解的第二层、第三层等小波系数仍采用硬阈值函数进行降噪处理。这种小波系数处理策略继承了硬阈值函数与强制去噪的计算简便、快速的优点，且基本保留了原始信号中低频段的动态特性。 0190 如图 1 所示：第二代小波变换重构信号，采用提升方案的第二代小波变换实施信号的重构： 0191 步骤（1）：恢复更新说明书 CN 103557856 A 12 10/11 页 13 0192 给定某层尺度系数 cn 。

41、和小波系数 dn，可恢复出偶序列： 0193 xen cn-U(dn) 0194 或 0195 步骤（2）：恢复预测 0196 由步骤（1）得出的 xe(n) 和给定的 dn 预测出奇序列： 0197 xon dn+P(xen) 0198 或 0199 步骤（3）：合并 0200 由步骤（1）和步骤（2）获得奇序列和偶序列，将其合并恢复出原信号（或上一层尺度信号）。 0201 0202 至此完成了第二代小波的一层信号的重构，重复步骤（1）至（3）可以继续进行重构，不同的是重复步骤（3）时，原信号的位置由上一层尺度系数 cj(n) 代替之。

42、。继续进行分解层数次的恢复更新、恢复预测、合并之后，可以复原出原始信号 xn。作为实时降噪的应用，此时只输出重构信号的最后一个值，即只更新最新的采样点。 0203 至此，实时降噪处理完毕，可将实时输出的最新采样点作为姿态结算的输入陀螺信号，进行实时姿态结算。 0204 如图 3 所示：本发明较传统基于一代小波的实时降噪方案，去噪后的模拟信号偏离期望的程度，即方差小于传统方案去噪后的结果，且去噪后的信噪比较传统方案提高约 3db，较含噪信号提高约 6.7db，采用文中提到的改进方案的实时降噪效果优于传统一代小波的实时降噪方案。 0205 方案类型：染噪。

43、信号的方差为 1.4992、信噪比为 -3.0198 ；传统方案的方差为 0.9260、信噪比为 0.6642 ；改进方案的方差为 0.7117、信噪比为 37120. 0206 如图 4 所示：在以下条件下，为比较滑动数据窗宽度的影响，利用该方案对一组染噪正弦信号进行模拟在线实时降噪实验： 0207 正弦信号频率： 0.1Hz ；幅值： 1 ；采样点数： 10000 ；采样频率： 100Hz ；噪声特征：每个采样点叠加一个标准正态分布的随机数；小波分解层数： 3 层；滑动数据窗宽度分别为： 128、 512、 1024 ； 0208 结。

44、果如图 4 所示，可以看出，滑动数据窗对本发明的实时降噪效果影响很小。 0209 在以下条件下，利用该方案对光纤陀螺输出信号进行模拟在线实时降噪实验： 0210 光纤陀螺捷联航姿系统放置于单轴转台上，转台无基准；采集静态数据 5 小时，采样频率： 98Hz ；小波分解层数： 3 层；滑动数据窗宽度： 1024 ；数据处理长度： 10000 ； 0211 利用本发明所述方法得到结果如图 5、图 6 所示，利用 Allan 方差分析法对本发明处理后的信号误差系数进行分析，比较了基于一代小波的实时降噪方案（传统方案）与本发明采用的基于二代小波变换的实时。

45、降噪方案（改进方案）的去噪性能：说明书 CN 103557856 A 13 11/11 页 14 0212 降噪方案：原始信号的方差为 4.7625、零偏稳定性 /(0) h-1为 0.1147392、量化噪声 /rad 为 1.5247672、速率斜坡 /(0)h-2) 为 0.2481738 ； 0213 传统方案的方差为 2.2447、零偏稳定性 /(0)h-1为 0.1028541、量化噪声 /rad 为 1.0404798、速率斜坡 /(0)h-2) 为 0.1753006 ； 0214 改进方案的方差为 0.8212、零偏稳定性 /(0)h-1为 0。

46、.0848103、量化噪声 /rad 为 0.6125475、速率斜坡 /(0)h-2) 为 0.1069661。 0215 在以下条件下，利用该方案对光线陀螺信号进行模拟实时降噪后，进行姿态解算实验： 0216 光纤陀螺捷联航姿系统放置于单轴转台上，转台无基准；采集静态数据 5 小时，采样频率： 98Hz ；小波分解层数： 3层滑动数据窗宽度： 1024 ；降噪起始位置： 3小时处；降噪数据长度： 10000（约 100s）； 0217 系统初始参数如下： 0218 载体初始位置：北纬 45.7796 ，东经 126.6705 ；地球自转角速度： 7.292115810-5rads ；重力加速度： 9.78049ms2；地球半径： 6378393.0m ； 0219 利用本发明方案从陀螺输出 3 小时处进行约 100s 的降噪处理，而后进行姿态解算，输出姿态如图 7 所示。虽然只进行了约 100s 的降噪处理，但。

展开阅读全文