2014, Vol. 57
2. 美国加州大学圣克鲁斯分校地球与行星科学学院, 美国加州 995064;
3. 中国科学技术大学地球与空间科学学院, 合肥 230026
2. Earth and Planetary Sciences Department, University of California, Santa Cruz, CA, 995064, USA;
3. School of Earth and Space Sciences, University of Sciences and Technology of China, Hefei 230026, China
反演理论和方法是依赖于正演模型的.不同的正演模型有着不同的用于反演的模型参数并最终可能获得完全不同的反演结果.通过对地震波场的走时的分析，传统的层析成像反演 主要用于重建速度模型中速度变化的低频分量部分.为了获得更精确的反演结果，非线性的全波形反演(Pratt and Worthington, 1988, 1990; Woodward, 1992; 杨文采等，1993；杨文采和杜剑渊， 1994；常旭和刘伊克，1998；Pratt et al., 1998; Pratt, 1999; 刘伊克和常旭，2000；杨文采，2002；Liu et al, 2005, 2007; Woodward et al., 2008) 已经发展了很多年.但好的线性反演方法，例如散射层析成像反演，仍然存在着重要的研究价值，因为它能为有效的非线性反演方法提供理论依据.
不同于基于均匀背景速度的传统的Born模型(Devaney，1982，1984；Slaney et al., 1984; Wu and Toksöz, 1987)，在新发展的广义的散射层析成像反演方法中(Gelius, 1991; Schlottmann, 2006; Wu, 2007; Cao, 2008)，非均匀背景介质中(Schultz and Jaggard, 1987; Gelius, 1991)和非均匀体(Taylor, 1972)的散射波场都可以通过利用局部扰动波场的Born近似的方法来求得.
Wu(2007)推导出非均匀介质中体散射的散射层析成像反演公式，该公式包含对波场的有限频率带宽的校正和采集系统的有限空间孔径的校正.该公式使用模型的前向散射中重新归一化的格林函数，而这种正演的模型是基于De Wolf 近似的(De Wolf, 1971,1985;Wu, 1996,2003;Wu et al., 2007).De Wolf 近似是一种对散射波场进行多次前向散射和单次后向散射的近似.
本文详细地给出了主要由Wu(2007)建立和发展的基于非均匀介质中的局部Born模型体散射的层析成像反演和多频叠加的方法来提高谱域信息的质量和对速度模型的低波数分量进行恢复( Zhu and Wu, 2009)，该方法为能恢复精确的地下速度模型提供基础的理论依据.本文的安排如下，首先详细地给出基于非均匀介质中的局部Born模型体散射的层析成像反演的基本理论，详细地介绍波场的成像条件和局部成像矩阵、地震成像的分辨率算子、基于Born模型体散射的分辨率核函数表达式和局部波数域中的反褶积滤波公式以及两种用来计算Born模型的振幅校正因子的方法；其次，通过利用广义层析成像的反演方法对数值模型的谱的恢复和速度模型的低频分量的重建进行测试，通过对结果的分析给出一些初步的结论；然后将本方法与最小平方反演方法进行比较和讨论；最后总结本文和对进一步研究提出展望.
成像过程从数学上可以描述为波场的反向传播(或逆时传播)与聚焦(或成像条件)(Claerbout, 1971, 1985).在实际的地震数据中，由于受到系统的采集孔径和震源频率带宽的限制，即使使用精确的波场传播算子和成像条件，实际得到的成像振幅与地下反射体的真实的反射振幅往往相差很大.为了获得接近于真实的反射振幅就必须利用采集系统的分辨率核函数对成像振幅进行采集孔径和波场传播过程中的吸收衰减进行真振幅校正，而由于在空间域中，分辨率核函数的计算和反褶积滤波难以实现，所以本文在局部波数域对波场进行真振幅校正.在本节首先介绍空间域的成像条件，然后由它导出局部波数域的成像条件并对局部成像矩阵作简要的介绍，且空间域和局部波数域的成像条件的推导以及局部成像矩阵的表达式大多基于频率域的.
在对地震数据进行处理的过程中，叠前偏移是一种利用双聚焦算子进行成像的过程，它将来自震源排列的波场和接收器排列的波场都聚焦到成像点.这种成像方法相当于将采集系统从地表向下延拓至成像深度并进行局部的散射成像处理.这种双 聚焦的成像方法包括波场的反向传播以用来消除 上覆的非均匀背景介质的影响和运用成像条件(Claerbout,1971;Claerbout,1985)获得成像点处的像.在空间域中基于互相关形式表示的标准的成像条件的表达式可写为
其中
在式(1)、(2)、(3)和(4)中，L(ω,x)是频率域中的源的波场和接收器接收的波场在成像点x处双聚焦的结果，将频率域中所有频率的像L(ω,x)与震源的谱的乘积进行积分求和便可得到成像点x处的最终的像，也就是时间域的双聚焦成像.“Re” 表示复函数的实部；D(ω,x;xs)是从地表上xs处的源激发的入射波(或下行波)波场；U(ω,x;xs)是成像点x处的散射波场，它是通过对地面上接收到波场做反向传播并在散射点处对波场进行Rayleigh 积分而得到的.psc(ω,xg;xs)由地表上xs处的源激发的波场被地表上xg处的接收器接收到的频率域的散射波场；‘’ 表示复共轭；GI(ω,x;xs)和GI(ω,x;xg)是成像过程中的格林函数，它们有别于正演过程中的格林函数；s(xs)是震源的空间分布且通常取为单位值；f(ω)是震源的谱函数，B(ω)是震源的频率域的带宽；A(xs)和A(xg)分别是震源和接收器的空间域的孔径.成像的振幅严重地受到系统采集孔径的影响，包括震源和接收器的孔径以及震源的带宽.
相对于经典的散射层析成像的成像条件来说，为了使得成像过程中来自震源的排列的波场聚焦和来自接收器的排列的波场聚集有着很好的对称性，本文将空间域中的基于互相关形式的成像条件(式(2))稍微地修改成下面的形式
式(5)是运用在散射层析成像中的成像条件(Wu and Toksöz, 1987; Wu, 2007)， 所以被称为层析成像的成像条件.若使用归一化的震源波场，层析成像的成像条件可以写为
为了得到成像场的局部角度谱的信息，传统的空间域的成像条件可以扩展到局部波数域(或小波束域)(Wu and Chen, 2002)，在局部波数域中的成像函数就不再是一个标量的值而是一个矩阵：局部成像矩阵 (Local Image Matrix, LIM).在此，本文定义局部波数域中源的波数向量为从成像点指向地表上的源的方向，它与波场的入射波数向量正好方向相反.因此，在频率和局部波数(或角度)域中，基于互相关形式的标准的成像条件可表示为
在式(7)和(8)中，
是在成像点x处的局部波数域的局部成像矩阵；
是成像点x处从地表上xs处的源激发的局部波数域的入射波(或下行波)的波场；
是成像点x处局部波数域的散射波(或上行波)的波场，它是通过Rayleigh积分将地表上接收器接收到的波场反向传播并聚焦到散射点而得到的.其中，ks=
分别是源和接收器的波数向量，且它们也可以表示成
分别是波数向量的水平和垂直分量，且有
分别是波数域中源和散射方向的单位向量；
是来自地面上xs处的点源入射到成像点x处的小波束域的入射波场的格林函数
是地面上xg处的接收到的成像点x处的小波束域的散射波场的格林函数.同理，可以得到局部波数域中以互相关形式表示的层析成像的成像条件如下：
在本文中, 将式(1)中的格林函数在深度z处的成像点 x′ =(x′,z′)附近分解到局部波数域，则该成像点附近任一点的波场可表示为
在式(10)和(11)中， r′ =(x′－x0,z′－z)，x0是分解位置处的坐标原点，通常选
xs)和
分别是不含相位因子的局部入射和反射的平面波系数，它们分别是由对格林函数
进行局部小波束域分解而得到的.求和变量‘s’和‘g’ 分别表示对局部入射和反射的平面波求和.除此以外，还可以利用局部倾角叠加(Local Slant Stack)的方法(Wu and Chen, 2002; Xie and Wu, 2002; Xie et al., 2003)来将格林函数分解到局部波数域.
将式(10)和(11)代入式(9)，在成像点 x′ 处的局部波数域的局部成像矩阵
可表示为
从式(12)可知，来自于从任一局部入射和散射角对的信息只能得到谱域的分量在
处的波场的局部谱信息.在图1中
是局部入射波数向量，由于本文定义源的波数向量
为从成像点指向地表上的源，所以有
.若设z′=z，则式(12)可以简化为
 | 图1 几种波数向量的定义和当源和接收器的排列均为无穷长时的Born模型的谱域的覆盖范围的示意图 Fig.1 Definitions of local wavenumbers and domain of spectral coverage for the Born model when the source and receiver arrays are both infinite |
如果利用多频率的数据对散射体进行成像，这时谱域的覆盖范围就可以得到扩展且成像点 x′ 处局部波数域的局部成像矩阵可表示为
如果在反演过程中所用的算子与正演(采集)过程中的格林函数完全一致，或它们的相位信息相同，那么成像点处的局部成像矩阵就是实的矩阵.然而，实际上由于采用近似的传播算子和速度模型的不精确，总体上，局部成像矩阵是复矩阵，可表示为
其中，
分别是成像点x处的局部波数域的局部成像矩阵的振幅和相位矩阵元素.若背景的速度模型存在误差，则对于不同的源和接收器的局部波数向量对
,来说它们的相位也会有差别，而这些相位的残差会导致最终成象结果的模糊，并存在一些假象和相关噪声.
本节主要详细地给出基于Born模型体散射的广义层析成像反演的理论和公式(Wu, 2007; Zhu and Wu, 2009; 朱小三，2010).
对模型的正演和成像的过程可以用算子的形式表示为
在式(16)和(17)中，F是模型的采集(或正演)算子，它是利用采集系统得到模型S的正演数据U；B是模型的成像算子，它将正演数据U反演成模型的像I；x0是模型空间的散射点，而x是模型空间的成像点.则模型的分辨率算子(或分辨率矩阵)可表示为
如果模型的成像算子是模型的正演算子的精确的逆算子，则分辨率算子R将是一个单位算子.但在通常情况下，分辨率矩阵并不是一个单位矩阵，且在矩阵中沿对角线周边分布的矩阵的各元素可以对模型的成像或反演的参数分辨率提供定量的描述.从式(17)和(18)可知，模型的像I可以看成分辨率算子R和模型S的褶积.
若使用归一化的源，则可以得到基于式(1)和(2)的互相关成像条件的成像算子的表达式
同理，也可以得到基于式(1)和(5)的层析成像的成像算子
模型的正演算子从总体上来说是非线性的，也就是模型中某一处的因模型参数的扰动而引起的响应将会受到其他位置处参数扰动的影响(如多次散射).在对模型的真实的反射振幅进行成像的方法中，总是利用优化的方法将问题线性化，这样处理不但能够恢复优化的地层反射率，而且还可以从基础理论上发展有效的迭代反演的方法.所以，本节中将介绍一种线性化的模型算子：Born模型.
在Born模型中，各个散射体都是独立的且不考虑它们之间的相互作用.这种近似只适用于模型中参数存在小的扰动或波场的传播距离比较短的情况.在Born模型中待反演的是相对于给定的背景介质的未知参数的体扰动，而模型中的格林函数是背景介质中的脉冲响应.对于标量 (声波和常密度)介质情况，反演的目标函数可用相对于背景速度c0(x)的速度分布c(x)来表示，其中x是模型中点的位置向量.从地表上xs处激发的源被地表上xg处接收器接收到的基于Born近似的散射波场可表示为
其中
在式(21)中，
是散射波场；GF(x;xs)和GF(xg;x)是对非均匀背景介质正演过程中的标量波场的格林函数，GF(x;xs)是在x处接收到的地表上xs处的源的波场对模型中非均匀体的响应，而GF(xg;x)是在地表上xg处接收到的由x处的波场对模型中非均匀体的响应；O(x)是模型的速度扰动函数，它是经典的散射层析成像的目标函数(Devaney, 1982, 1984; Slaney et al., 1984; Wu and Toksöz,1987)；k0=ω/c0(x)是背景介质的波数；V是非均匀介质积分的体积.式(21)可以直接从Lippmann-Schwinger方程中导出(Wu and Toksöz, 1987).因此，Born模型的正演算子可表示为
本文使用正演散射波场对模型中的格林函数进行重新归一化，且正演模型是基于De Wolf 近似的(De Wolf,1971;De Wolf,1985;Wu, 1996,2003;Wu et al., 2007)， 它通过对正演过程中的波场进行多次前向散射和单次后向散射的近似来更新格林函数.相对基于均匀背景介质的经典的散射层析成像反演的Born模型，基于式(21)的定义Born近似的散射波场和De Wolf近似的格林函数的模型可以称为De Wolf-Born模型(或局部Born模型).
由于经典的Born模型只适用于弱散射，即弱的参数扰动和积分体积比较小的非均匀介质.而对于强扰动的介质或大的非均匀体利用Born散射来进行反演是不合适的.但在De Wolf-Born模型中，入射波场和格林函数都是利用前向散射的波场重新归一化的，所以可以将它扩展到对参数的体扰动较强的模型进行反演，并且它能克服一些经典的Born模型的限制条件和不足之处.
模型中分辨率算子的核函数(或分辨率核函数) (Back us and Gilbert, 1967, 1968, 1970; Tarantola, 1984a, 1984b, 2005) 可以通过式(18)、(19)和(23)来获得，且基于互相关成像条件的Born模型的分辨率核函数可表示为
分辨率核函数是将模型空间中的散射点映射回模型空间的一种算子，也称为成像系统的点散射函数(Point Spreading Function,PSF)(Wu and Chen, 2001, 2006; Wu et al, 2004,2006,2008; Wu and Cao, 2005; Chen et al., 2006; Wu and Mao, 2007; Wu, 2007; Cao, 2008; Wu and Xie, 2009; 朱小三，2010)，它既受到模型中数据的采集(或正演)的影响也受到数据的反演(或成像) 过程的影响.系统的分辨率核函数的质量是由整个成像系统的分辨率来决定的.由于分辨率核函数模拟模型中散射点波场的采集和偏移成像的过程，所以，在模型空间x0的散射点处产生散射波场数据的正演格林函数对可表示为
除式(25)中的格林函数对以外，式(24)的其他部分的算子将波场数据反向传播回模型空间，并运用成像条件对散射能量进行重新聚焦.
同理，基于Born 模型的层析成像的成像条件的分辨率核函数可表示为
为了消除成像场中采集系统的影响，例如，采集孔径和接收带宽的限制，波场传播的路径效果和不精确的成像算子，可利用分辨率核函数对成像场进行反褶积滤波来处理.若可以计算分辨率核函数，则真实的模型可以通过利用分辨率核函数对成像场进行反褶积滤波(并做相关的归一化处理)处理而得到.但由于计算空间域的分辨率核函数和空间域的体反褶积滤波是比较费时且难以处理，而在局部角度域中的反褶积滤波容易实现且具有较高的效率.将分辨率核函数分解到局部角度域与对散射算子进行特征函数分解比较相近.对于位于地表的地震反射波的观测系统，由于受到采集孔径的限制，其照明信息中的受反射体的倾角的限制是其观测系统固有的缺陷.但在局部角度域，采集系统的照明分布是很清楚的.因此，本文将分辨率核函数分解到局部波数域(或小波束域)并在该域中利用分辨率核函数对成像波场进行反褶积滤波处理.本文对分辨率核函数R(x;x0)(见式(24))在以x0为坐标原点的空间域中进行局部的三维傅里叶变换(Wu et al., 2006)即可得到局部波数域中基于互相关形式的分辨率核函数表达式
其中
在式(27)和(28)中，
是通过对波场进行反向传播并积分而得到的，而该波场是利用有限孔径的采集系统作用于模型的正演的格林函数而得到的.对于地表上xs处的某一确定的源和局部波数域中某一确定的局部接收波数向量
，本文可以利用来衡量某一具体的小波束域波场的成像系统采集孔径的效果.同理，本文对空间域中基于层析成像的成像条件的分辨率核函数RT(x;x0)(见式(26))进行三维傅里叶变换得到其在局部波数域的表达式(Chen et al., 2006)
在式(29)中
的表达式与式(28)相同.
在局部波数(或角度)域中，反褶积滤波是成像矩阵
被分辨率核函数
相除的过程(Wu and Chen, 2001, 2006; Wu et al., 2004; 2008; Wu and Cao, 2005; Chen et al., 2006; Wu and Mao, 2007; Wu and Xie, 2009),可表示为
所以，可以通过对小波束域的波场进行反演来重建空间域的目标函数
式(30)和(31)是分别对最终的多频的局部成像矩阵进行局部波数域的反褶积滤波和空间域的目标函数的重建.当然，本文也可以在局部波数域中对单个的频率的局部成像矩阵进行反褶积滤波和空间域的目标函数的重建.二维空间域的单频的目标函数的重建表达式为
其中
在式(32)、(33)和(34)中，
是将坐标从
的Jacobian系数;
是层析成像的滤波算子，所以，在局部波数域中直接将层析成像的滤波算子运用于成像矩阵便可得到目标函数.因此，空间域的目标函数可通过对所有的单个频率的目标函数(式(32))进行频率域的求和而得到.除此之外，本文还可以对单个频率的成像波场在局部波数域中进行滤波，然后在谱域将所有单个的频率的目标函数的谱进行叠加平均从而得到平均的谱，再利用目标函数的平均的谱来重建空间域的目标函数.
一般情况下，成像矩阵和分辨率核函数都是复矩阵，且它们的剩余相位都是局部波数的函数.利用多频的数据对局部波数域中的谱进行恢复以及对成像矩阵进行反褶积滤波是利用层析成像进行反演的普适方法.若在反演过程中，使用真实的背景速度模型，且传播算子不存在相位误差，本方法仅需要校正成像场的局部角度谱的振幅便可对成像的振幅进行真振幅补偿.
对于大多数的震源子波(如Ricker子波)，波场所携带的能量主要分布在震源子波的主频附近.如果只利用主频处的振幅校正因子而不是所有频率的振幅校正因子来对波场进行校正，这对于完全的分辨率核函数来说是一个合理的近似，且可以节省大量的计算(Wu et al., 2004; 2008; Wu and Cao, 2005; Xie et al, 2005, 2006; Wu and Chen, 2006; Chen et al., 2006; Xu and Lambaré, 2006; 2009; Wu and Mao, 2007; Wu and Xie, 2009; Zhu and Wu, 2010).本文定义模型的振幅校正因子如下：
对于单一频率的分辨率核函数来说，对于一给定的源的局部波数向量，则局部波数向量
和接收器的局部波数向量
成为一一对应的，因此有
可以用
来衡量对于地表上xs处的源和主频为ω0的某一给定的局部波数向量 的地震数据的采集孔径.在理想的情况下，如果使用真实的速度模型，聚焦算子的格林函数(反向传播并积分)是正演格林函数的复共轭，则
应是一个实函数.然而，在通常的情况下，
是复函数.本文使用
并除掉其相位信息.同样，也对取绝对值，因此有
系统的采集孔径积分包含了对由点源激发的波场被地表接收器排列接收并将接收到的波场在小波束域反向传播回成像点处的过程.根据互易定理，该过程等于将小波束域波场(由接收器的局部波数向量定义的小波束域波场)传播到地表的接收器排列并将接收器接收到的波场利用有限的波前传播孔径进行反向传播的过程.对于反向传播算子，本文大多使用能保持能量守恒的格林函数.这样，在对波场重新聚焦的过程中，所有的能量的损失将被忽略，例如，波场的边界反射、P波和S波的相互转换、散射和黏弹性衰减等.所以，可以将接收器排列接收到的波场能量集中到最大程度，这也就是利用单程波场聚焦算子的特点.若利用
来代替
，除了一些边界附近的小波束域波场以外，则能量守恒可表示为
根据互易定理，同样可以通过在地表xg处的点源激发的波场在模型中x0处以同样的接收角和小波束宽度接收到的小波束域的波场来计算格林函数.式(38)中的近似避免了对波场的正传和反传所需的计算量，提高了计算效率.且式(38)中右边的表达式只包含模型的采集系统信息，类似于在直接照明分析中使用的计算方法.
所以，本文可以利用格林函数来计算模型的分辨率核函数.对于大多数情况，并不知道真实的格林函数或需要大量的计算才能得到，在上面的式(37)中，本文利用下面的近似
则基于标准互相关成像条件的振幅校正因子的表达式可简化为
同理, 基于层析成像的成像条件的振幅校正因子的表达式为
若利用局部波数变量
和局部角度变量
之间的相互关系
，本文可以得到基于标准互相关成像条件的局部角度域的振幅校正因子的表达式
以及基于层析成像的成像条件的局部角度域的振幅校正因子的表达式
在式(42)和(43)中，角是由局部波数向量
分别由源的局部波数向量和接收器的局部波数向量来定义,且有
通过应用上面的近似，对振幅进行校正的计算方法变得更加方便有效.
由于受到地震数据的采集孔径和频率带宽的限制，即使采用精确的波场传播算子和准确的成像条件，实际的成像振幅与地下反射目标体的真正的振幅仍相差很远.利用多频数据的一种方法就是估计积分角度，若利用关系式
将分辨率核函数的表达式从局部波数域转换到局部角度域，其中
是波场反向传播方向的单位向量(与偏移的倾角方向相反)，本文可以将同属于某一角度的所有频率的值进行叠加求和便可得到该角度束的分辨率核函数
同理，本文也可以将成像矩阵分解到局部角度域并得到以
表示的局部角度道集
在局部角度域中，利用分辨率核函数对成像矩阵进行反褶积滤波，从而得到成像点处更精确的角度谱.最终的像是将相应的角度域滤波之后的成像矩阵对角度进行积分求和而得到的
式中，ε是用于归一化的阻尼因子.
本文分别利用Born近似和有限差分的方法对选 定的模型产生炮集数据.目标函数O(x)是相对于背 景速度的扰动函数，而背景速度是只沿深度方向变化 的v(z)介质，可表示为c0(z)=3.5+0.12z (km/s). 在图2a中，模型中间的方块是速度异常体，异常体的尺寸为1 km×1 km，本文使用该异常体相对于背景速度的平均速度扰动分别为+23% (见图2a)和+30%两种速度模型.对这两种速度模型本文也利用有限差分的方法直接对二维的标量波方程(见式(47))求解来产生炮集数据
u(x,t)是地震波场；f(t)是位置在 xs处的震源时间函数；δ(x－ xs)是Dirac delta函数.
 | 图2 (a)用于产生正演炮集数据的含有方块异常体的速度模型.方块异常体相对于背景速度的平均速度扰动为 23%. (b)具有不同主频(如0.3, 0.5, 0.75, 1.13, 1.70, 2.56, 3.86和5.81 Hz)的Ricker子波的谱和它们的叠加的谱. (c) 当模型中速度的扰动是实函数时，多频的Born模型数据的谱域的覆盖范围示意图 Fig.2 (a) Velocity model of a square box with average velocity perturbations 23% higher than the background velocity. (b) Spectra of Ricker wavelets with different dominant frequencies (i.e., 0.3, 0.5, 0.75, 1.13, 1.70, 2.56, 3.86 and 5.81 Hz) and the final stacked spectrum. (c) Domain of multi-frequency spectral coverage of local image matrices for the Born model when the velocity perturbation is a real function in the model |
本文使用8个具有不同的主频的Ricker子波作为震源子波来产生炮集数据.这些子波的主频分别为 0.3, 0.5, 0.75, 1.13, 1.70, 2.56, 3.86 和 5.81 Hz (见图2b), 且这些子波的主频的选取是根据下面的公式(Sirgue and Pratt, 2004)
fn是第n个Ricker子波的主频；hmax/z是最大偏移距的一半与深度z的比率.
利用局部平面波分解方法将波场分解到局部波数域(或局部角度域)，并在局部波数域对来自源和接收器的局部平面波在成像点处利用成像条件，便可以获得模型中成像点x处的局部成像矩阵L(ω,x,)和局部分辨率矩阵R(ω,x,).在本文的数值测试中，只计算模型中每个点的单个频率的局部成像矩阵和相应的局部分辨率矩阵，而这些频率的选取是根据各炮集数据的主频来确定的.由于本文只计算8个不同Ricker子波(见图2b)的炮集数据，且可以节省大量的计算所以也只计算8个单个频率的局部成像矩阵和其相应的局部分辨率矩阵.
对于频率域中每个局部入射和接收角对
，只能得到目标函数的谱在一个分量
上的信息(见图1)，所以，来自一个局部入射角的所有散射角的谱只能覆盖谱域信息的一半区域.但当模型中速度的扰动为实函数时，则目标函数在谱域中的信息是共轭对称的(式(40)和(41))(Wu and Toksöz, 1987)，所以可以根据谱域中已知的那一半区域的信息来得到谱的另一半区域的信息, 即可以根据下面式子来分别恢复局部成像矩阵和局部分辨率矩阵中另一半区域谱的信息
然而，对于单个频率的局部成像矩阵来说，即便恢复了谱域的另一半区域的信息，在谱域中仍然有两个圆形的盲区.为了将这两个盲区的信息最大限度地填充起来，本文采用对局部成像矩阵的多频的谱进行平均来获得盲区的信息并提高谱域的信息的均一性(见图2c)，而这些具有多频谱信息的局部成像矩阵是由对不同主频的炮集数据处理而得到的.若采用更多的频率的局部成像矩阵的谱进行平均，则谱域的信息将会恢复得更完全.
本文选择图3a中用箭头指出的三个点利用声波数据(利用有限差分方法对二维的标量波动方程进行求解而得到的数据)来产生的局部成像矩阵的原始的谱和恢复校正后的谱进行比较(见图3b—3e).且在图3b—3e中显示的三个点谱的顺序与图3a中的对应的点的顺序一致.从图3b—3e可以看出，利用多频的谱叠加平均，局部成像矩阵的谱在盲区的信息得到了很大程度上的恢复，可以看出对于不同的平均速度扰动(23%和30%)，恢复后的所有的局部波数域的谱的质量无论从覆盖的面积范围还是幅值的均一性都有了很大的提高，特别是对于谱的低波数分量部分，即恢复校正后的局部成像矩阵的谱有着更好的覆盖范围，其幅值也有着更好的均一性.
 | 图3 局部成像矩阵的原始的谱和恢复校正后的谱，这些局部成像矩阵是利用基于局部余弦基的单程波传播算子运用于声波数据(利用有限差分方法对二维的标量波动方程进行求解而产生的数据)而得到的 (a)在速度模型中选取的显示谱恢复的三个点.从(b)到(e)中局部成像矩阵的原始的谱和恢复校正后的谱选取的顺序与(a)中的三个点相一致.(b) 利用模型中平均速度扰动为23%速度模型产生的局部成像矩阵的原始的谱；(c) 对(b)中的原始的谱利用相应的局部分辨率矩阵校正后的谱； (d) 利用模型中平均速度扰动为30%速度模型产生的局部成像矩阵的原始的谱；(e) 对(d)中的原始的谱利用相应的局部分辨率矩阵校正后的谱. Fig.3 Original spectra and their corresponding recovered spectra of local image matrices of the three selected points in the velocity model. Those local image matrices are generated using LCB one-way propagator to the shots data generated using finite difference method to the exact 2D scalar wave equation (a) Sampled three points in the velocity model for showing the spectral recovery. (b) Original spectra of LIM generated using the velocity model with average velocity perturbation 23%. (c) Recovered spectra using their corresponding local resolution matrices of (b).(d) Original spectra of LIM generated using the velocity model with average velocity perturbation 30%. (e) Recovered spectra using their corresponding local resolution matrices of (d). |
由于本文选取产生炮集数据的Ricker子波的主频是从0.3~5.81 Hz，所以本方法只能恢复速度模型中的低频分量.从图4可知，对于模型中不同的速度扰动(分别是23%和30%)，利用基于局部余弦基的单程波传播算子作用于Born数据进行速度模型重建的结果基本上一致，两者的垂向和水平方向的截面图都很相似.从利用声波数据(利用有限差分方法对二维的标量波动方程进行求解而得到的数据)产生的速度重建的结果(见图5)可以看出，除了图5d中垂向上的速度值有一些畸变之外，该结果对模型中速度异常体无论是在垂向还是水平方向都能很好地恢复.重建后的速度值在异常体的真实位置上下震荡可能是由于缺少足够的低波数分量信息和波场在高速的方块异常体内产生的多次波引起的.
 | 图4 基于局部余弦基的单程波传播算子作用于Born数据的速度重建的结果和结果的横截面图 (a) 对模型中相对于背景速度的平均速度扰动为23%速度重建的结果；(b) 对模型中相对于背景速度的平均速度扰动为30%速度重建的结果； (c)和(e) 分别是结果(a)的垂直和水平截面；(d)和(f) 分别是结果(b)的垂直和水平截面. Fig.4 Results for the velocity reconstruction using LCB one-way propagator to the Born data and cross-sections through those results (a) Velocity reconstruction for the model with average velocity perturbation 23% with respect to the background velocity; (b) Velocity reconstruction for the model with average velocity perturbation 30% with respect to the background velocity; (c) and (e) are the vertical and horizontal slices through the results of (a), respectively; (d) and (f) are the vertical and horizontal slices through the results of (b). |
 | 图5 基于局部余弦基的单程波传播算子作用于有限差分数据(利用有限差分方法对二维的 标量波动方程进行求解而产生的数据)的速度重建的结果和结果的横截面图 (a) 对模型中相对于背景速度的平均速度扰动为23%速度重建的结果；(b) 对模型中相对于背景速度的平均速度扰动为30%速度重建的结果；(c)和(e) 分别是结果(a)的垂直和水平截面；(d)和(f) 分别是结果(b)的垂直和水平截面. Fig.5 Results for the velocity reconstruction using LCB one-way propagator to the shots data generated using finite difference method to the exact 2D scalar wave equation and cross-sections through those results (a) Velocity reconstruction for the model with average velocity perturbation 23% with respect to the background velocity; (b) Velocity reconstruction for the model with average velocity perturbation 30% with respect to the background velocity; (c) and (e) are the vertical and horizontal slices through the results of (a), respectively; (d) and (f) are the vertical and horizontal slices through the results of (b). |
综上所述，在本测验中广义散射层析成像反演对介质模型中速度扰动较小的情况适用.该方法可以对相对于背景速度的平均速度扰动达到23%的速度异常体的低波数分量进行很好地恢复.
本文利用同图2a所示的速度模型同样的方法对Marmousi速度模型(见图6a)产生炮集数据.图6a中的Marmousi速度模型的长度为7450 m, 深度为3000 m. 本文采用8种主频分别是0.3、0.4、0.61、1.23、2.49、5.05、8.10和10.27 Hz的Ricker子波对Marmousi速度模型产生二维标量波动方程(见式(46))的有限差分炮集数据. 炮点和检波器均位于地表且从左自右排列，炮点间隔为50 m, 检波器间隔为 20 m.
 | 图6 以多尺度表示的Marmousi速度模型及其谱 (a)和(h)分别是用于正演二维标量波方程的炮集数据的Marmousi速度模型和它的谱；从(b)到(g)分别是速度模型的从第6到第1个尺度的表示；从(i)到(n)分别对应着从(b)到(g)的谱. Fig.6 Marmousi velocity model, its spectrum and their corresponding multi-scale representations (from 1 to 6 scales) generated by multi-resolution analysis using dyadic scaling (a) and (h) are the exact Marmousi velocity model for generating shots data using finite difference method to the exact 2D scalar wave equation and its spectrum, respectively. On the left column, from (b) to (g) are the corresponding multi-scale representations from 6 to 1 of the Marmousi velocity model, respectively; on the right column, from (i) to (n) are the corresponding spectra of those from (b) to (g) on the left column, respectively. |
本文利用基于局部余弦基的单程波传播算子的方法和两种不同的背景速度模型(见图7a和7c)对产生的炮集数据进行广义散射层析成像反演.为了详细地比较广义散射层析成像反演对Marmousi速度模型的低波数分量进行恢复的结果(见图7b和7d)与真实Marmousi速度模型(图6a)的符合程度，本文对Marmousi速度模型(见图6a)和两种层析成像反演的结果(见图7b和7d)都进行多尺度表示(见图6、8和9)(即谱域的二进制尺度的多分辨率分析)(Yu, et al., 2004). 在本文中使用曲波变换(Curvelet Transform) (Candès and Donoho, 2000; Candès and Demanet, 2005)的方法来对它们进行多尺度分解，都分别分解成6个尺度表示.图6b—6g和6i—6n是以6个尺度(即谱域的二进制尺度)表示的Marmousi速度模型及其对应的谱. 图7a是利用对速度只沿深度均匀增加的v(z)背景速度模型, 而图7c是对真实的速度模型进行严重的平滑而产生的背景速度模型.图7b和7d分别是利用图7a和7c这两种背景速度模型对Marmousi速度模型的低波数分量进行恢复的结果.
 | 图7 利用基于局部余弦基的单程波传播算子作用于炮集数据来进行速度模型重建 (a) 速度只沿深度方向递增的v(z)背景速度模型；(b) 利用背景速度模型(a)进行速度重建的结果；(c) 对真实速度模型进行严重平滑的背景速度模型；(d) 利用背景速度模型(c)进行速度重建的结果. Fig.7 Velocity reconstruction using LCB one-way propagator to the shots data generated using finite difference method to the exact 2D scalar wave equation (a) Background velocity model with velocity changes only in the depth direction (a v(z) velocity model); (b) Velocity reconstruction from the background velocity model (a). (c) Heavy smoothed exact Marmousi velocity model. (d) Velocity reconstruction from the background velocity model (c). |
 | 图8 以多尺度表示的Marmousi速度模型和利用速度只沿深度方向递增的v(z)背景速度模型进行层析成像反演结果从(a)到(f)分别是Marmousi速度模型的从第6到第1个尺度的表示；从(g)到(l)分别是利用速度只沿深度方向递增的v(z)背景速度模型进行层析成像反演结果的从第6到第1个尺度的表示. Fig.8 Compare the different scale representations (from 1 to 6 scales) between the exact Marmousi velocity model (Fig.7a) and the result generated using vertical gradient background velocity model (Fig.7b) On the left column, from (a) to (f) are the corresponding scale representations from 6 to 1 of the exact Marmousi velocity model, respectively; on the right column, from (g) to (l) are the corresponding scale representations from 6 to 1 of the result generated using vertical gradient background velocity model, respectively. |
 | 图9 以多尺度表示的Marmousi速度模型和利用对真实速度模型进行严重平滑的背景速度模型进行层析成像反演结果从(a)到(f)分别是Marmousi速度模型的从第6到第1个尺度的表示；从(g)到(l)分别是利用对真实速度模型进行严重平滑的背景速度模型进行层析成像反演结果的从第6到第1个尺度的表示. Fig.9 Compare the different scale representations (from 1 to 6 scales) between the exact Marmousi velocity model (Fig.7c) and the result generated using heavily smoothed Marmousi velocity model (Fig.7d) On the left column, from (a) to (f) are the corresponding scale representations from 6 to 1 of the exact Marmousi velocity model, respectively; on the right column, from (g) to (l) are the corresponding scale representations from 6 to 1 of the result generated using heavily smoothed Marmousi velocity model, respectively. |
图8是以6个尺度表示的Marmousi速度模型(图8a—8f)和利用速度只沿深度方向递增的v(z)背景速度模型进行层析成像反演结果(图8g—8l)；图9是以6尺度表示的Marmousi速度模型(图9a—9f)和利用对真实速度模型进行严重平滑的背景速度模型进行层析成像反演结果(图9g—9l).
从图6可以看出，速度模型的不同尺度的表示与其在谱域的带宽相联系，尺度越大的其对应谱域的能量越集中.从图7b可以看出，利用速度只沿深度方向递增的v(z)背景速度模型进行层析成像反演结果在除模型的底部以外对速度模型的其他部分的恢复结果还是很好，且受高频干扰较少；由图8可知，反演结果能很好地恢复Marmousi速度模型的第4，5和6个尺度的分量和部分地恢复第3个尺度的分量，在波数域有着比利用光滑层状背景速度模型反演的结果更好的带宽，故它也更接近于真实模型.从图7c—7d可以看出，利用对真实速度模型进行严重平滑的速度模型作为背景模型(图7c)进行反演的结果，除了模型的上部有一些畸变外，模型的其他部分都得到了较好的恢复，尤其是模型的下部的速度体得到比较清晰的恢复；由图9可知，反演结果能很好地恢复Marmousi速度模型的第3，4，5和6个尺度的分量和部分地恢复第2个尺度的分量，由于对模型谱域不同带宽的截断，所以在不同尺度 表示的边界附近存在着由吉布斯效应引起不同程 度的干扰，该结果明显好于图8所示的结果，对 Marmousi速度模型的低波数分量进行了近于完全地恢复.
当然，从图8和9可以看出，利用三种背景速度模型反演的结果都没有对Marmousi速度模型的高波数分量进行很好地恢复，这是由于本文中使用的广义散射层析成像反演的方法是基于Born模型的，而Born近似只适用于地震数据的低频分量，因此，本方法很难对速度模型的高波数分量进行恢复.若使用基于边界散射的广义散射层析成像反演的方法(Wu, 2007)便能够对速度模型的高波数分量进行很好地恢复.故结合这两种散射层析成像反演方法，本方法便有可能对速度模型的波数域的分量进行完全地恢复.
以上的对所有的速度模型的低波数分量的重建都是利用8个主频均为低频的Ricker子波产生的炮集数据，且只利用了8个不同的频率进行叠加对原始速度模型的低波数分量进行一次重建的结果，如果采用更多的频率、更宽的频带的数据进行叠加和迭代，本方法会得到更好的速度重建的结果.
最小平方反演也是一种能对成像矩阵进行真振幅校正的反演方法，而其中所用到的Hessian矩阵与本文中的分辨率核函数有着诸多的异同点.在最小平方反演中(Lailly, 1983; Beylkin, 1985; ten Kroode et al., 1994; Chavent and Plessix, 1999; Hu et al., 2001; Kuhl and Sacchi, 2003; Rickett, 2003; Clapp et al., 2005; Valenciano et al., 2006; Beylkin et al., 1986)，成像算子(或反向传播算子)是最小平方误差的函数，而 Hessian矩阵则是误差函数的二阶梯度，可表示为式中，J(o)是误差函数，o(x)是反演的目标函数，dm(xs,xg,o,ω)和d(xs,xg,ω)分别是某个确定模型的频率域的合成数据和真实的数据.Hessian矩阵包含了数据系统的采集孔径、速度模型、频率域 带宽以及透射能量损失等诸多信息.而假设式(51)中交叉项很小以至于可以将其忽略掉，则式(51)可写为
式中，f(ω)是震源的谱函数，GM(ω,x;xs)和GM(ω,x;xg)都是正演过程中的Green函数.通过比较Hessian矩阵的表达式和分辨率核函数的表达式，我们可知：它们有着相同的正演模型算子，但成像算子各异.对于Hessian 矩阵，其成像算子是模型 算子的伴随算子，而在分辨率核函数中，其成像算子是聚焦算子(或反向积分算子)，它可以利用任何近似的Green函数来进行计算.由于基于Born模型的模型算子中有k2因子，故模型算子是一个伴随算子，但聚焦算子并没有该因子.同理，若模型算子的Green 函数存在一个幂指数的衰减项，则其伴随算子必定有相同的衰减项，而对于聚焦算子则不一定要有衰减项，但它必须与模型算子在相位上是复 共轭的.而在最小平方反演中由于利用Hessian矩阵进行反演，这些共轭的相位信息则互相抵消，这就导致在强衰减的介质中，伴随算子的衰减项会在很大程度上衰减掉弱的散射信号，从而导致反演结果的不稳定.当系统的采集孔径为无限大时，Hessian矩阵一般呈对角化且能代表几乎所有的系统照明信息(Beylkin, 1985; ten Kroode et al., 1994; Chavent and Plessix, 1999; Hu et al., 2001; Kuhl and Sacchi, 2003; Rickett, 2003; Clapp et al., 2005)， 但它无法区分开各个倾角上的照明信息(Valenciano et al., 2006; Beylkin et al., 1986).而分辨率核函数则能很好地表示不同倾角上的照明信息.
由于本文的内容主要是吴如山和M. Nafi Toksöz的1987年的文章(Wu and Toksöz, 1987)的继承和发展，在此就该文章与本文的内容的几个方面的异同点进行讨论如下：首先，前者的散射层析成像公式的推导过程是基于均匀背景介质的，而本文由于利用De Wolf 近似对前向散射的Green函数进行了归一化校正，虽然在每个小的高斯窗内背景介质是看成均匀的， 但在整个模型内背景介质可以是非均匀的. 所以本文中的散射层析成像公式可以适应非均匀的背景介质；其次，前者利用了滤波算子对成像的波场进行滤波校正，但只利用了波场的全局性特征，而本文中则是在局部波数域(或局部角度域)利用了分辨率矩阵对模型中每个点的局部成像矩阵进行滤波校正，从而达到对成像波场的 真振幅校 正；最后，前者利用了SRP (Surface Reflection Profiling)和VSP (Vertical Seismic Profiling) 的数据对谱域的信息进行平均，来减小谱域的盲区. 但没有利用多频的SRP数据对谱域的信息进行平均，而本文中利用了多频的SRP地震数据对模型中的每个点的谱域信息进行平均，来减小谱域的盲区.
与传统的方法相比，本方法中最耗时的要数对模型中每个点处的波场进行局部成像矩阵分解，并将分解得到的成像矩阵作为波场的传播系数进行传播，这就导致了计算量的增加，且增加的计算量基本上与分解得到的局部成像矩阵的大小成正比.
广义散射层析成像反演是一种非均匀介质中基于局部Born模型的体散射的层析成像反演，且反演过程中的入射波场和格林函数都是利用前向散射的波场重新归一化的.通过对两种相对于背景速度不同的平均速度扰动(23%和30%)模型的Born数据和有限差分数据的测试结果可知，在对多个不同频率的局部成像矩阵的谱进行叠加平均之后谱域的盲区的信息大部分都得到恢复，校正后谱的质量无论从覆盖的面积范围和值的均一性上都有了很大的改善；在本文的模型的测验中，该反演方法既可以很好地恢复方块异常体速度模型中相对于背景速度的平均速度扰动达到23%的速度异常体的低波数分量，也可对于不同的背景速度模型都能很好地恢复Marmousi速度模型中的低频分量.所以广义散射层析成像反演能适应模型参数的体扰动较强的情况，能很好对局部成像矩阵的谱进行有限频率带宽和空间孔径的校正，且它能克服一些经典的Born模型的限制条件和不足之处，扩展了基于Born模型的适用范围.
作为后续进一步的研究可以从两个方面来着手.在利用非均匀广义散射层析成像反演来恢复速度模型的低波数分量的方法中，可以测试参数扰动更大的介质模型、不同的背景速度模型和定义一些影响速度模型长波长分量的恢复的因子来进一步提高广义散射层析成像反演的适用范围和更好地对速度模型进行恢复；利用基于Kirchhoff积分的边界散射的层析成像反演来恢复速度模型的高频分量(即模型的边界)，结合基于体散射的层析成像反演的速度模型的低频分量，本方法应该能够很好地重建速度模型.
致 谢 在本文的研究中得到了谢小碧教授、刘伊克研究员、宁杰远教授、曹军、郑应才、何耀峰、杨辉、贾晓峰副教授、毛剑、闫蕊、吴邦玉、耿瑜、任浩然和徐文君等人的帮助，在此表示诚挚的感谢！| [1] | Back us G E, Gilbert F. 1967. Numerical applications of a formalism for geophysical inverse problems. Geophys. J. Royal Astron. Soc. , 13(1-3): 247-246. |
| [2] | Backus G E, Gilbert F. 1968. The resolving power of gross earth data. Geophys. J. Royal Astron. Soc. , 16(2): 169-225. |
| [3] | Backus G E, Gilbert F. 1970. Uniqueness in the inversion of inaccurate gross earth data. Phil. Trans. Roy. Soc. London Ser. A, 266(1173): 123-192. |
| [4] | Beylkin G. 1985. Imaging of discontinuities in the inverse scattering problem by inversion of a causal generalized Radon transform. J. Math. Phys. , 26(1): 99-108. |
| [5] | Beylkin G, Otistaglio M, Miller D. 1986. Spatial resolution of migration algorithms. Acoustic Imaging, 14: 155-168. |
| [6] | Candès E J, Donoho D. 2000. Curvelets: a surprisingly effective nonadaptive representation for objects with edges: Curves and surfaces fitting. Cohen A, Rabut C, Schumarker L eds. Venderbilt University Press, Nashville, 105-120. |
| [7] | Candès E J, Demanet L. 2005. The curvelet representation of wave propagators is optimally sparse. Comm. Pure Appl. Math. , 58-11: 1472-1528. |
| [8] | Cao J. 2008. Toward a wave-equation based true-reflection imaging [Ph.D. thesis]. Santa Cruz: University of California. |
| [9] | Chang X, Liu Y K. 1998. Distinguishing seismic first break by means of Hausdorff fractal dimension. Chinese J. Geophys. (in Chinese), 41(6): 826-832. |
| [10] | Chavent G, Plessix R E. 1999. An optimal true-amplitude least squares prestack depth-migration operator. Geophysics, 64(2): 508-515. |
| [11] | Chen L, Wu R S, Chen Y. 2006. Target-oriented beamlet migration based on Gabor-Daubechies frame decomposition. Geophysics, 71(2): S37-S52. |
| [12] | Claerbout J F. 1971. Toward a unified theory of reflector mapping. Geophysics, 36(3): 467-481. |
| [13] | Claerbout J F. 1985. Imaging the Earth Interior. Oxford: Blackwell Scientific Publications. |
| [14] | Clapp M, Clapp R, Biondi B. 2005. Regularized least-squares inversion for 3-D subsalt imaging. 75th Annual International Meeting, SEG, Expanded Abstracts, 1814-1817. |
| [15] | De Wolf D A. 1971. Electromagnetic reflection from an extended turbulent medium. cumulative forward-scatter single-backscatter approximation. IEEE Trans. Ant. Propag. , AP-19(2): 254-262. |
| [16] | De Wolf D A. 1985. Renormalization of EM fields in application to large-angle scattering from randomly continuous media and sparse particle distributions. IEEE Trans. Ant. Propag. , AP-33(10): 608-615. |
| [17] | Devaney A J. 1982. A filtered back propagation algorithm for diffraction tomography. Ultrasonic Imaging, 4(4): 336-350. |
| [18] | Devaney A J. 1984. Geophysical diffraction tomography. Geoscience and Remote Sensing, IEEE Transactions on, GE-22(1): 3-13. |
| [19] | Gelius L J, Johansen I, Sponheim N, et al. 1991. A generalized diffraction tomography algorithm. J. Acoust. Soc. Am. , 89(2): 523-528. |
| [20] | Hu J, Schuster G T, Valasek P. 2001. Poststack migration deconvolution. Geophysics, 66(3): 939-952. |
| [21] | Kuhl H, Sacchi D. 2003. Least-squares wave-equation migration for AVP/AVA inversion. Geophysics, 68(1): 262-273. |
| [22] | Lailly P. 1983. The seismic inverse problem as a sequence of before stack Migrations. SIAM Conference on Inverse Scattering: Theory and Applications. |
| [23] | Liu Y K, Chang X. 2000. Quantitative assessment of inversion solution of seismic tomographys and its application. Chinese J. Geophys. (in Chinese), 43(2): 251-256. |
| [24] | Liu Y K, Sun H C, Chang X. 2005. Least-squares wave-path migration. Geophys. Prosp. , 53(6): 811-816. |
| [25] | Liu Y K, Sun H C, Chang X. 2007. Crosswell imaging by two-dimensional oriented wave path migration. Geophys. J. Int. , 169(1): 233-239. |
| [26] | Pratt R G. 1999. Seismic waveform inversion in the frequency domain, Part 1: Theory and verification in a physical scale model. Geophysics, 64(3): 888-901. |
| [27] | Pratt R G, Shin C S, Hicks G J. 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion. Geophys. J. Int. , 133(2): 341-262. |
| [28] | Pratt R G, Worthington M H. 1988. The application of diffraction tomography to cross-hole seismic data. Geophysics, 53(10): 1284-1294. |
| [29] | Pratt R G, Worthington M H. 1990. Inverse theory applied to multi-source cross-hole tomography I: Acoustic wave equation method. Geophys. Prosp. , 38(3): 287-310. |
| [30] | Rickett J E. 2003. Illumination-based normalization for wave-equation depth migration. Geophysics, 68(4): 1371-1379. |
| [31] | Schlottmann R B. 2006. Direct waveform inversion via iterative inverse propagation. 76th Annual International Meeting, SEG, Expanded Abstracts, 2146-2150. |
| [32] | Schultz K I, Jaggard D L. 1987. Microwave projection imaging for refractive objects: a new method. Journal of Optical Society of America, A4: 1773-1782. |
| [33] | Sirgue L, Pratt R G. 2004. Efficient waveform inversion and imaging: A strategy for selecting temporal frequencies. Geophysics, 69: 231-248. |
| [34] | Slaney M, Kak A C, Larsen L. 1984. Limitations of imaging with first-order diffraction tomography. IEEE Trans. Microwave Theory and Techniques, MTT-32(8): 860-873. |
| [35] | Tarantola A. 1984a. Inversion of seismic reflection data in the acoustic approximation. Geophysics, 49(8): 1259-1266. |
| [36] | Tarantola A. 1984b. Linearized inversion of seismic reflection data. Geophys. Prosp. , 32(6): 998-1015. |
| [37] | Tarantola A. 2005. Inverse problem theory and methods for model parameter estimation. SIAM (first edition, 1987). |
| [38] | Taylor J R. 1972. Scattering Theory: The Quantum Theory of Nonrelativistic Collisions. Wiley. |
| [39] | ten Kroode A P E, Smit D J, Verdel A R. 1994. A microlocal analysis of migration. Wave Motion, 28(2): 149-172. |
| [40] | Valenciano A, Biondi B, Guitton A. 2006. Target-oriented wave-equation inversion. Geophysics, 71(4): A35-A38. |
| [41] | Woodward M J. 1992. Wave-equation tomography. Geophysics, 57(1): 15-26. |
| [42] | Woodward M J, Nichols D, Zdraveva O, et al. 2008. A decade of tomography. Geophysics, 73: VE5-VE11. |
| [43] | Wu R S. 1996. Synthetic seismograms in heterogeneous media by one-return approximation. Pure and Applied Geophysics, 148(1-2): 155-173. |
| [44] | Wu R S. 2003. Wave propagation, scattering and imaging using Dual-domain oneway and One-return propagators. Pure and Applied Geophysics, 160(3): 509-539. |
| [45] | Wu R S. 2007. Generalized diffraction tomography in inhomogeneous background with finite data aperture. 77th Annual International Meeting, SEG, Expanded Abstracts, 2728-2732. |
| [46] | Wu R S, Cao J. 2005. WKBJ solution and transparent propagators. 67th Annual International Meeting, EAGE, Expanded Abstracts, 167-170. |
| [47] | Wu R S, Chen L. 2001. Beamlet migration using Gabor-Daubechies frame propagator. 63rd Annual International Meeting, EAGE, Expanded Abstracts, 74-78. |
| [48] | Wu R S, Chen L. 2002. Mapping directional illumination and acquisition-aperture efficacy by beamlet propagators. 72nd Annual International Meeting, SEG, Expanded Abstracts, 1352-1355. |
| [49] | Wu R S, Chen L. 2006. Directional illumination analysis using beamlet decomposition and propagation. Geophysics, 71(4): S147-S159. |
| [50] | Wu R S, Luo M, Chen S, Xie X B. 2004. Acquisition aperture correction in angle-domain and true-amplitude imaging for wave equation migration. 74th Annual International Meeting, SEG, Expanded Abstracts, 937-940. |
| [51] | Wu R S, Mao J. 2007. Beamlet migration using local harmonic bases. 77th Annual International Meeting, SEG, Expanded Abstracts, 2230-2234. |
| [52] | Wu R S, Toksöz M N. 1987. Diffraction tomography and multisource holography applied to seismic imaging. Geophysics, 52(1): 11-25. |
| [53] | Wu R S, Wang Y Z, Luo M Q. 2008. Beamlet migration using local cosine basis. Geophysics, 73(5): S207-S217. |
| [54] | Wu R S, Xie X B. 2009. Modeling primaries of acoustic/elastic waves by one-return approximation. The Leading Edge, 28(5): 576-581. |
| [55] | Wu R S, Xie X B, Wu X Y. 2006. One-way and one return approximations (De Wolf approximation) for fast elastic wave modeling in complex media. Advance in Geophysics, 48: 265-322. |
| [56] | Wu R S, Xie X B, Wu X Y. 2007. One-way and one-return approximations (De Wolf approximation) for fast elastic wave modeling in complex media.// Wu R S, Maupin V eds. Advances in Geophysics, 48: 265-322. |
| [57] | Xie X B, Jin S W, Wu R S. 2003. Three-dimensional illumination analysis using wave-equation based propagator. 73rd Annual International Meeting, SEG, Expanded Abstracts, 989-992. |
| [58] | Xie X B, Jin S W, Wu R S. 2006. Wave equation based seismic illumination analysis. Geophysics, 71(5): S169-S177. |
| [59] | Xie X B, Wu R S. 2002. Extracting angle domain information from migrated wave field. 72th Annual International Meeting, SEG, Expanded Abstracts, 21(2):1360-1363. doi:10.1190/1.1816910. |
| [60] | Xie X B, Wu R S, Fehler M, Huang L. 2005. Seismic resolution and illumination: A wave-equation-based analysis. 75th Annual International Meeting, SEG, Expanded Abstracts, 1862-1865. |
| [61] | Xu S, Lambaré G. 2006. True amplitude Kirchhoff prestack depth migration in complex media. Chinese J. Geophys.(in Chinese), 49(5): 1431-1444. |
| [62] |
Xu S, Lambaré G. 2009. True amplitude prestack depth migration in phase space. Geophysical Prospecting for Petroleum (in Chinese), 48(1): 34-39. |
| [63] | Yang W C. 2002. A perspective to development of geophysical inversions. Earth Science Frontiers (in Chinese), 9(4): 389-396. |
| [64] | Yang W C, Du J Y. 1994. A new algorithm of seismic tomography with application to engineering detections. Chinese J. Geophys. (in Chinese), 37(2): 239-244. |
| [65] | Yang W C, Li Y M. 1993. Applied Seismic Tomography (in Chinese). Beijing: Geological Publishing House. |
| [66] | Yu Z, MeMechan G A, Anno P D, et al. 2004. Wavelet transform-based prestack multiscale Kirchhoff migration. Geophysics, 69: 1505-1512. |
| [67] | Zhu X S. 2010. Diffraction points imaging and generalized diffraction tomography [Ph.D. thesis] (in Chinese). Beijing: Peking University. |
| [68] | Zhu X S, Wu R S. 2009. Numerical tests on generalized diffraction tomography. 79th Annual International Meeting, SEG, Expanded Abstracts, 2307-2311. |
| [69] | Zhu X S, Wu R S. 2010. Imaging diffraction points using the local image matrices generated in prestack migration. Geophysics, 75(1): S1-S9. |
| [70] | 常旭 ,刘伊克. 1998. Hausdorff分数维识别地震道初至走时. 地球物理学报, 41(6): 826-832. |
| [71] | 刘伊克, 常旭. 2000. 地震层析成像反演中解的定量评价及其应用. 地球物理学报, 43(2): 251-256. |
| [72] | 杨文采, 李幼铭. 1993. 应用地震层析成像. 北京: 地质出版社. |
| [73] | 杨文采, 杜剑渊. 1994. 层析成像新算法及其在工程检测上的应用. 地球物理学报, 37(2): 239-244. |
| [74] | 杨文采. 2002. 评地球物理反演的发展趋向. 地学前缘, 9(4): 389-396. |
| [75] | 徐升, Gilles Lambaré. 2009. 相空间中的保真振幅偏移. 石油物探,48(1): 34-39. |
| [76] | 徐升,Gilles Lambaré. 2006. 复杂介质下保真振幅Kirchhoff深度偏移. 地球物理学报, 49(5):1431-1444. |
| [77] | 朱小三. 2010. 地震散射点成像和非均匀介质中广义散射层析成像反演[博士论文]. 北京: 北京大学. |