MATLAB实现相移法提取面波频散曲线：从原理到实战避坑指南

简介：本资源是一套面向地球物理勘探专业师生及科研人员的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 多道叠加与“速度谱”的生成

相移法的巧妙之处在于，它利用了整个排列的所有道。我们不是两两对比，而是将所有道都向一个虚拟的“零偏移距”参考点进行相位移动。

1. 选择参考点：通常选择第一个检波器（最小偏移距）或排列中心作为参考点x_ref。
2. 遍历测试速度：对于一个给定的频率f，我们预设一个相速度v_test的扫描范围（比如从100 m/s到1000 m/s）。
3. 相位移动与叠加：对于排列中的第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速度传播的。
4. 相干叠加：如果实际的相速度恰好等于v_test，那么所有道经过上述相位移动后，它们在参考点处的“估计信号”的相位将会完全对齐。将所有道的这些“搬回来”的频谱在频率f处求和，其幅值将会达到最大。如果v_test不等于真实速度，各道相位参差不齐，叠加后会相互抵消，幅值较小。
5. 构建速度谱：对每个频率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/2+1, :); % 取单边谱 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_len+1: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(offset>0))); % 以最近的非零偏移距道为参考 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); end

4.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厚，Vs=200m/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_start=0.1s, t_window=0.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 关键参数选择策略

相移法的效果严重依赖以下几个参数，选择不当会导致失败。

1. 速度扫描范围[v_min, v_max]和间隔dv：

   • 策略：先宽后窄。第一次处理时，v_min可以设得很低（如100 m/s），v_max设得较高（如800-1000 m/s），dv可以大一些（如10 m/s），快速查看能量团的大致位置。
   • 依据：根据工区地质经验（如软土Vs一般150-300 m/s，硬土或风化岩300-500 m/s，完整岩石>500 m/s）。观察第一次生成的速度谱，能量主要集中在哪个速度区间，然后缩小范围，并减小dv（如2 m/s甚至1 m/s）以提高精度。
   • 注意：dv太小会急剧增加计算量，且可能引入速度谱的“锯齿状”噪声。需要在精度和效率间权衡。

2. 频率分析范围[f_min, f_max]：

   • 策略：基于数据频谱和勘探深度目标确定。
   • 如何确定：计算所有道的平均振幅谱。

     data_fft_all = fft(data_proc, NFFT); amp_spectrum = mean(abs(data_fft_all(1:NFFT/2+1, :)), 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不宜低于有效信号的最低频率。

3. 道间距dx与空间假频：

   • 核心问题：这是最容易被忽略的坑。相移法在波数域操作，必须满足空间采样定理，即道间距dx必须小于最小波长的一半：dx < λ_min / 2 = v_min / (2 * f_max)。
   • 举例：如果你的最高分析频率f_max=50Hz，预计最浅层（速度最低）的相速度v_min=150m/s，那么最小波长λ_min = 150/50 = 3m。要求dx < 1.5m。如果你的实际dx=2m，那么在50Hz附近就会出现空间假频，速度谱能量会模糊甚至出现虚假的高速度能量团。
   • 解决方案：如果dx不满足要求，要么降低f_max，要么在分析前对数据进行空间插值（需谨慎，会引入误差），要么接受高频段结果不可靠的事实。

5.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_min=180m/s, f_max=80Hz, 则需 dx < 180/(2*80)=1.125m。 % 实际 dx=2m，因此需要将 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('最终频散曲线');

通过这个迭代优化的过程，你就能从复杂的实测数据中提取出相对稳定、可靠的频散曲线，为后续的反演解释打下坚实基础。记住，没有一套参数能通吃所有数据，耐心调整和基于物理意义的判断，才是用好这个工具的关键。

本文还有配套的精品资源，点击获取