基于VSP数据的任意弱各向异性参数反演方法

陈占国 曾昭翰 杨心超

(中国石化石油物探技术研究院，江苏南京 211103)

0 引言

随着地震各向异性理论的发展和应用，工业界基于VTI介质的地震建模和成像技术日趋成熟； TTI各向异性建模和成像技术在实际中得到较广泛的应用[1]； 学术界已开始研究如基于正交各向异性理论等的更复杂应用技术。

地震各向异性参数的获取一直是各向异性应用的关键。与地面地震相比，VSP测量作为各向异性参数获取的一种途径，具有较高的精度和准确的深度优势。前人研究表明，应用 VTI模型可以更加精确地描述长排列上的地震波旅行时特性[2]； 基于 VTI假设，给出了不同形式的时差关系：Byun等[3-4]应用一种双参数速度分析对VSP合成数据进行速度分析； Hake等[5]导出了VTI介质时差的3项泰勒展开式； Tsvankin等[6-7]对Hake的公式进行了改进，并引入渐近拟合使新公式在大炮检距上更逼近真值； Alkhalifah等[8-9]提出用两个参数表达的非双曲时差关系描述 VTI介质的地震波旅行时，已在地面地震数据各向异性动校正等处理中得到了广泛的应用[10-12]； 白海军等[13]和陈沅忠等[14]将Alkhalifah的方法应用于VSP数据，均取得了较好的效果。

以上方法主要是基于VTI或TTI介质的纵波非双曲线时差进行速度分析和旅行时反演，难以获得更复杂的任意各向异性参数。VSP能够直接测量地震波的慢度和极化信息，且不受上覆地层影响。极化和慢度矢量是表征Christoffel方程的直接参数，理论上只要获得足够多的慢度和极化信息，就可以反演介质的所有各向异性参数。根据这一特点，Jech等[15]提出了各向异性介质的一阶扰动方法。Horne 等[16]和Gomes等[17]基于扰动理论，提出了应用慢度和极化矢量求解介质弱各向异性参数的方法，并得到了一些成功的应用。

在前人研究的基础上，本文系统地介绍了基于慢度和极化信息反演各向异性参数的方法原理，并给出了Walkaway VSP慢度和极化数据的获取方法，提出了基于加权迭代的任意弱各向异性参数反演算法的详细实现过程，应用一组TTI模型合成数据和一口实际井资料进行处理试验，取得了理想的效果。

1 方法原理

1.1 Christoffel 方程及扰动方法

介质中传播的平面波遵循Christoffel方程[18]

(1)

式中：Γjk=aijklninl为Christoffel矩阵，其中a是密度归一化弹性张量，n为波的传播方向矢量，i、j、k、l均取1～3，表示坐标方向；m=1、2、3，分别代表各向异性介质中传播的qS1、qS2、qP波；Gm是Christoffel矩阵的三个特征值；g(m)为三种波的极化矢量。式(1)中使用了爱因斯坦求和约定。

由于实际介质大多是弱各向异性[19-20]，基于扰动理论，Christoffel传播方程可近似表述为[21]

(2)

式中：线“～”表示在各向同性参考介质中的基本量； “Δ”表示在弱各向异性介质中的扰动量。

1.2 qP波慢度和极化矢量正反演公式

只考虑式(2)中的qP波(m=3)，忽略其中二阶项，并引入了三个相互垂直的单位矢量[22]

(3)

(4)

(5)

式中：p和g分别是qP波的慢度和极化矢量；α和β分别表示参考介质中qP波和qS波的相速度，K′=1、2；B为与弱各向异性参数相关的弱各向异性矩阵[25]，可表示为

(6)

其中Δa为密度归一化弹性张量的扰动量。

反过来，通过慢度矢量和极化矢量，可以求解弱各向异性介质的弹性参数[23-24]

(7)

B33=-2α4piΔpi

(8)

(9)

式(7)和式(8)是关于15个弱各向异性参数的线性公式。

1.3 qP波慢度矢量和极化矢量的计算

1.3.1 极化矢量

在数学上，质点的运动轨迹可以认为是空间点的集合。对于三分量的VSP记录，这种离散采样点的质点运动轨迹具有相关性，可以通过三分量记录的协方差矩阵统计其空间运动规律。根据矩阵的特征值和特征向量可定量地确定极化椭球(长、中、短)各轴的大小和椭球的空间位置及其标准椭球方程，并能求出质点运动轨迹的极化参数。协方差矩阵定义为

(10)

求解协方差矩阵的特征值和对应的特征向量，可以确定极化方向，其中最大特征值对应的特征向量就是qP波的极化矢量。

1.3.2 慢度矢量

一般情况下，可通过邻近两个检波器的时差dt和深度差dh求取慢度的垂向分量

(11)

慢度p是一个局部量，不依赖于上覆地层速度结构。根据地震波走时的互易性，如果介质是横向均匀的，则炮点到接收点的走时等于接收点到炮点的走时。当介质横向非均匀性比较弱时可以将接收点与炮点互易，利用共接收点的旅行时计算慢度矢量的水平分量。

旅行时曲线一般情况下都不光滑，其一阶导数不连续甚至不存在。为了得到更稳定的慢度矢量，本文采用多点平滑和相邻五点三次多项式求导的方法计算慢度矢量。

1.4 Walkaway VSP反演弱各向异性参数方法

通过将BKL代入式(7)和式(8)，可以建立弱各向异性参数与慢度及极化矢量的关系式。如果慢度矢量的扰动Δp的三个分量和极化矢量g(3)的三个分量都已知，可以单独使用式(7)或式(8)建立一个关于未知数为弱各向异性参数的线性方程组，或者联合式(7)和式(8)建立方程组，求解该线性方程组即可得弱各向异性参数。

(12)

式中：N是数据点个数；y表示式(7)和式(8)的观测参数和已知量；w为观测参数的权系数；x为待反演的弱各向异性参数，其维数为Npar；J=1,2,…,Npar；A为N×Npar阶系数矩阵。弱各向异性参数的估计可以通过奇异值分解求解，即

x=VΛ-1UTd

(13)

式中：U和V是系数矩阵ATwA两个正交的N×N奇异矩阵；Λ是由ATwA的奇异值构成N×N对角矩阵；d=ATwy。

权系数w由数据的初始权重wd和更新后的测量值与正演计算值残差的权重wr组成

(14)

数据的初始权重wd由两部分构成，一部分是参与运算的数据类型的质量估计权系数wq，另一部分权系数由慢度或极化数据的可靠性决定，即

(15)

权系数wr为反演迭代时慢度或者极化矢量残差的再加权[26]

(16)

式中：drI为第I个慢度或者极化矢量的观测值与正演计算值的残差；σmed=0.67449，为高斯噪声的中值；ω是dr的标准偏差截断因子，一般选择3～6间的数； drmed定义为

drmed=med(|drI-med(drI)|)

(17)

其中med(•)表示取数列的中值。

求解式(13)后，每个模型参数的最小二乘估计误差EJ为[26]

(18)

式中：TJJ为奇异值和奇异矩阵构成的协方差矩阵T=VΛ-2VT的对角元素；Ω为加权剩余量的变化

(19)

独立的弱各向异性参数共有21个，其中与qP波相关的有15个：εx、εy、εz、δx、δy、δz、χx、χy、χz、ε15、ε16、ε24、ε26、ε34和ε35。若要反演15个参数，需要布设5条以上Walkaway VSP测线； 如果对两个正交的测线进行联合反演，最多可得到9个参数：εx、εy、εz、δx、δy、ε15、ε24、ε34和ε35； 如果仅使用一条测线，只能反演5个参数，即εx、εz、δx、ε15和ε35[25]。

2 模型资料测试

模型是一个弹性参数随深度线性变化的TTI模型，相速度的各向异性程度(相速度的最大最小值之差与最大值的比值)为4%。地表处的TTI弹性参数(密度归一化的)矩阵Θ3和10km深处的弹性参数矩阵Θ4分别由相应的两个VTI弹性矩阵Θ1和Θ2经过绕y轴顺时针旋转40°、再绕z轴(向下为正)顺时针旋转30°后得到。四个弹性参数矩阵分别为

(20)

(21)

(22)

(23)

观测系统由6条等方位间隔的Walkaway VSP测线组成。每条测线布设32个炮点，井孔两边各16个； 井源距范围是100～3100m，炮间距为200m。井中1000～1550m深度等间距布设12个检波器，检波器间距为50m。

联合应用式(7)和式(8)，通过6条Walkaway VSP的慢度矢量和极化矢量的3个分量，反演了模型的15个弱各向异性参数。

以深度1250m的检波点为例，其他深度结果类似。图1展示了方位角为30°/210°、90°/270°和150°/330°测线测量的慢度和极化矢量随井源距的变化曲线，可见测量的慢度/极化矢量与真值完全吻合。图2为测量的极化矢量和慢度矢量的相对误差，可见极化矢量的相对误差非常小，慢度矢量的水平分量误差稍大一些，但基本都小于1%。

图3为6条测线联合反演的不同深度15个弱各向异性参数结果。反演的弱各向异性参数(散点)在不同深度上有所不同，但基本都在真值(直线) 的附近，其绝对误差(图4)最大为4.05×10-3，远小于最小的各向异性参数，说明模型的反演结果精度很高。

图5是利用反演的弱各向异性参数通过式(4)和式(5)正演计算的不同方位测线(30°/210°、 90°/270°和150°/330°)1250m深度的慢度和极化矢量与测量值的对比，可以看出二者高度吻合。图6是利用式(4)正演的不同方位、不同倾角在1250m深度处接收点的相速度曲线，可见：120°/300°测线的相速度曲线对称，说明这个方位是TTI对称平面的垂直方向，那么对称平面方位角则为30°/210°； 在对称平面上，-40°倾角时相速度最小，说明这是对称轴倾角，通过对称平面上的相速度两个极值计算的各向异性程度为4%，与TTI模型的设计值完全一致。表1给出了所有深度点计算的各向异性程度和TTI对称轴方向，与模型设计值高度一致。

图1 3条Walkaway测线测量的慢度、极化矢量与真值对比(a)30°/210°； (b)90°/270°； (c)150°/330°。上标“obs”表示测量值，“r”表示真值

图2 3条Walkaway测线测量的慢度、极化矢量的相对误差(a)30°/210°； (b)90°/270°； (c)150°/330°

图3 弱各向异性参数反演结果与真值的对比(a)εx、δx、χx、ε15、ε24； (b)εy、δy、χy、ε16、ε35； (c)εz、δz、χz、ε26、ε34。上标“i”表示反演值，“r”表示真值

图4 反演的弱各向异性参数的绝对误差条形图

图5 3条Walkaway测线正演计算的慢度、极化矢量与测量值的对比(a)30°/210°； (b)90°/270°； (c)150°/330°。上标“obs”表示测量值，“cal”表示计算值，图9、图10同

表1 不同深度接收点处反演的各向异性程度和TTI对称轴方向

图6 不同方位测线速度随倾角的变化曲线

3 实际资料应用

以A井3条Walkaway VSP采集数据为例。测线布置如图7，其中测线1布设450炮、测线2布设460炮、测线3布设470炮。接收点深度范围为4560～5670m，共75 级，间距为15m。

基于上述3条测线，可反演9个弱各向异性参数，其中5个参数(εx、εz、δx、ε15和ε35)结果比较稳定。

图7 实际资料测线布置示意图

图8a和图8b分别是反演的速度和5个弱各向异性参数随深度的变化曲线。在地层分界面附近，尤其是第二和第三界面处，相速度有较大的抖动(图8a)，反演估计的绝对误差也较大(图8c)，可能是由于在地层分界面上、下的接收点时差变化较大，造成垂直慢度的误差较大。在盐丘(第三条虚线以下)中各向异性参数有明显的变大趋势，说明盐丘各向异性很强。从图8c可发现，几个ε参数的反演误差都比δx小，这是因为ε只与qP波各向异性有关，δ不仅与qP波有关，还与qS波有关，而本文反演只用到了qP波信息，但反演误差的最大值为2.296×10-4，仍远小于弱各向异性参数值。

图9和图10分别是应用反演的弱各向异性参数正演计算的3条测线的第一个接收点的慢度和极化矢量与测量值的对比，可以看出，测线1的慢度和极化矢量的计算值与测量值都吻合的较好，测线2和测线3基本能吻合，说明反演的弱各向异性参数正确地反映了测量的慢度和极化矢量所携带的地层各向异性信息。

图11是计算的不同深度的各向异性程度及其各层的平均各向异性程度。如图中所示，测线1的各向异性程度较弱，测线2、测线3较强，其中测线2最强。从深度上看，深度5500m以下盐丘各向异性明显增强，其中测线2最大，超过18%。

图8 反演的速度(a)、弱各向异性参数(b)及其误差(c)随深度的变化曲线水平虚线表示地层界面

图9 应用反演的弱各向异性参数计算的慢度矢量与测量值的对比(a)测线1； (b)测线2； (c)测线3

图10 应用反演的弱各向异性参数计算的极化矢量与测量值的对比(a)测线1； (b)测线2； (c)测线3

图11 3条测线不同深度的各向异性程度曲线(a)与地层平均各向异性程度(b)的对比

4 结束语

本文系统地介绍了基于扰动理论的qP波的慢度和极化矢量正、反演方法，给出了精度高、稳定的慢度和极化矢量计算方法，提出了适用于慢度和极化矢量反演任意弱各向异性参数的加权迭代反演算法。

通过设计6条Walkaway VSP的TTI模型处理试验，反演的15个弱各向异性参数精度很高； 应用反演的弱各向异性参数计算的慢度和极化矢量与测量值完全匹配； 正演的相速度结果表明，模型各向异性程度和对称轴与理论值完全一致。

实际3条测线Walkaway VSP资料反演的5个弱各向异性参数取值范围和各向异性程度都符合地质认识，各向异性随深度变化趋势与地层分界面对应较好； 应用弱各向异性参数计算的慢度和极化矢量与测量值吻合很好； 反演误差远小于各向异性参数的数量级。

模型和实际应用表明，本文方法能够用于高精度的地震各向异性参数提取，具备较高的应用价值。