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 > nTracescompTrace = 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:numTracesToShowtraceIdx = 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-windowSizewindowData = 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 <= nSamplestrace(pos:pos+waveLength-1) = trace(pos:pos+waveLength-1) + reflectors(j)*baseWavelet;endend% 添加噪声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 > 50fcuts = [10, 20, 100, 120]; % 高频数据 elsefcuts = [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:nInlinescoherenceCube(:,:,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)