地球物理学报  2019, Vol. 62 Issue (7): 2645-2659   PDF    
基于归一化能量谱目标函数的全波形反演方法
邢贞贞1, 韩立国1, 胡勇1, 张晓培2, 尹语晨1     
1. 吉林大学地球探测科学与技术学院, 长春 130026;
2. 吉林大学建设工程学院, 长春 130026
摘要:基于L2范数的常规全波形反演目标函数是一个强非线性泛函,在反演过程中容易陷入局部极小值.本文提出归一化能量谱目标函数来缓解全波形反演过程中的强非线性问题,同时能够有效地缓解噪声和震源子波不准等因素的影响.能量谱目标函数是通过匹配观测数据与模拟数据随频率分布的能量信息来实现最小二乘反演的,其忽略了地震数据波形与相位变化的细节特征,这在反演的过程中能够有效缓解波形匹配错位等问题.数值测试结果表明,基于归一化能量谱目标函数在构建初始速度模型、抗噪性和缓解震源子波依赖等方面都优于归一化全波形反演目标函数.金属矿模型测试结果表明,即使地震数据缺失低频分量,基于归一化能量谱目标函数的全波形反演方法在像金属矿这样的强散射介质反演问题上同样具有一定的优势.
关键词: 全波形反演      多尺度      归一化      能量谱目标函数      金属矿勘探     
Full waveform inversion based on normalized energy spectrum objective function
XING ZhenZhen1, HAN LiGuo1, HU Yong1, ZHANG XiaoPei2, YIN YuChen1     
1. College of Geo-Exploration Science and Technology, Jilin University, Changchun 130026, China;
2. College of Construction Engineering, Jilin University, Changchun 130026, China
Abstract: The L2 norm-based conventional Full Waveform Inversion (FWI) method is a strong nonlinear functional, it is easy to fall into a local minimum during the inversion process. In this paper, the Normalized Energy Spectrum based Full Waveform Inversion (NESFWI) is proposed to alleviate the strong nonlinear problem and reduce the dependence of high Signal to Noise Ratio (SNR) seismic data and source wavelet. The objective function of NESFWI is to match the energy spectrum of the observed and synthetic data in the frequency domain, which ignores the waveform and phase oscillation characteristics of the seismic data. It can effectively mitigate the waveform mismatching problem during the inversion process, while reducing dependence of seismic source wavelet and high SNR seismic data. Numerical examples show that NESFWI results are better than the Normalized Full Waveform Inversion (NFWI) in building an initial velocity model, noise immunity and source wavelet dependence. Tests on a metallic deposit model show that even if the seismic data lack low frequency components, it still can get better inversion results than the conventional FWI method.
Keywords: Full waveform inversion    Multi-scale    Normalized    Energy spectrum objective function    Metallic deposit exploration    
0 引言

在地震勘探领域,全波形反演(Full Waveform Inversion,FWI)的核心思想是利用模拟数据来匹配观测数据中所有波形信息,并通过优化算法对给定的初始速度模型按照梯度方向不断地更新,最终获得地下模型参数的高分辨率反演结果.FWI随着伴随状态法在梯度求取过程中的应用(Lailly,1983Tarantola,1984),极大程度地缩短了每次更新迭代所需要的计算时间,使得该方法逐渐成为地球物理领域的热点问题.尽管FWI有很多优点,但是由于地震信号的复杂性,使得基于L2范数的常规FWI方法变成一个强非线性优化问题(Symes,2008Virieux and Operto, 2009).再者由于检波器的接收频带有限,实际采集到的地震数据中不含低频分量,导致了局部优化算法对初始速度模型进行更新的时候,目标函数容易陷入局部极小值,无法得到准确的反演结果.这给FWI实际生产应用带来了巨大的挑战(Virieux and Operto, 2009).近年来针对FWI方法的研究主要集中于如何构建一个弱非线性目标函数,来缓解波形反演过程中目标函数陷入局部极小值的问题.

为了解决基于L2范数FWI方法的强非线性优化问题,近年来已经提出了很多新的目标函数.不同的目标函数在反演过程中对应地震数据中不同信息(包络,走时,相位等),这便使得不同的目标函数具有不同的模型反演精度、抗噪能力和稳定性.例如:包络目标函数(Chi et al., 2014Wu et al., 2014董良国等,2015Luo and Wu, 2015Bharadwaj et al., 2016罗静蕊等,2016胡勇等,2017Hu et al., 2017Wu and Chen, 2017Ha et al., 2018)能够有效地重建地震数据中缺失的低频分量,来缓解在反演过程中出现的周波跳跃现象,可为常规FWI提供一个较好的初始速度模型.Luo和Schuster(1991)利用模拟数据与观测数据的走时差信息来构建速度反演的目标函数.由于走时信息与地下速度模型有着近似线性的对应关系,因此能够有效地避免目标函数陷入局部极小值等问题.相位信息与走时信息类似(Fichtner et al., 2008Bozda et al., 2011Luo et al., 2016Hu et al., 2018),相比于地震数据的波形信息,相位信息与地下介质具有更强的线性关系.互相关目标函数(Choi and Alkhalifah, 20122013)与相位反演类似,主要是强调观测数据与模拟数据的相似程度.后来Liu等(2016)Tao等(2018)在此基础上完善地分析了互相关目标函数在全波形反演应用中的优势与劣势.近年来能量运输目标函数(Engquist and Hamfeldt, 2013Métivier et al., 2016Yang and Engquist, 2017)得到了很好地应用.通过能量运输的思想来构建目标函数,其忽略了很多波形变化的细节信息,因此能够很好地恢复速度模型中的宏观尺度信息.Wu和Chen(2017)提出利用直接包络导数代替链式法则来解决强反差速度模型的反演问题,并对整个反演流程及理论方法做了重新定义(Chen et al., 2018Hu et al., 2018Wu and Chen, 2018Zhang et al., 2018Hu et al., 2019).Shin和Ha(2017)年提出利用累积能量目标函数来进行FWI测试,并在反演过程中利用扰动点理论来得到FWI的梯度.这从计算量的角度考虑,该方法暂时还不适合在实际生产中应用.

此外,震源子波和地震数据的信噪比对反演结果有着重大影响.其中通过构建目标函数来消除震源子波影响的方法有:Rickett等(2013)提出利用变量投影方法缓解了震源子波对频率域FWI的影响;Choi和Alkhalifah(2011)敖瑞德等(2015)利用交叉褶积目标函数消除震源子波对时间域FWI的影响;Zhang等(2017)利用褶积与反褶积之后的地震数据来构建目标函数,在消除了震源子波对反演结果的影响的同时,成功地重构地震数据缺失的低频分量.关于噪声对FWI的影响,一直是一个备受关注的话题.其中基于L2范数的FWI目标函数对高斯噪声并不敏感,即当在地震数据中加入微弱的高斯噪声后对反演结果影响不大.但是当信噪比较低的情况下,常规FWI目标函数很难得到较好的反演结果.基于L1范数FWI目标函数具有较好的抗噪性,但是L1范数目标函数在求解伴随震源的过程中会出现奇异值的问题(Liu et al., 2016),导致整个FWI流程不稳定.互相关目标函数能够很好地压制噪声干扰(Liu et al., 2016),这主要是因为观测数据中的噪声成分在模拟数据中没有与其匹配的波形(Liu et al., 2016).

为了解决以上方法存在的缺陷,本文提出归一化能量谱全波形反演方法(Normalized Energy Spectrum Full Waveform Inversion,NESFWI)来缓解常规FWI目标函数强非线性问题.能量谱目标函数主要强调的是观测数据与模拟数据对应的能量随着频率的分布特征,类似于包络目标函数中能量随着时间的分布特征.由于目标函数中只考虑了地震数据的能量信息,而忽略了地震波形变化的细节信息,因此对应的梯度具有较低的模型分辨率.文章首先介绍NESFWI目标函数构建方法,以及其对应的梯度和伴随震源公式推导;并从目标函数和伴随震源的角度,理论分析其不依赖震源子波和抗噪性等优势.此外本文还结合了低通滤波多尺度反演策略,将地震数据按照频率分成低频段和高频段,先利用低频段长波长的地震数据反演得到一个较好的初始速度模型,然后再利用高频段数据恢复速度模型中的细节信息.将NESFWI的结果与归一化FWI结果(Normalized Full Waveform Inversion,NFWI)进行对比分析,证明了其在缓解周波跳跃、抗噪性和不依赖震源子波等方面的优势.利用Marmousi模型数据测试了本文方法的可行性,并进一步利用金属矿强散射速度模型分别测试NESFWI和NFWI方法,结果表明NESFWI方法在强散射介质中同样具有很好的效果.

1 目标函数与梯度

全波形反演的目标函数决定着地震数据反演流程中所利用的信息,并将直接影响到最终的反演结果.为了缓解基于L2范数的FWI目标函数强非线性问题,本文利用随频率分布的地震数据能量谱信息来构建一个新的目标函数.

1.1 归一化FWI目标函数及其梯度

归一化全波形反演(NFWI)是利用归一化的模拟数据不断地去拟合归一化的观测数据来得到地下速度模型高分辨率的反演结果,关于时间域NFWI目标函数可以定义为

(1)

其中u表示时间域模拟数据;d表示时间域观测数据;‖u2、‖d2分别表示单炮模拟数据和观测数据的L2范数,则NFWI目标函数对速度参数的梯度有:

(2)

其中表示沿着检波器方向和时间轴方向的二重积分.

根据公式(2)时间域NFWI伴随震源可以表示为

(3)

当有了伴随震源,便可以采用伴随状态法来计算目标函数对速度参数的梯度,即正传波场与伴随震源反传波场的零延迟互相关(Pratt et al., 1998).从公式(3)可以看出,伴随震源中存在一个积分项,这一项可以理解为在归一化模拟数据上加一个权重因子,用以校正模拟数据和观测数据的能量差异,则时间域NFWI的目标函数的梯度可以表示为

(4)

其中(L-1)Tfs表示时间域NFWI伴随震源的反传波场,其中T表示矩阵的转置,则多源NFWI目标函数的梯度最终可以简写为

(5)

其中Pb表示由伴随震源向地下传播得到的反传波场,Pf表示由震源发出的正演模拟波场.则公式(5)即表示正传波场与伴随震源反传波场的零延迟互相关.

1.2 归一化能量谱目标函数及其梯度

本文提出的归一化能量谱全波形反演(NESFWI)的目标函数是利用模拟数据的归一化能量谱不断地去拟合观测数据的归一化能量谱,其忽略了地震数据波形剧烈振荡的特性,因此其相对于常规时间域FWI目标函数具有较低的模型分辨率,适合用于构建初始速度模型.则NESFWI的目标函数可以定义为

(6)

其中表示频率域模拟数据和观测数据;表示模拟数据和观测数据在频率域对应的能量谱.为了便于目标函数的梯度公式推导,可将观测数据和模拟数据的L2范数改写成积分表示形式.并利用链式法则对模型速度参数求梯度,则NESFWI目标函数对应梯度可以表示为

(7)

其中Re[·]表示取数据的实部,且表示沿着检波器方向和频率轴方向的二重积分.

对公式(7)继续对速度参数求导并简化得到:

(8)

其中,将其代入到公式(8)中,可以得到目标函数对地下速度参数的梯度为

(9)

其中*表示复数的共轭.进一步简化可以得到NESFWI的目标函数对地下速度参数的梯度为

(10)

其中F-1[·]表示傅里叶反变换.

由公式(10)可以得到归一化能量谱目标函数的伴随震源表达式,即:

(11)

经傅里叶反变换后,NESFWI目标函数的伴随震源计算公式类似于常规时间域FWI的伴随震源.则多源NESFWI目标函数对应的梯度公式同样可以简化为

(12)

其中Pb表示由NESFWI的伴随震源向地下传播的反传波场,Pf表示由震源发出的正演模拟波场.从公式(11)与可以看出,NESFWI的伴随震源在频率域存在模拟数据与观测数据的乘积运算,类似于一个滤波的作用.当观测数据和模拟数据在能量谱上不相关时,这部分不相关的信息在伴随震源中将受到一定的抑制.可见该目标函数更加强调观测数据与模拟数据在能量谱方面的相似性程度.

2 归一化能量谱目标函数优势分析 2.1 目标函数特性分析

本节通过均匀速度模型数值算例(Biondi and Almomin, 2013)来分析NESFWI和NFWI目标函数的非线性程度.图 1a是均匀速度模型,速度值为v=2.5 km·s-1,模型纵横大小为2 km×7 km,网格间距为25 m,并将该均匀速度模型称为真实速度模型.在真实速度模型上利用10 Hz主频的雷克子波进行正演模拟,得到的地震数据称为观测数据.背景速度模型同样是均匀介质,其变化范围在v0=[2 km·s-1, 3 km·s-1].图 1b是模拟数据与观测数据的残差,其中模拟数据是在背景速度(v0=2 km·s-1)上正演模拟得到.当真实速度模型与背景速度模型的速度参数相差过大时,观测数据和模拟数据的波形在时间上会完全分开,这便产生了反演过程中波形不匹配的问题,即周波跳跃现象.因此当真实速度模型与背景速度模型相差过大时,利用常规FWI目标函数方法无法收敛到全局最小点(图 1b).

图 1 均匀速度模型与数据残差 (a)均匀速度模型;(b)观测数据与模拟数据的残差. Fig. 1 Homogeneous velocity model and seismic data residuals (a) Homogeneous velocity model; (b) Residuals between synthetic and recorded data.

取第195道模拟数据和观测数据来分析NESFWI和NFWI目标函数的变化特性(图 2).当背景速度与真实速度相等的时候(v0=v=2.5 km·s-1),两个目标函数同时取得全局最小值.但当背景速度模型的速度值小于2.43 km·s-1或者大于2.58 km·s-1时(图 2点画线对应处),NFWI目标函数会分别陷入左右两侧的两个极小值中.而对于NESFWI方法在背景速度值在v0=[2 km·s-1, 3 km·s-1]变化范围内,目标函数值都可以收敛到2.5 km·s-1处的全局最小点.所以为了避免常规FWI目标函数陷入局部极小值的情况,可以选择让背景速度接近真实速度模型(起速度值应选在点画线之间v0=[2.43 km·s-1, 2.58 km·s-1]),但是这一点是很难做到的.而当反演速度模型越来越接近真实速度时,NESFWI目标函数值下降的速度越来越慢,同阶段NFWI目标函数下降速度较快.所以本文将NESFWI方法、NFWI方法和分频多尺度反演策略结合使用,即在反演的起始阶段先利用低频段NESFWI构建一个较好的初始速度模型,然后再利用高频段NFWI来得到最终高精度的反演结果.

图 2 目标函数随着背景速度参数的变化曲线 Fig. 2 Objective function versus background velocity parameter
2.2 震源子波对反演结果的影响分析

NESFWI目标函数不仅在缓解非线性问题上有着独特的优势,同时还具有不依赖震源子波的优点.本节将从理论的角度来分析NESFWI目标函数在不依赖震源子波方面的优势.观测数据和模拟数据在频率域可以表示为格林函数与震源子波的乘积形式,即:

(13)

(14)

其中G0为模拟数据对应的格林函数;G为真实地震数据对应的格林函数;s0是预测震源子波在频率域的表示形式;s是真实震源子波在频率域的表示形式.

将公式(13)与(14)代入到NESFWI目标函数中,可以得到:

(15)

从公式(15)中可以看出,该目标函数只与震源子波的能量谱的分布有关,与震源子波相位谱无关.再经过归一化以后,NESFWI目标函数与震源子波能量强弱,延时多少,以及相位差异等参数均无关.

同样采用2.1节所示参数来测试不同震源子波对应的目标函数变化特征,并从理论的角度来分析NESFWI目标函数在不依赖震源子波方面的优势.从图 3中可以看出,携带不同能量的震源子波对NFWI和NESFWI目标函数影响较小.这主要是因为利用归一化的目标函数对整个单炮数据的能量进行了校正.从图 4a中可以看出,真实震源子波与预测震源子波不仅相位存在着较大的差异(相位差为-135°),同时还存在着160 ms的走时差.利用图 4a所示子波进行测试,对应的目标函数结果如图 4b所示,可以看出NESFWI目标函数仍然是在2.5 km·s-1处取得全局的最小值.但是此时NFWI目标函数不仅在真实速度模型处取不到目标函数的最小值,且目标函数的全局最小点远远偏离真实速度参数.这主要是因为预测震源子波与真实震源子波在相位和走时信息上差异过大,使得观测数据与模拟数据波形不匹配,进而无法更新初始速度模型.

图 3 不同能量的震源子波及其对应的目标函数曲线分析 (a)震源子波波形图;(b)目标函数随着背景速度参数的变化曲线. Fig. 3 Objective function analysis for seismic wavelets with different energy (a) Seismic wavelet waveforms; (b) Objective function varying with background velocity parameter.
图 4 不同走时和相位震源子波及其对应的目标函数曲线分析 (a)震源子波波形图;(b)目标函数随着背景速度参数的变化曲线. Fig. 4 Objective function analysis for seismic wavelet with different traveltimes and phases (a) Seismic wavelet waveforms; (b) Objective function varying with background velocity parameter.
2.3 噪声对反演结果影响分析

实际数据中常常含有各种形式的噪声,当噪声达到一定级别以后,会直接影响FWI的最终反演结果.对于公式(11)NESFWI的伴随震源可以简写为

(16)

其中是频率域观测数据和模拟数据泛函. 表示模拟数据与泛函f(·)在频率域乘积,即在时间域相当于一种滤波作用.由于在FWI的过程中噪声成分只含在观测数据中,所以噪声成分在频率域模拟数据中没有与其对应的能量.这样经过乘积运算以后,噪声能量将会在伴随震源中得到很好地压制,如公式(16)所示.

FWI方法的抗噪声能力可以从以下两个角度来考虑: (1)目标函数自身对噪声的敏感程度.这只能从反演结果中验证,将在3.3节详细讨论;(2)目标函数的去噪能力.这可以从伴随震源的角度来进行详细分析.若在伴随震源中的噪声得到了压制,那么在梯度中噪声干扰也会得到相应的缓解.本节利用高斯噪声和坏道噪声分别来对NFWI和NESFWI目标函数进行测试.NFWI的伴随震源可以理解为归一化的模拟数据减去归一化的观测数据,并在模拟数据上加入一个权重因子.由于噪声只含在观测数据中,使得观测数据中的噪声会原封不动地传递到NFWI伴随震源中,如图 5a所示.从公式(16)分析可知,NESFWI目标函数利用变换域的优势,并能够有效地衰减伴随震源中的噪声成分.图 5a图 5b对比可以看出,原先较强的高斯噪声干扰,在NESFWI伴随震源中得到了有效的衰减,且远偏移距处的噪声衰减效果尤为明显.

图 5 含高斯噪声的伴随震源 (a) NFWI伴随震源;(b) NESFWI伴随震源. Fig. 5 Adjoint sources with Gaussian noise (a) NFWI adjoint source; (b) NESFWI adjoint source.

图 6设计了实际数据采集过程中可能会出现的坏道噪声干扰测试.每一道中的噪声为高斯分布的随机噪声,并且在每一道的噪声中乘以[0, 1]随机分布的常数,使得每一道的噪声干扰具有不同的能量强度.然后再将噪声随机地加到15%的地震道中. 由于NFWI目标函数不具有衰减噪声的能力(图 6a),所以坏道噪声会被完整的传递到伴随震源中.对比图 6a图 6b可以发现,在NESFWI伴随震源中坏道噪声得到了有效的衰减.

图 6 含坏道噪声的伴随震源 (a) NFWI伴随震源;(b) NESFWI伴随震源. Fig. 6 Adjoint sources with bad trace noise (a) NFWI adjoint source; (b) NESFWI adjoint source.

图 5图 6的震源位置位于47道处,由于NESFWI目标函数需要变换到频率域来构建目标函数,在这过程中由于震源处能量较强,对伴随震源有着很强的干扰作用.因此在本文所有数值测试中都对震源左右10道地震数据采取了切除处理.从图 6b中可以发现,在伴随震源的直达波以上存在多处假象,这主要是因为在频率域运算的过程中出现了吉普斯效应.这些假象在地震数据正传波场中没有与之对应的波形信息,因此在利用零延迟互相关成像条件计算梯度时,新引入的噪声并不影响最终的反演结果.

3 模型测试

利用重采样的Marmousi速度模型对本文新方法进行测试,如图 7a所示.Marmousi速度模型构造较为复杂,能够很好地测试出新方法的优越性.图 7b是反演过程中利用的初始速度模型,其速度值随着深度递增,呈现出一维变化趋势,基本不含有Marmousi速度模型任何构造细节.在速度模型表层含有深度为125 m的水层,水层的速度设置为1.5 km·s-1,整个模型的速度值变化范围从1.5 km·s-1到4 km·s-1.测试模型网格纵横数目为80×280,模型网格间距为25 m.整个模型的横向距离为7 km,纵向距离为2 km.在水中均匀设置了36个震源,且震源沉放深度为75 m.检波器沉放深度为25 m,共280个检波器均匀横向分布在模型的上表面.震源选用10 Hz主频的雷克子波,采样时间间隔为0.002 s,检波器接收时间总长度为4 s.

图 7 Marmousi速度模型 (a)真实速度模型;(b)线性初始速度模型. Fig. 7 Marmousi velocity models (a) True velocity model; (b) Linear initial velocity model.

为了充分利用低频段的有效信号,本文利用低通滤波多尺度策略将地震数据按照频率尺度分开来分步反演(Bunks et al., 1995).先利用低频分量反演速度模型的宏观构造,再加入地震数据的高频分量,得到最终高精度的反演结果.本章将地震数据分为两个频段来进行测试,即0~8 Hz定义为低频段,0~13 Hz定义为高频段.

3.1 归一化能量谱全波形反演测试

为了验证本文提出的基于归一化能量谱目标函数在缓解FWI周波跳跃方面的优势,利用图 7所示的真实速度模型和初始速度模型进行测试.NFWI在低频段的反演结果如图 8a所示,在Marmousi速度模型的左侧出现了明显的错误.再加入高频段地震数据后,利用图 8a的反演结果作为初始速度模型继续进行反演.最终的反演结果如图 8b所示,可以看出高频段的反演结果只能对那些已经恢复宏观构造的区域进一步更新模型的细节信息.然而对那些在低频段未能得到宏观构造的模型区域,在加入高频的地震数据后,依然无法反演出该区域正确的速度参数.

图 8 NFWI结果 (a)低频段反演结果(0~8 Hz);(b)以(a)作为初始速度模型的高频段反演结果(0~13 Hz). Fig. 8 NFWI results (a) Low frequency band inversion result (0~8 Hz); (b) High frequency band inversion result using (a) as an initial model (0~13 Hz).

图 8a图 9a对比可以看出,NESFWI结果得到了明显的改善,基本恢复了Marmousi速度模型的宏观信息,尤其针对Marmosi速度模型左侧的近水平层位置.这主要是因为NESFWI方法降低了目标函数的非线性程度,且忽略了地震波形振荡的细节信息.因此相对于NFWI目标函数,NESFWI目标函数能够有效地缓解反演结果陷入局部极小值的问题.但是从图 2目标函数变化曲线中可以看出,NESFWI目标函数在接近真实速度模型的时候,目标函数曲线的下降速度较为缓慢.为了快速得到高精度的反演结果,可以利用图 9a反演结果作为初始速度模型,并利用NFWI目标函数继续进行高频段地震数据反演,以快速获得最终高精度的反演结果(图 9b).

图 9 反演结果 (a)低频段NESFWI结果(0~8 Hz);(b)以(a)作为初始速度模型的高频段NFWI结果(0~13 Hz). Fig. 9 Inversion result (a) Low frequency band NESFWI result (0~8 Hz); (b) High frequency band NFWI result using (a) as an initial model (0~13 Hz).
3.2 震源子波不准确测试

本节利用图 7所示的真实速度模型和初始速度模型来测试NESFWI目标函数在不依赖震源子波方面的优势.在测试的过程中将低通滤波截断频率降低至4 Hz,主要是为了避免因地震数据主频过高影响NESFWI和NFWI震源子波测试的对比效果.将震源子波分成两种情况来考虑(如图 3a图 4a):一种是真实震源子波与预测震源子波之间只存在能量差异;另一种是真实震源子波与预测震源子波之间存在相位和走时差异.从图 4a中可以看出,两个震源子波存在160 ms的时间延迟,同时相位为-135°,使得这两个震源子波在波形上基本没有任何重叠的部分.

图 10的反演结果中可以看出,震源子波能量差异对归一化类目标函数的影响较小,两者都能很好地获得Marmousi速度模型的宏观信息.图 10b是NESFWI结果,相对来说分辨率较低,这也侧面反映了NESFWI方法在构建初始速度模型等方面的优势.

图 10 利用图 3a震源子波反演结果 (a) NFWI结果;(b) NESFWI结果. Fig. 10 Inversion results with Fig. 3a wavelet (a) NFWI result; (b) NESFWI result.

图 11是利用图 7所示速度模型和图 4a所示的震源子波得到的反演结果.当真实震源子波与预测震源子波在相位(相差-135°)与走时(相差160 ms)方面差异过大时,使得观测数据与模拟数据之间没有任何波形上的交集,导致NFWI目标函数从反演的起始阶段便无法下降,同时初始速度模型也无法得到更新(如图 11a所示).从公式(15)和图 11b的反演结果可以看出,NESFWI目标函数对震源子波的依赖性较小,即使在相位与走时信息完全错误的情况下,NESFWI也能够很好地反演出速度模型对应的宏观信息.图 11b图 10b对比发现,Marmousi速度模型深部构造反演结果基本相同,有力地证明了NESFWI目标函数在不依赖震源子波方面的优势.

图 11 利用图 4a震源子波反演结果 (a) NFWI结果;(b) NESFWI结果. Fig. 11 Inversion results with Fig. 4a wavelet (a) NFWI result; (b) NESFWI result.
3.3 噪声测试

野外采集的地震数据中常含有各种各样的噪声干扰,这些噪声的存在严重影响了FWI最终结果. B本文通过讨论高斯噪声以及坏道噪声对反演结果的影响,来分析对应目标函数的抗噪能力.为了避免地震数据主频过高而影响到噪声测试的对比效果,本节同样将低通滤波的截断频率降低至4 Hz.此外为了体现NESFWI目标函数对高斯噪声不敏感的同时,还具有压制强噪声的优势,便将地震数据的信噪比降低至SNR=-5.32,即噪声强度远远大于有效信号的强度.

图 12a中可以看出,由于高斯噪声能量强于有效信号的能量,使得NFWI结果中含有较强随机分布的干扰点.参照NFWI伴随震源(图 5a),在低信噪比的情况下,只能有部分直达波可见,其余有效波的信息完全淹没在了强高斯噪声中.同时图 10a图 12a对比可以发现,强高斯噪声严重影响NFWI目标函数最终的反演结果,无法获得正确的地质构造.相比于图 12a图 12b利用NESFWI目标函数在计算伴随震源的过程中还具有较强衰减噪声的能力(图 5b所示),这使得NESFWI结果(图 12b)受噪声影响较小,且抗噪性优于NFWI结果.

图 12 高斯噪声测试 (a) NFWI结果;(b) NESFWI结果. Fig. 12 Gaussian noise tests (a) NFWI result; (b) NESFWI result.

图 13是利用坏道噪声(图 6)进行测试的结果.在地震实际数据采集过程中,因地表扰动干扰,可能使得整道地震数据出现随机振荡的现象,像这种坏道噪声严重影响反演的精度.从图 13a可以看出,L2范数对此类噪声较为敏感(Liu et al., 2016).当地震数据中仅含有15%坏道噪声时,反演结果与图 10a相比差了很多,完全无法恢复Marmousi速度模型左侧及深部区域地质构造.从图 13b可以看出,NESFWI结果基本恢复了Marmousi速度模型的宏观构造,且有效缓解了左侧近水平构造反演陷入局部极小值等问题.图 13a图 13b对比发现,NESFWI目标函数相比于NFWI目标函数在压制坏道噪声方面具有更强的优势.

图 13 坏道噪声测试 (a) NFWI结果;(b) NESFWI结果. Fig. 13 Bad trace noise tests (a) NFWI result; (b) NESFWI result.
4 能量谱目标函数在金属矿速度模型中的应用

为了验证本文提出的NESFWI方法在金属矿地震勘探中应用的可行性,根据安徽省庐枞盆地某矿床的真实产出形态来重建复杂地球物理纵波速度模型(吕庆田等,2010廉玉广等,2011),并将NESFWI结果与NFWI结果进行对比分析.图 14a是对原始纵波速度模型进行重采样后得到的真实速度模型,模型网格点数为100×280,网格间距为25 m,速度模型横向长7 km,纵向深2.5 km,其中主矿体埋深距地表 1.2~1.9 km.该矿体总体呈复杂的似层状、平缓透镜状,空间上表现为穹隆状,中心以侵染状磁铁矿为主,富、厚矿多环于四周(孙宏宇等,2015).矿床由多个矿体组成,其中规模较大的有3个矿体,其余均为小矿体(吕庆田等,2010廉玉广等,2011吕庆田等,2014).图 14b是测试过程中利用的初始速度模型,且速度参数随着深度逐渐递增,速度变化范围从3 km·s-1到4.5 km·s-1,呈现出一维变化趋势.然而在真实速度模型中金属矿所在区域对应的速度值是5.6 km·s-1,目标矿区速度参数与初始速度参数存在着较大的差异,这也给整个反演流程带来了巨大挑战.在金属矿速度模型表面均匀分布着24个震源,每个震源对应着280个检波器.震源选用15 Hz主频的雷克子波,采样时间间隔为0.002 s,检波器接收时间总长度为3 s.

图 14 金属矿速度模型 (a)真实速度模型;(b)线性初始速度模型. Fig. 14 Metallic deposit velocity models (a) True velocity model; (b) Linear initial velocity model.

为了更好地模拟实际情况,同时证明本文提出的NESFWI方法在缓解周波跳跃方面的优势,首先对地震数据利用高通滤波器滤除5 Hz以下的低频分量,利用缺失低频分量的地震数据分别进行NFWI和NESFWI测试.此外在反演过程中还结合了低通滤波多尺度策略,沿用第3节频率分组策略,并将5~12 Hz定义为低频段数据,将5~50 Hz定义为高频段数据.NFWI方法在低频段的结果如图 15a所示,在金属矿速度模型的左侧和中间区域出现了大块假低速异常值,并且无法确定目标矿体的空间产出状态.将低频段的反演结果(图 15a)作为初始速度模型继续进行高频段NFWI,最终结果如图 15b所示.随着高频分量的加入,在金属矿浅部区域的分辨率有了明显的改善,但是假的低速异常区域依然无法恢复其对应位置真实的速度参数.同时在图 15b中,只能隐约的看到金属矿的上边界,这样的反演结果无法直接用于地质解释,甚至可能会对该地区的整个矿产普查产生错误的理解.在初始速度模型较差的情况下,NFWI方法无法得到金属矿模型精确的反演结果,主要是因为:金属矿区域速度参数明显高于围岩速度值,同时矿体呈透镜状,近似层状,具有较大的厚度,且分布在模型的中深层,这类矿体属于典型的强散射介质.Wu和Chen(2018)指出FWI方法是基于弱散射假设条件推导的梯度表达式,当在反演的过程中遇到强散射介质(强反差速度模型),且当初始速度模型与真实速度模型相差较大时,弱散射假设条件不再适用,导致常规FWI方法无法恢复强散射介质对应的物性参数.

图 15 NFWI结果 (a)低频段反演结果(5~12 Hz);(b)以(a)作为初始速度模型的高频段NFWI结果(5~50 Hz). Fig. 15 NFWI results (a) Low frequency band inversion result (5~12 Hz); (b) High frequency band NFWI result using (a) as an initial model (5~50 Hz).

图 16a图 15a对比可以看出,NESFWI的结果得到了明显的改善,并且很好地恢复了金属矿上部区域对应的速度参数.利用低频段NESFWI的结果作为高频段NFWI的初始速度模型,最终反演结果如图 16b所示.图 15b图 16b对比可以看出,NESFWI方法不仅改善了金属矿区域的速度反演结果,同时还能够有效避免在浅部弱散射介质中出现虚假低速异常等问题.但是图 14a所示的金属矿及其产出形态,属于典型的强散射介质,这在地震勘探中会对其下方的矿体存在着明显的“高速屏蔽”效应,即在初始速度模型较差的情况下,NESFWI方法也无法完整地将整个金属矿所在区域的速度参数精确地反演出来.

图 16 反演结果 (a)低频段NESFWI结果(5~12 Hz);(b)以(a)作为初始速度模型的高频段NFWI目标函数反演结果(5~50 Hz). Fig. 16 Inversion result (a) Low frequency band NESFWI result (5~12 Hz); (b) High frequency band NFWI result using (a) as an initial model (5~50 Hz).

图 17是单道速度参数对比结果,在反演模型的2.95 km处抽取对应的速度参数(图 15b16b红线所示).从图 17a中可以看出,NFWI结果只能真实的反映出浅部速度值变化较为平缓的区域.随着深度逐渐增加,尤其是金属矿所在区域,NFWI方法无法得到正确的反演结果.相比于NFWI结果,NESFWI结果(图 17b)得到了明显的改善,其不仅很好地恢复了浅层弱散射介质对应的速度参数,同时还能很好地反演出金属矿体区域顶层的速度参数.

图 17 速度参数剖面对比(速度模型2.95 km处) (a) NFWI结果;(b) NESFWI+NFWI结果. Fig. 17 Comparison of velocity parameter (at 2.95 km of velocity model) (a) NFWI result; (b) NESFWI+NFWI result.

图 18是单道速度参数对比结果,在反演模型的4.4 km处抽取对应的速度参数(图 15b16b红线所示).NFWI结果在深度为1 km处出现了一块假低速异常体,并且这个假低速异常区域严重影响了周围及其深部的反演结果,使其与真实地质层位出现空间位置偏差.然而对于NESFWI的结果,不仅在浅部没有出现低速异常的假象,同时还能很好地反演出深部层位的空间位置及其对应的物性参数.

图 18 速度参数剖面对比(速度模型4.4 km处) (a) NFWI结果;(b) NESFWI结果. Fig. 18 Comparison of velocity parameter (at 4.4 km of velocity model) (a) NFWI result; (b) NESFWI result.

本章对比了NFWI和NESFWI目标函数在金属矿产方向的应用.无论从反演速度模型的空间产出状态还是物性参数等角度来说,NESFWI结果都优于NFWI结果.这主要是因为NESFWI的目标函数降低了非线性程度,减弱了反演过程对初始速度模型的依赖,因此能够有效地缓解目标函数陷入局部极小值等问题,最终得到高精度的反演结果.

5 结论与讨论

(1) 本文提出NESFWI目标函数,其忽略了地震波形变化的细节信息,有效降低了地震反演过程中的非线性程度.因此该方法可为NFWI方法提供一个较好的初始速度模型,有效缓解初始速度模型不准确导致目标函数陷入局部极小值的问题.

(2) NESFWI目标函数是利用观测数据与模拟数据随着频率分布的能量信息来构建的,其对震源子波的相位和走时信息差异不敏感.当真实震源子波与预测震源子波在激发时间上有一个很大的延迟时,NESFWI的结果基本不受影响.因此在震源子波不准确以及初始速度模型较差的情况下,该方法依然能够有效缓解周波跳跃现象.

(3) 由于NESFWI的伴随震源在频率域涉及到观测数据与模拟数据乘积过程,这类似一个滤波作用.当模拟数据能量谱中没有与噪声对应的成分时,观测数据中的噪声将在伴随震源中得到很好的压制.强高斯噪声和坏道噪声测试结果表明,NESFWI目标函数能够很好地压制伴随震源中的多种噪声.即使在信噪比较低的情况下,利用NESFWI目标函数也能够得到正确的反演结果.

(4) 金属矿模型测试结果表明,NESFWI方法相比于NFWI方法更适用于强散射介质.但当金属矿体积较大时,且由于地震波存在“高速屏蔽”现象,NESFWI方法只能得到金属矿顶部的物性参数.

本文提出NESFWI目标函数,借助变换域的优点,即在频率域构建NESFWI目标函数,在时间域更新速度模型.因此可以将频率域的优势策略应用到时间域反演中,并为FWI方法实际应用发展一套综合的处理流程.将时间域和频率域的新理论与新策略相结合,有望能够解决像金属矿这类强反差速度模型的高精度反演问题.

References
Ao R D, Dong L G, Chi B X. 2015. Source-independent envelope-based FWI to build an initial model. Chinese Journal of Geophysics (in Chinese), 58(6): 1998-2010. DOI:10.6038/cjg20150615
Bharadwaj P, Mulder W, Drijkoningen G. 2016. Full waveform inversion with an auxiliary bump functional. Geophysical Journal International, 206(2): 1076-1092. DOI:10.1093/gji/ggw129
Biondi B, Almomin A. 2013. Tomographic full-waveform inversion (TFWI) by combining FWI and wave-equation migration velocity analysis. The Leading Edge, 32(9): 1074-1080. DOI:10.1190/tle32091074.1
Bozda E, Trampert J, Tromp J. 2011. Misfit functions for full waveform inversion based on instantaneous phase and envelope measurements. Geophysical Journal International, 185(2): 845-870. DOI:10.1111/gji.2011.185.issue-2
Bunks C, Saleck F M, Zaleski S, et al. 1995. Multiscale seismic waveform inversion. Geophysics, 60(5): 1457-1473. DOI:10.1190/1.1443880
Chen G X, Wu R S, Chen S C. 2018. Reflection multi-scale envelope inversion. Geophysical Prospecting, 66(7): 1258-1271. DOI:10.1111/gpr.2018.66.issue-7
Chi B X, Dong L G, Liu Y Z. 2014. Full waveform inversion method using envelope objective function without low frequency data. Journal of Applied Geophysics, 109: 36-46. DOI:10.1016/j.jappgeo.2014.07.010
Choi Y, Alkhalifah T. 2011. Source-independent time-domain waveform inversion using convolved wavefields:Application to the encoded multisource waveform inversion. Geophysics, 76(5): R125-R134. DOI:10.1190/geo2010-0210.1
Choi Y, Alkhalifah T. 2012. Application of multi-source waveform inversion to marine streamer data using the global correlation norm. Geophysical Prospecting, 60(4): 748-758. DOI:10.1111/gpr.2012.60.issue-4
Choi Y, Alkhalifah T. 2013. Frequency-domain waveform inversion using the phase derivative. Geophysical Journal International, 195(3): 1904-1916. DOI:10.1093/gji/ggt351
Dong L G, Huang C, Chi B X, et al. 2015. Strategy and application of waveform inversion based on seismic data subset. Chinese Journal of Geophysics (in Chinese), 58(10): 3735-3745. DOI:10.6038/cjg20151024
Engquist B, Hamfeldt B F. 2013. Application of the Wasserstein metric to seismic signals. Communications in Mathematical Sciences, 12(5): 979-988.
Fichtner A, Kennett B L N, Igel H, et al. 2008. Theoretical background for continental-and global-scale full-waveform inversion in the time-frequency domain. Geophysical Journal International, 175(2): 665-685. DOI:10.1111/gji.2008.175.issue-2
Ha J, Kim S, Koo N, et al. 2018. Full waveform inversion using a decomposed single frequency component from a spectrogram. Journal of Applied Geophysics, 153: 154-167. DOI:10.1016/j.jappgeo.2018.04.010
Hu Y, Han L G, Xu Z, et al. 2017. Adaptive multi-step full waveform inversion based on waveform mode decomposition. Journal of Applied Geophysics, 139: 195-210. DOI:10.1016/j.jappgeo.2017.02.017
Hu Y, Han L G, Xu Z, et al. 2017. Demodulation envelope multi-scale full waveform inversion based on precise seismic source function. Chinese Journal of Geophysics (in Chinese), 60(3): 1088-1105. DOI:10.6038/cjg20170321
Hu Y, Han L G, Yu J L, et al. 2018. Time-frequency domain multi-scale full waveform inversion based on adaptive non-stationary phase correction. Chinese Journal of Geophysics (in Chinese), 61(7): 2969-2988. DOI:10.6038/cjg2018L0421
Hu Y, Wu R S, Han L G, et al. 2018. Multi-Scale energy weighted envelope-phase inversion using direct envelope Fréchet derivative in local time-frequency.//88th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts.
Hu Y, Wu R S, Han L G, et al. 2019. Joint multiscale direct envelope inversion of phase and amplitude in the time-frequency domain. IEEE Transactions on Geoscience and Remote Sensing, in Press.
Lailly P. 1983. The seismic inverse problem as a sequence of before stack migrations.//Conference on Inverse Scattering: Theory and Application. Philadelphia, PA: Society for Industrial and Applied Mathematics, 206-220.
Lian Y G, Lü Q T, Han L G, et al. 2011. Research of seismic modeling of complex metal ore body-examples of the Luzong Luohe, Nihe and Dabaozhuang deposits. Acta Geologica Sinica (in Chinese), 85(5): 887-899.
Liu Y S, Teng J W, Xu T, et al. 2016. Robust time-domain full waveform inversion with normalized zero-lag cross-correlation objective function. Geophysical Journal International, 209(1): 106-122.
Luo J R, Wu R S. 2015. Seismic envelope inversion:reduction of local minima and noise resistance. Geophysical Prospecting, 63(3): 597-614. DOI:10.1111/1365-2478.12208
Luo J R, Wu R S, Gao F C. 2016. Time-domain full-waveform inversion using instantaneous phase with damping.//86th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 1472-1476.
Luo J R, Wu R S, Gao J H. 2016. Local minima reduction of seismic envelope inversion. Chinese Journal of Geophysics (in Chinese), 59(7): 2510-2518. DOI:10.6038/cjg20160716
Luo Y, Schuster G T. 1991. Wave-equation traveltime inversion. Geophysics, 56(5): 645-653. DOI:10.1190/1.1443081
Lü Q T, Dong S W, Shi D N, et al. 2014. Lithosphere architecture and geodynamic model of Middle and Lower Reaches of Yangtze Metallogenic Belt:A review from SinoProbe. Acta Petrologica Sinica (in Chinese), 30(4): 889-906.
Lü Q T, Han L G, Yan J Y, et al. 2010. Seismic imaging of volcanic hydrothermal iron-sulfur deposits and its hosting structure in Luzong ore district. Acta Petrologica Sinica (in Chinese), 26(9): 2598-2612.
Métivier L, Brossier R, Mérigot Q, et al. 2016. Measuring the misfit between seismograms using an optimal transport distance:application to full waveform inversion. Geophysical Journal International, 205(1): 345-377. DOI:10.1093/gji/ggw014
Pratt R G, Shin C, Hick G J. 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion. Geophysical Journal International, 133(2): 341-362. DOI:10.1046/j.1365-246X.1998.00498.x
Rickett J. 2013. The variable projection method for waveform inversion with an unknown source function. Geophysical Prospecting, 61(4): 874-881. DOI:10.1111/gpr.2013.61.issue-4
Shin C, Ha W. 2017. Accumulated energy norm for full waveform inversion of marine data. Journal of Applied Geophysics, 147: 91-101. DOI:10.1016/j.jappgeo.2017.10.005
Sun H Y, Han L G, Han M, et al. 2015. Elastic full waveform inversion based on visibility analysis and energy compensation for metallic deposit exploration. Chinese Journal of Geophysics (in Chinese), 58(12): 4605-4616. DOI:10.6038/cjg20151222
Symes W W. 2008. Migration velocity analysis and waveform inversion. Geophysical Prospecting, 56(6): 765-790. DOI:10.1111/gpr.2008.56.issue-6
Tao K, Grand S P, Niu F L. 2018. Full-waveform inversion of triplicated data using a normalized-correlation-coefficient-based misfit function. Geophysical Journal International, 210(3): 1517-1524.
Tarantola A. 1984. Inversion of seismic reflection data in the acoustic approximation. Geophysics, 49(8): 1259-1266. DOI:10.1190/1.1441754
Virieux J, Operto S. 2009. An overview of full-waveform inversion in exploration geophysics. Geophysics, 74(6): WCC1-WCC26. DOI:10.1190/1.3238367
Wu R S, Chen G X. 2017. New frechet derivative for envelope data and multi-scale envelope inversion.//87th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts.
Wu R S, Chen G X. 2018. Multi-scale seismic envelope inversion using a direct envelope Frechet derivative for strong-nonlinear full waveform inversion. arXiv preprint arXiv: 1808.05275.
Wu R S, Luo J R, Wu B Y. 2014. Seismic envelope inversion and modulation signal model. Geophysics, 79(3): WA13-WA24. DOI:10.1190/geo2013-0294.1
Yang Y N, Engquist B. 2017. Analysis of optimal transport and related misfit functions in full-waveform inversion.//87th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 1291-1296.
Zhang P, Han L G, Xu Z, et al. 2017. Sparse blind deconvolution based low-frequency seismic data reconstruction for multiscale full waveform inversion. Journal of Applied Geophysics, 139: 91-108. DOI:10.1016/j.jappgeo.2017.02.021
Zhang P, Wu R S, Han L G. 2018. Source-independent seismic envelope inversion based on the direct envelope Fréchet derivative. Geophysics, 83(6): R581-R595. DOI:10.1190/geo2017-0360.1
敖瑞德, 董良国, 迟本鑫. 2015. 不依赖子波、基于包络的FWI初始模型建立方法研究. 地球物理学报, 58(6): 1998-2010. DOI:10.6038/cjg20150615
董良国, 黄超, 迟本鑫, 等. 2015. 基于地震数据子集的波形反演思路、方法与应用. 地球物理学报, 58(10): 3735-3745. DOI:10.6038/cjg20151024
胡勇, 韩立国, 许卓, 等. 2017. 基于精确震源函数的解调包络多尺度全波形反演. 地球物理学报, 60(3): 1088-1105. DOI:10.6038/cjg20170321
胡勇, 韩立国, 于江龙, 等. 2018. 基于自适应非稳态相位校正的时频域多尺度全波形反演. 地球物理学报, 61(7): 2969-2988. DOI:10.6038/cjg2018L0421
廉玉广, 吕庆田, 韩立国, 等. 2011. 复杂金属矿体地震波正演模拟研究-以庐枞盆地罗河、泥河和大包庄矿床为例. 地质学报, 85(5): 887-899.
罗静蕊, 吴如山, 高静怀. 2016. 地震包络反演对局部极小值的抑制特性. 地球物理学报, 59(7): 2510-2518. DOI:10.6038/cjg20160716
吕庆田, 董树文, 史大年, 等. 2014. 长江中下游成矿带岩石圈结构与成矿动力学模型-深部探测(SinoProbe)综述. 岩石学报, 30(4): 889-906.
吕庆田, 韩立国, 严加永, 等. 2010. 庐枞矿集区火山气液型铁、硫矿床及控矿构造的反射地震成像. 岩石学报, 26(9): 2598-2612.
孙宏宇, 韩立国, 韩淼, 等. 2015. 基于可视性分析与能量补偿的金属矿弹性波全波形反演. 地球物理学报, 58(12): 4605-4616. DOI:10.6038/cjg20151222