MATLAB震源震动模拟：从点源模型到合成地震图

简介面向地震动力学与MATLAB仿真初学者的经典案例完整演示震源震动模拟从弹性波方程建立、有限差分/有限元离散到时间步进求解与结果可视化的流程对理解地震波传播规律、评估灾害风险和优化监测网络很有帮助。压缩包内共2个文件包括1个mp4操作录屏和1个.m源码脚本整体约22.49MB。录屏适合对照操作逐步理解建模思路脚本则集中实现了介质参数定义、震源触发、边界条件处理和波动场输出等核心代码可直接运行修改参看不同地层响应。通过该案例可掌握线性动力学方程、Lax-Wendroff与Leap-frog等时间步进算法、自由与无反射边界设定以及surf/slice等绘图工具的使用。资源已有169人学习案例代码结构清晰、注释完整适合地球物理、土木工程等专业学生作为MATLAB建模仿真的入门实践与教学演示辅助资料。1. 震源震动模拟到底在模拟什么把一次地震从发震时刻到台站记录笔尖的完整路径搬进MATLAB很多人第一反应是解波动方程、上有限差分实际上一份能跑的震源震动模拟只做三件事生成震源时间函数、叠加传播路径衰减、合成地面运动记录。这正是地震学里合成地震图的最小闭环也是这类编号案例反复训练的核心能力。对测震学、地震工程方向的师生以及准备MATLAB课设的工程师来说它把矩震级、震源深度、子波主频这些物理参数直接对应到波形幅值和频谱形态上。下面按点源理论、MATLAB实现、参数标定、结果验证这条路线展开每一步都给出可直接复制的脚本。2. 震源模型怎么选点源与有限断层的适用边界一次模拟从选模型开始。震源震动模拟里最容易被跳过的不是代码问题而是点源和有限断层模型各管什么尺度。模型选错后面的参数标定和波形解读都会跟着偏所以先把理论立住再写脚本。2.1 点源模型远场记录的第一近似点源模型的核心假设是断层破裂尺度远小于台站到震源的距离。这种情况可以把整个破裂面等效成一个位错点远场位移的数学形式是震源时间函数的导数与辐射图型因子的乘积再乘上传播路径上的衰减。也就是说只要给定震源时间函数、震级和距离就能算出台站处波形的相对形态。区域地震和近震场景里台站距通常有几十到几百公里而中小地震的破裂尺度只有几百米到几公里点源假设基本成立。有限断层模型则把断层面划分成许多子断层每个子断层有独立的滑移量、破裂时刻和上升时间再通过数值积分把所有子源的贡献叠加起来。它适合震级较大、近场观测丰富的场景计算量比点源大一个量级参数标定也复杂得多。课程案例和绝大多数单台合成记录用点源就够了。判断标准就一条观测距离是否大于断层面最大尺度的三倍满足就点源不满足才考虑有限断层。2.2 雷克子波默认的震源时间函数选择震源时间函数描述断层位错随时间增长的过程。MATLAB相关案例里最常用的震源时间函数是雷克子波表达式为w(t) (1 - 2π²f₀²t²) · exp(-π²f₀²t²)它本质上是高斯函数的二阶导数波形对称、零相位频谱只在主频 f₀ 附近有一个峰。带限特性让合成记录不会凭空出现高频毛刺。主频 f₀ 对应震源辐射能量的优势频率区域地震一般取 15 Hz近场强震动模拟取 510 Hz。子波的有效频带范围大致落在 0.3f₀2.5f₀超出这个范围的能量可以忽略这个近似对后面设置带通滤波很有用。雷克子波的优点是参数少、频带可控缺点是没描述震源破裂的随机性。需要更接近实测特征时换成布龙脉冲它的震源谱在高频按 ω⁻² 衰减更贴近 omega-square 模型。两种时间函数的对比如下特征雷克子波布龙脉冲时间域形态对称、零相位快速上升、缓慢衰减频谱特征单峰带限高频段按 ω⁻² 衰减参数个数仅主频 f0拐角频率与应力降适用场景教学案例、确定性合成记录随机振动、强地面运动合成初学阶段先跑通雷克子波再换布龙脉冲对比频谱差异能直观理解震源时间函数对最终波形的影响。2.3 震源震动模拟涉及的物理量写代码前先列一张参数表避免单位混乱。常见做法是全程以SI单位为主距离用公里、时间用秒、频率用赫兹只在最后输出时统一参数符号单位常用范围子波主频f0Hz110矩震级Mw无量纲3.07.0震源深度hkm550台站距离rkm10300S波速度vskm/s2.53.8介质密度ρkg/m³25002900品质因子Q无量纲100600最容易忽略的是品质因子 Q它控制高频衰减速率Q 越小高频衰减越快。远台记录模拟时 Q 取 150 上下比较常规近场强震动模拟通常取 300 以上。表中的 Mw 与地动位移幅值之间要经过地震矩换算不是直接相乘这一步放到第 4 章专门说明。3. 用MATLAB跑通震源震动模拟的最小流程理论定了开始把脚本跑起来。实现路线是先生成震源时间函数再对时间序列施加几何扩散和介质衰减最后做数值微分得到速度、加速度记录。全程不依赖额外工具箱R2023b 及更早版本都能直接运行。习惯用模块搭建的也可以把震源时间函数放进 Simulink 的 From Workspace 模块再接衰减环节但脚本模式在做参数扫描时更方便后续改动一行就够。3.1 生成雷克子波并落成震源时间函数先写一个生成雷克子波的函数单独保存为 ricker_wavelet.m后面所有脚本都调用它function [t, w] ricker_wavelet(f0, fs, tlen) % 生成雷克子波作为震源时间函数 % 输入: f0 子波主频(Hz), fs 采样率(Hz), tlen 时间窗长度(s) % 输出: t 时间轴(s), w 归一化幅值 t -tlen/2 : 1/fs : tlen/2 - 1/fs; w (1 - 2*(pi*f0*t).^2) .* exp(-(pi*f0*t).^2); w w / max(abs(w)); % 归一化到峰值1幅值交给震级换算环节 end时间轴从 -tlen/2 到 tlen/2是为了让零相位子波的峰值落在 t0这样后续做互相关和走时分析时零点语义清晰。归一化这步把子波峰值压到 1真实记录的幅值大小全部由震级换算环节决定参数职责分离排查问题时不用两头猜。调用方式如下fs 100; % 采样率100 Hz对应奈奎斯特频率50 Hz f0 2; % 主频2 Hz接近区域地震记录的特征频率 tlen 20; % 20秒时间窗覆盖主要能量 [t, w] ricker_wavelet(f0, fs, tlen); figure; plot(t, w, LineWidth, 1.5); xlabel(时间 (s)); ylabel(归一化幅值); title(雷克子波震源时间函数 f02Hz); grid on;运行后应当看到对称的钟形波形峰值在 0 时刻两端平滑趋近于零这就是震源辐射位移的雏形。出图时把背景设为白色用 exportgraphics 导出 pdf 或 eps写论文插图就够用了不需要再截屏。3.2 施加几何扩散与品质因子衰减震源发出的波向外传播振幅不会恒定。体波的几何扩散使振幅按 1/r 衰减介质非弹性又引入随频率变化的衰减 exp(-π f r / (Q v))。把两者叠到子波上就得到台站处的位移记录r 30; % 台站距 30 km vs 3.2; % S波速度 km/s Q 200; % 品质因子 f 0 : 1/tlen : fs - 1/tlen; % 频率轴长度与 w 一致 W fft(w); % 变换到频域施加衰减 att exp(-pi * f * r / (Q * vs)); % 介质衰减因子 W_atten W .* att(:); % 频域逐点相乘 w_geo real(ifft(W_atten)) / r; % 除以距离体现几何扩散几何扩散在时间域里直接除以距离就行介质衰减依赖频率必须在频域逐频点乘衰减因子。att 里频率越高衰减越快这解释了远台记录为什么看起来更“钝”高频成分被优先吃掉。att(:) 用来显式匹配 W 的维度避免 MATLAB 隐式扩展造成维度广播错误。这段代码里的 r 是标量实际做多台站模拟时把 r 换成向量就能一次算出一条测线的记录。3.3 用数值微分从位移推速度和加速度并落盘台站记录的是位移、速度还是加速度取决于仪器类型。有了位移序列用数值微分推速度和加速度vel diff(w_geo) * fs; % 位移对时间求导得到速度 vel [0; vel]; % diff 长度减1用0对齐时间轴 acc diff(vel) * fs; % 速度求导得到加速度 acc [0; acc]; acc bandpass(acc, [0.1 15], fs); % 滤掉数值微分放大出的高频噪声 T table(t(:), w_geo(:), vel(:), acc(:), ... VariableNames, {time_s,disp,vel,acc}); writetable(T, synthetic_station.csv);diff 是差分运算必须乘以采样率 fs 才是物理导数因为离散序列的采样间隔是 1/fs。连续两次求导会把高频噪声成倍放大所以最后用 bandpass 压掉频带外成分0.115 Hz 这个窗口适合区域地震记录。落盘成 CSV 之后后续用 readmatrix 重新读进来做 FFT 频谱分析或者交给其他语言二次处理都不需要再依赖原始脚本。4. 震源模拟的参数标定与三个高频踩坑点代码能跑不等于结果可信。震源震动模拟里参数设错造成的偏差往往比代码 bug 更大而且这类问题书上很少写明。下面三个坑按出现频率排序每一个都会直接污染合成波形。4.1 采样率与子波主频的匹配关系最常踩的坑是采样率设太低。奈奎斯特频率是 fs/2只有高于子波最高有效频率才能保留波形形态。采样定理要求 fs 至少是最高频率的两倍但雷克子波频谱是连续谱最高有效频率约为主频的 3 倍所以工程经验取 fs ≥ 10·f0。也就是说 f02 Hz 的震源至少用 20 Hz 采样稳妥起见直接取 100 Hz这也是前文代码的取值依据。fs 不足时子波峰值会被削平波形出现类似拍频的伪振荡。诊断方法是用 fft 看合成记录的高频端是否有能量被折叠回低频段出现异常尖峰就优先检查采样率。注意这个折算只针对雷克子波这类带限信号如果换用布龙脉冲其高频段衰减较慢需要的采样率还要更高取 fs ≥ 20·f_c 起步。4.2 矩震级与地震矩的换算幅值不能拍脑袋另一个高频错误是直接用 Mw 数值乘子波当作幅值。Mw 是对数标度不能参与线性运算必须先换算成地震矩 M0。SI 单位下的关系是Mw (2/3)·log10(M0) - 6.07其中 M0 单位为 N·m常数项依赖单位制换成 dyn·cm 时要改成 10.7。反算 M0 并估算远场位移量级的代码如下Mw 5.0; M0 10^(1.5 * Mw 9.1); % 由矩震级反算地震矩单位 N·m rho 2700; % 地壳密度 kg/m^3 c 3500; % P波速度 m/s R 30e3; % 台站距这里必须换算成米 d0 M0 / (4 * pi * rho * c^3 * R); % 远场位移峰值量级估算提示同一个距离在几何扩散里用公里、在峰值估算里用米是这类脚本最常见的隐蔽错误建议在文件头部统一声明单位换算关系。远场位移估算来自点源位错理论ρc³ 是介质阻抗贡献。Mw5.0、距离 30 km 代入后d0 在 10⁻³ m 量级附近和实际近震记录到的峰值位移量级一致。注意代码里 R 已经是米而前面的几何扩散用的是公里混用会让结果差三个数量级。4.3 时间窗截断与频谱泄漏第三个坑出现在时间窗不够长、窗两端不归零时。雷克子波在时间窗两端如果仍有非零值fft 会把截断当成周期信号处理频谱上出现旁瓣泄漏。解决方法是加锥形窗tw tukeywin(length(w_geo), 0.15); % 15%的锥形边沿 w_taper w_geo .* tw;tukeywin 第二参数为 1 时退化为完全汉宁窗为 0 时不加窗。0.15 是折中值能压掉两端泄漏又不至于损坏有效信号。另一件常被混淆的事是零填充把序列补零到 2 的幂次长度只提高 fft 频谱的插值密度不增加真实频率分辨率分辨率由时间窗实际长度决定。频谱峰变宽时应该加长时间窗而不是单纯补零。5. 用FFT频谱和双台走时验证模拟结果模拟做完要证明结果物理自洽。这里给一个不依赖外部数据的验证套路单台用频谱核对主频双台用走时核对波速。5.1 用 pwelch 核对主频是否对准震源子波合成记录的峰值频率应当与输入的子波主频一致。用 pwelch 估计功率谱并提取峰值[pxx, fvec] pwelch(w_geo, hann(256), 128, 512, fs); [~, idx] max(pxx); fprintf(峰值频率 %.2f Hz输入主频 %.2f Hz\n, fvec(idx), f0);如果两者偏差超过 20%先检查采样率是否满足 10 倍主频的下限再检查带通滤波是否把主频附近的能量削掉。这个检查脚本在每次调参后反复运行比肉眼看图判断可靠也方便写进批处理循环里做参数扫描。5.2 双台延迟验证波速一致性更进一步模拟两个距离已知的台站用互相关系数求到时差反推波速[c, lag] xcorr(rec_st1, rec_st2, coeff); [~, imax] max(c); dt lag(imax) / fs; % 互相关到时差单位秒 v_est (r2 - r1) / dt % 反推视速度单位 km/s反推出的视速度应当落在介质波速附近偏差超过 15% 时回头检查衰减因子里的波速和双台距离定义单位不统一是这类偏差最常见的根因。本文还有配套的精品资源点击获取