简介本资源是一套面向地球物理勘探专业师生及科研人员的MATLAB实操工具聚焦多道面波分析中相移法频散曲线提取这一核心任务解决野外地震记录中面波速度—频率关系建模难、相位解缠易出错等实际问题。压缩包共2个文件1个主程序PhaseShift.m与1个示例数据seis.mat总大小119KB其中MATLAB脚本完整实现数据预处理、FFT相位提取、unwrap相位解缠、相邻道相位差计算及频散曲线自动绘制全流程.mat文件提供真实单炮面波记录供即开即用验证。已有772人学习下载适用于课程实验、毕业设计或科研初期快速复现经典面波反演方法。用户可直接运行脚本通过调整采样率、频率范围等参数适配不同采集系统输出结果包含清晰频散图像与速度-频率数值表便于后续联合反演或地质解释。1. 项目缘起从一道面试题到一套实用工具几年前我还在处理工程物探数据当时面试一位新人我随口问了个问题“给你一条多道面波记录不用商业软件你怎么最快地把频散曲线给提出来” 他愣了一下然后开始讲各种变换和手动拾取。我告诉他其实有个很巧妙的方法叫相移法用MATLAB几十行代码就能实现而且抗噪性不错。后来我把这个思路整理成了程序成了团队内部的一个小工具。今天我就把这个“多道面波分析相移法频散曲线提取方法”的MATLAB实现从原理到代码再到实际处理中的各种坑完整地分享出来。对于搞浅层地震勘探、工程物探或者研究地震波传播的朋友来说面波频散曲线是反演地下横波速度结构的关键输入。传统方法比如f-k谱分析或τ-p变换要么对道间距要求苛刻要么计算量不小。相移法Phase Shift提供了一种相对直观、在频率-波数域直接计算相速度谱的思路特别适合处理常规排列采集的多道数据。这个方法不新鲜但网上能找到的、真正能跑通、附带详细注释和实用技巧的MATLAB代码并不多。本文将手把手带你理解相移法的核心并给你一套可以直接运行、修改的MATLAB程序同时会重点聊聊我在处理实际数据时遇到的波形畸变、噪声干扰和参数选择问题。2. 相移法核心原理为什么是“移相”而不是“变换”理解相移法关键在于跳出“变换”的思维定式。我们最终目标是得到每个频率成分对应的相速度。想象一下对于某个特定的频率f如果有一个平面波以速度v在这个频率下传播那么相邻两个检波器记录到的该频率信号的相位差应该是固定的。2.1 从波动方程到相位移动我们从简谐平面波的表达式出发。一个沿x方向传播、角频率为ω的单频平面波在位置x处的振动可以表示为u(x, t) A * exp(i * (kx - ωt))其中k是波数k ω / v 2πf / v。这里v就是我们要求的相速度。现在假设我们在位置x0处有一个记录道u(x0, t)。如果我们想“猜测”一个测试相速度v_test并计算按照这个速度传播信号在另一个位置x1处“应该”是什么样子我们可以对u(x0, t)进行一个相位移动操作。这个操作在频率域进行极其方便。具体来说对u(x0, t)做傅里叶变换到频率域得到U(x0, f)。那么根据上述平面波公式在位置x1处的波场U(x1, f)理论上应该是U(x1, f) U(x0, f) * exp(i * k * Δx) U(x0, f) * exp(i * 2πf * Δx / v_test)这里的exp(i * 2πf * Δx / v_test)就是一个相位移动因子。它把参考道x0处的频谱移动到了x1处。2.2 多道叠加与“速度谱”的生成相移法的巧妙之处在于它利用了整个排列的所有道。我们不是两两对比而是将所有道都向一个虚拟的“零偏移距”参考点进行相位移动。选择参考点通常选择第一个检波器最小偏移距或排列中心作为参考点x_ref。遍历测试速度对于一个给定的频率f我们预设一个相速度v_test的扫描范围比如从100 m/s到1000 m/s。相位移动与叠加对于排列中的第j个检波器其位置为x_j。我们计算它相对于参考点的距离Δx_j x_j - x_ref。然后将该道在频率f处的频谱U(x_j, f)乘以一个反向的相位移动因子U(x_j, f) * exp(-i * 2πf * Δx_j / v_test)。这个操作相当于把该道的信号“搬回”到参考点位置前提是信号确实是以v_test速度传播的。相干叠加如果实际的相速度恰好等于v_test那么所有道经过上述相位移动后它们在参考点处的“估计信号”的相位将会完全对齐。将所有道的这些“搬回来”的频谱在频率f处求和其幅值将会达到最大。如果v_test不等于真实速度各道相位参差不齐叠加后会相互抵消幅值较小。构建速度谱对每个频率f重复步骤2-4遍历所有v_test计算每个(f, v_test)组合下的叠加幅值。这个二维矩阵频率×速度的幅值就是相速度谱或称频散能量谱。对于每个频率f在速度轴上寻找幅值最大的点其对应的速度v就是该频率的相速度估计值。连接这些点就得到了提取的频散曲线。注意这里的“移相”是概念核心。exp(-i * 2πf * Δx / v)中的负号很关键它表示“将信号从当前位置移回参考点”。如果符号弄反结果将完全错误。3. MATLAB程序实现逐行拆解与关键函数理论清晰后我们来看代码实现。我将程序分为几个核心函数方便理解和调用。3.1 主函数dispersion_phase_shift.m这是程序的入口负责流程控制。function [disp_curve, velocity_spectrum, f_axis, v_axis] dispersion_phase_shift(data, dt, dx, v_min, v_max, dv, f_min, f_max) % 相移法提取面波频散曲线 % 输入 % data - 地震数据矩阵每一列是一个道行是时间采样点 (nt x nr) % dt - 时间采样间隔 (秒) % dx - 道间距 (米) % v_min - 扫描最小相速度 (m/s) % v_max - 扫描最大相速度 (m/s) % dv - 扫描速度间隔 (m/s) % f_min - 分析最小频率 (Hz) % f_max - 分析最大频率 (Hz) % 输出 % disp_curve - 提取的频散曲线两列矩阵 [频率, 相速度] % velocity_spectrum - 相速度谱矩阵 (频率轴 x 速度轴) % f_axis - 频率轴向量 % v_axis - 速度轴向量 [nt, nr] size(data); % nt: 时间点数 nr: 道数 t_axis (0:nt-1)*dt; % 时间轴 % 1. 计算频率轴 NFFT 2^nextpow2(nt); % 使用2的幂次以提高FFT效率 f_axis_full (0:NFFT/2) / (NFFT*dt); % 单边频谱频率轴 % 选取感兴趣的频率范围 f_idx find(f_axis_full f_min f_axis_full f_max); f_axis f_axis_full(f_idx); nf length(f_axis); % 2. 生成速度轴 v_axis v_min:dv:v_max; nv length(v_axis); % 3. 对每一道数据进行FFT并只取感兴趣频率部分 data_fft zeros(NFFT, nr); for i 1:nr data_fft(:, i) fft(data(:, i), NFFT); end data_fft data_fft(1:NFFT/21, :); % 取单边谱 data_fft data_fft(f_idx, :); % 截取频率范围 % 4. 定义检波器位置以第一道为参考点 x_pos (0:nr-1) * dx; % 检波器位置坐标 x_ref x_pos(1); % 参考点设为第一道 delta_x x_pos - x_ref; % 各道相对于参考点的距离 % 5. 初始化速度谱矩阵 velocity_spectrum zeros(nf, nv); % 6. 核心相移计算循环 % 为了提高计算效率我们逐频率计算并对向量化操作进行优化 for f_idx 1:nf f f_axis(f_idx); % 当前频率 U_f data_fft(f_idx, :); % 所有道在当前频率下的频谱值 (1 x nr 向量) for v_idx 1:nv v_test v_axis(v_idx); % 当前测试速度 % 计算相位移动因子向量 (1 x nr) % 注意这里使用矩阵运算避免内层循环 phase_shift_factor exp(-1i * 2 * pi * f * delta_x / v_test); % 将各道频谱移相后叠加 stacked_amplitude abs(sum(U_f .* phase_shift_factor)); % 存储到速度谱中 velocity_spectrum(f_idx, v_idx) stacked_amplitude; end % 可选显示进度对于大数据量很实用 if mod(f_idx, 10) 0 fprintf(Processing frequency %d / %d...\n, f_idx, nf); end end % 7. 从速度谱中提取频散曲线寻找每个频率下的能量峰值 disp_curve zeros(nf, 2); for f_idx 1:nf [~, max_idx] max(velocity_spectrum(f_idx, :)); disp_curve(f_idx, 1) f_axis(f_idx); disp_curve(f_idx, 2) v_axis(max_idx); end % 8. 可选简单的后处理去除明显异常的孤立点 % 例如可以基于速度的局部中值滤波 window_size 5; for i 1:nf start_idx max(1, i - floor(window_size/2)); end_idx min(nf, i floor(window_size/2)); median_v median(disp_curve(start_idx:end_idx, 2)); % 如果当前点速度与局部中值相差过大则用中值替代 if abs(disp_curve(i, 2) - median_v) 0.3 * median_v disp_curve(i, 2) median_v; end end end关键点解析FFT长度使用nextpow2确定FFT长度能显著提升计算速度尤其是当nt不是2的幂时。参考点选择代码中以第一道为参考点x_ref x_pos(1)。你也可以改为排列中心x_ref mean(x_pos)这有时能减少因波前非平面性引起的误差。循环优化最内层循环是对速度v_test的遍历。这里我选择在频率循环内嵌套速度循环结构清晰。对于nr道数很大的情况U_f .* phase_shift_factor这行利用MATLAB的广播机制进行向量化乘法比在道数上再套一层循环快得多。后处理直接取最大值得到的频散曲线可能包含“毛刺”。第8步提供了一个简单的基于局部中值的去噪方法这在处理低信噪比数据时非常有效。3.2 可视化函数plot_dispersion_results.m频散曲线和速度谱的可视化至关重要。function plot_dispersion_results(velocity_spectrum, f_axis, v_axis, disp_curve, fig_title) % 绘制相速度谱和提取的频散曲线 % 输入 % velocity_spectrum, f_axis, v_axis - 来自主函数 % disp_curve - 提取的频散曲线 [频率, 速度] % fig_title - 图标题 figure(Position, [100, 100, 900, 500]); % 子图1相速度谱能量谱 subplot(1, 2, 1); imagesc(v_axis, f_axis, velocity_spectrum); set(gca, YDir, normal); % 确保频率轴从低到高 xlabel(相速度 (m/s)); ylabel(频率 (Hz)); title([fig_title, - 相速度谱]); colorbar; colormap(jet); % 使用jet色图能量高亮显示更明显 axis tight; hold on; % 在速度谱上叠加提取的频散曲线 plot(disp_curve(:,2), disp_curve(:,1), w-, LineWidth, 2.5); plot(disp_curve(:,2), disp_curve(:,1), k--, LineWidth, 1.5); % 黑白双线使其在任何背景下都清晰 % 子图2单独的频散曲线 subplot(1, 2, 2); plot(disp_curve(:,1), disp_curve(:,2), b-o, LineWidth, 2, MarkerSize, 4, MarkerFaceColor, b); xlabel(频率 (Hz)); ylabel(相速度 (m/s)); title([fig_title, - 提取的频散曲线]); grid on; axis tight; % 通常频散曲线随频率升高速度降低可设置Y轴范围 ylim([min(v_axis), max(v_axis)]); end绘图技巧set(gca, YDir, normal)这是关键imagesc默认的Y轴方向是反的原点在左上角这个命令将其纠正使低频在下高频在上符合我们的阅读习惯。叠加曲线在速度谱上用黑白双线叠加频散曲线确保了无论在哪种颜色映射下曲线都清晰可见。子图布局并排显示速度谱和频散曲线方便对比检查提取结果是否合理地位于能量团的主轴上。3.3 数据预处理函数preprocess_sw_data.m原始数据通常不能直接使用预处理能极大提升效果。function data_proc preprocess_sw_data(data_raw, dt, t_start, t_window, taper_ratio, filter_low, filter_high) % 面波数据预处理 % 输入 % data_raw - 原始数据矩阵 % dt - 时间采样率 % t_start - 面波窗起始时间 (秒)相对于记录开始 % t_window - 面波窗长度 (秒) % taper_ratio - 时域两端taper的比例 (0~0.5)用于减少截断效应 % filter_low, filter_high - 带通滤波器的低、高截止频率 (Hz)设为0或[]则不滤波 % 输出 % data_proc - 预处理后的数据 [nt_raw, nr] size(data_raw); t_axis_raw (0:nt_raw-1)*dt; % 1. 截取面波时间窗 start_idx max(1, round(t_start/dt) 1); end_idx min(nt_raw, round((t_start t_window)/dt) 1); data_win data_raw(start_idx:end_idx, :); nt_win size(data_win, 1); % 2. 去除各道直流分量 (减去均值) data_win data_win - mean(data_win, 1); % 3. 时域加窗 (Taper) 以减少频谱泄漏 if taper_ratio 0 taper_len round(nt_win * taper_ratio); taper_win tukeywin(nt_win, 2*taper_ratio); % 使用Tukey窗taper_ratio控制平顶和锥化部分比例 % 如果信号处理工具箱没有tukeywin可以用汉宁窗部分替代 % taper_win hanning(nt_win); % taper_win(1:taper_len) linspace(0,1,taper_len); % taper_win(end-taper_len1:end) linspace(1,0,taper_len); data_win data_win .* taper_win; end % 4. 带通滤波 (保留面波有效频段) if ~isempty(filter_low) filter_low 0 ~isempty(filter_high) filter_high filter_low fs 1/dt; % 设计一个巴特沃斯带通滤波器 [b, a] butter(4, [filter_low, filter_high]/(fs/2), bandpass); % 使用filtfilt进行零相位滤波避免波形畸变 data_win filtfilt(b, a, data_win); end % 5. 可选能量均衡对各道数据乘以一个增益因子补偿几何扩散 % 这里采用简单的偏移距相关增益假设能量随1/sqrt(x)衰减 offset (0:nr-1) * mean(diff(x_pos)); % 需要传入x_pos这里假设已知 gain sqrt(offset / min(offset(offset0))); % 以最近的非零偏移距道为参考 gain(isinf(gain)|isnan(gain)) 1; % 处理第一道可能为0的情况 data_win data_win .* gain; data_proc data_win; end预处理要点时间窗截取只保留包含主要面波能量的时间段能有效压制体波和噪声干扰。t_start需要根据初至时间手动估算或通过其他方式确定。零相位滤波filtfilt函数进行前向-后向滤波避免了普通滤波filter引起的相位失真这对于依赖相位信息的相移法至关重要。能量均衡远道信号弱近道信号强均衡处理可以避免远道信号在叠加中被“淹没”。这里使用sqrt(offset)是一种近似实际中可能需要根据数据情况调整增益函数。4. 实战演练用合成数据测试与验证程序在处理真实数据前用合成数据验证程序是必不可少的一步。它能帮你确认程序逻辑正确并理解参数的影响。4.1 生成合成面波记录我们模拟一个简单的层状介质模型生成理论频散曲线然后合成多道记录。function [syn_data, t, x, true_curve] generate_synthetic_sw(dt, nt, dx, nr, v_model, f_range) % 生成合成面波记录基于频散曲线和简谐波叠加 % 输入 % dt, nt, dx, nr - 时间采样间隔、点数、道间距、道数 % v_model - 层状模型参数矩阵 [厚度(m), Vs(m/s), Vp(m/s), 密度(g/cm3)]最后一行是半空间 % f_range - 要合成的频率范围 [f_min, f_max] 和点数 nf % 输出 % syn_data - 合成数据矩阵 (nt x nr) % t, x - 时间轴和偏移距轴 % true_curve - 用于合成的理论频散曲线 [频率 相速度] % 1. 计算理论频散曲线 (这里调用一个外部函数例如基于Haskell-Thomson矩阵法的程序) % 假设已有函数 calc_dispersion_curve 返回频率和相速度 % [f_theory, v_theory] calc_dispersion_curve(v_model, f_range); % 为演示我们简单假设一个频散曲线速度随频率升高线性降低 f_theory linspace(f_range(1), f_range(2), f_range(3)); v_theory 500 - 100 * (f_theory - f_theory(1)) / (f_theory(end) - f_theory(1)); % 从500m/s降到400m/s true_curve [f_theory(:), v_theory(:)]; % 2. 生成时间和空间轴 t (0:nt-1)*dt; x (0:nr-1)*dx; % 3. 合成记录对每个频率成分生成一个以该频率对应相速度传播的平面波 syn_data zeros(nt, nr); for i 1:length(f_theory) f f_theory(i); v v_theory(i); % 该频率成分的波数 k 2 * pi * f / v; % 生成一个随机的初始相位和振幅模拟实际信号的随机性 A 1.0 / sqrt(f); % 振幅随频率衰减粗略模拟源频谱 phi0 2*pi*rand(); % 随机初始相位 % 为所有时间和空间点生成该频率的波场并叠加 % 这里使用向量化操作提高速度 [T, X] meshgrid(t, x); wave_component A * sin(2*pi*f*T - k*X phi0); syn_data syn_data wave_component; end % 4. 添加高斯白噪声 signal_power mean(syn_data(:).^2); snr_db 20; % 信噪比单位dB noise_power signal_power / (10^(snr_db/10)); noise sqrt(noise_power) * randn(size(syn_data)); syn_data syn_data noise; % 5. 简单滤波去除过高过低频率 fs 1/dt; [b, a] butter(4, [f_range(1)*0.8, f_range(2)*1.2]/(fs/2), bandpass); syn_data filtfilt(b, a, syn_data); end4.2 运行测试与结果分析现在我们将整个流程串起来测试。% 测试脚本 test_phase_shift.m clear; close all; clc; % 1. 生成合成数据参数 dt 0.001; % 1ms采样 nt 1024; % 1024个时间点 dx 2.0; % 2米道间距 nr 48; % 48道 v_model [5, 200, 600, 1.8; % 第一层5m厚Vs200m/s 10, 300, 900, 1.9; % 第二层 inf, 500, 1500, 2.0]; % 半空间 f_range [5, 50, 50]; % 频率从5Hz到50Hz共50个点 [syn_data, t, x, true_curve] generate_synthetic_sw(dt, nt, dx, nr, v_model, f_range); % 2. 预处理数据这里简单处理主要做滤波 data_proc preprocess_sw_data(syn_data, dt, 0.1, 0.5, 0.05, 5, 60); % t_start0.1s, t_window0.5s, taper 5%, 带通5-60Hz % 3. 设置相移法参数并运行 v_min 150; v_max 600; dv 2; % 速度扫描间隔2m/s精度高但计算量稍大 f_min 5; f_max 50; [disp_curve, velocity_spectrum, f_axis, v_axis] ... dispersion_phase_shift(data_proc, dt, dx, v_min, v_max, dv, f_min, f_max); % 4. 可视化结果 plot_dispersion_results(velocity_spectrum, f_axis, v_axis, disp_curve, 合成数据测试); hold on; % 在频散曲线子图上叠加理论曲线 subplot(1,2,2); plot(true_curve(:,1), true_curve(:,2), r--, LineWidth, 2); legend(提取曲线, 理论曲线, Location, best);运行这个脚本你应该能看到速度谱上有一条清晰的能量带提取的频散曲线蓝色实线与理论曲线红色虚线基本吻合。这验证了程序的基本正确性。5. 处理实测数据参数调优与常见问题排查合成数据很理想但实测数据充满挑战。下面结合我处理城市背景噪声或主动源面波数据的经验分享关键步骤和避坑指南。5.1 实测数据准备与初步观察假设你有一个SEG-Y格式的野外数据field_data.sgy。第一步是读入并观察。% 使用开源工具箱如 read_segy 或MATLAB自带函数需Signal Processing Toolbox % 这里假设数据已读入为矩阵 data_raw并获得了 dt 和 dx。 % 绘制原始单炮记录 figure; imagesc(1:nr, t, data_raw); set(gca, YDir, reverse); % 地震数据显示通常时间向下增加 xlabel(道号); ylabel(时间 (s)); title(原始单炮记录); colorbar; colormap(gray);观察记录识别出直达波、折射波、反射波和面波通常是最强、延续时间最长、呈扫帚状散开的能量团。确定面波的主要时间窗口。5.2 关键参数选择策略相移法的效果严重依赖以下几个参数选择不当会导致失败。速度扫描范围[v_min, v_max]和间隔dv策略先宽后窄。第一次处理时v_min可以设得很低如100 m/sv_max设得较高如800-1000 m/sdv可以大一些如10 m/s快速查看能量团的大致位置。依据根据工区地质经验如软土Vs一般150-300 m/s硬土或风化岩300-500 m/s完整岩石500 m/s。观察第一次生成的速度谱能量主要集中在哪个速度区间然后缩小范围并减小dv如2 m/s甚至1 m/s以提高精度。注意dv太小会急剧增加计算量且可能引入速度谱的“锯齿状”噪声。需要在精度和效率间权衡。频率分析范围[f_min, f_max]策略基于数据频谱和勘探深度目标确定。如何确定计算所有道的平均振幅谱。data_fft_all fft(data_proc, NFFT); amp_spectrum mean(abs(data_fft_all(1:NFFT/21, :)), 2); f_axis_full (0:NFFT/2)/(NFFT*dt); figure; plot(f_axis_full, amp_spectrum); xlabel(频率(Hz)); ylabel(平均振幅);从振幅谱上找到面波能量占优的频带。通常主动源面波有效频带在几Hz到几十Hz。f_max不宜超过尼奎斯特频率1/(2*dt)的一半。f_min不宜低于有效信号的最低频率。道间距dx与空间假频核心问题这是最容易被忽略的坑。相移法在波数域操作必须满足空间采样定理即道间距dx必须小于最小波长的一半dx λ_min / 2 v_min / (2 * f_max)。举例如果你的最高分析频率f_max50Hz预计最浅层速度最低的相速度v_min150m/s那么最小波长λ_min 150/50 3m。要求dx 1.5m。如果你的实际dx2m那么在50Hz附近就会出现空间假频速度谱能量会模糊甚至出现虚假的高速度能量团。解决方案如果dx不满足要求要么降低f_max要么在分析前对数据进行空间插值需谨慎会引入误差要么接受高频段结果不可靠的事实。5.3 常见问题与诊断当你得到的速度谱看起来不对劲时可以按以下流程排查问题现象可能原因诊断与解决方案速度谱能量分散没有清晰的能量团1. 信噪比太低。2. 面波窗选取不准包含了太多非面波能量。3. 波前非平面波假设不成立近场效应。1.检查原始记录增强显示增益看面波是否清晰。尝试叠加或滤波提升信噪比。2.调整时间窗尝试不同的t_start和t_window确保只截取最“干净”的面波段。3.检查偏移距如果最小偏移距太大可能已进入波前曲率明显的区域。可尝试以排列中心为参考点重新计算。提取的频散曲线在某个频率发生剧烈跳变1. 该频率处存在较强的干扰波如声波、车噪或反射波。2. 空间假频在该频率出现。3. 速度扫描间隔dv太大错过了真实的能量峰值。1.频谱分析查看该频率成分的单道频谱和所有道的相位关系是否有异常道。2.验证空间采样计算v_min/(2*f_max)与dx对比。如果接近或小于dx则高频跳变很可能是假频所致。3.减小dv在跳变频率附近缩小速度扫描范围并减小dv重新计算。高频段30Hz速度谱能量很弱曲线提取困难1. 高频信号本身衰减快能量弱。2. 检波器耦合或仪器响应在高频段不佳。3. 预处理滤波时不小心滤掉了高频。1.检查预处理确认带通滤波的f_high设置正确没有过早截断。2.能量均衡应用或调整preprocess_sw_data中的增益函数适当提升远道对高频敏感的权重。3.接受现实对于浅层勘探高频信号可能确实很弱频散曲线在高频段不连续是正常的。速度谱出现多条平行的能量带多模式频散。面波尤其是瑞雷波通常存在基阶和高阶模式。这是正常现象相移法将不同模式的能量都成像出来了。你需要判断哪一条是基阶模式通常速度最低、能量最强的那一条。在反演时可能需要分别提取基阶和高阶模式曲线。5.4 一个完整的实测数据处理示例% 假设 field_data, dt, dx 已加载 % 步骤1初步观察与预处理 figure(1); % ... 绘制原始记录确定面波窗大致为 0.2s 到 0.8s ... t_start 0.2; t_window 0.6; data_proc preprocess_sw_data(field_data, dt, t_start, t_window, 0.05, 5, 80); % 步骤2参数设置基于初步分析和工区经验 % 检查空间假频假设期望 v_min180m/s, f_max80Hz, 则需 dx 180/(2*80)1.125m。 % 实际 dx2m因此需要将 f_max 限制在 180/(2*2)45Hz 以下以保证无假频。 f_max_safe 180 / (2 * 2); % 约45Hz f_max_used min(80, f_max_safe); % 取保守值 v_min 150; v_max 600; dv 3; % 第一次用稍大的间隔 f_min 5; f_max f_max_used; % 使用安全的最大频率 % 步骤3运行相移法 [disp_curve1, vel_spec1, f_axis1, v_axis1] ... dispersion_phase_shift(data_proc, dt, dx, v_min, v_max, dv, f_min, f_max); plot_dispersion_results(vel_spec1, f_axis1, v_axis1, disp_curve1, 实测数据-初版); % 步骤4根据初版结果优化参数 % 假设从 vel_spec1 看到能量主要集中在 180-400 m/s v_min2 170; v_max2 450; dv2 1.5; % 缩小范围提高精度 % 同时发现10Hz以下和40Hz以上能量很弱可以调整频率范围 f_min2 8; f_max2 40; [disp_curve2, vel_spec2, f_axis2, v_axis2] ... dispersion_phase_shift(data_proc, dt, dx, v_min2, v_max2, dv2, f_min2, f_max2); plot_dispersion_results(vel_spec2, f_axis2, v_axis2, disp_curve2, 实测数据-优化后); % 步骤5结果后处理与输出 % 对提取的曲线进行平滑处理例如移动平均 windowSize 7; disp_curve_smooth disp_curve2; disp_curve_smooth(:,2) movmedian(disp_curve2(:,2), windowSize); % 使用中值滤波抗野值 % 绘制最终对比图 figure; plot(disp_curve2(:,1), disp_curve2(:,2), b., MarkerSize, 10); hold on; plot(disp_curve_smooth(:,1), disp_curve_smooth(:,2), r-, LineWidth, 2); xlabel(频率 (Hz)); ylabel(相速度 (m/s)); legend(原始提取点, 平滑后曲线, Location, best); grid on; title(最终频散曲线);通过这个迭代优化的过程你就能从复杂的实测数据中提取出相对稳定、可靠的频散曲线为后续的反演解释打下坚实基础。记住没有一套参数能通吃所有数据耐心调整和基于物理意义的判断才是用好这个工具的关键。本文还有配套的精品资源点击获取