CN113253350A - 一种基于联合字典的孔隙度反演方法 - Google Patents

Abstract

本发明公开了一种基于联合字典的孔隙度反演方法，包括以下步骤：S1、数据预处理：将测井数据进行分块处理；S2、将分块处理后的测井数据进行字典学习，得到阻抗字典和孔隙度字典；S3、将地震数据进行基于字典的孔隙度反演，得到地震数据的阻抗和孔隙度。本发明从测井数据出发，引入字典学习和稀疏表示理论，实现了孔隙度字典反演方法的构建，解决了传统方法对于复杂地质构造求解孔隙度适用性差的缺点，同时能够反演求解更多的测井物性参数。为提高测井数据在反演中的利用率提供了一种可行性方案，降低反演成本和反演周期，提高反演过程中测井数据的利用率。

Description

一种基于联合字典的孔隙度反演方法

技术领域

本发明属于地质统计学技术领域，特别涉及一种基于联合字典的孔隙度反演方法。

背景技术

储层物性参数是进行储层预测与评价的重要基础，孔隙度作为最基本的储层物性参数之一，其准确估计可以为储层预测提供可靠的参考依据。获取储层孔隙度的方法主要分为直接测量法和间接解释法两种，直接测量法是通过岩心测量获取孔隙度，后者是通过测井资料解释获得。钻井取心而后进行岩石物理分析所获得的孔隙度最为准确(潘新朋等，2018)，但取样和测试成本高，而且随着钻井深度的增加，获取孔隙度的难度也随之增加，故直接测量法仅限于某些层段，难以获取整个工区的完整孔隙度。间接解释法是据地质信息和测井资料之间的关系确定储层孔隙度，在油气资源勘探中变得越来越重要。但由于影响储层孔隙度的因素众多，且关系复杂，为储层孔隙度的准确预测带来了巨大困难。前人通过建立经验公式或简化地质模型来计算和预测未知层段的储层孔隙度(Khaksar et al.,1998；张显文等，2018)，这些方法对于解决一般地质储层问题能取得较好的结果。但这些方法不能很好地挖掘出孔隙度与测井资料之间复杂的非线性关系及空间的连续性。因此，探索发展高效、低成本的孔隙度预测方法，以实现精确储层预测与评价具有重要的意义。

通常，测井数据包含丰富的物性参数信息，如何对它们进行求解是地球物理勘探具有重要战略意义的研究方向。经研究表明，地下弹性参数之间存在很强的相关性。虽然测井数据较少但是其所包含的物性参数却是非常丰富的，在这些参数具有强相关性的前提下如何提高求解的参数的数量是一个非常有前景的研究方向。

地质统计学预测孔隙度所采取的的插值方法，需要依赖测井数据，即它的结果与测井数据关系密切。通常在工区钻井所得到的测井数据是比较少的，这时如果通过插值方法进行孔隙度预测的话，对于物理空间距离较远的井之间就会出现大面积孔隙度预测值都相近的情况；另外如果相邻井之间对应的值极性相反的话，就会在插值时候中间出现左右极性反向的情况，这些然是不符合实际情况的；一般的，地下的地质结构通常是复杂多变的，插值方法对于大多数的孔隙度预测情况可以说是都不适用的，它不能够刻画孔隙度的变化趋势，只能说是对与孔隙度的一个宏观大致预测，所得到的结果参考价值不大。

线性拟合的孔隙度预测方法，依赖于其它物性参数的测井数据，通过对其他物性参数数据与孔隙度数据利用线性回归的方法，拟合出一个描述它们之间联系的曲线关系。对于工区地下地质情况单一，各弹性参数联系单一的情况，该方法很好地解决了地质统计学通过插值方式预测孔隙度的问题，具有良好的效果；但是对于工区地下地质情况复杂的时候，该方法只能描述用于线性拟合的井曲线的单一关系，只具有局部单一性，并不能描述整体的关系，所以，该方法具有很大的局限性。

通过神经网络进行储层孔隙度预测的方法，主要分为数据处理、模型构建、模型训练、孔隙度预测、模型评价几个阶段。该方法相比较于上述地质统计学和线性拟合的方法，可以很好地应对复杂地质构造的情况；不同学者在应用神经网络方法对孔隙度进行预测时都使用了自己改进的网络模型，这说明针对于不同的地质情况来说，不同网络的表现性能有好有坏。说明通过神经网络预测的结果只对于特定的结构有效，更换数据后就要重新训练网络，这是因为地下地质情况多变，训练得到强大的可以预测所有地下地质情况的神经网络是不现实的。为了使训练出来的网络准确性高，通常是需要耗费很长时间来进行的。

地震反演问题本身就是一个病态反问题，其带来的后果就是直接求解反问题往往不能得到稳定且唯一的解。充分利用测井数据和地质构造信息作为约束，不仅可以使地震反演问题的解趋于稳定，而且也可以从单参数或部分参数求解实现多参数求解，使得反演结果更加符合实际地质情况。

受勘探开发成本的影响，往往只能采集得到少量的测井数据，并且现有的地震反演求解方法只是使用并求解单一或部分参数。虽然字典反演方法对于实际的复杂地质构造可以有很好的自适应性，但是现有的技术只能够计算部分弹性参数。

发明内容

本发明的目的在于克服现有技术的不足，提供一种从测井数据出发，引入字典学习和稀疏表示理论，实现了孔隙度字典反演方法的构建，能够反演求解更多的测井物性参数，能够降低反演成本和反演周期，提高反演过程中测井数据的利用率的基于联合字典的孔隙度反演方法。

本发明的目的是通过以下技术方案来实现的：一种基于联合字典的孔隙度反演方法，包括以下步骤：

S1、数据预处理：将测井数据进行分块处理，去除地震数据中的噪声信号；

S2、将分块处理后的测井数据进行字典学习，得到阻抗字典和孔隙度字典；

S3、将地震数据进行基于字典的孔隙度反演，得到地震数据的阻抗和孔隙度。

进一步地，测井数据进行分块处理的具体实现方法为：将测井曲线分成N个相同长度的数据块，每个数据块的长度为L，得到训练样本集

Y_Ip表示阻抗数据块集合，作为阻抗字典D_Ip学习的训练样本集；Y_por表示孔隙度数据块集合，作为孔隙度字典D_por学习的训练样本集。

进一步地，所述步骤S2具体实现方法为：将字典的学习过程描述为下述问题的优化过程：

其中，矩阵

为期望学习得到的稀疏字典；矩阵A＝[α₁ α₂ …α_N]是稀疏表示矩阵，向量α_i为阻抗字典D_Ip和孔隙度字典D_por所共有的稀疏系数；λ是正则化参数；

表示的是向量的F范数，即Frobenius范数，该式表示矩阵元素绝对值的平方和；||·||₀表示的是向量的0范数，表示向量中非零元素的个数；K表示的是预设的稀疏度最大值，即非零元素的个数；

具体字典学习的过程为：

S21、设置正则化参数λ、最大迭代次数T，初始化字典D和稀疏系数矩阵A，当前迭代次数置为t＝1；

S22、进行稀疏编码：假设字典D固定，使用匹配正交追踪算法对训练样本集Y进行稀疏编码，得到稀疏系数矩阵A；

S23、进行字典更新：采用逐个更新策略对字典D中的原子进行更新：k＝1,2,3,…,K，对于第k个原子，进行如下步骤：

S231、

为当前要更新的原子，记I_k＝{i|α_i(k)≠0,1≤i≤N}，α_i(k)为α_i中的第k个元素，则I_k表示所有样本中使用到第k列原子d_k的索引集；

S232、用

表示去除字典中第k个原子后的字典，

表示去除第k行系数以后的系数矩阵，然后计算新的字典与系数矩阵与原信号之间的残差矩阵

E^k表示去除第k个原子后对所有样本造成的误差；

S233、根据I_k中的索引值选取E_k中相应的列向量，构成新的误差矩阵

并对

进行奇异值分解，得到

S234、选取矩阵U中的第一列，即最大特征值对应的特征向量，作为更新后的原子，并更新相应的系数；

S24、令t＝t+1，如果t≤T，则返回步骤S22；否则停止迭代，输出字典

进一步地，所述步骤S3具体实现方法为：

S31、提取地震子波ω(t)；

S32、进行井震标定；

S33、建立初始模型：反演方法是建立在地震褶积模型的基础上的，反演模型为：

d(t)＝ω(t)*r(t) (2)

其中，d(t)是地震记录，ω(t)是地震子波，*是褶积运算操作符，r(t)是地层反射系数序列；

当地震波沿着地层的垂直方向入射时，地层反射系数与地层波阻抗具有如下关系：

令地震子波序列为ω(t)＝[ω₁ ω₂ …ω_m]^T，反射系数序列为r(t)＝[r₁ r₂ …r_n-1]^T，地震记录为d(t)＝[d₁ d₂ …d_n-1]^T，其中m为子波序列长度，将式(2)改写成如下的矩阵形式：

简写为：

d＝Wr (5)

其中W为子波褶积矩阵，由下式进行定义：

对于叠后波阻抗反演，根据Russell波阻抗近似公式，当波阻抗变化量相对比波阻抗的绝对数值来说非常小时，将式(3)改写为：

根据式(7)，将反射系数序列r表示为如下的矩阵形式：

将式(8)带入式(4)，得到：

令

式(9)表示为：

d＝WDm (10)

令G＝WD，式(10)进一步简写成：

d＝Gm (11)

其中d为(n-1)×1维向量，代表单道地震叠后记录；G为(n-1)×n维矩阵，代表正演模型矩阵；m为n×1维向量，代表单道波阻抗系数的对数；

S34、将

作为正则化约束项加入到反演目标函数中，得到最终目标求解函数：

其中，m_Ip表示阻抗，

表示通过目标函数得到的阻抗最优解，R_i是一个由0和1形成的对阻抗进行取块的矩阵，

表示矩阵元素绝对值的平方和；

S35、通过L曲线法选取正则化参数λ和μ；并设置稀疏度K；

S36、利用最优化相关算法最小化目标函数，求解公式(12)。

本发明的有益效果是：本发明对于现有技术预测孔隙度需要人为预先设定模型参数服从某种特定数学模型以及测井数据较少时出现的预测结果不全面等问题，从测井数据出发，引入字典学习和稀疏表示理论，实现了孔隙度字典反演方法的构建，解决了传统方法对于复杂地质构造求解孔隙度适用性差的缺点，同时能够反演求解更多的测井物性参数。为提高测井数据在反演中的利用率提供了一种可行性方案，降低反演成本和反演周期，提高反演过程中测井数据的利用率。

附图说明

图1为本发明基于联合字典的孔隙度反演方法的流程图；

图2为本实施例的Marmousi阻抗模型图；

图3为实施例1正线性关系的Marmousi孔隙度模型示意图；

图4为实施例1负线性关系的Marmousi孔隙度模型示意图；

图5为实施例1得到的阻抗字典和孔隙度字典；

图6为实施例1正线性拟合孔隙度预测结果；

图7为实施例1负线性拟合孔隙度预测结果；

图8为实施例1中，采用回归方法得到的孔隙度剖面结果图；

图9为实施例1中，采用本发明所提出方法的孔隙度剖面结果；

图10为实施例2得到的阻抗字典和孔隙度字典；

图11为实施例2得到的盲井孔隙度预测结果；

图12为实施例2孔隙度预测结果。

具体实施方式

本发明以字典学习和稀疏表示等理论为基础，考虑测井曲线不同物性参数之间的内在联系，以及反演求解过程中速度、密度、阻抗与孔隙度等之间的联系等因素，发明适用于参数求解的字典反演方法，实现除速度、密度、阻抗外其他测井曲线物性参数(孔隙度)的求解。本发明所使用的联合字典方法，字典训练时间短，且高效快速。本发明的字典原子也可以很好地表示地下沉积微相的构造特征。在这里，就可以不用像线性拟合方法那样事先考虑不同测井曲线之间的关系，通过联合字典的方法，从微观就可学习得到对于不同岩性特征的不同测井曲线之间的关系。而且只需要少量的测井数据，便可充分进行这种关系的学习。

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

如图1所示，本发明的一种基于联合字典的孔隙度反演方法，包括以下步骤：

S1、将测井数据和地震数据分别进行数据预处理：将测井数据进行分块处理；去除地震数据中的噪声信号；地震数据预处理主要是在专业软件上对采集到的地震数据进行动静校正、噪声压制、振幅步长、反褶积、抽道集、速度分析、偏移等。目的是得到能真实反映地层情况的地震数据，处理的目的是去除噪声干扰、增强地震信号等，便于后续反演方法的进行。

测井数据进行分块处理的具体实现方法为：将测井曲线分成N个相同长度的数据块，每个数据块的长度为L，得到训练样本集

Y_Ip表示阻抗数据块集合，作为阻抗字典D_Ip学习的训练样本集；Y_por表示孔隙度数据块集合，作为孔隙度字典D_por学习的训练样本集；L的取值范围可根据地震数据和测井曲线给定。

S2、将分块处理后的测井数据进行字典学习，得到阻抗字典和孔隙度字典；

测井数据物性参数能够精细刻画地下介质的实际情况，充分利用测井数据能够丰富字典学习的样本特征，为字典反演刻画复杂地质构造油气藏相关特征提供了更加完备的信息。常规字典反演方法并没有充分利用测井数据，本发明能够在一定程度上提高实际测井数据的使用率。

将字典的学习过程描述为下述问题的优化过程：

其中，矩阵

为期望学习得到的稀疏字典；矩阵A＝[α₁ α₂ …α_N]是稀疏表示矩阵，向量α_i为阻抗字典D_Ip和孔隙度字典D_por所共有的稀疏系数

λ是正则化参数；

表示的是向量的F范数，即Frobenius范数，该式表示矩阵元素绝对值的平方和；||·||₀表示的是向量的0范数(即稀疏度)，表示向量中非零元素的个数；K表示的是预设的稀疏度最大值，即非零元素的个数；

具体字典学习的过程为：

S21、针对字典学习的目标求解函数

设置正则化参数λ、最大迭代次数T，初始化字典D和稀疏系数矩阵A，当前迭代次数置为t＝1；

S22、进行稀疏编码：假设字典D

固定，使用匹配正交追踪算法对训练样本集Y进行稀疏编码，得到稀疏系数矩阵A；

S23、进行字典更新：采用逐个更新策略对字典D中的原子进行更新：k＝1,2,3,…,K，对于第k个原子，进行如下步骤：

S231、

为当前要更新的原子，记I_k＝{i|α_i(k)≠0,1≤i≤N}，α_i(k)为α_i中的第k个元素，则I_k表示所有样本中使用到第k列原子d_k的索引集；

S232、用

表示去除字典中第k个原子后的字典，

表示去除第k行系数以后的系数矩阵，然后计算新的字典与系数矩阵与原信号之间的残差矩阵

E^k表示去除第k个原子后对所有样本造成的误差；

S233、根据I_k中的索引值选取E_k中相应的列向量，构成新的误差矩阵

并对

进行奇异值分解，得到

假设

是m×n阶矩阵，则其中U为m×m酉矩阵，Λ为半正定m×n对角矩阵，V^T是n×n阶酉矩阵；

S234、选取矩阵U中的第一列，即最大特征值对应的特征向量，作为更新后的原子，并更新相应的系数；

S24、令t＝t+1，如果t≤T，则返回步骤S22；否则停止迭代，输出字典

S3、将地震数据进行基于字典的孔隙度反演，得到地震数据的阻抗和孔隙度；

具体实现方法为：

S31、提取地震子波ω(t)，可以采用确定性子波提取方法或统计性子波提取方法等常用的方法进行提取；

S32、进行井震标定，联系地震数据和测井数据；

S33、建立初始模型：根据地震层位解释信息，按照一定地质沉积规律在解释层位间细分出多个小层，搭建出相对精细的地质框架结构。然后基于该框架，采用数学插值方法，在单井震标定的基础上，将测井数据沿地震解释层位进行内插和外推，构建初始模型。

地震反演，是相对于地震正演来说的。本发明所提出的反演方法是建立在Robinson(1985)地震褶积模型(正演过程)的基础上的，反演模型为：

d(t)＝ω(t)*r(t) (2)

其中，d(t)是地震记录，ω(t)是地震子波，*是褶积运算操作符，r(t)是地层反射系数序列；

当地震波沿着地层的垂直方向入射时，地层反射系数r_l与地层波阻抗Z_l+1、Z_l具有如下关系：

令地震子波序列为ω(t)＝[ω₁ ω₂…ω_m]^T，反射系数序列为r(t)＝[r₁ r₂…r_n-1]^T，地震记录为d(t)＝[d₁ d₂ …d_n-1]^T，其中m为子波序列长度，将式(2)改写成如下的矩阵形式：

简写为：

d＝Wr (5)

其中W为子波褶积矩阵，由下式进行定义：

对于叠后波阻抗反演，根据Russell波阻抗近似公式，当波阻抗变化量相对比波阻抗的绝对数值来说非常小时，将式(3)改写为：

根据式(7)，将反射系数序列r表示为如下的矩阵形式：

将式(8)带入式(4)，得到：

令

式(9)表示为：

d＝WDm (10)

令G＝WD，式(10)进一步简写成：

d＝Gm (11)

其中d为(n-1)×1维向量，代表单道地震叠后记录；G为(n-1)×n维矩阵，代表正演模型矩阵；m为n×1维向量，代表单道波阻抗系数的对数；

S34、根据地质学相关理论可知，地下的物理属性具有一定的相似性和连续性，可以近似认为一定范围内的地下地层参数具有相同的稀疏表示基，这些基可以通过字典学习的方式从测井数据中获得。本发明在反演过程中进行改进，对于字典学习的稀疏表示正则化项进行改进。将步骤S2得到的字典稀疏表示

作为正则化约束项加入到反演目标函数中，得到最终目标求解函数：

其中，m_Ip表示阻抗，

表示通过目标函数得到的阻抗最优解，R_i是一个由0和1形成的对阻抗进行取块的矩阵，

表示矩阵元素绝对值的平方和；

S35、通过L曲线法选取正则化参数λ和μ；并设置稀疏度K；

S36、利用最优化相关算法最小化目标函数，求解公式(12)。

实施例1

将本发明的方法应用于Marmousi模型数据来验证方法的有效性。首先根据图2所示的Marmousi阻抗模型，利用数理统计中的回归分析方法生成两个Marmousi孔隙度模型：一个是由Marmousi阻抗模型正相关生成的孔隙度模型，如图3所示。一个是由Marmousi阻抗模型负相关生成的孔隙度模型，如图4所示。图中，纵坐标(Time)表示时间，横坐标(Tracenumbers)表示道集。并分别为两个模型添加30dB的高斯白噪声。

对于Marmousi波阻抗模型，本发明用主频为35Hz的雷克子波，使用褶积模型生成理论地震数据。与其他方法不同，本发明的方法需要学习一种不同参数的联合字典

因此，本发明方法从模型波阻抗和孔隙度数据中分别等间隔的抽取10道井数据作为训练数据。此次的训练大小为20×2000，即每个字典有2000个原子，每个原子的采样点数为20。对抽取的阻抗和孔隙度数据分别进行原子大小为20步长为1的滑块处理，得到联合字典训练集，通过联合字典学习的方式得到阻抗与孔隙度的联合字典集。图5为本实施例训练得到的字典集(部分)，其中，a)Ip_Dic为Marmousi阻抗模型学习得到的字典原子；b)Por_Dic为由Marmousi阻抗模型负相关生成的孔隙度模型学习得到的字典原子c)；Por_dic为由Marmousi阻抗模型正相关生成的孔隙度模型学习得到的字典原子。

可以很清楚的观察出不同拟合关系的阻抗和孔隙度所学习到的字典原子呈现出的不同关系。表1展示了图5中对应原子的线性拟合系数，可以看出他们之间的线性关系的微弱变化也可以表示出来。而传统方法的线性拟合方法就只能表示一种关系。说明了通过联合字典学习的到的不同物性参数之间的对应字典原子可以表示这种微相关系。

表1.图5中不同字典原子的拟合系数

接下来，我们致力于对式(12)的求解。经过以下的求解步骤：

(1)通过测井阻抗数据建立初始模型，给定最大迭代次数K和允许误差精度Errorlimit；

(2)给定阻抗的初值m₀；

(3)通过传统的正则化方法生成新的m_k；

(4)计算稀疏系数，进行稀疏字典重构得到m_k+1；

(5)判断相邻两次结果的误差值是否超出Error limit或者迭代次数超过K；

(5)不满足则跳到(3)，满足则停止。

本发明需要专注于那些参与训练的道。以道#43为例。图6展示了正线性拟合关系孔隙度使用线性回归方法和本发明的孔隙度反演方法的孔隙度预测输出结果，a)为线性回归方法预测结果；b)为本发明方法预测结果；图中，纵坐标(Time)表示时间，横坐标(Porosity)表示孔隙度。图7展示了负线性拟合关系孔隙度使用线性回归方法和本发明的孔隙度反演方法的孔隙度预测输出结果，a)为线性回归方法预测结果；b)为本发明方法预测结果。从图6a)和图7a)中可以看出圈出的位置线性回归方法并不能很好的预测出来。从图3和图4的孔隙度模型中我们可以看到，模型的上部构造较为单一且地层变化不明显，而模型的下部构造则较为复杂。对于模型中构造较为简单的位置线性拟合方法可以进行较为准确的预测。对于模型中构造较为复杂的位置，也就是图6a)中圈出的位置，线性拟合方法则不能够进行孔隙度的准确预测。拟合方法的结果在图中圈出的位置出现了异常值，这应该是由于拟合方法单一关系并不能准确描述数据整体趋势变化所导致的。而本发明所提出的方法能够将曲线不同位置的响应特征很好地刻画出来，从而得到与实际测井数据吻合度更高的结果。在字典反演领域并没有学者应用字典学习和稀疏表示方法来反演孔隙度。本发明方法可以吻合地下地质体的孔隙度发育趋势。

接下来不只是专注于某一道的预测结果，而是专注于整个剖面上的所有道，并对于反演结果进行研究。这个部位只进行负线性拟合孔隙度的求解。图8展示了线性拟合预测孔隙度剖面结果，图9展示了联合字典方法预测孔隙度剖面结果，图中，纵坐标(Time)表示时间，横坐标(Trace numbers)表示道集。通过对比结果图8和原始孔隙度剖面图4，可以很清楚的看到通过回归方法预测得到的孔隙度降低了分辨率。通过对比图8和图9，可以看到本文所提出的方法分辨率与原始模型数据更加吻合。由于本测试的数据物性关系比较单一，所以两种方法的预测结果都叫好，但是通过单井测试结合剖面说明本文所提出的方法与原始的真实数据吻合度更高，特征响应更加明显。

实施例2

为了进一步说明本发明方法的实际应用能力，将本发明方法应用于中国川中某气田地区工区。对测井的抗和孔隙度数据分别进行原子大小为20步长为1的滑块处理，得到联合字典训练集，通过联合字典学习的方式得到阻抗与孔隙度的联合字典集。图10为本实施例训练得到的字典集(部分)，其中，a)Ip_Dic为学习得到的阻抗字典、b)Por_Dic为学习得到的孔隙度字典。可以很清楚的观察出阻抗和孔隙度所学习到的字典原子呈现出的不同关系。表2展示了图10中对应原子的线性拟合系数。

表2图5中不同字典原子的拟合系数

线性拟合的拟合系数通过多井交汇的方法对阻抗数据和孔隙度数据进行交会图分析，从而确定线性拟合的系数。我们仍然专注于某一道的结果，以没有参与字典训练的井(盲井)为例，图11为盲井孔隙度预测结果，其中，a)为线性拟合结果，b)为本发明方法反演结果。图中，纵坐标(Time)表示时间，横坐标(Porosity)表示孔隙度。该地区下方阻抗数据与孔隙度数据关系较为复杂，可以看到回归方法的结果已经不能够准确描述这种非线性关系，而本文所提出的方法通过联合字典的学习方法能够将他们的非线性特征较为准确的描述出来，从而得到与实际测井数据吻合度更高的结果。

接下来不只是专注于某一道的预测结果，而是专注于整个剖面上的所有道，并对于反演结果进行研究。得到本实施例的孔隙度预测结果如图12所示，(a)为线性拟合结果，(b)为本发明方法反演结果；图中，纵坐标(Time)表示时间，横坐标(Trace numbers)表示道集。两个图中的黑色曲线为层位线。

根据实际资料显示，从古地理学角度对该地区地质时期自然环境的形成和发展演变进行追溯，结合该地质时期的地质构造进行分析。在地壳运动变化过程中，岩石相对高部位抬升出露地表，受到了长期的风化剥蚀作用，形成了较为良好的储集空间。经过其它地方岩石的搬运、沉积、压实覆盖在储层上方。可知该地区目的层上方孔隙度较低，为高阻对低孔，而目的层位下方则大致呈现出高阻对高孔的情况。地质统计学方法的孔隙度预测结果只能部分表现出这个对应关系，而其该方法很容易得到与实际地质情况不符的结果；看图12(a)所示的线性拟合结果，在目的层下方的结果确实比较准确的刻画了阻抗与孔隙度的对应关系，但是由于线性拟合方法的关系单一特性，目的层上面的孔隙度预测结果出现了与实际地质情况相反的情况；最后，图12(b)所展示的是本发明所使用的联合字典方法的孔隙度结果，可以看到，本发明的方法所得到的孔隙度结果也更加符合实际的地质规律，能为后面进一步的油气勘探提供更好的指导。

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

Claims (4)

1.一种基于联合字典的孔隙度反演方法，其特征在于，包括以下步骤：

S1、数据预处理：将测井数据进行分块处理，去除地震数据中的噪声信号；

S2、将分块处理后的测井数据进行字典学习，得到阻抗字典和孔隙度字典；

S3、将地震数据进行基于字典的孔隙度反演，得到地震数据的阻抗和孔隙度。

2.根据权利要求1所述的一种基于联合字典的孔隙度反演方法，其特征在于，测井数据进行分块处理的具体实现方法为：将测井曲线分成N个相同长度的数据块，每个数据块的长度为L，得到训练样本集

Y_Ip表示阻抗数据块集合，作为阻抗字典D_Ip学习的训练样本集；Y_por表示孔隙度数据块集合，作为孔隙度字典D_por学习的训练样本集。

3.根据权利要求2所述的一种基于联合字典的孔隙度反演方法，其特征在于，所述步骤S2具体实现方法为：将字典的学习过程描述为下述问题的优化过程：

其中，矩阵

为期望学习得到的稀疏字典；矩阵A＝[α₁ α₂ … α_N]是稀疏表示矩阵，向量α_i为阻抗字典D_Ip和孔隙度字典D_por所共有的稀疏系数；λ是正则化参数；

表示的是向量的F范数，即Frobenius范数，该式表示矩阵元素绝对值的平方和；||·||₀表示的是向量的0范数，K表示的是预设的稀疏度最大值；

具体字典学习的过程为：

S21、设置正则化参数λ、最大迭代次数T，初始化字典D和稀疏系数矩阵A，当前迭代次数置为t＝1；

S22、进行稀疏编码：假设字典D固定，使用匹配正交追踪算法对训练样本集Y进行稀疏编码，得到稀疏系数矩阵A；

S23、进行字典更新：采用逐个更新策略对字典D中的原子进行更新：k＝1,2,3,…,K，对于第k个原子，进行如下步骤：

S231、

为当前要更新的原子，记I_k＝{i|α_i(k)≠0,1≤i≤N}，α_i(k)为α_i中的第k个元素，则I_k表示所有样本中使用到第k列原子d_k的索引集；

S232、用

表示去除字典中第k个原子后的字典，

表示去除第k行系数以后的系数矩阵，然后计算新的字典与系数矩阵与原信号之间的残差矩阵

E^k表示去除第k个原子后对所有样本造成的误差；

S233、根据I_k中的索引值选取E_k中相应的列向量，构成新的误差矩阵

，并对

进行奇异值分解，得到

S234、选取矩阵U中的第一列，即最大特征值对应的特征向量，作为更新后的原子，并更新相应的系数；

S24、令t＝t+1，如果t≤T，则返回步骤S22；否则停止迭代，输出字典

4.根据权利要求3所述的一种基于联合字典的孔隙度反演方法，其特征在于，所述步骤S3具体实现方法为：

S31、提取地震子波ω(t)；

S32、进行井震标定；

S33、建立初始模型：反演方法是建立在地震褶积模型的基础上的，反演模型为：

d(t)＝ω(t)*r(t) (2)

其中，d(t)是地震记录，ω(t)是地震子波，*是褶积运算操作符，r(t)是地层反射系数序列；

当地震波沿着地层的垂直方向入射时，地层反射系数与地层波阻抗具有如下关系：

令地震子波序列为ω(t)＝[ω₁ ω₂ … ω_m]^T，反射系数序列为r(t)＝[r₁ r₂ … r_n-1]^T，地震记录为d(t)＝[d₁ d₂ … d_n-1]^T，其中m为子波序列长度，将式(2)改写成如下的矩阵形式：

简写为：

d＝Wr (5)

其中W为子波褶积矩阵，由下式进行定义：

对于叠后波阻抗反演，根据Russell波阻抗近似公式，当波阻抗变化量相对比波阻抗的绝对数值来说非常小时，将式(3)改写为：

根据式(7)，将反射系数序列r表示为如下的矩阵形式：

将式(8)带入式(4)，得到：

令

式(9)表示为：

d＝WDm (10)

令G＝WD，式(10)进一步简写成：

d＝Gm (11)

其中d为(n-1)×1维向量，代表单道地震叠后记录；G为(n-1)×n维矩阵，代表正演模型矩阵；m为n×1维向量，代表单道波阻抗系数的对数；

S34、将

作为正则化约束项加入到反演目标函数中，得到最终目标求解函数：

其中，m_Ip表示阻抗，

表示通过目标函数得到的阻抗最优解，R_i是一个由0和1形成的对阻抗进行取块的矩阵，

表示矩阵元素绝对值的平方和；

S35、通过L曲线法选取正则化参数λ和μ；并设置稀疏度K；

S36、利用最优化相关算法最小化目标函数，求解公式(12)。

Images