有限差分法在matlab中绘制0-100km内地震传播动图

有限差分法是一种常用的地震波传播数值模拟方法，可以用来模拟地震波在地球内部的传播情况。下面介绍如何在 Matlab 中使用有限差分法绘制地震波传播动图。

1. 准备数据

首先需要准备地震波传播的速度模型数据，可以使用公开的数据集，例如 IRIS 的 Global Velocity Model (GVM) 数据集，或者自行采集或创建速度模型数据。

2. 编写有限差分法程序

根据地震波传播的数学模型，编写有限差分法程序。有限差分法的主要思路是将地震波传播过程离散化，将地球内部划分为网格点，在每个网格点处求解波场的变化，从而模拟地震波的传播过程。

一般而言，有限差分法程序的主要步骤包括：

- 定义网格点和时间步长
- 初始化速度模型和波场
- 计算波场的一阶和二阶时间导数
- 更新波场
- 边界处理

这里给出一个简单的有限差分法程序示例，仅供参考：

% 定义模拟参数
nx = 200; % 网格点数
ny = 200;
nt = 1000; % 时间步数
dx = 10e3; % 网格间距
dy = 10e3;
dt = 0.005; % 时间步长

% 初始化速度模型和波场
vp = ones(nx, ny) * 2.0; % P波速度模型
vs = ones(nx, ny) * 1.0; % S波速度模型
rho = ones(nx, ny) * 2.0; % 密度模型
p = zeros(nx, ny); % P波波场
s = zeros(nx, ny); % S波波场
pold = zeros(nx, ny); % P波波场上一时刻的值
sold = zeros(nx, ny); % S波波场上一时刻的值

% 计算波场的一阶和二阶时间导数
for it = 1:nt
    % 计算P波波场的一阶和二阶时间导数
    ptt = (vp.^2 .* (p(3:end, 2:end-1) - 2*p(2:end-1, 2:end-1) + p(1:end-2, 2:end-1)) + ...
           vp.^2 .* (p(2:end-1, 3:end) - 2*p(2:end-1, 2:end-1) + p(2:end-1, 1:end-2))) / dx^2;
    ptx = (vp.^2 .* (p(2:end-1, 3:end) - p(2:end-1, 2:end-1)) - ...
           vp.^2 .* (p(2:end-1, 2:end-1) - p(2:end-1, 1:end-2))) / dx;
    pty = (vp.^2 .* (p(3:end, 2:end-1) - p(2:end-1, 2:end-1)) - ...
           vp.^2 .* (p(2:end-1, 2:end-1) - p(1:end-2, 2:end-1))) / dy;
    ptt2 = (vp.^2 .* (p(3:end, 2:end-1) + p(1:end-2, 2:end-1) - 2*p(2:end-1, 2:end-1)) + ...
            vp.^2 .* (p(2:end-1, 3:end) + p(2:end-1, 1:end-2) - 2*p(2:end-1, 2:end-1))) / dx^2;
    pxy = (vp.^2 .* (p(3:end, 3:end) - p(3:end, 1:end-2) - p(1:end-2, 3:end) + p(1:end-2, 1:end-2))) / (4*dx*dy);

    % 计算S波波场的一阶和二阶时间导数
    stt = (vs.^2 .* (s(3:end, 2:end-1) - 2*s(2:end-1, 2:end-1) + s(1:end-2, 2:end-1)) + ...
           vs.^2 .* (s(2:end-1, 3:end) - 2*s(2:end-1, 2:end-1) + s(2:end-1, 1:end-2))) / dx^2;
    stx = (vs.^2 .* (s(2:end-1, 3:end) - s(2:end-1, 2:end-1)) - ...
           vs.^2 .* (s(2:end-1, 2:end-1) - s(2:end-1, 1:end-2))) / dx;
    sty = (vs.^2 .* (s(3:end, 2:end-1) - s(2:end-1, 2:end-1)) - ...
           vs.^2 .* (s(2:end-1, 2:end-1) - s(1:end-2, 2:end-1))) / dy;
    stt2 = (vs.^2 .* (s(3:end, 2:end-1) + s(1:end-2, 2:end-1) - 2*s(2:end-1, 2:end-1)) + ...
            vs.^2 .* (s(2:end-1, 3:end) + s(2:end-1, 1:end-2) - 2*s(2:end-1, 2:end-1))) / dx^2;
    sxy = (vs.^2 .* (s(3:end, 3:end) - s(3:end, 1:end-2) - s(1:end-2, 3:end) + s(1:end-2, 1:end-2))) / (4*dx*dy);

    % 更新波场
    pnew = 2*p - pold + dt^2 * (ptt + ptt2 .* rho + ptx + pty + pxy .* rho);
    snew = 2*s - sold + dt^2 * (stt + stt2 .* rho + stx + sty + sxy .* rho);

    % 边界处理（此处采用完全吸收边界条件）
    pnew(1,:) = 0; pnew(end,:) = 0; pnew(:,1) = 0; pnew(:,end) = 0;
    snew(1,:) = 0; snew(end,:) = 0; snew(:,1) = 0; snew(:,end) = 0;

    % 更新波场变量
    pold = p; p = pnew; sold = s; s = snew;

    % 绘制波场动图
    imagesc(p); colorbar; caxis([-0.1 0.1]); drawnow;
end

3. 绘制动图

通过上述有限差分法程序，可以得到地震波在每个时间步长时的波场变化情况，从而可以绘制地震波传播的动图。可以使用 Matlab 的 `imagesc` 函数绘制波场图像，然后使用 `drawnow` 函数实时绘制波场动图。

下面给出一个简单的绘制动图的示例代码：

% 绘制波场动图
for it = 1:nt
    imagesc(p); colorbar; caxis([-0.1 0.1]); drawnow;
end

在绘制动图时，可以调整颜色条的范围、时间步长等参数，以适应不同的需求。