技术领域
[0001] 本发明涉及岩石物理参数智能预测技术领域，特别地涉及一种基于结构张量和联合字典的物性参数预测方法。
相关背景技术
[0002] 储层物性参数，如孔隙度，含水饱和度，泥质含量等，是直接进行地下油气藏定位的重要参数。地震勘探得到的地震数据中蕴含着丰富的地震弹性参数和储层物性参数信息。一般储层物性参数的预测方式分为两步进行：先通过地震反演从地震数据中获取弹性参数，然后将弹性参数进行岩石物理反演得到物性参数。
[0003] 通常先进行地震反演再进行岩石物理反演的物性参数方法，通常会存在误差累积的问题，即地震反演的弹性参数结果可能就存在误差，这一误差结果又被代入岩石物理反演中进行，从而造成误差进一步增大。同时由于实际的地下地层的情况复杂，在进行岩石物理反演之前通常需要选择一个合适的岩石物理模型，并设置合适的参数，且最终也不能保证岩石物理模型是最优的，从而导致方法的自适应性差，且由于岩石物理模型的强非线性，通常进行岩石物理反演时，需要选择一些概率统计学方法，或神经网络等方法，导致岩石物理反演的计算量大。
具体实施方式
[0068] 下面将结合附图对本发明作进一步说明。本发明的完整流程图如图2所示。
[0069] 本发明提供了一种基于结构张量和联合字典的物性参数预测方法，包括以下步骤：
[0070] S1、建立叠前地震正演模型，叠前地震正演模型基于以波阻抗反射系数为变量的近似方程以及地震褶积模型所构建；
[0071] S2、通过已知地震数据，对叠前地震正演模型反演得到弹性参数，使用联合稀疏表征技术，从测井曲线出发，使用联合字典方法提取得到测井曲线的特征，同时学习得到地下弹性参数与物性参数间的关系；将学习得到的字典划分为弹性参数字典和物性参数字典后，使用弹性参数字典约束叠前地震反演，同时将地震反演过程中得到的稀疏系数与物性参数字典相乘；
[0072] S3、求解得到物性参数后，对物性参数进行张量扩散。
[0073] 在一个实施例中，对二维图像I中的每一个像素点I(p)(p＝[x,y]代表像素点的坐标)，定义其梯度为：
[0074]
[0075] 则像素点p的初始结构张量T0为一个2×2的对称矩阵：
[0076]
[0077] 为了抑制图像中可能存在的噪声，并得到更清晰的图像边缘信息和梯度方向信息，
[0078] 通常会对初始结构张量场T0(p)进行平滑滤波处理：
[0079]
[0080] 其中 表示滤波器。若采用不同的滤波器(如线性的或非线性的)，便可以得到不同类型的结构张量。
[0081] 对于一个给定的图像，各像素点处的梯度向量代表了该像素点处的结构强度和方向信息，但不能充分地体现图像的局部结构特性。通过结构张量T(p)则可获得更丰富的局部结构信息，并且能够度量像素点p的领域内的各向异性结构特征。对T(p)进行特征值分解，有：
[0082] T(p)＝λuuuT+λvvvT
[0083] 其中，λu和λv为T(p)的特征值且λu≥v；u和v为对应的长度为2的特征向量，且u⊥v。针对结构张量 特征值和特征向量的计算公式分别为：
[0084]
[0085]
[0086] 和
[0087]
[0088] v⊥u
[0089] 其中，u对应图像边缘的法向，v对应图像边缘的切向。
[0090] 图1(a)‑图1(c)从几何角度直观的解释了不同情况下的结构张量所代表的图像局部结构特征。总的来说，可以分为如下三种情况：
[0091] (1)λu≈λv≈0，如图1(a)，表示图像在该点附近沿任何方向的变化都很小，对应图像平坦区域的特征。
[0092] (2)λu＞＞λv≈0，如图1(b)，表示图像某一方向的变化率远大于垂直于此方向的变化率，对应图像的边缘结构特征或流线状(Flow‑Like)结构特征。
[0093] (3)λu≈λv＞＞0，如图1(c)，表示图像在两个相互垂直的方向上均有较大的变化，对应图像拐角区域的特征。
[0094] 在一个实施例中，叠前地震正演模型主要是基于以波阻抗反射系数为变量的近似方程以及地震褶积模型所构建。近似方程形式如下所示：
[0095] RPP(θ)＝c1RP+c2RS+c3RD2 2 2
[0096] 其中 RPP(θ)为随角度变化的PP波反射系数；c1＝1+tanθ，c2＝‑8γtan2 2 2
θ，γ＝VS/VP；c3＝‑0.5tanθ+2γsinθ；
[0097] RP、RS、RD分别为纵波阻抗反射系数、横波阻抗反射系数、横波阻抗反射系数以及密度反射系数，他们与纵横波速度及密度有如下关系：
[0098]
[0099]
[0100]
[0101] 其中ΔVP、ΔVS和Δρ分别对应分界面两侧的纵波速度、横波速度和密度的差；Vp、VS和ρ为两侧对应属性的平均值。
[0102] 在一个实施例中，地震波在地下传播产生的反射波和透射波的过程通常可以认为是离散的过程，也即为地震正演的过程，从而我们可以将地下地层分为Nl个的分层，相邻的分层之间即为一个分界面，令i＝1,2,…,Nl‑1表示层i与层i+1的分界面。通常相邻的两个地层之间的物理性质差别不大，因此有相关研究认为有如下近似：
[0103]
[0104]
[0105] 其中 且 为第i层的纵波阻抗。针对Nl个采样点，上式可以写成如下矩阵形式：
[0106]
[0107] 上式可以简写为：
[0108]
[0109] 其中， 为纵波阻抗反射系数向量；
[0110] 表示纵波阻抗向量的自然对数；
[0111] 为一阶 差分 矩阵。同理 ，令 与分别表示横波阻抗和密度的自然对数；
[0112]
[0113] RD＝TLD
[0114] 其中， 和 分别为横波阻抗反射系数向量和横波阻抗向量的自然对数；
[0115] 和 分别表示密度反射系数向量和密度向量的自然对数。
[0116]
[0117] 假设地层岩石的纵横波阻抗的对数LP、LS以及密度的对数LD三者之间存在一种弱线性关系：
[0118] LS＝kLP+kc+ΔLS
[0119] LD＝mLP+mc+ΔLD
[0120] 其中k、kc、mc以及mc为所拟合出的线性关系的系数常量。可得：
[0121]
[0122] 其中 且 根据褶积模型，从反射系数到地震信号的正演模型即可表示为如下公式：
[0123] d(θ)＝W(θ)RPP(θ)
[0124] 其中 就是入射角度为θ时的合成地震信号， 表示长为Nd‑Nl+2(Nd＞Nl)的子波构造而成的卷积矩阵。令地震信号的入射角的个数为J，可将所有入射角对应的方程纳入到一个方程中，即得到如下叠前三参数正演模型：
[0125]
[0126] 由以上即得到了叠前正演模型，用于反演Vp、Vs、ρ三参数，简记为：
[0127] d＝Gm
[0128] 其中d为同一角道集的J个不同入射角的地震信号列向量话后组成的地震信号，G即为由子波、纵横波背景速度以及拟合系数构建而成的正演算子，m为模型参数由纵横波阻抗和密度组建得到。
[0129] 在一个实施例中，通过构建得到的正演模型，已知地震数据，我们通过反演得到弹性参数，然而反演求解常常为一个超定方程的求解问题，且存在有噪声的影响，因此从确定性角度出发，我们选取l2范数作为适配函数，使用最小二乘算法，构建下面的反演目标函数：
[0130]
[0131] 通常在高维地震信号反演过程中，存在病态性的问题，因此为了降低求解时的多解性，同时为了使得m的求解更加的稳定，且满足某些已知的条件，对于方程中m的估计，通常还要使用正则化技术，即在上式中再加入正则化项：
[0132]
[0133] 其中λ为正则化系数，用于控制正则化项的约束强度，Φ(m)即为正则化项，通过这一正则化项，我们即可加入一些先验信息，从而使得所估计的m符合某些特征，降低反演求解的多解性。然而通常大多数传统方法的正则项是简单的，经验化的，且不适应的，例如最简单的1阶光滑Tikhonov正则化方法便令该函数为(1th‑Tikh)，通过这一正则化项即可使得反演结果出现光滑的特性，然而实际地下地质情况往往是十分复杂的，因此这些方法并不适用具有复杂地质结构的境况。
[0134] 通常地下地层的参数之间存在有一定的横向连续性和一致性，并且地下弹性参数和物性参数之间存在很强的相关性，因此我们选择使用联合稀疏表征技术，从测井曲线出发，使用联合字典方法提取得到他们的特征，同时学习得到地下弹性参数与物性参数间的关系。将学习得到的字典划分为弹性参数字典和物性参数字典后，使用弹性参数字典约束叠前地震反演，从而提高地震反演的精度，同时将地震反演过程中得到的稀疏系数与物性参数字典相乘，得到物性参数，从而实现弹性参数与物性参数同步预测。
[0135] 从测井数据中学习得到弹性参数和物性参数的联合字典，首先需要构建他们的联合样本。由于测井数据的数量是很有限的，因此为了能够更加充分的利用测井数据，我们需要将测井数据通过滑动时窗进行分块的处理来构建联合字典的训练样本，具体操作如图3所示：
[0136] 经过以上的分块操作之后即可得到字典的训练样本Y＝[Y1,Y2,…,Yn]，为滑动时窗得到的分块个数，每一个训练样本Yi＝[VPi,VSi,ρi,φi,Swi,Ci]。同时由于不同的弹性参数、物性参数具有不同的尺度和量级，因此我们采用了最大最小值归一化方法，将弹性参数与物性参数放在同一个量级从而可以训练得到一个联合字典。
[0137] 字典学习的过程便是学习得到一个字典，通过选取这个字典中的少量原子组合即可重构训练样本，如果记此时稀疏系数集为A＝[α1,α2,…,αn]，则这一过程就可以表示为如下能量最小化的问题：
[0138]
[0139]
[0140] 其中，每一个长度为L的列向量αi为yi在字典D下的稀疏表示,||·||0是l0范数，用以计算某一个向量的含零元素的个数，K是表示稀疏程度的常量。我们通常通过K‑SVD算法求解上式，即可学习得到弹性参数和物性参数的联合字典 其中 为弹性参数的字典， 为物性参数的字典。训练得到联合字典后我们将其划分为弹性参数字典和物性参数字典，即可将其与叠前地震反演目标函数结合，构建一个叠前弹性参数与物性参数同步反演的目标函数：
[0141]
[0142] 并满足有下面的约束：
[0143]
[0144] 其中d为地震数据，G为叠前正演算子，m为模型参数， 为归一化后的弹性参i数， 为物性参数的集合，R 为分块矩阵，DE为弹性参数的字典，DP为物性参数的字i
典，α是由联合字典得到的联合稀疏系数，ψ(·)为模型参数与弹性参数之间的转换公式，‑1
ζ (·)为反归一化算子。
[0145] 然而，求解叠前弹性参数与物性参数同步反演的目标函数依然是极其困难的，主要是两方面的问题。第一，需要考虑的变量太多了；第二，由于ψ是一个包含了自然对数计算的元素级的算子，等式约束显得极其复杂。对于前者而言，一个可选方案是，采用坐标下降策略每次只更新一个变量，该策略也被广泛应用到其他问题中。因此，叠前弹性参数与物性参数同步反演的目标函数可以被分解为如下三个不同的子问题：
[0146] (1)模型参数更新问题
[0147] 当只考虑模型参数m为变量时，叠前弹性参数与物性参数同步反演的目标函数可简化为：
[0148]
[0149] 其中x即为VP、VS、ρ。对于约束条件，我们可以令 为惩罚项，从而将上式变为下面的问题：
[0150]
[0151] 可以发现上式中模型参数的更新会同时受到地震正演模型的约束，同时由于 主要是由上一轮迭代中稀疏重构结果而来，因此也会受到字典学习后稀疏表征的重构结果的约束。求解上式，可以通过共轭梯度或者拟牛顿共轭梯度优化算法进行求解。
[0152] (2)稀疏编码问题
[0153] 当只考虑稀疏表征系数α为变量时，叠前弹性参数与物性参数同步反演的目标函数可简化为多个稀疏编码问题：
[0154]
[0155] s.t.||αi||0≤K
[0156] 其中i＝(1,2,...,n)，这一问题直接通过OMP算法求解即可。同时在这一步中，根据对弹性参数进行稀疏编码得到的αi，我们与物性参数的字典DP相乘，即得到了物性参数的预测结果。
[0157] (3)弹性参数更新问题
[0158] 当只考虑弹性参数变量时，叠前弹性参数与物性参数同步反演的目标函数可写为下面的形式：
[0159]
[0160]
[0161] 由于模型参数为固定的，因此我们可以令弹性参数三分量分别为 和 它们‑1通过ψ (m)求逆得来，因此我们可以将上式转化为下面的三个子问题：
[0162]
[0163]
[0164]
[0165] 为了方便求解上述方程，本节通过添加惩罚项的方式，并将它们设置一个相同的权重，通过拉格朗日乘子法将上式近似地改写为如下无约束问题：
[0166]
[0167]
[0168]
[0169] 同理，上式所描述的弹性参数更新问题也具有解析解，因此根据上式，可以得到下面的弹性参数x的三个分量的更新公式：
[0170]
[0171]
[0172]
[0173] 其中：
[0174]
[0175] 和分别为第t次迭代中第i个小块对应的纵波速度、横波速度和密度的稀疏编码。γ系数主要用于控制稀疏表征结果的影响程度，当γ由零至无穷大，表示反演所受字典稀疏表征约束的影响越强。为了方便确定γ系数，我们将上式写为下面的形式：
[0176]
[0177]
[0178]
[0179] 其中三个分量 和 分别对应所有稀疏重构小块 和平均的结果，此时γ的取值范围为[0,1]。上式可化简为
[0180]
[0181] 其中
[0182] 在一个实施例中，3D地震数据的空间结构特征一定程度上代表了地下介质参数的空间结构特征，因此我们使用空间结构张量技术从叠后地震数据中学习得到地下的空间结构特征。
[0183] 对3D地震数据体中的每一个点I(p)(p＝[x,y,z])，定义其梯度为：
[0184]
[0185] 则点p的初始结构张量T0为一个3×3的对称正定矩阵：
[0186]
[0187] 其中，Ix,Iy,Iz则表示3维数据体在点p上关于x,y,z方向的梯度向量。同时为了得到更鲁棒的结构张量场，我们还在上式的基础上，加入了高斯滤波算子，从而得到新的矩阵如下：
[0188]
[0189] 对上式进行特征值分解，可得
[0190] TG(P)＝λuuuT+λvvvT+λwwwT
[0191] 其中u、v和w为结构张量场的特征向量，两两之间正交，它们对应的特征值为λu≥λv≥λw。特征向量u代表3D数据体内某局部平面的法向，特征向量v和w的方向则平行于该局部平面。
[0192] 下面我们给出了某工区地震数据体的一个沿层切片，俯视观察特征向量与数据空间结构的关系。如(图4(a))所示，向量场u有许多的小黑点，这是因为它们都是垂直与平面的关系，向量场v(图4(b))和向量场w(图4(c))的小黑点则较少，因为它们均平行于地震数据的切片平面。
[0193] 接下来对输入3D地震数据体某点I(p)和相应的结构张量场T(p)，各向异性扩散过程可以描述为方程
[0194]
[0195] 在上式中，div代表着散度算子，它的初始条件是g(p；0)＝I(p)。对于任何时间t＞0，我们可以将解g(p；0)视为输入数据I(p)的平滑版本。这里，扩散张量场E(p)由结构张量场T(p)通过共享T(p)的特征向量来构造：
[0196] E(p)＝μuuuT+μvvvT+μwwwT
[0197] 构造得到扩散张量场后E(p)后我们即提取到了地下的结构特征，下一步即将其引入我们反演的过程当中。我们将智能联合稀疏表征约束的多参数反演的目标函数分为了三个子问题，在第一个子问题中我们可以采用共轭梯度或者拟牛顿共轭梯度优化算法进行求解，因此第一个子问题的求解更新将会有下面的形式：
[0198] mt,k+1＝mt,k+αt,kpt,k
[0199] 其中，k表示求解子问题时的内循环迭代索引，mt,k表示外循环第t次迭代且内循环第k次迭代时的模型参数，αt,k为搜索步长，pt,k为搜索方向。为了引入模型参数m的空间结构张量约束，我们借鉴张量扩散的方式，将上式改写为如下的形式：
[0200]
[0201] 其中τ是控制扩散速度的权重参数，E为计算好的扩散矩阵，为元素级乘法(即在每一个模型参数点执行矩阵乘法)，通过这一式子即可在三个子问题迭代求解时引入空间结构的约束。同时对于物性参数的预测同样可以通过这一方式引入空间结构约束，即在第二个子问题求解得到物性参数后，同样使用下式进行张量扩散即可：
[0202]
[0203] 在一个实施例中，首先将本发明应用在Marmousi模型数据来验证方法的有效性。先根据Marmousi纵横波阻抗、密度模型，利用数理统计中的回归分析方法生成一个Marmousi孔隙度模型，并添加30dB的高斯白噪声，生成的孔隙度模型如图5所示：
[0204] 然后从Marmousi纵横波速度模型、密度模型和孔隙度模型中，每隔25道抽取出一个井数据，共15口井的阻抗数据和孔隙度数据作为已知的井数据。设置字典原子长度为40，字典原子个数为1000，测井数据分块的滑动时窗滑动步长则选择为1得到训练集，进行训练得到弹性参数和物性参数的联合字典，最终得到反演结果如下图6(a)‑6(b)所示：
[0205] 由上图的反演结果可以看出基于联合字典的方法和基于结构张量与联合字典的两种方法均能较好的反演还原出原始数据，但在联合字典的方法中，存在着有空间连续性不强，‘挂面条’现象，而在结构张量与联合字典的方法中，则呈现出更好的空间连续性，在局部特征上与原数据更加贴合。为进一步验证联合字典方法有效性，计算各方法的均方误差进行定量分析，逐道计算两种方法与原始数据的均方误差，得到的计算结果如下图7所示：
[0206] RMSE计算结果一般越小越好，由上图我们可以看出基于联合字典的两种方法均方误差都是比较小的，同时可以发现基于结构张量与联合字典的方法的均方误差要略低于仅使用联合字典的方法，从而证明了，引入结构张量后的联合字典反演方法，不仅能提高反演结果的空间连续性，同时能够进一步的提高物性参数预测的精度。
[0207] 将本发明的方法应用到四川盆地某气田地区工区中，以检验方法的有效性。工区中包括有采样得到的叠前地震数据，以及弹性参数(纵横波速度，密度)，物性参数(孔渗饱)的测井数据，共有5口井数据，采样时间为1ms。
[0208] 选择该工区中的三口井进行测试，分别比较联合字典预测方法和基于结构张量和联合字典方法的预测效果。两个方法中联合字典训练的部分，都采用相同的参数，经过测试选择原子长度为40，原子的个数为1600，训练样本选择三口井的纵横波速度、密度、孔隙度、泥质含量、含水饱和度，得到的训练的联合字典。预测结果中孔隙度与泥质含量的反演效果对比对如图8(a)‑8(b)所示。
[0209] 由图9(a)‑9(b)两种方法孔隙度与泥质含量的剖面反演结果可以看出，联合字典进行预测得到的结果，破碎感较强，空间连续性差，而基于结构张量和联合字典的方法的预测结果，则展现出更好的空间连续性。
[0210] 本发明从现有物性参数预测方法的不足出发，从测井数据中学习得到能够代表地下弹性参数和物性参数特征与关系的联合字典，将这一字典通过稀疏表征的方式加入叠前地震反演过程中，同时从叠后地震数据中学习得到地下地层的空间展布规律的信息，通过张量扩散的方式也加入到叠前地震反演过程中，实现了弹性参数与物性参数的同步预测，同时提高了预测结果的空间连续性。
