CN112526599A - 基于加权l1范数稀疏准则的子波相位估计方法及系统 - Google Patents

Abstract

本发明提供一种基于加权L1范数稀疏准则的子波相位估计方法及系统，该方法包括：S1、获取地震数据，并给定初始地震子波；S2、构建不同相位的子波群；S3、利用构建的子波群进行基于稀疏贝叶斯学习的稀疏反演，得到反演结果；S4、采用加权L1范数稀疏准则量化评价所述反演结果的稀疏性，获取地震子波相位和对应的反射系数序列中的一项。该系统包括：初始子波给定单元、子波群构建单元、反演单元、评价单元。本发明建立了加权L1范数稀疏准则与地震子波相位之间的关系，通过加权L1范数和稀疏贝叶斯学习的双重稀疏表示，分离了子波相位和反射系数的干涉耦合效应。

Description

基于加权L1范数稀疏准则的子波相位估计方法及系统

技术领域

本发明涉及地震勘探反演技术领域，尤其涉及一种基于加权L1范数稀疏准则的子波相位估计方法及系统。

背景技术

对反射波地震数据进行反褶积处理可以得到具有超高分辨率的反射系数或波阻抗结果从而获得尽可能清晰的油气藏图像，而子波正确与否对反褶积结果起着决定性作用。然而经过一系列处理后的地震剖面通常会有残留的子波相位，而残留的子波相位会造成时域主峰的偏移和旁瓣数目的变化，并有可能彻底改变子波的相位特性及因果特性，导致后续的反褶积处理结果误差变大并且产生假象进而影响后续的地质解释和勘探开发任务。

对于子波相位校正的约束条件一般分为两类：井资料约束与准则函数约束。当需要处理的数据中有大量井资料的时候，那么利用测井数据对井旁地震记录进行相位估计，然后将估计的结果应用于整个剖面，对井旁道地震数据进行确定性子波估计的方法是最准确的。然而实际上由于测井数据较为稀少，确定性子波估计方法通常难以应用在实际场景当中。此外，大量的学者研究发现不同的测井数据估计的子波通常会有相位上的差异。因此，对于大部分地震数据都需要从数据本身出发进行估计子波即采用统计性子波估计方法。

最早的子波常相位估计校正技术是在Wiggins的最小熵反褶积的基础上，Levy和Oldenburg等人对地震道进行相位旋转，利用最大方差模准则评价不同相位时地震道偏离高斯性的程度，提出一种自动相位校正的方法 (APC)直接从地震数据中估计子波相位。但是，该方法并未考虑到反射系数的干涉对地震数据相位的影响，而且当地震数据缺少低频信息时该方法也很难取得较好的效果。

考虑到峰度准则具有放大强反射和压制弱反射等性质，Lu和Liu等人通过对归一化的振幅取对数的处理方式对峰度准则进行改进，提出了一种改进的峰度准则。L1范数作为一种良好的稀疏准则被广泛应用于地震数据处理的各个环节当中，随后Scargle等将L1范数用来评价相位提取的结果。 Sacchi研究了Sech准则和Cauchy准则，并分别将其作为稀疏约束引入到稀疏脉冲反褶积中。在测井资料丰富时，采用零相位地震道作为标准地震道，可以采用相似系数准则估计相位，Taner等人将改进的相似系数准则应用于实际地震数据的相位校正中且取得了较好的效果。刘玉金和李振春等人利用局部相似度准则估计出井旁道的相位属相，再将其插值到整个地震剖面，从而实现对非稳态数据进行了相位校正处理。周兴元对常相位校正后的结果的包络进行了分析比较，且分析了子波剩余相位不是常数时进行常相位校正的意义。刘金俊等人阐述了子波零相位化、反褶积和地震分辨率之间的关系。刘俊州、马铭等基于偏斜度准则进行子波相位估计，分析了带宽、噪声以及反射系数模型等对结果的影响。随后Yu和Yuan等人分析对比了峰度、Parsimony、柯西和改进柯西等准则在子波常相位估计当中的应用。

值的注意的是，反射系数的干涉也会引起地震数据的相位变化，这种干涉效应通常与地震子波相位耦合在一块难以分离，进而影响子波真实相位的估计。上述现有技术中的子波相位估计方法，是根据地震记录或测井数据直接估计子波相位，并未考虑到反射系数的干涉效应对子波相位估计的影响，影响了子波真实相位的估计。

当我们应用构造的不同相位的子波进行反演或反褶积处理时，在地震信号波形匹配的前提下，会改变反演的反射系数序列的相位谱，会导致反演的反射系数其在时域的能量分布发散。换句话说，当使用不准确相位的地震子波去反演时，相位的差异导致无法采用单一的反射系数脉冲刻画该波形，而是需要多个脉冲才能够匹配上已知的信号，因此得到的反射系数反演结果会呈现不稀疏的现象。地震记录表现出来的相位是反射系数干涉与真实子波相位相互耦合的‘伪’相位。因此，如何去解耦这种干涉效应进而准确地估计真实子波相位是获得地下真实反射脉冲图像的前提条件，也是获得高分辨率处理结果中不可或缺的一步。

发明内容

本发明的特征和优点在下文的描述中部分地陈述，或者可从该描述显而易见，或者可通过实践本发明而学习。

本发明提供一种基于加权L1范数稀疏准则的子波相位估计方法，包括：

S1、获取地震数据，并给定初始地震子波；

S2、构建不同相位的子波群；

S3、利用构建的子波群进行基于稀疏贝叶斯学习的稀疏反演，得到反演结果；

S4、采用加权L1范数稀疏准则量化评价所述反演结果的稀疏性，获取地震子波相位和对应的反射系数序列中至少一项。

在本发明的一个实施例中，所述步骤S3中，包括：假设需反演的反射系数序列服从分等级高斯分布，则对应的反射系数先验概率密度分布函数为：

其中，r表示反射系数，h＝[h₁,h₂,L,h_m]^T代表M个相互独立的超参数，每一个超参数分别控制其对应反射系数大小的先验信息；

表示反射系数r_m服从均值为0，方差为

的高斯分布，其中，m＝1，2，…，M。

在本发明的一个实施例中，所述步骤S3中，包括：假定噪声服从零均值、方差为

的正态分布(高斯分布)，且独立，则可以推导对应的似然函数为：

其中，d表示地震数据，W为对应相位的子波褶积矩阵。

在本发明的一个实施例中，所述步骤S4具体包括：

将反演结果进行能量归一化后代入到加权L1范数函数，得到相对应的扫描曲线；

根据所述扫描曲线求所述似然函数的最小值对应的相位，即可得到地震子波相位和对应的反射系数序列。

在本发明的一个实施例中，所述加权L1范数函数为：

其中，

表示对反射系数序列进行能量归一化，

表示反射系数序列的“能量”，

表示权重，σ＞0是一个稳定因子，0＜p≤2。

本发明提供一种基于加权L1范数稀疏准则的子波相位估计系统，包括：

初始子波给定单元，用于获取地震数据，并给定初始地震子波；

子波群构建单元，用于构建不同相位的子波群；

反演单元，用于利用构建的子波群进行基于稀疏贝叶斯学习的稀疏反演，得到反演结果；

评价单元，用于采用加权L1范数稀疏准则量化评价所述反演结果的稀疏性，获取地震子波相位和对应的反射系数序列中至少一项。

在本发明的一个实施例中，所述反演单元，用于：假定噪声服从零均值、方差为

的正态分布(高斯分布)，且独立，则可以推导对应的似然函数为：

其中，d表示地震数据，W为对应相位的子波褶积矩阵。

在本发明的一个实施例中，所述评价单元，具体用于：

将反演结果进行能量归一化后代入到加权L1范数函数，得到相对应的扫描曲线；

根据所述扫描曲线求所述似然函数的最小值对应的相位，即可得到地震子波相位和/或对应的反射系数序列。

在本发明的一个实施例中，所述加权L1范数函数为：

其中，

表示对反射系数序列进行能量归一化，

表示反射系数序列的“能量”，

表示权重，σ＞0是一个稳定因子，0＜p≤2。

本发明提供一种计算机可读存储介质，所述计算机可读存储介质存储有计算机可执行的至少一个程序，其特征在于，所述至少一个程序被所述计算机执行时使所述计算机执行本发明任一实施例所述的方法中的步骤。

本发明基于反射系数稀疏假设，通过加权L1范数稀疏准则来评价稀疏反褶积结果的稀疏程度，进而估计得到地震子波相位的同时得到高分辨率图像结果，分离反射系数的干涉效应与真实子波相位相互耦合的影响。

本发明提供的基于加权L1范数稀疏准则的子波相位估计方法及系统，建立了加权L1范数稀疏准则与地震子波相位之间的关系，通过加权L1范数和稀疏贝叶斯学习的双重稀疏表示，分离了子波相位和反射系数的干涉耦合效应，进而可以提取出了高分辨率反射系数图像结果。

通过阅读说明书，本领域普通技术人员将更好地了解这些技术方案的特征和内容。

附图说明

下面通过参考附图并结合实例具体地描述本发明，本发明的优点和实现方式将会更加明显，其中附图所示内容仅用于对本发明的解释说明，而不构成对本发明的任何意义上的限制，在附图中：

图1A为两种地震记录。

图1B为图1A所示的两种地震记录对应的两种反射系数模型；

图1C为图1A所示的两种地震记录对应的不同相位的子波。

图2为本发明实施例的基于加权L1范数稀疏准则的子波相位估计方法的流程示意图。

图3为本发明实施例的基于加权L1范数稀疏准则的子波相位估计系统的结构示意图。

图4为加权L1范数准则的子波相位扫描曲线；其中，纵坐标的值都是对原函数值取对数后的结果。

图5为厚与低频成分对加权L1范数准则的影响；其中，横坐标为层厚；纵坐标表示对应层厚扫描准确的概率。

图6A为Marmousi模型测试中原始Marmousi褶积数据；

图6B为Marmousi模型测试中相位扫描曲线；

图6C为本发明的处理结果。

图7A为实际地震资料测试的原始地震数据；

图7B为相位扫描曲线；

图7C为本发明处理结果。

具体实施方式

褶积模型是对许多条件进行合理的假设后得到的符合实际情况的地震记录模型，其时间域的表达式为：

其中，s表示地震记录，r表示反射系数，

表示子波的相位。在方程 (1)中，

可以由一个初始零相位子波经过旋转任意相位得到，在频率域可以由下式实现：

其中，f表示频率，W(f)表示时域初始子波的傅里叶变换，

表示旋转相位后得到的子波的傅里叶变换结果，对其进行反傅里叶变换可以得到时间域的常相位子波

在时间域

可以由下式获得：

其中，w₀(t)表示初始子波，H[g]表示希尔伯特变换，imag表示对希尔伯特变换结果取虚部。

当我们在进行反褶积或反射系数反演处理时，地震信号在处理前后是匹配的。从方程(1)中可以看出，当使用不准确相位的地震子波去反演时，这种相位的差异会导致反演得到的反射系数结果会出现连续同极性脉冲。换句话说，当有剩余相位存在时，反演结果需要比正确的反射系数序列更多的脉冲才能够匹配上已知的信号，也就是会导致反演的反射系数其在时域的能量分布发散，导致不稀疏的反演结果；另一方面，反射系数的干涉同时也会引起如图1A至1C所示的地震数据的振幅和相位变化，这种反射系数的干涉势必会影响到子波相位的估计。因此想要得到更加精准的子波相位就需要解耦反射系数干涉与子波真实相位这两者对地震数据的影响。

本发明中，将方程(1)写成矩阵-向量形式为：

其中，

表示对应相位的子波褶积矩阵。

本发明通过相位扫描的策略，对方程(4)我们采用SBL(稀疏贝叶斯学习)进行稀疏反演，将不同相位子波对应的反演结果看作子波相位校正结果的集合，提出采用加权L1范数准则来量化反演结果的稀疏度。

如图2所示，本发明实施例提供一种双重稀疏约束盲反褶积提高分辨率方法，包括步骤：

S1、获取地震数据，并给定一个初始地震子波；

S2、构建不同相位的子波群；

在具体实施时，可以给定相位扫描范围，采用相位扫描的策略完成不同相位子波群的建立；

S3、利用构建的子波群进行基于稀疏贝叶斯学习的稀疏反演，得到反演结果；

S4、采用加权L1范数稀疏准则量化评价所述反演结果的稀疏性，获取地震子波相位和对应的反射系数序列中至少一项。

更具体地，先将反演结果进行能量归一化，然后代入到加权L1范数函数，得到相对应的扫描曲线。根据该扫描曲线求稀疏反演对应的似然函数的最小值对应的相位，即可得到地震子波相位和对应的反射系数序列。

在实际处理中可以对扫描曲线进行光滑处理以减少奇异值的影响和单道反演不稳定的影响。

具体地，在进行步骤S3时：

由于实际资料通常含有噪声，因此对方程(4)加入噪声，有：

Wr+n＝d (5)

其中，n代表噪声，d表示地震数据。由于地震子波是带限的，子波矩阵W是个病态矩阵，因此直接求解方程(5)是不适定的。为此，我们应用稀疏贝叶斯学习反演方法进行稀疏反演，本实施例中稀疏反演为盲反褶积处理。

在进行稀疏贝叶斯学习反褶积时，假定噪声服从零均值、方差为

的正态分布(高斯分布)，且独立，则可以推导对应的似然函数为

为了获得与地质假设一致的解，我们假设需反演的反射系数序列服从分等级高斯分布，则对应的反射系数先验概率密度分布函数为：

其中，h＝[h₁,h₂,L,h_m]^T代表M个相互独立的超参数，每一个超参数分别控制其对应反射系数大小的先验信息。

表示反射系数r_m服从均值为0，方差为

的高斯分布，其中，m＝1，2，…，M。

根据贝叶斯理论，在d，h和σ²已知的条件下，联立公式(6)和(7)，则有反射系数r的后验概率密度分布为：

其中，C表示一个常数；Σ表示协方差；μ表示均值；H＝diag(h₁,h₂,L,h_M)。为了获得盲反褶积结果，首先要估计超参数h和σ²的最佳值。根据贝叶斯框架，超参数的边缘分布通过下式计算：

p(d|h,σ²)＝-2∫p(d|r,σ²)p(r|h)dr＝(2π)^M|Q|exp(d^TQ^-1d) (10)

其中，

对公式(10)求最小即可求出超参数h。我们使用了一种快速的迭代算法-相关向量机来获得盲反褶积结果。实质上，通过公式(10)每次估计出来的超参数h为稀疏脉冲的位置，而通过公式(9)求出μ则对应于稀疏脉冲的大小。

在步骤S4中，可以通过最小化不同相位子波对应反演结果的加权L1 范数函数值，可以得到真实的子波相位以及对应的反射系数序列，加权L1 范数稀疏准则如下式所示：

其中，

表示对反射系数序列进行能量归一化，

表示反射系数序列的“能量”，

表示权重，σ＞0是一个稳定因子，0＜p≤2。虽然反演前后数据能量一致以及不同相位子波的能量是一致的，但由于地震记录的带限性质、反射系数反演的过程又是一个非因果过程，而且子波褶积矩阵是一个病态矩阵，导致反演过程本质上就有多解性以及算法的计算误差等因素会导致反演结果的能量稍有误差。因此，需要对反射系数反演结果进行能量归一化处理以去除反演结果能量不一致对我们应用稀疏准则函数的影响。

本发明提供一种计算机可读存储介质，所述计算机可读存储介质存储有计算机可执行的至少一个程序，其特征在于，所述至少一个程序被所述计算机执行时使所述计算机执行本发明任一实施例所述的方法中的步骤。

如图3所示，本发明提供一种基于加权L1范数稀疏准则的子波相位估计系统，包括：初始子波给定单元10、子波群构建单元20、反演单元30、评价单元40。

初始子波给定单元10用于获取地震数据，并给定初始地震子波；子波群构建单元20与所述初始子波给定单元10相连，子波群构建单元20用于构建不同相位的子波群；反演单元30与子波群构建单元20相连，用于利用构建的子波群进行基于稀疏贝叶斯学习的稀疏反演，得到反演结果；评价单元40与反演单元30相连，评价单元40用于采用加权L1范数稀疏准则量化评价所述反演结果的稀疏性，获取地震子波相位和对应的反射系数序列中的至少一项。

在具体实施时，反演单元用于实现上述步骤S3，在此不再赘述，评价单元用于实现上述步骤S4，在此不再赘述。

本发明将不同子波相位的稀疏反演或反褶积结果看作子波相位校正结果的集合，采用加权L1范数稀疏准则量化评价稀疏反演或反褶积结果的稀疏性，可辨识出真实或准确的地震子波相位

同时得到对应的反射系数序列。

本发明采用加权L1范数来量化反演结果的稀疏性与子波真实相位的关系，该准则函数能够平衡稀疏准则函数中L1和L0范数在幅值依赖性上的差异同时兼顾弱反射和强反射对函数值的贡献。对于加权L1范数而言，通过自适应调节权重值来增强弱信号对稀疏表示的贡献同时抑制强信号的贡献，从而加强弱反射对函数值的贡献。这一特性可以提高本方法在相位估计方面的稳定性，同时减少参数、反射系数模型等的影响，实现了子波相位估计和稀疏盲反褶积的同步处理。

下面辅以具体的实例进一步阐述本发明：

为了去除偶然性与测试方法的健壮性，我们设计了1000个随机反射系数模型，每个模型有30个非零脉冲，非零脉冲的幅值保持在-0.3-0.3之间，非零脉冲的间隔不超过20ms，真实子波是主频为30Hz的零相位Richer子波，采样间隔为1ms，总时间长度为1s。我们展示了基于加权L1范数准则的子波相位扫描结果(如图4所示)。图4中1000次扫描结果准确的概率为99.2％(扫描到准确相位的次数除以总的模型数)。从图中也可以看出，在无噪情况下利用加权L1范数准则估计的子波相位稳定且准确率高，而且扫描曲线的凸性也较好，不易受反射系数模型的影响。

为了测试了在不同层厚以及数据缺失低频的情况时该方法的性能，我们重新设置了100个随机反射系数模型，每个模型有30个非零脉冲，脉冲的幅值控制在-0.3-0.3之间，其中最小脉冲幅值为0.01(相当于模型包含相对的强反射与弱反射)，子波是峰值频率为30Hz的零相位子波，采样间隔为1ms，模型总时间长度为1s。随后，我们分别将100个反射系数模型非零脉冲的间隔从1ms一直增大到30ms。从图5中可以看出，当非零脉冲间隔很小时(即反射层过薄，例如小于12ms)，相位估计准确的概率是较低的，但随着时间间隔的增加准确的概率将达到100％；而且当缺失10Hz低频数据时，依然可以保持很稳定准确的估计结果；当缺失低频数据过多整体曲线趋势会往右移(如图5黑色箭头所示)，也就是说需要更厚的地层才能有好的估计效果。

为了进一步测试方法的有效性，我们采用MArmousi模型数据进行测试 (如图6A所示)。其中原始子波为峰值频率为30Hz的RiCker子波，采样间隔为1ms。从图6B可以看出，对于较为复杂的模型，该方法仍然可以较为准确地指示出正确相位。对比图6A和6C我们可以看出，经过处理以后的结果同相轴更加清晰，而且刻画出了更多的细节信息(如圆圈处所示)。

最后我们将该发明应用于实际地震资料(图7A)中。从图7B中可以看出，扫描曲线具有良好的凸性，而且每一道的拐点都较为一致，横向一致性较好。图7C为采用扫描出来的子波相位进行反演的结果。对比图7A和 7C，我们可以看出，经过处理以后的地震剖面不仅同相轴变细了而且细节信息更加丰富了。

本发明基于稀疏假设，利用加权L1范数来量化稀疏贝叶斯反射系数反演结果的稀疏性与子波相位之间的关系，通过双重稀疏约束进行盲反褶积处理，同时估计得到子波的相位以及反褶积结果。其优势主要有：发现了反射波干涉引起的相位变化与子波真实相位之间的内禀规律，设计了加权 L1范数函数去量化了这种内禀关系以达到分离两者的目的，从而更准确地估计得到子波相位。而且提高了子波相位估计的稳定性和准确性，同时减少了参数以及反射系数模型等的影响。

以上参照附图说明了本发明的优选实施例，本领域技术人员不脱离本发明的范围和实质，可以有多种变型方案实现本发明。举例而言，作为一个实施例的部分示出或描述的特征可用于另一实施例以得到又一实施例。以上仅为本发明较佳可行的实施例而已，并非因此局限本发明的权利范围，凡运用本发明说明书及附图内容所作的等效变化，均包含于本发明的权利范围之内。

Claims (10)

1.一种基于加权L1范数稀疏准则的子波相位估计方法，其特征在于，包括：

S1、获取地震数据，并给定初始地震子波；

S2、构建不同相位的子波群；

S3、利用构建的子波群进行基于稀疏贝叶斯学习的稀疏反演，得到反演结果；

S4、采用加权L1范数稀疏准则量化评价所述反演结果的稀疏性，获取地震子波相位和对应的反射系数序列中至少一项。

2.根据权利要求1所述基于加权L1范数稀疏准则的子波相位估计方法，其特征在于，所述步骤S3中，包括：假设需反演的反射系数序列服从分等级高斯分布，则对应的反射系数先验概率密度分布函数为：

其中，r表示反射系数，h＝[h₁,h₂,L,h_m]^T代表M个相互独立的超参数，每一个超参数分别控制其对应反射系数大小的先验信息；

表示反射系数r_m服从均值为0，方差为

的高斯分布，其中，m＝1，2，…，M。

3.根据权利要求1所述基于加权L1范数稀疏准则的子波相位估计方法，其特征在于，所述步骤S3中，包括：假定噪声服从零均值、方差为

的正态分布(高斯分布)，且独立，则可以推导对应的似然函数为：

其中，d表示地震数据，W为对应相位的子波褶积矩阵。

4.根据权利要求3所述基于加权L1范数稀疏准则的子波相位估计方法，其特征在于，所述步骤S4具体包括：

将反演结果进行能量归一化后代入到加权L1范数函数，得到相对应的扫描曲线；

根据所述扫描曲线求所述似然函数的最小值对应的相位，即可得到地震子波相位和对应的反射系数序列。

5.根据权利要求1所述基于加权L1范数稀疏准则的子波相位估计方法，其特征在于，所述加权L1范数稀疏准则为：

其中，

表示对反射系数序列进行能量归一化，

表示反射系数序列的“能量”，

表示权重，σ＞0是一个稳定因子，0＜p≤2。

6.一种基于加权L1范数稀疏准则的子波相位估计系统，其特征在于，包括：

初始子波给定单元，用于获取地震数据，并给定初始地震子波；

子波群构建单元，用于构建不同相位的子波群；

反演单元，用于利用构建的子波群进行基于稀疏贝叶斯学习的稀疏反演，得到反演结果；

评价单元，用于采用加权L1范数稀疏准则量化评价所述反演结果的稀疏性，获取地震子波相位和对应的反射系数序列中至少一项。

7.根据权利要求6所述基于加权L1范数稀疏准则的子波相位估计系统，其特征在于，所述反演单元，用于：假定噪声服从零均值、方差为σ_n2的正态分布(高斯分布)，且独立，则可以推导对应的似然函数为：

其中，d表示地震数据，W为对应相位的子波褶积矩阵。

8.根据权利要求7所述基于加权L1范数稀疏准则的子波相位估计系统，其特征在于，所述评价单元，具体用于：

将反演结果进行能量归一化后代入到加权L1范数函数，得到相对应的扫描曲线；

根据所述扫描曲线求所述似然函数的最小值对应的相位，即可得到地震子波相位和对应的反射系数序列。

9.根据权利要求1所述基于加权L1范数稀疏准则的子波相位估计系统，其特征在于，所述加权L1范数稀疏准则为：

其中，

表示对反射系数序列进行能量归一化，

表示反射系数序列的“能量”，

表示权重，σ＞0是一个稳定因子，0＜p≤2。

10.一种计算机可读存储介质，所述计算机可读存储介质存储有计算机可执行的至少一个程序，其特征在于，所述至少一个程序被所述计算机执行时使所述计算机执行上述权利要求1～5任一项所述的方法中的步骤。

Images