一种多波束通讯卫星的可见性快速预报算法

1.本发明提供一种高精度的多波束通讯卫星的可见性快速预报算法，它涉及一种适合地面终端的快速轨道预报模型，基于星下点轨迹的地面倾角的可见范围确定方法，基于分段线性插值的多波束可见性确定方法。属于卫星通讯技术领域。


背景技术：

2.低轨卫星一般指轨道高度在200km-2000km的卫星，低轨卫星主要用于地面目标侦察和地面终端通信。其中对于通信卫星，较低的轨道高度使其具有传输延迟短，路径损耗小的优点。低轨道卫星星座一般是指多个低轨卫星构成的可以合作完成特定任务的大型卫星系统，低轨道卫星星座是目前最具发展前景的卫星移动通信系统。高精度的可见性预报是卫星与地面终端建立通讯的关键技术之一。集成了高精度的可见性预报算法的地面终端在无人值守的情况下与卫星建立长期、稳定的数据传输通信协作任务，并且在非通讯时间通过待机实现低功耗长续航的优点。
3.目前对于对地成像的可见性预报已经有了很多研究，但这些研究主要集中在空间相机的成像区域建模，预报的时间精度达到秒级，对于通信卫星的分波束通信可见性预报还没有相关研究。此外，对于低轨卫星，其单波束的通讯时长只有10s-200s，为了充分利用通信资源，对于多波束通信的预报精度提出了更高的要求。由于目前的可见性预报模型并没有针对低轨卫星建立一套高精度快速轨道预报模型，限制了其可见性预报精度。
4.综上，基于低轨卫星通信系统具有巨大的研究与应用价值，本发明提出了一种多波束通讯卫星的高精度、快速可见性预报算法。


技术实现要素：

5.(一)发明目的
6.本发明首先推导了相对于地心地固系的轨道预报算法，然后介绍了基于星下点轨迹的地面倾角的可见范围确定方法，进一步基于分段线性插值的多波束可见性确定方法。目的在于提出了一种多波束通讯卫星的高精度、快速可见性预报算法，为构建智能化的卫星通讯终端构建提供技术支撑。
7.(二)技术方案
8.本发明所述的多波束通讯卫星的高精度、快速可见性预报算法实施步骤如下：
9.步骤一：准备工作
10.首先推导卫星相对于地心地固系的高精度轨道预报算法：
11.高精度的轨道预报通常是在惯性系下对空间位置速度进行积分，并且需要使用非常复杂的力学模型，包括高阶重力场模型，三体引力模型，潮汐引力模型，大气阻力，太阳光压等。由于模型误差的存在，轨道预报的精度会随着预报时间的增加而变得越来越差，这将影响过顶时间和通讯波束预报的精度。如果使用复杂的动力学模型进行轨道预报，虽然可以降低轨道预报误差，但会占用大量的计算资源。目前很多星上自主的轨道预报器使用的
是由norad开发的简化常规摄动模型(sgp4)，常用的sgp4 模型预报的位置精度在10km，对于低轨卫星的可见性预报误差达到1-2s，针对更高的预报时间精度要求，本发明利用轨道摄动cowell方法作为预报器，并给出了在地心地固系下轨道动力学方程。
12.低轨通信卫星运动微分方程：
[0013][0014]
式(1)中，r是卫星位置，f是卫星受到的中心引力和各类摄动力。
[0015]
低轨通信卫星在地心地固坐标系s
e
运动微分方程：
[0016][0017]
式(2)中，上标e表示矢量在地心地固坐标系s
e
下的分量列阵。r
e
，分别是卫星的位置，相对s
e
的速度，相对s
e
的加速度。
[0018]
分析不同类型的摄动加速度对于轨道预报的影响，确定最适合低轨卫星快速轨道预报模型：对于轨道高度大于500km的低轨卫星，可以采用9*9阶重力场模型，24 小时内的轨道预报位置误差小于1km。
[0019]
采用归一化球谐项系数描述的9
×
9阶次的地球引力势：
[0020][0021]
式(3)中，r，φ，λ分别表示地心距、地心纬度和地心经度，μ为地球引力常数， r
e
为地球赤道半径，和为归一化的n阶m次球谐项系数，为归一化n阶 m次缔合勒让德多项式，归一化系数
[0022]
根据勒让德多项式的递归关系以及三角函数的加法定理，将上述地球引力势中的递归项结合起来(cunningham)，地球引力势模型可以写为：
[0023][0024]
式(4)中，p
nm
(sinφ)cos(mλ)，p
nm
(sinφ)sin(mλ)
[0025]
低轨卫星所受的引力加速度：
[0026][0027]
由于地球引力势球谐项模型的参考坐标系为地心地固系，上述地球引力加速度实际上是卫星所受的绝对加速度在地心地固坐标系下的投影，因此用于轨道动力学数值积分的相对ecef的相对加速度可以表示成：
[0028][0029]
式(6)中，上标e表示在地心地固坐标系下的分量列阵，为绝对加速度，为牵连
加速度，为科氏加速度。
[0030]
简化的地球自转模型中不考虑极移和章动，认为地心地固坐标系是由地心赤道惯性系绕z轴旋转一定的角度(格林威治赤经)得到，因此可以写出牵连加速度和科氏加速度的计算公式：
[0031][0032]
a
c
＝2ω
e
×
v
r
＝[-2ω
e
v
y 2ω
e
v
x z]
ꢀꢀꢀꢀꢀꢀꢀ
(8)
[0033]
式(7)、(8)中，ω
e
为地球自转角速度。
[0034]
然后介绍了基于星下点轨迹地面倾角的终端临界通讯条件确定方法：
[0035]
假设地球是球形，低轨卫星是圆轨道，那么卫星的通讯范围是一个圆锥，在地球上的投影是一个圆，接收范围是指该投影圆的球面半径，在球面三角形中，可以用球心角表示。
[0036]
如说明书附图1所示：s表示测站所在位置，m表示卫星所在位置，o表示地球球心，n为卫星m的星下点。
[0037]
在三角形osm中，由正弦定理：
[0038][0039]
式(9)中，l
om
＝|r
m
|，l
ms
＝|r
s-r
m
|，l
os
＝r
e
+h，∠oms＝η。
[0040]
接收范围的半径：
[0041][0042]
式(10)中，η是本体系下满足通讯条件的最大俯仰角。
[0043]
最佳观测(通信)时刻是指卫星在终端所在的当地水平面上投影距离终端最近的一点。该点满足星下点速度方向与基线方向垂直。
[0044]
如说明书附图2所示：在地心地固坐标系oxyz中，图中n表示地球北极，s表示终端所在位置，n1表示轨道平面p1星下点轨迹上一点，o1表示n1所在纬圈平面与z轴的交点，且满足球面大圆弧n1处的地面倾角为i
e
，且满足即平面o1n1n2与赤道面的二面角为i
e
；平面o1s n1与平面o1n1n2垂直。因此n1表示图示终端的临界通讯位置。
[0045]
由上述几何关系，平面o1s n1与赤道面夹角为π/2-i
e
，得到下式方程组：
[0046][0047]
式(11)中，表示平面o1s n1的法向量，j表示z轴的单位向量。
[0048]
设s的经纬度为(λ
s
，φ
s
)，n1的经纬度为(λ1，φ1)，，由球面三角形余弦定理，可以表示成：
[0049][0050]
平面o1s n1的法向量其中：
[0051][0052]
s代表sin，c代表cos。
[0053]
代入式(11)计算可得：
[0054][0055]
联立求解式(12)中关于sinφ1的四次方程：
[0056]-s2i
e
s4φ1+2s2i
e
sφ
s
s3φ
1-(s2i
e
s2φ
s-1)s2φ
1-(2cγc2i
e
sφ
s
+2s2i
e
sφ
s
)sφ1+s2i
e
s2γ+c2γ-c2φ
s
＝0
[0057][0058]
求解后得到临界通讯位置的星下点经纬度(λ1，φ1)。
[0059]
接下来推导星下点轨迹地面倾角的计算方法。
[0060]
如说明书附图3所示：在球面三角形bmt中，其中b是升交点，m是卫星位置，t是通过m的子午面与赤道的交点。是卫星轨道平面p{ω}上的一段圆弧，角b为轨道倾角i，边是m的纬度幅角u，边是m的赤纬δ
m
，边等于α
m-ω
m
，其中m的赤经α
m
，轨道的升交点赤经ω
m
。
[0061]
由球面直角三角形的正弦定理和余弦定理：
[0062][0063]
当卫星m处在临界通讯位置，即对应着附图2中的星下点n1时，m的赤经赤纬 (α
m
，δ
m
)与星下点n1的经纬度(λ1，φ1)有如下关系：
[0064][0065]
式(14)中：α
g
(t)为当前卫星位置m对应时刻的greenwich赤经。
[0066]
在仅考虑地球扁率的长周期影响下，卫星的轨道根数仅有升交点赤经会发生变化，对于这一类轨道族p{ω}(具有相同半长轴、轨道倾角、偏心率、近地点幅角的卫星轨道)，考虑简化的地球模型，具有下列性质：(1)族中各轨道相同的赤纬处地面轨迹的倾角相等；(2)族中各轨道内从同一赤纬飞行经过相同纬度幅角所用的时间相等。
[0067]
由上述轨道族性质(2)，我们可以使用升交点赤经设计值对应的轨道来计算实际轨道的临界通讯位置n1所对应的星下点轨迹地面倾角。
[0068]
已知δ
m
＝φ1，令α
′
m-ω0＝α
m-ω
m
，代入式(14)中，得到：
[0069][0070]
考虑j2的平均摄动效果，得到：
[0071][0072]
式(16)中，w
ω
为考虑j2摄动的圆轨道升交点平均变化率：w
u
为考虑j2摄动的圆轨道纬度幅角的平均变化率：
[0073]
进一步，得到临界通讯位置n1的星下点轨迹地面倾角：
[0074][0075]
通过联立上述式(12)、(17)可以计算出终端临界通讯位置n1的经纬度(λ1，φ1)和对应的地面轨迹倾角i
e
。
[0076]
进一步，由上述轨道族的性质(2)推导基于星下点轨迹地面倾角的终端临界通讯条件，该条件用卫星穿越用终端纬线上的经度范围表示。
[0077]
如说明书附图4所示：s表示终端所在位置，m(n1)表示临界通讯位置，m2为该轨道平面p1与终端纬线的交点，n2为m2对应的星下点。虚线mm2对应的是卫星的轨道平面，实线nn1对应的是该轨道平面的星下点轨迹。其中m2，m的赤经赤纬分别为(α
m
，δ
m
)，n2，n1的经纬度分别为(λ2，φ
s
)，(λ1，φ1)。
[0078]
暂时不考虑非球形引力摄动，由式(12)得：
[0079][0080]
式(18)联立求解得：
[0081][0082]
考虑j2平均摄动影响，升交点在一个周期内的存在进动，考虑其平均变化率，对λ2进行修正，得：
[0083][0084]
由式(12)的解算结果可知终端临界通讯位置n1对应有两个点，因此n2也有两个对应点因此基于星下点轨迹地面倾角的终端临界通讯条件用卫星穿越终端纬线上的经度范围表示，对满足该条件时卫星位置继续进行可见性精确搜索。
[0085]
最后，构造分波束可见性判断函数。
[0086]
对于通讯卫星，一般认为卫星本体系相对于卫星轨道系定向，即x轴指向卫星飞行方向，z轴指向地心，y轴指向轨道角动量的反方向。其波束定义由发射天线在本体系的方位决定，可以用卫星-终端基线在本体系的俯仰角和方位角来判断终端位于某一波束的通讯范围。同时，在一个通讯窗口，卫星的位置信息具有很弱的非线性，可以采用分段插值拟合的方法表示出卫星位置的连续函数，进一步通过求解分波束可见性的判断函数得到精确地
可见性范围。
[0087]
已知卫星m在时刻t满足了基于星下点轨迹地面倾角的终端临界通讯条件，其在地心地固系中的位置为通过上述给出的轨道预报模型预报一段时间的轨道，并根据俯仰角判断条件判断是否终止轨道预报，预报的时间间隔使用10s，得到一系列的卫星位置
[0088]
对上述采用二次多项式拟合的方法，得到卫星m的位置关于时间的连续函数根据俯仰角信息构造全波束可见性判断函数：
[0089][0090]
式(20)中，r
s
是终端所在位置的地心地固系坐标。
[0091]
当f(t)＞cosη时，η是本体系下满足通讯条件的最大俯仰角。令f(t)＝cosη根为 (t1，t2)，则卫星全波束可见性范围为t∈(t1，t2)。
[0092]
进一步，根据不同波束对应的俯仰角和方位角信息构造分波束可见性判断函数。
[0093]
上述卫星-终端基线是表示在地心地固系下的分量列阵，通讯卫星的本体系相对于地心地固系的姿态转换是位置的函数，因此卫星-终端基线在本体系的分量列阵可以写成：
[0094][0095]
式(21)中，上标b表示卫星本体系，上标e表示地心地固系。
[0096]
进一步，卫星-终端基线的方位角信息与高度角信息可以表示成：
[0097][0098]
式(22)中，p(t)和a(t)分别表示俯仰角和方位角。
[0099]
分波束可见性判断函数：
[0100][0101]
式(23)中，beam
i
(η)和beam
i
(ψ)分别对应着第i个分波束的俯仰角和方位角边界。
[0102]
可以通过求解式(23)的根得到分波束可见性范围，但这种方法涉及姿态转换矩阵的求解和反三角函数的根的解算。为进一步减小计算量，可以将上述采样点代入式(22)中，得到一系列p
i
和a
i
，对新的俯仰角和方位角采样点进行拟合，由于俯仰角和方位角表现出较高的非线性，采用分段线性插值的方法，得到其连续函数和式(23)可以改写成：
[0103]
[0104]
求解式(24)的根得到通讯卫星分波束可见性范围。
[0105]
至此，多波束通讯卫星的高精度、快速可见性预报算法的准备工作已经完成。
[0106]
步骤二：低轨通信卫星的轨道参数初始化
[0107]
在当前通讯时刻，根据低轨卫星携带的导航接收机计算出卫星在地球固连系下的位置速度，并将其播发给地面终端。终端完成与低轨通讯卫星的时钟同步。
[0108]
步骤三：通讯波束、可见性预报时长、终端位置初始化
[0109]
根据通讯卫星的设计参数，完成多波束的俯仰角和方位角边界的初始化，并确定多波束可见性预报的总时长，默认设置成24小时。此外，由先验信息得到终端的经纬高。
[0110]
步骤四：计算卫星的密切轨道根数
[0111]
根据终端收到的低轨通信卫星在地心地固系中的位置，计算出卫星的密切轨道根数，进一步得到轨道周期与考虑j2摄动影响的升交点赤经进动角速度，得到在一个轨道周期卫星的星下点轨迹的西进经度。并计算基于星下点轨迹地面倾角的终端临界通讯条件。
[0112]
步骤五：预报以当前时刻为起点的第一个轨道周期内卫星位置
[0113]
使用考虑j2摄动影响的轨道预报方法，采用50s～200s的计算步长，得到一系列的卫星位置。并将卫星在地心地固系中的坐标转换成大地坐标系下星下点对应经纬度。
[0114]
步骤六：找出第一周期内最接近终端纬线的卫星位置所对应的时间
[0115]
通过步骤五种得到一系列纬度信息，找出距离终端所在纬度最近卫星位置，并得到该位置对应的时间，记为终端纬线穿越时间t0。
[0116]
步骤七：根据西进经度计算满足临界通讯条件的终端纬线穿越时间
[0117]
由步骤四、六中终端纬线穿越时间和西进经度，得到可见性预报时长内满足临界通讯条件的终端纬线穿越时刻t。
[0118][0119]
式(25)中，t0表示第一周期内终端纬线穿越时刻，t表示一个轨道周期，δλ表示当前卫星经度与临界通讯条件对应的经度的差，λ
w
表示西进经度，[x]
z
表示取整，具体是向上取整还是向下取整取决于临界通讯条件的东西边界。
[0120]
步骤八：高精度的轨道预报
[0121]
根据步骤二中的低轨卫星轨道参数，使用步骤一中推导的高精度轨道预报器，积分步长使用100s，求解出满足临界通讯条件的终端纬线穿越时刻的卫星位置速度。
[0122]
步骤九：短步长采样与判断函数拟合
[0123]
由步骤八中的卫星位置速度，采用高精度轨道预报器，采用小步长(5s-20s)，积分计算出一系列的卫星位置速度，计算卫星-终端基线在卫星本体系下的俯仰角和方位角，其中轨道预报的终止条件设置成步骤一中全波束可见性判断函数f(t)。
[0124]
进一步，通过一系列的卫星-终端基线在卫星本体系下的俯仰角和方位角进行分段线性拟合，得到步骤一中的俯仰角、方位角连续函数和
[0125]
步骤十：分波束可见性求解
[0126]
将步骤三中多通讯波束的俯仰角方位角判断条件代入到步骤一式(23)所述的方程，可以求解出分波束可见性区间。
[0127]
通过上述步骤，提出了一种多波束通讯卫星的可见性快速预报算法。该方法使用
了适合低轨卫星在轨轨道预报的精密轨道预报器，将一天内的轨道预报精度提高到分波束可见精度0.1s量级。同时采用基于星下点轨迹地面倾角的终端临界通讯条件求解方法，该方法满足了考虑j2平均摄动的精确搜索范围，使得步骤八-步骤十在更准确的终端临界通讯条件求解，避免了工程上采用过大的经验比例系数扩大终端临界通讯条件导致的无效求解过程。此外，针对通讯卫星提出了一种基于俯仰角和方位角的多波束可见性判断函数，并使用多项式拟合的方法避免了迭代求解，在不影响精度的前提下加快了计算速度。
[0128]
(三)优点
[0129]
本发明提供的提出了一种多波束通讯卫星的可见性快速预报算法的优点在于：
[0130]
①
本发明提出的算法涉及一种基于俯仰角和方位角的多波束可见性判断函数，可以针对不同通讯状态给出各波束的可见性预报结果。
[0131]
②
本发明提出的算法涉及一种地心地固系的轨道预报器，其在短期内的可见性预报达到0.1s，相比于常用的sgp4轨道预报器的可见性预报结果高1-2个量级，同时经过优化的轨道预报器在解算速度上也能适应终端的计算能力。
[0132]
③
本发明提出的算法涉及一种基于星下点轨迹地面倾角的终端临界通讯条件求解方法，避免了工程上采用过大的经验比例系数扩大终端临界通讯条件导致的无效求解过程。减少了计算量。
附图说明
[0133]
图1是本发明中接收范围的半径求解示意图。
[0134]
图2、图3、图4是基于星下点轨迹地面倾角的终端临界通讯条件求解示意图。
[0135]
图5是本发明的实施步骤流程图。
具体实施方式
[0136]
下面将结合附图5和技术方案对本发明的具体实施过程做进一步的详细说明。
[0137]
步骤一：输入参数初始化
[0138]
在当前通讯时刻，根据低轨卫星携带的导航接收机计算出卫星在地球固连系下的位置速度，并将其播发给地面终端。终端完成与低轨通讯卫星的时钟同步。
[0139]
根据通讯卫星的设计参数，完成多波束的俯仰角和方位角边界的初始化，并确定多波束可见性预报的总时长，默认设置成24小时。此外，由先验信息得到终端的经纬高。
[0140]
该步骤对应附图1中的第一个方框。
[0141]
步骤二：计算卫星的密切轨道根数与相关轨道参数
[0142]
根据终端收到的低轨通信卫星在地心地固系中的位置，计算出卫星的密切轨道根数，进一步得到轨道周期与考虑j2摄动影响的升交点赤经进动角速度，得到在一个轨道周期卫星的星下点轨迹的西进经度。并计算基于星下点轨迹地面倾角的终端临界通讯条件。
[0143]
该步骤对应附图1中的第二个方框。
[0144]
步骤三：预报第一周期终端纬线穿越时间t0[0145]
使用考虑j2摄动影响的轨道预报方法，采用50s～200s的计算步长，得到一系列的卫星位置。并将卫星在地心地固系中的坐标转换成大地坐标系下星下点对应经纬度。
[0146]
通过得到一系列纬度信息，找出距离终端所在纬度最近卫星位置，并得到该位置
对应的时间，记为终端纬线穿越时间t0。
[0147]
该步骤对应附图1中的第三个方框。
[0148]
步骤四：计算满足临界通讯条件的终端纬线穿越时间
[0149]
由步骤二、三中终端纬线穿越时间和西进经度，得到满足临界通讯条件的终端纬线穿越时刻t。
[0150]
该步骤对应附图1中的第四个方框。
[0151]
步骤五：高精度的轨道预报
[0152]
根据步骤二中的低轨卫星轨道参数，使用上述技术方案中高精度轨道预报器，积分步长使用100s，求解出满足临界通讯条件的终端纬线穿越时刻的卫星位置速度。
[0153]
该步骤对应附图1中的第五个方框。
[0154]
步骤六：基线俯仰角和方位角短步长采样
[0155]
由步骤五中的卫星位置速度，采用高精度轨道预报器，采用小步长(5s-20s)，积分计算出下一时刻的卫星位置速度，计算卫星-终端基线在卫星本体系下的俯仰角和方位角。
[0156]
该步骤对应附图1中的第六个方框。
[0157]
步骤七：精密轨道预报的终止条件
[0158]
轨道预报的终止条件设置成上述技术方案中全波束可见性判断函数f(t)＞cosη。
[0159]
该步骤对应附图1中的第七个方框。
[0160]
步骤八：判断函数拟合
[0161]
对步骤六、步骤七循环产生的卫星-终端基线在卫星本体系下的俯仰角和方位角进行分段线性拟合，得到上述技术方案中的俯仰角、方位角连续函数和
[0162]
该步骤对应附图1中的第八个方框。
[0163]
步骤九：分波束可见性求解
[0164]
将步骤一中多通讯波束的俯仰角方位角边界条件代入到上述技术方案中式(23) 所述的方程，可以求解出分波束可见性区间。
[0165]
该步骤对应附图1中的第九个方框。
[0166]
综上所述，本发明的算法流程如附图5所示。提出了一种多波束通讯卫星的可见性快速预报算法。该方法使用了适合低轨卫星在轨轨道预报的精密轨道预报器，将一天内的轨道预报精度提高到分波束可见精度0.1s量级。同时采用基于星下点轨迹地面倾角的终端临界通讯条件求解方法，该方法满足了考虑j2平均摄动的精确搜索范围，避免了工程上采用过大的经验比例系数扩大终端临界通讯条件导致的无效求解过程。最后，针对通讯卫星提出了一种基于俯仰角和方位角的多波束可见性判断函数，并使用多项式拟合的方法避免了迭代求解，在不影响精度的前提下加快了计算速度。