基于MATLAB的TTI介质弹性波正演模拟与震源实现解析

简介这是MATLAB建模仿真32个经典案例中的震源震动模拟项目面向地质学、地球物理学、土木工程等领域研究人员与学生目标是帮助读者掌握从地震波动方程到数值仿真结果的全链路建模方法。资源包总共2个文件包含1个mp4演示视频和1个m源码脚本压缩包大小22.49MB。视频展示完整的操作流程便于跟着复现m脚本基于弹性波方程与有限差分思想给出可复用的地震波传播仿真程序其中涉及线性动力学方程离散化、时间步进算法、无反射边界条件等核心知识点也方便用户修改介质参数或震源位置进一步拓展自己的模型。已有169人学习。通过学习该案例可掌握如何把地震波传播问题转化为可运行的MATLAB程序学会用surf、slice等函数呈现波场演化并理解模拟结果对地震灾害评估、监测台网布局的支撑价值。1. 从psv2dTTI.m说起震源震动模拟为什么难在介质描述解压这个案例压缩包里面最核心的东西其实只有一个psv2dTTI.m和一个演示视频。文件名里的PSV和TTI两个缩写基本已经把这次模拟的技术路线写死了波场是P波与SV波耦合传播介质是倾斜横向各向同性介质。做过各向异性正演的人都知道从各向同性跨到TTI真正麻烦的不是多写几行差分而是弹性刚度矩阵要随对称轴旋转震源的激发方式也得重新约定处理不当就会在波场里看到不该出现的伪S波。这个案例的价值在于它把震源、介质旋转、差分内核三个环节串在一起。适合正在做地震波正演、微震监测或各向异性介质建模的工程师也适合打算从各向同性模拟转向各向异性模拟的研究生。读代码时我建议重点关注两个地方震源是怎么注入的旋转后的弹性常数又是如何进入应力更新的。这两个位置看懂整个脚本的骨架也就清楚了。2. 震源模型的数学约定与Ricker子波的MATLAB实现2.1 Ricker子波参数主频、延迟与采样间隔的标定震源时间函数在弹性波正演里用Ricker子波是普遍做法原因不只是波形简单而是它的频谱集中在主频f0附近且没有零频分量方便预先控制波长与网格步长的关系。Ricker子波的常用形式为w(t) (1 - 2τ²)exp(-τ²)其中 τ πf₀(t - td)td是子波延迟时间通常取1/f0到1.5/f0。td取小了t0时刻子波导数不为零会在波场起始处引入一个不自然的冲击取大了则白白增加前段计算时间。实现时先把整个波形预生成主循环里按时间步索引取值代码如下function w ricker_wavelet(f0, dt, nt) % Ricker子波返回长度nt的行向量 % f0: 主频(Hz), dt: 时间步长(s), nt: 总时间步数 t (0:nt-1) * dt; td 1.2 / f0; % 延迟时间子波峰值出现在td附近 tau pi * f0 * (t - td); w (1 - 2 * tau.^2) .* exp(-tau.^2); end主频f0是第一个需要定下来的参数。它决定了最小波长 λ_min v_min / f0而空间步长必须满足 dx ≤ v_min / (8 * f0)否则波场里会出现肉眼可见的高频频散。反过来f0越高分辨效果越好但网格和时间步长的成本同步上升。二维教学模型里f0常取20到30Hz三维问题则普遍降到10到15Hz。2.2 TTI介质里的震源注入为什么不能只用单点压力各向同性介质中一个同时加在σxx和σzz上的等压源就能激发纯P波。到了TTI介质情况变了爆炸源产生的应力状态不再与介质本征方向对齐一部分能量会耦合进SV模式波场快照里出现一个不该存在的近似椭圆伪S波前。这不是差分格式的误差而是源项与介质本构不匹配导致的物理伪影。处理思路分两层。在弱各向异性介质里可以继续用压力源但接受少量伪波更稳妥的做法是只向对角线应力注入子波τxz始终保持不激励。工程上更常见的是在应力时间导数右端直接加源项等价于给介质一个受迫扰动。代码片段如下% 时间主循环内应力更新后注入震源 tau_xx(2:end-1, 2:end-1) tau_xx(2:end-1, 2:end-1) ... dt .* w(it) .* sxx_spatial; tau_zz(2:end-1, 2:end-1) tau_zz(2:end-1, 2:end-1) ... dt .* w(it) .* szz_spatial; % tau_xz不做源注入避免直接激励剪切模式w(it)是第it步的子波振幅sxx_spatial和szz_spatial是震源在空间网格上的展布矩阵。对τxx和τzz注入同幅值源等价于在介质里制造一个体积膨胀趋势如果介质各向异性较强可以把源向量替换为对应qP波本征方向的权重分布这一步留给后续二次开发。2.3 空间展布与网格采样高斯源半径怎么设单网格点注入在有限差分里会产生严重的空间高频振荡因为点源在波数域近似平坦数值格式会把网格尺度上的寄生分量一并激励起来。常见做法是把震源展布成二维高斯等效于对子波做一次空间低通滤波。高斯半径σ取1.5到2倍网格步长比较合适σ过大源尺度变大短波长信息会被模糊掉。% 构造高斯展布震源中心位于索引(cx, cz)半径sigma [X, Z] meshgrid((1:Nx) * dx, (1:Nz) * dz); gauss exp(-((X - cx * dx).^2 (Z - cz * dz).^2) / (2 * sigma^2)); sxx_spatial gauss / sum(gauss(:)); % 归一化到单位积分 szz_spatial sxx_spatial;X和Z矩阵同时用于震源位置与后续坐标相关参数避免索引与物理坐标错位。表1给出一组推荐起点新手可以按这个组合先跑通再根据波场质量调整。表1常用震源与网格参数起点参数推荐起点调节依据f020 Hz频率越高dx和dt成本越大dx dz10 m保证最小波长内不少于8个网格sigma2*dx抑制点源寄生振荡td1.2/f0保证起始时刻子波趋近于零dt1 ms由第3章的CFL条件约束3. TTI介质弹性波方程的离散化与交错网格差分内核3.1 弹性刚度矩阵旋转Bond变换的MATLAB实现TTI介质可以理解成VTI介质的对称轴在x-z剖面内旋转了一个角度θ。VTI介质本身只需要5个独立弹性常数C11、C13、C33、C44、C66旋转之后刚度矩阵不再保持VTI的稀疏结构C15、C35这些耦合项变为非零。它们直接连接正应力与剪应变是PSV耦合模拟里不能省掉的项。Bond变换是把四阶弹性张量从晶格坐标系投影到全局坐标系的标准做法。二维x-z剖面内绕y轴旋转θ角的6x6形式为C_TTI M C_VTI Mᵀ其中M包含cosθ、sinθ和sin2θ的组合项。注意θ必须使用弧度单位sin(2theta)不要拆成2sin(theta)*cos(theta)再带入矩阵直接写反而减少出错机会function C_tti bond_tilt(C_vti, theta) % 将VTI弹性矩阵绕y轴旋转theta弧度返回TTI弹性矩阵 c cos(theta); s sin(theta); M [c^2 s^2 0 0 0 -sin(2*theta); s^2 c^2 0 0 0 sin(2*theta); 0 0 1 0 0 0; 0 0 0 c s 0; 0 0 0 -s c 0; s*c -s*c 0 0 0 c^2 - s^2]; C_tti M * C_vti * M; end实际调用时先装配VTI矩阵。注意C11是水平方向模量C33是垂直方向模量两者差越大各向异性越强C11 32e9; C13 9e9; C33 26e9; C44 7e9; C66 9e9; rho 2600; % 密度 kg/m^3 C_vti [C11 C13 C13 0 0 0; C13 C11 C13 0 0 0; C13 C13 C33 0 0 0; 0 0 0 C44 0 0; 0 0 0 0 C44 0; 0 0 0 0 0 C66]; theta 30 * pi / 180; C bond_tilt(C_vti, theta); % 提取PSV模拟需要的6个弹性分量 C11c C(1,1); C13c C(1,3); C15c C(1,5); C33c C(3,3); C35c C(3,5); C55c C(5,5);检查旋转结果有个简单办法把θ设为0C应当与C_vti完全一致再把θ设为90度C33c应等于原C11C11c应等于原C33。这条互换关系不成立时问题多半在Bond矩阵中sin2θ项的符号。提示调试TTI代码时建议在Bond变换返回后立即检查theta0时C15与C35是否为零。这一步能省下大量排查时间。3.2 一阶速度-应力方程与中心差分离散二维TTI介质中的PSV波场用一阶速度-应力方程组描述未知量为质点速度vx、vz和应力τxx、τzz、τxz。TTI与VTI在方程形式上的差别集中在应力时间导数旋转后的刚度矩阵里C15、C35非零导致正应力更新必须同时用到∂vx/∂z和∂vz/∂x。对教学型代码把波场放在同一网格上做中心差分是比较直观的写法速度为% 速度更新rho(i,j)为密度场 vx(i,j) vx(i,j) dt / rho(i,j) * ... ((tau_xx(i1,j) - tau_xx(i-1,j)) / (2*dx) ... (tau_xz(i,j1) - tau_xz(i,j-1)) / (2*dz)); vz(i,j) vz(i,j) dt / rho(i,j) * ... ((tau_xz(i1,j) - tau_xz(i-1,j)) / (2*dx) ... (tau_zz(i,j1) - tau_zz(i,j-1)) / (2*dz));应力更新的关键是先把四个梯度一次性算出来再更新三个应力分量避免重复差分% 应力更新TTI介质C15c与C35c耦合项必须保留 dudx (vx(i1,j) - vx(i-1,j)) / (2*dx); dvdz (vz(i,j1) - vz(i,j-1)) / (2*dz); dvdx (vz(i1,j) - vz(i-1,j)) / (2*dx); dudz (vx(i,j1) - vx(i,j-1)) / (2*dz); tau_xx(i,j) tau_xx(i,j) dt * ... (C11c*dudx C13c*dvdz C15c*(dudz dvdx)); tau_zz(i,j) tau_zz(i,j) dt * ... (C13c*dudx C33c*dvdz C35c*(dudz dvdx)); tau_xz(i,j) tau_xz(i,j) dt * ... (C15c*dudx C35c*dvdz C55c*(dudz dvdx));dudx和dvdz对应正应力对正应变的贡献dudzdvdx对应工程剪应变。注意C15c与C35c把剪应变耦合进正应力这正是TTI介质与VTI介质的根本差异。调试时把theta设为0C15c和C35c必然为0应力更新退化为VTI格式观察这一点可以确认Bond变换部分没有写错。3.3 网格步长与CFL条畋稳定性条件是有限差分正演的第一步。对中心差分的二维速度-应力方程时间步长约满足dt ≤ dx / (v_max · √2)v_max取整个模型中可能出现的最大相速度通常用 sqrt(max(C11c, C33c) / rho_min) 估计。频散控制则要求每个最小波长内至少8个网格点即 dx ≤ v_min / (8·f0)v_min取最小剪切波速 sqrt(C55c / rho_max)。两个约束都要满足表2给出一组可直接运行的参数组合。表2均匀TTI模型模拟参数参数数值说明rho2600 kg/m³均匀密度v_p(0°)3162 m/s由sqrt(C33/rho)计算v_p(90°)3508 m/s由sqrt(C11/rho)计算f020 HzRicker子波主频dx dz10 m约等于最小波长/8dt1 msCFL上限约2.02 ms这组参数下模拟600×600网格、2000时间步大概耗时几十秒到几分钟取决于CPU。如果换成交错网格加高阶差分同样精度下可以把空间步长放宽到15m左右但实现复杂度会上升一个台阶。4. 边界处理与波场可视化让仿真结果可解释4.1 吸收衰减带简单可靠的边界处理波场到达模型边界时如果不做处理会反射回计算区域干扰后续波场。PML是完全匹配层效果好但实现复杂、参数多对教学型代码衰减带是更常见的选择在边界内部设置一段吸收带每一步对波场乘一个随空间变化的衰减系数。衰减系数要满足两点从内部到边界逐渐增大避免阻抗突变形成新的反射在边界处足够强让波场在离开计算域前衰减到可忽略。二次曲线剖面是常用形式function damp damping_profile(n, width, strength) % 生成一维衰减系数剖面返回长度n的列向量 % width: 吸收带厚度(网格数), strength: 边界最大衰减强度 damp zeros(n, 1); for k 1 : width r (width - k 1) / width; % 内部为0边界为1 damp(k) strength * r^2; damp(n - k 1) strength * r^2; end end把x和z两个方向的剖面扩成二维衰减场并在每次波场更新后乘一次nbound 30; % 吸收带厚度 strength 5.0; % 边界最大衰减强度 damp_x damping_profile(Nx, nbound, strength); damp_z damping_profile(Nz, nbound, strength); damp2d damp_x damp_z; % Nx x Nz % 主循环内波场更新完毕后统一衰减 vx vx .* exp(-damp2d); vz vz .* exp(-damp2d); tau_xx tau_xx .* exp(-damp2d); tau_zz tau_zz .* exp(-damp2d); tau_xz tau_xz .* exp(-damp2d);strength取4到6宽度取30到50。强度过大或宽度过窄衰减带本身会变成强散射体强度过小则反射波仍能传回中心区域。调试时可以先用各向同性模型检查四个边界的反射是否肉眼可见再逐步调整参数。4.2 波场快照与检波器记录模拟过程中需要定期保存波场快照。常见做法是每隔snap步保存一次vx、vz到.mat文件后处理时再读回。600×600网格的双精度浮点数组单次约2.9MB对现代机器没有压力。保存代码if mod(it, snap) 0 snapfile sprintf(snap_it%05d.mat, it); save(snapfile, vx, vz, it, dt, dx, dz); end检波器记录则是在固定位置持续采样输出模拟地震记录。在时间循环内把接收点位置的vz写入矩阵就得到常说的炮集记录rec_x 200 : 50 : 500; % 检波器水平位置(m) rec_x_idx round(rec_x / dx); rec_z_idx 5; % 距顶边界5个网格 seis(it, :) vz(rec_z_idx, rec_x_idx);记录时机放在速度更新之后、吸收衰减之前保证拿到的是物理波场而不是被边界系数削弱的版本。得到的seis矩阵用imagesc绘制即可查看同相轴形态。4.3 波场可视化与对称色标vz分量的波场快照通常用蓝色-白色-红色对称色标显示零振幅为白色正负振幅对比明显。MATLAB自带的jet色标做出版图层次不够可以自己构造一个对称色标cmap [linspace(0,1,128) linspace(0,1,128) ones(128,1); ones(128,1) linspace(1,0,128) linspace(1,0,128)]; figure; imagesc((0:Nx-1)*dx, (0:Nz-1)*dz, vz); axis equal tight; clim max(abs(vz(:))); caxis([-clim clim]); colormap(cmap); colorbar; xlabel(x (m)); ylabel(z (m)); title(sprintf(vz wavefield at t%.3f s, it*dt));绘图时要注意坐标方向与模型z向下的约定一致。如果看到波前轮廓关于震源中心明显不对称先检查坐标轴是否翻转再检查C15、C35的符号。表3给出边界与可视化的参数起点适合作为排错基准。表3边界与可视化参数起点参数推荐起点检查项nbound3050边界反射可见则增大strength46过大则边界处出现散射snap间隔50步约每50ms一张快照rec_z_idx510避开顶边界奇异区5. 参数标定与算例校核避免TTI模拟常见翻车5.1 用各向同性退化测试验证代码正确性拿到震源震动模拟代码第一件事不是调倾角而是先把算例退化到各向同性。取C11C33并令C13C11-2*C44theta设为0此时波场应当退化为标准P波和SV波两个同心圆。检查两个正交方向的P波到时模拟速度与sqrt(C33/rho)的偏差若超过1%问题多半在差分循环或弹性常数装配上。这一项通过之后再打开TTI介质的倾角才有意义。5.2 高频噪点、边界反射、伪S波的成因与参数联调现场调试讲究从现象反推参数。表4列出三种最常见的翻车现象和处置顺序。表4常见异常现象与调整策略现象优先检查调整策略棋盘状高频噪点dx、sigma加密网格或加大高斯半径边界强反射nbound吸收带加宽到50个网格异常伪S波震源注入去掉tau_xz源改用应力源参数联调顺序建议是先固定f0和rho调整dx让频散消失再缩小dt满足CFL最后处理边界。每一步验证一个小目标不要同时改多个参数否则出问题无法定位。5.3 从均匀到分层TTI参数按位置展开把均匀模型扩展成分层模型时不需要改差分内核只要把C11、C33、C13、C44等标量换成二维数组每个网格点取对应层的值。theta同样可以随空间变化用来模拟倾斜地层界面的局部倾角。建议在主循环外预分配C11m、C33m等矩阵循环内直接索引取值避免每个时间步重复做层位判断。这一步做完这个震源震动模拟就从概念验证变成了可以接实际地层参数的可用工具。本文还有配套的精品资源点击获取