地球物理学报  2020, Vol. 63 Issue (6): 2375-2385   PDF    
六分量弹性波场数值模拟与分析
张智1,2, 孙丽霞1, 唐国彬1,2, 徐涛3, 王赟4, 王敏玲1,2, 郭希1     
1. 桂林理工大学地球科学学院, 桂林 541004;
2. 广西隐伏金属矿产勘查重点实验室, 桂林 541004;
3. 中国科学院地质与地球物理研究所, 北京 100029;
4. 中国地质大学地球物理与信息技术学院, 北京 100083
摘要:地震波传播过程中,质点的振动不仅包括三个独立的平移部分,还包括三个独立的旋转部分.本文基于一阶速度-应力弹性波方程,采用分裂完全匹配层(SPML)的吸收边界条件,推导了时间导数二阶精度和空间导数高阶精度的交错网格有限差分格式的弹性波速度与应力各分量计算公式,模拟了各向同性介质中均匀模型和层状模型下的六分量波场,并对二维各向同性层状模型下的三个分量地震记录做高分辨率线性拉东变换得到各自的频散能谱.数值模拟分析结果表明:(1)旋转分量的能量要比平动分量弱的多;(2)在平动分量上,面波能量强,频率低,反射P波能量较强,反射S波能量稍弱;在旋转分量上,反射P波能量很弱,S波能量强;(3)与平动分量相比,旋转分量的频散能谱效果更好,能看到基阶和完整的高阶面波,即旋转分量能反映更多的地下介质信息.
关键词: 平动分量      转动分量      交错网格      数值模拟     
Numerical simulation of the six-component elastic-wave field
ZHANG Zhi1,2, SUN LiXia1, TANG GuoBin1,2, XU Tao3, WANG Yun4, WANG MinLing1,2, GUO Xi1     
1. College of Earth Sciences, Guilin University of Technology, Guilin 541004, China;
2. Guangxi Key Laboratory of Hidden Metallic Ore Deposits Exploration, Guilin 541004, China;
3. Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China;
4. School of Geophysics and Information Technology, China University of Geosciences, Beijing 100083, China
Abstract: In the seismic wave propagation, the particle motion includes not only three independent translations, but also three independent rotations. In this paper, we simulate seismic wave propagation in the elastic medium, including body and surface waves, and analyze the wave characteristics of all the six components. Based on the theory of first-order velocity-stress elastic wave equation, we use the absorbing boundary condition of the Splitting Complete Matching Layer (SPML) in the simulation. The calculation formulas of the velocity and stress in the staggered grid finite difference scheme are derived to simulate the six-component wave field in the isotropic medium. And then the three-component dispersive-energy spectra in a two-dimensional isotropic two-layer model are obtained by a high-resolution linear Radon transform. The simulation results show that the energy of rotational components is much weaker than that of translational components. In the translational components, surface waves have stronger energy and lower frequency, while the P-wave energy is stronger and the S-wave energy is weaker. On the contrary, in the rotational components, the P-wave energy is weaker and the S-wave energy is stronger. The dispersion energy spectra of rotational components contain more information on fundamental and higher modes, namely rotational components carry more underground information.
Keywords: Translational component    Rotational component    Staggered grid finite difference scheme    Simulation    
0 引言

地震波在介质中传播时,地面质点运动不仅包括线性运动的三个分量,而且还包括转动的三个分量.传统地震学以观测和模拟地表和建筑物的三分量平移为主,一直忽视了地面旋转运动的测量,然而,转动可以导致地面扭转和结构破坏,有时会成为建筑物破坏的主要原因(李宏男和孙立晔, 2001),这在大型结构和桥梁中更为严重(何超等, 2011).近年来强地震频发,造成重大人员伤亡和经济损失,地震预测研究工作面临空前的机遇与挑战.将旋转与平动结合能更好地得到地球内部结构信息,使无需走时参数的层析成像成为可能,也能为破裂过程提供更详细的线索,并将推动地震预报的深入研究.

近几十年来,在各国研究工作者的共同努力下,关于转动的理论和仪器研究都呈现出可喜的发展势头,对其今后的研究方向也已达成了广泛共识.自19世纪60年代以来,发展了很多获取转动的方法,如利用旋转地震仪直接测量、利用密集台阵测量的平动分量来间接获得、根据经典力学与地震学理论由地震平动记录直接推算;相应地转动地震仪的研究也取得了很大进展.目前被业界认可的且确实能够记录转动分量的转动地震仪主要分为以下两类(Jaroszewicz et al., 2016):

(1) 间接测量转动的转动地震仪.这种转动地震仪基于有限差分原理,并且由一对标准地震计组成.理论研究表明,旋转速度数值等于平动速度旋转值的一半.将两个距离足够近的地震台的水平方向地震记录之差除以水平距离,即得到地震台位置近似的绕铅垂轴的转动分量;同样地,两个地震台的垂向地震记录之差除以水平距离,即得到地震台位置近似的绕水平轴的转动分量.因此,通过记录平移运动,可以获得低灵敏度的旋转运动.这种方法原理简单、效果良好,但是要求地震台间距足够小,从而限制了其实用性和数据精度,且这类转动地震仪适合强震观测.

(2) 直接测量转动的转动地震仪.目前常见的三种直接测量转动的转动地震仪分别是基于电化学、机械和光学的(Ernauer et al., 2012).其中基于MEMS(微电子机械系统)的转动地震仪是目前应用最广泛的高精度地震仪.而利用流体作为惯性质量的电化学装置则是通过测量内部流体运动来记录转动,但是容易受环境的影响.光学转动地震仪则使用一个圆形、三角形或方形的充满气体的腔,且腔中有两个反向激光束的闭路传输,这种转动地震仪通过测量两个激光束的震动频率来记录转动,尽管非常精确,但价格昂贵,一般只适用于空间导航和科学研究.

基于经典力学和地震学的理论,地震旋转效应的重要原因是不对称的剪切应变,转动分量可以通过平动分量来计算(孙士军和陈国兴, 1998; 李宏男, 1991; 王君杰和胡聿贤, 1991),并且有许多初步成熟的方法,如两点差分法、行波法、频率域法(洪钟和崔太成, 2012)、有限差分法等.利用上述方法可以建立转动分量与平动分量之间的关系,并基于地震的平动分量记录得到转动分量.两点差分法简单、清晰、可靠,但需要建立成本高昂的密集台阵,难以广泛应用,因此该方法更适合局部地区.行波法更为可靠,其精度与震中距有关,震中距越大精度越高.频率域法是对行波法的改进,考虑了地震波入射时视速度的频散效应,平动分量不再被视为定型波,该方法可以从平动分量记录中获得高精度的转动分量,并且据频域法的分析结果,对体波(包含P波、SV波和SH波)而言,只有SH型体波能在地表产生旋转加速度,因而具有广泛的实用性.

由于介质模型和边界条件的复杂性,我们通常不能获得理论地震图的解析解,用来评估旋转运动重要性的一个可行方案就是利用数值方法和数学模型来模拟地震现象.近年来,对六分量波场的研究越来越多,为了解六分量波场的特征,需要研究地震波场的传播,这可以看做是求解波动方程(裴正林和牟永光, 2004).常用的求解波动方程的方法有行波法、分离变量法、格林函数法和数值解法.其中数值方法使用更加广泛,包括牛顿法、梯度法、傅里叶变换法、卷积微分算子法、辛几何算法、李群算法、伪谱法(刘鲁波等, 2007)、谱元法(Patera, 1984; 王童奎等, 2007)、有限差分法(王祥春和刘学伟, 2007)和有限元法(刘少林等, 2015),所有上述方法都可以用来模拟地震弹性波场.有限差分法主要有交错网格有限差分法(Virieux, 1984)、旋转交错网格有限差分法(Saenger et al., 2000)等,此外还发展了许多其它有限差分方法(Fornberg, 1998; Podgornova and Lisitsa, 2010).

需要指出的是,在选择数值计算方法时,需要考虑数值方法计算结果的精确性、对自由表面和边界条件的适应性以及是否易于并行.本文基于弹性波的一阶速度应力方程,在稳定条件下,利用交错网格有限差分法和完全匹配层吸收边界,模拟了不同介质中具有二阶差分精度时间导数和高阶差分精度空间导数的六分量地震波场.在交错网格有限差分格式中,给出了弹性波一阶速度应力方程的计算公式,并得到了三个平动位移和三个转动位移,进而比较了六分量地震波场中平动分量和转动分量的差异,有助于处理实际地震资料和加深对六分量波场的认识.

1 理论 1.1 弹性波方程

根据弹性波的一阶速度-应力方程,弹性波方程可以简介的表示为

(1)

(2)

其中,vxvyvz分别是粒子在xyz方向的振动速度,同时物理强度等于0.公式(1)和(2)是一阶速度应力方程的一般形式,弹性波场可以在此基础上进行模拟.

1.2 转动分量

利用泰勒展开法可以求出物体在外力作用下的位移场变化,任意粒子QP点附近的位移可以表示为

(3)

其中:

(4)

(5)

eji是对称张量,rji是反对称张量,分别表示平动和转动分量.其结果是,运动完全可以由三个平动分量和三个转动分量来描述.

1.3 数值方法

Madariaga(1976)首次提出了交错网格有限差分法,在不占用更多存储空间和保证计算复杂度的情况下,它相比传统的有限差分法具有4倍高的精度以及更快的收敛速度.在波动方程中,分别采用二阶差分精度的时间导数和六阶差分精度的空间导数,这有效提高了数值模拟的精度.差分系数am可以表示为(裴正林, 2004):

(6)

1.4 边界条件

当波传播到边界时,在数值模拟中会产生边界反射,干扰真实波场.因此,有必要处理边界并削弱边界反射,即边界条件的定义是不可避免的.对于边界的处理有多种方法,如透明人工边界条件(ABC)、傍轴近似人工边界条件和完全匹配吸收层人工边界条件(PML)(董良国等, 2000).PML包括完全匹配吸收层边界条件的分裂形式(SPML)和非分裂形式(NPML)两种(李佩笑等, 2015),两种方法均能很好地处理边界反射问题.本文左、右和下边界均采用SPML,上边界采用自由边界条件,具体为(王周等, 2012).

(7)

图 1展示了两种边界条件的影响,可以看到SPML很好地吸收了边界反射波,而自由边界条件则产生了面波,表明这两种方法能很好地处理边界问题.

图 1 SPML和自由边界条件下的波场 Fig. 1 Wave fields with SPML and free boundary conditions
2 二维合成数据 2.1 各向同性介质半空间均匀模型

基于二维模型,我们使用了一个简单的均匀半空间模型model 1来比较平动和转动之间的差异.震源子波为中心频率20 Hz的雷克子波,被加载在XZ方向,并作为膨胀源产生P波.炮点在模型中间,埋深2 m,模型尺寸为300 m(宽)×150 m(深).检波器位于地表,间隔1 m.介质参数为VP=3000 m·s-1VS=1800 m·s-1ρ=2600 kg·m-3.地震记录如图 2所示.

图 2 各向同性介质半空间均匀模型model 1的地震记录 (a)X;(b) Z;(c)Ry. Fig. 2 Seismic records of isotropic medium half space homogeneous model 1

直达P波在所有分量中都能看到,能量强且表现为一条直线.SV波在Ry分量中的能量强于平动分量,与P波特征相反.进一步提取每个分量的第80和第130道以观察波形的差异,如图 3所示.

图 3 model 1每个分量的第80和第130道波形 Fig. 3 Waveforms of 80 and 130 channels of each component of model 1

结果显示,第130道由于更接近炮点,相比第80道其初至波到时较短且能量更强.比较不同分量的波形差异,可以发现Z分量的直达P波和面波能量最强,比转动分量强四个数量级.同一道的所有分量具有相同的初至波到时、相同的初始相位和不同的面波频率.为了更好地理解不同分量的波场,我们绘制了三个分量在不同时间的波场快照,如图 4所示.

图 4 model 1的三分量波场快照 Fig. 4 Wave field snapshots of three components of model 1

在三分量的波场快照中,可以清晰地看到直达P波、S波和面波(图 4).对比不同时刻的波场快照,发现P波和S波的波前随时间的推移逐渐增大,到达边界时能量被吸收.同时,对比同一时刻不同分量的波场快照,发现P波能量在Z分量中最强,在转动分量中最弱,而S波能量在转动分量中强于平动分量.

对各向同性均匀半空间介质的数值模拟表明,自由边界条件的存在会导致面波的产生,从而使波场发生变化.埋在地下2 m深处的膨胀源产生P波,P波向上传播并在地表反射,产生转换S波.这表明自由边界对自由表面的模拟有重要影响.

2.2 二维各向同性介质的两层模型

为了进一步理解平动分量和转动分量之间的差异,我们使用了一个二维各向同性层状介质模型model 2,其尺寸为300 m(宽)×150 m(深).震源为中心频率20 Hz的雷克子波,位于模型中间,埋深2 m,被加载在XZ方向上用来模拟膨胀源.检波器间隔为1 m.记录长度为0.2 s,具有二阶时间精度和六阶空间精度.模型有两层,上层厚度为50 m,参数为VP=3000 m·s-1VS=1800 m·s-1ρ=2600 kg·m-3;下层厚度为100 m,参数为VP=4000 m·s-1VS= 2400 m·s-1ρ=2900 kg·m-3.获得的地震记录如图 5所示.

图 5 两层模型model2的地震记录 (a) X;(b) Z;(c) Ry. Fig. 5 Seismic records of two-layer model 2

结果显示三分量均可以清晰地看到直线形状的直达波、面波以及双曲线形状的反射P波、PS波和S波.由于P波速度最快,导致反射P波的走时最短,而S波速度相对较低,导致反射S波走时最长.比较三个分量的差异,可以发现在Z分量上P波的能量明显强于S波,而在Ry分量上则结果相反.此外,面波在三个分量上的能量均很强,并具有明显的频散现象.

通过对比平动和转动分量(图 6),可以发现直达P波在0.03 s到达,其次是面波和反射波,并且Ry分量的能量比XZ分量弱4个数量级,与model 1的结果一致.此外,由于地层界面的存在,波形明显更加复杂了.

图 6 model 2每个分量的第80和第130道波形 Fig. 6 Waveforms of 80 and 130 channels of each component of model 2

为进一步比较平动和转动分量,我们绘制了三个分量的波场快照,如图 7所示.可以看到随着时间的推移,地震波在50 m深度的界面处发生了反射和透射,并产生了反射PP波、PS波、SP波、SS波和透射PP波、PS波、SP波、SS波.另外还可以看到,沿地表传播且能量较强的面波发生了物理频散,并且Ry分量上的频散比在XZ分量上的频散更加明显,这有助于面波和体波的分离.

图 7 两层模型model 2的波场快照 Fig. 7 Wave field snapshots of two-layer model 2

各向同性介质两层模型的数值模拟发现,反射界面的存在对波场有很大影响,可以在界面处引起波的模式转换并产生波的反射和透射,从而使波场更加复杂.此外,地震波场的三分量具有差异性,比如在转动分量中S波的能量要强于P波,而这有助于P波和S波的分离.

2.3 含异常体的两层模型

基于2.2节提到的二维各向同性介质两层模型model 2,我们在第二层增加了一个圆形的低速异常体(model 3),并观察了异常体对平动分量和转动分量的影响,同时比较了低速异常体产生的散射波在三分量中的响应差异.异常体的参数为VP=2000 m·s-1VS=1200 m·s-1ρ=2000 kg·m-3,如图 8所示.

图 8 含异常体的两层模型model 3 Fig. 8 Two-layer with anomaly body model 3

图 9显示,散射波在三个分量上都非常显著,其能量比直达波和面波弱.另外,Z分量上散射波的能量比Ry分量强.随着时间的推移,散射波的能量逐渐降低,如图 10所示.

图 9 model 3的地震记录 (a) X;(b) Z;(c) Ry. Fig. 9 Seismic records of model 3
图 10 model 3每个分量的第80和第130道波形 Fig. 10 Waveforms of 80 and 130 channels of each component of model 3

不同时刻不同分量的波场快照(图 11)表明model 3的三分量波场更加复杂,可以清晰地看到在50 m深度处有一个地层界面,并且在x=150 m、z=100 m处有一个明显的异常体.由于转动分量中的P波能量弱,因此平动分量的波场比转动分量的波场更加复杂.

图 11 model 3的波场快照 Fig. 11 Wave field snapshots of model 3

为了进一步分析异常体造成的差异,我们用model 2和model 3的结果之差来表示异常体的波场,结果如图 12图 13所示.

图 12 第80和第130道的波形差异 Fig. 12 Waveform differences of 80 and 130 channels
图 13 波场差异 Fig. 13 Differences of wave fields

结果显示,低速异常体能够造成大量的散射波,对三分量波场特别是Z分量有很大影响.其中,Z分量散射波场的能量明显强于Ry分量,而在转动分量上低速异常体产生的波场则相对简单.此外,三个分量在界面处的二次反射和透射波也非常显著.

通过对比有异常体和无异常体的两层模型地震波场,我们发现异常体对波场的影响很大,因此记录多分量地震波场对于研究低速异常体的散射波和界面的二次反射和透射波有重要意义.

3 结论

基于一阶弹性波速度应力方程,我们采用交错网格有限差分法模拟了二维各向同性介质中的三分量地震波场.结果表明,在各向同性介质中,平动分量的能量比转动分量能量大得多;在平动分量上,面波能量强,频率低,反射P波能量较强,反射S波能量稍弱;转动分量中,S波和面波的能量明显强于P波,这对波场分离是有利的;与平动分量相比,旋转分量的频散能谱效果更好,能看到基阶和完整的高阶面波,即旋转分量能反映更多的地下介质信息.低速异常体对地震波场影响很大,会产生散射波和二次反射;不同分量对低速异常体的敏感性不同,Z分量受低速异常体影响较大,而转动分量则相反.从数值模拟的角度分析不同地质/地球物理模型对地震旋转运动的影响,以及同一地质/地球物理模型中,旋转运动随时空的变化特征,可以用来指导野外观测旋转运动的地震台站的布设;同时对强地面运动地震学、宽频带地震学、地震工程学、地震物理学、地震仪表设备、地震灾害、地震构造、大地测量学的研究具有重要的指导意义.

References
Dong L G, Ma Z T, Cao J Z, et al. 2000. A staggered-grid High-order difference method of one-order elastic wave equation. Chinese Journal of Geophysics (in Chinese), 43(3): 411-419.
Ernauer F, Wassermann J, Igel H. 2012. Rotational sensors-a comparison of different sensor types. Journal of Seismology, 16(4): 595-602. DOI:10.1007/s10950-012-9286-7
Fornberg B. 1998. A Practical Guide to Pseudospectral Methods. Cambridge: Cambridge University Press.
He C, Luo Q F, Hong Z. 2011. Brief discussion on the study of the seismic rotational components. Journal of Seismological Research (in Chinese), 34(1): 81-87.
Hong Z, Cui T C. 2012. Research of rotation component of ground motiona. Shanxi Architecture (in Chinese), 38(8): 39-40.
Jaroszewicz L R, Kurzych A, Krajewski Z, et al. 2016. Review of the usefulness of various rotational seismometers with laboratory results of Fibre-Optic ones tested for engineering applications. Sensors, 16(12): 2161. DOI:10.3390/s16122161
Li H N. 1991. Study on rotational components of ground motion. Journal of Shenyang Architectural and Civil Engineering Institute (in Chinese), 7(1): 88-93.
Li H N, Sun L Y. 2001. Rotational components of earthquake ground motions derived from surface waves. Earthquake Engineering and Engineering Vibration (in Chinese), 21(1): 15-23.
Li P X, Lin W J, Zhang X M, et al. 2015. Comparisons for regular splitting and non-splitting perfectly matched layer absorbing boundary conditions. Acta Acustica (in Chinese), 40(1): 44-53.
Liu L B, Chen X F, Wang Y B. 2007. Chebyshev pseudospectral method for seismic wavefield modeling. Northwestern Seismological Journal (in Chinese), 29(1): 18-25.
Liu S L, Li X F, Wang W S, et al. 2015. A Lax-Wendroff lumped mass finite element method for seismic simulations. Oil Geophysical Prospecting (in Chinese), 50(5): 905-911, 924.
Madariaga R. 1976. Dynamics of an expanding circular fault. Bulletin of the Seismological Society of America, 66(3): 639-666.
Patera A T. 1984. A spectral element method for fluid dynamics:laminar flow in a channel expansion. Journal of Computational Physics, 54(3): 468-488. DOI:10.1016/0021-9991(84)90128-1
Pei Z L. 2004. Numerical modeling using staggered-grid high-order finite-difference of elastic wave equation on arbitrary relief surface. Oil Geophysical Prospecting (in Chinese), 39(6): 629-634.
Pei Z L, Mu Y G. 2004. Numerical simulation of seismic wave propagation. Progress in Geophysics (in Chinese), 19(4): 933-941.
Podgornova O, Lisitsa V. 2010. Accuracy analysis of finite-difference staggered-grid numerical schemes for elastic-elastic and fluid-elastic interfaces.//80thAnn. Internat Mtg., Soc. Expi. Geophys..Expanded Abstracts, 3087-3091.
Saenger E H, Gold N, Shapiro S A. 2000. Modeling the propagation of elastic waves using a modified finite-difference grid. Wave Motion, 31(1): 77-92. DOI:10.1016/S0165-2125(99)00023-2
Sun S J, Chen G X. 1998. Synthesis method for estimation of rotation components of ground motion. Journal of Seismology (in Chinese), (1): 19-24.
Virieux J. 1984. SH wave propagation in heterogeneous media; Velocity-stress finite-difference method. Geophysics, 49(11): 1933-1942. DOI:10.1190/1.1441605
Wang J J, Hu Y X. 1991. A study on rotational components of surface ground motion. Earthquake Engineering and Engineering Vibration (in Chinese), 11(2): 1-10.
Wang T K, Li R H, Li X F, et al. 2007. Numerical spectral-element modeling for seismic wave propagation in transversely isotropic medium. Progress in Geophysics (in Chinese), 22(3): 778-784.
Wang X C, Liu X W. 2007. Simulation and analysis of seismic wavefield on relief surface by 2-D acoustic wave equation. Oil Geophysical Prospecting (in Chinese), 42(3): 268-276.
Wang Z, Li Z H, Long G H, et al. 2012. Comparison among implementations of free-surface boundary in elastic wave simulation using the finite-difference method. Engineering Mechanics (in Chinese), 29(4): 77-83.
董良国, 马在田, 曹景忠, 等. 2000. 一阶弹性波方程交错网格高阶差分解法. 地球物理学报, 43(3): 411-419. DOI:10.3321/j.issn:0001-5733.2000.03.015
何超, 罗奇峰, 洪钟. 2011. 关于地震动转动分量的研究. 地震研究, 34(1): 81-87. DOI:10.3969/j.issn.1000-0666.2011.01.013
洪钟, 崔太成. 2012. 地震动转动分量研究. 山西建筑, 38(8): 39-40. DOI:10.3969/j.issn.1009-6825.2012.08.023
李宏男. 1991. 关于地震动转动分量的研究. 沈阳建筑工程学院学报, 7(1): 88-93.
李宏男, 孙立晔. 2001. 地震面波产生的地震动转动分量研究. 地震工程与工程振动, 21(1): 15-23. DOI:10.3969/j.issn.1000-1301.2001.01.003
李佩笑, 林伟军, 张秀梅, 等. 2015. 常规分裂和非分裂完全匹配层吸收边界比较研究. 声学学报, 40(1): 44-53.
刘鲁波, 陈晓非, 王彦宾. 2007. 切比雪夫伪谱法模拟地震波场. 西北地震学报, 29(1): 18-25.
刘少林, 李小凡, 汪文帅, 等. 2015. Lax-Wendroff集中质量有限元法地震波场模拟. 石油地球物理勘探, 50(5): 905-911, 924.
裴正林. 2004. 任意起伏地表弹性波方程交错网格高阶有限差分法数值模拟. 石油地球物理勘探, 39(6): 629-634. DOI:10.3321/j.issn:1000-7210.2004.06.002
裴正林, 牟永光. 2004. 地震波传播数值模拟. 地球物理学进展, 19(4): 933-941. DOI:10.3969/j.issn.1004-2903.2004.04.038
孙士军, 陈国兴. 1998. 地面运动转动分量的合成方法. 地震学刊, (1): 19-24.
王君杰, 胡聿贤. 1991. 地震动旋转分量的研究. 地震工程与工程振动, 11(2): 1-10.
王童奎, 李瑞华, 李小凡, 等. 2007. 横向各向同性介质中地震波场谱元法数值模拟. 地球物理学进展, 22(3): 778-784. DOI:10.3969/j.issn.1004-2903.2007.03.018
王祥春, 刘学伟. 2007. 起伏地表二维声波方程地震波场模拟与分析. 石油地球物理勘探, 42(3): 268-276. DOI:10.3321/j.issn:1000-7210.2007.03.007
王周, 李朝晖, 龙桂华, 等. 2012. 求解弹性波有限差分法中自由边界处理方法的对比. 工程力学, 29(4): 77-83.