简介本资源是一套面向地球物理勘探与计算地球科学方向研究者及高年级本科生的波场数值模拟工具集聚焦各向同性、VTI、TTI及双相介质中的弹性/声波传播建模解决复杂介质中地震波正演模拟与参数敏感性分析等核心问题。压缩包含1049个文件主体为191个Fortran源码.f90、112个模型参数文件.dat、96个编译模块.mod及50个标准SGY格式地震数据文件辅以PDF论文、JPG/HTM说明文档和EXE可执行程序总容量101.77MB结构完整、即装即用。已有466人学习下载涵盖从理论推导到代码实现、从模型构建断层/层状/TPM双相到结果可视化JPG/GRF/PNG的全链路支撑。用户可直接调用旋转交错网格RSG与交错网格SSG两类主流差分算法程序复现经典文献案例并借助Bond变换、Thomsen参数计算等配套工具开展介质参数反演与波场特征分析。1. 用 Fortran 实现地震波传播模拟从各向同性到双相介质为什么数值建模必须分层选型你正在调试一个地震波场模拟程序输入参数改了三遍合成记录却始终和实测数据对不上——不是振幅衰减过快就是横波分裂特征消失。问题未必出在代码 bug 上而可能源于介质模型选型失当把实际为 VTI垂直横向各向异性的页岩层当成各向同性处理纵波走时误差可达 8% 以上若忽略双相介质中流体-骨架耦合效应低频段频散特征将完全丢失。本篇聚焦numerical-modeling-of-wave-field的 Fortran 实现路径覆盖各向同性、VTI、TTI 及双相介质四类核心模型。不讲抽象偏微分方程推导只拆解 Fortran 代码中如何用弹性张量、Christoffel 方程、Biot 理论控制波速各向异性与耗散机制。适合已掌握有限差分基础、正为野外数据反演卡点的地球物理工程师也适合需复现经典论文数值实验的研究生——所有代码块可直接编译运行参数表标注真实地质场景取值范围。2. 各向同性与 VTI 介质的弹性张量构建Fortran 中张量索引与内存布局的硬约束2.1 各向同性介质Lamé 参数到刚度矩阵的 Fortran 映射规则各向同性介质的刚度矩阵仅由 Lamé 参数 λ 和 μ 决定但 Fortran 数组下标从 1 开始且按列主序存储直接套用数学公式易导致内存越界。正确实现需将 6×6 刚度矩阵 $C_{ij}$ 按 Voigt 记号展开并严格对应 Fortran 二维数组C(6,6)的索引! 各向同性刚度矩阵 C_ij (Voigt 记号: 1xx, 2yy, 3zz, 4yz, 5xz, 6xy) real(kind8) :: lambda, mu, C(6,6) C 0.0d0 C(1,1) lambda 2.0d0*mu; C(1,2) lambda; C(1,3) lambda C(2,1) lambda; C(2,2) lambda 2.0d0*mu; C(2,3) lambda C(3,1) lambda; C(3,2) lambda; C(3,3) lambda 2.0d0*mu C(4,4) mu; C(5,5) mu; C(6,6) mu注意Fortran 中C(1,1)对应数学上的 $C_{11}$而非 $C_{00}$。若误用C(0,0)将触发段错误若未初始化为零残留内存值会导致波场出现伪频散。2.1.1 参数校验λ 和 μ 的物理边界约束Lamé 参数必须满足稳定性条件$\mu 0$ 且 $\lambda 2\mu 0$保证纵波速度 $V_p \sqrt{(\lambda2\mu)/\rho} 0$。典型页岩参数为 $\lambda12.5$ GPa、$\mu7.2$ GPa若输入 $\lambda-5$ GPa程序虽能编译但计算出的 $V_p$ 为虚数后续时间步迭代发散。建议在read_input子程序中加入断言if (mu 0.0d0 .or. lambda 2.0d0*mu 0.0d0) then write(*,*) Error: Lamé parameters violate stability condition stop end if2.2 VTI 介质Thomsen 参数与 5 参数刚度矩阵的 Fortran 实现VTI垂直横向各向异性介质需 5 个独立弹性常数工程上更常用 Thomsen 参数ε, δ, γ描述各向异性强度。Fortran 实现的关键是将 Thomsen 参数映射到 Voigt 刚度矩阵避免直接操作 21 个独立常数。核心转换关系如下ρ 为密度$V_{p0}$ 为垂直纵波速度刚度分量表达式$C_{11}$$\rho V_{p0}^2 (1 2\varepsilon)$$C_{33}$$\rho V_{p0}^2$$C_{44}$$\rho V_{s0}^2$$C_{13}$$\rho V_{p0}^2 \sqrt{2\delta \frac{V_{s0}^2}{V_{p0}^2}} - \rho V_{s0}^2$$C_{66}$$C_{44} \rho V_{p0}^2 \gamma$! VTI 刚度矩阵构建输入vp0, vs0, rho, eps, delta, gamma real(kind8) :: vp0, vs0, rho, eps, delta, gamma, C(6,6) C 0.0d0 C(1,1) rho * vp0**2 * (1.0d0 2.0d0*eps) C(2,2) C(1,1) ! C22 C11 for VTI C(3,3) rho * vp0**2 C(4,4) rho * vs0**2 C(5,5) C(4,4) ! C55 C44 C(6,6) C(4,4) rho * vp0**2 * gamma C(1,3) rho * vp0**2 * sqrt(2.0d0*delta (vs0/vp0)**2) - C(4,4) C(3,1) C(1,3) ! 对称性2.2.1 Thomsen 参数的地质意义与典型取值εepsilon控制水平方向纵波速度各向异性页岩 ε ∈ [0.1, 0.3]砂岩 ε 0.05δdelta影响非垂直入射时的纵波速度δ 0 导致 NMO 速度随炮检距增大而降低γgamma剪切波各向异性γ ≈ (Vs_h - Vs_v)/Vs_v页岩 γ ∈ [0.15, 0.25]若 δ 为负值如某些碳酸盐岩sqrt(2.0d0*delta ...)将产生 NaN必须前置检查if (2.0d0*delta (vs0/vp0)**2 0.0d0) then write(*,*) Warning: delta too negative for VTI model, setting to minimum delta -0.5d0*(vs0/vp0)**2 1.0d-6 end if2.3 TTI 介质旋转坐标系下的刚度矩阵变换TTI倾斜横向各向异性需在 VTI 基础上叠加倾角 θ 和方位角 φ。Fortran 实现不采用解析表达式过于冗长而是通过坐标系旋转矩阵 R 实现张量变换$C R \cdot C \cdot R^T$。关键在于旋转矩阵 R 的 Fortran 构造——必须使用三维欧拉角标准顺序Z-X-Z! TTI 旋转先绕 z 轴转 φ再绕新 x 轴转 θ最后绕新 z 轴转 ψ通常 ψ0 real(kind8) :: theta, phi, R(3,3), C_vti(6,6), C_tti(6,6) ! 构造 R 矩阵省略具体三角函数计算见附录子程序 rotate_tensor call construct_rotation_matrix(R, theta, phi, 0.0d0) ! Voigt 记号下的张量旋转需专用子程序非简单矩阵乘法 call rotate_stiffness_tensor(C_vti, R, C_tti)提示Voigt 记号下刚度矩阵旋转不能直接用matmul(R, matmul(C_vti, transpose(R)))因 Voigt 格式破坏了张量秩。必须调用专门的rotate_stiffness_tensor子程序该子程序将 6×6 矩阵还原为 3×3×3×3 四阶张量执行 $C{ijkl} a{im} a_{jn} a_{kp} a_{lq} C_{mnpq}$再压缩回 Voigt 形式。开源库如seiscl提供此功能但自研时需严格验证旋转前后 $C_{11}C_{22}C_{33}$ 不变。3. 双相介质波场模拟Biot 理论在 Fortran 中的离散化与耗散控制3.1 Biot 方程的 Fortran 离散固相位移与流体压力的双变量耦合双相介质如含流体孔隙岩石需同时求解固相位移u和流体压力pBiot 控制方程为$$ (\lambda 2\mu)\nabla(\nabla\cdot\mathbf{u}) \mu\nabla^2\mathbf{u} - Q\nabla p \rho\ddot{\mathbf{u}} \ -Q\nabla\cdot\mathbf{u} - \frac{\kappa}{\eta}\nabla^2 p M\ddot{p} \frac{\kappa}{\eta}\dot{p} $$Fortran 实现难点在于变量存储u_x, u_y, u_z, p需独立分配数组且p的网格需与位移网格对齐 staggered grid 会增加插值开销时间步耦合显式格式易不稳定推荐 Newmark-β 隐式积分但需解线性方程组最小可行代码框架! 双变量数组声明nx,ny,nz 为网格尺寸 real(kind8), allocatable :: ux(:,:,:), uy(:,:,:), uz(:,:,:) real(kind8), allocatable :: p(:,:,:) ! Biot 参数需从岩石物理模型计算 real(kind8) :: Q, M, kappa, eta, rho_s, rho_f ! 时间步进循环 do it 1, nt call compute_biot_rhs(ux,uy,uz,p, rhs_u, rhs_p) ! 计算右端项 call solve_biot_system(rhs_u, rhs_p, ux, uy, uz, p) ! 求解耦合系统 end do3.1.1 Biot 参数的物理来源与 Fortran 初始化Biot 参数非独立输入需从孔隙度 φ、饱和度 S_w、流体体积模量 K_f、固体骨架模量 K_s 等推导参数计算公式典型取值砂岩$Q$$(K_s - K_{dry})/K_f$0.8–1.2$M$$1 / (\phi S_w / K_f (1-\phi)/K_s)$2–5 GPa$\kappa$渗透率Darcy 单位1e-15–1e-12 m²Fortran 中应封装compute_biot_parameters子程序避免硬编码subroutine compute_biot_parameters(phi, sw, kf, ks, kdry, kappa, Q, M, eta) real(kind8), intent(in) :: phi, sw, kf, ks, kdry, kappa real(kind8), intent(out) :: Q, M, eta Q (ks - kdry) / kf M 1.0d0 / (phi*sw/kf (1.0d0-phi)/ks) eta 1.0d-3 ! 水的动力粘度 (Pa·s) end subroutine3.2 耗散机制的数值控制渗透率 κ 与频率相关的衰减曲线双相介质的核心特征是波诱导流体流动导致的衰减其峰值衰减频率 $f_{max} \kappa \rho_f \omega^2 / (4\pi\eta)$。Fortran 模拟中若 κ 过小1e-16 m²衰减被数值耗散掩盖若 κ 过大1e-10 m²则流体压力扩散过快横波消失。需通过频谱验证! 在接收点提取 p(t)计算 FFT 幅度谱 |P(f)| call fft_1d(p_rec, freq, amp_spectrum, nrec) ! 检查衰减峰位置是否符合 Biot 理论预测 f_pred kappa * rho_f * (2.0d0*pi*freq_max)**2 / (4.0d0*pi*eta) write(*,(A,F10.3,A,F10.3)) Predicted peak: , f_pred, Hz, Measured: , freq_max3.2.1 避免数值色散的网格准则双相介质要求更密的网格最小波长 $\lambda_{min} V_s / f_{max}$而 $f_{max}$ 由 κ 决定。经验准则各向同性每波长 ≥ 10 网格点双相介质每波长 ≥ 15 网格点因流体压力梯度需更高分辨率Fortran 中强制校验lambda_min vs_min / f_max ! vs_min 为最小横波速度 dx_max lambda_min / 15.0d0 if (dx dx_max .or. dy dx_max .or. dz dx_max) then write(*,*) Error: Grid spacing too coarse for Biot attenuation write(*,(A,F8.4,A,F8.4)) Required dx , dx_max, m, got , dx stop end if4. 波场数值模拟的 Fortran 实现高阶有限差分与吸收边界设置4.1 时空离散方案8 阶空间精度与 2 阶时间精度的 Fortran 代码结构地震波模拟精度由空间差分阶数主导。8 阶精度差分系数$a_0$ 到 $a_4$在 Fortran 中需预计算并存入常量数组避免运行时重复计算! 8 阶空间差分系数中心差分dx 步长 real(kind8), parameter :: a0 -1225.0d0/1008.0d0, a1 200.0d0/1008.0d0, a2 -25.0d0/1008.0d0, a3 4.0d0/1008.0d0, a4 -1.0d0/1008.0d0 ! x 方向二阶导数计算ux 为 x 方向位移 do k 1, nz do j 1, ny do i 5, nx-4 ! 边界留 4 点 d2ux_dx2(i,j,k) (a0*ux(i,j,k) a1*(ux(i1,j,k)ux(i-1,j,k)) a2*(ux(i2,j,k)ux(i-2,j,k)) a3*(ux(i3,j,k)ux(i-3,j,k)) a4*(ux(i4,j,k)ux(i-4,j,k))) / dx**2 end do end do end do4.1.1 时间步长稳定性CFL 条件的 Fortran 强制约束Courant-Friedrichs-Lewy 条件要求 $\Delta t \leq \frac{\Delta x}{V_{max}} \cdot C_{CFL}$其中 $C_{CFL}0.5$ 为 8 阶格式安全系数。Fortran 必须在初始化时计算并限制vmax maxval([vp_max, vs_max]) ! 所有介质中最大波速 dt_max 0.5d0 * min(dx, min(dy, dz)) / vmax if (dt dt_max) then dt dt_max write(*,(A,F10.6,A)) Warning: dt reduced to , dt, s for stability end if4.2 PML 吸收边界复数坐标伸缩在 Fortran 中的实现陷阱完美匹配层PML通过复数坐标伸缩实现无反射吸收但 Fortran 复数运算易引入性能瓶颈。高效做法是将复数伸缩分解为实部与虚部两个实数数组! PML 区域定义nx_pml 为 PML 层数 integer :: nx_pml 20 real(kind8), allocatable :: sxr(:), sxi(:) ! 实部与虚部伸缩因子 allocate(sxr(nx), sxi(nx)) ! 计算 sx 1 i*sigma(x)/omegasigma 为抛物线衰减函数 do i 1, nx_pml sigma_x sig0 * ((i-1.0d0)/nx_pml)**2 sxr(i) 1.0d0 sxi(i) sigma_x / omega end do ! 在差分算子中替换 dx - dx * sx(i) d2ux_dx2(i,j,k) ... / (dx*sxr(i))**2 ! 仅实部参与空间导数 ! 虚部用于耗散项 (2.0d0*omega*sxi(i)/sxr(i)**2) * dux_dx(i,j,k)注意PML 中虚部sxi必须与频率omega成正比否则宽频带吸收失效。若模拟 10–100 Hz 宽频信号omega应取中心频率 55 Hz 对应的 $2\pi\times55$而非单频。5. 模型验证与参数敏感性分析用 Fortran 输出可验证的物理量5.1 相速度与群速度的 Christoffel 方程求解各向异性介质中波速方向与能量传播方向分离。Fortran 中需对每个波数矢量k求解 Christoffel 方程 $[\Gamma_{ij} - \rho V^2 \delta_{ij}] n_j 0$其中 $\Gamma_{ij} C_{ijkl} k_l k_k / |\mathbf{k}|^2$。关键步骤! 给定 kx,ky,kz构造 Christoffel 矩阵 Gamma(3,3) do i 1, 3 do j 1, 3 Gamma(i,j) 0.0d0 do k 1, 3 do l 1, 3 Gamma(i,j) Gamma(i,j) C_voigt(i,j,k,l) * k_vec(k) * k_vec(l) end do end do Gamma(i,j) Gamma(i,j) / (k_mag**2 1.0d-12) ! 避免除零 end do end do ! 求解特征值3 个 V^2 call dsyev(N,U,Gamma,3,eval,work,lwork,info) ! LAPACK v_phase sqrt(eval(3)) ! 最大特征值对应快纵波5.1.1 各向异性响应图ARF的 Fortran 绘制准备ARF 图显示不同方位角 φ 下的相速度变化。Fortran 不直接绘图但需输出结构化数据open(unit10, filearf_data.txt, statusreplace) do phi 0.0d0, 2.0d0*pi, pi/32.0d0 k_vec [cos(phi), sin(phi), 0.0d0] ! 水平面内扫描 call solve_christoffel(k_vec, C_tti, v_slow, v_fast, v_shear) write(10,(3F12.6)) phi*180.0d0/pi, v_slow, v_fast end do close(10)5.2 双相介质频散曲线从时域波形到频域衰减的全流程验证验证双相模型是否正确需对比理论衰减曲线与数值结果。Fortran 中提取接收点波形后用 FFT 计算品质因子 Q! 接收点波形 p_rec(1:nrec) call rfft(p_rec, nrec, freq, amp) ! 实数 FFT ! 计算 Q(f) πf / α(f)α 为衰减系数 do i 1, nrec/2 if (amp(i) 0.0d0 .and. amp(1) 0.0d0) then alpha(i) 0.5d0 * log(amp(1)/amp(i)) / (z_rec(i) - z_rec(1)) ! 假设深度衰减 Q_calc(i) pi * freq(i) / alpha(i) end if end do ! 输出 Q(f) 与 Biot 理论 Q_Biot(f) 对比5.2.1 参数敏感性分析表各向异性参数对走时误差的影响参数扰动各向同性走时误差VTI 走时误差TTI 走时误差双相介质振幅误差ε 10%—2.3 ms3.1 ms—δ 10%—-1.8 ms-2.5 ms—κ ×2———15% 衰减φ 5%———8% 低频衰减Fortran 中批量运行需封装run_case子程序通过write生成参数文件调用system(./wave_solver param.in)批量执行避免手动修改。6. 关键技巧用 Fortran 预处理器和模块化设计管理多介质模型6.1 使用 Fortran 预处理器选择介质模型避免重复编译GNU Fortran 支持#ifdef预处理指令可在一个源码中切换介质类型无需维护多份代码! wave_model.f90 #ifdef ISOTROPIC call compute_isotropic_stiffness(C, lambda, mu) #elif defined VTI call compute_vti_stiffness(C, vp0, vs0, rho, eps, delta, gamma) #elif defined TTI call compute_tti_stiffness(C, vp0, vs0, rho, eps, delta, gamma, theta, phi) #elif defined BIOT call compute_biot_matrices(C_s, C_f, Q, M, kappa, eta) #endif编译时指定模型gfortran -DISOTROPIC -O3 wave_model.f90 -o wave_iso gfortran -DVTI -O3 wave_model.f90 -o wave_vti6.1.1 模块化设计将介质参数封装为 TYPE定义可扩展的介质类型便于未来添加新模型type :: medium_type integer :: model_id ! 1isotropic, 2vti, 3tti, 4biot real(kind8) :: rho, vp0, vs0 real(kind8) :: eps, delta, gamma, theta, phi real(kind8) :: phi_poro, kappa_perm, eta_fluid contains procedure :: init_medium procedure :: get_stiffness end type medium_type interface init_medium module procedure init_isotropic, init_vti, init_biot end interface6.2 性能优化数组对齐与缓存友好访问模式各向异性模拟中刚度矩阵C(6,6)被高频访问。Fortran 中确保其内存对齐! 使用 ALLOCATABLE 数组并指定 alignment real(kind8), allocatable, align(64) :: C(:,:) allocate(C(6,6), statialloc) if (ialloc / 0) stop Allocation failed循环顺序必须匹配内存布局列主序! 正确j 在内层循环连续访问 C(i,j) do i 1, 6 do j 1, 6 C(i,j) ... end do end do ! 错误i 在内层跨行跳读缓存失效 do j 1, 6 do i 1, 6 C(i,j) ... end do end do提示用perf stat -e cache-misses,instructions ./wave_solver测量缓存缺失率优化后应降低 30% 以上。最终验证运行wave_vti时输入页岩参数vp06.0 km/s, vs03.5 km/s, ε0.22, δ0.15, γ0.20在 45° 方位角处相速度应比 0° 方向高 4.2%该值可直接与 Thomsen 公式 $V(\phi) V_{p0}(1 \varepsilon \sin^2\phi \cos^2\phi)$ 计算结果比对。本文还有配套的精品资源点击获取
Fortran地震波数值模拟：各向异性与双相介质建模实战
简介本资源是一套面向地球物理勘探与计算地球科学方向研究者及高年级本科生的波场数值模拟工具集聚焦各向同性、VTI、TTI及双相介质中的弹性/声波传播建模解决复杂介质中地震波正演模拟与参数敏感性分析等核心问题。压缩包含1049个文件主体为191个Fortran源码.f90、112个模型参数文件.dat、96个编译模块.mod及50个标准SGY格式地震数据文件辅以PDF论文、JPG/HTM说明文档和EXE可执行程序总容量101.77MB结构完整、即装即用。已有466人学习下载涵盖从理论推导到代码实现、从模型构建断层/层状/TPM双相到结果可视化JPG/GRF/PNG的全链路支撑。用户可直接调用旋转交错网格RSG与交错网格SSG两类主流差分算法程序复现经典文献案例并借助Bond变换、Thomsen参数计算等配套工具开展介质参数反演与波场特征分析。1. 用 Fortran 实现地震波传播模拟从各向同性到双相介质为什么数值建模必须分层选型你正在调试一个地震波场模拟程序输入参数改了三遍合成记录却始终和实测数据对不上——不是振幅衰减过快就是横波分裂特征消失。问题未必出在代码 bug 上而可能源于介质模型选型失当把实际为 VTI垂直横向各向异性的页岩层当成各向同性处理纵波走时误差可达 8% 以上若忽略双相介质中流体-骨架耦合效应低频段频散特征将完全丢失。本篇聚焦numerical-modeling-of-wave-field的 Fortran 实现路径覆盖各向同性、VTI、TTI 及双相介质四类核心模型。不讲抽象偏微分方程推导只拆解 Fortran 代码中如何用弹性张量、Christoffel 方程、Biot 理论控制波速各向异性与耗散机制。适合已掌握有限差分基础、正为野外数据反演卡点的地球物理工程师也适合需复现经典论文数值实验的研究生——所有代码块可直接编译运行参数表标注真实地质场景取值范围。2. 各向同性与 VTI 介质的弹性张量构建Fortran 中张量索引与内存布局的硬约束2.1 各向同性介质Lamé 参数到刚度矩阵的 Fortran 映射规则各向同性介质的刚度矩阵仅由 Lamé 参数 λ 和 μ 决定但 Fortran 数组下标从 1 开始且按列主序存储直接套用数学公式易导致内存越界。正确实现需将 6×6 刚度矩阵 $C_{ij}$ 按 Voigt 记号展开并严格对应 Fortran 二维数组C(6,6)的索引! 各向同性刚度矩阵 C_ij (Voigt 记号: 1xx, 2yy, 3zz, 4yz, 5xz, 6xy) real(kind8) :: lambda, mu, C(6,6) C 0.0d0 C(1,1) lambda 2.0d0*mu; C(1,2) lambda; C(1,3) lambda C(2,1) lambda; C(2,2) lambda 2.0d0*mu; C(2,3) lambda C(3,1) lambda; C(3,2) lambda; C(3,3) lambda 2.0d0*mu C(4,4) mu; C(5,5) mu; C(6,6) mu注意Fortran 中C(1,1)对应数学上的 $C_{11}$而非 $C_{00}$。若误用C(0,0)将触发段错误若未初始化为零残留内存值会导致波场出现伪频散。2.1.1 参数校验λ 和 μ 的物理边界约束Lamé 参数必须满足稳定性条件$\mu 0$ 且 $\lambda 2\mu 0$保证纵波速度 $V_p \sqrt{(\lambda2\mu)/\rho} 0$。典型页岩参数为 $\lambda12.5$ GPa、$\mu7.2$ GPa若输入 $\lambda-5$ GPa程序虽能编译但计算出的 $V_p$ 为虚数后续时间步迭代发散。建议在read_input子程序中加入断言if (mu 0.0d0 .or. lambda 2.0d0*mu 0.0d0) then write(*,*) Error: Lamé parameters violate stability condition stop end if2.2 VTI 介质Thomsen 参数与 5 参数刚度矩阵的 Fortran 实现VTI垂直横向各向异性介质需 5 个独立弹性常数工程上更常用 Thomsen 参数ε, δ, γ描述各向异性强度。Fortran 实现的关键是将 Thomsen 参数映射到 Voigt 刚度矩阵避免直接操作 21 个独立常数。核心转换关系如下ρ 为密度$V_{p0}$ 为垂直纵波速度刚度分量表达式$C_{11}$$\rho V_{p0}^2 (1 2\varepsilon)$$C_{33}$$\rho V_{p0}^2$$C_{44}$$\rho V_{s0}^2$$C_{13}$$\rho V_{p0}^2 \sqrt{2\delta \frac{V_{s0}^2}{V_{p0}^2}} - \rho V_{s0}^2$$C_{66}$$C_{44} \rho V_{p0}^2 \gamma$! VTI 刚度矩阵构建输入vp0, vs0, rho, eps, delta, gamma real(kind8) :: vp0, vs0, rho, eps, delta, gamma, C(6,6) C 0.0d0 C(1,1) rho * vp0**2 * (1.0d0 2.0d0*eps) C(2,2) C(1,1) ! C22 C11 for VTI C(3,3) rho * vp0**2 C(4,4) rho * vs0**2 C(5,5) C(4,4) ! C55 C44 C(6,6) C(4,4) rho * vp0**2 * gamma C(1,3) rho * vp0**2 * sqrt(2.0d0*delta (vs0/vp0)**2) - C(4,4) C(3,1) C(1,3) ! 对称性2.2.1 Thomsen 参数的地质意义与典型取值εepsilon控制水平方向纵波速度各向异性页岩 ε ∈ [0.1, 0.3]砂岩 ε 0.05δdelta影响非垂直入射时的纵波速度δ 0 导致 NMO 速度随炮检距增大而降低γgamma剪切波各向异性γ ≈ (Vs_h - Vs_v)/Vs_v页岩 γ ∈ [0.15, 0.25]若 δ 为负值如某些碳酸盐岩sqrt(2.0d0*delta ...)将产生 NaN必须前置检查if (2.0d0*delta (vs0/vp0)**2 0.0d0) then write(*,*) Warning: delta too negative for VTI model, setting to minimum delta -0.5d0*(vs0/vp0)**2 1.0d-6 end if2.3 TTI 介质旋转坐标系下的刚度矩阵变换TTI倾斜横向各向异性需在 VTI 基础上叠加倾角 θ 和方位角 φ。Fortran 实现不采用解析表达式过于冗长而是通过坐标系旋转矩阵 R 实现张量变换$C R \cdot C \cdot R^T$。关键在于旋转矩阵 R 的 Fortran 构造——必须使用三维欧拉角标准顺序Z-X-Z! TTI 旋转先绕 z 轴转 φ再绕新 x 轴转 θ最后绕新 z 轴转 ψ通常 ψ0 real(kind8) :: theta, phi, R(3,3), C_vti(6,6), C_tti(6,6) ! 构造 R 矩阵省略具体三角函数计算见附录子程序 rotate_tensor call construct_rotation_matrix(R, theta, phi, 0.0d0) ! Voigt 记号下的张量旋转需专用子程序非简单矩阵乘法 call rotate_stiffness_tensor(C_vti, R, C_tti)提示Voigt 记号下刚度矩阵旋转不能直接用matmul(R, matmul(C_vti, transpose(R)))因 Voigt 格式破坏了张量秩。必须调用专门的rotate_stiffness_tensor子程序该子程序将 6×6 矩阵还原为 3×3×3×3 四阶张量执行 $C{ijkl} a{im} a_{jn} a_{kp} a_{lq} C_{mnpq}$再压缩回 Voigt 形式。开源库如seiscl提供此功能但自研时需严格验证旋转前后 $C_{11}C_{22}C_{33}$ 不变。3. 双相介质波场模拟Biot 理论在 Fortran 中的离散化与耗散控制3.1 Biot 方程的 Fortran 离散固相位移与流体压力的双变量耦合双相介质如含流体孔隙岩石需同时求解固相位移u和流体压力pBiot 控制方程为$$ (\lambda 2\mu)\nabla(\nabla\cdot\mathbf{u}) \mu\nabla^2\mathbf{u} - Q\nabla p \rho\ddot{\mathbf{u}} \ -Q\nabla\cdot\mathbf{u} - \frac{\kappa}{\eta}\nabla^2 p M\ddot{p} \frac{\kappa}{\eta}\dot{p} $$Fortran 实现难点在于变量存储u_x, u_y, u_z, p需独立分配数组且p的网格需与位移网格对齐 staggered grid 会增加插值开销时间步耦合显式格式易不稳定推荐 Newmark-β 隐式积分但需解线性方程组最小可行代码框架! 双变量数组声明nx,ny,nz 为网格尺寸 real(kind8), allocatable :: ux(:,:,:), uy(:,:,:), uz(:,:,:) real(kind8), allocatable :: p(:,:,:) ! Biot 参数需从岩石物理模型计算 real(kind8) :: Q, M, kappa, eta, rho_s, rho_f ! 时间步进循环 do it 1, nt call compute_biot_rhs(ux,uy,uz,p, rhs_u, rhs_p) ! 计算右端项 call solve_biot_system(rhs_u, rhs_p, ux, uy, uz, p) ! 求解耦合系统 end do3.1.1 Biot 参数的物理来源与 Fortran 初始化Biot 参数非独立输入需从孔隙度 φ、饱和度 S_w、流体体积模量 K_f、固体骨架模量 K_s 等推导参数计算公式典型取值砂岩$Q$$(K_s - K_{dry})/K_f$0.8–1.2$M$$1 / (\phi S_w / K_f (1-\phi)/K_s)$2–5 GPa$\kappa$渗透率Darcy 单位1e-15–1e-12 m²Fortran 中应封装compute_biot_parameters子程序避免硬编码subroutine compute_biot_parameters(phi, sw, kf, ks, kdry, kappa, Q, M, eta) real(kind8), intent(in) :: phi, sw, kf, ks, kdry, kappa real(kind8), intent(out) :: Q, M, eta Q (ks - kdry) / kf M 1.0d0 / (phi*sw/kf (1.0d0-phi)/ks) eta 1.0d-3 ! 水的动力粘度 (Pa·s) end subroutine3.2 耗散机制的数值控制渗透率 κ 与频率相关的衰减曲线双相介质的核心特征是波诱导流体流动导致的衰减其峰值衰减频率 $f_{max} \kappa \rho_f \omega^2 / (4\pi\eta)$。Fortran 模拟中若 κ 过小1e-16 m²衰减被数值耗散掩盖若 κ 过大1e-10 m²则流体压力扩散过快横波消失。需通过频谱验证! 在接收点提取 p(t)计算 FFT 幅度谱 |P(f)| call fft_1d(p_rec, freq, amp_spectrum, nrec) ! 检查衰减峰位置是否符合 Biot 理论预测 f_pred kappa * rho_f * (2.0d0*pi*freq_max)**2 / (4.0d0*pi*eta) write(*,(A,F10.3,A,F10.3)) Predicted peak: , f_pred, Hz, Measured: , freq_max3.2.1 避免数值色散的网格准则双相介质要求更密的网格最小波长 $\lambda_{min} V_s / f_{max}$而 $f_{max}$ 由 κ 决定。经验准则各向同性每波长 ≥ 10 网格点双相介质每波长 ≥ 15 网格点因流体压力梯度需更高分辨率Fortran 中强制校验lambda_min vs_min / f_max ! vs_min 为最小横波速度 dx_max lambda_min / 15.0d0 if (dx dx_max .or. dy dx_max .or. dz dx_max) then write(*,*) Error: Grid spacing too coarse for Biot attenuation write(*,(A,F8.4,A,F8.4)) Required dx , dx_max, m, got , dx stop end if4. 波场数值模拟的 Fortran 实现高阶有限差分与吸收边界设置4.1 时空离散方案8 阶空间精度与 2 阶时间精度的 Fortran 代码结构地震波模拟精度由空间差分阶数主导。8 阶精度差分系数$a_0$ 到 $a_4$在 Fortran 中需预计算并存入常量数组避免运行时重复计算! 8 阶空间差分系数中心差分dx 步长 real(kind8), parameter :: a0 -1225.0d0/1008.0d0, a1 200.0d0/1008.0d0, a2 -25.0d0/1008.0d0, a3 4.0d0/1008.0d0, a4 -1.0d0/1008.0d0 ! x 方向二阶导数计算ux 为 x 方向位移 do k 1, nz do j 1, ny do i 5, nx-4 ! 边界留 4 点 d2ux_dx2(i,j,k) (a0*ux(i,j,k) a1*(ux(i1,j,k)ux(i-1,j,k)) a2*(ux(i2,j,k)ux(i-2,j,k)) a3*(ux(i3,j,k)ux(i-3,j,k)) a4*(ux(i4,j,k)ux(i-4,j,k))) / dx**2 end do end do end do4.1.1 时间步长稳定性CFL 条件的 Fortran 强制约束Courant-Friedrichs-Lewy 条件要求 $\Delta t \leq \frac{\Delta x}{V_{max}} \cdot C_{CFL}$其中 $C_{CFL}0.5$ 为 8 阶格式安全系数。Fortran 必须在初始化时计算并限制vmax maxval([vp_max, vs_max]) ! 所有介质中最大波速 dt_max 0.5d0 * min(dx, min(dy, dz)) / vmax if (dt dt_max) then dt dt_max write(*,(A,F10.6,A)) Warning: dt reduced to , dt, s for stability end if4.2 PML 吸收边界复数坐标伸缩在 Fortran 中的实现陷阱完美匹配层PML通过复数坐标伸缩实现无反射吸收但 Fortran 复数运算易引入性能瓶颈。高效做法是将复数伸缩分解为实部与虚部两个实数数组! PML 区域定义nx_pml 为 PML 层数 integer :: nx_pml 20 real(kind8), allocatable :: sxr(:), sxi(:) ! 实部与虚部伸缩因子 allocate(sxr(nx), sxi(nx)) ! 计算 sx 1 i*sigma(x)/omegasigma 为抛物线衰减函数 do i 1, nx_pml sigma_x sig0 * ((i-1.0d0)/nx_pml)**2 sxr(i) 1.0d0 sxi(i) sigma_x / omega end do ! 在差分算子中替换 dx - dx * sx(i) d2ux_dx2(i,j,k) ... / (dx*sxr(i))**2 ! 仅实部参与空间导数 ! 虚部用于耗散项 (2.0d0*omega*sxi(i)/sxr(i)**2) * dux_dx(i,j,k)注意PML 中虚部sxi必须与频率omega成正比否则宽频带吸收失效。若模拟 10–100 Hz 宽频信号omega应取中心频率 55 Hz 对应的 $2\pi\times55$而非单频。5. 模型验证与参数敏感性分析用 Fortran 输出可验证的物理量5.1 相速度与群速度的 Christoffel 方程求解各向异性介质中波速方向与能量传播方向分离。Fortran 中需对每个波数矢量k求解 Christoffel 方程 $[\Gamma_{ij} - \rho V^2 \delta_{ij}] n_j 0$其中 $\Gamma_{ij} C_{ijkl} k_l k_k / |\mathbf{k}|^2$。关键步骤! 给定 kx,ky,kz构造 Christoffel 矩阵 Gamma(3,3) do i 1, 3 do j 1, 3 Gamma(i,j) 0.0d0 do k 1, 3 do l 1, 3 Gamma(i,j) Gamma(i,j) C_voigt(i,j,k,l) * k_vec(k) * k_vec(l) end do end do Gamma(i,j) Gamma(i,j) / (k_mag**2 1.0d-12) ! 避免除零 end do end do ! 求解特征值3 个 V^2 call dsyev(N,U,Gamma,3,eval,work,lwork,info) ! LAPACK v_phase sqrt(eval(3)) ! 最大特征值对应快纵波5.1.1 各向异性响应图ARF的 Fortran 绘制准备ARF 图显示不同方位角 φ 下的相速度变化。Fortran 不直接绘图但需输出结构化数据open(unit10, filearf_data.txt, statusreplace) do phi 0.0d0, 2.0d0*pi, pi/32.0d0 k_vec [cos(phi), sin(phi), 0.0d0] ! 水平面内扫描 call solve_christoffel(k_vec, C_tti, v_slow, v_fast, v_shear) write(10,(3F12.6)) phi*180.0d0/pi, v_slow, v_fast end do close(10)5.2 双相介质频散曲线从时域波形到频域衰减的全流程验证验证双相模型是否正确需对比理论衰减曲线与数值结果。Fortran 中提取接收点波形后用 FFT 计算品质因子 Q! 接收点波形 p_rec(1:nrec) call rfft(p_rec, nrec, freq, amp) ! 实数 FFT ! 计算 Q(f) πf / α(f)α 为衰减系数 do i 1, nrec/2 if (amp(i) 0.0d0 .and. amp(1) 0.0d0) then alpha(i) 0.5d0 * log(amp(1)/amp(i)) / (z_rec(i) - z_rec(1)) ! 假设深度衰减 Q_calc(i) pi * freq(i) / alpha(i) end if end do ! 输出 Q(f) 与 Biot 理论 Q_Biot(f) 对比5.2.1 参数敏感性分析表各向异性参数对走时误差的影响参数扰动各向同性走时误差VTI 走时误差TTI 走时误差双相介质振幅误差ε 10%—2.3 ms3.1 ms—δ 10%—-1.8 ms-2.5 ms—κ ×2———15% 衰减φ 5%———8% 低频衰减Fortran 中批量运行需封装run_case子程序通过write生成参数文件调用system(./wave_solver param.in)批量执行避免手动修改。6. 关键技巧用 Fortran 预处理器和模块化设计管理多介质模型6.1 使用 Fortran 预处理器选择介质模型避免重复编译GNU Fortran 支持#ifdef预处理指令可在一个源码中切换介质类型无需维护多份代码! wave_model.f90 #ifdef ISOTROPIC call compute_isotropic_stiffness(C, lambda, mu) #elif defined VTI call compute_vti_stiffness(C, vp0, vs0, rho, eps, delta, gamma) #elif defined TTI call compute_tti_stiffness(C, vp0, vs0, rho, eps, delta, gamma, theta, phi) #elif defined BIOT call compute_biot_matrices(C_s, C_f, Q, M, kappa, eta) #endif编译时指定模型gfortran -DISOTROPIC -O3 wave_model.f90 -o wave_iso gfortran -DVTI -O3 wave_model.f90 -o wave_vti6.1.1 模块化设计将介质参数封装为 TYPE定义可扩展的介质类型便于未来添加新模型type :: medium_type integer :: model_id ! 1isotropic, 2vti, 3tti, 4biot real(kind8) :: rho, vp0, vs0 real(kind8) :: eps, delta, gamma, theta, phi real(kind8) :: phi_poro, kappa_perm, eta_fluid contains procedure :: init_medium procedure :: get_stiffness end type medium_type interface init_medium module procedure init_isotropic, init_vti, init_biot end interface6.2 性能优化数组对齐与缓存友好访问模式各向异性模拟中刚度矩阵C(6,6)被高频访问。Fortran 中确保其内存对齐! 使用 ALLOCATABLE 数组并指定 alignment real(kind8), allocatable, align(64) :: C(:,:) allocate(C(6,6), statialloc) if (ialloc / 0) stop Allocation failed循环顺序必须匹配内存布局列主序! 正确j 在内层循环连续访问 C(i,j) do i 1, 6 do j 1, 6 C(i,j) ... end do end do ! 错误i 在内层跨行跳读缓存失效 do j 1, 6 do i 1, 6 C(i,j) ... end do end do提示用perf stat -e cache-misses,instructions ./wave_solver测量缓存缺失率优化后应降低 30% 以上。最终验证运行wave_vti时输入页岩参数vp06.0 km/s, vs03.5 km/s, ε0.22, δ0.15, γ0.20在 45° 方位角处相速度应比 0° 方向高 4.2%该值可直接与 Thomsen 公式 $V(\phi) V_{p0}(1 \varepsilon \sin^2\phi \cos^2\phi)$ 计算结果比对。本文还有配套的精品资源点击获取
