Deepwave 声波正演和弹性波正演

Deepwave

Deepwave 调用 scalar 方法实现声波和弹性波正演。

######## 声波正演 ######################
import torch
import numpy as np
import deepwave
from deepwave import scalardevice = torch.device('cuda' if torch.cuda.is_available()else 'cpu')## Set observation system paramaters
ny = 2301	# Y方向的网格数
nx = 751	# X方向的网格数
dx = 4.0	# 网格间距（m）
v = torch.from_file('marmousi_vp.bin',size=ny*nx).reshape(ny, nx).to(device)	# (2301, 751)	# 速度模型# source parameters
n_shots = 115  # 总共的地震源数
n_sources_per_shot = 1	# 每个地震记录中的源个数
d_source = 20		# 源的间距，20 * 4 = 80 m
first_source = 10	# 第一个源的位置: 10 m
source_depth = 2	# 源的深度: 2 m# receiver parameters
n_receivers_per_shot = 384	# 每个记录中有多少个接收器
d_receiver = 6				# 接收器之间的间距
first_receiver = 0			# 第一个接收器的位置
receiver_depth = 2			# 接收器的深度freq = 25		# 震源频率 Hz
nt = 750		# 总时间步数
dt = 0.004		# 采样间隔 s
peak_time = 1.5 / freq	# 震源持续波长# source locations 
# 三维: 一维->一个地震数据包含几炮; 二维->一炮包含几个震源; 三维->一个震源的位置 (x, y)
source_locations = torch.zeros(n_shots, n_sources_per_shot, 2, dtype=torch.long, device=device)
source_locations[..., 1] = source_depth		
source_locations[:, 0, 0] = (torch.arange(n_shots) * d_source + first_source)# receiver locations
# 三维: 一维->一个地震数据包含几炮; 二维->一炮包含几个接收器; 三维->一个接收器的位置 (x, y)
receiver_locations = torch.zeros(n_shots, n_receivers_per_shot, 2, dtype=torch.long, device=device)
receiver_locations[..., 1] = receiver_depth
receiver_locations[:, :, 0] = ((torch.arange(n_receivers_per_shot) * d_receiver + first_receiver).repeat(n_shots, 1))# source_amplitudes
source_amplitudes = (deepwave.wavelets.ricker(freq, nt, dt, peak_time)	# (震源频率, 接收器采样次数，接收器采样间隔，雷克子波波长).repeat(n_shots, n_sources_per_shot, 1).to(device)
)	# (20, 1, 750)## 调用 scalar 进行正演
## 声波
## 返回7个对象，前2个为波长快照的最后时刻，中间四个为对应的PML边界图，最后一个是接收到的地震记录。
out = scalar(v, dx, dt, source_amplitudes=source_amplitudes,source_locations=source_locations,receiver_locations=receiver_locations,accuracy=8,	# 空间有限差分的阶数 (2, 4, 6, 8)pml_freq=freq)	# 指定 PML 的主频率, 有助于最大限度地减少边缘反射.
receiver_amplitudes = out[-1] 	# 接收到的地震记录     

######## 弹性波正演 ######################
import torch
import numpy as np
import deepwave
from deepwave import elasticdevice = torch.device('cuda' if torch.cuda.is_available()else 'cpu')## Set observation system paramaters
ny = 2301	# Y方向的网格数
nx = 751	# X方向的网格数
dx = 4.0	# 网格间距（m）## create vp, vs, rho
vp_background = torch.ones(ny, nx, device=device) * 1500
vs_background = torch.ones(ny, nx, device=device) * 1000
rho_background = torch.ones(ny, nx, device=device) * 2200vp_true = vp_background.clone()
vp_true[10:20, 30:40] = 1600
vs_true = vs_background.clone()
vs_true[10:20, 45:55] = 1100
rho_true = rho_background.clone()
rho_true[10:20, 60:70] = 2300# source parameters
n_shots = 8 # 总共的地震源数
n_sources_per_shot = 1	# 每个地震记录中的源个数
d_source = 20 	# 源的间距，20 * 4 = 80 m
first_source = 10	# 第一个源的位置: 10 m
source_depth = 2	# 源的深度: 2 m# receiver parameters
n_receivers_per_shot = 384	# 每个记录中有多少个接收器
d_receiver = 6				# 接收器之间的间距
first_receiver = 0			# 第一个接收器的位置
receiver_depth = 2			# 接收器的深度freq = 25		# 震源频率 Hz
nt = 750		# 总时间步数
dt = 0.004		# 采样间隔 s
peak_time = 1.5 / freq	# 震源持续波长# source locations 
# 三维: 一维->一个地震数据包含几炮; 二维->一炮包含几个震源; 三维->一个震源的位置 (x, y)
source_locations = torch.zeros(n_shots, n_sources_per_shot, 2, dtype=torch.long, device=device)
source_locations[..., 1] = source_depth		
source_locations[:, 0, 0] = (torch.arange(n_shots) * d_source + first_source)# receiver locations
# 三维: 一维->一个地震数据包含几炮; 二维->一炮包含几个接收器; 三维->一个接收器的位置 (x, y)
receiver_locations = torch.zeros(n_shots, n_receivers_per_shot, 2, dtype=torch.long, device=device)
receiver_locations[..., 1] = receiver_depth
receiver_locations[:, :, 0] = ((torch.arange(n_receivers_per_shot) * d_receiver + first_receiver).repeat(n_shots, 1))# source_amplitudes
source_amplitudes = (deepwave.wavelets.ricker(freq, nt, dt, peak_time)	# (震源频率, 接收器采样次数，接收器采样间隔，雷克子波波长).repeat(n_shots, n_sources_per_shot, 1).to(device)
)	# (20, 1, 750)## 调用 elastic 进行弹性波正演
out = elastic(*deepwave.common.vpvsrho_to_lambmubuoyancy(vp_true, vs_true,rho_true),dx, dt,source_amplitudes_y=source_amplitudes,source_locations_y=source_locations,receiver_locations_y=receiver_locations,pml_freq=freq,
)[-2] 	# 接收到的地震记录  (8, 384, 750)