地壳介质普遍存在着不同程度的地震各向异性(Thomsen,1986).地震波在各向异性介质中传播,其相速度随传播方向不同而变化;且在一般情况下,波的偏振矢量既不平行也不垂直于波的传播方向.在各向异性介质中,除了少数简单介质(如横向各向同性VTI介质)外,Christoffel方程的特征根和特征矢量的表达式一般都不是弹性参数的解析函数.由于地壳岩石大多是弱各向异性介质(Thomsen, 1986),所以可以运用微扰理论对Christoffel方程进行求解,从而获得特征根和特征矢量的近似解析表达式(Jech and Pšenčík,1989;Pšenčík and Gajewski, 1998;Zheng and Pšenčík,2002;Zheng,2004).特征根和特征矢量的线性化不仅使极为复杂的qP波,qS1和qS2波相互解耦,也为各向异性参数反演奠定了广泛的应用基础.
研究地壳介质的各向异性对于油气田勘探和地震震源区应力场的监测都具有重要意义(Crampin, 1985;Crampin et al., 1999).然而获得地壳介质的各向异性并不简单,因为要获得介质各向异性需要全方位的观测.利用在地表布设的各种地震观测系统(包括二维剖面和三维台阵)对地壳介质的各向异性进行研究具有一定的局限性,利用这种地表观测系统只能给出地壳介质部分各向异性参数的粗略估计.对于浅层地壳(深度通常为几千米),walkaway VSP观测可以为地壳各向异性研究提供最有效的方法.walkaway VSP不仅能够准确地确定介质的P波和S波速度,而且还能确定介质的各向异性参数.利用慢度矢量的3个分量可以完全确定各向异性介质的弹性参数,在这种情况下,偏振矢量可作为辅助信息参与计算.然而,在许多情况下,慢度矢量的所有3个分量未知,为了确定各向异性介质的弹性参数,偏振矢量就成了必不可少的信息(Zheng and Pšenčík,2002;Zheng,2004).当只有一个钻井存在的情况下,只能获得慢度矢量沿钻井方向的一个分量.Gaiser(1990)提出钻孔中慢度矢量的垂直分量由在钻孔中记录到的走时确定,假设介质是横向均匀的,则可以利用走时的互易性确定钻孔中慢度矢量的水平分量.为了得到Java Sea的一条walkaway VSP各向异性参数,Horne和Leaney(2000)通过计算两千多个VTI模型并与观测资料进行拟合得到了一个误差最小的VTI模型.与Gaiser和Leaney使用的方法不同,Gomes等(2004)、刘宪彬和郑需要(2013)使用弱各向异性介质中的线性化反演方法对同一条Java Sea walkaway VSP数据进行了反演,得到了更为精确的各向异性参数.田鑫等(2015)使用qP波各向异性坐标变换方法和最小二乘求解方法,得到了与一般各向异性介质最接近的正交各向异性及其对称轴方向参数的计算公式.使用这些方法,他们对2条互相正交的walkaway VSP剖面数据进行了各向异性反演,获得了钻井中多个检波器处介质的WA参数.
因为VTI介质中传播的地震波的相速度和偏振矢量具有解析表达式,使得VTI介质在地震学和勘探地震学中被广泛研究.然而即使是这种简单的各向异性介质,仅仅使用地表观测系统也很难获得其准确的参数.为此,我们首先基于微扰理论和qP波坐标变换(Zheng,2007)推导了仅仅使用qP波反演TTI介质各向异性参数的公式.然后设计一个由多条剖面组成的walkaway VSP观测系统,分别对由2条相互正交的剖面,3条和6条等方位的剖面组成的3个数据集进行数值模拟和反演计算.最后对一个实际的walkaway VSP观测系统数据进行反演计算.
1 弱各向异性正反演公式如果介质的各向异性强度比较弱,可以利用微扰理论对Christoffel方程进行线性化近似求解,获得地震波在弱各向异性介质中传播时的解析表达式(Zheng and Pšenčík, 2002; Zheng,2004).对于qP波, 相速度和偏振矢量可以分别表达为
(1) |
(2) |
其中,方程(2)中使用了爱因斯坦求约定.在方程(1)和(2)中,v3和gj(3)分别代表qP波的相速度和偏振矢量,pi是各向异性介质中的慢度矢量,pkpk=v3-2.矢量ei(1), ei(2)和ni=ei(3)是定义在接收点的各向同性介质中沿P波射线的3个相互垂直的单位矢量.参考各向同性介质中的P波和S波的速度分别为α和β.矢量ni是射线的切线方向,单位矢量ei(1)和ei(2)可以在垂直于ni的平面内任意选取,但是下列选择将对公式的推导和计算带来很大的方便(Pšenčík and Gajewski, 1998; Zheng and Pšenčík, 2002):
(3) |
其中,
在波的传播方向上,n=(cosϕnsinθn, sinϕnsinθn, cosθn),其中,ϕn是方位角,而θn是波传播方向与z轴的夹角(0≤ϕn≤2π, 0≤θn≤π).所以,D=sinθn,而且e(1), e(2)可写为
(4) |
为了避免D=0或者非常接近于0的时候(3)式中的分母为零的问题,我们选择使用(4)式.
方程(1)和(2)中的B13, B23和B33被称为弱各向异性矩阵元素,它们是15个WA参数的函数(见附录),这15个WA参数可以完全描述弱各向异性介质中qP波的性质(Pšenčík and Gajewski, 1998; Zheng and Pšenčík, 2002).重写方程(1)和(2),使包含WA参数弱各向异性矩阵元素B13, B23和B33出现在方程的左边,它们当中包含着待求解线性方程组的未知参数:
(5) |
(6) |
方程(5)和(6)代表由15个WA参数作为未知参数的线性方程组.方程组等号右端的偏振矢量gk和慢度矢量pk可以从观测中得到.
对于VTI介质,qP波有3个独立的WA参数,它们分别是εx, εz和δx,这时B13, B23和B33可以分别表示为
(7) |
如果VTI介质的对称轴指向空间任意方向,那么介质就变成了TTI介质.描述TTI介质有5个独立参数,除了3个WA参数外,还有2个描述对称轴方向的角度参数,用ϕ和θ表示.它们分别代表TTI介质对称轴的方位角和与垂直轴的夹角.利用Zheng(2007)得到的关于qP波各向异性参数的坐标变换公式对B13和B33进行变换,可以得到
(8) |
(9) |
其中,系数h1, i, h3, i, ci, 1, ci, 2, ci, 3(i=1, …, 15)以及WA参数的系数在附录中给出.联立公式(5)、(6)和(8)、(9),得到TTI各向异性介质的反演公式.公式左边的各向异性矩阵元素B13和B33里包含3个WA参数εx,εz,δx和2个描述对称轴的方向参数ϕ和θ,这5个参数是待反演的参数,公式的右边是观测量.如果仅仅知道偏振矢量的3个分量和慢度矢量的垂直分量,可以通过消去qP波慢度矢量的2个水平分量获得下列反演公式:
(10) |
式中η代表P波慢度矢量的垂直分量,Δη代表qP波慢度垂直分量的扰动(Zheng,2004).
2 数值模拟实验考虑一个由6条剖面组成的walkaway VSP实验,其相邻剖面方位间隔为30°.沿每条剖面布设32个炮点,位于井孔的两边各有16个炮点,距离井孔最近炮点为0.1 km,最远炮点距离为3.1 km,相邻炮间距为0.2 km.井中等间距布设12个检波器,第1个检波器距地表 1.00 km,第12个检波器深度为1.55 km,相邻检波器间距为0.05 km(图 1),图中字母L后面的数字为剖面的方位角, 单位为°.
TTI介质模型是一个弹性参数随深度线性变化的梯度模型,给定地表的TTI(密度归一化的)弹性参数矩阵和5 km深处的弹性参数矩阵,这2个矩阵是分别由相应的2个VTI弹性矩阵经过2次坐标系旋转得到.
地表和地下5 km处的VTI弹性(对称)矩阵分别为
首先对这2个弹性矩阵绕坐标轴y顺时针旋转40°,然后再围绕z轴顺时针旋转30°.2次旋转后得到的地表和地下5 km处的弹性矩阵分别为
经过2次坐标变换之后,矩阵的所有元素都不为零.我们可以把这种介质看作是一般各向异性介质,其中qP波的各向异性大小为4%.
在这一观测系统中,每一个接收点有192条射线.使用各向异性介质中的运动学和动力学射线追踪Anray软件(Gajewski and Pšenčík, 1990),可以计算出震源到接收点之间的射线和检波器记录到的理论地震图.在模型的上部,射线覆盖相当好,最浅的接收点不仅被近于水平的射线覆盖,也被近垂直的射线覆盖.最深的接收点主要是被近垂直的射线覆盖,没有近水平的射线照射在这些接收点上.因此最深的接收点的射线覆盖不如浅部的接收点覆盖的那么好,我们并不期望在深的接收点处得到非常理想的反演结果.
从“观测的”合成地震图上可以获得qP波的初至和偏振矢量.在每一个接收点,我们使用近垂直的射线走时决定P波的参考速度α,S波的速度为β=
对于一般各向异性介质,完全恢复qP波的15个WA参数需要5条剖面(Pšenčík and Gajewski, 1998).对于具有高度对称性的TTI介质,2条互相正交的剖面可以决定3个qP波的WA参数和2个方向参数(田鑫等,2015).图 2给出了使用2条互相正交剖面(L0和L90)反演得到的检波器1, 3, 5, 7, 9, 11处qP波相速度分布曲线,偶数检波器处的相速度分布曲线与奇数检波器的相速度分布曲线类似.如前所述,由于正演模型和反演公式都使用TTI介质,所以line 30,line 60,line 120和line 150所在切面的qP波的相速度与line 0和line 90所在切面的相速度具有同样的精度.表 1给出了使用2条剖面反演得到的每一个检波器处qP波相速度的最大相对误差,最大误差出现在检波器12处,其值为0.81%.大部分检波器处TTI对称轴方向误差小于或等于4°,除了2个最深处的检波器TTI对称轴方向误差大于4°.对于每一个检波器,我们给出了6个方位的垂直于地表的切面内qP波相速度随偏振角(偏振矢量与z轴的夹角)的变化曲线.水平轴代表偏振角,沿剖面正方向偏振角为正,沿剖面反方向偏振角为负.由于对称性,偏振角的取值范围为-90°至90°.纵轴代表qP波的相速度,单位为km·s-1.在所有切面图中实线代表qP波相速度的精确值(由正演计算得到),虚线代表反演结果.对于每一个检波器的6个方位切面图,位于左边的3个切面图中的一个切面(line 0)与参与反演的L0剖面重合,右边的3个切面图中的一个切面(line 90)与参与反演的L90剖面重合.由于正演模型和反演公式都使用TTI介质,所以line 30, line 60, line 120, line 150四个切面图给出了与line 0和line 90两个切面图同样精度的结果.从图中可以看出,反演结果与精确解存在微小的差别,反演误差主要来源是因为处于较深的检波器没有很好的覆盖,求取慢度矢量分量时由于所采用的光滑函数存在缺陷,不能得到准确的慢度矢量分量,还有一个很重要的原因是舍弃了慢度矢量的切向分量,因为在仅有一口钻井存在的情况下,无法获得慢度矢量的切向分量信息.
由于实际野外观测系统有3条剖面(见第3节),所以我们对基于3条剖面(L0,L60,L120)的数据集进行了反演计算,得到了TTI介质的3个WA参数和对称轴的方向.其中12个检波器处qP波的相速度精度高于使用2条正交剖面得到的结果,最大相对误差出现在检波器7和8处,其值为0.47%.反演得到了TTI介质对称轴的两个方向参数,其中对称轴的方位角与精确解(ϕ=30°)完全相同,除了最深的检波器处的对称轴方向θ有3°的误差外,对称轴与z轴的夹角与精确解(θ=40°)有1°~2°的误差.
我们也得到了使用6条剖面数据反演TTI介质WA参数的结果.相速度随方向的变化图像与使用2条和3条剖面的结果类似,但相速度误差和TTI对称轴的误差进一步减小.反演得到的qP波相速度的最大相对误差为0.38%,对称轴方向角ϕ与理论值完全相同,θ的最大误差为3°.
3 实际观测数据反演实际观测地震资料来自Java Sea地区布设的由3条剖面组成的一个walkaway VSP观测系统.其中剖面L1方位角为25°,剖面L2方位角为84°,剖面L3方位角为146°.沿3条剖面分别布设了315个,310个和312个震源,它们分布在钻孔两侧-3.1~3.2 km的范围内,相邻震源间隔为20 m.在钻孔2.30~2.59 km深度范围内等距离布设30个三分量检波器,相邻检波器间距为10 m, 其中2.40~2.44 km为勘探目标层, 检波器序号依次为11, 12, 13, 14和15.数据由下行(直达波)和上行(反射)qP波组成.由于上行波含有较大的噪声而且不是初至波,我们仅研究下行波的反演问题.
walkaway VSP走时曲线可分为两类:一类是共炮点走时曲线,即地表炮点为震源,井中不同深度的检波器为接收点,通常横坐标表示深度,纵坐标表示地震波走时;另一类是共接收点(或检波器)走时曲线,根据地震波走时互易原理,可以假定接收点为震源,炮点接收地震波.通常取炮点到井口的连线方向为横坐标,走时为纵坐标.首先在地震图中拾取P波初至,将初至分别按共炮点和共接收点组成走时曲线.沿不同剖面获得多条共接收点走时曲线和共炮点走时曲线.接收点处慢度矢量的垂直分量可以通过计算共炮点走时曲线的斜率得到,但接收点处的水平分量并不能直接获得.因为地震波的走时具有互易性,如果介质是横向均匀的,则炮点到接收点的走时等于接收点到炮点的走时.在这种情况下,可以把接收点当作炮点,利用共接收点的走时计算慢度矢量的水平分量.当介质横向非均匀性比较弱时,所得结果可作为一个近似;如果介质的横向非均匀性很强,就不能得到接收点处的水平分量.若想得到水平分量,必须在邻近的井中放入检波器进行观测.共接收点和共炮点走时曲线一般情况下都不光滑,它们的一阶导数均不连续甚至不存在.为了得到可靠的慢度矢量分量,必须对走时曲线进行光滑处理.使用最小二乘三次样条函数对走时数据进行拟合,三次样条函数的一阶导数即为接收点处慢度矢量的垂直分量和水平分量.
从检波器记录到的来自所有震源的三分量地震数据中获得了偏振矢量和慢度矢量的垂直分量和水平分量.在对这些数据进行分析后,剔除了相对误差大于30%的数据.每个接收点处的参考速度是综合考虑测井资料和观测到的慢度矢量得到的.
本文使用TTI公式和一般各向异性公式对所得到的数据进行了反演.图 3给出了勘探目标层中3个检波器观测到的慢度矢量的垂直分量(×)以及使用TTI模型反演得到的慢度矢量的垂直分量(实线).在每一组图中, 纵坐标为慢度矢量的垂直分量,横坐标为偏振角.偏振角大于零代表剖面的正方向,偏振角小于零代表剖面的负方向.偏振角的绝对值小于90°的慢度矢量的垂直分量对应于下行波(直达波),偏振角的绝对值大于90°的慢度矢量的垂直分量对应于上行波(反射波).容易看出,上行波很离散并包含有较大的噪声.
图 4给出了反演得到的勘探目标层中3个检波器所在深度处的相速度分布.横轴是P波传播的方向与z轴正方向的夹角,纵轴为相速度大小.横轴正方向与剖面方位方向一致,横轴负方向指向离开井口相反的方向.图中虚线代表使用TTI模型反演得到的结果,实线代表由一般各向异性模型反演得到的结果.如果介质是各向同性介质,反演得到的每个接收点不同方位剖面下的介质的相速度分布应该是一条直线,不随波传播的方向而变化.显然,从图 4可以看出,所有不同深度接收点处的介质都呈现不同程度的各向异性.地震波的相速度不仅与剖面的方位有关,而且与波传播方向与垂直轴的夹角有关.表 2给出了使用3条剖面和TTI模型反演得到的目标层中每一个检波器处qP波的各向异性参数和对称轴的方向.不同检波器处的3个各向异性参数变化不大,对称轴与z轴的最大偏离为10°,对称轴的平均方位角为54.4°.表 3给出了使用3条剖面和TTI模型反演得到的目标层中每一个检波器处qP波的相速度最小值、最大值和各向异性大小.可以看出,不同检波器处和不同剖面的各向异性强度(an)变化并不显著,平均值为9.27%.不同剖面相速度的一致性主要是使用TTI模型所导致.从反演结果中可以看出一般各向异性模型反演结果与TTI模型反演结果存在较大的差异(图 4).尽管两种模型反演得到的对称轴方向比较接近,但是不同剖面不同接收点处的相速度和各向异性值有较大的不同.也可以发现两种模型反演得到的相速度最小值、最大值和各向异性值存在着差别.沿剖面L1不同接收点处的WA参数变化较大,最小值接近于TTI模型的各向异性值,最大值是TTI模型的2倍.沿剖面L2不同接收点处的各向异性最小值为13.8%,最大值为17.4%.沿剖面L3不同接收点处的各向异性值接近于TTI模型的各向异性值.反演结果还给出了用TTI模型和一般各向异性模型反演结果的误差(表 3最后一列),这个误差是反演方程等号左边的物理量(计算结果)和右边物理量(观测结果)之差的平方和.对于每一个检波器,使用TTI模型进行反演的误差明显大于使用一般各向异性模型进行反演的误差.
本文在地壳介质是弱各向异性介质的假设下,从Zheng和Pšenčík(2002)的一般弱各向异性反演公式, Zheng(2004, 2007)关于qP波各向异性参数坐标变换公式中得到了适用于TTI介质的WA参数反演公式.本文理论和数值实验均证明,在walkaway VSP观测中,2条正交的剖面完全可以决定TTI介质的3个独立的关于qP波的WA参数和对称轴方向.如果利用更多的观测剖面(≥3)进行反演,可以提高WA参数和对称轴方向的精度.由于大多数地壳介质并非简单的TTI介质,使用任意各向异性介质模型对单条或多条剖面进行反演更有实际意义.
我们知道,研究介质各向异性有实验室岩石样品测量方法和野外观测和勘探等方法.实验室测量通常使用高频超声波,能够得到非常复杂的介质的各向异性参数(Arts et al., 1991; Zheng and Tomaskova, 2004).然而,由于岩石样品尺寸的限制和使用高频信号,所获得的介质各向异性参数与野外大尺度和低频观测得到的各向异性参数有一定的差别.野外观测和勘探研究介质各向异性有许多方法(滕吉文等, 1992;张中杰, 2002;高原和滕吉文, 2005),例如体波观测方法、面波观测方法、二维和三维地震勘探等等.由于大部分观测台站位于地表,通过这些方法很难得到对地下介质各向异性全面和准确的描述.
在油气勘探中,如果没有考虑地下介质存在各向异性,将会导致计算得到的界面深度存在较大的误差.如果进行偏移成像,所得地下反射界面的图像也可能是比较模糊的.随着地震勘探的快速发展和勘探区域地质构造复杂性的增大,获得精确的地下介质各向异性参数对于提高勘探效率和降低成本已成为当务之急,walkaway VSP是一种获得地下介质各向异性参数的可靠和精确的方法,所以在油气田勘探,特别是页岩气和页岩油勘探中具有广阔的应用前景.
附录
(1) qP波弱各向异性矩阵元素
(A1) |
(A2) |
(A3) |
(2) qP波弱各向异性参数
(A4) |
其中Aij是密度归一化的弹性常数.
(3) 公式(8)和(9)中WA参数的系数
(A5) |
(A6) |
(A7) |
Arts R J, Helbig K, Rasolofosaon P N J.
1991. General anisotropic elastic tensor in rocks:Approximation, invariants, and particular directions. //61st Annual International Meeting, SEG. Expanded Abstracts: 1534-1537.
|
|
Crampin S.
1985. Evaluation of anisotropy by shear-wave splitting. Geophysics, 50(1): 142-152.
DOI:10.1190/1.1441824 |
|
Crampin S, Volti T, Stefánsson R.
1999. A successfully stress-forecast earthquake. Geophysical Journal International, 138(1): F1-F5.
DOI:10.1046/j.1365-246x.1999.00891.x |
|
Gaiser J E.
1990. Transversely isotropic phase velocity analysis from slowness estimates. Journal of Geophysical Research:Solid Earth, 95(B7): 11241-11254.
DOI:10.1029/JB095iB07p11241 |
|
Gajewski D, Pšenčík I.
1990. Vertical seismic profile synthetics by dynamic ray tracing in laterally varying layered anisotropic structures. Journal of Geophysical Research:Solid Earth, 95(B7): 11301-11315.
DOI:10.1029/JB095iB07p11301 |
|
Gao Y, Teng J W.
2005. Studies on seismic anisotropy in the crust and mantle on Chinese mainland. Progress in Geophysics (in Chinese), 20(1): 180-185.
|
|
Gomes E, Zheng X, Pšenčík I, et al.
2004. Local determination of weak anisotropy parameters from walkaway VSP qP-wave data in the Java sea region. Studia Geophysica et Geodaetica, 48(1): 215-231.
DOI:10.1023/B:SGEG.0000015593.29477.88 |
|
Horne S A, Leaney W S.
2000. Polarization and slowness component inversion for TI anisotropy. Geophysical Prospecting, 48(4): 779-788.
DOI:10.1046/j.1365-2478.2000.00205.x |
|
Jech J, Pšenčík I.
1989. First-order perturbation method for anisotropic media. Geophysical Journal International, 99(2): 369-376.
DOI:10.1111/gji.1989.99.issue-2 |
|
Liu X B, Zheng X Y.
2013. Calculation of weakly anisotropic parameters by slowness and polarization vectors of qP wave. Acta Seismologica Sinica (in Chinese), 35(2): 184-198.
|
|
Pšenčík I, Gajewski D.
1998. Polarization, phase velocity and NMO velocity of qP waves in arbitrary weakly anisotropic media. Geophysics, 63(5): 1754-1766.
DOI:10.1190/1.1444470 |
|
Teng J W, Zhang Z J, Wang A W, et al.
1992. The study of anisotropy in elastic medium:evolution, present situation and questions. Progress in Geophysics (in Chinese), 7(4): 14-28.
|
|
Thomsen L.
1986. Weak elastic anisotropy. Geophysics, 51(10): 1954-1966.
DOI:10.1190/1.1442051 |
|
Tian X, Hong Q Y, Zheng X Y.
2015. Seismic anisotropic inversion of multiple walkaway VSPs. Acta Seismologica Sinica (in Chinese), 37(2): 266-277.
|
|
Zhang Z J.
2002. A review of the seismic anisotropy and its applications. Progress in Geophysics (in Chinese), 17(2): 281-293.
|
|
Zheng X Y.
2004. Inversion for elastic parameters in weakly anisotropic media. Geophysical Journal International, 159(3): 1077-1089.
DOI:10.1111/gji.2004.159.issue-3 |
|
Zheng X Y.
2007. Identification of the symmetry planes in weakly anisotropic elastic media. Geophysical Journal International, 170(1): 468-478.
DOI:10.1111/gji.2007.170.issue-1 |
|
Zheng X, Pšenčík I.
2002. Local determination of weak anisotropy parameters from qP-wave slowness and particle motion measurements. Pure and Applied Geophysics, 159(7-8): 1881-1905.
DOI:10.1007/s00024-002-8713-z |
|
Zheng X Y, Tomaskova A.
2004. Inversion of travel times in weakly anisotropic rock samples. Pure and Applied Geophysics, 161(8): 1709-1724.
DOI:10.1007/s00024-004-2527-0 |
|
高原, 滕吉文.
2005. 中国大陆地壳与上地幔地震各向异性研究. 地球物理学进展, 20(1): 180–185.
|
|
刘宪彬, 郑需要.
2013. 利用qP波慢度和偏振矢量计算弱各向异性介质参数. 地震学报, 35(2): 184–198.
|
|
滕吉文, 张中杰, 王爱武, 等.
1992. 弹性介质各向异性研究沿革、现状与问题. 地球物理学进展, 7(4): 14–28.
|
|
田鑫, 洪启宇, 郑需要.
2015. 多测线变偏移距VSP地震各向异性反演. 地震学报, 37(2): 266–277.
|
|
张中杰.
2002. 地震各向异性研究进展. 地球物理学进展, 17(2): 281–293.
|
|