一阶弹性波方程的变网格高阶有限差分数值模拟Word格式.docx
《一阶弹性波方程的变网格高阶有限差分数值模拟Word格式.docx》由会员分享，可在线阅读，更多相关《一阶弹性波方程的变网格高阶有限差分数值模拟Word格式.docx（13页珍藏版）》请在冰豆网上搜索。
这对保持模型计算的灵活性非常不利,还会对模型的其他部分造成过采样,增加惊人的计算量。
所以,设计一种更合理的离散化方式进行正演模拟是非常必要的。
从空间采样的角度考虑,既要提高模拟精度、又要降低计算机内存需求的最有效方法就是在模型的不同区域采用不同步长的网格,即变网格（也称为不规则网格）。
为此Jastram等提出了基于二维声波方程针对某一深度变网格步长的算法[2];
Jastram等
[3];
[4Pitarka实现了低阶精度的不[5]。
本文在前人研究的基础上提出一种高阶变网格数值模拟方法,实现了高精度的局部网格加密技术,另外还改进了前人对变化网格过渡带采用的波场插值算法,避免了在过渡带区域进行大量的插值计算,也不需要平滑函数进行平滑,同时还保持了算法的稳定性。
2 方法原理
文中所提出的算法是基于Virieux[1]和Graves[6]
提出的一阶速度—应力方程,对于各向同性介质,有
τxx+5τ5tvx=b（5xzxz）+fx
τxz+5τ5tvz=b（5xzzz）+fzλ+2μ）5xvx+λ5τ5zvztxx=（
λ+2μ）5zvz+λ5τ5xvxtzz=（
（5zvx+5xvz）5τtxz=μ
（1）
ττ式中:
vx,vz分别表示速度;
τxx,zz,xz分别表示应力向量;
5x、5z、5t分别表示一阶偏微分算子5/5x、5/5z、5/5t;
b表示密度的倒数;
fx,fz为体力向量;
λ及μ为拉梅常数。
3山东省东营市中国石油大学地球资源与信息学院,257061
本文于2019年1月8日收到,修改稿于同年7月30日收到。
本课题由国家973专题（2019CB209605）、国家863专题（2019AA06Z206）和CNPC物探重点实验室中国石油大学（华东）研究室资助。
712 石油地球物理勘探2019年
通过对时间导数采用二阶近似,得到弹性波一阶速度—应力方程的离散形式为
vxv
n+2
i,j
=vx
n-2
τxx+Dτ+Δtbx（Dxzxz）|n
i,j
n
Δ 1=dxi/2,Δ2=dxi/2
n-1
n-1k=1
z i+2,j+2
=v
τ+Δtbz（Dxxz+Dτzzz）|i+,j+22
n+2i+2,jn+2i+2,j
Δn-1=
k=1
dxi+k+dxi/2,Δn=∑dxi-k+dxi/2∑
ττ
n+1
xxi+2,jn+1
zzi+2,j
=τ
xxi+2,jn
λ+2μ）Dxvx+λ+Δt[（Dzvz]
图1 计算对称点为i+
的差分系数所需的变网格节点2
=τ）Dzvz+λ+Δt[（λ+2μDxvx] 同理,以i为对称点的差分算子为Δn-1=
n+1n
τμxz（Dzvx+Dxvz）]xzi,j+2=τxzi,j+2+Δt[（
i,j+2
dx∑
i+（k-1）
+dxi+n-1/2,Δn=
i-k
+dxi-n/2。
（2）
式中:
Dx,Dz分别表示对x,z的一阶微分算子;
bx=（bi,j+bi+1,j）/2,bz=（bi,j+bi,j+1）/2;
μxz=（1/μi,j+1/μi+1,j+1/μi,j+1+1/μi+1,j+1）/4。
用变量g来表示速度vx、vz或应力向量τxx、ττzz、xz,则在常步长的交错网格技术中,用下式计算式
（2）中的g关于x的2n阶空间偏导数[7]
Dxg（x,z）=
C[g（x+Δx
Δx∑
ni
i=1
下面给出变网格算法的任意偶数阶精度差分近似式及差分系数的求取方法,即
i-1
[c∑
2i-1
g（x+Δ2i-1,z）+
（4）
+c2i-Δ2i,z）]
）-
i;
i是空间差分算子,它
dg关于x。
令g（x,z）=gzexp（ikx）,则式（4）可写成
（3） -g（x-Δxi-1zx2n式中,系数Cin由x（x2i-1,z）和g（x-Δx2i-1,z）。
与常网格不同,变网格的差分算子是空间变化
ik=
Δ2i-1）+c2iexp（-ikΔ2i）]exp（+ik
（5）
的。
给定地质模型,在确定网格剖分后就可求出各
网格节点的差分算子供递推计算时调用。
在交错网
格技术中差分算子有两个对称点,即i+或i,不
2
同的对称点对应的Δi也不同。
下面以对称点i+
为例求取Δi,如图1所示。
对式（5）中的指数项进行泰勒展开,它的2n阶泰勒展开式为
Δi）≈1+ikΔi+（ik）2Δexp（iki+…+
2n-12n
（kΔi）2n-1+o（Δ（6） +i）
（2n-1）!
把式（6）代入式（5）并整理,可以得到
2222
ΔΔ）（）（c1ΔΔ2 ik=（c1+c2+…+c2n）+ik（c1Δ-c+…-c+i1222n2n1+c2Δ2+…+c2n2n）+
+（i）3
3
6
332n-1
（c1ΔΔ31+c2Δ2+…-c2n2n）+…+（i）
n-12n-12n-1
（c1Δ-c2Δ+…-c2
nΔ2n）12
（7）
2n-1
写成矩阵的形式
1iΔ1
1
-iΔ2
…
…ω
1iΔ2n-1
-iΔ2n
c1 c2 …=
2n-22n-1
0i
（i）2n-2Δ1
2n-2
（i）2n-1Δ1
（i）2n-2Δ22n-12n-1
Δ-（i）2
……
（i）2n-2Δ2n
Δ-（i）2n
…（i）2n-2Δ2n-12n-1…（i）2n-1Δ2n-1
c2n-1c2n
00
（8）
第43卷 第6期李振春等:
一阶弹性波方程的变网格高阶有限差分数值模拟 713
通过解方程组（8）,可以得到有限差分算子Dx的系数ci。
同样地,差分算子Dz的差分系数也是通过这种方法求得。
对于交错的网格,由于变量所定义的网格点的位置不同,对于每一个差分算子Dx或Dz都有两组不同的差分系数。
在使用变网格时,网格步长的变化可能会导致数值反射的出现。
究其原因是因为波场离散后,相速度是网格步长的函数。
当相速度梯度较大时,即使速度和密度都没有变化,入射波的能量也会部分反射回来,导致了数值反射现象。
因为频散程度与频率有关,所以当入射波的频率较低时,产生的数值反射较高频时弱。
在网格变化的区域对波场进行插值计算虽可以在某种程度上压制数值反射,但是插值算法计算量大,且效果并不令人满意,所以本文没有采用插值计算。
本文提出的算法是根据各点所对应的步长的变化计算其差分系数,在一个循环内完成所有点的波场计算,不会产生数值反射,并且减小了计算量。
h=max（dxi,dzi）表示最大的空间网格步长;
vmin是地震波速的最小值。
本文所提算法需要满足的稳定性条件为
≤
min（dxi,dzi）ni=1
∑
（-1）i-1Ci
Vmax为最大速度值;
ρ为密度;
Ci为式（3）中做2n阶空间展开时的系数。
4 数值模拟实例
为了验证提出的变网格差分算法的稳定性及有效性,设计两个模型正演实例来说明文中所述方法的实际效果。
4.1 。
模×
3000m,纵波速度为3500m/s,横0125来计算,密度为常数。
P波震源激发,震源主频为30Hz,位于（1500m,1500m）处。
变网格步长为10m,5m。
与其进行对比的常规交错网格步长为5m。
数值模拟中,差分解法的精度为o（Δt2,Δx10
）。
弹性波波场快照和单炮记录分别记录于图2和图3中。
从上面的弹性波波场快照和单炮记录中可以看出,变网格的过渡区没有产生数值边界反射,说明本文提出的算法成功地消除了空间网格变化造成的数值反射,进而说明了该算法的合理性和良好的稳定性。
3 法较好地解决了边界反射问题。
使用透射边界能够吸收绝大部分边界反射的能量,剩余能量通过衰减边界再次减小,这样计算量增加不大,吸收也较彻底,可以较好地消除人为边界反射。
稳定性条件是有限差分数值模拟
升级会员 