旋转交错网格TTI介质波场模拟与波场分解

引用本文

李长江, 李庆春, 王亚妮. 旋转交错网格TTI介质波场模拟与波场分解[J]. 物探与化探, 2015,39(3): 553-557.
LI Chang-Jiang, LI Qing-Chun, WANG Ya-Ni. The wave-field simulation and decomposition in the TTI media by the rotated staggered grid[J]. Geophysical and Geochemical Exploration, 2015,39(3): 553-557.
Doi:10.11720/wtyht.2015.3.19 复制到剪切板

Permissions

《物探与化探》编辑部

旋转交错网格TTI介质波场模拟与波场分解

李长江, 李庆春, 王亚妮

长安大学地质工程与测绘学院,陕西西安 710054

作者简介: 李长江(1987-),男,在读研究生,长安大学地质工程与测绘学院,从事地震数据的波场模拟及分离分解。

收稿日期: 2014-11-08

基金: 国家自然科学基金项目(41374145,41304105); 教育部博士点基金项目(20120205130002)

摘要

为了得到TTI介质中传播的纯纵波和纯横波波场,采用旋转交错网格高阶有限差分法对TTI介质中传播的地震波场进行了数值模拟,并利用qP波、qSV波、SH波的极化方向相互垂直的特性,来实现qP波、qSV波、SH波的分解。通过对均匀介质模型和复杂层状介质模型分解前后波场快照和地震记录的对比,验证了本文方法的正确性。

关键词: 旋转交错网格; TTI介质; 有限差分; 波场模拟; 波场分解

中图分类号:P631.4 文献标志码:A 文章编号:1000-8918(2015)03-0553-05

The wave-field simulation and decomposition in the TTI media by the rotated staggered grid

LI Chang-Jiang, LI Qing-Chun, WANG Ya-Ni

School of Geological Engineering and Geomatics,Chang'an University,Xi'an 710054,China

Abstract

In order to obtain the pure P and SV and SH waves in the TTI media.The rotated staggered grid finite-difference method was used to simulate the waves propagate in the TTI media,and the properties of the P、SV、SH polarization directions perpendicular to each other was used to decompose the coupled waves.This method is testified by thecomparison of the original wave-field and the decomposed wave-fields of a homogeneous media and a complex 6 layers media.

Keyword: rotated staggered grid; TTI media; finite difference method; wave-field modeling; wave-field decomposition

Show Figures

完全弹性波动方程数值模拟中, 纵、横波耦合的问题是研究纯纵波和纯横波的传播规律以及利用纯纵、横波逆时偏移成像的难题, 因此波场分解和纯纵、横波波场模拟越来越受到重视^{[1, 2, 3, 4, 5]}。马德堂等^[6]、李振春等^[7]在散度旋度的理论上提出等价方程, 模拟了纯P波和纯S波的传播特征。Dellinger^{[8, 9]}提出通过将地震波场投影到纵、横波的极化向量上来分离各向异性介质中的纵波, 然后用全波场减去纵波波场来得到非纵波场。Yan 和 Sava^[10]将波数域中的波场分离算子进行反傅立叶变换, 得到空间域随空间变化的分离算子, 该算子能适应物性参数缓慢变化的复杂介质, 计算量大, 且精度会受到物性参数剧烈变化的影响。Qunshan Zhang^[11]引入了波场分解算子, 将VTI介质中的全波波场分解为P波、SV波波场, 该算子克服了散度、旋度波场分解算子引起的相移, 但其在得到最终的分解结果时, 采用了平均策略。杜启振^[12]指出平均策略降低了波场分解的精度, 所以采用在波数域对波场进行插值处理来提高精度。

在旋转交错网格中, 所有的速度分量被定义在一个位置, 所有的应力分量被定义在和速度分量相对的角点上^[13]。在进行各向异性介质波场模拟时, 不需要对弹性常数进行平滑或插值处理, 相较于普通交错网格, 旋转交错网格在TTI介质波场中的模拟精度更高。同理, 在进行波场分解时不需要采用任何平均策略和插值处理, 提高了分解精度和计算效率。

1 一阶速度— 应力方程

由本构方程和运动微分方程可以得到二维三分量TTI介质波动方程

$\begin{matrix} \begin{matrix} ρ \frac{\partial u}{\partial t} = \frac{\partial τ_{xx}}{\partial x} + \frac{\partial τ_{xz}}{\partial z}, ρ \frac{\partial v}{\partial t} = \frac{\partial τ_{xy}}{\partial x} + \frac{\partial τ_{yz}}{\partial z}, ρ \frac{\partial w}{\partial t} = \frac{\partial τ_{xz}}{\partial x} + \frac{\partial τ_{zz}}{\partial z} 。 (1) \\ \begin{matrix} \partial τ_{xx} / \partial t = C_{11} \partial u / \partial x + C_{13} \partial w / \partial z + C_{14} \partial v / \partial z + C_{15} (\partial u / \partial z + \partial w / \partial x) + C_{16} \partial v / \partial x, \\ \partial τ_{zz} / \partial t = C_{13} \partial u / \partial x + C_{33} \partial w / \partial z + C_{34} \partial v / \partial z + C_{35} (\partial u / \partial z + \partial w / \partial x) + C_{36} \partial v / \partial x, \\ \partial τ_{yz} / \partial t = C_{14} \partial u / \partial x + C_{34} \partial w / \partial z + C_{44} \partial v / \partial z + C_{45} (\partial u / \partial z + \partial w / \partial x) + C_{46} \partial v / \partial x, \\ \partial τ_{xz} / \partial t = C_{15} \partial u / \partial x + C_{35} \partial w / \partial z + C_{45} \partial v / \partial z + C_{55} (\partial u / \partial z + \partial w / \partial x) + C_{56} \partial v / \partial x, \\ \partial τ_{xy} / \partial t = C_{16} \partial u / \partial x + C_{36} \partial w / \partial z + C_{46} \partial v / \partial z + C_{56} (\partial u / \partial z + \partial w / \partial x) + C_{66} \partial v / \partial x 。 \end{matrix} (2) \end{matrix} \end{matrix}$

其中:u、v、w为质点的x、y、z速度分量, τ _xx、τ _zz、τ _yz、τ _xz、τ _xy为应力分量, c_ij为弹性系数, ρ 为介质密度。

TTI介质的弹性系数矩阵是通过Bond变换将VTI介质的弹性系数矩阵旋转到观测坐标系下得到的, 由于受到极化角和方位角的影响, TTI介质弹性矩阵几乎全不为零, 使得TTI介质的弹性波动方程变得复杂, 也更接近真实地球模型, 因此研究TTI介质中传播的地震波的传播规律具有重要的理论价值和实际意义。

2 TTI介质波场分解理论

在各向同性介质中可以利用P波和S波的极化方向相互垂直的特性来分解纵横波场, 同理可将各向异性介质中P波和S波的极化向量用于各向异性介质中P波和S波的波场分离与分解。设P波、SV波和SH波的极化向量分别为A_P、A_SV和A_SH, 通过求解Christoffel方程

$\begin{matrix} Γ - ρ V^{2} I = 0 (3) \end{matrix}$

来得到A_P、A_SV和A_SH。其中:Γ 是Christoffel矩阵; Γ _ij=c_ijkln_jn_l, c_ijkl为弹性张量, n_j、n_l为归一化之后的j、l方向的波数, (i, j, k, l)=1, 2, 3; V是相速度, 与地震波传播方向有关; I是单位向量。Christoffel方程的求解是一个标准的3× 3的特征值与特征向量问题, 三个特征值分别代表P波、SV波和SH波的相速度, 与之对应的特征向量分别代表P波、SV波和SH波的极化向量。

将Christoffel方程展开为

$\begin{matrix} [\begin{matrix} Γ_{11} - ρ V^{2} & Γ_{12} & Γ_{13} \\ Γ_{21} & Γ_{22} - ρ V^{2} & Γ_{23} \\ Γ_{31} & Γ_{32} & Γ_{33} - ρ V^{2} \end{matrix}] [\begin{matrix} A_{x} \\ A_{y} \\ A_{z} \end{matrix}] = 0, (4) \end{matrix}$

其中:Γ _ij是介质的弹性参数和地震波传播方向的函数。通过Christoffel方程可以求得P波、SV波和SH波的极化向量, 分别为:

P波, A_P=[A_P_x, A_P_y, A_P_z];

SV波, A_SV=[A_SV_x, A_SV_y, A_SV_z];

SH波, A_SH=[A_SH_x, A_SH_y, A_SH_z]。

各向异性介质中P波、SV波、SH波的波场分解公式为

$\begin{matrix} \begin{matrix} {\dot{W}}_{P} = A_{P} (A_{P} \cdot \dot{W}), \\ {\dot{W}}_{SV} = A_{SV} (A_{SV} \cdot \dot{W}), (5) \\ {\dot{W}}_{SH} = A_{SH} (A_{SH} \cdot \dot{W}) 。 \end{matrix} \end{matrix}$

对 $\begin{matrix} {\dot{W}}_{P} \end{matrix}$ 、 $\begin{matrix} {\dot{W}}_{SV} \end{matrix}$ 、 $\begin{matrix} {\dot{W}}_{SH} \end{matrix}$ 作三维或二维傅立叶反变换可得到各向异性介质中的波场分解之后的P波、SV波和SH波波场, $\begin{matrix} \dot{W} \end{matrix}$ 为质点速度分量在波数域中的表示。

3 旋转交错网格有限差分

旋转交错网格因其在模拟各向异性介质时不需要对弹性参数作任何特殊处理, 只需要对密度进行平均, 提高了TTI介质波场模拟的精度。图1为旋转交错网格和普通交错网格的对比。

	Figure Option View Download New Window
	图1 普通交错网格(左)和旋转交错网格(右)定义对比(据Bohlen T^[14]修改)

对密度进行平均

$\begin{matrix} \begin{matrix} \overset{̅}{ρ} (i + 1 / 2, j + 1 / 2) = 1 / 4 [ρ (i, j) + ρ (i + 1, j) + \\ ρ (i + 1, j + 1) + ρ (i, j + 1)] 。 (6) \end{matrix} \end{matrix}$

在二维情况下, 旋转交错网格有限差分中沿坐标轴的差分算子可以表示为

$\begin{matrix} \begin{matrix} \frac{\partial}{\partial \dot{z}} = \frac{Δx}{Δr} \frac{\partial}{\partial x} + \frac{Δz}{Δr} \frac{\partial}{\partial z}, \\ \frac{\partial}{\partial \dot{x}} = \frac{Δx}{Δr} \frac{\partial}{\partial x} - \frac{Δz}{Δr} \frac{\partial}{\partial z} 。 \end{matrix} (7) \end{matrix}$

通过求解式(7)可得

$\begin{matrix} \begin{matrix} \frac{\partial}{\partial x} = \frac{Δr}{2 Δx} (\frac{\partial}{\partial \dot{z}} + \frac{\partial}{\partial \dot{x}}), \\ \frac{\partial}{\partial z} = \frac{Δr}{2 Δz} (\frac{\partial}{\partial \dot{z}} - \frac{\partial}{\partial \dot{x}}) 。 \end{matrix} (8) \end{matrix}$

式中 $\begin{matrix} \dot{x} \end{matrix}$ 和 $\begin{matrix} \dot{z} \end{matrix}$ 分别代表两个对角线的方向导数, Δ x和Δ z分别是沿坐标轴x、z方向的空间步长, Δ r是对角线的长度, 且Δ r= $\begin{matrix} \sqrt[]{Δ x^{2} + Δ z^{2}} \end{matrix}$ 。 $\begin{matrix} \dot{x} \end{matrix}$ 和 $\begin{matrix} \dot{z} \end{matrix}$ 离散差分算子的计算与普通交错网格有限差分中差分算子的计算相同。

4 模型试算

4.1 均匀介质波场模拟及分解

为了验证本文分解算法的正确性, 对均匀TTI介质进行了波场模拟及波场分解处理, 其物性参数如表1所示, 极化角和方位角均为45° , 震源采用z方向集中力震源, 雷克子波主频为30 Hz, 网格点数为500× 500, 网格间距为10 m, 模拟时间步长为1 ms。鉴于篇幅本文仅给出x分量波场分解前后对比图(图2)和单道能量及相位对比图(图3), 图3中x分量为分解波场分解前的能量图, x_P、x_SV、x_SH分别为波场分解之后对应的P、SV、SH波的能量图。

表1 TTI介质模型各向异性参数

	Figure Option View Download New Window
	图2 x分量分解前后波场快照对比

	Figure Option View Download New Window
	图3 x分量分解前和分解后单道能量及相位对比

图3为抽取图2中x分量垂向网格200, 水平网格(0~500)处的全波、P波、S波波场的能量对比, 可以看出x_P和x分量中的P波部分振幅相同且相位不变, x_SV和x分量中的SV波部分振幅相同且相位不变, x_SH和x分量中的SH波部分振幅相同且相位不变。分解前后的能量和相位完全相同, 因此将分解之后的x分量的P波与S波波场相加得到的叠加全波波场和模拟所得全波波场x分量相同。同理, y、z分量也具有这样的特征。

4.2 复杂介质模型

设计了一个复杂模型, 网格点数为700× 400, 网格间距为5 m, 模拟时间步长为0.5 ms, 采样间隔为1 ms, 模拟时长为2.5 s, 震源采用30 Hz雷克子波垂向集中力源, 震源在模型(200 m, 1 750 m)处激发, 检波器布设在地表。对极化方位各向异性TTI介质(极化角45° , 方位角60° )进行波场模拟及分解处理。模型的物性参数见表2, 纵波速度分布如图4。

表2 复杂层状VTI介质各向异性参数

	Figure Option View Download New Window
	图4 复杂层状介质纵波速度

图5为x、y、z各分量分解前后波场快照, 图6为波场分解得到的地表地震记录。图5从左至右依次为全波、P波、SV波、SH波波场快照, 从图中可知, 与均匀TTI介质下所得的结论相同, x、y和z分量中的P波、SV波、SH波被完全分开, 且分解之后的波场能量和相位没有变化, 将波场分解之后的地震波场用于逆时偏移将得到更精确的地下界面成像。

	Figure Option View Download New Window
	图5 x、y、z分量分解前后波场快照

	Figure Option View Download New Window
	图6 x、y、z分量分解前后地震记录

为了能够对比分解前后地震道能量, 文中将分解前后的地震记录放在同一图中绘制(图6), 从左至右依次为全波、P波、SV波、SH波。从x、y、z三分量地震记录对比可以看出, 分解前后地震记录中同相轴的能量及位置没有变化, 验证了本文采用的波场分解方法是保幅的、无相变的。

5 结论

通过均匀和复杂TTI介质模型波场分解前后的快照和地震记录的对比, 验证了将TTI介质中传播的地震波映射到P、SV、SH极化向量方向进行波场分解方法的正确性, 分解得到的地震波的振幅、相位及物理意义和原始输入的地震波相同。

采用旋转交错网格进行TTI介质的波场模拟和波场分解, 相较于普通交错网格具有以下几点优势:①不需要对各向异性参数进行平滑处理, 提高了TTI介质波场模拟的精度; ②速度分量在相同位置, 不需要做平均处理提高了波场分解的精度; ③不需要对速度分量在空间域作平均处理或在波数域做插值处理提高了运算效率。

The authors have declared that no competing interests exist.

参考文献

View Option

[1]	Amundsen L, Reitan A. Decomposition of multicomponent sea-floor data into upgoing and downgoing P- and S-waves[J]. Geophysics, 1995, 60(2): 563-572. [本文引用:1]
[2]	Sun R, McMechan G A, Chuang H H. Amplitude balancing in separating P- and S-waves in 2D and 3D elastic seismic data[J]. Geophysics, 2011, 76(3): S103-S113. [本文引用:1]
[3]	van Der Baan M. PP/PS Wavefield separation by independent component analysis[J]. Geophysical Journal International, 2006, 166(1): 339-348. [本文引用:1]
[4]	Cheng J, Kang W. Propagating Pure Wave Modes in 3D General Anisotropic Media: PART I—P-Wave Propagator[C]//Expand ed Abstracts of 82^nd SEG Annual International Meeting, 2012. [本文引用:1]
[5]	Kang W, Cheng J. Propagating Pure Wave Modes in 3D General Anisotropic Media: PART II—SV and SH Wave[C]//Expand ed Abstracts of 82^nd SEG Annual International Meeting, 2012. [本文引用:1]
[6]	马德堂, 朱光明. 弹性波波场P波和S波分解的数值模拟[J]. 石油地球物理勘探, 2003, 38(5): 482-486. [本文引用:1]
[7]	李振春, 张华, 刘庆敏, 等. 弹性波交错网格高阶有限差分法波场分离数值模拟[J]. 石油地球物理勘探, 2007, 42(5): 510-515. [本文引用:1]
[8]	Dellinger J A. Anisotropic seismic wave propagation[D]. Stanford University, 1991. [本文引用:1]
[9]	Dellinger J, Etgen J. Wave-type separation in 3-D anisotropic media[C]//Expand ed Abstracts of 59^th SEG Annual International Meeting, 1989: 977-979. [本文引用:1]
[10]	Yan J, Sava P. Elastic wave-mode separation for VTI media[J]. Geophysics, 2009, 74(5): WB19-WB32. [本文引用:1]
[11]	Zhang Q, McMechan G A. 2D and 3D elastic wavefield vector decomposition in the wavenumber domain for VTI media[J]. Geophysics, 2010, 75(3): D13-D26. [本文引用:1]
[12]	Du Q Z, Zhang M Q. Wave-number domain interpolatin of staggered grid for true-amplitude wavefield separation[C]//Expand ed Abstracts of 83^rd SEG Annual International Meeting, 2013. [本文引用:1]
[13]	Thomsen L. Elastic anisotropy due to aligned cracks in porous rock[J]. Geophysical Prospecting, 1995, 43(6): 805-829 [本文引用:1]
[14]	Bohlen T, Saenger E H. Accuracy of heterogeneous staggered-grid finite-difference modeling of Rayleigh waves[J]. Geophysics, 2006, 71(4): T109-T115. [本文引用:1]

1995

0.0

... 完全弹性波动方程数值模拟中,纵、横波耦合的问题是研究纯纵波和纯横波的传播规律以及利用纯纵、横波逆时偏移成像的难题,因此波场分解和纯纵、横波波场模拟越来越受到重视^[1,2,3,4,5] ...

2011

0.0

2006

0.0

2012

0.0

2012

0.0

2003

0.0

马德堂, 朱光明. 弹性波波场P波和S波分解的数值模拟[J]. 石油地球物理勘探, 2003, 38(5): 482-486.

Common-used the foil elastic wave equation-to carry out the numerical modeling of elastic wavefield only obtained a hybrid wavefield of P-wave and S-wave.When you want to obtain P-wave and S-wave wavefield,generally the wavefield separation is carried out after hybrid wavefield being obtained or using acoustic equation to model the P-wave and S-wave wavefield,but these methods are very difficult to ensure the amplitude of P-wave and S-wave do not disorder.The paper gave an equivalent equation of full elastic wave equation, which contains variations of both the hybrid wavefield and the pure P-wave and S-wave wevefield. Using pseudo-spectrum method to solve the wave equation,obtained totally separated wavefield of pure P-wave and S-wave while obtaining hybrid wavefield,and preserved the information of mutual transformation of P-wave and S-wave energy.The analysis of modeling wavefield showed using nu-merical modeling of elastic wave separation has lmportant meaning in respect of knowing the propagating rule of elastic wave.

Department of Basic Course Chang'an University, Xi'an City, Shanxi Province 710064, China

通常采用完全弹性波波动方程进行弹性波波场数值模拟,只能得到P波和S波的混合波场。若要进一步获得P波和S波的波场,一般是在得到混合波场之后再进行波场分离,或用声波方程单独模拟P波和S波波场,但是这些方法往往很难确保P波和S波的振幅不出现畸变。本文给出了完全弹性波波动方程的一种等价方程,该方程既含混合波场变量,也含纯P波和纯S波波场变量。用伪谱法求解该波动方程,在得到混合波场的同时,也得到了完全分离的纯P波或纯S波波场,并保留了P波和S波能量相互转换的信息。模拟波场分析表明,采用这种弹性波波场分解的数值模拟方法对认识弹性波的传播规律有重要意义。

... 马德堂等^[6]、李振春等^[7]在散度旋度的理论上提出等价方程,模拟了纯P波和纯S波的传播特征 ...

2007

0.0

李振春, 张华, 刘庆敏, 等. 弹性波交错网格高阶有限差分法波场分离数值模拟[J]. 石油地球物理勘探, 2007, 42(5): 510-515.

地震波中的弹性波传播会产生纵波和横波，所以用完全弹性波波动方程进行弹性波波场数值模拟时，只能得到纵横波耦合的混合波场。本文从P波波场为无旋场，而S波波场为无散场的思路出发，推出满足此条件的一阶速度—应力弹性波波场分离方程，并利用交错网格高阶有限差分法对波场分离方程进行数值模拟。模拟实例表明，此方法不但成功地将P波波场和S波波场从混合波场中分离出来，而且这种方法的稳定性好、模拟精度高，可用于弹性波传播规律研究及地震资料处理。

... 马德堂等^[6]、李振春等^[7]在散度旋度的理论上提出等价方程,模拟了纯P波和纯S波的传播特征 ...

1991

0.0

... Dellinger^[8,9]提出通过将地震波场投影到纵、横波的极化向量上来分离各向异性介质中的纵波,然后用全波场减去纵波波场来得到非纵波场 ...

1989

0.0

... Dellinger^[8,9]提出通过将地震波场投影到纵、横波的极化向量上来分离各向异性介质中的纵波,然后用全波场减去纵波波场来得到非纵波场 ...

2009

0.0

... Yan 和 Sava^[10]将波数域中的波场分离算子进行反傅立叶变换,得到空间域随空间变化的分离算子,该算子能适应物性参数缓慢变化的复杂介质,计算量大,且精度会受到物性参数剧烈变化的影响 ...

2010

0.0

... Qunshan Zhang^[11]引入了波场分解算子,将VTI介质中的全波波场分解为P波、SV波波场,该算子克服了散度、旋度波场分解算子引起的相移,但其在得到最终的分解结果时,采用了平均策略 ...

2013

0.0

... 杜启振^[12]指出平均策略降低了波场分解的精度,所以采用在波数域对波场进行插值处理来提高精度 ...

1995

0.0

... 在旋转交错网格中,所有的速度分量被定义在一个位置,所有的应力分量被定义在和速度分量相对的角点上^[13] ...

2006

0.0