CN110618459B - 地震数据处理方法和装置 - Google Patents
地震数据处理方法和装置 Download PDFInfo
- Publication number
- CN110618459B CN110618459B CN201911061919.0A CN201911061919A CN110618459B CN 110618459 B CN110618459 B CN 110618459B CN 201911061919 A CN201911061919 A CN 201911061919A CN 110618459 B CN110618459 B CN 110618459B
- Authority
- CN
- China
- Prior art keywords
- wave
- wave field
- target
- seismic data
- wavefield
- 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
- 238000003672 processing method Methods 0.000 title claims abstract description 14
- 239000013598 vector Substances 0.000 claims abstract description 107
- 238000000034 method Methods 0.000 claims abstract description 60
- 230000010287 polarization Effects 0.000 claims abstract description 59
- 238000012545 processing Methods 0.000 claims abstract description 35
- 230000008569 process Effects 0.000 claims abstract description 28
- 238000003384 imaging method Methods 0.000 claims description 63
- 230000002441 reversible effect Effects 0.000 claims description 31
- 238000003860 storage Methods 0.000 claims description 7
- 238000004590 computer program Methods 0.000 claims description 3
- 238000006073 displacement reaction Methods 0.000 description 28
- 239000002245 particle Substances 0.000 description 27
- 238000000926 separation method Methods 0.000 description 21
- 230000005012 migration Effects 0.000 description 17
- 238000013508 migration Methods 0.000 description 17
- 238000006243 chemical reaction Methods 0.000 description 15
- 238000010586 diagram Methods 0.000 description 11
- 238000007781 pre-processing Methods 0.000 description 11
- 230000000694 effects Effects 0.000 description 10
- 230000006870 function Effects 0.000 description 9
- 238000004364 calculation method Methods 0.000 description 8
- 230000015572 biosynthetic process Effects 0.000 description 7
- 238000005516 engineering process Methods 0.000 description 7
- 239000000126 substance Substances 0.000 description 7
- 238000000605 extraction Methods 0.000 description 6
- 230000009466 transformation Effects 0.000 description 5
- 230000008859 change Effects 0.000 description 4
- 238000004458 analytical method Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 239000011159 matrix material Substances 0.000 description 2
- 230000003287 optical effect Effects 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- 230000009471 action Effects 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 239000013078 crystal Substances 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 230000004907 flux Effects 0.000 description 1
- 238000009776 industrial production Methods 0.000 description 1
- 229910052500 inorganic mineral Inorganic materials 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 239000011707 mineral Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 230000010363 phase shift Effects 0.000 description 1
- 230000000644 propagated effect Effects 0.000 description 1
- 230000003252 repetitive effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000001629 suppression 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/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
- G01V1/362—Effecting static or dynamic corrections; Stacking
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
技术领域
本申请实施例涉及地震技术领域，尤其涉及一种地震数据处理方法和装置。
背景技术
随着油气和矿产资源的不断开发，地震勘探逐渐向多元化、精细化发展。近些年来，伴随着多分量采集技术的推广和高性能计算(Highperformance computing，HPC)技术的成熟，弹性波逆时偏移技术越来越受到工业界和科研界的重视。弹性波逆时偏移技术可以处理包括反射波，折射波，多次波，回转波等在内的全波场信息，在复杂地质构造上成像效果突出，不受倾角限制。此外，弹性波偏移处理可以利用转换横波的信息，对气云，岩丘和断层下部构造进行良好照明和成像，为岩性识别和油藏检测等提供了直接的依据。
各向异性弹性波逆时偏移的必要步骤包括：震源波场正向延拓，检波器波场反向延拓，震源波场和检波器波场的波场分离(得到准P波和准S波)，纯波型波场的互相关成像以及角道集的抽取。目前，已有的弹性波逆时偏移技术利用的方向矢量是代表弹性波能流密度的波印廷矢量，该方法主要通过对各向异性介质传播速度参数以及各向异性参数的设置，建立群角与极化角之间的转换关系；应用弹性波方程对成像区域内的多分量地震数据进行波场延拓，得到震源正传波场和炮记录反传波场，以及相应的波印廷矢量；利用波印廷矢量和群角与极化角之间的转换关系，进行波场分离；最后利用多分量地震数据中对应的波场分离结果进行各向异性介质中弹性波矢量成像。
在上述现有的技术方案中，在波场分离时进行的群角到极化角的转化需要先通过群角转化为相角，再通过相角获得极化角，过程繁琐，效率低下。
发明内容
本申请实施例提供一种地震数据处理方法和装置，用于解决波场分离时极化角获取过程繁琐，效率低下的问题。
第一方面，本申请实施例提供一种地震数据处理方法，包括：
获取目标区域的地震数据；
根据所述目标区域的地震数据，对震源子波进行正向延拓，获得正向延拓后的边界波场值；
根据所述边界波场值进行反向延拓，以获得震源波场，以及根据所述地震数据进行反向延拓，以获得检波器波场；
在获得目标波场的过程中，在预设时间间隔到达时，根据目标波场的波场分量、所述波场分量的时间一阶导数和空间一阶导数，获得波矢量，以及根据所述波矢量，获得所述目标波场的极化方向；
其中，所述波矢量用于表示所述目标波场的相速度方向，所述目标波场为所述震源波场或检波器波场。
可选地，所述波场分量包括：水平方向上的波场分量和垂直方向上的波场分量。
可选地，根据所述目标波场的波场分量、所述波场分量的时间一阶导数和空间一阶导数，获得波矢量，包括：
根据如下公式，获得波矢量；
其中，vx表示波场在水平方向上的分量，vz表示波场在垂直方向上的分量，p表示波矢量。
可选地，所述方法还包括：
根据所述目标波场的极化方向，对波场分量分离，获得所述目标波场的准P波和准S波。
可选地，所述方法还包括：
根据所述波矢量，获得张角或地层倾角。
可选地，所述方法还包括：
根据每个预设时间间隔获得的所述张角或地层倾角，对每个预设时间间隔的所述震源波场的准P波、准S波与所述检波器波场的准P波、准S波进行互相关运算；
根据所有时间间隔获得的互相关运算结果，获得单炮偏移成像剖面。
可选地，所述方法还包括：
根据所有预设时间间隔获得的所述张角，提取出共成像点张角道集；或者，
根据所有预设时间间隔获得的地层倾角，提取出共成像点倾角道集。
可选地，所述地震数据为多炮地震数据中任一炮地震数据；
所述方法还包括：在获得所述多炮震数据分别对应的单炮偏移成像剖面后，将所有单炮偏移成像剖面叠加以形成所述目标区域的偏移成像剖面。
可选地，所述获取目标区域的地震数据，包括：
获取目标区域的原始地震数据；
对所述地震数据进行预处理，获得所述目标区域的地震数据；
所述预处理包括如下至少一项：去除直达波，面波，随机噪音以及去除子波效应。
可选地，所述波场分量是速度在水平方向上的分量和垂直方向上的分量，或者位移在水平方向上的分量和垂直方向上的分量。
第二方面，本申请实施例提供一种地震数据处理装置，包括
第一获取模块，用于获取目标区域的地震数据；
第一处理模块，用于根据所述目标区域的地震数据，对震源子波进行正向延拓，获得正向延拓后的边界波场值；根据所述边界波场值进行反向延拓，以获得震源波场，以及根据所述地震数据进行反向延拓，以获得检波器波场；
第二获取模块，用于在获得目标波场的过程中，在预设时间间隔到达时，根据目标波场的波场分量、所述波场分量的时间一阶导数和空间一阶导数，获得波矢量；
第二处理模块，用于根据所述波矢量，获得所述目标波场的极化方向；
其中，所述波矢量用于表示所述目标波场的相速度方向，所述目标波场为所述震源波场或检波器波场。
可选地，所述波场分量包括：水平方向上的波场分量和垂直方向上的波场分量。
可选地，所述第二获取模块，具体用于：
其中，vx表示波场在水平方向上的分量，vz表示波场在垂直方向上的分量，p表示波矢量。
可选地，所述第二处理模块，还用于：根据所述目标波场的极化方向，对波场分量分离，获得所述目标波场的准P波和准S波。
可选地，所述第二处理模块，还用于：
根据所述波矢量，获得张角或地层倾角。
可选地，所述第二处理模块，还用于：
根据每个预设时间间隔获得的所述张角或地层倾角，对每个预设时间间隔的所述震源波场的准P波、准S波与所述检波器波场的准P波、准S波进行互相关运算；
根据所有时间间隔获得的互相关运算结果，获得单炮偏移成像剖面。
可选地，所述第二处理模块，还用于：
根据所有预设时间间隔获得的所述张角，提取出共成像点张角道集；或者，
根据所有预设时间间隔获得的地层倾角，提取出共成像点倾角道集。
可选地，所述地震数据为多炮地震数据中任一炮地震数据；
所述第二处理模块，还用于：在获得所述多炮地震数据分别对应的单炮偏移成像剖面后，将所有单炮偏移成像剖面叠加以形成所述目标区域的偏移成像剖面。
可选地，所述第一获取模块，具体用于：
获取目标区域的原始地震数据；
对所述地震数据进行预处理，获得所述目标区域的地震数据；
所述预处理包括如下至少一项：去除直达波，面波，随机噪音以及去除子波效应。
可选地，所述波场分量是速度在水平方向上的分量和垂直方向上的分量，或者位移在水平方向上的分量和垂直方向上的分量。
第三方面，本申请实施例提供一种电子设备，包括：
存储器，用于存储程序指令；
处理器，用于调用并执行所述存储器中的程序指令，执行如本申请第一方面所述的地震数据处理方法。
第四方面，本申请实施例提供一种计算机可读存储介质，所述计算机存储介质存储有计算机程序，所述计算机程序被处理器执行时实现本申请第一方面所述的地震数据处理方法。
本申请实施例提供一种地震数据处理方法和装置。通过获取目标区域的地震数据；再根据所述目标区域的地震数据，获取目标波场；在获得目标波场的过程中，在预设时间间隔到达时，根据目标波场的波场分量、所述波场分量的时间一阶导数和空间一阶导数，获得波矢量，通过所述波矢量获得所述目标波场的极化方向，简化了获得极化方向的过程，提高了地震数据处理效率。
附图说明
图1为本申请一实施例提供的群角、相角以及极化角的转化关系的示意图；
图2为本申请一实施例提供的一种地震数据处理方法流程示意图；
图3为本申请一实施例提供的地层模型及其背景速度模型的参数示意图；
图4为本申请一实施例提供的张角或地层倾角的关系示意图；
图5为本申请一实施例提供的地震数据处理成像剖面示意图；
图6为本申请一实施例提供的地震数据处理装置的结构示意图；
图7为本申请一实施例提供的电子设备的结构示意图。
具体实施方式
各向异性广泛存在于地球介质中，形成机制包括：晶体各向异性，直接的应力作用以及岩性各向异性，其在地震波上的表现主要有传播速度(与传播方向有关)、体波间相互耦合及横波发生分裂等现象。忽略各向异性，将导致成像结果同相轴的错位，对进一步的处理和解释造成负面影响。研发各向异性下的弹性波逆时偏移技术，无论是对能源工业生产还是地球内部结构成像都有着十分重要的意义。
各向异性弹性波逆时偏移的必要步骤包括：震源波场正向延拓，检波器波场反向延拓，震源波场和检波器波场的波场分离(得到准P波和准S波)，纯波型波场的互相关成像以及角道集的抽取。其中，为了减少存储量，可以先进行震源波场正向延拓并将边界值存储下来，利用完备的边界条件反传来重构震源波场。这样可以同时反向延拓震源波场和检波器波场，方便了同一时刻纯波型波场互相关的计算。为了减少计算量，可以在满足奈奎斯特定理的采样间隔下等间隔地进行波场分离和互相关运算，例如每隔7个时间采样点进行一次分离和成像。
各向异性弹性波逆时偏移面临的主要问题是波场分离，压制低频噪音和共成像点角道集的抽取。关于各向异性弹性波场的分离，传统方法往往涉及波数域内的计算，需要大量的傅里叶变换，因此计算成本很高且占用内存巨大。或者利用拟散度旋度算子，一定程度上可以减少运算量和存储，但还是会破坏地震波的振幅和相位，需要额外的矫正。基于极化角的空间分离公式是一种简单高效的方案，但是极化方向的求取需要利用其他方向矢量转化。低频噪音主要是由强反射界面的反射波引起，在频率域表现为低频特性，一种流行做法是使用拉普拉斯算子对叠加的成像结果进行滤波，去除低频信号从而压制噪音，但是这样会改变成像轴的相位。利用低频噪音在张角域表现出的大角度特征，可以通过限定入射波和出射波之间的张角范围在成像过程中消减噪音的影响。逆时偏移角道集的提取有诸多方法，包括方向矢量法，局部平面波分解方法，局部相移方法等。在所有的策略中，使用方向矢量是成本较低，效果较好的一类方法。由以上讨论可知，最优的一种解决方案是利用波的传播方向信息来进行弹性波逆时偏移。
目前，已有的方向矢量类弹性波逆时偏移技术利用的方向矢量是代表弹性波能流密度的波印廷矢量。这一矢量对应的是群速度的方向，也就是地震波射线路径的切线方向，它与垂直方向的夹角被称为群角。空间域波场分离方案是基于极化方向的，极化方向与垂直方向的夹角被称为极化角，所以在波场分离之前要进行群角到极化角的转化。群角到极化角的关系是通过相角建立的，相角是地震波同相位面外法线的方向与垂直方向的夹角。因此转化实质上是通过这样的关系建立的：从群角到相角，从相角到极化角。从相角到极化角，考虑到分离准S波所用的极化方向垂直于准S波的质点振动方向，准P波和准S波都采用准P波的转换公式。但是从群角到相角，准P波和准S波的转化关系并不一致，即同一个群角，对应于准P波和准S波有不同的相角。其次，不同于从相角到群角有显式的函数关系，从群角到相角的过程是很难用函数直接表达出来的，因此需要数值求解，建表一一对应，这就需要函数满足单调变化的条件，但是对于准S波在某些各向异性参数条件下并不满足单调变化的条件。因此，这个分离方案有三个缺陷：1.需要两次分离，先按照准P波求取的极化方向分离，将准P波的能量提取出来后，再将其他区域按照准S波的极化方向分离；2.从群角到极化角的计算是经过查表和插值实现的，这个过程同样可能引起误差；3.建表查表占用内存，并且需要一定的输入输出操作。
正确的张角或者地层倾角计算是基于相角信息完成的，因此基于传统波印廷矢量的方法直接利用群角求取入射波与出射波之间的张角或者地层倾角是不合理的。一般需要将群角到相角的转化表也储存下来，一一查表转化完成后再进行角道集的提取，或者在提取之后对结果进行校正。
通常情况下，弹性波介质中的波印廷矢量指的是波传播的能流密度，是应力张量和速度矢量矩阵乘积的相反数，数学上可以表示成公式一：
s＝-τ·ν公式一
其中，
表示波印廷矢量，τ表示应力张量，v表示速度矢量。图1为本申请一实施例提供的相角到极化角的转化关系的示意图，如图1所示，在这一定义下，波印廷矢量与波的射线路径相切，对应的是群角。但是在地震数据处理过程中，无论是波场分离，还是角道集的提取，所使用的都是波矢量或者是相角。波矢量沿同相位面的外法线方向，其与垂直方向的夹角被称作相角。对于各向同性介质，能流密度矢量与波矢量沿同一方向，可以相互代替。但是对于各向异性介质，二者并不重合且从能流密度方向转化到波矢量的方向并不容易。首先，转化的过程十分繁琐。二维情况下，群角和相角之间满足下列关系，如公式二：
其中，ψ表示群角，θ表示相角，Vn表示准P或者准S的相速度，是关于模型参数和相角的复杂函数。上述公式二不易直接表达的反函数对应于群角到相角的转换，一般做法是采用建表查表插值的方法加以解决。因为描述这一关系的参数较多，包括模型的纵横波速度比，Thmosen参数ε，δ，以及群角本身。显然对于四维表格，无论是存储还是读取都降低了程序的效率。其次，方案具有多解性。对于准P波还有准S波二者的公式是不一样的，同一群角可以对应两个不同的相角。
本申请可以解决上述至少一个问题，具体方案如下所述。
图2为本申请一实施例提供的地震数据处理方法的流程示意图，如图2所示，本实施例的方法可以包括：
S210、获取目标区域的地震数据。
本实施例中，先获取目标区域的地震数据，其中，可以是实施本方案的设备接收其它设备发送的地震数据，或者，也可以是实施本方案的设备通过存储设备(例如U盘等)获取的地震数据。
可选地，获取目标区域的地震数据的一种可能的实现方式为：获取目标区域的原始地震数据；然后对获取的目标区域的原始地震数据进行预处理，所述原始地震数据预处理包括去除直达波，面波，随机噪音以及去除子波效应等。
可选地，所述地震数据为多炮地震数据中任一炮地震数据。例如：获取通过测线采集到的多分量地震数据，并对该多分量地震数据按炮号规则排列，获得多炮原始地震数据，针对每炮原始地震数据，通过上述预处理可以获得相应的每炮地震数据。其中，炮号规则可以参见相关技术中的描述，此处不再赘述。
S220、根据所述目标区域的地震数据，对震源子波进行正向延拓，获得正向延拓后的边界波场值。
本实施例中，在获得目标区域的地震数据之后，根据该地震数据获得震源子波，例如：对所述的目标区域的地震数据互相关运算，获得震源子波。然后根据地震数据，对震源子波进行正向延拓，获得正向延拓后的边界波场值，例如：根据背景速度模型，对震源子波进行正向延拓，获得正向延拓后的边界波场值，其中，背景速度模型是根据地震资料速度建模得到的地震波的传播速度和各向异性参数。其中，速度建模可以参见相关技术的描述，此处不再赘述。图3为本申请一实施例提供的地层模型及其背景速度模型的参数示意图，如图3所示，其中，VP表示纵波传播速度，VS表示横波传播速度，ε，δ用来衡量介质的各向异性参数。
其中，所述的边界波场值是指以质点振动的形式传播出的位移或者速度值。
需要说明的是，正向延拓的实现过程可以参见相关技术中的描述，此处不再赘述。
S230、根据所述边界波场值进行反向延拓，以获得震源波场，以及根据所述地震数据进行反向延拓，以获得检波器波场。
本实施例中，根据S220中获得的所述边界波场值进行反向延拓，以获得震源波场，以及根据所述目标区域的地震数据进行反向延拓，以获得检波器波场。
其中，反向延拓的实现过程可以参见相关技术中的描述，此处不再赘述。
S240、在获得目标波场的过程中，在预设时间间隔到达时，根据目标波场的波场分量、所述波场分量的时间一阶导数和空间一阶导数，获得波矢量，以及根据所述波矢量，获得所述目标波场的极化方向。
其中，所述波矢量用于表示所述目标波场的相速度方向，所述目标波场为所述震源波场或检波器波场。
本实施例中，在获得所述震源波场的过程中，在预设时间间隔到达时，根据所述震源波场更新的波场分量、所述震源波场分量的时间一阶导数和空间一阶导数，获得代表所述震源波场的相速度方向的波矢量，以及根据所述波矢量，获得所述震源波场的极化方向。
在获得所述检波器波场的过程中，在预设时间间隔到达时，根据所述检波器波场更新的波场分量、所述检波器波场分量的时间一阶导数和空间一阶导数，获得代表所述检波器波场的相速度方向的波矢量，以及根据所述波矢量，获得所述检波器波场的极化方向。
本实施例中，通过获取的所述目标区域的地震数据，进行相应的处理之后以获取震源波场或检测器波场，在获取震源波场或检测器波场的过程中，按照预设时间间隔，根据所述震源波场或检测器波场的波场分量、所述震源波场或检测器波场的所述波场分量的时间一阶导数和空间一阶导数，获得震源波场或检测器波场的波矢量，通过波矢量获得所述震源波场或检测器波场的极化方向，简化了获得极化方向的过程，提高了地震数据处理效率。
在一些实施例中，所述波场分量包括：水平方向上的波场分量和垂直方向上的波场分量。可选地，所述波场分量是速度在水平方向上的分量和垂直方向上的分量，或者位移在水平方向上的分量和垂直方向上的坐标分量。
在一些实施例中，上述根据所述目标波场的波场分量、所述波场分量的时间一阶导数和空间一阶导数，获得波矢量的一种可能的实现方式为：
根据如下公式三，获得波矢量；
其中，p表示所述目标波场的波矢量，vx表示质点速度或位移在水平方向上的坐标分量，vz表示质点速度或位移在垂直方向上的坐标分量，
分别表示质点速度或位移在对应空间方向上的一阶偏导数，
分别表示质点速度或位移在对应时间方向上的一阶偏导数。
其中，可以根据平面波理论和光流问题的求解方法，获得上述公式三，具体过程如下：
在获得目标波场的过程中(弹性波数值模拟速度会更新)，利用获取的所述目标波场的波场分量的时间一阶导数计算乘积并相加，获得所述目标波场的波场分量时间偏导数
如公式四所示。
其中，
表示所述目标波场的坐标分量时间偏导数，
分别表示质点速度或位移在对应时间方向上的一阶偏导数，vx表示质点速度或位移在水平方向上的坐标分量，vz表示质点速度或位移在垂直方向上的坐标分量。
在获得目标波场的过程中(弹性波数值模拟应力会更新)，利用获取的所述目标波场的波场分量的空间一阶导数计算乘积并相加，获得目标波场的波场分量空间梯度
如公式五所示。
其中，
分别表示质点速度或位移在对应空间方向上的一阶偏导数，
表示所述质点速度或位移模平方的空间梯度，vx表示质点速度或位移在水平方向上的坐标分量，vz表示质点速度或位移在垂直方向上的坐标分量。
获取上述
和
的乘积的相反数得到波矢量，如公式六所示：
其中，P表示所述目标波场的某一原始物理量，在声波方程中，P一般指代声压。我们选用的质点速度或位移绝对值的平方
也可以选用质点速度或位移的坐标分量vx，vz，或者是质点速度或位移的绝对值||v||，或者是原始波场经过散度或者旋度算子处理后得到的初始分离拟纵波波场uqP0和拟横波波场uqS0。选用质点速度模绝对值的平方，其包含的物理意义是弹性波传播时的质点动能。
从公式六可以看出，求取所述震源波场和检波器波场波矢量所需要的是：所述震源波场和检波器波场分量本身vx，vz，它们对时间的一阶偏导数
对空间的一阶偏导数
而上述这些参数都是波场迭代更新已知的，减少了震源波场和检波器波场波矢量求取的计算量，并且这种方法既适应于基于规则网格的有限差分方法，也适用于基于非规则，非结构化网格的有限元法、谱元法以及格子法。
可选地，为了提高计算结果的稳定性，对计算的结果求取局部加权平均意义下的波矢量，如公式七所示。
其中，Ω表示待求点的邻域，包括待求点本身及距离它最近的八个点，i点是该邻域内的一点，i的范围是1-9，αi表示该点的权重，可以设置成与距离有关的高斯函数，也可以直接设置为1，
表示所述目标波场的波矢量局部加权平均。
其中，波矢量表示相速度方向的方式为：对于逆时偏移而言，因为用到波矢量的方向信息，所以一般归一化之后用相角θ表示，如公式八所示：
其中，θ表示相角，
表示波矢量的水平分量，
表示波矢量的垂直分量。
在一些实施例中，根据所述波矢量，获得所述目标波场的极化方向的一种可能的实现方式，如下所述：
图1为本申请一实施例提供的相角到极化角的转化关系的示意图，如图1所示，根据所述的目标区域模型参数(纵波和横波速度参数以及各向异性参数)完成从所述波矢量方向到极化方向的转化。对于各向同性介质，所述波矢量所代表的相速度方向与极化方向重合，因此转化操作可以省略。在各向异性中，相角θ和极化角υ满足Christoffel方程，如公式九：
其中，
代表Christoffel张量，ρ代表介质的密度，Vn(n＝qP,qSV)表示不同波形的相速度，
是极化矢量，c11，c13，c44代表介质的弹性刚度系数，表示介质应力随应变的变化，I是单位矩阵，θ表示相角。
在垂直对称横向各项异性(Vertical Transverse Isotropic，VTI)介质中，四个刚度系数由公式十给出：
其中，ρ是介质的密度，VP0和VS0分别是准P波和准S波沿介质对称轴方向的传播速度，ε是度量准P波的各向异性强度参数，ε越大，介质的纵波各向异性越大，ε＝0，纵波无各向异性。δ是连接VP0和准P波是沿垂直于对称轴方向的速度VP90的一种过渡性参数，f是中间参数，可以方便计算。
解公式九，就可以推导出极化角υ和相角θ之间的关系：
其中，υ表示极化角，sgn()是符号函数，返回参数的符号，G11，G33和G13被称为Christoffel元。
公式十一中G11，G33和G13可由公式十二求解，
其中，G11为Christoffel元，c11，c44，c13，c33代表介质的弹性刚度系数，表示介质应力随应变的变化，θ表示相角。
显然，求出(c11-c44)sin2θ-(c33-c44)cos2θ]和(c13+c44)sinθcosθ后就可以避免许多重复运算，使得从相角到极化角的转化过程简单而高效。
对于对称轴不沿垂直方向的各向异性介质(Tilted Transversely Isotropic，TTI)，需要在求取到相角的基础上再加上倾斜角，在求取到极化角之后再减去倾斜角得到最终所用的极化角。
在一些实施例中，如图2所示，还可以包括：
S250：根据所述目标波场的极化方向，对所述目标波场的波场分量分离，获得所述目标波场的准P波和准S波。
本实施例中，根据所述震源波场的极化方向，对所述震源波场的波场分量分离，获得所述震源波场的准P波和准S波。根据所述检波器波场的极化方向，对所述检波器波场的波场分量分离，获得所述检波器波场的准P波和准S波。
其中，目标波场分离如公式十三所示：
其中，vx，vz是波场的坐标分量，UqP，UqS是波场的分离结果，υ是极化角。公式十三可以视作将初始波场投影到极化矢量方向和极化矢量的垂直方向，从而得到分离结果。
因此，通过上述方式获得的目标波场的准P波和准S波，不仅快速而且准确。
在一些实施例中，如图2所示，本实施例还可以包括：
S260：根据所述波矢量，获得张角或地层倾角。
如图2所示，利用所述震源波场和检波器波场的波矢量，可以得到所述震源波场的入射波的入射方向和所述检波器波场的反射波的出射方向之间的张角。所述震源波场正向延拓可以得到入射波的传播方向
所述检波器波场反向延拓可以得到反射波的出射方向
则入射波与反射波之间的张角α满足公式十四：
其中，
代表所述震源波场正向延拓得到的入射波的传播方向的波矢量，
代表所述检波器波场反向延拓得到的反射波的出射方向的波矢量，α代表入射波与反射波之间的张角。
地层倾角β满足公式十五：
其中，β代表地层倾角，
代表所述震源波场波矢量在X轴方向的坐标分量，
代表所述震源波场波矢量在Z轴方向的坐标分量，
代表所述检波器波场波矢量在水平方向上的坐标分量，
代表所述检波器波场波矢量在垂直方向上的坐标分量。图4为本申请一实施例提供的张角或地层倾角的关系示意图，如图4所示。
正确的张角或者地层倾角计算是基于相角信息完成的，因此使用波矢量的方向求取所述震源波场的入射波的入射方向和所述检波器波场的反射波的出射方向之间的张角或者地层倾角相较于群速度的方向更加合理。
在一些实施例中，如图2所示，本实施例的方法还可以包括：
S270、根据每个预设时间间隔获得的所述张角或地层倾角，对每个预设时间间隔的所述震源波场的准P波、准S波与所述检波器波场的准P波、准S波进行互相关运算；
S280、根据所有时间间隔获得的互相关运算结果，获得单炮偏移成像剖面。
上述步骤S270中，根据所述震源波场的波矢量(也就是入射波的波矢量)和所述检波器波场的波矢量(也就是出射波的波矢量)，计算二者之间的所述张角或地层倾角的大小，所述震源波场分离出的纯波型波场(准P波、准S波)和所述检波器波场分离出的纯波型波场(准P波、准S波)进行互相关运算，将所有时间间隔获得的计算值叠加在一起，得到最终的单炮弹性波逆时偏移剖面。可选地，在计算值叠加的时候，若某一时间间隔获得所述张角或地层倾角不满足预设角度范围，则针对张角或地层倾角的计算值不用于叠加处理，使得获得的单炮偏移成像剖面的成像效果更好。
在一些实施例中，本实施例的方法还可以包括：根据所有预设时间间隔获得的所述张角，提取出共成像点张角道集；或者，根据所有预设时间间隔获得的地层倾角，提取出共成像点倾角道集。
本实施例中，根据所述震源波场的入射波波矢量和所述检波器波场的出射波的波矢量计算二者之间的所述张角或地层倾角的大小，所述震源波场分离出的纯波型波场(准P波、准S波)和所述检波器波场分离出的纯波型波场(准P波、准S波)进行互相关运算，将结果放置在成像点的张角或地层倾角的位置，最后将所有时刻的计算值叠加在一起，得到最终的单炮弹性波逆时偏移张角道集或者倾角道集。
共成像点张角道集的抽取用公式十六表示：
其中，Psrc表示所述震源波场的入射P波，Prec表示所述检波器波场的出射P波，Srec表示所述检波器波场的出射S波，
表示成像点，t是成像时刻，Ipp表示震源波场分离出的准P波与检波器波场分离出的准P波成像得到的张角道集，Ips表示震源波场分离出的拟P波与检波器波场分离出的准S波成像得到的张角道集，Mpp表示对应两种角道集的角度投影函数，目的是根据求得的张角将互相关值投影到指定的位置上。
在公式十六中，Mpp可以由公式十七求出。
其中，
代表所述震源波场正向延拓得到的入射波的传播方向的波矢量，
代表所述检波器波场反向延拓得到的反射波的出射方向的波矢量，α代表入射波与反射波之间的张角，δ表示脉冲函数只有在
时，δ函数值为1，其他情况下其值为0。
类似地，可以得到共成像点倾角道集。
逆时偏移成像剖面则是共成像点角道集所有角度上值的叠加，即在每一个成像点上将一定角度范围内的值求和，可以表示为：
其中，
表示成像点，Ipp表示震源波场分离出的准P波与检波器波场分离出的准P波成像得到的张角道集，Ips表示震源波场分离出的准P波与检波器波场分离出的准S波成像得到的张角道集，
表示震源波场分离出的准P波与检波器波场分离出的准P波成像得到的张角道集在所有角度上值的叠加，即震源波场分离出的准P波与检波器波场分离出的准P波成像剖面，
表示震源波场分离出的准P波与检波器波场分离出的准S波成像得到的张角道集在所有角度上值的叠加，即震源波场分离出的准P波与检波器波场分离出的准S波成像剖面。
在一些实施例中，所述地震数据处理方法还包括：在获得所述多炮震数据分别对应的单炮偏移成像剖面后，将所有单炮偏移成像剖面叠加以形成所述目标区域的偏移成像剖面。图5为本申请一实施例提供的地震数据处理成像剖面示意图，如图5所示，(a)为PP成像剖面(震源波场P波与检波器波场P波互相关运算所得)；(b)为PS成像剖面(震源波场P波与检波器波场S波互相关运算所得)。
使用所述目标波场的波矢量方向作为张角或者地层倾角提取出的共成像点角道集更适合振幅随角度变化(Amplitude Versus Angle，AVA)分析，偏移速度分析(MigrationVelocity Analysis，MVA)等，并且不需要后续的转换，简化了过程。
图6为本申请一实施例提供的地震数据处理装置的结构示意图，如图6所示，本实施例的装置400可以包括：第一获取模块410，第一处理模块420，第二获取模块430和第二处理模块440。
第一获取模,410，用于获取目标区域的地震数据；
第一处理模块420，用于根据所述目标区域的地震数据，对震源子波进行正向延拓，获得正向延拓后的边界波场值；根据所述边界波场值进行反向延拓，以获得震源波场，以及根据所述地震数据进行反向延拓，以获得检波器波场；
第二获取模块430，用于在获得目标波场的过程中，在预设时间间隔到达时，根据目标波场的波场分量、所述波场分量的时间一阶导数和空间一阶导数，获得波矢量；
第二处理模块440，用于根据所述波矢量，获得所述目标波场的极化方向；
其中，所述波矢量用于表示所述目标波场的相速度方向，所述目标波场为所述震源波场或检波器波场。
在一些实施例中，所述波场分量为：水平方向上的波场分量和垂直方向上的波场分量。
在一些实施例中，所述第二获取模块430，具体用于：
其中，p表示所述目标波场的波矢量，vx表示质点速度或位移在水平方向上的坐标分量，vz表示质点速度或位移在垂直方向上的坐标分量，
分别表示质点速度或位移在对应空间方向上的一阶偏导数，
分别表示质点速度或位移在对应时间方向上的一阶偏导数。
在一些实施例中，所述第二处理模块440，还用于：根据所述目标波场的极化方向，对波场分量分离，获得所述目标波场的准P波和准S波。
在一些实施例中，所述第二处理模块440，还用于：
根据所述波矢量，获得张角或地层倾角。
在一些实施例中，所述第二处理模块440，还用于：
根据每个预设时间间隔获得的所述张角或地层倾角，对每个预设时间间隔的所述震源波场的准P波、准S波与所述检波器波场的准P波、准S波进行互相关运算；
根据所有时间间隔获得的互相关运算结果，获得单炮偏移成像剖面。
在一些实施例中，所述第二处理模块440，还用于：
根据所有预设时间间隔获得的所述张角，提取出共成像点张角道集；或者，
根据所有预设时间间隔获得的地层倾角，提取出共成像点倾角道集。
在一些实施例中，所述地震数据为多炮地震数据中任一炮地震数据；
所述第二处理模块440，还用于：在获得所述多炮震数据分别对应的单炮偏移成像剖面后，将所有单炮偏移成像剖面叠加以形成所述目标区域的偏移成像剖面。
在一些实施例中，所述第一获取模块410，具体用于：
获取目标区域的原始地震数据；
对所述地震数据进行预处理，获得所述目标区域的地震数据；
所述预处理包括如下至少一项：去除直达波，面波，随机噪音以及去除子波效应。
在一些实施例中，所述波场分量是速度在水平方向上的分量和垂直方向上的分量，或者位移在水平方向上的分量和垂直方向上的分量。
本实施例的装置，可以用于执行上述各方法实施例的技术方案，其实现原理和技术效果类似，此处不再赘述。
图7为本申请一实施例提供的电子设备的结构示意图，如图7所示，本实施例的电子设备500可以包括：存储器510、处理器520。
存储器510，用于存储程序指令；
处理器520，用于调用并执行所述存储器中的程序指令，执行：
获取目标区域的地震数据；
根据所述目标区域的地震数据，对震源子波进行正向延拓，获得正向延拓后的边界波场值；
根据所述边界波场值进行反向延拓，以获得震源波场，以及根据所述地震数据进行反向延拓，以获得检波器波场；
在获得目标波场的过程中，在预设时间间隔到达时，根据目标波场的波场分量、所述波场分量的时间一阶导数和空间一阶导数，获得波矢量，以及根据所述波矢量，获得所述目标波场的极化方向；
其中，所述波矢量用于表示所述目标波场的相速度方向，所述目标波场为所述震源波场或检波器波场。
在一些实施例中，所述波场分量为：水平方向上的波场分量和垂直方向上的波场分量。
在一些实施例中，所述处理器520，在根据所述目标波场的波场分量、所述波场分量的时间一阶导数和空间一阶导数，获得波矢量时，具体用于：
其中，p表示所述目标波场的波矢量，vx表示质点速度或位移在水平方向上的坐标分量，vz表示质点速度或位移在垂直方向上的坐标分量，
分别表示质点速度或位移在对应空间方向上的一阶偏导数，
分别表示质点速度或位移在对应时间方向上的一阶偏导数。
在一些实施例中，所述处理器520还用于：根据所述目标波场的极化方向，对波场分量分离，获得所述目标波场的准P波和准S波。
在一些实施例中，所述处理器520还用于：根据所述目标波场的极化方向，对波场分量分离，获得所述目标波场的准P波和准S波。
在一些实施例中，所述处理器520还用于：根据每个预设时间间隔获得的所述张角或地层倾角，对每个预设时间间隔的所述震源波场的准P波、准S波与所述检波器波场的准P波、准S波进行互相关运算；
根据所有时间间隔获得的互相关运算结果，获得单炮偏移成像剖面。
在一些实施例中，所述处理器520还用于：根据所有预设时间间隔获得的所述张角，提取出共成像点张角道集；或者，
根据所有预设时间间隔获得的地层倾角，提取出共成像点倾角道集。
在一些实施例中，所述地震数据为多炮地震数据中任一炮地震数据；
所述处理器520还用于：在获得所述多炮震数据分别对应的单炮偏移成像剖面后，将所有单炮偏移成像剖面叠加以形成所述目标区域的偏移成像剖面。
在一些实施例中，所述处理器520具体用于：
获取目标区域的原始地震数据；
对所述地震数据进行预处理，获得所述目标区域的地震数据；
所述预处理包括如下至少一项：去除直达波，面波，随机噪音以及去除子波效应。
在一些实施例中，所述波场分量是速度在水平方向上的分量和垂直方向上的分量，或者位移在水平方向上的分量和垂直方向上的分量。
本实施例的电子设备，可以用于执行上述各方法实施例的技术方案，其实现原理和技术效果类似，此处不再赘述。
本领域普通技术人员可以理解：实现上述各方法实施例的全部或部分步骤可以通过程序指令相关的硬件来完成。前述的程序可以存储于一计算机可读取存储介质中。该程序在执行时，执行包括上述各方法实施例的步骤；而前述的存储介质包括：只读内存(Read-Only Memory，ROM)、随机存取存储器(Random Access Memory，RAM)、磁碟或者光盘等各种可以存储程序代码的介质。
最后应说明的是：以上各实施例仅用以说明本申请的技术方案，而非对其限制；尽管参照前述各实施例对本申请进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本申请各实施例技术方案的范围。
Claims (7)
1.一种地震数据处理方法，其特征在于，包括：
获取目标区域的地震数据；
根据所述目标区域的地震数据，对震源子波进行正向延拓，获得正向延拓后的边界波场值；
根据所述边界波场值进行反向延拓，以获得震源波场，以及根据所述地震数据进行反向延拓，以获得检波器波场；
在获得目标波场的过程中，在预设时间间隔到达时，根据目标波场的波场分量、所述波场分量的时间一阶导数和空间一阶导数，获得波矢量，以及根据所述波矢量，获得所述目标波场的极化方向；
其中，所述波矢量用于表示所述目标波场的相速度方向，所述目标波场为所述震源波场或检波器波场；
根据所述目标波场的极化方向，对波场分量分离，获得所述目标波场的准P波和准S波；
根据所述目标波场的波场分量、所述波场分量的时间一阶导数和空间一阶导数，获得波矢量，包括：
根据如下公式，获得波矢量；
其中，vx表示波场在水平方向上的分量，vz表示波场在垂直方向的分量，p表示波矢量。
2.根据权利要求1所述的方法，其特征在于，还包括：
根据所述波矢量，获得张角或地层倾角。
3.根据权利要求2所述的方法，其特征在于，还包括：
根据每个预设时间间隔获得的所述张角或地层倾角，对每个预设时间间隔的所述震源波场的准P波、准S波与所述检波器波场的准P波、准S波进行互相关运算；
根据所有时间间隔获得的互相关运算结果，获得单炮偏移成像剖面。
4.根据权利要求2所述的方法，其特征在于，还包括：
根据所有预设时间间隔获得的所述张角，提取出共成像点张角道集；或者，
根据所有预设时间间隔获得的地层倾角，提取出共成像点倾角道集。
5.一种地震数据处理装置，其特征在于，包括：
第一获取模块，用于获取目标区域的地震数据；
第一处理模块，用于根据所述目标区域的地震数据，对震源子波进行正向延拓，获得正向延拓后的边界波场值；根据所述边界波场值进行反向延拓，以获得震源波场，以及根据所述地震数据进行反向延拓，以获得检波器波场；
第二获取模块，用于在获得目标波场的过程中，在预设时间间隔到达时，根据目标波场的波场分量、所述波场分量的时间一阶导数和空间一阶导数，获得波矢量；
第二处理模块，用于根据所述波矢量，获得所述目标波场的极化方向；
其中，所述波矢量用于表示所述目标波场的相速度方向，所述目标波场为所述震源波场或检波器波场；
第二处理模块，还用于根据所述目标波场的极化方向，对波场分量分离，获得所述目标波场的准P波和准S波；
所述第二获取模块，具体用于：
根据如下公式，获得波矢量；
其中，vx表示波场在水平方向上的分量，vz表示波场在垂直方向的分量，p表示波矢量。
6.一种电子设备，其特征在于，包括：
存储器，用于存储程序指令；
处理器，用于调用并执行所述存储器中的程序指令，执行如权利要求1至4中任一项所述的地震数据处理方法。
7.一种计算机可读存储介质，其特征在于，所述计算机存储介质存储有计算机程序，所述计算机程序被处理器执行时实现如权利要求1至4任一项所述的地震数据处理方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911061919.0A CN110618459B (zh) | 2019-11-01 | 2019-11-01 | 地震数据处理方法和装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911061919.0A CN110618459B (zh) | 2019-11-01 | 2019-11-01 | 地震数据处理方法和装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110618459A CN110618459A (zh) | 2019-12-27 |
| CN110618459B true CN110618459B (zh) | 2020-11-10 |
Family
ID=68927402
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201911061919.0A Active CN110618459B (zh) | 2019-11-01 | 2019-11-01 | 地震数据处理方法和装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110618459B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115469362B (zh) * | 2022-09-15 | 2023-10-10 | 中山大学 | 一种地震勘探中的能流密度矢量计算方法 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107894613A (zh) * | 2017-10-26 | 2018-04-10 | 中国石油天然气集团公司 | 弹性波矢量成像方法、装置、存储介质及设备 |
| CN108181653A (zh) * | 2018-01-16 | 2018-06-19 | 东北石油大学 | 针对vti介质逆时偏移方法、设备及介质 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US10295685B2 (en) * | 2017-04-06 | 2019-05-21 | Saudi Arabian Oil Company | Generating common image gather using wave-field separation |
| US20190018156A1 (en) * | 2017-07-11 | 2019-01-17 | The United State of America, as represented by the Secretary of the Department of the Interior | Highly accurate focal mechanism for microseismic envents |
-
2019
- 2019-11-01 CN CN201911061919.0A patent/CN110618459B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107894613A (zh) * | 2017-10-26 | 2018-04-10 | 中国石油天然气集团公司 | 弹性波矢量成像方法、装置、存储介质及设备 |
| CN108181653A (zh) * | 2018-01-16 | 2018-06-19 | 东北石油大学 | 针对vti介质逆时偏移方法、设备及介质 |
Non-Patent Citations (2)
| Title |
|---|
| TTI介质矢量波成像及波型分解;王娟 等;《地球物理学进展》;20141231;第29卷(第6期);第2754-2760页 * |
| 高效逆时偏移角度道集生成方法研究;吴成梁 等;《石油物探》;20180531;第57卷(第3期);第404-418页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110618459A (zh) | 2019-12-27 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Popovici et al. | 3-D imaging using higher order fast marching traveltimes | |
| EP3063562B1 (en) | Methods of subsurface exploration, computer program product and computer-readable storage medium | |
| CN107894613B (zh) | 弹性波矢量成像方法、装置、存储介质及设备 | |
| US20190072686A1 (en) | De-aliased source separation method | |
| WO2009066057A2 (en) | Seismic data processing method for surface related multiple attenuation | |
| US10215869B2 (en) | System and method of estimating anisotropy properties of geological formations using a self-adjoint pseudoacoustic wave propagator | |
| EP3062128B1 (en) | Method and apparatus for deblending seismic data using a non-blended dataset | |
| CN111025386B (zh) | 一种无分离假象的纵横波分离方法 | |
| US5530679A (en) | Method for migrating seismic data | |
| KR20180067650A (ko) | 진폭 보존을 갖는 fwi 모델 도메인 각도 스택들 | |
| EA032186B1 (ru) | Сейсмическая адаптивная фокусировка | |
| CN110895348B (zh) | 一种地震弹性阻抗低频信息提取方法、系统及存储介质 | |
| US20180059276A1 (en) | System and method for focusing seismic images | |
| WO2018102813A2 (en) | Seismic acquisition geometry full-waveform inversion | |
| WO2021127382A1 (en) | Full waveform inversion in the midpoint-offset domain | |
| CN111781635B (zh) | 海底四分量弹性波高斯束深度偏移方法和装置 | |
| EP3153890B1 (en) | Device and method for constrained seismic wave-field separation | |
| CN110618459B (zh) | 地震数据处理方法和装置 | |
| Roy Chowdhury | Seismic data acquisition and processing | |
| CN106597535A (zh) | 一种提高弹性波逆时偏移计算率和空间分辨率的方法 | |
| Sun et al. | Multiple attenuation using λ-f domain high-order and high-resolution Radon transform based on SL0 norm | |
| GB2425355A (en) | Residual move out processing of seismic data | |
| Amundsen et al. | Broadband seismic over/under sources and their designature-deghosting | |
| Xie et al. | Interpolation and regularization with the 3D CRS operator | |
| CN110703332A (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 |