有限差分法波动方程正演模型制作方法研究.ppt

高精度傅里叶有限差分法模型正演刘文革 3 ① 贺振华 ①② 黄德济 ① 杜增利 ③( ①成都理工大学信息工程学院 ; ②成都理工大学“油气藏地质与开发工程”国家重点实验室 ; ③西南石油大学)刘文革 ,贺振华 , 黄 德 济 , 杜 增 利 . 高 精 度 傅 里 叶 有 限 差 分 法 模 型 正 演 . 石 油 地 球 物 理 勘 探 , 2007 , 42( 6 ) : 629～633摘要 波动方程数值模拟方法分为有限差分法和频率 —波数域法两类 ,其中有限差分法的计算精度取决于波场外推算子的近似程度 、离散网格间距及差分方程阶数 ,它能适应速度任意横向变化 ,但在大倾角处易出现频 散现象及背景噪声 。

频率 —波数域法算法简单 、精度高 、噪声小 ,能适应任意地层倾角情况 ,但不适于速度场的 任意横向变化 。

文中结合有限差分法和频率 —波数域法的优点 ,应用傅里叶有限差分法 ( F F D) 实现在多域用高 精度延拓算子对模型进行地震记录的数值模拟 ,其波场外推算子由相移项 、折射项 ( 时移项) 和有限差分补偿项 组成 。

对 F F D 法进行了理论与误差分析 ,并用单程声波方程分别进行了层状模型和 SE G / EA GE 盐丘模型的数 值模拟试验 。

数值试验的对比分析表明 , F F D 法适用于速度场横向剧烈变化情形 ,且具有精度高 、无频散 、背景 噪声弱等优点 ,模拟结果反射特征清楚 ,能对复杂地质构造进行准确的地震数值模拟 。

关键词 波动方程 地震数值模拟 有限差分法 频率 —波数域法 傅里叶有限差分法提高微分方程对单平方根算子的逼近程度 ,改善了 成像效果[ 3 ]; C lae r b o u t 在研究时 —空域和频率 —空间域中的有限差分法后 ,认为有限差分法可以适应横向速度的任意变化 ,但是单平方根算子展开引起 了高波数近似 ,使得对大倾角构造的成像误差较大 , 易出现频散现象及背景噪声[ 4 ] 。

有限差分法地震正演模拟程序一．二阶公式推导1．二维的弹性波动方程22222222x x x ZU U U U t x z x z λμμλμρρρ∂∂∂∂++=++∂∂∂∂∂ 22222222xz z z U U U U t z x x zλμμλμρρρ∂∂∂∂++=++∂∂∂∂∂ 2．对方程进行中间差分（1）首先对时间进行二阶差分()()()()()112,,,22222n n n xi k xi k xi kx x x x U U U U t t U t U t t U t t t +--++-+-∂==∂ ()()()()()112,,,22222n n n zi k zi k zi kz z z z U U U U t t U t U t t U t t t +--++-+-∂==∂ （2）对x 方向进行二阶差分()()()()()21,,1,22222n n n xi k xi k xi kx x x x U U U U x x U x U x x U x x x +--++-+-∂==∂ ()()()()()21,,1,22222n n n zi k zi k zi kz z z z U U U U x x U x U x x U x x x +--++-+-∂==∂ （3）对z 方向进行二阶差分()()()()()2,1,,122222n n nxi k xi k xi k x x x x U U U U z z U z U z z U z z z +--++-+-∂==∂ ()()()()()2,1,,122222n n nzi k zi k zi k z z z z U U U U z z U z U z z U z z z +--++-+-∂==∂ （4）对xz 进行差分2,1,11,11,11,11,111,11,11,11,122224n nzi k zi k z zn n n nzi k zi k zi k zi k n n n n zi k zi k zi k zi k U U U U x z x zx z U U U U U U U U x x z z x+-++-++--++-++---⎛⎫-∂∂∂∂⎛⎫== ⎪ ⎪ ⎪∂∂∂∂∂⎝⎭⎝⎭-----+==2,1,11,11,11,11,111,11,11,11,122224n nxi k xi k x x n n n n xi k xi k xi k xi k n n n n xi k xi k xi k xi k U U U U x z x z x z U U U U U U U U x x z z x+-++-++--++-++---⎛⎫-∂∂∂∂⎛⎫== ⎪ ⎪ ⎪∂∂∂∂∂⎝⎭⎝⎭-----+==（5）把（1）（2）（3）（4）得到的结果带入波动方程3．写出差分方程()()()11,,,1,,1,22,1,,11,11,11,11,12222=2++4n n n n n nxi k xi k xi kxi k xi k xi k nnnnnnnxi k xi k xi k zi k zi k zi k zi k U U U U U U t x U U UU UUUz xz λμρμλμρρ+-+-+-++-++----+-++-+--++()()()11,,,1,,1,22,1,,11,11,11,11,12222=2++4n n n nnnzi k zi k zi kzi k zi k zi knn n nnnnzi k zi k zi k xi k xi k xi k xi k U U U U U U t x U UUU UUUz xz λμρμλμρρ+-+-+-++-++----+-++-+--++即得到()()1,,1,,1,,122112,,,1,11,11,11,1222+=2-++4n n n n n nxi k xi k xi kxi k xi k xi k n n n xi k xi k xi k n n n nzi k zi k zi k zi k U U U U U U x z U U U t U U U U z x λμμρρλμρ+-+-+-++-++---⎛⎫-+-++ ⎪ ⎪ ⎪--++ ⎪ ⎪⎝⎭()()1,,1,,1,,122112,,,1,11,11,11,1222+=2-++4n n n n n nzi k zi k zi k zi k zi k zi k n n n zi k zi k zi k n n n nxi k xi k xi k xi k U U U U U U x z U U U t U U U U z x λμμρρλμρ+-+-+-++-++---⎛⎫-+-++ ⎪ ⎪ ⎪--++ ⎪ ⎪⎝⎭二．MATLAB 程序 clear; clc;Nx=200; Nz=300;fi=30;%%%主频t_step=300;%%%%时间采样点dx=10.0;%空间采样间隔——x 方向 dz=10.0;%空间采样间隔——z 方向 dt=0.001;%时间采样间隔——1mslambda=66*1e9;%砂岩拉梅常数lamdamu=44*1e9;%砂岩拉梅常数murho=2650;%砂岩密度%%%%%%模型扩展%%%%%vp=zeros(Nz+2,Nx+2);%纵波速度vs=zeros(Nz+2,Nx+2);%横波速度c=zeros(Nz+2,Nx+2);%交叉项系数for i=1:Nz+2for j=1:Nx+2vp(i,j)=sqrt((lambda+2*mu)/rho);vs(i,j)=sqrt((mu)/rho);c(i,j)=sqrt((lambda+mu)/rho);endend%%%%%% 震源%%%%%%%t_wavelet=[1:t_step]*dt-1.0/fi;source=((1-2*(pi*fi*t_wavelet).^2).*exp(-(pi*fi*t_wavelet).^2));% 雷克子波amp=sqrt(2.0);% 振幅source_x=floor(Nx/2)+1;% 震源位置——x坐标source_z=2;% 震源位置——z坐标source_amp=zeros(Nz+2,Nx+2);% 震源振幅初始化（所有点处均为0）source_amp(source_z,source_x)=amp;% 震源振幅，在位置source_z,source_x处振幅为amp，其它位置为0 %%%%%%%%%%%%%%%%%%%%%%%%%%%U=zeros(Nz,Nx);% 弹性波x分量V=zeros(Nz,Nx);% 弹性波z分量U0=zeros(Nz+2,Nx+2);% 前一时刻的UU1=zeros(Nz+2,Nx+2);% 当前时刻的UU2=zeros(Nz+2,Nx+2);% 下一时刻的UV0=zeros(Nz+2,Nx+2);% 前一时刻的VV1=zeros(Nz+2,Nx+2);% 当前时刻的VV2=zeros(Nz+2,Nx+2);% 下一时刻的Vrecord_u=zeros(t_step,Nx);% x方向接收到的地震记录——Urecord_v=zeros(t_step,Nx);% x方向接收到的地震记录——V%%%%%% 有限差分计算U V %%%%%%for it=1:t_stepfor x=2:Nx+1for z=2:Nz+1U2(z,x)=2*U1(z,x)-U0(z,x)+(dt*dt)*(vp(z,x)^2*(1.0/(dx*dx))*(U1(z,x+1)-2*U1(z,x) +U1(z,x-1))+vs(z,x)^2*(1.0/(dz*dz))*(U1(z+1,x)-2*U1(z,x)+U1(z-1,x))+c(z,x)^2*(1. 0/(4*dx*dz))*(V1(z+1,x+1)-V1(z+1,x-1)-V1(z-1,x+1)+V1(z-1,x-1)))+source(it)*sour ce_amp(z,x)*dt*dt;V2(z,x)=2*V1(z,x)-V0(z,x)+(dt*dt)*(vs(z,x)^2*(1.0/(dx*dx))*(V1(z,x+1)-2*V1(z,x) +V1(z,x-1))+vp(z,x)^2*(1.0/(dz*dz))*(V1(z+1,x)-2*V1(z,x)+V1(z-1,x))+c(z,x)^2*(1. 0/(4*dx*dz))*(U1(z+1,x+1)-U1(z+1,x-1)-U1(z-1,x+1)+U1(z-1,x-1)));endendU0=U1;U1=U2;V0=V1;V1=V2;record_u(it,:)=U2(2,[2:201]);record_v(it,:)=V2(2,[2:201]);U=U2(2:301,2:201);V=V2(2:301,2:201);end%%%%% 绘制U 和V 的波场图%%%%%%figureimagesc(U(1:300,1:200));title('U');xlabel('x');ylabel('z');figureimagesc(V(1:300,1:200));title('V');xlabel('x');ylabel('z');%%%%% 绘制U 和V 的地震记录%%%%%%figureimagesc(record_u(1:300,1:200));title('U');xlabel('x');ylabel('t_step');figureimagesc(record_v(1:300,1:200));title('V');xlabel('x');ylabel('t_step');U分量波场图V分量波场图U分量地震记录V分量地震记录。

波动方程有限差分一、引言波动方程是自然界中许多现象的数学模型，如声波、地震波等。

为了解决波动方程的数值解，有限差分方法是一种常用的数值计算方法。

本文将详细介绍波动方程有限差分的原理、方法和应用。

二、波动方程波动方程描述了介质中物理量随时间和空间变化的规律。

具体来说，假设介质中某个物理量为u(x, t)，其中x表示空间坐标，t表示时间，则波动方程可以表示为：∂²u/∂t² = c²∇²u其中c表示介质中的传播速度，∇²表示拉普拉斯算子。

该方程描述了一个在介质中传播的二阶偏微分方程。

三、有限差分方法有限差分方法是一种常用的数值计算方法，其基本思想是将连续函数离散化为离散点上的函数值，并通过差商逼近导数或偏导数，从而得到原问题的近似解。

对于波动方程，在空间上进行网格剖分，并在每个网格点处离散化u(x, t)和其导数，可以得到如下形式的差分格式：(u(i, j+1) - 2u(i, j) + u(i, j-1)) / Δt² = c²((u(i+1, j) - 2u(i, j) + u(i-1, j)) / Δx² + (u(i, j+1) - 2u(i, j) + u(i, j-1)) / Δy²)其中i表示空间网格点的横坐标，j表示纵坐标，Δt、Δx和Δy分别为时间和空间上的步长。

这个差分方程可以通过迭代求解得到波动方程的数值解。

具体来说，可以使用显式差分法或隐式差分法进行求解。

四、应用波动方程有限差分方法在地震勘探、声学建模等领域得到广泛应用。

例如，在地震勘探中，可以通过模拟地震波传播过程得到地下岩层的结构信息；在声学建模中，可以计算音场传播过程，并预测噪声污染等问题。

五、总结本文介绍了波动方程有限差分方法的原理、方法和应用。

有限差分方法是一种常用的数值计算方法，在许多领域都有广泛应用。

对于波动方程这类偏微分方程，有限差分方法是一种有效的求解方法。

面向波场分析的波动方程有限差分正演方法随着勘探对象的日益复杂化,使得研究的模型介质也越来越复杂,为了更详细地了解和认识复杂介质波场,进行地震波场数值模拟是非常必要的。

波动方程有限差分法是近几年最流行和最受欢迎的数值模拟方法之一,其优点就是模拟波场信息丰富、实现简单和运算速度较快等。

波动方程有限差分法虽然有很多优点,但是也存在一些不足之处,其中稳定性条件就是一个很重要且急需解决的问题。

稳定性条件直接决定正演模拟结果的成败,在粘弹介质波动方程的有限差分解中,介质的粘弹特性会增加数值解的不稳定性。

本文从粘弹性理论入手,研究了Kelvin-Voigt和Maxwell在不同精度差分格式下的稳定性条件,而对于标准线性粘弹模型,只给出了稳定性条件与速度和品质因子的数值解。

当空间网格步长一定时正演模拟中的时间步长是由模型中的最大速度决定的,速度越大时间步长越小,所需的运算量也越大,为了在不影响结果的前提下提高波动方程有限差分法模拟的效率,本文提出了波动方程变时间步长有限差分数值模拟。

近地表的主要特点就是低速度、低密度和低品质因子。

随着山地、沙漠等地表复杂区地震勘探的发展,近地表地震数值模拟技术受到了地球物理学家的广泛关注。

结合上述得到的结果进行地震波场近地表效应的模拟及分析,对这一部分的模拟采用自由边界条件,对近地表的瑞雷面波进行了模拟和定量分析,确定模拟出的是瑞雷面波。

粘弹性介质是更符合实际介质,进行了含有低速层的粘弹性介质在自由边界条件下的数值模拟。

还从Q值大小、震源子波频率等方面定量讨论粘弹介质地震波衰减影响因素。

实际介质(和油气藏有关的介质)往往是非均匀的,非均匀介质中存在很多微小异常且分布极不规则(通过测井数据和岩心样品就可以发现)。

对于上述情况实现了随机介质的构造和模拟,得到非均匀介质中的波场。

最后对几个实际速度模型进行数值模拟,得到了复杂储层(含有孔洞、裂缝等)下的全波场和储层在不同含水饱和度下的波场,并且通过频谱和时频分析对结果进行了分析,还研究了低速层和近地表对波场模拟的影响。

三维波动方程有限差分正演方法摘要：本文主要介绍了三维波动方程有限差分正演方法的基本原理和实现过程，并对其优缺点进行了分析。

通过数值实验验证了该方法的可行性和准确性，为地震勘探、地下水文学等领域的研究提供了参考。

关键词：三维波动方程、有限差分、正演方法、数值实验一、引言三维波动方程是地震勘探、地下水文学等领域的基础理论，其求解方法对于相关领域的研究具有重要意义。

有限差分正演方法是求解三维波动方程常用的数值方法之一，其具有计算速度快、精度高等优点，在实际应用中得到了广泛的应用。

本文将介绍三维波动方程有限差分正演方法的基本原理和实现过程，并对其优缺点进行分析，同时通过数值实验验证该方法的可行性和准确性。

二、三维波动方程有限差分正演方法的基本原理三维波动方程可以表示为：其中，u(x,y,z,t)为波动场，c(x,y,z)为介质速度，ρ(x,y,z)为介质密度，t为时间。

有限差分正演方法通过将空间和时间离散化，将三维波动方程转化为差分方程，进而求解波动场在不同时刻的数值解。

具体而言，有限差分正演方法将空间和时间分别离散化，将空间网格点和时间网格点相结合，构成一个四维网格空间，其中每个网格点对应一个波动场的数值解。

有限差分正演方法的求解过程可以分为以下几个步骤：1. 空间离散化对于三维空间，可以将其分为x、y、z三个方向，分别进行离散化。

假设在x方向上，空间步长为Δx，在y方向上，空间步长为Δy，在z方向上，空间步长为Δz。

则可以将空间网格点表示为(xi,yj,zk)，其中i=1,2,...,nx，j=1,2,...,ny，k=1,2,...,nz，nx、ny、nz分别为x、y、z方向的网格数。

2. 时间离散化假设时间步长为Δt，则在t时刻波动场的数值解可以表示为ui,j,k^n，其中n表示时间步数，i、j、k表示空间网格点的索引。

3. 有限差分近似将三维波动方程中的导数项用有限差分近似表示，例如：其中，ui,j,k^n表示在t时刻，在(xi,yj,zk)处的波动场数值解，c(xi,yj,zk)表示在(xi,yj,zk)处的介质速度，ρ(xi,yj,zk)表示在(xi,yj,zk)处的介质密度，Δx、Δy、Δz、Δt分别为空间和时间步长。

波动方程正演模型及应用吴清岭　张　平　施泽龙3(大庆石油管理局勘探开发研究院)摘　要　地震资料解释经常用到正演模型。

常规的褶积模型不能模拟地震波的动力学特征。

本文采用声波方程,通过四阶有限差分近似,实现了复杂地质构造零炮检距的数值模拟。

文中同时展示了实际应用效果。

　主题词　正演模型 有限差分 零炮检距剖面作者简介　吴清岭,男,1962年生,1983年毕业于华东石油学院勘探系,硕士,高级工程师,现从事地震方法研究工作。

地址:(163712)黑龙江省大庆市让胡路区勘探开发研究院。

3　参加本工作的还有杨有林同志。

在地震资料解释中,人们力图得到能够保持地震波的运动学与动力学特征的波动方程正演模型,以达到精确模拟地震波传播特性的目的。

在求解波动方程的2种数值解法(有限差分法和有限元法)中,有限差分法是一种快速有效的方法,并且地质模型的复杂程度不影响运算速度。

本文介绍了对声波方程采用四阶有限差分近似制作零炮检距剖面的基本过程及应用效果。

一、基本原理1,计算公式在二维空间域内,二维声波方程为1C 292u 9t 2=92u 9x 2+92u9z 2式中　C ———声学介质下地震波的纵波速度;u ———声压。

设Δh 为空间采样步长;Δt 为时间采样步长;m 、n 、l 分别为正整数;则有x =m ・Δh z =n ・Δh t =l ・Δt 对时间域采用二阶有限差分;对空间域采用四阶有限差分(推导过程略),其数值计算公式为u (m ,n ,l +1)=(A 2/12){16[u (m +1,n ,l )+u (m -1,n ,l )+u (m ,n +1,l )+u (m ,n -1,l )]-[u (m +2,n ,l )+u (m -2,n ,l )+u (m ,n +2,l )+u (m ,n -2,l )]}+(2-5A 2)[u (m ,n ,l )-u (m ,n ,l -1)]其中 A 2=C 2(m ,n )Δt 2/Δh 2式中　C (m ,n )———介质速度的空间离散值;Δt ———时间离散步长;Δh ———空间离散步长。