[发明专利]一种用于Allan方差分析技术的噪声系数计算方法

摘要

本发明公开了一种用于Allan方差分析技术的噪声系数估算方法，解决了现有技术在估算5种典型随机噪声系数过程中存在的矩阵求逆奇异、噪声系数有负值、噪声系数不是最优的问题。其主要包括以下步骤：1)给定Allan方差序列和相应的相关时间序列；2)计算相关时间矩阵；3)设定停机条件并对目标函数线性化；4)设置循环迭代相关参数的初值；5)矩阵变换；6)球体变换，计算无约束最优增量；7)计算超平面非负约束最优增量；8)噪声系数更新；9)计算停机指数；10)判断是否满足停机准则。该方法估计出的随机噪声系数可用于惯性器件性能评定或随机测量误差的精确建模，对惯性基组合导航系统导航精度的提高有重要意义。

主权项

1.一种用于Allan方差分析技术的噪声系数计算方法，其特征是，包括以下步骤：1)给定待分析仪器静态时输出数据的Allan方差序列及对应的相关时间序列：将给定的Allan方差序列记为相关时间序列{τ₁,τ₂,…,τ_n}记为τ；2)计算相关时间矩阵：利用相关时间序列τ，根据Allan方差理论值（下标i表示序列的第i个元素，如无特殊说明，下同）与5种典型随机噪声系数C的关系计算相关时间矩阵A；所述与C的关系为：σi2=3τi2Q2+1τiN2+2ln2πB2+τi3K2+τi22R2---(1)]]>式(1)中，Q²为量化噪声系数，N²为白噪声系数，B²为零偏不稳定性系数，K²为速率随机游走系数，R²为速率斜坡系数，为相关时间为τ_i时的Allan方差理论值，则随机噪声系数C为：C＝[Q²,N²,B²,K²,R²]^T上标T表示向量或矩阵的转置，记则对于相关时间序列{τ₁,τ₂,…,τ_n}，相关时间矩阵A为：A=a(τ1)a(τ2)···a(τn)---(2)]]>记Allan方差理论值序列为σ²，根据式(1)，Allan方差理论值序列σ²可表示为：σ²＝A·C                                (3)3)设定停机条件并对目标函数线性化：根据最优化目标函数设定循环迭代的停机条件，并对目标函数进行线性化获取随机噪声系数增量dC的系数矩阵E和无约束条件下最优增量的中间向量b的表达式；所述的最优化目标函数f(C)为：f(C)=[ln(σ2)-ln(σ~2)]T·[ln(σ2)-ln(σ~2)]---(4)]]>式(4)中，σ²为根据式(3)计算的Allan方差理论值序列，为步骤1)中输入的Allan方差序列；定义停机指数e为：e=|f(C)-f(C_)|f(C_)---(5)]]>式(5)中，f(C)为本次循环得到的噪声系数利用式(4)计算的目标函数值，f(C_)为上次循环得到的噪声系数利用式(4)计算的目标函数值；取循环迭代的停机准则为：e≤e₀                                      (6)式(6)中，e为停机指数，e₀为停机条件，可根据所需的噪声系数的精度确定；将式(3)代入式(4)并进行线性化为：f(C+dC)＝dC^T·E^T·E·dC+2·b^T·dC+f(C)       (7)式(7)中，增量dC的系数矩阵E的元素E(i,j)计算方法为：E(i,j)=A(i,j)σi2---(8)]]>E(i,j)表示系数矩阵E的第i行第j列的元素，A(i,j)表示相关时间矩阵A的第i行第j列的元素，表示序列σ²中的第i个元素；式(7)中，无约束条件下最优增量的中间向量b的计算方法为：b=ET·[ln(σ2)-ln(σ~2)]---(9)]]>4)设置循环迭代相关参数的初值：包括5种典型随机噪声系数C的初值C＝[0,0,0,0,0]^T，计算的Allan方差理论值序列σ²的初值目标函数初值和停机指数的初值，目标函数f(C_)的初值应使得首次循环能够进行，停机指数e初值取为一个大于e₀的正数即可；5)矩阵变换：利用相关矩阵A和计算的Allan方差理论值序列σ²计算矩阵E，将矩阵E进行矩阵变换为矩阵以防止后续矩阵求逆出现奇异，并将增量dC、系数C和向量b进行相应变换为和以维持目标函数的平衡；所述矩阵变换操作按如下方法进行：记矩阵E的第j列向量的模为l_j，即：l_j＝|E(:,j)|其中E(:,j)表示矩阵E的第j列向量，|·|表示对向量的取模运算，则经过矩阵变换操作得到的矩阵的元素为：E^(i,j)=E(i,j)lj---(10)]]>增量dC、系数C的元素与变换后的的元素之间的关系为：dC^(j)=lj·dC(j),C^(j)=lj·C(j)]]>其中，dC(j)和分别表示dC和的第j个元素，C(j)和分别表示C和的第j个元素；根据式(9)和式(10)可知的表达式为：b^=E^T·[ln(σ2)-ln(σ~2)]---(11)]]>则式(7)可变换为：f(C+dC)=dC^T·(E^T·E^)·dC^+2·b^T·dC^+f(C)---(12)]]>根据式(12)可知由目标函数等值的所构成的集合是一个与同维的椭球体；6)球体变换，计算无约束最优增量：根据矩阵计算球体变换矩阵H，将增量进行球体变换为x，计算球体变换后无约束条件下的最优增量r；所述的球体变换操作按如下方法进行：由于是正定矩阵，必然存在可逆矩阵H使得：HT·H=E^T·E^---(13)]]>因此可对进行矩阵分解以得到球体变换矩阵H，将进行球体变换为x：x=H·dC^]]>令：r=(H-1)T·b^---(14)]]>则式(12)可重写为：f(C+dC)＝x^T·x+2·r^T·x+f(C)            (15)根据式(15)可知由目标函数等值的x所构成的集合是一个与x同维的球体，则使得目标函数值最小的x应该距离该球体球心r最近，在没有约束的条件下，x和r的最小距离为0，故r即为无约束最优增量；7)计算超平面非负约束最优增量：将球体变换矩阵H的逆矩阵H^-1的各行视为超平面的法向量，确定由H^-1和决定的超平面所约束的非负区域，在非负区域中计算距离无约束最优增量r最近的x；所述的超平面非负约束最优增量按如下方法进行：由于随机噪声系数C的每一个元素都是非负的，即：C≥0故对于所求的x应保证：H-1·x+C^≥0---(16)]]>式(16)所表示的约束可以这样理解：记矩阵H^-1的第i行向量为n_i，的第i个元素与n_i的模长之比为d_i，即：n_i＝H^-1(i,:)di=C^(i)|ni|]]>则n_i可理解为超平面i的法向量，d_i可理解为原点到超平面i的距离，易知由H^-1和所确定的超平面共有5个；若x在超平面i中，则ni·x+C^(i)=0;]]>若x在超平面i之上，则ni·x+C^(i)>0;]]>若x在超平面i之下，则ni·x+C^(i)<0;]]>故式(16)表示x在由H^-1和所确定的5个超平面之中或之上，也就是说由满足式(16)的x所构成的区域即为超平面约束的非负区域；目标函数值最小等价于x距离r最近，即要求(x-r)^T·(x-r)取最小值，表示为（式中min指最小值）：(x-r)^T·(x-r)＝min          (17)若r不在式(16)所表示的约束区域中，则满足式(16)和式(17)的x必然是r在约束区域边界上的投影，即超平面非负约束最优增量x就是r本身，或者x是r在非负区域边界上的投影；8)噪声系数更新：将获得的有约束最优增量x依次进行球体逆变换和矩阵逆变换后得到系数增量dC，计算新的随机噪声系数C；所述的对x进行球体逆变换，其计算方法为：dC^=H-1·x---(18)]]>对进行矩阵逆变换，其计算方法为：dC(j)=dC^(j)lj---(19)]]>则C+dC的值作为更新后的新的随机噪声系数C，即：C＝C+dC           (20)9)计算停机指数：利用新的随机噪声系数C和相关时间矩阵A根据式(3)计算Allan方差理论值序列σ²，根据式(4)计算目标函数值f(C)，再根据式(5)计算停机指数e；10)判断是否满足停机准则：若式(6)成立则满足停机准则，新的随机噪声系数C即为非负最优的；若式(6)不成立，则不满足停机准则，返回步骤5)继续计算。