2011, Vol. 54
2. 中国科学院研究生院,北京 100049
2. Graduate School of Chinese Academy of Sciences, Beijing 100049,China
在我国，目前的油气勘探重点已转移到了西部、西南部地区，其剧烈的地形起伏给地震勘探工作提出了严峻的挑战，传统的基于平缓构造勘探的地震数据采集、处理和解释方法在这类复杂地表地区已不再适用.研究起伏地表地震波场数值模拟，对于认识起伏地表下地震波传播的规律及特征，进而研发基于起伏地表的地震资料处理、采集和解释方法均有非常重要的意义[1~8].
为了研究起伏地表条件下的地震波传播特征，近年来提出了一些有针对性的数值模拟方法:包括有限元法[9~11]，谱元法[12， 13]，伪谱法[14~16]，边界元和边界积分法[17~19]，有限差分法[20~27]等.谱元法和有限元都能自然满足自由表面边界条件，可以使用曲边单元模拟二维自由表面和界面起伏.传统的有限元法计算代价很高；在谱元法中，需要使用比数值频散所要求的小得多的空间离散步长来刻画大曲率的地表面和界面[12]；伪谱法中使用的全局基使模拟结果变得不准确[28]；边界元方法的半解析性质决定了该方法不能适用于地表速度变化较大的情况[29]，而实际情况是，浅部地层受后期地质作用速度变化甚为剧烈.有限差分法是数值模拟中最常用的一种方法(与其他方法相比，有限差分法简单灵活)，用有限差分法求解二阶的弹性动力学方程已广泛应用于波场的数值模拟中[30~34]，但当处理自由边界条件且纵横波速度比较大时往往出现不稳定问题[35， 36].Ilan[37]提出了一种改进的方法，但当处理强非均匀介质时依然不稳定.Vidale和Clayton[38]提出了一种隐式的算法来处理自由边界.为解决不稳定问题，人们探索用一阶速度应力方程来处理自由边界，并采用交错网格来离散化.目前应用于地震波场数值模拟的差分算法大部分是基于交错网格技术的[39~42]，然而交错网格方法处理起伏地表比较困难，因而基于非结构网格的一类方法如谱元法等近年来得到了迅速发展.
最近，Nilsson等[43]提出了一种稳定的、可以广泛应用于任意纵横波速度比、且基于二阶波动方程的显式的处理自由边界的差分方法.在边界上，采用一阶差分算子来离散波动方程中混合导数项(∂x∂z和∂z∂x)中的法向导数(∂z)，而在内部网格点上，用二阶中心差分来离散.正如Nilsson 等[43]所证明的，方程中在边界上的一阶离散代入边界条件方程后精度就变成二阶了，这种处理自由地表边界条件的方法不仅不会降低数值计算的精度，反而使处理容易简便，且结果稳定[36].此后，Appelo和Petersson[44]通过坐标转换，把这种方法推广到曲线坐标系，提出了一种稳定的、可模拟不规则区域各向同性介质中地震波传播的方法.实际的地球介质常常是呈各向异性的[45， 46].这可能源于裂隙或孔隙诱导的各向异性[47~50]，薄各向同性层产生的各向异性[34， 51]以及组成岩石的矿物微粒在排列上有择优取向时形成岩石的固有各向异性[52] 等.本文中将Appelo 和 Petersson提出的方法[44]加以应用并推广以模拟复杂地表非均匀各向异性介质中地震波的传播并对其波场特征进行分析.
2 笛卡尔坐标和曲线坐标的转换对于复杂地表的介质，离散网格边界需与起伏的地表吻合以避免人为的产生边界散射.这种网格被称作贴体网格[53， 54]，且广泛应用于数值模拟中[25， 55~59].贴体网格可以通过由计算空间到物理空间的坐标变换来获得(图 1).通过坐标变换曲线坐标系的q，r映射到物理空间的笛卡尔坐标系，这两个坐标系中垂直轴都取向下方向为正.
贴体网格可分为规则网格和不规则网格两种.规则网格在每个坐标轴方向都有固定的单元数，在二维情况下这些单元一般为四边形.本研究中采用了规则的贴体网格.生成贴体曲线网格的方法主要有代数法、保角变换法及偏微分方程法等[53， 54， 60， 61].其中偏微分方程法因其能方便地控制所生成网格点的疏密，且所生成的网格在边界处有良好的正交性而得到广泛的应用.本文中采用偏微分方程法生成贴体网格.
贴体网格生成之后，笛卡尔坐标系中的网格点与曲线坐标系中的网格点就一一对应，即
|
(1) |
由链锁规则，我们有
|
(2a) |
|
(2b) |
|
(3a) |
|
(3b) |
这里qx表示∂q(x，z)/∂x，qz，rx，rz的意义也类似.这些导数叫做度量导数，把式(2a, b)分别代入式(3a, b)，经过化简，可得
|
(4) |
J是变换的雅克比矩阵，它可以表示为J=xqzr-xrzq，详细的推导过程请参见附录A.值得注意的是，即使映射关系式(1)有解析的形式，度量导数依然通过数值的方式来计算以避免当使用守恒形式的动量方程时系数偏导数引起的假源项[54].本文中所有例子的度量导数均采用二阶差分计算.
3 笛卡尔坐标系和曲线坐标系中的 VTI介质弹性波方程对VTI介质，不考虑外力时，笛卡尔坐标系中的弹性波方程为
|
(5a) |
|
(5b) |
其中，cij(x，z)表示介质的弹性参数；u，v分别是沿x，z轴方向的位移分量；ρ(x，z)为介质密度.
利用关系式(2a, b)，方程(5a, b)可变换为
|
(6) |
|
(7) |
笛卡尔坐标系中的自由边界条件为
|
(8) |
式中，[nx，nz] T 为自由表面的内法线.利用关系式(2a, b)，自由边界条件可变换为
|
(9) |
|
(10) |
式中，用法向度量来表示法向
|
(11) |
为了离散方程(6)~(7)，我们用以下网格点来离散一矩形区域:
|
式中，l，w分别为计算区域的长和宽，hq，hr>0且分别表示沿q，r方向的网格间距.两个波场分量为
|
为简化起见，仿照Appelo 和Petersson的做法引入下列差分算子表示符:
|
(12) |
方程(6)~(7)的右边部分的空间导数包含了四种基本的类型，可以分别按照以下方式来离散:
|
(13) |
式中，ω 表示u或v；a，b，c，d为度量导数和弹性参数的组合；E为求平均算子:
|
(14) |
在边界上(j=1)，采用一阶差分算子来离散波动方程中混合导数项(∂q∂r和∂r∂q)中的法向导数(∂r)，而在内部网格点上(j≥2)，用二阶中心差分来离散，即
|
(15) |
采用(13)式近似计算方程(6)~(7)的空间导数项，为了便于阅读，去掉了角标(i，j)，可得:
|
(16) |
|
(17) |
其中M为
|
(18) |
式中，k，l表示q或r；kx，kz，lx，lz由公式(4)计算.
在时间域我们用二阶精度的中心差分进行离散，完整形式的离散方程形式如下:
|
(19) |
|
(20) |
边界条件方程(9)~(10)可离散为
|
(21) |
|
(22) |
满足边界条件(21)~ (22)的波动方程(19)~(20)的具体的差分格式详见附录B.
在求∂q∂r和∂r∂q项中的法向导数(∂r)时巧妙地运用一阶算子
模拟中的吸收边界条件采用Cerjan 等[62]建议的吸收边界:
|
(23) |
其中I为给定的吸收边界带宽的节点数，i为吸收边界内的节点号，α 是衰减系数，其值与I的大小有关.
6 数值模拟实例 6.1 水平地表模型首先设计了一个均匀VTI介质模型，该模型尺寸为6000m×3000 m.模型介质弹性参数为c11=25.5GPa, c13 =14.0 GPa, c33 =18.4 GPa, c44 =5.6GPa, 密度为2.59g/cm3.模拟计算中，纵横向空间网格大小均为10 m, 时间采样间隔取1ms.垂直方向点力源位于地表(500m, 0m)处，力源子波函数为
|
(24) |
其中，t0=0.5s, 截止频率f0=10Hz.
把数值解与二维兰姆问题的解析解进行对比考查本文方法处理自由地表边界条件的正确性.将 VTI介质自由半空间格林函数和子波函数卷积即可得到VTI介质兰姆问题的解析解[45， 63].从地震记录中分别取炮检距为120m 和990m 的地震道，和解析解对比结果见图 2，数值解和解析解吻合得很好，说明本文采用的数值方法对自由边界条件的实现是正确的.
|
图 2 水平自由表面均匀半空间模型中偏移距分别为120m (a)和990 m(b)的地震道解析解和数值解对比.实线为解析解，虚线为数值解 Fig. 2 Waveform comparisons of the vertical component between analytical and numerical solutions at the receivers with the offsets of 120 and 990 m, respectively, for the homogeneous planar free surface model.Solid and dotted lines are analytical and numerical solutions, respectively. |
图 3(a, b)分别为模拟所得到的水平和垂直分量的地震记录.图中，直达qP 波、qSV 波和瑞雷面波(由于瑞雷波和横波的速度差异比较小，在波场记录中很难区分)清晰可见.由弹性波数值模拟传播1.8s后的波场快照(图 4)可见:各向异性介质中弹性波波前面具有明显的方向性；qSV 波波前初至出现三叉区现象[32， 64， 65]；自由界面产生较强能量的瑞雷面波，并随深度增加迅速衰减；发现连接qSV波和沿地表传播的qP波的首波H.
|
图 3 水平自由表面均勻半空间模型中波场位移x分量(a)和z分量记录(b) Fig. 3 Seismograms for the homogeneous planar free surface model |
|
图 4 二维水平自由表面均勻半空间模型波场位移1.83时工分量(a)和z分量(b)快照 Fig. 4 Snapshots of the wavefield at 1.8 s for the homogeneous planar free surface model |
起伏地表对地震波的传播有重要影响.选用了三个起伏地表模型加以研究.首先是一个边界光滑的小山丘模型，接着我们考虑一个自由表面有半圆形凹陷的半空间模型，对这个模型的研究有助于测试本文的方法处理不光滑地表地形的能力，为了简单起见，假设半空间是均匀的，介质的弹性参数与上例相同，震源位于地表距左侧边界500 m 处.最后是一个地表地形为正弦形变化的不均匀介质模型.三个模型中震源均与水平地表模型的相同.
6.2.1 小山丘地表模型模型大小为6km×3km, 在模型表面的中央有一个小山丘，它的高程可以用高斯函数表示:
|
(25) |
纵横向采样间隔均为10m, 时间采样间隔为0.8ms, 图 5为小山丘模型结构化网格剖分图.图 6(a, b)分别为模拟所得到的水平和垂直分量的地震记录.从图中可以看出，由于自由地表小山丘的影响(与图 3的地震记录比较)，小山丘右侧qP 波和瑞雷面波的振幅均有所减小，在直达qP 波之后观测到由瑞雷面波散射转换成的次生qP 波(RqPf)，同理，在瑞雷面波之前也可以观察到由于qP散射转换的次生瑞雷面波(qPRf)的存在，另外，观察到反射的瑞雷面波(qPRb, RR)和qP波(qPqP，RqPb)，这充分地体现了起伏的小山丘对地震波传播的影响.
|
图 5 小山丘模型结构化网格剖分为清楚起见，网格密度减少为原来的1/3. Fig. 5 The grids in the hill topography model.For clarity, the grids are displayed with a reducing density factor of 3 |
|
图 6 小山丘模型中波场位移x分量(a)和z分量(b)记录 图中qPd, Rd分别表示qP波衍射波和瑞雷面波衍射波；qPRf qPRb分别表示qP波发生散射，转换为Rayleigh面波向前和向后传播； RqPfRqPb分别表示Rayleigh面波发生散射，转换为qP波向前和向后传播；qPqP表示qP波反射波；RR表示Rayleigh波反射波. Fig. 6 Seismograms for the hill topography model Symbols mean the following (qPd) qP wave diffracts to qP wave； (Rd) Rayleigh wave diffracts to Rayleigh wave； (qPRf) qP wave scatters to Rayleigh wave and propagates forward； (qPRb) qP wave scatters to Rayleigh wave and propagates backward； (qPqP) qP wave reflectes to qP wave； (RqPf) Rayleigh wave scatters to qP wave and propagates forward； (RqPb) Rayleigh wave scatters to qP wave and propagates backward； (RR) Rayleigh wave reflectes to Rayleigh wave. |
图 7(a~h)是弹性波数值模拟的垂直分量在不同时刻的波场快照.从快照图 7(a, b)看，此刻波场中存在的波为qP 波，qSV 波，瑞雷面波以及连接 qSV波和沿地表传播的qP波的首波，在2.3s时(图 7c)瑞雷面波刚刚到达小山丘表面，开始产生反射波和转换波.与图 7(a, b)不同的是在快照图 7(e~h)上均可发现有比较明显的反射瑞雷面波的存在，由于在小凸面上存在着波的类型的转换等现象，故反射瑞雷面波的能量大为减弱.同样可发现由小凸面产生的反射qSV 波存在.
|
图 7 小山丘模型中波场位移2分量不同时刻的快照 Fig. 7 Snapshots of the vertical component of the wavefield at different propagation times for the hill topography model |
第一个模型实例的地表起伏是光滑连续的，即地表处处都有连续且有限的斜率.半圆形凹陷模型在凹陷边缘的斜率是无穷大的，因此模拟这样一个极端的起伏模型具有重要意义.半圆凹陷位于模型的中央，半径为150m.
数值模型沿x、z分别含有601、291 个网格点，沿x方向的网格间距均为10m, 垂直方向的网格间距随着深度是变化的，靠近自由表面时间距较小，靠近模型底部时间距较大，最小和最大间距分别为6m和12m, 沿垂向的网格间距平均约为10.3m (图 8)，时间采样间隔为0.8ms.
|
图 8 半圆形凹陷模型结构化网格剖分 为清楚起见，网格密度减少为原来的1/3. Fig. 8 The grids in the semi-circular shape depression topography model For clarity, the grids are displayed with a reducing density factor of 3. |
由于凹陷的存在，qP波和瑞雷面波的振幅在凹陷右侧都明显的减小.另外，在直达qP 波之后也可以观察到由瑞雷面波散射转换成的次生qP 波(RqPf)，同理，在瑞雷面波之前也可以观察到由于 qP散射转换的次生瑞雷面波(qPRf)的存在，还有一些反射的瑞雷面波(qPRb, RR)和qP 波(qPqP，RqPb).凹陷的边缘处(x=2850和x=3150m)，由于地形突变而引起的体波和瑞雷面波发生强烈的散射，这些都可以在地震记录图 9中清楚地观察到.由于面波的波长较小，因此瑞雷面波产生的散射要明显强于体波，这表明这样的陡凹陷模型能明显地阻碍瑞雷面波的传播.图 10为波场位移垂直分量的波场快照，与小山丘模型的波场快照(图 7)相比，我们可以看到由于瑞雷面波在半圆凹陷的边缘就发生散射，因此反射的瑞雷面波波前较小山丘模型中的(反射瑞雷面波波前)超前，而且能量也比后者强，这也表明了陡凹陷模型更能阻碍瑞雷面波的传播.
|
图 9 半圆形凹陷模型中波场位移x分量(a)和z分量(b)记录(图中符号的含义同图 6) Fig. 9 Seismograms for the semi-circular shape depression topography model (The meanings of the symbols are the same as in Fig.6) |
|
图 10 半圆形凹陷模型中波场位移z分量的快照 Fig. 10 Snapshots of the vertical component of the wavefield at different propagation times for the semi-circular shape depression topography model |
最后，我们选取一个地表高程可以用下列正弦函数表述的不均匀介质模型:
|
(26) |
模型分两层，模型长6000 m, 宽3000 m, 上、下层各宽1500m.选取上层介质的弹性参数与前面例子中的(弹性参数)相同，模型下层的弹性参数为:c11=71.8GPa, c13=1.2GPa, c33=53.4GPa, c44=26.1GPa, 密度为2.81g/cm3.图 11是正弦形地表地形模型结构化网格剖分图，纵横向网格大小分别为6和10m.图 12(a~f)分别显示传播0.8s, 1.3，1.8，2.3，3.3s和4.3s后垂直分量波场快照，图 13为地表观测的波场位移记录.从中可见，任意起伏地表使得波场成分复杂:单炮记录中直达波同相轴不再是直线，反射波同相轴变为非双曲线，并显示出与地形起伏强相关的特点；不规则起伏地表成为强散射源，导致地表记录和波场快照上散射波发育；散射波、直达波和反射波混叠使得识别地下的有效信息变得异常困难.
|
图 11 正弦形地表地形模型结构化网格剖分. 为清楚起见，网格密度减少为原来的1/5. Fig. 11 The grids in the sinusoidal topography model. For clarity, the grids are displayed with a reducing density factor of 5. |
|
图 12 正弦形地表地形模型中波场位移z分量的快照 Fig. 12 Snapshots of the vertical component of the wavefield for the sinusoidal topography model |
|
图 13 正弦形地表地形模型中波场位移x分量(a)和z分量(b)记录 Fig. 13 Synthetic seismograms coming from the sinusoidal surface topography model |
本文应用了Appelo 和Petersson的方法，实现了起伏地表下不均匀各向异性介质中地震波场传播的显式有限差分法模拟，获得如下基本认识与结论:
(1) 通过把矩形网格映射到曲线网格上处理了起伏地表问题，采用在边界上修饰的差分算子离散波动方程中含垂向导数的混合导数项，有效地处理了自由地表的问题.
(2) 复杂地形的自由地表模拟实例充分地展现了真实地球表面附近传播的地震波场的复杂.合成地震记录和波场快照显示:地表起伏会使直达波和反射波同相轴弯曲，并显示出与地形起伏强相关的特点；产生散射波、衍射波、多次反射波与转换波等复杂震相；在各向异性介质的波场快照中也观察到 qS波三叉区的存在.
附录A 度量导数把式2(a, b)分别代入式3(a, b)，可得
|
(A1) |
|
(A2) |
经过整理可得
|
(A3) |
|
(A4) |
对比偏导数两边的系数，有
|
(A5) |
|
(A6) |
|
(A7) |
|
(A8) |
求解上述方程可得
|
(A9) |
其中，J为转换的Jacobi矩阵J=xqzr-xrzq.
附录B VTI介质波动方程(曲线坐标系)的差分格式在边界上(j= 1)，波动方程需满足自由地表边界条件，把边界条件方程(21)代入波动方程(19)，可得到自由地表上水平向的位移分量u满足的差分格式:
|
(B1) |
方程中各变量的意义同文章中(19)~(20)式相同.把边界条件方程(22)代入波动方程(20)，可得到自由地表上垂直向的位移分量v满足的差分格式:
|
(B2) |
在内部网格点上(j≥2)，用二阶中心差分来离散，可分别得u、v所满足的差分格式:
|
(B3) |
|
(B4) |
通过观察上面的差分格式(需把M中的qx，qz，rx，rz用(4)式替换)不难发现，在实际计算中，差分格式中系数的hq，hr(曲线坐标系中的网格间距)，同M中的hq，hr抵消掉了，因此差分格式实际是与hq，hr无关的(这只是个虚拟步长，它的值不会影响计算的结果).
致谢本研究工作得到了中国科学院地质与地球物理研究所张中杰研究员的悉心指导和帮助．笔者在此表示衷心的感谢．
| [1] | 董良国, 郭晓玲, 吴晓丰, 等. 起伏地表弹性波传播有限差分法数值模拟. 天然气工业 , 2007, 27(10): 38–41. Dong L G, Guo X L, Wu X F, et al. Finite difference numerical simulation for the elastic wave propagation in rugged topography. Natural Gas Industry (in Chinese) , 2007, 27(10): 38-41. |
| [2] | 孙建国. 复杂地表条件下地球物理场数值模拟方法评述. 世界地质 , 2007, 26(3): 345–362. Sun J G. Methods for numerical modeling of geophysical fields under complex topographica1 conditions: a critical review. Global Geology (in Chinese) , 2007, 26(3): 345-362. |
| [3] | 张华, 李振春, 韩文功. 起伏地表条件下地震波数值模拟方法综述. 勘探地球物理进展 , 2007, 30(5): 334–339. Zhang H, Li Z C, Han W G. Review of seismic wave numerical simulation from irregular topography. Progress in Exploration Geophysics (in Chinese) , 2007, 30(5): 334-339. |
| [4] | 阎世信, 刘怀山, 姚雪根. 山地地球物理勘探技术. 北京: 石油工业出版社, 2000 . Yan S X, Liu H S, Yao X G. Geophysical Exploration Technology in the Mountainous Area (in Chinese). Beijing: Petroleum Industry Press, 2000 . |
| [5] | 邓志文. 复杂山地地震勘探. 北京: 石油工业出版社, 2006 . Deng Z W. Exploration in Complex Mountainous Area (in Chinese). Beijing: Petroleum Industry Press, 2006 . |
| [6] | 张永刚. 复杂介质地震波场模拟分析与应用. 北京: 石油工业出版社, 2007 . Zhang Y G. Seismic Wave Field Simulation Analysis and Application in Complex Media (in Chinese). Beijing: Petroleum Industry Press, 2007 . |
| [7] | 郑鸿明, 吕焕通, 娄兵, 等. 地震勘探近地表异常校正. 北京: 石油工业出版社, 2009 . Zheng H M, Lv H T, Lou B, et al. Near-Surface Anomaly Correction in Seismic Exploration (in Chinese). Beijing: Petroleum Industry Press, 2009 . |
| [8] | 裴正林, 何光明, 谢芳. 复杂地表复杂构造模型的弹性波方程正演模拟. 石油地球物理勘探 , 2010, 45(6): 807–818. Pei Z L, He H M, Xie F. Elastic wave equation forward modeling for complex surface and complex structure model. Oil Geophysical Prospecting (in Chinese) , 2010, 45(6): 807-818. |
| [9] | Huang Z P, Zhang M, Wu W Q, et al. A domain decomposition method for numerical simulation of the elastic wave propagation. Chinese J. Geophys. , 2004, 47(6): 1094-1100. |
| [10] | Rial J A, Saltzman N G, Ling H. Earthquake-induced resonance in sedimentary basins. American Scientist , 1992, 80(6): 566-578. |
| [11] | Toshinawa T, Ohmachi T. Love-wave propagation in a three-dimensional sedimentary basin. Bull. Seism. Soc. Am. , 1992, 82(4): 1661-1667. |
| [12] | Komatitsch D, Tromp J. Introduction to the spectral element method for three-dimensional seismic wave propagation. Geophysical Journal International , 1999, 139(3): 806-822. DOI:10.1046/j.1365-246x.1999.00967.x |
| [13] | Komatitsch D, Vilotte J. The spectral element method: An efficient tool to simulate the seismic response of 2D and 3D geological structures. Bulletin of the Seismological Society of America , 1998, 88(2): 368-392. |
| [14] | Nielsen P, If F, Berg P, et al. Using the pseudospectral technique on curved grids for 2D acoustic forward modelling. Geophysical Prospecting , 1994, 42(4): 321-342. DOI:10.1111/gpr.1994.42.issue-4 |
| [15] | Tessmer E, Kosloff D. 3-D elastic modeling with surface topography by a Chebychev spectral method. Geophysics , 1994, 59(3): 464-473. DOI:10.1190/1.1443608 |
| [16] | Tessmer E, Kosloff D, Behle A. Elastic wave propagation simulation in the presence of surface topography. Geophys. J. Int. , 1992, 108(2): 621-632. DOI:10.1111/gji.1992.108.issue-2 |
| [17] | Bouchon M, Campillo M, Gaffet S. A boundary integral equation-discrete wavenumber representation method to study wave propagation in multilayered media having irregular interfaces. Geophysics , 1989, 54(9): 1134-1140. DOI:10.1190/1.1442748 |
| [18] | Liu E R, Zhang Z J, Yue J H, et al. Boundary integral modelling of elastic wave propagation in multi-layered 2D medium with irregular interfaces. Commun. Comput. Phys. , 2008, 3(1): 52-62. |
| [19] | 符力耘, 牟永光. 弹性波边界元法正演模拟. 地球物理学报 , 1994, 37(4): 521–529. Fu L Y, Mou Y G. Boundary element method for elastic wave forward modeling. Chinese J. Geophys. (Acta Geophysica Sinica) (in Chinese) , 1994, 37(4): 521-529. |
| [20] | Gao H W, Zhang J F. Parallel 3-D simulation of seismic wave propagation in heterogeneous anisotropic media: a grid method approach. Geophysical Journal International , 2006, 165(3): 875-888. DOI:10.1111/gji.2006.165.issue-3 |
| [21] | Hestholm S, Ruud B. 2D finite-difference elastic wave modelling including surface topography. Geophys. Prosp , 1994, 42(5): 371-390. DOI:10.1111/gpr.1994.42.issue-5 |
| [22] | Jih R S, McLaughlin K L, Der Z A. Free-boundary conditions of arbitrary polygonal topography in a two-dimensional explicit elastic finite-difference scheme. Geophysics , 1988, 53(8): 1045-1055. DOI:10.1190/1.1442541 |
| [23] | Lombard B, Piraux J, Gélis C, et al. Free and smooth boundaries in 2-D finite-difference schemes for transient elastic waves. Geophysical Journal International , 2008, 172(1): 252-261. DOI:10.1111/gji.2008.172.issue-1 |
| [24] | Robertsson J O A. A numerical free-surface condition for elastic/viscoelastic finite-difference modeling in the presence of topography. Geophysics , 1996, 61(6): 1921-1934. DOI:10.1190/1.1444107 |
| [25] | Zhang W, Chen X F. Traction image method for irregular free surface boundaries in finite difference seismic wave simulation. Geophysical Journal International , 2006, 167(1): 337-353. DOI:10.1111/gji.2006.167.issue-1 |
| [26] | 张剑锋. 弹性波数值模拟的非规则网格差分法. 地球物理学报 , 1998, 41(S1): 357–366. Zhang J F. Non-orthogonal grid finite-difference method for numerical simulation of elastic wave propagation. Chinese J. Geophys. (in Chinese) , 1998, 41(S1): 357-366. |
| [27] | 张金海, 王卫民, 赵连锋, 等. 傅里叶有限差分法三维波动方程正演模拟. 地球物理学报 , 2007, 50(6): 1854–1862. Zhang J H, Wang W M, Zhao L F, et al. Modeling 3-D scalar waves using the Fourier finite-difference method. Chinese J. Geophys. (in Chinese) , 2007, 50(6): 1854-1862. |
| [28] | Komatitsch D, Barnes C, Tromp J. Simulation of anisotropic wave propagation based upon a spectral element method. Geophysics , 2000, 65(4): 1251-1260. DOI:10.1190/1.1444816 |
| [29] | Bouchon M, Schultz C A, Nafi Toks?z M. A fast implementation of boundary integral equation methods to calculate the propagation of seismic waves in laterally varying layered media. Bulletin of the Seismological Society of America , 1995, 85(6): 1679-1687. |
| [30] | Alterman Z S, Karal F C Jr. Propagation of elastic waves in layered media by finite difference methods. Bulletin of the Seismological Society of America , 1968, 58(1): 367-398. |
| [31] | Alterman Z S, Rotenberg A. Seismic waves in a quarter plane. Bulletin of the Seismological Society of America , 1969, 59(1): 347-368. |
| [32] | Yang D H, Liu E, Zhang Z J, et al. Finite-difference modelling in two-dimensional anisotropic media using a flux-corrected transport technique. Geophys. J. Int. , 2002, 148(2): 320-328. |
| [33] | Zhang Z J, He Q D, Teng J W. Forward modeling of 3-component seismic records in 2-D transversely isotropic media with finite difference method. Can. J. Exp. Geophys. , 1993, 29(1): 51-58. |
| [34] | Zhang Z J, Wang G J, Harris J M. Multi-component wavefield simulation in viscous extensively dilatancy anisotropic media. Physics of the Earth and Planetary Interiors , 1999, 114(1-2): 25-38. DOI:10.1016/S0031-9201(99)00043-6 |
| [35] | Ilan A, Loewenthal D. Instability of finite difference schemes due to boundary conditions in Elastic Media. Geophysical Prospecting , 1976, 24(3): 431-453. DOI:10.1111/gpr.1976.24.issue-3 |
| [36] | Lan H Q, Zhang Z J. Comparative study of the free-surface boundary condition in two-dimensional finite-difference elastic wave field simulation. Journal of Geophysics and Engineering , 2011, 8(2): 275-286. DOI:10.1088/1742-2132/8/2/012 |
| [37] | Ilan A. Stability of finite difference schemes for the problem of elastic wave propagation in a quarter plane. Journal of Computational Physics , 1978, 29(3): 389-403. DOI:10.1016/0021-9991(78)90141-9 |
| [38] | Vidale J E, Clayton R W. A stable free-surface boundary condition for two-dimensional elastic finite-difference wave simulation. Geophysics , 1986, 51(12): 2247. DOI:10.1190/1.1442078 |
| [39] | Levander A R. Fourth-order finite-difference P-SV seismograms. Geophysics , 1988, 53(11): 1425. DOI:10.1190/1.1442422 |
| [40] | 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 |
| [41] | 董良国, 马在田, 曹景忠, 等. 一阶弹性波方程交错网格高阶差分解法. 地球物理学报 , 2000, 43(3): 411–419. Dong L G, Ma Z T, Cao J Z, et al. A staggered-grid high-order difference method of one-order elastic wave equation. Chinese J. Geophys. (in Chinese) , 2000, 43(3): 411-419. |
| [42] | 牟永光, 裴正林. 三维复杂介质地震数值模拟. 北京: 石油工业出版社, 2005 . Mou Y G, Pei Z L. Seismic Numerical for 3D Complex Media (in Chinese). Beijing: Petroleum Industry Press, 2005 . |
| [43] | Nilsson S, Petersson N A, Sjgreen B, et al. Stable difference approximations for the elastic wave equation in second order formulation. SIAM Journal on Numerical Analysis , 2007, 45(5): 1902-1936. DOI:10.1137/060663520 |
| [44] | Appelo D, Petersson N. A stable finite difference method for the elastic wave equation on complex geometries with free surfaces. Communications in Computational Physics , 2009, 5(1): 84-107. |
| [45] | 何樵登, 张中杰. 横向各向同性介质中地震波及其数值模拟. 长春: 吉林大学出版社, 1996 . He Q D, Zhang Z J. Seismic Waves in Transversely Isotropie Medium and Numerical Modelling (in Chinese). Changchun: Jilin University Press, 1996 . |
| [46] | 张中杰. 多分量地震资料的各向异性处理与解释方法. 哈尔滨: 黑龙江教育出版社, 2002 . Zhang Z J. Processing and Interpretation Methods of Multi-Component Seismic Data for Anisotropy (in Chinese). Harbin: Heilongjiang Education Publisher, 2002 . |
| [47] | Crampin S. A review of the effects of anisotropic layering on the propagation of seismic waves. Geophysical Journal International , 1977, 49(1): 9-27. DOI:10.1111/gji.1977.49.issue-1 |
| [48] | Hudson J A. Wave speeds and attenuation of elastic waves in material containing cracks. Geophysical Journal International , 1981, 64(1): 133-150. DOI:10.1111/j.1365-246X.1981.tb02662.x |
| [49] | Liu E R, Crampin S, Queen J H, et al. Velocity and attenuation anisotropy caused by microcracks and microfractures in a multiazimuth reverse VSP. Canadian Journal of Exploration Geophysics , 1993, 29(1): 177-188. |
| [50] | 刘恩儒, 岳建华, 刘彦. 具有离散裂缝空间分布的二维固体中地震波传播的有限差分模拟. 地球物理学报 , 2006, 49(1): 180–188. Liu E R, Yue J H, L Y. Finite difference simulation of seismic wave propagation in 2-D solids with spatial distribution of discrete fractures. Chinese J. Geophys. (in Chinese) , 2006, 49(1): 180-188. |
| [51] | Helbig K. Anisotropy and dispersion in periodically layered media. Geophysics , 1984, 49(4): 364-373. DOI:10.1190/1.1441672 |
| [52] | Forsyth D W. The early structural evolution and anisotropy of the oceanic upper mantle. Geophysical Journal International , 1975, 43(1): 103-162. DOI:10.1111/j.1365-246X.1975.tb00630.x |
| [53] | Hvid S L. Three dimensional algebraic grid generation. 1994 . |
| [54] | Thompson J F, Warsi Z U A, Mastin C W. Numerical Grid Generation: Foundations and Applications. Amsterdam: North-Holland, 1985 . |
| [55] | Fornberg B. The pseudospectral method: Accurate representation of interfaces in elastic wave calculations. Geophysics , 1988, 53(5): 625-637. DOI:10.1190/1.1442497 |
| [56] | 孙章庆, 孙建国, 张东良. 2.5维起伏地表条件下坐标变换法直流电场数值模拟. 吉林大学学报 (地球科学版) , 2010, 40(2): 425–431. Sun Z Q, Sun J G, Zhang D L. 2. 5-D DC electric field numerical modeling including surface topography based on coordinate transformation method. Journal of Jilin University (Earth Science Edition) (in Chinese) , 2010, 40(2): 425-431. |
| [57] | 孙章庆, 孙建国, 张东良. 二维起伏地表条件下坐标变换法直流电场数值模拟. 吉林大学学报 (地球科学版) , 2009, 39(3): 528–534. Sun Z Q, Sun J G, Zhang D L. 2D DC electric field numerical modeling including surface topography using coordinate transformation method. Journal of Jilin Universit y (Earth Science Edition) (in Chinese) , 2009, 39(3): 528-534. |
| [58] | 张东良, 孙建国, 孙章庆. 二维起伏地表直流电场插值法数值模拟. 世界地质 , 2009, 28(2): 242–248. Zhang D L, Sun J G, Sun Z Q. Numerical modeling of DC electrical field interpolation on two dimension undulate topography. Global Geology (in Chinese) , 2009, 28(2): 242-248. |
| [59] | 张东良, 孙建国, 孙章庆. 2维和2. 5维起伏地表直流电法有限差分数值模拟. 地球物理学报 , 2011, 54(1): 234–244. Zhang D L, Sun J G, Sun Z Q. Finite-difference DC electrical field modeling on 2D and 2.5D undulate topography. Chinese J. Geophys. (in Chincsc) (in Chinese) , 2011, 54(1): 234-244. |
| [60] | 蒋丽丽, 孙建国. 基于Poisson方程的曲网格生成技术. 世界地质 , 2008, 27(3): 298–305. Jiang L L, Sun J G. Source terms of elliptic system in grid generation. Global Geology (in Chinese) , 2008, 27(3): 298-305. |
| [61] | 孙建国, 蒋丽丽. 用于起伏地表条件下地球物理场数值模拟的正交曲网格生成技术. 石油地球物理勘探 , 2009, 44(4): 494–500. Sun J G, Jiang L L. Orthogonal curvilinear grid generation technique used for numeric simulation of geophysical fields in undulating surface condition. Oil Geo physical Prospecting (in Chinese) , 2009, 44(4): 494-500. |
| [62] | Cerjan C, Kosloff D, Kosloff R, et al. A nonreflecting boundary condition for discrete acoustic and elastic wave equations. Geophysics , 1985, 50(4): 705-708. DOI:10.1190/1.1441945 |
| [63] | Payton R C. Elastic Wave Propagation in Transversely Isotropic Media. Springer , 1983. |
| [64] | Zheng H S, Zhang Z J, Liu E R. Non-linear seismic wave propagation in anisotropic media using the flux-corrected transport technique. Geophys. J. Int. , 2006, 165(3): 943-956. DOI:10.1111/gji.2006.165.issue-3 |
| [65] | 张中杰, 滕吉文, 万志超. 二维各向异性介质中地震波前面参数方程的建立. 科学通报 , 1994, 39(15): 1399–1402. Zhang Z J, Teng J W, Wan Z C. The establishment of the parametric equations of the seismic wavefronts in two-dimensional anisotropic media. Chinese Science Bulletin (in Chinese) , 1994, 39(15): 1399-1402. |