复杂近地表速度广义近似反演方法研究

李 哲,程丹丹

(中国石油大学(华东)，山东青岛266580)

随着我国油气勘探的重点目标由东部地区向西部和南方探区转移，复杂近地表成为制约地震勘探效果的关键因素。在以胜利油田和大庆油田为主体的东部探区，近地表相对比较简单，具有典型的低速带、降速带及正常速度带3层结构特征。此类近地表结构中地震波速度向下逐层递增，一般采用折射波法[1]来获取近地表结构。在复杂近地表探区，地形起伏剧烈，低速带与降速带不明显，同时表层横向岩性与速度变化快，常用的近地表调查和速度建模方法效果较差，不能满足地震勘探的要求。

地震波旅行时层析是一种利用地震波在介质中传播时间来反演地震波所经过路径上介质速度的方法。该方法对地下介质的速度分层结构没有特殊要求，并能适用于横向速度变化剧烈的地下介质，因而成为最有效的近地表速度反演方法。国内外众多学者对旅行时层析反演方法进行了研究。自1984年第54届SEG年会上设置了地震层析成像研究的专题之后，地震层析成像研究在地震勘探领域得以快速推进。经过近30年的发展，国内外关于地震层析成像方法的研究取得了丰硕的成果，主要有透射波层析[2-4]、反射波层析[5-8]、折射波层析[9-10]、井间层析[11-13]、多波联合层析[14-16]、初至波层析[17-19]、波形层析[20-21]、起伏地表层析[22-23]等。这些工作也极大地推动了我国地震层析成像方法研究和实际应用的步伐。

对于复杂近地表探区，旅行时速度反演在地震勘探中主要有3个用途：一是静校正量的计算，只要静校正量合适，不要求速度值与地层厚度正确，相当于求解半确定性问题;二是用于观测系统参数的确定，要求大尺度范围内一致，主要涉及震源激发深度的确定;三是偏移速度模型的建立，要求小尺度上的一致性。本文针对复杂近地表旅行时层析成像的数据特征，设计和推导了基于广义逆的层析反演近似算法，并提出了4种针对性的约束方法。通过约束广义近似旅行时层析速度反演建立合适的近地表速度模型，用于复杂近地表区地震资料的层析静校正，恢复了炮集记录上反射波同相轴的双曲线形态。

1 复杂近地表探区的特点

1.1 复杂近地表探区的地质特点

复杂近地表是指地形起伏大、速度变化迅速、岩性变化剧烈、厚度横向差异大的浅表结构。其主要地质特点体现在以下方面。

1) 地表高程差别大。例如南方海相碳酸盐岩探区，一般为大山区，地表高差往往在四五百米以上，鄂西、镇巴和黑池梁等探区局部高程差可达1000多米[24];在西部黄土塬探区，地形起伏高差也有100多米。

2) 表层速度变化迅速。有的地表速度只有几百米每秒，而不远处就能达到2000m/s以上;在一些灰岩出露地区，已经不存在传统意义上的低、降速带，表层速度已经达到3000～4000m/s。

3) 表层岩性变化剧烈。复杂近地表受到了大自然环境变化的影响，是大气圈渗入岩石圈而形成的特殊地质结构，其岩性特征、层理结构、胶结程度等受大自然环境作用而发生了改变，横向均质性极差。例如我国西部典型的山地、沙漠、戈壁、砾石覆盖区、黄土覆盖区等探区，横向上岩性极不稳定。正是这种横向上岩性的变化，结合差异风化作用，使得低、降速带厚度横向差异明显。

1.2 复杂近地表条件下地震资料特点

复杂近地表条件严重制约地震资料质量，一方面是通过震源激发特性影响地震资料的幅频特性和信噪比;另一方面是接收到的地震波先后两次经过复杂近地表传播，其运动学与动力学特性受到明显影响。由于复杂近地表介质的非完全弹性性质与横向非均质性，地震波发生吸收、衰减、散射和漫射等现象，不仅使地震波场变得复杂，也严重衰减有效波地震能量，降低了信噪比。

复杂近地表条件下，在共炮点道集上直达波同相轴不再是直线，反射波时距曲线也明显偏离双曲线特征(图1)。这类地震记录的初至波拾取存在2个问题：①存在极低信噪比的地震道，初至波无法识别与拾取，减少了可用于旅行时层析的数据;②由于地震资料信噪比低，以及波形发生畸变等原因，拾取到的初至时间存在较大的误差，这将严重影响表层旅行时层析速度反演的精度。因此，研究设计层析反演方法要考虑实际旅行时数据误差。

图1 起伏地表区典型炮集记录

为此，在复杂近地表旅行时层析反演中，需要通过旅行时数据的统计分析确保用于速度反演的旅行时数据正确，并充分利用已知信息来约束速度反演，以弥补旅行时数据的不足。

2 复杂近地表旅行时层析广义近似方法

旅行时层析成像的基本数学原理是Radon变换[25]。在旅行时层析中将地下介质按一定大小进行网格化，用vj(j=1,2,…,M)表示第j个网格的速度。根据地震波传播时间与速度、地震波传播路径的关系，得到的表达式为

(1)

式中：ti为地震波沿第i(i=1,2，…，N)条射线传播的时间;lij为第i条射线在第j个网格中的长度。公式(1)中速度vj与旅行时ti之间的关系不是线性关系，用速度vj的倒数sj来替换公式(1)中的1/vj，从而得到用慢度表示的方程：

(2)

公式(2)形式上旅行时ti与未知量慢度sj已经成为线性关系。但旅行时层析反演依然不是一个简单的线性反演，因为这个反演中地震波传播的路径也不是已知的。实际上，这个反演不仅要反演出速度分布，而且必须同时反演出射线分布。由于地震波传播路径由速度分布与费马原理共同决定，旅行时层析速度反演实际上是一个复杂的非线性反演问题。

2.1 广义反演及其近似表达式

用列向量d=(t1,t2，…,tN)T表示有效的旅行时数据，列向量s=(s1,s2，…，sM)T表示需要求取的慢度;将每条射线在网格中的长度用下面的射线矩阵G来表示

(3)

则公式(2)变为

(4)

射线矩阵G不是方阵，公式(4)不能直接求解。考虑到射线矩阵G与s的依赖关系，假设初始慢度s0，用于获取射线矩阵G，并计算出地震波在慢度s0下传播的时间d0，即

(5)

进而得到

(6)

式中：Δs表示s与s0的差值；Δd表示d与d0的差值。将公式(5)代入公式(6)得到

(7)

求解公式(7)得到修正的慢度。由于射线参数矩阵G是根据初始慢度s0计算的，与实际射线矩阵有差异，因而求取的慢度和实际慢度有差异。为此，需要用修正后的慢度作为初始慢度，重新进行计算。这样就形成了旅行时层析速度反演的迭代算法。

由于射线矩阵G不是正定矩阵，公式(7)不能直接求解。可以利用最小二乘法[27]得到慢度修正量Δs的广义反演形式：

(8)

公式(8)的求解需要计算M×M阶的逆矩阵，M为慢度修正量的未知数数目，通常比较大。公式(8)难以直接计算，同时参数矩阵G与实际射线矩阵有差异，公式(8)精确计算的结果肯定是错误的，为此考虑近似求解。

用J表示正定矩阵GTG，其元素Jkm是经过第k网格和第m网格的射线长度数据的互相关，即

(9)

公式(9)右边被求和项只有当一条射线同时通过第k网格和第m网格才不为零，其它项均为零。即：当k≠m时，Jkm接近于零;当k=m时，被求和项为网格的自相关，当一个网格没有射线经过时，Jkk=0。剔除没有射线经过的网格，则可以得到一个主对角线元素绝对占优的对角阵J。因此反演中只考虑J的主对角线元素，得到公式(8)的近似公式

(10)

公式(10)不需要求解大型矩阵，在近似计算中只需计算主对角线元素和射线矩阵的乘法，因此算法是稳定的。

2.2 复杂近地表初始模型建立

建立初始模型是旅行时层析反演的关键。常用的建模方法有2种。

1) 建立与已知信息相一致的初始模型。在地震勘探中往往利用微测井、小折射和岩性探测技术等获取探区内某些关键点上的实测表层速度与岩性信息。利用这些已知点的地层速度分布，结合地形起伏情况，可以建立起一个比较接近实际的速度模型。基于该初始模型，只需要反演出较小的慢度修正量。

2) 按旅行时层析反演算法的要求来建立初始模型。旅行时反演中只有那些有射线经过的网格才能被反演出来。对于近地表来说，如果一开始建立的模型浅层速度比深层速度大，那么地震波射线就只会在浅表，下部介质里没有射线，自然就反演不出来。因此，应建立速度随深度逐渐增加的模型，让射线能够下传到下部介质，符合大多数情况，即表层速度基本结构是速度随深度增加的。

无论采用哪种建模方法，都应在反演中利用微测井、小折射和岩性探测技术等提供的实测速度信息作为反演的约束条件，以减少反演的多解性。

2.3 旅行时拾取误差的影响

第一个地震波到达接收点，接收点开始振动的时刻为起跳时间，这个时间就是旅行时反演所需要的初至时间。由于存在背景噪声，当地震波引起的振动小于背景噪声时，无法确定哪个时刻是起跳时间。通常在技术上总是拾取地震记录上第一个能量明显高于背景振动的到达时间，作为初至波时间。

在复杂近地表条件下，地震记录变得十分复杂，无论是直达波还是反射波，其同相轴均偏离了正常的形态且伴随着低信噪比，进一步增加了旅行时拾取的难度和可信度。

在近地表旅行时层析反演中炮检距相对较小，因而旅行时拾取的相关误差很大。由于讨论相对误差对速度反演的影响比较困难，假设所有旅行时均具有相同的绝对误差σ，则根据公式(10)进行反演，估算慢度误差列向量为

(11)

公式(11)是旅行时数据误差传递方程。假设射线在每个网格中的长度相同，则速度反演结果将被放大N/L倍(L为射线的平均长度)。因此，反演中应该考虑旅行时数据的误差，避免将旅行时误差反演到速度模型中。

2.4 低、降速带的确定

旅行时层析反演获得近地表速度分布后，需要进一步确定低、降速带。根据公式(10)可知，没有射线经过的网格不参加反演，即这些网格的速度不正确。实际应用时首先根据射线分布确定速度反演的可靠范围，在该范围内分析低、降速带的分布。

利用微测井、小折射以及岩性探测等技术提供的局部近地表可靠速度分层信息，对近地表旅行时层析反演资料进行一致性检测，进而标定速度范围，同时结合层析射线的密度进行低、降速带划分。

3 实际地震资料试处理

对南方海相碳酸盐岩某复杂地表区实际二维地震资料进行旅行时层析广义近似反演试处理。该区地形起伏达到500m，岩石出露地表，最低速度在1400m/s左右，属于典型的复杂近地表结构。根据工区的基本情况，考虑到地表速度较高，建立如图2 所示的速度随深度逐渐增加的初始模型。

对该二维测线的所有炮集数据进行初至波拾取，并绘制初至波旅行时与炮检距的交会图，结果如图3a所示。观察图3a可知，炮检距在4000m以上时存在散点，即可能存在误差。近地表初始模型的深度为1000m，而大炮检距的地震波向下的穿透深度超过了1000m，为此切除了炮检距大于2500m的初至时间，结果如图3b所示。

图2 南方复杂地表区实际二维测线初始速度模型

图3 实际二维地震资料拾取的初至波旅行时随炮检距的分布a 所有初至波旅行时; b 小炮检距初至波旅行时

图4为直接利用图3b所示小炮检距初至波旅行时数据进行无约束广义近似层析反演给出的射线分布和速度分布，可以看出剖面两边由于边界原因，速度与初始模型接近，反演不充分;中部射线比较松散，浅表射线比较集中;从速度分布上看，没有清晰地反演出低、降速带。

图4 无约束广义近似层析反演的射线分布(a)与速度分布(b)

通过对实际资料的分析，根据工区特点引入4种约束条件：①速度区间的控制，根据已知的先验信息设定一个速度区间，防止速度过小和过大;②利用小折射和微测井资料提供的先验信息进行约束;③速度奇异点约束，利用中值滤波去除速度奇异值;④利用附近网格修正量均值修正无射线通过网格。图5为仍然采用图3b中的初至波旅行时数据进行增加约束后的广义近似层析反演给出的射线分布与速度模型。从射线分布(图5a)上看，对低、降速带有较好的覆盖，近地表速度模型(图5b)的层状结构非常明显。

利用约束广义近似层析反演出的合理的近地表速度模型(图5)，以3000m/s为低降速带门槛值，得到低、降速带底界，进而计算静校正量，对实际二维地震资料的炮集记录进行静校正处理。图6a 为原始炮集记录，直达波和反射波同相轴畸变明显;图6b是应用静校正量后的结果，反射波同相轴基本上恢复了双曲线的特性。

图5 约束广义近似层析反演的射线分布(a)与速度分布(b)

图6 复杂地表区实际二维地震资料单炮记录静校正前(a)、后(b)对比

4 结束语

复杂近地表条件下地震波场复杂，地震资料信噪比低，拾取的旅行时误差比较大。我们提出的广义近似层析反演方法是一种考虑旅行时拾取误差的层析成像方法，能够获得较稳定的近地表速度模型。针对旅行时误差与速度模型误差的关系，推导了慢度反演误差的计算公式，可进一步研究将其应用于反演精度控制的方法。

针对复杂地表区实际地震资料的无约束广义近似层析反演结果存在对近地表的分层不明显、射线完全集中在浅表等问题，研究并提出了4种约束方法，使得约束后反演的射线能较好地覆盖低降速带，速度分布清晰地展示了近地表地层结构。将约束广义近似层析反演出的速度模型用于复杂地表区实际地震资料的静校正，恢复了反射波同相轴的双曲线形态。

致谢:中国石油大学(华东)地球科学与技术学院地球物理系莫言工作组为本方法的研究提供了所有资料，课题组老师和同学对本文的研究提出了宝贵意见，在此表示衷心感谢！

参 考 文 献

[1] 徐凌,崔兴福,齐莉,等.复杂近地表层析反演与波场延拓联合基准面校正[J].石油物探,2007,46(3)：226-230

Xu L,Cui X F,Qi L,et al.Combined datum correction of tomography inversion with wavefield continuation under complex near surface condition[J].Geophysical Prospecting for Petroleum,2007,46(3):226-230

[2] Bording R P,Gersztenkorn A,Lines L R,et al.Applications of seismic travel-time tomography[J].Geophysical Journal of the Royal Astronomical Society,1987,90(2):285-303

[3] Gelbke C,Miranda F,Sattel G.Results of a seismic transmission tomography survey at the Grimsel Rock Laboratory[J].Log Analysist,1989,30(4):243-260

[4] 张建中,丁兴号.一种2.5维井间透射波层析成像方法[J].声学学报,2007,32(1):91-96

Zhang J Z,Ding X H.An approach to 2.5-D crosswell seismic transmission tomography[J].Acta Acustica,2007,32(1):91-96

[5] Zhou C,Crawley S,Whitmore D,et al.Seismic reflection tomography with 3D RTM angle gathers[EB/OL].(2013-06-10)[2013-12-30].http:∥earthdoc.eage.org/publication/publicationdetails/?publication=69508

[6] Stork C.Reflection tomography in the postmigrated domain[J].Geophysics,1992,57(5):680-692

[7] 井西利,杨长春,王世清.一种改进的地震反射层析成像方法[J].地球物理学报,2007,50(6):1831-1836

Jing X L,Yang C C,Wang S Q.A improved seismic reflection tomographic method[J].Chinese Journal of Geophysics,2007,50(6):1831-1836

[8] Brzostowski M A,McMechan G A.3-D tomographic imaging of near-surface seismic velocity and attenuation[J].Geophysics,1992,57(3):396-403

[9] Mari J L O,Mendes M.3D Near surface imaging thanks to refraction tomography and delay method[EB/OL].(2010-06-14)[2013-12-30].http:∥www.earthdoc.org/publication/publicationdetails/?publication=39624

[10] Afanasiev M,Pratt R G,Kamei R,et al.Crosshole Tomography,Waveform Inversion and Anisotropy-A Combined Approach Making Use of Simulated Annealing[EB/OL].(2012-06-04)[2013-12-30].http:∥www.earthdoc.org/publication/publicationdetails/?publication=59614

[11] 左建军,林松辉,孔庆丰,等.井间地震直达波和反射波联合层析成像及应用[J].石油地球物理勘探,2011,46(2):226-231

Zuo J J,Lin S H,Kong Q F,et al.Crosswell seismic direct wave and reflected wave joint tomographic imaging and Application[J].Oil Geophysical Prospecting,2011,46(2):226-231

[12] 瞿辰,杨文采,于常青.高分辨率过井地震层析成像的级联算法及其应用[J].物探与化探,2010,34(6):814-820

Huo C,Yang W C,Yu C Q.High resolution over the cascade algorithm well seismic tomography and its application[J].Geophysical and Geochemical Exploration,2010,34(6):814-820

[13] 李培明,闫智慧,郭明杰,等.联合应用初至和浅层反射的二维可形变层析静校正[J].石油地球物理勘探,2010,45(5):647-654

Li P M,Yan Z H,Guo M J.The combined application of first break and shallow reflection of 2D deformable tomographic static correction[J].Oil Geophysical Prospecting,2010,45(5):647-654

[14] 张树林,朱介寿,贺振华.井间地震和逆VSP 联合层析成像[J].石油地球物理勘探,1993,28(5):577-583

Zhang S L,Zhu J S,He Z H.Joint tomography of interborehole seismic data and inverse VSP data[J].Oil Geophysical Prospecting,1993,28(5):577-583

[15] Meléndez A,Korenaga J,Sallarǒs V,et al.TOMO3D-A new 3-d joint refraction and reflection travel-time tomography code for active-source seismic data[EB/OL].(2013-06-04)[2013-12-30].http:∥www.researchgate.net/publication/258621526_tomo3d_a_new_3-D_joint_refraction_and_reflection_travel-time_tomography_code_for_active-source_seismic_data

[16] Noble M,Thierry P,Taillandier C,et al.High-performance 3D first-arrival traveltime tomography[J].The Leading Edge,2010,29(1):86-93

[17] Li S,Fomel S,Vladimirsky A.Prestack first-break traveltime tomography using the double-square-root eikonal equation[J].Expanded Abstracts of 82ndAnnual Internat SEG Mtg,2012,1-6

[18] 袁刚,冯心远,蒋波,等.约束层析反演及其在地震速度计算中的应用[J].石油物探,2013,52(1):55-59

Yuan G,Feng X Y,Jiang B,et al.Constrained tomography inversion and its application in seismic velocity computation[J].Geophysical Prospecting for Petroleum,2013,52(1):55-59

[19] Smithyman B,Pratt R G,Hayles J,et al.Detecting near-surface objects with seismic waveform tomography[J].Geophysics,2009,74(6):WCC119-WCC127

[20] Sheng J,Leeds A,Buddensiek M,et al.Early arrival waveform tomography on near-surface refraction data[J].Geophysics,2006,71(4):U47-U57

[21] Zhang K,Li Z C,Zeng T S,et al.Tomographic velocity inversion for ADCIGs in areas with a rugged surface[J].Applied Geophysics,2012,9(3):313-318

[22] 张兵,徐兆涛,王华忠,等.山前带地震数据共成像点道集层析速度反演建模方法研究[J].石油物探,2012,51(6):590-597

Zhang B,Xu Z T,Wang H Z,et al.Common imaging gathers tomography velocity inversion and model building in foothill area[J].Geophysical Prospecting for Petroleum,2012,51(6):590-597

[23] 董良国,郭晓玲,吴晓丰,等.起伏地表弹性波传播有限差分法数值模拟[J].天然气工业,2007,27(10):38-41

Dong L G,Guo X L,Wu X F,et al.Finite difference numerical simulation for the elastic wave propagation in rugged topography[J].Natural Gas Industry,2007,27(10):38-41

[24] Helgason S.The Radon transform[M].2nd ed.USA:Birkhauser Boston Inc,1999:1-200

[26] Fitzgibbon A,Pilu M,Fisher R B.Direct least square fitting of ellipses[J].Pattern Analysis and Machine Intelligence,IEEE Transactions on,1999,21(5):476-480