地震正演培训

已于 2023-12-04 17:30:49 修改
阅读量116 收藏 1
点赞数
文章标签： 机器学习
于 2023-07-26 16:26:35 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/Zhang_ing/article/details/131940193
版权

目录

正演的大致流程

首先有一个观测数据，假设一个速度模型velocity mode，由velocity mode正演得到seismic数据，将seismic数据与观测数据进行比较，并由误差更改速度模型

二维声波方程

研究平面声波的传播行为，即 t t t与 x x x、 y y y之间的关系。

常规网格

连续偏导形式

二位声波方程的连续偏导形式如下：
1 v 2 ∂ 2 U ∂ t 2 = ∂ 2 U ∂ x 2 + ∂ 2 U ∂ y 2 \frac{1}{v^2}\frac{\partial^2 U}{\partial t^2}=\frac{\partial^2U}{\partial x^2} + \frac{\partial^2 U}{\partial y^2} v21​∂t2∂2U​=∂x2∂2U​+∂y2∂2U​
其中 U U U是声压（声波的压力变化）， v v v是声速

二维差分形式

在这里插入图片描述
对连续的偏导形式进行差分处理，实现离散化

泰勒级数

f ( x ) − f ( x − h ) = ∂ f ( x ) ∂ x h − 1 2 ∂ 2 f ( x ) ∂ x 2 + o ( h 2 ) f(x)-f(x-h)=\frac{\partial f(x)}{\partial x}h-\frac{1}{2}\frac{\partial^2 f(x)}{\partial x^2}+o(h^2) f(x)−f(x−h)=∂x∂f(x)​h−21​∂x2∂2f(x)​+o(h2)

f ( x + h ) − f ( x ) = ∂ f ( x ) ∂ x h + 1 2 ∂ 2 f ( x ) ∂ x 2 + o ( h 2 ) f(x+h)-f(x)=\frac{\partial f(x)}{\partial x}h+\frac{1}{2}\frac{\partial^2 f(x)}{\partial x^2}+o(h^2) f(x+h)−f(x)=∂x∂f(x)​h+21​∂x2∂2f(x)​+o(h2)

联立得

f ( x + h ) + f ( x − h ) + 2 f ( x ) h 2 = ∂ 2 f ( x ) ∂ x 2 + o ( h 2 ) \frac{f(x+h)+f(x-h)+2f(x)}{h^2}=\frac{\partial^2 f(x)}{\partial x^2}+o(h^2) h2f(x+h)+f(x−h)+2f(x)​=∂x2∂2f(x)​+o(h2)
其中 h h h就是差分单元，那么再结合上面的二位声波方程就可以得到如下的差分后方程
在这里插入图片描述

代码实现

% 声波二维均匀模型-中心差分格式
tic
clc
close all
clear all

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

Unow=zeros(CNX,CNZ);
Uprev=zeros(CNX,CNZ);
Unext=zeros(CNX,CNZ);

% 主程序
for time=0:delta_t:Endtime
  for i=2:CNX-1
      for j=2:CNZ-1
      A=(-2*Unow(i,j)+Unow(i+1,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*f0*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

以上是有限差分实现代码，图像如下
在这里插入图片描述

边界网格

地震波传播到介质中断层处时，界面的波阻抗出现差异导致波发生反射和透射，是一种自然的物理现象。在自然传播情况中，地球介质可被视为无限大的区域，但在数值模拟中，受到计算量的限制，可模拟的地震波传播区域大小受限。因而模拟计算区域必定存在人工边界。

代码

% 声波二阶-添加边界-常规网格
tic
clc
close all
clear all

Endtime=0.4;  %模拟时长
delta_t=0.001;%以秒为单位
delta_x=0.006;%以千米为单位   空间步长
delta_z=0.006;
f0=30;% 10~40HZ  震源频率
iteration=1;
nk=1;
n=Endtime/delta_t;

CNX=301; % 
CNZ=301; %
BNX=20;
BNZ=20;
XMAX=CNX+2*BNX;
ZMAX=CNZ+2*BNZ;
Sx=(XMAX+1)/2; %震源位置
Sz=(ZMAX+1)/2;
c=0;

%zeros(XMAX,ZMAX);%波速
v(1:100,1:ZMAX)=2.0;
v(101:200,1:ZMAX)=2.0;
v(201:XMAX,1:ZMAX)=2;

a=-0.15;
S=0.001;

Wave=zeros(n+1,1);
Fnow=zeros(n+1,CNZ);
Unow=zeros(XMAX,ZMAX);
Uprev=zeros(XMAX,ZMAX);
Unext=zeros(XMAX,ZMAX);

Uxnow=zeros(XMAX,ZMAX);%n
Uxprev=zeros(XMAX,ZMAX);%n-1
Uxnext=zeros(XMAX,ZMAX);%n+1
Uznow=zeros(XMAX,ZMAX);%n
Uzprev=zeros(XMAX,ZMAX);%n-1
Uznext=zeros(XMAX,ZMAX);%n+1

Psi2x=zeros(2*BNZ,XMAX);
Uxxprev=zeros(2*BNZ,XMAX);
Psi2xprev=zeros(2*BNZ,XMAX);
Eta2xprev=zeros(2*BNZ,XMAX);
Eta2x=zeros(2*BNZ,XMAX);
Psixprev=zeros(2*BNZ,XMAX);
Psix=zeros(2*BNZ,XMAX);
Etaxprev=zeros(2*BNZ,XMAX);
Etax=zeros(2*BNZ,XMAX);
Thetaxprev=zeros(2*BNZ,XMAX);
Thetax=zeros(2*BNZ,XMAX);


Uzzprev=zeros(ZMAX,2*BNZ);
Psi2zprev=zeros(ZMAX,2*BNZ);
Psi2z=zeros(ZMAX,2*BNZ);
Eta2zprev=zeros(ZMAX,2*BNZ);
Eta2z=zeros(ZMAX,2*BNZ);
Psizprev=zeros(ZMAX,2*BNZ);
Psiz=zeros(ZMAX,2*BNZ);
Etazprev=zeros(ZMAX,2*BNZ);
Etaz=zeros(ZMAX,2*BNZ);
Thetazprev=zeros(ZMAX,2*BNZ);
Thetaz=zeros(ZMAX,2*BNZ);

% 主程序
for time=0:delta_t:Endtime
    for m=4:XMAX-3
      for n=4:ZMAX-3
          c=v(m,n);
     %A=(-2*Unow(m,n)+Unow(m+1,n)+Unow(m-1,n))/delta_x^2;
     %B=(-2*Unow(m,n)+Unow(m,n+1)+Unow(m,n-1))/delta_z^2;
     A=(-2.8194*Unow(m,n)+1.575*(Unow(m+1,n)+Unow(m-1,n))-0.1827*(Unow(m+2,n)+Unow(m-2,n))+0.0174*(Unow(m+3,n)+Unow(m-3,n)))/delta_x^2;
     B=(-2.8194*Unow(m,n)+1.575*(Unow(m,n+1)+Unow(m,n-1))-0.1827*(Unow(m,n+2)+Unow(m,n-2))+0.0174*(Unow(m,n+3)+Unow(m,n-3)))/delta_z^2;
     Unext(m,n)=2*Unow(m,n)-Uprev(m,n)+c^2*(A+B)*delta_t^2;

     %fprintf('their is %8.5f\n',c);
     
     %***************添加边界***************
        if m<=BNX+1 || m>=BNX+CNX || n<=BNZ+1 || n>=BNZ+CNZ
            
              Uxx=A;
              Uzz=B;
  
            %U的X方向****************************
            if m<=BNX+1 || m>=CNX+BNX
                if m<=BNX+1
                    XBPos=m;
                    XBDistance=(BNX+1-m)*delta_x;
                else
                    XBPos=m-CNX;
                    XBDistance=(m-(BNX+CNX))*delta_x;
                end
                ZBPos=n;
                
                XBThick=BNX*delta_x ;
                XDecay=-3*c*log(S)/(2*XBThick)*(XBDistance/XBThick)^2;
                XDecay_order=-3*c*log(S)*XBDistance/XBThick^3;
                XAlpha=pi*f0*(1-XBDistance/XBThick);
                XAlpha_order=-pi*f0/XBThick;
                bX=exp(-(XDecay+XAlpha)*delta_t);
                aX=(1-bX)/(XDecay+XAlpha);    
            
                %傅里叶逆变换各部分表达式
                Psi2x(XBPos,ZBPos)=bX*Psi2xprev(XBPos,ZBPos)+aX*0.5*(Uxx+Uxxprev(XBPos,ZBPos));
                Eta2x(XBPos,ZBPos)=bX*Eta2xprev(XBPos,ZBPos)+aX*0.5*(Psi2x(XBPos,ZBPos)+Psi2xprev(XBPos,ZBPos));
                Psix(XBPos,ZBPos)=bX*Psixprev(XBPos,ZBPos)+aX*0.5*(Uxnow(m,n)+Uxprev(m,n));
                Etax(XBPos,ZBPos)=bX*Etaxprev(XBPos,ZBPos)+aX*0.5*(Psix(XBPos,ZBPos)+Psixprev(XBPos,ZBPos))    ;
                Thetax(XBPos,ZBPos)=bX*Thetaxprev(XBPos,ZBPos)+aX*0.5*(Etax(XBPos,ZBPos)+Etaxprev(XBPos,ZBPos)) ;
                
                %指标变换——迭代更新
                Uxxprev(XBPos,ZBPos)=Uxx;
                Psi2xprev(XBPos,ZBPos)=Psi2x(XBPos,ZBPos);
                Eta2xprev(XBPos,ZBPos)=Eta2x(XBPos,ZBPos);
                Psixprev(XBPos,ZBPos)=Psix(XBPos,ZBPos);
                Etaxprev(XBPos,ZBPos)=Etax(XBPos,ZBPos);
                Thetaxprev(XBPos,ZBPos)=Thetax(XBPos,ZBPos);

                %U的X方向二阶偏导数的傅里叶逆变换表达式
                Uxx=Uxx-2*XDecay*Psi2x(XBPos,ZBPos)+XDecay^2*Eta2x(XBPos,ZBPos)-XDecay_order*Psix(XBPos,ZBPos)...
                    +XDecay*(2*XDecay_order+XAlpha_order)*Etax(XBPos,ZBPos)-XDecay^2*(XDecay_order+XAlpha_order)*Thetax(XBPos,ZBPos);
            end
%   **************U的Z方向****************************
            if n<=BNZ+1 || n>=BNZ+CNZ
                if n<=BNZ+1
                    ZBPos=n;
                    ZBDistance=(BNZ+1-n)*delta_z;
                else
                    ZBPos=n-CNZ;
                    ZBDistance=(n-(BNZ+CNZ))*delta_z;
                end
                XBPos=m;
                
                ZBThick=BNZ*delta_z;
                ZDecay=-3*c/(2*ZBThick)*log(S)*(ZBDistance/ZBThick)^2;
                ZDecay_order=-3*c*log(S)*ZBDistance/ZBThick^3;
                ZAlpha=pi*f0*(1-ZBDistance/ZBThick);
                ZAlpha_order=-pi*f0/ZBThick;
                bZ=exp(-(ZDecay+ZAlpha)*delta_t);
                aZ=(1-bZ)/(ZDecay+ZAlpha);   
                
                %傅里叶逆变换各部分表达式
                Psi2z(XBPos,ZBPos)=bZ*Psi2zprev(XBPos,ZBPos)+aZ*0.5*(Uzz+Uzzprev(XBPos,ZBPos));
                Eta2z(XBPos,ZBPos)=bZ*Eta2zprev(XBPos,ZBPos)+aZ*0.5*(Psi2z(XBPos,ZBPos)+Psi2zprev(XBPos,ZBPos));
                Psiz(XBPos,ZBPos)=bZ*Psizprev(XBPos,ZBPos)+aZ*0.5*(Uznow(m,n)+Uzprev(m,n));
                Etaz(XBPos,ZBPos)=bZ*Etazprev(XBPos,ZBPos)+aZ*0.5*(Psiz(XBPos,ZBPos)+Psizprev(XBPos,ZBPos)) ;
                Thetaz(XBPos,ZBPos)=bZ*Thetazprev(XBPos,ZBPos)+aZ*0.5*(Etaz(XBPos,ZBPos)+Etazprev(XBPos,ZBPos)) ;
                
                %指标变换——迭代更新
                Uzzprev(XBPos,ZBPos)=Uzz;
                Psi2zprev(XBPos,ZBPos)=Psi2z(XBPos,ZBPos);
                Eta2zprev(XBPos,ZBPos)=Eta2z(XBPos,ZBPos);
                Psizprev(XBPos,ZBPos)=Psiz(XBPos,ZBPos);
                Etazprev(XBPos,ZBPos)=Etaz(XBPos,ZBPos);
                Thetazprev(XBPos,ZBPos)=Thetaz(XBPos,ZBPos);
                
                %Z方向二阶偏导数的傅里叶逆变换表达式
                Uzz=Uzz-2*ZDecay*Psi2z(XBPos,ZBPos)+ZDecay^2*Eta2z(XBPos,ZBPos)-ZDecay_order*Psiz(XBPos,ZBPos)...
                    +ZDecay*(2*ZDecay_order+ZAlpha_order)*Etaz(XBPos,ZBPos)-ZDecay^2*(ZDecay_order+ZAlpha_order)*Thetaz(XBPos,ZBPos);
            end
            Unext(m,n)=2*Unow(m,n)-Uprev(m,n)+(c*delta_t)^2*(Uxx+Uzz);
        end %边界循环结束
            if m == Sx
                Fnow(nk,n) = Unext(Sx,n);%地震记录
            end
              Wave(iteration) = Unext(101,151); %波形图
              iteration = iteration + 1;
      end
   end
  Unext(Sx,Sz)=delta_t^2*5.76*f0^2*(1-16*(0.6*f0*time-1)^2)*exp(-8*(0.6*f0*time-1)^2);

  Uprev=Unow;  %波场更新
  Unow=Unext;
  nk=nk+1;
end

%绘图
%plot(Wave)
%surf(Unow)
surf(Fnow)

shading interp;
view(2);%view(90,90)
colormap(gray);

toc

得到三个结果：
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

交错网格

网格划分

在这里插入图片描述

一阶速度应力格式

在这里插入图片描述
按照上面提到过的泰勒级数处理，得到
在这里插入图片描述

代码

% 声波二维_一阶速度-应力格式_交错网格
tic
clc
close all
clear all

Endtime=0.4;  %模拟时长
delta_t=0.0005;%以秒为单位
delta_x=0.008;%以千米为单位   空间步长
delta_z=0.008;
CNX=301; %x方向的网格数
CNZ=301; %
rho=1.2;%
v=2; %波速
Sx=(CNX+1)/2; %震源位置
Sz=(CNZ+1)/2;
f0=30;% 10~40HZ  震源频率

Unow=zeros(CNX,CNZ);
Unext=zeros(CNX,CNZ);
Vnow=zeros(CNX,CNZ);
Vnext=zeros(CNX,CNZ);
Pnow=zeros(CNX,CNZ);
Pnext=zeros(CNX,CNZ);


% 主程序
for time=0:delta_t:Endtime
    for i=2:CNX-1
        for j=2:CNZ-1
            M=(Pnow(i,j)-Pnow(i-1,j))/(rho*delta_x);
            N=(Pnow(i,j)-Pnow(i,j-1))/(rho*delta_z);
            Unext(i,j) = Unow(i,j) + delta_t*M;
            Vnext(i,j) = Vnow(i,j) + delta_t*N;
        end
    end
    
    for i=1:CNX-1
        for j=1:CNZ-1
            X = (Unext(i+1,j) - Unext(i,j))/delta_x;
            Y = (Vnext(i,j+1) - Vnext(i,j))/delta_z;
            Pnext(i,j)=Pnow(i,j) + rho*v*v*(X+Y)*delta_t;
        end
    end

  Pnext(Sx,Sz)=5.76*f0^2*(1-16*(0.6*f0*time-1)^2)*exp(-8*(0.6*f0*time-1)^2);
  
  %波场更新
  Unow=Unext;  %波场更新
  Vnow=Vnext; %波场更新
  Pnow=Pnext;
end
% max(max(Vnext))
%绘图
surf(Pnow)
%surf(Unow)%x方向速度分量
%surf(Vnow)%z方向速度分量
shading interp;

view(2);%view(90,90)
colormap(gray);
toc

结果：
在这里插入图片描述

弹性波

outragee
关注
  • 0
    点赞
  • 踩
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
Norsar 3D v5.44 Linux 1CD(功能强大的正演软件)
cadworx2016的博客
05-11 1714
JKBench v1.15 1CD（爆破设计与分析） Norsar 3D v5.44 Linux 1CD(功能强大的正演软件) Norsar 2D v5.16 Linux 1CD Norsar 培训教程 Petrosys v16.8-ISO 1CD Physprops v1.6.1 1CD Pipedrop v1.2.6 1CD Powerlog v3.3 FULL 1CD（专
研究报告（2）--基于有限差分的地震正演培训总结
qq_52245889的博客
07-14 442
对地震正演方法的总结
地震正演基础知识
fulishafulisha的博客
07-15 1689
地震正演基础知识
地震正演培训-记录篇
qq_52245889的博客
06-27 376
5天培训的会议记录及作业答案
正演培训视频学习总结6.25-6.26
shengMio的博客
10-16 109
地震正演
地震勘探原理自问自答
10-16 242
2、以陆地勘探为例，炮点激发，将地下介质近似为弹性介质，由激发信号经过地下介质的传播，在地表上采集到的信号称为地震数据，也就是我们常说的正演过程。地震数据包含：体波（反射比、直达波、折射波等）、面波（勒夫波、瑞利波/地滚波）、噪声。比如，是否含有面波？1、二维/三维（陆地/海洋）地震勘探，炮点（激发点）和检波点（接收点）的排布位置如何？2）利用合成数据去做反演，和利用真实数据去做反演，他们的反演过程/结果会有什么差异？3、请标出真实数据中的直达波、反射波、折射波、面波、初至波、噪声？
全波形反演的深度学习方法: 第 2 章 正演
minfanphd的博客
11-27 1368
本章介绍正演的基础知识. 仅供内部培训. 持续修改中.
全波形反演的深度学习方法: 第 4 章 基于正演的 FWI
minfanphd的博客
11-30 2291
本章论述经典的 FWI, 它基于正演方法. 本贴仅供内部培训.
第三周：tesseral 2D软件的使用
qq_41554955的博客
10-18 3496
参考资料： Tesseral中文用户手册.doc , https://max.book118.com/html/2021/0111/5343040031003110.shtm tesseral软件教学-数据处理到成剖面图 , https://www.bilibili.com/video/BV1f64y1U7gj?from=search&seid=4001933724500544312&spm_id_from=333.337.0.0 tesseral软件教学-观测系统的建立 , http
基于matlab 的地震正演模型实现
06-07
地震软件开发中，matlab应用十分广泛，主要作用是编写程序与计算 基于matlab 的地震正演模型实现
地震波动方程正演模拟-程序
04-29
地震正演模拟相关程序，可以用来进行正演得到炮记录和波场快照，非常方便。
地震资料正演软件
12-20
地震正演软件！
tessel地震正演软件
06-04
地震正演程序，俄罗斯人的软件，好用的很啊。大家支持
基于matlab的地震正演模型实现
08-20
内涵两个matlab程序和论文，模型为三层，含油模型和不含油模型...根据褶积理论 ,结合地质模型 ,在 MATLAB 环境中编写程序实现了地震正演。为了最大可 能的与实际情况相符合 ,地震模型中使用的子波是稳 定可实现的子波。
人工智能和机器学习的应用日益广泛，在医疗健康领域的具体应用是什么？
xiao77224li的博客
06-16 908
相较于传统的群体药动学建模方法，基于机器学习构建的模型能更精准地预测血药浓度和给药剂量，从而提高临床精准用药水平，减少不良反应的发生[51]。这些应用展示了AI和ML技术在医疗健康领域的巨大潜力，不仅能够提高医疗服务的质量和效率，还能推动医疗行业向更加智能化、高效化的方向发展。、分析患者相似度以准确预测药物组合的方法，能够提升药物推荐任务的安全性和准确性[57]。人工智能和机器学习在医疗健康领域中的应用非常广泛，不同的疾病预测模型根据其特定的应用场景和目标疾病类型，采用了不同的算法和技术。
机器学习课程复习——主成分分析
最新发布
Winkyyyyyy的博客
06-19 386
介绍了一种典型的无监督降维算法——主成分分析
Python数据分析与机器学习在金融风控中的应用
热门推荐
屿小夏.的知识博客
06-17 1万+
金融风控是金融机构确保其业务健康运行、减少损失的重要手段。随着大数据和人工智能技术的发展，利用Python进行数据分析和机器学习可以为金融风控提供强有力的支持。本文将探讨Python在金融风控中的应用，详细介绍如何利用Python进行数据收集、预处理、机器学习建模和评估，以提升金融风控的准确性和效率。
【机器学习】MSE和MAE损失函数有什么区别？回归问题的常用损失函数
q742971636的博客
06-17 1007
误差处理方式MSE 使用平方处理误差，对大误差更加敏感。MAE 使用绝对值处理误差，对所有误差均一对待。异常值影响MSE 受到异常值的显著影响。MAE 对异常值的影响较小，更加鲁棒。优化和求导MSE 的平方性质使其导数在很多优化算法中计算更为简便。MAE 的绝对值使其在某些优化过程中可能不如 MSE 简单，但在一些鲁棒优化问题中更加有效。
Python地震正演
10-08
关于Python地震正演，我可以为您提供一些信息。 Python地震正演是指使用Python语言编写地震波动传播的数值模拟程序。该程序可以通过计算机模拟地震波的传播过程，帮助地震学家预测地震的发生和地震波在地球内部的...

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

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

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值