地球物理学报  2019, Vol. 62 Issue (1): 260-275   PDF    
正交各向异性介质中qP波入射的二阶近似反射系数与透射系数
许茜茹1,2, 毛伟建1     
1. 中国科学院测量与地球物理研究所计算与勘探地球物理中心, 大地测量与地球动力学国家实验室, 武汉 430077;
2. 中国科学院大学, 北京 100049
摘要:地震波在各向异性介质中以一个准P波(qP)和两个准S波(qS1和qS2)的形式传播.研究三种波的相速度、群速度以及偏振方向等传播性质能够为各向异性介质中的正反演问题提供有效支撑.具有比横向各向同性(TI)介质更一般对称性的正交各向异性介质通常需要9个独立参数对其进行描述,这使得对传播特征的计算更为复杂.当两个准S波速度相近时具有耦合性,从而令慢度的计算产生奇异性.因此,奇异点(慢度面的鞍点和交叉点)附近的反射与透射(R/T)系数的求解不稳定,会导致波场振幅不准确.本文首次通过结合耦合S波射线理论和基于迭代的各向异性相速度与偏振矢量的高阶近似解,得到了适用于正交各向异性介质以qP波入射所产生的二阶R/T系数的计算方法.与基于一阶近似的结果相比,基于二阶近似的方法提高了qP波R/T系数的精度,能得到一阶耦合近似无法表达的准确的qP-qS转换波的R/T系数解,且方法适用于较强的各向异性介质.
关键词: 各向异性      R/T系数      耦合S波      二阶近似     
Reflection and transmission coefficients based on second-order approximation in orthorhombic anisotropic media
XU QianRu1,2, MAO WeiJian1     
1. Chinese Academy of Sciences, Institute of Geodesy and Geophysics, Center for Computational and Exploration Geophysics and State Key Laboratory of Geodesy and Earth's Dynamics, Wuhan 430077, China;
2. University of Chinese Academy of Sciences, Beijing 100049, China
Abstract: Seismic wave propagation in the anisotropic media with one quasi-P wave (qP) and two quasi-S waves (qS1 and qS2). Studying of the propagation properties, such as phase velocity, group velocity and polarization, provides effective support for seismic forward modeling and inverse problems in anisotropic media. To describe the orthorhombic anisotropic media with more general symmetry than transverse isotropic (TI) media usually requires 9 independent parameters, which makes the calculation of propagation properties very complicated. The two quasi-S waves with similar velocities are coupled, causing the singularities in slowness calculation. As a result, the reflection and transmission (R/T) coefficients are unstable in the neighborhood of singularities (saddle points and cross points on slowness surface), which makes the wavefield amplitude calculation inaccurate. In this paper we propose a method to obtain the second-order R/T coefficients generated by incident qP wave in orthorhombic media by using coupled ray theory and high-order iterative solution of phase velocity and polarization vector in anisotropic media. Compared with the first-order results, the proposed method improves the accuracy of qP wave R/T coefficients, and obtains accurate qP-qS converted wave R/T coefficients which is not available in the first-order coupling approximation method. Furthermore, our method is applicable to strong anisotropic media.
Keywords: Anisotropy    R/T coefficients    Couple-S wave    Second-order approximation    
0 引言

射线理论作为波动方程的高阶近似解,被广泛应用于地震学研究中.射线在给定介质中的传播特征作为射线类正反演理论的基础,已有大量的研究成果.

针对地球介质,各向异性假设更符合真实的地层介质特征,不同的各向异性介质需要使用相应的参数化方法进行描述.Thomsen(1986)假设具有单个对称轴的TI介质(VTI, HTI, TTI)能够描述地层的结构属性,给出了弱各向异性VTI介质的参数,并据此推导了相速度和群速度的近似表达式.而当介质中存在与地层走向垂直的裂隙时,使用正交各向异性模型能够对其更合理地描述(卢明辉等,2005).Tsvankin(1997)基于Thomsen的TI介质参数,通过分别在三个对称面内求解Christoffel方程,推导了正交各向异性介质的参数化方法,并给出了弱正交各向异性介质中qP波相速度近似表达式.Xu和Stovas(2017)结合椭圆背景场和非椭圆系数提出了面向声波正交各向异性介质的参数化方法,并且基于Shanks变换得到了较精确的走时近似.

在各向异性介质中,地震波以准P波和两个准S波的形式进行传播,传播特征分别与介质特征决定的Christoffel方程的三组特征值和特征向量相对应,但仅在部分具有简单对称性的各向异性介质中存在解析解,大部分情况下需要使用数值方法进行求解.吴国忱等(2010)推导了TTI介质的相速度与偏振方向的解析表达式;利用TI介质的传播特征信息,可以进行相应的正演模拟和逆时偏移计算(黄金强和李振春, 2017; 李振春等, 2017);Farra(2001)使用迭代方法推导了各向异性介质中高阶qP波和qS波相速度以及偏振方向的表达式,在未增加过多计算代价的情况下,可以得到相对精确的高阶解.Sripanich和Fomel(2015)使用了椭圆近似得到TI介质和正交各向异性介质中qP波的相速度和群速度解,并使用其进行波向量分解(Sripanich等,2017),Hao和Stovas(2016)利用广义时差近似(GMA)也给出了正交各向异性介质的qP波高精度近似相速度,达到了很高的精度,但是后两种方法仅涉及qP波的近似表达.

在各向异性介质中,当两个准S波相近时存在耦合现象,即在弱各向异性介质中或者一般各向异性介质的奇异点(慢度面的鞍点和交叉点)附近,常规的射线追踪方法通常无法得到稳定的解(Vavryčuk, 2003; Grechka, 2015).Coates和Chapman(1990)探讨了qS波在弱各向异性介质中的耦合现象,并使用耦合体积分对射线走时和偏振方向进行改正.Klimeš(2006)Farra和Pšenčík(2008)提出了耦合S波射线理论,并且推导了光滑各向异性弹性介质中的耦合S波的一阶射线追踪表达式.耦合射线理论不再区分两个速度相近的准S波,而是将其看作按两个偏振方向进行传播的一组波,某种程度上解决了准S波场在奇异点附近的解不稳定的问题.

R/T系数反映了波场振幅在界面处的变化情况.计算与位移相关的R/T系数需要分别对广义Snell定理和Zoeppritz方程进行求解,即令入射波和反射、透射波在界面切平面内的慢度分量相等,并且要求位移、应力连续.Slawinski等(2000)基于Thomsen的VTI介质中的相速度以及相速度面和偏振方向的几何关系,推导了适用于VTI介质的广义Snell定理.Daley和Hron (1977)Chapman(1994)对半空间条件下的两个均匀各向异性介质相关的R/T系数进行了研究. Červený(2001)使用射线理论计算了3维层状各向异性介质间的R/T系数.Jílek (2002)使用扰动方法推导了界面上下介质弱差异的弱各向异性相关的PS波反射系数.Farra和Pšenčík(2010)基于一阶射线追踪和耦合S波射线理论给出了弱各向异性介质中qP-qP非转换波的一阶近似R/T系数,但并没有给出qP-qS转换波的R/T系数解.梁锴等(2011)利用弱各向异性近似推导了TTI介质qP波入射精确和近似R/T系数.

本文通过Tsvankin的参数化方法和Bond矩阵定义正交各向异性介质、TI介质和各向同性介质,结合耦合S波射线理论与各向异性相速度和偏振矢量的高阶近似解得到了基于二阶近似的Snell定理,使用牛顿法对二阶qP波和耦合S波的反射、透射波慢度进行迭代求解,根据Zoeppritz方程得到基于二阶近似的以qP(P)波入射的转换波和非转换波的R/T系数并且将二阶近似解与对应的精确解、一阶近似解进行比较,二阶近似解的非转换波R/T系数精度高于一阶近似,并且能够得到两qS波相关的R/T系数解.

1 理论方法 1.1 各向异性介质的刚度矩阵参数化

为了讨论基于正交各向异性的R/T系数,本文使用Tsvankin(1997)提出的9参数对正交各向异性介质进行参数化,参数及其物理意义见附录A.此参数化方法以正交各向异性介质的对称轴为坐标系(即各向异性参考坐标系),能够退化表达TI介质和各向同性介质.

当各向异性参考坐标系与全局坐标系(z轴沿竖直方向的正交坐标系)不一致时,为了保持文中坐标系的统一,需要进行矩阵旋转操作.假设沿各向异性参考坐标系三轴正方向的单位向量在全局坐标系中表示为(e1, e2, e3),全局坐标系的三轴单位向量为(e1, e2, e3),可定义作用在刚度矩阵上的旋转矩阵为

(1)

与旋转矩阵R对应的Bond矩阵M

其中,rij为旋转矩阵R的元素.全局坐标系下的刚度矩阵Aglobal可以通过Bond矩阵对在各向异性参考坐标系下的刚度矩阵Alocal进行变换,表示为

(2)

其中,Aglobal, Alocal均为Voigt格式的6×6刚度矩阵,不做特殊说明的情况下,下文中出现的刚度矩阵均是基于全局坐标系的.

1.2 各向异性介质的二阶近似相速度

与慢度方向矢量n相关的Christoffel矩阵Γ(x, n)为

(3)

其中,弹性常量的四阶张量Cijkl经密度ρ正则化为aijkl.由于其具有对称性,与式(2)中Voigt格式表示的Aglobal, Alocal等价,具有映射关系(Thomsen, 1986).在各向异性介质中,式(3)的三个不同的特征值和特征向量分别对应qP、qS1、qS2波的速度平方和偏振方向.

S波耦合现象经常出现于弱各向异性介质中或者强各向异性介质的两个qS波波速相似的奇异点区域.针对这种情况,我们引入耦合射线的概念(Klimeš, 2006).耦合射线理论假设两个速度相近的qS波按照“同一”路径传播,分别计算其偏振属性和动力学参数,能够有效处理两个qS波,不会产生奇异解,避免了常规各向异性射线理论无法处理qS波耦合的问题,且形式简洁,可进一步简化用以处理各向同性介质的S波.对于一阶近似,各向异性介质中的qP波和耦合S波的特征值可以分别表示为(Farra和Pšenčík, 2008)

(4)

下文中,令qP(P)波相关的变量下标为P,耦合S波相关的变量下标为SM.式(4)中Bjl为对称矩阵B(x, n)的元素,

(5)

其中,ei[j]为一组三维正交基向量e[1], e[2], e[3]中的第j个向量的i分量,其中,e[3]=ne[1], e[2]为一组与e[3]正交的向量.需要说明的是,对任意e[1], e[2]来说,B11+B22B132+B232为定值,因此,其选择不影响一阶及下述的二阶特征值的大小.由式(4),qP波和耦合S波对应的相速度分别为

对于各向异性介质中qP波和耦合S波的二阶近似,我们先利用Farra(2001)给出的基于扰动迭代的高阶qP波和qS1、qS2波形式推导出最终形式为

(6)

其中,

(7)

令两个qS波沿耦合射线路径传播,即GS1(1)=GS2(1)=GSM(1),结合式(4),式(6),可以得到qP波和耦合S波的二阶近似

(8)

将耦合射线理论与二阶迭代近似相速度计算结合,得到二阶近似的qP波和耦合S波波速值.与Farra(2001)中二阶扰动解不同,此二阶近似式(8)的分母中,我们使用了一阶近似的qP波和耦合S波特征值GP(1)(x, n), GSM(1)(x, n),而非原公式中迭代二阶近似中使用的参考各向同性介质的P波和S波特征值,将迭代二阶近似向迭代三阶近似逼近,增加了二阶近似相速度的精度,形式简单,并且没有显著提升所需的计算代价.二阶近似qP波和耦合S波的相速度大小分别为慢度矢量为p(2)=n/v(2).

1.3 广义Snell定理和基于二阶近似的Snell定理

与式(3)类似,坐标x处的与慢度相关的Christoffel矩阵Γ

(9)

其中,p为特定波的慢度矢量,p=n/vv表示其相速度大小.各向异性介质中某种波的存在条件为

(10)

x位于界面上,则入射波与反射、透射波的慢度遵循Snell定理,即在界面切平面的投影分量a相等,写作向量形式为

(11)

其中,pinc为入射波的慢度矢量,Nx处的界面法向量,pR/T为反射、透射波的慢度矢量.由式(11)可知,

(12)

其中,σpR/T在界面法线方向上的投影.结合式(12)与(10),可以得到关于σ的基于广义Snell定理的六阶方程,忽略物理上不合理的三个解,从而可求得对应于qP、qS1、qS2波的三个解(Červený, 2001,公式2.3.53,2.3.54),但是由于复杂各向异性存在奇异点,当两个准S波对应的特征值相似时,会产生奇异解.Farra和Pšenčík(2010)提出的基于一阶近似的R/T系数计算虽然能够处理弱各向异性的qP波折射和反射,但是其无法适应强各向异性介质,并且其并未给出转换耦合S波相关的R/T系数.因此,我们结合二阶近似的相速度式(8)推导基于二阶近似的Snell定理.与式(10)对应的二阶qP波和耦合S波特征值可以表示为

(13)

其中,pP, pSM分别为两者的慢度矢量,同样遵循(12)式.结合式(12)与(13)可以得到关于σ的最高为12阶的高次方程.方程的阶数由各向异性复杂程度、全局坐标系与各向异性参考坐标系的差异决定.例如,各向同性介质对应的两个方程均为8阶,解包含两组三重根在内,共有4个,与一阶近似相同;VTI介质对应的两个方程也为8阶,其8个解可能含有重根;对称轴为全局坐标系三轴的正交各向异性介质也对应了8阶方程.而界面法线方向不沿竖直向上,且对称轴与全局坐标系三轴均不一致的正交各向异性介质,或更为复杂的各向异性介质则对应了最高阶为12阶的方程.虽然12阶的情况可以通过旋转全局坐标系为各向异性参考坐标系从而将方程降阶为一般的8阶,但是由于高阶方程系数的复杂度随着各向异性及地层的复杂程度增加,方程求解和筛选的计算代价也会超出可接受的范畴.因此,令以R/T系数相关的变量σ为未知数的目标函数为

(14)

使用牛顿法来迭代求解,第j次迭代的迭代公式为

(15)

其中,与∂G(2)/∂pk相关变量的完整形式见附录B.若解为复数,此方法仍能够有效求解出相应的σ.使用两组迭代求解可以分别得到二阶qP波和耦合S波的相速度慢度矢量.迭代求解避免了高阶方程求解后解的筛选过程,在有限的迭代次数中就能够得到合理的结果,比一般方程求解效率更高,且耦合S波保证了奇异点附近解的稳定性.

1.4 基于二阶近似的qP入射波的Zoeppritz方程

偏振方向的求取基于Farra(2001),根据qP波和耦合S波的二阶慢度矢量确定的单位方向向量,可以得到其对应的二阶偏振方向向量g(2),具体表达形式为

(16)

其中,Bij参数分别为对应的qP波慢度或耦合S波慢度.GSk(2)(p)为替换式(6)中单位方向向量参数为对应耦合S波慢度矢量得到,其对应的ei[j]与式(5)中相同即可.当慢度为复数矢量时,对应的偏振矢量亦为复数矢量形式.

我们通过如下正交各向异性模型给出关于慢度和偏振方向计算的算例.

(17)

首先,使用式(17)参数化的正交各向异性模型(Schoenberg和Helbig, 1997),按照式(3)进行特征值求解,计算此模型在三个对称面对应的慢度分布情况,如图 1,可见在正交各向异性中的两个qS波对应部分入射角的速度相似,即在奇异点(慢度曲线交叉点)附近存在耦合现象,此时若按照传统的各向异性射线追踪方法进行求解则会出现不稳定的问题.下面我们将通过结合各向异性高阶近似解与耦合射线理论对正交各向异性介质进行处理.不做特殊说明的情况下,将对包括奇异点附近范围在内的所有入射角度进行计算,检验此近似对一般情况的整体适应性.

图 1 使用式(17)参数化的正交各向异性模型在三个对称面上的慢度分布 (a)入射向量(n1, n2, n3=0)确定的水平平面内的慢度分布; (b)入射向量(n1=0, n2, n3)确定的竖直平面内的慢度分布; (c)入射向量(n1, n2=0, n3)确定的竖直平面内的慢度分布. Fig. 1 Slowness on three symmetric surfaces for orthorhombic model parameterized by (17) (a) Slowness on horizontal surface with incident vector (n1, n2, n3=0); (b) Slowness on vertical surface with incident vector (n1=0, n2, n3); (c) Slowness on vertical surface with incident vector (n1, n2=0, n3).

由于耦合S波理论的S波与常规的qS1、qS2波波速不直接可比,我们将(8)式中的二阶近似qP波、Farra的二阶、三阶近似解与精确解进行比较,得到如图 2所示的误差剖面,qP波相速度的精度随着阶数的增加而提高.其中,对于高阶近似qP波波速,图 2b使用了零阶近似P、S波波速,图 2c使用了qP波及耦合S波的一阶近似波速,图 2d使用了一阶近似qP、qS1、qS2波速.图 2c的精度介于图 2b图 2d之间,与一阶近似的图 2a相比,提高了一个数量级的精度,此二阶近似对R/T系数的改正作用在下文的实例中也能够体现出来.

图 2 近似qP波波速与精确解的误差(百分比) (a)一阶近似的误差; (b)二阶近似(Farra)的误差; (c)二阶近似(本文)的误差; (d)三阶近似(Farra)的误差. Fig. 2 Error (in percentage) of approximated qP wave phase velocity compared with exact one (a) Error of 1st order approximation; (b) Error of 2nd order approximation (Farra′s); (c) Error of 2nd order approximation (ours); (d) Error of 3rd order approximation (Farra′s).

接下来考察使用二阶近似计算的三个波对应的偏振矢量,按照式(16)计算偏振方向并绘制在相应入射方向nP, nSM的单位球上,如图 3,正交各向异性的qP波(对应nP)偏振方向与入射方向有一定偏差,但是仍较接近于入射方向;qS1波偏振方向和qS2波偏振方向(对应nSM)与入射方向的夹角不为90°,与Sripanich等(2017)Fig. 1的解趋势一致.

图 3 单位入射方向球上的正交各向异性模型的(a) qP, (b) qS1, (c) qS2偏振矢量 Fig. 3 Polarization vectors of (a) qP, (b) qS1, and (c) qS2 waves in the orthorhombic model plotted on the unit incident direction sphere

在界面法线单位矢量为N的坐标x处,二阶慢度矢量为p(2)的波对应的应力可以表示为

(18)

若在x点以qP(P)波入射,入射波的偏振矢量和应力矢量为gPinc, XPinc,由位移和应力连续,可得到反射\透射qP、耦合S波的偏振矢量和应力矢量gP, gS1, gS2XP, XS1, XS2,即Zoeppritz方程如下,

(19)

其中,

上标中的rt分别对应上层(入射)介质中的反射波以及下层(透射)介质中的透射波.CPr, CS1r, CS2r对应反射qP波、耦合S波的反射系数,CPt, CS1t, CS2t对应透射qP波、耦合S波的透射系数.此处及下文中的反射、透射系数均为以位移表示的反射、透射系数.方程(19)是一个较稳定的线性系统,可通过常规方法求解.

2 数值实验

为了验证本文中二阶近似R/T系数计算方法的有效性,我们使用三组各向异性模型来进行验证.在所有模型中,当上半空间介质为各向同性时,入射波为P波,反射波为P波,S波;当其介质为各向异性时,入射波为qP波,反射波为qP波和耦合S波.透射波在下半空间的各向异性介质中传播,透射波为qP波和耦合S波.为了表达和比较的统一和方便,后文也使用CS1r, CS2r来分别表示各向同性介质中的反射波的S1和S2分量.令入射倾角和方位角均从0°变化至90°,给出R/T系数及相关参数的等值线图.R/T系数的精确解是通过直接求解广义Snell定理((10),(12)联立的六阶方程)和对应的Zoeppritz方程得到,由于奇异点附近直接求解不稳定,需要额外计算和选择合理的qS波慢度、偏振及应力矢量.另外,由于使用式(19)求解的R/T系数可能为复数,我们在例子中均使用系数的模对各变量进行表达和比较.

2.1 P波从各向同性介质入射至弱各向异性HTI介质

首先,将二阶近似的R/T系数结果与已有的基于一阶近似的R/T系数结果进行对比,用于检测二阶近似的正确性和对弱各向异性的改正作用.下面的模型A为Farra和Pšenčík(2010)文中使用的第二组模型,入射P波所在的上层为各向同性介质,P波波速vP、S波波速vS以及密度ρ分别为

(20)

透射qP波和耦合S波所在的下层为弱HTI介质,可以将刚度系数矩阵换算为Tsvankin参数和旋转矩阵,即

(21)

由于其上层的各向同性介质的速度在一定倾角、方位角范围内小于下层介质,因此可能存在部分复R/T系数.图 4为模型A的qP(P)波R/T系数的精确解、二阶近似解与精确解的差值,以及一阶近似解与精确解的差值.因入射倾角和方位角分别在0到90°范围内变化,三维空间内的临界角可组成一条分界线,即图 4a所示的AB或图 4d所示的CD连线.在临界角分界线附近区域,|CPr|和|CPt|由实数变为复数,产生剧烈变化.比较图 4b4c4e4f在此区域的值,二阶近似解与精确解的差值要小于一阶近似解与精确解的差值.在整体上,对于弱各向异性介质,二阶近似解与一阶近似解的qP(P)波R/T系数的精度相近,但是在误差分布上,令误差的期望为0,二阶近似解、一阶近似解与精确解的差值的标准差分别为0.0086、0.0426(|CPr|)和0.0119、0.0506(|CPt|),即二阶近似解略优于一阶近似解,两者都具有较高的精度.

图 4 模型A的P波入射对应的P波反射系数和qP波透射系数 (a) P波反射系数的精确解; (b) P波反射系数的二阶近似解与精确解的差值; (c) P波反射系数的一阶近似解与精确解的差值; (d)(e)(f)与(a)(b)(c)类似,是关于qP波透射系数的值. Fig. 4 P wave reflection coefficient and qP wave transmission coefficient corresponding to incident P wave for model A (a) Exact P wave reflection coefficient; (b) Difference between second-order approximate and exact P wave reflection coefficient; (c) Difference between first-order approximate and exact P wave reflection coefficient; (d) (e) (f) are similar to (a) (b) (c) for qP wave transmission coefficients.

由于Farra计算的耦合S波R/T系数与精确解不可比(Farra和Pšenčík, 2010),此处及下文中,对于准S波(或S波),我们仅对二阶近似结果与精确解的差进行分析和比较.图 5为耦合S波的两个偏振方向(各向同性介质为S1、S2分量)对应的R/T系数的精确解以及二阶近似解与精确解的差值.由图 5可以看出,二阶近似解与精确解很接近,其标准差均小于0.0034(见表 1).

图 5 模型A的以P波入射对应的S波反射系数、qS波透射系数的精确解和二阶近似解与精确解的差值 (a) S波的S1分量反射系数精确解; (b) qS1波的透射系数精确解; (c) S波的S2分量反射系数精确解; (d) qS2波透射系数精确解; (e)(f)(g)(h)为二阶近似解与相应精确解(a)(b)(c)(d)的差值. Fig. 5 S wave reflection coefficients and qS wave transmission coefficients corresponding to incident P wave for model A (a) Exact S1 component of S wave reflection coefficient; (b) Exact qS1 wave transmission coefficient; (c) Exact S2 component of S wave reflection coefficient; (d) Exact qS2 wave transmission coefficient; (e) (f) (g) (h) are difference between second-order approximate and exact coefficients corresponding to (a)(b)(c)(d).
表 1 模型A中qS波(S波)二阶近似R/T系数的标准差 Table 1 Standard deviation of second-order approximation of qS wave (S wave) R/T coefficients for model A

为了更直观地分析qS波(S波)二阶近似R/T系数的计算精度,我们分别给出了固定入射倾角和固定方位角为45°时,qS(S)波R/T系数值的截面(见图 6a6b).由图 6可清楚地看到二阶近似S波的S1、S2分量反射系数、耦合S波两偏振相关的透射系数与其精确值相当接近.

图 6 当(a)固定入射倾角为45°,(b)固定入射方位角为45°时,模型A的qS(S)波R/T系数精确解(实线)和二阶近似解(虚线) Fig. 6 The exact (solid line) and 2nd-order (dashed line) qS(S) wave R/T coefficients for model A when (a) the incidence is fixed as 45°, (b) the azimuth is fixed as 45°

基于对模型A的上述比较和分析,二阶近似的qP(P)波R/T系数,与Farra的一阶近似结果精度近似,整体上优于一阶近似解,这是由于对于各向同性介质和弱各向异性介质来说,一阶近似已能够达到所需的精度要求.而二阶近似的qS(S)波R/T系数精度也较高,能够较准确地处理弱各向异性介质.

2.2 P波从各向同性介质入射至强正交各向异性介质

为了验证基于二阶近似的R/T系数计算方法同样适用于更为一般的强正交各向异性介质,我们选用Tsvankin(1997)提供的强正交各向异性参数作为模型的一部分来说明二阶近似算法的适用性.

模型B中的入射P波所在的上层为各向同性介质,P波波速vP、S波波速vS以及密度ρ分别为

(22)

透射波所在的下层为强正交各向异性介质,其定义在全局坐标系下的各向异性参数和旋转矩阵为

(23)

图 7可知,对于强正交各向异性,一阶近似与二阶近似的R/T系数差别较大.一阶近似的透射系数|CPt|的最大相对误差为12%,反射系数|CPr|的最大相对误差则高达40%,一阶近似解在陡倾角区域精度无法达到要求.与之相应的,二阶近似的透射系数|CPt|和反射系数的最大相对误差分别为3%和6%,二阶近似的精确度整体上比一阶近似解要高约一位有效数字.在误差分布上,二阶近似和一阶近似与精确解的标准差分别为0.0030、0.0207(|CPr|)和0.0043、0.0310(|CPt|),即二阶近似的精度更高.

图 7 模型B的P波入射对应的P波反射系数和qP波透射系数,其它与图 4类似 Fig. 7 P wave reflection coefficient and qP wave transmission coefficient corresponding to incident P wave for model B. Others are similar to Fig. 4, but for model B

图 5类似,比较图 8中二阶近似的CS1r, CS1t, CS2r, CS2t与精确解的误差,其标准差如表 2.图 9a9b分别给出了固定入射倾角和方位角为45°时,qS(S)波的R/T系数大小,与上述分析一致,二阶近似与精确值匹配程度很高.因此,对于强各向异性介质,基于二阶近似的qS波(S波)R/T系数在图示入射角度分布范围内仍然同样具有较高的精度.

图 8 模型B的以P波入射对应的S波反射系数、qS波透射系数的精确解和二阶近似解与精确解的差值,其他与图 5类似 Fig. 8 S wave reflection coefficients and qS wave transmission coefficients corresponding to incident P wave for model B. Others are similar to Fig. 5, but for model B
表 2 模型B中qS波(S波)二阶近似R/T系数的标准差 Table 2 Standard deviation of second-order approximation of qS wave (S wave) R/T coefficients for model B
图 9 当(a)固定入射倾角为45°,(b)固定入射方位角为45°时,模型B的qS(S)波R/T系数精确解(实线)和二阶近似解(虚线) Fig. 9 The exact (solid line) and 2nd-order (dashed line) qS (S) wave R/T coefficients for model B when (a) the incidence is fixed as 45°, (b) the azimuth is fixed as 45°

对于强各向异性介质来说,一阶近似已经无法准确计算R/T系数,而考虑了速度和偏振方向二阶修正项的R/T系数能够较好地逼近准确值,因此,基于二阶近似R/T系数计算方法同样能够适应强各向异性介质.

2.3 qP波从弱各向异性TTI介质入射至倾斜的强正交各向异性介质

以上的例子均使用含有较简单的倾角(90°或0°)信息的各向异性旋转矩阵,下面给出一个例子来说明当各向异性坐标系的各对称轴均与界面法线方向不符时,基于二阶近似的R/T系数计算仍旧有效.

模型C的入射qP波所在的上层为TTI介质,参数分别为

(24)

同时,透射波所在的下层为强正交各向异性介质,参数为(Tsvankin, 1997)

(25)

定义上下两层的旋转矩阵分别为

(26)

通过对比图 10中qP波R/T系数的二阶近似、一阶近似与准确解的差,可知,与模型B类似,当下层介质为强各向异性时,一阶近似的qP波R/T系数精度与弱各向异性的情况相比,有降低趋势.由于模型C下层的各向异性比模型B稍弱,其一阶近似解相较模型B的精度略有提高,但反射系数和折射系数的最大相对误差仍然达到近25%和10%,而与之对应的二阶近似解的最大相对误差则约为10%和5%,比一阶提高了一倍左右, 仍具有相对较高的精度.二阶近似和一阶近似的标准差分别为0.0080、0.0141(|CPr|)和0.0085、0.0544(|CPt|)也佐证了上述结论,即二阶近似解更接近精确解.

图 10 模型C的qP波入射对应的qP波反射/透射系数,其他与图 4类似 Fig. 10 qP wave R/T coefficients corresponding to incident qP wave for model C. Others are similar to Fig. 4, just for model C

接下来考察qS波相关的二阶近似解,如图 11图 12CS1r, CS1t, CS2r, CS2t的标准差见表 3.四组R/T系数的精度对大多数入射角度来说较高,对波从弱各向异性介质入射至强各向异性介质的情况,二阶近似解也能较好地对R/T系数进行表达.

图 11 模型C的以qP波入射对应的qS波反射/透射系数的精确解和二阶近似解与精确解的差值,其他与图 5类似 Fig. 11 qS wave R/T coefficients corresponding to incident qP wave for model C. Others are similar to Fig. 5, just for model C
图 12 当(a)固定入射倾角为45°,(b)固定入射方位角为45°时,模型C的qS波R/T系数精确解(实线)和二阶近似解(虚线) Fig. 12 The exact (solid line) and 2nd-order (dashed line) qS waves R/T coefficients for model C when (a) the incidence is fixed as 45°, (b) the azimuth is fixed as 45°
表 3 模型C中qS波二阶近似R/T系数的标准差 Table 3 Standard deviation of second-order approximation of qS waves R/T coefficients for model C

当界面法线方向与强各向异性定义坐标系的基向量不符时,处理各向异性的Snell定理表达式更为复杂,但是此方法得到的二阶近似的R/T系数仍然能与准确解高度近似,因此,使用与旋转矩阵相关的Bond矩阵计算二阶近似的R/T系数是有效的.

3 结论

复杂各向异性介质的描述需要同时考虑地质构造和介质的各向异性特征,使用合理的各向异性参数化方法能够更简便地表达介质信息.在正交各向异性中,相速度、群速度及偏振矢量等传播特征的计算比TI介质更为复杂,在特定的入射角范围内具有奇异性,qS波在速度相近时的耦合性,且一阶近似仅针对弱各向异性介质.为了解决上述问题,我们通过结合基于迭代的相速度和偏振矢量的高阶近似与耦合S波射线理论,提出基于二阶近似的R/T系数计算方法,对所有小于临界角的入射角均获得了较精确的解.奇异点附近的二阶近似解的误差主要由近似慢度、偏振矢量与真实解的误差传递产生,由于其避免了奇异点异常,此误差在可接受的范畴内.根据此方法计算的二阶近似的qP波(P波) R/T系数与Farra和Pšenčík(2010)的一阶近似结果相比,在精度上有一定的提升,并且由Farra方法无法计算的qP-qS转换波的R/T系数也能够通过我们的二阶近似表达,在不增加过多计算量的前提下,具有较高的精度.同时,对含有强各向异性介质的模型也有一定的适应性,不需要界面上下介质弱差异的条件.此方法可以配合常规的R/T系数计算方法,针对强各向异性奇异点附近区域进行特殊处理.由于我们除参数化方法外,给出的表达式均为针对刚度矩阵的函数,因此,能够便于推广到比正交各向异性更复杂的各向异性介质正反演计算中.

附录 附录A Tsvankin各向异性参数表示及物理意义

设矩阵A为Voigt形式密度正则化的3×3的刚度矩阵,Tsvankin的9参数可以表示为

(A1)

其中,vP为P波在竖直方向的速度分量大小,vS为在x1方向偏振的S波竖直方向上的分量,ε1, γ1, δ1相当于[x1, x3]对称面上的VTI的ε, γ, δ参数,ε2, γ2, δ2相当于[x2, x3]对称面上的VTI的ε, γ, δ参数,δ3相当于[x1, x2]对称面上的VTI的δ参数.使用(A1)中所示的9个参数表达VTI介质时,令ε1=ε2, γ1=γ2, δ1=δ2δ3=0,另外,在表达各向同性的介质时,ε1=ε2=γ1=γ2=δ1=δ2=δ3=0.HTI,TTI以及对称轴倾斜的正交各向异性介质可以通过增加一个关于坐标系的旋转矩阵(1)进行定义.

附录B 二阶Snell定理相关变量表达式

为了辅助基于二阶近似的Snell定理的表示和推导,我们给出一组三维正交基向量:

(B1)

与式(5)类似,定义Bjl为对称矩阵B(x, p)的元素

(B2)

由此,根据式(8)和(B2),在x处与对应慢度矢量p相关的二阶近似的特征值及其关于p的一阶导数为

(B3)

其中,各符号分别为

(B4)

由于bij可直接通过(B2)和(B4)求得,下面只列出偏导数相关项:

(B5)

(B6)

(B7)

(B8)

(B9)

当介质为定义在全局坐标系下的正交各向异性时,(B6)—(B9)中除了A11, A22, A33, A44, A55, A66, A12, A13, A23项(即(A1)中定义项),其余项均为0.

References
Červený. 2001. Seismic Ray Theory. Cambridge: Cambridge University Press.
Chapman C H. 1994. Reflection/transmission coefficient reciprocities in anisotropic media. Geophysical Journal International, 116(2): 498-501. DOI:10.1111/j.1365-246X.1994.tb01811.x
Coates R T, Chapman C H. 1990. Quasi-shear wave coupling in weakly anisotropic 3-D media. Geophysical Journal International, 103(2): 301-320. DOI:10.1111/j.1365-246X.1990.tb01773.x
Daley P F, Hron F. 1977. Reflection and transmission coefficients for transversely isotropic media. Bulletin of the Seismological Society of America, 67(8): 661-675.
Farra V. 2001. High-order perturbations of the phase velocity and polarization of qP and qS waves in anisotropic media. Geophysical Journal International, 147(1): 90-104. DOI:10.1046/j.1365-246X.2001.00510.x
Farra V, Pšeník I. 2008. First-order ray computations of coupled S waves in inhomogeneous weakly anisotropic media. Geophysical Journal International, 173(3): 979-989. DOI:10.1111/j.1365-246X.2008.03778.x
Farra V, Pšeník I. 2010. First-order reflection/transmission coefficients for unconverted plane P waves in weakly anisotropic media. Geophysical Journal International, 183(3): 1443-1454. DOI:10.1111/j.1365-246X.2010.04794.x
Grechka V. 2015. Shear-wave group-velocity surfaces in low-symmetry anisotropic media. Geophysics, 80(1): C1-C7. DOI:10.1190/GEO2014-0156.1
Hao Q, Stovas A. 2016. Analytic calculation of phase and group velocities of P-waves in orthorhombic media. Geophysics, 81(3): C79-C97. DOI:10.1190/geo2015-0156.1
Huang J Q, Li Z C. 2017. Modeling and reverse time migration of pure quasi-P-waves in complex TI media with a low-rank decomposition. Chinese Journal of Geophysics (in Chinese), 60(2): 704-721. DOI:10.6038/cjg20170223
Jílek P. 2002. Converted PS-wave reflection coefficients in weakly anisotropic media. Pure and Applied Geophysics, 159(7-8): 1527-1562. DOI:10.1007/s00024-002-8696-9
Klimeš L. 2006. Common-ray tracing and dynamic ray tracing for S waves in a smooth elastic anisotropic medium. Studia Geophysica et Geodaetica, 50(3): 449-461. DOI:10.1007/s11200-006-0028-6
Li Z C, Huang J Q, Huang J P, et al. 2017. Fast least-squares reverse time migration based on plane-wave encoding for VTI media. Chinese Journal of Geophysics (in Chinese), 60(1): 240-257. DOI:10.6038/cjg20170120
Liang K, Yin X Y, Wu G C. 2011. Exact and approximate reflection and transmission coefficient for incident qP wave in TTI media. Chinese Journal of Geophysics (in Chinese), 54(1): 208-217. DOI:10.3969/j.issn.0001-5733.2011.01.022
Lu M H, Tang J H, Yang H Z, et al. 2005. P-wave traveltime analysis and Thomsen parameters inversion in orthorhombic media. Chinese Journal of Geophysics (in Chinese), 48(5): 1167-1171. DOI:10.3321/j.issn:0001-5733.2005.05.026
Schoenberg M, Helbig K. 1997. Orthorhombic media: Modeling elastic wave behavior in a vertically fractured earth. Geophysics, 62(6): 1954-1974. DOI:10.1190/1.1444297
Slawinski M A, Slawinski R A, Brown R J, et al. 2000. A generalized form of Snell′s law in anisotropic media. Geophysics, 65(2): 632-637. DOI:10.1190/1.1444759
Sripanich Y, Fomel S. 2015. On anelliptic approximations for qP velocities in transversely isotropic and orthorhombic media. Geophysics, 80(5): C89-C105. DOI:10.1190/geo2014-0534.1
Sripanich Y, Fomel S, Sun J Z, et al. 2017. Elastic wave-vector decomposition in heterogeneous anisotropic media. Geophysical Prospecting, 65(5): 1231-1245. DOI:10.1111/1365-2478.12482
Thomsen L. 1986. Weak elastic anisotropy. Geophysics, 51(10): 1954-1966. DOI:10.1190/1.1442051
Tsvankin I. 1997. Anisotropic parameters and P-wave velocity for orthorhombic media. Geophysics, 62(4): 1292-1309. DOI:10.1190/1.1444231
Vavryčuk V. 2003. Behavior of rays near singularities in anisotropic media. Physical Review B, 67(5): 054105. DOI:10.1103/PhysRevB.67.054105
Wu G C, Liang K, Yin X Y. 2010. The analysis of phase velocity and polarization feature for elastic wave in TTI media. Chinese Journal of Geophysics (in Chinese), 53(8): 1914-1923. DOI:10.3969/j.issn.0001-5733.2010.08.017
Xu S B, Stovas A. 2017. A new parameterization for acoustic orthorhombic media. Geophysics, 82(6): C229-C240. DOI:10.1190/GEO2017-0215.1
黄金强, 李振春. 2017. 基于Low-rank分解的复杂TI介质纯qP波正演模拟与逆时偏移. 地球物理学报, 60(2): 704-721. DOI:10.6038/cjg20170223
李振春, 黄金强, 黄建平, 等. 2017. 基于平面波加速的VTI介质最小二乘逆时偏移. 地球物理学报, 60(1): 240-257. DOI:10.6038/cjg2017012
梁锴, 印兴耀, 吴国忱. 2011. TTI介质qP波入射精确和近似反射透射系数. 地球物理学报, 54(1): 208-217. DOI:10.3969/j.issn.0001-5733.2011.01.022
卢明辉, 唐建侯, 杨慧珠, 等. 2005. 正交各向异性介质P波走时分析及Thomsen参数反演. 地球物理学报, 48(5): 1167-1171. DOI:10.3321/j.issn:0001-5733.2005.05.026
吴国忱, 梁锴, 印兴耀. 2010. TTI介质弹性波相速度与偏振特征分析. 地球物理学报, 53(8): 1914-1923. DOI:10.3969/j.issn.0001-5733.2010.08.017