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

摘要

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

主权项

一种基于时间序列分析消噪的光纤陀螺温度补偿方法，其特征在于：该方法具体步骤如下：步骤一：设计温度特性测试实验方案，对光纤陀螺进行定点温度实验，采集陀螺数据；温度实验采用静态测试，其中，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}{\Sigma }{V}_k^2}=\sqrt{\frac{1}{n-1}[{\left({\Sigma }_k\right)}^2-{\left(\underset{k}{\Sigma \bar{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时用正态分布来近似，此时选用统计量为：其中：${\mu }_r=\frac{{2N}_1{N}_2}{N}+1;{\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}的三阶矩即偏态系数ξ与四阶矩即峰态系数ν是否满足正态随机变量的特性；偏态系数ξ与峰态系数ν的定义为：$\xi =E{\left[\frac{x_k-{\mu }_x}{{\sigma }_x}\right]}^3;v=E{\left[\frac{x_k-{\mu }_x}{{\sigma }_x}\right]}^4---(5-7)$偏态系数ξ反映了随机变量的概率密度函数峰的对称性，峰态系数ν反映了随机变量的概率密度函数峰的状态；对于正态随机变量：ξ＝0；ν＝3；对{x_k}计算ξ和ν的估计值：$\left\{\begin{array}{c}\hat{\xi }=\frac{1}{{n\hat{\sigma }}_x^3}\underset{k=1}{\overset{n}{\Sigma }}{(x_k-{\hat{\mu }}_x)}^3\\ \hat{v}=\frac{1}{{n\hat{\sigma }}_x^4}\underset{k=1}{\overset{n}{\Sigma }}{(x_k-{\hat{\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为模型的阶数，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 \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[{\phi }_i^T{\phi }_{i0}\right]}^{-1}---(5-15)$${\hat{\alpha }}_0={\hat{\alpha }}_{i0}=P_{i0}{\phi }_{i0}^T{Y}_{i0}---(5-15)$其中$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估计得到$\hat{\partial }={({\hat{\partial }}_1,{\hat{\partial }}_2,...,{\hat{\partial }}_{\mathrm{pN}})}^T={\left(X_{\mathrm{pN}}^T{X}_{\mathrm{pN}}\right)}^{-1}{X}_{\mathrm{pN}}^T{Y}_{\mathrm{nP}}---(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计算残差${\hat{\partial }}_t=x_t-\underset{j=1}{\overset{P_N}{\Sigma }}{\hat{\partial }}_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)模型的诸参数值：$\hat{\beta }=\left[\begin{array}{c}\hat{\alpha }\\ \hat{\theta }\end{array}\right]={\left(X_{\mathrm{pq}}^T{X}_{\mathrm{pq}}\right)}^{-1}{X}_{\mathrm{pq}}^T{Y}_{\mathrm{pq}}---(5-19)$其中，$\hat{\beta }={({\hat{a}}_1,{\hat{a}}_2,...,{\hat{a}}_p,-{\hat{\theta }}_1,-{\hat{\theta }}_2,...,-{\hat{\theta }}_q)}^T,$Y_pq＝(x_pN+1,x_pN+2,…,x_N‑1,x_N)^T，$X_{\mathrm{pq}}=\begin{array}{}\left[\begin{array}{cccccccc}{x}_{\mathrm{pN}}& x_{\mathrm{pN}-1}& ...& x_{\mathrm{pN}+1-p}& {\hat{\partial }}_{\mathrm{pN}}& {\hat{\partial }}_{\mathrm{pN}-1}& ...& {\hat{\partial }}_{\mathrm{pN}+1-q}\\ x_{\mathrm{pN}+1}& x_{\mathrm{pN}}& ...& x_{\mathrm{pN}+2-p}& {\hat{\partial }}_{\mathrm{pN}+1}& {\hat{\partial }}_{\mathrm{pN}}& ...& {\hat{\partial }}_{\mathrm{pN}+2-q}\\ .& & & & & & & .\\ .& & & & & & & .\\ .& & & & & & & .\\ x_{N-2}& x_{N-3}& ...& x_{N-1-p}& {\hat{\partial }}_{N-2}& {\hat{\partial }}_{N-3}& ...& {\hat{\partial }}_{N-1-q}\\ x_{N-1}& x_{N-2}& ...& x_{N-p}& {\hat{\partial }}_{N-1}& {\hat{\partial }}_{N-2}& ...& {\hat{\partial }}_{N-q}\end{array}\right]\end{array}$上式中已在步骤②由长自回归估计得到，所以该算法不涉及非线性求解问题，只用到线性最小二乘估计，建模过程简单，便于计算机实现；(3)模型的适用性检验；模型适用性的最基本检验是检验模型残差是否为白噪声；模型的适用性检验准则，根据检验形式的不同，分为三类：第一，检验残差是否为白噪声，称为白噪声检验准则；第二，检验残差平方和S或残差方差是否显著减小，称为残差平方和即残差方差检验准则；第三，Akaike信息准则，即考虑残差方差的下降和模型阶次的升高带来的利弊；采用残差平方和检验准则中的F‑准则进行模型的适用性检验，按前述方法计算出残差后，按下式计算残差平方和S：$S=\underset{t=n+1}{\overset{N}{\Sigma }}{\partial }_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自相关函数互相关函数卡尔曼滤波的递推算式为：状态一步预测${\hat{X}}_{k|k-1}=A{\hat{X}}_{k-1}$协方差阵一步预测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}{\sigma }_a^2& 0& ...& 0\\ 0& {\sigma }_a^2& ...& 0\\ .& & & .\\ .& & & .\\ .& & & .\\ 0& 0& ...& {\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}0& -{\theta }_1& -{\theta }_2& ...& -{\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)与其算术平均值y的差值平方和称为离差平方和，记作$S=\underset{i=1}{\overset{N}{\Sigma }}{(y_i-\bar{y})}^2=l_{\mathrm{yy}}---(5-24)$根据上式得：$S=\underset{i=1}{\overset{N}{\Sigma }}{[(y_i-{\hat{y}}_i)+({\hat{y}}_i-\bar{y})]}^2=\underset{i=1}{\overset{N}{\Sigma }}{({\hat{y}}_i-\bar{y})}^2+\underset{i=1}{\overset{N}{\Sigma }}{(y_i-{\hat{y}}_i)}^2+2\underset{i=1}{\overset{N}{\Sigma }}[(y_i-{\hat{y}}_i)({\hat{y}}_i-\bar{y})---(5-25)$证明交叉项因此总的离差平方和分解为两部分，即：$S=\underset{i=1}{\overset{N}{\Sigma }}{({\hat{y}}_i-\bar{y})}^2+\underset{i=1}{\overset{N}{\Sigma }}{(y_i-{\hat{y}}_i)}^2---(5-26)$上式简写为：S＝U+Q   (5‑27)其中，$Q=\underset{i=1}{\overset{N}{\Sigma }}{(y_i-{\hat{y}}_i)}^2,U=\underset{i=1}{\overset{N}{\Sigma }}{({\hat{y}}_i-\bar{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不存在；为了进一步提高模型参数辨识的精度，采用平衡法进行解决；平衡法指：计算$s_i=\underset{1\le j\le n}{\max }\left|T_{\mathrm{ij}}\right|(i=\mathrm{1,2},...,n),$令$D=\mathrm{diag}(\frac{1}{s_1},\frac{1}{s_2},...,\frac{1}{s_n}),$得到与原方程同解的方程组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)，则得到所求模型及参数。