利用qP波慢度和偏振矢量计算弱各向异性介质参数

刘宪彬 郑需要

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

引言

本文提出了一个模型试验来验证Zheng和Pšenˇcík（2002）给出的反演公式.与Gaiser（1990）和Horne等（1998）使用的方法相比，本试验更具有普遍意义.在实验中，使用多震源垂向地震剖面（vertical seismic profile，简写为VSP）数据反演介质的局部各向异性参数，三分量地震记录数据不受自由界面的影响.所研究的介质是一种垂向非均匀弱各向异性介质，具有任意的对称性.对于这样的各向异性介质，我们寻找一种参考各向同性介质，各向异性介质与参考各向同性介质略有不同.偏振矢量可以很直接地从三分量地震图中得到，然而，慢度矢量的确定却比较复杂.在此，我们使用Gaiser（1990）提出的方法，钻孔中慢度矢量的垂直分量由在钻孔中记录到的走时确定，因为假设介质是横向均匀的，所以能够利用走时的互易性确定钻孔中慢度矢量的水平分量.应该指出，Gaiser（1990）仅仅考虑具有垂直对称轴的横向各向同性介质（transversely isotropic media with vertical axis of symmetry，简写为VTI），慢度矢量总是被限制在有震源和钻孔所决定的垂直平面内，而在我们的模型中介质可以具有任意对称性，慢度矢量可以背离垂直平面.如果qP波的慢度矢量和偏振矢量是已知的，那么，慢度公式和偏振公式形成一套线性代数方程组，使用这套方程组可以确定钻孔中检波器周围介质的各向异性参数.我们进一步把垂向非均匀各向异性介质扩展到横向非均匀情况，从慢度公式和偏振公式中消掉慢度矢量的水平分量，得到了适用于仅有一个慢度分量时的反演公式，并用一条实际的VSP观测资料进行验证.

1 基本方法

使用一个右手笛卡尔坐标系，x轴和y轴在水平面内，z轴垂直向下为正方向.在弱各向异性（weakly anisotropic，简写为WA）介质中，可以把P波的慢度c－1和偏振矢量gi表示为（Zheng，Pšenˇcík，2002；Zheng，2004；郑需要等，2010）

式（1）和式（2）中，pk是WA介质中的慢度矢量；pkpk＝c－2；矢量e（1）i，e（2）i，ni＝e（3）i分别为定义在接收点的各向同性介质中沿P波射线的3个相互垂直的单位矢量；α和β分别为各向同性介质中的P波和S波的速度；矢量ni是射线的切线方向；单位矢量e（1）i和e（2）i可以在垂直于ni的平面内任意选取.但下面的选择为公式的推导和计算带来很大的方便（Pšenˇcík，Gajewski，1998；Zheng，Pšenˇcík，2002）

式中，D＝（n21＋n22）1／2，n21＋n22＋n23＝1.

在波的传播方向，n＝（cosφsinθ，sinφsinθ，cosθ），其中，φ是方位角，而θ是波的传播方向与z轴的夹角（0≤φ≤2π，0≤θ≤π）.所以，D＝sinθ，而且e（1），e（2）可写为

为了避免D＝0或者非常接近于零的时候式（3）中的分母为零的问题，我们选择使用式（4）.

式（1）和式（2）中的B13，B23和B33被称为弱各向异性矩阵元素，它们是15个 WA参数的函数.这15个WA参数可以完全描述弱各向异性介质中qP波的性质（Pšenˇcík，Gajewski，1998；Zheng，Pšenˇcík，2002）.重写式（1）和式（2），使包含 WA参数弱各向异性矩阵元素B13，B23和B33出现在方程的左边，它们是待求解线性方程组中的未知参数

式（5）和式（6）代表由15个WA参数作为未知参数的线性方程组.方程组等号右端的偏振矢量gk和慢度矢量pk可以从观测中得到.为方便起见，将式（5）称作慢度公式，式（6）称作偏振公式.

2 合成数据实验

考虑一个合成的多方位、多震源VSP实验.使用右手笛卡尔坐标系，x轴和y轴在水平面内，z轴正方向朝下.对于钻孔中的每一个接收点，如果知道了偏振矢量和慢度矢量，便可得到3 Ns×Np个方程.其中，Ns为沿一条剖面的震源数目，Np为具有不同方位的剖面数目.要想由qP波唯一地确定15个各向异性参数，至少需要5条剖面（Pšenˇcík，Gajevski，1998）.模型位于立方体中，其尺度参数为－1.0≤x≤1.0，－1.0≤y≤1.0，－1.0≤z≤1.0（单位：km）.介质为VTI介质，由下面的两个密度归一化的弹性矩阵说明.在地表z＝0km处，弹性参数（单位：km2／s2）为

在z＝1km处，弹性参数为VTI介质的对称轴先绕y轴（离开垂直轴）转动80°，然后再绕x轴转动25°，最后得到的两个矩阵的所有元素都不为零，可以把它们看作是一般的弹性矩阵.

钻孔坐落在坐标系的原点，在钻孔里等间隔地布设13个三分量地震仪，间距为0.05km，深度为0.1—0.7km，忽略自由地表效应.在地表跨过钻井口布设6条剖面，方位（从x轴开始逆时针方向）依次为0°，30°，60°，90°，120°和150°.在钻孔一边每一条剖面上有9个震源，最近的震源距钻孔0.1km，相邻震源的距离为0.1 km（图1）.在这一观测系统中，每一个接收点有108条射线.使用各向异性介质中的射线追踪方程，可以计算出震源到接收点之间的射线.在模型的上部，射线覆盖相当好，最浅的接收点不仅被近水平的射线照亮，也被近垂直的射线照亮.最深的接收点的射线覆盖不如浅部的接收点覆盖得那么好.它们主要是被近垂直的射线照亮，没有近水平的射线照在这些接收点上.因此，我们不能期望在深的接收点处得到好的反演结果.“观测的”三分量合成地震图是由修改版的 Anray软件包计算得到的（Gajewski，Pšenˇcík，1990）.

从“观测的”合成地震图上可以获得qP波的到时和偏振矢量.在每一个接收点，使用到时决定P波的参考速度α，S波的速度为β＝.利用相邻接收点的到时决定慢度矢量的垂直分量，慢度矢量的径向分量根据走时的互易性决定，切向分量很小而被忽略.根据每一接收点的参考速度，可以建立背景各向同性模型，它是一个一维垂向非均匀速度模型.在这个模型里进行射线追踪，可以得到每个接收点处的，和ni；也可以用接收点处的偏振矢量近似地确定这些矢量.后者与上覆介质无关，可以避免复杂的射线追踪.由此便可得到式（5）和式（6）左边所有的系数和右边的观测值.因为式（5）和式（6）是线性方程组，可以使用奇异值分解的方法进行求解反演.

基于6条剖面的数据，得到了图2所示的反演结果.图2a，b，c各选择了3个深度的接收点的反演结果，其深度值依次为0.25，0.45和0.65km.在每一个接收点，左边三列给出了精确的对称轴的投影（实线）和反演得到的结果（虚线）；最右边一列给出了包含对称轴的相速度曲线，水平坐标表示波的传播方向与对称轴的夹角，从0°（对称轴方向）到90°（各向同性面），垂直坐标表示相速度.可以看出，反演结果与精确解非常接近，慢度公式与偏振公式的联合反演结果好于使用单一公式的反演结果.正如所期望的那样，最深的接收点的反演结果不如最浅的接收点的反演结果好，因为最深的接收点没有很好地被射线照亮.

当使用3条剖面时，反演结果如图3所示.3条剖面的方位角分别为0°，60°和120°.如前所述，3条剖面少于完全恢复所有各向异性参数所需的剖面的数量（5条）.可以看出，反演结果与使用6条剖面的结果相比稍差一些，特别是对称轴在（x，y）平面内的投影和相速度的变化形式，然而相速度与真速度的基本特征仍然得以保留.慢度公式反演结果比较好，联合反演也产生了较好的结果.

图1 多方位多离源距的VSP实验示意图Fig.1 Schematic diagram of multi－－azimuthal and multi－－source offset VSP experiment

图2 使用6条剖面资料反演各向异性参数的结果（a）由慢度公式（5）反演得到；（b）由偏振公式（6）反演得到；（c）由慢度公式和偏振公式联合反演得到.各子图的上、中、下3行对应的接收点深度分别为0.25，0.45和0.65km.左边三列给出精确的对称轴投影（实线）和反演结果（虚线）；最右边一列给出包含对称轴的相速度曲线，水平坐标表示波的传播方向与对称轴的夹角，垂直坐标表示相速度，实线为相速度的准确值，虚线为反演结果Fig.2 Inversion for WA parameters by using data from 6profiles（a）Obtained with slowness equation（5）；（b）Obtained with polarization equation（6）；（c）Obtained by combining equations（5）and（6）.The receiver depths for the upper，middle and lower rows in each plot are 0.25，0.45and 0.65km，respectively.The first three frames from the left show the exact（solid line）and inverted（dashed line）projections of axis of symmetry；the frame on the right shows phase velocity sections within the symmetry plane，in which the horizontal coordinate represents angles between wave propagation direction and axis of symmetry，vertical coordinate represents phase velocity，solid line denotes the exact value of phase velocity and dashed line denotes inversed result

图3 使用3条剖面资料反演各向异性参数的结果.（a）由慢度公式（5）反演得到；（b）由偏振公式（6）反演得到；（c）由慢度公式和偏振公式联合反演得到.图注说明同图2Fig.3 Inversion for WA parameters by using data from 3profiles.（a）Obtained with slowness equation（5）；（b）Obtained with polarization equation（6）；（c）Obtained by combining equations（5）and（6）.The illustration is the same as for Fig.2

下面的结果和图形显示了不同水平的噪声对数据及反演结果的影响.首先在观测数据中加入10%的随机噪声.图4和图5分别显示了使用6条剖面和3条剖面的反演结果.在图4中反演得到的对称轴与精确的对称轴有一点点偏离，其它的与图2所示结果没有明显差别.图5显示噪声对于对称轴在（x，y）平面的投影有较强的影响，反演得到的相速度剖面与精确解仍然比较吻合.

图4 使用6条剖面资料、数据中加入10%随机噪声后的反演结果.（a）由慢度公式（5）反演得到；（b）由偏振公式（6）反演得到；（c）由慢度公式和偏振公式联合反演得到.图注说明同图2Fig.4 Inversion for WA parameters by using data with 10%random noise from 6profiles.（a）Obtained with slowness equation（5）；（b）Obtained with polarization equation（6）；（c）Obtained by combining equations（5）and（6）.The illustration is the same as for Fig.2

图5 使用3条剖面资料、数据中加入10%随机噪声后的反演结果.（a）由慢度公式（5）反演得到；（b）由偏振公式（6）反演得到；（c）由慢度公式和偏振公式联合反演得到.图注说明同图2Fig.5 Inversion for WA parameters by using data with 10%random noise from 3profiles.（a）Obtained with slowness equation（5）；（b）Obtained with polarization equation（6）；（c）Obtained by combining equations（5）and（6）.The illustration is the same as for Fig.2

图6 使用6条剖面资料、数据中加入20%随机噪声后的反演结果.（a）由慢度公式（5）反演得到；（b）由偏振公式（6）反演得到；（c）由慢度公式和偏振公式联合反演得到.图注说明同图2Fig.6 Inversion for WA parameters by using data with 20%random noise from 6profiles.（a）Obtained with slowness equation（5）；（b）Obtained with polarization equation（6）；（c）Obtained by combining equations（5）and（6）.The illustration is the same as for Fig.2

当观测数据中加入了20%的随机噪声时，从6条剖面得到的反演结果仍然很好（图6）.这种水平的噪声对使用3条剖面的反演结果有一些影响（图7）.如果噪声水平达到30%，可以看到使用6条剖面的反演结果（图8）和使用3条剖面的反演结果（图9）都与精确解有较大的偏离.使用慢度公式得到的结果好于使用偏振公式的结果，偏振公式对联合反演也有较大的影响.使用6条剖面和慢度公式得到的结果仍然是可以接受的，使用偏振公式和联合反演的结果有非常大的畸变.进一步增大噪声水平将导致更严重的畸变.

图7 使用3条剖面资料、数据中加入20%随机噪声后的反演结果.（a）由慢度公式（5）反演得到；（b）由偏振公式（6）反演得到；（c）由慢度公式和偏振公式联合反演得到.图注说明同图2Fig.7 Inversion for WA parameters by using data with 20%random noise from 3profiles.（a）Obtained with slowness equation（5）；（b）Obtained with polarization equation（6）；（c）Obtained by combining equations（5）and（6）.The illustration is the same as for Fig.2

图8 使用6条剖面资料、数据中加入30%随机噪声后的反演结果.（a）由慢度公式（5）反演得到；（b）由偏振公式（6）反演得到；（c）由慢度公式和偏振公式联合反演得到.图注说明同图2Fig.8 Inversion for WA parameters by using data with 30%random noise from 6profiles.（a）Obtained with slowness equation（5）；（b）Obtained with polarization equation（6）；（c）Obtained by combining equations（5）and（6）.The illustration is the same as for Fig.2

图9 使用3条剖面资料、数据中加入30%随机噪声后的反演结果.（a）由慢度公式（5）反演得到；（b）由偏振公式（6）反演得到；（c）由慢度公式和偏振公式联合反演得到.图注说明同图2Fig.9 Inversion for WA parameters by using data with 30%random noise from 3profiles.（a）Obtained with slowness equation（5）；（b）Obtained with polarization equation（6）；（c）Obtained by combining equations（5）and（6）.The illustration is the same as for Fig.2

3 实际数据反演

把上面的反演方法应用到在爪哇海（Java Sea）地区获得的一条walkaway VSP剖面中（Horne，Leaney，2000）.剖面上布设了228个震源，分布在钻孔两侧－2.5—2.5km的范围内，相邻震源间隔约为25m.以钻孔1.63km为中心布设5个三分量检波器，相邻检波器间距为5m.数据由下行和上行（反射）qP波和qS波组成，由Esmersoy（1990）和Leaney（1990）提出的参数反演方法提取.本文仅仅研究qP波问题，并且假定地下介质是横向非均匀、弱各向异性介质.在这种情况下只能从观测记录中得到慢度矢量的垂向分量和偏振矢量.通过在式（5）和式（6）中消去慢度矢量的水平分量，可以得到由慢度矢量的垂向分量和偏振矢量表示的反演公式（Zheng，Pšenˇcík，2002）.因为只有一条剖面，所以仅能反演5个 WA参数，即εx，εz，δx，ε15和ε35（Zheng，Pšenˇcík，2002）.这5个 WA参数完全决定介质在（x，z）平面内的性质.

从walkaway VSP剖面中得到了228个下行qP波和228个上行qP波数据（图10）（Horne，Leaney，2000）.图10中纵轴表示慢度矢量的垂直分量，横轴表示偏振角.偏振角在－90°―90°之间的数据为下行波，偏振角在±（90°―180°）之间的数据为上行波（反射波）.从图10中可以看出，上行波比下行波有较大的离散性.图中实线是联合使用上行波和下行波反演得到的理论结果.在反演中假定介质是一般各向异性介质.

图10 来自爪哇海地区的一条变井源距观测剖面Fig.10 A multi-source offset VSP experiment from the Java Sea region

在qP波偏振角0°—75°范围内可用的下行波和上行波数目分别为133个和170个.偏振角大于75°的数据误差较大，因为在这个范围内，慢度矢量和偏振矢量的垂直分量较小，信噪比较低.可以单独使用下行波组成133个线性方程组，也可以联合使用下行波和上行波组成303个线性方程组.对前面描述的由5个WA参数组成的方程进行求解，可以获得相关的WA参数.使用上面公式得到的参考介质的P波速度为α＝2.75km／s.图11给出了介质为一般各向异性介质模型时得到的（x，z）平面内的相速度曲线.图中也给出了Horne和Leaney（2000）使用VTI介质模型和不同反演方法得到的结果.当介质为各向同性时，相速度为常数，在图11中显示为一条直线.Horne和Leaney（2000）假定介质是VTI介质，单独使用慢度数据和联合使用慢度与偏振数据得到了两个稍有差别的结果，相速度在水平方向较为一致，但是在垂直方向有较大的区别.图11中的细实线和粗实线分别表示使用下行波和联合使用下行波和上行波得到的相速度曲线.与Horne和Leaney（2000）结果的主要区别在于：他们得到介质的各向异性为14%（单独使用慢度）和19%（联合使用慢度和偏振）；而本文通过求解前面描述的反演方程得到的介质各向异性为8%.本文得到的相速度的最小值也不在垂直方向，而是偏离垂直方向大约20°，这一方向可能是对称轴在（x，z）平面内的投影.

图11 用不同方法得到的反演结果Fig.11 Results from inversion using different methods

4 讨论与结论

合成数据实验表明，本文提出的反演方法能够确定具有任意对称性的弱各向异性参数，尽管在反演时假定慢度矢量的切向分量为零.使用较少的剖面表明联合使用慢度公式和偏振公式可以得到较好的结果.噪声实验表明，慢度公式是确定各向异性参数的较为可靠的工具.

显然，肯定和证实上述结论需要进一步的实验.笔者希望继续进行具有不规则震源剖面分布和减少震源数目的实验，探讨介质横向非均匀性和低对称性介质对反演的影响.我们计划研究那些使慢度矢量严重偏离传播面的地震各向异性介质.

上面的实验可以很容易地推广到钻孔为非垂直的和震源高程不同的情况.本文已经将介质垂向非均匀性的条件扩展到了介质的横向非均匀性，由于使用了慢度分量和偏振矢量，使得在反演中不用考虑上覆介质的影响，避免了复杂的射线追踪计算.除了介质是弱各向异性的假设条件外，本文提出的正演公式和反演公式具有普遍性.

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

Esmersoy C.1990.Inversion of P and SV waves from multicomponent offset vertical seismic profiles［J］.Geophysics，55（1）：39－－50.

Gaiser J E.1990.Transversely isotropic phase velocity analysis from slowness estimates［J］.J Geophys Res，95（B7）：11241－－11254.

Gajewski D，Pšenˇcík I.1990.Vertical seismic profile synthetics by dynamic ray tracing in laterally varying layered anisotropic structures［J］.J Geophys Res，95（B7）：11301－－11315.

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

Horne S A，McGarrity J P，Sayers C M，Simith R L，Wijnands F.1998.Fractured reservoir characterization using multiazimuthal walkaway VSPs［C］∥Expanded Abstracts of the 68th Annual International Meeting of the Society of Exploration Geophysics.Tulsa：1640－－1643.

Leaney W S.1990.Parametric wavefield decomposition and application［C］∥Expanded Abstracts of the 68th Annual International Meeting of the Society of Exploration Geophysics.Tulsa：1097－－1100.

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.

Zheng X Y，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 Y.2004.Inversion for elastic parameters in weakly anisotropic media［J］.Geophys J Int，159（3）：1077－－1089.