CN108490486B - 一种地震资料反演的方法及装置、设备 - Google Patents

Abstract

本发明涉及一种地震资料反演的方法及装置、设备。该方法建立基于l₂范数约束的弹性阻抗反演最优化模型。提出先借助阈值滤波子和L‑曲线准则确定最优正则参数，然后利用Tikhonov滤波子求解最优化模型的极小解的方法，能有效地提高计算精度以及解的稳定性。通过低信噪比理论模型和实际地震资料对算法进行测试。低信噪比理论模型的反演结果与模型正演情况基本吻合；而基于实际地震资料的弹性阻抗反演结果与该地区岩石物理分析结果相符合，进一步获得的脆性指数分布的预测结果也与钻井吻合。该算法具有良好的稳定性和抗噪性，有实际应用推广的前景。

Description

一种地震资料反演的方法及装置、设备

技术领域

本发明涉及领域地震勘探领域，具体涉及一种应用于地震波弹性阻抗反演中的双滤子正则化方法。

背景技术

页岩气作为一种清洁、高效的非常规能源，它的开发和利用对发展低碳经济、优化能源结构具有重要意义。岩石脆性作为遴选高品质页岩的重要评价指标，是页岩气储层描述的主要内容之一，而脆性指数是岩石脆性的主要评价指标。通常，脆性指数的选取是基于模型正演敏感性分析，利用纵、横波阻抗、拉梅系数、杨氏模量、泊松比等弹性参数，获得不同方法的脆性指数结果，最终优选出敏感的脆性指数。

目前，叠后声波阻抗反演、AVO反演和弹性阻抗反演是提取上述弹性参数的主要方法。叠后声波阻抗反演方法结果单一、不能充分地利用叠前信息。AVO反演方法对数据信噪比要求高，数据信噪比低而导致地震子波随偏移距变化的情况会降低AVO反演结果的准确性。与叠后声波阻抗反演和AVO反演相比，Connolly(1999)提出的弹性阻抗反演方法，考虑AVO效应，利用部分叠加来提高数据信燥比，能对近、远偏移距资料进行有效标定并获得更为丰富、稳定、可靠的弹性参数反演结果，有利于页岩敏感性脆性指数的选取。

弹性阻抗反演是通过求解离散线性算子方程相关的最小二乘问题来实现的。线性反问题可用如下形式描述：

Az＝u z∈F，u∈U (1)

式中，A是已知算子，u是实测数据，z是待求的模型参数，F、U可称为“空间”。弹性阻抗反演是Hadamard意义下的不适定问题，即不能同时满足解的存在性、唯一性和稳定性的要求。因此直接求解离散线性算子方程(公式1)相关的最小二乘问题会带来不稳定的计算结果，而正则化方法是求解不适定问题的有效方法。正则化方法是用一簇与原问题相“邻近”的适定问题的解去逼近原问题的真实解，即通过构造正则算子得到原问题稳定的近似解，如公式(2)所示。

Az+αz＝uα＞0 (2)

式中，α是正则参数，控制与原问题的“临近程度”。

正则化方法包括Tikhonov正则化方法、Landweber迭代法、Newton型方法、共轭梯度法等。其中，由苏联科学院院士Tikhonov所提出的Tikhonov正则化方法是通过引入光滑泛函来构造正则算子，进而利用其极小解去逼近原问题的真实解，应用最为广泛。

最优正则参数的选取是正则化方法的核心问题之一。若正则参数α取值过大会导致公式(2)与原问题(公式1)相差太远；反之，α取值太小会导致公式(2)过多的“继承”原问题的不适定性而导致解的不稳定。正则参数的选取方法通常分为先验和后验两大类。选取正则参数的先验方法具有理论分析的意义，但是在实际中常常难以实现。后验方法有如下几种形式：偏差原理、广义偏差原理、误差极小化准则、Arcangeli准则等，这些方法需要预先估计原始数据的误差水平。针对原始资料误差水平未知的情况，Tikhonov等提出拟最优准则、Hansen等提出L-曲线准则、Golub等提出广义交叉校验准则。其中，利用L-曲线准则确定最优正则参数的方法在图像恢复、生命科学、遥感技术等领域已经得到广泛应用。

发明内容

根据地震资料数据量大，信噪比差，资料误差水平难以准确估计的情况，为克服地震波弹性阻抗反演问题的不适定性，将解的l₂范数最小作为约束条件加入到原问题；提出用L-曲线准则选择最优正则化参数，平衡目标函数和约束项；在数值实现上，把Tikhonov正则化与截断奇异值正则化方法结合起来，利用双滤子技巧，获得较好的弹性阻抗反演结果。

步骤1：对投影模型空间到观测空间的算子矩阵W进行奇异值分解，[U，∑，V]＝svd(W)，其中∑＝diag(σ₁，σ₂，…σ_n)是矩阵奇异值，呈降序排列；U、V分别是m×m和n×n阶正交矩阵，U＝(u₁，u₂，…，u_m)，V＝(v₁，v₂，…，v_n)，

即秩为r的m*n阶矩阵，r是W的数值秩，W的奇异***可以记为{σ_i，u_i，v_j}；

步骤2：计算解的范数

和残差的范数

其中，σ.是矩阵W的奇异值，α是正则参数，d＝WR+N，式中，d＝(d₁，d₂，…，d_m)^T，R＝(R₁，R₂，…，R_n)^T，N＝(N₁，N₂，…N_m)^T；

步骤3：根据步骤2中解的范数和残差的范数绘制关系曲线(log||R_α||₂，log||WR_α-d||₂)，计算最优正则解。

所述双滤子分别为Tikhonov滤波子和阈值滤波子。

所述双滤子分别为Tikhonov滤波子和阈值滤波子的获取方法如下：采用Tikhonov正则化方法将弹性阻抗反演问题表达式简化为||WR-d||₂→min形式，建立基于带l₂范数约束的最优化模型：

其中，K是规范化矩阵，为正定或半正定矩阵，R⁰为初始解，令其为0；α＞0是正则化参数，最优化模型的极小解根据Picard定理可以写成下面的形式：

其中，f_i(α)是滤波子，Tikhonov滤波子在K为单位矩阵，即K＝I的条件下，f_i(α)选取如下形式：

阈值滤波子则为：

所述步骤3中计算最优正则解的过程中需要先确定最优正则化参数：绘制关系曲线(log||R_α||₂，log||WR_α-d||_２)，该曲线在log-log尺度下会存在一个明显的拐角，最优化正则参数所对应的值位于该曲线曲率最大值附近，通过寻找该曲线最大曲率并依据硬阈值滤波子来确定最优正则化参数。

所述曲线的曲率函数构建方法如下：令ρ(α)＝log||WR_α-d||₂，θ(α)＝log||R_α||_２，有关正则参数α的非线性函数如下

其中，ρ′，ρ″和θ′，θ″分别为函数ρ(α)和θ(α)的一阶和二阶导数，

)称作参数目的曲率函数。

所述最优正则解的计算方法如下：利用Picard定理，采用Tikhonov滤波子

获得正则解

一种地震资料反演的装置，其特征在于，所述装置包括：

奇异值分解模块：对投影模型空间到观测空间的算子矩阵W进行奇异值分解，[U，∑，V]＝svd(W)，其中∑＝diag(σ₁，σ₂，…σ_n)是矩阵奇异值，呈降序排列；U、V分别是m×m和n×n阶正交矩阵，U＝(u₁，u₂，…，u_m)，V＝(v₁，v₂，…，v_n)，

即秩为r的m*n阶矩阵，r是W的数值秩，W的奇异***可以记为{σ_i，u_i，v_j}；

范数计算模块：计算解的范数

和残差的范数

其中，σ.是矩阵W的奇异值，α是正则参数，d＝WR+N，式中，d＝(d₁，d_２，…，d_m)^T，R＝(R₁，R₂，…，R_n)^T，N＝(N₁，N₂，…N_m)^T；

最优正则解获取模块：根据解的范数和残差的范数(log||R_α||₂，log||WR_α-d||₂)绘制关系曲线，计算最优正则解。

一种地震资料反演的设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述程序时可实现如上所述的方法。

本发明的有益技术效果如下：

(1)建立基于l₂范数约束的弹性阻抗反演最优化模型。

(2)提出先借助阈值滤波子和L-曲线准则确定最优正则参数，然后利用Tikhonov滤波子求解最优化模型的极小解的方法，能有效地提高计算精度以及解的稳定性。

(3)通过低信噪比理论模型和实际地震资料对算法进行测试。低信噪比理论模型的反演结果与模型正演情况基本吻合；而基于实际地震资料的弹性阻抗反演结果与该地区岩石物理分析结果相符合，进一步获得的脆性指数分布的预测结果也与钻井吻合。该算法具有良好的稳定性和抗噪性，有实际应用推广的前景。

附图说明

图1地层的弹性阻抗随入射角的变化曲线(0-30°范围内)

图2模型正演的弹性阻抗剖面

图3模型正演的地震剖面(无噪音)

图4模型正演的地震剖面(s/n＝5)

图5 L-曲线选取正则参数

图6弹性阻抗反演结果剖面

图7反演结果的误差(％)

图8 L80主测线25°入射角的地震剖面(1600-1840ms)

图9页岩和灰岩弹性阻抗特征分析

图10 L-曲线选取正则参数

图11弹性阻抗反演结果

图12基于弹性反演结果所计算的敏感脆性指数(Brit4)

具体实施方式

实测地震记录可视为地层反射系数与震源子波的卷积，即

式中，t_i表示时间延迟，ω表示震源子波，d(t)表示地震信号，N(t)表示可加性噪音干扰。

公式(3)可以写成矩阵-向量方程的形式：

d＝WR+N (4)

式中，d＝(d₁，d₂，…，d_m)^T，R＝(R₁，R₂，…，R_n)^T，N＝(N₁，N₂，…N_m)^T以及

W是将模型空间投影到观测空间的算子。线性反问题(公式4)的求解是弹性阻抗反演的一个关键。

Tikhonov正则化方法是根据Lagrange乘子法或优化理论中构造泛函的思想，通过引入光滑泛函来构造正则算子，进而对不适定性反问题进行求解。本文利用上述Tikhonov正则化思想，将线性反问题(4)的求解简化为数据拟合问题，即||WR-d||₂→min的形式，并建立基于l₂范数约束的最优化模型：

式中，K是规范化矩阵，通常是一个(半)正定矩阵，在本文中，令K为一个恒等算子；R⁰是初始解(提供了解空间的先验信息)，在本文中，令R⁰＝0；α＞0是正则化参数，平衡残差项

与附加约束项

以避免过正则化或欠正则化。

在算子广义逆和谱分析的理论框架下对最优化模型(公式6)进行求解。对于公式(6)中的矩阵算子

即秩为r的m*n阶矩阵，进行奇异值分解：

W＝U∑V (7)

式中，U、V为酉矩阵，U＝(u₁，u₂，…，u_m)，V＝(v₁，v_２，…，v_n)，i＝1，2，…m，j＝1，2，…n；∑形式如下：

其中，∑_r＝diag(σ₁，σ₂，…，σ_r)，σ₁≥σ₂≥…≥σ_r。矩阵算子W的奇异***可以记为{σ_i，u_i，v_j]。基于Moore-Penrose广义逆，根据Picard定理可将最优化问题(6)的极小解写成下面的形式：

当数据存在噪音干扰，σ_i→0(i→r)时，因式

噪音干扰所导致的误差将会被放大。此时，由上述公式(8)求得的近似解会不稳定。通过引入合适的滤波器可以将1/σ_i的放大作用衰减下来：

其中，f_i(α)是滤波子。为了获得原问题的稳定近似解，需要选择合适的滤波子。本文针对弹性阻抗反演问题考虑两种滤波子。

Tikhonov滤波子。在K是单位矩阵(K＝1)的情况下，Tikhonov滤波子具体有如下形式：

(11)此时，极小解(公式10)可以表示为

由于矩阵W的病态性，不适定问题(4)的解可以利用W的数值秩r_δ来表示。设B_δ为W的扰动矩阵，则数值秩r_δ满足r_δ＝min{rank(B_δ+A)：B_δ∈R^m×n，||B_δ||₂≤δ} (12)

注意到σ_i是矩阵W的奇异值，则

因此，存在参数

使得

阈值滤波子。根据以上讨论，公式(10)的滤波函数f_i(α)也可以取截断奇异值分解的形式：

在数值实现上，本文把上述Tikhonov滤波子和阈值滤波子结合起来，形成双滤波子正则化。由于引起弹性阻抗反演问题求解不稳定的因素主要是正演算子小奇异值以及地震数据中的噪音干扰，因此，在确定最优正则参数α时，借助阈值滤波子把小奇异值滤掉，以提高算法的稳定性。然后，再利用Tikhonov滤波子，求解最优化模型(公式6)的极小解。

在正则化方法中，最优的正则参数α起到平衡残差项和有关解约束条件的作用。Hansen等人正是基于这种考虑，提出L-曲线准则，即在双对数坐标***下借助残差项||WR_α-d||_２和约束条件||K(R_α-R⁰)||₂随正则参数α变化的曲线来确定最优正则参数的方法。该方法因曲线在对数尺度中会呈现“L”形状而得名；其关键参数是L-曲线隅角。Hanke等^[22]定义L-曲线隅角为其在对数尺度下的最大曲率。

令ρ(α)＝log||WR_α-d||_２，θ(α)＝log||R_α||_２，定义一个有关正则参数α的非线性函数：

其中，ρ′，ρ″和θ′，θ″分别为函数ρ(α)和θ(α)的一阶和二阶导数，

称作参数α的曲率函数。

下面给出利用Tikhonov方法和截断奇异值方法相结合，以及L-曲线准则求取弹性阻抗反演问题的双滤子正则化算法。

双滤子正则化算法

(1)算子矩阵W作奇异值分解：

[U，∑，V]＝svd(W) (15)

其中，∑＝diag(σ₁，σ₂，…σ_n)是矩阵奇异值，呈降序排列；U、V分别是m×m和n×n阶正交矩阵。

(2)对σ_r≤α≤σ₁(r≤n)，r是W的数值秩，借助阈值滤波子，计算反问题逼近解的范数：

以及残差的范数：

(3)以图形形式(log||R_α||_２，log||WR_α-d||_２)，绘出L曲线确定曲线隅角，隅角处对应的正则参数α就是最优正则参数。

(4)借助Tikhonov滤波子，获得反问题的正则解

首先，利用Castagna和Smith(1994)提出的第三类经典AVO模型^[27]进行试算。如图2所示，构造一个三层的楔状构造模型，时间范围为0-250ms，共140个地震道，地层岩性从上至下分别为页岩、含气砂岩、含水砂岩，各地层的弹性参数分布情况如表1所示。图1是三个地层弹性阻抗值随入射角的变化曲线(0-30°)。由图可见，页岩层和含水砂岩层的弹性阻抗值随入射角而增大，含气砂岩层的弹性阻抗值随入射角而减小。也就是说，当入射角为30°时，各反射界面上、下地层的弹性阻抗差异最大。

表1 Castagna和Smith(1994)的第三类经典AVO模型

为了方便试算结果的分析对比，本文主要研究模型30°入射角的情况。图2是模型正演的弹性阻抗剖面。正演过程采用零相位雷克子波：

其中，f_m是子波主频，f_m＝30Hz；t是时间采样间隔，t＝1ms；子波样点个数是101。图3、图4分别是无噪声、信噪比为5的情况下模型正演的地震剖面。由地震剖面可见，页岩和含气砂岩之间的反射同相轴为负极性弱反射；当数据中存在噪音干扰时，难以在地震剖面上准确定位页岩和含气砂岩地层的分界面。利用双滤子正则化算法对图4的地震数据进行反演。图5是该反演过程的L-曲线示意图，横、纵坐标分别是对数尺度下反问题逼近解的范数log||R_α||₂和残差项的范数log||WR_α-d||_２。根据Hanke等的定义，利用公式(10)确定L-曲线的隅角(曲线红色标记处)以及它所对应的正则参数α＝0。在图4中，最优正则参数α＝0.8846。图6是弹性阻抗反演结果剖面。为了验证双滤子正则化算法的精度，图7给出反演结果与真实模型的误差分布情况。从图7可以看出，页岩地层反演结果的误差在7％左右；而含气砂岩、含水砂岩的反演结果误差相对较高，分别为10％和13％。

利用某油田的实际三维地震资料对双滤子正则化算法进行测试。本次反演目的层为泥盆系页岩气层，构造较为平缓，顶面埋深2500-4500m，沉积于开阔海域环境，最大海进盆地快速充填期，上覆和下伏地层均为泥灰岩。从图8所示的80主测线25度入射角的地震剖面可以看出，研究目的层地层为一低频弱振幅波谷反射特征，反射特征清楚，横向连续性好。

通过岩石物理分析，页岩和灰岩阻抗差异明显，并且随入射角度增大差异增大。如图9所示，在25度入射角情况下，页岩弹性阻抗值范围是9x10³-1.1x10⁴g/cm³*m/s，灰岩弹性阻抗值范围是1.1x10⁴-1.7x10⁴g/cm³*m/s，说明弹性阻抗能有效地区分页岩和泥岩。工区地震数据主频28Hz，频带范围是8-100Hz。目的层页岩储层有效厚度20-50m。上述3个条件满足弹性阻抗反演的要求。图10是实际数据反演过程中的L-曲线法示意图。横、纵坐标分别是对数尺度下反问题逼近解的范数log||R_α||₂和残差项的范数log||WR_α-d||₂。根据Hanke等^[11]的定义，利用公式(10)确定L-曲线隅角(曲线红色标记处)以及它所对应的正则参数α＝0.5469。

图11是弹性阻抗反演结果剖面。借助弹性阻抗剖面可以明显区分页岩和灰岩。目的层页岩是剖面1700ms-1750ms处的深蓝色部分，表现为明显的低阻抗特征，厚度大，具备较好的追踪解释条件，厚度在15-20米。剖面1650ms-1700ms的浅蓝色区域为灰岩，弹性阻抗值范围是1.15x10⁴-1.3x10⁴g/cm³*m/s。弹性阻抗反演结果与图9所示的该地区岩石物理分析结果相符合。基于弹性阻抗反演结果可获得纵、横波阻抗、杨氏模量、泊松比等弹性参数。对比不同方法的脆性指数结果

最终优选出此数据的敏感脆性指数是Brit₄＝Eρ。图12是敏感脆性指数Brit₄的计算结果，如图12所示，在页岩层段内部脆性指数较高的区域即黄色到红色区域为页岩甜点区，内部细节较为丰富具有很好的指示作用。根据反演结果与测线附近井资料的对比分析，目的层段弹性阻抗均值与井资料的误差为11％，脆性指数分布也与钻井吻合，有效地满足了实际生产的需求。

以上所述仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内所作的任何修改、等同替换、改进等，均包含在本发明的保护范围内。

Claims (6)

1.一种地震资料反演的双滤子正则化方法，具体步骤如下：

步骤1：对投影模型空间到观测空间的算子矩阵W进行奇异值分解，[U，∑，V]＝svd(W)，其中Σ＝diag(σ₁，σ₂，…σ_n)是矩阵奇异值，呈降序排列；U、V分别是m×m和n×n阶正交矩阵，U=(u₁，u₂，…，u_m)，V＝(v₁，v₂，…，v_n)，

即秩为r的m*n阶矩阵，r是W的数值秩，W的奇异***可以记为{σ_i，u_i，v_j}，其中i＝1，2，…m，j＝1，2，…n；

步骤2：计算解的范数

和残差的范数

其中，σ是矩阵W的奇异值，α是正则参数，d＝WR+N，式中，d＝(d₁，d₂，…，d_m)^T，R＝(R₁，R₂，…，R_n)^T，N＝(N₁，N₂，…N_m)^T；

步骤3：根据步骤2中解的范数和残差的范数绘制(log||R_a||₂，log||WR_a-d||₂)关系曲线，计算最优正则解；

所述双滤子分别为Tikhonov滤波子和阈值滤波子，其获取方法如下：采用Tikhonov正则化方法将弹性阻抗反演问题表达式简化为||WR-d||₂→min形式，建立基于带l₂范数约束的最优化模型：

其中，K是规范化矩阵，为正定或半正定矩阵，R⁰为初始解，令其为0；α＞0是正则化参数，最优化模型的极小解根据Picard定理可以写成下面的形式:

其中，f_i(α)是滤波子，Tikhonov滤波子在K为单位矩阵，即K＝I的条件下，f_i(α)选取如下形式：

阈值滤波子则为：

2.根据权利要求1所述的方法，其特征在于，步骤3中计算最优正则解的过程中需要先确定最优正则化参数：绘制(log||R_α||₂，log||WR_α-d||₂)关系曲线，该曲线在log-log尺度下会存在一个明显的拐角，最优化正则参数所对应的值位于该曲线曲率最大值附近，通过寻找该曲线最大曲率并依据硬阈值滤波子来确定最优正则化参数。

3.根据权利要求2所述的方法，其特征在于，曲线的曲率函数构建方法如下：令ρ(α)＝log||WR_α-d||₂，θ(α)＝log||R_α||₂，有关正则参数α的非线性函数如下

其中，ρ′，ρ″和θ′，θ″分别为函数ρ(α)和θ(α)的一阶和二阶导数，

称作参数α的曲率函数。

4.根据权利要求3所述的方法，其特征在于，所述最优正则解的计算方法如下：利用Picard定理，采用Tikhonov滤波子

获得正则解

5.一种地震资料反演的装置，其特征在于，所述装置包括：

奇异值分解模块：对投影模型空间到观测空间的算子矩阵W进行奇异值分解，[U，∑，V]＝svd(W)，其中Σ＝diag(σ₁，σ₂，…σ_n)是矩阵奇异值，呈降序排列；U、V分别是m×m和n×n阶正交矩阵，U＝(u₁，u₂，…，u_m)，V＝(v₁，v₂，…，v_n)，

即秩为r的m*n阶矩阵，r是W的数值秩，W的奇异***可以记为{σ_i，u_i，v_j}，其中i＝1，2，…m，j＝1，2，…n；

范数计算模块：计算解的范数

和残差的范数

其中，σ是矩阵W的奇异值，α是正则参数，d＝WR+N，式中，d＝(d₁，d₂，…，d_m)^T，R＝(R₁，R₂，…，R_n)^T，N＝(N₁，N₂，…N_m)^T；

最优正则解获取模块：根据上述解的范数和残差的范数绘制(log||R_a||₂，log||WR_a-d||₂)关系曲线，计算最优正则解；

双滤子分别为Tikhonov滤波子和阈值滤波子，其获取方法如下：采用Tikhonov正则化方法将弹性阻抗反演问题表达式简化为||WR-d||₂→min形式，建立基于带l₂范数约束的最优化模型：

其中，K是规范化矩阵，为正定或半正定矩阵，R⁰为初始解，令其为0；α＞0是正则化参数，最优化模型的极小解根据Picard定理可以写成下面的形式:

其中，f_i(α)是滤波子，Tikhonov滤波子在K为单位矩阵，即K＝I的条件下，f_i(α)选取如下形式：

阈值滤波子则为：

6.一种地震资料反演的设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述程序时可实现如权利要求1-4中任一项的方法。

Images