三维VSP地震初至波走时层析成像研究

黄光南， 邓居智， 李红星， 李泽林， 张 华， 王安东

(1.东华理工大学放射性地质与勘探技术国防重点学科实验室，江西 南昌 330013；2.中国地质大学(武汉)地球内部多尺度成像湖北省重点实验室，湖北 武汉 430074；3.中国石油大学(北京)油气资源与探测国家重点实验室，北京 102249)



三维VSP地震初至波走时层析成像研究

黄光南1,2,3， 邓居智1， 李红星1， 李泽林1， 张 华1， 王安东1

(1.东华理工大学放射性地质与勘探技术国防重点学科实验室，江西 南昌 330013；2.中国地质大学(武汉)地球内部多尺度成像湖北省重点实验室，湖北 武汉 430074；3.中国石油大学(北京)油气资源与探测国家重点实验室，北京 102249)

在油气田勘探开发阶段，VSP(垂直地震剖面)作为一种精细勘探技术发挥着重要作用，三维VSP地震初至波走时层析成像可以直接得到地层速度模型。它结合了三维VSP地震射线追踪算法和数学反演算法，它能够实现对地下地层速度模型进行重构，速度成像结果对认识地质结构提供了直观的依据。为了测试算法的稳定性和抗干扰能力，对三维VSP正演走时加入了20%的高斯噪声，然后利用三维地震初至波走时层析成像进行速度反演，得到了三维地层速度模型、垂直切片图和水平切片图。射线方向的交叉性程度决定了反演结果的质量和可靠性，由于射线在浅部地层的交叉性比深部地层的交叉性更好，浅部地层的反演效果更加清晰；相反，深部地层的速度反演效果较差。算法对含噪声数据取得了良好的反演结果，验证了该算法具有较好的抗干扰能力。

垂直地震剖面；高斯噪声；射线追踪；地震层析成像

黄光南，邓居智，李红星，等.2016.三维VSP地震初至波走时层析成像研究[J].东华理工大学学报：自然科学版，39(3)：266-272.

Huang Guang-nan, Deng Ju-zhi, Li Hong-xing,et al.2016.3D VSP seismic first arrival traveltime tomography[J].Journal of East China University of Technology (Natural Science), 39(3)：266-272.

VSP是一种井-地地震观测技术。与地表采集数据相比,VSP观测方式的采集数据具有更高的信噪比和分辨率。VSP观测方式能够提供地下地层结构和地面测量参数之间的直接对应关系，这种关系有助于地震资料处理实施时深转换，并提供精确的地层速度模型。地震走时层析成像是利用大量地震观测走时数据反演地层物性参数(纵波速度、横波速度、品质因子和各向异性参数等)分布的一种地球物理反演方法。它的反演结果能够以图像的方式非常清晰、直观地显示地层物性参数分布，所以VSP技术一出现就备受地球物理学者的重视。

地震走时层析成像在天然地震学和地球动力学等重大科学领域具有广泛应用价值。在天然地震学领域，地震台站不间断地接收不同的地震事件，利用这些地震事件可以反演出地球区域范围内的地层速度分布。根据地层速度反演结果，可以进一步反演天然地震震源位置信息和分析地震的孕震机制(Aki et al.，1976；徐涛等，2014)。在地震发生的短时间内，震源位置信息对于防震减灾工作具有非常重要的意义。利用地震事件反演得到的地层速度分布信息，还可以用于发现和调查大洋洋中脊分布、大陆板块消减带、地壳内部的火山分布、大型断裂带、地幔与地壳的分界面等重大科学问题(Calvert et al.，2001；Lees，1992；Rawlinson et al.，2012；Zhao et al.，2013)。另外，全球地层速度分布信息可以进一步解释地球内部的动力学机制，地质学家认为低速、低密度的新物质产生于洋中脊，而高速、高密度的物质消失于洋壳与陆壳的接触带，这种动力学机制进一步为板块构造学说提供了解释依据。地震走时层析成像在石油、煤田和矿产等资源勘探领域也具有实际应用价值。随着全球经济的快速发展，社会对能源的需求急剧增加，这种社会背景巨大地促进了地球物理勘探技术的发展，地震走时层析成像技术也得到了前所未有的发展(和锐等，2007；孙章庆等，2012；Zhang et al.，1998；Zhou，2006)。在油气地震勘探时，经常将地震资料偏移结果、层析成像剖面和测井资料相结合，对油气藏等地下异常体的特征进行综合判断(常旭等，1996；李录明等，2013；秦宁等，2011)。煤炭等矿产资源与围岩的速度存在较大差异，它们在地震层析成像剖面上通常会呈现出速度异常体特征。地震走时层析成像在水文地质、工程与环境地质等方面也具有广泛的应用实例，例如水库大坝的防灾调查、建筑与工程设施的地基勘查和地下水污染调查等(Zelt et al.，2006；黄光南等，2013)。近年来VSP技术得到了较快发展，常规地面地震方法在VSP领域取得拓展性应用，基于VSP的速度和Q值层析成像也得到了一定程度的发展。但是大多数VSP地层速度反演是采用零偏观测方式情形。李文杰等(2004)利用零偏VSP资料先计算地层的反射系数，然后利用地层反射系数与速度的关系式求取层速度，数值模拟实验表明该方法具有较高的分辨率。刘冰等(2007)利用质心频移法对VSP资料进行地层Q值反演，数值实例表明这种方法对厚层反演效果较为准确，但是对薄层反演效果不稳定。乔玉雷(2010)将质心频移法运用于零偏VSP资料的地层Q值反演，胜利油田垦71井的地层Q值反演结果表明观测走时和地层速度对Q值反演效果具有较大影响。零偏VSP速度和Q值反演具有较高的精度，但是它的观测角度范围非常小。非零偏VSP能够扩大射线照射角度，因此可以反演出井周围的地层速度。王仰华(1988)结合非零偏VSP地震初至和广义线性反演算法反演地层速度，数值模拟和实际数据均取得了非常好的效果。Huang等(2012)将变阻尼约束算法加入层析成像方程，并对非零偏VSP资料进行地震走时层析成像，结果表明该算法可以克服射线不均匀覆盖带来的影响，从而提高速度反演结果的质量。周珺等(2012)利用逐层递推速度反演算法对非零偏VSP资料进行速度反演，数值模拟结果证明直线法适用于小井源距，而折线法适用于不同井源距情形。

上述研究表明地震走时层析成像在天然地震学和地球动力学等重大学科领域，以及石油、煤炭等资源勘探领域均具有实际应用价值。但是VSP地震走时层析成像大部分是针对一维和二维观测系统情形，三维VSP地震走时层析成像的研究非常缺乏，本文研究三维VSP地震初至波走时层析成像，旨在为三维VSP资料的速度反演结果提供指导作用。本文包括以下几部分研究内容：首先，介绍有限差分解“程函方程”算法计算速度模型的走时分布，以及射线路径的求取方法。然后，引出地震走时层析成像的基本原理和实现基本流程。最后，利用地震初至波走时层析成像对三维VSP速度模型进行速度反演，并分析反演结果与真实模型的差异，总结三维VSP速度反演结果的特点为将来实际资料的应用提供指导。

1 方法原理

1.1 三维有限差分走时射线追踪算法

有限差分走时计算方法就是采用离散的矩形网格对程函方程进行差分近似，结合每个网格节点的已知速度值，求取网格模型走时分布的一种方法。程函方程的一般数学表达式如下

(1)

式中，T是网格模型的走时，S是对应网格模型的慢度值，它是速度的倒数，X是表示网格模型的空间坐标向量，三维网格模型X=(x,y,z)，Xs是表示震源点坐标的向量。

三维模型的程函方程表达式如下

(2)

式中，t为走时；x,y和z分别为直角坐标轴；s为慢度(速度的倒数)。

本文采用了有限差分走时计算方法，在局部的有限差分走时算计过程中考虑了首波、体波和散射波的存在，同时采用反向延拓算法，使得算法能够适应复杂的地层速度模型，也保留了有限差分方法快速和精确的优点。当计算完三维速度模型的走时之后，可以得到三维网格模型的走时分布，根据网格模型的三维走时分布可以求取走时梯度。

x轴方向的走时梯度计算表达式为

gradt[1]=(t(i+1,j,k)+t(i+1,j+1,k)+t(i+1,j,k+1)+t(i+1,j+1,k+1)-t(i,j,k)-t(i,j+1,k)-t(i,j,k+1)-t(i,j+1,k+1))/(4×h)

(3)

y轴方向的走时梯度计算表达式为

gradt[2]=(t(i,j+1,k)+t(i+1,j+1,k)+t(i,j+1,k+1)+t(i+1,j+1,k+1)-t(i,j,k)-t(i+1,j,k)-t(i,j,k+1)-t(i+1,j,k+1))/(4×h)

(4)

z轴方向的走时梯度计算表达式为

gradt[3]=(t(i,j,k+1)+t(i+1,j,k+1)+t(i,j+1,k+1)+t(i+1,j+1,k+1)-t(i,j,k)-t(i+1,j,k)-t(i,j+1,k)-t(i+1,j+1,k))/(4×h)

(5)

射线路径是沿走时梯度最快下降的方向，从检波点至炮点进行反向追踪计算得到。在整个网格模型的走时场，射线路径垂直于走时梯度的等值线。射线路径可以从模型范围内的任意一个检波点开始，利用相应炮点的走时分布进行反向追踪求得。每一炮求取射线路径的运算量线性正比于射线路径的条数与射线路径的长度。为了提高射线追踪的计算效率，可以根据走时互易原理，如果炮点数远大于检波点数，可以将震源点设置于在检波点处，在炮点位置进行接收走时信息(Hole,1992)。

1.2 地震走时层析成像方程

地震层析成像通常采用线性化和离散化方法，将速度模型划分成网格模型，并假设每个网格单元的慢度值为常数。地震走时层析成像的走时扰动可以用慢度扰动沿射线路径的线性积分来表示。整个反演问题可以用一个大型稀疏线性方程组来表示

Ax=b

(6)

通常地震层析成像的阻尼最小二乘表达式写成如下形式

(7)

通常利用估计误差的倒数对每个观测点的数据进行加权，并在正则化过程中加入平滑约束算子，这样可以得到新的反演表达式

(8)

(9)

式中，λ参数的作用是平衡数据拟合度与平滑约束大小。Lees等(1989)描述了平滑约束算子的基本概念和作用，它是最小化模型粗糙度的一种约束方法，它是在层析反演的迭代过程中压制异常值，以避免反演陷入不正确的奇异值解。以三维速度模型为计算实例，第个网格的平滑约束方程为

6xj-(xj+1+xj-1+xj+n+xj-n+xj-m+xj+m)=0

(10)

式中，xj为网格j的慢度扰动值，m为网格模型一个平面内的网格数，n为网格模型一行的网格数(Tryggvason et al.,2002)。在实际计算过程中，除了采用平滑约束算子外，还可以通过设置速度反演的上下限数值大小，对不合理的速度值进行剔除处理。Zhou等(1992)综合运用了平滑约束算子和速度上下限算法对速度差异很大的非均匀速度模型进行了非线性走时反演。

1.3 三维VSP地震走时层析成像算法流程

图1 三维VSP地震初至波走时层析成像算法流程图Fig.1 Flowchart of 3D VSP first arrival traveltime tomography

三维VSP地震初至波走时层析成像结合了三维VSP地震射线追踪算法、地震走时层析成像方程与数学反演算法。根据地震走时层析成像的基本原理，这里制定了三维VSP地震初至波走时层析成像算法流程(图1)。实现流程包括以下几个步骤：(1)从采集的地震记录拾取观测走时。(2)建立三维初始速度模型，该初始速度模型应尽可能地接近真实模型，以便快速反演得到正确速度模型。(3)网格化初始速度模型，并利用有限差分算法计算各个网格节点的走时，根据走时分布求取地震波的射线路径坐标和接收点走时。(4)计算观测走时与正演走时的走时残差和，判断其是否满足给定的精度要求，如果满足精度要求，则输出速度模型，如果不满足精度误差，则进入下一环节。(5)建立地震走时层析方程组，采用线性或者非线性反演算法求解方程组，由于LSQR算法节省内存和计算量，也更适合求解不适定性问题，所以本文采用LSQR求解地震层析成像方程。(6)将求得的速度模型扰动量对初始速度模型进行更新得到新速度模型。(7)将其作为新的速度模型进行新一轮的正、反演与精度误差判断，经过一定的迭代次数之后，可以得到满足给定精度误差要求的速度模型，它代表了地下地层的真实速度分布情况。

2 数值模拟

为了得到VSP观测方式下三维地质模型的速度反演结果，并对速度反演结果进行分析。本文建立三维检测板速度模型，三维速度模型的大小为204 m×204 m×104 m。速度异常体的大小为30 m×30 m×15 m，低速异常体的速度值为1 800 m/s，高速异常体的速度值为2 200 m/s。采用VSP观测方式，炮井在地表的投影坐标为(105 m, 90 m, 0 m)，炮点数为51，第一个炮点的位于深度0m处，炮点在井中等间距(ΔS=2 m)分布。检波器分布于模型表面，在x与y轴方向均为等间距(ΔR=4 m)分布，共有2500个检波器。反演模型在x、y和z轴方向的网格间距均为Δ=2 m，三维速度模型的初始速度为2 000 m/s。本文首先对检测板速度模型进行正演运算得到理论走时(图2)。为了检验三维初至走时反演算法的稳定性和抗干扰能力，然后将理论走时加入了20%的高斯噪声(Zhou et al.,2008)。初至走时个数共有51×2 500个，为方便和清晰地对比理论初至与含高斯噪声初至之间的差别，因此图2仅仅展示了第51炮的一条测线的初至走时。

本文对理论初至和含高斯噪声初至分别进行了速度反演。为了得到合理的速度反演结果，在反演迭代过程中加入了平滑约束因子、阻尼因子与阿尔法滤波器。所有的速度反演结果都经过了20次迭代反演过程。图3为检测板速度模型、理论初至反演结果和含高斯噪声初至反演结果，由于采用了三维VSP观测方式，可以看出速度反演结果类似于“倒锥形”形状，浅部的反演效果比深部的反演效果更好，这是因为浅部的射线交叉性比深部更良好。为了更加方便地观察检测板速度模型的反演效果，本文选取检测板速度模型与反演结果在x=105 m处的垂直切片和x=8 m处的水平切片分别进行比较。对比检测板速度模型和反演结果的垂直切片(图4)，可以看出初至走时层析成像取得了较好的反演效果，井附近异常体的位置在反演前后的对应性较好，理论初至反演结果与含高斯噪声初至的反演结果没有明显区别。由于采用VSP观测方式，射线在浅部地层的交叉性较好，在深部地层趋向于单一射线方向，所以浅部地层的速度反演效果相比深部地层的速度反演效果更好。由于受地表检波点分布的约束，浅部异常体的速度反演结果“拉伸”变形作用较弱，反演结果的数值与真实结果相接近。然而深部异常体的速度反演结果“拉伸”变形作用较强，其反演结果的数值与真实结果相差较大。对比检测板速度模型和反演结果的水平切片(图5)，可以看出反演得到的异常体位置与检测板速度模型的异常体位置相对应。相比含高斯噪声初至反演得到的异常体形状，理论初至反演得到的异常体形状更加接近检测板速度模型的异常体形状(圆形)，这是因为大量高斯噪声使反演结果在一定程度上产生了畸变，如果继续提高噪声的百分含量会进一步使得异常体产生畸变。

图2 正演初至与添加20%高斯噪声的初至Fig.2 First arrival of forward modeling and 20% Gaussian noise added first arrivals

图3 三维VSP初至走时层析成像结果Fig.3 3D VSP tomogram by use of first arrival traveltimesa.三维检测板模型； b.理论初至反演三维 结果；c.含噪声初至反演三维结果

图4 三维检测板模型和速度反演结果的垂直切片Fig.4 Vertical slices of the 3D checkerboard model and the inverted velocity modelsa.检测板模型的垂直切片；b.理论初至反演结果的垂直切片；c.含噪声初至反演结果的垂直切片

3 结论

三维VSP勘探技术不仅可以获取井壁的地质信息，而且还可以获取井附近一定范围内的地质信息。三维VSP地震初至波走时层析成像作为一项VSP成像技术，它可以获取地层附近的速度信息，这对VSP偏移成像和地质建模具有重要作用。本文对数值实例进行研究和评价，这对VSP实际资料的速度反演具有指导意义。本文对三维VSP初至波走时层析成像的反演结果做出总结：

(1)受VSP观测方式的约束，三维VSP速度反演结果的形状类似于“倒圆锥形”。

(2)射线的交叉现象也更有利于提高反演质量，所以浅部地层比深部地层的反演效果较好。

图5 三维检测板模型和速度反演结果的水平切片Fig.5 Horizontal slices of the 3D checkerboard model and the inverted velocity modelsa.检测板模型的水平切片；b.理论初至反演结果的水平切片；c.含噪声初至反演结果的水平切片

(3)随着深度的增加，射线方向趋于单一射线方向——趋近于井深方向，射线的交叉现象越少，所以井的底部速度反演效果较差。

(4)通过对含高斯噪声的初至进行反演可知该算法具有较好的稳定性和抗干扰能力。根据以上几点结论，如果需要改善VSP初至波走时层析成像质量，就需要尽可能地扩大走时观测范围，使射线的交叉性程度得到进一步增强。

致谢 本文的研究工作得到了瑞典乌普沙拉大学地球科学系Ari Tryggvason研究员和斯伦贝谢公司Artem Kashubin博士的指导，在此感谢他们的热情帮助。

常旭, 刘伊克, 王志君, 等. 1996. 用波速层析成像方法提高储层横向预测精度 [J]. 地球物理学报, 6(6): 813-822.

和锐, 杨建思, 张翼. 2007. 地震层析成像方法综述 [J]. CT理论与应用研究, 16(1): 35-48.

黄光南, 刘洋, Tryggvason A, 等. 2013. 变网格间距速度层析成像方法 [J]. 石油地球物理勘探, 48(3): 379-389.

李录明, 贺玉山, 罗省贤. 2013. 高效高精度初至波层析静校正方法及应用 [J]. 成都理工大学学报:自然科学版, 40(2): 113-119.

李文杰, 魏修成, 刘洋. 2004. 利用VSP资料反演地层层速度的一种新途径 [J]. 石油物探, 43(2): 126-129.

刘冰, 王德利, 巩向博,等. 2007. 基于VSP资料进行的吸收系数Q值反演 [J]. 吉林大学学报, 37(增刊): 49-52.

乔玉雷. 2010. 零偏VSP资料质心频移法在胜利油田Q值反演中的应用研究 [J]. 山东科技大学学报, 29(5): 8-12.

秦宁, 李振春, 杨晓东. 2011. 基于角道集的井约束层析速度反演 [J]. 石油地球物理勘探, 46(5): 725-731.

孙章庆, 孙建国, 韩复兴. 2012. 针对复杂地形的三种地震波走时算法及对比 [J]. 地球物理学报, 55(2): 560-568.

王仰华. 1988. 根据VSP初至走时反演地层速度 [J]. 石油物探, 27(3): 16-25.

徐涛, 张明辉, 田小波, 等. 2014. 丽江-清镇剖面上地壳速度结构及其与鲁甸Ms6.5级地震孕震环境的关系 [J]. 地球物理学报, 57(9): 3069-3079.

周珺, 谢春辉, 杨鹏. 2012. VSP初至逐层递推反演层速度[J]. 物探与化探, 36(2): 242-245.

Aki K，Lee W H K. 1976. Determination of three dimensional velocity anomalies under a seismic array using first P arrival times from local earthquakes, 1. a homogeneous initial model [J]. Journal of Geophysical Research, 81(23):4381-4399.

Calvert A J, Fisher M A. 2001. Imaging the Seattle Fault Zone with high-resolution seismic tomography [J]. Geophysical Research Letters, 28(12): 2337-2340.

Hole J A. 1992. Nonlinear High-Resolution Three-Dimensional Seismic Travel Time Tomography [J]. Journal of Geophysical Research, 97(B5): 6553-6562.

Huang G N，Liu Y. 2012. Variable damping constraint tomography and its application in VSP data [J]. Applied Geophysics, 9(2): 177-185.

Lees J M. 1992. The magma system of Mount St. Helens non-linear high-resolution P-wave tomography [J]. Journal of Volcanology and Geothermal Research, 53(1): 103-116.

Lees J M，Crosson R S. 1989. Tomographic inversion for three-dimensional velocity structure at Mount St. Helens using earthquake data [J]. Journal of Geophysical Research, 94(B5): 5716-5728.

Phillips W S，Fehler M C. 1991. Traveltime tomography: A comparison of popular methods [J]. Geophysics, 56(10):1639-1649.

Rawlinson N，Fishwick S. 2012. Seismic structure of the southeast Australian lithosphere from surface and body wave tomography [J]. Tectonophysics, 572(1): 111-122.

Tryggvason A, Rognvaldsson S T，Flovenz O G. 2002. Three dimensional imaging of the P- and S-wave velocity structure and earthquake locations beneath southwest Iceland [J]. Geophysical Journal International, 151(1): 848-866.

Zelt C A, Azaria A, Levander A. 2006. 3D seismic refraction traveltime tomography at a groundwater contamination site [J]. Geophysics, 71(5): H67-H78.

Zhang J，Toksoz M N. 1998. Nonlinear refraction traveltime tomography [J]. Geophysics, 63(5): 1726-1737.

Zhao L F, Xie X B, He J K， et al. 2013. Crustal flow pattern beneath the Tibetan Plateau constrained by regional Lg-wave Q tomography [J]. Earth and Planetary Science Letters, 383(1): 113-122.

Zhou B, Sinadinovski C, Greenhalgh S A. 1992. Non-linear inversion travel-time tomography: imaging high-contrast inhomogeneities [J]. Explorational Geophysics, 23(1): 459-464.

Zhou B, Greenhalgh S, Green A. 2008. Nonlinear traveltime inversion scheme for crosshole seismic tomography in tilted transversely isotropic media [J]. Geophysics, 73(4): D17-D33.

Zhou H W. 2006. Multiscale deformable-layer tomography [J]. Geophysics, 71(3): R11-R19.

3D VSP Seismic First Arrival Traveltime Tomography

HUANG Guang-nan1,2,3, DENG Ju-zhi1, LI Hong-xing1, LI Ze-lin1, ZHANG Hua1, WANG An-dong1

(1. Fundamental Science on Radioactive Geology and Exploration Technology Laboratory, East China University of Technology, Nanchang，JX 330013, China; 2. Hubei Subsurface Multi-scale Imaging Key Laboratory (SMIL), China University of Geosciences, Wuhan, HB 430074, China; 3. State Key Laboratory of Petroleum Resources and Prospecting, China University of Petroleum, Beijing 102249, China)

VSP (Vertical seismic profile) is a high resolution exploration technology. It plays a significant role in the exploration and development stage of oil and gas. 3D VSP first arrival traveltime tomography can invert the underground velocity model directly. This method combines 3D VSP seismic ray tracing algorithm and mathematical inversion algorithm. It can reconstruct the underground velocity models directly. These models can provide apparent evidence to know the underground geological structure better. In order to test the stability and anti-noise ability of the algorithm, there are 20% Gaussian noises were added into the 3D VSP forward traveltimes. Then, the 3D VSP first arrival traveltime tomography is used to invert the velocity models. We have obtained 3D velocity models for the checkerboard model, and its vertical slices and horizontal slices. The more rays cross with each other, the better the reliability and effect of inversion result are. Therefore, the effect of shallow strata inversion result is better than the deep strata, because there are more rays cross with other in shallow strata than the deep strata. On the contrary, the velocity effect in deep strata is much worse. The inversion algorithm has inverted good results by use of the noise contained traveltimes. This confirms that the robust algorithm has strong ability to anti-noise.

vertical seismic profile; Gaussian noise; ray tracing; seismic tomography

2015-09-16

国家自然科学基金(41504095, 41104073, 41004048, 41364004, 41104074, 41304097); 国家科技支撑计划(2011BAB04B03);国家科技重大专项(2011ZX05024-001-02); 中国地质大学(武汉)地球内部多尺度成像湖北省重点实验室基金(SMIL-2015-10);江西省教育厅基金(GJJ14476); 东华理工大学博士科研启动基金(DHBK2013212)

黄光南(1983—)，男，讲师，主要从事地震层析成像方法研究。E-mail: bobking2@126.com

10.3969/j.issn.1674-3504.2016.03.010

P631.4

A

1674-3504(2016)03-0266-07