基于压缩感知的快速最小二乘逆时偏移

屠 宁

(同济大学海洋与地球科学学院,上海200092)

压缩感知理论的出现[6-7]为解决规模庞大的地震反演问题提供了新思路。在地震反演偏移成像问题中,对于一个压缩的物理观测模型,压缩感知理论提出了:①精确求解原始观测模型所对应解的充分条件,从而指导我们最优化设计压缩观测系统;②如何利用稀疏优化的方法从压缩观测的数据中恢复真实的地震像。本文将介绍作者与HERRMANN等在这个领域共同研究取得的一些重要成果[8-11]:①最小二乘逆时偏移的快速计算;②消除逆时偏移中未知的子波参数对反演结果的影响;③包含表面多次波的地震数据的成像。同时展示了合成和实际数据的实验结果,以验证本文方法的可行性和有效性。

1 方法和合成数据实验

1.1 快速最小二乘偏移

HERRMANN等[8]提出了基于压缩感知的最小二乘偏移方法,TU等[9]在此基础上,分析了波动方程线性化过程中引入的线性误差对该方法的影响。与传统的炮编码一样,基于压缩感知的最小二乘偏移也是通过对数据降维来实现波动方程模拟量的减少。与传统的最小二乘偏移不同,求解压缩感知问题的Basis Pursuit De-Noise(BPDN)模型[6]引入了稀疏约束,其在频率域求解的形式如下:

(1)

式中:Ω表示对数据降维时选择的频率子集;Σ表示选择的编码炮集;i和j分别表示频率和炮集的编号;d表示随机降采样的观测数据(在本文中,根据上下文可判断出是仅包含一次波的观测数据,还是同时包含一次波和多次波的观测数据);m0表示光滑的背景模型;wi表示地震子波w在某个频率上的频谱值;s表示不包含子波效应的随机降采样后的震源波场(即主要包含与观测数据相对应的震源的位置信息,子波效应包含在w中,降采样方式和观测数据d一致);C表示曲波变换(C*表示曲波反变换);x是模型扰动δm在曲波域的系数(δm=C*x);σ是可调节参数,用以控制匹配数据的精度,其选择往往反映数据的信噪比。(1)式的意义是,在所有可能的曲波系数中,寻找一个最稀疏的系数x,经过曲波反变换后得到模型扰动的一个估计,再经过反偏移,可以在由σ指定的误差范围内匹配观测数据d。引入曲波变换的意义在于,由于地质模型在空间域并非稀疏,曲波变换可将其稀疏化[12-13]。直观上解释,上述问题的含义是:在所有可能的曲波系数中,寻找最稀疏(一范数最小)的解,这个解经过反曲波变换(变换到空间域中),得到地下介质的扰动模型,该模型经过反偏移算子的作用,能在σ指定的误差范围内匹配观测到的地震数据。HERRMANN等[8]利用谱梯度投影算法(SPGL1)求解上述问题[14]。

对比传统的最小二乘偏移算法,压缩感知引入了稀疏约束(公式(1)第2行)。然而压缩感知对最小二乘偏移的影响还不止于此:压缩感知理论对数据如何进行降采样(同时对应对炮集如何进行编码)也有指导作用。压缩感知理论指出,在采样过程中引入随机因素,是最终恢复原始信号的关键条件。信号恢复的充分条件由Restricted Isometry Property(RIP)条件给出[6];在地震信号的处理中,HENNENFENT等[15]对随机采样的重要性亦作了直观的阐述。因此,在对观测数据进行降采样时,应随机选取频率的子集;在编码炮集时,也应引入足够的随机因素。详细的降采样方案请参考文献[8]和文献[11]。

另一个问题就是波动方程正演的线性化(即用反偏移代替正演)带来的误差。在数据观测较多,成像系统是超定系统时,最小二乘的方法可以自动地将正交于反偏移算子的列空间的噪声排除在解空间之外;然而,随着降采样带来的系统欠定性增强,噪声可能会在解空间投影,从而给真实解带来干扰。TU等[9]详细分析了线性误差对成像系统带来的影响,并提出,根据特定的规则在不同的迭代之间引入不同的降采样算子(即对数据和炮集重新进行采样)会显著提高成像系统对线性噪声的鲁棒性。

利用SEG/EAGE盐丘模型的一个二维剖面,对上述算法的有效性作了验证[10]。实验模型的参数如下:模型3.9km深,15.7km长,网格间距约为24m。共设置323个炮点和检波点,利用互易原理,将数据设置为每一炮都有共同的检波器来接收。真实速度模型及其平滑之后的背景模型如图1所示。图2展示了一系列实验结果。首先,利用线性数据(即利用反偏移算子作用到真实模型扰动上生成的地震数据)分析无任何噪声的状况下所提算法的有效性。图2a为对线性数据做单次逆时偏移的结果;图2b为对线性数据做快速最小二乘偏移的结果。在快速最小二乘偏移中,所用合成炮的数量约为原炮数的5%,所用频率子集的大小约为总频率集合的15%,我们运行60个SPGL1的迭代,在这样的反演参数设置下,总的波动方程的求解数量相当于利用所有数据进行单次逆时偏移,即图2a和图2b所展示的结果包含几乎相同的波动方程求解。接着,验证存在相干噪声的情况下该算法的效果。用时间域的iWave波场正演引擎来模拟观测波场[16],用我们的频率域波场线性反偏移引擎来进行成像,结果如图2c所示。这两种波场模拟引擎间存在数值算法上的差别,同时,iWave正演的观测波场含有高阶散射。在这种情况下,快速最小二乘成像方法依然取得了较好的效果。和图2a相比,图2b和图2c均展现出了较高的分辨率。

图1 快速最小二乘逆时偏移合成数据实验所用的真实速度模型(a)和光滑速度模型(b)

图2 线性数据的单次逆时偏移结果(a)和线性数据(b)、iWave正演数据(c)的快速最小二乘偏移结果

1.2 震源子波估计

无论是单次逆时偏移还是最小二乘逆时偏移,都需要预先知道地震子波。如果地震子波估算不准,那么正传波场就会有误,从而影响地震成像以及最小二乘偏移中梯度的求取。求取可靠的地震子波在技术上或者计算上面临较大的挑战:依赖统计信息[17]或者大量的计算[18]。最小二乘偏移本身是一种反演算法,因此我们将子波的反演融入其中,在不给出子波先验信息的情况下,得到可靠的地震成像结果。

借助在非线性最小二乘问题中变量投影(Variable Projection)的方法[19-20],并将其理论推广到包含稀疏约束的非线性最小二乘问题中,提出了不依赖于子波信息的最小二乘逆时偏移算法[10]。对于子波w来说,一旦x确定,存在下述解析表达:

(2)

使用交替求解的方式,每更新一次x,就更新一次w的估计,即可在未给定子波的情况下求解最小二乘偏移的梯度[10]。

继续利用前文所用的模型进行合成数据的实验,主要结果如图3所示。图3a展示了利用线性化数据进行最小二乘逆时偏移的结果,但正演中使用的地震子波和真实的子波存在一个相位差,其影响是反射层发生错位(如箭头所示)和整体地震像散焦。图3b展示了利用相同的线性化数据,在未给定子波(第一个迭代中的子波起始估计是白谱子波，即时间域的脉冲)的情况下,利用本文包含子波估计的成像方法得到的地震像。图3c展示了利用iWave正演的数据进行实验的结果(其它参数与图3b所用的一致)。从结果来看,本文方法可以在未知真实子波的情况下,得到可靠的地震像。对比图3b,图3c和图2b,图2c可以看出,采用本文方法得到的地震像与知道真实子波情况下得到的地震像具有相似的高精度和分辨率。

1.3 表面多次波成像

海洋表面多次波是海洋地震资料处理时面临的强干扰。逆时偏移本身基于单次散射的假设,并不能对多次波进行有效成像。我们从描述自由表面反射的自由表面多次波消除公式出发,研究出对多次波进行成像以及一次波和多次波联合成像的方法[10-11]。海洋表面多次波的正确成像依赖于:①正确识别表面多次波对应的震源波场。一次波的成像中,正传波场为点震源波场;多次波的成像中,正传波场是海洋表面的整体下行波场,可以看作一个面震源[11,21-23]。②新的成像条件。对于含有多阶海洋表面反射的多次波来说,广泛应用于一次波成像的互相关成像条件会造成各阶次反射之间产生非物理的串扰噪声[24-25],形成对正确地震像的干扰。应用反褶积成像条件可以一定程度缓解这个问题[24,26],而真正衰减串扰噪声则需要利用最小二乘反演的方法[11,27]。

基于我们提出的快速最小二乘逆时偏移的框架,实现对表面多次波的正确成像。通过对正传震源波场做修正,用自由表面的下行波场代替点震源波场,可实现对表面多次波的正确成像[11]以及一次波和多次波的联合成像[10,28]。数据的预处理细节参考文献[10]。本文主要介绍一次波和多次波联合成像的方法。这时,数据匹配的目标变为:

(3)

需要注意的是,此处的观测数据d同时包含了上行的一次波波场和多次波波场。对于预测数据,可以看到,在反偏移算子中,震源波场wisj-di,j不只是包含点震源s(经过波场传播后预测一次波),同时包含广义的面震源波场即观测数据d(前面的“-”表示水面的反射系数为-1,上行的观测数据经过水面反射,再经过和一次波相同的波场传播后,预测表面多次波)。除去震源波场的变化,其它成像过程与前文所述一次波的反演成像过程相同。另外,这种情况下,由于震源波场的变化,子波的估计方法则为:

(4)

选用专为验证与表面多次波有关的算法而设计的Sigsbee 2B模型,对包含海洋表面多次波的数据进行了测试。真实和背景速度模型见图4。数据正演使用iWave,使用自由表面边界条件来正演包含多次波和鬼波的观测数据。假设真实地震子波已知,利用本文算法对一次波和多次波进行联合快速最小二乘偏移成像,结果如图5a所示,其整体的波动方程模拟次数依然控制在单次逆时偏移的水平上。图5b 展示了通过子波估计对一次波和多次波进行联合成像的结果。可见,图5a和图5b结果具有很高的一致性,与真实模型对比可以看出,这两个结果中多次波带来的串扰噪声均得到了很好的压制。文献[10]中详细讨论了本文方法对子波初始估计的鲁棒性,本文方法在子波具有错误的初始振幅和相位的情况下,均能得到高精度的地震像。

图4 包含多次波数据偏移实验所用真实的Sigsbee 2B速度模型(a)和光滑后的速度模型(b)

图5 使用真实子波(a)和子波估计方法(b)得到的一次波与多次波联合快速最小二乘偏移成像结果

2 实际数据实验

利用英国北海海域某工区实际数据,对文中所述算法进行验证。图6a为利用互相关成像条件对多次波进行偏移成像的结果(成像过程中,震源为面震源,由表面位置的下行波场构成,反传数据为多次波),可

以看到不同阶次多次波导致的串扰噪声。图6b为利用本文快速最小二乘偏移算法对多次波成像的结果。图7和图8分别为图6中区域A和区域B放大后的结果。可以看出,图6a中出现的极强的串扰噪声在图6b中得到了很好的压制。

图9为不同数据的偏移成像结果。图9a是一次波成像结果;图9b是多次波单独成像结果;图9c是一次波与多次波联合成像的结果。可以看出,多次波单独成像和一次波成像存在较好的可比性,较强的连续反射层(红色框中的反射层)均能得到较高精度的成像,但是整体信噪比较低,有些反射层(黑色箭头所指)没有成像;对于某些一次波成像没能揭示出的结构(白色箭头所指),多次波成像能较好的刻画出来。从图9c可以看出,一次波与多次波的联合成像,不仅在数据处理环节避免了多次波和一次波的分离,而且改进了多次波成像本身信噪比偏低的缺点。

图6 利用互相关成像条件(a)和快速最小二乘逆时偏移(b)对多次波进行成像的结果

图7 图6a(a)和图6b(b)中区域A放大后的结果

图8 图6a(a)和图6b(b)中区域B放大后的结果

图9 不同数据的偏移成像结果a 仅用一次波; b 仅用多次波;c 一次波和多次波联合成像

3 结论

本文将压缩感知应用于地震成像中,通过人工压缩观测系统和观测数据,大幅降低最小二乘逆时偏移的成本。快速最小二乘逆时偏移是个灵活的框架,在这个框架下,介绍了如何利用变量投影技术解决地震偏移中未知子波的问题,并利用广义面震源实现了海洋多次波的成像以及一次波与多次波的联合成像。模型数据以及海上实际数据应用结果表明,本文方法在未知地震子波信息的情况下,利用相当于单次逆时偏移的计算量,可以得到高分辨率的最小二乘偏移成像结果,并对多次波进行正确的归位和成像,为今后更加有效和多样化地利用多次波信息提供了思路。

[1] LAILLY P.The seismic inverse problem as a sequence of before stack migrations[C]∥BEDNAR J B.Conference on inverse scattering:theory and application,Philadelphia:SIAM,1983:206-220

[2] NEMETH T,WU C,SCHUSTER G T.Least-squares migration of incomplete reflection data[J].Geophysics,1999,64(1):208-221

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

[4] ROMERO L,GHIGLIA D,OBER C,et al.Phase encoding of shot records in prestack migration[J].Geophysics,2000,65(2):426-436

[5] SCHUSTER G T,WANG X,HUANG Y,et al.Theory of multisource crosstalk reduction by phase-encoded statics[J].Geophysical Journal International,2011,184(3):1289-1303

[6] DONOHO D.Compressed sensing[J].IEEE Transactions on Information Theory,2006,52(4):1289-1306

[8] HERRMANN F J,LI X.Efficient least-squares imaging with sparsity promotion and compressive sensing[J].Geophysical Prospecting,2012,60(4):696-712

[9] TU N,LI X,HERRMANN F J.Controlling linearization errors in l_1 regularized inversion by rerandomization[J].Expanded Abstracts of 83rdAnnual Internat SEG Mtg,2013:4640-4644

[10] TU N,ARAVKIN A Y,LEEUWEN T,et al.Source estimation with surface-related multiples-fast ambiguity-resolved seismic imaging[J].Geophysical Journal International,2016,205(3):1492-1511

[11] TU N,HERRMANN F J.Fast imaging with surface-related multiples by sparse inversion[J].Geophysical Journal International,2015,201(1):304-317

[13] HERRMANN F J,MOGHADDAM P P,STOLK C.Sparsity- and continuity-promoting seismic image recovery with curvelet frames[J].Applied and Computational Harmonic Analysis,2008,24(2):150-173

[14] VAN DEN BERG E,FRIEDLANDER M P.Probing the Pareto frontier for basis pursuit solutions[J].SIAM Journal on Scientific Computing,2008,31(2):890-912

[15] HENNENFENT G,HERRMANN F J.Simply denoise:wavefield reconstruction via jittered undersampling[J].Geophysics,2008,73(3):V19-V28

[16] TERENTYEV I S,VDOVINA T,SYMES W W,et al.iWave:a framework for wave simulation[EB/OL].[2017-01-10].http:∥www.trip.caam.rice.edu/software/iwave/doc/html/index.html

[17] ULRYCH T J,VELIS D R,SACCHI M D.Wavelet estimation revisited[J].The Leading Edge,1995,14(11):1139-1143

[18] VIRIEUX J,OPERTO S.An overview of full-waveform inversion in exploration geophysics[J].Geophysics,2009,74(6):WCC1-WCC26

[19] GOLUB G,PEREYRA V.The differentiation of pseudo-inverses and nonlinear least squares problems whose variables separate[J].SIAM Journal on Numerical Analysis,1973,10(2):413-432

[20] GOLUB G,PEREYRA V.Separable nonlinear least squares:the variable projection method and its applications[J].Inverse Problems,2003,19(2):R1-R26

[21] BERKHOUT A J.Migration of multiple reflections[J].Expanded Abstracts of 63rdAnnual Internat SEG Mtg,1993:1022-1025

[22] GUITTON A.Shot-profile migration of multiple reflections[J].Expanded Abstracts of 72ndAnnual Internat SEG Mtg,2002:1296-1299

[23] ZUBERI A,Alkhalifah T.Imaging by forward propagating the data:theory and application[J].Geophysical Prospecting,2013,61(S1):248-267

[24] MUIJS R,ROBERTSSON J O A,HOLLIGER K.Prestack depth migration of primary and surface-related multiple reflections:part Ⅰ—imaging[J].Geophysics,2007,72(2):S59-S69

[25] LIU Y,CHANG X,JIN D,et al.Reverse time migration of multiples for subsalt imaging[J].Geophysics,2011,76(5):WB209-WB216

[26] WHITMORE N D,VALENCIANO A A,SOLLNER W.Imaging of primaries and multiples using a dual-sensor towed streamer[J].Expanded Abstracts of 80thAnnual Internat SEG Mtg,2010:3187-3192

[27] WONG M,BIONDI B,RONEN S.Imaging with multiples using least-squares reverse time migration[J].The Leading Edge,2014,33(9):970-976

[28] TU N,HERRMANN F J.Fast least-squares imaging with surface-related multiples:application to a North-Sea data set[J].The Leading Edge,2015,34(7):788-794