keywords:
误差分配研究得比较早了，18年的时候学习过一段时间，当时一直用严老师提出的误差分配与可观测度分析方法。这套误差分配方法还用在了一些项目上，算法理论虽然简单，但是推导过程有点复杂。
这个方法实验室一直有人用，主要用在初始对准、误差传递分析和旋转惯导研究上，经常一直讨论这个方法，发现了一些问题，一直没有整理。恰好想投CCC2022的会议，就重新整理了一下，主要修改了算法的推导思路和过程，我觉得修改后更加容易理解了。
<!--more-->
目录：
1 什么是SINS误差分配？
1.1 设计一个导航系统
1.2 选取合适的组合导航状态量
1.3 Mont-Carlo仿真分析
2 协方差分解
2.1 时间更新分析
2.2 按照状态分解
2.3 协方差分解流程
3 SINS误差分配
3.1 SINS误差模型
3.2 误差分配流程
3.3 误差分配软件
4 开源项目与附件
1 什么是SINS误差分配？
1.1 选取合适的组合导航状态量
组合导航中往往需要根据不同的应用条件选取合适的状态向量，选择的状态向量如果可观测性较低，则有可能影响组合导航精度，增加组合导航的计算量。
常见的15维组合导航状态向量如下:
\boldsymbol{X} = \left[ \begin{matrix} \boldsymbol{\phi} & \delta\boldsymbol{v} & \delta\boldsymbol{p} & \boldsymbol{\varepsilon}^{b} & \boldsymbol{\nabla}^{b} \end{matrix} \right]
在不同的使用场景，可以对状态向量进行修改，增加影响较大的误差项作为状态，去除对导航结果没有影响的状态，增加系统可靠性。
那么要如何对每个状态量，特别是导航误差以外的误差项进行分析，选择出合适的导航向量呢？
如果有一张表可以表示各个误差源在导航误差中的占比，那就可以轻松选择状态向量。
| 导航误差 | 陀螺零偏 | 加计零偏 | 角度随机游走 | 加计随机游走 |
|---|
| 东向失准角 | 30% | 20% | 10% | 9% |
| 水平位置误差 | 50% | 30% | 6% | 3% |
1.2 状态可观测度分析
在SINS系统级标定中，我们希望所有的IMU误差项都可以被辨识出来，因此需要设计一套转动方案，使得完成要求的转动后，所有误差项全部可观测。状态的可观测度就与相应误差项在量测信息中的占比有关，如果在量测中占比较小，甚至引起的量测变化小于量测噪声，则很难将其辨识出来。
转台的转动种类有限，无非就是绕三轴的转动以及在6个位置的静止状态，以此为基础设计可以辨识全部IMU误差参数的转动方案就需要先对各种基础转动下导航误差中各个误差源的占比进行分析，保证经过一系列转动后，所有的误差项在量测的导航误差中均有较大影响，就能通过Kalman滤波估计出IMU误差参数。
我在硕士论文中也给出了如何设计系统级标定转动方案，第一步就是分析绕各轴转动不同角度，不同误差项的可观测度，而可观测度就与误差占比有关。
1.3 Monte-Carlo仿真分析
最常见的误差传递分析方法是Monte-Carlo仿真，因为惯导随机误差的存在，想要获得导航误差规律，必然需要进行多次仿真统计导航结果。文献都说Monte-Carlo仿真可以用于误差分配，但是不太理解如果研究各个误差源在最终导航误差中的占比，一直也没有找到相关文献资料。
我的想法是这样的，首先设计一条长航时轨迹，并设置完整的误差参数；然后利用Monte-Carlo仿真得到导航误差；最后将被分析误差项置零，利用Monte-Carlo仿真得到导航误差，与之前完成误差下导航误差的差值就可以认为是该误差源造成的导航误差。
看着很有道理，也不知道对不对，也没有仿真过。
2 协方差分解
这里简单列一下推导的结论，具体的过程可以看论文，网站上写了应当也没人看
2.1 时间更新分析
系统状态方程的连续形式和离散形式如下：
\boldsymbol{\dot{X}} = \boldsymbol{FX}+\boldsymbol{GU}
\boldsymbol{X}_{k} = \boldsymbol{\Phi}_{k/k-1} \boldsymbol{X}_{k-1} + \boldsymbol{\Gamma}_{k-1} \boldsymbol{U}_{k-1}
按照时间进行更新的过程为：
\left\{ \begin{aligned}
&\boldsymbol{X}_1 = \boldsymbol{X}_{1/0} \boldsymbol{X}_0 + \boldsymbol{\Gamma}_0 \boldsymbol{U}_0 \\
&\boldsymbol{X}_2 = \boldsymbol{X}_{2/1} \boldsymbol{X}_1 + \boldsymbol{\Gamma}_1 \boldsymbol{U}_1 \\
\vdots \\
&\boldsymbol{X}_n = \boldsymbol{X}_{n/n-1} \boldsymbol{X}_{n-1} + \boldsymbol{\Gamma}_{n-1} \boldsymbol{U}_{n-1}
\end{aligned} \right.
整理后就得到:
\begin{aligned}
\boldsymbol{X}_{n} =&\boldsymbol{\Phi}_{n \backslash n-1} \boldsymbol{X}_{n-1}+\boldsymbol{\Gamma}_{n-1} \boldsymbol{U}_{n-1} \\
=&\boldsymbol{\Phi}_{n \backslash n-1}\left(\boldsymbol{\Phi}_{n-1 \backslash n-2} \boldsymbol{X}_{n-2}+\boldsymbol{\Gamma}_{n-2} \boldsymbol{U}_{n-2}\right)+ \\
&\boldsymbol{\Gamma}_{n-1} \boldsymbol{U}_{n-1} \\
=&\boldsymbol{\Phi}_{n \backslash n-1} \boldsymbol{\Phi}_{n-1 \backslash n-2} \boldsymbol{X}_{n-2}+ \\
&\boldsymbol{\Phi}_{n \backslash n-1} \boldsymbol{\Gamma}_{n-2} \boldsymbol{U}_{n-2}+ \boldsymbol{\Gamma}_{n-1} \boldsymbol{U}_{n-1} \\
& \cdots \\
=&\prod_{i=1}^{n} \boldsymbol{\Phi}_{i \backslash i-1} \boldsymbol{X}_{0}+\sum_{j=0}\left(\prod_{k=j}^{n-2} \boldsymbol{\Phi}_{k+2 \backslash k+1} \boldsymbol{\Gamma}_{k} \boldsymbol{U}_{k}\right) \\
&+\boldsymbol{\Gamma}_{n-1} \boldsymbol{U}_{n-1}
\end{aligned}
如果输入向量是随机误差，就需要用协方差更新公式表示，即：
\begin{aligned}
\boldsymbol{P}_{n} =& \boldsymbol{\Phi}_{n/n-1} \boldsymbol{P}_{n-1} \boldsymbol{P}_{n-1}^{\mathrm{T}} + \boldsymbol{\Gamma}_{n-1} \boldsymbol{Q}_{n-1} \boldsymbol{\Gamma}_{n-1} \\
=& \prod_{i=1}^{n}\boldsymbol{\Phi}_{i/i-1}\boldsymbol{P}_{0}\boldsymbol{\Phi}_{i/i-1}^{\mathrm{T}} + \\
& \sum_{j=0}^{n-2} \left( \prod_{k=j}^{n-2} \boldsymbol{\Phi}_{k+2/k+1} \boldsymbol{\Gamma}_{k} \boldsymbol{Q}_{k} \boldsymbol{\Gamma}_{k}^{\mathrm{T}} \boldsymbol{\Phi}_{k+2/k+1}^{\mathrm{T}} \right)+ \\
& \boldsymbol{\Gamma}_{n-1} \boldsymbol{Q}_{n-1} \boldsymbol{\Gamma}_{n-1}^{\mathrm{T}}
\end{aligned}
这几步公式非常简单，就是为了说明任意时刻的导航结果都是初始时刻的状态和之前的全部输入确定的。这其实是不用分析的，但是因为误差分配需要这样分解，得到新的状态更新形式：
\left\{ \begin{aligned}
& \boldsymbol{\bar{X}}_{k+1} = \boldsymbol{\Phi}_{k+1/k}\boldsymbol{\bar{X}}_{k} \\
& \boldsymbol{\bar{U}}_{k+1} = \boldsymbol{\Phi}_{k+1/k}\boldsymbol{\bar{U}}_{k} + \boldsymbol{\Gamma}_{k} \boldsymbol{\bar{U}}_k \\
& \boldsymbol{X}_{k+1} = \boldsymbol{\bar{X}}_{k+1} + \boldsymbol{\bar{U}}_{k+1}
\end{aligned}\right. \space \space
\left( \begin{gathered}
\boldsymbol{\bar{X}}_0 = \boldsymbol{X}_0 \\
\boldsymbol{\bar{U}}_0 = \boldsymbol{0} \\
k=0,1,\cdots,n-1
\end{gathered}\right)
\begin{gathered}
\left\{ \begin{aligned}
&\bar{\boldsymbol{P}}_{k+1} = \boldsymbol{\Phi}_{k+1/k} \bar{\boldsymbol{P}}_{k} \boldsymbol{\Phi}_{k+1/k}^{\mathrm{T}} \\
&\bar{\boldsymbol{Q}}_{k+1} = \boldsymbol{\Phi}_{k/k-1} \bar{\boldsymbol{Q}}_{k} \boldsymbol{\Phi}_{k/k-1}^{\mathrm{T}} + \boldsymbol{\Gamma}_{k} \boldsymbol{Q}_{k} \boldsymbol{\Gamma}_{k}^{\mathrm{T}} \\
&\boldsymbol{P}_{k+1} = \bar{\boldsymbol{P}}_{k+1}+\bar{\boldsymbol{Q}}_{k+1}
\end{aligned}\right. \space,
\left( \begin{gathered}
\bar{\boldsymbol{P}}_{0} = \boldsymbol{P}_{0} \\
\bar{\boldsymbol{Q}}_{0} = \boldsymbol{0}
\end{gathered} \right) \\
\left(k=0,1,\cdots,n-1\right)
\end{gathered}
2.2 按照状态分解
对于状态向量，按照状态分解过程为：
\boldsymbol{\bar{X}}_{k} = \boldsymbol{\bar{X}}_{k,(1)} + \boldsymbol{\bar{X}}_{k,(2)} + \cdots + \boldsymbol{\bar{X}}_{k,(n)}
\left\{ \begin{aligned}
&\boldsymbol{\bar{X}}_{0,(1)} = \left[\begin{matrix} X_{0,(1)} & \boldsymbol{0}_{n-1\times1} \end{matrix}\right]^{\mathrm{T}} \\
\vdots \\
&\boldsymbol{\bar{X}}_{0,(i)} = \left[\begin{matrix} \boldsymbol{0}_{i-1\times1} & X_{0,(i)} & \boldsymbol{0}_{n-i\times1} \end{matrix}\right]^{\mathrm{T}} \\
\vdots \\
&\boldsymbol{\bar{X}}_{0,(n)} = \left[\begin{matrix} \boldsymbol{0}_{n-1\times1} & X_{0,(n)} \end{matrix}\right]^{\mathrm{T}}
\end{aligned} \right.
同理，协方差也可以分解。类似地，可以对协方差按照对角元素进行分类，即通过将原来的一个误差协方差矩阵按照误差源分解为多个矩阵求和的形式。
\bar{\boldsymbol{P}}_{k} = \bar{\boldsymbol{P}}_{k,(1)} + \bar{\boldsymbol{P}}_{k,(2)} + \cdots + \bar{\boldsymbol{P}}_{k,(p)}
这样就可以得到按照分解后状态进行更新的系统方程：
\begin{aligned}
&\left\{ \begin{aligned}
& \boldsymbol{\bar{X}}_{k+1,(i)} = \boldsymbol{\Phi}_{k+1/k}\boldsymbol{\bar{X}}_{k,(i)} \\
& \boldsymbol{\bar{U}}_{k+1,(j)} = \boldsymbol{\Phi}_{k+1/k}\boldsymbol{\bar{U}}_{k},(j) + \boldsymbol{\Gamma}_{k} \boldsymbol{\bar{U}}_{k,(j)} \\
& \boldsymbol{X}_{k+1} = \sum_{i=1}^{n}\boldsymbol{\bar{X}}_{k+1,(i)} + \sum_{j=1}^{m}\boldsymbol{\bar{U}}_{k+1,(j)}
\end{aligned}\right.
\left( \begin{gathered}
\boldsymbol{\bar{U}}_{0,(j)} = \boldsymbol{0}
\end{gathered}\right) \\
&\left( \begin{gathered}
i=1,2,\cdots,n; j=1,2,\cdots,m; k=0,1,\cdots,n-1
\end{gathered}\right)
\end{aligned}
\begin{gathered}
\left\{ \begin{aligned}
&\bar{\boldsymbol{P}}_{k+1,(i)} = \boldsymbol{\Phi}_{k+1/k} \bar{\boldsymbol{P}}_{k,(i)} \boldsymbol{\Phi}_{k+1/k}^{\mathrm{T}} \\
&\bar{\boldsymbol{Q}}_{k+1,(j)} = \boldsymbol{\Phi}_{k/k-1} \bar{\boldsymbol{Q}}_{k,(j)} \boldsymbol{\Phi}_{k/k-1}^{\mathrm{T}} + \boldsymbol{\Gamma}_{k} \boldsymbol{Q}_{k,(j)} \boldsymbol{\Gamma}_{k}^{\mathrm{T}} \\
&\boldsymbol{P}_{k+1} = \sum_{i=1}^{p}\bar{\boldsymbol{P}}_{k+1,(j)} + \sum_{j=1}^{m}\bar{\boldsymbol{Q}}_{k+1,(j)}
\end{aligned}\right. \\
(i=1,\cdots,p; j=1,\cdots,m; k=0,1,\cdots, n-1)
\end{gathered}
$\bar{\boldsymbol{P}}{k+1,(j)} $和$\bar{\boldsymbol{Q}}{k+1,(j)} $的每一个元素都包含丰富的SINS误差传递信息，后文会进行数据重组，得到误差分配结果。
以上分解原理很简单，但是非常繁琐，其主要目的就是将原有的误差传递信息重新按照误差源类型、导航误差类型和时间进行信息重组。通过对新的分解后的协方差矩阵进行数据重构，挖掘出导航误差分配信息。 为了方便理解，重构后的数据就如下图所示：
每一个分解后的协方差矩阵表示一个（类）误差源对导航误差的影响，其对角线元素表示其造成的导航误差方差。按照论文中的分类方法，取前7个对角元素作为误差分配分析的数据基础。
提取出每一个数字都对应某一时刻某误差源对某导航误差的作用。
完成对角元素提取后，就可以将数据整理，得到一个三维数据集，三个维度分别表示：时间、导航误差和误差源。同一时间下，一个导航误差中不同误差源的占比就是误差分配清单，可以将其表示成百分比的形式。
3 SINS误差分配
3.1 SINS误差模型
捷联惯导误差模型没有什么需要讲的，为了使误差分配的结果比较丰富，不要仅仅只分析常值零偏的影响，在论文中还加入了陀螺和加表的刻度系数误差以及安装偏角误差。仿真使用的误差参数表如下：
3.2 误差分配流程
这里的流程指的是仿真测试流程，具体的步骤如下流程图：
除了最后一步是仿真分析的步骤，其他都是误差分配的计算过程。
3.3 误差分配软件
误差分配的基本原理非常简单，就是过程繁琐，如果每个研究者都去完成矩阵拆分、重组、捷联惯导误差建模是没有意义的。误差分配更多地是进行快速定性分析，所以我就趁自己闲下来的时间（不想写文章的时间），写了一个误差分配的MatlabApp，同时也是丰富一下我的gitee内容。
这个软件非常简单，导入轨迹数据，设置误差参数后就可以进行误差分配计算，得到三张误差分配图，分别对应姿态误差、速度误差和位置误差。轨迹数据当然是按照严老师PSINS工具箱的trj数据格式，这样分析的时候，直接调用以前的仿真数据就可以了，省去了跑轨迹仿真的时间。软件界面如下：
说到轨迹仿真，其实放程序的gitee仓最开始就是为轨迹仿真程序设置的，后来导航相关的程序多了又不想创建太多的仓，就都放在了一起，用文件夹区分。
5 开源项目与附件
论文下载：
[arXiv]Analysis Method of Strapdown Inertial Navigation Error Distribution Based on Covariance Matrix Decomposition
点击下载PDF
论文报告ppt:
【石墨文档-Error Distribution Analysis Method】
误差分类程序：
[Gitee]SINS误差分配分析程序
这个ppt做的非常不认真，主要原因就是这届CCC办得太差了，提前半个月我就放弃了ppt制作，本来准备好好做20多页。最后很多都删掉了。
