2014, Vol. 57
叠前地震成像是通过地震波所携带的信息认识地下构造的重要方法.随着油气勘探的深入及勘探目标复杂程度的增加，构造成像已不能满足勘探的需求，对岩性成像的需求越来越迫切.即在保证成像位置准确的情况下，还需要达到真振幅成像(李振春等，2008).大部分真振幅成像算法是建立在完全弹性介质条件的假设下，然而实际地下介质，特别是在近地表及油气储层，多为一种黏弹性介质.地震波在黏弹性介质中的传播主要表现为速度频散与振幅衰减.不考虑黏滞性的叠前成像算法不仅会使成像位置发生偏离，而且还会引起成像振幅的欠估计，严重影响甚至误导随后的地震数据处理、解释等工作.
对由黏滞性引起的地震波吸收衰减的补偿，常规的方法是进行反Q滤波以增强地震波特别是高频成分的能量，达到提高分辨率的效果(Wang，2002，2006；李振春等，2007).但常规的反Q滤波方法大部分基于层状介质假设，虽能够对能量进行一定补偿，但在复杂介质条件下，该方法并不符合地震波的传播规律，不是一种精确的补偿方法(Zhang et al.，2010).另一种方法是在成像的过程中对地震波进行能量补偿及相位校正，即反Q偏移方法.目前已有的方法主要有三类，即基于射线理论(Ribodetti et al., 1998)、基于单程波方程与基于双程波方程的反Q偏移方法.基于射线理论的成像方法由于基于高频近似，较难处理多波至问题，在复杂介质模型中的应用受到限制.基于波动理论的偏移方法在复杂介质条件下的成像精度更高.常规反Q偏移的思路较为简单，主要是波场外推中的不稳定问题限制了方法的应用与发展.Mittet et al.(1995)首先提出了黏声介质中单程波偏移的振幅补偿和频散校正策略，随后Zhang et al.(2002)、Wang et al.(2004)、Zhou et al.(2012)借鉴Mittet的思想应用不同的外推算子对相应的理论进行了发展.在处理不稳定问题时，采用增益截频或对传播角度进行限制的方法以使得算法稳定.Wang,2008借鉴反Q滤波中的思想，提出了基于稳定因子法的反Q滤波偏移方法并着重讨论了稳定性问题.张立彬等(2010)将两种方法进行结合得到了一种新的稳定性控制方法.稳定性的问题在双程波逆时偏移中同样存在，且由于时间采样点数较多，不稳定问题更为严重.Causse et al.(1999，2000)在黏声介质全波形反演中通过分离频散项与振幅项推导了稳定的振幅增加介质，并将其应用到黏声介质的逆时偏移中.该方法虽较好地解决了稳定性问题，但较难推广到其他吸收介质模型中，并且计算量显著增加.Deng et al.(2007,2008)分别实现了黏声介质及黏弹性介质中的真振幅逆时偏移，但文中对稳定性问题并未进行过多讨论.Zhang Y et al.(2010)基于频散关系推导了新的黏声介质中的波动方程，并指出通过加入正则化项可使得外推稳定；Suh et al.(2012)将Zhang Y et al.(2010)的方法推广到黏声各向异性介质，并指出将高切滤波器作用于耗散项使得外推稳定.综上所述无论是单程波方法还是双程波方法，都存在稳定性的问题.虽然通过一些方法技术可使得算法稳定，但这在一定程度上也损害了波场的真实信息、影响成像精度.
最小平方偏移的理论基础最早由Tarantola(1984)提出，但由于计算能力的限制，发展缓慢.之前的最小平方成像算法主要基于Kirchhoff算子与单程波算子(Nemeth et al., 1999; 王彦飞等，2009；Tang et al, 2009；刘玉金等，2013)，近几年发展了基于双程波动方程的最小平方逆时偏移方法(LSRTM).Dai et al.(2010)通过相位编码技术有效地降低了LSRTM的计算成本并将其应用到3D数据；(Yao et al.,2012a，2012b)提出了线性及非线性的LSRTM；Dong et al.(2012)通过模型试算及野外实际资料处理结果验证了LSRTM具有振幅保持、高成像分辨率等特性.在黏声介质中的最小平方偏移方面研究较少，Zhang C et al.(2010)基于最小平方偏移方法推导了时变反褶积算子，并将该算子作用于成像剖面以得到对吸收衰减的校正；Mulder et al.(2009)与Hak et al.(2010，2011)在频率-空间域理论分析了黏声介质成像的多解性问题，并将其 应用到黏声介质的线性及非线性反演.虽然在频率-空间域处理黏滞性问题较为方便，但数值实现需要求解大型线性方程组，在处理大尺度3D模型时受到限制.
本文将地下介质的黏滞性引入LSRTM，在避免了不稳定性的同时对吸收衰减进行补偿，得到地下介质构造的高精度真振幅成像.本文所述方法在时间域求解问题，相对于Mulder的方法较易扩展到三维大尺度模型.虽然基于广义标准线性固体(GSLS)模型进行相应算法的推导，但可较易将其扩展到其他黏滞性模型.首先给出了基于GSLS模型的黏声介质波动方程并推导了其线性化形式，然后借助伴随状态法推导了黏声介质LSRTM的迭代算法，并引入动态相位编码技术进一步提高计算效率，最后通过模型试算验证了方法的正确性及有效性.
基于GSLS模型的黏声介质波动方程(Robertsson，et al.，1994)为
其中,MR为松弛模量；τεi与τσi为松弛时间；H(t)为单位阶跃函数；I为标准线性固体的个数；
为梯度算子；· 为散度算子,*为时间上的卷积算子.由于地下介质密度与地层品质因子是渐变的，反射波可认为主要是由速度突变产生的，因此本文只对慢度扰动进行成像.将松弛模量用速度与密度表示并进行整理可得
其中,
.波场p及震源项f为空间位置与时间的函数，介质参数v、ρ、τεi、τσi均为空间位置的函数.为了公式表述的简洁，在公式推导中会忽略一部分函数的自变量.
由场的叠加原理，可以将波场分解为在背景波场和扰动波场的叠加，即
其中,p0为背景波场；ps为扰动波场.背景波场是由震源激发并在背景介质中传播的波场，满足方程
定义慢度扰动量为
将式(3)—(5)代入式(2)，并应用Born近似可得到扰动波场的控制方程
由(6)式可知，扰动波场是背景波场与慢度扰动的相互作用作为二次震源产生的波场，与慢度扰动呈线性关系.在实际地震资料处理中，通过偏移速度分析、层析等速度分析方法(秦宁等，2013)获取的宏观速度模型可视为此处的背景介质；通过偏移成像等方法得到的成像结果即为此处的慢度扰动.由此可知，上述基于场的叠加及Born近似的波动方程线性化是合理的.
为了随后公式推导的简洁，可以将波动方程及其线性化写成算子的形式，即
其中,
.为方便波动方程的数值求解，引入辅助变量:
通过整理，原始方程可化简为
通过引入辅助变量ri，将原始波动方程中的时间褶积转化为相应的辅助微分方程，减少了直接计算褶积对内存的消耗.此外通过引入辅助变量ux、uz将原有的空间二阶求导转化为一阶求导，可在空间交错网格上数值求解波动方程，提高了正演的精 度，也方便了PML边界条件的引入.本文采用未分裂形式的CFS-PML(Complex-Frequency-Shifted Perfectly Matched Layer, 复频移完全匹配层)边界条件对边界反射进行吸收，通过辅助微分方程进行实现.
基于反演的成像方法寻求最优的地下介质模型以使得正演波场与观测波场残差的模最小，是一个最小范数问题.为保持公式推导的简洁，此处以单炮成像为例进行相应理论的推导，对于多炮成像为单炮情况下的叠加.为推导伴随算子，定义数据空间的标量乘(Tanrantola, 2005)为
模型空间的标量乘为
其中,V( x )为空间范围；t0、t1为初始及终止时间.
定义LSRTM的目标函数为
其中,第一项为数据匹配项，第二项为模型约束项.Wd为数据加权函数，主要控制数据在目标函数中的权重；Wm为模型加权函数，可通过迭代再加权方法实现对模型的L1模稀疏约束.
对于梯度类的局部寻优算法，首先需求取梯 度，即目标函数相对于模型参数的导数.通过推导可得
其中,
为偏导数波场，反应扰动波场相对于模型参数的敏感性.偏导数波场可以通过对正演方程(7b)两边相对于模型参数求导得到
由(14)式可以看到，地下 x′ 处的偏导数波场为将该点视为绕射点得到的散射波场，具有明确的物理意义.结合式(13)可进一步对梯度进行物理解释：梯度存在于那些偏导数波场与数据剩余量最佳匹配的点上.通过在这些点上对模型进行更新，可进一步减少数据剩余量，以达到最终数据匹配的结果.
从计算的角度看，对于地下任一点的偏导数波场都需要一次正演过程，若利用式(13)计算梯度，计算量巨大，目前计算设备无法承受.借助伴随状态法，可以有效的提高计算效率.
将偏导数波场的显式形式代入式(13)，并利用正演算子的伴随算子可得
关于伴随算子的详细推导见附录A.定义伴随波场为
则梯度向量可表示为
由上式可知，梯度向量为背景波场的时间二阶导数与伴随波场的零延迟互相关，其中伴随波场为在伴随算子及终止条件控制下的波场剩余量的反向传播.可以证明，与常规的黏声介质逆时偏移方法不同，当满足CFL稳定性条件时，伴随算子在反向外推的过程中是稳定的，不存在不稳定性的问题.从计算的角度，借助伴随状态法将梯度向量的计算简化为一次背景波场的正向传播与一次伴随波场的反向传播，相对于式(13)极大的提高了计算效率.
利用最速下降法，可对模型进行迭代更新，直至收敛到终止条件.模型更新过程为
其中, P 为预处理算子.最优的预处理算子应为Hessian矩阵的逆，但计算量巨大，直接求取并不现实.本文用背景波场的能量来近似Hessian矩阵的对角元素，在减少计算量的同时加速了收敛速度.
在黏声介质中，由于吸收衰减的影响，深层反射波能量很弱，使得这种弱信号在迭代过程中权重很小，在一定程度上降低了收敛效率.可通过修改数据加权矩阵Wd增加相应数据的权重，在一定程度起到增加收敛速率的作用.此外，由于数据并不能完全 描述模型的各个分量，在迭代过程中需要加入模型的先验约束，本文采用对模型的L1模稀疏约束，通过迭代再加权方法实现.
对于LSRTM方法，每次迭代至少需要3次波场外推.由于线性化波动方程的正演相当于两次常规正演模拟，因此对于N次迭代，其计算成本大致是常规逆时偏移方法的2.5N倍.庞大的计算量限制了方法的应用.
不管LSRTM还是常规的逆时偏移方法，其计算量都与所用的炮数成线性关系.与面炮偏移思想类似(杨敬磊等，2008)， 通过相位编码技术将多个炮集编码成一个超道集，可有效地提高计算效率；利用动态相位编码技术可在迭代过程中压制由于相位编码引 起的串扰噪声(Schuster et al.，2011).编码之后LSRTM 的计算成本可与常规逆时偏移在同一个数量级上.
超道集的构建可用下式表示:
其中,s为炮号；e(s)为编码函数.本文所用编码函数为随机的1，-1，编码函数在每次迭代中随机生成.由于正演波场与震源的线性关系，一个超道集只需一次正演计算就可得到.
将上式带入梯度公式并做相应整理可得：
由式(20)可知在编码后的超道集上求取的梯度不仅有正确的梯度(式中第一项)，还引入了不同炮之间的串扰噪声(式中最后一项).
由于在迭代过程中动态的更新编码函数，并且不同迭代之间串扰噪声是随机的，通过多次迭代可进一步增强真实的成像结果，压制随机的串扰噪声.此外可通过增加每次迭代过程中所用超道集的个数，降低串扰噪声的干扰，进一步提高信噪比.极限情况，当所用超道集的个数与原始炮数相同时，退化为原始的最小平方偏移算法.
为验证方法的正确性，分别用简单平层模型、二维SEG/EAGE盐丘模型对算法进行测试.其中通过平层模型主要分析黏滞性的引入对叠前数据及成像结果的影响；通过盐丘模型主要测试方法对振幅的恢复、偏移噪声的压制、盐下弱照明区域成像的能力，同时对比分析相位编码的引入对效率的提高.
所用平层模型为如图1a所示的三层介质模型， 横纵向分别为150个采样点，空间采样率为10 m.从上到下三层介质的速度分别为1500，2000，2500 m/s； 品质因子分别为30，60，100.在地表以30 m间隔均匀分布50个炮点，主频20 Hz的雷克子波作为震源子波进行正演模拟，150个在地表以10 m间隔均匀分布的检波点进行接收.
 | 图1 平层模型 (a)速度模型；(d)慢度扰动. Fig.1 Three layers model (a) Velocity model; (b) Slowness perturbation. |
将原始模型平滑之后的速度场作为背景速度场，进而可由式(5)得到慢度扰动，该慢度扰动可视为理想的成像结果，如图1b所示.基于线性化的波动方程可正演得到没有多次波的正演记录，并将此记录作为成像的输入道集.为对比分析黏滞性对波场的影响，分别正演得到黏声及声波数据.为定量对比反射波场的振幅、波形及走时差异，将黏声数据炮点左侧数据与声波数据炮点右侧数据进行拼接后成图，如图2a所示，图中曲线为对应位置处的波形.由图中可以看到，相对于声波数据，黏声数据中反射波由于吸收导致振幅严重衰减，同时由于速度频散导致走时差异及波形的变化.图2b为两种数据第一个反射轴的振幅谱对比，可以看到由于吸收衰减的影响导致反射波总体能量降低且主频往低频端移动，高频吸收严重.若在成像过程中不对这些效应进行消除，则会引起成像位置错误、成像振幅偏低、成像分辨率降低等问题.
 | 图2 黏声数据与声波数据对比 (a)炮集对比，左侧为黏声介质正演记录，右侧为声波介质正演记录；(b)第一个反射轴振幅谱对比，第一道为黏声数据振幅谱，第二道为声波数据振幅谱. Fig.2 Comparison of visco-acoustic and acoustic data (a) Comparison of shot gathers and waveform, left is visco-acoustic data, right is acoustic data; (b) Comparison of amplitude spectrum, trace 1 is visco-acoustic data, trace 2 is acoustic data. |
利用本文所述方法，对黏声介质下正演数据分 别进行黏声介质及常规声波LSRTM，结果如图3所示.由于在成像过程中利用近似的Hessian矩阵对角元素对梯度进行预处理，模型中心位置在第一次迭代的成像结果已近似保振幅，如图3a及图3d中第三条曲线所示.在模型两侧，由于本身的照明不均衡及预处理方法不是精确的Hessian矩阵，振幅明显欠估计.迭代80次之后的成像结果如图3b所示.对比图3a与图3b可知，随着迭代次数的增加，模型的高波数成分得到恢复，低波数噪声得到有效压制，同时有效地对由于照明不均衡引起的振幅横向变化进行了校正.对比图3d中第3、4条曲线可以看到，多次迭代之后的成像结果在振幅及波形上更加逼近真实的慢度扰动.
 | 图3 成像结果对比 (a)黏声介质LSRTM第1次迭代结果；(b)黏声介质LSRTM第80次迭代结果；(c)对黏声数据应用常规声波LSRTM成像结果；(d)水平 位置750 m处成像结果与真实慢度扰动模型对比，其中第1条曲线为真实慢度扰动，第2、3、4条曲线中实线分别对应图c、a、b在750 m 处成像曲线，图中虚线为真实慢度扰动. Fig.3 Comparison of image results (a) Image result at the 1st iteration of visco-acoustic LSRTM; (b) Image result at the 80th iteration of visco-acoustic LSRTM; (c) Image result of acoustic LSRTM with visco-acoustic data; (d) Comparison of image results at x=750 m. The 1st line is the true slowness perturbation, while the 2nd, 3rd ,4th line is image result from figure c, a, b respectively. The dot line is also the true slowness perturbation. |
常规声波LSRTM的成像结果如图3c所示.对比图3b可知，基于声波介质的成像方法不能有效地对由黏滞性引起的速度频散进行校正，导致了成像界面位置错误，并引入部分成像假象.由于振幅衰减未能得到补偿，成像振幅远小于真实慢度扰动，且分辨率较低，如图3d所示.
原始二维盐丘模型主要用来测试存在高速异常体时成像算法的优劣，特别是盐下弱照明区域的成像，速度模型如图4a所示.本文在原有速度模型的基础上利用李氏经验公式(李庆忠，1993)计算品质因子来构建黏声介质模型，以此来测试本文所提算法在复杂介质模型下的应用效果.首先通过对比初次迭代与最终反演结果说明该方法在成像噪声压制、振幅恢复等方面的优势，并通过成像结果与理论 成像结果的对比验证该方法真振幅成像的能力；然后通过与常规声波介质LSRTM结果的对比，进一步验证在成像过程中考虑黏滞性对成像分辨率提高、振幅恢复等方面的必要性；最后对比分析了相位编码技术的引入对计算效率的提高及对成像结果的影响.
 | 图4 盐丘模型 (a)速度模型；(b)品质因子模型；(c)背景速度模型；(d)慢度扰动模型. Fig.4 Salt model (a) Velocity model; (b) Quality factor model; (c) Background velocity model; (d) Slowness perturbation. |
所用模型横向采样点数为780，纵向采样点数为209，空间采样率为20 m.表层最小速度为1500 m/s，盐丘最大速度为4450 m/s，最小品质因子为34，最大品质因子为373.该模型在表层吸收最为严重，盐丘对应位置吸收最弱，模型其余部分品质因子在50到100之间.图4c所示为背景速度模型，图4d为对应的慢度扰动，可视为理想的成像结果.在地表以160 m间隔均匀分布97个炮点进行正演模拟，震源子波是主频10 Hz的雷克子波，每炮780个检波点进行接收，检波点在地表每个网格点上都有定义.随后的数值测试利用该正演数据进行.
黏声介质LSRTM第一次迭代结果如图5a所示，与常规逆时偏移方法类似，成像结果被严重的低频噪声覆盖，特别是盐丘顶部最为明显，如图5a中椭圆虚线框所示.虽然利用背景波场能量的倒数进行预处理使得中深层能量得到了相对均衡，但吸收衰减及成像噪声严重影响成像结果，该成像结果与真振幅成像差距很大.对图5a中结果进行高通滤波可在一定程度上消除原有的低频噪声，但该类成像后处理方法并不是保幅的处理，在一定程度上破坏了原有的振幅相对强弱关系，如在震源附近引入了更强的噪声(如图5b中矩形虚线框中所示)、盐丘顶部能量也受到了一定的损失(如图5b中椭圆虚线框所示)等.此外，由于吸收衰减的影响，成像振幅及成像分辨率均偏低，特别是在深层尤为明显.迭代50次之后的结果如图5c所示.对比图5a与图5c可知，通过迭代反演，不仅有效的压制了低频成像噪声，同时还对由于波场的吸收衰减及弱照明等因素 引起的成像振幅欠估计进行了补偿，成像分辨率也得到了极大的提高.在图6中对成像结果与真实慢度扰动进行了精确的曲线对比，从中也可看到成像曲线与真实慢度扰动基本吻合，达到了真振幅成像的目的.
 | 图5 盐丘模型成像结果对比分析 (a)黏声介质LSRTM第一次迭代成像结果；(b)图a进行高通滤波之后的结果；(c)黏声介质LSRTM 50次迭代之后成像结果；(d)常规声波LSRTM 50次迭代之后成像结果； (e)深度方向波数谱对比， 其中标号1,2,3分别为常规声波LSRTM、黏声介质LSRTM、理论成像结果的波数谱；(f)图5d中A、B位置处成像结果曲线对比，实线为常规声波LSRTM，虚线为黏声介质LSRTM. Fig.5 Comparison of image results of salt model (a) Image result at 1st iteration; (b) applying high-pass filtering to Figure 5a; (c) Image result at 50th iteration of visco-acoustic LSRTM; (d) Image result at 50th iteration of acoustic LSRTM; (e) Comparison of wavenumber spectrums, where the trace numbered with 1, 2, 3 is the wavenumber spectrum of the image result of acoustic LSRTM, visco-acoustic LSRTM, ideal image result respectively; (f) Comparison of image result in Figure 5c (dash line) and 5d (real line) at location A and B. |
 | 图6 成像结果与真实慢度扰动对比图 其中背景图片为真实慢度扰动，红色虚线为在该CDP点的慢度扰动曲线，蓝色实线为成像结果曲线. Fig.6 Comparison of image result and true slowness perturbation The background picture is the image of true slowness perturbation; the red dot line is exact value of slowness perturbation at this CDP; the blue line is the image result at this CDP. |
由图5c及图6可知，除图5c中矩形虚线框所示区域外，其余区域成像结果均达到了较高精度.通过照明分析等研究表明，该区域成像效果不理想的主要原因是照明极弱.对于盐丘上表面方框所示区域，由于此处倾角较陡，在此处的反射波基本无法在检波点处接收到；对于盐丘以下方框所示区域，由于盐丘的阻挡由震源到达此处的地震波能量本身就很弱，而反射波在盐丘下侧表现为很强的折射效应，透射波能量同样很弱，这也导致了该处几乎无法成像.由此可见，基于迭代反演的成像方法可对弱照明区域进行一定的振幅补偿，但是当照明极弱近乎盲区时，该方法依然不能准确成像.若加入额外的人为约束条件，可望得到此处的像，但成像结果的准确性依赖约束条件的合适与否.
为了进一步说明本文所述方法在振幅补偿及分辨率提高等方面的优势，对前述正演数据利用常规声波LSRTM进行成像以作对比.将黏滞性参数的影响消除，本文所述的黏声介质LSRTM便退化为常规声波LSRTM算法.保持其余参数不变，50次迭代反演之后的结果如图5d所示.对比图5c与5d可以看到，忽略黏滞性的影响使得成像噪声增加、成像分辨率降低、成像振幅欠估计，相对于图5c该成像结果在信噪比、保真度等方面明显降低.
成像噪声的增加主要是由于地震波在黏滞性介质与完全弹性介质中传播规律的不耦合引起的，因为很难找到一种声波介质模型使得其正演记录与黏声介质下的正演记录完全匹配.这种不耦合对于浅层大偏移距的情况最为严重，因此浅层的成像噪声明显强于深层.为定量对比成像分辨率，分别求取常 规声波LSRTM、黏声介质LSRTM成像结果及真实慢度扰动在深度方向的波数谱并在横向上进行叠加得到平均波数谱，分别如图5e中1、2、3所示.虽然该波数谱仅反映深度方向上的分辨率，但在一定程度也能够对比成像分辨率的高低.通过对比可知，黏声介质LSRTM成像结果最接近真实慢度扰动，成像分辨率更高，有效的恢复了模型的高波数成分；常规声波LSRTM的波数谱相对来说有效宽度更窄，对应成像分辨率最低，未能有效的恢复模型的高波数成分.此外常规声波LSRTM波数谱总体能量较弱，对应空间域成像振幅的欠估计.为更详细的对比成像振幅的强弱，分别从图5c、图5d中A、B位置处抽取成像曲线进行对比，如图5f所示.通过成像曲线的对比可以看到，声波介质LSRTM成像振幅普遍欠估计且随着黏滞性增强振幅欠估计越严重，此外成像分辨率较低.
通过以上的对比分析，说明了在成像过程中考虑介质黏滞性必要性的同时，也验证了本文所述方法具有压制成像噪声、对吸收衰减进行补偿、提高成像分辨率等优势，能够实现对黏声介质的真振幅成像.该方法一个主要的缺点是计算量庞大，通过引入相位编码技术可有效减少计算量.本节后半部分主要分析相位编码的引入对成像结果的影响及计算效率的提高.
为测试相位编码技术的引入对成像结果的影响及对计算效率的提高，将炮集分别编码成1、4、8个超道集进行模型试算.图7为将炮集编码成一个超道集的成像结果.图7a为第10次迭代的成像结果，在其中虽能看到大体的构造，但是成像剖面被严重的串扰噪声所覆盖，信噪比很低.将图7a中结果与常规LSRTM第10次迭代结果做差，可进一步验证串扰噪声是由相位编码引入的，如图7b所示.随着迭代次数的增加，在100次迭代时，串扰噪声得到有效地压制，成像剖面信噪比得到极大的提高，如图7c所示.将图7c中所示成像结果与常规LSRTM第50次迭代结果进行求差得到两者差剖面，如图7d所示.绘图所用差剖面为真实差值的十倍，图7c与图7d成图时所用色标相同.由差剖面可以清楚的看到，将炮集编码成一个超道集的LSRTM方法100次迭代后的结果与常规LSRTM 50次迭代后的结果相当，但计算成本是原先算法的0.0206倍，大约为原先的1/50，极大的提高了计算效率.
 | 图7 基于相位编码LSRTM成像结果对比分析 (a)第10次迭代结果；(b)第10次迭代结果与常规LSRTM第10次迭代结果的差；(c)第100次迭代结果； (d)第100次迭代结果与常规LSRTM第50次迭代结果的差的10倍. Fig.7 Comparison and analysis of image result using LSRTM with phase-encoding (a) Image result at 10th iteration; (b) Difference between image result using phase-encoding and conventional LSRTM at 10th iteration; (c) Image result at 100th iteration using phase-encoding; (d) Ten times the difference between image result in figure c and the one using conventional LSRTM at 50th iteration. |
由于利用当前PC-cluster等高性能计算设备较易实现不同炮集之间的并行处理，可以通过增加超道集的数目来进一步提高收敛的速度，降低串扰噪声的影响.极限情况下，当超道集的数目与原始炮集数目相等时，基于相位编码的LSRTM退化为常规的LSRTM.图8为分别应用4、8个超道集在第10次迭代时的成像结果.对比图7a与图8a、b可以看到，随着超道集数目的增加，由相位编码引起的串扰噪声得到更好的压制.在实际处理中可根据现有的计算设备进行优化选取.
 | 图8 应用多个超道集的基于相位编码LSRTM第10次迭代成像结果 (a)基于4个超道集的成像结果；(b)基于8个超道集的成像结果. Fig.8 Image result at 10th iteration using multiple supergathers (a) The image result using 4 supergathers; (b) The image result using 8 supergathers. |
针对黏滞性吸收对成像振幅及分辨率的重要性，本文提出了一种时间域黏声介质条件下的LSRTM方法.该方法相对于常规反Q偏移方法主要有以下优势：(1)在避开不稳定性问题的同时对黏滞性吸收进行补偿；(2)能够实现对地下介质的真振幅成像；(3)能够有效地压制成像噪声并对弱照明区域进行照明补偿；(4)基于动态相位编码的LSRTM可有效地降低计算成本，将其计算量降至与常规逆时偏移方法同一数量级.
本文虽然在二维情况下对基于GSLS模型的黏声LSRTM方法进行了相关的理论推导，但可较方便的推广到三维及其他吸收介质模型.此外由于在时间域进行算法的推导，相对于频率域方法在处理大尺度模型时，该方法计算量更小且更适用于当前的PC-Cluster集群等高性能计算设备.
虽然本文所述方法在理论上是完善的，且通过模型试算验证了方法的正确性及有效性，但若将其应用到实际资料还有一系列的难题需要克服.这些困难主要包括：(1)保振幅的前期处理与评价工作，主要涉及去噪、震源能量一致性校正、面波等弹性效应去除等叠前数据处理及保幅评价方法研究；(2)高精度的速度模型估计与品质因子提取.对于第一点，这是所有真振幅成像及波形反演等方法所共有的问题，本文对其并不做过多评价.对于第二点主要是该方法对精确背景介质参数模型的依赖性问题，作者认为在实际资料的处理中主要可从以下几个方面考虑：
(1)建立合理的地下介质参数成像流程分别对不同波数成分进行成像，如通过层析成像、基于成像道集的偏移速度分析方法、Laplace域波形反演等技术手段可获取地下介质模型的低波数成分；通过频率域或时间域波形反演方法可获取模型的中波数成分；最后利用本文所述方法获取地下介质参数的高波数成分.通过不同方法逐步对地下介质参数细节部分进行成像，最终获取精确的地下介质成像.总的来说，即将LSRTM视为一种反演模型高波数成分的线性化波形反演.
(2)精细的近地表模型，特别是近地表品质因子模型的建立，对陆上数据LSRTM的应用至关重要.复杂的近地表结构严重影响地震波的传播并对随后的成像结果产生重要影响，目前较常用的初至走时反演、有限频反演方法等虽能建立较精确的速度模型，但是对品质因子的估计较为困难，而常规的品质因子估计方法一般精度有限.通过理论分析与比较，作者认为可将针对初至波的波形反演方法作为近地表刻画的主要反演工具，同时反演近地表速度及品质因子.
(3)将已知的地下介质分布信息作为约束条件加入成像算法中，如井数据约束、地层倾角约束等，通过约束条件的加入可在一定程度上减弱算法对背景模型的依赖性，增加算法的稳定性.总的来说，基于波形匹配的反演方法在理论上具有不可替代的优势，但在实际资料的处理过程中，还需要发展相应配套的方法技术及更稳定的LSRTM方法.针对上述问题完善算法本身并发展相应配套的新方法技术将是下一步的工作重点.
致 谢 作者感谢两位匿名审稿人的宝贵意见.感谢Madagascar软件(http://www.reproducibility.org)提供的绘图支持.借助式(10)关于数据空间标量乘算子的定义，伴随算子满足：
其中，L为波动方程算子，p为波动方程算子控制下的正演波场，
为波动方程算子L的伴随算子,
为伴随算子控制下的正演波场.将波动方程算子代入可得
在推导伴随算子的表达形式之前，先定义原始波动方程的初边值条件，初始条件满足静止初始条件即
对于自由地表边界条件，波场满足
对于吸收边界条件，即假设研究区域为无限大，引入Sommerfeld辐射条件，可得
对于式(A2)中第一项与时间偏导有关的项，通过分步积分可得
假设伴随波场满足终止条件
则式(A6)可整理为
对于第二项与空间求导有关的推导，为了表示方便，用松弛率函数Ψ(t)代替与黏滞性有关的项，即
将时间褶积与空间求导展开，并只考虑z方向，通过分步积分可得
假设伴随波场满足:
式(A10)可进一步整理为
当模拟无限大介质时，式(A11)对应吸收边界条件；当存在自由地表时，式(A11)对应自由地表边界条件.
对式(A12)利用
变换积分顺序并对符号做相应替换可得
若定义反因果函数
满足:
则式(A13)可整理为
同理对于x方向可得
将式(A8)、(A15)、(A16)代入式(A2)可得原始波动方程的伴随方程为
其中最右端项为震源项，将
显示写出，则上式可整理为
可以看到伴随波动方程与原始波动方程具有近似的形式，可以利用与原始波动方程近似的数值解法进行数值求解.
可以证明，在伴随波动方程控制下的波传播问题，沿时间正向传播的波场是增加的，数值求解是不稳定的；但是在满足伴随波场终止条件(A7)的条件下，沿时间逆向传播的波场是衰减的，数值求解是稳定的.通过引入伴随波动方程有效地避开了波场反向传播中的不稳定问题.
| [1] | Causse E, Mittet R, Ursin B. 1999. Preconditioning of full-waveform inversion in viscoacoustic media. Geophysics, 64(1): 130-145, doi: 10.1190/1.1444510. |
| [2] | Causse E, Ursin B. 2000. Viscoacoustic reverse-time migration. Journal of Seismic Exploration, 9: 165-184. |
| [3] | Dai W, Fowler P, Schuster G T. 2012. Multi-source least-squares reverse time migration. Geophysical Prospecting, 60(4): 681-695, doi: 10.1111/j. 1365-2478.2012.01092.x. |
| [4] | Deng F, McMechan G A. 2007. True-amplitude prestack depth migration. Geophysics, 72(3): S155-S166, doi: 10.1190/1.2714334. |
| [5] | Deng F, McMechan G A. 2008. Viscoelastic true-amplitude prestack reverse-time depth migration. Geophysics, 73(4): S143-S155, doi: 10.1190/1.2938083. |
| [6] | Dong S, Cai J, Guo M, et al. 2012. Least-squares reverse time migration: towards true amplitude imaging and improving the resolution. 82ed Ann. Internat Mtg., Soc. Expi. Gephys. Expanded Abstracts, doi: 10. 1190/segam2012-1488.1. |
| [7] | Hak B, Mulder W A. 2010. Migration for velocity and attenuation perturbations. Geophysical Prospecting, 58(6): 939-951, doi: 10.1111/j. 1365-2478. 2010. 00866.x. |
| [8] | Hak B, Mulder W A. 2011. Seismic attenuation imaging with causality. Geophys. J. Int., 184(1): 439-451, doi: 10.1111/j. 1365-246X. 2010. 04848.x. |
| [9] | Li Q Z. 1993. The Way to the Accurate Seismic Exploration (in Chinese). Beijing: Petroleum Industry Press. |
| [10] | Li Z C, Wang Q Z. 2007. A review of research on mechanism of seismic attenuation and energy compensation. Progress in Geophys. (in Chinese), 22(4): 1147-1152, doi: 10.3969/j. issn. 1004-2903. 2007.04.021. |
| [11] | Li Z C, Zhu X F, Han W G. 2008. Review of true-amplitude migration methods. Progress in Exploration Geophysics (in Chinese), 31(1): 10-15, 64. |
| [12] | Liu Y J, Li Z C, Wu D, et al. 2013. The research on local slop constrained least-squares migration. Chinese J. Geophys. (in Chinese), 56(3): 1003-1011, doi: 10.6038/cjg20130328. |
| [13] | Mittet R, Sollie R, Hokstad K. 1995. Prestack depth migration with compensation for absorption and dispersion. Geophysics, 60(5): 1485-1494, doi: 10.1190/1.1443882. |
| [14] | Mulder W A, Hak B. 2009. An ambiguity in attenuation scattering imaging. Geophys. J. Int., 178(3): 1614-1624, doi: 10. 1111/j. 1365-246X. 2009. 04253.x. |
| [15] | Nemeth T, Wu C, Schuster G T. 1999. Least-squares migration of incomplete reflection data. Geophysics, 64(1): 208-221, doi: 10.1190/1.1444517. |
| [16] | Qin N, Li Z C, Yang X D, et al. 2013. Pre-stack joint migration velocity modeling with multi-stage optimization. Progress in Geophys. (in Chinese), 28(1): 320-328, doi: 10.6038/pg20130135. |
| [17] | Ribodetti A, Virieux J. 1998. Asymptotic theory for imaging the attenuation factor Q. Geophysics, 63(5): 1767-1778, doi: 10. 1190/1.1444471. |
| [18] | Robertsson J O A, Blanch J O, Symes W W. 1994. Viscoelastic finite-difference modeling. Geophysics, 59(9): 1444-1456, doi: 10.1190/1.1443701. |
| [19] | Schuster G T, Wang X, Hang Y, et al. 2011. Theory of multisource crosstalk reduction by phase-encoded statics. Geophys. J. Int., 184(3): 1289-1303, doi: 10. 1111/j. 1365-246X. 2010. 04906.x. |
| [20] | Suh S, Yoon K, Cai J, et al. 2012. Compensating visco-acoustic effects in anisotropic reverse-time migration. 82ed Ann. Internat Mtg., Soc. Expi. Gephys. Expanded Abstracts, doi: 10.1190/segam2012-1297. 1. |
| [21] | Tang Y X. 2009. Target-oriented wave-equation least-squares migration/inversion with phase-encoded Hessian. Geophysics, 74(6): WCA95-WCA107, doi: 10.1190/1.3204768. |
| [22] | Tarantola A. 1984. Linearized inversion of seismic reflection data. Geophysical Prospecting, 32(6): 998-1015, doi: 10.1111/j. 1365-2478.1984. tb00751.x. |
| [23] | Tanrantola A. 2005. Inverse problem theory and methods for model parameter estimation. Society for Industrial and Applied Mathematics. |
| [24] | Wang H Z, Zhang L B, Ma Z T. 2004. Seismic wave imaging in visco-acoustic media. Science in China, 47(1): 146-154, doi: 10.1360/04za0013. |
| [25] | Wang Y F, Yang C C, Duan Q L. 2009. On iterative regularization methods for migration deconvolution and inversion in seismic imaging. Chinese J. Geophys. (in Chinese), 52(6): 1615-1624, doi: 10.3969/j. issn. 0001-5733. 2009. 06.024. |
| [26] | Wang Y H. 2002. A stable and efficient approach of inverse Q filtering. Geophysics, 67(2): 657-663, doi: 10.1190/1.1468627. |
| [27] | Wang Y H. 2006. Inverse Q-filter for seismic resolution enhancement. Geophysics, 71(3): V51-V60, doi: 10. 1190/1.2192912. |
| [28] | Wang Y H. 2008. Inverse-Q filtered migration. Geophysics, 73(1): S1-S6, doi: 10.1190/1.2806924. |
| [29] | Yang J L, Li Z C, Ye Y M, et al. 2008. Review of seismic illumination pre-stack depth migration methods. Progress in Geophys. (in Chinese), 23(1): 146-152 |
| [30] | Yao G, Jakubowicz H. 2012a. Least-squares Reverse-Time migration. 82ed Ann. Internat Mtg., Soc. Expi. Gephys.. Expanded Abstracts, doi: 10. 1190/segam2012-1425.1. |
| [31] | Yao G, Jakubowicz H. 2012b. Non-linear least-squares reverse-time migration. 82ed Ann. Internat Mtg., Soc. Expi. Gephys.. Expanded Abstracts, doi: 10.1190/segam2012-1435.1. |
| [32] | Zhang C J, Ulrych T J. 2010. Refocusing migrated seismic images in absorptive media. Geophysics, 75(3): S103-S110, doi: 10.1190/1.3374434. |
| [33] | Zhang J F, Wapenaar K. 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. |
| [34] | Zhang L B, Wang H Z. 2010. A stable inverse Q migration method. Geophysical Prospecting for Petroleum (in Chinese), 49(2): 115-120, doi: 10.3969/j. issn. 1000-1441.2010.02 002. |
| [35] | Zhang Y, Zhang P, Zhang H Z. 2010. Compensating for visco-acoustic effects with reverse time migration. 80ed Ann. Internat Mtg., Soc. Expi. Gephys. Expanded Abstracts, 3160-3164, doi: 10.1190/1.3513503. |
| [36] | Zhou H, Lin H, Sheng S B, et al. 2012. High angle prestack depth migration with absorption compensation. Applied Geophysics, 9(3): 293-300, doi: 10.1007/s11770-012-0339-z. 李庆忠. 1993. 走向精确勘探的道路. 北京: 石油工业出版社. |
| [37] | 李振春, 王清振. 2007. 地震波衰减机理及能量补偿研究综述.地球物理学进展, 22(4): 1147-1152, doi: 10. 3969/j. issn. 1004-2903. 2007. 04.021. |
| [38] | 李振春, 朱绪峰, 韩文功. 2008. 真振幅偏移方法综述. 勘探地球物理进展, 31(1): 10-15, 64. |
| [39] | 刘玉金, 李振春, 吴丹等. 2013. 局部倾角约束最小二乘偏移方法研究. 地球物理学报, 56(3): 1003-1011, doi: 10. 6038/cjg20130328. |
| [40] | 秦宁, 李振春, 杨晓东等. 2013. 叠前多级优化联合偏移速度建模. 地球物理学进展, 28(1): 320-328, doi: 10. 6038/pg20130135. |
| [41] | 王彦飞, 杨长春, 段秋梁. 2009. 地震偏移反演成像的迭代正则化方法研究. 地球物理学报, 52(6): 1615-1624,doi: 10.3969/j. issn. 0001-5733. 2009.06.024. |
| [42] | 杨敬磊, 李振春, 叶月明等. 2008. 地震照明叠前深度偏移方法综述. 地球物理学进展, 23(1): 146-152. |
| [43] | 张立彬, 王华忠. 2010. 稳定的反Q偏移方法研究. 石油物探, 49(2): 115-120, doi: 10.3969/j.issn. 1000-1441. 2010. 02.002. |