简介基于随机几何的低轨星座下行链路仿真与分析完整资料面向具备Python编程基础的卫星通信研究人员与系统工程师重点解决低轨星座建模、星地链路损耗计算及多星干扰评估问题。其中以二项点过程(BPP)构建星座模型深入分析单星与多星场景下的路径损耗和同频干扰提出干扰期望计算方法可用于巨型低轨星座性能分析与星地链路设计参考。配套代码覆盖参数设置、卫星与地面站随机生成、自由空间及大气损耗、接收功率与噪声功率计算、SINR分析、结果可视化等完整环节代码可运行且附带逐段解释方便读者对照学习。压缩包内为单个docx文档共1个文件大小约51KB轻量便携目前已有105人学习/下载适合作为低轨卫星通信仿真的入门与实践参考。内容还展望了天地一体化网络、大规模星座干扰建模等未来方向为后续研究提供延伸切入点。1. 为什么低轨星座下行链路仿真要改用随机几何低轨星座规模越做越大几十上百颗卫星只是起步千星级别的系统已经开始规划。传统仿真做法信赖星历文件加轨道积分器逐颗卫星计算星下点、仰角、路径损耗精度确实高但代价是算力需求随卫星数量膨胀——分析一条下行链路要遍历全部可见卫星换一个用户位置就得重算一轮。随机几何的思路完全不同它把卫星在观察时刻的空间分布当作随机点过程的一次实现用强度函数近似星座的整体几何直接推导SINR、覆盖概率、中断概率的表达式或半解析结果。这套方法牺牲掉单颗卫星的精确轨迹换来参数尺度上的可解释性轨道高度、卫星数量、发射功率、频率复用因子如何影响覆盖率几分钟内就能跑完。本文就用这套思路从星座建模开始逐步实现完整的低轨卫星下行链路仿真。2. 低轨星座建模从Walker构型到随机点过程的映射2.1 Walker星座参数与初始轨道生成低轨星座最常见的构型是Walker Delta星座四个参数就能描述T是卫星总数P是轨道面数F是相位因子i是轨道倾角。例如常见的参考方案T288、P18、F1、i53°每面16颗卫星近圆形轨道高度550km。相位因子F控制相邻轨道面间卫星的相位偏移量大小为2πF/T这个参数直接影响面间相邻卫星的最小间距进而决定覆盖均匀度和同频干扰的分布特性。生成星座初始ECI位置的最小Python实现如下import numpy as np def generate_walker_constellation(T288, P18, F1, i_deg53.0, h_km550.0): 生成Walker Delta星座的ECI位置向量圆轨道假设 T: 卫星总数 P: 轨道面数 F: 相位因子 i_deg: 轨道倾角度 h_km: 轨道高度km 返回卫星ECI坐标数组形状 (T, 3) R_E 6371.0e3 # 地球平均半径m r h_km * 1e3 R_E # 轨道半径m S T // P # 每面卫星数 sats np.zeros((T, 3)) idx 0 for p in range(P): raan 2 * np.pi * p / P # 升交点赤经均分 for s in range(S): # 面内相角等间隔 相位因子偏移 nu 2 * np.pi * s / S 2 * np.pi * F * p / T # 轨道平面内的位置圆轨道 x_orb r * np.cos(nu) y_orb r * np.sin(nu) # 1) 绕X轴旋转倾角 cos_i, sin_i np.cos(np.radians(i_deg)), np.sin(np.radians(i_deg)) x1, y1, z1 x_orb, y_orb * cos_i, y_orb * sin_i # 2) 绕Z轴旋转升交点赤经 cos_r, sin_r np.cos(raan), np.sin(raan) sats[idx, 0] x1 * cos_r - y1 * sin_r sats[idx, 1] x1 * sin_r y1 * cos_r sats[idx, 2] z1 idx 1 return sats逻辑说明先按轨道面遍历每个面的升交点赤经RAAN在360°内等间隔分布再在面内按相角等间隔放置S颗卫星并叠加相位因子偏移。圆轨道假设下真近点角就等于相角省去了开普勒方程求解。旋转矩阵分两步执行先绕X轴旋转轨道倾角把轨道平面从赤道平面倾斜到指定角度再绕Z轴旋转RAAN把升交点对准正确方向。参数说明轨道高度h_km决定轨道周期和覆盖半径550km对应的轨道周期约95分钟覆盖圆半径约2400km按最小仰角10°计算。T和P的比例决定了面内卫星间距P越多单面卫星越少面间相对运动越复杂。F取0到P-1之间的整数工程上倾向选F1或F2使相邻轨道面的卫星相位错开避免形成规则的干扰纹理。这里的地球半径用平均半径6371km如果要做特定真实星座的精确复现需要换成WGS84椭球并计入J2摄动引起的RAAN长期漂移。提示Walker参数T288、P18、F1适用于系统级覆盖预研。如果你要复现某个具体星座的论文结果先确认论文用的地面覆盖最小仰角这个值对可见卫星数量影响极大。2.2 从离散卫星到随机点过程强度函数的建模逻辑有了离散卫星的位置为什么还要引入随机几何核心原因是计算复杂度的量级差异。假设T288颗卫星地面用户在任意时刻平均可见约14颗卫星550km高度、10°仰角约束下遍历计算完全可以接受。但当星座规模到10000颗时平均可见卫星数接近500颗逐颗精确计算距离、天线增益和仰角的时间开销随T线性增长而蒙特卡洛仿真又需要成千上万次重复整体开销迅速变得难以承受。随机几何的做法是用泊松点过程PPP近似卫星的空间分布把问题从「枚举每颗卫星」转成「对空间积分」。核心思想是用户可见的球冠区域内卫星位置被建模为独立均匀分布的随机点密度λ T / (4πR²)R是轨道半径。可见球冠面积由轨道高度和最小仰角决定平均可见卫星数为λ乘以球冠面积。这个近似在卫星数量大、最小仰角高、轨道高度低的时候最准因为可见区域收缩后边界效应减弱区域内分布更接近均匀独立。下表给出了不同参数组合下的近似适用性判断直接决定你后边仿真用确定性方法还是概率方法参数典型范围对PPP近似的影响仿真建议卫星总数T100~10000T越大泊松近似越准T300可逐颗精确计算最小仰角10°~40°仰角越高可见区域越小近似越可靠10°~20°适合PPP验证轨道高度340~1200km高度越低覆盖半径越小边界效应越弱340km处PPP误差最小轨道面数P12~24P越少轨道纹理越明显偏离PPP越大P12建议用确定性模型关键是要意识到PPP不是对真实星座的精确描述而是一种把「卫星在可见区域里以某密度随机出现」当作等效模型的工程近似。低轨星座的轨道面如果只有几个卫星分布会呈现明显的带状纹理此时PPP会低估干扰的空间聚集性。一般做法是先在确定性Walker星座上验证一次再切到PPP模型做参数扫描两种方法互为校验。2.3 可见性判断与用户-卫星几何计算无论用哪种模型仿真里都要计算用户到一颗卫星的仰角。定义用户位置为经纬高卫星位置为ECI坐标需要把卫星位置转到用户本地东北天ENU坐标系才能算出仰角。下面的函数实现了这个几何过程def compute_elevation_angle(user_lla, sat_eci): 计算用户到卫星的仰角 user_lla: (纬度deg, 经度deg, 高度m) sat_eci: 卫星ECI坐标 (x, y, z) m 返回仰角度 lat np.radians(user_lla[0]) lon np.radians(user_lla[1]) alt user_lla[2] # WGS84椭球参数 a 6378137.0 # 长半轴m e2 6.69437999014e-3 # 偏心率平方 N a / np.sqrt(1 - e2 * np.sin(lat)**2) # 用户ECEF坐标 ux (N alt) * np.cos(lat) * np.cos(lon) uy (N alt) * np.cos(lat) * np.sin(lon) uz (N * (1 - e2) alt) * np.sin(lat) # 卫星到用户的向量 dx sat_eci[0] - ux dy sat_eci[1] - uy dz sat_eci[2] - uz # ECEF到ENU旋转矩阵 sin_lat, cos_lat np.sin(lat), np.cos(lat) sin_lon, cos_lon np.sin(lon), np.cos(lon) e -sin_lon * dx cos_lon * dy n -sin_lat * cos_lon * dx - sin_lat * sin_lon * dy cos_lat * dz u cos_lat * cos_lon * dx cos_lat * sin_lon * dy sin_lat * dz return np.degrees(np.arctan2(u, np.sqrt(e**2 n**2)))逻辑说明先把用户经纬高通过WGS84椭球转成ECEF直角坐标然后用ECEF到ENU的方向余弦矩阵把卫星相对用户的向量分解到东、北、天三个方向。仰角是天向分量与水平分量的反正切水平分量由东向和北向的平方根给出。参数说明最小仰角阈值通常取10°或25°对应不同的覆盖概率要求。10°仰角可以看到更远的低轨卫星但信号穿过的对流层路径更长大气衰减和闪烁更严重25°仰角服务质量更稳定但可见卫星数量减少约30%。在干扰分析里这个选择直接影响干扰源的密度建议在参数表中同时输出两个仰角阈值下的结果。3. 下行链路损耗计算从自由空间到大气衰减3.1 自由空间路径损耗与多普勒频移低轨卫星下行链路损耗中自由空间路径损耗占绝对主导。公式为PL 20log10(4πd/λ)其中d是用户到卫星的直线距离。对20GHz载频、550km轨道高度、仰角10°时的最大斜距约2400km来算路径损耗约为186dB左右仰角90°卫星在头顶时斜距等于轨道高度550km路径损耗约173dB。也就是说仅仰角变化就能引入13dB的链路预算波动。多普勒频移方面低轨卫星轨道速度约7.6km/s径向相对速度在地面用户看来的极端情况可达7km/s量级。20GHz承载频率下最大多普勒频移可接近500kHz远超地面移动通信系统通常面临的规模。仿真相干时间内如果只取瞬时快照多普勒不会影响SINR但若仿真时间窗口拉长必须计算每颗卫星相对用户的径向速度并折算到频偏否则接收机频率补偿误差会变成实际损耗。def free_space_path_loss(d_m, freq_hz): 自由空间路径损耗d_m为斜距米freq_hz为载波频率 c 2.99792458e8 return 20 * np.log10(4 * np.pi * d_m * freq_hz / c) def doppler_shift(sat_vel, user_pos, sat_pos, freq_hz): 计算多普勒频移sat_vel为卫星速度矢量(m/s) c 2.99792458e8 r_vec sat_pos - user_pos r_dist np.linalg.norm(r_vec) radial_v np.dot(sat_vel, r_vec) / r_dist return radial_v * freq_hz / c逻辑说明多普勒频移的代码里关键一步是把卫星速度矢量投影到用户-卫星连线上得到径向速度。径向速度为正表示卫星正在靠近用户频移为正为负表示远离频移为负。LEO卫星过顶过程中多普勒由正到负的过零变化是常见特征。参数说明载波频率越高多普勒越显著这也是Ka频段比L频段更难做低轨通信的原因之一。如果你仿真的是20GHz下行建议把子载波间隔和频偏补偿算法纳入链路级仿真但在系统级信噪比分析中只需记录最近服务卫星的多普勒值用于判断是否需要引入频偏校正开销。3.2 大气衰减与雨衰模型自由空间损耗之上电波穿过对流层时还有气体吸收和雨衰。低频段L/S的大气衰减可以忽略到了Ka/Ku频段氧气和水蒸气的分子吸收不可再忽略。ITU-R P.676给出了天顶方向大气衰减的标准计算方法。在20GHz、天顶方向氧气衰减约0.2dB水蒸气约0.1dB但仰角降到10°时电波在大气层中的路径长度增加约5倍衰减合计可超过1.5dB。雨衰是另一个大头用ITU-R P.618推荐的模型估算。核心输入是区域降雨强度R_0.01每年超过0.01%时间的降雨率单位为mm/h以及仰角、极化倾角和频率。一个粗略的参考值在K频段中纬度地区0.01%时间降雨强度约30mm/h对应20GHz、30°仰角下雨衰约3-6dB。频段频率范围天顶大气衰减0.1%时间雨衰参考值适用场景L1-2GHz可忽略几乎可忽略星间链路、物联网S2-4GHz可忽略0.5dB移动通信备份Ku12-18GHz0.3-0.8dB1-3dB传统广播Ka20-30GHz0.5-1.5dB3-10dB宽带低轨星座仿真里处理雨衰有一个常见选择直接用概率分布注入而不是精确计算每次降雨的时空分布。最简单的方式是假设雨衰服从对数正态分布在每次蒙特卡洛迭代里独立抽样叠加。3.3 接收功率与链路预算表把上述损耗项组合起来就得到下行链路预算。接收到的载波功率可以写成P_r P_t G_t - L_fs - L_atm - L_rain G_r其中P_t为卫星发射功率G_t为卫星天线增益L_fs为自由空间损耗L_atm为大气气体衰减L_rain为雨衰G_r为用户终端天线增益。以下是一组典型参数下的链路预算参数数值计算说明卫星发射功率10W40dBm下行通常比上行功率充裕卫星天线增益30dBi多波束相控阵单波束增益EIRP70dBm40 30自由空间损耗173dB550km垂直高度、20GHz大气衰减0.8dB仰角30°、中纬度雨衰0.1%时间2.5dB20GHz中纬度中等降雨用户天线增益25dBi固定终端接收功率-81.3dBm70 - 173 - 0.8 - 2.5 25噪声功率-104dBm500MHz带宽噪声系数1.5dBSNR22.7dB接收功率减噪声功率这个表里SNR的余量很大因为系统级干扰分析的主要矛盾在于同频干扰而不是噪声。低轨星座通常全频复用服务卫星以外还有成百上千颗卫星落在用户天线波束内它们造成的集总干扰要比底噪高一个量级。所以下一章要解决的核心问题是在那么多同频卫星里SINR还剩下多少。4. 干扰分析与下行链路SINR仿真4.1 同频干扰的随机几何建模考虑一个下行场景一颗服务卫星向目标用户发送信号其他同频卫星的信号同时到达用户天线形成集总干扰。在传统确定性仿真里你需要得到每一颗干扰卫星的精确位置、天线增益、发射功率然后枚举求和。随机几何的建模则把干扰卫星集合视为可见球冠区域内的PPP密度为λ_int α·λ其中α是频率复用系数全复用1。为什么PPP在这里是合理的近似低轨星座干扰源数量大、距离散布广、各干扰链路间的相互依赖较弱。用户天线通常有指向性只有落在主波束范围内的干扰卫星才显著贡献干扰功率在这个波束覆盖的天空区域内卫星是否出现、出现几颗本质上接近一个随机事件。用PPP就可以直接写出干扰累积分布对应的强度测度而不需要明确每颗干扰卫星的编号和位置。仿真时生成干扰卫星位置的关键代码def generate_interfering_satellites(lambda_int, area_visible, n_realize): 在可见球冠区域内生成PPP干扰源位置 lambda_int: 干扰卫星密度每平方米 area_visible: 可见球冠面积平方米 n_realize: 需要生成的随机实现数 返回每个实现中干扰源数量列表与位置列表 import numpy as np reals [] for _ in range(n_realize): n np.random.poisson(lambda_int * area_visible) # 泊松抽样 # 在球冠区域内均匀抽样简化用均匀球面角 cos_theta 1 - np.random.rand(n) * (1 - np.cos(np.radians(30))) # 30°半锥角 theta np.arccos(cos_theta) phi 2 * np.pi * np.random.rand(n) reals.append((n, theta, phi)) return reals逻辑说明泊松抽样决定了每次实现中干扰源的数量期望值等于λ_int乘以可见面积。球冠内的均匀抽样采用逆变换法先对cosθ均匀抽样再对方位角均匀抽样这比直接在θ上均匀抽样更能反映球面面积的分布。30°半锥角对应的是用户天线主波束宽度代表只有这个锥角内的干扰卫星会被接收机看到。4.2 蒙特卡洛SINR仿真框架有了服务卫星和干扰卫星的几何SINR定义如下SINR P_s · G_s(θ_s) / [Σ_i P_i · G_i(θ_i) N_0 · B]P_s和P_i分别为服务与干扰卫星的发射功率G_s和G_i为卫星天线在用户方向上的增益θ是离轴角。这里假设所有卫星发射功率相同采用全频复用即所有可见卫星都是潜在干扰源。一个简化的单快照SINR计算函数如下直接整合了前两章的损耗计算def compute_sinr(user_lla, sat_position, sat_velocity, freq_hz, p_tx_dbm, g_tx_db, interferer_positions, bandwidth_hz): 计算下行链路SINR user_lla: 用户位置 (lat, lon, alt) sat_position: 服务卫星位置 interferer_positions: 干扰卫星位置数组 (N_interf, 3) # 服务链路距离和路径损耗 d_s np.linalg.norm(sat_position - user_ecef(user_lla)) pl_s free_space_path_loss(d_s, freq_hz) # 接收信号功率 p_rx p_tx_dbm g_tx_db - pl_s - atm_loss(d_s) # 集总干扰功率 p_interf 0.0 for pos in interferer_positions: d_i np.linalg.norm(pos - user_ecef(user_lla)) pl_i free_space_path_loss(d_i, freq_hz) p_interf 10**((p_tx_dbm g_tx_db - pl_i) / 10) # 噪声功率 k 1.38e-23 noise_w k * 290 * bandwidth_hz p_noise_dbm 10 * np.log10(noise_w) 30 sinr_linear 10**(p_rx / 10) / (p_interf 10**(p_noise_dbm / 10)) return 10 * np.log10(sinr_linear)逻辑说明先算服务链路的接收功率再遍历所有干扰卫星把每颗干扰卫星的接收功率线性累加最后与噪声功率一起求SINR。这里要注意干扰功率必须在线性域累加不能用dB值直接相加这是最容易犯错的地方。参数说明天线增益G_t在简化模型里被当作常数但实际卫星相控阵天线对用户方向的增益随离轴角变化如果仿真要求精度高应该把天线方向图函数加进来——通常用cos^q(θ)模型q值决定了波束尖锐程度。干扰卫星的发射功率在频分多址或极化复用场景下应该乘以复用因子否则会过大估计干扰。4.3 覆盖概率与中断概率统计分析蒙特卡洛仿真的输出是一组SINR采样值覆盖概率定义为SINR超过阈值γ的概率即P_cov(γ) P(SINR γ)。这个指标直接对应系统设计中的应用场景γ0dB时覆盖概率高低反映基本通信能力γ10dB时反映高质量链路的可用性。def monte_carlo_coverage(n_iter10000, sinr_threshold_db0): 蒙特卡洛仿真覆盖概率 返回覆盖概率、平均SINR、SINR分布 sinr_samples [] for _ in range(n_iter): # 生成服务卫星取可见区域内距离用户最近的一颗 sats generate_walker_constellation(T288, P18) elevs [compute_elevation_angle(USER_POS, s) for s in sats] visible [s for s, e in zip(sats, elevs) if e 10] if not visible: sinr_samples.append(-np.inf) continue # 选距离最近通常仰角最高的卫星作为服务星 serv min(visible, keylambda s: np.linalg.norm(s - USER_ECEF)) # 生成干扰PPP area_vis visible_area() interf generate_interfering_satellites(LAMBDA_INT, area_vis, 1) sinr compute_sinr(USER_POS, serv, np.zeros(3), 20e9, 40, 30, interf, 500e6) sinr_samples.append(sinr) cov_prob np.mean(np.array(sinr_samples) sinr_threshold_db) return cov_prob, np.mean(sinr_samples), np.percentile(sinr_samples, [10, 50, 90])逻辑说明每次迭代先做一次星座遍历找到仰角大于10°的可见卫星集合选距离用户最近的一颗作为服务卫星随后在可见区域内生成PPP干扰源计算SINR。无效迭代没有可见卫星计数为无穷大的中断计入中断概率。最终覆盖概率就是SINR超过阈值的样本占比。参数说明n_iter取10000时统计波动在±1%以内如果仿真速度慢可以降到2000但覆盖概率的置信区间会变宽。阈值γ建议同时输出0dB、5dB、10dB三个值对应不同业务需求。这里用T288颗完整Walker星座逐颗遍历验证PPP模型得到的结果与论文参考值对比误差通常在1dB以内。5. 仿真加速与参数验证的三个技巧5.1 预计算可见卫星表代替逐次遍历蒙特卡洛每轮迭代都重新生成整个Walker星座再遍历一遍数据上浪费很大。常见做法是仿真前先预计算一次用户位置对应的可见卫星集合和距离、仰角然后只从集合里抽样服务卫星和干扰卫星。对于静止终端用户到卫星的几何关系只随时间变化因此在快照式仿真中预计算可以节省接近一半的运行时间。# 预计算可见卫星只做一次 all_sats generate_walker_constellation(...) visible_sats [s for s in all_sats if compute_elevation_angle(user_pos, s) 10] # 蒙特卡洛主循环里直接用预计算结果 for i in range(n_iter): serv np.random.choice(visible_sats.shape[0]) # 干扰源从剩余可见卫星中抽取密度按频段复用比例5.2 三个必做敏感性参数仿真做完别直接写结论先跑三个参数敏感性实验最小仰角、频率复用系数和星座相位因子。低轨星座全频复用场景下干扰受限是常态频率复用系数从1降到0.5覆盖概率可能提升20%以上相位因子F从1调到2部分纬度带上覆盖概率会抖动2-3dB。这些规律在确定性Walker模型中很难直接看出来PPP模型通过密度参数的变化让因果链条变得清晰。5.3 用解析表达式交叉验证仿真器PPP模型有一个可解析的覆盖概率近似公式P_cov ≈ exp(-λ · C · γ^(-2/α))其中C是取决于路径损耗指数和星座几何的常数α是路径损耗指数。这个公式虽然粗糙但可以直接用仿真参数代入把结果跟你写出的蒙特卡洛结果对比。如果两者差距超过20%先检查代码里干扰累加是不是用了对数域、可见区域面积的计算是不是用了球面近似这些是低轨卫星干扰仿真里最常被忽略的两个误差来源。本文还有配套的精品资源点击获取
简介基于随机几何的低轨星座下行链路仿真与分析完整资料面向具备Python编程基础的卫星通信研究人员与系统工程师重点解决低轨星座建模、星地链路损耗计算及多星干扰评估问题。其中以二项点过程(BPP)构建星座模型深入分析单星与多星场景下的路径损耗和同频干扰提出干扰期望计算方法可用于巨型低轨星座性能分析与星地链路设计参考。配套代码覆盖参数设置、卫星与地面站随机生成、自由空间及大气损耗、接收功率与噪声功率计算、SINR分析、结果可视化等完整环节代码可运行且附带逐段解释方便读者对照学习。压缩包内为单个docx文档共1个文件大小约51KB轻量便携目前已有105人学习/下载适合作为低轨卫星通信仿真的入门与实践参考。内容还展望了天地一体化网络、大规模星座干扰建模等未来方向为后续研究提供延伸切入点。1. 为什么低轨星座下行链路仿真要改用随机几何低轨星座规模越做越大几十上百颗卫星只是起步千星级别的系统已经开始规划。传统仿真做法信赖星历文件加轨道积分器逐颗卫星计算星下点、仰角、路径损耗精度确实高但代价是算力需求随卫星数量膨胀——分析一条下行链路要遍历全部可见卫星换一个用户位置就得重算一轮。随机几何的思路完全不同它把卫星在观察时刻的空间分布当作随机点过程的一次实现用强度函数近似星座的整体几何直接推导SINR、覆盖概率、中断概率的表达式或半解析结果。这套方法牺牲掉单颗卫星的精确轨迹换来参数尺度上的可解释性轨道高度、卫星数量、发射功率、频率复用因子如何影响覆盖率几分钟内就能跑完。本文就用这套思路从星座建模开始逐步实现完整的低轨卫星下行链路仿真。2. 低轨星座建模从Walker构型到随机点过程的映射2.1 Walker星座参数与初始轨道生成低轨星座最常见的构型是Walker Delta星座四个参数就能描述T是卫星总数P是轨道面数F是相位因子i是轨道倾角。例如常见的参考方案T288、P18、F1、i53°每面16颗卫星近圆形轨道高度550km。相位因子F控制相邻轨道面间卫星的相位偏移量大小为2πF/T这个参数直接影响面间相邻卫星的最小间距进而决定覆盖均匀度和同频干扰的分布特性。生成星座初始ECI位置的最小Python实现如下import numpy as np def generate_walker_constellation(T288, P18, F1, i_deg53.0, h_km550.0): 生成Walker Delta星座的ECI位置向量圆轨道假设 T: 卫星总数 P: 轨道面数 F: 相位因子 i_deg: 轨道倾角度 h_km: 轨道高度km 返回卫星ECI坐标数组形状 (T, 3) R_E 6371.0e3 # 地球平均半径m r h_km * 1e3 R_E # 轨道半径m S T // P # 每面卫星数 sats np.zeros((T, 3)) idx 0 for p in range(P): raan 2 * np.pi * p / P # 升交点赤经均分 for s in range(S): # 面内相角等间隔 相位因子偏移 nu 2 * np.pi * s / S 2 * np.pi * F * p / T # 轨道平面内的位置圆轨道 x_orb r * np.cos(nu) y_orb r * np.sin(nu) # 1) 绕X轴旋转倾角 cos_i, sin_i np.cos(np.radians(i_deg)), np.sin(np.radians(i_deg)) x1, y1, z1 x_orb, y_orb * cos_i, y_orb * sin_i # 2) 绕Z轴旋转升交点赤经 cos_r, sin_r np.cos(raan), np.sin(raan) sats[idx, 0] x1 * cos_r - y1 * sin_r sats[idx, 1] x1 * sin_r y1 * cos_r sats[idx, 2] z1 idx 1 return sats逻辑说明先按轨道面遍历每个面的升交点赤经RAAN在360°内等间隔分布再在面内按相角等间隔放置S颗卫星并叠加相位因子偏移。圆轨道假设下真近点角就等于相角省去了开普勒方程求解。旋转矩阵分两步执行先绕X轴旋转轨道倾角把轨道平面从赤道平面倾斜到指定角度再绕Z轴旋转RAAN把升交点对准正确方向。参数说明轨道高度h_km决定轨道周期和覆盖半径550km对应的轨道周期约95分钟覆盖圆半径约2400km按最小仰角10°计算。T和P的比例决定了面内卫星间距P越多单面卫星越少面间相对运动越复杂。F取0到P-1之间的整数工程上倾向选F1或F2使相邻轨道面的卫星相位错开避免形成规则的干扰纹理。这里的地球半径用平均半径6371km如果要做特定真实星座的精确复现需要换成WGS84椭球并计入J2摄动引起的RAAN长期漂移。提示Walker参数T288、P18、F1适用于系统级覆盖预研。如果你要复现某个具体星座的论文结果先确认论文用的地面覆盖最小仰角这个值对可见卫星数量影响极大。2.2 从离散卫星到随机点过程强度函数的建模逻辑有了离散卫星的位置为什么还要引入随机几何核心原因是计算复杂度的量级差异。假设T288颗卫星地面用户在任意时刻平均可见约14颗卫星550km高度、10°仰角约束下遍历计算完全可以接受。但当星座规模到10000颗时平均可见卫星数接近500颗逐颗精确计算距离、天线增益和仰角的时间开销随T线性增长而蒙特卡洛仿真又需要成千上万次重复整体开销迅速变得难以承受。随机几何的做法是用泊松点过程PPP近似卫星的空间分布把问题从「枚举每颗卫星」转成「对空间积分」。核心思想是用户可见的球冠区域内卫星位置被建模为独立均匀分布的随机点密度λ T / (4πR²)R是轨道半径。可见球冠面积由轨道高度和最小仰角决定平均可见卫星数为λ乘以球冠面积。这个近似在卫星数量大、最小仰角高、轨道高度低的时候最准因为可见区域收缩后边界效应减弱区域内分布更接近均匀独立。下表给出了不同参数组合下的近似适用性判断直接决定你后边仿真用确定性方法还是概率方法参数典型范围对PPP近似的影响仿真建议卫星总数T100~10000T越大泊松近似越准T300可逐颗精确计算最小仰角10°~40°仰角越高可见区域越小近似越可靠10°~20°适合PPP验证轨道高度340~1200km高度越低覆盖半径越小边界效应越弱340km处PPP误差最小轨道面数P12~24P越少轨道纹理越明显偏离PPP越大P12建议用确定性模型关键是要意识到PPP不是对真实星座的精确描述而是一种把「卫星在可见区域里以某密度随机出现」当作等效模型的工程近似。低轨星座的轨道面如果只有几个卫星分布会呈现明显的带状纹理此时PPP会低估干扰的空间聚集性。一般做法是先在确定性Walker星座上验证一次再切到PPP模型做参数扫描两种方法互为校验。2.3 可见性判断与用户-卫星几何计算无论用哪种模型仿真里都要计算用户到一颗卫星的仰角。定义用户位置为经纬高卫星位置为ECI坐标需要把卫星位置转到用户本地东北天ENU坐标系才能算出仰角。下面的函数实现了这个几何过程def compute_elevation_angle(user_lla, sat_eci): 计算用户到卫星的仰角 user_lla: (纬度deg, 经度deg, 高度m) sat_eci: 卫星ECI坐标 (x, y, z) m 返回仰角度 lat np.radians(user_lla[0]) lon np.radians(user_lla[1]) alt user_lla[2] # WGS84椭球参数 a 6378137.0 # 长半轴m e2 6.69437999014e-3 # 偏心率平方 N a / np.sqrt(1 - e2 * np.sin(lat)**2) # 用户ECEF坐标 ux (N alt) * np.cos(lat) * np.cos(lon) uy (N alt) * np.cos(lat) * np.sin(lon) uz (N * (1 - e2) alt) * np.sin(lat) # 卫星到用户的向量 dx sat_eci[0] - ux dy sat_eci[1] - uy dz sat_eci[2] - uz # ECEF到ENU旋转矩阵 sin_lat, cos_lat np.sin(lat), np.cos(lat) sin_lon, cos_lon np.sin(lon), np.cos(lon) e -sin_lon * dx cos_lon * dy n -sin_lat * cos_lon * dx - sin_lat * sin_lon * dy cos_lat * dz u cos_lat * cos_lon * dx cos_lat * sin_lon * dy sin_lat * dz return np.degrees(np.arctan2(u, np.sqrt(e**2 n**2)))逻辑说明先把用户经纬高通过WGS84椭球转成ECEF直角坐标然后用ECEF到ENU的方向余弦矩阵把卫星相对用户的向量分解到东、北、天三个方向。仰角是天向分量与水平分量的反正切水平分量由东向和北向的平方根给出。参数说明最小仰角阈值通常取10°或25°对应不同的覆盖概率要求。10°仰角可以看到更远的低轨卫星但信号穿过的对流层路径更长大气衰减和闪烁更严重25°仰角服务质量更稳定但可见卫星数量减少约30%。在干扰分析里这个选择直接影响干扰源的密度建议在参数表中同时输出两个仰角阈值下的结果。3. 下行链路损耗计算从自由空间到大气衰减3.1 自由空间路径损耗与多普勒频移低轨卫星下行链路损耗中自由空间路径损耗占绝对主导。公式为PL 20log10(4πd/λ)其中d是用户到卫星的直线距离。对20GHz载频、550km轨道高度、仰角10°时的最大斜距约2400km来算路径损耗约为186dB左右仰角90°卫星在头顶时斜距等于轨道高度550km路径损耗约173dB。也就是说仅仰角变化就能引入13dB的链路预算波动。多普勒频移方面低轨卫星轨道速度约7.6km/s径向相对速度在地面用户看来的极端情况可达7km/s量级。20GHz承载频率下最大多普勒频移可接近500kHz远超地面移动通信系统通常面临的规模。仿真相干时间内如果只取瞬时快照多普勒不会影响SINR但若仿真时间窗口拉长必须计算每颗卫星相对用户的径向速度并折算到频偏否则接收机频率补偿误差会变成实际损耗。def free_space_path_loss(d_m, freq_hz): 自由空间路径损耗d_m为斜距米freq_hz为载波频率 c 2.99792458e8 return 20 * np.log10(4 * np.pi * d_m * freq_hz / c) def doppler_shift(sat_vel, user_pos, sat_pos, freq_hz): 计算多普勒频移sat_vel为卫星速度矢量(m/s) c 2.99792458e8 r_vec sat_pos - user_pos r_dist np.linalg.norm(r_vec) radial_v np.dot(sat_vel, r_vec) / r_dist return radial_v * freq_hz / c逻辑说明多普勒频移的代码里关键一步是把卫星速度矢量投影到用户-卫星连线上得到径向速度。径向速度为正表示卫星正在靠近用户频移为正为负表示远离频移为负。LEO卫星过顶过程中多普勒由正到负的过零变化是常见特征。参数说明载波频率越高多普勒越显著这也是Ka频段比L频段更难做低轨通信的原因之一。如果你仿真的是20GHz下行建议把子载波间隔和频偏补偿算法纳入链路级仿真但在系统级信噪比分析中只需记录最近服务卫星的多普勒值用于判断是否需要引入频偏校正开销。3.2 大气衰减与雨衰模型自由空间损耗之上电波穿过对流层时还有气体吸收和雨衰。低频段L/S的大气衰减可以忽略到了Ka/Ku频段氧气和水蒸气的分子吸收不可再忽略。ITU-R P.676给出了天顶方向大气衰减的标准计算方法。在20GHz、天顶方向氧气衰减约0.2dB水蒸气约0.1dB但仰角降到10°时电波在大气层中的路径长度增加约5倍衰减合计可超过1.5dB。雨衰是另一个大头用ITU-R P.618推荐的模型估算。核心输入是区域降雨强度R_0.01每年超过0.01%时间的降雨率单位为mm/h以及仰角、极化倾角和频率。一个粗略的参考值在K频段中纬度地区0.01%时间降雨强度约30mm/h对应20GHz、30°仰角下雨衰约3-6dB。频段频率范围天顶大气衰减0.1%时间雨衰参考值适用场景L1-2GHz可忽略几乎可忽略星间链路、物联网S2-4GHz可忽略0.5dB移动通信备份Ku12-18GHz0.3-0.8dB1-3dB传统广播Ka20-30GHz0.5-1.5dB3-10dB宽带低轨星座仿真里处理雨衰有一个常见选择直接用概率分布注入而不是精确计算每次降雨的时空分布。最简单的方式是假设雨衰服从对数正态分布在每次蒙特卡洛迭代里独立抽样叠加。3.3 接收功率与链路预算表把上述损耗项组合起来就得到下行链路预算。接收到的载波功率可以写成P_r P_t G_t - L_fs - L_atm - L_rain G_r其中P_t为卫星发射功率G_t为卫星天线增益L_fs为自由空间损耗L_atm为大气气体衰减L_rain为雨衰G_r为用户终端天线增益。以下是一组典型参数下的链路预算参数数值计算说明卫星发射功率10W40dBm下行通常比上行功率充裕卫星天线增益30dBi多波束相控阵单波束增益EIRP70dBm40 30自由空间损耗173dB550km垂直高度、20GHz大气衰减0.8dB仰角30°、中纬度雨衰0.1%时间2.5dB20GHz中纬度中等降雨用户天线增益25dBi固定终端接收功率-81.3dBm70 - 173 - 0.8 - 2.5 25噪声功率-104dBm500MHz带宽噪声系数1.5dBSNR22.7dB接收功率减噪声功率这个表里SNR的余量很大因为系统级干扰分析的主要矛盾在于同频干扰而不是噪声。低轨星座通常全频复用服务卫星以外还有成百上千颗卫星落在用户天线波束内它们造成的集总干扰要比底噪高一个量级。所以下一章要解决的核心问题是在那么多同频卫星里SINR还剩下多少。4. 干扰分析与下行链路SINR仿真4.1 同频干扰的随机几何建模考虑一个下行场景一颗服务卫星向目标用户发送信号其他同频卫星的信号同时到达用户天线形成集总干扰。在传统确定性仿真里你需要得到每一颗干扰卫星的精确位置、天线增益、发射功率然后枚举求和。随机几何的建模则把干扰卫星集合视为可见球冠区域内的PPP密度为λ_int α·λ其中α是频率复用系数全复用1。为什么PPP在这里是合理的近似低轨星座干扰源数量大、距离散布广、各干扰链路间的相互依赖较弱。用户天线通常有指向性只有落在主波束范围内的干扰卫星才显著贡献干扰功率在这个波束覆盖的天空区域内卫星是否出现、出现几颗本质上接近一个随机事件。用PPP就可以直接写出干扰累积分布对应的强度测度而不需要明确每颗干扰卫星的编号和位置。仿真时生成干扰卫星位置的关键代码def generate_interfering_satellites(lambda_int, area_visible, n_realize): 在可见球冠区域内生成PPP干扰源位置 lambda_int: 干扰卫星密度每平方米 area_visible: 可见球冠面积平方米 n_realize: 需要生成的随机实现数 返回每个实现中干扰源数量列表与位置列表 import numpy as np reals [] for _ in range(n_realize): n np.random.poisson(lambda_int * area_visible) # 泊松抽样 # 在球冠区域内均匀抽样简化用均匀球面角 cos_theta 1 - np.random.rand(n) * (1 - np.cos(np.radians(30))) # 30°半锥角 theta np.arccos(cos_theta) phi 2 * np.pi * np.random.rand(n) reals.append((n, theta, phi)) return reals逻辑说明泊松抽样决定了每次实现中干扰源的数量期望值等于λ_int乘以可见面积。球冠内的均匀抽样采用逆变换法先对cosθ均匀抽样再对方位角均匀抽样这比直接在θ上均匀抽样更能反映球面面积的分布。30°半锥角对应的是用户天线主波束宽度代表只有这个锥角内的干扰卫星会被接收机看到。4.2 蒙特卡洛SINR仿真框架有了服务卫星和干扰卫星的几何SINR定义如下SINR P_s · G_s(θ_s) / [Σ_i P_i · G_i(θ_i) N_0 · B]P_s和P_i分别为服务与干扰卫星的发射功率G_s和G_i为卫星天线在用户方向上的增益θ是离轴角。这里假设所有卫星发射功率相同采用全频复用即所有可见卫星都是潜在干扰源。一个简化的单快照SINR计算函数如下直接整合了前两章的损耗计算def compute_sinr(user_lla, sat_position, sat_velocity, freq_hz, p_tx_dbm, g_tx_db, interferer_positions, bandwidth_hz): 计算下行链路SINR user_lla: 用户位置 (lat, lon, alt) sat_position: 服务卫星位置 interferer_positions: 干扰卫星位置数组 (N_interf, 3) # 服务链路距离和路径损耗 d_s np.linalg.norm(sat_position - user_ecef(user_lla)) pl_s free_space_path_loss(d_s, freq_hz) # 接收信号功率 p_rx p_tx_dbm g_tx_db - pl_s - atm_loss(d_s) # 集总干扰功率 p_interf 0.0 for pos in interferer_positions: d_i np.linalg.norm(pos - user_ecef(user_lla)) pl_i free_space_path_loss(d_i, freq_hz) p_interf 10**((p_tx_dbm g_tx_db - pl_i) / 10) # 噪声功率 k 1.38e-23 noise_w k * 290 * bandwidth_hz p_noise_dbm 10 * np.log10(noise_w) 30 sinr_linear 10**(p_rx / 10) / (p_interf 10**(p_noise_dbm / 10)) return 10 * np.log10(sinr_linear)逻辑说明先算服务链路的接收功率再遍历所有干扰卫星把每颗干扰卫星的接收功率线性累加最后与噪声功率一起求SINR。这里要注意干扰功率必须在线性域累加不能用dB值直接相加这是最容易犯错的地方。参数说明天线增益G_t在简化模型里被当作常数但实际卫星相控阵天线对用户方向的增益随离轴角变化如果仿真要求精度高应该把天线方向图函数加进来——通常用cos^q(θ)模型q值决定了波束尖锐程度。干扰卫星的发射功率在频分多址或极化复用场景下应该乘以复用因子否则会过大估计干扰。4.3 覆盖概率与中断概率统计分析蒙特卡洛仿真的输出是一组SINR采样值覆盖概率定义为SINR超过阈值γ的概率即P_cov(γ) P(SINR γ)。这个指标直接对应系统设计中的应用场景γ0dB时覆盖概率高低反映基本通信能力γ10dB时反映高质量链路的可用性。def monte_carlo_coverage(n_iter10000, sinr_threshold_db0): 蒙特卡洛仿真覆盖概率 返回覆盖概率、平均SINR、SINR分布 sinr_samples [] for _ in range(n_iter): # 生成服务卫星取可见区域内距离用户最近的一颗 sats generate_walker_constellation(T288, P18) elevs [compute_elevation_angle(USER_POS, s) for s in sats] visible [s for s, e in zip(sats, elevs) if e 10] if not visible: sinr_samples.append(-np.inf) continue # 选距离最近通常仰角最高的卫星作为服务星 serv min(visible, keylambda s: np.linalg.norm(s - USER_ECEF)) # 生成干扰PPP area_vis visible_area() interf generate_interfering_satellites(LAMBDA_INT, area_vis, 1) sinr compute_sinr(USER_POS, serv, np.zeros(3), 20e9, 40, 30, interf, 500e6) sinr_samples.append(sinr) cov_prob np.mean(np.array(sinr_samples) sinr_threshold_db) return cov_prob, np.mean(sinr_samples), np.percentile(sinr_samples, [10, 50, 90])逻辑说明每次迭代先做一次星座遍历找到仰角大于10°的可见卫星集合选距离用户最近的一颗作为服务卫星随后在可见区域内生成PPP干扰源计算SINR。无效迭代没有可见卫星计数为无穷大的中断计入中断概率。最终覆盖概率就是SINR超过阈值的样本占比。参数说明n_iter取10000时统计波动在±1%以内如果仿真速度慢可以降到2000但覆盖概率的置信区间会变宽。阈值γ建议同时输出0dB、5dB、10dB三个值对应不同业务需求。这里用T288颗完整Walker星座逐颗遍历验证PPP模型得到的结果与论文参考值对比误差通常在1dB以内。5. 仿真加速与参数验证的三个技巧5.1 预计算可见卫星表代替逐次遍历蒙特卡洛每轮迭代都重新生成整个Walker星座再遍历一遍数据上浪费很大。常见做法是仿真前先预计算一次用户位置对应的可见卫星集合和距离、仰角然后只从集合里抽样服务卫星和干扰卫星。对于静止终端用户到卫星的几何关系只随时间变化因此在快照式仿真中预计算可以节省接近一半的运行时间。# 预计算可见卫星只做一次 all_sats generate_walker_constellation(...) visible_sats [s for s in all_sats if compute_elevation_angle(user_pos, s) 10] # 蒙特卡洛主循环里直接用预计算结果 for i in range(n_iter): serv np.random.choice(visible_sats.shape[0]) # 干扰源从剩余可见卫星中抽取密度按频段复用比例5.2 三个必做敏感性参数仿真做完别直接写结论先跑三个参数敏感性实验最小仰角、频率复用系数和星座相位因子。低轨星座全频复用场景下干扰受限是常态频率复用系数从1降到0.5覆盖概率可能提升20%以上相位因子F从1调到2部分纬度带上覆盖概率会抖动2-3dB。这些规律在确定性Walker模型中很难直接看出来PPP模型通过密度参数的变化让因果链条变得清晰。5.3 用解析表达式交叉验证仿真器PPP模型有一个可解析的覆盖概率近似公式P_cov ≈ exp(-λ · C · γ^(-2/α))其中C是取决于路径损耗指数和星座几何的常数α是路径损耗指数。这个公式虽然粗糙但可以直接用仿真参数代入把结果跟你写出的蒙特卡洛结果对比。如果两者差距超过20%先检查代码里干扰累加是不是用了对数域、可见区域面积的计算是不是用了球面近似这些是低轨卫星干扰仿真里最常被忽略的两个误差来源。本文还有配套的精品资源点击获取
