一阶弹性波方程的变网格高阶有限差分数值模拟.docx

1、一阶弹性波方程的变网格高阶有限差分数值模拟一阶弹性波方程的变网格高阶有限差分数值模拟 2019年12 月第43卷第6期 正演技术 一阶弹性波方程的变网格高阶有限差分数值模拟 李振春3张慧张华 (中国石油大学地球资源与信息学院, 山东东营257061) 李振春, 张慧, 张华. 一阶弹性波方程的变网格高阶有限差分数值模拟. 石油地球物理勘探, 2019, 43(6) : 711716 摘要使用可变网格的有限差分法进行地震模拟有许多独特的优点, 主要表现为对地质模型的离散化更为合理, 在低速带和复杂构造区域, 可将局部网格划分得相对精细些, 不仅提高了模拟精度, 消除了因采样不足导致的频散现象,

2、而且可以减少计算机内存需求, 保持模型计算的灵活性。本文提出一种新的基于高阶交错网格技术的弹性波数值模拟方法, 通过改变网格的空间步长实现了局部网格加密技术, 弥补了常规网格的缺陷和不足。试算结果表明, 本文提出的算法稳定性较好, 且能够提高模拟精度, 减少计算时间, 提高计算效率。关键词变网格交错网格高阶有限差分弹性波 1引言 法技术(地震采集。为适应复杂勘探条件的需要, 、高效率的地震数值模拟技术变得越来越重要。 自Virieux 1率先采用交错网格有限差分法模拟P 2SV 波在非均匀介质中的传播以来, 交错网格成为了有限差分法正演模拟中普遍使用的网格剖分技术。同二阶位移方程相比, 使用交

3、错网格的一阶速度应力方程, 无论是模拟精度还是模拟的稳定性都有了较大提高。但是如果局部地区构造较为复杂, 或者存在低速层时, 为了提高复杂构造的正演精度, 通常对整个模型使用较小的采样网格, 以保证对频散现象的抑制和对稳定性的保障。这对保持模型计算的灵活性非常不利, 还会对模型的其他部分造成过采样, 增加惊人的计算量。所以, 设计一种更合理的离散化方式进行正演模拟是非常必要的。 从空间采样的角度考虑, 既要提高模拟精度、又要降低计算机内存需求的最有效方法就是在模型的不同区域采用不同步长的网格, 即变网格(也称为不规则网格) 。为此J ast ram 等提出了基于二维声波方程针对某一深度变网格步

4、长的算法2;J ast ram 等 3; 4Pitarka 实现了低阶精度的不5。 本文在前人研究的基础上提出一种高阶变网格数值模拟方法, 实现了高精度的局部网格加密技术, 另外还改进了前人对变化网格过渡带采用的波场插值算法, 避免了在过渡带区域进行大量的插值计算, 也不需要平滑函数进行平滑, 同时还保持了算法的稳定性。 2方法原理 文中所提出的算法是基于Virieux 1和Graves 6 提出的一阶速度应力方程, 对于各向同性介质, 有 xx +55t v x =b (5x z xz ) +f x xz +55t v z =b (5x z zz ) +f z +2) 5x v x +55z

5、 v z t xx =( +2) 5z v z +55x v x t zz =( (5z v x +5x v z ) 5t xz = (1) 式中:v x , v z 分别表示速度; xx , zz , xz 分别表示应力向量; 5x 、5z 、5t 分别表示一阶偏微分算子5/5x 、5/5z 、5/5t ; b 表示密度的倒数; f x , f z 为体力向量; 及为拉梅常数。 3山东省东营市中国石油大学地球资源与信息学院,257061 本文于2019年1月8日收到, 修改稿于同年7月30日收到。本课题由国家973专题(2019CB209605) 、国家863专题(2019AA06Z206)

6、 和CNPC 物探重点实验室中国石油大学(华东) 研究室资助。 712石油地球物理勘探2019年 通过对时间导数采用二阶近似, 得到弹性波一阶速度应力方程的离散形式为 v x v n+2 i , j =v x n-2 i , j xx +D +tb x (D x z xz ) |n i , j n 1=d x i /2, 2=d x i /2 n-1 n-1k =1 n+2 z i+2, j+2 =v n-2 z i+2, j+2 +tb z (D x xz +D z zz ) |i+, j+22 n+2i+2, j n+2i+2, j n-1= k =1 d x i+k+d x i /2,

7、n =d x i-k +d x i /2 n+1 xx i +2, j n+1 zz i +2, j = n xx i +2, j n zz i +2, j +2) D x v x +t(D z v z 图1计算对称点为i + 的差分系数所需的变网格节点2 n-1 =) D z v z +t(+2D x v x 同理, 以i 为对称点的差分算子为n-1= n-1k =1 n+1n xz (D z v x +D x v z ) xz i , j +2=xz i , j +2+t( n+2 i , j +2 d x i+(k-1) +d x i+n-1/2, n = k =1 d x i-k +d

8、 x i-n /2。 (2) 式中:D x , D z 分别表示对x , z 的一阶微分算子; b x =(b i , j +b i +1, j ) /2, b z =(b i , j +b i , j +1) /2; xz =(1/i , j +1/i +1, j +1/i , j +1+1/i +1, j +1) /4。 用变量g 来表示速度v x 、v z 或应力向量xx 、zz 、xz , 则在常步长的交错网格技术中, 用下式计算式 (2) 中的g 关于x 的2n 阶空间偏导数7 D x g (x , z ) = C g (x +x x n i i =1 下面给出变网格算法的任意偶数阶

9、精度差分近似式及差分系数的求取方法, 即 n D x g (x , z ) = i -1 c 2i-1 g (x +2i-1, z ) + (4) +c 2i -2i , z ) n , ) - i ; i 是空间差分算子, 它 d g 关于x 。令g (x , z ) =g z exp (i k x ) , 则式(4) 可写成 n (3) -g (x -x i-1z x 2n 式中, 系数C i n 由x (x 2i -1, z ) 和g (x -x 2i -1, z ) 。 与常网格不同, 变网格的差分算子是空间变化 i k = i -1 c 2i-1 2i-1) +c 2i exp (-

10、i k 2i ) exp (+i k (5) 的。给定地质模型, 在确定网格剖分后就可求出各 网格节点的差分算子供递推计算时调用。在交错网 格技术中差分算子有两个对称点, 即i +或i , 不 2 同的对称点对应的i 也不同。下面以对称点i + 2 为例求取i , 如图1所示。 对式(5) 中的指数项进行泰勒展开, 它的2n 阶泰勒展开式为 2 i ) 1+i k i +(i k ) 2exp (i k i + 2 2n-12n (k i ) 2n-1+o (6) +i ) (2n -1) ! 把式(6) 代入式(5) 并整理, 可以得到 2222 ) () (c 12i k =(c 1+c

11、2+c 2n ) +i k (c 1-c +-c +i 1222n 2n 1+c 22+c 2n 2n ) + 2 +(i ) 3 3 6 332n-1 (c 131+c 22+-c 2n 2n ) +(i ) 2 n-12n-12n-1 (c 1-c 2+-c 2 n 2n ) 12 (2n -1) ! (7) 2n-1 写成矩阵的形式 1i 1 1 -i 2 1i 2n-1 1 -i 2n c 1c 2= 2n-22n-1 0i (i ) 2n-21 2n-2 2n-1 (i ) 2n-11 2n-2 (i ) 2n-222n-12n-1 -(i ) 2 (i ) 2n-22n 2n-1

12、-(i ) 2n 2n-2 (i ) 2n-22n-12n-1(i ) 2n-12n-1 c 2n-1c 2n 00 (8) 第43卷第6期李振春等:一阶弹性波方程的变网格高阶有限差分数值模拟713 通过解方程组(8) , 可以得到有限差分算子D x 的系数c i 。同样地, 差分算子D z 的差分系数也是通过这种方法求得。对于交错的网格, 由于变量所定义的网格点的位置不同, 对于每一个差分算子D x 或D z 都有两组不同的差分系数。 在使用变网格时, 网格步长的变化可能会导致数值反射的出现。究其原因是因为波场离散后, 相速度是网格步长的函数。当相速度梯度较大时, 即使速度和密度都没有变化,

13、 入射波的能量也会部分反射回来, 导致了数值反射现象。因为频散程度与频率有关, 所以当入射波的频率较低时, 产生的数值反射较高频时弱。在网格变化的区域对波场进行插值计算虽可以在某种程度上压制数值反射, 但是插值算法计算量大, 且效果并不令人满意, 所以本文没有采用插值计算。本文提出的算法是根据各点所对应的步长的变化计算其差分系数, 在一个循环内完成所有点的波场计算, 不会产生数值反射, 并且减小了计算量。 式中:h =max (d x i , d z i ) 表示最大的空间网格步长; v min 是地震波速的最小值。 本文所提算法需要满足的稳定性条件为 min (d x i , d z i )

14、 n i =1 (-1) i-1C i 式中:V max 为最大速度值; 为密度; C i 为式(3) 中做2n 阶空间展开时的系数。 4数值模拟实例 为了验证提出的变网格差分算法的稳定性及有效性, 设计两个模型正演实例来说明文中所述方法的实际效果。4. 1。模3000m , 纵波速度为3500m/s , 横0125来计算, 密度为常数。P 波震源激发, 震源主频为30Hz , 位于(1500m ,1500m ) 处。变网格步长为10m ,5m 。与其进行对比的常规交错网格步长为5m 。数值模拟中, 差分解法的精度为o (t 2, x 10 ) 。弹性波波场快照和单炮记录分别记录于图2和图3中。 从上面的弹性波波场快照和单炮记录中可以看出, 变网格的过渡区没有产生数值边界反射, 说明本文提出的算法成功地消除了空间网格变化造成的数值反射, 进而说明了该算法的合理性和良好的稳定性。 3法较好地解决了边界反射问题。使用透射边界能够吸收绝大部分边界反射的能量, 剩余能量通过衰减边界再次减小, 这样计算量增加不大, 吸收也较彻底, 可以较好地消除人为边界反射。 稳定性条件是有限差分数值模拟