【经典文献翻译】IMU Error Modeling Tutorial

Catalogue
  1. 1. IMU Error Modeling Tutorial (INS state estimation with real-time sensor calibration)
  2. 2. 介绍
  3. 3. 应用简介
  4. 4. State Estimation Error Model
  5. 5. 线性状态空间系统的功率谱密度
  6. 6. Available Error Specification Information
  7. 7. 问题描述
  8. 8. Allan Variance
  9. 9. Power Spectral Density and Allan Variance
  10. 10. Modeling via Independent Noise Sources
  11. 11. Continuous-Time State-Space Models
    1. 11.1. 随机游走误差\(z_N(t)\):角速度和速度
    2. 11.2. 随机游走误差\(z_K(t)\):速率和加速度
    3. 11.3. 累积误差模型:N,K
    4. 11.4. 偏差不稳定性\(z_B(t)\)
    5. 11.5. Gauss-Markov Error Model
    6. 11.6. Cumulative Error Model: N , B, K

IMU Error Modeling Tutorial (INS state estimation with real-time sensor calibration)

Clipboard 2023年3月11日 10.14

介绍

自动驾驶汽车技术正在快速发展。关键的使能因素是计算和传感系统能力的提高和成本的降低,这些系统使传感融合能够感知车辆的状态和环境。出于控制目的,必须以足够高的采样率和足够高的带宽准确、可靠地估计车辆状态。对于具有高带宽的系统,这些要求通常通过辅助惯性导航系统(INS)[1]、[2]、[3]、[4]、[5]、[6]来实现。

NS通过IMU高采样率的运动学模型集成来自惯性测量单元(IMU)的6个数据,以计算状态估计。辅助INS使用从辅助传感器(如视觉、激光雷达、雷达和全局导航卫星系统(GNSS))数据纠正这一状态估计。传感器融合状态估计可以通过以下任何一种方法来完成:卡尔曼滤波器(KF)[7],[8],[9],[10],[11],扩展卡尔曼滤波器(EKF)[12],[13],[14],[15],无迹卡尔曼滤波器(UKF)[16],[17],[18],粒子滤波(PF)[19],[20],[21],12和最大后验(MAP)优化[22],[23],[24],[25],[26],[27]。

如果将IMU和辅助传感器数据结合在状态空间形式中,则结合IMU和辅助传感器数据的数据融合系统将能够通过实时校准实现更好的性能。IMU制造商提供了一个描述预期IMU性能的数据表。根据规范标准[28]、[29]、[30],这种性能通常用艾伦方差(AV)来表示。

然而,目前还不清楚如何将这些数据表中的AV信息转换成合适的状态空间模型。多种模型已经被知晓并使用了几十年[31],[32],[33],[34],[35],20[36],[37],[38],[39],[40]。尽管它们很重要,但现有文献中没有对潜在思想、问题和权衡进行清晰的教程阐述。本文的目的是提供这样一个指导性的讨论。“辅助INS历史”中讨论了这些想法的悠久历史以及与成功应用相关的问题。

本教程讨论与示例导航系统设计方法相关的问题和权衡:

  • (1)使用AV信息来指定IMU状态空间随机误差模型的连续时间参数(例如,\(p, n_z, A_z, B_z, C_z, S_z\), and \(S_\eta\));
  • (2)将该连续时间模型转换为实现状态估计器所需的离散时间、状态空间误差模型参数(\(\Phi, Q_{z_d}, H, Q_{\eta_d}\))
  • (3)相对于Allan信息验证IMU模型。以前的一些文章讨论了上述专题中的一些[32]、[34]、[36]、[39]、[43]、[44]、[46]、[47]。
  • 本文的目标是以教程的方式清楚而全面地介绍背景和主要思想,使用符合仪器规范标准[28]、[29]的符号和术语。为了澄清问题,通篇都包括了一些例子。

应用简介

INS设计是基于车辆运动学模型:

Clipboard 2023年3月11日 10.27

其中,\(\vec{x}_v\)表示车辆的状态\(\vec{u} \in \Re^6\)表示系统输入(即比力和角速度向量)。典型的车辆状态向量可能包括位置、速度和姿态的子向量。导航系统基于信号\(\vec{u}(t)\)的测量来求解(1):

Clipboard 2023年3月11日 10.30

其中 \(\hat{\vec{u}}(t)\)是使用实时估计的校准因子(从IMU 测量\(\tilde{\vec{u}}(t)\)计算出来的)。“"Simplified INS Example.”中介绍了简化的二维惯性导航示例。

对于标量信号,将传感器测量̃\(\tilde{u}(t)\)与期望真值信号\(u(t)\)相关联的模型是(参见[28]、[29]):

Clipboard 2023年3月11日 10.32

期望信号\(u(t)\)的测量\(\tilde{u}(t)\)被确定性误差\(d(\vec{u}(t))\)和累积随机误差\(z(t)\)破坏。确定性误差是具有未知确定性系数的分析模型足够的传感器缺陷。其中一些系数在仪器的生命周期内仅表现出微小的变化。

这些可以在工厂校准过程中估计和补偿。其他确定性错误,如启动偏差,可以通过状态增广来实时估计。确定性错误的形式可能包括比例因子误差、非线性、对陀螺的 g 敏感性、非正交轴和交叉耦合的轴[5]。请参阅“问题和权衡讨论”中的“确定性错误”小节。

本文的重点是由\(z(t)\)表示的随机误差,这些误差可能来自各种物理现象(参见“背景”部分)。每次仪器打开时,随机误差是不同的,随时间变化的函数,不能根据传感器测量\(\tilde{u}(t)\)进行预测。

为清楚起见和简单性,本教程的大部分将\(z(t)\)视为标量信号。基本思想适用于六自由度IMU随机误差模型。

State Estimation Error Model

当导航系统通过积分式(2)的非线性模型来传播车辆状态向量时,误差表示为\(\delta \vec{x}_v(t)=\vec{x}_v(t)-\hat{\vec{x}}_v(t)\),无论选择哪种姿态表示法(例如方向余弦矩阵或四元数),姿态误差都可以由具有三个分量的向量来表示。因此,车辆误差状态包含位置、速度和姿态误差的子矢量,每个是具有三个分量的矢量,即\(\delta \vec{x}_v(t) \in \Re^{n_v}, n_v=9\)

车辆状态误差向量可以使用辅助传感器[10]、[13]、[12]、[31]、[43]、[44]的测量结果实时估计。估计算法包含误差状态的线性化状态空间模型:

Clipboard 2023年3月11日 10.46

其中,

  • \(F(t)=\left.\frac{\partial f(\vec{x}, \vec{u})}{\partial \vec{x}}\right|_{\hat{\vec{x}}(t), \hat{\vec{u}}(t)}\)
  • \(G(t)=\left.\frac{\partial f(\vec{x}, \vec{u})}{\partial \vec{u}}\right|_{\hat{\vec{x}}(t), \hat{\vec{u}}(t)}\)
  • \(\delta \vec{u}(t)=\vec{u}(t)-\hat{\vec{u}}(t)\)
  • 参见“Simplified INS Example.”中的示例F(t)和G(t)矩阵

这个状态空间模型并不完整,除非指定\(\delta \vec{u}(t)\)的IMU误差模型。对于状态估计,确定状态空间形式的IMU误差模型。

实时状态估计过程用于估计增广的状态向量:

Clipboard 2023年3月11日 10.49

由车辆误差状态向量\(\delta \vec{x}_v(t) \in \Re^{n_v}\)组成,其中,

  • \(\vec{x}_d(t) \in \Re^{n_d}\)表示imu确定性误差
  • \(\vec{x}_z(t) \in \Re^{n_z}\)表示imu随机误差
  • 错误状态的总维度\(n_x=n_v+n_d+n_2\)

状态增广的过程,在“State Augmentation”中进行了讨论。这里的介绍将完全集中在随机IMU误差上\(z(t)\),其在式(3)中的累积用\(\vec{z}(t)\)表示。

标量\(z(t)\)的状态空间模型:

Clipboard 2023年3月11日 10.53

其中,

  • \(A_z \in \Re^{n_z \times n_z}, B_z \in \Re^{n_z \times p}\), and \(C_z \in \Re^{1 \times n_z}\)
  • 参数p表示IMU误差模型的微分方程式部分中不同的独立的噪声过程的数目
  • 参数\(n_z\)表示IMU随机误差模型中的状态数

随机信号\(\vec{\omega}_z(t)\) and \(\eta_z(t)\)式相互独立的具有功率谱密度\((S_{\omega_z} \in \Re^{p \times p}\) and \(S_{\eta_z} \in \Re)\)的高斯白噪声过程。假设\(\vec{\omega}_z(t)\)的元素是相互独立的,因此得到\(S_{\omega_z}\)是对角线的。

设计者在选择模型结构(特别是\(n_z\)\(P\))时必须谨慎,因为整个模型将有6\(n_z\)状态和6\(P\)个独立噪声源。

线性状态空间系统的功率谱密度

对应于(6)-(7)中的状态空间模型,频域模型为

Clipboard 2023年3月11日 11.12

其中,S是拉普拉斯变量。从\(\omega_z(t)\)\(z(t)\)的传递函数模型为:

Clipboard 2023年3月11日 11.13

它有1行和p列,符号\(Z(s), \Omega_z(s)\), and \(\eta_z(s)\)分别表示为信号\(z(t), \omega_z(t)\), and \(\eta_z(t)\)的拉普拉斯变换,因此,对应于信号\(z(t)\)的功率谱密度(PSD):

Clipboard 2023年3月11日 11.15

假设驱动噪声向量\(\omega_z(t)\)和输出噪声\(\eta_z(t)\)的所有元素是相互独立的并且是白的,这简化为

Clipboard 2023年3月11日 11.16

其中,

Clipboard 2023年3月11日 11.17

\(T_i(s)\)\(\vec{\omega}_z(t)\)的第i个元素到\(z(t)\)的(标量)传递函数,\(B_{z_i} \in \Re^{n_z \times 1}\)\(B_z\)的第i列。每个\(T_i(s)\)\(s\)中分子与分母之比,\(S_z\)\(\omega\)的正实函数。因此,每个\(T_i(j \omega) T_i(-j \omega)\)\(S_z(\omega)\)将始终是拉普拉斯变量\(s\)的偶次幂的多项式函数的比率。证明这一事实的例子在《Finite Dimensional Linear StateSpace Systems have Even Power Spectra.》中给出

\(S_z\)是正实函数的事实导致了IMU误差建模方法中的主要挑战之一,因为一些IMU随机误差分量的PSD不能用式(10)的总和中的项精确地拟合。因此,设计者必须在近似状态空间模型中做出明智的选择,以达到满意的折衷。这个问题将在"Modeling via Independent Noise Sources”一节中进一步讨论。

Available Error Specification Information

为了表征IMU的质量(根据IEEE规范[28]、[29]、[30]),制造商提供了AV曲线图、Allan标准偏差(ASD)曲线图或从中提取的参数。ASD图表既可以帮助仪器设计者理解和改进他们的传感器,也可以将预期性能传达给预期的用户。在本文的上下文中,主要主题是INS设计人员如何使用来自ASD图的信息来指定式(6)-(7)中的状态空间模型的参数。

图1和图2显示了示例ASD曲线图。图1显示了Crossbow μNav IMU陀螺仪的ASD曲线图[45]。

Clipboard 2023年3月11日 11.28

这些ASD是根据[36]的作者提供的数据计算得出的。蓝色、黑色和绿色星号(‘*’)标记ASD数据点。

图2中的ASD10绘图是根据制造商提供的数据计算得出的,该数据来自安装在隔振系统顶部的大大理石平板上的IMU。图2中的每个蓝色‘x’标记一个ASD数据点。每个ASD图的水平轴是集群时间(或大小),符号τ以秒为单位。

Clipboard 2023年3月11日 11.28

请注意,这两张图中的ASD图既有相似之处,也有不同之处。对于较小的集群大小,两者都以-1/2的特征斜率减小,如每个图中的红色虚线切线所示。然后,在虚线青色切线所指示的值处,两者都变平到零的斜率。

对于较大的集群大小,图2中的ASD以1/2的斜率增加,如用黑色虚线绘制的切线所示。对于像τ=1000这样大的集群时间,图1中的ASD曲线图并没有(强烈地)表现出斜率为1/2的这种增加。

这些τ的值和发生变化的ASD值对于每个仪器都是不同的。它们指定了某些参数,这些参数对于比较惯性仪器的性能和评估权衡构建IMU随机误差模型都很有用。

问题描述

本文的目的是讨论与开发模型相关的方法、问题和权衡,该模型将IMU误差的随机部分的性质定量地传达给状态估计算法。该误差模型的输入将是随机信号。

制造商提供的可用于模型开发的信息是AV或ASD特性。设计者无法获得产生这些特性的IMU真值输出数据;因此,系统辨识方法不适用。

由于误差状态估计算法是状态空间形式的,所以IMU误差模型也是状态空间形式的。输入被建模为独立的高斯白噪声过程。挑战在于构建这些模型,使它们具有与IMU相同的输出统计特性(即,AV方差)。

请注意,这种随机状态空间模型并不是唯一的。事实上,本文并不打算提出一个特定的模型;尽管本文讨论了一个特定的ASD情节,并给出了一个模型作为教程示例。相反,本文的目标是清楚地展示各种文章、书籍[9]和标准16[30]、[28]、[29]中使用或暗示的方法(使用各种特定模型)。作者认为,这种方法的许多方面都是行业标准的。不幸的是,大多数描述这种方法的出版物都不公开。

从高层次上,该方法的概要如下:

  • (1)利用来自AV/ASD图的信息构建连续时间状态空间模型。
  • (2)将连续时间为IMU的随机模型转化为等价的离散时间模型。
  • (3)在模拟中使用离散时间模型产生数据,计算ASD图与仪器的ASD图进行比较。
  • (4)当设计者对IMU误差模型满意时,将其附加到车辆状态误差模型中,并用于INS误差状态估值器的设计。

Allan Variance

AV是一种众所周知的时域分析技术,最初开发该技术是为了表征和研究振荡器的频率稳定性[48]、[49]、[50]。由于其相对简单,已被成功地用于传递IMU性能规范和来表征它们的随机误差[28]、[29]、[32]、[34]、[36]、[39]、[43]、[44]、[46]、[47]

给定一组数据,计算AV的过程如下。设\(D=\left\{\tilde{u}_i\right\}_{i=1}^L\)是一组(去趋势的)比力(或角速率)数据(由静止的IMU以恒定采样间隔T测量)。对于每个\(n \in[1, L / 2]\),以\(\tau=n T\)为变量,以\(\tau\)的范围为从T到LT/2的簇时间的值分别计算方差。

对于给定的n,在每个时刻\(t_i \in[T, 2 T, \ldots,(L-n) T]\),一组n个连续的数据点(从\(t_i\)开始)形成一个数据集\(\left\{\tilde{u}_j\right\}_{j=i}^{i+n-1}\),为每个这样的n点聚类计算平均值\(\bar{u}_i(\tau)=\frac{1}{n} \sum_{j=0}^{n-1} \tilde{u}_{i+j}\)。然后,以\((L-2 n)\)平方簇差[28]、[49]的平均值来计算持续时间τ的AV:

Clipboard 2023年3月11日 16.37

由于一些IMU(尤其是那些高级IMU)提供比力(或角速率)的积分,表示为\(\tilde{\theta}_i\),因此\(\bar{u}_i(\tau)\)可替换为:

Clipboard 2023年3月11日 16.38

将式(13)带入式(12)得到

Clipboard 2023年3月11日 16.40

这是计算AV的另一种公式[28]。

对于图形分析,AV的平方根\(\hat{\sigma}_u(\tau)\),称为ASD,被称为ASD,通常以对数-对数比例绘制,由于数据集D的长度有限,每组数据簇的数据量会随着\(\tau\)的增加而减少;因此,ASD的标准差(注意,不是ASD本身)计算为:

Clipboard 2023年3月11日 16.48

其中,\(\kappa\)是经验常数,对于imu分析,通常取\(\kappa \approx 1 / \sqrt{2}\),[53], [28], [29], [34], [54].

Power Spectral Density and Allan Variance

AV通过以下方式与双边PSD相关:

Clipboard 2023年3月11日 16.56

IEEE标准952-1997[28]中(C.1)后面的文字将此等式解释为,当通过传递函数时,Allan Variance与信号\(u\)中的总噪声功率成比例,该传递函数由用于创建和操作簇的方法确定。式(16)的推导可以在[55]第79页找到,(16)没有求逆公式(见[53])。在这个表达式中,\(S_u(f)=\left.S_u(s)\right|_{s=j 2 \pi f}\),其中\(s \in \mathbb{C}\)是拉普拉斯变量,\(j=\sqrt{-1}\)是复数、\(f \in \Re\)的单位是hz。

Modeling via Independent Noise Sources

当功率谱表示为频率f的幂函数级数时,它的形式为:

Clipboard 2023年3月11日 17.02

这种形式的PSD很方便。根据叠加原理,它与信号的功率谱相对应:

Clipboard 2023年3月11日 17.02

其中,信号\(z_N(t), z_B(t)\), and \(z_K(t)\)是相互独立的零均值噪声过程。

在此假设下,对式(16)到(17)应用,将产生具有以下形式的AV:

Clipboard 2023年3月11日 17.04

其中每个AV项的具体功能形式可以计算,并可在各种来源中获得[28]、[29]、[30]、[34]、[55]。每个AV项的函数形式很容易与ASD图的一部分相关联。

“Continuous-Time State-Space Models”一节描述了建立这种关联并为每个术语定义连续时间状态空间模型的方法。每项的状态空间模型由它自己的独立的高斯白驱动噪声驱动,从而导致每个信号\(z_N(t), z_B(t)\), and \(z_K(t)\)相互独立。

  • 这些状态空间模型对于对应于式(17)中f的偶函数的项(如\(\frac{K^2}{(2 \pi f)^2}\))可以是精确的
  • 然而(如前所述并在侧栏中举例说明)是频率f的奇函数的功率谱项(如\(\frac{B^2}{2 \pi f}\))不能由任何有限维、线性、状态空间模型精确建模,因此,必须对这些术语进行近似建模,仔细权衡。

在式(17)的幂级数表示中可以包括任意数量的项。这导致式(18)的信号模型和式(19)的AV模型中的项的数量相同。每个项代表来自独立来源的不同类型的噪声。ASD图的典型形状如图3所示,有五个独立的噪声源(另见[28]中的图C.8)

Clipboard 2023年3月11日 17.13

在ASD图中,每种噪声类型都与一个特征斜率相关联,该特征斜率有助于识别该噪声类型及其模型参数。并不是所有的噪声类型在每个设备都很明显。当存在时,噪声项占主导地位的\(\tau\)的模型参数和范围对于每个仪器可能是不同的。

在商业级IMU中,N、B和K项通常占主导地位(例如,参见图1和图2)。对于较小的\(\tau\),仪器设计的选择(例如量化方法和10个采样周期)会导致随机误差表现为白噪声。这种白噪声在随机游走噪声项中被考虑(也就是N)。然而,随机误差并不是完全真正的白噪声。随着集群时间\(\tau\)的增加,ASD曲线图可能表现出偏置不稳定(B)、速率随机游动(K)和其他噪声类型。要使ASD曲线图显示这些其他噪声类型,用于生成ASD曲线图的IMU数据集必须非常长。

ASD图绘制到几分钟(例如,几百秒),对于在观测不可用的间隔期间分析性能是有意义的。然而,对于非常大的\(\tau\),ASD曲线的具体形状通常是不确定的,并不重要。

表1的第1列和第2列包括了陀螺仪和加速度计的 N、B 和 K 个噪声项的具体名称。第3列和第4列总结了这些噪声类型的AV和PSD之间的关系,[28]、[34]

Clipboard 2023年3月11日 17.25

其中,N、B 和 K 项将是“"Continuous-Time State-Space Models”部分讨论的焦点。

Continuous-Time State-Space Models

本节考虑开发近似再现ASD图和式(17)PSD的连续时间状态空间模型。使用图2中的示例ASD来说明整个思路。图1将仅在参考B和K幂级数项时才讨论。整体模型将具有以下形式:

Clipboard 2023年3月11日 17.30

其中, \(z_N(t), z_B(t)\), and \(z_K(t)\)分别是与系数N、B和K相关的IMU随机误差信号。

随机游走误差\(z_N(t)\):角速度和速度

式(17)中的PSD项\(N^2\)相对于频率f是恒定的,这对应于白噪声的功率谱[56]。因此,

Clipboard 2023年3月11日 17.32

其中,\(\omega_N(t)\)是PSD为\(S_N=N^2\)的高斯白噪声。

在文献中和制造商规范中,这种类型的误差称为陀螺角速度随机游走误差加速度计的速度随机游走误差

将式(16)中的变换应用于\(S_{z_N}(f)=N^2\),将得到[28],如表1中所示:

Clipboard 2023年3月11日 17.34

这表明,在ASD图上,角/速度随机游走将用斜率为\(-\frac{1}{2}\)的直线表示,如图4所示。

Clipboard 2023年3月11日 17.36

随机游走参数N的值可以从提供ASD绘图的制造商中近似确定,这是通过识别ASD图上斜率为\(-\frac{1}{2}\)\(\tau\)范围,并绘制其切线来得到的。

在图2中,在\(\tau \in[0.01,30]\)范围内绘制了红色虚线切线。从式(23)可以看出\(\left.\sigma_{z_N}(\tau)\right|_{\tau=1}=N\),因此,可以从ASD图中(\(\tau=1 \mathrm{~s}\),斜率为\(-\frac{1}{2}\)的切线)提取N的值,对于图2,结果是\(N \approx 0.0033 \mathrm{~m} /\mathrm{s}/ \mathrm{s}^{1 / 2}\)

随机游走误差\(z_K(t)\):速率和加速度

式(17)中的项\(\frac{K^2}{(2 \pi f)^2}=\left.\frac{K}{s} \frac{K}{s^*}\right|_{s=j 2 \pi f}\)对应于具有传递函数为\(T(s)=\frac{1}{s}\)的线性系统,因此,其状态空间模型为:

Clipboard 2023年3月11日 17.49

其中\(z_K(t)\)是输出,输入\(\omega_K(t)\)是带有PSD为\(S_K=K^2\)的高斯白噪声。

在文献中和制造商规范中,这种类型的误差称为陀螺的速率随机游走误差和加速度计的加速随机游走误差.

给定式(24)和(25),\(z_K(t)\)的PSD是:

Clipboard 2023年3月11日 17.57

其所需形式对应于式(17)中的第三个项。使用式(16)代换\(S_{z_K}(f)\),得到:

Clipboard 2023年3月11日 17.58

总结在表1的相应行中,等式 (27) 显示,在ASD图上,速率/加速随机游走误差将由斜率为\(+\frac{1}{2}\)的线表示,如图5所示。

Clipboard 2023年3月11日 18.00

速率/加速度随机游走参数K可以从制造商提供的ASD图中近似确定。

第一步是识别ASD图上斜率为\(+\frac{1}{2}\)\(\tau\)范围(如果存在),并绘制与其切线。在图2中,绘制的虚线黑线为\(\tau \in[3,500]\)秒,斜率为\(+\frac{1}{2}\),在\(\tau \geq 100\)秒时与ASD曲线大致相切。因为\(\tau\)很大,这部分ASD通常具有更高程度的不确定性[如相对于 式(15) 所讨论的那样]。

第二步使用切线来估计K。从式(27)可以看出\(\left.\sigma_{z_K}(\tau)\right|_{\tau=3}=K\),因此,从ASD图中估计K值的简单方法是在\(\tau = 3 s\)时找到直线近似的值,在图2的例子中,\(K \approx 0.00014 \mathrm{~m} / \mathrm{s}^2 / \mathrm{s}^{1 / 2}\)

基于图1中的ASD图,μNav单元中的陀螺可能不需要包含角速率随机游走噪声,即使聚类时间可达1000秒。

累积误差模型:N,K

由于角随机游走和速率随机游走误差(或速度随机游走和加速度随机游走误差)都具有功率谱,因此可以直接建立状态空间模型来再现功率谱中相应的项,以及它们的ASD图部分。

基于前面的两个部分,状态空间模型将是:

Clipboard 2023年3月11日 18.09

其中,\(z_N(t)=\omega_N(t)\)是PSD为\(N^2\)的白噪声,\(\omega_K(t)\)是PSD为\(K^2\)的白噪声。

随机信号\(\omega_N(t)\)\(\omega_K(t)\)是独立的,因此,\(z_N(t)\)\(z_K(t)\)也是独立的。

该模型的ASD为:

Clipboard 2023年3月11日 18.12

该模型的ASD图使用值\(N \approx 0.0033 \mathrm{~m} /\mathrm{s}/ \mathrm{s}^{1 / 2}\)\(K \approx 0.00014 \mathrm{~m} / \mathrm{s}^2 / \mathrm{s}^{1 / 2}\),如图2所示的绿色实线。

偏差不稳定性\(z_B(t)\)

一些ASD图(如图1中的图)没有表现出与\(z_K(t)\)相关的\(+\frac{1}{2}\)斜率,但对于较大的\(\tau\)值,确实有一个较宽的平坦区域。该平面区域不能由N或K值项来很好地建模。在这种情况下,式(30)中的NK模型的AV(和 ASD)图中,\(\tau\)在N和K这个中间范围内太小。在任何一种情况下,都有足够的偏差不稳定性,这样可以通过在模型中考虑其性能来改进。

\(S_{z_B}(f)=\frac{B^2}{2 \pi f}\)对应的误差项\(z_B(t)\)通常被称为偏差不稳定性(或闪烁噪声)[28], [29], [34], [44], [39].

将 式(16) 应用于\(S_{z_B}(f)=\frac{B^2}{2 \pi f}\) for \(f \leq f_0\)(并且排除0的情况),将会得到:

Clipboard 2023年3月11日 23.09

其中,

  • \(x=\pi f_0 \tau\)
  • \(C i\)是[56]的cosine积分函数
  • 参数\(f_0\)是截断频率[28]

偏差不稳定性ASD图如图6所示。

Clipboard 2023年3月12日 09.42

该图显示,对于小\(\tau\)\(\sigma_{z_B}(\tau)\)增长,直到\(\tau>\frac{1}{f_0}\)的平台阶段。因此,\(\tau \approx \frac{1}{f_0}\)的值定义了偏差不稳定性(或闪烁噪声)对其最大值贡献的部分。在该区域内,可以证明式(31)中的正弦项和余弦项接近于零,因此在平坦区域,有:

Clipboard 2023年3月12日 09.47

这些方程提供了一种从ASD图中提取B近似值的简单方法。在这种方法中(例如,参见[35]中的p.6,[44]中的p.21, [54]中的p.10和[55] 中的p.114),可以从[28]中的B.4.5 节中推断,ASD图平坦的部分所对应的\(\sigma_z^2(\tau)\)\(\frac{2 B^2 \ln (2)}{\pi}\)近似,来推断参数B的值。

对于图1中的ASD图,青色水平线近似于9.5e-3和1.40e-2deg/s的最小ASD值,对应的B的值为1.43e-2和2.11e-2 deg/s (9.5e-3/0.664, 1.40e-2/0.664),对于图2中的ASD图,\(7.4 e-4 \mathrm{~m} / \mathrm{s}^2\)的最小ASD值对应于\(B=1.11 e-3 \mathrm{~m} / \mathrm{s}^2\)

因为偏差不稳定性项的功率谱项(即\(\frac{B^2}{2 \pi f}\))不是\(s=j 2 \pi f\)的偶数幂,所以没有完全拟合它的有限阶线性状态空间模型。因此,导航系统设计者必须选择一个状态空间模型来逼近偏差不稳定性误差效应。

人们提出了各种方法来近似地解释偏置不稳定性。其中包括一阶高斯-马尔科夫[4],[32],[36],[39],[44],[58]和高阶自回归模型[43],[46],[47]。一个重要的权衡是,随着状态空间模型维数的增加,近似值的保真度可能会增加,但状态估计算法所需的实时计算负载也会增加。此外,更精细的模型可能对未建模的动力学和非线性不鲁棒,特别是当一些添加的状态是弱可观的。这些主题将在“Discussion of Issues and Tradeoffs”一节中进一步分析。

为了举例说明这个想法,下一节将考虑一阶高斯-马尔可夫模型,该模型使用指数相关噪声来模拟偏差不稳定性误差。

Gauss-Markov Error Model

一阶连续时间高斯-马尔可夫模型为[9],[56]:

Clipboard 2023年3月12日 10.02

且:

Clipboard 2023年3月12日 10.03

其中,

  • 符号\(T_B\)表示进程的相关时间。
  • 符号\(\omega_B(t)\)表示带有PSD为\(S_B\)的白色驱动噪声

式(33)对应的传递函数为:

Clipboard 2023年3月12日 10.04

产生的PSD为:

Clipboard 2023年3月12日 10.05

将式(16)应用于\(S_{z_G}(s)\)将得到:

Clipboard 2023年3月12日 10.07
Clipboard 2023年3月12日 10.33

\(\sigma_{z_G}(\tau)\)的曲线如图7所示。一些特殊情况值得注意:

(1) 对于较小的集群时间(\(\tau<<T_B\)

Clipboard 2023年3月12日 10.10

使ASD图对于小的\(\tau\)阶段的斜率为\(+\frac{1}{2}\)

(2) 当\(\tau=1.89 T_B\),曲线平滑且:

Clipboard 2023年3月12日 10.12

\(\sigma_{z_G}\left(1.89 T_B\right)=0.4365 \sqrt{S_B T_B}\)

此处公式要注意,\(\sigma_{z_G}\left(1.89 T_B\right)\)是一个值,不是两个值相乘,所以\(\sigma_{z_G}^2\left(1.89 T_B\right)\)应该看作\(\sigma_{z_G}^2(\tau)\)

(1) 对于较大的集群时间(\(\tau>>T_B\)

Clipboard 2023年3月12日 10.14

使ASD图对于大的\(\tau\)阶段的斜率为\(-\frac{1}{2}\)

一阶标量高斯-马尔可夫过程可用于(近似地)对ASD图的平坦部分(即偏置不稳定性)建模。

为什么长这样,与直觉相反?因为为了产生平坦区域,前面部分斜率为\(+\frac{1}{2}\)刚好可以跟角度/速度随机游走N的斜率\(-\frac{1}{2}\)抵消。后面部分也同样道理,但是这样就需要有截断的频率。

如果制造商只提供了\(B\)\(T_B\)的值,那么\(\mu_B\)的值可以用(34)计算。通过使式(32)中的\(\sigma_{z_B}^2\)等于式(37)中的\(\sigma_{z_G}^2\),可以求解\(S_B\)

Clipboard 2023年3月12日 10.18

\(\mu_B\)\(S_B\)已知时,式(33)的状态空间模型完全指定。

相反,如果制造商提供ASD图,并且偏差不稳定性足够大,足以保证包含在模型中,那么可以首先根据\(1.89T_B\)靠近ASD图的平坦部分来得到\(T_B\)参数,然后使式(37)中定义的\(\sigma_{z_G}\left(1.89 T_B\right)\)的值接近ASD图平坦区域的值,从而得到\(S_B\)

对于图2中的ASD图,在\(\tau = 60 s\)时,最小值为\(7.4 e-4 \mathrm{~m} / \mathrm{s}^2\)。对应\(T_B= 60/1.89 = 31.7\), \(S_B=9.0 e-8 \mathrm{~m}^2 / \mathrm{s}^5\), \(B=1.11 e-3 \mathrm{~m} / \mathrm{s}^2\)

Cumulative Error Model: N , B, K