CN109946728B

CN109946728B - 一种适用于卫星用户站数字跟踪接收机的程序跟踪方法

Info

Publication number: CN109946728B
Application number: CN201910188102.3A
Authority: CN
Inventors: 孙永磊; 邹世源
Original assignee: Beijing Institute of Remote Sensing Equipment
Current assignee: Beijing Institute of Remote Sensing Equipment
Priority date: 2019-03-13
Filing date: 2019-03-13
Publication date: 2021-10-26
Anticipated expiration: 2039-03-13
Also published as: CN109946728A

Abstract

本发明公开了一种适用于卫星用户面数字跟踪接收机的程序跟踪方法，包括星历轨道根数到WGS‑84的变换、用户站定姿定位、天线执行角计算。所述包括星历轨道根数到WGS‑84的变换是将星历预报模型给出的星历轨道根数转换到WGS‑84坐标系下的表示；用户站定姿定位是利用载体的姿态测量组合计算出载体在WGS‑84坐标系下的表示；最后将天线方向矢量由WGS‑84坐标系转换到载体坐标系，计算出天线执行角度或波束切换角度。本发明通过卫星对准和跟踪模块能精确跟踪卫星运动轨迹，具有可靠，有效，实用的特点，功能稳定，效果优良。

Description

一种适用于卫星用户站数字跟踪接收机的程序跟踪方法

技术领域

本发明涉及一种卫星跟踪技术领域，特别是一种适用于低轨卫星地面用户站数字跟踪机的程序跟踪方法。

背景技术

在移动互联网时代，高速率、低时延、广域无差别的互联网接入能力是迫切需求。在

卫星通信更容易发挥其广域覆盖的优势，尤其是低轨卫星通信，其覆盖性不受地形约束，其通信时延与地面光纤网络处于相同的数量级，其通信系统容量可以通过增加卫星数量来提升。在卫星技术快速发展的今天，低轨小卫星星座以其成本和能力优势能满足建设偏远地区空间信息走廊的需求。

结合我国地区开发的战略部署，提出了低轨卫星商业航天项目，以提供随时随地、基于“热点”模式的互联网接入能力的目的，基于低轨卫星星座以及地面系统构建天地一体化信息系统。

系统预计发射154颗低轨卫星，由于处于前期技术验证阶段，先发射一颗技术验证星，且卫星绕地球运转一周时间约为100分钟，相对于地球自转速度较快，所以地面用户站需要设计跟踪接收机，使天线波束时刻跟踪对准卫星，保证地面用户站接收到的卫星信号强度最大。

发明内容

为了解决上述技术问题，本发明提出了一种适用于卫星用户站数字跟踪接收机的程序跟踪方法，该方法结合卫星运行轨道模型和用户站当前姿态角度，实时计算出天线执行角度或者波束切换方向。

对此本发明提出一种适用于卫星用户站数字跟踪接收机的程序跟踪方法，其特征在于，所述方法的步骤如下：第一步，将星历轨道根数(a,e,i,Ω,ω,θ)转换到WGS-84坐标系下的表示；第二步，将来自姿态测量组合的经、纬、高

转换到WGS-84坐标系下的表示；第三步，将站星方向矢量转换为n系下的表示；第四步，对用户站进行定姿态，即将用户站由n系转换到b系；第五步，计算天线执行角。

其中，所述第一步，将星历轨道根数(a,e,i,Ω,ω,θ)转换到WGS-84坐标系下的表示，包括如下转换步骤：

(1)、星历轨道根数到J2000惯性坐标系的转换

f)

g)

h)第一向量

i)第二向量

j)卫星在J2000坐标系下的位置向量

其中星历轨道根数(a,e,i,Ω,ω,θ)定义如下：

a为半长轴；e为偏心率；i为轨道倾角；Ω为升交点赤经；ω为近地点幅角；θ为真近点角，

μ为地球引力常量，其值为3.986004418×10¹⁴；

(2)、J2000惯性坐标系到瞬时平赤道坐标系MOD的转换

g)将协调世界时转换成儒略日：

其中fix()为取整函数；

Year,Month,Day,Hour,Minute,Second为协调世界时UTC时间的年，月，日，时，分，秒；

h)太阳系质心动力学时：

其中LeapSecond为跳秒数，值为常数37；

i)由太阳系质心动力学时对应的儒略儒略世纪数：

j)赤道面进动欧拉角：

ζ_A＝2306.2181*Δτ+0.30188*Δτ²+0.017998*Δτ³ (9)

Z_A＝2306.2181*Δτ+1.09468*Δτ²+0.018203*Δτ³ (10)

θ_A＝2004.3109*Δτ-0.42665*Δτ²-0.041833*Δτ³ (11)

k)转换矩阵：

A＝M₃(-Z_A)M₂(θ_A)M₃(-ζ_A) (12)

其中

其中:M₃(x)，M₂(x)表示函数，x为函数自变量；

l)位置向量

(3)、瞬时平赤道坐标系MOD到瞬时真赤道坐标系TOD的转换

d)平黄赤交角：

e)章动转换矩阵：

N＝M₁(-ε-Δε)M₃(-Δψ)M₁(ε) (17)

其中

f)位置向量

(4)、瞬时真赤道坐标系TOD到准固定地球坐标系ECPF的转换

g)世界时：

UT1＝UTC+UT1UTC+DUT1 (20)

其中：UT1UTC为当前转换时刻的UTC与对应的世界时UT1之间的差值，其值可通过EopDE405星历表查到，DUT1为UT1UTC的误差，定为0；

h)由转换时刻的UTC时间对应的世界时UT1，从2000年1月1日12时起算的儒略世纪数：

i)格林尼治平恒星时角：

j)格林尼治真恒星时角：

其中：

为步骤4中c中的格林尼治平恒星时角，Δψ为黄经章动；Δε为交角章动；其值可通过NutationDE405星历表查到，ε为步骤(3)中求得的平黄赤交角；

k)转换矩阵：B＝M₃(θ_g) (24)

l)位置向量：

(5)、准固定地球坐标系ECPF到WGS-84坐标系的转换

c)极移修正转换矩阵：

其中：X_p、Y_p为地球自转极移角；

d)位置向量

其中：

为卫星在WGS-84坐标系下的位置向量，W为极移修正转换矩阵，B为步骤4中的转换矩阵，N为步骤3中的章动转换矩阵，A为步骤2中的转换矩阵，

为步骤1中卫星在J2000坐标系下的位置向量。

其中，所述第二步，将来自姿态测量组合的经、纬、高

转换到WGS-84坐标系下的表示，其转换步骤如下：

c)卯酉圈曲率半径：

d)位置向量

其中a为地球半长轴，其数值为6378137.0米，b为地球半短轴，其数值为6356752.3142米。

其中，所述第三步，将站星方向矢量转换为n系下的表示，实施方式如下：

(1)、首先将站星方向矢量转换为WGS-84坐标系下的表示，即将星历轨道根数在WGS-84坐标系下的表示

减去姿态测量组合的经、纬、高在WGS-84坐标系下的表示

即：

其中：

为站星方向在WGS-84坐标系下的表示，

分别为向量

对应的三个坐标表示。

(2)、将站星方向矢量在WGS-84坐标系下的表示，转换到n系下的表示：

d)转换矩阵：

e)

其中

为站星方向在n系下的三个坐标表示；

f)对

进行归一化处理。

其中，所述第四步，对用户站进行定姿态，将用户站由n系转换到b系，具体实施方式如下：

d)由姿态测量组合给出的旋转角测量值

构建由n系向b系的旋转变换矩阵

e)转换矩阵

f)坐标转换：

其中

为用户站在b系下的三个坐标表示。

其中，所述第五步，计算天线执行角如下：

天线执行角

在b系的含义如下：

为离轴角，其是归一化天线方向矢量

与oz_b轴之间的夹角，范围0°-57°；

为旋转角，其是归一化天线方向矢量

在x_boy_b平面内的投影

与ox_b轴之间的夹角，顺时针方向为正；

天线执行角与归一化天线方向矢量之间的换算关系如下：

其中，载体坐标系ox_by_bz_b b系中，原点o位于天线中心位置；ox_b轴指向载体正前方；oy_b轴指向载体正左方；oz_b轴指向载体正上方。

本发明采用了上述技术方案，由于卫星相对于地球运动速度较快，对信号处理的速度和精度要求较高，在工程实现过程中采用FPGA或者专用DSP芯片对上述计算过程进行运算，提高了数据处理速率，本发明具有可靠，有效，实用的特点，功能稳定，效果优良。

附图说明

图1为本发明的一个用户站整机系统结构示意图；

图2为本发明的一个数字跟踪接收机系统结构示意图；

图3为本发明的一个程序跟踪模型示意图；

图4为本发明的一个程序跟踪方法示意图。

具体实施方式

下面结合附图说明根据本发明的具体实施方式。数字跟踪接收机中程序跟踪模型如附图3所示：

在程序跟踪状态下，根据跟踪预报模型给出的对星角

和姿态测量组合给出的由n系向b系的旋转变换矩阵

通过计算确定天线执行角

实现将天线波束对准卫星的功能。

其实现步骤如下：

a)根据预报模型给出的对星角

计算归一化方向矢量

作为n系基准方向矢量；

b)根据姿态测量组合给出的旋转角测量值

构建由n系向b系的旋转变换矩阵

得到b系预计方向矢量

c)由式(34)计算天线预计执行角

d)相控阵天线对波束方向的控制采用离散化的方式，因此需要将天线预计执行角

对应至临近的可实现波束控制方向，对应的天线执行角为

e)调转天线波束至指定位置；

f)等待下一周期的对星角

和旋转角测量值

重复a)-e)。

实施例

本发明一种适用于卫星地面数字跟踪接收机的程序跟踪方法，其实现步骤如图4所示，具体实施例的步骤如下：

第一步，将星历轨道根数(a,e,i,Ω,ω,θ)转换到WGS-84坐标系下的表示；

第二步，将来自姿态测量组合的经、纬、高

转换到WGS-84坐标系下的表示；

第三步，将站星方向矢量转换为n系下的表示；

第四步，对用户站进行定姿态，即将用户站由n系转换到b系；

第五步，计算天线执行角。

下面对各个步骤的具体实施方式做出详细说明。

第一步中，将星历轨道根数(a,e,i,Ω,ω,θ)转换到WGS-84坐标系下的表示，包括如下转换步骤：

(1)、星历轨道根数到J2000惯性坐标系的转换

k)

l)

m)第一向量

n)第二向量

o)卫星在J2000坐标系下的位置向量

其中星历轨道根数(a,e,i,Ω,ω,θ)定义如下：

a为半长轴；e为偏心率；i为轨道倾角；Ω为升交点赤经；ω为近地点幅角；θ为真近点角，μ为地球引力常量，其值为3.986004418×10¹⁴。

(2)、J2000惯性坐标系到瞬时平赤道坐标系MOD的转换

m)将协调世界时转换成儒略日：

其中fix()为取整函数；

n)太阳系质心动力学时：

其中LeapSecond为跳秒数，值为常数37；

o)由太阳系质心动力学时对应的儒略儒略世纪数：

p)赤道面进动欧拉角：

ζ_A＝2306.2181*Δτ+0.30188*Δτ²+0.017998*Δτ³ (9)

Z_A＝2306.2181*Δτ+1.09468*Δτ²+0.018203*Δτ³ (10)

θ_A＝2004.3109*Δτ-0.42665*Δτ²-0.041833*Δτ³ (11)

q)转换矩阵：

A＝M₃(-Z_A)M₂(θ_A)M₃(-ζ_A) (12)

其中

其中:M₃(x)，M₂(x)表示函数，x为函数自变量；

r)位置向量

(3)、瞬时平赤道坐标系MOD到瞬时真赤道坐标系TOD的转换

g)平黄赤交角：

h)章动转换矩阵：

N＝M₁(-ε-Δε)M₃(-Δψ)M₁(ε) (17)

其中

i)位置向量

(4)、瞬时真赤道坐标系TOD到准固定地球坐标系ECPF的转换

m)世界时：

UT1＝UTC+UT1UTC+DUT1 (20)

其中：UT1UTC为当前转换时刻的UTC与对应的世界时UT1之间的差值，其值可通过EopDE405星历表查到，DUT1为UT1UTC的误差，暂定为0；

n)由转换时刻的UTC时间对应的世界时UT1，从2000年1月1日12时起算的儒略世纪数：

o)格林尼治平恒星时角：

p)格林尼治真恒星时角：

其中：

q)转换矩阵：B＝M₃(θ_g) (24)

r)位置向量：

(5)、准固定地球坐标系ECPF到WGS-84坐标系的转换

e)极移修正转换矩阵：

其中：X_p、Y_p为地球自转极移角；

f)位置向量

其中：

为步骤1中卫星在J2000坐标系下的位置向量。

第二步中，将来自姿态测量组合的经、纬、高

转换到WGS-84坐标系下的表示，其转换步骤如下：

e)卯酉圈曲率半径：

f)位置向量

第三步中，将站星方向矢量转换为n系下的表示，实施方式如下：

减去姿态测量组合的经、纬、高在WGS-84坐标系下的表示

即：

其中：

为站星方向在WGS-84坐标系下的表示，

分别为向量

对应的三个坐标表示。

g)转换矩阵：

h)

其中

为站星方向在n系下的三个坐标表示；

i)对

进行归一化处理。

第四步中，对用户站进行定姿态，将用户站由n系转换到b系，具体实施方式如下：

g)由姿态测量组合给出的旋转角测量值

构建由n系向b系的旋转变换矩阵

h)转换矩阵

i)坐标转换：

其中

为用户站在b系下的三个坐标表示。

第五步中，计算天线执行角：

天线执行角

在b系的含义如下：

为离轴角，其是归一化天线方向矢量

与oz_b轴之间的夹角，范围0°-57°；

为旋转角，其是归一化天线方向矢量

在x_boy_b平面内的投影

与ox_b轴之间的夹角，顺时针方向为正；

天线执行角与归一化天线方向矢量之间的换算关系如下：

其中，载体坐标系ox_by_bz_b b系中，

e)原点o位于天线中心位置；

f)ox_b轴指向载体正前方；

g)oy_b轴指向载体正左方；

h)oz_b轴指向载体正上方。

Claims

1.一种适用于卫星用户站数字跟踪接收机的程序跟踪方法，其特征在于，所述方法的步骤如下：

第一步，将星历轨道根数(a,e,i,Ω,ω,θ)转换到WGS-84坐标系下的表示；将星历轨道根数(a,e,i,Ω,ω,θ)转换到WGS-84坐标系下的表示，包括如下转换步骤：

(1)、星历轨道根数到J2000惯性坐标系的转换

a)

b)

c)第一向量

d)第二向量

e)卫星在J2000坐标系下的位置向量

其中，星历轨道根数(a,e,i,Ω,ω,θ)定义如下：a为半长轴；e为偏心率；i为轨道倾角；Ω为升交点赤经；ω为近地点幅角；θ为真近点角，μ为地球引力常量，其值为3.986004418×10¹⁴；

(2)、J2000惯性坐标系到瞬时平赤道坐标系MOD的转换

a)将协调世界时转换成儒略日：

其中fix()为取整函数；

b)太阳系质心动力学时：

其中LeapSecond为跳秒数，值为常数37；

c)由太阳系质心动力学时确定对应的儒略世纪数：

d)赤道面进动欧拉角：

ζ_A＝2306.2181*Δτ+0.30188*Δτ²+0.017998*Δτ³ (9)

Z_A＝2306.2181*Δτ+1.09468*Δτ²+0.018203*Δτ³ (10)

θ_A＝2004.3109*Δτ-0.42665*Δτ²-0.041833*Δτ³ (11)

e)转换矩阵：

A＝M₃(-Z_A)M₂(θ_A)M₃(-ζ_A) (12)

其中

其中:M₃(x)，M₂(x)表示函数，x为函数自变量；

f)位置向量

(3)、瞬时平赤道坐标系MOD到瞬时真赤道坐标系TOD的转换

a)平黄赤交角：

b)章动转换矩阵：

N＝M₁(-ε-Δε)M₃(-Δψ)M₁(ε) (17)

其中

c)位置向量

(4)、瞬时真赤道坐标系TOD到准固定地球坐标系ECPF的转换

a)世界时：

UT1＝UTC+UT1UTC+DUT1 (20)

其中：UT1UTC为当前转换时刻的UTC与对应的世界时UT1之间的差值，UT1UTC值可通过EopDE405星历表查到，DUT1为UT1UTC的误差，定为0；

b)由转换时刻的UTC时间对应的世界时UT1，从2000年1月1日12时起算的儒略世纪数：

c)格林尼治平恒星时角：

d)格林尼治真恒星时角：

其中：

为步骤4中c中的格林尼治平恒星时角，Δψ为黄经章动；Δε为交角章动；

其值可通过NutationDE405星历表查到，ε为步骤(3)中求得的平黄赤交角；

e)转换矩阵：B＝M₃(θ_g) (24)

f)位置向量：

(5)、准固定地球坐标系ECPF到WGS-84坐标系的转换

a)极移修正转换矩阵：

其中：X_p、Y_p为地球自转极移角；

b)位置向量

其中：

为卫星在WGS-84坐标系下的位置向量，W为极移修正转换矩阵，B为步骤(4)中的转换矩阵，N为步骤(3)中的章动转换矩阵，A为步骤(2)中的转换矩阵，

为步骤(1)中卫星在J2000坐标系下的位置向量；

第二步，将来自姿态测量组合的经、纬、高

转换到WGS-84坐标系下的表示；

其转换步骤如下：

a)卯酉圈曲率半径：

b)位置向量

其中a为地球半长轴，其数值为6378137.0米，b为地球半短轴，其数值为6356752.3142米；

第三步，将站星方向矢量转换为n系下的表示；

实施方式如下：

减去姿态测量组合的经、纬、高在WGS-84坐标系下的表示

即：

其中：

为站星方向在WGS-84坐标系下的表示，

分别为向量

对应的三个坐标表示；

a)转换矩阵：

b)

其中

为站星方向在n系下的三个坐标表示；

c)对

进行归一化处理；

第四步，对用户站进行定姿态，即将用户站坐标由n系转换到b系；

实施方式如下：

a)由姿态测量组合给出的旋转角测量值

构建由n系向b系的旋转变换矩阵

b)转换矩阵

c)坐标转换：

其中

为用户站在b系下的三个坐标表示；

第五步，计算天线执行角；

天线执行角

在b系的含义如下：

为离轴角，其是归一化天线方向矢量

与oz_b轴之间的夹角，范围0°-57°；

为旋转角，其是归一化天线方向矢量

在x_boy_b平面内的投影

与ox_b轴之间的夹角，顺时针方向为正；

天线执行角与归一化天线方向矢量之间的换算关系如下：

其中，载体坐标系ox_by_bz_bb系中，

a)原点o位于天线中心位置；

b)ox_b轴指向载体正前方；

c)oy_b轴指向载体正左方；

d)oz_b轴指向载体正上方。