全波形反演的一个新目标函数:数据域中的微分相似优化

高福春,Paul WILLIAMSON,R. Gerhard PRATT

(1.道达尔勘探和生产研究技术中心,休士顿77002,美国;2.西安大略大学,伦敦N6A 3K7,加拿大)

全波形反演的一个新目标函数:数据域中的微分相似优化

高福春1,Paul WILLIAMSON1,R. Gerhard PRATT2

(1.道达尔勘探和生产研究技术中心,休士顿77002,美国;2.西安大略大学,伦敦N6A 3K7,加拿大)

全波形反演(full waveform inversion,FWI)目前已有广泛的工业实践,但因其本质上的非线性,不如走时层析成像等传统速度建模技术稳健,非线性程度也因目标函数不同而不同。研究分析了FWI中几种不同目标函数的性质,基于定义在数据域中的微分相似概念,提出了一种新的目标函数。初步试验表明,这种目标函数对于比较大范围的数据残差都有凸状性质,基于梯度优化法时使用该目标函数的FWI比传统FWI更稳健,而且波形反演的良好分辨率基本得以保留。

全波形反演;目标函数;非线性;微分相似优化

三维全波形反演(FWI)在20年前还被认为不可实现,但是目前已经在工业上得到应用,这主要得益于计算机硬件和数值算法的持续发展[1-2]。标准的FWI可以得到比反射层析成像等常规速度建模方法分辨率更高的速度模型,但因其高度非线性,不如常规速度建模方法稳健。为了确保FWI实现物理意义上的收敛,基于二范数目标函数的常规FWI需要初始模型足够接近真实模型,或者数据中含有可靠的低频信息[3]。通常,人们通过折射波走时层析成像(如PRATT[4],GAO等[5])或深度域成像建模流程(包括反射波层析成像,如SIRGUE等[6])来得到接近真实模型的初始模型,然而,当介质比较复杂时,这些方法不一定能产生一个足够接近真实模型的初始模型,从而确保FWI在使用典型地震数据带宽时总是获得成功。另一方面,不同类型的目标函数对相同模型的扰动以及数据中噪声的敏感度不同。为替代传统的二范数目标函数,人们提出多种新的目标函数,主要有以下几种:①人工合成数据和实际观测数据之间基于其它范数的差异性(如TARANTOLA[7],PRATT等[8])的衡量。HUGUES等[9]使用一范数目标函数反演了一个墨西哥湾的二维数据。GUITTON等[10]推导出基于休伯范数的反演算法,并证明该算法对数据中的异常值不太敏感。KIM等[11]设计出一种FWI算法,旨在频率域里减少人工合成数据和观测数据之间基于二范数的相位差,而不是振幅和相位的整体差。ARAVKIN等[12]推导出基于学生-T范数(休伯范数的一个变种)的FWI算法,该算法对集中在某一区域的异常值不很敏感。②人工合成数据和实际观测数据相似性的度量。ROUTH等[13]设计了基于相关性目标函数的FWI算法,通过最大限度地提高计算数据和观测数据之间的相似性来更新速度模型,而不是最小化它们之间的差别。这两类替代目标函数的提出大都是为了达到某一个特定的目标。③人工合成数据和观测数据之间经过相同变换或处理以后的差异衡量。HA等[14]推导出拉普拉斯域中的FWI算法,指出该方法通过最小化拉普拉斯域中的差异,使其比传统FWI对数据中低频成分的依赖性相对减弱。LUO等[15]推导并实现了地震波波形包络(而不是波形本身)反演方法,可重建速度模型的大尺度成分。BAEK等[16]推导出地震波事件记录引导下的最小二乘反演方法,以减轻传统FWI在全域最小附近收敛半径较窄的问题。这类方法的共性是先进行数据变换或处理使得信号在其它域中具有相对低频成分,然后再比较人工合成数据和观测数据之间的差异使其发生周跳的可能性降低,减弱对初始模型的敏感性。上述方法都有其各自的缺点。定义在像域中的微分相似优化(differential semblance optimization,简称DSO)概念最初由SYMES等[17]提出,作为一种在像域中为更新速度模型而建立的目标函数,它具有凸状特性优势,受到格外重视。PRATT等[18]分析了一系列基于相似性的目标函数,初步证明DSO概念可以扩展到数据域中,并指出这些函数可作为替代目标函数,使FWI更加稳健。虽然PRATT等指出了这些目标函数对于模型的扰动比常规目标函数更线性,然而他们没有实现用这些函数作为目标函数的FWI算法,所以这些函数在理论上拥有的优点还有待进一步证实。本文提出了与PRATT等相似但更加线性化的目标函数,分析并选择了其中一种作为FWI的目标函数并予以实现。该FWI算法更加稳健,分辨率较传统的全波形反演有所降低,但较常规射线类反演方法有明显提高。将这种目标函数与传统的目标函数联合使用能得到分辨率较高的模型。

1 方法原理

1.1 不同的目标函数

传统FWI定义在频率和空间域中基于二范数的目标函数(如PRATT等[8])为:

(1)

式中:u和d分别是人工正演模拟和观察到的波场;“*”表示复共轭。PRATT等[18]提出了如下相干性的衡量:

其中,s被定义为震源因子:

(4)

式中:“T”表示转置。公式(4)中的震源因子可用m个地震道集来计算,1≤m≤N,N为共震源道集中所有的地震道数量。公式(3)中的导数可以是对震源空间位置或检波器空间位置的导数。

在本次研究工作中,我们提出了另外一种新的目标函数,称之为窗口式相似性目标函数(简称WS目标函数)。该函数在原理上和公式(3)中描述的目标函数没有本质的区别,但是更便于实现,在实际应用中更具有灵活性。

为了比较由方程(1)—(3)所定义的目标函数和WS目标函数对于模型扰动的不同敏感程度,利用PRATT等[18]使用的二维人工合成模型,采用井间数据采集方式对4种目标函数进行了测试。两口井水平距离150m,井深300m,分别布置51个震源和51个检波器。模型的速度只是深度的函数,即一维结构模型。其速度分布是一个常背景速度加上一恒定的梯度(v=1000+5z),所以速度随深度增加,从1000m/s上升到底部的2000m/s,见图1a。模型中使用的梯度扰动了20次,每次改变5%,从0一直变化到正确梯度的95%,便得到了20个不同的扰动模型,参见图1b。

图2显示了4种不同目标函数随速度扰动的演化过程，可以看到：除了常规的二范数目标函数,其它3种目标函数的变化中都不存在区域最小。常规目标函数的变化中存在区域最小,而且不止一个,如图2a 所示。采用梯度法时,只有当初始模型的梯度值是正确值的85%时,常规FWI才会收敛到全域最小。从图2b至图2d中3条曲线可以看到,如果它们被当作目标函数,即使速度模型离正确值很远,也可以收敛到全域最小。但这3个函数的演化过程不同:图2b中由公式(2)定义的函数变化一直缓慢,直到离正确模型很近时才加快收敛,其它两个函数的收敛则比较均匀。由于WS函数相对比较容易实施,我们这里选择WS函数作为目标函数。

图1 速度模型(a)和各种扰动后的速度随深度的线性变化(b)

图2 四种目标函数收敛性随速度扰动的变化a 公式(1); b 公式(2); c 公式(3); d WS目标函数

1.2 使用WS目标函数的FWI

与主流FWI的实施过程一样,我们也是用迭代法使得目标函数变小来更新模型,即:

(5)

2 应用实例

为了验证WS目标函数的有效性,我们用Marmousi模型数据进行了一系列测试。该数据为典型的浅海采集数据,具有复杂的沉积构造和很强的横向不均匀性(图3a)。在2.0km深、9.2km宽的区域内,速度变化范围为1500～5500m/s,共有93个共炮点道集,检波器在地表均匀分布,位置固定不变,间距为50m,共185个,炮点间距100m。

首先使用一个较好的初始模型进行了一系列测试。该初始模型为真实模型的平滑结果,见图3b。依次反演了3个频率分量,分别为3.0,3.5和4.0Hz。

图3 Marmousi模型数据测试结果a 真实模型; b 初始模型; c 常规FWI反演结果; d 基于WS目标函数(窗口宽度3)的FWI反演结果; e 基于WS目标函数(窗口宽度21)的FWI反演结果

采用WS目标函数时,我们选择了两个不同的窗口宽度,分别为3和21。图3c为常规FWI反演的最终速度模型,图3d和图3e均是用WS目标函数FWI反演的最终结果,前者窗口宽度为3,后者窗口宽度为21。可以看到,用WS目标函数时,不管窗口长度如何,在反演频率允许的分辨率范围内,都可以成功地重建真实模型,但它们的分辨率,尤其是采用的窗口宽度比较大时得到的反演模型,要稍低于传统FWI得到的速度模型。从沉积层和模型边缘尤其能看出分辨率的不同。沉积层分辨率的差异可归因于目标函数的不同,边缘分辨率的差异可归因于窗口宽度的不同。由上述分析可以看出,WS目标函数对模型的扰动敏感性降低了,所以FWI反演的结果分辨率会低一些。窗口宽度越大,边缘效应就越强,所以模型的边缘部分受到的影响较大。

以上测试只是确定了WS函数可以作为FWI的目标函数,但是还没有证实它的收敛半径较大。为了进一步测试该目标函数的性质,我们对不同目标函数做了一系列测试(图4,图5)。采用的两个初始模型(模型1和模型2)都是一维结构模型,即速度只是深度的函数,离真实模型比较远。这两个模型虽然都是一维模型,但是差异比较大。模型1速度从地表的1500m/s以常梯度变化到底部的2625m/s,而模型2 则从地表1500m/s变化到底部的2000m/s。显而易见,模型2偏离真实模型更远,几乎是常速度模型。用模型1(图4b)反演了2～6Hz的数据,常规FWI得到的最终模型(图4c)显然陷入了区域最小,而用WS目标函数得到的最终模型(图4d)则比较合理地重建了真实模型。将图4d对应的速度模型作为初始模型,我们可以通过常规FWI得到分辨率较高的速度模型(图4e)。

图4 初始模型用模型1的测试结果a 真实模型; b 初始模型; c 常规FWI反演的最终模型; d 用WS目标函数FWI得到的最终模型; e 常规FWI用WS反演模型作为初始模型得到的最终模型(2～6Hz)

图5 初始模型用模型2的测试结果a 初始模型; b 常规FWI反演的最终模型; c 用WS目标函数FWI得到的最终模型; d 常规FWI用WS反演模型作为初始模型得到的最终模型(2～6Hz)

图5为用模型2测试的结果。因为模型2离真实模型更远,可以看到常规FWI结果误差很大(图5b),用WS目标函数的FWI则可以在一定区域内重建真实模型(图5c)。这里是指左右下角之外的区域,原因在于这两个区域受到地震波照明的限制。用WS目标函数的FWI结果作为初始模型,我们又能通过常规FWI得到分辨率比较高的模型(图5d)。

3 讨论和结论

本文展示并实现了一个FWI常规目标函数的替代函数,即WS目标函数。它的提出受到了PRATT等于2002年提出的几个目标函数的启发。WS目标函数和PRATT等提出的目标函数有相似的性质,但相对比较容易实现,而且可以通过窗口的宽度调节目标函数的性质。Marmousi模型测试证实WS目标函数具有比较大的收敛半径,较少的局部最小值,能使FWI变得稳健。下一步研究工作是对该算法作进一步测试,并将其用于实际地震数据。

WS目标函数会使FWI的分辨率稍微变低,这可以理解而且不是问题,因为该目标函数可以和常规目标函数结合起来使用,使得FWI既稳健又具有很好的分辨率。

致谢:感谢道达尔勘探研究和生产技术中心允许发表这项研究成果,感谢Laurent Sirgue为本研究中的测试提供了所需要的初始模型。

[1] SIRGUE L,ETGEN J,ALBERTIN U.3D full waveform inversion:wide-versus narrow-azimuth acquisitions[J].Expanded Abstracts of 77thAnnual Internat SEG Mtg,2007:1760-1763

[2] PLESSIX R E,GUIDO B,DE MAAG J W,et al.Application of acoustic full waveform inversion to a low-frequency large-offset land data set[J].Expanded Abstracts of 80thAnnual Internat SEG Mtg,2010:930-934

[3] PRATT G R,GAO F,ZEL T C,et al.A comparison of ray-based and waveform tomography:implications for migration[R].Florence,Italy:64thEAGE Conference and Exhibition,2002

[4] PRATT G R.Velocity models from frequency-domain waveform tomography:past,present and future[R].Paris,France:66thEAGE Conference and Exhibition,2004

[5] GAO F C,LEVANDER A R,PRATT R G,et al.Waveform tomography at a ground water contamination site:VSP-surface dataset[J].Geophysics,2006,71(1):H1-H11

[6] SIRGUE L,DENEL B,GAO F.Integrating 3D full waveform inversion into depth imaging projects[J].Expanded Abstracts of 81stAnnual Internat SEG Mtg,2011:2354-2358

[7] TARANTOLA A.Inversion of seismic reflection data in the acoustic approximation[J].Geophysics,1984,49(8):1259-1266

[8] PRATT G R,SHIN C,HICKS G J.Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[J].Geophysical Journal International,1998,133(2):341-362

[9] HUGUES A D,TARANTOLA A.Multiparameter l1 norm waveform fitting:interpretation of Gulf of Mexico reflection seismograms[J].Geophysics,1999,64(4):1023-1035

[10] GUITTON A,SYMES W W.Robust inversion of seismic data using the Huber norm[J].Geophysics,2003,68(4):1310-1319

[11] KIM W,SHIN C.Phase inversion of seismic data with unknown source wavelet:synthetic examples[J].Expanded Abstracts of 75thAnnual Internat SEG Mtg,2005:1685-1689

[12] ARAVKIN A,VAN LEEUWEN T,HERMANN F.Robust full-waveform inversion using the student’s t-distribution[J].Expanded Abstracts of 81stAnnual Internat SEG Mtg,2011:2669-2673

[13] ROUTH P,KREBS J,LAZARATOS S,et al.Encoded simultaneous source full-wavefield inversion for spectrally shaped marine streamer data[J].Expanded Abstracts of 81stAnnual Internat SEG Mtg,2011:2433-2437

[14] HA W,SHIN C.Laplace-domain full-waveform inversion of seismic data lacking low frequency information[J].Geophysics,2012,77(5):R199-R206

[15] LUO J,WU R.Envelope inversion—some application issues[J].Expanded Abstracts of 83rdAnnual Internat SEG Mtg,2013:1042-1047

[16] BAEK H,CALANDRA H,DEMANET L.Velocity estimation via registration-guided least-squares inversion[J].Geophysics,2014,79(2):R79-R89

[17] SYMES W W,CARAZZONE J J.Velocity inversion by differential semblance optimization[J].Geophysics,1991,56(5):654-663

[18] PRATT G R,SYMES W W.Semblance and differential semblance optimization for waveform tomography:a frequency domain implementation[J].Journal of Conference Abstracts,2002,7(2):183-184

(编辑：戴春秋)

A new objective function for full waveform inversion:differential semblance optimization in data domain

GAO Fuchun1,Paul WILLIAMSON1,R. Gerhard PRATT2

(1.TotalE&PResearchandTechnology,Houston77002,USA;2.UniversityofWesternOntario,LondonN6A3K7,Canada)

Full Waveform Inversion (FWI),while now widely practiced industrially,is less robust than many conventional velocity model building techniques,such as travel time tomography,due to its high non-linearity.Different objective functions in FWI have different degrees of non-linearity.In this study,we investigate the behavior of FWI with different objective functions and propose a new objective function based on differential semblance defined in the data domain.Preliminary tests suggest that this objective function is convex for a large range of data residuals.Gradient-based optimization schemes are therefore more robust than for the standard least-squares formulation;however,the good resolving power of waveform inversion is mostly retained.

full waveform inversion,objective function,non-linearity,differential semblance optimization

2016-10-08;改回日期：2016-11-18。

高福春(1968—),男,道达尔勘探和生产研究技术中心地球物理研究员,研究方向为反演及其在地震波成像中的应用、地震波偏移等。

P631

A

1000-1441(2017)01-0026-05

10.3969/j.issn.1000-1441.2017.01.003