CN102650527A - 一种基于时间序列分析消噪的光纤陀螺温度补偿方法 - Google Patents

Abstract

一种基于时间序列分析消噪的光纤陀螺温度补偿方法，它有四个步骤：步骤一：设计实验方案，对光纤陀螺进行定点高低温测试实验，利用采集软件进行数据采集。步骤二：对陀螺零偏数据进行时间序列分析，建立光纤陀螺随机误差的数学模型。步骤三：采用卡尔曼滤波算法滤除光纤陀螺零偏数据中的随机噪声。步骤四：利用卡尔曼滤波消噪后的数据对光纤陀螺温度漂移误差模型结构进行辨识，并对辨识后模型进行解算参数。本发明经过时间序列分析和卡尔曼滤波消噪处理，并进行温度漂移误差模型结构和参数辨识，建立光纤陀螺静态温度漂移误差多项式模型。该方法完全符合工程上的实时补偿要求，在航空航天导航技术领域里具有较好的实用价值和广阔的应用前景。

Description

一种基于时间序列分析消噪的光纤陀螺温度补偿方法

（一）技术领域

本发明涉及一种基于时间序列分析消噪的光纤陀螺温度补偿方法，它详细阐述了时间序列分析消除噪声以及温度补偿的方法，这种方法应用于光纤陀螺输出和温度关系的辨识中，可以有效地辨识出这种关系，建立精确的模型，提高温度补偿的有效性。该技术属于航空、航天导航技术领域。

（二）背景技术

光纤陀螺作为一种发展极为迅速的新型惯性角速度传感器，以其特有的技术和性能优势，如全固态结构、可靠性高、寿命长；启动速度快，响应时间短；测量范围大，动态范围宽；抗冲击、振动，耐化学腐蚀；体积小、重量轻、成本低；适合大批量生产等，已经广泛用于各工程技术领域。工程化的光纤陀螺，为适应不同领域的应用，一般要具有较宽的工作温度范围（通常-40℃～60℃）。构成光纤陀螺的主要器件如光纤线圈、集成光学器件、光源、耦合器等对温度较为敏感，当环境温度发生变化时，在陀螺的输出信号中将产生非互易相位误差，导致光纤陀螺零偏和标度因数的不稳定，并最终影响光纤陀螺在不同温度下的检测精度。因此，进行光纤陀螺温度特性的研究并对其实施温度漂移误差补偿，以提高陀螺在不同温度条件下的适用性，是光纤陀螺走向实用化的必要环节。

国内外近年来对光纤陀螺温度漂移误差补偿进行的各项工作，主要包括光纤陀螺机理结构的改善、硬件温控措施及软件补偿等内容。对光纤陀螺进行结构改进或对其进行硬件温控受到成本及技术等多方面制约，带有硬件温控的光纤陀螺结构复杂、功耗大，启动时间长，不适应快速启动的要求，而软件补偿可改善这一缺陷，更易于实现。补偿光纤陀螺温度漂移误差的方法主要包括下述几个方面：

（1）光纤陀螺结构及组成部件的改进。可选择热膨胀系数和折射率温度系数低的石英材料制作纤芯，采用热传导系数高的光线外涂复合材料使各光纤层之间形成“热短路”等。

（2）改善光纤环绕制技术。为减少热致非互易性，同时考虑到微弯损耗，光纤环通常采用如下四种绕制方法：1）柱形绕法；。2）对称绕法；3）双极对称绕法；4）四级对称绕法。

（3）对光纤陀螺进行温度控制。光纤陀螺中的光学元器件性能随温度的变化，及电路部分各器件发热不同造成的陀螺内部温度场的非线性变化，都会影响陀螺工作的热平衡状态，直接造成光纤陀螺的零偏漂移。使光纤陀螺迅速进入并稳定工作在热平衡状态，会极大地改善陀螺的输出特性。

（4）误差建模及补偿

对于给定的光纤陀螺，在一定的温度变化情况下，产生的非互易相位噪声是确定的。因此，对光纤陀螺温度特性进行实验研究，建立温度漂移误差模型并实施温度补偿是解决光纤陀螺温度漂移问题的有效手段之一。目前常用的模型主要有如下几种：1）线性模型或多项式模型；2）神经网络模型；3）小波网络；4）模糊逻辑；5）受控马氏链模型；6）支持向量机。其中多项式模型及神经网络模型应用最为广泛。

时间序列分析是一种时域分析方法，不仅研究过程的确定性变化，而且更着重于研究过程的随机性变化，是一种对有序的随机数据进行分析和研究，将现在时刻的值用过去时刻的值表示，并预测系统未来值的数据处理方法。时间序列分析法建模就是对所观测到的时间序列数据y₁,y₂,…,y_n拟合出适用的时间序列模型。建模内容包括数据的采集、数据的统计分析（平稳性分析和相关函数分析）与预处理、模型形式的选取、模型参数的估计、模型适用性检验等问题，其中模型参数的选择最关键，其次是模型的适用性检验。对于光纤陀螺的建模一般采用自回归模型（AR模型）、自回归-滑动平均混合模型（ARMA模型）。当输出序列是非平稳时间序列时，一般采用自回归差分滑动平均模型（ARIMA模型）。在时间序列分析中，用相邻观测数据或一定间隔的观测数据的差分方法，能方便地除去序列中的非周期性或周期性趋势项，利于对光纤输出数据进行模型参数辨识和噪声统计特性计算。在求出模型参数和噪声统计特性的基础上，通过卡尔曼滤波/平滑技术可达到良好的消噪效果。

本发明提出了通过大量实验得到的数据，经过时间序列分析方法建立陀螺零偏随机误差模型，通过卡尔曼滤波技术进行消噪处理；并在此基础上进行温度模型结构和参数辨识，建立光纤陀螺温度漂移误差模型，实现温度补偿。该方法能够更有效的辨识出陀螺零偏温度特性，建立精确的模型进行实时补偿，满足工程化要求。

（三）发明内容

本发明提供了一种基于时间序列分析消噪的光纤陀螺温度补偿方法。该方法利用时间序列分析建立光纤陀螺零偏的随机误差模型，在求出模型参数和噪声统计特性的基础上，采用卡尔曼滤波实现光纤陀螺零偏数据中噪声的降低，有效辨识出陀螺零偏温度特性，最后对消噪后的陀螺数据建立温度漂移误差的多项式模型，用最小二乘法求得模型参数，实现光纤陀螺的温度漂移误差补偿。

本发明一种基于时间序列分析消噪的光纤陀螺温度补偿方法，该方法具体步骤如下：

步骤一：设计温度特性测试实验方案，对光纤陀螺进行定点温度实验，采集陀螺数据。温度实验采用静态测试，其中，X，Y，Z轴光纤陀螺分别指向东、北、天。为了研究温度对光纤陀螺零偏的影响，分别在-30℃、-20℃、-10℃、0℃、10℃、20℃、30℃、40℃、50℃、60℃环境温度下（或根据工作温度确定所选温度点），对光纤陀螺进行高低温测试。在每一个温度点保温两小时后测量一小时，并通过采集软件记录光纤陀螺自身的温度和相应光纤陀螺零偏值。

步骤二：对陀螺零偏数据进行时间序列分析，建立光纤陀螺随机误差的数学模型。如图1所示，具体的建模步骤为：

（1）对陀螺测试样本数据进行统计检验和预处理。首先剔除数据中的异常值，其次进行平稳性检验，如发现为非平稳的，应剔除其中确定性的趋势项，再次进行周期性检验，如发现潜周期分量，应剔除其中能量较大的潜周期分量，最后对去除了趋势项和潜周期分量的残差序列进行正态性检验。如果经过检验的陀螺测试数据的残差序列是非正态时间序列，应进行差分处理使之成为正态时间序列，然后建立随机误差模型。具体方法如下。

1）数据异常值剔除

本发明采用拉依达准则判别异常值。假定数据的总体是正态分布，则：

P(|x-μ|>3σ)<0.003    （5-1）

其中：x为随机变量，μ和σ为数据总体的数学期望和标准差，P为满足括号内条件的概率。

设数据为x₁,x₂,…,x_n，均值为

残差为V_k(k=1,2,…,n)，标准差为：

$\mathrm{\&sigma;}=\sqrt{\frac{1}{n-1}\underset{k}{\mathrm{\&Sigma;}}{V}_k^2}=\sqrt{\frac{1}{n-1}[\underset{k}{\mathrm{\&Sigma;}}{\left(x_k\right)}^2-{\left(\underset{k}{\mathrm{\&Sigma;}}\stackrel{\&OverBar;}{x}\right)}^2/n]}---(5-2)$

若某个测量值x_d的残差V_d(1<d<n)满足|V_d|>3σ，则认为x_d是异常值，予以剔除。

2）平稳性检验

时间序列的平稳性是建模的重要前提。设观测序列的值为x_k(k=1,2,…,n)，均值为

用符号“+”表示

“-”表示

在保持数据原有顺序的情况下，游程定义为具有相同符号的序列，这种符号可以把观测值分成两个互相排斥的类。游程检验假设：样本数据出现的顺序没有明显的趋势，就是平稳的。游程太多或太少都被认为存在非平稳性趋势。设：N₁为一种符号出现的总数；N₂为另一种符号出现的总数；r为游程的总数。当N₁或N₂超过15时可以用正态分布来近似，此时选用统计量为：

其中： ${\mathrm{\&mu;}}_r=\frac{2N_1{N}_2}{N}+1;$ ${\mathrm{\&sigma;}}_r=\sqrt{2N_1{N}_2(2N_1{N}_2-N)/N^2(N-1)};$ N=N₁+N₁。

对于显著水平α=0.05的情况，如果|Z|≤1.96(按2σ原则)，可接受原假设，认为序列是平稳的，否则拒绝。

当平稳性检验中判断{x_k}是非平稳时间序列时，应提取{x_k}中所含的非平稳部分，称为趋势项d_k，一般表示为时间k的多项式：

d_k=a₀+a₁k+a₂k²+…+a_mk^m    （5-4）

式中，a₀,a₁,a₂,…,a_m是多项式数学模型的系数，可由回归分析求出。

光纤陀螺输出的非平稳时间序列中，较常见的趋势项为：

d_k=a₀+a₁k    （5-5）

它表明光纤陀螺输出的均值随时间呈线性变化，但一般是极为缓慢的。在光纤陀螺输出非平稳时间序列{x_k}中去除出了趋势项d_k以后，得到时间序列{y_k}：

y_k=x_k-d_k    （5-6）

再对时间序列{y_k}建模。

3）周期性检验

周期性检验用来识别光纤陀螺输出数据中是否包含随机量以外的周期性分量。周期性检验的方法是直接考察由输出数据得到的概率密度函数、自相关函数或功率谱密度函数的图形。

输出数据中如果存在周期性分量，会在概率密度函数图形中反映出来。随机量和周期性分量的概率密度函数图形的不同之处：随机量的概率密度函数图形呈钟形曲线，周期性分量的概率密度函数图形呈盆形曲线；如果输出数据中同时含有随机量和周期性分量，会出现双峰形曲线，峰值处对应的横坐标值为周期性分量的振幅。

输出数据中如果存在周期性分量，也会在自相关函数图形中反映出来。随机量和周期性分量的自相关函数图形的不同之处：随机量的自相关函数图形在时间间隔增大时，是一条衰减的曲线，周期性分量的自相关函数图形不管时间间隔怎样增大，总是一条不衰减的振荡曲线；如果输出数据中同时含有随机量和周期性分量，其自相关函数图形在一定的时间间隔内呈衰减趋势，然后维持不衰减的振荡。

输出数据中如果存在周期性分量，还会在功率谱密度函数图形中反映出来。当输出数据中同时含有随机量和周期性分量时，功率谱密度函数曲线中会有明显的尖峰，尖峰处所对应的横坐标值为周期性分量的频率。

为了寻找隐含周期，识别并提取随机序列中的周期函数项，通常采用周期图分析法。分析出周期性分量后，应予以剔除。剔除方法为：根据前述方法分析出的周期性分量对应的频率，并对输出数据序列进行傅里叶级数展开，找出周期性分量对应的傅里叶系数，由系数构造出周期性分量序列，在原始数据序列中减去周期性分量序列，实现周期性分量的剔除。

4）正态性检验

正态性检验用来判断光纤陀螺输出数据序列是否具有正态分布的特性。对输出序列{x_k}正态性的检验，最基本的是检验{x_k}的三阶矩(偏态系数ξ)与四阶矩(峰态系数v)是否满足正态随机变量的特性。偏态系数ξ与峰态系数v的定义为：

$\mathrm{\&xi;}=E{\left[\frac{x_k-{\mathrm{\&mu;}}_x}{{\mathrm{\&sigma;}}_x}\right]}^3;$ $v=E{\left[\frac{x_k-{\mathrm{\&mu;}}_x}{{\mathrm{\&sigma;}}_x}\right]}^4---(5-7)$

偏态系数ξ反映了随机变量的概率密度函数峰的对称性，峰态系数v反映了随机变量的概率密度函数峰的状态。对于正态随机变量：ξ=0；v=3。

对{x_k}计算ξ和v的估计值：

$\left\{\begin{array}{c}\widehat{\mathrm{\&xi;}}=\frac{1}{n{\widehat{\mathrm{\&sigma;}}}_x^3}\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{(x_k-{\widehat{\mathrm{\&mu;}}}_x)}^3\\ \hat{v}=\frac{1}{n{\widehat{\mathrm{\&sigma;}}}_x^4}\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{(x_k-{\widehat{\mathrm{\&mu;}}}_x)}^4\end{array}\right.---(5-8)$

式中，和

分别是{x_k}的均值和方差的估计值。当算得

时，认为{x_k}是正态序列。如果序列不是平稳的，那么采用若干次差分都可以使序列变成平稳的。

（2）模型形式的选取及参数估计。其具体实现过程如下：

1）选择模型的形式

光纤陀螺的建模采用自回归模型（AR模型）、自回归-滑动平均混合模型（ARMA模型）。自回归模型指任何一个时刻k上的数值y_k可以表示为过去p个时刻上数值的线性组合加上k时刻的白噪声，表示为：

y_k=a₁y_k-1+…+a_py_k-p+r_k    （5-9）

其中：常数p(正整数)为模型的阶数，常系数a₁,a₂,…,a_p为模型参数，{r_k}为均值为0、方差为σ²的白噪声，p阶模型简记为AR(p)。自回归-滑动平均混合模型是在自回归模型的基础上减去过去q个时刻上白噪声的线性组合，表示为：

y_k=a₁y_k-1+…+a_py_k-p+r_k-θ₁r_k-1-θ₂r_k-2-…-θ_qr_k-q    （5-10）

其中p和q为模型阶数，a₁,a₂,…,a_p、θ₁,θ₂,…,θ_q为ARMA模型的参数，模型简记为ARMA(p，q)。

一般首先对数据序列进行AR(p)建模，如果找不到适用性模型，再进行ARMA(p，q)建模。在估计模型参数之前，需要定义模型的结构，即确定模型的阶数。本发明采用Akaike的AIC准则确定模型阶数。AIC定阶准则的定义：

$\mathrm{AIC}(p,q)=\log \mathrm{\&delta;}+\frac{2(p+q)}{N}---(5-11)$

其中δ是拟合残差的方差，p、q分别是AR和MA部分的阶数，N是参与估计的样本个数。其定阶思想是从低阶到高阶寻优，先设定模型的阶数，利用最小二乘法等算法估计出模型的参数和残差。再利用上式求出每个模型的AIC值，AIC值最小的模型结构是最优的。

2）AR模型的参数估计

本发明采用下面的快速算法-RLS法进行AR模型的参数估计。

考虑式（5-9）所示的AR(p)模型，阶次p已知，参数a_i和σ²未知。问题是基于已知观测值（y_k，y_k-1，…，y₀，…，y_1-p）求估值和

式（5-10）可以写成如下向量形式：

式中T为矩阵的转置，

α^T=(a₁,…,a_p)。

定义

AR(p)模型参数α的估值公式如下：

初值

和P₀是利用少量观测数据（y₁,…,y_i0），通过以下两式来求得：

$P_0=P_{i0}={\left[{\mathrm{\&phi;}}_{i0}^T{\mathrm{\&phi;}}_{i0}\right]}^{-1}---(5-15)$

${\widehat{\mathrm{\&alpha;}}}_0={\widehat{\mathrm{\&alpha;}}}_{i0}=P_{i0}{\mathrm{\&phi;}}_{i0}^T{Y}_{i0}---(5-16)$

其中 $Y_{i0}=\left(\begin{array}{c}{y}_1\\ .\\ .\\ .\\ y_{i0}\end{array}\right).$

3）ARMA模型的参数估计

本发明采用长自回归白噪声估计法建立式（5-10）所示ARMA模型。其主要步骤为：

①建立长自回归模型AR(p_N)。阶数取lg N的适当倍数，即p_N=(lgN)^1+δ，选δ为一正数：0≤δ≤1。AR(p_N)的自回归参数为

由线性最小二乘估计(LS估计)得到

$\widehat{\&PartialD;}={({\widehat{\&PartialD;}}_1,{\widehat{\&PartialD;}}_2,...,{\widehat{\&PartialD;}}_{\mathrm{pN}})}^T={\left(X_{\mathrm{pN}}^T{X}_{\mathrm{pN}}\right)}^{-1}{X}_{\mathrm{pN}}^T{Y}_{\mathrm{pN}}---(5-17)$

其中 $X_{\mathrm{pN}}=\left[\begin{array}{cccc}{x}_{\mathrm{pN}}& x_{\mathrm{pN}-1}& ...& x_1\\ x_{\mathrm{pN}+1}& x_{\mathrm{pN}}& ...& x_2\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ x_{N-1}& x_{N-2}& ...& x_{N-\mathrm{pN}}\end{array}\right],$ Y_pN＝(x_pN+1，x_pN+2，…，x_N)^T。

②求长自回归模型残差，检验其独立性。用步骤①所得的及样本值x=(x₁,x₂,…,x_N)^T计算残差

${\widehat{\&PartialD;}}_t=x_t-\underset{j=1}{\overset{P_N}{\mathrm{\&Sigma;}}}{\widehat{\&PartialD;}}_j{x}_{t-j},t=P_N+1,P_N+2,...,N---(5-18)$

并检验

的独立性。若不独立，则增大p_N，再重新进行①、②两步，否则进行下一步。

③估计ARMA模型参数，联合白噪声估计值

t＝P_N+1,P_N+2,…,N和样本值x_t，t=1,2,…,N，按LS估计ARMA(p，q)模型的诸参数值：

$\widehat{\mathrm{\&beta;}}=\left[\begin{array}{c}\widehat{\mathrm{\&alpha;}}\\ \widehat{\mathrm{\&theta;}}\end{array}\right]={\left(X_{\mathrm{pq}}^T{X}_{\mathrm{pq}}\right)}^{-1}{X}_{\mathrm{pq}}^T{Y}_{\mathrm{pq}}---(5-19)$

其中 $\widehat{\mathrm{\&beta;}}={({\hat{a}}_1,{\hat{a}}_2,...,{\hat{a}}_p,-{\widehat{\mathrm{\&theta;}}}_1,-{\widehat{\mathrm{\&theta;}}}_2,...,-{\widehat{\mathrm{\&theta;}}}_q)}^T,$ Y_pq=(x_pN+1,x_pN+2,…,x_N-1，x_N)^T，

$X_{\mathrm{pq}}=\left[\begin{array}{cccccccc}{x}_{\mathrm{pN}}& x_{\mathrm{pN}-1}& ...& x_{\mathrm{pN}+1-p}& {\widehat{\&PartialD;}}_{\mathrm{pN}}& {\widehat{\&PartialD;}}_{\mathrm{pN}-1}& ...& {\widehat{\&PartialD;}}_{\mathrm{pN}+1-q}\\ x_{\mathrm{pN}+1}& x_{\mathrm{pN}}& ...& x_{\mathrm{pN}+2-p}& {\widehat{\&PartialD;}}_{\mathrm{pN}+1}& {\widehat{\&PartialD;}}_{\mathrm{pN}}& ...& {\widehat{\&PartialD;}}_{\mathrm{pN}+2-q}\\ .& & & & & & & .\\ .& & & & & & & .\\ .& & & & & & & .\\ x_{N-2}& x_{N-3}& ...& x_{N-1-p}& {\widehat{\&PartialD;}}_{N-2}& {\widehat{\&PartialD;}}_{N-3}& ...& {\widehat{\&PartialD;}}_{N-1-q}\\ x_{N-1}& x_{N-2}& ...& x_{N-p}& {\widehat{\&PartialD;}}_{N-1}& {\widehat{\&PartialD;}}_{N-2}& ...& {\widehat{\&PartialD;}}_{N-q}\end{array}\right].$

上式中

已在步骤②由长自回归估计得到，所以本算法不涉及非线性求解问题，只用到线性最小二乘估计，建模过程简单，便于计算机实现。

（3）模型的适用性检验。模型适用性的最基本检验是检验模型残差

是否为白噪声。

模型的适用性检验准则，根据检验形式的不同，一般分为三类：第一，检验残差是否为白噪声，称为白噪声检验准则；第二，检验残差平方和S或残差方差

是否显著减小，称为残差平方和(残差方差)检验准则；第三，Akaike信息准则，即考虑残差方差

的下降和模型阶次的升高带来的利弊。本发明采用残差平方和检验准则中的F-准则进行模型的适用性检验。按前述方法计算出残差

后，可按下式计算残差平方和S：

$S=\underset{t=n+1}{\overset{N}{\mathrm{\&Sigma;}}}{\&PartialD;}_t^2---(5-20)$

若高阶模型的残差平方和S_h与低阶模型的残差平方和S_l相比有显著性下降，则接受模型，认为其是适用的，否则拒绝。

步骤三：采用卡尔曼滤波算法滤除光纤陀螺零偏数据中的随机噪声。

对于式（5-9）所示的AR(p)模型，设系统的状态为X_k=[y_k,y_k-1，…,y_kp+1]^T，过程噪声为V_k=[r_k,0,…,0]_1×p ^T，状态方程为：

X_k=AX_k-1+BV_k    （5-21）

其中， $A={\left(\begin{array}{ccccc}{a}_1& a_2& ...& a_{p-1}& a_p\\ 1& 0& ...& 0& 0\\ 0& 1& ...& 0& 0\\ .& .& & .& .\\ .& .& & .& .\\ .& .& & .& .\\ 0& 0& ...& 1& 0\end{array}\right)}_{p\&times;p},$ $B={\left(\begin{array}{cccc}1& 0& ...& 0\\ 0& 0& ...& 0\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ 0& 0& ...& 0\end{array}\right)}_{p\&times;p}.$

设光纤陀螺静态输出为测量量Z_k，则系统测量方程为：

Z_k=CX_k+W_k    （5-22）

其中，C=[1 0…0]_1×p，W_k为测量误差。V_k、W_k的统计特性为：

均值E(V_k)＝E(W_k)＝0

自相关函数

互相关函数

卡尔曼滤波的递推算式为：

状态一步预测

协方差阵一步预测P_k|k-1＝AP_k-1A^T+BOB^T

滤波增益K_k＝P_k|k-1C^T(CP_k|k-1C^T+R)^-1

协方差阵估计P_k＝(I-K_kC)P_k|k-1

状态估计 ${\hat{X}}_k={\hat{X}}_{k|k-1}+K_k(Z_k-C{\hat{X}}_{k|k-1})$

其中，R为系统量测噪声方差，值为陀螺的零偏稳定性

Q值取为 ${\left(\begin{array}{cccc}{\mathrm{\&sigma;}}_a^2& 0& ...& 0\\ 0& {\mathrm{\&sigma;}}_a^2& ...& 0\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ 0& 0& ...& {\mathrm{\&sigma;}}_a^2\end{array}\right)}_{p\&times;p},$ P的初值选为单位阵，I为单位阵，X的初值选为零阵。

对于式(5-10)所示的ARMA(p,q)模型，需将上述过程中的V_k改为[r_k,r_k-1，…,r_k-q]^T，矩阵B改为 ${\left(\begin{array}{ccccc}1& -{\mathrm{\&theta;}}_1& -{\mathrm{\&theta;}}_2& ...& -{\mathrm{\&theta;}}_q\\ 0& 0& 0& ...& 0\\ 0& 0& 0& ...& 0\\ .& .& .& & .\\ .& .& .& & .\\ .& .& .& & .\\ 0& 0& 0& 0& 0\end{array}\right)}_{p\&times;q},$ 其他变量不变。将光纤陀螺零偏数据作为卡尔曼滤波器的输入，经计算实现光纤陀螺随机噪声的滤除。

步骤四：光纤陀螺温度漂移误差模型结构及参数辨识。本发明采用统计检验的方法进行光纤陀螺温度漂移误差模型定阶，采用最小二乘方法计算模型参数。

本发明采用多项式拟合光纤陀螺零偏和温度的关系，模型如下：

y=a₀+a₁T+a₂T²+a₃T³+...+a_mT^m    (5-23)

上式中，y为光纤陀螺零偏，单位为°/h；T为光纤陀螺温度，单位为℃。

（1）采用统计检验的方法进行模型定阶。

设总的观测值个数为N，观测值y_i(i＝1,...,N)与其算术平均值

的差值平方和称为离差平方和，记作

$S=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{(y_i-\stackrel{\&OverBar;}{y})}^2=l_{\mathrm{yy}}---(5-24)$

根据上式可得：

$S=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{[(y_i-{\hat{y}}_i)+({\hat{y}}_i-\stackrel{\&OverBar;}{y})]}^2=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{({\hat{y}}_i-\stackrel{\&OverBar;}{y})}^2+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{(y_i-{\hat{y}}_i)}^2+2\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}[(y_i-{\hat{y}}_i)({\hat{y}}_i-\stackrel{\&OverBar;}{y})---(5-25)$

可以证明交叉项

因此总的离差平方和可以分解为两部分，即：

$S=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{({\hat{y}}_i-\stackrel{\&OverBar;}{y})}^2+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{(y_i-{\hat{y}}_i)}^2---(5-26)$

上式可简写为：

S=U+Q    （5-27）

其中， $Q=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{(y_i-{\hat{y}}_i)}^2,$ $U=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{({\hat{y}}_i-\stackrel{\&OverBar;}{y})}^2$ 为回归平方和。

引入方差比F和复相关系数R，即：

$F=\frac{U/m}{Q/N-m-1}---(5-28)$

$R^2=\frac{U}{S}=1-\frac{Q}{S}---(5-29)$

如果在拟合曲线中去掉式（5-23）中第m次项，仅用m-1次多项式拟合，这时对应Q′将增大，T^m的偏回归平方和为U-U'=Q'-Q；偏回归平方和越大，说明第m次项越重要；反之，则不重要。本发明通过显著性检验确定。由式（5-29），可得

$F_m=\frac{Q^{\&prime;}-Q/m}{Q/N-m-1}---(5-30)$

F_m是符合自由度为1和Q/N-m-1的F分布。对于给定的α，由统计表查出F分布的理论临界值F_α(1,N-m-1)，若F_m>F_α(1,N-m-1)，则说明m次项重要，须引入拟合曲线中；反之，则不引入。引入m次项后，还需要考查m+1次项是否显著，把m+1视为m，m视为m-1，重复上述步骤，直到F_m不大于F_α(1,N-m-1)为止。

（2）采用最小二乘方法求解模型参数。

式（5-23）的测量方程为：

$\left.\begin{array}{c}{Y}_1=a_0+a_1{T}_1+a_2{T_1}^2+a_3{T_1}^3+...+a_m{T_1}^m\\ Y_2=a_0+a_1{T}_2+a_2{T_2}^2+a_3{T_2}^3+...+a_m{T_2}^m\\ ...\\ Y_N=a_0+a_1{T}_N+a_2{T_N}^2+a_3{T_N}^3+...+a_m{T_N}^m\end{array}\right\}---(5-31)$

相应的估计量如下：

$\left.\begin{array}{c}{y}_1=a_0+a_1{t}_1+a_2{t_1}^2+a_3{t_1}^3+...+a_m{t_1}^m\\ y_2=a_0+a_1{t}_2+a_2{t_2}^2+a_3{t_2}^3+...+a_m{t_2}^m\\ ...\\ y_N=a_0+a_1{t}_N=a_2{t_N}^2+a_3{t_N}^3+...+a_m{t_N}^m\end{array}\right\}---(5-32)$

其误差方程为：

V=L-Y，Y=TA    （5-33）

其中L为实际测量值， $y=\left[\begin{array}{c}{y}_1\\ y_2\\ .\\ .\\ .\\ y_N\end{array}\right],$ $A=\left[\begin{array}{c}{A}_1\\ A_2\\ .\\ .\\ .\\ A_N\end{array}\right],$ $T=\left[\begin{array}{cccc}1& t_1& ...& {t_1}^m\\ 1& t_2& ...& {t_2}^m\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ 1& t_N& ...& {t_N}^m\end{array}\right].$

根据最小二乘法理论可知，V^TV=最小时，可求解出模型系数A，即A=(T^TT)^-1T^Ty。

在实际工程中，利用最小二乘法时会出现信息矩阵的病态问题，即(T^TT)^-1不存在。为了进一步提高模型参数辨识的精度，本发明采用平衡法进行解决。平衡法指：计算令

得到与原方程同解的方程组DY=DTA，然后根据上述最小二乘法进行计算。

（3）温度补偿具体算法如下：

1）模型为m=1阶，并给定置信度(1-α)；

2）利用卡尔曼滤波消噪后数据，通过最小二乘法解算出m阶模型系数；

3）通过统计方法检查模型最高阶次的显著性；

4）若F_m>F_α(1,N-m-1)，则m阶需引入模型。并将m+1阶引入模型，把m+1视为m，m视为m-1，返回2）步骤；

5）若F_m<F_α(1,N-m-1)，则得到所求模型及参数。

本发明的优点在于：

（1）方法适用于所有惯导系统的光纤陀螺温度补偿，具有通用性；

（2）时间序列分析不仅研究过程的确定性变化，而且更着重于研究过程的随机性变化。采用相邻观测数据的差分方法能方便地去除序列中的非周期项或周期性趋势项，便于对陀螺随机误差模型进行简单确切地描述，具有较高的建模精度。

（3）采用卡尔曼滤波算法对光纤陀螺的输出数据进行处理，能够有效地滤除随机噪声。通过卡尔曼滤波的消噪处理，使光纤陀螺零偏和温度的关系辨识变得容易，更有效。

（4）建立温度漂移误差多项式模型进行温度补偿，结构简单，成本低，易于实现，且实时性好，能够满足高精度实时测量要求。

（四）附图说明

图1为时间序列分析建模流程图；

图2为本发明的温度补偿方法流程图。

图2中符号说明如下：

m为光纤陀螺温度漂移误差多项式模型的阶数，a为显著性水平，1-a为置信度，N为观测值个数，F_m是符合自由度为1的F分布，F_a(1,N-m-1)为F分布表中对应显著性水平a的数值。

（五）具体实施方式

下面将结合附图对本发明作进一步的详细说明。

见图2，本发明一种基于时间序列分析消噪的光纤陀螺温度补偿方法，该方法具体实施步骤如下：

步骤一：设计实验方案，对光纤陀螺进行定点温度实验，采集数据。温度实验采用静态测试，其中，X，Y，Z轴光纤陀螺分别指向东、北、天。为了研究温度对光纤陀螺零偏的影响，分别在-30℃、-20℃、-10℃、0℃、10℃、20℃、30℃、40℃、50℃、60℃环境温度下（或根据工作温度确定所选温度点），对光纤陀螺进行高低温测试。在每一个温度点保温两小时后测量一小时，并通过采集软件记录光纤陀螺自身的温度和相应光纤陀螺零偏值。

步骤二：对陀螺零偏数据进行时间序列分析，建立光纤陀螺随机误差的数学模型。见图1，具体的建模步骤为：

（1）对陀螺测试所得到的样本数据序列进行统计检验和预处理。首先剔除数据中的异常值，其次进行陀螺测试数据的平稳性检验，如发现为非平稳的随机时间序列，应提取其中确定性的趋势项，再次进行周期性检验，如发现潜周期分量，应提取其中能量较大的潜周期分量，最后对去除了趋势项和潜周期分量的残差序列进行正态性检验。如果经过检验的陀螺样本数据的残差序列是非正态时间序列，则进行差分处理使之成为正态时间序列，然后建立随机误差模型。具体方法如下。

1）数据异常值剔除

本发明采用拉依达准则判别异常值。假定数据的总体是正态分布，则：

P(|x-μ|>3σ)<0.003    （5-1）

其中：x为随机变量，μ和σ为数据总体的数学期望和标准差，P为满足括号内条件的概率。

设数据为x₁,x₂,…,x_n，均值为

残差为V_k(k=1,2,…,n)，标准差为：

$\mathrm{\&sigma;}=\sqrt{\frac{1}{n-1}\underset{k}{\mathrm{\&Sigma;}}{V}_k^2}=\sqrt{\frac{1}{n-1}[\underset{k}{\mathrm{\&Sigma;}}{\left(x_k\right)}^2-{\left(\underset{k}{\mathrm{\&Sigma;}}\stackrel{\&OverBar;}{x}\right)}^2/n]}---(5-2)$

若某个测量值x_d的残差V_d(1<d<n)满足|V_d|>3σ，则认为x_d是异常值，予以剔除。

2）平稳性检验

时间序列的平稳性是建模的重要前提。设观测序列的值为x_k(k=1,2,…,n)，均值为

用符号“+”表示

“-”表示

在保持数据原有顺序的情况下，游程定义为具有相同符号的序列，这种符号可以把观测值分成两个互相排斥的类。游程检验假设：样本数据出现的顺序没有明显的趋势，就是平稳的。游程太多或太少都被认为存在非平稳性趋势。设：N₁为一种符号出现的总数；N₂为另一种符号出现的总数；r为游程的总数。当N₁或N₂超过15时可以用正态分布来近似，此时选用统计量为：

其中： ${\mathrm{\&mu;}}_r=\frac{2N_1{N}_2}{N}+1;$ ${\mathrm{\&sigma;}}_r=\sqrt{2N_1{N}_2(2N_1{N}_2-N)/N^2(N-1)};$ N=N₁+N₂。

对于显著水平a=0.05的情况，如果|Z|≤1.96(按2σ原则)，可接受原假设，认为序列是平稳的，否则拒绝。

当平稳性检验中判断{x_k}是非平稳时间序列时，应提取{x_k}中所含的非平稳部分，称为趋势项d_k，一般表示为时间k的多项式：

d_k=a₀+a₁k+a₂k²+…+a_mk^m    （5-4）

式中，a₀,a₁,a₂,…,a_m是多项式数学模型的系数，可由回归分析求出。

光纤陀螺输出的非平稳时间序列中，较常见的趋势项为：

d_k=a₀+a₁k    （5-5）

它表明光纤陀螺输出的均值随时间呈线性变化，但一般是极为缓慢的。在光纤陀螺输出非平稳时间序列{x_k}中去除出了趋势项d_k以后，得到时间序列{y_k}：

y_k=x_k-d_k    （5-6）

再对时间序列{y_k}建模。

3）周期性检验

周期性检验用来识别光纤陀螺输出数据中是否包含随机量以外的周期性分量。周期性检验的方法是直接考察由输出数据得到的概率密度函数、自相关函数或功率谱密度函数的图形。

输出数据中如果存在周期性分量，会在概率密度函数图形中反映出来。随机量和周期性分量的概率密度函数图形的不同之处：随机量的概率密度函数图形呈钟形曲线，周期性分量的概率密度函数图形呈盆形曲线；如果输出数据中同时含有随机量和周期性分量，会出现双峰形曲线，峰值处对应的横坐标值为周期性分量的振幅。

输出数据中如果存在周期性分量，也会在自相关函数图形中反映出来。随机量和周期性分量的自相关函数图形的不同之处：随机量的自相关函数图形在时间间隔增大时，是一条衰减的曲线，周期性分量的自相关函数图形不管时间间隔怎样增大，总是一条不衰减的振荡曲线；如果输出数据中同时含有随机量和周期性分量，其自相关函数图形在一定的时间间隔内呈衰减趋势，然后维持不衰减的振荡。

输出数据中如果存在周期性分量，还会在功率谱密度函数图形中反映出来。当输出数据中同时含有随机量和周期性分量时，功率谱密度函数曲线中会有明显的尖峰，尖峰处所对应的横坐标值为周期性分量的频率。

为了寻找隐含周期，识别并提取随机序列中的周期函数项，通常采用周期图分析法。分析出周期性分量后，应予以剔除。剔除方法为：根据前述方法分析出的周期性分量对应的频率，并对输出数据序列进行傅里叶级数展开，找出周期性分量对应的傅里叶系数，由系数构造出周期性分量序列，在原始数据序列中减去周期性分量序列，实现周期性分量的剔除。

4）正态性检验

正态性检验用来判断光纤陀螺输出数据序列是否具有正态分布的特性。对输出序列{x_k}正态性的检验，最基本的是检验{x_k}的三阶矩(偏态系数ξ)与四阶矩(峰态系数v)是否满足正态随机变量的特性。偏态系数ξ与峰态系数v的定义为：

$\mathrm{\&xi;}=E{\left[\frac{x_k-{\mathrm{\&mu;}}_x}{{\mathrm{\&sigma;}}_x}\right]}^3;$ $v=E{\left[\frac{x_k-{\mathrm{\&mu;}}_x}{{\mathrm{\&sigma;}}_x}\right]}^4---(5-7)$

偏态系数ξ反映了随机变量的概率密度函数峰的对称性，峰态系数v反映了随机变量的概率密度函数峰的状态。对于正态随机变量：ξ=0；v=3。

对{x_k}计算ξ和v的估计值：

$\left\{\begin{array}{c}\widehat{\mathrm{\&xi;}}=\frac{1}{n{\widehat{\mathrm{\&sigma;}}}_x^3}\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{(x_k-{\widehat{\mathrm{\&mu;}}}_x)}^3\\ \hat{v}=\frac{1}{n{\widehat{\mathrm{\&sigma;}}}_x^4}\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{(x_k-{\widehat{\mathrm{\&mu;}}}_x)}^4\end{array}\right.---(5-8)$

式中，

和

分别是{x_k}的均值和方差的估计值。当算得

时，认为{x_k}是正态序列。如果序列不是平稳的，那么采用若干次差分都可以使序列变成平稳的。

（2）模型形式的选取及参数估计。

1）选择模型的形式

光纤陀螺的建模采用自回归模型（AR模型）、自回归-滑动平均混合模型（ARMA模型）。自回归模型指任何一个时刻k上的数值y_k可以表示为过去p个时刻上数值的线性组合加上k时刻的白噪声，表示为：

y_k=a₁y_k-1+…+a_py_k-p+r_k    （5-9）

其中：常数p(正整数)为模型的阶数，常系数a₁,a₂,…,a_p为模型参数，{r_k}为均值为0、方差为σ²的白噪声，p阶模型简记为AR(p)。自回归-滑动平均混合模型是在自回归模型的基础上减去过去q个时刻上白噪声的线性组合，表示为：

y_k=a₁y_k-1+…+a_py_k-p+r_k-θ₁r_k-1-θ₂r_k-2-…-θ_qr_k-q    （5-10）

其中p和q为模型阶数，a₁,a₂,…,a_p、θ₁,θ₂,…,θ_q为ARMA模型的参数，模型简记为ARMA(p，q)。

一般首先对数据序列进行AR(p)建模，如果找不到适用性模型，再进行ARMA(p，q)建模。在估计模型参数之前，需要定义模型的结构，即确定模型的阶数。本发明采用Akaike的AIC准则确定模型阶数。AIC定阶准则的定义：

$\mathrm{AIC}(p,q)=\log \mathrm{\&delta;}+\frac{2(p+q)}{N}---(5-11)$

其中δ是拟合残差的方差，p、q分别是AR和MA部分的阶数，N是参与估计的样本个数。其定阶思想是从低阶到高阶寻优，先设定模型的阶数，利用最小二乘法等算法估计出模型的参数和残差。再利用上式求出每个模型的AIC值，AIC值最小的模型结构是最优的。

2）AR模型的参数估计

本发明采用下面的快速算法-RLS法进行AR模型的参数估计。

考虑式（5-9）所示的AR(p)模型，阶次p已知，参数a_i和σ²未知。问题是基于已知观测值（y_k，y_k-1，…，y₀，…，y_1-p）求估值

和

式（5-10）可以写成如下向量形式：

式中T为矩阵的转置，α^T=(a₁，…,a_p)。

定义

AR(p)模型参数α的估值公式如下：

初值

和P₀是利用少量观测数据（y₁，…,y_i0），通过以下两式来求得：

$P_0=P_{i0}={\left[{\mathrm{\&phi;}}_{i0}^T{\mathrm{\&phi;}}_{i0}\right]}^{-1}---(5-15)$

${\widehat{\mathrm{\&alpha;}}}_0={\widehat{\mathrm{\&alpha;}}}_{i0}=P_{i0}{\mathrm{\&phi;}}_{i0}^T{Y}_{i0}---(5-16)$

其中 $Y_{i0}=\left(\begin{array}{c}{y}_1\\ .\\ .\\ .\\ y_{i0}\end{array}\right).$

3）ARMA模型的参数估计

本发明采用长自回归白噪声估计法建立式（5-10）所示ARMA模型。其主要步骤为：

①建立长自回归模型AR(p_N)。阶数取lgN的适当倍数，即p_N=(lgN)^1+δ，选δ为一正数：0≤δ≤1。AR(p_N)的自回归参数为

由线性最小二乘估计(LS估计)得到

$\widehat{\&PartialD;}={({\widehat{\&PartialD;}}_1,{\widehat{\&PartialD;}}_2,...,{\widehat{\&PartialD;}}_{\mathrm{pN}})}^T={\left(X_{\mathrm{pN}}^T{X}_{\mathrm{pN}}\right)}^{-1}{X}_{\mathrm{pN}}^T{Y}_{\mathrm{pN}}---(5-17)$

其中 $X_{\mathrm{pN}}=\left[\begin{array}{cccc}{x}_{\mathrm{pN}}& x_{\mathrm{pN}-1}& ...& x_1\\ x_{\mathrm{pN}+1}& x_{\mathrm{pN}}& ...& x_2\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ x_{N-1}& x_{N-2}& ...& x_{N-\mathrm{pN}}\end{array}\right],$ Y_pN=(x_pN+1,x_pN+2,…,x_N)^T。

②求长自回归模型残差，检验其独立性。用步骤①所得的

及样本值x=(x₁,x₂,…,x_N)^T计算残差

${\widehat{\&PartialD;}}_t=x_t-\underset{j=1}{\overset{P_N}{\mathrm{\&Sigma;}}}{\widehat{\&PartialD;}}_j{x}_{t-j},t=P_N+1,P_N+2,...,N---(5-18)$

并检验

的独立性。若不独立，则增大p_N，再重新进行①、②两步，否则进行下一步。

③估计ARMA模型参数，联合白噪声估计值

t=P_N+1,P_N+2,…,N和样本值x_t，t=1,2,…,N，按LS估计ARMA(p，q)模型的诸参数值：

$\widehat{\mathrm{\&beta;}}=\left[\begin{array}{c}\widehat{\mathrm{\&alpha;}}\\ \widehat{\mathrm{\&theta;}}\end{array}\right]={\left(X_{\mathrm{pq}}^T{X}_{\mathrm{pq}}\right)}^{-1}{X}_{\mathrm{pq}}^T{Y}_{\mathrm{pq}}---(5-19)$

其中 $\widehat{\mathrm{\&beta;}}={({\hat{a}}_1,{\hat{a}}_2,...,{\hat{a}}_p,-{\widehat{\mathrm{\&theta;}}}_1,-{\widehat{\mathrm{\&theta;}}}_2,...,-{\widehat{\mathrm{\&theta;}}}_q)}^T,$ Y_pq=(x_pN+1，x_pN+2,…,x_N-1，x_N)^T，

$X_{\mathrm{pq}}=\left[\begin{array}{cccccccc}{x}_{\mathrm{pN}}& x_{\mathrm{pN}-1}& ...& x_{\mathrm{pN}+1-p}& {\widehat{\&PartialD;}}_{\mathrm{pN}}& {\widehat{\&PartialD;}}_{\mathrm{pN}-1}& ...& {\widehat{\&PartialD;}}_{\mathrm{pN}+1-q}\\ x_{\mathrm{pN}+1}& x_{\mathrm{pN}}& ...& x_{\mathrm{pN}+2-p}& {\widehat{\&PartialD;}}_{\mathrm{pN}+1}& {\widehat{\&PartialD;}}_{\mathrm{pN}}& ...& {\widehat{\&PartialD;}}_{\mathrm{pN}+2-q}\\ .& & & & & & & .\\ .& & & & & & & .\\ .& & & & & & & .\\ x_{N-2}& x_{N-3}& ...& x_{N-1-p}& {\widehat{\&PartialD;}}_{N-2}& {\widehat{\&PartialD;}}_{N-3}& ...& {\widehat{\&PartialD;}}_{N-1-q}\\ x_{N-1}& x_{N-2}& ...& x_{N-p}& {\widehat{\&PartialD;}}_{N-1}& {\widehat{\&PartialD;}}_{N-2}& ...& {\widehat{\&PartialD;}}_{N-q}\end{array}\right].$

上式中

已在步骤②由长自回归估计得到，所以本算法不涉及非线性求解问题，只用到线性最小二乘估计，建模过程简单，便于计算机实现。

（3）模型的适用性检验。

模型适用性的最基本检验是检验模型残差

是否为白噪声。模型的适用性检验准则，根据检验形式的不同，一般分为三类：第一，检验残差

是否为白噪声，称为白噪声检验准则；第二，检验残差平方和S或残差方差

是否显著减小，称为残差平方和(残差方差)检验准则；第三，Akaike信息准则，考虑残差方差

的下降和模型阶次的升高带来的利弊。本发明采用残差平方和检验准则中的F-准则进行模型的适用性检验。按前述方法计算出残差

后，可按下式计算残差平方和S：

$S=\underset{t=n+1}{\overset{N}{\mathrm{\&Sigma;}}}{\&PartialD;}_t^2---(5-20)$

若高阶模型的残差平方和S_h与低阶模型的残差平方和S_l相比有显著性下降，则接受模型，认为其是适用的，否则拒绝。

步骤三：采用卡尔曼滤波算法滤除光纤陀螺零偏数据中的随机噪声。

对于式（5-9）所示的AR(p)模型，设系统的状态为X_k=[y_k,y_k-1,…,y_k-p+1]^T，过程噪声为V_k=[r_k,0,…,0]_1×p ^T，状态方程为：

X_k=AX_k-1+BV_k    （5-21）

其中， $A={\left(\begin{array}{ccccc}{a}_1& a_2& ...& a_{p-1}& a_p\\ 1& 0& ...& 0& 0\\ 0& 1& ...& 0& 0\\ .& .& & .& .\\ .& .& & .& .\\ .& .& & .& .\\ 0& 0& ...& 1& 0\end{array}\right)}_{p\&times;p},$ $B={\left(\begin{array}{cccc}1& 0& ...& 0\\ 0& 0& ...& 0\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ 0& 0& ...& 0\end{array}\right)}_{p\&times;p}.$

设光纤陀螺静态输出为测量量Z_k，则系统测量方程为：

Z_k=CX_k+W_k    （5-22）

其中，C=[1 0…0]_1×p，W_k为测量误差。V_k、W_k的统计特性为：

均值E(V_k)=E(W_k)=0

自相关函数

互相关函数

卡尔曼滤波的递推算式为：

状态一步预测

协方差阵一步预测P_k|k-1=AP_k-1A^T+BQB^T

滤波增益K_k=P_k|k-1C^T(CP_k|k-1C^T+R)^-1

协方差阵估计P_k=(I-K_kC)P_k|k-1

状态估计 ${\hat{X}}_k={\hat{X}}_{k|k-1}+K_k(Z_k-C{\hat{X}}_{k|k-1})$

其中，R为系统量测噪声方差，值为陀螺的零偏稳定性

Q值取为 ${\left(\begin{array}{cccc}{\mathrm{\&sigma;}}_a^2& 0& ...& 0\\ 0& {\mathrm{\&sigma;}}_a^2& ...& 0\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ 0& 0& ...& {\mathrm{\&sigma;}}_a^2\end{array}\right)}_{p\&times;p},$ P的初值选为单位阵，I为单位阵，X的初值选为零阵。

对于式（5-10）所示的ARMA(p,q)模型，需将上述过程中的V_k改为[r_k,r_k-1，…,r_k-q]^T，矩阵B改为 ${\left(\begin{array}{ccccc}1& -{\mathrm{\&theta;}}_1& -{\mathrm{\&theta;}}_2& ...& -{\mathrm{\&theta;}}_q\\ 0& 0& 0& ...& 0\\ 0& 0& 0& ...& 0\\ .& .& .& & .\\ .& .& .& & .\\ .& .& .& & .\\ 0& 0& 0& 0& 0\end{array}\right)}_{p\&times;q},$ 其他变量不变。将光纤陀螺零偏数据作为卡尔曼滤波器的输入，经计算实现光纤陀螺随机噪声的滤除。

步骤四：光纤陀螺温度漂移误差多项式模型结构辨识及参数辨识。本发明采用统计检验的方法进行模型定阶，用最小二乘方法计算模型参数。

本发明采用多项式拟合光纤陀螺零偏和温度的关系，模型如下：

y=a₀+a₁T+a₂T²+a₃T³+...+a_mT^m    （5-23）

上式中，y为光纤陀螺零偏，单位为°/h；T为光纤陀螺温度，单位为℃。

（1）采用统计检验的方法进行模型定阶。

设总的观测值个数为N，观测值y_i(i=1,...,N)与其算术平均值

的差值平方和称为离差平方和，记作

$S=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{(y_i-\stackrel{\&OverBar;}{y})}^2=l_{\mathrm{yy}}---(5-24)$

根据上式可得：

$S=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{[(y_i-{\hat{y}}_i)+({\hat{y}}_i-\stackrel{\&OverBar;}{y})]}^2=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{({\hat{y}}_i-\stackrel{\&OverBar;}{y})}^2+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{(y_i-{\hat{y}}_i)}^2+2\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}[(y_i-{\hat{y}}_i)({\hat{y}}_i-\stackrel{\&OverBar;}{y})---(5-25)$

可以证明交叉项

因此总的离差平方和可以分解为两部分，即：

$S=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{({\hat{y}}_i-\stackrel{\&OverBar;}{y})}^2+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{(y_i-{\hat{y}}_i)}^2---(5-26)$

上式可简写为：

S=U+Q    （5-27）

其中， $Q=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{(y_i-{\hat{y}}_i)}^2,$ $U=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{({\hat{y}}_i-\stackrel{\&OverBar;}{y})}^2$ 为回归平方和。

引入方差比F和复相关系数R，即：

$F=\frac{U/m}{Q/N-m-1}---(5-28)$

$R^2=\frac{U}{S}=1-\frac{Q}{S}---(5-29)$

如果在拟合曲线中去掉式（5-23）中第m次项，仅用m-1次多项式拟合，这时对应Q′将增大，T^m的偏回归平方和为U-U'=Q'-Q；偏回归平方和越大，说明第m次项越重要；反之，则不重要。本发明通过显著性检验确定。由式（5-29），可得

$F_m=\frac{Q^{\&prime;}-Q/m}{Q/N-m-1}---(5-30)$

F_m是符合自由度为1和Q/N-m-1的F分布。对于给定的α，由统计表查出F分布的理论临界值F_α(1,N-m-1)，若F_m>F_α(1,N-m-1)，则说明m次项重要，须引入拟合曲线中；反之，则不引入。引入m次项后，还需要考查m+1次项是否显著，把m+1视为m，m视为m-1，重复上述步骤，直到F_m不大于F_α(1,N-m-1)为止。

（2）采用最小二乘方法求解模型参数。

式（5-23）的测量方程为：

$\left.\begin{array}{c}{Y}_1=a_0+a_1{T}_1+a_2{T_1}^2+a_3{T_1}^3+...+a_m{T_1}^m\\ Y_2=a_0+a_1{T}_2+a_2{T_2}^2+a_3{T_2}^3+...+a_m{T_2}^m\\ ...\\ Y_N=a_0+a_1{T}_N+a_2{T_N}^2+a_3{T_N}^3+...+a_m{T_N}^m\end{array}\right\}---(5-31)$

相应的估计量如下：

$\left.\begin{array}{c}{y}_1=a_0+a_1{t}_1+a_2{t_1}^2+a_3{t_1}^3+...+a_m{t_1}^m\\ y_2=a_0+a_1{t}_2+a_2{t_2}^2+a_3{t_2}^3+...+a_m{t_2}^m\\ ...\\ y_N=a_0+a_1{t}_N=a_2{t_N}^2+a_3{t_N}^3+...+a_m{t_N}^m\end{array}\right\}---(5-32)$

其误差方程为：

V=L-Y，Y=TA    （5-33）

其中L为实际测量值， $y=\left[\begin{array}{c}{y}_1\\ y_2\\ .\\ .\\ .\\ y_N\end{array}\right],$ $A=\left[\begin{array}{c}{A}_1\\ A_2\\ .\\ .\\ .\\ A_N\end{array}\right],$ $T=\left[\begin{array}{cccc}1& t_1& ...& {t_1}^m\\ 1& t_2& ...& {t_2}^m\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ 1& t_N& ...& {t_N}^m\end{array}\right].$

根据最小二乘法理论可知，V^TV=最小时，可求解出模型系数A，即A=(T^TT)^-1T^Ty。

在实际工程中，利用最小二乘法时会出现信息矩阵的病态问题，即(T^TT)^-1不存在。为了进一步提高模型参数辨识的精度，本发明采用平衡法进行解决。平衡法指：计算

令

得到与原方程同解的方程组DY=DTA，

然后根据上述最小二乘法进行计算。

（3）温度补偿具体算法如下：

1）模型为m=1阶，并给定置信度(1-α)；

2）利用卡尔曼滤波消噪后数据，通过最小二乘法解算出m阶模型系数；

3）通过统计方法检查模型最高阶次的显著性；

4）若F_m>F_α(1,N-m-1)，则m阶需引入模型。并将m+1阶引入模型，把m+1视为m，m视为m-1，返回2）步骤；

5）若F_m<F_α(1,N-m-1)，则得到所求模型及参数。

Claims (1)

1.一种基于时间序列分析消噪的光纤陀螺温度补偿方法，其特征在于：该方法具体步骤如下：

步骤一：设计温度特性测试实验方案，对光纤陀螺进行定点温度实验，采集陀螺数据；温度实验采用静态测试，其中，X，Y，Z轴光纤陀螺分别指向东、北、天；为了研究温度对光纤陀螺零偏的影响，分别在-30℃、-20℃、-10℃、0℃、10℃、20℃、30℃、40℃、50℃、60℃环境温度下，对光纤陀螺进行高低温测试；在每一个温度点保温两小时后测量一小时，并通过采集软件记录光纤陀螺自身的温度和相应光纤陀螺零偏值；

步骤二：对陀螺零偏数据进行时间序列分析，建立光纤陀螺随机误差的数学模型；其具体的建模步骤为：

（1）对陀螺测试样本数据进行统计检验和预处理，首先剔除数据中的异常值，其次进行平稳性检验，如发现为非平稳的，应剔除其中确定性的趋势项，再次进行周期性检验，如发现潜周期分量，应剔除其中能量较大的潜周期分量，最后对去除了趋势项和潜周期分量的残差序列进行正态性检验；如果经过检验的陀螺测试数据的残差序列是非正态时间序列，应进行差分处理使之成为正态时间序列，然后建立随机误差模型；具体方法如下：

1）数据异常值剔除

采用拉依达准则判别异常值；假定数据的总体是正态分布，则：

P(|x-μ|>3σ)<0.003    （5-1）

其中：x为随机变量，μ和σ为数据总体的数学期望和标准差，P为满足括号内条件的概率；设数据为x₁,x₂,…,x_n，均值为

残差为V_k(k=1,2,…,n)，标准差为：

$\mathrm{\&sigma;}=\sqrt{\frac{1}{n-1}\underset{k}{\mathrm{\&Sigma;}}{V}_k^2}=\sqrt{\frac{1}{n-1}[\underset{k}{\mathrm{\&Sigma;}}{\left(x_k\right)}^2-{\left(\underset{k}{\mathrm{\&Sigma;}}\stackrel{\&OverBar;}{x}\right)}^2/n]}---(5-2)$

若某个测量值x_d的残差V_d(1<d<n)满足|V_d|>3σ，则认为x_d是异常值，予以剔除；

2）平稳性检验

时间序列的平稳性是建模的重要前提，设观测序列的值为x_k(k=1,2,…,n)，均值为

用符号“+”表示

“-”表示

在保持数据原有顺序的情况下，游程定义为具有相同符号的序列，这种符号可以把观测值分成两个互相排斥的类；游程检验假设：样本数据出现的顺序没有明显的趋势，就是平稳的，游程太多或太少都被认为存在非平稳性趋势；设：N₁为一种符号出现的总数，N₂为另一种符号出现的总数，r为游程的总数，当N₁或N₂超过15时用正态分布来近似，此时选用统计量为：

其中： ${\mathrm{\&mu;}}_r=\frac{2N_1{N}_2}{N}+1;$ ${\mathrm{\&sigma;}}_r=\sqrt{2N_1{N}_2(2N_1{N}_2-N)/N^2(N-1)};$ N=N₁+N₂；

对于显著水平α=0.05的情况，如果|Z|≤1.96(按2σ原则)，接受原假设，认为序列是平稳的，否则拒绝；

当平稳性检验中判断{x_k}是非平稳时间序列时，应提取{x_k}中所含的非平稳部分，称为趋势项d_k，表示为时间k的多项式：

d_k=a₀+a₁k+a₂k²+…+a_mk^m    （5-4）

式中，a₀,a₁,a₂,…,a_m是多项式数学模型的系数，可由回归分析求出；

光纤陀螺输出的非平稳时间序列中，较常见的趋势项为：

d_k=a₀+a₁k    （5-5）

它表明光纤陀螺输出的均值随时间呈线性变化，但是极为缓慢的；在光纤陀螺输出非平稳时间序列{x_k}中去除出了趋势项d_k以后，得到时间序列{y_k}：

y_k=x_k-d_k    （5-6）

再对时间序列{y_k}建模；

3）周期性检验

周期性检验用来识别光纤陀螺输出数据中是否包含随机量以外的周期性分量，周期性检验的方法是直接考察由输出数据得到的概率密度函数、自相关函数或功率谱密度函数的图形；

输出数据中如果存在周期性分量，会在概率密度函数图形中反映出来，随机量和周期性分量的概率密度函数图形的不同之处：随机量的概率密度函数图形呈钟形曲线，周期性分量的概率密度函数图形呈盆形曲线；如果输出数据中同时含有随机量和周期性分量，会出现双峰形曲线，峰值处对应的横坐标值为周期性分量的振幅；

输出数据中如果存在周期性分量，也会在自相关函数图形中反映出来，随机量和周期性分量的自相关函数图形的不同之处：随机量的自相关函数图形在时间间隔增大时，是一条衰减的曲线，周期性分量的自相关函数图形不管时间间隔怎样增大，总是一条不衰减的振荡曲线；如果输出数据中同时含有随机量和周期性分量，其自相关函数图形在一定的时间间隔内呈衰减趋势，然后维持不衰减的振荡；

输出数据中如果存在周期性分量，还会在功率谱密度函数图形中反映出来，当输出数据中同时含有随机量和周期性分量时，功率谱密度函数曲线中会有明显的尖峰，尖峰处所对应的横坐标值为周期性分量的频率；

为了寻找隐含周期，识别并提取随机序列中的周期函数项，通常采用周期图分析法，分析出周期性分量后，应予以剔除；剔除方法为：根据前述方法分析出的周期性分量对应的频率，并对输出数据序列进行傅里叶级数展开，找出周期性分量对应的傅里叶系数，由系数构造出周期性分量序列，在原始数据序列中减去周期性分量序列，实现周期性分量的剔除；

4）正态性检验

正态性检验用来判断光纤陀螺输出数据序列是否具有正态分布的特性，对输出序列{x_k}正态性的检验，最基本的是检验{x_k}的三阶矩即偏态系数ξ与四阶矩即峰态系数v是否满足正态随机变量的特性；偏态系数ξ与峰态系数v的定义为：

$\mathrm{\&xi;}=E{\left[\frac{x_k-{\mathrm{\&mu;}}_x}{{\mathrm{\&sigma;}}_x}\right]}^3;$ $v=E{\left[\frac{x_k-{\mathrm{\&mu;}}_x}{{\mathrm{\&sigma;}}_x}\right]}^4---(5-7)$

偏态系数ξ反映了随机变量的概率密度函数峰的对称性，峰态系数v反映了随机变量的概率密度函数峰的状态；对于正态随机变量：ξ=0；v=3；

对{x_k}计算ξ和v的估计值：

$\left\{\begin{array}{c}\widehat{\mathrm{\&xi;}}=\frac{1}{n{\widehat{\mathrm{\&sigma;}}}_x^3}\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{(x_k-{\widehat{\mathrm{\&mu;}}}_x)}^3\\ \hat{v}=\frac{1}{n{\widehat{\mathrm{\&sigma;}}}_x^4}\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{(x_k-{\widehat{\mathrm{\&mu;}}}_x)}^4\end{array}\right.---(5-8)$

式中，

和

分别是{x_k}的均值和方差的估计值；当算得

时，认为{x_k}是正态序列；如果序列不是平稳的，那么采用若干次差分都使序列变成平稳的；

（2）模型形式的选取及参数估计；其具体实现过程如下：

1）选择模型的形式

光纤陀螺的建模采用自回归模型即AR模型、自回归-滑动平均混合模型即ARMA模型；自回归模型指任何一个时刻k上的数值y_k表示为过去p个时刻上数值的线性组合加上k时刻的白噪声，表示为：

y_k=a₁y_k-1+…+a_py_k-p+r_k    （5-9）

其中：常数p(正整数)为模型的阶数，常系数a₁,a₂,…,a_p为模型参数，{r_k}为均值为0、方差为σ²的白噪声，p阶模型简记为AR(p)；自回归-滑动平均混合模型是在自回归模型的基础上减去过去q个时刻上白噪声的线性组合，表示为：

y_k=a₁y_k-1+…+a_py_k-p+r_k-θ₁r_k-1-θ₂r_k-2-…-θ_qr_k-q    （5-10）

其中，p和q为模型阶数，a₁,a₂,…,a_p、θ₁,θ₂,…,θ_q为ARMA模型的参数，模型简记为ARMA(p，q)；

一般首先对数据序列进行AR(p)建模，如果找不到适用性模型，再进行ARMA(p，q)建模；在估计模型参数之前，需要定义模型的结构，即确定模型的阶数；采用Akaike的AIC准则确定模型阶数，AIC定阶准则的定义：

$\mathrm{AIC}(p,q)=\log \mathrm{\&delta;}+\frac{2(p+q)}{N}---(5-11)$

其中，δ是拟合残差的方差，p、q分别是AR和MA部分的阶数，N是参与估计的样本个数；其定阶思想是从低阶到高阶寻优，先设定模型的阶数，利用最小二乘法等算法估计出模型的参数和残差，再利用上式求出每个模型的AIC值，AIC值最小的模型结构是最优的；

2）AR模型的参数估计

采用下面的快速算法-RLS法进行AR模型的参数估计；

考虑式（5-9）所示的AR(p)模型，阶次p已知，参数a_i和σ²未知，问题是基于已知观测值（y_k，y_k-1，…，y₀，…，y_1-p）求估值和

式（5-10）写成如下向量形式：

式中T为矩阵的转置，

α^T=(a_T,…,a_p)；

定义

AR(p)模型参数α的估值公式如下：

初值和P₀是利用少量观测数据（y₁,…,y_i0），通过以下两式来求得：

$P_0=P_{i0}={\left[{\mathrm{\&phi;}}_{i0}^T{\mathrm{\&phi;}}_{i0}\right]}^{-1}---(5-15)$

${\widehat{\mathrm{\&alpha;}}}_0={\widehat{\mathrm{\&alpha;}}}_{i0}=P_{i0}{\mathrm{\&phi;}}_{i0}^T{Y}_{i0}---(5-16)$

其中 $Y_{i0}=\left(\begin{array}{c}{y}_1\\ .\\ .\\ .\\ y_{i0}\end{array}\right);$

3）ARMA模型的参数估计

采用长自回归白噪声估计法建立式（5-10）所示ARMA模型，其主要步骤为：

①建立长自回归模型AR(p_N)；阶数取lgN的适当倍数，即p_N=(lgN)^1+δ，选δ为一正数：0≤δ≤1；AR(p_N)的自回归参数为

由线性最小二乘估计即LS估计得到

$\widehat{\&PartialD;}={({\widehat{\&PartialD;}}_1,{\widehat{\&PartialD;}}_2,...,{\widehat{\&PartialD;}}_{\mathrm{pN}})}^T={\left(X_{\mathrm{pN}}^T{X}_{\mathrm{pN}}\right)}^{-1}{X}_{\mathrm{pN}}^T{Y}_{\mathrm{pN}}---(5-17)$

其中， $X_{\mathrm{pN}}=\left[\begin{array}{cccc}{x}_{\mathrm{pN}}& x_{\mathrm{pN}-1}& ...& x_1\\ x_{\mathrm{pN}+1}& x_{\mathrm{pN}}& ...& x_2\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ x_{N-1}& x_{N-2}& ...& x_{N-\mathrm{pN}}\end{array}\right],$ Y_pN=(x_pN+1,x_pN+2,…,x_N)^T；

②求长自回归模型残差，检验其独立性；用步骤①所得的

及样本值x=(x₁,x₂,…,x_N)^T计算残差

${\widehat{\&PartialD;}}_t=x_t-\underset{j=1}{\overset{P_N}{\mathrm{\&Sigma;}}}{\widehat{\&PartialD;}}_j{x}_{t-j},t=P_N+1,P_N+2,...,N---(5-18)$

并检验

的独立性；若不独立，则增大p_N，再重新进行①、②两步，否则进行下一步；

③估计ARMA模型参数，联合白噪声估计值

t=P_N+1,P_N+2,…,N和样本值x_t，t=1,2,…,N，按LS估计ARMA(p，q)模型的诸参数值：

$\widehat{\mathrm{\&beta;}}=\left[\begin{array}{c}\widehat{\mathrm{\&alpha;}}\\ \widehat{\mathrm{\&theta;}}\end{array}\right]={\left(X_{\mathrm{pq}}^T{X}_{\mathrm{pq}}\right)}^{-1}{X}_{\mathrm{pq}}^T{Y}_{\mathrm{pq}}---(5-19)$

其中， $\widehat{\mathrm{\&beta;}}={({\hat{a}}_1,{\hat{a}}_2,...,{\hat{a}}_p,-{\widehat{\mathrm{\&theta;}}}_1,-{\widehat{\mathrm{\&theta;}}}_2,...,-{\widehat{\mathrm{\&theta;}}}_q)}^T,$ Y_pq=(x_pN+1,x_pN+2，…,x_N-1，x_N)^T，

$X_{\mathrm{pq}}=\left[\begin{array}{cccccccc}{x}_{\mathrm{pN}}& x_{\mathrm{pN}-1}& ...& x_{\mathrm{pN}+1-p}& {\widehat{\&PartialD;}}_{\mathrm{pN}}& {\widehat{\&PartialD;}}_{\mathrm{pN}-1}& ...& {\widehat{\&PartialD;}}_{\mathrm{pN}+1-q}\\ x_{\mathrm{pN}+1}& x_{\mathrm{pN}}& ...& x_{\mathrm{pN}+2-p}& {\widehat{\&PartialD;}}_{\mathrm{pN}+1}& {\widehat{\&PartialD;}}_{\mathrm{pN}}& ...& {\widehat{\&PartialD;}}_{\mathrm{pN}+2-q}\\ .& & & & & & & .\\ .& & & & & & & .\\ .& & & & & & & .\\ x_{N-2}& x_{N-3}& ...& x_{N-1-p}& {\widehat{\&PartialD;}}_{N-2}& {\widehat{\&PartialD;}}_{N-3}& ...& {\widehat{\&PartialD;}}_{N-1-q}\\ x_{N-1}& x_{N-2}& ...& x_{N-p}& {\widehat{\&PartialD;}}_{N-1}& {\widehat{\&PartialD;}}_{N-2}& ...& {\widehat{\&PartialD;}}_{N-q}\end{array}\right]$

上式中

已在步骤②由长自回归估计得到，所以该算法不涉及非线性求解问题，只用到线性最小二乘估计，建模过程简单，便于计算机实现；

（3）模型的适用性检验；模型适用性的最基本检验是检验模型残差是否为白噪声；

模型的适用性检验准则，根据检验形式的不同，分为三类：第一，检验残差

是否为白噪声，称为白噪声检验准则；第二，检验残差平方和S或残差方差

是否显著减小，称为残差平方和即残差方差检验准则；第三，Akaike信息准则，即考虑残差方差

的下降和模型阶次的升高带来的利弊；采用残差平方和检验准则中的F-准则进行模型的适用性检验，按前述方法计算出残差后，按下式计算残差平方和S：

$S=\underset{t=n+1}{\overset{N}{\mathrm{\&Sigma;}}}{\&PartialD;}_t^2---(5-20)$

若高阶模型的残差平方和S_h与低阶模型的残差平方和S_l相比有显著性下降，则接受模型，认为其是适用的，否则拒绝；

步骤三：采用卡尔曼滤波算法滤除光纤陀螺零偏数据中的随机噪声；

对于式（5-9）所示的AR(p)模型，设系统的状态为X_k=[y_k,y_k-1，…,y_k-p+1]^T，过程噪声为V_k=[r_k,0,…,0]_1×p ^T，状态方程为：

X_k=AX_k-1+BV_k    （5-21）

其中， $A={\left(\begin{array}{ccccc}{a}_1& a_2& ...& a_{p-1}& a_p\\ 1& 0& ...& 0& 0\\ 0& 1& ...& 0& 0\\ .& .& & .& .\\ .& .& & .& .\\ .& .& & .& .\\ 0& 0& ...& 1& 0\end{array}\right)}_{p\&times;p},$ $B={\left(\begin{array}{cccc}1& 0& ...& 0\\ 0& 0& ...& 0\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ 0& 0& ...& 0\end{array}\right)}_{p\&times;p};$

设光纤陀螺静态输出为测量量Z_k，则系统测量方程为：

Z_k=CX_k+W_k    （5-22）

其中，C=[10…0]_1×p，W_k为测量误差；V_k、W_k的统计特性为：

均值E(V_k)=E(W_k)=0

自相关函数

互相关函数

卡尔曼滤波的递推算式为：

状态一步预测

协方差阵一步预测P_k|k-1=AP_k-1A^T+BQB^T

滤波增益K_k=P_k|k-1C^T(CP_k|k-1C^T+R)^-1

协方差阵估计P_k=(I-K_kC)P_k|k-1

状态估计 ${\hat{X}}_k={\hat{X}}_{k|k-1}+K_k(Z_k-C{\hat{X}}_{k|k-1})$

其中，R为系统量测噪声方差，值为陀螺的零偏稳定性

Q值取为 ${\left(\begin{array}{cccc}{\mathrm{\&sigma;}}_a^2& 0& ...& 0\\ 0& {\mathrm{\&sigma;}}_a^2& ...& 0\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ 0& 0& ...& {\mathrm{\&sigma;}}_a^2\end{array}\right)}_{p\&times;p},$ P的初值选为单位阵，I为单位阵，X的初值选为零阵；

对于式（5-10）所示的ARMA(p,q)模型，需将上述过程中的V_k改为[r_k,r_k-1，…,r_k-q]^T，矩阵B改为 ${\left(\begin{array}{ccccc}1& -{\mathrm{\&theta;}}_1& -{\mathrm{\&theta;}}_2& ...& -{\mathrm{\&theta;}}_q\\ 0& 0& 0& ...& 0\\ 0& 0& 0& ...& 0\\ .& .& .& & .\\ .& .& .& & .\\ .& .& .& & .\\ 0& 0& 0& 0& 0\end{array}\right)}_{p\&times;q},$ 其它变量不变；将光纤陀螺零偏数据作为卡尔曼滤波器的输入，经计算实现光纤陀螺随机噪声的滤除；

步骤四：光纤陀螺温度漂移误差模型结构及参数辨识；采用统计检验的方法进行光纤陀螺温度漂移误差模型定阶，采用最小二乘方法计算模型参数；

采用多项式拟合光纤陀螺零偏和温度的关系，模型如下：

y＝a₀+a₁T+a₂T²+a₃T³+...+a_mT^m    (5-23)

上式中，y为光纤陀螺零偏，单位为°/h；T为光纤陀螺温度，单位为℃；

(1)采用统计检验的方法进行模型定阶；

设总的观测值个数为N，观测值y_i(i＝1,...,N)与其算术平均值的差值平方和称为离差平方和，记作

$S=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{(y_i-\stackrel{\&OverBar;}{y})}^2=l_{\mathrm{yy}}---(5-24)$

根据上式得：

$S=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{[(y_i-{\hat{y}}_i)+({\hat{y}}_i-\stackrel{\&OverBar;}{y})]}^2=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{({\hat{y}}_i-\stackrel{\&OverBar;}{y})}^2+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{(y_i-{\hat{y}}_i)}^2+2\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}[(y_i-{\hat{y}}_i)({\hat{y}}_i-\stackrel{\&OverBar;}{y})---(5-25)$

证明交叉项

因此总的离差平方和分解为两部分，即：

$S=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{({\hat{y}}_i-\stackrel{\&OverBar;}{y})}^2+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{(y_i-{\hat{y}}_i)}^2---(5-26)$

上式简写为：

S＝U+Q    (5-27)

其中， $Q=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{(y_i-{\hat{y}}_i)}^2,$ $U=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{({\hat{y}}_i-\stackrel{\&OverBar;}{y})}^2$ 为回归平方和；

引入方差比F和复相关系数R，即：

$F=\frac{U/m}{Q/N-m-1}---(5-28)$

$R^2=\frac{U}{S}=1-\frac{Q}{S}---(5-29)$

如果在拟合曲线中去掉式(5-23)中第m次项，仅用m-1次多项式拟合，这时对应Q′将增大，T^m的偏回归平方和为U-U′＝Q′-Q；偏回归平方和越大，说明第m次项越重要；反之，则不重要；通过显著性检验确定，由式（5-29），得

$F_m=\frac{Q^{\&prime;}-Q/m}{Q/N-m-1}---(5-30)$

F_m是符合自由度为1和Q/N-m-1的F分布；对于给定的α，由统计表查出F分布的理论临界值F_α(1,N-m-1)，若F_m>F_α(1,N-m-1)，则说明m次项重要，须引入拟合曲线中；反之，则不引入；引入m次项后，还需要考查m+1次项是否显著，把m+1视为m，m视为m-1，重复上述步骤，直到F_m不大于F_α(1,N-m-1)为止；

（2）采用最小二乘方法求解模型参数；

式（5-23）的测量方程为：

$\left.\begin{array}{c}{Y}_1=a_0+a_1{T}_1+a_2{T_1}^2+a_3{T_1}^3+...+a_m{T_1}^m\\ Y_2=a_0+a_1{T}_2+a_2{T_2}^2+a_3{T_2}^3+...+a_m{T_2}^m\\ ...\\ Y_N=a_0+a_1{T}_N+a_2{T_N}^2+a_3{T_N}^3+...+a_m{T_N}^m\end{array}\right\}---(5-31)$

相应的估计量如下：

$\left.\begin{array}{c}{y}_1=a_0+a_1{t}_1+a_2{t_1}^2+a_3{t_1}^3+...+a_m{t_1}^m\\ y_2=a_0+a_1{t}_2+a_2{t_2}^2+a_3{t_2}^3+...+a_m{t_2}^m\\ ...\\ y_N=a_0+a_1{t}_N=a_2{t_N}^2+a_3{t_N}^3+...+a_m{t_N}^m\end{array}\right\}---(5-32)$

其误差方程为：

V=L-Y，Y=TA    （5-33）

其中，L为实际测量值， $y=\left[\begin{array}{c}{y}_1\\ y_2\\ .\\ .\\ .\\ y_N\end{array}\right],$ $A=\left[\begin{array}{c}{A}_1\\ A_2\\ .\\ .\\ .\\ A_N\end{array}\right],$ $T=\left[\begin{array}{cccc}1& t_1& ...& {t_1}^m\\ 1& t_2& ...& {t_2}^m\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ 1& t_N& ...& {t_N}^m\end{array}\right];$

根据最小二乘法理论知，V^TV=最小时，求解出模型系数A，即A=(T^TT)^-1T^Ty；

在实际工程中，利用最小二乘法时会出现信息矩阵的病态问题，即(T^TT)^-1不存在；为了进一步提高模型参数辨识的精度，采用平衡法进行解决；平衡法指：计算

令

得到与原方程同解的方程组DY=DTA，

然后根据上述最小二乘法进行计算；

（3）温度补偿具体算法如下：

1）模型为m=1阶，并给定置信度(1-α)；

2）利用卡尔曼滤波消噪后数据，通过最小二乘法解算出m阶模型系数；

3）通过统计方法检查模型最高阶次的显著性；

4）若F_m>F_α(1,N-m-1)，则m阶需引入模型；并将m+1阶引入模型，把m+1视为m，m视为m-1，返回2）步骤；

5）若F_m<F_α(1,N-m-1)，则得到所求模型及参数。

Images