已知走时场和炮点检波点的位置做射线追踪的方法，注意各向同性和各向异性介质有所不同

一、各向同性介质

在各向同性介质中，群速度方向等于相速度方向（走时场的梯度方向），因此，可以简单地沿走时场的梯度方向追踪到射线路径。

function ray_tracing_example()
    % 参数设置
    Lx = 10;  % x方向长度 (km)
    Lz = 10;  % z方向长度 (km)
    Nx = 101; % x方向网格点数
    Nz = 101; % z方向网格点数
    dx = Lx / (Nx - 1);  % x方向网格间距 (km)
    dz = Lz / (Nz - 1);  % z方向网格间距 (km)

    % 生成网格
    x = linspace(0, Lx, Nx);
    z = linspace(0, Lz, Nz);
    [X, Z] = meshgrid(x, z);

    % 生成示例走时场（假设波速均匀）
    v = 3;  % 波速 (km/s)
    T = sqrt(X.^2 + Z.^2) / v;

    % 炮点和多个检波点位置
    source = [0, 0];   % 炮点位置 (x, z)
    receivers = [7, 7; 3, 7; 5, 8; 9, 2; 6, 9; 8, 4]; % 多个检波点位置 (x, z)

    % 进行射线追踪
    figure;
    contourf(X, Z, T, 20); hold on;
    plot(source(1), source(2), 'go', 'MarkerFaceColor', 'g'); % 炮点
    for i = 1:size(receivers, 1)
        receiver = receivers(i, :);
        [ray_x, ray_z] = ray_tracing(T, source, receiver, dx, dz);
        plot(ray_x, ray_z, 'r', 'LineWidth', 2); % 绘制射线路径
        plot(receiver(1), receiver(2), 'bo', 'MarkerFaceColor', 'b'); % 检波点
    end
    xlabel('距离 (km)');
    ylabel('深度 (km)');
    title('射线追踪');
    colorbar;
    hold off;
end

function [ray_x, ray_z] = ray_tracing(T, source, receiver, dx, dz)
    % 初始化
    current_pos = receiver;
    ray_x = current_pos(1);
    ray_z = current_pos(2);
    
    % 梯度计算
    [dTdx, dTdz] = gradient(T, dx, dz);

    while norm(current_pos - source) > max(dx, dz)
        % 获取当前点的索引
        ix = round(current_pos(1) / dx) + 1;
        iz = round(current_pos(2) / dz) + 1;

        % 获取梯度
        grad_T = [dTdx(iz, ix), dTdz(iz, ix)];

        % 计算步进方向（沿梯度的反方向）
        step_direction = -grad_T / norm(grad_T);

        % 更新位置
        current_pos = current_pos + step_direction * min(dx, dz);

        % 记录射线位置
        ray_x(end + 1) = current_pos(1);
        ray_z(end + 1) = current_pos(2);
    end
end

三维代码

function ray_tracing_example_3D()
    % 参数设置
    Lx = 10;  % x方向长度 (km)
    Ly = 10;  % y方向长度 (km)
    Lz = 10;  % z方向长度 (km)
    Nx = 101; % x方向网格点数
    Ny = 101; % y方向网格点数
    Nz = 101; % z方向网格点数
    dx = Lx / (Nx - 1);  % x方向网格间距 (km)
    dy = Ly / (Ny - 1);  % y方向网格间距 (km)
    dz = Lz / (Nz - 1);  % z方向网格间距 (km)

    % 生成网格
    x = linspace(0, Lx, Nx);
    y = linspace(0, Ly, Ny);
    z = linspace(0, Lz, Nz);
    [X, Y, Z] = meshgrid(x, y, z);

    % 生成示例走时场（假设波速均匀）
    v = 3;  % 波速 (km/s)
    T = sqrt(X.^2 + Y.^2 + Z.^2) / v;

    % 炮点和多个检波点位置
    source = [0, 0, 0];   % 炮点位置 (x, y, z)
    receivers = [7, 7, 7; 3, 5, 5; 5, 8, 8; 9, 2, 2; 6, 9, 9; 8, 4, 4]; % 多个检波点位置 (x, y, z)

    % 进行射线追踪
    figure;
    hold on;

    % 先绘制走时场切片图，并设置透明度
    h = slice(X, Y, Z, T, [5], [5], [3, 6]);
    set(h, 'FaceAlpha', 0.2); % 设置透明度
    shading interp;

    % 绘制射线路径和点
    for i = 1:size(receivers, 1)
        receiver = receivers(i, :);
        [ray_x, ray_y, ray_z] = ray_tracing_3D(T, source, receiver, dx, dy, dz);
        plot3(ray_x, ray_y, ray_z, 'r', 'LineWidth', 2); % 绘制射线路径
        plot3(receiver(1), receiver(2), receiver(3), 'bo', 'MarkerFaceColor', 'b'); % 检波点
    end
    plot3(source(1), source(2), source(3), 'go', 'MarkerFaceColor', 'g'); % 炮点

    xlabel('距离 (km)');
    ylabel('距离 (km)');
    zlabel('深度 (km)');
    xlim([0, Lx]);
    ylim([0, Ly]);
    zlim([0, Lz]);
    title('三维射线追踪');
    colorbar;
    hold off;
end

function [ray_x, ray_y, ray_z] = ray_tracing_3D(T, source, receiver, dx, dy, dz)
    % 初始化
    current_pos = receiver;
    ray_x = current_pos(1);
    ray_y = current_pos(2);
    ray_z = current_pos(3);
    
    % 梯度计算
    [dTdx, dTdy, dTdz] = gradient(T, dx, dy, dz);

    while norm(current_pos - source) > max([dx, dy, dz])
        % 获取当前点的索引
        ix = round(current_pos(1) / dx) + 1;
        iy = round(current_pos(2) / dy) + 1;
        iz = round(current_pos(3) / dz) + 1;

        % 边界检查
        if ix < 1 || ix > size(T, 2) || iy < 1 || iy > size(T, 1) || iz < 1 || iz > size(T, 3)
            break;
        end

        % 获取梯度
        grad_T = [dTdx(iy, ix, iz), dTdy(iy, ix, iz), dTdz(iy, ix, iz)];

        % 计算步进方向（沿梯度的反方向）
        step_direction = -grad_T / norm(grad_T);

        % 更新位置
        current_pos = current_pos + step_direction * min([dx, dy, dz]);

        % 记录射线位置
        ray_x(end + 1) = current_pos(1);
        ray_y(end + 1) = current_pos(2);
        ray_z(end + 1) = current_pos(3);
    end
end

二、各向异性介质

而在各向异性介质中进行射线追踪时，需要沿着群速度方向，而不是简单地沿走时场的梯度方向。我们需要计算波的群速度，并沿群速度方向进行追踪。

function anisotropic_ray_tracing_example()
    % 参数设置
    Lx = 10;  % x方向长度 (km)
    Lz = 10;  % z方向长度 (km)
    Nx = 101; % x方向网格点数
    Nz = 101; % z方向网格点数
    dx = Lx / (Nx - 1);  % x方向网格间距 (km)
    dz = Lz / (Nz - 1);  % z方向网格间距 (km)

    % 生成网格
    x = linspace(0, Lx, Nx);
    z = linspace(0, Lz, Nz);
    [X, Z] = meshgrid(x, z);

    % 生成示例走时场（假设各向异性介质）
    epsilon = 0.4;
    delta = 0.1;
    v0 = 3;  % 波速 (km/s)
    T = generate_anisotropic_traveltime(X, Z, v0, epsilon, delta);

    % 炮点和多个检波点位置
    source = [0, 0];   % 炮点位置 (x, z)
    receivers = [7, 7; 3, 7; 5, 8; 9, 2; 6, 9; 8, 4]; % 多个检波点位置 (x, z)


    % 进行射线追踪
    figure;
    contourf(X, Z, T, 20); hold on;
    plot(source(1), source(2), 'go', 'MarkerFaceColor', 'g'); % 炮点
    for i = 1:size(receivers, 1)
        receiver = receivers(i, :);
        [ray_x, ray_z] = anisotropic_ray_tracing(T, source, receiver, v0, epsilon, delta, dx, dz);
        plot(ray_x, ray_z, 'r', 'LineWidth', 2); % 绘制射线路径
        plot(receiver(1), receiver(2), 'bo', 'MarkerFaceColor', 'b'); % 检波点
    end
    xlabel('距离 (km)');
    ylabel('深度 (km)');
    title('各向异性介质中的射线追踪');
    colorbar;
    hold off;
end

function T = generate_anisotropic_traveltime(X, Z, v0, epsilon, delta)
    % 生成各向异性介质中的走时场
    [Nz, Nx] = size(X);
    T = zeros(Nz, Nx);
    for i = 1:Nx
        for j = 1:Nz
            theta = atan2(Z(j, i), X(j, i));
            v = v0 * sqrt(1 + 2 * epsilon * sin(theta)^2 + 2 * delta * cos(theta)^2);
            T(j, i) = sqrt(X(j, i)^2 + Z(j, i)^2) / v;
        end
    end
end

function [ray_x, ray_z] = anisotropic_ray_tracing(T, source, receiver, v0, epsilon, delta, dx, dz)
    % 初始化
    current_pos = receiver;
    ray_x = current_pos(1);
    ray_z = current_pos(2);

    max_iter = 1000;
    tol = max(dx, dz);

    for iter = 1:max_iter
        % 获取当前点的索引
        ix = round(current_pos(1) / dx) + 1;
        iz = round(current_pos(2) / dz) + 1;

        % 边界检查
        if ix < 1 || ix > size(T, 2) || iz < 1 || iz > size(T, 1)
            break;
        end

        % 计算当前点的速度和群速度方向
        theta = atan2(current_pos(2), current_pos(1));
        v = v0 * sqrt(1 + 2 * epsilon * sin(theta)^2 + 2 * delta * cos(theta)^2);
        vg_x = v * cos(theta);
        vg_z = v * sin(theta);

        % 更新位置
        step_length = tol;
        step_direction = [vg_x, vg_z] / norm([vg_x, vg_z]);
        current_pos = current_pos - step_length * step_direction;

        % 记录射线位置
        ray_x(end + 1) = current_pos(1);
        ray_z(end + 1) = current_pos(2);

        % 终止条件
        if norm(current_pos - source) < tol
            break;
        end
    end
end

三维各向异性

function anisotropic_ray_tracing_example_3D()
    % 参数设置
    Lx = 10;  % x方向长度 (km)
    Ly = 10;  % y方向长度 (km)
    Lz = 10;  % z方向长度 (km)
    Nx = 51;  % x方向网格点数
    Ny = 51;  % y方向网格点数
    Nz = 51;  % z方向网格点数
    dx = Lx / (Nx - 1);  % x方向网格间距 (km)
    dy = Ly / (Ny - 1);  % y方向网格间距 (km)
    dz = Lz / (Nz - 1);  % z方向网格间距 (km)

    % 生成网格
    x = linspace(0, Lx, Nx);
    y = linspace(0, Ly, Ny);
    z = linspace(0, Lz, Nz);
    [X, Y, Z] = meshgrid(x, y, z);

    % 生成示例走时场（假设各向异性介质）
    epsilon = 0.4;
    delta = 0.1;
    v0 = 3;  % 波速 (km/s)
    T = generate_anisotropic_traveltime_3D(X, Y, Z, v0, epsilon, delta);

    % 炮点和多个检波点位置
    source = [0, 0, 0];   % 炮点位置 (x, y, z)
   receivers = [7, 7, 7; 3, 5, 5; 5, 8, 8; 9, 2, 2; 6, 9, 9; 8, 4, 4]; % 多个检波点位置 (x, y, z)

    % 进行射线追踪
    figure;
    hold on;

    % 绘制走时场切片
    xslice = [0]; % x方向的切片位置
    yslice = [Ly/2, 4*Ly/4]; % y方向的切片位置
    zslice = [Lz/2, 0]; % z方向的切片位置
    hslice = slice(X, Y, Z, T, xslice, yslice, zslice);
    set(hslice, 'EdgeColor', 'none', 'FaceAlpha', 0.5); % 设置透明度并移除网格

    % 添加比色卡和标题
    colorbar;
    colormap();
    c = colorbar;
    c.Label.String = '走时 (s)';

    plot3(source(1), source(2), source(3), 'go', 'MarkerFaceColor', 'g'); % 炮点
    for i = 1:size(receivers, 1)
        receiver = receivers(i, :);
        [ray_x, ray_y, ray_z] = anisotropic_ray_tracing_3D(T, source, receiver, v0, epsilon, delta, dx, dy, dz);
        plot3(ray_x, ray_y, ray_z, 'r', 'LineWidth', 2); % 绘制射线路径
        plot3(receiver(1), receiver(2), receiver(3), 'bo', 'MarkerFaceColor', 'b'); % 检波点
    end

    xlabel('X (km)');
    ylabel('Y (km)');
    zlabel('Z (km)');
    title('各向异性介质中的3D射线追踪');
    axis([0 Lx 0 Ly 0 Lz]);
    grid on;
    hold off;
end

function T = generate_anisotropic_traveltime_3D(X, Y, Z, v0, epsilon, delta)
    % 生成各向异性介质中的走时场
    [Ny, Nx, Nz] = size(X);
    T = zeros(Ny, Nx, Nz);
    for i = 1:Nx
        for j = 1:Ny
            for k = 1:Nz
                theta = atan2(sqrt(X(j, i, k)^2 + Y(j, i, k)^2), Z(j, i, k));
                phi = atan2(Y(j, i, k), X(j, i, k));
                v = v0 * sqrt(1 + 2 * epsilon * sin(theta)^2 * cos(phi)^2 + 2 * delta * cos(theta)^2);
                T(j, i, k) = sqrt(X(j, i, k)^2 + Y(j, i, k)^2 + Z(j, i, k)^2) / v;
            end
        end
    end
end

function [ray_x, ray_y, ray_z] = anisotropic_ray_tracing_3D(T, source, receiver, v0, epsilon, delta, dx, dy, dz)
    % 初始化
    current_pos = receiver;
    ray_x = current_pos(1);
    ray_y = current_pos(2);
    ray_z = current_pos(3);

    max_iter = 1000;
    tol = min([dx, dy, dz]);

    for iter = 1:max_iter
        % 获取当前点的索引
        ix = round(current_pos(1) / dx) + 1;
        iy = round(current_pos(2) / dy) + 1;
        iz = round(current_pos(3) / dz) + 1;

        % 边界检查
        if ix < 1 || ix > size(T, 2) || iy < 1 || iy > size(T, 1) || iz < 1 || iz > size(T, 3)
            break;
        end

        % 计算当前点的速度和群速度方向
        theta = atan2(sqrt(current_pos(1)^2 + current_pos(2)^2), current_pos(3));
        phi = atan2(current_pos(2), current_pos(1));
        v = v0 * sqrt(1 + 2 * epsilon * sin(theta)^2 * cos(phi)^2 + 2 * delta * cos(theta)^2);
        vg_x = v * sin(theta) * cos(phi);
        vg_y = v * sin(theta) * sin(phi);
        vg_z = v * cos(theta);

        % 更新位置
        step_length = tol;
        step_direction = [vg_x, vg_y, vg_z] / norm([vg_x, vg_y, vg_z]);
        current_pos = current_pos - step_length * step_direction;

        % 记录射线位置
        ray_x(end + 1) = current_pos(1);
        ray_y(end + 1) = current_pos(2);
        ray_z(end + 1) = current_pos(3);

        % 终止条件
        if norm(current_pos - source) < tol
            break;
        end
    end
end