面向反射/折射初至走时的 2D 射线追踪对标井下/地面石油勘探场景多层介质 层速度模型 射线方程程函/哈密顿形式 界面 Snell 透射/反射 初至走时。一、问题建模先对齐术语石油勘探语境下射线追踪通常是速度建模 / 走时层析 / Kirchhoff 偏移的前置要素含义速度模型v(x,z)v(x,z)v(x,z)层状 / 线性梯度 / 网格射线方程drdsp, dpds∇(ln⁡v) ⇒\frac{d\boldsymbol{r}}{ds}\boldsymbol{p},\; \frac{d\boldsymbol{p}}{ds}\nabla(\ln v)\;\Rightarrowdsdr​p,dsdp​∇(lnv)⇒哈密顿形式界面条件Snell:sin⁡θ1v1sin⁡θ2v2\frac{\sin\theta_1}{v_1}\frac{\sin\theta_2}{v_2}v1​sinθ1​​v2​sinθ2​​透射反射同理目标炮-检走时ttt射线路径覆盖分析两套实现A 版层状模型 解析 SnellB 版网格速度模型 哈密顿 ODE RK4二、A 版层状模型 Snell 透射/反射2.1 速度模型定义velocity_model.m%% 层状速度模型石油勘探常用% 每层: [顶深(m), 底深(m), v_top(m/s), v_bottom(m/s)]% 层内速度线性插值 v(z) v_top (v_bottom-v_top)*(z-top)/(bottom-top)functionvlayered_velocity(z,layers)% z: 深度(m)可标量或向量% layers: N×4 矩阵vzeros(size(z));fori1:size(layers,1)toplayers(i,1);bottomlayers(i,2);v_toplayers(i,3);v_bottomlayers(i,4);mask(ztop)(zbottom);ifisize(layers,1)mask(ztop);% 最底层无下底endv(mask)v_top(v_bottom-v_top)*(z(mask)-top)/(bottom-topeps);endend%% 示例模型浅层沉积基底% 地表 z0深度layers[...0,800,1800,2200;% 表层泥岩800,1800,2200,2800;% 砂岩1800,3000,2800,3500;% 灰岩];2.2 单层内射线传播常梯度解析近似function[x2,z2,t12,pz2]ray_segment_layer(z1,z2_target,x1,pz1,vtop,vbot,layer_top,layer_bottom)% 单层内射线追踪线性梯度 v(z)v0 g*z% 输入: 起点(z1,x1), 目标深度 z2_target, 入射慢度 pz1, 层参数% 输出: 出点(x2,z2), 走时 t12, 出射慢度 pz2g(vbot-vtop)/(layer_bottom-layer_top);% 梯度v0vtop-g*layer_top;% v(z)v0g*z% 慢度 p sinθ/v, 水平慢度 px 守恒% 由 Snell: px pz1 * (??) —— 更稳的做法是用 Hamilton:% dz/ds pz*v, dpz/ds -dv/dz -g% 常梯度层内可解析:% 令 pz(z) pz1 - g*(z-z1)/v(z)? 不直接用慢度形式:% 改用慢度 u1/v 的形式更稳% 程函: |∇τ|² u², 射线: dx/dt u*p 之类% 工程上石油界常用梯形法则Snell这里给数值段通用pxpz1*0;% 水平慢度——等等得从入射角来% 重新入参: 改成给入射角更直观end上面这个函数写得太教科书实际石油界层内常用两段做法层内常速 → 直线 Snell层内线性梯度 → 圆弧形射线解析x(z)x(z)x(z)是sin⁡\sinsin/cos⁡\coscos形式2.3 主追踪器层状 透射/反射raytrace_2d.m%% 2D 射线追踪 —— 层状模型shooting 法% 炮 (xs, zs), 检波器 (xr, zr)多层常速层界面 Snell 透射% 这里做: 直达 折射临界 反射 三类function[rays]raytrace_2d_shot(xs,zs,xr_arr,zr_arr,layers)% xs,zs: 炮坐标% xr_arr,zr_arr: 检波器数组% layers: N×4 层定义nclength(xr_arr);rayscell(nc,1);forir1:nc xrxr_arr(ir);zrzr_arr(ir);% 先试: 直达波同一层内iffloor(zs/layers(1,2))floor(zr/layers(1,2))% 粗糙同层判vlayered_velocity((zszr)/2,layers);Lsqrt((xr-xs)^2(zr-zs)^2);t_directL/v;ray.typedirect;ray.path[[xs,zs];[xr,zr]];ray.ttt_direct;rays{ir}ray;continue;end% 透射波射线穿过界面% shooting: 从炮出发猜入射角追到检波器深度看水平偏差% 这里给简化版: 假设垂直层状射线只在界面处偏折层内直线% 找出炮→检经过哪些界面z_interface[layers(:,1);layers(end,2)];% 所有界面深度z_ptsunique(sort([zs;zr;z_interface]));z_ptsz_pts(z_ptsmin(zs,zr)z_ptsmax(zs,zr));% 每段: 常速, Snell 在界面处% 入射角 θ1: sinθ1/v1 sinθ2/v2 p (射线参数, 守恒)% shooting: 搜 p (水平慢度)p_candidateslinspace(1e-6,1/min(layers(:,3))*0.99,200);% p 1/v_maxbest_errinf;best_pNaN;best_path[];best_ttinf;forip1:length(p_candidates)pp_candidates(ip);% 水平慢度, 守恒% 从炮出发x_curxs;z_curzs;tt0;path[x_cur,z_cur];oktrue;% 向下追界面foriz2:length(z_pts)z_nextz_pts(iz);% 当前所在层ilayerfind(z_curlayers(:,1)z_cur[layers(2:end,1);inf],1);ifisempty(ilayer),ilayersize(layers,1);endv_layerlayered_velocity((z_curz_next)/2,layers);% 由 p 求垂直慢度: pz sqrt(u² - p²), u1/vu1/v_layer;ifpu okfalse;break;% 全反射/临界endpzsqrt(u^2-p^2);% 层内 dz, dxdzabs(z_next-z_cur);dsdz/(v_layer*pz*v_layer);% 错——ds dz/(v*cosθ), cosθ pz*u? 重算:% 慢度矢量 p⃗ (px, pz) (p, pz)% dx/ds v * px v*p, dz/ds v*pz% 所以 dz v*pz * ds → ds dz/(v*pz)% 但 pz √(u²-p²)dsdz/(v_layer*pz);dxv_layer*p*ds;ifz_nextz_cur,x_curx_curdx;elsex_curx_cur-dx;endz_curz_next;ttttds;path[path;x_cur,z_cur];endif~ok,continue;end% 看水平偏差errabs(x_cur-xr);iferrbest_err best_errerr;best_pp;best_pathpath;best_tttt;endend% 收敛判据ifbest_err5% 5m 以内算追到ray.typetransmitted;ray.pathbest_path;ray.ttbest_tt;ray.pbest_p;rays{ir}ray;elserays{ir}[];% 没追上可能是反射/临界这里略endendend2.4 主脚本seismic_raytrace_demo.m%% 石油勘探 2D 射线追踪 democlear;clc;close all;%% 1. 速度模型layers[...0,800,1800,2200;% 表层800,1800,2200,2800;% 砂岩1800,3000,2800,3500;% 灰岩];%% 2. 观测系统xs0;zs0;% 地表炮xr0:50:2000;zrzeros(size(xr));% 地表检波器% 也可以加井下: zr 1000*ones(size(xr));%% 3. 追踪raysraytrace_2d_shot(xs,zs,xr,zr,layers);%% 4. 走时曲线ttzeros(size(xr));fori1:length(xr)if~isempty(rays{i})tt(i)rays{i}.tt;elsett(i)NaN;endend%% 5. 可视化figure(Color,w,Position,[1001001200500]);% 速度剖面subplot(1,2,1);hold on;grid on;box on;% 画层fori1:size(layers,1)plot([02500],[layers(i,1)layers(i,1)],k-,LineWidth,1.5);plot([02500],[layers(i,2)layers(i,2)],k--,LineWidth,1);% 速度填充示意zzlinspace(layers(i,1),layers(i,2),100);vvlayered_velocity(zz,layers);yyvv*ones(size(zz))? 不对——换思路用 contourfend% 画几条射线colorsparula(length(rays));cnt0;fori1:5:length(rays)% 稀疏画if~isempty(rays{i})cntcnt1;rrays{i};plot(r.path(:,1),r.path(:,2),-,Color,colors(cnt,:),LineWidth,1);plot(r.path(:,1),r.path(:,2),.k,MarkerSize,3);endend% 炮检plot(xs,zs,r*,MarkerSize,14,LineWidth,2);plot(xr,zr,kv,MarkerSize,8,LineWidth,1);xlabel(x (m));ylabel(z (m));title(射线追踪层状模型);axis ij;ylim([30000]);xlim([-1002100]);% 走时曲线subplot(1,2,2);hold on;grid on;plot(xr,tt,bo-,LineWidth,1.5);xlabel(偏移距 (m));ylabel(走时 (s));title(初至走时曲线 (t-x));% 画 v(z) 示意axes(Position,[0.70.150.10.6]);hold on;zzlinspace(0,3000,200);vvlayered_velocity(zz,layers);plot(vv,zz,r-,LineWidth,2);set(gca,YDir,reverse);ylim([30000]);xlabel(v (m/s));ylabel();title(v(z));sgtitle(石油勘探 2D 射线追踪 Demo,FontSize,14,FontWeight,bold);参考代码 raytrace石油勘探的射线追踪程序www.youwenfan.com/contentcsw/82208.html三、B 版网格速度模型 哈密顿 ODEA 版够教学但真做走时层析一般用网格v(x,z)v(x,z)v(x,z) 最短路径法SPM或 Runge-Kutta 追射线。给你 RK 版骨架%% 哈密顿射线追踪网格速度模型% 射线方程慢度形式:% dx/dt v² * px% dz/dt v² * pz% dpx/dt -v * ∂v/∂x% dpz/dt -v * ∂v/∂z% 其中 p ∇τ慢度矢量, v 1/ufunction[x_ray,z_ray,t_ray]hamilton_raytrace(x0,z0,px0,pz0,v_grid,xvec,zvec,dt)% v_grid: 网格速度 (Nz×Nx)% xvec, zvec: 网格坐标% dt: RK 步长% 初始xx0;zz0;pxpx0;pzpz0;t0;x_rayx;z_rayz;t_rayt;whilet1e4% 追到检波器逻辑外面包% 插值 v, dv/dx, dv/dzvinterp2(xvec,zvec,v_grid,x,z,linear,2000);[dvdx,dvdz]gradient(v_grid,mean(diff(xvec)),mean(diff(zvec)));dvdx_iinterp2(xvec,zvec,dvdx,x,z,linear,0);dvdz_iinterp2(xvec,zvec,dvdz,x,z,linear,0);% RK4% k1k1xv^2*px;k1zv^2*pz;k1px-v*dvdx_i;k1pz-v*dvdz_i;% k2, k3, k4 略模板化写...% 这里给欧拉示意想稳就用 RK4xxdt*k1x;zzdt*k1z;pxpxdt*k1px;pzpzdt*k1pz;ttdt;x_ray[x_ray;x];z_ray[z_ray;z];t_ray[t_ray;t];% 出口判追到检波器ifzzvec(end)-10,break;endendend工程上更稳的是最短路径法SPM/ Fast Marching网格上直接解程函方程得 τ(x,z)再插射线∇τ 方向。石油界商业 codeHampson-Russell, Omega底层都是这个思路。