CN110378012B - 一种考虑高阶重力场的严格回归轨道设计方法、系统及介质 - Google Patents

一种考虑高阶重力场的严格回归轨道设计方法、系统及介质 Download PDF

Info

Publication number
CN110378012B
CN110378012B CN201910640736.8A CN201910640736A CN110378012B CN 110378012 B CN110378012 B CN 110378012B CN 201910640736 A CN201910640736 A CN 201910640736A CN 110378012 B CN110378012 B CN 110378012B
Authority
CN
China
Prior art keywords
regression
orbit
satellite
semi
design
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.)
Active
Application number
CN201910640736.8A
Other languages
English (en)
Other versions
CN110378012A (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 Jiaotong University
Original Assignee
Shanghai Jiaotong University
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 Jiaotong University filed Critical Shanghai Jiaotong University
Priority to CN201910640736.8A priority Critical patent/CN110378012B/zh
Publication of CN110378012A publication Critical patent/CN110378012A/zh
Application granted granted Critical
Publication of CN110378012B publication Critical patent/CN110378012B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/06Power analysis or power optimisation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

本发明提供了一种考虑高阶重力场的严格回归轨道设计方法、系统及介质,包括:卫星轨道初设计步骤:依据卫星轨道设计总体输入,包括回归天数和回归圈数,进行考虑J2摄动影响的回归轨道参数解析设计,获得解析设计结果,通过牛顿迭代法对解析设计结果进行微分修正,获得J2项摄动卫星轨道初设计初值;严格回归轨道设计步骤:根据获得的严格回归轨道优化设计初值,将牛顿迭代法与多目标优化算法NSGA‑II相结合,通过高精度重力场下的回归轨道混合优化设计算法,实现严格回归轨道的快速优化设计。本发明相比现有方法考虑了更精确的动力学环境,并能够大大提升回归精度,满足合成孔径雷达干涉测量系统的应用需求。

Description

一种考虑高阶重力场的严格回归轨道设计方法、系统及介质
技术领域
本发明涉及航天技术领域,具体的涉及一种考虑地球高阶重力场的严格回归轨道设计方法。
背景技术
回归轨道是指卫星星下点轨迹经过一段时间后重复出现,且在同一纬度圈上相邻轨迹间距相同的一类特殊轨道。回归轨道可满足对特定地面目标进行周期性动态观测的任务要求,目前,已在对地成像观测和地球科学探测等任务中得到广泛应用。特别的,卫星轨道的严格回归特性保障了合成孔径雷达干涉测量(InSAR)这一类新型空间技术的应用和实施。InSAR卫星系统需要获取对同一地面目标两次观测所产生的相位差和空间距离差,通过解析其相互之间的关系以提取地表三维信息和形变信息,因此,针对高精度回归轨道设计方法研究具有重要意义。
目前针对回归轨道设计的研究大都只考虑了地球非球形引力J2项的摄动影响。通过简化地球非球形引力摄动影响,以解析求解回归轨道参数(主要是半长轴和偏心率)的设计思路,可满足对卫星轨道回归精度要求不高的航天任务需求。由于仅考虑J2摄动影响,轨道的设计回归精度与实际回归精度尚存在较大差异,并不能满足合成孔径雷达干涉测量系统的应用需求:一方面,合成孔径雷达干涉测量技术对卫星轨道的回归精度提出了更高要求;另一方面,考虑实际工程应用,为了在轨道设计阶段尽可能精确模拟卫星在轨运行的真实情况,回归轨道设计过程中应考虑高阶地球非球形摄动的影响。
专利文献CN104484493B(申请号:201410608513.0)公开了一种飞船返回预定落点回归轨道设计方法,通过选择轨道平近点角实现星下点准确经过预定落点,通过调整轨道半长轴使轨道满足回归特性,建立轨道半长轴和轨道平近点角2个参数双层迭代求解流程,得到轨道倾角、轨道半长轴和轨道平近点角相互匹配的设计参数,保证了飞船星下点轨迹每回归周期准确经过预定落点。
发明内容
针对现有技术中的缺陷,本发明的目的是提供一种考虑高阶重力场的严格回归轨道设计方法、系统及介质。
根据本发明提供的一种考虑高阶重力场的严格回归轨道设计方法,包括:
卫星轨道初设计步骤:依据卫星轨道设计总体输入,包括回归天数和回归圈数,进行考虑J2摄动影响的回归轨道参数解析设计,获得解析设计结果,通过牛顿迭代法对解析设计结果进行微分修正,获得J2项摄动卫星轨道初设计初值;
严格回归轨道设计步骤:根据获得的严格回归轨道优化设计初值,将牛顿迭代法与多目标优化算法NSGA-II相结合,通过高精度重力场下的回归轨道混合优化设计算法,实现严格回归轨道的快速优化设计。
优选地,所述卫星轨道初设计步骤:
轨道半长轴初步确定步骤:根据回归天数N和回归轨道数R,在二体运动模型下,初步确定轨道半长轴
Figure GDA0003090287480000021
J2摄动下的半长轴计算步骤:根据获得的初步确定的轨道半长轴
Figure GDA0003090287480000022
引入J2项摄动,得到J2摄动下的卫星半长轴
Figure GDA0003090287480000023
轨道倾角获取步骤:根据获得的J2摄动下的半长轴
Figure GDA0003090287480000024
再根据太阳同步约束,求得轨道倾角
Figure GDA0003090287480000025
升交点赤经获取步骤:根据分离点的经纬度
Figure GDA0003090287480000026
和分离时刻TS,可求得升交点赤经Ω。
优选地,所述轨道半长轴初步确定步骤:
在二体运动模型下基于回归天数N和回归轨道数R初步确定轨道半长轴
Figure GDA0003090287480000027
Figure GDA0003090287480000028
Figure GDA0003090287480000029
其中,
Figure GDA00030902874800000210
表示初步确定的轨道半长轴;
Figure GDA00030902874800000211
为地球引力常数;
P表示回归系数的倒数;
N表示回归天数;
R表示回归轨道数;
所述J2摄动下的半长轴计算步骤:
Figure GDA0003090287480000031
其中,
Figure GDA0003090287480000032
表示J2摄动下的卫星半长轴;
J2表示J2摄动项;
Figure GDA0003090287480000033
表示升交点赤经的进动值;
Figure GDA0003090287480000034
为地球赤道半径;
所述轨道倾角获取步骤:
Figure GDA0003090287480000035
Figure GDA0003090287480000036
其中,
J2表示J2摄动项;
year表示一年时间对应的秒数;
Figure GDA0003090287480000037
表示轨道倾角;
所述升交点赤经获取步骤:
假设S为分离点,为参考轨道上的某一特定时间TS,即分离时刻,其余的卫星轨道参数可由球面几何计算得到,给定纬度
Figure GDA0003090287480000038
和飞行方向,北面为正:
Figure GDA0003090287480000039
其中,
u表示纬度幅角;
Figure GDA00030902874800000310
表示纬度;
Figure GDA00030902874800000311
表示轨道倾角;
sign表示符号函数;
从S点到下一个升交点的飞行时间可表示为轨道周期的一部分,基于卫星自转,升交点赤经沿经度方向的偏移为:
Figure GDA0003090287480000041
其中,
ΔλN表示升交点赤经沿经度方向的偏移;
Tday表示一天时间对应的秒数,数值为86400;
基于球面几何,升交点赤经与S点的经度差可表示为:
Figure GDA0003090287480000042
其中,
Δλ表示升交点赤经与S点的经度差;
给定S点的经度值λ,轨道节点经度可由下式计算
λN=λ-Δλ+ΔλN
其中,
λ表示分离点S的经度值;
λN表示轨道节点经度;
最后,针对给定的时刻,升交点赤经可表示为
Ω=GMST(tUT)+λN
其中,
Ω表示升交点赤经;
GMST表示格林尼治时角;
tUT1表示平太阳时。
优选地,所述严格回归轨道设计步骤:
基于牛顿迭代法的回归轨道回归特性修正步骤:
根据获得的J2项摄动卫星轨道初设计初值半长轴参数a、轨道倾角i以及升交点赤经Ω,设定其余轨道参数,为优化初设计轨道的回归特性,提出基于牛顿迭代法的严格回归重力场修正方法:
设回归后的纬度
Figure GDA0003090287480000043
和回归后的经度λ为半长轴a和倾角i的函数:
Figure GDA0003090287480000051
λ=f(a,i)
对上式求导数,并且在整个回归周期内积分,从回归周期的开始时刻T0到回归周期的结束时刻Tf=T0+Tperiod,可得到下式
Figure GDA0003090287480000052
其中Δ代表在整个回归周期的差分运算符,Tperiod为回归周期。
Figure GDA0003090287480000053
为相对于半长轴a和倾角i的偏导数。
优选地,所述严格回归轨道设计步骤还包括:
基于微分修正的严格回归轨道修正步骤:
依据牛顿迭代算法计算原理,以卫星轨道根数为变量,构建基于PV=[px,Py,pz,vx,vy,vz]状态向量函数为:
Figure GDA0003090287480000054
其余,
px,py,pz,vx,vy,vz分别表示x方向位置、y方向位置、z方向位置、x方向速度、y方向速度、z方向速度;
a,e,i,w,Ω,M分别表示轨道半长轴、偏心率、轨道倾角、近地点幅角、升交点赤经、平近点角;
fpx表示x方向位置关于轨道根数的函数;
fpx(a,e,i,w,Ω,M)表示x方向位置关于轨道根数的函数;
fpy表示y方向位置关于轨道根数的函数;
fpz表示z方向位置关于轨道根数的函数;
fvx表示x方向速度关于轨道根数的函数;
fvy表示y方向速度关于轨道根数的函数;
fvz表示z方向速度关于轨道根数的函数;
依据上述状态变量函数,通过微分处理,可得高阶重力场环境下卫星轨道六要素牛顿迭代方程为
Figure GDA0003090287480000061
式中:
Δ、
Figure GDA0003090287480000062
为回归误差计算算子;
Δpx、Δpy、Δpz、Δvx、Δvy、Δvz为在当前轨道根数下,卫星状态分量的回归误差;
Figure GDA0003090287480000063
为依次给定的六个轨道根数微分量;
Figure GDA0003090287480000064
为在新轨道根数卫星状态分量的回归误差;
Δa,Δe,Δi,Δω,ΔΩ,ΔM分别表示轨道半长轴、偏心率、轨道倾角、近地点幅角、升交点赤经、平近点角等轨道6要素的微元;
基于NSGA-II的严格回归轨道多目标优化设计步骤:
将严格回归轨道设计问题描述为多目标优化问题,该问题可描述为:
优化目标:min(Δxp),min(xv)
优化变量:Δα=[Δa,Δe,Δi,Δω,ΔΩ,ΔM]
约束条件:
x′=f(t,x);
Δxp=xp(Tf)-xp(T0);
Δxv=xv(Tf)-xv(T0);
其中,
x′表示卫星状态变量的导数;
Δxp表示回归周期始末的位置差;
Δxv表示回归周期始末的速度差;
f为高阶重力场影响下的卫星轨道状态动力学方程;xp(T0),xp(Tf)分别为回归周期始末地固系下的位置矢量;xv(T0),xv(Tf)分别为回归周期始末地固系下的速度矢量;
为求解上述多目标优化问题,选择NSGA-II多目标优化算法进行求解,以求取最终的卫星轨道根数。
根据本发明提供的一种考虑高阶重力场的严格回归轨道设计系统,包括:
卫星轨道初设计模块:依据卫星轨道设计总体输入,包括回归天数和回归圈数,进行考虑J2摄动影响的回归轨道参数解析设计,获得解析设计结果,通过牛顿迭代法对解析设计结果进行微分修正,获得J2项摄动卫星轨道初设计初值;
严格回归轨道设计模块:根据获得的严格回归轨道优化设计初值,将牛顿迭代法与多目标优化算法NSGA-II相结合,通过高精度重力场下的回归轨道混合优化设计算法,实现严格回归轨道的快速优化设计。
优选地,所述卫星轨道初设计模块:
轨道半长轴初步确定模块:根据回归天数N和回归轨道数R,在二体运动模型下,初步确定轨道半长轴
Figure GDA0003090287480000071
J2摄动下的半长轴计算模块:根据获得的初步确定的轨道半长轴
Figure GDA0003090287480000072
引入J2项摄动,得到J2摄动下的卫星半长轴
Figure GDA0003090287480000073
轨道倾角获取模块:根据获得的J2摄动下的半长轴
Figure GDA0003090287480000074
再根据太阳同步约束,求得轨道倾角
Figure GDA0003090287480000075
升交点赤经获取模块:根据分离点的经纬度
Figure GDA0003090287480000076
和分离时刻TS,可求得升交点赤经Ω。
优选地,所述轨道半长轴初步确定模块:
在二体运动模型下基于回归天数N和回归轨道数R初步确定轨道半长轴
Figure GDA0003090287480000077
Figure GDA0003090287480000078
Figure GDA0003090287480000081
其中,
Figure GDA0003090287480000082
表示初步确定的轨道半长轴;
Figure GDA0003090287480000083
为地球引力常数;
P表示回归系数的倒数;
N表示回归天数;
R表示回归轨道数;
所述J2摄动下的半长轴计算模块:
Figure GDA0003090287480000084
其中,
Figure GDA0003090287480000085
表示J2摄动下的卫星半长轴;
J2表示J2摄动项;
Figure GDA0003090287480000086
表示升交点赤经的进动值;
Figure GDA0003090287480000087
为地球赤道半径;
所述轨道倾角获取模块:
Figure GDA0003090287480000088
Figure GDA0003090287480000089
其中,
J2表示J2摄动项;
year表示一年时间对应的秒数;
Figure GDA00030902874800000810
表示轨道倾角;
所述升交点赤经获取模块:
假设S为分离点,为参考轨道上的某一特定时间TS,即分离时刻,其余的卫星轨道参数可由球面几何计算得到,给定纬度
Figure GDA00030902874800000811
和飞行方向,北面为正:
Figure GDA0003090287480000091
其中,
u表示纬度幅角;
Figure GDA0003090287480000092
表示纬度;
Figure GDA0003090287480000093
表示轨道倾角;
sign表示符号函数;
从S点到下一个升交点的飞行时间可表示为轨道周期的一部分,基于卫星自转,升交点赤经沿经度方向的偏移为:
Figure GDA0003090287480000094
其中,
ΔλN表示升交点赤经沿经度方向的偏移;
Tday表示一天时间对应的秒数,数值为86400;
基于球面几何,升交点赤经与S点的经度差可表示为:
Figure GDA0003090287480000095
其中,
Δλ表示升交点赤经与S点的经度差;
给定S点的经度值λ,轨道节点经度可由下式计算
λN=λ-Δλ+ΔλN
其中,
λ表示分离点S的经度值;
λN表示轨道节点经度;
最后,针对给定的时刻,升交点赤经可表示为
Ω=GMST(tUT)+λN
其中,
Ω表示升交点赤经;
GMST表示格林尼治时角;
tUT1表示平太阳时。
优选地,所述严格回归轨道设计模块:
基于牛顿迭代法的回归轨道回归特性修正模块:
根据获得的J2项摄动卫星轨道初设计初值半长轴参数a、轨道倾角i以及升交点赤经Ω,设定其余轨道参数,为优化初设计轨道的回归特性,提出基于牛顿迭代法的严格回归重力场修正方法:
设回归后的纬度
Figure GDA0003090287480000105
和回归后的经度λ为半长轴a和倾角i的函数:
Figure GDA0003090287480000101
λ=f(a,i)
对上式求导数,并且在整个回归周期内积分,从回归周期的开始时刻T0到回归周期的结束时刻Tf=T0+Tperiod,可得到下式
Figure GDA0003090287480000102
其中Δ代表在整个回归周期的差分运算符,Tperiod为回归周期。
Figure GDA0003090287480000103
为相对于半长轴a和倾角i的偏导数;
基于微分修正的严格回归轨道修正模块:
依据牛顿迭代算法计算原理,以卫星轨道根数为变量,构建基于PV=[px,py,pz,vx,vy,vz]状态向量函数为:
Figure GDA0003090287480000104
其余,
px,py,pz,vx,vy,vz分别表示x方向位置、y方向位置、z方向位置、x方向速度、y方向速度、z方向速度;
a,e,i,w,Ω,M分别表示轨道半长轴、偏心率、轨道倾角、近地点幅角、升交点赤经、平近点角;
fpx表示x方向位置关于轨道根数的函数;
fpx(a,e,i,w,Ω,M)表示x方向位置关于轨道根数的函数;
fpy表示y方向位置关于轨道根数的函数;
fpz表示z方向位置关于轨道根数的函数;
fvx表示x方向速度关于轨道根数的函数;
fvy表示y方向速度关于轨道根数的函数;
fvz表示z方向速度关于轨道根数的函数;
依据上述状态变量函数,通过微分处理,可得高阶重力场环境下卫星轨道六要素牛顿迭代方程为
Figure GDA0003090287480000111
式中:
Δ、
Figure GDA0003090287480000112
为回归误差计算算子;
Δpx、Δpy、Δpz、Δvx、Δvy、Δvz为在当前轨道根数下,卫星状态分量的回归误差;
Figure GDA0003090287480000113
为依次给定的六个轨道根数微分量;
Figure GDA0003090287480000114
为在新轨道根数卫星状态分量的回归误差;
Δa,Δe,Δi,Δω,ΔΩ,ΔM分别表示轨道半长轴、偏心率、轨道倾角、近地点幅角、升交点赤经、平近点角等轨道6要素的微元;
基于NSGA-II的严格回归轨道多目标优化设计模块:
将严格回归轨道设计问题描述为多目标优化问题,该问题可描述为:
优化目标:min(Δxp),min(xv)
优化变量:Δα=[Δa,Δe,Δi,Δω,ΔΩ,ΔM]
约束条件:
x′=f(t,x);
Δxp=xp(Tf)-xp(T0);
Δxv=xv(Tf)-xv(T0);
其中,
x′表示卫星状态变量的导数;
Δxp表示回归周期始末的位置差;
Δxv表示回归周期始末的速度差;
f为高阶重力场影响下的卫星轨道状态动力学方程;xp(T0),xp(Tf)分别为回归周期始末地固系下的位置矢量;xv(T0),xv(Tf)分别为回归周期始末地固系下的速度矢量;
为求解上述多目标优化问题,选择NSGA-II多目标优化算法进行求解,以求取最终的卫星轨道根数。
根据本发明提供的一种存储有计算机程序的计算机可读存储介质,其特征在于,所述计算机程序被处理器执行时实现上述中任一项所述的考虑高阶重力场的严格回归轨道设计方法的步骤。
与现有技术相比,本发明具有如下的有益效果:
本发明相比现有方法考虑了更精确的动力学环境,并能够大大提升回归精度,满足合成孔径雷达干涉测量系统的应用需求。
具体实施方式
下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明,但不以任何形式限制本发明。应当指出的是,对本领域的普通技术人员来说,在不脱离本发明构思的前提下,还可以做出若干变化和改进。这些都属于本发明的保护范围。
根据本发明提供的一种考虑高阶重力场的严格回归轨道设计方法,包括:
卫星轨道初设计步骤:依据卫星轨道设计总体输入,包括回归天数和回归圈数,进行考虑J2摄动影响的回归轨道参数解析设计,获得解析设计结果,通过牛顿迭代法对解析设计结果进行微分修正,获得J2项摄动卫星轨道初设计初值;
严格回归轨道设计步骤:根据获得的严格回归轨道优化设计初值,将牛顿迭代法与多目标优化算法NSGA-II相结合,通过高精度重力场下的回归轨道混合优化设计算法,实现严格回归轨道的快速优化设计。
具体地,所述卫星轨道初设计步骤:
轨道半长轴初步确定步骤:根据回归天数N和回归轨道数R,在二体运动模型下,初步确定轨道半长轴
Figure GDA0003090287480000131
J2摄动下的半长轴计算步骤:根据获得的初步确定的轨道半长轴
Figure GDA0003090287480000132
引入J2项摄动,得到J2摄动下的卫星半长轴
Figure GDA0003090287480000133
轨道倾角获取步骤:根据获得的J2摄动下的半长轴
Figure GDA0003090287480000134
再根据太阳同步约束,求得轨道倾角
Figure GDA0003090287480000135
升交点赤经获取步骤:根据分离点的经纬度
Figure GDA00030902874800001312
和分离时刻TS,可求得升交点赤经Ω。
具体地,所述轨道半长轴初步确定步骤:
在二体运动模型下基于回归天数N和回归轨道数R初步确定轨道半长轴
Figure GDA0003090287480000136
Figure GDA0003090287480000137
Figure GDA0003090287480000138
其中,
Figure GDA0003090287480000139
表示初步确定的轨道半长轴;
Figure GDA00030902874800001310
为地球引力常数;
P表示回归系数的倒数;
N表示回归天数;
R表示回归轨道数;
所述J2摄动下的半长轴计算步骤:
Figure GDA00030902874800001311
其中,
Figure GDA0003090287480000141
表示J2摄动下的卫星半长轴;
J2表示J2摄动项;
Figure GDA0003090287480000142
表示升交点赤经的进动值;
Figure GDA0003090287480000143
为地球赤道半径;
所述轨道倾角获取步骤:
Figure GDA0003090287480000144
Figure GDA0003090287480000145
其中,
J2表示J2摄动项;
year表示一年时间对应的秒数;
Figure GDA0003090287480000146
表示轨道倾角;
所述升交点赤经获取步骤:
假设S为分离点,为参考轨道上的某一特定时间TS,即分离时刻,其余的卫星轨道参数可由球面几何计算得到,给定纬度
Figure GDA00030902874800001411
和飞行方向,北面为正:
Figure GDA0003090287480000147
其中,
u表示纬度幅角;
Figure GDA0003090287480000148
表示纬度;
Figure GDA0003090287480000149
表示轨道倾角;
sign表示符号函数;
从S点到下一个升交点的飞行时间可表示为轨道周期的一部分,基于卫星自转,升交点赤经沿经度方向的偏移为:
Figure GDA00030902874800001410
其中,
ΔλN表示升交点赤经沿经度方向的偏移;
Tday表示一天时间对应的秒数,数值为86400;
基于球面几何,升交点赤经与S点的经度差可表示为:
Figure GDA0003090287480000151
其中,
Δλ表示升交点赤经与S点的经度差;
给定S点的经度值λ,轨道节点经度可由下式计算
λN=λ-Δλ+ΔλN
其中,
λ表示分离点S的经度值;
λN表示轨道节点经度;
最后,针对给定的时刻,升交点赤经可表示为
Ω=GMST(tUT1)+λN
其中,
Ω表示升交点赤经;
GMST表示格林尼治时角;
tUT1表示平太阳时。
具体地,所述严格回归轨道设计步骤:
基于牛顿迭代法的回归轨道回归特性修正步骤:
根据获得的J2项摄动卫星轨道初设计初值半长轴参数a、轨道倾角i以及升交点赤经Ω,设定其余轨道参数,为优化初设计轨道的回归特性,提出基于牛顿迭代法的严格回归重力场修正方法:
设回归后的纬度
Figure GDA0003090287480000153
和回归后的经度λ为半长轴a和倾角i的函数:
Figure GDA0003090287480000152
λ=f(a,i)
对上式求导数,并且在整个回归周期内积分,从回归周期的开始时刻T0到回归周期的结束时刻Tf=T0+Tperiod,可得到下式
Figure GDA0003090287480000161
其中Δ代表在整个回归周期的差分运算符,Tperiod为回归周期。
Figure GDA0003090287480000162
为相对于半长轴a和倾角i的偏导数。
具体地,所述严格回归轨道设计步骤还包括:
基于微分修正的严格回归轨道修正步骤:
依据牛顿迭代算法计算原理,以卫星轨道根数为变量,构建基于PV=[px,py,pz,vx,vy,vz]状态向量函数为:
Figure GDA0003090287480000163
其余,
px,py,pz,vx,vy,vz分别表示x方向位置、y方向位置、z方向位置、x方向速度、y方向速度、z方向速度;
a,e,i,w,Ω,M分别表示轨道半长轴、偏心率、轨道倾角、近地点幅角、升交点赤经、平近点角;
fpx表示x方向位置关于轨道根数的函数;
fpx(a,e,i,w,Ω,M)表示x方向位置关于轨道根数的函数;
fpy表示y方向位置关于轨道根数的函数;
fpz表示z方向位置关于轨道根数的函数;
fvx表示x方向速度关于轨道根数的函数;
fvy表示y方向速度关于轨道根数的函数;
fvz表示z方向速度关于轨道根数的函数;
依据上述状态变量函数,通过微分处理,可得高阶重力场环境下卫星轨道六要素牛顿迭代方程为
Figure GDA0003090287480000171
式中:
Δ、
Figure GDA0003090287480000172
为回归误差计算算子;
Δpx、Δpy、Δpz、Δvx、Δvy、Δvz为在当前轨道根数下,卫星状态分量的回归误差;
Figure GDA0003090287480000173
为依次给定的六个轨道根数微分量;
Figure GDA0003090287480000174
为在新轨道根数卫星状态分量的回归误差;
Δa,Δe,Δi,Δω,ΔΩ,ΔM分别表示轨道半长轴、偏心率、轨道倾角、近地点幅角、升交点赤经、平近点角等轨道6要素的微元;
基于NSGA-II的严格回归轨道多目标优化设计步骤:
将严格回归轨道设计问题描述为多目标优化问题,该问题可描述为:
优化目标:min(Δxp),min(xv)
优化变量:Δα=[Δa,Δe,Δi,Δω,ΔΩ,ΔM]
约束条件:
x′=f(t,x);
Δxp=xp(Tf)-xp(T0);
Δxv=xv(Tf)-xv(T0);
其中,
x′表示卫星状态变量的导数;
Δxp表示回归周期始末的位置差;
Δxv表示回归周期始末的速度差;
f为高阶重力场影响下的卫星轨道状态动力学方程;xp(T0),xp(Tf)分别为回归周期始末地固系下的位置矢量;xv(T0),xv(Tf)分别为回归周期始末地固系下的速度矢量;
为求解上述多目标优化问题,选择NSGA-II多目标优化算法进行求解,以求取最终的卫星轨道根数。
本发明提供的考虑高阶重力场的严格回归轨道设计系统,可以通过本发明给的考虑高阶重力场的严格回归轨道设计方法的步骤流程实现。本领域技术人员可以将所述考虑高阶重力场的严格回归轨道设计方法,理解为所述考虑高阶重力场的严格回归轨道设计系统的一个优选例。
根据本发明提供的一种考虑高阶重力场的严格回归轨道设计系统,包括:
卫星轨道初设计模块:依据卫星轨道设计总体输入,包括回归天数和回归圈数,进行考虑J2摄动影响的回归轨道参数解析设计,获得解析设计结果,通过牛顿迭代法对解析设计结果进行微分修正,获得J2项摄动卫星轨道初设计初值;
严格回归轨道设计模块:根据获得的严格回归轨道优化设计初值,将牛顿迭代法与多目标优化算法NSGA-II相结合,通过高精度重力场下的回归轨道混合优化设计算法,实现严格回归轨道的快速优化设计。
具体地,所述卫星轨道初设计模块:
轨道半长轴初步确定模块:根据回归天数N和回归轨道数R,在二体运动模型下,初步确定轨道半长轴
Figure GDA0003090287480000181
J2摄动下的半长轴计算模块:根据获得的初步确定的轨道半长轴
Figure GDA0003090287480000182
引入J2项摄动,得到J2摄动下的卫星半长轴
Figure GDA0003090287480000183
轨道倾角获取模块:根据获得的J2摄动下的半长轴
Figure GDA0003090287480000184
再根据太阳同步约束,求得轨道倾角
Figure GDA0003090287480000185
升交点赤经获取模块:根据分离点的经纬度
Figure GDA0003090287480000187
和分离时刻TS,可求得升交点赤经Ω。
具体地,所述轨道半长轴初步确定模块:
在二体运动模型下基于回归天数N和回归轨道数R初步确定轨道半长轴
Figure GDA0003090287480000186
Figure GDA0003090287480000191
Figure GDA0003090287480000192
其中,
Figure GDA0003090287480000193
表示初步确定的轨道半长轴;
Figure GDA0003090287480000194
为地球引力常数;
P表示回归系数的倒数;
N表示回归天数;
R表示回归轨道数;
所述J2摄动下的半长轴计算模块:
Figure GDA0003090287480000195
其中,
Figure GDA0003090287480000196
表示J2摄动下的卫星半长轴;
J2表示J2摄动项;
Figure GDA0003090287480000197
表示升交点赤经的进动值;
Figure GDA0003090287480000198
为地球赤道半径;
所述轨道倾角获取模块:
Figure GDA0003090287480000199
Figure GDA00030902874800001910
其中,
J2表示J2摄动项;
year表示一年时间对应的秒数;
Figure GDA00030902874800001911
表示轨道倾角;
所述升交点赤经获取模块:
假设S为分离点,为参考轨道上的某一特定时间TS,即分离时刻,其余的卫星轨道参数可由球面几何计算得到,给定纬度
Figure GDA0003090287480000206
和飞行方向,北面为正:
Figure GDA0003090287480000201
其中,
u表示纬度幅角;
Figure GDA0003090287480000202
表示纬度;
Figure GDA0003090287480000203
表示轨道倾角;
sign表示符号函数;
从S点到下一个升交点的飞行时间可表示为轨道周期的一部分,基于卫星自转,升交点赤经沿经度方向的偏移为:
Figure GDA0003090287480000204
其中,
ΔλN表示升交点赤经沿经度方向的偏移;
Tday表示一天时间对应的秒数,数值为86400;
基于球面几何,升交点赤经与S点的经度差可表示为:
Figure GDA0003090287480000205
其中,
Δλ表示升交点赤经与S点的经度差;
给定S点的经度值λ,轨道节点经度可由下式计算
λN=λ-Δλ+ΔλN
其中,
λ表示分离点S的经度值;
λN表示轨道节点经度;
最后,针对给定的时刻,升交点赤经可表示为
Ω=GMST(tUT1)+λN
其中,
Ω表示升交点赤经;
GMST表示格林尼治时角;
tUT1表示平太阳时。
具体地,所述严格回归轨道设计模块:
基于牛顿迭代法的回归轨道回归特性修正模块:
根据获得的J2项摄动卫星轨道初设计初值半长轴参数a、轨道倾角i以及升交点赤经Ω,设定其余轨道参数,为优化初设计轨道的回归特性,提出基于牛顿迭代法的严格回归重力场修正方法:
设回归后的纬度
Figure GDA0003090287480000215
和回归后的经度λ为半长轴a和倾角i的函数:
Figure GDA0003090287480000211
λ=f(a,i)
对上式求导数,并且在整个回归周期内积分,从回归周期的开始时刻T0到回归周期的结束时刻Tf=T0+Tperiod,可得到下式
Figure GDA0003090287480000212
其中Δ代表在整个回归周期的差分运算符,Tperiod为回归周期。
Figure GDA0003090287480000213
为相对于半长轴a和倾角i的偏导数;
基于微分修正的严格回归轨道修正模块:
依据牛顿迭代算法计算原理,以卫星轨道根数为变量,构建基于PV=[px,py,pz,vx,vy,vz]状态向量函数为:
Figure GDA0003090287480000214
其余,
px,py,pz,vx,vy,vz分别表示x方向位置、y方向位置、z方向位置、x方向速度、y方向速度、z方向速度;
a,e,i,w,Ω,M分别表示轨道半长轴、偏心率、轨道倾角、近地点幅角、升交点赤经、平近点角;
fpx表示x方向位置关于轨道根数的函数;
fpx(a,e,i,w,Ω,M)表示x方向位置关于轨道根数的函数;
fpy表示y方向位置关于轨道根数的函数;
fpz表示z方向位置关于轨道根数的函数;
fvx表示x方向速度关于轨道根数的函数;
fvy表示y方向速度关于轨道根数的函数;
fvz表示z方向速度关于轨道根数的函数;
依据上述状态变量函数,通过微分处理,可得高阶重力场环境下卫星轨道六要素牛顿迭代方程为
Figure GDA0003090287480000221
式中:
Δ、
Figure GDA0003090287480000222
为回归误差计算算子;
Δpx、Δpy、Δpz、Δvx、Δvy、Δvz为在当前轨道根数下,卫星状态分量的回归误差;
Figure GDA0003090287480000223
为依次给定的六个轨道根数微分量;
Figure GDA0003090287480000224
为在新轨道根数卫星状态分量的回归误差;
Δa,Δe,Δi,Δω,ΔΩ,ΔM分别表示轨道半长轴、偏心率、轨道倾角、近地点幅角、升交点赤经、平近点角等轨道6要素的微元;
基于NSGA-II的严格回归轨道多目标优化设计模块:
将严格回归轨道设计问题描述为多目标优化问题,该问题可描述为:
优化目标:min(Δxp),min(xv)
优化变量:Δα=[Δa,Δe,Δi,Δω,ΔΩ,ΔM]
约束条件:
x′=f(t,x);
Δxp=xp(Tf)-xp(T0);
Δxv=xv(Tf)-xv(T0);
其中,
x′表示卫星状态变量的导数;
Δxp表示回归周期始末的位置差;
Δxv表示回归周期始末的速度差;
f为高阶重力场影响下的卫星轨道状态动力学方程;xp(T0),xp(Tf)分别为回归周期始末地固系下的位置矢量;xv(T0),xv(Tf)分别为回归周期始末地固系下的速度矢量;
为求解上述多目标优化问题,选择NSGA-II多目标优化算法进行求解,以求取最终的卫星轨道根数。
根据本发明提供的一种存储有计算机程序的计算机可读存储介质,其特征在于,所述计算机程序被处理器执行时实现上述中任一项所述的考虑高阶重力场的严格回归轨道设计方法的步骤。
下面通过优选例,对本发明进行更为具体地说明。
优选例1:
首先,依据卫星轨道设计总体输入,包括回归天数和回归圈数,给出了考虑J2摄动影响的回归轨道参数解析设计方法,并通过牛顿迭代法对解析设计结果进行微分修正以作为严格回归轨道优化设计初值;接着,通过将牛顿迭代法与多目标优化算法NSGA-II相结合,研究了一种高精度重力场下的回归轨道混合优化设计算法,以实现严格回归轨道的快速优化设计。主要内容为:
1)考虑J2项影响下的卫星轨道初设计
在仅考虑J2重力场下,回归轨道解析求解思路可描述为:首先,由回归天数N和回归轨道数R,在两体运动假设条件下,初步确定轨道半长轴
Figure GDA0003090287480000231
其次,引入J2项摄动,得到J2项摄动下的半长轴
Figure GDA0003090287480000232
再次,由太阳同步约束,求得轨道倾角
Figure GDA0003090287480000233
最后依据分离点的经纬度
Figure GDA00030902874800002413
和分离时刻TS,可求得升交点赤经Ω。上述回归轨道解析求解思路的具体流程为:
考虑两体运动模型,卫星半长轴可表示为
Figure GDA0003090287480000241
Figure GDA0003090287480000242
其中,
Figure GDA0003090287480000243
表示卫星半长轴;
Figure GDA0003090287480000244
为地球引力常数;
P表示单个卫星轨道周期;
N表示回归天数;
R表示回归轨道圈数;
考虑J2摄动,忽略偏心率影响,可得到
Figure GDA0003090287480000245
其中,
Figure GDA0003090287480000246
表示J2摄动下的卫星半长轴;
J2表示J2摄动项;
Figure GDA0003090287480000247
表示升交点赤经的进动值;
Figure GDA0003090287480000248
为地球赤道半径且
Figure GDA0003090287480000249
其中,
J2表示J2摄动项;
Figure GDA00030902874800002410
表示轨道倾角;
Figure GDA00030902874800002411
为升交点赤经的进动值,由太阳同步特性决定。由式(1)-(3)可求得,
Figure GDA00030902874800002412
Figure GDA0003090287480000251
假设S为分离点,为参考轨道上的某一特定时间TS(分离时刻),其余的卫星轨道参数可由球面几何计算得到。给定纬度
Figure GDA0003090287480000256
和飞行方向,北面为正:
Figure GDA0003090287480000252
为给定点的纬度幅角值。
其中,
Figure GDA0003090287480000253
表示纬度;
i表示轨道倾角;
sign表示符号函数;
从S点到下一个升交点的飞行时间可表示为轨道周期的一部分。基于卫星自转,升交点赤经沿经度方向的偏移为:
Figure GDA0003090287480000254
其中,
ΔλN表示升交点赤经沿经度方向的偏移;
u表示纬度幅角;
Tday表示一天时间对应的秒数,数值为86400;
基于球面几何,升交点赤经与S点的经度差可表示为
Figure GDA0003090287480000255
其中,
Δλ表示升交点赤经与S点的经度差;
给定S点的经度值λ,轨道节点经度可由下式计算
λN=λ-Δλ+ΔλN (8)
其中,
λ表示分离点S的经度值;
λN表示轨道节点经度;
最后,针对给定的时刻,升交点赤经可表示为
Ω=GMST(tUT1)+λN (9)
其中,
Ω表示升交点赤经;
GMST为格林尼治时角,tUT1为平太阳时。
2)基于牛顿迭代法的回归轨道回归特性修正方法
J2项摄动下的卫星轨道初设计结果为严格回归轨道设计奠定了基础,但因仅考虑了地球非球谐摄动中的第2项,在高阶重力场环境下的回归误差高达几十千米。为缩小高阶重力场环境下精密回归轨道优化求解空间,提升优化求解速度,以J2项摄动卫星轨道初设计输出的卫星轨道半长轴参数a、轨道倾角i和升交点赤经Ω为输入,在设定其余轨道参数的条件下,为优化初设计轨道的回归特性,提出了基于牛顿迭代法的严格回归重力场修正方法。
该方法的研究思路为,假设回归后的纬度
Figure GDA0003090287480000265
和经度λ为半长轴a和倾角i的函数:
Figure GDA0003090287480000261
对上式求导数,并且在整个回归周期内积分,从回归周期的开始时刻T0到回归周期的结束时刻Tf=T0+Tperiod,可得到下式
Figure GDA0003090287480000262
其中Δ代表在整个回归周期的差分运算符,Tperiod为回归周期。
Figure GDA0003090287480000263
Figure GDA0003090287480000264
为相对于半长轴a和倾角i的偏导数。
3)基于微分修正的严格回归轨道修正方法
高阶重力场环境下,受高阶重力场模型强非线性特性,以及轨道回归递推周期过长的因素影响。利用单一的牛顿迭代方法已无法完成高精度回归轨道优化设计。考虑到低价重力场回归轨道优化结果在高阶重力场环境起始和末状态依然存在较大偏差,如直接选用现有的多目标优化算法进行寻优求解,寻优空间大,迭代次数多,优化时间长。通过均衡考虑回归轨道寻优精度和速度,结合已有的低阶重力场回归轨道修正处理经验,提出基于牛顿迭代法和多目标优化求解的高阶重力场回归轨道快速、精确设计方法。首先通过建立高阶重力场环境下的牛顿迭代方程,对低阶设计结果进行修正,大幅减小精密回归轨道优化设计寻优搜索空间;其次,引入计算量小、收敛速度较快的NSGA-II多目标优化算法,完成精密回归轨道优化设计。
依据牛顿迭代算法计算原理,以卫星轨道根数为变量,构建基于PV=[px,py,pz,vx,vy,vz]状态向量函数为
Figure GDA0003090287480000271
其余,
px,py,pz,vx,vy,vz分别表示x方向位置、y方向位置、z方向位置、x方向速度、y方向速度、z方向速度;
a,e,i,w,Ω,M分别表示轨道半长轴、偏心率、轨道倾角、近地点幅角、升交点赤经、平近点角;
fpx表示x方向位置关于轨道根数的函数;
fpx(a,e,i,w,Ω,M)表示x方向位置关于轨道根数的函数;
fpy表示y方向位置关于轨道根数的函数;
fpz表示z方向位置关于轨道根数的函数;
fvx表示x方向速度关于轨道根数的函数;
fvy表示y方向速度关于轨道根数的函数;
fvz表示z方向速度关于轨道根数的函数;
依据上述状态变量函数,通过微分处理,可得高阶重力场环境下卫星轨道六要素牛顿迭代方程为
Figure GDA0003090287480000272
式中:
Δ、
Figure GDA0003090287480000281
为回归误差计算算子;
Δpx、Δpy、Δpz、Δvx、Δvy、Δvz为在当前轨道根数下,卫星状态分量的回归误差;
Figure GDA0003090287480000282
为依次给定的六个轨道根数微分量;
Figure GDA0003090287480000283
为在新轨道根数卫星状态分量的回归误差;
Δa,Δe,Δi,Δω,ΔΩ,ΔM分别表示轨道半长轴、偏心率、轨道倾角、近地点幅角、升交点赤经、平近点角等轨道6要素的微元;
4)基于NSGA-II的严格回归轨道多目标优化设计方法
在高阶重力场作用下,基于上述的牛顿迭代法已无法进一步优化卫星轨道的回归特性,故为进一步提高所设计轨道的严格回归特性,需将严格回归轨道设计问题描述为多目标优化问题。该问题可描述为:
优化目标:min(Δxp),min(xv)
优化变量:Δα=[Δa,Δe,Δi,Δω,ΔΩ,ΔM]
约束条件:
x′=f(t,x);
Δxp=xp(Tf)-xp(T0);
Δxv=xv(Tf)-xv(T0);
其中,
x′表示卫星状态变量的导数;
Δxp表示回归周期始末的位置差;
Δxv表示回归周期始末的速度差;
f为高阶重力场影响下的卫星轨道状态动力学方程;xp(T0),xp(Tf)分别为回归周期始末地固系下的位置矢量;xv(T0),xv(Tf)分别为回归周期始末地固系下的速度矢量。
为求解上述多目标优化问题,这里选择NSGA-II多目标优化算法进行求解,以求取最终的卫星轨道根数。
优选例2:
步骤1:在二体运动模型下基于回归天数N和回归轨道数R初步确定轨道半长轴
Figure GDA0003090287480000291
Figure GDA0003090287480000292
其中
Figure GDA0003090287480000293
为地球引力常数,
Figure GDA0003090287480000294
P表示回归系数的倒数,N表示回归天数,R表示回归轨道数;
步骤2:求取J2项摄动下的半长轴
Figure GDA0003090287480000295
和轨道倾角
Figure GDA0003090287480000296
Figure GDA0003090287480000297
其中
Figure GDA0003090287480000298
为地球赤道半径且
Figure GDA0003090287480000299
为升交点赤经的进动值,由太阳同步特性决定。
J2项摄动下的轨道倾角
Figure GDA00030902874800002910
Figure GDA00030902874800002911
步骤3:求取J2项摄动下的升交点赤经
Figure GDA00030902874800002912
假设S为分离点,为参考轨道上的某一特定时间TS(分离时刻),其余的卫星轨道参数可由球面几何计算得到。给定纬度
Figure GDA00030902874800002913
和飞行方向,北面为正:
Figure GDA00030902874800002914
为给定点的纬度幅角值。从S点到下一个升交点的飞行时间可表示为轨道周期的一部分。基于卫星自转,升交点赤经沿经度方向的偏移为:
Figure GDA00030902874800002915
基于球面几何,升交点赤经与S点的经度差可表示为
Figure GDA0003090287480000301
给定S点的经度值λ,轨道节点经度可由下式计算
λNN=λ-Δλ+ΔλN
最后,针对给定的时刻,升交点赤经可表示为
Ωj2=GMST(tUT1)+λN
其中GMST为格林尼治时角,tUT1为平太阳时。
步骤4:确定回归轨道六要素在J2项摄动下的预估值:
上述步骤确定了半长轴、轨道倾角和升交点赤经,此外,对于近圆轨道,偏心率可近似看为0,即
Figure GDA0003090287480000302
考虑轨道的冻结特性,近地点幅角
Figure GDA0003090287480000303
可取为90°。给定卫星在轨道的初始位置(由平近点角
Figure GDA0003090287480000304
表示),由以上推导可得回归轨道六要素预估设计值
Figure GDA0003090287480000305
其中
Figure GDA0003090287480000306
步骤5:基于J2项摄动的回归轨道参数微分迭代求解:
以步骤4确定的J2项摄动轨道预估值为回归轨道参数迭代求解初值,假定卫星的地心经度λ和地心纬度
Figure GDA00030902874800003012
与轨道半长轴am和轨道倾角im存在非线性函数关系:
Figure GDA0003090287480000307
依据公式,将卫星在A、B两点对应的地心地固坐标系下位置矢量
Figure GDA0003090287480000308
分别转为对应的球坐标系下位置量
Figure GDA0003090287480000309
依据非线性方程线性近似原理,构造卫星在A、B两点的经纬度差分关于半长轴修正量Δam和轨道倾角修正量Δim的线性方程,从而得到每一步迭代的修正量:
Figure GDA00030902874800003010
其中,
Figure GDA00030902874800003011
基于J2摄动项的微分迭代法计算过程如下所示:
1)初始化变量k=0;
2)输入轨道预报的开始时刻ts和结束时刻te,轨道参数预估值
Figure GDA0003090287480000311
半长轴微元δam、轨道倾角微元δim
3)令初始时刻卫星轨道参数为
Figure GDA0003090287480000312
考虑J2摄动影响,利用解析法预报轨道,得到初始时刻和回归时刻星下点在地固坐标系下的地心经纬度
Figure GDA0003090287480000313
求得经纬度差分量Δλ、
Figure GDA0003090287480000314
4)依据多元函数展开原理,用数值方法近似求解[f;g]关于变量am和im的偏导数,构造状态量对设计变量的雅克比矩阵
Figure GDA0003090287480000315
(以
Figure GDA0003090287480000316
的求解为例):
2-1)保持轨道倾角及其他轨道要素不变,令
Figure GDA0003090287480000317
2-2)输入轨道参数
Figure GDA0003090287480000318
利用解析法预报轨道,得到始末时刻星下点的经纬度差分量δλ和
Figure GDA0003090287480000319
2-3)求得
Figure GDA00030902874800003110
5)由雅克比矩阵和经纬度差分量Δλ、
Figure GDA00030902874800003111
求得设计变量的迭代修正量Δam、Δim
Figure GDA00030902874800003112
6)利用Δam、Δim修正
Figure GDA00030902874800003113
Figure GDA00030902874800003114
得到新的半长轴
Figure GDA00030902874800003115
和轨道倾角
Figure GDA00030902874800003116
更新初始时刻卫星轨道参数
Figure GDA00030902874800003117
7)判断迭代算法是否收敛,即经纬度差分量Δλ、
Figure GDA00030902874800003118
是否足够小:若成立,则令
Figure GDA00030902874800003119
并终止计算;若不成立,令m=m+1并返回至步骤3)。
通过微分迭代算法对轨道参数预估值
Figure GDA00030902874800003120
进行修正,最终可得到J2摄动影响下的严格回归轨道参数
Figure GDA00030902874800003121
步骤6:通过平均轨道要素和密切轨道要素转换关系将
Figure GDA00030902874800003122
转换为对应的密切轨道要素
Figure GDA00030902874800003123
步骤7:基于高阶地球重力场的回归轨道参数微分迭代求解
建立高阶重力场环境下卫星轨道六要素牛顿迭代方程为
Figure GDA0003090287480000321
式中:
Δ、
Figure GDA0003090287480000322
为回归误差计算算子;
Δpx、Δpy、Δpz、Δvx、Δvy、Δvz为在当前轨道根数下,卫星状态分量的回归误差;
Figure GDA0003090287480000323
为依次给定的六个轨道根数微分量;
Figure GDA0003090287480000324
为在新轨道根数卫星状态分量的回归误差。
其中
Figure GDA0003090287480000325
为上述回归轨道参数微分迭代法的设计初值。基于高阶地球重力场的微分迭代法计算过程与基于J2摄动项的微分迭代法近似,不同的是前者需要计算状态矢量的差分以及六个轨道要素的微分修正量,而得到微分迭代收敛结果,记为
Figure GDA0003090287480000326
步骤8:基于高阶地球重力场的回归轨道参数多目标求解
将卫星的位置回归误差ΔrE和速度回归误差
Figure GDA0003090287480000327
描述为两个优化目标,对
Figure GDA0003090287480000328
进行进一步优化设计。建立优化问题模型为:
优化目标:
Figure GDA0003090287480000329
优化变量:
x=(ao,eo,iooo,Mo)T
约束条件:
Figure GDA00030902874800003210
Figure GDA00030902874800003211
式中:
H:在地球非球引力摄动影响下,卫星位置和速度满足的非线性函数关系;
Figure GDA0003090287480000331
轨道要素的寻优范围的上、下限。
G(x)表示:优化目标函数;
ΔrE(x)表示:卫星的位置回归误差;
步骤9:基于NGSA-II算法构建高阶地球重力场下回归轨道参数求解
采用NSGA-II算法对步骤8中的多目标优化问题进行求解,具体过程如下:
1)初始化变量k=0;
2)输入仿真参数:轨道预报的开始时刻ts和结束时刻te,轨道预报积分步长h0,轨道六要素优化范围[xl,xu],遗传算法参数(包括最大寻优代数NG,种群大小NP,交叉概率PC,变异概率PM);
3)随机生成包含NP个个体的初始化种群P0
4)计算初始种群P0中每个个体的目标函数,对P0进行非支配排序,初始化每个个体的等级值和拥挤距离;
5)基于二元锦标赛选择策略从第k代种群Pk中随机选择个体,进行交叉和变异操作,产生后代种群Qk
6)合并Pk和Qk,产生出组合种群Rk=Pk∪Qk
7)对Rk进行快速非支配排序,基于等级值和拥挤距离选出NP个个体,组成新一代种群Pk+1
8)更新遗传代数k=k+1,判断是否满足终止条件:若k>NG,则令
Figure GDA0003090287480000332
并结束循环;否则,跳转至步骤5)。
设计出在高阶地球重力场下满足设计目标的回归轨道参数,记为
Figure GDA0003090287480000333
本领域技术人员知道,除了以纯计算机可读程序代码方式实现本发明提供的系统、装置及其各个模块以外,完全可以通过将方法步骤进行逻辑编程来使得本发明提供的系统、装置及其各个模块以逻辑门、开关、专用集成电路、可编程逻辑控制器以及嵌入式微控制器等的形式来实现相同程序。所以,本发明提供的系统、装置及其各个模块可以被认为是一种硬件部件,而对其内包括的用于实现各种程序的模块也可以视为硬件部件内的结构;也可以将用于实现各种功能的模块视为既可以是实现方法的软件程序又可以是硬件部件内的结构。
以上对本发明的具体实施例进行了描述。需要理解的是,本发明并不局限于上述特定实施方式,本领域技术人员可以在权利要求的范围内做出各种变化或修改,这并不影响本发明的实质内容。在不冲突的情况下,本申请的实施例和实施例中的特征可以任意相互组合。

Claims (5)

1.一种考虑高阶重力场的严格回归轨道设计方法,其特征在于,包括:
卫星轨道初设计步骤:依据卫星轨道设计总体输入,包括回归天数和回归圈数,进行考虑J2摄动影响的回归轨道参数解析设计,获得解析设计结果,通过牛顿迭代法对解析设计结果进行微分修正,获得J2项摄动卫星轨道初设计初值;
严格回归轨道设计步骤:根据获得的严格回归轨道优化设计初值,将牛顿迭代法与多目标优化算法NSGA-II相结合,通过高精度重力场下的回归轨道混合优化设计算法,实现严格回归轨道的快速优化设计;
所述卫星轨道初设计步骤:
轨道半长轴初步确定步骤:根据回归天数N和回归轨道数R,在二体运动模型下,初步确定轨道半长轴
Figure FDA0003090287470000011
J2摄动下的卫星半长轴计算步骤:根据获得的初步确定的轨道半长轴
Figure FDA0003090287470000012
引入J2项摄动,得到J2摄动下的卫星半长轴
Figure FDA0003090287470000013
轨道倾角获取步骤:根据获得的J2摄动下的卫星半长轴
Figure FDA0003090287470000014
再根据太阳同步约束,求得轨道倾角
Figure FDA0003090287470000015
升交点赤经获取步骤:根据分离点的经纬度
Figure FDA0003090287470000016
和分离时刻TS,可求得升交点赤经Ω;
所述严格回归轨道设计步骤:
基于牛顿迭代法的回归轨道回归特性修正步骤:
根据获得的J2项摄动卫星轨道初设计初值半长轴参数a、轨道倾角i以及升交点赤经Ω,设定其余轨道参数,为优化初设计轨道的回归特性,提出基于牛顿迭代法的严格回归重力场修正方法:
设回归后的纬度
Figure FDA0003090287470000017
和回归后的经度λ为半长轴a和轨道倾角i的函数:
Figure FDA0003090287470000018
λ=f(a,i)
对上式求导数,并且在整个回归周期内积分,从回归周期的开始时刻T0到回归周期的结束时刻Tf=T0+Tperiod,可得到下式
Figure FDA0003090287470000021
其中Δ代表在整个回归周期的差分运算符,Tperiod为回归周期;
Figure FDA0003090287470000022
Figure FDA0003090287470000023
为相对于半长轴a和轨道倾角i的偏导数;
所述严格回归轨道设计步骤还包括:
基于微分修正的严格回归轨道修正步骤:
依据牛顿迭代算法计算原理,以卫星轨道根数为变量,构建基于PV=[px,py,pz,vx,vy,vz]状态向量函数为:
Figure FDA0003090287470000024
其中,
px,py,pz,vx,vy,vz分别表示x方向位置、y方向位置、z方向位置、x方向速度、y方向速度、z方向速度;
a,e,i,w,Ω,M分别表示轨道半长轴、偏心率、轨道倾角、近地点幅角、升交点赤经、平近点角;
fpx表示x方向位置关于轨道根数的函数;
fpy表示y方向位置关于轨道根数的函数;
fpz表示z方向位置关于轨道根数的函数;
fvx表示x方向速度关于轨道根数的函数;
fvy表示y方向速度关于轨道根数的函数;
fvz表示z方向速度关于轨道根数的函数;
依据上述状态向量函数,通过微分处理,可得高阶重力场环境下卫星轨道六要素牛顿迭代方程为
Figure FDA0003090287470000031
式中:
Δ、
Figure FDA0003090287470000032
为回归误差计算算子;
Δpx、Δpy、Δpz、Δvx、Δvy、Δvz为在当前轨道根数下,卫星状态分量的回归误差;
Figure FDA0003090287470000033
为依次给定的六个轨道根数微分量;
Figure FDA0003090287470000034
为在新轨道根数卫星状态分量的回归误差;
Δa,Δe,Δi,Δω,ΔΩ,ΔM分别表示轨道半长轴、偏心率、轨道倾角、近地点幅角、升交点赤经、平近点角轨道6要素的微元;
基于NSGA-II的严格回归轨道多目标优化设计步骤:
将严格回归轨道设计问题描述为多目标优化问题,该问题可描述为:
优化目标:min(ΔXp),min(ΔXv),
优化变量:Δα=[Δa,Δe,Δi,Δω,ΔΩ,ΔM]
约束条件:
x′=f(t,x);
Δxp=xp(Tf)-xp(T0);
Δxv=xv(Tf)-xv(T0);
其中,
x′表示卫星状态变量的导数;
Δxp表示回归周期始末的位置差;
Δxv表示回归周期始末的速度差;
f为高阶重力场影响下的卫星轨道状态动力学方程;xp(T0),xp(Tf)分别为回归周期始末地固系下的位置矢量;xv(T0),xv(Tf)分别为回归周期始末地固系下的速度矢量;
为求解上述多目标优化问题,选择NSGA-II多目标优化算法进行求解,以求取最终的卫星轨道根数。
2.根据权利要求1所述的考虑高阶重力场的严格回归轨道设计方法,其特征在于,所述轨道半长轴初步确定步骤:
在二体运动模型下基于回归天数N和回归轨道数R初步确定轨道半长轴
Figure FDA0003090287470000041
Figure FDA0003090287470000042
Figure FDA0003090287470000043
其中,
Figure FDA0003090287470000044
表示初步确定的轨道半长轴;
Figure FDA0003090287470000045
为地球引力常数;
P表示回归系数的倒数;
N表示回归天数;
R表示回归轨道数;
所述J2摄动下的卫星半长轴计算步骤:
Figure FDA0003090287470000046
其中,
Figure FDA0003090287470000047
表示J2摄动下的卫星半长轴;
J2表示J2摄动项;
Figure FDA0003090287470000051
表示升交点赤经的进动值;
Figure FDA0003090287470000052
为地球赤道半径;
所述轨道倾角获取步骤:
Figure FDA0003090287470000053
Figure FDA0003090287470000054
其中,
J2表示J2摄动项;
year表示一年时间对应的秒数;
Figure FDA0003090287470000055
表示轨道倾角;
所述升交点赤经获取步骤:
假设S为分离点,为参考轨道上的某一特定时间TS,即分离时刻,其余的卫星轨道参数可由球面几何计算得到,给定纬度
Figure FDA0003090287470000056
和飞行方向,北面为正:
Figure FDA0003090287470000057
其中,
u表示纬度幅角;
Figure FDA0003090287470000058
表示纬度;
Figure FDA0003090287470000059
表示轨道倾角;
sign表示符号函数;
从S点到下一个升交点的飞行时间可表示为轨道周期的一部分,基于卫星自转,升交点赤经沿经度方向的偏移为:
Figure FDA00030902874700000510
其中,
ΔλN表示升交点赤经沿经度方向的偏移;
Tday表示一天时间对应的秒数,数值为86400;
基于球面几何,升交点赤经与S点的经度差可表示为:
Figure FDA0003090287470000061
其中,
Δλ表示升交点赤经与S点的经度差;
i表示轨道倾角;
给定S点的经度值λ,轨道节点经度可由下式计算
λN=λ-Δλ+ΔλN
其中,
λ表示分离点S的经度值;
λN表示轨道节点经度;
最后,针对给定的时刻,升交点赤经可表示为
Ω=GMST(tUT1)+λN
其中,
Ω表示升交点赤经;
GMST表示格林尼治时角;
tUT1表示平太阳时。
3.一种考虑高阶重力场的严格回归轨道设计系统,其特征在于,包括:
卫星轨道初设计模块:依据卫星轨道设计总体输入,包括回归天数和回归圈数,进行考虑J2摄动影响的回归轨道参数解析设计,获得解析设计结果,通过牛顿迭代法对解析设计结果进行微分修正,获得J2项摄动卫星轨道初设计初值;
严格回归轨道设计模块:根据获得的严格回归轨道优化设计初值,将牛顿迭代法与多目标优化算法NSGA-II相结合,通过高精度重力场下的回归轨道混合优化设计算法,实现严格回归轨道的快速优化设计;
所述卫星轨道初设计模块:
轨道半长轴初步确定模块:根据回归天数N和回归轨道数R,在二体运动模型下,初步确定轨道半长轴
Figure FDA0003090287470000071
J2摄动下的卫星半长轴计算模块:根据获得的初步确定的轨道半长轴
Figure FDA0003090287470000072
引入J2项摄动,得到J2摄动下的卫星半长轴
Figure FDA0003090287470000073
轨道倾角获取模块:根据获得的J2摄动下的卫星半长轴
Figure FDA0003090287470000074
再根据太阳同步约束,求得轨道倾角
Figure FDA0003090287470000075
升交点赤经获取模块:根据分离点的经纬度
Figure FDA0003090287470000076
和分离时刻TS,可求得升交点赤经Ω;
所述严格回归轨道设计模块:
基于牛顿迭代法的回归轨道回归特性修正模块:
根据获得的J2项摄动卫星轨道初设计初值半长轴参数a、轨道倾角i以及升交点赤经Ω,设定其余轨道参数,为优化初设计轨道的回归特性,提出基于牛顿迭代法的严格回归重力场修正方法:
设回归后的纬度
Figure FDA0003090287470000077
和回归后的经度λ为半长轴a和轨道倾角i的函数:
Figure FDA0003090287470000078
λ=f(a,i)
对上式求导数,并且在整个回归周期内积分,从回归周期的开始时刻T0到回归周期的结束时刻Tf=T0+Tperiod,可得到下式
Figure FDA0003090287470000079
其中Δ代表在整个回归周期的差分运算符,Tperiod为回归周期;
Figure FDA00030902874700000710
Figure FDA00030902874700000711
为相对于半长轴a和轨道倾角i的偏导数;
基于微分修正的严格回归轨道修正模块:
依据牛顿迭代算法计算原理,以卫星轨道根数为变量,构建基于PV=[px,py,pz,vx,vy,vz]状态向量函数为:
Figure FDA0003090287470000081
其中,
px,py,pz,vx,vy,vz分别表示x方向位置、y方向位置、z方向位置、x方向速度、y方向速度、z方向速度;
a,e,i,w,Ω,M分别表示轨道半长轴、偏心率、轨道倾角、近地点幅角、升交点赤经、平近点角;
fpx表示x方向位置关于轨道根数的函数;
fpy表示y方向位置关于轨道根数的函数;
fpz表示z方向位置关于轨道根数的函数;
fvx表示x方向速度关于轨道根数的函数;
fvy表示y方向速度关于轨道根数的函数;
fvz表示z方向速度关于轨道根数的函数;
依据上述状态向量函数,通过微分处理,可得高阶重力场环境下卫星轨道六要素牛顿迭代方程为
Figure FDA0003090287470000082
式中:
Δ、
Figure FDA0003090287470000083
为回归误差计算算子;
Δpx、Δpy、Δpz、Δvx、Δvy、Δvz为在当前轨道根数下,卫星状态分量的回归误差;
Figure FDA0003090287470000091
为依次给定的六个轨道根数微分量;
Figure FDA0003090287470000092
为在新轨道根数卫星状态分量的回归误差;
Δa,Δe,Δi,Δω,ΔΩ,ΔM分别表示轨道半长轴、偏心率、轨道倾角、近地点幅角、升交点赤经、平近点角轨道6要素的微元;
基于NSGA-II的严格回归轨道多目标优化设计模块:
将严格回归轨道设计问题描述为多目标优化问题,该问题可描述为:
优化目标:min(ΔXp),min(ΔXv)
优化变量:Δα=[Δa,Δe,Δi,Δω,ΔΩ,ΔM]
约束条件:
x′=f(t,x);
Δxp=xp(Tf)-xp(T0);
Δxv=xv(Tf)-xv(T0);
其中,
x′表示卫星状态变量的导数;
Δxp表示回归周期始末的位置差;
Δxv表示回归周期始末的速度差;
f为高阶重力场影响下的卫星轨道状态动力学方程;xp(T0),xp(Tf)分别为回归周期始末地固系下的位置矢量;xv(T0),xv(Tf)分别为回归周期始末地固系下的速度矢量;
为求解上述多目标优化问题,选择NSGA-II多目标优化算法进行求解,以求取最终的卫星轨道根数。
4.根据权利要求3所述的考虑高阶重力场的严格回归轨道设计系统,其特征在于,所述轨道半长轴初步确定模块:
在二体运动模型下基于回归天数N和回归轨道数R初步确定轨道半长轴
Figure FDA0003090287470000101
Figure FDA0003090287470000102
Figure FDA0003090287470000103
其中,
Figure FDA0003090287470000104
表示初步确定的轨道半长轴;
Figure FDA0003090287470000105
为地球引力常数;
P表示回归系数的倒数;
N表示回归天数;
R表示回归轨道数;
所述J2摄动下的卫星半长轴计算模块:
Figure FDA0003090287470000106
其中,
Figure FDA0003090287470000107
表示J2摄动下的卫星半长轴;
J2表示J2摄动项;
Figure FDA0003090287470000108
表示升交点赤经的进动值;
Figure FDA0003090287470000109
为地球赤道半径;
所述轨道倾角获取模块:
Figure FDA00030902874700001010
Figure FDA00030902874700001011
其中,
J2表示J2摄动项;
year表示一年时间对应的秒数;
Figure FDA00030902874700001012
表示轨道倾角;
所述升交点赤经获取模块:
假设S为分离点,为参考轨道上的某一特定时间TS,即分离时刻,其余的卫星轨道参数可由球面几何计算得到,给定纬度
Figure FDA0003090287470000111
和飞行方向,北面为正:
Figure FDA0003090287470000112
其中,
u表示纬度幅角;
Figure FDA0003090287470000113
表示纬度;
Figure FDA0003090287470000114
表示轨道倾角;
sign表示符号函数;
从S点到下一个升交点的飞行时间可表示为轨道周期的一部分,基于卫星自转,升交点赤经沿经度方向的偏移为:
Figure FDA0003090287470000115
其中,
ΔλN表示升交点赤经沿经度方向的偏移;
Tday表示一天时间对应的秒数,数值为86400;
基于球面几何,升交点赤经与S点的经度差可表示为:
Figure FDA0003090287470000116
其中,
Δλ表示升交点赤经与S点的经度差;
i表示轨道倾角;
给定S点的经度值λ,轨道节点经度可由下式计算
λN=λ-Δλ+ΔλN
其中,
λ表示分离点S的经度值;
λN表示轨道节点经度;
最后,针对给定的时刻,升交点赤经可表示为
Ω=GMST(tUT1)+λN
其中,
Ω表示升交点赤经;
GMST表示格林尼治时角;
tUT1表示平太阳时。
5.一种存储有计算机程序的计算机可读存储介质,其特征在于,所述计算机程序被处理器执行时实现权利要求1至2中任一项所述的考虑高阶重力场的严格回归轨道设计方法的步骤。
CN201910640736.8A 2019-07-16 2019-07-16 一种考虑高阶重力场的严格回归轨道设计方法、系统及介质 Active CN110378012B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910640736.8A CN110378012B (zh) 2019-07-16 2019-07-16 一种考虑高阶重力场的严格回归轨道设计方法、系统及介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910640736.8A CN110378012B (zh) 2019-07-16 2019-07-16 一种考虑高阶重力场的严格回归轨道设计方法、系统及介质

Publications (2)

Publication Number Publication Date
CN110378012A CN110378012A (zh) 2019-10-25
CN110378012B true CN110378012B (zh) 2021-07-16

Family

ID=68253430

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910640736.8A Active CN110378012B (zh) 2019-07-16 2019-07-16 一种考虑高阶重力场的严格回归轨道设计方法、系统及介质

Country Status (1)

Country Link
CN (1) CN110378012B (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111381256B (zh) * 2020-03-10 2022-07-26 上海卫星工程研究所 主动遥感卫星天线相位中心偏移误差计算的方法和系统
CN111551182B (zh) * 2020-06-08 2021-09-17 中国人民解放军战略支援部队信息工程大学 一种可见光室内终端定位方法及可见光定位系统
CN111731513B (zh) * 2020-06-15 2022-03-04 航天东方红卫星有限公司 一种基于单脉冲轨控的高精度引力场中回归轨道维持方法
CN112257343B (zh) * 2020-10-22 2023-03-17 上海卫星工程研究所 一种高精度地面轨迹重复轨道优化方法及系统
CN113734468B (zh) * 2021-08-30 2023-02-03 北京宇航系统工程研究所 一种基于迭代制导的轨道面精确控制方法
CN116562067B (zh) * 2023-07-12 2023-09-19 银河航天(北京)网络技术有限公司 用于构建严格回归目标轨道库的方法、设备和存储介质

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2002039616A2 (en) * 2000-11-04 2002-05-16 Space Resource International Corporation Virtual geostationary satellite constellation and method of satellite communications
CN104484493A (zh) * 2014-10-29 2015-04-01 中国人民解放军63920部队 一种飞船返回预定落点回归轨道设计方法
CN106021784A (zh) * 2016-05-31 2016-10-12 北京航空航天大学 一种基于两层优化策略的全轨迹优化设计方法
CN106092105A (zh) * 2016-06-03 2016-11-09 上海航天控制技术研究所 一种近地卫星严格回归轨道的确定方法
CN107065025A (zh) * 2017-01-13 2017-08-18 北京航空航天大学 一种基于重力梯度不变量的轨道要素估计方法
CN107065930A (zh) * 2017-06-01 2017-08-18 上海航天控制技术研究所 一种复杂约束严格回归轨道控制方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105160417B (zh) * 2015-08-04 2019-05-14 大连大学 基于改进nsga-ii算法的航天器任务规划求解方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2002039616A2 (en) * 2000-11-04 2002-05-16 Space Resource International Corporation Virtual geostationary satellite constellation and method of satellite communications
CN104484493A (zh) * 2014-10-29 2015-04-01 中国人民解放军63920部队 一种飞船返回预定落点回归轨道设计方法
CN106021784A (zh) * 2016-05-31 2016-10-12 北京航空航天大学 一种基于两层优化策略的全轨迹优化设计方法
CN106092105A (zh) * 2016-06-03 2016-11-09 上海航天控制技术研究所 一种近地卫星严格回归轨道的确定方法
CN107065025A (zh) * 2017-01-13 2017-08-18 北京航空航天大学 一种基于重力梯度不变量的轨道要素估计方法
CN107065930A (zh) * 2017-06-01 2017-08-18 上海航天控制技术研究所 一种复杂约束严格回归轨道控制方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Multi-objective Optimization Method For Repeat Ground-track Orbit Design Considering the Orbit Injection Error;Mingjun Pu 等;《Journal of Aerospace Technology and Management》;20180610;第7-8页和图3 *
卫星轨道扫描覆盖仿真程序功能设计;谭振华;《中国优秀硕士学位论文全文数据库 工程科技Ⅱ辑》;20131215(第12期);全文 *
基于迭代修正方法的严格回归轨道设计;杨盛庆 等;《宇航学报》;20160430;第37卷(第4期);全文 *

Also Published As

Publication number Publication date
CN110378012A (zh) 2019-10-25

Similar Documents

Publication Publication Date Title
CN110378012B (zh) 一种考虑高阶重力场的严格回归轨道设计方法、系统及介质
CN112257343B (zh) 一种高精度地面轨迹重复轨道优化方法及系统
Song et al. Solar-sail trajectory design for multiple near-Earth asteroid exploration based on deep neural networks
CN104298647B (zh) 基于低轨道地球卫星的地影时刻预报的星上确定方法
Zuo et al. A case learning-based differential evolution algorithm for global optimization of interplanetary trajectory design
CN108875244B (zh) 一种基于随机森林的轨道预报精度改进方法
CN110553653B (zh) 基于多源数据驱动的航天器轨道确定方法
CN109211246B (zh) 不确定环境下行星着陆轨迹规划方法
CN109032176A (zh) 一种基于微分代数的地球同步轨道确定和参数确定方法
CN110816896B (zh) 一种卫星星上简易轨道外推方法
CN107526368A (zh) 一种考虑误差的多脉冲环月卫星编队初始化方法
CN111731513B (zh) 一种基于单脉冲轨控的高精度引力场中回归轨道维持方法
Lee et al. Robust position and attitude control for spacecraft formation flying
CN115081343B (zh) 一种基于神经网络结合遗传算法的天基无源探测定轨方法
He et al. Solution set calculation of the Sun-perturbed optimal two-impulse trans-lunar orbits using continuation theory
Samsam et al. Nonlinear Model Predictive Control of J2-perturbed impulsive transfer trajectories in long-range rendezvous missions
Yan et al. ANN-based method for fast optimization of Jovian-moon gravity-assisted trajectories in CR3BP
CN111814313A (zh) 一种高精度引力场中回归轨道设计方法
Zhai et al. Improvement of orbit prediction accuracy using extreme gradient boosting and principal component analysis
Ren et al. Research on satellite orbit prediction based on neural network algorithm
Cao et al. HE2LM-AD: Hierarchical and efficient attitude determination framework with adaptive error compensation module based on ELM network
Ananthasayanam Tuning of the Kalman filter using constant gains
CN115392540A (zh) 一种用于月球轨道交会制导的快速预报方法
CN113128032B (zh) 一种基于轨道解析摄动解的交点时刻及位置求解算法
CN110579784B (zh) 基于卫星组合导航系统的卫星自主导航方法

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