一种机抖激光陀螺惯性测量单元数字滤波器设计方法

摘要

一种机抖激光陀螺惯性测量单元数字滤波器设计方法，其特征在于滤波器包括带阻部分和低通部分；首先根据机抖激光陀螺机抖噪声特点设计线性相位IIR数字滤波器带阻部分；然后在时域内采用带遗忘因子的最小二乘参数估计使线性相位IIR数字滤波器低通部分逼近线性相位FIR低通滤波器，获得近似线性相位；最后建立频域内线性相位IIR数字滤波器低通部分目标函数，利用DFP法优化目标函数，修正时域设计结果，使线性相位IIR数字滤波器低通部分具有所期望的频率特性。其优点是相位近似线性，运算量低，群延时小；降低了滤波器对IMU硬件的要求，提高了IMU测量实时性，保证了时延补偿精度，为IMU高精度实时测量奠定了基础。

说明书

技术领域

本发明属于惯性技术领域，涉及一种机抖激光陀螺惯性测量单元(IMU)数字滤波器设计方法，可以应用于位置姿态测量系统(Position and OrientationSystem，POS)和捷联惯性导航系统(Strap down Navigation System，SINS)，以及SINS/GPS(Global Position System，全球定位系统)组合导航系统的数据预处理。

背景技术

惯性测量单元(Inertial Measurement Unit，IMU)是POS、SINS、SINS/GPS等惯性测量系统的核心部件，惯性测量系统是以牛顿力学定律为基础，利用IMU测量载体的线运动和角运动参数，在给定初始条件下，由计算机推算得到载体的位置速度姿态参数。它是一种自主式的位置速度姿态测量单元，完全依靠机载设备自主地完成测量任务，和外界不发生任何光、电联系，具有隐蔽性好，工作不受气象条件限制的优点，因此在航空、航天和航海领域得到广泛的使用。IMU主要包含陀螺仪模块、加速度计模块、数据采集模块和数据处理模块，陀螺仪的精度很大程度上决定了IMU的测量精度。由于激光陀螺具有高动态、低成本、高可靠、性能价格比高的特点，在军用、民用方面被广泛应用。目前全世界应用的激光陀螺大部分是美国霍尼尔公司和利顿公司的机抖激光陀螺，而在国内激光陀螺主要是国防科技大学和航空618所的机抖激光陀螺，由于机抖激光陀螺的突出优点与广泛应用，机抖激光陀螺IMU已经成主要发展方向之一。

激光陀螺由于“闭锁效应”的存在，通常采用机械抖动偏频，这种激光陀螺即为机抖激光陀螺，机抖抖动偏频是在激光陀螺敏感轴方向输入偏置角速率，使其工作点全部或大部分从锁区偏置出来，因此机抖激光陀螺输出的信号包含抖动角速率分量，需要进行抖动解调，剥离由机械抖动引起的抖动角速率分量，保留IMU载体的运动角速率分量。抖动解调的方法有光学补偿法，整周期计数法，误差信号修正法和高速采样低通滤波方法。其中高速采样低通滤波方法解调精度高、设计简单，是目前主要的解调方式之一，该方法的关键是低通数字滤波器的设计。从数字滤波器的冲激响应来看，滤波器可分有限长单位冲激响应(FIR)滤波器和无限长单位冲激响应(IIR)滤波器。机抖激光陀螺IMU信号滤波采用FIR滤波器能够获得严格的线性相位，但由于机抖激光陀螺机械抖动噪声强度远远超出有用信号强度，要想达到理想的幅频特性，FIR滤波器阶数将会很高，运算量大，对IMU硬件系统提出了苛刻的要求，而且IMU输出的群延时也与FIR滤波器阶数成比例，因此FIR滤波器使得IMU输出数据延时大。机抖激光陀螺IMU采用IIR滤波器可以获得良好的噪声抑制效果，而且IIR滤波器阶数低，IIR滤波器包括巴特沃斯滤波器、切比雪夫I型滤波器、切比雪夫II型滤波器和椭圆滤波器，按照这些准则设计的IIR滤波器，不仅相位非线性严重，使得不同频率的IMU测量信号输出延时差别大，导致IMU测量精度下降，而且其群延时也很大，降低了IMU测量的实时性。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提供一种机抖激光陀螺IMU数字滤波器设计方法，该方法设计的数字滤波器具有相位近似线性，运算量小，群延时小，群延时波动小的优点，降低了滤波器对IMU硬件系统的要求，提高了IMU测量实时性和时延补偿精度，为IMU高精度实时测量打下基础。

本发明技术解决方案为：一种机抖激光陀螺IMU数字滤波器设计方法，其特点在于包括下列步骤：

(1)根据机抖激光陀螺IMU性能要求，确定线性相位IIR数字滤波器带阻部分和低通部分的指标参数：带阻部分阻带衰减参数A_stop1，中心频率f₀；低通部分阻带衰减参数A_stop2，低通部分通带纹波A_pass2，低通部分群延时波动G_dw；线性相位IIR数字滤波器整体阻带衰减参数A_stop。并确定机抖激光陀螺抖动噪声峰值频率f_n和强度A_n以及机抖激光陀螺抖频波动范围f_range；

(2)采用通带最平坦准则，以中心频率f₀＝f_n，设计线性相位IIR数字滤波器带阻部分，获得其传递函数分子系数a_s1，a_s2，…，a_sn和分母系数b_s0，b_s1，b_s2，…，b_sn，其传递函数$>H_s\left(z\right)=\frac{b_{s0}+b_{s1}{z}^{-1}+b_{s2}{z}^{-2}+...+b_{\mathrm{sn}}{z}^{-n}}{1+a_{s1}{z}^{-1}+a_{s2}{z}^{-2}+...+a_{\mathrm{sn}}{z}^{-n}};$>

(3)在机抖激光陀螺的抖频调节范围为3f_range～5f_range的情况下测试H_s(z)，采集经H_s(z)滤波后的信号进行功率谱分析，如果抖动噪声峰值低于0.7A_n，调节H_s(z)的中心频率f₀重新进行步骤(2)、步骤(3)；

(4)在时域内，以线性相位FIR低通滤波器与线性相位IIR数字滤波器低通部分H_l(z)的误差为量测误差，利用最小二乘估计使H_l(z)相位线性化，获得H_l(z)初始分子系数a₁(0)，a₂(0)，…，a_n(0)和分母系数b₀(0)，b₁(0)…b_n(0)，即初始$>H_l\left(z\right)=\frac{b_0\left(0\right)+b_1\left(0\right)z^{-1}+b_2\left(0\right)z^{-2}+...+b_n\left(0\right)z^{-n}}{1+a_1\left(0\right)z^{-1}+a_2\left(0\right)z^{-2}+...+a_n\left(0\right)z^{-n}};$>

(5)在频域内，建立优化H_l(z)的目标函数J，J的相频误差评价函数权系数为K_p，幅频误差评价函数过渡带和阻带权系数为K_M2，K_M3；以a₁(0)，a₂(0)，…，a_n(0)，b₀(0)，b₁(0)，b₂(0)…b_n(0)为初值，采用多变量无约束非线性极小值问题求解方法DFP(Daridon-Fletcher-Powell)法对J进行优化，获得频域优化后的H_l(z)的分子和分母系数a_l1(n)，a_l2(n)，…，a_ln(n)，b_l0(n)，b_l1(n)…b_ln(n)；

(6)通过单位冲激响应测试法，测试频域优化后的H_l(z)，并采集经H_l(z)滤波后机抖激光陀螺信号，进行功率谱分析，如果H_l(z)的通带纹波大于A_pass2，同时减小K_p，K_M2，K_M3，如果H_l(z)的群延时波动大于G_dw，增大权系数K_p，如果机抖激光陀螺抖动噪声衰减低于A_stop1，增大权系数K_M3，重新进行步骤(6)；

(7)利用公式H(z)＝H_s(z)H_l(z)求线性相位IIR数字滤波器传递函数H(z)，并采集滤波后信号进行功率谱分析，如果H(z)阻带内噪声衰减低于A_stop，增大A_stop1、A_stop2，重新进行，求得H(z)的分子系数b₀(n)，b₁(n)，b₂(n)…b_n(n)和分母系数a₁(n)，a₂(n)，…，a_n(n)，最终获得线性相位IIR数字滤波器传递函数$>H\left(z\right)=\frac{b_0\left(n\right)+b_1\left(n\right)z^{-1}+b_2\left(n\right)z^{-2}+...+b_n\left(n\right)z^{-n}}{1+a_1\left(n\right)z^{-1}+a_2\left(n\right)z^{-2}+...+a_n\left(n\right)z^{-n}},$>完成线性相位IIR数字滤波器设计。

本发明的原理是：机抖激光陀螺信号包含被测角速率信息，机械抖动噪声和其它噪声。其中机械抖动噪声具有频率集中，噪声强度大。一般单独采用FIR低通滤波器阶数会很高，其主要用于抑制抖动噪声，但抖动噪声频带很窄，高阶的FIR滤波器却将整个阻带信号全部衰减，FIR滤波器实际利用率低。因此本发明采用二阶或三阶线性相位IIR数字滤波器带阻部分，抑制机抖激光陀螺抖动噪声，该线性相位IIR数字滤波器带阻部分在低频段幅值波动和相位非线性比线性相位IIR数字滤波器低通部分小一个数量级以上，可以忽略，对整个数字滤波器通带内的幅相特性没有影响。为进一步保证整个数字滤波器的滤波器效果，设计一个线性相位FIR低通滤波器，由于IMU不要求数字滤波器具有严格的线性相位，滤波器非线性只要对IMU测量精度影响可以忽略即可，因此考虑到如果能够设计一个线性相位的IIR滤波器低通部分，整个滤波器运算量会进一步降低。但是求解IIR滤波器最优问题是一个多变量无约束优化问题，这种问题在参数较多时很难求得优化解。本发明采用时域设计与频域设计的优化方法，首先在时域内，利用最小二乘法IIR滤波器低通部分优化跟踪所期望频率特性的线性相位FIR低通滤波器的单位冲激响应，从而获得这个多变量无约束优化问题的初值，优化范围已经缩小，然后在频域中采用DFP法求解多变量无约束优化问题。通过这种时域与频域结合设计的线性相位IIR数字滤波器低通部分相位接近线性，而且相位线性度可以通过频域设计时权重系数调节，同时该滤波器阶数相对于具有同等幅频特性效果的FIR滤波器阶数降低超过一半，这种近似线性相位滤波器进一步降低了整个数字滤波器的运算量，而且保证了整个频段的幅频特性。

本发明与现有技术相比的优点在于：本发明兼具了常规IIR滤波和FIR滤波器不能同时具有的低阶小运算量和近似线性相位的优点，而且群延时小；在机抖激光陀螺IMU中应用，降低了数字滤波器对IMU硬件的要求，提高了IMU测量实时性，保证了IMU时延补偿精度，为IMU高精度实时测量奠定了基础。

附图说明

图1为本发明的具体实施的POS系统组成框图；

图2为本发明的算法应用的硬件组成框图；

图3为本发明的线性IIR滤波器设计的流程框图；

图4为本发明的线性相位IIR数字滤波器低通部分时域设计原理图；

图5为本发明的递推最小二乘估计的解算流程图。

具体实施方式

如图1所示，为本发明方法实施对象IMU所在的POS系统，IMU包括了激光陀螺仪模块、石英加速度计模块、数据采集电路模块和数据处理模块。本发明设计的滤波器是数据预处理模块的主要组成部分。如图2所示为滤波器所在硬件系统平台，图中标明了信号流向，滤波器程序将在图2中DSP最小子系统中实现。如图3所示，具体实施方法如下：

(1)通过4KHz采集机抖激光陀螺输出信号，采用功率谱密度分析方法，确定机抖激光陀螺抖动噪声峰值频率f_n和强度A_n以及机抖激光陀螺抖频波动范围f_range。线性相位IIR数字滤波器带阻部分和低通部分的指标参数：带阻部分阻带衰减参数A_stop1＝0.7A_n，中心频率f₀＝f_n；低通部分阻带衰减参数A_stop2＝60dB，低通部分通带纹波A_pass2＝0.01dB，以IMU时延要求的波动范围的0.9倍为滤器低通部分群延时波动范围G_dw；线性相位IIR数字滤波器整体阻带衰减参数A_stop应满足通带内信号强度超过阻带内噪声强度一个数量级以上；

(2)采用通带最平坦准则设计IIR滤波器带阻部分H_s(z)，通过选择阶数为二阶或三阶的带阻部分，并选用对称于f₀的带阻参数，使带阻部分的相位非线性相对于低通部分低一个数量级以上；

(3)在机抖激光陀螺的抖频调节范围为3f_range～5f_range的情况下测试H_s(z)。以4KHz采样频率采集H_s(z)滤波后的机抖激光陀螺2s信号，采集三组分析其中功率谱密度，求其均值，如果H_s(z)阻带衰减不到0.7A_n，将f₀上移1Hz重新进行(2)(3)；再将H_s(z)中心频率移动3f_range～5f_range，如果H_s(z)的阻带衰减不到0.7A_s，将f₀上移0.2Hz重新进行步骤(2)、步骤(3)；

(4)线性相位IIR数字滤波器低通部分时域设计。如图4所示线性相位IIR数字滤波器低通部分时域设计原理图，根据等纹波设计准则设计线性相位FIR低通滤波器H_f(z)；以线性相位IIR数字滤波器低通部分H(z)输入输出组成状态方程量测矩阵，建立最小二乘参数估计器；然后在高斯白噪声激励下，利用最小二乘估计器使线性相位IIR数字滤波器低通部分优化跟踪上述线性相位FIR低通滤波器，获得时域设计的线性相位IIR数字滤波器低通部分，具体步骤如下：

(a)以待求系数a₁，a₂，…，a_n，b₀，b₁…b_n建立H_l(z)的模型：

$>H\left(z\right)=\frac{b_0+b_1{z}^{-1}+b_2{z}^{-2}+...+b_n{z}^{-n}}{1+a_1{z}^{-1}+a_2{z}^{-2}+...+a_n{z}^{-n}}---\left(1\right)$>

采用等纹波设计准则设计线性相位FIR低通滤波器其中h₀，h₁，…，h_N-1为线性相位FIR低通滤波器传递函数系数，FIR低通滤波器阶数约为H_l(z)阶数的一倍，通带纹波应比H_l(z)大一个数量级以上，阻带衰减应比H_l(z)少一个数量级；

(b)建立最小二乘估计方程。以a₁，a₂，…，a_n，b₀，b₁，b₂…b_n组成最小二乘估计的估计状态X(k)＝[a₁，a₂，…，a_n，b₀，b₁，…，b_n]^T，以H_l(z)输入x(k)和输出y(k)组成最小二乘估计量测矩阵H(k)＝[-y(k-1)，…，-y(k-n)，x(k)，…，x(k-n)]，建立最小二乘估计方程：

Z(k)＝H(k)X(k)+V(k)    (2)

其中Z(k)为k时刻量测量，H(k)为k时刻状态转移矩阵，X(k)为k时刻系统状态量，V(k)为k时刻量测噪声矩阵；

采用带遗忘因子的最小二乘递推方法，算法如下：

$>\hat{X}\left(k\right)=\hat{X}(k-1)+K\left(k\right)[Z\left(k\right)-H\left(k\right)\hat{X}(k-1)]$>

K(k)＝P(k-1)H^T(k)[H(k)P(k-1)H^T(k)+r(k)]^-1    (3)

$>P\left(k\right)=\frac{1}{r\left(k\right)}[I-K\left(k\right)H\left(k\right)]P(k-1)$>

其中为k时刻系统状态估计量，为k-1时刻系统状态估计量，K(k)为k时刻增益矩阵，Z(k)为k时刻量测量，H(k)为k时刻状态转移矩阵，P(k)为k时刻系统估计量的协方差矩阵，P(k-1)为k-1时刻系统估计量的协方差矩阵，r(k)为k时刻遗忘因子，I为单位矩阵，最小二乘估计遗忘因子采用常值遗忘因子r(k)＝λ，为兼顾最小二乘收敛速度和灵敏度取0.980＜λ＜0.998。最小二乘解算流程如图5所示；

(c)利用(b)中的最小二乘算法在高斯白噪声激励下，取最小二乘估计初值P(0)＝rI，r＞0，使H_l(z)逼近(a)中的线性相位FIR低通滤波器H_f(z)，以获得具有线性相位的H_l(z)，其系数为a₁(0)，a₂(0)，…，a_n(0)，b₀(0)，b₁(0)，b₂(0)…b_n(0)；

(5)线性相位IIR数字滤波器低通部分频域设计。建立线性相位IIR数字滤波器低通部分与两个待逼近FIR低通滤波器单位冲激响应的幅值和相位误差的评价函数，以时域设计结果为初值，采用DFP法对J进行优化，以获得具有所期望频率特性的线性相位IIR数字滤波器；然后将优化后的IIR滤波器用于机抖激光陀螺抖频抑制，测试线性相位IIR数字滤波器低通部分滤波效果。具体步骤如下：

(a)在频域内建立两个线性相位FIR低通滤波器与H_l(z)的幅值误差和相位误差的评价函数：

$>J_{M1}=\underset{i=l_{m0}}{\overset{l_{m1}}{\Sigma }}{\left[\right|H_{\mathrm{dm}}\left(e^{{\mathrm{j\omega }}_i}\right)|-|H_l\left(e^{{\mathrm{j\omega }}_i}\right)\left|\right]}^2---\left(4\right)$>

$>J_{M2}=\underset{i=l_{m1}}{\overset{l_{m2}}{\Sigma }}{\left[\right|H_{\mathrm{dm}}\left(e^{{\mathrm{j\omega }}_i}\right)|-|H_l\left(e^{{\mathrm{j\omega }}_i}\right)\left|\right]}^2---\left(5\right)$>

$>J_{M3}=\underset{i=l_{m2}}{\overset{l_{m3}}{\Sigma }}{\left[\right|H_{\mathrm{dm}}\left(e^{{\mathrm{j\omega }}_i}\right)|-|H_l\left(e^{{\mathrm{j\omega }}_i}\right)\left|\right]}^2---\left(6\right)$>

$>J_P=\underset{i=l_{p0}}{\overset{l_{p1}}{\Sigma }}{[\angle H_{\mathrm{dp}}\left(e^{j{\omega }_i}\right)-\angle H_l\left(e^{j{\omega }_i}\right)]}^2---\left(7\right)$>

其中为线性相位IIR数字滤波器低通部分的传递函数，H_dm(e^jω)和H_dp(e^jω)为两个线性相位FIR低通滤波器传递函数；J_M1，J_M2，J_M3分别为通带，过渡带和阻带的幅值误差的评价函数，J_P为相位误差评价函数；l_m0，l_m1，l_m2，l_m3为幅频误差评价区间的分界点，l_p0，l_p1为通带相频特性的评价区间点；

(b)建立线性相位IIR数字滤波器低通部分的最优准则：

J＝aJ_P1+β₁J_M1+β₂J_M2+β₃J_M3    (8)

评价函数J_M1，J_M2，J_M3和J_P有相当大差距，为使权系数直观反映各部分的相对权重，取

$>\alpha =K_p\frac{J_{M1}\left({\hat{\theta }}_0\right)}{J_{P1}\left({\hat{\theta }}_0\right)},$>β₁＝1

                                (9)

$>{\beta }_2=K_{M2}\frac{J_{M1}\left({\hat{\theta }}_0\right)}{J_{M2}\left({\hat{\theta }}_0\right)},$>$>{\beta }_3=K_{M3}\frac{J_{M1}\left({\hat{\theta }}_0\right)}{J_{M3}\left({\hat{\theta }}_0\right)}$>

其中，X(n)＝[a₁(0)，a₂(0)，…，b₀(0)，b₁(0)，b₂(0)…b_n(0)]^T为时域设计得到的H_l(z)的时域初值，K_p，K_M2，K_M3是直观反映IIR滤波器幅相特性各部分权重的权系数；

(c)设计待逼近的两个FIR低通滤波器H_dm(e^jω)和H_dp(e^jω)。FIR低通滤波H_dm(e^jω)具有待求线性相位IIR数字滤波器低通部分期望的幅频特性，其阶数为H_l(z)阶数的3～5倍，FIR低通滤波H_dp(e^jω)具有所期望的相频特性，其阶数不超过H_l(z)阶数的2倍；

(d)根据H_l(z)幅值和相位权重，确定J权系数K_p，K_M2，K_M3，以a₁(0)，a₂(0)，…，a_n(0)，b₀(0)，b₁(0)，b₂(0)…b_n(0)为初值利用DFP法求目标函数J极小值，以获得具有所期望频率特性的线性相位IIR数字滤波器低通部分分子系数a_l1(n)，a_l2(n)，…，a_ln(n)和分母系数b_l0(n)，b_l1(n)…b_ln(n)，传递函数$>H_l\left(z\right)=\frac{b_{l0}\left(n\right)+b_{l1}\left(n\right)z^{-1}+b_{l2}\left(n\right)z^{-2}+...+b_{\ln }\left(n\right)z^{-n}}{1+a_{l1}\left(n\right)z^{-1}+a_{l2}\left(n\right)z^{-2}+...+a_{\ln }\left(n\right)z^{-n}};$>

(e)通过单位冲激响应测试法，测试优化后的H_l(z)，并将优化后的H_l(z)用于机抖激光陀螺抖频抑制，以4KHz采集2s滤波后信号，进行功率谱密度，如果H_l(z)的通带纹波大于A_pass2，同时减小K_p，K_M2，K_M3至原来的0.9倍，如果H_l(z)的群延时波动大于G_dw，增大权系数K_p至原来的1.05倍，如果机抖激光陀螺抖动噪声衰减低于A_stop1，增大权系数K_M3至原来的1.05倍，重新进行(6)；

(7)利用公式(10)求线性相位IIR数字滤波器H(z)：

H(z)＝H_s(z)H_l(z)    (10)

将H(z)用于机抖激光陀螺噪声抑制，并以4KHz采集2s滤波后信号进行功率谱分析，如果H(z)阻带内噪声衰减低于A_stop，将A_stop1和A_stop2增大5dB，重新进行步骤(2)-步骤(7)，求得H(z)的分子系数b₀(n)，b₁(n)，b₂(n)…b_n(n)和分母系数a₁(n)，a₂(n)，…，a_n(n)，最终获得线性相位IIR数字滤波器传递函数$>H\left(z\right)=\frac{b_0\left(n\right)+b_1\left(n\right)z^{-1}+b_2\left(n\right)z^{-2}+...+b_n\left(n\right)z^{-n}}{1+a_1\left(n\right)z^{-1}+a_2\left(n\right)z^{-2}+...+a_n\left(n\right)z^{-n}},$>完成线性相位IIR数字滤波器设计；在具有相同机抖激光陀螺噪声抑制效果的情况下本发明设计的线性相位IIR数字滤波器比线性相位FIR低通滤波器阶数低，运算量减小，群延时可减少50％以上，为IMU高精度实时测量奠定了良好的基础。

本发明未详细阐述部分属于本领域公知技术。

显然，对于本领域的普通技术人员来说，参照上文所述的实施例还可能做出其它的实施方式。本发明中的实施例都只是示例性的、而不是局限性的。所有的在本发明的权利要求技术方案的本质之内的修改都属于其所要求保护的范围。