CN107589448B - 一种多道地震记录反射系数序列同时反演方法 - Google Patents

Abstract

本发明公开了一种多道地震记录反射系数序列同时反演方法。该方法充分考虑地层反射系数序列的稀疏特征以及地震道的空间横向连续性，利用参数化指数函数来度量反射系数序列的稀疏性，建立相应的约束优化目标函数。首先，基于褶积模型的特点和反射系数的稀疏性假设，在不丢失任何信息的前提下，设计降维预处理矩阵，降低约束目标函数的维度；然后，利用信赖域序列二次规划算法对优化问题进行高效求解。该算法具有全局收敛性和超线性收敛率，计算效率高，且具有较高反演精度。相比于传统的单道地震记录反褶积算法，本发明算法由于考虑了地震道的空间相关性，反演结果更稳定、可靠。

Description

一种多道地震记录反射系数序列同时反演方法

技术领域

本发明属于地球物理勘探技术领域，涉及一种地震记录反射系数序列反演方法，特别涉及一种用于多道地震记录反射系数序列同时反演方法。

背景技术

在地层－岩性油气藏中，薄互层型储层是主要储层类型之一。由于激发地震所用人工震源频带范围有限，薄层产生的反射波走时差别小，从而使得反射波相互叠加干涉，接收到的地震记录无法直接分辨出薄互层结构，为薄互层型储层识别带来诸多困难。因此，有必要提高反射地震记录的分辨率，实现对薄互层地层结构的高精度刻画。

反褶积类方法通过压缩地震子波，进而获取地层反射系数序列，是地震资料处理中提高纵向分辨率、精细刻画地层结构常用方法之一。反褶积的理论基础一般是建立在满足若干基本假设的褶积模型之上。根据对地震资料的不同假设，地震记录反褶积方法种类较多，比如维纳反褶积、最小熵反褶积、最小平方反褶积、预测反褶积、互信息率反褶积、统计反褶积和稀疏脉冲反褶积等。

稀疏脉冲反褶积是建立在地层反射系数满足稀疏性假设的条件之上的，即大部分反射系数幅值趋近于零，剩余的大幅值反射系数对应着主要的岩性变化界面。基于对大量的实际测井数据分析，人们认识到在一定的地质沉积背景下，地层反射系数可用稀疏脉冲序列进行表示。统计上，可用 Bernoulli-Gaussian过程来描述服从稀疏分布的地层反射系数序列。在此基础之上，已经发展了一系列稀疏脉冲反褶积方法。Taylor等于1979年率先探讨了利用l_p范数来实现地震稀疏脉冲反褶积方法。Oldenbug等利用l₁范数反演层状地质结构的反射系数序列，取得较好的效果。Kormylo和Mendel在稀疏脉冲假设条件下，提出最大似然反褶积技术。Lavielle提出了基于 Bernoulli-Gaussian过程的Bayesian反褶积方法。Sachhi提出利用重加权策略实现地震反褶积。Kaaresen基于后验概率极大准则，提出迭代窗口校正函数极大化的稀疏脉冲反褶积方法。Velis提出统计稀疏脉冲反褶积方法，利用模拟退火法来确定反射层所在的位置，然后利用线性最小二乘方法估计反射系数的振幅。Gholami和Sacchi将基追踪去噪模型引入地震反褶积，提出一种从含噪地震记录中估计稀疏脉冲反射系数序列的方法。

以上现有技术具有如下缺点：

(1)常规的反射系数序列反演算法通常是基于梯度或共轭梯度构造下降方向，其收敛率是线性的，收敛速度较慢；

(2)常规的反射系数序列反演通常是逐道反演的，忽略了地震记录的横向连续性，进而使得后续叠加阻抗低频噪声严重。

发明内容

本发明的目的在于克服上述现有技术的缺点，提供一种新的用于多道地震记录反射系数序列同时反演方法。该方法基于地震记录褶积模型，充分考虑临近地震道之间的相关性，并利用近似零范数度量反射系数序列的稀疏特征，建立一种新的多道地震记录反射系数同时反演模型。在数值实现上，基于褶积模型特点和反射系数稀疏性特征，利用预处理技术降低求解维度，然后采用信赖域二次规划算法对本发明建立的反演模型进行求解，具有较高的计算效率和全局收敛性。该技术方案易于实现，可操作性强。

本发明的目的是通过以下技术方案来解决的：

该种用于多道地震记录反射系数序列同时反演方法，包括以下步骤：

1)采集原始地震资料，然后对采集到的地震资料进行预处理，包括静校正处理、噪音压制处理、真振幅恢复处理，得到叠后地震数据，称这个地震数据为观测地震数据，记为S(x,t)，其中x表示CDP号，t表示时间轴；

2)基于统计性方法提取地震子波w(t)；

3)利用快速衰减的指数函数f_σ(·)来逐渐逼近零范数，并给定初始值σ；

4)确定正则化参数μ，估计地震记录噪声能量ε²，记反射系数剖面为R，记优化目标函数为h(R,σ,μ)，记约束条件为c(R)，建立多道地震记录反射系数序列同时反演优化目标函数；

5)利用阻尼最小二乘法估计反射系数初值R⁰；

6)构造预处理降维矩阵A；

7)确定当前迭代步的零范数逼近参数σ_k，构造优化问题对应的Lagrange 函数L(R^(k),λ^(k))，计算目标函数的梯度g_k，以及Lagrange函数的海森阵的正定近似矩阵B_k；

8)利用信赖域技术求解二次规划子问题Q(d,λ)，寻求当前迭代点R^(k)的修正量d^(k)并改进相应的Lagrange乘子λ_k；

9)进行迭代k＝k+1；判断是否满足终止条件，若不满足，则返回步骤7)；否则，执行步骤10)；

10)输出最终的反射系数剖面R^*。

上述步骤3)中，利用快速衰减的指数函数来逐渐逼近零范数，该函数的表达式如下

式中，σ为高斯函数的标准差，当σ趋近于0时，f_σ(x)趋近于零范数。

上述步骤4)中，优化目标函数h(R,σ,μ)的表达式为

式中，M为地震道总道数，N为时间采样点数，i为时间离散指标，j为道离散指标，μ为正则化参数，由GCV准则确定，ε²为噪声能量水平，R为反射系数剖面，R_i,j为第j道第i个时间采样点反射系数。

约束条件c(R)的表达式为

式中，W为地震子波w(t)张成的Toepolitz矩阵，S为叠后地震记录。

进而，多道地震记录反射系数序列同时反演优化目标函数可简写为如下形式

上述步骤6)中，基于高斯随机分布构造降维矩阵A，该矩阵的行数远小于列数，然后对其进行谱范数归一化。记降维处理后的地震记录和Toeplitz 矩阵分别为

上述步骤7)中，经过降维处理后，约束条件可改写为

此时，优化问题对应的Lagrange函数L(R^(k),λ^(k))的表达式为

L(R^(k),λ^(k))＝h(R,σ_k,μ)-λc(R)

式中，k表示迭代步数。目标函数的梯度B_k由BFGS校正公式生成。

上述步骤8)中，二次规划子问题Q(d,λ)如下所示

式中，d为当前迭代点R^(k)的修正量，ρ_k为信赖域半径。

本发明具有以下有益效果：

本发明从地震记录褶积模型出发，充分考虑了临近地震道之间的空间相关性和反射系数的稀疏特征，建立了一种新的多道地震记录反射系数同时反演的约束优化目标函数。利用参数化的光滑指数函数逐渐逼近零范数，以度量反射系数的稀疏特征。基于褶积模型的特点和反射系数的稀疏性假设，在不丢失任何信息的前提下，设计降维预处理矩阵，降低约束目标函数的维度。在算法实现上，利用信赖域序列二次规划算法对优化问题进行高效求解。该算法具有全局收敛性和超线性收敛率，计算效率高。相比于传统的单道地震记录反褶积算法，本发明算法考虑了地震道的空间相关性，反演结果更稳定、可靠。

附图说明

图1是本发明流程示意图；

图2是地震子波估计示意图；

图3是某陆上实际叠后地震资料。

图4是地震子波估计结果。

图5是本发明反射系数反演剖面。

具体实施方式

下面结合附图对本发明做进一步详细描述：

反射系数序列反演方法通过压缩地震子波，进而获取地层反射系数序列，是地震资料处理中提高纵向分辨率、精细刻画地层结构常用方法之一。本发明基于反射系数稀疏性假设，充分考虑地震道的空间相关性，建立多道地震记录反射系数同时反演模型。在算法实现上，利用信赖域二次规划算法对该模型进行高效求解，反演结果稳定、可靠。

本发明的物质基础是通过野外地震采集设备采集到的大量地震数据。本发明的理论基础是地震记录褶积模型。

在一定的假设条件之下，地震记录可以看做是地震子波与地层反射系数序列褶积结果。若假设地震子波在一定时空区域内是平稳的，即地震子波与空间和时间无关，此时，离散形式的多道褶积模型可写为

S＝WR+E

式中，R∈R^N×M为反射系数序列，共M道，N个时间采样点；若地震子波的长度为l，则W∈R^L×N为地震子波构成的Toeplitz矩阵，L＝N+l-1；E∈R^L×M为加性噪声；S∈R^L×M为观测到的叠后地震记录。

本发明的多道地震记录反射系数序列同时反演框架，具体步骤分别为：

1)采集原始地震资料，然后对采集到的地震资料做常规预处理，包括静校正处理、噪音压制处理、真振幅恢复处理。处理后得到叠前共中心点道集或叠后地震数据，称这个地震数据为观测地震数据，记为S(x,t)，其中x表示CDP号，t表示时间轴。

2)基于统计性方法提取地震子波w(t)。先对叠后地震记录S(x,t)逐道做傅里叶变换，然后计算多道平均振幅谱，利用多项式拟合地震子波的振幅谱，然后利用峭度最大化准则估计地震子波的相位，结合估计得到的地震子波振幅谱和相位，通过逆傅里叶变换得到地震子波w(t)。提取地震子波的流程如图2所示。

3)构造逼近零范数的光滑函数f_σ(·)，并给定初始值σ。一个序列的零范数定义为该序列中非零元素的个数。尽管零范数是度量反射系数序列稀疏特征最有力的工具，但是零范数度量非凸非光滑，求解难度大，计算量高。为了后续迭代算法设计的需要，需要选取合适的函数对零范数进行近似，这里选用如下所示的高斯函数

式中，σ为高斯函数的标准差，当σ趋近于0时，f_σ(x)趋近于零范数。σ在后续迭代中可动态选取，这里设定其初值为0.5。

4)首先，通过广义交叉确认法，估计正则化参数μ。地震记录噪声能量ε²可通过对地震记录做小波变换分析得到。记反射系数剖面为R，为了后续叙述方便，记优化目标函数为h(R,σ,μ)，其表达式，

式中，i为时间离散指标，j为道离散指标，μ为正则化参数，由GCV准则确定，ε²为噪声能量水平，R_i,j为第j道第i个时间采样点反射系数。

类似地，为了后续叙述方便，记约束条件为c(R)，其表达式为

最终，本发明建立多道地震记录反射系数序列同时反演约束优化目标函数如下所示

5)利用阻尼最小二乘法估计反射系数R⁰，作为后续迭代的初值。阻尼最小二乘法估计公式为

R⁰＝(W^TW+δI)^-1W^TS

式中，T为矩阵转置符号，I为单位矩阵，δ为阻尼参数，为了保证初值估计的稳定性，可取为0.05。

6)构造预处理降维矩阵A。由于褶积模型是一个超定矩阵，即方程个数大于未知量的维度，有必要对褶积模型进行降维处理。此外，由于反射系数序列满足稀疏性假设，在不丢失任何信息的前提下，可对褶积模型 S＝WR+E两边同时乘以预处理矩阵A，其中，A∈R^K×N是行满秩的，其维度满足K＜N＜L，且||A||₂＝1。记降维处理后的地震记录和Toeplitz矩阵分别为则褶积模型的维度得到有效降低。此外，由于||A||₂＝1，则||AE||₂≤||A||₂||E||₂≤||E||₂。

7)确定当前迭代步的零范数逼近参数σ_k，构造优化问题对应的Lagrange 函数L(R^(k),λ^(k))，计算目标函数的梯度g_k，以及Lagrange函数的海森阵的正定近似矩阵B_k。经过降维处理后，约束条件可改写为

此时，优化问题对应的Lagrange函数L(R^(k),λ^(k))的表达式为

L(R^(k),λ^(k))＝h(R,σ_k,μ)-λc(R)

式中，k表示迭代步数。目标函数的梯度B_k由BFGS校正公式生成

式中，y_k＝g_k-g_k-1，η_k＝R^(k)-R^(k-1)，trc(·)表示矩阵的迹。当k＝1时，B₁＝I 为单位矩阵。

8)利用信赖域技术求解二次规划子问题Q(d,λ)，寻求当前迭代点R^(k)的修正量d^(k)并改进相应的Lagrange乘子λ_k；

二次规划子问题Q(d,λ)如下所示

式中，d为当前迭代点R^(k)的修正量，ρ_k为信赖域半径。利用预处理共轭梯度法求取信赖域子问题Q(d,λ)，获取最优的修正量d^(k)。

9)进行迭代R^(k+1)＝R^(k)+d^(k)，k＝k+1；判断是否满足终止条件

式中，ξ为误差容许阈值，可取为0.0001。

若不满足，则返回步骤7)；否则，R^*＝R^(k+1)，执行步骤10)；

10)输出最终的反射系数剖面R^*。

实际地震资料测试

本节利用实际地震资料来验证本发明提出的多道地震记录反射系数同时反演算法的有效性。图3为某陆上实际叠后地震资料中的一条测线。该测线共1751 道，道间距20m，截取时间范围为1400～2200ms，时间采样间隔为2ms。测井位于CDP 1020处。可见，地震资料中采集脚印干扰严重，信噪比较差。利用图2给出流程估计地震子波，估计结果如图4所示。图5为本发明提出的多道地震记录反射系数同时反演算法反演得到的反射系数剖面。可见，反射系数剖面结构清晰、空间连续性好，且具有较高的纵向分辨率，实现对地质结构得到精细刻画。

Claims (6)

1.一种多道地震记录反射系数序列同时反演方法，其特征在于，包括以下步骤：

1)采集原始地震资料，对采集到的地震资料进行预处理，包括静校正处理、噪音压制处理、真振幅恢复处理，得到叠后地震数据，称这个地震数据为观测地震数据，记为S(x,t)，其中x表示CDP号，t表示时间轴；

2)基于统计性方法提取地震子波w(t)；

3)利用快速衰减的指数函数f_σ(·)来逐渐逼近零范数，并给定初始值σ；

4)确定正则化参数μ，估计地震记录噪声能量ε²，记反射系数剖面为R，记优化目标函数为h(R,σ,μ)，记约束条件为c(R)，建立多道地震记录反射系数序列同时反演的优化目标函数；

5)利用阻尼最小二乘法估计反射系数初值R⁰；

6)构造预处理降维矩阵A；

7)确定当前迭代步的零范数逼近参数σ_k，构造优化问题对应的Lagrange函数L(R^(k),λ^(k))，计算目标函数的梯度g_k，以及Lagrange函数的海森阵的正定近似矩阵B_k；

8)利用信赖域技术求解二次规划子问题Q(d,λ)，寻求当前迭代点R^(k)的修正量d^(k)并改进相应的Lagrange乘子λ_k；

9)进行迭代k＝k+1；判断是否满足终止条件，若不满足，则返回步骤7)；否则，执行步骤10)；

10)输出最终的反射系数剖面R^*。

2.根据权利要求1所述的一种多道地震记录反射系数序列同时反演方法，其特征在于，所述步骤3)中，利用快速衰减的指数函数来逐渐逼近零范数，该函数的表达式如下

式中，σ为高斯函数的标准差，当σ趋近于0时，f_σ(x)趋近于零范数，σ的初值可取较大数值，可取值为0.5，以保证目标函数的光滑性。

3.根据权利要求1所述的一种多道地震记录反射系数序列同时反演方法，其特征在于，上述步骤4)中，优化目标函数h(R,σ,μ)的表达式为

式中，M为地震道总道数，N为时间采样点数，i为时间离散指标，j为道离散指标，μ为正则化参数，由GCV准则确定，ε²为噪声能量水平，R为反射系数剖面，R_i,j为第j道第i个时间采样点反射系数；

约束条件c(R)的表达式为

式中，W为地震子波w(t)张成的Toepolitz矩阵，S为叠后地震记录；

进而，多道地震记录反射系数序列同时反演优化目标函数可简写为如下形式

4.根据权利要求1所述的一种多道地震记录反射系数序列同时反演方法，其特征在于，上述步骤6)中，基于高斯随机分布构造降维矩阵A，该矩阵的行数远小于列数，然后对其进行谱范数归一化，记降维处理后的地震记录和Toeplitz矩阵分别为其中，W为地震子波w(t)张成的Toepolitz矩阵，S为叠后地震记录。

5.根据权利要求1所述的一种多道地震记录反射系数序列同时反演方法，其特征在于，上述步骤7)中，经过降维处理后，约束条件可改写为

记降维处理后的地震记录，为降维后Toeplitz矩阵，此时，优化问题对应的Lagrange函数L(R^(k),λ(^k))的表达式为

L(R^(k),λ^(k))＝h(R,σ_k,μ)-λc(R)

式中，k表示迭代步数，目标函数的梯度B_k由BFGS校正公式生成。

6.根据权利要求1所述的一种多道地震记录反射系数序列同时反演方法，其特征在于，上述步骤8)中，二次规划子问题Q(d,λ)如下所示

式中，d为当前迭代点R(^k)的修正量，ρ_k为信赖域半径。