2018, Vol. 61
2. 中国科学院大学, 北京 100049;
3. 中国科学院油气资源研究重点实验室, 北京 100029;
4. 中国科学院地球科学研究院, 北京 100029;
5. 中国科学院青藏高原地球科学卓越创新中心, 北京 100101
2. University of Chinese Academy of Sciences, Beijing 100049, China;
3. Key Laboratory of Petroleum Resources Research, Chinese Academy of Sciences, Beijing 100029, China;
4. Institutions of Earth Science, Chinese Academy of Sciences, Beijing 100029, China;
5. CAS Center for Excellence in Tibetan Plateau Earth Sciences, Beijing 100101, China
在石油勘探领域，横向各向同性介质(Transversely Isotropy，TI)是最常用来描述实际地质体的各向异性介质(Fedorov, 1968; Auld, 1973; Winterstein, 1990; Helbig, 1994).TI介质是一种关于对称轴无限对称的各向异性介质，它是地壳中最常见的各向异性介质类型.上地壳定向排列的微裂隙被认为是上地壳各向异性的起源(Crampin, 1984; Barruol and Kern, 1996; Crampin and Zatsepin, 1997; Zatsepin and Crampin, 1997)，它表现出的各向异性性质为HTI介质，即对称轴平行于地面的TI介质.下地壳及沉积盆地的各向异性被认为是由水平叠加的薄层引起的(Backus, 1962; Wild and Crampin, 1991)，它表现出的各向异性性质为VTI介质，即对称轴垂直于地面的TI介质.TTI介质(倾斜对称轴)是TI介质的一般形式，它比各向同性介质及VTI介质更加接近真实的地层结构.因此三维TTI介质射线追踪一直都是国内外地震领域中的热点研究方向之一.
射线追踪方法可以追踪模型中的地震波射线路径及走时信息，如试射法(Červený, 1972; Langan et al., 1985; Rawlinson et al., 2001; Xu et al., 2008; 徐涛等，2004)、弯曲法(Julian and Gubbins, 1977; Um and Thurber, 1987; Xu et al., 2006, 2010, 2014)、最短路径法(Nakanishi and Yamaguchi, 1986; Moser, 1991; Bai et al., 2009)、波前法(Vinje et al., 1993)等，以及基于程函方程数值解获得走时场，继而追踪路径法，如有限差分法(Vidale, 1988)、快速推进法(Sethian and Popovici, 1999; Rawlinson and Sambridge, 2004; De Kool et al., 2006)、快速扫描法(Zhao, 2005, 2007)、线性走时插值法(Asakawa and Kawanaka, 1993；马德堂和朱光明，2006；张东等，2009)等.
高斯运动学射线追踪(Červený, 1972, 2001)可以追踪任意弹性介质中的体波射线路径及走时，且在临界区及阴影区等特殊情况下算法稳定性好，国内外对其进行了许多研究(Hanyga, 1986；Zhu et al., 2005; 邓飞和刘超颖, 2009；孙成禹等，2011).因为高斯运动学射线追踪快速精确的特点，它被广泛应用于地震波层析成像(Chapman and Pratt, 1992; Pratt and Chapman, 1992)、高斯偏移成像(Červený, 1985; Hill, 1990; Han and Wang, 2015)、理论地震图计算(Thomson et al., 1992; Guest and Kendall, 1993; Alkhalifah, 1995)及震源定位(曹雷等, 2015)等领域中.三维TTI介质中传统的运动学射线追踪通常需要利用Bond变换或者四阶张量变换求取TTI介质中的弹性参数，由于TTI介质的21个弹性参数绝大部分情况下都非零，所以只能使用最复杂的各向异性射线追踪方程.这样的处理方式相当于将TTI介质当成最复杂的三斜晶系各向异性介质对待，从而产生了许多额外的计算量.
TTI介质与VTI介质唯一的区别在于对称轴方向有所不同，因此实际野外对称轴方向统一的复杂TTI介质模型可以看成VTI介质模型旋转一定角度获得.为此，如果首先在本构坐标系中进行VTI射线追踪，然后再通过坐标旋转将射线路径旋转至观测坐标系中，获得TTI介质射线路径，将会显著提高射线追踪的效率.因为这样避免了使用Bond变换及最复杂的各向异性追踪方程.本文拟用该方法研究单层及多层TTI介质中的P波，SV波及SH波的波前面及射线路径.
走时曲线是研究地震波射线传播规律非常重要的手段之一.目前在许多地震勘探中仍然使用水平层状VTI介质走时曲线来解释倾斜地层的TTI介质，这样会使得最终地质解释与实际的地层结构产生很大的偏差.因此本文还简单地分析了多层VTI与TTI介质中P波走时曲线的差别.
1 理论方法 1.1 TTI介质传统射线追踪方法 1.1.1 射线追踪方程根据爱因斯坦求和约定，任意弹性介质中的射线追踪方程可以表述如下(Červený, 1972，2001)：
|
(1) |
其中gjgk=Djk/D，Djk及D的表达式为：
|
(2) |
其中Γik=aijklpjpl，Γ为Christoffel矩阵.xi为射线当前位置在i方向的分量，pj为相慢度在j方向的分量，aijkl=cijkl/ρ，即密度归一化的弹性参数.该射线追踪方程对于qP波射线追踪始终是正确的，qS波射线追踪的情况要复杂许多.在临界区等特殊情况下，qS1与qS2会部分耦合，此时会导致射线追踪方程失效.此外，如果D=0时，有Djk/D=0/0，此时射线追踪方程也会失效(Shearer and Chapman, 1989; Pšeník and Dellinger, 2001).对于临界区的qS波射线追踪，需要加入一些控制条件来保证其正确性(Vavrycuk, 2000).
1.1.2 Bond变换观测系统xyz与本构系统x0y0z0(x0, y0, z0分别对应射线追踪方程中的x1, x2, x3)的关系如图 1所示，本构坐标系中对称轴为z0轴.因为TI介质是关于对称轴无限对称的介质，所以我们将y0轴设定在xoy平面内.相慢度矢量p与z0轴的夹角为Φ，即相角为Φ.相慢度矢量p与z轴的夹角为α，在xoy平面内的投影与x轴的夹角为β.从x0y0z0坐标系到xyz坐标系经历了两次坐标旋转.第一次旋转是以y0轴为旋转轴旋转了θ，x0轴变为了x′轴，z0轴变为了z轴.第二次是以z轴为旋转轴旋转了φ，x′轴变为了x轴，y0轴变为了y轴.从xyz坐标系到x0y0z0坐标系是其逆过程.
|
图 1 观测系统与本构系统的关系 Fig. 1 Relationship between the observation system and constitutive system |
由x0y0z0坐标系到xyz坐标系的旋转矩阵为：
|
(3) |
由xyz坐标系到x0y0z0坐标系的旋转矩阵为：
|
(4) |
对于TI介质，以z0轴为对称轴的本构坐标系下密度归一化弹性矩阵的Voigt形式为：
|
(5) |
本构坐标系下密度归一化弹性矩阵A到观测坐标系下TTI介质的密度归一化弹性矩阵Axyz可由Bond变换(Bond，1943)得到：Axyz=αAαT.利用四阶张量变换(姚陈和蔡明刚，2009)同样可以得到Axyz.由(3)式可求得Bond变换矩阵α的具体表达式：
|
(6) |
将α, αT代入Axyz=αAαT中即可求得TTI介质(φ, θ都取任意值)中的密度归一化弹性矩阵.因为TTI介质大部分情况21个弹性参数都非零，因此只能使用最复杂的各向异性追踪方程(1)进行射线追踪.
1.2 TTI介质高效射线追踪方法本文推导了三维VTI介质射线追踪方程，提出先在本构坐标系中进行VTI射线追踪，然后再通过坐标旋转将射线路径旋转至观测坐标系中，获得TTI介质射线路径.
1.2.1 三维VTI射线追踪方程Červený(2001)给出了二维VTI射线追踪方程，但是该方程无法应用到三维VTI介质中.我们将密度归一化的VTI弹性矩阵代入Christoffel矩阵及各向异性追踪方程(1)中，得出三维VTI射线追踪方程：
|
(7) |
|
(8) |
Christoffel矩阵分量Γik同样可以得到简化：
|
(9) |
我们在三维VTI射线追踪方程的推导过程中没有使用任何近似，同时旋转矩阵也是精确表达式.因此本构系统中的三维VTI射线追踪与观测系统中的TTI射线追踪在理论上结果应该是完全一致，不存在任何误差.
1.2.2 三维TTI高效射线追踪相较于传统方法，本文选择的方法是先进行VTI射线追踪，然后再通过VTI介质射线路径获得TTI介质射线路径.该方法的基本原理是：
(ⅰ)利用旋转矩阵(4)式将观测坐标系xyz中的相慢度矢量变换到本构坐标系x0y0z0中；
(ⅱ)在本构坐标系x0y0z0中使用VTI介质射线追踪方程(7)—(9)进行射线追踪获得本构坐标系中的射线路径以及走时信息；
(ⅲ)为了方便进行坐标变换，将震源点位置平移至坐标原点处，追踪所得的射线路径同样进行相应地平移；
(ⅳ)利用对称轴方向(φ, θ)以及旋转矩阵(3)式将本构系统中的VTI介质射线路径旋转成为观测系统中的TTI介质射线路径；
(ⅴ)重新将震源点平移回原来的位置，将旋转后所得的TTI介质射线路径进行相应地平移.因为走时不受坐标变换的影响，所以观测坐标系中的走时与本构坐标系中的走时保持一致.
在射线追踪过程中，射线每走一步都需要重复上述过程，这样可以确保结果的正确性.利用VTI介质射线追踪来获取TTI介质射线路径的方法可以极大地提高射线追踪的效率.我们从理论角度来分析对比传统方法与本文方法计算量的大小：①传统方法需要利用Bond变换或者四阶张量变换将所有采样点的VTI弹性矩阵变换为TTI弹性矩阵.相比较而言，本文的方法避免了使用Bond变换，即省略了这一步骤；②本文利用VTI射线追踪方程，大大简化了任意各向异性介质射线追踪方程：以x方向为例，任意各向异性介质运动学射线追踪方程中每一步∂x/∂T的值需要计算27项，而使用VTI射线追踪方程每一步∂x/∂T的值只需要计算7项.从公式的简化程度上可以发现，本文方法在计算∂x/∂T的效率上比传统方法提高了约3.86倍左右.同理，传统方法中∂px/∂T需要计算81项，本文方法中∂px/∂T只需要计算21项，也就是说本文方法在计算∂px/∂T的效率上比传统方法也提高约3.86倍左右.另外，考虑到Christoffel矩阵的对称性，Γik的计算量由原来的54项减为15项，效率提高了3.6倍；③本文方法多出的步骤是将本构坐标系中的VTI射线路径变换为观测坐标系中的TTI射线路径，但是这一步的计算量很小.
1.2.3 初始条件对于各向异性介质的射线追踪，给定的初始条件为：震源点的位置xi=xi0及初始相慢度pi=pi0，它决定了波前面的法线方向，而非射线出射方向.利用Thomsen参数VP0, VS0, δ, ε, γ, δ*, D*及相角Φ，可以将各向异性介质中的精确相速度公式(Thomsen, 1986)表达为：
|
(10) |
其中Thomsen参数与密度归一化弹性参数之间的转换关系为：
|
(11) |
参数δ*, D*与VP0, VS0, δ, ε, γ的关系为：
|
(12) |
如图 1所示，利用精确相速度公式(10)及相角Φ可以求得相速度V.相慢度分量pi可由相速度V以及初始相慢度方向(α, β)求出：
|
(13) |
虽然射线追踪方程使用的初始条件是相慢度，方向为波前面的法线方向.但是追踪过程中已经将相慢度分量转换为群速度分量∂xi/∂T=aijklplgjgk，所以最后追踪出来的xi是射线的路径.从震源点出发的各个方向的射线在同一时刻到达的位置是一个等相位面，即波前面.因此波前面的形状应该由射线的速度所决定，即波前面形状应该与群速度面形状一致，而非相速度面.
1.2.4 微分方程解法及参数插值本文求解微分方程组(7)式与(8)式所采用的方法是四阶龙格库塔法.令∂xi/∂T=f(xi, pi), ∂pi/∂T=φ(xi, pi)，则四阶龙格库塔的斜率可以表示为：
|
(14) |
其中dt为时间采样间隔.利用已知的(xi, pi)通过(7)式与(8)式可以依次求出对应的四阶龙格库塔的斜率.然后通过求得的四阶龙格库塔斜率可以将微分方程组转换为差分方程组：
|
(15) |
利用差分方程组(15)式以及初始条件，可以从震源点开始逐步进行射线追踪，获取射线路径以及走时信息.
本文建立速度模型时，采样点位置处的各向异性参数是初始给定的.采样点之外其他位置的各向异性参数是通过对包围该点的立方体八个顶点进行拉格朗日插值获得.假设插值点的坐标为(x, y, z)，需要进行插值的参数为δ.插值点坐标(x, y, z)处于采样点号(xi, yi, zi)与(xi+1, yi+1, zi+1)之间(采样点号乘以对应方向的采样间隔才是该采样点的坐标)，(x, y, z)方向的采样间隔分别为(dx, dy, dz).那么通过包围该点的立方体八个顶点处的空间坐标，可以求出八个顶点处的插值基函数：
|
(16) |
通过求得的八个顶点处插值基函数Ri以及对应的δi参数可以插值出(x, y, z)处的δ参数：δ= 
本文采用的模型参数化方法为常规的三维节点法(Thurber, 1983).如果对于同一介质，本文方法与传统方法计算结果相同.那么这就说明本文的方法是可行的.在数值模拟实验中，我们建立了大量的速度模型，并发现两种方法计算所得的结果均完全一致，证明了两种方法理论上没有误差的结论.这里只给出一个简单单层模型中P波的计算结果作为示例.模型的各向异性参数δ, ε, γ均取0.2，VP0, VS0分别取6 km·s-1与3 km·s-1，对称轴的φ取0°，θ取45°.首先，我们将xyz三个方向的采样点都设为10个，即模型总采样点数为103个，采样间隔为1 km.两种方法计算所得0.6 s时刻的P波波前面如图 2所示.
|
图 2 (a) 传统方法P波波前面, (b)本文方法P波波前面 Fig. 2 (a) Wavefronts of P waves calculated by traditional method, (b) Wavefronts of P waves calculated by our method |
由图 2可以看出两种方法计算的波前面完全一致.其中图 2a是使用传统方法计算所得的波前面，除去读入文件等预处理时间，CPU(Intel Core i7 2.50 GHz)耗时6.49 s.图 2b使用的是本文方法计算所得的波前面，CPU耗时1.66 s.从数值计算角度，可以发现本文方法的效率是传统方法的3.91倍左右.该结果略高于理论公式分析角度得出的3.86倍.这是因为传统方法在Bond变换或者四阶张量变换上需要消耗额外的时间.当模型采样点数增加时，Bond变换耗时同样会增加，而本文方法的坐标变换耗时增加相对较少.为此，我们将模型网格采样点加密，将xyz三个方向的网格采样点都增加到100个，即模型总采样点数为1003个，采样间隔相应地降为0.1 km.对于该模型而言，因为各向异性参数简单并且统一，所以两种网格疏密程度不同的模型都为精确描述，波前面结果也完全一致.但是密网格模型的传统方法CPU耗时变为12.07 s，本文方法的CPU耗时变为2.64 s，此时本文方法的效率是传统方法的4.57倍.因此对于一些复杂密网格模型而言，本文方法相对于传统方法具有更高的效率.
2.2 模型试验 2.2.1 单层模型为检验本文方法的效果，我们建立了多个单层TTI介质模型(表 1)来进行数值模拟分析.首先，通过VTI射线追踪方程获得本构坐标系下的P波，SV波及SH波的波前面及射线路径坐标.然后通过坐标旋转得到观测系统中的TTI介质射线路径.
|
表 1 单层TTI模型各向异性参数 Table 1 Anisotropic parameters of single-layer TTI models |
射线路径是同一射线在不同时间到达不同位置的轨迹，波前面是不同方向出射的射线在同一时间到达的位置组成的等相位面.因此类似于波前面的定义，我们刻画波前面的方法是每间隔1°从震源(5 km, 5 km, -5 km)追踪一条射线.即初始相慢度的方向α从0°到180°每隔1°取一个值，β从0°到359°每隔1°取一个值，t时刻的射线位置就是t时刻的波前面.
图 3a—3c分别给出了M1模型在本构坐标系中的P波0.5 s时刻以及SV波，SH波1 s时刻的波前面.图 3d—3f分别给出了P波波前面在xoz, yoz, xoy面的投影.图 3a—3c中可以看出，P波波前面近似一个椭球，SV波波前面复杂一些，而SH波波前面是一个标准椭球.本文所有图中的波前面颜色都设定为随深度的变化而变化，这是为了使得三维图像看起来更有立体感.图中蓝色的粗线是对称轴，黑色的细线就是射线，我们每隔30°相慢度方向追踪一条射线路径.虽然初始相慢度方向是均匀给定的(每30°出射一条射线)，但是图 3d—3f中可以看出xoz, yoz面射线出射方向的分布并不均匀.这是因为TI介质在纵向上群速度与相速度不一致引起的.因为TI介质在横向上是各向同性的，所以xoy面的射线是均匀分布的.P波波前面在xoz, yoz面的投影形状完全一致，在xoy面的投影是一个圆.该结果与TI介质关于对称轴无限对称的结论一致.
|
图 3 M1模型本构坐标系中射线路径及波前面 (a) P波; (b) SV波; (c) SH波. P波射线路径及波前面在(d)xoz面的投影; (e) yoz面的投影; (f) xoy面的投影. Fig. 3 Ray paths and wavefronts of M1 model in constitutive system (a) P waves; (b) SV waves; (c) SH waves. Ray paths and wavefronts of P waves projecting onto (d) xoz plane; (e) yoz plane; (f) xoy plane. |
利用坐标旋转，我们得到了M1模型观测系统中的P波、SV波、SH波的波前面以及P波波前面在xoz, yoz, xoy面的投影(图 4).从图 3a—3c与图 4a—4c的对比可以看出射线路径与波前面的确从本构坐标系中的VTI图形变为了观测坐标系中的TTI图形.利用对称轴取向(φ, θ)，可以求得对称轴在xoz面的投影与z轴的夹角为arctan(cosφtanθ)，在yoz面的投影与z轴的夹角为arctan(sinφtanθ)，在xoy面的投影与x轴的夹角为φ.M1模型对称轴的方向(φ, θ)取值为(0°, 45°)，因此对称轴在xoz面及yoz面的投影与z轴的夹角分别为45°(图 4d)及0°(图 4e)，对称轴在xoy面的投影与x轴的夹角为0°(图 4f).如图 4所示，随着对称轴方向的变化，波前面的形状也发生了相应的变化，但始终相对于对称轴呈对称分布.
|
图 4 M1模型观测坐标系中射线路径及波前面 (a) P波, (b) SV波, (c) SH波. P波射线路径及波前面在(d)xoz面的投影; (e)yoz面的投影; (f)xoy面的投影. Fig. 4 Ray paths and wavefronts of M1 model in observation system (a) P waves, (b) SV waves, (c) SH waves. Ray paths and wavefronts of P waves projecting onto (d) xoz plane; (e) yoz plane; (f) xoy plane. |
为了避免重复，后续其他模型都只给出观测坐标系下的波前面结果图.M2模型的特点是δ=ε，此时P波的波前面是一个标准的椭球面，即椭圆各向异性(Daley and Hron, 1979)，SV波的相速度恒等于VS0，因此SV波的波前面类似于各向同性介质中的形状，是一个标准的球面(Thomsen, 1986).图 5a—5c分别给出了P波0.5 s时刻及SV与SH波1 s时刻的波前面图.从图中可以看出P波及SV波波前面与Thomsen (1986)的结论一致，一个是标准椭球面(图 5a)，一个是标准球面(图 5b)，SH波的波前面始终是一个标准椭球面(图 5c).图 5d—5f分别给出了P波波前面在xoz, yoz, xoy面的投影.从图中可以看出，虽然对称轴在各方向的投影不同，但是波前面投影一定相对于对称轴投影呈对称分布.
|
图 5 M2模型观测坐标系中射线路径及波前面 (a) P波, (b) SV波, (c) SH波. P波射线路径及波前面在(d)xoz面的投影; (e) yoz面的投影; (f) xoy面的投影. Fig. 5 Ray paths and wavefronts of M2 model in observation system (a) P waves, (b) SV waves, (c) SH waves. Ray paths and wavefronts of P waves projecting onto (d) xoz plane; (e) yoz plane; (f) xoy plane. |
M3与M4是两个比较特殊的模型，M3的特点是δ为负数，M4的特点是各向异性强度非常大.图 6给出了M3中P波0.5 s时刻、SV波0.8 s时刻及SH波1 s时刻的波前面.图 7给出了M4中P波0.5 s时刻、SV波1.2 s时刻及SH波1 s时刻的波前面.在这两种复杂各向异性参数模型中，P波的波前面不再类似于椭球面，不同方向的群速度变化相对于M1与M2模型更加复杂(图 6a, 7a).SV波的波前面最为复杂，形成了TI介质中特有的线形奇异值(图 6b, 7b)，从而引起了多值走时(Crampin and Yedlin, 1981; Crampin, 1991).SH波的波前面最为简单，始终是标准椭球面，相对于M1与M2模型不同的是椭球扁率变的非常大(图 6c, 7c).本文方法射线追踪所得的波前面与梁锴(2009)计算所得理论各向异性介质中的群速度面完全一致.这说明本文方法在复杂各向异性参数模型中的计算结果仍然是完全精确的，与传统方法不存在任何误差.
|
图 6 M3模型观测坐标系中射线路径及波前面 (a) P波, (b) SV波, (c) SH波. Fig. 6 Ray paths and wavefronts of M3 model in observation system (a) P waves; (b) SV waves; (c) SH waves. |
|
图 7 M4模型观测坐标系中射线路径及波前面 (a) P波, (b) SV波, (c) SH波. Fig. 7 Ray paths and wavefronts of M4 model in observation system (a) P waves; (b) SV waves; (c) SH waves. |
对于多层介质而言，分界面的形状不会影响坐标变换，因此本文方法可以应用到多层弯曲分界面介质中.为检验本文方法适应多层弯曲分界面模型的能力，我们建立了一个多层且具有倾斜弯曲分界面的TTI介质模型M5(表 2)，模型中分界面使用二维三次B样条函数定义.如表 2所示，M5由三层不同各向异性参数的介质所组成，对称轴取向(φ, θ)统一设定为(90°, 45°).震源位置设置在(5 km, 5 km, -5 km)，位于第二层介质中.
|
|
表 2 M5模型各向异性参数 Table 2 Anisotropic parameters of model M5 |
图 8给出了M5中P波0.6 s时刻、SV波及SH波1 s时刻的波前面，该结果与传统方法结果仍然一致.需要注意的是，射线经过分界面时，通过TI介质中的Snell定律获取折射波(Daley and Hron, 1977; Slawinski et al., 2000).因此相比于单层介质，多层介质需要消耗额外的CPU时间来搜索折射点及折射角.这个过程对于传统方法及本文方法的耗时是一致的，这样就会导致多层介质中本文方法的效率相对于单层介质稍微有所减少.如图 8所示，P、SV及SH波从中间层传播到第一层及第三层时，波前面有明显的折射现象.但是由于分界面形状变化比较平缓，所以第一层及第三层的波前面比较规则.对于具有复杂分界面的TTI介质模型，其处理过程与简单分界面的介质模型一致，并且在CPU耗时上没有明显区别.
|
图 8 M5模型观测坐标系中射线路径及波前面 (a) P波; (b) SV波; (c) SH波. Fig. 8 Ray paths and wavefronts of M5 model in observation system (a) P waves; (b) SV waves; (c) SH waves. |
走时曲线是研究地震波射线传播规律非常重要的手段之一.目前在许多地震勘探中仍然使用水平层状VTI介质走时曲线来解释倾斜地层的TTI介质，这样会使得最终地质解释与实际的地层结构产生很大的偏差.因此研究VTI介质走时曲线与TTI介质走时曲线的异同点有很大的应用价值.因为P波在各向异性介质中不像S波一样发生分裂，并且从M1—M4模型的对比可以发现P波波前面相对于SV波波前面要简单得多，所以在地震勘探中主要使用的是P波.因此本文简单地分析了多层VTI与TTI介质中P波走时曲线的差别.
走时曲线不仅会受到模型参数的影响，还会受到分界面形状的影响.为了减少分界面形状的影响，突出对称轴方向及地层倾斜程度对走时曲线的影响，我们将分界面设为平面.建立了多层模型M6 (表 3)，在M6中每层的各向异性参数都不一致，并且沿着对称轴方向的P波与S波速度随着深度的增加而增加.这里没有给出对称轴方向的取值，因为本文会对不同对称轴方向的走时曲线进行对比.M5中震源位置设置在(15 km, 15 km, -27 km)，始终位于模型的最底层.
|
|
表 3 M6模型各向异性参数 Table 3 Anisotropic parameters of M6 model |
如图 9所示，本文给出了M6模型TTI介质(此处取φ=90°, θ=15°作为示例)中P波在xoz截面与yoz截面的波前面及射线路径图.图中蓝色的点代表波前面，蓝色实线代表对称轴，红色的实线代表分界面，黑色的实线为射线.由M1—M4模型的波前面分析可知，单层介质中波前面在任意方向的投影仍然会相对于对称轴呈对称分布.在图 9中可以发现，该规律对于简单的多层TTI介质同样适用.VTI介质的对称轴垂直于地表，所以沿着任意方向测量所得到的走时曲线都一致，并且相对于震源点的投影呈对称分布.但是TTI介质的对称轴是倾斜的，其不同切面方向的走时曲线会有所差别.其中平行于对称轴地面投影方向(图 9b)的走时曲线受到对称轴倾斜角的影响最大.而在垂直于对称轴地面投影方向(图 9a)的走时曲线应该与VTI介质的走时曲线一致.
|
图 9 M6模型P波射线路径及波前面 (a) xoz截面; (b) yoz截面. Fig. 9 Ray paths and wavefronts of P waves in M6 model (a) Projection onto xoz plane; (b) Projection onto yoz plane. |
因为平行于对称轴地面投影方向的走时曲线受到对称轴倾斜角的影响最大，所以我们只给出了该方向(yoz截面)在不同倾斜程度下(θ取值0°(VTI)，15°, 30°, 45°, 60°，75°)的走时曲线(图 10).从图中可以看出，VTI介质中的走时曲线的确是关于震源点在地面的投影呈对称分布的.但是当VTI对称轴开始倾斜成TTI介质时，走时曲线开始慢慢变得不对称，并且变得复杂起来.因为垂直于对称轴方向的P波群速度达到最大，所以随着对称轴倾斜程度越大，射线到达地面所需要的走时就会越小.这在图 10中表现为θ越大，走时最小值越小，并且最小值点不再处于震源的地面投影处.从图 10中的走时曲线对比可知，使用VTI介质的走时曲线去研究倾斜地层的TTI介质会产生很大的偏差.
|
图 10 M6模型平行于对称轴投影方向走时曲线随倾斜角变化 Fig. 10 Traveltime curves versus varied dips of M6 model in the direction of y axis |
本文利用VTI介质弹性矩阵对任意介质运动学射线追踪方程进行简化，得到了三维VTI介质运动学射线追踪方程.然后利用三维VTI介质运动学射线追踪方程研究了三维TTI介质中的体波射线路径及波前面.相比于基于Bond变换及任意各向异性介质运动学射线追踪方程的传统方法，本文提出先在本构坐标系中进行VTI介质射线追踪，再通过坐标旋转获得观测系统中的TTI介质射线路径.通过理论分析及数值计算，证实了本文方法的正确性及高效性.本文方法要求多层介质的对称轴必须一致.如果对称轴不一致，不同采样点之间的各向异性参数插值结果会不正确，未来工作中，我们将进一步改进该方法的不足.
本文通过建模分析发现：(1)体波的波前面在强各向异性及各种特殊情况下都与理论群速度面一致；(2)波前面在某个面的投影始终相对于对称轴的投影呈对称分布；(3)给定的初始条件为相慢度的方向，而非射线的方向.因此虽然相慢度方向是均匀给定的，但是追踪出来的射线方向分布不是均匀的；(4) VTI介质的走时曲线相对于震源位置呈对称分布且相对较简单.TTI介质的走时曲线不再相对于震源位置呈对称分布，并且比VTI介质的走时曲线复杂；(5)因为垂直于对称轴方向的群速度达到最大，所以随着TTI介质对称轴倾斜程度越大，射线到达地面所需要的走时就会越小.
Alkhalifah T.
1995. Efficient synthetic-seismogram generation in transversely isotropic, inhomogeneous media. Geophysics, 60(4): 1139-1150.
DOI:10.1190/1.1443842 |
|
Asakawa E, Kawanaka T.
1993. Seismic ray tracing using linear traveltime interpolation. Geophysical Prospecting, 41(1): 99-111.
DOI:10.1111/gpr.1993.41.issue-1 |
|
Auld B A. 1973. Acoustic resonators. //Acoustic Fields and Waves in Solids. 2nd ed. New York: John Wiley & Sons.
|
|
Backus G E.
1962. Long-wave elastic anisotropy produced by horizontal layering. Journal of Geophysical Research, 67(11): 4427-4440.
DOI:10.1029/JZ067i011p04427 |
|
Bai C Y, Tang X P, Zhao R.
2009. 2-D/3-D multiply transmitted, converted and reflected arrivals in complex layered media with the modified shortest path method. Geophysical Journal International, 179(1): 201-214.
DOI:10.1111/gji.2009.179.issue-1 |
|
Barruol G, Kern H.
1996. Seismic anisotropy and shear-wave splitting in lower-crustal and upper-mantle rocks from the Ivrea Zone-experimental and calculated data. Physics of the Earth and Planetary Interiors, 95(3-4): 175-194.
DOI:10.1016/0031-9201(95)03124-3 |
|
Bond W L.
1943. The mathematics of the physical properties of crystals. Bell System Technical Journal, 22(1): 1-72.
DOI:10.1002/bltj.1943.22.issue-1 |
|
Cao L, Zhang J H, Yao Z X.
2015. Source location using 3-D Gaussian beam migration imaging. Chinese Journal of Geophysics (in Chinese), 58(2): 481-494.
DOI:10.6038/cjg20150212 |
|
Červený V.
1972. Seismic rays and ray intensities in inhomogeneous anisotropic media. Geophysical Journal International, 29(1): 1-13.
DOI:10.1111/j.1365-246X.1972.tb06147.x |
|
Červený V.
1985. Gaussian beam synthetic seismograms. Journal of Geophysics Zeitschrift Fur Geophysik, 58: 44-72.
|
|
Červený V. 2001.
Seismic Ray Theory. Cambridge: Cambridge University Press.
|
|
Chapman C H, Pratt R G.
1992. Traveltime tomography in anisotropic media-I. Theory. Geophysical Journal International, 109(1): 1-19.
DOI:10.1111/gji.1992.109.issue-1 |
|
Crampin S, Yedlin M.
1981. Shear-wave singularities of wave propagation in anisotropic media. Journal of Geophysics, 49: 43-46.
|
|
Crampin S.
1984. Effective anisotropic elastic constants for wave propagation through cracked solids. Geophysical Journal International, 76(1): 135-145.
DOI:10.1111/j.1365-246X.1984.tb05029.x |
|
Crampin S.
1991. Effects of point singularities on shear-wave propagation in sedimentary basins. Geophysical Journal International, 107(3): 531-543.
DOI:10.1111/gji.1991.107.issue-3 |
|
Crampin S, Zatsepin S V.
1997. Modelling the compliance of crustal rock:Ⅱ-response to temporal changes before earthquakes. Geophysical Journal International, 129(3): 495-506.
DOI:10.1111/gji.1997.129.issue-3 |
|
Daley P F, Hron F.
1977. Reflection and transmission coefficients for transversely isotropic media. Bulletin of the Seismological Society of America, 67(3): 661-675.
|
|
Daley P F, Hron F.
1979. Reflection and transmission coefficients for seismic waves in ellipsoidally anisotropic media. Geophysics, 44(1): 27-38.
DOI:10.1190/1.1440920 |
|
De Kool M, Rawlinson N, Sambridge M.
2006. A practical grid-based method for tracking multiple refraction and reflection phases in three-dimensional heterogeneous media. Geophysical Journal International, 167(1): 253-270.
DOI:10.1111/gji.2006.167.issue-1 |
|
Deng F, Liu C Y.
2009. 3-D rapid ray-tracing and Gaussian ray-beam forward simulation. Oil Geophysical Prospecting (in Chinese), 44(2): 158-165.
|
|
Fedorov F I. 1968.
Theory of Elastic Waves in Crystals. New York: Plenum Press.
|
|
Guest W S, Kendall J M.
1993. Modelling seismic waveforms in anisotropic inhomogeneous media using ray and Maslov asymptotic theory:applications to exploration seismology. Canadian Journal of Exploration Geophysics, 29: 78-92.
|
|
Han J G, Wang Y.
2015. Gaussian beam prestack depth migration in heterogeneous transversely isotropic media. Exploration Geophysics, 46(2): 153-158.
DOI:10.1071/EG13061 |
|
Hanyga A.
1986. Gaussian beams in anisotropic elastic media. Geophysical Journal International, 85(3): 473-504.
DOI:10.1111/j.1365-246X.1986.tb04528.x |
|
Helbig K. 1994.
Foundations of Anisotropy for Exploration Seismics. Oxford: Elsevier Science Ltd.
|
|
Hill N R.
1990. Gaussian beam migration. Geophysics, 55(11): 1416-1428.
DOI:10.1190/1.1442788 |
|
Julian B R, Gubbins D.
1977. Three-dimensional seismic ray tracing. Journal of Geophysics, 43: 95-113.
|
|
Langan R T, Lerche I, Cutler R T.
1985. Tracing of rays through heterogeneous media:an accurate and efficient procedure. Geophysics, 50(9): 1456-1465.
DOI:10.1190/1.1442013 |
|
Liang K. 2009. The study on propagation feature and forward modeling of seismic wave in TI media (in Chinese). Beijing: China University of Petroleum.
|
|
Ma D T, Zhu G M.
2006. Computation of traveltimes of seismic first breaks in transversely isotropic medium. Oil Geophysical Prospecting (in Chinese), 41(1): 26-31.
|
|
Moser T J.
1991. Shortest path calculation of seismic rays. Geophysics, 56(1): 59-67.
DOI:10.1190/1.1442958 |
|
Nakanishi I, Yamaguchi K.
1986. A numerical experiment on nonlinear image reconstruction from first-arrival times for two-dimensional island arc structure. Journal of Physics of the Earth, 34(2): 195-201.
DOI:10.4294/jpe1952.34.195 |
|
Pratt R G, Chapman C H.
1992. Traveltime tomography in anisotropic media-Ⅱ. Application. Geophysical Journal International, 109(1): 20-37.
DOI:10.1111/gji.1992.109.issue-1 |
|
Pšeník I, Dellinger J A.
2001. Quasi-shear waves in inhomogeneous weakly anisotropic media by the quasi-isotropic approach:A model study. Geophysics, 66(1): 308-319.
DOI:10.1190/1.1444909 |
|
Rawlinson N, Houseman G A, Collins C D N.
2001. Inversion of seismic refraction and wide-angle reflection traveltimes for three-dimensional layered crustal structure. Geophysical Journal International, 145(2): 381-400.
DOI:10.1046/j.1365-246x.2001.01383.x |
|
Rawlinson N, Sambridge M.
2004. Multiple reflection and transmission phases in complex layered media using a multistage fast marching method. Geophysics, 69(5): 1338-1350.
DOI:10.1190/1.1801950 |
|
Sethian J A, Popovici A M.
1999. 3-D traveltime computation using the fast marching method. Geophysics, 64(2): 516-523.
DOI:10.1190/1.1444558 |
|
Shearer P M, Chapman C H.
1989. Ray tracing in azimuthally anisotropic media-I. results for models of aligned cracks in the upper crust. Geophysical Journal International, 96(1): 51-64.
DOI:10.1111/gji.1989.96.issue-1 |
|
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 |
|
Sun C Y, Zhang W Y, Ni C K, et al.
2011. Seismic wave field modeling by energy controlled Gaussian beam method. Oil Geophysical Prospecting (in Chinese), 46(6): 856-861.
|
|
Thomsen L.
1986. Weak elastic anisotropy. Geophysics, 51(10): 1954-1966.
DOI:10.1190/1.1442051 |
|
Thomson C J, Kendall J M, Guest W S.
1992. Geometrical theory of shear-wave splitting:Corrections to ray theory for interference in isotropic/anisotropic transitions. Geophysical Journal International, 108(1): 339-363.
DOI:10.1111/gji.1992.108.issue-1 |
|
Thurber C H.
1983. Earthquake locations and three-dimensional crustal structure in the Coyote Lake Area, central California. Journal of Geophysical Research:Solid Earth, 88(B10): 8226-8236.
DOI:10.1029/JB088iB10p08226 |
|
Um J, Thurber C H.
1987. A fast algorithm for two-point seismic ray tracing. Bulletin of the Seismological Society of America, 77(3): 972-986.
|
|
Vavrycuk V. 2000. Ray tracing for anisotropic media in singular directions of S waves. //Seismic Waves in Complex 3-D Structures. Prague: Faculty of Mathematics and Physics, Charles University, 161-183.
http://sw3d.cz/papers/r10vv1.htm |
|
Vidale J.
1988. Finite difference calculation of travel times. Bulletin of the Seismological Society of America, 78(6): 2062-2076.
|
|
Vinje V, Iverson E, Gjøystdal H.
1993. Traveltime and amplitude estimation using wavefront construction. Geophysics, 58(8): 1157-1166.
DOI:10.1190/1.1443499 |
|
Wild P, Crampin S.
1991. The range of effects of azimuthal isotropy and EDA anisotropy in sedimentary basins. Geophysical Journal International, 107(3): 513-529.
DOI:10.1111/gji.1991.107.issue-3 |
|
Winterstein D F.
1990. Velocity anisotropy terminology for geophysicists. Geophysics, 55(8): 1070-1088.
DOI:10.1190/1.1442919 |
|
Xu T, Xu G M, Gao E G, et al.
2004. Block modeling and shooting ray tracing in complex 3-D media. Chinese Journal of Geophysics (in Chinese), 47(6): 1118-1126.
|
|
Xu T, Xu G M, Gao E G, et al.
2006. Block modeling and segmentally iterative ray tracing in complex 3D media. Geophysics, 71(3): T41-T51.
DOI:10.1190/1.2192948 |
|
Xu T, Zhang Z, Zhao A, et al.
2008. Sub-triangle shooting ray-tracing in complex 3D VTI media. Journal of Seismic Exploration, 17(2): 133-146.
|
|
Xu T, Zhang Z J, Gao E G, et al.
2010. Segmentally iterative ray tracing in complex 2D and 3D heterogeneous block models. Bulletin of the Seismological Society of America, 100(2): 841-850.
DOI:10.1785/0120090155 |
|
Xu T, Li F, Wu Z B, et al.
2014. A successive three-point perturbation method for fast ray tracing in complex 2D and 3D geological models. Tectonophysics, 627: 72-81.
DOI:10.1016/j.tecto.2014.02.012 |
|
Yao C, Cai M G.
2009. Analytical expression of TI elastic tensor with arbitrary orientation. Chinese Journal of Geophysics (in Chinese), 52(9): 2345-2348.
DOI:10.3969/j.issn.0001-5733.2009.09.019 |
|
Zatsepin S V, Crampin S.
1997. Modelling the compliance of crustal rock:I-response of shear-wave splitting to differential stress. Geophysical Journal International, 129(3): 477-494.
DOI:10.1111/gji.1997.129.issue-3 |
|
Zhang D, Fu X R, Yang Y, et al.
2009. 3-D seismic ray tracing algorithm based on LTI and partition of grid interface. Chinese Journal of Geophysics (in Chinese), 52(9): 2370-2376.
DOI:10.3969/j.issn.0001-5733.2009.09.023 |
|
Zhao H K.
2005. A fast sweeping method for Eikonal equations. Mathematics of Computation, 74(250): 603-627.
|
|
Zhao H.
2007. Parallel implementations of the fast sweeping method. Journal of Computational Mathematics, 25(4): 421-429.
|
|
Zhu T F, Gray S, Wang D L. 2005. Kinematic and dynamic raytracing in anisotropic media: theory and application. //75th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts.
http://scitation.aip.org/getabs/servlet/GetabsServlet?prog=normal&id=SEGEAB000024000001000096000001&idtype=cvips&gifs=Yes |
|
曹雷, 张金海, 姚振兴.
2015. 利用三维高斯射线束成像进行地震定位. 地球物理学报, 58(2): 481–494.
DOI:10.6038/cjg20150212 |
|
邓飞, 刘超颖.
2009. 三维射线快速追踪及高斯射线束正演. 石油地球物理勘探, 44(2): 158–165.
|
|
梁锴. 2009. TI介质地震波传播特征与正演方法研究. 北京: 中国石油大学.
http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=Y1543669 |
|
马德堂, 朱光明.
2006. 横向各向同性介质中的初至波旅行时计算. 石油地球物理勘探, 41(1): 26–31.
|
|
孙成禹, 张文颖, 倪长宽, 等.
2011. 能量约束下的高斯射线束法地震波场正演. 石油地球物理勘探, 46(6): 856–861.
|
|
徐涛, 徐果明, 高尔根, 等.
2004. 三维复杂介质的块状建模和试射射线追踪. 地球物理学报, 47(6): 1118–1126.
|
|
姚陈, 蔡明刚.
2009. 任意空间取向TI弹性张量解析表述. 地球物理学报, 52(9): 2345–2348.
DOI:10.3969/j.issn.0001-5733.2009.09.019 |
|
张东, 傅相如, 杨艳, 等.
2009. 基于LTI和网格界面剖分的三维地震射线追踪算法. 地球物理学报, 52(9): 2370–2376.
DOI:10.3969/j.issn.0001-5733.2009.09.023 |
|