当前位置：首页 > news >正文

地震数据处理基本程序（MATLAB实现）

news 2026/5/11 20:20:15

一、数据读取与预处理

%% 1. 读取SEGY格式地震数据
% 使用内置函数读取SEGY文件（示例文件：input.sgy）
segyStruct = segyread('input.sgy');% 提取数据体与道头信息
data = segyStruct.data;      % 数据体 [nTraces × nSamples]
nTraces = size(data, 1);     % 道数
nSamples = size(data, 2);    % 采样点数
dt = segyStruct.b;           % 采样间隔（秒）%% 2. 数据预处理
% 去直流分量（均值偏移）
data = detrend(data, 'constant');% 时间轴构建
time = (0:nSamples-1) * dt;  % 时间向量 [0, dt, 2*dt, ..., (nSamples-1)*dt]

二、去噪处理

1. 小波阈值去噪

%% 小波去噪（db4小波基，5层分解）
[c, l] = wavedec(data(:, 1), 5, 'db4');  % 对单道信号分解
sigma = median(abs(c)) / 0.6745;         % 估计噪声标准差
thr = wthrmngr('dw1ddenoLVL', 'penalhi', c, sigma);  % 阈值计算
denoised = wdenoise(data(:, 1), 5, 'Wavelet', 'db4', 'ThresholdRule', 'soft');  % 去噪

2. 变分模态分解（VMD）去噪

%% VMD去噪（参数设置）
alpha = 2000;   % 带宽约束
tau = 0.1;      % 噪声容忍
K = 4;          % 分解模态数
tol = 1e-6;     % 收敛容差% 执行VMD分解
[u, ~, ~] = VMD2D(data, alpha, tau, K, 0, 1, tol);% 保留前K-1个IMF分量（去噪）
denoised_vmd = sum(u(1:K-1, :), 1);

三、滤波处理

1. 带通滤波（FIR滤波器）

%% 带通滤波（5-80 Hz，4阶巴特沃斯）
fs = 1/dt;      % 采样率（Hz）
fcuts = [5, 80]; % 截止频率
order = 4;      % 滤波器阶数% 设计滤波器
[b, a] = butter(order, fcuts/(fs/2), 'bandpass');% 应用滤波器
filtered_data = filtfilt(b, a, data');  % 对每道信号滤波

2. F-K滤波（频散分析）

%% F-K滤波（压制线性干扰）
dx = 25;        % 空间采样间隔（米）
fk = fftfilt(data, dx);  % 频散分析
f_k_filtered = fk .* (abs(fft(fk)) > 0.1*max(abs(fft(fk))));  % 保留主频成分

四、叠加处理

1. 动校正与CMP叠加

%% 动校正（假设叠加速度为2000 m/s）
vel = 2000;     % 叠加速度（m/s）
time_shift = (0:nSamples-1) * dx / vel;  % 动校正时移量% 对每道应用时移
corrected_data = zeros(size(data));
for i = 1:nTracescorrected_data(i, :) = interp1(time, data(i, :), time - time_shift(i));
end%% CMP叠加
stacked_data = sum(corrected_data, 1);  % 简单叠加

2. 加权叠加（反距离加权）

%% 加权叠加（道间距dx=25m）
weights = 1 ./ (1 + (0:nTraces-1)*dx);  % 反距离权重
weighted_stack = sum(corrected_data .* weights', 1);

五、可视化与结果分析

%% 绘制处理流程对比图
figure;% 原始信号
subplot(3, 2, 1);
plot(time, data(1, :));
title('原始信号');
xlabel('时间 (s)');
ylabel('振幅');% 小波去噪
subplot(3, 2, 2);
plot(time, denoised);
title('小波去噪');
xlabel('时间 (s)');
ylabel('振幅');% 带通滤波
subplot(3, 2, 3);
plot(time, filtered_data(1, :));
title('带通滤波');
xlabel('时间 (s)');
ylabel('振幅');% CMP叠加
subplot(3, 2, 4);
plot(time, stacked_data);
title('CMP叠加');
xlabel('时间 (s)');
ylabel('振幅');% 频谱对比
subplot(3, 2, 5);
pwelch(data(1,:), [], [], [], fs);
hold on;
pwelch(filtered_data(1,:), [], [], [], fs);
legend('原始信号', '滤波后信号');
title('频谱对比');

参考代码地震数据处理基本程序，保存去噪，滤波，叠加，适合学习理解地震数据处理的基本原理 www.youwenfan.com/contentcnr/99629.html