技术领域
[0001] 本
发明涉及
地震数据处理技术领域，特别涉及一种三维无射线追踪回折波层析成像方法及装置。
背景技术
[0002] 在地震数据处理中，旅行时层析成像方法是一种重要的速度成像方法。旅行时层析成像可以按照计算流程分为正演旅行时、反演慢度两个主要步骤。正演旅行时包括旅行时计算，通常也含有射线追踪。依据正演旅行时后是否进行射线追踪，可以简单将旅行时层析成像划分为两类情形：
[0003] (1)传统的射线追踪旅行时层析成像方法。本类方法设定速度分布在离散网格
节点上，节点间的面元速度为常数或渐变。对旅行时正演后的网格可以用诸如梯度法、互换法求取射线路径，从而得到每个接收点正演旅行时和拾取初至的误差，然后进行
迭代反演更新速度场。
[0004] 现有的三维回折波层析成像主要采用传统射线追踪旅行时层析成像方法，此时约束速度横向不变，垂向渐变递增，从而满足回折波传播的条件，然后迭代进行三维射线追踪和反演以更新慢度。
[0005] 本类方法优点是对地下速度分布无任何假定，可以适应任意复杂或简单的速度场。缺点是复杂的速度场对射线路径的计算方法有较高要求，需要处理旅行时多至、射线焦散、速度梯度不连续等复杂问题。另外，射线路径一般为曲线，使得反演呈现强非线性，求解困难并存在多解性。当约定速度仅仅沿垂向渐变时，波沿垂向渐变递增的速度场以回折
波形式传播并回折至地表，得到传统的回折波旅行时层析成像方法。
[0006] (2)无射线追踪回折波旅行时层析成像方法。本类方法同样约定速度仅仅沿垂向渐变。但是，依据回折波的传播规律，可以推导公式解析计算回折波回折点的深度和速度，从而避免复杂的正演旅行时和射线追踪、反演更新慢度过程。
[0007] 现有的三维回折波层析成像主要采用无射线追踪回折波旅行时层析成像方法，该方法本质上是一种二维
算法，无法进行三维速度场成像。
[0008] 综上，现有方法是对三维地震数据的初至时间进行三维分解(Osypov,2000)，然后对分解后的旅行时应用二维算法进行计算。但是分解的算法复杂且需要多次迭代。
发明内容
[0009] 为了解决
现有技术的问题，本发明提出一种三维无射线追踪回折波层析成像方法及装置，扩展二维无射线追踪回折波层析成像技术至三维地震资料处理，避免了复杂的三维旅行时迭代分解过程。
[0010] 为实现上述目的，本发明提供了一种三维无射线追踪回折波层析成像方法，该方法包括：
[0011] 拾取单炮初至，对拾取的每个单炮初至按照方位
角度间隔
选定至少两个方位角的初至，并对每一个方位角初至按照偏移距递增方式排序，形成0～x范围内的初至曲线；其中，x为最大偏移距；
[0012] 以偏移距为横坐标、初至时间为纵坐标，对0～x范围内的初至曲线进行低阶多项式最小二乘法整体拟合，得到每一单炮初至的方位角对应地多项式系数；
[0013] 根据所述多项式系数获得对应地每一单炮初至的方位角慢度曲线；
[0014] 根据所述每一单炮初至的方位角慢度曲线获得拾取的每个单炮初至所选定的方位角对应地回折波回折点的速度和深度；
[0015] 根据拾取的每个单炮初至所选定的方位角对应地回折波回折点的速度和深度插值成三维速度场。
[0016] 优选地，所述拾取的每个单炮初至所选定的方位角对应地回折波回折点的速度获取方法为：
[0017] 对拾取的每一单炮初至的方位角慢度曲线在0～x范围内任取M点为回折波回折点；
[0018] 求出点i对应
位置的慢度曲线值的倒数，得到点i的速度；其中，i＝1,2，…，M。
[0019] 优选地，所述拾取的每个单炮初至所选定的方位角对应地回折波回折点的深度获取方法为：
[0020] 对点i的深度采用积分公式获得；其中，所述积分公式为：
[0021]
[0022] 其中，zi为回折点i对应的深度，xi为点i对应的偏移距。
[0023] 优选地，所述每一单炮初至的方位角慢度曲线获取的方法为：
[0024] 根据所述多项式系数获得多项式；
[0025] 对所述多项式进行求导，获得每一单炮初至的方位角慢度曲线。
[0026] 优选地，所述方位角度间隔为相当或不相等的方位角度间隔。
[0027] 为实现上述目的，本发明还提供了一种三维无射线追踪回折波层析成像装置，该装置包括：
[0028] 初至曲线获取单元，用于拾取单炮初至，对拾取的每个单炮初至按照方位角度间隔选定至少两个方位角的初至，并对每一个方位角初至按照偏移距递增方式排序，形成0～x范围内的初至曲线；其中，x为最大偏移距；
[0029] 拟合单元，用于以偏移距为横坐标、初至时间为纵坐标，对0～x范围内的初至曲线进行低阶多项式最小二乘法整体拟合，得到每一单炮初至的方位角对应地多项式系数；
[0030] 方位角慢度曲线获取单元，用于根据所述多项式系数获得对应地每一单炮初至的方位角慢度曲线；
[0031] 回折波回折点的速度和深度获取单元，用于根据所述每一单炮初至的方位角慢度曲线获得拾取的每个单炮初至所选定的方位角对应地回折波回折点的速度和深度；
[0032] 插值单元，用于根据拾取的每个单炮初至所选定的方位角对应地回折波回折点的速度和深度插值成三维速度场。
[0033] 优选地，所述回折波回折点的速度和深度获取单元包括：回折波回折点获取模
块和第一计算模块；其中，
[0034] 所述回折波回折点获取模块，用于对拾取的每一单炮初至的方位角慢度曲线在0～x范围内任取M点为回折波回折点；
[0035] 所述第一计算模块，用于求出点i对应位置的慢度曲线值的倒数，得到点i的速度；其中，i＝1,2，…，M。
[0036] 优选地，所述回折波回折点的速度和深度获取单元还包括第二计算模块；其中，[0037] 所述第二计算模块，用于对点i的深度采用积分公式获得；其中，所述积分公式为：
[0038]
[0039] 其中，zi为回折点i对应的深度，xi为点i对应的偏移距。
[0040] 优选地，所述方位角慢度曲线获取单元包括：多项式获取模块和求导模块；其中，[0041] 所述多项式获取模块，用于根据所述多项式系数获得多项式；
[0042] 所述求导模块，用于对所述多项式进行求导，获得每一单炮初至的方位角慢度曲线。
[0043] 优选地，所述初至曲线获取单元采用的方位角度间隔为相当或不相等的方位角度间隔。
[0044] 上述技术方案具有如下有益效果：
[0045] 将传统无射线追踪的回折波旅行时二维层析成像方法沿着三维地震数据的多个方位角进行分方位角应用，从而得到三维计算结果。本技术方案避免了复杂的三维旅行时迭代分解过程，因此易于实现，当选择方位角个数较少，例如少于40个，计算速度也很快。
附图说明
[0046] 为了更清楚地说明本发明
实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0047] 图1为本发明提出的一种三维无射线追踪回折波层析成像方法
流程图；
[0048] 图2为本发明提出的一种三维无射线追踪回折波层析成像装置
框图；
[0049] 图3为本实施例的三维无射线追踪回折波层析成像方法流程图；
[0050] 图4为本实施的层析成像效果图。
具体实施方式
[0051] 下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0052] 本技术方案的工作原理是：对二维的单炮初至进行
抽取，获得多个方位角的二维线，运用二维无射线追踪回折波旅行时层析成像算法进行计算，依据多方位角的组合达到三维计算的效果，避免复杂的迭代三维旅行时分解计算。
[0053] 如图1所示，为本发明提出的一种三维无射线追踪回折波层析成像方法流程图。该方法包括：
[0054] 步骤101)拾取单炮初至，对拾取的每个单炮初至按照方位角度间隔选定至少两个方位角的初至，并对每一个方位角初至按照偏移距递增方式排序，形成0～x范围内的初至曲线；其中，x为最大偏移距；
[0055] 在步骤101)中，所述单炮初至拾取由处理
软件通过自动或者交互拾取的方式获得，获得的初至数据中必须含有本道的方位角信息。对拾取的单炮初至数据，按照相等或不相等的方位角度间隔，提取多个方位角的初至，并对每一个方位角的初至按偏移距递增排序，形成0～x范围内初至曲线。
[0056] 单炮初至的各个方位角，可以是一个，也可以是多个方位角。只有一个方位角时，等同二维层析计算的效果。当规定多个方位角时，可以按相等或不相等的方位角度间隔，灵活选择需要的角度。
[0057] 步骤102)以偏移距为横坐标、初至时间为纵坐标，对0～x范围内的初至曲线进行低阶多项式最小二乘法整体拟合，得到每一单炮初至的方位角对应地多项式系数；
[0058] 在步骤102)中，所述初至曲线整体拟合基于选定的单方位角的初至曲线：以偏移距为横坐标，初至时间为纵坐标，对0～x范围的初至进行低阶多项式最小二乘法整体拟合。拟合时选择能用多项式根加权求和进行快速积分的多项式，如勒让德多项式，契比
雪夫多项式。多项式用5阶以下的低阶，通常选择2阶。拟合完毕后得到初至数据对应的多项式系数。
[0059] 步骤103)根据所述多项式系数获得对应地每一单炮初至的方位角慢度曲线；
[0060] 在步骤103)中，步骤102)中的多项式系数获得多项式，所述的慢度曲线由对多项式求导数得到。计算方法是：用步骤102)中计算的多项式的导数公式，计算0～x范围内导数p＝dt/dx的值，即可得到慢度曲线。其中，t为初至时间，x为偏移距，d为求导运算符号。由于慢度曲线计算直接使用了多项式导数公式，因此，计算采用的坐标可以是规则排列的，也可以是不规则排列的。
[0061] 步骤104)根据所述每一单炮初至的方位角慢度曲线获得拾取的每个单炮初至所选定的方位角对应地回折波回折点的速度和深度；
[0062] 在步骤104)中，所述的速度计算方法是：对步骤103)求得的
指定方位角慢度曲线，在0～x范围内，任取M个点，求出点i位置的慢度曲线值的倒数1/(dt/dx)i(i＝1,2,…,M)，得到该点回折点的速度。
[0063] 在步骤104)中，点i的深度采用如下积分公式：
[0064]
[0065] 其中，zi为回折点i对应的深度，xi为偏移距。深度计算是在0～xi范围内，逐点计算积分核函数，然后，对0～xi积分得到zi。由于步骤103)中多项式选择了能用多项式根加权求和进行快速积分的多项式，因此积分采用高斯积分进行快速积分，对于M个积分计算点，可以获得2M-1阶多项式的精确解。由于步骤102)中拟合初至曲线采用了低阶多项式，如2阶，因此M可以取得很小，如M＝2，即可获得3阶多项式精确解。所以本发明在回折点深度计算时，不但积分
精度高，而且积分计算量很小，计算效率高。
[0066] 在步骤104)计算时，选取的M点坐标可以等间距，也可以不等间距。M个点可以
覆盖全部偏移距范围，也可以只覆盖一部分偏移距范围。
[0067] 步骤105)根据拾取的每个单炮初至所选定的方位角对应地回折波回折点的速度和深度插值成三维速度场。
[0068] 如图2所示，为本发明提出的一种三维无射线追踪回折波层析成像装置框图。该装置包括：
[0069] 初至曲线获取单元201，用于拾取单炮初至，对拾取的每个单炮初至按照方位角度间隔选定至少两个方位角的初至，并对每一个方位角初至按照偏移距递增方式排序，形成0～x范围内的初至曲线；其中，x为最大偏移距；
[0070] 拟合单元202，用于以偏移距为横坐标、初至时间为纵坐标，对0～x范围内的初至曲线进行低阶多项式最小二乘法整体拟合，得到每一单炮初至的方位角对应地多项式系数；
[0071] 方位角慢度曲线获取单元203，用于根据所述多项式系数获得对应地每一单炮初至的方位角慢度曲线；
[0072] 回折波回折点的速度和深度获取单元204，用于根据所述每一单炮初至的方位角慢度曲线获得拾取的每个单炮初至所选定的方位角对应地回折波回折点的速度和深度；
[0073] 插值单元205，用于根据拾取的每个单炮初至所选定的方位角对应地回折波回折点的速度和深度插值成三维速度场。
[0074] 优选地，所述回折波回折点的速度和深度获取单元204包括：回折波回折点获取模块和第一计算模块；其中，
[0075] 所述回折波回折点获取模块，用于对拾取的每一单炮初至的方位角慢度曲线在0～x范围内任取M点为回折波回折点；
[0076] 所述第一计算模块，用于求出点i对应位置的慢度曲线值的倒数，得到点i的速度；其中，i＝1,2，…，M。
[0077] 优选地，所述回折波回折点的速度和深度获取单元204还包括第二计算模块；其中，
[0078] 所述第二计算模块，用于对点i的深度采用积分公式获得；其中，所述积分公式为：
[0079]
[0080] 其中，zi为回折点i对应的深度，xi为点i对应的偏移距。
[0081] 优选地，所述方位角慢度曲线获取单元203包括：多项式获取模块和求导模块；其中，
[0082] 所述多项式获取模块，用于根据所述多项式系数获得多项式；
[0083] 所述求导模块，用于对所述多项式进行求导，获得每一单炮初至的方位角慢度曲线。
[0084] 优选地，所述初至曲线获取单元201采用的方位角度间隔为相当或不相等的方位角度间隔。
[0085] 为了深入了解本发明，下面结合附图和应用实例对本发明的实施例作进一步详细描述。
[0086] 如图3所示，为本实施例的三维无射线追踪回折波层析成像方法流程图。具体实施方式为：
[0087] 步骤1)拾取单炮初至并按方位角分选：首先拾取初至数据，并依据工区范围和炮、检点地表高程，建立网格速度模型，初始化速度为0。对拾取的单炮初至，按方位角分选，并对分选后的初至按偏移距递增排序，形成0～3000米范围内初至曲线，其中3000米为选择的最大偏移距。
[0088] 步骤2)对指定方位角的初至曲线进行整体拟合：在单炮内，对指定方位角的初至，采用2阶勒让德多项式，以偏移距为横坐标，初至时间为纵坐标，进行最小二乘法整体拟合，解出本方位角对应的多项式系数。
[0089] 步骤3)对指定方位角用拟合的多项式曲线计算慢度曲线：利用步骤2)中得到的多项式系数，依据2阶勒让德多项式导数公式，以步骤2)中的偏移距为横坐标，计算0～3000米范围内慢度曲线p＝dt/dx。
[0090] 步骤4)计算该方位角的回折波回折点的对应速度和深度：首先，在0～3000范围内，逐点计算步骤3)得到的慢度曲线p的倒数1/(dt/dx)i(i＝1,2,…,M)，得到第i点对应的回折点的速度；其次，采用3点高斯-勒让德积分，计算点i的回折点深度。利用深度计算公式：
[0091]
[0092] 首先对0～i点逐点取导数值并除以常数(dt/dx)i，计算cosh-1函数值；其次，查数学手册得到3点高斯-勒让德多项式根和权系数，在0～xi范围加权求和求得精确积分，计算出回折点深度zi。
[0093] 高斯-勒让德积分可以对等距节点进行积分，也可以对不等距节点进行积分，因此取原初至对应的偏移距为计算坐标点，可以适应检波点规则排列和非规则排列的情形。计算时如只需要对部分偏移距计算，则限制积分范围即可。
[0094] 步骤5)对单炮的各个方位角，重复步骤2)～步骤4)，得到相应的回折点深度和对应速度：此时选择36个方位角，以10度相等间隔分布。
[0095] 步骤6)对所有单炮重复步骤2)～步骤5)，得到全部回折点深度和对应速度。
[0096] 步骤7)将速度点插值成三维速度场。
[0097] 如图4所示，为本实施的层析成像效果图。对于图4展示的山地资料，本发明得到的三维速度场速度自地表向深层逐渐变化，在空间上向各个方向展布均匀，变化自然，无明显异常体出现，且在地表附近也比较均匀，没有出现常规射线层析由于地表变化剧烈引起的计算结果不均匀、速度分布异常。另外，本发明在计算图4结果时，仅用集群上1个节点1个线程，耗时3分钟，而相同集群上常规射线层析8节点并行，每节点采用10线程，耗时30分钟才得到相近结果。因此本发明的效率优势非常明显。
[0098] 以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何
修改、等同替换、改进等，均应包含在本发明的保护范围之内。
