CN113077386A - 基于字典学习和稀疏表征的地震资料高分辨率处理方法 - Google Patents

Abstract

本发明公开了一种基于字典学习和稀疏表征的地震资料高分辨率处理方法，包括以下步骤：S1、进行井震标定，提取地震子波；S2、建立联合样本；S3、利用K‑SVD算法训练得到联合字典；S4、将实际采集到的地震数据做归一化处理；S5、将归一化后的地震数据在地震数据字典下进行稀疏表示，得到稀疏表示系数；S6、利用反射系数字典和稀疏表示系数进行稀疏重构，得到归一化的反射系数；S7、对归一化的反射系数进行反归一化处理，得到地震数据对应的反射系数。本发明将测井反射系数和井旁地震数据结合起来进行字典的训练，然后进行基于联合字典的反射系数反演，不仅可以提高地震资料分辨率，并且可以通过字典学习的方法引入测井反射系数特征。

Description

基于字典学习和稀疏表征的地震资料高分辨率处理方法

技术领域

本发明属于地球物理勘探解释技术领域，是一种通过测井数据(反射系数)与井旁道地震数据联合训练字典，实现井数据反射系数和地震数据联合特征的获取。在联合字典对地震数据稀疏表征基础上，获取反射系数反演结果，以实现井数据特征约束下的地震数据稀疏反褶积。

背景技术

地震勘探一直是地球物理勘探中探测油气藏最重要、最有效的一种方法。所谓的地震勘探技术就是利用地震波的响应推测地下介质的属性和地层构造，其主要工作流程为：通过工业爆破或者打桩等人工方式引发地震从而产生地壳振动形成地震波，地震波沿着地层向下传播，当遇到弹性介质分界面时会产生反射，此时地面上的精密仪器会记录下地上各采样点的震动情况，通过分析记录的采样点的特征，可以帮助我们寻找油气藏的位置、预测油气藏的种类、判断地层岩性特征等。

地震反演技术是地震勘探的重要手段之一。地震反演是利用地表记录地震数据，通过前人经验构建目标函数，将钻井资料和测井数据作为约束，反映地下空间构造和地下介质的物理特性并进行成像的过程。地震反演促进了地震岩性分析技术和地震流体分析技术的发展，并为油气藏的勘探开发提供了重要的依据。

地震资料的分辨率一直是地震勘探领域研究的重点。无论是叠前还是叠后数据，采用反褶积或者反射系数反演来提高地震资料的垂向分辨率，是目前最常用的方法。在反褶积方面主要有线性反褶积和非线性反褶积。线性反褶积在原有地震数据频带范围内通过扩展高频分量来提高分辨率，并不能创造出地震数据频带以外的高频分量，因此其提高分辨率的能力有限；非线性反褶积通过加入各种稀疏约束项，来创造出地震数据原始频带以外的高频，可以较大幅度的提高地震资料的分辨率。

稀疏脉冲反褶积是目前最成熟也是最常用的非线性反褶积方法，其基本思路是：对反射系数施加L₁范数稀疏约束。稀疏脉冲反褶积的基本假设是地层的强反射系数是稀疏的，即地层反射系数由一系列叠加于高斯背景的强轴组成。但是，真实的地层反射系数并非完全满足稀疏脉冲反褶积的假设条件。比如，对大套强反射层周围存在薄层弱反射时，稀疏脉冲反褶积只能重构出强反射系数，对弱反射有压制作用，并且很难识别薄互层，这显然不符合实际地层特性。稀疏脉冲反褶积对地震子波的依赖性很强，然而地震子波的准确提取一直是地震反演的难点。

为了将全频带的测井数据特征引入地震反演中，我们考虑将测井反射系数与井旁道地震数据联合起来，通过字典学习的方法来获取联合特征。在联合字典对地震数据稀疏表征基础上，获取反射系数反演结果，以实现井数据特征约束下的地震数据稀疏反褶积。

Robinson褶积模型建立了地震数据s(t)，反射系数r(t)和地震子波w(t)三者之间的关系。为了模拟真实地震数据在采集过程中伴有的随机噪声的干扰，可以在褶积模型中加入噪声n(t)。上述过程可以表示为以下形式：

s(t)＝w(t)*r(t)+n(t)

式中，*表示卷积运算符。由于在数字信号处理中，我们习惯用矩阵-向量的形式简化卷积运算，因此上式可以表示为以下矩阵-向量形式：

S＝WR+N

其中，S是地震数据，R是反射系数，N是随机加性噪声，W表示子波褶积矩阵，其具体形式如下：

k表示子波长度。

上述褶积过程可以看作高频反射系数经过滤波算子滤波之后变为低频的地震数据过程，该滤波算子就是震源子波。稀疏脉冲反褶积过程是褶积过程的逆过程，但是由于该问题的多解性，所以在求解反问题时，在传统的最小二乘反褶积的基础上加上一个稀疏先验约束项，通过求解L₂范数和L₁范数之和形式的目标函数得到高分辨率的反演结果。稀疏脉冲反褶积的目标函数为：

基于稀疏脉冲反褶积的基本原理，学者们提出了各种稀疏脉冲反褶积方法：最小熵反褶积通过搜索数量最少的强脉冲使得其正演结果与地震记录保持一致，Debeye和Reil于1990年首次提出了基于p范数的稀疏脉冲反褶积，揭开了研究约束稀疏脉冲反褶积方法的新篇章，Sacchi分别在1997年和1999年使用Huber分布、Cauchy分布和Sech分布来恢复稀疏反射序列，Meng在2005年使用Cauchy准则作为约束项，并利用预条件共轭梯度法估计地震子波和稀疏反射系数序列，Danilo RVelis在2008年提出了随机稀疏脉冲反褶积，并指出稀疏脉冲反褶积本质上是从含噪记录中寻找一些脉冲，并确定它们的位置和振幅，快速模拟退火算法被应用于寻找脉冲的位置，而脉冲的振幅可以由线性最小二乘法来确定。

稀疏字典学习算法源于图像处理领域。人们尝试着将图像压缩，即通过某种方法去除掉图像的冗余信息从而来简化图像的表达方式，这有利于提升存储空间、提高处理效率以及降低处理成本等。这种用尽可能少的资源来表示尽可能多的信息的方法，通常被称为稀疏表示方法。

最早的稀疏表示方法是将图像表示成一组正交基向量的线性组合，不过这些正交基都是固定的，不能自适应性地表示图像。Meyer提出了自适应多尺度几何分析方法，该方法提出的字典是可以随图像的变化而变化的，克服了正交基有限方向的缺点。20世纪90年代，Mallat和Zhang发现当原子个数远远大于信号维数时，所训练出的字典将包含更多有效的信息，并将这样的字典称之为过完备字典，同时还提出了匹配追踪算法以解决过完备字典稀疏表示的问题。匹配追踪算法仅与最新的原子正交，不能保证与所有选出的原子正交，所以该算法是一种不收敛的贪婪算法。Pati等在匹配追踪算法的基础上提出了正交匹配追踪算法，该方法在每一步的分解中都会对已选择的原子进行正交化处理，在提升了收敛速度的同时保证每一步都选择了不同的原子。在字典训练方面，Aharon等提出了非常著名的快速奇异值分解算法，该方法在更新字典时保证每次只更新一个原子和其对应的稀疏系数，且更新一个原子时其他原子保持不变，大大加快了迭代的收敛速度。传统的字典训练方法对样本数据进行集中训练，不适合处理经常更新或数量过多的样本数据。Mairal等提出了能解决上述问题的在线字典学习方法，该方法每次只处理部分样本，能够在新样本加入时动态地更新字典，提升了训练字典的速度，增强了字典的适应性。

由于地下地层的不连续性，地震反射系数序列往往表现出稀疏特征，同时地下弹性参数表现出块状特征，于是地球物理学家考虑将稀疏表示的思想引入到地震勘探中。基于稀疏表示的地震反演主要是将测井数据作为训练样本学习得到含有地层信息的稀疏字典，并将稀疏字典作为反演目标函数的先验信息约束，最后利用最优化方法和稀疏求解算法反解出需要的地层弹性参数。

发明内容

本发明的目的在于克服现有技术的不足，提供一种将测井反射系数和井旁地震数据结合起来，进行字典的训练；然后进行基于联合字典的反射系数反演，以实现井数据特征约束下的地震数据稀疏反褶积的基于字典学习和稀疏表征的地震资料高分辨率处理方法。本发明不仅可以提高地震资料分辨率，并且可以通过字典学习的方法引入测井反射系数特征，克服了稀疏脉冲反褶积压制弱反射以及难以识别薄互层的缺点，减弱了反演结果对子波的依赖程度。

本发明的目的是通过以下技术方案来实现的：基于字典学习和稀疏表征的地震资料高分辨率处理方法，包括以下步骤：

S1、进行井震标定，提取地震子波；

S2、将所有井的测井反射系数与地震子波进行合成，得到合成地震数据；将每口井的合成地震数据与实际采集到的井旁道地震数据进行比较，将二者相似度达到85％以上的井的井旁道地震数据作为地震数据样本，建立联合样本{Y_S,Y_R}，其中Y_S和Y_D分别是井旁道地震数据样本和测井反射系数样本；

S3、利用K-SVD算法训练得到联合字典{D_S,D_R}，D_S为地震数据字典，D_R为反射系数字典；

S4、将实际采集到的地震数据S做归一化处理；

S5、将归一化后的地震数据在字典D_S下进行稀疏表示，得到稀疏表示系数α；

S6、利用反射系数字典D_R和稀疏表示系数α进行稀疏重构，得到归一化的反射系数；

S7、对归一化的反射系数进行反归一化处理，得到地震数据对应的反射系数。

进一步地，所述步骤S3中，联合字典{D_S,D_R}的训练方法为：

λ是正则化参数；A代表是字典学习过程中的表示系数；|| ||₁表示L₁范数；

联合上两式得到目标函数如下：

其中，N和M分别是地震数据样本和反射系数样本的维度，在联合字典中，取N＝M。

进一步地，所述步骤S5具体实现方法为：稀疏表示的模型表示为以下形式：

r＝D_Rα，s＝D_Sα

其中，r和s分别表示归一化后未知的反射系数和地震数据；α是一个大部分元素为零或者接近于零的向量，称为稀疏表示系数；本发明需要从未知反射系数的地震数据s中恢复出相对应的反射系数r，首先需要通过求解以下目标函数

来得到稀疏表示系数α：

||α||₀为L₀范数，表示向量α中非零元素的个数；在稀疏编码中，α是足够稀疏的，因此用L₁范数近似L₀范数：

通过标准线性规划的方法来对上式进行求解，得到最优解α。

进一步地，所述步骤S6具体实现方法为：利用反射系数字典D_R和步骤S5得到的稀疏表示系数α，对未知的井旁道地震数据进行稀疏重构，得到归一化反射系数：

r＝D_Rα。

本发明的有益效果是：本发明将测井反射系数和井旁地震数据结合起来，进行字典的训练；然后进行基于联合字典的反射系数反演，以实现井数据特征约束下的地震数据稀疏反褶积。本发明不仅可以提高地震资料分辨率，并且可以通过字典学习的方法引入测井反射系数特征，克服了稀疏脉冲反褶积压制弱反射以及难以识别薄互层的缺点，减弱了反演结果对子波的依赖程度。

附图说明

图1为本发明的基于字典学习和稀疏表征的地震资料高分辨率处理方法的流程图；

图2为30Hz Ricker子波；

图3为单道合成数据的联合字典；

图4为采用本发明出的方法和稀疏脉冲反褶积分别对地震数据做处理，并与真实反射系数的对比图；

图5为多道合成数据的联合字典；

图6为本发明处理效果与传统的稀疏脉冲反褶积方法作对比图；

图7为抽取任意三道不参与字典训练的地震道进行结果对比的对比图；

图8为实际资料处理频谱分析图；

图9为实际资料的联合字典；

图10为不同方法对同一剖面的数据处理的对比图。

具体实施方式

联合字典首先应用于图像超分辨率问题中，通过建立低分辨率图像和高分辨率图像的联合字典，然后将低分辨率图像在低分辨率字典下进行稀疏表示，得到表示系数，然后利用该表示系数与高分辨率字典可以重构得到高分辨率图像。由于低分辨率图像可以看做对高分辨率图像做一个模糊操作，所以地震数据和反射系数之间关系与高低分辨率图像之间的这种关系很相似。受该思想的启发，我们设想采用联合字典的思想从地震数据中重构得到反射系数。

想要利用本发明方法获得理想反演结果，最关键是训练得到能够用于稀疏表示的联合字典。

假设训练样本数据集为Y＝{y₁,y₂,…,y_n},通过字典学习的方法，比如KSVD，通过求解以下目标函数，训练得到过完备字典D：

在本发明中，我们需要训练两个字典D_S和D_R。

下面结合附图进一步说明本发明的技术方案。

如图1所示，本发明的基于字典学习和稀疏表征的地震资料高分辨率处理方法，包括以下步骤：

S1、进行井震标定，提取地震子波；

S2、将所有井的测井反射系数与地震子波进行合成，得到合成地震数据；将每口井的合成地震数据与实际采集到的井旁道地震数据进行比较，将二者相似度达到85％以上的井的井旁道地震数据作为地震数据样本，建立联合样本{Y_S,Y_R}，其中Y_S和Y_D分别是井旁道地震数据样本和测井反射系数样本；

S3、利用K-SVD算法训练得到联合字典{D_S,D_R}，D_S为地震数据字典，D_R为反射系数字典；联合字典{D_S,D_R}的训练方法为：

λ是正则化参数；A代表是字典学习过程中的表示系数；|| ||₁表示L₁范数；

联合上两式得到目标函数如下：

其中，N和M分别是地震数据样本和反射系数样本的维度，在联合字典中，取N＝M。

S4、将实际采集到的地震数据S做归一化处理；

S5、将归一化后的地震数据在字典D_S下进行稀疏表示，得到稀疏表示系数α；

在地震勘探中，同一工区的地震记录有着相似的形态结构特征，所以可以认为测井资料在一定程度上代表着整个工区的地层特性，而稀疏表示的主要思想是以尽可能少的资源来表示尽可能多的信息，因此我们考虑利用较少的测井资料来表征整个工区的地层参数，将测井反射系数和相应的井旁道地震数据作为训练联合字典的样本，通过字典学习的方法得到可以表示整个工区地层特征的字典。具体实现方法为：稀疏表示的模型表示为以下形式：

r＝D_Rα，s＝D_Sα

其中，r和s分别表示归一化后未知的反射系数和地震数据；α是一个大部分元素为零或者接近于零的向量，称为稀疏表示系数；本发明需要从未知反射系数的地震数据s中恢复出相对应的反射系数r，首先需要通过求解以下目标函数

来得到稀疏表示系数α：

||α||₀为L₀范数，表示向量α中非零元素的个数；在稀疏编码中，α是足够稀疏的，因此用L₁范数近似L₀范数：

通过标准线性规划的方法来对上式进行求解，得到最优解α。

S6、利用反射系数字典D_R和稀疏表示系数α进行稀疏重构，得到归一化的反射系数：

r＝D_Rα。

S7、对归一化的反射系数进行反归一化处理，得到地震数据对应的反射系数。

通过合成单道地震数据以及利用经典的Marmousi模型合成多道地震数据，来验证本发明提出方法的有效性。

(1)单道合成数据

随机生成反射系数序列，然后用30Hz的Ricker子波(如图2)来生成地震记录，为了模拟实际地震数据，加入20dB的随机噪声。利用字典学习方法训练联合字典，如图3所示，图3(a)为地震数据字典部分原子；图3(b)为反射系数字典部分原子。然后采用本文提出的方法和稀疏脉冲反褶积分别对地震数据做处理，并与真实反射系数做对比，其结果如图4所示，图中，(a)为真实反射系数；(b)为含噪地震数据；(c)为本发明得到的反射系数；(d)为真实反射系数与本发明结果的误差；(e)为稀疏脉冲反褶积得到的结果；(f)为真实反射系数与稀疏脉冲反褶积结果的误差。

从图4的中可以看出，采用本发明的方法可以重构出弱反射系数，而稀疏脉冲反褶积对弱反射有压制作用；从蓝色虚线椭圆和黑色椭圆中可以看出，本文方法对薄层反射系数能够很好地重构。从两者的误差图中可以看出，采用本发明的方法误差明显小于稀疏脉冲反褶积的误差。

(2)多道合成数据

为了将本发明提出的联合字典方法与稀疏脉冲反褶积方法作进一步的对比，我们采用经典的Marmousi模型来做实验。我们利用30Hz的Ricker子波(图2)与Marmousi反射系数模型来合成地震数据，为了模拟真实地震数据采集过程中的由于采集设备、信号传输以及数据预处理等环节中的噪声干扰，在合成地震数据中加入信噪比为20dB的随机噪声，然后对该地震数据进行高分辨率处理。

在Marmousi模型中，共有377道地震数据，为了使训练得到的字典原子包含有整个工区的地震数据特征和反射系数特征，在本实验中我们抽取样本数据的取样间隔是10，也就是说均匀取出该剖面中的37道地震数据以及相对应的反射系数作为训练字典的样本。为了使训练字典的样本尽可能的多，训练得到的字典原子特征尽可能的丰富，我们对该样本数据做取小块操作。

接下来我们利用KSVD方法来训练联合字典，然后根据矩阵行分块的思想，分别得到含有该工区地震数据特征和反射系数特征的地震数据字典D_S和反射系数字典D_R，结果如图5所示，图中，(a)为地震数据字典部分原子；(b)为反射系数字典部分原子。

在字典训练完成之后，将该工区地震数据S，利用式

进行稀疏表示得到稀疏表示系数α，接着在反射系数字典D_R下，根据公式

重构反射系数。为了验证本发明的处理效果，与传统的稀疏脉冲反褶积方法作对比，处理结果如下图6所示，图中，(a)为含噪地震数据；(b)为Marmousi反射系数剖面；(c)为本发明处理结果；(d)为稀疏脉冲反褶积方法处理结果。从图6可以看出，本发明处理后的结果能够提升地震数据分辨率，而且比稀疏脉冲反褶积方法的处理结果更接近真实反射系数。从红色矩形和蓝色椭圆中可以看出，采用本发明处理的结果具有更多的细节信息，可以更真实的反映地下地层结构，以方便地质解释人员进行地层地质解释。

为了进一步比较两者的处理效果，我们抽取任意三道不参与字典训练的地震道进行结果对比。在本实验中，我们抽取第100道、200道和第300道数据进行对比分析，其结果如图7所示，图中，(a)为第100道数据对比；(b)为第200道数据对比；(c)为第300道数据对比。从图7中可以看出，稀疏脉冲反褶积方法只能反演出强反射系数，对弱反射有压制作用，稀疏脉冲反褶积结果与真实反射系数相差较大。图7中，True Reflectivity为真实反射系数，Reflectivity of JD为联合字典(Joint Dictionary)得到的反射系数，Reflectivity ofSSD为稀疏脉冲反褶积(Sparse-spike Deconvolution)得到的反射系数。

我们将两种方法的处理结果与真实反射系数做误差分析以及相关性分析，结果如下表1所示：

表1 RMSE和相关性分析

	本文方法	稀疏脉冲反褶积方法
			均方根误差	0.0087	0.0121
相似性	0.8087	0.4303

从表1结果对比中可以看出，不管是均方根误差，还是相似度，本文提出的联合字典的方法均比稀疏脉冲反褶积方法效果好。

图8为实际资料处理频谱分析，图中，Seismic Data为地震数据频带，TrueReflectivity为真实反射系数频带，Reflectivity of JD为联合字典(Joint Dictionary)得到的反射系数频带，Reflectivity of SSD为稀疏脉冲反褶积(Sparse-spikeDeconvolution)得到的反射系数频带。从图中可以看出，利用本发明处理后的结果高频分量增加，频带比地震数据频带宽，分辨率增强。而稀疏脉冲反褶积方法的结果明显偏离真实反射系数频谱，更多脉冲的存在，使得频谱中出现了假(高)频现象。

为了验证本文方法的实际应用效果，采用东北某工区的实际叠后地震数据进行测试。该工区共有测井104口，由于本文方法对测井数据要求较高，所以我们选取测井数据与地震数据相关系数大于85％的19口井作为训练联合字典的样本数据。由于样本数据较少，为了尽可能多地提取样本数据的特征，我们对样本数据分小块操作，取小块长度为30，每隔1个采样点取下一个小块，通过联合字典的方法训练得到联合字典，如图9所示，图中，(a)为地震数据字典部分原子；(b)为反射系数字典部分原子。然后划分子字典并分别做列向量归一化处理。

由于训练字典的19口井数据来自整个工区，故我们可以认为这19口井训练得到的字典包含了该工区所有地震数据和反射系数的特征，所以训练得到的联合字典可以适用于整个工区任一剖面。在本实验中，我们选取某一剖面来验证方法的有效性。为了与稀疏脉冲反褶积方法作对比，我们采用稀疏脉冲反褶积方法处理同一剖面的数据，结果图10所示，图中，(a)为地震数据剖面；(b)为联合字典得到的反射系数剖面；(c)为稀疏脉冲反褶积得到的反射系数剖面。从图10可以看出，通过联合字典得到的反射系数剖面比稀疏脉冲反褶积的结果包含有更多的地层分布信息，而稀疏脉冲反褶积压制了弱反射，而且不能识别薄互层，只能识别出大套层，所以其反褶积效果将会影响地质解释人员对地层结构的判断。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (4)

1.基于字典学习和稀疏表征的地震资料高分辨率处理方法，其特征在于，包括以下步骤：

S1、进行井震标定，提取地震子波；

S2、将所有井的测井反射系数与地震子波进行合成，得到合成地震数据；将每口井的合成地震数据与实际采集到的井旁道地震数据进行比较，将二者相似度达到85％以上的井的井旁道地震数据作为地震数据样本，建立联合样本{Y_S,Y_R}，其中Y_S和Y_D分别是井旁道地震数据样本和测井反射系数样本；

S3、利用K-SVD算法训练得到联合字典{D_S,D_R}，D_S为地震数据字典，D_R为反射系数字典；

S4、将实际采集到的地震数据S做归一化处理；

S5、将归一化后的地震数据在字典D_S下进行稀疏表示，得到稀疏表示系数α；

S6、利用反射系数字典D_R和稀疏表示系数α进行稀疏重构，得到归一化的反射系数；

S7、对归一化的反射系数进行反归一化处理，得到地震数据对应的反射系数。

2.根据权利要求1所述的基于字典学习和稀疏表征的地震资料高分辨率处理方法，其特征在于，所述步骤S3中，联合字典{D_S,D_R}的训练方法为：

λ是正则化参数；A代表是字典学习过程中的表示系数；|| ||1₁表示L₁范数；

联合上两式得到目标函数如下：

其中，N和M分别是地震数据样本和反射系数样本的维度，在联合字典中，取N＝M。

3.根据权利要求1所述的基于字典学习和稀疏表征的地震资料高分辨率处理方法，其特征在于，所述步骤S5具体实现方法为：稀疏表示的模型表示为以下形式：

r＝D_Rα，s＝D_Sα

其中，r和s分别表示归一化后未知的反射系数和地震数据；α是一个大部分元素为零或者接近于零的向量，称为稀疏表示系数；本发明需要从未知反射系数的地震数据s中恢复出相对应的反射系数r，首先需要通过求解以下目标函数

来得到稀疏表示系数α：

||α||₀为L₀范数，表示向量α中非零元素的个数；在稀疏编码中，α是足够稀疏的，因此用L₁范数近似L₀范数：

通过标准线性规划的方法来对上式进行求解，得到最优解α。

4.根据权利要求1所述的基于字典学习和稀疏表征的地震资料高分辨率处理方法，其特征在于，所述步骤S6具体实现方法为：利用反射系数字典D_R和步骤S5得到的稀疏表示系数α，对未知的地震数据进行稀疏重构，得到归一化反射系数：

r＝D_Rα。

Images