CN115061171A - 一种卫星大气摄动测算方法及装置 - Google Patents
一种卫星大气摄动测算方法及装置 Download PDFInfo
- Publication number
- CN115061171A CN115061171A CN202210494763.0A CN202210494763A CN115061171A CN 115061171 A CN115061171 A CN 115061171A CN 202210494763 A CN202210494763 A CN 202210494763A CN 115061171 A CN115061171 A CN 115061171A
- Authority
- CN
- China
- Prior art keywords
- satellite
- perturbation
- calculating
- gravity
- coordinate
- 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.)
- Pending
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO 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
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
- G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
- G01S19/42—Determining position
- G01S19/45—Determining position by combining measurements of signals from the satellite radio beacon positioning system with a supplementary measurement
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
- G01C21/24—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 specially adapted for cosmonautical navigation
-
- Y—GENERAL 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
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
- Y02E60/00—Enabling technologies; Technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Computer Networks & Wireless Communication (AREA)
- Astronomy & Astrophysics (AREA)
- Automation & Control Theory (AREA)
- Navigation (AREA)
Abstract
本发明公开了一种卫星大气摄动测算方法及装置。包括:基于北斗定位和授时功能,获取卫星的卫星坐标、卫星速度和当前时刻;获取卫星的卫星测量姿态信息和太阳帆板转角信息;根据卫星坐标、卫星速度和当前时刻,确定卫星在惯性坐标系下的摄动总和;根据卫星位置和第三方引力场模型,确定卫星对应的非球形引力摄动;根据在当前时刻的太阳坐标和月球坐标,确定卫星对应的日引力摄动和月引力摄动;根据当前时刻、卫星坐标、太阳坐标、卫星测量姿态信息和太阳帆板转角信息,确定卫星对应的太阳光压摄动;根据摄动总和、非球形引力摄动、日引力摄动、月引力摄动、太阳光压摄动,确定卫星的卫星大气摄动。本发明可以提高大气摄动量的计算准确度。
Description
技术领域
本发明涉及卫星大气摄动测算技术领域,特别是一种卫星大气摄动测算方法及装置。
背景技术
低轨道卫星(轨道高度100-1500km)在空间中的轨道摄动主要有非球形引力摄动和大气摄动,其次受到日月引力摄动和太阳光压力摄动。轨道预报需要对各种摄动进行准确预估,而非球形引力、日月引力、太阳光压力均能较好的求解:对于非球形引力,利用高阶地球引力场模型可以精确求解;对于日月引力,利用历表或日月轨道参数也容易确定;对于太阳光压力,只需确定光照条件、卫星面积-质量比、卫星表面反射系数即可求解。大气摄动最为复杂,与摄动加速度直接相关的大气密度、大气阻力系数、大气运动速度均存在极大的不确定性:
1)大气密度ρ受多种空间环境(如地磁活动、太阳活动)影响剧烈,还存在周日项、半周日项等因素,只能通过经验公式估计,难以得到准确值。
2)大气阻力系数Cd大小难以准确估计,受航天器形状、表面材料影响,受大气成分影响而实时变化,但工程上通常取常值2.2,给大气阻力摄动计算造成了误差。
3)高空(100km以上)大气运动规律复杂,机制尚不清楚。轨道摄动分析中通常予以简化,认为大气跟随地球同步旋转,也给卫星和大气的相对运动速度计算带来误差。
发明内容
本发明解决的技术问题是:克服现有技术的不足,提供了一种卫星大气摄动测算方法及装置。
本发明的技术解决方案是:
第一方面,本发明实施例提供了一种卫星大气摄动测算方法,包括:
基于北斗定位和授时功能,获取卫星的卫星坐标、卫星速度和当前时刻;
获取所述卫星的卫星测量姿态信息和太阳帆板转角信息;
根据所述卫星坐标、所述卫星速度和所述当前时刻,确定所述卫星在惯性坐标系下的摄动总和;
根据所述卫星位置和第三方引力场模型,确定所述卫星对应的非球形引力摄动;
根据在所述当前时刻的太阳坐标和月球坐标,确定所述卫星对应的日引力摄动和月引力摄动;
根据所述当前时刻、所述卫星坐标、所述太阳坐标、所述卫星测量姿态信息和所述太阳帆板转角信息,确定所述卫星对应的太阳光压摄动;
根据所述摄动总和、所述非球形引力摄动、所述日引力摄动、所述月引力摄动、所述太阳光压摄动,确定所述卫星的卫星大气摄动。
可选地,所述根据所述卫星坐标、所述卫星速度和所述当前时刻,确定所述卫星在惯性坐标系下的摄动总和,包括:
获取所述当前时刻的J2000历元;
根据所述当前时刻的J2000历元,计算得到地固系与J2000地心平赤道之间的坐标系转移矩阵;
根据所述坐标系转移矩阵、所述卫星坐标对应的位置矢量和所述卫星速度对应的速度矢量,计算得到所述卫星在J2000地心平赤道坐标系下的目标位置矢量和目标速度矢量;
根据所述目标位置矢量和所述目标速度矢量,计算得到所述卫星的卫星轨道六根数;
根据所述卫星轨道六根数和所述当前时刻,计算得到所述卫星在二体运动下的卫星理论位置;
根据所述卫星理论位置、所述卫星在地固系下的卫星位置和所述当前时刻,计算得到所述卫星在惯性坐标系下的摄动总和。
可选地,所述根据所述卫星位置和第三方引力场模型,确定所述卫星对应的非球形引力摄动,包括:
根据所述卫星位置,计算得到所述卫星在二体运动下的地球引力、及所述卫星对应的卫星经度、卫星纬度和地心距;
根据所述卫星经度、所述卫星纬度、所述地心距和所述第三方引力场模型,计算得到所述卫星对应的基准地球引力;
根据所述地球引力和所述基准地球引力,计算得到所述卫星对应的非球形引力摄动。
可选地,所述根据在所述当前时刻的太阳坐标和月球坐标,确定所述卫星对应的日引力摄动和月引力摄动,包括:
根据所述当前时刻,计算得到太阳在当前时刻的太阳坐标,及月球在当前时刻的月球坐标;
根据所述卫星在二体运动下的卫星理论位置、所述太阳坐标和所述月球坐标,计算得到所述卫星对应的日引力摄动和月引力摄动。
可选地,所述根据所述摄动总和、所述非球形引力摄动、所述日引力摄动、所述月引力摄动、所述太阳光压摄动,确定所述卫星的卫星大气摄动,包括:
计算得到所述日引力摄动、所述月引力摄动和所述太阳光压摄动之间的摄动和值;
计算所述摄动总和与所述摄动和值之间的差值,并将该差值作为所述卫星大气摄动。
第二方面,本发明实施例提供了一种卫星大气摄动测算装置,包括:
卫星信息获取模块,用于基于北斗定位和授时功能,获取卫星的卫星坐标、卫星速度和当前时刻;
卫星姿态获取模块,用于获取所述卫星的卫星测量姿态信息和太阳帆板转角信息;
摄动总和确定模块,用于根据所述卫星坐标、所述卫星速度和所述当前时刻,确定所述卫星在惯性坐标系下的摄动总和;
非球形引力摄动确定模块,用于根据所述卫星位置和第三方引力场模型,确定所述卫星对应的非球形引力摄动;
日月引力摄动确定模块,用于根据在所述当前时刻的太阳坐标和月球坐标,确定所述卫星对应的日引力摄动和月引力摄动;
太阳光压摄动确定模块,用于根据所述当前时刻、所述卫星坐标、所述太阳坐标、所述卫星测量姿态信息和所述太阳帆板转角信息,确定所述卫星对应的太阳光压摄动;
卫星大气摄动确定模块,用于根据所述摄动总和、所述非球形引力摄动、所述日引力摄动、所述月引力摄动、所述太阳光压摄动,确定所述卫星的卫星大气摄动。
可选地,所述摄动总和确定模块包括:
J2000历元获取单元,用于获取所述当前时刻的J2000历元;
坐标系转移矩阵计算单元,用于根据所述当前时刻的J2000历元,计算得到地固系与J2000地心平赤道之间的坐标系转移矩阵;
目标位置速度矢量计算单元,用于根据所述坐标系转移矩阵、所述卫星坐标对应的位置矢量和所述卫星速度对应的速度矢量,计算得到所述卫星在 J2000地心平赤道坐标系下的目标位置矢量和目标速度矢量;
卫星轨道六根数计算单元,用于根据所述目标位置矢量和所述目标速度矢量,计算得到所述卫星的卫星轨道六根数;
卫星理论位置计算单元,用于根据所述卫星轨道六根数和所述当前时刻,计算得到所述卫星在二体运动下的卫星理论位置;
摄动总和计算单元,用于根据所述卫星理论位置、所述卫星在地固系下的卫星位置和所述当前时刻,计算得到所述卫星在惯性坐标系下的摄动总和。
可选地,所述非球形引力摄动确定模块包括:
地球引力计算单元,用于根据所述卫星位置,计算得到所述卫星在二体运动下的地球引力、及所述卫星对应的卫星经度、卫星纬度和地心距;
基准地球引力计算单元,用于根据所述卫星经度、所述卫星纬度、所述地心距和所述第三方引力场模型,计算得到所述卫星对应的基准地球引力;
非球形引力摄动计算单元,用于根据所述地球引力和所述基准地球引力,计算得到所述卫星对应的非球形引力摄动。
可选地,所述日月引力摄动确定模块包括:
日月坐标计算单元,用于根据所述当前时刻,计算得到太阳在当前时刻的太阳坐标,及月球在当前时刻的月球坐标;
日月引力摄动计算单元,用于根据所述卫星在二体运动下的卫星理论位置、所述太阳坐标和所述月球坐标,计算得到所述卫星对应的日引力摄动和月引力摄动。
可选地,所述卫星大气摄动确定模块包括:
摄动和值计算单元,用于计算得到所述日引力摄动、所述月引力摄动和所述太阳光压摄动之间的摄动和值;
卫星大气摄动获取单元,用于计算所述摄动总和与所述摄动和值之间的差值,并将该差值作为所述卫星大气摄动。
本发明与现有技术相比的优点在于:
1、利用北斗系统的授时、导航功能,通过联合求解多种摄动要素,避免了大气不确定性造成的计算误差,实现了大气摄动量的准确计算;
2、利用北斗系统的短报文通信功能,在数据获取环节摆脱了对卫星遥测通道的依赖,实现全弧段、全天候数据获取,实现摄动量的实时计算;
3、通过设计运算平台软件,实现了轨道大气摄动量的全流程自动求解。
附图说明
图1为本发明实施例提供的一种卫星大气摄动测算方法的步骤流程图;
图2为本发明实施例提供的一种卫星大气摄动测算系统的示意图;
图3为本发明实施例提供的一种运算平台内大气摄动计算流程的示意图;
图4为本发明实施例提供的一种运算平台的示意图;
图5为本发明实施例提供的一种赤道惯性坐标系下卫星轨道要素的示意图;
图6为本发明实施例提供的一种卫星大气摄动测算装置的结构示意图。
具体实施方式
实施例一
参照图1,示出了本发明实施例提供的一种卫星大气测算方法的步骤流程图,如图1所示,该卫星大气测算方法可以包括以下步骤:
步骤101:基于北斗定位和授时功能,获取卫星的卫星坐标、卫星速度和当前时刻。
本发明实施例以卫星通过北斗导航获取的坐标、速度、时刻以及卫星自身测量姿态等信息驱动,同时利用北斗系统的短报文通信功能保证信息的实时回传,构建软件平台对卫星所受的大气摄动量进行全流程计算。
本发明的测试方法原理如图2所示,共包含6条数据流,测试过程按照数据流的编号顺序进行,具体过程为:
数据流1为通过北斗终端和卫星姿态控制系统获取的四项初始信息:卫星坐标、卫星速度、当前时刻、卫星测量姿态与太阳帆板转角,该四项信息是数据流2至6的输入。星载北斗导航接收机获取卫星坐标、卫星速度、当前时刻三项信息,卫星姿态控制系统获取卫星测量姿态与太阳帆板转角信息,星载北斗短报文终端将该四项信息传输给地面北斗短报文终端。
数据流2输入为卫星坐标、卫星速度、当前时刻,输出为计算卫星所受的摄动总和,作为数据流6的输入。
数据流3输入为卫星坐标,输出为卫星所受的非球形引力摄动,作为数据流6的输入。
数据流4输入为卫星坐标、当前时刻,输出为卫星所受的日、月引力摄动,作为数据流6的输入。
数据流5输入为卫星坐标、当前时刻、卫星测量姿态与太阳帆板转角,输出为卫星所受的太阳光压力摄动,作为数据流6的输入。
数据流6输入为卫星摄动总和、非球形引力摄动、日月引力摄动、太阳光压力摄动(即数据流2、3、4、5的输出),输出为卫星的大气摄动。
本实施例提供的测试方法原理可以如图2所示,分为数据获取和运算平台两部分,其中,数据获取部分可以包括:星载北斗接收机和卫星姿态控制系统,其中,星载北斗接收机可以利用北斗定位和授时功能,获取卫星的卫星坐标、卫星速度和当前时刻。
步骤102:获取所述卫星的卫星测量姿态信息和太阳帆板转角信息。
卫星姿态控制系统可以获取卫星测量姿态信息和太阳帆板转角信息。
如图2所示,在通过星载北斗接收机获取卫星的卫星坐标、卫星速度和当前时刻,并通过卫星姿态控制系统获取卫星测量姿态信息和太阳帆板转角信息之后,可以由星载北斗接收机将卫星坐标、卫星速度和当前时刻发送给星载北斗短报文终端,并由卫星姿态控制系统将卫星测量姿态信息和太阳帆板转角信息发送给星载北斗短报文终端,进而,可以由星载北斗短报文终端通过短报文发送给地面北斗短报文终端,该地面北斗短报文终端可以将这些信息发送给运算平台。
在本示例中,卫星坐标、速度的描述均为3个变量,姿态四元数为4个变量,当前时刻与帆板转角均为1个变量,为保证足够的精度,可取卫星坐标、速度、姿态、帆板转角均占用4字节,当前时刻占用8字节,共占用52字节。而北斗短报文的区域通信能力为每次1750字节,全球通信能力为每次70字节,单次短报文容量完全满足数据传输需求。
本发明实施例利用北斗系统的短报文通信功能保证信息的实时回传,可以实现卫星大气摄动的实时测算。
在获取到卫星的卫星坐标、卫星速度和当前时刻,以及卫星的卫星测量姿态信息和太阳帆板转角信息之后,执行步骤103。
步骤103:根据所述卫星坐标、所述卫星速度和所述当前时刻,确定所述卫星在惯性坐标系下的摄动总和。
在获取到卫星的卫星坐标、卫星速度和当前时刻之后,可以根据卫星的卫星坐标、卫星速度和当前时刻确定出卫星在惯性坐标系下的摄动总和,具体地,可以结合下述具体实现方式进行详细描述。
首先,对本发明实施例涉及的一些公式进行描述。
如图5所示,Υ-春分点,N-卫星轨道的升交点,S-卫星位置,P-卫星轨道的近地点,Ω-升交点赤经,ω-近地点幅角,f—真近点角,i-轨道倾角,e- 偏心率矢量,从地心指向近地点,长度为e,W-轨道平面法线的单位矢量,r- 卫星向径。
在本发明的一种具体实现方式中,上述步骤103可以包括:
子步骤A1:获取所述当前时刻的J2000历元。
在本发明实施例中,在运算平台接收到卫星的上述信息之后,可以进行相应的计算,该计算流程可以如图2右侧部分、图3和图4所示。
可以获取当前时刻的J2000历元,具体地,自标准历元J2000起算,当前时刻的儒略世纪数,其计算公式为:
其中JD(t0)为t0时刻对应的儒略日。
子步骤A2:根据所述当前时刻的J2000历元,计算得到地固系与J2000 地心平赤道之间的坐标系转移矩阵。
在获取到当前时刻的J2000历元之后,可以根据当前时刻的J2000历元计算得到地固系与J2000地心平赤道之间的坐标系转移矩阵,具体地,可以参照下述公式的计算过程。
地固系下的矢量或坐标p到惯性系p’的表达式为:
p'=(PR)T(NR)T(ER)T(EP)Tp=(IR)p (2-2)
上式中,EP为极移,ER地球自转方向,NR为章动,PR为岁差,各个矩阵的表达式
如下:
EP=Ry(-xp)Rx(-yp) (2-3)
ER=Rz(SG) (2-4)
NR=Rx(-Δε)Ry(Δθ)Rz(-Δμ) (2-5)
PR=Rz(-zA)Ry(θA)Rz(-ξA) (2-6)
其中Rx(θ)、Ry(θ)、Rz(θ)为坐标转移矩阵,其表达式为:
由于极移一般不超过0".8,计算时忽略其影响,式(5)中,SG为格林尼治真恒星时角,其表达式为:
岁差模型采用IAU1976,式(2-6)中,ξA、zA、θA为岁差角,其表达式如下:
章动包括黄经章动ΔΨ、赤经章动Δμ、赤纬章动Δθ、交角章动Δε,章动序列采用IAU1980,与日月引力有关的项共106项,其表达式如下:
式(2-11)序列中,αi(t)所表示的五项基本变量分别为月球平近点角l,太阳平近点角l’,月球平升交角据F,月球平角距D和月球轨道升交点黄经Ω,其表达式如下:
相应的赤经和赤纬章动Δμ和Δθ为:
Δμ=Δψcosε (2-13)
Δθ=Δψsinε (2-14)
其中黄赤交角的计算公式如下:
ε=84381.448”-46.8150”t-0.00059”t2+0.001813”t3 (2-15)
式(2-9)至式(2-15)中,t由(2-1)所得。
子步骤A3:根据所述坐标系转移矩阵、所述卫星坐标对应的位置矢量和所述卫星速度对应的速度矢量,计算得到所述卫星在J2000地心平赤道坐标系下的目标位置矢量和目标速度矢量。
目标位置矢量是指卫星在J2000地心平赤道坐标系下的位置矢量。
目标速度矢量是指卫星在J2000地心平赤道坐标系下的速度矢量。
在计算得到坐标系转移矩阵之后,可以根据坐标系转移矩阵、卫星坐标对应的位置矢量和卫星速度对应的速度矢量,计算得到卫星在J2000地心平赤道坐标系下的目标位置矢量和目标速度矢量,具体地计算过程可以参照下述公式的描述。
vi=(IR)g(ve+ve0) (2-16)
ri=(IR)·re (2-17)
ve0为地固系随地球自转形成的速度分量,有
子步骤A4:根据所述目标位置矢量和所述目标速度矢量,计算得到所述卫星的卫星轨道六根数。
在计算得到目标位置矢量和目标速度矢量之后,可以根据目标位置矢量和目标速度矢量计算得到卫星的卫星轨道六根数,具体地计算过程可以参照下述公式的描述。
卫星轨道由六根数a、e、i、Ω、ω、f来描述,六根数定义如图3所示,需要2-c计算得到的位置速度矢量求得轨道根数。
在J2000地心平赤道坐标系下,卫星动量矩为h=ri×vi。
轨道升交点赤经也由与动量矩决定,为Ω=arccos(Ngux),式中 N=(uz×h)/|uz×h|为轨道节线单位矢量,ux=(0,0,1)为+X坐标轴的单位向量。
近点角也可用偏近点角或平近点角表达,真近点角与偏近点角E的关系为
由偏近点角可得平近点角
M=E-esinE (2-20)
子步骤A5:根据所述卫星轨道六根数和所述当前时刻,计算得到所述卫星在二体运动下的卫星理论位置。
在计算到卫星轨道六根数之后,可以根据卫星轨道六根数、当前时刻计算得到卫星在二体运动下的卫星理论位置(即卫星的实际位置),具体地计算过程可以参照下述公式的描述。
二体运动将地球和卫星均当做质点。设一帧输入数据对应的时刻为t0,下帧数据对应的时刻为t1。以t0时刻的卫星定位数据为起始点(输入),不计一切摄动,按二体运动规律推算t1时刻的理论卫星位置。
按上述步骤计算的卫星在t0时刻的轨道参数,记为,a、e、i、Ω、ω、f0, E0,M0。则二体运动下t1时刻轨道参数为a、e、i、Ω、ω、fb,Eb,Mb,可见五个轨道参数保持不变,近点角计算方法如下。
平近点角Mb为均匀变化:
偏近点角Eb由平近点角多次迭代求得:
E1=Mb+esinE0
E2=Mb+esinE1
L
En=Mb+esinEn-1
Eb=Mb+esinEn (2-22)
真近点角fb由偏近点角按公式(2-19)推导,得到真近点角后,可推导出卫星在J2000地心平赤道坐标系的理论坐标rib为:
子步骤A6:根据所述卫星理论位置、所述卫星在地固系下的卫星位置和所述当前时刻,计算得到所述卫星在惯性坐标系下的摄动总和。
在计算得到卫星的卫星理论位置之后,可以根据卫星理论位置、卫星在地固系下的卫星位置和当前时刻计算得到卫星在惯性坐标系下的摄动总和,具体地计算过程可以参照下述公式的描述。
记t1时刻由北斗系统获取的卫星实际位置为ri1=[xi1,yi1,zi1]T
记惯性坐标系下的摄动加速度总和为Wall=[xall,yall,zall]T,则有
步骤104:根据所述卫星位置和第三方引力场模型,确定所述卫星对应的非球形引力摄动。
在获取到卫星位置之后,可以根据卫星位置和第三方引力场模型确定出卫星对应的非球形引力摄动,具体地,可以结合下述具体实现方式的描述。
在本发明的另一种具体实现方式中,上述步骤104可以包括:
子步骤B1:根据所述卫星位置,计算得到所述卫星在二体运动下的地球引力、及所述卫星对应的卫星经度、卫星纬度和地心距。
在本发明实施例中,在获取到卫星的卫星位置之后,可以根据卫星位置计算得到卫星在二体运动下的地球引力、及卫星对应的卫星经度、卫星纬度和地心距,具体地,可以参照下述公式的描述。
1、计算二体运动下地球引力
2、计算卫星经度、纬度、地心距
子步骤B2:根据所述卫星经度、所述卫星纬度、所述地心距和所述第三方引力场模型,计算得到所述卫星对应的基准地球引力。
在计算得到卫星经度、卫星纬度、地心距和地球引力之后,可以根据卫星经度、卫星纬度、地心距和地球引力,结合第三方引力场模型,计算得到卫星对应的基准地球引力,即卫星的精确地球引力,具体地计算过程可以参照下述公式的描述。
根据归一化的天体引力场位函数、引力场模型,卫星受到的精确地球引力可写成下列形式:
子步骤B3:根据所述地球引力和所述基准地球引力,计算得到所述卫星对应的非球形引力摄动。
在获取到地球引力和基准地球引力之后,可以根据地球引力和基准地球引力计算得到卫星对应的非球形引力摄动,具体地计算过程可以参照下述公式的描述。
根据上述步骤中的引力场函数求偏导数,得到卫星所受精确地球引力Wa在径向、东向、北向的分量:
涉及勒让德多项式的导数,其递推公式为:
卫星所受地球引力摄动为Wga=Wa-Ga,其径向、东向、北向分量为:
步骤105:根据在所述当前时刻的太阳坐标和月球坐标,确定所述卫星对应的日引力摄动和月引力摄动。
在获取到当前时刻之后,可以获取当前时刻的太阳坐标和月球坐标,并根据太阳坐标和月球坐标确定出卫星对应的日引力摄动和月引力摄动,具体地,可以结合下述具体实现方式进行详细描述。
在本发明的另一种具体实现方式中,上述步骤105可以包括:
子步骤C1:根据所述当前时刻,计算得到太阳在当前时刻的太阳坐标,及月球在当前时刻的月球坐标。
子步骤C2:根据所述卫星在二体运动下的卫星理论位置、所述太阳坐标和所述月球坐标,计算得到所述卫星对应的日引力摄动和月引力摄动。
在本发明实施例中,计算日月坐标有两种方法:用月球、地球的平均轨道参数求取近似坐标或读取天文历表。前者优点是软件实现时占用的空间极小,后者的优点是更为精确。此处选取平均轨道参数法。
采用的月球平均轨道根数中,半长轴、偏心率、倾角、升交点黄经计算公式为:
al=383397.7916+0.0038t
el=0.055545526-0.000000016t
il=5°.15668983-0″.00008t+0″.02966t2-0″.000042t3-0″.00000013t4
Ωl=125°.04455504-6967919″.3622t+6″.3622t2+0″.007625t3-0″.00003586t4
月球近地点幅角、平近点角计算公式为:
λl=218°.31665436+1732559343″.73604t-5″.8883t2+0″.006604t3-0″.00003169t4
参考系为J2000地心黄道坐标系,由轨道根数求解月球在该坐标系下的坐标参照公式(2-19)(2-22)(2-23),得月球坐标为rl0=(xl0,yl0,zl0)T。
J2000地心平赤道坐标系下月球坐标为rl=(xl,yl,zl)T=rl0Rx(-ε0),式中ε0为J2000.0时刻平黄赤交角,有ε0=23°26′21″.448,Rx函数见式2-7。
地球绕日平均轨道根数中,坐标系为J2000日心黄道坐标系,采用的半长轴、偏心率、倾角、升交点黄经计算公式为:
ae=1.0000010178au
ee=0.0167086342-0.0004203654t-0.0000126734t2+144410-10t3-210-10t4+310- 10t5
ie=469″.97289t-3″.35053t2-0″.12374t3+0″.00027t4-0″.00001t5+0″.00001t6
Ωe=174°.87317577-8679″.27034t+15″.34191t2+0″.00532t3-0″.03734t4-0″.00073t5+0″.00004t6
其中,au为天文单位,1au=1.495978707108km。
近地点幅角、平近点角计算公式为:
λe=100°.46645683+1295977422″.83429t-2″.04411t2-0″.00523t3
由轨道根数求解地球在J2000日心黄道坐标系中的坐标参照公式(2-19) (2-22)(2-23),得地球坐标为re0=(xe0,ye0,ze0)T。J2000地心黄道坐标系下太阳坐标为rs0=(xs0,ys0,zs0)T=-re0。J2000地心平赤道坐标系下太阳坐标为
rs=(xs,ys,zs)T=rs0Rx(-ε0)。
根据上述步骤,有矢量为:
Vrs=ri-rs,Vrl=ri-rl (2-30)
太阳引力摄动为:
月球引力摄动为:
G为万有引力常数,ms为太阳质量,ml为月球质量。日月引力摄动为:
Wsl=Ws+Wl (2-33)
步骤106:根据所述当前时刻、所述卫星坐标、所述太阳坐标、所述卫星测量姿态信息和所述太阳帆板转角信息,确定所述卫星对应的太阳光压摄动。
在获取到卫星的相关参数之后,可以根据当前时刻、卫星坐标、太阳坐标、卫星测量姿态信息和太阳帆板转角信息确定出卫星对应的太阳光压摄动,加推的:
卫星相对J2000地心平赤道坐标系的姿态四元数为q=[q1,q2,q3,q4]T,太阳帆板转角为ap为通过北斗短报文获取的卫星信息。太阳帆板的面积为Sp,反射系数为ηp,为已知量。设帆板转动轴与星体+Y平行,且帆板转角为0时,其法向与星体-Z轴平行。
卫星的姿态矩阵为:
J2000地心平赤道坐标系下,卫星太阳帆板法向的单位向量为:
太阳光源方向的单位矢量为:
帆板法向与太阳光源方向夹角(阳光入射角)为:
太阳辐射压强取ρf=4.5×10-6N/m。
帆板所受太阳光压力为:
星体的近似等效截面积为:
卫星太阳帆板受照面积大且保持对准太阳,吸收和反射的太阳能量较高,因此作为太阳光压摄动的主要部分。而星体的照射情况受自身形状、帆板遮挡、自身姿态等影响较大,且受照面积通常明显小于帆板,对其光压进行简化估算。低轨卫星太阳光压摄动量级通常为10-8,远低于非球形引力(10-3)、日月引力 (10-7)、大气阻力(10-6),且星体受到的光压在总光压里又不占主要部分,因此这种估算是可以接受的。
如果已知卫星质量为ms,则卫星收到的太阳光压摄动为
在具体实现中,仅在卫星受到光照时进行此项计算,如果卫星处于地球或月球阴影中,则太阳光压摄动Wf=0。
步骤107:根据所述摄动总和、所述非球形引力摄动、所述日引力摄动、所述月引力摄动、所述太阳光压摄动,确定所述卫星的卫星大气摄动。
在通过上述步骤得到摄动总和、非球形引力摄动、日引力摄动、月引力摄动、太阳光压摄动之后,可以根据摄动总和、非球形引力摄动、日引力摄动、月引力摄动和太阳光压摄动确定出卫星的卫星大气摄动,具体地,可以计算得到日引力摄动、月引力摄动和太阳光压摄动之间的摄动和值,然后计算得到摄动总和与摄动和值之间的差值,并将该差值作为卫星大气摄动。具体地计算过程可以下述公式所示:
WA=Wall-Wga-Wsl-Wf (2-37)
本发明实施例利用北斗系统的授时、导航功能,通过联合求解多种摄动要素,避免了大气不确定性造成的计算误差,实现了大气摄动量的准确计算。
实施例二
参照图6,示出了本发明实施例提供的一种卫星大气测算装置的结构示意图,如图6所示,该卫星大气测算装置可以包括以下模块:
卫星信息获取模块610,用于基于北斗定位和授时功能,获取卫星的卫星坐标、卫星速度和当前时刻;
卫星姿态获取模块620,用于获取所述卫星的卫星测量姿态信息和太阳帆板转角信息;
摄动总和确定模块630,用于根据所述卫星坐标、所述卫星速度和所述当前时刻,确定所述卫星在惯性坐标系下的摄动总和;
非球形引力摄动确定模块640,用于根据所述卫星位置和第三方引力场模型,确定所述卫星对应的非球形引力摄动;
日月引力摄动确定模块650,用于根据在所述当前时刻的太阳坐标和月球坐标,确定所述卫星对应的日引力摄动和月引力摄动;
太阳光压摄动确定模块660,用于根据所述当前时刻、所述卫星坐标、所述太阳坐标、所述卫星测量姿态信息和所述太阳帆板转角信息,确定所述卫星对应的太阳光压摄动;
卫星大气摄动确定模块670,用于根据所述摄动总和、所述非球形引力摄动、所述日引力摄动、所述月引力摄动、所述太阳光压摄动,确定所述卫星的卫星大气摄动。
可选地,所述摄动总和确定模块包括:
J2000历元获取单元,用于获取所述当前时刻的J2000历元;
坐标系转移矩阵计算单元,用于根据所述当前时刻的J2000历元,计算得到地固系与J2000地心平赤道之间的坐标系转移矩阵;
目标位置速度矢量计算单元,用于根据所述坐标系转移矩阵、所述卫星坐标对应的位置矢量和所述卫星速度对应的速度矢量,计算得到所述卫星在 J2000地心平赤道坐标系下的目标位置矢量和目标速度矢量;
卫星轨道六根数计算单元,用于根据所述目标位置矢量和所述目标速度矢量,计算得到所述卫星的卫星轨道六根数;
卫星理论位置计算单元,用于根据所述卫星轨道六根数和所述当前时刻,计算得到所述卫星在二体运动下的卫星理论位置;
摄动总和计算单元,用于根据所述卫星理论位置、所述卫星在地固系下的卫星位置和所述当前时刻,计算得到所述卫星在惯性坐标系下的摄动总和。
可选地,所述非球形引力摄动确定模块包括:
地球引力计算单元,用于根据所述卫星位置,计算得到所述卫星在二体运动下的地球引力、及所述卫星对应的卫星经度、卫星纬度和地心距;
基准地球引力计算单元,用于根据所述卫星经度、所述卫星纬度、所述地心距和所述第三方引力场模型,计算得到所述卫星对应的基准地球引力;
非球形引力摄动计算单元,用于根据所述地球引力和所述基准地球引力,计算得到所述卫星对应的非球形引力摄动。
可选地,所述日月引力摄动确定模块包括:
日月坐标计算单元,用于根据所述当前时刻,计算得到太阳在当前时刻的太阳坐标,及月球在当前时刻的月球坐标;
日月引力摄动计算单元,用于根据所述卫星在二体运动下的卫星理论位置、所述太阳坐标和所述月球坐标,计算得到所述卫星对应的日引力摄动和月引力摄动。
可选地,所述卫星大气摄动确定模块包括:
摄动和值计算单元,用于计算得到所述日引力摄动、所述月引力摄动和所述太阳光压摄动之间的摄动和值;
卫星大气摄动获取单元,用于计算所述摄动总和与所述摄动和值之间的差值,并将该差值作为所述卫星大气摄动。
本申请所述具体实施方式可以使本领域的技术人员更全面地理解本申请,但不以任何方式限制本申请。因此,本领域技术人员应当理解,仍然对本申请进行修改或者等同替换;而一切不脱离本申请的精神和技术实质的技术方案及其改进,均应涵盖在本申请专利的保护范围中。
本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。
Claims (10)
1.一种卫星大气摄动测算方法,其特征在于,包括:
基于北斗定位和授时功能,获取卫星的卫星坐标、卫星速度和当前时刻;
获取所述卫星的卫星测量姿态信息和太阳帆板转角信息;
根据所述卫星坐标、所述卫星速度和所述当前时刻,确定所述卫星在惯性坐标系下的摄动总和;
根据所述卫星位置和第三方引力场模型,确定所述卫星对应的非球形引力摄动;
根据在所述当前时刻的太阳坐标和月球坐标,确定所述卫星对应的日引力摄动和月引力摄动;
根据所述当前时刻、所述卫星坐标、所述太阳坐标、所述卫星测量姿态信息和所述太阳帆板转角信息,确定所述卫星对应的太阳光压摄动;
根据所述摄动总和、所述非球形引力摄动、所述日引力摄动、所述月引力摄动、所述太阳光压摄动,确定所述卫星的卫星大气摄动。
2.根据权利要求1所述的方法,其特征在于,所述根据所述卫星坐标、所述卫星速度和所述当前时刻,确定所述卫星在惯性坐标系下的摄动总和,包括:
获取所述当前时刻的J2000历元;
根据所述当前时刻的J2000历元,计算得到地固系与J2000地心平赤道之间的坐标系转移矩阵;
根据所述坐标系转移矩阵、所述卫星坐标对应的位置矢量和所述卫星速度对应的速度矢量,计算得到所述卫星在J2000地心平赤道坐标系下的目标位置矢量和目标速度矢量;
根据所述目标位置矢量和所述目标速度矢量,计算得到所述卫星的卫星轨道六根数;
根据所述卫星轨道六根数和所述当前时刻,计算得到所述卫星在二体运动下的卫星理论位置;
根据所述卫星理论位置、所述卫星在地固系下的卫星位置和所述当前时刻,计算得到所述卫星在惯性坐标系下的摄动总和。
3.根据权利要求1所述的方法,其特征在于,所述根据所述卫星位置和第三方引力场模型,确定所述卫星对应的非球形引力摄动,包括:
根据所述卫星位置,计算得到所述卫星在二体运动下的地球引力、及所述卫星对应的卫星经度、卫星纬度和地心距;
根据所述卫星经度、所述卫星纬度、所述地心距和所述第三方引力场模型,计算得到所述卫星对应的基准地球引力;
根据所述地球引力和所述基准地球引力,计算得到所述卫星对应的非球形引力摄动。
4.根据权利要求2所述的方法,其特征在于,所述根据在所述当前时刻的太阳坐标和月球坐标,确定所述卫星对应的日引力摄动和月引力摄动,包括:
根据所述当前时刻,计算得到太阳在当前时刻的太阳坐标,及月球在当前时刻的月球坐标;
根据所述卫星在二体运动下的卫星理论位置、所述太阳坐标和所述月球坐标,计算得到所述卫星对应的日引力摄动和月引力摄动。
5.根据权利要求1所述的方法,其特征在于,所述根据所述摄动总和、所述非球形引力摄动、所述日引力摄动、所述月引力摄动、所述太阳光压摄动,确定所述卫星的卫星大气摄动,包括:
计算得到所述日引力摄动、所述月引力摄动和所述太阳光压摄动之间的摄动和值;
计算所述摄动总和与所述摄动和值之间的差值,并将该差值作为所述卫星大气摄动。
6.一种卫星大气摄动测算装置,其特征在于,包括:
卫星信息获取模块,用于基于北斗定位和授时功能,获取卫星的卫星坐标、卫星速度和当前时刻;
卫星姿态获取模块,用于获取所述卫星的卫星测量姿态信息和太阳帆板转角信息;
摄动总和确定模块,用于根据所述卫星坐标、所述卫星速度和所述当前时刻,确定所述卫星在惯性坐标系下的摄动总和;
非球形引力摄动确定模块,用于根据所述卫星位置和第三方引力场模型,确定所述卫星对应的非球形引力摄动;
日月引力摄动确定模块,用于根据在所述当前时刻的太阳坐标和月球坐标,确定所述卫星对应的日引力摄动和月引力摄动;
太阳光压摄动确定模块,用于根据所述当前时刻、所述卫星坐标、所述太阳坐标、所述卫星测量姿态信息和所述太阳帆板转角信息,确定所述卫星对应的太阳光压摄动;
卫星大气摄动确定模块,用于根据所述摄动总和、所述非球形引力摄动、所述日引力摄动、所述月引力摄动、所述太阳光压摄动,确定所述卫星的卫星大气摄动。
7.根据权利要求6所述的装置,其特征在于,所述摄动总和确定模块包括:
J2000历元获取单元,用于获取所述当前时刻的J2000历元;
坐标系转移矩阵计算单元,用于根据所述当前时刻的J2000历元,计算得到地固系与J2000地心平赤道之间的坐标系转移矩阵;
目标位置速度矢量计算单元,用于根据所述坐标系转移矩阵、所述卫星坐标对应的位置矢量和所述卫星速度对应的速度矢量,计算得到所述卫星在J2000地心平赤道坐标系下的目标位置矢量和目标速度矢量;
卫星轨道六根数计算单元,用于根据所述目标位置矢量和所述目标速度矢量,计算得到所述卫星的卫星轨道六根数;
卫星理论位置计算单元,用于根据所述卫星轨道六根数和所述当前时刻,计算得到所述卫星在二体运动下的卫星理论位置;
摄动总和计算单元,用于根据所述卫星理论位置、所述卫星在地固系下的卫星位置和所述当前时刻,计算得到所述卫星在惯性坐标系下的摄动总和。
8.根据权利要求6所述的装置,其特征在于,所述非球形引力摄动确定模块包括:
地球引力计算单元,用于根据所述卫星位置,计算得到所述卫星在二体运动下的地球引力、及所述卫星对应的卫星经度、卫星纬度和地心距;
基准地球引力计算单元,用于根据所述卫星经度、所述卫星纬度、所述地心距和所述第三方引力场模型,计算得到所述卫星对应的基准地球引力;
非球形引力摄动计算单元,用于根据所述地球引力和所述基准地球引力,计算得到所述卫星对应的非球形引力摄动。
9.根据权利要求7所述的装置,其特征在于,所述日月引力摄动确定模块包括:
日月坐标计算单元,用于根据所述当前时刻,计算得到太阳在当前时刻的太阳坐标,及月球在当前时刻的月球坐标;
日月引力摄动计算单元,用于根据所述卫星在二体运动下的卫星理论位置、所述太阳坐标和所述月球坐标,计算得到所述卫星对应的日引力摄动和月引力摄动。
10.根据权利要求6所述的装置,其特征在于,所述卫星大气摄动确定模块包括:
摄动和值计算单元,用于计算得到所述日引力摄动、所述月引力摄动和所述太阳光压摄动之间的摄动和值;
卫星大气摄动获取单元,用于计算所述摄动总和与所述摄动和值之间的差值,并将该差值作为所述卫星大气摄动。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210494763.0A CN115061171A (zh) | 2022-05-07 | 2022-05-07 | 一种卫星大气摄动测算方法及装置 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210494763.0A CN115061171A (zh) | 2022-05-07 | 2022-05-07 | 一种卫星大气摄动测算方法及装置 |
Publications (1)
Publication Number | Publication Date |
---|---|
CN115061171A true CN115061171A (zh) | 2022-09-16 |
Family
ID=83197423
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202210494763.0A Pending CN115061171A (zh) | 2022-05-07 | 2022-05-07 | 一种卫星大气摄动测算方法及装置 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN115061171A (zh) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN116659543A (zh) * | 2023-06-21 | 2023-08-29 | 中国人民解放军61540部队 | 基于遥感卫星轨道根数的卫星位置姿态估计方法和装置 |
-
2022
- 2022-05-07 CN CN202210494763.0A patent/CN115061171A/zh active Pending
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN116659543A (zh) * | 2023-06-21 | 2023-08-29 | 中国人民解放军61540部队 | 基于遥感卫星轨道根数的卫星位置姿态估计方法和装置 |
CN116659543B (zh) * | 2023-06-21 | 2024-05-07 | 中国人民解放军61540部队 | 基于遥感卫星轨道根数的卫星位置姿态估计方法和装置 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Poghosyan et al. | CubeSat evolution: Analyzing CubeSat capabilities for conducting science missions | |
Williams et al. | OSIRIS-REx flight dynamics and navigation design | |
Van Patten et al. | A possible experiment with two counter-orbiting drag-free satellites to obtain a new test of einstein's general theory of relativity and improved measurements in geodesy | |
CN100533065C (zh) | 基于多天体路标的星际巡航自主导航方法 | |
CN111427002B (zh) | 地面测控天线指向卫星的方位角计算方法 | |
Bradley et al. | Cislunar navigation accuracy using optical observations of natural and artificial targets | |
Slavinskis et al. | Nanospacecraft fleet for multi-asteroid touring with electric solar wind sails | |
Jacobson | The orbits of the main Saturnian satellites, the Saturnian system gravity field, and the orientation of Saturn’s pole | |
Hugentobler et al. | Satellite orbits and attitude | |
Hesar et al. | Small body gravity field estimation using LIAISON supplemented optical navigation | |
CN115061171A (zh) | 一种卫星大气摄动测算方法及装置 | |
CN113589832B (zh) | 对地表固定区域目标稳定观测覆盖的星座快速设计方法 | |
Zeitlhöfler | Nominal and observation-based attitude realization for precise orbit determination of the Jason satellites | |
Langley | The orbits of GPS satellites | |
CN116125503A (zh) | 一种高精度卫星轨道确定及预报算法 | |
CN111141278B (zh) | 一种星下点定时回归的赤道轨道半长轴确定方法 | |
Wang et al. | Autonomous orbit determination using pulsars and inter-satellite ranging for Mars orbiters | |
Lee | CubeSat constellation implementation and management using differential drag | |
Straka et al. | Navigation of satellite measurements without ground control points | |
CN114894199B (zh) | 一种地月空间航天器的天基测定轨方法 | |
CN111366126A (zh) | 地面测站天线对卫星指向的视向量计算系统 | |
Karimi et al. | A Performance Based Comparison of Deep-space Navigation using Optical Communication and Conventional Navigation Techniques: Small Body Missions | |
Ibrahim | Attitude and orbit control of small satellites for autonomous terrestrial target tracking | |
CN112379398B (zh) | 一种地月空间卫星导航定位方法 | |
Class | SCHOOL OF ENGINEERING AND ARCHITECTURE |
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 |