读取SEG2格式地震勘探数据的MATLAB实现
读取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 == -1
error('无法打开文件: %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_traces
fprintf('已读取 %d/%d 道数据...\n', trace_idx, file_header.num_traces);
end
end
headers.trace_headers = trace_headers;
%% 3. 关闭文件
fclose(fid);
fprintf('成功读取 %d 道地震数据\n', file_header.num_traces);
catch ME
% 发生错误时关闭文件
fclose(fid);
rethrow(ME);
end
end
function 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;
end
file_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});
end
end
break; % 文件描述块后是道数据块
else
% 跳过未知块
fseek(fid, block_size, 'cof');
end
end
end
function trace_header = read_trace_header(fid, trace_idx)
% 读取道头信息
% 道头起始标识 (类型48)
block_type = fread(fid, 1, 'uint8');
if block_type ~= 48
error('道 %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;
end
trace_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_type
case 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;
otherwise
error('道 %d: 未知的数据类型: %d', trace_idx, trace_header.data_type);
end
trace_header.bytes_per_sample = bytes_per_sample;
trace_header.data_size = trace_header.num_samples * bytes_per_sample;
end
function trace_data = read_trace_data(fid, trace_header)
% 读取地震道数据
% 数据块起始标识 (类型64)
block_type = fread(fid, 1, 'uint8');
if block_type ~= 64
error('道 %d: 无效的数据块起始标识', trace_header.trace_index);
end
% 数据块大小
block_size = fread(fid, 1, 'uint16');
expected_size = trace_header.data_size;
if block_size ~= expected_size
warning('道 %d: 数据块大小(%d)与预期(%d)不符', ...
trace_header.trace_index, block_size, expected_size);
end
% 读取二进制数据
switch trace_header.data_type
case 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);
otherwise
error('道 %d: 不支持的数据类型: %d', trace_header.trace_index, trace_header.data_type);
end
% 确保数据方向正确
if size(trace_data, 1) == 1
trace_data = trace_data';
end
end
function 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 > 0
value = value * 16^(exponent - 64);
else
value = value / 16^(64 - exponent);
end
if sign_bit
value = -value;
end
ieee_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:3
idx = 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(:)];
end
function 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
end
function 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 > 0
agc_trace(j) = trace(j) / rms_val;
else
agc_trace(j) = trace(j);
end
end
balanced_data(i, :) = agc_trace;
end
end
SEG2文件格式解析
SEG2格式是地震勘探中常用的数据存储格式，包含以下主要部分：
1. 文件头结构
- 文件标识符：4字节的"SEG2"
- 描述符块：包含全局信息
- 道数 (num_traces)
- 每道采样点数 (samples_per_trace)
- 采样间隔 (sample_interval)
- 扫描类型 (scan_type)
2. 道数据块
每个地震道包含：
- 道头：
- 块类型标识 (48)
- 道头大小
- 描述符字符串（包含道号、位置、传感器类型等信息）
- 数据块：
- 块类型标识 (64)
- 数据大小
- 地震数据（多种数据格式）
3. 支持的数据类型
- IBM浮点数（32位）
- 32位整数
- 16位整数（最常见）
- IEEE浮点数
- 24位整数
主要函数
read_seg2：主函数，协调文件读取过程read_file_header：解析SEG2文件头read_trace_header：读取单道头信息read_trace_data：读取单道地震数据ibm2ieee：IBM浮点数转IEEE浮点数
参考代码 读取 地震勘探数据seg2 www.youwenfan.com/contentcne/66097.html
数据处理流程
- 文件验证：检查文件存在性和SEG2标识
- 文件头解析：提取全局参数（道数、采样率等）
- 道数据读取：
- 解析道头信息
- 根据数据类型读取二进制数据
- 特殊格式转换（如IBM浮点数）
- 数据组织：将数据组织为矩阵（道 × 时间）
数据可视化
- 单道波形显示：显示三道典型波形
- 地震剖面：二维剖面显示（使用地震专用色标）
- 振幅统计：显示各道RMS振幅变化
- 频谱分析：单道频谱特性
数据预处理
- 带通滤波：
- 使用Butterworth滤波器
- 可调通带频率（默认5-10Hz至80-100Hz）
- 自动增益控制(AGC)：
- 时间域振幅均衡
- 可调窗口大小（默认100ms）
此代码提供了从SEG2文件读取到预处理的全流程解决方案，适用于陆地地震勘探、海洋地震勘探等多种地震数据采集场景。
浙公网安备 33010602011771号