在勘探地震的偏移成像中,近地表速度(P波速度和S波速度)结构(刘伊克等, 2001)的精确度制约着深层地质结构的构建.目前为止,大多数面波反演(Xia et al., 1999; 石耀霖和金文, 1995; 胡家富等, 1999; 李明明和何玉梅, 2011)基于层状模型假设,这一假设在实际复杂的地质环境中并不适用.为了克服这个层状模型假设,许多研究人员尝试采用波形反演的思路对面波进行反演.Schäfer等(2013)采用全波形反演方法对两个实际数据进行测试,结果表明反演过程收敛速度极慢以及反演结果往往陷入局部最小值.Solano等(2013)提出了一种振幅波形反演(Amplitude Waveform Inversion)的方法,其目标函数为观测数据振幅谱与预测数据振幅谱的最小二范数,而传统的全波形反演(Tarantola, 1984; Pratt et al., 1998; Virieux and Operto, 2009)的目标函数为观测数据与预测数据的最小二范数,结果表明振幅波形反演具有更好的收敛性质.这是因为传统的全波形反演匹配的是所有波形,而振幅波形反演只匹配振幅谱,振幅谱仅仅是所有波形信息中的一部分,这种简化的数据(Luo and Schuster, 1991a, 1991b)使得振幅波形反演具有更快的收敛速度和更好的收敛效果.Solano等(2014)提出了一种加窗振幅波形反演方法(Windowed-Amplitude Waveform Inversion),与振幅波形反演相比,该方法更加稳定.然而与传统全波形反演相比,振幅波形反演是以牺牲分辨率为代价的.为了部分恢复振幅波形反演丢失的分辨率,本文提出了一种初至P波与瑞雷面波的交叉梯度联合反演(于鹏等, 2009; 彭淼等, 2013; 李兆祥等, 2015)方法,采用交叉梯度同时作为初至P波反演和瑞雷面波加窗振幅波形反演的约束项,将分辨率较低的瑞雷面波振幅波形反演与高分辨率的P波波形反演结合起来,可以获得更高分辨率的S波速度结构.
地球物理数据的联合反演可以分为两大类.第一类是基于不同参数之间的岩石物理关系,例如饱和度和孔隙度(Hoversten et al., 2006),电导率和介电常数(Kowalsky et al., 2005).例如,P波和S波速度之间的关系是依赖于介质属性的.Castagna等(1985)指出,在砂岩中P波与S波比约为1.5~1.7,在页岩中该比值大于2.这就意味着在联合反演中,采用不同参数之间错误的经验关系可能会带来严重的耦合误差.
第二类联合反演利用不同参数之间的结构相似性,同一地质模型可以用不同的物性参数来描述,例如密度、P波速度、电导率等,这些参数存在着结构相似性.Gallardo和Meju (2003)采用交叉梯度联合反演方法成功地对直流电法和地震走时数据进行了反演.Colombo和De Stefano (2007)等利用交叉梯度联合反演方法同时对地震、重力和大地电磁数据进行了反演.本文采用第二类联合反演策略对初至地震P波和瑞雷面波进行交叉梯度联合反演.
对于初至纵波,本文采用全波形反演获得地下P波速度结构,采用加窗振幅波形反演方法对瑞雷面波进行反演,同时采用P波速度和S波速度的交叉梯度作为二者的约束项.结果显示,与单独的瑞雷面波振幅波形反演相比,交叉梯度联合反演的S波速度具有更高的精度.本文第2节和第3节将分别给出交叉梯度联合反演的方法原理以及数值模拟结果.之后,将该联合反演方法应用于实际数据.最后是结论.
2 方法原理对于地震数据,初至P波可以通过求解二维常密度声波方程来进行模拟,
(1) |
其中v是纵波速度,p(x, t)是压力场,f(x, t)是震源项.
另一方面,瑞雷面波可以通过求解弹性波方程进行模拟获得,
(2) |
其中ρ是密度,ui是第i个位移分量,cijkl是刚度矩阵,si是震源项.
2.1 交叉梯度交叉梯度的定义为两个模型参数梯度的叉积:
(3) |
其中mP代表P波速度,而mS代表S波速度.交叉梯度正则项是基于这样的事实,当两个模型参数的交叉梯度达到最小化时,它们的结构相似度将达到最大化.如果在P波和S波速度中存在相同的边界,且该边界的位置和方向一致,那么它们的交叉梯度将等于零.
对于二维问题,方程(3)可以改写为:
(4) |
如图 1a为P波速度模型,波速从左到右逐渐增大;图 1b为S波速度模型,一个低速方形异常体、一个高速圆形异常体及一个高速方形异常体嵌于一均匀背景速度模型;图 1c为二者的交叉梯度值.以左侧低速异常体为例,虽然在S波速度中四条边界处的梯度值均不为零;然而在P波速度中低速方形所在位置的左、右边界处的梯度与S波速度中对应的梯度平行,故而交叉梯度值为零,仅在上下边界处其交叉梯度值不为零.
波动方程波形反演的目标函数为:
(5) |
其中Δdrs=dcal(xr, xs)-dobs(xr, xs)是初至P波残差,它表示为预测数据与观测数据之差.对于P波速度,它可以通过任何一种梯度类方法进行迭代更新,例如通过共轭梯度(Hestenes and Stiefel, 1952; 胡祖志, 2006)方法对P波速度进行更新:
(6) |
其中vkP为第k次P波速度模型,αkP为迭代步长,dkP为迭代更新方向:
(7) |
(8) |
当前的迭代搜索方向dkP取决于目标函数对模型参数的梯度以及上一次的搜索方向dk-1P.目标函数对模型参数的梯度表示为:
(9) |
其中,
(10) |
符号*代表时间卷积,
加窗振幅波形反演(Solano et al., 2014)是在全波形反演的基础上衍生出来的,它的目标函数为:
(11) |
(12) |
(13) |
其中, w代表不同窗口,Dwcal和Dwobs分别为加窗预测数据和观测数据的二维傅里叶变换,dcal为预测数据,dobs为观测数据,wwT为窗函数.与P波速度类似,S波速度同样可以采用共轭梯度方法进行迭代更新:
(14) |
其中vkS代表第k次的S波速度模型;αkS是迭代步长;dkS代表迭代更新方向,它取决于上一次更新迭代方向dk-1S及目标函数的梯度:
(15) |
其中
(16) |
uz(x, t| xs)cal为预测的垂直位移分量,GS(x, -t| xr, 0)为方程(2)格林函数.
(17) |
其中
为了提高反演结果的精度,本文将初至P波与S波速度进行联合同步反演,对于初至P波,其全波形联合反演目标函数为:
(18) |
同理对于瑞雷面波,其加窗振幅波形联合同步反演的目标函数则为:
(19) |
其中λ1和λ2分别为各自的正则化系数,当这两个正则化系数同时置为零时,方程(18)和(19)退化为无耦合的单独反演方程.对于交叉梯度联合反演,关键问题是计算目标函数(18)和(19)分别对P波速度和S波速度的偏导,即获得各自的梯度.注意到交叉梯度联合反演目标函数均包含两项,原始目标函数项和正则约束项,而原始目标函数对应的梯度已知,故关键问题是求取正则约束项分别对P波速度和S波速度的梯度,根据Hu等(2009),经过数学推导,正则项的梯度分别表示为:
(20) |
(21) |
全波形反演是病态的非线性反演问题,非线性、跳周、子波估计不准等现象往往使得波形反演问题陷入局部最小值.本文选取加窗振幅波形反演对瑞雷面波进行反演,是因为与全波形反演相比,它更容易收敛到全局最小值.为了比较全波形反演、振幅波形反演及加窗振幅波形反演三种目标函数的形态,首先构建一简单模型如图 2所示,图 2为S波速度模型,一条带状异常体速度vS1嵌于均匀背景速度vS0中,真实速度模型和初始速度模型同时取决于vS1和vS0.假设真实背景速度大小为800 m·s-1,真实异常体速度大小为960 m·s-1;通过改变vS1和vS0获得不同的初始模型,其中vS1的变化范围为400~1400 m·s-1,vS0的变化范围为600~1400 m·s-1.同时P波速度设为2000 m·s-1,密度设为1000 kg·m-3.图中红色星号代表震源所在位置,采用中心频率为20 Hz的雷克子波作为震源;检波器均匀分布于地表.图 3为三种不同目标函数的形态,十字虚白线交叉点为真模型;图 3a为全波形反演目标函数形态;图 3b为振幅波形反演目标函数形态;图 3c为加窗振幅波形反演目标函数形态,可见加窗目标函数具有最宽的全局最小范围,最容易收敛到全局最小.相比于初至P波,瑞雷面波具有更强的非线性性质,故而本文对瑞雷面波采用加窗振幅波形反演方法;而对初至P波采用全波形反演方法.
首先将交叉梯度联合反演方法应用于简单数值模型图 4.其中图 4a是P波速度模型,图 4b为S波速度模型.该模型包括三层岩层,从地表到深部其速度逐渐递增,且不同岩层之间界面为不规则界面.模型被离散为30×120的等体积立方体网格,网格间距为3 m.图中红色X表示震源所在位置,其间距为6 m;地表均匀布设120个检波器,检波器间距为3 m.源子波为中心频率为20 Hz的雷克子波.初始速度模型见图 5,它是将真实速度模型进行50次高斯滤波圆滑获得.
图 6a是经过10次迭代后单独反演的P波速度结果,由图可知,P波反演结果中基本可以识别出三层岩层,其中浅层形态与真实模型吻合较好;底层高速层的形态及反演结果与真实模型基本一致;中间层的形态与真实模型基本一致,然而其反演结果与真实模型存在一定误差.图 6b是经过10次迭代后交叉梯度联合反演P波速度结果,其结果同样并没有得到提高.
图 7所示为经过10次迭代反演S波速度结果.图 7a为单独反演S波速度结果,由图可知,浅层速度结构刻画得较好,浅层与中间层的界面与真实模型基本吻合,但是中间层及底层的反演结果与真实模型存在较大误差.这是因为瑞雷面波的穿透深度仅为其波长的1/3至1/2,而这里浅层的速度为800 m·s-1,对应的波长为40 m,故而其探测深度应在20 m范围内.图 7b为交叉梯度联合反演S波速度结果,可以看到:交叉梯度联合反演S波速度结果中,浅层界面的连续性好于单独反演结果.
此二维地震勘探的目的是为了探明地下某隐伏断层在地表的位置.该二维地震剖面包括240个垂直分量检波器,检波器间距为5 m.时间采样率为1 ms,记录长度为3 s,每个地震记录进行32次叠加,总共240炮地震记录.图 8a,8b,8c,8d分别为偏移距为10 m,30 m,50 m,70 m的共偏移距地震剖面.通过对比不同共偏移距剖面,我们认为该隐伏断层可能位于第55炮所在位置,即位于水平位置275 m处,即如图中红色虚线所示.故本文选取第21炮至第81炮地震记录进行反演.
图 9a为第21炮原始地震记录,图中红色虚线所示为初至P波,而蓝色虚线所示为瑞雷面波.通过对两条虚线的斜率进行估算,可以得到近地表的P波速度和S波速度分别约为2900 m·s-1和460 m·s-1.在数据预处理之前,首先手动提取初至P波和瑞雷面波走时,提取的走时将用于后续提取初至P波和瑞雷面波波形.
数据预处理的第一步是选取3~50 Hz滤波器对原始数据进行带通滤波;其次对数据进行三维到二维转换及几何扩散校正;第三步去除坏道并采用POCS (Projection onto Convex Sets)(Abma and Kabir, 2006)对缺失道进行数据插值;第四步在时间域对数据进行插值,使得时间采样率满足正演模拟需求;最后切除近偏移距地震道以及根据提取的走时将初至P波和瑞雷面波提取出.经过预处理后的数据如图 9b所示.
4.2 反演结果通过对瑞雷面波频散曲线进行反演(Zhang et al., 2015)获得的S波速度可以作为联合反演的初始S波速度模型,如图 10b所示,该初始S波速度模型被离散为30×121等体积立方体网格,网格间距为2.5 m.从图可知,在30 m深度范围内,从左到右,S波速度从420 m·s-1递增到530 m·s-1,低速异常体位于水平方向100~150 m,这意味着断层地表水平位置可能位于150~170 m;深层为一高速带,故在45 m以下的区域设置为650 m·s-1常速度.基于此S波速度模型,以及从图 9a估算的P波速度值(2900 m·s-1),我们可以构建初始P波速度模型,如图 10a所示.
图 11a是经过15次迭代后单独反演P波速度结果,图 11b是经过15次迭代后交叉梯度联合反演P波速度结果.由图可知,相比于单独反演,交叉梯度联合反演P波速度结果并没有得到提高.然而与初始P波速度相比,反演结果确有很大提高,其中低速异常体面积扩大,从中可以推断隐伏断层水平位置可能位于230~280 m,这与共偏移距结果指示的隐伏断层位置更加吻合.
图 12所示为经过15次迭代反演S波速度结果.图 12a为单独反演S波速度结果,可见在水平位置100~200 m处仍然存在一低速异常体,但是该异常体的形态与P波反演结果形态略有不同.并且在浅部水平位置250~300 m处存在一高速异常体,这与P波反演结果不同.图 12b为交叉梯度联合反演S波速度结果,该结果中低速异常体的位置、形态与P波速度反演结果非常相似,而且浅部也未见局部高速异常体.单独反演中,S波速度反演结果与P波速度反演结果相似度不高;而在交叉梯度联合反演中,S波速度与P波速度反演结果相似度极高,其反演结果更真实可靠.
本文提出了一种初至P波和瑞雷面波联合反演的方法,对初至P波进行全波形反演,以及对瑞雷面波采用加窗振幅波形反演,同时将P波速度和S波速度的交叉梯度作为二者反演的正则化约束项.数值模拟结果表明,与单独反演相比,联合反演可以消除S波反演结果中的伪异常,提高S波反演精度;对于P波速度,联合反演并不能获得更好的反演结果.实际资料反演结果表明,交叉梯度联合反演的S波速度与P波速度反演结果相似度更高,其反演结果更可信.
交叉梯度联合反演的缺点是,当真实的P波速度模型与S波速度模型在结构上相似度不高时,采用交叉梯度作为反演约束项会得到错误的反演结果.
Abma R, Kabir N. 2006. 3D interpolation of irregular data with a POCS algorithm. Geophysics, 71(6): E91-E97. DOI:10.1190/1.2356088 | |
Castagna J P, Batzle M L, Eastwood R L. 1985. Relationships between compressional-wave and shear-wave velocities in clastic silicate rocks. Geophysics, 50(4): 571-581. DOI:10.1190/1.1441933 | |
Colombo D, De Stefano M. 2007. Geophysical modeling via simultaneous joint inversion of seismic, gravity, and electromagnetic data: Application to prestack depth imaging. The Leading Edge, 26(3): 326-331. DOI:10.1190/1.2715057 | |
Gallardo L A, Meju M A. 2003. Characterization of heterogeneous near-surface materials by joint 2D inversion of DC resistivity and seismic data. Geophysical Research Letters, 30(13): 1658. | |
Hestenes M R, Stiefel E. 1952. Methods of conjugate gradients for solving linear systems. Journal of Research of the National Bureau of Standards, 49(6): 409-436. DOI:10.6028/jres.049.044 | |
Hoversten G M, Cassassuce F, Gasperikova E, et al. 2006. Direct reservoir parameter estimation using joint inversion of marine seismic AVA and CSEM data. Geophysics, 71(3): C1-C13. DOI:10.1190/1.2194510 | |
Hu J F, Duan Y K, Hu Y L, et al. 1999. Inversion of Shear-wave velocity structure in shallow soil from Rayleigh waves. Chinese J. Geophys. (in Chinese), 42(3): 393-400. | |
Hu W Y, Abubakar A, Habashy T M. 2009. Joint electromagnetic and seismic inversion using structural constraints. Geophysics, 74(6): R99-R109. DOI:10.1190/1.3246586 | |
Hu Z Z, Hu X Y, He Z X. 2006. Pseudo-Three-Dimensional magnetotelluric inversion using nonlinear conjugate gradients. Chinese J. Geophys. (in Chinese), 49(4): 1226-1234. | |
Kowalsky M B, Finsterle S, Peterson J, et al. 2005. Estimation of field-scale soil hydraulic and dielectric parameters through joint inversion of GPR and hydrological data. Water Resources Research, 41(11): W11425. DOI:10.1029/2005WR004237 | |
Li M M, He Y M. 2011. Lithospheric structure beneath northeastern boundary region of the north China craton from Rayleigh wave dispersion inversion. Acta Seismologica Sinica (in Chinese), 33(2): 143-155. | |
Li Z X, Tan H D, Fu S S, et al. 2015. Two-dimensional synchronous inversion of TDIP with cross-gradient constraint. Chinese J. Geophys. (in Chinese), 58(12): 4718-4726. DOI:10.6038/cjg20151232 | |
Liu Y K, Chang X, Wang H, et al. 2001. Estimation of near-surface velocity and seismic tomographic static corrections. Chinese J. Geophys. (in Chinese), 44(2): 272-278. | |
Luo Y, Schuster G T. 1991a. Wave equation inversion of skeletalized geophysical data. Geophysical Journal International, 105(2): 289-294. DOI:10.1111/gji.1991.105.issue-2 | |
Luo Y, Schuster G T. 1991b. Wave-equation traveltime inversion. Geophysics, 56(5): 645-653. DOI:10.1190/1.1443081 | |
Peng M, Tan H D, Jiang M, et al. 2013. Three-dimensional joint inversion of magnetotelluric and seismic travel time data with cross-gradient constraints. Chinese J. Geophys. (in Chinese), 56(8): 2728-2738. DOI:10.6038/cjg20130821 | |
Pratt R G, Shin C, Hick G J. 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion. Geophysical Journal International, 133(2): 341-362. DOI:10.1046/j.1365-246X.1998.00498.x | |
Schäfer M, Groos L, Forbriger T, et al. 2013. 2D Full waveform inversion of recorded shallow seismic Rayleigh waves on a significantly 2D structure.//Near Surface Geoscience 2013-19th EAGE European Meeting of Environmental and Engineering Geophysics. EAGE. | |
Shi Y L, Jin W. 1995. Genetic algorithms inversion of lithospheric structure from surface wave dispersion. Chinese J. Geophys. (Acta Geophysica Sinica) (in Chinese), 38(2): 189-198. | |
Solano C A P, Donno D, Chauris H. 2013. 2D surface wave inversion in the F-K domain.//75th EAGE Conference & Exhibition incorporating SPE EUROPEC. EAGE. | |
Solano C A P, Donno D, Chauris H. 2014. Alternative waveform inversion for surface wave analysis in 2-D media. Geophysical Journal International, 198(3): 1359-1372. DOI:10.1093/gji/ggu211 | |
Tarantola A. 1984. Inversion of seismic reflection data in the acoustic approximation. Geophysics, 49(8): 1259-1266. DOI:10.1190/1.1441754 | |
Virieux J, Operto S. 2009. An overview of full-waveform inversion in exploration geophysics. Geophysics, 74(6): WCC1-WCC26. DOI:10.1190/1.3238367 | |
Xia J H, Miller R D, Park C B. 1999. Estimation of near-surface shear-wave velocity by inversion of Rayleigh waves. Geophysics, 64(3): 691-700. DOI:10.1190/1.1444578 | |
Yu P, Dai M G, Wang J L, et al. 2009. Joint inversion of magnetotelluric and seismic data based on random resistivity and velocity distributions. Chinese J. Geophys. (in Chinese), 52(4): 1089-1097. | |
Zhang Z D, Liu Y K, Schuster G. 2015. Wave equation inversion of skeletonized surface waves.//2015 SEG Annual Meeting. New Orleans, Louisiana: Society of Exploration Geophysicists. | |
胡家富, 段永康, 胡毅力, 等. 1999. 利用Rayleigh波反演浅土层的剪切波速度结构. 地球物理学报, 42(3): 393–400. | |
胡祖志, 胡祥云, 何展翔. 2006. 大地电磁非线性共轭梯度拟三维反演. 地球物理学报, 49(4): 1226–1234. | |
李明明, 何玉梅. 2011. 利用瑞雷面波反演华北克拉通东北部边界的岩石圈结构. 地震学报, 33(2): 143–155. | |
李兆祥, 谭捍东, 付少帅, 等. 2015. 基于交叉梯度约束的时间域激发极化法二维同步反演. 地球物理学报, 58(12): 4718–4726. DOI:10.6038/cjg20151232 | |
刘伊克, 常旭, 王辉, 等. 2001. 三维复杂地形近地表速度估算及地震层析静校正. 地球物理学报, 44(2): 272–278. | |
彭淼, 谭捍东, 姜枚, 等. 2013. 基于交叉梯度耦合的大地电磁与地震走时资料三维联合反演. 地球物理学报, 56(8): 2728–2738. DOI:10.6038/cjg20130821 | |
石耀霖, 金文. 1995. 面波频散反演地球内部构造的遗传算法. 地球物理学报, 38(2): 189–198. | |
于鹏, 戴明刚, 王家林, 等. 2009. 电阻率和速度随机分布的MT与地震联合反演. 地球物理学报, 52(4): 1089–1097. | |