CN113483765B - 一种卫星自主姿态确定方法 - Google Patents
一种卫星自主姿态确定方法 Download PDFInfo
- Publication number
- CN113483765B CN113483765B CN202110565164.9A CN202110565164A CN113483765B CN 113483765 B CN113483765 B CN 113483765B CN 202110565164 A CN202110565164 A CN 202110565164A CN 113483765 B CN113483765 B CN 113483765B
- Authority
- CN
- China
- Prior art keywords
- attitude
- ekf
- sensitive
- calculation
- sensor
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 42
- 238000005259 measurement Methods 0.000 claims abstract description 96
- 238000004364 calculation method Methods 0.000 claims abstract description 89
- 238000001914 filtration Methods 0.000 claims abstract description 68
- 239000013598 vector Substances 0.000 claims abstract description 42
- 230000000295 complement effect Effects 0.000 claims abstract description 31
- 230000008569 process Effects 0.000 claims abstract description 8
- 230000035945 sensitivity Effects 0.000 claims description 63
- 239000011159 matrix material Substances 0.000 claims description 27
- 230000005389 magnetism Effects 0.000 claims description 22
- 238000004422 calculation algorithm Methods 0.000 claims description 9
- 238000004590 computer program Methods 0.000 claims description 6
- 238000010586 diagram Methods 0.000 description 12
- 230000003287 optical effect Effects 0.000 description 10
- 238000009434 installation Methods 0.000 description 9
- 239000000758 substrate Substances 0.000 description 6
- 238000011161 development Methods 0.000 description 5
- 230000009466 transformation Effects 0.000 description 5
- RTZKZFJDLAIYFH-UHFFFAOYSA-N Diethyl ether Chemical compound CCOCC RTZKZFJDLAIYFH-UHFFFAOYSA-N 0.000 description 4
- 150000001875 compounds Chemical class 0.000 description 4
- 230000006870 function Effects 0.000 description 4
- 230000036544 posture Effects 0.000 description 4
- 238000012545 processing Methods 0.000 description 4
- 230000004044 response Effects 0.000 description 4
- 238000005070 sampling Methods 0.000 description 4
- 230000009897 systematic effect Effects 0.000 description 4
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 3
- 230000004927 fusion Effects 0.000 description 3
- 238000007499 fusion processing Methods 0.000 description 3
- 206010020751 Hypersensitivity Diseases 0.000 description 2
- 238000012937 correction Methods 0.000 description 2
- 230000009977 dual effect Effects 0.000 description 2
- 238000005286 illumination Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 239000013307 optical fiber Substances 0.000 description 2
- 230000002093 peripheral effect Effects 0.000 description 2
- 230000000644 propagated effect Effects 0.000 description 2
- XKVYZLLWKHGKMT-BEJOYRPXSA-N Gemin D Natural products O([C@@H]([C@@H](O)C=O)[C@@H]1[C@@H](O)COC(=O)c2c(c(O)c(O)c(O)c2)-c2c(O)c(O)c(O)cc2C(=O)O1)C(=O)c1cc(O)c(O)c(O)c1 XKVYZLLWKHGKMT-BEJOYRPXSA-N 0.000 description 1
- 238000003491 array Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 229930192479 gemin Natural products 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000006855 networking Effects 0.000 description 1
- 239000004065 semiconductor Substances 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
Images
Classifications
-
- 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
-
- 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/005—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 with correlation of navigation data from several sources, e.g. map or contour matching
-
- 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/10—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration
- G01C21/12—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning
- G01C21/16—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation
- G01C21/165—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation combined with non-inertial navigation instruments
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Automation & Control Theory (AREA)
- General Physics & Mathematics (AREA)
- Astronomy & Astrophysics (AREA)
- Navigation (AREA)
- Gyroscopes (AREA)
Abstract
本发明一个实施例公开了一种卫星自主姿态确定方法,包括:S100、通过各敏感器测量信息进行敏感器姿态解算得到姿态信息,其中,各敏感器姿态解算包括陀螺递推姿态解算、双星敏姿态解算、单星敏姿态解算、太敏和地敏双矢量姿态解算、太敏和磁双矢量姿态解算以及地敏和磁双矢量解算;S102、将步骤S100中经陀螺和其他敏感器姿态解算得到的姿态信息进行有效组合完成滤波姿态解算,进一步得到连续的姿态信息,其中,滤波姿态解算包括互补滤波姿态解算和扩展卡尔曼滤波姿态解算即EKF姿态解算;S104、制定数据判断选用规则,考虑各敏感器及滤波定姿的精度,设计数据判断流程,在步骤S102的互补滤波姿态解算和EKF姿态解算的结果中选取最优化姿态信息。
Description
技术领域
本发明涉及航天技术领域。更具体地,涉及一种卫星自主姿态确定方法,计算机设备以及计算机可读存储介质。
背景技术
随着航天事业的推进,航天任务呈现灵活多样性的发展势态,小型卫星具有功能密度大、研制和发射成本低、体积小重量轻、研制周期短、发射手段灵活和易于组网等特点,其空间应用需求日渐增加,市场化及小型化成为卫星未来发展的主流趋势。控制系统作为卫星可靠运行的核心部分,获得较为精准的姿态输入极为重要。
为确保卫星任务的可靠执行,卫星通常配备有多种测量敏感器,如星敏感器、太阳敏感器和红外地平仪等。目前,通常采用星敏感器作为姿态确定的主要敏感器,陀螺仪作为基础角速度测量器件,然而星敏不能提供连续的姿态信息,存在由于丢失星图而造成量测不稳定等问题;高精度陀螺则存在结构复杂且重量、体积较大的问题。受小型卫星自身成本、重量等条件的限制,其配备的敏感器均为低成本的中低精度器件。因此,借助中低精度敏感器通过设计实现高性能高可靠性的姿态确定目标,成为小型卫星快速发展的迫切需求。
发明内容
本发明的目的在于提供一种卫星自主姿态确定方法。以解决现有技术存在的问题中的至少一个。
为达到上述目的,本发明采用下述技术方案:
第一方面,本发明提供了一种卫星自主姿态确定方法,包括:
S100、通过各敏感器测量信息进行敏感器姿态解算得到姿态信息,其中,各敏感器姿态解算包括陀螺递推姿态解算、双星敏姿态解算、单星敏姿态解算、太敏和地敏双矢量姿态解算、太敏和磁双矢量姿态解算以及地敏和磁双矢量解算;
S102、将步骤S100中经陀螺和其他敏感器姿态解算得到的姿态信息进行有效组合完成滤波姿态解算,进一步得到连续的姿态信息,其中,所述滤波姿态解算包括互补滤波姿态解算和扩展卡尔曼滤波姿态解算即EKF姿态解算;
S104、制定数据判断选用规则,考虑各敏感器及滤波定姿的精度,设计数据判断流程,在步骤S102的互补滤波姿态解算和EKF姿态解算的结果中选取最优化姿态信息。
在一个具体实施方式中,
在进行所述敏感器姿态解算时,除陀螺递推姿态解算外,依据优先级仅采用一种姿态解算方法,敏感器优先级依次为双星敏、单星敏、太敏和地敏、太敏和磁,或地敏和磁。
在一个具体实施方式中,
所述互补滤波算法的算式为:
qerr=q-1×qt
q=q×q_tmp
式中,q为陀螺递推解算得到的姿态四元数;qt为星敏、太敏和地敏、太敏和磁、或地敏和磁解算得到的姿态四元数;qerr为两者姿态四元数的差值;k为互补滤波增益系数;q_tmp为qerr的归一化值。
在一个具体实施方式中,
在进行所述EKF姿态解算时,分为四种EKF模式,包括:星敏+陀螺EKF、太敏和地敏+陀螺EKF、太敏和磁+陀螺EKF、地敏和磁+陀螺EKF,根据滤波有效性判断依据,选择最优滤波方案进行姿态解算到姿态信息。
在一个具体实施方式中,所述进行EKF姿态解算的步骤包括:
S 1020、确定各敏感器测量模型;
S1021、确定各种EKF模式计算方法;
S1022、进行EKF有效性判断,根据预定规则判断四种EKF模式的有效性,选取最优项。
在一个具体实施方式中,
所述预定规则为:判断敏感器失效时长是否超限;若敏感器失效时长未超限,判断EKF是否收敛,如果在一定时间内系统估计均方差阵的范数稳定且小于设定阈值则EKF收敛;若EKF已收敛且星敏感器定姿有效,将EKF结果与星敏感器定姿结果进行比较,误差不超限时认为EKF有效,优先使用滤波结果。
在一个具体实施方式中,
根据所述数据判断选用规则,经过综合计算得到双星敏定姿精度最优,其次为星敏EKF定姿,之后依次为太敏和地敏EKF定姿、单星敏定姿、太敏和地敏定姿、太敏和磁EKF定姿、地敏和磁EKF定姿、太敏和磁定姿或地敏和磁定姿。
在一个具体实施方式中,
若判定选用EKF定姿,则将EKF解算得到的角速度信息、姿态信息赋值到系统相应信息,若采用敏感器定姿,则无需使用EKF解算得到的数据;当敏感器定姿无效且EKF尚未收敛时,仍采用EKF解算的姿态信息。
第二方面,本发明还提供一种计算机可读存储介质,其上存储有计算机程序,该程序被处理器执行时实现如本申请第一方面提供的方法。
第三方面,本发明还提供一种计算机设备,包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序,所述处理器执行所述程序时实现如本申请第一方面提供的方法。
本发明的有益效果如下:
本发明提供的一种卫星自主姿态确定方法,能有效避免因单一敏感器故障导致的姿态信息紊乱,保证姿态的稳定性和可靠性。本发明的一种卫星自主姿态确定方法,借助中低精度敏感器,通过数据融合处理解算得到高精度姿态信息,保证卫星业务运行的需求。
附图说明
为了更清楚地说明本发明实施例中的技术方案,下面将对实施例描述中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其他的附图。
图1示出根据本发明一个实施例的一种卫星自主姿态确定方法流程图。
图2示出根据本发明一个实施例的一种卫星自主姿态确定方法示意图。
图3示出根据本发明一个实施例的卫星轨道系及本体系示意图。
图4示出根据本发明一个实施例的互补滤波姿态解算示意图。
图5示出根据本发明一个实施例的EKF姿态解算流程图。
图6示出根据本发明一个实施例的EKF姿态解算原理图。
图7示出根据本发明一个实施例的EKF有效性判断流程图。
图8示出根据本发明一个实施例的星敏1姿态解算姿态角误差示意图。
图9示出根据本发明一个实施例的陀螺+星敏1互补滤波姿态角误差示意图。
图10示出根据本发明一个实施例的四种EKF形式的姿态角误差示意图
图11示出根据本发明一个实施例的收敛后四种EKF形式的姿态角误差示意图。
图12示出根据本发明一个实施例的姿态数据融合选用流程图。
图13示出适于用来实现本申请实施例的计算机设备的结构示意图。
图14示出根据本发明一个实施例的星敏EKF形式的姿态角误差示意图。
具体实施方式
为使本发明的技术方案和优点更加清楚,下面将结合附图对本发明实施方式作进一步地详细描述。附图中相似的部件以相同的附图标记进行表示。本领域技术人员应当理解,下面所具体描述的内容是说明性的而非限制性的,不应以此限制本发明的保护范围。
第一实施例
如图1和图2所示,本发明的一个实施例公开了一种卫星自主姿态确定方法,包括:
S100、通过各敏感器测量信息进行敏感器姿态解算得到姿态信息,其中,各敏感器姿态解算包括陀螺递推姿态解算、双星敏姿态解算、单星敏姿态解算、太敏和地敏双矢量姿态解算、太敏和磁双矢量姿态解算以及地敏和磁双矢量解算。
具体的,卫星装备的测量敏感器包括1套三轴陀螺仪(以下简称陀螺)、2套星敏感器(以下简称星敏)、1套三轴太阳敏感器(以下简称太敏)、1套地球敏感器(以下简称地敏)和1套三轴磁强计(以下简称磁),各敏感器均属于中低精度等级。设定陀螺测量精度为0.5°/h,星敏1和星敏2测量精度分别为5角秒、10角秒,太敏和地敏测量精度分别为0.05°、0.15°,磁测量精度为500nT。
在本实施例中,定义姿态信息为本体系相对轨道系的姿态;惯性系为赤道惯性坐标系;轨道系OcXoYoZo,取卫星质心为坐标系原点Oc,OcZo轴指向卫星地心矢径的反方向,OcXo轴过Oc点,在与OcZo轴垂直的平面内指向速度方向,OcYo轴与它们按照右手法则螺旋。本体系坐标原点Oc,OcZb轴指向与对地板外法向方向,OcYb过Oc点,在与OcZb垂直的平面内指向卫星标准飞行姿态下的轨道平面负法线方向,OcXb轴与它们按照右手法则螺旋;安装系由各敏感器在卫星本体的安装位置坐标决定。卫星轨道系及本体系示意图如图3所示。
在本实施例中,首先对双矢量定姿原理进行说明。已知在参考坐标系中有两个互不平行的参考矢量V1,V2,它们在本体系中观测矢量为U1,U2。现求姿态转换矩阵A满足U1=AV1,U2=AV2,利用基准矢量的不平行性,分别在参考坐标系和本体系中建立新的正交坐标系R、S,各坐标轴的单位矢量是
在本实施例中,对陀螺递推姿态解算进行说明,陀螺仪输出卫星相对惯性空间的角速度,通过上一时刻惯性系到本体系的四元数,一步递推得到当前时刻惯性系到本体系的四元数qi,结合卫星轨道信息可以得到惯性系到轨道系的四元数qoi,于是可以得到卫星本体在轨道系下的姿态四元数q,进而解算得到姿态信息。
在本实施例中,对双星敏姿态解算、单星敏姿态解算进行说明。
星敏输出其安装坐标系相对惯性空间的姿态四元数,当两个星敏均有效工作时,采用双星敏姿态解算,否则采用单星敏姿态解算。
在双星敏姿态解算方法中,根据星敏输出的姿态四元数可以得到星敏安装系相对惯性系的坐标转换矩阵,即
V1=Fsi1·[0 0 1]T,V2=Fsi2·[0 0 1]T
U1=Ms1·[0 0 1]T,U2=Ms2·[0 0 1]T
之后采用双矢量定姿原理可以得到惯性系到本体系的姿态转换矩阵,进而得到星敏在惯性系的姿态四元数,结合轨道数据最终解算得到姿态信息。
在单星敏姿态解算方法中,借助星敏输出的惯性系的姿态四元数,结合轨道数据即可解算得到姿态信息。
在本实施例中,对太敏和地敏双矢量、太敏和磁双矢量、地敏和磁双矢量姿态解算进行说明。
太敏输出相对其安装坐标系下的太阳矢量Sb;地敏输出相对其安装坐标系下的地心矢量Eb,磁强计输出相对其安装坐标系下的地磁矢量Bb。结合轨道数据中太阳矢量、地心矢量和磁矢量在轨道系下的分量So,Eo,Bo,
U11=Sb,V12=So
U21=Eb,V22=Eo
U31=Bb,V32=Bo
采用双矢量定姿原理即可分别得到轨道系大盘本体系的姿态转换矩阵,进而得到相应姿态四元数,解算姿态信息。
在一个具体实施例中,在进行所述敏感器姿态解算时,除陀螺递推姿态解算外,依据优先级仅采用一种姿态解算方法,敏感器优先级依次为双星敏、单星敏、太敏和地敏、太敏和磁,或地敏和磁。
S102、将步骤S100中经陀螺和其他敏感器姿态解算得到的姿态信息进行有效组合完成滤波姿态解算,进一步得到高精度连续的姿态信息,其中,所述滤波姿态解算包括互补滤波姿态解算和扩展卡尔曼滤波姿态解算即EKF姿态解算;
具体的,互补滤波是一种利用两种或两种以上的输入得到最优化输出的滤波器技术,在卫星各敏感器姿态解算中,陀螺仪动态响应特性良好,但会产生累积误差;星敏、太敏和地敏、太敏和磁、地敏和磁计算姿态没有累积误差,但动态响应特性相对较差。因此,两者在频域上特性互补,可以采用互补滤波算法将两者得到的姿态数据进行融合,提高解算姿态的精度和动态可靠性。
在一个具体实施例中,如图4所示,进行互补滤波姿态解算,互补滤波是一种利用两种或两种以上的输入得到最优化输出的滤波器技术,其输入通常具有不同的频率特性,即频域上具有互补的关系。
输入信号x(t),系统中包含两种频域上互补的观测量,其中一种在高频区域可信度较高,另一种在低频区域可信度高。互补滤波其通过对第一种观测量进行高通滤波,即1-G(t),另一种观测量进行低通滤波,即G(t),结合两种输入得到优于任意输入的输出值y(t),
y(t)=[x(t)+δ(t)][1-G(t)]+[x(t)+η(t)]G(t)
=x(t)+δ(t)[1-G(t)]+η(t)G(t)
式中,δ(t),η(t)为信号噪声。
在卫星各敏感器姿态解算中,陀螺仪动态响应特性良好,但会产生累积误差;星敏、太敏和地敏、太敏和磁,或地敏和磁计算姿态没有累积误差,但动态响应特性相对较差。因此,两者在频域上特性互补,可以采用互补滤波算法将两者得到的姿态数据进行融合,提高解算姿态的精度和动态可靠性。
在本实施例中,所述互补滤波算法的算式为:
qerr=q-1×qt
q=q×q_tmp
式中,q为陀螺递推解算得到的姿态四元数;qt为星敏、太敏和地敏、太敏和磁、或地敏和磁解算得到的姿态四元数;qerr为两者姿态四元数的差值;k为互补滤波增益系数;q_tmp为qerr的归一化值。
具体的,EKF姿态解算是一种状态估计法,通过建立状态量变化的状态方程及观测方程,采用估计算法根据观测信息估计出状态量,并成为一定准则下的最优估计。考虑实际应用中的非线性因素,采用扩展卡尔曼滤波算法进行卫星姿态确定,该算法的核心是利用泰勒展开将非线性的状态方程和量测方程进行局部线性化,之后利用卡尔曼滤波原理进行计算。
如图5所示,EKF算法大体上可以分为两个阶段,即预测阶段和校正阶段。其核心思想为先根据上一时刻系统状态,预测当前时刻系统状态,再结合当前时刻系统状态的观测值修正得到当前时刻“真实”状态。预估过程建立对当前状态的先验估计,校正过程负责反馈,利用当前测量变量校正估计状态,建立起对状态改进的后验估计。
在一个具体实施例中,在EKF姿态解算中,考虑各敏感器使用特性,设计四种EKF模式,所述四种EKF模式包括:星敏+陀螺EKF、太敏和地敏+陀螺EKF、太敏和磁+陀螺EKF、地敏和磁+陀螺EKF,根据滤波有效性判断依据,选择最优滤波方案进行姿态解算到姿态信息。图6为EKF姿态解算原理图。
在一个具体实施例中,所述进行EKF姿态解算的步骤包括:
S1020、确定各敏感器测量模型;
陀螺测量模型
ω=ωb+b+vg
星敏测量模型
星敏测量输出模型为
qbm=qb+vm
式中,qbm为星敏测量值,qb为真实值,vm为星敏测量噪声。记星敏测量残差为星敏测量输出四元数与估计值之间的差值。
Δqb=HmΔq+vm
Hm=I3×3
式中,Hm为星敏量测矩阵。
太敏测量模型
本体系下太阳矢量Sb=[Sbx Sby Sbz],太敏测量角η,ξ,设
Dη=tanη=-Sbx/Sby
Dξ=tanξ=Sbz/Sby
太敏测量输出模型为
式中,vs为太敏测量噪声,记太敏测量残差为太敏测量输出与估计值之间的差值,
式中,Hs为太敏量测矩阵。
地敏测量模型
式中,He为地敏量测矩阵。
磁测量模型
为方便计算,简化磁测模型,
Bbm=Bb+vc
式中,Bbm为本体系下磁测量值,Bb为真实磁矢量,vc为磁测量噪声。记磁测量残差为磁测量输出值与估计值之间的差值。
ΔBb=2[Bbm×]Δq+vc
=HcΔq+vc
式中,Hc为磁测量测矩阵。
S1021、确定各种EKF模式计算方法;
星敏+陀螺EKF
在星敏+陀螺EKF姿态解算中,结合卫星姿态运动学模型,以星敏测量数据作为量测矢量,以姿态四元数作为系统状态量,并将陀螺测量数据引入系统状态方程,可建立对姿态的滤波观测。滤波输出值除姿态的最优估计外,还包括陀螺的常值漂移误差的最优估计。
i.系统状态方程
系统状态量选取为姿态误差四元数和陀螺漂移,即
XT=[Δq1 Δq2 Δq3 Δbx Δby Δbz]T
根据四元数运动方程和陀螺的输出模型,可得误差四元数表示的方程,
Δq0=0
系统状态方程为
式中,F(t),G(t)分别为状态方差阵和系统噪声方差阵,Wk为系统激励噪声序列。本文中后续定义均与此一致。
将系统状态方程离散化
X(k+1)=Φ(k+1,k)·X(k)+Γ(k+1,k)·W
式中,Φ(k)=I6+F(k)·Δt,Γ(k)=G(t)·Δt,Δt为采样步长,Φ(k+1,k)为tk时刻至tk+1时刻的一步转移阵;Γ(k+1,k)为系统噪声驱动阵。本文中后续定义均与此一致。
ii.系统量测方程
系统量测量采用四元数表示,星敏测量残差矢量部分为两个星敏测量的四元数与陀螺积分四元数的差值,记为[Δq1 Δq2]T。
系统量测方程为
Z(t)=H·X(t)+V
太敏和地敏+陀螺EKF
在太敏和地敏+陀螺EKF姿态解算中,结合卫星姿态运动学模型,以太敏和地敏测量数据作为量测矢量,以姿态四元数作为系统状态量,并将陀螺和地敏测量数据引入系统状态方程,可建立对姿态的滤波观测。滤波输出值除姿态的最优估计外,还包括陀螺和地敏常值测量误差的最优估计。
i.系统状态方程
系统状态量选取为姿态误差四元数、陀螺漂移和地敏测量误差,即
XT=[Δq1 Δq2 Δq3 Δbx Δby Δbz Δγbias Δθbias]T
系统滤波状态方程为
将滤波状态方程离散化,
X(k+1)=Φ(k+1,k)·X(k)+Γ(k+1,k)·W
式中,Φ(k)=I8+F(k)·Δt,Γ(k)=G(t)·Δt,Δt为采样步长。
ii.系统量测方程
若卫星处于光照区,太敏和地敏均正常输出,系统量测量采用地敏和太敏的测量残差,记为ZT=[ΔDη ΔDξ Δγ Δθ]T。
系统量测方程为
Z(t)=H·X(t)+V
若卫星处于阴影区,太敏失效,系统量测量采用地敏的测量残差,记为ZT=[0 0Δγ Δθ]T,量测矩阵中Hs=0,vs=0;若地敏失效,系统量测量采用太敏的测量残差,ZT=[ΔDη ΔDξ 0 0]T,量测矩阵中He=0,ve=0。
太敏和磁+陀螺EKF
在太敏和磁+陀螺EKF姿态解算中,结合卫星姿态运动学模型,以太敏和磁测量数据作为量测矢量,以姿态四元数作为系统状态量,并将陀螺测量数据引入系统状态方程,可建立对姿态的滤波观测。滤波输出值除姿态的最优估计外,还包括陀螺常值测量误差的最优估计。
i.系统状态方程
系统状态量选取为姿态误差四元数和陀螺误差,即
XT=[Δq1 Δq2 Δq3 Δbx Δby Δbz]T
系统滤波状态方程为
将滤波状态方程离散化,
X(k+1)=Φ(k+1,k)·X(k)+Γ(k+1,k)·W
式中,Φ(k)=I6+F(k)·Δt,Γ(k)=G(t)·Δt,Δt为采样步长。
ii.系统量测方程
若卫星处于光照区,太敏和磁均正常输出,系统量测量采用太敏和磁的测量残差,记为ZT=[ΔDη ΔDξ ΔBbx ΔBby ΔBbz]T。
系统量测方程为Z(t)=H·X(t)+V
若卫星处于阴影区,太敏失效,系统量测量采用磁测量残差,记为ZT=[0 0 ΔBbxΔBby ΔBbz]T,量测矩阵中Hs=0,vs=0;若磁失效,系统量测量采用太敏的测量残差,ZT=[ΔDη ΔDξ 0 0 0]T,量测矩阵中Hc=0,vc=0。
地敏和磁+陀螺EKF
在地敏和磁+陀螺EKF姿态解算中,结合卫星姿态运动学模型,以地敏和磁测量数据作为量测矢量,以姿态四元数作为系统状态量,并将陀螺和地敏测量数据引入系统状态方程,可建立对姿态的滤波观测。滤波输出值除姿态的最优估计外,还包括陀螺和地敏常值测量误差的最优估计。
i.系统状态方程
系统状态量选取为姿态误差四元数、陀螺漂移和地敏测量误差,即
XT=[Δq1 Δq2 Δq3 Δbx Δby Δbz Δγbias Δθbias]T
系统滤波状态方程为
将滤波状态方程离散化,
X(k+1)=Φ(k+1,k)·X(k)+Γ(k+1,k)·W
式中,Φ(k)=I8+F(k)·Δt,Γ(k)=G(t)·Δt,Δt为采样步长。
ii.系统量测方程
若地敏和磁均有效,系统量测量采用地敏和磁的测量残差,记为ZT=[Δγ ΔθΔBbx ΔBby ΔBbz]T。
系统量测方程为Z(t)=H·X(t)+V
若地敏失效,系统量测量采用磁测量残差,记为ZT=[0 0 ΔBbx ΔBby ΔBbz]T,量测矩阵中Hs=0,vs=0;若磁失效,系统量测量采用地敏的测量残差,ZT=[Δγ Δθ 0 0 0]T,量测矩阵中Hc=0,vc=0。
S1022、进行EKF有效性判断,根据预定规则判断四种EKF模式的有效性,选取最优项。
在滤波姿态解算中,四种EKF模式均可以解算得到姿态信息,需要制定相应规则判断其有效性,选取最优项。在本实施例中,如图7所示,所述预定规则为:判断敏感器失效时长是否超限;若敏感器失效时长未超限,判断EKF是否收敛,如果在一定时间内系统估计均方差阵的范数稳定且小于设定阈值则EKF收敛;若EKF已收敛且星敏感器定姿有效,将EKF结果与星敏感器定姿结果进行比较,误差不超限时认为EKF有效,优先使用滤波结果。
S104、制定数据判断选用规则,考虑各敏感器及滤波定姿的精度,设计数据判断流程,在步骤S102的互补滤波姿态解算和EKF姿态解算的结果中选取最优化姿态信息。姿态数据融合选用流程图如图12所示。
对于星敏而言,其光轴误差相对较大,敏感器双星敏定姿中采取了光轴误差补偿的措施;而为保证星敏滤波中数据最大化使用,滤波观测量为星敏原始数据,并未消除光轴误差,在一个具体实施例中,根据所述数据判断选用规则,经过综合计算得到双星敏定姿精度最优,其次为星敏EKF定姿,之后依次为太敏和地敏EKF定姿、单星敏定姿、太敏和地敏定姿、太敏和磁EKF定姿、地敏和磁EKF定姿、太敏和磁定姿或地敏和磁定姿。
在一个具体实施例中,
若判定选用EKF定姿,则将EKF解算得到的角速度信息、姿态信息赋值到系统相应信息,若采用敏感器定姿,则无需使用EKF解算得到的数据;当敏感器定姿无效且EKF尚未收敛时,仍采用EKF解算的姿态信息。
下面将结合具体示例对本发明进行进一步说明。
在本示例中,卫星处于倾斜圆轨道,轨道高度1125km、倾角86.5°、周期108.8min。卫星装备中低精度的1套陀螺、2套星敏、1套太敏、1套地敏和1套磁强计作为姿态测量敏感器。仿真时长为一个轨道周期6528s,计算周期0.25s。在三轴对地稳定状态下,设初始姿态角为[5°;5°;5°],角速度为0,除星敏2外各敏感器初始设置均为有效。
S100、通过各敏感器测量信息进行敏感器姿态解算得到姿态信息。
在本示例中,在敏感器姿态解算部分采用星敏1姿态解算。星敏1输出安装系相对惯性系的四元数,通过星敏1在本体系下的安装矩阵可以得到本体系相对惯性系的四元数,结合轨道数据即可解算得到姿态信息。
从图8可以看出,星敏1姿态解算的误差在0.03°以内。
S102、将步骤S100中经陀螺和其他敏感器姿态解算得到的姿态信息进行有效组合完成滤波姿态解算,进一步得到高精度连续的姿态信息,其中,所述滤波姿态解算包括互补滤波姿态解算和扩展卡尔曼滤波姿态解算即EKF姿态解算;
具体的,在本示例中考虑敏感器姿态确定精度,在滤波姿态解算部分采用陀螺+星敏1互补滤波姿态解算及星敏EKF滤波姿态解算。
进行互补滤波姿态解算,陀螺+星敏1互补滤波,陀螺递推解算得到的姿态四元数q,星敏1解算得到的姿态四元数qt,增益系数k设置为0.1,
qerr=q-1×qt
q=q×q_tmp
式中,q为陀螺递推解算得到的姿态四元数;qt为星敏、太敏和地敏、太敏和磁、或地敏和磁解算得到的姿态四元数;qerr为两者姿态四元数的差值;k为互补滤波增益系数;q_tmp为qerr的归一化值。
因此,如图9所示,可以得到互补滤波数据融合处理之后的姿态四元数,结合轨道数据最终解算得到姿态信息,误差在0.005°以内。
分别进行四种EKF形式的姿态角解算,结果如图10所示。如图11所示,收敛后,星敏EKF姿态误差在0.001°范围内,其余三项滤波姿态误差在0.01°范围内。
S104、制定数据判断选用规则,考虑各敏感器及滤波定姿的精度,设计数据判断流程,在步骤S102的互补滤波姿态解算和EKF姿态解算的结果中选取最优化姿态信息。
在本示例中,根据图12所示的姿态数据融合选用流程图,最终选择星敏EKF作为姿态解算方式,得到高精度姿态角信息,星敏EKF形式的姿态角误差如图14所示。
本发明提供的一种小型卫星多模冗余自主姿态确定方法,能有效避免因单一敏感器故障导致的姿态信息紊乱,保证姿态的稳定性和可靠性。本发明的一种小型卫星多模冗余自主姿态确定方法,借助中低精度敏感器,通过数据融合处理解算得到高精度姿态信息,保证卫星业务运行的需求。
第二实施例
图13示出了本申请的另一个实施例提供的一种计算机设备的结构示意图。图13显示的计算机设备50仅仅是一个示例,不应对本申请实施例的功能和使用范围带来任何限制。如图13所示,计算机设备50以通用计算设备的形式表现。计算机设备50的组件可以包括但不限于:一个或者多个处理器或者处理单元500,系统存储器516,连接不同系统组件(包括系统存储器516和处理单元500)的总线501。
总线501表示几类总线结构中的一种或多种,包括存储器总线或者存储器控制器,外围总线,图形加速端口,处理器或者使用多种总线结构中的任意总线结构的局域总线。举例来说,这些体系结构包括但不限于工业标准体系结构(ISA)总线,微通道体系结构(MAC)总线,增强型ISA总线、视频电子标准协会(VESA)局域总线以及外围组件互连(PCI)总线。
计算机设50典型地包括多种计算机系统可读介质。这些介质可以是任何能够被计算机设备50访问的可用介质,包括易失性和非易失性介质,可移动的和不可移动的介质。
系统存储器516可以包括易失性存储器形式的计算机系统可读介质,例如随机存取存储器(RAM)504和/或高速缓存存储器506。计算机设备50可以进一步包括其它可移动/不可移动的、易失性/非易失性计算机系统存储介质。仅作为举例,存储系统508可以用于读写不可移动的、非易失性磁介质(图13未显示,通常称为“硬盘驱动器”)。尽管图13中未示出,可以提供用于对可移动非易失性磁盘(例如“软盘”)读写的磁盘驱动器,以及对可移动非易失性光盘(例如CD-ROM,DVD-ROM或者其它光介质)读写的光盘驱动器。在这些情况下,每个驱动器可以通过一个或者多个数据介质接口与总线501相连。存储器516可以包括至少一个程序产品,该程序产品具有一组(例如至少一个)程序模块,这些程序模块被配置以执行实施例一的功能。
具有一组(至少一个)程序模块512的程序/实用工具510,可以存储在例如存储器516中,这样的程序模块512包括但不限于操作系统、一个或者多个应用程序、其它程序模块以及程序数据,这些示例中的每一个或某种组合中可能包括网络环境的实现。程序模块512通常执行本申请所描述的实施例中的功能和/或方法。
计算机设备50也可以与一个或多个外部设备70(例如键盘、指向设备、显示器60等)通信,还可与一个或者多个使得用户能与该计算机设备50交互的设备通信,和/或与使得该计算机设备50能与一个或多个其它计算设备进行通信的任何设备(例如网卡,调制解调器等等)通信。这种通信可以通过输入/输出(I/O)接口502进行。并且,计算机设备50还可以通过网络适配器514与一个或者多个网络(例如局域网(LAN),广域网(WAN)和/或公共网络,例如因特网)通信。如图13所示,网络适配器514通过总线501与计算机设备50的其它模块通信。应当明白,尽管图13中未示出,可以结合计算机设备50使用其它硬件和/或软件模块,包括但不限于:微代码、设备驱动器、冗余处理单元、外部磁盘驱动阵列、RAID系统、磁带驱动器以及数据备份存储系统等。
处理器单元500通过运行存储在系统存储器516中的程序,从而执行各种功能应用以及数据处理,例如实现本申请实施例一所提供的一种卫星自主姿态确定方法。
本申请针对目前现有的问题,制定一种适用于卫星自主姿态确定方法的计算机设备,有效避免因单一敏感器故障导致的姿态信息紊乱,保证姿态的稳定性和可靠性。
第三实施例
本申请的另一个实施例提供了一种计算机可读存储介质,其上存储有计算机程序,该程序被处理器执行时实现上述实施例一所提供的方法。在实际应用中,所述计算机可读存储介质可以采用一个或多个计算机可读的介质的任意组合。计算机可读介质可以是计算机可读信号介质或者计算机可读存储介质。
计算机可读存储介质例如可以是但不限于电、磁、光、电磁、红外线、或半导体的系统、装置或器件,或者任意以上的组合。计算机可读存储介质的更具体的例子(非穷举的列表)包括:具有一个或多个导线的电连接、便携式计算机磁盘、硬盘、随机存取存储器(RAM)、只读存储器(ROM)、可擦式可编程只读存储器(EPROM或闪存)、光纤、便携式紧凑磁盘只读存储器(CD-ROM)、光存储器件、磁存储器件、或者上述的任意合适的组合。在本实施例中,计算机可读存储介质可以是任何包含或存储程序的有形介质,该程序可以被指令执行系统、装置或者器件使用或者与其结合使用。
计算机可读的信号介质可以包括在基带中或者作为载波一部分传播的数据信号,其中承载了计算机可读的程序代码。这种传播的数据信号可以采用多种形式,包括但不限于电磁信号、光信号或上述的任意合适的组合。计算机可读的信号介质还可以是计算机可读存储介质以外的任何计算机可读介质,该计算机可读介质可以发送、传播或者传输用于由指令执行系统、装置或者器件使用或者与其结合使用的程序。
计算机可读介质上包含的程序代码可以用任何适当的介质传输,包括但不限于无线、电线、光缆、RF等等,或者上述的任意合适的组合。可以以一种或多种程序设计语言或其组合来编写用于执行本申请操作的计算机程序代码,所述程序设计语言包括面向对象的程序设计语言-诸如Java、Smalltalk、C++,还包括常规的过程式程序设计语言-诸如“C”语言或类似的程序设计语言。程序代码可以完全地在用户计算机上执行、部分地在用户计算机上执行、作为一个独立的软件包执行、部分在用户计算机上部分在远程计算机上执行、或者完全在远程计算机或服务器上执行。在涉及远程计算机的情形中,远程计算机可以通过任意种类的网络——包括局域网(LAN)或广域网(WAN)-连接到用户计算机,或者,可以连接到外部计算机(例如利用因特网服务提供商来通过因特网连接)。
显然,本发明的上述实施例仅仅是为清楚地说明本发明所作的举例,而并非是对本发明的实施方式的限定,对于所属领域的普通技术人员来说,在上述说明的基础上还可以做出其它不同形式的变化或变动,这里无法对所有的实施方式予以穷举,凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。
Claims (9)
1.一种卫星自主姿态确定方法,其特征在于,包括:
S100、通过各敏感器测量信息进行敏感器姿态解算得到姿态信息,其中,各敏感器姿态解算包括陀螺递推姿态解算、双星敏姿态解算、单星敏姿态解算、太敏和地敏双矢量姿态解算、太敏和磁双矢量姿态解算以及地敏和磁双矢量解算;
S102、将步骤S100中经陀螺和其他敏感器姿态解算得到的姿态信息进行有效组合完成滤波姿态解算,进一步得到连续的姿态信息,其中,所述滤波姿态解算包括互补滤波姿态解算和扩展卡尔曼滤波姿态解算即EKF姿态解算;
S104、制定数据判断选用规则,考虑各敏感器及滤波定姿的精度,设计数据判断流程,在步骤S102的互补滤波姿态解算和EKF姿态解算的结果中选取最优化姿态信息;
所述互补滤波算法的算式为:
qerr=q-1×qt
q=q×q_tmp
式中,q为陀螺递推解算得到的姿态四元数;qt为星敏、太敏和地敏、太敏和磁、或地敏和磁解算得到的姿态四元数;qerr为两者姿态四元数的差值;k为互补滤波增益系数;q_tmp为qerr的归一化值。
2.根据权利要求1所述的方法,其特征在于,
在进行所述敏感器姿态解算时,除陀螺递推姿态解算外,依据优先级仅采用一种姿态解算方法,敏感器优先级依次为双星敏、单星敏、太敏和地敏、太敏和磁,或地敏和磁。
3.根据权利要求1所述的方法,其特征在于,
在进行所述EKF姿态解算时,分为四种EKF模式,包括:星敏+陀螺EKF、太敏和地敏+陀螺EKF、太敏和磁+陀螺EKF、地敏和磁+陀螺EKF,根据滤波有效性判断依据,选择最优滤波方案进行姿态解算到姿态信息。
4.根据权利要求3所述的方法,其特征在于,所述进行EKF姿态解算的步骤包括:
S1020、确定各敏感器测量模型;
S1021、确定各种EKF模式计算方法;
S1022、进行EKF有效性判断,根据预定规则判断四种EKF模式的有效性,选取最优项。
5.根据权利要求4所述的方法,其特征在于,
所述预定规则为:判断敏感器失效时长是否超限;若敏感器失效时长未超限,判断EKF是否收敛,如果在一定时间内系统估计均方差阵的范数稳定且小于设定阈值则EKF收敛;若EKF已收敛且星敏感器定姿有效,将EKF结果与星敏感器定姿结果进行比较,误差不超限时认为EKF有效,优先使用滤波结果。
6.根据权利要求1所述的方法,其特征在于,
根据所述数据判断选用规则,经过综合计算得到双星敏定姿精度最优,其次为星敏EKF定姿,之后依次为太敏和地敏EKF定姿、单星敏定姿、太敏和地敏定姿、太敏和磁EKF定姿、地敏和磁EKF定姿、太敏和磁定姿或地敏和磁定姿。
7.根据权利要求6所述的方法,其特征在于,
若判定选用EKF定姿,则将EKF解算得到的角速度信息、姿态信息赋值到系统相应信息,若采用敏感器定姿,则无需使用EKF解算得到的数据;当敏感器定姿无效且EKF尚未收敛时,仍采用EKF解算的姿态信息。
8.一种计算机可读存储介质,其上存储有计算机程序,其特征在于,该程序被处理器执行时实现如权利要求1-7中任一项所述的方法。
9.一种计算机设备,包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序,其特征在于,所述处理器执行所述程序时实现如权利要求1-7中任一项所述的方法。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110565164.9A CN113483765B (zh) | 2021-05-24 | 2021-05-24 | 一种卫星自主姿态确定方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110565164.9A CN113483765B (zh) | 2021-05-24 | 2021-05-24 | 一种卫星自主姿态确定方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN113483765A CN113483765A (zh) | 2021-10-08 |
CN113483765B true CN113483765B (zh) | 2023-03-24 |
Family
ID=77933027
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202110565164.9A Active CN113483765B (zh) | 2021-05-24 | 2021-05-24 | 一种卫星自主姿态确定方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN113483765B (zh) |
Families Citing this family (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113932802B (zh) * | 2021-10-12 | 2024-05-14 | 中国科学院微小卫星创新研究院 | 多个星敏感器的优先级变更方法及系统 |
CN114001602A (zh) * | 2021-10-26 | 2022-02-01 | 东北大学秦皇岛分校 | 基于四元数卡尔曼滤波去噪融合的火箭炮扰动检测方法 |
Family Cites Families (10)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US6463366B2 (en) * | 2000-03-10 | 2002-10-08 | Schafer Corp | Attitude determination and alignment using electro-optical sensors and global navigation satellites |
US8825399B2 (en) * | 2008-07-24 | 2014-09-02 | Raytheon Company | System and method of passive and autonomous navigation of space vehicles using an extended Kalman filter |
CN105158784A (zh) * | 2015-07-07 | 2015-12-16 | 中国人民解放军第二炮兵工程大学 | 动中通卫星通信系统级联卡尔曼滤波载体姿态估计方法 |
EP3182067A1 (en) * | 2015-12-18 | 2017-06-21 | Universite De Montpellier | Method and apparatus for determining spacecraft attitude by tracking stars |
CN107747953B (zh) * | 2017-10-25 | 2020-04-24 | 上海航天控制技术研究所 | 一种多敏感器数据与轨道信息时间同步方法 |
CN108279010A (zh) * | 2017-12-18 | 2018-07-13 | 北京时代民芯科技有限公司 | 一种基于多传感器的微小卫星姿态确定方法 |
CN108168549B (zh) * | 2018-03-24 | 2018-10-30 | 陕西博实测控科技有限公司 | 一种卫星动中通姿态检测方法 |
CN108827299B (zh) * | 2018-03-29 | 2022-04-12 | 南京航空航天大学 | 一种基于改进四元数二阶互补滤波的飞行器姿态解算方法 |
CN110411438B (zh) * | 2019-07-12 | 2021-02-09 | 北京控制工程研究所 | 一种基于多星敏感器的自适应组合确定卫星姿态角的方法 |
CN110608741A (zh) * | 2019-09-25 | 2019-12-24 | 北京安达维尔科技股份有限公司 | 一种提高飞机航姿参考系统姿态解算精度的方法 |
-
2021
- 2021-05-24 CN CN202110565164.9A patent/CN113483765B/zh active Active
Also Published As
Publication number | Publication date |
---|---|
CN113483765A (zh) | 2021-10-08 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN109001787B (zh) | 一种姿态角解算与定位的方法及其融合传感器 | |
CN102692225B (zh) | 一种用于低成本小型无人机的姿态航向参考系统 | |
Gebre-Egziabher et al. | A gyro-free quaternion-based attitude determination system suitable for implementation using low cost sensors | |
CN113483765B (zh) | 一种卫星自主姿态确定方法 | |
CN111102978A (zh) | 一种车辆运动状态确定的方法、装置及电子设备 | |
US20050240347A1 (en) | Method and apparatus for adaptive filter based attitude updating | |
CN110954102B (zh) | 用于机器人定位的磁力计辅助惯性导航系统及方法 | |
WO2004063669A2 (en) | Attitude change kalman filter measurement apparatus and method | |
CN106767767A (zh) | 一种微纳多模星敏感器系统及其数据融合方法 | |
Soken | A survey of calibration algorithms for small satellite magnetometers | |
CN110285815A (zh) | 一种可在轨全程应用的微纳卫星多源信息姿态确定方法 | |
CN114485641A (zh) | 一种基于惯导卫导方位融合的姿态解算方法及装置 | |
CN116817896A (zh) | 一种基于扩展卡尔曼滤波的姿态解算方法 | |
CN115033844A (zh) | 一种无人机状态估计方法、系统、设备及可读存储介质 | |
CN111141285B (zh) | 一种航空重力测量装置 | |
CN115727842B (zh) | 一种无人机快速对准方法、系统、计算机设备及存储介质 | |
Somov et al. | Guidance, navigation and control of a surveying satellite when an area imagery for disaster management. | |
CN115878939A (zh) | 基于飞行器舵面偏转的高精度动态测量方法 | |
CN111141283A (zh) | 一种通过地磁数据判断行进方向的方法 | |
Candan et al. | Estimation of attitude using robust adaptive Kalman filter | |
CN113932803B (zh) | 适用于高动态飞行器的惯性/地磁/卫星组合导航系统 | |
CN112882118B (zh) | 地固坐标系下动基座重力矢量估计方法、系统及存储介质 | |
JPH11248456A (ja) | 3軸姿勢検出装置 | |
CN114526729A (zh) | 一种基于冗余技术的mems惯性定位系统航向优化方法 | |
CN117508642B (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 |