CN106153046B - 一种基于自适应Kalman滤波的陀螺随机噪声AR建模方法 - Google Patents

Abstract

本发明公开了一种基于自适应Kalman滤波的陀螺随机噪声AR建模方法，包括步骤如下：(1)对陀螺仪随机噪声进行数据预处理；通过数字微分保证随机噪声数据满足平稳随机过程的要求；(2)将AR模型的参数做为系统状态量；(3)建立系统的观测方程；(4)建立系统的状态方程；(5)采用未知定常的观测噪声统计估值器估计观测噪声的均值和方差；(6)在估值器估计的观测噪声的均值和方差基础上运用自适应Kalman滤波器对AR模型的参数进行估计。本发明可以有效减少采样样本数和采样时间，建模速度快；当新的噪声样本到来时，自适应Kalman滤波器能够保证建立的AR模型参数可以及时进行更新，具有快速实时建模的优点。

Description

一种基于自适应Kalman滤波的陀螺随机噪声AR建模方法

技术领域

本发明涉及一种基于自适应Kalman滤波的陀螺随机噪声AR建模方法，属于捷联惯性导航技术领域。

背景技术

陀螺仪输出误差的建模和补偿是提高捷联惯性导航系统精度的重要方法。陀螺仪的输出误差可以分为确定性误差和随机误差。确定性误差可以通过实验标定的方法予以大部分消除。对陀螺仪随机误差的处理一般是从信号处理角度，通过数学工具进行建模和补偿，其中时间序列建模方法得到了广泛的应用。

常用的时间序列模型有ARMA(Auto Regressive Moving Average，自回归滑动平均模型)、AR(Auto Regressive，自回归模型)和MA(MovingAverage，滑动平均模型)三种。由Wold定理和Kolmogorov-Szego定理可知，这三种模型可以相互转化，都具有普遍适用性。但如果选择的模型合适，就可以用低阶的模型获得相对高的建模精度，从而使系统更加简单。

实际工程中不同型号和批次的陀螺仪其随机噪声特性并不相同。有些型号批次陀螺仪的噪声数据其自相关和偏相关特性采用AR模型建模更加适合。传统陀螺随机噪声AR建模方法存在着算法收敛速度慢，所需样本多等缺点。

发明内容

为解决传统陀螺随机噪声建模方法的不足，本发明提出一种基于自适应Kalman滤波技术的陀螺随机噪声AR模型建模方法。相比传统方法，本发明方法建模速度快，所需样本少，建模精度高。

本发明采用的技术方案为：

一种基于自适应Kalman滤波技术的陀螺随机噪声建模方法，包括以下步骤：

(1)首先对陀螺仪随机噪声进行平稳随机检验。采用数字微分方法保证随机噪声数据满足平稳随机过程的要求。

(2)将AR模型的参数做为系统状态量。以AR(2)模型为例：

z(k)＝a₁z(k-1)+a₂z(k-2)+ε(k)；

式中a₁、a₂是AR模型的自回归系数；z(k-1)是一阶自回归项；z(k-2)是二阶自回归项；ε(k)是均值为0零，方差未知(常数)的白噪声；a₁，a₂是待定的模型参数。将a₁，a₂做为系统的状态量X：

X＝[a₁(k)，a₂(k)]^T；

(3)建立系统观测方程：

z(k)＝H(k)X(k)+v(k)；

式中H(k)是系统观测阵：H(k)＝[z(k-1)，z(k-2)]；v(k)是系统观测噪声，它主要由AR(2)模型的白噪声ε(k)构成。

(4)建立系统状态方程：

X(k+1)＝X(k)；

(5)采用未知定常的观测噪声统计估值器估计观测噪声的均值和方差，在此基础上运用自适应Kalman滤波器对系统状态量X，也即AR模型的参数进行估计。

进一步，上述步骤(1)所述数据预处理包括以下步骤：

(a)周期性数据检验和滤除；

(b)平稳随机检验；

(c)分析数据的自相关和偏相关特性；当自相关系数函数呈现拖尾性而偏相关系数函数呈现截尾性时，则选择AR模型建模。且偏相关系数函数截断处就是AR模型的阶数。

且步骤(b)所述的平稳随机检验还包括差分处理过程。

进一步，用下标“k”表示Kalman滤波的迭代次数，则步骤(5)所述的自适应Kalman滤波器的迭代过程如下所示：

状态一步预测：

式中

是Kalman滤波器根据：k-1迭代时得到的系统状态量X的估计值

对第k次迭代时系统状态量X估计值的预测，称为状态一步预测值；Φ_k|k-1是系统状态转移阵，当对AR(2)进行参数估计时，它为2×2维的单位矩阵；

是k-1次迭代时得到的状态量X的估计值；

一步预测误差方差阵：

式中P_k|k-1是一步预测误差方差阵；Q是系统激励噪声的方差，Q＝0；P_k-1|k-1是k-1次迭代时得到的估计误差方差阵；

滤波增益矩阵：

式中K_k是Kalman滤波器第k次迭代时的滤波增益矩阵；P_k|k-1是一步预测误差方差阵；H_k是第k次迭代时的系统观测矩阵；

是第k-1次迭代时得到的系统观测噪声方差的估计；

新息：

式中

是第k次迭代时的新息；z_k是第k次迭代时的系统观测量，也即陀螺随机噪声的采样；H_k是第k次迭代时的系统观测矩阵；

是状态一步预测；

是第k-1次迭代时得到的系统观测噪声均值的估计；

状态估计：

式中

是Kalman滤波在迭代k次时得到的状态量X的估计值；

是状态一步预测；K_k是第k次迭代时的滤波增益矩阵；

是第k次迭代时的新息；

估计误差方差阵：

P_k|k＝[I_n-K_kH_k]P_k|k-1；

式中P_k|k是第k次迭代时的估计误差方差阵；I_n是n×n维的单位矩阵，当对AR(2)进行参数估计时，n＝2；K_k是第k次迭代时的滤波增益矩阵；H_k是第k次迭代时的系统观测矩阵；P_k|k-1是一步预测误差方差阵；

未知定常的观测噪声均值估计：

式中

为滤波器第k次迭代时观测噪声均值的估计值；z_k是第k次迭代时的系统观测量；

是第k次迭代时状态量的一步预测，H_k是第k次迭代时的系统观测阵；

未知定常的观测噪声方差估计：

式中

为滤波器第k次迭代时观测噪声方差的估计值；

是第k次迭代时的新息；P_k|k-1是第k次迭代时的一步预测误差方差阵，H_k是第k次迭代时的系统观测阵。

只要AR模型阶数选取合适，通过自适应Kalman滤波器迭代，则可迅速得到AR模型参数的最优估计。本发明的有益之处在于：本发明可以有效减少采样样本数和采样时间，建模速度快。当新的噪声样本到来时，自适应Kalman滤波器能够保证建立的AR模型参数可以及时进行更新，具有快速实时建模的优点。

附图说明

图1为某国产光纤陀螺的原始噪声数据；

图2为该陀螺原始噪声经一阶差分后满足平稳随机过程要求的数据；

图3为该陀螺随机噪声数据的自相关特性分析；

图4为该陀螺随机噪声数据的偏相关特性分析；

图5为发明方法对该陀螺随机噪声进行AR(2)建模结果a₁；

图6为发明方法对该陀螺随机噪声进行AR(2)建模结果a₂；

图7是为本发明方法步骤过程示意图。

具体实施方式

以下通过具体实施例对本发明作具体的介绍。

下面以对某国产光纤陀螺仪输出随机噪声的AR(2)建模为例，首先对陀螺随机噪声的采样数据进行预处理，保证其满足平稳随机过程的要求。数据预处理包括3步：1)周期性数据检验和滤除；2)平稳随机检验；3)根据自相关和偏相关特性选择合适的时间序列模型建模。

(a)周期性数据检验和滤除

如果信号的采样含有周期性或准周期性的数据，则在功率谱上会出现尖峰，尖峰处对应的角频率也即周期项的角频率，而不含周期项的随机数据的功率谱没有明显尖峰。信号功率谱和自相关函数的定义如下：

如果陀螺的噪声数据含有周期为T的周期项，则使用T步差分是去除周期项的一种有效办法。

(b)平稳随机检验

经过趋势项和周期项检验和滤除后的陀螺噪声数据可以认为是随机序列。但还需检验其平稳性。如不满足平稳随机序列的要求，也需通过差分处理的方法使之平稳化。实际工程应用中一般用非参数法，例如游程检验法来检验时间序列的平稳性。游程是指在保持随机序列原有顺序的情况下，随机序列中具有相同符号的序列。这种检验法里把观测值分成两类：例如

和

是序列均值。分别用“+”和“-”表示，称为正负符号。游程检验法是以正负符号的数目N₁和N₂做为统计检验量，根据正态分布表来确定检验的接受域。例如，以置信度95％为可接受，即假设检验的显著性水平α＝0.05。按照双边假设检验(2σ)查找标准正态分布表，则检验统计量U在[-1.96，1.96]范围时为假设接受域。可以认为该序列是平稳的。检验统计量的计算公式为：

式中U是检验统计量；γ是游程数；μ_γ是游程的期望数；σ_γ是游程的标准差；计算方法如下：

式中N₁、N₂分别是游程中正、负符号的数目；μ_γ是游程的期望数；σ_γ是游程的标准差。不管多复杂的非平稳随机过程，经过一定次数的差分后，总可以变成平稳随机过程。而对于陀螺的随机噪声数据来说，一般经过一阶或二阶差分处理，即可基本满足平稳性要求。

(c)根据现代时间序列分析理论，时间序列模型选择的依据在于平稳随机数据的自相关系数函数和偏相关系数函数的截尾和拖尾特性的不同，具体区分如下表所示：

表1三种时间序列模型的自相关和偏相关特性

其中，自相关系数函数p(h)和偏自相关系数函数

的数学定义如下：

式中x(k)、x(k-h)分别代表第k和第k-h个样本数据；“Cov”代表自协方差函数；“Var”代表方差函数；

和

分别是利用x(n+1)，x(n+2)，…，x(n+k-1)对x(n)和x(n+k)做的最佳线性估计。

附图1是地面静态实验测得的某国产光纤陀螺的原始噪声数据(x轴，共3600s)。采样频率100Hz，单位：°/s，陀螺坐标轴指向：x-东，y-北，z-天。首先按对其进行数据预处理，通过一阶差分，最终得到满足平稳随机过程要求的数据(见附图2)。分析其自相关系数函数p(h)和偏相关系数函数

定义如式(4)所示。

从附图3、附图4中可以看出，该陀螺仪随机噪声数据的自相关系数函数明显呈现拖尾性，而偏自相关系数函数呈现截尾特性，且在“2”后截断。从表1可知应优先选用AR模型建模。且偏相关系数函数的截断处“2”就是AR模型的阶数。因此选用二阶AR模型进行建模。

传统AR模型建模方法有最小二乘法、Yule-Walker法等。但这些传统方法需要的样本数较多，算法收敛速度较慢。为了解决传统陀螺随机噪声AR建模方法中存在的问题，本发明提出了一种新的基于自适应Kalman滤波的陀螺随机噪声AR模型建模方法，如图7所示。

首先对陀螺仪随机噪声进行平稳随机检验，通过数字微分保证随机噪声数据满足平稳随机过程的要求后，将AR模型的参数做为系统状态量；

假设陀螺仪随机噪声的AR模型阶数为AR(2)：

z(k)＝a₁z(k-1)+a₂z(k-2)+ε(k) (5)

式中的ε(k)是均值为0方差未知(常数)的白噪声。取状态量：X(k)＝[a₁(k)，a₂(k)^T，则可以把上式转化为系统观测方程：

z(k)＝H(k)X(k)+v(k) (6)

式中的H是系统观测阵：H(k)＝[z(k-1)，z(k-2)]，v(k)＝ε(k)。考虑到当陀螺仪的随机噪声特性相对稳定时，当满足一定的样本数后，所建立的AR模型的参数估计值应当稳定收敛于真值，不再随样本数量的增加而改变。即

式中上标“^”表示参数的估计值；可以建立系统的状态方程：

X(k+1)＝X(k) (7)

从上式可以看出，系统噪声w(k)为0，状态转移阵Φ＝I。观测噪声v(k)主要由AR模型中的白噪声ε(k)引起。有：

式中“E[]”代表均值函数；r_k、R_k是系统观测噪声v(k)的均值阵和方差阵；δ_kj是克罗内克函数。

实际系统由于不可避免存在建模误差，所以r_k不一定为0。又由于AR模型中白噪声ε(k)的方差特性是未知定常的，因此R_k也是未知定常型的。对于未知定常型的观测噪声，通常利用Sage-Husa的次优无偏MAP(Maximum A Posteriori极大后验)噪声统计估值器估计未知定常观测噪声的统计特性。对应的Kalman滤波器称为自适应Kalman滤波。未知定常观测噪声统计估值器如下式：

对该光纤陀螺随机噪声数据采用发明提出的自适应Kalman滤波方法进行建模，滤 波初值选取

建模过程如附图5-6中 的实曲线所示。图中虚直线是采用传统“Y-W”法(1小时采样数据360000样本点)对该陀螺随 机噪声数据进行AR(2)建模的结果：

z(k)＝-0.671z(k-1)-0.334z(k-2)+ε(k) (10)

在滤波迭代时，可以根据系统的建模精度要求设定阈值θ。如有连续10次运行结果中各参数估计值

的最大值和最小值之间的差值都小于阈值θ，则做为迭代的退出条件。这里选取θ_a1＝θ_b1＝θ_b2＝0.001。程序在运行980次迭代后退出，状态量估计

收敛于：

对应的AR(2)模型为：

z(k)＝-0.669z(k-1)-0.338z(k-2)+ε(k) (12)

如果仍采用传统Y-W法，但只对前980个陀螺随机噪声数据进行建模，其建模结果为：

z(k)＝-0.682z(k-1)-0.349z(k-2)+ε(k) (13)

分别将发明的基于自适应Kalman滤波的陀螺随机噪声建模结果(980点数据)，与传统Y-W法建模的结果(360000点数据和980点数据)列于下表：

表2发明方法与传统Y-W法的建模结果对比

从表2可以看出，只用980点陀螺随机噪声数据进行建模时，发明的基于自适应Kalman滤波的陀螺随机噪声AR建模方法的建模结果(式(12))与传统的Y-W法建模结果(式(13))相比，更接近1小时完整数据的建模结果(式(10))。说明发明方法在快速建模方面具有优势。

此外，当有新的陀螺随机噪声样本数据到来时，传统的AR建模方法必须对全部历史数据重新计算一次才能得到新的模型参数。而自适应Kalman滤波器则可以利用新的噪声观测数据对AR模型的参数估值进行修正。既解决了重复计算问题，也可以使建立的模型能够实时跟踪陀螺随机噪声特性的变化，建模的精度高。

全部历史数据重新计算一次才能得到新的模型参数。而自适应Kalman滤波器则可以利用新的噪声观测数据对AR模型的参数估值进行修正。既解决了重复计算问题，也可以使建立的模型能够实时跟踪陀螺随机噪声特性的变化，建模的精度高。

发明方法在实际应用时，选用的AR模型的阶数并不唯一，本发明方法同样适用于对陀螺随机噪声其它阶AR模型的建模，例如AR(1)，AR(3)等，前提是该陀螺随机噪声数据的自相关系数函数和偏相关系数函数的特性满足该模型的要求；但不同阶数的AR模型在进行自适应Kalman滤波建模时，滤波参数的维数需要随之改变；例如如果陀螺随机噪声数据适合AR(3)模型建模，则系统状态量将变为：X＝[a₁(k)，a₂(k)，a₃(k)]^T；系统观测阵将变为：H(k)＝[z(k-1)，z(k-2)，z(k-3)]；一步预测误差方差阵P_k|k-1、估计误差方差阵P_k|k和滤波增益矩阵K_k也将由原来的2阶变为3阶。

上述实施例不以任何形式限制本发明，凡采用等同替换或等效变换的方式所获得的技术方案，均落在本发明的保护范围内。

Claims (3)

1.一种基于自适应Kalman滤波的陀螺随机噪声AR建模方法，其特征在于，包括步骤如下：

(1)首先对陀螺仪随机噪声进行数据预处理，通过数字微分保证随机噪声数据满足平稳随机过程的要求；所述数据预处理包括以下步骤：

(a)周期性数据检验和滤除，包括定义信号功率谱和自相关函数如下：

使用T步差分是去除周期项；

(b)平稳随机检验；以正负符号的数目N₁和和N₂做为统计检验量，根据正态分布表来确定检验的接受域，以置信度95％为可接受，即假设检验的显著性水平α＝0.05，按照双边假设检验(2σ)查找标准正态分布表，则检验统计量U在[-1.96,1.96]范围时为假设接受域，认为该序列是平稳的，且检验统计量的计算公式为：

式中U是检验统计量；γ是游程数；μ_γ是游程的期望数：σ_γ是游程的标准差；计算方法如下：

(c)分析数据的自相关和偏相关特性；当自相关系数函数呈现拖尾性而偏相关系数函数呈现截尾性时，则选择AR模型建模，且偏相关系数函数截断处就是AR模型的阶数；

定义自相关系数函数ρ_h和偏自相关系数函数

的数学定义如下：

式中X(k)、X(k-h)分别代表第k和第k-h个样本数据；“Cov”代表自协方差函数；“Var”代表方差函数；

和

分别是利用X(n+1),x(n+2),…,X(n+k-1)对X(n)和X(n+k)做的最佳线性估计；

(2)将AR模型的参数做为系统状态量，以AR(2)模型为例：z(k)＝a₁z(k-1)+a₂z(k-2)+ε(k)；

式中a₁、a₂是AR模型的自回归系数；z(k-1)是一阶自回归项；z(k-2)是二阶自回归项；ε(k)是均值为0，方差未知的白噪声；a₁、a₂是待定的模型参数，将a₁、a₂做为系统的状态量X：

X＝[a₁(k)，a₂(k)]^T；

(3)建立系统观测方程：

z(k)＝H(k)X(k)+v(k)；

式中H(k)是系统观测阵：H(k)＝[z(k-1)，z(k-2)]；v(k)是系统观测噪声，它由AR(2)模型的白噪声ε(k)构成；

(4)建立系统的状态方程：

X(k+1)＝X(k)；

从上式可以看出，当系统噪声w(k)＝0，状态转移阵Φ＝1，观测噪声V(k)由AR模型中的白噪声ε(k)引起，有：

式中“E[]”代表均值函数、r_k、R_k是系统观测噪声V(k)的均值阵和方差阵，δ_kj是克罗内克函数；

(5)采用均值、方差未知且定常的系统观测噪声统计估值器估计观测噪声的均值和方差；在此基础上运用自适应Kalman滤波器对系统状态量X，也即AR模型的参数进行估计；用下标“k”表示Kalman滤波的迭代次数，则步骤(5)所述的自适应Kalman滤波器的迭代过程如下所示：

状态一步预测：

x^|_k|k-1＝Φ_k|k-1x^_k-1|k-1＝x^_k-1|k-1;

式中x^_k-1|k-1是Kalman滤波器根据：k-1迭代时得到的系统状态量X的估计值对第k次迭代时系统状态量X估计值的预测，称为状态一步预测值；Φ_k|k-1是系统状态转移阵，当对AR(2)进行参数估计时，它为2×2维的单位矩阵；是k-1次迭代时得到的状态量X的估计值；

一步预测误差方差阵：

P_k|k-1＝Φ_k|k-1P_k-1|k-1Φ_k|k-1 ^T+Γ_k|k-1Q_k-1Γ_k|k-1 ^T＝P_k-1|k-1；

式中P_k|k-1是一步预测误差方差阵；Q是系统激励噪声的方差，Q＝0；P_k-1|k-1是k-1次迭代时得到的估计误差方差阵；

滤波增益矩阵：

K_k＝P_k|k-1H_k ^T[H_kP_k|k-1H_k ^T+R^_k-1]^-1

式中K_k是Kalman滤波器第k次迭代时的滤波增益矩阵；P_k|k-1是一步预测误差方差阵；H_k是第k次迭代时的系统观测矩阵；R^_k-1是第k-1次迭代时得到的系统观测噪声方差的估计；

新息：

ε^_k＝z_k-H_kx^_k|k-1-r^_k-1；

式中ε^_k是第k次迭代时的新息；z_k是第k次迭代时的系统观测量，也即陀螺随机噪声的采样；H_k是第k次迭代时的系统观测矩阵；x^_k|k-1是状态一步预测；r^_k-1是第k-1次迭代时得到的系统观测噪声均值的估计；

状态估计：

x^_k|K＝x^_k|k-1+K_kε ^_k ；

式中x^_k|K是Kalman滤波在迭代k次时得到的状态量X的估计值；x^_k|k-1是状态一步预测；K是第k次迭代时的滤波增益矩阵；ε^_k是第k次迭代时的新息；

估计误差方差阵：

P_k|k＝[I_n-K_kH_k]P_k|k-1；

式中P_k|k是第k次迭代时的估计误差方差阵；I_n是n×n维的单位矩阵，当对AR(2)进行参数估计时，n＝2；K_k是第k次迭代时的滤波增益矩阵；H_k是第k次迭代时的系统观测矩阵；P_k|k-1是一步预测误差方差阵；

未知定常的观测噪声均值估计：

r^_k＝(1-1/k)r^_k-1+1/k(z_k-H_kx^_k|k-1)；

式中r^_k为滤波器第k次迭代时观测噪声均值的估计值；z_k是第k次迭代时的系统观测量；x^_k|k-1是第k次迭代时状态量的一步预测，H_k是第k次迭代时的系统观测阵；

未知定常的观测噪声方差估计：

R^_k＝(1-1/k)R^_k-1+1/k(ε ^_kε ^_k ^TH_kP_k|k-1H_k ^T)；

式中R^_k为滤波器第k次迭代时观测噪声方差的估计值；ε^_k是第k次迭代时的新息；P_k|k-1是第k次迭代时的一步预测误差方差阵，H_k是第k次迭代时的系统观测阵。

2.根据权利要求1所述的一种基于自适应Kalman滤波的陀螺随机噪声AR建模方法，其特征在于，步骤(b)所述的平稳随机检验还包括差分处理过程。

3.根据权利要求1或2所述的一种基于自适应Kalman滤波的陀螺随机噪声AR建模方法，其特征在于，所述的自适应Kalman滤波器的迭代过程中还根据系统的建模精度要求设定有阈值θ。

Images