追逃博弈在航天领域并不新鲜但把“不完全信息”和“Epsilon纳什均衡”组合进末端追逃问题再叠加上EKF在线参数估计就值得认真复现一次了。这篇文章我想完整拆解一下这一类基于 Matlab 的复现工作到底在做什么从建模、EKF的“引线”作用、均衡求解到代码实现和调试时踩过的坑一次性讲透希望能给正在搞航天制导控制、空间博弈算法或者想复现类似期刊工作的朋友提供一些参考。很多刚接触这个方向的人会被“Epsilon纳什均衡”这几个字吓到觉得理论门槛很高。实际上我们在末端追逃场景里提到的Epsilon纳什均衡落地到算法层面就是一句话受限于在线估计精度和模型简化博弈双方做不到绝对的、精准的最优响应只能做到“误差有界”的最优响应这个误差界的宽窄就对应Epsilon。理解了这层物理含义后续推导和Matlab代码就会顺很多。本文涉及的内容包括相对运动建模、扩展卡尔曼滤波EKF对对手参数的在线估计、基于估计值求解近似纳什均衡策略以及整套流程的Matlab仿真实测。1. 末端追逃博弈的第一个问题动力学模型的尺度选择1.1 为什么博弈框架在这个尺度上成立末端追逃的典型场景是两个航天器在最后几十公里甚至几公里的相对距离上进行接近与规避对抗。追击方想最小化末端脱靶量逃逸方想最大化末端脱靶量双方同时决策、同时机动这天然是一个双方零和微分博弈问题。和传统单边最优制导最大的区别在于博弈框架下追击方的策略必须考虑逃逸方的“理性反击”而不是简单地把对手当作固定目标或已知轨迹。有人会问直接用鲁棒控制或者滑模制导不也能处理外界扰动吗从工程结果上看某些鲁棒制导的确能扛住一定程度的机动不确定性但鲁棒控制不关心对手是否也在“优化自己的收益”它只是把所有不确定性当作最坏情况的干扰。而博弈策略的核心优势是能够对齐双方的“决策交互”对手机动策略变化时追击方的策略可以自适应地跟着调整。在EKF估计的配合下这种调整不仅是基于状态的反馈还是基于对手“底牌”的反馈。1.2 相对运动方程选择C-W方程还是非线性轨道方程末端追逃问题里最常见的做法是选用C-W方程Hill方程描述相对运动因为它形式简单、解析性质好方便构造线性二次型博弈。线性化的C-W方程如下相对位置向量 r [x; y; z]相对速度向量 v [vx; vy; vz]在圆轨道参考系下相对运动近似满足 \dot r v\dot v_x 3 n^2 x 2 n vy u_x d_x\dot v_y -2 n vx u_y d_y\dot v_z -n^2 z u_z d_z其中 n 表示参考轨道角速度u 为控制加速度d 为机动加速度扰动如果逃逸方的机动能力用有界控制项表示整个方程可以整理成状态空间形式\dot X A X B_p u_p B_e u_e。A和B矩阵的具体形式和参考轨道参数有关复现时只要把轨道高度、偏心率这些参数代入即可。但要注意C-W方程只适用于近距离、短时间、近圆轨道相对运动。如果初始相对距离超过几十公里或者双方机动加速度很大线性误差会明显影响结果的合理性。我在复现时对比过线性模型和非线性轨道相对运动模型的差异发现当末端时间在30秒以内、相对距离小于10公里时C-W方程和数值轨道模型的相对位置误差能控制在5%以内这个精度对博弈策略验证是够用的但如果末端时间拉到上百秒线性简化带来的偏差会导致末端脱靶量出现“假性收敛”也就是仿真看着效果不错实际场景中却跑偏。所以复现这类期刊工作时建议先明确末端时间和初始距离的量级再决定模型复杂度。对于末端追逃这种“末端”属性很强的问题短时间、近距离场景下用C-W方程完全合理。1.3 不完全信息的数学描述把“不知道”翻译成参数化形式不完全信息在博弈论里有很多种表达方式在追逃问题里最常见的是对手的成本函数参数、机动能力上限、或者终端偏好权重未知。举个具体的例子追击方知道逃逸方会最小化自己被捕获的概率但不知道逃逸方对控制能量消耗的“心疼程度”——是愿意不惜一切代价机动还是倾向于省燃料只做小幅规避。这个偏好权重就是逃逸方的“底牌”。在数学上我们把逃逸方的控制加权矩阵 R_e 写成已知结构 未知参数的形式比如 R_e(θ) diag([θ1, θ2, θ3])θ 就是要估计的未知参数向量。这样不完全信息问题就转换为一个参数辨识问题。此时“不完全信息”不再是抽象意义上的“敌方意图”而是具体到某些能被EKF实时估计的参数。在这个转换中一个很重要的前提是要保证系统可辨识性。也就是说需要确认测量数据能够把 θ 和状态量区分开。复现时我专门看了量测矩阵对扩展状态的敏感性如果只有视线角测量θ 的估计收敛就比较慢如果加上相对距离和相对速度测量θ 的收敛速度会明显加快。这也是很多期刊复现代码里把测量模型设置为“相对位置相对速度”的原因。2. EKF参数估计把对手的“底牌”从数据里滤出来2.1 为什么纳什均衡需要参数估计如果假设完全信息博弈双方的参数都已知那么整个问题就是一个标准的二次型微分博弈可以直接求解Riccati微分方程得到反馈控制策略。但现实里追击方并不掌握逃逸方的真实偏好参数那么计算出的“最优”策略很可能不是真正针对对手的最优响应甚至可能被对手的机动方式牵着走。EKF在这里起的作用就是“在线采参数”通过量测数据估计逃逸方的未知参数随着滤波收敛把原本不确定的模型逐步逼近真实模型然后再基于这个逼近后的模型求解博弈控制。这样一来整个闭环系统就有了“感知-估计-决策-执行”的完整链条也就是标题里“自适应博弈”的含义。有朋友可能会问直接用强跟踪滤波STF或者无迹卡尔曼UKF是否更好对这个问题我自己的体会是EKF在计算效率和实现简洁性上最平衡尤其当你已经把状态方程线性化模板写好之后切换到不同场景只需要改雅可比矩阵。而UKF虽然对强非线性更友好但Sigma点传播会引入额外计算量在机器人或航天嵌入式平台上不太划算。期刊复现里大多以EKF为主也是因为它的计算量小方便快速验证博弈算法的核心逻辑。2.2 扩展状态构造与EKF预测更新流程构造扩展状态向量时我们把原来的相对运动状态和未知参数拼接在一起X_aug [x; y; z; vx; vy; vz; θ1; θ2; θ3]对应的状态方程为前6个状态遵循C-W方程后3个参数状态取 \dot θ 0表示在滤波时间尺度内对手偏好不突变量测方程取Z H X_aug v_noise其中H为线性量测矩阵v_noise为零均值高斯白噪声。因为EKF的状态方程本身带有控制输入项而这个控制输入项在估计时是已知的追击方自己发出的控制指令所以EKF预测阶段要把 u_p 代进去。逃逸方的控制则依赖当前估计参数在滤波里用一个“当前最优猜测”近似形成“EKF估计-博弈求解-实际执行”的循环。EKF每个滤波周期的过程大致是计算雅可比矩阵 F ∂f/∂X_augΦ expm(F * T_sample)预测协方差 P_pred Φ * P * Φ Q计算卡尔曼增益 K P_pred * H * (H * P_pred * H R)^{-1}状态更新 X_new X_pred K * (Z_real - Z_pred)协方差更新 P_new (I - K * H) * P_pred这里最容易出事的是雅可比矩阵的推导。如果状态方程里同时包含追击方控制 u_p 和逃逸方估计控制 \hat u_e那么F要对状态求偏导而 \hat u_e 本身又是估计参数 θ 的函数这一步需要把链式法则写对。我自己在复现时就在这里犯过错误F矩阵少了一项导致EKF在末端十几秒时估计值突然发散。后来把F矩阵的每个元素都做了符号求导验证才定位到问题。2.3 可观性分析与测量配置的工程考量EKF能不能收敛很大程度上取决于扩展状态的可观性。我们做了简化推导如果量测只是相对位置的三个分量那么系统对 θ 的可观性矩阵往往不满秩尤其当相对加速度方向不变时θ的某些分量无法被有效激励。实际问题里逃逸方机动通常是在末端才开始明显变化所以前段飞行时 θ 估计波动大末端机动激发后才逐渐收敛。工程上的对策是增加量测维度把相对速度甚至视线角速率引入量测。这样虽然EKF的矩阵维度变大但可观性明显改善。我在代码里把量测设置为相对速度和相对位置同时测量对比实验表明 θ 估计的均方根误差从20%以上降到了8%左右。另外有一个细节EKF的噪声矩阵R不能设得太小否则滤波增益过高量测噪声会被放大到状态估计中导致博弈策略抖动也不能设得太大否则EKF跟踪不上末端快速变化的相对运动。一个可行的调参方案是从名义值出发通过蒙特卡洛仿真对比估计误差和末端脱靶量的趋势选定一个在中段和末端都能接受的折中值。3. Epsilon纳什均衡在不完美估计下定义“次优但够用”的策略3.1 不完全信息均衡的定义回顾纳什均衡的核心思想是没有任何一个参与者可以单方面改变策略获得更好收益。完全信息下这个定义是严格的但在不完全信息下由于存在估计误差策略不是建立在“真实参数”上而是建立在“估计参数”上因此很难达到精确的纳什均衡。Epsilon纳什均衡的定义很简洁对每个参与者如果其策略相对于任意单边偏离所能带来的收益增益不超过一个很小的正数ε则称该策略组合构成ε-纳什均衡。放到追逃里就是说追击方和逃逸方各自按照估计模型采用近似最优策略双方都不存在能使自己多获得超过ε的收益的“单方面聪明改变”。这个定义既保留了纳什均衡的核心精神又给在线估计和模型简化留了余地。复现中判断ε是否可接受通常直接看性能指标差异拿估计策略和真实完全信息策略分别计算末端脱靶量如果差异足够小就认为ε在工程可接受范围内。这个方法虽然有点“事后验证”的味道但对仿真验证来说是最直观有效的。3.2 末端追逃博弈策略的求解路径在短时间末端场景下线性二次型微分博弈是最常用的可解析求解框架。目标函数一般写成J 1/2 * x(T) * Q_T * x(T) 1/2 * ∫_0^T [u_p * R_p * u_p - γ^2 * u_e * R_e(θ) * u_e] dt这里追击方最小化J逃逸方最大化J。Q_T是末端相对位置误差的惩罚矩阵R_p和R_e是控制能量惩罚。γ是逃逸方的“能量上限系数”反映了逃逸方的机动能力边界。对于这个目标函数标准结果是鞍点策略可以通过求解微分Riccati方程得到\dot S S A A S - S B_p R_p^{-1} B_p S γ^{-2} S B_e R_e^{-1} B_e S Q 0末端条件 S(T) Q_T最优控制u_p -R_p^{-1} B_p S xu_e R_e^{-1} B_e S x / γ^2从形式上看追击方和逃逸方的控制增益都由同一个S矩阵决定只是角色不同。也就是说博弈双方的控制策略本质上是在抢同一个“Riccati状态”的不同分量。这个形式在代码里非常好实现先用数值积分从T到0逆时间积分得到S(t)再在正向仿真中根据当前时刻的S(t)和状态x计算出双方的控制指令。需要提醒的是Riccati方程求解时要特别注意数值稳定性尤其是当S矩阵在末端附近变化剧烈时。我在Matlab里习惯用ode15s求解这类刚性方程而不建议用ode45硬积分——后者在末端几秒内容易出现矩阵发散导致控制增益出现巨大尖峰。3.3 自适应博弈闭环估计与决策的耦合方式整个闭环的运行逻辑可以这样理解EKF不断输出当前时刻对θ的最优估计把θ估计值代入R_e(θ)重新生成Riccati方程中的参数矩阵重新求解S(t)更新追击方的控制增益逃逸方也按照当前估计采用同类策略双方按照更新后的策略执行控制产生新的相对运动量测再反馈给EKF这个“估计-重构-决策-执行”循环不是同时发生的而是一个带时序的串行过程。工程实现上通常把EKF更新频率设为0.1秒甚至0.05秒而Riccati方程的重解频率可以设低一些比如1秒一次。因为Riccati方程重解的计算量相对较大没必要每个滤波周期都重算而估计参数的变化相对平缓1秒重解一次已经足够跟上末端信息更新的节奏。我在复现中发现如果每一个EKF周期都重解Riccati方程仿真速度会明显下降而且控制曲线抖动增多反而导致末端脱靶量略微变差。这是因为参数估计本身有噪声频繁重解会把估计噪声直接传导到控制端。把重解频率降下来相当于给决策回路增加了一个低通特性实际胜率反而提升了。4. Matlab代码实现从模块划分到关键函数解析4.1 代码整体架构与文件组织复现代码最好模块化否则参数多了根本没法维护。以下是我推荐的目录结构分别对应建模、滤波、均衡求解和仿真主循环main.m主脚本负责设置仿真参数、调用模块、记录日志dynamics_CW.mC-W方程相对运动函数输入状态、控制输出状态导数jacobian_ekf.mEKF状态方程雅可比矩阵ekf_predict.m/ekf_update.mEKF预测与更新game_riccati.m求解博弈Riccati微分方程game_control.m根据S(t)和当前状态计算双双方控制指令plot_results.m绘图与结果统计这个划分最核心的好处是你可以独立替换dynamics_CW.m为非线性轨道模型或加入J2摄动只需要保证状态维度和雅可比接口一致你也可以独立替换game_riccati.m里的求解器而不用碰EKF代码。对复现工作而言这种解耦能够极大地方便后续做敏感性分析。4.2 EKF核心代码的关键片段EKF部分的Matlab代码并不长但每一行都可能影响结果。放一段核心更新代码结合注释说明关键设置% EKF update for augmented state function [X_aug_new, P_new] ekf_update(X_aug_prior, P_prior, Z_meas, Hk, Rk) % innovation Z_pred Hk * X_aug_prior; innov Z_meas - Z_pred; S Hk * P_prior * Hk Rk; K P_prior * Hk * inv(S); X_aug_new X_aug_prior K * innov; P_new (eye(length(X_aug_prior)) - K * Hk) * P_prior; end这里有两个地方值得专门说。第一inv(S)虽然读起来方便但在尺寸较大时建议用S \ (Hk * P_prior)代替数值稳定性更好也避免求逆的显式矩阵。第二协方差矩阵P_new的更新我用了最基础的Joseph形式下的简化版本在增益和Hk足够精确时没问题但如果你发现协方差变得不对称或非正定建议补上Joseph形式的更新并定期做对称化处理P_new (P_new P_new) / 2。在调用这个更新之前我们需要从预测模块拿到X_aug_prior和P_prior而预测模块最关键的是雅可比矩阵和连续状态方程向离散EKF的转换。预测代码大致如下function [X_aug_pred, P_pred] ekf_predict(X_aug_est, P_est, u_p, Ts, Q_d) x X_aug_est(1:6); theta_est X_aug_est(7:end); % 当前估计的逃逸方控制 [u_e_est, ~] game_control_escape(x, theta_est); f dynamics_CW_aug(X_aug_est, u_p, u_e_est); X_aug_pred X_aug_est f * Ts; F jacobian_ekf(X_aug_est, u_p, theta_est); Phi expm(F * Ts); P_pred Phi * P_est * Phi Q_d; end注意这里有一个隐含的处理逃逸方控制u_e_est是基于当前状态和θ估计生成的最优响应预测。这相当于假设EKF内部有一个“内模”在模拟对方决策。很多复现代码没有显式做这一步直接把逃逸方控制当扰动但那样EKF的预测方程与实际动力学偏差增大估计精度会明显下降。把这个内模环节做进EKF预测是让滤波跟随对方机动变化的关键。4.3 博弈Riccati方程求解的代码实现求解Riccati微分方程我建议用函数句柄加ode15s这样代码清晰又稳定。核心代码片段如下function [t_grid, S_grid] solve_game_riccati(A, Bp, Be, Q_T, Rp, Re_theta, gamma, t_end, n_time) % 反向积分Riccati方程 S0 Q_T; Svec0 S0(:); % 6x6矩阵展成36维向量 opt odeset(RelTol,1e-6,AbsTol,1e-8); [t_out, Svec_out] ode15s((t,Svec) riccati_rhs(t,Svec,A,Bp,Be,Q_T,Rp,Re_theta,gamma), ... [t_end 0], Svec0, opt); % 反向积分拿到的t是递减的后续插值需要转成递增 t_grid flipud(t_out); S_grid flipud(reshape(Svec_out, size(A))); endriccati_rhs里的核心就是前面提到的微分方程这里不再重复。需要注意Re_theta是随θ估计值变化的所以求解器每次调用都要传入当前θ估计值对应的矩阵。在仿真循环里每1秒重解一次Riccati然后把S_grid插值到控制律计算所需的每个时刻。还有一个小技巧Riccati方程在末端附近S矩阵可能存在快速变化。如果在某些时间点上S(t)变成非对称或者出现明显尖峰不要只调求解器先检查一下Re(theta)是否退化——比如θ估计出接近0或非常大的值导致矩阵病态。在我做参数敏感性分析时θ估计值偶尔会收敛到负值此时R_e(θ)变成负定矩阵Riccati方程完全失去博弈意义末端脱靶量直接爆表。我的解决方式是在EKF更新后对θ做一次投影theta_est max(theta_est, theta_min)限制在物理允许范围内。4.4 仿真主循环设计与结果后处理仿真主循环的结构大致如下初始化状态x0、协方差P0、真实逃逸方参数theta_true、估计初值theta0对于每个采样时刻用当前S_grid插值得到控制增益计算追击方控制u_p真实逃逸方控制用theta_true计算我们能看到“上帝视角”但滤波器看不到一步数值积分得到真实相对运动状态生成带噪声量测EKF预测-更新更新theta估计判断是否到Riccati重解时刻如果是重新求解S(t)记录每个时刻的状态、控制、估计值、协方差迹绘制相对轨迹、控制曲线、θ估计收敛曲线、脱靶量变化等在结果后处理中最值得关注的三张图分别是相对轨迹图看末端是否有明显螺旋收敛或发散θ估计收敛曲线看参数估计是否稳定收敛到真值附近脱靶量随末端时间变化曲线看策略是否真的达到性能要求这三张图是论文里最常出现的三类图复现时要确保它们能对应上。如果仿真跑完发现轨迹很漂亮但θ估计不收敛那说明EKF的测量配置或噪声设置有问题先去调滤波参数不要急着调博弈部分。反过来如果θ估计收敛但脱靶量表现不佳那多半是Riccati求解参数或末端权重设置不合理。5. 复现过程中遇到的典型问题与调试实录5.1 EKF发散从数值稳定性到模型失配我在调EKF时遇到过最典型的现象是滤波跑几秒后协方差矩阵的迹突然急剧增大估计值飞到物理不可信的区间。排查下来常见的原因有三个。第一个是雅可比矩阵错误。尤其当状态方程里同时包含追击控制输入项和逃逸方控制内模项时链式法则非常容易漏项。排查办法很简单用数值差分法对比自动求导得到的雅可比矩阵把二者差值打印出来通常几分钟就能定位到问题。第二个是协方差矩阵非正定。上面提到过当量测更新频繁且增益较大时简单形式的协方差更新可能导致对称性丢失。解决方式是改用Joseph形式并在每个EKF周期末尾强制对称化。第三个是过程噪声矩阵Q设置太小导致滤波器过于相信模型预测。如果模型本身与实际动力学有偏差Q太小时估计会僵化跟不上测量变化。我一般先用一段没有博弈的匀速相对运动数据做EKF预热把Q和R调到比较合理的初值再接上完整的博弈策略仿真。5.2 Riccati方程数值发散刚性问题与病态参数Riccati方程本质上是一个非线性矩阵微分方程末端边界条件给定反向积分时容易出现刚性。如果积分器报错或者输出NaN优先检查以下三项矩阵A是否稳定C-W方程的A矩阵虽然是线性的但在某些参考轨道参数下包含正实部特征值Riccati方程反向积分时可能不稳定Re(θ)是否正定θ估计值掉到负区间会导致完全病态γ是否设置合理γ太小意味着逃逸方能量上限被卡得太死Riccati方程容易出现数值奇异性我在调试时习惯把Riccati方程分解成两个阶段来处理前80%时间段用较粗的相对误差容限后20%时间段用更严格的容限。这样既保证了末端附近的高精度又不会让整段时间的积分变得太慢。5.3 θ估计收敛慢或不收敛信息激励不足末段追逃的难点在于逃逸方的机动可能要到末端十几秒内才真正激发出来前段飞行过程中逃逸方几乎不怎么做大机动此时量测对θ的敏感性很低θ估计基本停留在初值附近。这是不完全信息博弈的天然难点不是单纯调滤波参数能解决的。我做了几组实验后得到一个有效的处理思路在EKF里把θ的过程噪声Q_theta设成相对较大的值这样允许θ在后段快速变化。尽管这意味着稳态时θ的波动会大一些但换来的是末端机动段里θ能够快速跟上。如果Q_theta太小θ会一直“赖”在初值附近末端即使量测信息丰富也来不及收敛。另一个实用技巧是给量测增加末端时段更高权重的处理。因为末端数据信息含量最高可以让滤波在末端时段更激进而中段更平缓。这种自适应噪声调度实现起来很简单把R矩阵按剩余时间做一个递减函数。实测下来对θ收敛速度的提升比较明显而且实现成本很低强烈推荐尝试。5.4 末端脱靶量表现不理想先看控制曲线再下结论如果最终脱靶量比你预期的大不要急着改权重矩阵。先看追击方和逃逸方的控制加速度曲线检查是否存在控制饱和或者振荡。博弈策略的最优控制公式只是无约束情况下的数学解一旦触及推力上限实际控制指令会被截断等于破坏了最优性条件此时末端脱靶量的恶化几乎是必然的。我的做法是在仿真里显式加入控制限幅然后在限幅约束下重新校准R_p和γ的取值。具体来说通过蒙特卡洛扫描R_p的权重观察控制量是否贴近限幅边界选一个既不过饱和又能保证脱靶量满意的折中值。这里没有万能参数只能针对具体场景做标定。还有一个经常被忽略的点末端时间T的选取直接影响博弈效果。太短的T意味着双方几乎没有时间响应太长的T则可能让线性化C-W方程的误差积累过大。我做过一组T从10秒到60秒的扫描发现对于10公里级初始距离T在20到30秒时博弈策略的优势最明显。这给复现时设置仿真场景提供了参考依据。6. 一些个人的实操心得复现这类不完全信息博弈策略整体难度并不在于某个单点算法而在于把EKF、博弈求解、动力学仿真这三个模块咬合成一个能稳定收敛的闭环。EKF输出一个不完美但可用的θ估计Riccati求解器基于这个估计生成控制指令控制指令又反过来影响量测和估计任何一个环节的误差都会在闭环中被放大。调试时建议一步步来先做“完全信息”基线仿真把Riccati博弈控制单独验证到脱靶量满意再冻结博弈策略只调EKF收敛最后才合并起来做联合仿真。这样一旦联合仿真出问题你能快速定位是哪个环节破坏的。另外一个值得留意的点是控制器重解频率和下行的数值插值。如果重解频率过低S(t)插值出来的控制增益和真实最优解偏差大过高又会让估计噪声传导到控制端。我最终选择的是1秒重解一次并且在每个周期内用三次样条插值效果比线性插值更平滑控制曲线几乎没毛刺。这套自适应博弈的思路除了追逃拦截之外在航天器绕飞规避、空间碎片主动清除的接近操作、以及多航天器编队对抗中都可能用到。换个场景时只需要替换动力学方程和量测方程博弈求解框架和EKF架构基本不用大动。从工程复现的角度看这种“感知-估计-均衡-执行”的闭环结构才是整个工作的核心资产。 
长河君 · 2026/10/20 11:37:59 · 阅读 1.2k 不完全信息末端追逃博弈的EKF参数估计与Epsilon纳什均衡Matlab复现
追逃博弈在航天领域并不新鲜但把“不完全信息”和“Epsilon纳什均衡”组合进末端追逃问题再叠加上EKF在线参数估计就值得认真复现一次了。这篇文章我想完整拆解一下这一类基于 Matlab 的复现工作到底在做什么从建模、EKF的“引线”作用、均衡求解到代码实现和调试时踩过的坑一次性讲透希望能给正在搞航天制导控制、空间博弈算法或者想复现类似期刊工作的朋友提供一些参考。很多刚接触这个方向的人会被“Epsilon纳什均衡”这几个字吓到觉得理论门槛很高。实际上我们在末端追逃场景里提到的Epsilon纳什均衡落地到算法层面就是一句话受限于在线估计精度和模型简化博弈双方做不到绝对的、精准的最优响应只能做到“误差有界”的最优响应这个误差界的宽窄就对应Epsilon。理解了这层物理含义后续推导和Matlab代码就会顺很多。本文涉及的内容包括相对运动建模、扩展卡尔曼滤波EKF对对手参数的在线估计、基于估计值求解近似纳什均衡策略以及整套流程的Matlab仿真实测。1. 末端追逃博弈的第一个问题动力学模型的尺度选择1.1 为什么博弈框架在这个尺度上成立末端追逃的典型场景是两个航天器在最后几十公里甚至几公里的相对距离上进行接近与规避对抗。追击方想最小化末端脱靶量逃逸方想最大化末端脱靶量双方同时决策、同时机动这天然是一个双方零和微分博弈问题。和传统单边最优制导最大的区别在于博弈框架下追击方的策略必须考虑逃逸方的“理性反击”而不是简单地把对手当作固定目标或已知轨迹。有人会问直接用鲁棒控制或者滑模制导不也能处理外界扰动吗从工程结果上看某些鲁棒制导的确能扛住一定程度的机动不确定性但鲁棒控制不关心对手是否也在“优化自己的收益”它只是把所有不确定性当作最坏情况的干扰。而博弈策略的核心优势是能够对齐双方的“决策交互”对手机动策略变化时追击方的策略可以自适应地跟着调整。在EKF估计的配合下这种调整不仅是基于状态的反馈还是基于对手“底牌”的反馈。1.2 相对运动方程选择C-W方程还是非线性轨道方程末端追逃问题里最常见的做法是选用C-W方程Hill方程描述相对运动因为它形式简单、解析性质好方便构造线性二次型博弈。线性化的C-W方程如下相对位置向量 r [x; y; z]相对速度向量 v [vx; vy; vz]在圆轨道参考系下相对运动近似满足 \dot r v\dot v_x 3 n^2 x 2 n vy u_x d_x\dot v_y -2 n vx u_y d_y\dot v_z -n^2 z u_z d_z其中 n 表示参考轨道角速度u 为控制加速度d 为机动加速度扰动如果逃逸方的机动能力用有界控制项表示整个方程可以整理成状态空间形式\dot X A X B_p u_p B_e u_e。A和B矩阵的具体形式和参考轨道参数有关复现时只要把轨道高度、偏心率这些参数代入即可。但要注意C-W方程只适用于近距离、短时间、近圆轨道相对运动。如果初始相对距离超过几十公里或者双方机动加速度很大线性误差会明显影响结果的合理性。我在复现时对比过线性模型和非线性轨道相对运动模型的差异发现当末端时间在30秒以内、相对距离小于10公里时C-W方程和数值轨道模型的相对位置误差能控制在5%以内这个精度对博弈策略验证是够用的但如果末端时间拉到上百秒线性简化带来的偏差会导致末端脱靶量出现“假性收敛”也就是仿真看着效果不错实际场景中却跑偏。所以复现这类期刊工作时建议先明确末端时间和初始距离的量级再决定模型复杂度。对于末端追逃这种“末端”属性很强的问题短时间、近距离场景下用C-W方程完全合理。1.3 不完全信息的数学描述把“不知道”翻译成参数化形式不完全信息在博弈论里有很多种表达方式在追逃问题里最常见的是对手的成本函数参数、机动能力上限、或者终端偏好权重未知。举个具体的例子追击方知道逃逸方会最小化自己被捕获的概率但不知道逃逸方对控制能量消耗的“心疼程度”——是愿意不惜一切代价机动还是倾向于省燃料只做小幅规避。这个偏好权重就是逃逸方的“底牌”。在数学上我们把逃逸方的控制加权矩阵 R_e 写成已知结构 未知参数的形式比如 R_e(θ) diag([θ1, θ2, θ3])θ 就是要估计的未知参数向量。这样不完全信息问题就转换为一个参数辨识问题。此时“不完全信息”不再是抽象意义上的“敌方意图”而是具体到某些能被EKF实时估计的参数。在这个转换中一个很重要的前提是要保证系统可辨识性。也就是说需要确认测量数据能够把 θ 和状态量区分开。复现时我专门看了量测矩阵对扩展状态的敏感性如果只有视线角测量θ 的估计收敛就比较慢如果加上相对距离和相对速度测量θ 的收敛速度会明显加快。这也是很多期刊复现代码里把测量模型设置为“相对位置相对速度”的原因。2. EKF参数估计把对手的“底牌”从数据里滤出来2.1 为什么纳什均衡需要参数估计如果假设完全信息博弈双方的参数都已知那么整个问题就是一个标准的二次型微分博弈可以直接求解Riccati微分方程得到反馈控制策略。但现实里追击方并不掌握逃逸方的真实偏好参数那么计算出的“最优”策略很可能不是真正针对对手的最优响应甚至可能被对手的机动方式牵着走。EKF在这里起的作用就是“在线采参数”通过量测数据估计逃逸方的未知参数随着滤波收敛把原本不确定的模型逐步逼近真实模型然后再基于这个逼近后的模型求解博弈控制。这样一来整个闭环系统就有了“感知-估计-决策-执行”的完整链条也就是标题里“自适应博弈”的含义。有朋友可能会问直接用强跟踪滤波STF或者无迹卡尔曼UKF是否更好对这个问题我自己的体会是EKF在计算效率和实现简洁性上最平衡尤其当你已经把状态方程线性化模板写好之后切换到不同场景只需要改雅可比矩阵。而UKF虽然对强非线性更友好但Sigma点传播会引入额外计算量在机器人或航天嵌入式平台上不太划算。期刊复现里大多以EKF为主也是因为它的计算量小方便快速验证博弈算法的核心逻辑。2.2 扩展状态构造与EKF预测更新流程构造扩展状态向量时我们把原来的相对运动状态和未知参数拼接在一起X_aug [x; y; z; vx; vy; vz; θ1; θ2; θ3]对应的状态方程为前6个状态遵循C-W方程后3个参数状态取 \dot θ 0表示在滤波时间尺度内对手偏好不突变量测方程取Z H X_aug v_noise其中H为线性量测矩阵v_noise为零均值高斯白噪声。因为EKF的状态方程本身带有控制输入项而这个控制输入项在估计时是已知的追击方自己发出的控制指令所以EKF预测阶段要把 u_p 代进去。逃逸方的控制则依赖当前估计参数在滤波里用一个“当前最优猜测”近似形成“EKF估计-博弈求解-实际执行”的循环。EKF每个滤波周期的过程大致是计算雅可比矩阵 F ∂f/∂X_augΦ expm(F * T_sample)预测协方差 P_pred Φ * P * Φ Q计算卡尔曼增益 K P_pred * H * (H * P_pred * H R)^{-1}状态更新 X_new X_pred K * (Z_real - Z_pred)协方差更新 P_new (I - K * H) * P_pred这里最容易出事的是雅可比矩阵的推导。如果状态方程里同时包含追击方控制 u_p 和逃逸方估计控制 \hat u_e那么F要对状态求偏导而 \hat u_e 本身又是估计参数 θ 的函数这一步需要把链式法则写对。我自己在复现时就在这里犯过错误F矩阵少了一项导致EKF在末端十几秒时估计值突然发散。后来把F矩阵的每个元素都做了符号求导验证才定位到问题。2.3 可观性分析与测量配置的工程考量EKF能不能收敛很大程度上取决于扩展状态的可观性。我们做了简化推导如果量测只是相对位置的三个分量那么系统对 θ 的可观性矩阵往往不满秩尤其当相对加速度方向不变时θ的某些分量无法被有效激励。实际问题里逃逸方机动通常是在末端才开始明显变化所以前段飞行时 θ 估计波动大末端机动激发后才逐渐收敛。工程上的对策是增加量测维度把相对速度甚至视线角速率引入量测。这样虽然EKF的矩阵维度变大但可观性明显改善。我在代码里把量测设置为相对速度和相对位置同时测量对比实验表明 θ 估计的均方根误差从20%以上降到了8%左右。另外有一个细节EKF的噪声矩阵R不能设得太小否则滤波增益过高量测噪声会被放大到状态估计中导致博弈策略抖动也不能设得太大否则EKF跟踪不上末端快速变化的相对运动。一个可行的调参方案是从名义值出发通过蒙特卡洛仿真对比估计误差和末端脱靶量的趋势选定一个在中段和末端都能接受的折中值。3. Epsilon纳什均衡在不完美估计下定义“次优但够用”的策略3.1 不完全信息均衡的定义回顾纳什均衡的核心思想是没有任何一个参与者可以单方面改变策略获得更好收益。完全信息下这个定义是严格的但在不完全信息下由于存在估计误差策略不是建立在“真实参数”上而是建立在“估计参数”上因此很难达到精确的纳什均衡。Epsilon纳什均衡的定义很简洁对每个参与者如果其策略相对于任意单边偏离所能带来的收益增益不超过一个很小的正数ε则称该策略组合构成ε-纳什均衡。放到追逃里就是说追击方和逃逸方各自按照估计模型采用近似最优策略双方都不存在能使自己多获得超过ε的收益的“单方面聪明改变”。这个定义既保留了纳什均衡的核心精神又给在线估计和模型简化留了余地。复现中判断ε是否可接受通常直接看性能指标差异拿估计策略和真实完全信息策略分别计算末端脱靶量如果差异足够小就认为ε在工程可接受范围内。这个方法虽然有点“事后验证”的味道但对仿真验证来说是最直观有效的。3.2 末端追逃博弈策略的求解路径在短时间末端场景下线性二次型微分博弈是最常用的可解析求解框架。目标函数一般写成J 1/2 * x(T) * Q_T * x(T) 1/2 * ∫_0^T [u_p * R_p * u_p - γ^2 * u_e * R_e(θ) * u_e] dt这里追击方最小化J逃逸方最大化J。Q_T是末端相对位置误差的惩罚矩阵R_p和R_e是控制能量惩罚。γ是逃逸方的“能量上限系数”反映了逃逸方的机动能力边界。对于这个目标函数标准结果是鞍点策略可以通过求解微分Riccati方程得到\dot S S A A S - S B_p R_p^{-1} B_p S γ^{-2} S B_e R_e^{-1} B_e S Q 0末端条件 S(T) Q_T最优控制u_p -R_p^{-1} B_p S xu_e R_e^{-1} B_e S x / γ^2从形式上看追击方和逃逸方的控制增益都由同一个S矩阵决定只是角色不同。也就是说博弈双方的控制策略本质上是在抢同一个“Riccati状态”的不同分量。这个形式在代码里非常好实现先用数值积分从T到0逆时间积分得到S(t)再在正向仿真中根据当前时刻的S(t)和状态x计算出双方的控制指令。需要提醒的是Riccati方程求解时要特别注意数值稳定性尤其是当S矩阵在末端附近变化剧烈时。我在Matlab里习惯用ode15s求解这类刚性方程而不建议用ode45硬积分——后者在末端几秒内容易出现矩阵发散导致控制增益出现巨大尖峰。3.3 自适应博弈闭环估计与决策的耦合方式整个闭环的运行逻辑可以这样理解EKF不断输出当前时刻对θ的最优估计把θ估计值代入R_e(θ)重新生成Riccati方程中的参数矩阵重新求解S(t)更新追击方的控制增益逃逸方也按照当前估计采用同类策略双方按照更新后的策略执行控制产生新的相对运动量测再反馈给EKF这个“估计-重构-决策-执行”循环不是同时发生的而是一个带时序的串行过程。工程实现上通常把EKF更新频率设为0.1秒甚至0.05秒而Riccati方程的重解频率可以设低一些比如1秒一次。因为Riccati方程重解的计算量相对较大没必要每个滤波周期都重算而估计参数的变化相对平缓1秒重解一次已经足够跟上末端信息更新的节奏。我在复现中发现如果每一个EKF周期都重解Riccati方程仿真速度会明显下降而且控制曲线抖动增多反而导致末端脱靶量略微变差。这是因为参数估计本身有噪声频繁重解会把估计噪声直接传导到控制端。把重解频率降下来相当于给决策回路增加了一个低通特性实际胜率反而提升了。4. Matlab代码实现从模块划分到关键函数解析4.1 代码整体架构与文件组织复现代码最好模块化否则参数多了根本没法维护。以下是我推荐的目录结构分别对应建模、滤波、均衡求解和仿真主循环main.m主脚本负责设置仿真参数、调用模块、记录日志dynamics_CW.mC-W方程相对运动函数输入状态、控制输出状态导数jacobian_ekf.mEKF状态方程雅可比矩阵ekf_predict.m/ekf_update.mEKF预测与更新game_riccati.m求解博弈Riccati微分方程game_control.m根据S(t)和当前状态计算双双方控制指令plot_results.m绘图与结果统计这个划分最核心的好处是你可以独立替换dynamics_CW.m为非线性轨道模型或加入J2摄动只需要保证状态维度和雅可比接口一致你也可以独立替换game_riccati.m里的求解器而不用碰EKF代码。对复现工作而言这种解耦能够极大地方便后续做敏感性分析。4.2 EKF核心代码的关键片段EKF部分的Matlab代码并不长但每一行都可能影响结果。放一段核心更新代码结合注释说明关键设置% EKF update for augmented state function [X_aug_new, P_new] ekf_update(X_aug_prior, P_prior, Z_meas, Hk, Rk) % innovation Z_pred Hk * X_aug_prior; innov Z_meas - Z_pred; S Hk * P_prior * Hk Rk; K P_prior * Hk * inv(S); X_aug_new X_aug_prior K * innov; P_new (eye(length(X_aug_prior)) - K * Hk) * P_prior; end这里有两个地方值得专门说。第一inv(S)虽然读起来方便但在尺寸较大时建议用S \ (Hk * P_prior)代替数值稳定性更好也避免求逆的显式矩阵。第二协方差矩阵P_new的更新我用了最基础的Joseph形式下的简化版本在增益和Hk足够精确时没问题但如果你发现协方差变得不对称或非正定建议补上Joseph形式的更新并定期做对称化处理P_new (P_new P_new) / 2。在调用这个更新之前我们需要从预测模块拿到X_aug_prior和P_prior而预测模块最关键的是雅可比矩阵和连续状态方程向离散EKF的转换。预测代码大致如下function [X_aug_pred, P_pred] ekf_predict(X_aug_est, P_est, u_p, Ts, Q_d) x X_aug_est(1:6); theta_est X_aug_est(7:end); % 当前估计的逃逸方控制 [u_e_est, ~] game_control_escape(x, theta_est); f dynamics_CW_aug(X_aug_est, u_p, u_e_est); X_aug_pred X_aug_est f * Ts; F jacobian_ekf(X_aug_est, u_p, theta_est); Phi expm(F * Ts); P_pred Phi * P_est * Phi Q_d; end注意这里有一个隐含的处理逃逸方控制u_e_est是基于当前状态和θ估计生成的最优响应预测。这相当于假设EKF内部有一个“内模”在模拟对方决策。很多复现代码没有显式做这一步直接把逃逸方控制当扰动但那样EKF的预测方程与实际动力学偏差增大估计精度会明显下降。把这个内模环节做进EKF预测是让滤波跟随对方机动变化的关键。4.3 博弈Riccati方程求解的代码实现求解Riccati微分方程我建议用函数句柄加ode15s这样代码清晰又稳定。核心代码片段如下function [t_grid, S_grid] solve_game_riccati(A, Bp, Be, Q_T, Rp, Re_theta, gamma, t_end, n_time) % 反向积分Riccati方程 S0 Q_T; Svec0 S0(:); % 6x6矩阵展成36维向量 opt odeset(RelTol,1e-6,AbsTol,1e-8); [t_out, Svec_out] ode15s((t,Svec) riccati_rhs(t,Svec,A,Bp,Be,Q_T,Rp,Re_theta,gamma), ... [t_end 0], Svec0, opt); % 反向积分拿到的t是递减的后续插值需要转成递增 t_grid flipud(t_out); S_grid flipud(reshape(Svec_out, size(A))); endriccati_rhs里的核心就是前面提到的微分方程这里不再重复。需要注意Re_theta是随θ估计值变化的所以求解器每次调用都要传入当前θ估计值对应的矩阵。在仿真循环里每1秒重解一次Riccati然后把S_grid插值到控制律计算所需的每个时刻。还有一个小技巧Riccati方程在末端附近S矩阵可能存在快速变化。如果在某些时间点上S(t)变成非对称或者出现明显尖峰不要只调求解器先检查一下Re(theta)是否退化——比如θ估计出接近0或非常大的值导致矩阵病态。在我做参数敏感性分析时θ估计值偶尔会收敛到负值此时R_e(θ)变成负定矩阵Riccati方程完全失去博弈意义末端脱靶量直接爆表。我的解决方式是在EKF更新后对θ做一次投影theta_est max(theta_est, theta_min)限制在物理允许范围内。4.4 仿真主循环设计与结果后处理仿真主循环的结构大致如下初始化状态x0、协方差P0、真实逃逸方参数theta_true、估计初值theta0对于每个采样时刻用当前S_grid插值得到控制增益计算追击方控制u_p真实逃逸方控制用theta_true计算我们能看到“上帝视角”但滤波器看不到一步数值积分得到真实相对运动状态生成带噪声量测EKF预测-更新更新theta估计判断是否到Riccati重解时刻如果是重新求解S(t)记录每个时刻的状态、控制、估计值、协方差迹绘制相对轨迹、控制曲线、θ估计收敛曲线、脱靶量变化等在结果后处理中最值得关注的三张图分别是相对轨迹图看末端是否有明显螺旋收敛或发散θ估计收敛曲线看参数估计是否稳定收敛到真值附近脱靶量随末端时间变化曲线看策略是否真的达到性能要求这三张图是论文里最常出现的三类图复现时要确保它们能对应上。如果仿真跑完发现轨迹很漂亮但θ估计不收敛那说明EKF的测量配置或噪声设置有问题先去调滤波参数不要急着调博弈部分。反过来如果θ估计收敛但脱靶量表现不佳那多半是Riccati求解参数或末端权重设置不合理。5. 复现过程中遇到的典型问题与调试实录5.1 EKF发散从数值稳定性到模型失配我在调EKF时遇到过最典型的现象是滤波跑几秒后协方差矩阵的迹突然急剧增大估计值飞到物理不可信的区间。排查下来常见的原因有三个。第一个是雅可比矩阵错误。尤其当状态方程里同时包含追击控制输入项和逃逸方控制内模项时链式法则非常容易漏项。排查办法很简单用数值差分法对比自动求导得到的雅可比矩阵把二者差值打印出来通常几分钟就能定位到问题。第二个是协方差矩阵非正定。上面提到过当量测更新频繁且增益较大时简单形式的协方差更新可能导致对称性丢失。解决方式是改用Joseph形式并在每个EKF周期末尾强制对称化。第三个是过程噪声矩阵Q设置太小导致滤波器过于相信模型预测。如果模型本身与实际动力学有偏差Q太小时估计会僵化跟不上测量变化。我一般先用一段没有博弈的匀速相对运动数据做EKF预热把Q和R调到比较合理的初值再接上完整的博弈策略仿真。5.2 Riccati方程数值发散刚性问题与病态参数Riccati方程本质上是一个非线性矩阵微分方程末端边界条件给定反向积分时容易出现刚性。如果积分器报错或者输出NaN优先检查以下三项矩阵A是否稳定C-W方程的A矩阵虽然是线性的但在某些参考轨道参数下包含正实部特征值Riccati方程反向积分时可能不稳定Re(θ)是否正定θ估计值掉到负区间会导致完全病态γ是否设置合理γ太小意味着逃逸方能量上限被卡得太死Riccati方程容易出现数值奇异性我在调试时习惯把Riccati方程分解成两个阶段来处理前80%时间段用较粗的相对误差容限后20%时间段用更严格的容限。这样既保证了末端附近的高精度又不会让整段时间的积分变得太慢。5.3 θ估计收敛慢或不收敛信息激励不足末段追逃的难点在于逃逸方的机动可能要到末端十几秒内才真正激发出来前段飞行过程中逃逸方几乎不怎么做大机动此时量测对θ的敏感性很低θ估计基本停留在初值附近。这是不完全信息博弈的天然难点不是单纯调滤波参数能解决的。我做了几组实验后得到一个有效的处理思路在EKF里把θ的过程噪声Q_theta设成相对较大的值这样允许θ在后段快速变化。尽管这意味着稳态时θ的波动会大一些但换来的是末端机动段里θ能够快速跟上。如果Q_theta太小θ会一直“赖”在初值附近末端即使量测信息丰富也来不及收敛。另一个实用技巧是给量测增加末端时段更高权重的处理。因为末端数据信息含量最高可以让滤波在末端时段更激进而中段更平缓。这种自适应噪声调度实现起来很简单把R矩阵按剩余时间做一个递减函数。实测下来对θ收敛速度的提升比较明显而且实现成本很低强烈推荐尝试。5.4 末端脱靶量表现不理想先看控制曲线再下结论如果最终脱靶量比你预期的大不要急着改权重矩阵。先看追击方和逃逸方的控制加速度曲线检查是否存在控制饱和或者振荡。博弈策略的最优控制公式只是无约束情况下的数学解一旦触及推力上限实际控制指令会被截断等于破坏了最优性条件此时末端脱靶量的恶化几乎是必然的。我的做法是在仿真里显式加入控制限幅然后在限幅约束下重新校准R_p和γ的取值。具体来说通过蒙特卡洛扫描R_p的权重观察控制量是否贴近限幅边界选一个既不过饱和又能保证脱靶量满意的折中值。这里没有万能参数只能针对具体场景做标定。还有一个经常被忽略的点末端时间T的选取直接影响博弈效果。太短的T意味着双方几乎没有时间响应太长的T则可能让线性化C-W方程的误差积累过大。我做过一组T从10秒到60秒的扫描发现对于10公里级初始距离T在20到30秒时博弈策略的优势最明显。这给复现时设置仿真场景提供了参考依据。6. 一些个人的实操心得复现这类不完全信息博弈策略整体难度并不在于某个单点算法而在于把EKF、博弈求解、动力学仿真这三个模块咬合成一个能稳定收敛的闭环。EKF输出一个不完美但可用的θ估计Riccati求解器基于这个估计生成控制指令控制指令又反过来影响量测和估计任何一个环节的误差都会在闭环中被放大。调试时建议一步步来先做“完全信息”基线仿真把Riccati博弈控制单独验证到脱靶量满意再冻结博弈策略只调EKF收敛最后才合并起来做联合仿真。这样一旦联合仿真出问题你能快速定位是哪个环节破坏的。另外一个值得留意的点是控制器重解频率和下行的数值插值。如果重解频率过低S(t)插值出来的控制增益和真实最优解偏差大过高又会让估计噪声传导到控制端。我最终选择的是1秒重解一次并且在每个周期内用三次样条插值效果比线性插值更平滑控制曲线几乎没毛刺。这套自适应博弈的思路除了追逃拦截之外在航天器绕飞规避、空间碎片主动清除的接近操作、以及多航天器编队对抗中都可能用到。换个场景时只需要替换动力学方程和量测方程博弈求解框架和EKF架构基本不用大动。从工程复现的角度看这种“感知-估计-均衡-执行”的闭环结构才是整个工作的核心资产。