地震散射波的高精度数值模拟与振幅分析

刘铁华

中铁第四勘察设计院集团有限公司，武汉 430063

地震散射波的高精度数值模拟与振幅分析

刘铁华

中铁第四勘察设计院集团有限公司，武汉 430063

随着地震勘探要求和技术的不断提高，散射波地震勘探方法逐渐引起人们的重视，尤其在构造复杂地区的高精度勘探中，散射波理论具有广泛的应用前景．本文在讨论了现阶段几种散射波数值模拟方法的利弊后，设计了一种基于微扰论的FK域积分法，在散射场的二次震源和空间能量衰减处理两方面做了改进，其计算精度和效率均得到了提高．此外，采用该算法针对散射波动力学进行了正演模拟试算，从扰动介质空间分布和扰动量两个角度定量地分析了散射波振幅动力学方面的特征．在此基础上得到了五个经验常数，并给出了对应的计算公式和物理意义，通过这五个常数可以划分散射波不同形态的动力学特征．

地震散射波，高精度，数值模拟，FK域积分，微扰论，动力学特征，定量分析

1 引 言

地震勘探已逐步从地震地质条件相对优越的地区向复杂地区推移，人们对地震勘探的精度提出了更高的要求，基于水平层状介质理论的反射法地震勘探技术遇到了前所未有的挑战，在介质适应能力和信息提取能力等方面均难以达到勘探精度的要求［1］．

事实上对于石油、煤田和工程类而言，水平层状介质的假设并不完全适用，非均匀介质理论更具普适性，可以弥补基于水平介质假设的反射法地震勘探的这一缺陷．散射波是由地下非均匀体扰动引起的波场畸变现象，一般情况下，把能够用几何光学（射线理论）近似处理的情况剔除在外，这种散射波被称为狭义的地震散射波［2］．传统上所讨论的各种波动形式（如反射波、折射波、透射波等）均是散射波在不同约束条件下的表现，因此对散射波理论的研究具有深远意义，近十年来，散射波理论在基础理论和正反演技术思路探索方面得到了长足的发展［3－7］．

早期，在能量衰减机制方面［1，3］涉及到散射波动力学方面的研究，在后期，散射波研究主要集中在运动学理论和应用方面［8－15］，关于宏观散射波动力学研究则比较少．散射波理论中包括单次散射、多次散射等诸多重要的领域［1，2，4］，本文从单次散射波理论入手，重点讨论其振幅这一动力学因素的特征，定量地分析它与扰动介质之间的关系．

2 散射场计算方法

在散射波数值模拟方面，1995年Wu和Huang［8］提出了相位屏算子，计算了二维情况下的散射波场；1998年，符力耘等［15］采用配置法来对体积分方程进行了求解，将非均匀介质等效地看作为某一均匀背景介质扰动情况下的扰动；同年，David W．Eaton采用Born近似和射线理论近似的方法进行了散射场的计算；2001年，孙明［16］采用高斯射线束的方法进行了简单块状构造的模拟，取得了一定的成果；2003年，黄雪继［17］基于微扰论，采用相移法求解波动方程计算一次场，再用FK域积分公式进行了散射场的正演模拟；2007年，秦雪霏［18］采用六阶有限差分进行了散射波的数值模拟；2009年，郭继亮［19］采用有限差分法针对随机介质进行了模拟并对波场进行了分析．

2．1 数值模拟计算方法

在基于微扰论的现有各种散射波数值模拟方法中，几种是基于常规波动方程方法的，这些方法的数学模型还不够完备，只能从某些方面反映散射场的形态．有些方法有了严密散射理论的方法（Wu和黄雪继），数学模型比较完备，他们都是基于微扰论的，均将扰动对入射波的响应当成二次的信号源（二次震源）．基本思想如下：

根据微扰理论，可将介质分为背景介质和扰动介质两部分，介质Q可以描述为背景介质Q0（x）和扰动介质δQ（x）：

基于声波方程，对应的波场可写为

其中Xs表示炮点位置，Xr表示检波点位置，U0表示背景介质中从炮点直接传波到检波点的波场，也称为一次场，Us表示经过地下介质扰动后形成的散射场．

对（2）式采用格林函数法可以得到散射波的计算公式［8］：

其中

n为介质扰动量，G为均匀介质中的格林函数，Y可以看作是求取散射场的二次震源场．微扰论认为，在扰动量微弱的前提下一次场等于总场，散射场只是个微小量．式（4）中的U＊本质意义是总场，在微扰论范围类可以用一次场U0代替．对公式（3）在横向做傅氏变换，采用频率波数域的格林函数［8］代入（3）式得到下式：

其中

Wu和黄雪继在处理二次震源Y时采用了不同的计算方法，且都存在着二次震源近似和误差累计的情况．首先，前人将二次震源表达式（4）中的U＊用背景场U0代替，在背景场微小的情况下是可以的，但是，扰动量增大后就存在较大的误差；其次，在计算一次场U0时采用的都是FK域的方法，存在与计算散射场时同样的FK域数值误差．散射场本身就是一种扰动场，如果误差累计增大就必然影响波场的识别；最后，在计算一次场时都是基于单程波方程，由于计算散射场前要得到每个深度的全时间数据，而相位屏法或相移法一次只能计算一个深度处的记录，虽然在计算一个深度的时间花费不高，但采样加密深度增大后就需要很长的时间，所以计算效率不是很理想．

针对以上三点，本文设计了一种基于总场的FK域散射场计算方法，将（4）式中的U＊不做近似，采用常规高阶有限差分法计算，采用了两个修正项消除FK的数值干扰和衰减问题保证了能量的合理性．进行算法改进后，一方面，理论上散射场的源U＊不存在误差，而且不存在FK域的误差累计现象；另一方面，总场不需要额外计算，且总场中不包含FK域积分计算的数值误差，便于波场的比较分析．在二维情况下（5）式改写为［20］

其中

其中Nz为纵向网格个数，（7）式为修正后的散射波计算公式，（8）式是计算散射的源场（二次场）计算公式，（11）、（12）和（13）是用于进行频散和能量处理的修正公式．这里给出的表达形式多为形式上的解析形式，因入射波场由其他数值方法解得，所以这里的解为解析—数值混合形式，该方法的精度与计算入射总场U时所使用的数值方法的精度相关．

在程序实现过程中，采用了两步法进行计算，先采用常规波动方程正演方法（程序中采用的是高阶有限差分法）求取波场U，该波场是包含了一次场和散射场信息的全波场U，本文设计的计算方法中以全波场U代替前人采用的一次场U0，然后由公式（7）和（8）计算散射场Us．从散射场的计算过程可以看出，散射场Us是常规地震波场U的再次提取，是地震波场与介质异常耦合后的积分结果．在计算U0过程中采用均匀背景介质中的格林函数G0，存在一定的近似，在背景介质变化较大（超过15%）时误差将会比较明显，所以这种方法在背景介质变化率小于15%时较为理想．

2．2 模型试算

为了体现算法的效果，现结合两个模型进行计算测试，分别是点模型和一个复合模型．

模型一 如图1所示，背景介质速度为2500m／s，三个深度分别为300m、400m和500m的扰动点速度为2800m／s，震源位于扰动点正上方地表处，震源子波主频为40Hz，检波器均匀分布于地表．

图1 点模型Fig．1 Model of point

经过试算得到如图2所示的单炮散射波记录，三个双曲线同相轴代表着三个扰动点的点散射，不论从走时还是组合关系来看都比较合理．与常规方法相比有个明显的不同点是：没有直达波，这是因为散射是扰动场的反映，直达波没有经过扰动改造，所以在散射场中体现不出来．

图2 散射波单炮记录Fig．2 Scatter record of single shot

图3 复合模型Fig．3 Model of complex

模型二 如图3所示，该模型除了考虑盐丘构造还一并分析了向斜构造和小角度斜层构造．从上往下共有三个主要分界面，在第一层水平界面上有一强扰动点，该扰动点坐标为（700m，200m），速度为3000m／s，一个顶部带有背斜构造的盐丘体穿透下面的两个地层．震源位（1000m，20m）处，震源子波主频为40Hz，接收点均匀分布于地表．

图4是通过本文的改进算法计算出的散射场单炮记录，可以发现：首先，散射场中是没有直达波信息的；其次，在散射场中出现了由于盐丘顶部的背斜构造形成的“蝴蝶结”；再次，可以清楚地发现由于界面1上扰动点形成的点散射．从整个记录上可以发现，波场的细节部分比较明显．

散射场是由于介质的扰动对波场的改造形成的波场，所以散射能量与介质扰动量的分布有着密切的关系．如果将介质中每个点所有时刻中最大的散射振幅抽取出来就可以得到散射最大振幅剖面，图5就是复合模型对应的散射最大振幅剖面．

散射最大振幅剖面（图5）中的虚线是盐丘的实际轮廓，与散射能量的分布在顶部存在较好的对应关系，界面2和界面3也吻合较好．在盐丘深部，由于照明不充分轮廓没有很好地反映出来．此外，由于模型底部的扰动量相对集中且较高，从而界面1在图5中反映的不明显．从散射最大振幅剖面与模型的吻合度可以发现，散射理论的介质信息提取能比较突出，再加上对细节的描述能力较强，从而分辨率也是比较高的，这也从另一方面证实了本算法的正确性．

图4 散射场单炮记录Fig．4 Scatter record of single shot

图5 散射最大振幅剖面Fig．5 Scattered wave maximum amplitude profile

3 散射波振幅特征

由（3）式可以看出，散射波的能量是由扰动区域内各点处的一次场延拓累积而成，只与扰动介质的尺寸和扰动程度有关，与观测系统没有本质的关系．所以本文采用分步单次散射波计算法［13］，重点围绕扰动体的尺寸和扰动量两个因素进行分析．

根据Wu［14］和田丽花［15］的相关文献，散射波场是在FK域计算的，散射场横向的动力学特征相对于纵向更加丰富，所以设计了如图6所示的模型，500个接受点均匀设置于地表，250m深度上存在一个厚为h宽L位于中央的扰动体，介质背景速度为2500m／s，震源位于扰动体的正上方地表处，震源主频40Hz，峰值为1个单位．

3．1 散射振幅与横向尺度的特征

图6 扰动模型Fig．6 Model of perturbing body

在所有接受点中，位于扰动体正上方的第250道因为来自地下介质的能量最为均匀，其动力学特征最具代表性、对波场的变化最敏感、接收到的能量最大，所以在分析时就以该道数据为主要研究对象．下图给出了该道随扰动体横向尺寸变化时的能量曲线（虚线）以及其能量关于扰动体横向尺寸L的梯度关系图（实线）．其中，扰动体的纵向厚度h为2m（远小于子波的波长），扰动体的扰动量为10%．不过从分析的合理性来考虑，讨论尺度L—波数k0变化时的能量特性更为合理，这在Wu等［8，14］前人的工作中有所体现，本文还是想单独地讨论空间异常随空间尺度变化时的特征情况．

图7 尺度L变化时的能量特性Fig．7 The energy characteristics changes over perturbation scale

由图7可以看出，随着扰动体横向尺度L的增大，能量是先快速增加，后经过少许降低再略为抬高后趋于稳定．在梯度曲线上可以看到三个关键分界点，分别是：

L1：首次振幅递增速度拐点．结合实际合成地震记录可以发现，扰动体横向尺度L小于该值时，来自地下各点的散射能量相互加强；在L大于L1后出现相干衰弱现象，此时由断点处形成的点散射能量逐渐凸显出来．总体而言，L1是狭义散射波和绕射波的分界点，其大小等于对应的第一菲涅尔半径；

L2：首次振幅极大值点，随着各扰动点散射波之间干涉作用的逐渐加强，在L达到L2时能量不再增加，达到最大值．超过L2后，内部的各扰动点的散射干涉逐渐演变为反射波，这种现象会随着尺度的变大而更加明显．L2是反射形态波场出现的起始点，其间端点绕射波始终存在；

L3：首次振幅阶段性极小点．尽管由于扰动体横向尺度的继续增加，反射波的能量逐渐增加，但是当L小于L3时，这种能量的增加量没有由于断点绕射远离导致的能量降低量大，所以出现降低的现象；当L达到L3后，反射波的能量逐渐占据绝对比重，能量表现为增加的趋势．所以说，L3是绕射和反射主导地位分界点，当扰动体横向尺度大于L3后记录中主要表现为反射波．

到此，讨论了单道振幅能量的情况，如果将所有接受道振幅的最大和最小能量统计后就出现如图8所示的结果，可以发现L1和L3对应的最大振幅基本一致，L2处出现了最大和最小振幅．如果我们将最终的最大振幅值称为稳定最大振幅，那么L1就是首次达到平衡最大振幅的L值，L3是最终趋于平衡最大振幅的起始点．经过模型试验计算后统计发现L1、L2和L3之间存在着密切的关系，这里给出通过本文算法计算得出的经验公式：

其中R为第一菲涅尔半径，e自然对数的底数，取2．71828183．

3．2 散射动力学特征与纵向尺度的关系

保持观测系统不变，扰动体横向尺寸L固定为2m，扰动体的扰动量为10%，通过改变纵向的尺寸h，统计各道的极值振幅如图9所示．

图8 横向尺度L与振幅极值的关系Fig．8 The relationship between horizontal scale and amplitude extreme

图9 纵向尺度h与振幅极值的关系Fig．9 The relationship between vertical scale and amplitude extreme

由纵向尺度h与振幅极值的关系可以看出，振幅正负最大值的表现特征不太一样，这是由于负极值受到上下两个界面的干涉作用引起的，上界面的负相位和下界面的正相位发生了叠加现象．在图9中有两个关键的分界点，如图中的L4和L5所示：

L4：最高纵向散射分辨精度．该值表示只要扰动介质纵向的尺度小于L4都可以看作是一个散射单元，当L＝L4时正负最大振幅值是一致的．

L5：上下界面波形完全分离的最小尺度．通过实际合成资料的波形分析发现，当L介于L4和L5之间时，上下界面的波形能量是连续的，没有完全分离；当L大于L5后上下波形完全分离，能够直观地定量地划定界面的分界点．

通过试验数据的统计分析发现，由纵向尺度h的散射动力学特征得到的两个关键参数L4和L5的大小可以表示为

其中R为第一菲涅尔半径．

3．3 散射动力学特征与扰动量的关系

维持原有的观测系统不变，扰动体的横纵尺寸均为2m，通过改变扰动量n得到了如图10所示的能量特征曲线，其中实线表示第250道记录的能量走势，虚线反映的是各道平均能量的走势情况．

图10 扰动量变化时的能量特性Fig．10 The energy characteristics changes over perturbation quantity

图10中两条曲线基本一致，记录的能量都是随着扰动量的增加而不断地提高，并且前期的增速明显高于后期，这也是公式（3）的数值表现形式．

4 结 论

本文设计了一种基于总场的散射波的FK域数值模拟方法，该方法具有误差小且计算效率高的特点．由于误差较小，对于介质扰动的信息提取精度得到了提高，由散射场提取出来的最大散射振幅剖面与模型剖面具有较高的吻合度，从而达到了地下扰动信息的保真提取．此外，针对于单次散射波进行的散射波振幅动力学的特征分析，主要讨论了扰动体的空间尺度和扰动量这两个因素对散射波振幅能量的定量关系．通过模型试算得到了五个经验参数，分别是L1、L2、L3、L4和L5，并给出了它们的经验公式和对应的物理意义．

散射波的动力学特征除了振幅信息以外还应包含相位、频率等属性本文没有涉及，此外，在扰动量较大或者扰动介质过大时应该采用多次散射波的计算公式来取代单次散射波，这样的计算精度才能达到动力学分析的要求，这两方面是进一步研究的方向．

致 谢在本文的数据计算过程中得到长安大学朱光明教授和刘永华副教授的点拨，在此表示诚挚的感谢．

（References）

［1］ 尹军杰，刘学伟，李文慧．地震波散射理论及应用研究综述．地球物理学进展，2005，20（1）：123－134．Yin J J，Liu X W，Li W H．The view of seismic wave scattering theory and its applications．Progress in Geophysics（in Chinese），2005，20（1）：123－134．

［2］ 吴如山，［美］安艺敬一．地震波的散射与衰减．李裕澈，卢寿德等译．北京：地震出版社，1993．Wu R S，Aki K．Scattering and Attenuation of Seismic Waves（in Chinese）．Translated by Li Y C，Lu S D，et al．Beijing：Seismological Press，1993．

［3］ Hudson J A，Knopoff L．Predicting the overall properties of composite materials with small－scale inclusions or cracks．Pure and Applied Geophysics，1989，131（4）：551－576．

［4］ Tournat V，Pagneux V，Lafarge D，et al．Multiple scattering of acoustic waves and porous absorbing media．Physical Review E，2004，70（2）：1－10．

［5］ van Wijk K，Komatitsch D，Scales J A，et al．Analysis of strong scattering at the micro－scale．J．Acoust．Soc．Am．，2004，115（3）：1006－1011．

［6］ 贾豫葛，李小凡，张美根等．地震波散射研究的若干重要进展．地球物理学进展，2005，20（4）：939－944．Jia Y G，Li X F，Zhang M G，et al．Some important progress in research upon seismic wave scattering．Progress in Geophysics（in Chinese），2005，20（4）：939－944．

［7］ 高伟，王宁．基于背向散射脉冲传播时间统计特性的随机界面参数反演研究．地球物理学报，2008，51（1）：260－265．Gao W，Wang L．Inversion of rough surface parameters based on travel－time statistics of backscattered pulses．Chinese J．Geophys．（in Chinese），2008，51（1）：260－265．

［8］ Wu R S，Huang L J，Xie X B．Backscattered wave calculation using the De wolf approximation and a phase－screen propagator：Expanded Abstracts of the Technical Program．SEG′65th Annual Meeting，1995：1293－1296．

［9］ 李振春，王清振．地震波衰减机理及能量补偿研究综述．地球物理学进展，2007，22（4）：1147－1152．Li Z C，Wang Q Z．A review of research on mechanism of seismic attenuation and energy compensation．Progress in Geophysics（in Chinese），2007，22（4）：1147－1152．

［10］ 刘洪，袁江华，勾永峰等．地震逆散射波场和算子的谱分解．地球物理学报，2007，50（1）：240－247．Liu H，Yuan J H，Gou Y F，et al．Spectral factorization of wavefield and operator in seismic inverse scattering．Chinese J．Geophys．（in Chinese），2007，50（1）：240－247．

［11］ 陶智勇．弹性波逆散射微扰论．地球物理学报，2005，48（4）：918－923．Tao Z Y．Elastic wave inverse scattering perturbation theory．Chinese J．Geophys．（in Chinese），2005，48（4）：918－923．

［12］ 丁科，宋守根，谢忠球．逆散射理论的发展及应用前景．地球物理学进展，2005，20（3）：661－666．Ding K，Song S G，Xie Z Q．Development and future application of inverse scattering theory．Progress in Geophys．（in Chinese），2005，20（3）：661－666．

［13］ 高伟，王宁．基于背向散射脉冲传播时间统计特性的随机界面参数反演研究．地球物理学报，2008，51（1）：260－265．Gao W，Wang N．Inversion of rough surface parameters based on travel－time statistics of backscattered pulses．Chinese J．Geophys．（in Chinese），2008，51（1）：260－265．

［14］ Wu R S，Aki K．Elastic wave scattering by a random medium and the small scale inhomogeneities in the lithosphere．Journal of Geophysical Research，1985，90（B12）：10261－10273．

［15］ 符力耘，杨慧珠，牟永光．非均匀介质散射问题的体积分方程数值解法．力学学报，1998，30（4）：488－494．Fu L Y，Yang H J，Mu Y G．Volume integral equation method for scattering problems of inhomogeneous media．Journal of Theoretical and Applied Mechanics（in Chinese），1998，30（4）：488－494．

［16］ 孙明，林君．金属矿地震散射波场的数值模拟研究．地质与勘探，2001，37（4）：68－70．Sun M，Lin J．Study of seismic scattering wave field numerical modelfor metallic ore exploration．Geology and Prospecting（in Chinese），2001，37（4）：68－70．

［17］ 田丽花．地震波散射场正演数值模拟．北京：中国地质大学（北京），2007．Tian L H．Seismic Wave Numerical Value Forward Modeling of the Scattered Field（in Chinese）．Beijing：China University of Geosciences（Beijing），2007．

［18］ 秦雪霏．地震波多次散射波场属性分析．长春：吉林大学，2007．Qin X F．The Characters Analysis about the Multiply Scatter of Seismic Wave（in Chinese）．Changchun：Jilin University，2007．

［19］ 郭继亮．随机介质散射波场数值模拟与特征分析．西安：长安大学，2009．Guo J L．The Numerical Simulation of Scattering in Random Media and the Analysis of the Wavefield（in Chinese）．Xi′an：Chang′an University，2009．

［20］ 刘铁华．地震散射波理论与数值模拟．西安：长安大学，2010．Liu T H．The Theory and Numerical Simulation of Scattering（in Chinese）．Xi′an：Chang′an University，2010．

High precision numerical modeling and dynamic characteristic analysis of seismic scattering wave

LIU Tie－Hua
China Railway Siyuan Survey and Design Group CO．，LTD，Wuhan 430063，China

The scattered wave seismic method is now well appreciated and developed with the advancement of technology and requirement of seismic survey．Particularly in high－precision prospecting of complicated structure area，it has broad application prospect．Different current numerical modeling of scattered wave has their advantages and disadvantages．This paper has designed a method of Integration in FKdomain based on perturbation theory to solve scattered field，which could obtain scattered field by integration in FKdomain after get the total wave field．It has improved the data source of scattered field and the processing of space energy．The computational accuracy and efficiency of scattered field have got better．In addition，based on the algorithm a forward modeling test is carried out for the dynamics of scattered field．The features of kinetics of scattered wave amplitude are analyzed quantitatively from the perspectives of space distribution of perturbing body and disturbance quantity．Based on this，it obtains five experiential constants and provides the corresponding formula and physical significance．Scattered wave could be classified into different dynamical characteristics by the constants．

Scattered wave，High precision，Numerical modeling，Integration in FKdomain，Perturbation theory，Dynamic characteristic，Quantitative analysis

P631收修定稿2011－11－08，2012－02－27收修定稿

刘铁华．地震散射波的高精度数值模拟与振幅分析．地球物理学报，2012，55（4）：1318－1324，

10．6038／j．issn．0001－5733．2012．04．027．

Liu T H．High precision numerical modeling and dynamic characteristic analysis of seismic scattering wave．Chinese J．Geophys．（in Chinese），2012，55（4）：1318－1324，doi：10．6038／j．issn．0001－5733．2012．04．027．

10．6038／j．issn．0001－5733．2012．04．027

刘铁华，男，1983年生，工学硕士，助理工程师，主要从事散射波理论和工程地球物理综合方法的研究与应用工作．E－mial：lthrwx＠163．com

（本文编辑 胡素芳）