CN115079264A - 一种不依赖子波的粘声最小二乘逆时偏移成像方法 - Google Patents

Abstract

本发明属于油藏开发技术领域，具体涉及一种不依赖子波粘声最小二乘逆时偏移成像方法，在方法中采用常分数阶解耦拉普拉斯粘声波动方程模拟地震波在粘声介质中的传播，建立卷积目标函数，构建具有相同震源的模拟数据和观测数据间的匹配，从而降低对震源子波的敏感性，推导正演算子与伴随算子，在迭代过程中补偿衰减效应，实现不依赖子波的粘声最小二乘逆时偏移，在震源子波不准确的情况下，获取可靠的能够补偿深层能量的成像结果。本方法更能满足油气勘探开发对复杂地下构造成像的需求，提高对油气藏的识别精度，提升我国油气勘探开发核心竞争力。

Description

一种不依赖子波的粘声最小二乘逆时偏移成像方法

技术领域

本发明属于油藏开发技术领域，具体涉及一种不依赖子波的粘声最小二乘逆时偏移成像方法。

背景技术

油气资源是现代工业的血液，是生存发展不可或缺的战略资源之一。地震勘探通过人工激发的地震波来定位地下油气藏，是寻找油气资源的重要技术。做为地震勘探技术的重要环节之一，地震数据成像会影响对地下构造识别的精度。近年来，随着油气勘探领域的不断深入，勘探对象日趋复杂化，由常规油藏转变为非常规油气、复杂致密介质，中浅层勘探过渡为深层勘探，现有的地震成像技术难以满足油气勘探开发的需求，需要地震成像精度大幅度提升。

实际上，地下介质具有粘滞性，包含衰减和频散效应，致使地震波主频降低，振幅减小，尤其在高衰减区域，衰减现象尤为明显，振幅衰减与相位频散积累，会降低地震成像的分辨率，令地下构造定位不准确，影响地震解释的可靠性。为获取正确的地下构造成像结果，应在逆时偏移成像计算过程中考虑地层吸收衰减作用，并进行补偿。粘性介质吸收衰减补偿方法可以在地震波传播过程中进行振幅与相位的补偿，例如Q补偿逆时偏移(Q-RTM)，然而Q-RTM在振幅补偿时，会导致波场记录中的高频信息，尤其是高频噪声被放大，引起数值不稳定。针对该现象，可采用低通滤波器压制高频噪音，Suh et al.(2012)提出应用高切滤波器控制高频成分的放大。Sun和Zhu(2015)提出一种稳定衰减补偿方法，通过求取速度频散波场和粘声波场的比值，构造稳定的补偿算子，然而，此方法计算成本大约是传统RTM的1.5–2倍，而且还有除以零的不稳定性问题；田坤等(2017)对非相位频散粘声方程添加正则化项，实现波场稳定传播；Wang et al.(2018)基于稳定反Q滤波，提出一种具有时变性和Q依赖性的自适应稳定Q-RTM方法；Yang et al.(2021)提出一种基于具有鲁棒性的空间波数域衰减补偿算子，该算子在波场传播过程中能够自动考虑与波数相关的衰减效应并进行校正，可用于叠前与叠后粘声逆时偏移成像。然而，这些方法在成像过程中均存在高频噪音的影响，需要引入滤波器或修改方程进行压制。针对此现象，陈汉明等(2020)提出一种基于常分数阶拉普拉斯算子粘声波方程的最小二乘逆时偏移，该方法正演和伴随算子均呈现能量衰减特征，不存在数值不稳定问题，且最小二乘逆时偏移可以提升成像分辨率。

基于反演理论，最小二乘逆时偏移通过最大限度地缩小模拟数据和观测数据之间的差距来更新反射率成像，能够有效压制偏移串扰，增加照明，提高空间分辨率。研究人员对最小二乘展开广泛研究，Dutta和Schuster(2014)提出了一种时域粘声波最小二乘逆时偏移成像方法，该方法从含记忆变量的SLS粘声波方程出发，有效补偿深层能量，改善成像剖面质量；Yang et al.(2018)提出一种基于时域复数波动方程的粘声最小二乘逆时偏移方法，引入TV正则化增强鲁棒性，并使用对角Hessian作为预条件子来加速收敛；Yang etal.(2022)设计空间波数滤波器来逼近高斯-牛顿Hessian(GNH)，并将其用作最小二乘梯度的预条件算子，以加速收敛速率。

传统基于L2范数的粘声最小二乘逆时偏移(VLSRTM)虽然可以在迭代过程中对衰减效应进行有效补偿，但它对子波具有一定的依赖性，鉴于实际数据中很难得到精确的子波，该方法通常无法获得令人满意的补偿成像效果。

为了满足油气工业对地下粘滞性介质准确成像的需求，急需一种可以解决震源子波不准确导致成像结果精度降低的方法。基于此，本专利提出了一种不依赖子波粘声最小二乘逆时偏移成像方法。

发明内容

为了解决上述技术问题，提出了一种不依赖子波粘声最小二乘逆时偏移成像方法，在方法中采用常分数阶解耦拉普拉斯粘声波动方程模拟地震波在粘声介质中的传播，建立卷积目标函数，构建具有相同震源的模拟数据和观测数据间的匹配，从而降低对震源子波的敏感性，推导正演算子与伴随算子，在迭代过程中补偿衰减效应，实现不依赖子波的粘声最小二乘逆时偏移，在震源子波不准确的情况下，获取可靠的能够补偿深层能量的成像结果。本方法更能满足油气勘探开发对复杂地下构造成像的需求，提高对油气藏的识别精度，提升我国油气勘探开发核心竞争力。

本发明采用的技术方案为：一种不依赖子波的粘声最小二乘逆时偏移成像方法，主要包含两个重点计算部分，其一为正演算子与伴随算子的推导，其二为卷积目标函数的构建：

首先，常分数阶拉普拉斯粘滞声波方程为：

式中：Q为品质因子，

为拉普拉斯算子，c₀为地震波速度，ω₀为参考角频率，ω_d＝2πf_d为主角频率，f_d为子波主频，s(t；x_s)为震源子波，x_s为震源位置，该方程能够解耦相位频散与振幅衰减效应，并将变分数阶拉普拉斯算子合理的近似为不受混合域影响的常分数阶算子，在Q≥15的情况下，能够精确近似变分数阶拉普拉斯粘声波方程；

根据扰动理论，将模型参数分解为背景模型c₀和扰动模型δc两部分，分别对应背景波场u₀和扰动波场δu，采用泰勒公式展开：

代入公式(1)，可得：

并与公式(1)相减，忽略高阶项，可得Born正演方程：

Aδu＝mDu₀ (5)

其中，

为反射率；

常规的粘声最小二乘逆时偏移基于L2范数建立目标函数：

式中，U_syn为模拟数据，U_obs为观测数据，根据Tarantola(1984)，地震数据可表达为：

U＝G*_s (7)

式中，*为卷积运算，U为地震数据，G为格林函数，将式(7)代入式(6)，可得

式中，s_syn为模拟数据子波，s_obs为观测数据子波，当两者不一致时，会影响成像结果的计算；

为降低子波对成像结果的影响，建立卷积型目标函数：

式中，

和

分别为模拟数据参考道和观测数据参考道，当地震数据信噪比较高时，选择距震源较近的地震道做为参考道，反之，则选择平均道以压制噪音，将式(7)代入式(9)可得

式中，

s'＝s_obs*s_syn，因此，式(10)看作是具有相同震源s'的模拟数据与观测数据间的匹配，能够降低真实子波与偏移子波不一致带来的影响；

利用拉格朗日乘子法求解约束最优化问题，

式中，H为计算区域，μ为伴随波场，令

可得伴随方程

式中，

代表互相关算子，由于参考道中直达波后的反射经卷积与互相关计算会引入噪音，因此需要对参考道应用时窗，保留直达波部分，所以式(12)可以重写为：

式中，W为时窗，

式中，t_c为时窗长度，n为时窗阶数；

进而得到梯度：

采用共轭梯度法迭代更新反射率成像，公式如下：

式中，g_k代表梯度，d_k代表共轭梯度方向，β_k满足：

反射率成像迭代更新：

m_k+1＝m_k+α_kd_k (19)

式中，k为迭代次数，α_k为迭代更新的步长。

本发明的有益效果：提供了一种不依赖子波粘声最小二乘逆时偏移成像方法，在方法中采用常分数阶解耦拉普拉斯粘声波动方程模拟地震波在粘声介质中的传播，建立卷积目标函数，构建具有相同震源的模拟数据和观测数据间的匹配，从而降低对震源子波的敏感性，推导正演算子与伴随算子，在迭代过程中补偿衰减效应，实现不依赖子波的粘声最小二乘逆时偏移，在震源子波不准确的情况下，获取可靠的能够补偿深层能量的成像结果。本方法更能满足油气勘探开发对复杂地下构造成像的需求，提高对油气藏的识别精度，提升我国油气勘探开发核心竞争力。其主要优点如下：

(1)、本方法正演算子和伴随算子均呈现衰减特征，避免数值不稳定现象，可在迭代过程中实现衰减效应补偿，提高成像结果分辨率。

(2)、基于卷积型目标函数，构建具有相同震源子波的模拟数据与观测数据间的匹配，解决模拟数据子波与观测数据子波不一致带来的影响。

(3)、卷积目标函数具有一定的抗噪性，在观测地震数据含噪时，本方法可以恢复较为清晰的成像结果。

附图说明

图1为不依赖子波粘声最小二乘逆时偏移成像的流程图；

图2为测试一中的纵波速度模型图；

图3为测试一中的Q模型图；

图4为测试一中的扰动模型图；

图5为测试一中ricker子波图；

图6为测试一中高斯一阶导子波图；

图7为测试一中ricker子波向右平移1/f的子波图；

图8为测试一中ricker子波向右平移2/f的子波图；

图9为测试一中子波正确的常规VLSRTM成像图；

图10为测试一中子波正确的不依赖子波VLSRTM成像图；

图11为测试一中子波错误的常规VLSRTM成像图；

图12为测试一中子波错误的不依赖子波VLSRTM成像图；

图13为测试一中子波正确的常规VLSRTM单道信息图；

图14为测试一中子波正确的不依赖子波VLSRTM单道信息图；

图15为测试一中子波错误的常规VLSRTM单道信息图；

图16为测试一中子波错误的不依赖子波VLSRTM单道信息图；

图17为测试一中ricker子波向右平移1/f常规VLSRTM成像图；

图18为测试一中ricker子波向右平移1/f不依赖子波VLSRTM成像图；

图19为测试一中ricker子波向右平移2/f常规VLSRTM成像图；

图20为测试一中ricker子波向右平移2/f不依赖子波VLSRTM成像图；

图21为测试一中ricker子波向右平移1/f常规VLSRTM单道信息图；

图22为测试一中ricker子波向右平移1/f不依赖子波VLSRTM单道信息图；

图23为测试一中ricker子波向右平移2/f常规VLSRTM单道信息图；

图24为测试一中ricker子波向右平移2/f不依赖子波VLSRTM单道信息图；

图25为测试二中纵波速度模型图；

图26为测试二中Q模型图；

图27为测试二中扰动模型图；

图28为测试二中ricker子波图；

图29为测试二中高斯一阶导子波图；

图30为测试二中ricker子波向左平移1/f的子波图；

图31为测试二中声波地震记录图；

图32为测试二中粘声波地震记录图；

图33为测试二中单道振幅对比图；

图34为测试二中频谱对比图；

图35为测试二中声波LSRTM成像图；

图36为测试二中基于L2范数的VLSRTM成像图；

图37为测试二中不依赖子波VLSRTM成像图；

图38为测试二中声波LSRTM单道信息图；

图39为测试二中基于L2范数的VLSRTM单道信息图；

图40为测试二中不依赖子波VLSRTM单道信息图；

图41为测试二中声波LSRTM的波数谱图；

图42为测试二中基于L2范数的VLSRTM的波数谱图；

图43为测试二中不依赖子波VLSRTM的波数谱图；

图44为测试三中子波错误的常规VLSRTM成像图；

图45为测试三中子波错误的不依赖子波VLSRTM成像图；

图46为测试三中子波错误的常规VLSRTM单道信息图；

图47为测试三中子波错误的不依赖子波VLSRTM单道信息图；

图48为测试三中ricker子波向左平移的常规VLSRTM成像图；

图49为测试三中ricker子波向左平移的不依赖子波VLSRTM成像图；

图50为测试三中ricker子波向左平移的常规VLSRTM单道信息图；

图51为测试三中ricker子波向左平移的不依赖子波VLSRTM单道信息图；

图52为测试四中原始观测地震数据图；

图53为测试四中含噪观测地震数据图；

图54为测试四中基于含噪观测数据的常规VLSRTM成像图；

图55为测试四中基于含噪观测数据的不依赖子波VLSRTM成像图；

图56为图36的局部放大图；

图57为图37的局部放大图；

图58为图54的局部放大图；

图59为图55的局部放大图。

具体实施方式

实施例一

一种不依赖子波的粘声最小二乘逆时偏移成像方法，主要包含两个重点计算部分，其一为正演算子与伴随算子的推导，其二为卷积目标函数的构建：

首先，常分数阶拉普拉斯粘滞声波方程为：

式中：Q为品质因子，

为拉普拉斯算子，c₀为地震波速度，ω₀为参考角频率，ω_d＝2πf_d为主角频率，f_d为子波主频，s(t；x_s)为震源子波，x_s为震源位置，该方程能够解耦相位频散与振幅衰减效应，并将变分数阶拉普拉斯算子合理的近似为不受混合域影响的常分数阶算子，在Q≥15的情况下，能够精确近似变分数阶拉普拉斯粘声波方程；

根据扰动理论，将模型参数分解为背景模型c₀和扰动模型δc两部分，分别对应背景波场u₀和扰动波场δu，采用泰勒公式展开：

代入公式(1)，可得：

并与公式(1)相减，忽略高阶项，可得Born正演方程：

Aδu＝mDu₀ (5)

其中，

为反射率；

常规的粘声最小二乘逆时偏移基于L2范数建立目标函数：

式中，U_syn为模拟数据，U_obs为观测数据，根据Tarantola(1984)，地震数据可表达为：

U＝G*s (7)

式中，*为卷积运算，U为地震数据，G为格林函数，将式(7)代入式(6)，可得

式中，s_syn为模拟数据子波，s_obs为观测数据子波，当两者不一致时，会影响成像结果的计算；

为降低子波对成像结果的影响，建立卷积型目标函数：

式中，

和

分别为模拟数据参考道和观测数据参考道，当地震数据信噪比较高时，选择距震源较近的地震道做为参考道，反之，则选择平均道以压制噪音，将式(7)代入式(9)可得

式中，

s'＝s_obs*s_syn，因此，式(10)可看作是具有相同震源s'的模拟数据与观测数据间的匹配，能够降低真实子波与偏移子波不一致带来的影响；

利用拉格朗日乘子法求解约束最优化问题，

式中，H为计算区域，μ为伴随波场，令

可得伴随方程

式中，

代表互相关算子，由于参考道中直达波后的反射经卷积与互相关计算会引入噪音，因此需要对参考道应用时窗，保留直达波部分，所以式(12)可以重写为：

式中，W为时窗，

式中，t_c为时窗长度，n为时窗阶数；

进而得到梯度：

采用共轭梯度法迭代更新反射率成像，公式如下：

式中，g_k代表梯度，d_k代表共轭梯度方向，β_k满足：

反射率成像迭代更新：

m_k+1＝m_k+α_kd_k (19)

式中，k为迭代次数，α_k为迭代更新的步长。

实施例二

测试一：

图2-图4分别为纵波速度模型、Q模型和扰动模型图，如图5-图8分别为ricker子波图、高斯一阶导子波图、ricker子波向右平移1/f的子波图和ricker子波向右平移2/f的子波图。在第一个数值算例中，首先对层状模型进行测试来验证所提出方法在子波不准确时的有效性；所使用的层状介质模型大小为301×301个网格点，空间采样间隔为10m，时间采样率为1.0ms，采样点数为2500；共设置15炮震源激发，均匀分布于地表，炮间距为200m，每炮布置201个检波器接收，均匀分布于炮点两侧；采用主频为25Hz的ricker子波做为正确子波，如图5所示，将扰动模型做为参考解。

为测试子波波形对成像结果的影响，采用主频为25Hz的高斯一阶导作为偏移子波，如图6所示。图9和图10为子波正确的常规VLSRTM成像结果与不依赖子波的VLSRTM成像结果，图11和图12为子波错误的常规VLSRTM成像结果与不依赖子波VLSRTM成像结果。对比扰动解，可以看出子波正确的VLSRTM、不依赖子波VLSRTM以及子波错误的不依赖子波VLSRTM均可清晰成像，反射界面正确归位，与扰动解相位一致，而子波错误的常规VLSRTM的相位与扰动解有明显差异，按图9白色虚线所在位置抽取水平方向1.5km处单道信息，如图13-16所示，图中，reference为扰动解，LSRTM为粘声最小二乘逆时偏移成像，可以看出子波错误的常规VLSRTM结果与扰动解明显不匹配。

进一步测试子波相位对成像结果的影响，对25Hz的ricker子波分别向右平移1/f、2/f做为偏移子波，f为子波主频，如图7和图8所示。如图17和图19为常规VLSRTM成像，图18和图20为不依赖子波VLSRTM成像结果；如图17-图20中白色虚线框内所示，可以看出随着子波相位向右平移，常规VLSRTM成像结果反射界面向上移动，而不依赖子波VLSRTM成像结果反射界面位置不变，抽取单道信息如图21-图24所示，对比可知，随着子波相位平移幅度增大，常规VLSRTM结果相位与振幅均与扰动解逐渐偏离，无法准确识别构造位置，而不依赖子波VLSRTM成像在子波错误的情况下仍可获取正确的成像结果。

测试二：

对复杂Marmousi模型进行测试，验证所提出方法对衰减效应的补偿效果。如图25所示为速度模型，该模型尺寸为767×330个网格点，Q模型、扰动模型分别如图26和图27所示，空间网格间距为10m，时间采样率为1.0ms，记录时长为4s。共设置20炮震源激发，均匀分布于地表，炮间距为380m，每炮设置767个检波器接收。采用如图28所示的主频为25Hz的ricker子波做为正确子波生成观测数据，并将扰动模型做为参考解。

如图31为声波地震数据，图32为粘声波地震数据，分别抽取第385道数据，输出如图33的单道振幅与如图34的频谱，图中，ac为声波地震数据，vis为粘声波地震数据，可以看出，相比于声波地震数据，粘声波地震数据深层振幅与高频成分均发生衰减现象。采用粘声地震数据进行成像，如图35为声波LSRTM成像结果，图36和图37分别为基于L2范数的VLSRTM和不依赖子波的VLSRTM成像结果，对比扰动解，可以看出声波LSRTM由于忽略地层粘滞性，浅层噪音严重，深层能量弱，而两种VLSRTM成像能够清晰刻画复杂地下地质构造，按图35中白色曲线所在位置抽取水平方向4.0km处单道信息，如图38-40所示，并输出三种LSRTM成像结果的波数谱，如图41-43所示，可以看出基于L2范数的VLSRTM和不依赖子波的VLSRTM均可有效补偿衰减效应，增强深层能量，提高成像分辨率。

测试三：

对Marmousi模型进行测试，验证所提出方法在地下构造复杂时对子波的不敏感性。模型参数与测试二相同，同样采用如图28所示的主频为25Hz的ricker子波做为正确子波生成观测数据。

测试子波波形对成像结果的影响，采用如图29所示的主频25Hz高斯一阶导做为偏移子波。图44和图45为在子波错误情况下的常规VLSRTM成像结果与不依赖子波VLSRTM成像结果。对比扰动解，可以看出子波错误的不依赖子波VLSRTM空间归位正确，背斜、不整合面等地质结构清晰，能够对复杂构造准确成像，而子波错误的常规VLSRTM相位与扰动解有明显差异。按图35中白色虚线所在位置抽取单道信息，如图46和图47所示，可以看出子波错误的常规VLSRTM成像相位、振幅与扰动解不匹配，并且如箭头指示，存在层位错动问题，而不依赖子波VLSRTM与扰动解较为一致。

进一步测试子波相位对成像结果的影响，将图28中的ricker子波向左平移做为偏移子波，如图30所示。图48和图49分别为子波错误的常规VLSRTM成像和子波错误的不依赖子波VLSRTM成像。对比扰动解，可以看出当子波相位向左平移时，常规VLSRTM成像结果模糊，且如图48和图49中白色虚线框内所示，层位向下移动，而不依赖子波VLSRTM成像结果与扰动解较为匹配。按图35中白色虚线所在位置抽取单道信息，如图50和图51所示，对比可知，不依赖子波VLSRTM成像可降低对震源子波的敏感性。

测试四：

对Marmousi模型进行测试，验证所提出方法的抗噪性。模型参数与测试二相同，对地震数据加入高斯随机噪音进行测试，信噪比S/N＝5，图52和图53分别为无噪与含噪的第10炮地震数据，采用正确子波进行成像。图54和图55为采用含噪观测数据进行计算的常规VLSRTM成像与不依赖子波VLSRTM成像，两者均可刻画主要地质结构，但成像剖面质量降低，为清晰对比，按图54中白色虚线框所在位置抽取无噪数据成像(图36和图37)与含噪数据成像(图54和图55)的部分成像结果，如图56-图59所示。经对比，可以看出，不依赖子波VLSRTM可以恢复更清晰的成像结果，这是因为卷积目标函数作为带通滤波器能够抑制一些频带超出子波频带的噪声。

不依赖子波的粘声最小二乘逆时偏移成像方法能够有效补偿地层吸收衰减作用，提高深层能量，正演算子和伴随算子均为能量衰减型，可以避免数值不稳定现象。采用卷积目标函数，可以降低对子波的敏感性。与现有技术相比具有的优点：比如产率、质量、精度和效率的提高，能耗、原材料、工序的节省；加工、操作、控制、使用的简便，环境污染的治理或根治等。

Claims (1)

1.一种不依赖子波的粘声最小二乘逆时偏移成像方法，其特征在于：不依赖子波的粘声最小二乘逆时偏移成像方法为：

首先，常分数阶拉普拉斯粘滞声波方程为：

式中：Q为品质因子，

为拉普拉斯算子，c₀为地震波速度，ω₀为参考角频率，ω_d＝2πf_d为主角频率，f_d为子波主频，s(t；x_s)为震源子波，x_s为震源位置，该方程能够解耦相位频散与振幅衰减效应，并将变分数阶拉普拉斯算子合理的近似为不受混合域影响的常分数阶算子，在Q≥15的情况下，能够精确近似变分数阶拉普拉斯粘声波方程；

根据扰动理论，将模型参数分解为背景模型c₀和扰动模型δc两部分，分别对应背景波场u₀和扰动波场δu，采用泰勒公式展开：

代入公式(1)，可得：

并与公式(1)相减，可得Born正演方程：

Aδu＝mDu₀ (5)

其中，

为反射率；

常规的粘声最小二乘逆时偏移基于L2范数建立目标函数：

式中，U_syn为模拟数据，U_obs为观测数据，地震数据可表达为：

U＝G*s (7)

式中，*为卷积运算，U为地震数据，G为格林函数，将式(7)代入式(6)，可得

式中，s_syn为模拟数据子波，s_obs为观测数据子波，当两者不一致时，会影响成像结果的计算；

为降低子波对成像结果的影响，建立卷积型目标函数：

式中，

和

分别为模拟数据参考道和观测数据参考道，将式(7)代入式(9)可得

式中，

s'＝s_obs*s_syn，因此，式(10)看作是具有相同震源s'的模拟数据与观测数据间的匹配，能够降低真实子波与偏移子波不一致带来的影响；

利用拉格朗日乘子法求解约束最优化问题，

式中，H为计算区域，μ为伴随波场，令

可得伴随方程

式中，

代表互相关算子，由于参考道中直达波后的反射经卷积与互相关计算会引入噪音，因此需要对参考道应用时窗，保留直达波部分，所以式(12)可以重写为：

式中，W为时窗，

式中，t_c为时窗长度，n为时窗阶数；

进而得到梯度：

采用共轭梯度法迭代更新反射率成像，公式如下：

式中，g_k代表梯度，d_k代表共轭梯度方向，β_k满足：

反射率成像迭代更新：

m_k+1＝m_k+α_kd_k (19)

式中，k为迭代次数，α_k为迭代更新的步长。

Images