阅读说明：本技术 一种基于波阻抗低秩正则化的地震波阻抗反演方法 (Seismic wave impedance inversion method based on wave impedance low-rank regularization ) 是由 陈雷平 李曙 贺达江 丁黎明 于 2019-12-02 设计创作，主要内容包括：本发明公开了一种基于波阻抗低秩正则化的地震波阻抗反演方法。本发明直接针对地震波阻抗建立低秩正则化项,并在此基础上建立具有波阻抗低秩约束的目标函数,最后巧妙地利用奇异值分解进行求解,利用波阻抗的低秩性进行地震反演,充分利用了地震波阻抗数据的低秩性先验信息且提高了地震波阻抗反演的准确性。(The invention discloses a seismic wave impedance inversion method based on wave impedance low-rank regularization. The method directly establishes a low-rank regularization item aiming at the seismic wave impedance, establishes a target function with low-rank wave impedance constraint on the basis, finally solves the problem by skillfully utilizing singular value decomposition, performs seismic inversion by utilizing the low-rank property of the wave impedance, fully utilizes the low-rank prior information of seismic wave impedance data and improves the accuracy of seismic wave impedance inversion.)

一种基于波阻抗低秩正则化的地震波阻抗反演方法

技术领域

本发明属于地震勘探领域，尤其涉及一种基于波阻抗低秩正则化的地震波阻抗反演方法。

背景技术

地震波阻抗是一种重要的地球介质参数，它在数值上等于地下岩石的密度与地震波在岩石中传播时的速度的乘积。它与岩石的岩性、孔隙度等紧密相关。因此，如何有效且准确地获得地震波阻抗一直是地震勘探中的重要研究内容。地震波阻抗反演技术根据观测得到的地震、测井等数据利用先验信息反向推演地下介质的波阻抗，是目前获得地震波阻抗的最主要途径。随着地震采集技术和信号处理技术的发展，地震数据的质量越来越高。地震反演结果的好坏越来越明显地受制于先验信息的利用。如何挖掘且有效地使用先验信息成为了当前地震反演领域的热点问题。

数据的低秩性是近些年来被广泛关注和使用的重要先验信息。如：语音和图像信号去噪、地震数据去噪/重构等领域都有大量基于数据低秩性先验的文献报导。在地震勘探领域，低秩性主要用于地震数据去噪、缺失地震道补全、地震数据重构等方面。在这些研究中，所利用的均是地震数据的低秩性。即：地震数据中波形相似的不同子区块的地震数据组合成的矩阵具有低秩的特性。当使用低秩性后，数据的去噪/重构效果比传统的方法有明显的提升。因此，这些方法在实验上证明了低秩性作为一种先验信息在地震数据处理中的应用优势。需要指出的是，这些地震数据处理方法都只是利用了地震数据在数据域的低秩性。它们处理完后的输出结果仍然只是地震数据。

既然地震数据的不同子区块之间具有低秩性，地震数据与地下岩层的性质和结构息息相关的。前者是后者的信号化表现形式。如果前者的不同子区块之间具有相似性，后者的对应位置也应该具有这种相似性。这时，这些子区块的波阻抗、密度、泊松比等波阻抗张成的矩阵就具有低秩性。低秩性是一种非局部化的先验信息。在电磁参数反演中，如ManuelAmaya等人(2016)发表于Geophysics期刊上的论文《A low-rank approximation forlarge-scale 3D controlled-source electromagnetic Gauss-Newton inversion》提到了使用低秩近似的思想。但是，他们的思路是对反演过程中的Hessian矩阵使用其低秩近似予以替代，以降低计算复杂度，并未直接利用要反演的模型参数的低秩性。目前，还未见利用波阻抗的低秩性，根据地震数据，使用地震反演方法得到波阻抗的文献报导。

现有的地震波阻抗反演方法主要是利用了地震数据或波阻抗的局部化先验信息。例如，稀疏反演所利用的稀疏性先验信息即是一种局部化先验信息。它是利用地震数据或波阻抗的局部区域在L0、L1或Lp(0<p<1)范数下的稀疏性，缺乏对不同子区域之间的先验信息的利用。近年来，备受关注的基于全变分正则化的地震波阻抗反演方法利用了波阻抗的梯度稀疏性先验信息，但是这种先验信息也是局部化的，同样没有利用非局部化先验信息。

综合上述，低秩性作为一种重要的非局部先验信息在信号/数据去噪、重构等领域具有突出优势，因而在地震数据处理方面发挥了重要作用。地震波阻抗也应该具有低秩性。然而，现有的地震波阻抗反演方法均未利用这一先验信息。

发明内容

为解决上述问题，本发明公开了一种基于波阻抗低秩正则化的地震波阻抗反演方法。本发明直接针对地震波阻抗建立低秩正则化项，并在此基础上建立具有波阻抗低秩约束的目标函数，最后巧妙地利用奇异值分解进行求解，利用波阻抗的低秩性进行地震反演，充分利用了地震波阻抗数据的低秩性且提高了地震波阻抗反演的准确性。

为实现上述目的，本发明的技术方案为：

一种基于波阻抗低秩正则化的地震波阻抗反演方法，包括如下步骤：

步骤一、输入地震数据、波阻抗初始模型，构建波阻抗的初始低秩近似矩阵：采用基于相似性搜索的方法构建初始低秩近似矩阵；

步骤二、建立波阻抗低秩正则化目标函数

设整个波阻抗剖面可以分成K个相互之间有交叠的子数据块，其中的第i个子数据块记为d_i，则将波阻抗最优低秩近似问题表述为：

Φ_i是一个相似性度量与选择算子，以d_i为参考块，用于从地震波阻抗m中选择出与d_i的欧氏距离最小的Q个子数据块，并排列成一个大小为P×Q的矩阵，P表示子数据块d_i的长度；Ο_i为待求的最优低秩近似矩阵；|| ||_F表示矩阵的Frobenius范数，rank()为求矩阵的秩的运算；λ表示低秩约束的正则化参数；

在式(1)的基础上，构建波阻抗低秩正则化地震波阻抗反演问题的目标函数：

其中，G表示地震波阻抗反演时的正演算子，S表示观测得到的地震数据；η为波阻抗低秩正则项的正则化系数，用于平衡数据保真项

和低秩正则项在反演中所起的作用的大小，值越大说明低秩正则项起的作用越大；

步骤三、目标函数的求解

对于式(2)的目标函数，分两步进行求解，首先求解最优低秩近似，即求解(1)的最优值，式(1)的最优低秩近似写成：

使用奇异值分解和硬阈值收缩求解式(3)，令Φ_im的奇异值分解为：

Φ_im＝Αdiag(ω)Β^H (4)

其中，A和B分别表示奇异值分解后产生的两个酉矩阵；ω是奇异值分解得到的奇异值向量；diag(ω)表示由奇异值分解得到的奇异值向量ω张成的对角矩阵，对角线上的元素为ω；H表示共轭算子；最优低秩近似的解写成：

其中λ表示硬阈值收缩时使用的阈值，HT_λ()表示硬阈值函数，定义式为：

其中，Γ表示输入HT_λ()函数中的数据；

在求得O_i的最优解后，在式(2)目标函数中固定O_i，得到波阻抗反演结果的更新公式为：

式(7)是一个最小二乘问题，将式(7)整理为关于

的方程：

其中，G^*和Φ_i ^*分别为G和Φ_i的伴随算子，使用共轭梯度法快速求解得到反演结果

步骤四、将反演结果

按照与构成长的列向量m相反的过程，以初始波阻抗模型的尺寸，恢复成2维剖面或3维数据体，即完成了整个反演过程。

进一步的改进，η通过交叉验证的方法确定，即：根据反演结果与真实值之间的差异而调整，选择出使得反演误差最小的值即η的值；所述地震数据为由地震检波器记录下来的地震信号。

进一步的改进，所述步骤一中，构建波阻抗的初始低秩近似矩阵步骤如下：先将初始的波阻抗模型按列排列成一个长的列向量m；然后确定一个作为参考的子数据块d，假设子数据块d长度为P；再依次从m中选取与d的欧氏距离最小的Q个子数据块；将Q个子数据块排列成一个矩阵，则完成初始低秩近似矩阵的构建；初始的波阻抗为2维的矩阵或3维的数据体。

附图说明

图1为初始低秩近似矩阵的构建流程图；图中(a)表示波阻抗剖面，(b)表示欧氏距离相近的相似块，(c)初始低秩近似矩阵；

图2为本发明的流程图；

图3为取自Marmousi2模型的部分数据作为波阻抗真实模型；

图4为波阻抗初始模型；

图5为用于反演的地震数据；

图6为主频为40Hz的雷克子波；

图7为反演结果图。

具体实施方式

以下结合附图及实施例对本发明做进一步说明。

实施例1

为了解决现有的地震波阻抗反演方法存在的没有利用波阻抗的非局部先验信息、无法从全局的角度把握地震波阻抗的本质特征等问题，提出了一种基于波阻抗低秩正则化的地震波阻抗反演方法。该方法的核心是：直接针对地震波阻抗建立低秩正则化项，并在此基础上建立具有波阻抗低秩约束的目标函数，最后巧妙地利用奇异值分解进行求解。目前，尚无利用波阻抗的低秩性进行地震反演的研究报导。

(1)构建波阻抗的低秩近似矩阵

采用基于相似性搜索的方法构建低秩近似矩阵。具体实现方法是：先将初始波阻抗模型按列排列成一个长的列向量m。然后，确定一个作为参考的子数据块d，假设其长度为P。再依次从m中选取与d的欧氏距离最小的Q个子数据块。最后，将这Q个子数据块排列成一个矩阵，则得到的将是一个低秩矩阵。构造过程如图1所示。需要说明的是，这里为了便于观察子数据块之间的相似性，没有事先将图1中，(a)中的波阻抗做向量化。

(2)建立波阻抗低秩正则化目标函数

假设整个波阻抗剖面可以分成K个相互之间有交叠的子数据块，其中的第i个子数据块记为d_i，则可以将波阻抗最优低秩近似问题表述为：

Φ_i是一个相似性度量与选择算子，它以d_i为参考块，用于从地震波阻抗m中选择出与d_i的欧氏距离最小的Q个子数据块，并排列成一个大小为P×Q的矩阵。O_i为待求的最优低秩近似矩阵。|| ||_F表示矩阵的Frobenius范数，rank()为求矩阵的秩的运算。λ表示低秩约束的正则化参数。在此基础上，可以构建波阻抗低秩正则化地震波阻抗反演问题的目标函数：

其中，G表示地震波阻抗反演时的正演算子，S表示观测得到的地震数据。η为波阻抗低秩正则项的正则化系数，用于平衡数据保真项

和低秩正则项在反演中所起的作用的大小，其值越大说明低秩正则项起的作用越大。

(3)目标函数的求解

对于上一步构建的目标函数，可以分两步进行求解。第一步是求解最优低秩近似，即求解下式的最优值：

该式的最优化表达式可写成：

使用奇异值分解和硬阈值收缩求解该问题。令Φ_im的奇异值分解为：

Φ_im＝Αdiag(ω)Β^H

其中，A和B分别表示奇异值分解后产生的两个酉矩阵。ω是奇异值分解得到的奇异值向量。diag(ω)表示由奇异值分解得到的奇异值向量ω张成的对角矩阵，对角线上的元素为ω。鉴于此，可以将最优低秩近似

的解写成：

其中λ表示硬阈值收缩时使用的阈值，HT_λ()表示硬阈值函数，其定义式为：

其中，Γ表示输入HT_λ()函数中的数据。

在求得O_i的最优解后，在目标函数中固定O_i，可得到波阻抗反演结果的更新公式为：

这是一个最小二乘问题，可以将其整理为关于

的方程：

其中，G^*和Φ_i ^*分别为G和Φ_i的伴随算子。此问题可以使用共轭梯度法快速求解。

最后，将反演结果按照与构成长的列向量m相反的过程，以初始波阻抗模型的尺寸，恢复成2维剖面或3维数据体，即完成了整个反演过程。

本发明的具体算例如下：取自Marmousi2模型的部分数据作为波阻抗真实模型如图3所示，深度方向有452个采样点，距离方向有512道。波阻抗初始模型如图4所示是对真实模型的高斯低通滤波；用于反演的地震数据如图5所示是用波阻抗真实模型求得反射系数，然后用反射系数与主频为40Hz的雷克子波褶积，再添加20％的高斯白噪声而得到。子波为主频40Hz的雷克子波，如图6所示；结果如图7所示。

可以看出，采用本发明方法得到的波阻抗反演结果与真实的波阻抗模型十分相近，且可以在反演结果中看到地层的分界面刻画得十分清晰。

尽管本发明的实施方案已公开如上，但并不仅仅限于说明书和实施方案中所列运用，它完全可以被适用于各种适合本发明的领域，对于熟悉本领域的人员而言，可容易地实现另外的修改，因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里所示。