SGY地震数据自相关与互相关分析
MATLAB实现方案
%% SGY地震数据自相关与互相关分析
clc; clear; close all;
%% 1. 读取SGY地震数据
% 注意：实际使用时需替换为真实SGY文件路径
% 这里使用模拟数据代替
[sgyData, t, traceHeaders] = readSGY('sample.sgy');
% 显示基本信息
[nTraces, nSamples] = size(sgyData);
fprintf('地震数据信息:\n');
fprintf('道数: %d\n', nTraces);
fprintf('采样点数: %d\n', nSamples);
fprintf('采样间隔: %.2f ms\n', t(2)-t(1));
fprintf('总时间: %.2f s\n', t(end)/1000);
%% 2. 数据预处理
% 去除直流分量
sgyData = detrend(sgyData')';
% 带通滤波 (可选)
fs = 1/(t(2)-t(1)); % 采样频率 (Hz)
fcuts = [5, 10, 80, 100]; % 通带范围 (Hz)
order = 4; % 滤波器阶数
[b,a] = butter(order, fcuts/(fs/2), 'bandpass');
filteredData = filtfilt(b, a, sgyData')';
% 选择分析区域
startSample = 500; % 起始采样点
endSample = 1500; % 结束采样点
analysisData = filteredData(:, startSample:endSample);
analysisTime = t(startSample:endSample);
%% 3. 自相关分析
% 选择参考道 (中间道)
refTrace = ceil(nTraces/2);
selectedTrace = analysisData(refTrace, :);
% 计算自相关
[acor, lags] = xcorr(selectedTrace, 'coeff');
lags_ms = lags * (t(2)-t(1)); % 转换为毫秒
% 归一化自相关
acor = acor / max(acor);
%% 4. 互相关分析
% 选择另一道进行对比
compTrace = refTrace + 5;
if compTrace > nTraces
compTrace = refTrace - 5;
end
compareTrace = analysisData(compTrace, :);
% 计算互相关
[xcor, lagsX] = xcorr(selectedTrace, compareTrace, 'coeff');
lagsX_ms = lagsX * (t(2)-t(1));
%% 5. 生成图件
% 设置图形参数
set(0, 'DefaultAxesFontSize', 10);
set(0, 'DefaultTextFontSize', 10);
set(0, 'DefaultLineLineWidth', 1.5);
% 图件1: 原始地震道
figure('Name', '地震道波形', 'Position', [100, 100, 1200, 800]);
subplot(2,1,1);
plot(analysisTime, selectedTrace, 'b', 'LineWidth', 1.5);
hold on;
plot(analysisTime, compareTrace, 'r', 'LineWidth', 1.5);
xlabel('时间 (ms)');
ylabel('振幅');
title('参考道与对比道波形');
legend(['道 ', num2str(refTrace)], ['道 ', num2str(compTrace)]);
grid on;
% 添加地质标记
hold on;
yLimits = ylim;
plot([300, 300], yLimits, 'k--', 'LineWidth', 1);
text(310, yLimits(2)*0.9, '浅层反射', 'FontSize', 10);
plot([700, 700], yLimits, 'k--', 'LineWidth', 1);
text(710, yLimits(2)*0.9, '中层反射', 'FontSize', 10);
plot([1200, 1200], yLimits, 'k--', 'LineWidth', 1);
text(1210, yLimits(2)*0.9, '深层反射', 'FontSize', 10);
% 图件2: 自相关分析
subplot(2,1,2);
plot(lags_ms, acor, 'b', 'LineWidth', 1.5);
xlabel('延迟时间 (ms)');
ylabel('归一化自相关');
title('参考道自相关分析');
grid on;
xlim([-200, 200]);
% 标记主要周期
[peaks, locs] = findpeaks(acor, 'MinPeakHeight', 0.3);
hold on;
for i = 1:min(3, length(locs))
plot(lags_ms(locs(i)), peaks(i), 'ro', 'MarkerSize', 8);
text(lags_ms(locs(i)), peaks(i)+0.05, ...
sprintf('%.1f ms', lags_ms(locs(i))), ...
'HorizontalAlignment', 'center');
end
% 图件3: 互相关分析
figure('Name', '互相关分析', 'Position', [100, 100, 1200, 600]);
subplot(1,2,1);
plot(analysisTime, selectedTrace, 'b', 'LineWidth', 1.5);
hold on;
plot(analysisTime, compareTrace, 'r', 'LineWidth', 1.5);
xlabel('时间 (ms)');
ylabel('振幅');
title('参考道与对比道波形');
legend(['道 ', num2str(refTrace)], ['道 ', num2str(compTrace)]);
grid on;
subplot(1,2,2);
plot(lagsX_ms, xcor, 'm', 'LineWidth', 1.5);
xlabel('延迟时间 (ms)');
ylabel('归一化互相关系数');
title('参考道与对比道互相关分析');
grid on;
xlim([-200, 200]);
% 标记最大相关位置
[maxCorr, maxIdx] = max(xcor);
hold on;
plot(lagsX_ms(maxIdx), maxCorr, 'ko', 'MarkerSize', 10, 'MarkerFaceColor', 'y');
text(lagsX_ms(maxIdx), maxCorr+0.05, ...
sprintf('最大相关: %.2f\n延迟: %.1f ms', maxCorr, lagsX_ms(maxIdx)), ...
'HorizontalAlignment', 'center');
% 图件4: 多道自相关分析
figure('Name', '多道自相关分析', 'Position', [100, 100, 1200, 800]);
numTracesToShow = min(9, nTraces);
cols = 3;
rows = ceil(numTracesToShow/cols);
for i = 1:numTracesToShow
traceIdx = floor((i-1)*nTraces/numTracesToShow) + 1;
traceData = analysisData(traceIdx, :);
[acorTrace, lagsTrace] = xcorr(traceData, 'coeff');
lagsTrace_ms = lagsTrace * (t(2)-t(1));
subplot(rows, cols, i);
plot(lagsTrace_ms, acorTrace, 'b', 'LineWidth', 1.2);
title(['道 ', num2str(traceIdx)]);
xlabel('延迟 (ms)');
ylabel('自相关');
grid on;
xlim([-150, 150]);
end
% 图件5: 相干分析 (多道互相关)
figure('Name', '相干分析', 'Position', [100, 100, 1200, 800]);
% 计算平均相干系数
coherence = zeros(1, length(analysisTime));
windowSize = 50; % 相干计算窗口大小
for i = 1:nTraces-windowSize
windowData = analysisData(i:i+windowSize-1, :);
meanTrace = mean(windowData, 1);
centeredData = windowData - repmat(meanTrace, windowSize, 1);
covMatrix = centeredData * centeredData' / (nSamples-1);
coherence = coherence + diag(covMatrix, 0);
end
coherence = coherence / (nTraces-windowSize);
% 绘制相干图
subplot(2,1,1);
imagesc(analysisTime, 1:nTraces, analysisData);
colormap(gray);
colorbar;
xlabel('时间 (ms)');
ylabel('道号');
title('地震数据剖面');
axis tight;
subplot(2,1,2);
plot(analysisTime, coherence, 'b', 'LineWidth', 2);
xlabel('时间 (ms)');
ylabel('相干系数');
title('地震数据相干性分析');
grid on;
ylim([0, 1]);
%% 6. 保存分析结果
save('seismic_correlation_results.mat', 'acor', 'xcor', 'lags_ms', 'lagsX_ms', 'coherence');
%% 辅助函数: 读取SGY文件 (简化版)
function [data, t, traceHeaders] = readSGY(filename)
% 实际实现需要完整的SGY解析器
% 这里使用模拟数据代替
% 模拟参数
nTraces = 100; % 道数
nSamples = 2000; % 每道采样点数
dt = 0.002; % 采样间隔 (秒)
% 创建时间向量 (毫秒)
t = (0:nSamples-1)*dt*1000;
% 创建模拟地震数据
data = zeros(nTraces, nSamples);
for i = 1:nTraces
% 基础波形
baseWavelet = ricker(30, dt, 1.5);
waveLength = length(baseWavelet);
% 随机反射系数序列
reflectors = randn(1, 15)*0.1;
reflectorTimes = sort(randperm(nSamples-500, 15) + 250);
% 构建地震道
trace = zeros(1, nSamples);
for j = 1:length(reflectors)
pos = reflectorTimes(j);
if pos + waveLength <= nSamples
trace(pos:pos+waveLength-1) = trace(pos:pos+waveLength-1) + reflectors(j)*baseWavelet;
end
end
% 添加噪声
noise = randn(1, nSamples)*0.05;
data(i, :) = trace + noise;
% 添加道间差异
data(i, :) = data(i, :) * (0.9 + 0.2*rand());
end
% 模拟道头信息
traceHeaders = struct();
traceHeaders.TRACE_SAMPLE_COUNT = repmat(nSamples, nTraces, 1);
traceHeaders.SAMPLE_INTERVAL = repmat(dt*1000000, nTraces, 1); % 微秒
traceHeaders.TraceNumber = (1:nTraces)';
end
%% Ricker小波生成函数
function wavelet = ricker(freq, dt, duration)
t = -duration/2:dt:duration/2;
arg = pi*freq*t;
wavelet = (1 - 2*arg.^2) .* exp(-arg.^2);
wavelet = wavelet / max(abs(wavelet)); % 归一化
end
分析说明
1. 自相关分析
- 目的：识别地震信号中的重复模式，检测周期性特征
- 数学表达： \(R_{xx}(τ)=∫_{−∞}^∞x(t)x(t+τ)dt\)
- 应用： 确定地震子波主频 检测多次波 估算地层厚度
2. 互相关分析
- 目的：测量两个地震道之间的相似性，检测相对时移
- 数学表达： \(R_{xy}(τ)=∫_{−∞}^∞x(t)y(t+τ)dt\)
- 应用： 速度分析 构造平移校正 相似断面提取
3. 相干分析
- 目的：量化多道数据间的相似性
- 计算方法： \(γ_{ij}^2=\frac{∣Pij∣^2}{P_{ii}P_{jj}}\)
- 应用： 断层检测 裂缝预测 岩性变化识别
图件解释
图1: 地震道波形
- 蓝色曲线：参考道（通常为中间道）
- 红色曲线：对比道（相邻道）
- 黑色虚线：典型地质界面标记
- 用途：直观比较不同道的波形特征
图2: 自相关分析
- 横轴：延迟时间（毫秒）
- 纵轴：归一化自相关系数
- 红圈标记：主要周期成分
- 用途：识别地震子波主频和地层周期特性
图3: 互相关分析
- 左图：原始波形对比
- 右图：互相关函数
- 黄色标记：最大相关位置和相关系数
- 用途：确定两道间的时移和相关性强度
图4: 多道自相关分析
- 网格布局：展示多条测线的自相关特征
- 用途：比较不同位置的地震信号特征变化
图5: 相干分析
- 上图：地震数据剖面（灰度表示振幅）
- 下图：沿测线相干系数变化
- 低相干区：对应断层、裂缝或岩性变化
- 用途：地质构造解释和异常检测
实际应用建议
-
数据选择： 选择信噪比较高的数据段 避免使用强能量干扰区域 考虑不同偏移距道集的分析
-
参数优化：
% 自适应参数调整 if dominant_freq > 50 fcuts = [10, 20, 100, 120]; % 高频数据 else fcuts = [2, 5, 40, 50]; % 低频数据 end -
高级分析扩展：
% 倾角扫描分析 dips = -5:0.5:5; % 倾角范围 dipCorrelations = zeros(length(dips), nSamples); for i = 1:length(dips) shiftedTrace = shift_trace(compareTrace, dips(i), t); [corr, ~] = xcorr(selectedTrace, shiftedTrace, 'coeff'); dipCorrelations(i, :) = corr; end -
三维可视化：
% 创建三维相干体 coherenceCube = zeros(nTraces, nSamples, nInlines); for inline = 1:nInlines coherenceCube(:,:,inline) = compute_inline_coherence(data(:,:,inline)); end % 可视化 volumeViewer(coherenceCube);
参考代码 对SGY地震数据进行自相关和互相关分析并生成图件 www.youwenfan.com/contentcno/82810.html
常见问题解决方案
-
低频噪音干扰：
- 增加高通滤波
- 使用小波去噪
% 小波去噪 [C, L] = wavedec(selectedTrace, 5, 'db4'); thr = wthrmngr('dw1ddenoLVL','penalhi',C,L,1); cleanedTrace = wdencmp('gbl',C,L,'db4',5,thr,'s'); -
计算效率优化：
- 使用并行计算
parfor i = 1:nTraces [acor(i,:), ~] = xcorr(analysisData(i,:), 'coeff'); end- 降采样处理
decimation_factor = 4; downsampledData = analysisData(:, 1:decimation_factor:end); -
结果验证方法： 合成数据测试 已知地质模型验证 井震标定对比
地质应用实例
- 断层检测： 相干切片中低值区对应断层 互相关时移量反映断层面倾角
- 薄层识别： 自相关峰值间距对应地层厚度 公式：h=2v⋅Δt
- 岩性预测： 高相干区：均匀沉积层 低相干区：裂缝发育带或岩性突变带
- 速度分析： 互相关最大峰值位置对应时差 叠加速度计算：θ=tan−1(ΔxΔt)
浙公网安备 33010602011771号