CN109471164A - 基于Ho-RPCA的地震断层增强方法 - Google Patents

Abstract

本发明提供一种基于Ho‑RPCA的地震断层增强方法，该基于Ho‑RPCA的地震断层增强方法包括：步骤1，对输入三维相干体数据进行预处理,对低秩的断层部分和稀疏的噪声部分进行分离；步骤2，构建基于Ho‑RPCA的目标函数；步骤3，通过交替方向乘子法算法框架来求解该最优化问题；步骤4，对张量核范数进行最小化，采用基于张量随机奇异值分解的核范数近似求解。该基于Ho‑RPCA的地震断层增强方法从张量模型的角度来对相干体进行增强，更好地保护了相干属性体的高维特性，可以较好地去除噪声，增强后的断层连续性得到提升，能够较好地增强相干体中的断层信息，总体效果比较理想。

Description

基于Ho-RPCA的地震断层增强方法

技术领域

本发明涉及地球物理勘探领域技术领域，特别是涉及到一种基于 Ho-RPCA的地震断层增强方法。

背景技术

地震勘探采集到的地震资料通常是三维或者四维的数据体。地震资料 数据是由多道对应工区地表坐标的地震反射波组成，每道波形都会随着反 射层的地质特性差异而产生不同的波形幅值。地震资料能够为地质分析与 解释提供充足的物理依据。地震数据的纹理属性提取方法是地震资料解释 的重要手段。在众多的纹理属性提取方法中，最具代表性是相干体算法， 它能够通过突出地震资料中的不连续性来识别断层，从而提高断层解释的 效率。由于断层在油气勘探和开发中有着重要的控圈控储控藏作用，因此 发展先进的相干体算法与断层增强技术，实现对断层的准确刻画，一直是 地质领域中研究的热点。

地震资料在信号采集以及资料解释的过程中，会受到高斯噪声、相干 噪声以及其他地质噪声的干扰。现有的纹理属性提取方法仅突出了地震资 料中特定的纹理信息，并不能从根本上消除这些噪声干扰。对原始的振幅 数据体进行相干计算后，得到的相干体信噪比仍较低，反映在相干切片上 断层的横向分辨率低、可辨识度差以及断层的连续性差。对于这样的相干 体，不管是进行断层的人工解释还是自动解释，都具有较高的难度。因此，通常还需要对相干体进行适当的断层增强处理，以削弱各类噪声的影响， 突出断层的线型特征，提升相干体的横向分辨率，以便于后续的构造解释 或是断层自动识别。

地震纹理属性提取方法属于多维图像处理领域，主要研究的是从多维 地震资料中提取出感兴趣的地质信息的方法。在计算机图像学中，纹理指 的是在一个区域中成周期性的均匀变化的结构，在现实世界中，比较典型 的例子是树木纹理，而类似的纹理构造广泛存在于地震资料数据中。断层 提取是地震纹理属性提取方法的一个重要研究方向。断层是地壳中岩层连 续性受到破坏后，沿断裂面发生明显相对移动的一种构造现象，矿床的形成、矿体产状及其分布等，都要受到断层构造的控制；断层的形态和分布 决定了油气田分布，断层提取在油气勘探开发中有着及其重要的作用。研 究人员针对断层识别提取问题发展出了一系列断层纹理属性提取算法，其 中最具代表性的是相干体算法。对原始的振幅数据体进行相干计算后，得 到的相干体信噪比仍较低，反映在相干切片上断层的横向分辨率低、可辨 识度差以及断层的连续性差。对于这样的相干体，不管是进行断层的人工 解释还是自动解释，都具有较高的难度。因此，通常还需要对相干体进行 适当的断层增强处理，以削弱各类噪声的影响，突出断层的线型特征，提 升相干体的横向分辨率，以便于后续的构造解释或是断层自动识别。

对于相干体三维数据而言，其中有用信息是断层部分，在空间上具有 很强的相关性和冗余性，这样的性质可以通过张量的低秩性来描述；而噪 声部分则在三维数据体中呈近似稀疏分布。因此，采用低秩张量恢复模型 (Ho-RPCA)对相干体进行建模，从相干体数据中分离出低秩的断层部分， 以达到断层增强的目的。为此我们发明了一种新的基于Ho-RPCA的地震断 层增强方法，解决了以上技术问题。

发明内容

本发明的目的是提供一种基于Ho-RPCA的地震断层增强方法。

本发明的目的可通过如下技术措施来实现：基于Ho-RPCA的地震断层 增强方法，该基于Ho-RPCA的地震断层增强方法包括：步骤1，对输入三维 相干体数据进行预处理,对低秩的断层部分和稀疏的噪声部分进行分离； 步骤2，构建基于Ho-RPCA的目标函数；步骤3，通过交替方向乘子法算法 框架来求解该最优化问题；步骤4，对张量核范数进行最小化，采用基于 张量随机奇异值分解的核范数近似求解。

本发明的目的还可通过如下技术措施来实现：

在步骤1中，采用鲁棒主成分分析法对输入三维相干体数据进行预处 理，鲁棒主成分分析法的数学模型可表示为：

D＝A+E (1-1)

如果将相干体图像用D来表示，由于断层在相干体数据中是规律的平 面结构，所以它对应的是低秩部分A；在相干体图像中，所有噪声的和是 满足近似稀疏性的，即噪声的和可以对应鲁棒主成分分析法模型中的稀疏 部分E，这样，就将三维相干体数据D用低秩部分A和稀疏部分E表示出 来。

在步骤2中，采用低秩张量恢复模型Ho-RPCA来对三维相干数据体进行 断层增强；

根据Ho-RPCA对相干体三维图像建模，类似于鲁棒主成分分析法的图 像去噪模型，有如下公式：

其中，rank(A)表示对矩阵A取秩；||E||₀是E的0范数，表示矩阵E中非 零元素的个数；参数λ是用来平衡低秩项和稀疏项的权重参数，一般λ>0； D是待去噪增强的相干体数据，A对应了低秩的断 层部分；E对应了稀疏的噪声部分；

当低秩矩阵的奇异向量分布合理且稀疏矩阵的非零元素均匀分布, 该目标函数可转换为如下问题的求解；

其中，||A||_*是矩阵A的核范数，||E||₁表示矩阵E的1范 数，即σ_k(A)表示矩阵A的奇异值。

在步骤3中，对于公式(1-3)中的两个凸函数，套用ADMM原始的公式， 得到：

f(E)＝λ||E||₁ (1-4)

g(A)＝||A||_* (1-5)

其中，||·||_*表示张量核范数，采用张量奇异值分解T-SVD，写成如下形 式：

g(A)＝||A||_TNN (1-6)

利用增广拉格朗日方法求解，将约束项放入到目标函数中；目标函数 变为：

对于公式(1-7)，<W,D-A-E>表示张量的内积，ADMM算法首先固定其 他变量求解A；然后固定A和W，将E作为变量；最后更新拉格朗日乘子 W；β是惩罚参数，作用是平衡残差的权重；增广是指在加入了二次惩罚 性

对于最小化目标函数J的ADMM具体迭代求解方法为：

针对A变量：

针对E变量：

针对W变量：

引理1：假设矩阵那么秩为r的矩阵A的奇异值分解SVD表 示为：

A＝UE_rV^T，E_r＝diag({σ_i}_1≤i≤r) (1-14)

其中U和V各列是单位正交的，E_r是对角阵，对角元素是A的奇异值 σ_i；

对于目标函数它的奇异值软阈值算子如 公式(1-15)所示：

S_δ(W)＝U*S_δ(E_r)*V^T (1-15)

S_δ(E_r)＝diag(max((σ_i-δ),0)) (1-16)

在对RPCA模型的目标函数进行优化求解时，其ADMM算法子问题对其矩 阵核范数进行最小化，是采用了引理1中给出的矩阵的奇异值近似解；同 理，在Ho-RPCA模型的目标函数求解时，也需要对公式(1-11)中的张量 核范数求解最小化。

在步骤4中，采用基于TR-SVD分解的核范数近似求解对公式(1-11) 中的张量核范数进行最小化优化处理；由于在时域进行TR-SVD分解计算 相对较复杂，充分利用这种张量分解经过DFT变换所具有的性质，本次 在频域中执行TR-SVD计算。

步骤5包括：

(1)设定TR-SVD分解的截断项k，过采样参数p，并构建高斯随机张 量且使得G满足标准的正态分布；

(2)对于待处理的和均进行傅里叶变换处理获得和

(3)分别对不同的频率进行计算并进行正交三角(QR)分解， 获得正交酉矩阵并对正交酉矩阵做运算

(4)对做奇异值分解，获得

(5)根据设置的截断参数获得

(6)对每一个频率执行(3)、(4)、(5)步骤，并执行反傅里叶变换获得并输出结果： 和

将算法应用到实际地震数据中，并对得到的结果进行展示、分析和评 价，通过现实采集的数据来验证基于Ho-RPCA的断层增强方法。

本发明中的基于Ho-RPCA的地震断层增强方法，在矩阵稀疏低秩分解 模型的基础上，通过T-SVD张量分解进行高维拓展，将相干体数据中的断 层部分和噪声部分分离，由此达到相干体增强的目的。本发明从张量模型 的角度来对相干体进行增强，更好地保护了相干属性体的高维特性，可以 较好地去除噪声，增强后的断层连续性得到提升。对实际工区数据进行测 试，并与现有算法进行比较，本发明能够较好地增强相干体中的断层信息， 总体效果比较理想。

附图说明

图1为本发明的基于Ho-RPCA的地震断层增强方法的一具体实施例的 流程图。

图2为本发明的一具体实施例中GSM工区三维地震数据三维显示图；

图3为本发明的一具体实施例中GSM工区相干体切片的示意图；

图4为本发明的一具体实施例中本发明增强后的切片的示意图；

图5为本发明的一具体实施例中商业软件增强后的切片的示意图；

图6为本发明的一具体实施例中F3工区相干体切片的示意图；

图7为本发明的一具体实施例中增强后的切片的示意图；

图8为本发明的一具体实施例中商业软件增强后的切片的示意图。

具体实施方式

在步骤101，对输入三维相干体数据进行预处理,对低秩的断层部分和 稀疏的噪声部分进行分离。鲁棒主成分分析(RPCA)是将压缩感知理论在 矩阵上面进行的推广。RPCA的数学模型可以表示为，

D＝A+E (1-1)

如果将相干体图像用D来表示，由于断层在相干体数据中是规律的平 面结构，所以它对应的是低秩部分A。相干体中的噪声成分比较复杂，里 面主要包含了高斯噪声、相干噪声。这里对相干噪声进行简要说明：它的 形成与地震数据的采集以及数据前期的预处理有关；相干噪声往往分布在 断层面附近成絮状，它的图像灰度值要低于断层的图像灰度值；与高斯噪 声相比，相干噪声的分布具有规律性，高斯噪声在相干体图像中是满秩且 稀疏的，而相干噪声是满秩且近似稀疏。在相干体图像中，所有噪声的和 是满足近似稀疏性的，即噪声的和可以对应RPCA模型中的稀疏部分E。这 样，就将三维相干体数据D用低秩部分A和稀疏部分E表示出来。

在步骤102，构建基于Ho-RPCA的目标函数。张量作为向量和矩阵的高 阶推广，能更好地表达高阶数据的复杂本质结构。目前，Ho-RPCA在一些 视频去噪恢复的问题上已经取得了很好的效果，视频去噪问题与相干体数 据的断层增强问题是可以相类比，所以本发明采用低秩张量恢复模型 (Ho-RPCA)来对三维相干数据体进行断层增强在理论上是可行的。

根据Ho-RPCA对相干体三维图像建模，类似于RPCA的图像去噪模型， 有如下公式：

其中，rank(A)表示对矩阵A取秩；||E||₀是E的0范数，表示矩阵E中非 零元素的个数。参数λ是用来平衡低秩项和稀疏项的权重参数，一般λ>0。 D是待去噪增强的相干体数据，A对应了低秩的断 层部分；E对应了稀疏的噪声部分。然而，由于目标函数中rank(A)和||E||₀都是非线性非凸的组合优化函数，对上述问题的求解是十分困难的。

当低秩矩阵的奇异向量分布合理且稀疏矩阵的非零元素均匀分布, 该目标函数可以转换为如下问题的求解。

其中，||A||_*是矩阵A的核范数，||E||₁表示矩阵E的1范 数，即σ_k(A)表示矩阵A的奇异值。

在步骤103，通过交替方向乘子法(Alternating Direction Method ofMultipliers,ADMM)框架来求解该最优化问题。对于公式(1-3)中的 两个凸函数，套用ADMM原始的公式，得到：

f(E)＝λ||E||₁ (1-4)

g(A)＝||A||_* (1-5)

其中，||·||_*表示张量核范数(tensor nuclear norm TNN)。因为本发 明采用的是张量奇异值分解(T-SVD)分解，所以可以写成如下形式：

g(A)＝||A||_TNN (1-6)

利用增广拉格朗日方法求解，将约束项放入到目标函数中。目标函数 变为：

对于公式(1-7)，<W,D-A-E>表示张量的内积，ADMM算法首先固定其 他变量求解A；然后固定A和W，将E作为变量；最后更新拉格朗日乘子 W。β是惩罚参数，作用是平衡残差的权重。“增广”是指在加入了二次 惩罚性该项的作用是使得目标函数更加鲁棒和收敛。

对于最小化目标函数J的ADMM具体迭代求解方法为：

针对A变量：

针对E变量：

针对W变量：

引理1：假设矩阵那么秩为r的矩阵A的奇异值分解(SVD) 表示为：

A＝UE_rV^T，E_r＝diag({σ_i}_1≤i≤r) (1-14)

其中U和V各列是单位正交的，E_r是对角阵，对角元元素是A的奇异 值σ_i。

对于目标函数它的奇异值软阈值算子如 公式(1-15)所示。

S_δ(W)＝U*S_δ(E_r)*V^T (1-15)

S_δ(E_r)＝diag(max((σ_i-δ),0)) (1-16)

在对RPCA模型的目标函数进行优化求解时，其ADMM算法子问题对其矩 阵核范数进行最小化，是采用了引理1中给出的矩阵的奇异值近似解。同 理，在Ho-RPCA模型的目标函数求解时，也需要对公式(1-11)中的张量 核范数求解最小化。由于待分解的相干体数据是一个三维的数据体，在对 其进行张量奇异值近似时，不能使用矩阵的SVD(奇异值分解)。本发明引 入T-SVD(张量奇异值分解)分解来求取三阶张量的奇异值近似,并在T-SVD 分解中加入高斯随机张量来提升奇异值分解速度。

在步骤104，对张量核范数进行最小化，采用基于张量随机奇异值 (TR-SVD)分解的核范数近似求解。

本发明引入T-SVD分解方法来求解如下的问题：

但当输入数据D很大时,T-SVD分解速度很慢，于是本发明引入 TR-SVD(张量随机奇异值分解)方法来代替T-SVD分解。TR-SVD是将R-SVD 扩张到T-SVD的过程中，TR-SVD的方法是为了找到张量分解的 最优近似。TR-SVD求解过程分为两步。(1)寻找张量Q使得A≈Q*Q^T*A；(2) 利用低秩表示来进行TR-SVD分解。

在采用TR-SVD求解张量核范数最小化时，当计算在时域上进行，会产 生大量的t-product运算，计算量较大。而张量在时域上的t-product运算， 可以通过离散傅里叶变换转换到频域上，此时计算复杂度高的t-product 运算可以由频域上的循环卷积来替换，这将大大减少张量TR-SVD分解的的 计算量。

TR-SVD频域形式的计算流程为：

(1)设定TR-SVD分解的截断项k，过采样参数p，并构建高斯随机张 且使得G满足标准的正态分布；

(2)对于待处理的和均进行傅里叶变换处理获得和

(3)分别对不同的频率进行计算并进行正交三角(QR)分解， 获得正交酉矩阵并对正交酉矩阵做运算

(4)对做奇异值分解，获得

(5)根据设置的截断参数获得

(6)对每一个频率执行(3)、(4)、(5)步骤，并执行反傅里叶变换获得并输出结果： 和

将算法应用到实际地震数据中，并对得到的结果进行展示、分析和评 价，通过现实采集的数据来验证基于Ho-RPCA的断层增强方法。

图2为本发明的一具体实施例中GSM工区三维地震数据三维显示图；

在应用本发明的一具体实例中，选用的测试工区是一组是GSM工区地 震数据，另一组是位于荷兰北海的F3工区数据。本发明处理的数据是三维 相干体数据，且是对整个数据体进行处理。本节对整个工区范围内进行基 于Ho-RPCA的断层增强方法。一组是GSM工区地震数据，GSM工区的三维地 震数据三维显示图如图2所示。首先，在工区的三维振幅数据体中，截取 范围在Inline7600～8400、Crossline3500～4300、Time500～3000ms区域的振幅数据 体，然后对该区域的数据体用叠后C3算法计算相干体；再对该工区的相干 体进行断层增强。

图3为GSM工区的相干数据体在时间t＝1500ms处的切片。在该切片图 像上可以大致看出断层面貌分布的趋势：其中，白色像素区域对应的是背 景，代表了均匀分布的地层结构；而黑色像素区域则是断层和少部分其他 的地质构造(如河道等)，在断层附近的絮状阴影就是相干噪声。

图4为GSM工区数据相干数据体在时间t＝1500ms处的相干体时间切片 断层增强效果图，可以看到断层线条附近的成絮状的相干噪声已明显减 少，且基本消除了背景中的高斯噪声；图5为商业软件AASPI增强后的切片 的示意图。比较图4和图5，可以发现在视觉效果上，本发明在断层线条与 背景区域的对比度方面效果提升显著，西北区域的裂缝发育特征也更加清 晰。

图6是F3工区在t＝1200ms处的时间切片图，在图7所示AASPI软件增 强后的相干体切片中，图像黑圈区域里的断层仍是难以区分，在断层线条 附近依然残留有一些的相干噪声。而在图8所示的本节方法增强后的相干 体切片中，较好的抑制了图像黑圈区域里的相干噪声，断层线条与背景已 经可以较好的区分开，本方法取得了更为理想的断层增强效果。

本发明将低秩张量恢复理论应用于相干体三维图像的断层增强中，提 出了基于Ho-RPCA的断层增强方法。将RPCA模型拓展到高维的 Ho-RPCA模型，构建了Ho-RPCA模型的目标函数，并通过ADMM算法框 架来求解该凸优化问题，在求解过程中，采用了TR-SVD来对张量范数进 行求解。最后给出了基于Ho-RPCA的断层增强方法的流程。从张量模型 的角度来对相干体进行增强，更好地保护了相干属性体的高维特性，可以 较好地去除噪声，增强后的断层连续性得到提升。对实际工区数据进行测 试，并与现有算法进行比较，本发明能够较好地增强相干体中的断层，总 体效果比较理想。

Claims (6)

1.基于Ho-RPCA的地震断层增强方法，其特征在于，该基于Ho-RPCA的地震断层增强方法包括：

步骤1，对输入三维相干体数据进行预处理,对低秩的断层部分和稀疏的噪声部分进行分离；

步骤2，构建基于Ho-RPCA的目标函数；

步骤3，通过交替方向乘子法算法框架来求解该最优化问题；

步骤4，对张量核范数进行最小化，采用基于张量随机奇异值分解的核范数近似求解。

2.根据权利要求1所述的基于Ho-RPCA的地震断层增强方法，其特征在于，在步骤1中，采用鲁棒主成分分析法对输入三维相干体数据进行预处理，鲁棒主成分分析法的数学模型可表示为：

D＝A+E (1-1)

如果将相干体图像用D来表示，由于断层在相干体数据中是规律的平面结构，所以它对应的是低秩部分A；在相干体图像中，所有噪声的和是满足近似稀疏性的，即噪声的和可以对应鲁棒主成分分析法模型中的稀疏部分E，这样，就将三维相干体数据D用低秩部分A和稀疏部分E表示出来。

3.根据权利要求1所述的基于Ho-RPCA的地震断层增强方法，其特征在于，在步骤2中，采用低秩张量恢复模型Ho-RPCA来对三维相干数据体进行断层增强；

根据Ho-RPCA对相干体三维图像建模，类似于鲁棒主成分分析法的图像去噪模型，有如下公式：

其中，rank(A)表示对矩阵A取秩；||E||₀是E的0范数，表示矩阵E中非零元素的个数；参数λ是用来平衡低秩项和稀疏项的权重参数，一般λ>0；D,A,D是待去噪增强的相干体数据，A对应了低秩的断层部分；E对应了稀疏的噪声部分；

当低秩矩阵的奇异向量分布合理且稀疏矩阵的非零元素均匀分布,该目标函数可转换为如下问题的求解；

其中，||A||_*是矩阵A的核范数，||E||₁表示矩阵E的1范数，即σ_k(A)表示矩阵A的奇异值。

4.根据权利要求3所述的基于Ho-RPCA的地震断层增强方法，其特征在于，在步骤3中，对于公式(1-3)中的两个凸函数，套用ADMM原始的公式，得到：

f(E)＝λ||E||₁ (1-4)

g(A)＝||A||_* (1-5)

其中，||·||_*表示张量核范数，采用张量奇异值分解T-SVD，写成如下形式：

g(A)＝||A||_TNN (1-6)

利用增广拉格朗日方法求解，将约束项放入到目标函数中；目标函数变为：

对于公式(1-7)，<W,D-A-E>表示张量的内积，ADMM算法首先固定其他变量求解A；然后固定A和W，将E作为变量；最后更新拉格朗日乘子W；β是惩罚参数，作用是平衡残差的权重；增广是指在加入了二次惩罚性

对于最小化目标函数J的ADMM具体迭代求解方法为：

针对A变量：

针对E变量：

针对W变量：

引理1：假设矩阵那么秩为r的矩阵A的奇异值分解SVD表示为：

A＝UE_rV^T，E_r＝diag({σ_i}_1≤i≤r) (1-14)

其中U和V各列是单位正交的，E_r是对角阵，对角元素是A的奇异值σ_i；

对于目标函数它的奇异值软阈值算子如公式(1-15)所示：

S_δ(W)＝U*S_δ(E_r)*V^T (1-15) S_δ(E_r)＝diag(max((σ_i-δ),0)) (1-16)

在对RPCA模型的目标函数进行优化求解时，其ADMM算法子问题对其矩阵核范数进行最小化，是采用了引理1中给出的矩阵的奇异值近似解；同理，在Ho-RPCA模型的目标函数求解时，也需要对公式(1-11)中的张量核范数求解最小化。

5.根据权利要求4所述的基于Ho-RPCA的地震断层增强方法，其特征在于，在步骤4中，采用基于TR-SVD分解的核范数近似求解对公式(1-11)中的张量核范数进行最小化优化处理；由于在时域进行TR-SVD分解计算相对较复杂，充分利用这种张量分解经过DFT变换所具有的性质，本次在频域中执行TR-SVD计算。

6.根据权利要求5所述的基于Ho-RPCA的地震断层增强方法，其特征在于，步骤5包括：

(1)设定TR-SVD分解的截断项k，过采样参数p，并构建高斯随机张量且使得G满足标准的正态分布；

(2)对于待处理的和均进行傅里叶变换处理获得和

(3)分别对不同的频率进行计算并进行正交三角分解，获得正交酉矩阵并对正交酉矩阵做运算

(4)对做奇异值分解，获得

(5)根据设置的截断参数获得

(6)对每一个频率执行(3)、(4)、(5)步骤，并执行反傅里叶变换获得并输出结果： 和

将算法应用到实际地震数据中，并对得到的结果进行展示、分析和评价，通过现实采集的数据来验证基于Ho-RPCA的断层增强方法。