自组织是自然界的普遍现象，反映的是系统如何自动地由无序走向有序.如河流的沉积，首先沉积的是大颗粒，并逐渐沉积小颗粒.地球介质内部，46亿年的演化导致介质分布逐渐从混沌无序，到产生地核、地幔和地壳等几个有序区域，并逐渐从低级有序走向高级有序.随机介质表征的地球内部介质的自组织性，反映的是大尺度背景场上附加的小尺度非均匀性和统计性特征，代表了如构造活动的强烈程度和物质交换程度，有丰富的地球动力学含义.大量的观测和研究表明，无论从空间分布，还是尺度变化上看，地球介质内部的非均匀性是普遍存在的(Aki，1969；Wu，1982；Zhang et al., 1993；Yang et al., 2002；张中杰，2002；张智等， 2005，2013；杨顶辉，2002；滕吉文等，2004；Zhang and Klemperer， 2005；吴国忱等，2010；梁锴等，2011；兰海强等， 2011，2012；刘有山等，2013；Crampin and Gao， 2013；Liu et al., 2014)，小尺度扰动比大尺度背景场对反向散射波场的特征影响更大(Levander and Holliger， 1992；Ikelle et al., 1993；Baig et al., 2003).如测井曲线和岩芯样本显示出的小尺度非均匀性(Wu et al., 1994; 吴何珍等，2008)，地震地壳尺度的CMP叠加剖面上常常存在大量的横向上不连续，范围不大的较短的同相轴的反射段等等(Frankel and Clayton， 1986；Zhang et al., 2004).根据复杂介质中非均匀地质体(或非均匀异常)随机分布的特点，可借助基于统计学方法的随机介质模型进行研究和描述(Ikelle et al., 1993；奚先和姚姚，2005；徐涛等，2007；刘永霞等， 2007a，2007b；Xu et al., 2010，2014；郭乃川等，2012).通过随机介质表征的自组织介质模型上的波场模拟，使我们可能更加充分地利用实际地震记录中的信息，尤其某些以往认为是干扰的信息，并研究介质在小尺度上的非均匀特征.

不同类型、不同尺度的非均匀性会引起不同形式的地震波散射.对于实际的地震波，介质中非均匀异常引起的波前扰动，会随着传播距离的增大而减小，即波前愈合效应.而射线理论是基于对地震波无限高频近似而建立起来的，认为在射线路径上的速度异常体造成的震相初至扰动与震源距无关(Ojo and Mereu， 1986；徐涛等，2004；Xu et al., 2006；赵烽帆等，2014).因此，利用射线理论计算的走时与真实波场有一定的差异，这种差异会随着地震波长和震源距的增大而增大(Hung et al., 2001；Baig et al., 2003；Yang and Hung， 2005；江燕和陈晓非，2011).

对于大尺度的地球介质波场模拟，受计算量和计算时间的制约，通常只能进行较低精度的波场模拟.为了同时满足计算量和计算精度的需求，近年来，国内外学者提出了多尺度数值模拟方法，即在有限的范围内采取高密度的网格分布，获得局部区域更精细的波场特征，在保证全区域精度要求的情况下，能够节省大量的计算时间.目前的主要方法包括多尺度有限体积法(Jemry et al., 2003)、多尺度有限元法(Babuska and Osborn， 1983)、均匀化方法(Hou et al., 1999)、非均匀化多尺度方法(E et al., 2003; Ren and E， 2005)以及小波数值均匀化方法(Dorobantu and Engquist， 1998)等.

为研究非均匀介质中波前愈合效应，本文在地球自组织模型基础上进行了非均匀化多尺度的波场模拟，研究背景速度场与附加自组织介质速度场中的波场走时差，并和由射线方法获得的走时差进行比较，探讨波前愈合效应引起的射线理论适应性条件.

自组织介质通常表征为在大尺度速度背景场上叠加小尺度的随机扰动(Frankel and Clayton， 1986). 在自组织介质中，空间坐标点 x 处的速度可表示为

描述介质自组织特征的参数包括自相关函数、相关长度、均值、方差(或标准差)等(Ikelle et al., 1993).其中，自相关函数及其对应的功率谱密度是体现自组织介质分布特征的重要参数，常见的自相关函数主要有高斯型、指数型、von Karman型(零阶)和Kummer型等(Levander and Holliger， 1992； 刘永霞等，2007a).表 1为两种常见的各向同性自相关函数及其功率谱密度(Klimeš，2002).其中a为相关长度，k_r为波数.自组织模型是通过滤波因子在频率域或者波数域对白噪声频谱进行滤波再进行傅里叶反变换得到的.不同的滤波因子产生不同的自组织模型.

表 1(Table 1)

表 1 各向同性自相关函数及其功率谱密度Table 1 Isotropic autocorrelation functions and associated power spectral density 类型 自相关函数 功率谱密度 高斯型 er2/a2 e-a2kr2/4 指数型 er/a (1+a2kr2)-3/2

表 1 各向同性自相关函数及其功率谱密度Table 1 Isotropic autocorrelation functions and associated power spectral density

综上所述，给定自相关函数等系列参数，就可以进行自组织介质的建模.具体过程(Ikelle et al., 1993；Klimeš，2002；徐涛等，2007；刘永霞等，2007b)如下:(1)生成二维白噪声矩阵序列w(x，z)，并计 算其二维傅里叶变换W(k_x，k_z);(2)选择各向异性自

相似滤波因子 (k_x，k_z)，计算U(k_x，k_z)= (k_x，k_z)×W(k_x，k_z)；(3)对U(k_x，k_z)进行反傅里叶变换，得到空间域自组织模型u(x，z);(4)修正随机扰动u(x，z)的均值和方差; 即可得到特定均值、方差和满足指定特征分布的二维随机序列.图 1为采用相同的白噪声生成的各向同性高斯型和指数型介质.其中，横向和纵向相关长度均为10 km.

图 1

Fig. 1

图 1 相同白噪声生成的两种自组织模型 (a)高斯型；(b)指数型.Fig. 1 Two self-organized media generated from same white noises (a)Gaussian；(b)Exponent.

波场数值模型的网格化通常为均匀的，而非均匀化多尺度方法的核心思想为：首先建立包括宏观尺度和微观尺度的两套网格体系，微观尺度在每一个单元内体现，然后分别在宏观尺度和微观尺度上建立波动方程的差分格式，由微观尺度上的解来计算宏观尺度上的数据，从而有效地实现不同尺度间数据的耦合，并最终通过宏观格式的求解得到原问题的宏观解(E et al., 2007).

非均匀化多尺度方法包括两个主要部分，一是根据微观模型的解在粗网格(感兴趣的网格)是实现宏观递推，二是在稀疏(非均匀)空间域在微观求解初始方程.下面以VTI 介质SH 波波动方程为例，介绍非均匀化多尺度方法的实现过程.

VTI 介质SH 波波动方程为

用粗网格(x_i，z_j)将模型离散，其中i，j=1，2，…N，Δx=x_i+1-x_i，Δz=z_j+1-z_j，并且远大于ε. 这样做是为了形成一个宏观模型，其通量形式方程为

假设在t^k时间，我们有粗网格(x_i，z_j)处方程(6)的数值解，表示为U_ij^k.为了找到t^k+1时刻的解U_ij^k+1，我们分三个步骤进行：微观尺度波动方程求解、宏观通量计算以及宏观解的递推.

(1)微观尺度波动方程求解

根据精细网格的差分格式以及构造的初值条件和边界条件，用有限差分方法求解得到每一个单元内波动方程微观尺度的解.

在每个点(x_i，z_j)周围，定义4个微观尺度的ε单元I_i±1/2，j^ε，I_i，j±1/2^ε(如图 2). 定义d_x-=(Δx-ε)/2，d_x+=(Δx+ε)/2，d_z-=(Δz-ε)/2，d_z+=(Δz+ε)/2. 则围绕每个点(x_i，z_j)的4个ε单元可表示为：

图 2

Fig. 2

图 2 粗网格点(xi，zj)与4个ε单元示意图Fig. 2 Schematic diagram of 4 ε-cells at coarse point(xi，zj)

将在微观网格求解初始方程，即在ε域求解ε尺度.微观网格的定义为

对区域I_i+1/2，j^ε采用如下网格剖分(区域I_i-1/2，j^ε、I_i，j+1/2^ε、I_i，j-1/2^ε与之类似)：

原方程在细网格上的离散格式为

(2)微观尺度单元的宏观通量计算

每个单元的总通量为

(3)用通量计算得到宏观解

将波动方程变换成通量和位移的方程，在此基础上，根据单元的通量利用差分得到下一个时刻宏观网格上的解，实现宏观解的递推.

由所得宏观通量J_{xx i+1/2，j}^m和J_{zz i，j+1/2}^m，对宏观方程(6)进行离散求解，得到tⁿ⁺¹时刻的解U_i，jⁿ⁺¹.由方程(6)可得

因为标准有限差分方法需要用微观模型对整个域离散并且求解方程，所以它会产生大量的方程组，难以数值求解.而有限差分非均匀化多尺度方法，其数值的主要工作在微观模型求解，不过这只在初始域的小子域进行.由于微观单元问题是独立的，可以实现并行求解.

我们选取VTI介质来进行数值模拟试验.模型 参数，γ=0.2，ε=10 m，模型尺寸大小为10×10 km，震源位于模型中心位置，震源子波为主频10 Hz的Ricker子波，时间采样间隔为1 ms.介质速度为

图 3

Fig. 3

图 3 非均匀VTI介质波场快照 (a)粗网格方法；(b)精细网格方法；(c)非均匀化多尺度方法.Fig. 3 Snapshots of wave fields in heterogeneous VTI media (a)Coarse grid method；(b)Fine grid method；(c)Heterogeneous multi-scale method.

我们分别通过射线理论方法和波场模拟的方法，比较背景速度场和附加自组织介质速度场的走时差，来探讨波前愈合效应.

对于任意分布的非均匀介质，射线轨迹可以近似离散为直线段的集合，因此射线轨迹的走时T可以近似表示为

建立高斯型和指数型两种各向同性自组织介质，γ=0.0.模型尺寸大小为168×168 km，震源s位于地表一侧，震源子波为主频3 Hz的Ricker子波，时间采样间隔为6 ms，空间采样间隔dx=dz=150 m，背景速度为6.3 km·s^-1，波长λ=2.1 km.介质的扰动标准差ε=0.01，0.02，0.03，相关长度a=0.75λ，1.5λ，2.25λ，3λ.以震源 s 为中心，在半径L=|r-s|=4λ，8λ，12λ，16λ，20λ，40λ，60λ，80λ处的1/4圆弧上等间隔布设检波器 r，间隔为1°，这样得到一组呈辐射状的检波器阵列.

一般情况下，尾波的波动比波前面的波长长、频率低.这是由于在小尺度情况下，低频波向后散射，高频波向前散射造成的(Morse and Ingard， 1986).图 4为t=10 s时高斯型和指数型自组织介质的波场快照.可以看到当取相同的ε和a/λ时，相对于高斯型介质，指数型介质的波场更加丰富.在指数型介质中，与a=3λ时的情况相比，当a=0.75λ时可以清楚的看到小尺度体持续的向后散射(图 4b1和b3).在左侧的高斯型介质中，波前面的走时波纹随着非均匀程度的增大(ε=0.03)，波场响应以波前面的走时波纹为主，并伴随尾波衍射.而右侧的指数型自组织介质的波前面形态非常相似，并有更明显的不连续的尾波.

图 4

Fig. 4

图 4 高斯型和指数型自组织介质在t=10 s时的波场快照 (a)高斯型介质；(b)指数型介质.Fig. 4 Snapshots in Gaussian and exponential self-organized media at t=10 s (a)Gaussian media;(b)Exponential media.

介质的非均匀程度对波前面震荡和尾波的振幅也有着明显的影响.在高斯型自组织介质中，当a=3λ时，ε=0.03的波场的扰动要比ε=0.01时的显著很多.在所有ε=0.03的模拟中，可以清楚的看到波前面的扰动，其中a=3λ时最为明显，具有大尺度的非均匀性.小尺度的非均匀性在波前面产生了小的扰动，并在衍射的过程中不断的被修复(Hung et al., 2001).

我们利用互相关方法(Robinson and Treitel， 1980)计算存在随机扰动的自组织介质与均匀无扰动介质的相关走时差.对于检波器阵列中的每一个接收点处的自组织介质的地震记录和无扰动均匀介质的地震记录做互相关计算，得到互相关走时差为

传播距离L变化时，地震记录也会相应的变化.图 5为当ε=0.01，a=0.75λ时，取角度为5°时的一条射线上的传播距离L分别为4λ，8λ，12λ，16λ时的高斯型和指数型自组织介质地震记录.随着传播距离的增加，在两种介质中，波形变形的程度会加重，在尾波中会出现低振幅震荡，并且互相关方法计算所得的到时差δT会随着传播距离L的增加而增加.

图 5

Fig. 5

图 5 高斯型(a)和指数型(b)自组织介质在ε=0.01，a=0.75λ时地震记录 其中点划线代表无扰动的背景场波场记录，实线代表附加自组织介质的波场记录.Fig. 5 Synthetic seismograms in Gaussian(a) and exponential(b)media，with a heterogeneity strength ε=0.01 and correlation scalelength a=0.75λ Dash-dotted lines denote homogeneous-media seismograms and solid lines denote heterogeneous-media seismograms.

图 6为高斯型和指数型两种自组织介质，当L=40λ时，选取某一固定接收点，且ε=0.01时，相关长度a变化时(0.75λ，1.5λ，2.25λ，3λ)的地震记录.可以看到，随着a的增大，走时差也逐渐增大.而相关长度a通常代表介质非均匀尺度大小，a越小，非均匀尺度越小.表明当非均匀体尺度较小时，波前面的愈合效果要更好，并且指数型自组织介质的地震记录中尾波震荡更加明显.

图 6

Fig. 6

图 6 高斯型(a)和指数型(b)自组织介质在L=40λ，ε=0.01时地震记录 其中点划线代表无扰动的背景场波场记录，实线代表附加自组织介质的波场记录.Fig. 6 Synthetic seismograms in Gaussian(a) and exponential(b)media，with a heterogeneity strength ε=0.01 and source-receiver distance L=40λ Dash-dotted lines denote the homogeneous-media seismograms and solid lines denote heterogeneous-media seismograms.

为了研究波前愈合效应，比较射线理论和波场模拟的走时差差异，利用非均匀化多尺度方法对不同相关长度和扰动程度(标准差)的高斯型自组织介质和无扰动均匀介质进行波场模拟，采用辐射状的检波器阵列记录不同传播方向和不同震源距的自组织介质和无扰动均匀介质的地震记录.利用互相关方法求取二者的互相关走时差.同时利用射线理论也可以得到高斯型自组织介质和无扰动速度背景场的走时差.将两种方法得到的走时差进行比较时，以射线理论的结果为横坐标，波场模拟的结果作为纵 坐标画出走时差散点图(图 7a和图 8a).同样的，可 以画出指数型自组织介质的走时差散点(图 7b和图 8b).

图 7

Fig. 7

图 7 高斯型(a)和指数型(b)自组织介质互相关走时差与射线理论走时差散点图(ε=0.01，L=40λ) 横坐标为射线理论计算得到的走时差，纵坐标为互相关计算得到的走时差；实线为走时差散点图拟合的斜率，虚线斜率为1.0.Fig. 7 Scatter plots of cross-correlation traveltime shift(vertical axis)with the corresponding ray-theoretical prediction(horizontal axis)in Gaussian(a) and exponential(b)self-organized media with ε=0.01 and L=40λ X-coordinate denotes traveltime shifts derived from ray theory.Z-coordinate denotes traveltime shifts derived from cross-correlation calculation.A least-squares criterion is used to find the best-fitting line(solid); dashed lines denote slope of 1.0.

图 8

Fig. 8

图 8 高斯型(a)和指数型(b)自组织介质互相关走时差与射线理论走时差散点图(ε=0.01，a=0.75λ) 横坐标为射线理论计算得到的走时差，纵坐标为互相关计算得到的走时差；实线为走时差散点图拟合的斜率，虚线斜率为1.0.Fig. 8 Scatter plots of cross-correlation traveltime shift(vertical axis)with the corresponding ray-theoretical prediction(horizontal axis)in Gaussian(a) and exponential(b)self-organized media with ε=0.01 and a=0.75λ. X-coordinate denotes traveltime shifts derived from ray theory. Z-coordinate denotes traveltime shifts derived from cross-correlation calculation. A least-squares criterion hase been used to find the best-fitting line(solid); dashed lines denote slope of 1.0.

由于波场模拟能模拟直达波、反射波、散射波、折射波、转换波、多次波等全波场信息，能够模拟较真实的地震波场特征.而高频近似的射线理论，通常只考虑部分波场信息.在本文中，在速度扰动为一级近似的情况下，炮点与接收点之间用直线连接，并直接计算速度扰动导致的直达波走时差.如果射线理论的结果精度足够高，两种方法得到的走时差就基本相等，则以这两个结果为横纵坐标的点就会沿着过原点且斜率为1.0的直线分布，表明波前愈合效应越不明显.图 7是高斯型和指数型自组织介质当传播距离L=40λ，相关长度a分别等于0.75λ，1.5λ，2.25λ，3λ时的走时差散点图.图 8是高斯型和指数型自组织介质对应当a=0.75λ，L分别等于20λ，40λ，60λ，80λ时的走时差散点图.这里四个图中自组织介质的非均匀扰动标准差都为ε=0.01，实线是对散点做最小二乘拟合的结果.

从图 7和图 8可以看出，当a/λ较小、L/λ较大时，走时差散点图拟合斜率(实线)和1.0(虚线)相差较大，存在更明显的波前愈合效应：波场模拟得到的走时差要比射线理论计算得到的走时差小.这是由于射线理论认为在传播的射线路径上，波前面可以清楚的“记得”每一个小尺度异常，且与传播的射线路径之外的速度扰动场无关.而实际波场传播过程中，由于多种波场的叠加，全空间的小尺度异常扰动均能对接收点的走时产生影响，即产生波前愈合效应.

波场模拟结果与射线理论计算结果的拟合斜率是散射波波前愈合重要性的度量.基于模拟的所有结果，图 9为对于不同的ε时，每一个a/(λL)^1/2对应的斜率值，左边是高斯型自组织介质，右边为指数型自组织介质.从图中可以看到，两种自组织介质的散点拟合斜率值在a/(λL)^1/2≤0.5时变化剧烈，表明波前愈合效应俞明显，射线理论计算的走时误差较大；当a/(λL)^1/2>0.5大于时，拟合的斜率逐渐平稳，并趋近于1，且ε越小，这种趋近程度越明显，误差也越小，表明波前愈合效应不明显，此时射线理论计算的走时误差较小.对比两种介质，相对于指数型自组织介质，高斯型自组织介质在该范围内的散点拟合斜率值变化范围更大.

图 9

Fig. 9

图 9 高斯型(a)和指数型(b)自组织介质互相关走时差与射线理论走时差散点拟合斜率图 横坐标为参数a/(λL)1/2，垂直虚线代表参数等于0.5；纵坐标为走时差散点图拟合斜率，水平虚线代表斜率为1.0；a1—a3代表不同标准差的高斯型自组织介质，b1—b3代表不同标准差的指数型自组织介质.Fig. 9 Best fitting scatterplot slope for the full suite of Gaussian(a) and exponential(b)self-organized media X-coordinate denotes parameter a/(λL)1/2， vertical dashed line denotes the parameter equal to 0.5; Z-coordinate denotes slope of best-fitting，horizontal dashed line denotes the slope equal to 1.0.

本文利用声波非均匀化多尺度方法对不同扰动强度及相关长度的高斯型和指数型自组织介质进行 了波场模拟，采用互相关方法计算背景速度场和附加自组织介质速度场之间波形的走时差，并与射线理论方法获得的走时差进行比较.结果表明，非均匀化多尺度方法是解决非均匀大尺度模型的有效波场模拟方法，节省计算时间的同时又可保持计算精度.介质相关长度a越小、波长λ越长，且传播距离L越远时，波前愈合效应在增强.对两种方法得到的走时差结果进行散点拟合与散点斜率统计表明，当a/(λL)^1/2≤0.5时，两种方法得到的走时差差别明显，随着该值逐渐增大，散点斜率趋近于1，且ε越小，这种趋近程度越明显，拟合的误差越小.而且高 斯型自组织介质在该值小于0.5时的散点拟合斜率值变化范围大于指数型自组织介质.即当a≤0.5(λL)^1/2时，波前愈合效应显著，射线理论的误差较大.