用Python和SGP4计算低轨卫星可见性与通信延迟全流程-北京尧图网络

简介一份基于 Python 的低轨道卫星通信模拟项目面向通信、计算机、自动化等相关专业的学生与开发者适合作为课程设计、毕业设计或初期项目演示项目围绕低轨卫星轨道运动模拟、星地可见性判断与通信延迟计算展开同时模拟数据在卫星和地面站间的传输过程。压缩包共 28 个文件含 14 个 Python 脚本、11 个 pyc 编译文件、2 个文本说明和 1 个 CSV 仿真数据文件整体大小 708KB脚本按功能拆分覆盖轨道数值积分、姿态解算、磁强计模拟、通信链路仿真与结果可视化等模块并附带文本说明与仿真记录便于对照学习和二次开发。目前已有 369 人浏览学习。代码结构清晰、可直接运行既能帮助初学者快速建立卫星通信链路模拟的完整认知也可供进阶者修改和扩展实现自定义轨道与通信场景模块化设计还方便继续接入其他传感器模型或通信协议整体具有较好的学习借鉴价值。1. 低轨道卫星通信里可见性和延迟才是先要算清的两笔账做低轨道LEO卫星通信仿真很多人第一反应是去画星座图、算链路预算结果发现连“这颗卫星现在能不能跟地面站通信”都没法回答。实际上LEO 卫星的高速运动导致星地链路每隔几分钟就要切换一次卫星从地平线升起到落下地面站能用的时间窗口往往只有 515 分钟。这个窗口内仰角决定了链路质量星地距离决定了传播延迟。低轨道的高度范围通常在 4002000 km单程传播延迟大约 1.36.7 ms看似不大但考虑到 LEO 星座组网需要多跳中继端到端延迟会叠加直接影响 TCP 拥塞控制行为和实时业务的调度策略。标题里的“可见性”和“通信延迟”本质是两个独立但耦合的计算问题可见性需要先算卫星在惯性系中的位置再转换到地面站所在地心地固系中求仰角和方位角通信延迟则在可见的前提下用星地距离除以光速得到传播时延再叠加大气、处理和排队时延。Python 生态里有一套成熟的方案核心是 SGP4 轨道传播器和 TLE 轨道根数本文就用这两个工具把从轨道计算到可见性判断、再到延迟估算的完整链路跑通。2. 用 Python 和 SGP4 把 TLE 轨道根数变成卫星坐标2.1 为什么选 SGP4 而不自己去解轨道方程低轨道卫星的轨道计算有两类主流方法一类是直接数值积分摄动方程精度高但需要高精度的初始状态向量和力模型另一类是 SGP4Simplified General Perturbations Model 4它吃的是北美防空司令部发布的 TLE 双行轨道根数把地球非球形引力、大气阻力、日月引力等摄动以经验公式内置在模型中计算速度快精度对绝大多数通信仿真场景足够。TLE 数据在 CelesTrak 等平台公开获取每天更新是工程上做 LEO 卫星可见性预报的事实标准。SGP4 在 Python 里的实现有两个常用选择sgp4库提供底层接口skyfield在它之上封装了时间系统和坐标系转换更适合做地面站相关的计算。两者的关系是sgp4算出卫星在 TEMETrue Equator, Mean Equinox坐标系下的位置和速度skyfield负责把 TEME 坐标转到地心地固坐标系再和地面站的经纬度做几何解算。下面先看底层sgp4的用法。2.2 安装环境和解析 TLE 的最小代码先装依赖建议在虚拟环境里操作pip install sgp4 skyfield numpy然后拿一颗典型 LEO 卫星的 TLE 来解析。这里以国际空间站ISS的 TLE 为例格式如下ISS (ZARYA) 1 25544U 98067A 24050.55555556 .00016717 00000-0 10270-3 0 9999 2 25544 51.6400 247.4753 0007350 75.6057 60.1642 15.49418968456788解析并传播到指定时刻的代码from sgp4.api import Satrec from sgp4.api import jday from datetime import datetime, timezone # 两行TLE字符串来自CelesTrak或Space-Track line1 1 25544U 98067A 24050.55555556 .00016717 00000-0 10270-3 0 9999 line2 2 25544 51.6400 247.4753 0007350 75.6057 60.1642 15.49418968456788 sat Satrec.twoline2rv(line1, line2) # 目标时刻UTC 2024年2月19日 12:00:00 t datetime(2024, 2, 19, 12, 0, 0, tzinfotimezone.utc) jd, fr jday(t.year, t.month, t.day, t.hour, t.minute, t.second) # 传播得到TEME坐标系下的位置(km)和速度(km/s) e, r, v sat.sgp4(jd, fr) if e ! 0: print(fSGP4 传播错误错误码: {e}) else: print(f位置 (TEME, km): x{r[0]:.2f}, y{r[1]:.2f}, z{r[2]:.2f}) print(f速度 (TEME, km/s): vx{v[0]:.4f}, vy{v[1]:.4f}, vz{v[2]:.4f})这段代码里有几个关键点。Satrec.twoline2rv把两行 TLE 解析成内部轨道状态jday函数把 UTC 时间拆成年月日时分秒并转换成儒略日Julian Date注意返回值有两个jd是儒略日的整数部分fr是小数部分两者相加才是一个完整的儒略日时刻。sgp4(jd, fr)返回三个值e是错误码0 表示正常r是 TEME 坐标系下的位置向量单位是 kmv是速度向量单位是 km/s。这里的坑在于坐标系。sgp4库直接输出的是 TEME 坐标系它相对惯性空间有缓慢的极移和岁差章动修正不能直接拿来和地面站的经纬度做几何计算。要么手动做 TEME 到 ECEF 的转换需要读取地球定向参数要么直接换用skyfield库它内部封装了这套转换逻辑。工程上更常见的做法是直接用skyfield。2.3 用 Skyfield 直接拿到地心地固坐标skyfield的load.tle()可以从字符串或文件加载 TLEEarthSatellite对象可以直接通过at(time)方法传播from skyfield.api import load, wgs84 ts load.timescale() satellite load.tle(https://celestrak.org/NORAD/elements/gp.php?CATNR25544) # 用一个时间序列从当前时刻开始每隔10秒采样一次共50个点 t0 ts.now() t_array ts.utc(t0.utc_datetime() [__import__(datetime).timedelta(secondsi*10) for i in range(50)]) # 传播得到GCRS坐标系下的位置然后转成地面经纬度 geocentric satellite.at(t_array) # 地心地固系下的经纬度和高度 subpoint wgs84.subpoint(geocentric) lat subpoint.latitude.degrees lon subpoint.longitude.degrees alt subpoint.elevation.km print(f第一个时刻: 纬度 {lat[0]:.3f}°, 经度 {lon[0]:.3f}°, 高度 {alt[0]:.3f} km)wgs84.subpoint()返回的是卫星星下点subpoint的纬度和经度这对应卫星在地表正下方的投影点。注意区分两个概念星下点经纬度描述的是卫星自己在哪里而我们真正要算的是地面站看卫星的仰角和方位角这两者之间还差一个地面站位置参与的空间几何关系。3. 地面站可见性判定仰角计算和最小仰角阈值3.1 可见性的几何定义地面站收发天线指向卫星的视线方向与地面站水平面的夹角叫仰角elevation范围从 -90° 到 90°。卫星在地平线以下时仰角为负此时信号被地球本身遮挡无法通信卫星刚冒出地平线时仰角为 0°随着卫星升高仰角增大到卫星过顶时仰角接近 90°。实际系统中天线不可能贴着地平线工作——城市环境中有建筑物遮挡低仰角时大气衰减和地面噪声温度都会恶化链路质量所以工程上会设一个最小工作仰角典型值是 5° 到 10°。地面站与卫星之间满足“卫星仰角 ≥ 最小仰角”的时间区间就是一次可见窗口access window。可见性计算本质上是一个三维几何问题。已知卫星在 ECEF 坐标系下的位置向量r_sat和地面站的位置向量r_station先求出地面站指向卫星的向量d r_sat - r_station再把这个向量和地面站水平面的法向量做点积取反正弦就得到仰角。难点在于地面站位置向量本身依赖经纬度和海拔需要从大地坐标转成 ECEF 直角坐标。3.2 用 Skyfield 计算仰角和方位角skyfield对这整套几何做了一层高层封装altaz()方法直接返回给定地面站视角下的仰角、方位角和距离from skyfield.api import load, wgs84 from datetime import timedelta # 加载TLE satellite load.tle(https://celestrak.org/NORAD/elements/gp.php?CATNR25544) # 地面站北京某坐标海拔50米 station wgs84.latlon(39.9042, 116.4074, elevation_m50) ts load.timescale() # 以当前时刻为起点往后扫2小时每5秒一个采样点 now ts.now() times ts.utc([now.utc_datetime() timedelta(secondsi*5) for i in range(1440)]) # 计算整个时间序列的仰角/方位角/距离 difference satellite - station topocentric difference.at(times) alt, az, distance topocentric.altaz() # 打印第一个可见窗口 min_elev 10 # 最小工作仰角10度 visible alt.degrees min_elev for i in range(1, len(visible)): if visible[i] and not visible[i-1]: print(f进入可见窗口: {times[i].utc_strftime(%H:%M:%S)}, 仰角 {alt.degrees[i]:.2f}°) if not visible[i] and visible[i-1]: print(f离开可见窗口: {times[i].utc_strftime(%H:%M:%S)}, 仰角 {alt.degrees[i]:.2f}°) # 当探测到可见时输出距离 if any(visible): idx [i for i, v in enumerate(visible) if v] mid idx[len(idx)//2] print(f窗口内中间时刻: 仰角 {alt.degrees[mid]:.2f}°, 方位角 {az.degrees[mid]:.2f}°, 星地距离 {distance.km[mid]:.2f} km)这里difference satellite - station构建了一个“从地面站看卫星”的观测几何对象at(times)沿时间序列传播altaz()一次性给出仰角、方位角和距离。alt.degrees是 NumPy 数组直接用布尔比较生成可见性掩码再检测掩码的上升沿和下降沿就能定位窗口边界。需要提一个精度问题skyfield的altaz()默认不做大气折射修正因为大气折射对仰角的影响在低仰角时才显著仰角 0° 时折射约 0.5°而我们的最小仰角通常设在 5° 以上工程上直接忽略折射是合理的。如果要更精确altaz()有一个refractionTrue参数但注意它依赖气温气压没有气象数据时反而引入额外误差。3.3 手动计算 ECEF 坐标理解底层才能调参用库固然省事但只靠altaz()你没法理解距离和延迟的关系。这里给出手动从大地坐标转 ECEF 的实现方便你在没有skyfield的环境里也能算import math import numpy as np def geodetic_to_ecef(lat_deg, lon_deg, alt_m): a 6378.137 # 地球长半轴 km f 1 / 298.257223563 # WGS-84扁率 e2 f * (2 - f) lat math.radians(lat_deg) lon math.radians(lon_deg) N a / math.sqrt(1 - e2 * math.sin(lat)**2) x (N alt_m/1000) * math.cos(lat) * math.cos(lon) y (N alt_m/1000) * math.cos(lat) * math.sin(lon) z (N * (1 - e2) alt_m/1000) * math.sin(lat) return np.array([x, y, z]) def elevation_angle(r_sat_ecef, r_station_ecef): d r_sat_ecef - r_station_ecef dist np.linalg.norm(d) # 地面站位置向量的单位向量近似作为当地垂线方向 up r_station_ecef / np.linalg.norm(r_station_ecef) # 仰角 90° - 视线与垂线的夹角 cos_zenith np.dot(d, up) / dist elev 90.0 - math.degrees(math.acos(np.clip(cos_zenith, -1, 1))) return elev, dist station_ecef geodetic_to_ecef(39.9042, 116.4074, 50) # 假设从skyfield或sgp4拿到了ECEF坐标 sat_ecef np.array([-3478.2, 4752.3, 3456.7]) # 示例值 elev, dist elevation_angle(sat_ecef, station_ecef) print(f仰角: {elev:.2f}°, 距离: {dist:.2f} km)这段手动代码用了 WGS-84 椭球模型geodetic_to_ecef先把经纬度和海拔换成地心直角坐标。仰角计算把地面站的 ECEF 位置方向的单位向量近似当作当地垂线方向这在工程精度范围内是成立的因为地面站海拔远小于地球半径。注意sat_ecef需要是真正的 ECEF 坐标TEME 坐标直接喂进来会算错这是最容易踩的坑。4. 通信延迟从星地距离到传播时延的完整链路4.1 传播延迟拆解不是只有距离除以光速通信延迟的计算看起来简单单程传播延迟星地距离 / 光速。但真实系统里这个延迟要拆成四层来看第一层是自由空间传播延迟LEO 卫星在仰角 10° 时星地距离通常在 15002500 km 之间取决于轨道高度单程延迟约 5.08.3 ms卫星过顶时距离最短约等于轨道高度减地球遮蔽的部分。第二层是电离层和对流层的附加延迟电离层对 L 频段信号造成的群延迟可达 110 ns 量级对纳秒级要求的同步协议不能忽略。第三层是地面站和卫星端的处理时延包括调制解调、编解码、协议处理通常在几毫秒到几十毫秒。第四层是排队和调度时延由网络负载决定。对大多数可见性仿真场景后三层要么是常量要么依赖具体设备能在仿真中精确建模的核心就是第一层——自由空间传播延迟它直接由轨道动力学决定。4.2 在可见窗口内计算并绘制通信延迟曲线接着第 3 节的代码直接复用altaz()返回的distance数组用光速做除法import numpy as np import matplotlib.pyplot as plt # 复用前面的 satellite / station / times / alt / az / distance c 299792.458 # 光速 km/s # 单程传播延迟 one_way_delay distance.km / c # 单位: 秒 min_elev 10 visible_mask alt.degrees min_elev # 只保留可见区间的延迟数据 vis_times times.utc_datetime() vis_delay one_way_delay * 1000 # 转成毫秒 # 找出每个可见窗口的边界索引 window_starts [] window_ends [] in_window False for i, v in enumerate(visible_mask): if v and not in_window: window_starts.append(i) in_window True elif not v and in_window: window_ends.append(i) in_window False if in_window: window_ends.append(len(visible_mask)) for idx, (s, e) in enumerate(zip(window_starts, window_ends)): # 找窗口内延迟最小即距离最短的时刻 min_delay_idx s np.argmin(one_way_delay[s:e]) print(f窗口{idx1}: {vis_times[s].strftime(%H:%M:%S)} ~ {vis_times[e-1].strftime(%H:%M:%S)}, f时长 {(e-s)*5}秒, 最小延迟 {vis_delay[min_delay_idx]:.3f} ms) plt.figure(figsize(10, 5)) plt.plot(vis_times, vis_delay, b-, linewidth1) plt.axhline(yone_way_delay.min()*1000, colorr, linestyle--, alpha0.7) plt.xlabel(UTC 时间) plt.ylabel(单程传播延迟 (ms)) plt.title(ISS 相对北京地面站的传播延迟随时间变化) plt.grid(True, alpha0.3) plt.tight_layout() plt.show()这段代码把可见性和延迟数据合到一起做了两件事一是定位每个可见窗口输出窗口起止时间和窗口内的最小延迟二是画出延迟随时间变化的曲线。曲线形状是典型的“U 形”——卫星从地平线升起时延迟最大过顶时达到最小值然后再次增大。这个 U 形曲线直接决定了链路预算中最差工况和最优工况的差距。注意distance.km是skyfield返回的星地直线距离它包含了地球曲率的影响。如果把距离简单想成卫星高度会低估实际延迟仰角 10° 时400 km 轨道的卫星实际距离约 1400 km对应延迟约 4.7 ms而不是 1.3 ms。这就是为什么延迟计算必须和可见性计算耦合在一起做不能分开估算。4.3 往返延迟对协议设计的意义上面算的是单程延迟实际通信协议关心的是往返时间RTT。LEO 场景的 RTT 包括上行单程延迟加下行单程延迟如果上下行都走同一颗卫星RTT 就是两倍单程延迟。以 10° 最小仰角、550 km 轨道高度为例单程延迟约 5.4 msRTT 约 10.8 ms这远低于 GEO 卫星的约 500 ms RTT是 LEO 星座在低延迟场景的核心卖点。但你做仿真时要注意实际端到端延迟还会经过地面骨干网和星座星间链路。星间链路ISL在极地轨道星座里通常跨越几百到几千公里每跳增加 1.36.7 ms 延迟。如果一颗卫星不能直接连接到目标地面站需要经过 23 跳 ISL端到端延迟可能逼近 3040 ms这对实时语音仍有明显感知。所以延迟计算要区分“星地传播延迟”和“端到端延迟”别拿一个数覆盖所有场景。5. 完整模拟示例找出延迟最小的可用窗口把前面三部分串起来做一个完整的 LEO 卫星过境模拟。场景设定地面站位于北京39.9042°N, 116.4074°E目标卫星是 ISS轨道高度约 420 km倾角 51.6°最小工作仰角 10°模拟时长 24 小时采样间隔 5 秒。目标是从所有可见窗口里找出延迟最小的一次过境并输出该窗口内延迟随时间的变化和对应的星下点轨迹。from skyfield.api import load, wgs84 from datetime import timedelta import numpy as np # 1. 加载 TLE 和地面站 satellite load.tle(https://celestrak.org/NORAD/elements/gp.php?CATNR25544) station wgs84.latlon(39.9042, 116.4074, elevation_m50) ts load.timescale() # 2. 生成24小时采样点 start ts.now() sample_count 24 * 3600 // 5 # 24小时每5秒一个采样 times ts.utc([start.utc_datetime() timedelta(secondsi*5) for i in range(sample_count)]) # 3. 计算仰角和距离 difference satellite - station topocentric difference.at(times) alt, az, distance topocentric.altaz() c 299792.458 one_way_delay_ms distance.km / c * 1000 min_elev 10.0 visible alt.degrees min_elev # 4. 找到所有可见窗口并评估 windows [] in_win False for i in range(len(visible)): if visible[i] and not in_win: win_start i in_win True elif not visible[i] and in_win: win_end i - 1 win_mid win_start np.argmin(one_way_delay_ms[win_start:win_end1]) windows.append({ start: times[win_start].utc_strftime(%Y-%m-%d %H:%M:%S), end: times[win_end].utc_strftime(%H:%M:%S), duration_s: (win_end - win_start) * 5, min_delay_ms: one_way_delay_ms[win_mid], max_elev_deg: alt.degrees[win_start:win_end1].max(), mid_time: times[win_mid].utc_strftime(%H:%M:%S) }) in_win False if in_win: win_end len(visible) - 1 win_mid win_start np.argmin(one_way_delay_ms[win_start:win_end1]) windows.append({ start: times[win_start].utc_strftime(%Y-%m-%d %H:%M:%S), end: times[win_end].utc_strftime(%H:%M:%S), duration_s: (win_end - win_start) * 5, min_delay_ms: one_way_delay_ms[win_mid], max_elev_deg: alt.degrees[win_start:win_end1].max(), mid_time: times[win_mid].utc_strftime(%H:%M:%S) }) # 5. 排序并输出 windows.sort(keylambda w: w[min_delay_ms]) print(f24小时内共检测到 {len(windows)} 个可见窗口) for w in windows[:5]: print(f开始 {w[start]} | 结束 {w[end]} | 时长 {w[duration_s]:4d}s | f最小延迟 {w[min_delay_ms]:.2f}ms | 最大仰角 {w[max_elev_deg]:.1f}° | f延迟最小点 {w[mid_time]})运行时注意两点。第一load.tle()需要联网下载 TLE 数据如果网络不可用先前把 TLE 文本存成.txt文件用load.tle_file(iss.txt)读取。第二load.tle()默认会缓存数据但 TLE 数据本身有时效性——超过一周的 TLE 预报误差会显著增大如果仿真结果里窗口起止时间和实际不符先检查 TLE 的 epoch 日期是否太老。运行结果通常会发现几个规律可见窗口持续时间在 310 分钟之间最大仰角越高的窗口持续越长、延迟越小24 小时内同一地面站对单颗 LEO 卫星的可见窗口大约 46 次这是因为 LEO 卫星轨道周期约 90 分钟但受地球自转和轨道面进动的影响不是每次过境都恰好经过地面站上空。6. 采样间隔和 TLE 时效两个直接影响仿真可信度的参数最后一个值得单独拿出来说的参数是时间采样间隔。前面所有计算都基于离散采样采样间隔决定了两个关键指标的上限窗口边界定位精度和最小延迟估计精度。如果采样间隔是 30 秒而一次可见窗口只有 5 分钟那么窗口内最多只有 10 个采样点窗口起止时刻的误差可能高达 ±30 秒对应的仰角误差可能达到几度直接影响“最小工作仰角是否满足”的判断。我一般会按轨道周期和仰角变化率来选采样间隔。LEO 卫星在过顶时仰角变化率最快大约每秒 0.5°1°如果要求仰角判定误差小于 0.5°采样间隔应控制在 1 秒以内如果只是粗筛窗口510 秒间隔足够。延迟方面过顶时距离变化率约 0.51 km/s5 秒采样间隔下最小延迟的误差约 0.01 ms对工程判断无影响所以延迟精度不构成采样间隔的约束条件——仰角精度才是。另一个容易忽视的是 TLE 的预报误差边界。SGP4 模型的误差随预报时长非线性增长第 1 天误差约 13 km第 3 天可能到 510 km第 7 天可能超过 20 km。对可见性计算10 km 的位置误差在仰角 30° 以上时只会带来零点几度的误差但在仰角接近最小阈值如 10°时误差会被放大——因为此时视线接近地球切面同样的位置误差对应更大的仰角变化。所以做高精度过境预报时TLE 的年龄不要超过 3 天超过之后窗口起止时间的预测误差可能达到分钟级。验证计算正确性可以从两个方向入手一是对比 CelesTrak 网站上的过境预报二是用已知卫星做自检——比如以中国空间站或某颗铱星为对象手动检查一个窗口内最大仰角时刻是否和星下点最接近地面站的时刻重合。如果想把 SGP4 输出的 TEME 坐标精度再往上提可以用sgp4库的WGS72和WGS84重力场常数切换功能对于低轨卫星两个模型下的位置差异通常在百米量级只有当你的仰角阈值正好落在可见性边界附近时这种细微差异才可能改变窗口的判定结果。本文还有配套的精品资源点击获取

用Python和SGP4计算低轨卫星可见性与通信延迟全流程

相关推荐

YOLOv8玉米叶病害检测实战：数据集训练与PyQt界面部署

C51单片机控制ISD1820PY语音录放模块：原理图、接线与代码详解

PSO-DV-Hop无线传感器网络定位算法详解与MATLAB实现

更多资讯

Unity消消乐开发：菜单界面与游戏界面场景解耦实战

yolov5垃圾分类识别检测：数据、训练到部署全流程实战

玩Steam游戏用什么开黑语音软件？2026兼容全平台的语音工具盘点

硕士论文高效写作四步法：从选题到终稿全流程解析

Python机器学习股市情感分析全流程：从词向量到量化因子

CopilotKit + A2A + A2UI 实战：用 AG-UI 协议构建可动态渲染 UI 的餐厅预订 Agent

今日更新

Linux 内核 ARM 虚拟内存布局全解析：从 4GB 地址空间划分到源码级验证

AIGC降重工具测评与教育应用指南

数据容灾核心指标与实战方案解析

本周精选

超人会飞不算本事：系统稳定依赖清晰规则与边界设计

超人VS蜘蛛侠：拆解超级IP的影响力与传播方法论

基于CNN的调制信号识别：MATLAB实现时频图分类实战

本月热门

基于YOLO的焊缝缺陷检测：从数据标注到工业部署全流程实践

分布式服务的决策记录

创业产品的首版范围

读完文章，想聊聊您的网站？