VTI介质地震波传播模拟（PML边界条件）
用于模拟VTI（垂直横向各向同性）介质中的地震波传播，使用PML（完美匹配层）边界条件，并生成波场快照。
% VTI介质地震波传播模拟 - PML边界条件
% 使用各向异性弹性波方程和有限差分法
clear all; close all; clc;
%% 参数设置
% 模型参数
nx = 200; % x方向网格点数
nz = 200; % z方向网格点数
dx = 10; % x方向网格间距(m)
dz = 10; % z方向网格间距(m)
dt = 0.5e-3; % 时间步长(s)
nt = 1000; % 时间步数
f0 = 25; % 震源主频(Hz)
% VTI介质参数 (Thomsen参数)
rho = 2500; % 密度(kg/m^3)
vp0 = 3000; % 垂直P波速度(m/s)
vs0 = 1500; % 垂直S波速度(m/s)
epsilon = 0.2; % Thomsen参数ε
delta = 0.1; % Thomsen参数δ
gamma = 0.15; % Thomsen参数γ
% 计算VTI刚度矩阵元素
C33 = rho * vp0^2;
C44 = rho * vs0^2;
C11 = C33 * (1 + 2*epsilon);
C66 = rho * vs0^2 * (1 + 2*gamma);
C13 = sqrt(C33^2 + C11^2 - 2*C33*C11 - 4*C44^2)/2 - C44;
% PML参数
pml_width = 20; % PML层宽度(网格点数)
R = 1e-6; % 理论反射系数
damp_max = 3; % 最大阻尼系数
% 震源位置
src_x = nx/2;
src_z = nz/2;
% 接收点位置
rec_x = src_x + 50;
rec_z = src_z;
%% 初始化模型
% 创建空间坐标
x = (0:nx-1)*dx;
z = (0:nz-1)*dz;
% 初始化场变量
vx = zeros(nx, nz); % x方向速度分量
vz = zeros(nx, nz); % z方向速度分量
sxx = zeros(nx, nz); % xx应力分量
szz = zeros(nx, nz); % zz应力分量
sxz = zeros(nx, nz); % xz剪切应力分量
% 初始化PML参数
[pml_alpha_x, pml_alpha_z] = setup_pml(nx, nz, pml_width, damp_max);
%% 震源函数 - Ricker子波
t = (0:nt-1)*dt;
src_func = (1 - 2*(pi*f0*t).^2) .* exp(-(pi*f0*t).^2);
%% 主循环 - 时间推进
for it = 1:nt
% 添加震源 (力源在中心点)
if it <= length(src_func)
sxx(src_x, src_z) = sxx(src_x, src_z) + src_func(it) * dt;
end
% 更新速度分量 (x方向)
for i = 2:nx-1
for j = 2:nz-1
% 计算空间导数
dsxx_dx = (sxx(i+1,j) - sxx(i-1,j))/(2*dx);
dsxz_dz = (sxz(i,j+1) - sxz(i,j-1))/(2*dz);
% 更新速度 (考虑PML阻尼)
vx(i,j) = vx(i,j) + dt/rho * (dsxx_dx + dsxz_dz) - pml_alpha_x(i,j)*vx(i,j)*dt;
end
end
% 更新速度分量 (z方向)
for i = 2:nx-1
for j = 2:nz-1
% 计算空间导数
dsxz_dx = (sxz(i+1,j) - sxz(i-1,j))/(2*dx);
dszz_dz = (szz(i,j+1) - szz(i,j-1))/(2*dz);
% 更新速度 (考虑PML阻尼)
vz(i,j) = vz(i,j) + dt/rho * (dsxz_dx + dszz_dz) - pml_alpha_z(i,j)*vz(i,j)*dt;
end
end
% 更新应力分量 (xx)
for i = 2:nx-1
for j = 2:nz-1
% 计算空间导数
dvx_dx = (vx(i+1,j) - vx(i-1,j))/(2*dx);
dvz_dz = (vz(i,j+1) - vz(i,j-1))/(2*dz);
% 更新应力 (VTI介质)
sxx(i,j) = sxx(i,j) + dt * (C11*dvx_dx + C13*dvz_dz);
end
end
% 更新应力分量 (zz)
for i = 2:nx-1
for j = 2:nz-1
% 计算空间导数
dvx_dx = (vx(i+1,j) - vx(i-1,j))/(2*dx);
dvz_dz = (vz(i,j+1) - vz(i,j-1))/(2*dz);
% 更新应力 (VTI介质)
szz(i,j) = szz(i,j) + dt * (C13*dvx_dx + C33*dvz_dz);
end
end
% 更新应力分量 (xz)
for i = 2:nx-1
for j = 2:nz-1
% 计算空间导数
dvx_dz = (vx(i,j+1) - vx(i,j-1))/(2*dz);
dvz_dx = (vz(i+1,j) - vz(i-1,j))/(2*dx);
% 更新应力 (VTI介质)
sxz(i,j) = sxz(i,j) + dt * C44 * (dvx_dz + dvz_dx);
end
end
% 边界条件 - 自由表面 (顶部边界)
j = 1; % 顶部边界
szz(:,j) = 0;
sxz(:,j) = 0;
% 保存波场快照 (每隔一定步数)
if mod(it, 50) == 0
snap_vx(:,:,it/50) = vx;
snap_vz(:,:,it/50) = vz;
snap_sxx(:,:,it/50) = sxx;
snap_szz(:,:,it/50) = szz;
snap_sxz(:,:,it/50) = sxz;
time_snap(it/50) = it*dt;
end
end
%% 波场快照可视化
plot_wavefield_snapshots(time_snap, snap_vx, snap_vz, snap_sxx, snap_szz, snap_sxz, nx, nz, dx, dz);
%% PML设置函数
function [pml_alpha_x, pml_alpha_z] = setup_pml(nx, nz, width, damp_max)
% 初始化PML阻尼系数
pml_alpha_x = zeros(nx, nz);
pml_alpha_z = zeros(nx, nz);
% x方向PML (左右边界)
for i = 1:width
% 左边界
coeff = damp_max * (1 - cos(pi*i/(2*width))) / 2;
pml_alpha_x(i, :) = coeff;
pml_alpha_x(nx-i+1, :) = coeff;
% z方向PML (上下边界)
pml_alpha_z(:, i) = coeff;
pml_alpha_z(:, nz-i+1) = coeff;
end
end
%% 波场快照可视化函数
function plot_wavefield_snapshots(time_snap, snap_vx, snap_vz, snap_sxx, snap_szz, snap_sxz, nx, nz, dx, dz)
% 创建网格
[X, Z] = meshgrid(0:dx:(nx-1)*dx, 0:dz:(nz-1)*dz);
% 选择要显示的快照
snapshot_indices = [2, 4, 6, 8, 10]; % 对应时间步200, 400, 600, 800, 1000
% 创建图形
figure('Position', [100, 100, 1200, 900]);
for i = 1:length(snapshot_indices)
idx = snapshot_indices(i);
t = time_snap(idx);
% 绘制垂直速度分量
subplot(5, 3, (i-1)*3+1);
imagesc(X(1,:), Z(:,1), snap_vz(:,:,idx)');
colormap(seismic); colorbar;
title(sprintf('Vz (t=%.1fs)', t));
xlabel('X (m)'); ylabel('Z (m)');
axis equal tight;
caxis([-1e-6, 1e-6]);
% 绘制水平速度分量
subplot(5, 3, (i-1)*3+2);
imagesc(X(1,:), Z(:,1), snap_vx(:,:,idx)');
colormap(seismic); colorbar;
title(sprintf('Vx (t=%.1fs)', t));
xlabel('X (m)'); ylabel('Z (m)');
axis equal tight;
caxis([-1e-6, 1e-6]);
% 绘制压力场 (P波)
pressure = (snap_sxx(:,:,idx) + snap_szz(:,:,idx))/2;
subplot(5, 3, (i-1)*3+3);
imagesc(X(1,:), Z(:,1), pressure');
colormap(seismic); colorbar;
title(sprintf('Pressure (t=%.1fs)', t));
xlabel('X (m)'); ylabel('Z (m)');
axis equal tight;
caxis([-1e7, 1e7]);
end
sgtitle('VTI介质地震波传播波场快照 (PML边界条件)');
% 单独绘制剪切波场
figure('Position', [100, 100, 800, 600]);
shear_stress = snap_sxz(:,:,end);
imagesc(X(1,:), Z(:,1), shear_stress');
colormap(seismic); colorbar;
title(sprintf('剪切应力 Sxz (t=%.1fs)', time_snap(end)));
xlabel('X (m)'); ylabel('Z (m)');
axis equal tight;
caxis([-1e7, 1e7]);
end
程序功能说明
1. 物理模型
- VTI介质：使用Thomsen参数(ε, δ, γ)描述垂直横向各向同性
- 弹性波方程：实现各向异性介质中的速度-应力形式方程
- PML边界：完美匹配层吸收边界条件，有效减少人工反射
2. 数值方法
- 有限差分法：二阶中心差分格式
- 交错网格：速度和应力变量错开存储
- 时间推进：显式蛙跳格式
3. 关键组件
- 震源函数：Ricker子波作为地震激发源
- 自由表面：顶部边界设置为自由表面条件
- 波场快照：定期保存并可视化波场演化
4. 可视化功能
- 多时刻波场快照显示
- 速度分量(Vx, Vz)和应力分量(Sxx, Szz, Sxz)可视化
- 压力场(P波)和剪切波场分离显示
- 使用seismic色标增强对比度
运行结果
程序运行后会生成两组图形：
-
主波场快照图：5行3列共15个子图，显示5个不同时间点的：
- 垂直速度分量(Vz)
- 水平速度分量(Vx)
- 压力场(平均应力)
-
剪切波场图：显示最终时刻的剪切应力分量(Sxz)
参数调整建议
-
模型尺寸：
- 增加
nx和nz提高空间分辨率 - 减小
dx和dz提高精度（需相应减小dt）
- 增加
-
震源参数：
- 修改
f0改变主频（低频穿透深，高频分辨率高） - 调整震源位置
src_x,src_z
- 修改
-
介质参数：
- 修改Thomsen参数
epsilon,delta,gamma模拟不同VTI特性 - 调整密度
rho和速度vp0,vs0
- 修改Thomsen参数
-
PML设置：
- 增加
pml_width提高吸收效果 - 调整
damp_max控制阻尼强度
- 增加
参考代码 通过Matlab编程，给定初始条件及PML边界条件，模拟VTI介质中地震波传播的波场快照 www.youwenfan.com/contentcss/79964.html
扩展功能
如需进一步扩展，可考虑：
- 添加各向异性参数渐变层
- 实现更复杂的震源机制
- 添加接收器记录地震道
- 实现三维模拟
- 添加地层界面反射
