CN112083487A - 宽频带频散曲线提取方法、装置 - Google Patents
宽频带频散曲线提取方法、装置 Download PDFInfo
- Publication number
- CN112083487A CN112083487A CN202010978546.XA CN202010978546A CN112083487A CN 112083487 A CN112083487 A CN 112083487A CN 202010978546 A CN202010978546 A CN 202010978546A CN 112083487 A CN112083487 A CN 112083487A
- Authority
- CN
- China
- Prior art keywords
- array
- frequency
- dispersion
- stations
- energy
- 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
Links
- 239000006185 dispersion Substances 0.000 title claims abstract description 400
- 238000000034 method Methods 0.000 title claims abstract description 59
- 230000010363 phase shift Effects 0.000 claims abstract description 137
- 238000004364 calculation method Methods 0.000 claims abstract description 27
- 238000005314 correlation function Methods 0.000 claims description 31
- 238000000605 extraction Methods 0.000 claims description 24
- 239000011159 matrix material Substances 0.000 claims description 13
- 238000001514 detection method Methods 0.000 claims description 11
- 230000004927 fusion Effects 0.000 claims description 10
- 230000009466 transformation Effects 0.000 claims description 5
- 238000007781 pre-processing Methods 0.000 claims description 3
- 230000002349 favourable effect Effects 0.000 abstract 1
- 238000010586 diagram Methods 0.000 description 15
- 230000000694 effects Effects 0.000 description 8
- 230000000737 periodic effect Effects 0.000 description 7
- 238000003491 array Methods 0.000 description 6
- 238000012545 processing Methods 0.000 description 5
- 238000004590 computer program Methods 0.000 description 4
- 239000000284 extract Substances 0.000 description 4
- 238000000638 solvent extraction Methods 0.000 description 4
- 238000004458 analytical method Methods 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 3
- 238000003384 imaging method Methods 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 230000006870 function Effects 0.000 description 2
- 238000005065 mining Methods 0.000 description 2
- 230000003287 optical effect Effects 0.000 description 2
- 238000012952 Resampling Methods 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- PCHJSUWPFVWCPO-UHFFFAOYSA-N gold Chemical compound [Au] PCHJSUWPFVWCPO-UHFFFAOYSA-N 0.000 description 1
- 239000010931 gold Substances 0.000 description 1
- 229910052737 gold Inorganic materials 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 238000004806 packaging method and process Methods 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 239000004065 semiconductor Substances 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 239000000758 substrate Substances 0.000 description 1
- 239000013589 supplement Substances 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 230000002087 whitening effect Effects 0.000 description 1
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/30—Analysis
- G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
-
- 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
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)
- Geophysics And Detection Of Objects (AREA)
Abstract
一种宽频带频散曲线提取方法,包括:利用密集线性台阵内的各台站采集地震面波数据,分别以密集台阵内的各个台站为中心台站,将密集台阵按一定孔径划分为多个阵内台站和多个阵外台站,基于与各中心台站对应的多个阵内台站采集的地震面波数据,计算阵内相移频散能量,以及,基于全部台站采集的地震面波数据,计算阵外相移频散能量,提取出对应的阵内频散曲线和阵外频散曲线,给各中心台站对应的阵内频散曲线和阵外频散曲线在各频率点分别分配权重,并进行加权叠加计算,使阵内频散曲线和阵外频散曲线融合为一条宽频带频散曲线。本公开提供的方法有效改善了地震背景噪声中低频面波的频散能量的质量,为更容易的提取出宽频带频散曲线创造了有利条件。
Description
技术领域
本公开涉及地震背景噪声成像技术领域,尤其涉及一种宽频带频散曲线提取方法、装置。
背景技术
地震背景噪声可用来获取地下结构的速度信息,其中面波能量占据地表地震记录的70%以上,通过提取台阵下方的面波频散曲线可以反演得到横波速度结构。频散曲线提取是背景噪声面波勘探中最基础、最重要的环节,频散曲线中的信息量直接影响勘探的最终结果。其中,频散曲线中的高频频散信息能够用于反演更精细的浅层地下结构,而频散曲线中的低频频散信息能够用于反演更深部的地下结构。然而,高频和低频频散曲线的获取之间存在着矛盾,同时提取低频频散曲线和高频频散曲线提取较为困难。
Park等人(Imaging dispersion curves of surface waves on multi-channelrecord,SEG Technical Program Expanded Abstracts,1998)提出的主动源相移法被广泛应用于提取十几Hz以上的高频频散曲线,通常只能约束到距离地表30m以浅的范围内;Park等人(Multichannel Analysis of Passive Surface Waves-Modeling and ProcessingSchemes,Site Characterization and Modeling,2005)后来对主动源相移法进行了拓展,通过布设“十”字正交型台阵使其具备处理被动源地震记录的能力,可提取出更低的频率,约数Hz;Cheng等人(Multichannel analysis of passive surface waves based oncross-correlations,Geophysics,2016)在该方法中引入互相关函数,使表达式和计算更为简练,并且不再依赖于“十”字台阵布设,提高了相移法的实际应用范围。但目前被动源相移法想要突破到1Hz以下依然是非常困难的。然而浅地表勘探对有效探测深度的要求越来越高,甚至达到了公里级,传统的浅地表面波勘探技术不再能满足新的发展需要。因此,需要一种新的频散提取方法获取更宽频带的频散曲线,以实现对地下结构由浅至深的精细刻画。
发明内容
鉴于上述问题,本发明提供了一种宽频带频散曲线提取方法,应用于密集线性台阵,包括:利用所述密集线性台阵内的各台站采集地震面波数据;分别以所述密集台阵内的各个台站为中心台站,将所述密集台阵按一定半径划分为多个阵内台站和多个阵外台站;基于各所述中心台站对应的多个阵内台站采集的地震面波数据,计算各所述中心台站对应的阵内相移频散能量,以及,基于全部台站采集的地震面波数据,计算各所述中心台站对应的阵外相移频散能量;根据各所述中心台站对应的所述阵内相移频散能量和所述阵外相移频散能量,提取出对应的阵内频散曲线和阵外频散曲线;给各所述中心台站对应的所述阵内频散曲线和所述阵外频散曲线在各频率点分别分配权重,并进行加权叠加计算,使所述阵内频散曲线和所述阵外频散曲线融合为一条宽频带频散曲线。
可选的,所述利用所述密集线性台阵内的各台站采集地震面波数据之后,方法还包括:对所述地震面波数据分别进行单道预处理。
可选的,所述分别以所述密集台阵内的各个台站为中心台站,将所述密集台阵按一定半径划分为多个阵内台站和多个阵外台站包括:若所述密集线性台阵为台间距为d的N个台站组成,当以第n个台站为中心台站时,以r=md为半径作圆,所述圆内部的多个台站为所述多个阵内台站,所述圆外部的多个台站成为所述多个阵外台站,n=1,…,N。
可选的,所述基于各所述中心台站对应的多个阵内台站采集的地震面波数据,计算各所述中心台站对应的阵内相移频散能量,以及,基于全部台站采集的地震面波数据,计算各所述中心台站对应的阵外相移频散能量,包括:基于各所述中心台站对应的所述多个阵内台站采集的所述地震面波数据,利用互相关函数,计算与各所述中心台站对应的第一正支频散能量、第一负支频散能量和第一对称叠加频散能量;从所述第一正支频散能量、第一负支频散能量和第一对称叠加频散能量中选出一种频散能量,作为所述中心台站的阵内相移频散能量;基于全部台站采集的所述地震面波数据,利用互相关函数,计算与各所述中心台站对应的第二正支频散能量、第二负支频散能量和第二对称叠加频散能量;从所述第二正支频散能量、第二负支频散能量和第二对称叠加频散能量中选出一种频散能量,作为所述中心台站的阵外相移频散能量。
可选的,所述基于各所述中心台站对应的所述多个阵内台站采集的所述地震面波数据,利用互相关函数,计算与各所述中心台站对应的第一正支频散能量、第一负支频散能量和第一对称叠加频散能量包括:令表示所述第一正支频散能量,表示所述第一负支频散能量,表示所述第一对称叠加频散能量,f表示地震面波的频率,CT表示所述地震面波的估计相速度,xkl表示第k个阵内台站和第l个阵内台站之间的距离,分别为所述k,l两台站记录的地震面波数据的互相关函数正支和负支的傅里叶变换,所述多个阵内台站为所述密集线性台阵内的第n-m至n+m个台站,则:
可选的,所述基于各所述中心台站对应的所述多个阵外台站采集的所述地震面波数据,利用互相关函数,计算与各所述中心台站对应的第二正支频散能量、第二负支频散能量和第二对称叠加频散能量包括:令表示所述第二正支频散能量,表示所述第二负支频散能量,表示所述第二对称叠加频散能量,f表示所述地震面波的频率,CT表示所述地震面波的估计相速度,xkl表示第k个所述阵外台站和第l个阵内台站之间的距离,分别为k,l两台站记录的地震面波数据的互相关函数正支和负支的傅里叶变换,所述密集线性台阵包括N个台站,所述多个阵内台站为所述密集线性台阵内的第n-m至n+m个台站,所述多个阵外台站为所述密集线性台阵内的除所述多个阵内台站以外的台站,则:
可选的,所述根据各所述中心台站对应的所述阵内相移频散能量和所述阵外相移频散能量,提取出对应的阵内频散曲线和阵外频散曲线包括:针对各所述中心台站,获得其采集到的不同频率的所述地震面波对应的阵内相移频散能量的最大值,以估计各频率的所述地震面波对应的第一真实相速度;基于所述第一真实相速度和对应的所述地震面波的频率,得到各所述中心台站对应的所述阵内频散曲线;针对各所述中心台站,获得不同频率的所述地震面波对应的阵外相移频散能量的最大值,以估计各频率的所述地震面波对应的第二真实相速度;基于所述第二真实相速度和对应的所述地震面波的频率,得到各所述中心台站对应的所述阵外频散曲线。
可选的,所述给各所述中心台站对应的所述阵内频散曲线和所述阵外频散曲线在各频率点分别分配权重,并进行加权叠加计算,使所述阵内频散曲线和所述阵外频散曲线融合为一条宽频带频散曲线包括:比较所述中心台站对应的阵内频散曲线和阵外频散曲线,选取阵内频散曲线和阵外频散曲线最接近的频段;将所述频段内的所述阵内相移频散曲线和所述阵外相移频散曲线进行加权叠加,融合为一条宽频带频散曲线。
可选的,所述将所述频段内的所述阵内相移频散曲线和所述阵外相移频散曲线进行加权叠加,融合为一条宽频带频散曲线包括:令C(f)、Cint(f)、Cext(f)分别表示在所述宽频带频散曲线、阵内相移频散曲线、阵外相移频散曲线中频率为f的所述地震面波对应的相速度,ωint(f)表示所述阵内相移频散曲线中频率点f对应的相速度的权重,ωext(f)表示所述阵外相移频散曲线中频率点f对应的相速度权重,则:
C(f)=ωint(f)·Cint(f)+ωext(f)·Cext(f)。
本公开还提供了一种宽频带频散曲线提取装置,包括:探测模块,用于利用所述密集线性台阵内的各台站采集地震面波数据;台阵划分模块,用于分别以所述密集台阵内的各个台站为中心台站,将所述密集台阵按一定半径划分为多个阵内台站和多个阵外台站;频散能量计算模块,用于基于各所述中心台站对应的多个阵内台站采集的地震面波数据,计算各所述中心台站对应的阵内相移频散能量,以及,基于全部台站采集的地震面波数据,计算各所述中心台站对应的阵外相移频散能量;频散曲线提取模块,用于根据各所述中心台站对应的所述阵内相移频散能量和所述阵外相移频散能量,提取出对应的阵内频散曲线和阵外频散曲线;曲线融合模块,用于给各所述中心台站对应的所述阵内频散曲线和所述阵外频散曲线在各频率点分别分配权重,并进行加权叠加计算,使所述阵内频散曲线和所述阵外频散曲线融合为一条宽频带频散曲线。
在本公开实施例采用的上述至少一个技术方案能够达到以下有益效果:
本公开提供的一种宽频带频散曲线提取方法,应用于密集线性台阵,将密集线性台阵分为分为阵内台阵和阵外台阵,分别计算阵内相移频散能量和阵外相移频散能量。阵内频散能量用来提取中高频频散曲线,面波频散高频信号的保留保证了浅层结构的分辨率;阵外频散能量用来提取中低频频散曲线,改善了面波频散的低频信号质量,极大的增加了有效探测深度。通过该方法提取的频散曲线,低频段由数赫兹降到了约0.55赫兹,信噪比更高的可达到0.5Hz,拓宽了提取频散曲线的频带。该方法在保证对浅地表结构具有较高分辨率同时增大有效探测深度,从而有助于对地下结构进行从深到浅的整体研究。
附图说明
为了更完整地理解本公开及其优势,现在将参考结合附图的以下描述,其中:
图1示意性示出了本公开实施例提供的一种宽频带频散曲线提取方法的示意图;
图2示意性示出了本公开实施例提供的一种划分阵内台站和阵外台站的示意图;
图3示意性示出了本公开实施例提供的一种阵内频散曲线和阵外频散曲线融合的示意图;
图4示意性示出了本公开实施例提供的一种面波的互相关信号的示意图;
图5示意性示出了本公开实施例提供的一种台站分布图;
图6示意性示出了本公开实施例提供的一条SmartSolo测线的互相关结果示意图;
图7示意性示出了本公开实施例提供的一种经典相移法和阵外相移频散能量对比图;
图8示意性示出了本公开实施例提供的一种SmartSolo测线频散曲线融合示意图;
图9示意性示出了本公开实施例提供的对SmartSolo测线经典相移法提取频散曲线结果示意图;
图10示意性示出了本公开实施例提供的一种SmartSolo和Zland 2条测线利用拓距相移法提取的全部频散的示意图。
具体实施方式
以下,将参照附图来描述本公开的实施例。但是应该理解,这些描述只是示例性的,而并非要限制本公开的范围。在下面的详细描述中,为便于解释,阐述了许多具体的细节以提供对本公开实施例的全面理解。然而,明显地,一个或多个实施例在没有这些具体细节的情况下也可以被实施。此外,在以下说明中,省略了对公知结构和技术的描述,以避免不必要地混淆本公开的概念。
在此使用的术语仅仅是为了描述具体实施例,而并非意在限制本公开。在此使用的术语“包括”、“包含”等表明了所述特征、步骤、操作和/或部件的存在,但是并不排除存在或添加一个或多个其他特征、步骤、操作或部件。
在此使用的所有术语(包括技术和科学术语)具有本领域技术人员通常所理解的含义,除非另外定义。应注意,这里使用的术语应解释为具有与本说明书的上下文相一致的含义,而不应以理想化或过于刻板的方式来解释。
附图中示出了一些方框图和/或流程图。应理解,方框图和/或流程图中的一些方框或其组合可以由计算机程序指令来实现。这些计算机程序指令可以提供给通用计算机、专用计算机或其他可编程数据处理装置的处理器,从而这些指令在由该处理器执行时可以创建用于实现这些方框图和/或流程图中所说明的功能/操作的装置。
因此,本公开的技术可以硬件和/或软件(包括固件、微代码等)的形式来实现。另外,本公开的技术可以采取存储有指令的计算机可读介质上的计算机程序产品的形式,该计算机程序产品可供指令执行系统使用或者结合指令执行系统使用。在本公开的上下文中,计算机可读介质可以是能够包含、存储、传送、传播或传输指令的任意介质。例如,计算机可读介质可以包括但不限于电、磁、光、电磁、红外或半导体系统、装置、器件或传播介质。计算机可读介质的具体示例包括:磁存储装置,如磁带或硬盘(HDD);光存储装置,如光盘(CD-ROM);存储器,如随机存取存储器(RAM)或闪存;和/或有线/无线通信链路。
经典相移法是用于提取面波频散曲线的有效方法,该方法对面波中高频部分的频散曲线提取效果较好,但对面波中的低频部分的频散曲线的提取效果较差。
在本公开实施例中,为保证对地下浅部结构的探测具有较高分辨率同时增大有效探测深度,本公开提供了一种宽频带频散曲线提取方法,包括:利用经典相移法计算阵内频散能量提取中高频频散曲线;利用阵外相移法计算阵外频散能量提取中低频频散曲线。本公开提供的该宽频带频散曲线提取方法在不损失高频频散信息的同时,有效改善了低频频散曲线的信号质量。
参阅图1,本公开提供了一种宽频带频散曲线提取方法,应用于密集线性台阵,包括步骤S110~S150。
S110,利用密集线性台阵内的各台站采集地震面波数据。
在本公开实施例中,采集到地震面波数据之后,对地震面波数据分别进行单道预处理,以保证采集到的数据可靠性,包括数据切割、去除仪器响应、重采样、去均值、去线性、带通滤波、时域归一化、谱白化。
S120,分别以密集台阵内的各个台站为中心台站,将密集台阵按一定半径划分为多个阵内台站和多个阵外台站。
参阅图2,在本公开实施例中,若密集线性台阵为台间距为d的N个台站组成,当以第n个台站为中心台站时,以r=md为半径作圆,圆内部的多个台站为多个阵内台站,圆外部的多个台站成为多个阵外台站,n=1,…,N。
S130,基于各中心台站对应的多个阵内台站采集的地震面波数据,计算各中心台站对应的阵内相移频散能量,以及,基于全部台站采集的地震面波数据,计算各中心台站对应的阵外相移频散能量。
根据Cheng等人(Multichannel analysis of passive surface waves based oncross-correlations.Geophysics,2016)提出的被动源面波频散曲线提取方法,计算面波的各频率谐波的阵内相移频散能量的公式为:
研究人员在实际研究中发现,由于被动源的噪声源分布并不理想,如图4所示,被动源面波的互相关信号常常出现较强的单支性,因此,将台站对地震面波数据的互相关函数正支和负支对称叠加可能会导致有效信号被无效信号“污染”,降低数据质量。因此,在本公开实施例中,对台站对地震记录的互相关函数正支、负支、对称叠加频散能量分别计算,从中选取质量最好的进行提取。
在本公开实施例中,具体的,步骤S130包括S131~S134。
S131,基于各中心台站对应的多个阵内台站采集的地震面波数据,利用互相关函数,计算与各中心台站对应的第一正支频散能量、第一负支频散能量和第一对称叠加频散能量。
令表示第一正支频散能量,表示第一负支频散能量,表示第一对称叠加频散能量,f表示地震面波的频率,CT表示地震面波的估计相速度,xkl表示第k个阵内台站和第l个阵内台站之间的距离,分别为k,l两台站记录的地震面波数据的互相关函数正支和负支的傅里叶变换,多个阵内台站为密集线性台阵内的第n-m至n+m个台站,则:
S132,从第一正支频散能量、第一负支频散能量和第一对称叠加频散能量中选出一种频散能量,作为中心台站的阵内相移频散能量。
在本公开实施例中,从第一正支频散能量、第一负支频散能量和第一对称叠加频散能量中选出质量最好的一种频散能量作为对应的中心台站的阵内相移频散能量。参阅图7,图7中的小图a)、c)、e)分别表示了利用阵内相移法对互相关函数正支、负支和对称叠加计算,得到阵内相移的第一正支频散能量、第一负支频散能量和第一对称叠加频散能量示意图,比较图a)、c)、e),可看出图a)表示的第一正支频散能量的图形效果比图c)表示的第一负支频散能量的图形效果好,图e)表示的第一对称叠加频散能量的图形效果与图a)相当,说明此时互相关函数的正支计算得到的频散能量质量更高,因此,选择第一正支频散能量作为此时的阵内相移频散能量。
S133,基于全部台站采集的地震面波数据,利用互相关函数,计算与各中心台站对应的第二正支频散能量、第二负支频散能量和第二对称叠加频散能量。
令表示第二正支频散能量,表示第二负支频散能量,表示第二对称叠加频散能量,f表示地震面波的频率,CT表示地震面波的估计相速度,xkl表示第k个阵外台站和第l个阵内台站之间的距离,分别为k,l两台站记录的地震面波数据的互相关函数正支和负支的傅里叶变换,密集线性台阵包括N个台站,多个阵内台站为密集线性台阵内的第n-m至n+m个台站,多个阵外台站为密集线性台阵内的除多个阵内台站以外的台站,则:
S134,从第二正支频散能量、第二负支频散能量和第二对称叠加频散能量中选出一种频散能量,作为中心台站的阵外相移频散能量。
在本公开实施例中,从第二正支频散能量、第二负支频散能量和第二对称叠加频散能量中选出质量最好的一种频散能量作为对应的中心台站的阵内相移频散能量。参阅图7,图7中的小图b)、d)、f)分别表示了利用阵外相移法对互相关函数正支、负支和对称叠加计算得到阵外相移的第二正支频散能量、第二负支频散能量和第二对称叠加频散能量示意图,比较图b)、d)、f),可看出图b)表示的第二正支频散能量的图形效果比图d)、f)表示的第二负支频散能量和第二对称叠加频散能量的图形效果好,说明互相关函数正支计算得到的频散能量质量更高,互相关函数负支计算得到的频散能量很差,将互相关函数对称叠加时由于负支的加入降低了频散能量的质量,低频频散能量受互相关函数负支的干扰明显,因此,选择第二正支频散能量作为此时的阵外相移频散能量。
S140,根据各中心台站对应的阵内相移频散能量和阵外相移频散能量,提取出对应的阵内频散曲线和阵外频散曲线。
在本公开实施例中,步骤S140包括步骤S141~S144。
步骤S141,针对各中心台站,获得不同频率的地震面波对应的阵内相移频散能量的最大值,以估计各频率的地震面波对应的第一真实相速度。
步骤S142,基于第一真实相速度和对应的地震面波的频率,得到各中心台站对应的阵内频散曲线。
具体的,步骤S141、S142是指,基于阵内相移频散能量,在保证连续性和合理性的前提下,提取频散质量较高的中高频段频散曲线,针对每一个目标频率点,选取该频率下频散能量最大值对应的相速度估计该频率的真实相速度,不同频率的相速度共同构成了频域内的阵内频散曲线,将其在周期域等间隔插值得到周期域的阵内频散曲线。
步骤S143,针对各中心台站,获得不同频率的地震面波对应的阵外相移频散能量的最大值,以估计各频率的地震面波对应的第二真实相速度;
步骤S144,基于第二真实相速度和对应的地震面波的频率,得到各中心台站对应的阵外频散曲线。
具体的,步骤S143、S144是指,基于阵外相移频散能量,在保证连续性和合理性的前提下,提取频散质量较高的中低频段频散曲线,针对每一个目标频率点选取该频率下频散能量最大值对应的相速度估计该频率的真实相速度,不同频率的相速度共同构成了频域内的阵外频散曲线,将其在周期域等间隔插值得到周期域的阵外频散曲线。
S150,给各中心台站对应的阵内频散曲线和阵外频散曲线在各频率点分别分配权重,并进行加权叠加计算,使阵内频散曲线和阵外频散曲线融合为一条宽频带频散曲线。
在本公开实施例中,阵内频散曲线为中高频频散曲线,阵外频散曲线为中低频频散曲线,在重叠频段分别施加给阵内频散曲线和阵外频散曲线不同的权重进行加权叠加得到频带更宽的频散曲线。具体步骤包括S151~S152。
步骤S151,比较中心台站对应的阵内频散曲线和阵外频散曲线,选取阵内频散曲线和阵外频散曲线最接近的频段。
参阅图3,图3示出了阵内频散曲线和阵内频散曲线,周期为0.2s~0.4s(即频率为2.5Hz~5Hz)频段内的阵内频散曲线和阵外频散曲线最接近,选择该频率范围内的阵内频散曲线和阵外频散曲线进行融合。
步骤S152,将频段内的阵内相移频散曲线和阵外相移频散曲线进行加权叠加,融合为一条宽频带频散曲线。其中,在高频段内,阵内频散曲线的权重更大,在低频段内,阵外频散曲线的权重更大。
令C(f)、Cint(f)、Cext(f)分别表示在宽频带频散曲线、阵内相移频散曲线、阵外相移频散曲线中频率为f的地震面波对应的相速度,ωint(f)表示阵内相移频散曲线中频率点f对应的相速度的权重,ωext(f)表示阵外相移频散曲线中频率点f对应的相速度权重,则:
C(f)=ωint(f)·Cint(f)+ωext(f)·Cext(f)。
在本公开实施例中,基于步骤S140,结合阵内相移频散能量和阵外相移频散能量各自的特点和优势,分别利用阵内频散能量、阵外相移频散能量提取中高频、中低频频散曲线,并选择了适当的重叠频域施以不同的权重进行加权叠加,得到一条连续的、频带更宽的频散曲线。这样,对地下结构进行反演时,能够保证在浅层具有较高分辨率的同时对深部结构具有更强的约束。
需要说明的是,在本公开实施例中,阵内相移法计算中面波信号的叠加次数Nint为:
阵外相移法计算中面波信号的叠加次数Next为:
Next=(N-2m-1)(2m+1);
假设Nint=Next,有:
m(2m+1)&=(N-2m-1)(2m+1);
则:
N=3m+1;
即当台站总数N=3m+1时,面波在阵内台站之间、阵外台阵与阵内台阵之间的相移的叠加次数相同。那么,当N>3m+1时,阵外相移叠加次数更高。在处理实际数据时,该条件容易得到满足,即在一般情况下阵外相移叠加次数高于阵内相移。如此,可在不降低对地下结构探测的分辨率的前提下,增加面波信号在台站之间的叠加次数,并且提高对相关面波信号的提取质量,使频散能量图表现出更高的稳定性和收敛性,更有利于低频信息的提取。
实施例1
在本实施例中,将本公开提供的一种宽频带频散曲线提取方法应用于湖南沃溪金矿成像的频散曲线提取工作中。为了研究沃溪矿区浅地表地质构造,技术人员在沃溪矿区东部布设了2条被动源密集测线,测线方位角为160°,在2019年6月到2019年7月进行了为期一个月的持续观测。台阵布设分布见图5,SmartSolo测线长约8km,台间距为100m,采集仪器为75个SmartSolo-3C(5Hz)节点式地震仪;Zland测线长约12km,台间距为100m,采集仪器为118个Zland-3C(5Hz)节点式地震仪。
计算采集获得地震面波数据的互相关函数,以SmartSolo测线的互相关结果为例进行展示,见图6。与图4类似,互相关函数对称性很差,正支信号质量较高,但是负支信号质量较差,表现出强烈的单边性。
为对比测试阵内相移(经典相移法)和阵外相移的低频频散曲线提取能力,以1500m为半径分别计算频散能量,见图7。图7中反映了很多重要信息,对比图7的a)、c)、e)可以看出,在此例中经典相移法利用互相关函数正支的计算结果明显优于负支,这与互相关信号(图6)表现出的特征一致,对称叠加的频散能量与正支相当,负支的加入并没有产生明显影响;图7的b)、d)、f)同样反映互相关函数正支优于负支的特点,但对称叠加时由于负支的加入降低了数据质量,低频频散信号受其干扰明显而得不到准确的频散信息。对比图7的a)和7b)可以发现,传统相移法低频只能提取到1.5Hz左右,而阵外相移包含有效信息的低频信号达到了约0.5Hz,转换到周期域内由0.7s扩展到了2s,这对频带扩展效果非常突出。此外,在高频段很容易发现经典相移法(阵内相移)的频散信号质量优于阵外相移,验证了经典相移法在高频信号处理时具有独到的优势,因此,需要计算阵内相移来补充高频信号。
分别利用阵内相移和阵外相移计算频散能量,在计算时为了更充分发挥各自的优势周期,阵内相移计算半径为600m,阵外相移计算半径为1500m。利用阵内相移提取中高频的频散曲线,用阵外相移提取中低频的频散曲线,并在周期域内等间隔插值获得等周期间隔的频散曲线。以图3所示的方案进行融合,图8以SmartSolo测线为例进行展示融合前后的频散曲线。图8的a)中黑线、灰线分别为阵内相移、阵外相移提取出的频散曲线,可以发现阵内相移(经典相移法)随周期增大相速度急剧上升,出现不真实的频散特性,反映出了传统相移法难以提取低频信号的问题。选择在周期为0.3s附近的一个频带内进行融合,高频部分以阵内相移频散曲线为主,低频部分以阵外相移频散曲线为主,融合后的频散曲线如图8的b)所示,可以看到提取的大多数低频信号都达到了周期为1.8s,部分质量较高的频散曲线能够进一步提取到周期为2s(频率为0.5Hz)的信号。
利用经典相移法(阵内相移)以1500m为半径对SmartSolo测线所有台站计算了频散能量并提取,来验证本发明提出的拓距相移法其中的阵外相移方法的有效性。图9为利用经典相移法对SmartSolo测线提取出的频散曲线,可以看到绝大多数频散曲线在0.8s左右截至,提取不出超过1s的频散信息。由此可见,与图8的b)用拓距相移法提取的频散曲线相比,传统相移法(阵内相移)缺乏提取低频频散曲线的能力。通过拓距相移法提取的频散曲线,低频端由数赫兹降到了约1.8s(0.55Hz),信噪比更高的可达到2s(0.5Hz;参见图10),将阵内相移和阵外相移相结合很好的拓宽了提取频散曲线的频带,可为相关技术人员更多的地下介质信息。
通过以上实例和相关对比结果证明,本发明提出的拓距相移法在保留高频频散信息的同时,大大提升了低频频散曲线的提取能力,拓宽了被动源面波勘探的频散曲线带宽,可以为技术人员提供了更多关于地下介质的信息,从而增大了探测的有效深度,对被动源面波勘探具有重大意义。
本公开还提供了一种宽频带频散曲线提取装置,包括:探测模块、台阵划分模块、频散能量计算模块、频散曲线提取模块、曲线融合模块。
探测模块,用于利用密集线性台阵内的各台站采集地震面波数据。
台阵划分模块,用于分别以密集台阵内的各个台站为中心台站,将密集台阵按一定半径划分为多个阵内台站和多个阵外台站。
频散能量计算模块,用于基于各中心台站对应的多个阵内台站采集的地震面波数据,计算各中心台站对应的阵内相移频散能量,以及,基于全部台站采集的地震面波数据,计算各中心台站对应的阵外相移频散能量。
频散曲线提取模块,用于根据各中心台站对应的阵内相移频散能量和阵外相移频散能量,提取出对应的阵内频散曲线和阵外频散曲线。
曲线融合模块,用于给各中心台站对应的阵内频散曲线和阵外频散曲线在各频率点分别分配权重,并进行加权叠加计算,使阵内频散曲线和阵外频散曲线融合为一条宽频带频散曲线。
本公开提供的一种宽频带频散曲线提取装置具有与本公开提供的一种宽频带频散曲线提取方法相同的作用和有益效果,在此不做赘述。
可以理解的是,台阵划分模块、探测模块、频散能量计算模块、频散能量融合模块、曲线生成模块可以合并在一个模块中实现,或者其中的任意一个模块可以被拆分成多个模块。或者,这些模块中的一个或多个模块的至少部分功能可以与其他模块的至少部分功能相结合,并在一个模块中实现。根据本发明的实施例,台阵划分模块、探测模块、频散能量计算模块、频散能量融合模块、曲线生成模块中的至少一个可以至少被部分地实现为硬件电路,例如现场可编程门阵列(FPGA)、可编程逻辑阵列(PLA)、片上系统、基板上的系统、封装上的系统、专用集成电路(ASIC),或可以以对电路进行集成或封装的任何其他的合理方式等硬件或固件来实现,或以软件、硬件以及固件三种实现方式的适当组合来实现。或者,台阵划分模块、探测模块、频散能量计算模块、频散能量融合模块、曲线生成模块中的至少一个可以至少被部分地实现为计算机程序模块,当该程序被计算机运行时,可以执行相应模块的功能。
本领域技术人员可以理解,本公开的各个实施例和/或权利要求中记载的特征可以进行多种组合或/或结合,即使这样的组合或结合没有明确记载于本公开中。特别地,在不脱离本公开精神和教导的情况下,本公开的各个实施例和/或权利要求中记载的特征可以进行多种组合和/或结合。所有这些组合和/或结合均落入本公开的范围。
尽管已经参照本公开的特定示例性实施例示出并描述了本公开,但是本领域技术人员应该理解,在不背离所附权利要求及其等同物限定的本公开的精神和范围的情况下,可以对本公开进行形式和细节上的多种改变。因此,本公开的范围不应该限于上述实施例,而是应该不仅由所附权利要求来进行确定,还由所附权利要求的等同物来进行限定。
Claims (10)
1.一种宽频带频散曲线提取方法,应用于密集线性台阵,其特征在于,包括:
利用所述密集线性台阵内的各台站采集地震面波数据;
分别以所述密集台阵内的各个台站为中心台站,将所述密集台阵按一定半径划分为多个阵内台站和多个阵外台站;
基于各所述中心台站对应的多个阵内台站采集的地震面波数据,计算各所述中心台站对应的阵内相移频散能量,以及,基于全部台站采集的地震面波数据,计算各所述中心台站对应的阵外相移频散能量;
根据各所述中心台站对应的所述阵内相移频散能量和所述阵外相移频散能量,提取出对应的阵内频散曲线和阵外频散曲线;
给各所述中心台站对应的所述阵内频散曲线和所述阵外频散曲线在各频率点分别分配权重,并进行加权叠加计算,使所述阵内频散曲线和所述阵外频散曲线融合为一条宽频带频散曲线。
2.根据权利要求1所述的方法,其特征在于,所述利用所述密集线性台阵内的各台站采集地震面波数据之后,方法还包括:
对所述地震面波数据分别进行单道预处理。
3.根据权利要求1所述的方法,其特征在于,所述分别以所述密集台阵内的各个台站为中心台站,将所述密集台阵按一定半径划分为多个阵内台站和多个阵外台站包括:
若所述密集线性台阵为台间距为d的N个台站组成,当以第n个台站为中心台站时,以r=md为半径作圆,所述圆内部的多个台站为所述多个阵内台站,所述圆外部的多个台站成为所述多个阵外台站,n=1,…,N。
4.根据权利要求1所述的方法,其特征在于,所述基于各所述中心台站对应的多个阵内台站采集的地震面波数据,计算各所述中心台站对应的阵内相移频散能量,以及,基于全部台站采集的地震面波数据,计算各所述中心台站对应的阵外相移频散能量,包括:
基于各所述中心台站对应的所述多个阵内台站采集的所述地震面波数据,利用互相关函数,计算与各所述中心台站对应的第一正支频散能量、第一负支频散能量和第一对称叠加频散能量;
从所述第一正支频散能量、第一负支频散能量和第一对称叠加频散能量中选出一种频散能量,作为所述中心台站的阵内相移频散能量;
基于全部台站采集的所述地震面波数据,利用互相关函数,计算与各所述中心台站对应的第二正支频散能量、第二负支频散能量和第二对称叠加频散能量;
从所述第二正支频散能量、第二负支频散能量和第二对称叠加频散能量中选出一种频散能量,作为所述中心台站的阵外相移频散能量。
6.根据权利要求4所述的方法,其特征在于,所述基于全部台站采集的所述地震面波数据,利用互相关函数,计算与各所述中心台站对应的第二正支频散能量、第二负支频散能量和第二对称叠加频散能量包括:
令表示所述第二正支频散能量,表示所述第二负支频散能量,表示所述第二对称叠加频散能量,f表示所述地震面波的频率,CT表示所述地震面波的估计相速度,xkl表示第k个所述阵外台站和第l个阵内台站之间的距离,分别为k,l两台站记录的地震面波数据的互相关函数正支和负支的傅里叶变换,所述密集线性台阵包括N个台站,所述多个阵内台站为所述密集线性台阵内的第n-m至n+m个台站,所述多个阵外台站为所述密集线性台阵内的除所述多个阵内台站以外的台站,则:
7.根据权利要求1所述的方法,其特征在于,所述根据各所述中心台站对应的所述阵内相移频散能量和所述阵外相移频散能量,提取出对应的阵内频散曲线和阵外频散曲线包括:
针对各所述中心台站,获得其采集到的不同频率的所述地震面波对应的所述阵内相移频散能量的最大值,以估计各频率的所述地震面波对应的第一真实相速度;
基于所述第一真实相速度和对应的所述地震面波的频率,得到各所述中心台站对应的所述阵内频散曲线;
针对各所述中心台站,获得其采集到的不同频率的所述地震面波对应的所述阵外相移频散能量的最大值,以估计各频率的所述地震面波对应的第二真实相速度;
基于所述第二真实相速度和对应的所述地震面波的频率,得到各所述中心台站对应的所述阵外频散曲线。
8.根据权利要求1所述的方法,其特征在于,所述给各所述中心台站对应的所述阵内频散曲线和所述阵外频散曲线在各频率点分别分配权重,并进行加权叠加计算,使所述阵内频散曲线和所述阵外频散曲线融合为一条宽频带频散曲线包括:
比较所述中心台站对应的阵内频散曲线和阵外频散曲线,选取阵内频散曲线和阵外频散曲线最接近的频段;
将所述频段内的所述阵内相移频散曲线和所述阵外相移频散曲线进行加权叠加,融合为一条宽频带频散曲线。
9.根据权利要求8所述的方法,其特征在于,所述将所述频段内的所述阵内相移频散曲线和所述阵外相移频散曲线进行加权叠加,融合为一条宽频带频散曲线包括:
令C(f)、Cint(f)、Cext(f)分别表示在所述宽频带频散曲线、阵内相移频散曲线、阵外相移频散曲线中频率为f的所述地震面波对应的相速度,ωint(f)表示所述阵内相移频散曲线中频率点f对应的相速度的权重,ωext(f)表示所述阵外相移频散曲线中频率点f对应的相速度权重,则:
C(f)=ωint(f)·Cint(f)+ωext(f)·Cext(f)。
10.一种宽频带频散曲线提取装置,其特征在于,包括:
探测模块,用于利用所述密集线性台阵内的各台站采集地震面波数据;
台阵划分模块,用于分别以所述密集台阵内的各个台站为中心台站,将所述密集台阵按一定半径划分为多个阵内台站和多个阵外台站;
频散能量计算模块,用于基于各所述中心台站对应的多个阵内台站采集的地震面波数据,计算各所述中心台站对应的阵内相移频散能量,以及,基于全部台站采集的地震面波数据,计算各所述中心台站对应的阵外相移频散能量;
频散曲线提取模块,用于根据各所述中心台站对应的所述阵内相移频散能量和所述阵外相移频散能量,提取出对应的阵内频散曲线和阵外频散曲线;
曲线融合模块,用于给各所述中心台站对应的所述阵内频散曲线和所述阵外频散曲线在各频率点分别分配权重,并进行加权叠加计算,使所述阵内频散曲线和所述阵外频散曲线融合为一条宽频带频散曲线。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010978546.XA CN112083487B (zh) | 2020-09-16 | 2020-09-16 | 宽频带频散曲线提取方法、装置 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010978546.XA CN112083487B (zh) | 2020-09-16 | 2020-09-16 | 宽频带频散曲线提取方法、装置 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN112083487A true CN112083487A (zh) | 2020-12-15 |
CN112083487B CN112083487B (zh) | 2021-12-14 |
Family
ID=73736779
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202010978546.XA Active CN112083487B (zh) | 2020-09-16 | 2020-09-16 | 宽频带频散曲线提取方法、装置 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN112083487B (zh) |
Cited By (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN116360369A (zh) * | 2023-03-29 | 2023-06-30 | 山东农业工程学院 | 一种陶瓷刀片智能配料控制方法及系统 |
CN116400406A (zh) * | 2023-04-21 | 2023-07-07 | 中国地震局地球物理研究所 | 一种基于阵列的被动源多模式面波频散曲线提取方法 |
CN116577829A (zh) * | 2023-05-15 | 2023-08-11 | 中国矿业大学(北京) | 一种基于背景噪音频散曲线自动化提取方法 |
CN116577829B (zh) * | 2023-05-15 | 2024-06-04 | 中国矿业大学(北京) | 一种基于背景噪音频散曲线自动化提取方法 |
Citations (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105785440A (zh) * | 2016-02-29 | 2016-07-20 | 河南理工大学 | 一种矿井槽波双分量地震信号频散曲线提取方法 |
CN106094020A (zh) * | 2016-05-31 | 2016-11-09 | 中国石油天然气集团公司 | 一种地震反演方法及装置 |
WO2016187252A1 (en) * | 2015-05-20 | 2016-11-24 | Conocophillips Company | Surface wave tomography using sparse data acquisition |
CN109477904A (zh) * | 2016-06-22 | 2019-03-15 | 休斯敦大学系统 | 地震或声波频散的非线性信号比较和高分辨率度量 |
CN109884709A (zh) * | 2019-04-01 | 2019-06-14 | 西安石油大学 | 一种基于面波旅行时层析的转换波静校正方法 |
CN109923440A (zh) * | 2017-10-12 | 2019-06-21 | 南方科技大学 | 面波勘探方法及终端设备 |
CN110568495A (zh) * | 2019-09-24 | 2019-12-13 | 中南大学 | 基于广义目标函数的瑞雷波多模式频散曲线反演方法 |
CN111103621A (zh) * | 2019-12-09 | 2020-05-05 | 东华理工大学 | 一种主动源共成像点叠加多道面波分析方法 |
-
2020
- 2020-09-16 CN CN202010978546.XA patent/CN112083487B/zh active Active
Patent Citations (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2016187252A1 (en) * | 2015-05-20 | 2016-11-24 | Conocophillips Company | Surface wave tomography using sparse data acquisition |
CN105785440A (zh) * | 2016-02-29 | 2016-07-20 | 河南理工大学 | 一种矿井槽波双分量地震信号频散曲线提取方法 |
CN106094020A (zh) * | 2016-05-31 | 2016-11-09 | 中国石油天然气集团公司 | 一种地震反演方法及装置 |
CN109477904A (zh) * | 2016-06-22 | 2019-03-15 | 休斯敦大学系统 | 地震或声波频散的非线性信号比较和高分辨率度量 |
CN109923440A (zh) * | 2017-10-12 | 2019-06-21 | 南方科技大学 | 面波勘探方法及终端设备 |
CN109884709A (zh) * | 2019-04-01 | 2019-06-14 | 西安石油大学 | 一种基于面波旅行时层析的转换波静校正方法 |
CN110568495A (zh) * | 2019-09-24 | 2019-12-13 | 中南大学 | 基于广义目标函数的瑞雷波多模式频散曲线反演方法 |
CN111103621A (zh) * | 2019-12-09 | 2020-05-05 | 东华理工大学 | 一种主动源共成像点叠加多道面波分析方法 |
Non-Patent Citations (2)
Title |
---|
孙天为 等: "利用密集台阵研究宾川沉积层及地壳结构", 《中国优秀硕士学位论文全文数据库》 * |
宋先海: "基于模式识别算法的高频瑞雷波频散曲线非线性反演研究", 《中国博士学位论文全文数据库》 * |
Cited By (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN116360369A (zh) * | 2023-03-29 | 2023-06-30 | 山东农业工程学院 | 一种陶瓷刀片智能配料控制方法及系统 |
CN116360369B (zh) * | 2023-03-29 | 2024-04-19 | 山东农业工程学院 | 一种陶瓷刀片智能配料控制方法及系统 |
CN116400406A (zh) * | 2023-04-21 | 2023-07-07 | 中国地震局地球物理研究所 | 一种基于阵列的被动源多模式面波频散曲线提取方法 |
CN116400406B (zh) * | 2023-04-21 | 2023-12-19 | 中国地震局地球物理研究所 | 一种基于阵列的被动源多模式面波频散曲线提取方法 |
CN116577829A (zh) * | 2023-05-15 | 2023-08-11 | 中国矿业大学(北京) | 一种基于背景噪音频散曲线自动化提取方法 |
CN116577829B (zh) * | 2023-05-15 | 2024-06-04 | 中国矿业大学(北京) | 一种基于背景噪音频散曲线自动化提取方法 |
Also Published As
Publication number | Publication date |
---|---|
CN112083487B (zh) | 2021-12-14 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN106094029B (zh) | 利用偏移距矢量片地震数据预测储层裂缝的方法 | |
US7725265B2 (en) | VH signal integration measure for seismic data | |
Mi et al. | Near-surface imaging from traffic-induced surface waves with dense linear arrays: An application in the urban area of Hangzhou, China | |
CN104280775B (zh) | 一种基于全波形矢量偏移叠加的微地震监测定位方法 | |
US20080288173A1 (en) | Seismic attributes for reservoir localization | |
CN112083487B (zh) | 宽频带频散曲线提取方法、装置 | |
Cornou et al. | Contribution of dense array analysis to the identification and quantification of basin-edge-induced waves, Part II: Application to Grenoble basin (French Alps) | |
EP2419761A1 (en) | Interferometric seismic data processing | |
CN112883564B (zh) | 一种基于随机森林的水体温度预测方法及预测系统 | |
CN112285767B (zh) | 海底地震仪四分量海洋面波多阶频散能量成像装置及方法 | |
CN106970417B (zh) | 椭圆展开转换波速度分析方法与系统 | |
CN107728214A (zh) | 一种裂缝预测方法 | |
CN104570116A (zh) | 基于地质标志层的时差分析校正方法 | |
Ning et al. | High-frequency surface-wave imaging from traffic-induced noise by selecting in-line sources | |
CN105137479B (zh) | 一种面元覆盖次数的计算方法及装置 | |
Dewangan et al. | Nature of the ambient noise, site response, and orientation of ocean‐bottom seismometers (OBSs): Scientific results of a passive seismic experiment in the Andaman Sea | |
Boaga et al. | Shear wave structural models of Venice Plain, Italy, from time cross correlation of seismic noise | |
Xu et al. | Optimized workflows for high-frequency seismic interferometry using dense arrays | |
CN107918142A (zh) | 一种地震勘探方法 | |
CN110850474A (zh) | 一种海洋地震资料震源鬼波压制方法和系统 | |
CN114415234B (zh) | 基于主动源面波频散和h/v确定浅地表横波速度的方法 | |
CN111691876B (zh) | 一种利用声波测井对邻井成像的方法、装置及存储介质 | |
CN110967751B (zh) | 基于地面浅井监测的微地震事件的定位方法及存储介质 | |
CN108802822B (zh) | 方位各向异性介质中的保幅直接叠前时间偏移方法及装置 | |
Kuyuk | On the use of Stockwell transform in structural dynamic analysis |
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 |