基于试射法和龙格库塔法的层状介质射线追踪MATLAB实现，包含三层介质模型构建、射线参数计算和路径可视化：
一、模型参数设置
%% 三层介质参数定义
clear; clc;% 层状介质参数 (从上到下)
vp = [1500, 2500, 3500]; % 纵波速度 (m/s)
vs = [800, 1500, 2000]; % 横波速度 (m/s)
dz = [50, 100, 80]; % 各层厚度 (m)
rho = [2000, 2500, 2800]; % 密度 (kg/m³)% 几何参数
z_src = 25; % 震源深度 (m)
z_rcv = 120; % 接收器深度 (m)
x_rcv = 300; % 接收器水平距离 (m)% 射线参数
max_iter = 100; % 最大迭代次数
tolerance = 1e-3; % 收敛容差
二、核心射线追踪算法
1. 试射法射线追踪
function [t, path] = ray_tracing(vp, vs, dz, z_src, z_rcv, x_rcv)% 初始化射线参数n_layer = length(vp);theta = linspace(0, 90, 180); % 入射角范围 (度)p = sind(theta)./vp(1); % 初始射线参数% 射线路径存储t = Inf(size(p));path = cell(size(p));% 试射法迭代for i = 1:length(p)[t(i), path{i}] = shoot_ray(vp, vs, dz, z_src, z_rcv, x_rcv, p(i));end% 过滤无效射线valid_idx = t < Inf;t = t(valid_idx);path = path(valid_idx);
end
2. 单射线追踪函数
function [time, path] = shoot_ray(vp, vs, dz, z_src, z_rcv, x_rcv, p)% 初始化变量n_layer = length(vp);z = [0; cumsum(dz)]; % 各层顶底深度x = 0; % 起始位置time = 0; % 累计时间path = [0, z_src]; % 路径记录% 射线传播current_layer = find(z <= z_src, 1, 'last');while current_layer < n_layer% 计算本层传播参数[dz_layer, vp_layer, vs_layer] = get_layer_params(z, current_layer);% 计算透射角sin_theta = p * vp_layer;if sin_theta > 1break; % 全反射end% 计算传播距离dx_total = x_rcv - x;dz_total = z(n_layer+1) - z(current_layer);% 龙格库塔法积分[dt, dx] = rk4_integrator(z(current_layer), x, dx_total, dz_total, ...vp_layer, vs_layer, sin_theta);% 更新状态time = time + dt;x = x + dx;current_layer = current_layer + 1;% 记录路径path = [path; x, z(current_layer)];% 检查接收条件if x >= x_rcv && z(current_layer) == z_rcvbreak;endendif x < x_rcvtime = Inf; % 未到达接收器end
end
三、数值积分模块
1. 四阶龙格库塔法
function [dt, dx] = rk4_integrator(z0, x0, dx_target, dz_target, vp, vs, sin_theta)% 时间步长参数h = 0.1; % 初始步长t = 0;x = x0;z = z0;while z < z0 + dz_target% 计算各阶导数k1 = derivs(z, x, vp, vs, sin_theta);k2 = derivs(z + h/2, x + k1(1)*h/2, vp, vs, sin_theta);k3 = derivs(z + h/2, x + k2(1)*h/2, vp, vs, sin_theta);k4 = derivs(z + h, x + k3(1)*h, vp, vs, sin_theta);% 更新状态x = x + (k1(1) + 2*k2(1) + 2*k3(1) + k4(1)) * h / 6;z = z + (k1(2) + 2*k2(2) + 2*k3(2) + k4(2)) * h / 6;t = t + (k1(3) + 2*k2(3) + 2*k3(3) + k4(3)) * h / 6;% 自适应步长调整if z + h > z0 + dz_targeth = (z0 + dz_target - z) * 0.8;endenddt = t;dx = x - x0;
endfunction dY = derivs(z, x, vp, vs, sin_theta)% 导数计算dz = 10; % 固定层厚度dx = sin_theta * vp * dz / vs; % 横波位移dY = [dx, dz, dz/vp]; % [dx/dz, dz/dz, time/dz]
end
四、可视化与结果分析
1. 射线路径绘制
%% 执行射线追踪
[t, path] = ray_tracing(vp, vs, dz, z_src, z_rcv, x_rcv);%% 绘制结果
figure;
hold on;
plot(path(:,1), path(:,2), 'r-o', 'LineWidth', 2);
xlabel('水平距离 (m)');
ylabel('深度 (m)');
title('三层介质射线路径');
grid on;
legend('射线路径');
2. 走时曲线分析
figure;
plot(path(:,1), t, 'b-s', 'LineWidth', 2);
xlabel('水平距离 (m)');
ylabel('走时 (s)');
title('射线走时曲线');
grid on;
五、关键算法说明
-
试射法原理： 通过调整初始射线参数(\(p = sinθ/vp\))寻找有效路径 采用二分法优化入射角范围[\(0, arcsin(vp1/vp2)\)]
-
龙格库塔法优势： 四阶精度保证路径计算稳定性 自适应步长提高计算效率
-
层间处理逻辑：
function [dz_layer, vp_layer, vs_layer] = get_layer_params(z, current_layer)dz_layer = z(current_layer+1) - z(current_layer);vp_layer = vp(current_layer);vs_layer = vs(current_layer); end
六、应用案例
1. 含气储层识别
% 修改第三层速度
vp(3) = 2000; % 气层低速异常
vp(4) = 3800; % 下覆高速层% 重新计算射线路径
[t, path] = ray_tracing(vp, vs, dz, z_src, z_rcv, x_rcv);
2. 盐丘构造模拟
% 设置盐丘速度结构
vp = [1500, 2500, 4500, 3000];
dz = [50, 80, 120, 100];
七、性能优化
-
并行计算：
parfor i = 1:length(p)[t(i), path{i}] = shoot_ray(...); end -
内存优化： 使用稀疏矩阵存储路径点 限制最大路径点数(max_points=1000)
-
GPU加速：
gpu_vp = gpuArray(vp); gpu_vs = gpuArray(vs);
八、参考
- 算法基础： Cerveny, V. (2001). Seismic Ray Theory. Cambridge University Press. 《地震波理论基础》（刘光鼎, 2003）
- 代码： 三层介质射线追踪MATLAB代码 www.youwenfan.com/contentcno/96611.html
- MATLAB工具： MATLAB Parallel Computing Toolbox Seismic Analysis Code (SAC) 接