射线追踪与地震处理完整实现

1. 射线追踪实现

ray_tracing.m – 射线追踪核心算法

function [travel_times, ray_paths, incidence_angles] = ray_tracing(v_model, depths, offsets, method)
% 射线追踪算法
% 输入:
%   v_model - 速度模型 [v1, v2, ..., vn] (m/s)
%   depths - 界面深度 [z1, z2, ..., zn] (m)
%   offsets - 偏移距数组 (m)
%   method - 追踪方法 ('dix', 'kirchhoff', 'snell')
% 输出:
%   travel_times - 旅行时矩阵
%   ray_paths - 射线路径
%   incidence_angles - 入射角

    fprintf('开始射线追踪...\n');
    
    n_layers = length(v_model);
    n_offsets = length(offsets);
    
    % 初始化输出
    travel_times = zeros(n_offsets, n_layers);
    ray_paths = cell(n_offsets, n_layers);
    incidence_angles = zeros(n_offsets, n_layers);
    
    switch method
        case 'dix'
            % Dix方法射线追踪
            for i = 1:n_offsets
                for j = 1:n_layers
                    [travel_times(i,j), ray_paths{i,j}, incidence_angles(i,j)] = ...
                        dix_ray_tracing(v_model(1:j), depths(1:j), offsets(i));
                end
            end
            
        case 'snell'
            % Snell定律射线追踪
            for i = 1:n_offsets
                for j = 1:n_layers
                    [travel_times(i,j), ray_paths{i,j}, incidence_angles(i,j)] = ...
                        snell_ray_tracing(v_model(1:j), depths(1:j), offsets(i));
                end
            end
            
        otherwise
            error('不支持的射线追踪方法: %s', method);
    end
    
    fprintf('射线追踪完成\n');
end

function [t, ray_path, incidence_angle] = dix_ray_tracing(v, z, offset)
% Dix方法射线追踪
    n_layers = length(v);
    
    % 计算层厚度
    thickness = zeros(1, n_layers);
    thickness(1) = z(1);
    for i = 2:n_layers
        thickness(i) = z(i) - z(i-1);
    end
    
    % 计算均方根速度
    v_rms = zeros(1, n_layers);
    for i = 1:n_layers
        numerator = 0;
        denominator = 0;
        for j = 1:i
            numerator = numerator + v(j)^2 * thickness(j);
            denominator = denominator + thickness(j);
        end
        v_rms(i) = sqrt(numerator / denominator);
    end
    
    % 计算旅行时
    t = sqrt((offset / v_rms(end))^2 + (2 * sum(thickness) / v_rms(end))^2);
    
    % 简化射线路径
    ray_path.x = [0, offset/2, offset];
    ray_path.z = [0, sum(thickness)/2, 0];
    
    % 计算入射角
    incidence_angle = atan(offset / (2 * sum(thickness)));
end

function [t, ray_path, incidence_angle] = snell_ray_tracing(v, z, offset)
% Snell定律射线追踪
    n_layers = length(v);
    
    % 计算层厚度
    thickness = zeros(1, n_layers);
    thickness(1) = z(1);
    for i = 2:n_layers
        thickness(i) = z(i) - z(i-1);
    end
    
    % 使用迭代方法求解射线参数p
    p_min = 0;
    p_max = 1 / min(v);
    tolerance = 1e-6;
    max_iter = 100;
    
    p = (p_min + p_max) / 2;
    
    for iter = 1:max_iter
        % 计算当前p对应的偏移距
        x_calc = 0;
        for i = 1:n_layers
            denominator = sqrt(1/v(i)^2 - p^2);
            if imag(denominator) ~= 0
                denominator = 0;
            end
            x_calc = x_calc + thickness(i) * p * v(i) / denominator;
        end
        
        x_calc = 2 * x_calc;  % 往返路径
        
        if abs(x_calc - offset) < tolerance
            break;
        elseif x_calc < offset
            p_min = p;
        else
            p_max = p;
        end
        
        p = (p_min + p_max) / 2;
    end
    
    % 计算旅行时
    t = 0;
    for i = 1:n_layers
        denominator = sqrt(1/v(i)^2 - p^2);
        if imag(denominator) ~= 0
            denominator = 0;
        end
        t = t + 2 * thickness(i) / (v(i)^2 * denominator);
    end
    
    % 计算射线路径
    ray_path.x = [0];
    ray_path.z = [0];
    
    current_x = 0;
    current_z = 0;
    
    for i = 1:n_layers
        angle = asin(p * v(i));
        dx = thickness(i) * tan(angle);
        current_x = current_x + dx;
        current_z = current_z + thickness(i);
        
        ray_path.x = [ray_path.x, current_x];
        ray_path.z = [ray_path.z, current_z];
    end
    
    % 对称路径
    ray_path.x = [ray_path.x, 2*current_x];
    ray_path.z = [ray_path.z, 0];
    
    incidence_angle = asin(p * v(1));
end

2. 佐普利兹方程求解

zoeppritz_solver.m – 反射系数计算

function [Rpp, Rps, Tpp, Tps] = zoeppritz_solver(vp1, vs1, rho1, vp2, vs2, rho2, theta)
% 佐普利兹方程求解器
% 输入:
%   vp1, vs1, rho1 - 上层纵波速度、横波速度、密度
%   vp2, vs2, rho2 - 下层纵波速度、横波速度、密度
%   theta - 入射角 (度)
% 输出:
%   Rpp - PP反射系数
%   Rps - PS反射系数
%   Tpp - PP透射系数
%   Tps - PS透射系数

    % 转换为弧度
    theta = deg2rad(theta);
    
    % 计算透射角
    theta2 = asin(vp2/vp1 * sin(theta));  % 纵波透射角
    phi1 = asin(vs1/vp1 * sin(theta));    % 上层横波反射角
    phi2 = asin(vs2/vp1 * sin(theta));    % 下层横波透射角
    
    % 计算矩阵元素
    a = rho2 * (1 - 2*vs2^2/vp2^2 * sin(theta2)^2) - rho1 * (1 - 2*vs1^2/vp1^2 * sin(theta)^2);
    b = rho2 * (1 - 2*vs2^2/vp2^2 * sin(theta2)^2) + 2 * rho1 * vs1^2/vp1^2 * sin(theta)^2;
    c = rho1 * (1 - 2*vs1^2/vp1^2 * sin(theta)^2) + 2 * rho2 * vs2^2/vp2^2 * sin(theta2)^2;
    d = 2 * (rho2 * vs2^2 - rho1 * vs1^2) / vp1^2;
    
    E = (b * cos(theta)/vp1) + (c * cos(phi1)/vs1);
    F = (b * cos(phi2)/vs2) - (c * cos(theta2)/vp2);
    G = a - d * cos(theta)/vp1 * cos(phi2)/vs2;
    H = a - d * cos(phi1)/vs1 * cos(theta2)/vp2;
    
    D = E*F + G*H;
    
    % 计算反射和透射系数
    Rpp = (b * (E - F) - c * (G - H)) / D;
    Rps = -2 * cos(theta)/vp1 * (a*b + c*d * cos(theta2)/vp1 * cos(phi2)/vs2) / (vs1 * D);
    Tpp = 2 * rho1 * cos(theta)/vp1 * (F + H) / (vp1 * D);
    Tps = 2 * rho1 * cos(theta)/vp1 * (E - G) / (vs1 * D);
    
    % 确保结果为实数
    Rpp = real(Rpp);
    Rps = real(Rps);
    Tpp = real(Tpp);
    Tps = real(Tps);
end

function [Rpp_curve, angles] = zoeppritz_AVO(vp1, vs1, rho1, vp2, vs2, rho2, max_angle)
% 计算AVO响应曲线
% 输入:
%   max_angle - 最大入射角 (度)
    
    angles = 0:0.5:max_angle;
    Rpp_curve = zeros(size(angles));
    
    for i = 1:length(angles)
        [Rpp, ~, ~, ~] = zoeppritz_solver(vp1, vs1, rho1, vp2, vs2, rho2, angles(i));
        Rpp_curve(i) = Rpp;
    end
end

function plot_avo_response(Rpp_curve, angles, layer_info)
% 绘制AVO响应曲线
    
    figure('Position', [100, 100, 800, 600]);
    plot(angles, Rpp_curve, 'b-', 'LineWidth', 2);
    grid on;
    xlabel('入射角 (度)');
    ylabel('PP反射系数');
    title(sprintf('AVO响应曲线 - %s', layer_info));
    
    % 添加参考线
    hold on;
    plot([min(angles), max(angles)], [0, 0], 'k--', 'LineWidth', 1);
    
    % 标记特殊角度
    [min_val, min_idx] = min(Rpp_curve);
    [max_val, max_idx] = max(Rpp_curve);
    
    plot(angles(min_idx), min_val, 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'red');
    plot(angles(max_idx), max_val, 'go', 'MarkerSize', 8, 'MarkerFaceColor', 'green');
    
    legend('AVO曲线', '零线', '最小值', '最大值', 'Location', 'best');
end

3. 合成地震记录

synthetic_seismogram.m – 地震记录合成

function [seismic_data, time_axis] = synthetic_seismogram(v_model, density_model, depth_model, offsets, wavelet, dt, t_max)
% 合成地震记录
% 输入:
%   v_model - 速度模型
%   density_model - 密度模型
%   depth_model - 深度模型
%   offsets - 偏移距数组
%   wavelet - 子波
%   dt - 时间采样率
%   t_max - 最大时间
% 输出:
%   seismic_data - 地震数据矩阵
%   time_axis - 时间轴

    fprintf('开始合成地震记录...\n');
    
    n_offsets = length(offsets);
    n_times = round(t_max / dt) + 1;
    time_axis = 0:dt:t_max;
    
    seismic_data = zeros(n_times, n_offsets);
    
    % 射线追踪获取旅行时和入射角
    [travel_times, ~, incidence_angles] = ray_tracing(v_model, depth_model, offsets, 'snell');
    
    n_layers = length(v_model);
    
    % 为每个界面合成地震道
    for offset_idx = 1:n_offsets
        fprintf('处理偏移距 %.1f m...\n', offsets(offset_idx));
        
        for layer_idx = 2:n_layers  % 从第一个界面开始
            % 获取当前界面的物性参数
            vp1 = v_model(layer_idx-1);
            vs1 = vp1 / 1.7;  % 估算横波速度
            rho1 = density_model(layer_idx-1);
            
            vp2 = v_model(layer_idx);
            vs2 = vp2 / 1.7;
            rho2 = density_model(layer_idx);
            
            % 计算反射系数
            incidence_angle_deg = rad2deg(incidence_angles(offset_idx, layer_idx));
            Rpp = zoeppritz_solver(vp1, vs1, rho1, vp2, vs2, rho2, incidence_angle_deg);
            
            % 获取旅行时
            t_arrival = travel_times(offset_idx, layer_idx);
            
            % 将子波放置在到达时间
            if t_arrival <= t_max
                sample_idx = round(t_arrival / dt) + 1;
                wavelet_length = length(wavelet);
                
                start_idx = max(1, sample_idx - floor(wavelet_length/2));
                end_idx = min(n_times, sample_idx + floor(wavelet_length/2));
                
                wavelet_start = max(1, floor(wavelet_length/2) - (sample_idx - start_idx) + 1);
                wavelet_end = min(wavelet_length, floor(wavelet_length/2) + (end_idx - sample_idx) + 1);
                
                % 添加反射信号
                seismic_data(start_idx:end_idx, offset_idx) = seismic_data(start_idx:end_idx, offset_idx) + ...
                    Rpp * wavelet(wavelet_start:wavelet_end)';
            end
        end
    end
    
    fprintf('地震记录合成完成\n');
end

function wavelet = ricker_wavelet(freq, dt, length_ms)
% 生成Ricker子波
% 输入:
%   freq - 主频 (Hz)
%   dt - 时间采样率 (s)
%   length_ms - 子波长度 (ms)

    length_s = length_ms / 1000;
    t = -length_s/2:dt:length_s/2;
    
    % Ricker子波公式
    wavelet = (1 - 2 * (pi * freq * t).^2) .* exp(-(pi * freq * t).^2);
    
    % 归一化
    wavelet = wavelet / max(abs(wavelet));
end

function plot_seismic_section(seismic_data, time_axis, offsets, title_str)
% 绘制地震剖面
    
    figure('Position', [100, 100, 1000, 800]);
    
    % 地震数据显示
    subplot(2,1,1);
    imagesc(offsets, time_axis, seismic_data);
    colormap(seismic_colormap());
    colorbar;
    xlabel('偏移距 (m)');
    ylabel('时间 (s)');
    title(title_str);
    set(gca, 'YDir', 'reverse');
    
    % 叠加道显示
    subplot(2,1,2);
    stack_trace = mean(seismic_data, 2);
    plot(stack_trace, time_axis, 'b-', 'LineWidth', 1.5);
    set(gca, 'YDir', 'reverse');
    xlabel('振幅');
    ylabel('时间 (s)');
    title('叠加道');
    grid on;
end

function cmap = seismic_colormap()
% 地震数据显示颜色映射
    n_colors = 64;
    cmap = zeros(n_colors, 3);
    
    % 红色到白色到蓝色
    for i = 1:n_colors
        if i <= n_colors/2
            % 红色到白色
            cmap(i,:) = [1, i/(n_colors/2), i/(n_colors/2)];
        else
            % 白色到蓝色
            idx = i - n_colors/2;
            cmap(i,:) = [1 - idx/(n_colors/2), 1 - idx/(n_colors/2), 1];
        end
    end
end

4. 常规动校正(NMO)

nmo_correction.m – 常规动校正

function [nmo_corrected, stretch_factors] = nmo_correction(seismic_data, time_axis, offsets, v_nmo)
% 常规动校正
% 输入:
%   seismic_data - 输入地震数据
%   time_axis - 时间轴
%   offsets - 偏移距数组
%   v_nmo - NMO速度函数
% 输出:
%   nmo_corrected - NMO校正后的数据
%   stretch_factors - 拉伸因子

    fprintf('开始常规动校正...\n');
    
    [n_times, n_offsets] = size(seismic_data);
    nmo_corrected = zeros(size(seismic_data));
    stretch_factors = zeros(size(seismic_data));
    
    dt = time_axis(2) - time_axis(1);
    
    for offset_idx = 1:n_offsets
        offset = offsets(offset_idx);
        
        for time_idx = 1:n_times
            t0 = time_axis(time_idx);
            
            % 获取当前时间的NMO速度
            if isa(v_nmo, 'function_handle')
                v_current = v_nmo(t0);
            else
                v_current = v_nmo;
            end
            
            % 计算NMO校正时间
            t_nmo = sqrt(t0^2 + (offset / v_current)^2);
            
            % 找到最近的时间采样点
            nmo_idx = round(t_nmo / dt) + 1;
            
            if nmo_idx <= n_times
                nmo_corrected(time_idx, offset_idx) = seismic_data(nmo_idx, offset_idx);
                
                % 计算拉伸因子
                stretch_factors(time_idx, offset_idx) = t_nmo / t0;
            end
        end
        
        if mod(offset_idx, 10) == 0
            fprintf('  完成偏移距 %d/%d\n', offset_idx, n_offsets);
        end
    end
    
    fprintf('常规动校正完成\n');
end

function v_nmo_func = estimate_nmo_velocity(v_model, depth_model, time_axis)
% 估算NMO速度函数
% 使用Dix公式计算层状介质的NMO速度

    n_times = length(time_axis);
    v_nmo_func = zeros(1, n_times);
    
    % 计算均方根速度作为NMO速度
    for i = 1:n_times
        t_current = time_axis(i);
        
        % 找到当前时间对应的层位
        layer_idx = 1;
        t_accumulated = 0;
        
        while layer_idx <= length(v_model) && t_accumulated < t_current/2
            if layer_idx == 1
                layer_time = depth_model(1) / v_model(1);
            else
                layer_time = (depth_model(layer_idx) - depth_model(layer_idx-1)) / v_model(layer_idx);
            end
            
            if t_accumulated + layer_time > t_current/2
                % 部分层位
                remaining_time = t_current/2 - t_accumulated;
                break;
            else
                t_accumulated = t_accumulated + layer_time;
                layer_idx = layer_idx + 1;
            end
        end
        
        % 计算均方根速度
        if layer_idx > 1
            v_rms = 0;
            for j = 1:layer_idx-1
                if j == 1
                    layer_thickness = depth_model(1);
                else
                    layer_thickness = depth_model(j) - depth_model(j-1);
                end
                v_rms = v_rms + v_model(j)^2 * (layer_thickness / v_model(j));
            end
            v_rms = sqrt(v_rms / t_accumulated);
        else
            v_rms = v_model(1);
        end
        
        v_nmo_func(i) = v_rms;
    end
end

5. 无拉伸动校正

stretch_free_nmo.m – 无拉伸动校正

function [nmo_stretch_free, correction_quality] = stretch_free_nmo(seismic_data, time_axis, offsets, v_nmo, max_stretch)
% 无拉伸动校正
% 输入:
%   seismic_data - 输入地震数据
%   time_axis - 时间轴
%   offsets - 偏移距数组
%   v_nmo - NMO速度函数
%   max_stretch - 最大允许拉伸因子
% 输出:
%   nmo_stretch_free - 无拉伸NMO校正后的数据
%   correction_quality - 校正质量指标

    fprintf('开始无拉伸动校正...\n');
    
    [n_times, n_offsets] = size(seismic_data);
    nmo_stretch_free = zeros(size(seismic_data));
    correction_quality = zeros(size(seismic_data));
    
    dt = time_axis(2) - time_axis(1);
    
    for offset_idx = 1:n_offsets
        offset = offsets(offset_idx);
        
        for time_idx = 1:n_times
            t0 = time_axis(time_idx);
            
            % 获取当前时间的NMO速度
            if isa(v_nmo, 'function_handle')
                v_current = v_nmo(t0);
            else
                v_current = v_nmo;
            end
            
            % 计算NMO校正时间
            t_nmo = sqrt(t0^2 + (offset / v_current)^2);
            
            % 计算拉伸因子
            stretch_factor = t_nmo / t0;
            
            if stretch_factor <= max_stretch
                % 拉伸在允许范围内，使用常规NMO
                nmo_idx = round(t_nmo / dt) + 1;
                
                if nmo_idx <= n_times
                    nmo_stretch_free(time_idx, offset_idx) = seismic_data(nmo_idx, offset_idx);
                    correction_quality(time_idx, offset_idx) = 1;  % 优质校正
                end
            else
                % 拉伸超过限制，使用插值或拒绝该采样
                % 这里使用线性插值来减少拉伸影响
                t_nmo_sample = t_nmo / dt + 1;
                
                if t_nmo_sample < n_times
                    idx_floor = floor(t_nmo_sample);
                    idx_ceil = ceil(t_nmo_sample);
                    
                    weight_ceil = t_nmo_sample - idx_floor;
                    weight_floor = 1 - weight_ceil;
                    
                    nmo_stretch_free(time_idx, offset_idx) = ...
                        weight_floor * seismic_data(idx_floor, offset_idx) + ...
                        weight_ceil * seismic_data(idx_ceil, offset_idx);
                    
                    % 质量指标：拉伸因子的倒数
                    correction_quality(time_idx, offset_idx) = 1 / stretch_factor;
                end
            end
        end
        
        if mod(offset_idx, 10) == 0
            fprintf('  完成偏移距 %d/%d\n', offset_idx, n_offsets);
        end
    end
    
    fprintf('无拉伸动校正完成\n');
end

function [optimized_vnmo, residuals] = optimize_nmo_velocity(seismic_data, time_axis, offsets, v_initial)
% 优化NMO速度
% 使用扫描方法找到最优NMO速度
    
    fprintf('开始NMO速度优化...\n');
    
    v_range = 0.8 * v_initial : 50 : 1.2 * v_initial;
    n_velocities = length(v_range);
    
    stack_quality = zeros(1, n_velocities);
    
    for v_idx = 1:n_velocities
        v_current = v_range(v_idx);
        
        % 应用NMO校正
        nmo_corrected = nmo_correction(seismic_data, time_axis, offsets, v_current);
        
        % 计算叠加道
        stack_trace = mean(nmo_corrected, 2);
        
        % 使用叠加道能量作为质量指标
        stack_quality(v_idx) = sum(stack_trace.^2);
    end
    
    % 找到最佳速度
    [~, best_idx] = max(stack_quality);
    optimized_vnmo = v_range(best_idx);
    
    residuals = stack_quality;
    
    fprintf('NMO速度优化完成，最佳速度: %.1f m/s\n', optimized_vnmo);
    
    % 绘制速度扫描结果
    figure('Position', [100, 100, 800, 400]);
    plot(v_range, stack_quality, 'b-', 'LineWidth', 2);
    hold on;
    plot(optimized_vnmo, max(stack_quality), 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'red');
    xlabel('NMO速度 (m/s)');
    ylabel('叠加质量');
    title('NMO速度扫描');
    grid on;
    legend('质量指标', '最佳速度', 'Location', 'best');
end

6. 主程序示例

main_seismic_processing.m – 主程序

%% 地震数据处理完整流程演示
clear; clc; close all;

fprintf('=== 地震数据处理完整流程 ===\n\n');

%% 1. 设置模型参数
fprintf('1. 设置模型参数...\n');

% 速度模型 (m/s)
v_model = [1500, 2000, 2500, 3000, 3500];

% 密度模型 (kg/m³)
density_model = [2000, 2100, 2200, 2300, 2400];

% 界面深度 (m)
depth_model = [500, 1000, 1500, 2000, 2500];

% 偏移距范围 (m)
offsets = 0:50:3000;

fprintf('   速度模型: %s m/s\n', mat2str(v_model));
fprintf('   界面深度: %s m\n', mat2str(depth_model));
fprintf('   偏移距范围: 0 - %d m, %d 道\n', max(offsets), length(offsets));

%% 2. 射线追踪
fprintf('\n2. 射线追踪...\n');

[travel_times, ray_paths, incidence_angles] = ray_tracing(...
    v_model, depth_model, offsets, 'snell');

% 绘制射线路径
figure('Position', [100, 100, 1200, 800]);
subplot(2,2,1);
for i = 1:min(10, length(offsets))
    plot(ray_paths{i,end}.x, ray_paths{i,end}.z, 'b-', 'LineWidth', 1);
    hold on;
end
set(gca, 'YDir', 'reverse');
xlabel('水平距离 (m)');
ylabel('深度 (m)');
title('射线路径');
grid on;

% 绘制旅行时曲线
subplot(2,2,2);
for i = 1:length(v_model)
    plot(offsets, travel_times(:,i), 'LineWidth', 2);
    hold on;
end
xlabel('偏移距 (m)');
ylabel('旅行时 (s)');
title('旅行时曲线');
legend('第1层', '第2层', '第3层', '第4层', '第5层', 'Location', 'best');
grid on;

%% 3. AVO分析
fprintf('\n3. AVO分析...\n');

% 选择第一个界面进行分析
vp1 = v_model(1);
vs1 = vp1 / 1.7;
rho1 = density_model(1);

vp2 = v_model(2);
vs2 = vp2 / 1.7;
rho2 = density_model(2);

% 计算AVO响应
max_angle = 50;
[Rpp_curve, angles] = zoeppritz_AVO(vp1, vs1, rho1, vp2, vs2, rho2, max_angle);

% 绘制AVO曲线
subplot(2,2,3);
plot_avo_response(Rpp_curve, angles, '第一界面');

%% 4. 合成地震记录
fprintf('\n4. 合成地震记录...\n');

% 生成子波
dt = 0.002;  % 2ms采样
t_max = 4;   % 4秒记录长度
freq = 30;   % 30Hz主频

wavelet = ricker_wavelet(freq, dt, 100);

% 合成地震记录
[seismic_data, time_axis] = synthetic_seismogram(...
    v_model, density_model, depth_model, offsets, wavelet, dt, t_max);

% 显示合成记录
subplot(2,2,4);
plot_seismic_section(seismic_data, time_axis, offsets, '合成地震记录');

%% 5. 常规动校正
fprintf('\n5. 常规动校正...\n');

% 估算NMO速度
v_nmo_func = estimate_nmo_velocity(v_model, depth_model, time_axis);

% 应用常规NMO
[nmo_corrected, stretch_factors] = nmo_correction(...
    seismic_data, time_axis, offsets, v_nmo_func);

% 显示NMO校正结果
figure('Position', [100, 100, 1200, 800]);
subplot(2,2,1);
plot_seismic_section(seismic_data, time_axis, offsets, '原始CMP道集');

subplot(2,2,2);
plot_seismic_section(nmo_corrected, time_axis, offsets, '常规NMO校正后');

% 显示拉伸因子
subplot(2,2,3);
imagesc(offsets, time_axis, stretch_factors);
colorbar;
xlabel('偏移距 (m)');
ylabel('时间 (s)');
title('拉伸因子');
set(gca, 'YDir', 'reverse');
caxis([1, 2]);

%% 6. 无拉伸动校正
fprintf('\n6. 无拉伸动校正...\n');

max_stretch = 1.3;  % 最大允许拉伸30%

[nmo_stretch_free, correction_quality] = stretch_free_nmo(...
    seismic_data, time_axis, offsets, v_nmo_func, max_stretch);

% 显示无拉伸NMO结果
subplot(2,2,4);
plot_seismic_section(nmo_stretch_free, time_axis, offsets, '无拉伸NMO校正后');

%% 7. 叠加效果比较
fprintf('\n7. 叠加效果比较...\n');

figure('Position', [100, 100, 1000, 600]);

% 计算叠加道
stack_original = mean(seismic_data, 2);
stack_nmo = mean(nmo_corrected, 2);
stack_stretch_free = mean(nmo_stretch_free, 2);

% 显示叠加道比较
subplot(1,2,1);
plot(stack_original, time_axis, 'k-', 'LineWidth', 1, 'DisplayName', '原始叠加');
hold on;
plot(stack_nmo, time_axis, 'r-', 'LineWidth', 2, 'DisplayName', '常规NMO叠加');
plot(stack_stretch_free, time_axis, 'b-', 'LineWidth', 2, 'DisplayName', '无拉伸NMO叠加');
set(gca, 'YDir', 'reverse');
xlabel('振幅');
ylabel('时间 (s)');
title('叠加道比较');
legend('Location', 'best');
grid on;

% 计算叠加质量
quality_nmo = sum(stack_nmo.^2) / sum(stack_original.^2);
quality_stretch_free = sum(stack_stretch_free.^2) / sum(stack_original.^2);

subplot(1,2,2);
bar([1, 2], [quality_nmo, quality_stretch_free]);
set(gca, 'XTickLabel', {'常规NMO', '无拉伸NMO'});
ylabel('相对能量');
title('叠加质量比较');
grid on;

fprintf('   常规NMO叠加质量: %.4f\n', quality_nmo);
fprintf('   无拉伸NMO叠加质量: %.4f\n', quality_stretch_free);

%% 8. NMO速度优化
fprintf('\n8. NMO速度优化...\n');

% 使用第一个反射层的时间进行速度优化
target_time = travel_times(1, 2);  % 第一个界面的零偏移距旅行时
v_initial = v_nmo_func(find(time_axis >= target_time, 1));

[optimized_vnmo, residuals] = optimize_nmo_velocity(...
    seismic_data, time_axis, offsets, v_initial);

fprintf('\n=== 处理完成 ===\n');

% 保存结果
save('seismic_processing_results.mat', ...
    'seismic_data', 'nmo_corrected', 'nmo_stretch_free', ...
    'time_axis', 'offsets', 'v_model', 'travel_times');

7. 高级功能扩展

advanced_nmo.m – 高级NMO功能

function [nmo_advanced, velocity_field] = advanced_nmo(seismic_data, time_axis, offsets, method, parameters)
% 高级NMO校正方法
% 支持多种高级NMO技术
    
    switch method
        case 'anisotropic'
            % 各向异性NMO
            [nmo_advanced, velocity_field] = anisotropic_nmo(...
                seismic_data, time_axis, offsets, parameters);
            
        case 'multiples'
            % 多次波压制NMO
            [nmo_advanced, velocity_field] = multiple_suppression_nmo(...
                seismic_data, time_axis, offsets, parameters);
            
        case 'high_resolution'
            % 高分辨率NMO
            [nmo_advanced, velocity_field] = high_resolution_nmo(...
                seismic_data, time_axis, offsets, parameters);
            
        otherwise
            error('不支持的NMO方法: %s', method);
    end
end

function [nmo_aniso, velocity_field] = anisotropic_nmo(seismic_data, time_axis, offsets, parameters)
% 各向异性NMO校正
    
    v_nmo = parameters.v_nmo;
    eta = parameters.eta;  % 各向异性参数
    
    [n_times, n_offsets] = size(seismic_data);
    nmo_aniso = zeros(size(seismic_data));
    velocity_field = zeros(size(seismic_data));
    
    dt = time_axis(2) - time_axis(1);
    
    for offset_idx = 1:n_offsets
        offset = offsets(offset_idx);
        
        for time_idx = 1:n_times
            t0 = time_axis(time_idx);
            
            % 各向异性NMO公式
            t_nmo = t0 * sqrt(1 + (offset / (v_nmo * t0))^2 - ...
                2 * eta * (offset / (v_nmo * t0))^4 / (1 + (1 + 2*eta) * (offset / (v_nmo * t0))^2));
            
            nmo_idx = round(t_nmo / dt) + 1;
            
            if nmo_idx <= n_times
                nmo_aniso(time_idx, offset_idx) = seismic_data(nmo_idx, offset_idx);
            end
            
            velocity_field(time_idx, offset_idx) = v_nmo;
        end
    end
end

参考代码 射线追踪 www.youwenfan.com/contentzhf/63651.html

这个完整的实现包含了从射线追踪到动校正的完整地震处理流程，特别强调了无拉伸动校正技术来保持波形特征。所有代码都经过模块化设计，便于理解和扩展。