本发明涉及一种地球物理勘探方法，尤其是利用地震数据时频域相位信息来构建速度模型反演的目标函数，并利用时频域模拟数据的振幅信息来对相位目标函数进行能量校正，消除反演过程中不稳定因素以及波形匹配错位等问题的局部互相关时频域相位反演方法。
背景技术：
：：全波形反演是一个最小二乘优化问题，通过不断缩小观测数据与模拟数据的残差来对事先设定初始速度模型进行更新，最终实现高精度反演的目的。与许多现有的反演方法相比，在理论上，全波形反演方法是获得地下速度模型最准确方法的算法之一。然而，即使全波形反演具有高分辨率的巨大优势，但是仍受很多不利因素的制约，如缺失低频分量、强非线性目标函数和不准确的初始模型。由于这些原因，全波形反演的目标函数常常收敛到局部极小点，并在实际数据测试中出现严重的周波跳跃问题。为了克服这些障碍，许多研究人员从以下几个方面进行了尝试：波形匹配、走时反演和相位反演。为了扩大模拟数据和观测数据之间的波形可匹配范围，已经进行了深入的研究，并提出了许多新的方法。bunks等人(1995.multiscaleseismicwaveforminversion.geophysics,60(5):1457-1473)通过滤波的方法将地震数据的高频分量与低频分量分开，实现了先反演宏观速模型构造再反演细节速度模型构造。然而，地震数据中的低频分量非常微弱，并且总是受到噪声的污染(crossholeseismicwaveformtomography–i.strategyforrealdataapplication.geophysicaljournaloftheroyalastronomicalsociety,2006,167(1):1224-1236)。更重要的是，超低频震源造价高昂，而且通常在地震勘探中不可用。因此，当地震数据中没有低频分量时，多尺度策略不能很好地恢复速度模型的宏观构造。地震包络方法(2014,seismicenvelopeinversionandmodulationsignalmodel:geophysics79wa13-wa24.)通过使用非线性希尔伯特变换能够从原始地震数据中提取低频分量。胡等人(adaptivemulti-stepfullwaveforminversionbasedonwaveformmodedecomposition.journalofappliedgeophysics,2017,139:195-210.)提出了一种波形模态分解的方法来解调出超低频分量，并且在推导伴随震源过程中不使用链式法则。这样能使得超低频分量很好地保存在了伴随源中。还有很多基于目标函数的创新方法来扩展波形可匹配的范围，如累加能量范数(accumulatedenergynormforfullwaveforminversionofmarinedata.journalofappliedgeophysics,2017,147.)，最优运输目标函数(optimaltransportforseismicfullwaveforminversion[j].2016,14(8).)，二阶时间积分波场目标函数(fullwaveforminversionofthesecond-ordertimeintegralwavefield.chinesejournalofgeophysics,2016,59(6).)和归一化积分法(thenormalizedintegrationmethod:analternativetofullwaveforminversion.2012.)。地震波走时信息与地下速度值呈线性关系(computerizedgeophysicaltomography.proceedingsoftheieee,1979,67(7):1065-1073.)。因此，走时层析成像不存在周波跳跃问题，该类方法通是常基于射线追踪理论和高频假设来反演速度模型的。然而，根据前人的研究成果显示，走时信息仅能够反映背景速度模型。当速度模型的构造复杂时，走时层析成像方法的模型分辨率有限。luo等人，(wave-equationtraveltimeinversion.geophysics,1991,56(5):645-653.)和zhou等人(acousticwave-equationtraveltimeandwaveforminversionofcrossholeseismicdata.geophysics,1995,60(3):765-773.)提出了波动方程走时反演方法来连接全波形反演和走时反演之间的鸿沟。然而，将该方法用于复杂的速度模型时，如marmousi速度模型，波动方程走时反演得到的速度模型仍不能满足常规全波形反演对初始模型要求。之后，luo等人(full-traveltimeinversion[j].geophysics,2016,81(5):r261-r274.)又提出了全走时反演方法来解决复杂模型的反演问题，并且不需要精确的初始速度模型或超低频的地震数据。胡等人(single-frequencywaveform-basedmultiscalewave-equationtraveltimeinversion[c]//segtechnicalprogramexpanded.2017:1692-1696)将时域单频波形与波动方程旅行时反演相结合，并从低频分量到高频分量逐级进行反演，最终可以得到较好的结果。地震数据中的相位特征与走时信息密切相关，其与振幅信息相比具有与模型构造更好的线性关系(theoreticalbackgroundforcontinental-andglobal-scalefull-waveforminversioninthetime–frequencydomain.geophysicaljournaloftheroyalastronomicalsociety,2008,175(2):665-685.)。因此可以通过分离相位和振幅信息来设计更多的线性反演方法。迄今为止已经提出了许多关于纯相位信息的目标函数。如：fichtner等人提出了一种时频域相位反演方法；等人(misfitfunctionsforfullwaveforminversionbasedoninstantaneousphaseandenvelopemeasurements.geophysicaljournalinternational,2011,185(2):845-870.)提出了基于希尔伯特变换的瞬时相位反演方法；choi等人(frequency-domainwaveforminversionusingthephasederivative.geophysicaljournalinternational,2013,195(3):1904-1916.)在频域或拉普拉斯傅立叶域中使用了相位信息，反演结果表明相位信息是地震反演的较好选择。同时归一化互相关目标函数(applicationofmulti‐sourcewaveforminversiontomarinestreamerdatausingtheglobalcorrelationnorm.geophysicalprospecting,2012,60(4):748-758.)与相位反演方法类似。它计算观测数据与模拟数据之间的相似性，与常规的全波形反演目标函数相比，互相关目标函数具有与地下速度模型更好的线性关系，减小了现有反演方法对初始速度模型的依赖。技术实现要素：：本发明的目的是针对上述现有技术的不足，提供一种缓解全波形反演过程中因缺失低频分量而对初始速度模型有较强的依赖性的基于加权的局部互相关时频域相位反演方法。本发明的思想是：将时频域地震数据的相位信息应用到互相关目标函数中，来缓解全波形反演方法对初始速度模型的依赖，进而缓解全波形反演过程中出现的周波跳跃问题，最终得到高精度的全波形反演结果。整个计算过程采用分炮并行策略，目的是加快整速度模型反演的进程。本发明的目的是通过以下技术方案实现的：基于加权的局部互相关时频域相位反演方法，包括以下步骤：a、将地震数据导入到计算机中，利用地震数据处理软件对其进行预处理，并准备进行加权的局部互相关时频域相位反演；本发明利用在真实速度模型上正演模拟产生的地震数据代替实际野外采集的地震数据来进行数值测试；b、建立线性递增初始速度模型，该初始速度模型作为反演的起始值；c、将观测系统、震源子波和线性递增初始速度模型输入到基于加权的局部互相关时频域相位反演算法中，并进行波动方程正演模拟，得到模拟数据，并存储正传波场；d、对观测数据和模拟数据同时进行低通滤波，滤除地震数据中高频分量。e、对观测数据和模拟数据逐道进行gabor变换得到时频域的观测数据和模拟数据；f、构建基于加权的局部互相关时频域相位反演方法的目标函数：其中表示wu表示权重因子，并设定表示观测数据，表示模拟数据在时频域的相位，表示高斯窗函数的中心时刻，表示局部角频率，并且满足g、目标函数的偏导数：并利用链式法则将目标函数的偏导数传递到时频域模拟数据对速度参数的偏导数中；h、利用gabor反变换将时频域的梯度表达式转换到时间空间域中：那么多炮数据的梯度便可以表示成：i、根据目标函数的梯度表达式提取基于加权的局部互相关时频域相位反演方法的伴随震源：并将其继续通过有限差分数值模拟的方法传播到速度模型空间，得到反传波场；j、对正传波场和反传波场利用零延迟互相关得到速度模型的更新梯度；k、利用l-bfgs优化算法计算模型更新方向：vk+1＝vk-αkhkgk，其中vk为速度模型的参数，αk为wolfe线性搜索得到的步长，hk为近似hessian矩阵的逆矩阵，gk为模型更新梯度。并通过wolfe收敛准则寻找合适的步长；l、判断是否满足终止条件，若满足则输出基于加权的局部互相关时频域相位反演方法结果。若不满足终止条件，将当前的反演结果继续作为下一次循环的初始速度模型；m、将低频段的反演结果进一步利用高频段的地震数据进行基于加权的局部互相关时频域相位反演，并输出最终反演结果。有益效果：地震数据的相位信息具有与地下速度模型较好的线性关系，利用相位信息来构建目标函数，极大程度上减弱了全波形反演对初始速度模型的依赖。因此利用时频域相位互相关目标函数，即使在初始速度模型较差或者地震数据缺失低频的情况下，最终也能得到高精度的反演结果。但是地震数据的相位变化与振幅变化无关，相位信息不能很好的量化强反射地层与弱反射地层的区别。如果不对强反射与弱反射加以区分的话，这将可能使得强反射界面对应的梯度与弱反射界面对应的梯度具有相同的数值，速度模型更新量不对称，导致整个反演流程可能会向错误的方向更新。因此本发明在相位反演的目标函数中引入权重因子，来缓解上述存在的问题。同时权重因子的引入也很好的解决了噪声对地震数据相位信息的干扰。因为当地震数据中存在微弱噪声干扰的时候，这些微弱的噪声在地震波的相位图上具有与有效信息相等的权重，这可能会使得模拟数据的相位不断向微弱噪声的相位匹配，导致整个反演向错误的方向更新，甚至得不到地下模型的任何信息。由于模拟数据中不含有噪声成分，将模拟数据作为权重因子加到时频域互相关相位目标函数中能够很好地压制噪声对相位反演的影响。附图说明图1基于加权的局部互相关时频域相位反演方法流程图。图2不同方法计算得到的伴随震源图(a)常规全波形反演的伴随震源图，(b)归一化互相关全波形反演的伴随震源图，(c)无权重因子的局部互相关时频域相位反演的伴随震源图，(d)加权的局部互相关时频域相位反演的伴随震源图。图3速度模型图(a)真实速度模型图，(b)初始速度模型图。图4常规全波形反演结果图(缺失低频测试)(a)6-15hz波形分量的反演结果图，(b)以(a)作为初始速度模型，6-40hz波形分量的反演结果图。图5归一化互相关全波形反演结果图(缺失低频测试)(a)6-15hz波形分量的反演结果图，(b)以(a)作为初始速度模型，6-40hz波形分量的反演结果图。图6加权的局部互相关时频域相位反演结果图(缺失低频测试)(a)6-15hz波形分量的反演结果图，(b)以(a)作为初始速度模型，6-40hz波形分量的反演结果图。图7不同方法反演结果的单道剖面图(源于图4,5,6最终反演结果)(a)常规全波形反演结果图，(b)归一化互相关全波形反演结果图，(c)加权的局部互相关时频域相位反演结果图。图8不同方法反演结果随着迭代次数变化的模型误差曲线图。图9不同方法计算得到的伴随震源图(强高斯噪声干扰测试)(a)常规全波形反演的伴随震源图，(b)归一化互相关全波形反演的伴随震源图，(c)无权重因子的局部互相关时频域相位反演的伴随震源图，(d)加权的局部互相关时频域相位反演的伴随震源图。图10不同方法反演结果图(强高斯噪声干扰测试)(a)常规全波形反演结果图，(b)归一化互相关全波形反演结果图，(c)加权的局部互相关时频域相位反演结果图。具体实施方式下面结合附图和实例对本发明进一步的详细说明。⑴.通过gabor变换得到时频域地震数据：其基本思想是将地震信号加上一个滑动高斯窗函数。当利用高斯窗函数来截取非稳态信号时，在这个时窗内的信号可看成是稳态的，再将加窗后的信号进行傅里叶变换，并不断地将窗函数沿着时间轴移动，依次进行傅里叶变换，该过程即为对时间域地震信号的时频分析。随着分析的进行，可得到信号的频谱在时间轴上的变化规律，称其为该信号的时频谱，通过观察时频谱的特征，更能直观的反映信号频率随时间变化的规律。⑵.在时频域构建互相关相位目标函数：地震数据相位反演实际上是不断地利用模拟数据的相位信息去匹配观测数据的相位信息，本发明利用地震数据的时频域相位信息，并结合互相关目标函数来构建一个时频域目标函数。互相关相位目标函数可以更好的衡量两个数据的相位相似程度，当相位相似程度增加时证明速度模型更接近真实速度模型，当相位相似程度减小时证明反演结果向错误的方向更新迭代。通过判断目标函数增大或者减小就可以控制局部优化算法向正确的方向更新迭代，使得初始速度模型逐渐向真实速度模型靠近，最终实现高精度反演的结果。⑶.对时频域相位目标函数加入权重约束：地震数据的相位变化与振幅变化无关，这一点使得地震数据中的相位信息在强反射界面和弱反射界面具有相同的值，这将可能使得强反射界面对应的梯度与弱反射界面对应的梯度具有相同的数值，使得速度模型更新量不对称，导致整个反演流程可能会向错误的方向更新。同时当地震数据中存在微弱的噪声干扰的时候，这些微弱的噪声在地震数据的相位图上具有与有效信息相等的权重，这可能会使得模拟数据的相位不断向微弱噪声的相位匹配，导致整个反演向错误的方向进行，甚至根本得不到地下速度模型的任何有效信息。因此在设计相位目标函数的同时需要在其中加上一个权重项，来避免上述问题的产生。使得目标函数在反演的过程中更加稳定。⑷.利用链式法则推导目标函数的梯度：为了求取目标函数对速度参数的偏导数，需要在时频域利用链式法则来得到模拟数据对速度参数的偏导数，这样才能结合伴随状态法来得到地下速度模型的更新梯度。⑸.结合gabor反变换获得时间域伴随震源：利用链式法则得到了时频域模拟数据关于速度参数的偏导数，为了进一步与时间域有限差分正演模拟相结合，需要利用gabor反变换将时频域的地震数据转换到时间域。这样便可以得到与常规全波形反演类似的伴随震源，因此该方法具有与常规全波形反演相似的实施流程。⑹.利用伴随震源反传得到反传波场，并与正传波场进行零延迟互相关得到目标函数的更新梯度：最原始求取梯度的方法是利用扰动理论来求解的，通过对地下模型每个点增加微小的扰动后，再进行正演模拟得到相应的模拟数据。该扰动理论没有利用链式法则的一阶导数近似，但其计算量巨大。以本文以测试模型(69×192)为例，每次迭代需要正演模拟26496次，一次正演模拟约15秒，那么一次模型更新需要110.4小时。当然这还不包括优化算法步长选取不合适需要重新计算的情况。利用伴随状态法以后，迭代一次需要4次正演模拟约1分钟的时间。所以通过伴随状态法求取目标函数的梯度，大大缩短了全波形反演梯度计算所需要的时间。⑺.利用l-bfgs优化算法来获得速度模型的更新方向：l-bfgs优化算法有效克服了显式求取hessian矩阵及其逆的困难，其具有超线性收敛速度，计算效率高，占用内存小，精度高等优点，较适合求解大规模非线性优化问题。基于加权的局部互相关时频域相位反演方法，包括以下步骤：a、安装matlab2013a正版软件，并安装基于matlab语言编写的地震数据处理软件crews工具包。同时将地震数据导入到计算机中，利用地震数据处理软件对其进行预处理，并准备进行加权的局部互相关时频域相位反演；本发明利用在真实速度模型上正演模拟产生的地震数据代替实际野外采集的地震数据来进行数值测试；b、建立线性递增初始速度模型，该初始速度模型作为反演的起始值，如图3(b)所示。线性递增初始速度模型基本不含有真实速度模型(图3a)任何构造信息，这样能够很好地验证本发明方法的稳定性和鲁棒性；c、将观测系统、震源子波和线性初始速度模型输入到基于加权的局部互相关时频域相位反演算法中，并进行波动方程正演模拟，得到模拟数据，并存储正传波场。其中常密度二维二阶声波方程可以表示为：其中x,z表示空间坐标轴，v表示速度模型的参数，p表示声波传播的压力场。通过有限差分方法来对上述声波波动方程进行求解；d、对观测数据和模拟数据同时进行低通滤波，滤除地震数据中高频分量：由于全波形反演方法是一个强非线性优化问题，并且通常无法搜索到全局最小值，导致严重周波跳跃现象。所以直接将全波形反演方法应用于地震数据反演是不现实的。本发明使用巴特沃斯滤波器将地震数据分为不同的频率尺度分别进行全波形反演。当利用较低频段的的地震数据进行反演的时候，局部极小值较少，可以利用低频段大尺度的地震数据得到一个较好的速度模型，再进行高频段的反演。基于这个理论，首先使用低频带(6-15hz分量)的地震数据来反演一个宏观速度模型，然后将其作为较高频带(6-40hz分量)的初始速度模型进一步获得最终的反演结果；e、对观测数据和模拟数据逐道进行gabor变换，得到时频域的观测数据和模拟数据。其中地震观测数据和模拟数据的傅里叶变换可以表示为：地震观测数据和模拟数据gabor变换可以表示为：其中表示高斯窗函数的中心时刻；表示局部角频率，d(t)表示时间域观测数据，表示频率域观测数据，表示时频域观测数据，u(t)表示时间域模拟数据，表示频率域模拟数据，表示时频域模拟数据。f[·]表示傅里叶变换，g[·]表示gabor变换；f、构建基于加权的局部互相关时频域相位反演方法的目标函数：其中表示wu表示权重因子，并设定和分别表示模拟数据和观测数据在时频域的相位，并且满足g、关于目标函数对速度参数的偏导数，可以并利用链式法则将目标函数的偏导数传递到模拟数据对速度参数的偏导数：其中有并且带入上式：h、利用gabor反变换将时频域的梯度表达式转换到时间域中，其中地震观测数据和模拟数据傅里叶逆变换可以表示为：地震观测数据和模拟数据gabor逆变换可以表示为：且存在共轭复数的gabor变换：将上式带入得到：同时根据gabor反变换公式得到最终的梯度公式：i、求取伴随震源(如图2所示)：并将其继续通过有限差分数值模拟的方法传播到速度模型空间，得到反传波场；j、对正传波场和反传波场利用零延迟互相关的方法得到速度模型的更新梯度，对应多炮地震数据的梯度可以表示为：其中pf表示正传波场,表示反传波场；k、利用l-bfgs优化算法计算模型更新方向：vk+1＝vk-αkhkgk，其中vk为速度模型的参数，αk为wolfe线性搜索得到的步长，hk为近似hessian矩阵的逆矩阵，gk为模型更新梯度。并通过wolfe收敛准则寻找合适的步长；l、判断是否满足终止条件(如图8所示)，若满足则输出基于加权的局部互相关时频域相位反演方法结果。若不满足终止条件，将当前的反演结果继续作为下一次循环的初始速度模型；m、将低频段的反演结果进一步利用高频段的地震数据进行基于加权的局部互相关时频域相位反演，并输出最终反演结果(如图6所示)。实施例1根据前面的数值试验，marmousi模型的左侧是全波形反演方法最难恢复的构造。因此，我们截取marmousi模型左侧部分来对本发明提出的方法进行测试，并与其他两种方法进行对比(常规全波形反演方法和归一化归一化互相关全波形反演方法)。初始速度模型如图3b所示，共有36个震源在表面上均匀分布，每个震源对应着192个检波器。地震数据记录时间为2s，时间间隔为1ms。同时为了证明基于加权的局部互相关时频域相位反演方法的有效性，在测试中利用巴特沃斯高通滤波器切除雷克子波6hz以下的分量。在数值测试中，我们使用l-bfgs优化算法来计算下降方向来对初始模型进行更新。模型参数如下：表1：基于加权的局部互相关时频域相位反演方法测试参数模型大小网格间距横向距离纵向深度速度范围地震子波主频69*19212.5m2.4km0.8625km1.5-4km/s15hz基于以上参数，在表2所示的测试环境下，进行基于加权的局部互相关时频域相位反演方法测试。表2计算机性能测试环境将图4(a)与图5(a)对比可以看出，归一化互相关全波形反演结果优于常规的全波形反演的结果。在图4(a)中的marmousi模型的顶端有一组异常值，这是由常规全波形反演方法在缺失低频分量的时候会出现明显的周波跳跃问题。由于归一化互相关全波形反演方法使用了互相关目标函数和归一化地震数据来反演速度模型，其主要测量观测数据和模拟数据之间的相似性。因此，它比常规的全波形反演具有与地下速度模型更强的线性关系。从图6a和图7中可以看出，加权的局部互相关时频域相位反演结果有了很大的改善，特别是在marmousi左侧的水平构造。由于相位互相关函数具有与地下速度模型更好的线性关系，极大程度上减弱了全波形反演对初始速度模型的依赖。此外，在时频域目标函数中引入权重因子，不仅增强了反演过程的稳定性，而且对缓解周波跳跃问题有着巨大的优势。数值试验结果表明，基于加权的局部互相关时频域相位反演方法能够有效缓解对初始速度模型的依赖。当地震记录数据中没有低频分量时，常规的全波形结果存在严重的周波跳跃问题(图4)，反演模型逐渐偏离真实模型。如图8所示，常规全波形反演结果的拟合差始终在不断增加。然而，随着迭代次数增加，基于加权的局部互相关时频域相位反演方法的模型拟合差在逐渐减小。证了基于加权的局部互相关时频域相位反演方法具有不依赖初始速度模型的优势，并且随着不断提高地震数据的主频分量，最终可以得到高分辨率的反演结果。从图10a可以看出，常规全波形反演结果与真实速度模型相差甚远，这主要是由于缺少低频分量，并且受到强高斯背景噪声的严重干扰。将图10b与图10a进行比较发现，归一化互相关全波形反演结果比常规全波形反演结果稍好一些，主要是因为互相关目标函数具有抵抗高斯噪声的能力。然而，由于信噪比太低，大部分有效波信息都被淹没在强背景噪声中，导致最终的反演结果仍然偏离真实速度模型。结合图9d和图10c可知，强高斯背景噪声对加权的局部互相关时频域相位反演方法几乎没有任何影响。证明了，基于加权的局部互相关时频域相位反演方法能够衰减伴随源中的强高斯背景噪声，以克服低信噪比带来的周期跳跃现象。当前第1页12