基于解耦EKF-DSF估计的强鲁棒性在线MPC跟踪时变高超音速轨迹方法

《Acta Astronautica》：Robust online MPC tracking of tychastic hypersonic trajectories via decoupled EKF-DSF estimation

　　摘要：高超音速再入飞行器面临着严峻的空气热效应限制以及复杂的环境不确定性。本文针对在双重不确定性条件下实现稳健在线跟踪的问题展开研究，这两种不确定性包括较大的初始状态误差以及极端的大气密度偏差。首先从数学上证明了，由于瞬时不可观测性及外推误差的放大效应，传统的耦合估计方法并不适用于这项任务。为克服这些理论限制，本文提出了一种解耦架构，该架构包含一个运动学扩展卡尔曼滤波器以及专用的密度尺度因子估计器。这种多速率估计方案为降阶模型预测控制器提供输入，后者则依据离线生成的稳健参考轨迹进行跟踪。通过利用参考轨迹所具有的概率安全裕度——这一裕度已通过坎泰利不等式得到验证——所提出的模型预测控制器能够在无需进行复杂的在线约束管理的情况下保证热安全。大量的蒙特卡洛模拟结果证明了该架构的有效性：与基于标准大气模型的基准控制器相比，后者在终端阶段的偏差较大，平均高度误差超过4公里，而所提出的架构则实现了更好的跟踪性能，平均误差约为-63.8米。更为重要的是，其性能与具备完美大气信息理论的系统相比并无显著差异（平均误差约-56.2米）。这些结果表明，解耦框架能够将运动学噪声与参数偏差分开处理，为安全、自主的高超音速飞行提供了可行的解决方案。

1. 引言
在大气层再入引导方面，由于空气动力学、大气不确定性以及飞行器控制能力之间的紧密关联，它始终是航天任务设计中最具挑战性的问题之一[1]、[2]。历史上，再入引导依赖于增益调定的跟踪法则，其中最著名的当属阿波罗计划和航天飞机所使用的阻力能量算法，这类算法依靠离线生成的预计算参考轨迹来引导飞行[3]。随着对精确着陆、更强鲁棒性以及热安全要求的提升，模型预测控制逐渐成为该领域的主流方法。通过在飞行器上解决有限时间最优控制问题，模型预测控制为处理非线性动力学、强烈的状态-输入耦合关系以及严格约束条件提供了有力的理论框架[3]、[4]、[5]。

在基于模型预测控制的再入引导技术方面，近期取得了两大进展：一是提高了计算的可行性，二是增强了抗干扰能力。为减轻 onboard优化带来的计算负担，有研究提出了模型预测引导框架，用全局收敛的固定点算法替代了复杂的二次规划求解器[6]、[7]。在此基础上，又有研究将这一理念拓展到大气进入阶段，提出了一种两阶段策略，用以解决倾斜入轨飞行器中存在的控制力不足问题。与此同时，人们也致力于降低环境不确定性带来的影响。早期的干扰观测器虽然能够提供有效的前馈补偿，但并未考虑执行器饱和问题。为解决这一缺陷，有研究将状态空间观测器引入非线性模型预测控制框架中，从而实现对大气不确定性以及执行器限制条件的联合建模。最近，还有研究通过基于预定时间控制律和障碍李雅普诺夫函数的方法，确保在受限不确定性条件下仍能实现有限时间内的跟踪性能。

尽管在计算能力和反馈鲁棒性方面取得了诸多进展，但仍然存在一个根本性的问题：通常情况下，参考轨迹被视为固定的基准。这些方法隐含地假设在线控制器能够充分补偿环境偏差。然而，当名义轨迹接近诸如加热率或动压极限等关键安全边界时，即便模型预测控制器运行正常，也可能无法在非正常工况下保持系统的可行性。正是这一认识推动了稳健轨迹生成技术的发展。其中，多项式混沌展开法作为一种强大的工具，被用于通过非线性动力学传递不确定性。有研究证明了这种无侵入式的多项式混沌展开法在轨道不确定性传播方面的有效性，还有研究指出，将多项式混沌展开法与凸优化相结合，能够显著提高计算效率。与谱方法相辅相成的是，基于离散化的概率方法也被开发出来，用于处理航天器轨迹规划中的随机约束问题。有研究通过平滑似然近似和轨迹集合，将随机轨迹优化问题转化为确定的数值优化问题。

虽然这些技术都非常有效，但罗斯提出的时变最优控制框架提供了一种更具吸引力的替代方案，因为它无需使用伊藤微积分和随机微分方程。这一框架源自黎曼-斯蒂尔杰斯最优控制理论[17]、[18]，后来发展成了无迹轨迹优化方法[19]。该方法通过无迹变换将概率分布映射为确定的西格玛点[20]、[21]、[22]，从而能够在标准的数值优化算法中直接设置均值-方差安全区间。鉴于其在控制轨迹分散方面的优异表现[23]、[24]，本文采用了这一框架来生成所提出的引导架构所需的稳健参考轨迹。

尽管在规划与跟踪技术方面都已取得进展，但两者之间仍存在重要的集成缺口。现有的基于模型预测控制的再入引导架构往往是在假定下降过程中不确定性逐渐减小的扰动模型下进行验证的。例如，有研究中使用的大气扩散模型就假设在较低高度时密度偏差会显著减小，从而减轻了在最敏感阶段的在线估计负担。然而，实际运行数据表明，真实的大气变化要严重得多。正如有研究指出的那样，高层大气对太阳辐射通量和地磁活动极为敏感，导致每日密度不确定性超过50%[27]。这类变化通常并非高频波动，而是持续存在的较大规模偏差，往往难以通过简单的过滤策略加以消除[28]。当如此大的大气偏差与显著的初始状态误差叠加时，就会出现所谓的“双重不确定性”情况，此时传统的耦合状态-参数估计器会变得不稳定，而线性化控制器则有可能出现灾难性的发散现象。

为克服这些限制，本文提出了一种专为在严重不确定性条件下实现稳健大气层再入而设计的解耦估计-控制架构。其主要贡献如下：
1. 从数学上证明了，由于实际观测限制，通过耦合状态-参数方案直接在线估计大气参数是不稳定的。这一理论发现为采用降阶参数跟踪方法提供了依据。
2. 提出了一种多速率估计框架，该框架将运动学状态估计与环境识别严格分离。通过将导航解与空气动力学建模误差隔离开来，这种设计能够防止在极端大气偏差情况下估计器出现发散现象。
3. 建立了一种用于约束管理的责任分离策略。通过利用离线无迹轨迹生成技术预先设定概率安全裕度，在线模型预测控制器无需承担处理随机约束优化所带来的高昂成本。大量仿真结果表明，该架构能够将热违规率降至可接受的工程水平，同时其跟踪精度与具备完美大气信息的理想控制器相当。

2. 问题描述
整个问题分为两个阶段：第一阶段是在存在较大大气不确定性的情况下离线生成稳健的参考轨迹；第二阶段则是使用模型预测控制框架实时跟踪该轨迹，而该框架必须能够估算出真实的大气条件。

2.1 全三维时变模型
为准确处理高超音速飞行中固有的不确定性，本文采用了时变最优控制理论框架，该框架不同于基于伊藤微积分的传统随机最优控制方法[29]。传统的随机方法通常使用由维纳过程驱动的随机微分方程来描述不确定性，这需要运用复杂的伊藤微积分，并且还需要诸如欧拉-马鲁亚马法之类的特殊数值积分方法[30]。然而，如有的研究所指出的，高超音速飞行轨迹上大气密度的变化实际上属于参数不确定性，而非维纳过程所具有的那种高频、剧烈的噪声特性。因此，这里的不确定性是通过时变随机变量p来描述的，这些随机变量定义在某个支撑集supp(p)上，而非动态噪声过程[17]、[19]。这样的建模方式使得系统在特定p取值情况下具有确定性：对于p的任何具体取值，系统的动态行为都是确定且平滑的，这样就可以使用标准的龙格-库塔求解器以及成熟的伪谱优化技术，同时还能有效地处理统计分散性问题。该系统由以下参数依赖的常微分方程描述：

x?3dof(t)=f3dof(x3dof(t),u(t),p)

其中，全三维状态向量x3dof=[r,θ,?,V,γ,ψ]T∈R6，用于描述飞行器在以行星为中心的旋转坐标系中的位置和速度，各分量分别表示径向距离、经度、纬度、相对于行星的速度、飞行路径角度以及航向角。控制向量u=[α,σ]T，则代表攻角和侧倾角。

核心的不确定性来自于大气密度模型，该模型由时变向量p∈R2参数化：

p=[ρ0,β]T～N(pμ,Σp)

忽略行星自转Ω（尽管它在系统模型中有所体现），三维动态方程可表示为[32]：

r?=Vsinγ

θ?=Vcosγsinψrcos?

??=Vcosγcosψr

V?=?Dm?g(r)sinγ

γ?=LcosσmV?g(r)V?Vrcosγ

ψ?=LsinσmVcosγ?Vrcosγsin?tan?

其中的空气动力系数L和D是未知参数向量p的函数，它们与密度ρ有关：

L=12ρ(h;p)V2CL(α)Sref

D=12ρ(h;p)V2CD(α)Sref

式中，h=r?Re，CL(α)=CL0+CL1α，CD(α)=CD0+CD1α+CD2α2，而密度则由以下指数模型给出：

ρ(h;p)=p1exp(?p2?h)=ρ0exp(?βh)

注1：在该框架下，准确的密度建模至关重要，因为空气动力系数与密度成正比。虽然像Mars-GRAM或Earth-GRAM这样的高精度工程模型能够捕捉复杂的大气现象[33]，而基于LSTM网络的数据驱动方法也具备很强的重建能力[34]、[35]，但它们要么计算成本过高，不适合用于在线轨迹优化，要么缺乏梯度型求解器所需的解析导数。诚然，指数假设会引入结构偏差，但这种形式能够体现轨迹规划所必需的主要尺度高度特征。更重要的是，该模型将密度分布的无限维不确定性降到了低维参数空间p=[ρ0,β]T中，从而使所提出的EKF-DSF架构下的在线估计问题变为可观测问题，同时也使得稳健轨迹优化在计算上变得可行。

2.2 降阶模型
对于实时模型预测控制而言，使用全三维时变模型不仅计算成本极高，而且由于倾斜入轨飞行器的固有特性，其动态特性也不稳定。正如有研究指出的，通过侧倾控制实现大气层进入属于控制力严重不足的问题，因为一个主要的控制变量（侧倾角σ）需要同时调控六个相互关联的状态量。如果试图同时跟踪纵向（距离、能量）和横向（横向距离、航向）动态，往往会引发控制需求之间的冲突——为纠正距离误差而需要的侧倾角会与保持航向的要求相矛盾，从而导致跟踪性能下降甚至系统无法正常工作[36]。

为解决这一问题，本文采用了解耦架构。在线模型预测控制器使用降阶的纵向模型，该模型仅关注平面内的动态行为（能量管理和距离控制）。这种设计利用了这样一个事实：纵向轨迹的变化主要取决于垂直升力分量Lcosσ的大小，而与侧倾方向无关[37]。横向动态则单独处理，用以控制σ的正负值。

由此产生的纵向状态向量x=[r,V,γ]T由以下方程描述：

r?=Vsinγ

V?=?D(x,α,p)m?g(r)sinγ

γ?=L(x,α,p)cosσmV?g(r)V?Vrcosγ

其中，空气动力系数L和D仍然明确依赖于时变参数向量p，这使得纵向模型预测控制器能够对由EKF-DSF估计器得到的大气不确定性保持稳健性。

注2：本文仅聚焦于纵向生存问题，即在极端密度不确定性条件下防止热失控以及飞行器跳出轨道或坠毁。横向引导则属于常规的几何问题，通常通过独立的侧倾反转逻辑来处理[38]。由于横向机动并不会从根本上改变能量-密度估计流程，因此为了突出所提出的估计架构的作用，本文在模型预测控制的公式推导和仿真分析中均忽略了横向逻辑。

2.3 在线估计问题与测量模型
虽然离线优化能够在概率分布p～N(pμ,Σp)下得到理论上较为稳健的参考策略uref(t)，但在线控制问题则是确定性的，不过存在不确定性。飞行器面临的是大气参数ptrue的特定且未知的取值，该取值很可能与先验均值有很大偏差，同时还会存在不可避免的初始状态误差x(t0)≠xref(t0)。因此，在线控制的核心目标就是设计一种在计算上可行的估计架构，能够实时重建无法直接观测到的大气密度，从而让模型预测控制器能够在真实的动态条件下跟踪参考轨迹xref(t)。

为实现这一目标，估计器利用两类不同频率的 onboard测量数据为低频率的模型预测控制器提供输入。至关重要的是，估计器所能获取的信息被分为与空气动力学相关的部分和与之无关的部分。惯性测量单元能够提供感测到的特定力值，这些力值实际上反映了飞行器对未知大气的空气动力学响应：

zimu=aD,measaL,meas=1mD(x,u,ptrue)L(x,u,ptrue)+wimu

其中，wimu～N(0,Rimu)代表传感器噪声。而导航系统（如GPS/INS）则提供与大气参数无关的运动学状态测量值：

znav=rmeasVmeas=rtrueVtrue+wnav

其中，wnav～N(0,Rnav)。将测量数据分为与参数相关的力值（zimu）和与参数无关的运动学数据（znav），这一划分构成了后续章节中提出的解耦EKF-DSF架构的物理基础。

3. 离线稳健轨迹生成
为使随机最优控制问题在计算上可行，本文采用了无迹变换技术，将大气参数p=[ρ0,β]T的连续概率分布映射为一组有限的确定性样本点（即西格玛点）[26]。这一过程将单一的随机轨迹优化问题转化为了确定的超轨迹优化问题。

3.1增广确定性问题中采用了单纯形西格玛点集，该方法能够用Nσ=d+1或2d个点来最小化表示d维随机变量的前两个统计矩——均值和协方差。对于二维不确定性情况（d=2），基于归一化的单纯形集合χ∈R2×4来确定4个西格玛点，该集合由公式(17)χ=200?202?20定义，每个点对应的权重wi为1/4。物理参数实值Pi=[ρ0,i,βi]T是通过将这些归一化点通过假设的均值pμ和协方差Σp进行变换而得到的[39]，公式为(18)[P1,…,PNσ]=pμ11×Nσ+Σpχ。无迹框架的核心创新在于，在单一的共同控制策略u(t)下，同时传播这Nσ个不同的大气状态实值。定义了一个增广状态向量X(t)∈Rnx×Nσ，它将所有西格玛点的状态轨迹拼接在一起(19)X(t)=[x1(t)T,x2(t)T,…,xNσ(t)T]，这样动力学问题就被扩展成了Nσ个相互独立的常微分方程系统，这些方程仅通过共享的控制量u(t)相互关联(20)x?i(t)=f6dof(xi(t),u(t),Pi)，i=1,…,Nσ。

3.2 统计约束转换与理论风险界限随机加热率约束被转换为基于集合统计的确定性代数路径约束。对于每一个西格玛点i，瞬时加热率可通过公式(21)q?i(t)=C(α)ρ(ri,Pi)0.5Vi3.07计算。然后利用无迹权重来计算整个集合中加热率的瞬时均值μq?和标准差σq?，公式分别为(22)μq?(t)=∑i=1Nσwiq?i(t)(23)σq?(t)=∑i=1Nσwiq?i(t)?μq?(t)2。这些通过统一处理技术计算得到的矩值可以近似真实加热率分布的矩值。由于加热率模型对p是非线性的，所以在高斯输入假设下，西格玛点加权求和只能近似得到精确矩值，而非封闭形式的解析解。因此，所得到的约束实际上是基于统一处理技术的、无需考虑分布的坎泰利安全准则[40]的实现形式。随机加热率约束被转化为了确定性代数路径约束。按照NASA热防护系统设计的标准做法，会区分操作规划限制值(q?plan)和最终结构极限值(q?limit)[41]、[42]。根据火星2020任务中的气动热裕度设计原则，会在预测的对流加热环境基础上乘以1.1到1.23之间的系数[41]，因此这里采用了安全因子ηsafe=1.23。这样一来，离线优化就会保守地以q?plan=794.5kW/m2作为目标值，而结构失效阈值则为q?limit=ηsafe?q?plan≈977kW/m2。因此，所采用的鲁棒替代约束为(24)μq?(t)+κσq?(t)≤q?plan，?t∈[t0,tf]。虽然通常情况下κ=3的选择是基于高斯假设的经验值（意味着风险约约为0.13%），但由于高超音速动力学的非线性特性，实际分布并非高斯分布。为确保绝对安全，这一选择是通过坎泰利不等式（单边切比雪夫不等式）[40]来证明其合理性的。对于任何具有有限均值和方差的加热率分布Y=q?(t;p)，只要满足约束方程(24)，其实际违规概率就会被严格限制在(25)Pr(Y>q?plan)≤11+κ2之内。这一界限与具体分布无关。当κ=3时，根据公式(25)可得出违规概率θ≤1/(1+32)=10%。在当前基于统一处理技术的实现方案中，这一数值就被作为名义上的安全界限。通过1000次离线蒙特卡洛模拟，得到的实际违规率约为2.5%，这一数值与前述安全准则的保守性是一致的。在第五章第3节所描述的在线闭环蒙特卡洛模拟中，情景B的热违规率为9.4%，即便存在在线跟踪动态带来的额外扰动，这一数值依然低于10%的名义安全界限。这些结果表明，基于离线统一处理技术设定的安全裕度在在线跟踪阶段依然有效。为了对上文推导出的理论界限进行实证验证，研究人员开展了包含1000次分散模拟的高精度蒙特卡洛模拟。大气参数是从完整的3σ不确定性范围内抽样的，而飞行器动力学则依据所得到的控制策略来进行推演。图1展示了不同轨迹下的加热率对比曲线，从中可以看出两种方法在安全性能方面存在显著差异。如图1(a)所示，确定性基准控制策略会紧紧跟随约束边界（0.79MW/m2）以尽可能提升性能。然而，这种缺乏安全裕度的设计在非正常工况下会带来致命后果，几乎在所有模拟场景中都会出现严重的约束违反现象，使得实际飞行轨迹无法实现。相反，鲁棒的无迹优化方案（图1(b)）则能有效将加热率分布拉离约束极限。通过强制施加替代约束（μ+3σ≤q?plan），优化器会自动形成必要的安全缓冲区。因此，蒙特卡洛分析证实，违规率被大幅降低到了2.5%。这一实际违规率不仅与[26]的研究结果相符，也严格处于公式(25)所保证的理论安全界限（≤10%）之内。当优化过程收敛后，就会得到一个统一的鲁棒控制策略u?(t)，以及与每一个大气状态实值相对应的Nσ条独立状态轨迹。传递给在线跟踪阶段的名义参考轨迹xref(t)则是这些集合轨迹的加权平均值(26)xref(t)=∑i=1Nσwixi?(t)，uref(t)=u?(t)。这样的设计确保了模型预测控制算法能够追踪飞行器路径的统计期望值，而该期望值恰好位于由西格玛点所界定的安全飞行走廊的中心位置。备注3：严格的离线验证表明，xref(t)具备与实际飞行风险相匹配的固有概率安全裕度。这一重要发现为在线控制策略的变革提供了依据：无需让机载计算机承担计算成本高昂的带机会约束优化任务，在线控制器只需专注于精确跟踪即可。只要在线估计器能够收敛到真实的大气状态（详见第4节），严格追踪这一鲁棒参考值就足以确保热安全。下载：下载高分辨率图像（317KB）下载：下载全尺寸图像图1. 不同轨迹下的蒙特卡洛加热率曲线（N=1000）。4. 在线模型预测控制与估计架构在确立了鲁棒的参考轨迹xref(t)之后，接下来就要介绍用于追踪该轨迹的在线架构。为应对双重不确定性问题，本文提出了一种分离式的估计与控制架构。这种设计并非随意为之，而是基于对传统耦合式估计方案为何会在这种特殊的高超音速环境下失效的分析而严格推导出来的。下载：下载高分辨率图像（381KB）下载：下载全尺寸图像图2. 系统架构图，展示了两阶段处理方式：离线阶段的鲁棒轨迹生成（上方）以及在线阶段的多元速率分离式跟踪（下方）。4.1 理论分析与提出的估计架构对于大气不确定性的常规处理方法是，向状态向量x中加入未知参数p=[ρ0,β]T，然后再采用递归卡尔曼滤波器或批量加权最小二乘法估计器。然而，仔细分析后会发现，对于入轨飞行时的空气动力学问题而言，这种处理方式存在根本性的条件不适问题。具体来说，首先，卡尔曼滤波器存在结构上不可观测的问题；其次，当用于反解ρ0时，对数线性加权最小二乘法估计器会变得数值上不稳定。以下将对这些问题进行正式阐述。定义大气密度模型ρ(h)=ρ0exp(?βh)，并考虑由惯性测量单元获得的空气动力学加速度测量值zimu。下面的结论表明，ρ0和β是无法同时被估计出来的。引理1：空气动力学观测量的梯度共线性对于模型ρ(h)=ρ0e?βh，其测量灵敏度满足(27)?zimu?β=(?ρ0h)?zimu?ρ0，因此对于所有h>0的情况，这些灵敏度都是完全共线的。这样一来，关于(ρ0,β)的观测雅可比矩阵的秩就为1。证明：设zimu(h)=g(h)ρ(h)，其中g(h)包含了所有已知的空气动力学项。那么有(28)?ρ?ρ0=e?βh，(29)?ρ?β=ρ0e?βh?h=(?ρ0h)?ρ?ρ0。将g(h)乘到等式两边后，共线性依然保持不变，至此证明完毕。□命题1：大气参数的不可观测性由于引理1中所提到的雅可比矩阵的各列之间存在线性依赖关系，因此仅凭空气动力学测量数据是无法观测到参数对(ρ0,β)的。如果试图用卡尔曼滤波器同时估计这两个参数，那么其协方差矩阵就会出现发散现象。这就解释了为什么在短距入轨飞行过程中，无法区分表面密度ρ0和尺度高度β的变化：从数学上来看，滤波器实际上只能“感知”到一个有效的单一尺度因子。备注4：引理1和命题1都表明，在任意某一特定高度h处，关于(ρ0,β)的空气动力学观测雅可比矩阵的秩均为1；也就是说，仅凭单次测量数据就无法立即观测到这两个参数。在整个入轨飞行过程中，飞行器会经过从大约80公里到24公里不等的不同高度，而共线性因子(?ρ0h)的值会随着高度的变化而变化。原则上，如果在足够不同的高度处进行测量，那么就能够获得独立的梯度方向，从而在一定程度上恢复信息矩阵的秩。不过，有两个实际因素限制了这种跨全程的信息获取能力：第一，相对于整个飞行轨迹而言，实时估计所涉及的短暂时间窗口内的高度变化幅度相对较小，这就限制了相邻雅可比矩阵行之间的瞬时角度差；第二，指数型的密度分布会导致对数线性反演过程将即便是很小的测量噪声也会放大成ρ?0方面的极大误差（见引理2）。正是这些因素共同作用，使得即便有来自多个高度的测量数据，耦合式(ρ0,β)参数估计仍然存在条件不适的问题，这也促使本研究采用了降阶的密度尺度因子公式。虽然使用对数线性化方法的批量加权最小二乘法[43]看似很有吸引力，但这种方法存在一种隐性的误差放大效应。对lnρ进行线性化处理后可以得到(30)ln(ρ(h))=ln(ρ0)?βh，这样一来，估计误差就会满足(31)ln(ρ?0)≈intercept+β?hmean。由于在入轨飞行过程中hmean的值通常在20–30公里左右，哪怕是β?方面的微小噪声，也会导致ρ?0出现极大的误差。具体来说：引理2：对数线性反演的不稳定性设δβ表示对数线性加权最小二乘法带来的估计误差，那么ρ?0的误差将会满足(32)δρ?0ρ?0∝exp(hmeanδβ)，随着hmean的增大，这一误差会无限制地增长，这就违背了加权最小二乘法所要求的残差为零的前提条件，进而导致数值计算失败。备注5：引理1和引理2共同表明，当(ρ0,β)这两个参数都具有不确定性时，无论是采用卡尔曼滤波器增强型方法还是对数线性化批量估计方法，都不适合用于高超音速入轨飞行场景。鉴于上述存在的问题，本研究没有采用直接的参数估计方法，而是采用了单一的、随时间变化的密度尺度因子（DSF）Kρ(t)，这一做法与航空捕获技术中的做法是一致的[44]。这一替换是通过一种分离式架构来实现的，该架构打破了参数误差与状态估计之间的反馈循环。第一个组成部分是一个高频运行的纯运动学卡尔曼滤波器。为确保状态估计值x?k不会因为大气建模误差而出现偏差，该滤波器的校正步骤仅限于使用那些与空气动力学无关的导航测量数据（即公式(16)中所给出的znav）。这样一来，测量雅可比矩阵H就变成了一个静态的选择矩阵[I2×2,02×1]，它将运动学解与动态模型中的不确定性分隔开了。第二个组成部分则是DSF跟踪器，它利用卡尔曼滤波器得到的纯净运动学状态来反解阻力方程。通过将实际测量的阻力加速度aD,meas与模型预测值进行比较，就可以确定出瞬时的密度偏差(33)ρ?est(k)=2maD,meas(k)CD(αk)SrefV?k2 (34)Kest(k)=ρ?est(k)ρnom(h?k)由于原始的Kest值存在噪声，因此需要通过一阶低通滤波器对其进行处理，以便提取出背后的大气趋势(35)K?ρ(k)=(1?Gdsf)Kest(k)+GdsfK?ρ(k?1)这种仅依赖单个参数的估计方法完全避免了共线性问题。得到的估计值会被送入模型预测控制算法中，该算法会利用ρpredict(h)=K?ρ?ρnom(h)来调整自身的内部模型，从而让控制器与真实的环境状态保持同步。4.2 省略在线加热率约束的理由在所提出的架构中，一个至关重要的设计决策就是刻意不在在线模型预测控制公式中加入明确的加热率路径约束条件（q?≤q?plan）。这一决定绝非疏忽，而是基于对加热率模型相对于大气参数而言具有严重非线性特性的严格风险分析后做出的。传统的实时模型预测控制实现方法通常依赖于对约束条件的线性化处理，以此来保证问题的凸性并提升计算效率。为了评估在这种极端大气不确定性条件下采用此类方法的合理性，研究人员开展了大规模的蒙特卡洛模拟分析（N=1000）。通过将完整非线性模拟得到的真实峰值加热率（q?true）与基于名义大气状态pμ的一阶泰勒展开式得到的线性化预测值（q?lin）进行对比，从而分析了参数不确定性所带来的影响(36)q?lin(t)=q?ref(t)+?pq?(t)|pμ?(ptrue?pμ)。图3展示了分析结果，从中可以明显看出这种线性近似方法的根本性缺陷。图3(a)显示，在靠近约束边界的高风险场景中，线性模型会系统性地低估热负荷。这种不保守的偏差会给人一种危险的虚假安全感，使得模型预测控制算法认为某个飞行轨迹是可行的，而实际上该轨迹已经超出了热安全限制。从定量角度来看，线性化误差的直方图（图3(b)）显示，在35.0%的模拟案例中，仅线性化误差就已经超过了0.04MW/m2的5%标准安全裕度。这项分析证实，加热率约束在参数空间中属于高度非凸的函数。因此，如果试图通过在线线性化的方式来强制执行这一约束，就会带来不可接受的安全风险。为了解决这个问题，本研究采用了责任分离的策略：将处理非线性热安全问题这一计算成本较高的任务交由离线鲁棒优化算法来完成（详见第3节），该算法会生成一个位于概率安全走廊中心的参考轨迹xref。这样一来，在线模型预测控制算法就无需再承担这一复杂的非线性任务，它的功能就被严格限制在精确追踪这一鲁棒基准轨迹上。下载：下载高分辨率图像（359KB）下载：下载全尺寸图像图3. 参数不确定性带来的影响（N=1000）。虽然在极端参数条件下严格追踪平均参考值可能会导致暂时的热负荷波动（详见第5.3节），但相比而言，这种策略能够避免因在线约束条件线性化处理错误而导致的控制器发散以及控制权丧失等更为严重的风险。4.3 多速率模型预测控制公式在线控制器以较慢的计算周期（Tmpc=1.0s）运行，用以解决有限时间范围内的最优跟踪问题。在每个时间步tk中，MPC利用当前状态估计值x?k以及第4节中给出的大气参数估计值K?ρ,k来计算控制序列，从而最小化与离线参考轨迹xref的偏差。该最优控制问题被构建为一个在预测范围N内的二次规划问题：minδX,δU‖δxk+N‖P2+∑i=0N?1‖δxk+i‖Q2+‖δuk+i‖R2+‖Δuk+i‖Rδ2s.t.δxk+i+1=Ak+iδxk+i+Bk+iδuk+i+dk+iumin≤uref,k+i+δuk+i≤umaxδxk=x?k?xref(tk)其中δx和δu表示与参考状态及控制的偏差，而Δuk+i=δuk+i?δuk+i?1则表示控制增量。为确保数值稳定性，尤其是当状态变量的量级差异极大时（例如，r～106米与γ～10?1弧度），权重矩阵会经过严格归一化处理。状态惩罚矩阵Q是通过物理比例因子（Sr=Re,SV=V0,Sγ=1.0弧度）来构建的：Q=diagwrSr2,wVSV2,wγSγ2这种归一化处理确保了高度上1%的相对误差与速度上1%的误差对成本函数的贡献相当，从而避免优化器偏向于处理量级较大的状态。控制惩罚项R=diag(rα,rσ)则用于对偏离离线参考值uref的情况进行惩罚。尽管在受扰动的大气条件下，参考控制可能已不再是最优的，但它仍可作为稳定的基准，确保在线解决方案仅在必要时偏离已验证的飞行轨迹。此外，还设置了控制速率惩罚项Rδ，以抑制高频振荡并保证控制的平滑性。最后，终端权重矩阵P的选择与阶段成本相匹配（即P≈Q）。这样的配置能够提供足够的终端惩罚力，引导预测轨迹在边界处趋近于参考轨迹，从而避免因处理高度非线性的高超音速动力学问题而产生的数值异常[45]。为确保在高度非线性的高超音速动力学环境中，能够在相对较长的控制时间内保持预测精度，这里采用了四阶龙格-库塔积分方法而非简单的欧拉近似法。降阶状态的离散时间演化公式为：xk+1=FRK4xk,uk,K?ρ,k其中FRK4(?)表示对连续纵向动力学方程（公式(12)–(14)）在一个时间步长内的数值积分。关键在于，此积分过程中的空气动力作用力是依据估计得到的大气密度ρ(h)=K?ρ,k?ρnom(h)来计算的，这样内部模型就能适应当前的环境状况。在构建该二次规划问题时，这个非线性离散映射会在离线参考轨迹(xref,uref)附近进行线性化处理。采用的是“先离散化再线性化”的方法，通过RK4积分器的中心差分法来计算雅可比矩阵：Ak=?FRK4?x|ref,k,Bk=?FRK4?u|ref,k公式(37)中的仿射项dk用于考虑参考轨迹本身的离散化误差。这种高精度的线性化处理确保了MPC预测出的未来状态与实际的非线性系统行为一致。完整的在线系统架构被实现为一个多速率反馈系统，如图2所示。该架构严格区分了不同的时间尺度：MPC优化在较慢的引导周期内运行，以控制计算负荷，而解耦的EKF和DSF估计器则在较快的周期内运行，以实现高带宽的状态跟踪。这一循环的详细执行流程见算法1。下载：下载高分辨率图像（235KB）下载：下载全尺寸图像4.4. 解耦架构的稳定性特性所提出的架构的一个显著特点就是严格分离了估计与控制的时间尺度：运动学EKF和DSF滤波器以快速的传感器更新频率Tf=0.1秒运行，而MPC则在较慢的引导周期Tmpc=1.0秒内解决最优跟踪问题。这种多速率结构，再加上将运动学测量结果与空气动力模型误差相分离，使得我们可以对其进行分析，分别确定各个子系统的有界性，进而整体了解整个系统的行为特性。接下来的讨论将从内部的估计循环开始，逐步过渡到外部的控制循环，最终为整个闭环架构提供实用的稳定性分析。这里的分析是在给定假设和特定进入窗口条件下的EKF-DSF MPC级联系统的局部有界性分析，并非旨在为任意飞行条件下的全非线性高超音速动力学系统提供通用的闭环稳定性保证。首先分析内部的估计循环。运动学EKF的设计使得修正步骤完全依赖于导航测量值znav（公式(16)），这些值是对运动状态（r,V）的直接测量结果，虽然存在噪声干扰，但与大气参数p完全无关。因此，测量矩阵H=[I2×2,02×1]∈R2×3实际上是一个恒定的、满秩的选择矩阵。设x?k=x?k?xk表示EKF的估计误差。以下命题阐述了由这种结构所带来的有界性特性。命题2 仅基于运动学信息的EKF的随机有界性考虑这样一个离散时间系统，其状态xk∈R3在纵向动力学方程（公式(12)–(14)）的作用下演变，而大气密度则是通过某个任意的、有界的比例因子K?ρ(k)>0来估算的。EKF的预测步骤使用相同的动力学模型，而修正步骤则采用线性测量模型znav,k=Hxk+vk，其中H=[I2×2,02×1]，且vk～N(0,Rnav)。假设满足以下条件：(A1)线性化后的动力学模型Fk=?f/?x|x?k是有界的，并且对于所有的k来说，组合体(Fk,H)都是可观测的；(A2)非线性动力学函数f(x,u,K?ρ)在状态x上具有利普希茨连续性，且在飞行范围内其利普希茨常数是有界的；(A3)初始估计误差x?0以及噪声协方差矩阵Qk、Rnav都是有界的。那么，估计误差的均方值将呈指数级有界：存在常数η>0和0<ν<1，使得对于所有的k≥0，都有(41)e[‖x?k‖2]≤ηνke[‖x?0‖2]+ρ，其中ρ的值取决于噪声协方差以及函数f的利普希茨常数。重要的是，由于h是一个恒定的选择矩阵，因此测量更新过程和可观测性结构都与大气参数p无关。虽然有界的大气参数偏差k?ρ(k)≠ktrue(tk)会通过空气动力作用影响预测步骤，但它只是以每个快速时间步长tf=0.1秒内的有界建模误差形式出现，并会被纳入有效的过程噪声中，从而不会破坏误差的有界性结构。证明在条件(a1)–(a3)下，参考文献[46]中的定理3.1可以保证通用离散时间ekf的指数级有界性。其核心原理在于测量模型是线性的且与参数无关（h为恒定且满秩）。在所考虑的进入窗口以及不断变化的速度-高度分布条件下，降阶后的运动学模型在模拟轨迹上仍然是局部可观测的；因此，可观测子空间对大气参数变化的敏感度极低。由于(a2)以及有限飞行范围内的空气动力作用力的有界性，单个快速时间步长tf=0.1秒内由大气参数偏差k?ρ≠ktrue带来的预测误差也是有界的；这类误差会被纳入ekf协方差传播过程中的有效过程噪声项中。因此，与参数无关的修正步骤主导了预测误差，避免了大气参数偏差与状态估计错误之间的正反馈效应，进而防止了耦合型估计器出现不稳定现象（见命题1）。□命题2的实际意义在于，即便大气模型存在严重误差，导航解依然能够保持可靠性；这一点在图4中通过数值实验得到了验证，图中显示，稳态下的ekf误差与极端的大气条件下的情况并无关联。之后，dsf跟踪器会使用经过过滤的运动学状态x?k来反演阻力测量值（公式(33)–(35)）。设ktrue(t)=ρtrue(h(t))/ρnom(h(t))表示真实的、随时间变化的密度比例因子，而k?ρ(k)=k?ρ(k)?ktrue(tk)则表示dsf的估计误差。dsf滤波器实际上是一个一阶线性时不变低通滤波器，其增益为标量gdsf∈(0,1)。该滤波器的极点位于z=gdsf处，且严格处于单位圆内部，因此能够保证有界输入下的有界输出稳定性。由于ekf的估计误差是均方有界的（见命题2），加之imu加速度测量值的噪声方差也是有界的，因此原始的密度比例因子估计值kest(k)本质上是一个有界的随机序列。经过过滤后的估计值会以指数级收敛到真实密度比例因子的附近区域：(42)??|k?ρ(k)?ktrue(tk)|≤gdsfk|k?ρ(0)|+11?gdsf(?snr+?drift)，其中??snr受传感器噪声水平以及瞬时的空气动力信噪比的影响，而??drift则反映了由于真实密度比例因子ktrue(t)在滤波器窗口内的缓慢变化所导致的跟踪滞后现象。这一特性可以解释图10(a)中所观察到的“宽-窄-宽”方差分布特征：在高空动态压力较大的飞行阶段，强烈的阻力信号会导致较低的??snr值，从而使置信区间较为狭窄；而在高空及末端阶段，由于空气动力作用较弱，方差则会有所增大。值得关注的是，整个过程中的偏差始终保持在零点附近，这正符合有界输入有界输出稳定性特性。接下来分析外部控制循环。需要明确的是，这里的讨论仅针对有限进入窗口条件下的局部实际稳定性问题，并非针对整个高超音速动力学系统的全局非线性稳定性证明。设δxk=xk?xref(tk)表示跟踪误差。降阶后的mpc模型（公式(37)）实际上是一个带有终端状态惩罚项p的有限时间范围跟踪问题。在无约束情况下，其解与有限时间范围的lqr问题一致；如果p满足离散代数riccati方程，那么最优成本函数v?(δx)=‖δx‖p2就可以作为闭环线性化系统的李亚普诺夫函数，从而保证跟踪误差的局部渐近稳定性[5]，[45]。在存在控制约束的情况下，此处采用的终端惩罚项p≈q实际上是对无限时间范围成本函数的二次近似。虽然这种选择并不能构成原始非线性高超音速动力学系统的严格控制李亚普诺夫函数，但它能够提供足够的终端成本，避免因有限时间范围的截断而导致控制策略出现短视或不稳定现象[5]。此外，由于在每个mpc迭代步骤中，ak、bk这两个矩阵都是通过rk4积分器的中心差分法重新计算的，因此线性预测模型能够准确捕捉沿参考轨迹的局部动态行为，从而减少可能降低系统稳定性的线性化误差。整个闭环系统实际上形成了一个级联结构：ekf和dsf估计器负责提供状态和密度估计值给mpc，而mpc则根据这些信息生成控制指令。由于估计循环的运行速度比控制循环快一个数量级（tf=0.1秒 vs. tmpc=1.0秒），因此由初始化或大气条件变化所引发的暂时性估计误差会在连续的mpc迭代之间迅速衰减。从mpc的角度来看，估计残差x?k和k?ρ(k)被视为有界的附加干扰。关于跟踪mpc的现有稳健性研究[5]表明，有界的输入干扰只会导致状态偏离名义轨迹一定的范围；这一特性被正式定义为输入到状态的实际稳定性isps[47]。因此，整个复合系统具备实际稳定性：运动状态、密度比例因子估计值以及轨迹跟踪误差都保持在由传感器噪声水平、dsf滤波器带宽以及mpc预测范围所决定的有界区域内。第5.3节中进行的蒙特卡洛仿真进一步验证了这一理论结论：在采用所提出架构的500次模拟中，所有结果都处于可恢复的飞行范围内，且终端误差几乎都集中在零点附近（平均高度误差为-63.8米，标准差为346.9米）。总而言之，在给定的假设条件以及特定的进入窗口范围内，各个子系统都会对整个系统产生有界的误差传递效应，而整个复合架构也确实具备上述定义的实际稳定性。5. 数值仿真本节旨在评估在严重大气不确定性条件下的所提出的ekf-dsf mpc架构的性能。所有实验都有意将闭环系统置于超出其常规工作范围的极端条件下，以此来考察系统的稳健性、故障机制以及估计器的精确度。实验的主要目标有两个：其一，验证解耦式的估计组件能够在初始条件存在极大不确定性的情况下，依然能够可靠地恢复状态和密度参数；其二，证明将mpc与估计器相结合的完整在线系统，确实能够高精度地跟踪参考轨迹，从而避免在传统控制方案中出现的性能下降现象。5.1. 仿真设置与参数数值实验是在一台配备有intel core i5-12500h处理器（基础频率2.5吉赫兹，拥有12个性能核心）以及16gb ddr5-4800速率内存的笔记本电脑上进行的，实验过程分为两个与所提出架构相对应的不同计算阶段。在离线的稳健规划阶段，第3节中提到的时变最优控制问题被转化为一个非线性规划问题，这一转换是通过通用的最优控制软件gpops-ii [48]来实现的。之后，这个大规模的稀疏非线性规划问题又通过顺序二次规划求解器snopt [49]来处理，该求解器采用主动集方法来严格约束路径约束条件。进入在线实时跟踪阶段后，则通过高精度的积分步骤来模拟多速率的估计与控制循环（见第4节）。对于慢速循环的mpc部分，其凸形二次规划问题是通过yalmip软件（版本r20230622）[50]来建模的，随后则由mosek内点优化器（版本10.1）[51]来求解。该二次规划求解器设定的相对容差为10^-6，且每进行一次mpc迭代最多允许200次内点运算。 ρnom(h(t))表示真实的、随时间变化的密度比例因子，而k?ρ(k)=K?ρ(k)?Ktrue(tk)则表示DSF的估计误差。DSF滤波器实际上是一个一阶线性时不变低通滤波器，其增益为标量Gdsf∈(0,1)。该滤波器的极点位于z=Gdsf处，且严格处于单位圆内部，因此能够保证有界输入下的有界输出稳定性。由于EKF的估计误差是均方有界的（见命题2），加之IMU加速度测量值的噪声方差也是有界的，因此原始的密度比例因子估计值Kest(k)本质上是一个有界的随机序列。经过过滤后的估计值会以指数级收敛到真实密度比例因子的附近区域：(42)??|K?ρ(k)?Ktrue(tk)|≤Gdsfk|K?ρ(0)|+11?Gdsf(?snr+?drift)，其中??snr受传感器噪声水平以及瞬时的空气动力信噪比的影响，而??drift则反映了由于真实密度比例因子Ktrue(t)在滤波器窗口内的缓慢变化所导致的跟踪滞后现象。这一特性可以解释图10(a)中所观察到的“宽-窄-宽”方差分布特征：在高空动态压力较大的飞行阶段，强烈的阻力信号会导致较低的??snr值，从而使置信区间较为狭窄；而在高空及末端阶段，由于空气动力作用较弱，方差则会有所增大。值得关注的是，整个过程中的偏差始终保持在零点附近，这正符合有界输入有界输出稳定性特性。接下来分析外部控制循环。需要明确的是，这里的讨论仅针对有限进入窗口条件下的局部实际稳定性问题，并非针对整个高超音速动力学系统的全局非线性稳定性证明。设δxk=xk?xref(tk)表示跟踪误差。降阶后的MPC模型（公式(37)）实际上是一个带有终端状态惩罚项P的有限时间范围跟踪问题。在无约束情况下，其解与有限时间范围的LQR问题一致；如果P满足离散代数Riccati方程，那么最优成本函数V?(δx)=‖δx‖P2就可以作为闭环线性化系统的李亚普诺夫函数，从而保证跟踪误差的局部渐近稳定性[5]，[45]。在存在控制约束的情况下，此处采用的终端惩罚项P≈Q实际上是对无限时间范围成本函数的二次近似。虽然这种选择并不能构成原始非线性高超音速动力学系统的严格控制李亚普诺夫函数，但它能够提供足够的终端成本，避免因有限时间范围的截断而导致控制策略出现短视或不稳定现象[5]。此外，由于在每个MPC迭代步骤中，Ak、Bk这两个矩阵都是通过RK4积分器的中心差分法重新计算的，因此线性预测模型能够准确捕捉沿参考轨迹的局部动态行为，从而减少可能降低系统稳定性的线性化误差。整个闭环系统实际上形成了一个级联结构：EKF和DSF估计器负责提供状态和密度估计值给MPC，而MPC则根据这些信息生成控制指令。由于估计循环的运行速度比控制循环快一个数量级（Tf=0.1秒 vs. tmpc=1.0秒），因此由初始化或大气条件变化所引发的暂时性估计误差会在连续的MPC迭代之间迅速衰减。从MPC的角度来看，估计残差x?k和K?ρ(k)被视为有界的附加干扰。关于跟踪MPC的现有稳健性研究[5]表明，有界的输入干扰只会导致状态偏离名义轨迹一定的范围；这一特性被正式定义为输入到状态的实际稳定性ISpS[47]。因此，整个复合系统具备实际稳定性：运动状态、密度比例因子估计值以及轨迹跟踪误差都保持在由传感器噪声水平、DSF滤波器带宽以及MPC预测范围所决定的有界区域内。第5.3节中进行的蒙特卡洛仿真进一步验证了这一理论结论：在采用所提出架构的500次模拟中，所有结果都处于可恢复的飞行范围内，且终端误差几乎都集中在零点附近（平均高度误差为-63.8米，标准差为346.9米）。总而言之，在给定的假设条件以及特定的进入窗口范围内，各个子系统都会对整个系统产生有界的误差传递效应，而整个复合架构也确实具备上述定义的实际稳定性。5. 数值仿真本节旨在评估在严重大气不确定性条件下的所提出的ekf-dsf mpc架构的性能。所有实验都有意将闭环系统置于超出其常规工作范围的极端条件下，以此来考察系统的稳健性、故障机制以及估计器的精确度。实验的主要目标有两个：其一，验证解耦式的估计组件能够在初始条件存在极大不确定性的情况下，依然能够可靠地恢复状态和密度参数；其二，证明将mpc与估计器相结合的完整在线系统，确实能够高精度地跟踪参考轨迹，从而避免在传统控制方案中出现的性能下降现象。5.1. 仿真设置与参数数值实验是在一台配备有intel core i5-12500h处理器（基础频率2.5吉赫兹，拥有12个性能核心）以及16gb ddr5-4800速率内存的笔记本电脑上进行的，实验过程分为两个与所提出架构相对应的不同计算阶段。在离线的稳健规划阶段，第3节中提到的时变最优控制问题被转化为一个非线性规划问题，这一转换是通过通用的最优控制软件gpops-ii [48]来实现的。之后，这个大规模的稀疏非线性规划问题又通过顺序二次规划求解器snopt>该配置确保优化过程在1.0秒的设定时间范围内顺利收敛。停滞点加热速率是通过修正了攻角相关修正因子的广义查普曼方程来建模的，如[32]中所描述：(43) q?=kQ?Cq(α)?ρn1/Vn2，其中Cq(α)=Cq0+Cq1α+Cq2α2+Cq3α3用于描述加热系数随飞行器姿态的变化。具体的系数值以及其他飞行器物理参数详见表1。为严格评估系统的鲁棒性，研究中引入了显著的随机不确定性。依据[26]中的随机不确定性建模方法，假设真实的大气参数向量遵循联合高斯分布。大气协方差Σp被构建为对角矩阵diag(σρ02,σβ2)，前提是各参数之间不存在先验相关性。包括这些大气波动、初始状态误差以及传感器噪声特性在内的所有仿真参数汇总于表2中。为评估在这些不确定性条件下的性能，定义了三种主要测试场景：表1展示了飞行器物理参数及任务边界条件。参数符号及其SI单位数值如下：飞行器属性方面，质量m为102,285千克；参考面积Sref为249.9平方米；升力系数多项式CL0…3的取值为[?0.207,1.676,0,0]；阻力系数多项式CD0…3的取值为[0.079,?0.353,2.040,0]。气动热模型方面，基础热系数kQ为1.066×10?4（SI单位）；加热率指数n1、n2分别为0.5和3.07；攻角相关加热系数多项式Cq0…3的取值为[1.067,?1.101,0.699,?0.190]；规划阶段的加热率限制为q?plan=794.5千瓦/平方米；结构层面的加热率限制为q?limit=1.23×q?plan。初始条件方面，t0时刻的高度h0为79,248米，速度V0为7334米/秒，飞行路径角度γ0为?1.06°。目标条件方面，tf时刻的高度hf为24,384米，速度Vf为762米/秒，飞行路径角度γf为?5.0°。控制约束方面，攻角α的取值范围为[αmin,αmax]，即[?90.0°,+90.0°]；侧倾角σ的取值范围为[σmin,σmax]，即[?90.0°,+1.0°]。表2列出了仿真参数及其不确定性分布（3σ值）。参数符号、对应数值或分散度（3σ值）如下：定时方面，滤波器步长Tf为0.1秒，MPC步长Tmpc为1.0秒；初始状态误差方面，径向位置δr的误差范围为±1000米，速度δV的误差范围为±3.3米/秒，飞行路径角度δγ的误差范围为±0.1°；大气不确定性方面，平均密度ρ0的均值μ为1.225千克/立方米，平均尺度高度β的均值μ为1.378×10?4米^-1，密度的分散度（3σ）为σρ0±15%×ρ0、μ，β的分散度（3σ）为σβ±5%×βμ；传感器噪声方面，GPS位置的噪声标准差wr为10米，GPS速度的噪声标准差wV为1.0米/秒，IMU加速度的噪声标准差wa为0.01米/秒^2。估计器调优方面，EKF的初始协方差矩阵P0为diag(102平方米,12米/秒^2,12度^2)；DSF的增益Gdsf为0.2（低通滤波器）。三种测试场景如下：1.场景A（“朴素”基准）：该场景下MPC运行时不使用DSF估计器，仅依赖离线版本的标称大气模型（p=pμ）进行预测，实际上相当于假设密度尺度因子Kρ≡1.0，代表了不进行在线参数调整时的基准性能。2.场景B（“自适应”方案）：所提出的EKF-DSF架构处于完全激活状态，MPC会利用DSF估算出的实时密度尺度因子K?ρ(t)来修正其内部模型。3.场景C（“理想”极限）：MPC能够获得系统实际使用到的真实大气参数（ptrue），这一场景可作为理论上的性能上限，用于量化由估计误差带来的性能损失。4.场景D（“对比”基准）：该场景基于[10]中的状态空间观测器方法，采用扰动观测器辅助的MPC进行对比。扰动观测器通过三个等效的扰动模式η=[ηD,ηL,ηρ]T来补充运动学状态，这些扰动模式遵循一阶高斯-马尔可夫过程演变。气动力被建模为D=(1+ηρ+ηD)Dnom，L=(1+ηρ+ηL)Lnom，观测器则通过实测与预测气动加速度之间的差异来修正这些扰动模式。这些扰动估计值会被作为力级修正量注入MPC预测模型中。观测器的过程噪声及初始协方差经过保守调整，以在抑制扰动影响与减少噪声放大之间取得平衡，同时扰动幅度也被限制在物理上合理的范围内。为确保比较的公平性，四种场景共享相同的条件：相同的500组大气参数样本，相同的初始状态误差范围，相同的控制约束，相同的MPC预测时间范围（N=15）及成本权重矩阵，以及相同的传感器噪声特性。各场景的区别在于估计策略：场景B采用显式的密度剖面重建方式，而场景D则采用等效的力级扰动补偿方式。5.2 极端大气条件下的估计器验证为严格测试系统在无需考虑物理上不太可能出现的组合情况（例如多个无关参数同时出现3σ偏差）下的性能极限，测试案例是根据3σ联合置信椭圆[52]来设定的：(44) δρ0/σρ02 + δβ/σβ2 = 3。在联合置信椭圆与其对角线的交点处选取了四个战略性的测试点，对应于±3σ/2的偏差范围，这一范围大约涵盖了98.9%的所有可能大气条件。这些场景包括了两种极端的物理状况：一种是“稀薄/罕见”的大气条件（情况1：?3σρ0/2,+3σβ/2），另一种是“密集/厚重”的大气条件（情况2：+3σρ0/2,?3σβ/2），此外还有两种混合条件（情况3：?/?和情况4：+/+），以此确保全面覆盖参数不确定性范围。图4展示了仅使用GPS的EKF的估计误差情况。尽管初始时刻的初始化误差很大（δr0=1000米，δV0=3.3米/秒，δγ0=0.1度），但该滤波器仍在前50秒内实现了快速收敛。关键的是，稳态误差保持在可控范围内，且与极端大气条件无关，这证明了分离式设计的有效性：通过将EKF与气动模型分开，即便密度模型存在严重误差，导航结果依然保持无偏。下载：下载高分辨率图像（517KB）下载：下载完整尺寸图像图4. 在3σ初始误差条件下，仅使用GPS的EKF状态估计误差的时间变化曲线。图5(a)展示了DSF即Kρ(t)的跟踪性能，该估计器能够正确识别各种场景下的物理状况。例如在密集大气条件下的情况2中，DSF正确识别出大于1.0的较大尺度因子Kρ，以此解释额外的阻力现象。估算得到的K?ρ与理论上的真实比值Ktrue相当，这表明低通滤波器能够有效捕捉由β参数不匹配引起的时变偏差。下载：下载高分辨率图像（491KB）下载：下载完整尺寸图像图5. 大气不确定性条件下的DSF/MPC性能对比。最关键的测试是在大气不确定性条件下的闭环跟踪能力。图5(b)比较了在三种不同场景下，朴素MPC（场景A）、所提出的自适应EKF-DSF MPC（场景B）以及扰动观测器辅助的MPC对比方案（场景D）相对于参考轨迹的高度跟踪误差情况。图5(b)中的红色虚线显示了忽视参数不确定性所带来的严重后果，尤其是如敏感性分析所示（见引理1），β参数偏差会导致误差呈指数级放大。如图5(a)所示，在情况1和情况4这类稀薄大气条件下，实时密度尺度因子Kρ会低于1，意味着大气密度低于正常水平。在这种低密度环境下，朴素MPC会高估阻力和升力，导致早期无法有效消耗能量，进而出现持续下降的高度，最终在到达目标之前就导致飞行器坠毁（h<0）。而在密集大气条件下的情况2和情况3中，MPC则会因为过大的升力而导致飞行器在中期出现高度上升现象，随后又因持续的高阻力而耗尽能量，最终导致终端阶段高度不足。与基准方案的失败相比，自适应EKF-DSF MPC（蓝色实线）通过利用实时估算的K?ρ值动态调整侧倾角，成功避免了这些任务失败的情况。在要求极高的密集大气条件下的情况2中，朴素MPC根本无法维持飞行高度，而自适应MPC则能够完成整个飞行轨迹，体现了在极限状况下优先保证飞行器存活而非追求最终精度的设计思路。尽管其终端高度误差达到了2.67公里，但这一误差是由于飞行器升力能力达到物理极限所致，而且飞行器始终处于可恢复的飞行范围内。另外，在稀薄大气条件下的情况1中，自适应控制器也能正确识别出低密度环境，最终实现的终端高度误差仅为?64.81米，速度误差为?28.56米/秒。场景D（绿色点划线）表明，扰动观测器辅助的MPC也能从在线补偿中受益，其性能明显优于朴素控制器，不过其在四种不同大气条件下的性能表现并不一致。由于扰动观测器的补偿是在等效的阻力/升力层面进行的，而非通过显式重建密度剖面来实现，因此局部的气动残差并不一定能够反映MPC预测范围内的整体密度偏差。所以，虽然场景D在有利条件下也能达到类似的精度，但其最严重的终端高度误差却高达约?9.14公里，这一数值明显大于所提出的EKF-DSF MPC的误差水平。总体而言，尽管不同极端条件下的终端误差存在差异（高度误差范围从64.8米到2665米，速度误差范围从3.05米/秒到28.56米/秒），但该系统在整个不确定性范围内始终保持着稳定性。通过避免基准方案中出现的任务失败情况，所提出的架构成功地将极端不确定性转化为可控的跟踪误差。5.3 蒙特卡洛鲁棒性分析虽然前一节中的单案例分析已经展示了所提出架构的工作原理，但要想获得更具代表性的验证结果，还需要通过统计分析来进行。为此，研究人员开展了大规模的蒙特卡洛模拟，共进行了Nmc=500次模拟运行。在介绍跟踪统计结果之前，先概述一下所提出的EKF-DSF MPC循环的在线计算成本。在模拟过程中，每次MPC更新所需的实际时间都是在调用在线控制程序时记录下来的。所记录的时间包括时间范围线性化、数据整理以及QP求解器的运行时间，但不包括植物模型的传播、测量数据的生成、图表绘制以及控制器的一次性构建时间。表3给出了三种不同预测时间范围内的平均计时统计结果。第二行数据既包含了整个飞行过程中的慢周期MPC更新次数，也包含了快周期的EKF-DSF估计更新次数，所有数值都是以1789.8秒的标称飞行时间为基准进行归一化的。在这三种时间范围设置下，所有模拟均顺利完成在线计时评估，没有出现求解器故障或计时违规的情况。正如预期，MPC的平均更新时间会随着预测时间范围的延长而增加，但总的在线计算时间仍然仅占实际飞行时间的1%左右。对于所采用的N=15的时间范围，MPC的平均更新时间为7.34毫秒，最大更新时间为100.95毫秒，这两个数值都远低于Tmpc=1.0秒的设定时间要求，即便在最不利的情况下也有充足的实时冗余。每次EKF更新以及DSF滤波器的计算时间都不到0.1毫秒，对整体计算时间的影响可以忽略不计。表3. 不同预测时间范围内所提出的EKF-DSF MPC循环的计算时间。指标N=15时N=20时N=25时MPC平均更新时间（毫秒）7.349.2410.90第95百分位数的MPC更新时间（毫秒）9.2911.7013.80MPC最大更新时间（毫秒）100.95128.08149.91在线计算时间/飞行时间百分比0.800.991.15需要指出的是，这些计时结果并非作为与硬件无关的下限值。所测得的计算速度是基于第5节中提到的i5-12500H处理器所得出的。所提出的公式中的在线优化问题属于带约束的QP问题，若使用更高效的专用求解器、预热策略、固定点或主动集实现方法，或是用C/C++语言进行底层实现，都可以进一步缩短计算时间。因此，当前的结果应被视为在高层次研究实现条件下的实际可行性证明，而非该架构的理论计算极限。任务失败是指在飞行过程中出现以下任意一种情况：(1)高度h<0（撞击地面）；(2)在到达目标之前速度降至50米/秒以下；或者(3)瞬时加热率q?(t)>q?limit。相反，任务成功则是指能够以Vf≥50米/秒且q?(t)≤q?limit的条件进入终端阶段。在每次模拟运行中，都会引入双重不确定性。首先，真实参数ptrue是从联合高斯分布N(pμ,Σp)中随机抽取的；其次，飞行器的初始状态x0会同时受到高斯导航误差的影响，这些误差的标准差是根据表2中的3σ限值确定的。这样的设置能够模拟任务中可能遇到的极端操作条件。上文提到的四种场景会被相互比较，为了进一步体现所提出架构相较于现有的在线扰动观测器辅助导航策略的优势，还会将场景B与场景D进行对比。图6展示了这两种方案下终端高度误差和速度误差的概率密度函数（PDF）直方图。扰动观测器辅助的场景D虽然并非刻意设计的性能较弱的控制器，但它仍能提供良好的对比基准——其任务成功率为89.2%，与场景B的90.6%相当。不过，其终端误差的分布范围更广。如表4所示，以及图6(a)中的高度误差PDF曲线所示，场景D的平均高度误差为?270.1米，标准差为1192.0米，而场景B的高度误差则更为集中，平均值为?63.8米，标准差为346.9米。尽管场景D在原点附近也出现了许多误差较小的情况，但其较大的误差尾部表明大误差出现的频率更高。图6(b)中的终端速度PDF曲线也呈现出类似的趋势，场景D的速度标准差为72.7米/秒，而所提出架构对应的速度标准差仅为9.2米/秒。这说明，虽然等效的扰动补偿能够提升系统的基本鲁棒性，但它无法像精确的密度剖面重建那样，在整个MPC预测时间内保持一致的密度分布。下载：下载高分辨率图像（294KB）下载：下载完整尺寸图像图6. 蒙特卡洛误差分布：场景D（对比组）与场景B（自适应组）。下载：下载高分辨率图像（244KB）下载：下载完整尺寸图像图7. 蒙特卡洛误差分布：场景B（自适应组）与场景C（理想组）。在确认了系统的鲁棒性之后，接下来通过将自适应场景B与完美信息场景C进行对比，来评估估计器的最优性。图7展示了这两种场景的误差分布情况，两者几乎完全重合，这一结论在图8(b)中也有清晰体现。如表4所示，估算解的速度标准差为9.2米/秒，与完美信息基准值8.0米/秒相当。这说明，与由系统动态和离散控制引起的固有引导误差相比，EKF-DSF架构引入的估算误差可以忽略不计。对场景B（9.4%）、场景C（10.6%）和场景D（10.8%）中的热故障分析表明存在根本的“跟踪与安全冲突”。在密度较大的大气中，物理上最优的路径自然需要更强的减速。然而，MPC的任务是跟踪平均参考值xref。当飞行器处于密度更高的流体中（ρ>ρref）时，控制器仍需迫使其保持较高的名义速度（V→Vref），这就不得不承受会消耗1.23安全裕度的临时热负荷。场景D中的热故障次数略少，但其终端能量误差较大，这表明其跟踪行为较为松散，并非精度更高。这些结果证实，故障是由刚性跟踪的物理特性导致的，而非估算误差。场景B的9.4%在线热故障率仍低于通过Cantelli不等式（方程(25)）离线得出的10%最坏情况阈值，这表明在线跟踪动态并未从根本上否定离线概率安全裕度。所提出的架构成功地将致命碰撞转化为可控制的临时负荷，这些负荷大多仍处于既定的工程安全缓冲范围内。为进一步评估自适应架构与理想基准之间的性能差距是否显著，进行了500次蒙特卡洛模拟，并采用了配对Wilcoxon符号秩检验[53]。最终得到的终端高度误差的p值为0.4848，速度误差的p值为0.5542，均远高于0.05的显著性阈值。因此，无法拒绝两种分布相同的零假设，这意味着在Wilcoxon符号秩检验下，性能差距并无统计显著性。这些结果表明，所提出的架构能够有效恢复具备完美大气知识系统的性能。此外，为可视化系统在不确定性条件下的动态演变，图8(a)对比了传统MPC（红色）与所提出的自适应MPC（蓝色）的高度和速度轨迹。表4为蒙特卡洛统计结果（N=500）。架构成功率高度误差（米）速度误差（米/秒）空值比率均值标准差（σ）均值标准差（σ）场景A：传统MPC58.8%?4470.96102.1?450.1267.4场景B：EKF-DSF MPC90.6%?63.8346.9?4.39.2场景C：真实MPC89.4%?56.2353.0?4.38.0场景D：DOB辅助MPC89.2%?270.11192.0?14.772.7下载：下载高分辨率图像（710KB）下载：下载全尺寸图像图8：蒙特卡洛轨迹簇（N=500）。下载：下载高分辨率图像（428KB）下载：下载全尺寸图像图9：两个典型故障场景的控制与状态历史详细对比。下载：下载高分辨率图像（314KB）下载：下载全尺寸图像图10：Kρ的统计数据（N=500）。场景A的轨迹表现出明显的扩散现象。由于没有参数自适应，MPC依赖于错误的密度模型。在低密度区域，由于气动制动不足，飞行器保持过高的速度，从而导致在到达目标之前就过度下坠并发生碰撞（h<0）。在高密度区域，由于过大的阻力，飞行器过度减速，出现中期高度超调，随后能量耗尽，最终无法精确到达目标。相比之下，场景B的轨迹紧密围绕参考值（黑色）分布。尽管存在随机的大气变化，EKF-DSF架构仍能将飞行器保持在接近参考值的轨道上，有效抵消了未知参数的影响。为了解释为何传统控制器会失效，需研究图9中所示的两个典型碰撞场景的内部控制动态。在图9(a)所示的早期故障案例中，飞行器遇到了密度较低的大气环境（Kρ<1）。传统MPC没有意识到这一差异，初期纠正力度不足，导致高度误差逐渐累积。大约在t=600秒时，飞行器过早地进入密度更高的低层大气，从而受到突然的力矩冲击。为应对这一问题，优化器在t=800?900秒左右指令将俯仰角σ迅速提升至上限，最终在t=1600秒左右将攻角α提升至其物理极限值。这种“饱和”状态意味着控制能力完全丧失：由于飞行器的配置限制，它无法应对突然增加的力矩，从而导致在t≈850秒时速度急剧下降并发生碰撞。在图9(b)所示的晚期故障案例中，故障发展更为缓慢，飞行器面临的是密度高于名义模型的大气环境（Kρ>1）。传统MPC不断通过调整俯仰角σ来对抗过大的阻力，大约在t=1000秒时逐渐将其提升至上限，以产生更多垂直向上的升力。虽然这样能使飞行时间更长，但轨迹在第一个热负荷峰值之后会出现不稳定的振荡，反弹高度超过参考值，后续下降速度虽慢但更为明显。未经校正的能量耗尽最终导致攻角α在波动中下降，然后在t=1600秒左右趋于饱和，从而导致飞行器突然急速下降，最终偏离目标超过10公里。至关重要的是，自适应MPC（蓝线）完全避免了这些异常模式。它通过实时估算密度比例因子K?ρ，能够及早发现密度不匹配问题。它还会主动调整俯仰角曲线以保存能量，从而保持平稳的控制裕度，确保飞行器安全抵达目标位置。最后，为进一步验证所提架构的内部可靠性，分析了密度比例因子估算误差的统计特性，该误差定义为eK(t)=K?ρ(t)?Ktrue(t)。图10(a)展示了所有蒙特卡洛模拟中均值误差（偏差）及3σ置信区间（经过平滑处理）的时间演变情况，结果表明，估算精度主要受空气动力学测量值的时变信噪比影响。方差呈现出“宽-窄-宽”的规律：在高空阶段（0–900秒）和终端阶段（>1600秒），由于大气较稀薄或速度较低，阻力信号较弱，接近传感器噪声底限，此时置信区间较宽；而在轨迹中间阶段（900–1600秒），较高的动压会产生较强的阻力信号，使得估算器能够快速收敛，方差较小。尽管瞬时信号质量存在波动，但整个轨迹过程中的均值误差（偏差）始终在统计上接近零，这表明DSF低通滤波器能够有效滤除传感器噪声，而不引入系统性偏差。此外，500次蒙特卡洛模拟的统计性能表现为总均方根误差为0.0085，均值偏差仅为0.00012，这充分证明了该架构能够有效将大气参数与随机传感器噪声分离。作为时间历史分析的补充，图10(b)展示了最终估算误差的概率分布。直方图呈现出严格以零为中心的高斯分布特征，这进一步证实了估算器的无偏性。误差残差紧密聚集在±0.004范围内，这体现了系统的高精度。6. 结论本文研究了在初始状态误差与极端大气密度不匹配的共同不确定性条件下，自主高超音速再入的问题。理论分析表明，传统引导方案之所以失败，是由于雅可比矩阵的共线性导致耦合估算存在根本性的病态问题，而非控制能力不足。为克服这些限制，本文提出了一种解耦的估算-控制架构。内环将与大气无关的运动状态估算与环境参数识别分开，从而能够可靠地获得有效的密度修正值。外环则采用责任分离策略，通过离线风险约束优化预先设定热安全裕度。这使得可以采用计算效率更高的降阶跟踪MPC，从而无需让机载计算机承担主动的热约束线性化工作。大量的蒙特卡洛模拟验证了该方法的有效性。所提出的框架避免了传统控制器中出现的巨大轨迹偏差，其跟踪精度与具备完美大气知识的理论系统相比并无统计上的显著差异。这些结果表明，对于具有状态与环境不确定性高度耦合的高非线性高超音速系统，通过在估算、规划与控制之间实现架构解耦，而非仅仅依赖越来越复杂的在线优化，有助于提升系统的自主性。未来的工作将整合预测性横向引导技术，并研究针对非指数型大气结构的数据驱动残差学习方法。作者贡献说明李宇堂：撰写原始稿件。李松洲：监督指导。周迪：资金筹集。何震：资源提供。李宇堂|李松洲|周迪|何震中国黑龙江省哈尔滨市哈尔滨工业大学航天学院，邮编150000