CN109738952A - 基于全波形反演驱动的被动源直接偏移成像方法 - Google Patents

Abstract

本发明涉及一种基于全波形反演驱动的被动源直接偏移成像方法。首先采用全波形反演，将地震干涉法得到的被动源虚拟炮集作为观测记录，之后采用卷积目标函数来消除全波形反演过程中存在的子波依赖问题。针对被动源数据的特点，对不依赖震源子波被动源全波形反演的震源编码方法，参考道选择方法以及梯度加权等方法做了改进；之后利用被动源全波形反演得到的结果作为被动源偏移的速度模型，采用直接成像条件的进行逆时偏移。本方法对速度模型精度依赖小，在一个低精度的速度模型下，对噪声被动源数据成像，得到了高质量的偏移结果，不需要主动源数据的参与。

Description

基于全波形反演驱动的被动源直接偏移成像方法

技术领域：

本发明涉及一种地震勘探中的被动源直接偏移成像方法。利用不依赖子波 的震源编码全波形反演求得速度模型，之后利用直接成像的边界条件进行被动 源逆时偏移。

背景技术：

地震勘探数据可以分为主动源数据和被动源数据，对于日常勘探生产中则 主要以主动源数据为主，被动源数据则视为背景噪声而被去除。主动源勘探工 程量大，成本高，不适合大范围勘察，被动源数据不需要人工震源，仅需要放 置检波器接收即可，而背景噪声中往往携带有大量的地下地质体信息，特别是 震源位置在底层深部地区，携带大量低频信息和深部位置信息，如能利用好被 动源数据，对地下信息位置工区进行预勘探，或者辅助主动源数据，都是一个 成本低廉的勘探手段。

有价值的数据信号可能随机分布在记录时间中以及一些如井下噪声、地表 环境噪声、潮汐、固体潮等噪声之内，由于不能预测这些噪声发生的时间，震 源位置未知，并且被动源震动不是连续发生的，如何合理利用被动源数据，是 提取出地下地质体信息的关键。目前常规的被动源方法都是基于地震干涉法的 数据重构或者数据提取，由于被动源重构数据低信噪比，仅基于被动源的反演 或者成像没有得到一个高精度的结果，目前大多数被动源数据的研究局限于利 用被动源数据提高主动源反演或者成像质量。

地震干涉法中的互相关重构法，通过接受到的各个检波道之间互相关，重 构出以其中一个检波点为震源，其余检波点为接收点的虚拟炮集。Artman (2006,Imagingpassive seismic data.Geophysics,71(4):SI177– SI187)提出了直接成像的方法，直接采用噪声源数据本身作为边界条件，在 频率域采用单程波偏移，得到了地下地质体构造图像，对于地震干涉法噪声源 成像，免去了重构虚拟炮集时产生的二次噪声，提高了成像的准确性，但偏移 过程依赖于真实的速度模型，且成像质量一般。

对于偏移过程需要地下的速度模型，可以通过被动源全波形反演(FWI， Fullwaveforminversion)提供。被动源中含有的地下地质体信息，可以对 工区地下构造模型进行矫正，常规主动源全波形反演和逆时偏移中，采用观测 数据和模拟数据构造目标函数，用最小二乘迭代数据模型，取得了不错的效 果。但是被动源生成的虚拟炮集中含有大量的噪声，子波也难以估计，直接使 用难以构造目标函数。Choi和Alkhalifah(2011.Source-independent time-domain waveforminversion using convolved wavefields:Application to the encoded multisource waveform inversion. Geophysics,76(5):R125-R134)提出了一种不依赖子波的时间域全波形反 演，消除了子波估计对反演的影响。不依赖子波的方法不仅可以消除错误子波 估计的影响，还可以降低被动源虚拟炮集中噪声的影响。被动源数据由于信噪 比低，需要多次迭代才能得到相应结果，导致计算效率低下。

对于被动源信号低信噪比的特点，全波形反演需要多次迭代以提高结果的 准确性，震源编码是一种提高计算效率的方法，Kerbs等(2009.Fast full- wavefield seismicinversion using encoded sources.Geophysics,74(6): WCC177-WCC188)将震源编码技术用于全波形反演，对单炮记录采用+1 和-1进行随机编码，压制了串扰噪声，Moghaddam等(2013.A new optimization approach for source-encoding full-waveforminversion.Geophysics,78(3):R125-R132)发现基于动态随机编码的全波形反演，其 梯度是无偏的随机估计，采用指数加权的方式，提高了收敛速度。由于重构后 的被动源数据信噪比低，常规的震源编码方法可能导致被动源全波形反演收敛 速度慢，反演精度不够等问题。

发明内容：

本发明的目的是针对上述现有技术的不足，尤其是针对被动源数据的特性 对FWI进行了改进，并采用直接成像的边界条件，利用FWI求得的速度进行 被动源直接逆时偏移，提供一种基于全波形反演驱动的被动源直接偏移成像方 法。

本发明的目的是通过以下技术方案实现的：

发明思想：首先采用地震干涉法重构虚拟炮集，作为全波形反演的观测记 录，采用褶积构造的目标函数，消除子波的影响，针对被动源数据的特性进行 多尺度全波形反演，改进了震源编码方法以及参考道选择方法，提高了反演效 率；之后将直接成像的边界条件应用于时间域被动源逆时偏移，排除常规逆时 偏移子波选取的干扰，且对速度依赖程度低，以被动源全波形反演得到的速度 模型为基础，采用被动源数据得到浅层正确的偏移结果。

基于全波形反演驱动的被动源直接偏移成像方法是通过MATLAB平台实 现的。

基于全波形反演驱动的被动源直接偏移成像方法，包括以下步骤：

a、地表采集获得噪声被动源数据；

b、利用地震干涉法，重构采集的噪声数据，重构为虚拟炮集：

C(x₂,x₁,t)＝T(x₁,x_S,t)*T(x₂,x_S,-t)

式中，T(x₁,x_S,t)为检波器x1接收到的传动波场数据，C(x₂,x₁,t)为通过互 相关干涉得到的序列；

c、互相关后的序列中包含因果和非因果部分，将因果非因果部分加和得到 虚拟震源炮集,由格林函数法可得：

u_vir(x₂,x₁,t)＝G(x₂,x₁,t)*S(t)

S(t)为虚拟震源的子波；

d、把重构的单炮记录做频谱分析，确定反演子波范围；

e、将得到的虚拟炮集采用全新的方式进行编码，编码方式如下所示：

为第i个编码炮，编码炮数量为nc。将各单炮数据采用震源位置随 机，得到一组随机序列，依次从随机序列中取m个值，得到nc个序列s，每个 序列s即为每个编码炮的震源位置编码，公式中p是随机生成的+1和-1序列， 长度为ns；

f、从构成该编码炮的单炮中随机选择一炮，令其炮点位置为参考道，正演 模拟记录与观测记录采用相同的参考道位置；

g、将编码炮作为观测数据，采用褶积的方法构造目标函数：

通过此目标函数消除了虚拟震源子波的影响；

h、同过伴随源法求取梯度，褶积构造的目标函数其伴随源为：

r'是第一互相关残差，作为除参考道之外的道集伴随震源，r”是第二互相 关残差，作为参考道处的伴随震源，梯度计算公式为：

r为将第一第二互相关残差合并后的总残差,v是迭代的速度模型，ns为震 源数，nr为检波器数，t0为记录时间；

i、梯度通过改进的指数加权算法求得：

这里α_k为线性搜索后的步长，α_k∈(0,1)，a为基于信噪比的收敛因子，k 是参与加权的迭代数。且这里的步长是通过每一编码炮残差求得的步长取算数 平均后得到的；

j、重复迭代步骤e到步骤i，迭代速度模型:

v_k+1＝v_k-α_kg_k

v_k为第k次迭代的速度模型，α_k为步长，g_k为梯度。

求得速度模型v，每次随机生成新的编码，且每ni次提高正演时子波主 频，ni大于每个编码炮包含的单炮集数m；

k、平滑速度模型；

l、通过直接被动源逆时偏移，利用平滑速度模型v成像，这里采用震源归 一化的互相关成像条件：

这里I(x,z)为逆时偏移结果，x,z为积分空间，T_obs(x,z,t_i)为接收到传动波 场，nm为叠加次数，m＝K/nt，K为被动源数据总采集时间，nt为单次叠加 时间。

有益效果：

本发明将直接偏移成像条件应用于被动源逆时偏移，对被动源FWI进行了 针对性的改进，保证反演正确的同时提高了计算效率，并且将直接偏移的边界 条件进入被动源偏移成像，取得了不错的效果。有以下优点：①改进的震源编 码全波形反演提高了计算的效率；②改进的参考道选择方式，适应于被动源褶 积构造的目标函数全波形反演，得到了正确的结果；③直接成像条件的被动源 逆时偏移对速度模型依赖程度低，可以在低精度模型下求得正确的偏移结果； ④被动源数据接收成本低，相对于主动源数据不需要人工震源的激发，通过本 方法可以得到高质量的偏移结果，有效提高了生产效率。

附图说明：

图1基于全波形反演驱动的被动源直接偏移成像方法流程图。

图2被动源地震干涉法示意图。

(a)一次直达波；

(b)一次反射响应；

(c)重构的虚拟炮记录；

图3试算重构记录采用的模型以及震源。

(a)模型与震源点分布；

(b)截取的一段噪声震源；

图4一组炮集记录。

(a)接收到的被动源传动波场数据；

(b)选取了一炮互相关法重构的记录；

(c)相同模型下主动源正演记录；

图5直接成像原理示意图。

(a)一次上传直达波场；

(b)一次下传反射波场；

(c)二次上传波场。

图6虚拟炮频谱分析以及震源编码后的炮记录。

(a)虚拟单炮记录的频谱分析；

(b)随机选取的一组震源编码后的编码炮；

图7震源编码的被动源全波形反演数值算例。

(a)初始线性梯度模型；

(b)震源为8Hz子波时的迭代结果；

(c)震源为20Hz子波时的迭代结果；

(d)震源为32Hz子波时的迭代结果；

图8平滑结果数值算例。

(a)平滑后的被动源全波形反演结果；

(b)采用平滑FWI速度模型直接逆时偏移结果；

(c)采用平滑真实速度模型的直接成像逆时偏移结果；

具体实施方式：

下面结合附图和实例对本发明进一步的详细说明。

基于全波形反演驱动的被动源直接偏移成像方法，包括以下步骤：

a、地表采集获得噪声被动源数据；

b、利用地震干涉法，重构采集的噪声数据，重构为虚拟炮集,一维反射波 干涉基本原理为，被动源震源在地下，检波器x1接收来自地下震源的数据， 其中包含了来自各个震源的不同延迟的信号数据(图2(a))。检波器x1接 收的波被地面反射，后被地下的地质体散射，波上行至检波器x2，被检波器 x2接收(图2(b))。由于经过检波器x1到x2的传播过程，此时x2接收 的数据中包含了新的延迟。对比图2a和2b发现，x1和x2接收的数据中同时 包含了从震源到检波器x1的传播延迟的数据。此相同的传播路径延迟可以通 过将两检波器的数据做互相关加以消除。从而生成以x1为震源，经地下地质 体散射至地表，被x2所接收的地震数据(图2(c)所示)：

C(x₂,x₁,t)＝T(x₁,x_S,t)*T(x₂,x_S,-t)

式中，T(x₁,x_S,t)为检波器x1接收到的传动波场数据，C(x₂,x₁,t)为通过互 相关干涉得到的序列；

c、互相关后的序列中包含因果和非因果部分，将因果非因果部分加和得到 虚拟震源炮集,时间域中，通过格林函数法可知

G(x₂,x₁,t)＝G(x₁,x_s,t)*G(x₂,x_s,-t)

其中，G(x₂,x₁,t)代表震源在点x₁处，在x₂处接收到的格林函数。

对于噪声源，检波器接收的数据可以表示为格林函数和噪声序列的褶积， 则上式可表示为：

C(x₂,x₁,t)＝T(x₁,x_s,t)*T(x₂,x_s,-t)

其中C(x₂,x₁,t)可表示为：

C(x₂,x₁,t)＝(G(x₂,x₁,t)+G(x₂,x₁,-t))*S(t)

S(t)为虚拟震源的子波，T(x₁,x_S,t)为检波器x1接收到的传动波场数据， C(x₂,x₁,t)为通过互相关干涉得到的序列。互相关后的序列中包含因果和非因果 部分，本方法中将因果非因果部分加和得到虚拟震源炮集：

u_vir(x₂,x₁,t)＝G(x₂,x₁,t)*S(t)

d、把重构的单炮记录做频谱分析，确定反演子波范围；

e、将得到的虚拟炮集采用全新的方式进行编码，编码方式如下所示：

为第i个编码炮，编码炮数量为nc。将各单炮数据采用震源位置随 机，得到一组随机序列，依次从随机序列中取m个值，得到nc个序列s，每个 序列s即为每个编码炮的震源位置编码，公式中p是随机生成的+1和-1序列， 长度为ns；

f、从构成该编码炮的单炮中随机选择一炮，令其炮点位置为参考道，正演 模拟记录与观测记录采用相同的参考道位置，相应的参考道为：

u_scal为编码后的正演模拟数据；

g、将编码炮作为观测数据，采用褶积的方法构造目标函数：

通过此目标函数消除了虚拟震源子波的影响；

h、同过伴随源法求取梯度，褶积构造的目标函数其伴随源为：

r'是第一互相关残差，作为除参考道之外的道集伴随震源，r”是第二互相 关残差，作为参考道处的伴随震源，梯度计算公式为：

r为将第一第二互相关残差合并后的总残差,v是迭代的速度模型，ns为震 源数，nr为检波器数，t0为记录时间；

i、梯度通过改进的指数加权算法求得：

这里α_k为线性搜索后的步长，α_k∈(0,1)，a为基于信噪比的收敛因子，k 是参与加权的迭代数。且这里的步长是通过每一编码炮残差求得的步长取算数 平均后得到的；

j、重复迭代步骤e到步骤i，迭代速度模型:

v_k+1＝v_k-α_kg_k

v_k为第k次迭代的速度模型，α_k为步长，g_k为梯度。

求得速度模型v，每次随机生成新的编码，且每ni次提高正演时子波主 频，ni大于每个编码炮包含的单炮集数m；

k、平滑速度模型；

l、通过直接被动源逆时偏移，利用求得的速度模型v成像，直接成像的原 理同互相关虚拟炮集原理基本一致，如图5所示，对于每一检波器接收到的被 动源数据，包含了(图5(a))中大量的地下震源传播到地表的一次上传直达 波波场U0，一次直达波经自由表面反射后，产生(图5(b))一次下传反射 波场D1，经地下地质体反射产生(图5(c))上传至地表U1波场，这里D1 和U1构成一对震源波场和检波点波场，波场U1经过自由表面反射，再次作 为震源向下传播后被检波器接收，构成二阶震源检波点波场。依次类推，产生 n阶震源检波波场。由于被动源数据中本身包含了作为震源的数据，直接用接 收到的被动源数据本身替代逆时偏移中的正演数据，同时被动源数据中也有相 应的检波点信息，将数据反传，之后使用成像条件，便得到了地下地质体偏移 数据。注意到，由于地下噪声源是位置随机分布，且发生时间不可预知，无法 准确提取出(图5(a))一次直达波U0，波场U0只作为第一次反射中的震 源。

这里采用震源归一化的互相关成像条件：

这里I(x,z)为逆时偏移结果，x,z为积分空间，T_obs(x,z,t_i)为接收到传动波 场，nm为叠加次数，m＝K/nt，K为被动源数据总采集时间，nt为单次叠加 时间。

实施例：

a、输入噪声记录以及地下速度模型；采用Marmousi模型来测试(图3 (a))，模型大小121×365，网格间距10m，1000个震源随机分布在地 下1000m-1200m，水平20m-3600m之间，震源位置随机，震源震动时长 随机，均为无序噪声，地表均匀分布92个检波器，间距40m，图中黑框部分 是随机震源位置，截取一段噪声源(图3(b))，为了模拟真实环境中噪声源 振幅较弱的特性，这里每个噪声源振幅数量级为10^-5；

b、利用地震干涉法，重构采集的噪声数据，重构为虚拟炮集：

C(x₂,x₁,t)＝T(x₁,x_S,t)*T(x₂,x_S,-t)

式中，T(x₁,x_S,t)为检波器x1接收到的传动波场数据，C(x₂,x₁,t)为通过互 相关干涉得到的序列；

c、互相关后的序列中包含因果和非因果部分，将因果非因果部分加和得到 虚拟震源炮集,由格林函数法可得：

u_vir(x₂,x₁,t)＝G(x₂,x₁,t)*S(t)

S(t)为虚拟震源的子波。共生成92炮虚拟炮，选取地震干涉法得到的第 46炮虚拟炮集(图4(b))，对比真实模型正演得到的46炮单炮记录，可 以看到与主动源记录相比，虚拟炮集的数据中含有大量的噪声，这是由于被动 源数据本身含有有噪声以及地震干涉法导致的二次噪声；

d、把重构的单炮记录做频谱分析，确定反演子波范围；这里抽取了虚拟炮 中的第46炮为参考(图4(b))，通过频谱分析(图6(a))可知，其频 率主要分布在8Hz～32Hz之间，采用多尺度FWI反演策略，由于其卷积构造 的目标函数特性，在分频率多尺度反演的过程中不需要对反演的观测数据进行 频率域滤波，仅需要更改相应的正演子波主频即可。

e、将得到的虚拟炮集采用全新的方式进行编码，编码方式如下所示：

为第i个编码炮，编码炮数量为nc。将各单炮数据采用震源位置随 机，得到一组随机序列，依次从随机序列中取m个值，得到nc个序列s，每个 序列s即为每个编码炮的震源位置编码，公式中p是随机生成的+1和-1序列， 长度为ns；图6(b)为随机选取的，采用本方法提出的编码和参考道策略的 一个震源编码的编码炮。这里将92个虚拟单炮分为12炮编码炮。

f、从构成该编码炮的单炮中随机选择一炮，令其炮点位置为参考道，正演 模拟记录与观测记录采用相同的参考道位置；

g、将编码炮作为观测数据，采用褶积的方法构造目标函数：

通过此目标函数消除了虚拟震源子波的影响；

h、同过伴随源法求取梯度，褶积构造的目标函数其伴随源为：

r'是第一互相关残差，作为除参考道之外的道集伴随震源，r”是第二互相 关残差，作为参考道处的伴随震源，梯度计算公式为：

r为将第一第二互相关残差合并后的总残差,v是迭代的速度模型，ns为震 源数，nr为检波器数，t0为记录时间；

i、梯度通过改进的指数加权算法求得：

式中α_k为线性搜索后的步长，α_k∈(0,1)，a为基于信噪比的收敛因子，k 是参与加权的迭代数。且这里的步长是通过每一编码炮残差求得的步长取算数 平均后得到的；

j、重复迭代步骤e到步骤i，迭代速度模型:

v_k+1＝v_k-α_kg_k

v_k为第k次迭代的速度模型，α_k为步长，g_k为梯度。

求得速度模型v，每次随机生成新的编码，且每ni次提高正演时子波主 频，ni大于每个编码炮包含的单炮集数m；这里令迭代次数为20次，以线性 梯度模型为初始模型，分别采用了8Hz、20Hz、32Hz频率的雷克子波，每个 频率迭代20次，每个频率迭代后的结果为下一个频率的初始模型，分别得到 了不同频率下的结果(图7)。这里对于加权梯度的k值，我们选择10为合适 的参数，a为0.5。可以看到从8Hz子波(图7(b))到32Hz子波(图7 (c))精度逐步提高；

k、平滑速度模型；这里采用高通滤波器滤波，得到平滑后的速度(图8 (a))；

l、通过直接被动源逆时偏移，利用平滑速度模型v成像，这里采用震源归 一化的互相关成像条件：

式中I(x,z)为逆时偏移结果，x,z为积分空间，T_obs(x,z,t_i)为接收到传动波 场，nm为叠加次数，m＝K/nt，K为被动源数据总采集时间，nt为单次叠加 时间。基于FWI最终得到的速度模型，采用直接成像的边界条件，通过逆时偏 移进行成像(图8(b))，初步反应了地下构造，平滑速度模型结果(图8 (c))成像质量更稳定。

Claims (1)

1.一种基于全波形反演驱动的被动源直接偏移成像方法，其特征在于包括以下步骤：

a、地表采集获得噪声被动源数据；

b、利用地震干涉法，重构采集的噪声数据，重构为虚拟炮集：

C(x₂,x₁,t)＝T(x₁,x_S,t)*T(x₂,x_S,-t)

式中，T(x₁,x_S,t)为检波器x1接收到的传动波场数据，C(x₂,x₁,t)为通过互相关干涉得到的序列；

c、互相关后的序列中包含因果和非因果部分，将因果非因果部分加和得到虚拟震源炮集，由格林函数法得：

u_vir(x₂,x₁,t)＝G(x₂,x₁,t)*S(t)

S(t)为虚拟震源的子波；

d、把重构的单炮记录做频谱分析，确定反演子波范围；

e、将得到的虚拟炮集采用全新的方式进行编码，编码方式如下所示：

为第i个编码炮，编码炮数量为nc，将各单炮数据采用震源位置随机，得到一组随机序列，依次从随机序列中取m个值，得到nc个序列s，每个序列s即为每个编码炮的震源位置编码，公式中p是随机生成的+1和-1序列，长度为ns；

f、从构成该编码炮的单炮中随机选择一炮，令其炮点位置为参考道，正演模拟记录与观测记录采用相同的参考道位置；

g、将编码炮作为观测数据，采用褶积的方法构造目标函数：

通过该目标函数消除了虚拟震源子波的影响；

h、同过伴随源法求取梯度，褶积构造的目标函数其伴随源为：

r'是第一互相关残差，作为除参考道之外的道集伴随震源，r”是第二互相关残差，作为参考道处的伴随震源，梯度计算公式为：

r为将第一第二互相关残差合并后的总残差,v是迭代的速度模型，ns为震源数，nr为检波器数，t0为记录时间；

i、梯度通过改进的指数加权算法求得：

式中α_k为线性搜索后的步长，α_k∈(0,1)，a为基于信噪比的收敛因子，k是参与加权的迭代数，且这里的步长是通过每一编码炮残差求得的步长取算数平均后得到的；

j、重复迭代步骤e到步骤i，迭代速度模型:

v_k+1＝v_k-α_kg_k

v_k为第k次迭代的速度模型，α_k为步长，g_k为梯度；

求得速度模型v，每次随机生成新的编码，且每ni次提高正演时子波主频，ni大于每个编码炮包含的单炮集数m；

k、平滑速度模型；

l、通过直接被动源逆时偏移，利用平滑速度模型v成像，采用震源归一化的互相关成像条件：

式中I(x,z)为逆时偏移结果，x,z为积分空间，T_obs(x,z,t_i)为接收到传动波场，nm为叠加次数，m＝K/nt，K为被动源数据总采集时间，nt为单次叠加时间。