2. 同济大学海洋与地球科学学院, 波现象与反演成像研究组, 上海 200092
2. School of Ocean and Earth Science, Wave Phenomena and Inversion Imaging Group(WPI), Tongji University, Shanghai 200092, China
叠加速度分析技术(Taner and Koehler, 1969;Davis,1972)是经典的初始速度建模方法,主要由速度谱的计算和拾取两部分构成.最经典且直观的速度谱计算方法是叠加能量法,对所有t0时刻,用一系列速度计算反射同相轴的走时,并将反射走时所对应的能量(振幅平方)叠加得到叠加能量速度谱.然而叠加能量法分辨率较低,因此可以用相似系数(Neidell and Taner, 1971)代替能量叠加准则,得到相似系数速度谱(简称相似系数谱).为了进一步提高速度谱的分辨率以及抗噪声能力,涌现了各种加权相似系数谱(Luo and Hale, 2012;Chen et al., 2015).然而,当CMP道集中存在AVO效应,尤其是极性反转时,这些加权相似系数方法的效果会降低甚至失效.为此,Sarkar等(2002),Fomel(2009)针对振幅极性的空间变化提出了相应的改进措施(AB相似系数等).此外,图像的特征结构(Biondi and Kostov, 1989;Kirlin,1992)、MUSIC(Multiple Signal Classification)方法(Barros et al., 2015)、bootstrap策略(Sacchi,1998;Abbad and Ursin, 2012)等也常用于构造速度谱.
另一方面,速度谱的拾取也是一个重要问题,且需要较好的地质与地球物理背景.速度谱的拾取通常有两种方式:人工拾取及自动拾取.人工拾取方法受人的主观影响大,通常拾取效果较好但需要较长的处理周期.而自动拾取方法通常借助如图像识别(Beveridge et al., 2002)、神经网络等方法(Schmidt et al., 1992;Fish and Kusuma, 1994;Huang and Yang, 2015)实现自动速度谱拾取,但难以与速度谱的计算过程自动融合,仍然需要一定程度的人工干预及质量控制.
本文从模型参数估计的角度重新审视叠加速度估计问题,将叠加速度谱的估计与Radon变换相类比.传统的Radon谱可以在L2范数意义下反演(Trad et al., 2003),也可以在L1或者L0范数意义下反演稀疏的Radon谱(王雄文和王华忠,2014).基于稀疏反演的思想,本文将重点讨论如何将叠加速度分析转变为模型参数估计问题,从而实现叠加速度的自动估计.
2 理论 2.1 正问题:CMP道集预测模型CMP道集dCMP(t, x)可以用不同反射同相轴的叠加近似为:
(1) |
其中,n为同相轴的个数,ui(t, x)为第i个同相轴,t为时间,x为偏移距.
在层状介质中,CMP道集中的反射同相轴的走时平方可以展开为偏移距的偶数次多项式求和(Taner and Koehler, 1969).在中小偏移距的假设下,时距关系可以进一步简化为双曲线.若仅考虑CMP道集的运动学特征,反射同相轴ui(t, x)可以近似为时移的零偏移距反射子波:
(2) |
其中,wi(t)为第i个同相轴对应的零偏移距反射子波,ΔTi(x)为第i个同相轴在偏移距x处相对于零偏移距子波的时移量(即正常时差校正量,NMO时差).然而,(2) 式无法考虑随空间位置变化的反射波振幅以及子波的相位空间变化.因此,可以引入振幅的空间变化预测函数,此时同相轴表示为:
(3) |
结合(1) 和(3) 式,CMP道集预测模型可以近似为:
(4) |
(4) 式的物理解释为:反射同相轴可以用时移ΔTi(x)后的反射子波wi(t)表达,且空变振幅用Ai(x)近似.
对于层状介质,反射波的时距关系存在解析表达式.然而对于横向变速以及倾斜界面,难以给出解析的反射波走时关系.因此,可以将NMO时移量ΔT分解为解析部分ΔTa和扰动部分Δtp:
(5) |
通过合理选择反射走时的近似表达式(ΔTa(x)存在解析表达),如经典的双曲速度关系,或考虑各向异性的4阶走时预测函数(Alkhalifah and Tsvankin, 1995)等,可以使得Δtp(x)远小于ΔTa(x).
本文将(4) 式给出的正问题定义为CMP道集预测模型,可以用算子形式简记为:d=Lw.其中算子L表示(4) 式中的时移、振幅和走时校正以及同相轴叠加过程.其中,子波wi(t)、预测误差Δti(x)、空变振幅Ai(x)为模型参数,而时差ΔTi(x)为预先给定的走时预测模型,是关于叠加速度v、垂向双程走时t0以及偏移距x的函数.
2.2 反问题:叠加速度反演根据(5) 式所定义的正问题,反问题可以描述为:已知时距关系T(x)的具体表达形式(如双曲关系)及观测的CMP道集,反演叠加速度v与对应的t0时间(模型参数).常规的速度分析方法通过计算速度谱(如相似系数谱)并拾取其能量团的中心得到叠加速度.本文将速度分析问题同Radon变换方法相类比.为此,首先定义CMP道集所在的空间为数据空间,速度谱所在的空间为模型空间,并用广义Radon变换(沿时距关系定义的曲线积分)作为数据空间向模型空间的投影算子.此时,速度谱可以用Radon谱表征,而速度谱拾取可以转化为模型参数估计问题(此时叠加速度即Radon谱的下标:t0,v).传统的Radon谱反演可以在L2或Lp(p=0, 1) 范数意义下求解.利用Radon谱的稀疏性作为约束条件,(5) 式对应的反问题的数学表达为:
(6) |
其中,ε为噪声水平,p取0或1.
考虑到向量的L0范数表示向量中非零元素的个数,用L0范数可以最佳描述解的稀疏性.因此本文令(6) 式中p=0.通常而言,(6) 式是一个NP难(non-deterministic polynomial-time hard)的问题,可以用贪婪算法近似求解,如匹配追踪(matching pursuit,简记为MP)算法(Mallat and Zhang, 1993)及其各种变种.考虑到(5) 式中的模型参数除了反射子波w(t)及叠加速度v、垂向双程走时t0外,还包含预测误差Δti(x)及空变振幅Ai(x).因此难以直接利用传统的匹配追踪算法求解.为此,本文在标准匹配追踪方法的基础之上,引入预测误差校正的思想,给出(6) 式的一种近似求解方法.
为简化算法描述,以最简单的双曲时距关系为例:
令迭代次数k=k+1,并对数据残差进行离散双曲Radon变换,将CMP道集投影到模型空间:
(7) |
其中,Ntr为CMP道集中地震道的数目,Nt为离散时间采样点数,Nv为扫描速度的个数,mk=mk(t, v)为基于双曲Radon变换的速度谱.
(7) 式实现了从数据空间向模型空间的投影,可以用算子形式简记为:mk=Rdk-1,R为双曲Radon变换算子.由于Radon变换算子的非正交性(RHR≠I),导致模型空间中存在能量泄漏.为克服变换算子的非正交性限制,本文给出如下近似正交化方法:
首先将最优叠加速度定义为使得Radon谱能量最大的时间-速度对:
(8) |
相应的,模型空间中同相轴的投影结果wk(t)(零偏移距反射子波)可以从Radon谱中提取:
(9) |
其中,h(t; ti)为以ti为中心的窗函数,如汉明(Hamming)窗等.
基于(8) 式中得到的叠加速度以及(9) 式的子波估计,数据空间中对应的反射同相轴可以用Radon变换算子的伴随变换RH预测为:
(10) |
与(7) 式类比,(10) 式可以用算子形式简记为:uk=RHwk=RHPmk,P为由窗函数构成的滤波算子.
由于真实CMP道集中存在振幅空间变化和反射时距关系的近似表达,以及离散Radon变换对叠加速度的离散采样等因素的影响,利用伴随变换(10) 式预测的同相轴与真实CMP道集中对应的同相轴通常存在走时和振幅预测误差.根据(4) 式中定义的正问题,时差ΔTi(x)及空变振幅Ai(x)亦为模型参数.因此,本文给出如下预测误差校正方法,克服伴随Radon变换算子RH引入的预测误差.假定走时扰动Δtp小于半个子波周期,(7) 式定义的数据空间向模型空间的投影过程可以近似成立.此时,走时预测误差(即走时扰动Δtp)可用互相关函数直接测量:
(11) |
假定走时预测误差校正后的预测结果与观测数据中对应的同相轴仅相差常振幅比例因子,此时空变振幅Ai(x)可以用二者的相似系数表示为:
(12) |
最后,考虑了振幅空变以及走时预测误差后的预测结果表示为:
(13) |
若将走时扰动与振幅函数分别用算子T和A
表示,则(13) 式可以表示为:
相应的,数据残差更新为:
(14) |
(14) 式消除了稀疏反演算法选定的叠加速度(ti, vj)k对应的反射同相轴在CMP道集中的贡献,从而有助于反演逐步收敛.
同时,反演的叠加速度更新为:
(15) |
迭代执行公式(7) —(15) 可以逐步减去CMP道集中符合双曲时距关系的同相轴,直至反演终止条件即可得到离散的叠加速度S.可选的迭代终止条件有:数据残差的相对误差能量小于预先给定的阈值或反演迭代次数达到叠加速度的最大数目等.具体实现中可以在反演过程中灵活地施加单一或多种终止条件的组合.
3 数值试验 3.1 二维合成数据测试首先,我们用公式(16) 合成一个CMP道集(如图 1所示)来测试本文方法的反演精度与抗噪能力.
(16) |
其中,Ri为反射系数,w(t; t0, fm)为时延t0的Ricker子波,主频为fm,同相轴的走时采用双曲时距关系,
在反演中,设定扫描速度范围为1000~3000 m·s-1,扫描速度间隔为10 m·s-1.数据空间中第k次迭代的残差能量定义为:
为测试本文方法对噪声的容忍能力,我们将合成CMP道集加入高斯随机噪声进行测试.图 4为加入s/n=5高斯随机噪声的CMP道集(用SeismicUnix软件中suaddnoise模块实现),反演参数同第一个试验.图 5为图 4中CMP道集的双曲Radon谱.每次迭代过程中的数据残差如图 6所示.同第一个试验结果相类似,反演结束时残差能量较小.
为进一步测试随机噪声对本文方法的影响,本文将信噪比设定为s/n=1,并再次测试反演算法.然后将无噪声、s/n=5与s/n=1的反演结果用列表对比(表 2).由于本文合成的CMP道集中反射同相轴的主频依次降低,导致Radon谱的分辨率也逐渐下降(积分孔径不变).当输入数据的噪声逐渐增强时,Radon谱的信噪比也相应降低.尤其对于低频和/或能量弱的同相轴的影响更加严重.从表 2中可以看出,当s/n=1时,最后一个低主频同相轴对应的叠加速度的误差较大.
接着,我们设计了一个包含三个反射同相轴的CMP道集,如图 7a所示,第一个同相轴存在极性反转,第二个同相轴的振幅随偏移距逐渐增大,第三个同相轴为常振幅.三个同相轴的波峰振幅随偏移距变化曲线如图 8所示.图 9为图 7a中CMP道集的双曲Radon谱.即使在CMP道集存在AVO效应时,反演结果仍然较好.图 7b为反演结束时本文方法输出的预测数据,图 7c为观测数据与预测数据的残差.从图 7c中可以看出,残差相对于观测数据可以忽略不计,仅在同相轴交叉处有少量残余,从另一方面验证了本文方法的有效性.
我们将采用不同信噪比的实际资料测试本文方法的适应能力.首先,中国川西某探区的一个CMP道集如图 10所示,该CMP信噪比较高.偏移距展布范围为0.95~7.5 km.图 11为基于加权相似系数算法的常规速度谱.在反演中,取ε=0.5,并设定最大叠加次数n=50.同时,预先设置速度谱的时间精度为0.1 s(即反演的叠加速度的时间间隔最小为0.1 s).反演总共输出了11个时间-速度对,在图 11中用红色三角形表示.与常规算法得到的速度谱相比,本文方法自动反演的叠加速度恰好位于速度谱主要能量团的中心.图 11中的绿色曲线为某商业软件输出的叠加速度(人工拾取+插值).由于商业软件的处理结果已经进行了插值,生成了光滑的三维叠加速度体,因此绿色曲线不完全与本文结果相重合,但本文方法输出的叠加速度的变化趋势与商业软件处理结果基本相符.图 12为本文方法预测的CMP道集,可以看出,主要反射同相轴的运动学和动力学信息与原始数据吻合较好.图 13为反演结束时CMP道集的残差,仅含有极少量的有效信号残留.这也从另一方面证明了本文方法的稳健性.
接着,用来自中国川东北某探区的一个信噪比较低的CMP道集进一步测试本文方法.CMP道集如图 14所示.偏移距展布范围为0.19~5.9 km.图 15为基于加权相似系数算法的常规速度谱.由于数据质量不高,相似系数谱的能量团的聚焦特性较差.在反演中,设定最大叠加次数n=50,且速度谱的时间精度为0.1 s.在图 15中,红色三角形为本文的反演结果,绿色曲线为某商业软件输出的叠加速度,在T<2.5 s与本文方法结果吻合较好,而在T>2.5 s二者有一定的误差.针对这种信噪比较低,且速度谱能量团不够聚焦的数据,本文方法仍然需要进一步完善.
最后,用某浅海资料测试本文方法.CMP道集如图 16所示,偏移距展布范围为0.16~3.6 km.该CMP道集中反射同相轴较多,且无明显可识别的多次波.在反演中,设定最大叠加次数n=20,且速度谱的时间精度为0.1 s.图 17为基于加权相似系数算法的常规速度谱,红色三角形为本文的反演结果.对比速度谱可以看出,本文算法仅反演出了0.3~1.8 s区间内的叠加速度.2.0~3.0 s区间内多次波的速度谱由于能量较弱,本算法未能自动识别.若增加迭代次数,这些多次波速度谱亦能自动拾取.然而理论上,本文算法无法区分多次波与一次波的叠加速度.另一方面,在速度谱的人工拾取过程中,有经验的处理人员可以自动识别多次波速度谱从而拾取结果更加准确.
在反演精度方面,本文采用广义Radon变换实现数据空间向模型空间的投影,其速度敏感性低于基于相似系数类的方法.因此,当反射同相轴主频较低且能量较弱时(尤其在低信噪比情况下),速度谱的聚焦较差导致反演精度降低.这种情况下本文方法仍然需要进一步研究.
此外,反演终止条件的选择也是一个重要问题.由于输入数据的信噪比不同,难以给出一个普适的残差能量阈值作为反演终止条件,因而需要预先测试并估算数据信噪比.另外,可以引入额外的约束条件帮助反演收敛,如最大迭代次数.同时,还可以对反演结果进行一定的后处理,如根据模型空间中的能量以及时间分辨率约束条件等对输出结果过滤,即可输出主要能量反射同相轴对应的叠加速度,无需人工干预.
本文方法也存在一些局限性,如文中定义的正问题(CMP道集预测模型)无法考虑反射子波相位的空间变化(仅可以考虑极性反转).对于信噪比较低或波形畸变严重的实际资料,会导致速度谱的反演精度降低甚至导致算法不收敛.同时,本文方法无法区分多次波与一次反射波的速度谱,因而当CMP道集中存在明显的多次波干扰时,本文的建议是首先进行多次波压制再进行自动叠加速度估计.
5 结论本文从参数估计的角度出发,提出了一种具有明确的数学物理意义的自动叠加速度反演方法.本文方法的核心是将传统的叠加速度分析和拾取问题转化为稀疏反演框架下的模型参数估计问题,并提出了一种基于预测校正思想的匹配追踪算法自动反演模型参数(离散的叠加速度).
当地震数据的信噪比较高时,本文方法的反演精度较高且收敛性较好,可以提高叠加速度建模的效率并降低人力成本,有望实现自动化叠加速度建模并为后续的高精度速度反演方法提供较好的初始模型.
致谢感谢两位匿名审稿人提出中肯的修改建议.
Abbad B, Ursin B. 2012. High-resolution bootstrapped differential semblance. Geophysics, 77(3): U39-U47. DOI:10.1190/geo2011-0231.1 | |
Alkhalifah T, Tsvankin I. 1995. Velocity analysis for transversely isotropic media. Geophysics, 60(5): 1550-1566. DOI:10.1190/1.1443888 | |
Barros T, Lopes R, Tygel M. 2015. Implementation aspects of eigendecomposition-based high-resolution velocity spectra. Geophysical Prospecting, 63(1): 99-115. DOI:10.1111/gpr.2015.63.issue-1 | |
Beveridge J R, Ross C, Whitley D, et al. 2002. Augmented geophysical data interpretation through automated velocity picking in semblance velocity images. Machine Vision and Applications, 13(3): 141-148. DOI:10.1007/s001380100068 | |
Biondi B L, Kostov C. 1989. High-resolution velocity spectra using eigenstructure methods. Geophysics, 54(7): 832-842. DOI:10.1190/1.1442712 | |
Chen Y K, Liu T T, Chen X H. 2015. Velocity analysis using similarity-weighted semblance. Geophysics, 80(4): A75-A82. DOI:10.1190/geo2014-0618.1 | |
Davis J M. 1972. Interpretation of velocity spectra through an adaptive modeling strategy. Geophysics, 37(6): 953-962. DOI:10.1190/1.1440318 | |
Fish B C, Kusuma T. 1994. A neural network approach to automate velocity picking.//64th Annual International Meeting, SEG. Expanded Abstracts, 185-188. | |
Fomel S. 2009. Velocity analysis using AB semblance. Geophysical Prospecting, 57(3): 311-321. DOI:10.1111/gpr.2009.57.issue-3 | |
Huang K, Yang J. 2015. Seismic velocity picking using Hopfield neural network.//85th Annual International Meeting, SEG. Expanded Abstracts, 5317-5321. | |
Kirlin R L. 1992. The relationship between semblance and eigenstructure velocity estimators. Geophysics, 57(8): 1027-1033. DOI:10.1190/1.1443314 | |
Luo S, Hale D. 2012. Velocity analysis using weighted semblance. Geophysics, 77(2): U15-U22. DOI:10.1190/geo2011-0034.1 | |
Mallat S, Zhang Z. 1993. Matching pursuits with time-frequency dictionaries. IEEE Transactions on Signal Processing, 41(12): 3397-3415. DOI:10.1109/78.258082 | |
Neidell N S, Taner M T. 1971. Semblance and other coherency measures for multichannel data. Geophysics, 36(3): 482-497. DOI:10.1190/1.1440186 | |
Sacchi M D. 1998. A bootstrap procedure for high-resolution velocity analysis. Geophysics, 63(5): 1716-1725. DOI:10.1190/1.1444467 | |
Sarkar D, Baumel R T, Larner K L. 2002. Velocity analysis in the presence of amplitude variation. Geophysics, 67(5): 1664-1672. DOI:10.1190/1.1512814 | |
Schmidt J, Hadsell F A. 1992. Neural network stacking velocity picking.//62nd Annual International Meeting, SEG. Expanded Abstracts, 18-21. | |
Taner M T, Koehler F. 1969. Velocity spectra:Digital computer derivation and applications of velocity functions. Geophysics, 34(6): 859-881. DOI:10.1190/1.1440058 | |
Trad D, Ulrych T, Sacchi M. 2003. Latest views of the sparse Radon transform. Geophysics, 68(1): 386-399. DOI:10.1190/1.1543224 | |
Wang X W, Wang H Z. 2014. A research of high-resolution plane-wave decomposition based on compressed sensing. Chinese Journal of Geophysics, 57(9): 2946-2960. DOI:10.6038/cjg20140920 | |
王雄文, 王华忠. 2014. 基于压缩感知的高分辨率平面波分解方法研究. 地球物理学报, 57(9): 2946–2960. DOI:10.6038/cjg20140920 | |