用 Python + TLE 数据,打造交互式卫星轨迹可视化

from astropy.utils import iersiers.conf.auto_download = Trueiers.conf.iers_auto_url = 'https://datacenter.iers.org/data/9/finals2000A.all'iers.conf.auto_max_age = 30

sat = Satrec.twoline2rv(line1, line2)jd, fr = jday(dt.year, dt.month, dt.day, dt.hour, dt.minute, dt.second)error, pos, vel = sat.sgp4(jd, fr)

teme_cart = CartesianRepresentation(pos_teme * u.km)teme = TEME(teme_cart, obstime=obstime)itrs = teme.transform_to(ITRS(obstime=obstime))

def compute_elevation(r_ecef, station):    r_sta = np.array([station.x.to(u.km).value, ...])    d = r_ecef - r_sta                           # 卫星相对站的矢量    up = r_sta / np.linalg.norm(r_sta)           # 当地"上"方向    cos_zenith = np.dot(d / np.linalg.norm(d), up)    return 90.0 - np.rad2deg(np.arccos(cos_zenith))

ecef_arr = np.zeros((num_frames, 3))elevation_dict = {sname: np.zeros(num_frames) for sname in stations}

frame_data.append(dict(    type='scatter3d',   # 必须声明，否则 Plotly 默认按 scatter 处理    x=[pos[0]], y=[pos[1]], z=[pos[2]],    marker=dict(...)))

import numpy as npimport pandas as pdfrom datetime import datetime, timedeltafrom astropy.time import Timefrom astropy.coordinates import (    ITRS, TEME, EarthLocation, CartesianRepresentation)import astropy.units as ufrom sgp4.api import Satrec, jdayimport plotly.graph_objects as gofrom plotly.subplots import make_subplots# ==================== 地球定向参数（必须在其他 astropy 导入前配置） ====================from astropy.utils import iersfrom astropy.utils.iers import IERS_Auto, IERS_B# 启用自动下载，保留高精度iers.conf.auto_download = Trueiers.conf.iers_auto_url = 'https://datacenter.iers.org/data/9/finals2000A.all'iers.conf.auto_max_age = 30  # 缓存有效期 30 天# 显式初始化：下载 IERS_A，失败则回退 IERS_Btry:    iers.IERS.iers_table = IERS_Auto.open()    print('IERS-A 高精度表格已加载')except Exception as e:    print(f'IERS-A 下载失败，回退内置 IERS-B: {e}')    iers.IERS.iers_table = IERS_B.open()# ==================== 常数 ====================R_EARTH = 6371.0                     # kmELEVATION_THRESHOLD = 5.0            # 可视仰角阈值 (°)# ==================== 多地面站定义 ====================stations = {    'Beijing': EarthLocation(lat=39.9 * u.deg, lon=116.4 * u.deg, height=0 * u.m),    'Shanghai': EarthLocation(lat=31.2 * u.deg, lon=121.4 * u.deg, height=0 * u.m),    'NewYork': EarthLocation(lat=40.7 * u.deg, lon=-74.0 * u.deg, height=0 * u.m),}# ==================== TLE 数据（国际空间站 + 示例卫星） ====================# 注意：可替换为最新的 TLE，从 celestrak.org 获取TLE_SETS = [    {        'name': 'ISS (ZARYA)',        'line1': '1 25544U 98067A   24001.50000000  .00001234  00000+0  23456-4 0  9991',        'line2': '2 25544  51.6420 123.4567 0001234  45.6789 314.5678 15.50123456123456',        'color': 'red'    },    {        'name': 'CSS (TIANHE)',        'line1': '1 48274U 21035A   24001.50000000  .00000567  00000+0  12345-4 0  9992',        'line2': '2 48274  41.4690  78.1234 0001567  89.1234 270.9876 15.61234567890123',        'color': 'cyan'    }]# ==================== SGP4 传播与坐标转换 ====================def tle_to_eci(sat, dt):    """通过 SGP4 计算卫星在 dt 时刻的 ECI (TEME) 位置 (km)"""    jd, fr = jday(dt.year, dt.month, dt.day, dt.hour, dt.minute, dt.second)    error, pos, vel = sat.sgp4(jd, fr)    if error != 0:        raise RuntimeError(f"SGP4 error at {dt}: {error}")    return np.array(pos)   # TEME 坐标系def teme_to_itrs(pos_teme, obstime):    """TEME -> ITRS (ECEF) 高精度转换"""    teme_cart = CartesianRepresentation(pos_teme * u.km)    teme = TEME(teme_cart, obstime=obstime)    itrs = teme.transform_to(ITRS(obstime=obstime))    return np.array([itrs.x.to(u.km).value,                     itrs.y.to(u.km).value,                     itrs.z.to(u.km).value])def compute_elevation(r_ecef, station):    """简单向量几何计算仰角（度），无需 astropy 变换链，速度极快"""    # 地面站 ECEF 坐标 (km)    r_sta = np.array([station.x.to(u.km).value,                       station.y.to(u.km).value,                       station.z.to(u.km).value])    # 卫星相对地面站的矢量    d = r_ecef - r_sta    # 地面站当地"上"方向（简化为地球球心径向）    up = r_sta / np.linalg.norm(r_sta)    # 仰角 = 90° − 天顶角    d_norm = np.linalg.norm(d)    if d_norm < 1e-6:        return 90.0    cos_zenith = np.dot(d / d_norm, up)    # clamp 防止浮点溢出    cos_zenith = max(-1.0, min(1.0, cos_zenith))    return 90.0 - np.rad2deg(np.arccos(cos_zenith))# ==================== 时间轴构建 ====================# 传播 24 小时，200 帧start_time = datetime(2023, 12, 1, 0, 0, 0)num_frames = 200duration = timedelta(hours=24)time_steps = [start_time + duration * i / num_frames for i in range(num_frames)]# ==================== 主计算循环 ====================all_sat_data = []   # 存放每颗卫星的完整数组for tle in TLE_SETS:    # 初始化 SGP4 卫星对象    sat = Satrec.twoline2rv(tle['line1'], tle['line2'])    name = tle['name']    color = tle['color']    # 预分配数组    ecef_arr = np.zeros((num_frames, 3))    lons = np.zeros(num_frames)    lats = np.zeros(num_frames)    # 为每个地面站存储仰角数组    elevation_dict = {sname: np.zeros(num_frames) for sname in stations}    times_out = []    for idx, t in enumerate(time_steps):        # SGP4 得到 TEME 位置        pos_teme = tle_to_eci(sat, t)        obstime = Time(t)        # TEME -> ITRS (ECEF)        r_ecef = teme_to_itrs(pos_teme, obstime)        ecef_arr[idx] = r_ecef        # 经纬度        x, y, z = r_ecef        lat = np.rad2deg(np.arctan2(z, np.sqrt(x**2 + y**2)))        lon = np.rad2deg(np.arctan2(y, x))        lats[idx] = lat        lons[idx] = lon        # 所有地面站仰角        for sname, sta in stations.items():            elev = compute_elevation(r_ecef, sta)            elevation_dict[sname][idx] = elev        times_out.append(t)    # 保存这颗卫星的所有数据    sat_data = {        'name': name,        'color': color,        'times': times_out,        'ecef': ecef_arr,        'lons': lons,        'lats': lats,        'elevations': elevation_dict    }    all_sat_data.append(sat_data)    # ---- 输出 CSV（每颗卫星一个文件） ----    df = pd.DataFrame({        'Time(UTC)': [t.strftime('%Y-%m-%d %H:%M:%S') for t in times_out],        'ECEF_X_km': ecef_arr[:, 0],        'ECEF_Y_km': ecef_arr[:, 1],        'ECEF_Z_km': ecef_arr[:, 2],        'Longitude_deg': lons,        'Latitude_deg': lats    })    for sname in stations:        df[f'Elevation_{sname}_deg'] = elevation_dict[sname]    csv_name = f"{name.replace(' ', '_')}_trajectory.csv"    df.to_csv(csv_name, index=False, float_format='%.6f')    print(f"已保存 {csv_name} (含{len(stations)}个地面站仰角)")# ==================== 地球网格 ====================def sphere_mesh(radius, steps=60):    phi = np.linspace(0, np.pi, steps)    theta = np.linspace(0, 2*np.pi, steps)    phi, theta = np.meshgrid(phi, theta)    x = radius * np.sin(phi) * np.cos(theta)    y = radius * np.sin(phi) * np.sin(theta)    z = radius * np.cos(phi)    return x, y, zx_e, y_e, z_e = sphere_mesh(R_EARTH)# ==================== 构建 Plotly 图形 ====================fig = make_subplots(    rows=1, cols=2,    specs=[[{'type': 'scatter3d'}, {'type': 'scattergeo'}]],    subplot_titles=('真实卫星三维空间轨迹 (TLE + SGP4)', '星下点轨迹与地面站'),    column_widths=[0.65, 0.35])# ---- 左图：3D 场景 ----fig.add_trace(go.Surface(    x=x_e, y=y_e, z=z_e,    colorscale='Blues', opacity=0.4,    showscale=False, name='地球'), row=1, col=1)scatter_indices = []   # 记录每颗卫星当前位置 trace 的索引（用于动画更新）for sat_data in all_sat_data:    color = sat_data['color']    name = sat_data['name']    # 空间轨迹线    fig.add_trace(go.Scatter3d(        x=sat_data['ecef'][:, 0], y=sat_data['ecef'][:, 1],        z=sat_data['ecef'][:, 2],        mode='lines', line=dict(color=color, width=2),        name=f"{name} 空间轨迹"    ), row=1, col=1)    # 星下点投影线    norms = np.linalg.norm(sat_data['ecef'], axis=1)    ground = R_EARTH * sat_data['ecef'] / norms[:, np.newaxis]    fig.add_trace(go.Scatter3d(        x=ground[:, 0], y=ground[:, 1], z=ground[:, 2],        mode='lines', line=dict(color=color, width=1.5, dash='dot'),        name=f"{name} 星下点投影"    ), row=1, col=1)    # 卫星当前位置（初始第一帧）    trace_idx = len(fig.data)    fig.add_trace(go.Scatter3d(        x=[sat_data['ecef'][0, 0]], y=[sat_data['ecef'][0, 1]],        z=[sat_data['ecef'][0, 2]],        mode='markers', marker=dict(size=8, color=color, symbol='circle'),        name=f"{name} 当前位置"    ), row=1, col=1)    scatter_indices.append(trace_idx)# ---- 右图：2D 地图 + 地面站 ----# 1) 每颗卫星的完整星下点轨迹（实线）for sat_data in all_sat_data:    fig.add_trace(go.Scattergeo(        lon=sat_data['lons'], lat=sat_data['lats'],        mode='lines', line=dict(width=2, color=sat_data['color']),        name=f"{sat_data['name']} 星下点"    ), row=1, col=2)# 2) 地面站标记station_colors = ['gold', 'orange', 'magenta']  # 每个站一个颜色for (sname, sta), sc in zip(stations.items(), station_colors):    fig.add_trace(go.Scattergeo(        lon=[sta.lon.deg], lat=[sta.lat.deg],        mode='markers+text',        marker=dict(size=10, color=sc, symbol='star', line=dict(width=1, color='black')),        text=[sname],        textposition='top center',        name=f'地面站 {sname}'    ), row=1, col=2)# 地图风格（显式配置 geo2 子图，确保可旋转缩放）fig.update_layout(    geo2=dict(        projection=dict(            type='orthographic',            rotation=dict(lon=100, lat=30, roll=0)        ),        showcoastlines=True, coastlinecolor='Black',        showland=True, landcolor='lightgray',        showocean=True, oceancolor='lightblue',    ))# ==================== 动画帧 ====================frames = []for frame_idx in range(num_frames):    frame_data = []    for sat_idx, sat_data in enumerate(all_sat_data):        pos = sat_data['ecef'][frame_idx]        # 使用对 Beijing 站的仰角决定标记形状（也可以选择一个参考站）        elev_beijing = sat_data['elevations']['Beijing'][frame_idx]        marker_symbol = 'circle' if elev_beijing >= ELEVATION_THRESHOLD else 'x'        frame_data.append(dict(            type='scatter3d',            x=[pos[0]],            y=[pos[1]],            z=[pos[2]],            marker=dict(size=10, color=sat_data['color'], symbol=marker_symbol)        ))    frames.append(go.Frame(        data=frame_data,        traces=scatter_indices,        name=f"t{frame_idx}"    ))fig.frames = frames# 播放控件与时间滑块fig.update_layout(    updatemenus=[dict(        type='buttons',        buttons=[            dict(label='播放', method='animate',                 args=[None, {'frame': {'duration': 80, 'redraw': True},                              'fromcurrent': True}]),            dict(label='暂停', method='animate',                 args=[[None], {'frame': {'duration': 0, 'redraw': True},                                'mode': 'immediate', 'transition': {'duration': 0}}])        ],        showactive=False,        x=0.1, y=0, xanchor='right', yanchor='top'    )],    sliders=[dict(        active=0,        yanchor='top', xanchor='left',        currentvalue={'prefix': 'UTC: '},        pad={'b': 10, 't': 50},        len=0.9, x=0.1,        steps=[dict(method='animate',                    args=[[f't{k}'], {'frame': {'duration': 80, 'redraw': True},                                      'mode': 'immediate', 'transition': {'duration': 0}}],                    label=time_steps[k].strftime('%H:%M'))               for k in range(0, num_frames, 5)]    )],    scene=dict(aspectmode='data',               xaxis_title='X ECEF (km)',               yaxis_title='Y ECEF (km)',               zaxis_title='Z ECEF (km)',               dragmode='turntable',               camera=dict(eye=dict(x=1.5, y=1.5, z=1.2))),    title='TLE 驱动多卫星轨迹与地面站可见性',    height=750)# ==================== 保存为 HTML ====================# output_html = 'satellite_animation.html'# fig.write_html(output_html)# print(f"交互式动画已保存至 {output_html}，可直接用浏览器打开。")fig.show(config=dict(scrollZoom=True, displayModeBar=True, displaylogo=False))# 若仍希望在 notebook 中显示，可取消下一行注释# fig.show()

pip install numpy pandas astropy sgp4 plotlypython tLETo3D-2D.py