求解声波方程的径向基函数无网格法_蒋婵君.pdf

第 44 卷 第 5 期 煤田地质与勘探 Vol. 44 No.5 2016 年 10 月 COAL GEOLOGY EXPLORATION Oct. 2016 收稿日期 2016-01-12 基金项目 国家自然科学基金项目（41174077）；大陆构造与动力学国家重点实验室开放基金项目（k201505） Foundation itemNational Natural Sciences Foundation of China（41174077）；Open Foundation of State Key Laboratory of Continental Tectonics and Dynamics of China（k201505） 作者简介 蒋婵君（1984），女，博士研究生，讲师，从事地震波场正演与面波层析成像工作. E-mail554083077 通讯作者王有学（1961），男，博士，教授，从事地球物理与地球内部结构及动力学工作. E-mailuxue.wang 引用格式 蒋婵君，王有学，李川，等. 求解声波方程的径向基函数无网格法［J］. 煤田地质与勘探，2016，44（5）136-141. JIANG Chanjun，WANG Youxue，LI Chuan，et al. Numerical solution of the acoustic wave equation using meshless with radial basis functions［J］. Coal Geology Exploration，2016，44（5）136-141. 文章编号 1001-1986（2016）05-0136-06 求解声波方程的径向基函数无网格法 蒋婵君 1，2，王有学1，2，李 川1，曾高福3，贠 鹏1 （1. 桂林理工大学地球科学学院，广西 桂林 541006； 2. 大陆构造与动力学国家重点实验室，北京 100037； 3. 中国有色桂林矿产地质研究院有限公司，广西 桂林 541004） 摘要 地震勘探广泛应用于油气、煤田勘探。地震波场数值模拟是整个地震勘探数据处理技术的 基石。将径向基函数（RBF）引入地震声波波场数值模拟中，在空间上用径向基函数无网格法来构造 二阶导数，而在时间上采用简单的二阶差分公式，并重点讨论了形状参数 c 对该方法精度的影响， 总结 c 经验取值范围为 2～4 倍平均数据点间距。设计不同模型，利用径向基函数无网格法进行声 波波场模拟，并与空间四阶时间二阶的有限差分计算结果进行对比，结果表明同样精度下，径 向基函数每个波长所取的数据点数远小于空间四阶矩形网格有限差分每个波长所取的网格点数， 即径向基函数的空间采样率更低，这表明径向基函数具有更小的数值频散。 关 键 词径向基函数；声波方程；波场模拟；无网格法 中图分类号P631 文献标识码A DOI 10.3969/j.issn.1001-1986.2016.05.026 Numerical solution of the acoustic wave equation using meshless with radial basis functions JIANG Chanjun1，2， WANG Youxue1，2， LI Chuan1， ZENG Gaofu3， YUN Peng1 （1. College of Earth Science ， Guilin University of Technology， Guilin 541006， China； 2. State Key Laboratory of Continental Tectonics and Dynamics， Beijing 100037， China； 3. China Nonferrous Metals（Guilin） Geology and Mining Co.， Ltd.，Guilin 541004， China） Abstract Seismic exploration has been widely used in oil and gas as well as coal exploration. Numeric modeling of the seismic wave field is the footstone of the overal seismic exploration technique. The paper， introducing the radial basis function（RDF） into the numeric modeling of seismic wave field， constructiong in space the second de- rivative by using the meshless of the radial basis function， and in time domain using the simple second or- der difference ula， discussed emphatically the influence of the shape parameters（c） on the accuracy of the ， summed up the range of empirical value of c of 2～4 times of the avarage interval of data points. Different models were designed， sonic wave field was simulated by using meshless of radial basis function， and compared with the calculation results of the finite difference in four order space and two order time. The results show that at the same accuracy the number of data points taken by each wave length of the radial basis function was much less than the mesh points taken by each wave length of finite difference of rectangular mesh of four order space， that is， the spatial sampling rate of the radial basis function was lower， suggesting that the radial basis func- tion had smaller numerical dispersion. Key words radial basis function； sonic wave equation； simulation of wave field； meshless 地震勘探作为最重要的物探方法之一，广泛应 用于在油气、煤田勘探中。地震波场数值模拟是整 个地震勘探数据处理技术的基石。目前，基于波动 方程理论的地震波场数值模拟方法主要有有限差分 ChaoXing 第 5 期 蒋婵君等 求解声波方程的径向基函数无网格法 137 法（Finite Difference ， 简称 FDM）［1］、伪谱法［2］、 有限元法［3］、谱元法［4-5］、边界元法［6］等。通常意义 下，上述模数值拟方法各有优缺点［7］有限差分法 计算速度快，能够获得全波场信息，但是其频散较 严重，而且难以解决地表起伏问题；伪谱法计算精 度高，占用内存相对较小，并不利于并行计算，也 不适用复杂的地质体模型；有限元法处理地表起伏 问题的能力强，但计算量较大；谱元法和边界元法 便于处理边界起伏问题，但也有计算量大的不足。 1990 年，Kansa 首次将数学界研究的径向基函数 （Radial Basis Function，简称 RBF）用于求解偏微分 方程组，为求解偏微分方程提供了一种新兴的数值 方法［8-9］。径向基函数无网格法作为一种纯粹的无网 格法， 克服了传统数值方法对网格的依赖性， 适合于 复杂求解域问题或高维问题的计算， 因此该方法很快 得到了推广和发展［10］求解椭圆型偏微分方程［11］、 流体力学中的稳态解［12］、波动方程的频率域解［13］、 抛物型和双曲型偏微分方程以及瞬态问题和时域问 题的求解［14-15］等。本文将尝试用 Kansa 提出的径向 基函数无网格法［8-9］进行地震声波波场数值模拟。 1 基本原理 地震声波方程可以表示为 2 22 2 0 u vu t ∂ -∇ ∂ （1） 式中 v为波速； 2 ∇为拉普拉斯算子；u 为压力（或 标量位）。 在时间上采用简单的二阶差分方法进行计算， 式（1）写为 12221 2 kkkk UUvtUU - Δ ∇- （2） 本文在空间上采用径向基函数无网格法来构造 二阶导数 2u ∇，其实质是应用径向基函数构造插值 函数来逼近待求函数u ，再进行二阶求导。 以一维声波方程为例，给定数据 1 ｛，｝ j N jj xu n RR∈，选取径向基函数 RRφ →，利用平移构 造基函数系 1 ｛ （）｝N jj xxφ -，并寻找插值函数（ ）ux * 来逼近未知函数（ ）u x 1 （ ）（ ）（） N jj j u xuxxxλ φ * ≈- ∑ （3） 且满足插值条件 11 （）（） NN ijijjji jj uxxxλ φλ φ - ∑∑ （4） 式中 j λ为待求系数。 式（4）可写成 UΦΛ （5） 式中 12 ［，，，］T N u uuU， 12 ［，，，］T N λ λλΛ 11211 12222 12 （）（）（） （）（）（） （）（）（） N N NNNN xxx xxx xxx φφφ φφφ φφφ ■■ ■■ ■■ ■■ ■■ ■■ ⋮⋮⋮⋮ Φ 将式（5）代入式（2）中，可得 1222 （2） kkkk vtΛ -11 Δ ∇-Λ ΛΦ ΛΦΦΦ ΛΦΦ ΛΦΦΦ ΛΦ （6） 将初始条件代入式（6）， 即可沿时间向前推进求解。 利用径向基函数求解地震声波方程的过程中， 只使用求解空间的离散节点来进行计算，既不要求 节点间的连接信息，也没有形成单元，属于无网格 方法。 常见的径向基函数类型主要有 Kriging 的 Gauss 分布函数 2 2 （ ）e c r rφ - 、 Hardy 的 Multi-Quadric 函数 22 （ ）（）rcr β φ和 逆 Multi-Quadric 函 数（ ） rφ 22 （）cr β- ，其中 r 为数据点之间的距离范数，c 为 形状参数，β为正实数。 2 误差分析 基于径向基函数求解偏微分方程是一种新的 数值方法，尚未形成完整的误差分析理论。径向 基函数的类型、阶次和形状参数 c的取值，数据点 个数以及分布，都会影响径向基函数无网格法的 精度。 本文采用 Kansa 推荐的 Multi-Quadric 函数 22 （ ） rcrφ（β0.5），可以克服径向基函数插值 的不稳定性，从而提高计算精度［8-9］。数据点个数越 多，精度越高，但时间计算成本越高，数据点个数 选取应同时考虑精度和时间成本。数据点分布应尽 量均匀充满 R，本文将数据点规则分布，但需要指 出的是，数据点虽然呈网格状规则分布，但数据点 之间并无关联，这正是无网格法的基本思想。 如何确定形状参数 c 还没有很明确的办法， 一般 只能采取试算。 为了研究形状参数 c 的取值对精度的 影响，设计一个小模型，模型大小为 300 m300 m， 震源位于模型中心，波速 v2 000 m/s，震源选用 Ricker 子波，主频 0 f分别取 20 Hz、40 Hz 和 80 Hz， 时间步长为 Δt0.000 5 s，数据点间距 Δh5 m，数 据点 6161 个。根据数值模拟结果（图 1）可知，c 的 取值对精度影响较大，随着 c 的取值增大，精度变 高，但是当 c 超过 20 时，精度又迅速降低，c 经验 取值范围为 2～4 倍数据点间距。数据点个数不变， 当数据点随机分布时， 根据数值模拟结果（图 2）可知 ChaoXing 138 煤田地质与勘探 第 44 卷 c 经验取值范围为 2～4 倍平均数据点间距（图 2）。对 比图 1a 和图 2c，同样数据点个数，数据点均匀规 则分布比随机分布数值模拟效果更好。 3 波场模拟 对于一个二维均匀介质模型， 其大小为 2 400 m 2 400 m，波速 v2 000 m/s。震源位于模型中心 图 1 不同频率不同形状参数 c 取值下 T0.055 s 波场快照图 Fig.1 Snapshots of acoustic wave fields at time 0.055 sec under different values of shape parameters c and different frequency 图 2 不同数据点间距时不同形状参数 c 取值下 T0.055 s 波场快照图 Fig.2 Snapshots of acoustic wave fields at time 0.055 sec under different data point intervals and different values of shape parameters c ChaoXing 第 5 期 蒋婵君等 求解声波方程的径向基函数无网格法 139 （1200，1 200），选用 Ricker 子波，主频 0 f取 20 Hz。 时间步长 Δt0.002 s，数据点间距为 Δh30 m，形 状参数 c70，采用径向基函数无网格法（RBF）在 T0.5 s 时波场快照见图 3a；时间步长 Δt0.001 s， 网格间距 Δh5 m，空间四阶、时间二阶矩形网格有 限差分方法（FDM）在同一时刻的波场快照见图 3b。 这两种方法在（1 800，1 200）处的地震记录见图 3c。 比较可知，径向基函数方法和有限差分法均得到了 清晰的波场快照和地震记录。径向基函数方法在每 个波长取 3 个数据点左右就能达到空间四阶矩形网 格有限差分在每个波长取 20 个网格点的精度， 这表 明径向基函数具有更小的数值频散。 对于一个三层介质模型， 其大小为2 400 m2 400 m， 第一层厚度为 780 m，波速 2 500 m/s；第二层厚度为 720 m，波速为 2 000 m/s；第三层厚度为 900 m，波速 为 3 000 m/s。震源位于（1 200，1 140），子波亦选用 Ricker 子波。时间步长为 Δt0.002 s，数据点间距为 Δh30 m， 形状参数 c70， 径向基函数无网格法（RBF） 在 T0.5 s 时波场快照见图 4a；时间步长 Δt0.001 s， 网格间距 Δh5 m，空间四阶、时间二阶矩形网格有 限差分法（FDM）在相同时刻的波场快照见图 4b。径向 基函数方法和有限差分法在同一时刻波场快照图一 致，均能够清晰呈现反射波和透射波。这两种方法在 （1 400，1 140）处的地震记录重合， 均能能够准确记录直 达波和反射波波形，见图 4c。径向基函数方法在每个 波长取 3 个数据点左右就能达到空间四阶矩形网格有 限差分在每个波长取 20 个网格点的精度， 这表明径向 基函数具有更小的数值频散。 对于包含高速异常体的模型 2 400 m 2 400 m， 声波在背景区域中速度为 2 000 m/s，在异常体中的 速度为 6 000 m/s，异常体是边长为 20 m 的正方形， 异常体中心位于（1 450，1 450）。Ricker 子波主频为 40 Hz，震源位于模型中心（1 200，1 200）。时间步长 Δt0.001 s，数据点间距为 Δh20 m，形状参数 c80，采用径向基函数无网格法（RBF）在 T0.5 s 时波场快照见 5a；时间步长 Δt0.000 5 s，网格间 距 Δh5 m，空间四阶、时间二阶有限差分方法 （FDM）在同一时刻的波场快照见图 5b。 这两种方法 在（1 440，1 200）处的地震记录见图 5c。比较可知， 径向基函数无网格法可以清晰地呈现出直达波和异 常体产生的绕射波。在相同精度下，径向基函数方 法数据点总数是空间四阶矩形网格有限差分网格点 总数的 1/4，即径向基函数的空间采样率更低，这表 明径向基函数具有更小的数值频散。 图 3 均匀介质模型 RBF 和 FDM 模拟结果对比 Fig.3 Comparison of results generated by RBF and FDM simulation of models for isotropic medium ChaoXing 第 5 期 蒋婵君等 求解声波方程的径向基函数无网格法 141 4 结 语 本文尝试将径向基函数无网格法引入声波方程 的求解中，在空间上用径向基函数无网格法来构造 二阶导数，而在时间上采用简单的二阶差分公式， 并重点讨论了形状参数 c 对该方法精度的影响，认 为 c 经验取值范围应为为 2～4 倍平均数据点间距。 通过对均匀介质、 层状均匀介质模型及异常体模型， 利用径向基函数无网格法进行声波波场模拟，并与 空间四阶、时间二阶的有限差分计算结果进行了对 比，结果表明同样精度下，径向基函数每个波长 所取的数据点数远小于空间四阶矩形网格有限差分 每个波长所取的网格点数，即径向基函数的空间采 样率更低， 这表明径向基函数具有更小的数值频散。 这是因为径向基函数法属于无网格法，数据点之间 没有关联，不像传统的数值方法强烈依赖网格的选 取，不会产生有限差分法或者有限元因网格节点选 取导致求解失败等问题。径向基函数法是地震波场 数值模拟的一种新尝试，值得进一步探讨研究。 参考文献 ［1］ VIRIEUX J. P-SV wave propagation in heterogeneous media velocity-stress finite-difference ［J］. Geophysics，1986， 51（4）889-901. ［2］ KOLSOFF D D，BAYSAL E. Forward modeling by a Fourier ［J］. Geophysics，1982，47（10）1402-1412. ［3］ 张美根，王妙月，李小凡，等．时间域全波场各向异性弹性 参数反演［J］．地球物理学报，2003，46（1）94-100． ZHANG Meigen，WANG Miaoyue，LI Xiaofan，et al. Full wavefield inversion of anisotropic elastic parameters in the time domain［J］. Chinese Journal of Geophysics， 2003， 46（1） 94-100. ［4］ SERIANI G，PRIOLO E. Spectral element for acoustic wave simulation in heterogeneous media［J］. Finite Elements in Analysis and Design，1994，16（3-4）337-348. ［5］ KOMATITSCH D，BARNES C，TROMP J. Wave propagation near a fluid-solid interfacea spectral-element approach［J］. Geophysics，2000，65（2）623-631. ［6］ BOUCHON M，COUTANT O. Calculation of synthetic seis- mograms in a laterally varying medium by the boundary ele- ment-discrete wave number ［J］. Bulletin of the Seismol- ogical Society of America，1994，84（6）1869-1881. ［7］ 朱多林，白超英．基于波动方程理论的地震波场数值模拟方 法综述［J］．地球物理学进展，2011，26（5）1588-1599． ZHU Duolin，BAI Chaoying. Review on the seismic wavefield forward modelling［J］. Progress in Geophysics，2011，26（5） 1588-1599. ［8］ KANSA E J. MultiquadricsA scattered data approximation scheme with applications to computational fluid dynamics-I sur- face approximations and partial derivative estimates［J］. Com- puters and Mathematics with Applications，1990，19（8/9） 127-145. ［9］ KANSA E J. MultiquadricsA scattered data approximation scheme with applications to computational fluid dynamics-II so- lutions to parabolic， hyperbolic and elliptic partial differential equations［J］. Computers and Mathematics with Applications， 1990，19（8/9）147-161. ［10］ 李坤，黄其柏，林立广，等．二阶时域波动方程的无网格方 法求解［J］．华中科技大学学报（自然科学版），2011，39（3） 26-29． LI Kun，HUANG Qibai，LIN Liguang，et al. Solving of sec- ond-order time domain wave equations by meshless ［J］. Journal of Huazhong University of Science of Technol- ogy（Natural Science Edition），2011，39（3）26-29. ［11］ FAIRWEATHER G，KARAGEORGHIS A. The of fundamental solution for elliptic boundary problems［J］. Ad- vances in Computational Mathematics，1998，969-95. ［12］ CHEN W， HON Y C. Numerical investigation on convergence of boundary knot in the analysis of Helmholtz，modified Helmholtz and convection-diffusion problems［J］. Computer s in Applied Mechanics and Engineering， 2003， 192（15） 1859-1875. ［13］ 李鹏，彭伟才，李志江，等．加权最小二乘无网格法求解亥 姆霍兹方程［J］． 华中科技大学学报（自然科学版）， 2010， 38（7） 40-43． LI Peng，PENG Weicai，LI Zhijiang，et al. Meshless least- squares for solving Helmholtz equation［J］. Journal of Huazhong University of Science of Technology（Natural Science Edition），2010，38（7）40-43. ［14］ DEHGHAN M，SHOKRI A. A meshless for numerial solution of the one-dimensional wave equation with an integral condition using radial basis functions［J］. Numerical Algorithms， 2009，52461-477. ［15］ YOUNG D L，GU M H，FAN C M. The time-marching of fundamental solutions for wave equations［J］. Engineering Analysis with Boundary Elements，2009，33（12），1411-1425. ［16］ 赖生建．计算电磁学中的径向基无网格法［D］．西安西安电 子科技大学，2010. （责任编辑 聂爱兰） ChaoXing