简介这份基于MATLAB/Simulink的卫星避碰方案代码包面向航天工程、轨道动力学及仿真技术相关方向的开发者与研究者围绕卫星在轨碰撞风险评估与规避机动策略展开。包内包括8个文件其中5个.m脚本用于构建卫星运动模型、避碰检测算法与机动决策核心逻辑2个txt文件可提供参数配置与运行说明1个md文件则对方案整体结构进行梳理压缩包整体仅4KB便于快速下载与本地部署。方案以开普勒定律和牛顿运动定律为理论基础兼顾地球非球形引力、日月摄动等影响因素结合Simulink图形化仿真环境搭建完整的避碰验证链路可直接用于学习避碰算法设计思路、修改仿真参数并观察规避效果。目前已有53人学习下载对于准备课程设计、航天仿真入门或避碰策略研究的读者这份压缩包提供了紧凑而完整的可运行示例。1. 卫星避碰不是“躲导弹”先看清Simulink方案能解决哪一段在轨卫星收到碰撞预警时工程上最纠结的不是“要不要躲”而是“往哪躲、躲多少、什么时候出手”。基于MATLAB/Simulink的卫星避碰方案解决的就是这一段决策链把相对运动外推、最接近时刻TCA估计、告警判断、机动量计算和执行验证放进同一个仿真环境让设计人员在出方案阶段就能反复调参、跑Monte Carlo。它不负责预报碰撞本身也不替代地面测控的碰撞风险评估而是把碰撞预警之后的“避让动作”变成可量化、可复现的控制逻辑。适合做星座碰撞规避策略论证、卫星控制专业课题以及想用Simulink把轨道力学落成工程模型的同行。2. 避碰算法选型先定运动模型和判据Simulink才有东西可搭拿到避碰需求最怕一上来就拖模块连信号。先问自己三个问题用什么模型外推相对位置用什么判据决定告警机动沿哪个方向给这三件事不定清楚Simulink模型搭得再漂亮结果也没法用来做决策。2.1 相对运动模型HCW方程在什么条件下够用常见做法是用HCW方程Clohessy-Wiltshire描述目标星圆轨道下的线性化相对运动。它把两星相对位置分解到轨道坐标系x为径向y为迹向z为法向运动方程为x 3n²x 2ny y -2nx z -n²z其中n为目标轨道角速度。这组方程成立的前提很苛刻目标轨道是圆轨道、两星相对距离远小于轨道半径、忽略摄动。在LEO 500 km圆轨道、相对距离几十公里的场景下HCW在几个轨道周期的预报窗口内精度够用但目标在椭圆轨道、或者预报时间超过半天就必须考虑J2摄动甚至切换成数值积分。Simulink里实现相对运动外推有两条路。一条是在MATLAB Function Block里写HCW解析状态转移矩阵速度快、可重复性好适合后面要跑几百次Monte Carlo的方案论证阶段。另一条是用Aerospace Blockset里的轨道外推模块做数值积分精度高、能带J2但仿真速度慢调一次参数可能等上几分钟。我一般建议首版用解析解等算法逻辑稳定了再上数值模型做交叉验证。2.2 碰撞判据两级告警距离阈值和概率判据分工碰撞预警不能只看相对距离这一个数。工程上常用两级判据先做粗筛选相对距离进入预设门限通常取5~10 km再进入精细判断在TCA附近计算碰撞概率。一级判据是距离门限二级判据是概率门限。Simulink里这两级可以做成一个触发式状态机也可以用MATLAB Function Block里的if-else实现。碰撞概率需要两星位置误差协方差。工程上常用定轨系统给出的3σ位置误差合成联合协方差再在TCA时刻对相对位置做高斯积分。如果拿不到协方差数据只能退回到距离阈值判据。这里有一条血泪经验阈值不能为了“安全”取得太大否则虚警率会高到逼着操作员把所有告警都忽略掉告警系统就废了。定轨误差是公里级时阈值取5 km以下就是天天报警精密定轨误差是十米级阈值才能收到1 km附近。2.3 机动方向同面追越和异面接近要分开处理同轨道面内的前后追越是LEO星座最常见场景也是最容易算的场景。这种几何下机动以迹向脉冲为主沿迹向给一个小的速度增量直接改变半长轴进而改变轨道周期让两星相位差在TCA之前被拉开。可以做一个快速估算。对圆轨道沿迹向dv导致的周期变化率约为dT/T ≈ 3dv/v相位漂移率则为3ndv/v。等待时间T_wait后相位累积量大约为Δθ ≈ 3n·T_wait·dv / v假设500 km轨道高度轨道速度约7.6 km/s角速度约0.0011 rad/s想在1000秒内拉开30 km的尾迹距离对应相位约0.0044 rad需要的dv只有3.5 m/s量级。这个量级对Monoprop推力器完全可行所以同面避碰优先选沿迹向脉冲。异面接近倾角差或升交点差主导则不同平面内机动效果不明显更有效的做法是法向脉冲直接改变轨道面方向让两星在空间上错开。但法向机动燃料代价高、点火时刻约束多属于最后手段。在Simulink里做机动决策时需要先把相对速度分解到轨道坐标系判断主导分量在迹向还是法向再决定机动方向。3. 在Simulink里搭避碰闭环模块划分与两段关键代码避碰仿真在Simulink里不是单点计算而是带反馈的决策闭环。整个模型可以按下述结构搭输入层提供两星轨道状态核心算法层做外推与决策输出层记录结果供后处理。下面给出我惯用的模块划分和两段可以直接搬进模型的代码。3.1 顶层架构输入、核心算法、输出三层分开输入层用From Workspace或From Spreadsheet载入两星的轨道状态位置、速度。如果只有TLE数据建议先在外面用MATLAB脚本做SGP4轨道外推生成一段星历文件再导入模型不要在Simulink里实时解TLE既慢又难调。核心算法层全部用MATLAB Function Block实现不依赖额外库方便跨机器拷贝。算法层包含三个子块相对运动外推块、告警判断块、机动决策块。信号流为两星绝对状态 → 转换为相对状态 → 外推得到相对距离时间序列 → 判定TCA与距离 → 输出机动标志与dv。输出层用To Workspace记录时间、相对距离、机动标志。仿真跑完后回MATLAB工作区画图比在Simulink里调Scope高效得多。顶层模型里我还会加一个脉冲发生器驱动Triggered Subsystem让决策逻辑每隔60秒评估一次而不是每个仿真步长都重算——这对后面调参和实时化都很重要。3.2 相对运动外推代码HCW解析状态转移把下面这段写进MATLAB Function Block就能实现相对状态的一步外推。输入是当前相对状态x0、步长dt和目标轨道角速度n输出是dt秒后的相对状态。function x_rel hcw_propagate(x0, dt, n) % x0: 相对状态 [x; y; z; vx; vy; vz]轨道坐标系x径向/y迹向/z法向 % dt: 外推步长, s % n: 目标轨道角速度, rad/s tau n * dt; Phi zeros(6,6); Phi(1,1) 4 - 3*cos(tau); Phi(1,2) sin(tau); Phi(1,4) sin(tau)/n; Phi(1,5) 2*(1-cos(tau))/n; Phi(2,1) 6*(sin(tau)-tau); Phi(2,2) 1; Phi(2,4) 2*(cos(tau)-1)/n; Phi(2,5) (4*sin(tau)-3*tau)/n; Phi(3,3) cos(tau); Phi(3,6) sin(tau)/n; Phi(4,1) 3*n*sin(tau); Phi(4,2) n*cos(tau); Phi(4,4) cos(tau); Phi(4,5) 2*sin(tau); Phi(5,1) 6*n*(cos(tau)-1); Phi(5,2) -n*sin(tau); Phi(5,4) -2*sin(tau); Phi(5,5) 4*cos(tau)-3; Phi(6,3) -n*sin(tau); Phi(6,6) cos(tau); x_rel Phi * x0; end这段代码的核心是用状态转移矩阵Phi直接做一步解析外推。tau n*dt是归一化时间把轨道周期折叠进去这样不管轨道高度是300 km还是800 km步长的相对含义一致。MATLAB Function Block里建议把输入输出都声明成double型避免默认的int32截断把亚米级位置误差吃掉。注意HCW模型的前提是近距离、圆轨道如果相对距离超过100 km或者轨道偏心率大于0.01这个外推结果就该谨慎使用了。3.3 告警和机动决策代码输出TCA、最小距离和机动矢量外推只能给出相对运动轨迹真正做决策还需要一个函数来扫描整个预报窗口找到最接近时刻并计算机动量。下面这段代码可以直接放进另一个MATLAB Function Block它会返回TCA、最小相对距离、告警标志和需要的dv。function [TCA, rmin, flag, dv] collision_assessment(x0, t0, tmax, n, r_safe, t_maneuver) % x0: 当前相对状态 [x; y; z; vx; vy; vz], 单位 m, m/s % t0: 当前仿真时刻, s % tmax: 预报长度, s % n: 目标轨道角速度, rad/s % r_safe: 安全距离阈值, m % t_maneuver: 预期机动点到TCA之间的时间, s dt 1; % 外推步长固定1秒LEO场景足够 t t0:dt:t0tmax; x zeros(6, length(t)); x(:,1) x0; for k 1:length(t)-1 x(:,k1) hcw_propagate(x(:,k), dt, n); end r sqrt(sum(x(1:3,:).^2, 1)); [rmin, idx] min(r); TCA t(idx); flag rmin r_safe; % 只有在TCA还有足够执行时间时才考虑机动 dv [0; 0; 0]; if flag (TCA - t0) 120 dv plan_avoidance(x0, n, r_safe, t_maneuver); end end function dv plan_avoidance(x0, n, r_safe, t_maneuver) % 简化的沿迹向脉冲估计通过改变周期拉开相位 v_orbit 7600; % 500km轨道速度近似值, m/s phase_req max(2*r_safe / (norm(x0(1:3)) eps), 0.01); dv_y v_orbit * phase_req / (3 * n * t_maneuver); dv [0; dv_y; 0]; % 沿迹向脉冲 end这个函数的执行逻辑是先用HCW解析外推扫出整个预报窗口内的相对距离序列找到最小距离和对应的TCA再判断最小距离是否突破安全阈值最后在TCA还有执行余量时计算机动量。plan_avoidance里的公式是2.3节的相位估算反推出来的属于工程快速估算不考虑推力器最小脉冲约束和姿态机动时间。实际接入卫星模型时dv还要经过推力器脉宽调制和姿态控制闭环这里只给参考输入。在Simulink里这段逻辑不要每个步长都调用。常见做法是用Triggered Subsystem包起来用脉冲发生器每60秒触发一次评估。决策输出通过Unit Delay保持到下一次触发这样既不会漏掉TCA变化也不会因为高频重算导致抖动。3.4 输出与可视化跑完数据回MATLAB画图模型里加两个To Workspace模块分别记录仿真时间t和相对距离r采样间隔设成1秒。仿真结束后回MATLAB命令窗口执行simOut sim(collision_avoidance_model.slx); r simOut.r.Data; t simOut.t.Data; plot(t/60, r/1000); xlabel(时间 (min)); ylabel(相对距离 (km)); grid on;这段脚本会把整段仿真时间内两星相对距离画成曲线。看到曲线最低点落在安全阈值线下方、且TCA前后距离呈典型抛物线形说明外推逻辑基本正常。如果曲线出现锯齿或突变优先检查相对速度的坐标系是否一致。4. 三个必调参数机动提前量、速度增量和安全阈值的边界参数标定是Simulink避碰方案里最耗时间的部分。三个参数决定整个决策逻辑的行为机动提前量、速度增量、安全距离阈值。每个都有明确的边界越界就会让结果完全失真。4.1 机动提前量留足推力器执行时间机动提前量指从决策时刻到TCA之间的时间。太短则推力器来不及完成点火太长则轨道预报误差累积机动后的轨迹可能已经偏离预期。工程上最低要求是留出姿态机动时间加推力器最长点火时间。一个典型的LEO卫星姿态机动约30~60秒推力器连续点火按任务约束可能限制在几十秒到几分钟所以TCA前10分钟锁定决策是比较稳妥的底限。更保险的做法是把评估窗口设在TCA前30分钟开始、每60秒更新一次TCA前10分钟停止更新。Simulink里用Triggered Subsystem配合Clock模块很容易模拟这种时间约束。4.2 速度增量不是越大越好dv计算不能只看“拉开距离”这一个目标还要受燃料预算和推力器能力约束。第3章的plan_avoidance函数给的是理论最小dv实际使用时需要加上至少1.5倍裕度以吸收轨道预报误差。但dv也不是越大越好过大的迹向脉冲会把卫星推到明显不同的轨道面后续恢复星座构型要花更多燃料。以500 km轨道为例迹向dv每增加1 m/s大约造成每天数十公里的相位漂移。想要在1000秒内拉开30 km3.5 m/s够用但若把dv加到10 m/s虽然拉开了近百公里卫星的轨道周期改变量也大了近3倍事后恢复构型需要反方向二次点火。这就是常见的“躲过了碰撞、丢了构型”的尴尬局面。Simulink里做参数扫描时把dv范围限制在0.5~10 m/s按0.5 m/s步进扫描观察TCA距离和轨道周期变化两条曲线。4.3 安全距离阈值跟着导航误差走安全阈值是告警逻辑的触发门限也是最容易被拍脑袋定错的值。它不应该是“感觉上安全的距离”而应该由定轨误差决定。定轨方式位置误差量级阈值建议初值TLE两行根数公里级10 kmGPS自主定轨10~30 m1~2 km星间相对测量米级300~500 m如果拿不到定轨误差数据保守做法是先取10 km跑通流程再逐步收紧观察虚警率变化。Simulink里把阈值设成Constant模块或Model Workspace变量方便Monte Carlo扫描时批量修改。5. 避坑排查Simulink卫星避碰仿真常见的5个踩坑现场这个主题的坑很多不是算法问题而是Simulink本身的操作问题。下面按“现象 → 原因 → 解决”列出5个最常见的问题每一条都是实测中真实遇到过的。5.1 Bus Selector选不出信号现象双击Bus Selector弹出的列表里一片空白或者明明连接了总线却只能看到个别信号。原因最常见的原因是总线信号没有定义Bus对象Simulink无法识别信号层级另一种情况是上游模块输出的其实是向量Mux而不是总线Bus Selector自然选不到。新版MATLAB里还有一个坑模型里用了Bus对象但没保存到数据字典换个机器打开就丢。解决如果只是简单把几路信号传进决策模块直接用Mux/Demux替代总线。如果有十几路信号需要组合传递就在模型资源管理器里显式创建Bus对象并保存到数据字典。遇到换机器丢信号的问题打开模型前先确认数据字典已加载别省这一步。5.2 MATLAB Function Block里中文注释乱码现象代码逻辑全对但中文注释变成一片问号或者乱码保存再打开依旧如此。原因MATLAB从2023版开始对UTF-8编码支持变严格而很多旧工程文件是以GBK保存的。Simulink的MATLAB Function Block内嵌代码源文件默认按系统编码读取两边不一致就出现乱码。解决在MATLAB主页预设里把字符编码设为UTF-8重新保存模型已乱码的注释直接删掉重写别试图用编码转换“救”回来转换过程会二次污染。写Simulink里的MATLAB Function代码时我后来一律用英文注释虽然难看但省事——这也是不少做仿真的人的共同习惯。5.3 多个S-Function Builder编译互相干扰现象有两个S-Function Builder模块修改其中一个并重新编译后另一个在仿真时报“找不到符号”或加载失败。原因S-Function Builder会把生成的C代码写进同一个build目录后编译的模块覆盖了先编译的产物但Simulink缓存里还记着旧符号两个模块的源文件混在一起就冲突了。解决这是Simulink的老毛病不是用法错误。简单方案是彻底清空模型相关的build目录后全量重建根治方案是放弃S-Function Builder逻辑不复杂时改用MATLAB Function Block复杂时用Simulink Coder生成独立S-Function。实测后者最稳妥虽然前期要多写一点配置。5.4 HCW外推发散距离曲线出现异常振荡现象相对距离曲线随时间变成锯齿状甚至越来越大完全不符合理论上的周期振荡规律。原因两类原因。一是步长过大HCW状态转移矩阵里的三角函数在tau接近π的倍数时数值敏感固定步长超过10秒就可能在数个轨道周期后积累出明显误差二是相对速度的单位或坐标系给错了比如把惯性系速度当成轨道系速度输入HCW方程直接失效。解决先把步长压到5秒以内跑一次如果曲线立刻恢复正常就是步长问题。如果还是振荡检查初始相对状态的符号和单位。常见错误是把径向速度写成迹向速度数值会差好几个量级。5.5 改了参数但仿真结果不变现象在模型对话框里把r_safe从10 km改成1 km重新运行告警时刻和dv输出完全不变。原因参数没有真正进模型。Simulink的MATLAB Function Block里如果引用了工作区变量变量必须存在基工作区或Model Workspace直接改模型里显示的数值改的只是默认值的副本。解决打开Model Explorer把参数定义在Model Workspace里并勾选“在模型加载时从数据字典读取”。同时检查模型的InitFcn回调确认没有初始化语句把参数重置回默认值。调完参数跑一次前在命令行执行clear all清理旧变量避免工作区残留干扰。6. 从Simulink验证到STK校核把避碰方案从“能跑”做到“可信”Simulink模型跑通只是第一步离“可信”还差两级验证。第一级是跟STK等高精度轨道工具对照第二级是Monte Carlo统计。STK对照的常见做法是把Simulink算出的机动dv写成STK的Maneuver命令在STK里用高精度轨道预报HPOP重新推演一遍对比两边的TCA时刻和最小距离。通过标准可以按相对精度定TCA时刻偏差小于1分钟最小距离偏差小于10%。如果偏差超限优先怀疑HCW线性化误差而不是STK的问题。Monte Carlo验证是把导航误差注入初值批量跑200次以上统计漏警率和虚警率。Simulink里用parfor并行跑for i 1:200 x0_noisy x0 randn(6,1) .* sigma; simOut sim(collision_avoidance_model.slx); % 记录是否告警结束 end校验对象方法建议通过标准相对运动模型STK对照HCW外推最小距离偏差10%机动策略STK重放机动指令TCA距离安全阈值闭环整体Monte Carlo 200次漏警率0虚警率5%我习惯在交付任何避碰策略前先把Monte Carlo跑一遍再挑一个“几乎要触发告警”的场景手动复盘。有一次就是Monte Carlo里有个样本在TCA前3分钟才突破阈值虽然没漏警但留给推力器的时间只剩3分钟工程上根本来不及执行。后来我把决策锁定时间从10分钟改到15分钟才把这个隐患堵上。这类时间余量问题纯看单次仿真根本发现不了必须靠批量统计暴露。希望帮到你。本文还有配套的精品资源点击获取 卫星避碰仿真方案：基于MATLAB/Simulink的HCW建模与机动决策
简介这份基于MATLAB/Simulink的卫星避碰方案代码包面向航天工程、轨道动力学及仿真技术相关方向的开发者与研究者围绕卫星在轨碰撞风险评估与规避机动策略展开。包内包括8个文件其中5个.m脚本用于构建卫星运动模型、避碰检测算法与机动决策核心逻辑2个txt文件可提供参数配置与运行说明1个md文件则对方案整体结构进行梳理压缩包整体仅4KB便于快速下载与本地部署。方案以开普勒定律和牛顿运动定律为理论基础兼顾地球非球形引力、日月摄动等影响因素结合Simulink图形化仿真环境搭建完整的避碰验证链路可直接用于学习避碰算法设计思路、修改仿真参数并观察规避效果。目前已有53人学习下载对于准备课程设计、航天仿真入门或避碰策略研究的读者这份压缩包提供了紧凑而完整的可运行示例。1. 卫星避碰不是“躲导弹”先看清Simulink方案能解决哪一段在轨卫星收到碰撞预警时工程上最纠结的不是“要不要躲”而是“往哪躲、躲多少、什么时候出手”。基于MATLAB/Simulink的卫星避碰方案解决的就是这一段决策链把相对运动外推、最接近时刻TCA估计、告警判断、机动量计算和执行验证放进同一个仿真环境让设计人员在出方案阶段就能反复调参、跑Monte Carlo。它不负责预报碰撞本身也不替代地面测控的碰撞风险评估而是把碰撞预警之后的“避让动作”变成可量化、可复现的控制逻辑。适合做星座碰撞规避策略论证、卫星控制专业课题以及想用Simulink把轨道力学落成工程模型的同行。2. 避碰算法选型先定运动模型和判据Simulink才有东西可搭拿到避碰需求最怕一上来就拖模块连信号。先问自己三个问题用什么模型外推相对位置用什么判据决定告警机动沿哪个方向给这三件事不定清楚Simulink模型搭得再漂亮结果也没法用来做决策。2.1 相对运动模型HCW方程在什么条件下够用常见做法是用HCW方程Clohessy-Wiltshire描述目标星圆轨道下的线性化相对运动。它把两星相对位置分解到轨道坐标系x为径向y为迹向z为法向运动方程为x 3n²x 2ny y -2nx z -n²z其中n为目标轨道角速度。这组方程成立的前提很苛刻目标轨道是圆轨道、两星相对距离远小于轨道半径、忽略摄动。在LEO 500 km圆轨道、相对距离几十公里的场景下HCW在几个轨道周期的预报窗口内精度够用但目标在椭圆轨道、或者预报时间超过半天就必须考虑J2摄动甚至切换成数值积分。Simulink里实现相对运动外推有两条路。一条是在MATLAB Function Block里写HCW解析状态转移矩阵速度快、可重复性好适合后面要跑几百次Monte Carlo的方案论证阶段。另一条是用Aerospace Blockset里的轨道外推模块做数值积分精度高、能带J2但仿真速度慢调一次参数可能等上几分钟。我一般建议首版用解析解等算法逻辑稳定了再上数值模型做交叉验证。2.2 碰撞判据两级告警距离阈值和概率判据分工碰撞预警不能只看相对距离这一个数。工程上常用两级判据先做粗筛选相对距离进入预设门限通常取5~10 km再进入精细判断在TCA附近计算碰撞概率。一级判据是距离门限二级判据是概率门限。Simulink里这两级可以做成一个触发式状态机也可以用MATLAB Function Block里的if-else实现。碰撞概率需要两星位置误差协方差。工程上常用定轨系统给出的3σ位置误差合成联合协方差再在TCA时刻对相对位置做高斯积分。如果拿不到协方差数据只能退回到距离阈值判据。这里有一条血泪经验阈值不能为了“安全”取得太大否则虚警率会高到逼着操作员把所有告警都忽略掉告警系统就废了。定轨误差是公里级时阈值取5 km以下就是天天报警精密定轨误差是十米级阈值才能收到1 km附近。2.3 机动方向同面追越和异面接近要分开处理同轨道面内的前后追越是LEO星座最常见场景也是最容易算的场景。这种几何下机动以迹向脉冲为主沿迹向给一个小的速度增量直接改变半长轴进而改变轨道周期让两星相位差在TCA之前被拉开。可以做一个快速估算。对圆轨道沿迹向dv导致的周期变化率约为dT/T ≈ 3dv/v相位漂移率则为3ndv/v。等待时间T_wait后相位累积量大约为Δθ ≈ 3n·T_wait·dv / v假设500 km轨道高度轨道速度约7.6 km/s角速度约0.0011 rad/s想在1000秒内拉开30 km的尾迹距离对应相位约0.0044 rad需要的dv只有3.5 m/s量级。这个量级对Monoprop推力器完全可行所以同面避碰优先选沿迹向脉冲。异面接近倾角差或升交点差主导则不同平面内机动效果不明显更有效的做法是法向脉冲直接改变轨道面方向让两星在空间上错开。但法向机动燃料代价高、点火时刻约束多属于最后手段。在Simulink里做机动决策时需要先把相对速度分解到轨道坐标系判断主导分量在迹向还是法向再决定机动方向。3. 在Simulink里搭避碰闭环模块划分与两段关键代码避碰仿真在Simulink里不是单点计算而是带反馈的决策闭环。整个模型可以按下述结构搭输入层提供两星轨道状态核心算法层做外推与决策输出层记录结果供后处理。下面给出我惯用的模块划分和两段可以直接搬进模型的代码。3.1 顶层架构输入、核心算法、输出三层分开输入层用From Workspace或From Spreadsheet载入两星的轨道状态位置、速度。如果只有TLE数据建议先在外面用MATLAB脚本做SGP4轨道外推生成一段星历文件再导入模型不要在Simulink里实时解TLE既慢又难调。核心算法层全部用MATLAB Function Block实现不依赖额外库方便跨机器拷贝。算法层包含三个子块相对运动外推块、告警判断块、机动决策块。信号流为两星绝对状态 → 转换为相对状态 → 外推得到相对距离时间序列 → 判定TCA与距离 → 输出机动标志与dv。输出层用To Workspace记录时间、相对距离、机动标志。仿真跑完后回MATLAB工作区画图比在Simulink里调Scope高效得多。顶层模型里我还会加一个脉冲发生器驱动Triggered Subsystem让决策逻辑每隔60秒评估一次而不是每个仿真步长都重算——这对后面调参和实时化都很重要。3.2 相对运动外推代码HCW解析状态转移把下面这段写进MATLAB Function Block就能实现相对状态的一步外推。输入是当前相对状态x0、步长dt和目标轨道角速度n输出是dt秒后的相对状态。function x_rel hcw_propagate(x0, dt, n) % x0: 相对状态 [x; y; z; vx; vy; vz]轨道坐标系x径向/y迹向/z法向 % dt: 外推步长, s % n: 目标轨道角速度, rad/s tau n * dt; Phi zeros(6,6); Phi(1,1) 4 - 3*cos(tau); Phi(1,2) sin(tau); Phi(1,4) sin(tau)/n; Phi(1,5) 2*(1-cos(tau))/n; Phi(2,1) 6*(sin(tau)-tau); Phi(2,2) 1; Phi(2,4) 2*(cos(tau)-1)/n; Phi(2,5) (4*sin(tau)-3*tau)/n; Phi(3,3) cos(tau); Phi(3,6) sin(tau)/n; Phi(4,1) 3*n*sin(tau); Phi(4,2) n*cos(tau); Phi(4,4) cos(tau); Phi(4,5) 2*sin(tau); Phi(5,1) 6*n*(cos(tau)-1); Phi(5,2) -n*sin(tau); Phi(5,4) -2*sin(tau); Phi(5,5) 4*cos(tau)-3; Phi(6,3) -n*sin(tau); Phi(6,6) cos(tau); x_rel Phi * x0; end这段代码的核心是用状态转移矩阵Phi直接做一步解析外推。tau n*dt是归一化时间把轨道周期折叠进去这样不管轨道高度是300 km还是800 km步长的相对含义一致。MATLAB Function Block里建议把输入输出都声明成double型避免默认的int32截断把亚米级位置误差吃掉。注意HCW模型的前提是近距离、圆轨道如果相对距离超过100 km或者轨道偏心率大于0.01这个外推结果就该谨慎使用了。3.3 告警和机动决策代码输出TCA、最小距离和机动矢量外推只能给出相对运动轨迹真正做决策还需要一个函数来扫描整个预报窗口找到最接近时刻并计算机动量。下面这段代码可以直接放进另一个MATLAB Function Block它会返回TCA、最小相对距离、告警标志和需要的dv。function [TCA, rmin, flag, dv] collision_assessment(x0, t0, tmax, n, r_safe, t_maneuver) % x0: 当前相对状态 [x; y; z; vx; vy; vz], 单位 m, m/s % t0: 当前仿真时刻, s % tmax: 预报长度, s % n: 目标轨道角速度, rad/s % r_safe: 安全距离阈值, m % t_maneuver: 预期机动点到TCA之间的时间, s dt 1; % 外推步长固定1秒LEO场景足够 t t0:dt:t0tmax; x zeros(6, length(t)); x(:,1) x0; for k 1:length(t)-1 x(:,k1) hcw_propagate(x(:,k), dt, n); end r sqrt(sum(x(1:3,:).^2, 1)); [rmin, idx] min(r); TCA t(idx); flag rmin r_safe; % 只有在TCA还有足够执行时间时才考虑机动 dv [0; 0; 0]; if flag (TCA - t0) 120 dv plan_avoidance(x0, n, r_safe, t_maneuver); end end function dv plan_avoidance(x0, n, r_safe, t_maneuver) % 简化的沿迹向脉冲估计通过改变周期拉开相位 v_orbit 7600; % 500km轨道速度近似值, m/s phase_req max(2*r_safe / (norm(x0(1:3)) eps), 0.01); dv_y v_orbit * phase_req / (3 * n * t_maneuver); dv [0; dv_y; 0]; % 沿迹向脉冲 end这个函数的执行逻辑是先用HCW解析外推扫出整个预报窗口内的相对距离序列找到最小距离和对应的TCA再判断最小距离是否突破安全阈值最后在TCA还有执行余量时计算机动量。plan_avoidance里的公式是2.3节的相位估算反推出来的属于工程快速估算不考虑推力器最小脉冲约束和姿态机动时间。实际接入卫星模型时dv还要经过推力器脉宽调制和姿态控制闭环这里只给参考输入。在Simulink里这段逻辑不要每个步长都调用。常见做法是用Triggered Subsystem包起来用脉冲发生器每60秒触发一次评估。决策输出通过Unit Delay保持到下一次触发这样既不会漏掉TCA变化也不会因为高频重算导致抖动。3.4 输出与可视化跑完数据回MATLAB画图模型里加两个To Workspace模块分别记录仿真时间t和相对距离r采样间隔设成1秒。仿真结束后回MATLAB命令窗口执行simOut sim(collision_avoidance_model.slx); r simOut.r.Data; t simOut.t.Data; plot(t/60, r/1000); xlabel(时间 (min)); ylabel(相对距离 (km)); grid on;这段脚本会把整段仿真时间内两星相对距离画成曲线。看到曲线最低点落在安全阈值线下方、且TCA前后距离呈典型抛物线形说明外推逻辑基本正常。如果曲线出现锯齿或突变优先检查相对速度的坐标系是否一致。4. 三个必调参数机动提前量、速度增量和安全阈值的边界参数标定是Simulink避碰方案里最耗时间的部分。三个参数决定整个决策逻辑的行为机动提前量、速度增量、安全距离阈值。每个都有明确的边界越界就会让结果完全失真。4.1 机动提前量留足推力器执行时间机动提前量指从决策时刻到TCA之间的时间。太短则推力器来不及完成点火太长则轨道预报误差累积机动后的轨迹可能已经偏离预期。工程上最低要求是留出姿态机动时间加推力器最长点火时间。一个典型的LEO卫星姿态机动约30~60秒推力器连续点火按任务约束可能限制在几十秒到几分钟所以TCA前10分钟锁定决策是比较稳妥的底限。更保险的做法是把评估窗口设在TCA前30分钟开始、每60秒更新一次TCA前10分钟停止更新。Simulink里用Triggered Subsystem配合Clock模块很容易模拟这种时间约束。4.2 速度增量不是越大越好dv计算不能只看“拉开距离”这一个目标还要受燃料预算和推力器能力约束。第3章的plan_avoidance函数给的是理论最小dv实际使用时需要加上至少1.5倍裕度以吸收轨道预报误差。但dv也不是越大越好过大的迹向脉冲会把卫星推到明显不同的轨道面后续恢复星座构型要花更多燃料。以500 km轨道为例迹向dv每增加1 m/s大约造成每天数十公里的相位漂移。想要在1000秒内拉开30 km3.5 m/s够用但若把dv加到10 m/s虽然拉开了近百公里卫星的轨道周期改变量也大了近3倍事后恢复构型需要反方向二次点火。这就是常见的“躲过了碰撞、丢了构型”的尴尬局面。Simulink里做参数扫描时把dv范围限制在0.5~10 m/s按0.5 m/s步进扫描观察TCA距离和轨道周期变化两条曲线。4.3 安全距离阈值跟着导航误差走安全阈值是告警逻辑的触发门限也是最容易被拍脑袋定错的值。它不应该是“感觉上安全的距离”而应该由定轨误差决定。定轨方式位置误差量级阈值建议初值TLE两行根数公里级10 kmGPS自主定轨10~30 m1~2 km星间相对测量米级300~500 m如果拿不到定轨误差数据保守做法是先取10 km跑通流程再逐步收紧观察虚警率变化。Simulink里把阈值设成Constant模块或Model Workspace变量方便Monte Carlo扫描时批量修改。5. 避坑排查Simulink卫星避碰仿真常见的5个踩坑现场这个主题的坑很多不是算法问题而是Simulink本身的操作问题。下面按“现象 → 原因 → 解决”列出5个最常见的问题每一条都是实测中真实遇到过的。5.1 Bus Selector选不出信号现象双击Bus Selector弹出的列表里一片空白或者明明连接了总线却只能看到个别信号。原因最常见的原因是总线信号没有定义Bus对象Simulink无法识别信号层级另一种情况是上游模块输出的其实是向量Mux而不是总线Bus Selector自然选不到。新版MATLAB里还有一个坑模型里用了Bus对象但没保存到数据字典换个机器打开就丢。解决如果只是简单把几路信号传进决策模块直接用Mux/Demux替代总线。如果有十几路信号需要组合传递就在模型资源管理器里显式创建Bus对象并保存到数据字典。遇到换机器丢信号的问题打开模型前先确认数据字典已加载别省这一步。5.2 MATLAB Function Block里中文注释乱码现象代码逻辑全对但中文注释变成一片问号或者乱码保存再打开依旧如此。原因MATLAB从2023版开始对UTF-8编码支持变严格而很多旧工程文件是以GBK保存的。Simulink的MATLAB Function Block内嵌代码源文件默认按系统编码读取两边不一致就出现乱码。解决在MATLAB主页预设里把字符编码设为UTF-8重新保存模型已乱码的注释直接删掉重写别试图用编码转换“救”回来转换过程会二次污染。写Simulink里的MATLAB Function代码时我后来一律用英文注释虽然难看但省事——这也是不少做仿真的人的共同习惯。5.3 多个S-Function Builder编译互相干扰现象有两个S-Function Builder模块修改其中一个并重新编译后另一个在仿真时报“找不到符号”或加载失败。原因S-Function Builder会把生成的C代码写进同一个build目录后编译的模块覆盖了先编译的产物但Simulink缓存里还记着旧符号两个模块的源文件混在一起就冲突了。解决这是Simulink的老毛病不是用法错误。简单方案是彻底清空模型相关的build目录后全量重建根治方案是放弃S-Function Builder逻辑不复杂时改用MATLAB Function Block复杂时用Simulink Coder生成独立S-Function。实测后者最稳妥虽然前期要多写一点配置。5.4 HCW外推发散距离曲线出现异常振荡现象相对距离曲线随时间变成锯齿状甚至越来越大完全不符合理论上的周期振荡规律。原因两类原因。一是步长过大HCW状态转移矩阵里的三角函数在tau接近π的倍数时数值敏感固定步长超过10秒就可能在数个轨道周期后积累出明显误差二是相对速度的单位或坐标系给错了比如把惯性系速度当成轨道系速度输入HCW方程直接失效。解决先把步长压到5秒以内跑一次如果曲线立刻恢复正常就是步长问题。如果还是振荡检查初始相对状态的符号和单位。常见错误是把径向速度写成迹向速度数值会差好几个量级。5.5 改了参数但仿真结果不变现象在模型对话框里把r_safe从10 km改成1 km重新运行告警时刻和dv输出完全不变。原因参数没有真正进模型。Simulink的MATLAB Function Block里如果引用了工作区变量变量必须存在基工作区或Model Workspace直接改模型里显示的数值改的只是默认值的副本。解决打开Model Explorer把参数定义在Model Workspace里并勾选“在模型加载时从数据字典读取”。同时检查模型的InitFcn回调确认没有初始化语句把参数重置回默认值。调完参数跑一次前在命令行执行clear all清理旧变量避免工作区残留干扰。6. 从Simulink验证到STK校核把避碰方案从“能跑”做到“可信”Simulink模型跑通只是第一步离“可信”还差两级验证。第一级是跟STK等高精度轨道工具对照第二级是Monte Carlo统计。STK对照的常见做法是把Simulink算出的机动dv写成STK的Maneuver命令在STK里用高精度轨道预报HPOP重新推演一遍对比两边的TCA时刻和最小距离。通过标准可以按相对精度定TCA时刻偏差小于1分钟最小距离偏差小于10%。如果偏差超限优先怀疑HCW线性化误差而不是STK的问题。Monte Carlo验证是把导航误差注入初值批量跑200次以上统计漏警率和虚警率。Simulink里用parfor并行跑for i 1:200 x0_noisy x0 randn(6,1) .* sigma; simOut sim(collision_avoidance_model.slx); % 记录是否告警结束 end校验对象方法建议通过标准相对运动模型STK对照HCW外推最小距离偏差10%机动策略STK重放机动指令TCA距离安全阈值闭环整体Monte Carlo 200次漏警率0虚警率5%我习惯在交付任何避碰策略前先把Monte Carlo跑一遍再挑一个“几乎要触发告警”的场景手动复盘。有一次就是Monte Carlo里有个样本在TCA前3分钟才突破阈值虽然没漏警但留给推力器的时间只剩3分钟工程上根本来不及执行。后来我把决策锁定时间从10分钟改到15分钟才把这个隐患堵上。这类时间余量问题纯看单次仿真根本发现不了必须靠批量统计暴露。希望帮到你。本文还有配套的精品资源点击获取