CN102679983B - 一种余度惯性测量单元测量噪声估计和数据融合方法 - Google Patents

Abstract

本发明提出的一种余度惯性测量单元测量噪声估计和数据融合方法，属于惯性导航领域，具体包括获得载体轴的空间合成矢量、将空间合成矢量向余度传感器的各个测量轴投影并得到空间投影数据、得到空间投影数据与惯性测量数据间的投影残差、对余度传感器各个轴向的惯性测量数据加权、利用加权最小二乘法合成载体的空间合成矢量等步骤。本发明通过对余度系统的分析，实现了对余度系统的测量误差的估计，并充分利用测量噪声的估计结果，对测量数据进行了合理的加权，有效的提高了合成矢量的计算精度，并提高了余度惯性导航系统的精度；且本发明流程简单，易于操作。

Description

一种余度惯性测量单元测量噪声估计和数据融合方法

技术领域

本发明涉及惯性导航领域，具体来说，是一种基于加权最小二乘的余度惯性测量单元测量噪声估计方法和数据融合方法。

背景技术

余度惯性测量单元是一种为了增加系统的可靠性和准确度，而采用多个传感器器件，并将其配置成余度方式的技术。余度惯性测量单元相对普通惯性测量单元使用更多的传感器，以减少测量噪声并保证其中部分传感器出现故障时系统仍然能够正常工作。余度惯性测量单元的测量数据需要进行数据融合后才能使用，其目的是将多个测量轴的惯性测量信息融合为三个正交轴的信息。一般情况下，测量单元的传感器具有不同精度，由此导致传感器测量数据的误差也不相同。目前，通常使用基于常规最小二乘的数据融合方法来处理这一问题。由于常规最小二乘法是一种等权重数据加权融合方法，在处理从不同测量精度传感器获得的测量数据时，常规最小二乘数据融合方法不能实现准确的估计。常规最小二乘数据融合方法难以有效利用测量数据获得准确的结果，因此造成了数据融合精度低，对不同噪声分布的适应能力差等现状。

发明内容

针对现有技术中存在的问题，本发明提出了一种基于加权最小二乘的传感器测量误差估计方法和数据融合方法，能够有效的估计传感器测量噪声，提高数据融合精度，提高余度惯性导航系统的精度，本发明不仅仅适用于余度导航系统，也适合于其他余度系统如余度地震波检测系统、余度风速检测系统、余度温度测量系统等，具有实用性和广泛的应用前景。

本发明提出一种余度惯性测量单元测量噪声估计和数据融合方法，具体包括以下几个步骤：

步骤一：获得载体轴的空间合成矢量ω；

通过由载体轴与余度传感器测量轴之间的空间关系，根据最小二乘法对惯性测量数据进行数据融合，得到沿载体轴的空间合成矢量ω；

余度传感器采用六单轴余度传感器对称配置，6个余度传感器的测量轴分别沿正十二面体的六个面法向配置，余度传感器惯性测量单元沿载体轴的空间合成矢量ω为：ω＝[ω_xω_yω_z]^T，其中ω_x、ω_y、ω_z分别表示空间合成矢量ω沿载体的x轴、y轴、z轴的分量；

根据空间几何关系得到六个余度传感器所测的惯性测量数据m为m＝Hω+Δm；

其中，m＝[m₁m₂m₃m₄m₅m₆]^T，m₁～m₆分别表示各个测量轴的惯性测量数据；

Δm＝[Δm₁Δm₂Δm₃Δm₄Δm₅Δm₆]^T，Δm₁～Δm₆分别表示各个测量轴惯性测量数据的测量噪声； $H={\left[\begin{array}{cccccc}\sin \mathrm{\α}& -\sin \mathrm{\α}& \cos \mathrm{\α}& \cos \mathrm{\α}& 0& 0\\ 0& 0& \sin \mathrm{\α}& -\sin \mathrm{\α}& \cos \mathrm{\α}& \cos \mathrm{\α}\\ \cos \mathrm{\α}& \cos \mathrm{\α}& 0& 0& \sin \mathrm{\α}& -\sin \mathrm{\α}\end{array}\right]}^T,$ 表示配置矩阵；通过最小二乘法对余度传感器的惯性测量数据m进行数据融合，得到ω＝(H^TH)^-1H^Tm；

步骤二：根据余度传感器载体轴与测量轴的空间关系，将沿载体轴的空间合成矢量向余度传感器的各个测量轴分别进行投影，得到空间投影数据ω_m为ω_m＝Hω，H为配置矩阵；

步骤三：根据步骤二中得到的空间投影数据ω_m与步骤一中的惯性测量数据m，得到空间投影数据与惯性测量数据间的投影残差δm，选择统计区间，对投影残差进行方差统计，得到余度传感器各轴上的投影残差的方差统计结果

1)获得空间投影数据ω_m与惯性测量数据m间的投影残差δm为投影残差：

δm＝ω_m-m

2)误差方差统计

采用窗口为n的连续滑动窗口法统计随机变量即投影残差δm的方差，对六个轴的投影残差δm_i分别进行方差统计，具体为：

$\begin{array}{cc}E\left({\mathrm{\δm}}_{\mathrm{ji}}^2\right)\stackrel{\mathrm{\Δ}}{=}\frac{1}{n}\underset{k=0}{\overset{n-1}{\mathrm{\Σ}}}{\left({\mathrm{\δm}}_{j(i-k)}\right)}^2& i\≥n,j=\mathrm{1,2},...,6;\end{array}$

其中，表示第j个测量轴、第i次测量的投影残差δm投影残差统计方差值；δm_ji表示第j个测量轴、第i次测量的投影残差，k表示累加的上下界；

令6个测量轴第i次测量的投影残差统计方差向量为：

$E\left({\mathrm{\δm}}_i^2\right)={\left[\begin{array}{cccccc}E\left({\mathrm{\δm}}_{1i}^2\right)& E\left({\mathrm{\δm}}_{2i}^2\right)& E\left({\mathrm{\δm}}_{3i}^2\right)& E\left({\mathrm{\δm}}_{4i}^2\right)& E\left({\mathrm{\δm}}_{5i}^2\right)& E\left({\mathrm{\δ}}_{6i}^2\right)\end{array}\right]}^T,$

6个测量轴第i次测量的真实测量误差方差向量为：

$E\left({\mathrm{\Δm}}_i^2\right)={\left[\begin{array}{cccccc}E\left({\mathrm{\Δm}}_{1i}^2\right)& E\left({\mathrm{\Δm}}_{2i}^2\right)& E\left({\mathrm{\Δm}}_{3i}^2\right)& E\left({\mathrm{\Δm}}_{4i}^2\right)& E\left({\mathrm{\Δm}}_{5i}^2\right)& E\left({\mathrm{\Δ}}_{6i}^2\right)\end{array}\right]}^T,$

根据测量轴与惯性空间的配置关系得到：

E(δm²)＝TE(Δm²)

其中，T表示真实测量误差与统计方差之间的传递关系，获得真实误差方差估计E(Δm²)：

E(Δm²)＝T^-1E(δm²)

实现对余度传感器各测量轴的真实噪声方差的估计；

步骤四：根据步骤三中得到的余度传感器各轴的真实测量误差方差的估计结果，对余度传感器各个轴向的惯性测量数据加权，求出加权阵，再利用加权最小二乘法合成载体的空间合成矢量ω_ji；

采用加权最小二乘法对余度传感器各轴测量误差的真实方差估计结果进行数据融合，得到空间合成矢量ω_ji为ω_ji＝(H^TW_iH)^-1H^TW_im_i；

其中，ω_ji表示j个载体轴、第i次测量的加权最小二乘融合的空间合成矢量；H表示配置矩阵；m_i表示第i次的测量值，

$W=\mathrm{diag}\{\frac{1}{E\left({\mathrm{\Δm}}_{1i}^2\right)},\frac{1}{E\left({\mathrm{\Δm}}_{2i}^2\right)},\frac{1}{E\left({\mathrm{\Δm}}_{3i}^2\right)},\frac{1}{E\left({\mathrm{\Δm}}_{4i}^2\right)},\frac{1}{E\left({\mathrm{\Δm}}_{5i}^2\right)},\frac{1}{E\left({\mathrm{\Δm}}_{6i}^2\right)}\},$ 表示加权阵。

本发明的优点在于：

1、本发明提出的一种余度惯性测量单元测量噪声估计和数据融合方法，通过对余度系统的分析，实现了对余度系统的测量误差的估计；

2、本发明提出的一种余度惯性测量单元测量噪声估计和数据融合方法，充分利用测量噪声的估计结果，对测量数据进行了合理的加权，有效的提高了合成矢量的计算精度，并提高了余度惯性导航系统的精度；

3、本发明提出的一种余度惯性测量单元测量噪声估计和数据融合方法，流程简单，易于操作。

附图说明

图1：本发明提出的一种余度惯性测量单元测量噪声估计和数据融合方法流程图；

图2：本发明中对称六单轴传感器余度配置示意图；

图3：本发明中误差估计方法矢量分析示意图；

图4-A：分别使用最小二乘与加权最小二乘数据融合方法融合的数据进行捷联解算时东向速度解算结果；

图4-B：分别使用最小二乘与加权最小二乘数据融合方法融合的数据进行捷联解算时北向速度解算结果。

具体实施方式

下面结合附图与实施例对本发明做进一步的详细说明。

本发明提出的一种余度惯性测量单元测量噪声估计和数据融合方法，具体包括以下几个步骤：

步骤一：获得载体轴的空间合成矢量ω；

余度传感器中惯性测量单元测量得到沿余度传感器测量轴相对于惯性空间的角速率或加速度分量信息，这些信息称为惯性测量数据。通过由载体轴与余度传感器测量轴之间的空间关系，根据最小二乘法对惯性测量数据进行数据融合，得到沿载体轴的空间合成矢量ω；

所述余度传感器有两种配置方式：分别为测量轴正交配置和斜置配置。在两种配置具有相同容错能力时，正交配置所需要的传感器数要多于斜置配置所需的数目。因此，本发明中余度传感器采用斜置配置方式，具体为六单轴余度传感器对称配置，是指6个余度传感器的测量轴分别沿正十二面体的六个面法向配置，如图2所示，其中，O、X、Y、Z表示参考正交坐标系。这种结构使每个余度传感器的测量轴都彼此相距一个球面角2α＝63°26′5.8″，每一对余度传感器的测量轴都位于参考正交坐标系的一个平面内，并且与正交轴之间的夹角为α。

设余度传感器中陀螺仪或加速度计的随机常值漂移得到补偿，剩余测量噪声为高斯分布、零均值的白噪声，不同余度传感器的测量噪声不相关，且噪声方差在短时间内不发生剧烈的变化和波动，余度传感器惯性测量单元沿载体轴的空间合成矢量ω为：ω＝[ω_xω_yω_z]^T，其中ω_x、ω_y、ω_z分别表示空间合成矢量ω沿载体的x轴、y轴、z轴的分量。当使用陀螺仪作为测量单元时，ω表示空间合成角速度矢量；当使用加速度计作为测量单元时，ω表示空间合成加速度矢量。

根据空间几何关系得到六个余度传感器所测的惯性测量数据m为m＝Hω+Δm，

其中，m＝[m₁m₂m₃m₄m₅m₆]^T，m₁～m₆分别表示各个测量轴的惯性测量数据；

Δm＝[Δm₁Δm₂Δm₃Δm₄Δm₅Δm₆]^T，Δm₁～Δm₆分别表示各个测量轴惯性测量数据的测量噪声；

H为配置矩阵，由余度传感器安装位置与参考坐标系之间的关系得到，由图2所示的余度传感器安装位置与参考坐标系之间的关系得到：

$H={\left[\begin{array}{cccccc}\sin \mathrm{\α}& -\sin \mathrm{\α}& \cos \mathrm{\α}& \cos \mathrm{\α}& 0& 0\\ 0& 0& \sin \mathrm{\α}& -\sin \mathrm{\α}& \cos \mathrm{\α}& \cos \mathrm{\α}\\ \cos \mathrm{\α}& \cos \mathrm{\α}& 0& 0& \sin \mathrm{\α}& -\sin \mathrm{\α}\end{array}\right]}^T;$

通过最小二乘法对余度传感器的惯性测量数据m进行数据融合，得到ω＝(H^TH)^-1H^Tm

步骤二：根据余度传感器载体轴与测量轴的空间关系，将沿载体轴的空间合成矢量向余度传感器的各个测量轴分别进行投影，得到空间投影数据ω_m为ω_m＝Hω；其中，H为配置矩阵，由余度传感器安装位置与参考坐标系之间的关系得到；

步骤三：根据步骤二中得到的空间投影数据ω_m与步骤一中的惯性测量数据m，得到空间投影数据与惯性测量数据间的投影残差δm，选择统计区间，对投影残差进行方差统计，得到余度传感器各轴上的投影残差的方差统计结果

1)获得空间投影数据ω_m与惯性测量数据m间的投影残差δm为投影残差：

δm＝ω_m-m

2)误差方差统计

采用窗口为n的连续滑动窗口法统计随机变量即投影残差δm的方差，对六个轴的投影残差δm_i分别进行方差统计，具体如下：

$\begin{array}{cc}E\left({\mathrm{\δm}}_{\mathrm{ji}}^2\right)\stackrel{\mathrm{\Δ}}{=}\frac{1}{n}\underset{k=0}{\overset{n-1}{\mathrm{\Σ}}}{\left({\mathrm{\δm}}_{j(i-k)}\right)}^2& i\≥n,j=\mathrm{1,2},...,6\end{array}$

其中，表示第j个测量轴、第i次测量的投影残差δm投影残差统计方差值；δm_ji表示第j个测量轴、第i次测量的投影残差。k表示累加的上下界。

令6个测量轴第i次测量的投影残差统计方差向量为：

$E\left({\mathrm{\δm}}_i^2\right)={\left[\begin{array}{cccccc}E\left({\mathrm{\δm}}_{1i}^2\right)& E\left({\mathrm{\δm}}_{2i}^2\right)& E\left({\mathrm{\δm}}_{3i}^2\right)& E\left({\mathrm{\δm}}_{4i}^2\right)& E\left({\mathrm{\δm}}_{5i}^2\right)& E\left({\mathrm{\δ}}_{6i}^2\right)\end{array}\right]}^T,$

6个测量轴第i次测量的真实测量误差方差向量为：

$E\left({\mathrm{\Δm}}_i^2\right)={\left[\begin{array}{cccccc}E\left({\mathrm{\Δm}}_{1i}^2\right)& E\left({\mathrm{\Δm}}_{2i}^2\right)& E\left({\mathrm{\Δm}}_{3i}^2\right)& E\left({\mathrm{\Δm}}_{4i}^2\right)& E\left({\mathrm{\Δm}}_{5i}^2\right)& E\left({\mathrm{\Δ}}_{6i}^2\right)\end{array}\right]}^T,$

δm＝ω_m-m＝H(H^TH)^-1H^Tm-m＝H(H^TH)^-1H^T(m₀+Δm)-(m₀+Δm)，

可以得到：δm＝H(H^TH)^-1H^TΔm-Δm＝[H(H^TH)^-1H^T-I]Δm，其中，I表示六阶单位阵。将矩阵H的值代入即可以得到投影残差δm和测量噪声Δm的具体转移关系：

$\mathrm{\δm}=\left[\begin{array}{cccccc}-0.5000& 0.2236& 0.2236& 0.2236& 0.2236& -0.2236\\ 0.2236& -0.5000& -0.2236& -0.2236& 0.2236& -0.2236\\ 0.2236& -0.2236& -0.5000& 0.2236& 0.2236& 0.2236\\ 0.2236& -0.2236& 0.2236& -0.5000& -0.2236& -0.2236\\ 0.2236& 0.2236& 0.2236& -0.2236& -0.5000& 0.2236\\ -0.2236& -0.2236& 0.2236& -0.2236& 0.2236& -0.5000\end{array}\right]\mathrm{\Δm}$

上式实质为六个等式，将每一个等式两边平方后再取其数学期望，利用不同轴向的测量噪声不相关的假设条件可得到：

$\left[\begin{array}{c}E\left({\mathrm{\δm}}_1^2\right)\\ E\left({\mathrm{\δm}}_2^2\right)\\ E\left({\mathrm{\δm}}_3^2\right)\\ E\left({\mathrm{\δm}}_4^2\right)\\ E\left({\mathrm{\δm}}_5^2\right)\\ E\left({\mathrm{\δm}}_6^2\right)\end{array}\right]=\left[\begin{array}{cccccc}0.25& 0.05& 0.05& 0.05& 0.05& 0.05\\ 0.05& 0.25& 0.05& 0.05& 0.05& 0.05\\ 0.05& 0.05& 0.25& 0.05& 0.05& 0.05\\ 0.05& 0.05& 0.05& 0.25& 0.05& 0.05\\ 0.05& 0.05& 0.05& 0.05& 0.25& 0.05\\ 0.05& 0.05& 0.05& 0.05& 0.05& 0.25\end{array}\right]\left[\begin{array}{c}E\left({\mathrm{\Δm}}_1^2\right)\\ E\left({\mathrm{\Δm}}_2^2\right)\\ E\left({\mathrm{\Δm}}_3^2\right)\\ E\left(\mathrm{\Δ}{m}_4^2\right)\\ E\left({\mathrm{\Δm}}_5^2\right)\\ E\left({\mathrm{\Δm}}_6^2\right)\end{array}\right]\stackrel{\mathrm{\Δ}}{=}T\left[\begin{array}{c}E\left({\mathrm{\Δm}}_1^2\right)\\ E\left({\mathrm{\Δm}}_2^2\right)\\ E\left({\mathrm{\Δm}}_3^2\right)\\ E\left({\mathrm{\Δm}}_4^2\right)\\ E\left({\mathrm{\Δm}}_5^2\right)\\ E\left({\mathrm{\Δm}}_6^2\right)\end{array}\right]$

根据测量轴与惯性空间的配置关系得到：

E(δm²)＝TE(Δm²)

其中， $T=\left[\begin{array}{cccccc}0.25& 0.05& 0.05& 0.05& 0.05& 0.05\\ 0.05& 0.25& 0.05& 0.05& 0.05& 0.05\\ 0.05& 0.05& 0.25& 0.05& 0.05& 0.05\\ 0.05& 0.05& 0.05& 0.25& 0.05& 0.05\\ 0.05& 0.05& 0.05& 0.05& 0.25& 0.05\\ 0.05& 0.05& 0.05& 0.05& 0.05& 0.25\end{array}\right],$ 表示真实测量误差与统计方差之间的传递关系，由配置矩阵H决定。通过如下公式解出真实误差方差估计E(Δm²)：

E(Δm²)＝T^-1E(δm²)

由此实现对余度传感器各测量轴的真实噪声方差的估计。

步骤四：根据步骤三中得到的余度传感器各轴的真实测量误差方差的估计结果，对余度传感器各个轴向的惯性测量数据加权，求出加权阵，再利用加权最小二乘法合成载体的空间合成矢量ω_ji；

当各个余度传感器测量精度不同时，必须区别对待不同精度的测量值，才能实现准确估计，简单的区别对待方法就是加权处理。本发明采用加权最小二乘法对余度传感器各轴测量误差的真实方差估计结果进行数据融合，得到：

ω_ji＝(H^TW_iH)^-1H^TW_im_i

其中，ω_ji表示j个载体轴、第i次测量的加权最小二乘融合的空间合成矢量；H表示配置矩阵；m_i表示第i次的测量值，W表示加权阵，理论上最优的加权阵W为余度传感器各轴测量误差的真实方差E(Δm²)的倒数，即加权阵W表达式为：

$W=\mathrm{diag}\{\frac{1}{E\left({\mathrm{\Δm}}_{1i}^2\right)},\frac{1}{E\left({\mathrm{\Δm}}_{2i}^2\right)},\frac{1}{E\left({\mathrm{\Δm}}_{3i}^2\right)},\frac{1}{E\left({\mathrm{\Δm}}_{4i}^2\right)},\frac{1}{E\left({\mathrm{\Δm}}_{5i}^2\right)},\frac{1}{E\left({\mathrm{\Δm}}_{6i}^2\right)}\}.$

实施例：

本实施例中载体的初始位置为东经116。7°，北纬39°，高度800m，速度257m/s，航向角35.2°，飞行时间200s。陀螺随机漂移0.01°/h，加速度计随机漂移1μg，数据测量周期为0.01s，惯导解算周期为0.01s。应用本发明提供的一种基于最小二乘的余度传感器测量误差估计方法和数据融合方法，如图1所示流程图，通过下列步骤实现：

步骤一：获得载体轴的空间合成矢量；

余度传感器中惯性测量单元测量得到沿余度传感器测量轴相对于惯性空间的角速率或加速度分量信息，这些信息称为惯性测量数据。利用惯性测量数据，通过由载体轴与余度传感器测量轴之间的空间关系对惯性测量数据根据最小二乘法进行数据融合，可得到沿载体轴的空间合成矢量。本实施例以典型的非正交配置(六单轴传感器对称配置)进行说明。

六单轴传感器对称配置是指6个传感器的测量轴分别沿正十二面体的六个面的法向配置，如图2所示，其中，O、X、Y、Z表示参考正交坐标系。这种结构使每个余度传感器的测量轴都彼此相距一个球面角2α＝63°26′5.8″，每一对余度传感器的测量轴都位于参考正交坐标系的一个平面内，并且与正交轴之间的夹角为α。

由于最小二乘法考虑的是数据融合精度问题，因此设余度传感器中陀螺仪、加速度计的随机常值漂移得到补偿，以零均值白噪声模拟各余度传感器测量噪声中的非零分量。为了获取准确的噪声值，从而对噪声的估计结果进行评估，在准确的加速度计和陀螺惯性测量数据中添加零均值白噪声，添加噪声的方差如表1所示。

表1噪声方差

传感器编号	加速度计噪声	陀螺噪声
			h1	0.1g0	0.1*0.01*pi/180
h2	0.1g0	0.1*0.01*pi/180
			h3	0.01g0	10*0.01*pi/180
h4	0.01g0	0.1*0.01*pi/180
			h5	0.1g0	10*0.01*pi/180
h6	0.01g0	0.1*0.01*pi/180

其中h₁～h₆为余度传感器，g₀为重力加速度常量，pi为圆周率常量；

令余度传感器惯性测量单元的状态合成矢量为：ω＝[ω_xω_yω_z]^T

当使用陀螺仪作为测量单元时，ω表示合成空间的角速度矢量；当使用加速度计作为测量单元时，ω表示合成空间的加速度矢量。

如图2所示的空间几何关系可得6个余度传感器所测的惯性测量数据为：

m＝Hω+Δm

其中：

m＝[m₁m₂m₃m₄m₅m₆]^T

H为配置矩阵，由余度传感器安装位置与参考坐标系之间的关系可以得到，由图1余度传感器安装位置与参考坐标系之间的关系可知：

$H={\left[\begin{array}{cccccc}\sin \mathrm{\α}& -\sin \mathrm{\α}& \cos \mathrm{\α}& \cos \mathrm{\α}& 0& 0\\ 0& 0& \sin \mathrm{\α}& -\sin \mathrm{\α}& \cos \mathrm{\α}& \cos \mathrm{\α}\\ \cos \mathrm{\α}& \cos \mathrm{\α}& 0& 0& \sin \mathrm{\α}& -\sin \mathrm{\α}\end{array}\right]}^T$

通过最小二乘法对测量数据进行数据融合，得到惯性信息沿载体轴的合成矢量为：

ω＝(H^TH)^-1H^Tm

步骤二：根据余度传感器载体轴与测量轴的空间关系，将载体轴空间合成矢量向余度传感器各个测量轴分别进行投影，可以得到空间合成矢量在余度传感器各个测量轴上的投影信息，这些信息称为空间投影数据，计算空间投影数据ω_m为ω_m＝Hω。

步骤三：因为通过空间投影数据ω与惯性测量数据m之间的投影残差可以估计测量轴的测量误差，所以对空间投影数据ω与惯性测量数据m之间的投影残差进行方差统计，就可以估计出余度传感器各轴的测量误差；

为了直观表达，下面通过3载体轴余度传感器系统来说明6载体轴余度传感器系统的测量误差估计原理。

三个载体轴的空间配置关系如图3所示：

其中，m＝[m₁m₂m₃]^T表示余度传感器测得的惯性测量数据；ω表示空间合成矢量；ω_m＝[ω_m1ω_m2ω_m3]^T表示合成矢量在余度传感器测量轴上的投影。则H，m、ω、ω_m之间转换关系如下：

m＝Hω

ω＝(H^TH)^-1H^Tm

ω_m＝Hω

由于最小二乘数据融合方法是一种最小化误差平方和的数据匹配方法，合成矢量在余度传感器各测量轴上的投影ω_m就会与余度传感器测得的惯性测量数据m之间存在一定的误差δm，其中：

δm＝ω_m-m

将δm作为统计量，令δm的误差统计区间为n；在这个区间内分别对各个测量轴的Δm进行统计。其中，n为误差统计区间样本数量，若n取值太小，统计特性不明显，容易造成误统计；若n取值太大，虽然统计特性比较可靠，但不利于实时解算分析。对于n的确定，要根据实际采样频率与系统实时性要求来确定。高采样频率下的区间样本数量n应大于低采样频率下的样本数量。

对误差统计项δm进行方差统计，具体统计公式如下所示：

$\begin{array}{cc}E\left({\mathrm{\δm}}_{\mathrm{ji}}^2\right)\stackrel{\mathrm{\Δ}}{=}\frac{1}{n}\underset{k=0}{\overset{n-1}{\mathrm{\Σ}}}{\left({\mathrm{\δm}}_{j(i-k)}\right)}^2& i\≥n,j=\mathrm{1,2},...,6\end{array}$

其中，E(δm_ji ²)表示第j个测量轴、第i次测量的统计方差值；δm_ji表示第j个测量轴、第i次测量的投影残差。

令6个测量轴第i次测量的统计方差向量为：

$E\left({\mathrm{\δm}}_i^2\right)={\left[\begin{array}{cccccc}E\left({\mathrm{\δm}}_{1i}^2\right)& E\left({\mathrm{\δm}}_{2i}^2\right)& E\left({\mathrm{\δm}}_{3i}^2\right)& E\left({\mathrm{\δm}}_{4i}^2\right)& E\left({\mathrm{\δm}}_{5i}^2\right)& E\left({\mathrm{\δ}}_{6i}^2\right)\end{array}\right]}^T;$

6个测量周第i次测量的真实测量误差方差向量为：

$E\left({\mathrm{\Δm}}_i^2\right)={\left[\begin{array}{cccccc}E\left({\mathrm{\Δm}}_{1i}^2\right)& E\left({\mathrm{\Δm}}_{2i}^2\right)& E\left({\mathrm{\Δm}}_{3i}^2\right)& E\left({\mathrm{\Δm}}_{4i}^2\right)& E\left({\mathrm{\Δm}}_{5i}^2\right)& E\left({\mathrm{\Δ}}_{6i}^2\right)\end{array}\right]}^T,$ 则根据测量轴与惯性空间的配置关系得到：

E(δm²)＝TE(Δm²)；

其中， $T=\left[\begin{array}{cccccc}0.25& 0.05& 0.05& 0.05& 0.05& 0.05\\ 0.05& 0.25& 0.05& 0.05& 0.05& 0.05\\ 0.05& 0.05& 0.25& 0.05& 0.05& 0.05\\ 0.05& 0.05& 0.05& 0.25& 0.05& 0.05\\ 0.05& 0.05& 0.05& 0.05& 0.25& 0.05\\ 0.05& 0.05& 0.05& 0.05& 0.05& 0.25\end{array}\right],$ 表示真实测量误差与统计方差之间的传递关系，由配置矩阵H决定。通过如下公式解出真实误差方差向量D_i：

E(Δm²)＝T^-1E(δm²)

由此实现对余度传感器各测量轴的误差方差的估计。

本实施例取统计区间样本数量n＝200(2秒时长)。

步骤四：根据步骤三中得到的余度传感器各轴测量误差的真实方差估计结果，对余度传感器各个轴向的测量数据加权，求出加权阵，再利用加权最小二乘法合成空间矢量；

当各个余度传感器测量精度不同时，必须区别对待不同精度的测量值，才能实现准确估计，简单的区别对待方法就是加权处理。本发明采用加权最小二乘法对余度传感器各轴测量误差的方差估计结果进行数据融合，得到：

ω_jq＝(H^TWH)^-1H^TWm

其中，ω_ji表示j个载体轴、第i次测量的加权最小二乘融合的空间合成矢量；H表示配置矩阵；m_i表示第i次的测量值，W表示加权阵。

理论上最优的加权阵W为真实误差项方差E(Δm²)的倒数，即加权阵W表达式为：

$W=\mathrm{diag}\{\frac{1}{E\left({\mathrm{\Δm}}_{1i}^2\right)},\frac{1}{E\left({\mathrm{\Δm}}_{2i}^2\right)},\frac{1}{E\left({\mathrm{\Δm}}_{3i}^2\right)},\frac{1}{E\left({\mathrm{\Δm}}_{4i}^2\right)},\frac{1}{E\left({\mathrm{\Δm}}_{5i}^2\right)},\frac{1}{E\left({\mathrm{\Δm}}_{6i}^2\right)}\}$

对加噪的数据，分别采用最小二乘与加权最小二乘方法进行数据融合。加速度计a₁、a₂、a₃数据融合仿真结果的数据统计表如表2所示。

表2数据融合统计表

由表2可知，由于传感器的测量误差方差是一段时间的统计值，与真实测量误差方差之间存在误差，所以加权最小二乘数据融合结果要差于理论最优加权最小二乘数据的融合结果，但总的来说，加权最小二乘数据融合方法要优于最小二乘数据融合方法。

由误差分析及大量实验仿真结果可知，当6余度传感器高、低精度不同配置时，本发明的数据融合算法达到最优。将本发明的加权最小二乘数据融合的空间合成矢量结果送入捷联解算过程中使用，进行正常的捷联惯性导航解算过程。

加噪的数据，分别使用最小二乘与加权最小二乘数据融合方法融合的数据进行捷联解算。东向速度和北向速度解算结果如图4所示，数据统计结果如表3所示：

表3融合数据捷联解算对比表

由图4-A、图4-B、表3可知，采用加权最小二乘融合的数据进行捷联解算的误差，要明显小于常规最小二乘融合的数据进行捷联解算的误差。以上结果说明，本发明的数据融合方法能够有效地估计各测量轴的测量误差、提高余度惯性导航系统的精度。

Claims (1)

1.一种余度惯性测量单元测量噪声估计和数据融合方法，其特征在于：具体包括以下几个步骤：

步骤一：获得载体轴的空间合成矢量ω；

通过由载体轴与余度传感器测量轴之间的空间关系，根据最小二乘法对惯性测量数据进行数据融合，得到沿载体轴的空间合成矢量ω；

余度传感器采用六单轴余度传感器对称配置，6个余度传感器的测量轴分别沿正十二面体的六个面法向配置，余度传感器惯性测量单元沿载体轴的空间合成矢量ω为：ω＝[ω_x ω_y ω_z]^T，其中ω_x、ω_y、ω_z分别表示空间合成矢量ω沿载体的x轴、y轴、z轴的分量；

根据空间几何关系得到六个余度传感器所测的惯性测量数据m为m＝Hω+Δm；

其中，m＝[m₁ m₂ m₃ m₄ m₅ m₆]^T，m₁～m₆分别表示各个测量轴的惯性测量数据；

Δm＝[Δm₁ Δm₂ Δm₃ Δm₄ Δm₅ Δm₆]^T，Δm₁～Δm₆分别表示各个测量轴惯性测量数据的测量噪声； $H={\left[\begin{array}{cccccc}\sin \mathrm{\α}& -\sin \mathrm{\α}& \cos \mathrm{\α}& \cos \mathrm{\α}& 0& 0\\ 0& 0& \sin \mathrm{\α}& -\sin \mathrm{\α}& \cos \mathrm{\α}& \cos \mathrm{\α}\\ \cos \mathrm{\α}& \cos \mathrm{\α}& 0& 0& \sin \mathrm{\α}& -\sin \mathrm{\α}\end{array}\right]}^T,$ 表示配置矩阵；通过最小二乘法对余度传感器的惯性测量数据m进行数据融合，得到ω＝(H^TH)^-1H^Tm；每一对余度传感器的测量轴都位于参考正交坐标系的一个平面内，并且与正交轴之间的夹角为α；

步骤二：根据余度传感器载体轴与测量轴的空间关系，将沿载体轴的空间合成矢量向余度传感器的各个测量轴分别进行投影，得到空间投影数据ω_m为ω_m＝Hω，H为配置矩阵；

步骤三：根据步骤二中得到的空间投影数据ω_m与步骤一中的惯性测量数据m，得到空间投影数据与惯性测量数据间的投影残差δm，选择统计区间，对投影残差进行方差统计，得到余度传感器各轴上的投影残差的方差统计结果

1)获得空间投影数据ω_m与惯性测量数据m间的投影残差δm为投影残差：

δm＝ω_m-m

2)误差方差统计

采用窗口为n的连续滑动窗口法统计随机变量即投影残差δm的方差，对六个轴的投影残差δm_i分别进行方差统计，具体为：

$E\left(\mathrm{\δ}{m}_{\mathrm{ji}}^2\right)\stackrel{\mathrm{\Δ}}{=}\frac{1}{n}\underset{k=0}{\overset{n-1}{\mathrm{\Σ}}}{\left(\mathrm{\δ}{m}_{j(i-k)}\right)}^2,i\≥n,j=\mathrm{1,2},...,6;$

其中，表示第j个测量轴、第i次测量的投影残差δm投影残差统计方差值；δm_ji表示第j个测量轴、第i次测量的投影残差，k表示累加的上下界；

令6个测量轴第i次测量的投影残差统计方差向量为：

$E\left(\mathrm{\δ}{m}_i^2\right)={\left[\begin{array}{cccccc}E\left(\mathrm{\δ}{m}_{1i}^2\right)& E\left(\mathrm{\δ}{m}_{2i}^2\right)& E\left(\mathrm{\δ}{m}_{3i}^2\right)& E\left(\mathrm{\δ}{m}_{4i}^2\right)& E\left(\mathrm{\δ}{m}_{5i}^2\right)& E\left(\mathrm{\δ}{m}_{6i}^2\right)\end{array}\right]}^T,$

6个测量轴第i次测量的真实测量误差方差向量为：

$E\left(\mathrm{\Δ}{m}_i^2\right)={\left[\begin{array}{cccccc}E\left(\mathrm{\Δ}{m}_{1i}^2\right)& E\left(\mathrm{\Δ}{m}_{2i}^2\right)& E\left(\mathrm{\Δ}{m}_{3i}^2\right)& E\left(\mathrm{\Δ}{m}_{4i}^2\right)& E\left(\mathrm{\Δ}{m}_{5i}^2\right)& E\left(\mathrm{\Δ}{m}_{6i}^2\right)\end{array}\right]}^T,$

根据测量轴与惯性空间的配置关系：

ω＝(H^TH)^-1H^Tm

得到：

δm＝[H(H^TH)^-1H^T-1]Δm

将矩阵H的值代入，得到投影残差δm和测量噪声Δm的具体转移关系：

$\mathrm{\δm}=\left[\begin{array}{cccccc}-0.5000& 0.2236& 0.2236& 0.2236& 0.2236& -0.2236\\ 0.2236& -0.5000& -0.2236& -0.2236& 0.2236& -0.2236\\ 0.2236& -0.2236& -0.5000& 0.2236& 0.2236& 0.2236\\ 0.2236& -0.2236& 0.2236& -0.5000& -0.2236& -0.2236\\ 0.2236& 0.2236& 0.2236& -0.2236& -0.5000& 0.2236\\ -0.2236& -0.2236& 0.2236& -0.2236& 0.2236& -0.5000\end{array}\right]\mathrm{\Δm}$

上式实质为六个等式，将每一个等式两边平方后再取其数学期望，利用不同轴向的测量噪声不相关的假设条件得到：

$\left[\begin{array}{c}E\left(\mathrm{\δ}{m}_1^2\right)\\ E\left(\mathrm{\δ}{m}_2^2\right)\\ E\left(\mathrm{\δ}{m}_3^2\right)\\ E\left(\mathrm{\δ}{m}_4^2\right)\\ E\left(\mathrm{\δ}{m}_5^2\right)\\ E\left(\mathrm{\δ}{m}_6^2\right)\end{array}\right]=\left[\begin{array}{cccccc}0.25& 0.05& 0.05& 0.05& 0.05& 0.05\\ 0.05& 0.25& 0.05& 0.05& 0.05& 0.05\\ 0.05& 0.05& 0.25& 0.05& 0.05& 0.05\\ 0.05& 0.05& 0.05& 0.25& 0.05& 0.05\\ 0.05& 0.05& 0.05& 0.05& 0.25& 0.05\\ 0.05& 0.05& 0.05& 0.05& 0.05& 0.25\end{array}\right]\left[\begin{array}{c}E\left(\mathrm{\Δ}{m}_1^2\right)\\ E\left(\mathrm{\Δ}{m}_2^2\right)\\ E\left(\mathrm{\Δ}{m}_3^2\right)\\ E\left(\mathrm{\Δ}{m}_4^2\right)\\ E\left(\mathrm{\Δ}{m}_5^2\right)\\ E\left(\mathrm{\Δ}{m}_6^2\right)\end{array}\right]\stackrel{\mathrm{\Δ}}{=}T\left[\begin{array}{c}E\left(\mathrm{\Δ}{m}_1^2\right)\\ E\left(\mathrm{\Δ}{m}_2^2\right)\\ E\left(\mathrm{\Δ}{m}_3^2\right)\\ E\left(\mathrm{\Δ}{m}_4^2\right)\\ E\left(\mathrm{\Δ}{m}_5^2\right)\\ E\left(\mathrm{\Δ}{m}_6^2\right)\end{array}\right]$

根据测量轴与惯性空间的配置关系得到：

E(δm²)＝TE(Δm²)

其中，T表示真实测量误差与统计方差之间的传递关系，获得真实误差方差估计E(Δm²)：

E(Δm²)＝T^-1E(δm²)

实现对余度传感器各测量轴的真实噪声方差的估计；

所述的真实测量误差与统计方差之间的传递关系T为：

$T=\left[\begin{array}{cccccc}0.25& 0.05& 0.05& 0.05& 0.05& 0.05\\ 0.05& 0.25& 0.05& 0.05& 0.05& 0.05\\ 0.05& 0.05& 0.25& 0.05& 0.05& 0.05\\ 0.05& 0.05& 0.05& 0.25& 0.05& 0.05\\ 0.05& 0.05& 0.05& 0.05& 0.25& 0.05\\ 0.05& 0.05& 0.05& 0.05& 0.05& 0.25\end{array}\right];$

步骤四：根据步骤三中得到的余度传感器各轴的真实测量误差方差的估计结果，对余度传感器各个轴向的惯性测量数据加权，求出加权阵，再利用加权最小二乘法合成载体的空间合成矢量ω_ji；

采用加权最小二乘法对余度传感器各轴测量误差的真实方差估计结果进行数据融合，得到空间合成矢量ω_ji为ω_ji＝(H^TW_iH)^-1H^TW_im_i；

其中，ω_ji表示j个载体轴、第i次测量的加权最小二乘融合的空间合成矢量；H表示配置矩阵；m_i表示第i次的测量值，

$W_i=\mathrm{diag}\{\frac{1}{E\left(\mathrm{\Δ}{m}_{1i}^2\right)},\frac{1}{E\left(\mathrm{\Δ}{m}_{2i}^2\right)},\frac{1}{E\left(\mathrm{\Δ}{m}_{3i}^2\right)},\frac{1}{E\left(\mathrm{\Δ}{m}_{4i}^2\right)},\frac{1}{E\left(\mathrm{\Δ}{m}_{5i}^2\right)},\frac{1}{E\left(\mathrm{\Δ}{m}_{6i}^2\right)}\},$ 表示加权阵。