1. 项目背景与核心问题
航天器末端追逃博弈是空间对抗领域的关键课题,其本质是追踪方与逃逸方在有限时间内的动态策略对抗。传统研究通常假设双方完全掌握对方的动力学参数和控制策略,这种理想假设在实际任务中往往难以成立。当逃逸方通过机动变轨或信息隐藏手段使追踪方无法获取真实参数时,博弈就进入了不完全信息状态——这正是本项目要解决的核心难题。
在近地轨道环境中,航天器的相对运动通常采用Clohessy-Wiltshire(C-W)方程描述。该模型将三维空间运动解耦为径向、迹向和法向三个独立通道,其状态方程可表示为:
A = [zeros(3,3) eye(3); 3*Omega^2 0 0 0 2*Omega 0; 0 0 0 -2*Omega 0 0; 0 0 -Omega^2 0 0 0]; % 系统矩阵 B = [zeros(3,3); eye(3)]; % 控制矩阵其中Omega为轨道角速度。当逃逸方的控制矩阵B存在未知参数时,基于错误参数计算的追踪策略将导致拦截性能显著下降——仿真显示参数误差20%可使拦截时间延长50%,最终相对距离增大至15米。
2. 关键技术方案设计
2.1 Epsilon纳什均衡理论框架
针对不完全信息问题,本项目引入Epsilon纳什均衡作为性能评价标准。定义策略组合(u*,v*)满足Epsilon均衡的条件是:
|J(u*,v) - J(u*,v*)| ≤ ε, ∀v∈V |J(u,v*) - J(u*,v*)| ≤ ε, ∀u∈U其中J表示收益函数,ε为预设阈值。这意味着在有限时间内,任何单方策略偏离带来的收益提升不超过ε,实现了对完全信息纳什均衡的逼近。
2.2 基于EKF的参数估计
扩展卡尔曼滤波(EKF)是本项目的核心估计算法。通过将逃逸方的未知控制参数r扩展为系统状态,构建增广状态空间模型:
% 增广系统状态定义 X_hat = [x_pos; y_pos; z_pos; x_vel; y_vel; z_vel; r_est]; % EKF预测步骤 F = [A B*r_est; zeros(1,7)]; % 状态转移雅可比 P_pred = F*P*F' + Cov_W; % 协方差预测 % 更新步骤 K = P_pred*H'/(H*P_pred*H' + Cov_V); X_hat = X_hat + K*(z_meas - H*X_hat); P = (eye(7) - K*H)*P_pred;关键参数设置:
- 过程噪声协方差Cov_W = diag([1e-6,1e-6,1e-6,0.25e-6,0.25e-6,0.25e-6,1e10])/2
- 测量噪声协方差Cov_V = diag([1e-8,1e-8,1e-8,0.25e-8,0.25e-8,0.25e-8])/2
2.3 自适应博弈策略
基于实时参数估计,追踪方动态调整其最优控制策略:
% 实时策略计算 u = -inv(R_P)*B_est'*P_t(:,:,i)*X_hat(1:6); v = inv(R_E)*B'*P_t(:,:,i)*(X_hat(1:6)-X_E);其中P_t通过逆向求解Riccati微分方程获得,反映了当前时刻的最优控制增益。
3. MATLAB实现详解
3.1 仿真环境配置
clear all; close all; clc; Omega = 0.001; % 轨道角速度(rad/s) T = 500; % 博弈总时长(s) dt = 1; % 时间步长(s) % 初始状态设置(单位:km, km/s) X0_P = [1.5; 0.5; 0; 0; 0; 0]; % 追踪航天器 X0_E = [0; 0; 0; -0.05; 0; 0.05]; % 逃逸航天器 % 权重矩阵 R_P = eye(3)*1e6; % 追踪方控制权重 R_E = 1.5*eye(3)*1e6; % 逃逸方控制权重3.2 四阶龙格库塔求解
采用ode45求解Riccati微分方程:
P_T = Q_T; % 终端条件 options = odeset('RelTol',1e-6,'AbsTol',1e-8); sol = ode45(@(t,P) riccati_ode(t,P,A,B,R_P,R_E,Q), [T 0], P_T, options); P_t = deval(sol, 0:dt:T); % 存储时变增益其中Riccati微分方程定义为:
function dP = riccati_ode(t,P,A,B,R_P,R_E,Q) S = B*(inv(R_E) - inv(R_P))*B'; dP = -(A'*P + P*A - P*S*P + Q); end3.3 主仿真循环
for k = 1:T % EKF参数估计 [X_hat, P] = ekf_update(X_hat, P, z_meas, A, B_est, Cov_W, Cov_V); % 控制策略计算 u = -inv(R_P)*B_est'*P_t(:,:,k)*X_hat(1:6); v = inv(R_E)*B'*P_t(:,:,k)*(X_hat(1:6)-X_E); % 动力学更新 X_P = dyn_update(X_P, u, A, B, dt); X_E = dyn_update(X_E, v, A, B, dt); % 相对状态记录 X_rel(:,k) = X_P - X_E; end4. 结果分析与验证
4.1 参数估计性能
仿真显示EKF能在200秒内将参数估计误差收敛至5%以内(图1)。收敛速度受以下因素影响:
- 过程噪声协方差Cov_W:取值过大会导致估计波动,过小则降低适应性
- 测量噪声协方差Cov_V:需要与实际传感器精度匹配
- 初始协方差P:建议设置为对角阵,对角元素与各状态量级相当
关键技巧:可通过"调整遗忘因子"增强参数跟踪能力,在系统矩阵F中引入时变衰减系数0.95~0.99
4.2 拦截性能对比
| 场景 | 拦截时间(s) | 最终误差(m) | 参数误差收敛性 |
|---|---|---|---|
| 完全信息 | 320 | 0 | - |
| 无参数估计 | 480 | 15 | 持续20%误差 |
| EKF自适应策略 | 350 | 2 | 200s内收敛至5% |
4.3 典型问题排查
估计结果发散:
- 检查系统矩阵F的雅可比计算是否正确
- 验证过程噪声协方差是否过小(建议不低于1e-6)
- 确认测量更新频率与动力学步长匹配
拦截轨迹振荡:
- 调整控制权重矩阵R_P/R_E的比值(通常1:1.5~1:2)
- 检查Riccati方程求解的数值稳定性(RelTol建议1e-6)
实时性不足:
- 将龙格库塔求解离线预处理存储P_t
- 简化EKF更新步骤(可尝试UKF或无迹变换)
5. 工程实践建议
传感器配置优化:
- 相对测量至少需要位置+速度信息
- 角度测量需配合星敏感器提高精度
计算资源分配:
- EKF计算复杂度O(n^3),n为状态维数
- 7维状态单次更新约需0.3ms(i7-11800H)
鲁棒性增强:
% 添加参数估计饱和限制 r_est = max(0.5, min(2.0, r_est)); % 控制量限幅 u = sign(u).*min(abs(u), 2.0); % 加速度限制2m/s²扩展应用场景:
- 多航天器协同追逃:需设计分布式EKF架构
- 非线性动力学模型:改用无迹卡尔曼滤波(UKF)
- 通信延迟补偿:引入状态预测缓冲区