CN104820758A - 一种基于奇异值分解的传递对准精度评估可观测性分析方法 - Google Patents

Abstract

一种基于奇异值分解的传递对准精度评估可观测性分析方法。本发明的目的在于，提供一种能够有效分析精度评估系统各状态的可观测度的基于奇异值分解的传递对准精度评估可观测性分析方法。按照以下步骤：选取精度评估系统的第一个时间段，并记为j＝1，计算A_j及H_j，并计算该时间段所对应的可观测性矩阵Q_j，选取此时的SOM，并记为Q_s(j)，由系统的外观测量的大小及精度，计算该时间段的外观测量Y_i，求取该时间段SOM的奇异值根据式计算出对应于每个状态变量X₀的奇异值；求取的奇异值即为状态变量可观测性的判断标准。本发明有效分析精度评估系统各状态的可观测度，且可观测度越好，卡尔曼平滑器的性能越好。

Description

一种基于奇异值分解的传递对准精度评估可观测性分析方法

技术领域

本发明涉及的是一种基于奇异值分解的传递对准精度评估可观测性分析方法。

背景技术

传递对准的精度评估是基于卡尔曼滤波和平滑方法进行的，而卡尔曼滤波和平滑与系统的可观测性有很大关系。若系统状态变量的可观测性好，则卡尔曼状态估计收敛速度快、精度高。通过对精度评估系统进行可观测性分析，可以确定所设计的评估方案与最终平滑结果的关系，因此，可利用可观测性分析方法对平滑器的性能进行表征，为后续设计传递对准精度评估方案提供依据。

进行系统可观测性分析一般分为两个内容：一是确定系统是否完全可观测；二是对不完全可观测系统确定哪些状态变量可观测以及哪些状态变量不可观测。目前无论线性定常系统还是分段线性定常系统，其可观测性分析方法都只能对系统状态变量的可观测情况进行定性的分析，尤其对不完全可观测系统,只能知道那些状态变量可观测,以及那些状态变量不可观测，而无法知道状态变量的可观测程度。然而状态变量的可观测程度才真正反映卡尔曼滤波器进行状态估计时的收敛速度和收敛精度。国内学者房建成等提出了一种依据系统可观测性矩阵的奇异值分解来确定状态可观测度的简单有效的方法。本发明将基于奇异值分解的可观测性分析方法用于传递对准精度评估，可以有效分析精度评估系统各状态的可观测度。

发明内容

本发明的目的在于，提供一种能够有效分析精度评估系统各状态的可观测度的基于奇异值分解的传递对准精度评估可观测性分析方法。

为实现上述目的，本发明采用的技术方案包括下列步骤：

按照以下步骤：

1)、选取精度评估系统的第一个时间段，并记为j＝1；

2)、计算A_j及H_j，并计算该时间段所对应的可观测性矩阵Q_j，计算公式为：

$Q_j=\left[\begin{array}{c}{H}_j\\ H_j{\mathrm{\Φ}}_j\\ .\\ .\\ .\\ H_j{\mathrm{\Φ}}_j^{n-1}\end{array}\right];$

3)、选取此时的SOM，并记为Q_s(j)；4)、由系统的外观测量的大小及精度，计算该时间段的外观测量Y_i；

5)、求取该时间段SOM的奇异值其中λ_i为负定矩阵A^HA的i个特征值；

6)、根据式计算出对应于每个状态变量X₀的奇异值；求取的奇异值即为状态变量可观测性的判断标准；

7)、若当前得时间段为非末时间段，则进入下一时间段精度评估系统的可观测性分析；记j＝2，重复进行步骤2)-7)，直至全部时间段的可观测性分析结束；

奇异值分解如下：

对于A∈C^m×n,非负定矩阵A^HA的i个特征值λ_i(≥0)的算数根称为A的奇异值；

令A∈R^m×n，则存在酉阵U∈R^m×n,V∈R^m×n，使得

A＝UΣV^T

其中，U＝[u₁,u₂,…,u_m]，V＝[v₁,v₂,…,v_m]为正交阵， $\mathrm{\Σ}=\left[\begin{array}{cc}S& 0\\ 0& 0\end{array}\right]$ 为m×r阶矩阵，且S＝diag(σ₁,σ₂,…σ_r)为对角阵，对角线元素σ₁≥σ₂≥…≥σ_r＞0，且r＝rank(A)；σ₁,σ₂,…,σ_r和σ_r+1＝…＝σ_n＝0称为矩阵A的奇异值；

则矩阵A的奇异值分解形式可表示为：

$A=\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{\mathrm{\σ}}_i{u}_i{v}_i^T$

所涉及的传递对准精度评估系统，其离散线性其次方程为：

$\left\{\begin{array}{c}X(k+1)=F\left(k\right)X\left(k\right)\\ Y\left(k\right)=\mathrm{HX}\left(k\right)\end{array}\right.$

将初始状态X(0)表示成一组观测值Y(0),Y(1),…,Y(k)的函数：

Y(0)＝HX(0)

Y(1)＝HF(0)X(0)

Y(2)＝HF(1)F(0)X(0)

…

$Y\left(k\right)=H\underset{j=0}{\overset{k-1}{\mathrm{\Π}}}F\left(j\right)X\left(0\right)$

$R_k=\left[\begin{array}{c}H\\ \mathrm{HF}\left(0\right)\\ \mathrm{HF}\left(1\right)F\left(0\right)\\ ...\\ H\underset{j=0}{\overset{k-1}{\mathrm{\Π}}}\left[{\mathrm{WTBX}}_{j=1}^{k-1}\right]\end{array}\right],Y=\left[\begin{array}{c}Y\left(0\right)\\ Y\left(1\right)\\ Y\left(2\right)\\ ...\\ Y\left(K\right)\end{array}\right]$

则

R_kX(0)＝Y

式中，R_k为系统的可观测性阵；

根据奇异值分解法，R_k表示为：

R_k＝UΣV^T

根据定义，式R_kX(0)＝Y改写为：

$Y=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\σ}}_i\left(v_i^T{X}_0\right)u_i$

可知，观测量Y可表示成初始状态X(0)在由[σ₁v₁,σ₂v₂,…,σ_rv_r]生成的子空间上的投影；

若σ_n＞0，初始状态X₀可由m×n个观测量估计而得

$X_0={\left({\mathrm{U\ΣV}}^T\right)}^{-1}Y=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left(\frac{u_i^{T}y}{{\mathrm{\σ}}_i}\right)v_i$

若σ₁≥σ₂≥…≥σ_r＞0，σ_r+1＝…＝σ_n＝0，则V可由两个子空间表示

V＝[V₁,V₂]

式中，V₁＝[v₁,v₂,…,v_r]，V₂＝[v_r+1,v_r+2,…,v_n]，V₂即为R_k阵的零子空间；则初始状态X₀可表示为

$X_0=\underset{i=1}{\overset{r}{\mathrm{\Σ}}}\left(\frac{u_i^{T}y}{{\mathrm{\σ}}_i}\right)v_i+\underset{j=r+1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_j{v}_j$

式中，α_j(j＝r+1,…,n)为零子空间中的任意系数；

状态变量的可观测度定义如下：

${\mathrm{\η}}_k=\frac{{\mathrm{\σ}}_i}{{\mathrm{\σ}}_0},i=\mathrm{1,2},...,n$

${\mathrm{\σ}}_i~\max \left[\frac{u_i^{T}y{v}_i}{{\mathrm{\σ}}_i}\right]$

式中，η_k为k时刻状态变量的可观测度；σ₀为外观测量相对应的奇异值；σ₁为当 $\left[\frac{u_i^{T}y{v}_i}{{\mathrm{\σ}}_i}\right](i=\mathrm{1,2},...,n)$ 取得最大时对应的奇异值；

以下重新定义σ₀：

根据最小二乘法定义σ₀为：

${\mathrm{\σ}}_0=\frac{\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{\mathrm{\σ}}_i}{4}$

使外观测量所对应的σ_i与标准σ₀之间的差值的平方和达到最小，即满足

$J\left({\mathrm{\σ}}_0\right)=\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{({\mathrm{\σ}}_i-{\mathrm{\σ}}_0)}^2=\min$

奇异值与状态变量之间的对应关系不变，因此定义：

${\tilde{X}}_0=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left(\frac{u_i^T\tilde{y}}{{\mathrm{\σ}}_i}\right)v_i$

式中，其维数与y相同。

外部基准信息包括：差分GPS速度和位置信息，主惯导水平姿态辅助差分GPS速度和位置信息，主惯导水平、方位姿态辅助差分GPS速度和位置信息；选用二通道惯导系统数学模型，其状态方程如下式：

$\left\{\begin{array}{c}\stackrel{.}{\mathrm{\φ}}=-{\mathrm{\ω}}_{\mathrm{in}}^n\×\mathrm{\φ}+\mathrm{\δ}{\mathrm{\ω}}_{\mathrm{in}}^n-C_b^n{\mathrm{\ϵ}}^b\\ \mathrm{\δ}\stackrel{.}{V}=f^n\×\mathrm{\φ}+C_b^n{\▿}^b-(2{\mathrm{\ω}}_{\mathrm{ie}}^n+{\mathrm{\ω}}_{\mathrm{en}}^n)\×\mathrm{\δV}\\ \mathrm{\δ}\stackrel{.}{P}=\mathrm{\δV}\\ \stackrel{.}{\mathrm{\ϵ}}=0\\ \stackrel{.}{\▿}=0\end{array}\right.$

式中，φ为子惯导姿态误差；为地理坐标系旋转角速率；为地理坐标系旋转角速率误差；为子惯导姿态矩阵；δV为子惯导速度误差；δP为子惯导位置误差；fⁿ为沿地理系的加速度计测量信息；ε为陀螺漂移；为加速度计零偏；为地球坐标系旋转角速率；为地理坐标系相对地球坐标系旋转角速率；机动方式包括匀速运动和加速运动。

在选择差分GPS速度和位置信息作为外部基准信息时，以子惯导与差分GPS的速度差、位置差作为滤波观测量，其量测方程为如下式：

式中，和分别为差分GPS的东向速度和北向速度；和分别为差分GPS的东向位置和北向位置；和分别为捷联捷联惯导系统的东向速度和北向速度；和分别表示捷联惯导系统的东向位置和北向位置。

在选择主惯导水平姿态辅助差分GPS速度和位置信息作为外部基准信息时，以子惯导与差分GPS系统的速度差、位置差及子主惯导水平姿态差作为滤波观测量，其量测方程为：

式中，

Z_φ'为子惯导系统的水平姿态误差量测量；和分别为捷联惯导系统的东向速度和北向速度；和分别为差分GPS的东向速度和北向速度；和分别为捷联惯导系统的东向位置和北向位置；和分别为差分GPS的东向位置和北向位置；和分别为子惯导系统的东向失准角和北向失准角；和分别为主惯导系统的东向失准角和北向失准角。

在选择主惯导水平姿态辅助差分GPS速度和位置信息作为外部基准信息时，将子惯导与差分GPS系统的速度差、位置差及子主惯导姿态差作为滤波观测量，其量测方程为：

式中，

其中：和分别为捷联惯导系统的东向速度和北向速度；和分别为差分GPS的东向速度和北向速度；和分别为捷联惯导系统的东向位置和北向位置；和分别为差分GPS的东向位置和北向位置；分别为子惯导系统的东向失准角、北向失准角、航向失准角；分别为主惯导系统的东向失准角、北向失准角、航向失准角

本发明的有益效果在于：利用奇异值分解法，在四种不同机动方式，同时引入不同参考信息的方案下，对传递对准精度评估系统进行可观测性分析，通过仿真结果明确了不同机动方式和引入不同参考信息对各姿态失准角平滑精度的影响，由此引出设计一种无需舰船做出强机动即能实现姿态失准角精确估计的精度评估方法的必要性，具有一定的工程应用价值。该方法可有效地分析精度评估系统各状态的可观测度，且可观测度越好，卡尔曼平滑器的性能越好，即姿态误差角的收敛速度越快，精度越高，可为设计传递对准精度评估方案提供有效依据。

附图说明

图1为方案(a)的可观测性直方图(第一时间段)，

图2为方案(a)的可观测性直方图(第二时间段)，

图3为方案(b)的可观测性直方图(第一时间段)，

图4为方案(b)的可观测性直方图(第二时间段)，

图5为方案(c)的可观测性直方图(第一时间段)，

图6为方案(c)的可观测性直方图(第二时间段)，

图7为方案(d)的可观测性直方图(第一时间段)，

图8为方案(d)的可观测性直方图(第二时间段)；

具体实施方式

下面结合附图1-8对本发明作进一步的详细描述。

本发明提出了一种基于奇异值分解的传递对准精度评估可观测性分析方法，该方法的按照以下步骤：

1)、选取精度评估系统的第一个时间段，并记为j＝1；

2)、计算A_j及H_j，并计算该时间段所对应的可观测性矩阵Q_j，计算公式为：

$Q_j=\left[\begin{array}{c}{H}_j\\ H_j{\mathrm{\Φ}}_j\\ .\\ .\\ .\\ H_j{\mathrm{\Φ}}_j^{n-1}\end{array}\right];$

3)、选取此时的SOM(StrippedObservabilityMatrix，提取可观测性矩阵,缩写为SOM)，并记为Q_s(j)；

4)、由系统的外观测量的大小及精度，计算该时间段的外观测量Y_i；

5)、求取该时间段SOM的奇异值其中λ_i为负定矩阵A^HA的i个特征值；

6)、根据式计算出对应于每个状态变量X₀的奇异值；求取的奇异值即为状态变量可观测性的判断标准；

7)、若当前得时间段为非末时间段，则进入下一时间段精度评估系统的可观测性分析；记j＝2，重复进行步骤2)-7)，直至全部时间段的可观测性分析结束。

奇异值分解理论如下：

对于A∈C^m×n,非负定矩阵A^HA的i个特征值λ_i(≥0)的算数根称为A的奇异值。

令A∈R^m×n，则存在酉阵U∈R^m×n,V∈R^m×n，使得

A＝UΣV^T

其中，U＝[u₁,u₂,…,u_m]，V＝[v₁,v₂,…,v_m]为正交阵， $\mathrm{\Σ}=\left[\begin{array}{cc}S& 0\\ 0& 0\end{array}\right]$ 为m×r阶矩阵，且S＝diag(σ₁,σ₂,…σ_r)为对角阵，对角线元素σ₁≥σ₂≥…≥σ_r＞0，且r＝rank(A)。σ₁,σ₂,…,σ_r和σ_r+1＝…＝σ_n＝0称为矩阵A的奇异值。

则矩阵A的奇异值分解形式可表示为：

$A=\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{\mathrm{\σ}}_i{u}_i{v}_i^T$

所涉及的传递对准精度评估系统，其离散线性其次方程为：

$\left\{\begin{array}{c}X(k+1)=F\left(k\right)X\left(k\right)\\ Y\left(k\right)=\mathrm{HX}\left(k\right)\end{array}\right.$

将初始状态X(0)表示成一组观测值Y(0),Y(1),…,Y(k)的函数：

Y(0)＝HX(0)

Y(1)＝HF(0)X(0)

Y(2)＝HF(1)F(0)X(0)

…

$Y\left(k\right)=H\underset{j=0}{\overset{k-1}{\mathrm{\Π}}}F\left(j\right)X\left(0\right)$

令

$R_k=\left[\begin{array}{c}H\\ \mathrm{HF}\left(0\right)\\ \mathrm{HF}\left(1\right)F\left(0\right)\\ ...\\ H\underset{j=0}{\overset{k-1}{\mathrm{\Π}}}\left[{\mathrm{WTBX}}_{j=1}^{k-1}\right]\end{array}\right],Y=\left[\begin{array}{c}Y\left(0\right)\\ Y\left(1\right)\\ Y\left(2\right)\\ ...\\ Y\left(K\right)\end{array}\right]$

则

R_kX(0)＝Y

式中，R_k为系统的可观测性阵。初始状态X(0)的估计是由系统SOM决定的。根据奇异值分解法，R_k表示为：

R_k＝UΣV^T

根据定义，式R_kX(0)＝Y改写为：

$Y=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\σ}}_i\left(v_i^T{X}_0\right)u_i$

可知，观测量Y可表示成初始状态X(0)在由[σ₁v₁,σ₂v₂,…,σ_rv_r]生成的子空间上的投影。因此，若想确定唯一的状态X(0)，至少需要r个观测量。

若σ_n＞0，初始状态X₀可由m×n个观测量估计而得

$X_0={\left({\mathrm{U\ΣV}}^T\right)}^{-1}Y=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left(\frac{u_i^{T}y}{{\mathrm{\σ}}_i}\right)v_i$

若σ₁≥σ₂≥…≥σ_r＞0，σ_r+1＝…＝σ_n＝0，则V可由两个子空间表示

V＝[V₁,V₂]

式中，V₁＝[v₁,v₂,…,v_r]，V₂＝[v_r+1,v_r+2,…,v_n]，V₂即为R_k阵的零子空间。则初始状态X₀可表示为

$X_0=\underset{i=1}{\overset{r}{\mathrm{\Σ}}}\left(\frac{u_i^{T}y}{{\mathrm{\σ}}_i}\right)v_i+\underset{j=r+1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_j{v}_j$

式中，α_j(j＝r+1,…,n)为零子空间中的任意系数。α_j显然具有无数解，因此初始状态X₀中的某些状态不能由m×n次观测量Y估计得到。根据工程具体情况，某些奇异值几乎为零的观测量，可将其视为零；估计效果较好的量测量则对应其相应的最大的奇异值。

由以上分析可知，系统状态X序列与可观测性阵的奇异值序列σ存在一一对应关系。其可观测性的强弱可由可观测度进行表征，以下介绍可观测度的概念。

状态变量的可观测度定义如下：

${\mathrm{\η}}_k=\frac{{\mathrm{\σ}}_i}{{\mathrm{\σ}}_0},i=\mathrm{1,2},...,n$

${\mathrm{\σ}}_i~\max \left[\frac{u_i^{T}y{v}_i}{{\mathrm{\σ}}_i}\right]$

式中，η_k为k时刻状态变量的可观测度；σ₀为外观测量相对应的奇异值；σ₁为当 $\left[\frac{u_i^{T}y{v}_i}{{\mathrm{\σ}}_i}\right](i=\mathrm{1,2},...,n)$ 取得最大时对应的奇异值。由式 ${\mathrm{\η}}_k=\frac{{\mathrm{\σ}}_i}{{\mathrm{\σ}}_0},i=\mathrm{1,2},...,n$ 可知，任意变量的可观测度，可表示为当其奇异值取最大值时与外观测量对应的奇异值之比。精度评估系统中，外观测量为速度误差δV与位置误差δP。由于要求参考系统提供的参考信息至少比子惯导系统高一个数量级，这将导致其对应的奇异值相差较大，因此外观测量σ₀的选取不唯一。针对这种情况，以下重新定义σ₀。

根据最小二乘法定义σ₀为：

${\mathrm{\σ}}_0=\frac{\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{\mathrm{\σ}}_i}{4}$

使外观测量所对应的σ_i与标准σ₀之间的差值的平方和达到最小，即满足

$J\left({\mathrm{\σ}}_0\right)=\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{({\mathrm{\σ}}_i-{\mathrm{\σ}}_0)}^2=\min$

也可省略外观测量y，直接利用矩阵A的右奇异向量，机V阵的列向量来分析可观测性。式(3-21)中，等号右边含有外观测量y。其维数随时间段的增加而增大，在相同时间段内保持不变。因此，当分析某一时间段内所有状态的可观测度时，可将外观测量y省略。所求的和与初始状态X₀虽不相等，但其奇异值与状态变量之间的对应关系不变，因此定义：

${\tilde{X}}_0=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left(\frac{u_i^T\tilde{y}}{{\mathrm{\σ}}_i}\right)v_i$

式中，其维数与y相同。该方法避免了外观测量的求取，有效简化了分析步骤，降低了计算量。

根据研究结果可知，σ_i的右奇异向量v_i，其中元素的最大绝对值分别对应于系统状态X_j。因此，仅分析v_i的列元素就可以得到相应状态的可观测度。该方法不但避免了对奇异值标准的定义，同时简化了分析步骤。

因为外部基准信息和机动方式的不同会对传递对准精度评估产生各种各样的误差，所以针对外部基准信息以及机动方式的不同，设计了四种仿真方案，其中外部基准信息包括：差分GPS速度和位置信息，主惯导水平姿态辅助差分GPS速度和位置信息，主惯导水平、方位姿态辅助差分GPS速度和位置信息。针对于舰船领域，选用二通道惯导系统数学模型，其状态方程如下式：

$\left\{\begin{array}{c}\stackrel{.}{\mathrm{\φ}}=-{\mathrm{\ω}}_{\mathrm{in}}^n\×\mathrm{\φ}+\mathrm{\δ}{\mathrm{\ω}}_{\mathrm{in}}^n-C_b^n{\mathrm{\ϵ}}^b\\ \mathrm{\δ}\stackrel{.}{V}=f^n\×\mathrm{\φ}+C_b^n{\▿}^b-(2{\mathrm{\ω}}_{\mathrm{ie}}^n+{\mathrm{\ω}}_{\mathrm{en}}^n)\×\mathrm{\δV}\\ \mathrm{\δ}\stackrel{.}{P}=\mathrm{\δV}\\ \stackrel{.}{\mathrm{\ϵ}}=0\\ \stackrel{.}{\▿}=0\end{array}\right.$

式中，φ为子惯导姿态误差；为地理坐标系旋转角速率；为地理坐标系旋转角速率误差；为子惯导姿态矩阵；δV为子惯导速度误差；δP为子惯导位置误差；fⁿ为沿地理系的加速度计测量信息；ε为陀螺漂移；为加速度计零偏；为地球坐标系旋转角速率；为地理坐标系相对地球坐标系旋转角速率。

在选择差分GPS速度和位置信息作为外部基准信息时，以子惯导与差分GPS的速度差、位置差作为滤波观测量，其量测方程为如下式：

式中，和分别为差分GPS的东向速度和北向速度；和分别为差分GPS的东向位置和北向位置；和分别为捷联捷联惯导系统的东向速度和北向速度；和分别表示捷联惯导系统的东向位置和北向位置。

在选择主惯导水平姿态辅助差分GPS速度和位置信息作为外部基准信息时，以子惯导与差分GPS系统的速度差、位置差及子主惯导水平姿态差作为滤波观测量，其量测方程为：

式中，

Z_φ'为子惯导系统的水平姿态误差量测量；和分别为捷联惯导系统的东向速度和北向速度；和分别为差分GPS的东向速度和北向速度；和分别为捷联惯导系统的东向位置和北向位置；和分别为差分GPS的东向位置和北向位置；和分别为子惯导系统的东向失准角和北向失准角；和分别为主惯导系统的东向失准角和北向失准角。

在选择主惯导水平姿态辅助差分GPS速度和位置信息作为外部基准信息时，将子惯导与差分GPS系统的速度差、位置差及子主惯导姿态差作为滤波观测量，其量测方程为：

式中，

其中：和分别为捷联惯导系统的东向速度和北向速度；和分别为差分GPS的东向速度和北向速度；和分别为捷联惯导系统的东向位置和北向位置；和分别为差分GPS的东向位置和北向位置；分别为子惯导系统的东向失准角、北向失准角、航向失准角；分别为主惯导系统的东向失准角、北向失准角、航向失准角。

仿真方案设置：

通常在进行传递对准的精度评估时，机动方式有两种：匀速运动和加速运动。因此，以下将系统按机动方式的不同分为两个时间段。鉴于现有舰船领域惯导系统，其传递对准精度评估方法大都采用基于主惯导姿态信息的方法，解决其评估性能受限于舰船有限机动能力的问题，现针对引入机动与基于主惯导姿态信息的精度评估系统可观测性展开分析。

具体的方案描述为如表1所示

表1 可观测分析的方案描述

Matlab仿真条件设置：

滤波器初始值设置

Q＝diag{(0.05°/h)²,(0.05°/h)²,(0.05°/h)²,

(50μg)²,(50μg)²,(50μg)²,0,0,0,0,0,0,0}

传统机动方案的量测噪声矩阵：

R＝diag{(0.01m/s)²,(0.01m/s)²,(6m)²,(6m)²}

基于主惯导“水平”姿态信息的量测噪声矩阵：

R＝diag{(0.01m/s)²,(0.01m/s)²,(6m)²,(6m)²,(6')²,(6')²,(6')²}

基于主惯导“水平+方位”姿态信息的量测噪声矩阵：

R＝diag{(0.01m/s)²,(0.01m/s)²,(0.01m/s)²,

(6m)²,(6m)²,(6')²,(6')²,(6')²}

仿真结果：

为简化分析，利用系统的SOM代替TOM来进行可观测性分析。根据上述PWCS和奇异值分解理论，由所建立的精度评估数学模型分别设计各方案的MATLAB仿真程序，得到四种方案下系统13个状态变量的可观测度直方图。其中方案(a)在第一时间段和第二时间段的可观测度直方图分别如图1、2所示。值得注意的是，为简化处理，本文将忽略小于10^-4的奇异值，将其视为零。

在每个直方图中，状态变量的序列为：

${\mathrm{\φ}}_x,{\mathrm{\φ}}_y,{\mathrm{\φ}}_z,{\mathrm{\δV}}_x,{\mathrm{\δV}}_y,{\mathrm{\δP}}_x,{\mathrm{\δP}}_y,{\mathrm{\ϵ}}_{\mathrm{bx}}^b,{\mathrm{\ϵ}}_{\mathrm{by}}^b,{\mathrm{\ϵ}}_{\mathrm{bz}}^b,{\tilde{N}}_{\mathrm{bx}}^b,{\tilde{N}}_{\mathrm{by}}^b,{\tilde{N}}_{\mathrm{bz}}^b.$

由图1可知，系统在第一时间段SOM的秩为9，系统为不完全可观测的，其中有4个不可观测变量。为简化分析，令外观测量所对应的奇异值为1。因此，每个状态变量的可观测度等于其奇异值。直方图的读取方式为：每个直方图中，柱状条绝对值为最大时，其数字对应于状态变量序列中相应的状态变量。那么这个状态变量的可观测度即为柱状条最大时对应的奇异值。因此，φ_x、φ_y的可观测度为13.9417；的可观测度为13.8697；δV_x、δV_y的可观测度为1.4142；δP_x、δP_y的可观测度为1；的可观测度为0.0007；的可观测度为0；因此，的可观测度相对较低；而四个状态完全不可观测。

在第二时间段，系统SOM的秩为11，有2个不可观测变量； 的可观测性有所提高；而由不可观测变为可观测，但的可观测性相对较低；而状态变量仍不可观测。

方案(b)在第一时间段和第二时间段的可观测性直方图分别如图3、4所示。

方案(b)在第一时间段内，不可观测的变量有4个。在第二时间段内，方位失准角φ_z由不可观测变为可观测，但其可观测度相对较低；不可观测。

方案(c)在第一时间段和第二时间段的可观测性直方图分别如图5、6所示。

方案(c)在第一时间段时，不可观测的状态变量有2个，分别为φ_z和且的可观测度依然很低。在第二时间段内，所有13个状态都可观测；但φ_z的可观测度为0.0008，几乎不可观测。

方案(d)在第一时间段和第二时间段的可观测性直方图分别如图7、8所示。

方案(d)在第一时间段时，系统SOM秩为12，其中有一个不可观测变量为在第二时间段内，系统13个状态变量都可观测；且方位失准角φ_z的可观测度为1.4142，说明可观测性较好，系统可以对方位失准角进行精确估计。

本文为突出研究重点，仅对水平姿态误差和方位姿态误差进行可观测度的对比研究。由于东向失准角与北向失准角的可观测度相近，因此以下仅列出东向失准角的可观测度，各方案的可观测性结果如表2所示。

表2 各方案的可观测性结果

由表2可知，四种方案中φ_x的奇异值均较大，说明φ_x的可观测性较好，可以得到精确的估计；而在第一时间段内，对于方位失准角φ_z，只有方案(d)中所设计的引入主惯导水平和方位信息的方法可以实现对φ_z的估计。在整个时间段内，方案(a)及方案(c)的可观测度几乎都为零，无法完成对方位失准角φ_z的估计；方案(b)由于进行加速机动，使得φ_z的奇异值增大，可观测性好，但由于加速度的绝对值相对较小，所以对方位失准角的估计精度依然很差；方案(d)通过引入主惯导水平及方位信息，使得方位失准角的可观测度上升为1.4142，说明在外界参考信息高精度的情况下，可以精确地对方位失准角进行估计。

结合上述分析，得到如下的分析结果：通过本发明所提出的基于奇异值分解的传递对准精度评估可观测性分析方法，可有效分析精度评估系统各状态的可观测度。同时，可观测度越好，卡尔曼平滑器的性能越好，即姿态误差角的收敛速度越快、精度越高，可为设计传递对准精度评估方案提供有效依据。应理解，这些实施例仅用于说明本发明而不用于限制本发明的范围。此外应理解，在阅读了本发明讲授的内容之后，本领域技术人员可以对本发明作各种改动或修改，这些等价形式同样落于本申请所附权利要求书所限定的范围。

Claims (5)

1.一种基于奇异值分解的传递对准精度评估可观测性分析方法，其特征在于，按照以下步骤：

1)、选取精度评估系统的第一个时间段，并记为j＝1；

2)、计算A_j及H_j，并计算该时间段所对应的可观测性矩阵Q_j，计算公式为：

$Q_j=\left[\begin{array}{c}{H}_j\\ H_j{\mathrm{\Φ}}_j\\ \·\\ \·\\ \·\\ H_j{\mathrm{\Φ}}_j^{n-1}\end{array}\right];$

3)、选取此时的SOM，并记为Q_s(j)；4)、由系统的外观测量的大小及精度，计算该时间段的外观测量Y_i；

5)、求取该时间段SOM的奇异值其中λ_i为负定矩阵A^HA的i个特征值；

6)、根据式i＝1,2,...,n，计算出对应于每个状态变量X₀的奇异值；求取的奇异值即为状态变量可观测性的判断标准；

7)、若当前得时间段为非末时间段，则进入下一时间段精度评估系统的可观测性分析；记j＝2，重复进行步骤2)-7)，直至全部时间段的可观测性分析结束；

奇异值分解如下：

对于A∈C^m×n,非负定矩阵A^HA的i个特征值λ_i(≥0)的算数根称为A的奇异值；

令A∈R^m×n，则存在酉阵U∈R^m×n,V∈R^m×n，使得

A＝UΣV^T

其中，U＝[u₁,u₂,...,u_m]，V＝[v₁,v₂,...,v_m]为正交阵， $\mathrm{\Σ}=\left[\begin{array}{cc}S& 0\\ 0& 0\end{array}\right]$ 为m×r阶矩阵，且S＝diag(σ₁,σ₂,...σ_r)为对角阵，对角线元素σ₁≥σ₂≥...≥σ_r＞0，且r＝rank(A)；σ₁,σ₂,...,σ_r和σ_r+1＝...＝σ_n＝0称为矩阵A的奇异值；

则矩阵A的奇异值分解形式可表示为：

$A=\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{\mathrm{\σ}}_i{u}_i{v}_i^T$

所涉及的传递对准精度评估系统，其离散线性其次方程为：

$\left\{\begin{array}{c}X(k+1)=F\left(k\right)X\left(k\right)\\ Y\left(k\right)=\mathrm{HX}\left(k\right)\end{array}\right.$

将初始状态X(0)表示成一组观测值Y(0),Y(1),...,Y(k)的函数：

$\begin{array}{c}Y\left(0\right)=\mathrm{HX}\left(0\right)\\ Y\left(1\right)=\mathrm{HF}\left(0\right)X\left(0\right)\\ Y\left(2\right)=\mathrm{HF}\left(1\right)F\left(0\right)X\left(0\right)\\ \·\·\·\\ Y\left(k\right)=H\underset{j=0}{\overset{k-1}{\mathrm{\Π}}}F\left(j\right)X\left(0\right)\end{array}$

令

$R_k=\left[\begin{array}{c}H\\ \mathrm{HF}\left(0\right)\\ \mathrm{HF}\left(1\right)F\left(0\right)\\ \·\·\·\\ H\underset{j=0}{\overset{k-1}{\mathrm{\Π}}}\left[{\mathrm{WTBX}}_{j=1}^{k-1}\right]\end{array}\right],Y=\left[\begin{array}{c}Y\left(0\right)\\ Y\left(1\right)\\ Y\left(2\right)\\ \·\·\·\\ Y\left(K\right)\end{array}\right]$

则

R_kX(0)＝Y

式中，R_k为系统的可观测性阵；

根据奇异值分解法，R_k表示为：

R_k＝UΣV^T

根据定义，式R_kX(0)＝Y改写为：

$Y=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\σ}}_i\left(v_i^T{X}_0\right)u_i$

可知，观测量Y可表示成初始状态X(0)在由[σ₁v₁,σ₂v₂,...,σ_rv_r]生成的子空间上的投影；

若σ_n＞0，初始状态X₀可由m×n个观测量估计而得

$X_0={\left({\mathrm{U\ΣV}}^T\right)}^{-1}Y=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left(\frac{u_i^{T}y}{{\mathrm{\σ}}_i}\right)v_i$

若σ₁≥σ₂≥...≥σ_r＞0，σ_r+1＝...＝σ_n＝0，则V可由两个子空间表示

V＝[V₁,V₂]

式中，V₁＝[v₁,v₂,...,v_r]，V₂＝[v_r+1,v_r+2,...,v_n]，V₂即为R_k阵的零子空间；则初始状态X₀可表示为

$X_0=\underset{i=1}{\overset{r}{\mathrm{\Σ}}}\left(\frac{u_i^{T}y}{{\mathrm{\σ}}_i}\right)v_i+\underset{j=r+1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_j{v}_j$

式中，α_j(j＝r+1,...,n)为零子空间中的任意系数；

状态变量的可观测度定义如下：

${\mathrm{\η}}_k=\frac{{\mathrm{\σ}}_i}{{\mathrm{\σ}}_0},$ i＝1,2,...,n

${\mathrm{\σ}}_i~\max \left[\frac{u_i^{T}y{v}_i}{{\mathrm{\σ}}_i}\right]$

式中，η_k为k时刻状态变量的可观测度；σ₀为外观测量相对应的奇异值；σ₁为当(i＝1,2,...,n)取得最大时对应的奇异值；

以下重新定义σ₀：

根据最小二乘法定义σ₀为：

${\mathrm{\σ}}_0=\frac{\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{\mathrm{\σ}}_i}{4}$

使外观测量所对应的σ_i与标准σ₀之间的差值的平方和达到最小，即满足

$J\left({\mathrm{\σ}}_0\right)=\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{({\mathrm{\σ}}_i-{\mathrm{\σ}}_0)}^2=\min$

奇异值与状态变量之间的对应关系不变，因此定义：

${\tilde{X}}_0=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left(\frac{u_i^T\tilde{y}}{{\mathrm{\σ}}_i}\right)v_i$

式中，其维数与y相同。

2.根据权利要求1所述的一种基于奇异值分解的传递对准精度评估可观测性分析方法，其特征在于，外部基准信息包括：差分GPS速度和位置信息，主惯导水平姿态辅助差分GPS速度和位置信息，主惯导水平、方位姿态辅助差分GPS速度和位置信息；选用二通道惯导系统数学模型，其状态方程如下式：

$\left\{\begin{array}{c}\stackrel{\·}{\mathrm{\φ}}=-{\mathrm{\ω}}_{\mathrm{in}}^n\×\mathrm{\φ}+\mathrm{\δ}{\mathrm{\ω}}_{\mathrm{in}}^n-C_b^n{\mathrm{\ϵ}}^b\\ \mathrm{\δ}\stackrel{\·}{V}=f^n\×\mathrm{\φ}+C_b^n{\▿}^b-(2{\mathrm{\ω}}_{\mathrm{in}}^n+{\mathrm{\ω}}_{\mathrm{in}}^n)\×\mathrm{\δV}\\ \mathrm{\δ}\stackrel{\·}{P}=\mathrm{\δV}\\ \stackrel{\·}{\mathrm{\ϵ}}=0\\ \stackrel{\·}{\▿}=0\end{array}\right.$

式中，φ为子惯导姿态误差；为地理坐标系旋转角速率；为地理坐标系旋转角速率误差；为子惯导姿态矩阵；δV为子惯导速度误差；δP为子惯导位置误差；fⁿ为沿地理系的加速度计测量信息；ε为陀螺漂移；为加速度计零偏；为地球坐标系旋转角速率；为地理坐标系相对地球坐标系旋转

角速率；机动方式包括匀速运动和加速运动。

3.根据权利要求2所述的一种基于奇异值分解的传递对准精度评估可观测性分析方法，其特征在于，在选择差分GPS速度和位置信息作为外部基准信息时，以子惯导与差分GPS的速度差、位置差作为滤波观测量，其量测方程为如下式：

式中，和分别为差分GPS的东向速度和北向速度；和分别为差分GPS的东向位置和北向位置；和分别为捷联捷联惯导系统的东向速度和北向速度；和分别表示捷联惯导系统的东向位置和北向位置。

4.根据权利要求2所述的一种基于奇异值分解的传递对准精度评估可观测性分析方法，其特征在于，在选择主惯导水平姿态辅助差分GPS速度和位置信息作为外部基准信息时，以子惯导与差分GPS系统的速度差、位置差及子主惯导水平姿态差作为滤波观测量，其量测方程为：

式中，

Z_φ'为子惯导系统的水平姿态误差量测量；和分别为捷联惯导系统的东向速度和北向速度；和分别为差分GPS的东向速度和北向速度；和分别为捷联惯导系统的东向位置和北向位置；和分别为差分GPS的东向位置和北向位置；和分别为子惯导系统的东向失准角和北向失准角；和分别为主惯导系统的东向失准角和北向失准角。

5.根据权利要求2所述的一种基于奇异值分解的传递对准精度评估可观测性分析方法，其特征在于，在选择主惯导水平姿态辅助差分GPS速度和位置信息作为外部基准信息时，将子惯导与差分GPS系统的速度差、位置差及子主惯导姿态差作为滤波观测量，其量测方程为：

式中，

其中：和分别为捷联惯导系统的东向速度和北向速度；和分别为差分GPS的东向速度和北向速度；和分别为捷联惯导系统的东向位置和北向位置；和分别为差分GPS的东向位置和北向位置；分别为子惯导系统的东向失准角、北向失准角、航向失准角；分别为主惯导系统的东向失准角、北向失准角、航向失准角。