1. 项目到底在解决什么问题不完全信息下的末端追逃博弈先直接说结论这个项目复现的是一类带有不完全信息假设的航天器末端追逃博弈策略参与双方都装备了扩展卡尔曼滤波器EKF来在线估计对方的机动参数并在Epsilon纳什均衡ε-纳什均衡的意义下求解自己的控制策略最终用一套Matlab仿真代码把整个“估计—博弈求解—控制闭环”跑通。追逃博弈这个概念听起来很“军事”但它在学术上属于微分博弈differential game和轨道机动控制交叉的经典问题一个追踪航天器Pursuer试图靠近目标而一个逃逸航天器Evader试图拉开距离。末端段意味着双方已经进入相对距离几十公里甚至几公里的量级此时引力差、J2摄动等长期轨道效应相对次要主导方程是相对运动动力学Clohessy-Wiltshire方程或其非线性版本控制力来自双方自身的推力加速度。以前教材里的追逃博弈总是先假设“完全信息”双方都知道对方的动力学模型、控制约束和收益函数。但工程实际立刻打脸——追踪航天器不可能精确知道逃逸航天器的最大推力、发动机响应延迟、真实导航误差水平逃逸方同样也只知道追踪方的大致性能包线。这就是“不完全信息”。我们复现的期刊策略就是在双方互相“猜参数”的前提下用EKF把未知参数估计出来再基于估计结果做自适应博弈求解。适合看这篇文章的读者有三类第一类是研究航天器相对运动控制的研究生想找一个能跑通的仿真框架第二类是研究博弈论与控制结合的人想理解ε-纳什均衡怎么落到工程实现第三类是刚接触EKF参数辨识的工程师希望看到一个完整的、把状态估计和决策耦合起来的闭环例子。2. 从“完全理性”到“有限理性”为什么必须引入Epsilon纳什均衡2.1 完全信息纳什均衡的工程困境先回顾一下教科书定义。在追逃博弈中设追踪方策略为 (u_p)逃逸方策略为 (u_e)追踪方收益函数为 (J_p(u_p, u_e))越小越好逃逸方收益为 (J_e(u_p, u_e))越大越好。常规纳什均衡 ((u_p^, u_e^)) 满足[ J_p(u_p^, u_e^) \le J_p(u_p, u_e^*), \quad \forall u_p \in \mathcal{U}_p ][ J_e(u_p^, u_e^) \ge J_e(u_p^*, u_e), \quad \forall u_e \in \mathcal{U}_e ]换句话说任何一方单独改变策略都不会让自己更好。这个定义很优美但落到航天器末端追逃场景立即遇到三个现实问题。第一计算代价太高。严格纳什均衡通常涉及求解一对Hamilton-Jacobi-IsaacsHJI偏微分方程或者进行连续策略空间上的不动点迭代。在末端博弈这种时间以秒到分钟计算的场合星载计算机根本扛不住高维PDE数值求解。第二收益估计本身有噪声。因为我们用了EKF去估计对方参数估计误差必然传播到收益函数 (J_p)、(J_e) 的计算上。你辛辛苦苦求出的“严格均衡”可能是基于一个有偏差的收益函数求出来的那就谈不上严格了。第三双方理性程度是有限的。真实航天器控制系统不会真的去追求绝对最优策略它只需要“足够好”——收益损失不超过某个可接受的小量即可。所以项目里选择ε-纳什均衡本质上是把“绝对最优”放松成“近似最优”如果双方收益损失都不超过一个小常数 (\varepsilon)就认为达到了可接受的博弈解。这个放松是完全合理的工程折中。2.2 Epsilon纳什均衡的数学定义给定 (\varepsilon 0)策略组合 ((u_p^, u_e^)) 是ε-纳什均衡当且仅当[ J_p(u_p^, u_e^) \le J_p(u_p, u_e^*) \varepsilon_p, \quad \forall u_p \in \mathcal{U}_p ][ J_e(u_p^, u_e^) \ge J_e(u_p^*, u_e) - \varepsilon_e, \quad \forall u_e \in \mathcal{U}_e ]直观理解追踪方如果单方面换策略最多只能再减少 (\varepsilon_p) 的收益逃逸方最多只能再多拿 (\varepsilon_e) 的收益。只要 (\varepsilon) 取得足够小这个“差不多最优”的策略组合在工程上完全可用。让我用一个生活化类比来解释考试时大家都在有限时间内交卷你不需要保证自己的答案全校最优只需要保证“无论怎么改选某道题总分提升不超过2分”——这个2分就是你的ε。ε越小你和“绝对最优”越接近但解题时间越长ε取得太大策略可能明显偏离最优。实际复现时的难点在于ε到底取多少这个没有普适答案要结合动力学时间常数、EKF收敛速度和控制能力一起来试。2.3 项目中的自适应博弈求解思路在复现期刊代码时我注意到它的求解思路不是一次性求全局均衡而是在每个控制周期内根据当前EKF参数估计局部迭代逼近ε-纳什均衡。具体来说做了三件事用EKF给出对手参数的当前估计值 (\hat{\theta})把这个估计当作“当前时刻的近似真实参数”。在当前参数下构建线性化或局部二次近似的收益函数。用迭代最优响应法Iterative Best Response或梯度投影法迭代若干轮后得到一个满足ε条件的策略组合。这个思路非常实用因为它在每个周期只做有限次迭代计算量可控同时由于参数估计持续更新策略也会自适应调整。仿真里我最关注的就是ε的收敛曲线——如果策略迭代在后期稳定在一个很小的范围说明算法确实收敛到了近似均衡。3. EKF参数辨识与状态估计追逃双方如何“看穿”对方3.1 状态-参数联合估计问题建模不完全信息是分层次的。最底层的不完全信息是“我不知道你的状态”这个靠相对导航传感器测距、测速、视线角解决更深一层是“我不知道你的动力学模型参数”比如逃逸航天器的最大推力加速度 (a_{max})、控制响应时间常数 (\tau)或者追踪方导航系统量测噪声的真实方差 (R)。在Matlab复现里最常用的做法是增广状态向量把未知参数放进滤波状态里一起估计。设原始状态 (x \in \mathbb{R}^6)相对位置和相对速度未知参数为 (\theta \in \mathbb{R}^m)增广状态就是 (X [x; \theta])。例如期望估计追逃双方的“推力上限”则增广状态为[ X [r_x, r_y, r_z, v_x, v_y, v_z, a_{p,max}, a_{e,max}]^T ]为什么可以用一个公共滤波器同时估计状态和参数因为参数 (\theta) 虽然是常数或缓变但它会影响系统的演化方程所以通过观测残差可以反推参数。这就是EKF做在线参数辨识的基本原理。3.2 扩展卡尔曼滤波的核心流程EKF处理的是非线性系统[ X_{k1} f(X_k, u_k, w_k), \quad w_k \sim \mathcal{N}(0, Q) ][ z_k h(X_k) v_k, \quad v_k \sim \mathcal{N}(0, R) ]EKF与传统KF的区别在于它对非线性函数 (f) 和 (h) 在当前估计 (\hat{X}_k) 处做一阶Taylor展开得到雅可比矩阵 (F_k) 和 (H_k)然后套用线性卡尔曼滤波的预测-更新框架。预测步[ \hat{X}{k1|k} f(\hat{X}{k|k}, u_k) ][ P_{k1|k} F_k P_{k|k} F_k^T Q ]更新步[ K_k P_{k1|k} H_k^T \left( H_k P_{k1|k} H_k^T R \right)^{-1} ][ \hat{X}{k1|k1} \hat{X}{k1|k} K_k \left( z_{k1} - h(\hat{X}_{k1|k}) \right) ][ P_{k1|k1} (I - K_k H_k) P_{k1|k} ]这套流程在Matlab里就是几个矩阵操作的事十来行代码就能写出来。真正麻烦的是两个细节。第一个细节是雅可比矩阵的推导。相对运动方程如果是CW方程那对状态部分的雅可比是常值矩阵但如果用了非线性Clohessy-Wiltshire修正项或考虑了J2项就必须手动求偏导或使用数值雅可比。我复现时先用解析雅可比保证精度调试通过后再尝试数值差分雅可比发现只要步长选得合适两者结果几乎一致但解析式更快更稳。第二个细节是参数可观测性。EKF不是万能的。如果观测方程 (h(X)) 完全不包含某个参数的信息那EKF对那个参数的估计值就只由初始猜测决定根本不会修正。典型例子你试图从相对位置观测中估计逃逸航天器的质量但如果观测信息里质量完全不出现在动态方程中比如你用的是线性CW方程那就估计不出来。务必先做可观测性分析或者保底的办法是给参数估计加上小的过程噪声扰动让协方差不至于迅速收缩到零。3.3 噪声方差整定的经验Q和R的选取直接决定EKF收敛速度和稳态精度。我复现时用了两个阶段阶段一前几秒把过程噪声Q设得稍大让滤波器“敢信”测量、快速收敛。比如相对状态部分取 (Q_{state} 10^{-6} \times I) 量级参数部分可以取 (Q_{\theta} 10^{-4})因为参数初值不确定度大。阶段二估计收敛后适当减小Q让滤波器更信任动力学模型降低稳态振荡。参数部分的Q可以降一个数量级但要保留一个小的值比如 (10^{-6})防止EKF过度自信导致参数锁死。R的选取相对直接按传感器噪声方差标定。我假设测距噪声标准差5米、测速噪声标准差0.1米/秒因此R取对角阵 (\text{diag}(25, 25, 25, 0.01, 0.01, 0.01))。如果你用视线角观测R量纲就变成rad²需要按实际传感器重标。4. 自适应博弈策略设计从估计值到控制指令的闭环4.1 闭环架构估计、博弈、控制三步循环整个仿真在每个时间步内走三个环节测量与估计读取相对导航传感器的带噪量测运行EKF得到当前增广状态估计 (\hat{X})。博弈求解以 (\hat{X}) 为初值结合估计出的对手参数利用ε-纳什均衡迭代求解双方控制指令 ((u_p, u_e))。动力学推进将控制指令代入真实仿真动力学推进一个步长生成下一周期的观测数据循环往复。这个架构里博弈求解器是核心。我建议把博弈求解单独封装成一个函数输入是当前状态估计、双方估计参数和当前时间输出是双方控制指令。这样容易调试也方便后续替换不同的求解算法对比效果。4.2 收益函数与控制策略设计项目中最常见的一类收益函数是二次型终端代价加控制代价的混合形式。追踪方想最小化终端相对距离同时不想过度消耗燃料其收益可以写成[ J_p |r(t_f)|^2 \int_{t_0}^{t_f} u_p^T W_p u_p , dt ]逃逸方收益则是类似形式的相反数或另一种组合。在这个框架下博弈问题转化为双方在各自的控制集合 (\mathcal{U}_p)、(\mathcal{U}_e) 内选择控制律使终端指标达到ε-均衡。求解时不能直接解HJI方程而是用迭代最优响应固定逃逸方当前策略 (u_e^{(k)})求解追踪方的最优响应 (u_p^{(k1)} \arg\min J_p(u_p, u_e^{(k)}))固定追踪方策略 (u_p^{(k1)})求解逃逸方最优响应 (u_e^{(k2)} \arg\max J_e(u_p^{(k1)}, u_e))循环迭代直到连续两次迭代的收益差都小于 (\varepsilon_p)、(\varepsilon_e)。每一次“最优响应”求解本质上是一个单方最优控制问题可以用直接法配点法、伪谱法或间接法协态方程打靶求解。但为了满足实时性我复现时采用的是在线线性二次型近似把每个时刻的博弈问题近似成双线性二次博弈最优响应闭式解就能出来迭代几轮后近似得到ε均衡。这个简化在末端较短、状态变化平稳的工况下实测效果足够稳定。4.3 自适应机制如何体现自适应体现在两个层面。第一个层面是对对手参数的在线修正。EKF把当前对 (a_{p,max})、(a_{e,max}) 的估计喂给博弈求解器。比如逃逸方一旦探测到追踪方实际推力更强(a_{p,max}) 估计值上升它就会在博弈求解中更激进地拉开距离、减少接近风险。第二个层面是策略更新频率随估计置信度变化。当EKF协方差较大、估计不可靠时博弈求解器的步长可以适当保守跳变幅度受限当协方差收缩、估计可信时再让策略更快响应。这个机制类似强化学习里的自适应学习率用起来效果很好能显著抑制初期震荡。我个人的一个经验不要把EKF估计值直接当作真值来更新控制最好做一个简单的低通滤波或速率限制否则估计噪声引起的控制抖动会反过来污染EKF的测量输入造成正反馈振荡。这个坑我在初版代码里踩得很深后面在“常见问题”部分细说。5. Matlab仿真实现完整复现流程与关键代码5.1 仿真结构与初始化我复现时的仿真结构分四个层级顶层run_simulation.m设置参数调用循环汇总结果绘图中层sim_step.m单步仿真调度串联EKF、博弈求解、动力学推进底层ekf_update.m和solve_epsilonNE.m分别完成滤波和博弈求解辅助函数dynamics.m、measurement.m生成真实状态演化和带噪观测。初始化参数我建议单列一个脚本init_params.m方便反复调整% init_params.m dt 0.1; % 控制周期/仿真步长 T_end 60; % 总仿真时间 t 0:dt:T_end; % 初始相对状态追踪方在逃逸方后方1km速度差50m/s r0 [1000; 0; 0]; v0 [50; 0; 0]; X0 [r0; v0; 20; 15]; % 增广状态后两位是双方最大推力加速度初值 % 传感器噪声标准差 sigma_r 5; % 位置量测噪声m sigma_v 0.1; % 速度量测噪声m/s R diag([sigma_r^2*ones(1,3), sigma_v^2*ones(1,3)]); % 过程噪声 Q 1e-6 * diag([1 1 1 1 1 1 1e2 1e2]); % 参数部分噪声大一些 % 博弈参数 epsilon_p 0.5; epsilon_e 0.5; max_iter_gd 100; % 博弈求解最大迭代次数5.2 EKF更新函数实现核心的EKF函数大概是这样只列主体逻辑完整代码建议自己补齐function [X_est, P] ekf_update(X_est, P, z, u_p, u_e, dt, Q, R) % 状态转移与雅可比 [X_pred, F] dynamics_full(X_est, u_p, u_e, dt); % 预测协方差 P_pred F * P * F Q; % 观测预测与雅可比 [z_pred, H] measurement_full(X_pred); % 新息 y z - z_pred; S H * P_pred * H R; K P_pred * H / S; % 更新 X_est X_pred K * y; P (eye(length(X_est)) - K * H) * P_pred; % 协方差对称化与正定修正 P 0.5 * (P P); P (P P) / 2; if min(eig(P)) 0 [V, D] eig(P); D max(D, 1e-10); P V * D * V; end end这个函数有几个细节值得注意。第一dynamics_full不仅要返回预测状态还要返回状态雅可比矩阵 (F)。如果你在dynamics.m里只更新状态而不算雅可比那EKF就无从进行协方差预测。我建议dynamics_full的输出格式固定为[X_pred, F]即使暂时用数值雅可比也要保证接口一致。第二协方差矩阵在数值运算中容易失去对称性甚至非正定。我加了对称化处理和特征值截断这是一种“外科手术式”的补救措施。更体系化的做法是使用Joseph形式更新协方差但在仿真里特征值截断已经足够而且速度更快。第三观测雅可比 (H) 的计算。如果观测是直接相对位置和相对速度那么 (H) 就是两个稀疏块矩阵如果观测是视线角加距离则要仔细算偏导。我建议先从全状态直接量测开始调试等滤波收敛正常后再换更真实的观测模型。5.3 博弈求解函数实现博弈求解我采用最简单的梯度迭代法。基本思想是每个控制向量沿自身收益梯度方向调整但双方交替进行近似逼近均衡点。function [u_p, u_e] solve_epsilonNE(X_est, params) % 初始化策略取上一时刻策略为初值 u_p params.u_p_prev; u_e params.u_e_prev; for k 1:params.max_iter_gd % 追踪方梯度下降min J_p, 所以沿负梯度方向更新 grad_p compute_grad_Jp(u_p, u_e, X_est, params); u_p_new u_p - params.alpha_p * grad_p; u_p_new project_control(u_p_new, params.a_p_max_est); % 逃逸方梯度上升max J_e, 所以沿正梯度方向更新 grad_e compute_grad_Je(u_p_new, u_e, X_est, params); u_e_new u_e params.alpha_e * grad_e; u_e_new project_control(u_e_new, params.a_e_max_est); % 判断收敛连续两次策略差足够小 if norm(u_p_new - u_p) norm(u_e_new - u_e) params.tol_gd u_p u_p_new; u_e u_e_new; break; end u_p u_p_new; u_e u_e_new; end end这个函数结构清晰但也需要说明几个点。一是梯度计算。compute_grad_Jp里如果把收益函数看作一次单步预测的形式可以用解析梯度直接算如果用多步前向Rollout就要用伴随法或数值差分。末端长度10秒以内时我直接对控制轨迹做数值差分也能接受末端长度超过30秒时数值差分的计算量会爆炸必须换伴随法。二是投影算子project_control。它的作用是约束推力加速度不超过当前EKF估计的最大值。这里特别要注意投影后的控制量必须再用于下一轮博弈迭代不能把投影放到迭代外只做最后一步否则迭代可能来回震荡不收敛。我是在迭代循环内部每次更新后都做投影。三是学习率 (\alpha_p)、(\alpha_e) 的选取。这里没有通用公式我建议从0.1开始逐步调整。如果策略发散就减小到0.05如果收敛太慢就试0.2。还有一个更稳妥的做法是用自适应步长比如梯度方向的步长取收益改进量的线性搜索结果但代码复杂度会高不少。5.4 主循环与结果输出主循环写起来很直白% run_simulation.m 主循环片段 for k 1:length(t)-1 % 1. 生成真实带噪测量 z measurement_true(X_true(:,k)) mvnrnd(zeros(6,1), R); % 2. EKF估计 [X_est(:,k1), P] ekf_update(X_est(:,k), P, z, u_p(:,k), u_e(:,k), dt, Q, R); % 3. 利用估计值求解博弈 params.X_est X_est(:,k1); params.u_p_prev u_p(:,k); params.u_e_prev u_e(:,k); [u_p(:,k1), u_e(:,k1)] solve_epsilonNE(params); % 4. 真实动力学推进 X_true(:,k1) dynamics_true(X_true(:,k), u_p(:,k1), u_e(:,k1), dt); end仿真结果重点看四个图相对轨迹图看追踪方是否成功逼近逃逸方是否成功逃开距离-时间曲线看末端最小距离是否达到设计目标EKF参数估计曲线看 (a_{p,max}) 和 (a_{e,max}) 的估计值是否收敛到真实值附近ε收敛曲线看博弈迭代次数和残差是否随时间趋于稳定。我的复现经验是如果EKF参数估计收敛良好距离曲线也会比较平滑如果距离曲线频繁出现异常跳变十有八九是博弈求解器振荡先检查步长和投影逻辑大概率能解决。6. 常见问题与排查技巧实录6.1 问题速查表现象可能原因解决方法EKF估计发散状态直接飞出合理范围初始协方差P0设置过小滤波器过度信任初值把P0设置到合理物理量级比如位置方差取100²速度方差取10²参数估计收敛到错误值参数不可观测或激励不足调整观测模型或在博弈轨迹里加入小幅探测激励博弈迭代不收敛策略来回震荡学习率过大、投影顺序不对降低学习率确保每次迭代内投影后再参与下一步迭代仿真前期控制量抖动剧烈EKF估计初期协方差大估计噪声直接进入博弈求解对博弈求解输出的控制指令做速率限制或低通滤波后期控制量长期不变化估计过程噪声Q设得过小参数锁死适当增大参数对应Q给估计留一点“呼吸空间”相对距离曲线出现毛刺传感器野值未处理添加野值剔除逻辑比如残差超过3倍标准差时降权处理6.2 三个值得展开的避坑经验第一个经验是关于EKF参数可观测性的。我初版仿真里试图从相对位置观测中估计双方最大推力加速度但跑了500组参数发现逃逸方参数估计经常跑偏。后来一查原来逃逸方只对“拉距离”有影响而追踪方一直紧追使得相对距离变化缓慢逃逸方推力始终未被充分激发参数自然不可观测。解决方案是在博弈收益函数里加入一个小的探针激励项让逃逸方在博弈框架下“被迫”间歇性拉开距离从而让观测序列包含参数信息。这个过程听起来有点绕但本质上就是系统辨识里的“持续激励”条件。第二个经验是关于控制抖动与滤波正反馈的。有段时间仿真总在末端10秒左右出现距离曲线剧烈波动我最初以为是博弈求解器出问题。后来把EKF估计的控制输入和真实控制输入画在同一张图上才发现问题博弈策略中包含了估计误差带来的高频噪声这些噪声又通过动力学进入下一次量测滤波器误以为状态真的高频变化于是估计协方差增大、参数估计也随之偏离。这个“控制-估计相互污染”的问题在工程上非常典型。解决起来也不难我在博弈求解器输出端加了一个一阶惯性滤波[ u_{applied} \beta u_{desired} (1 - \beta) u_{applied, prev} ](\beta) 取0.3到0.5之间即可既抑制了高频抖动又不至于让控制响应太迟钝。第三个经验是关于仿真步长的。项目里我一开始用定步长 (dt 0.1) 秒跑博弈迭代每步大概70次左右。后来想加快仿真速度把步长改成0.5秒结果EKF经常在剧烈机动段发散。原因是动力学模型在这段时间内明显非线性和强耦合线性化假设在步长过大时不再成立。我后来改成自适应步长状态变化率小时用0.2秒接近博弈末端时自动切到0.02秒。虽然代码多了一点分支逻辑但仿真稳定性肉眼可见地上升。如果你不想写自适应步长至少也要保证在距离小于300米的“末端末端”把步长手动调小。6.3 关于复现期刊代码的心态建议最后说点软性的东西。这类“期刊复现”项目最大的陷阱不是算法看不懂而是你以为自己复现成功了其实没有。我调试过程中反复遇到“这次结果比论文好”“下次结果比论文差”的现象后来明白这纯粹是参数初值不同造成的差异而不是算法有本质优劣。正确做法是把论文里的核心结论作为一个边界条件只要你的仿真在合理参数范围内能复现出“估计收敛、博弈接近ε均衡、控制闭环稳定”这三个特征就可以认为复现成功。至于绝对值是否和论文图表完全一致反而不重要。7. 一点个人实践体会这个项目我前后折腾了两周最深的感受是真正困难的不是EKF也不是博弈求解而是把两者在同一个闭环里稳定地跑起来。单独测试EKF收敛很容易——给一组标准轨迹滤波效果一目了然单独测试博弈求解也很容易——固定参数下迭代几十轮看结果就行。但它们一结合各种正反馈、参数耦合、时序关系问题就全都冒出来了。如果让我给后来者排一个合理的复现顺序我会建议这样走先做纯完全信息的博弈仿真确认策略本身合理再加入EKF做状态估计但不做参数辨识确认滤波和控制可以共存最后才把参数辨识加进来实现完整的不完全信息博弈。这个顺序虽然多花一点时间但每一步的变量都很少出问题能快速定位。项目里还有一个值得继续扩展的方向目前EKF假设噪声都是高斯白噪声但真实相对导航在末端可能会遇到非高斯厚尾噪声或通信断续。后续可以考虑用鲁棒卡尔曼滤波或粒子滤波替代EKF看看在这种更恶劣的噪声环境下ε-纳什均衡策略还能保持多少性能。这个方向如果再发一篇论文应该不难产出而且工作量集中在滤波部分博弈部分基本可以复用。如果时间有限优先把ε的取值敏感性分析做一遍这个结果对你理解整体策略的鲁棒性格外有价值。