CN113108782B - 一种海空旋转调制惯导/天文组合导航方法 - Google Patents

一种海空旋转调制惯导/天文组合导航方法 Download PDF

Info

Publication number
CN113108782B
CN113108782B CN202110502728.4A CN202110502728A CN113108782B CN 113108782 B CN113108782 B CN 113108782B CN 202110502728 A CN202110502728 A CN 202110502728A CN 113108782 B CN113108782 B CN 113108782B
Authority
CN
China
Prior art keywords
coordinate system
grid
error
navigation
representing
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
CN202110502728.4A
Other languages
English (en)
Other versions
CN113108782A (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.)
National University of Defense Technology
Original Assignee
National University of Defense Technology
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 National University of Defense Technology filed Critical National University of Defense Technology
Priority to CN202110502728.4A priority Critical patent/CN113108782B/zh
Publication of CN113108782A publication Critical patent/CN113108782A/zh
Application granted granted Critical
Publication of CN113108782B publication Critical patent/CN113108782B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01CMEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/10Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration
    • G01C21/12Navigation; 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/16Navigation; 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/165Navigation; 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
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01CMEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/02Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by astronomical means

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)

Abstract

本发明属于导航技术领域,公开了一种海空旋转调制惯导/天文组合导航方法,适用于航海、航空领域的全球导航航行。本发明以地理坐标系、格网坐标系下的组合导航滤波器为基础,建立了系统误差状态及其协方差矩阵在两个导航坐标系之间的转换关系,能够实现滤波状态的稳定平滑过渡,避免坐标系转换过程中的滤波状态跳变问题;以CNS提供的姿态信息为观测量,对RINS的姿态误差进行估计校正,并能保证导航过程自主性。此外,设计的滤波器采用了开环结构,能够保证RINS、CNS系统导航信息的独立性,满足了大船、大飞机的全球安全可靠航行。

Description

一种海空旋转调制惯导/天文组合导航方法
技术领域
本发明属于导航技术领域,涉及惯性/天文组合导航方法,特别涉及一种海空旋转调制惯导/天文组合导航方法,适用于航海、航空领域的全球导航航行。
背景技术
近年来,航海、航空制造业进步迅速,大船、大飞机逐步在军民领域得到了越来越多的应用。对于大船、大飞机而言,安全可靠的全球导航能力是衡量其先进性的重要指标。为此,在大船、大飞机导航系统设计时,需要对导航系统的设计问题进行针对性研究。
大船、大飞机对导航的精度指标、可靠性有更高的要求,需要着重考虑这两方面的因素。近年来,旋转调制惯导(Rotational Inertial Navigation System,RINS)在航海、航空等长航时领域得到了越来越多的应用,旋转调制惯导能够抵消惯性器件的确定性误差,进而提高导航精度,相较于一般的纯捷联惯导,其导航精度可以提高一个数量级,但其存在长期定位误差发散的问题,需要外界其它导航系统的辅助才能抑制误差发散问题;天文导航系统(Celestial Navigaiton System,CNS)能够提供载体的姿态信息,具有长期稳定性高、误差不累积的特点,并且天文导航系统具备导航自主性,不受电磁环境影响,但天文导航系统存在易受天气因素影响的缺点,并且不能提供全部的导航参数;旋转调制惯导与天文导航系统结合起来构成RINS/CNS组合导航系统是一种十分理想的方案,具有完全的自主性、优势明显,但须针对可靠、全球导航进行特别设计。
现有针对大船、大飞机全球导航的研究主要关注的仅仅是在高纬度地区的区域导航能力,对大船、大飞机在不同纬度、不同区域之间的连续航行过程缺少足够认识。目前,在中低纬度地区,RINS/CNS组合导航算法一般在当地水平地理坐标系下设计,在高纬度地区一般在格网坐标系下设计。当大船、大飞机在两个地区之间连续航行时,组合导航算法需要在不同坐标系之间转换,以实现组合导航滤波器的一致估计,避免滤波状态震荡,而这正是现有技术忽视的地方。另一方面,为了保证导航信息的可靠性,在设计组合导航算法时,要保证RINS与CNS之间信息的独立性,传统的闭环反馈滤波方法难以适用。
本发明针对目前存在的问题,提出一种海空旋转调制惯导/天文组合导航方法,以地理坐标系、格网坐标系下的组合导航滤波器为基础,建立了系统误差状态及其协方差矩阵在两个导航坐标系之间的转换关系,能够实现滤波状态的稳定平滑过渡,避免坐标系转换过程中的滤波状态跳变问题;以CNS提供的姿态信息为观测量,对RINS的姿态误差进行估计校正,并能保证导航过程自主性。此外,设计的滤波器采用了开环结构,能够保证RINS、CNS系统导航信息的独立性,满足了大船、大飞机的全球安全可靠航行。
发明内容
本发明要解决的技术问题就在于:提供完全自主的全球导航方案,解决大船、大飞机全球航行过程中导航坐标系转换导致的滤波不稳定问题,实现系统误差状态的平滑过渡,提高导航精度,并且同时保证旋转调制惯导系统与天文导航系统的独立性,为大船、大飞机安全可靠航行提供更加准确的导航信息。
为解决上述技术问题,本发明提出的解决方案为:
一种海空旋转调制惯导/天文组合导航方法,包括以下步骤:
(1)确定大船、大飞机在高纬度航行区域的导航坐标系及位置表示方式,包括如下步骤:
(1.1)确定大船、大飞机高纬度地区航行时的导航坐标系,高纬度地区导航坐标系确定为格网坐标系,其中,格网坐标系的定义为:格网平面平行于格林尼治子午面,其与大船、大飞机位置点处切平面的交线为格网北向,地理北向与格网北向的夹角为格网角,以顺时针为正;格网天向与当地地理坐标系天向相同,其与格网东向、北向一起构成右手直角坐标系;将格网角σ表示为
Figure BDA0003057067400000021
Figure BDA0003057067400000022
其中,表示L当地纬度,λ表示当地经度;
(1.2)确定大船、大飞机在格网坐标系下的位置矩阵
Figure BDA0003057067400000023
与高度h,其中位置矩阵
Figure BDA0003057067400000024
定义为格网坐标系G与地球坐标系e之间的方向余弦矩阵,高度h即大船、大飞机相对于水平面的高度,
Figure BDA0003057067400000025
表示如下:
Figure BDA0003057067400000026
其中,
Figure BDA0003057067400000027
表示格网坐标系相对于地理坐标系n的方向余弦矩阵,
Figure BDA0003057067400000028
表示地理坐标系n相对于地球坐标系e的方向余弦矩阵;
(2)确定大船、大飞机在格网坐标系下的更新方程,包括姿态更新方程、速度更新方程、位置更新方程,具体实施如下:
(2.1)确定格网坐标系下的姿态更新方程为:
Figure BDA0003057067400000029
其中,
Figure BDA00030570674000000210
表示格网坐标系相对于载体坐标系b的方向余弦矩阵,
Figure BDA00030570674000000211
表示载体坐标系相对于惯性坐标系i的旋转角速度,
Figure BDA00030570674000000212
表示格网坐标系相对于惯性坐标系的旋转角速度;
(2.2)确定格网坐标系下的速度vG的更新方程为:
Figure BDA00030570674000000213
其中,
Figure BDA00030570674000000214
Figure BDA0003057067400000031
式中,fb表示载体坐标系下表示的比力,gG表示格网坐标系下表示的重力矢量,
Figure BDA0003057067400000032
表示地球坐标系相对于惯性坐标系的旋转角速度在格网坐标系下的投影,
Figure BDA0003057067400000033
表示格网坐标系相对于地球坐标系的旋转角速度在格网坐标系下的投影,
Figure BDA0003057067400000034
表示地球坐标系相对于惯性坐标系的旋转角速度在地球坐标系下的投影,ωie表示地球旋转角速度,Rx为格网东向的曲率半径,Ry为格网北向的曲率半径,τf为扭曲半径,
Figure BDA0003057067400000035
表示格网东向速度,
Figure BDA0003057067400000036
表示格网北向速度;
(2.3)确定格网坐标系下的位置更新方程为:
Figure BDA0003057067400000037
Figure BDA0003057067400000038
其中,位置更新包括位置矩阵
Figure BDA0003057067400000039
的更新与高度h的更新,
Figure BDA00030570674000000310
表示格网垂向速度;
(3)确定大船、大飞机在格网坐标系下的姿态误差方程、速度误差方程、位置误差方程,具体实施如下:
确定姿态误差φG的方程如下:
Figure BDA00030570674000000311
其中,
Figure BDA00030570674000000312
表示格网坐标系相对于惯性坐标系的旋转角速度误差,
Figure BDA00030570674000000313
表示载体坐标系相对于惯性坐标系的旋转角速度误差;
确定速度误差δvG的方程如下:
Figure BDA00030570674000000314
其中,
Figure BDA00030570674000000315
表示地球坐标系相对于惯性坐标系的旋转角速度误差,
Figure BDA00030570674000000316
表示格网坐标系相对于地球坐标系的旋转角速度误差,δfb表示比力误差;
确定位置误差方程,位置误差包括位置矩阵误差θG与高度误差δh,且位置矩阵误差方程采用位置误差角的微分方程表示:
Figure BDA00030570674000000317
高度误差方程为:
Figure BDA00030570674000000318
式中,
Figure BDA00030570674000000319
表示格网垂向速度误差;
(4)分别确定RINS/CNS组合导航滤波器在地理坐标系与格网坐标系下的观测方程;
在地里坐标系下的观测方程为:
Figure BDA0003057067400000041
且地理坐标系下表示的东向、北向、垂向姿态误差
Figure BDA0003057067400000042
分别为0.5(Π3223)、0.5(Π1331)、0.5(Π2112);
在格网坐标系下的观测方程为:
Figure BDA0003057067400000043
且格网坐标系下表示的东向、北向、垂向姿态误差
Figure BDA0003057067400000044
分别为0.5(Δ3223)、0.5(Δ1331)、0.5(Δ2112);
其中,
Figure BDA0003057067400000045
分别表示
Figure BDA0003057067400000046
的解算值,φn、φG分别表示地理坐标系下姿态误差与格网坐标系下姿态误差,
Figure BDA0003057067400000047
表示惯性坐标系与载体坐标系之间的方向余弦矩阵,
Figure BDA0003057067400000048
由天文导航系统提供,
Figure BDA0003057067400000049
表示地球坐标系与惯性坐标系之间的方向余弦矩阵;
(5)确定大船、大飞机导航参数在地理坐标系与格网坐标系之间的转换关系并进行转换,导航参数的转换包括姿态转换、速度转换、位置转换;
其中,大船、大飞机姿态参数在地理坐标系、格网坐标系之间的转换关系为:
Figure BDA00030570674000000410
Figure BDA00030570674000000411
式中,
Figure BDA00030570674000000412
表示地理坐标系n与载体坐标系b之间的方向余弦矩阵,
Figure BDA00030570674000000413
表示地理坐标系与格网坐标系之间的方向余弦矩阵;
大船、大飞机速度参数在地理坐标系、格网坐标系之间的转换关系为:
Figure BDA00030570674000000414
Figure BDA00030570674000000415
式中,vn表示地理坐标系下表示的速度;
大船、大飞机位置参数在地理坐标系、格网坐标系之间的转换关系为:
Figure BDA00030570674000000416
离开高纬度地区时,纬度、经度通过位置矩阵
Figure BDA00030570674000000417
的元素c31,c32,c33通过三角函数运算获得,其中c31,c32,c33分别为
Figure BDA00030570674000000418
的第3行第1-3列元素;
高度h在两个坐标系下保持不变;
(6)完成RINS/CNS组合导航滤波器在地理坐标系与格网坐标系之间的转换,其中RINS/CNS组合导航滤波器采用开环反馈校正方式,具体实施如下:
(6.1)分别确定地理坐标系与格网坐标系下的系统误差状态为:
地理坐标系下的系统误差状态xn(t)为
Figure BDA0003057067400000051
格网坐标系下的系统误差状态xG(t)为
Figure BDA0003057067400000052
其中,
Figure BDA0003057067400000053
分别表示地理坐标系下表示的东向、北向、垂向姿态误差,
Figure BDA0003057067400000054
分别表示格网坐标系下表示的东向、北向、垂向姿态误差,
Figure BDA0003057067400000055
分别表示地理坐标系下表示的东向、北向、垂向速度误差,
Figure BDA0003057067400000056
分别表示格网坐标系下表示的东向、北向、垂向速度误差,δL,δλ分别表示纬度、经度误差,
Figure BDA0003057067400000057
分别表示位置误差角东向、北向误差,
Figure BDA0003057067400000058
分别表示x、y、z轴向陀螺常值零偏,
Figure BDA0003057067400000059
分别表示x、y、z轴向加表常值零偏;
(6.2)分别确定姿态误差、速度误差、位置误差在地理坐标系与格网坐标系下间的转换关系为:
首先确定地理坐标系下姿态误差φn与格网坐标系下姿态误差φG之间的转换关系
Figure BDA00030570674000000510
式中,
Figure BDA00030570674000000511
Figure BDA00030570674000000512
其次确定地理坐标系下速度误差δvn与格网坐标系下速度误差δvG之间的转换关系
Figure BDA00030570674000000513
式中,
Figure BDA00030570674000000514
表示格网坐标系相对于地理坐标系方向余弦矩阵的误差;
进而确定纬度误差δL、经度误差δλ与示位置误差角东向误差
Figure BDA00030570674000000515
北向误差
Figure BDA00030570674000000516
之间的转换关系
Figure BDA00030570674000000517
高度误差δh、陀螺常值零偏
Figure BDA00030570674000000518
加表常值零偏
Figure BDA00030570674000000519
在地理坐标系与格网坐标系下保持不变;
确定格网坐标系下的系统误差状态xG(t)与地理坐标系下的系统误差状态xn(t)之间的转换关系如下:
xG(t)=Φxn(t),xn(t)=Φ-1xG(t)
其中,Φ为转换系数矩阵,并且根据φG与φn之间的转换关系,δvG与δvn之间的转换关系,
Figure BDA00030570674000000520
与δL、δλ之间的转换关系,并考虑高度误差δh、陀螺常值零偏
Figure BDA00030570674000000521
加表常值零偏
Figure BDA0003057067400000061
在地理坐标系与格网坐标系下的不变性进行确定;
(6.3)根据步骤(6.2),确定地理坐标系下系统误差状态协方差矩阵Pn(t)与格网坐标系下系统误差状态协方差矩阵PG(t)的转换关系:
Figure BDA0003057067400000062
Pn(t)=Φ-1PG(t)Φ-T
式中,
Figure BDA0003057067400000063
表示格网坐标系下表示的系统误差状态估计值,
Figure BDA0003057067400000064
表示地理坐标系下表示的系统误差状态估计值;
(6.4)根据步骤(4),当大船、大飞机处于中低纬度时,采用地理坐标系下的观测方程,当大船、大飞机处于高纬度时,采用格网坐标系下的观测方程,并且观测方程与系统状态方程对应,系统状态方程确定后,观测方程相应确定;
(6.5)当大船、大飞机在中纬度、高纬度地区连续航行时,开环反馈RINS/CNS组合导航滤波器完成在地理坐标系与格网坐标系之间的系统误差状态、协方差矩阵转换,转换方式按照步骤(6.2)、步骤(6.3)所述,转换前后xn(t)、Pn(t),xG(t)、PG(t)按照如下方式进行更新:
Figure BDA0003057067400000065
Figure BDA0003057067400000066
Figure BDA0003057067400000067
Figure BDA0003057067400000068
Figure BDA0003057067400000069
式中,上标+、-分别表示更新后时刻、更新前时刻,下标k+1、k分别表示离散化k+1、k时刻,K、P、H、R、Q、F、Υ分别表示增益矩阵、协方差矩阵、观测矩阵、观测噪声强度矩阵、系统噪声强度矩阵、状态转移矩阵、系统噪声矩阵,x、z分别表示系统状态向量、观测向量,I为单位矩阵;
(7)采用输出校正方式对RINS导航参数信息进行校正,在地理坐标系、格网坐标系下的导航参数校正方式分别如下:
Figure BDA00030570674000000610
Figure BDA00030570674000000611
Figure BDA00030570674000000612
式中,
Figure BDA00030570674000000613
分别表示
Figure BDA00030570674000000614
的解算值,
Figure BDA00030570674000000615
分别表示vn、vG的解算值,
Figure BDA00030570674000000616
Figure BDA00030570674000000617
分别表示L、λ、h的解算值,
Figure BDA00030570674000000618
Figure BDA00030570674000000619
的解算值。
进一步的,若大船、大飞机接收到GNSS定位信息时,利用GNSS位置点信息完成对所述步骤(4)中
Figure BDA0003057067400000071
Figure BDA0003057067400000072
的装订更新。
进一步的,所述步骤(5)中导航参数在地理坐标系与格网坐标系之间转换时基于转换时刻的纬度阈值判断,且地理坐标系转换到格网坐标系、格网坐标系转换到地理坐标系两种情况下的阈值设定不同。
进一步的,当CNS受环境因素影响无法输出姿态参考信息时,所述步骤(6.5)只进行RINS/CNS组合导航滤波的时间更新过程。
进一步的,所述的RINS为单轴旋转调制惯导、或双轴旋转调制惯导、或三轴旋转调制惯导。
通过以上步骤可以实现大船、大飞机RINS/CNS全球组合导航方法,实现全球范围内的全自主导航定位,不会出现导航滤波器震荡问题,并且能够保证RINS、CNS导航信息的独立性。
与现有技术相比,本发明的优点在于:
(1)本发明解决了大船、大飞机全球航行过程中,由于导航坐标系转换带来的组合导航滤波器震荡问题,有效提高导航精度。
(2)本发明设计滤波器能够保证RINS、CNS导航信息的独立性,实现全自主导航定位。
附图说明
图1为本发明方法的流程示意图。
具体实施方式
以下将结合说明书附图和具体实施例对本发明作进一步详细说明。
如图1所示,一种海空旋转调制惯导/天文组合导航方法,包括以下步骤:
(1)确定大船、大飞机在高纬度航行区域的导航坐标系及位置表示方式,包括如下步骤:
(1.1)确定大船、大飞机高纬度地区航行时的导航坐标系,高纬度地区导航坐标系确定为格网坐标系,其中,格网坐标系的定义为:格网平面平行于格林尼治子午面,其与大船、大飞机位置点处切平面的交线为格网北向,地理北向与格网北向的夹角为格网角,以顺时针为正;格网天向与当地地理坐标系天向相同,其与格网东向、北向一起构成右手直角坐标系;将格网角σ表示为
Figure BDA0003057067400000073
Figure BDA0003057067400000074
其中,表示L当地纬度,λ表示当地经度;
(1.2)确定大船、大飞机在格网坐标系下的位置矩阵
Figure BDA0003057067400000075
与高度h,其中位置矩阵
Figure BDA0003057067400000076
定义为格网坐标系G与地球坐标系e之间的方向余弦矩阵,高度h即大船、大飞机相对于水平面的高度,
Figure BDA0003057067400000077
表示如下:
Figure BDA0003057067400000081
其中,
Figure BDA0003057067400000082
表示格网坐标系相对于地理坐标系n的方向余弦矩阵,
Figure BDA0003057067400000083
表示地理坐标系n相对于地球坐标系e的方向余弦矩阵;
(2)确定大船、大飞机在格网坐标系下的更新方程,包括姿态更新方程、速度更新方程、位置更新方程,具体实施如下:
(2.1)确定格网坐标系下的姿态更新方程为:
Figure BDA0003057067400000084
其中,
Figure BDA0003057067400000085
表示格网坐标系相对于载体坐标系b的方向余弦矩阵,
Figure BDA0003057067400000086
表示载体坐标系相对于惯性坐标系i的旋转角速度,
Figure BDA0003057067400000087
表示格网坐标系相对于惯性坐标系的旋转角速度;
(2.2)确定格网坐标系下的速度vG的更新方程为:
Figure BDA0003057067400000088
其中,
Figure BDA0003057067400000089
Figure BDA00030570674000000810
式中,fb表示载体坐标系下表示的比力,gG表示格网坐标系下表示的重力矢量,
Figure BDA00030570674000000811
表示地球坐标系相对于惯性坐标系的旋转角速度在格网坐标系下的投影,
Figure BDA00030570674000000812
表示格网坐标系相对于地球坐标系的旋转角速度在格网坐标系下的投影,
Figure BDA00030570674000000813
表示地球坐标系相对于惯性坐标系的旋转角速度在地球坐标系下的投影,ωie表示地球旋转角速度,Rx为格网东向的曲率半径,Ry为格网北向的曲率半径,τf为扭曲半径,
Figure BDA00030570674000000814
表示格网东向速度,
Figure BDA00030570674000000815
表示格网北向速度;
(2.3)确定格网坐标系下的位置更新方程为:
Figure BDA00030570674000000816
Figure BDA00030570674000000817
其中,位置更新包括位置矩阵
Figure BDA00030570674000000818
的更新与高度h的更新,
Figure BDA00030570674000000819
表示格网垂向速度;
(3)确定大船、大飞机在格网坐标系下的姿态误差方程、速度误差方程、位置误差方程,具体实施如下:
确定姿态误差φG的方程如下:
Figure BDA00030570674000000820
其中,
Figure BDA0003057067400000091
表示格网坐标系相对于惯性坐标系的旋转角速度误差,
Figure BDA0003057067400000092
表示载体坐标系相对于惯性坐标系的旋转角速度误差;
确定速度误差δvG的方程如下:
Figure BDA0003057067400000093
其中,
Figure BDA0003057067400000094
表示地球坐标系相对于惯性坐标系的旋转角速度误差,
Figure BDA0003057067400000095
表示格网坐标系相对于地球坐标系的旋转角速度误差,δfb表示比力误差;
确定位置误差方程,位置误差包括位置矩阵误差θG与高度误差δh,且位置矩阵误差方程采用位置误差角的微分方程表示:
Figure BDA0003057067400000096
高度误差方程为:
Figure BDA0003057067400000097
式中,
Figure BDA0003057067400000098
表示格网垂向速度误差;
(4)分别确定RINS/CNS组合导航滤波器在地理坐标系与格网坐标系下的观测方程;
在地里坐标系下的观测方程为:
Figure BDA0003057067400000099
且地理坐标系下表示的东向、北向、垂向姿态误差
Figure BDA00030570674000000910
分别为0.5(Π3223)、0.5(Π1331)、0.5(Π2112);
在格网坐标系下的观测方程为:
Figure BDA00030570674000000911
且格网坐标系下表示的东向、北向、垂向姿态误差
Figure BDA00030570674000000912
分别为0.5(Δ3223)、0.5(Δ1331)、0.5(Δ2112);
其中,
Figure BDA00030570674000000913
分别表示
Figure BDA00030570674000000914
的解算值,φn、φG分别表示地理坐标系下姿态误差与格网坐标系下姿态误差,
Figure BDA00030570674000000915
表示惯性坐标系与载体坐标系之间的方向余弦矩阵,
Figure BDA00030570674000000916
由天文导航系统提供,
Figure BDA00030570674000000917
表示地球坐标系与惯性坐标系之间的方向余弦矩阵;
(5)确定大船、大飞机导航参数在地理坐标系与格网坐标系之间的转换关系并进行转换,导航参数的转换包括姿态转换、速度转换、位置转换;
其中,大船、大飞机姿态参数在地理坐标系、格网坐标系之间的转换关系为:
Figure BDA00030570674000000918
Figure BDA00030570674000000919
式中,
Figure BDA00030570674000000920
表示地理坐标系n与载体坐标系b之间的方向余弦矩阵,
Figure BDA00030570674000000921
表示地理坐标系与格网坐标系之间的方向余弦矩阵;
大船、大飞机速度参数在地理坐标系、格网坐标系之间的转换关系为:
Figure BDA0003057067400000101
Figure BDA0003057067400000102
式中,vn表示地理坐标系下表示的速度;
大船、大飞机位置参数在地理坐标系、格网坐标系之间的转换关系为:
Figure BDA0003057067400000103
离开高纬度地区时,纬度、经度通过位置矩阵
Figure BDA0003057067400000104
的元素c31,c32,c33通过三角函数运算获得,其中c31,c32,c33分别为
Figure BDA0003057067400000105
的第3行第1-3列元素;
高度h在两个坐标系下保持不变;
(6)完成RINS/CNS组合导航滤波器在地理坐标系与格网坐标系之间的转换,其中RINS/CNS组合导航滤波器采用开环反馈校正方式,具体实施如下:
(6.1)分别确定地理坐标系与格网坐标系下的系统误差状态为:
地理坐标系下的系统误差状态xn(t)为
Figure BDA0003057067400000106
格网坐标系下的系统误差状态xG(t)为
Figure BDA0003057067400000107
其中,
Figure BDA0003057067400000108
分别表示地理坐标系下表示的东向、北向、垂向姿态误差,
Figure BDA0003057067400000109
分别表示格网坐标系下表示的东向、北向、垂向姿态误差,
Figure BDA00030570674000001010
分别表示地理坐标系下表示的东向、北向、垂向速度误差,
Figure BDA00030570674000001011
分别表示格网坐标系下表示的东向、北向、垂向速度误差,δL,δλ分别表示纬度、经度误差,
Figure BDA00030570674000001012
分别表示位置误差角东向、北向误差,
Figure BDA00030570674000001013
分别表示x、y、z轴向陀螺常值零偏,
Figure BDA00030570674000001014
分别表示x、y、z轴向加表常值零偏;
(6.2)分别确定姿态误差、速度误差、位置误差在地理坐标系与格网坐标系下间的转换关系为:
首先确定地理坐标系下姿态误差φn与格网坐标系下姿态误差φG之间的转换关系
Figure BDA00030570674000001015
式中,
Figure BDA00030570674000001016
Figure BDA00030570674000001017
其次确定地理坐标系下速度误差δvn与格网坐标系下速度误差δvG之间的转换关系
Figure BDA0003057067400000111
式中,
Figure BDA0003057067400000112
表示格网坐标系相对于地理坐标系方向余弦矩阵的误差;
进而确定纬度误差δL、经度误差δλ与示位置误差角东向误差
Figure BDA0003057067400000113
北向误差
Figure BDA0003057067400000114
之间的转换关系
Figure BDA0003057067400000115
高度误差δh、陀螺常值零偏
Figure BDA0003057067400000116
加表常值零偏
Figure BDA0003057067400000117
在地理坐标系与格网坐标系下保持不变;
确定格网坐标系下的系统误差状态xG(t)与地理坐标系下的系统误差状态xn(t)之间的转换关系如下:
xG(t)=Φxn(t),xn(t)=Φ-1xG(t)
其中,Φ为转换系数矩阵,并且根据φG与φn之间的转换关系,δvG与δvn之间的转换关系,
Figure BDA0003057067400000118
与δL、δλ之间的转换关系,并考虑高度误差δh、陀螺常值零偏
Figure BDA0003057067400000119
加表常值零偏
Figure BDA00030570674000001110
在地理坐标系与格网坐标系下的不变性进行确定;
(6.3)根据步骤(6.2),确定地理坐标系下系统误差状态协方差矩阵Pn(t)与格网坐标系下系统误差状态协方差矩阵PG(t)的转换关系:
Figure BDA00030570674000001111
Pn(t)=Φ-1PG(t)Φ-T
式中,
Figure BDA00030570674000001112
表示格网坐标系下表示的系统误差状态估计值,
Figure BDA00030570674000001113
表示地理坐标系下表示的系统误差状态估计值;
(6.4)根据步骤(4),当大船、大飞机处于中低纬度时,采用地理坐标系下的观测方程,当大船、大飞机处于高纬度时,采用格网坐标系下的观测方程,并且观测方程与系统状态方程对应,系统状态方程确定后,观测方程相应确定;
(6.5)当大船、大飞机在中纬度、高纬度地区连续航行时,开环反馈RINS/CNS组合导航滤波器完成在地理坐标系与格网坐标系之间的系统误差状态、协方差矩阵转换,转换方式按照步骤(6.2)、步骤(6.3)所述,转换前后xn(t)、Pn(t),xG(t)、PG(t)按照如下方式进行更新:
Figure BDA0003057067400000121
Figure BDA0003057067400000122
Figure BDA0003057067400000123
Figure BDA0003057067400000124
Figure BDA0003057067400000125
式中,上标+、-分别表示更新后时刻、更新前时刻,下标k+1、k分别表示离散化k+1、k时刻,K、P、H、R、Q、F、Υ分别表示增益矩阵、协方差矩阵、观测矩阵、观测噪声强度矩阵、系统噪声强度矩阵、状态转移矩阵、系统噪声矩阵,x、z分别表示系统状态向量、观测向量,I为单位矩阵;
(7)采用输出校正方式对RINS导航参数信息进行校正,在地理坐标系、格网坐标系下的导航参数校正方式分别如下:
Figure BDA0003057067400000126
Figure BDA0003057067400000127
Figure BDA0003057067400000128
式中,
Figure BDA0003057067400000129
分别表示
Figure BDA00030570674000001210
的解算值,
Figure BDA00030570674000001211
分别表示vn、vG的解算值,
Figure BDA00030570674000001212
Figure BDA00030570674000001213
分别表示L、λ、h的解算值,
Figure BDA00030570674000001214
Figure BDA00030570674000001215
的解算值。
进一步的,若大船、大飞机接收到GNSS定位信息时,利用GNSS位置点信息完成对所述步骤(4)中
Figure BDA00030570674000001216
Figure BDA00030570674000001217
的装订更新。
进一步的,所述步骤(5)中导航参数在地理坐标系与格网坐标系之间转换时基于转换时刻的纬度阈值判断,且地理坐标系转换到格网坐标系、格网坐标系转换到地理坐标系两种情况下的阈值设定不同。
进一步的,当CNS受环境因素影响无法输出姿态参考信息时,所述步骤(6.5)只进行RINS/CNS组合导航滤波的时间更新过程。
进一步的,所述的RINS为单轴旋转调制惯导、或双轴旋转调制惯导、或三轴旋转调制惯导。
以上所述仅是本发明的优选实施方式,本发明的保护范围并不仅局限于上述实施例,凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出,对于本技术领域的普通技术人员来说,在不脱离本发明原理前提下的若干改进和润饰,这些改进和润饰也应视为本发明的保护范围。

Claims (5)

1.一种海空旋转调制惯导/天文组合导航方法,其特征在于,包括以下步骤:
(1)确定大船、大飞机在高纬度航行区域的导航坐标系及位置表示方式,包括如下步骤:
(1.1)确定大船、大飞机高纬度地区航行时的导航坐标系,高纬度地区导航坐标系确定为格网坐标系,其中,格网坐标系的定义为:格网平面平行于格林尼治子午面,其与大船、大飞机位置点处切平面的交线为格网北向,地理北向与格网北向的夹角为格网角,以顺时针为正;格网天向与当地地理坐标系天向相同,其与格网东向、北向一起构成右手直角坐标系;将格网角σ表示为
Figure FDA0003621956300000011
Figure FDA0003621956300000012
其中,表示L当地纬度,λ表示当地经度;
(1.2)确定大船、大飞机在格网坐标系下的位置矩阵
Figure FDA0003621956300000013
与高度h,其中位置矩阵
Figure FDA0003621956300000014
定义为格网坐标系G与地球坐标系e之间的方向余弦矩阵,高度h即大船、大飞机相对于水平面的高度,
Figure FDA0003621956300000015
表示如下:
Figure FDA0003621956300000016
其中,
Figure FDA0003621956300000017
表示格网坐标系相对于地理坐标系n的方向余弦矩阵,
Figure FDA0003621956300000018
表示地理坐标系n相对于地球坐标系e的方向余弦矩阵;
(2)确定大船、大飞机在格网坐标系下的更新方程,包括姿态更新方程、速度更新方程、位置更新方程,具体实施如下:
(2.1)确定格网坐标系下的姿态更新方程为:
Figure FDA0003621956300000019
其中,
Figure FDA00036219563000000110
表示格网坐标系相对于载体坐标系b的方向余弦矩阵,
Figure FDA00036219563000000111
表示载体坐标系相对于惯性坐标系i的旋转角速度,
Figure FDA00036219563000000112
表示格网坐标系相对于惯性坐标系的旋转角速度;
(2.2)确定格网坐标系下的速度vG的更新方程为:
Figure FDA00036219563000000113
其中,
Figure FDA00036219563000000114
Figure FDA00036219563000000115
式中,fb表示载体坐标系下表示的比力,gG表示格网坐标系下表示的重力矢量,
Figure FDA0003621956300000021
表示地球坐标系相对于惯性坐标系的旋转角速度在格网坐标系下的投影,
Figure FDA0003621956300000022
表示格网坐标系相对于地球坐标系的旋转角速度在格网坐标系下的投影,
Figure FDA0003621956300000023
表示地球坐标系相对于惯性坐标系的旋转角速度在地球坐标系下的投影,ωie表示地球旋转角速度,Rx为格网东向的曲率半径,Ry为格网北向的曲率半径,τf为扭曲半径,
Figure FDA0003621956300000024
表示格网东向速度,
Figure FDA0003621956300000025
表示格网北向速度;
(2.3)确定格网坐标系下的位置更新方程为:
Figure FDA0003621956300000026
Figure FDA0003621956300000027
其中,位置更新包括位置矩阵
Figure FDA0003621956300000028
的更新与高度h的更新,
Figure FDA0003621956300000029
表示格网垂向速度;
(3)确定大船、大飞机在格网坐标系下的姿态误差方程、速度误差方程、位置误差方程,具体实施如下:
确定姿态误差φG的方程如下:
Figure FDA00036219563000000210
其中,
Figure FDA00036219563000000211
表示格网坐标系相对于惯性坐标系的旋转角速度误差,
Figure FDA00036219563000000212
表示载体坐标系相对于惯性坐标系的旋转角速度误差;
确定速度误差δvG的方程如下:
Figure FDA00036219563000000213
其中,
Figure FDA00036219563000000214
表示地球坐标系相对于惯性坐标系的旋转角速度误差,
Figure FDA00036219563000000215
表示格网坐标系相对于地球坐标系的旋转角速度误差,δfb表示比力误差;
确定位置误差方程,位置误差包括位置矩阵误差θG与高度误差δh,且位置矩阵误差方程采用位置误差角的微分方程表示:
Figure FDA00036219563000000216
高度误差方程为:
Figure FDA00036219563000000217
式中,
Figure FDA00036219563000000218
表示格网垂向速度误差;
(4)分别确定RINS/CNS组合导航滤波器在地理坐标系与格网坐标系下的观测方程;
在地理坐标系下的观测方程为:
Figure FDA00036219563000000219
且地理坐标系下表示的东向、北向、垂向姿态误差
Figure FDA00036219563000000220
分别为0.5(Π3223)、0.5(Π1331)、0.5(Π2112);
在格网坐标系下的观测方程为:
Figure FDA0003621956300000031
且格网坐标系下表示的东向、北向、垂向姿态误差
Figure FDA0003621956300000032
分别为0.5(Δ3223)、0.5(Δ1331)、0.5(Δ2112);
其中,
Figure FDA0003621956300000033
分别表示
Figure FDA0003621956300000034
的解算值,φn、φG分别表示地理坐标系下姿态误差与格网坐标系下姿态误差,
Figure FDA0003621956300000035
表示惯性坐标系与载体坐标系之间的方向余弦矩阵,
Figure FDA0003621956300000036
由天文导航系统提供,
Figure FDA0003621956300000037
表示地球坐标系与惯性坐标系之间的方向余弦矩阵;
(5)确定大船、大飞机导航参数在地理坐标系与格网坐标系之间的转换关系并进行转换,导航参数的转换包括姿态转换、速度转换、位置转换;
其中,大船、大飞机姿态参数在地理坐标系、格网坐标系之间的转换关系为:
Figure FDA0003621956300000038
Figure FDA0003621956300000039
式中,
Figure FDA00036219563000000310
表示地理坐标系n与载体坐标系b之间的方向余弦矩阵,
Figure FDA00036219563000000311
表示地理坐标系与格网坐标系之间的方向余弦矩阵;
大船、大飞机速度参数在地理坐标系、格网坐标系之间的转换关系为:
Figure FDA00036219563000000312
Figure FDA00036219563000000313
式中,vn表示地理坐标系下表示的速度;
大船、大飞机位置参数在地理坐标系、格网坐标系之间的转换关系为:
Figure FDA00036219563000000314
离开高纬度地区时,纬度、经度通过位置矩阵
Figure FDA00036219563000000315
的元素c31,c32,c33通过三角函数运算获得,其中c31,c32,c33分别为
Figure FDA00036219563000000316
的第3行第1-3列元素;
高度h在两个坐标系下保持不变;
(6)完成RINS/CNS组合导航滤波器在地理坐标系与格网坐标系之间的转换,其中RINS/CNS组合导航滤波器采用开环反馈校正方式,具体实施如下:
(6.1)分别确定地理坐标系与格网坐标系下的系统误差状态为:
地理坐标系下的系统误差状态xn(t)为
Figure FDA00036219563000000317
格网坐标系下的系统误差状态xG(t)为
Figure FDA00036219563000000318
其中,
Figure FDA0003621956300000041
分别表示地理坐标系下表示的东向、北向、垂向姿态误差,
Figure FDA0003621956300000042
分别表示格网坐标系下表示的东向、北向、垂向姿态误差,
Figure FDA0003621956300000043
分别表示地理坐标系下表示的东向、北向、垂向速度误差,
Figure FDA0003621956300000044
分别表示格网坐标系下表示的东向、北向、垂向速度误差,δL,δλ分别表示纬度、经度误差,
Figure FDA0003621956300000045
分别表示位置误差角东向、北向误差,
Figure FDA0003621956300000046
分别表示x、y、z轴向陀螺常值零偏,
Figure FDA0003621956300000047
分别表示x、y、z轴向加表常值零偏;
(6.2)分别确定姿态误差、速度误差、位置误差在地理坐标系与格网坐标系下间的转换关系为:
首先确定地理坐标系下姿态误差φn与格网坐标系下姿态误差φG之间的转换关系
Figure FDA0003621956300000048
式中,
Figure FDA0003621956300000049
Figure FDA00036219563000000410
其次确定地理坐标系下速度误差δvn与格网坐标系下速度误差δvG之间的转换关系
Figure FDA00036219563000000411
式中,
Figure FDA00036219563000000412
表示格网坐标系相对于地理坐标系方向余弦矩阵的误差;
进而确定纬度误差δL、经度误差δλ与示位置误差角东向误差
Figure FDA00036219563000000413
北向误差
Figure FDA00036219563000000414
之间的转换关系
Figure FDA00036219563000000415
高度误差δh、陀螺常值零偏
Figure FDA00036219563000000416
加表常值零偏
Figure FDA00036219563000000417
在地理坐标系与格网坐标系下保持不变;
确定格网坐标系下的系统误差状态xG(t)与地理坐标系下的系统误差状态xn(t)之间的转换关系如下:
xG(t)=Φxn(t),xn(t)=Φ-1xG(t)
其中,Φ为转换系数矩阵,并且根据φG与φn之间的转换关系,δvG与δvn之间的转换关系,
Figure FDA00036219563000000418
与δL、δλ之间的转换关系,并考虑高度误差δh、陀螺常值零偏
Figure FDA00036219563000000419
加表常值零偏
Figure FDA00036219563000000420
在地理坐标系与格网坐标系下的不变性进行确定;
(6.3)根据步骤(6.2),确定地理坐标系下系统误差状态协方差矩阵Pn(t)与格网坐标系下系统误差状态协方差矩阵PG(t)的转换关系:
Figure FDA0003621956300000051
Pn(t)=Φ-1PG(t)Φ-T
式中,
Figure FDA0003621956300000052
表示格网坐标系下表示的系统误差状态估计值,
Figure FDA0003621956300000053
表示地理坐标系下表示的系统误差状态估计值;
(6.4)根据步骤(4),当大船、大飞机处于中低纬度时,采用地理坐标系下的观测方程,当大船、大飞机处于高纬度时,采用格网坐标系下的观测方程,并且观测方程与系统状态方程对应,系统状态方程确定后,观测方程相应确定;
(6.5)当大船、大飞机在中纬度、高纬度地区连续航行时,开环反馈RINS/CNS组合导航滤波器完成在地理坐标系与格网坐标系之间的系统误差状态、协方差矩阵转换,转换方式按照步骤(6.2)、步骤(6.3)所述,转换前后xn(t)、Pn(t),xG(t)、PG(t)按照如下方式进行更新:
Figure FDA0003621956300000054
Figure FDA0003621956300000055
Figure FDA0003621956300000056
Figure FDA0003621956300000057
Figure FDA0003621956300000058
式中,上标+、-分别表示更新后时刻、更新前时刻,下标k+1、k分别表示离散化k+1、k时刻,K、P、H、R、Q、F、Υ分别表示增益矩阵、协方差矩阵、观测矩阵、观测噪声强度矩阵、系统噪声强度矩阵、状态转移矩阵、系统噪声矩阵,x、z分别表示系统状态向量、观测向量,I为单位矩阵;
(7)采用输出校正方式对RINS导航参数信息进行校正,在地理坐标系、格网坐标系下的导航参数校正方式分别如下:
Figure FDA0003621956300000059
式中,
Figure FDA00036219563000000510
分别表示
Figure FDA00036219563000000511
的解算值,
Figure FDA00036219563000000512
分别表示vn、vG的解算值,
Figure FDA00036219563000000513
Figure FDA00036219563000000514
分别表示L、λ、h的解算值,
Figure FDA00036219563000000515
Figure FDA00036219563000000516
的解算值。
2.如权利要求1所述的一种海空旋转调制惯导/天文组合导航方法,其特征在于,若大船、大飞机接收到GNSS定位信息时,利用GNSS位置点信息完成对所述步骤(4)中
Figure FDA00036219563000000517
Figure FDA00036219563000000518
的装订更新。
3.如权利要求1所述的一种海空旋转调制惯导/天文组合导航方法,其特征在于,所述步骤(5)中导航参数在地理坐标系与格网坐标系之间转换时基于转换时刻的纬度阈值判断,且地理坐标系转换到格网坐标系、格网坐标系转换到地理坐标系两种情况下的阈值设定不同。
4.如权利要求1所述的一种海空旋转调制惯导/天文组合导航方法,其特征在于,当CNS受环境因素影响无法输出姿态参考信息时,所述步骤(6.5)只进行RINS/CNS组合导航滤波的时间更新过程。
5.如权利要求1至4中任一项所述的一种海空旋转调制惯导/天文组合导航方法,其特征在于,所述的RINS为单轴旋转调制惯导、或双轴旋转调制惯导、或三轴旋转调制惯导。
CN202110502728.4A 2021-05-09 2021-05-09 一种海空旋转调制惯导/天文组合导航方法 Active CN113108782B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110502728.4A CN113108782B (zh) 2021-05-09 2021-05-09 一种海空旋转调制惯导/天文组合导航方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110502728.4A CN113108782B (zh) 2021-05-09 2021-05-09 一种海空旋转调制惯导/天文组合导航方法

Publications (2)

Publication Number Publication Date
CN113108782A CN113108782A (zh) 2021-07-13
CN113108782B true CN113108782B (zh) 2022-06-14

Family

ID=76721731

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110502728.4A Active CN113108782B (zh) 2021-05-09 2021-05-09 一种海空旋转调制惯导/天文组合导航方法

Country Status (1)

Country Link
CN (1) CN113108782B (zh)

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107543545A (zh) * 2017-10-30 2018-01-05 中国人民解放军国防科技大学 极区双航海惯性导航系统定位信息融合方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9341718B2 (en) * 2012-09-07 2016-05-17 Honeywell International Inc. Method and system for providing integrity for hybrid attitude and true heading

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107543545A (zh) * 2017-10-30 2018-01-05 中国人民解放军国防科技大学 极区双航海惯性导航系统定位信息融合方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Transfer alignment method for SINS based on reverse navigation solution and data fusion;Sun Jin 等;《Journal of Chinese Inertial Technology》;20151231;第23卷(第6期);第727-732页 *
航海多惯导协同定位与误差参数估计;王林;《中国博士学位论文全文数据库工程科技Ⅱ辑》;20200215(第02期);第C032-36页 *

Also Published As

Publication number Publication date
CN113108782A (zh) 2021-07-13

Similar Documents

Publication Publication Date Title
CN103245360B (zh) 晃动基座下的舰载机旋转式捷联惯导系统自对准方法
CN102692225B (zh) 一种用于低成本小型无人机的姿态航向参考系统
CN102830414B (zh) 一种基于sins/gps的组合导航方法
CN104344837B (zh) 一种基于速度观测的冗余惯导系统加速度计系统级标定方法
CN110954102B (zh) 用于机器人定位的磁力计辅助惯性导航系统及方法
CN101949703A (zh) 一种捷联惯性/卫星组合导航滤波方法
CN103217174B (zh) 一种基于低精度微机电系统的捷联惯导系统初始对准方法
CN103363992A (zh) 基于梯度下降的四旋翼无人机姿态航向参考系统解算方法
CN102519470A (zh) 多级嵌入式组合导航系统及导航方法
CN103900608A (zh) 一种基于四元数ckf的低精度惯导初始对准方法
CN113108783B (zh) 一种无人潜航器惯性/多普勒组合导航方法
CN109489661B (zh) 一种卫星初始入轨时陀螺组合常值漂移估计方法
CN110058288A (zh) 无人机ins/gnss组合导航系统航向误差修正方法及系统
CN104236586A (zh) 基于量测失准角的动基座传递对准方法
CN103674059A (zh) 一种基于外测速度信息的sins水平姿态误差修正方法
CN102393204B (zh) 一种基于sins/cns的组合导航信息融合方法
CN110398242A (zh) 一种高旋高过载条件飞行器的姿态角确定方法
CN116222551A (zh) 一种融合多种数据的水下导航方法及装置
CN113108788B (zh) 一种长航时惯导/天文全球组合导航方法
CN106643726B (zh) 一种统一惯性导航解算方法
CN113108787B (zh) 一种长航时惯导/卫星全球组合导航方法
CN115574817B (zh) 一种基于三轴旋转式惯导系统的导航方法及导航系统
CN113108782B (zh) 一种海空旋转调制惯导/天文组合导航方法
CN105180928B (zh) 一种基于惯性系重力特性的船载星敏感器定位方法
CN113155125B (zh) 一种大飞机ins/gnss全球组合导航方法

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