CN111591469B - 一种低轨星座***相位保持方法、***、设备和存储介质 - Google Patents

一种低轨星座***相位保持方法、***、设备和存储介质 Download PDF

Info

Publication number
CN111591469B
CN111591469B CN202010140847.5A CN202010140847A CN111591469B CN 111591469 B CN111591469 B CN 111591469B CN 202010140847 A CN202010140847 A CN 202010140847A CN 111591469 B CN111591469 B CN 111591469B
Authority
CN
China
Prior art keywords
satellite
phase angle
time
deviation
semi
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
CN202010140847.5A
Other languages
English (en)
Other versions
CN111591469A (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.)
CASIC Space Engineering Development Co Ltd
Original Assignee
CASIC Space Engineering Development Co Ltd
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 CASIC Space Engineering Development Co Ltd filed Critical CASIC Space Engineering Development Co Ltd
Priority to CN202010140847.5A priority Critical patent/CN111591469B/zh
Publication of CN111591469A publication Critical patent/CN111591469A/zh
Application granted granted Critical
Publication of CN111591469B publication Critical patent/CN111591469B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • BPERFORMING OPERATIONS; TRANSPORTING
    • B64AIRCRAFT; AVIATION; COSMONAUTICS
    • B64GCOSMONAUTICS; VEHICLES OR EQUIPMENT THEREFOR
    • B64G1/00Cosmonautic vehicles
    • B64G1/22Parts of, or equipment specially adapted for fitting in or to, cosmonautic vehicles
    • B64G1/24Guiding or controlling apparatus, e.g. for attitude control
    • B64G1/242Orbits and trajectories

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Chemical & Material Sciences (AREA)
  • Combustion & Propulsion (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)
  • Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)

Abstract

本发明公开了一种低轨星座***相位保持方法、***、设备和存储介质,其中,该方法的步骤包括:根据卫星的参考平相位角和实际平相位角,确定卫星实际半长轴和参考半长轴的偏差以及偏差随时间的变化量;根据卫星实际半长轴和参考半长轴的偏差以及偏差随时间的变化量,确定卫星轨道机动时间和半长轴改变量;利用卫星轨道机动时间和半长轴改变量,在卫星相位保持任务周期内调整低轨星座***的相位。本申请所述技术方案可以用于实现低轨星座***高精度相位保持,能够在卫星位置速度确定精度不高、而且定位数据采样率较低的情况下,实现高精度相位保持(精度要求例如±0.1°),降低了高精度相位保持对卫星定轨精度以及定轨数据采样存储要求。

Description

一种低轨星座***相位保持方法、***、设备和存储介质
技术领域
本发明涉及卫星控制领域。更具体地,涉及一种用于低轨星座***的高精度相位保持方法、***、设备和存储介质。
背景技术
近年来,多个国家正在加速研制和部署大规模低轨星座***。随着大规模低轨星座的部署运行,星座***将面临***内部卫星间以及***内部卫星与空间碎片的碰撞威胁。由于大气阻力等空间摄动的影响,卫星的相位将相对于参考相位发生长期变化,从而对空间轨道安全产生影响,为使低轨星座满足空间轨道安全要求,低轨星座的卫星相位需要保持在一定的较小区间内(例如参考相位的±0.1°以内)。此外,低轨星座的覆盖范围约束也要求卫星相位尽量保持在一定的较小区间内。
目前国内对卫星星座相位保持的研究主要针对中高轨星座,低轨卫星相位保持研究较少。虽然,在低轨卫星相位保持研究方面,国内一些团队利用大气阻力实现低轨卫星星座的相对相位保持,仿真表明空气阻力进行相对相位保持可以实现稳定性很好的星座相对相位保持。通过对目标卫星实施一次相对基准卫星的轨道高度抬升/降低,可以维持星间的相对位置关系,在实际工程应用也表明了此策略的有效性。但是,以上无论是采用大气阻力被动控制,还是采用推力进行主动控制的低轨星座相位保持策略,位置保持的精度量级都超过了±1°,不符合低轨星座***高精度相位保持要求(例如±0.1°)。
发明内容
本发明的一个目的在于提供一种用于低轨星座***的高精度相位保持方法、***、设备和存储介质。
为达到上述目的,本方案采用下述技术方案:
第一方面,本发明提供一种低轨星座***相位保持方法,该方法的步骤包括:
根据卫星的参考平相位角和实际平相位角,确定卫星实际半长轴和参考半长轴的偏差以及偏差随时间的变化量;
根据卫星实际半长轴和参考半长轴的偏差以及偏差随时间的变化量,确定卫星轨道机动时间和半长轴改变量;
利用卫星轨道机动时间和半长轴改变量,在卫星相位保持任务周期内调整低轨星座***的相位。
在一种优选地实施例中,所述根据卫星的参考平相位角和实际平相位角,确定卫星实际半长轴和参考半长轴的偏差以及偏差随时间的变化量的步骤包括:
将卫星同一时刻的实际平相位角和参考相位角作差,获得卫星的平相位角偏差;
根据卫星半长轴随时间的变化趋势和平相位角偏差随时间的变化趋势,确定平相位角偏差随时间变化的因子参数初值;所述相位角偏差随时间变化的因子参数包括:初始时刻的平相位角偏差、初始时刻的半长轴偏差导致平相位角漂移一阶项和半长轴随时间线性变化导致平相位角漂移二阶项;
若某一时刻平相位角偏差预报值超出卫星的漂移边界,则以平相位角偏差随时间变化的因子参数初值,作为确定卫星轨道机动时间和半长轴改变量的依据。
在一种优选地实施例中,所述根据卫星的参考平相位角和实际平相位角,确定卫星实际半长轴和参考半长轴的偏差以及偏差随时间的变化量的步骤还包括:
若某一时平相位角偏差预报值未超出卫星的漂移边界,则重新确定该时刻的实际平相位角和该时刻的卫星平相位角偏差,并根据重新确定的该时刻的卫星平相位角偏差更新平相位角偏差随时间变化的因子参数,直至该时刻的平相位角度偏差超出漂移边界。
在一种优选地实施例中,所述参考平相位角的确定方式为:
根据地球中心引力场和地球非球形摄动,确定卫星的参考平相位角随时间的一阶变化率;
根据卫星在卫星相位保持任务周期的初始时刻时的参考相位角和卫星的参考平相位角随时间的一阶变化率,确定卫星在某一时刻的参考相位角。
在一种优选地实施例中,所述实际平相位角的确定方式为:从卫星相位保持任务周期的初始时刻开始,在某一段时间内每隔固定时长采样一次卫星的位置和速度,并将其转化为实际平相位角。
在一种优选地实施例中,根据卫星实际半长轴和参考半长轴的偏差以及偏差随时间的变化量,确定卫星轨道机动时间和半长轴改变量的步骤包括:
根据平相位角偏差随时间变化的因子参数,确定当前相位保持任务周期的初始时刻的实际半长轴与参考半长轴的偏差;
根据轨道半长轴随时间的一阶线性变化率,确定当前相位保持任务周期的自由漂移结束时刻的实际半长轴和参考半长轴的偏差;
根据平相位角偏差区间和半长轴随时间线性变化导致平相位角漂移二阶项,确定卫星下一相位保持两次轨道高度抬升之间的时长;
根据两次轨道高度抬升之间的半长轴漂移量,确定轨道机动的半长轴改变量。
在一种优选地实施例中,该方法的步骤还包括:
将当前卫星相位保持任务周期的初始时为上一卫星相位保持任务周期的结束时刻。
第二方面,本方案提供一种低轨星座***相位保持***,该***包括:
偏差计算模块,根据卫星的参考平相位角和实际平相位角,确定卫星实际半长轴和参考半长轴的偏差以及偏差随时间的变化量;
控制量计算模块,根据卫星实际半长轴和参考半长轴的偏差以及偏差随时间的变化量,确定卫星轨道机动时间和半长轴改变量;
调整模块,利用卫星轨道机动时间和半长轴改变量,在卫星相位保持任务周期内调整低轨星座***的相位。
第三方面,本方案提供一种设备,包括:存储器,一个或多个处理器;存储器与处理器通过通信总线相连;处理器被配置为执行存储器中的指令;所述存储器中存储有用于执行如上所述方法中各个步骤的指令。
第四方案,本方案提供一种计算机可读存储介质,其上存储有计算机程序,该程序被处理器执行时实现如上所述方法的步骤。
本发明的有益效果如下:
本申请所述技术方案可以用于实现低轨星座***高精度相位保持,能够在卫星位置速度确定精度不高、而且定位数据采样率较低的情况下,实现高精度相位保持(精度要求例如±0.1°),降低了高精度相位保持对卫星定轨精度以及定轨数据采样存储要求。
附图说明
为了更清楚地说明本发明实施例中的技术方案,下面将对实施例描述中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其他的附图。
图1示出本方案所述低轨星座***相位保持方法的示意图。
图2示出本方案实施例所述不施加控制的平相位角偏差的时间历程:2-a为仿真时间10天的示意图;2-b为仿真时间内10个轨道周期。
图3示出本方案实施例所述相位保持任务周期内,平相位角偏差的瞬时轨道值和平均轨道值的时间历程:3-a为仿真时间365天的结果;3-b为仿真60天的局部结果。
图4示出本方案实施例所述相位保持任务周期内,平相位角偏差的瞬时轨道值的时间历程:4-a为仿真时间365天的结果;4-b为仿真时间60天的局部结果。
图5示出本方案实施例所述相位保持任务周期内,平相位角偏差的平均轨道值的时间历程:5-a为仿真时间365天的结果;5-b为仿真时间60天的局部结果。
具体实施方式
为使本发明的、技术方案和优点更加清楚,下面将结合附图对本发明实施方式作进一步地详细描述。显然,所描述的实施例仅是本申请的一部分实施例,而不是所有实施例的穷举。需要说明的是,在不冲突的情况下,本申请中的实施例及实施例中的特征可以相互组合。
经过对现有技术的分析和研究,低轨星座***的相位漂移主要是由大气阻力导致的卫星轨道高度降低造成的,为了实现相位保持,低轨卫星的半长轴保持在参考半长轴附近,因此需要尽量高精度地确定实际半长轴与参考半长轴的偏差。采用间隔较短的长时间地面测控数据和GPS数据进行精密定轨,可以提高实际半长轴的定轨精度,但计算量较大,对于大规模低轨星座***,实现难度较大。因此,本方案提出一种方案:通过卫星实际平相位角与参考平相位角偏差随时间的变化趋势,精确确定实际半长轴与参考半长轴偏差的初始值以及偏差随时间的变化趋势,用于确定相位保持轨道高度改变量,实现低轨星座***高精度相位保持。本方案提出的方法对卫星的位置、速度确定精度要求不高,典型GPS的定位精度即可满足要求;而且不需要对卫星位置、速度数据进行高频率采样,仅仅需要间隔一定时间的位置、速度数据作为相位保持的输入,例如每隔一天的位置、速度数据。
以下,对本方案提出的一种用于低轨星座***的高精度相位保持方法进行详细说明。如图1所示,该方法主要包括如下三个步骤:
根据卫星的参考平相位角和实际平相位角,确定卫星实际半长轴和参考半长轴的偏差以及偏差随时间的变化量中;
根据卫星实际半长轴和参考半长轴的偏差以及偏差随时间的变化量,确定卫星轨道机动时间和半长轴改变量;
利用卫星轨道机动时间和半长轴改变量,在卫星相位保持任务周期内调整低轨星座***的相位。
在第一个实施步骤中,可以先将卫星同一时刻的实际平相位角和参考相位角作差,获得卫星的平相位角偏差;
根据卫星半长轴随时间的变化趋势和平相位角偏差随时间的变化趋势,确定平相位角偏差随时间变化的因子参数初值;所述相位角偏差随时间变化的因子参数包括:初始时刻的平相位角偏差、初始时刻的半长轴偏差导致平相位角漂移一阶项和半长轴随时间线性变化导致平相位角漂移二阶项;
一个可能的实现方式是,同一时刻卫星的实际平相位角与其相应的参考平相位角作差,获得时刻t=ti+l×Δtgc的平相位角偏差,并将卫星的平相位角偏差记为
Δν(ti+l×Δtgc)=νreal(ti+l×Δtgc)-νref(ti+l×Δtgc) (1)
在大气阻力的作用下,半长轴随时间的长期变化趋势为一阶线性变化,平相位角偏差随时间的变化趋势是接近二阶非线性的,将平相位角偏差Δν(t)随时间t的变化趋势近似为
Δv(t)≈Δv0+Δv1×(t-ti)+Δν2×(t-ti)2 (2)
其中:Δv0是初始时刻ti的平相位角偏差,Δν1是初始时刻ti的半长轴偏差导致平相位角漂移一阶项,Δv2是半长轴随时间线性变化导致平相位角漂移二阶项,这三项即为平相位角偏差随时间变化的因子参数。
根据时刻t=ti+l×Δtgc,卫星的平相位角偏差Δν(ti+l×Δtgc)数值,以及表达式(3),利用最小二乘法计算因子参数Δν0、Δν1和Δν2的初值。
在第一个实施步骤中,所述参考平相位角的确定方式为:根据地球中心引力场和地球非球形摄动,确定卫星的参考平相位角随时间的一阶变化率;根据卫星在卫星相位保持任务周期的初始时刻时的参考相位角和卫星的参考平相位角随时间的一阶变化率,确定卫星在某一时刻的参考相位角。
所在第一个实施步骤中,确定参考平相位角的一种可能实现方式为:
1、确定卫星的参考平相位角随时间的一阶变化率
本专利中的平相位角ν,为近地点幅角ω与平近点角M之和,即v=ω+M。通过本专利的控制策略将参考平相位角与实际平相位角的偏差Δv控制在区间 [Δνmin Δvmax]内。在大气阻力的作用下卫星轨道高度逐渐降低,轨道角速度逐渐变大,在平相位角偏差接近Δνmax时,实施轨道抬高机动,将卫星的平相位角偏差Δν控制在区间[Δνmin Δνmax]内。
参考平相位角随时间的变化仅考虑地球中心引力场和地球非球形摄动低阶项导致的长期变化。将地球中心引力场和地球非球形摄动低阶项导致的平相位角v随时间的一阶变化率记为nref。一阶变化率nref可以根据考虑的地球中心引力场和地球非球形摄动低阶项直接解析获得,本专利仅给出考虑地球中心引力场和地球非球形摄动J2 项的一阶变化率
Figure GDA0003304977800000061
其中:轨道角速度
Figure GDA0003304977800000062
a为半长轴;半通径p=a(1-e2);e为轨道偏心率i为轨道倾角。μe是地球的引力常数,μe=3.986×1014m3/s2;Req是地球的平均赤道半径,Req=6378.140km;J2=1.08263×10-3
不妨假设任务开始时刻,卫星的平相位角已经调整到保持区间内,将整个任务的开始时刻作为第一个相位保持周期的初始时刻。
2、确定相位保持周期内参考平相位角随时间的变化关系
将卫星在初始时刻ti的参考平相位角记为vref(ti),则卫星在时刻t的参考平相位角为
vref(t)=νref(ti)+nref×(t-ti) (4)
在第一个实施步骤中,所述实际平相位角的确定方式为:从卫星相位保持任务周期的初始时刻开始,在某一段时间内每隔固定时长采样一次卫星的位置和速度,并将其转化为实际平相位角。
所在第一个实施步骤中,确定实际平相位角的一种可能实现方式为:
从初始时刻开始,在某段时间长内每隔固定时长采样一次卫星的位置和速度,并将其转化为实际平相位角
将位置和速度采样值的获取间隔记为Δtgc,从初始时刻开始获取L+1个采样点的数据,从而总采样时长为Δtp=L×Δtgc。从初始时刻开始,每隔时长Δtgc采样一次卫星的位置矢量r和速度矢量v的值,并通过位置、速度到轨道要素的转化方法,计算卫星在时刻t=ti+l×Δtgc(l∈[0,1,2,…,L])的实际平相位角vreal(ti+l×Δtgc)。
在第一个实施步骤中,进一步需要判断某一时刻平相位角偏差预报值超出卫星的漂移边界;若某一时刻平相位角偏差预报值超出卫星的漂移边界,则以平相位角偏差随时间变化的因子参数初值,作为确定卫星轨道机动时间和半长轴改变量的依据;若某一时平相位角偏差预报值未超出卫星的漂移边界,则重新确定该时刻的实际平相位角和该时刻的卫星平相位角偏差,并根据重新确定的该时刻的卫星平相位角偏差更新平相位角偏差随时间变化的因子参数,直至该时刻的平相位角度偏差超出漂移边界。
该步骤的一种可能实现方式为:
根据平相位角偏差随时间变化的因子参数Δν0、Δν1和Δν2的初值,可以预报时刻ti+(L+1)×Δtgc的平相位角偏差。如果时刻ti+(L+1)×Δtgc的平相位角偏差预报值超出漂移边界Δνmax,转到0确定相位保持轨道控制量;否则,在时间段 [ti+(L)×Δtgc ti+(L+1)×Δtgc]内不施加相位保持轨道控制,采用0的方法获得时刻ti+(L+1)×Δtgc的实际平相位角,并利用0的方法确定时刻ti+(L+1)×Δtgc平相位角偏差,在考虑时刻ti+(L+1)×Δtgc平相位角偏差基础上,对平相位角偏差随时间变化的因子参数进行更新。
采用上述步骤依次预报时刻ti+(L+2)×Δtgc、ti+(L+3)×Δtgc、……的平相位角偏差,直到相应时刻的平相位角偏差超出漂移边界Δνmax,记录相应的从初始时刻ti开始的卫星无控制自由漂移时间(L+H)×Δtgc,以及该段时间内的平相位角偏差随时间变化的因子参数。
在第二个实施步骤中,可以根据卫星实际半长轴和参考半长轴的偏差以及偏差随时间的变化量,确定卫星轨道机动时间和半长轴改变量的步骤包括:
根据平相位角偏差随时间变化的因子参数,确定当前相位保持任务周期的初始时刻的实际半长轴与参考半长轴的偏差;
根据轨道半长轴随时间的一阶线性变化率,确定当前相位保持任务周期的自由漂移结束时刻的实际半长轴和参考半长轴的偏差;
根据平相位角偏差区间和半长轴随时间线性变化导致平相位角漂移二阶项,确定卫星下一相位保持两次轨道高度抬升之间的时长;
根据两次轨道高度抬升之间的半长轴漂移量,确定轨道机动的半长轴改变量。
该步骤的一种可能实现方式为:
根据卫星从时刻ti到时刻ti+(L+H)×Δtgc的平相位角偏差随时间变化的因子参数Δν0、Δν1和Δν2。通过下式可以计算初始时刻ti实际半长轴与参考半长轴的偏差Δa(ti)。
Figure GDA0003304977800000081
通过下式可以计算这段时间段内卫星的轨道半长轴随时间的一阶线性变化率
Figure GDA0003304977800000084
Figure GDA0003304977800000082
从而可以确定ti+(L+H)·Δtgc的实际半长轴与参考半长轴的偏差为
Figure GDA0003304977800000083
根据平相位角偏差区间[Δνmin Δvmax]以及平相位角漂移二阶项Δv2,可以估算下一相位保持两次轨道高度抬升之间的时长为
Figure GDA0003304977800000091
在两次轨道高度抬升之间的半长轴漂移量为
Figure GDA0003304977800000092
所以轨道机动的半长轴改变量(即轨道高度改变量)为
Figure GDA0003304977800000093
对于近圆轨道,以上半长轴改变量需要的横向速度增量为
Figure GDA0003304977800000094
在接下来的时间长Δtgc内,即时刻ti+(L+H+1)·Δtgc前,施加切向推力,获得横向速度增量ΔvT
相位保持周期的总时长定义为(L+H+1)·Δtgc,其中在前时长(L+H)·Δtgc内,卫星进行自由漂移;在最后一个时长Δtgc内,卫星完成轨道抬高机动。
为了配合本方案所述低轨星座***相位保持方法的应用,本方案进一步提供了一种低轨星座***相位保持***,该***包括:偏差计算模块、控制量计算模块和调整模块。偏差计算模块根据卫星的参考平相位角和实际平相位角,确定卫星实际半长轴和参考半长轴的偏差以及偏差随时间的变化量;再利用控制量计算模块根据卫星实际半长轴和参考半长轴的偏差以及偏差随时间的变化量,确定卫星轨道机动时间和半长轴改变量;最后,通过调整模块利用卫星轨道机动时间和半长轴改变量,在卫星相位保持任务周期内调整低轨星座***的相位。
此外,本方案所述低轨星座***相位保持方法还可以通过电子设备来实现相应的功能。其中,该设备包括:存储器,一个或多个处理器;存储器与处理器通过通信总线相连;处理器被配置为执行存储器中的指令;所述存储器中存储有用于执行如上所述低轨星座***相位保持方法中各个步骤的指令。此外,为了方便工作人员使用和操作,也可以进一步为设备配置显示器,以实现人机交互。
另外,本方案所述低轨星座***相位保持方法还可以通过将方法编写为计算机程序存储计算机可读存储介质中;该程序被处理器执行时实现如上所述低轨星座***相位保持方法的步骤。
实施例
实施例1
本实施例提供了一种用于低轨星座***的高精度相位保持方法,该方法对卫星的位置、速度确定精度要求不高,典型GPS的定位精度即可满足要求;而且不需要对卫星位置、速度数据进行高频率采样,仅仅需要间隔一定时间的位置、速度数据作为相位保持的输入,例如每隔一天的位置、速度数据。对于星座内的每颗卫星,相位保持策略相同,仅仅参考相位不同,通过本专利的控制策略,使每颗卫星的实际平相位角跟踪其参考平相位角。由于每颗卫星的相位保持策略相同,后续不再针对星座展开描述,直接针对卫星展开描述。此外,本方案对于星座内的每颗卫星,相位保持策略相同,仅仅参考相位不同,通过本专利的控制策略,使每颗卫星的实际平相位角跟踪其参考平相位角。由于每颗卫星的相位保持策略相同,后续不再针对星座展开描述,直接针对卫星展开描述。具体的,该方法的实施过程如下:
步骤一、确定卫星的参考平相位角随时间的一阶变化率
本专利中的平相位角v,为近地点幅角ω与平近点角M之和,即v=ω+M。通过本专利的控制策略将参考平相位角与实际平相位角的偏差Δv控制在区间 [Δνmin Δvmax]内。在大气阻力的作用下卫星轨道高度逐渐降低,轨道角速度逐渐变大,在平相位角偏差接近Δvmax时,实施轨道抬高机动,将卫星的平相位角偏差Δν控制在区间[Δνmin Δvmax]内。
参考平相位角随时间的变化仅考虑地球中心引力场和地球非球形摄动低阶项导致的长期变化。将地球中心引力场和地球非球形摄动低阶项导致的平相位角v随时间的一阶变化率记为nref。一阶变化率nref可以根据考虑的地球中心引力场和地球非球形摄动低阶项直接解析获得,本专利仅给出考虑地球中心引力场和地球非球形摄动J2 项的一阶变化率
Figure GDA0003304977800000111
其中:轨道角速度
Figure GDA0003304977800000112
a为半长轴;半通径p=a(1-e2);e为轨道偏心率i为轨道倾角。μe是地球的引力常数,μe=3.986×1014m3/s2;Req是地球的平均赤道半径,Req=6378.140km;J2=1.08263×10-3
不妨假设任务开始时刻,卫星的平相位角已经调整到保持区间内,将整个任务的开始时刻作为第一个相位保持周期的初始时刻。
步骤二、确定相位保持周期内参考平相位角随时间的变化关系
将卫星在初始时刻ti的参考平相位角记为vref(ti),则卫星在时刻t的参考平相位角为
νref(t)=νref(ti)+nref×(t-ti) (13)
步骤三、从初始时刻开始,在某段时间长内每隔固定时长采样一次卫星的位置和速度,并将其转化为实际平相位角
将位置和速度采样值的获取间隔记为Δtgc,从初始时刻开始获取L+1个采样点的数据,从而总采样时长为Δtp=L×Δtgc。从初始时刻开始,每隔时长Δtgc采样一次卫星的位置矢量r和速度矢量v的值,并通过位置、速度到轨道要素的转化方法,计算卫星在时刻t=ti+l×Δtgc(l∈[0,1,2,…,L])的实际平相位角νreal(ti+l×Δtgc)。
步骤四、计算平相位角偏差随时间变化的因子参数初值
将以上相应时刻卫星的实际平相位角与其相应的参考平相位角作差,获得时刻 t=ti+l×Δtgc的平相位角偏差,并将卫星的平相位角偏差记为
Δv(ti+l×Δtgc)=νreal(ti+l×Δtgc)-vref(ti+l×Δtgc) (14)
在大气阻力的作用下,半长轴随时间的长期变化趋势为一阶线性变化,平相位角偏差随时间的变化趋势是接近二阶非线性的,将平相位角偏差Δv(t)随时间t的变化趋势近似为
Δv(t)≈Δv0+Δv1×(t-ti)+Δv2×(t-ti)2 (15)
其中:Δv0是初始时刻ti的平相位角偏差,Δv1是初始时刻ti的半长轴偏差导致平相位角漂移一阶项,Δv2是半长轴随时间线性变化导致平相位角漂移二阶项,这三项即为平相位角偏差随时间变化的因子参数。
根据时刻t=ti+l×Δtgc,卫星的平相位角偏差Δv(ti+l×Δtgc)数值,以及表达式(3),利用最小二乘法计算因子参数Δv0、Δv1和Δv2的初值。
步骤五、确定平相位角偏差值超出漂移边界对应的时刻以及卫星无控时间段内的平相位角偏差随时间变化的因子参数
以上确定了平相位角偏差随时间变化的因子参数Δν0、Δv1和Δv2的初值,从而可以预报时刻ti+(L+1)×Δtgc的平相位角偏差。如果时刻ti+(L+1)×Δtgc的平相位角偏差预报值超出漂移边界Δvmax,转到0确定相位保持轨道控制量;否则,在时间段 [ti+(L)×Δtgcti+(L+1)×Δtgc]内不施加相位保持轨道控制,采用0的方法获得时刻ti+(L+1)×Δtgc的实际平相位角,并利用0的方法确定时刻ti+(L+1)×Δtgc平相位角偏差,在考虑时刻ti+(L+1)×Δtgc平相位角偏差基础上,对平相位角偏差随时间变化的因子参数进行更新。
采用上述步骤依次预报时刻ti+(L+2)×Δtgc、ti+(L+3)×Δtgc、……的平相位角偏差,直到相应时刻的平相位角偏差超出漂移边界Δvmax,记录相应的从初始时刻ti开始的卫星无控制自由漂移时间(L+H)×Δtgc,以及该段时间内的平相位角偏差随时间变化的因子参数。
步骤六、确定卫星轨道机动时间和半长轴改变量。
以上确定了卫星从时刻ti到时刻ti+(L+H)×Δtgc的平相位角偏差随时间变化的因子参数Δv0、Δv1和Δν2。通过下式可以计算初始时刻ti实际半长轴与参考半长轴的偏差Δa(ti)。
Figure GDA0003304977800000121
通过下式可以计算这段时间段内卫星的轨道半长轴随时间的一阶线性变化率
Figure GDA0003304977800000123
Figure GDA0003304977800000122
从而可以确定ti+(L+H)×Δtgc的实际半长轴与参考半长轴的偏差为
Figure GDA0003304977800000131
根据平相位角偏差区间[Δνmin Δνmax]以及平相位角漂移二阶项Δν2,可以估算下一相位保持两次轨道高度抬升之间的时长为
Figure GDA0003304977800000132
在两次轨道高度抬升之间的半长轴漂移量为
Figure GDA0003304977800000133
所以轨道机动的半长轴改变量(即轨道高度改变量)为
Figure GDA0003304977800000134
对于近圆轨道,以上半长轴改变量需要的横向速度增量为
Figure GDA0003304977800000135
在接下来的时间长Δtgc内,即时刻ti+(L+H+1)×Δtgc前,施加切向推力,获得横向速度增量ΔvT
相位保持周期的总时长定义为(L+H+1)×Δtgc,其中在前时长(L+H)×Δtgc内,卫星进行自由漂移;在最后一个时长Δtgc内,卫星完成轨道抬高机动。
步骤七、将时刻ti+(L+H+1)×Δtgc作为下一控制周期的初始时刻,并重复步骤 0至0,直到完成卫星整个相位保持任务周期的相位保持。
实施例2
本实施例中,根据仿真实验对本方案提出一种用于低轨星座***的高精度相位保持方法作进一步说明。
首先分析轨道要素的选择,图2为不施加控制的平相位角偏差的时间历程,其中,图2-a为仿真时间10天的示意图;图2-b为仿真时间内10个轨道周期。图中平均轨道要素指的是在瞬时轨道要素的基础上,通过解析方法,去除地球非球形摄动J2导致的短周期项。
可以看出:平相位角偏差的瞬时轨道要素的振荡幅度约为±0.05°,对应的振荡周期约为0.5个轨道周期;平相位角偏差的瞬时轨道值的振荡幅度约为±0.05°,对应的振荡周期约为0.5个轨道周期;平相位角偏差的拟平均轨道值(仅考虑去掉J2项导致的短周期项,可以通过单点数据解析获得)的振荡幅度约为±0.005°,对应的振荡周期约为0.5天。为了减小短周期项对轨道控制的影响,将平均轨道要素作为相位保持控制的输入。
假设任务要求将平相位角偏差的瞬时轨道值保持在±0.1°的区间内,由于平相位角偏差的瞬时轨道值的振荡幅度约为±0.05°,所以考虑该项振荡后,平相位角偏差的拟平均轨道值应保持在±0.05°的区间内。由于平相位角偏差的拟平均轨道值的振荡幅度约为±0.005°,而平相位角偏差的拟平均轨道值作为相位保持的控制时,不确定采样点处于最高点还是最低点,考虑该项振荡后,平相位角偏差的拟平均轨道值应保持在±0.04°的区间内。本专利的相位保持策略没有考虑日月第三体引力以及太阳光压的影响,这些摄动力也会造成卫星平相位角的漂移,综合考虑未建模的摄动力,将实施例采样点处的平相位角偏差的拟平均轨道值保持在区间[-0.02° 0.04°]内,即Δvmin=-0.02°,Δvmax=0.04°。
除了特别说明,后续实施例中的平相位角偏差都是指的拟平均轨道值。星座***的半长轴、偏心率、轨道倾角的标称值分别为a=7528137m、e1=0.001、i=88°。将位置和速度采样值的获取间隔取为Δtgc=1day,卫星的初始质量取为200kg,太阳帆板面积取为3m2,L取为10,位置三维分量误差均取为30m,速度的三维分量误差均取为0.01m/s,相位保持任务周期取为365天。以下根据发明内容给出星座内一颗卫星的相位保持仿真结果,其他卫星可以根据发明内容进行相应计算。
步骤一、确定卫星的参考平相位角随时间的一阶变化率
在半长轴、偏心率、轨道倾角给定的情况下,一阶变化率nref可以根据考虑的地球中心引力场和地球非球形摄动低阶项直接解析获得
nref=5.532°×10-2[sec-1]
在任务开始时刻t=0,平相位角的参考值为99.98°,不妨设卫星的平相位角已经调整到控制区间内,符合这一条件的一组参数瞬时轨道参数如下表所示,将整个任务的开始时刻作为第一个相位保持周期的初始时刻,即ti=0。
表1卫星在任务开始时刻的初始瞬时轨道要素
瞬时轨道要素 卫星
半长轴,m 7519972
偏心率 5.9181×10<sup>-4</sup>
轨道倾角,deg 86.9983
升交点赤经,deg 19.9994
近地点幅角,deg 311.1214
平近点角,deg 148.8416
步骤二、确定相位保持周期内参考平相位角随时间的变化关系
卫星在初始时刻ti的参考平相位角为vref(ti)=99.98°,则卫星在时刻t的参考平相位角为
vref(t)=99.98°+5.532°×10-2×(t-ti)
步骤三、从初始时刻开始,在某段时间长内每隔固定时长采样一次卫星的位置和速度,并将其转化为实际平相位角。
从初始时刻开始,每隔时长Δtgc=1day采样一次卫星的位置矢量r和速度矢量v的值,并通过位置、速度到平均轨道要素的转化方法,计算卫星在时刻t=l×86400sec (l∈[0,1,2,…,11])的实际平相位角νreal(l×86400sec)。
步骤四、计算平相位角偏差随时间变化的因子参数初值
将以上相应时刻卫星的实际平相位角与其相应的参考平相位角作差,获得不同仿真时刻的平相位角偏差如下表所示。
表2不同仿真时刻的实际平相位角与参考平相位角之差
Figure GDA0003304977800000151
Figure GDA0003304977800000161
根据平相位角偏差随时间t的变化规律,并结合上表的数据可以计算平相位角偏差随时间变化的因子参数初值分别为:
Figure GDA0003304977800000162
步骤五、确定平相位角偏差值超出漂移边界对应的时刻以及卫星无控时间段内的平相位角偏差随时间变化的因子参数
根据平相位角偏差随时间变化的因子参数,可以预测t=11×86400sec的平相位角为-0.009395°,没有超出边界值0.04°,在时间段[(10)×86400(11)×86400]内不施加相位保持轨道控制,继续获得时刻(11)×86400sec的实际平相位角,并确定时刻 (11)×86400sec平相位角偏差以及考虑时刻(11)×86400sec平相位角偏差的平相位角偏差随时间变化的因子参数。
采用上述步骤依次预报时刻(12)×86400sec、(13)×86400sec的平相位角偏差,直到相应时刻的平相位角偏差超出漂移边界Δνlimit,确定从时刻ti开始的卫星无控制自由漂移时间为(10+12)×86400sec,即H=12,以及该段时间内的平相位角偏差随时间变化的因子参数
Figure GDA0003304977800000163
步骤六、确定卫星轨道机动时间和半长轴改变量。
计算初始时刻实际半长轴与参考半长轴的偏差。
Δa(0)≈2.416m
计算这段时间段内卫星的轨道半长轴随时间的一阶线性变化率
Figure GDA0003304977800000172
Figure GDA0003304977800000171
从而可以确定时刻(10+12)×86400sec的实际半长轴与参考半长轴的偏差为
Δa((10+12)×86400sec)≈-5.725m (23)
根据平相位角偏差区间[-0.02° 0.04°]以及平相位角漂移二阶项Δν2=4.979×10-6deg/day2,可以估算下一相位保持两次轨道高度抬升之间的时长为
Δt=37.713day
在两次轨道高度抬升之间的半长轴漂移量为
Δad=-13.349m
所以轨道机动的半长轴改变量(即轨道高度改变量)为
δamane=12.400m
对于近圆轨道,以上半长轴改变量需要的横向速度增量为
ΔvT=5.993×10-3m/s
在时刻23×86400sec前,施加切向推力,获得以上横向速度增量。
步骤七、将时刻23×86400sec作为下一控制周期的初始时刻,并重复步骤上述步骤,直到完成卫星整个相位保持任务周期的相位保持。
通过图3中相位保持任务周期内,平相位角偏差的瞬时轨道值和平均轨道值的时间历程:3-a为仿真时间365天的结果;3-b为仿真60天的局部结果。图4中相位保持任务周期内,平相位角偏差的瞬时轨道值的时间历程:4-a为仿真时间365天的结果;4-b为仿真时间60天的局部结果。图5中相位保持任务周期内,平相位角偏差的平均轨道值的时间历程:5-a为仿真时间365天的结果;5-b为仿真时间60天的局部结果。可以看出平相位角偏差的瞬时轨道值保持在±0.1°以内,平相位角偏差的平均轨道值保持在±0.05°以内,通过本方案的控制策略实现了高精度的相位保持。
显然,本发明的上述实施例仅仅是为清楚地说明本发明所作的举例,而并非是对本发明的实施方式的限定,对于所属领域的普通技术人员来说,在上述说明的基础上还可以做出其它不同形式的变化或变动,这里无法对所有的实施方式予以穷举,凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。

Claims (8)

1.一种低轨星座***相位保持方法,其特征在于,该方法的步骤包括:
根据卫星的参考平相位角和实际平相位角,确定卫星实际半长轴和参考半长轴的偏差以及偏差随时间的变化量,步骤包括:
将卫星同一时刻的实际平相位角和参考平相位角作差,获得卫星的平相位角偏差,
根据卫星半长轴随时间的变化趋势和平相位角偏差随时间的变化趋势,确定平相位角偏差随时间变化的因子参数初值,
所述相位角偏差随时间变化的因子参数包括:初始时刻的平相位角偏差、初始时刻的半长轴偏差导致平相位角漂移一阶项和半长轴随时间线性变化导致平相位角漂移二阶项,
若某一时刻平相位角偏差预报值超出卫星的漂移边界,则以平相位角偏差随时间变化的因子参数初值,作为确定卫星轨道机动时间和半长轴改变量的依据;
根据卫星实际半长轴和参考半长轴的偏差以及偏差随时间的变化量,确定卫星轨道机动时间和半长轴改变量,步骤包括:
根据平相位角偏差随时间变化的因子参数,确定当前相位保持任务周期的初始时刻的实际半长轴与参考半长轴的偏差,
根据轨道半长轴随时间的一阶线性变化率,确定当前相位保持任务周期的自由漂移结束时刻的实际半长轴和参考半长轴的偏差,
根据平相位角偏差区间和半长轴随时间线性变化导致平相位角漂移二阶项,确定卫星下一相位保持两次轨道高度抬升之间的时长,
根据两次轨道高度抬升之间的半长轴漂移量,确定轨道机动的半长轴改变量;
利用卫星轨道机动时间和半长轴改变量,在卫星相位保持任务周期内调整低轨星座***的相位。
2.根据权利要求1所述的方法,其特征在于,所述根据卫星的参考平相位角和实际平相位角,确定卫星实际半长轴和参考半长轴的偏差以及偏差随时间的变化量的步骤还包括:
若某一时平相位角偏差预报值未超出卫星的漂移边界,则重新确定该时刻的实际平相位角和该时刻的卫星平相位角偏差,并根据重新确定的该时刻的卫星平相位角偏差更新平相位角偏差随时间变化的因子参数,直至该时刻的平相位角度偏差超出漂移边界。
3.根据权利要求1所述的方法,其特征在于,所述参考平相位角的确定方式为:
根据地球中心引力场和地球非球形摄动,确定卫星的参考平相位角随时间的一阶变化率;
根据卫星在卫星相位保持任务周期的初始时刻时的参考相位角和卫星的参考平相位角随时间的一阶变化率,确定卫星在某一时刻的参考相位角。
4.根据权利要求1所述的方法,其特征在于,所述实际平相位角的确定方式为:从卫星相位保持任务周期的初始时刻开始,在某一段时间内每隔固定时长采样一次卫星的位置和速度,并将其转化为实际平相位角。
5.根据权利要求1至4任一项所述的方法,其特征在于,该方法的步骤还包括:
将当前卫星相位保持任务周期的初始时为上一卫星相位保持任务周期的结束时刻。
6.一种低轨星座***相位保持***,其特征在于,该***包括:
偏差计算模块,根据卫星的参考平相位角和实际平相位角,确定卫星实际半长轴和参考半长轴的偏差以及偏差随时间的变化量;
控制量计算模块,根据卫星实际半长轴和参考半长轴的偏差以及偏差随时间的变化量,确定卫星轨道机动时间和半长轴改变量;
调整模块,利用卫星轨道机动时间和半长轴改变量,在卫星相位保持任务周期内调整低轨星座***的相位。
7.一种低轨星座***相位保持设备,其特征在于,包括:存储器,一个或多个处理器;存储器与处理器通过通信总线相连;处理器被配置为执行存储器中的指令;所述存储器中存储有用于执行如权利要求1至5任一项所述方法中各个步骤的指令。
8.一种计算机可读存储介质,其上存储有计算机程序,其特征在于,该程序被处理器执行时实现如权利要求1至5任一项所述方法的步骤。
CN202010140847.5A 2020-03-03 2020-03-03 一种低轨星座***相位保持方法、***、设备和存储介质 Active CN111591469B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010140847.5A CN111591469B (zh) 2020-03-03 2020-03-03 一种低轨星座***相位保持方法、***、设备和存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010140847.5A CN111591469B (zh) 2020-03-03 2020-03-03 一种低轨星座***相位保持方法、***、设备和存储介质

Publications (2)

Publication Number Publication Date
CN111591469A CN111591469A (zh) 2020-08-28
CN111591469B true CN111591469B (zh) 2021-12-07

Family

ID=72181316

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010140847.5A Active CN111591469B (zh) 2020-03-03 2020-03-03 一种低轨星座***相位保持方法、***、设备和存储介质

Country Status (1)

Country Link
CN (1) CN111591469B (zh)

Families Citing this family (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112124626B (zh) * 2020-08-27 2022-02-15 中国人民解放军战略支援部队航天工程大学 一种Walker星座构型维持方法和终端设备
CN112769466B (zh) * 2020-12-22 2022-08-12 火眼位置数智科技服务有限公司 低轨卫星星座构型保持方法
CN113189619B (zh) * 2021-04-02 2023-05-09 中国空间技术研究院 一种低轨星座相位保持参数估计方法
CN113447957B (zh) * 2021-06-23 2022-11-25 北京航空航天大学 一种高精度实时卫星轨道瞬时半长轴的去振方法
CN113602534B (zh) * 2021-06-26 2023-02-28 山东航天电子技术研究所 一种微型电推进推力大小的在轨标定方法
CN113697129B (zh) * 2021-07-16 2023-06-06 北京航空航天大学 一种星座相对相位控制的方法及计算机设备
CN113859587B (zh) * 2021-09-26 2023-06-02 中国西安卫星测控中心 一种巨型低轨星座卫星绝对站位保持方法
CN114248949B (zh) * 2021-11-22 2022-12-27 浙江大学 异面质比星座的部署方法及装置、电子设备、存储介质
CN115806060B (zh) * 2022-11-10 2023-05-19 北京航天驭星科技有限公司 卫星相对相位保持策略模型的建模方法、模型、获取方法
CN115892516B (zh) * 2022-11-10 2023-06-20 北京航天驭星科技有限公司 卫星相对相位保持策略模型的建模方法、模型、获取方法
CN115806061B (zh) * 2022-11-10 2023-05-09 北京航天驭星科技有限公司 卫星相对相位保持策略模型的建模方法、模型、获取方法
CN116873225B (zh) * 2023-08-07 2024-01-23 北京航天驭星科技有限公司 采用电推方式维持相位差的方法、***、电子设备和介质
CN117311394B (zh) * 2023-11-28 2024-03-12 中国西安卫星测控中心 一种非严格多星串联编队***中星间相位超越控制方法
CN117508648B (zh) * 2024-01-05 2024-04-26 北京航天驭星科技有限公司 同一轨道面星座卫星的轨控方法、装置、设备和介质

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102880184A (zh) * 2012-10-24 2013-01-16 北京控制工程研究所 一种静止轨道卫星自主轨道控制方法
CN107031868A (zh) * 2017-03-23 2017-08-11 北京空间飞行器总体设计部 一种低轨遥感卫星自主轨道控制方法
CN110053788A (zh) * 2019-03-15 2019-07-26 中国西安卫星测控中心 一种考虑复杂摄动的星座长期保持控制频次估计方法
EP3587282A1 (en) * 2018-06-21 2020-01-01 Airbus Oneweb Satellites SAS Satellite control apparatuses and methods

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US3588000A (en) * 1968-07-29 1971-06-28 Westinghouse Electric Corp Orbit position control of planar reflector satellites
US6453220B1 (en) * 2000-01-31 2002-09-17 Space Systems/Loral, Inc. Low earth orbit satellite constellation stationkeeping algorithm with absolute altitude control
US6672542B2 (en) * 2002-06-03 2004-01-06 The Aerospace Corporation Method and system for controlling the eccentricity of a near-circular orbit
CN101794154A (zh) * 2009-11-25 2010-08-04 哈尔滨工业大学 一种编队卫星相对轨道与姿态解耦控制方法
CN102591343B (zh) * 2012-02-09 2013-10-16 航天东方红卫星有限公司 基于两行根数的卫星轨道维持控制方法
FR3020348B1 (fr) * 2014-04-24 2016-05-13 Snecma Procede de deploiement d'une constellation de satellites
CN106628257B (zh) * 2016-09-28 2018-01-19 西北工业大学 地球摄动引力场中近地航天器相对运动轨道的保持方法
CN107797130B (zh) * 2017-10-16 2021-01-05 中国西安卫星测控中心 低轨航天器多点多参数轨道上行数据计算方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102880184A (zh) * 2012-10-24 2013-01-16 北京控制工程研究所 一种静止轨道卫星自主轨道控制方法
CN107031868A (zh) * 2017-03-23 2017-08-11 北京空间飞行器总体设计部 一种低轨遥感卫星自主轨道控制方法
EP3587282A1 (en) * 2018-06-21 2020-01-01 Airbus Oneweb Satellites SAS Satellite control apparatuses and methods
CN110053788A (zh) * 2019-03-15 2019-07-26 中国西安卫星测控中心 一种考虑复杂摄动的星座长期保持控制频次估计方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
一种考虑摄动影响的星座构型稳定性设计方法;陈长春等;《上海航天(中英文)》;20200225(第01期);全文 *
载人航天器与空间站共轨飞行轨道维持策略;刘欢等;《航天器工程》;20190415(第02期);全文 *

Also Published As

Publication number Publication date
CN111591469A (zh) 2020-08-28

Similar Documents

Publication Publication Date Title
CN111591469B (zh) 一种低轨星座***相位保持方法、***、设备和存储介质
CN109625323B (zh) 一种卫星化学推进变轨方法及***
CN112257343B (zh) 一种高精度地面轨迹重复轨道优化方法及***
CN104142686B (zh) 一种卫星自主编队飞行控制方法
CN103557871B (zh) 一种浮空飞行器捷联惯导空中初始对准方法
CN107797130B (zh) 低轨航天器多点多参数轨道上行数据计算方法
US10046869B2 (en) Inertial sensing augmentation for navigation of spacecraft
CN109080854B (zh) 航天器返回预定落点的大椭圆轨道变轨规划方法
CN110104219B (zh) 一种控制探测器着陆地外天体的方法及装置
CN109344449B (zh) 航天器月地转移轨道逆向设计方法
CN104015938A (zh) 一种电推进静止轨道卫星的位置保持方法
CN108388135B (zh) 一种基于凸优化的火星着陆轨迹优化控制方法
CN106114910B (zh) 一种航天器飞行轨道滚动时域控制方法
CN115639580A (zh) 卫星自主导航方法
CN110789739A (zh) 一种j2摄动下长时间轨道交会最优速度增量快速估计方法
CN112607065B (zh) 一种基于电推进***的高精度相位控制方法
CN110053788B (zh) 一种考虑复杂摄动的星座长期保持控制频次估计方法
CN111301715A (zh) 基于霍曼变轨的同轨道特定相位分布的星座布局与轨道调整方法、装置及计算机存储介质
CN102878997A (zh) 一种大偏心率轨道的星上快速高精度外推方法
CN113525721A (zh) 卫星轨道变换方法、装置、电子设备及存储介质
CN108082538B (zh) 一种考虑始末约束的多体***低能量捕获轨道方法
CN104020678B (zh) 一种基于月表地形的动力下降初始点参数优化方法
CN111060111A (zh) 一种低轨卫星入轨初期定轨方法
CN113311854A (zh) 一种月球采样返回任务中的定点着陆轨道设计方法
CN111428912B (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