快速行进法（FMM）求解程函方程的基本原理，以及VTI介质和TTI介质中的走时计算

function fast_marching_method
    % 参数设置
    Lx = 2000;  % x方向长度 (m)
    Lz = 1000;  % z方向长度 (m)
    Nx = 101;   % x方向网格点数
    Nz = 51;    % z方向网格点数
    dx = Lx / (Nx - 1);  % x方向网格间距 (m)
    dz = Lz / (Nz - 1);  % z方向网格间距 (m)

    x = linspace(0, Lx, Nx);
    z = linspace(0, Lz, Nz);
    [X, Z] = meshgrid(x, z);

    % 定义双层速度模型
    v1 = 1500;  % 上层速度 (m/s)
    v2 = 3000;  % 下层速度 (m/s)
    depth = 500;  % 两层分界深度 (m)

    v = v1 * ones(Nz, Nx);
    v(Z > depth) = v2;

    % 初始化时间场
    T = inf(Nz, Nx);
    T(1, floor(Nx/2)) = 0;  % 在地表中心放置震源

    % 快速行进法求解程函方程
    T = fast_marching(T, v, dx, dz);

    % 绘制时间场
    figure;
    imagesc(x, z, T);
    colorbar;
    title('时间场');
    xlabel('距离 (m)');
    ylabel('深度 (m)');
end

function T = fast_marching(T, v, dx, dz)
    [Nz, Nx] = size(T);
    known = false(Nz, Nx);
    trial = false(Nz, Nx);

    % 初始化已知和候选点
    known(T == 0) = true;
    [zi, xi] = find(T == 0);
    for i = 1:length(zi)
        neighbors = [zi(i)-1, xi(i); zi(i)+1, xi(i); zi(i), xi(i)-1; zi(i), xi(i)+1];
        for k = 1:size(neighbors, 1)
            nz = neighbors(k, 1);
            nx = neighbors(k, 2);
            if nz > 0 && nz <= Nz && nx > 0 && nx <= Nx && ~known(nz, nx)
                T(nz, nx) = min(T(nz, nx), solve_quadratic(T, v, nz, nx, dx, dz));
                trial(nz, nx) = true;
            end
        end
    end

    % 主循环
    while any(trial(:))
        [~, idx] = min(T(trial));
        [zi, xi] = find(trial);
        zi = zi(idx);
        xi = xi(idx);

        known(zi, xi) = true;
        trial(zi, xi) = false;

        neighbors = [zi-1, xi; zi+1, xi; zi, xi-1; zi, xi+1];
        for k = 1:size(neighbors, 1)
            nz = neighbors(k, 1);
            nx = neighbors(k, 2);
            if nz > 0 && nz <= Nz && nx > 0 && nx <= Nx && ~known(nz, nx)
                T(nz, nx) = min(T(nz, nx), solve_quadratic(T, v, nz, nx, dx, dz));
                trial(nz, nx) = true;
            end
        end
    end
end

function Tn = solve_quadratic(T, v, nz, nx, dx, dz)
    % 使用有限差分法求解离散程函方程
    T_left = inf;
    T_right = inf;
    T_up = inf;
    T_down = inf;

    if nz > 1
        T_up = T(nz-1, nx);
    end
    if nz < size(T, 1)
        T_down = T(nz+1, nx);
    end
    if nx > 1
        T_left = T(nz, nx-1);
    end
    if nx < size(T, 2)
        T_right = T(nz, nx+1);
    end

    T_x = min(T_left, T_right);
    T_z = min(T_up, T_down);

    if isinf(T_x)
        Tn = T_z + dz / v(nz, nx);
    elseif isinf(T_z)
        Tn = T_x + dx / v(nz, nx);
    else
        a = 1 / dx^2 + 1 / dz^2;
        b = -2 * (T_x / dx^2 + T_z / dz^2);
        c = T_x^2 / dx^2 + T_z^2 / dz^2 - 1 / v(nz, nx)^2;
        delta = b^2 - 4 * a * c;
        if delta >= 0
            Tn = (-b + sqrt(delta)) / (2 * a);
        else
            Tn = inf;
        end
    end
end

function fast_marching_vtiqq
    % 参数设置
    Lx = 2000;  % x方向长度 (m)
    Lz = 1000;  % z方向长度 (m)
    Nx = 101;   % x方向网格点数
    Nz = 51;    % z方向网格点数
    dx = Lx / (Nx - 1);  % x方向网格间距 (m)
    dz = Lz / (Nz - 1);  % z方向网格间距 (m)

    x = linspace(0, Lx, Nx);
    z = linspace(0, Lz, Nz);
    [X, Z] = meshgrid(x, z);

    % 定义VTI介质参数
    vp1 = 1500;  % 上层P波速度 (m/s)
    vp2 = 3000;  % 下层P波速度 (m/s)
    epsilon1 = 0.1;  % 上层Thomsen参数
    epsilon2 = 0.2;  % 下层Thomsen参数
    delta1 = 0.05;   % 上层Thomsen参数
    delta2 = 0.1;    % 下层Thomsen参数
    depth = 500;     % 两层分界深度 (m)

    vp = vp1 * ones(Nz, Nx);
    epsilon = epsilon1 * ones(Nz, Nx);
    delta = delta1 * ones(Nz, Nx);
    vp(Z > depth) = vp2;
    epsilon(Z > depth) = epsilon2;
    delta(Z > depth) = delta2;

    % 初始化时间场
    T = inf(Nz, Nx);
    T(1, floor(Nx/2)) = 0;  % 在地表中心放置震源

    % 快速行进法求解程函方程
    T = fast_marching_vti(T, vp, epsilon, delta, dx, dz);

    % 绘制时间场
    figure;
    imagesc(x, z, T);
    colorbar;
    title('时间场');
    xlabel('距离 (m)');
    ylabel('深度 (m)');
end

function T = fast_marching_vti(T, vp, epsilon, delta, dx, dz)
    [Nz, Nx] = size(T);
    known = false(Nz, Nx);
    trial = false(Nz, Nx);

    % 初始化已知和候选点
    known(T == 0) = true;
    [zi, xi] = find(T == 0);
    for i = 1:length(zi)
        neighbors = [zi(i)-1, xi(i); zi(i)+1, xi(i); zi(i), xi(i)-1; zi(i), xi(i)+1];
        for k = 1:size(neighbors, 1)
            nz = neighbors(k, 1);
            nx = neighbors(k, 2);
            if nz > 0 && nz <= Nz && nx > 0 && nx <= Nx && ~known(nz, nx)
                T(nz, nx) = min(T(nz, nx), solve_vti(T, vp, epsilon, delta, nz, nx, dx, dz));
                trial(nz, nx) = true;
            end
        end
    end

    % 主循环
    while any(trial(:))
        [~, idx] = min(T(trial));
        [zi, xi] = find(trial);
        zi = zi(idx);
        xi = xi(idx);

        known(zi, xi) = true;
        trial(zi, xi) = false;

        neighbors = [zi-1, xi; zi+1, xi; zi, xi-1; zi, xi+1];
        for k = 1:size(neighbors, 1)
            nz = neighbors(k, 1);
            nx = neighbors(k, 2);
            if nz > 0 && nz <= Nz && nx > 0 && nx <= Nx && ~known(nz, nx)
                T(nz, nx) = min(T(nz, nx), solve_vti(T, vp, epsilon, delta, nz, nx, dx, dz));
                trial(nz, nx) = true;
            end
        end
    end
end

function Tn = solve_vti(T, vp, epsilon, delta, nz, nx, dx, dz)
    % 使用有限差分法求解离散程函方程，考虑VTI介质的影响
    T_left = inf;
    T_right = inf;
    T_up = inf;
    T_down = inf;

    if nz > 1
        T_up = T(nz-1, nx);
    end
    if nz < size(T, 1)
        T_down = T(nz+1, nx);
    end
    if nx > 1
        T_left = T(nz, nx-1);
    end
    if nx < size(T, 2)
        T_right = T(nz, nx+1);
    end

    T_x = min(T_left, T_right);
    T_z = min(T_up, T_down);

    if isinf(T_x)
        Tn = T_z + dz / (vp(nz, nx) * sqrt(1 + 2 * epsilon(nz, nx)));
    elseif isinf(T_z)
        Tn = T_x + dx / (vp(nz, nx) * sqrt(1 + 2 * epsilon(nz, nx)));
    else
        a = 1 / dx^2 + 1 / dz^2;
        b = -2 * (T_x / dx^2 + T_z / dz^2);
        c = T_x^2 / dx^2 + T_z^2 / dz^2 - (1 / vp(nz, nx)^2) * (1 + 2 * epsilon(nz, nx));
        delta = b^2 - 4 * a * c;
        if delta >= 0
            Tn = (-b + sqrt(delta)) / (2 * a);
        else
            Tn = inf;
        end
    end
end

function fast_marching_tti1
    % 参数设置
    Lx = 2000;  % x方向长度 (m)
    Lz = 1000;  % z方向长度 (m)
    Nx = 101;   % x方向网格点数
    Nz = 51;    % z方向网格点数
    dx = Lx / (Nx - 1);  % x方向网格间距 (m)
    dz = Lz / (Nz - 1);  % z方向网格间距 (m)

    x = linspace(0, Lx, Nx);
    z = linspace(0, Lz, Nz);
    [X, Z] = meshgrid(x, z);

    % 定义TTI介质参数
    vp1 = 1500;  % 上层P波速度 (m/s)
    vp2 = 3000;  % 下层P波速度 (m/s)
    epsilon1 = 0.1;  % 上层Thomsen参数
    epsilon2 = 0.2;  % 下层Thomsen参数
    delta1 = 0.05;   % 上层Thomsen参数
    delta2 = 0.1;    % 下层Thomsen参数
    theta1 = 10;     % 上层倾斜角度 (度)
    theta2 = 20;     % 下层倾斜角度 (度)
    depth = 500;     % 两层分界深度 (m)

    vp = vp1 * ones(Nz, Nx);
    epsilon = epsilon1 * ones(Nz, Nx);
    delta = delta1 * ones(Nz, Nx);
    theta = theta1 * ones(Nz, Nx);
    vp(Z > depth) = vp2;
    epsilon(Z > depth) = epsilon2;
    delta(Z > depth) = delta2;
    theta(Z > depth) = theta2;

    % 将倾斜角度转换为弧度
    theta = deg2rad(theta);

    % 初始化时间场
    T = inf(Nz, Nx);
    T(1, floor(Nx/2)) = 0;  % 在地表中心放置震源

    % 快速行进法求解程函方程
    T = fast_marching_tti(T, vp, epsilon, delta, theta, dx, dz);

    % 绘制时间场
    figure;
    imagesc(x, z, T);
    colorbar;
    title('时间场');
    xlabel('距离 (m)');
    ylabel('深度 (m)');
end

function T = fast_marching_tti(T, vp, epsilon, delta, theta, dx, dz)
    [Nz, Nx] = size(T);
    known = false(Nz, Nx);
    trial = false(Nz, Nx);

    % 初始化已知和候选点
    known(T == 0) = true;
    [zi, xi] = find(T == 0);
    for i = 1:length(zi)
        neighbors = [zi(i)-1, xi(i); zi(i)+1, xi(i); zi(i), xi(i)-1; zi(i), xi(i)+1];
        for k = 1:size(neighbors, 1)
            nz = neighbors(k, 1);
            nx = neighbors(k, 2);
            if nz > 0 && nz <= Nz && nx > 0 && nx <= Nx && ~known(nz, nx)
                T(nz, nx) = min(T(nz, nx), solve_tti(T, vp, epsilon, delta, theta, nz, nx, dx, dz));
                trial(nz, nx) = true;
            end
        end
    end

    % 主循环
    while any(trial(:))
        [~, idx] = min(T(trial));
        [zi, xi] = find(trial);
        zi = zi(idx);
        xi = xi(idx);

        known(zi, xi) = true;
        trial(zi, xi) = false;

        neighbors = [zi-1, xi; zi+1, xi; zi, xi-1; zi, xi+1];
        for k = 1:size(neighbors, 1)
            nz = neighbors(k, 1);
            nx = neighbors(k, 2);
            if nz > 0 && nz <= Nz && nx > 0 && nx <= Nx && ~known(nz, nx)
                T(nz, nx) = min(T(nz, nx), solve_tti(T, vp, epsilon, delta, theta, nz, nx, dx, dz));
                trial(nz, nx) = true;
            end
        end
    end
end

function Tn = solve_tti(T, vp, epsilon, delta, theta, nz, nx, dx, dz)
    % 使用有限差分法求解离散程函方程，考虑TTI介质的影响
    T_left = inf;
    T_right = inf;
    T_up = inf;
    T_down = inf;

    if nz > 1
        T_up = T(nz-1, nx);
    end
    if nz < size(T, 1)
        T_down = T(nz+1, nx);
    end
    if nx > 1
        T_left = T(nz, nx-1);
    end
    if nx < size(T, 2)
        T_right = T(nz, nx+1);
    end

    T_x = min(T_left, T_right);
    T_z = min(T_up, T_down);

    if isinf(T_x)
        Tn = T_z + dz / (vp(nz, nx) * sqrt(1 + 2 * epsilon(nz, nx) * sin(theta(nz, nx))^2));
    elseif isinf(T_z)
        Tn = T_x + dx / (vp(nz, nx) * sqrt(1 + 2 * epsilon(nz, nx) * sin(theta(nz, nx))^2));
    else
        a = 1 / dx^2 + 1 / dz^2;
        b = -2 * (T_x / dx^2 + T_z / dz^2);
        c = T_x^2 / dx^2 + T_z^2 / dz^2 - (1 / vp(nz, nx)^2) * (1 + 2 * epsilon(nz, nx) * sin(theta(nz, nx))^2);
        delta = b^2 - 4 * a * c;
        if delta >= 0
            Tn = (-b + sqrt(delta)) / (2 * a);
        else
            Tn = inf;
        end
    end
end