地球物理学报  2019, Vol. 62 Issue (2): 619-633   PDF    
初至波相位走时层析
刘玉柱, 伍正, 耿志成     
同济大学海洋地质国家重点实验室, 上海 200092
摘要:近地表速度结构通常是利用射线走时层析或菲涅尔体走时层析等反演方法得到的,但它们的目标函数仍利用射线走时残差构建,导致反演精度不高.为此,本文提出了基于散射积分算法的初至波相位走时层析成像方法.该方法的核心是:①提出了依赖于频率的相位走时概念;②利用依赖于频率的相位走时信息,而非单一的无限频率射线走时;③发展了一种改进的相位展开方法,即通过监测相位不连续性和2π周期判定来消除相位折叠现象;④考虑了地震波传播的有限频特征,即基于波动理论而非传统的射线路径或有限空间的菲涅尔体构建核函数.通过利用Overthrust模型的数值实验及与传统射线走时层析和菲涅尔体走时层析的对比表明:本文提出的方法是一种有效的初至波走时反演方法.同时,基于Overthrust模型的数值试验还证明了下列结论,即通过挖掘更多的走时信息的确可以获得更高的反演精度和分辨率.
关键词: 地震层析成像      初至波      相位走时      散射积分算法      核函数     
First-arrival phase-traveltime tomography
LIU YuZhu, WU Zheng, GENG ZhiCheng     
State Key Laboratory of Marine Geology, Tongji University, Shanghai 200092, China
Abstract: Near-surface velocity structures are generally built by the ray-based traveltime tomography or the Fresnel-volume tomography in which, however, the objective function is usually constructed with the traveltime residual of raypaths. As a result, the inversion accuracy is not high. To obtain results with higher accuracy and resolution, we propose a first-arrival phase-traveltime tomographic method with an improved scattering-integral algorithm. The core of this method is:①putting forward the concept of frequency-dependent phase traveltime; ②inverting the frequency-dependent traveltime instead of the asymptotic ray-based traveltime; ③developing an improved phase-unwrapping method, e.g. unwrapping the phase by means of monitoring the phase discontinuity and judging the 2π jumps; ④taking the finite-frequency effect of the seismic wave propagation into account, using the wave equation but not raypaths or Fresnel volumes as the sensitivity kernel. The numerical experiments on an Overthrust model and the comparisons between this method with the conventional ray-based traveltime tomography and Fresnel-volume traveltime tomography prove that this method is an effective first-arrival traveltime inversion method. Meanwhile, the numerical experiments also prove that the inversion accuracy and resolution can be apparently improved by digging more traveltime information contained in the first-arrival waveform.
Keywords: Seismic tomography    First arrival    Phase-traveltime    Scattering-integral algorithm    Sensitivity kernel    
0 引言

准确估算表层速度模型,借此消除复杂表层因素在地震资料上产生的静校正问题,是野外勘探和室内资料处理的共同任务.通常,近地表速度结构是通过射线走时层析得到的(Bishop et al., 1985; Zhu et al., 1992),由此会带来一些问题.事实上,地震波是有限频带时域信号,不仅几何射线路径上的介质性质会影响地震波的传播,而且其周围介质的物理性质同样会影响地震波的走时.此外,波的散射、衍射、波前分裂和波前弥合(Liu and Dong, 2012)等都将导致地震波的走时不仅仅依赖于几何射线路径上的介质性质.因此,如果直接利用地震波的高频渐近理论,即几何光学近似理论来实施速度反演,会使最终结果或多或少地偏离真实的速度分布,尤其是难以分辨出小尺度不均匀体,即难以考虑尺度小于优势波长的不均匀体对速度不均匀分布的贡献.

相对于传统的基于几何射线理论的层析成像方法,有限频层析以“波路径”替代传统走时层析成像中的“几何射线路径”,更准确地描述了地震信号对介质速度结构的灵敏度.此外,“波路径”对低速区域具有同高速区域一样的敏感性,从而改善了高频射线路径对高、低速区域采样的差异所产生的问题(Wielandt, 1987).Spetzler和Snieder (2004)刘玉柱等(2009)基于有限频理论提出了菲涅尔体地震层析成像方法,进一步改善了射线走时层析的精度.然而,该方法只在有限的空间区域被认为对观测数据是敏感的,反演精度仍有待提高.Sheng等(2006)提出了一种反演近地表速度结构的初至波波形层析成像方法.相对于传统的射线走时层析,初至波波形层析利用了更多的信息,因此能够得到更精确的结果,但它的计算成本与常规的全波形反演(Full Waveform Inversion, FWI)相同,且同样面临FWI中的诸多挑战.

除了几何射线路径和波路径以外,利用波场相位信息同样可以反演地下速度结构(Shin and Min, 2006; Bednar et al., 2007).这是因为在Rytov近似框架下介质扰动引起的波长残差仅仅来自于相位.然而,绝大多数相位反演的应用都仅仅提取了相位而没有消除相位折叠现象,这样会使反演陷入局部极值,严重破坏了反演的稳定性.Oppenheim和Schaefer (1975)首次提出了对波场取对数和对波场取导数两种提取相位的方法,但没有解决相位折叠问题.Poggiagliolmi等(1982)通过监测相位的不连续性进行相位展开,取得了一定成效,但对高频情况下低信噪比数据的应用很不理想.Bednar等(2007)将相位反演不稳定现象归因于进行相位展开操作时丢失了某些信息.Gélis等(2007)在基于Rytov波路径的弹性波反演实验中证明了相位展开的重要性,但却没有实际应用.Choi和Alkhalifah (2013)提出了瞬时走时的概念,其实质是对波场取导数求取频率依赖的走时信息,并将相位展开操作应用到实际反演中,得到了高分辨率的反演结果,然而其实验中频带上限只有5 Hz,对于更高的频率瞬时走时的求取会出现较大误差.

针对上述在相位反演中所存在的问题,本文对相位反演进行了改进并将其基本思想用于初至波走时层析中,提出了初至波相位走时层析成像方法.一方面,通过对波场取对数提取频率依赖的走时(为了与其他走时区分,本文称其为相位走时),丰富了对走时信息的利用;另一方面,通过对波场施加强阻尼得到初至波形,从而压制噪声,解决了Poggiagliolmi等(1982)相位展开方法中对高频情况下低信噪比数据应用不理想的问题,避免了相位提取中的周期跳跃问题.相比于瞬时走时,本文方法实验中频带上限也能做到更高.同时,该方法利用Rytov近似计算相位走时层析的全空间敏感核函数,通过使用Liu等(2015)提出的改进的散射积分算法,很好地解决大规模矩阵的存储问题,提高了反演的稳定性.本文最后将Overthrust模型的初至波相位走时层析反演结果与传统射线走时层析和菲涅尔体走时层析的反演结果进行对比.

1 理论方法

指数衰减波场u(t)可以表示为原波场u(t)与指数衰减项e-αt的乘积:

(1)

其中,α为阻尼因子.对上述阻尼波场作傅里叶变换得

(2)

其中,ω为角频率.在(2)式的最右项中,阻尼因子α位于复角频率的虚部,这意味着频率域正演时使用复角频率将会生成指数衰减波场,其中复角频率的虚部为阻尼因子.

当阻尼因子很大时,可以压制初至波之后的波形,从而得到初至波形(Shin et al., 2002),近似表达为:

(3)

其中,ξτ分别为初至波的振幅和走时.对(3)式作傅里叶变换可以写成

(4)

在(4)式中,相位为初至走时与复角频率的乘积.对该式左右取对数,将振幅与相位分离,再用相位项除以对应的角频率ω

(5)

便能够获得初至时刻τ.该方法计算得到的初至时刻是依赖于频率的,相当于提取了地震记录的相位走时信息.

本文将通过监测相位的不连续性和2π周期判定来展开相位.当某点的相位跳跃大于某个特定值时,该点即为一个不连续点,在该点将相位翻折上去,即进行了一个相位展开操作.对所有不连续点都进行该操作,即可得到真实的频率依赖的相位走时.

反演中使用的频率相关的相位走时层析核函数为(Woodward,1992)

(6)

式中,x表示空间任意一点,xr表示检波点位置,xs表示炮点位置,G0(ω, x, xr)或G0(ω, x, xs)为当前模型c0xrxsx点的格林函数,u0(ω, xr, xs)为炮点xs激发在x点的波场,u0(ω, xr, xs)为炮点xs激发在xr点的波场.

同时,本文使用Liu等(2015)提出的改进的散射积分算法,将梯度计算中大规模的核函数-向量乘运算表示为具有明确物理含义的向量-标量乘的累加运算,大大减少内存占用量的同时也能够方便地实现预条件最速下降方向(7),提高收敛效率.

(7)

式中, Δτ为提取的相位走时残差.

2 初至波相位信息的提取

本节展示初至波相位信息的提取过程,使用Overthrust模型进行数值实验.理论模型如图 1所示,模型包括801×187个网格,网格间距均为25 m,模型总大小为20 km×4.65 km.共设置199个炮点,均匀分布在地表以下25 m的水平线上,炮间距为100 m,第一炮的位置在横向100 m处.每一炮包含200个检波点,均匀分布于地表,检波点间距为100 m,第一个检波点的位置在横向50 m处.利用声波方程正演得到地震记录,震源函数为Ricker子波,主频为4 Hz,采样长度为8000 ms,采样间隔为1 ms.

图 1 真实Overthrust速度模型 Fig. 1 The true Overthrust velocity model

利用施加最佳匹配层(Perfectly Matched Layer,PML)边界条件的声波方程正演合成地震记录.图 2a为第101炮地震记录,对起跳时间后的波形施加阻尼得到初至波形如图 2b所示.

图 2 (a) 第101炮的地震记录(截取至6000 ms);(b)阻尼后得到的第101炮地震记录初至波形,阻尼因子α=30.0,并做了道间归一化处理 Fig. 2 (a) Seismogram of the 101 shot (within 6000 ms); (b) Seismogram contains only the first-arrival waveform by damping figure (a) with the damping factor of 30.0 followed by inter-trace normalization

对得到的初至波形作傅里叶变换,再对其取对数便得到初至相位信息.对第101炮的地震记录进行上述操作,图 3a展示了该地震记录不同频率的初至相位走时,随着频率的增大,周期跳跃的现象越来越严重,且折叠的值域范围也在不断收敛.通过监测相位不连续性和2π周期判定进行相位展开,图 3b为不同频率的非折叠相位(Choi and Alkhalifah, 2013)走时.为了更直观地显示,将不同频率的非折叠相位走时放在一起进行对比.图 4a分别展示了第1炮、第101炮和第151炮所有频率的非折叠相位走时,图 4b图 4a中对应最大偏移距处的放大细节,可见不同频率的非折叠相位走时是有差别的,且偏移距越大差别越大.在勘探频带内,20 km偏移距不同频率的走时差异约为20 ms,但这种差别在极低频时则非常明显,如图 5所示.

图 3 第101炮地震记录1 Hz、4 Hz、8 Hz的(a)折叠初至相位走时与(b)非折叠初至相位走时 Fig. 3 The 1, 4 and 8 Hz (a) wrapped phase traveltime and (b) unwrapped phase traveltime
图 4 (a) 第1炮、第101炮和第151炮所有频率的非折叠相位走时;(b)最大偏移距处的细节放大 Fig. 4 (a) The unwrapped phase traveltimes of 1, 101 and 151 shot gathers for all frequencies; (b) Zoomed views of the traveltimes at the maximum offsets
图 5 极低频时的非折叠相位走时,不同频率会有很大差别 Fig. 5 The unwrapped phase traveltimes for extra-low frequencies. The differences are very apparent
3 阻尼因子的影响

对地震记录施加不同的阻尼因子会对提取的相位走时影响很大.图 6展示了不同阻尼因子下第1炮和第101炮的折叠相位走时,其与频率和偏移距关系明显.可以看到,当阻尼因子很小时(图 6a, c),折叠相位走时的形态十分复杂且不规则,而强阻尼因子下的折叠相位走时(图 6b, d)则比较简单规则. 图 7图 6中的折叠相位走时在偏移距为5 km、15 km(第1炮)和2 km、8 km (第101炮)处的剖面.当阻尼因子很小时(图 7a, c, e, g),折叠相位走时随频率变化很大,而强阻尼因子下的折叠相位走时(图 7b, d, f, h)随频率的变化总体控制在一定范围内.图 8图 6中的折叠相位走时在频率为2 Hz和6 Hz处的剖面.当阻尼因子很小时(图 8a, c, e, g),折叠相位走时在部分位置不平整会出现跳跃,而强阻尼因子下的折叠相位走时(图 8b, d, f, h)十分平整规律.这是因为强阻尼因子可以压制初至波后的“噪声”,从而得到初至波形.

图 6 (a, b)第1炮和(c, d)第101炮记录在(a, c) 0.1 s-1与(b, d) 30 s-1的阻尼因子下获得的折叠相位走时与频率和偏移距的关系 Fig. 6 The wrapped phase traveltimes as a function of frequency and offset of the (a, b) 1 and (c, d) 101 shot gathers for damping factors of (a, c) 0.1 s-1 and (b, d) 30 s-1
图 7 图 6中的折叠相位走时在(a, c, e, g) 0.1 s-1和(b, d, f, h) 30 s-1的阻尼因子下,偏移距分别为(a, b, c, d)第1炮:5 km、15 km和(e, f, g, h)第101炮:2 km、8 km处的剖面 Fig. 7 The profiles of wrapped phase traveltimes map in Fig. 6 at offsets of (a, b, c, d) shot gather 1 (5 km, 15 km) and (e, f, g, h) shot gather 101 (2 km, 8 km) for damping factors: (a, c, e, g) 0.1 s-1, (b, d, f, h) 30 s-1
图 8 图 6中第1炮(a, b, c, d)和第101炮(e, f, g, h)的折叠相位走时在(a, c, e, g) 0.1 s-1和(b, d, f, h) 30 s-1阻尼因子下频率为(a, b, e, f)2 Hz和(c, d, g, h)6 Hz的剖面 Fig. 8 The profiles of wrapped phase traveltimes of the 1 shot (a, b, c, d) and 101 shot (e, f, g, h) map in Fig. 6 at frequencies of (a, b, e, f) 2 Hz and (c, d, g, h) 6 Hz for damping factors of (a, c, e, g) 0.1 s-1 and (b, d, f, h) 30 s-1

通过监测相位不连续性和2π周期判定将相位展开,图 9展示了图 6中的折叠相位展开后的结果.可以看到,阻尼因子很小时(图 9a, c)的非折叠相位走时形态仍是复杂不规则的,而强阻尼因子(图 9b, d)下的非折叠相位走时十分规则,很好地体现了走时信息的变化.图 10图 9中的非折叠相位走时在偏移距为5 km、15 km(第1炮)和2 km、8 km(第101炮)处的剖面.可以发现,强阻尼因子下的非折叠相位走时(图 10b, d, f, h)随频率的变化具有单调性,且不同频率的相位走时差异与偏移距有关,偏移距越大差异越大.图 11图 9中的非折叠相位走时在频率为2 Hz和6 Hz处的剖面.可以看到,阻尼因子很小时(图 11a, c, e, g)未能完全对相位进行展开,仍存在一定的跳跃,而强阻尼因子下(图 11b, d, f, h)的相位展开操作能够得到真实的与频率有关的相位走时信息.

图 9 (a, b)第1炮和(c, d)第101炮在(a, c) 0.1 s-1和(b, d) 30 s-1阻尼因子下与频率、偏移距有关的非折叠相位走时 Fig. 9 The unwrapped phase traveltimes as a function of frequency and offset of (a, b) 1 and (c, d) 101 shot gathers for damping factors of (a, c) 0.1 s-1 and (b, d) 30 s-1
图 10 图 9中的非折叠相位走时在(a, c, e, g) 0.1 s-1和(b, d, f, h) 30 s-1不同阻尼因子下偏移距为:(a, b, c, d)第1炮:5 km、15 km和(e, f, g, h)第101炮:2 km、8 km处的剖面 Fig. 10 The profiles of unwrapped phase traveltimes map in Fig. 9 at offsets of (a, b, c, d) shot gather 1 (5 km, 15 km) and (e, f, g, h) shot gather 101 (2 km, 8 km) for damping factors: (a, c, e, g) 0.1 s-1, (b, d, f, h) 30 s-1
图 11 图 9中的非折叠相位走时在(a, c, e, g)0.1 s-1和(b, d, f, h)30 s-1不同阻尼因子下在频率为(a, b, e, f)2 Hz和(c, d, g, h)6 Hz处的剖面 Fig. 11 The profiles of unwrapped phase traveltimes map in Fig. 9 at frequencies of (a, b, e, f) 2 Hz and (c, d, g, h) 6 Hz for damping factors: (a, c, e, g) 0.1 s-1, (b, d, f, h) 30 s-1
4 层析核函数分析

本文使用的核函数为Rytov近似下的频率域表达式(6).为了显示不同频率下核函数的区别,本节做了下述数值实验.图 12av=2000 m·s-1的均匀速度模型,包含201×201个网格,网格大小为10 m×10 m,左边蓝星表示激发点位置,右边蓝星表示接收点位置.图 12b12d展示了不同频率的核函数.可以看到,在低频时,菲涅尔体胖度较大,能量分布在整个菲涅尔体范围内,随着频率的增高,菲涅尔体胖度逐渐减小,能量也逐渐向中心射线集中.对2~100 Hz之间每2 Hz离散采样的单频层析核函数进行加权叠加,加权函数为Ricker子波的振幅谱,得到带限层析核函数KT如下所示:

(8)

图 12 (a) 简单均匀速度模型,v=2000 m·s-1;(b)—(d)不同频率的单频层析核函数;(e) 0~100 Hz之间每2 Hz离散采样的单频层析核函数进行加权叠加得到的带限层析核函数(a)左边蓝星为激发点,位置在(250,1000)m处,右边蓝星为接收点,位置在(1750,1000)m处;(b) 4 Hz核函数;(c) 10 Hz核函数;(d) 50 Hz核函数;(e)带限层析核函数 Fig. 12 (a) A homogeneous velocity model with v=2000 m·s-1; (b)—(d) are different monochromatic sensitivity kernels; (e) Band-limited sensitivity kernel by weighted-stacking the 0~100 Hz monochromatic kernels (a) The left blue star indicates the shot position at (250, 1000) m, and the right one denotes the receiver position at (1750, 1000) m; (b) 4 Hz; (c) 10 Hz; (d) 50 Hz; (e) Band-limited

其中,KT(r, ω)为单频层析核函数.加权函数P(ω)为

(9)

其中,ω0=2πf0f0为Ricker子波的主频.

最后得到的带限层析核函数剖面如图 12e所示.图 13x=1000 m处的垂向剖面,可以发现单频层析核函数在中心射线的一定领域范围内基本能够进行同相叠加,而在该范围之外则因异相叠加而相互抵消.由于在主频对应的菲涅尔体范围内,绿线和红线能够很好地吻合,所以在常规的菲涅尔体层析中常利用主频对应的单频菲涅尔体边界近似代替带限菲涅尔体边界(Liu et al., 2009).然而,不同频率的层析核函数不尽相同,利用主频对应的单频菲涅尔体边界约束带限菲涅尔体会导致一定的反演误差.本文提出的相位走时层析方法使用所有频率的全空间(而非限于菲涅尔体内)层析核函数信息,且利用频率依赖的相位信息,从低频到高频,不同频率下都进行了迭代反演,从而保证可以得到更精确的反演结果.

图 13 (a) 根据公式(6)计算得到的单频层析核函数x=1000 m处垂直剖面;(b)红线为对所有单频层析核函数进行振幅谱加权叠加得到的带限层析核函数x=1000 m处垂直剖面,绿线为30 Hz主频对应的单频层析核函数x=1000 m处垂直剖面,均做了归一化处理 Fig. 13 Extracted sections of (a) the monochromatic sensitivity kernels according to formula (6), and (b) the red curve is the weighted-stacking result of the monochromatic sensitivity kernels, and the green curve is the monochromatic kernel of the dominant frequency 30 Hz, at x=1000 m. Both of them are normalized
5 Overthrust模型的初至波相位走时层析实验

使用Overthrust模型来验证初至波相位走时层析方法的有效性.理论模型如图 1所示,初始模型如图 14所示,为v=v0v·z的等梯度模型,其中v0=2300 m·s-1,Δv=22.581 m·s-1z为模型深度.

图 14 初始梯度模型 Fig. 14 The constant-gradient initial velocity model

反演选取的频带范围为1~8 Hz,每个频率迭代10次.图 15图 17分别展示了1 Hz、4 Hz和8 Hz的反演结果和对应频率最后一次迭代的梯度.分析可以发现,随着频率的增大,迭代次数的增多,反演精度也在不断提高,特别是地表浅部的速度结构尤为明显;从梯度上看,低频时梯度更新量集中在速度模型的浅部,而随着频率的增大,梯度更新量逐渐往模型深部发展,梯度值逐渐减小,且对模型局部速度的更新越来越精细.图 18为抽取不同频率反演结果在地下200 m、300 m处的速度剖面,更直观地显示了在近地表区域随着频率的增大,反演精度与分辨率也在不断提升.

图 15 (a) 1 Hz初至波相位走时层析反演结果;(b)最后一次迭代对应梯度 Fig. 15 The first-arrival phase-traveltime tomographic result after 1 Hz inversion, and (b) the corresponding gradient in the last iteration
图 16 (a) 4 Hz初至波相位走时层析反演结果;(b)最后一次迭代对应梯度 Fig. 16 (a) The first-arrival phase-traveltime tomographic result after 4 Hz inversion, and (b) the corresponding gradient in the last iteration
图 17 (a) 8 Hz初至波相位走时层析反演结果;(b)最后一次迭代对应梯度 Fig. 17 (a) The first-arrival phase-traveltime tomographic result after 8 Hz inversion, and (b) the corresponding gradient in the last iteration
图 18 (a) 深度200 m处的速度剖面;(b)深度300 m处的速度剖面 Fig. 18 Velocity sections at depth of (a) 200 m and (b) 300 m

为了验证初至波相位走时层析的反演能力,本文将其与传统射线走时层析和菲涅尔体走时层析(Liu et al., 2009)作对比.图 19a为传统射线走时层析的反演结果,图 19b为菲涅尔体走时层析的反演结果.从反演结果可以看出,初至波相位走时层析与菲涅尔体走时层析都准确地反演出了模型中的中低波数扰动信息,对该模型近地表的速度结构都有一定程度的更新.相比于传统射线走时层析,相位走时层析(图 17a)浅层低速异常反演更加准确.尽管随着深度的加深,初至波相位走时层析方法的反演精度逐渐降低,但总体上还是优于射线走时层析的,能大致反映真实速度的变化趋势.图 20为三种方法的目标函数收敛曲线,可以看出,初至波相位走时层析方法收敛得更快,且在不同频率会出现阶梯状下降的现象.尽管在高频时目标函数值会有轻微上升,但整体趋于稳定,且最终收敛值更小.为了进行细节对比,分别抽取地表以下100 m、200 m、300 m处的速度剖面,如图 21abc所示.可以看到,传统射线走时层析的反演速度变化极为平滑,且容易受深部高速体的影响,导致在部分浅部反演速度过大;相位走时层析和菲涅尔体走时层析的反演结果也出现该现象,但总体控制得更好.在横向速度剧烈变化的部分区域,菲涅尔体走时层析的反演速度变化较为平缓,而相位走时层析反演的速度在这些位置能体现出相应的变化.因此,本文提出的相位走时层析方法更加稳定,反演结果更精确,特别是对于速度模型的局部细节,处理得更好.另一方面,初至波相位走时层析与菲涅尔体走时层析在表层的反演结果基本相近,一定程度上说明了在地震勘探的频带范围内,利用射线走时构建目标函数仍是较为精确的一种方法,这样可以避免时间域模拟,进而大幅降低计算量.

图 19 (a) 射线走时层析反演结果;(b)菲涅尔体走时层析反演结果 Fig. 19 Inversion results of (a) ray-based traveltime tomography and (b) Fresnel volume traveltime tomography
图 20 三种方法的目标函数收敛曲线 绿线表示射线走时层析,蓝线表示菲涅尔体走时层析,红线表示初至波相位走时层析 Fig. 20 Convergence curves of the three methods Green line indicates the ray-based traveltime tomography, blue line denotes the Fresnel volume traveltime tomography and the red line represents the phase-traveltime tomography
图 21 (a)、(b)、(c)为三种方法分别在地表以下100 m、200 m、300 m处的速度剖面对比结果 黄线:真实模型;黑线:初始模型;红线:相位走时层析结果;蓝线:菲涅尔体走时层析结果;绿线:射线走时层析结果 Fig. 21 Velocity sections at depths of (a) 100 m, (b) 200 m and (c) 300 m of the three methods Yellow line: the true model; Black line: the initial model; Red line: the phase-traveltime tomographic result; Blue line: the Fresnel volume traveltime tomographic result; Green line: the ray-based traveltime tomographic result
6 结论

本文提出了一种基于相位反演的初至波相位走时层析成像方法,该方法的核心是:①利用依赖于频率的相位走时信息,而非单一的无限频率射线走时;②通过监测相位不连续性和2π周期判定来消除相位折叠现象.通过对原始数据进行阻尼变换和取对数,定义并提取了与频率有关的初至相位信息.实践表明:随着频率的增大,相位周期跳跃现象越来越严重,而通过监测相位的不连续性和2π周期判定可以对相位进行有效展开.同时,实验证明:强阻尼应用下的相位展开操作能够得到真实的与频率有关的相位走时信息.相比于射线走时层析,相位走时层析利用基于波动方程导出的核函数而非射线路径,考虑了波传播的有限频特征,具有更高的反演精度;相比于菲涅尔体走时层析,初至波相位走时层析利用频率依赖的无折叠走时,利用了更多的初至波信息,并利用全空间的核函数信息,改善了反演质量.它们在表层的反演结果相近,一定程度上说明了在勘探地震的频带范围内,利用射线走时残差构建波动方程走时层析目标函数仍是较为精确的一种方法,且可以避免时间域模拟,大幅降低计算量.

射线走时和瞬时走时是目前文献中利用走时信息构建目标函数的两种方法,它们使用的核函数不尽相同.与这两种走时相比,本文提出了相位走时的概念实际上是一种频率相关的初至走时.关于射线走时、瞬时走时和相位走时这三者之间的区别和联系,是今后需要深入探讨的问题.

References
Bednar J B, Shin C, Pyun S. 2007. Comparison of waveform inversion, Part 2:Phase approach. Geophysical Prospecting, 55(4): 465-475. DOI:10.1111/gpr.2007.55.issue-4
Bishop T N, Bube K P, Cutler R T, et al. 1985. Tomographic determination of velocity and depth in laterally varying media. Geophysics, 50(4): 903-923.
Choi Y, Alkhalifah T. 2013. Frequency-domain waveform inversion using the phase derivative. Geophysical Journal International, 195(3): 1904-1916. DOI:10.1093/gji/ggt351
Gélis C, Virieux J, Grandjean G. 2007. Two-dimensional elastic full waveform inversion using Born and Rytov formulations in the frequency domain. Geophysical Journal International, 168(2): 605-633. DOI:10.1111/gji.2007.168.issue-2
Liu Y Z, Dong L G, Li P M, et al. 2009. Fresnel volume tomography based on the first arrival of the seismic wave. Chinese Journal of Geophysics (in Chinese), 52(9): 2310-2320. DOI:10.3969/j.issn.0001-5733.2009.09.015
Liu Y Z, Dong L G, Wang Y W, et al. 2009. Sensitivity kernels for seismic Fresnel volume tomography. Geophysics, 74(5): U35-U46. DOI:10.1190/1.3169600
Liu Y Z, Dong L G. 2012. Influence of wave front healing on seismic tomography. Science China Earth Sciences, 55(11): 1891-1900. DOI:10.1007/s11430-012-4441-0
Liu Y Z, Yang J Z, Chi B X, et al. 2015. An improved scattering-integral approach for frequency-domain full waveform inversion. Geophysical Journal International, 202(3): 1827-1842. DOI:10.1093/gji/ggv254
Oppenheim A V, Schafer R W. 1975. Digital Signal Processing. Inglewood Gliffs, N J: Prentice Hall.
Poggiagliolmi E, Berkhout A J, Boone M M. 1982. Phase unwrapping, possibilities and limitations. Geophysical Prospecting, 30(3): 281-291. DOI:10.1111/gpr.1982.30.issue-3
Sheng J M, Leed A, Buddensiek M, et al. 2006. Early arrival waveform tomography on near-surface refraction data. Geophysics, 71(4): 47-57.
Shin C, Min D J, Marfurt K J, et al. 2002. Traveltime and amplitude calculations using the damped wave solution. Geophysics, 67(5): 1637-1647. DOI:10.1190/1.1512811
Shin C, Min D J. 2006. Waveform inversion using a logarithmic wavefield. Geophysics, 71(3): R31-R42. DOI:10.1190/1.2194523
Spetzler G, Snieder R. 2004. The Fresnel volume and transmitted waves. Geophysics, 69(3): 653-663. DOI:10.1190/1.1759451
Wielandt E. 1987. On the validity of the ray approximation for interpreting delay times.//Nolet G ed. Seismic Tomography. Dordrecht: Springer, 85-98.
Woodward M J. 1992. Wave-equation tomography. Geophysics, 57(1): 15-26. DOI:10.1190/1.1443179
Zhu X H, Sixta D P, Angstman B G. 1992. Tomostatics: Turning-ray tomography+static corrections.//62nd Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts.
刘玉柱, 董良国, 李培明, 等. 2009. 初至波菲涅尔体地震层析成像. 地球物理学报, 52(9): 2310-2320. DOI:10.3969/j.issn.0001-5733.2009.09.015