CN103217544B - 用星敏感器星点位置变化估算星体角速度的方法及系统 - Google Patents

用星敏感器星点位置变化估算星体角速度的方法及系统 Download PDF

Info

Publication number
CN103217544B
CN103217544B CN201310091394.1A CN201310091394A CN103217544B CN 103217544 B CN103217544 B CN 103217544B CN 201310091394 A CN201310091394 A CN 201310091394A CN 103217544 B CN103217544 B CN 103217544B
Authority
CN
China
Prior art keywords
star
sigma
asterism
moment
alpha
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
CN201310091394.1A
Other languages
English (en)
Other versions
CN103217544A (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 Xinyue Instrument Factory
Original Assignee
Shanghai Xinyue Instrument Factory
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 Xinyue Instrument Factory filed Critical Shanghai Xinyue Instrument Factory
Priority to CN201310091394.1A priority Critical patent/CN103217544B/zh
Publication of CN103217544A publication Critical patent/CN103217544A/zh
Application granted granted Critical
Publication of CN103217544B publication Critical patent/CN103217544B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Landscapes

  • Navigation (AREA)

Abstract

本发明公开了一种用星敏感器星点位置变化估算星体角速度的方法,包括步骤:计算星图转角α;计算星图平移量(Δx,Δy);根据星图转角α和星图平移量(Δx,Δy)计算惯性角速度。本发明还提供相应的系统。本发明可以在敏感器完成星图识别前提供角速度信息,解决了以往算法在没有连续有效姿态四元数的条件下不能完成角速度输出功能,大幅度放宽了限制条件。角速度的算法比较简单,总体上计算量很小,对于星敏感器处理器速度,各项功能的正常完成影响很小。在控制系统实际使用方面,可以在飞行器发射入轨阶段,能够快速、并可靠地输出角速度信息,用于速率阻尼、姿态快速机动等各种控制模式,从提高系统可靠性等工程应用方面,有着很大实际意义。

Description

用星敏感器星点位置变化估算星体角速度的方法及系统
技术领域
本发明涉及一种利用星敏感器星点位置变化估算星体角速度的技术。
背景技术
星敏感器姿态测量精度高、产品性能可靠,在空间飞行器上的应用广泛。部分星敏感器产品根据前后拍的姿态的变化估算并输出惯性角速度,完成这种方法的条件是星敏感器需要运行在正常跟踪、并姿态确定数据有连续有效。
发明内容
针对现有技术存在的不足,本发明解决的技术问题是提出一种采用星敏感器中星点位置变化来估算角速度的方法,通过星点位置变化与角度变化之间的关系估算角速度,可消除星敏感器须有连续姿态确定数据的比较苛刻的条件,以改善产品使用适应性。
为解决上述技术问题,本发明是通过以下的技术方案实现的。
根据本发明的一个方面,提供一种用星敏感器星点位置变化估算星体角速度的方法,包括如下步骤:
步骤1:计算星图转角α;
步骤2:计算星图平移量(Δx,Δy);
步骤3:根据所述星图转角α和星图平移量(Δx,Δy)计算惯性角速度。
优选地,所述步骤1,具体为,
第i和第j个星点间相对位置矢量的计算公式为:
q → ij ( k ) = ( x i ( k ) - x j ( k ) , y i ( k ) - y j ( k ) )
式中:(xi(k),yi(k))为k时刻第i个星点的标定位置,i=1,2,...,n,yi(k)-yj(k))为k时刻第j个星点的标定位置,j=1,2,...,n;
k时刻的星点相对位置矢量和k+1时刻的星点相对位置矢量之间的夹角αij(k+1)的计算公式如下:
α ij ( k + 1 ) = arccos ( q → ij ( k ) q → ij ( k + 1 ) | q → ij ( k ) | | q → ij ( k + 1 ) | )
得到个角度值,对这个角度值进行平均处理,得到第k时刻到k+1时刻星图的转动角度α(k+1)的计算公式如下:
α ( k + 1 ) = 1 C n 2 Σ i = 1 n - 1 Σ j = i + 1 n α ij ( k + 1 ) .
优选地,所述步骤2,具体为,
第k时刻各星点在空间中的位置(xi′(k),yi′(k))计算公式如下:
xi′(k)=xi(k)·cosα(k+1)-yi(k)·sinα(k+1)
yi′(k)=xi(k)·sinα(k+1)+yi(k)·cosα(k+1)
星点位置平移量计算公式如下:
Δxi(k+1)=xi(k+1)-xi′(k)
Δyi(k+1)=yi(k+1)-yi′(k)
对n个平移值求平均,得到k+1时刻星图的位置平移为:
Δx ( k + 1 ) = 1 n ( Σ i = 1 n x i ( k + 1 ) - Σ i = 1 n x i ( k ) · cos α ( k + 1 ) + Σ i = 1 n y i ( k ) · sin α ( k + 1 ) )
Δy ( k + 1 ) = 1 n ( Σ i = 1 n y i ( k + 1 ) - Σ i = 1 n x i ( k ) · sin α ( k + 1 ) - Σ i = 1 n y i ( k ) · cos α ( k + 1 ) ) .
优选地,所述步骤3,具体为,
惯性角速度wi(k+1)计算公式如下:
w i ( k + 1 ) = 1 Tb ( Δy ( k + 1 ) , - Δx ( k + 1 ) , - α ( k + 1 ) ) T
其中,Tb为第k到k+1之间的时间间隔。
根据本发明的另一个方面,还提供一种用星敏感器星点位置变化估算星体角速度的系统,包括如下装置:
第一计算装置,用于计算星图转角α;
第二计算装置,用于计算星图平移量(Δx,Δy);
第三计算装置,用于根据所述星图转角α和星图平移量(Δx,Δy)计算惯性角速度。
优选地,所述第一计算装置,具体为,
第i和第j个星点间相对位置矢量的计算公式为:
q → ij ( k ) = ( x i ( k ) - x j ( k ) , y i ( k ) - y j ( k ) )
式中:(xi(k),yi(k))为k时刻第i个星点的标定位置,i=1,2,...,n,yi(k)-yj(k))为k时刻第j个星点的标定位置,j=1,2,...,n;
k时刻的星点相对位置矢量和k+1时刻的星点相对位置矢量之间的夹角αij(k+1)的计算公式如下:
α ij ( k + 1 ) = arccos ( q → ij ( k ) q → ij ( k + 1 ) | q → ij ( k ) | | q → ij ( k + 1 ) | )
得到个角度值,对这个角度值进行平均处理,得到第k时刻到k+1时刻星图的转动角度α(k+1)的计算公式如下:
α ( k + 1 ) = 1 C n 2 Σ i = 1 n - 1 Σ j = i + 1 n α ij ( k + 1 ) .
优选地,所述第二计算装置,具体为,
第k时刻各星点在空间中的位置(xi′(k),yi′(k))计算公式如下:
xi′(k)=xi(k)·cosα(k+1)-yi(k)·sinα(k+1)
yi′(k)=xi(k)·sinα(k+1)+yi(k)·cosα(k+1)
星点位置平移量计算公式如下:
Δxi(k+1)=xi(k+1)-xi′(k)
Δyi(k+1)=yi(k+1)-yi′(k)
对n个平移值求平均,得到k+1时刻星图的位置平移为:
Δx ( k + 1 ) = 1 n ( Σ i = 1 n x i ( k + 1 ) - Σ i = 1 n x i ( k ) · cos α ( k + 1 ) + Σ i = 1 n y i ( k ) · sin α ( k + 1 ) )
Δy ( k + 1 ) = 1 n ( Σ i = 1 n y i ( k + 1 ) - Σ i = 1 n x i ( k ) · sin α ( k + 1 ) - Σ i = 1 n y i ( k ) · cos α ( k + 1 ) ) .
优选地,所述第三计算装置,具体为,
惯性角速度wi(k+1)计算公式如下:
w i ( k + 1 ) = 1 Tb ( Δy ( k + 1 ) , - Δx ( k + 1 ) , - α ( k + 1 ) ) T
其中,Tb为第k到k+1之间的时间间隔。
星敏感器采用星点位移算法计算惯性角速度的方法,可以在敏感器完成星图识别前提供角速度信息,解决了以往算法在没有连续有效姿态四元数的条件下不能完成角速度输出功能,大幅度放宽了限制条件。
角速度的算法比较简单,比通常的角速度算法计算量略为增加,但总体上计算量还是很小,对于星敏感器处理器速度,各项功能的正常完成影响很小。
在控制系统实际使用方面,可以在飞行器发射入轨阶段,能够快速、并可靠地输出角速度信息,用于速率阻尼、姿态快速机动等各种控制模式,对提高系统可靠性等工程应用方面,有着很大实际意义。
附图说明
通过阅读参照以下附图对非限制性实施例所作的详细描述,本发明的其它特征、目的和优点将会变得更明显:
图1为星点位置变化示意图。
具体实施方式
下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明,但不以任何形式限制本发明。应当指出的是,对本领域的普通技术人员来说,在不脱离本发明构思的前提下,还可以做出若干变形和改进。这些都属于本发明的保护范围。
1、星图转角α的计算
星敏感器CCD上像点位置变化,是由于敏感器坐标系绕某个空间轴转动引起。根据星点相对位置关系保持固定,敏感器坐标系绕某空间轴旋转一个角度,敏感器CCD测得的星点图形形状不变。星点标定位置变化可以分解为绕Zs轴转角α,与平面内的位置平动(Δx,Δy)二部分,如图1所示,设某星点开始位置(x1,y1),绕Os转动质心位置到(x′,y′),再通过平移到(x2,y2)。
根据敏感器坐标系运动关系,星图平移量Δx、Δy和转动角度α,是以3、1、2转序下得到的欧拉转角,分别对应于俯仰角、滚动角和偏航角。
图1中可以看出,线段A1B1绕Os转动α得到A′B′,因此只要计算线段A1B1与A′B′间的夹角就可得到转角α。点A1和B1为移动前的位置,敏感器可以测量得到,而点A′和B′只是星点位置移动前后理论分解的坐标,不能直接测量得到。
由于A′B′通过平移得到A2B2,A′B′与A2B2平行,因此A1B1与A′B′间的夹角,等于A1B1与A2B2之间的夹角,而点A2、B2位置可以测量得到,因此通过计算星点相对位置矢量A1B1和A2B2的夹角得到转动角度。
通常情况下敏感器中有多个恒星点,当经过处理判断敏感器视场中有效的星点数目为n,k时刻测量得到的标定位置分别为(xi(k),yi(k)),i=1,2,...,n。通过矢量计算,第i和第j个星点间相对位置矢量为
q → ij ( k ) = ( x i ( k ) - x j ( k ) , y i ( k ) - y j ( k ) ) - - - ( 1 )
其中,i=1,2,...,n-1,j=i+1,2,...,n。通过(1)式计算k和k+1时刻的星点相对位置矢量两个矢量之间的夹角可以用余弦定理计算得到:
α ij ( k + 1 ) = arccos ( q → ij ( k ) q → ij ( k + 1 ) | q → ij ( k ) | | q → ij ( k + 1 ) | ) - - - ( 2 )
可以得到个角度值,对这些值的进行平均处理以提高估算精度,得到第k时刻到k+1时刻星图的转动角度α(k+1):
α ( k + 1 ) = 1 C n 2 Σ i = 1 n - 1 Σ j = i + 1 n α ij ( k + 1 ) - - - ( 3 )
2、星图平移量(Δx,Δy)的计算
第k时刻的n个星点标定位置(xi(k),yi(k)),i=1,2,...,n,因星图绕光轴Os转动角度为α(k+1),由此各星点的标定位置同样绕Os转动α(k+1)角度,得到各星点在空间中的位置(xi′(k),yi′(k)):
xi′(k)=xi(k)·cosα(k+1)-yi(k)·sinα(k+1)    (4)
yi′(k)=xi(k)·sinα(k+1)+yi(k)·cosα(k+1)    (5)
敏感器在k+1时刻测得的星点标定位置为(xi(k+1),yi(k+1)),与第k时刻转动α(k+1)得到各星点在空间中的位置比较,可以得到n个星点位置平移量:
Δxi(k+1)=xi(k+1)-xi′(k)    (6)
Δyi(k+1)=yi(k+1)-yi′(k)    (7)
对n个平移值求平均,得到k+1时刻星图的位置平移为:
Δx ( k + 1 ) = 1 n ( Σ i = 1 n x i ( k + 1 ) - Σ i = 1 n x i ( k ) · cos α ( k + 1 ) + Σ i = 1 n y i ( k ) · sin α ( k + 1 ) ) - - - ( 8 )
Δy ( k + 1 ) = 1 n ( Σ i = 1 n y i ( k + 1 ) - Σ i = 1 n x i ( k ) · sin α ( k + 1 ) - Σ i = 1 n y i ( k ) · cos α ( k + 1 ) ) - - - ( 9 )
3、惯性角速度的估算
根据转动运动,敏感器绕OsXs轴正向转动,像点的纵坐标将变大,Δy为正值;敏感器绕OsYs正向转动,像点的横坐标将变小,Δx为负值;敏感器绕OsZs正向转动,星图将反相转动,α为负值。
若星敏感器曝光时刻第k到k+1之间的时间间隔为Tb,敏感器中的星图位置变化分解为星图位置绕Zs轴转角α(k+1)和平移(Δx(k+1),Δy(k+1))。
因此相对于敏感器坐标系,飞行器的惯性角速度wi(k+1)可以用下式进行测算:
w i ( k + 1 ) = 1 Tb ( Δy ( k + 1 ) , - Δx ( k + 1 ) , - α ( k + 1 ) ) T - - - ( 10 )
以上对本发明的具体实施例进行了描述。需要理解的是,本发明并不局限于上述特定实施方式,本领域技术人员可以在权利要求的范围内做出各种变形或修改,这并不影响本发明的实质内容。

Claims (2)

1.一种用星敏感器星点位置变化估算星体角速度的方法,其特征在于,包括如下步骤:
步骤1:计算星图转角α;
步骤2:计算星图平移量(Δx,Δy);
步骤3:根据所述星图转角α和星图平移量(Δx,Δy)计算惯性角速度;
所述步骤1,具体为,
第i和第j个星点间相对位置矢量的计算公式为:
式中:(xi(k),yi(k))为k时刻第i个星点的标定位置,i=1,2,...,n,(xj(k),yj(k))为k时刻第j个星点的标定位置,j=1,2,...,n;
k时刻的星点相对位置矢量和k+1时刻的星点相对位置矢量之间的夹角αij(k+1)的计算公式如下:
得到个角度值,对这个角度值进行平均处理,得到第k时刻到k+1时刻星图的转动角度α(k+1)的计算公式如下:
α ( k + 1 ) = 1 C n 2 Σ i = 1 n - 1 Σ j = i + 1 n α ij ( k + 1 ) ;
所述步骤2,具体为,
第k时刻各星点在空间中的位置(x′i(k),y′i(k))计算公式如下:
x′z(k)=xi(k)·cosα(k+1)-yi(k)·sinα(k+1)
y′i(k)=xi(k)·sinα(k+1)+yi(k)·cosα(k+1)
星点位置平移量计算公式如下:
Δxi(k+1)=xi(k+1)-x′i(k)
Δyi(k+1)=yi(k+1)-y′i(k)
对n个平移值求平均,得到k+1时刻星图的位置平移为:
Δx ( k + 1 ) = 1 n ( Σ i = 1 n x i ( k + 1 ) - Σ i = 1 n x i ( k ) · cos α ( k + 1 ) + Σ i = 1 n y i ( k ) · sin α ( k + 1 ) )
Δy ( k + 1 ) = 1 n ( Σ i = 1 n y i ( k + 1 ) - Σ i = 1 n x i ( k ) · sin α ( k + 1 ) - Σ i = 1 n y i ( k ) · cos α ( k + 1 ) ) ;
所述步骤3,具体为,
惯性角速度wi(k+1)计算公式如下:
w i ( k + 1 ) = 1 Tb ( Δy ( k + 1 ) , - Δx ( k + 1 ) , - α ( k + 1 ) ) T
其中,Tb为第k到k+1之间的时间间隔。
2.一种用星敏感器星点位置变化估算星体角速度的系统,其特征在于,包括如下装置:
第一计算装置,用于计算星图转角α;
第二计算装置,用于计算星图平移量(Δx,Δy);
第三计算装置,用于根据所述星图转角α和星图平移量(Δx,Δy)计算惯性角速度;
所述第一计算装置,具体为,
第i和第j个星点间相对位置矢量的计算公式为:
式中:(xi(k),yi(k))为k时刻第i个星点的标定位置,i=1,2,...,n,为k时刻第j个星点的标定位置,j=1,2,...,n;
k时刻的星点相对位置矢量和k+1时刻的星点相对位置矢量之间的夹角αij(k+1)的计算公式如下:
得到个角度值,对这个角度值进行平均处理,得到第k时刻到k+1时刻星图的转动角度α(k+1)的计算公式如下:
α ( k + 1 ) = 1 C n 2 Σ i = 1 n - 1 Σ j = i + 1 n α ij ( k + 1 ) ;
所述第二计算装置,具体为,
第k时刻各星点在空间中的位置(x′i(k),y′i(k))计算公式如下:
x′i(k)=xi(k)·cosα(k+1)-yi(k)·sinα(k+1)
y′i(k)=xi(k)·sinα(k+1)+yi(k)·cosα(k+1)
星点位置平移量计算公式如下:
Δxi(k+1)=xi(k+1)-x′i(k)
Δyi(k+1)=yi(k+1)-y′i(k)
对n个平移值求平均,得到k+1时刻星图的位置平移为:
Δx ( k + 1 ) = 1 n ( Σ i = 1 n x i ( k + 1 ) - Σ i = 1 n x i ( k ) · cos α ( k + 1 ) + Σ i = 1 n y i ( k ) · sin α ( k + 1 ) )
Δy ( k + 1 ) = 1 n ( Σ i = 1 n y i ( k + 1 ) - Σ i = 1 n x i ( k ) · sin α ( k + 1 ) - Σ i = 1 n y i ( k ) · cos α ( k + 1 ) ) ;
所述第三计算装置,具体为,
惯性角速度wi(k+1)计算公式如下:
w i ( k + 1 ) = 1 Tb ( Δy ( k + 1 ) , - Δx ( k + 1 ) , - α ( k + 1 ) ) T
其中,Tb为第k到k+1之间的时间间隔。
CN201310091394.1A 2013-03-21 2013-03-21 用星敏感器星点位置变化估算星体角速度的方法及系统 Active CN103217544B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201310091394.1A CN103217544B (zh) 2013-03-21 2013-03-21 用星敏感器星点位置变化估算星体角速度的方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201310091394.1A CN103217544B (zh) 2013-03-21 2013-03-21 用星敏感器星点位置变化估算星体角速度的方法及系统

Publications (2)

Publication Number Publication Date
CN103217544A CN103217544A (zh) 2013-07-24
CN103217544B true CN103217544B (zh) 2015-06-03

Family

ID=48815510

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201310091394.1A Active CN103217544B (zh) 2013-03-21 2013-03-21 用星敏感器星点位置变化估算星体角速度的方法及系统

Country Status (1)

Country Link
CN (1) CN103217544B (zh)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105547327B (zh) * 2015-12-10 2018-10-30 上海新跃仪表厂 一种基于空间转换的星敏感器精度测试方法
CN106441282B (zh) * 2016-09-19 2019-07-09 上海航天控制技术研究所 一种星敏感器星跟踪方法
CN107588768B (zh) * 2017-08-21 2020-07-07 中国科学院长春光学精密机械与物理研究所 基于星图的帧间角速度计算方法
CN107607737B (zh) * 2017-08-21 2019-07-05 中国科学院长春光学精密机械与物理研究所 基于星光矢量的无陀螺角速度测量方法
CN107942090B (zh) * 2017-12-28 2019-10-29 北京航空航天大学 一种基于模糊星图提取光流信息的航天器角速度估计方法
CN108508905B (zh) * 2018-04-08 2020-12-25 上海微小卫星工程中心 一种基于最短空间轴的姿态机动控制和导引律规划方法
CN113532445B (zh) * 2021-07-20 2022-07-05 北京控制工程研究所 一种卷帘曝光星敏感器的高动态快速自主捕获方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1229307A1 (de) * 2001-01-31 2002-08-07 Philips Corporate Intellectual Property GmbH Anordnung zum Erfassen des Drehwinkels eines drehbaren Elements
CN102435763A (zh) * 2011-09-16 2012-05-02 中国人民解放军国防科学技术大学 一种基于星敏感器的航天器姿态角速度测量方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1229307A1 (de) * 2001-01-31 2002-08-07 Philips Corporate Intellectual Property GmbH Anordnung zum Erfassen des Drehwinkels eines drehbaren Elements
CN102435763A (zh) * 2011-09-16 2012-05-02 中国人民解放军国防科学技术大学 一种基于星敏感器的航天器姿态角速度测量方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
星敏感器在轨标定算法研究;袁彦红;《中国优秀硕士学位论文全文数据库 工程科技Ⅱ辑》;20090215(第2(2009年)期);第36页至第48页 *

Also Published As

Publication number Publication date
CN103217544A (zh) 2013-07-24

Similar Documents

Publication Publication Date Title
CN103217544B (zh) 用星敏感器星点位置变化估算星体角速度的方法及系统
CN112649016B (zh) 一种基于点线初始化的视觉惯性里程计方法
Phuong et al. A DCM based orientation estimation algorithm with an inertial measurement unit and a magnetic compass
CN102435763B (zh) 一种基于星敏感器的航天器姿态角速度测量方法
CN109625335B (zh) 一种基于角速度估计信息和太阳敏感器的捕获太阳方法
CN106709222B (zh) 基于单目视觉的imu漂移补偿方法
CN108534784B (zh) 一种基于空间圆轨迹的非合作航天器自旋角速度估计方法
Shijie et al. Monocular vision-based two-stage iterative algorithm for relative position and attitude estimation of docking spacecraft
CN103268067B (zh) 一种基于拟四元数与拟四元数运动学方程的卫星指向跟踪控制方法
CN111609868A (zh) 一种基于改进光流法的视觉惯性里程计方法
CN102607558A (zh) 一种基于惯性测量单元的相对变形测量方法
CN108426584B (zh) 汽车多传感器的校准方法
CN108645400B (zh) 用于空间非合作目标相对导航的惯性参数辨识方法及系统
EP2879011B1 (en) On-board estimation of the nadir attitude of an Earth orbiting spacecraft
CN115540854A (zh) 一种基于uwb辅助的主动定位方法、设备和介质
Dhahbane et al. Attitude determination and attitude estimation in aircraft and spacecraft navigation. A survey
CN103466103B (zh) 一种太阳敏感器故障情况下的对日定向姿态控制方法
CN111649738B (zh) 微重力场下的加速度计初始姿态解算方法
US20220205788A1 (en) Ground vehicle monocular visual-inertial odometry via locally flat constraints
CN114964259A (zh) 基于红外焦平面阵列时域信号与红外图像的姿态解算方法及系统
CN117015719A (zh) 用于确定刚体的运动状态的方法
CN107168350B (zh) 一种针对服务航天器定轴转动时的最佳自转角速度的计算方法
Guo et al. Plane based visual odometry for structural and low-texture environments using rgb-d sensors
Yu et al. Stereo vision based relative state estimation for non-cooperative spacecraft with outliers
CN116909305B (zh) 基于时域变换的预设时间预设精度姿轨一体化跟踪控制方法及其控制系统

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant