本发明涉及地震勘探领域，尤其涉及一种基于两点射线追踪的地震走时层析反演方法。
背景技术：
地震层析成像是指利用地震观测数据反演研究区域速度结构的一种方法。地震层析成像的原理类似于医学ct技术，是根据弹性波理论及其在地层介质中的传播规律，对观测得到的弹性波在岩土体介质中的走时或波形进行反演计算，重建被测范围内岩体弹性波速度分布规律的图像，从而达到确定地层结构或圈定地质异常体的一种物探反演解释方法。基于射线理论的射线追踪法是正演算法的一种，射线追踪法基于高频近似，只计算源点到接收点的射线路径，计算效率高。反演计算一般是采用最优化算法，例如最速下降法、共轭梯度法、牛顿迭代法、随机搜索等，求解满足给定最小误差准则的最优解。
当地层速度只在深度方向上变化，水平方向上无变化时，地震层析成像可用一维速度模型描述研究区域的速度结构，即速度只是深度的函数。在近地表浅层的层析成像问题中，小区域范围内的地层速度变化一般可用一维模型近似，在地震数据采集密度不高的情况下，工程物探一般会采用基于一维模型的地震层析成像方法。现有的一维走时层析速度反演假设每层层内速度都是均匀的，通过多次迭代来得到一个层状的速度模型，一般需要划分比较多的层数才能比较准确地描述实际地层速度变化的特征。实际地层的速度变化多为非均匀且呈现出连续变化的特征，因此利用现有技术得到的结果存在固有近似性。并且层数越多需要反演的模型参数就越多，因而地震层析成像的计算量越大。另外，在大入射角的情况下(一般震源与接收点的距离比地震信号反射深度大很多的时候会出现此情况)，利用现有技术中的两点走时射线追踪牛顿迭代方法存在收敛慢或不收敛的问题。
技术实现要素：
本发明的目的在于提供一种基于两点射线追踪的地震走时层析反演方法，旨在解决现有技术中地震层析成像的计算量繁重和现有的两点走时射线追踪牛顿迭代方法存在的收敛慢或不收敛的问题。
本发明是这样实现的，一种基于两点射线追踪的地震走时层析反演方法，包括下述步骤:
在研究区域进行地震数据采集，获取直达波走时数据和反射波走时数据；
对研究区域进行模型参数化，建立层内速度可连续变化的初始一维连续层状模型；
根据所述层内速度可连续变化的一维连续层状模型，将射线参数p用变量q表示，震源距x表示为变量q的函数x＝f(q)，利用牛顿迭代法求解x＝f(q)，进而得到射线参数p，射线路径由射线参数p唯一确定，得到射线参数p后，计算得到理论的直达波走时和反射波走时；
比较计算得到的理论直达波走时和反射波走时与所述地震数据采集获得的实际的直达波走时数据和反射波走时数据，判断理论直达波走时和反射波走时与所述地震数据采集获得的实际的直达波走时数据和反射波走时数据的差别是否满足给定的误差标准，是则进行模型输出，否则进行下一步；
用最优化算法调整所述层内速度可连续变化的一维连续层状模型，直到计算得到的理论的直达波走时和反射波走时和实际的直达波走时数据和反射波走时数据的差别满足给定的误差标准为止，并进行模型输出。
本发明提供的基于两点射线追踪的地震走时层析反演方法与现有技术相比存在以下有益效果：通过建立地层速度可连续变化的一维连续层状模型，使基于此模型的反演变量数可以大幅度减少，可以更加准确地描述实际地层速度结构，显著提高反演计算效率；通过将射线参数p用变量q表示，间接利用变量q求解射线参数p，使得迭代求解过程稳定，收敛快速，有效避免在大入射角的情况下迭代不收敛的问题。
附图说明
图1是本发明实施例提供的基于两点射线追踪的地震走时层析反演方法流程图；
图2是本发明实施例提供的基于两点射线追踪的地震走时层析反演方法连续层状模型的参数定义；
图3是本发明实施例提供的基于两点射线追踪的地震走时层析反演方法地层模型示意图；
图4是本发明实施例提供的基于两点射线追踪的地震走时层析反演方法地震勘探地表数据采集示意图；
图5是本发明实施例提供的基于两点射线追踪的地震走时层析反演方法地表激发垂直地震剖面示意图；
图6是本发明实施例提供的基于两点射线追踪的地震走时层析反演方法井中激发垂直地震剖面示意图；
图7是本发明实施例提供的基于两点射线追踪的地震走时层析反演方法井间成像示意图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
请参阅图1至图3，本发明实施例提供了一种基于两点射线追踪的地震走时层析反演方法，该方法包括下述步骤：
步骤s1：采集实际走时数据。在研究区域进行地震数据采集，获取直达波走时数据和反射波走时数据；
步骤s2：模型参数化。对研究区域进行模型参数化，建立层内速度可连续变化的初始一维连续层状模型；
具体地，该层内速度可以是纵波或横波速度，随深度递增或递减；
步骤s3：计算射线参数。根据层内速度可连续变化的一维连续层状模型，将射线参数p用变量q表示，震源距x表示为变量q的函数x＝f(q)，利用牛顿迭代法求解x＝f(q)，进而得到射线参数p，射线路径由射线参数p唯一确定；
步骤s4：计算理论走时。根据射线参数p计算直达波走时和反射波走时；
步骤s5：比较理论走时与实际走时。比较计算得到的理论直达波走时和反射波走时与所述地震数据采集获得的实际的直达波走时数据和反射波走时数据，判断理论直达波走时和反射波走时与所述地震数据采集获得的实际的直达波走时数据和反射波走时数据的差别是否满足给定的误差标准，是则进行步骤s7，否则进行步骤s6；
步骤s6：优化模型。用最优化算法调整层内速度可连续变化的一维连续层状模型，直到计算得到的理论的直达波走时和反射波走时和实际的直达波走时数据和反射波走时数据的差别满足给定的误差标准为止。
具体地，使用的最优化算法可以为最速下降法、共轭梯度法、牛顿迭代法、随机搜索等方法；
步骤s7：模型输出。计算得到的理论到时和实际数据到时的差别满足给定的误差标准后，进行模型输出。
本发明实施例提供的基于两点射线追踪的地震走时层析反演方法建立了层内速度可连续变化的一维连续层状模型，将层内速度表示为深度的函数，此模型允许层内速度连续变化，大大减少了划分层数，从而大幅度减少了反演变量数，可以更加准确地描述实际地层速度结构，显著提高反演计算效率；并且通过将射线参数p用变量q表示，间接利用变量q求解射线参数p，使得迭代求解过程稳定，收敛快速，有效避免了在大入射角的情况下迭代不收敛的问题。
进一步地，请参考图2，在层内速度可连续变化的一维连续层状模型中，层内速度vk为深度z的函数，第k层速度函数表示为：
vk＝akz+bk，
其中，下标k表示第k层，ak和bk是需要反演的模型参数，分别为第k层速度函数的梯度和截距，当反演地层纵波速度模型时，层内速度vk为纵波速度；当反演地层横波速度模型时，则层内速度vk为横波速度；当反演地层转换波速度模型时，则层内速度vk由层内转换波属性决定其为纵波速度或横波速度，这里的转换波为纵波到横波转换波或横波到纵波转换波，该层内转换波属性为纵波属性或横波属性。
当ak＝0时，表示在第k层为层内速度不变的均匀层。
进一步地，根据斯涅尔定律，震源距x表示为射线参数p的函数，震源距x表示为：
其中，各项定义为：
其中，εk，ωk，hk，μk和δk为中间参数，下标s表示震源所在的层，可取值范围为1至n，n为反射波反射发生所在层的标号，zs为震源深度，z(k)表示第k层的深度，下标k表示第k层，k＝0项是关于震源位置的一个修正项。
进一步地，为了避免在大入射角情况下迭代不收敛的问题，将射线参数p用变量q表示为：
进一步地，将震源距x表示为变量q的函数，震源距x表示为：
其中，各项定义为：
给定震源距x，利用牛顿迭代法求解x＝f(q)，得到参数q的值，将参数q代回到与射线参数p的关系式中，即可得到射线参数p的值。
进一步地，牛顿迭代法求解方程x＝f(q)时，q的初值通过以下方法获得：将震源距x和变量q用有理函数公式逼近，有理函数公式为：
其中，系数α1，α2，β1和β2通过有理函数公式和震源x的函数公式的泰勒展开式获得，利用上述有理函数公式，获得q的初值估算公式：
其中，各项系数α1，α2，β1和β2的表达式分别为：
利用所述q的初值估算公式获得初值进行迭代计算，经过迭代得到q的精确解。利用此方法获得的初值精度一般能达到95％以上，只需要1至3次迭代就可以得到q的精确解。
进一步地，得到射线参数p的值后，直达波走时和反射波走时的计算公式为：
进一步地，请参考图4至图7，在研究区域进行地震数据采集的步骤具体为：
如图4所示的地表地震勘探，震源和接收检波器都布置在地表；
或者为：如图5所示的垂直地震剖面勘探，震源在地表，接收器在井中；
或者为：如图6所示的垂直地震剖面勘探，震源在井中，接收器在地表；
或者为：如图7所示的井间成像勘探，震源和接收检波器在两口不同的井中。
进一步地，地震数据采集采用的地震信号是纵波、或横波、或纵波到横波的转换波，或横波到纵波的转换波。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。