含震源二维波动方程求解

已于 2024-06-09 13:48:42 修改
阅读量577 收藏 17
点赞数 6
文章标签： 机器学习 算法 人工智能 python pytorch 深度学习
于 2024-06-09 13:46:03 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_58509609/article/details/139556527
版权

目录

前言

上次发了一个二维无震源的波动方程求解。为了模拟一下地震波场的传播，这次做一个含震源的均匀介质波动方程求解。
在这里插入图片描述

一、控制方程、初始条件及边界条件

我们考虑在x∈[0,8];y∈[0,8]的二维空间中，存在一个沿Z轴方向振动且波速为的波场。选取雷克子波函数作为空间中的震源，则其在空间中振动时间-位移关系满足下式：
在这里插入图片描述
设置边界条件和初始条件如下：
在这里插入图片描述
在这里插入图片描述

二、采样点设置

对于此二维波动方程，我们在x ∈ [0,8]; y ∈ [0,8]; t ∈ [0,2]构成的数值空间中，共计采样 10000 个采样点；在左右上下四个边界上各获取 5000 个采样点；在两个初始条件下获取 5000 个采样点。采样情况如下图所示。
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

三、训练

我们设置的网络主要有五个区块组成，每个区块里包含一个全连接层，以及对应的激活函数，全连接层包含 20 个节点，使用双曲正切函数为激活函数；训练期间向模型输入训练集和验证集，训练集为 3-3-2 节提到的采样数据，验证集为同样使用拉丁超立方采样算法进行采样的 1000 个波场内部采样数据；训练的学习率为0.001，预设置训练轮次为 10000 次。
在这里插入图片描述

四、结果展示

先看看结果吧，反正这个东西看了也就图一乐，初始条件的拟合结果还可以，波场的传播图像也符合常识。
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

五、代码详解

话不多上，上代码，反正参数咱也不咋会调，只能保证跑个结果出来。

1.预前准备

######################################################################
import torch
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import torch.nn as nn
from pyDOE import lhs
import imageio
import os
######################################################################
#自动更换CPU,GPU，并设置浮点型位数
default_device = "cuda" if torch.cuda.is_available() else "cpu"
dtype=torch.float32
torch.set_default_dtype(dtype)
#####################################################################
# 设置随机数种子
def setup_seed(seed):
    torch.manual_seed(seed)
    torch.cuda.manual_seed_all(seed)
    torch.backends.cudnn.deterministic = True

setup_seed(42)
######################################################################
#转换张量函数，存储损失的时候用到
def Tran_fun(tensor):
    a = []
    for i in tensor:
        i = i.to("cpu")
        i = i.detach().numpy()
        a.append(i)
    return(np.array(a))
######################################################################
#设置采样策略
def LHS_Sampling(N,be,en):
    return be + (en-be)*lhs(3, N) 
######################################################################
#定义雷克子波，在激发的情况下使用
def rec_wave(x,y,t,x0,y0,t0,f):
    a = (1-2*(torch.pi*f*(t-t0))**2)*torch.exp(-(np.pi*f*(t-t0))**2)
    b = torch.exp(-0.5*((x-x0)**2+(y-y0)**2)/(0.8)**2)
    return a*b
######################################################################
#求导函数
def gradients(u, x, order=1):
    u.requires_grad_(True)
    x.requires_grad_(True)
    if order == 1:
        return torch.autograd.grad(u, x, grad_outputs=torch.ones_like(u),
                                   create_graph=True,
                                   only_inputs=True, )[0]
    else:
        return gradients(gradients(u, x), x, order=order - 1)
######################################################################

2.网络编译

######################################################################
#网络结构的编译
class PINN(nn.Module):
    def __init__(self, NN): 
        super(PINN, self).__init__()

        self.input_layer = nn.Linear(3, NN)
        self.hidden_layer1 = nn.Linear(NN,NN) 
        self.hidden_layer2 = nn.Linear(NN, NN)  
        self.hidden_layer3 = nn.Linear(NN, NN)  
        self.hidden_layer4 = nn.Linear(NN, NN)  
        self.output_layer = nn.Linear(NN, 1)

    def forward(self, x): # 一种特殊的方法 __call__() 回调
        out = torch.tanh(self.input_layer(x))
        out = torch.tanh(self.hidden_layer1(out))
        out = torch.tanh(self.hidden_layer2(out))
        out = torch.tanh(self.hidden_layer3(out))
        out = torch.tanh(self.hidden_layer4(out))
        #out = torch.tanh(self.hidden_layer5(out))
        out_final = self.output_layer(out)
        return out_final

######################################################################
#损失函数，网络，优化器对象化
loss = torch.nn.MSELoss()
loss = loss.to(device = default_device)
#pinn = PINN(20)
#训练一部分保存的网络，在这里继续训练
pinn = torch.load("D:/trl/波场结构/add_resoure_wave_model2d_523.pth",map_location=torch.device('cpu'))
pinn = pinn.to(device = default_device)
opt = torch.optim.Adam(pinn.parameters(), lr=1e-3)
scheduler = torch.optim.lr_scheduler.StepLR(opt, step_size=5000, gamma=0.9)
######################################################################

3.损失函数设置

######################################################################
#pde损失函数设置
def interior_loss(Mat):
    # 内点
    Mat = torch.tensor(Mat)
    x = Mat[:, 0:1].float()
    y = Mat[:, 1:2].float()
    t = Mat[:, 2:3].float()
    x.requires_grad_(True)
    y.requires_grad_(True)
    t.requires_grad_(True)
    x = x.to(device=default_device)
    y = y.to(device = default_device)
    t = t.to(device = default_device)
    u = pinn(torch.cat([t, x, y], dim=1))
    cond1 = gradients(u, t, order=2)
    cond2 = (gradients(u, x, order=2)+gradients(u, y, order=2))*c**2 + rec_wave(x,y,t,x0=4,y0=4,t0=0,f=0.75)
    return loss(cond1,cond2)

def left_loss(Mat):
    # 左边界点
    Mat = torch.tensor(Mat)
    Mat[:,0:1] = torch.zeros_like(Mat[:,0:1])
    x = Mat[:, 0:1].float()
    y = Mat[:, 1:2].float()
    t = Mat[:, 2:3].float()
    zero = torch.zeros_like(t).float()
    x = x.to(device=default_device)
    y = y.to(device = default_device)
    t = t.to(device = default_device)
    zero = zero.to(device = default_device)
    u = pinn(torch.cat([t, x, y], dim=1))
    cond1 = u
    cond2 = zero
    return loss(cond1,cond2)

def right_loss(Mat):
    # 右边界点
    Mat = torch.tensor(Mat)
    Mat[:,0:1] = torch.ones_like(Mat[:,0:1])*8
    x = Mat[:, 0:1].float()
    y = Mat[:, 1:2].float()
    t = Mat[:, 2:3].float()
    zero = torch.zeros_like(t).float()
    x.requires_grad_(True)
    y.requires_grad_(True)
    t.requires_grad_(True)
    x = x.to(device=default_device)
    y = y.to(device = default_device)
    t = t.to(device = default_device)
    zero = zero.to(device = default_device)
    u = pinn(torch.cat([t, x, y], dim=1))
    zero = torch.zeros_like(u).requires_grad_(True)
    zero = zero.to(device = default_device)
    cond1 = u
    cond2 = zero
    return loss(cond1,cond2)

def down_loss(Mat):
    # 下边界点
    Mat = torch.tensor(Mat)
    Mat[:,1:2] = torch.zeros_like(Mat[:,1:2])
    x = Mat[:, 0:1].float()
    y = Mat[:, 1:2].float()
    t = Mat[:, 2:3].float()
    zero = torch.zeros_like(t).float()
    x.requires_grad_(True)
    y.requires_grad_(True)
    t.requires_grad_(True)
    zero.requires_grad_(True)
    x = x.to(device=default_device)
    y = y.to(device = default_device)
    t = t.to(device = default_device)
    zero = zero.to(device = default_device)
    u = pinn(torch.cat([t, x, y], dim=1))
    cond1 = u
    cond2 = zero
    return loss(cond1,cond2)

def up_loss(Mat):
    # 上边界点
    Mat = torch.tensor(Mat)
    Mat[:,1:2] = torch.ones_like(Mat[:,1:2])*8
    x = Mat[:, 0:1].float()
    y = Mat[:, 1:2].float()
    t = Mat[:, 2:3].float()
    zero = torch.zeros_like(t).float()
    x.requires_grad_(True)
    y.requires_grad_(True)
    t.requires_grad_(True)
    x = x.to(device=default_device)
    y = y.to(device = default_device)
    t = t.to(device = default_device)
    zero = zero.to(device = default_device)
    u = pinn(torch.cat([t, x, y], dim=1))
    cond1 = u
    cond2 = zero
    return loss(cond1,cond2)

def ini1_loss(Mat):
    # 初始条件1
    Mat = torch.tensor(Mat)
    Mat[:,2:3] = torch.zeros_like(Mat[:,2:3])
    x = Mat[:, 0:1].float()
    y = Mat[:, 1:2].float()
    t = Mat[:, 2:3].float()
    zero = torch.zeros_like(t).float()
    x.requires_grad_(True)
    y.requires_grad_(True)
    t.requires_grad_(True)
    zero.requires_grad_(True)
    x = x.to(device=default_device)
    y = y.to(device = default_device)
    t = t.to(device = default_device)
    u = pinn(torch.cat([t, x, y], dim=1))
    cond1 = u
    cond2 = rec_wave(x,y,t,x0=4,y0=4,t0=0,f=0.75)
    return loss(cond1,cond2)

def ini2_loss(Mat):
    # 初始条件2
    Mat = torch.tensor(Mat)
    Mat[:,2:3] = torch.zeros_like(Mat[:,2:3])
    x = Mat[:, 0:1].float()
    y = Mat[:, 1:2].float()
    t = Mat[:, 2:3].float()
    zero = torch.zeros_like(t).float()
    x.requires_grad_(True)
    y.requires_grad_(True)
    t.requires_grad_(True)
    zero.requires_grad_(True)
    x = x.to(device=default_device)
    y = y.to(device = default_device)
    t = t.to(device = default_device)
    u = pinn(torch.cat([t, x, y], dim=1))
    cond1 = gradients(u, t, order=1)
    cond2 = torch.exp(-0.5*((x-4)**2+(y-4)**2)/(0.8)**2)
    return loss(cond1,cond2)
######################################################################

4.网络训练

#####################################################################
#参数设置
c,len_x,len_y,len_t = 2, 2, 2,6 
t_b,t_e,t_num = 0.0,2,21
x_b,x_e,x_num = 0.0,8,40
y_b,y_e,y_num = 0.0,8,40
sample_num = 10000
######################################################################
#范围设置
t_range = [t_b, t_e]
x_range = [x_b, x_e]
y_range = [y_b, y_e]
lb = np.asarray([x_range[0], y_range[0], t_range[0]])
ub = np.asarray([x_range[1], y_range[1], t_range[1]])
#内点，边界，初始条件采样
xyt_in = LHS_Sampling(sample_num,lb,ub)
xyt_left = LHS_Sampling(sample_num,lb,ub)
xyt_right = LHS_Sampling(sample_num,lb,ub)
xyt_down = LHS_Sampling(sample_num,lb,ub)
xyt_up = LHS_Sampling(sample_num,lb,ub)
xyt_ini1 = LHS_Sampling(sample_num,lb,ub)
xyt_ini2 = LHS_Sampling(sample_num,lb,ub)
######################################################################
######################################################################
#设置训练次数，设置接收损失值的列表
epochs = 10000
loss_i,loss_le,loss_ri,loss_d,loss_u,loss_i1,loss_i2,loss_all = [],[],[],[],[],[],[],[]
######################################################################
#训练
for i in range(epochs):
    opt.zero_grad()
    a1= interior_loss(xyt_in)
    a2 = left_loss(xyt_left)
    a3 = right_loss(xyt_right)
    a4 = down_loss(xyt_down)
    a5 = up_loss(xyt_up)
    a6 = ini1_loss(xyt_ini1)
    a7 = ini2_loss(xyt_ini2)
    l = a1+a2+a3+a4+a5+a6+a6
    loss_all.append(l)
    l.backward()
    opt.step()
    if i % 100 == 0:
        print(i,l)
c8 = Tran_fun(loss_all)

5.网络保存

##################################################################
model = pinn
#torch.save(pinn,"D:/trl/波场结构/add_resoure_wave_model2d.pth")
##################################################################

6.绘图函数编写

这里给了两个绘图函数，可以画平面图，可以画3D图

######################################################################
#编写绘图
class Output_Fun():
    def __init__(self, t_b,t_e,t_num,x_b,x_e,x_num,y_b,y_e,y_num):
        self.t_b = t_b
        self.t_e = t_e
        self.t_num = t_num
        self.x_b = x_b
        self.x_e = x_e
        self.x_num = x_num
        self.y_b = y_b
        self.y_e = y_e
        self.y_num = y_num
        
    def U_Pred(self,model_path):
        model = torch.load(model_path,map_location=torch.device('cpu'))
        tc = torch.linspace(self.t_b, self.t_e, self.t_num)
        xc = torch.linspace(self.x_b, self.x_e, self.x_num)
        yc = torch.linspace(self.y_b, self.y_e, self.y_num)
        tm, xm, ym = torch.meshgrid(tc, xc, yc)
        xx = xm.reshape(-1, 1)
        yy = ym.reshape(-1, 1)
        tt = tm.reshape(-1, 1)
        xy = torch.cat([tt,xx,yy], dim=1)
        u_pred = model(xy)
        return u_pred

    def U_Draw(self,pred_mat):
        plt.rcParams['font.sans-serif']=["Times New Roman"] # 用来正常显示中文标签
        plt.rcParams['axes.unicode_minus']=False # 用来正常显示
        my_font_label = {"family":"Microsoft JhengHei", "size":40}
        my_font_title = {"family":"Microsoft JhengHei", "size":50}
        my_font_ax = {"family":"Times New Roman", "size":40}
        x = torch.linspace(self.x_b,self.x_e,self.x_num)
        y = torch.linspace(self.y_b,self.y_e,self.y_num)
        X,Y = torch.meshgrid(x,y)
        dt = (self.t_e-self.t_b)/int(self.t_num-1)
        pred_mat = pred_mat.detach().numpy()
        lenspace = int(self.x_num**2)
        frames_pre = []
        for i in range(0,int(self.t_num),1):
            u_i = pred_mat[i*lenspace:(i+1)*lenspace,:].reshape((self.x_num,self.y_num))
            fig = plt.figure(figsize=(10,10),tight_layout=True)
            t = round(dt*i,3)
            
            ax = fig.add_subplot(111,projection='3d')
            ax.set_zlim3d(zmin=-3, zmax=3)
            ax.set_xlabel('X',fontdict=my_font_ax,labelpad=25)
            ax.set_ylabel('Y',fontdict=my_font_ax,labelpad=25)
            ax.set_zlabel('U',fontdict=my_font_ax,labelpad=5)
            #ax.set_title("第"+str(t)+"秒波场预测图像",fontdict=my_font_title)
            ax.tick_params(labelsize=30)
            
            surf = ax.plot_trisurf(X.reshape(-1), Y.reshape(-1), u_i.reshape(-1),cmap='RdYlBu_r',vmin =-1,vmax=1)
            cb=fig.colorbar(surf,shrink=0.6,aspect=20,ticks=np.linspace(-1,1,5))
            cb.ax.tick_params(labelsize=30)  #设置色标刻度字体大小。
            cb.set_label('波场值',fontdict=my_font_label) #设置colorbar的标签字体及其大小
            
            
            
            plt.savefig("D:/trl/图像波场/含震源二维预测波场/第"+str(t)+"秒波场图像.jpg")
            #image_name = "D:/trl/图像波场/含震源二维预测波场/第"+str(t)+"秒波场图像.jpg"
            #frames_pre.append(imageio.imread(image_name))
        #imageio.mimsave("D:/trl/图像波场/含震源二维预测波场/pred.gif", frames_pre, 'GIF', duration=0.2) # 保存在当前文件夹
######################################################################
#编写绘图
class Output_Fun_plane():
    def __init__(self, t_b,t_e,t_num,x_b,x_e,x_num,y_b,y_e,y_num):
        self.t_b = t_b
        self.t_e = t_e
        self.t_num = t_num
        self.x_b = x_b
        self.x_e = x_e
        self.x_num = x_num
        self.y_b = y_b
        self.y_e = y_e
        self.y_num = y_num
        
    def U_Pred(self,model_path):
        model = torch.load(model_path,map_location=torch.device('cpu'))
        tc = torch.linspace(self.t_b, self.t_e, self.t_num)
        xc = torch.linspace(self.x_b, self.x_e, self.x_num)
        yc = torch.linspace(self.y_b, self.y_e, self.y_num)
        tm, xm, ym = torch.meshgrid(tc, xc, yc)
        xx = xm.reshape(-1, 1)
        yy = ym.reshape(-1, 1)
        tt = tm.reshape(-1, 1)
        xy = torch.cat([tt,xx,yy], dim=1)
        u_pred = model(xy)
        return u_pred

    def U_Draw(self,pred_mat):
        plt.rcParams['font.sans-serif']=['Times New Roman'] # 用来正常显示中文标签
        plt.rcParams['axes.unicode_minus']=False # 用来正常显示
        my_font_label = {"family":"Microsoft JhengHei", "size":40}
        my_font_title = {"family":"Microsoft JhengHei", "size":50}
        my_font_ax = {"family":"Times New Roman", "size":40}
        x = torch.linspace(self.x_b,self.x_e,self.x_num)
        y = torch.linspace(self.y_b,self.y_e,self.y_num)
        X,Y = torch.meshgrid(x,y)
        dt = (self.t_e-self.t_b)/int(self.t_num-1)
        pred_mat = pred_mat.detach().numpy()
        lenspace = int(self.x_num**2)
        frames_pre = []
        for i in range(0,int(self.t_num),1):
            u_i = pred_mat[i*lenspace:(i+1)*lenspace,:].reshape((self.x_num,self.y_num))
            fig = plt.figure(figsize=(12,10))
            t = round(dt*i,3)
            
            ax = fig.add_subplot(111)
            ax.set_xlabel('X',fontdict=my_font_ax,labelpad=25)
            ax.set_ylabel('Y',fontdict=my_font_ax,labelpad=25)
            ax.set_title("第"+str(t)+"秒波场预测图像",fontdict=my_font_title)
            ax.tick_params(labelsize=30)
            
            surf = ax.pcolor(X, Y, u_i,cmap='RdYlBu_r',vmin =-0.3,vmax=0.3)
            cb=fig.colorbar(surf,shrink=1,aspect=20,ticks=np.linspace(-0.3,0.3,7))
            cb.ax.tick_params(labelsize=30)  #设置色标刻度字体大小。
            cb.set_label('波场值',fontdict=my_font_label) #设置colorbar的标签字体及其大小
            
            
            plt.savefig("D:/trl/图像波场/含震源二维预测波场/第"+str(t)+"秒波场图像二维.jpg")

5.网络保存

#求解方程并绘制3D图像
t_b,t_e,t_num = 0.0,2,21
x_b,x_e,x_num = 0.0,8,40
y_b,y_e,y_num = 0.0,8,40
sol = Output_Fun(t_b,t_e,t_num,x_b,x_e,x_num,y_b,y_e,y_num)
mat1 = sol.U_Pred("D:/trl/波场结构/add_resoure_wave_model2d.pth")
sol.U_Draw(mat1)
######################################################################
#求解方程并绘制图像
#sol = Output_Fun_plane(t_b,t_e,t_num,x_b,x_e,x_num,y_b,y_e,y_num)
#mat1 = sol.U_Pred("D:/trl/波场结构/add_resoure_wave_model2d.pth")
#sol.U_Draw(mat1)
######################################################################

最后的话

这是一个均匀层状介质模型的实现，但是我一直没有调到非常满意的结果，可以看看这两张gif图片，效果就是这样。下一期会给出一个层状介质的模拟结果，也是这个系列的最后一期了。
在这里插入图片描述
在这里插入图片描述

邓肯黑火
关注
  • 6
    点赞
  • 踩
  • 17
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
基于Python的二维有限元声波方程正演计算
fortunemouse的博客
11-29 1238
基于Python的二维有限元声波方程正演计算
声波正演c语言程序,二维频率域声波方程正演模拟
weixin_33343700的博客
05-21 1774
1. 概述频率域波场正演相对于时间域数值模拟来说，有其自身的优势。首先，在多炮数值模拟情况下，频率域相对于时间域效率更高，每个频率成分的阻抗矩阵只需要计算一次，加入并行计算后可以极大地提高计算效率；其次，就方程本身而言，频率域的衰减项更容易添加；另外，频率域模拟方法是计算区域所有网格点上的全部时间上的频率切片解，误差将分配到所有参与计算的每一个网格点上，因此不存在类似时间域模拟的累计误差，可以考虑...
二维波动方程数值模拟(非均匀介质)
Sumbrella_的博客
12-09 5841
文章用有限差分法模拟并可视化了非均匀介质波场，可以观测到反射波、折射波等。
一个简单的二维波动方程有限差分正演示例（均匀介质，无吸收边界条件及震源）
Reading8的博客
09-05 700
均匀介质中无震源及吸收边界条件的二维波动方程有限差分正演
二维声波传播方程的有限差分模拟
DSPSTACK.COM
09-09 3419
二维声波传播方程的有限差分解法 二维声波方程在Oxz平面表示： 有限差分表示： 其中f(t)表示源函数，我们用Ricker作为激发源。 离散化的二维声波方程 matlab示例 x,z向共201个节点，节点间隔h=8m,时间采样点位400，采样间隔为0.001s。假设声音传播速度为3km/s，激发源在i=100,j=100处。Ricker主频为20Hz,频带控制参数r=3. clc; c...
二维声波方程的有限差分法数值模拟
Neverlevsun的博客
08-11 5761
二维声波方程的有限差分法数值模拟 文章目录二维声波方程的有限差分法数值模拟一、实现效果二、代码分享 一、实现效果 二、代码分享 close all clear clc % 此程序是有限差分法实现声波方程数值模拟 %% 参数设置 delta_t = 0.001; % s delta_s = 10; % 空间差分：delta_s = delta_x = delta_y (m) nx = 800; ny = 800; nt = 1000; fmain = 12.5; %loop:按照10000s为一次大
波动方程的PINN方法求解
qq_58509609的博客
04-27 1419
波动方程的求解一直是地球物理领域中的一大重要问题。波动方程的求解有多种方法，传统的方法有差分法、谱元法、伪谱法等，此外利用深度学习算法求解偏微分方程的方法在近几年也出了很多文章。这段时间刚好在做毕业设计，内容是基于深度学习算法求解Biot方程，其中简单二维波动方程的求解是至关重要的一步，因此将计算的流程和代码分享出来，供给大家参考。
基于Python的二维谱元法声波方程正演计算
fortunemouse的博客
12-04 1287
基于Python的二维谱元法声波方程正演计算
modeling_pre_2D_y (2)_波动方程正演_pre_二维正演程序_
10-01
标题中的"modeling_pre_2D_y (2)_波动方程正演_pre_二维正演程序_"表明这是一个关于二维波动方程正演模拟的程序，主要用于地质勘探领域，特别是地震成像技术。在这个程序中，"pre"通常代表"前处理"，意味着它可能...
FD_规则网格_二维声波方程规则网格有限差分正演_声波有限差分_声波差分_声波规则_
10-04
在IT领域，尤其是在地球物理勘探、信号处理以及计算科学中，有限差分法（Finite Difference Method, FDM）是一种常用的技术，用于数值求解偏微分方程，如本问题中的二维声波方程。标题"FD_规则网格_二维声波方程规则...
基于matlab二维射线追踪程序地震声波正演源程序
04-18
程序需要一个二维的速度模型作为输入，这通常由地震测井数据或地震波的逆向问题求解得到。 4. **源项设置**：在正演模拟中，需要定义震源的性质，如震源位置、震级、辐射模式等。这些参数会影响模拟出的地震波形。 ...
【项目代码】二维射线追踪程序！matlab编写！！ 地震声波正演源程序.rar
08-19
二维射线追踪是一种在地球物理学领域中用于模拟地震波传播的重要技术，特别是在地震勘探和地壳结构分析中占据核心地位。Matlab作为一种强大的数学计算和可视化工具，是编写这种复杂算法的理想选择，因为它提供了丰富...
二维声波正演 自激自收，二维声波方程有限差分正演，PML吸收边界，
09-02
声波方程描述了压力波动如何在介质中传播，而有限差分法是一种数值计算方法，用于将偏微分方程离散化，进而转化为可以计算机求解的代数方程组。在二维情况下，这意味着我们需要考虑水平和垂直两个方向的波传播。通过...
无监督语义分割综述
qq_46981910的博客
07-16 314
语义分割是计算机视觉领域的一个重要任务，旨在将图像中的每个像素分配给特定的语义类别。然而，传统的语义分割方法通常依赖大量标注数据，获取这些数据既费时又昂贵。无监督语义分割方法旨在通过不依赖标注数据或仅需少量标注数据的情况下，实现对图像的精确分割。这些方法通过利用数据的内在结构、图像特征和自监督学习技术，从大量未标注的数据中提取有用的信息，显著降低了对标注数据的依赖。本文综述了无监督语义分割的主要方法，包括基于聚类的方法、基于图割的方法、基于生成模型的方法、基于自监督学习的方法以及基于领域自适应性的方法。
机器学习 第五章
zxsdss的博客
07-16 647
BP算法的目标是最小化训练集D上的累积误差，基于累积误差最小化的更新规则，就可以得到累积误差逆传播算法，即每次读取全部的数据集一遍，进行一轮学习，从而基于当前的累积误差进行权值调整，其参数更新的频率相比标准BP算法低了很多。基于梯度的搜索是使用最广泛的参数寻优方法，梯度下降法的主要思想就是沿着负梯度方向去搜索最优解，负梯度方向是函数值下降最快的方向，若迭代到某处的梯度为0，则表示达到一个局部最小，参数更新停止。机器学习中的神经网络指的是神经网络学习，即机器学习与神经网络两个学科领域的交叉部分。
YOLO--置信度（超详细解读）
m0_71212744的博客
07-16 745
数值范围：置信度是一个介于0和1之间的数值。意义：它表示模型对预测框内存在目标对象的确信程度。如果置信度接近1，表示模型非常确信预测框内包含了目标对象；如果置信度接近0，则表示模型认为预测框内可能不包含目标对象。
蒙特卡洛采样
最新发布
weixin_46713695的博客
07-16 397
蒙特卡洛采样
写一个水平分层介质的二维波动方程正演模拟matlab代码，包含边界条件，震源项，收敛条件，并解释每一行代码的意义
05-13
以下是一个二维水平分层介质的波动方程正演模拟的 MATLAB 代码。本代码考虑了各向异性介质，使用了有限差分时间域（FDTD）方法。 ```matlab % 声波二维正演模拟 clc; clear; % 设定模拟区域大小 xmax = 1000; % x方向最大值（单位：m） ymax = 1000; % y方向最大值（单位：m） dx = 10; % 网格间距（单位：m） dy = 10; % 网格间距（单位：m） nx = floor(xmax/dx); % x方向网格数 ny = floor(ymax/dy); % y方向网格数 % 设定模拟时间参数 tmax = 2; % 最大模拟时间（单位：s） dt = 0.001; % 时间步长（单位：s） nt = floor(tmax/dt); % 时间步数 % 设定介质参数 vp = zeros(nx, ny); % 声波速度（单位：m/s） vs = zeros(nx, ny); % 纵波速度（单位：m/s） rho = zeros(nx, ny); % 密度（单位：kg/m^3） epsilon = zeros(nx, ny); % 纵波各向异性系数 delta = zeros(nx, ny); % 横波各向异性系数 % 设定震源参数 sx = nx/2; % 震源位置x坐标 sy = ny/2; % 震源位置y坐标 f0 = 30; % 震源频率 t0 = 0.05; % 震源持续时间 src = zeros(nt,1); % 震源波形 % 生成Ricker子波作为震源波形 for i = 1:nt t = (i-1)*dt; src(i) = (1-2*pi^2*f0^2*(t-t0)^2)*exp(-pi^2*f0^2*(t-t0)^2); end % 初始化波场和边界条件 p = zeros(nx,ny); % 压力波场 vx = zeros(nx,ny); % x方向速度波场 vy = zeros(nx,ny); % y方向速度波场 pml_len = 40; % PML吸收边界层数 pml_coef = 0.001; % PML吸收系数 pml_xmin = zeros(pml_len,ny); % PML吸收边界x方向左侧 pml_xmax = zeros(pml_len,ny); % PML吸收边界x方向右侧 pml_ymin = zeros(nx,pml_len); % PML吸收边界y方向下侧 pml_ymax = zeros(nx,pml_len); % PML吸收边界y方向上侧 % 循环进行时间步进计算 for it = 1:nt % 计算PML吸收边界 if pml_len > 0 for i = 1:pml_len x = (pml_len-i+0.5)*dx; sigmax = pml_coef*(i/pml_len)^2; for j = 1:ny pml_xmin(i,j) = (1-sigmax)*pml_xmin(i,j) + sigmax*p(i,j); pml_xmax(i,j) = (1-sigmax)*pml_xmax(i,j) + sigmax*p(nx-i+1,j); end end for i = 1:nx y = (pml_len-i+0.5)*dy; sigmay = pml_coef*(i/pml_len)^2; for j = 1:pml_len pml_ymin(i,j) = (1-sigmay)*pml_ymin(i,j) + sigmay*p(i,j); pml_ymax(i,j) = (1-sigmay)*pml_ymax(i,j) + sigmay*p(i,ny-j+1); end end end % 计算应力源 if it <= length(src) p(sx,sy) = p(sx,sy) + src(it); end % 计算速度波场 for i = 2:nx-1 for j = 2:ny-1 % 计算x方向速度 vx(i,j) = vx(i,j) + dt/rho(i,j)*(epsilon(i,j)*... (p(i,j+1)-p(i,j))-delta(i,j)*(p(i+1,j+1)-p(i+1,j-1))/4/dy); % 计算y方向速度 vy(i,j) = vy(i,j) + dt/rho(i,j)*(epsilon(i,j)*... (p(i+1,j)-p(i,j))-delta(i,j)*(p(i+1,j+1)-p(i-1,j+1))/4/dx); end end % 计算应力波场 for i = 2:nx-1 for j = 2:ny-1 % 计算x方向应力 p(i,j) = p(i,j) + dt*vp(i,j)^2*rho(i,j)*... (epsilon(i,j)*(vx(i,j)-vx(i,j-1))/dy+delta(i,j)*... (vy(i,j)-vy(i-1,j))/dx); end end % 边界处理 if pml_len > 0 p(1:pml_len,:) = pml_xmin; p(nx-pml_len+1:nx,:) = pml_xmax; p(:,1:pml_len) = pml_ymin; p(:,ny-pml_len+1:ny) = pml_ymax; end % 显示进度条 fprintf('Progress: %.2f%%\n', it/nt*100); end ``` 代码解释如下： 1. 第1-3行：清空命令窗口，清除已经定义的变量。 2. 第6-10行：设定模拟区域大小和网格间距。 3. 第13-15行：设定模拟时间参数。 4. 第18-22行：设定介质参数，包括声波速度、纵波速度、密度、纵波各向异性系数和横波各向异性系数。 5. 第25-28行：设定震源参数，包括震源位置、震源频率和震源持续时间。 6. 第31-36行：生成Ricker子波作为震源波形。 7. 第39-48行：初始化波场和边界条件。波场包括压力波场和速度波场，边界条件包括PML吸收边界。 8. 第51-93行：循环进行时间步进计算。首先计算PML吸收边界，然后计算应力源，接着计算速度波场，最后计算应力波场。边界处理使用PML吸收边界。在每个时间步长结束后，显示进度条。 以上就是一个二维水平分层介质的波动方程正演模拟的 MATLAB 代码。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值