基于MSFM算法与最速下降法的射线追踪技术

基于MSFM算法与最速下降法的射线追踪技术

1、射线追踪效果

为了方便展示射线追踪的过程，用matlab将射线追踪过程做成了gif动图的形式，直观的展示一下跨孔方式的射线追踪效果，如下图所示，左图是波前走时等值线图（波前走时场），右图是射线路径图。

上图对应的速度场如下图

（上图用matlab画的较粗糙，没做调整，但速度场形态和位置都是正确的）

2、技术原理

具体的原理可以看文献，比如‘基于HAFMM的无射线追踪跨孔雷达走时层析成像’及其相关的文献，讲的很详细，我也是看的这些文献学习的。

3、代码

close all
clear 
clc

format long;

% 此程序是运行MSFM算法
% 用最速下降法得到接收点的走时与路径

disp('开始MSFM算法计算程函方程得到走时场...')
pause(1);

% 保存数据的文件名
outfile = '模型正演数据v03.dat';

% 初始化参数
% 网格单元个数
nex = 201;
ney = 301;

% 网格间距,单位为米。
dx = 0.1;
dy = 0.1;

% 总的网格单元个数
nexy = nex*ney;

% 给定初始速度模型
vp2d = zeros(ney,nex);

% 背景的速度
vp2d(1:ney,1:nex) = 2000;
% 异常体的速度
vp2d(120:200,50:90) = 4500;

vp2d(20:50,150:170) = 500;

% 显示速度模型
figure(1),pcolor(vp2d);shading flat;axis equal

% 发射点的坐标
% 发射点占横向网格单元的个数为X坐标。
% 发射点占纵向网格单元的个数为Y坐标。
% 钻孔两边设置30个发射点，T代表发射点，R代表接受点。

nPoint = 30;
xt = 1*ones(1,nPoint);
yt = linspace(6,296,nPoint);

xr = 201*ones(1,nPoint);
yr = linspace(6,296,nPoint);

tic
nRay = nPoint*nPoint;
Vmean = vp2d/dx;

% 发射点组成的矩阵
SourcePoints = [yt;xt];
% 接收点组成的矩阵
ReceivePoints = [yr;xr];

% 运行msfm2d程序，得到走时场，并提取接收点处的走时值。
% 每条射线的网格坐标及走时
TRpoints_firstTime = zeros(nRay,5);

tbegin = cputime;
for i = 1:30
    % i = 20;

    fprintf('\n');
    pause(0.5);  % 等待一秒，观察图像

    SourcePoint = SourcePoints(:,i);
    Traveltime = msfm2d(Vmean,SourcePoint,true,true);
    figure(2);
    set(gcf,'position',[100 100 1150 600]);
    subplot(1,2,1),contourf(Traveltime*1e3,20);
    colormap(jet),h = colorbar;
    h.Label.String  = 'Time(ms)';
    h.Label.FontWeight = 'bold';
    axis([-.5 202.5 -.5 302.5]);
    xlabel('Distance(m)','FontWeight','bold');
    ylabel('Depth(m)','FontWeight','bold');
    nx = linspace(0,201,5);
    ny = linspace(0,301,7);
    set(gca,'ydir','reverse','ytick',ny);
    set(gca,'yticklabel',{'0','-5','-10','-15','-20','-25','-30'});
    set(gca,'xtick',nx,'xticklabel',{'0','5','10','15','20'});
    set(gca,'xaxislocation','top');
    set(gca,'fontsize',16,'fontname','times new roman');
    firstTime(1+(i-1)*nPoint:i*nPoint) = Traveltime(yr,nex);
    jieshou = [xr',yr'];
    fashe_reverse = fliplr(SourcePoint');
    fashe = repmat(fashe_reverse,nPoint,1);
    TRtime = Traveltime(yr,nex);
    TRpoints_firstTime(1+(i-1)*nPoint:i*nPoint,:) = [fashe,jieshou,TRtime];
    tran = SourcePoint;
    % rec = [301,201];
    T = Traveltime;
    x_begin = 0;
    y_begin = -0.05;
    di = 0.1;
    for j = 1:30
        aaa = ReceivePoints(:,j);
        rec = aaa';
        [raylength,raypath] =steepest_descent_raytracing(T,tran,rec,x_begin,y_begin,di);
        subplot(1,2,2),hold on
        plot(raypath(1,:),raypath(2,:),'r-');
        xlabel('Distance(m)','FontWeight','bold');
        ylabel('Depth(m)','FontWeight','bold');
        set(gca,'xaxislocation','top');
        axis([-.3 20.3 -30.1 .1]);
        box on
        grid on
        set(gca,'fontsize',16,'fontname','times new roman');
    end
    
%     outputfile = strcat('CT',num2str(i),'.png');
% %     disp(outputfile);
%     print(outputfile,'-dpng');
    
    tend = cputime - tbegin;
    
%     if i == 1
%         tend = tend1;
%     else
%         tend = tend - tend1;
%     end
    fprintf('第 %d 个发射点计算时间为 %g 秒\n',i,tend);
%     fprintf('运行 %d 个发射点计算时间为 %g 秒\n',i,tend1);
end

% 将网格数据转化为坐标数据
TRpoints_firstTime(:,1:4) = TRpoints_firstTime(:,1:4)*dx - dx;
% 保存数据
dlmwrite(outfile,TRpoints_firstTime,...
    'delimiter',' ','precision','%.6f');
toc
fprintf('\n');
tend = cputime;
stime = tend - tbegin;
sout = sprintf('总耗时 %g 秒',stime);
disp(sout);

这段程序计算的结果不是很准确，但是我还未发现问题出在哪里，希望感兴趣的读者帮忙指出问题，我再修改一下。

上面的代码是是主程序代码，包含的两个子函数在网上可以找到，代码较长，我整理成文件夹上传到优快云资源，欢迎感兴趣的读者下载学习交流。

搬砖不易，走过路过的，点个赞支持一下，谢谢！