CN107390267A - 一种同步挤压变换域的地震资料衰减补偿方法 - Google Patents
一种同步挤压变换域的地震资料衰减补偿方法 Download PDFInfo
- Publication number
- CN107390267A CN107390267A CN201710626472.1A CN201710626472A CN107390267A CN 107390267 A CN107390267 A CN 107390267A CN 201710626472 A CN201710626472 A CN 201710626472A CN 107390267 A CN107390267 A CN 107390267A
- Authority
- CN
- China
- Prior art keywords
- mrow
- msub
- msup
- frequency
- time
- 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
- 238000000034 method Methods 0.000 title claims abstract description 60
- 230000001360 synchronised effect Effects 0.000 title claims abstract description 48
- 230000008569 process Effects 0.000 claims abstract description 23
- 239000011159 matrix material Substances 0.000 claims description 29
- 230000009466 transformation Effects 0.000 claims description 25
- 238000006243 chemical reaction Methods 0.000 claims description 24
- 238000001228 spectrum Methods 0.000 claims description 10
- 238000013507 mapping Methods 0.000 claims description 6
- 230000002123 temporal effect Effects 0.000 claims description 6
- 230000008859 change Effects 0.000 claims description 4
- 230000009471 action Effects 0.000 claims description 3
- 238000004364 calculation method Methods 0.000 claims 1
- 230000000694 effects Effects 0.000 abstract description 8
- 238000009792 diffusion process Methods 0.000 abstract description 4
- 230000004807 localization Effects 0.000 abstract description 4
- 238000013341 scale-up Methods 0.000 abstract description 2
- 230000002463 transducing effect Effects 0.000 abstract 1
- 239000006185 dispersion Substances 0.000 description 6
- 238000001914 filtration Methods 0.000 description 6
- 230000015572 biosynthetic process Effects 0.000 description 4
- 238000013213 extrapolation Methods 0.000 description 4
- 238000003786 synthesis reaction Methods 0.000 description 3
- 230000006854 communication Effects 0.000 description 2
- 238000009826 distribution Methods 0.000 description 2
- 235000013399 edible fruits Nutrition 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 230000006641 stabilisation Effects 0.000 description 2
- 230000005428 wave function Effects 0.000 description 2
- 230000003044 adaptive effect Effects 0.000 description 1
- 230000002238 attenuated effect Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 239000012141 concentrate Substances 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000010363 phase shift Effects 0.000 description 1
- 230000021715 photosynthesis, light harvesting Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
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
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/40—Transforming data representation
- G01V2210/48—Other transforms
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
一种同步挤压变换域的地震资料衰减补偿方法，该方法在同步挤压变换域中，将衰减补偿过程表示为一个带有L1范数正则化的反问题，并通过迭代重加权算法进行求解，从而实现稳定化的衰减补偿。其中，同步挤压变换能压制时频能量的扩散，得到更为稀疏和高度局域化的时频表示，从而有利于对非平稳的地震信号实现衰减的补偿。该方法利用了信号在同步挤压变换域的稀疏性质，引入了L1范数正则化方法来稳定化衰减补偿过程，从而避免噪声放大问题。利用提出的方法实现地震衰减效应的补偿，可以增强地震资料的分辨率，为后续的地震反演和储层描述带来便利。
Description
技术领域
本发明属于地震勘探技术领域，涉及一种地震资料的时频域衰减补偿方法，尤其是一种同步挤压变换域的地震资料衰减补偿方法。
背景技术
地震波在地下介质的传播过程中，由于介质的衰减和频散效应，会降低地震波的幅度并使得子波相位畸变，从而显著的降低地震数据的分辨率。衰减和频散一般用与频率无关的品质因子Q来定量的表示，在估计得到此参数的情况下，可以对衰减和频散效应进行补偿，从而提高地震资料的分辨率以有利于后续的解释。
目前，衰减和频散的补偿可以描述为反褶积问题或者波场外推问题。反褶积是基于地震数据提高资料纵向分辨率的主要方法之一，传统的褶积模型要求地震是平稳的，也即地震子波在传播过程中保持不变。然而由于Q值效应的影响，地震子波在地下传播过程中是非平稳的。Margrave等将传统褶积模型推广为非平稳褶积模型，提出了Gabor域反褶积方法校正衰减效应。Wang等则利用自适应的窗口划分，提出了改进的Gabor域非平稳反褶积方案。此类方法中并不需要Q值的先验信息，但需假设子波满足最小相位。
反Q滤波是一种基于波场外推的地震数据衰减补偿方法，其需要预先估计得到Q值模型。反Q滤波算子包含相位(频散)补偿算子和幅度(衰减)补偿算子。相位补偿算子是一个无条件稳定算子，容易处理。但对于幅度补偿算子，其作用是放大被衰减信号的能量，从而在含噪声情况下会导致不稳定问题。因此，在反Q滤波算法中，必须对衰减补偿算子作稳定化处理以避免不稳定。
Wang基于Kolsky的衰减频散模型，利用波场外推相移算子，给出了一种Gabor域的稳定化反Q滤波算法。其中，幅度补偿算子的稳定化，从反演的角度来看，等价于求解一个带有L2范数正则化的反问题。Braga和Moraes则认为，连续小波变换相对Gabor变换，具有更好的时频局域化性，而且具有多分辨率特性，从而推广得到小波域的反Q滤波方法。
在时频域实现反Q滤波算法时，需要对地震道数据进行时频分解，如上述的Gabor变换和小波变换。利用Gabor变换或小波变换进行时频分析，是通过将待分析的信号与给定的一族模板函数作内积实现，Gabor变换的模板函数是时移和调制的窗函数，而小波变换的模板函数是时移和伸缩的小波基。一般地，模板函数的选取不可避免的影响了时频表示，同时不确定性原理也限制了表示的时频分辨率。在得到的时频表示中，存在明显的时频能量的扩散，将不利于对每个频率成分进行准确的衰减补偿。
发明内容
本发明目的在于克服现有技术的不足，提供一种同步挤压变换域的地震资料衰减补偿方法，该方法在同步挤压变换域中，利用了信号在同步挤压变换域的稀疏性质，将衰减补偿过程表示为一个带有L1范数正则化的反问题，并通过迭代重加权算法进行求解，从而实现稳定化的衰减补偿。
为达到上述目的，本发明采用以下的技术方案：
一种同步挤压变换域的地震资料衰减补偿方法，该方法在同步挤压变换域中，利用了信号在同步挤压变换域的稀疏性质，将衰减补偿过程表示为一个带有L1范数正则化的反问题，并通过迭代重加权算法进行求解，从而实现稳定化的衰减补偿。
本发明进一步的改进在于，具体包括以下步骤：
(1)对待补偿的时域地震道进行同步挤压变换；
(2)将同步挤压变换域补偿过程表示为带L1范数正则化的反问题，利用迭代重加权算法求解，得到矩阵形式下的时频表示S(ω,τ)，其中ω是频率，τ是时间；
(3)利用同步挤压变换的重构公式，重建补偿后的时域信号波形。
本发明进一步的改进在于，步骤(1)的同步挤压变换具体过程包括如下步骤：
步骤a：对待补偿的地震信号s0(t)作连续小波变换，记为W(a,τ)
其中，t是时间变量，a和τ分别是小波变换尺度伸缩和时间平移因子，ψ(t)是母小波函数，*表示复共轭运算；对于母小波函数ψ(t)，其傅里叶变换为满足如下的解析小波性质
根据帕斯瓦尔定理，小波变换的频率形式为
其中，是信号s0(t)的傅里叶变换，ξ是频率变量；
步骤b：计算瞬时频率
对连续小波变换结果，计算瞬时频率
其中W(a,τ)≠0；是虚数单位；
步骤c：时间-尺度平面映射到时间-频率平面
利用计算得到的瞬时频率，将信号的时频表示从时间尺度平面映射到同步挤压变换域的时间频率平面，也即考虑到频率变量ω和尺度变量a在实现时都是离散的，也就是小波变换需要在离散的尺度点ak进行计算，其中ak-ak-1＝(Δa)k；对同步挤压变换，则是在以ωl为中心的小箱中操作，其中ωl-ωl-1＝Δω；将小箱中的能量进行求和，则在离散的频率点ωl处，信号s0(t)的同步挤压变换S0(ω,τ)表示为
本发明进一步的改进在于，步骤(2)的具体过程包括如下步骤：
步骤a：时频表示的向量化，正向衰减算子G(ω,τ)表示为对角矩阵，得到反问题形式；
设待补偿信号s0(t)的同步挤压变换域时频表示S0(ω,τ)有L个频率点，M个时间点，则S0(ω,τ)看成是一个L行M列的矩阵S0，同样地，补偿后信号的时频表示S(ω,τ)看成是矩阵S，正向衰减算子G(ω,τ)看成是矩阵G；注意到信号衰减的过程其关于时间频率点是一对一点乘的关系；为了将其表示成矩阵向量形式，按照列优先的方式重排列，被衰减的信号的时频表示S0(ω,τ)重排为列向量s0＝vec(S0)，补偿后信号的时频表示S(ω,τ)重排为列向量s＝vec(S)；对于正向衰减算子G(ω,τ)，首先将其列向量化为g＝vec(G)，然后向量g的元素作为对角线元素，得到对角矩阵F＝diag(g)；从而，用矩阵乘法来表示时频域的衰减过程
s0＝Fs+n
其中，n是噪声的时频表示对应的列向量；从而，衰减补偿过程表示为如下的具有L1范数正则化的反问题，目标函数Ψ(s,γ)为：
其中，||·||2表示L2范数，是L1范数，sj是向量s的第j个元素，γ是正则化因子；L1范数正则化表示对反问题施加了稀疏性约束，这符合同步挤压变换域时频表示的稀疏特性；
步骤b：利用迭代重加权算法求解反问题
利用近似反问题改写为
式中，ε是一个给定的非负的小量；对此目标函数关于sj的导数为零得到
FT(Fs-s0)+γWs＝0
其中，W是一个对角矩阵
设已有解sk，对下一次迭代的解sk+1，其满足
FT(Fsk+1-s0)+γWk+1sk+1＝0
假设Wk+1＝Wk，从而
FT(Fsk+1-s0)+γWksk+1＝0
则得到迭代重加权格式为
sk+1＝(FTF+γWk)-1FTs0
迭代的停机准则设置为：
其中，δ是停机参数；或者当达到预先设定的最大迭代次数时，停止迭代；
步骤c：将反演得到的列向量s重排为矩阵形式下的时频表示S(ω,τ)。
本发明进一步的改进在于，步骤(3)中利用下面公式重建补偿后的时域信号波形；
式中是一个与母小波相关的常数，Re[·]是取实部运算，S(ωl,τ)是S(ω,τ)在离散频率点ωl处的取值。
与现有技术相比，本发明具有以下有益效果：
本发明在同步挤压变换域中，将衰减补偿过程表示为一个带有L1范数正则化的反问题，并通过迭代重加权算法进行求解，从而实现稳定化的衰减补偿。本发明中利用同步挤压变换压制时频能量的扩散，得到更为稀疏和高度局域化的时频表示，从而有利于对非平稳的地震信号实现衰减的补偿。本发明利用了信号在同步挤压变换域的稀疏性质，引入了L1范数正则化方法来稳定化衰减补偿过程，从而避免噪声放大问题。本发明提出的方法实现地震衰减效应的补偿，可以增强地震资料的分辨率，为后续的地震反演和储层描述带来便利。
附图说明
图1(a)为单频余弦信号的波形图；
图1(b)为单频余弦信号的小波变换时频图；
图1(c)为单频余弦信号的同步挤压变换时频图；
图2(a)为实际地震信号图；
图2(b)连续小波变换时频图；
图2(c)为同步挤压变换时频图；
图3(a)为合成的无衰减地震信号图；
图3(b)为无衰减信号的同步挤压变换时频图；
图3(c)为合成的含衰减地震信号图；
图3(d)为含衰减信号的同步挤压变换时频图；
图4(a)为本发明方法实现衰减补偿后的补偿后的地震信号图；
图4(b)为本发明方法实现衰减补偿后的同步挤压变换时频图；
图5(a)为含噪情况下无正则化直接补偿后的地震信号图；
图5(b)为本发明方法补偿后的地震信号图；
图6(a)为衰减补偿前的实际地震记录剖面图；
图6(b)为衰减补偿后的实际地震记录剖面图。
具体实施方式
下面结合附图对本发明做进一步详细描述：
本发明在同步挤压变换域中，利用了信号在同步挤压变换域的稀疏性质，将衰减补偿过程表示为一个带有L1范数正则化的反问题，并通过迭代重加权算法进行求解，从而实现稳定化的衰减补偿，具体如下：
根据波场外推理论，对含衰减的地震道s0(t)，其时频域表示为S0(ω,τ)，在时频域的衰减补偿过程可以描述为：
其中，S(ω,τ)是衰减补偿后信号的时频表示，ω表示频率，τ表示时间，Λ(ω,τ)是衰减补偿算子，Q是预先给定的介质品质因子。
衰减补偿算子Λ(ω,τ)是一个关于时间和频率的增函数，如果直接进行补偿，如果信号中存在噪声，此补偿算子会在放大有用信号的同时也放大信号噪声，从而使得衰减补偿过程不稳定。将上式可以改写为信号衰减的过程：
其中，G(ω,τ)称为正向衰减算子，此算子是关于时间和频率的减函数，不会导致不稳定问题。因此，本发明在同步挤压变换域，从信号在时频域的衰减过程出发，将待补偿的地震道数据的时频表示S0(ω,τ)作为观测数据，将补偿后信号的时频表示S(ω,τ)作为模型参数，通过建立反问题来求解S(ω,τ)，从而在时频域实现稳定化的衰减补偿。
本发明在同步挤压变换域实现地震资料衰减补偿，包括以下步骤：
(1)对待补偿的时域地震道进行同步挤压变换，其中同步挤压变换具体包括如下步骤：
步骤a：对待补偿的地震信号s0(t)作连续小波变换，记为W(a,τ)
其中，t是时间变量，a和τ分别是小波变换尺度伸缩和时间平移因子，ψ(t)是母小波函数，*表示复共轭运算。对于母小波函数ψ(t)，其傅里叶变换为满足如下的解析小波性质
根据帕斯瓦尔定理，小波变换的频率形式为
其中，是信号s0(t)的傅里叶变换，ξ是频率变量。
考虑一个纯谐波信号s0(t)＝Acos(ωt)，A是幅度，其小波变换为
由此式可知，如果小波集中在ξ＝ω0处，那么小波变换的结果聚焦在a＝ω0/ω。然而，实际的小波函数在频率域总是有一定的支撑范围，也就意味着纯谐波信号的小波变换在时间尺度平面上是在尺度a＝ω0/ω附件能量扩散的。为了减少这种能量扩散，得到更为准确更为局域化的时频分布，需将扩散的能量挤压回到对应的瞬时频率处。
步骤b：计算瞬时频率
对连续小波变换结果，计算瞬时频率
其中W(a,τ)≠0；是虚数单位。
步骤c：时间-尺度平面映射到时间-频率平面
利用计算得到的瞬时频率，可以将信号的时频表示从时间尺度平面映射到同步挤压变换域的时间频率平面，也即考虑到频率变量ω和尺度变量a在实现时都是离散的，也就是小波变换需要在离散的尺度点ak进行计算，其中ak-ak-1＝(Δa)k；对同步挤压变换，则是在以ωl为中心的小箱中操作，其中ωl-ωl-1＝Δω。将小箱中的能量进行求和，则在离散的频率点ωl处，信号s0(t)的同步挤压变换S0(ω,τ)表示为
如图1(a)、图1(b)、图1(c)所示，一个单频余弦信号s(t)＝cos(40πt)，其瞬时频率为20Hz，对此信号进行小波变换，并将时间尺度图转换为时间-频率图，可以明显地观察到信号的时频能量沿着频率轴扩散，而同步挤压变换的时频表示聚焦在20Hz。图2还给出了一道实际地震信号的波形，及其小波变换与同步挤压变换的时频图。从图2(a)、2(b)2(c)可以看到，同步挤压变换具有更好的时频聚集性和稀疏性，能够更准确的刻画信号的时频分布。
(2)将同步挤压变换域补偿过程表示为带L1范数正则化的反问题，利用迭代重加权算法求解，具体包括如下步骤：
步骤a：时频表示的向量化，正向衰减算子G(ω,τ)表示为对角矩阵，得到反问题形式；
设待补偿信号s0(t)的同步挤压变换域时频表示S0(ω,τ)有L个频率点，M个时间点，则S0(ω,τ)可以看成是一个L行M列的矩阵S0，同样地，补偿后信号的时频表示S(ω,τ)可以看成是矩阵S，正向衰减算子G(ω,τ)可以看成是矩阵G。注意到信号衰减的过程其关于时间频率点是一对一点乘的关系。为了将其表示成矩阵向量形式，按照列优先的方式重排列，被衰减的信号的时频表示S0(ω,τ)重排为列向量s0＝vec(S0)，补偿后信号的时频表示S(ω,τ)重排为列向量s＝vec(S)。对于正向衰减算子G(ω,τ)，首先将其列向量化为g＝vec(G)，然后向量g的元素作为对角线元素，得到对角矩阵F＝diag(g)。从而，可以用矩阵乘法来表示时频域的衰减过程
s0＝Fs+n
其中，n可以看成是噪声的时频表示对应的列向量。从而，衰减补偿过程表示为如下的具有L1范数正则化的反问题，目标函数Ψ(s,γ)为：
其中，||·||2表示L2范数，是L1范数，sj是向量s的第j个元素，γ是正则化因子。L1范数正则化表示对反问题施加了稀疏性约束，这符合同步挤压变换域时频表示的稀疏特性。
步骤b：利用迭代重加权算法求解反问题
利用近似反问题可以改写为
式中，ε是一个给定的非负的小量。对此目标函数关于sj的导数为零得到
FT(Fs-s0)+γWs＝0
其中，W是一个对角矩阵
设已有解sk，对下一次迭代的解sk+1，其满足
FT(Fsk+1-s0)+γWk+1sk+1＝0
假设Wk+1＝Wk，从而
FT(Fsk+1-s0)+γWksk+1＝0
则得到迭代重加权格式为
sk+1＝(FTF+γWk)-1FTs0
迭代的停机准则设置为：
其中，δ是停机参数。或者当达到预先设定的最大迭代次数时，也停止迭代。
步骤c：将反演得到的列向量s重排为矩阵形式下的时频表示S(ω,τ)
(3)利用同步挤压变换的重构公式，重建补偿后的时域信号波形
式中是一个与母小波相关的常数，Re[·]是取实部运算，S(ωl,τ)是S(ω,τ)在离散频率点ωl处的取值。
数值仿真结果如下：
首先将提出的方法应用于合成记录的衰减补偿。
震源选取峰值频率为35Hz的Ricker子波，品质因子设置为Q＝60，时间采样间隔为1毫秒。合成的含衰减的非平稳地震记录如图3(a)、3(b)、3(c)、3(d)所示，并利用同步挤压变换对此信号作时频分析。为了比较，图3(a)、3(b)、3(c)、3(d)也给出了无衰减情况下的平稳地震道及其相应的时频图。对比可知，由于衰减效应的作用，在时间波形上，地震信号的幅度整体是下降趋势，波形展宽，且深层数据受到的影响大。在时频图中，信号的高频能量损失明显，尤其在深层。
对于衰减的地震道，利用本发明提出的方法，选取正则化参数γ＝10-7，ε＝10-10，停机参数δ＝10-3，最大迭代次数为50次。通过迭代反演得到补偿后信号的时频表示，然后利用同步挤压变换的逆变换重建得到补偿后的时间信号，结果如图4(a)和4(b)所示。经过补偿，恢复了信号的幅度和时频域能量，提高了信号的分辨率。
对衰减后的信号加入15dB的高斯噪声，如果不引入L1范数正则化方法对衰减补偿过程进行稳定化，也就是γ＝0，会导致出现不稳定性问题，如图5(a)所示。利用本发明提出的稳定化补偿方法，设置γ＝5×10-7，ε＝10-10，停机参数δ＝10-3，可以得到稳定的补偿结果，如图5(b)所示。
本发明的方法用于实际地震数据的衰减补偿。
如图6(a)所示，此叠后地震数据中有100道，时间从1.3秒到2.0秒，时间采样间隔2毫秒。在提出的方法中正则化参数γ＝10-7，ε＝10-10，预先估计得到的介质品质因子问为Q＝60。对剖面中的每一道数据分别进行补偿，最终获得的补偿后剖面如图6(b)所示。结果显示，经过本发明提出的方法实现衰减补偿，地震资料的分辨率得到了增强，比如图中箭头和椭圆框标记的位置，这将有利于解释更精细的地层结构。
Claims (5)
1.一种同步挤压变换域的地震资料衰减补偿方法，其特征在于，该方法在同步挤压变换域中，利用了信号在同步挤压变换域的稀疏性质，将衰减补偿过程表示为一个带有L1范数正则化的反问题，并通过迭代重加权算法进行求解，从而实现稳定化的衰减补偿。
2.根据权利要求1所述的一种同步挤压变换域的地震资料衰减补偿方法，其特征在于，具体包括以下步骤：
(1)对待补偿的时域地震道进行同步挤压变换；
(2)将同步挤压变换域补偿过程表示为带L1范数正则化的反问题，利用迭代重加权算法求解，得到矩阵形式下的时频表示S(ω,τ)，其中ω是频率，τ是时间；
(3)利用同步挤压变换的重构公式，重建补偿后的时域信号波形。
3.根据权利要求2所述的一种同步挤压变换域的地震资料衰减补偿方法，其特征在于，步骤(1)的同步挤压变换具体过程包括如下步骤：
步骤a：对待补偿的地震信号s0(t)作连续小波变换，记为W(a,τ)
<mrow>
<mi>W</mi>
<mrow>
<mo>(</mo>
<mi>a</mi>
<mo>,</mo>
<mi>&tau;</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mo>&Integral;</mo>
<msub>
<mi>s</mi>
<mn>0</mn>
</msub>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<msup>
<mi>a</mi>
<mrow>
<mo>-</mo>
<mn>1</mn>
<mo>/</mo>
<mn>2</mn>
</mrow>
</msup>
<msup>
<mi>&psi;</mi>
<mo>*</mo>
</msup>
<mrow>
<mo>(</mo>
<mfrac>
<mrow>
<mi>t</mi>
<mo>-</mo>
<mi>&tau;</mi>
</mrow>
<mi>a</mi>
</mfrac>
<mo>)</mo>
</mrow>
<mi>d</mi>
<mi>t</mi>
</mrow>
其中，t是时间变量，a和τ分别是小波变换尺度伸缩和时间平移因子，ψ(t)是母小波函数，*表示复共轭运算；对于母小波函数ψ(t)，其傅里叶变换为满足如下的解析小波性质
<mrow>
<mover>
<mi>&psi;</mi>
<mo>^</mo>
</mover>
<mrow>
<mo>(</mo>
<mi>&omega;</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mn>0</mn>
<mo>,</mo>
<mi>&omega;</mi>
<mo><</mo>
<mn>0</mn>
</mrow>
根据帕斯瓦尔定理，小波变换的频率形式为
<mrow>
<mi>W</mi>
<mrow>
<mo>(</mo>
<mi>a</mi>
<mo>,</mo>
<mi>&tau;</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mrow>
<mn>2</mn>
<mi>&pi;</mi>
</mrow>
</mfrac>
<mo>&Integral;</mo>
<msub>
<mover>
<mi>s</mi>
<mo>^</mo>
</mover>
<mn>0</mn>
</msub>
<mrow>
<mo>(</mo>
<mi>&xi;</mi>
<mo>)</mo>
</mrow>
<msup>
<mi>a</mi>
<mrow>
<mn>1</mn>
<mo>/</mo>
<mn>2</mn>
</mrow>
</msup>
<msup>
<mover>
<mi>&psi;</mi>
<mo>^</mo>
</mover>
<mo>*</mo>
</msup>
<mrow>
<mo>(</mo>
<mi>a</mi>
<mi>&xi;</mi>
<mo>)</mo>
</mrow>
<msup>
<mi>e</mi>
<mrow>
<mi>i</mi>
<mi>&tau;</mi>
<mi>&xi;</mi>
</mrow>
</msup>
<mi>d</mi>
<mi>&xi;</mi>
</mrow>
其中，是信号s0(t)的傅里叶变换，ξ是频率变量；
步骤b：计算瞬时频率
对连续小波变换结果，计算瞬时频率
<mrow>
<msub>
<mi>&omega;</mi>
<msub>
<mi>s</mi>
<mn>0</mn>
</msub>
</msub>
<mrow>
<mo>(</mo>
<mi>a</mi>
<mo>,</mo>
<mi>&tau;</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mo>-</mo>
<mi>i</mi>
<msup>
<mrow>
<mo>(</mo>
<mi>W</mi>
<mo>(</mo>
<mrow>
<mi>a</mi>
<mo>,</mo>
<mi>&tau;</mi>
</mrow>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mfrac>
<mrow>
<mo>&part;</mo>
<mi>W</mi>
<mrow>
<mo>(</mo>
<mi>a</mi>
<mo>,</mo>
<mi>&tau;</mi>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>&tau;</mi>
</mrow>
</mfrac>
</mrow>
其中W(a,τ)≠0；是虚数单位；
步骤c：时间-尺度平面映射到时间-频率平面
利用计算得到的瞬时频率，将信号的时频表示从时间尺度平面映射到同步挤压变换域的时间频率平面，也即考虑到频率变量ω和尺度变量a在实现时都是离散的，也就是小波变换需要在离散的尺度点ak进行计算，其中ak-ak-1＝(Δa)k；对同步挤压变换，则是在以ωl为中心的小箱中操作，其中ωl-ωl-1＝Δω；将小箱中的能量进行求和，则在离散的频率点ωl处，信号s0(t)的同步挤压变换S0(ω,τ)表示为
<mrow>
<msub>
<mi>S</mi>
<mn>0</mn>
</msub>
<mrow>
<mo>(</mo>
<msub>
<mi>&omega;</mi>
<mi>l</mi>
</msub>
<mo>,</mo>
<mi>&tau;</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msup>
<mrow>
<mo>(</mo>
<mi>&Delta;</mi>
<mi>&omega;</mi>
<mo>)</mo>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<munder>
<mo>&Sigma;</mo>
<mrow>
<msub>
<mi>a</mi>
<mi>k</mi>
</msub>
<mo>:</mo>
<mo>|</mo>
<mi>&omega;</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>a</mi>
<mi>k</mi>
</msub>
<mo>,</mo>
<mi>b</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<msub>
<mi>&omega;</mi>
<mi>l</mi>
</msub>
<mo>|</mo>
<mo>&le;</mo>
<mi>&Delta;</mi>
<mi>&omega;</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
</munder>
<mi>W</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>a</mi>
<mi>k</mi>
</msub>
<mo>,</mo>
<mi>&tau;</mi>
<mo>)</mo>
</mrow>
<msubsup>
<mi>a</mi>
<mi>k</mi>
<mrow>
<mo>-</mo>
<mn>3</mn>
<mo>/</mo>
<mn>2</mn>
</mrow>
</msubsup>
<msub>
<mrow>
<mo>(</mo>
<mi>&Delta;</mi>
<mi>a</mi>
<mo>)</mo>
</mrow>
<mi>k</mi>
</msub>
<mo>.</mo>
</mrow>
4.根据权利要求3所述的一种同步挤压变换域的地震资料衰减补偿方法，其特征在于，步骤(2)的具体过程包括如下步骤：
步骤a：时频表示的向量化，正向衰减算子G(ω,τ)表示为对角矩阵，得到反问题形式；
设待补偿信号s0(t)的同步挤压变换域时频表示S0(ω,τ)有L个频率点，M个时间点，则S0(ω,τ)看成是一个L行M列的矩阵S0，同样地，补偿后信号的时频表示S(ω,τ)看成是矩阵S，正向衰减算子G(ω,τ)看成是矩阵G；注意到信号衰减的过程其关于时间频率点是一对一点乘的关系；为了将其表示成矩阵向量形式，按照列优先的方式重排列，被衰减的信号的时频表示S0(ω,τ)重排为列向量s0＝vec(S0)，补偿后信号的时频表示S(ω,τ)重排为列向量s＝vec(S)；对于正向衰减算子G(ω,τ)，首先将其列向量化为g＝vec(G)，然后向量g的元素作为对角线元素，得到对角矩阵F＝diag(g)；从而，用矩阵乘法来表示时频域的衰减过程
s0＝Fs+n
其中，n是噪声的时频表示对应的列向量；从而，衰减补偿过程表示为如下的具有L1范数正则化的反问题，目标函数Ψ(s,γ)为：
<mrow>
<mi>&Psi;</mi>
<mrow>
<mo>(</mo>
<mi>s</mi>
<mo>,</mo>
<mi>&gamma;</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mo>|</mo>
<mo>|</mo>
<msub>
<mi>s</mi>
<mn>0</mn>
</msub>
<mo>-</mo>
<mi>F</mi>
<mi>s</mi>
<mo>|</mo>
<msubsup>
<mo>|</mo>
<mn>2</mn>
<mn>2</mn>
</msubsup>
<mo>+</mo>
<mi>&gamma;</mi>
<mo>|</mo>
<mo>|</mo>
<mi>s</mi>
<mo>|</mo>
<msub>
<mo>|</mo>
<mn>1</mn>
</msub>
</mrow>
其中，||·||2表示L2范数，是L1范数，sj是向量s的第j个元素，γ是正则化因子；L1范数正则化表示对反问题施加了稀疏性约束，这符合同步挤压变换域时频表示的稀疏特性；
步骤b：利用迭代重加权算法求解反问题
利用近似反问题改写为
<mrow>
<mi>&Psi;</mi>
<mrow>
<mo>(</mo>
<mi>s</mi>
<mo>,</mo>
<mi>&gamma;</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mo>|</mo>
<mo>|</mo>
<msub>
<mi>s</mi>
<mn>0</mn>
</msub>
<mo>-</mo>
<mi>F</mi>
<mi>s</mi>
<mo>|</mo>
<msubsup>
<mo>|</mo>
<mn>2</mn>
<mn>2</mn>
</msubsup>
<mo>+</mo>
<mi>&gamma;</mi>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mrow>
<mi>L</mi>
<mo>&times;</mo>
<mi>M</mi>
</mrow>
</munderover>
<msqrt>
<mrow>
<msubsup>
<mi>s</mi>
<mi>j</mi>
<mn>2</mn>
</msubsup>
<mo>+</mo>
<msup>
<mi>&epsiv;</mi>
<mn>2</mn>
</msup>
</mrow>
</msqrt>
</mrow>
式中，ε是一个给定的非负的小量；对此目标函数关于sj的导数为零得到
FT(Fs-s0)+γWs＝0
其中，W是一个对角矩阵
<mrow>
<mi>W</mi>
<mo>=</mo>
<mi>d</mi>
<mi>i</mi>
<mi>a</mi>
<mi>g</mi>
<mrow>
<mo>(</mo>
<mrow>
<mn>1</mn>
<mo>/</mo>
<mrow>
<mo>(</mo>
<mrow>
<mn>2</mn>
<msqrt>
<mrow>
<msubsup>
<mi>s</mi>
<mn>1</mn>
<mn>2</mn>
</msubsup>
<mo>+</mo>
<msup>
<mi>&epsiv;</mi>
<mn>2</mn>
</msup>
</mrow>
</msqrt>
</mrow>
<mo>)</mo>
</mrow>
<mo>,</mo>
<mn>...</mn>
<mo>,</mo>
<mn>1</mn>
<mo>/</mo>
<mrow>
<mo>(</mo>
<mrow>
<mn>2</mn>
<msqrt>
<mrow>
<msubsup>
<mi>s</mi>
<mi>j</mi>
<mn>2</mn>
</msubsup>
<mo>+</mo>
<msup>
<mi>&epsiv;</mi>
<mn>2</mn>
</msup>
</mrow>
</msqrt>
</mrow>
<mo>)</mo>
</mrow>
<mo>,</mo>
<mn>...</mn>
<mo>,</mo>
<mn>1</mn>
<mo>/</mo>
<mrow>
<mo>(</mo>
<mrow>
<mn>2</mn>
<msqrt>
<mrow>
<msubsup>
<mi>s</mi>
<mrow>
<mi>L</mi>
<mo>&times;</mo>
<mi>M</mi>
</mrow>
<mn>2</mn>
</msubsup>
<mo>+</mo>
<msup>
<mi>&epsiv;</mi>
<mn>2</mn>
</msup>
</mrow>
</msqrt>
</mrow>
<mo>)</mo>
</mrow>
</mrow>
<mo>)</mo>
</mrow>
</mrow>
设已有解sk，对下一次迭代的解sk+1，其满足
FT(Fsk+1-s0)+γWk+1sk+1＝0
假设Wk+1＝Wk，从而
FT(Fsk+1-s0)+γWksk+1＝0
则得到迭代重加权格式为
sk+1＝(FTF+γWk)-1FTs0
迭代的停机准则设置为：
<mrow>
<mfrac>
<mrow>
<mo>|</mo>
<mo>|</mo>
<msup>
<mi>s</mi>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>-</mo>
<msup>
<mi>s</mi>
<mi>k</mi>
</msup>
<mo>|</mo>
<msub>
<mo>|</mo>
<mn>2</mn>
</msub>
</mrow>
<mrow>
<mo>|</mo>
<mo>|</mo>
<msup>
<mi>s</mi>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>|</mo>
<msub>
<mo>|</mo>
<mn>2</mn>
</msub>
</mrow>
</mfrac>
<mo><</mo>
<mi>&delta;</mi>
</mrow>
其中，δ是停机参数；或者当达到预先设定的最大迭代次数时，停止迭代；
步骤c：将反演得到的列向量s重排为矩阵形式下的时频表示S(ω,τ)。
5.根据权利要求4所述的一种同步挤压变换域的地震资料衰减补偿方法，其特征在于，步骤(3)中利用下面公式重建补偿后的时域信号波形；
<mrow>
<mi>s</mi>
<mrow>
<mo>(</mo>
<mi>&tau;</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mi>Re</mi>
<mo>&lsqb;</mo>
<msubsup>
<mi>C</mi>
<mi>&psi;</mi>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msubsup>
<munder>
<mo>&Sigma;</mo>
<mi>l</mi>
</munder>
<mi>S</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>&omega;</mi>
<mi>l</mi>
</msub>
<mo>,</mo>
<mi>&tau;</mi>
<mo>)</mo>
</mrow>
<mrow>
<mo>(</mo>
<mi>&Delta;</mi>
<mi>&omega;</mi>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
</mrow>
式中是一个与母小波相关的常数，Re[·]是取实部运算，S(ωl,τ)是S(ω,τ)在离散频率点ωl处的取值。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710626472.1A CN107390267B (zh) | 2017-07-27 | 2017-07-27 | 一种同步挤压变换域的地震资料衰减补偿方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710626472.1A CN107390267B (zh) | 2017-07-27 | 2017-07-27 | 一种同步挤压变换域的地震资料衰减补偿方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107390267A true CN107390267A (zh) | 2017-11-24 |
| CN107390267B CN107390267B (zh) | 2019-03-01 |
Family
ID=60341953
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710626472.1A Active CN107390267B (zh) | 2017-07-27 | 2017-07-27 | 一种同步挤压变换域的地震资料衰减补偿方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107390267B (zh) |
Cited By (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108710851A (zh) * | 2018-05-21 | 2018-10-26 | 闽南师范大学 | 地震信号随机噪声衰减方法、终端设备及存储介质 |
| CN108845357A (zh) * | 2018-06-13 | 2018-11-20 | 成都信息工程大学 | 一种基于同步挤压小波变换估计地层等效品质因子的方法 |
| CN109143370A (zh) * | 2018-07-25 | 2019-01-04 | 中国地震局地球物理研究所 | 地震动加速度记录基线漂移的校正方法 |
| CN109521421A (zh) * | 2018-01-27 | 2019-03-26 | 河南工业大学 | 一种探地雷达薄层目标识别定位方法 |
| CN110579800A (zh) * | 2019-10-19 | 2019-12-17 | 西南石油大学 | 一种基于高精度同步挤压变换的地震数据数字处理方法 |
| CN110873900A (zh) * | 2018-09-04 | 2020-03-10 | 中国石油化工股份有限公司 | 频率域叠前地震道q补偿方法及系统 |
| CN111708081A (zh) * | 2020-05-29 | 2020-09-25 | 成都理工大学 | 考虑衰减频散的深度域地震记录合成方法 |
| CN112394402A (zh) * | 2019-08-19 | 2021-02-23 | 中国石油化工股份有限公司 | 基于同步挤压小波变换检测微地震信号的方法和系统 |
| CN112904416A (zh) * | 2021-01-20 | 2021-06-04 | 中国石油大学(北京) | 基于反射结构约束的多道吸收补偿处理方法 |
| CN113777650A (zh) * | 2021-08-12 | 2021-12-10 | 西安交通大学 | 一种基于混合范数和小波变换的稀疏时频谱分解方法、装置、设备及存储介质 |
| CN114842800A (zh) * | 2022-05-19 | 2022-08-02 | 姜英 | 采用离线标定减弱amoled显示屏退化的补偿方法 |
| CN114842800B (zh) * | 2022-05-19 | 2024-05-31 | 姜英 | 采用离线标定减弱amoled显示屏退化的补偿方法 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103645502A (zh) * | 2013-12-11 | 2014-03-19 | 中国海洋石油总公司 | 一种曲波域中地震波衰减补偿方法 |
| WO2014191011A1 (en) * | 2013-05-27 | 2014-12-04 | Statoil Petroleum As | High resolution estimation of attenuation from vertical seismic profiles |
| CN104880730A (zh) * | 2015-03-27 | 2015-09-02 | 西安交通大学 | 基于Synchrosqueezing变换的地震资料时频分析和衰减估计方法 |
| CN106291700A (zh) * | 2016-09-28 | 2017-01-04 | 西安交通大学 | 基于同步挤压变换的地震加权平均瞬时频率提取方法 |
| CN106443785A (zh) * | 2016-11-03 | 2017-02-22 | 中国矿业大学(北京) | 绕射波场获取方法和装置 |
-
2017
- 2017-07-27 CN CN201710626472.1A patent/CN107390267B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2014191011A1 (en) * | 2013-05-27 | 2014-12-04 | Statoil Petroleum As | High resolution estimation of attenuation from vertical seismic profiles |
| CN103645502A (zh) * | 2013-12-11 | 2014-03-19 | 中国海洋石油总公司 | 一种曲波域中地震波衰减补偿方法 |
| CN104880730A (zh) * | 2015-03-27 | 2015-09-02 | 西安交通大学 | 基于Synchrosqueezing变换的地震资料时频分析和衰减估计方法 |
| CN106291700A (zh) * | 2016-09-28 | 2017-01-04 | 西安交通大学 | 基于同步挤压变换的地震加权平均瞬时频率提取方法 |
| CN106443785A (zh) * | 2016-11-03 | 2017-02-22 | 中国矿业大学(北京) | 绕射波场获取方法和装置 |
Non-Patent Citations (1)
| Title |
|---|
| 鲁才: "《多维地震信号正则化处理方法研究》", 《电子科技大学硕士学位论文》 * |
Cited By (17)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109521421A (zh) * | 2018-01-27 | 2019-03-26 | 河南工业大学 | 一种探地雷达薄层目标识别定位方法 |
| CN108710851B (zh) * | 2018-05-21 | 2021-04-16 | 闽南师范大学 | 地震信号随机噪声衰减方法、终端设备及存储介质 |
| CN108710851A (zh) * | 2018-05-21 | 2018-10-26 | 闽南师范大学 | 地震信号随机噪声衰减方法、终端设备及存储介质 |
| CN108845357A (zh) * | 2018-06-13 | 2018-11-20 | 成都信息工程大学 | 一种基于同步挤压小波变换估计地层等效品质因子的方法 |
| CN108845357B (zh) * | 2018-06-13 | 2020-12-22 | 成都信息工程大学 | 一种基于同步挤压小波变换估计地层等效品质因子的方法 |
| CN109143370A (zh) * | 2018-07-25 | 2019-01-04 | 中国地震局地球物理研究所 | 地震动加速度记录基线漂移的校正方法 |
| CN110873900A (zh) * | 2018-09-04 | 2020-03-10 | 中国石油化工股份有限公司 | 频率域叠前地震道q补偿方法及系统 |
| CN110873900B (zh) * | 2018-09-04 | 2021-07-27 | 中国石油化工股份有限公司 | 频率域叠前地震道q补偿方法及系统 |
| CN112394402A (zh) * | 2019-08-19 | 2021-02-23 | 中国石油化工股份有限公司 | 基于同步挤压小波变换检测微地震信号的方法和系统 |
| CN110579800A (zh) * | 2019-10-19 | 2019-12-17 | 西南石油大学 | 一种基于高精度同步挤压变换的地震数据数字处理方法 |
| CN111708081A (zh) * | 2020-05-29 | 2020-09-25 | 成都理工大学 | 考虑衰减频散的深度域地震记录合成方法 |
| CN111708081B (zh) * | 2020-05-29 | 2022-04-15 | 成都理工大学 | 考虑衰减频散的深度域地震记录合成方法 |
| CN112904416A (zh) * | 2021-01-20 | 2021-06-04 | 中国石油大学(北京) | 基于反射结构约束的多道吸收补偿处理方法 |
| CN113777650A (zh) * | 2021-08-12 | 2021-12-10 | 西安交通大学 | 一种基于混合范数和小波变换的稀疏时频谱分解方法、装置、设备及存储介质 |
| CN113777650B (zh) * | 2021-08-12 | 2022-10-25 | 西安交通大学 | 一种基于混合范数和小波变换的稀疏时频谱分解方法、装置、设备及存储介质 |
| CN114842800A (zh) * | 2022-05-19 | 2022-08-02 | 姜英 | 采用离线标定减弱amoled显示屏退化的补偿方法 |
| CN114842800B (zh) * | 2022-05-19 | 2024-05-31 | 姜英 | 采用离线标定减弱amoled显示屏退化的补偿方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107390267B (zh) | 2019-03-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107390267A (zh) | 一种同步挤压变换域的地震资料衰减补偿方法 | |
| US11333783B2 (en) | Integrated method for estimation of seismic wavelets and synthesis of seismic records in depth domain | |
| CN104849756B (zh) | 一种提高地震数据分辨率增强有效弱信号能量的方法 | |
| CN106707334B (zh) | 一种提高地震资料分辨率的方法 | |
| CN106680874A (zh) | 基于波形形态特征稀疏化建模的谐波噪声压制方法 | |
| CN106646303A (zh) | 一种欠采样磁共振波谱的快速重建方法 | |
| CN107132579A (zh) | 一种保地层结构的地震波衰减补偿方法 | |
| CN102053273A (zh) | 一种对地震波信号进行反q滤波的方法 | |
| CN103728660A (zh) | 基于地震数据的多道匹配追踪方法 | |
| CN112327358B (zh) | 一种粘滞性介质中声波地震数据正演模拟方法 | |
| CN107894612A (zh) | 一种q吸收衰减补偿的声波阻抗反演方法及系统 | |
| Zoukaneri et al. | A combined Wigner-Ville and maximum entropy method for high-resolution time-frequency analysis of seismic data | |
| CN104932018A (zh) | 补偿变分辨率因子s变换的复时-频谱提高地震剖面分辨率的方法 | |
| CN103777238A (zh) | 一种纯纵波各向异性波场模拟方法 | |
| CN110967749A (zh) | 一种vsp地震资料频变q值估计与反q滤波方法 | |
| CN104422956B (zh) | 一种基于稀疏脉冲反演的高精度地震谱分解方法 | |
| Chen et al. | A statistical instantaneous frequency estimator for high-concentration time-frequency representation | |
| Wang et al. | Data-driven pre-stack AVO inversion method based on fast orthogonal dictionary | |
| CN104932008B (zh) | 补偿j变换的复时‑频谱提高地震剖面分辨率的方法 | |
| Cheng et al. | Structural geosteering constrained multi-trace sparse reflectivity inversion based on mixed norms | |
| CN104932009B (zh) | 补偿Morlet小波变换的复时‑频谱提高地震剖面分辨率的方法 | |
| Xin et al. | Seismic high-resolution processing method based on spectral simulation and total variation regularization constraints | |
| CN106443771B (zh) | 提高转换波地震数据分辨率的方法及其速度反演方法 | |
| CN112327354A (zh) | 一种提高低频弱信号的方法、装置、电子设备和可读介质 | |
| CN112925013A (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 |