1997年，为方便在与同步轨道卫星业务共用的某些Ku和Ka频段内引入非对地静止固定业务系统，WRC-97通过了等效功率通量密度EPFD的概念^[1]。EPFD的定义为：非对地静止卫星系统范围内，所有发射电台在地球表面或在对地静止轨道中的对地静止卫星系统接收电台产生功率通量密度的总和。在WRC-97中，EPFD为对非静止轨道(NGSO)星座系统的发射采取的“硬”限值，即绝不能超过此限制条件^[1-2]。随着时间的推移，后续WRC大会加入了更多频段的EPFD限制^[3]。

传统计算EPFD的方法主要通过卫星轨道外推，即在某一段时间内外推星座卫星的轨道位置，计算不同时刻的EPFD，最后统计不同EPFD出现时间占比以获得最终的概率分布和累计概率分布。通过此结果判断是否会对静止轨道(GSO)星座系统产生干扰。在实际应用中，目前通常将ITU-1503建议书^[4]提供的参数输入仿真软件工具中，通过下行功率通量密度(pfd)或有效等向辐射功率(e.i.r.p)掩码计算EPFD。本质上，ITU-1503建议书中所采用的软件工具依旧使用外推方式(建议书中称为时间仿真法)，在每个时间步长上对干扰电平进行评估以获得EPFD的概率分布。

近几年，随着以Starlink为代表的巨型星座计划相继被提出，FCC授权通过的星座卫星数量越来越多：柯伊伯星座计划发射3 236颗卫星，Starlink计划发射41 943颗卫星，OneWeb星座则计划发射48 560颗卫星^[5-9]。巨型星座卫星数量远大于传统星座、分布密集，不仅增加了对原有GSO星座系统造成同频干扰的风险，也给干扰分析带来了挑战。由于GSO卫星相对地球站位置固定，且传统NGSO星座的卫星数量较少，采用外推法即可完成干扰分析的仿真工作，因此目前简化EPFD计算的研究较少。在巨型星座场景中，万颗量级的卫星将导致外推仿真时间过长，计算效率较低，甚至出现普通PC机无法完成仿真等情况^[10-12]。

文[13-18]分析了卫星出现概率、干扰场景建模等研究方法，为基于概率模型的EPFD计算方法提供了研究思路。在此基础上，本文提出了一种计算EPFD及其概率分布的方法。通过不同的星座快照计算EPFD；同时，以基准卫星的卫星出现概率表征快照概率，进而获得EPFD的概率分布。此方法不仅适用于传统卫星数量较少的NGSO星座，也可用于巨型星座的计算。仿真结果表明，本文提出的计算方法和外推法具有相同的准确度，但计算效率更高。

1 模型与方法 1.1 系统模型与EPFD计算方法

考虑NGSO星座系统与GSO星座系统同频共存的场景。图 1中，GSO系统由1个接入的GSO卫星和1个地球站组成，定义为受扰系统；NGSO星座系统为施扰系统，由S颗NGSO卫星及其地球站组成。

若某时刻GSO地球站可见s颗NGSO卫星，则这s颗可见卫星均会对GSO产生干扰。EPFD计算公式为

$ \mathrm{EPFD}=10 \lg \left[\sum\limits_{i=1}^{s} 10^{\frac{P_{i}}{10}} \frac{G_{\mathrm{t}}\left(\theta_{i}\right)}{4 {\rm{ \mathsf{ π} }} d_{i}^{2}} \frac{G_{\mathrm{r}}\left(\varphi_{i}\right)}{G_{\mathrm{r}, \max }}\right]. $

(1)

其中：P_i为NGSO卫星i的发射功率，G_t(θ_i)为NGSO卫星i星上波束在受扰地球站方向上的发射增益，θ_i为卫星i波束指向方向与受扰地球站和该卫星连线的夹角，G_r(φ_i)为受扰地球站在卫星方向上的接收增益，φ_i为GSO卫星和受扰地球站连线与受扰地球站和卫星i连线的夹角，d_i为卫星i到受扰地球站的距离，G_{r, max}为受扰地球站接收机最大接收增益。

在ITU-1503建议书中，pfd或e.i.r.p掩码定义了NGSO空间站(空间电台)和NGSO地球站辐射的功率包络，主要用于计算施扰系统卫星在不同位置的功率通量密度。在每个时间步长内，计算不同可视施扰卫星的位置向量，将其位置向量等信息作为输入参数，通过pfd或e.i.r.p掩码以获取不同施扰卫星的功率通量密度，同时，位置向量信息也可用于计算NGSO卫星的离轴角，进而获得接收增益。通过上述两部分结果计算EPFD，重复此过程，在不同时间步长内计算EPFD，可获取最终概率分布。通过上述分析可知，为获得不同施扰卫星位置向量，外推卫星轨道获取卫星位置仍是本方法的核心。

由式(1)可知，若通信参数(如增益、功率)确定，则EPFD值取决于施扰星座可视卫星数量及不同卫星与GSO通信链路的相对位置即θ_i、φ_i的大小。由于GSO位置与地球站相对位置恒定，若可确定施扰NGSO星座所有可视卫星位置，则EPFD即可计算。基于此，本文提出了一种通过概率与卫星快照的方式快速获取NGSO星座卫星位置的方法，进而可计算EPFD及其概率分布。

1.2 NGSO星座卫星位置获取方法

Walker星座因其具有良好的纬度和全球覆盖特性，是目前NGSO星座系统最常用的星座构型。具体地，Walker星座又可细分为星形星座、δ星座(玫瑰星座)等。本文所述NGSO星座均为Walker星座。

在Walker星座中，当星座中任意一颗卫星的轨道位置参数确定，便可通过此卫星参数与Walker星座构型参数推算其他所有卫星的位置，称此颗卫星为基准卫星(显然，星座内任意一颗卫星均可作为基准卫星)。δ星座中其他卫星状态推算公式如下：

$ \varOmega_{i, j}=\varOmega_{0}+i(2 {\rm{ \mathsf{ π} }} / P) ; $

(2a)

$ M_{i, j}=M_{0}+2 {\rm{ \mathsf{ π} }}\left(\frac{F}{N} i+\frac{P}{N} j\right); $

(2b)

$ i=0,1, \cdots, P-1 ; j=0,1, \cdots, \frac{N}{P}-1 . $

(2c)

其中：Ω₀代表基准卫星的升交点赤经，M₀代表基准卫星的平近点角。当星座为星形星座时，星座中其他卫星的升交点赤经推算公式如下：

$ \varOmega_{i, j}=\varOmega_{0}+i({\rm{ \mathsf{ π} }} / P), i=0,1, \cdots, P-1 . $

(3)

因此，一旦确定任意一颗卫星位置与星座构型参数，即可获得星座内所有卫星的分布位置，进而可计算EPFD。对某一时刻的卫星分布，一般称为星座快照。为计算所有可能的EPFD，需获取不同的星座快照。

由式(2)和(3)可知，Walker星座内卫星分布具有规律性：对任意时刻，任意相邻的4颗卫星围成的区域内(如图 2中的灰色区域)，在概率上有且仅存在1颗卫星。基于此，可在此区域内不同位置设置基准卫星，生成星座快照，进而可得到所有可能的星座快照。称此区域为设置基准卫星区域，简称设置区域。

图 3中，设置区域的范围由同轨道相邻卫星、相邻轨道相邻卫星组成。设同轨道相位差为F_d，相邻轨道相邻卫星的升交点赤经差为Ω_d，则设置区域范围为：

$ F_{\mathrm{d}}=2 {\rm{ \mathsf{ π} }} /(N / P) ; $

(4)

$ \varOmega_{\mathrm{d}}= \begin{cases}2 {\rm{ \mathsf{ π} }} / P, & \delta \text { 星座; } \\ {\rm{ \mathsf{ π} }} / P, & \text {星形星座. }\end{cases} $

(5)

传统外推方法是以某时长作为仿真步长，将连续的时间离散化，外推获取不同时刻点卫星位置以计算EPFD。本文方法需在设置区域内不同位置设置基准卫星，以获得所有可能的星座快照，因此同样需要对设置区域离散化处理，即在设置区域内划分子区域，划分子区域的面积应与仿真步长t精度相同。设子区域的面积为l₁×l₂，则l₁与l₂需满足：

$ l_{1}=l_{2}=\frac{t}{T} \cdot 2 {\rm{ \mathsf{ π} }}. $

(6)

其中T为卫星的轨道周期。令卫星设置区域为R_N1×N2，其中N₁为平近点角间隔数，N₂为升交点赤经间隔数。

$ N_{1}=F_{\mathrm{d}} / l_{1}, $

(7)

$ N_{2}=\varOmega_{\mathrm{d}} / l_{2}. $

(8)

对每一个子区域r_{l, k}∈R_N1×N2(l=1, 2, …, N₁; k=1, 2, …, N₂)，由于其满足仿真的精度要求且子区域r_{l, k}面积很小，将基准卫星放置在r_{l, k}的任何位置均可认为产生的EPFD相近。为方便计算，后文中统一将基准卫星放置在r_{l, k}中心。

联立式(4)—(8)可得：

$ N_{1}=\frac{T}{t(N / P)}, $

(9)

$ N_{2}= \begin{cases}\frac{T}{t P}, & \delta \text { 星座; } \\ \frac{T}{2 t P}, & \text {星形星座. }\end{cases} $

(10)

随着星座规模的扩大，同轨道卫星数N/P、轨道面数P将大幅增加，式(9)和(10)中N₁和N₂将减小，导致子区域r_{l, k}数量即N₁N₂减小。

在巨型星座场景下，卫星数量很多，则4颗相邻卫星围成的设置区域面积很小，较小的设置区域将导致设置的基准卫星较少，进而导致快照数量较少。在后文计算中，较少的快照数将导致较少的计算次数。因此，此方法在巨型星座场景下可有效减少计算次数，提高计算效率。

1.3 EPFD概率计算方法

节1.2给出了获取星座中不同卫星位置的方法，进而可计算EPFD。然而，为计算其概率分布，还需确定不同EPFD的概率。由于EPFD唯一取决于卫星位置，即EPFD概率唯一取决于快照概率，而一旦基准卫星确定，星座快照即可确定，因此，快照概率又唯一取决于基准卫星的出现概率。综上所述，若确定了基准卫星的出现概率，则可获得EPFD的概率。

文[19]中给出卫星出现概率的计算公式。当星座卫星总数为n，对某一空域D，在此空域中星座内任意一颗卫星的出现概率为

$ p=n \frac{A}{2 {\rm{ \mathsf{ π} }}^{2} \sin \alpha \cos L}. $

(11)

其中：A为空域D的球面积(sterad)，单位为rad²；α为卫星轨迹与D中心纬度线的夹角；L为计算空域D中心的纬度。

$ A=\Delta \theta_{\varepsilon} \cdot \Delta \theta_{\beta}, $

(12)

$ \alpha=\arccos \frac{\cos \gamma}{\cos L}. $

(13)

其中：Δθ_ε、Δθ_β为空域D边长对应的地心角度，由于本文中子区域均为矩形，因此式(12)仅给出矩形面积的计算方法；γ为卫星轨道的倾角。

在设置区域内，不同子区域中均放置一颗基准卫星，由于子区域很小，子区域中卫星的出现概率可等价为子区域中基准卫星的出现概率。式(6)给出了子区域的球面积，结合不同子区域所在纬度位置，代入式(12)中即可获得不同基准卫星的出现概率。由前述分析，基准卫星的出现概率一旦确定，则EPFD的概率即可确定。

1.4 EPFD概率分布计算方法

如图 4所示，在设置区域内依据仿真步长精度需求划分子区域，得到设置区域R_N1×N2，依节1.3方法计算得到不同子区域的概率，记为P_N1×N2。若选取基准卫星RS_{l, k, 1}(l=1, 2, …, N₁; k=1, 2, …, N₂)，其卫星出现概率为p_{l, k}，生成星座快照，则快照概率为p_{l, k}。此时，施扰NGSO系统中所有卫星位置已知，即可通过式(1)计算EPFD。设此时产生的EPFD为e_{l, k}，则其对应的概率为p_{l, k}。

遍历计算N₁N₂个子区域中基准卫星生成的星座快照产生的EPFD，得到EPFD矩阵E_N1×N2，其对应的概率矩阵为P_N1×N2。为得到最终概率分布，需对计算结果进一步处理：将不同EPFD归并至相应EPFD区间，不同EPFD所对应的概率将被累加为当前EPFD区间的概率值，设以x dB为区间粒度，统计EPFD在区间[e-x, e)的概率，其中e为统计的EPFD值，则统计方法为

$ p(e)=\sum\limits_{l=1}^{N_{1}} \sum\limits_{k=1}^{N_{2}} p\left(e_{l, k}\right), e-x \leqslant e_{l, k}<e . $

(14)

对所有干扰值与概率值重复此过程，最终得到EPFD的概率分布。

2 仿真分析

为对比本文方法与实际外推的效果，选用表 1和2中的参数完成相应的仿真，分析NGSO星座对GSO星座下行EPFD的概率分布情况。为方便计算，采用相同通信频点相同带宽计算干扰最严重的情况。GSO卫星的天线方向图为ITU-R S.465-5；NGSO星座卫星使用固定波束，采用对地定向方式，天线方向图为ITU-R S.1528。

表 1中NGSO星座采用了传统星座构型，δ为星座，卫星数量较少。仿真结果如图 5和6所示。可以看出，本文方法的仿真结果在概率分布、累积概率分布曲线上，与传统外推方法结果一致。图 3选取了部分EPFD值及其概率，外推结果与本文方法的误差最大在0.03%左右。

图 7给出了仿真1、30 d与本文方法的对比结果。需要说明，业界的仿真模拟时长通常为数月，以获取较为准确的仿真结果，本文最长仿真模拟时长选取一个月以对比分析。可以看出，仿真模拟时长为1 d时，曲线抖动较大，且最大EPFD(放大区域)无法仿真到。在仿真时长上，外推方法需要325 s，而本文方法需要433 s，在卫星数量较少的场景，本文方法和外推方法的仿真时长相差不大。

表 4为一组NGSO巨型星座构型参数，为星形星座，卫星数量很多。仿真结果如图 8和9所示，可以看出，在概率分布曲线上，本文方法存在一定抖动，而在累计概率分布曲线上与外推方法几乎一致。概率分布曲线存在抖动的主要原因是：在巨型星座场景下，相邻4颗卫星距离较近，因此设置区域面积较小，导致通过式(8)和(9)计算得到的子区域数量较少(百量级)，因此计算结果存在一定的抖动，但误差很小(最大0.015)，且对累积概率分布曲线几乎没有影响。此时，若减小仿真步长/粒度，即增加子区域数(计算点)，可进一步降低由于计算点较少导致的抖动误差。图 10展示了本文方法在不同仿真步长/粒度下干扰概率分布曲线的变化，随着步长/粒度不断缩小，本文方法曲线的抖动也逐渐减小。

需要说明，由于计算机内存限制，巨型星座场景下，本仿真使用的PC机仅可仿真处理外推1 d的情况，最终外推30 d仿真结果由多次仿真的合并处理得到。具体地，本文在合并处理时，将上一次仿真的场景终点作为下一次仿真的场景起点，即在上一次仿真结束时，保留当前的时间点、轨道位置、波束指向等一系列仿真相关的参数，作为下一次仿真起点的初始参数，以保证合并后仿真的连续性。本场景仿真30 d，即重复此过程30次，因此30次仿真完整覆盖了30 d的场景变化，具有完备性。最终，将所有仿真结果合并处理后可获得与一次性仿真30 d相同的效果。

本仿真中，外推法仿真时长为13 019 s，而本文方法仿真时长为325 s，显而易见，本文方法在巨型星座场景下仿真时间缩短了至少一个数量级，计算效率大幅提升。与δ星座的结果对比可知，本文方法的仿真时间不随星座规模扩大而增加；相反，由于巨型星座卫星更密，设置区域面积变小，本文方法在处理巨型星座场景时具有更高效率。

综上所述，本文方法不论在传统星座或者巨型星座，均可获得较为准确的EPFD概率分布曲线。在传统星座场景中，由于卫星数量较少，本文方法计算效率与外推法处于同一数量级，相差不大；随着卫星数量的增加，外推法计算时长大幅上升，而本文方法可以快速获得EPFD的概率分布。

3 结论

本文提出了一种计算EPFD概率分布的方法。通过在设置区域内不同位置放置基准卫星生成星座快照，获取所有可能的卫星分布，以此计算所有可能的EPFD；同时，通过基准卫星的卫星出现概率表征快照概率，以获取不同EPFD的概率，进而可获得EPFD的概率分布。仿真结果表明，本文方法可同时适用于传统星座和巨型星座，不仅与传统方法具有相同准确性，且计算效率更高，尤其在巨型星座场景中计算效率提高了一个数量级，可用于快速评估星座间的EPFD。