2018, Vol. 61
2. 美国地球物理先进技术公司, 休斯顿 77478
2. Advanced Geophysical Technology Inc., Houston 77478, U. S. A
地震勘探中经常遇到有强衰减介质(介质品质因子Q值较低)的区域(Müller et al., 2010).地震波在其中传播会耗散大量能量.反映在地震记录上，表现为振幅下降和波形畸变.如果直接利用这些地震记录进行逆时偏移成像，会降低成像质量(Yu et al., 2002; Zhou et al., 2011)，具体表现为衰减区域及其下方照明严重不足，分辨率降低，偏移深度错误(Zhang et al., 2010).因此，在偏移成像过程中，对强衰减介质区域进行衰减补偿是十分必要的.
为了提高强衰减介质区域的成像质量和分辨率，前人发展了诸多在偏移成像中进行衰减补偿的方法.早期反Q滤波方法(Bickel and Natarajan, 1985; Hargreaves and Calvert, 1991)被广泛采用.但是此方法只适合于一维的Q值模型，难以处理Q值有横向变化的情况.基于射线理论的叠前深度偏移衰减补偿方法(Zhang et al., 2002, 2013)、变换域衰减补偿方法(刘喜武等，2006)和带衰减补偿的Kirchhoff偏移(Traynin et al., 2008)都可以考虑横向变化的Q，可以更加精确地补偿衰减效应.但是受限于高频近似，在复杂介质尤其是有多散射体或尖锐反射面存在的情况下会导致成像精度下降.基于单程波动方程的衰减补偿偏移(Dai and West, 1994; Yu et al., 2002；崔建军和何继善，2001；李金丽等，2015；白敏等，2016)能处理更复杂的速度模型包括倾斜地层、断层等强散射体.然而在地下速度结构具有强烈的横向非均匀性时，也难于得到高分辨率图像.带衰减补偿的最小平方偏移方法(Dutta and Schuster, 2014; Sun et al., 2015a; 邓文志等，2015)适用于更复杂的地下介质，但计算量非常大.
由于逆时偏移(RTM)技术(Baysal et al., 1983)利用双程波动方程建立了全波场正演模拟，对复杂地质结构具有良好适应性，近年来得到广泛应用(Etgen et al., 2009).同时，随着Q值层析成像方法的逐步发展(Xin et al., 2010; Clark et al., 2009; Hu et al., 2011; Cavalca et al., 2011; Shen and Zhu, 2015; Dutta and Schuster, 2016；张固澜等，2014)，人们可以得到更加精确的Q值模型，使得在逆时偏移中进行衰减补偿(Q-RTM)成为可能.不过，虽然物理衰减的正演模拟很早就被发展起来(Carcione et al., 1988; Robertsson et al., 1994)，Q-RTM中关于衰减的补偿仍然具有一定难度.一方面，衰减补偿是一个将波场能量指数放大的过程，不容易稳定；另一方面，为了同时纠正物理衰减带来的照明不足和波形畸变问题，我们必须同时补偿振幅和相位(Zhou et al., 2018).
随着Q-RTM技术的发展，出现了两大类时间域的衰减补偿方法.其一是基于分数阶微积分的常数Q模型(FDCQ)衰减补偿方法(Zhang et al., 2010; Zhu et al., 2014a, 2014b, 2017; Sun et al., 2015b; 吴玉等，2017)，其二是基于广义标准线性体模型(GSLS)的衰减补偿方法(Causse and Ursin, 2000; Deng and McMechan, 2007; Deng et al., 2008; Zhou et al., 2018).FDCQ方法利用分数阶微积分本构关系得到了常数Q(不随频率变化)的模型，并实现了衰减部分和频散部分的分解(Zhu et al., 2014a)，能够在保持相速度不变的情况下将Q反号，从而实现同时补偿振幅和相位.但由于该方法使用基于分数阶微分算子描述的本构关系，只能用全局求解波场的伪谱法实现，导致FDCQ方法计算费用较高.同时也难以并行化，难以利用GPU进行加速.所以，在三维大区域计算时，全局方法不具有优势(Ayala and Wang, 2013).
基于GSLS的衰减补偿方法利用数个标准线性体串联起来，在一定频率范围内拟合一个常数Q(Blanch et al., 1995)，可以利用具有局域性的有限差分法实现，即数值上局域地求解波场，只用相邻几个点的值计算导数，易于灵活地进行计算区域分解从而易于细粒度并行和GPU加速(Abubakar et al., 2009, 2011; Abdelkhalek et al., 2012).但进行振幅和频散的分离是其难点，无法在保持相速度频散关系不变的情况下补偿振幅(赵岩等，2014; Guo et al., 2016).早期基于GSLS的Q-RTM方法(Causse and Ursin, 2000; Deng and McMechan, 2007; Deng et al., 2008)只做到了补偿振幅，而不能保证相速度频散关系不变.为了克服该方法的这个弱点，我们之前提出了一种基于多级优化来修正相速度频散关系的方法(Hu et al., 2016; Zhou et al., 2018).利用多项式优化广义标准线性体的复模量实部.这样在保持微分算子局域性的基础上，既可以补偿振幅，又可以在一定频率范围内保持相速度频散关系基本不变.我们称之为MORC方法(Multi-stage Optimization on Real part of Complex moduli).这个方法使得利用有限差分实现准确的Q-RTM成为了可能.但是，MORC方法在低频范围(5~10 Hz)精度不高.此时由于相速度与真实模型的误差较大，会发生一定的波形的畸变和成像深度的偏差.
本文提出一种新的衰减补偿方法，即保持广义标准线性体的复模量实部不变，通过多项式优化复模量的虚部，在保持正确的相速度不变的同时实现振幅的补偿.通过合理设计优化项，我们在更宽的频率范围(5~70 Hz)内实现了准确的振幅补偿.由于同样能够应用有限差分进行计算，本文方法同样适合进行高效的并行计算，也可以容易地与现有的RTM方法进行合成，不需要进行过多的代码重构(Hu et al., 2016; Zhou et al., 2018).
1 基于线性黏弹性体复模量虚部拟合的声波振幅衰减补偿 1.1 基于广义标准线性体(GSLS)的常数Q近似方法在勘探地球物理中涉及的频率范围内，常数Q模型(不随频率变化)是对弹性能衰减的合理表达(Kolsky, 1956; Lomnitz, 1957; Kjartansson, 1979).其中既有物理衰减的贡献，也有散射的贡献.因为二者均表现为主要地震震相弹性能的损失，所以形式上可以归于物理衰减.Robertsson等(1994)利用几个标准线性体串联形成的广义标准线性体(GSLS)来模拟一定频率范围内的常数Q.这一方法在国内外波动方程模拟和偏移工作中得到了广泛应用(Deng and McMechan, 2007; 徐文才等，2016).由Blanch等(1995)进一步给出简化的黏声波波动方程如下:
|
(1) |
|
(2) |
|
(3) |
这是一组速度-应力形式的波动方程.其中，p是应力，v是质点速度, MR是弹性模量, ρ是密度，l是串联的标准线性体个数, τσl是每个标准线性体的应力弛豫时间.在近似条件τ≪1下(Q＞10)，τ是一个只与Q相关的无量纲参数.rl是记忆变量.记忆变量的作用是避免直接求解GSLS模型本构方程中的卷积(Carcione et al., 1988).
Q值和相速度随频率的变化关系如下(Robertsson et al., 1994; Blanch et al., 1995)：
|
(4) |
|
(5) |
其中，M(ω)为作为角频率ω函数的复模量.
在一定频率范围内近似常数的Q是通过几个标准线性体叠加形成的.在本文的研究中，我们选用5个标准线性体来得到近似常数Q，拟合方法同Blanch等(1995).图 1中实线所示的是一个Q=60的例子，包括Q值随频率变化和相速度频散关系.
|
图 1 GSLS模型和负τ方法的Q值和相速度频散关系 (a) Q值随频率变化关系.实线为GSLS模型拟合的常数Q值.虚线为理想常数Q值(Q=60)，点线为负τ方法的Q值; (b)相速度频散关系.实线为GSLS模型的相速度频散，点线为负τ方法的相速度频散. Fig. 1 Relation between Q value and normalized phase velocity dispersion for the GSLS model and negative τ-method (a) Q value varying with frequency. Solid line is GSLS model, dashed line is ideal constant Q (Q=60), dotted line is negative τ method. (b) Phase velocity dispersion curve. Solid line is GSLS model, dotted line is negative τ method. |
类似于RTM，Q-RTM也是一个三步成像过程.第一步是得到震源侧的波场；第二步是将数据时间域反转求解波动方程，得到接收器侧的反传波场；最后利用互相关成像条件得到偏移结果(Claerbout，1971).震源侧波场求解以及接收器侧波场求解都需要在求解波动方程的时候对衰减进行补偿.为了同时补偿振幅衰减和相位扰动，要求新的波动方程(Q补偿方程)对应的Q值应等于GSLS波动方程Q值的相反数，而其对应的相速度频散关系与GSLS波动方程相速度频散关系一致(Hu et al., 2016; Zhou et al., 2018).
修改方程组(1)—(3)建立Q补偿方程可以自然保留局域化算子.根据公式(4)，直接更改τσl的符号似乎可以同时满足Q补偿方程对应的Q值反号和相速度不变.然而这种方法在数值上不稳定(Xie et al., 2015).在数值测试中，这种方法甚至在计算开始的最初几步就出现了不稳定现象.因此，早期的Q-RTM工作(Causse and Ursin, 2000; Deng et al., 2007, 2008)采用更改τ或记忆变量rl的符号的方法补偿振幅.为了表述方便，我们简称这类方法为负τ方法.负τ方法在使Q反号(并非仅仅反号，Q的绝对值也发生了变化，如图 1a中点线所示)的同时，也改变了相速度.如图 1b中点线所示，这类方法使相速度频散关系发生了大的变化，导致无法补偿由衰减造成的相位偏离，甚至使之加剧.负τ方法的成像结果中会出现偏移深度错误，波形畸变等严重问题(Guo et al., 2016; Zhou et al., 2018).
在负τ方法中，根据公式(5)，如果将参数τ反号，就会严重改变相速度频散关系.这是由于负的τ值导致了复模量的实部发生改变.我们观察到，如果不改变τ的符号，则相速度就不会被改变，这对相位补偿是很好的性质.本文提出利用多项式多级优化方法去调整复模量的虚部，使Q值在一定频率范围内反号，实现了对振幅的补偿.
根据波动方程(1)—(3)，在GSLS模型中，复模量为：
|
(6) |
在公式(6)的虚部上加入多项式补偿项：
|
(7) |
其中, αnω2n-1是多项式优化项，n=0, 1, 2, ….由于在变换域中，(iω)n对应时间域中的n阶时间偏导数∂n/∂tn，而整数阶导数可以用局域方法求解，故多项式的补偿项可以保持微分算子的局域性.补偿项只加在公式(6)的虚部上，所以我们只利用ω的奇数次多项式来进行复模量虚部的优化补偿，称之为n级优化.这是关于频率ω的2n-1次多项式.特别地，n=0时即零级优化时，优化项是关于频率ω的一次分式α0ω-1.
MOIC方法可以表示为如下的多项式优化问题：
|
(8) |
其中，“+”号表示优化的目标为GSLS模型复模量虚部Im(M(ω))的相反数.由于Q值可以写成复模量的实部除以虚部，而复模量的实部不变.为了方便表示，我们考察Q值的倒数.1/Q值直接表示了一个周期内衰减(补偿)的大小.
|
(9) |
其中分母为复模量的实部，在优化过程中保持不变.分子为(8)式求解出来的优化虚部.因此，对复模量虚部的优化相当于优化1/Q.根据数值测试结果，类似于一般的多项式拟合.如图 2所示.
|
图 2 2-4级优化虚部修正方法的1/Q曲线 (a) 2级优化；(b) 3级优化；(c) 4级优化；(d) 4级优化但无分式优化项.各图中实线表示理想的1/Q值随频率变化关系，虚线表示优化的1/Q值随频率变化关系. Fig. 2 2-4 stage 1/Q optimization result for MOIC method Solid line is the ideal 1/Q value, dashed line is the optimized 1/Q value. (a) 2-stage; (b) 3-stage; (c) 4-stage; (d) 4-stage without fractional optimization term. |
图 2显示的是Q=60情况下从2级优化到4级优化的1/Q随频率的关系.可以看出，从2级优化到4级优化的过程中，能拟合的频带范围不断加宽，并且拟合程度逐渐提升，误差逐渐减小.2级优化可以较好地拟合5~25 Hz，3级优化能较好地拟合5~40 Hz，而4级优化可以在5~70 Hz的频率范围内对1/Q实现15%以内误差的逼近.与MORC方法(Hu et al., 2016; Zhou et al., 2018)相比，虽然MOIC方法在振幅补偿上有更大的误差，但其能在保持全频带正确的相速度不变的同时实现振幅的补偿，并能拟合更宽的频带范围，从而避免了MORC方法在低频范围(5~10 Hz)波形畸变和成像深度偏差大的问题.
另外，分式优化项α0ω-1很重要.比较图 2c—2d，可以发现，MOIC方法低频适应性很好的主要原因就是加入了分式优化项α0ω-1.分式的加入可以将低频段1/Q频散曲线的趋势调整过来，对多项式优化十分有利.如果不加分式优化，低频段将很难调整成负值.而且分式优化项也可以通过记忆变量的形式转化为整数阶偏微分算子，保持整个系统偏微分算子的局域性.基于这些良好的性质，我们最终实现了在低频到高频的一个较宽频带范围内，振幅和相位能同时准确补偿的效果.
下面推导MOIC方法的波动方程.首先给出复速度：
|
(10) |
然后用压强p和质点速度v改写方程(10)，使之变成频率域速度-应力形式的方程：
|
(11) |
|
(12) |
然后定义记忆变量：
|
(13) |
其中rl是对应于标准线性体的记忆变量，而s是对应于分式优化项的记忆变量.这里我们也借助记忆变量将分式优化项转化成独立的微分方程，避免了直接变换回时间域带来的卷积.用公式(13)替换公式(11)中相关的项，剩余的优化项我们将其归入速度中.最后将替换后的公式(11)和(12)都反变换回时间域，就得到MOIC方法的速度-应力形式波动方程：
|
(14) |
|
(15) |
|
(16) |
|
(17) |
|
(18) |
方程组(14)—(18)是MOIC方法的Q补偿波动方程.求解(14)—(18)可以获得一个衰减补偿的波场.它既可以精确地补偿相位，也可以在一个设定的频率范围内(如4级优化：5~70 Hz)准确地补偿振幅的衰减.方程(16)是优化补偿项.这一项随着优化级数的不同，含有不同的高阶时间导数.一般而言，计算高阶时间导数如用高阶算法保证精度，则需要存储额外的历史波场值，增大计算量及内存开销.为了数值计算的方便，我们可以近似地利用空间导数替代时间导数.由式(15)知道，在τ≪1时记忆变量rl也远小于1，并且由优化结果得知优化项同样远小于1.故空间导数和时间导数近似地具有如下关系：
|
(19) |
公式(19)的近似条件τ≪1也是GSLS模型的近似条件.一般情况下，Q＞10，满足这一近似条件(Blanch et al., 1995).
2 数值测试 2.1 MOIC方法的振幅补偿和波形恢复能力测试
采用时间2阶精度，空间8阶精度的有限差分方法求解方程组(14)—(18).方程组(14)—(18)与GSLS方程组(1)—(3)有相同的稳定性条件(Robertsson et al., 1994)，即cmΔt/(Δx, Δy, Δz)min＜1，其中
|
表 1 Q=60时4级优化参数表 Table 1 4-stage optimization parameters for Q=60 |
利用一系列数值模型验证了本文方法的有效性和准确性.首先，为了分析MOIC方法对不同主频的信号之振幅补偿和波形恢复的效果，我们对一个Q=60，cp=2200 m·s-1的均匀模型进行一维测试，使用主频分别为10、25 Hz的Ricker子波分别传播4000、2000 m到达接收点，随后再将接收点记录到的波形时间上反转，作为震源输入，求解方程组(14)—(18)，传播相同的距离到达震源处，得到补偿结果，如图 3和图 4所示.
|
图 3 10 Hz波形补偿测试 (a) MOIC补偿结果和参考波形对比；(b) MORC补偿结果和参考波形对比；以上两图中实线表示理想波形，虚线和点线分别表示通过MOIC和MORC方法补偿过的波形；(c)两种方法的频谱.灰色实线表示理想频谱，点虚线为MORC补偿的频谱，虚线为MOIC补偿的频谱，灰色点线为不补偿的频谱. Fig. 3 Waveform recovery test for a 10 Hz Ricker wavelet (a) MOIC method; (b) MORC method; In (a) (b), solid, dashed, and dotted lines are respectively the ideal, MOIC and MORC compensated waveforms; (c) Amplitude spectrum of the two methods. Solid line in light gray is the ideal spectrum, dot-dashed line is the compensated spectrum for MORC method, dashed line is the compensated spectrum for MOIC method, and gray dotted line is the uncompensated spectrum. |
|
图 4 25 Hz波形补偿测试 (a) MOIC补偿结果和参考波形对比；(b) MORC补偿结果和参考波形对比；以上两图中实线表示理想波形，点线表示补偿过的波形; (c)两种方法的频谱.灰色实线表示理想频谱，点虚线为MORC补偿的频谱，虚线为MOIC补偿的频谱，灰色点线为不补偿的频谱. Fig. 4 Waveform recovery test for a 25 Hz Ricker wavelet (a) MOIC method; (b) MORC method; In (a) (b), solid and dotted lines are respectively the ideal and compensated waveforms; (c) Amplitude spectrum of two methods. Solid line in light gray is the ideal spectrum, dot-dashed line is the compensated spectrum for MORC method, dashed line is the compensated spectrum for MOIC method, gray dotted line is the attenuated spectrum. |
比较图 3a和图 3b，我们新提出的复模量虚部优化(MOIC)方法对于主频10 Hz的低频信号可以很好地恢复衰减.补偿后的波形与参考的Ricker子波十分接近.然而之前提出的MORC方法由于补偿后的相速度和理想相速度差距过大，导致波形的对称性无法恢复.图 3c显示了两种方法的频谱.可见MORC方法对频谱的恢复是比较好的.因为MORC方法在Q值拟合上更加精确.MOIC方法对频谱的恢复不如MORC方法，主要原因是Q值拟合的误差大于MORC方法.但是由于MORC方法对低频信号波形恢复不佳，在低频分量占主导时只能使用MOIC方法.
在主频为25 Hz的补偿测试中(如图 4所示)，可以看到Ricker子波的频率范围更宽，覆盖了0~70 Hz.我们新提出的MOIC方法同样能很好地恢复波形，并有良好的对称性.如图 4a所示.但补偿的波形出现了一些旁瓣假象.由图 4c可以发现这些旁瓣假象来源于不同频率补偿量的不同.回顾图 2c，我们发现在25 Hz附近，1/Q值约小10%，而在45 Hz附近，1/Q值约大10%.这对应了图 4c中25 Hz频率附近的过补偿和45 Hz附近的欠补偿.反应在波形上就是旁瓣假象.但由于MOIC方法不修改相速度频散关系，相速度是准确的.所以补偿后的波形仍然具有非常好的对称性.图 4b是MORC方法的补偿结果.可以看到它是没有旁瓣假象的，但波形的前部有一些高频假象，造成波形对称性稍微弱于MOIC方法.这些高频假象是由于40~50 Hz频率范围内的相速度偏快的拟合误差造成高频分量传播稍快，堆积在波形前部造成的(Zhou et al., 2018).但由于相速度拟合误差在高频区间十分小，所以在高频范围内MORC方法也能够获得比较准确的相位补偿结果.图 4c则显示了MORC方法能获得精确的振幅谱恢复效果.MOIC方法的优势主要在于相位补偿更加精确.由此带来更准确的偏移深度恢复和波形形状.
2.2 MOIC方法的Q-RTM成像测试我们利用MOIC的Q补偿方程建立了Q-RTM方法，并首先在一个简单的三层模型上进行成像测试.其中第一层为一水平层，第二、三层为倾斜层.模型的大小为3200 m×3200 m，网格间距为8 m，总共400×400个网格.区域四周加入了近似完全匹配层(NPML, Nearly Perfect Matched Layer)边界条件(Hu and Cummer, 2004; Hu et al., 2007)以消除边界反射波.模型参数如图 5所示.选用10 Hz的Ricker子波，以突出MOIC方法在低频范围内成像上的有效性.
|
图 5 三层Q-RTM模型示意图 Fig. 5 Sketch of three-layer Q-RTM model |
我们给出了当数据有衰减时，用传统的RTM，以及MOIC Q-RTM及MORC Q-RTM等三种方法进行成像的结果，同时给出了数据无衰减时的RTM图像作为参考，结果如图 6所示.图 6a是根据衰减资料得到的RTM结果.可以看到反射面的像不清晰，照明不足且较实际位置偏上(是由于相位没得到补偿导致的).图 6b和图 6c分别是MOIC Q-RTM和MORC Q-RTM，分别为虚部优化和实部优化.图 6d是根据未衰减资料得到的RTM作为参考.对比图 6b，6c可以看出，MOIC可以很好地补偿振幅和相位，大幅度提高成像的清晰度.并且，偏小的偏移深度也归位了.虽然看上去MORC和MOIC的补偿效果接近，但MOIC更接近于参考RTM图像.图 6e是水平距离x=1600 m处的垂直剖面.比较四种图像可以发现，确实在低主频情况下，MOIC比MORC更接近参考模型.此时，MOIC几乎完全与参考模型一致，波形的恢复相当好，并且没有偏移深度误差.而MORC发生了较明显的波形畸变和补偿不足现象，并且在2000 m深度上出现偏移深度偏差.这是低频范围实部优化的相速度误差过大造成的.
|
图 6 三层模型Q-RTM测试，使用10 Hz Ricker子波 (a)衰减数据的RTM结果；(b)衰减数据的MOIC Q-RTM结果；(c)衰减数据的MORC Q-RTM结果；(d)无衰减数据RTM结果(参考)；(e)四种图的x=1600 m处垂直道比较，深灰色实线为无衰减数据RTM结果，作为参考.黑色虚线为MOIC Q-RTM结果，灰色点线为MORC Q-RTM结果，浅灰色点虚线为衰减数据的RTM结果(无任何补偿). Fig. 6 Q-RTM tests on 3-layer model with 10 Hz Ricker wavelets (a) Conventional RTM image with attenuated data; (b) Q-RTM image with MOIC method; (c) Q-RTM image with MORC method; (d) Normal RTM with unattenuated data; (e) Vertical traces at x=1600 m. Solid line in dark gray is the RTM result with unattenuated data (reference), dashed line in black is MOIC Q-RTM result, dotted line in gray is MORC Q-RTM result, and dot-dashed line in light gray is RTM with attenuated data (uncompensated). |
为了检验基于MOIC的Q-RTM方法在接近实际的地质模型中的实用性，我们利用Sigsbee2A模型进行了进一步测试.模型大小为10000 m×3750 m，网格间距为6.25 m，共使用1600×600个网格.但区域边界50个网格为NPML边界条件区域.速度模型和Q模型如图 7所示.震源采用20 Hz Ricker子波，并在地表以50 m等间隔分布；接收器在地表以25 m等间隔分布.
|
图 7 Sigsbee 2A模型 (a)速度模型；(b) Q模型. Fig. 7 Sigsbee 2A model (a) Velocity model; (b) Q model. |
我们分别计算了衰减数据利用MOIC和MORC方法做Q补偿得到的Q-RTM结果.作为对比，还展示了衰减数据传统RTM(不进行Q补偿)，无衰减数据RTM以及基于伪谱法的Q-RTM(Zhu et al., 2014b)的结果.如图 8所示.
|
图 8 Sigsbee 2A模型的RTM结果 (a)衰减数据RTM(无补偿)；(b)无衰减数据RTM(参考)；(c) MORC方法RTM；(d) MOIC方法RTM；(e)伪谱法RTM(Zhu et al., 2014b)；(f)在x=2187 m处垂直剖面的比较.浅灰色实线为衰减数据RTM(无补偿)，深灰色实线为正常无衰减数据RTM，灰色点线为MORC Q-RTM，黑色虚线为MOIC Q-RTM，深灰色粗点线为伪谱法Q-RTM. Fig. 8 RTM results of Sigsbee 2A model (a) Conventional RTM with attenuated data (without compensation); (b) Normal acoustic RTM (without attenuation); (c) Q-RTM with MORC method; (d) Q-RTM with MOIC method; (e) Q-RTM with pseudo-spectral method (Zhu et al., 2014b); (f) A trace comparison of (a)—(e) at x=2187 m. Solid line in light gray is RTM with attenuated data (uncompensated), solid line in dark gray is RTM with unattenuated data, dotted line in gray is MORC Q-RTM, dashed line in black is MOIC Q-RTM, and thick dotted line in dark gray is pseudo-spectral Q-RTM. |
图 8a是用RTM做衰减数据(不进行Q补偿)的结果，图 8b和8c分别是MOIC Q-RTM方法和MORC Q-RTM方法得到的偏移成像结果.可以看出两种方法都很好地补偿了强衰减区域带来的分辨率下降和照明减弱问题.尤其是在Q区域下方，图 8a几乎没有成像.图 8b和8c都恢复出了强衰减区域下方的结构.图 8d是没有衰减区域的传统RTM图像，作为对比.可以看出图 8b—8d都十分接近，基本上完整恢复了由于衰减而缺失的信息，并没有增加额外的假象.图 8e是伪谱法(Zhu et al., 2014b)的Q-RTM图像.对比图 8b—8e, 这几种Q-RTM方法和参考图都非常接近，我们用较少的计算量达到了伪谱法的计算精度.图 8f是各个方法在x=2187 m处的垂直剖面比较.在垂直道上可以更明显地看出这些方法都很接近参考图像，MOIC方法稍微有些过补偿，但波形的对称性稍好一些.通过这个复杂模型我们进一步验证了MOIC方法在更宽频带范围内的补偿有效性.
3 讨论本文提出的MOIC方法通过多项式优化复模量的虚部，将Q值强行修正到负Q，避免了相速度的改变.但相应的，振幅补偿精度就会下降.在振幅方面，两种方法的Q值拟合误差起决定性作用.在5~70 Hz频率范围内，MORC方法的Q值拟合误差随频率变化约为0.5%~1%，这种误差一般对波形补偿图像和RTM图像都只产生非常小的影响.而MOIC方法的Q值拟合误差约为10%~15%.MOIC方法会造成一定程度上的旁瓣假象(图 3, 4)，如Q=60情况下补偿2 km时，旁瓣假象的振幅相当于主瓣的6%.当然，一般情况下，这样的旁瓣误差在RTM图像中对图像形状的影响也是有限的(图 6, 8).相位方面，两种方法的相速度拟合误差起决定性作用.在较高频范围(如15~70 Hz)，2级MORC优化的相速度误差大概为0.5%，3级优化中则缩小到0.3%(Zhou et al., 2018)，一般情况对波形补偿图像和RTM图像也都只产生非常小的影响.在较低频范围(如5~10 Hz)，2级、3级MORC优化的相速度误差快速上升到约10%和5%.这种程度的相速度误差已经对相位补偿造成了非常大的影响.如Q=60情况下补偿2 km时，就已经恢复不了波形.然而，MOIC方法不改变相速度，在全频率范围都不存在相速度拟合误差.
在勘探常用的频率范围，综合振幅和相位的补偿效果，MORC方法比MOIC方法更准确.在较低频率范围内更适合采用MOIC方法以更加准确地恢复子波形状，避免对解释造成干扰.在某些特殊情况下，如Q值很低且衰减区域非常大，两种方法都会导致误差累积.此时我们更推荐MOIC方法，因为在地球物理问题中往往相位对成像更为关键，保证相位准确而牺牲振幅在某种程度上更加可以接受.另外，在Q-FWI(即同时反演Q值和速度结构)中，我们也可以利用MOIC方法做低频范围的衰减补偿，以在相位不变的情况下提高振幅信息.在反演问题中，这有利于减轻病态问题，使收敛速度加快(Sun et al., 2015a).
在之前的MORC方法中，需要采用一个低通滤波器来防止高频噪声被不断放大.因为Q补偿的过程本身是一个振幅指数放大的过程，而常数Q又意味着高频分量放大的速率要远快于低频分量.如果不加限制，高频噪声将迅速增大以至于超过信号振幅，造成严重的不稳定现象.在MOIC方法中，一般情况下也需要采用低通滤波器.然而某些级的优化理论上是不需要低通滤波器的.单纯考虑Q补偿的影响，2级优化和4级优化理论上就不需要滤波器，因为它们在高频段1/Q值迅速增大，由负变正(图 2a, 2c).这样就使得高频分量变成了衰减.因此高频噪声不会被放大，可以稳定地传播.但3级优化由于高频段1/Q值迅速下降，趋于负无穷，导致补偿速率快速增大，需要通过很强的低通滤波器滤波才能使算法稳定.故我们一般不用3级优化.在实际计算中2级优化确实不需要滤波器，但4级优化仍然需要低通滤波器.这是由于4级优化中修正项出现了5阶和7阶导数.高阶空间导数数值非常大，并且对噪声等扰动十分敏感.为了保证求导的过程的稳定性，在5阶和7阶导数上加滤波器是必要的.我们选用频率采样法和Kaiser窗来设计低通滤波器，可以实现低阶数的高效线性相位滤波.
由多项式优化理论，通过提高优化级数，MORC方法和MOIC方法都可以更精确地逼近理想情况.然而优化级数的提高带来高阶导数.比如3级MORC方法用到4阶导数，而4级MOIC方法需要使用7阶导数.波动方程传播过程中，高阶空间导数的数值非常大，而其系数又非常小，很容易受噪声或数值误差影响引起不稳定.为了数值计算上的方便，我们需要尽量避免高阶导数的出现.数值测试也表明，采用更高阶的空间导数，就需要截止频率更低一些的滤波器来保证稳定性.而滤波器总会增加计算量以及损失一部分有用的高频信息.选择较低阶空间导数也有利于选用更高截止频率的滤波器，以保留更多高频信息.
4 结论我们发展了基于GSLS模型复模量虚部多级优化的Q-RTM方法.不同于之前的复模量实部优化，我们不采用负τ方案，而是保持复模量实部不变并对复模量的虚部进行多项式多级优化，同时引入负一次项来改善低频区域的拟合效果.最终得到可用的虚部优化Q-RTM方法.它适用于更宽的频率范围：5~70 Hz.RTM数值测试表明它可以准确补偿振幅衰减和相位扰动，从而提高了成像质量.作为基于GSLS模型的多项式优化Q-RTM方法，它能够利用有限差分实现，从而很好地保持了算子的局域性，在并行计算和GPU加速方面具有高度灵活性.此外，GSLS模型还可以通过调整标准线性体的数量和参数来模拟任意随频率变化的Q值，因此有可能将来应用于被动源成像.
致谢在本文数值测试过程中我们利用了一部分GPU加速算法，研究中受到了北京大学地球与空间科学学院硕士生崔丛越的启发和帮助，匿名审稿人提供了建设性的修改意见，极大地提高了论文质量.本研究得到了中国地震科学台阵探测——华北地区中部项目(DQJB16A0307)和地震行业专项(201408013)的资助，也得到了美国地球物理先进技术公司的支持.在此一并表示感谢!
Abdelkhalek R, Calandra H, Coulaud O, et al. 2012. Fast seismic modeling and reverse time migration on a GPU cluster. //International Conference on High Performance Computing & Simulation. Leipzig, Germany: IEEE, 36-43.
|
Abubakar A, Hu W Y, Habashy T M, et al. 2009. Application of the finite-difference contrast-source inversion algorithm to seismic full-waveform data. Geophysics, 74(6): WCC47-WCC58. DOI:10.1190/1.3250203 |
Abubakar A, Pan G, Li M, et al. 2011. Three-dimensional seismic full-waveform inversion using the finite-difference contrast source inversion method. Geophysical Prospecting, 59(5): 874-888. |
Ayala O, Wang L P. 2013. Parallel implementation and scalability analysis of 3D fast Fourier transform using 2D domain decomposition. Parallel Computing, 39(1): 58-77. DOI:10.1016/j.parco.2012.12.002 |
Bai M, Chen X H, Wu J, et al. 2016. Multiple-component Gaussian beam reverse-time migration based on attenuation compensation. Chinese Journal of Geophysics, 59(9): 3379-3393. DOI:10.6038/cjg20160921 |
Baysal E, Kosloff D D, Sherwood J W C. 1983. Reverse time migration. Geophysics, 48(11): 1514-1524. DOI:10.1190/1.1441434 |
Bickel S H, Natarajan R R. 1985. Plane-wave Q deconvolution. Geophysics, 50(9): 1426-1539. DOI:10.1190/1.1442011 |
Blanch J O, Robertsson J O A, Symes W W. 1995. Modeling of a constant Q:Methodology and algorithm for an efficient and optimally inexpensive viscoelastic technique. Geophysics, 60(1): 176-184. DOI:10.1190/1.1443744 |
Carcione J M, Kosloff D, Kosloff R. 1988. Wave propagation simulation in a linear viscoacoustic medium. Geophysical Journal International, 93(2): 393-401. DOI:10.1111/j.1365-246X.1988.tb02010.x |
Causse E, Ursin B. 2000. Viscoacoustic reverse-time migration. Journal of Seismic Exploration, 9(1): 165-184. |
Cavalca M, Moore I, Zhang L, et al. 2011. Ray-based tomography for Q estimation and Q compensation in complex media. //71st Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 3989-3993.
|
Claerbout J F. 1971. Toward a united theory of reflector mapping. Geophysics, 36(3): 467-481. DOI:10.1190/1.1440185 |
Clark R A, Benson P M, Carter A J, et al. 2009. Anisotropic P-wave attenuation measured from a multi-azimuth surface seismic reflection survey. Geophysical Prospecting, 57(5): 835-845. DOI:10.1111/gpr.2009.57.issue-5 |
Cui J J, He J S. 2001. Viscoelastic wave equation forward modeling and migration. Journal of Central South University of Technology (Natural Science), 32(5): 441-444. DOI:10.3969/j.issn.1672-7207.2001.05.001 |
Dai N X, West G F. 1994. Inverse Q migration. //64th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 1418-1421.
|
Deng F, McMechanG A. 2007. True-amplitude prestack depth migration. Geophysics, 72(3): S155-S166. DOI:10.1190/1.2714334 |
Deng F, McMechanG A. 2008. Viscoelastic true-amplitude prestack reverse-time depth migration. Geophysics, 73(4): S143-S155. DOI:10.1190/1.2938083 |
Deng W Z, Li Z C, Wang Y G, et al. 2015. The least-squares reverse time migration for visco-acoustic medium based on a stable reverse-time propagator. Geophysical and Geochemical Exploration, 39(4): 791-796. DOI:10.11720/wtyht.2015.4.22 |
Dutta G, Schuster G T. 2014. Attenuation compensation for least-squares reverse time migration using the viscoacoustic-wave equation. Geophysics, 79(6): S251-S262. DOI:10.1190/geo2013-0414.1 |
Dutta G, Schuster G T. 2016. Wave-equation Q tomography. Geophysics, 81(6): R471-R484. DOI:10.1190/geo2016-0081.1 |
Etgen J, Gray S H, Zhang Y. 2009. An overview of depth imaging in exploration geophysics. Geophysics, 74(6): WCA5-WCA17. DOI:10.1190/1.3223188 |
Guo P, McMechan G A, Guan H M. 2016. Comparison of two viscoacoustic propagators for Q-compensated reverse time migration. Geophysics, 81(5): S281-S297. DOI:10.1190/geo2015-0557.1 |
Hargreaves N D, Calvert A J. 1991. Inverse Q filtering by Fourier transform. Geophysics, 56(4): 519-527. DOI:10.1190/1.1443067 |
Hu W, Cummer S A. 2004. The nearly perfectly matched layer is a perfectly matched layer. IEEE Antennas and Wireless Propagation Letters, 3: 137-140. DOI:10.1109/LAWP.2004.831077 |
Hu W Y, Abubakar A, Habashy T M. 2007. Application of the nearly perfectly matched layer in acoustic wave modeling. Geophysics, 72(5): SM169-SM175. DOI:10.1190/1.2738553 |
Hu W Y, Liu J, Bear L, et al. 2011. A robust and accurate seismic attenuation tomography algorithm. //71st Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts.
|
Hu W Y, Zhou T, Ning J Y. 2016. An efficient Q-RTM algorithm based on local differentiation operators. //86th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 4183-4187.
|
Kjartansson E. 1979. Constant Q-wave propagation and attenuation. Journal of Geophysical Research:Solid Earth, 84(B9): 4737-4748. DOI:10.1029/JB084iB09p04737 |
Kolsky H. 1956. LXXI. The propagation of stress pulses in viscoelastic solids. Philosophical Magazine, 1(8): 693-710. DOI:10.1080/14786435608238144 |
Li J L, Li Z C, Guan L P, et al. 2015. The method of seismic attenuation and energy compensation. Geophysical and Geochemical Exploration, 39(3): 456-465. DOI:10.11720/wtyht.2015.3.04 |
Liu X W, Nian J B, Liu H. 2006. Generalized S-transform based compensation for stratigraphic absorption of seismic attenuation. Geophysical Prospecting for Petroleum, 45(1): 9-14. DOI:10.3969/j.issn.1000-1441.2006.01.002 |
Lomnitz C. 1957. Linear dissipation in solids. Journal of Applied Physics, 28(2): 201-205. DOI:10.1063/1.1722707 |
Müller T M, Gurevich B, Lebedev M. 2010. Seismic wave attenuation and dispersion resulting from wave-induced flow in porous rocks-a review. Geophysics, 75(5): 75A147-75A164. DOI:10.1190/1.3463417 |
Robertsson J O, Blanch J O, Symes W W. 1994. Viscoelastic finite-difference modeling. Geophysics, 59(9): 1444-1456. DOI:10.1190/1.1443701 |
Shen Y, Zhu T Y. 2015. Image-based Q tomography using reverse time Q migration. //85th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 3694-3698.
|
Sun J Z, Fomel S, Zhu T Y. 2015a. Preconditioning least-squares RTM in viscoacoustic media by Q-compensated RTM. //85th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 3959-3965.
|
Sun J Z, Zhu T Y, Fomel S. 2015b. Viscoacoustic modeling and imaging using low-rank approximation. Geophysics, 80(5): A103-A108. DOI:10.1190/geo2015-0083.1 |
Traynin P, Liu J, Reilly J M. 2008. Amplitude and bandwidth recovery beneath gas zones using Kirchhoff prestack depth Q-migration. //78th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 2412-2416.
|
Wu Y, Fu L Y, Chen G X. 2017. Forward modeling and reverse time migration of viscoacoustic media using decoupled fractional Laplacians. Chinese Journal of Geophysics, 60(4): 1527-1537. DOI:10.6038/cjg20170425 |
Xie Y, Sun J, Zhang Y, et al. 2015. Compensating for visco-acoustic effects in TTI reverse time migration. //85th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 3996-4001.
|
Xin K, Hung B. 2010. 3D Q tomography for compensating frequency dependent attenuation and dispersion. //72nd EAGE Conference and Exhibition Incorporating SPE EUROPEC 2010, Expanded Abstracts.
|
Xu W C, Yang G Q, Li Z C, et al. 2016. Pseudo acoustic equation for TI medium attenuation based on the GSLS model. Chinese Journal of Geophysics, 59(6): 2232-2244. DOI:10.6038/cjg20160626 |
Yu Y, Lu R S, Deal M M. 2002. Compensation for the effects of shallow gas attenuation with viscoacoustic wave-equation migration. //72nd Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts:, 2062-2065.
|
Zhang G L, Wang X M, Ge Z H, et al. 2014. Interval Q inversion based on zero-offset VSP data and applications. Applied Geophysics, 11(2): 235-244, 255. DOI:10.1007/s11770-014-0416-6 |
Zhang J F, Wapenaar C P A. 2002. Wavefield extrapolation and prestack depth migration in anelastic inhomogeneous media. Geophysical Prospecting, 50(6): 629-643. DOI:10.1046/j.1365-2478.2002.00342.x |
Zhang J F, Wu J Z, Li X Y. 2013. Compensation for absorption and dispersion in prestack migration:An effective Q approach. Geophysics, 78(1): S1-S14. DOI:10.1190/geo2012-0128.1 |
Zhang Y, Zhang P, Zhang H Z. 2010. Compensating for visco-acoustic effects in reverse-time migration. //80th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 3160-3164.
|
Zhao Y, Liu Y, Ren Z M. 2014. Viscoacoustic prestack reverse time migration based on the optimal time-space domain high-order finite-difference method. Applied Geophysics, 11(1): 50-62. DOI:10.1007/s11770-014-0414-8 |
Zhou J, Birdus S, Hung B, et al. 2011. Compensating attenuation due to shallow gas through Q tomography and Q-PSDM, a case study in Brazil. //81st Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 3332-3336.
|
Zhou T, Hu W Y, Ning J Y. 2018. An efficient local operator-based Q-compensated reverse time migration algorithm with multistage optimization. Geophysics, 83(3): S249-S259. DOI:10.1190/geo2017-0026.1 |
Zhu T Y, Harris J M. 2014a. Modeling acoustic wave propagation in heterogeneous attenuating media using decoupled fractional Laplacians. Geophysics, 79(3): T105-T116. DOI:10.1190/geo2013-0245.1 |
Zhu T Y, Harris J M, Biondi B. 2014b. Q-compensated reverse-time migration. Geophysics, 79(3): S77-S87. DOI:10.1190/geo2013-0344.1 |
Zhu T Y, Sun J Z. 2017. Viscoelastic reverse time migration with attenuation compensation. Geophysics, 82(2): S61-S73. DOI:10.1190/geo2016-0239.1 |
白敏, 陈小宏, 吴娟, 等. 2016. 基于吸收衰减补偿的多分量高斯束逆时偏移. 地球物理学报, 59(9): 3379-3393. DOI:10.6038/cjg20160921 |
崔建军, 何继善. 2001. 粘弹性波动方程正演和偏移. 中南大学学报(自然科学版), 32(5): 441-444. DOI:10.3969/j.issn.1672-7207.2001.05.001 |
邓文志, 李振春, 王延光, 等. 2015. 基于稳定逆时传播算子的黏声介质最小二乘逆时偏移. 物探与化探, 39(4): 791-796. DOI:10.11720/wtyht.2015.4.22 |
李金丽, 李振春, 管路平, 等. 2015. 地震波衰减及补偿方法. 物探与化探, 39(3): 456-465. DOI:10.11720/wtyht.2015.3.04 |
刘喜武, 年静波, 刘洪. 2006. 基于广义S变换的吸收衰减补偿方法. 石油物探, 45(1): 9-14. DOI:10.3969/j.issn.1000-1441.2006.01.002 |
吴玉, 符力耘, 陈高祥. 2017. 基于分数阶拉普拉斯算子解耦的黏声介质地震正演模拟与逆时偏移. 地球物理学报, 60(4): 1527-1537. DOI:10.6038/cjg20170425 |
徐文才, 杨国权, 李振春, 等. 2016. 基于GSLS模型TI介质衰减拟声波方程. 地球物理学报, 59(6): 2232-2244. DOI:10.6038/cjg20160626 |
张固澜, 王熙明, 贺振华, 等. 2014. 零偏移距VSP资料层Q反演及其应用研究. 应用地球物理, (2): 235-244, 255. DOI:10.1007/s11770-014-0416-6 |
赵岩, 刘洋, 任志明. 2014. 基于优化时空域高阶有限差分方法的粘滞声波叠前逆时偏移. 应用地球物理, 11(1): 50-62. DOI:10.1007/s11770-014-0414-8 |