尧图网页设计 
简介基于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(idx1, 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_relTCA 时刻相对位置轨道外推输出v_relTCA 时刻相对速度轨道外推输出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.53 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.11 s平衡精度与仿真时长外推时长2472 h覆盖典型避碰窗口结果采样1030 s控制落盘数据量停止时间86400 s24 小时模型参数不要散落在各个模块里。推荐用模型回调函数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.53 km。模型跑完输出和这个理论值对得上说明 CW 方程的系数、单位换算和模块接线基本没问题对不上就要回头查代码别急着往下跑蒙特卡洛。我做这类仿真有个习惯会强制在模型里留一个validate_mode开关打开后只跑理论算例并自动对比结果对不上就报警。这个习惯救过我很多次尤其是隔了几个月再拿旧模型时模型里的单位换算和坐标系很可能被不小心改动过理论算例能第一时间暴露问题。希望这套基于 MATLAB/Simulink 的避碰方案骨架也能帮你少踩几个同样的坑。本文还有配套的精品资源点击获取 基于MATLAB/Simulink的卫星避碰仿真：轨道外推、碰撞概率与CW机动决策
简介基于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(idx1, 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_relTCA 时刻相对位置轨道外推输出v_relTCA 时刻相对速度轨道外推输出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.53 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.11 s平衡精度与仿真时长外推时长2472 h覆盖典型避碰窗口结果采样1030 s控制落盘数据量停止时间86400 s24 小时模型参数不要散落在各个模块里。推荐用模型回调函数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.53 km。模型跑完输出和这个理论值对得上说明 CW 方程的系数、单位换算和模块接线基本没问题对不上就要回头查代码别急着往下跑蒙特卡洛。我做这类仿真有个习惯会强制在模型里留一个validate_mode开关打开后只跑理论算例并自动对比结果对不上就报警。这个习惯救过我很多次尤其是隔了几个月再拿旧模型时模型里的单位换算和坐标系很可能被不小心改动过理论算例能第一时间暴露问题。希望这套基于 MATLAB/Simulink 的避碰方案骨架也能帮你少踩几个同样的坑。本文还有配套的精品资源点击获取