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

读取SEG2格式地震勘探数据的MATLAB实现

news 2025/10/1 1:33:34

读取SEG2格式的地震勘探数据，并对数据进行可视化和初步处理。

function [seismic_data, headers] = read_seg2(filename)% 读取SEG2格式地震勘探数据% 输入：%   filename - SEG2文件路径% 输出：%   seismic_data - 地震数据矩阵（道 x 时间）%   headers - 包含文件头和道头的元胞数组% 检查文件是否存在if ~exist(filename, 'file')error('文件不存在: %s', filename);end% 打开文件fid = fopen(filename, 'r', 'ieee-be'); % SEG2通常使用大端字节序if fid == -1error('无法打开文件: %s', filename);end% 初始化输出seismic_data = [];headers = struct('file_header', [], 'trace_headers', []);try%% 1. 读取文件头file_header = read_file_header(fid);headers.file_header = file_header;fprintf('读取SEG2文件: %s\n', filename);fprintf('文件标识符: %s\n', file_header.file_id);fprintf('格式版本: %s\n', file_header.format_version);fprintf('扫描类型: %s\n', file_header.scan_type);fprintf('道数: %d\n', file_header.num_traces);fprintf('采样点数/道: %d\n', file_header.samples_per_trace);fprintf('采样间隔: %.6f ms\n', file_header.sample_interval*1000);%% 2. 读取道头和数据trace_headers = cell(file_header.num_traces, 1);seismic_data = zeros(file_header.num_traces, file_header.samples_per_trace);for trace_idx = 1:file_header.num_traces% 读取道头trace_header = read_trace_header(fid, trace_idx);trace_headers{trace_idx} = trace_header;% 读取地震数据seismic_data(trace_idx, :) = read_trace_data(fid, trace_header);% 显示进度if mod(trace_idx, 10) == 0 || trace_idx == file_header.num_tracesfprintf('已读取 %d/%d 道数据...\n', trace_idx, file_header.num_traces);endendheaders.trace_headers = trace_headers;%% 3. 关闭文件fclose(fid);fprintf('成功读取 %d 道地震数据\n', file_header.num_traces);catch ME% 发生错误时关闭文件fclose(fid);rethrow(ME);end
endfunction file_header = read_file_header(fid)% 读取SEG2文件头% SEG2文件头起始标识file_id = fread(fid, 4, 'char')';if ~isequal(char(file_id), 'SEG2')error('无效的SEG2文件: 缺少SEG2标识');end% 创建文件头结构file_header = struct();file_header.file_id = 'SEG2';% 读取文件描述块while true% 读取块类型和大小block_type = fread(fid, 1, 'uint8');block_size = fread(fid, 1, 'uint16');if isempty(block_type) || isempty(block_size)error('文件头读取不完整');end% 文件描述块 (类型1)if block_type == 1% 读取描述符块descriptor_block = fread(fid, block_size, 'char')';desc_str = char(descriptor_block);% 解析描述符字符串tokens = regexp(desc_str, '(\w+)\s*=\s*([^;]+);', 'tokens');for i = 1:length(tokens)key = tokens{i}{1};value = tokens{i}{2};% 转换数值类型num_value = str2double(value);if ~isnan(num_value)value = num_value;endfile_header.(key) = value;end% 确保必要字段存在required_fields = {'num_traces', 'samples_per_trace', 'sample_interval'};for i = 1:length(required_fields)if ~isfield(file_header, required_fields{i})error('缺少必要字段: %s', required_fields{i});endendbreak; % 文件描述块后是道数据块else% 跳过未知块fseek(fid, block_size, 'cof');endend
endfunction trace_header = read_trace_header(fid, trace_idx)% 读取道头信息% 道头起始标识 (类型48)block_type = fread(fid, 1, 'uint8');if block_type ~= 48error('道 %d: 无效的道头起始标识', trace_idx);end% 道头大小block_size = fread(fid, 1, 'uint16');% 读取道描述块descriptor_block = fread(fid, block_size, 'char')';desc_str = char(descriptor_block);% 初始化道头结构trace_header = struct();trace_header.trace_index = trace_idx;trace_header.block_size = block_size;% 解析描述符字符串tokens = regexp(desc_str, '(\w+)\s*=\s*([^;]+);', 'tokens');for i = 1:length(tokens)key = tokens{i}{1};value = tokens{i}{2};% 转换数值类型num_value = str2double(value);if ~isnan(num_value)value = num_value;endtrace_header.(key) = value;end% 确保必要字段存在if ~isfield(trace_header, 'data_type')trace_header.data_type = 1; % 默认为IBM浮点数end% 数据大小和采样点数if ~isfield(trace_header, 'num_samples')error('道 %d: 缺少采样点数信息', trace_idx);end% 数据大小 (字节)switch trace_header.data_typecase 1 % IBM浮点数 (4字节)bytes_per_sample = 4;case 2 % 32位整数 (4字节)bytes_per_sample = 4;case 3 % 16位整数 (2字节)bytes_per_sample = 2;case 4 % IEEE浮点数 (4字节)bytes_per_sample = 4;case 5 % 未使用error('道 %d: 不支持的数据类型: %d', trace_idx, trace_header.data_type);case 6 % 未使用error('道 %d: 不支持的数据类型: %d', trace_idx, trace_header.data_type);case 7 % 24位整数 (3字节)bytes_per_sample = 3;otherwiseerror('道 %d: 未知的数据类型: %d', trace_idx, trace_header.data_type);endtrace_header.bytes_per_sample = bytes_per_sample;trace_header.data_size = trace_header.num_samples * bytes_per_sample;
endfunction trace_data = read_trace_data(fid, trace_header)% 读取地震道数据% 数据块起始标识 (类型64)block_type = fread(fid, 1, 'uint8');if block_type ~= 64error('道 %d: 无效的数据块起始标识', trace_header.trace_index);end% 数据块大小block_size = fread(fid, 1, 'uint16');expected_size = trace_header.data_size;if block_size ~= expected_sizewarning('道 %d: 数据块大小(%d)与预期(%d)不符', ...trace_header.trace_index, block_size, expected_size);end% 读取二进制数据switch trace_header.data_typecase 1 % IBM浮点数% 转换为IEEE浮点数 (需要转换函数)raw_data = fread(fid, trace_header.num_samples, 'uint32');trace_data = ibm2ieee(raw_data);case 2 % 32位整数trace_data = fread(fid, trace_header.num_samples, 'int32');case 3 % 16位整数trace_data = fread(fid, trace_header.num_samples, 'int16');case 4 % IEEE浮点数trace_data = fread(fid, trace_header.num_samples, 'float32');case 7 % 24位整数% 读取24位数据 (3字节)raw_data = fread(fid, trace_header.num_samples*3, 'uint8');% 重组为32位整数trace_data = reshape(raw_data, 3, []);trace_data = double(typecast(trace_data(:), 'int32')) / (2^8);otherwiseerror('道 %d: 不支持的数据类型: %d', trace_header.trace_index, trace_header.data_type);end% 确保数据方向正确if size(trace_data, 1) == 1trace_data = trace_data';end
endfunction ieee_data = ibm2ieee(ibm_data)% 将IBM浮点数转换为IEEE浮点数% IBM浮点数格式: 1位符号, 7位指数, 24位尾数ieee_data = zeros(size(ibm_data));for i = 1:length(ibm_data)% 提取符号位sign_bit = bitand(ibm_data(i), hex2dec('80000000'));% 提取指数exponent = double(bitshift(bitand(ibm_data(i), hex2dec('7F000000')), -24));% 提取尾数fraction = double(bitand(ibm_data(i), hex2dec('00FFFFFF')));% 计算值 (IBM格式: (-1)^s * 16^(e-64) * 0.f)value = fraction / 2^24;if exponent > 0value = value * 16^(exponent - 64);elsevalue = value / 16^(64 - exponent);endif sign_bitvalue = -value;endieee_data(i) = value;end
end%% 主程序：读取和可视化SEG2数据
clc;
clear;
close all;% 设置文件路径
filename = 'seismic_data.seg2'; % 替换为你的SEG2文件路径% 读取SEG2文件
[seismic_data, headers] = read_seg2(filename);% 获取基本参数
fs = 1 / (headers.file_header.sample_interval * 1000); % 采样率 (Hz)
num_traces = size(seismic_data, 1);
num_samples = size(seismic_data, 2);
time_axis = (0:num_samples-1) / fs; % 时间轴 (秒)%% 可视化地震数据
% 1. 单道波形显示
figure('Position', [100, 100, 1200, 800], 'Name', '地震数据可视化');% 选择三道进行显示
trace_indices = [1, floor(num_traces/2), num_traces];subplot(2, 2, 1);
hold on;
colors = lines(3);
for i = 1:3idx = trace_indices(i);plot(time_axis, seismic_data(idx, :) + i*0.5, 'Color', colors(i, :), 'LineWidth', 1.5);
end
hold off;
xlabel('时间 (s)');
ylabel('振幅 (偏移)');
title('三道地震波形');
legend(sprintf('道 %d', trace_indices(1)), sprintf('道 %d', trace_indices(2)), sprintf('道 %d', trace_indices(3)));
grid on;% 2. 地震剖面显示
subplot(2, 2, 2);
imagesc(1:num_traces, time_axis, seismic_data');
colormap(seismic_colormap); % 使用地震数据专用色标
colorbar;
xlabel('道号');
ylabel('时间 (s)');
title('地震剖面');
axis tight;% 3. 振幅统计
subplot(2, 2, 3);
trace_amplitudes = rms(seismic_data, 2);
plot(1:num_traces, trace_amplitudes, 'b-o', 'LineWidth', 1.5);
xlabel('道号');
ylabel('RMS 振幅');
title('道振幅统计');
grid on;% 4. 频谱分析
subplot(2, 2, 4);
trace_idx = floor(num_traces/2);
[pxx, f] = pwelch(seismic_data(trace_idx, :), [], [], [], fs);
semilogy(f, pxx, 'r', 'LineWidth', 1.5);
xlabel('频率 (Hz)');
ylabel('功率谱密度');
title(sprintf('道 %d 频谱分析', trace_idx));
grid on;
xlim([0, fs/2]);%% 显示道头信息
% 选择一道显示详细信息
trace_idx = 1;
trace_header = headers.trace_headers{trace_idx};fprintf('\n===== 道头信息 (道 %d) =====\n', trace_idx);
disp(trace_header);%% 数据预处理
% 1. 增益恢复 (可选)
% 2. 带通滤波
filtered_data = apply_bandpass_filter(seismic_data, fs, 5, 10, 80, 100);% 3. 振幅均衡
balanced_data = apply_agc(filtered_data, 0.1*fs); % 100ms窗口%% 显示处理后的数据
figure('Position', [100, 100, 1000, 600], 'Name', '处理后的地震数据');% 原始数据
subplot(2, 1, 1);
imagesc(1:num_traces, time_axis, seismic_data');
colormap(seismic_colormap);
caxis(prctile(seismic_data(:), [5, 95]));
colorbar;
title('原始地震剖面');
xlabel('道号');
ylabel('时间 (s)');% 处理后的数据
subplot(2, 1, 2);
imagesc(1:num_traces, time_axis, balanced_data');
colormap(seismic_colormap);
caxis(prctile(balanced_data(:), [5, 95]));
colorbar;
title('处理后地震剖面 (带通滤波 + AGC)');
xlabel('道号');
ylabel('时间 (s)');%% 保存处理后的数据
output_filename = 'processed_seismic_data.mat';
save(output_filename, 'balanced_data', 'headers', 'fs');
fprintf('\n处理后的数据已保存至: %s\n', output_filename);%% 辅助函数
function cmap = seismic_colormap()% 创建地震数据专用色标n = 64;r = [zeros(1, n/2), linspace(0, 1, n/2)];g = [linspace(0, 1, n/2), linspace(1, 0, n/2)];b = [linspace(0.5, 1, n/2), linspace(1, 0, n/2)];cmap = [r(:), g(:), b(:)];
endfunction filtered_data = apply_bandpass_filter(data, fs, f1, f2, f3, f4)% 应用带通滤波器 (Butterworth)% f1, f2: 通带下限频率 (Hz)% f3, f4: 通带上限频率 (Hz)nyquist = fs / 2;Wn = [f2, f3] / nyquist; % 通带频率归一化[b, a] = butter(4, Wn, 'bandpass');% 对每道数据进行滤波filtered_data = zeros(size(data));for i = 1:size(data, 1)filtered_data(i, :) = filtfilt(b, a, double(data(i, :)));end
endfunction balanced_data = apply_agc(data, window_size)% 应用自动增益控制 (AGC)% window_size: 窗口大小 (采样点数)balanced_data = zeros(size(data));num_samples = size(data, 2);for i = 1:size(data, 1)trace = data(i, :);agc_trace = zeros(1, num_samples);for j = 1:num_samples% 计算窗口范围start_idx = max(1, j - floor(window_size/2));end_idx = min(num_samples, j + floor(window_size/2));% 计算窗口内的RMS值window = trace(start_idx:end_idx);rms_val = sqrt(mean(window.^2));% 避免除以零if rms_val > 0agc_trace(j) = trace(j) / rms_val;elseagc_trace(j) = trace(j);endendbalanced_data(i, :) = agc_trace;end
end