CN102445203A - 一种刚体空间运动状态的埃米特输出方法 - Google Patents

一种刚体空间运动状态的埃米特输出方法 Download PDF

Info

Publication number
CN102445203A
CN102445203A CN201110280617XA CN201110280617A CN102445203A CN 102445203 A CN102445203 A CN 102445203A CN 201110280617X A CN201110280617X A CN 201110280617XA CN 201110280617 A CN201110280617 A CN 201110280617A CN 102445203 A CN102445203 A CN 102445203A
Authority
CN
China
Prior art keywords
centerdot
phi
rigid body
state
output
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
Application number
CN201110280617XA
Other languages
English (en)
Other versions
CN102445203B (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.)
Xian Feisida Automation Engineering Co Ltd
Original Assignee
Xian Feisida Automation Engineering Co Ltd
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 Xian Feisida Automation Engineering Co Ltd filed Critical Xian Feisida Automation Engineering Co Ltd
Priority to CN 201110280617 priority Critical patent/CN102445203B/zh
Publication of CN102445203A publication Critical patent/CN102445203A/zh
Application granted granted Critical
Publication of CN102445203B publication Critical patent/CN102445203B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Landscapes

  • Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)

Abstract

本发明公开了一种刚体空间运动状态的埃米特输出方法,该方法通过定义三元数,使得机体轴系三个速度分量和三元数构成线性微分方程组,并采用埃米特正交多项式对滚转、俯仰、偏航角速度p,q,r进行近似逼近描述,可以按照任意阶保持器的方式求解系统的状态转移矩阵,进而得到刚体运动离散状态方程的表达式,避免了姿态方程奇异问题,从而得到刚体主要运动状态;本发明通过引入三元数使得状态转移矩阵为分块上三角形式,可以降阶求解状态转移矩阵,大大简化了计算复杂度,便于工程使用。

Description

一种刚体空间运动状态的埃米特输出方法
技术领域
本发明涉及空间运动刚体模型,特别涉及飞行器大机动飞行状态输出问题。
背景技术
机体轴系刚体运动微分方程是描述飞行器、鱼雷、航天器等空间运动的基本方程。通常,在数据处理等应用中,体轴系的状态变量主要包含3个速度分量、三个欧拉角、以及地面坐标系的XE,YE,ZE等,由于ZE定义为垂直地面指向地球中心,因此ZE实际为负的飞行高度;XE,YE通常主要依赖GPS、GNSS、北斗等直接给出;欧拉角表示刚体空间运动姿态,而刻画刚体姿态的微分方程又是其中的核心,通常以三个欧拉角即俯仰、滚转和偏航角来描述。当刚体的俯仰角为±90°时,滚转角和偏航角无法定值,同时临近该奇点的区域求解误差过大,导致工程上不可容忍的误差而不能使用;为了避免这一问题,人们首先采用限制俯仰角取值范围的方法,这使得方程式退化,不能全姿态工作,因而难以广泛用于工程实践。随着对飞行器极限飞行的研究,人们又相继采用了方向余弦法、等效转动矢量法、四元数法等推算刚体运动姿态。
方向余弦法避免了欧拉角描述方法的“奇异”现象,用方向余弦法计算姿态矩阵没有方程退化问题,可以全姿态工作,但需要求解9个微分方程,计算量较大,实时性较差,无法满足工程实践要求。等效转动矢量法如单子样递推、双子样转动矢量、三子样转动矢量和四子样旋转矢量法以及在此基础上的各种修正算法和递推算法等。文献中研究旋转矢量时,都是基于速率陀螺输出为角增量的算法。然而在实际工程中,一些陀螺的输出是角速率信号,如光纤陀螺、动力调谐陀螺等。当速率陀螺输出为角速率信号时,旋转矢量法的算法误差明显增大。四元数法是定义4个欧拉角的函数来计算航姿,能够有效弥补欧拉角描述方法的奇异性,只要解4个一阶微分方程式组即可,比方向余弦姿态矩阵微分方程式计算量有明显的减少,能满足工程实践中对实时性的要求。其常用的计算方法有毕卡逼近法、二阶、四阶龙格-库塔法和三阶泰勒展开法等。毕卡逼近法实质是单子样算法,对有限转动引起的不可交换误差没有补偿,在高动态情况下姿态解算中的算法漂移会十分严重。采用四阶龙格-库塔法求解四元数微分方程时,随着积分误差的不断积累,会出现三角函数取值超出±1的现象,从而导致计算发散;泰勒展开法也因计算精度的不足而受到制约。当刚体大机动时,角速率较大导致上述方法的误差更大;不仅如此,姿态估计的误差常常会导致速度4个分量、高度输出的误差急剧增大。
发明内容
为了克服现有刚体运动模型输出误差大的问题,本发明提供一种刚体空间运动状态的埃米特输出方法,该方法通过定义三元数,使得机体轴系三个速度分量和三元数构成线性微分方程组,并采用埃米特正交多项式对滚转、俯仰、偏航角速度p,q,r进行近似逼近描述,可以按照任意阶保持器的方式求解系统的状态转移矩阵,进而得到刚体运动离散状态方程的表达式,避免了姿态方程奇异问题,从而得到刚体主要运动状态。
本发明解决其技术问题采用的技术方案是,一种刚体空间运动状态的埃米特输出方法,其特征包括以下步骤:
1、机体轴系三个速度分量输出为:
u ( t ) v ( t ) w ( t ) t = ( k + 1 ) T = Φ v [ ( k + 1 ) T , kT ] u ( t ) v ( t ) w ( t ) t = kT + g Φ v [ ( k + 1 ) T , kT ] Φ s [ ( k + 1 ) T , kT ] s 1 ( t ) s 2 ( t ) s 3 ( t ) t = kT
+ g ∫ kT ( k + 1 ) T Φ v [ ( k + 1 ) T , τ ] n x n y n z dτ
其中:u,v,w分别为沿刚体机体轴系x,y,z轴的速度分量,nx,ny,nz分别为沿x,y,z轴的过载,g为重力加速度,s1、s2、s3为定义的三元数,且
s 1 ( t ) s 2 ( t ) s 3 ( t ) t = ( k + 1 ) T = Φ s [ ( k + 1 ) T , kT ] s 1 ( t ) s 2 ( t ) s 3 ( t ) t = kT
Φ v [ ( k + 1 ) T , kT ] ≈ I + Π v Hξ ( t ) | kT ( k + 1 ) T + Π v ∫ kT ( k + 1 ) T [ ξ ( t ) ξ T ( t ) ] dt H T Π v T - Π v Hξ ( t ) | kT ( k + 1 ) T Π v Hξ ( kT )
Φ s [ ( k + 1 ) T , kT ] ≈ I + Π s Hξ ( t ) | kT ( k + 1 ) T + Π s ∫ kT ( k + 1 ) T [ ξ ( t ) ξ T ( t ) ] dt H T Π s T - Π s Hξ ( t ) | kT ( k + 1 ) T Π s Hξ ( kT )
p,q,r分别为滚转、俯仰、偏航角速度,T为采样周期;全文参数定义相同;
I = 1 0 0 0 1 0 0 0 1 , ξ(t)=[ξ0(t)ξ1(t)…ξn-1(t) ξn(t)]T
ξ 0 ( t ) = 1 ξ 1 ( t ) = 2 t ξ 2 ( t ) = 4 t 2 - 2 ξ 3 ( t ) = 8 t 3 - 12 t ξ 4 ( t ) = 16 t 4 - 48 t 2 + 12 ξ 5 ( t ) = 32 t 5 - 160 t 3 + 120 t ξ 6 ( t ) = 64 t 6 - 480 t 4 + 720 t - 120 · · · ξ i + 1 ( t ) = 2 t ξ i ( t ) - 2 i ξ i - 1 ( t ) i = 2,3 , · · · , n - 1
为埃米特正交多项式的递推形式,滚转、俯仰、偏航角速度p,q,r的展开式分别为
p(t)=[p0 p1…pn-1 pn][ξ0(t)ξ1(t)…ξn-1(t)ξn(t)]T
q(t)=[q0 q1…qn-1 qn][ξ0(t)ξ1(t)…ξn-1(t)ξn(t)]T
r(t)=[r0 r1…rn-1 rn][ξ0(t)ξ1(t)…ξn-1(t)ξn(t)]T
Π v = 0 0 0 0 0 1 0 - 1 0 p 0 p 1 · · · p n - 1 p n
+ 0 0 - 1 0 0 0 1 0 0 q 0 q 1 · · · q n - 1 q n + 0 1 0 - 1 0 0 0 0 0 r 0 r 1 · · · r n - 1 r n
Π s = 0 0 0 0 0 1 0 - 1 0 p 0 p 1 · · · p n - 1 p n
+ 0 0 1 0 0 0 - 1 0 0 q 0 q 1 · · · q n - 1 q n + 0 - 1 0 1 0 0 0 0 0 r 0 r 1 · · · r n - 1 r n
Figure BSA00000577641100035
h12=0.5,h23=0.25, h i ( i + 1 ) = h ( i - 1 ) i 1 + 2 h ( i - 1 ) i , · · ·
h21=0.5,h41=-1.5, h ( 2 i ) 1 = 4 ih ( 2 i - 2 ) 1 1 + 2 h ( 2 i ) ( 2 i + 1 ) , · · ·
其余hij=0;
2、高度输出为:
h · = u v w s 1 s 2 s 3
其中:h为高度;
3、姿态角的输出为:
Figure BSA00000577641100041
Figure BSA00000577641100042
ψ ( t ) = ψ ( kT ) + ∫ kT t qs 2 ( t ) + rs 3 ( t ) s 2 2 ( t ) + s 3 2 ( t ) dt
其中:
Figure BSA00000577641100044
θ,ψ分别表示滚转、俯仰、偏航角, s 1 ( t ) s 2 ( t ) s 3 ( t ) = Φ s ( t , kT ) s 1 ( t ) s 2 ( t ) s 3 ( t ) t = kT .
本发明的有益效果是:通过引入三元数使得状态转移矩阵为分块上三角形式,可以降阶求解状态转移矩阵,大大简化了计算复杂度,便于工程使用。
下面结合实施例对本发明作详细说明。
具体实施方式
1、机体轴系三个速度分量输出为:
u ( t ) v ( t ) w ( t ) t = ( k + 1 ) T = Φ v [ ( k + 1 ) T , kT ] u ( t ) v ( t ) w ( t ) t = kT + g Φ v [ ( k + 1 ) T , kT ] Φ s [ ( k + 1 ) T , kT ] s 1 ( t ) s 2 ( t ) s 3 ( t ) t = kT
+ g ∫ kT ( k + 1 ) T Φ v [ ( k + 1 ) T , τ ] n x n y n z dτ
其中:u,v,w分别为沿刚体机体轴系x,y,z轴的速度分量,nx,ny,nz分别为沿x,y,z轴的过载,g为重力加速度,s1、s2、s3为定义的三元数,且
s 1 ( t ) s 2 ( t ) s 3 ( t ) t = ( k + 1 ) T = Φ s [ ( k + 1 ) T , kT ] s 1 ( t ) s 2 ( t ) s 3 ( t ) t = kT
Φ v [ ( k + 1 ) T , kT ] ≈ I + Π v Hξ ( t ) | kT ( k + 1 ) T + Π v ∫ kT ( k + 1 ) T [ ξ ( t ) ξ T ( t ) ] dt H T Π v T - Π v Hξ ( t ) | kT ( k + 1 ) T Π v Hξ ( kT )
Φ s [ ( k + 1 ) T , kT ] ≈ I + Π s Hξ ( t ) | kT ( k + 1 ) T + Π s ∫ kT ( k + 1 ) T [ ξ ( t ) ξ T ( t ) ] dt H T Π s T - Π s Hξ ( t ) | kT ( k + 1 ) T Π s Hξ ( kT )
p,q,r分别为滚转、俯仰、偏航角速度,T为采样周期;全文参数定义相同;
I = 1 0 0 0 1 0 0 0 1 , ξ(t)=[ξ0(t)ξ1(t)…ξn-1(t)ξn(t)]T
ξ 0 ( t ) = 1 ξ 1 ( t ) = 2 t ξ 2 ( t ) = 4 t 2 - 2 ξ 3 ( t ) = 8 t 3 - 12 t ξ 4 ( t ) = 16 t 4 - 48 t 2 + 12 ξ 5 ( t ) = 32 t 5 - 160 t 3 + 120 t ξ 6 ( t ) = 64 t 6 - 480 t 4 + 720 t - 120 · · · ξ i + 1 ( t ) = 2 t ξ i ( t ) - 2 i ξ i - 1 ( t ) i = 2,3 , · · · , n - 1
为埃米特正交多项式的递推形式,滚转、俯仰、偏航角速度p,q,r的展开式分别为
p(t)=[p0 p1…pn-1 pn][ξ0(t)ξ1(t)…ξn-1(t)ξn(t)]T
q(t)=[q0 q1…qn-1 qn][ξ0(t)ξ1(t)…ξn-1(t)ξn(t)]T
r(t)=[r0 r1…rn-1 rn][ξ0(t)ξ1(t)…ξn-1(t)ξn(t)]T
Π v = 0 0 0 0 0 1 0 - 1 0 p 0 p 1 · · · p n - 1 p n
+ 0 0 - 1 0 0 0 1 0 0 q 0 q 1 · · · q n - 1 q n + 0 1 0 - 1 0 0 0 0 0 r 0 r 1 · · · r n - 1 r n
Π s = 0 0 0 0 0 1 0 - 1 0 p 0 p 1 · · · p n - 1 p n
+ 0 0 1 0 0 0 - 1 0 0 q 0 q 1 · · · q n - 1 q n + 0 - 1 0 1 0 0 0 0 0 r 0 r 1 · · · r n - 1 r n
Figure BSA00000577641100061
h12=0.5,h23=0.25, h i ( i + 1 ) = h ( i - 1 ) i 1 + 2 h ( i - 1 ) i , · · ·
h21=0.5,h41=-1.5, h ( 2 i ) 1 = 4 ih ( 2 i - 2 ) 1 1 + 2 h ( 2 i ) ( 2 i + 1 ) , · · ·
其余hij=0;
2、高度输出为:
h · = u v w s 1 s 2 s 3
其中:h为高度;
3、姿态角的输出为:。
Figure BSA00000577641100066
ψ [ ( k + 1 ) T ] = ψ ( kT ) + ∫ kT ( k + 1 ) T qs 2 ( t ) + rs 3 ( t ) s 2 2 ( t ) + s 3 2 ( t ) dt
其中:
Figure BSA00000577641100068
θ,ψ分别表示滚转、俯仰、偏航角, s 1 ( t ) s 2 ( t ) s 3 ( t ) t = ( k + 1 ) T = Φ s [ ( k + 1 ) T , kT ) s 1 ( t ) s 2 ( t ) s 3 ( t ) t = kT .

Claims (1)

1.一种刚体空间运动状态的埃米特输出方法,其特征包括以下步骤:
机体轴系三个速度分量输出为:
u ( t ) v ( t ) w ( t ) t = ( k + 1 ) T = Φ v [ ( k + 1 ) T , kT ] u ( t ) v ( t ) w ( t ) t = kT + g Φ v [ ( k + 1 ) T , kT ] Φ s [ ( k + 1 ) T , kT ] s 1 ( t ) s 2 ( t ) s 3 ( t ) t = kT
+ g ∫ kT ( k + 1 ) T Φ v [ ( k + 1 ) T , τ ] n x n y n z dτ
其中:u,v,w分别为沿刚体机体轴系x,y,z轴的速度分量,nx,ny,nz分别为沿x,y,z轴的过载,g为重力加速度,s1、s2、s3为定义的三元数,且
s 1 ( t ) s 2 ( t ) s 3 ( t ) t = ( k + 1 ) T = Φ s [ ( k + 1 ) T , kT ] s 1 ( t ) s 2 ( t ) s 3 ( t ) t = kT
Φ v [ ( k + 1 ) T , kT ] ≈ I + Π v Hξ ( t ) | kT ( k + 1 ) T + Π v ∫ kT ( k + 1 ) T [ ξ ( t ) ξ T ( t ) ] dt H T Π v T - Π v Hξ ( t ) | kT ( k + 1 ) T Π v Hξ ( kT )
Φ s [ ( k + 1 ) T , kT ] ≈ I + Π s Hξ ( t ) | kT ( k + 1 ) T + Π s ∫ kT ( k + 1 ) T [ ξ ( t ) ξ T ( t ) ] dt H T Π s T - Π s Hξ ( t ) | kT ( k + 1 ) T Π s Hξ ( kT )
p,q,r分别为滚转、俯仰、偏航角速度,T为采样周期;全文参数定义相同;
I = 1 0 0 0 1 0 0 0 1 , ξ(t)=[ξ0(t)ξ1(t)…ξn-1(t)ξn(t)T
ξ 0 ( t ) = 1 ξ 1 ( t ) = 2 t ξ 2 ( t ) = 4 t 2 - 2 ξ 3 ( t ) = 8 t 3 - 12 t ξ 4 ( t ) = 16 t 4 - 48 t 2 + 12 ξ 5 ( t ) = 32 t 5 - 160 t 3 + 120 t ξ 6 ( t ) = 64 t 6 - 480 t 4 + 720 t - 120 · · · ξ i + 1 ( t ) = 2 t ξ i ( t ) - 2 i ξ i - 1 ( t ) i = 2,3 , · · · , n - 1
为埃米特正交多项式的递推形式,滚转、俯仰、偏航角速度p,q,r的展开式分别为
p(t)=[p0 p1…pn-1 pn][ξ0(t)ξ1(t)…ξn-1(t)ξn(t)]T
q(t)=[q0 q1…qn-1 qn][ξ0(t)ξ1(t)…ξn-1(t)ξn(t)]T
r(t)=[r0 r1…rn-1 rn][ξ0(t)ξ1(t)…ξn-1(t)ξn(t)]T
Π v = 0 0 0 0 0 1 0 - 1 0 p 0 p 1 · · · p n - 1 p n
+ 0 0 - 1 0 0 0 1 0 0 q 0 q 1 · · · q n - 1 q n + 0 1 0 - 1 0 0 0 0 0 r 0 r 1 · · · r n - 1 r n
Π s = 0 0 0 0 0 1 0 - 1 0 p 0 p 1 · · · p n - 1 p n
+ 0 0 1 0 0 0 - 1 0 0 q 0 q 1 · · · q n - 1 q n + 0 - 1 0 1 0 0 0 0 0 r 0 r 1 · · · r n - 1 r n
h12=0.5,h23=0.25, h i ( i + 1 ) = h ( i - 1 ) i 1 + 2 h ( i - 1 ) i , · · ·
h21=0.5,h41=-1.5, h ( 2 i ) 1 = 4 ih ( 2 i - 2 ) 1 1 + 2 h ( 2 i ) ( 2 i + 1 ) , · · ·
其余hij=0;
高度输出为:
h · = u v w s 1 s 2 s 3
其中:h为高度;
姿态角的输出为:
Figure FSA00000577641000029
Figure FSA000005776410000210
ψ ( t ) = ψ ( kT ) + ∫ kT t qs 2 ( t ) + rs 3 ( t ) s 2 2 ( t ) + s 3 2 ( t ) dt
其中:
Figure FSA00000577641000031
θ,ψ分别表示滚转、俯仰、偏航角, s 1 ( t ) s 2 ( t ) s 3 ( t ) = Φ s ( t , kT ) s 1 ( t ) s 2 ( t ) s 3 ( t ) t = kT .
CN 201110280617 2011-09-20 2011-09-20 一种刚体空间运动状态的埃米特输出方法 Expired - Fee Related CN102445203B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN 201110280617 CN102445203B (zh) 2011-09-20 2011-09-20 一种刚体空间运动状态的埃米特输出方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN 201110280617 CN102445203B (zh) 2011-09-20 2011-09-20 一种刚体空间运动状态的埃米特输出方法

Publications (2)

Publication Number Publication Date
CN102445203A true CN102445203A (zh) 2012-05-09
CN102445203B CN102445203B (zh) 2013-05-08

Family

ID=46007929

Family Applications (1)

Application Number Title Priority Date Filing Date
CN 201110280617 Expired - Fee Related CN102445203B (zh) 2011-09-20 2011-09-20 一种刚体空间运动状态的埃米特输出方法

Country Status (1)

Country Link
CN (1) CN102445203B (zh)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1422380A (zh) * 2000-02-03 2003-06-04 独立技术有限责任公司 用变更四元数数据表示对可倾斜物体中的方位角估算
US20070213889A1 (en) * 2004-02-27 2007-09-13 Instituto Nacional De Tecnica Aeroespacial "Esteba Terradas" Sensor Fusion System and Method for Estimating Position, Speed and Orientation of a Vehicle, in Particular an Aircraft
US20110190964A1 (en) * 2010-01-29 2011-08-04 Eurocopter Turning-stabilized estimation of the attitude angles of an aircraft

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1422380A (zh) * 2000-02-03 2003-06-04 独立技术有限责任公司 用变更四元数数据表示对可倾斜物体中的方位角估算
US20070213889A1 (en) * 2004-02-27 2007-09-13 Instituto Nacional De Tecnica Aeroespacial "Esteba Terradas" Sensor Fusion System and Method for Estimating Position, Speed and Orientation of a Vehicle, in Particular an Aircraft
US20110190964A1 (en) * 2010-01-29 2011-08-04 Eurocopter Turning-stabilized estimation of the attitude angles of an aircraft

Also Published As

Publication number Publication date
CN102445203B (zh) 2013-05-08

Similar Documents

Publication Publication Date Title
CN102589553B (zh) 建立飞行器运动模型的切换方法
CN102436437B (zh) 基于角速度的飞行器极限飞行时四元数傅里埃近似输出方法
CN102445202B (zh) 一种刚体空间运动状态的拉盖尔输出方法
CN102384747A (zh) 一种刚体空间运动状态的Hartley输出方法
CN102359789B (zh) 一种刚体空间运动状态的任意阶输出方法
CN102323990A (zh) 一种刚体空间运动的气动模型
CN102346729B (zh) 一种刚体空间运动状态的勒让德输出方法
CN102384746A (zh) 一种刚体空间运动状态的切比雪夫输出方法
CN102508818A (zh) 一种刚体空间运动状态的任意步长正交级数输出方法
CN102323992A (zh) 一种刚体空间运动状态的多项式类输出方法
CN102508819B (zh) 基于角速度的飞行器极限飞行时四元数勒让德近似输出方法
CN102445203B (zh) 一种刚体空间运动状态的埃米特输出方法
CN102495831B (zh) 基于角速度的飞行器极限飞行时四元数埃米特近似输出方法
CN102495830B (zh) 基于角速度的飞行器极限飞行时四元数Hartley近似输出方法
CN102495829B (zh) 基于角速度的飞行器极限飞行时四元数沃尔什近似输出方法
CN102359790B (zh) 一种刚体空间运动状态的傅里埃输出方法
CN102323991B (zh) 一种刚体空间运动状态沃尔什输出模型的建模方法
CN102508821A (zh) 一种刚体空间运动的状态输出方法
CN102506866B (zh) 基于角速度的飞行器极限飞行时四元数切比雪夫近似输出方法
CN102375803B (zh) 一种刚体空间运动的气流轴系模型的建立方法
CN102346727A (zh) 一种刚体空间运动的机体轴系模型
CN102506864B (zh) 飞行器极限飞行时四元数任意步长正交级数近似输出方法
CN102506865B (zh) 基于角速度的飞行器极限飞行时四元数多项式类近似输出方法
CN102495825B (zh) 基于角速度的飞行器极限飞行时四元数超线性输出方法
CN102494688B (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
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20130508