CN108872973B - 一种弹道导弹目标定轨的ekf滤波方法 - Google Patents

一种弹道导弹目标定轨的ekf滤波方法 Download PDF

Info

Publication number
CN108872973B
CN108872973B CN201811000966.XA CN201811000966A CN108872973B CN 108872973 B CN108872973 B CN 108872973B CN 201811000966 A CN201811000966 A CN 201811000966A CN 108872973 B CN108872973 B CN 108872973B
Authority
CN
China
Prior art keywords
value
moment
current time
state
current
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
CN201811000966.XA
Other languages
English (en)
Other versions
CN108872973A (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.)
Beijing Institute of Electronic System Engineering
Original Assignee
Beijing Institute of Electronic System Engineering
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 Beijing Institute of Electronic System Engineering filed Critical Beijing Institute of Electronic System Engineering
Priority to CN201811000966.XA priority Critical patent/CN108872973B/zh
Publication of CN108872973A publication Critical patent/CN108872973A/zh
Application granted granted Critical
Publication of CN108872973B publication Critical patent/CN108872973B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/006—Theoretical aspects

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Radar Systems Or Details Thereof (AREA)
  • Aiming, Guidance, Guns With A Light Source, Armor, Camouflage, And Targets (AREA)

Abstract

一种弹道导弹目标定轨的EKF滤波方法，其特征在于，方法包括：确定EKF算法的初始值；计算当前时刻的目标定轨的系统误差值；建立EKF算法的系统模型，并在初始值与系统误差值的基础上进行滤波，以计算当前时刻的状态估计值；以及得到当前时刻的目标定轨评估范围。

Description

一种弹道导弹目标定轨的EKF滤波方法
技术领域
本发明涉及目标定轨领域。更具体地，涉及一种弹道导弹目标定轨的EKF(Extended Kalman Filter，扩展卡尔曼滤波器)滤波方法。
背景技术
目前适用于弹道导弹目标定轨滤波的算法很多，常用的滤波方法有EKF，QUKF、HIF、ESO和SVF。EKF、QUKF、HIF和SVF均需要系统的模型，ESO则不需要系统的精确模型。目前，EKF、QUKF、HIF和SVF的滤波精度比较接近，ESO的精度稍低一些。
以上方法均能实现对弹道导弹目标进行定轨滤波，但使用范围具有局限性(使用时不能够同时满足初始值不稳定、数据断片等情况)，且现有技术的这些方法仅考虑了随机性定轨误差而忽略系统性定轨误差，而在实际应用中系统性定轨误差常常是大于随机性定轨误差的。此外，现有的方法均不能对定轨滤波结果进行实时在线估计。
因此，需要提供一种能够保证在各种条件下均可实现算法收敛性和稳定性、考虑探测设备的系统误差与随机误差的定轨误差实时在线估计算法。
发明内容
本发明的目的在于提供一种能够对定轨结果进行实时在线估计的弹道导弹目标定轨的EKF滤波方法。
为达到上述目的，本发明采用下述技术方案：
一种弹道导弹目标定轨的EKF滤波方法，方法包括：确定EKF算法的初始值；计算当前时刻的目标定轨的系统误差值；建立所述EKF算法的系统模型，并在所述初始值与所述系统误差值的基础上进行滤波，以计算当前时刻的状态估计值；以及得到当前时刻的目标定轨误差评估范围。
优选地，所述计算当前时刻的状态估计值进一步包括：利用上一时刻的状态预报值或状态估计值、上一时刻的P矩阵和当前时刻的量测值计算当前时刻的状态估计值。
优选地，所述计算当前时刻的状态估计值进一步包括：当未能获得上一时刻的状态估计值时，利用所述上一时刻的状态预报值、所述上一时刻的P矩阵和所述当前时刻的量测值计算当前时刻的状态估计值；以及当能够获得上一时刻的状态估计值时，利用所述上一时刻的状态估计值、所述上一时刻的P矩阵和所述当前时刻的量测值计算当前时刻的状态估计值。
优选地，计算所述当前时刻的状态估计值的方程为
Figure BDA0001782951070000021
其中，k为上一时刻，k+1为当前时刻，
Figure BDA0001782951070000022
表示对当前时刻状态的预报值，
Figure BDA0001782951070000023
表示对上一时刻状态的估计值，
Figure BDA0001782951070000024
Figure BDA0001782951070000025
为上一时刻状态转移矩阵的预报值，Bs为增益矩阵，Kk+1为当前时刻的新息增益矩阵，C＝[I 0]，I为单位矩阵，Pk+1为当前时刻的协方差矩阵，
Figure BDA0001782951070000026
表示雷达当前时刻量测在地球系下的方差阵，
Figure BDA0001782951070000027
为当前时刻量测在地球系下的投影，Qk为上一时刻系统误差值的中间参量。
优选地，按照下述步骤计算当前时刻的系统误差值：设当前时刻的状态估计值为
Figure BDA0001782951070000028
雷达系到地球系的转换矩阵为
Figure BDA0001782951070000029
雷达系下的量测值分量为r,b,e，计算得到雷达系下的量测值分量r,b,e的导数估计值为：
Figure BDA0001782951070000031
其中
Figure BDA0001782951070000032
Figure BDA0001782951070000033
表示上一时刻位置估值在地心系的矢量，
Figure BDA0001782951070000034
表示上一时刻速度估值在地心系的矢量，rooR为雷达点位在地心系的值；
利用
Figure BDA0001782951070000035
进一步得到变量的估计值如下：
Figure BDA0001782951070000036
以及
计算得到系统性定轨误差值：
Figure BDA0001782951070000037
其中，Δr为测距的系统误差指标，Δb为方位角的系统误差指标，Δe为高低角的系统误差指标。
本发明的有益效果如下：
本发明所述技术方案提供了一种能够保证在各种条件下均可实现算法收敛性和稳定性、考虑探测设备的系统误差与随机误差的定轨误差实时在线估计算法。
附图说明
下面结合附图对本发明的具体实施方式作进一步详细的说明；
图1为示出根据本申请的示例性EKF滤波方法的定轨计算流程图；
图2为示出根据本申请的示例性EKF滤波方法的滤波初始化函数计算流程图；以及
图3为示出根据本申请的EKF滤波方法的定轨计算性能验证图。
具体实施方式
为了更清楚地说明本发明，下面结合优选实施例和附图对本发明做进一步的说明。附图中相似的部件以相同的附图标记进行表示。本领域技术人员应当理解，下面所具体描述的内容是说明性的而非限制性的，不应以此限制本发明的保护范围。
下面结合图1和图2描述本申请的方法。图1为示出根据本申请的示例性EKF滤波方法的定轨计算流程图；图2为示出根据本申请的示例性EKF滤波方法的滤波初始化函数计算流程图。
根据本申请的实施例，当对例如来袭的弹道导弹进行拦截时，需要根据探测设备探测到的量测数据预测被拦截的目标导弹的轨道，从而实施准确拦截。在本申请的实施例中，探测设备可以是雷达，其既可以是单雷达探测也可以是多雷达融合探测的情形。即在本申请中，既可以进行单雷达定轨滤波也可以进行多雷达融合定轨滤波，通过由雷达得到量测值利用本申请的滤波方法进行实时估计得准确的定轨误差评估范围，从而准确得目标导弹的轨道的范围估计。以下参照图2的流程图并结合图1详细描述本申请的示例性EKF滤波方法的定轨滤波方法。
如图1所示，在本申请的实施例中，如步骤S101，首先对计算EKF算法的初始值并对系统进行建模。
对于初始值的计算，在图2中以流程图的方式示出了计算过程。
具体地，在步骤S1011中，得到初始的雷达测量值r0,b0,e0即初始测量分量，该初始量测值是在雷达系下。
接下来在步骤S1013中，计算
Figure BDA0001782951070000041
和从雷达系到地球系的转换矩阵
Figure BDA0001782951070000042
其中图1中还示出了需要计算的雷达点位在地心系的值rooR以及地心到雷达系转换矩阵
Figure BDA0001782951070000043
以用于以下步骤的中间计算过程。在步骤S1015和S1017中计算并利用得到的雷达量测值分量r,b,e和
Figure BDA0001782951070000044
计算EKF滤波方法的初始值：
Figure BDA0001782951070000045
Figure BDA0001782951070000051
其中，
Figure BDA0001782951070000052
为初始状态估计值，
Figure BDA0001782951070000053
为雷达初始量测在地球系下的投影，R0为雷达初始量测在地球系下的方差阵，P0为初始矩阵。
接下来回到图1，在步骤S103中，计算当前时刻的系统误差值。具体地，设当前获得的状态估计值为
Figure BDA0001782951070000054
利用雷达系下的量测值分量r,b,e，计算它们的导数估计值
Figure BDA0001782951070000055
为：
Figure BDA0001782951070000056
应注意，多得到的分量的导数估计值是地球系下的，且其中，
Figure BDA0001782951070000057
Figure BDA0001782951070000058
表示上一时刻位置估值在地心系的矢量，
Figure BDA0001782951070000059
表示上一时刻速度估值在地心系的矢量。
利用
Figure BDA00017829510700000510
进一步得到中间变量的估计值如下：
Figure BDA00017829510700000511
最后，系统误差值计算如下：
Figure BDA00017829510700000512
其中，Δr为测距的系统误差指标，Δb为方位角的系统误差指标，Δe为高低角的系统误差指标。
下面进入步骤S105，计算当前时刻的状态估计值。值得注意的是，在该步中需要利用模型弹道导弹目标运动学方程对系统进行EKF建模得到系统的离散化模型。
首先根据模型弹道导弹目标运动学方程建立如下模型：
Figure BDA0001782951070000061
进而得到表达式：
Figure BDA0001782951070000062
其中，
Figure BDA0001782951070000063
利用式(7)对进行离散化处理得到如下离散化模型：
Figure BDA0001782951070000064
在步骤S105中，建模后，基于已经得到的当前时刻的系统误差值表达式和初始值，利用EKF方法进行滤波处理，得到当前时刻的状态估计值。
具体地，在本申请的实施例中，在进行滤波处理计算得到当前时刻的状态估计值时，需要用到上一时刻的状态预报值或状态估计值、上一时刻的P矩阵和当前时刻的量测值。而当前时刻的量测值对于当前时刻的状态估计值是必须的。
具体地，在本申请的实施例中，充分考虑到雷达探测时的各种情况。比如，在雷达探测过程中，在上一时刻，可能由于外界干扰或导弹位置探测不易而出现量测数据的断片。这样的情况下，对于上一时刻而言是不能够得到量测值的，无法计算状态估计值。因此，计算状态预报值，在当前时刻利用状态预报值充当上一时刻的状态估计值进行计算。当上一时刻没有存在数据断片，即能够得到上一时刻的量测值的情况下，利用上一时刻的状态估计值来计算当前时刻的状态估计值。
具体地，下面详细描述通过滤波计算得到当前时刻的状态估计值的过程，注意到，为了表述简便，当上一时刻量测数据断片时，以上一时刻的状态预报值充当上一时刻的状态估计值进行计算，在下面则没有专门在计算式中列出，仅以
Figure BDA0001782951070000071
表示上一时刻的状态估计值或上一时刻的状态预报值。
基于式(8)得到的离散化模型，利用EKF的上一时刻(k时刻)的状态估计值或状态预报值
Figure BDA0001782951070000072
上一时刻的P阵和当前时刻(k+1)时刻的量测值，生成当前时刻的状态估计值和当前时刻的P阵。
设上一时刻的状态估计值或上一时刻的状态预报值为
Figure BDA0001782951070000073
P阵为Pk，当前时刻的的量测值为Yk+1，则
Figure BDA0001782951070000074
的Pk+1计算如下：
进行一致EKF滤波算法：
Figure BDA0001782951070000075
其中，
Figure BDA0001782951070000076
表示当前时刻的状态预报值，
Figure BDA0001782951070000077
表示所述上一时刻的状态估计值或上一时刻的状态预报值，
Figure BDA0001782951070000078
Figure BDA0001782951070000079
为为上一时刻状态转移矩阵的预报值，Bs为增益矩阵，Kk+1为当前时刻的新息增益矩阵，C＝[I 0]，I为单位矩阵，Pk+1为当前时刻的协方差矩阵，
Figure BDA00017829510700000711
表示雷达当前时刻量测在地球系下的方差阵，
Figure BDA00017829510700000710
为当前时刻量测在地球系下的投影，Qk为上一时刻的系统误差值的中间参量，Qk通过如下算法得到：
Figure BDA0001782951070000081
其中a＞0是需要调节的参数。
此外，应注意，如果当前时刻雷达未能探测到量测值，则可以通过式(9)中的前两个等式来计算当前时刻的状态预报值，将该时刻的状态预报值作为下一时刻的基础计算下一时刻的状态估计值。即：
Figure BDA0001782951070000082
经过上述计算，在步骤S107中，可以得到当前时刻的目标定轨误差评估范围：
Figure BDA0001782951070000083
得到的定轨误差评估范围形成了当前时刻目标弹道导弹的定轨带，从而精确预测了当前时刻目标弹道导弹的位置所处的范围，以进行准确拦截。
本领域技术人员应理解，利用上述方法得到的状态估计值可以逐步递推，从而预测目标弹道导弹的整个轨迹，进而根据估计进行实行拦截。
下面参照图3对本申请的滤波方法进行验证。图3为示出根据本申请的EKF滤波方法的定轨计算性能验证图。
图3中，实线为已知目标实际弹道情况下计算的真实定轨误差的模值，两条虚线为不知该弹道情况下利用本申请的方法得到的实时估计范围。可以看出，本申请的方法能够实现对定轨精度的实时估计，效果良好。
显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定，对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动，这里无法对所有的实施方式予以穷举，凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。

Claims (4)

1.一种弹道导弹目标定轨的EKF滤波方法，其特征在于，所述方法包括：
确定EKF算法的初始值；
计算当前时刻的目标定轨的系统误差值；
建立所述EKF算法的系统模型，并在所述初始值与所述系统误差值的基础上进行滤波，以计算当前时刻的状态估计值；以及
得到当前时刻的目标定轨误差评估范围，
其中，所述计算当前时刻的系统误差值进一步包括：
设当前时刻的状态估计值为
Figure FDA0003590867600000011
雷达系到地球系的转换矩阵为
Figure FDA0003590867600000012
雷达系下的量测值分量为r,b,e，计算得到所述雷达系下的量测值分量r,b,e的导数估计值为：
Figure FDA0003590867600000013
其中
Figure FDA0003590867600000014
Figure FDA0003590867600000015
表示上一时刻位置估值在地心系的矢量，
Figure FDA0003590867600000016
表示上一时刻速度估值在地心系的矢量，rooR为雷达点位在地心系的值；
利用
Figure FDA0003590867600000017
进一步得到变量的估计值如下：
Figure FDA0003590867600000021
以及
计算得到所述当前时刻的系统性误差值：
Figure FDA0003590867600000022
其中，Δr为测距的系统误差指标，Δb为方位角的系统误差指标，Δe为高低角的系统误差指标。
2.如权利要求1所述的弹道导弹目标定轨的EKF滤波方法，其特征在于，所述计算当前时刻的状态估计值进一步包括：
利用上一时刻的状态预报值或状态估计值、上一时刻的P矩阵和当前时刻的量测值计算当前时刻的状态估计值。
3.如权利要求2所述的弹道导弹目标定轨的EKF滤波方法，其特征在于，所述计算当前时刻的状态估计值进一步包括：
当未能获得上一时刻的状态估计值时，利用所述上一时刻的状态预报值、所述上一时刻的P矩阵和所述当前时刻的量测值计算当前时刻的状态估计值；以及
当能够获得上一时刻的状态估计值时，利用所述上一时刻的状态估计值、所述上一时刻的P矩阵和所述当前时刻的量测值计算当前时刻的状态估计值。
4.如权利要求1-3中的任意一项所述的弹道导弹目标定轨的EKF滤波方法，其特征在于，计算所述当前时刻的状态估计值的方程为：
Figure FDA0003590867600000023
其中，k为上一时刻，k+1为当前时刻，
Figure FDA0003590867600000031
表示对当前时刻状态的预报值，
Figure FDA0003590867600000032
表示对上一时刻状态的估计值，
Figure FDA0003590867600000033
Figure FDA0003590867600000034
为上一时刻状态转移矩阵的预报值，Bs为增益矩阵，Kk+1为当前时刻的新息增益矩阵，C＝[I 0]，I为单位矩阵，Pk+1为当前时刻的协方差矩阵，
Figure FDA0003590867600000035
表示雷达当前时刻量测在地球系下的方差阵，
Figure FDA0003590867600000036
为当前时刻量测在地球系下的投影，Qk为上一时刻系统误差值的中间参量。
CN201811000966.XA 2018-08-30 2018-08-30 一种弹道导弹目标定轨的ekf滤波方法 Active CN108872973B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811000966.XA CN108872973B (zh) 2018-08-30 2018-08-30 一种弹道导弹目标定轨的ekf滤波方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811000966.XA CN108872973B (zh) 2018-08-30 2018-08-30 一种弹道导弹目标定轨的ekf滤波方法

Publications (2)

Publication Number Publication Date
CN108872973A CN108872973A (zh) 2018-11-23
CN108872973B true CN108872973B (zh) 2022-07-29

Family

ID=64322410

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811000966.XA Active CN108872973B (zh) 2018-08-30 2018-08-30 一种弹道导弹目标定轨的ekf滤波方法

Country Status (1)

Country Link
CN (1) CN108872973B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109990789A (zh) * 2019-03-27 2019-07-09 广东工业大学 一种飞行导航方法、装置及相关设备

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102338870A (zh) * 2011-08-25 2012-02-01 北京理工大学 一种采用前向散射雷达的三维目标跟踪方法
CN103438892A (zh) * 2013-09-16 2013-12-11 哈尔滨工业大学 一种改进的基于ekf的天文自主定轨算法
CN105022035A (zh) * 2015-07-31 2015-11-04 中国电子科技集团公司第三十八研究所 一种基于模型修正的弹道目标发射点估计装置及其方法
CN107315171A (zh) * 2017-07-02 2017-11-03 中国航空工业集团公司雷华电子技术研究所 一种雷达组网目标状态与系统误差联合估计算法
CN107765242A (zh) * 2017-09-16 2018-03-06 太原理工大学 基于状态增广迭代扩展卡尔曼滤波的系统状态估计方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10175348B2 (en) * 2014-10-08 2019-01-08 Src, Inc. Use of range-rate measurements in a fusion tracking system via projections

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102338870A (zh) * 2011-08-25 2012-02-01 北京理工大学 一种采用前向散射雷达的三维目标跟踪方法
CN103438892A (zh) * 2013-09-16 2013-12-11 哈尔滨工业大学 一种改进的基于ekf的天文自主定轨算法
CN105022035A (zh) * 2015-07-31 2015-11-04 中国电子科技集团公司第三十八研究所 一种基于模型修正的弹道目标发射点估计装置及其方法
CN107315171A (zh) * 2017-07-02 2017-11-03 中国航空工业集团公司雷华电子技术研究所 一种雷达组网目标状态与系统误差联合估计算法
CN107765242A (zh) * 2017-09-16 2018-03-06 太原理工大学 基于状态增广迭代扩展卡尔曼滤波的系统状态估计方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
A Direct Quadrature Based Nonlinear Filtering with Extended Kalman Filter Update for Orbit Determination;Jangho Yoon et al.;《Proceedings of the 2010 American Control Conference》;20100729;第2855—2860页 *
一种基于多模型算法的纯弹道式弹道落点预报方法;陈映等;《宇航学报》;20100731;第31卷(第7期);第1825—1831页 *
利用最优定轨算法鉴别弹道有源假目标;饶彬等;《系统工程与电子技术》;20100630;第32卷(第6期);第1195—1199页 *
基于EKF和多信息源融合的空间交会对接过程实时定轨方法;马鹏斌等;《载人航天》;20111231;第17卷(第3期);第27—30页 *

Also Published As

Publication number Publication date
CN108872973A (zh) 2018-11-23

Similar Documents

Publication Publication Date Title
EP2972475B1 (en) Iterative kalman filtering
CN108535720B (zh) 用于改进的卡尔曼滤波目标跟踪的自适应过程噪声描述
CN110889862B (zh) 一种网络传输攻击环境中多目标跟踪的组合测量方法
Lammas et al. A 6 DoF navigation algorithm for autonomous underwater vehicles
KR101908534B1 (ko) 이동체의 위치 및 자세 결정 장치 및 방법
CN110231620B (zh) 一种噪声相关系统跟踪滤波方法
CN106969767B (zh) 一种动平台传感器系统偏差的估计方法
KR102038049B1 (ko) 예상 궤적을 이용한 위치 결정 방법 및 장치
CN108872973B (zh) 一种弹道导弹目标定轨的ekf滤波方法
CN106569180B (zh) 一种基于Prony方法的方位估计算法
Ra et al. Robust horizontal line-of-sight rate estimator for sea skimming anti-ship missile with two-axis gimballed seeker
CN114779642A (zh) 一种基于新息抗差估计的gnss/ins紧组合欺骗检测方法
CN107390166B (zh) 一种自适应干扰源定位飞行校验方法
Wu et al. A sequential converted measurement Kalman filter with Doppler measurements in ECEF coordinate system
CN114488104B (zh) 基于交互一致性的天波超视距雷达目标跟踪方法
Janczak et al. Data fusion for ballistic targets tracking using least squares
CN111679270B (zh) 一种反射点不确定场景下多路径融合目标检测算法
Li et al. Tracking an underwater maneuvering target using an adaptive Kalman filter
KR102252823B1 (ko) 표적을 추적하고 탄두를 방출하는 장치 및 그 방법
CN113933798A (zh) 一种基于相似性原理的全局传感器系统误差分区配准算法
Kumru et al. Performance evaluation of nonlinear filters on impact point prediction of ballistic targets
Fathi et al. Adaptive Fusion of Inertial Navigation System and Tracking Radar Data
US6232914B1 (en) Method of and apparatus for determining the relative weight and weapon class of battlefield projectiles insensitive to errors in meteorological data and radar measurements
US20220334238A1 (en) Method and Device for Estimating a Velocity of an Object
CN113514823B (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