波动方程有限差分数值模拟研究_龚良钢.pdf

2020年第12期西部探矿工程 * 收稿日期 2020-03-22 作者简介 龚良钢 （1995-） ， 男 （汉族） ， 江西南昌人， 东华理工大学在读硕士研究生， 研究方向 测井资料解释评价。 波动方程有限差分数值模拟研究 龚良钢* 东华理工大学地球物理与测控技术学院， 江西 南昌330013 摘要 地震波波场数值模拟是地震学领域中的一种重要的方法， 在油气田勘探中起到重要作用。 有限差分法具有精度高、 计算速度快等优点， 在数值模拟研究中受到广泛的应用。以波动理论和泰 勒级数展开公式为基础， 通过对声波波动方程进行推导， 最终得出二阶波动方程的有限差分格式。 并对稳定性及收敛性、 震源函数、 频散问题、 吸收边界条件等几个关键问题进行分析。建立在理论的 基础上， 利用Matlab程序设计得到二维地震波有限差分数值模拟， 主要对两层水平层状介质进行模 拟， 得到相对应的地震记录。研究表明， 波动方程有限差分数值模拟可以很方便的计算出不同地质 体的空间位置， 形态的地震记录， 具有高效、 精确等特点， 为进一步了解地下介质中地震波的传播规 律提供了理论依据。 关键词 波动方程； 有限差分； 数值模拟 中图分类号 P631.44 文献标识码 A 文章编号 1004-5716202012-0085-05 1概述 地震波场数值模拟的实质 假设在地下的介质构造 模型和相应的物理参数已经知道的情况下， 模拟地震波 在地下的传播规律， 计算在地下和地面观测点观测到数 据记录， 最终画出地震波的快照、 记录等[1-7]。其理论基 础就是地震波在地下各种介质中传播， 根据波在不同介 质中的传播速度来判别地下地质的形态特征等。因此 波场数值模拟在地震勘探中起到非常重要的作用。 地震波场的数值模拟可以分为三大类 波动方程 法、 射线追踪法、 积分方程法。而波动方程法中几种常 用的逼近方法又可以分为 伪普法、 谱元法和有限差分 法[8]。利用有限差分法的算法较为稳定， 适应性相对来 说较好， 能够得到在各种介质及复杂的地层的传播规律 与理论结果较为匹配的模拟， 是地震波场数值模拟中常 用的方法之一。本文基于Matlab程序， 采用波动方程有 限差分法对水平层状介质进行数值模拟， 模拟声波在水 平层状介质中的传播过程， 并对不同时间的波场快照进 行比较和对最终得到的地震记录进行分析， 可以更好地 理解地震波在分层介质中的传播规律。综上所述， 为了 深入地了解地震波在地下不同介质中的传播规律和准 确地解释实际勘探中的地震资料， 地震波场的数值模拟 何尝不是一种行之有效的方法之一。 2波动方程数值模拟原理 二维均匀介质声波波动方程为 1 v2 ∂2u ∂t2 ∂2u ∂x2 ∂2u ∂z2 fx,z,t（1） 式中 v声波波速； fx,z,t震源函数； ux,z,t声波波场值。 令ukm,numΔx,nΔz,kΔt， 其中Δx、Δz是X、 Z 轴方向上的间隔，Δt是时间间隔。用k表示时间方向 的离散网格， m表示x方向的离散网格， n表示z方向的 离散网格。利用泰勒级数将uk1 m,n在u k m,n展开， 有 uk1 m,nux,z,tΔt ukm,nΔt∂u ∂t km,n 1 2Δt 2∂2u ∂t2 km,n 1 6Δt 3∂3u ∂t4 km,nOΔt4 （2） 同理， 利用泰勒级数将ukm,n在uk-1 m,n展开有 uk-1 m,nu k m,n-Δt∂u ∂t km,n 1 2Δt 2∂2u ∂t2 km,n- 1 6Δt 3∂3u ∂t4 km,nOΔt4 （3） 地质与矿业工程 85 2020年第12期西部探矿工程 将 （2） 式与 （3） 式两式相减得 ∂u ∂t k-1 m,n uk1 m,n-u k-1 m,n 2Δt 上式即为关于t的一介中心差分。 从而进一步可以得到关于t的二阶差分 ∂ 2u ∂t2 km,n uk1 m,n-2u k m,nu k-1 m,n Δt2 （4） 同理， 可以得到关于x、 z的二阶中心差分 关于x方向 ∂ 2u ∂x2 k m,n uxΔx-2uxux-Δx Δx2 ukm1,n-2ukm,nukm-1,n Δx2 （5） 关于z方向 ∂ 2u ∂z2 km,n uzΔz-2uzuz-Δz Δz2 ukm,n1-2ukm,nukm,n-1 Δz2 （6） 将 （4） 式、（5） 式、（6） 式代入 （1） 式得 1 v2 ∂ 2u ∂t2 km,n∂ 2u ∂x2 k m,n∂ 2u ∂z2 km,n ukm1,n-2ukm,nukm-1,n Δx2 ukm,n1-2ukm,nukm,n-1 Δz2 1 v2 uk1 m,n-2u k m,nu k-1 m,n Δt2 令ΔxΔzh（即网格相等） ， 则 1 v2 uk1 m,n-2u k m,nu k-1 m,n Δt2 ukm1,n-2ukm,nukm,n1-2ukm,nukm,n-1 Δh2 将上式继续推导得 uk1 m,n2u k m,n-u k-1 m,n v2Δt2 h2 [ukm1,nukm-1,nukm,n1ukm,n-1-4ukm,n] 21-2M2ukm,n-uk-1 m,nM 2uk m1,nu k m-1,nu k m,n1u k m,n-1 其中M vΔt h ， 上式就是二阶波动方程的有限差 分格式。 3波动方程有限差分的几个关键问题 3.1边界条件问题 由于在计算机上进行数值模拟时受到计算机的内 存和运算速度的因素影响， 因此， 波动方程数值模拟只 能模拟一个有限的区域， 那么在这个有限的区域四周 就会有边界， 当波传播到人工边界时， 它将发生反射， 这将导致最终模拟效果受到干扰[9]。显然在如图1中， 在波传播到该区域的边上时， 发生了反射情况， 从而影 响最终的模拟效果。所以这种情况是我们需要避免 的， 因此， 在我们需要对人工边界进行处理， 消除或减 弱这种边界反射条件。 图1未加入边界条件的波场快照 而为了消除或减弱人工边界条件， 反映波在真实 地下地质的传播情况， 从而达到准确地反映地下地质 情况。我们通常是对人工边界进行处理 即加入人工 吸收边界条件。当我们加入了人工吸收边界条件后， 波运动到区域边界后被吸收， 即得到了无限区域的正 演结果。我们在处理边界反射常用的处理方法有 透 明边界条件和吸收边界条件等。从图2中可以看出， 在 我们加入了吸收边界条件后， 边界所发生的反射效果 大大减弱。 图2加入边界条件后的波场快照 86 2020年第12期西部探矿工程 3.2稳定性及收敛性 在波动方程有限差分数值模拟中， 首先要考虑的 问题就是其稳定条件[10]。因为利用有限差分法对地震 波场数值模拟时， 是随时间递推计算的， 所以计算下一 个时刻的波场值会受到上一时刻的波场值影响。这就 使得误差会随时间的推移而累加， 导致最终的模拟效 果受到影响。如果差分方程的解的误差不随时间推进 而累计， 则差分方程的解是稳定的。稳定性反映了差 分格式在计算过程中控制误差传递的能力[11]。有限差 分稳定性的条件也根据差分阶数而变化， 差分阶数越 高， 稳定性越差。根据Lax等价定理， 在保证稳定性， 同时也保证了差分格式的收敛性[12]。 对于均匀介质， 差分方程的解的稳定性是确定的， 即稳定条件为 vΔt Δx ≤ 1 2,即Δt≤ Δx v 2 其中v是波在均匀介质中的传播速度，Δx表示 （Δx，Δz） 的最小值。 对于非均匀介质， 就与波在不同介质中的传播速 度有关， 其稳定条件为 vmaxΔt Δx ≤ 1 2,即△t≤ Δx vmax2 其中vmax是各种均匀介质中的最大的波速值。这 里Δx是表示 （Δx，Δz） 的最小值。 3.3震源函数 在有限差分波动方程正演模拟中， 震源项通常表 示为震源空间函数与震源时间函数乘积的形式。时间 函数描述震源力作用时间的延续性， 通常用地震子波 函数表示;空间函数描述震源震源力的空间作用范围， 通常用指数衰减函数表示[13]。对地震子波而言， 子波 延续的时间越短， 频带越宽， 子波的垂直分辨率也就越 高。在做有限差分计算时避免不了会发生数值频散现 象， 模拟结果受到震源函数选择的影响， 因此我们需要 合理地选择子波主频。所以想要更加深入地了解和完 善地震波的传播知识， 在进行波动方程数值模拟前， 必 须要先选定合适的地震子波。目前普遍认为雷克提出 的地震子波数学模型具有广泛的代表性， 其数学表达 式为 bt1-2πft2exp-πft2 其中f为主频， t为时间。 图3就是选用主频为30Hz的雷克子波 （a） 及其子 波振幅谱 （b） 。根据其波形不难看出， 利用雷克子波作 为一个震源函数效果较好。 3.4频散问题 波动方程有限差分法的最大缺点就是频散问题， 在有限差分数值求解的过程中， 是利用离散化的有限 差分方程去逼近波动方程， 当波场内空间的采样点太 少时， 进而产生较大的误差， 即产生数值频散现象， 从 而降低数值模拟结果分辨率， 严重的话会使得结果产 生错误。这种频散现象它是有限差分法求解波动方程 数值解的本质特性， 是其内部因素引起的， 频散不可避 免， 只能尽量降低[14]。所以， 一定要在数值模拟的时 候， 尽力地降低频散的程度。一般情况下， 在考虑到计 算机的运行内存和运行速度的前提下， 适当的较少时 间和空间网格步长， 和提高差分的阶数， 又或者降低地 震子波的频率。这些方法都可以有效地降低频散现 象。 4模型分析 设计双层介质模型， 网格大小200200， 第一层为 0～70m， 速度为2400m/s,第二层为70～200m， 速度为 3000m/s， 采样时间间隔为0.001， 震源位置 （100,50） ， 子波选用雷克子波， 子波频率为30Hz。见图4。 正演结果 图5 （a） 反映的是波动刚产生不久时的波场状态， 图5 （b） 反映的是波场在第二层产生的反射。图5 （c） 反 映的是第二层反射波传到地面被检波器记录的波场状 态。很明显当地震波入射到两种不同的介质时， 会产 生透射波和反射波。在地震记录中， 我们可以看到激 发震源后产生的直达波， 以及在70m分层处发生的反 射， 即说明了这是一个双层水平层状介质的结构， 与设 图330Hz雷克子波和子波振幅谱 87 2020年第12期西部探矿工程 计的模型基本一致， 达到所预期的效果。同时从图5 （a） 与5 （b） 中可以看出， 数值模拟计算所得的结果与理 论基本一致， 反射和透射的状态清晰明了， 利用有限差 分法能够较好地完成模拟地震波场的计算。 5结论 （1） 利用计算机对地震波场数值模拟需要考虑到 边界反射的问题， 即对人工边界进行处理， 加入人工吸 收边界条件来改善数值模拟效果。 （2） 选用合适的震源子波是一个很重要的问题， 因 为子波频率对成像结果有很大的影响， 虽然提高子波 频率有助于提高分辨率， 但同时也会发生能量衰减及 频散问题， 影响成像效果。 （3） 既要提高计算效率又要降低频散现象， 需要选 取合理的时间和空间网格步长。进一步提高对波动方 程有限差分数值模拟出现频散问题的认识。 图4双层介质模型 图5正演结果 （下转第94页） 88 2020年第12期西部探矿工程 （2） 结合已知钻孔及矿体， 初步建立了物探异常模 型。即 AMT显示为低电阻率异常 （ρ＜3200Ωm） ； SIP 显示低电阻率 （ρ＜1000Ωm） 、 高充电率 （m＞3） 、 高时 间常数 （τ＞3.5S） 、 低频率相关系数 （c＜0.45） 的异常组 合特征。 （3） 综合物探方法研究表明 音频大地电磁方法 （AMT） 测深在该区取得了良好的应用效果， 与频谱激 电法 （SIP） 测深资料， 相互验证、 相互补充， 并经异常验 证见矿， 验证了电磁方法在寻找深部隐伏铅锌矿体中 的有效性， 该方法可以为同模式铅锌矿攻深找盲提供 技术支撑。 参考文献 [1]金中国.黔西北地区铅锌矿控矿因素、 成矿规律与找矿预测 研究[D].中南大学,2006. [2]唐骄.贵州赫章县三万坪铅锌矿矿石特征[J].中国西部科技, 2014,131018,36. [3]秦小军.赫章某铅锌矿床地质特征及矿床成因浅析[J].企业 技术开发,2014,331963-64,67. [4]贵州省地矿局.贵州省区域地质志[M].北京地质出版社, 1987501-505. [5]卢卯.瞬变电磁法在黔西北找矿成果的研究[C]//第九届中国 国际地球电磁学术讨论会论文集.中国地球物理学会 （地球 电磁专业委员会） [Geo-Electromagnetic Special Committee of CGS]中国地球物理学会,2009152-158. [6]马为,李世斌,徐新学,郑军,赵洪鹏,李华强.重磁电联合反演 方法在天津基岩构造研究中的应用[J].地质调查与研究, 2014,373212-216. Application of Audio Magnetotelluric and Spectrum Induced Polarization in the Detection of a Lead-Zinc Mine in Northwest Guizhou HAN Yao-fei， YANG Bing-nan， ZHANG De-quan， LI Shu-lin The 103 Geological Brigade of Guizhou Bureau of Geology and Mineral Resources Exploration and Development，Tongren Gui- zhou 554300, China Abstract The lead-zinc deposit in northwestern Guizhou has the potential to find medium and large deposits. According to the geo- logical background of the mining area and the electrical character- istics of the rock and ore, audio frequency magnetotelluric and spectral polarization soundings were carried out, combined with relevant geological data, to determine the electrical character- istics of the rock and ore body in the deep space and the abnormal response of the ore body Characteristics, delineate the structural pattern, find deep hidden ore bodies, and break through the bottle- neck of prospecting. The research results show that the two geo- physical s have achieved good working results in the area, and the mines have been verified by drilling to verify the effective- ness of the electromagnetic in finding deep hidden lead- zinc ore bodies. The can be the same model of lead Zinc mines provide technical support by finding the blind. Key words lead-zinc ore；audio magnetotelluric ；spec- trum induced polarization； Northwest Guizhou （上接第88页） 参考文献 [1]裴正林,牟永光.地震波传播数值模拟[J].地球物理学进展, 20044213-221. [2]牟永光,裴正林.三维复杂介质地震数值模拟[M].北京石油 工业出版社,2005. [3]李信富,李小凡,张美根.地震波数值模拟方法研究综述[J].防 灾减灾工程学报,2007,272241-248. [4]杨莹.二维地震波场有限差分法数值模拟研究[D].中国地质 大学 （北京） ,2009. [5]刘庆敏.高阶差分数值模拟方法研究[D].中国石油大学华 东,2007. [6]王修敏.神通正演模拟子系统研发与应用以二维声波 正演模块的应用为例[J].油气地球物理,2016,16224-26,31. [7]杨晓柳,竺俊,姚铭,等.二维声波方程交错网格有限差分数值 模拟研究[J].中国锰业,2018. [8]冯英杰,杨长春,吴萍.地震波有限差分模拟综述[J].地球物理 学进展,2007,222487-491. [9]关东.二维声波波动方程波场正演的计算机数值模拟的研 究[D].成都理工大学,2013. [10]袁崇鑫.波动方程有限差分正演技术研究[D].2015. [11]李国平,姚逢昌,石玉梅， 等.有限差分法地震波数值模拟的 几个关键问题[J].地球物理学进展,2011,262469-476. [12]倪长宽.复杂介质地震波场正演模拟方法研究[D].中国石 油大学,2008. [13]黎殿来.波动方程时域有限差分地震正演建模方法研究 [D].电子科技大学,2012. [14]杨顶辉,滕吉文.各向异性介质中三分量地震记录的FCT有 限差分模拟[J].石油地球物理勘探,19972181-190. 94