CN110716172B - 基于频率挑选的矢量水听器包络谱估计方法 - Google Patents

基于频率挑选的矢量水听器包络谱估计方法 Download PDF

Info

Publication number
CN110716172B
CN110716172B CN201911005544.6A CN201911005544A CN110716172B CN 110716172 B CN110716172 B CN 110716172B CN 201911005544 A CN201911005544 A CN 201911005544A CN 110716172 B CN110716172 B CN 110716172B
Authority
CN
China
Prior art keywords
signal
spectrum
frequency
self
signals
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
CN201911005544.6A
Other languages
English (en)
Other versions
CN110716172A (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.)
Harbin Engineering University
Original Assignee
Harbin Engineering University
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 Harbin Engineering University filed Critical Harbin Engineering University
Priority to CN201911005544.6A priority Critical patent/CN110716172B/zh
Publication of CN110716172A publication Critical patent/CN110716172A/zh
Application granted granted Critical
Publication of CN110716172B publication Critical patent/CN110716172B/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
    • G01S3/00—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received
    • G01S3/80—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received using ultrasonic, sonic or infrasonic waves
    • G01S3/86—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received using ultrasonic, sonic or infrasonic waves with means for eliminating undesired waves, e.g. disturbing noises
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01H—MEASUREMENT OF MECHANICAL VIBRATIONS OR ULTRASONIC, SONIC OR INFRASONIC WAVES
    • G01H3/00—Measuring characteristics of vibrations by using a detector in a fluid

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)

Abstract

基于频率挑选的矢量水听器包络谱估计方法，涉及舰船目标识别领域。本发明是为了解决舰船辐射噪声的包络谱提取不准的问题。本发明所述的本发明所述的基于频率挑选的矢量水听器包络谱估计方法，对船辐射噪声信号进行频带的优选，利用窗函数法构造有限冲激响应滤波器，从而得到高信噪比的目标信号进行包络谱分析，获得准确的包络谱。

Description

基于频率挑选的矢量水听器包络谱估计方法
技术领域
本发明属于舰船目标识别领域，尤其涉及舰船包络谱的估计方法。
背景技术
包络谱特征是舰船对目标进行识别的主要依据。现有技术中，舰船包络谱特征的提取方法主要有绝对值法、平方法等。包络谱特征提取的性能主要取决于分析频带内舰船目标辐射噪声的信噪比。因此，如何对船辐射噪声信号进行频带的优选、并设计相应滤波器得到高信噪比的目标信号进行包络谱分析成为获得清晰包络谱的关键。此外，当接收到的船辐射噪声存在其他方向的船噪声干扰时，包络谱特征提取性能急剧下降或者难以与目标对应。
针对上述问题，目前常用的一种解决办法就是从空域上对舰船目标辐射噪声频带进行优选，使其从噪声和其余目标中分离出来。但是，该方法在舰船辐射噪声的包络谱提取并不准确。
发明内容
本发明是为了解决舰船辐射噪声的包络谱提取不准的问题，现提供基于频率挑选的矢量水听器包络谱估计方法。
基于频率挑选的矢量水听器包络谱估计方法，包括以下步骤：
步骤一：利用矢量水听器在同一点采集待测目标辐射的声信号，所述声信号包括一路声压信号p(n)和两路振速信号vx(n)和vy(n)，0<n<Q，Q为每路信号的样本数，所述两路振速信号的方向位于同一水平面、且相互垂直；
步骤二：对两路振速信号vx(n)和vy(n)进行方向补偿，获得补偿后指向正北方向的振速信号
Figure GDA0003010778140000011
和补偿后指向正东方向的振速信号
Figure GDA0003010778140000012
步骤三：分别在声压信号p(n)和补偿后的两路振速信号
Figure GDA0003010778140000013
和
Figure GDA0003010778140000014
中截取一段包括N个采样点子信号、并分别进行N点傅里叶变换，获得声压信号的频谱值P(k)、两路振速信号的频谱值Vx(k)和Vy(k)，其中，N为每秒钟的信号采样数，N<Q，
利用P(k)分别对Vx(k)和Vy(k)求解互谱值和自谱值；
步骤四：利用步骤三获得的两个互谱值的实部求解不同信号频率k的方位θ(k)；
步骤五：根据不同信号频率k的方位θ(k)对各信号频率的自谱进行直方图统计，获得同一方位上所有信号频率的能量和，将最大能量和所在方位作为估计目标方位；
步骤六：从直方图中提取估计目标方位上信号频率的自谱值，根据该自谱值对应的待测目标的信号通频带、利用窗函数法构造有限冲激响应滤波器；
步骤七：利用步骤六获得的有限冲激响应滤波器对声压信号p(n)进行滤波、并对滤波结果取绝对值，利用韦尔奇法计算取绝对值后滤波结果的功率谱，将低于50Hz的低频段线谱作为目标信号的包络谱。
进一步的，步骤二中根据以下公式采用罗经对两路振速信号进行补偿，获得补偿后的两路振速信号
Figure GDA0003010778140000021
和
Figure GDA0003010778140000022
Figure GDA0003010778140000023
Figure GDA0003010778140000024
所述罗经与矢量水听器同步转动、且罗经的北向指向振速信号vx(n)的方向，φ(n)为罗经采集的偏转角度。
进一步的，步骤三中根据下式获得声压信号的频谱值P(k)：
Figure GDA0003010778140000025
根据下式分别获得两路振速信号的频谱值Vx(k)和Vy(k)：
Figure GDA0003010778140000026
Figure GDA0003010778140000027
其中，k为每秒采集的信号数、即信号频率，0≤k≤N-1。
进一步的，步骤三中根据下式分别求解Vx(k)和Vy(k)的互谱值和自谱值：
Figure GDA0003010778140000028
Figure GDA0003010778140000029
SPP(k)＝P(k)P(k)H 公式8
其中，
Figure GDA0003010778140000031
为p(n)与
Figure GDA0003010778140000032
的互谱值，
Figure GDA0003010778140000033
为p(n)与
Figure GDA0003010778140000034
的互谱值，SPP(k)为p(n)的自谱值。
进一步的，步骤四中利用下式获得不同信号频率k的方位θ(k)：
Figure GDA0003010778140000035
其中，k1<k<k2，
Figure GDA0003010778140000036
Figure GDA0003010778140000037
其中，fs为信号采样频率，fL和fH分别为工作频带的下限频率和上限频率。
进一步的，步骤五中对各信号频率的自谱进行直方图统计的具体方法为：
利用下式将同一方位上所有信号频率所对应的自谱值进行累加获得同一方位上所有信号频率的能量和S(α)：
S(α)＝∑Pα(k) 公式12
其中，Pα(k)表示α方位上信号频率k的自谱值，满足：
Figure GDA0003010778140000038
Δα为角度间隔；
所述估计目标方位
Figure GDA0003010778140000039
的表达式如下：
Figure GDA00030107781400000310
进一步的，步骤六中从直方图中提取估计目标方位
Figure GDA00030107781400000311
上信号频率k的自谱值
Figure GDA00030107781400000312
根据
Figure GDA00030107781400000313
采用窗函数法设计阶数为N的有限冲激响应滤波器，该滤波器的冲激响应为h(n)。
进一步的，步骤七中通过下式对声压信号p(n)进行滤波：
Figure GDA00030107781400000314
其中，pBP(n)为有限冲激响应滤波器的滤波结果，
Figure GDA00030107781400000315
表示卷积。
进一步的，步骤七中获得目标信号的包络谱的具体方法为：
利用窗函数法设计上限频率为50Hz、阶数为32的低通有限冲激响应滤波器，该滤波器的冲激响应为hLP(n)，
将滤波结果pBP(n)取绝对值输入上述低通有限冲激响应滤波器中，获得输出结果y(n)：
Figure GDA0003010778140000041
对y(n)进行N点傅里叶变换，获得变换结果Y(k)：
Figure GDA0003010778140000042
根据下式获得目标信号方向的包络谱E(k)：
Figure GDA0003010778140000043
本发明所述的基于频率挑选的矢量水听器包络谱估计方法，对船辐射噪声信号进行频带的优选，利用窗函数法构造有限冲激响应滤波器，从而得到高信噪比的目标信号进行包络谱分析，获得准确的包络谱。
附图说明
图1为本发明所述矢量水听器包络谱估计方法的流程图；
图2为矢量水听器x轴、y轴和罗经方向的示意图；
图3为矢量水听器方位估计流程图。
具体实施方式
具体实施方式一：参照图1至3具体说明本实施方式，本实施方式所述的基于频率挑选的矢量水听器包络谱估计方法，包括以下步骤：
步骤一：利用矢量水听器在同一点采集待测目标辐射的声信号，所述矢量水听器以潜标的方式锚泊在水下约40米深的位置，矢量水听器所在海域深度为80米；
所述声信号包括一路声压信号p(n)和两路振速信号vx(n)和vy(n)，0<n<Q，Q为每路信号的样本数，所述两路振速信号的方向位于同一水平面、且相互垂直。
步骤二：在矢量水听器的正下方安装罗经，所述罗经与矢量水听器同步转动、且罗经的北向指向振速信号vx(n)的方向，
利用罗经采集的偏转角度φ(n)分别对两路振速信号vx(n)和vy(n)进行方向补偿，获得补偿后指向正北方向的振速信号
Figure GDA0003010778140000051
和补偿后指向正东方向的振速信号
Figure GDA0003010778140000052
Figure GDA0003010778140000053
Figure GDA0003010778140000054
步骤三：分别在声压信号p(n)和补偿后的两路振速信号
Figure GDA0003010778140000055
和
Figure GDA0003010778140000056
中截取一段包括N个采样点的子信号、并分别进行N点傅里叶变换，获得声压信号的频谱值P(k)、两路振速信号的频谱值Vx(k)和Vy(k)，N为每秒钟的信号采样数，N<Q，
Figure GDA0003010778140000057
Figure GDA0003010778140000058
Figure GDA0003010778140000059
其中，k为每秒采集的信号数、即信号频率，0≤k≤N-1；
利用P(k)分别对Vx(k)和Vy(k)求解互谱值和自谱值：
Figure GDA00030107781400000510
Figure GDA00030107781400000511
SPP(k)＝P(k)P(k)H 公式8
其中，
Figure GDA00030107781400000512
为p(n)与
Figure GDA00030107781400000513
的互谱值，
Figure GDA00030107781400000514
为p(n)与
Figure GDA00030107781400000515
的互谱值，SPP(k)为p(n)的自谱值。
步骤四：利用步骤三获得的两个互谱值的实部求解不同信号频率k的方位θ(k)：
Figure GDA00030107781400000516
其中，k1<k<k2，
Figure GDA0003010778140000061
Figure GDA0003010778140000062
其中，fs为信号采样频率，fL和fH分别为工作频带的下限频率和上限频率。
步骤五：根据不同信号频率k的方位θ(k)对各信号频率的自谱进行直方图统计：
利用下式将同一方位上所有信号频率所对应的自谱值进行累加获得同一方位上所有信号频率的能量和S(α)：
S(α)＝∑Pα(k) 公式12
其中，Pα(k)表示α方位上信号频率k的自谱值，满足：
Figure GDA0003010778140000063
Δα为角度间隔；
将最大能量和所在方位作为估计目标方位
Figure GDA0003010778140000064
Figure GDA0003010778140000065
步骤六：从直方图中提取估计目标方位
Figure GDA0003010778140000066
上信号频率k的自谱值
Figure GDA0003010778140000067
根据
Figure GDA0003010778140000068
利用窗函数法构造阶数为N的有限冲激响应滤波器，该滤波器的冲激响应为h(n)。
步骤七：利用步骤六获得的有限冲激响应滤波器对声压信号p(n)进行滤波：
Figure GDA0003010778140000069
其中，pBP(n)为有限冲激响应滤波器的滤波结果，
Figure GDA00030107781400000610
表示卷积；
利用窗函数法设计上限频率为50Hz、阶数为32的低通有限冲激响应滤波器，该滤波器的冲激响应为hLP(n)，
将滤波结果pBP(n)取绝对值输入上述低通有限冲激响应滤波器中，获得输出结果y(n)：
Figure GDA00030107781400000611
对y(n)进行N点傅里叶变换，获得变换结果Y(k)：
Figure GDA0003010778140000071
根据下式获得目标信号方向的包络谱E(k)：
Figure GDA0003010778140000072

Claims (10)

1.基于频率挑选的矢量水听器包络谱估计方法，其特征在于，包括以下步骤：
步骤一：利用矢量水听器在同一点采集待测目标辐射的声信号，所述声信号包括一路声压信号p(n)和两路振速信号vx(n)和vy(n)，0<n<Q，Q为每路信号的样本数，所述两路振速信号的方向位于同一水平面、且相互垂直；
步骤二：对两路振速信号vx(n)和vy(n)进行方向补偿，获得补偿后指向正北方向的振速信号
Figure FDA0003010778130000011
和补偿后指向正东方向的振速信号
Figure FDA0003010778130000012
步骤三：分别在声压信号p(n)和补偿后的两路振速信号
Figure FDA0003010778130000013
和
Figure FDA0003010778130000014
中截取一段包括N个采样点的子信号、并分别进行N点傅里叶变换，获得声压信号的频谱值P(k)、两路振速信号的频谱值Vx(k)和Vy(k)，其中，N为每秒钟的信号采样数，N<Q，
利用P(k)分别对Vx(k)和Vy(k)求解互谱值和自谱值；
步骤四：利用步骤三获得的两个互谱值的实部求解不同信号频率k的方位θ(k)；
步骤五：根据不同信号频率k的方位θ(k)对各信号频率的自谱进行直方图统计，获得同一方位上所有信号频率的能量和，将最大能量和所在方位作为估计目标方位；
步骤六：从直方图中提取估计目标方位上信号频率的自谱值，根据该自谱值对应的待测目标的信号通频带、利用窗函数法构造有限冲激响应滤波器；
步骤七：利用步骤六获得的有限冲激响应滤波器对声压信号p(n)进行滤波、并对滤波结果取绝对值，利用韦尔奇法计算取绝对值后滤波结果的功率谱，将低于50Hz的低频段线谱作为目标信号的包络谱。
2.根据权利要求1所述的基于频率挑选的矢量水听器包络谱估计方法，其特征在于，利用矢量水听器采集声信号，所述矢量水听器以潜标的方式锚泊在水下约40米深的位置，矢量水听器所在海域深度为80米。
3.根据权利要求1所述的基于频率挑选的矢量水听器包络谱估计方法，其特征在于，步骤二中，根据以下公式采用罗经对两路振速信号进行补偿，获得补偿后的两路振速信号
Figure FDA0003010778130000015
和
Figure FDA0003010778130000016
Figure FDA0003010778130000017
Figure FDA0003010778130000018
所述罗经与矢量水听器同步转动、且罗经的北向指向振速信号vx(n)的方向，φ(n)为罗经采集的偏转角度。
4.根据权利要求1所述的基于频率挑选的矢量水听器包络谱估计方法，其特征在于，步骤三中，根据下式获得声压信号的频谱值P(k)：
Figure FDA0003010778130000021
根据下式分别获得两路振速信号的频谱值Vx(k)和Vy(k)：
Figure FDA0003010778130000022
Figure FDA0003010778130000023
其中，k为每秒采集的信号数、即信号频率，0≤k≤N-1。
5.根据权利要求4所述的基于频率挑选的矢量水听器包络谱估计方法，其特征在于，步骤三中，根据下式分别求解Vx(k)和Vy(k)的互谱值和自谱值：
Figure FDA0003010778130000024
Figure FDA0003010778130000025
SPP(k)＝P(k)P(k)H 公式8
其中，
Figure FDA00030107781300000214
为p(n)与
Figure FDA0003010778130000026
的互谱值，
Figure FDA0003010778130000027
为p(n)与
Figure FDA0003010778130000028
的互谱值，SPP(k)为p(n)的自谱值。
6.根据权利要求1所述的基于频率挑选的矢量水听器包络谱估计方法，其特征在于，步骤四中，利用下式获得不同信号频率k的方位θ(k)：
Figure FDA0003010778130000029
其中，k1<k<k2，
Figure FDA00030107781300000215
为p(n)与
Figure FDA00030107781300000210
的互谱值，
Figure FDA00030107781300000211
为p(n)与
Figure FDA00030107781300000212
的互谱值，
Figure FDA00030107781300000213
Figure FDA0003010778130000031
其中，fs为信号采样频率，fL和fH分别为工作频带的下限频率和上限频率。
7.根据权利要求1所述的基于频率挑选的矢量水听器包络谱估计方法，其特征在于，步骤五中，对各信号频率的自谱进行直方图统计的具体方法为：
利用下式将同一方位上所有信号频率所对应的自谱值进行累加获得同一方位上所有信号频率的能量和S(α)：
S(α)＝∑Pα(k) 公式12
其中，Pα(k)表示α方位上信号频率k的自谱值，满足：
Figure FDA0003010778130000032
Δα为角度间隔，SPP(k)为p(n)的自谱值；
所述估计目标方位
Figure FDA0003010778130000033
的表达式如下：
Figure FDA0003010778130000034
8.根据权利要求7所述的基于频率挑选的矢量水听器包络谱估计方法，其特征在于，步骤六中，从直方图中提取估计目标方位
Figure FDA0003010778130000035
上信号频率k的自谱值
Figure FDA0003010778130000036
根据
Figure FDA0003010778130000037
采用窗函数法设计阶数为N的有限冲激响应滤波器，该滤波器的冲激响应为h(n)。
9.根据权利要求8所述的基于频率挑选的矢量水听器包络谱估计方法，其特征在于，步骤七中，通过下式对声压信号p(n)进行滤波：
Figure FDA0003010778130000038
其中，pBP(n)为有限冲激响应滤波器的滤波结果，
Figure FDA0003010778130000039
表示卷积。
10.根据权利要求9所述的基于频率挑选的矢量水听器包络谱估计方法，其特征在于，步骤七中，获得目标信号的包络谱的具体方法为：
利用窗函数法设计上限频率为50Hz、阶数为32的低通有限冲激响应滤波器，该滤波器的冲激响应为hLP(n)，
将滤波结果pBP(n)取绝对值输入上述低通有限冲激响应滤波器中，获得输出结果y(n)：
Figure FDA0003010778130000041
对y(n)进行N点傅里叶变换，获得变换结果Y(k)：
Figure FDA0003010778130000042
根据下式获得目标信号方向的包络谱E(k)：
Figure FDA0003010778130000043
CN201911005544.6A 2019-10-22 2019-10-22 基于频率挑选的矢量水听器包络谱估计方法 Active CN110716172B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911005544.6A CN110716172B (zh) 2019-10-22 2019-10-22 基于频率挑选的矢量水听器包络谱估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911005544.6A CN110716172B (zh) 2019-10-22 2019-10-22 基于频率挑选的矢量水听器包络谱估计方法

Publications (2)

Publication Number Publication Date
CN110716172A CN110716172A (zh) 2020-01-21
CN110716172B true CN110716172B (zh) 2021-07-09

Family

ID=69213067

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911005544.6A Active CN110716172B (zh) 2019-10-22 2019-10-22 基于频率挑选的矢量水听器包络谱估计方法

Country Status (1)

Country Link
CN (1) CN110716172B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111735525B (zh) * 2020-05-28 2023-03-31 哈尔滨工程大学 一种适用于无人声纳的demon谱特征提取方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2001051039A (ja) * 1999-08-09 2001-02-23 Nec Corp 信号処理方式
CN101271155A (zh) * 2008-04-25 2008-09-24 哈尔滨工程大学 用浅海低频相干矢量声场干涉结构探测并判别目标的方法
CN101539597B (zh) * 2009-04-29 2011-01-26 哈尔滨工程大学 一种分辨具有相同频带辐射噪声的多目标方法
CN104237871B (zh) * 2013-06-08 2017-01-11 中国科学院声学研究所 一种基于相位补偿的时延差估计方法
CN105223576B (zh) * 2015-09-24 2017-11-21 哈尔滨工程大学 一种基于单矢量潜标的线谱信号目标自动检测方法
CN105182345B (zh) * 2015-09-26 2017-10-03 哈尔滨工程大学 基于单矢量潜标的连续谱信号目标自动检测方法

Also Published As

Publication number Publication date
CN110716172A (zh) 2020-01-21

Similar Documents

Publication Publication Date Title
US20080270033A1 (en) Methods of hydrocarbon detection using spectral energy analysis
CN112285767B (zh) 海底地震仪四分量海洋面波多阶频散能量成像装置及方法
CN103592628A (zh) 一种基于球谐域实值权重波束形成的多声源定位方法
CN111024209B (zh) 一种适用于矢量水听器的线谱检测方法
CN101539597A (zh) 一种分辨具有相同频带辐射噪声的多目标方法
CN102590856A (zh) 基于小波频谱分析的位场异常分离方法
CN110716172B (zh) 基于频率挑选的矢量水听器包络谱估计方法
CN111856401A (zh) 一种基于互谱相位拟合的时延估计方法
CN102323618B (zh) 基于分数阶傅里叶变换的相干噪声抑制方法
CN103916199B (zh) 一种天线信号的时延和相位调整装置与方法
CN102323616B (zh) 提高灰岩出露区地震数据分辨率的分频匹配方法
CN115015942B (zh) 一种自适应水下目标声学测速装置及方法
CN110703198B (zh) 基于频率挑选的四元十字阵包络谱估计方法
CN115236592B (zh) 一种基于单阵元时频曲线匹配的冰声定位方法
CN102338884B (zh) 物探中的椭圆窗方向带通保幅滤波数据处理方法
GB2104218A (en) Detecting harmonically-rich acoustic sources
CN102338890B (zh) 一种地球物理勘探中的圆形窗带通保幅滤波数据处理方法
CN107957571A (zh) 水听器测向方法、装置、计算机可读存储介质及计算机设备
CN112464855A (zh) 基于eemd的随钻泥浆正脉冲信号处理方法及装置
CN109682459B (zh) 一种水下管口外辐射噪声的测量装置和方法
CN102368094B (zh) 一种地球物理勘探中的矩形窗带通保幅滤波数据处理方法
CN112162235A (zh) 一种光滑分段随机共振增强的声矢量信号定向方法
CN103344990A (zh) 基于matlab电磁式可控震源输出信号的分析方法
CN109358361B (zh) 一种提取地震射线矢量的方法
CN117169816B (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