CN108226996B - 基于能量频带分布的自适应各向异性分频分区滤波方法 - Google Patents

基于能量频带分布的自适应各向异性分频分区滤波方法 Download PDF

Info

Publication number
CN108226996B
CN108226996B CN201711453771.6A CN201711453771A CN108226996B CN 108226996 B CN108226996 B CN 108226996B CN 201711453771 A CN201711453771 A CN 201711453771A CN 108226996 B CN108226996 B CN 108226996B
Authority
CN
China
Prior art keywords
imf
frequency
seismic data
filtering
adaptive
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
CN201711453771.6A
Other languages
English (en)
Other versions
CN108226996A (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.)
Chengdu Univeristy of Technology
Original Assignee
Chengdu Univeristy of 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 Chengdu Univeristy of Technology filed Critical Chengdu Univeristy of Technology
Priority to CN201711453771.6A priority Critical patent/CN108226996B/zh
Publication of CN108226996A publication Critical patent/CN108226996A/zh
Application granted granted Critical
Publication of CN108226996B publication Critical patent/CN108226996B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/364—Seismic filtering
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/20—Trace signal pre-filtering to select, remove or transform specific events or signal components, i.e. trace-in/trace-out
    • G01V2210/21—Frequency-domain filtering, e.g. band pass
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/30—Noise handling
    • G01V2210/32—Noise reduction
    • G01V2210/324—Filtering

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种基于能量频带分布的自适应各向异性分频分区滤波方法，包括以下步骤：(1)输入二维叠后地震数据u；(2)利用VMD对二维叠后地震数据进行频域分解，得到不同频率范围的IMF剖面uk(k＝1,2,...,n)，n为整数；(3)将分解后的每个IMF剖面uk(k＝1,2,...,n)分别经自适应阈值各向异性滤波算法进行多次迭代处理，每次迭代处理后得到一重构的最终结果
Figure DDA0001528901430000011
(4)将多次迭代得到的多个最终结果
Figure DDA0001528901430000012
分别计算其信噪比SNR和相似度SSIM；(5)选取最优的SNR和SSIM对应的最终结果进行输出。本发明结合信号在频率域和时间域的特征进行去噪，更好的在滤除噪声的同时保护地震资料纹理的局部特征和主要结构信息，同时提高了地震资料的质量。

Description

基于能量频带分布的自适应各向异性分频分区滤波方法
技术领域
本发明涉及一种地震资料分析处理方法，尤其涉及一种基于能量频带分布的自适应各向异性分频分区滤波方法。
背景技术
野外采集的地震资料中包含着与地下结构和岩性有关的有用信息，而地下结构的不连续性和地层边缘信息往往预示着该处的油气藏可能性。但在地震资料采集的过程中，往往都会受到外界噪声的干扰，使得边界和断层信息被扭曲甚至掩盖，导致地震资料的质量降低，不宜直接用来进行地质解释，因此需要先对地震数据进行滤波处理。在不断的研究中发现，噪声在不同的频率范围内分布的强度是不同的，传统的全频带单一处理方法如Chambolle和Lions提出的CL模型，它容易造成某些频率范围内频率的过分压制，使得一些有效信息被滤除(比如边界、断层等信息)。因此，分频滤波成了目前研究的热点。
目前主流的信号分解方法主要有经验模态分解(简称EMD)，它能很好的对非线性和非稳态的数据进行自适应和多尺度分析。但是EMD是一种经验模态分解法，缺乏一定的数学理论基础支撑，并且存在模态混叠现象。集合经验模态分解(Ensemble Empirical ModeDecomposition，简称EEMD)和完备集经验模态分解(Complete Ensemble Empirical ModeDecomposition，简称CEEMD)都对EMD的模态混叠现象进行了改进，但两者存在模态冗余问题，并且在分解信号时容易受到强噪声的干扰，使分解的模态不稳定。针对上述信号分解存在的问题，Dragomiretskiy和Zosso提出了变分模态分解(Variational ModeDecomposition，简称VMD)。VMD不仅很好的避免了模态混叠和冗余的现象，同时对高频噪声有一定的抑制作用，适合于对多尺度地震信号进行分析。
传统的CL模型采用一个常量梯度阈值β(β＞0)来区分高低特征区域。但在一些复杂情况下，β的选取比较困难，并且针对多尺度噪声，全局采用一个梯度阈值进行选择容易造成过度滤波导致边界模糊等情况，影响后期地震资料的解释。基于偏微分方程的各向异性扩散滤波的改进方法很多，基于阈值改进的也有，但都没有考虑到噪声在不同频率范围的分布强度是不同的，依然是在全频带进行处理。
发明内容
本发明的目的就在于提供一种解决上述问题，克服了CL模型阈值选取困难以及在全频带范围内处理的不足的问题，而且提高了各向异性滤波在多尺度噪声环境下的去噪和边界恢复能力的基于能量频带分布的自适应各向异性分频分区滤波方法。
为了实现上述目的，本发明采用的技术方案是这样的：一种基于能量频带分布的自适应各向异性分频分区滤波方法，包括以下步骤：
(1)输入二维叠后地震数据u；
(2)利用VMD对二维叠后地震数据进行频域分解，得到不同频率范围的IMF剖面uk(k＝1,2,...,n)，n为整数且为IMF剖面总数；
(3)将分解后的每个IMF剖面uk(k＝1,2,...,n)分别经自适应阈值各向异性滤波算法进行多次迭代处理，每次迭代处理后得到一重构的最终结果
Figure BDA00015289014100000315
具体步骤如下：
(31)先计算第一个IMF剖面的自适应阈值
Figure BDA00015289014100000316
梯度值
Figure BDA0001528901410000031
和保真项εc；
(32)根据下式将梯度值
Figure BDA0001528901410000032
与自适应阈值
Figure BDA00015289014100000317
相比较，当
Figure BDA0001528901410000033
时，采取L2-norm滤波；当
Figure BDA0001528901410000034
时，采用TV方法滤波，得到滤波后的IMF剖面
Figure BDA0001528901410000035
Figure BDA0001528901410000036
(33)再以
Figure BDA0001528901410000037
为基础，重复步骤(31)(32)迭代m次，得到滤波后的IMF剖面
Figure BDA0001528901410000038
j为迭代次数；
(34)重复步骤(31)-(33)，得到每个IMF剖面对应的
Figure BDA0001528901410000039
(35)将每个IMF剖面中迭代次数相同的ukj进行叠加重构，得到重构后的最终结果
Figure BDA00015289014100000310
(4)将多次迭代得到的多个最终结果
Figure BDA00015289014100000311
分别计算其信噪比SNR和相似度SSIM；
(5)判定信噪比SNR和相似度SSIM，选取最优的SNR和SSIM对应的最终结果
Figure BDA00015289014100000312
作为最终结果
Figure BDA00015289014100000313
进行输出。
作为优选：步骤(2)中，VMD对二维叠后地震数据进行频域分解的方法为，
根据下式计算uk(k＝1,2,...,n)
Figure BDA00015289014100000314
其中，δ(t)是狄拉克函数，ωk表示每个频带范围的中心频率。
作为优选：步骤(31)中自适应阈值
Figure BDA0001528901410000041
梯度值
Figure BDA0001528901410000042
和保真项εc的计算方法如下：
Figure BDA0001528901410000043
Figure BDA0001528901410000044
式中，
Figure BDA0001528901410000045
为梯度算子，λ,μ,ρ∈(0,1)，所述ρ由IMF剖面的中心频率决定，中心频率越高，ρ值越大。
作为优选：步骤(4)中计算信噪比SNR和相似度SSIM的方法为：
Figure BDA0001528901410000046
Figure BDA0001528901410000047
式中，μi是平均强度，δi表示标准偏差，这里的i∈x,y，δxy表示图像x和y的协方差，Cj＝(KjL)2，其中L是像素值的自适应范围，Kj＜＜1为常数。一般K选用很小的常数。
与现有技术相比，本发明的优点在于：它将图形去噪方法运用到地震资料去噪中，结合基于CL模型的阈值选取问题进行改进得到自适应阈值各向异性滤波方法，考虑信号在频率域与时空域的局部特征，针对不同IMF剖面进行自适应阈值调整，同时结合不同频带内的能量强度比不同，区分不同强度噪声和特征属性。该方法不仅克服了CL模型阈值选取困难以及在全频带范围内处理的不足的问题，而且提高了各向异性滤波在多尺度噪声环境下的去噪和边界恢复能力。具体如下：
1.结合信号在频率域和时间域的特征进行去噪；
2.每个IMF剖面蕴含了对应频率范围纹理的局部特征和主要结构信息，根据其特征进行自适应阈值分区去噪，更好的在滤除噪声的同时保护纹理的局部特征和主要结构信息，同时提高了地震资料的质量。
附图说明
图1为本发明流程图；
图2为合成的地堑模型剖面；
图3为图2加入随机高斯噪声后的模型剖面图；
图4为将图3经为本发明方法处理后的剖面图；
图5为局部原始实际资料剖面图；
图6为图5经为本发明方法处理后的剖面图。
具体实施方式
下面将结合附图对本发明作进一步说明。
实施例1：参见图1，一种基于能量频带分布的自适应各向异性分频分区滤波方法，包括以下步骤：
(1)输入二维叠后地震数据u；
(2)利用VMD对二维叠后地震数据进行频域分解，得到不同频率范围的IMF剖面uk(k＝1,2,...,n)，n为整数；步骤(2)中，VMD对二维叠后地震数据进行频域分解的方法为，
根据下式计算uk(k＝1,2,...,n)
Figure BDA0001528901410000061
其中，δ(t)是狄拉克函数，ωk表示每个频带范围的中心频率；
(3)将分解后的每个IMF剖面uk(k＝1,2,...,n)分别经自适应阈值各向异性滤波算法进行多次迭代处理，每次迭代处理后得到一重构的最终结果u；
具体步骤如下：
(31)先计算第一个IMF剖面的自适应阈值
Figure BDA0001528901410000062
梯度值
Figure BDA0001528901410000063
和保真项εc；；步骤(31)中自适应阈值
Figure BDA0001528901410000064
梯度值
Figure BDA0001528901410000065
和保真项εc的计算方法如下：
Figure BDA0001528901410000066
Figure BDA0001528901410000067
式中，
Figure BDA0001528901410000068
为梯度算子，λ,μ,ρ∈(0,1)，所述ρ由IMF剖面的中心频率决定，中心频率越高，ρ值越大；
(32)根据下式将梯度值
Figure BDA0001528901410000069
与自适应阈值
Figure BDA00015289014100000610
相比较，当
Figure BDA00015289014100000611
时，采取L2-norm滤波；当
Figure BDA00015289014100000612
时，采用TV方法滤波，得到滤波后的IMF剖面
Figure BDA00015289014100000613
Figure BDA00015289014100000614
(33)再以
Figure BDA00015289014100000615
为基础，重复步骤(31)(32)迭代m次，得到滤波后的IMF剖面
Figure BDA00015289014100000616
j为迭代次数；
(34)重复步骤(31)-(33)，得到每个IMF剖面对应的
Figure BDA00015289014100000617
(35)将每个IMF剖面中迭代次数相同的
Figure BDA0001528901410000071
进行叠加重构，得到重构后的最终结果
Figure BDA0001528901410000072
(4)将多次迭代得到的多个最终结果
Figure BDA0001528901410000073
分别计算其信噪比SNR和相似度SSIM；步骤(4)中计算信噪比SNR和相似度SSIM的方法为：
Figure BDA0001528901410000074
Figure BDA0001528901410000075
式中，μi是平均强度，δi表示标准偏差，这里的i∈x,y，δxy表示图像x和y的协方差，Cj＝(KjL)2，其中L是像素值的自适应范围，Kj＜＜1为常数。
(5)判定信噪比SNR和相似度SSIM，选取最优的SNR和SSIM对应的最终结果
Figure BDA0001528901410000077
作为最终结果
Figure BDA0001528901410000078
进行输出。
本发明第一步输入地震数据，第二步利用VMD对二维叠后地震数据进行频率分解，最优分解IMF个数由数据本身确定，一般为4～6个。
第三步采用自适应阈值各向异性滤波算法进行迭代处理，具体来说，就是利用阈值和梯度值对我们的图像进行了划分，针对不同的分类结果选择不同的更适合的滤波方法，这样我们就可以在消除噪声的时候更好的保护纹理的局部特征；
具体方法可通过以下演示来讲解：
为了便于描述，在这里引入一个概念
Figure BDA0001528901410000076
n为分解的IMF剖面个数，m为迭代次数；
首先，输入二维叠后地震数据u并进行VMD分解，假设分解得到4个不同频率范围的IMF剖面uk(k＝1,2,...,n)，为了便于理解，我们标记为：IMF1＝u1、IMF2＝u2、IMF3＝u3、IMF4＝u4；
先按照步骤(31)计算IMF1＝u1的自适应阈值
Figure BDA00015289014100000814
梯度值|▽uk|和保真项εc；
再根据步骤(32)进行滤波，得到
Figure BDA0001528901410000081
此时为IMF1剖面第一次滤波得到的值。
再以
Figure BDA0001528901410000082
为基础，重复步骤(31)(32)迭代。依次迭代19次，得到滤波后的IMF剖面
Figure BDA0001528901410000083
此处对应步骤(33)
然后，根据步骤(34，)得到
IMF2的20个滤波后的IMF剖面
Figure BDA0001528901410000084
IMF3的20个滤波后的IMF剖面
Figure BDA0001528901410000085
IMF4的20个滤波后的IMF剖面
Figure BDA0001528901410000086
再根据步骤(35)进行重构，得到：
Figure BDA0001528901410000087
Figure BDA0001528901410000088
......
Figure BDA0001528901410000089
根据步骤(4)(5)，计算
Figure BDA00015289014100000810
的SNR和SSIM；假设
Figure BDA00015289014100000811
的SNR和SSIM最优，则将
Figure BDA00015289014100000812
作为最终结果
Figure BDA00015289014100000813
进行输出。
实施例2：
如图2所示的人工合成的二维地堑地质模型，其中模型包含250道，每道含有490个采样点，采样频率为1ms。图中设计了两个主要断层，每个断层包含了不同断距的断点。
如图3所示，图3是图2加入信噪比为-1db的随机高斯噪声后的地堑模型剖面。由于噪声的影响，断距断点都变得模糊，难以清楚分辨。
图4为将图3经为本发明方法处理后的剖面图；从图4上部和地层可见，断点恢复清晰，横向边界清楚，弱振幅加强，信噪比也从原来的-1dB提升到13.4575dB，相似度从0.6651提升到0.9775。
实施例3：
如图5、图6所示的是局部原始实际资料剖面及其处理后的剖面。图5可见，实际资料往往含有不同种类的噪声，由于噪声的存在，很多弱振幅被噪声掩盖，小断距消失，地层横向连续性差，从图5与图6对比可见，经过本发明方法处理后，被噪声掩盖的弱振幅得到恢复加强，断距清晰，地层横向连续性加强。
上述实例仅用于说明本发明，其中方法的各实施步骤等同具体实施步骤描述相同/都是可以有所改变的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排出在本发明的保护范围之外。

Claims (2)

1.一种基于能量频带分布的自适应各向异性分频分区滤波方法，其特征在于：包括以下步骤：
(1)输入二维叠后地震数据u；
(2)利用VMD对二维叠后地震数据进行频域分解，得到不同频率范围的IMF剖面uk，k＝1,2,...,n，n为整数且为IMF剖面总数；
VMD对二维叠后地震数据进行频域分解的方法为，
根据下式计算uk，k＝1,2,...,n，
Figure FDA0002404018960000011
其中，δ(t)是狄拉克函数，ωk表示每个频带范围的中心频率；
(3)将分解后的每个IMF剖面uk，k＝1,2,...,n分别经自适应阈值各向异性滤波算法进行多次迭代处理，每次迭代处理后得到一重构的最终结果
Figure FDA0002404018960000012
具体步骤如下：
(31)先计算第一个IMF剖面的自适应阈值
Figure FDA0002404018960000013
梯度值
Figure FDA0002404018960000014
和保真项εc，自适应阈值
Figure FDA0002404018960000015
梯度值
Figure FDA0002404018960000016
和保真项εc的计算方法如下：
Figure FDA0002404018960000017
Figure FDA0002404018960000018
式中，
Figure FDA0002404018960000019
为梯度算子，λ,μ,ρ∈(0,1)，所述ρ由IMF剖面的中心频率决定，中心频率越高，ρ值越大；
(32)根据下式将梯度值
Figure FDA00024040189600000110
与自适应阈值
Figure FDA00024040189600000111
相比较，当
Figure FDA00024040189600000112
时，采取L2-norm滤波；当
Figure FDA0002404018960000021
时，采用TV方法滤波，得到滤波后的IMF剖面
Figure FDA0002404018960000022
Figure FDA0002404018960000023
(33)再以
Figure FDA0002404018960000024
为基础，重复步骤(31)(32)迭代m次，得到滤波后的IMF剖面
Figure FDA0002404018960000025
p为迭代次数；
(34)重复步骤(31)-(33)，得到每个IMF剖面对应的
Figure FDA0002404018960000026
(35)将每个IMF剖面中迭代次数相同的
Figure FDA0002404018960000027
进行叠加重构，得到重构后的最终结果
Figure FDA0002404018960000028
(4)将多次迭代得到的多个最终结果
Figure FDA0002404018960000029
分别计算其信噪比SNR和相似度SSIM；
(5)判定信噪比SNR和相似度SSIM，选取最优的SNR和SSIM对应的最终结果
Figure FDA00024040189600000210
作为最终结果
Figure FDA00024040189600000211
进行输出。
2.根据权利要求1所述的基于能量频带分布的自适应各向异性分频分区滤波方法，其特征在于：步骤(4)中计算信噪比SNR和相似度SSIM的方法为：
Figure FDA00024040189600000212
Figure FDA00024040189600000213
式中，μi是平均强度，δi表示标准偏差，这里的i∈x,y，δxy表示图像x和y的协方差，Cq＝(KqL)2，q∈1,2，其中L是像素值的自适应范围，Kq＜＜1为常数。
CN201711453771.6A 2017-12-28 2017-12-28 基于能量频带分布的自适应各向异性分频分区滤波方法 Active CN108226996B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711453771.6A CN108226996B (zh) 2017-12-28 2017-12-28 基于能量频带分布的自适应各向异性分频分区滤波方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711453771.6A CN108226996B (zh) 2017-12-28 2017-12-28 基于能量频带分布的自适应各向异性分频分区滤波方法

Publications (2)

Publication Number Publication Date
CN108226996A CN108226996A (zh) 2018-06-29
CN108226996B true CN108226996B (zh) 2020-05-22

Family

ID=62648213

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711453771.6A Active CN108226996B (zh) 2017-12-28 2017-12-28 基于能量频带分布的自适应各向异性分频分区滤波方法

Country Status (1)

Country Link
CN (1) CN108226996B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110533602A (zh) * 2019-07-19 2019-12-03 中国石油天然气集团有限公司 基于信噪比场动态约束的潜山内幕成像增强方法和装置
CN113625164A (zh) * 2021-08-02 2021-11-09 南京航空航天大学 航空发电机故障特征提取方法、系统、介质及计算设备

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105549076B (zh) * 2015-12-08 2017-06-09 中国石油天然气股份有限公司 一种基于交替方向法和全变分理论的地震数据处理方法

Also Published As

Publication number Publication date
CN108226996A (zh) 2018-06-29

Similar Documents

Publication Publication Date Title
US8330464B2 (en) Data acquisition method with a three dimensional small bin electromagnetic consecutive array
Górszczyk et al. Application of curvelet denoising to 2D and 3D seismic data—Practical considerations
Wu et al. Noise attenuation for 2-D seismic data by radial-trace time-frequency peak filtering
CN102854533A (zh) 一种基于波场分离原理提高地震资料信噪比的去噪方法
Zhang et al. Signal preserving and seismic random noise attenuation by Hurst exponent based time–frequency peak filtering
CN107179550B (zh) 一种数据驱动的地震信号零相位反褶积方法
CN108226996B (zh) 基于能量频带分布的自适应各向异性分频分区滤波方法
CN104849757A (zh) 消除地震信号中随机噪声系统及方法
CN113238190A (zh) 一种基于emd联合小波阈值的探地雷达回波信号去噪方法
Lv et al. Noise removal for semi-airborne data using wavelet threshold and singular value decomposition
CN114779343A (zh) 一种基于曲波变换-联合双边滤波的地震数据去噪方法
CN112764099B (zh) 一种基于地震几何学信息的地震资料拓频方法
CN112255690B (zh) 基于地震相位分解的自适应围岩强反射分离方法
Zhang et al. Seismic random noise attenuation and signal-preserving by multiple directional time-frequency peak filtering
CN111257931B (zh) 一种去除海洋地震勘探过船干扰噪音的方法
CN109782344B (zh) 沉积层序边界识别方法及装置
CN116068619A (zh) 一种自适应的多阶频散面波压制方法、装置及设备
CN111929726B (zh) 地震相干数据体处理方法及装置
CN111239803B (zh) 基于微分变分模态分解与重构的地震资料处理方法
CN104122583A (zh) 一种拓宽地震数据频谱的方法和装置
CN112764108B (zh) 一种基于改进经验小波变换的新型地震资料噪声压制算法
CN109884705B (zh) 双重约束时频域子波谱提高地震分辨率处理方法
van der Baan et al. Recognition and reconstruction of coherent energy with application to deep seismic reflection data
CN104865603A (zh) 一种针对倾斜层的svd滤波方法及装置
CN112200069A (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