正演培训视频学习总结6.25-6.26

shengMio

已于 2023-10-16 23:11:17 修改

阅读量120

点赞数

分类专栏：视频学习总结文章标签：学习

于 2023-10-16 22:27:54 首次发布

本文链接：https://blog.csdn.net/m0_56382522/article/details/133846231

版权

视频学习总结专栏收录该内容

5 篇文章 0 订阅

订阅专栏

1. 地震正演

地震正演模拟是一种以地震波传播理论为基础，在已知地下介质结构的情况下，模拟地震波在介质中的传播过程的研究方法。地震波场正演模拟主要分为物理模拟和数值两种方法。

物理模拟通过进行大量物理模拟实验来研究波的传播。
已发展的数值模拟方法大致可被分为射线追踪法（快）以及基于波动方程的微分方程法和积分方程法。微分方程法主要包含有限差分法（应用广）、有限元法、谱元法等。

2. 科学计算

科学与工程计算最重要的内容就是求解各种偏微分方程或方程组。例如核反应和核爆炸的过程：巨大能量在极短时间内爆发，核装置内部的细致反应过程不能用仪器测量出来，核试验只提供综合的数据。而描述核反应和爆炸物理过程的数学模型是一个很复杂的非线性偏微分方程组，没办法得到该方程组的精确解，所以发展核武器的国家都在计算机上对核反应过程进行数值模拟。采用数值风洞求解空气动力学的偏微分方程可减少风洞实验次数。
解析解：模型简单
数值解：非精确解

3. 波动方程

3.1 弦振动方程

$\frac{\partial ^2u }{\partial t^2} = a^{2}\frac{\partial ^2u }{\partial x^2} + F(x,t)$

数学物理方程讲义（姜礼尚）弦振动方程-CSDN博客

偏微分1 | 弦振动方程 - 知乎 (zhihu.com)

3.2 二维声波方程（声波只有纵波）

$\frac{1}{v^{2}}\frac{\partial^2U }{\partial t^2} = \frac{\partial^2U }{\partial x^2} + \frac{\partial^2U }{\partial z^2}$

注释：U——声压即压力；x、z——方向

上式 = 一阶速度-应力格式：

$\left\{\begin{matrix} \frac{\partial v_{x}}{\partial t} = \frac{1}{\varrho }\frac{\partial P}{\partial x}, \\ \frac{\partial v_{z}}{\partial t} = \frac{1}{\varrho }\frac{\partial P}{\partial y}, \\ \frac{\partial P}{\partial t} = \rho v^{2}(\frac{\partial v_{x}}{\partial x} + \frac{\partial v_{z}}{\partial y}), \end{matrix}\right.$

注释：P——应力； $v_{x}$ ——速度

3.3 二维弹性波方程

弹性波传播中，能量在弹性物体的受力、变形与平衡过程中传递。因此建立弹性波波动方程的关键是研究弹性体的应力、应变和位移三者间关系。

地球实际上是弹性介质，可按传播方式被分为纵波、横波部分，传播过程较为复杂。弹性理论来源最早可以追溯到1660年的Hook定律，在1821年 Navier用于平衡方程和振动方程的研究。而后在1822年，Cauchy建立了6个应力分析及6个应变分量的理论。Poisson根据牛顿吻体内分子间作用力理论，发现了Р波及S波。

$\left\{\begin{matrix} \rho \frac{\partial v_{x}}{\partial t} = \frac{\partial \tau _{xx}}{\partial x} + \frac{\partial \tau _{xz}}{\partial z} \\ \rho \frac{\partial v_{z}}{\partial t} = \frac{\partial \tau _{xz}}{\partial x} + \frac{\partial \tau _{zz}}{\partial z} \\ \frac{\partial \tau _{xx}}{\partial t} = c_{11}\frac{\partial v_{x}}{\partial x} + c_{13}\frac{\partial v_{z}}{\partial z} \\ \frac{\partial \tau _{zz}}{\partial t} = c_{13}\frac{\partial v_{x}}{\partial x} + c_{33}\frac{\partial v_{z}}{\partial z} \\ \frac{\partial \tau _{xz}}{\partial t} = c_{44}\frac{\partial v_{z}}{\partial x} + c_{44}\frac{\partial v_{x}}{\partial z} \end{matrix}\right.$

其中， $v_{x}$ 、 $v_{z}$ 分别为:x、z方向速度分量， $\tau$ 为应力， $\rho$ 为密度， $c_{ij}$ 为介质的弹性常数，在各向同性介质中， $c_{11}=c_{33}=\lambda +2\mu ,c_{13}=\lambda ,c_{44}=\mu$ 。

补充：三维弹性波波动方程推导如下：

为什么各向异性后，为什么弹性系数矩阵中左下和右上部分会被省略？（未解决）

如果各个方向的测量结果是相同的，说明其物理性质与取向无关,就称为各向同性。

波场快照：在指定时刻波的信息，即波的传播情况，如每个空间点现在的压力、振幅是多少

Laplace算子：

$\Delta =\frac{\partial^2 }{\partial x_{1}^2} + \frac{\partial^2 }{\partial x_{2}^2} + ... + \frac{\partial^2 }{\partial x_{n}^2}$

$\bigtriangledown = e_{1}\frac{\partial }{\partial x_{1}} + e_{2}\frac{\partial }{\partial x_{2}} + ...+ e_{n}\frac{\partial }{\partial x_{n}}$ e：正向单位向量

4. 建立差分格式

4.1 网格剖分

4.2 泰勒级数展开

其中，k代表时间层，О表示同阶项，o表示高阶无穷小，f(x)看成U。忽略高阶项可得到二阶中心差分格式，同理对时间项也进行二阶中心差分，有:

$\frac{\partial^2 U^{k} }{\partial t^2} = \frac{U^{k+1}-2U^{k}+U^{k-1}}{\Delta t^{2}}$
其中， $\Delta x$ 是差分网格的宽度，被称为空间步长(h)。

在声波波动方程上进行近似差分可以得到：即把声波波动方程离散化

$\frac{1}{v^{2}} \frac{U_{i,j}^{k+1}-2U_{i,j}^{k}+U_{i,j}^{k-1}}{\Delta t^{2}} = \frac{U_{i+1,j}^{k}-2U_{i,j}^{k}+U_{i-1,j}^{k}}{\Delta x^{2}} + \frac{U_{i,j+1}^{k}-2U_{i,j}^{k}+U_{i,j-1}^{k}}{\Delta z^{2}}$

以上就是一种用于二维声波波动方程的规则网格差分格式。利用上式，可在已知k及k-1时间层各网格点的声压情况下，计算出k＋1时间层上各网格点的声压。可由前一时间层状态求出后一时间层状态，构成了时间上的正序迭代，因此能够通过计算机模拟声波的传播。

4.3 运行代码

Endtime=0.5;%模拟时长
delta_t=0.0005;%以秒为单位
delta_x=6;%以米为单位  空间步长
delta_z=6;
CNX=301;%x方向的网格数
CNZ=301;%z
v=1500;%波速，m/s
Sx=(CNX+1)/2;%震源位置
Sz=(CNZ+1)/2;
fO=30; % 10~40HZ,震源主频
Unow=zeros(CNX,CNZ);
Uprev=zeros(CNX,CNZ);
Unext=zeros(CNX,CNZ);


%主程序
for time=O:delta_t:Endtime
    for i=2:CNX-1
        for j=2:CNZ-1
        A=(-2*Unow(i,j)+Unow(i+l,j)+Unow(i-1,j))/delta_x^2;
        B=(-2*Unow(i,j)+Unow(i,j+1)+Unow(i,j-1))/delta_z^2;
        Unext(i,j)=2*Unow(i,j)-Uprev(i,j)+v^2*(A+B)*delta_t^2;
        end
    end
    Unext(Sx,Sz)=5.76*f0^2*(1-16*(0.6*fO*time-1)^2)*exp(-8*(0.6*f0*time-1)^2);%震源函数,给出的                          
    uprev=Unow;%波场更新
    Unow=Unext;
end

%绘图
surf(Unow)
shading interp;
view(2); %view (90,90)
colormap(gray);
toc

结果图：