CN111412887A - 一种基于卡尔曼滤波的攻角、侧滑角辨识方法 - Google Patents
一种基于卡尔曼滤波的攻角、侧滑角辨识方法 Download PDFInfo
- Publication number
- CN111412887A CN111412887A CN202010242132.0A CN202010242132A CN111412887A CN 111412887 A CN111412887 A CN 111412887A CN 202010242132 A CN202010242132 A CN 202010242132A CN 111412887 A CN111412887 A CN 111412887A
- Authority
- CN
- China
- Prior art keywords
- angle
- attack
- sideslip
- obtaining
- ppre
- 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.)
- Granted
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C1/00—Measuring angles
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
Abstract
本发明提供了一种基于卡尔曼滤波的攻角、侧滑角辨识方法,包括:建立攻角α与侧滑角β的递推模型和观测模型;得到攻角α和侧滑角β的一步预测值Xpre;得到一步预测均方差Ppre;得到测量矩阵H和测量量Z;得到滤波增益K;得到卡尔曼滤波后的攻角α(k+1)和卡尔曼滤波后的侧滑角β(k+1);得到更新后的误差方差阵P′;得到攻角α和侧滑角β更新后的一步预测值Xpre′,并返回S4基于更新后的误差方差阵P′和更新后的一步预测值Xpre′进行下一轮迭代。本发明通过建立攻角α与侧滑角β的递推模型和观测模型,并基于攻角α与侧滑角β的递推模型和观测模型进行卡尔曼滤波,可实现攻角α和侧滑角β的高精度在线辨识。本发明在无攻角、侧滑角测量装置的条件下使用,节省了一套传感采集设备,而且可实现对攻角、侧滑角的同步测量。
Description
技术领域
本发明涉及航空飞行器飞行过程中的攻角和侧滑角测量技术领域,尤其涉及一种基于卡尔曼滤波的攻角、侧滑角辨识方法。
背景技术
目前,对于飞航导弹或无人机等飞行器而言,普遍不装备攻角、侧滑角传感器,因此缺乏直接测量攻角、侧滑角的手段。而在飞行控制、发动机控制中,飞行器的攻角、侧滑角又是一个重要的参控量。
目前,控制专业常使用过载替代攻角、侧滑角。但是,影响过载的因素很多,除攻侧角外,还有动压、飞行速度等因素,因此精度较差,可靠性低。
发明内容
本发明提供了一种基于卡尔曼滤波的攻角、侧滑角辨识方法,能够解决现有技术中在无攻角、侧滑角测量装置时无法精确测量攻角和侧滑角的技术问题。
本发明提供了一种基于卡尔曼滤波的攻角、侧滑角辨识方法,所述方法包括:
S1、建立攻角α与侧滑角β的递推模型;
S2、建立攻角α与侧滑角β的观测模型;
S3、基于攻角α和侧滑角β各自的卡尔曼滤波初值和所述递推模型,得到攻角α和侧滑角β的一步预测值Xpre;
S4、基于误差方差阵P和系统方差阵Q,得到一步预测均方差Ppre;
S5、基于所述观测模型,得到测量矩阵H和测量量Z;
S6、基于一步预测均方差Ppre、测量矩阵H和测量方差阵R,得到滤波增益K;
S7、基于一步预测值Xpre、滤波增益K、测量矩阵H、测量量Z、当前时刻的攻角α(k)和当前时刻的侧滑角β(k),得到卡尔曼滤波后的攻角α(k+1)和卡尔曼滤波后的侧滑角β(k+1);
S8、基于滤波增益K、测量矩阵H和一步预测均方差Ppre,得到更新后的误差方差阵P′;
S9、基于卡尔曼滤波后的攻角α(k+1)、卡尔曼滤波后的侧滑角β(k+1)和所述递推模型,得到攻角α和侧滑角β更新后的一步预测值Xpre′,并返回S4基于更新后的误差方差阵P′和更新后的一步预测值Xpre′进行下一轮迭代。
优选的,建立的攻角α与侧滑角β的递推模型为:
式中,α(k-1)为上一时刻的攻角,β(k-1)为上一时刻的侧滑角。
优选的,建立的攻角α与侧滑角β的观测模型为:
优选的,通过下式得到一步预测均方差Ppre:
Ppre=A·P·AT+Q,
式中,A为模型方阵。
优选的,通过下式得到滤波增益K:
K=Ppre·HT(H·Ppre·HT+R)-1,
优选的,通过下式得到卡尔曼滤波后的攻角α(k+1)和卡尔曼滤波后的侧滑角β(k+1):
优选的,通过下式得到更新后的误差方差阵P′:
P′=(I-KH)Ppre,
式中,I为单位矩阵。
应用本发明的技术方案,通过建立攻角α与侧滑角β的递推模型和观测模型,并基于攻角α与侧滑角β的递推模型和观测模型进行卡尔曼滤波,可实现攻角α和侧滑角β的高精度在线辨识。本发明在无攻角、侧滑角测量装置的条件下使用,节省了一套传感采集设备,而且可实现对攻角、侧滑角的同步测量。本发明可广泛运用于各种巡航类飞行器上,参与飞控算法解算,达到提升系统可靠性目的。
附图说明
所包括的附图用来提供对本发明实施例的进一步的理解,其构成了说明书的一部分,用于例示本发明的实施例,并与文字描述一起来阐释本发明的原理。显而易见地,下面描述中的附图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其他的附图。
图1示出了根据本发明一种实施例的基于卡尔曼滤波的攻角、侧滑角辨识方法的流程图;
图2示出了根据本发明一种实施例的攻角辨识仿真结果;
图3示出了根据本发明一种实施例的侧滑角辨识仿真结果。
具体实施方式
需要说明的是,在不冲突的情况下,本申请中的实施例及实施例中的特征可以相互组合。下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅仅是本发明一部分实施例,而不是全部的实施例。以下对至少一个示例性实施例的描述实际上仅仅是说明性的,决不作为对本发明及其应用或使用的任何限制。基于本发明中的实施例,本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
需要注意的是,这里所使用的术语仅是为了描述具体实施方式,而非意图限制根据本申请的示例性实施方式。如在这里所使用的,除非上下文另外明确指出,否则单数形式也意图包括复数形式,此外,还应当理解的是,当在本说明书中使用术语“包含”和/或“包括”时,其指明存在特征、步骤、操作、器件、组件和/或它们的组合。
除非另外具体说明,否则在这些实施例中阐述的部件和步骤的相对布置、数字表达式和数值不限制本发明的范围。同时,应当明白,为了便于描述,附图中所示出的各个部分的尺寸并不是按照实际的比例关系绘制的。对于相关领域普通技术人员已知的技术、方法和设备可能不作详细讨论,但在适当情况下,所述技术、方法和设备应当被视为授权说明书的一部分。在这里示出和讨论的所有示例中,任何具体值应被解释为仅仅是示例性的,而不是作为限制。因此,示例性实施例的其它示例可以具有不同的值。应注意到:相似的标号和字母在下面的附图中表示类似项,因此,一旦某一项在一个附图中被定义,则在随后的附图中不需要对其进行进一步讨论。
如图1所示,示出了一种基于卡尔曼滤波的攻角、侧滑角辨识方法,所述方法包括:
S1、建立攻角α与侧滑角β的递推模型;
S2、建立攻角α与侧滑角β的观测模型;
S3、基于攻角α和侧滑角β各自的卡尔曼滤波初值和所述递推模型,得到攻角α和侧滑角β的一步预测值Xpre;
S4、基于误差方差阵P和系统方差阵Q,得到一步预测均方差Ppre;
S5、基于所述观测模型,得到测量矩阵H和测量量Z;
S6、基于一步预测均方差Ppre、测量矩阵H和测量方差阵R,得到滤波增益K;
S7、基于一步预测值Xpre、滤波增益K、测量矩阵H、测量量Z、当前时刻的攻角α(k)和当前时刻的侧滑角β(k),得到卡尔曼滤波后的攻角α(k+1)和卡尔曼滤波后的侧滑角β(k+1);
S8、基于滤波增益K、测量矩阵H和一步预测均方差Ppre,得到更新后的误差方差阵P′;
S9、基于卡尔曼滤波后的攻角α(k+1)、卡尔曼滤波后的侧滑角β(k+1)和所述递推模型,得到攻角α和侧滑角β更新后的一步预测值Xpre′,并返回S4基于更新后的误差方差阵P′和更新后的一步预测值Xpre′进行下一轮迭代。
本发明通过建立攻角α与侧滑角β的递推模型和观测模型,并基于攻角α与侧滑角β的递推模型和观测模型进行卡尔曼滤波,可实现攻角α和侧滑角β的高精度在线辨识。本发明在无攻角、侧滑角测量装置的条件下使用,节省了一套传感采集设备,而且可实现对攻角、侧滑角的同步测量。本发明可广泛运用于各种巡航类飞行器上,参与飞控算法解算,达到提升系统可靠性目的。
根据本发明的一种实施例,在飞行动力学中,可以默认在极短时间内,攻角α与侧滑角β不会发生变化,从而建立攻角α与侧滑角β的递推模型为:
式中,α(k)为当前时刻的攻角,β(k)为当前时刻的侧滑角,α(k-1)为上一时刻的攻角,β(k-1)为上一时刻的侧滑角。
在工程应用中,常使用纵向过载ny和侧向过载nz来估计攻角α和侧滑角β,以用于飞行控制及发动机控制中。目前使用的观测模型如下式所示:
式中,ny为纵向过载,nz为侧向过载,k1为第一飞行相关系数,k2为第二飞行相关系数。
式中,Fy为纵向气动力,Fz为侧向气动力,Cyα为纵向气动力系数,Cyβ为侧向气动力系数,S为飞行器特征面积,q为动压。其中,纵向气动力系数Cyα和侧向气动力系数Cyβ由风洞试验吹风获得,飞行器特征面积S是飞行器的固有状态变量,动压q通过下式获得:
q=ρv2/2,
式中,ρ为空气密度,v为空速。其中,空气密度ρ依据标准大气并通过惯导测量的高度获得。由于很多飞行器未装备空速测量装置,因此在空速不可测的情况下,则可使用惯导测量的地速代替,从而获得近似动压q。
显然,受气动参数误差、实际空气密度与标准大气偏移、风场变化等因素影响,仅使用纵向过载ny和侧向过载nz来估计攻角α和侧滑角β存在较大误差。
在本发明中,攻角α和侧滑角β的观测模型增加了额外的约束,以提高攻角α和侧滑角β测量的精度和可靠性。
在本发明中,攻角α、侧滑角β与欧拉动力学关系如下式所示:
在飞行过程中,由于攻角α、侧滑角β、航向角ψ与航迹角ψv之差一般为小量,上述的攻角α、侧滑角β与欧拉动力学关系可简化为下式:
将上式中的第一中间变量G消除,可得到下式:
根据本发明的一种实施例,对于卡尔曼滤波,精确的初值可以减短收敛时间。因此,选取攻角α和侧滑角β各自的卡尔曼滤波初值为:
根据本发明的一种实施例,通过下式得到攻角α和侧滑角β的一步预测值Xpre:
式中,αpre为攻角α的一步预测值,βpre为侧滑角β的一步预测值,α(k)为当前时刻的攻角,β(k)为当前时刻的侧滑角。
根据本发明的一种实施例,通过下式得到一步预测均方差Ppre:
Ppre=A·P·AT+Q,
根据本发明的一种实施例,通过下式得到滤波增益K:
K=Ppre·HT(H·Ppre·HT+R)-1,
根据本发明的一种实施例,通过下式得到卡尔曼滤波后的攻角α(k+1)和卡尔曼滤波后的侧滑角β(k+1):
根据本发明的一种实施例,通过下式得到更新后的误差方差阵P′:
P′=(I-KH)Ppre,
下面以某型号飞行器为例,对本发明提供的基于卡尔曼滤波的攻角、侧滑角辨识方法进行具体说明。
第一步,建立的攻角α与侧滑角β的递推模型为:
建立的攻角α与侧滑角β的观测模型为:
第二步,选取攻角α与侧滑角β的观测模型中的第一飞行相关系数k1和第二飞行相关系数k2。
由飞行动力学可知,第一飞行相关系数k1和第二飞行相关系数k2与气动系数、特征面积、动压、飞行器质量及当前飞行重力加速度相关。其中重力加速度、特征面积数值较为固定;动压与气动系数随飞行状态变化较大;飞行器质量是可估计的换变量。因此,若想提高辨识精度,可以依据气动数据、大气模型、耗油率,制定随飞行时间、高度、速度变化的差值表。以本实施例中的飞行器为例,设定k1=1170,k2=-0.1。
第三步,选取卡尔曼滤波的三个参数,误差方差阵P、系统方差阵Q和测量方差阵R。
卡尔曼滤波参数关系卡尔曼滤波器性能。误差方差阵P、系统方差阵Q和测量方差阵R的选取与飞行器性能等多因素相关。以本实施例中的飞行器为例,卡尔曼滤波的三个参数分别为:
第四步,基于攻角α与侧滑角β的递推模型和观测模型进行卡尔曼滤波,并得到攻角α和侧滑角β的仿真结果,如图2、图3所示。由图2、图3可知,若飞行器处于遇到弹道机动阶段,攻角α与侧滑角β的辨识精度稍差,攻角α误差为0.29°,侧滑角β误差为0.26°;若飞行器处于巡航平飞阶段,攻角α与侧滑角β的辨识精度较高,攻角α误差为0.0033°,侧滑角β误差为0.0077°,可以忽略不计;若飞行器处于倾斜转弯(BTT)阶段,攻角α误差为0.026°,侧滑角β误差为0.02°。
综上所述,飞行器在稳定巡航飞行过程中,攻角α与侧滑角β的辨识误差均不超过0.1°,在弹道机动的过程中,攻角α与侧滑角β的辨识误差均不超过0.3°,满足测量的精度要求。
在本发明的描述中,需要理解的是,方位词如“前、后、上、下、左、右”、“横向、竖向、垂直、水平”和“顶、底”等所指示的方位或位置关系通常是基于附图所示的方位或位置关系,仅是为了便于描述本发明和简化描述,在未作相反说明的情况下,这些方位词并不指示和暗示所指的装置或元件必须具有特定的方位或者以特定的方位构造和操作,因此不能理解为对本发明保护范围的限制;方位词“内、外”是指相对于各部件本身的轮廓的内外。
为了便于描述,在这里可以使用空间相对术语,如“在……之上”、“在……上方”、“在……上表面”、“上面的”等,用来描述如在图中所示的一个器件或特征与其他器件或特征的空间位置关系。应当理解的是,空间相对术语旨在包含除了器件在图中所描述的方位之外的在使用或操作中的不同方位。例如,如果附图中的器件被倒置,则描述为“在其他器件或构造上方”或“在其他器件或构造之上”的器件之后将被定位为“在其他器件或构造下方”或“在其他器件或构造之下”。因而,示例性术语“在……上方”可以包括“在……上方”和“在……下方”两种方位。该器件也可以其他不同方式定位(旋转90度或处于其他方位),并且对这里所使用的空间相对描述作出相应解释。
此外,需要说明的是,使用“第一”、“第二”等词语来限定零部件,仅仅是为了便于对相应零部件进行区别,如没有另行声明,上述词语并没有特殊含义,因此不能理解为对本发明保护范围的限制。
以上所述仅为本发明的优选实施例而已,并不用于限制本发明,对于本领域的技术人员来说,本发明可以有各种更改和变化。凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。
Claims (7)
1.一种基于卡尔曼滤波的攻角、侧滑角辨识方法,其特征在于,所述方法包括:
S1、建立攻角α与侧滑角β的递推模型;
S2、建立攻角α与侧滑角β的观测模型;
S3、基于攻角α和侧滑角β各自的卡尔曼滤波初值和所述递推模型,得到攻角α和侧滑角β的一步预测值Xpre;
S4、基于误差方差阵P和系统方差阵Q,得到一步预测均方差Ppre;
S5、基于所述观测模型,得到测量矩阵H和测量量Z;
S6、基于一步预测均方差Ppre、测量矩阵H和测量方差阵R,得到滤波增益K;
S7、基于一步预测值Xpre、滤波增益K、测量矩阵H、测量量Z、当前时刻的攻角α(k)和当前时刻的侧滑角β(k),得到卡尔曼滤波后的攻角α(k+1)和卡尔曼滤波后的侧滑角β(k+1);
S8、基于滤波增益K、测量矩阵H和一步预测均方差Ppre,得到更新后的误差方差阵P′;
S9、基于卡尔曼滤波后的攻角α(k+1)、卡尔曼滤波后的侧滑角β(k+1)和所述递推模型,得到攻角α和侧滑角β更新后的一步预测值Xpre′,并返回S4基于更新后的误差方差阵P′和更新后的一步预测值Xpre′进行下一轮迭代。
4.根据权利要求3所述的一种基于卡尔曼滤波的攻角、侧滑角辨识方法,其特征在于,通过下式得到一步预测均方差Ppre:
Ppre=A·P·AT+Q,
式中,A为模型方阵。
7.根据权利要求3所述的一种基于卡尔曼滤波的攻角、侧滑角辨识方法,其特征在于,通过下式得到更新后的误差方差阵P′:
P′=(I-KH)Ppre,
式中,I为单位矩阵。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010242132.0A CN111412887B (zh) | 2020-03-31 | 2020-03-31 | 一种基于卡尔曼滤波的攻角、侧滑角辨识方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010242132.0A CN111412887B (zh) | 2020-03-31 | 2020-03-31 | 一种基于卡尔曼滤波的攻角、侧滑角辨识方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN111412887A true CN111412887A (zh) | 2020-07-14 |
CN111412887B CN111412887B (zh) | 2021-12-10 |
Family
ID=71491577
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202010242132.0A Active CN111412887B (zh) | 2020-03-31 | 2020-03-31 | 一种基于卡尔曼滤波的攻角、侧滑角辨识方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN111412887B (zh) |
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112146561A (zh) * | 2020-09-09 | 2020-12-29 | 无锡卡尔曼导航技术有限公司 | 一种霍尔角度传感器安装角度偏置的估计方法 |
CN112799417A (zh) * | 2020-12-29 | 2021-05-14 | 中国航空工业集团公司西安飞机设计研究所 | 一种大侧滑状态迎角信号修正方法 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102520726A (zh) * | 2011-12-19 | 2012-06-27 | 南京航空航天大学 | 大攻角飞行状态下的大气攻角及侧滑角估计方法 |
CN105005099A (zh) * | 2015-07-03 | 2015-10-28 | 南京航空航天大学 | 一种基于捷联惯导与飞行控制系统的大气参数解算方法 |
WO2019089966A1 (en) * | 2017-11-02 | 2019-05-09 | The Government Of The United States Of America, As Represented By The Secretary Of The Navy | Improved atmospheric thermal location estimation |
-
2020
- 2020-03-31 CN CN202010242132.0A patent/CN111412887B/zh active Active
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102520726A (zh) * | 2011-12-19 | 2012-06-27 | 南京航空航天大学 | 大攻角飞行状态下的大气攻角及侧滑角估计方法 |
CN105005099A (zh) * | 2015-07-03 | 2015-10-28 | 南京航空航天大学 | 一种基于捷联惯导与飞行控制系统的大气参数解算方法 |
WO2019089966A1 (en) * | 2017-11-02 | 2019-05-09 | The Government Of The United States Of America, As Represented By The Secretary Of The Navy | Improved atmospheric thermal location estimation |
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112146561A (zh) * | 2020-09-09 | 2020-12-29 | 无锡卡尔曼导航技术有限公司 | 一种霍尔角度传感器安装角度偏置的估计方法 |
CN112799417A (zh) * | 2020-12-29 | 2021-05-14 | 中国航空工业集团公司西安飞机设计研究所 | 一种大侧滑状态迎角信号修正方法 |
CN112799417B (zh) * | 2020-12-29 | 2023-03-14 | 中国航空工业集团公司西安飞机设计研究所 | 一种大侧滑状态迎角信号修正方法 |
Also Published As
Publication number | Publication date |
---|---|
CN111412887B (zh) | 2021-12-10 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
RU2755843C2 (ru) | Система и способ оценки воздушной скорости летательного аппарата на основании модели накопления данных о погоде. | |
US10520389B2 (en) | Aerodynamic modeling using flight data | |
EP3719447B1 (en) | Deep neural network-based inertial measurement unit (imu) sensor compensation method | |
CN111412887B (zh) | 一种基于卡尔曼滤波的攻角、侧滑角辨识方法 | |
CN107544074B (zh) | 一种无人机识别虚假gps信号的方法 | |
Borup et al. | A machine learning approach for estimating air data parameters of small fixed-wing UAVs using distributed pressure sensors | |
Wenz et al. | Moving horizon estimation of air data parameters for UAVs | |
CN112257186B (zh) | 针对小型四旋翼飞行器气动参数的时域辨识方法 | |
CN105865446A (zh) | 基于大气辅助的惯性高度通道阻尼卡尔曼滤波方法 | |
Borup et al. | Kalman filters for air data system bias correction for a fixed-wing UAV | |
Brossard et al. | Tightly coupled navigation and wind estimation for mini UAVs | |
CN114329766A (zh) | 面向深度强化学习的飞行动力学模型可信度评估方法 | |
US10876920B1 (en) | Auxiliary aerial vehicles for flow characterization | |
Koehl et al. | Wind-disturbance and aerodynamic parameter estimation of an experimental launched micro air vehicle using an EKF-like observer | |
CN107063248A (zh) | 基于旋翼转速的动力学模型辅助惯导的导航方法 | |
CN113916222A (zh) | 基于卡尔曼滤波估计方差约束的组合导航方法 | |
Jurado et al. | Complete online algorithm for air data system calibration | |
Ariante et al. | Velocity and attitude estimation of a small unmanned aircraft with micro Pitot tube and Inertial Measurement Unit (IMU) | |
Woodbury et al. | Flight test results of Observer/Kalman Filter Identifi [|# 12#|] cation of the Pegasus unmanned vehicle | |
Sun et al. | A gnss/imu-based 5-hole pitot tube calibration algorithm | |
WO2021140491A1 (en) | Method and system for estimating aerodynamic angles of a flying body | |
Siu et al. | Flight test results of an angle of attack and angle of sideslip calibration method using Output-Error optimization | |
CN106871892B (zh) | 一种航空器组合导航方法和装置 | |
McCrink et al. | Aerodynamic Parameter Estimation for Derived Angle-of-Attack Systems | |
CN117521561B (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 |