基于弹性波速度-应力方程与交错网格有限差分的地震波场模拟器

简介本资源是一款面向地球物理专业研究人员、勘探工程师及高年级本科生的地震波场正演模拟软件基于弹性波速度-应力波动方程采用交错网格有限差分法实现P波与S波在复杂介质中的高精度传播模拟有效支撑地震成像、反演建模与教学实验等核心任务。压缩包共76个文件含7个核心cpp源码、4个h头文件、5个可执行exe程序、4个dat模型与记录数据、1个完整Visual Studio解决方案.sln及配套UI界面.ui、资源文件.qrc和详细使用手册.doc总大小477.28MB结构清晰支持模型导入、参数配置、波场快照动态可视化及单炮地震记录输出。已有736人学习下载用户可直接运行exe进行波场模拟亦可基于完整C工程源码深入理解交错网格离散策略、内存管理机制ManageMemory.h/cpp与速度-应力耦合更新逻辑具备良好的可扩展性与二次开发基础。 前阵子做一批模型数据被手头的射线法程序坑了一把断层夹在两层高速体中间射线追踪出来的反射同相轴怎么看怎么不对劲。跟同事一合计干脆直接上波场正演用弹性波速度-应力波动方程和交错网格有限差分自己写一个地震波场模拟软件。本来是应急用的结果越写越顺手最后整理成了一份带完整源代码的小工程从均匀介质验证到层状模型反射记录前后折腾了一个多星期。这篇就把这个模拟器从方程推导到网格设计、边界处理、源码组织、算例验证完整拆开讲一遍里面哪些坑我已经替你们踩平了也会一并交代。1. 为什么弃用位移方程而选速度-应力方程项目定方向的时候第一个要定的是用哪种方程形式。这个选择直接决定了后面所有代码怎么写所以值得单独拿出来聊。1.1 同一个波动两种数学描述弹性波传播的经典描述是位移形式的二阶偏微分方程。对各向同性介质Lamé常数λ和μ、密度ρ给定时位移u满足ρ ∂²u/∂t² (λμ)∇(∇·u) μ∇²u f这是课本上最常见的写法很多初学者上手就写这个。数值求解时一般把二阶时间导数化为两个一阶导数或者直接用中心差分同时处理时间和空间。但这么做有几个麻烦二阶空间导数需要较大的差分模板边界条件的物理意义不直观而且要把位移场离散后还要再求梯度才能得到应力而检波器记录和很多工程判断恰恰依赖的是应力或速度分量。另一种做法是引入质点速度 v ∂u/∂t 和应力张量 σ得到一阶速度-应力方程组。二维情况下弹性各向同性介质里有五个未知量水平速度 vx、垂直速度 vz、正应力 σxx、σzz、切应力 σxz。方程组写成ρ ∂vx/∂t ∂σxx/∂x ∂σxz/∂z ρ ∂vz/∂t ∂σxz/∂x ∂σzz/∂z ∂σxx/∂t (λ2μ)∂vx/∂x λ∂vz/∂z ∂σzz/∂t λ∂vx/∂x (λ2μ)∂vz/∂z ∂σxz/∂t μ(∂vx/∂z ∂vz/∂x)这组方程在数学上是典型的双曲型守恒律方程组非常适合用时间递推方式求解。速度分量在时间n1/2步更新应力分量在时间n1步更新两者交替推进就是常说的蛙跳格式。1.2 速度-应力形式对编程的三个实在好处第一个好处是时间方向只需一阶导数递推公式极简。速度的更新依赖应力的空间导数应力的更新依赖速度的空间导数时间步进不需要二阶导数因此每个时间步的计算量非常固定。第二个好处是能和交错网格天然配合。速度定义在网格半节点、应力定义在整节点空间导数在对方所在位置求值每一步差分都是两点之差模板紧凑数值频散特性明显优于非交错网格的同阶格式。这个我后面单独展开。第三个好处是扩展性强。粘弹性、各向异性、衰减介质、孔隙介质这类物理机制每加一个往往只需要在速度-应力方程右端增加一个辅助变量或附加项。如果当初用的是位移二阶方程改造成本会大很多。从我实际写代码的经验看速度-应力方程还有一个隐蔽优势稳定性和频散表现更直观。调试时如果波形出现高频振荡第一反应查时间步长是否超CFL如果出现拖尾频散第一反应查空间采样密度不需要在一堆二阶导数的截断误差里找原因。2. 交错网格不是错开一点那么简单交错网格staggered grid是这套模拟器性能的核心。很多初次接触的人以为只是把网格错开半个格距实际上这里面的物理量排布和差分系数推导有好几个容易踩坑的地方。2.1 五个场量、四套格点的排布二维情况下五个场量不是全都放在同一个网格点上的。经典做法是σxx、σzz、ρ 定义在整网格点 (i, j)vx 定义在水平半网格点 (i1/2, j)vz 定义在垂直半网格点 (i, j1/2)σxz 定义在对角半网格点 (i1/2, j1/2)为什么要这样放以运动方程第一式为例∂σxx/∂x 要在 vx 所在位置求值而 σxx 正好定义在整网格点(σxx_{i1,j} - σxx_{i,j})/dx 给出的就是半网格点 (i1/2, j) 处的一阶导数这个位置恰好也是 vx 所在的格点。再看 ∂σxz/∂zσxz 定义在对角半网格点沿 z 方向做差分同样落在 (i1/2, j)。两个导数在同一个位置求值相加后直接更新 vx完美的对齐。我第一版代码把 σxz 放错了位置放进整网格点结果更新 vx 和 vz 时总要费劲做四邻域平均频散也明显偏大。后来老老实实按对角半网格点存放代码简洁了波形也干净了。2.2 高阶差分模板的系数交错网格的空间差分可以用二阶、四阶、六阶甚至更高精度。二阶最简单∂σxx/∂x ≈ (σxx_{i1,j} - σxx_{i,j}) / dx这就是两点差分精度只有二阶。实际模拟中二阶格式频散严重每波长通常要20个以上网格点计算量偏大。所以多数生产级模拟器用四阶或更高阶。交错网格四阶差分模板的系数和普通中心差分不一样这点很多人会弄混。普通中心四阶差分系数是 2/3 和 -1/12但交错网格四阶格式的系数是∂f/∂x ≈ (9/8)(f_{i1} - f_i) - (1/24)(f_{i2} - f_{i-1}) 再除以 dx注意这里 f 是整网格点的量导数结果落在半网格点上。系数 9/8 和 -1/24 是按 Taylor 展开消去三阶截断误差推出来的。六阶可以用类似办法递推我在这套代码里默认用四阶稳定性和计算量平衡得最好。写代码时这些系数是硬编码还是动态算初期建议硬编码跑起来再看效果。等确实需要高阶格式了再写一个系数生成函数。2.3 时间递推的顺序问题时间上采用蛙跳格式已知应力场时间层 n和速度场时间层 n-1/2用应力场空间导数更新速度场到 n1/2 层用速度场空间导数更新应力场到 n1 层循环速度更新和应力更新在时间上错开半个步长这个顺序不能颠倒。我最早把更新顺序写反了算出来的波场整体错了一个时间步做解析解对比时相位对不上查了半天才反应过来。震源加载也要注意时间位置。力源加在速度方程右端严格说应该加在速度更新的半步层上爆炸源加在应力方程右端应该加在应力更新的整步层上。代码里可以一边更新一边把子波值叠加上去。3. 边界处理从衰减带换到PML的折腾有限计算区域必然产生边界反射。地震波场模拟对边界反射容忍度极低因为一次反射波和边界假象在观测记录上可能完全混在一起。边界处理是整个模拟器里最影响成像质量的部分。3.1 衰减带本文还有配套的精品资源点击获取