CN110967734B - 基于快速傅立叶变换的虚源重构方法及系统 - Google Patents
基于快速傅立叶变换的虚源重构方法及系统 Download PDFInfo
- Publication number
- CN110967734B CN110967734B CN201811139507.XA CN201811139507A CN110967734B CN 110967734 B CN110967734 B CN 110967734B CN 201811139507 A CN201811139507 A CN 201811139507A CN 110967734 B CN110967734 B CN 110967734B
- Authority
- CN
- China
- Prior art keywords
- virtual source
- fourier transform
- cross
- fast fourier
- correlation
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 22
- 230000005284 excitation Effects 0.000 claims abstract description 26
- 238000004364 calculation method Methods 0.000 claims abstract description 12
- 238000001228 spectrum Methods 0.000 claims description 36
- 238000005070 sampling Methods 0.000 claims description 24
- 125000004122 cyclic group Chemical group 0.000 claims description 13
- 238000010586 diagram Methods 0.000 description 14
- 230000003111 delayed effect Effects 0.000 description 3
- 238000002679 ablation Methods 0.000 description 1
- 150000001875 compounds Chemical class 0.000 description 1
- 230000001934 delay Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 238000005305 interferometry Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000004044 response 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
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
公开了一种基于快速傅立叶变换的虚源重构方法及系统。该方法可以包括:步骤1:输入地震数据,确定虚源A与虚源检波器B;步骤2:利用快速傅立叶变换,分别计算在每一个实际震源激发下对应的A与B的波场记录的互相关值;步骤3:将所有的互相关值进行叠加,得到A点激发B点接收的虚源波场;步骤4:重复进行步骤1‑3,完成所有虚源波场的计算。本发明通过快速傅立叶变换在频率域实现了互相关运算,极大提高了虚源重构的效率;频率域的互相关运算可以很方便地拓展为反褶积或互相干处理。
Description
技术领域
本发明涉及地震勘探技术领域,更具体地,涉及一种基于快速傅立叶变换的虚源重构方法及系统。
背景技术
虚源重构(地震干涉法)是一种通过互相关(反褶积、互相干)进行波场重构的方法,可以总结为“互相关叠加”处理:对两个接收器记录的波场进行互相关,然后将围绕检波器的所有震源的互相关结果进行叠加,就可以得到两个接收器之间的波场响应,其中一个接收器就称为虚源。
虚源重构的核心是波场之间的互相关,也是计算量最大的地方。通常是在时间域进行互相关运算,其原因是时间域进行互相关物理意义非常明确,易于理解。虽然虚源重构的原理是较为简单的,但在时间域实现的计算量是非常大的,远远不能满足理论研究与实际生产的要求。因此,有必要开发一种基于快速傅立叶变换的虚源重构方法及系统。
公开于本发明背景技术部分的信息仅仅旨在加深对本发明的一般背景技术的理解,而不应当被视为承认或以任何形式暗示该信息构成已为本领域技术人员所公知的现有技术。
发明内容
本发明提出了一种基于快速傅立叶变换的虚源重构方法及系统,其能够通过快速傅立叶变换在频率域实现了互相关运算,极大提高了虚源重构的效率;频率域的互相关运算可以很方便地拓展为反褶积或互相干处理。
根据本发明的一方面,提出了一种基于快速傅立叶变换的虚源重构方法。所述方法可以包括:步骤1:输入地震数据,确定虚源A与虚源检波器B;步骤2:利用快速傅立叶变换,分别计算在每一个实际震源激发下对应的A与B的波场记录的互相关值;步骤3:将所有的互相关值进行叠加,得到A点激发B点接收的虚源波场;步骤4:重复进行步骤1-3,完成所有虚源波场的计算。
优选地,所述步骤2包括:根据每一个实际震源激发下对应的A与B的波场记录,获得时间序列x与h;确定x与h的采样点数均为M,构造时间序列其中,采样点数均为N=2*M;对所述时间序列进行快速傅立叶变换得到频谱将的振幅谱相乘,相位谱相减,得到互相关的频谱对所述频谱进行反快速傅立叶变换得到循环相关结果根据所述循环相关结果计算A与B的波场记录的互相关值。
其中,xi为时间序列x中的第i个参量。
其中,hi为时间序列h中的第i个参量。
优选地,通过公式(3)计算所述A与B的波场记录的互相关值:
根据本发明的另一方面,提出了一种基于快速傅立叶变换的虚源重构系统,其特征在于,该系统包括:存储器,存储有计算机可执行指令;处理器,所述处理器运行所述存储器中的计算机可执行指令,执行以下步骤:步骤1:输入地震数据,确定虚源A与虚源检波器B;步骤2:利用快速傅立叶变换,分别计算在每一个实际震源激发下对应的A与B的波场记录的互相关值;步骤3:将所有的互相关值进行叠加,得到A点激发B点接收的虚源波场;步骤4:重复进行步骤1-3,完成所有虚源波场的计算。
优选地,所述步骤2包括:根据每一个实际震源激发下对应的A与B的波场记录,获得时间序列x与h;确定x与h的采样点数均为M,构造时间序列其中,采样点数均为N=2*M;对所述时间序列进行快速傅立叶变换得到频谱将的振幅谱相乘,相位谱相减,得到互相关的频谱对所述频谱进行反快速傅立叶变换得到循环相关结果根据所述循环相关结果计算A与B的波场记录的互相关值。
其中,xi为时间序列x中的第i个参量。
其中,hi为时间序列h中的第i个参量。
优选地,通过公式(3)计算所述A与B的波场记录的互相关值:
本发明的方法和装置具有其它的特性和优点,这些特性和优点从并入本文中的附图和随后的具体实施方式中将是显而易见的,或者将在并入本文中的附图和随后的具体实施方式中进行详细陈述,这些附图和具体实施方式共同用于解释本发明的特定原理。
附图说明
通过结合附图对本发明示例性实施例进行更详细的描述,本发明的上述以及其它目的、特征和优势将变得更加明显,其中,在本发明示例性实施例中,相同的参考标号通常代表相同部件。
图1示出了根据本发明的基于快速傅立叶变换的虚源重构方法的步骤的流程图。
图2示出了根据本发明的一个实施例的雷克子波及延迟雷克子波的示意图。
图3示出了根据本发明的一个实施例的两个时间序列在时间域与频率域的互相关结果的示意图。
图4示出了根据本发明的一个实施例的层状速度模型的示意图。
图5a和图5b分别示出了根据本发明的一个实施例的全波场的共炮点道集与下行波场的共炮点道集的示意图。
图6a和图6b分别示出了根据本发明的一个实施例的直达波的共接收点道集与自由表面下行波的共接收点道集的示意图。
图7a和图7b分别示出了根据本发明的一个实施例的时间域进行虚源重构的虚源炮集与频率域进行虚源重构的虚源炮集的示意图。
图8示出了根据本发明的一个实施例的地表x=200m处的正演模拟炮集的示意图。
具体实施方式
下面将参照附图更详细地描述本发明。虽然附图中显示了本发明的优选实施例,然而应该理解,可以以各种形式实现本发明而不应被这里阐述的实施例所限制。相反,提供这些实施例是为了使本发明更加透彻和完整,并且能够将本发明的范围完整地传达给本领域的技术人员。
图1示出了根据本发明的基于快速傅立叶变换的虚源重构方法的步骤的流程图。
在该实施例中,根据本发明的基于快速傅立叶变换的虚源重构方法可以包括:步骤1:输入地震数据,确定虚源A与虚源检波器B;步骤2:利用快速傅立叶变换,分别计算在每一个实际震源激发下对应的A与B的波场记录的互相关值;步骤3:将所有的互相关值进行叠加,得到A点激发B点接收的虚源波场;步骤4:重复进行步骤1-3,完成所有虚源波场的计算。
在一个示例中,步骤2包括:根据每一个实际震源激发下对应的A与B的波场记录,获得时间序列x与h;确定x与h的采样点数均为M,构造时间序列其中,采样点数均为N=2*M;对时间序列进行快速傅立叶变换得到频谱将的振幅谱相乘,相位谱相减,得到互相关的频谱对频谱进行反快速傅立叶变换得到循环相关结果根据循环相关结果计算A与B的波场记录的互相关值。
其中,xi为时间序列x中的第i个参量。
其中,hi为时间序列h中的第i个参量。
在一个示例中,通过公式(3)计算A与B的波场记录的互相关值:
具体地,虚源重构在时间域可以表示公式(4):
其中,ukA(t)表示A点记录的震源k激发的地震道,ukA(-t)表示地震道ukA(t)的逆时记录,叠加是沿N个震源进行的,DAB(t)表示虚源A激发而B点记录的地震道,*表示褶积,表示互相关。为了提高计算效率,将虚源重构在频率域进行表示为公式(5):
其中,UkA(f)表示ukA(t)的频谱,表示UkA(f)的共轭复数。在频率域进行虚源重构可以提高计算效率的本质原因是利用了快速傅立叶变换(FFT),采用这种算法使离散傅立叶变换的运算量减少了几个数量级,特别是采样点数越多,FFT算法节省的计算量就越大。
根据本发明的基于快速傅立叶变换的虚源重构方法可以包括:
步骤1:输入地震数据,确定虚源A与虚源检波器B。
步骤2:利用快速傅立叶变换,分别计算在每一个实际震源激发下对应的A与B的波场记录的互相关值:根据每一个实际震源激发下对应的A与B的波场记录,获得时间序列x与h;确定时间序列x与h的采样点数均为M0,如果M0≠2k,通过补零使得采样点数变为M,确保M=2k≥M0。通过公式(1)计算时间序列通过公式(2)计算时间序列其中,采样点数均为N=2*M;对时间序列进行快速傅立叶变换得到频谱将的振幅谱相乘,相位谱相减,得到互相关的频谱对频谱进行反快速傅立叶变换得到循环相关结果根据循环相关结果通过公式(3)计算A与B的波场记录的互相关值,包括正延迟与负延迟。
步骤3:将所有的互相关值进行叠加,得到A点激发B点接收的虚源波场。
步骤4:重复进行步骤1-3,完成所有虚源波场的计算。
本方法通过快速傅立叶变换在频率域实现了互相关运算,极大提高了虚源重构的效率;频率域的互相关运算可以很方便地拓展为反褶积或互相干处理。
为便于理解本发明实施例的方案及其效果,以下给出两个具体应用示例。本领域技术人员应理解,该示例仅为了便于理解本发明,其任何具体细节并非意在以任何方式限制本发明。
应用示例1
图2示出了根据本发明的一个实施例的雷克子波及延迟雷克子波的示意图。
图3示出了根据本发明的一个实施例的两个时间序列在时间域与频率域的互相关结果的示意图。
利用FFT实现两个时间序列的互相关。首先通过一个简单的算例验证利用FFT实现互相关的正确性,时间序列uA为雷克子波,序列uB为延迟的100ms的雷子子波,且最大振幅设为0.5,如图2所示。分别在时间域与频率域通过FFT进行时间序列uA与uB的互相关运算,结果如图3所示,可以看出利用FFT实现互相关的结果与时间域互相关结果完全一致,证实了利用FFT在频率域实现互相关的可行性。
应用示例2
以VSP波场的虚源重构为例。通过虚源重构的方法将VSP多次波重构为地表一次反射波,即常规地表激发地表接收的观测方式,可以极大拓展VSP波场的照明范围。
图4示出了根据本发明的一个实施例的层状速度模型的示意图。
图5a和图5b分别示出了根据本发明的一个实施例的全波场的共炮点道集与下行波场的共炮点道集的示意图。
考虑如图4所示的层状速度模型,模型大小为1500*1000m,井位于x=0m处。151个震源沿地表激发,间距10m;47个接收器沿井均匀布置,间距15m,第一个接收器的坐标为z=200m;采样点数2001,采样间隔0.5ms。图5a为共炮点道集记录,图5b为对应的下行波记录,在下行波记录中利用切除将直达波和自由地表下行波进行分离。
图6a和图6b分别示出了根据本发明的一个实施例的直达波的共接收点道集与自由表面下行波的共接收点道集的示意图。
图7a和图7b分别示出了根据本发明的一个实施例的时间域进行虚源重构的虚源炮集与频率域进行虚源重构的虚源炮集的示意图。
图8示出了根据本发明的一个实施例的地表x=200m处的正演模拟炮集的示意图。
将分离后的共炮点道集分选为共接收点道集(即IVSP的共炮点波场),图6a、6b分别为井中第一个接收器的直达波与自由地表下行波的共接收点道集,虚源位置为地表x=200m处。图7a、7b分别为利用时间域与频率域互相关进行虚源重构得到的虚源炮集,虚源位置为地表x=200m处,可以看出两者是完全一致的,效率大约提高了5倍,而且本例的采样点数只有2001,当采样点数越多时基于FFT进行虚源重构的计算效率也越高。图8为对应的地表x=200m处的正演模拟炮集,对比可以看出虚源炮集反射波的走时特征大部分得到准确重构,部分道由于孔径不足使得反射走时不能准确重构。
综上所述,本发明通过快速傅立叶变换在频率域实现了互相关运算,极大提高了虚源重构的效率;频率域的互相关运算可以很方便地拓展为反褶积或互相干处理。
本领域技术人员应理解,上面对本发明的实施例的描述的目的仅为了示例性地说明本发明的实施例的有益效果,并不意在将本发明的实施例限制于所给出的任何示例。
根据本发明的实施例,提供了一种基于快速傅立叶变换的虚源重构系统,其特征在于,该系统包括:存储器,存储有计算机可执行指令;处理器,所述处理器运行所述存储器中的计算机可执行指令,执行以下步骤:步骤1:输入地震数据,确定虚源A与虚源检波器B;步骤2:利用快速傅立叶变换,分别计算在每一个实际震源激发下对应的A与B的波场记录的互相关值;步骤3:将所有的互相关值进行叠加,得到A点激发B点接收的虚源波场;步骤4:重复进行步骤1-3,完成所有虚源波场的计算。
在一个示例中,步骤2包括:根据每一个实际震源激发下对应的A与B的波场记录,获得时间序列x与h;确定x与h的采样点数均为M,构造时间序列其中,采样点数均为N=2*M;对时间序列进行快速傅立叶变换得到频谱将的振幅谱相乘,相位谱相减,得到互相关的频谱对频谱进行反快速傅立叶变换得到循环相关结果根据循环相关结果计算A与B的波场记录的互相关值。
其中,xi为时间序列x中的第i个参量。
其中,hi为时间序列h中的第i个参量。
在一个示例中,通过公式(3)计算A与B的波场记录的互相关值:
本系统通过快速傅立叶变换在频率域实现了互相关运算,极大提高了虚源重构的效率;频率域的互相关运算可以很方便地拓展为反褶积或互相干处理。
本领域技术人员应理解,上面对本发明的实施例的描述的目的仅为了示例性地说明本发明的实施例的有益效果,并不意在将本发明的实施例限制于所给出的任何示例。
以上已经描述了本发明的各实施例,上述说明是示例性的,并非穷尽性的,并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下,对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。
Claims (8)
1.一种基于快速傅立叶变换的虚源重构方法,其特征在于,包括:
步骤1:输入地震数据,确定虚源A与虚源检波器B;
步骤2:利用快速傅立叶变换,分别计算在每一个实际震源激发下对应的A与B的波场记录的互相关值;
步骤3:将所有的互相关值进行叠加,得到A点激发B点接收的虚源波场;
步骤4:重复进行步骤1-3,完成所有虚源波场的计算;
其中,所述步骤2包括:
根据每一个实际震源激发下对应的A与B的波场记录,获得时间序列x与h;
5.一种基于快速傅立叶变换的虚源重构系统,其特征在于,该系统包括:
存储器,存储有计算机可执行指令;
处理器,所述处理器运行所述存储器中的计算机可执行指令,执行以下步骤:
步骤1:输入地震数据,确定虚源A与虚源检波器B;
步骤2:利用快速傅立叶变换,分别计算在每一个实际震源激发下对应的A与B的波场记录的互相关值;
步骤3:将所有的互相关值进行叠加,得到A点激发B点接收的虚源波场;
步骤4:重复进行步骤1-3,完成所有虚源波场的计算;
其中,所述步骤2包括:
根据每一个实际震源激发下对应的A与B的波场记录,获得时间序列x与h;
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811139507.XA CN110967734B (zh) | 2018-09-28 | 2018-09-28 | 基于快速傅立叶变换的虚源重构方法及系统 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811139507.XA CN110967734B (zh) | 2018-09-28 | 2018-09-28 | 基于快速傅立叶变换的虚源重构方法及系统 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN110967734A CN110967734A (zh) | 2020-04-07 |
CN110967734B true CN110967734B (zh) | 2022-03-08 |
Family
ID=70027719
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201811139507.XA Active CN110967734B (zh) | 2018-09-28 | 2018-09-28 | 基于快速傅立叶变换的虚源重构方法及系统 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN110967734B (zh) |
Families Citing this family (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111538081B (zh) * | 2020-06-05 | 2021-05-25 | 吉林大学 | 一种地震数据初至波的外推方法 |
CN112817047B (zh) * | 2020-12-31 | 2021-10-08 | 北京东方联创地球物理技术有限公司 | 海洋地震自适应去鬼波方法、装置、电子设备及介质 |
CN113740906A (zh) * | 2021-10-19 | 2021-12-03 | 中国地质大学(北京) | 一种水下垂直缆地震波干涉成像方法及装置 |
Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102141431A (zh) * | 2010-02-01 | 2011-08-03 | 鸿远亚太科技(北京)有限公司 | 双层介质空间中的声场测量与变换方法 |
CN103713312A (zh) * | 2012-10-09 | 2014-04-09 | 中国石油化工股份有限公司 | 一种虚源地震观测系统的设计方法 |
CN103728013A (zh) * | 2013-12-25 | 2014-04-16 | 广西科技大学 | 噪声源识别方法 |
CN104345343A (zh) * | 2014-12-01 | 2015-02-11 | 中国海洋石油总公司 | 一种复杂海底相关的层间多次波预测方法 |
CN106291684A (zh) * | 2015-06-06 | 2017-01-04 | 中国石油化工股份有限公司 | 一种盲源地震波场的地震响应恢复与虚源道集构建方法 |
US9697620B2 (en) * | 2013-03-28 | 2017-07-04 | Reeves Wireline Technologies Limited | Borehole log data processing methods |
CN107402405A (zh) * | 2016-05-18 | 2017-11-28 | 中国石油化工股份有限公司 | 静相位虚源道集构建方法 |
-
2018
- 2018-09-28 CN CN201811139507.XA patent/CN110967734B/zh active Active
Patent Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102141431A (zh) * | 2010-02-01 | 2011-08-03 | 鸿远亚太科技(北京)有限公司 | 双层介质空间中的声场测量与变换方法 |
CN103713312A (zh) * | 2012-10-09 | 2014-04-09 | 中国石油化工股份有限公司 | 一种虚源地震观测系统的设计方法 |
US9697620B2 (en) * | 2013-03-28 | 2017-07-04 | Reeves Wireline Technologies Limited | Borehole log data processing methods |
CN103728013A (zh) * | 2013-12-25 | 2014-04-16 | 广西科技大学 | 噪声源识别方法 |
CN104345343A (zh) * | 2014-12-01 | 2015-02-11 | 中国海洋石油总公司 | 一种复杂海底相关的层间多次波预测方法 |
CN106291684A (zh) * | 2015-06-06 | 2017-01-04 | 中国石油化工股份有限公司 | 一种盲源地震波场的地震响应恢复与虚源道集构建方法 |
CN107402405A (zh) * | 2016-05-18 | 2017-11-28 | 中国石油化工股份有限公司 | 静相位虚源道集构建方法 |
Non-Patent Citations (3)
Title |
---|
Surface Wave Scattering at Vertical Discontinuties;M. Munasinghe et al.;《IEEE》;20051205;第267-270页 * |
虚源法地震技术及数值模型试验;陈国金 等;《地球物理学进展》;20131031;第28卷(第5期);第2725-2733页 * |
虚源法相关道集噪声压制技术;雷朝阳 等;《地球物理学进展》;20170331;第32卷(第3期);第1154-1160页 * |
Also Published As
Publication number | Publication date |
---|---|
CN110967734A (zh) | 2020-04-07 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Zhang et al. | A stable and practical implementation of least-squares reverse time migration | |
US8437998B2 (en) | Hybrid method for full waveform inversion using simultaneous and sequential source method | |
US8688381B2 (en) | Simultaneous source inversion for marine streamer data with cross-correlation objective function | |
Boonyasiriwat et al. | An efficient multiscale method for time-domain waveform tomography | |
RU2612896C2 (ru) | Ортогональное кодирование источника и приемника | |
US9625593B2 (en) | Seismic data processing | |
US20120275267A1 (en) | Seismic Data Processing | |
Yang et al. | Time-domain least-squares migration using the Gaussian beam summation method | |
US20130311149A1 (en) | Tomographically Enhanced Full Wavefield Inversion | |
US9766357B2 (en) | Seismic image dip decomposition estimation and recomposition | |
US7639564B2 (en) | 3-D TAU-P interpolation | |
Kim et al. | Frequency-domain reverse-time migration with source estimation | |
US20120051176A1 (en) | Reverse time migration back-scattering noise removal using decomposed wavefield directivity | |
CN110967734B (zh) | 基于快速傅立叶变换的虚源重构方法及系统 | |
WO2007143355A2 (en) | Diplet-based seismic processing | |
US9348050B2 (en) | Near-surface noise prediction and removal for data recorded with simultaneous seismic sources | |
RU2570827C2 (ru) | Гибридный способ для полноволновой инверсии с использованием способа одновременных и последовательных источников | |
Araujo et al. | Symplectic scheme and the Poynting vector in reverse-time migration | |
Sun et al. | Alternating first-arrival traveltime tomography and waveform inversion for near-surface imaging | |
Vigh et al. | Comparisons for waveform inversion, time domain or frequency domain? | |
CN105572735B (zh) | 一种提高叠前深度偏移成像精度的方法及装置 | |
Yu et al. | Application of weighted early-arrival waveform inversion to shallow land data | |
EP2496968B1 (en) | System and method for seismic beam formation that accounts for equipment misalignment | |
US10338253B2 (en) | Method of suppressing spectral artefacts of wavefield decomposition caused by imperfect extrapolation | |
Witten et al. | Signal-to-noise estimates of time-reverse images |
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 |