CN110133715B - 一种基于初至时差和波形叠加的微地震震源定位方法 - Google Patents

一种基于初至时差和波形叠加的微地震震源定位方法 Download PDF

Info

Publication number
CN110133715B
CN110133715B CN201910458434.9A CN201910458434A CN110133715B CN 110133715 B CN110133715 B CN 110133715B CN 201910458434 A CN201910458434 A CN 201910458434A CN 110133715 B CN110133715 B CN 110133715B
Authority
CN
China
Prior art keywords
time
arrival time
arrival
function
detector
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
CN201910458434.9A
Other languages
English (en)
Other versions
CN110133715A (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.)
Yangtze University
Original Assignee
Yangtze 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 Yangtze University filed Critical Yangtze University
Priority to CN201910458434.9A priority Critical patent/CN110133715B/zh
Publication of CN110133715A publication Critical patent/CN110133715A/zh
Application granted granted Critical
Publication of CN110133715B publication Critical patent/CN110133715B/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
    • G01V1/288—Event detection in seismic signals, e.g. microseismics

Landscapes

  • Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Environmental & Geological Engineering (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Emergency Management (AREA)
  • Business, Economics & Management (AREA)
  • Acoustics & Sound (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

一种基于初至时差和波形叠加的微地震震源定位方法，步骤如下：S1：输入速度模型；S2：拾取并输入实际初至时间，以及读取地震数据；S3：依据S1中速度模型，计算可行解区域内所有网格点到每一个检波器的理论初至时间表；S4：构造走时残差函数Tr；S5：构造波形叠加函数Ews；S6：输入权重系数β，依据S4中走时残差函数Tr和S5中波形叠加函数Ews，构造改进目标函数；S7：通过网格搜索法寻找改进目标函数的最小值，对应的最优解即震源位置。本发明针对基于走时目标函数的定位方法对初至误差敏感的问题，结合初至时差和波形叠加信息构造改进的目标函数，可增强定位方法的抗噪能力，可提高反演方法的收敛性，从而可改善微地震震源定位的精度。

Description

一种基于初至时差和波形叠加的微地震震源定位方法
技术领域
本发明属于石油地球物理勘探开发技术领域，具体涉及一种基于初至时差和波形叠加的微地震震源定位方法。
背景技术
水力压裂微地震监测技术是监测致密含油气和页岩气等非常规储层改造及效果的油藏地球物理技术，其应用结果可评价压裂效果、调整压裂设计和井网布置，并对下一步开发提供有效指导，从而提高非常规油气藏的产能。
在水力压裂监测资料处理中，其核心部分是微地震震源定位。前人基本上是利用单一的走时信息或者波形信息来建立目标函数，可分为基于走时的射线追踪定位和基于波形的偏移定位方法。而在实际水力压裂实时监测作业中，采集的微地震资料可能存在信噪比不高、能量较弱和背景噪音较大等特点，极大地影响了初至拾取的精度。基于走时信息目标函数的定位方法对初至误差更为敏感，而基于波形信息目标函数的定位方法虽然不需要精确的初至时间，但是其计算量巨大，而导致定位效率不高。为此，亟需结合走时信息和波形信息建立目标函数，研发一种基于初至时差和波形叠加的微地震震源定位方法。
发明内容
本发明的目的在于克服上述背景技术的不足，而提供一种基于初至时差和波形叠加的微地震震源定位方法，该方法利用目标函数中微地震资料的走时和波形信息同时约束定位，并且具有一定抗噪能力，可提高反演方法的收敛性，以改善微地震震源定位的精度。
为实现上述目的，本发明提供的一种基于初至时差和波形叠加的微地震震源定位方法，针对任意一个微地震震源(事件)定位而言，包括以下步骤：
S1：输入速度模型；
S2：拾取并输入实际初至时间，以及读取地震数据；
S3：依据步骤S1中速度模型，计算可行解区域内所有网格点到每一个检波器的理论初至时间表；
S4：基于步骤S2中输入的实际初至时间和步骤S3中计算的理论初至时间表，依据公式(1)和(2)构造走时残差函数Tr；
Figure BDA0002077311450000021
Figure BDA0002077311450000022
上式中：M、N分别是观测到的横波和纵波初至的个数；
Figure BDA0002077311450000023
和
Figure BDA0002077311450000024
分别是第j个检波器记录的横波和纵波初至时间，
Figure BDA0002077311450000025
和
Figure BDA0002077311450000026
分别是与之对应的第j个检波器的横波及纵波理论计算时间；
Figure BDA0002077311450000027
和
Figure BDA0002077311450000028
分别是第i个检波器记录的横波和纵波初至时间，
Figure BDA0002077311450000031
和
Figure BDA0002077311450000032
分别是与之对应的第i个检波器的横波及纵波理论计算时间；Tshift是观测及理论初至之间的一个恒定的漂移量，以解决实际资料监测中微地震震源发震时间未知的难题；
S5：将步骤S2中读入的地震数据、走时目标函数对应微地震事件中第m个检波器记录的可识别纵波初至时间
Figure BDA0002077311450000033
和步骤S3中理论初至时间表结合，通过极性修正，然后依据公式(3)、(4)和(5)构造波形叠加函数Ews；
Figure BDA0002077311450000034
Figure BDA0002077311450000035
Figure BDA0002077311450000036
式中ΔP、ΔS分别是纵横波滑动时窗样点数；dt是采样间隔；n是检波器个数；
Figure BDA0002077311450000037
分别是假设震源点到第i个检波器的纵横波理论时间
Figure BDA0002077311450000038
与第m个检波器的纵波理论时间
Figure BDA0002077311450000039
之差，
Figure BDA00020773114500000310
分别是第i个检波器纵横波振幅的符号系数，ui是地震数据第i道振幅值；
S6：输入权重系数β，依据步骤S4中走时残差函数Tr和步骤S5中波形叠加函数Ews，按照公式(6)构造改进目标函数；
Q＝Tr-βEws (6)
S7：通过网格搜索法寻找改进目标函数的最小值，输出此时对应的最优解，即震源位置。
上述技术方案中，在所述步骤S1中，速度模型可通过声波测井资料和地质分层资料建立初始速度模型，然后利用射孔资料进行速度模型校正。
上述技术方案中，在所述步骤S3中，可行解区域内所有网格点是指以射孔点为中心，给定一个求解区域，然后对该区域进行网格化，每个网格点均为可能的震源点。
上述技术方案中，在所述步骤S5中，极性修正可通过符号系数
Figure BDA0002077311450000041
实现，即判断检波器记录中时窗内振幅符号，需要正负号保持一致，此处保证波形的同向叠加。
上述技术方案中，在所述步骤S6中，权重系数β>0，大小与微地震资料的信噪比情况成反比。
上述技术方案中，在所述步骤S7中，在网格搜索定位中，假设网格点是震源位置时，Tr为最小，Ews为最大(即-βEws为最小)，故改进目标函数Q为最小值。
与现有技术相比，本发明具有如下优点：
本发明利用改进目标函数中微地震资料的走时和波形信息同时约束定位，增强了定位方法的抗噪能力，提高了反演方法的收敛性，可改善微地震震源定位的精度。
附图说明
图1是本发明的实施流程示意图；
图2是基于均匀介质模型的二维微地震井中观测***，横、纵坐标分别表示水平距离x和深度z；
图3a是理论模型的z分量地震记录；
图3b是加入随机噪音后理论模型的z分量地震记录；
图4a是基于走时目标函数方法的定位结果图，横、纵坐标分别表示水平距离x和深度z；
图4b是本发明方法的定位结果图，横、纵坐标分别表示水平距离x和深度z。
具体实施方式
下面结合附图和实施例对本发明的实施方式做进一步的详细说明，但它们并不构成对本发明的限定，仅作举例而已。同时通过说明本发明的优点将变得更加清楚和容易理解。
本实施例的一种基于初至时差和波形叠加的微地震震源定位方法，如图1所示，包括以下步骤：
S1：输入速度模型；
S2：拾取并输入实际初至时间，以及读取地震数据；
S3：依据步骤S1中速度模型，计算可行解区域内所有网格点到每一个检波器的理论初至时间表；
S4：基于步骤S2中输入的实际初至时间和步骤S3中计算的理论初至时间表，依据公式(1)和(2)构造走时残差函数Tr；
Figure BDA0002077311450000051
Figure BDA0002077311450000061
上式中：M、N分别是观测到的横波和纵波初至的个数；
Figure BDA0002077311450000062
和
Figure BDA0002077311450000063
分别是第j个检波器记录的横波和纵波初至时间，
Figure BDA0002077311450000064
和
Figure BDA0002077311450000065
分别是与之对应的第j个检波器的横波及纵波理论计算时间；
Figure BDA0002077311450000066
和
Figure BDA0002077311450000067
分别是第i个检波器记录的横波和纵波初至时间，
Figure BDA0002077311450000068
和
Figure BDA0002077311450000069
分别是与之对应的第i个检波器的横波及纵波理论计算时间；Tshift是观测及理论初至之间的一个恒定的漂移量，以解决实际资料监测中微地震震源发震时间未知的难题；
S5：将步骤S2中读入的地震数据、走时目标函数对应微地震事件中第m个检波器记录的可识别纵波初至时间
Figure BDA00020773114500000610
和步骤S3中理论初至时间表结合，通过极性修正，然后依据公式(3)、(4)和(5)构造波形叠加函数Ews；
Figure BDA00020773114500000611
Figure BDA00020773114500000612
Figure BDA00020773114500000613
式中ΔP、ΔS分别是纵横波滑动时窗样点数；dt是采样间隔；n是检波器个数；
Figure BDA00020773114500000614
分别是假设震源点到第i个检波器的纵横波理论时间
Figure BDA00020773114500000615
与第m个检波器的纵波理论时间
Figure BDA00020773114500000616
之差，
Figure BDA00020773114500000617
分别是第i个检波器纵横波振幅的符号系数，ui是地震数据第i道振幅值；
S6：输入权重系数β，依据步骤S4中走时残差函数Tr和步骤S5中波形叠加函数Ews，按照公式(6)构造改进目标函数；
Q＝Tr-βEws (6)
S7：通过网格搜索法寻找改进目标函数的最小值，输出此时对应的最优解，即震源位置。
上述技术方案中，在所述步骤S1中，速度模型可通过声波测井资料和地质分层资料建立初始速度模型，然后利用射孔资料进行速度模型校正。
上述技术方案中，在所述步骤S3中，可行解区域内所有网格点是指以射孔点为中心，给定一个求解区域，然后对该区域进行网格化，每个网格点均为可能的震源点。
上述技术方案中，在所述步骤S5中，用大致准确的纵波初至时间
Figure BDA0002077311450000071
代替了常规波形叠加方法中未知的震源发震时间。
上述技术方案中，在所述步骤S5中，极性修正可通过符号系数
Figure BDA0002077311450000072
实现，即判断检波器记录中时窗内振幅符号，需要正负号保持一致，此处保证波形的同向叠加。
上述技术方案中，在所述步骤S6中，权重系数β大小(β>0)取决于微地震资料的信噪比情况，如果信噪比较低，则波形叠加部分权重要大些，故设置β大些，反之设置β小些。
上述技术方案中，在所述步骤S7中，在网格搜索定位中，假设网格点是震源位置时，Tr为最小，Ews为最大(即-βEws为最小)，故改进目标函数Q为最小值。
实施例：首先建立基于均匀介质模型的二维微地震井中观测***，纵波速度Vp＝4500m/s，横波速度Vs＝2500m/s，密度为2.425g/cm3。如图2所示，水平距离x＝2000m处监测井(直井)中10级检波器深度z分别位于1500-1590m，其道间距为10m。设置微地震震源发生的位置为(x，z)＝(2200，1570)m。然后采用100Hz的雷克子波进行二维弹性波动方程正演模拟，采样间隔为0.5ms，得到该理论模型的z分量地震记录(图3a)，图3b是对图3a加入随机噪音后理论模型的z分量地震记录。
为了测试本发明方法对初至误差的稳定性，对10个检波器微地震事件初至分别加入200次[-2，2]ms随机误差。然后依据公式(1)和(2)建立走时残差目标函数后，采用步骤7中网格搜索法对微地震震源进行定位，得到基于走时目标函数方法的定位结果(图4a)。依据公式(1)至(6)建立本发明的目标函数，采用步骤7中网格搜索法对微地震震源进行定位，得到本发明方法的定位结果(图4b)。
从图4a与图4b对比可以看出，在一定初至误差下，图4b的定位结果比图4a更加收敛，即本发明定位误差要比基于走时目标函数的定位方法误差小的多。
综上所述，本发明基于初至时差和波形叠加相结合的思路，对微地震震源进行定位，可以提高定位方法的抗噪能力及收敛性，可改善微地震震源定位的精度。
本说明书未作详细描述的内容属于本领域专业技术人员公知的现有技术。

Claims (5)

1.一种基于初至时差和波形叠加的微地震震源定位方法，其特征在于，包括以下步骤：
S1：输入速度模型；
S2：拾取并输入实际初至时间，以及读取地震数据；
S3：依据步骤S1中速度模型，计算可行解区域内所有网格点到每一个检波器的理论初至时间表；
所述可行解区域内所有网格点是指以射孔点为中心，给定一个求解区域，然后对该区域进行网格化，每个网格点均为可能的震源点；
S4：基于步骤S2中输入的实际初至时间和步骤S3中计算的理论初至时间表，依据公式(1)和(2)构造走时残差函数Tr；
Figure FDA0002728497220000011
Figure FDA0002728497220000012
上式中：M、N分别是观测到的横波和纵波初至的个数；
Figure FDA0002728497220000013
和
Figure FDA0002728497220000014
分别是第j个检波器记录的横波和纵波初至时间，
Figure FDA0002728497220000015
和
Figure FDA0002728497220000016
分别是与之对应的第j个检波器的横波及纵波理论计算时间；
Figure FDA0002728497220000017
和
Figure FDA0002728497220000018
分别是第i个检波器记录的横波和纵波初至时间，
Figure FDA0002728497220000019
和
Figure FDA00027284972200000110
分别是与之对应的第i个检波器的横波及纵波理论计算时间；Tshift是观测及理论初至之间的一个恒定的漂移量，以解决实际资料监测中微地震震源发震时间未知的难题；
S5：将步骤S2中读入的地震数据、步骤S4中的走时残差函数对应微地震事件中第m个检波器记录的可识别纵波初至时间
Figure FDA0002728497220000021
和步骤S3中理论初至时间表结合，通过极性修正，然后依据公式(3)、(4)和(5)构造波形叠加函数Ews；
Figure FDA0002728497220000022
Figure FDA0002728497220000023
Figure FDA0002728497220000024
式中ΔP、ΔS分别是纵、横波滑动时窗样点数；dt是采样间隔；n是检波器个数；
Figure FDA0002728497220000025
分别是假设震源点到第i个检波器的纵、横波理论时间
Figure FDA0002728497220000026
与第m个检波器的纵波理论时间
Figure FDA0002728497220000027
之差，
Figure FDA0002728497220000028
分别是第i个检波器纵、横波振幅的符号系数，ui是地震数据第i道振幅值；
S6：输入权重系数β，依据步骤S4中走时残差函数Tr和步骤S5中波形叠加函数Ews，按照公式(6)构造改进目标函数；
Q＝Tr-βEws (6)
S7：通过网格搜索法寻找改进目标函数的最小值，输出此时对应的最优解，即震源位置。
2.根据权利要求1所述的一种基于初至时差和波形叠加的微地震震源定位方法，其特征在于：在所述步骤S1中，速度模型可通过声波测井资料和地质分层资料建立初始速度模型，然后利用射孔资料进行速度模型校正。
3.根据权利要求1所述的一种基于初至时差和波形叠加的微地震震源定位方法，其特征在于：在所述步骤S5中，极性修正可通过符号系数
Figure FDA0002728497220000031
实现，即判断检波器记录中时窗内振幅符号，需要正负号保持一致，此处保证波形的同向叠加。
4.根据权利要求1所述的一种基于初至时差和波形叠加的微地震震源定位方法，其特征在于：在所述步骤S6中，权重系数β>0，大小与微地震资料的信噪比情况成反比。
5.根据权利要求1所述的一种基于初至时差和波形叠加的微地震震源定位方法，其特征在于：在所述步骤S7中，在网格搜索定位中，假设网格点是震源位置时，Tr为最小，Ews为最大，即-βEws为最小，故改进目标函数Q为最小值。
CN201910458434.9A 2019-05-29 2019-05-29 一种基于初至时差和波形叠加的微地震震源定位方法 Active CN110133715B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910458434.9A CN110133715B (zh) 2019-05-29 2019-05-29 一种基于初至时差和波形叠加的微地震震源定位方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910458434.9A CN110133715B (zh) 2019-05-29 2019-05-29 一种基于初至时差和波形叠加的微地震震源定位方法

Publications (2)

Publication Number Publication Date
CN110133715A CN110133715A (zh) 2019-08-16
CN110133715B true CN110133715B (zh) 2021-01-08

Family

ID=67582702

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910458434.9A Active CN110133715B (zh) 2019-05-29 2019-05-29 一种基于初至时差和波形叠加的微地震震源定位方法

Country Status (1)

Country Link
CN (1) CN110133715B (zh)

Families Citing this family (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110568499B (zh) * 2019-08-20 2021-06-04 长江大学 一种vsp地震资料的初至波时差校正方法及装置
CN111077569B (zh) * 2019-12-23 2022-05-06 中国石油天然气股份有限公司 全波形反演中分时窗提取数据的方法及装置
CN111221034B (zh) * 2020-01-20 2022-02-25 山东黄金矿业股份有限公司新城金矿 矿山微地震源定位方法及模拟检验***
CN111352160B (zh) * 2020-03-19 2020-11-10 中国科学院地质与地球物理研究所 一种海底地震仪自动重定位装置及方法
CN111580165A (zh) * 2020-05-27 2020-08-25 中国科学院地质与地球物理研究所 一种海底地震仪到时差定位装置及方法
CN111856581B (zh) * 2020-07-27 2022-02-22 广州海洋地质调查局 一种obs时钟漂移校正方法及处理终端
CN115327620B (zh) * 2021-05-11 2023-07-28 中国石油化工股份有限公司 微地震联合时差叠加定位方法
CN113608257A (zh) * 2021-07-07 2021-11-05 长江大学 一种基于改进波形叠加函数的微地震事件偏移定位方法
CN113960532B (zh) * 2021-10-20 2024-05-24 西北大学 一种基于假想源的二次定位计算的微地震定位方法
CN114167495B (zh) * 2021-11-30 2023-08-11 长江大学 一种用于减少纵波压制的叠加自相关滤波方法及装置
CN115373029B (zh) * 2022-10-25 2023-01-31 中国科学院地质与地球物理研究所 基于深度学习的实时微地震震源机制计算方法及***

Family Cites Families (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7978563B2 (en) * 2009-08-18 2011-07-12 Microseismic, Inc. Method for passive seismic emission tomography including polarization correction for source mechanism
CN102841373B (zh) * 2012-08-23 2015-02-04 中国石油集团川庆钻探工程有限公司地球物理勘探公司 基于方位角约束的微地震事件定位方法
CN103399300B (zh) * 2013-07-31 2015-07-08 中国石油集团川庆钻探工程有限公司地球物理勘探公司 波包叠加微地震地面定位方法
CN105510880A (zh) * 2014-09-23 2016-04-20 中国石油化工股份有限公司 一种基于双差法的微地震震源定位方法
CN106353821B (zh) * 2015-07-17 2020-06-30 中国石油化工股份有限公司 一种微地震事件定位方法
CN106353792B (zh) * 2015-07-17 2021-03-23 中国石油化工股份有限公司 一种适用于水力压裂微震震源定位的方法
CN105954795A (zh) * 2016-04-25 2016-09-21 吉林大学 一种用于微地震定位的网格逐次剖分方法
CN109085642B (zh) * 2017-06-14 2020-05-15 中国石油化工股份有限公司 一种各向异性介质微地震事件定位方法
CN109212594B (zh) * 2017-07-01 2020-04-07 中国石油化工股份有限公司 一种各向异性介质纵横波联合定位方法
CN109031419A (zh) * 2018-07-27 2018-12-18 长江大学 一种拾取微地震初至的方法及***

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
"Automated Microseismic Event Location Using Finite Difference Travel time Calculation and Enhanced Waveform Stacking";J.W. Huang等;《Conference Proceedings, 75th EAGE Conference & Exhibition incorporating SPE EUROPEC 2013》;20130731;第348-00943页 *
"Comparison of migration-based location and detection methods for microseismic events";Trojanowski J.等;《Geophysical Prospecting》;20171231;第65卷(第1期);第47-63页 *
"Improved relative locations of clustered earthquakes using con-strained multiple event location";Fehler M.等;《Bulletin of the Seismological Society of America》;20001231;第90卷(第3期);第775-780页 *

Also Published As

Publication number Publication date
CN110133715A (zh) 2019-08-16

Similar Documents

Publication Publication Date Title
CN110133715B (zh) 一种基于初至时差和波形叠加的微地震震源定位方法
CN111239802B (zh) 基于地震反射波形和速度谱的深度学习速度建模方法
CN109425896B (zh) 白云岩油气储层分布预测方法及装置
CN105510880A (zh) 一种基于双差法的微地震震源定位方法
CN109669212B (zh) 地震数据处理方法、地层品质因子估算方法与装置
CN102053263B (zh) 调查表层结构的方法
WO2012139082A1 (en) Event selection in the image domain
Wuestefeld et al. Benchmarking earthquake location algorithms: A synthetic comparison
CN111722284B (zh) 一种基于道集数据建立速度深度模型的方法
CN110687602A (zh) 浅层地震多波联合勘探方法
CN109991658B (zh) 一种基于“震源-台站”速度模型的微地震事件定位方法
CN112305591B (zh) 隧道超前地质预报方法、计算机可读存储介质
CN104570116A (zh) 基于地质标志层的时差分析校正方法
CN105607119B (zh) 近地表模型构建方法与静校正量求取方法
EP2321671A2 (en) Processing seismic data in common group-center gathers
CN110780341B (zh) 一种各向异性地震成像方法
CN103645506B (zh) 一种检测地层中裂缝发育程度的方法
RU2490677C2 (ru) Способ комплексной обработки геофизических данных и технологическая система "литоскан" для его осуществления
CN108375789B (zh) 联合采集地震数据的同步匹配方法
CN108693560B (zh) 一种基于互相关道的散射波成像方法及***
Strahser et al. Polarisation and slowness of seismoelectric signals: a case study
CN111352153A (zh) 一种基于瞬时相位互相关加权的微地震干涉定位方法
CN113534236B (zh) 一种基于检波器间距约束的微地震初至拾取方法
CN110764148B (zh) 一种各向异性矢量波场井地联合定位方法
CN110764136B (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