2017, Vol. 32
2. 中南大学地球科学与信息物理学院, 长沙 410083
3. 河北省国控矿业开发投资有限公司, 石家庄 050035
2. School of Geosciences and Info-Physics, Central South University, Changsha 410083, China
3. Hebei State-owned Assets Holding Mining Development & Investment Co., Ltd, Shijiazhuang 050035, China
近年来国家重视基础设施建设，工程勘察的需求强劲，地震勘察由于资料丰富、设备简单、解决问题能力强，在前期勘察中普遍使用.但现阶段勘察工作多遭遇复杂地形，对野外数据采集和资料处理有较大影响 (杨旭明，2002；邓志文，2006；田文辉，2006)，因而研究起伏地形下的地震波场特征显得十分必要.
起伏地形下地震波传播的解析解难以求出，一般通过数值模拟 (牟永光和裴正林，2005；冯英杰等，2007) 研究其波场特征和传播规律.有限差分法在处理起伏地形时的难点是自由边界问题.自由边界为介质与空气的物性突变面，垂直该界面的应力分量为零 (张华等，2007).目前，起伏地形下的有限差分法处理算法大致分两类：(1) 根据起伏地形条件，改进波动方程以适应地形模型，属全局法；(2) 直接对起伏地形处理，需要丰富的计算方法和技巧，属局部法.Tessmer等 (1992), Tessmer和Kosloff (Tessmer E et al., 1992; Tessmer E and Kosloff D, 1994) 将起伏地形模型通过坐标变换映射到新的水平地表坐标系下进行波动方程的模拟，在时间上用有限差分法，但在空间导数上用Chebyshev变换计算纵向波场，用Fourier变换处理横向波场；Hestholm和Rund (2000)借鉴其思想，通过坐标转换把起伏地表转换为水平面后，实现了在时间和空间上的有限差分法处理.在国内，董良国 (2005)、王祥春 (2007)等也使用纵坐标转换方法，研究了二维、三维起伏地表下的地震波数值模拟，以上均属于全局法.关于局部法，Jih等 (1988)采用直接离散方法处理局部起伏地表，计算量较大，影响了模拟结果的精度.裴正林 (2004)结合广义虚像法和真空法处理起伏地表，在自由界面采用4阶差分，很好地实现了任意起伏地形的弹性波波场模拟.董良国，吴晓丰 (2006)等通过在起伏自由界面上加一层虚拟网格点来实现自由边界条件，并根据不同起伏模型，将整个二维空间离散为24类不同的差分格式实现对起伏地形的处理.真空法也可以用来处理起伏自由边界，但是由于不满足自由边界条件，相对来说精度和稳定性不好.由应力镜像法发展起来的广义虚像法是根据虚像原理，对起伏地形分为不同的节点采用不同的镜像方法来实现自由边界条件.此外常用的局部法还有：横向各向同性介质近似代替自由边界法及声学-弹性介质边界代替自由边界法.汪利民 (2009)将XU的AEA在水平自由地表处理方法推广到起伏地表，并应用到三维情形中.总的来说，直接法处理自由界面不是无条件收敛的，要考虑收敛条件问题；间接法可以灵活处理界面，但要求地形函数连续、光滑、可导，难以满足实际地形的任意性.
由于有限差分计算方法理论比较成熟，算法快捷易实现，本文借鉴前人的思路 (Crase E, 1990；侯安宁，1994；Graves R W, 1996；张剑锋，2000)，选择交错网格有限差分法进行弹性波数值模拟，根据起伏地形，对自由边界条件进行了详细探讨和改进，实现了任意起伏地形下的地震波场模拟.通过地质模型的模拟，分析地震波传播规律及波场特征，旨在为实际工程地震勘探数据处理与资料解释做参考.
1 弹性波方程的高阶交错网格有限差分格式弹性波波动方程是研究地震正演的理论基础，也是研究地震波传播规律的出发点.根据弹性波理论，描述弹性体质点位移、应变和应力之间的关系可由运动平衡微分方程 (Navier方程)、本构方程 (Hook定律) 和几何方程 (Cauchy方程) 来概括.同时，为了避免对介质的弹性常数进行空间微分运算，利用位移与速度的关系，可将方程中的位移分量用速度分量表示.最终可建立起二维弹性各向同性介质中的一阶速度-应力弹性波动方程，公式为
|
(1) |
其中，vx、vz为质点在x，z方向的速度；σxx和σzz为质点在x，z方向的正应力，σxz为质点在xz平面内的剪应力；ρ为弹性介质的密度，t为时间，λ和μ表示介质的拉梅常数.
利用交错网格差分格式，对方程 (1)]进行时间二阶，空间八阶差分，得到差分方程组为
|
(2) |
其中：U、V、R、T和H分别为vx、vz、σxx、σzz和σxz的离散量，上标k为时间离散值，下标i, j为空间离散值，即x，z方向上的网格节点，Δx和Δz为x，z方向的空间采样间隔，Δt为时间步长，Cn为差分系数.
2 自由边界处理方法的改进 2.1 全局法最早由Tessmer (1992)等提出，将起伏地形模型通过坐标变换映射到新的水平坐标系下再进行波动方程模拟, Tessmer通过坐标变换在时间上用有限差分法，在空间导数上，利用Chebyshev变换计算纵向波场对空间的导数，用Fourier变换处理横向波场.该方法不能直接用在有限差分方法中，但是对起伏地形处理有一定指导意义.随后Hestholm、董良国、王祥春等借鉴了前者的思想并加以改进，使之能通过有限差分法求解方程进行数值模拟.全局法的思想就是通过纵向坐标的变换，以消除地形的影响，在狭义上也可认为是纵坐标转换法.基本思路是：
设 (x, z) 为包含起伏地形的曲线坐标系，(ξ, η) 为矩形网格坐标系，如图 (1) 所示.通过映射关系：
|
(3) |
其中，z0(ξ) 为地形函数.且对于 (ξ, η) 矩形网格坐标系的网格边界为ξ=0，ξ=ξmax，η=0和η=ηmax.这样就把起伏地表转换为水平的规则区域，曲线网格也相应地转换为矩形网格.考虑函数f(x, z) 为连续可导函数，则由上式可得：
|
(4) |
于是在原波动方程对x、z的空间导数就可以转换为对ξ和η的空间导数，则在原来 (x, z) 坐标系的波动方程在 (ξ, η) 坐标系下就转换为
|
(5) |
式中U′=(uξ, uη, uξξ, uηη, uξη)T为在 (ξ, η) 坐标系下速度-应力变量的列向量，p(ξ)=
|
图 1 网格映射示意图 Figure 1 Diagram of grids mapping |
相应的自由边界条件也需要进行相应的坐标转换.在起伏地表引入局部坐标 (x′, z′)，其中z′轴为自由地表的法线方向.在局部坐标 (x′, z′) 的边界条件为
|
(6) |
|
(7) |
然后通过 (8) 旋转关系转换到原来 (x, z) 坐标系下.公式为
|
(8) |
其中v和v′分别为 (x, z) 和 (x′, z′) 坐标系下的速度向量，A为旋转矩阵，公式为
|
(9) |
θ为x′与x的夹角.
该法虽然避免了对自由边界的特殊处理，按照常规的水平自由界面处理即可，但是要求地表必须是可导函数，连续光滑，这种条件在实际地形中是不可能的，缺乏实际的意义.另外，地表函数在局部变化过快，会引起吉布斯效应，甚至导致模拟结果的震荡.
2.2 局部法局部法：直接对起伏地表进行处理，不考虑介质的波动方程.主要方法有直接离散法、真空法、广义虚像法及组合方法.
Jih (1988)采用直接离散方法处理起伏自由边界条件，即根据局部地表，划分地表的每个边界点为上倾、下倾、陡坡、缓坡等6种不同情况，并给出了每一种情况具体的差分格式.但是，这种方法在地表处需要通过大量的插值来实现自由边界，这样就影响了模拟的精度.
真空法由于可以根据数据变量和弹性特性自动地识别起伏地形，不需要增加额外的边界计算，因此也适用于起伏地表的情况，但这种方法没有满足自由边界条件的要求，且稳定性不高，仅二阶算子稳定.
广义虚像法是对应力镜像法进行改进，能够适应起伏自由边界的情况.这种方法是根据起伏情况，按照地表附近差分网格上变量及参数分布确定其在网格的角内还是角外，根据其思路，在水平和垂直两个方向向外拓展网格，进行自由边界处理，效果较良好.
裴正林 (2004)采用广义虚像法和真空法相结合的方法处理自由边界，在自由界面采用4阶精度差分，进行了起伏地形的弹性波场模拟，取得了良好的结果.
2.3 改进型真空法针对任意起伏地形自由界面处理问题，本文提出改进型真空法，通过对自由界面上网格质点的放置和弹性参数的设置，来改善算法的精度和稳定性.首先在对自由界面的弹性参数的设置上，参考Mittet和XU对自由界面弹性参数的设置方法，对弹性参数进行平均，即
|
(10) |
其中，ρ、λ和μ为自由界面上的密度和拉梅常数，ρ1和μ1为自由界面下的介质的密度和拉梅常数，这样可以使真空法满足自由边界条件，并增加了稳定性.其次，在自由边界上加一虚拟层 (如图 2)，虚拟层的厚度为半个网格点，设置虚拟层以上为真空，真空区域内所有质点的弹性参数和物理量均设为零.自由界面在原来的位置上移半个网格距离，这样，剪切应力σxz置于自由边界上，正应力置于自由界面下不予考虑，只需要处理剪切应力σxz，令σxz=0，即水平界面、垂直界面上 (A和B点) 和内外角点处 (C和D点) 的切应力均为零.这种与水平自由界面类似的处理方法可使算法更好的满足自由边界条件.
|
图 2 起伏自由表面改进真空法节点网格分布图 Figure 2 Gridding schema of nodes on rugged surface by the improved vacuum method |
对自由界面上速度的处理，在加虚拟层后，水平质点速度vx和垂直质点速度vz也可以准确地置于自由边界上.对速度的计算遵循以下方法：自由界面上质点速度为虚拟层上的速度值与介质内部速度值的平均值.例如E+1点速度：vE+1=1/2(vE0+vE2)，其中，vE0、vE+1与vE2分别为质点E0、E+1与E+2的垂直方向上的速度.
由于是直接法，以折带曲地处理倾斜界面，必须考虑在阶梯状拐角点处的虚假反射.然而有限差分考虑消除数值频散，已经选择了足够精细的网格.研究表明1个最小波长所占网格数大于16时，虚假反射相对于反射波可以忽略不计 (Mitter R, 2002).Hayashi (2001)的研究则表明，虚假反射对模拟精度的影响远小于数值频散.
3 数值模拟结果(1) 倾斜界面模型：一个倾斜界面地质模型.模型大小800 m×800 m，震源主频30 Hz，其模型参数和震源位置见图 3a.
|
图 3 (a) 倾斜界面模型; (b) 粗网格下450 ms时刻波场快照; (c) 细网格下450 ms时刻波场快照 Figure 3 (a) Model with inclined surface; (b) Snapshot at 450 ms with coarse mesh; (c) Snapshot at 450 ms with fine mesh |
|
图 4 (a) 凸起地表模型; (b) 450 ms时刻波场快照图; (c) 单炮弹性波模拟记录 (水平分量); (d) 单炮弹性波模拟记录 (垂直分量) Figure 4 (a) Model with convex surface; (b) Snapshot at 450 ms; (c) Elastic simulation record (horizontal component) for a shot point; (d) Elastic simulation record (vertical component) for the shot point |
从波场快照中可以看到界面处的面波，向外扩展运动的纵波和横波.图 3b计算的参数为空间网格间距1 m，时间采样间隔0.5 ms，从模拟结果看在界面处产生了少许的虚假反射，这与以折带曲的网格呈现阶梯状有关.图 3c是采用网格间距0.25 m，时间采样间隔0.1 ms的模拟结果，可以看出，利用改进的真空法处理自由界面，取得效果不错，但须要精细的网格，这会增加计算成本，对算法的稳定性条件提出更高要求.在综合考虑精度和算法稳定性之后，本文采用网格间距1 m，时间间隔0.2 ms的参数进行计算.
(2) 凸起地表模型：一个地表随机凸起的地质模型.模型大小600 m×600 m，空间采样间隔1 m，时间间隔0.2 ms，震源主频30 Hz.
从波场快照和单炮记录可以清楚地看到P波和S波向外传播，地震波传至凸界面发生多次扰动，产生强散射波和多次反射波，压制能量较弱的散射体波，还可以发生面波和体波之间的转化，波前形态也发生变化，直达波扭曲，走时截断，同相轴不连续且形态变化.同时，侧反射的存在使得波场复杂化.由于凸起模型属于横向不均匀，对水平分量的影响相对较大.
(3) 凹陷地表模型：一个地表任意凹陷的地质模型.其模型参数与凸起地表模型相同.
从波场快照图和单炮模拟记录可以看出，地震波传播到凹陷地形时，面波和体波的相互转化，产生了大量的散射体波与散射面波，散射面波能量强，散射P、S波较弱，凹陷区域的角点处发生绕射.且当地震波从低洼处向上传播时，能量会变的很弱，几乎看不到波.直达波在凹陷区域走时发生截断，同相轴不连续，波形扭曲.
(4) 复杂起伏地表模型：一个任意起伏且存在地下反射界面的地质模型.模型大小为1000 m×1000 m，空间采用间隔1 m，时间采样间隔0.2 ms.震源主频30 Hz.
可见，复杂地形使地震波场变得异常复杂.直达波与地表的形态基本呈负相关性，反射波同相轴扭曲变形，在干扰强烈的地方甚至被淹没.在起伏地表上，可以产生二次或多次震源，生成多次散射波和多次反射波，还可以产生多次转换波及能量很强的侧面波，降低反射波信噪比，由于干扰波的存在，反射波和直达波都变得模糊，对比之下，图 5e中单炮声波模拟记录记录上，反射波、直达波则清晰可辨.
|
图 5 (a) 凹陷地表模型; (b) 500 ms时刻波场快照图; (c) 单炮弹性波模拟记录 (水平分量); (d) 单炮弹性波模拟记录 (垂直分量) Figure 5 (a) Model with concave surface; (b) Snapshot at 500 ms; (c) Elastic simulation record (horizontal component) for a shot point; (d) Elastic simulation record (vertical component) for the shot point |
|
图 6 (a) 复杂起伏地表模型; (b) 750 ms时刻波场快照图; (c) 单炮弹性波模拟记录 (水平分量); (d) 单炮弹性波模拟记录 (垂直分量); (e) 单炮声波模拟记录 Figure 6 (a) Model with complex undulating surface; (b) Snapshot at 750 ms; (c) Elastic simulation record (horizontal component) for a shot point; (d) Elastic simulation record (vertical component) for the shot point.; (e) Acoustic simulation record for the shot point |
(1) 在起伏地形自由界面处理时，采用了本文提出的改进型真空法，实现了任意起伏地面下波场模拟，并且只需要处理剪切应力，相比于真空法，增加了数值模拟的稳定性，同时保证了算法的精度.
(2) 起伏地形使波场复杂，产生强烈的散射体波，散射面波及转换波.同时对波的能量传播有调制作用，使面波能量强，反射波能量弱，同相轴扭曲、截断甚至淹没，信噪比降低.
4.2 建议(1) 对起伏地形采用的以折代曲模拟方法，增加了计算量，建议采用并行算法来提高计算效率.
(2) 起伏地表使弹性波场复杂化，亟需加强弹性波场的资料处理方法研究，尤其是叠前弹性波场成像方法的研究.
致谢 感谢审稿专家提出的修改意见和编辑部的大力支持！| [] | Crase E. 1990. High-order (space and time) finite-difference modeling of the elastic wave equation[C].//1990 SEG Technical Program. Expanded Abstracts, 987-995. |
| [] | Deng Z W. 2006. Complex Mountain Seismic Exploration (in Chinese)[M]. Beijing: Petroleum Industry Press. |
| [] | Dong L G. 2005. Numerical simulation of seismic wave propagation under complex near surface conditions[J]. Progress in Exploration Geophysics (in Chinese), 28(3): 187–194. |
| [] | Feng Y J, Yang C C, Wu P. 2007. The review of the finite-difference elastic wave motion modeling[J]. Progress in Geophysics (in Chinese), 22(2): 487–491. DOI:10.3969/j.issn.1004-2903.2007.02.021 |
| [] | Graves R W. 1996. Simulating seismic wave propagation in 3d elastic media using staggered-grid finite differences[J]. Bulletin of the Seismological Society of America, 86(4): 1091–1106. |
| [] | Hayashi K, Burns D R, Toksöz M N. 2001. Discontinuous-grid finite-difference seismic modeling including surface topography[J]. Bulletin of the Seismological Society of America, 91(6): 1750–1764. DOI:10.1785/0120000024 |
| [] | Hestholm S O, Ruud B. 2000. 2D finite-difference viscoelastic wave modelling including surface topography[J]. Geophysical Prospecting, 48(2): 341–373. DOI:10.1046/j.1365-2478.2000.00185.x |
| [] | Jih R S, McLaughlin K L, Der Z A. 1988. Free-boundary conditions of arbitrary polygonal topography in a two-dimensional explicit elastic finite-difference scheme[J]. Geophysics, 53(8): 1045–1055. DOI:10.1190/1.1442541 |
| [] | Mittet R. 2002. Free-surface boundary conditions for elastic staggered-grid modeling schemes[J]. Geophysics, 67(5): 1616–1623. DOI:10.1190/1.1512752 |
| [] | Moczo P. 1998. Introduction to modeling seismic wave propagation by the finite-difference method[R]. Kyoto:Kyoto University. |
| [] | Mou Y G, Pei Z L. 2005. Seismic Numerical Modeling for 3-D Complex Media (in Chinese)[M]. Beijing: Petroleum Industry Press. |
| [] | Pei Z L. 2004. Numerical modeling using staggered-grid high-order finite-difference of elastic wave equation on arbitrary relief surface[J]. Oil Geophysical Prospecting (in Chinese), 39(6): 629–634. |
| [] | Tessmer E, Kessler D, Kosloff D, et al. 1992. Multi-domain Chebyshev-Fourier method for the solution of the equations of motion of dynamic elasticity[J]. Journal of Computational Physics, 100(2): 355–363. DOI:10.1016/0021-9991(92)90241-P |
| [] | Tessmer E, Kosloff D. 1994. 3-D elastic modeling with surface topography by a Chebychev spectral method[J]. Geophysics, 59(3): 464–473. DOI:10.1190/1.1443608 |
| [] | Tian W H. 2006. Seismic imaging method study based on irregular topography[Ph. D. thesis] (in Chinese). Dongying:China University of Petroleum. |
| [] | Wang L M. 2009. Modeling of rayleigh-wave propagation in three-dimensional media with topography using staggered-grid finite difference scheme[MSc thesis] (in Chinese). Wuhan:China University of Geosciences. |
| [] | Wang X C, Liu X W. 2007. Simulation and analysis of seismic wavefield on relief surface by 2-D acoustic wave equation[J]. Oil Geophysical Prospecting (in Chinese), 42(3): 268–276. |
| [] | Wu X F. 2006. Numerical simulation of seismic wave propagation in the ground surface condition[MSc thesis] (in Chinese). Shanghai:Tongji University. |
| [] | Yang X M. 2002. The research and application of the subsurface seismic wave scatter[Ph. D. thesis] (in Chinese). Chengdu:Chengdu University of Technology. |
| [] | Zhang H, Li Z C, Han W G. 2007. Review of seismic wave numerical simulation from irregular topography[J]. Progress in Exploration Geophysics (in Chinese), 30(5): 334–339. |
| [] | Zhang J F, Liu T L. 2000. Numerical simulation of elastic wave propagation in anisotropic media[J]. Acta Mechanica Solida Sinica (in Chinese), 21(3): 234–242. |
| [] | 邓志文. 2006. 复杂山地地震勘探[M]. 北京: 石油工业出版社. |
| [] | 董良国. 2005. 复杂地表条件下地震波传播数值模拟[J]. 勘探地球物理进展, 28(3): 187–194. |
| [] | 冯英杰, 杨长春, 吴萍. 2007. 地震波有限差分模拟综述[J]. 地球物理学进展, 22(2): 487–491. DOI:10.3969/j.issn.1004-2903.2007.02.021 |
| [] | 侯安宁.1994.各向异性弹性波及其波动方程正反演研究[博士论文].长春:长春地质学院吉林大学. |
| [] | 牟永光, 裴正林. 2005. 三维复杂介质地震数值模拟[M]. 北京: 石油工业出版社. |
| [] | 裴正林. 2004. 任意起伏地表弹性波方程交错网格高阶有限差分法数值模拟[J]. 石油地球物理勘探, 39(6): 629–634. |
| [] | 田文辉. 2006.起伏地表条件下的地震成像方法研究[博士论文].东营:中国石油大学. |
| [] | 汪利民. 2009.三维带地形瑞雷面波交错网格有限差分法正演技术研究[硕士论文].武汉:中国地质大学. |
| [] | 王祥春, 刘学伟. 2007. 起伏地表二维声波方程地震波场模拟与分析[J]. 石油地球物理勘探, 42(3): 268–276. |
| [] | 吴晓丰. 2006.起伏地表条件下地震波传播数值模拟[硕士论文].上海:同济大学. |
| [] | 杨旭明. 2002.近地表地震波散射的研究及应用[博士论文].成都:成都理工大学. |
| [] | 张华, 李振春, 韩文功. 2007. 起伏地表条件下地震波数值模拟方法综述[J]. 勘探地球物理进展, 30(5): 334–339. |
| [] | 张剑锋, 刘铁林. 2000. 各向异性介质中弹性波的数值模拟[J]. 固体力学学报, 21(3): 234–242. |