非静止轨道卫星星座系统功率通量密度包络的计算方法

张晓燕 刘畅 李明明 李伟

（国家无线电监测中心，北京 100037）

引 言

近年来，随着公众对高清电视、交互式多媒体服务和宽带互联网接入等新兴应用的需求增加，不仅促使5G通信进入了商用，而且也使非静止轨道(non-geo stationary orbit, NGSO)卫星互联网星座系统进入了又一个迅猛发展的浪潮.这是由于卫星互联网星座系统采用了较高频段(Ku、Ka和Q/V频段)、中低轨道、星座组网、卫星模块化设计、频率复用、多波束天线、一箭多星等技术[1-3]，可以满足消费者以较高性价比随时随地无缝访问电信服务的愿望.典型的NGSO卫星互联网星座系统包括O3b、OneWeb和Starlink等.卫星互联网星座的飞速发展必然会为已有的无线电系统，尤其是同频或邻频电台带来额外的频率干扰问题.新兴的卫星星座系统动辄由成百上千颗NGSO卫星编队组成，在数字通信有效载荷、高级调制方案、频率复用和多波束天线等方面存在巨大差别，除了卫星设计操作或运营者外，其他人很难得到计算星座系统辐射能量的详细参数或方案，但开展NGSO卫星星座系统与其他无线电系统间的干扰协调以及对系统的最大辐射能量进行监管约束必不可少[4-6].因此本文主要解决如何计算NGSO卫星星座系统下行链路的最大辐射能量——功率通量密度(power-flux-density，PFD)包络问题.

NGSO卫星星座系统PFD包络的定义为从地球表面上任意一点观察到的，干扰NGSO系统中任意空间无线电台产生的最大集总PFD值[4].即NGSO卫星星座系统空间无线电台在其生存周期的不同阶段，不论使用哪种资源分配和切换策略等，辐射功率都能包括在计算出的PFD包络内，且又可满足不对NGSO系统施加太多的限制，同时能更好地适应各种系统设计，从而实现对PFD包络的恰当估计，兼顾计算效率，以达到减少不必要的协调和高效利用频谱效率的目的.PFD包络计算方法可参见文献[7-11]等，而卫星互联网星座系统的PFD包络计算均是参考文献[11]，然后由用户各自计算.但文献[11]中只给出了提交给国际电联的PFD包络的格式要求、定义与基本框架描述，并没有给出如何具体计算PFD包络的方法.因此，经过国内外调研，兼顾我国空间业务的管理现状，本文在文献[11]的计算框架上进行了改进和细化，如图1所示，主要包含将NGSO卫星天线实际发射增益加入计算，代替现有计算中的最大天线增益，明确最小隔离角的计算和波束配置与覆盖的数学表达，以及简化和完善计算步骤等.

图1 NGSO卫星星座系统PFD包络计算框图Fig.1 Calculation block of the PFD mask radiated by NGSO constellation system

改进计算框图1的计算步骤如下：

1)根据输入参数、干扰规避和减缓策略以及卫星建链的限制条件，计算所需轨道参数(详见1.1节)，并对星座系统对应轨道面的NGSO卫星进行可用波束配置和覆盖分析，得出对应的GSO可见弧范围.在此步骤中明确了波束配置与覆盖的数学表达等内容，详见1.2节.

2)在上述基础上，计算覆盖区域内每点的最小隔离角以及对应的最小经度差.在此步骤中明确了最小隔离角的计算等内容，详见1.3节.

3)计算上述隔离角到达地面的偏轴角及其对应的波束的卫星发射增益.此步骤主要用于计算NGSO卫星天线实际发射增益，详见1.2节.

4)计算覆盖区域内每点的同频波束集总PFD，详见1.4节.

5)按照步骤1～4循环遍历星座所有轨道面，计算轨道面覆盖纬度范围内所有点的集总PFD值，得出星座的PFD包络：PFD{(αminij，△Lonij)}，详见第2节.

1 算 法

对于NGSO卫星星座来说，相同轨道面上NGSO卫星的覆盖特性大致相同，可认为同一轨道面上的卫星具有相同的PFD包络，因此对任意轨道面的一颗NGSO卫星进行波束配置和覆盖分析后，按照下文步骤计算该NGSO卫星的集总PFD包络.由于PFD包络与NGSO卫星运行时刻无关，因此文中只需遍历星座系统中的所有轨道面和星座覆盖范围的几何关系，就可得出NGSO卫星星座的PFD包络.

1.1 NGSO轨道参数

国际电联无线电通信部门(International Telecommunications Union-Radio communications Sector，ITU-R)公布的NGSO卫星参数，包括星座中卫星的数目、轨道倾角i、周期T、远地点高度ra、近地点高度rp、升交点赤经 Ω、近地点幅角ω，卫星发射信号时离地最低高度h0、初始相位φ0、同频波束中心点间距离df.则轨道偏心率e=(ra−rp)/(2a)，长半轴a=RE+(ra+rp)/2，其中地球半径RE=6 378 km.

通过不同时刻的平近点角

式中，NGSO卫星重力常数µ≈398 600 km3/s2.按照若ME<π时E的初始值取ME+e/2，ME≥π时E的初始值取ME−e/2，由牛顿法迭代求开普勒方程[12-14]

直到达到所要求的精度范围，从而求出各时刻的偏近点角E.由偏近点角E求出真近点角

和NGSO到地心距离

再根据近地点幅角ω和初始相位φ0，可求出t=0时刻的真近点角θ0：

若NGSO卫星轨道为圆轨道，则θ=E=ME.接着通过球面余弦定理，求出NGSO卫星的星下点纬度L atS：

又有NGSO卫星在地心赤道坐标系中的坐标

则

卫星发射离地高度

NGSO卫星的可见GSO轨位弧范围为

式 中： ∆Lon=LonG−LonS， LonG和 LonS分 别 为GSO可见轨位弧采样点经度和NGSO卫星的经度；cosx1=RE/RG，RG=42 162 km为GSO轨道半径；cosx2=RE/RS.同时需兼顾NGSO卫星覆盖区到GSO可见轨位弧的仰角不小于5°，即两者间的距离范围为[35 786，41 124.624] km以及h≥h0等建链限制条件.

1.2 波束配置和覆盖

由诸如同频波束中心点间距离、频率及其复用方式、极化、覆盖区、发射功率、天线辐射方向图、干扰规避和减缓策略等确定NGSO卫星的波束配置和覆盖，如图2所示.栅格区为某时刻NGSO卫星某波束的覆盖区，栅格区的中心为该波束的中心点C、栅格的任意采样点为P.点C坐标为( LatC，LonC) ，点P坐标为( L atP， LonP)，即波束覆盖范围中的第i行第j列采样点，i和j为自然数.

为了得到星下点D、波束中心C与采样点P(用户)间的数学表达式，我们使用卫星体坐标系表示.如图2所示，即以NGSO卫星所在位置S为原点，S到星下点D的方向矢量为卫星体坐标系的YS轴，S指向覆盖范围内任意点P的方向矢量在XSZS面的投影间的夹角为相对于卫星的仰角L∈[0，π/2]，与YS轴正向的夹角为相对于卫星的方位角A∈[0，2π).与YS轴间的夹角为φ，在YS轴的垂足为M，则∠NMC为P相对于星下点的方位角θ.于是有：

对于±800 kV滇西北至广东特高压直流输电线路工程共塔段线路，若将接地极线的铁塔横担更换为复合横担，复合横担长度暂取4.2 m(暂取该工程各塔型铁塔横担长度中的最小值)，雷击闪络率如表2所示。

图2 波束配置与栅格采样示意图Fig.2 Schematic of beam configuration and grid sampling

波束覆盖一般用圆形表示，波束中心相对于卫星的位置可确定，即(θC，φC)已知，那么该波束天线的指向就确定了.

利用上述推导出的各点相对于卫星的仰角和方位角，可得出卫星、波束中心和用户位置的相对关系，便于进行后续的遍历和计算.即，波束中心的极坐标为三维直角坐标为[cosLCcosAC， sinLCcosAC， cosLC]T；用户的极坐标为A，π/2−L]T，三维直角坐标为[cosLcosA，sinLcosA，cosL]T.间的夹角β(即用户到波束中心的半视角)是该波束的偏轴角.

将β代入波束天线辐射方向图就可以求出NGSO卫星天线在位置P的实际发射增益.同理，可求出其他同频波束到点P的实际发射增益.

1.3 隔离角

为了更方便地计算出NGSO星座系统下行链路对GSO卫星系统的影响，引入NGSO卫星波束与GSO轨位弧间隔离角α(如图3)，即NGSO与GSO轨位弧在地面所成的隔离角.隔离角越小意味着NGSO卫星辐射到地面进入GSO卫星接收地球站的能量越多，GSO卫星下行链路更容易受到NGSO卫星的影响.

图3 隔离角几何示意图Fig.3 Geometry of the separation angle

图3中GSO轨位弧上的黑点为采样点，下面求隔离角.遍历式(8)的GSO可见轨位弧范围，设任一采样点的地心赤道三维坐标为RG[ cos(L onG) ，sin(L onG)，0]T.从图3中可知隔离角为

为了求出隔离角，将图2中以NGSO卫星为原点的坐标系XSYSZS转到地心赤道坐标系XYZ，通过式(16)进行变换，即先将NGSO卫星原点平移到地心，再进行坐标轴旋转，则任意用户E在地心赤道坐标中的三维坐标为

式 中：a11=−cos(LonS)；a12=−sin(LonS)；a13=0；a21=cos(LatS)sin(LonS)+sin(LatS)；a22=−cos(LatS)cos(LonS)；a23=−sin(LatS)；a31=−sin(LatS)sin(LonS)+cos(LatS)；a32=sin(LatS)cos(LonS)；a33=cos(LatS).

式中， LonG=LonS+k∆step，k=1,2,3,···， ∆step为可见轨位弧的采样间隔.将式(17)、(18)代入式(15)可求出最小隔离角αmin=min{α}.遍历所有波束中心点，计算对应的最小隔离角αmin，若最小隔离角小于对应位置的GSO卫星排他角，则相当于关闭此波束，后续计算集总PFD包络时，此波束发射功率为0.

1.4 集总PFD包络

NGSO卫星星座中每轨道面上卫星到达地面的PFD为

对每个纬度的所有波束覆盖区域进行栅格采样，计算每个采样点P的(αmin，△Lon)以及叠加计算到达该采样点的同频波束集总PFD.最后找出该纬度PFD序列中具有相同αmin的PFD值的最大值作为该纬度的PFD包络值.按照上述方法，对NGSO卫星星座的所有轨道面及其通信覆盖范围进行遍历，得出该NGSO卫星星座的集总PFD包络.

2 仿真计算

以卫星网络资料3ECOM-1进行仿真计算，如图4和图5所示.

图4 3ECOM-1卫星星座Fig.4 Satellite constellation of 3ECOM-1

图5 3ECOM-1卫星星座仿真配置示意图Fig.5 The simulation setup of 3ECOM-1 satellite constellation

2.1 仿真参数

计算该卫星星座PFD包络的主要输入参数如表1所示.

表1 仿真输入参数Tab.1 The simulation input parameters

2.2 仿真步骤

1)用仿真输入参数根据1.1节构造出图4的星座构型.

2)遍历星座中的每个轨道面.

3)针对每个轨道面，计算每个卫星辐射频率组的PFD值.将频率组中频率重叠的PFD包络值进行合并处理，即相同频段的PFD包络值取不同频率组中计算出的PFD的最大值.

4)遍历卫星网络中每个频率组，获取当前频率组的最低频率以及最高频率，确定参考带宽.参考带宽是由《无线电规则》第22条相应频段的限值参数确定的.

5)计算每个辐射频率组的PFD值.对该轨道的波束覆盖纬度范围，按照二分法逐级比对计算出的PFD值，直到两个步进计算出的PFD值很接近为止.然后选取PFD值很接近时，步进较大的那个作为纬度遍历的最终步进值.

①根据1.1节求解NGSO卫星位置

由得出的纬度步进和已知的卫星星座纬度范围以及轨道参数，可求出NGSO卫星的坐标位置.

②确定波束配置和覆盖

根据1.2节配置NGSO卫星波束参数，如图2所示，配置该卫星星座的波束形状、相对位置关系、频率复用情况、采用的天线辐射方向图、每个频率组对应的天线最大增益和发射功率、df、排他角等.注意将那些波束中心处的最小隔离角小于对应位置排他角的波束进行关闭处理，即该波束不参与计算.

③覆盖区栅格化

对覆盖区进行栅格化遍历，经度和纬度步进均以0.1°进行，由此得到一个点集：

将上述栅格化点集按照1.3节方法转化为对应的最小隔离角和经度差点集：

其中，GSO可见轨位弧按照步长为10−5rad的步进，进行遍历计算.

④计算各栅格点的集总PFD

按照1.4节计算式(20)中每点的集总PFD.首先根据栅格点和波束中心点的位置，计算各同频波束到达此栅格点的偏轴角，再根据对应波束的天线辐射方向图得出对应的波束到达此点的增益，然后结合对应波束的发射功率代入式(19)中计算出该栅格点的集总PFD.其他各点的PFD值同理计算.

注意：本次计算的星座中卫星的发射天线辐射方向图为ITU-R空间电台天线的标准辐射方向图即各个方向均匀辐射，天线增益为各频段对应的最大天线增益.另外，因未找到卫星发射天线的交叉极化覆盖范围及其天线赋形、交叉极化隔离度等相关信息，本次仿真计算中未考虑同频天线交叉极化的影响，即设交叉极化隔离度为无穷大.

⑤得出每个频率组的PFD包络

得出每个栅格点的PFD后，找出式(20)中具有相同最小隔离角的栅格点中PFD的最大值及其对应的(αminij，△Lonij)作为对应频率组下，此纬度处该栅格点处的PFD包络值.

遍历所有轨道后，可求出该卫星星座的PFD包络集，仿真结果见2.3节.

2.3 仿真结果

由上述仿真计算，可得出3ECOM-1卫星星座的PFD包络，结果如表2所示.

表2 PFD包络计算结果Tab.2 The mask of PFD

频率/GHz Lats/(°) αmin/(°) △Lon/(°) PFD(dBW·40 kHz)17.7～18.6 ±60−90.0 −80 −138−90.0 80 −138−39.0 −80 −138−30.0 −80 −138−0.5 −80 −168−0.5 80 −168 0.5 −80 −168 0.5 80 −168 30.0 −80 −138 30.0 80 −138 90.0 −80 −138 90.0 80 −138 19.7～20.2 ±60−90.0 −80 −151.5−90.0 80 −151.5−25.0 −80 −151.5−25.0 −80 −151.5−0.5 −80 −168.0−0.5 80 −168.0 0.5 −80 −168.0 0.5 80 −168.0 25.0 −80 −151.5 25.0 80 −151.5 90.0 −80 −151.5 90.0 80 −151.5

由于无法真实地得到3ECOM-1卫星星座的详细系统参数和运作参数，这些参数以及干扰减缓和规避策略等往往是卫星操作和运营者的商业秘密，很难获得.因此本文的仿真参数仅能依据国际电联国际频率信息周报(international frequency information circular, IFIC)第2 920期中的SRS数据库，详见仿真输入参数表1，而无法得知该卫星星座真实的波束复用、干扰规避方案等.本文的仿真结果可能会高于3ECOM-1卫星星座系统的实际包络值，但本文的计算方法为一种通用的方法，对所要仿真的系统的各项参数了解得越详细，得到的仿真结果就会越接近真实值.

3 结 论

本文的方法为计算NGSO卫星网络的PFD包络的一种通用方法.为了使仿真结果更接近NGSO卫星系统，尤其是星座系统的真实值，需要详细地了解所要仿真系统的各项参数，尤其是波束配置和覆盖等星座运作参数.本文计算出的PFD包络满足国际电联PFD掩码格式，可用于参与后续下行等效PFD的计算，从而用于判断NGSO卫星星座是否满足《国际电联》第22条的限值以及开展NGSO卫星星座系统与其他无线电系统间的干扰协调和监管约束等.

本文主要通过计算用户到波束中心的夹角β，将NGSO卫星天线实际发射增益加入计算PFD包络中，从而代替现有国际电联ITU-R S.1503-3计算方法中采用最大天线增益的全向天线参与计算的模式；补充明确用户点处最小隔离角的计算和波束配置与覆盖的数学表达；明确覆盖范围栅格遍历及其PFD包络的仿真计算步骤，以达到更精确估计PFD包络且兼顾计算效率的目的.以本文为基础，后续会开展NGSO卫星星座下行链路等效PFD的研究.