基于无分裂复频移卷积完全匹配层边界的黏弹介质勒夫波模拟

谢俊法 孙成禹 伍敦仕 乔志浩

1) 中国山东青岛266580中国石油大学(华东)地球科学与技术学院2) 中国兰州730020中国石油勘探开发研究院西北分院



基于无分裂复频移卷积完全匹配层边界的黏弹介质勒夫波模拟

谢俊法1,2)孙成禹1),伍敦仕1)乔志浩1)

1) 中国山东青岛266580中国石油大学(华东)地球科学与技术学院2) 中国兰州730020中国石油勘探开发研究院西北分院

本文建立了无分裂复频移卷积完全匹配层(CFS-CPML)吸收边界条件， 利用交错网格下的高精度有限差分格式对黏弹性介质中的勒夫波场进行了数值模拟； 分析了松弛机制个数对品质因子拟合精度的影响， 验证了CFS-CPML边界条件对大角度掠射波的吸收效果． 数值结果表明： 本文方法所使用的5个松弛机制和空间4阶差分精度， 即可在保证计算效率的前提下满足目前理论研究的需要； 随着品质因子的减小， 频散特征曲线的相速度逐渐向增高的方向偏离理论频散特征曲线的相速度， 且各模式的高频能量也随之减弱． 本文结果可为发展高精度的面波反演方法提供必要的理论依据．

勒夫波 正演模拟 复频移卷积完全匹配层(CFS-CPML) 吸收边界条件 频散特征

引言

勒夫波是在低速层分界面附近传播的一种SH型不均匀平面波， 包含了近地表的构造信息， 因此可以利用勒夫波获取近地表的横波速度、 品质因子等地层信息(Lee， Solomon， 1975， 1978； Xiaetal， 2002， 2013； Winsborrowetal， 2003， 2005； Safanietal， 2005， 2006； Dongetal， 2013)．

地震波在实际地层的传播过程中， 伴随着能量的吸收与衰减， 因而采用线性黏弹性的假设能够更准确地表征地震波尤其是面波在近地表衰减作用下的传播特征． 一般利用品质因子Q来表征介质的黏弹性， 且其在地震频带范围内可近似为常数(Kanamori， Anderson， 1977)． 时间域常Q黏弹性介质中的地震波模拟， 通常采用广义标准线性体(Genera-lized Zener body， 简写为GZB)模型(Carcioneetal， 1988； Carcione， 1993)或广义麦克斯韦尔(generalized Maxwell body define by Emmerich and Korn， 简写为GMB-EK)模型(Emmerich， Korn， 1987； Kristek， Moczo， 2003)实现常Q近似， Moczo和Kristek(2005)已证明这两种模型是等价的． 关于常Q黏弹性介质中体波模拟的研究较多(Liuetal， 1976； 孙成禹等， 2010)， 而针对面波的研究则相对较少. Carcione(1992)曾采用伪谱法对黏弹性半空间情形下的瑞雷波进行正演， 并采取广义标准线性体模型进行常Q拟合， 但其松弛函数遗漏了1/L(L为标准线性体个数或松弛机制个数)因子， 该疏漏后来在其专著中得以修正(Carcione， 2001)． Zhang等(2011)同样采用伪谱法结合广义标准线性体模型进一步分析了两层介质情况下瑞雷面波的频散特征变化， 但是由于其采用的模型较为简单， 高阶模式的变化规律不明显． 高静怀等(2012)基于位移-应力旋转交错网格有限差分法研究了黏弹性介质中瑞雷波的正演， 并与完全弹性情形进行了对比， 但其采用的松弛函数同样也没有1/L因子． 就勒夫波而言， Boore(1970)采用常规网格有限差分法实现了完全弹性横向非均匀介质中的勒夫波数值模拟； Luo等(2010)基于标准交错网格有限差分法研究了完全弹性介质中3种典型地层的勒夫波频散特征． 而关于黏弹性常Q介质中勒夫波数值模拟的研究， 特别是Q值对勒夫波频散特征的影响程度， 尚不多见．

波场数值模拟需要高效的边界反射处理技术． Bérenger(1994)在进行电磁波模拟时提出完全匹配层技术(perfectly matched layer， 简写为PML). 由于PML具有良好的吸收效果， 目前已成为地震波数值模拟中最流行的边界反射处理技术． 早期的PML多采用波场分裂方式， 称为分裂式完全匹配层(Split-PML)(Collino， Tsogka， 2011)， 相对于旁轴近似、 指数衰减等吸收边界处理技术(Clayton， Engquist， 1977； Cerjanetal， 1985)， 其吸收性能有了很大的提高， 但仍然存在某些局限性， 例如， 当震源离PML层很近或接收排列过长时， 不能有效地吸收以掠射角入射的波． 为了克服该问题， 无需波场分裂的PML技术得以发展， 其中比较典型的为复频移卷积完全匹配层技术(complex frequency shifted convolutional perfectly matched layer， 简写为CFS-CPML)． Komatitsch和Martin(2007)比较了完全弹性介质中CFS-CPML与Split-PML在吸收掠射角入射波时的特性， Martin和Komatitsch(2009)进一步比较了黏弹性介质中两种边界的吸收特性． 国内多位研究人员对吸收边界问题也进行了大量的研究工作(张显文等， 2009； 张鲁新等， 2010； 田坤等， 2013； Wangetal， 2013)， 结果表明， 采用CFS-CPML能够有效解决掠射情况下PML吸收效果差的问题， 而且实现过程中无需分裂变量， 提高了计算效率．

鉴于实际地层并非完全弹性， 尤其是近地表介质， 其黏弹性特征较为明显. 为了模拟黏弹性介质中的勒夫波， 本文将首先利用5个松弛机制的广义麦克斯韦尔模型来刻画黏弹性特征， 采取无需分裂的复频移卷积完全匹配层技术处理人工边界反射， 结合时间2阶、 空间4阶精度的标准交错网格有限差分法， 以实现黏弹性常Q介质中的勒夫波正演模拟． 其次， 分析不同Q值对勒夫面波频散特征的影响， 以及速度递增、 夹高速层和夹低速层的复杂地层的勒夫波频散特征， 为利用勒夫波获取近地表信息提供必要的理论依据．

1 勒夫波一阶速度-应力-记忆变量方程推导

在完全弹性条件下， 应力与应变之间是一种瞬时对应关系， 不随时间变化； 而在黏弹性条件下， 应力与应变均为时间的函数， 三维各向同性线性黏弹性介质的本构方程可表示为(Robertssonetal， 1994)：

(1)

式中， “*”表示时域卷积，ψπ(t)和ψμ(t)为应力松弛函数，σ为应力，ε为应变，t为时间， δ为克罗内克函数， 下标遵循爱因斯坦求和约定．

通常基于广义麦克斯韦尔模型(Emmerich， Korn， 1987； Kristek， Moczo， 2003)来实现常Q近似， 基于该模型得到的松弛函数为

(2)

(3)

式中:μu为未松弛模量;ξxyl和ξzyl为对应第l个松弛机制的记忆变量；v为速度分量， 单位为m/s;ρ为密度， 单位为kg/m3．

2 勒夫波黏弹方程的CFS-CPML吸收边界条件

高效的吸收边界是进行波场数值模拟必不可少的条件， 而采用无分裂波场的CFS-CPML吸收边界， 避免了非物理的波场分裂所引起的数值误差， 同时引入了复频移技术， 提高了对高角度掠射波场的吸收效果． 对式(3)进行傅里叶变换， 得到其频域形式：

(4)

式中 “^”表示对应的频域物理量． 经典完全匹配层(perfectly matched layer， 简写为PML)技术实际上是对PML层内的参数进行复坐标拉伸， 即

(5)

(6)

式中sx和sz为复拉伸函数， 具体表示为

(7)

经典PML技术通常采用波场分裂的方式， 对质点的速度分量和应力分量沿坐标轴方向进行分裂(Qinetal， 2009) ． 当震源离PML层很近或接收排列过长时， 不能有效地吸收以掠射角入射的波． 为解决该问题， CFS-CPML技术(Komatitsch， Martin， 2007； Martin， Komatitsch， 2009； 张显文等， 2009； Zhang， Yang， 2010； 张鲁新等， 2010； 田坤等， 2013； 王华等， 2013)引入两个频移因子αx，αz和两个收缩因子kx，kz， 将式(7)推广为

(8)

式中，kx≥1，kz≥1，αx≥0，αz≥0． 频移因子αx和αz的引入有效避免了衰减系数出现奇异值； 收缩因子kx和kz的引入使得地震波进入PML区域后， 传播方向沿该区域的法方向弯曲， 这样高角度入射波能够更加深入PML区域中并逐渐衰减(Zhang， Yang， 2010)． 进一步整理式(8)为

(9)

将式(9)代入式(4)， 经傅里叶反变换， 得到基于CFS-CPML技术实现的一阶速度-应力-记忆变量时域方程如下：

(10)

其中

(11)

为4个内部变量. 沿x方向的衰减系数dx， 收缩因子kx， 频移因子αx分别为(Zhang， Yang， 2010)

(12)

图1 不同物理参数在交错网格 单元中的布局方式 Fig.1 Distribution of different physical parameters within the staggered-grid cells

式中：x表示网格节点到PML层内界面的距离；L为PML层的厚度；d0，k0，α0为常数， 通常取d0=-3cSlogRc/2L，cS为横波速度，Rc为反射系数． 当PML层的厚度为10个网格节点时， 可取Rc=0.1%，k0=2，α0=πf0，f0为震源子波的主频(Collino， Tsogka， 2001)． 沿z方向的处理办法与x方向类似．

由于CFS-CPML技术处理边界不需要进行波场分裂， 降低了存储和求解的复杂程度， 其计算效率比使用Split-PML技术高．

3 标准交错网格离散差分方式

根据图1所示的网格交错方式， 对式(10)进行离散， 即可实现黏弹性SH波的数值模拟． 限于篇幅， 这里仅给出速度分量和内部变量的时间2阶、 空间4阶离散形式：

(13)

(14)

其中

(15)

式(10)中的其它变量， 如σxy，σzy，ξxyl，ξzyl， 可类似导出．

利用显式有限差分进行数值模拟时， 须考虑计算过程的稳定性． 以黏弹介质中的未松弛(高频极限)速度代替弹性速度(田坤， 2014)， 得到

(16)

可以将此作为黏弹性介质稳定性的判断依据． 式中Δx和Δz为横向和纵向的空间网格采样步长， Δt为时间采样间隔，cu为高频极限速度，an为相应的交错网格差分算子系数．

4 数值模拟算例与波场分析

4.1 松弛机制个数的选取

利用式(2)的松弛函数进行常Q近似时， 松弛机制个数对拟合精度具有较大的影响， 松弛机制个数越多， 拟合精度越高， 但相应的计算量也越大； 因此， 在保证拟合精度足够高的前提下， 需要尽量减少松弛机制的个数， 从而减小计算量． 图2给出了分别采用3个和5个松弛机制在1—10 Hz范围内的常Q拟合结果. 可以看出， 采用3个松弛机制的拟合曲线与Q=20的曲线相差较大， 而采用5个松弛机制的拟合曲线则与Q=20的曲线吻合得较好．

图2 松弛机制个数与常Q拟合精度的关系

number and constant Q fitting precision

图3给出了1—10 Hz范围内， 采用5个松弛机制进行常Q拟合的数值相速度与使用Futterman(1962)公式计算所得理论相速度的对比结果， 其中相速度的变化范围为2560—2680 m/s， 其最大相对误差为0.15%， 这充分说明采用5个松弛机制得到的拟合结果具有较高的精度， 故本文的正演模拟中选取5个松弛机制．

图4 均匀各向同性线性黏弹性介质模型示意图Fig.4 Schematic diagram of the homogeneous isotropic linear viscoelastic medium

4.2 边界条件的吸收效果对比

图5 基于黏弹性介质模型采用Split-PML吸收边界所获取的0.14 s (上)，0.3 s (中)， 0.54 s (下)时刻的垂直分量波场快照

图7给出了采用Split-PML和CFS-CPML两种吸收边界获得的炮记录对比， 两侧各40道为吸收边界， 考虑到边界反射较弱， 制图时对图7a和图7b均加了4%的增益． 可以看到： 图7a中炮检距较大的地震道存在顶部边界的反射(箭头所示区域)， 而炮检距较小的地震道则没有， 这是由于大炮检距处波入射到匹配层的入射角过大， 匹配层吸收不完全所致； 图7b中没有顶部边界的反射， 说明CFS-CPML吸收边界对大角度掠射到边界的地震波具有较好的吸收效果．

图6 基于黏弹介质模型采用CFS-CPML吸收边界所获取的0.14 s (上)，0.3 s (中)， 0.54 s (下)时刻的垂直分量波场快照

图7 炮记录的对比(a) 采用Split-PML技术， 箭头所指为顶部边界的反射； (b) 采用CFS-CPML技术

将图4所示的介质模型扩大成6000 m×6000 m， 即网格点数为3000×3000， 激发点设在(3000 m， 3000 m)点处， 接收线的深度设为3002 m． 在这种条件下， 当计算时间tcal≤0.54 s时， 直达波未到达边界， 属于无边界反射记录， 本文将该记录作为评价边界吸收效果的标准． 选取图7a和图7b中炮检距为零的地震数据进行波形显示， 如图8a所示， 可以看出， 采用Split-PML和CFS-CPML吸收边界的波形与无边界反射记录非常一致． 图8a下部放大图中所显示的0.05—0.1 s内振幅为10-3数量级， 存在边界反射的波形与无边界反射的波形几乎一致， 说明地震波以小角度入射到边界时， 采用Split-PML和CFS-CPML吸收边界均有较好的吸收效果． 选取图7a和图7b中炮检距为400 m的地震数据进行波形显示， 如图8b所示. 从图8b下图的直达波尾部放大图可以看出， 采用CFS-CPML吸收边界的波形更接近无边界反射的波形， 说明对于大角度掠射的地震波， 采用CFS-CPML吸收边界的处理效果较好．

图8 单道波形的对比

4.3 波场模拟与分析

4.3.1 双层介质模型的正演模拟与分析

由于勒夫波不存在于均匀半空间中， 其正演模拟需采用多层介质模型． 建立最简单的两层介质模型， 参数见表1. 由于密度对勒夫波的影响较小， 所以模型中两层介质的密度均设为2000 kg/m3． 整个模型由801×200个网格点构成， 网格大小为Δx=Δz=0.5 m， 时间步长为Δt=0.5 ms， CFS-CPML层由10个网格点构成； 震源函数采用主频为f0=40 Hz的雷克子波， 反射系数Rc=0.1%， 收缩因子k0=2， 频移因子α0=πf0， 采用5个松弛机制在1—100 Hz范围内进行常Q拟合． 鉴于Q>150时， 黏弹性衰减效应很小， 因此当地层Q≥200时， 该地层可视为完全弹性介质． 保持底层半空间介质为完全弹性， 依次使第一层的横波品质因子QS为200， 40， 20和10， 采用时间2阶、 空间4阶的有限差分法进行正演得到4个勒夫波记录， 通过倾斜叠加算法(Parketal， 1999)计算其所对应的频散特征(图9a--d)， 并将其与理论频散特征进行对比以考察表层Q值变化对勒夫波频散特征的影响．

表1 两层介质模型参数

当第一层介质的横波品质因子QS值为200时， 两层介质均为完全弹性， 其勒夫波地震记录及相应的频散特征如图9a所示． 从勒夫波记录中可以看出CFS-CPML边界的吸收效果较好． 从图9a给出的频散特征(彩色区域)可以看到： 该记录中包含基阶、 第一高阶、 第二高阶和第三高阶模式的勒夫波； 对每个模式而言， 随着频率的增大， 相速度逐渐减小并趋于表层横波速度300 m/s； 随着频率的减小， 相速度趋于底层半空间横波速度600 m/s． 图9a中的白色曲线为根据矩阵传递算法(Haskell， 1953)正演得到的理论勒夫波频散曲线， 很容易看到， 各模式的频散曲线与其吻合得较好， 说明正演算法能够准确地模拟该地层条件下的勒夫波． 此外， 理论上已证明， 在两层完全弹性介质的情况下， 勒夫波高阶模式存在截止频率， 即当频率低于截止频率时不存在高阶模式(Aki， Richards， 2002)， 这在图9a中也得以体现， 即当频率低于20 Hz时只存在基阶模式， 不存在高阶模式．

图9b给出了当第一层介质QS值为40时的勒夫波地震记录及其相应频散特征． 可以看出， 4个模式的高频成分损耗十分明显， 从约80 Hz降低到近似70 Hz， 这是由于吸收衰减作用所致． 此外， 由于黏弹性所导致的固有频散， 各模式频散曲线的最低相速度均高于完全弹性情况下的理论相速度． 当第一层的QS值进一步减小到20(图9c)甚至10(图9d)时， 从勒夫波记录中可以看出高频成分较完全弹性情况下的地震记录损失较严重， 特别是第601道之后的远偏移距部分更为严重． 从图9c频散图中也可以看到， 此时第三高阶模式高频能量已近乎完全衰减， 第二高阶模式的能量也仅集中在40—60 Hz之间的狭窄频带内， 基阶和第一高阶模式的最高频率降低至60 Hz附近． 若表层介质QS值减小到10(图9d)， 第三高阶模式高频能量则完全衰减， 第二高阶模式的最高频率降低至50 Hz左右， 基阶和第一高阶模式的最高频率也降低， 并且相速度在高频时趋于340 m/s， 与完全弹性时的相速度300 m/s的相对偏差可达13.3%．

图9 双层介质模型中的勒夫波记录(左)及其频散特征图(右)频散图中的白色曲线为根据矩阵传递算法(Haskell， 1953)正演所得的理论勒夫波频散曲线(a) 上层QS=200； (b) 上层QS=40； (c) 上层QS=20; (d) 上层QS=10

综上可见， 勒夫波频散曲线中各模式的最低相速度随第一层介质QS值的减小而增大． 然而， 实际地层并不是完全弹性的， 近地表介质尤其如此， 因此， 在利用面波获取近地表速度时， 第一层介质的速度受频散曲线最低相速度的影响较大， 导致利用面波所获取的第一层速度与真实速度存在一定的误差， 且QS值越小， 该误差越大．

4.3.2 4层介质模型的正演模拟与分析

近地表介质复杂多变， 垂向上的速度并非总是递增， 且近地表压实程度较低， 速度和Q值通常较小． 为便于对比， 研究多层黏弹介质的勒夫波特征时， 保持底层半空间介质为完全弹性， 对3种黏弹介质进行组合， 构成4层介质模型， 具体参数如表2所示． 横波速度为400， 600， 800 m/s和1000 m/s时， 横波品质因子QS分别为10， 20， 50， 200． 除模型参数设置不同外， 4层介质模型的数值模拟参数与表1所示的双层介质模型一致． 从图9a可知弹性介质的勒夫波频散曲线与理论曲线吻合， 因此仅分析黏弹介质的勒夫波频散曲线特征， 弹性介质的频散曲线特征可根据理论曲线作出判断． 下面分速度递增模型、 夹高速层模型和夹低速层模型等3种情况逐一讨论．

表2 4层介质模型参数

图10a--c分别给出了速度递增模型、 夹高速层模型和夹低速层模型的勒夫波记录及其相应的频散曲线， 其中， 彩色区域是从勒夫波记录中采用倾斜叠加算法所提取的频散曲线， 白色线条是理论频散曲线． 可以看出： 速度递增模型和夹高速层模型的频散曲线都只存在基阶模式和第一高阶模式， 基阶模式的能量最强， 且随着频率的增加， 每个模式的频散曲线与理论频散曲线的吻合程度逐渐变差； 与速度递增模型的频散曲线(图10a)相比， 夹高速层模型频散曲线(图10b)的基阶模式弯曲程度更小， 第一高阶模式的能量更弱； 夹低速层模型频散曲线(图10c)尚存在第二高阶模式， 但其能量较弱， 总体频率范围和各模式的能量较速度递增模型增加， 且基阶模式的频率范围更小， 与理论频散曲线的吻合程度更高， 此外， 第一高阶模式的频率范围更大．

由于黏弹性介质对高频能量的吸收较其对低频能量的吸收强， 因此频散曲线中各模式高频部分与理论曲线的吻合程度相对较差． 勒夫波是由一次波和多次反射的SH波相长干涉产生的， 基于几种介质不同组合形成的模型得到的勒夫波频散曲线， 其各阶模式存在差异， 即速度递增模型、 夹高速层模型和夹低速层模型的基阶模式和第一高阶在能量分布、 频率范围和弯曲程度上都存在差异． 反之， 勒夫波包含了地下构造信息， 其频散曲线中各模式的能量分布、 频率范围和弯曲程度是地下构造信息的一种反映， 采用适当的方法能够从勒夫波中获取地下的构造信息， 目前常用的做法是拾取频散特征曲线， 再结合最小二乘或者遗传算法等优化方法来反演得到地下横波速度等底层信息(Xiaetal， 2002， 2013； Safanietal， 2005， 2006)．

图10 4层介质中勒夫波记录(左)及其频散特征图(右)频散图中的白色曲线为根据矩阵传递算法(Haskell， 1953)正演所得的理论勒夫波频散曲线(a) 速度递增模型； (b) 夹高速层模型； (c) 夹低速层模型

5 讨论与结论

本文以无分裂复频移卷积完全匹配层为吸收边界, 利用高精度交错网格有限差分， 推导了黏弹性介质波动方程的一阶速度-应力-记忆变量方程， 并给出了数值解法， 实现了黏弹性勒夫波的数值模拟. 通过模型试算， 得到结论如下：

1) 相对于传统的分列式完全匹配层边界条件， CFS-CPML边界条件可以更加有效地吸收掠射到边界的地震波， 且无需进行波场分裂， 从而提高了计算效率．

2) 本文试验表明， 使用5个松弛机制和4阶差分精度， 能够在保证计算效率的前提下满足目前理论研究的需要．

3) 介质的黏弹性使频散特征曲线的相速度高于理论频散特征曲线的相速度， 各模式的高频成分能量减弱， 且Q值越小， 频散特征曲线与理论频散特征曲线的偏差越大， 各模式高频成分的能量越弱．

实际的近地表地层具有较强的黏弹性， 本文发展了基于无分裂复频移卷积完全匹配层边界的黏弹性勒夫波数值模拟方法， 模拟的频散特征更加接近真实地层情况， 可为面波反演提供理论依据．

高静怀， 何洋洋， 马逸尘. 2012. 黏弹性与弹性介质中Rayleigh面波特性对比研究[J]. 地球物理学报， 55(1): 207--218.

Gao J H， He Y Y， Ma Y C. 2012. Comparison of the Rayleigh wave in elastic and viscoelastic media[J].ChineseJournalofGeophysics， 55(1): 207--218 (in Chinese).

孙成禹， 肖云飞， 印兴耀， 彭洪超. 2010. 黏弹介质波动方程有限差分解的稳定性研究[J]. 地震学报， 32(2): 147--156.

Sun C Y， Xiao Y F， Yin X Y， Peng H C. 2010. Study on the stability of finite difference solution of visco-elastic wave equations[J].ActaSeismologicaSinica， 32(2): 147--156 (in Chinese).

田坤. 2014. 黏性介质正演及逆时偏移成像方法研究[D]. 青岛: 中国石油大学(华东): 35--39.

Tian K. 2014.StudyonMethodofForwardModelingandReverseTimeMigrationinViscousMedia[D]. Qingdao: China University of Petroleum: 35--39 (in Chinese).

田坤， 黄建平， 李振春， 李娜， 孔雪， 李庆洋. 2013. 实现复频移完全匹配层吸收边界条件的递推卷积方法[J]. 吉林大学学报： 地球科学版， 43(3): 1022--1032.

Tian K， Huang J P， Li Z C， Li N， Kong X， Li Q Y. 2013. Recursive convolution method for implementing complex frequency-shifted PML absorbing boundary condition[J].JournalofJilinUniversity：EarthScienceEdition， 43(3): 1022--1032 (in Chinese).

张鲁新， 符力耘， 裴正林. 2010. 不分裂卷积完全匹配层与旋转交错网格有限差分在孔隙弹性介质模拟中的应用[J]. 地球物理学报， 53(10): 2470--2483.

Zhang L X， Fu L Y， Pei Z L. 2010. Finite difference modeling of Biot’s poroelastic equations with unsplit convolutional PML and rotated staggered grid[J].ChineseJournalofGeophysics， 53(10): 2470--2483 (in Chinese).

张显文， 韩立国， 黄玲， 李翔， 刘前坤. 2009. 基于递归积分的复频移PML弹性波方程交错网格高阶差分法[J]. 地球物理学报， 52(7): 1800--1807.

Zhang X W， Han L G， Huang L， Li X， Liu Q K. 2009. A staggered-grid high-order difference method of complex frequency-shifted PML based on recursive integration for elastic wave equation[J].ChineseJournalofGeophysics， 52(7): 1800--1807 (in Chinese).

Aki K， Richards P G. 2002.QuantitativeSeismology[M]. 2nd ed. Sausalito: University Science Books: 251--254.

Bérenger J P. 1994. A perfectly matched layer for the absorption of electromagnetic waves[J].JComputPhys， 114(2): 185--200.

Boore D M. 1970. Love waves in nonuniform wave guides: Finite difference calculations[J].JGeophysRes， 75(8): 1512--1527.

Carcione J M， Kosloff D， Kosloff R. 1988. Viscoacoustic wave propagation simulation in the earth[J].Geophysics， 53(6): 769--777.

Carcione J M. 1992. Modeling anelastic singular surface waves in the earth[J].Geophysics， 57(6): 781--792.

Carcione J M. 1993. Seismic modeling in viscoelastic media[J].Geophysics， 58(1): 110--120.

Carcione J M. 2001.WaveFieldsinRealMedia：TheoryandNumericalSimulationofWavePropagationinAnisotropic，AnelasticandPorousMedia[M]. Oxford: Pergamon Press: 45--78.

Cerjan C， Kosloff D， Kosloff R， Reshef M. 1985. A nonreflecting boundary condition for discrete acoustic and elastic wave equations[J].Geophysics， 50(4): 705--708.

Clayton R， Engquist B. 1977. Absorbing boundary conditions for acoustic and elastic wave equations[J].BullSeismolSocAm， 67(6): 1529--1540.

Collino F， Tsogka C. 2001. Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[J].Geophysics， 66(1): 294--307.

Dong H F， Ke G P， Duffaut K. 2013. Seabed shear-wave velocity estimation from dispersion curves of Love waves[J].ActaAcustunitedAc， 99(6): 890--897.

Emmerich H， Korn M. 1987. Incorporation of attenuation into time-domain computations of seismic wave fields[J].Geophysics， 52(9): 1252--1264.

Futterman W I. 1962. Dispersive body waves[J].JGeophysRes， 67(13): 5279--5291.

Haskell N A. 1953. The dispersion of surface waves on multilayered media[J].BullSeismolSocAm， 43(1): 17--34.

Kanamori H， Anderson D L. 1977. Importance of physical dispersion in surface wave and free oscillation problems: Review[J].RevGeophysSpacePhys， 15(1): 105--112.

Komatitsch D， Martin R. 2007. An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation[J].Geophysics， 72(5): SM155--SM167.

Kristek J， Moczo P. 2003. Seismic-wave propagation in viscoelastic media with material discontinuities: A 3D fourth-order staggered-grid finite-difference modeling[J].BullSeismolSocAm， 93(5): 2273--2280.

Lee W B， Solomon S C. 1975. Inversion schemes for surface wave attenuation andQin the crust and the mantle[J].GeophysJRastrSoc， 43(1): 47--71.

Lee W B， Solomon S C. 1978. Simultaneous inversion of surface wave phase velocity and attenuation: Love waves in western North America[J].JGeophysRes， 83(B7): 3389--3400.

Liu H P， Anderson D L， Kanamori H. 1976. Velocity dispersion due to anelasticity: Implications for seismology and mantle composition[J].GeophysJInt， 47(1): 41--58.

Luo Y H， Xia J H， Xu Y X， Zeng C， Liu J P. 2010. Finite-difference modeling and dispersion analysis of high-frequency Love waves for near-surface applications[J].PureApplGeophys， 167(12): 1525--1536.

Martin R， Komatitsch D. 2009. An unsplit convolutional perfectly matched layer technique improved at grazing incidence for the viscoelastic wave equation[J].GeophysJInt， 179(1): 333--344.

Moczo P， Kristek J. 2005. On the rheological models used for time-domain methods of seismic wave propagation[J].GeophysResLett， 32(1): L01306.

Park C B， Miller R D， Xia J H. 1999. Multichannel analysis of surface waves (MASW)[J].Geophysics， 64(3): 800--808.

Qin Z， Lu M H， Zheng X D， Yao Y， Zhang C， Song J Y. 2009. The implementation of an improved NPML absorbing boundary condition in elastic wave modeling[J].ApplGeophys， 6(2): 113--121.

Robertsson J O A， Blanch J O， Symes W W. 1994. Viscoelastic finite-difference modeling[J].Geophysics， 59(9): 1444--1456.

Safani J， O′Neill A， Matsuoka T， Sanada Y. 2005. Applications of Love wave dispersion for improved shear-wave velocity imaging[J].JEnvironEngGeophys， 10(2): 135--150.

Safani J， O′Neill A， Matsuoka T. 2006. Full SH-wavefield modeling and multiple-mode Love wave inversion[J].ExplGeophys， 37(4): 307--321.

Wang H， Tao G， Shang X F， Fang X D， Burns D R. 2013. Stability of finite difference numerical simulations of acoustic logging-while-drilling with different perfectly matched layer schemes[J].ApplGeophys， 10(4): 384--396.

Winsborrow G， Huws D G， Muyzert E. 2003. Acquisition and inversion of Love wave data to measure the lateral variability of geo-acoustic properties of marine sediments[J].JApplGeophys， 54(1): 71--84.

Winsborrow G， Huws D G， Muyzert E. 2005. The estimation of shear-wave statics usinginsitumeasurements in marine near-surface sediments[J].GeophysProsp， 53(4): 557--577.

Xia J H， Miller R D， Park C B， Tian G. 2002. DeterminingQof near-surface materials from Rayleigh waves[J].JApplGeophys， 51(2/3/4): 121--129.

Xia J H， Yin X F， Xu Y X. 2013. Feasibility of determiningQof near-surface materials from Love waves[J].JApplGeophys， 95(8): 47--52.

Zhang K， Luo Y H， Xia J H， Chen C. 2011. Pseudospectral modeling and dispersion analysis of Rayleigh waves in viscoelastic media[J].SoilDynEarthqEng， 31(10): 1332--1337.

Zhang W， Yang S. 2010. Unsplit complex frequency-shifted PML implementation using auxiliary differential equations for seismic wave modeling[J].Geophysics， 75(4): T141--T154.

Love wave modeling in viscoelastic media with unsplit CFS-CPML conditions

Xie Junfa1,2)Sun Chengyu1),Wu Dunshi1)Qiao Zhihao1)

1)SchoolofGeosciences，ChinaUniversityofPetroleum，ShandongQingdao266580，China2)ResearchInstituteofPetroleumExploration&Development-Northwest,PetroChina,Lanzhou730020,China

Numerical simulation is an important way to study the characteristics of seismic wavefield. In this paper， we established the unsplit complex frequency shifted convolutional perfectly matched layer (CFS-CPML) absorbing boundary condition， and staggered grid and high precision finite difference were used to simulate Love wavefield in viscoelastic medium. And then we analyzed the influence of the number of relaxation mechanisms on theQvalue fitting and verified the absorbing effect of CFS-CPML boundary conditions on the waves with large incident angles. The numerical results show that the requirements of the present theoretical wavefield study can be met with higher calculation efficiency if five relaxation mechanisms and finite difference with spatial four-order accuracy are selected for this method. At the same time， the phase velocity of dispersion curve gradually deviates from that of the theoretical curve with the decrease of quality factor， and the high frequency energy of each pattern also weakens. The study provided certain theoretical basis for high precision surface wave inversion method.

Love wave; forward modeling; complex frequency shifted convolutional perfectly matched layer (CFS-CPML); absorbing boundary condition; dispersion characteristic

谢俊法， 孙成禹， 伍敦仕， 乔志浩. 2016. 基于无分裂复频移卷积完全匹配层边界的黏弹介质勒夫波模拟. 地震学报, 38(2): 244--258. doi:10.11939/jass.2016.02.009.

Xie J F， Sun C Y， Wu D S， Qiao Z H. 2016. Love wave modeling in viscoelastic media with unsplit CFS-CPML conditions.ActaSeismologicaSinica, 38(2): 244--258. doi:10.11939/jass.2016.02.009.

国家自然科学基金(41374123)资助.

2015-08-13收到初稿， 2015-12-07决定采用修改稿.

e-mail: suncy@upc.edu.cn

10.11939/jass.2016.02.009

P315.3+1

A