简介一个MATLAB工具包面向航天工程与轨道分析领域用于求解圆形或椭圆形地球轨道上两颗卫星之间的最近接近距离。脚本基于Brent一维最小化算法定位几何最近条件轨道传播采用Kozai方法并将地球半径增加2%以近似大气层阻力影响适合航天工程师、科研人员及空间技术学习者进行碰撞风险评估与轨道安全分析。压缩包共14个文件、230KB以12个m源码文件为主另有1个pdf说明文档和1个in输入数据示例各模块划分清晰如主控脚本、距离计算函数、最小化模块、Kozai传播函数、儒略日转换及结果打印等便于按需修改与复用调试。目前已有49人学习下载借助该工具可快速掌握卫星近距离接近的数值求解思路并依据输入参数模拟不同轨道场景为空间碎片规避、防碰撞预警与轨道设计提供实用参考。1. 两个地球卫星的最近距离先解决“时间”再谈“距离”手里拿着两颗卫星的两行根数想知道它们在接下来几个小时内会不会靠得太近这是空间态势感知里最常见的提问。两颗地球卫星之间的最近距离在英文里叫 Closest Approach防碰撞报告里经常简写成 CA。计算它不是把两个瞬时位置丢进欧氏距离公式就结束真正的难点在时间——每一颗卫星都在自己的开普勒轨道上漂移最近距离不是一个静态几何量而是一条随时间变化的相对距离曲线上的最小值。我会用 Matlab 从轨道根数重建 ECI 位置先粗扫再用优化工具箱精化把“最近距离、发生时刻、相对速度方向”一起算出来。这套流程适合写过轨道仿真但没有系统处理过角度换算和时间系统的人也适合想在任务规划里加一个独立校验脚本的工程师。2. 从轨道根数算卫星位置坐标系和开普勒方程必须先写对在把两颗卫星放到同一张图之前得先把各自的位置算对。这一步不涉及高级数学但单位、坐标系和旋转顺序三个地方错一个后面的距离计算就全白做。下面这套代码基于二体模型适合几分钟到几小时的接近筛查如果需要高精度交会评估第 4 章再补充时间系统和 J2 摄动。2.1 轨道六根数输入格式与单位是第一个坑经典轨道六根数是最常见的输入分别决定轨道的大小、形状、倾斜和当前相位。对“两颗卫星最近距离”这个问题两颗卫星必须在同一历元下定义时间轴从同一时刻起算否则算出来的 CA 时刻毫无意义。符号含义常用单位说明a半长轴km圆轨道下就是轨道半径LEO 约 6800 kme偏心率无量纲0 到 10.001 以下按小偏心率处理i轨道倾角deg相对赤道面0° 赤道轨道90° 极轨RAAN升交点赤经deg春分点到升交点的夹角argp近地点幅角deg升交点到近地点的夹角M0历元平近点角deg时间变量的起点替代真近点角作为输入单位坑集中在角度上。Matlab 的三角函数默认弧度而 TLE、星历表里给的都是角度。第一步就是把[i, RAAN, argp, M0]全部用deg2rad换算。偏心率也要检查必须满足0 e 1对偏心率接近 1 的细长轨道牛顿迭代初始值要另做处理对低轨卫星e 一般在 0.001 量级直接用平近点角当初始值就够了。2.2 用牛顿迭代解开普勒方程避免 atan2 奇异点平近点角 M 和偏近点角 E 的关系是开普勒方程M E - e * sin(E)给定 M 求 E 是超越方程工程上最常用牛顿法。迭代式是E_new E - (E - esin(E) - M) / (1 - ecos(E))分母在 e 接近 1 且 E 接近 2π 时会退化但近地轨道不会遇到。下面这个函数可以直接保存成kepler_E.mfunction E kepler_E(M, e, tol) % M 平近点角, rad % e 偏心率, 0 e 1 % tol 迭代容差, 默认 1e-10 if nargin 3, tol 1e-10; end E M; % 小偏心率轨道从平近点角出发收敛足够快 for k 1:100 dE (E - e*sin(E) - M) / (1 - e*cos(E)); E E - dE; if abs(dE) tol break; end end if k 100 warning(kepler_E:notConverged, 开普勒方程未在100步内收敛); end end逻辑说明迭代的核心是每一步计算修正量dE当修正量小于tol时停止。平近点角在传播时会持续增加调用前先做mod(M, 2*pi)避免大时间差下角度累积误差。tol1e-10对应偏近点角精度约 1e-10 rad换算成位置误差在毫米量级足够 CA 计算使用。2.3 从轨道面坐标转到 ECI旋转顺序不能记反偏近点角 E 算出后近焦点坐标系中的位置矢量为r_pf [a * (cos(E) - e); a * sqrt(1 - e^2) * sin(E); 0]要把这个矢量转到 ECI通常指 J2000 或 MOD需要依次绕 Z 轴旋转 -RAAN、绕 X 轴旋转 -i、绕 Z 轴旋转 -argp。这里的旋转方向约定来自 Curtis 的《Orbital Mechanics for Engineering Students》和 Matlab 自带的rotz方向不一定一致所以最好手写旋转矩阵function r_eci orbital2eci(a, e, i_deg, raan_deg, argp_deg, E) % 由半长轴、偏心率、倾角、升交点赤经、近地点幅角、偏近点角计算 ECI 位置 % 输入: a km, e 无量纲, 角度单位 deg, E rad % 输出: r_eci 3x1 km i deg2rad(i_deg); raan deg2rad(raan_deg); argp deg2rad(argp_deg); cosE cos(E); sinE sin(E); r_pf a * [cosE - e; sqrt(1 - e^2) * sinE; 0]; % 近焦点系 % 手写方向余弦矩阵Z-X-Z 顺序 R3 (th) [cos(th) sin(th) 0; -sin(th) cos(th) 0; 0 0 1]; R1 (th) [1 0 0; 0 cos(th) sin(th); 0 -sin(th) cos(th)]; r_eci R3(-raan) * R1(-i) * R3(-argp) * r_pf; end逻辑说明R3是绕 Z 轴旋转的矩阵R1是绕 X 轴旋转的矩阵。顺序不能调换因为轨道面的朝向是“先定升交点、再定倾角、最后定近地点”的物理次序。验证方法很简单当i0、RAAN0时结果应当等价于在赤道面内绕 Z 轴旋转-argp你可以用一组简单数自查。2.4 速度向量和轨道周期最近点校验要用的另一半后面校验最近点时要判断“相对位置与相对速度是否垂直”所以速度向量也要解析算出来不能用差分代替。近焦点系速度公式为v_pf sqrt(mu / p) * [-sin(E); sqrt(1 - e^2) * cos(E); 0]其中p a * (1 - e^2)mu 398600.4418 km^3/s^2。旋转矩阵与位置完全相同复用R3(-raan) * R1(-i) * R3(-argp)即可function v_eci orbital2eci_vel(a, e, i_deg, raan_deg, argp_deg, E) % 与 orbital2eci 相同的旋转顺序返回 3x1 km/s mu 398600.4418; p a * (1 - e^2); v_pf sqrt(mu / p) * [-sin(E); sqrt(1 - e^2) * cos(E); 0]; i deg2rad(i_deg); raan deg2rad(raan_deg); argp deg2rad(argp_deg); R3 (th) [cos(th) sin(th) 0; -sin(th) cos(th) 0; 0 0 1]; R1 (th) [1 0 0; 0 cos(th) sin(th); 0 -sin(th) cos(th)]; v_eci R3(-raan) * R1(-i) * R3(-argp) * v_pf; end轨道周期用于设置时间扫描范围T 2 * pi * sqrt(a^3 / mu)。两颗星周期不同相对距离曲线会出现拍频现象所以扫描窗口至少要覆盖其中较长轨道周期的一半否则可能漏掉最接近点。3. 用 Matlab 扫描相对距离曲线粗扫定位再精算到秒位置函数就绪后两颗卫星之间的距离变成单个标量函数d(t) |r1(t) - r2(t)|。看似简单但直接对 t 求导再求根容易被局部极值骗工程上我一般先均匀扫一遍再在最小值附近用fminbnd收敛。这正好也是 Matlab 优化工具箱的典型用法。3.1 把传播和距离封装成“输入时刻输出距离”的函数为了让fminbnd能反复调用需要把轨道传播和距离计算封装成“输入一个标量时间输出一个标量距离”的函数。两颗星的轨道根数放进结构体避免用全局变量污染工作区function d relative_distance(t, s1, s2) % t: 从历元起算的秒数 % s1, s2: 结构体字段为 a, e, i_deg, raan_deg, argp_deg, M0_deg mu 398600.4418; % WGS84 地球引力常数, km^3/s^2 n1 sqrt(mu / s1.a^3); n2 sqrt(mu / s2.a^3); % 平近点角随时间线性变化角度起点来自 M0 M1 deg2rad(s1.M0_deg) n1 * t; M2 deg2rad(s2.M0_deg) n2 * t; E1 kepler_E(mod(M1, 2*pi), s1.e, 1e-10); E2 kepler_E(mod(M2, 2*pi), s2.e, 1e-10); r1 orbital2eci(s1.a, s1.e, s1.i_deg, s1.raan_deg, s1.argp_deg, E1); r2 orbital2eci(s2.a, s2.e, s2.i_deg, s2.raan_deg, s2.argp_deg, E2); d norm(r1 - r2); end逻辑说明M0_deg是结构体里的平近点角单位度内部用deg2rad换成弧度再加n * t。mod把平近点角约束在 0 到 2π避免时间久后浮点误差把角度算飘。注意两颗星的历元必须一致如果 TLE 的 epoch 不同要先统一到同一个 UTC 时间点再换算各自的平近点角。3.2 粗扫描找到最小值区间再用 fminbnd 收敛单独调用fminbnd有可能收敛到局部极小值所以正确做法是先用固定步长算出完整距离曲线定位全局最小值落在哪个时间栅格里再在这个栅格附近用fminbnd做黄金分割搜索。示例代码% 两小时预报步长 5 秒 tspan 0:5:7200; dvals zeros(size(tspan)); for k 1:numel(tspan) dvals(k) relative_distance(tspan(k), s1, s2); end [dmin_grid, idx] min(dvals); t_coarse tspan(idx); % 在粗扫栅格附近精算 [t_close, d_close] fminbnd((tt) relative_distance(tt, s1, s2), ... t_coarse - 10, t_coarse 10, ... optimset(TolX, 1e-6)); fprintf(粗扫最小距离 %.3f km %.1f s\n, dmin_grid, t_coarse); fprintf(精算最近距离 %.6f km %.3f s\n, d_close, t_close);参数说明时间步长 5 秒是经验值。近地轨道相对接近速度通常 7 km/s 左右5 秒对应约 35 km 栅格足够把最近点定位到分钟级窗口如果两颗星相对速度更大或者轨道几何导致距离曲线出现尖锐凹谷需要把步长缩到 0.5 秒。fminbnd搜索区间取粗扫点前后 10 秒保证局部极小值落在区间内。optimset(TolX, 1e-6)控制时间收敛精度到微秒量级。不同步长对结果的影响可以用下面这张表快速估计时间步长栅格分辨率约适用场景60 s数百 km快速筛查一整天内的大尺度接近5 s约 35 km日常防碰撞接近筛查0.5 s约 3.5 km接近时刻的二次确认3.3 用相对速度垂直条件验证最近点最近距离处相对距离平方对时间的导数为零等价于相对位置向量和相对速度向量的点积为零。这个几何条件不依赖优化器适合独立验证fminbnd的搜索结果。校验代码如下function [r_rel, v_rel, dot_rv] relative_state(t, s1, s2) % 同时返回相对位置、相对速度和点积 % 位置、速度计算复用前面的 orbital2eci 和 orbital2eci_vel mu 398600.4418; n1 sqrt(mu / s1.a^3); n2 sqrt(mu / s2.a^3); M1 deg2rad(s1.M0_deg) n1 * t; M2 deg2rad(s2.M0_deg) n2 * t; E1 kepler_E(mod(M1, 2*pi), s1.e, 1e-10); E2 kepler_E(mod(M2, 2*pi), s2.e, 1e-10); r1 orbital2eci(s1.a, s1.e, s1.i_deg, s1.raan_deg, s1.argp_deg, E1); r2 orbital2eci(s2.a, s2.e, s2.i_deg, s2.raan_deg, s2.argp_deg, E2); v1 orbital2eci_vel(s1.a, s1.e, s1.i_deg, s1.raan_deg, s1.argp_deg, E1); v2 orbital2eci_vel(s2.a, s2.e, s2.i_deg, s2.raan_deg, s2.argp_deg, E2); r_rel r1 - r2; v_rel v1 - v2; dot_rv dot(r_rel, v_rel); end调用时只看dot_rv的数量级[~, ~, dot_rv] relative_state(t_close, s1, s2); fprintf(最近点处 dot(r_rel, v_rel) %.3e km^2/s\n, dot_rv);理论上最近点处应严格为零。由于fminbnd的容差限制实际输出可能在1e-4量级如果看到1e-1以上说明粗扫步长太大局部极小值没有真正被包住或者输入轨道根数的时间基准不一致。4. 时间系统、坐标转换与 J2 摄动别让结果差出几十公里二体模型在短时间窗口内已经够用但真实任务里还藏着三个额外的坑时间系统、坐标系混用、以及长时间预报下的摄动。下面逐个说清。4.1 用 datetime 管理历元时刻用“秒”统一轨道根数来自 TLE 或精密星历时历元通常是 UTC 日期字符串。在 Matlab 里用datetime保存历元计算时用seconds(t - t0)得到相对秒数。我一般不会用datenum做减法再乘 86400因为时区、闰秒和浮点误差会让时间差偏移 1 秒对近地轨道1 秒对应的沿轨位移约 7 km直接毁掉 CA 结果。t0 datetime(2026-01-01 00:00:00, InputFormat, yyyy-MM-dd HH:mm:ss); t_end t0 hours(2); t_sec seconds(t_end - t0); % 7200 秒 % 打印最近点时刻 ca_time_utc t0 seconds(t_close); disp(char(ca_time_utc, yyyy-MM-dd HH:mm:ss.SSS));这里的char转换就是常见的“datetime 转 string”需求。注意datetime默认的时区是本地时区处理 TLE 时应该显式指定TimeZone, UTC避免本地时区偏移污染时间差。如果你用datenum做测试至少要在脚本开头写清所有时间都是相对秒数不要混用绝对时刻和相对时刻。4.2 ECEF 和 ECI 不能混用加一次 GMST 旋转如果轨道根数来自站心观测或者某些广播星历输出的位置可能是 ECEF 而不是 ECI。直接用 ECEF 坐标算相对距离相当于忽略地球自转几分钟内就能造成几十千米误差。对于经典轨道六根数位置默认在惯性系 ECI 里不需要这步但当你拿到的是 ECEF 位置时必须转回 ECI。Matlab 的 Aerospace Toolbox 提供ecef2eci如果没装这个工具箱可以用简化 GMST 旋转自己写一个function r_eci ecef2eci_gmst(r_ecef, utc_time) % r_ecef: 3x1 km % utc_time: datetime, UTC 时区 jd juliandate(utc_time); T (jd - 2451545.0) / 36525; % 儒略世纪数 gmst mod(280.46061837 360.98564736629 * (jd - 2451545.0) ... 0.000387933 * T^2, 360); gmst_rad deg2rad(gmst); R [cos(gmst_rad) sin(gmst_rad) 0; -sin(gmst_rad) cos(gmst_rad) 0; 0 0 1]; r_eci R * r_ecef; end逻辑说明GMST 是格林尼治恒星时表示地球相对惯性系转过的角度。这个简化版本没有包含极移和章动对近地接近判断的精度在秒级以内。如果追求严格可以用ecef2eci_icrf等工具箱函数但日常筛查用 GMST 就够了。代码里必须把utc_time的TimeZone设为UTC否则juliandate会按本地时区解释差出 8 小时就是 120° 的旋转误差。4.3 长时间预报把 J2 造成的轨道面漂移加进去二体模型在几小时内误差有多大以 500 km 高度 LEO 卫星为例J2 摄动主要让 RAAN 和近地点幅角线性漂移RAAN 漂移率对非极轨卫星可达每天几度到十几度。两小时漂移不到 1°对最近距离的影响在高轨大偏心率场景下可能从几百米放大到几公里。常见做法是把 J2 的一阶长期项加到平均根数上。对近地卫星RAAN 和 argp 的漂移率近似为RAAN_dot -1.5 * J2 * (Re / p)^2 * n * cos(i) argp_dot 0.75 * J2 * (Re / p)^2 * n * (5 * cos(i)^2 - 1)写成 Matlab 函数function [raan, argp] j2_drift(t, a, e, i_deg, raan0, argp0) % t 秒, a km, e 无量纲, 角度单位 deg % 返回 t 时刻的 RAAN 和 argp, 单位 deg mu 398600.4418; Re 6378.137; J2 1.08262668e-3; p a * (1 - e^2); n sqrt(mu / a^3); % rad/s i deg2rad(i_deg); raan_dot -1.5 * J2 * (Re / p)^2 * n * cos(i); argp_dot 0.75 * J2 * (Re / p)^2 * n * (5 * cos(i)^2 - 1); raan raan0 raan_dot * t * 180 / pi; argp argp0 argp_dot * t * 180 / pi; end调用时把j2_drift返回的raan、argp传给orbital2eci其余流程不变。这里有一个很隐蔽的坑TLE 给的是平均根数本身已经包含了 J2 的一阶长期项再往上加等于重复计算。这套j2_drift只适用于你自己用密切根数做初值、又想快速评估摄动影响的场景。对接近预警二体加粗扫描已经能筛掉 99% 的无威胁事件真到了要报碰撞概率时再换成 SGP4 或数值积分。5. 把最近距离曲线画出来接近窗口自动标红5.1 距离曲线与时标算完 CA 只输出两个数还不够。距离曲线能立刻暴露三类问题粗扫步长是否大到漏谷、是否存在多个几乎持平的局部极小值、结束时刻相对距离是否重新拉高。下面画图脚本会标出距离低于阈值的所有窗口并把精算出的最近点用竖线固定下来。thresh 10; % 报警阈值, km tplot 0:1:7200; % 1 秒重采样画图即可不必用 0.1 秒 dplot zeros(size(tplot)); for k 1:numel(tplot) dplot(k) relative_distance(tplot(k), s1, s2); end figure(Color, w); plot(tplot / 60, dplot, LineWidth, 1.2); hold on; below dplot thresh; if any(below) area(tplot(below) / 60, dplot(below), ... FaceColor, r, FaceAlpha, 0.3); end yline(thresh, --k); xline(t_close / 60, -, sprintf(CA %.2f km, d_close)); xlabel(Time since epoch (min)); ylabel(Relative distance (km)); legend(distance, below threshold, threshold, closest approach, ... Location, northwest); grid on;逻辑说明area只填充离散点中低于阈值的部分如果 1 秒采样仍漏掉极窄的接近窗口先用find(diff(below) ~ 0)找出边界索引再在边界附近加密到 0.1 秒。xline的标签直接给出接近时刻和最近距离配合图例能一眼看出最小值的形态。5.2 输出接近报告表格把结果写成 CSV方便后续批量处理或多圈筛选T table(t_close, d_close, dot_rv, ... VariableNames, {Time_s, Distance_km, Dot_rv_km2_s}); writetable(T, close_approach.csv);这里dot_rv是上一节校验函数输出的最近点处点积。CSV 里建议把相对秒数换成 UTC 时间字符串用t0 seconds(t_close)生成datetime再char转 string 写入表格。这个脚本不依赖 Aerospace Toolbox只要有基础 Matlab 环境就能跑把s1、s2两个结构体换成从 TLE 解析出的轨道根数就可以直接扩展成批量接近筛查工具每次多圈扫描只需要改tspan和阈值。本文还有配套的精品资源点击获取 
长河君 · 2026/9/12 20:01:48 · 阅读 1.2k Matlab实现卫星最近距离计算：从轨道根数到CA实践
简介一个MATLAB工具包面向航天工程与轨道分析领域用于求解圆形或椭圆形地球轨道上两颗卫星之间的最近接近距离。脚本基于Brent一维最小化算法定位几何最近条件轨道传播采用Kozai方法并将地球半径增加2%以近似大气层阻力影响适合航天工程师、科研人员及空间技术学习者进行碰撞风险评估与轨道安全分析。压缩包共14个文件、230KB以12个m源码文件为主另有1个pdf说明文档和1个in输入数据示例各模块划分清晰如主控脚本、距离计算函数、最小化模块、Kozai传播函数、儒略日转换及结果打印等便于按需修改与复用调试。目前已有49人学习下载借助该工具可快速掌握卫星近距离接近的数值求解思路并依据输入参数模拟不同轨道场景为空间碎片规避、防碰撞预警与轨道设计提供实用参考。1. 两个地球卫星的最近距离先解决“时间”再谈“距离”手里拿着两颗卫星的两行根数想知道它们在接下来几个小时内会不会靠得太近这是空间态势感知里最常见的提问。两颗地球卫星之间的最近距离在英文里叫 Closest Approach防碰撞报告里经常简写成 CA。计算它不是把两个瞬时位置丢进欧氏距离公式就结束真正的难点在时间——每一颗卫星都在自己的开普勒轨道上漂移最近距离不是一个静态几何量而是一条随时间变化的相对距离曲线上的最小值。我会用 Matlab 从轨道根数重建 ECI 位置先粗扫再用优化工具箱精化把“最近距离、发生时刻、相对速度方向”一起算出来。这套流程适合写过轨道仿真但没有系统处理过角度换算和时间系统的人也适合想在任务规划里加一个独立校验脚本的工程师。2. 从轨道根数算卫星位置坐标系和开普勒方程必须先写对在把两颗卫星放到同一张图之前得先把各自的位置算对。这一步不涉及高级数学但单位、坐标系和旋转顺序三个地方错一个后面的距离计算就全白做。下面这套代码基于二体模型适合几分钟到几小时的接近筛查如果需要高精度交会评估第 4 章再补充时间系统和 J2 摄动。2.1 轨道六根数输入格式与单位是第一个坑经典轨道六根数是最常见的输入分别决定轨道的大小、形状、倾斜和当前相位。对“两颗卫星最近距离”这个问题两颗卫星必须在同一历元下定义时间轴从同一时刻起算否则算出来的 CA 时刻毫无意义。符号含义常用单位说明a半长轴km圆轨道下就是轨道半径LEO 约 6800 kme偏心率无量纲0 到 10.001 以下按小偏心率处理i轨道倾角deg相对赤道面0° 赤道轨道90° 极轨RAAN升交点赤经deg春分点到升交点的夹角argp近地点幅角deg升交点到近地点的夹角M0历元平近点角deg时间变量的起点替代真近点角作为输入单位坑集中在角度上。Matlab 的三角函数默认弧度而 TLE、星历表里给的都是角度。第一步就是把[i, RAAN, argp, M0]全部用deg2rad换算。偏心率也要检查必须满足0 e 1对偏心率接近 1 的细长轨道牛顿迭代初始值要另做处理对低轨卫星e 一般在 0.001 量级直接用平近点角当初始值就够了。2.2 用牛顿迭代解开普勒方程避免 atan2 奇异点平近点角 M 和偏近点角 E 的关系是开普勒方程M E - e * sin(E)给定 M 求 E 是超越方程工程上最常用牛顿法。迭代式是E_new E - (E - esin(E) - M) / (1 - ecos(E))分母在 e 接近 1 且 E 接近 2π 时会退化但近地轨道不会遇到。下面这个函数可以直接保存成kepler_E.mfunction E kepler_E(M, e, tol) % M 平近点角, rad % e 偏心率, 0 e 1 % tol 迭代容差, 默认 1e-10 if nargin 3, tol 1e-10; end E M; % 小偏心率轨道从平近点角出发收敛足够快 for k 1:100 dE (E - e*sin(E) - M) / (1 - e*cos(E)); E E - dE; if abs(dE) tol break; end end if k 100 warning(kepler_E:notConverged, 开普勒方程未在100步内收敛); end end逻辑说明迭代的核心是每一步计算修正量dE当修正量小于tol时停止。平近点角在传播时会持续增加调用前先做mod(M, 2*pi)避免大时间差下角度累积误差。tol1e-10对应偏近点角精度约 1e-10 rad换算成位置误差在毫米量级足够 CA 计算使用。2.3 从轨道面坐标转到 ECI旋转顺序不能记反偏近点角 E 算出后近焦点坐标系中的位置矢量为r_pf [a * (cos(E) - e); a * sqrt(1 - e^2) * sin(E); 0]要把这个矢量转到 ECI通常指 J2000 或 MOD需要依次绕 Z 轴旋转 -RAAN、绕 X 轴旋转 -i、绕 Z 轴旋转 -argp。这里的旋转方向约定来自 Curtis 的《Orbital Mechanics for Engineering Students》和 Matlab 自带的rotz方向不一定一致所以最好手写旋转矩阵function r_eci orbital2eci(a, e, i_deg, raan_deg, argp_deg, E) % 由半长轴、偏心率、倾角、升交点赤经、近地点幅角、偏近点角计算 ECI 位置 % 输入: a km, e 无量纲, 角度单位 deg, E rad % 输出: r_eci 3x1 km i deg2rad(i_deg); raan deg2rad(raan_deg); argp deg2rad(argp_deg); cosE cos(E); sinE sin(E); r_pf a * [cosE - e; sqrt(1 - e^2) * sinE; 0]; % 近焦点系 % 手写方向余弦矩阵Z-X-Z 顺序 R3 (th) [cos(th) sin(th) 0; -sin(th) cos(th) 0; 0 0 1]; R1 (th) [1 0 0; 0 cos(th) sin(th); 0 -sin(th) cos(th)]; r_eci R3(-raan) * R1(-i) * R3(-argp) * r_pf; end逻辑说明R3是绕 Z 轴旋转的矩阵R1是绕 X 轴旋转的矩阵。顺序不能调换因为轨道面的朝向是“先定升交点、再定倾角、最后定近地点”的物理次序。验证方法很简单当i0、RAAN0时结果应当等价于在赤道面内绕 Z 轴旋转-argp你可以用一组简单数自查。2.4 速度向量和轨道周期最近点校验要用的另一半后面校验最近点时要判断“相对位置与相对速度是否垂直”所以速度向量也要解析算出来不能用差分代替。近焦点系速度公式为v_pf sqrt(mu / p) * [-sin(E); sqrt(1 - e^2) * cos(E); 0]其中p a * (1 - e^2)mu 398600.4418 km^3/s^2。旋转矩阵与位置完全相同复用R3(-raan) * R1(-i) * R3(-argp)即可function v_eci orbital2eci_vel(a, e, i_deg, raan_deg, argp_deg, E) % 与 orbital2eci 相同的旋转顺序返回 3x1 km/s mu 398600.4418; p a * (1 - e^2); v_pf sqrt(mu / p) * [-sin(E); sqrt(1 - e^2) * cos(E); 0]; i deg2rad(i_deg); raan deg2rad(raan_deg); argp deg2rad(argp_deg); R3 (th) [cos(th) sin(th) 0; -sin(th) cos(th) 0; 0 0 1]; R1 (th) [1 0 0; 0 cos(th) sin(th); 0 -sin(th) cos(th)]; v_eci R3(-raan) * R1(-i) * R3(-argp) * v_pf; end轨道周期用于设置时间扫描范围T 2 * pi * sqrt(a^3 / mu)。两颗星周期不同相对距离曲线会出现拍频现象所以扫描窗口至少要覆盖其中较长轨道周期的一半否则可能漏掉最接近点。3. 用 Matlab 扫描相对距离曲线粗扫定位再精算到秒位置函数就绪后两颗卫星之间的距离变成单个标量函数d(t) |r1(t) - r2(t)|。看似简单但直接对 t 求导再求根容易被局部极值骗工程上我一般先均匀扫一遍再在最小值附近用fminbnd收敛。这正好也是 Matlab 优化工具箱的典型用法。3.1 把传播和距离封装成“输入时刻输出距离”的函数为了让fminbnd能反复调用需要把轨道传播和距离计算封装成“输入一个标量时间输出一个标量距离”的函数。两颗星的轨道根数放进结构体避免用全局变量污染工作区function d relative_distance(t, s1, s2) % t: 从历元起算的秒数 % s1, s2: 结构体字段为 a, e, i_deg, raan_deg, argp_deg, M0_deg mu 398600.4418; % WGS84 地球引力常数, km^3/s^2 n1 sqrt(mu / s1.a^3); n2 sqrt(mu / s2.a^3); % 平近点角随时间线性变化角度起点来自 M0 M1 deg2rad(s1.M0_deg) n1 * t; M2 deg2rad(s2.M0_deg) n2 * t; E1 kepler_E(mod(M1, 2*pi), s1.e, 1e-10); E2 kepler_E(mod(M2, 2*pi), s2.e, 1e-10); r1 orbital2eci(s1.a, s1.e, s1.i_deg, s1.raan_deg, s1.argp_deg, E1); r2 orbital2eci(s2.a, s2.e, s2.i_deg, s2.raan_deg, s2.argp_deg, E2); d norm(r1 - r2); end逻辑说明M0_deg是结构体里的平近点角单位度内部用deg2rad换成弧度再加n * t。mod把平近点角约束在 0 到 2π避免时间久后浮点误差把角度算飘。注意两颗星的历元必须一致如果 TLE 的 epoch 不同要先统一到同一个 UTC 时间点再换算各自的平近点角。3.2 粗扫描找到最小值区间再用 fminbnd 收敛单独调用fminbnd有可能收敛到局部极小值所以正确做法是先用固定步长算出完整距离曲线定位全局最小值落在哪个时间栅格里再在这个栅格附近用fminbnd做黄金分割搜索。示例代码% 两小时预报步长 5 秒 tspan 0:5:7200; dvals zeros(size(tspan)); for k 1:numel(tspan) dvals(k) relative_distance(tspan(k), s1, s2); end [dmin_grid, idx] min(dvals); t_coarse tspan(idx); % 在粗扫栅格附近精算 [t_close, d_close] fminbnd((tt) relative_distance(tt, s1, s2), ... t_coarse - 10, t_coarse 10, ... optimset(TolX, 1e-6)); fprintf(粗扫最小距离 %.3f km %.1f s\n, dmin_grid, t_coarse); fprintf(精算最近距离 %.6f km %.3f s\n, d_close, t_close);参数说明时间步长 5 秒是经验值。近地轨道相对接近速度通常 7 km/s 左右5 秒对应约 35 km 栅格足够把最近点定位到分钟级窗口如果两颗星相对速度更大或者轨道几何导致距离曲线出现尖锐凹谷需要把步长缩到 0.5 秒。fminbnd搜索区间取粗扫点前后 10 秒保证局部极小值落在区间内。optimset(TolX, 1e-6)控制时间收敛精度到微秒量级。不同步长对结果的影响可以用下面这张表快速估计时间步长栅格分辨率约适用场景60 s数百 km快速筛查一整天内的大尺度接近5 s约 35 km日常防碰撞接近筛查0.5 s约 3.5 km接近时刻的二次确认3.3 用相对速度垂直条件验证最近点最近距离处相对距离平方对时间的导数为零等价于相对位置向量和相对速度向量的点积为零。这个几何条件不依赖优化器适合独立验证fminbnd的搜索结果。校验代码如下function [r_rel, v_rel, dot_rv] relative_state(t, s1, s2) % 同时返回相对位置、相对速度和点积 % 位置、速度计算复用前面的 orbital2eci 和 orbital2eci_vel mu 398600.4418; n1 sqrt(mu / s1.a^3); n2 sqrt(mu / s2.a^3); M1 deg2rad(s1.M0_deg) n1 * t; M2 deg2rad(s2.M0_deg) n2 * t; E1 kepler_E(mod(M1, 2*pi), s1.e, 1e-10); E2 kepler_E(mod(M2, 2*pi), s2.e, 1e-10); r1 orbital2eci(s1.a, s1.e, s1.i_deg, s1.raan_deg, s1.argp_deg, E1); r2 orbital2eci(s2.a, s2.e, s2.i_deg, s2.raan_deg, s2.argp_deg, E2); v1 orbital2eci_vel(s1.a, s1.e, s1.i_deg, s1.raan_deg, s1.argp_deg, E1); v2 orbital2eci_vel(s2.a, s2.e, s2.i_deg, s2.raan_deg, s2.argp_deg, E2); r_rel r1 - r2; v_rel v1 - v2; dot_rv dot(r_rel, v_rel); end调用时只看dot_rv的数量级[~, ~, dot_rv] relative_state(t_close, s1, s2); fprintf(最近点处 dot(r_rel, v_rel) %.3e km^2/s\n, dot_rv);理论上最近点处应严格为零。由于fminbnd的容差限制实际输出可能在1e-4量级如果看到1e-1以上说明粗扫步长太大局部极小值没有真正被包住或者输入轨道根数的时间基准不一致。4. 时间系统、坐标转换与 J2 摄动别让结果差出几十公里二体模型在短时间窗口内已经够用但真实任务里还藏着三个额外的坑时间系统、坐标系混用、以及长时间预报下的摄动。下面逐个说清。4.1 用 datetime 管理历元时刻用“秒”统一轨道根数来自 TLE 或精密星历时历元通常是 UTC 日期字符串。在 Matlab 里用datetime保存历元计算时用seconds(t - t0)得到相对秒数。我一般不会用datenum做减法再乘 86400因为时区、闰秒和浮点误差会让时间差偏移 1 秒对近地轨道1 秒对应的沿轨位移约 7 km直接毁掉 CA 结果。t0 datetime(2026-01-01 00:00:00, InputFormat, yyyy-MM-dd HH:mm:ss); t_end t0 hours(2); t_sec seconds(t_end - t0); % 7200 秒 % 打印最近点时刻 ca_time_utc t0 seconds(t_close); disp(char(ca_time_utc, yyyy-MM-dd HH:mm:ss.SSS));这里的char转换就是常见的“datetime 转 string”需求。注意datetime默认的时区是本地时区处理 TLE 时应该显式指定TimeZone, UTC避免本地时区偏移污染时间差。如果你用datenum做测试至少要在脚本开头写清所有时间都是相对秒数不要混用绝对时刻和相对时刻。4.2 ECEF 和 ECI 不能混用加一次 GMST 旋转如果轨道根数来自站心观测或者某些广播星历输出的位置可能是 ECEF 而不是 ECI。直接用 ECEF 坐标算相对距离相当于忽略地球自转几分钟内就能造成几十千米误差。对于经典轨道六根数位置默认在惯性系 ECI 里不需要这步但当你拿到的是 ECEF 位置时必须转回 ECI。Matlab 的 Aerospace Toolbox 提供ecef2eci如果没装这个工具箱可以用简化 GMST 旋转自己写一个function r_eci ecef2eci_gmst(r_ecef, utc_time) % r_ecef: 3x1 km % utc_time: datetime, UTC 时区 jd juliandate(utc_time); T (jd - 2451545.0) / 36525; % 儒略世纪数 gmst mod(280.46061837 360.98564736629 * (jd - 2451545.0) ... 0.000387933 * T^2, 360); gmst_rad deg2rad(gmst); R [cos(gmst_rad) sin(gmst_rad) 0; -sin(gmst_rad) cos(gmst_rad) 0; 0 0 1]; r_eci R * r_ecef; end逻辑说明GMST 是格林尼治恒星时表示地球相对惯性系转过的角度。这个简化版本没有包含极移和章动对近地接近判断的精度在秒级以内。如果追求严格可以用ecef2eci_icrf等工具箱函数但日常筛查用 GMST 就够了。代码里必须把utc_time的TimeZone设为UTC否则juliandate会按本地时区解释差出 8 小时就是 120° 的旋转误差。4.3 长时间预报把 J2 造成的轨道面漂移加进去二体模型在几小时内误差有多大以 500 km 高度 LEO 卫星为例J2 摄动主要让 RAAN 和近地点幅角线性漂移RAAN 漂移率对非极轨卫星可达每天几度到十几度。两小时漂移不到 1°对最近距离的影响在高轨大偏心率场景下可能从几百米放大到几公里。常见做法是把 J2 的一阶长期项加到平均根数上。对近地卫星RAAN 和 argp 的漂移率近似为RAAN_dot -1.5 * J2 * (Re / p)^2 * n * cos(i) argp_dot 0.75 * J2 * (Re / p)^2 * n * (5 * cos(i)^2 - 1)写成 Matlab 函数function [raan, argp] j2_drift(t, a, e, i_deg, raan0, argp0) % t 秒, a km, e 无量纲, 角度单位 deg % 返回 t 时刻的 RAAN 和 argp, 单位 deg mu 398600.4418; Re 6378.137; J2 1.08262668e-3; p a * (1 - e^2); n sqrt(mu / a^3); % rad/s i deg2rad(i_deg); raan_dot -1.5 * J2 * (Re / p)^2 * n * cos(i); argp_dot 0.75 * J2 * (Re / p)^2 * n * (5 * cos(i)^2 - 1); raan raan0 raan_dot * t * 180 / pi; argp argp0 argp_dot * t * 180 / pi; end调用时把j2_drift返回的raan、argp传给orbital2eci其余流程不变。这里有一个很隐蔽的坑TLE 给的是平均根数本身已经包含了 J2 的一阶长期项再往上加等于重复计算。这套j2_drift只适用于你自己用密切根数做初值、又想快速评估摄动影响的场景。对接近预警二体加粗扫描已经能筛掉 99% 的无威胁事件真到了要报碰撞概率时再换成 SGP4 或数值积分。5. 把最近距离曲线画出来接近窗口自动标红5.1 距离曲线与时标算完 CA 只输出两个数还不够。距离曲线能立刻暴露三类问题粗扫步长是否大到漏谷、是否存在多个几乎持平的局部极小值、结束时刻相对距离是否重新拉高。下面画图脚本会标出距离低于阈值的所有窗口并把精算出的最近点用竖线固定下来。thresh 10; % 报警阈值, km tplot 0:1:7200; % 1 秒重采样画图即可不必用 0.1 秒 dplot zeros(size(tplot)); for k 1:numel(tplot) dplot(k) relative_distance(tplot(k), s1, s2); end figure(Color, w); plot(tplot / 60, dplot, LineWidth, 1.2); hold on; below dplot thresh; if any(below) area(tplot(below) / 60, dplot(below), ... FaceColor, r, FaceAlpha, 0.3); end yline(thresh, --k); xline(t_close / 60, -, sprintf(CA %.2f km, d_close)); xlabel(Time since epoch (min)); ylabel(Relative distance (km)); legend(distance, below threshold, threshold, closest approach, ... Location, northwest); grid on;逻辑说明area只填充离散点中低于阈值的部分如果 1 秒采样仍漏掉极窄的接近窗口先用find(diff(below) ~ 0)找出边界索引再在边界附近加密到 0.1 秒。xline的标签直接给出接近时刻和最近距离配合图例能一眼看出最小值的形态。5.2 输出接近报告表格把结果写成 CSV方便后续批量处理或多圈筛选T table(t_close, d_close, dot_rv, ... VariableNames, {Time_s, Distance_km, Dot_rv_km2_s}); writetable(T, close_approach.csv);这里dot_rv是上一节校验函数输出的最近点处点积。CSV 里建议把相对秒数换成 UTC 时间字符串用t0 seconds(t_close)生成datetime再char转 string 写入表格。这个脚本不依赖 Aerospace Toolbox只要有基础 Matlab 环境就能跑把s1、s2两个结构体换成从 TLE 解析出的轨道根数就可以直接扩展成批量接近筛查工具每次多圈扫描只需要改tspan和阈值。本文还有配套的精品资源点击获取