CN111427002A - 地面测控天线指向卫星的方位角计算方法 - Google Patents

地面测控天线指向卫星的方位角计算方法 Download PDF

Info

Publication number
CN111427002A
CN111427002A CN202010197399.2A CN202010197399A CN111427002A CN 111427002 A CN111427002 A CN 111427002A CN 202010197399 A CN202010197399 A CN 202010197399A CN 111427002 A CN111427002 A CN 111427002A
Authority
CN
China
Prior art keywords
satellite
antenna
ground
station
time
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202010197399.2A
Other languages
English (en)
Other versions
CN111427002B (zh
Inventor
俞航
易灵
郑渊
刘登岭
辛阁霖
宋效正
吕旺
边志强
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shanghai Institute of Satellite Engineering
Original Assignee
Shanghai Institute of Satellite Engineering
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shanghai Institute of Satellite Engineering filed Critical Shanghai Institute of Satellite Engineering
Priority to CN202010197399.2A priority Critical patent/CN111427002B/zh
Publication of CN111427002A publication Critical patent/CN111427002A/zh
Application granted granted Critical
Publication of CN111427002B publication Critical patent/CN111427002B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S1/00Beacons or beacon systems transmitting signals having a characteristic or characteristics capable of being detected by non-directional receivers and defining directions, positions, or position lines fixed relatively to the beacon transmitters; Receivers co-operating therewith
    • G01S1/02Beacons or beacon systems transmitting signals having a characteristic or characteristics capable of being detected by non-directional receivers and defining directions, positions, or position lines fixed relatively to the beacon transmitters; Receivers co-operating therewith using radio waves
    • G01S1/08Systems for determining direction or position line
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

本发明公开了一种地面测控天线指向卫星的方位角计算方法,包括:利用以星历时刻t0、t0时刻的轨道平根数[a0,e0,i000,M0]和地面测站天线的地理位置信息,经过卫星轨道的相关计算以及多个相关坐标系的转换计算,最终转换为t1时刻在站心系下的天线指向方位角度:高低角
Figure DDA0002418106560000011
水平角ψ,以完成地面测站天线对卫星指向方位的预报过程。本发明不依赖于仿真软件或过多假设内容,考虑卫星实际运行情况计算地面测站天线对卫星的指向,并且在计算中考虑地球引力场摄动中的J2至J4项摄动,适合长时间、高精度的卫星位置预测,可以适用于卫星在轨的多种任务需求,有效解决了地面站接收天线对卫星的指向控制问题,而且达到了比较高的指向精度。

Description

地面测控天线指向卫星的方位角计算方法
技术领域
本发明涉及轨道动力学中的卫星轨道计算领域,具体地,涉及地面测控天线指向卫星的方位角计算方法。
背景技术
雷达在我国科技建设中起着非常重要的作用,随着外空间目标探测、控制的需求,实用雷达快速发展,目前已应用到制导和超视距探测等多个重要领域。随着空间侦查技术的发展,对雷达天线的跟踪和搜索能力提出了越来越高的要求,由于卫星信号微弱而且方向性强,为了捕获运动卫星上的通讯信号,必须实时调整天线姿态与卫星的位置偏差以满足通讯的需求,因为星地指向的偏差会导致链路传输信息的信噪比下降,若是超过了最大站位容限,甚至会出现信号丢失的现象。这就要求雷达天线必须要根据指令调整指向,实时跟踪运动目标。因此,雷达天线指向过程的动态精度已经成为天线***功能的重要指标之一,设计指向精度高的指向计算方法具有普遍的实用意义。
地面测站天线对卫星的指向方位预报,主要是依据卫星的轨道信息、时间信息以及地面测站天线的位置信息,计算出目标时刻的天线指向方位角,并通过天线指向控制,实现在目标时刻对卫星的精确指向。近年来对于人造卫星在指定时刻位置的预测,已经成为越来越重要的问题。高精度的轨道预报作为航天技术中一项重要技术,在卫星轨道设计和轨道优化中起着至关重要的作用,同时可以为天线指向、跟踪定位等卫星在轨任务提供可靠的轨道信息参考。
国内现有的星地指向算法研究,大多集中于地面测站位置固定的情况下,在卫星本体坐标系下,进行卫星对地面测站指向的优化设计,而对于地面测站天线对卫星的指向的研究较少。在地面测站可移动的情况下,需要在已知自身地理经纬度及高程,被追踪在轨飞行器轨道信息的情况下,自主完成对飞行器的指向控制,实时计算天线指向方位角度。目前,应用于地面站的卫星轨道递推方法和高精度定轨算法一般由高性能计算机实现,其中包括了高精度积分算法和高精度动力学模型。本发明针对这一实际情况,提出一种用于地面测站天线的精度较高的地面测站实时定位方法,实现地面测站天线对卫星的自主指向方位的预报。
专利“一种深空探测器天线指向的设计方法”(专利号:CN104369877A)介绍了一种深空探测器天线对地心指向的方法,此方法是用于实现深空探测器天线对地心的定向。该专利针对天线对地心的定向而并非是地表的给定位置,而且直接给出了探测器天线对地心的指向矢量,没有通过轨道参数计算卫星位置的算法。本发明与其不同之处在于,设计了针对地表给定位置的地面站对卫星指向的计算方法,完成了地面站的定位计算,而且设计了通过给定时刻的卫星轨道参数计算地面站-卫星指向矢量的计算过程。
专利“一种数传天线指向角度的仿真分析方法”(专利号:CN105184002A)介绍了一种计算星载数传天线对地面站指向的方法,此方法通过已有的卫星轨道仿真软件STK对卫星的实际位置进行仿真求解,并计算数传天线的二维指向角。该专利的不足在于卫星位置计算依赖于卫星轨道仿真软件STK,没有具体的计算过程,而且对坐标系转换的描述比较简略,没有给出转换矩阵的算法。本发明的优势在于提出了一种无需依靠STK软件,根据指定时刻的卫星的轨道参数计算卫星实际位置的方法,以及设计了一套相关坐标系转换矩阵的详细计算流程。
专利“一种绕月卫星双轴天线对地指向控制方法”(专利号:CN101204994A)介绍了一种计算绕月卫星对地心指向的方法,根据地面上注星历数据推算卫星位置,计算卫星对地球的可见区域,并计算出双轴天线的指向角度。该专利为对地心指向,没有对地表位置进行定向,且主要结合月球相关坐标系进行计算。本发明与其不同之处在于主要结合地球及地表位置相关坐标系进行计算,完成了地面测站天线的位置计算及地面测站天线对卫星的定向计算,且天线指向角度的定义和计算方法不同。
李丹,于洋等人在“基于轨道根数的低轨卫星轨道预测算法”(见《光学精密工程》,2016年,10期)论文中提出了一种利用椭圆曲线来预测卫星轨道的方法,但是在求解过程中需要计算系数的偏微分。
中国发明专利“一种适用于圆轨道卫星的星上自主轨道外推方法”(专利号:CN103995800A)中,介绍了一种适用于圆轨道卫星的轨道递推的方法。但是该方法只考虑了摄动。
基于以上考虑,本发明公开的一种地面测控天线指向卫星的方位角计算方法,其优势在于,精度高,适用于长期预报。
发明内容
针对现有技术中的缺陷,本发明的目的是提供一种地面测控天线指向卫星的方位角计算方法。
根据本发明提供的一种地面测控天线指向卫星的方位角计算方法,包括:
利用星历时刻t0、t0时刻的轨道平根数[a0,e0,i000,M0]和地面测站天线的地理位置信息,经过卫星轨道的相关计算以及多个相关坐标系的转换计算,最终转换为t1时刻在站心系下的天线指向方位角度:高低角
Figure BDA0002418106540000031
水平角ψ,以完成地面测站天线对卫星指向方位的预报过程。进一步地,星历时刻t0、t0时刻的轨道平根数[a0,e0,i000,M0]和地面测站天线的地理位置信息这些参数为已知量,通过查询卫星的星历表直接获得。
优选地,包括:
输入参数无奇异处理步骤:
当轨道偏心率很小,即近圆轨道,e≈0时,为了避免计算中出现奇点,转化为3个无奇点变量,即:
ξ0=e0cos(ω0)
η0=-e0sin(ω0)
λ0=ω0+M0
优选地,还包括:
输入参数归一化处理步骤:
对于递推时间dt进行归一化处理:
dt=t1-t0
Figure BDA0002418106540000032
其中,
t0为轨道递推起始时刻;
t1表示轨道递推结束时刻;
dtn表示递推时间经过归一化处理的结果,下标n表示归一化;
对于半长轴a0进行归一化处理:
Figure BDA0002418106540000033
所述的时间的归一化单位为
Figure BDA0002418106540000041
其中,
Ge为地心引力常数;
Re表示地球赤道半径。
优选地,还包括:
摄动项计算步骤:
考虑地球引力场摄动中的J2至J4项摄动,包含计算一阶长期项、一阶短周期项和二阶长期项,各个摄动项的表达式如下:
一阶长期项计算:
Figure BDA0002418106540000042
Figure BDA0002418106540000043
Figure BDA0002418106540000044
其中,
Ω1表示卫星轨道升交点赤经的一阶长周期变化量;
ω1表示卫星轨道近地点幅角的一阶长周期变化量;
λ1表示本文计算所引入的无奇点变量λ的一阶长周期变化量;
所述p为半通径,表达式为p=a×(1-e0 2),轨道平均角速度
Figure BDA0002418106540000045
J2=1.624×10-3
一阶短周期项计算:
Figure BDA0002418106540000046
Figure BDA0002418106540000047
Figure BDA0002418106540000048
Figure BDA0002418106540000049
Figure BDA00024181065400000410
Figure BDA0002418106540000051
其中,
Figure BDA0002418106540000052
表示卫星轨道半长轴的一阶短周期变化量;
Figure BDA0002418106540000053
表示卫星轨道倾角的一阶短周期变化量;
Figure BDA0002418106540000054
表示卫星轨道升交点赤经的一阶短周期变化量;
Figure BDA0002418106540000055
分别表示本文计算所引入的三个无奇点变量ξ、η、λ的一阶短周期变化量;
所述u由一阶长期项计算得到,计算过程如下:
ξz1=ξ0 cos(ω1 dtn)+η0 sin(ω1 dtn)
ηz1=η0 cos(ω1 dtn)-ξ0 sin(ω1 dtn)
λz1=λ0+(n+λ1)dtn
Figure BDA0002418106540000056
ωz1=arc cos(ξz1/ez1)
Mz1=mod(λz1z1,2π)
Figure BDA0002418106540000057
u=fz1z1
其中,
ξz1、ηz1、λz1分别表示本文计算所引入的三个无奇点变量ξ、η、λ的根据一阶长期变化量的积分结果;
ez1表示卫星轨道偏心率以一阶长周期变化量为积分量的积分结果;
ωz1表示卫星轨道近地点幅角以一阶长周期变化量为积分量的积分结果;
Mz1表示卫星轨道平近点角以一阶长周期变化量为积分量的积分结果;
fz1表示卫星轨道真近点角以一阶长周期变化量为积分量的积分结果;
二阶长期项计算:
Figure BDA0002418106540000061
Figure BDA0002418106540000062
Figure BDA0002418106540000063
其中,
Ω2表示卫星轨道升交点赤经的二阶长周期变化量;
ξ2、λ2分别表示本文计算所引入的无奇点变量ξ、λ的二阶长周期变化量;
J3=2.5356×10-6,J4=7.1022×10-6
优选地,还包括:
递推主公式计算步骤:
对于引入的无奇点变量,结合长期项和短周期项摄动项进行解析解构造,忽略长周期项摄动,轨道递推主公式如下:
Figure BDA0002418106540000064
Figure BDA0002418106540000065
Figure BDA0002418106540000066
Figure BDA0002418106540000067
Figure BDA0002418106540000068
Figure BDA0002418106540000069
其中,
αt表示表示t1时刻,卫星轨道半长轴的归一化结果;
is表示t1时刻,卫星轨道倾角;
Ωs表示t1时刻,卫星轨道升交点赤经;
ξs、ηs、λs分别表示t1时刻,计算所引入的三个无奇点变量ξ、η、λ的值。
优选地,还包括:
无奇点变量还原步骤:
计算结束后将3个无奇点变量还原
Figure BDA0002418106540000071
ωs=arc cos(ξs/es)
Ms=mod(λss,2π)
优选地,还包括:
归一化变量还原步骤:
将卫星轨道长轴at还原为常规单位:
as=at×Re
所述的as单位:m。
优选地,还包括:
惯性系卫星位置计算步骤:
输入t1时刻卫星轨道瞬根数[as,es,isss,Ms],输出卫星的位置矢量在J2000.0坐标系下的分量RwECI,计算方法为:
RwECI=Q*rp
其中,旋转矩阵Q按照3-1-3旋转顺序进行描述:
Figure BDA0002418106540000072
矢量rp
Figure BDA0002418106540000073
其中,M1为真近点角:
Figure BDA0002418106540000074
优选地,还包括:
地固系卫星位置计算步骤:
输入目标时刻t1和所述计算出的卫星在惯性系下的位置RwECI,输出t1时刻卫星在地固系下的位置RwECF,具体方法如下:
根据给定的目标时刻t1,计算历元J2000.0(2000年1月1日12时)至给定目标时刻的秒计数值tc,输入t1时刻的年、月、日、时、分、秒,计算儒略日JD:
Figure BDA0002418106540000081
其中,floor()为向下取整运算;
根据儒略日JD计算历元J2000.0至给定目标时刻的秒计数值tc
tc=(JD-2455197.5)×86400+315547200
根据计算得到的历元J2000.0至给定目标时刻的秒计数值tc,计算地球自转矩阵ER、章动矩阵NR以及岁差矩阵PR,计算惯性系到地固系的转换矩阵MECI2ECF
MECI2ECF=ER*NR*PR
并根据所述计算得到的t1时刻卫星在惯性系下的位置RwECI,计算t1时刻卫星在地固系下的位置RwECF
RwECF=MECI2ECF*RwECI
优选地,还包括:
地固系测站天线位置计算步骤:
根据给定的地面测站天线的经纬度和高程,计算地面测站天线在地固系下的位置RtECF,具体方法如下:
输入地面测站天线的经度lon、纬度lat、高程h;
计算坐标分量G1、G2:
Figure BDA0002418106540000082
Figure BDA0002418106540000083
其中,f为地球椭球体几何扁率,f=1/298.257;
计算地面测站天线在地固系下的位置RtECF
Figure BDA0002418106540000091
优选地,还包括:
站心系卫星位置计算步骤:
输入计算出的t1时刻卫星在地固系下的位置RwECF和所述计算出的地面测站天线在地固系下的位置RtECF,输出刻卫星在站心系下的位置RwCT,具体方法如下:
在站心系下,计算地固系到站心系的转换矩阵MECF2CT,描述为一次绕地固系的Z轴的旋转和一次绕地固系的X轴的旋转:
MECF2CT=Rx(90°-lat)Rz(90°+lon)
其中,lon为地面测站天线的地理经度;lat为地面测站天线的地理纬度;
Figure BDA0002418106540000092
Figure BDA0002418106540000093
并根据所述t1时刻卫星在地固系下的位置RwECF、地面测站天线在地固系下的位置RtECF,将坐标原点由地心平移到地面测站天线处,并计算t1时刻卫星在站心系下的位置RwCT
RwCT=MECF2CT*(RwECF-RtECF)
优选地,还包括:
测站天线指向方位角计算步骤:
输入卫星在站心系下的位置,输出t1时刻地面测站天线指向方位角度:高低角
Figure BDA0002418106540000094
水平角ψ,具体方法如下:
高低角、水平角在站心系下定义,在站心系OCTXCTYCTZCT中,OCT表示坐标系原点,XCT表示坐标系的X轴,YCT表示坐标系的Y轴,OCTXCTYCTZCT表示该坐标系的X轴与Y轴所在的平面,高低角
Figure BDA0002418106540000101
为指向矢量RwCT与OCTXCTYCT平面的夹角,定义RwCT矢量与OCTZCT夹角小于90°为正;水平角ψ为指向矢量RwCT在OCTXCTYCT平面的投影与OCTXCT轴的夹角,定义绕OCTZCT轴从OCTXCT轴顺时针转向指向矢量RwCT在OCTXCTYCT面的投影为正,根据此定义求出天线指向方位角。假设地面测站的位置位于站心系的原点,则测站天线指向矢量在站心系下的投影为RwCT,记RwCT为:
Figure BDA0002418106540000102
其中,
xCT、yCT、zCT分别表示测站天线指向矢量RwCT,在站心系OCTXCTYCTZCT中,对应的X轴、Y轴、Z轴的坐标分量;
对天线指向方位角的计算高低角
Figure BDA0002418106540000103
水平角ψ为:
Figure BDA0002418106540000104
Figure BDA0002418106540000105
并根据xCT、yCT、zCT的正负,将高低角
Figure BDA0002418106540000106
水平角ψ划分到对应的角度范围内,完成天线指向方位预报过程:
若zCT≥0,则将高低角
Figure BDA0002418106540000107
划分到范围
Figure BDA0002418106540000108
若xCT≥0,yCT≥0,则将水平角ψ划分到范围
Figure BDA0002418106540000109
若xCT<0,yCT≥0,则将水平角ψ划分到范围
Figure BDA00024181065400001010
若xCT<0,yCT<0,则将水平角ψ划分到范围
Figure BDA00024181065400001011
若xCT≥0,yCT<0,则将水平角ψ划分到范围
Figure BDA00024181065400001012
与现有技术相比,本发明具有如下的有益效果:
本发明不依赖于仿真软件或过多假设内容,考虑卫星实际运行情况计算地面测站天线对卫星的指向,并且在计算中考虑地球引力场摄动中的J2至J4项摄动,适合长时间、高精度的卫星位置预测,可以适用于卫星在轨的多种任务需求,有效解决了地面站接收天线对卫星的指向控制问题,而且达到了比较高的指向精度。
附图说明
通过阅读参照以下附图对非限制性实施例所作的详细描述,本发明的其它特征、目的和优点将会变得更明显:
图1为一种地面测控天线指向卫星的方位角计算方法流程示意图。
图2为卫星的星上轨道位置自主预报示意图。
图3为地面测站天线对卫星指向方位示意图。
图4为站心系OCTXCTYCTZCT示意图。
图5为站心系下的高低角
Figure BDA0002418106540000112
与水平角ψ示意图。
具体实施方式
下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明,但不以任何形式限制本发明。应当指出的是,对本领域的普通技术人员来说,在不脱离本发明构思的前提下,还可以做出若干变化和改进。这些都属于本发明的保护范围。
根据本发明提供的一种地面测控天线指向卫星的方位角计算方法,包括:
利用星历时刻t0、t0时刻的轨道平根数[a0,e0,i000,M0]和地面测站天线的地理位置信息,经过卫星轨道的相关计算以及多个相关坐标系的转换计算,最终转换为t1时刻在站心系下的天线指向方位角度:高低角
Figure BDA0002418106540000111
水平角ψ,以完成地面测站天线对卫星指向方位的预报过程。进一步地,星历时刻t0、t0时刻的轨道平根数[a0,e0,i000,M0]和地面测站天线的地理位置信息这些参数为已知量,通过查询卫星的星历表直接获得。
具体地,包括:
输入参数无奇异处理步骤:
当轨道偏心率很小,即近圆轨道,e≈0时,为了避免计算中出现奇点,转化为3个无奇点变量,即:
ξ0=e0 cos(ω0)
η0=-e0 sin(ω0)
λ0=ω0+M0
具体地,还包括:
输入参数归一化处理步骤:
对于递推时间dt进行归一化处理:
dt=t1-t0
Figure BDA0002418106540000121
其中,
t0为轨道递推起始时刻;
t1表示轨道递推结束时刻;
dtn表示递推时间经过归一化处理的结果,下标n表示归一化;
对于半长轴a0进行归一化处理:
Figure BDA0002418106540000122
所述的时间的归一化单位为
Figure BDA0002418106540000123
其中,
Ge为地心引力常数;
Re表示地球赤道半径。
具体地,还包括:
摄动项计算步骤:
考虑地球引力场摄动中的J2至J4项摄动,包含计算一阶长期项、一阶短周期项和二阶长期项,各个摄动项的表达式如下:
一阶长期项计算:
Figure BDA0002418106540000131
Figure BDA0002418106540000132
Figure BDA0002418106540000133
其中,
Ω1表示卫星轨道升交点赤经的一阶长周期变化量;
ω1表示卫星轨道近地点幅角的一阶长周期变化量;
λ1表示本文计算所引入的无奇点变量λ的一阶长周期变化量;
所述p为半通径,表达式为p=a×(1-e0 2),轨道平均角速度
Figure BDA0002418106540000134
J2=1.624×10-3
一阶短周期项计算:
Figure BDA0002418106540000135
Figure BDA0002418106540000136
Figure BDA0002418106540000137
Figure BDA0002418106540000138
Figure BDA0002418106540000139
Figure BDA00024181065400001310
其中,
Figure BDA00024181065400001311
表示卫星轨道半长轴的一阶短周期变化量;
Figure BDA00024181065400001312
表示卫星轨道倾角的一阶短周期变化量;
Figure BDA00024181065400001313
表示卫星轨道升交点赤经的一阶短周期变化量;
Figure BDA00024181065400001314
分别表示本文计算所引入的三个无奇点变量ξ、η、λ的一阶短周期变化量;
所述u由一阶长期项计算得到,计算过程如下:
ξz1=ξ0 cos(ω1 dtn)+η0 sin(ω1 dtn)
ηz1=η0 cos(ω1 dtn)-ξ0 sin(ω1 dtn)
λz1=λ0+(n+λ1)dtn
Figure BDA0002418106540000141
ωz1=arc cos(ξz1/ez1)
Mz1=mod(λz1z1,2π)
Figure BDA0002418106540000142
u=fz1z1
其中,
ξz1、ηz1、λz1分别表示本文计算所引入的三个无奇点变量ξ、η、λ的根据一阶长期变化量的积分结果;
ez1表示卫星轨道偏心率以一阶长周期变化量为积分量的积分结果;
ωz1表示卫星轨道近地点幅角以一阶长周期变化量为积分量的积分结果;
Mz1表示卫星轨道平近点角以一阶长周期变化量为积分量的积分结果;
fz1表示卫星轨道真近点角以一阶长周期变化量为积分量的积分结果;
二阶长期项计算:
Figure BDA0002418106540000143
Figure BDA0002418106540000144
Figure BDA0002418106540000145
Figure BDA0002418106540000146
其中,
Ω2表示卫星轨道升交点赤经的二阶长周期变化量;
ξ2、λ2分别表示本文计算所引入的无奇点变量ξ、λ的二阶长周期变化量;
J3=2.5356×10-6,J4=7.1022×10-6
具体地,还包括:
递推主公式计算步骤:
对于引入的无奇点变量,结合长期项和短周期项摄动项进行解析解构造,忽略长周期项摄动,轨道递推主公式如下:
Figure BDA0002418106540000151
Figure BDA0002418106540000152
Figure BDA0002418106540000153
Figure BDA0002418106540000154
Figure BDA0002418106540000155
Figure BDA0002418106540000156
其中,
αt表示表示t1时刻,卫星轨道半长轴的归一化结果;
is表示t1时刻,卫星轨道倾角;
Ωs表示t1时刻,卫星轨道升交点赤经;
ξs、ηs、λs分别表示t1时刻,计算所引入的三个无奇点变量ξ、η、λ的值。
具体地,还包括:
无奇点变量还原步骤:
计算结束后将3个无奇点变量还原
Figure BDA0002418106540000157
ωs=arc cos(ξs/es)
Ms=mod(λss,2π)
具体地,还包括:
归一化变量还原步骤:
将卫星轨道长轴at还原为常规单位:
as=at×Re
所述的as单位:m。
具体地,还包括:
惯性系卫星位置计算步骤:
输入t1时刻卫星轨道瞬根数[as,es,isss,Ms],输出卫星的位置矢量在J2000.0坐标系下的分量RwECI,计算方法为:
RwECI=Q*rp
其中,旋转矩阵Q按照3-1-3旋转顺序进行描述:
Figure BDA0002418106540000161
矢量rp
Figure BDA0002418106540000162
其中,M1为真近点角:
Figure BDA0002418106540000163
具体地,还包括:
地固系卫星位置计算步骤:
输入目标时刻t1和所述计算出的卫星在惯性系下的位置RwECI,输出t1时刻卫星在地固系下的位置RwECF,具体方法如下:
根据给定的目标时刻t1,计算历元J2000.0(2000年1月1日12时)至给定目标时刻的秒计数值tc,输入t1时刻的年、月、日、时、分、秒,计算儒略日JD:
Figure BDA0002418106540000164
其中,floor()为向下取整运算;
根据儒略日JD计算历元J2000.0至给定目标时刻的秒计数值tc
tc=(JD-2455197.5)×86400+315547200
根据计算得到的历元J2000.0至给定目标时刻的秒计数值tc,计算地球自转矩阵ER、章动矩阵NR以及岁差矩阵PR,计算惯性系到地固系的转换矩阵MECI2ECF
MECI2ECF=ER*NR*PR
并根据所述计算得到的t1时刻卫星在惯性系下的位置RwECI,计算t1时刻卫星在地固系下的位置RwECF
RwECF=MECI2ECF*RwECI
具体地,还包括:
地固系测站天线位置计算步骤:
根据给定的地面测站天线的经纬度和高程,计算地面测站天线在地固系下的位置RtECF,具体方法如下:
输入地面测站天线的经度lon、纬度lat、高程h;
计算坐标分量G1、G2:
Figure BDA0002418106540000171
Figure BDA0002418106540000172
其中,f为地球椭球体几何扁率,f=1/298.257;
计算地面测站天线在地固系下的位置RtECF
Figure BDA0002418106540000173
具体地,还包括:
站心系卫星位置计算步骤:
输入计算出的t1时刻卫星在地固系下的位置RwECF和所述计算出的地面测站天线在地固系下的位置RtECF,输出刻卫星在站心系下的位置RwCT,具体方法如下:
在站心系下,计算地固系到站心系的转换矩阵MECF2CT,描述为一次绕地固系的Z轴的旋转和一次绕地固系的X轴的旋转:
MECF2CT=Rx(90°-lat)Rz(90°+lon)
其中,lon为地面测站天线的地理经度;lat为地面测站天线的地理纬度;
Figure BDA0002418106540000181
Figure BDA0002418106540000182
并根据所述t1时刻卫星在地固系下的位置RwECF、地面测站天线在地固系下的位置RtECF,将坐标原点由地心平移到地面测站天线处,并计算t1时刻卫星在站心系下的位置RwCT
RwCT=MECF2CT*(RwECF-RtECF)
具体地,还包括:
测站天线指向方位角计算步骤:
输入卫星在站心系下的位置,输出t1时刻地面测站天线指向方位角度:高低角
Figure BDA0002418106540000183
水平角ψ,具体方法如下:
高低角、水平角在站心系下定义,在站心系OCTXCTYCTZCT中,OCT表示坐标系原点,XCT表示坐标系的X轴,YCT表示坐标系的Y轴,OCTXCTYCT表示该坐标系的X轴与Y轴所在的平面,高低角
Figure BDA0002418106540000184
为指向矢量RwCT与OCTXCTYCT平面的夹角,定义RwCT矢量与OCTZCT夹角小于90°为正;水平角ψ为指向矢量RwCT在OCTXCTYCT平面的投影与OCTXCT轴的夹角,定义绕OCTZCT轴从OCTXCT轴顺时针转向指向矢量RwCT在OCTXCTYCT面的投影为正,根据此定义求出天线指向方位角。假设地面测站的位置位于站心系的原点,则测站天线指向矢量在站心系下的投影为RwCT,记RwCT为:
Figure BDA0002418106540000185
其中,
xCT、yCT、zCT分别表示测站天线指向矢量RwCT,在站心系OCTXCTYCTZCT中,对应的X轴、Y轴、Z轴的坐标分量;
对天线指向方位角的计算高低角
Figure BDA0002418106540000191
水平角ψ为:
Figure BDA0002418106540000192
Figure BDA0002418106540000193
并根据xCT、yCT、zCT的正负,将高低角
Figure BDA0002418106540000194
水平角ψ划分到对应的角度范围内,完成天线指向方位预报过程:
若zCT≥0,则将高低角
Figure BDA0002418106540000195
划分到范围
Figure BDA0002418106540000196
若xCT≥0,yCT≥0,则将水平角ψ划分到范围
Figure BDA0002418106540000197
若xCT<0,yCT≥0,则将水平角ψ划分到范围
Figure BDA0002418106540000198
若xCT<0,yCT<0,则将水平角ψ划分到范围
Figure BDA0002418106540000199
若xCT≥0,yCT<0,则将水平角ψ划分到范围
Figure BDA00024181065400001910
下面通过优选例,对本发明进行更为具体地说明。
优选例1:
本发明的技术解决问题是:克服现有技术的不足,提供一种精度较高的地面测控天线指向卫星的方位角计算方法,通过给定目标时刻的卫星星历数据和地面测站天线的地理位置信息,经过卫星轨道的相关计算以及多个相关坐标系的转换计算,最终转换为在站心系下的天线指向方位角以完成地面测站天线对卫星方位的预报过程。
本发明结合一种工程实际情况:在地面测站可移动的情况下,需要在已知自身地理经纬度及高程,被追踪在轨飞行器轨道信息的情况下,自主完成对飞行器的指向跟踪,实时计算天线指向方位角。提出一种用于地面测站天线的精度较高的地面测站实时定位方法,实现地面测站天线对卫星的自主指向方位的预报。
本发明提出的地面测控天线指向卫星的方位角计算方法,通过卫星星历信息实时计算卫星的位置,且对坐标系转换关系的影响因素考虑全面,计算精度较高,并提出了一种适用于地面测站天线的指向方位角定义,有效解决了地面测站天线实时对卫星指向方位的预报需求。
为了实现该方法,采用如下的技术方案:
一种地面测控天线指向卫星的方位角计算方法,其过程如下:
1.输入参数无奇异处理模块
轨道偏心率很小(近圆轨道,e≈0)时,为了避免计算中出现奇点,转化为3个无奇点变量,即
ξ0=e0 cos(ω0)
η0=-e0 sin(ω0)
λ0=ω0+M0
2.输入参数归一化处理模块
对于递推时间dt进行归一化处理,
dt=t1-t0
Figure BDA0002418106540000201
对于半长轴a0进行归一化处理,
Figure BDA0002418106540000202
所述的时间的归一化单位为,
Figure BDA0002418106540000203
其中,Ge为地心引力常数;
所述的长度的归一化单位为Re=6378140m
3.摄动项计算模块
考虑地球引力场摄动中的J2至J4项摄动,包含计算一阶长期项、一阶短周期项和二阶长期项,各个摄动项的表达式:
一阶长期项计算:
Figure BDA0002418106540000211
Figure BDA0002418106540000212
Figure BDA0002418106540000213
所述的p,为半通径,表达式为p=a×(1-e0 2),轨道平均角速度
Figure BDA0002418106540000214
J2=1.624×10-3
一阶短周期项计算:
Figure BDA00024181065400002113
Figure BDA0002418106540000215
Figure BDA0002418106540000216
Figure BDA0002418106540000217
Figure BDA0002418106540000218
Figure BDA0002418106540000219
所述的,u由一阶长期项计算得到。
Figure BDA00024181065400002110
Figure BDA00024181065400002111
ξz1=ξ0 cos(ωc1dtn)+η0 sin(ωc1dtn)
ηz1=η0 cos(ωc1dtn)-ξ0 sin(ωc1dtn)
λz1=λ0+(n+λc1)dtn
Figure BDA00024181065400002112
ωz1=arc cos(ξz1/ez1)
Mz1=mod(λz1z1,2π)
Figure BDA0002418106540000221
u=fz1z1
二阶长期项计算:
Figure BDA0002418106540000222
Figure BDA0002418106540000223
Figure BDA0002418106540000224
所述的,J3=2.5356×10-6,J4=7.1022×10-6
4.递推主公式计算模块
对于引入的无奇点变量,结合长期项和短周期项摄动项进行解析解构造,忽略长周期项摄动。轨道递推主公式如下:
Figure BDA0002418106540000225
Figure BDA0002418106540000226
Figure BDA0002418106540000227
Figure BDA0002418106540000228
Figure BDA0002418106540000229
Figure BDA00024181065400002210
5.无奇点变量还原模块
计算结束后将3个无奇点变量还原
Figure BDA00024181065400002211
ωs=arc cos(ξs/es)
Ms=mod(λss,2π)
6.归一化变量还原模块
将卫星轨道长轴at还原为常规单位:
as=at×Re
所述的as单位:m。
7.惯性系卫星位置模块
通过t1时刻卫星轨道瞬根数[as,es,isss,Ms]计算卫星的位置矢量在J2000.0坐标系下的分量RwECI,计算方法为:
RwECI=Q*rp
其中,旋转矩阵Q按照3-1-3旋转顺序进行描述:
Figure BDA0002418106540000231
矢量rp
Figure BDA0002418106540000232
其中,M1为真近点角:
Figure BDA0002418106540000233
8.地固系卫星位置模块
根据给定的目标时刻t1和所述步骤(7)计算出的卫星在惯性系下的位置RwECI,计算t1时刻卫星在地固系下的位置RwECF的方法如下:
根据给定的目标时刻t1,计算历元J2000.0(2000年1月1日12时)至给定目标时刻的秒计数值tc,输入t1时刻(UTC时间)的年(year)、月(month)、日(day)、时(hour)、分(min)、秒(sec),计算儒略日JD:
Figure BDA0002418106540000234
其中,floor()为向下取整运算。
根据儒略日JD计算历元J2000.0(2000年1月1日12时)至给定目标时刻的秒计数值tc
tc=(JD-2455197.5)×86400+315547200
根据计算得到的历元J2000.0(2000年1月1日12时)至给定目标时刻的秒计数值tc,计算地球自转矩阵ER,章动矩阵NR,岁差矩阵PR,由于极移对转换矩阵的计算影响很小,本发明中不考虑此项。计算惯性系到地固系的转换矩阵MECI2ECF
MECI2ECF=ER*NR*PR
并根据所述步骤(7)计算得到的t1时刻卫星在惯性系下的位置RwECI,计算t1时刻卫星在地固系下的位置RwECF
RwECF=MECI2ECF*RwECI
9.地固系测站天线位置模块
其特征在于,根据给定的地面测站天线的经纬度和高程,计算地面测站天线在地固系下的位置RtECF的方法如下:
输入地面测站天线的经度lon、纬度lat、高程h。
计算坐标分量G1、G2:
Figure BDA0002418106540000241
Figure BDA0002418106540000242
其中,f为地球椭球体几何扁率,f=1/298.257。
计算地面测站天线在地固系下的位置RtECF
Figure BDA0002418106540000243
10.站心系卫星位置模块
根据所述步骤(8)计算出的t1时刻卫星在地固系下的位置RwECF和所述步骤(9)计算出的地面测站天线在地固系下的位置RtECF,计算t1时刻卫星在站心系下的位置RwCT的方法如下:
在站心系下,计算地固系到站心系的转换矩阵MECF2CT,描述为一次绕地固系的Z轴的旋转和一次绕地固系的X轴的旋转:
MECF2CT=Rx(90°-lat)Rz(90°+lon)
其中,lon为地面测站天线的地理经度;lat为地面测站天线的地理纬度。
Figure BDA0002418106540000251
Figure BDA0002418106540000252
并根据所述步骤(8)计算得到的t1时刻卫星在地固系下的位置RwECF,所述步骤(9)计算得到的地面测站天线在地固系下的位置RtECF,将坐标原点由地心平移到地面测站天线处,并计算t1时刻卫星在站心系下的位置RwCT
RwCT=MECF2CT*(RwECF-RtECF)
11.测站天线指向方位角模块
根据所述步骤(10)计算出的卫星在站心系下的位置,计算t1时刻地面测站天线指向方位角度:高低角
Figure BDA0002418106540000253
水平角ψ的方法如下:
高低角、水平角。高低角、水平角在站心系下定义,高低角
Figure BDA0002418106540000254
为指向矢量RwCT与OCTXCTYCT平面的夹角,定义RwCT矢量与OCTZCT夹角小于90°为正;水平角ψ为指向矢量RwCT在OCTXCTYCT平面的投影与OCTXCT轴的夹角,定义绕OCTZCT轴从OCTXCT轴顺时针转向指向矢量RwCT在OCTXCTYCT面的投影为正,如图5所示,根据此定义求出天线指向方位角。假设地面测站的位置位于站心系的原点,则测站天线指向矢量在站心系下的投影为RwCT,记RwCT为:
Figure BDA0002418106540000255
对天线指向方位角的计算高低角
Figure BDA0002418106540000256
水平角ψ为:
Figure BDA0002418106540000261
Figure BDA0002418106540000262
并根据xCT、yCT、zCT的正负,将高低角
Figure BDA0002418106540000263
水平角ψ划分到对应的角度范围内,完成天线指向方位的预报过程:
Figure BDA0002418106540000264
优选例2:
本发明需要用到的坐标系:所述的惯性系为J2000.0惯性坐标系,所述的地固系为WGS-84坐标系。下面给出站心系的定义。
站心系OCTXCTYCTZCT
站心系的定义为,原点OCT为地面天线原点,基本平面OCTXCTYCT面为当地水平面,OCTXCT沿当地子午圈指向正北,OCTZCT垂直基本平面指向天顶,OCTYCT按右手法则确定,如图4所示。
下面详述本发明的计算过程:
利用MATLAB对此算法进行仿真[a0,e0,i000,M0]验证,地球相关参数以及站心系按上文所述设置,某型号卫星在UTC时间2019年1月5日4时的星历数据如下:
Figure BDA0002418106540000265
Figure BDA0002418106540000271
从2019年1月5日4时递推2019年1月6日4时卫星的轨道六根数,并计算2019年1月6日4时的天线指向方位角(高低角、水平角)。
1.输入参数无奇异处理模块
轨道偏心率很小(近圆轨道,e≈0)时,为了避免计算中出现奇点,转化为3个无奇点变量,即
ξ0=e0 cos(ω0)
η0=-e0 sin(ω0)
λ0=ω0+M0
2.输入参数归一化处理模块
对于递推时间dt进行归一化处理,
dt=t1-t0
Figure BDA0002418106540000272
对于半长轴a0进行归一化处理,
Figure BDA0002418106540000273
所述的时间的归一化单位为,
Figure BDA0002418106540000274
其中,Ge为地心引力常数;
所述的长度的归一化单位为Re=6378140m
3.摄动项计算模块
考虑地球引力场摄动中的J2至J4项摄动,包含计算一阶长期项、一阶短周期项和二阶长期项,各个摄动项的表达式:
一阶长期项计算:
Figure BDA0002418106540000281
Figure BDA0002418106540000282
Figure BDA0002418106540000283
所述的p,为半通径,表达式为p=a×(1-e0 2),轨道平均角速度
Figure BDA0002418106540000284
J2=1.624×10-3
一阶短周期项计算:
Figure BDA0002418106540000285
Figure BDA0002418106540000286
Figure BDA0002418106540000287
Figure BDA0002418106540000288
Figure BDA0002418106540000289
Figure BDA00024181065400002810
所述的,u由一阶长期项计算得到。
Figure BDA00024181065400002811
Figure BDA00024181065400002812
ξz1=ξ0 cos(ωc1dtn)+η0 sin(ωc1dtn)
ηz1=η0 cos(ωc1dtn)-ξ0 sin(ωc1dtn)
λz1=λ0+(n+λc1)dtn
Figure BDA00024181065400002813
ωz1=arc cos(ξz1/ez1)
Mz1=mod(λz1z1,2π)
Figure BDA0002418106540000291
u=fz1z1
二阶长期项计算:
Figure BDA0002418106540000292
Figure BDA0002418106540000293
Figure BDA0002418106540000294
所述的,J3=2.5356×10-6,J4=7.1022×10-6
4.递推主公式计算模块
对于引入的无奇点变量,结合长期项和短周期项摄动项进行解析解构造,忽略长周期项摄动。轨道递推主公式如下:
Figure BDA0002418106540000295
Figure BDA0002418106540000296
Figure BDA0002418106540000297
Figure BDA0002418106540000298
Figure BDA0002418106540000299
Figure BDA00024181065400002910
5.无奇点变量还原模块
计算结束后将3个无奇点变量还原
Figure BDA00024181065400002911
ωs=arc cos(ξs/es)
Ms=mod(λss,2π)
6.归一化变量还原模块
将卫星轨道长轴at还原为常规单位:
as=at×Re
所述的as单位:m。
上述步骤计算得到的结果为:
Figure BDA0002418106540000301
7.惯性系卫星位置模块
通过t1时刻卫星轨道瞬根数[as,es,isss,Ms]计算卫星的位置矢量在J2000.0坐标系下的分量RwECI,计算方法为:
RwECI=Q*rp
其中,旋转矩阵Q按照3-1-3旋转顺序进行描述:
Figure BDA0002418106540000302
矢量rp
Figure BDA0002418106540000303
其中,M1为真近点角:
Figure BDA0002418106540000304
计算结果为:
Figure BDA0002418106540000305
Figure BDA0002418106540000311
8.地固系卫星位置模块
根据给定的目标时刻t1和所述步骤(7)计算出的卫星在惯性系下的位置RwECI,计算t1时刻卫星在地固系下的位置RwECF的方法如下:
根据给定的目标时刻t1,计算历元J2000.0(2000年1月1日12时)至给定目标时刻的秒计数值tc,输入t1时刻(UTC时间)的年(year)、月(month)、日(day)、时(hour)、分(min)、秒(sec),计算儒略日JD:
Figure BDA0002418106540000312
其中,floor()为向下取整运算。
根据儒略日JD计算历元J2000.0(2000年1月1日12时)至给定目标时刻的秒计数值tc
tc=(JD-2455197.5)×86400+315547200
根据计算得到的历元J2000.0(2000年1月1日12时)至给定目标时刻的秒计数值tc,计算地球自转矩阵ER,章动矩阵NR,岁差矩阵PR,由于极移对转换矩阵的计算影响很小,本发明中不考虑此项。计算惯性系到地固系的转换矩阵MECI2ECF
MECI2ECF=ER*NR*PR
并根据所述步骤(7)计算得到的t1时刻卫星在惯性系下的位置RwECI,计算t1时刻卫星在地固系下的位置RwECF
RwECF=MECI2ECF*RwECI
计算结果为:
Figure BDA0002418106540000313
Figure BDA0002418106540000321
9.地固系测站天线位置模块
其特征在于,根据给定的地面测站天线的经纬度和高程,计算地面测站天线在地固系下的位置RtECF的方法如下:
输入地面测站天线的经度lon、纬度lat、高程h。
计算坐标分量G1、G2:
Figure BDA0002418106540000322
Figure BDA0002418106540000323
其中,f为地球椭球体几何扁率,f=1/298.257。
计算地面测站天线在地固系下的位置RtECF
Figure BDA0002418106540000324
计算结果为:
Figure BDA0002418106540000325
10.站心系卫星位置模块
根据所述步骤(8)计算出的t1时刻卫星在地固系下的位置RwECF和所述步骤(9)计算出的地面测站天线在地固系下的位置RtECF,计算t1时刻卫星在站心系下的位置RwCT的方法如下:
在站心系下,计算地固系到站心系的转换矩阵MECF2CT,描述为一次绕地固系的Z轴的旋转和一次绕地固系的X轴的旋转:
MECF2CT=Rx(90°-lat)Rz(90°+lon)
其中,lon为地面测站天线的地理经度;lat为地面测站天线的地理纬度。
Figure BDA0002418106540000331
Figure BDA0002418106540000332
并根据所述步骤(8)计算得到的t1时刻卫星在地固系下的位置RwECF,所述步骤(9)计算得到的地面测站天线在地固系下的位置RtECF,将坐标原点由地心平移到地面测站天线处,并计算t1时刻卫星在站心系下的位置RwCT
RwCT=MECF2CT*(RwECF-RtECF)
计算结果为:
Figure BDA0002418106540000333
11.测站天线指向方位角模块
根据所述步骤(10)计算出的卫星在站心系下的位置,计算t1时刻地面测站天线指向方位角度:高低角
Figure BDA0002418106540000334
水平角ψ的方法如下:
高低角、水平角。高低角、水平角在站心系下定义,高低角
Figure BDA0002418106540000335
为指向矢量RwCT与OCTXCTYCT平面的夹角,定义RwCT矢量与OCTZCT夹角小于90°为正;水平角ψ为指向矢量RwCT在OCTXCTYCT平面的投影与OCTXCT轴的夹角,定义绕OCTZCT轴从OCTXCT轴顺时针转向指向矢量RwCT在OCTXCTYCT面的投影为正,如图5所示,根据此定义求出天线指向方位角。假设地面测站的位置位于站心系的原点,则测站天线指向矢量在站心系下的投影为RwCT,记RwCT为:
Figure BDA0002418106540000341
对天线指向方位角的计算高低角
Figure BDA0002418106540000342
水平角ψ为:
Figure BDA0002418106540000343
Figure BDA0002418106540000344
并根据xCT、yCT、zCT的正负,将高低角
Figure BDA0002418106540000345
水平角ψ划分到对应的角度范围内,完成天线指向方位的预报过程:
Figure BDA0002418106540000346
计算结果为:
Figure BDA0002418106540000347
在本申请的描述中,需要理解的是,术语“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系,仅是为了便于描述本申请和简化描述,而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作,因此不能理解为对本申请的限制。
本领域技术人员知道,除了以纯计算机可读程序代码方式实现本发明提供的***、装置及其各个模块以外,完全可以通过将方法步骤进行逻辑编程来使得本发明提供的***、装置及其各个模块以逻辑门、开关、专用集成电路、可编程逻辑控制器以及嵌入式微控制器等的形式来实现相同程序。所以,本发明提供的***、装置及其各个模块可以被认为是一种硬件部件,而对其内包括的用于实现各种程序的模块也可以视为硬件部件内的结构;也可以将用于实现各种功能的模块视为既可以是实现方法的软件程序又可以是硬件部件内的结构。
以上对本发明的具体实施例进行了描述。需要理解的是,本发明并不局限于上述特定实施方式,本领域技术人员可以在权利要求的范围内做出各种变化或修改,这并不影响本发明的实质内容。在不冲突的情况下,本申请的实施例和实施例中的特征可以任意相互组合。

Claims (12)

1.一种地面测控天线指向卫星的方位角计算方法,其特征在于,包括:
利用星历时刻t0、t0时刻的轨道平根数[a0,e0,i000,M0]和地面测站天线的地理位置信息,经过卫星轨道的相关计算以及多个相关坐标系的转换计算,最终转换为t1时刻在站心系下的天线指向方位角度:高低角
Figure FDA0002418106530000011
水平角ψ,以完成地面测站天线对卫星指向方位的预报过程。
2.如权利要求1所述的一种地面测控天线指向卫星的方位角计算方法,其特征在于,包括:
输入参数无奇异处理步骤:
当轨道偏心率很小,即近圆轨道,e≈0时,为了避免计算中出现奇点,转化为3个无奇点变量,即:
ξ0=e0cos(ω0)
η0=-e0sin(ω0)
λ0=ω0+M0。
3.如权利要求2所述的一种地面测控天线指向卫星的方位角计算方法,其特征在于,还包括:
输入参数归一化处理步骤:
对于递推时间dt进行归一化处理:
dt=t1-t0
Figure FDA0002418106530000012
其中,
t0为轨道递推起始时刻;
t1表示轨道递推结束时刻;
dtn表示递推时间经过归一化处理的结果,下标n表示归一化;
对于半长轴a0进行归一化处理:
Figure FDA0002418106530000021
所述的时间的归一化单位为
Figure FDA0002418106530000022
其中,
Ge为地心引力常数;
Re表示地球赤道半径。
4.如权利要求3所述的一种地面测控天线指向卫星的方位角计算方法,其特征在于,还包括:
摄动项计算步骤:
考虑地球引力场摄动中的J2至J4项摄动,包含计算一阶长期项、一阶短周期项和二阶长期项,各个摄动项的表达式如下:
一阶长期项计算:
Figure FDA0002418106530000023
Figure FDA0002418106530000024
Figure FDA0002418106530000025
其中,
Ω1表示卫星轨道升交点赤经的一阶长周期变化量;
ω1表示卫星轨道近地点幅角的一阶长周期变化量;
λ1表示本文计算所引入的无奇点变量λ的一阶长周期变化量;
所述p为半通径,表达式为p=a×(1-e0 2),轨道平均角速度
Figure FDA0002418106530000026
J2=1.624×10-3
一阶短周期项计算:
Figure FDA0002418106530000027
Figure FDA0002418106530000028
Figure FDA0002418106530000031
Figure FDA0002418106530000032
Figure FDA0002418106530000033
Figure FDA0002418106530000034
其中,
Figure FDA0002418106530000035
表示卫星轨道半长轴的一阶短周期变化量;
Figure FDA0002418106530000036
表示卫星轨道倾角的一阶短周期变化量;
Figure FDA0002418106530000037
表示卫星轨道升交点赤经的一阶短周期变化量;
Figure FDA0002418106530000038
分别表示本文计算所引入的三个无奇点变量ξ、η、λ的一阶短周期变化量;
所述u由一阶长期项计算得到,计算过程如下:
ξz1=ξ0cos(ω1dtn)+η0sin(ω1dtn)
ηz1=η0cos(ω1dtn)-ξ0sin(ω1dtn)
λz1=λ0+(n+λ1)dtn
Figure FDA0002418106530000039
ωz1=arc cos(ξz1/ez1)
Mz1=mod(λz1z1,2π)
Figure FDA00024181065300000310
u=fz1z1
其中,
ξz1、ηz1、λz1分别表示本文计算所引入的三个无奇点变量ξ、η、λ的根据一阶长期变化量的积分结果;
ez1表示卫星轨道偏心率以一阶长周期变化量为积分量的积分结果;
ωz1表示卫星轨道近地点幅角以一阶长周期变化量为积分量的积分结果;
Mz1表示卫星轨道平近点角以一阶长周期变化量为积分量的积分结果;
fz1表示卫星轨道真近点角以一阶长周期变化量为积分量的积分结果;
二阶长期项计算:
Figure FDA0002418106530000041
Figure FDA0002418106530000042
Figure FDA0002418106530000043
其中,
Ω2表示卫星轨道升交点赤经的二阶长周期变化量;
ξ2、λ2分别表示本文计算所引入的无奇点变量ξ、λ的二阶长周期变化量;
J3=2.5356×10-6,J4=7.1022×10-6
5.如权利要求4所述的一种地面测控天线指向卫星的方位角计算方法,其特征在于,还包括:
递推主公式计算步骤:
对于引入的无奇点变量,结合长期项和短周期项摄动项进行解析解构造,忽略长周期项摄动,轨道递推主公式如下:
Figure FDA0002418106530000044
Figure FDA0002418106530000045
Figure FDA0002418106530000046
Figure FDA0002418106530000047
Figure FDA0002418106530000048
Figure FDA0002418106530000049
其中,
αt表示表示t1时刻,卫星轨道半长轴的归一化结果;
is表示t1时刻,卫星轨道倾角;
Ωs表示t1时刻,卫星轨道升交点赤经;
ξs、ηs、λs分别表示t1时刻,计算所引入的三个无奇点变量ξ、η、λ的值。
6.如权利要求5所述的一种地面测控天线指向卫星的方位角计算方法,其特征在于,还包括:
无奇点变量还原步骤:
计算结束后将3个无奇点变量还原
Figure FDA0002418106530000051
ωs=arccos(ξs/es)
Ms=mod(λss,2π) 。
7.如权利要求6所述的一种地面测控天线指向卫星的方位角计算方法,其特征在于,还包括:
归一化变量还原步骤:
将卫星轨道长轴at还原为常规单位:
as=at×Re
所述的as单位:m。
8.如权利要求7所述的一种地面测控天线指向卫星的方位角计算方法,其特征在于,还包括:
惯性系卫星位置计算步骤:
输入t1时刻卫星轨道瞬根数[as,es,isss,Ms],输出卫星的位置矢量在J2000.0坐标系下的分量RwECI,计算方法为:
RwECI=Q*rp
其中,旋转矩阵Q按照3-1-3旋转顺序进行描述:
Figure FDA0002418106530000052
矢量rp
Figure FDA0002418106530000053
其中,M1为真近点角:
Figure FDA0002418106530000061
9.如权利要求8所述的一种地面测控天线指向卫星的方位角计算方法,其特征在于,还包括:
地固系卫星位置计算步骤:
输入目标时刻t1和所述计算出的卫星在惯性系下的位置RwECI,输出t1时刻卫星在地固系下的位置RwECF,具体方法如下:
根据给定的目标时刻t1,计算历元J2000.0(2000年1月1日12时)至给定目标时刻的秒计数值tc,输入t1时刻的年、月、日、时、分、秒,计算儒略日JD:
Figure FDA0002418106530000062
其中,floor()为向下取整运算;
根据儒略日JD计算历元J2000.0至给定目标时刻的秒计数值tc
tc=(JD-2455197.5)×86400+315547200
根据计算得到的历元J2000.0至给定目标时刻的秒计数值tc,计算地球自转矩阵ER、章动矩阵NR以及岁差矩阵PR,计算惯性系到地固系的转换矩阵MECI2ECF
MECI2ECF=ER*NR*PR
并根据所述计算得到的t1时刻卫星在惯性系下的位置RwECI,计算t1时刻卫星在地固系下的位置RwECF
RwECF=MECI2ECF*RwECI。
10.如权利要求9所述的一种地面测控天线指向卫星的方位角计算方法,其特征在于,还包括:
地固系测站天线位置计算步骤:
根据给定的地面测站天线的经纬度和高程,计算地面测站天线在地固系下的位置RtECF,具体方法如下:
输入地面测站天线的经度lon、纬度lat、高程h;
计算坐标分量G1、G2:
Figure FDA0002418106530000071
Figure FDA0002418106530000072
其中,f为地球椭球体几何扁率,f=1/298.257;
计算地面测站天线在地固系下的位置RtECF
Figure FDA0002418106530000073
11.如权利要求10所述的一种地面测控天线指向卫星的方位角计算方法,其特征在于,还包括:
站心系卫星位置计算步骤:
输入计算出的t1时刻卫星在地固系下的位置RwECF和所述计算出的地面测站天线在地固系下的位置RtECF,输出刻卫星在站心系下的位置RwCT,具体方法如下:
在站心系下,计算地固系到站心系的转换矩阵MECF2CT,描述为一次绕地固系的Z轴的旋转和一次绕地固系的X轴的旋转:
MECF2CT=Rx(90°-lat)Rz(90°+lon)
其中,lon为地面测站天线的地理经度;lat为地面测站天线的地理纬度;
Figure FDA0002418106530000074
Figure FDA0002418106530000075
并根据所述t1时刻卫星在地固系下的位置RwECF、地面测站天线在地固系下的位置RtECF,将坐标原点由地心平移到地面测站天线处,并计算t1时刻卫星在站心系下的位置RwCT
RwCT=MECF2CT*(RwECF-RtECF) 。
12.如权利要求11所述的一种地面测控天线指向卫星的方位角计算方法,其特征在于,还包括:
测站天线指向方位角计算步骤:
输入卫星在站心系下的位置,输出t1时刻地面测站天线指向方位角度:高低角
Figure FDA0002418106530000081
水平角ψ,具体方法如下:
高低角、水平角在站心系下定义,在站心系OCTXCTYCTZCT中,OCT表示坐标系原点,XCT表示坐标系的X轴,YCT表示坐标系的Y轴,OCTXCTYCT表示该坐标系的X轴与Y轴所在的平面,高低角
Figure FDA0002418106530000082
为指向矢量RwCT与OCTXCTYCT平面的夹角,定义RwCT矢量与OCTZCT夹角小于90°为正;水平角ψ为指向矢量RwCT在OCTXCTYCT平面的投影与OCTXCT轴的夹角,定义绕OCTZCT轴从OCTXCT轴顺时针转向指向矢量RwCT在OCTXCTYCT面的投影为正,根据此定义求出天线指向方位角。假设地面测站的位置位于站心系的原点,则测站天线指向矢量在站心系下的投影为RwCT,记RwCT为:
Figure FDA0002418106530000083
其中,
xCT、yCT、zCT分别表示测站天线指向矢量RwCT,在站心系OCTXCTYCTZCT中,对应的X轴、Y轴、Z轴的坐标分量;
对天线指向方位角的计算高低角
Figure FDA0002418106530000084
水平角ψ为:
Figure FDA0002418106530000085
Figure FDA0002418106530000086
并根据xCT、yCT、zCT的正负,将高低角
Figure FDA0002418106530000087
水平角ψ划分到对应的角度范围内,完成天线指向方位预报过程:
若zCT≥0,则将高低角
Figure FDA0002418106530000091
划分到范围
Figure FDA0002418106530000092
若xCT≥0,yCT≥0,则将水平角ψ划分到范围
Figure FDA0002418106530000093
若xCT<0,yCT≥0,则将水平角ψ划分到范围
Figure FDA0002418106530000094
若xCT<0,yCT<0,则将水平角ψ划分到范围
Figure FDA0002418106530000095
若xCT≥0,yCT<0,则将水平角ψ划分到范围
Figure FDA0002418106530000096
CN202010197399.2A 2020-03-19 2020-03-19 地面测控天线指向卫星的方位角计算方法 Active CN111427002B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010197399.2A CN111427002B (zh) 2020-03-19 2020-03-19 地面测控天线指向卫星的方位角计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010197399.2A CN111427002B (zh) 2020-03-19 2020-03-19 地面测控天线指向卫星的方位角计算方法

Publications (2)

Publication Number Publication Date
CN111427002A true CN111427002A (zh) 2020-07-17
CN111427002B CN111427002B (zh) 2022-07-12

Family

ID=71548179

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010197399.2A Active CN111427002B (zh) 2020-03-19 2020-03-19 地面测控天线指向卫星的方位角计算方法

Country Status (1)

Country Link
CN (1) CN111427002B (zh)

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112329202A (zh) * 2020-09-30 2021-02-05 北京空间飞行器总体设计部 一种火星车对环绕器天线指向算法的优化实现方法
CN112948741A (zh) * 2021-02-04 2021-06-11 上海卫星工程研究所 一种深空探测器可见弧段计算方法及***
CN114002710A (zh) * 2021-10-20 2022-02-01 上海航天空间技术有限公司 小偏心率低轨卫星的星上轨道位置自主预报方法
CN114002713A (zh) * 2021-10-20 2022-02-01 上海航天空间技术有限公司 一种卫星轨道参数递推处理预报***
CN114137573A (zh) * 2021-09-16 2022-03-04 北京微纳星空科技有限公司 一种卫星测控站及卫星测控方法、设备及存储介质
CN114826438A (zh) * 2022-03-30 2022-07-29 中国空间技术研究院 一种地面测控天线安装位置的确定方法
CN116659543A (zh) * 2023-06-21 2023-08-29 中国人民解放军61540部队 基于遥感卫星轨道根数的卫星位置姿态估计方法和装置

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO1986000863A1 (en) * 1984-07-20 1986-02-13 Hughes Aircraft Company Spin-stabilized satellite with nutation control subsystem
CN102230969A (zh) * 2011-03-22 2011-11-02 航天恒星科技有限公司 一种卫星星座星间链路的长时间自主维持方法
CN104729457A (zh) * 2015-04-16 2015-06-24 哈尔滨工业大学 太阳相对近地轨道微小卫星位置的确定方法
CN105158777A (zh) * 2015-07-31 2015-12-16 上海卫星工程研究所 用于无源测向定位的数据源生成方法
CN107656293A (zh) * 2017-08-30 2018-02-02 清华大学 基于无奇点根数的14参数广播星历卫星位置估计方法
CN108974395A (zh) * 2018-06-21 2018-12-11 中国人民解放军战略支援部队航天工程大学 基于天基激光平台驱动的空间目标变轨计算方法及其装置
CN110816896A (zh) * 2019-10-28 2020-02-21 中国空间技术研究院 一种卫星星上简易轨道外推方法
CN110838864A (zh) * 2018-08-19 2020-02-25 南京理工大学 无人值守的卫星地面站跟踪控制***

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO1986000863A1 (en) * 1984-07-20 1986-02-13 Hughes Aircraft Company Spin-stabilized satellite with nutation control subsystem
CN102230969A (zh) * 2011-03-22 2011-11-02 航天恒星科技有限公司 一种卫星星座星间链路的长时间自主维持方法
CN104729457A (zh) * 2015-04-16 2015-06-24 哈尔滨工业大学 太阳相对近地轨道微小卫星位置的确定方法
CN105158777A (zh) * 2015-07-31 2015-12-16 上海卫星工程研究所 用于无源测向定位的数据源生成方法
CN107656293A (zh) * 2017-08-30 2018-02-02 清华大学 基于无奇点根数的14参数广播星历卫星位置估计方法
CN108974395A (zh) * 2018-06-21 2018-12-11 中国人民解放军战略支援部队航天工程大学 基于天基激光平台驱动的空间目标变轨计算方法及其装置
CN110838864A (zh) * 2018-08-19 2020-02-25 南京理工大学 无人值守的卫星地面站跟踪控制***
CN110816896A (zh) * 2019-10-28 2020-02-21 中国空间技术研究院 一种卫星星上简易轨道外推方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
HEERA SINGH R ET AL.: "Design of Drive Control System for High Efficiency Compact Ground Station Antenna for Remote Sensing Satellites", 《2012 5TH INTERNATIONAL CONFERENCE ON COMPUTERS AND DEVICES FOR COMMUNICATION(CODEC)》 *
晁宁 等: "一种基于卫星天线指向的地面站高度修正方法", 《空间电子技术》 *

Cited By (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112329202A (zh) * 2020-09-30 2021-02-05 北京空间飞行器总体设计部 一种火星车对环绕器天线指向算法的优化实现方法
CN112329202B (zh) * 2020-09-30 2023-08-15 北京空间飞行器总体设计部 一种火星车对环绕器天线指向算法的优化实现方法
CN112948741A (zh) * 2021-02-04 2021-06-11 上海卫星工程研究所 一种深空探测器可见弧段计算方法及***
CN112948741B (zh) * 2021-02-04 2023-02-28 上海卫星工程研究所 一种深空探测器可见弧段计算方法及***
CN114137573A (zh) * 2021-09-16 2022-03-04 北京微纳星空科技有限公司 一种卫星测控站及卫星测控方法、设备及存储介质
CN114002710A (zh) * 2021-10-20 2022-02-01 上海航天空间技术有限公司 小偏心率低轨卫星的星上轨道位置自主预报方法
CN114002713A (zh) * 2021-10-20 2022-02-01 上海航天空间技术有限公司 一种卫星轨道参数递推处理预报***
CN114826438A (zh) * 2022-03-30 2022-07-29 中国空间技术研究院 一种地面测控天线安装位置的确定方法
CN114826438B (zh) * 2022-03-30 2024-05-31 中国空间技术研究院 一种地面测控天线安装位置的确定方法
CN116659543A (zh) * 2023-06-21 2023-08-29 中国人民解放军61540部队 基于遥感卫星轨道根数的卫星位置姿态估计方法和装置
CN116659543B (zh) * 2023-06-21 2024-05-07 中国人民解放军61540部队 基于遥感卫星轨道根数的卫星位置姿态估计方法和装置

Also Published As

Publication number Publication date
CN111427002B (zh) 2022-07-12

Similar Documents

Publication Publication Date Title
CN111427002B (zh) 地面测控天线指向卫星的方位角计算方法
CN107450582B (zh) 一种基于星上实时规划的相控阵数传引导控制方法
US5957982A (en) Method and system for space navigation
CN112629543A (zh) 一种大椭圆轨道及小倾角圆轨道的轨道规划方法
CN111427003A (zh) 地面测站天线对卫星的指向导引***
CN112713922A (zh) 一种多波束通讯卫星的可见性快速预报算法
US20070080858A1 (en) Control segment-based lever-arm correction via curve fitting for high accuracy navigation
CN111427004A (zh) 适用于地面测站天线对卫星指向的坐标转换方法
CN112649006A (zh) 一种太阳同步圆轨道的轨道规划方法
CN111427001A (zh) 适用于地面测站天线对卫星指向的目标定位方法
CN116125503A (zh) 一种高精度卫星轨道确定及预报算法
CN112014869A (zh) 基于天文导航的星间链路自主导航方法及***
CN111679242A (zh) 适用于指向在轨航天器的地面天线导引方法
CN112607056B (zh) 雷达卫星对目标观测自主启动触发方法和***
CN112762925A (zh) 一种基于地磁计和陀螺仪的低轨卫星定姿方法
Barbee et al. Guidance and navigation for rendezvous and proximity operations with a non-cooperative spacecraft at geosynchronous orbit
CN116692028A (zh) 一种小卫星对地快速凝视指向跟踪控制方法及装置
Gustavsson Development of a MatLab based GPS constellation simulation for navigation algorithm developments
CN111366126A (zh) 地面测站天线对卫星指向的视向量计算***
CN112013834B (zh) 基于天文导航的星间链路自主恢复方法及***
CN112394381B (zh) 基于球卫星的全自主月面导航和数据通信方法
CN111427000A (zh) 适用于地面测站天线对卫星指向的目标视向量确定方法
Lincoln et al. Components of a vision assisted constrained autonomous satellite formation flying control system
CN114002710A (zh) 小偏心率低轨卫星的星上轨道位置自主预报方法
Keil Kalman filter implementation to determine orbit and attitude of a satellite in a molniya orbit

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant