2012, Vol. 55
2. 中国石油大学(北京)油气资源与探测国家重点实验室, 北京 102249
2. State Key Laboratory of Petroleum Resource and Prospecting, China University of Petroleum (Beijing), Beijing 102249, China
在地震勘探中，地质体通常被一连串大的均匀､各向同性层近似替代.然而高分辨率技术的出现及发展､剪切波越来越多地被利用､以及人们开始对地层学细节上加以了注意，都需要对所做的假设重新进行考察.
对于均匀性的假设条件而言———首先，这样的模型没有考虑到由测井曲线和岩芯样本显示出的小尺度非均匀性;其次，从这种模型得到的合成地震数据同野外地震资料间存在着巨大差异;另外，地震波在地质体中传播时，其本身受到了多种尺度非均匀性的影响.一个取而代之的模型是同时包含大尺度和小尺度非均匀性的模型.大尺度非均匀性是地质体的平均特征，小尺度非均匀性则是加在这些平均特征基础之上的扰动(这种扰动的思想的提出可追述到20世纪五､六十年代[1-2])，由于这些小尺度非均匀性数量既多分布又不规则，则所能期望的就是从地震数据中重建它们的统计信息.因此，可用空间随机分布过程来描述小尺度非均匀性，并从统计意义上对其进行研究.deHoop等[3]给出了1D 精细层状介质模型，该模型的成层厚度甚至可小于实际声波测井的分辨率;Kerner[4]创新式地给出了1D 周期式及过渡式随机水平层状介质模型的建立方法，并证实了该模型的有效性;Frankel和Clayton[5]采用具有各向同性自相关函数的2D 随机介质模型来描述小尺度非均匀性，该模型较适合于对比分析大范围的地质区域;Ikelle等[6]考虑了小尺度非均匀性在水平方向和垂直方向的择优取向问题，给出了一种依据椭圆自相关函数建立2D 随机介质模型的方法;我国学者奚先等人，于21 世纪初期首次较清晰地给出了一种随机介质的建模步骤，并在该基础上进行了一系列的拓展研究[7-9]，朱生旺等[10]借助随机介质模型对孔洞型介质进行了描述.然而，目前较为流行的建模方法对小尺度非均匀性择优取向问题的考虑，仅仅局限在水平方向和垂直方向，灵活性不足;并且国内有关从统计意义上描述非均匀性的文献中，均忽略了一个至关重要的问题，即在离散波数域中计算自相关函数时会产生误差，该误差严重影响建模的有效性，本文第一个目的就是解决上述两方面问题.
对于各向同性的假设条件而言———早在20 世纪50年代就有学者提出质疑，当时的一个公认长波长理论为:一个由薄层组成的水平层状介质(其各个组成薄层为弹性､均匀且各向同性的)，当地震波长足够长时，这个介质就表现出一个对称轴垂直于成层面､均匀､横向各向同性､非频散的各向异性固体介质性质.Postma[11]证实了在成层介质只包含2种岩石成分情况下该理论的正确性;Backus[12]将该理论推广到任意多层以及组成层可以为横向各向同性的情况，并给出了具有里程碑意义的平均化计算公式;Berryman[13]基于Postma和Backus的理论，对成层介质的各向异性进行了研究，并认为近海区长波长各向异性对纵波的影响可以忽略，但对剪切波的影响是很明显的;Daley和Hron[14]首次在地震扰动记录方面对比分析了两种在数学上等价的模型(以双层介质为单元形成的周期式介质和与其相对应的等效均匀横向各向同性介质)的数值响应，发现模拟结果惊人的相似;Melia等[15]第一次采用实验室测量的手段，用环氧树脂及玻璃组成的模型对该理论进行验证，发现测量结果同理论值没有统计上的显著差异，并且随着测量频率的增加，所得速度值是从长波长理论值逐渐向利用时间平均关系预测的射线理论值变化，指出当地震波长λ 同成层单元长度d之比为10~100 时，长波长理论近似成立;Helbig[16]通过对周期式层状介质频散方程的估算，认为对于SH 波而言λ/d值至少为3时，才可利用长波长理论;Carcione等[17]用数值模拟方法研究周期式层状介质中波的传播，发现长波长理论的适用条件，对于环氧树脂-玻璃式周期结构来说λ/d=8，对于砂岩-石灰岩式周期结构λ/d=5~6;Marion等[18]和Mukerji等[19]均研究了射线理论与长波长理论的过渡带问题，发现过渡区间对于由钢和塑料组成的周期式层状介质而言大约为λ/d=8~15，在2D 的随机介质中大约为λ/d=1~10，随机介质中的d是建立模型时所采用的自相关函数的相关长度，也叫做非均匀性尺度;Griffiths和Steinke[20]给出了一个周期性层状介质中地震波传播的一般性理论;Stovas和Ursin[21]基于Griffiths和Steinke的研究成果，导出周期层状介质中反射和透射波场的显式表达式，指出反射系数以及周期性层状介质两组成部分的单程走时之比是两个控制波在周期式介质中传播的重要参数，在小反射系数情况，长波长理论和射线理论的过渡带发生在λ/d=4附近.可见，多数对于长波长理论的研究均局限在由两种成分组成的周期式层状介质中，并且学者们得到的长波长理论适用条件几乎各不相同，这正是本文对长波长理论进行进一步研究的原因.
本文在前人研究的基础上，首先给出一种形式多样化的且可更全面考虑小尺度非均匀性择优取向的随机介质建模新方法，并对离散波数域中计算自相关函数时生成的误差进行了分析及有效处理，展示了典型建模成图.其次着重讨论了地震波在水平层状介质中传播时的长波长理论，指出对该理论进行研究的重要意义———利于测井和地震资料的匹配.接下来，基于1D 随机层状介质，从局部各向异性曲线的角度对长波长理论进行了研究，指出了地震波在这种介质中传播的特点，并讨论了密度及泊松比扰动对各向异性因子曲线形态的影响.然后，利用加载了完全匹配层(PML)边界的高阶交错网格有限差分法，在长波长理论不能完全满足的情况下，对1D 和2D 随机介质进行数值模拟研究，分析地震记录特征.为了便于对比分析，文中也计算了水平层状介质在不同平均化长度下，通过长波长理论得到的等效介质的理论走时(走时的计算是在弱各向异性假设条件下完成的)，最后分析了长波长理论对随机介质模型的适用性，给出了有指导意义的结论.
2 小尺度非均匀性的描述 2.1 基本假设首先，引入一个被用来描述各向同性弹性介质的有限集m(x)= {mi(x)}，i= 1，2，3，…，集合中的元素mi(x)用来描述介质所具有的弹性参数，x= (x，z)为坐标矢量.假设集合m(x)可被表示成
|
(1) |
其中m0 代表介质的大尺度非均匀性，在此假定为均匀的，Δm代表平稳随机分布在均匀介质中的小尺度非均匀性(亦可称为扰动)，下面将集合m(x)称为随机介质.
由于小尺度非均匀性数目巨大且分布极其无序，则可利用它们的统计信息对其进行描述.另外，最可能从地震资料中恢复出的随机介质的统计特征是其低阶统计矩，尤其是前两阶矩，那么就可将研究限定在二阶统计学范围.不失一般性，本文设随机介质m(x)由纵波速度vP ､纵横波速度比vP/vS 及密度ρ 三个分量构成.则有:
|
(2) |
根据Birch原理[22]，可假定纵波速度和横波速度具有相同的扰动性质，且密度扰动仅简单同速度扰动具有线性关系，那么只讨论随机介质中元素vP(x)的建立即可.
对于vP(x)=vP0(x)+ΔvP(x)，由于其一阶统计矩为vP0(x)，则可得第一个约束条件:
|
(3) |
其中〈〉表示对括号中的量做加权平均.
随机介质的二阶矩是由相关函数控制的，在假定随机介质中的元素是相互独立的前提下，只需选定因子vP(x)的自相关函数C(x1，x2)= 〈vP(x1)vP(x2)〉.目前可选择的自相关函数形式多样，Mukerji[19]选择了各向同性高斯型自相关函数，该函数所描述的非均匀性没有择优取向;Ikelle等[6]考虑了非均匀性在水平及垂直方向上的择优取向问题，选择了一种椭圆式的指数型自相关函数;奚先和姚姚[8]选择了一种混合型自相关函数，这时椭圆式的高斯型和指数型自相关函数只是给定了粗糙度因子r时的特例.本文，为了更全面地考虑小尺度非均匀性的择优取向问题，给出如下的新型自相关函数:
|
(4) |
其中a､b是自相关长度因子，r是粗糙度因子，φ 是用于表示小尺度非均匀性沿顺时针方向旋转角度的角度因子，其取值范围是[0，2π).这样，就可通过适当地选择自相关长度及角度因子去描述介质的非均匀性为各向同性的､在xz平面内的任意两个垂直方向上适度延长了､或者在某一方向无限延伸了的介质.该自相关函数当φ =0､r=0时，退化为高斯型自相关函数，当φ=0､r=1时，退化为指数型自相关函数，很明显，这种新型自相关函数既保持了原有自相关函数的特点又丰富了小尺度非均匀性的择优取向.
2.2 建模步骤具备上述一系列条件后，即可按照如下步骤进行随机介质的建模:
①选定好自相关函数的因子后，通过快速傅里叶变换将Φ(x，z)从空间域变换到波数域得珟Φ(kx，kz);
②在区间[0，2π)中生成服从独立､均匀分布的随机数ξ;
③利用概率论中的性质“自相关函数的傅里叶变换是随机介质的功率谱"，即可得到随机介质的波数域表达式:
|
(5) |
④利用反傅里叶变换将珟S(kx，kz)从波数域变换到空间域得S(x，z)，在确保方程(3)得以满足的前提下，将随机介质归一化到想要的方差σP 即可.
显然上述步骤理应是在连续域中进行计算的，然而实际上的计算是离散的，则在计算的过程中必然会引入误差，以至影响建模的可信度.据文献调研发现，国内从统计意义上描述非均匀性的文献中，没有对此问题进行考虑及适当处理，在这种情况下，倘若所建立的模型明显受到误差的影响而严重偏离其应具有的性质，那么后续的研究成果将不具可信性，这是极其危险的.因此，本文将在下一部分中给出相应的处理方式.
2.3 误差的处理及建模成图由于假定小尺度非均匀性是平稳随机分布在均匀介质中的，也就得到了方程(3)所示的约束条件，然而误差的存在会使得该约束条件不能得到保证.为了较为直观地展示离散计算误差对随机介质的影响，下面将在1D 情况下对其进行讨论.
首先，利用建立随机介质的4个基本步骤生成图 1a所示的具有如下参数的1D 随机介质:由2400个间距Δz=1m 的网格组成､平均纵波速度vP0(z)=3500m/s､自相关长度为b=3 m､归一化方差为σP=0.1､粗糙度因子r=0.然后，为考查介质的稳定性，计算该介质的移动平均值，选取窗口大小为20个网格，得到图 1b所示的结果.可见移动平均值同均值间存在显著的差异，这势必会严重影响建模的可信度.为了将移动平均值同均值间的差异尽量减小，一个有效的方法是在波数域计算随机介质时引入锥形函数T(k)[23].此时1D 和2D 情况下，在波数域计算随机介质的表达式分别为
|
图 1 (a)离散计算生成误差未被压制时的ID随机介质;(b)窗口大小为20个网格时,该介质的移动平均值 Fig. 1 (a) A 1D random medium whose errors due to discrete computation are not suppressed;(b)The moving average value of the medium when the window length is the same as 20 grids; |
|
(6) |
|
(7) |
T(k)在1D 和2D 情况下的表达式分别为
|
(8) |
|
(9) |
其中kzmax 表示锥形函数的长度，具体选择细节参考Marple所出专著[23].
图 2a及图 2b分别展示了不改变图 1a所示介质参数的前提下，引入锥形函数后的随机介质成图和移动平均值成图，此时移动平均后的扰动程度已明显降低，相对于随机介质的扰动已经达到可以忽略的程度.从而，通过引入锥形函数，离散计算的误差已被有效的压制，从而确保了随机介质模型的可信度.
|
图 2 (a)压制了离散计算生成误差时的ID随机介质;(b)窗口大小为20个网格时,该介质的移动平均值 Fig. 2 (a) A ID random medium whose errors due to discrete computation are suppressed;(b) The moving average value of the medium when the window length is the same as 20 grids; |
为了便于理解锥形函数的性质，图 3给出了在建立大小为500×500个网格的随机模型时，简单选取对应半径为100 个网格的2D 锥形函数，可以看出其作用实际是压制了随机介质中低波数的分量.在接下来的研究内容中，所建立的模型均引入了锥形函数T(k).
|
图 3 针对网格大小为500X 500的随机模型所选取锥 形函数的图示,图中锥形函数半径为0.2 m _1,对应的 空间范围是100个网格 Fig. 3 Diagram of the 2D tapering function used to establish a 500X500 random model. The tapering radius is 0.2 m—1 corresponding to 100 grids' space |
由于不可能穷举本文所示建模方法可生成的所有模型，并且较为常用的两种自相关函数为高斯型和指数型自相关函数，因而图 4 仅给出这两种自相关函数情况下，较具代表性的模型.为便于对比分析，假设各模型均有vP0(x，z)=3500m/s､σP=0.1，每个模型的其他参数情况为:(a)a=3 m､b=3 m､r=1､φ=0;(b)a=20 m､b=10 m､r=1､φ=π/4;(c)a=3 m､b=3 m､r=0､φ=0;(d)a=18 m､b=6m､r=0､φ=π/3;(e)a=∞ m､b=3 m､r=0､φ=0;(f)a=3m､b=∞ m､r=0､φ=0.其中(a)､(c)所描述的小尺度非均匀性是各向同性的，(b)､(d)描述的非均匀性是在互相垂直的两个方向上适当延长了的，(e)､(f)则描述了非均匀性在某一方向无限延伸时的情况.显然，该建模方法通过角度因子φ 可以方便地将小尺度非均匀性的择优取向拓展到任意的两个互相垂直的方向上(不再只局限于水平和垂直方向上).另外，通过观察图 4a—4d可发现，对于大小相当的自相关长度值，依据高斯型自相关函数而建立出的模型所描述的小尺度非均匀性具有更为良好的视觉效果，并且Mukerji等[19]研究长波长理论的过渡带问题时，只采用了各向同性的高斯型自相关函数建立随机介质，因而本文在接下来的对于长波长理论的研究中，所采用的模型均设定r=0.在利用本文所建模型对长波长理论进行研究前，对射线理论､介质的弱各向异性､各向异性和频散的相互关系进行介绍及适当分析是很有必要的.
|
图 4 速度均值为3500 m/s､归一化方差为0.1前提下的六个典型建模成图,各介质的其他参数(a､､r和分别为(a)3 m､ 3 m､1､0;(b)20 m､10 m､1､n/4;(c)3 m､3 m､0､0;(d)18 m､6 m､0､n/3;(e)⑵m､3 m､0､0;(f)3 m､∞m､0､0. Fig. 4 Six representative random models all of which have the same average velocity (3500 m/s) and normalized variance (0. 1). The other parameters (a^b^r and cp) of each random medium are (a)3 m,3 m､1,0;(b)20 m,10 m,1,π/4;(c)3 m,3 m,0,0; (d)18 m,6 m,0,π/3;(e)⑵ ∞m,3 m,0,0;(f)3 m,⑵ m,0,0,respectively |
当水平层状介质的各向同性组成层相对于地震波长很长时(短波长极限)，地震波会在地层分界面处发生散射，且以时间平均速度(亦称射线速度)传播，这就是通常所说的射线理论;而当地震波长足够长时(长波长极限)，地震波在层状介质中传播时好像经过了一个对称轴垂直于成层面､均匀､非频散的横向各向同性固体介质，这则是被广为讨论研究的长波长理论.下面将介绍在长波长极限下，一种可简单而有效地处理水平层状介质的方法.
假设所研究的地质体为线性弹性的，则其中的每一个应力分量τji线性依赖于每个应变分量ζkl.由于每个方向的指标都可假定值1､2､3(分别代表x､y､z方向)，这样就有9个关系式，每个包含1 个应力分量和9个应变分量，9个关系式可简写为
|
(10) |
其中有关cijkl的3×3×3×3弹性模量张量刻画着介质的弹性性质.由应力的对称性τji=τij和应变的对称性＜ζkl=ζlk，可以通过下面的下标变换法:
|
(11) |
有效地简化介质弹性性质的表示形式，将3×3×3×3张量表示为6×6的矩阵Cαβ.每个特定的对称类型有着特定的非零独立的Cαβ 表示形式，其中各向同性介质和最简单的具有广泛地球物理适用性的均匀横向各向同性介质的张量矩阵分别为(由于对称性，仅给出上三角组成部分):
|
(12) |
|
(13) |
Backus[12]经论证得到结论:对于一个由多层横向各向同性(每层的对称轴方向都垂直于层面)材料组成的层状介质，在长波长极限条件，可等效成对称轴垂直于成层面的横向各向同性介质，即VTI介质.相应的等效弹性模量可由每个组成层的弹性系数计算出，此时的等效张量矩阵及其等效模量分别为
|
(14) |
|
(15) |
Backus指出，方程(15)也适合于成层为各向同性情况，Levin[24]则给出当成层为各向同性时，用各层纵波速度vP､横波速度vS 和密度ρ 来表示等效弹性模量的表达式:
|
(16) |
另外，层状介质的等效密度表达式为
|
(17) |
那么，平行于和垂直于层状介质传播的剪切波和压缩波速度可以通过等效密度ρEM 和等效弹性模量C11､C33､C44 及C66 表示为
|
(18) |
其中下标字母v､h分别表示速度传播方向垂直于和平行于层状介质，SH 和SV 分别表明剪切波水平向极化和竖直向极化.
总体来讲，长波长(极限)条件下的剪切波和压缩波速度，平行于和垂直于由两种物质组成的层状介质传播时，理论上，可以利用各向同性层的密度和弹性波速度，以及它们在合成固体中的百分含量而方便地计算得到.
图 5显示了在长波长极限下，水平层状介质等效成VTI介质的示意图，其中层状介质的厚度L也叫做平均化长度.
|
图 5 水平层状介质a及其在长波长极限 条件下等效成的VTI介质(b) Fig. 5 A horizontal layered medium (a) and its equivalent VTI medium (b) when the long wavelength limit condition is satisfied |
然而，该理论的适用条件“地震波波长足够长"过于含糊，为了能够给出一个定性甚至定量的分析，许多学者从多种角度进行了相关研究.分歧较大､讨论较多的一个话题是用波长同非均匀性尺度之比λ/d来确定长波长极限.这里的非均匀性尺度d，对于周期成层介质来说是组成该介质的成层单元长度，对于随机介质来说通常被认为是描述介质时所采用的空间自相关长度.目前较普遍被认可的一种观点是“当波长与非均匀性尺度之比满足λ/d≥10时，长波长理论近似成立".在假定此观点正确的前提下，本文可容易地通过如下示例阐明对长波长理论进行细致研究的重要意义.
假设用单频震源来测量纵波在一个周期性层状介质(成层单元各成分密度相同)中的垂向传播走时，介质的成层单元由厚度d1 =5 m､纵波速度vP1 =3000m/s的水平层和厚度d2=5m､纵波速度vP2=4000m/s的水平层组成，则d=10 m.若采用地震勘探频段范围的震源进行测量，如f震=30 Hz, 则有λ震min=100m､λ震max=133 m;若震源的频率处于测井频段，如f井=3000 Hz, 则有λ井min =1.0 m､λ井max=1.3m.显然对于地震勘探频率有λ震min/d=10，波的传播满足长波长理论;对于测井勘探频率有λ井max/d=0.13，在认为λ井max/d=0.13 1 时，波的传播满足射线理论.由于射线理论速度VRT 大于等效介质理论速度VEMT[18](详细推导见附录A)，则面对同样的一个地质体，测井测量速度要高于地震测量速度，这就是井地资料匹配不准的原因之一.
可见，若能够解决长波长理论的适用条件问题，一个重要作用就是利于测井､地震资料的匹配，这也是本文对该理论进行讨论的根本原因.
3.2 弱各向异性已经知道，当长波长理论得以满足时，层状介质可以等效成一个对称轴垂直于成层面､均匀､非频散的VTI介质.对于VTI介质而言，地震波在其中的每个传播方向上均有三个解，即一个准纵波vqP，一个准横波vqSV 和一个横波vSH.Daley和Horn[25]给出了三个相速度的表达式:
|
(19a) |
|
(19b) |
|
(19c) |
|
(19d) |
其中ρ 为VTI介质密度，θ 为相位角(图 6 给出了VTI介质中相位角和群角之间的关系)，D(θ)是一个复杂表达式的缩记，正是它的复杂性，成为了利用横向各向同性模型分析地震数据的主要障碍.
|
图 6 VTI介质中相位角0和群角0的关系 Fig. 6 Relationship between phase angle 0 and group angle 0 in the VTI medium |
Thomsen[26]通过对比分析方程(12)､(13)，首先提出用三个各向异性因子ε､γ､δ*及垂向纵波速度vP0和垂向横波速度vS0来精确描述相速度的方法，此时相速度的表达式为
|
(20a) |
|
(20b) |
|
(20c) |
|
(20d) |
其中:
|
(21a) |
|
(21b) |
|
(21c) |
|
(21d) |
|
(21e) |
接下来，Thomsen通过对一系列沉积岩的各向异性数据进行观察，发现绝大多数岩石的各向异性因子取值均小于0.2，即岩石的各向异性普遍呈“弱"状态.从而Thomsen指出绝大多数弹性介质所具有的各向异性是弱的，并在此结论下对方程(20d)进行了一系列适当化简得:
|
(22) |
将方程(22)代入方程(20a)和(20b)，并进行线性化处理即得到一组适用于弱各向异性VTI介质的关于相速度的表达式:
|
(23a) |
|
(23b) |
|
(23c) |
其中:
|
(23d) |
可见，在弱各向异性条件下，VTI介质中的相速度同垂向纵波速度vP0､垂向横波速度vS0､各向异性因子ε､γ､δ 及相位角θ 的依赖关系就变得非常清晰，利用各向异性模型分析地震数据的难度被显著降低.由于Thomsen的结论是建立在大量的数据支撑下的，具有很高的可信度，因此本文在后面的研究中将会利用Thomsen所提出的各向异性参数对长波长理论进行讨论.为了便于对比分析，下面先把Thomsen对各向异性数据的观察结果概括如下:①对90%的样本，纵波各向异性因子ε 取值范围处于0到0.2之间，横波各向异性因子γ 取值范围处于0到0.3之间;②介质中横向速度大于垂向速度的只有很少的一部分，且对应的负因子的绝对值较小;③得到了或正或负的δ*和δ 因子，且大部分因子的绝对值小于0.3.
3.3 各向异性和频散的相互关系对于地震波在弹性介质中的传播而言，各向异性和频散是同一物理事件(相速度对波矢量k的依赖关系)的两种表现形式.如果波速只依赖于波矢量k的大小，则谈论频散;如果只依赖于波矢量k的方向，则提及各向异性.波速对波矢量具有依赖性的一个结果就是单频平面波(具有极小带宽的平面波)的速度是不同于能量传播速度的(由于假定介质为弹性，则能量速度等于群速度).在仅频散的情况下，这两个速度具有相同的方向，但有着不同的大小;在只有各向异性的情况下，它们有着不同的方向和不同的大小，此时单频平面波的速度为能量传播速度沿波法线方向上的投影.很明显，各向异性可以是频率依赖的，频散效果可以是方向依赖的，也就是说各向异性和频散可以按照不同的组合形式而共同发生.
随之即来的问题则是:长波长理论所得到的等效介质是一个非频散的各向异性介质，那么当长波长理论刚好不能满足时，层状介质会表现出什么性质呢?比如说，该介质能否依然是各向异性的?若是各向异性的，那么是否和频散同时发生?Helbig[16]首次从频散方程出发，研究了由两种物质组成的周期式层状介质中SH 波的频散曲线，指出当长波长理论不能完全满足时，各向异性和频散是共同发生的.虽然该结论是在研究周期式层状介质中得出的，但也颇具参考意义;Kerner[4]利用泊松型和褶积型随机层状介质讨论了该问题，对于剪切波的波场记录，得到了和Helbig相类似的结论.在下部分内容中，将在本文所建立的随机介质的基础上，对上述问题进行分析研究.
4 基于随机介质的长波长理论研究 4.1 随机层状介质的局部各向异性因子曲线研究长波长理论的方法有多种，本文先从随机层状介质的局部各向异性因子曲线的角度对其进行分析.若对于整个1D 随机介质的移动平均过程产生了一个几乎均匀的介质，长波长理论则为1D 随机介质生成了一个等效的均匀横向各向同性弹性介质(这里需要注意的一点是，利用该理论时，进行移动平均的目标长度也就是平均化长度L，一定要小于地震波信号的主波长).下面将计算图 2a所示的随机层状介质(即表 1中的media1)的局部各向异性因子曲线.若选取平均化长度L=60m, 则计算步骤为:根据方程(16)计算介质从深度0 m 到60 m区域内的等效VTI介质的弹性模量C11､C12､C33､C13､C44 及C66，代入方程(21a)—(21c)､(23d)计算出此深度范围内的等效各向异性参数ε､γ､δ* 及δ;然后再以同样的方式计算介质从深度1 m 到61 m范围内的等效各向异性参数;…;以1 m 为移动步长，直至计算完2340 m 到2400 m 的深度范围，即得局部各向异性因子曲线图.这种做法“模拟"了地震波在介质中传播时的动态平均化过程.
|
表 1 随机介质模型及其对应的参数 (各个模型的角度因子P均为0) Table 1 Random media and their parameters (each medium s angle factor 9 is equal to 0) |
图 7a和7b分别给出了平均化长度L为60 m和120m 两种情况下，随机层状介质media1 的局部各向异性曲线成图，其中“epsilon"曲线代表ε､“gamma"曲线代表γ､“delta1"曲线代表δ*､“delta2"曲线则代表δ.对于两种平均化长度，各向异性因子曲线均显示出随深度的变化而变化，并且即使当L=120m 时，曲线依旧存在明显扰动.由于模型的非均匀性尺度d=b=3 m, 这意味着若想得到整个1D 随机层状介质的等效均匀横向各向同性介质，平均化长度L要远大于非均匀性尺度，地震信号的主波长也就要远大于非均匀性尺度.如果对该模型进行数值模拟，即使波长相对非均匀性尺度而言很长，散射现象多多少少也会出现;倘若波长相对于非均匀性而言不太长，散射现象应明显发生.图 7c则给出了media1 在平均化长度L=1000 m, 即L/d=333情况下的局部各向异性参数曲线，可发现微弱扰动依旧存在.然而，尽管对于L为60m 和120m情况的平均化过程，没有产生均匀的各向异性等效介质，但生成的介质依旧是各向异性的.
|
图 7 media 1介质在平均化长度分别为60 m (a)､120 m (b)和1000 m (c)时的局部各向异性曲线以及 media 1'(d)､media l"(e)和media ]/〃(〇介质在平均化长度为60 m时的局部各向异性曲线 Fig. 7 Local anisotropy factors' curves of media 1 when the average length are 60 m (a)，120 m (b) and 1000 m (c) respectively.Local anisotropy factors' curves ofmedia 1 (d)、media 1〃（e) and media 1〃’ ⑴ when the average length is 60 m |
通过观察还可发现:因子γ 的曲线值是最高的(暗示横向各向同性等效介质对于剪切波来说，其非均匀程度更高)，其曲线形态和ε 类似､同δ* 成“对称"关系;δ* 值是最小的，其曲线的扰动幅度相对γ来说较小，但同ε相当;δ 值略高于δ* 值，但其曲线的扰动非常微弱(暗示横向各向同性等效介质对于δ 因子来讲，其非均匀程度几乎可以忽略).
同时，为了研究密度扰动和纵横波速度比的扰动对局部各向异性曲线的影响，在平均化长度为60m和介质media1的其他参数性质不变的前提下，图 7d给出了密度标准差σρ =0时，即media1′介质的各向异性因子曲线图;图 7e给出了纵横波速度比的标准差σP/S=0时，即media1″介质的各向异性曲线图;而图 7f则给出了σρ =0和σP/S=0共同发生时，即media1介质的各向异性因子曲线图.由于纵横波速度比的变化在一定程度上反映了介质泊松比的变化，那么经对比分析各图可得以下结论:层状介质的密度是否扰动，对δ 曲线几乎无影响，若层状介质的各层密度变得恒定，则ε､γ､δ* 曲线均会向0值靠近;若随机介质的泊松比突然变得恒定，除了ε､γ､δ* 曲线均会向0值靠近外，δ 曲线会变成一条直线且γ 值会降低到和ε 相同的水平，可见泊松比恒定的假设对曲线的影响是显著的.
值得注意的一点是，计算得到的各向异性参数值均在0.1以下，也就是说在给定平均化长度下计算出来的等效介质的各向异性是弱的，这同Thomsen的结论相吻合.
4.2 随机层状介质的数值模拟前面已经讨论过各向异性和频散的关系，并且计算了图 2a所示的随机介质模型在多种平均化长度下的局部各向异性曲线，发现平均化过程生成的介质是非均匀各向异性介质，那么当采取一定的波长(随机层状介质大小为主波长的若干倍)，对随机层状介质进行数值模拟时，是否会有各向异性､频散或散射现象发生呢?为给出解答，下面将对随机层状介质进行数值模拟研究.
时间及空间的离散化会引起数值频散以及数值各向异性，而本文研究的目标是物理频散及物理各向异性.从微扰理论推测，这些物理效果可能只达到几个百分点的影响程度[27]，这意味着，如果在模拟中采用低阶的有限差分方法，这些影响效果可能会被数值误差干扰，高阶有限差分法则可在不增加网格数目的情况下解决模拟精度问题;另外，随机介质中泊松比可以在很大一个范围内变化，且包含大泊松比值情况，而交错网格法[28]则可在大泊松比情况下，减少数值频散;Levander[29]采用四阶交错网格有限差分法进行数值模拟时，发现当每一个波长中的网格点数多于十个时，数值频散和数值各向异性均可忽略不计.综上，本文选取高阶交错网格法进行数值模拟研究，另外，在模拟过程中均采用了放射纵波的爆炸式震源，子波为雷克子波.
在假设介质是线性弹性､各向同性的､且在t=0时刻介质中的粒子速度和应力张量均为零，那么介质的弹性波方程组可写为
|
(24) |
其中vx､vz是速度向量，τxx､τzz､τxz为应力张量，λ､μ 是拉梅参数，ρ 为密度.
为更好地研究随机介质中的波场形态，本文在利用六阶交错网格有限差分法对方程(24)进行求解的同时，在模型四周加载了具有良好吸收效果的PML 吸收边界[30].
不失一般性，首先选取图 2a 所示介质的0~1000m 深度范围进行垂直地震剖面(VSP)数值模拟.由于随机层状介质的平均纵波速度vP0(z)=3500m/s, 为了使波长同非均匀性尺度之比较大，令震源主频fm=30Hz.对于纵波，则有λm0/d=38.9，相比于λ/d=10来说，主波长λm0=116.7m 已经很长，显然此时介质既不满足强散射理论(由于d/λm0=3/116.7\[\ll \]1)又不满足长波长理论(由于在平均化长度L=120m＞λm0 时，介质的局部各向异性曲线扰动依旧明显).
图 8给出了依据图 2a所示介质0~1000 m 范围，在偏移距为250m 情况下的声波VSP地震记录(a)､(b)和弹性波VSP地震记录(c)､(d)，其中指示因子X表示地震记录的水平方向分量，Z表示垂向分量.注:在此之后的地震记录成图的指示因子X､Z的含义均与图 8相同.通过观察可发现，虽然声波同弹性波地震记录的初至吻合良好，但对于仅激发出纵波的震源，弹性波模拟时观测到的转换波能量也是相当可观的，因此由于缺乏转换波信息，声波响应不能较完善地描述这种介质;另外，也是本文更为关心的现象，在这种强散射理论不能满足的情况下，无论纯声波响应还是弹性波响应，每道地震记录中都存在微弱散射，且能够明显观察到上行波，但是无论对于记录中的纵波还是横波均没有发现显著频散现象.
|
图 8 图2a所示介质0〜1000 m范围,在偏移距为250 m情况下的声波VSP地震记录(a)､(b)和弹性波VSP地震记录(c)､(d)指示因子X表示地震记录的水平方向分量､Z表示垂向分量 Fig. 8 Seismograms for acoustic wave VSP (a)， (b) and elastic wave VSP (c)， (d) (offset：250 m) recorded in the medium shown in Fig. 2a. The recording depth ranges from 0 to 1000 m The indicative factor X represents seismograms； component of horizontal direction and Z represents seismograms； component of vertical direction. |
图 9给出依据图 2a所示介质0~1000m 范围，在偏移距为500m 情况下的弹性波VSP地震记录.图中除了依旧可清楚观察到散射和上行波外，还可发现随着偏移距的增加，地震记录的垂向分量成图中，可明显记录到横波能量的深度范围有所增加，且上行波记录变得不明显.
|
图 9 图2a所示介质0〜1000 m范围,在偏移距为500 m情况下的弹性波VSP地震记录 Fig. 9 Seismograms for elastic wave VSP (offset：500 m) recorded in the medium showed in Fig. 2a. The recording depth ranges from 0 to 1000 m |
可见在λm0/d=38.9 的条件下，并没有观测到强散射及显著频散现象.为了对比分析，本文又建立了由2400 个间距Δz=0.25 m 的网格组成的随机层状介质(其具体参数参考表 1中的media2)进行数值模拟研究;模拟中令震源主频fm=300Hz, 则对于纵波有λm0/d=7.8，显然不满足λ/d≥10，此时波长相对于非均匀性尺度而言不是很长.那么，当波在介质中传播时，各向异性､频散和明显散射现象应同时发生.由于VSP地震记录不利于观测各向异性现象，因而设计了如图 10 所示的观测系统，其中震源位于x=z=100m 处，50个检波器均匀分布在以震源为中心，半径r=351 m 的1/4圆弧上，则震源到检波器的方向同直角坐标系垂直方向之间的夹角ψ(即图 6所示的群角)的变化范围是0°到90°.这样就可方便地对地震波在随机层状介质中传播时的各向异性､频散及散射现象进行观测.
|
图 10 观测系统示意图.其中黑点代表震源位置，星号代表检波器位置.检波器位于以震源为中心，351m为 半径的1/4圆弧上 Fig. 10 Diagram of the recording geometry，where the black spot denotes the source' s position and the asterisks denote the receivers' positions. Receivers are located in a one-fourth circle of which the center is the source and the radius is 351 m |
图 11首先给出图 10所示观测系统情况下，地震波在vP=3500m/s､vS=2333m/s､ρ=2.6g/cm3的均匀各向同性介质中传播时的地震记录，实线代表射线理论预测的走时，其结果同数值模拟结果有着良好的一致性;另外，可发现水平方向分量的强度随着入射角ψ 的增加而增加，垂向分量的强度随着ψ 的增加而减少，地震记录形态完全符合爆炸式震源的特点.注:从此以后，本文中的地震记录成图都是在震源主频为300 Hz､观测系统为图 10 所示的情况下得到的.
|
图 11 图10所示观测系统,震源主频为300 Hz条件下,地震波在纵波速度为3500 m/s､横波速度为2333 m/s､密度为2. 6 g/cm3的各向同性均勻介质中传播时的地震记录实线代表该均匀介质在射线理论预测下得到的走时. 实线代表该均匀介质在射线理论预测下得到的走时 Fig. 11 Seismograms for seismic wave propagating in a homogeneous isotropic elastic medium where the P-wave velocity is 3500 m/s,S-wave velocity is 2333 m/s and density is 2. 6 g/cm3. The main frequency of the source s 300 Hz The solid line represents traveltimes predicted by the ray theory. |
图 12 给出地震波在media2 介质中传播时的地震记录，记录的走时在一定程度上具有随入射角ψ 增加而减小的趋势，这种特点正是各向异性的表现;另外记录中展示出了介质的强散射现象，其中无论是地震记录的水平方向分量还是垂向分量，大入射角的尾波能量要显著高于小入射角情况;而且地震波在不同的传播方向上的初至脉冲长度是不同的，这表明了频散现象的发生.可见，对随机层状介质进行数值模拟研究，当主波长同非均匀性尺度之比降低到7.8(小于10)时，各向异性､散射及频散效应均会显著地显示在地震记录中.
|
图 12 地震波在media 2介质中传播时的地震记录,实线含义同图11 Fig. 12 Seismograms for seismic wave propagating in media 2. The solid line has the same meaning as in Fig.11 |
图 13给出了该随机层状介质以100 m 深度处为起点，平均化长度分别为L1=5 m､L2=50 m 以及L3=250m 情况下，依据长波长理论法得到的等效介质在图 10所示观测系统下的理论走时，走时的计算是在弱各向异性的假设条件下得到的，具体的计算步骤请参考附录B.可以发现，数值模拟结果同L1=5m 情况下的走时差异明显，同L3=250 m 情况下的差异最小，虽然这在一定程度上表明在给定起点的情况下，大平均化长度的理论值更接近数值模拟结果，然而并不说明可以用大平均化长度的理论值来描述地震波在介质中传播的实际情况.因为对于群角度在18°到72°范围内的实测走时都要显著低于长波长理论预测值.即可认为，在λm0/d=7.8的情况下，除了接近水平方向和垂直方向传播的走时外，大平均化长度下的长波长理论是不能准确描述地震波在介质中传播的走时的;在面对地震记录中初至脉冲的拉伸和能量显著的尾波时，长波长理论也不能给予解释.
|
图 13 media 2介质以100 m深度处为起点，在平均化长度分别为L1=5m、L2=50m、L3 =250 m情况下，依 据长波长理论法得到的等效介质在图10所示观测系统 下的理论走时，预测出的走时是在弱各向异性假设条件下得到的 Fig. 13 From the depth of 100 meters of media 2 , with the average length to be L1 = 5 m、L2 = 50 m and L3 = 250 m，respectively，we apply the long wavelength theory to obtain the effective medium whose theoretic traveltimes are recorded by using the recording geometry showed in Fig. 10. The predicted traveltimes are got under the hypothesis that the anisotropy is weak |
以上关于长波长理论的分析都是基于1D 随机介质的，下面将从2D 随机介质出发，对该理论进行研究.为了便于同非均匀性尺度为d=1.5 m 的media2介质的模拟结果进行对比研究，建立表 1中所示的media3—6四个2D 随机介质，介质网格间距Δx=Δz=0.25 m, 大小为2400×2400 个网格.这一系列2D 随机介质具有相同的自相关长度因子b=1.5m, 且另一个因子a的值从30m 逐渐降低至1.5m, 可见media6介质的非均匀性尺度也为d=1.5m.这样不仅可以对比具有不同维数､相同非均匀性尺度的media2和media6介质的地震记录差异，还可以研究从media2到media6介质间的过渡介质的地震记录情况.
图 14给出了图 10所示观测系统，震源主频为300Hz条件下，地震波在media3—6中传播时的地震记录.其中(a)､(b)对应于media3介质，(c)､(d)对应media4介质，(e)､(f)对应media5介质，(g)､(h)对应media6介质，其中实线的含义同图 11.通过对比分析图 11､12和14可以观察到，无论是1D还是2D 随机介质，地震记录中始终有能量很强的尾波存在，即存在散射;相对于图 11 中的初始脉冲而言，每个介质中的初始脉冲都已经拓宽了(由于多次反射或多次散射的作用);从media3介质逐渐过渡到media6 介质的过程中，地震记录的各向异性强度逐渐降低，media6 介质展示出了各向同性性质.
|
图 14 地震波在 media 3(a)(b)､media 4(c)(d)､media 5(e)⑴和 media 6(g)(h)介质中传播时的地震记录,实线含义同图11 Fig. 14 Seismograms for seismcc wave propagating in media 3 (a)(b)，media 4 (c)(d), media 5 (e)(f)and media 6 (g)(h). The solid line has the same meaning as in Fig.11 |
显然，在本文给出的数值模拟条件下，具有相同非均匀性尺度d的1D 和2D 随机介质所具有的性质是截然不同的:media2 介质具有显著各向异性性质，其尾波的能量主要集中在大入射角部分，且其不同传播方向上的初至脉冲长度几乎各不相同;而media6介质展示出各向同性性质，其尾波的能量均匀分布在各个传播角度，且其不同传播方向上的初至脉冲长度几乎相同;另外，值得注意的一点是，由media6 介质得到的走时要大于实线所示走时，并且要小于地震波在media2介质中以入射角ψ=0°传播时的走时.
由于尾波的强度几乎和初至脉冲相当，不可忽略，下面本文再对尾波特点进行讨论.
对于1D 的media2介质，多数尾波能量都处于大入射角范围，由于1D 介质的参数仅仅是深度依赖的，那么在大角度尾波的能量本质上是由超过临界角的隧道波(较早抵达)及导波(较晚抵达)引起的;对于2D 随机介质，尾波的能量本质上是由多次散射引起的，随着自相关长度因子a逐渐降低到和因子b相等，尾波能量从集中于大入射角范围变为均匀分布(由于a=b时，模型的非均匀性为各向同性的了).对于media6 介质，尾波能量变得均匀分布的现象也证明了本文建模方法的逼真性.
5 结 论本文主要研究了如下两方面内容:
首先，在前人的研究基础之上，为了更逼真地描述地质体中的小尺度非均匀性，给出一种形式多样化的且拓展了小尺度非均匀性择优取向的随机介质建模新方法，并指出了多年来国内有关随机介质建模文献中存在的隐患，即没有对离散计算随机介质模型时所产生的误差进行适当分析及处理.本文则通过引入一种锥形函数，很好地解决了建模误差的问题，使得随机介质模型的可信度更高.
其次，由于对长波长理论的细致研究有着利于测井与地震资料匹配的重要意义，并且多数对于该理论的研究都是简单的基于周期性层状介质的(同实际地质体情况的差异较大)，本文则以所建立的高斯型随机介质模型为基础，分别从1D 随机介质的局部各向异性曲线的角度以及对1D 和2D 随机介质进行数值模拟的角度来研究长波长理论的适用性.主要得到以下几点结论:若想得到整个1D 随机层状介质的一个等效的均匀横向各向同性介质，则需要一个远大于(至少相差两个数量级)非均匀性的尺度的平均化长度，若采用小于该平均化长度的波长对随机介质进行数值模拟，散射现象必然发生;1D 随机介质的局部各向异性因子曲线的绝对值均小于0.2，这同Thomsen对各向异性数据的观察结果相吻合;由于缺乏转换波信息，声波响应不能较完善地描述随机介质中的波动特征;对于1D 随机介质，当主波长同小尺度非均匀性之比较大时(λm0/d=38.9)地震记录中几乎观察不到频散现象，但依旧存在微弱散射，当比值较小时(λm0/d=7.8)各向异性､频散及散射均可被明显观测到，此时长波长理论既不能准确描述地震波在介质中传播的走时，又不能解释地震记录中的脉冲拓宽和尾波现象;当比值较小时(λm0/d=7.8)，具有相同非均匀性尺度d的1D和2D 随机介质所具有的性质是截然不同的，最显著的差异为1D 随机介质表现出显著各向异性性质，尾波能量集中在大入射角部分，而2D 随机介质表现出各向同性性质，并且其尾波能量均匀分布在各个入射角度(这种表现也说明了本文所建立模型的有效性)，并且群角度ψ=0°时，2D 介质中的群速度大于1D 介质中的群速度.
通过以上研究成果可见，如果随机介质模型可以真实地反映地质体性质，那么当主波长与小尺度非均匀性之比大于某一个给定值时(例如λm0/d=10)，就简单地认为长波长理论是成立的，这是有待商榷的，本文的研究结果表明该比值的差异要达到两个数量级，并且介质的维数也是一个显著的影响因素.
本文关于长波长理论问题的研究，主要是基于1D 随机介质模型的，后续的文章将会以2D 随机介质为主要出发点，对长波长理论进行进一步研究.
6 附录A:射线理论速度和等效介质理论速度的比较当考虑纵波垂直于由多层各向同性层组成的水平层状介质进行传播时，根据方程(16)､(17)和(18)可知，在长波长理论成立的条件下，即λ/d\[\gg \] 1 时，纵波垂直于介质进行传播的等效介质理论速度VEMT 等于VP，v, 即:
|
(A1) |
在射线理论成立的条件下，即λ/d\[\ll \]1时，纵波垂直于介质进行传播的走时是波在每一层介质中传播时间的总和，则可得关于射线理论速度VRT 的表达式为
|
(A2) |
经推导分析以上两式可以证明VRT >VEMT.
当所考虑的水平层状介质是一个由等体积分数的两种物质组成的周期性层状介质时，通过方程(A1)和(A2)可以推导出将两种速度联系起来的简单而又具有利用价值的关于VRT/VEMT 的表达式.假设两种物质的密度分别为ρ1 和ρ2､纵波速度分别为vP1 和vP2，并且令ρr =ρ1/ρ2，vPr =vP1/vP2.将上述符号代入方程(A1)和(A2)并联立可得:
|
(A3) |
并且当ρr ≈1时，进一步可得:
|
(A4) |
根据方程(16)计算media2介质从深度100 m为起点，平均化长度为L时的等效VTI介质的弹性模量C11､C12､C33､C13､C44､C66 及等效密度ρEM .由Thomsen[26]理论，在弱各向异性假设条件下，VTI介质的准纵波相速度的计算公式可写为
|
(B1) |
其中:
|
(B2) |
|
(B3) |
|
(B4) |
Thomsen进一步指出，对于一个给定的群角度ψqP，当利用(B5)式计算出相位角θqP 时:
|
(B5) |
可利用方程式(B1)及方程式
|
(B6) |
计算出VTI介质的群速度VqP(ψqP).则对于给定群角度ψqP 的走时为
|
(B7) |
其中r为炮检距，本文中群角度ψqP 的变化范围是0°~90°.
致谢感谢中国石油大学(北京)CNPC 物探重点实验室提供了良好的科研环境;特别感谢中国石油大学(北京)CNPC 物探重点实验室沈金松教授提供的帮助.
| [1] | Bergmann P G. Propagation of radiation in a medium with Random inhomogeneities. Physical Review , 1946, 70(7-8): 486-492. DOI:10.1103/PhysRev.70.486 |
| [2] | Aki K. Analysis of the seismic coda of local earthquakes as scattered waves. J. Geophys. Res. , 1969, 74(2): 615-631. DOI:10.1029/JB074i002p00615 |
| [3] | de Hoop M V, Burridge R, Chang H W. Wave propagation with tunneling in a highly discontinuous layered medium. Wave Motion , 1991, 13(4): 307-329. DOI:10.1016/0165-2125(91)90067-X |
| [4] | Kerner C. Anisotropy in sedimentary rocks modeled as random media. Geophysics , 1992, 57(4): 564-576. DOI:10.1190/1.1443270 |
| [5] | Frankel A, Clayton R W. A finite-difference simulation of wave propagation in two-dimensional random media. Bull. Seis. Soc. Am. , 1984, 74(6): 2467-2186. |
| [6] | Ikelle L T, Yung S K, Daube F. 2-D random media with ellipsoidal autocorrelation functions. Geophysics , 1993, 58(9): 1359-1372. DOI:10.1190/1.1443518 |
| [7] | 奚先, 姚姚. 二维随机介质及波动方程正演模拟. 石油地球物理勘探 , 2001, 36(5): 546–552. Xi X, Yao Y. 2-D random media and wave equation forward modeling. Oil Geophysical Prospecting (in Chinese) (in Chinese) , 2001, 36(5): 546-552. |
| [8] | 奚先, 姚姚. 随机介质模型的模拟与混合型随机介质. 地球科学--中国地质大学学报 , 2002, 27(1): 57–67. Xi X, Yao Y. Simulations of random medium model and intermixed random medium. Earth Science-Journal of China University of Geosciences (in Chinese) (in Chinese) , 2002, 27(1): 57-67. |
| [9] | 奚先, 姚姚. 二维横各向同性弹性随机介质中的波场特征. 地球物理学进展 , 2004, 19(4): 924–932. Xi X, Yao Y. The wave field characteristics of 2-D transversely isotropic elastic random medium. Progress in Geophys. (in Chinese) (in Chinese) , 2004, 19(4): 924-932. |
| [10] | 朱生旺, 魏修成, 曲寿利, 等. 用随机介质模型方法描述孔洞型油气储层. 地质学报 , 2008, 82(3): 420–427. Zhu S W, Wei X C, Qu L, et al. Description of the carbonate karst reservoir with random media model. Acta Geologica Sinica (in Chinese) (in Chinese) , 2008, 82(3): 420-427. |
| [11] | Postma G W. Wave propagation in a stratified medium. Geophysics , 1955, 20(4): 780-806. DOI:10.1190/1.1438187 |
| [12] | Backus G. Long-wave elastic anisotropy produced by horizontal layering. .J. Geophys. Res , 1962, 67(11): 4427-4440. DOI:10.1029/JZ067i011p04427 |
| [13] | Berryman J G. Long-wave elastic anisotropy in transversely isotropic media. Geophysics , 1979, 44(5): 896-917. DOI:10.1190/1.1440984 |
| [14] | Daley P F, Hron F. Computation of synthetic seismograms for a thin layered periodic structure. Canadian Journal of Earth Sciences , 1982, 19(7): 1449-1453. |
| [15] | Melia P J, Carlson R L. An experimental test of P-wave anisotropy in stratified media. Geophysics , 1984, 49(4): 374-378. DOI:10.1190/1.1441673 |
| [16] | Helbig K. Anisotropy and dispersion in periodically layered media. Geophysics , 1984, 49(4): 364-373. DOI:10.1190/1.1441672 |
| [17] | Carcione J M, Kosloff D, Behle A. Long-wave anisotropy in stratified media: A numerical test. Geophysics , 1991, 56(2): 245-254. DOI:10.1190/1.1443037 |
| [18] | Marion D, Mukerji T, Mavko G. Scale effects on velocity dispersion: From ray to effective medium theories in stratified media. Geophysics , 1994, 59(10): 1613-1619. DOI:10.1190/1.1443550 |
| [19] | Mukerji T, Mavko G, Mujica D, et al. Scale-dependent seismic velocity in heterogeneous media. Geophysics , 1995, 60(4): 1222-1233. DOI:10.1190/1.1443851 |
| [20] | Griffiths D J, Steinke C A. Waves in locally periodic media. American Journal of Physics , 2000, 69(2): 137-154. |
| [21] | Stovas A, Ursin B. Equivalent time-average and effective medium for periodic layers. Geophys. Prosp , 2007, 55(6): 871-882. DOI:10.1111/gpr.2007.55.issue-6 |
| [22] | Birch F. The velocity of compressional waves in rocks to 10 kilobars, Part 2. Geophys. Res , 1961, 66(7): 2199-2224. DOI:10.1029/JZ066i007p02199 |
| [23] | Marple S L. Digital Spectral Analysis with Applications. New Jersey: Prentice-Hall, 1987. http://cn.bing.com/academic/profile?id=2152783555&encoded=0&v=paper_preview&mkt=zh-cn |
| [24] | Levin F K. Seismic velocities in transversely isotropic media. Geophysics , 1979, 44(5): 918-936. DOI:10.1190/1.1440985 |
| [25] | Daley P F, Hron F. Reflection and transmission coefficients for transversely isotropic media. Bull. Seis. Soc. Am , 1977, 67(3): 661-675. |
| [26] | Thomsen L. Weak elastic anisotropy. Geophysics , 1986, 51(10): 1954-1966. DOI:10.1190/1.1442051 |
| [27] | Kneib G, Pestana R, Kerner C. Accurate and efficient seismic modelling in random media. The 53rd Ann. Internet. Mtg, EAGE, Expanded Abstracts, 1991: 628-629. http://cn.bing.com/academic/profile?id=2326248607&encoded=0&v=paper_preview&mkt=zh-cn |
| [28] | Virieux J. P-SV-wave propagation in heterogeneous media: Velocity-stress finite-difference method. Geophysics , 1986, 51(4): 889-901. DOI:10.1190/1.1442147 |
| [29] | Levander A R. Fourth-order finite-difference P-SV seismograms. Geophysics , 1988, 53(11): 1425-1436. DOI:10.1190/1.1442422 |
| [30] | Hastings F D, Schneider J B, Broschat S L. Application of the perfectly matched layer (PML) absorbing boundary condition to elastic wave propagation. J. Acoust. Soc. Am. , 1996, 100(5): 3062-3069. |