扩展卡尔曼滤波器与非线性演化的不稳定空间数据同化的C软件包实现

*stec=SoftwareX 7（2018）28原始软件出版物EKF-AUS-NL算法在没有线性正切模型和存在参数模型误差的情况下实现Luigi Palatellaa，*，Fabio GrassobaLiceo Scientific Statale ''C. Deb大气和气候科学研究所ISAC-CNR，Str. 箴Lecce-Colononi km 1，200，I-73100，Lecce（莱切），意大利ar t i cl e i nf o文章历史记录：接收8九月2017收到修订版2017年12月14日接受2018年关键词：卡尔曼滤波AUSLorenz96SLAMa b st ra ct在本文中，我们提出了一个C++-软件包实现算法EKF-AUS-NL（扩展卡尔曼滤波器与非线性演化的不稳定空间中的同化），旨在进行数据同化的不稳定空间时的微分方程的雅可比矩阵无法计算。我们还提出了一个简单的方法来考虑存在的模型误差的框架EKF-AUS-NL。该软件使用EKF-AUS-NL算法进行数据同化，并将动力系统定义为单独实现的通用时间演化例程我们提出了两个说明性的例子的基础上Lorenz96和SLAM系统。版权所有©2018作者.由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）中找到。代码元数据当前代码版本v1.0用于此代码版本的代码/存储库的永久链接https://github.com/ElsevierSoftwareX/SOFTX-D-17-00068法律代码许可证Mozilla公共许可证，2.0（MPL-2.0）https://mozilla.org/MPL/2.0/使用的代码版本控制系统无使用C++的软件代码语言、工具和服务编译要求，操作环境gcc 4.2.1，Eigen软件包http://eigen.tuxfamily.org[16]，MacOSX，Linux开发者手册链接https://github.com/ElsevierSoftwareX/SOFTX-D-17-00068问题支持电子邮件luigi. yahoo.it1. 动机和意义在最近的一些工作[1-这种方法也导致了一个相当自然的非线性扩展的算法，由于计算的自我和交叉非线性之间的相互作用的领先李雅普诺夫向量[3]。先前的工作处理系统提供的切线线性模型计算感谢的雅可比矩阵的解析表达式的确定性微分方程xi Fi（x）/xj。在本文中，我们提出了一个变化的算法EKF-AUS-NL设计执行数据同化过程时，通讯作者。电子邮件地址：luigi. yahoo.it，luigi. liceodegiorgi.gov.it（L. Palatella）。https://doi.org/10.1016/j.softx.2018.01.001F i（x）/此外，我们提出了一个简单的方法，以考虑到存在的参数模型误差的框架中的EKF-AUS-NL例程。我们提出了一个C++软件包，执行数据同化使用EKF-AUS-NL算法的动力系统定义为一个通用的时间演化例程。该软件包可以执行完美的模型情景测试（处理受模型误差影响的参数作为额外的状态变量）和实际的数据同化测试，使用一个程序设计，以获得实际的测量。时间演化与EKF-AUS-NL例程分开实现，因此它可以用任何语言编写，而动态系统的作为说明性的例子，我们展示了两个完美的模型测试应用于两个模型：Lorenz96模型[4]和同时定位和地图模型（SLAM）[52352-7110/©2018作者。由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）。可在ScienceDirect上获得目录列表SoftwareX期刊主页：www.elsevier.com/locate/softxL. Palatella，F.Grasso/SoftwareX 7（2018）2829DT=+++=FJJJr=1q=1在本文和代码中，我们使用统一符号LLJJK我LL我我a是新的（正交的）扰动。Lr=1q=1tJJ软件包是用C++面向对象语言实现的。为了保持软件体系结构和动 力 系 统 代 码 之 间 的 分 离 ， 所 提 供 的 包 包 括 两 个 例 程（PrepareForEvolution和PrepareForAnal- ysis），其从分析和预测误差扰动Xa，f开始，读取/写入由时间演化例程（其使用系统C++命令运行）使用的初始条件文件。分析/预报误差扰动Xa，f是适当数量的列向量，表示在同化过程之后和之前的误差的最佳估计。这个软件是，在我们看来，非常有价值的，因为它允许EKF-AUS依赖最近已被证明为线性系统的分析[12]。为了处理线性近似不完全成立的情况，我们使用[3]中描述的非线性1.1.1. EKF-AUS的非线性扩展：EKF-AUS-NL算法[3]的作者在他们的附录A中得到，如果考虑到直到二阶的非线性，则扰动，即X的列，从t=tk时的Xa演化到tk+1 =tk时的Xf。tk+τ（其中τ是分析区间），根据微分方程（爱因斯坦关于重复指数的约定，Fi;j 我的朋友，对于非常复杂的模式，使用EKF-AUS-NL技术进行数据同化，只需编写一个接口程序Fi;jk Fi）运行动力系统集合成员。此外，管理参数模型误差的可能性允许人们使用我们的算法，即使在某些模型参数受不确定性影响的情况这种情况发生，例如，对于SLAM算法，dX=Fi;jXjs，s≤mdX是（qr）=FijXjs（qr）+1α<$Fi;jkXjqXkr，（四）Rithm在下面介绍在这里，我们简要回顾EKF-AUS-NL算法，以便dt，2介绍正确的符号，并解释我们如何处理其中q≤r;q，r∈ [1，m l];s（q，r）= m+∑ml ∑q≤r1，其中M ≤模型误差lm是涉及到的线矢量的数量，非线性相互作用，α<$=3。人们当然可以考虑所有1.1. 理论背景EKF-AUS 滤 波 器 是 [9] 中 正 式 介 绍 的 扩 展 卡 尔 曼 滤 波 器（EKF）的一种特殊平方根实现[8]EKF-AUS算法是通过将近似限制在一个m维流形上而得到的。当m等于系统的自由度的数目N标准EKF方程。当mN+N0的简化形式，本文给出了非稳定子空间中的带同化（EKF-AUS），其中N+和N0分别为正和中性李雅普诺夫指数的个数注意，对于一般的动力系统，关系式N++N0<$N成立。m个扰动向量可能的非线性相互作用这意味着该算法应该涉及m（m1）/2更多向量，即所有可能的对获得的数量 m线性向量。在一些情况下，由于要考虑的矢量数量太大，无法采用这种方法。出于这个原因，我们限制非线性相互作用的领先（最不稳定）ml向量。关于索引s， ml， q， r的条件只不过是一种惯例，以保持向量根据与第一个m分离的非线性项进行EKF-AUS算法中已涉及N+N0个向量。通知在分析过程中涉及的向量的总数现在变为m+m（m+1）/2。[10]中：系统的总自由度为N，测度数为p。分析误差协方差矩阵表示为矩阵乘积Pa =XaXT，其中X为以下之一前面几段所述的背景是基于在[3]上。然而，我们需要理解如何处理非线性，因为它经常发生，微分方程xstec=the ‘‘square root’’ of一. 在预报期间，驱动动力系统的F（x）是已知的，但存在难以定义或计算一阶和二阶导数步骤，线性化的时间演化作用于Pa的列（扰动）作为Pf=MXa（MXa）T=XfXT.（一）在这里提出的算法中，时间演化是用整个动力系统来执行的，该动力系统生成轨迹xa的系综，每个轨迹由下式给出：f（x）。作为第一步，我们必须找到一个典型的尺度η（在相空间中，而不是在模型的物理空间中），在该尺度下，两个轨迹之间的差异的演化基本上是线性的。有不同的技术来估计η，我们建议读者参考文献，例如[13，14]。这里是[3]的附录B中建议的处理之间的轻微差异为了更好地再现线性化动态，xa（tk）=xa（tk）+δxa（二）线性切线模型是不可用的，我们意识到，其中δxa是矩阵Xa的第j列，xa（tk）是第k个分析时间tk的分析，即系统状态的最佳估计。经过时间演化，我们得到了最优状态xf（tk+1）的估计（预测）和Xf矩阵构建形式上正确的方法。与标准情况一样，我们有一个定义分析xa（tk）的轨迹。驱动从分析时间tk到预测时间的轨迹t k+1 =t k+τ是M（xa）。因此我们有xf（t k+τ）= M（xa（tk））。向上收集列向量δxf，δxf=xf（tk+1）−xf（tk+1），（3）给定在时刻t的m+m（m +1）/2个扰动Xa，让我们将Xa的第j列定义为δXa。我们得到第一个m+m l（m l+1）/2其中，xf（tk+1）是定义在中的状态xa（tk）的时间演化，xa（t）xa（t）Xa J1m m（m（1）第2项（五）ijfjK=k+η δj，∈[，+ll+/]当量（一）. 然后，我们对扰动δxi进行正交正规化，即Xf的列，然后我们遵循标准EKF-AUS方法则剩余的ml（ml+1）/2个轨迹被获得为如[9]中所述作为常规同化的最终结果，一a1a a得到分析状态xa（tk+1）和分析扰动由矩阵Xa给出 = [δxa，δxa，. . .，δxa]。矢量δxa，xs（tk）=x（tk）+2[δxr+δxq]，（6）X的列1 2m其中r，q∈ [1，m]，q≤r;s（q，r）=m+m（mLL+1）/2+当观测足够密集时，这些结果成立[11]这是一个很好的解释，也是一个很好的解释。∑ml ∑q≤r 1。在这一点上，所有的轨迹都进化到EKF正常工作而不受时间tk+1=k到分歧事件。收敛假设xf（tk+1）=M（xa（tk）），j∈[1，m+ ml（ml+1）].（七）在附录B中简单地提出并不完全正确。我们现在表明以这种方式进化的轨迹+τ使用全非线性系统我30L. Palatella，F.Grasso/SoftwareX 7（2018）28线性扰动模型误差非线性扰动而在进化之后，线性化和+2. 软件描述非线性非线性扰动模型+⎝δx一、二、∑∑+===J+==+=+α<$+=RQ一一⎠≤≤注意，在[3]的线性化处理中，ml（ml1）/2的非线性扰动的数量对应于ml（ml1）个新的系综成员。右边的Eq。（十二）、⎛0⎞为了在给定的时间演化之后获得扰动Xf，通过动力学模型，我们如下进行对于第一个M向量δxf（tk+1）=1[xf（tk+1）−xf（tk+1）]，k∈[1，m]（8）Xa=X二，一2，m⎟jηj↓M而其余的列则用来给出对应于方程2中α′项的非线性强迫。（四）：1η一一，一2002年2月1日XF =δx. . .δxam. . .δxam三、=0（十二）δxf（tk+1）= [xf（tk+1）−xf（tk+1）]. . .. . .. . .⎜ ⎟S s{1}下一页（九）尼日利亚，1. . .δxamM Lq ≤rr，q∈ [1，ml]，q≤r;s（q，r）=m+1，（10）r=1q= 1其中，δxf（τ）、δxf（τ）由等式给出（8）为r， q ml m.因此，最后一个方程的意义是简单地迫使最后m l（m l1）/2个非线性扰动的演化，其中线性和非线性时间演化之间的差乘以因子α<$。1.1.2. EKF-AUS-NL算法中参数模型误差的处理本文提出了一种EKF-AUS-NL算法中参数模型误差的处理方法众所周知，考虑模型误差的一种直接方法是将受这种误差影响的参数作为附加状态变量进行管理（例如，参见[15]）。处理模型误差的一个典型假设是认为它具有无偏差的高斯分布。此外，通常假设模式误差在连续同化时间tk之间没有相关性。为了更好地解释这里使用的方法，我们展示了SLAM示例的具体应用。在这种情况下，两个参数，即速度V和转向角γ，被“机械地"确定，但是分别受到具有标准偏差σ V和σ γ的高斯误差的V和γ被添加到系统的状态变量中，并且在每个时间演化间隔内被假定为恒定的，而它们可以根据车辆试图遵循的路径（在我们的示例中为正圆）在不同的演化间隔之间变化就在第k次分析之后的时间演化之前，在时间tk处，我们得到附加状态变量V和γ上的误差与分析误差的其他分量因此，要演化的不同系综成员从如在等式的左侧中布置的扰动矩阵Xa开始获得（十二）、这种选择是合理的，因为V和γ相对于其他状态变量的主要区别在于，就在时间演化之前，对应的扰动具有除对应于受模型误差影响的变量之外的全部零分量，即δxa=δij（ 11）其中指数j仅指作为参数模型误差处理的扰动（在这种情况下为V和γ）。在演化之后，误差显然与其他变量纠缠在一起，并被用于分析，但就在这一步之后，扰动被重置为方程的状态（十一）、因此，在时间tk+1处的分析是利用如图10中所示的扰动矩阵Xf来为了结束这一理论部分，显示一个明确的例子是有用的。为了解释不同的扰动排列在矩阵Xa或Xf中。假设，正和零李雅普诺夫指数m并且存 在受模型误差影响的两个参数（如在下面示出的SLAM示例中）。用户可以使用函数SetMod- elErrorVariable（long istart，long iend，MatrixXd error，MatrixXd Xa）来命名要作为模型错误处理的某些特定变量，在我们的情况下，应该将其用作SetModelErrorVariable（4，因此意味着总共6个线性扰动中的最后两个必须作为模型误差变量处理。6的线性扰动总数由正和零李雅普诺夫指数以及两个模型误差变量之和给出。假设我们还想考虑前两个李雅普诺夫向量之间的非线性相互作用的影响，因此我们设置非线性相互作用扰动的数量ml2从而在演化过程中得到m1（m1）6个非线性扰动。因此，演化过程中的扰动总数为6 6 12、在这个 道：非线性演化用于建立等于m l（m l）的分析的非线性扰动的数 量1）/2三是这样安排的：xaxa. . .Xa0 1 3xax a4 5xax ax a6 7 8代码的核心是EKF_AUS类，其正确方法如下所示。我们使用[10]的统一符号。我们还报告了用于非线性和模型误差项的名称。N（）：系统的自由度（包括受误差影响的参数）;P（）：尺寸的的测量linM（）：以线性化方式演化的扰动数（包括受模型误差影响的参数）;ml，NonlinM（）：非线性相互作用的扰动数;• NumberNonLinPert（）：非线性演化中涉及的扰动数等于ml（ml+1）;⎠δx···xf（τ）−xf（τ）−2[δxf（τ）+δxf（τ）]、SR=0Q=00⎛xaxa. . .Xa0 1 3xax a4 5xa. . .. . .Xa6 11δxa。. . δxa001，11，mδxa。. . δxa0. . .. . .. . .0δxa。. . δxaN，13，m000 0 0σV00 0 00σγL. Palatella，F.Grasso/SoftwareX 7（2018）2831=H··+•×•×•×==-=-·++•× [++]====-=≥Fig. 1. Lorenz96模型的误差均方根，N40，σo0。1，p20和τ0。125对于不同的m l值。很明显，对于m l3，4，滤波器是稳定的，并且不像线性（m l0）和m l<3情况那样发散。HalfNumberNonLinPert（）：这是在由ml（ml1）/2给出的分析步骤中使用的补充扰动的实际数目TotNumberPert（）：由总和m ml（ml1）给出的演化扰动总数η，Lfactor（）：η的值。所有这些参数都必须通过不同的构造函数进行设置（有关详细信息，请参阅doxygen文件主方法显然称为Assimilate，并且具有不同的输入，即：measure：测量值的p1矩阵; NonLinH：指向实现非线性测量运算符的例程的函数指针。注意，不需要线性化的测量算子;p1矩阵R包含测量误差协方差矩阵的对角项（我们作了简化假设，即不同测量值不相关），n1矩阵xf包含演化后的状态（预报）;包含分析前扰动的n m ml（ml1）/2矩阵Xf在输出上，xf被分析xa代替，而扰动Xf被分析扰动Xa代替。从系综成员传递到扰动Xf以及从扰动X f传递到系综成员所需的计算分别由例程PrepareForAnalysis和PrepareForEvolution执行。如果用户想将我们的算法应用到他们的系统中，他们不需要修改这些例程。相反，他们需要修改例程ReadFile，WriteFile和Evolve。前两个函数将系统的状态文件转换为一个单列状态向量（Eigen类的实例[16]），反之亦然。evolve例程通过使用运行动态系统软件的脚本调用系统2.1. 软件构架EKF-AUS-NL的这种实现可以应用于几个系统，由适当的动力学方程描述。我们在下面描述了应该实现的方法用户必须写下一个自定义类，其中包含接口类IAssimilate中声明的所有正确方法的定义这些方法考虑用户希望在同化测试中使用如前所述，这些例程是evolve、ReadFile、WriteFile、Dimensions（给出从模型的状态文件开始的状态向量的维数N与此例程一起，必须编写主过程例程，该主过程例程可以是例程完美或实际测试过程下的完美模型测试（如在所呈现的在第二种情况下，只需消除真值状态（显然是未知的），并写下一个能够从输入设备读取实际测量值的例程。在本文所介绍的软件中，我们展示了Lorenz96（L96）和SLAM算法的两个在接口IAssimilate中，声明的函数都是虚拟的，以便在编译后执行后期绑定。编译器生成vtable，在main（）例程中，编译后的软件根据所选的模型在运行时实例化正确的类我们参考doxygen文档的软件的其他实现细节和类层次结构。同化/进化循环的整个过程的完整方案如图所示。32.2. 软件功能该软件提供了两个主要功能。进行完美的模型情景测试，以验证测量和同化方案是否能够通过一组适当的实际测量来控制动力系统和实际同化实验的可能性这两种方法之间的选择必须在编译之前在主例程中完成。3. 说明性实例在本文中，我们展示了两个示例。第一个是Lorenz96模型[4]。在不使用线性切线模型的情况下实现我们的算法会导致与[3]中提出的完全相同的结果。此外，关于滤波器稳定性的结果与前一篇论文中获得的结果相同，因此我们没有再次报告它们，在这种情况下不是那么原始。在图1中，我们报告了Lorenz96模型分析误差的均方根（RMS）的行为。结果参考了针对不同ml值进行的完美模型试验。实验是用N40个自由度，σo0的情况。1，p20（我们每两个变量测量一个）τ 0。125（两次连续分析之间的时间），对于ml0、 1、 2、 3、 4的不同值。结果表明，在线性情况下，滤波器迅速发散。对于ml1， 2发生相同的情况，而对于ml3，滤波器相对于非线性效应是稳定的。有关模型的其他详细信息，请参见[3]。作为第二个例子，我们提出了SLAM算法，如在[5]中实现的，具有已知的和良好区分的参考地标。当然，在SLAM的实际应用中，必须应用合适的算法来识别在不同时间进行的什么测量属于同一地标（参见例如[5]的附录II在这里，我们想集中在应用于SLAM的EKF-AUS-NL算法的主要属性出于这个原因，我们认为已经建立了地标识别关于模型误差，我们认为两个变量，即速度（V）和转向角（γ），不是由驱动装置精确因此，我们增加了一些扰动，这些扰动等于受模型误差影响的变量的数量，然后将其视为两个额外的状态变量。在这里，我们显示作为一个例子的SLAM系统的运行与姿态的误差的时间演变。在提出的例子中，连续同化之间的时间间隔τ= 0。05秒，界标距离误差σd=20 cm，角度误差σθ=3cm。模型误差设定为σV=5cm/s和σγ=1cm/s。车辆的速度设置为V= 1。0米/秒。线性因子为η=0的情况。1，而扰动的数目被设置为m=5， 7， 10， 20·32L. Palatella，F.Grasso/SoftwareX 7（2018）28===-==图二. 上图：对于不同的扰动数量m值，RMS作为时间的函数。（下图）：使用SLAM的完美模型测试的示例。我们用真实（红线）、分析（黑色空点）、实际（空圆圈）和分析（点）标志显示了x-y这必须与系统N85（具有40个界标）的整个维度进行比较。 已经和我滤 波 器工作得非常好，这表明假设不稳定和中性的李雅普诺夫指数对应于变量x、y、V、φ（车辆相对于x轴的角度）加上对应于模型误差变量V、γ的两个变量的线性组合。该算法表现出很好的精度和稳定性，所以在这种情况下，没有必要考虑非线性相互作用。 在图2中，我们显示了实际和估计的轨迹以及40个地标的估计和位置。在该示例中，车辆遵循圆形轨迹。在上图中，我们将车辆位置的误差显示为函数在同化周期中测量的时间（两个周期之间的时间τ 0. 05 s）和扰动次数m.重要的是要强调的是，结果实际上是独立的扰动考虑M5， 7， 10，20的数量。此外，考虑到最小数量的矢量，滤波器能够控制系统。当处理具有多个地标的情况时（例如在自动驾驶应用中），此功能可能非常重要如果你有L个地标，L 1，我们有维数N因此，L图三. 软件包中代表进化-同化循环的图表。由于Anna Trevisan及其同事的开创性工作，几乎十年前[18]发表的这篇文章中，很少有在真实系统上的实际在我们看来，这是由于缺乏一个软件包，帮助不同的研究人员在使用这一重要的技术。特别是，我们认为这个软件是重要的，因为它提出了两个重要的特点：一个简单的软件结构，它允许在不知道算法细节（特别是那些与矩阵公式有关的细节）的情况下实现滤波器，而只知道该方法的主要思想;保持同化和动力系统代码之间完全分离的软件结构。有了这个软件包，动力系统软件通过调用命令系统来执行。因此，模型代码可以使用任何语言编写我们认为，这套软件包将主要在两个科学界推广。第一个是由参与大气和海洋模型数据同化的人员组成，第二个是面临自动导航，自动驾驶等主题的科学家。5. 结论本文介绍了一个包含C++程序的软件包，它实现了任意动力系统的EKF-AUS-NL算法，而不需要微分方程的雅可比矩阵。该软件还可以应用于具有参数误差的模型，并可以考虑非线性矩阵P和Pf中元素的个数2号 线，不同李雅普诺夫向量之间的相互作用，如[3]中所述对卡尔曼滤波器的性能的禁止性后果相反，在EKF-AUS应用于SLAM算法中，仅m5个扰动就足以控制系统，以这种方式极大地提高了滤波器的性能4. 影响本文提出的算法在不同的论文中使用[1最重要的是，以前版本的软件需要切线线性模型，而现在版本的代码不需要。在我们看来，这个软件的影响将非常广泛。尽管澳大利亚的数据同化方法是第一个我们认为，所提出的结果非常重要。事实上，在极少数情况下，当使用数据同化算法时，观测值足够密集和准确，以保持线性区域中的误差为了处理这一点，我们提出了文[3]中的非线性推广然而，所遵循的方法需要微分方程的解析表达式，这是方法几乎是不适用的。在这里，我们提出了一种方法来处理非线性没有切线线性模型。这一点有些微妙和不平凡。事实上，我们意识到[3]的附录B中遵循的方法是不正确的，我们修改它以恢复正确的结果。本文给出的两个例子可作为混沌大系统（如数值天气预报模型）和导航··L. Palatella，F.Grasso/SoftwareX 7（2018）2833算法类似于SLAM示例。我们希望这个软件将铺平道路的实际应用的AUS的想法，数据同化问题在几个动力系统。致谢谨以此文纪念安娜·特维桑.几年前，她推出了AUS方法在数据同化连同其他几个非常重要的课题在这一领域的研究。她总是用自己的才华、激情和对发现的热爱我们衷心感谢G. Lacorata提供有用的建议。引用[1] TrevisanA，Palatella L. 混沌与天气预报：不稳定子空间在可预报性和状态估计问题中的作用。Int J BifurcationChaos2011;21（12）：3389-415.[2] Palatella L，Carrassi A，Trevisan A.李雅普诺夫向量与不稳定子空间同化：理论与应用。 JPhys A 2013;46（25）：254020.[3] Palatella L，Trevisan A.李雅普诺夫向量在卡尔曼滤波器非线性扩展公式中的相互作用。 PhysRev E 2015;91（4）：042905。[4] Lorenz EN.可预测性：部分解决的问题。载于：《可预测性研讨会论文集》，第1卷。一九九六年。[5] 张文龙，张文龙.同时定位和地图构建（SLAM）问题的解决方案。 IEEETransaRobot Autom2001;17（3）：229-41.[6] AhmadH，Namerikawa T.基于扩展卡尔曼滤波器的移动机器人间歇测量定位。Syst Sci Control Eng An Open Access J2013;1（1）：113-26.[7] 黄S，迪萨纳亚克G.基于扩展卡尔曼滤波的SLAM算法的收敛性和一致性分析。IEEE Trans Robot2007;23（5）：1036-49.[8] 放大图片作者：J.离散卡尔曼滤波算法的数值比较：轨道确定案例研究，喷气推进实验室技术备忘录33-771。一九七六年[9] TrevisanA，Palatella L. 在卡尔曼滤波器中，误差协方差陷入不稳定子空间。非线性过程地球物理2011;18（2）：243-50.[10] Ide K，Courtier P，Ghil M，Lorenc AC.数据同化的统一符号：操作、顺序和变分（气象学和海洋学中的数据同化特刊：理论与实践）。J Meteorol SocJapan SerII1997;75（1B）：181-9.[11] Bocquet M，Carrassi A.四维系综变分资料同化与不稳定子空间。 TellusA：Dynamic Meteorology and Oceanography2017;69（1）：1304504.[12] Bocquet M，Gurumoorthy KS，Apte A，Carrassi A，Grudzien C，JonesCK. 退化 卡尔 曼滤 波误 差协 方差 及其收 敛于 不稳 定子 空间 。SIAM/ASA JUncertain Quantif2017;5（1）：304-33。[13] Toth Z，Kalnay E. NMC的包围预报：扰动的产生。《美国气象学通报》，1993年;74（12）：2317-30。[14] TothZ，Kalnay E.NCEP中的围扰预报及其育种方法。《星期一天气评论》1997年;125（12）：3297-319。[15] 雅兹温斯基随机过程与滤波理论。中国科学出版社.哈考特·布雷斯·约万诺维奇出版社的子[16] Guennebaud G等人，Eigen v3，http://eigen.tuxfamily.org。2010年。[17] [10]李文辉，李文辉.交通模型的非线性稳定性和利用李雅普诺夫向量估计交通状态。Phys Rev E2013;88（2）：022901。[18] 作者：A，M，O.不稳定子空间中的四维变分逼近与最优子空间维数。Q JRMeteorol Soc2010;136（647）：487-96.