简介：基于MATLAB/Simulink的卫星避碰方案仿真工程包，面向航天器轨道设计、任务规划与仿真验证方向的科研学习者及工程师，用于解决低轨卫星数量增长带来的碰撞风险建模与规避机动决策问题。压缩包共8个文件，其中5个.m脚本构成核心代码，覆盖轨道动力学建模、避碰检测算法与机动策略计算；2个.txt文档和1个README.md提供运行说明、参数设置与使用流程，帮助理解程序逻辑。整个资源包仅4KB，轻量紧凑，但设计链条完整，基于开普勒定律与牛顿运动方程，并考虑了地球非球形引力及日月摄动等影响因素，从实时监控相对位置、速度到生成规避策略均可仿真验证。内容预览显示项目结构包含Satellite-Collision-Avoidance主目录、scripts脚本子目录与data数据目录，便于在Simulink中直接打开、调整参数并观察效果。该资源已有53人学习下载，适合作为卫星避碰算法课程设计、科研预研或入门验证的参考模板。
1. 卫星避碰方案：Simulink能替你把关的每个决策节点
在轨航天器收到一份碰撞预警之后，真正要做的事只有两件：把未来 24 小时的接近弧段算清楚，在几个可行机动方案里选出代价最小的那一个。但这两件事中间隔着至少三道计算——轨道外推、TCA（最接近时间）求解、碰撞概率积分，任何一步参数错了，机动不做是风险，做了是燃料白烧。这套基于 MATLAB/Simulink 的卫星避碰仿真方案，就是把整条链路做成一个可追溯、能改参、能批量跑仿真的工程骨架。它适合做轨道力学课程设计或毕业设计的同学，也适合在轨任务工程师在做复算复核时当参照模板，核心价值是让你不用再从零攒模型。
2. 轨道外推与近接筛选：TCA 算不准，后面全是徒劳
避碰仿真的第一层不是 Simulink 模型，而是你喂给模型的轨道状态和力模型。很多第一次做避碰的人上来就搭模型，结果初始轨道根数写错、摄动力开关忘开，后面所有结论都站不住。这里先把两个底层问题讲透。
2.1 力模型怎么选：二体打底，J2 开关必须显式化
二体模型的好处是闭合解，速度快，适合粗扫和参数初始化；但它忽略了地球扁率的影响。J2 摄动会引起升交点赤经的长期漂移，在 500 km 高度太阳同步轨道上量级大约是一天漂移 0.5°～1°，24 小时外推下来终点位置差几十公里很正常。避碰场景本来就要求相对位置精度到公里级甚至百米级，所以模型里必须留一个 J2 开关，让用户能对比“开与不开”的差异。
工程上一般这样组织力模型：一个 MATLAB Function 接收卫星位置矢量，返回加速度矢量，加速度写成“中心引力 + J2 修正项”两部分。中心引力项就是经典的两体加速度，J2 修正项的表达式固定，但系数乘一个 mask 参数enable_J2。这样在快速方案筛选时可以关掉 J2 跑初值，在正式出结论时必须打开。
| 轨道高度 | J2 引起的交点漂移量级 | 24 h 位置误差影响 |
|---|---|---|
| 500 km | 约 0.8°～1°/天 | 几十 km 量级 |
| 800 km | 约 0.5°/天 | 十几 km 量级 |
| 1200 km | 约 0.3°/天 | 几 km 量级 |
这张表想说明的核心观点是：轨道越低，J2 开关越不能省。如果你做的是 LEO 避碰，把一个关掉 J2 的模型拿来出结论，基本等于白算。
模型里的轨道初值也不要手填六个不统一的量。常见做法是在 MATLAB 工作区里定义一套结构体，半长轴、偏心率、倾角、升交点赤经、近地点幅角、真近点角各一字段，由初始化脚本统一换算成 ECI 位置速度后，再喂给积分器。这个流程虽然多写几行，但能避免很多低级错误。
2.2 TCA 计算：最短接近时刻不是看曲线最低点
TCA 是避碰仿真的核心产出之一，多余的“最近距离”和“碰撞概率”都要挂在这个时刻上。最容易犯的错是用大仿真步长跑完，然后趴在曲线图上肉眼找最低点，再拿那个点当 TCA。这个做法有两个问题：一是步长 10 秒或 30 秒，找出来的点天然带半个步长的误差；二是相对距离变化快的时候，曲线最低点附近非常陡，肉眼判读误差能到几公里。
正确做法是两步走：先用大步长粗扫定位到 TCA 附近的小区间，再用二分法或牛顿迭代在小区间内精细收敛。下面这段代码是我在方案里实际使用的 TCA 求解函数。
function tca = find_tca(r1, v1, r2, v2, t0, t1, dt) % 输入：t0 时刻两星位置(m)与速度(m/s)，搜索窗口[t0,t1]与粗扫步长dt t = t0:dt:t1; d2 = zeros(size(t)); for k = 1:numel(t) d2(k) = dist2_at(r1, v1, r2, v2, t(k) - t0); end [~, idx] = min(d2); % 粗扫最小下标 % 以相邻两点夹住最小值，进入二分收敛 tL = t(max(idx-1, 1)); tR = t(min(idx+1, numel(t))); for k = 1:60 tm = 0.5 * (tL + tR); if dist2_at(r1, v1, r2, v2, tm - t0) > dist2_at(r1, v1, r2, v2, tL - t0) tR = tm; else tL = tm; end end tca = 0.5 * (tL + tR); end function d2 = dist2_at(r1, v1, r2, v2, tau) % 短弧段内使用匀速直线近似外推相对位置 dr = (r1 + v1 * tau) - (r2 + v2 * tau); d2 = dr' * dr; end这段代码的适用范围是短弧段场景。避碰预警窗口一般只有几分钟到几十分钟，相对速度在一两公里每秒量级，加速度引起的弯曲效应在这个时间尺度下很小，所以r + v * tau的线性近似足够用。二分 60 次可以让 TCA 收敛到毫秒级，dt 取 300 秒做粗扫也不会漏峰。
但必须说清楚边界：如果是大偏心率轨道，或者相对距离本身就是几十公里以上、接近弧段跨越多个轨道周期，线性外推的假设就不成立了。这时候要把dist2_at里的线性外推替换成数值积分结果，函数签名不变，只改内部实现即可。
2.3 碰撞概率：黑匣子里的那个数怎么算出来的
TCA 给了我们最近接近时刻，但“最近距离”本身不是完整判据，因为定轨误差的存在让位置带有不确定性。碰撞概率做的事是把两个星的联合位置误差投影到碰撞平面上，再对碰撞半径圆域做二维高斯积分。对很多人来说，这个概率输出像个黑匣子，但它的输入其实就四个东西：TCA 时刻相对位置、相对速度、联合协方差矩阵、碰撞半径。
联合协方差矩阵直接由两星定轨协方差相加得到，因为两星观测相互独立。碰撞半径取两星包络半径之和，工程上还会再加一个安全裕度，比如各加 50 米。概率公式本身不复杂：
[ P_c = \frac{1}{2\pi \sqrt{|C|}} \iint_{x^2 + y^2 \le R_c^2} \exp\left(-\frac{1}{2} (\mathbf{r}-\boldsymbol{\mu})^T C^{-1} (\mathbf{r}-\boldsymbol{\mu})\right) dx,dy ]
这里的 (\mathbf{r}) 是碰撞平面上的相对位置矢量，(C) 是联合协方差矩阵在碰撞平面上的二维投影。常见的简化处理是忽略积分区间内概率密度的剧烈变化，把二维高斯积分用网格求和逼近，网格边长为碰撞半径的 1/20 到 1/50，精度已经足够。
| 输入参数 | 含义 | 典型取值/来源 |
|---|---|---|
| r_rel | TCA 时刻相对位置 | 轨道外推输出 |
| v_rel | TCA 时刻相对速度 | 轨道外推输出 |
| C | 联合协方差矩阵投影 | 两星定轨协方差相加 |
| R_c | 碰撞半径 | 两星包络半径之和 + 裕度 |
这一步是“要不要机动”的最终依据。很多情况下两星距离只有几百米，但协方差很小，概率反而不高；另一些情况距离几公里，协方差很大，概率反而超过阈值。这也是为什么避碰决策不能只看距离，必须落到概率上。
3. 机动决策逻辑：CW 方程与 MATLAB Function 的落地边界
算完碰撞概率，模型进入决策环节。如果概率没超阈值，保持当前轨道即可；如果超了阈值，就要回答“往哪个方向推、推多少”。机动规划的动力学基础是 Clohessy-Wiltshire 方程，也就是常说的 CW 方程。
3.1 CW 方程的适用边界
CW 方程把目标星轨道作为圆参考轨道，在 Hill 坐标系下描述追踪星相对于目标星的线性化运动。x 轴沿径向向外，y 轴沿迹向，z 轴沿轨道面法向。这个方程有两个硬前提：参考轨道近圆，相对距离远小于轨道半径。避碰场景恰好落在这个范围内——轨道偏心率一般小于 0.01，接近距离从几百米到几十公里，相对轨道半径可以忽略。所以用 CW 方程做快速机动方案筛选是合理的，没必要在决策阶段就上完整非线性轨道积分。
CW 方程给出了三种典型机动的不同响应特征，这是选方向的核心依据。
| 机动方向 | 相对运动短期效果 | 24 小时后趋势 | 燃料效率 |
|---|---|---|---|
| 径向 | 改变相对位置径向分量 | 周期性调制，无长期漂移 | 中 |
| 迹向 | 改变相对半长轴，产生线性漂移 | 漂移持续累积，分离效果明显 | 最高 |
| 法向 | 改变相对轨道面 | 小幅周期振荡，几乎无长期分离 | 低 |
所以工程上最常用的是沿迹向机动。迹向脉冲产生的相对迹向漂移速度大约是 3 倍脉冲量级，一个 1 cm/s 的迹向脉冲，24 小时后能拉开约 2.5～3 km 的相对距离。径向和法向机动想要达到同等分离效果，需要明显更大的速度增量，而且径向机动的相位响应会随时间变化，实际使用中还要反复评估时序，复杂度高不少。
3.2 MATLAB Function 块里的决策逻辑封装
Simulink 里实现决策逻辑，最简洁的方式是 MATLAB Function 块，不需要额外装工具包。输入端接相对位置、相对速度、碰撞概率和阈值参数，输出端给机动脉冲矢量与模式标记。核心逻辑是先判断概率是否越线，再对候选方案做燃料最小选择。
function [dv_vec, mode] = maneuver_logic(r_rel, v_rel, Pc, P_th) % 输入端：相对位置(m)，相对速度(m/s)，当前碰撞概率，概率阈值 % 输出端：建议脉冲(m/s)，机动方向编码 0/1/2/3 dv_vec = zeros(3, 1); mode = 0; if Pc < P_th return; % 概率未越线，不机动 end % 三方向候选脉冲：迹向 1cm/s 作为基线，径向/法向按比例放大 dv_base = 0.01; candidates = [ 0.5*dv_base, 0, 0; % 径向候选 0, dv_base, 0; % 迹向候选 0, 0, 0.3*dv_base; % 法向候选 ]; % 按脉冲模长最小原则选方案 [~, mode] = min(sum(candidates.^2, 2)); dv_vec = candidates(mode, :)'; end这个函数有两点要说明。第一，P_th阈值不是随便拍的。工程上低轨避碰概率阈值通常取 (1\times10^{-4}) 到 (1\times10^{-6}) 之间，取决于任务风险等级和机动能力。课程设计建议用 (1\times10^{-4})，更容易看清整条链路的效果。第二，候选脉冲的幅值不是最终答案，它只是一个初值。真实流程是把候选脉冲带回完整轨道外推模型，看 TCA 时刻碰撞概率是否降到阈值以下，不满足就放大脉冲再来一轮。Simulink 模型里建议把决策函数和轨道积分器做成闭环，决策输出接到积分器输入，让模型自动迭代出满足条件的最小脉冲。
3.3 要不要上 Chart 状态机
如果避碰流程还要分阶段处理，比如“预警 → 评估 → 决策 → 执行 → 确认”五步，中间夹杂人工确认标志和超时计时，可以考虑用 Simulink 的 Chart 做状态机，把 MATLAB Function 当成里面的一个 action 函数调用。但大多数教学和验证场景用不到这么重的手段，一个 MATLAB Function 块就够了。状态机的额外收益是流程可视化，代价是增加状态同步的调试成本，我一般只在需要严格时序控制时才引入。
4. 可追溯的仿真工程：从模块接线到批量跑参
前面三章把单点算法讲清楚了，这一章说怎么把它们组装成一个能反复改参、能批量跑仿真的 Simulink 模型。这个工程骨架本身也是这套资源的核心交付物。
4.1 模型分层与模块职责
顶层模型不要拍平了画，按功能拆成四个子系统：轨道外推子系统、近接检测子系统、机动决策子系统、结果记录。轨道外推子系统内部是两个并行的积分器，分别积分两星状态，加速度输入来自力模型函数；近接检测子系统每个仿真步计算一次相对距离，并输出到工作区；决策子系统读入碰撞概率，输出脉冲；结果记录用 To Workspace 模块把距离序列、概率序列、速度增量序列统一落盘。
脉冲注入积分器的方式要特别注意。常见做法是把脉冲折算成有限推力弧段，在决策触发时刻往积分器的加速度输入端加一个持续几十秒的推力。这样做数值上是连续的，不会出现状态突变；直接在积分器输出端改状态值虽然简单，但容易触发代数环或者丢事件。
4.2 仿真参数配置参考
Simulink 求解器的选择直接决定结果可复现性。我的习惯是固定步长 ode4，步长 0.1 秒到 1 秒之间，外推 24 小时用 1 秒步长即可，50 万步在普通笔记本上几分钟跑完。如果要看概率积分的稳定性，可以把步长缩到 0.1 秒对比一下，结果差异在 1% 以内说明步长足够。
| 配置项 | 推荐值 | 原因 |
|---|---|---|
| 求解器 | ode4 固定步长 | 结果可复现，不随求解器版本漂移 |
| 步长 | 0.1～1 s | 平衡精度与仿真时长 |
| 外推时长 | 24～72 h | 覆盖典型避碰窗口 |
| 结果采样 | 10～30 s | 控制落盘数据量 |
| 停止时间 | 86400 s | 24 小时 |
模型参数不要散落在各个模块里。推荐用模型回调函数InitFcn统一初始化，所有轨道根数和力模型参数都从工作区结构体读取。
function InitFcn_Callback() % 模型初始化：统一维护轨道根数与力模型开关 global SC SC.mu = 3.986004418e14; % 地球引力常数 m^3/s^2 SC.Re = 6378.137e3; % 地球半径 m SC.J2 = 1.08262668e-3; % J2 摄动系数 SC.a = 6878.137e3; % 目标星半长轴，500km 圆轨道 SC.e = 1e-4; % 近圆轨道小偏心率 SC.i = 97.4 * pi / 180; % 倾角 SC.raan = 90 * pi / 180; % 升交点赤经 SC.argp = 0; % 近地点幅角 SC.nu0 = 0; % 初始真近点角 % 追踪星初始状态由轨道机动场景单独赋值 end这样写的收益很直接：改轨道初值不用进模型内部翻模块参数，改完工作区变量重新初始化即可。所有仿真配置集中在一个地方，后续出问题也好回溯。
4.3 用脚本批量跑参，不一个个点仿真
避碰方案验证不能只跑单场景。三个机动方向、三档脉冲幅值、两档阈值，组合下来就是十几个算例，手动点仿真按钮会把人耗死。正确做法是用 Simulink.SimulationInput 对象做批量离线仿真。
clear; load_system('sat_collision_avoid_model'); scenarios = struct(); scenarios(1).dv = [0; 0.01; 0]; % 迹向 +1cm/s scenarios(2).dv = [0.005; 0; 0]; % 径向 +0.5cm/s scenarios(3).dv = [0; 0; 0.005]; % 法向 +0.5cm/s for i = 1:numel(scenarios) simIn(i) = Simulink.SimulationInput('sat_collision_avoid_model'); simIn(i) = simIn(i).setVariable('dv_burn', scenarios(i).dv); simIn(i) = simIn(i).setVariable('burn_time', 30); % 推力持续30秒 end out = sim(simIn, 'StopTime', '86400'); for i = 1:numel(out) range_sig = out(i).logsout.get('range_m').Values.Data; [min_dist(i), ~] = min(range_sig); fprintf('方案 %d: 24h 最小距离 = %.3f km\n', i, min_dist(i) / 1000); end这段脚本的核心是setVariable在仿真前动态覆盖模型工作区变量，不需要手动打开模型改参数。sim函数支持向量化输入，多个方案一次性提交，结果按顺序返回。跑出来后从logsout里提取相对距离序列，统计 24 小时内最小值，就能横向比较三个方向的机动效果。
这里有个小提醒：dv_burn必须和模型里的变量名完全一致，大小写都不能错，否则setVariable不会报错但变量值不会被用到，最后结果全是一个方案的重复。检查办法是跑一个算例后把dv_burn打出来核对，别嫌这一步麻烦。
5. 避坑指南：卫星避碰仿真最容易翻车的五个点
资源我用过、也看别人翻车过，下面这五条是出现频率最高的坑，每一条都对应着实际仿真里的具体报错或异常结果。
5.1 S-Function Builder 编译后，其他 S-Function 也被带挂了
现象：模型里有多个 S-Function Builder 块，编译其中一个后，其他块的输出变成旧值或者直接报维度错误。 原因：S-Function Builder 生成的 C MEX 文件在 build 时可能把共享的模型头文件覆盖或重写，导致其他块的接口定义失效。 解决：不要在同一个模型里混用多个 S-Function Builder。能用 MATLAB Function 块实现的逻辑全部用 MATLAB Function 块替代，需要 C 代码的场景单独建一个模型验证，验证完再集成到主模型。如果确实必须用多个，每次重新生成后对所有 S-Function 块统一执行一次 rebuild。
5.2 机动脉冲加进模型了，但仿真结果毫无变化
现象：决策模块输出了明显的速度增量，目标星轨道却纹丝不动。 原因：脉冲是离散事件，但被当作普通连续信号接在了积分器输入端。如果推力持续时间为零，连续求解器根本不会在事件时刻做积分，脉冲被直接跳过。 解决：把脉冲折算成有限推力弧段，持续 30 秒以上再输入积分器。或者在模型中用 State 重置事件触发逻辑，让求解器在脉冲注入时刻强制停下。课程设计里用有限推力弧段最简单，既可以避免事件问题，还能顺便算燃料消耗。
5.3 TCA 迭代不收敛或收敛到窗口外
现象：二分法迭代出来的 TCA 落在搜索窗口边缘，或者概率计算结果反复跳动。 原因：粗扫步长太大，把真实最小值附近的两个采样点都漏掉了；或者相对距离函数在窗口内有两个峰值，粗扫只抓到了其中一个。 解决：第一步用 60 秒小步长做预扫，检查最小值是否在窗口中部；如果贴边，扩大窗口重新扫。第二步把二分初值改成预扫最小点前后各 2 个点，给足容差。TCA 求解这种事，初值差一点后面全白算，宁可多扫一轮。
5.4 碰撞概率异常大，大到像科幻片
现象：两星相距几十公里，碰撞概率却算出来接近 1，明显不合物理直觉。 原因：坐标系混用。TCA 外推用的 ECI 坐标，但协方差矩阵里混进了 ECEF 或其他地固系分量，或者碰撞平面法向量取的方向不对，导致概率积分把本不该重叠的区域算进去了。 解决：给所有坐标和协方差矩阵标注框架名，统一到同一惯性系；碰撞平面的法向量必须取相对速度方向，协方差矩阵投影时先做坐标旋转再截取二维分量。这条我几乎每次复核别人模型都会遇到，属于最高发坑位。
5.5 Bus Selector 没有任何可选信号
现象：连了 Bus Creator 之后，Bus Selector 的下拉列表是空的，选不到任何信号。 原因：Simulink 的总线信号是靠类型推断传播的。如果 Bus Creator 输出没有定义总线对象，或者模型编译顺序导致下游块先于上游块求类型，Bus Selector 就拿不到可用信号列表。 解决：在模型工作区显式定义 Bus 对象，把 Bus Creator 的 Output data type 设为该 Bus 对象，再连 Bus Selector。如果仍然选不到，检查是否把 Bus 信号接到了普通 Mux 或向量信号线上——这两个看起来像，实际完全不是一种东西。
6. 用蒙特卡洛外推验证：碰撞概率降下来才算数
机动方案算完，最后一步是验证。单次仿真只能说明“这个场景下有效”，但真实工程里初始定轨误差是有分布的，我们必须回答：在误差范围内，碰撞概率是否稳定降到了阈值以下。这一步最直接的做法是蒙特卡洛外推。
% 多次抽样初始位置误差并跑仿真，统计碰撞概率分布 N = 200; sigma_pos = 100; % 初始位置误差标准差 100m Pc_after = zeros(N, 1); for k = 1:N r1_err = sigma_pos * randn(3, 1); % 三轴独立抽样 simIn = Simulink.SimulationInput('sat_collision_avoid_model'); simIn = simIn.setVariable('r1_off', r1_err); out = sim(simIn, 'StopTime', '86400'); % 取仿真第四万多步处的碰撞概率 Pc_sig = out.logsout.get('Pc').Values.Data; Pc_after(k) = Pc_sig(end); end fprintf('机动后平均碰撞概率: %.2e\n安全率: %.1f%%\n', ... mean(Pc_after), 100 * sum(Pc_after < 1e-4) / N);这个脚本抽样的是追踪星初始位置误差，200 次采样后看结果分布。如果 95% 以上的算例碰撞概率都低于阈值，说明机动方案有足够鲁棒性；如果只有一半算例达标，说明脉冲幅值给小了，需要放大候选脉冲重新跑一轮。
还有一个更轻量的验证技巧：用手算一个基线算例。500 km 圆轨道上沿迹向给 1 cm/s 脉冲，理论上 24 小时后相对分离约 2.5～3 km。模型跑完输出和这个理论值对得上，说明 CW 方程的系数、单位换算和模块接线基本没问题；对不上就要回头查代码，别急着往下跑蒙特卡洛。
我做这类仿真有个习惯，会强制在模型里留一个validate_mode开关，打开后只跑理论算例并自动对比结果，对不上就报警。这个习惯救过我很多次，尤其是隔了几个月再拿旧模型时，模型里的单位换算和坐标系很可能被不小心改动过，理论算例能第一时间暴露问题。希望这套基于 MATLAB/Simulink 的避碰方案骨架也能帮你少踩几个同样的坑。
本文还有配套的精品资源，点击获取
