简介本资源是面向地球物理专业研究生、地震工程研究人员及勘探技术人员的MATLAB面波反演工具包聚焦于多道面波频散分析MASW与地下剪切波速结构反演这一核心任务。资源包含16个文件15个.m函数脚本1个.dat示例数据总大小仅57KB轻量但功能完整涵盖数据读取、刚度矩阵计算、频散曲线提取与成像、理论曲线生成、正演模拟、非线性反演及多维可视化等关键模块支持从原始地震记录到速度模型输出的全流程处理。内容预览显示其具备冰岛典型火山-沉积复合地层案例的完整测试流程Test_MASWaves_version1.m调用SampleData.dat并提供误差评估misfit、半空间/分层介质建模Ke_halfspace/Ke_layer等专业实现。目前已有429人学习下载适合需要快速部署面波反演、理解频散物理机制或开展教学实验的科研与工程实践者。1. 面波频散不是“画曲线”而是把地表振动信号翻译成地下剪切波速剖面你拿到一段24道地震检波器记录的面波数据用MATLAB跑完MASWaves_extract_dispersion_curve.m屏幕上跳出一条光滑的频散曲线——但这条线本身毫无地质意义。真正关键的是它背后隐含的是一组离散的、物理可解释的剪切波速Vs随深度变化的参数组合。MASWaves-version1-07-2017这个包本质不是绘图工具而是一套闭环的正演建模→频散计算→反演求解→模型验证链路。它强制你面对一个现实面波频散曲线是高度非唯一的同一组频散数据可能对应十几种完全不同的Vs剖面而MASWaves通过MASWaves_inversion.m中嵌入的最小二乘优化框架结合MASWaves_theoretical_dispersion_curve.m对半无限空间层状介质的精确频散正演把这种模糊性压缩到工程可接受的误差带内。适合两类人一是刚接触面波反演的地球物理研究生需要从源码级理解“为什么不能直接拟合频散曲线”二是已有野外采集经验的工程师想跳过商业软件黑箱用可控参数重跑冰岛火山岩区那种高梯度Vs跃变模型。它不处理原始地震计电压信号只接受预处理后的位移/速度时程如SampleData.dat格式这意味着你必须在前序环节完成去噪、道均衡、时间对齐——这点在Test_MASWaves_version1.m里被刻意省略但实际项目中80%的反演失败源于此处。2. 从原始数据到频散图像四步不可跳过的MATLAB预处理链2.1 数据格式与通道校验MASWaves_read_data.m的隐式约束MASWaves要求输入数据为列向量矩阵每列代表一道检波器记录行数为采样点数。SampleData.dat示例中24道×1024点采样率1000Hz道间距2m。关键约束在于时间轴必须严格等间隔MASWaves_read_data.m不进行重采样仅校验diff(t)是否恒定振幅单位需统一为位移m或速度m/s若用加速度需自行积分cumtrapz两次首道与末道的空间坐标必须能推导出线性阵列几何MASWaves_plot_data.m会据此生成距离-时间图。提示若你的数据来自SAC格式需先用rdseed转为ASCII再用load(-ascii)读入严禁用importdata——它会破坏矩阵维度对齐。% 正确加载示例假设SampleData.dat为24列 data load(SampleData.dat); % 直接生成24列矩阵 if size(data,2) ~ 24 error(通道数不匹配期望24道实际%d道, size(data,2)); end % 校验采样率一致性隐含在Test_MASWaves_version1.m中 dt 0.001; % 必须与实际采样间隔一致 t (0:size(data,1)-1) * dt;2.2 频散成像核心MASWaves_dispersion_imaging.m的参数博弈该函数将时域数据转换为频率-相速度二维能量图其质量直接决定后续反演收敛性。核心参数有三组参数名默认值物理意义调整逻辑fmin,fmax1, 50 Hz频率扫描范围冰岛案例需扩展至80Hz火山岩高频响应强但低于5Hz信噪比骤降cmin,cmax100, 500 m/s相速度搜索区间火山岩区设为300–1200 m/s否则漏掉高速基底nf,nc200, 100频率/速度网格密度过密300×150导致内存溢出过疏100×50丢失拐点% 执行频散成像以冰岛火山岩为例 [f_grid, c_grid, image] MASWaves_dispersion_imaging(... data, dt, 2, ... % data:24道矩阵, dt:0.001s, dx:2m 1, 80, ... % fmin1Hz, fmax80Hz 300, 1200, ... % cmin300m/s, cmax1200m/s 250, 120); % nf250, nc120 % 关键输出image为250×120矩阵每点(i,j)对应频率f_grid(i)与速度c_grid(j)的能量值2.2.1 能量计算原理相位差法 vs. 功率谱比法MASWaves_dispersion_imaging.m默认采用相位差法Phase Difference Method对每对相邻道计算互谱相位除以道间距得相速度。这比功率谱比法Power Spectrum Ratio抗噪性更强但要求道间相干性0.7。当image中出现大面积低能量区值0.1需检查是否存在某道振幅异常用MASWaves_plot_data.m逐道查看阵列是否弯曲dx参数应为平均道距非标称值fmin是否过低导致长周期噪声淹没信号。2.3 频散曲线提取MASWaves_extract_dispersion_curve.m的阈值陷阱该函数在image上搜索局部最大值生成(f,c)点集。但默认阈值thres0.3归一化能量在复杂地质中常失效火山岩区高频段能量衰减快thres0.3会截断有效高频点沉积层区低频段能量集中thres0.3可能合并多个模式。解决方案是分频段动态阈值% 分频段提取冰岛案例实测有效 f_low 1:2:20; % 1–20Hz步长2Hz f_high 22:4:80; % 22–80Hz步长4Hz f_all [f_low, f_high]; c_curve zeros(length(f_all), 1); for k 1:length(f_all) f_idx find(abs(f_grid - f_all(k)) min(abs(f_grid - f_all(k))), 1); % 在f_idx行找最大能量对应的速度索引 [~, c_idx] max(image(f_idx, :)); c_curve(k) c_grid(c_idx); end % 输出c_curve即为频散曲线长度与f_all一致注意MASWaves_extract_dispersion_curve.m内部使用imregionalmax检测峰值若image存在条纹噪声常见于仪器谐波需先用imgaussfilt(image, 2)平滑否则提取点呈锯齿状。3. 反演引擎拆解MASWaves_inversion.m中的三层参数控制3.1 正演模型选择半空间 vs. 层状介质的物理边界MASWaves_inversion.m调用MASWaves_theoretical_dispersion_curve.m计算理论频散而后者依赖两个核心子函数MASWaves_Ke_halfspace.m计算半无限空间无基底的频散适用于松散沉积层MASWaves_Ke_layer.m计算N层介质含基底的频散冰岛案例必须用此函数。关键区别在于半空间模型假设Vs随深度单调递增而层状模型允许Vs跃变如玄武岩盖层→安山岩基底。MASWaves_inversion.m通过nlayer参数切换当nlayer1时自动调用半空间nlayer1则强制层状。若误设nlayer1处理火山岩数据反演结果会出现虚假的渐变过渡层。% 冰岛案例反演配置3层模型 nlayer 3; % 必须≥2才能启用层状正演 initial_model [ % 列向量[厚度1; Vs1; 厚度2; Vs2; Vs3] 15; 350; % 第1层厚15mVs350m/s风化层 25; 720; % 第2层厚25mVs720m/s玄武岩 1200]; % 第3层半无限基底Vs1200m/s深部安山岩 % 注意Vs3无厚度参数由程序自动设为Inf3.2 反演算法参数options结构体的实战调优MASWaves_inversion.m接受options结构体控制优化过程其中三个参数决定成败字段默认值作用冰岛案例建议值MaxIter50最大迭代次数设为100火山岩收敛慢TolFun1e-4目标函数容差改为5e-5避免早停Jacobianfinite-difference雅可比矩阵计算方式保持默认解析雅可比在层状模型中未实现% 构建反演选项 options optimset(MaxIter, 100, TolFun, 5e-5, ... Display, iter, Algorithm, levenberg-marquardt); % 执行反演f_obs/c_obs为实测频散initial_model见上 [best_model, resnorm, residual] MASWaves_inversion(... f_obs, c_obs, nlayer, initial_model, options); % best_model为反演后参数向量需用reshape转为物理模型3.2.1 残差诊断residual向量揭示模型缺陷residual是每个频率点的理论vs实测相速度差单位m/s。若残差绝对值在高频段40Hz持续30m/s说明初始模型Vs3过低基底速度不足需提高initial_model(end)或cmax设置过小导致高频理论曲线被截断。此时应重新运行MASWaves_dispersion_imaging.m扩大cmax至1500m/s。3.3 刚度矩阵计算MASWaves_stiffness_matrix.m的数值稳定性该函数计算层状介质传递矩阵是正演的核心。其稳定性取决于泊松比ν固定为0.25代码中硬编码无法修改。这意味着所有层均按不可压缩介质处理对饱和粘土等ν≈0.45的介质会引入系统偏差厚度参数下限为0.1m若initial_model中某层厚度0.1m程序自动设为0.1m可能导致薄层被忽略。解决方案对含薄互层的沉积序列需手动合并厚度0.5m的层用等效Vs代替。4. 结果验证与可视化用MASWaves_plot_theor_exp_dispersion_curves.m做交叉检验4.1 理论-实测曲线叠置识别模式混淆的关键动作MASWaves_plot_theor_exp_dispersion_curves.m将反演得到的理论频散曲线蓝线与实测点红点绘制在同一图中。但仅看拟合优度R²是危险的——面波存在多模式fundamental mode与higher modes而MASWaves默认只反演基阶模式。若实测点在高频段明显高于理论线大概率是higher mode污染。此时需回查MASWaves_dispersion_imaging.m输出的image确认是否存在第二能量带若存在在MASWaves_extract_dispersion_curve.m中增加mode2参数提取高阶曲线用MASWaves_misfit.m计算双模式联合残差而非单模式。% 双模式验证提取基阶一阶高阶 [f_fund, c_fund] MASWaves_extract_dispersion_curve(image, f_grid, c_grid, 1); [f_h1, c_h1] MASWaves_extract_dispersion_curve(image, f_grid, c_grid, 2); % 分别反演后用同一函数绘制 figure; hold on; plot(f_fund, c_fund, ro, MarkerSize, 4); plot(f_h1, c_h1, go, MarkerSize, 4); % 绘制理论曲线需分别调用 c_theor_fund MASWaves_theoretical_dispersion_curve(f_fund, best_model_fund, nlayer); c_theor_h1 MASWaves_theoretical_dispersion_curve(f_h1, best_model_h1, nlayer); plot(f_fund, c_theor_fund, b-, LineWidth, 1.5); plot(f_h1, c_theor_h1, m-, LineWidth, 1.5); xlabel(Frequency (Hz)); ylabel(Phase Velocity (m/s)); legend(Fundamental Mode Obs, Higher Mode Obs, Fundamental Model, Higher Model);4.2 速度剖面图MASWaves_plot_dispersion_image_2D.m的深度标定技巧该函数生成速度-深度剖面图但纵坐标是层底深度而非中心深度。例如initial_model[15;350;25;720;1200]图中第一层显示为0–15m第二层为15–40m1525第三层为40m以下。若要对比钻孔数据需将钻孔Vs值插值到层底深度点对钻孔深度z_i找到其所在层k满足depth_{k-1} z_i ≤ depth_k用线性插值计算该层内Vs(z_i) Vs_{k-1} (Vs_k - Vs_{k-1}) × (z_i - depth_{k-1}) / thickness_k。提示MASWaves_plot_dispersion_image_2D.m默认y轴反转深度向下增大若需常规坐标执行set(gca,YDir,normal)。5. 冰岛火山岩区实战处理高梯度Vs跃变的三步修正法5.1 问题定位高频残差突增的物理根源在冰岛某火山口边缘采集的24道数据中反演后residual在35–60Hz区间出现45m/s突增理论值低于实测。检查image发现该频段存在两条平行能量带上带为基阶模式下带为一阶高阶模式但MASWaves_extract_dispersion_curve.m因阈值过高仅捕获上带。根本原因是火山岩Vs梯度达150m/s/m导致高阶模式能量显著增强。5.2 修正步骤一动态阈值提取双模式% 对35–60Hz频段单独处理 f_target 35:1:60; f_idx_target find(f_grid 35 f_grid 60); % 计算该频段内每行的最大能量值 row_max max(image(f_idx_target, :), [], 2); % 设定双阈值主模式用0.4高阶模式用0.25因能量弱 thres_fund 0.4; thres_h1 0.25; c_fund zeros(length(f_target), 1); c_h1 zeros(length(f_target), 1); for k 1:length(f_target) f_pos find(f_grid f_target(k), 1); if ~isempty(f_pos) f_pos length(f_grid) % 主模式找能量0.4×row_max(f_pos)的最高速度点 mask_fund image(f_pos, :) thres_fund * row_max(k); if any(mask_fund) [~, idx_fund] max(image(f_pos, mask_fund)); c_fund(k) c_grid(find(mask_fund, 1, first) idx_fund - 1); end % 高阶模式找能量0.25×row_max(f_pos)且速度0.8×c_fund(k)的点 mask_h1 image(f_pos, :) thres_h1 * row_max(k) c_grid 0.8*c_fund(k); if any(mask_h1) [~, idx_h1] max(image(f_pos, mask_h1)); c_h1(k) c_grid(find(mask_h1, 1, first) idx_h1 - 1); end end end5.3 修正步骤二分层反演约束Vs梯度将3层模型改为4层强制在25m深度处设置Vs跃变约束% 新增约束第2层底界深度25mVs从720→950跃变 nlayer 4; initial_model [15; 350; 10; 720; 15; 950; 1200]; % 注意厚度参数必须为正且总和覆盖目标深度 % 反演时固定第3层厚度15m对应25–40m仅优化Vs参数 fixed_params [0,0,0,0,1,0,0]; % 1表示该参数固定0表示可变 [best_model, ~, ~] MASWaves_inversion(f_obs, c_obs, nlayer, initial_model, options, fixed_params);5.4 修正步骤三用MASWaves_plot_dispersion_image_3D.m验证空间一致性最后将反演得到的Vs剖面导入MASWaves_plot_dispersion_image_3D.m生成三维频散曲面。若沿测线方向x轴的曲面在火山口位置出现陡峭褶皱而周边平缓则证实高梯度模型成功捕捉了地质突变。此时导出best_model中的Vs值即可作为后续地震危险性评估的输入参数——这才是MASWaves源码交付的终极价值不是一张图而是一组可嵌入工程模型的、经物理方程验证的地下参数。本文还有配套的精品资源点击获取 MASWaves面波反演原理与火山岩区高梯度Vs建模实战
简介本资源是面向地球物理专业研究生、地震工程研究人员及勘探技术人员的MATLAB面波反演工具包聚焦于多道面波频散分析MASW与地下剪切波速结构反演这一核心任务。资源包含16个文件15个.m函数脚本1个.dat示例数据总大小仅57KB轻量但功能完整涵盖数据读取、刚度矩阵计算、频散曲线提取与成像、理论曲线生成、正演模拟、非线性反演及多维可视化等关键模块支持从原始地震记录到速度模型输出的全流程处理。内容预览显示其具备冰岛典型火山-沉积复合地层案例的完整测试流程Test_MASWaves_version1.m调用SampleData.dat并提供误差评估misfit、半空间/分层介质建模Ke_halfspace/Ke_layer等专业实现。目前已有429人学习下载适合需要快速部署面波反演、理解频散物理机制或开展教学实验的科研与工程实践者。1. 面波频散不是“画曲线”而是把地表振动信号翻译成地下剪切波速剖面你拿到一段24道地震检波器记录的面波数据用MATLAB跑完MASWaves_extract_dispersion_curve.m屏幕上跳出一条光滑的频散曲线——但这条线本身毫无地质意义。真正关键的是它背后隐含的是一组离散的、物理可解释的剪切波速Vs随深度变化的参数组合。MASWaves-version1-07-2017这个包本质不是绘图工具而是一套闭环的正演建模→频散计算→反演求解→模型验证链路。它强制你面对一个现实面波频散曲线是高度非唯一的同一组频散数据可能对应十几种完全不同的Vs剖面而MASWaves通过MASWaves_inversion.m中嵌入的最小二乘优化框架结合MASWaves_theoretical_dispersion_curve.m对半无限空间层状介质的精确频散正演把这种模糊性压缩到工程可接受的误差带内。适合两类人一是刚接触面波反演的地球物理研究生需要从源码级理解“为什么不能直接拟合频散曲线”二是已有野外采集经验的工程师想跳过商业软件黑箱用可控参数重跑冰岛火山岩区那种高梯度Vs跃变模型。它不处理原始地震计电压信号只接受预处理后的位移/速度时程如SampleData.dat格式这意味着你必须在前序环节完成去噪、道均衡、时间对齐——这点在Test_MASWaves_version1.m里被刻意省略但实际项目中80%的反演失败源于此处。2. 从原始数据到频散图像四步不可跳过的MATLAB预处理链2.1 数据格式与通道校验MASWaves_read_data.m的隐式约束MASWaves要求输入数据为列向量矩阵每列代表一道检波器记录行数为采样点数。SampleData.dat示例中24道×1024点采样率1000Hz道间距2m。关键约束在于时间轴必须严格等间隔MASWaves_read_data.m不进行重采样仅校验diff(t)是否恒定振幅单位需统一为位移m或速度m/s若用加速度需自行积分cumtrapz两次首道与末道的空间坐标必须能推导出线性阵列几何MASWaves_plot_data.m会据此生成距离-时间图。提示若你的数据来自SAC格式需先用rdseed转为ASCII再用load(-ascii)读入严禁用importdata——它会破坏矩阵维度对齐。% 正确加载示例假设SampleData.dat为24列 data load(SampleData.dat); % 直接生成24列矩阵 if size(data,2) ~ 24 error(通道数不匹配期望24道实际%d道, size(data,2)); end % 校验采样率一致性隐含在Test_MASWaves_version1.m中 dt 0.001; % 必须与实际采样间隔一致 t (0:size(data,1)-1) * dt;2.2 频散成像核心MASWaves_dispersion_imaging.m的参数博弈该函数将时域数据转换为频率-相速度二维能量图其质量直接决定后续反演收敛性。核心参数有三组参数名默认值物理意义调整逻辑fmin,fmax1, 50 Hz频率扫描范围冰岛案例需扩展至80Hz火山岩高频响应强但低于5Hz信噪比骤降cmin,cmax100, 500 m/s相速度搜索区间火山岩区设为300–1200 m/s否则漏掉高速基底nf,nc200, 100频率/速度网格密度过密300×150导致内存溢出过疏100×50丢失拐点% 执行频散成像以冰岛火山岩为例 [f_grid, c_grid, image] MASWaves_dispersion_imaging(... data, dt, 2, ... % data:24道矩阵, dt:0.001s, dx:2m 1, 80, ... % fmin1Hz, fmax80Hz 300, 1200, ... % cmin300m/s, cmax1200m/s 250, 120); % nf250, nc120 % 关键输出image为250×120矩阵每点(i,j)对应频率f_grid(i)与速度c_grid(j)的能量值2.2.1 能量计算原理相位差法 vs. 功率谱比法MASWaves_dispersion_imaging.m默认采用相位差法Phase Difference Method对每对相邻道计算互谱相位除以道间距得相速度。这比功率谱比法Power Spectrum Ratio抗噪性更强但要求道间相干性0.7。当image中出现大面积低能量区值0.1需检查是否存在某道振幅异常用MASWaves_plot_data.m逐道查看阵列是否弯曲dx参数应为平均道距非标称值fmin是否过低导致长周期噪声淹没信号。2.3 频散曲线提取MASWaves_extract_dispersion_curve.m的阈值陷阱该函数在image上搜索局部最大值生成(f,c)点集。但默认阈值thres0.3归一化能量在复杂地质中常失效火山岩区高频段能量衰减快thres0.3会截断有效高频点沉积层区低频段能量集中thres0.3可能合并多个模式。解决方案是分频段动态阈值% 分频段提取冰岛案例实测有效 f_low 1:2:20; % 1–20Hz步长2Hz f_high 22:4:80; % 22–80Hz步长4Hz f_all [f_low, f_high]; c_curve zeros(length(f_all), 1); for k 1:length(f_all) f_idx find(abs(f_grid - f_all(k)) min(abs(f_grid - f_all(k))), 1); % 在f_idx行找最大能量对应的速度索引 [~, c_idx] max(image(f_idx, :)); c_curve(k) c_grid(c_idx); end % 输出c_curve即为频散曲线长度与f_all一致注意MASWaves_extract_dispersion_curve.m内部使用imregionalmax检测峰值若image存在条纹噪声常见于仪器谐波需先用imgaussfilt(image, 2)平滑否则提取点呈锯齿状。3. 反演引擎拆解MASWaves_inversion.m中的三层参数控制3.1 正演模型选择半空间 vs. 层状介质的物理边界MASWaves_inversion.m调用MASWaves_theoretical_dispersion_curve.m计算理论频散而后者依赖两个核心子函数MASWaves_Ke_halfspace.m计算半无限空间无基底的频散适用于松散沉积层MASWaves_Ke_layer.m计算N层介质含基底的频散冰岛案例必须用此函数。关键区别在于半空间模型假设Vs随深度单调递增而层状模型允许Vs跃变如玄武岩盖层→安山岩基底。MASWaves_inversion.m通过nlayer参数切换当nlayer1时自动调用半空间nlayer1则强制层状。若误设nlayer1处理火山岩数据反演结果会出现虚假的渐变过渡层。% 冰岛案例反演配置3层模型 nlayer 3; % 必须≥2才能启用层状正演 initial_model [ % 列向量[厚度1; Vs1; 厚度2; Vs2; Vs3] 15; 350; % 第1层厚15mVs350m/s风化层 25; 720; % 第2层厚25mVs720m/s玄武岩 1200]; % 第3层半无限基底Vs1200m/s深部安山岩 % 注意Vs3无厚度参数由程序自动设为Inf3.2 反演算法参数options结构体的实战调优MASWaves_inversion.m接受options结构体控制优化过程其中三个参数决定成败字段默认值作用冰岛案例建议值MaxIter50最大迭代次数设为100火山岩收敛慢TolFun1e-4目标函数容差改为5e-5避免早停Jacobianfinite-difference雅可比矩阵计算方式保持默认解析雅可比在层状模型中未实现% 构建反演选项 options optimset(MaxIter, 100, TolFun, 5e-5, ... Display, iter, Algorithm, levenberg-marquardt); % 执行反演f_obs/c_obs为实测频散initial_model见上 [best_model, resnorm, residual] MASWaves_inversion(... f_obs, c_obs, nlayer, initial_model, options); % best_model为反演后参数向量需用reshape转为物理模型3.2.1 残差诊断residual向量揭示模型缺陷residual是每个频率点的理论vs实测相速度差单位m/s。若残差绝对值在高频段40Hz持续30m/s说明初始模型Vs3过低基底速度不足需提高initial_model(end)或cmax设置过小导致高频理论曲线被截断。此时应重新运行MASWaves_dispersion_imaging.m扩大cmax至1500m/s。3.3 刚度矩阵计算MASWaves_stiffness_matrix.m的数值稳定性该函数计算层状介质传递矩阵是正演的核心。其稳定性取决于泊松比ν固定为0.25代码中硬编码无法修改。这意味着所有层均按不可压缩介质处理对饱和粘土等ν≈0.45的介质会引入系统偏差厚度参数下限为0.1m若initial_model中某层厚度0.1m程序自动设为0.1m可能导致薄层被忽略。解决方案对含薄互层的沉积序列需手动合并厚度0.5m的层用等效Vs代替。4. 结果验证与可视化用MASWaves_plot_theor_exp_dispersion_curves.m做交叉检验4.1 理论-实测曲线叠置识别模式混淆的关键动作MASWaves_plot_theor_exp_dispersion_curves.m将反演得到的理论频散曲线蓝线与实测点红点绘制在同一图中。但仅看拟合优度R²是危险的——面波存在多模式fundamental mode与higher modes而MASWaves默认只反演基阶模式。若实测点在高频段明显高于理论线大概率是higher mode污染。此时需回查MASWaves_dispersion_imaging.m输出的image确认是否存在第二能量带若存在在MASWaves_extract_dispersion_curve.m中增加mode2参数提取高阶曲线用MASWaves_misfit.m计算双模式联合残差而非单模式。% 双模式验证提取基阶一阶高阶 [f_fund, c_fund] MASWaves_extract_dispersion_curve(image, f_grid, c_grid, 1); [f_h1, c_h1] MASWaves_extract_dispersion_curve(image, f_grid, c_grid, 2); % 分别反演后用同一函数绘制 figure; hold on; plot(f_fund, c_fund, ro, MarkerSize, 4); plot(f_h1, c_h1, go, MarkerSize, 4); % 绘制理论曲线需分别调用 c_theor_fund MASWaves_theoretical_dispersion_curve(f_fund, best_model_fund, nlayer); c_theor_h1 MASWaves_theoretical_dispersion_curve(f_h1, best_model_h1, nlayer); plot(f_fund, c_theor_fund, b-, LineWidth, 1.5); plot(f_h1, c_theor_h1, m-, LineWidth, 1.5); xlabel(Frequency (Hz)); ylabel(Phase Velocity (m/s)); legend(Fundamental Mode Obs, Higher Mode Obs, Fundamental Model, Higher Model);4.2 速度剖面图MASWaves_plot_dispersion_image_2D.m的深度标定技巧该函数生成速度-深度剖面图但纵坐标是层底深度而非中心深度。例如initial_model[15;350;25;720;1200]图中第一层显示为0–15m第二层为15–40m1525第三层为40m以下。若要对比钻孔数据需将钻孔Vs值插值到层底深度点对钻孔深度z_i找到其所在层k满足depth_{k-1} z_i ≤ depth_k用线性插值计算该层内Vs(z_i) Vs_{k-1} (Vs_k - Vs_{k-1}) × (z_i - depth_{k-1}) / thickness_k。提示MASWaves_plot_dispersion_image_2D.m默认y轴反转深度向下增大若需常规坐标执行set(gca,YDir,normal)。5. 冰岛火山岩区实战处理高梯度Vs跃变的三步修正法5.1 问题定位高频残差突增的物理根源在冰岛某火山口边缘采集的24道数据中反演后residual在35–60Hz区间出现45m/s突增理论值低于实测。检查image发现该频段存在两条平行能量带上带为基阶模式下带为一阶高阶模式但MASWaves_extract_dispersion_curve.m因阈值过高仅捕获上带。根本原因是火山岩Vs梯度达150m/s/m导致高阶模式能量显著增强。5.2 修正步骤一动态阈值提取双模式% 对35–60Hz频段单独处理 f_target 35:1:60; f_idx_target find(f_grid 35 f_grid 60); % 计算该频段内每行的最大能量值 row_max max(image(f_idx_target, :), [], 2); % 设定双阈值主模式用0.4高阶模式用0.25因能量弱 thres_fund 0.4; thres_h1 0.25; c_fund zeros(length(f_target), 1); c_h1 zeros(length(f_target), 1); for k 1:length(f_target) f_pos find(f_grid f_target(k), 1); if ~isempty(f_pos) f_pos length(f_grid) % 主模式找能量0.4×row_max(f_pos)的最高速度点 mask_fund image(f_pos, :) thres_fund * row_max(k); if any(mask_fund) [~, idx_fund] max(image(f_pos, mask_fund)); c_fund(k) c_grid(find(mask_fund, 1, first) idx_fund - 1); end % 高阶模式找能量0.25×row_max(f_pos)且速度0.8×c_fund(k)的点 mask_h1 image(f_pos, :) thres_h1 * row_max(k) c_grid 0.8*c_fund(k); if any(mask_h1) [~, idx_h1] max(image(f_pos, mask_h1)); c_h1(k) c_grid(find(mask_h1, 1, first) idx_h1 - 1); end end end5.3 修正步骤二分层反演约束Vs梯度将3层模型改为4层强制在25m深度处设置Vs跃变约束% 新增约束第2层底界深度25mVs从720→950跃变 nlayer 4; initial_model [15; 350; 10; 720; 15; 950; 1200]; % 注意厚度参数必须为正且总和覆盖目标深度 % 反演时固定第3层厚度15m对应25–40m仅优化Vs参数 fixed_params [0,0,0,0,1,0,0]; % 1表示该参数固定0表示可变 [best_model, ~, ~] MASWaves_inversion(f_obs, c_obs, nlayer, initial_model, options, fixed_params);5.4 修正步骤三用MASWaves_plot_dispersion_image_3D.m验证空间一致性最后将反演得到的Vs剖面导入MASWaves_plot_dispersion_image_3D.m生成三维频散曲面。若沿测线方向x轴的曲面在火山口位置出现陡峭褶皱而周边平缓则证实高梯度模型成功捕捉了地质突变。此时导出best_model中的Vs值即可作为后续地震危险性评估的输入参数——这才是MASWaves源码交付的终极价值不是一张图而是一组可嵌入工程模型的、经物理方程验证的地下参数。本文还有配套的精品资源点击获取