多测线变偏移距VSP地震各向异性反演

田 鑫 洪启宇 郑需要

（中国北京100081中国地震局地球物理研究所）

引言

地壳介质普遍存在着不同程度的地震各向异性（Thomsen，1986）.这种各向异性可能是由于岩石裂隙和裂缝，薄层的周期排列，应力引起的岩石颗粒、空隙和晶体的定向排列等原因所引起.无论地震波在各向同性还是在各向异性介质中传播均遵守Christoffel方程（ˇCervený，2001）.地震波在各向同性介质中传播有P波和S波，两种波均为线性偏振波，P波的传播方向与偏振方向相同，S波的传播方向与偏振方向垂直，且两种波的相速度与波的传播方向无关.与波在各向同性介质中传播不同，在各向异性介质中，波的相速度和偏振矢量是传播方向的函数，一般情况下波的偏振矢量既不平行也不垂直于波的传播方向.同时，S波在各向异性介质中传播时还会发生分裂现象.在各向异性介质中，除了少数简单介质（如横向各向同性VTI介质）外，一般Christoffel方程的特征根和特征矢量的表达式都不是弹性参数的解析函数.由于多种地壳岩石都是弱各向异性介质（Thomsen，1986），所以可以运用微扰理论对Christoffel方程进行求解，以获得特征根和特征矢量的近似解析表达式（Jech，Pšenˇcík，1989；Zheng，2004）.其中，qP波的表达式是线性的，qS波的是非线性的.Zheng和Pšenˇcík（2002）同时考虑介质相对于参考介质的扰动和波传播方向的扰动，得到了更为普遍的关于qP波的解析表达式.基于微扰理论，Zheng（2004）得到了关于qP波和qS波的反演方程，并且给出了已知qP波和qS波慢度矢量的1，2或3个分量以及偏振矢量情况下反演介质WA参数的公式.

地震各向异性理论在地震学和地球物理学中有着广泛的应用.在地震勘探中，利用多条变偏VSP能够准确地确定介质的WA参数和地震波传播的相速度（郑需要等，2010；刘宪彬，郑需要，2013）.如果使用慢度矢量的全部3个分量，可以很容易地确定介质的各向异性参数，而偏振矢量可以作为辅助信息和约束条件参与WA参数的反演计算.然而在许多情况下，由于多种原因不能得到慢度矢量的所有3个分量.如果仍想确定介质的各向异性性质，就必须知道偏振矢量.在野外观测中，如果仅有一个钻井存在，严格地讲，只能获得慢度矢量的垂向分量和波的偏振矢量.在这种情况下，我们仍然希望获得介质的各向异性性质.本文将详细介绍如何利用一条变偏VSP剖面和两个相互正交的变偏VSP剖面计算介质各向异性参数的方法，并将此方法应用到两条相互正交变偏VSP测线实际观测数据的反演和解释中，同时对所获得的各向异性参数和相速度分布进行较为详尽的讨论.

1 各向异性正演公式

Pšenˇcík和 Gajewski（1998）以及 Zheng和 Pšenˇcík（2002）曾对地震波在各向异性介质中的传播理论进行详细描述.为方便阅读起见，本节和下面几节将给出常用的基本公式.

地震波在各向同性和各向异性介质中传播均满足Christoffel方程 （ˇCervený，2001），即

式中：¯Γjk＝aijklninl表示Christoffel矩阵元素，ni为波传播方向的单位矢量，aijkl为各向同性或各向异性介质中密度归一化的弹性参数张量；¯Gm＝v2m（m＝1，2，3）表示Christoffel矩阵的3个特征根；g（m）j为对应的特征矢量；vm为波传播的相速度，m＝1，2和3分别代表各向异性介质中3种不同的波（qS1波、qS2波和qP波）.如果介质的各向异性比较弱，就可以利用扰动理论对Christoffel方程进行线性化近似求解，以获得地震波在弱各向异性介质中传播时的解析表达式.

对于qP波可以得到相速度和偏振矢量的表达式为

式中：pi为慢度矢量；被称为弱各向异性矩阵元素（见附录），＝1，2，3）是沿射线坐标系中相互正交的3个单位矢量（Pšenˇcík，Gajewski，1998），其与波传播方向n的关系为；α和β分别表示地震波在各向同性参考介质中P波和S波的速度.

式（2）是Christoffel方程的一个特征根，偏振矢量公式（3）是与特征根对应的特征矢量.如果我们知道WA参数和波传播的方向，就可以使用公式（2）和（3）计算qP波在该方向传播的相速度和偏振矢量.这两个公式通常被称为慢度公式和偏振公式.公式（2）和（3）的右边都包含有与弹性参数相联系的WA参数，左边是相速度和偏振矢量（Zheng，2004）.

2 各向异性反演公式

地震波（qP波或两个qS波）在弱各向异性介质中传播，其慢度矢量可以表示为

式（4）中，如果忽略Δξ，Δζ和Δη的二次项，可以得到慢度平方的近似表达式为

比较式（2）与（5），可以得到

把式（3）的两边分别乘以e（K）i（K＝1，2）可以得到

使用式（6）和（7）可以对qP波进行线性化反演.公式（6）和（7）的左边是15个待求的WA参数，右边是野外观测到的物理量，表示qP波的偏振矢量和慢度矢量.这15个WA参数能够很好地确定qP波在弱各向异性介质中的传播过程.假如我们能够在观测中知道慢度矢量的3个分量，那么仅利用式（6）就可以完全确定一般弱各向异性介质与qP波相联系的15个WA参数，这时式（7）可以作为辅助或约束条件参加反演（Zheng，2004）.

多种原因使得野外观测受到限制.例如，由于经费原因，不可能为了得到慢度矢量的水平分量而同时钻两口相距很近的井.在许多实际地震观测中，我们经常不能得到慢度矢量的所有3个分量（Horne，Leaney，2000）.如果根据观测仅仅得到慢度矢量的一个或两个分量，那么就很难确定所有独立的WA参数.在这种情况下，利用偏振矢量信息就能使反演顺利进行.在公式（6）和（7）中，可以消去qP波慢度矢量的两个水平分量，得到一个与慢度矢量水平分量无关的方程：

式中，η（P）为参考介质中P波慢度矢量的垂直分量，Δη（qP）为观测到的qP波慢度矢量的垂直分量的扰动值.在偏振矢量和慢度矢量的一个分量已知的情况下，式（8）是反演介质WA参数的一般公式.

反演公式（6）、（7）和（8）中有15个相互独立的WA参数.若要唯一确定这15个参数，需布设5条变偏VSP剖面（Pšenˇcík，Gajewski，1998）.如果仅使用一条变偏 VSP剖面，则可得到5个相互独立的WA参数.假定剖面和井所构成的平面与（x，z）平面相同，qP波在这个平面内传播时，慢度矢量分量在n2方向（y轴）的投影等于零.这时式（8）中有5个独立的 WA参数，它们是εx，εz，δx，ε15和ε35，这5个参数完全决定了在（x，z）平面内传播的qP波的性质.利用qP波可以反演得到5个独立的WA参数.在两个相互正交的变偏VSP剖面下，取第一个位于（x，z）面，第二个位于（y，z）面.如果分别对两个剖面进行反演，对于每一个剖面都可以得到5个WA参数，其中有一个共同的参数是εz.如果对两个剖面进行联合反演，则可得到9个相互独立的WA参数.除了上面的5个WA参数外，剩余的4个参数为εy，δy，ε14和ε24.这9个 WA参数完全决定了（x，z）和（y，z）两个平面内介质的各向异性性质.εx，εy和εz分别代表地震波沿3个坐标轴方向传播时的各向异性大小.当地震波传播方向在xoz和yoz平面内偏离z轴一个较小角度时，δx和δy近似代表地震波沿该方向的各向异性（Thomsen，1986）.这两个参数在各向异性反射地震观测研究中有着重要应用，其它4个WA参数则与介质的对称性有关.

3 从一般弱各向异性介质中计算近似的TTI各向异性参数

VTI介质是一种最简单并且研究最广泛的各向异性介质.不仅是因为这种介质简单，更重要的是因为地震波在VTI介质中传播，其相速度和偏振矢量均有解析表达式.VTI介质有5个独立的参数，当其对称轴沿水平方向时，称为HTI介质；如果其对称轴指向空间任意方向，则称为TTI介质.尽管TTI介质可能有21个弹性参数，但独立的参数仍然是5个.虽然VTI介质的相速度和偏振矢量有解析表达式，但是非常复杂，它们是弹性参数的非线性函数.Thomsen（1986）在假定介质为弱各向异性的情况下，对VTI介质进行线性化解释，得到了用5个简单参数描述该介质的相速度，其中ε，δ和γ称为Thomsen参数.

在介质为弱各向异性的条件下，一般用21个弱各向异性参数进行描述（Pšenˇcík，Gajewski，1998；Zheng，2004）.对于qP波，需要15个WA参数进行描述（见附录）.如前所述，要完全决定这15个 WA参数，至少需要5条测线（Pšenˇcík，Gajewski，1998；Zheng，2004）.在两条测线的情况下，可以确定9个WA参数.为了更好地解释这9个参数，我们寻找最接近这种介质的TTI介质所对应的WA参数.对于qP波，TTI介质可以用两个代表方向的参数和3个 WA参数（或者Thomsen参数）进行描述，它们是φ，θ，εz，εx和δx.其中φ为对称轴的方位角，θ为对称轴与z轴正方向的夹角.

假定在坐标系（x，y，z）里VTI介质用5个WA参数表示（其中3个参数描述qP波），该介质的对称轴沿z方向.为方便起见，我们把该坐标系称为自然坐标系.首先将自然坐标系（x，y，z）绕y轴旋转角度θ，得到新的坐标系（x′，y′，z′）.在新的坐标系中，尽管VTI介质中的qP波可以用9个WA参数描述，但它们并不独立.它们与自然坐标系里的3个独立的WA参数之间的关系可以用附录中式（A5）表示（Zheng，2007）.

如果把刚得到的新的坐标系（x′，y′，z′）中的9个WA参数再围绕坐标轴z′旋转角度－φ，在新的坐标系（x″，y″，z″）中得到对称轴为任意指向的TTI介质，其WA参数由附录中式（A6）给出.在（x″，y″，z″）坐标系里，如果有两个相互正交的剖面，其中一个沿x″轴，另一个沿y″轴.如前所述，沿这两个坐标轴，对于qP波有9个独立的WA参数，即

如果将附录中式（A8）代入式（9），可以得到任意坐标系中TTI介质在两个相互正交的坐标平面（x，z）和（y，z）内的9个WA参数和自然坐标系中的3个WA参数的关系.

从附录中式（A6）可以得到沿x″轴和y″轴的TTI介质各向异性矩阵元素分别为

如果将式（9）和附录中式（A5）代入式（10）和式（11），可以得到弱各向异性矩阵元素为TTI介质中3个WA参数和对称轴的两个方向参数的函数.

利用式（9）、（10）、（11）和附录中式（A5），可以构造最小二乘公式为

为了得到自然坐标系（x，y，z）中的VTI介质的3个WA参数，需要给出VTI介质的对称轴方向，然后利用附录中式（A8）和式（9）求解式（12）的最小二乘解.通过寻找相速度最小误差的最小二乘解，可以得到最接近一般各向异性介质的VTI介质的对称轴方向和3个WA参数.同理，利用该方法也可以得到与一般弱各向异性介质最接近的正交各向异性介质参数.对于qP波，描述其正交各向异性介质的参数为对称轴的3个方向参数和6个WA参数，其数量正好与由两条相互正交的测线得到的9个WA参数数量相同.

4 慢度矢量分量与偏振矢量的计算

4.1 慢度矢量分量的计算

本文研究的变偏VSP地震资料取自瓜哇海地区布设的两条相互正交的变偏VSP测线.其中测线1方位角为40°，测线2方位角为130°.两条测线上分别布设了320个和324个震源，分布在钻孔两侧－3.24—3.24km的范围内，相邻震源间隔约为20m.在钻孔3.0—3.1km深度范围内等距离布设10个三分量检波器，相邻检波器间距为10m.数据由下行（直达波）和上行（反射）qP波和qS波组成.本文仅研究qP波的反演问题.

变偏VSP走时曲线可分为两类：一类是共炮点走时曲线，即地表炮点为震源，井中不同深度的检波器为接收点，通常横坐标表示深度（km），纵坐标表示地震波走时；另一类是共接收点（或检波器）走时曲线，根据地震波走时互易原理，可以假定接收点为震源，炮点接收地震波.通常取炮点到井口的连线方向为横坐标，走时为纵坐标.首先在地震图中拾取P波初至，将初至分别按共炮点和共接收点组成走时曲线.沿测线1获得了10条共接收点走时曲线，320条共炮点走时曲线；沿测线2获得了10条共接收点走时曲线，324条共炮点走时曲线.接收点处慢度矢量的垂直分量可以通过计算共炮点走时曲线的斜率得到，但接收点处的水平分量并不能直接获得.因为地震波的走时具有互易性，如果介质是横向均匀的，则炮点到接收点的走时等于接收点到炮点的走时.在这种情况下，可以把接收点当作炮点，利用共接收点的走时计算慢度矢量的水平分量.当介质横向非均匀性比较弱时，所得结果可作为一个近似；如果介质的横向非均匀性很强，就不能得到接收点处的水平分量.若想得到水平分量，必须在邻近的井中放入检波器进行观测.共接收点和共炮点走时曲线一般情况下都不光滑，它们的一阶导数均不连续甚至不存在.为了得到可靠的慢度矢量分量，必须对走时曲线进行光滑处理.本文使用最小二乘三次样条函数对走时数据进行了拟合，三次样条函数的一阶导数即为接收点处慢度矢量的垂直分量和水平分量.

4.2 偏振矢量的计算

检波器记录到的P波的偏振矢量是由三分量地震图的质点运动轨迹合成得到的.由于多种原因，P波的偏振矢量常常偏离线性偏振.为了得到P波的线性偏振矢量，需要使用空间三维最小二乘法拟合P波质点的运动轨迹.

已知空间直线的标准方程为

式中，xc，yc和zc是直线上已知点的坐标，对该方程进行整理得到直线的射影式方程为x＝X（z－zc）／Z＋xc，y＝Y（z－zc）／Z＋yc，即

式中，a＝X／Z，b＝Y／Z.式中的两个方程均为三维平面方程，如果这两个平面的法线不一致，它们必然相交于一条直线，所以可以分别对两个方程进行数据拟合（袭杨，2009），即

当Q取最小值时，a值和b值即为方程的系数.它们满足

时，Q值最小.由式（15）和（16）可以得到

进而可以得到P波偏振矢量的方向余弦为

对井中检波器记录到的P波在获得初至和振幅最大值的基础上，截取时间长度等于初至到最大振幅之间时间长度3倍的地震记录数据，利用上面推导出的公式计算P波的偏振矢量.计算过程中使用了质心公式，避免出现奇异解的可能性.

5 实际资料反演

Zheng和Pšenˇcík（2002）以及Zheng（2004）曾利用数值模拟方法分别对3个和5个变偏VSP剖面进行正反演计算，并对反演结果进行了理论分析和讨论.这里我们使用式（6）—（8）对在瓜哇海地区获得的两条相互正交的变偏VSP测线资料进行反演计算.

我们从检波器记录到的来自所有震源的三分量地震数据中获得了偏振矢量和慢度矢量的垂直分量和水平分量.在对这些数据进行分析后，发现有些数据包含较大的误差，对反演结果有很大的影响.因此，在反演中我们没有使用相对误差大于30%的数据.每个接收点处的参考速度是综合考虑测井资料和观测到的慢度矢量得到的.由于反演公式除了假设介质是弱各向异性外，没有其它任何假定条件，因此其完全适用于任意复杂介质，特别是不受上覆介质和介质的横向非均匀性的影响.

在反演中，可以单独使用下行波，也可以单独使用上行波（反射波），还可以联合使用下行波和上行波.因为上行波中存在较大的误差，我们没有单独使用上行波进行反演.下面主要给出联合使用下行波和上行波进行反演的结果.

图1 10个接收点所在深度处的介质各向异性相速度分布曲线绿色和黑色曲线分别表示测线1和测线2的相速度分布，红色和蓝色曲线分别表示测线1和测线2与一般各向异性最接近的TTI介质的相速度分布Fig.1 Phase velocity distribution of qP wave at ten receivers with different depth for line 1（green）and line 2（black）obtained by inversion in general anisotropic medium.Red curve represents phase velocity for line 1and blue curve for line 2corresponding to TTI model.The depth of ten receivers is 3.00，3.01，3.02，…，3.09km in order

图1 给出了反演得到的10个接收点所在深度处的介质各向异性相速度分布图.横轴表示P波传播方向与z轴正方向的夹角，纵轴表示相速度分布值.横轴正方向与测线方位方向一致，横轴负方向为离开井口与测线的负方向一致.如果介质是各向同性，则反演得到的每个接收点不同方位测线下的介质的相速度分布应该是一条直线，不随波传播的方向而变化.显然，从图1中可以看出，所有不同深度接收点处的介质均呈现不同程度的各向异性，介质的相速度不仅与测线的方位有关，而且与波传播方向与垂直轴的夹角有关.对两条测线所表示的曲线进行比较发现，大部分接收点处介质相速度分布比较一致，但是第6个接收点处的两条测线得到的相速度分布存在较大差别.表1给出了使用一般各向异性介质反演得到的9个 WA参数，其中εx，εy和εz的极大值分别是0.047，0.055和0.082，依次出现在第3个接收点和第4个接收点.特别应该指出的是，WA参数的最大值出现在第9个接收点，它们是δx和δy，其值分别为0.311和0.248，仍然在合理的范围内（Thomsen，1986）.

表1 使用一般各向异性介质反演得到的WA参数Table 1 WA parameters obtained by inversion based on the general anisotropic medium

使用两个相互正交的剖面，对于qP波我们得到了9个独立的WA参数.尽管9个独立参数不能完全确定一般任意各向异性介质（任意各向异性介质有15个独立的参数），但是它们能够确定具有高度对称性的各向异性介质.对于正交各向异性介质，如果仅使用qP波，则需要6个独立的WA参数和3个方向参数.对于TTI介质，使用qP波需要3个独立的WA参数和两个方向参数.因此可以看出，用两个相互正交的剖面可以确定介质是否为正交各向异性介质或TTI介质.我们使用第3节中给出的方法计算了与一般各向异性介质最接近的TTI介质的各向异性参数.从图1中可以看出代表TTI介质的相速度分布与其对应的一般各向异性介质相速度分布的差别.有的接收点处相速度分布差别比较大，有的差别比较小.相速度分布差别较大的说明一般各向异性介质偏离TTI各向异性介质，差别较小的意味着一般各向异性介质更接近于TTI各向异性介质.在所有的接收点中，第7个接收点处两条测线对应的相速度最为接近.如果一般各向异性介质与TTI各向异性介质的相速度分布很接近，则可以认为该接收点处的介质为TTI介质.

表2给出了从一般各向异性介质参数中求得的与其最接近的TTI介质的WA参数（εz，εx和δx）和TTI介质的对称轴方向参数（φ和θ）.其中第2列是各个接收点处P波的参考速度，其与一般各向异性反演时使用的参考速度相同（表1中第2列）.可以看出，第4个和第9个接收点处介质具有很强的各向异性，εz的绝对值分别为0.421和0.327.从对称轴的两个方向参数可以看出，所有接收点处的介质既不是VTI（对称轴沿垂直方向）介质也不是HTI（对称轴沿水平方向）介质，因为其对称轴与z轴的夹角既不等于0°也不等于90°.

表2 TTI介质的各向异性参数和对称轴方向Table 2 WA parameters and symmetry direction parameters of TTI media

6 讨论与结论

在地壳介质为弱各向异性介质的假设下，我们从Zheng和Pšenˇcík（2002）的一般弱各向异性反演公式（Zheng，2004）中得到了适用于两条相互正交测线的 WA参数反演公式.在这种情况下，使用qP波可以确定9个独立的WA参数，这9个参数可以完全地描述与两条测线对应的接收点处介质的各向异性性质.为了寻找与一般各向异性介质最接近的具有高对称性的正交各向异性介质和TTI各向异性介质，我们使用最小二乘方法和Zheng（2007）提出的qP波各向异性坐标变换方法，得到了与一般各向异性介质最接近的正交各向异性和TTI各向异性参数及其对称轴方向参数的计算公式.通过引入质心计算方法，获得了求取P波偏振矢量的计算公式.该方法计算简单、稳定、高效，而且不存在解的奇异问题.数值计算和实际资料处理结果表明，本文所使用的反演方法能够很好地获得井中接收点处介质各向异性参数，是地震勘探中研究地壳介质各向异性性质最直接和最可靠的方法.然而，我们尚未对qS波进行处理和解释，联合使用qP和qS波进行反演可以对介质的各向异性作出更全面的、完整的解释，这将是我们下一步研究的方向.

刘宪彬，郑需要.2013.利用qP波慢度和偏振矢量计算弱各向异性介质参数［J］.地震学报，35（2）：184-198.

Liu X B，Zheng X Y.2013.Calculation of weakly anisotropic parameters by slowness and polarization vectors of qP wave［J］.Acta Seismologica Sinica，35（2）：184-198（in Chinese）.

袭杨.2009.空间直线拟合的一种方法［J］.齐齐哈尔大学学报，25（2）：64-68.

Xi Y.2009.A method for fitting of a space straight line［J］.Journal of Qiqihar University，25（2）：64-68（in Chinese）.

郑需要，董天立，刘一霖.2010.变井源距垂直地震剖面各向异性参数反演［J］.地震学报，32（2）：157-166.

Zheng X Y，Dong T L，Liu Y L.2010.Inversion of seismic anisotropic parameters for walkaway VSP data［J］.Acta Seismologica Sinica，32（2）：157-166（in Chinese）.

ˇCervenýV.2001.Seismic Ray Theory［M］.Cambridge：Cambridge University Press：189-199.

Horne S A，Leaney W S.2000.Polarization and slowness component inversion for TI anisotropy［J］.Geophys Prospect，48（4）：779-788.

Jech J，Pšenˇcík I.1989.First-order perturbation method for anisotropic media［J］.Geophys J Int，99（2）：369-376.

Pšenˇcík I，Gajewski D.1998.Polarization，phase velocity and NMO velocity of qP waves in arbitrary weakly anisotropic media［J］.Geophysics，63（5）：1754-1766.

Thomsen L.1986.Weak elastic anisotropy［J］.Geophysics，51（10）：1954-1966.

Zheng X，Pšenˇcík I.2002.Local determination of weak anisotropy parameters from qP-wave slowness and particle motion measurements［J］.Pure Appl Geophys，159（7／8）：1881-1905.

Zheng X.2004.Inversion for elastic parameters in weakly anisotropic media［J］.Geophys J Int，159（3）：1077-1089.

Zheng X.2007.Identification of the symmetry planes in weakly anisotropic elastic media［J］.Geophys J Int，170（1）：468-478.

附录 qP波弱各向异性矩阵元素和弱各向异性参数坐标变换

1 qP波弱各向异性矩阵BMN的表达式

式中ni为波传播的方向.

15个WA参数与密度归一化的弹性参数的关系为

2 各向异性参数的坐标变换