CN112906482B - 基于随机共振与混沌协同的井下微弱特征信号检测方法 - Google Patents

基于随机共振与混沌协同的井下微弱特征信号检测方法 Download PDF

Info

Publication number
CN112906482B
CN112906482B CN202110095126.1A CN202110095126A CN112906482B CN 112906482 B CN112906482 B CN 112906482B CN 202110095126 A CN202110095126 A CN 202110095126A CN 112906482 B CN112906482 B CN 112906482B
Authority
CN
China
Prior art keywords
signal
stochastic resonance
duffing
chaos
bistable
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
CN202110095126.1A
Other languages
English (en)
Other versions
CN112906482A (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 Shiyou University
Original Assignee
Xian Shiyou 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 Xian Shiyou University filed Critical Xian Shiyou University
Priority to CN202110095126.1A priority Critical patent/CN112906482B/zh
Publication of CN112906482A publication Critical patent/CN112906482A/zh
Application granted granted Critical
Publication of CN112906482B publication Critical patent/CN112906482B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/08Feature extraction
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/11Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems

Abstract

本发明公开了基于随机共振与混沌协同的井下微弱特征信号检测方法,涉及信号检测技术领域,以双稳态Duffing系统为理论框架,利用其随机共振与混沌的协同作用获取信号完整的参数信息。不仅回避了混沌检测对待测信号频率的匹配要求;而且弥补了随机共振无法测定信号“大小”的不足。利用双稳态Duffing系统产生的随机共振与混沌效应,有效识别与噪声频谱混叠的微弱特征信号,并大幅提高待测信号的信噪比,从而为随钻动态测量弱信号检测提供解决途径。

Description

基于随机共振与混沌协同的井下微弱特征信号检测方法
技术领域
本发明涉及信号检测技术领域,特别涉及基于随机共振与混沌协同的井下微弱特征信号检测方法。
背景技术
受钻具振动、井下干扰磁场及高温高压等因素的影响,加速度计、磁通门和陀螺仪等传感器装置均会产生巨大的漂移误差,导致原始测量信号被严重干扰。因此,简单将两种或多种传感器机械组合,钻具姿态参数仍然存在不准确甚至不可测的问题。
FIR滤波是一种基于频域分析的噪声抑制方法,在处理随钻强振动噪声这种时变非平稳信号时,会导致部分有用信号的丢失;并且在抑制低频噪声时存在过度去噪、损害特征信号的问题。小波变换可以同时分析信号的时频信息,但该方法存在小波基参数选择难度大、选取依赖主观性等缺点;并且可处理信号的信噪比较高,无法满足随钻测量极低信噪比的要求。自适应滤波在进行参数控制时,算法的稳定性和收敛性无法保证,难以满足姿态测量信号处理对可靠性和实时性的要求。综上所述,基于噪声抑制或抵消的滤波方法在实际应用时,势必会受到特征信号信噪比及噪声频谱范围的限制。
在钻井工程中,由于底部钻具的强烈振动和快速旋转导致姿态测量信号中含有多频高幅值的干扰信号,造成有用信号幅度微弱并且信噪比极低。目前,随钻动态测量系统主要是借鉴惯性测量技术,将加速度计传感器与磁性传感器装置或陀螺仪联合,组成多传感器测量系统。常用的井下去噪方法包括FIR滤波、小波变换和自适应滤波等。现有的井下信号检测方法大都基于噪声抑制机制,在滤除噪声的同时,有用信号也不可避免的受到抑制或损害;甚至与噪声信号频谱混叠的有用信号很难被识别,导致姿态测量无效。
针对现有技术存在的问题,本申请提供了一种基于随机共振与混沌协同的井下微弱特征信号检测方法,以双稳态Duffing系统为理论框架,利用其产生的随机共振与混沌特性协同作用获取信号完整的参数信息。不仅回避了混沌检测对待测信号频率的匹配要求,而且弥补了随机共振无法测定信号“大小”的不足。
发明内容
本发明的目的在于提供一种基于随机共振与混沌协同的井下微弱特征信号检测方法,利用双稳态Duffing系统产生的随机共振与混沌效应,有效识别与噪声频谱混叠的微弱特征信号,并大幅提高待测信号的信噪比,从而为随钻动态测量弱信号检测提供解决途径。
本发明提供了一种基于随机共振与混沌协同的井下微弱特征信号检测方法,包括以下步骤:
步骤1:建立双稳态Duffing系统随机共振模型,如下所示:
Figure GDA0003639004580000021
其中,k是阻尼比,-ax+bx3是非线性恢复力,a和b是大于零的实系数,s(t)是特征信号;n(t)为背景噪声信号。
步骤2:对双稳态Duffing系统随机共振模型进行尺度变换;
步骤3:将尺度变换后的姿态测量信号输入双稳态Duffing系统随机共振模型,使系统输出处于混沌至大尺度周期的临界态;
步骤4:将两个双稳态Duffing系统随机共振模型的输出误差作为自适应控制器的输入,并根据两个双稳态Duffing系统随机共振模型方程建立控制目标;
步骤5:构造目标函数,即Lyapunov函数,如下所示:
Figure GDA0003639004580000031
步骤6:根据基于最速梯度下降的自适应控制算法,计算目标函数随系统输出变化的梯度,进而得到关于λ、
Figure GDA0003639004580000032
的两个伪梯度函数;
步骤7:对伪梯度函数进行迭代计算,直至控制目标收敛,得出系统稳定后的控制量u是钻具弱信号在Lyapunov稳定条件下关于幅值和相位的最优估计。
进一步地,所述步骤S1中公式(1)的特征信号为
Figure GDA0003639004580000033
其中f是测量信号的实际频率值,引入变尺度系数R,且R=100·f0,t’=Rt,此时,微弱特征信号表示为:
s(t’)=λcos(2πf0t)=λcos(2πf0·t’/R)=λcos(2π·0.01·t’) (3)。
进一步地,所述步骤2对数值计算的步长进行尺度变换的步骤如下:
步骤201:假设采样频率为fs,则数值计算的步长为dt=1/fs
步骤202:引入变尺度系数R,数值计算步长变为:
dt’=Rdt=R/fs (4);
步骤203:将信号输入双稳态Duffing系统随机共振模型,并令该信号的步长满足式(4),完成对随钻测量信号的频率重构。
进一步地,所述步骤4中的控制目标如下:
Figure GDA0003639004580000034
当t→∞ (7)
令控制量
u=-Ke1+3xxge11sinωt+θ2cosωt (9)
其中,K是反馈增益,Ke1
Figure GDA0003639004580000035
的线性近似;θ1、θ2分别作为
Figure GDA0003639004580000036
Figure GDA0003639004580000037
的估计;λ、
Figure GDA0003639004580000038
分别表示待恢复信号的幅值和初相角。
与现有技术相比,本发明具有如下显著优点:
本发明提出的一种基于随机共振与混沌协同的井下微弱特征信号检测方法,以双稳态Duffing系统为理论框架,利用其随机共振与混沌的协同作用获取信号完整的参数信息。不仅回避了混沌检测对待测信号频率的匹配要求;而且弥补了随机共振无法测定信号“大小”的不足。利用双稳态Duffing系统产生的随机共振与混沌效应,有效识别与噪声频谱混叠的微弱特征信号,并大幅提高待测信号的信噪比,从而为随钻动态测量弱信号检测提供解决途径。
附图说明
图1为本发明实施例提供的基于双耦合Duffing混沌系统的弱信号参数同步控制图;
图2为本发明实施例提供的输出信号信噪比随输入噪声的变化规律图;
图3为本发明实施例提供的混沌系统相态跃变图;
图4为本发明实施例提供的共振输出效应图。
具体实施方式
下面结合本发明中的附图,对本发明实施例的技术方案进行清楚、完整的描述,显然,所描述的实施例是本发明的一部分实施例,而不是全部实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例,都应属于本发明保护的范围。
目前,随钻动态测量系统主要是借鉴惯性测量技术,将加速度计传感器与磁性传感器装置或陀螺仪联合,组成多传感器测量系统。常用的井下去噪方法包括FIR滤波、小波变换和自适应滤波等。具体性能比较如表1所示。
表1常用去噪方式性能比较
Figure GDA0003639004580000041
Figure GDA0003639004580000051
近年来,非线性学科的不断发展,为微弱信号检测提供了新的思路,其中随机共振和混沌是两种备受关注的检测方法。
随机共振是指在输入信号和噪声的协同作用下,利用非线性系统使输出信号的信噪比随噪声的增大呈现先增后减的趋势,并在某一噪声强度下出现峰值,产生共振输出现象,如图2所示。研究表明,基于双稳态Duffing系统的随机共振模型可以有效识别强噪声背景下的微弱周期信号,例如随钻测量时加速度计与磁性传感器所输出的正弦周期信号。
然而,非线性系统在合适参数条件下所产生的随机共振效应,可以有效提高系统输出信号的信噪比,但不能对微弱特征信号的“大小”进行测定,即无法检测信号的幅值和相位参数。因而不能获取姿态测量信号的完整参数信息。
利用混沌相变对系统参数摄动的敏感性以及对噪声的免疫能力,可以在较少的测量数据和任意有色噪声背景下实现极低的检测信噪比。在前期研究中发现:将姿态测量传感器输出的微弱周期信号输入处于混沌临界状态的双稳态Duffing系统,系统输出相态会发生根本变化,如图3所示,依据这种变化可以实现随钻微弱特征信号的检测。
目前,将双稳态Duffing系统混沌特性应用于随钻动态测量微弱信号检测也具有一定的局限性:受近钻头强振动影响,井下姿态测量信号的频率是以钻具转速为基准不断波动变化,而非固定不变,因此不能直接将其作为特征信号与系统驱动项进行线性叠加,进而无法通过混沌相变对特征信号进行识别。
参照图1-4,本发明提供了一种基于随机共振与混沌协同的井下微弱特征信号检测方法,包括以下步骤:
步骤1:建立双稳态Duffing系统随机共振模型,如下所示:
Figure GDA0003639004580000061
其中,k是阻尼比,取0.5;-ax+bx3是非线性恢复力,a和b是大于零的实系数,取值1;s(t)是特征信号;n(t)为背景噪声信号。
如图4所示,在合适的参数条件下,双稳态Duffing系统在适当的参数条件下可以产生共振输出效应。对比输入和输出信号的频谱可以看出,在特征信号频率,即f0=0.01Hz处,信号幅值明显加强,为典型的随机共振效应。然而,该系统只能识别小参数信号的频率值(f0=0.01Hz),应用于井下姿态测量信号时严重受限。因此,针对传统随机共振理论的不足,本申请提出基于尺度变换Duffing系统的随机共振频率检测方法。
步骤2:为了使姿态测量信号满足随机共振对小频率参数的限制条件,对双稳态Duffing系统随机共振模型进行尺度变换;
步骤3:将尺度变换后的姿态测量信号输入双稳态Duffing系统随机共振模型,使系统输出处于混沌至大尺度周期的临界态;
步骤4:将两个双稳态Duffing系统随机共振模型的输出误差作为自适应控制器的输入,并根据两个双稳态Duffing系统随机共振模型方程建立控制目标;
步骤5:构造目标函数,即Lyapunov函数,如下所示:
Figure GDA0003639004580000062
步骤6:根据基于最速梯度下降的自适应控制算法,计算目标函数随系统输出变化的梯度,进而得到关于λ、
Figure GDA0003639004580000063
的两个伪梯度函数;
步骤7:对伪梯度函数进行迭代计算,直至控制目标收敛,得出系统稳定后的控制量u是钻具弱信号在Lyapunov稳定条件下关于幅值和相位的最优估计。
参照图1,其中,基于Duffing系统混沌相变的幅值、相位同步估计方案,从而获得完整的姿态测量信息。为此,设计基于双耦合Duffing混沌系统的弱信号参数同步控制方案,该控制方案包含两个参数相同的Duffing系统,且被同一内置信号驱动,处于混沌到大尺度周期的临界态。
实施例1
假设所述步骤S1中公式(1)的特征信号为
Figure GDA0003639004580000071
其中f是测量信号的实际频率值,保持Duffing系统其余参数条件不变,引入变尺度系数R,使s(t)在时间轴上放大R倍,即t’=Rt。此时,在数值上使变步长系数R=100·f0,则测量信号变为:
Figure GDA0003639004580000072
从公式(3)可以看出,随钻姿态测量信号的频率从f(约1~3Hz)被压缩为0.01Hz或接近于0.01Hz的值。此时,双稳态Duffing模型符合小参数条件,可以应用随机共振特性检测信号s(t’)是否存在。
所述步骤2对数值计算的步长进行尺度变换的步骤如下:
步骤201:假设采样频率为fs,则数值计算的步长为dt=1/fs
步骤202:引入变尺度系数R,数值计算步长变为:
dt’=Rdt=R/fs (4);
此时将信号的时间间隔增大了R倍,相应的信号频率被压缩了R倍。即一个采样频率为fs、特征频率为f的测量信号,通过步长变换变为采样频率为fs/R、特征频率为f/R的信号;
步骤203:将信号输入双稳态Duffing系统随机共振模型,并令该信号的步长满足式(4),完成对随钻测量信号的频率重构。进而通过公式(1)识别出微弱信号的频率值。
实施例2
对于图1中的Duffing主系统,根据公式(1)可建立如下系统方程:
Figure GDA0003639004580000081
其中,Acosωt是驱动信号,
Figure GDA0003639004580000082
是经由实施例1确定频率后的待测信号。
对于图1中的Duffing受控系统,建立系统方程如下:
Figure GDA0003639004580000083
所述步骤4中的控制目标如下:
Figure GDA0003639004580000084
将式(6)与(5)相减,可得
Figure GDA0003639004580000085
其中,e1=x-xd,e2=y-yd。令控制量
u=-Ke1+3xxde11sinωt+θ2cosωt (9)
其中,K是反馈增益,通常取大于1的常数;Ke1
Figure GDA0003639004580000086
的线性近似;θ1、θ2分别作为
Figure GDA0003639004580000087
Figure GDA0003639004580000088
的估计,λ、
Figure GDA0003639004580000089
分别表示待恢复信号的幅值和初相角。
为实现式(7)给出的控制目标,建立步骤5中的Lyapunov函数:
Figure GDA00036390045800000810
当K>1时,可得
Figure GDA00036390045800000811
公式(10)说明,选择公式(8)作为目标函数可以满足使控制系统稳定的可达性条件。因此,令控制目标函数Q(e)=V(e)。选择伪梯度函数如下:
Figure GDA00036390045800000812
由(11)可得快速梯度自适应算法如下:
Figure GDA0003639004580000091
其中,β、δ是自适应增益,取大于零的常数。
公式(9)和(12)满足快速梯度算法的正则、凸性、可达性以及有界性条件,能够实现控制目标(7),即式(6)的受控系统与式(5)的主系统达到同步。此时,可以认为输出控制量u是微弱待测信号
Figure GDA0003639004580000092
在Lyapunov稳定条件下的最优估计。
以上公开的仅为本发明的几个具体实施例,但是,本发明实施例并非局限于此,任何本领域的技术人员能思之的变化都应落入本发明的保护范围。

Claims (4)

1.一种基于随机共振与混沌协同的井下微弱特征信号检测方法,其特征在于,包括以下步骤:
步骤1:建立双稳态Duffing系统随机共振模型,如下所示:
Figure FDA0003667072500000011
其中,k是阻尼比,-ax+bx3是非线性恢复力,a和b是大于零的实系数,s(t)是特征信号,n(t)为背景噪声信号;
步骤2:对双稳态Duffing系统随机共振模型进行尺度变换;
步骤3:将尺度变换后的姿态测量信号输入双稳态Duffing系统随机共振模型,使系统输出处于混沌至大尺度周期的临界态;
步骤4:将两个双稳态Duffing系统随机共振模型的输出误差作为自适应控制器的输入,并根据两个双稳态Duffing系统随机共振模型方程建立控制目标;
其中,双稳态Duffing系统随机共振模型的输出误差作为自适应控制器的输入包括:
对于Duffing主系统,根据公式(1)可建立系统方程:
Figure FDA0003667072500000012
其中,Acosωt是驱动信号,
Figure FDA0003667072500000013
是待测信号;
对于Duffing受控系统,建立系统方程:
Figure FDA0003667072500000014
将(3)和(2)式相减,可得
Figure FDA0003667072500000015
其中,e1=x-xd,e2=y-yd,xd表示Duffing主系统中的变量x,u=-Ke1+3xxde11sin ωt+θ2cos ωt为控制量;
步骤5:构造目标函数,即Lyapunov函数,如下所示:
Figure FDA0003667072500000021
其中,V(e)表示Lyapunov函数,K表示反馈增益,取大于1的常数;
步骤6:根据基于最速梯度下降的自适应控制算法,计算目标函数随系统输出变化的梯度,进而得到关于λ、
Figure FDA0003667072500000022
的两个伪梯度函数,其中,伪梯度函数包括:
Figure FDA0003667072500000023
其中,β是自适应增益,通常取大于零的常数;
步骤7:对伪梯度函数进行迭代计算,直至控制目标收敛,得出系统稳定后的控制量u是钻具弱信号在Lyapunov稳定条件下关于幅值和相位的最优估计。
2.如权利要求1所述的一种基于随机共振与混沌协同的井下微弱特征信号检测方法,其特征在于,所述步骤1中公式(1)的特征信号为
Figure FDA0003667072500000024
其中f0是测量信号的实际频率值,引入变尺度系数R,且R=100·f0,t’=Rt,此时,微弱特征信号表示为:
Figure FDA0003667072500000025
3.如权利要求2所述的一种基于随机共振与混沌协同的井下微弱特征信号检测方法,其特征在于,所述步骤2对数值计算的步长进行尺度变换的步骤如下:
步骤201:假设采样频率为fs,则数值计算的步长为dt=1/fs
步骤202:引入变尺度系数R,数值计算步长变为:
dt’=Rdt=R/fs (7);
步骤203:将信号输入双稳态Duffing系统随机共振模型,并令该信号的步长满足式(4),完成对随钻测量信号的频率重构。
4.如权利要求1所述的一种基于随机共振与混沌协同的井下微弱特征信号检测方法,其特征在于,所述步骤4中的控制目标如下:
Figure FDA0003667072500000031
令控制量
u=-Ke1+3xxge11sinωt+θ2cosωt (9)
其中,K是反馈增益,Ke1
Figure FDA0003667072500000032
的线性近似;θ1、θ2分别作为
Figure FDA0003667072500000033
Figure FDA0003667072500000034
的估计;λ、
Figure FDA0003667072500000035
分别表示待恢复信号的幅值和初相角。
CN202110095126.1A 2021-01-25 2021-01-25 基于随机共振与混沌协同的井下微弱特征信号检测方法 Active CN112906482B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110095126.1A CN112906482B (zh) 2021-01-25 2021-01-25 基于随机共振与混沌协同的井下微弱特征信号检测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110095126.1A CN112906482B (zh) 2021-01-25 2021-01-25 基于随机共振与混沌协同的井下微弱特征信号检测方法

Publications (2)

Publication Number Publication Date
CN112906482A CN112906482A (zh) 2021-06-04
CN112906482B true CN112906482B (zh) 2022-07-29

Family

ID=76117618

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110095126.1A Active CN112906482B (zh) 2021-01-25 2021-01-25 基于随机共振与混沌协同的井下微弱特征信号检测方法

Country Status (1)

Country Link
CN (1) CN112906482B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113608021B (zh) * 2021-07-23 2022-05-20 华中科技大学 一种混沌振子电路和基于混沌理论的微弱信号检测系统

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104734824A (zh) * 2015-03-31 2015-06-24 华南理工大学 强噪声背景下混沌同步检测微弱周期信号的方法
JP2017184101A (ja) * 2016-03-31 2017-10-05 キヤノン株式会社 信号抽出処理装置および信号抽出処理方法
CN111125613A (zh) * 2019-12-06 2020-05-08 江苏大学 一种提高Duffing混沌振子检测微弱谐振信号抗噪声能力的方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104734824A (zh) * 2015-03-31 2015-06-24 华南理工大学 强噪声背景下混沌同步检测微弱周期信号的方法
JP2017184101A (ja) * 2016-03-31 2017-10-05 キヤノン株式会社 信号抽出処理装置および信号抽出処理方法
CN111125613A (zh) * 2019-12-06 2020-05-08 江苏大学 一种提高Duffing混沌振子检测微弱谐振信号抗噪声能力的方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
随钻弱SNR信号的Duffing振子混沌检测与恢复;杨一等;《仪器仪表学报》;20200324;第41卷(第02期);全文 *

Also Published As

Publication number Publication date
CN112906482A (zh) 2021-06-04

Similar Documents

Publication Publication Date Title
Zheng et al. Real-time dynamic displacement monitoring with double integration of acceleration based on recursive least squares method
CN105247432B (zh) 频率响应测定装置
CN107688554B (zh) 基于自适应傅里叶分解的频域辨识方法
CN104453857B (zh) 一种小井斜下井斜和工具面角动态测量方法及装置
WO2013047787A1 (ja) 加速度センサ回路
CN112906482B (zh) 基于随机共振与混沌协同的井下微弱特征信号检测方法
Henry et al. Prism signal processing for sensor condition monitoring
CN104516991A (zh) 一种伽马传感器全温度范围补偿方法
CN108680247A (zh) 基于振动烈度低频滤波修正的振动信号转换方法
Flammini et al. Least mean square method for LVDT signal processing
Li et al. Complex bandpass filtering for Coriolis mass flow meter signal processing
CN104500036B (zh) 一种适用于小井斜条件下井斜工具面角动态测量的方法及装置
Tyagi et al. Frequency estimation techniques in capacitance-to-frequency conversion measurement
CN111551785B (zh) 基于无迹卡尔曼滤波的频率与谐波检测方法
Thenozhi et al. A novel numerical integrator for structural health monitoring
CN109033646B (zh) 基于小量程mems加速度传感器的随钻测斜仪重力工具面角计算方法
Hewitson et al. Calibration of the power-recycled gravitational wave detector, GEO 600
CN107664499B (zh) 一种船用捷联惯导系统的加速度计在线降噪方法
Li et al. Modified wavelet filtering algorithm applied to gyro servo technology for the improvement of test-precision
Wielandt et al. Measuring seismometer nonlinearity on a shake table
Grecheneva et al. Compensation of the accelerometer errors in solving the problem of kinematic control of dynamic objects
Yang et al. Empirical analysis of the MIMU random-drift error based on wavelet soft-threshold and ARIMA model
CN112149291B (zh) 微弱谐波信号检测系统及方法
Cardona-Morales et al. Instantaneous frequency estimation and order tracking based on Kalman filters
US20220236308A1 (en) Method for improving emc robustness of integrated capacitive sensors

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