简介面向卫星通信研究者和工程师这份资源复现了基于随机几何的低轨星座下行链路仿真与分析工作。内容围绕低轨卫星网络移动性强、电磁环境复杂的特点基于二项点过程BPP构建星座模型覆盖单星与多星场景的路径损耗、接收功率、噪声及干扰计算并通过干扰期望计算方法评估星地链路性能。包体为1个docx文档压缩后约51KB文档内含详细理论推导和可直接运行的Python代码涵盖参数设置、卫星与地面站生成、SINR计算、干扰分析与可视化配有逐步解释方便对照学习。目前已有105人学习。借助该文档可理解随机几何在星座建模中的应用掌握星座参数对通信性能的影响为巨型低轨星座的干扰建模和链路设计提供可扩展的仿真参考。1. 低轨卫星通信的随机几何仿真先回答“为什么用PPP做下行链路”低轨卫星通信的下行链路仿真最花时间的往往不是轨道积分而是干扰说不清。真实低轨星座有几十个轨道面每颗卫星都在移动用户天线还要跟着波束切换如果只用精确星历逐秒跑结果只是某个具体时刻、某个纬度的特解。随机几何把卫星位置建模成球面上的泊松点过程PPP用密度和高度代替轨道参数覆盖概率和信干噪比分布这类指标第一次可以从统计上直接回答。这套方法适合正在复现低轨卫星通信论文、准备毕业设计或者做频率复用预研的工程师。下面按“星座建模 → 损耗计算 → 干扰分析 → 避坑 → 验证”的顺序给出可以直接改参数跑的 Python 代码和工程坑。2. 星座建模与可见性判断用球面泊松点过程生成低轨星座的代码实现2.1 为什么低轨星座适合用随机几何的PPP近似低轨卫星通信系统的轨道高度通常在 500km 到 1200km 之间。与静止轨道不同每颗低轨卫星相对地面用户移动很快几分钟内会从地平线附近升到头顶再降到另一侧地平线。对一个固定位置的地面用户来说卫星在哪个方向、距离多远是时刻变化的星座设计者更关心“任意时刻有至少一颗卫星可见”的概率而不是某颗卫星某秒的精确坐标。随机几何正是为这类统计问题准备的。常见的低轨星座建模路线有两种。确定性轨道建模需要输入轨道面、相位、TLE 或者开普勒根数适合做星间切换、波束指向的事件级仿真但一次只能代表一套具体星座构型。随机几何建模忽略卫星编号只关心单位面积上期望有多少颗卫星、在统计上服从什么分布输出的是覆盖概率、中断概率这类系统级指标。两者不矛盾确定性建模用来校核特定方案随机几何用来回答“如果星座密度提升一倍覆盖概率怎么变”这类趋势问题。实际仿真中使用 PPP 还有一个物理上的合理性。真实低轨星座虽然结构确定但经过一段时间后对地面用户而言卫星到达方向会随轨道运动和用户位置变化形成近似随机的模式。PPP 的独立性让干扰分析中的聚合干扰拉普拉斯变换可以分解成单卫星的乘积这是随机几何理论推导的基础。需要精确匹配卫星总数时可以用 BPP 固定每轮快照的卫星数这也被大量论文当作 Monte Carlo 验证方式。2.2 在球面上生成PPP卫星的代码先给出核心代码。注意这里的“卫星”指某时刻能与用户发生链路的可见卫星不是轨道上的全部卫星。我们用一个固定高度球面来近似所有卫星位置。import numpy as np def simulate_leo_sat_positions(mean_sat_num, height_km): 在球面上生成泊松点过程的卫星位置。 参数 - mean_sat_num : floatPPP的期望卫星数即密度乘以球面面积 - height_km : float轨道高度km 返回 - sat_pos : ndarray, shape (N,3)地心直角坐标单位 m R_earth 6371.0e3 # 地球半径m R_orbit R_earth height_km * 1e3 # 泊松点过程卫星数是随机变量先抽一个泊松数 N np.random.poisson(mean_sat_num) # 均匀抽样球面上的点。注意经度可以直接均匀极角必须用反余弦 lon 2.0 * np.pi * np.random.rand(N) cos_theta 1.0 - 2.0 * np.random.rand(N) # theta ∈ [0, pi] 的余弦值均匀在[-1,1] sin_theta np.sqrt(1.0 - cos_theta ** 2) # 球坐标 (r, theta, lon) 转地心直角坐标 x R_orbit * sin_theta * np.cos(lon) y R_orbit * sin_theta * np.sin(lon) z R_orbit * cos_theta return np.column_stack([x, y, z])这里最关键的是球面均匀抽样。很多人第一次写会直接让纬度和经度均匀分布结果高纬度地区卫星密度偏大造成干扰分析失真用cos_theta 1 - 2*rand等效于先均匀抽样单位球面切向量再做反变换能保证在每个相等面积上落入相同期望点数。mean_sat_num就是 PPP 的强度参数 λ 乘上球面面积 (4\pi R_{orbit}^2)。如果你要匹配一个具体星座比如“轨道面 15 个每个轨道面 40 颗”那期望卫星数就是 600但 PPP 的瞬时卫星数会随机波动这正是随机几何与确定性建模的差别。如果你想让卫星总数固定就把上面生成 N 的那一行换成N int(mean_sat_num)其余不用动这就变成了 BPP。做论文对比时我建议两边都跑用 PPP 做理论推导用 BPP 做蒙特卡洛验证两者差异就是随机几何近似的代价。2.3 地面用户到卫星的距离、仰角与可见性判断有了卫星位置下一步是判断地面用户能看到哪些卫星。这里需要把用户经纬度也变成地心直角坐标然后计算用户到卫星的距离矢量和仰角。def llh_to_ecef(lat_deg, lon_deg, alt_m): 把大地坐标转地心直角坐标简单球体模型 R 6371.0e3 alt_m lat np.deg2rad(lat_deg) lon np.deg2rad(lon_deg) cos_lat np.cos(lat) return np.array([ R * cos_lat * np.cos(lon), R * cos_lat * np.sin(lon), R * np.sin(lat) ]) def compute_distance_elevation(user_ecef, sat_pos): 计算每个卫星到用户的距离和仰角。 返回 distance_m, elev_deg d_vec sat_pos - user_ecef # 从用户指向卫星 dist_m np.linalg.norm(d_vec, axis1) # 用户位置指向地心的单位向量也就是用户本地天顶方向 zenith user_ecef / np.linalg.norm(user_ecef) # 仰角 sin 值 距离向量在天顶方向上的投影 / 距离 sin_elev np.sum(d_vec * zenith, axis1) / dist_m sin_elev np.clip(sin_elev, -1.0, 1.0) elev_deg np.degrees(np.arcsin(sin_elev)) return dist_m, elev_deg这段代码的几何关系可以这样理解用户站在球面上本地“竖直向上”的方向就是用户位置向量向外延长的那根线即user_ecef的单位向量。从用户指向卫星的距离矢量和竖直方向做点积得到的是距离在竖直方向的投影投影除以距离就是仰角的正弦。卫星在正头顶时仰角 90°在地平线附近仰角接近 0°跑在另一边就会算出负值应当丢弃。有了距离和仰角就可以根据最低仰角做筛选。工程上传统卫星通信通常要求最低通信仰角在 8° 到 25° 之间低轨卫星为了增加覆盖可能放宽到 5°但低于 5° 的链路大气损耗大、遮挡多不建议计入统计。筛选逻辑很简单def get_visible_satellites(user_ecef, sat_pos, min_elev_deg10.0): dist_m, elev_deg compute_distance_elevation(user_ecef, sat_pos) vis elev_deg min_elev_deg return dist_m[vis], elev_deg[vis]筛选完剩余卫星数通常在几颗到几十颗。如果平均剩下不到一颗说明星座密度太低或最低仰角设得太高结果会以“脱网”场景为主覆盖概率整体偏低。这个现象在第 5 章的坑里再展开。2.4 星座建模的边界PPP的前提与参数标定用 PPP 做低轨星座建模有一个物理前提需要交代清楚PPP 要求卫星位置在球面上均匀且独立。真实星座由有限个轨道面构成存在明显的结构相关性某些地区在某段时间可能“看见”多颗星另一些地区则稀疏。如果把真实星座的瞬时快照跑几千个时刻平均下来覆盖概率大概率落在 PPP 仿真结果的置信区间内但如果论文要比较“轨道面数对性能的影响”PPP 就表达不了轨道结构。因此这里有两个常见标定方法。一是把 PPP 的期望卫星数设为星座卫星总数用于评估系统容量上限二是把 PPP 期望卫星数设为“某个时刻用户可见或可干扰到的同频卫星数”用于评估干扰受限场景。前者适合覆盖分析后者适合干扰分析。实际项目中我更倾向于按轨道高度和天线波束宽度先算一个“干扰可见圈”再在圈内用 PPP 建模这样能减少全局 PPP 带来的过度保守干扰。3. 下行链路损耗计算自由空间、大气衰减与天线增益怎么进链路预算3.1 先算自由空间路径损耗斜距不是高度别拿轨道高度直接填低轨卫星与地面用户的距离随仰角变化非常大。以 550km 轨道高度为例卫星过顶时距离约 550km但仰角 10° 时斜距大约 2100km距离相差近 4 倍自由空间损耗相差约 12dB。很多初版仿真直接把h当距离填进弗里斯公式结果覆盖概率虚高干扰分析也失去意义。正确的做法是使用前面算出的dist_m先计算自由空间路径损耗[ FSPL(dB) 20\log_{10}(d) 20\log_{10}(f) - 147.55 ]其中d单位是米f单位是 Hz147.55来自 ( (4\pi/c)^2 ) 的对数变换。注意这个式子里的常数如果你用 km 和 GHz常数会变成 32.45单位搞错是最常见的低级错误。下面这段代码把距离、频率、仰角一并处理并保留逐链路的损耗供后续干扰求和def free_space_path_loss(distance_m, freq_hz): 自由空间路径损耗单位 dB return 20.0 * np.log10(distance_m) 20.0 * np.log10(freq_hz) - 147.553.2 大气衰减、雨衰和天线增益从单值补偿到按仰角插值LEO 下行链路的工作频率从 L 波段到 Ka 波段都有。频率越高大气衰减越明显。工程上最粗略的做法是在链路预算里加一个固定大气损耗余量比如 L 波段加 0.5dB、Ka 波段加 2dB但这是做系统预算用的不是做覆盖概率仿真用的。随机几何仿真关心不同仰角下的链路差异最好用一个按仰角变化的简单模型。对于晴空气体衰减可以采用“等效天顶衰减 1/sin(仰角)”的近似。物理依据是大气厚度在低仰角方向被拉长路径近似满足正割关系。下面是一个常用的近似大气衰减函数def atmospheric_loss(elev_deg, zenith_loss_db): 按仰角近似大气衰减。 zenith_loss_db: 天顶方向总大气衰减dB sin_elev np.sin(np.deg2rad(elev_deg)) if sin_elev 1e-6: return 100.0 # 仰角太低直接判为链路中断 return zenith_loss_db / sin_elev典型参考值如下晴空、中纬度地区仅作仿真示例频段天顶大气衰减dB10° 仰角近似衰减dBL / S2GHz左右0.10.6C / X8GHz左右0.21.2Ku12GHz左右0.52.9Ka20GHz左右1.05.8这里给的是“天顶大气衰减”不是雨衰。Ka 频段雨衰可以大到十几 dB要准确算雨衰需要 ITU-R P.618 建议书里的雨区模型输入降水概率、极化和路径仰角比较复杂。我在做预研仿真时一般先不给雨衰只把晴空大气算进去留一个“雨衰余量”的开关需要时再按 99.9% 可用度加固定余量。如果你在复现论文论文里写了“clear sky”就只算气体衰减写了“rain fading”才需要加雨衰模型。别两样混在一起。天线增益处理是另一个容易出错的地方。下行链路用户接收到的信号功率由发射功率 EIRP 和接收天线增益共同决定。站上天线如果是抛物面增益 ( G(\theta) ) 随离轴角变化低轨卫星因为频繁过顶地面站需要波束跟踪不能总按轴向最大增益来算。最简模型是“主瓣高斯近似 等旁瓣包络”def antenna_gain_deg(off_axis_deg, peak_gain_db, beamwidth_deg): 简化的离轴天线增益。 主瓣按二次衰减近似旁瓣设为比峰值低 30dB 的平层。 off np.maximum(np.abs(off_axis_deg), 1e-3) # 主瓣内增益随离轴角平方下降 main_lobe peak_gain_db - 12.0 * (off / beamwidth_deg) ** 2 # 旁瓣平层峰值-30dB side_lobe peak_gain_db - 30.0 return np.maximum(main_lobe, side_lobe)这个函数的思路是在 3dB 波束宽度内增益下降速度大致是 12dB/束宽倍数平方超过一定角度后不再无限下降而是进入旁瓣平台。这样做不会让干扰卫星的天线增益低到离谱干扰功率更接近真实系统。3.3 把损耗合到一起逐链路接收功率计算下面给出一个完整的逐链路功率计算函数。为了后续干扰分析方便返回值同时保留损耗 dB 和接收功率 dBm。def compute_link_power(tx_power_dbm, freq_hz, dist_m, elev_deg, gain_tx_dbi, gain_rx_dbi, zenith_loss_db): 计算单条下行链路的接收功率dBm。 tx_power_dbm: 卫星发射功率dBm gain_tx_dbi: 卫星发射天线增益dBi gain_rx_dbi: 用户接收天线增益dBi fsp free_space_path_loss(dist_m, freq_hz) atm atmospheric_loss(elev_deg, zenith_loss_db) total_loss_db fsp atm - gain_tx_dbi - gain_rx_dbi rx_power_dbm tx_power_dbm - total_loss_db return rx_power_dbm, total_loss_db注意公式里的符号天线增益作为“有效增益”放进预算它降低损耗所以是减去大气衰减增加损耗所以是加上。我在早期版本里写过一次 gain结果所有卫星接收功率变成 159dBm整个星座都在“发热”排查半天才发现是符号问题。后面第 5 章会专门列一条避坑记录。使用这个函数时卫星发射功率典型值可以从几瓦到几十瓦不等换算成 dBm 后大概在 33dBm 到 47dBm 之间。接收天线增益则需要和你用的波束带宽匹配如果仿真里的用户是手持终端增益 0dBi 合理如果是固定地球站按照天线口径和频率用标准公式算峰值增益。传播模型仿真和系统级仿真的区别在损耗这一步就显现了系统级仿真会逐秒更新距离、逐链路算雨衰随机几何仿真则是每个 Monte Carlo 快照中重新生成星座然后对每条链路重算这些损耗。只要保证距离、仰角、大气、天线增益四个量在同一个快照里统计结果才一致。这里特别提醒接收功率 dBm 在干扰分析时不能直接相减因为信号和干扰在线性域相加。后续第 4 章的代码中会用10**(rx_power_dbm/10)把 dBm 转成 mW求和后再转回 dB。这一步应该在链路预算之后做不要在 dB 域直接累加。还有一个注意如果tx_power_dbm是 dBW 而不是 dBm常数要改成 30论文里经常混用建议在仿真开头统一单位。我用一个全局变量UNIT dBm的方式避免后面越改越乱。4. 干扰分析与覆盖概率聚合干扰、SINR 和 Monte Carlo 仿真的完整流程4.1 干扰源集合可见卫星不等于干扰卫星但要先按最坏情况统计下行链路中用户锁定一颗卫星作为服务星其余同频可见卫星都会在用户接收机处产生干扰。什么算“同频”取决于频率复用设计。如果整个星座复用因子为 1所有卫星同频干扰最重如果星座把波束按 4 色或 7 色复用那么同频卫星只是其中的子集。在随机几何仿真里最简单的做法是把所有可见卫星当作干扰源这对应频率复用因子 1可以作为性能上界/下界干扰上界。若要模拟频率复用可以在生成星座后按经纬度或波束编号随机剔除一部分卫星让可见卫星里只有一部分干扰。更精细的做法是引入卫星天线指向的离轴角只有指向用户的波束或旁瓣才产生信号其他波束信号在用户接收天线处衰减更大。从工程角度我建议第一版先做“所有可见卫星均同频”的最坏情况因为这样跑出来的覆盖概率是下限如果下限都满足需求就不用担心复用方案。如果下限不满足再在代码里加一个同频比例开关frac_cochannel用泊松稀疏化将每颗卫星按概率保留为同频干扰源。4.2 SINR定义与从dB到线性的处理对单个用户SINR 的线性形式是[ \mathrm{SINR} \frac{S}{I N_0} ]其中 S 是服务卫星接收功率I 是各路干扰功率之和N0 是接收机噪声功率。算 dB 值时必须先在线性域求和最后再取 10log10。dB 域不能直接做“SdBm - (I_dBm N_dBm)”因为两个功率在 dB 域相加不等于功率和需要把每个功率转回 mW 再相加。噪声功率按经典公式[ N_0 k T B ]其中 k 是玻尔兹曼常数T 是等效噪声温度B 是信道带宽。工程上常用 -174dBm/Hz 表示 kT在 290K 时再加上接收机噪声系数 NF[ N_0_dBm -174 10\log_{10}(B_{Hz}) NF ]4.3 Monte Carlo仿真代码一次完整下行链路SINR统计下面这段代码串起第 2 章和第 3 章全部逻辑输出所有快照的 SINR dB 值。注意里面用np.random.default_rng(seed)保持结果可复现。import numpy as np def simulate_downlink_sinr( mean_sat_num200, height_km550, freq_hz2.0e9, tx_power_dbm40.0, gain_tx_dbi10.0, gain_rx_dbi0.0, min_elev_deg10.0, zenith_loss_db0.1, bandwidth_hz2e6, noise_figure_db2.0, iterations2000, seed42): 基于PPP的低轨星座下行链路SINR仿真。 返回所有快照的SINR dB数组。 rng np.random.default_rng(seed) # 噪声功率底 n0_dbm -174.0 10.0 * np.log10(bandwidth_hz) noise_figure_db n0_mw 10 ** (n0_dbm / 10.0) sinr_list [] for _ in range(iterations): # 1) 生成星座快照 N rng.poisson(mean_sat_num) theta np.arccos(1.0 - 2.0 * rng.random(N)) lon 2.0 * np.pi * rng.random(N) R_orb 6371e3 height_km * 1e3 sat_pos np.column_stack([ R_orb * np.sin(theta) * np.cos(lon), R_orb * np.sin(theta) * np.sin(lon), R_orb * np.cos(theta) ]) # 2) 用户放在赤道某点高度0 user_pos np.array([6371e3, 0.0, 0.0]) # 3) 距离/仰角/可见性 d_vec sat_pos - user_pos dist_m np.linalg.norm(d_vec, axis1) zenith user_pos / np.linalg.norm(user_pos) sin_elev np.sum(d_vec * zenith, axis1) / dist_m sin_elev np.clip(sin_elev, -1.0, 1.0) elev_deg np.degrees(np.arcsin(sin_elev)) vis elev_deg min_elev_deg if np.sum(vis) 0: # 一个可见卫星都没有SINR记作-50dB覆盖概率贡献0 sinr_list.append(-50.0) continue # 4) 计算每条可见链路的接收功率dBm fsp_db 20.0 * np.log10(dist_m[vis]) 20.0 * np.log10(freq_hz) - 147.55 atm_db zenith_loss_db / np.sin(np.deg2rad(elev_deg[vis])) # 简化天线增益这里先固定增益离轴角细化留给进阶 rx_power_dbm tx_power_dbm - fsp_db - atm_db gain_tx_dbi gain_rx_dbi rx_power_mw 10 ** (rx_power_dbm / 10.0) # 5) 服务卫星取接收功率最大的那颗 best_idx np.argmax(rx_power_mw) signal_mw rx_power_mw[best_idx] # 6) 干扰所有其他可见卫星功率之和 interfer_mw np.sum(rx_power_mw) - signal_mw sinr_lin signal_mw / (interfer_mw n0_mw) sinr_dB 10.0 * np.log10(sinr_lin) sinr_list.append(sinr_dB) return np.array(sinr_list)这段代码故意只用 numpy没有引入 scipy方便新手直接跑。逻辑上每一步都在完成一个独立任务第 1 步生成 PPP 星座第 2 步摆用户位置第 3 步筛选仰角第 4 步算损耗和接收功率第 5 步选服务星第 6 步把其他功率当干扰。几个参数需要解释freq_hz2.0e9对应 S 波段适合做低轨通信实验如果你复现 Ka 波段论文改成 20e9并同步调整zenith_loss_db和天线口径。min_elev_deg10.0但大气损耗仍按/sin(elev)计算仰角越小损耗越大这个模型中不需要再额外加遮挡。iterations2000通常够看趋势如果要画平滑 CDF建议 5000 以上。第 6 章会介绍量化怎么才算“够”。在选择服务卫星时我按“接收功率最大”而不是“距离最近”。因为损耗还包括大气和天线增益接收功率最大更符合实际只有在自由空间损耗占绝对主导时两者才等价。4.4 从SINR样本到覆盖概率/中断概率得到 SINR 数组后统计目标很简单def coverage_probability(sinr_db, threshold_db): return np.mean(sinr_db threshold_db) # 示例SINR 0dB 的覆盖概率 sinr_db simulate_downlink_sinr(iterations2000) print(P(SINR 0 dB) , coverage_probability(sinr_db, 0.0)) print(P(SINR 10 dB) , coverage_probability(sinr_db, 10.0))因为np.mean(bool数组)就是取 True 的比例。覆盖概率就是 CDF 的补( P_c(\tau) P(\mathrm{SINR} \tau) )中断概率等于 ( 1 - P_c(\tau) )。这里有一个容易被忽略的点阈值 label论文里可能用 dB可能用线性值写代码时要统一。我在做参数扫描时习惯把threshold_db存成列表[0, 5, 10, 15]一次性输出一条覆盖概率阶梯。为了检验干扰的贡献可以在同一组快照下分别跑“含干扰”和“不含干扰”两种模式。把第 6 行interfer_mw替换成 0.0得到的 SINR 就是信号只在噪声下的信噪比 SNR两条曲线之间的差距就是同频干扰带来的损失。这是论文图表里最常见的对比之一。如果要模拟频率复用可以把第 6 步改成对干扰卫星生成一个伯努利掩码cochannel rng.random(sum(vis)) frac_cochannel服务卫星强制保留其余只有掩码为 True 的卫星计入干扰。这样改既不破坏原有链路又能快速观察复用因子对覆盖概率的影响。5. 随机几何仿真避坑与排查5个反复出现的翻车现场下面这些坑是我把随机几何卫星仿真从论文公式搬到代码时真实踩过的。每一条都按“现象 → 原因 → 解决”写排列顺序基本按出现频率。如果你跑出来的结果和论文对不上先对照这一章排查。5.1 覆盖概率曲线剧烈波动仿真发散到没法看现象固定参数跑两次覆盖概率一个 0.35、一个 0.75迭代次数加到 10000 也只是稍微稳一点曲线还是锯齿状。有时候把脚本换一台机器重跑数值又变了。原因不是随机数种子问题而是低轨场景中可见卫星数量太少时PPP 生成的卫星位置方差太大。尤其在 min_elev 取得较高、卫星期望数较小时某个快照可能一颗可见卫星都没有下一个快照却有 5 颗SINR 样本分布呈重尾。这种“仿真发散”本质上是样本量不够覆盖极端事件。解决先检查“每个快照可见卫星数”的直方图。如果平均可见数小于 3说明星座密度和最低仰角不匹配仿真结果不稳定是正常的需要提高卫星期望数或降低最低仰角。另外不要对少数几个快照求均值应累计大量样本做 CDF再读取百分位。也可以将 PPP 换成 BPP 固定卫星总数减小卫星数波动带来的影响。我自己的习惯是在加iterations之前先跑一次可见卫星数统计确认均值在 5 颗以上再继续。5.2 天线增益符号反了接收功率高到离奇现象所有链路接收功率都是正的几十 dBm信号直接“放大”而不是衰减覆盖概率接近 1。甚至发射功率设成 0 也能收到正功率。原因链路预算是tx_power - loss gain但代码写成tx_power - (fsp atm gain)将天线增益当成了损耗的一部分。解决用单颗卫星自检。假设距离 1000km、频率 2GHz、发射功率 40dBm、发射和接收增益各 10dBi自由空间损耗约 158.5dB接收功率应在 -98.5dBm 量级。如果结果偏离超过几 dB先检查增益符号。另外一个自检方法把发射天线和接收天线增益都设为 0dBi、大气衰减设 0跑一次“无增益”仿真看接收功率是否等于tx_power - fsp。这样能快速定位是符号还是单位问题。5.3 低仰角卫星被选作服务星覆盖虚高现象覆盖概率在 10° 仰角门限附近异常高实际系统在低仰角很难建链仿真结果却比论文乐观很多。打开日志发现服务卫星经常是地平线附近那颗。原因只按自由空间损耗选星没有充分考虑大气和地形衰减。前面损耗函数里虽然用1/sin(elev)补偿大气但如果门限设置过低比如 min_elev5°10° 以下甚至 5° 链路仍会被保留在城市或山区场景低仰角还有地形遮挡几何可见不等价于电波可见。解决将 min_elev 提到 15° 或更高或在接收功率中增加一个低仰角衰落项比如在仰角小于 20° 时增加一个与仰角相关的附加损耗。对于论文复现先看一眼论文用的是哪个最低仰角很多低轨卫星论文用 10° 或 15°不要凭空改成 5° 去优化覆盖概率。记住“更低仰角 更高覆盖”在仿真里是数学结果工程上是陷阱。5.4 PPP边界效应用平面近似球面导致距离计算错现象用户在某个经纬度附近卫星干扰统计的 CDF 在远距离处出现异常平台覆盖概率呈现出和卫星密度不符的饱和。曲线形状在经度 90° 和赤道附近差异很大。原因如果直接用二维平面 PPP 在用户周围生成“局部卫星”没有考虑地球曲率那么远离用户的平面点在真实球面上可能已经在地平线下但仍然被当成干扰源或者说平面距离和球面大圆距离不一致导致损耗偏小。解决始终使用球面坐标生成卫星并按仰角过滤。不要在用户位置建立一个边长为几百公里的方形区域然后往里丢卫星。对低轨星座来说干扰源可以出现在地平线以外但球面 PPP 天然把位置限制在轨道球面上不需要额外设范围。检查方法把卫星到用户的最小仰角分布画出来如果有很多卫星仰角为负值且仍参与计算说明过滤没生效。注意用平面 PPP 也不是完全不能做但适用的场景是“用户被极大星座淹没干扰主要来自近处卫星”只要研究低轨星座轨道高度才几百公里地球曲率在小尺度内已经能造成明显距离误差别再省这一步。5.5 噪声带宽和dB单位混用SINR曲线整体平移现象SINR 分布曲线形状正确但整体向右或向左偏移几个 dB似乎所有卫星功率都变了。对比不同论文时同一组参数下差 10dB。原因带宽B用了 MHz 而不是 Hz发射功率用了 dBW 但代码按 dBm 处理噪声系数忘了加。单位偏差 10dB 很典型。解决在仿真文件开头集中定义单位换算DBM_TO_MW 1e-3全部用 dBm带宽必须换算成 Hz 后再取 log噪声功率用-174 10log10(BW_Hz) NF。可以做一次无干扰的 SNR 仿真将 SNR 与理论链路预算比对能整体平移说明单位问题局部偏差才是随机性问题。我在项目里会把所有功率变量名带上_dbm后缀方便肉眼检查。6. 结果验证与参数扫描让随机几何仿真结果可复现、可信6.1 一个轻量的收敛性检验批次平均而不是只看总均值Monte Carlo 仿真的覆盖概率是一个随机估计量如果只跑 2000 次就报数别人换了随机种子结果就对不上。我常用的验证方法是把样本分成 M 批每批 N 次计算每批的覆盖概率再看这些批次的均值和标准差。标准差就是估计量的标准误。下面是一个只依赖 numpy 的批次检验函数def batch_estimate(sinr_db, threshold_db, batch_size500): n_batch len(sinr_db) // batch_size batch_probs [] for i in range(n_batch): seg sinr_db[i*batch_size:(i1)*batch_size] batch_probs.append(np.mean(seg threshold_db)) return np.mean(batch_probs), np.std(batch_probs, ddof1)如果批次均值稳定而批次标准差小于 0.01这个覆盖概率的小数点后两位才是可以报的。如果标准差在 0.05 以上需要增加迭代次数或改换 BPP 固定卫星数否则任何参数对比都没有意义。这个检验动作应该是参数扫描之前的固定动作而不是最后补一个“置信区间”了事。6.2 参数扫描用同一个随机种子做“可复现对比”随机几何最常见的结果图是“覆盖概率 vs 星座密度”或“覆盖概率 vs 阈值”。这类对比必须保证所有曲线在哪一轮快照上使用同一套随机数才能让曲线差异只来自参数本身。最简单做法是在函数里加seed参数并且每次调用把 seed 偏移一个固定常数。for density_factor in [0.5, 1.0, 2.0]: base_mean 200 sinr simulate_downlink_sinr( mean_sat_numbase_mean * density_factor, iterations5000, seed100 int(density_factor * 100)) print(density_factor, coverage_probability(sinr, 10.0))这样跑出来的曲线是平滑且可复现的。更高阶的做法是把不同参数组合写成循环最后用 pandas 存成 CSV。我不建议在同一个 for 循环里反复修改全局变量很容易把上一轮的星座带到下一轮。参数扫描时我习惯先看几个典型方向的趋势对照物理直觉检查参数调大对覆盖概率的典型影响排查注意卫星密度干扰也增大覆盖不一定单调上升看干扰受限区间轨道高度单链路损耗增大但可见卫星更多综合影响是非线性最低仰角可见卫星减少覆盖下降别用来“优化”结果频率复用因子同频干扰减少覆盖上升需建模波束方向至少就我自己而言这类仿真项目从“能跑出数”到“能相信这个数”中间隔的不是代码量而是验证方法。先确认链路预算没有符号错误再查可见卫星统计最后才用批次标准误判断收敛。如果在每一步都把中间量打印出来看一眼翻车的概率会小很多。希望帮到你。本文还有配套的精品资源点击获取 低轨卫星通信随机几何仿真：PPP建模与下行链路覆盖分析
简介面向卫星通信研究者和工程师这份资源复现了基于随机几何的低轨星座下行链路仿真与分析工作。内容围绕低轨卫星网络移动性强、电磁环境复杂的特点基于二项点过程BPP构建星座模型覆盖单星与多星场景的路径损耗、接收功率、噪声及干扰计算并通过干扰期望计算方法评估星地链路性能。包体为1个docx文档压缩后约51KB文档内含详细理论推导和可直接运行的Python代码涵盖参数设置、卫星与地面站生成、SINR计算、干扰分析与可视化配有逐步解释方便对照学习。目前已有105人学习。借助该文档可理解随机几何在星座建模中的应用掌握星座参数对通信性能的影响为巨型低轨星座的干扰建模和链路设计提供可扩展的仿真参考。1. 低轨卫星通信的随机几何仿真先回答“为什么用PPP做下行链路”低轨卫星通信的下行链路仿真最花时间的往往不是轨道积分而是干扰说不清。真实低轨星座有几十个轨道面每颗卫星都在移动用户天线还要跟着波束切换如果只用精确星历逐秒跑结果只是某个具体时刻、某个纬度的特解。随机几何把卫星位置建模成球面上的泊松点过程PPP用密度和高度代替轨道参数覆盖概率和信干噪比分布这类指标第一次可以从统计上直接回答。这套方法适合正在复现低轨卫星通信论文、准备毕业设计或者做频率复用预研的工程师。下面按“星座建模 → 损耗计算 → 干扰分析 → 避坑 → 验证”的顺序给出可以直接改参数跑的 Python 代码和工程坑。2. 星座建模与可见性判断用球面泊松点过程生成低轨星座的代码实现2.1 为什么低轨星座适合用随机几何的PPP近似低轨卫星通信系统的轨道高度通常在 500km 到 1200km 之间。与静止轨道不同每颗低轨卫星相对地面用户移动很快几分钟内会从地平线附近升到头顶再降到另一侧地平线。对一个固定位置的地面用户来说卫星在哪个方向、距离多远是时刻变化的星座设计者更关心“任意时刻有至少一颗卫星可见”的概率而不是某颗卫星某秒的精确坐标。随机几何正是为这类统计问题准备的。常见的低轨星座建模路线有两种。确定性轨道建模需要输入轨道面、相位、TLE 或者开普勒根数适合做星间切换、波束指向的事件级仿真但一次只能代表一套具体星座构型。随机几何建模忽略卫星编号只关心单位面积上期望有多少颗卫星、在统计上服从什么分布输出的是覆盖概率、中断概率这类系统级指标。两者不矛盾确定性建模用来校核特定方案随机几何用来回答“如果星座密度提升一倍覆盖概率怎么变”这类趋势问题。实际仿真中使用 PPP 还有一个物理上的合理性。真实低轨星座虽然结构确定但经过一段时间后对地面用户而言卫星到达方向会随轨道运动和用户位置变化形成近似随机的模式。PPP 的独立性让干扰分析中的聚合干扰拉普拉斯变换可以分解成单卫星的乘积这是随机几何理论推导的基础。需要精确匹配卫星总数时可以用 BPP 固定每轮快照的卫星数这也被大量论文当作 Monte Carlo 验证方式。2.2 在球面上生成PPP卫星的代码先给出核心代码。注意这里的“卫星”指某时刻能与用户发生链路的可见卫星不是轨道上的全部卫星。我们用一个固定高度球面来近似所有卫星位置。import numpy as np def simulate_leo_sat_positions(mean_sat_num, height_km): 在球面上生成泊松点过程的卫星位置。 参数 - mean_sat_num : floatPPP的期望卫星数即密度乘以球面面积 - height_km : float轨道高度km 返回 - sat_pos : ndarray, shape (N,3)地心直角坐标单位 m R_earth 6371.0e3 # 地球半径m R_orbit R_earth height_km * 1e3 # 泊松点过程卫星数是随机变量先抽一个泊松数 N np.random.poisson(mean_sat_num) # 均匀抽样球面上的点。注意经度可以直接均匀极角必须用反余弦 lon 2.0 * np.pi * np.random.rand(N) cos_theta 1.0 - 2.0 * np.random.rand(N) # theta ∈ [0, pi] 的余弦值均匀在[-1,1] sin_theta np.sqrt(1.0 - cos_theta ** 2) # 球坐标 (r, theta, lon) 转地心直角坐标 x R_orbit * sin_theta * np.cos(lon) y R_orbit * sin_theta * np.sin(lon) z R_orbit * cos_theta return np.column_stack([x, y, z])这里最关键的是球面均匀抽样。很多人第一次写会直接让纬度和经度均匀分布结果高纬度地区卫星密度偏大造成干扰分析失真用cos_theta 1 - 2*rand等效于先均匀抽样单位球面切向量再做反变换能保证在每个相等面积上落入相同期望点数。mean_sat_num就是 PPP 的强度参数 λ 乘上球面面积 (4\pi R_{orbit}^2)。如果你要匹配一个具体星座比如“轨道面 15 个每个轨道面 40 颗”那期望卫星数就是 600但 PPP 的瞬时卫星数会随机波动这正是随机几何与确定性建模的差别。如果你想让卫星总数固定就把上面生成 N 的那一行换成N int(mean_sat_num)其余不用动这就变成了 BPP。做论文对比时我建议两边都跑用 PPP 做理论推导用 BPP 做蒙特卡洛验证两者差异就是随机几何近似的代价。2.3 地面用户到卫星的距离、仰角与可见性判断有了卫星位置下一步是判断地面用户能看到哪些卫星。这里需要把用户经纬度也变成地心直角坐标然后计算用户到卫星的距离矢量和仰角。def llh_to_ecef(lat_deg, lon_deg, alt_m): 把大地坐标转地心直角坐标简单球体模型 R 6371.0e3 alt_m lat np.deg2rad(lat_deg) lon np.deg2rad(lon_deg) cos_lat np.cos(lat) return np.array([ R * cos_lat * np.cos(lon), R * cos_lat * np.sin(lon), R * np.sin(lat) ]) def compute_distance_elevation(user_ecef, sat_pos): 计算每个卫星到用户的距离和仰角。 返回 distance_m, elev_deg d_vec sat_pos - user_ecef # 从用户指向卫星 dist_m np.linalg.norm(d_vec, axis1) # 用户位置指向地心的单位向量也就是用户本地天顶方向 zenith user_ecef / np.linalg.norm(user_ecef) # 仰角 sin 值 距离向量在天顶方向上的投影 / 距离 sin_elev np.sum(d_vec * zenith, axis1) / dist_m sin_elev np.clip(sin_elev, -1.0, 1.0) elev_deg np.degrees(np.arcsin(sin_elev)) return dist_m, elev_deg这段代码的几何关系可以这样理解用户站在球面上本地“竖直向上”的方向就是用户位置向量向外延长的那根线即user_ecef的单位向量。从用户指向卫星的距离矢量和竖直方向做点积得到的是距离在竖直方向的投影投影除以距离就是仰角的正弦。卫星在正头顶时仰角 90°在地平线附近仰角接近 0°跑在另一边就会算出负值应当丢弃。有了距离和仰角就可以根据最低仰角做筛选。工程上传统卫星通信通常要求最低通信仰角在 8° 到 25° 之间低轨卫星为了增加覆盖可能放宽到 5°但低于 5° 的链路大气损耗大、遮挡多不建议计入统计。筛选逻辑很简单def get_visible_satellites(user_ecef, sat_pos, min_elev_deg10.0): dist_m, elev_deg compute_distance_elevation(user_ecef, sat_pos) vis elev_deg min_elev_deg return dist_m[vis], elev_deg[vis]筛选完剩余卫星数通常在几颗到几十颗。如果平均剩下不到一颗说明星座密度太低或最低仰角设得太高结果会以“脱网”场景为主覆盖概率整体偏低。这个现象在第 5 章的坑里再展开。2.4 星座建模的边界PPP的前提与参数标定用 PPP 做低轨星座建模有一个物理前提需要交代清楚PPP 要求卫星位置在球面上均匀且独立。真实星座由有限个轨道面构成存在明显的结构相关性某些地区在某段时间可能“看见”多颗星另一些地区则稀疏。如果把真实星座的瞬时快照跑几千个时刻平均下来覆盖概率大概率落在 PPP 仿真结果的置信区间内但如果论文要比较“轨道面数对性能的影响”PPP 就表达不了轨道结构。因此这里有两个常见标定方法。一是把 PPP 的期望卫星数设为星座卫星总数用于评估系统容量上限二是把 PPP 期望卫星数设为“某个时刻用户可见或可干扰到的同频卫星数”用于评估干扰受限场景。前者适合覆盖分析后者适合干扰分析。实际项目中我更倾向于按轨道高度和天线波束宽度先算一个“干扰可见圈”再在圈内用 PPP 建模这样能减少全局 PPP 带来的过度保守干扰。3. 下行链路损耗计算自由空间、大气衰减与天线增益怎么进链路预算3.1 先算自由空间路径损耗斜距不是高度别拿轨道高度直接填低轨卫星与地面用户的距离随仰角变化非常大。以 550km 轨道高度为例卫星过顶时距离约 550km但仰角 10° 时斜距大约 2100km距离相差近 4 倍自由空间损耗相差约 12dB。很多初版仿真直接把h当距离填进弗里斯公式结果覆盖概率虚高干扰分析也失去意义。正确的做法是使用前面算出的dist_m先计算自由空间路径损耗[ FSPL(dB) 20\log_{10}(d) 20\log_{10}(f) - 147.55 ]其中d单位是米f单位是 Hz147.55来自 ( (4\pi/c)^2 ) 的对数变换。注意这个式子里的常数如果你用 km 和 GHz常数会变成 32.45单位搞错是最常见的低级错误。下面这段代码把距离、频率、仰角一并处理并保留逐链路的损耗供后续干扰求和def free_space_path_loss(distance_m, freq_hz): 自由空间路径损耗单位 dB return 20.0 * np.log10(distance_m) 20.0 * np.log10(freq_hz) - 147.553.2 大气衰减、雨衰和天线增益从单值补偿到按仰角插值LEO 下行链路的工作频率从 L 波段到 Ka 波段都有。频率越高大气衰减越明显。工程上最粗略的做法是在链路预算里加一个固定大气损耗余量比如 L 波段加 0.5dB、Ka 波段加 2dB但这是做系统预算用的不是做覆盖概率仿真用的。随机几何仿真关心不同仰角下的链路差异最好用一个按仰角变化的简单模型。对于晴空气体衰减可以采用“等效天顶衰减 1/sin(仰角)”的近似。物理依据是大气厚度在低仰角方向被拉长路径近似满足正割关系。下面是一个常用的近似大气衰减函数def atmospheric_loss(elev_deg, zenith_loss_db): 按仰角近似大气衰减。 zenith_loss_db: 天顶方向总大气衰减dB sin_elev np.sin(np.deg2rad(elev_deg)) if sin_elev 1e-6: return 100.0 # 仰角太低直接判为链路中断 return zenith_loss_db / sin_elev典型参考值如下晴空、中纬度地区仅作仿真示例频段天顶大气衰减dB10° 仰角近似衰减dBL / S2GHz左右0.10.6C / X8GHz左右0.21.2Ku12GHz左右0.52.9Ka20GHz左右1.05.8这里给的是“天顶大气衰减”不是雨衰。Ka 频段雨衰可以大到十几 dB要准确算雨衰需要 ITU-R P.618 建议书里的雨区模型输入降水概率、极化和路径仰角比较复杂。我在做预研仿真时一般先不给雨衰只把晴空大气算进去留一个“雨衰余量”的开关需要时再按 99.9% 可用度加固定余量。如果你在复现论文论文里写了“clear sky”就只算气体衰减写了“rain fading”才需要加雨衰模型。别两样混在一起。天线增益处理是另一个容易出错的地方。下行链路用户接收到的信号功率由发射功率 EIRP 和接收天线增益共同决定。站上天线如果是抛物面增益 ( G(\theta) ) 随离轴角变化低轨卫星因为频繁过顶地面站需要波束跟踪不能总按轴向最大增益来算。最简模型是“主瓣高斯近似 等旁瓣包络”def antenna_gain_deg(off_axis_deg, peak_gain_db, beamwidth_deg): 简化的离轴天线增益。 主瓣按二次衰减近似旁瓣设为比峰值低 30dB 的平层。 off np.maximum(np.abs(off_axis_deg), 1e-3) # 主瓣内增益随离轴角平方下降 main_lobe peak_gain_db - 12.0 * (off / beamwidth_deg) ** 2 # 旁瓣平层峰值-30dB side_lobe peak_gain_db - 30.0 return np.maximum(main_lobe, side_lobe)这个函数的思路是在 3dB 波束宽度内增益下降速度大致是 12dB/束宽倍数平方超过一定角度后不再无限下降而是进入旁瓣平台。这样做不会让干扰卫星的天线增益低到离谱干扰功率更接近真实系统。3.3 把损耗合到一起逐链路接收功率计算下面给出一个完整的逐链路功率计算函数。为了后续干扰分析方便返回值同时保留损耗 dB 和接收功率 dBm。def compute_link_power(tx_power_dbm, freq_hz, dist_m, elev_deg, gain_tx_dbi, gain_rx_dbi, zenith_loss_db): 计算单条下行链路的接收功率dBm。 tx_power_dbm: 卫星发射功率dBm gain_tx_dbi: 卫星发射天线增益dBi gain_rx_dbi: 用户接收天线增益dBi fsp free_space_path_loss(dist_m, freq_hz) atm atmospheric_loss(elev_deg, zenith_loss_db) total_loss_db fsp atm - gain_tx_dbi - gain_rx_dbi rx_power_dbm tx_power_dbm - total_loss_db return rx_power_dbm, total_loss_db注意公式里的符号天线增益作为“有效增益”放进预算它降低损耗所以是减去大气衰减增加损耗所以是加上。我在早期版本里写过一次 gain结果所有卫星接收功率变成 159dBm整个星座都在“发热”排查半天才发现是符号问题。后面第 5 章会专门列一条避坑记录。使用这个函数时卫星发射功率典型值可以从几瓦到几十瓦不等换算成 dBm 后大概在 33dBm 到 47dBm 之间。接收天线增益则需要和你用的波束带宽匹配如果仿真里的用户是手持终端增益 0dBi 合理如果是固定地球站按照天线口径和频率用标准公式算峰值增益。传播模型仿真和系统级仿真的区别在损耗这一步就显现了系统级仿真会逐秒更新距离、逐链路算雨衰随机几何仿真则是每个 Monte Carlo 快照中重新生成星座然后对每条链路重算这些损耗。只要保证距离、仰角、大气、天线增益四个量在同一个快照里统计结果才一致。这里特别提醒接收功率 dBm 在干扰分析时不能直接相减因为信号和干扰在线性域相加。后续第 4 章的代码中会用10**(rx_power_dbm/10)把 dBm 转成 mW求和后再转回 dB。这一步应该在链路预算之后做不要在 dB 域直接累加。还有一个注意如果tx_power_dbm是 dBW 而不是 dBm常数要改成 30论文里经常混用建议在仿真开头统一单位。我用一个全局变量UNIT dBm的方式避免后面越改越乱。4. 干扰分析与覆盖概率聚合干扰、SINR 和 Monte Carlo 仿真的完整流程4.1 干扰源集合可见卫星不等于干扰卫星但要先按最坏情况统计下行链路中用户锁定一颗卫星作为服务星其余同频可见卫星都会在用户接收机处产生干扰。什么算“同频”取决于频率复用设计。如果整个星座复用因子为 1所有卫星同频干扰最重如果星座把波束按 4 色或 7 色复用那么同频卫星只是其中的子集。在随机几何仿真里最简单的做法是把所有可见卫星当作干扰源这对应频率复用因子 1可以作为性能上界/下界干扰上界。若要模拟频率复用可以在生成星座后按经纬度或波束编号随机剔除一部分卫星让可见卫星里只有一部分干扰。更精细的做法是引入卫星天线指向的离轴角只有指向用户的波束或旁瓣才产生信号其他波束信号在用户接收天线处衰减更大。从工程角度我建议第一版先做“所有可见卫星均同频”的最坏情况因为这样跑出来的覆盖概率是下限如果下限都满足需求就不用担心复用方案。如果下限不满足再在代码里加一个同频比例开关frac_cochannel用泊松稀疏化将每颗卫星按概率保留为同频干扰源。4.2 SINR定义与从dB到线性的处理对单个用户SINR 的线性形式是[ \mathrm{SINR} \frac{S}{I N_0} ]其中 S 是服务卫星接收功率I 是各路干扰功率之和N0 是接收机噪声功率。算 dB 值时必须先在线性域求和最后再取 10log10。dB 域不能直接做“SdBm - (I_dBm N_dBm)”因为两个功率在 dB 域相加不等于功率和需要把每个功率转回 mW 再相加。噪声功率按经典公式[ N_0 k T B ]其中 k 是玻尔兹曼常数T 是等效噪声温度B 是信道带宽。工程上常用 -174dBm/Hz 表示 kT在 290K 时再加上接收机噪声系数 NF[ N_0_dBm -174 10\log_{10}(B_{Hz}) NF ]4.3 Monte Carlo仿真代码一次完整下行链路SINR统计下面这段代码串起第 2 章和第 3 章全部逻辑输出所有快照的 SINR dB 值。注意里面用np.random.default_rng(seed)保持结果可复现。import numpy as np def simulate_downlink_sinr( mean_sat_num200, height_km550, freq_hz2.0e9, tx_power_dbm40.0, gain_tx_dbi10.0, gain_rx_dbi0.0, min_elev_deg10.0, zenith_loss_db0.1, bandwidth_hz2e6, noise_figure_db2.0, iterations2000, seed42): 基于PPP的低轨星座下行链路SINR仿真。 返回所有快照的SINR dB数组。 rng np.random.default_rng(seed) # 噪声功率底 n0_dbm -174.0 10.0 * np.log10(bandwidth_hz) noise_figure_db n0_mw 10 ** (n0_dbm / 10.0) sinr_list [] for _ in range(iterations): # 1) 生成星座快照 N rng.poisson(mean_sat_num) theta np.arccos(1.0 - 2.0 * rng.random(N)) lon 2.0 * np.pi * rng.random(N) R_orb 6371e3 height_km * 1e3 sat_pos np.column_stack([ R_orb * np.sin(theta) * np.cos(lon), R_orb * np.sin(theta) * np.sin(lon), R_orb * np.cos(theta) ]) # 2) 用户放在赤道某点高度0 user_pos np.array([6371e3, 0.0, 0.0]) # 3) 距离/仰角/可见性 d_vec sat_pos - user_pos dist_m np.linalg.norm(d_vec, axis1) zenith user_pos / np.linalg.norm(user_pos) sin_elev np.sum(d_vec * zenith, axis1) / dist_m sin_elev np.clip(sin_elev, -1.0, 1.0) elev_deg np.degrees(np.arcsin(sin_elev)) vis elev_deg min_elev_deg if np.sum(vis) 0: # 一个可见卫星都没有SINR记作-50dB覆盖概率贡献0 sinr_list.append(-50.0) continue # 4) 计算每条可见链路的接收功率dBm fsp_db 20.0 * np.log10(dist_m[vis]) 20.0 * np.log10(freq_hz) - 147.55 atm_db zenith_loss_db / np.sin(np.deg2rad(elev_deg[vis])) # 简化天线增益这里先固定增益离轴角细化留给进阶 rx_power_dbm tx_power_dbm - fsp_db - atm_db gain_tx_dbi gain_rx_dbi rx_power_mw 10 ** (rx_power_dbm / 10.0) # 5) 服务卫星取接收功率最大的那颗 best_idx np.argmax(rx_power_mw) signal_mw rx_power_mw[best_idx] # 6) 干扰所有其他可见卫星功率之和 interfer_mw np.sum(rx_power_mw) - signal_mw sinr_lin signal_mw / (interfer_mw n0_mw) sinr_dB 10.0 * np.log10(sinr_lin) sinr_list.append(sinr_dB) return np.array(sinr_list)这段代码故意只用 numpy没有引入 scipy方便新手直接跑。逻辑上每一步都在完成一个独立任务第 1 步生成 PPP 星座第 2 步摆用户位置第 3 步筛选仰角第 4 步算损耗和接收功率第 5 步选服务星第 6 步把其他功率当干扰。几个参数需要解释freq_hz2.0e9对应 S 波段适合做低轨通信实验如果你复现 Ka 波段论文改成 20e9并同步调整zenith_loss_db和天线口径。min_elev_deg10.0但大气损耗仍按/sin(elev)计算仰角越小损耗越大这个模型中不需要再额外加遮挡。iterations2000通常够看趋势如果要画平滑 CDF建议 5000 以上。第 6 章会介绍量化怎么才算“够”。在选择服务卫星时我按“接收功率最大”而不是“距离最近”。因为损耗还包括大气和天线增益接收功率最大更符合实际只有在自由空间损耗占绝对主导时两者才等价。4.4 从SINR样本到覆盖概率/中断概率得到 SINR 数组后统计目标很简单def coverage_probability(sinr_db, threshold_db): return np.mean(sinr_db threshold_db) # 示例SINR 0dB 的覆盖概率 sinr_db simulate_downlink_sinr(iterations2000) print(P(SINR 0 dB) , coverage_probability(sinr_db, 0.0)) print(P(SINR 10 dB) , coverage_probability(sinr_db, 10.0))因为np.mean(bool数组)就是取 True 的比例。覆盖概率就是 CDF 的补( P_c(\tau) P(\mathrm{SINR} \tau) )中断概率等于 ( 1 - P_c(\tau) )。这里有一个容易被忽略的点阈值 label论文里可能用 dB可能用线性值写代码时要统一。我在做参数扫描时习惯把threshold_db存成列表[0, 5, 10, 15]一次性输出一条覆盖概率阶梯。为了检验干扰的贡献可以在同一组快照下分别跑“含干扰”和“不含干扰”两种模式。把第 6 行interfer_mw替换成 0.0得到的 SINR 就是信号只在噪声下的信噪比 SNR两条曲线之间的差距就是同频干扰带来的损失。这是论文图表里最常见的对比之一。如果要模拟频率复用可以把第 6 步改成对干扰卫星生成一个伯努利掩码cochannel rng.random(sum(vis)) frac_cochannel服务卫星强制保留其余只有掩码为 True 的卫星计入干扰。这样改既不破坏原有链路又能快速观察复用因子对覆盖概率的影响。5. 随机几何仿真避坑与排查5个反复出现的翻车现场下面这些坑是我把随机几何卫星仿真从论文公式搬到代码时真实踩过的。每一条都按“现象 → 原因 → 解决”写排列顺序基本按出现频率。如果你跑出来的结果和论文对不上先对照这一章排查。5.1 覆盖概率曲线剧烈波动仿真发散到没法看现象固定参数跑两次覆盖概率一个 0.35、一个 0.75迭代次数加到 10000 也只是稍微稳一点曲线还是锯齿状。有时候把脚本换一台机器重跑数值又变了。原因不是随机数种子问题而是低轨场景中可见卫星数量太少时PPP 生成的卫星位置方差太大。尤其在 min_elev 取得较高、卫星期望数较小时某个快照可能一颗可见卫星都没有下一个快照却有 5 颗SINR 样本分布呈重尾。这种“仿真发散”本质上是样本量不够覆盖极端事件。解决先检查“每个快照可见卫星数”的直方图。如果平均可见数小于 3说明星座密度和最低仰角不匹配仿真结果不稳定是正常的需要提高卫星期望数或降低最低仰角。另外不要对少数几个快照求均值应累计大量样本做 CDF再读取百分位。也可以将 PPP 换成 BPP 固定卫星总数减小卫星数波动带来的影响。我自己的习惯是在加iterations之前先跑一次可见卫星数统计确认均值在 5 颗以上再继续。5.2 天线增益符号反了接收功率高到离奇现象所有链路接收功率都是正的几十 dBm信号直接“放大”而不是衰减覆盖概率接近 1。甚至发射功率设成 0 也能收到正功率。原因链路预算是tx_power - loss gain但代码写成tx_power - (fsp atm gain)将天线增益当成了损耗的一部分。解决用单颗卫星自检。假设距离 1000km、频率 2GHz、发射功率 40dBm、发射和接收增益各 10dBi自由空间损耗约 158.5dB接收功率应在 -98.5dBm 量级。如果结果偏离超过几 dB先检查增益符号。另外一个自检方法把发射天线和接收天线增益都设为 0dBi、大气衰减设 0跑一次“无增益”仿真看接收功率是否等于tx_power - fsp。这样能快速定位是符号还是单位问题。5.3 低仰角卫星被选作服务星覆盖虚高现象覆盖概率在 10° 仰角门限附近异常高实际系统在低仰角很难建链仿真结果却比论文乐观很多。打开日志发现服务卫星经常是地平线附近那颗。原因只按自由空间损耗选星没有充分考虑大气和地形衰减。前面损耗函数里虽然用1/sin(elev)补偿大气但如果门限设置过低比如 min_elev5°10° 以下甚至 5° 链路仍会被保留在城市或山区场景低仰角还有地形遮挡几何可见不等价于电波可见。解决将 min_elev 提到 15° 或更高或在接收功率中增加一个低仰角衰落项比如在仰角小于 20° 时增加一个与仰角相关的附加损耗。对于论文复现先看一眼论文用的是哪个最低仰角很多低轨卫星论文用 10° 或 15°不要凭空改成 5° 去优化覆盖概率。记住“更低仰角 更高覆盖”在仿真里是数学结果工程上是陷阱。5.4 PPP边界效应用平面近似球面导致距离计算错现象用户在某个经纬度附近卫星干扰统计的 CDF 在远距离处出现异常平台覆盖概率呈现出和卫星密度不符的饱和。曲线形状在经度 90° 和赤道附近差异很大。原因如果直接用二维平面 PPP 在用户周围生成“局部卫星”没有考虑地球曲率那么远离用户的平面点在真实球面上可能已经在地平线下但仍然被当成干扰源或者说平面距离和球面大圆距离不一致导致损耗偏小。解决始终使用球面坐标生成卫星并按仰角过滤。不要在用户位置建立一个边长为几百公里的方形区域然后往里丢卫星。对低轨星座来说干扰源可以出现在地平线以外但球面 PPP 天然把位置限制在轨道球面上不需要额外设范围。检查方法把卫星到用户的最小仰角分布画出来如果有很多卫星仰角为负值且仍参与计算说明过滤没生效。注意用平面 PPP 也不是完全不能做但适用的场景是“用户被极大星座淹没干扰主要来自近处卫星”只要研究低轨星座轨道高度才几百公里地球曲率在小尺度内已经能造成明显距离误差别再省这一步。5.5 噪声带宽和dB单位混用SINR曲线整体平移现象SINR 分布曲线形状正确但整体向右或向左偏移几个 dB似乎所有卫星功率都变了。对比不同论文时同一组参数下差 10dB。原因带宽B用了 MHz 而不是 Hz发射功率用了 dBW 但代码按 dBm 处理噪声系数忘了加。单位偏差 10dB 很典型。解决在仿真文件开头集中定义单位换算DBM_TO_MW 1e-3全部用 dBm带宽必须换算成 Hz 后再取 log噪声功率用-174 10log10(BW_Hz) NF。可以做一次无干扰的 SNR 仿真将 SNR 与理论链路预算比对能整体平移说明单位问题局部偏差才是随机性问题。我在项目里会把所有功率变量名带上_dbm后缀方便肉眼检查。6. 结果验证与参数扫描让随机几何仿真结果可复现、可信6.1 一个轻量的收敛性检验批次平均而不是只看总均值Monte Carlo 仿真的覆盖概率是一个随机估计量如果只跑 2000 次就报数别人换了随机种子结果就对不上。我常用的验证方法是把样本分成 M 批每批 N 次计算每批的覆盖概率再看这些批次的均值和标准差。标准差就是估计量的标准误。下面是一个只依赖 numpy 的批次检验函数def batch_estimate(sinr_db, threshold_db, batch_size500): n_batch len(sinr_db) // batch_size batch_probs [] for i in range(n_batch): seg sinr_db[i*batch_size:(i1)*batch_size] batch_probs.append(np.mean(seg threshold_db)) return np.mean(batch_probs), np.std(batch_probs, ddof1)如果批次均值稳定而批次标准差小于 0.01这个覆盖概率的小数点后两位才是可以报的。如果标准差在 0.05 以上需要增加迭代次数或改换 BPP 固定卫星数否则任何参数对比都没有意义。这个检验动作应该是参数扫描之前的固定动作而不是最后补一个“置信区间”了事。6.2 参数扫描用同一个随机种子做“可复现对比”随机几何最常见的结果图是“覆盖概率 vs 星座密度”或“覆盖概率 vs 阈值”。这类对比必须保证所有曲线在哪一轮快照上使用同一套随机数才能让曲线差异只来自参数本身。最简单做法是在函数里加seed参数并且每次调用把 seed 偏移一个固定常数。for density_factor in [0.5, 1.0, 2.0]: base_mean 200 sinr simulate_downlink_sinr( mean_sat_numbase_mean * density_factor, iterations5000, seed100 int(density_factor * 100)) print(density_factor, coverage_probability(sinr, 10.0))这样跑出来的曲线是平滑且可复现的。更高阶的做法是把不同参数组合写成循环最后用 pandas 存成 CSV。我不建议在同一个 for 循环里反复修改全局变量很容易把上一轮的星座带到下一轮。参数扫描时我习惯先看几个典型方向的趋势对照物理直觉检查参数调大对覆盖概率的典型影响排查注意卫星密度干扰也增大覆盖不一定单调上升看干扰受限区间轨道高度单链路损耗增大但可见卫星更多综合影响是非线性最低仰角可见卫星减少覆盖下降别用来“优化”结果频率复用因子同频干扰减少覆盖上升需建模波束方向至少就我自己而言这类仿真项目从“能跑出数”到“能相信这个数”中间隔的不是代码量而是验证方法。先确认链路预算没有符号错误再查可见卫星统计最后才用批次标准误判断收敛。如果在每一步都把中间量打印出来看一眼翻车的概率会小很多。希望帮到你。本文还有配套的精品资源点击获取