井间及井地联合地震波CT成像方法及装置与流程

本发明涉及地震勘探及开发领域，更具体地，涉及一种井间及井地联合地震波ct成像方法。

背景技术：

地震层析成像是一种基于模型的数据拟合方法。与医学中的计算机辅助层析成像(computertomography，ct)类似并受其启发，地震层析成像通过分析台站记录到的由天然地震或者人工地震激发的地震波信号来获取地球内部结构的信息。

1976年，aki和lee从60个台站处提取了32个地震的走时信息，并利用这些走时信息成功地获得了美国加里福尼亚地区的地壳结构(akiandlee,1976)。他们开创性的工作开启了利用地震层析成像方法来研究地球内部结构的新时代。根据研究范围的大小，地震层析成像可以分为局部区域层析成像、地区层析成像和全球层析成像；根据所用地震波的类型，又可分为体波层析成像和面波层析成像。面波层析成像利用的数据是具有长周期的面波，因而分辨率较低，适用于全球或者地区的成像研究；而体波层析成像利用的数据是具有更短波长的体波，分辨率相对较高，适用于局部区域、地区以及全球等多尺度结构的研究。

在实际工程建设中隐伏的不良地质问题，如岩溶、陷落柱、裂隙、裂缝、断裂破碎带、软弱夹层、地下空洞和不明埋设物等，会造成工程质量隐患或施工困难、浪费，故必须进行精细的工程勘察。由于工程检测要求探测规模很小的异常体，并要求较高的分辨率，因此需要高精度的勘察手段。而在油气田资源勘探过程中，为了获取精确的速度模型，往往开展大量的钻井及测井工作。钻孔勘探的精度较高，但成本高、耗工期，不可能在工区内大量采用，层析成像方法可以弥补这一勘察手段的缺陷。

目前，地震ct从观测方式上来分析可以分为井地ct、井间ct、垂直地震剖面ct，其中井地ct中，由于大多采用井中激发地面接收，射线方向主要是垂直的，故垂向分辨率不如横向分辨率，存在垂向模糊性。而相反，井间ct地震记录上的初至波主要包括透射波、饶射波，由于波传播路径段，具有较高的能量及较宽的频带。波的传播方向基本上是水平的，横向分辨率不如垂向分辨率，存在横向模糊性。

因此，有必要开发一种地震波ct成像技术，能够兼具井间及井地ct的优点，以提高勘探分辨率，实现空间三维探测。

本发明背景技术部分的信息仅仅旨在加深对本发明的一般背景技术的理解，而不应当被视为承认或以任何形式暗示该信息构成已为本领域技术人员所公知的现有技术。

技术实现要素：

本发明提出了一种井间及井地联合地震波ct成像方法及装置，其能够通过野外单炮记录初至波拾取获得地震波旅行时间，并依据测井等资料建立初始速度模型，求解正演射线路径和理论旅行时，以理论和实际旅行时残差为目标函数建立反演方程，逐次反演迭代计算，最终获得高精度的介质速度场，解决了速度异常体勘探以及为地震资料处理提供精确的速度模型。

根据本发明的一方面，提出了一种井间及井地联合地震波ct成像方法，包括以下步骤：

步骤1：建立井间及井地联合观测系统，所述井间及井地联合观测系统包括设置于井间的井间检波器和炮点，以及设置于地面的地面检波器；

步骤2：通过所述井间检波器和地面检波器拾取实际初至波旅行时；

步骤3：根据测井资料建立初始速度模型，根据所述初始速度模型计算理论旅行时和射线路径；

步骤4：以所述实际初至波旅行时与所述理论旅行时的残差作为目标函数，基于sirt算法迭代反演速度修正值，对所述初始速度模型进行修正。

优选地，通过滑动时窗能量比法拾取所述实际初至波旅行时。

优选地，基于sirt算法迭代反演速度修正值，对所述初始速度模型进行修正包括：

步骤401：以sj表示初始模型中的第j个网格的慢度，第k个炮点对应nk个接收点，其中，j＝1,2,…,m，m为总的网格数量，k＝1,2,…,num，num为总的炮点数量；

步骤402：根据射线追踪方法计算第k个炮点对应的接收点的理论旅行时：

其中，表示在第i次迭代中第k个炮点对应的接收点的理论旅行时，i表示迭代次数，i的初始取值为1，anj表示第n条射线在第j个网格内的路径长度；

步骤403：计算第k个炮点对应的接收点的实际初至波旅行时tn和理论旅行时之间的残差：

其中，表示在第i次迭代中第k个炮点对应的接收点的实际初至波旅行时tn和理论旅行时之间的残差；

步骤404：计算每个网格的慢度修正值：

其中，表示针对第k个炮点在第i次迭代中第j个网格的慢度修正值，mj表示射线经过第j个网格的次数，即射线密度；

步骤405：循环执行步骤402至404，针对所有炮点求出第j个网格的慢度修正值，并计算所有炮点的第j个网格的慢度修正值的和

步骤406：对初始速度模型的每个网格的慢度值进行修正：

其中，表示在第i次迭代中第j个网格的慢度值；

步骤407：重复步骤402至406，对i进行迭代，直到残差小于预定值。

根据本发明的另一方面，提出了一种井间及井地联合地震波ct成像装置，包括：

井间及井地联合观测系统，所述井间及井地联合观测系统包括设置于井间的井间检波器和炮点，以及设置于地面的地面检波器；

实际初至波旅行时拾取单元，用于基于所述井间检波器和地面检波器拾取实际初至波旅行时；

初始速度模型建立单元，用于根据测井资料建立初始速度模型，根据所述初始速度模型计算理论旅行时和射线路径；

反演单元，用于以所述实际初至波旅行时与所述理论旅行时的残差作为目标函数，基于sirt算法迭代反演速度修正值，对所述初始速度模型进行修正。

优选地，通过滑动时窗能量比法拾取所述实际初至波旅行时单元。

优选地，基于sirt算法迭代反演速度修正值，对所述初始速度模型进行修正包括以下步骤：

步骤401：以sj表示初始模型中的第j个网格的慢度，第k个炮点对应nk个接收点，其中，j＝1,2,…,m，m为总的网格数量，k＝1,2,…,num，num为总的炮点数量；

步骤402：根据射线追踪方法计算第k个炮点对应的接收点的理论旅行时：

其中，表示在第i次迭代中第k个炮点对应的接收点的理论旅行时，i表示迭代次数，i的初始取值为1，anj表示第n条射线在第j个网格内的路径长度；

步骤403：计算第k个炮点对应的接收点的实际初至波旅行时tn和理论旅行时之间的残差：

其中，表示在第i次迭代中第k个炮点对应的接收点的实际初至波旅行时tn和理论旅行时之间的残差；

步骤404：计算每个网格的慢度修正值：

其中，表示针对第k个炮点在第i次迭代中第j个网格的慢度修正值，mj表示射线经过第j个网格的次数，即射线密度；

步骤405：循环执行步骤402至404，针对所有炮点求出第j个网格的慢度修正值，并计算所有炮点的第j个网格的慢度修正值的和

步骤406：对初始速度模型的每个网格的慢度值进行修正：

其中，表示在第i次迭代中第j个网格的慢度值；

步骤407：重复步骤402至406，对i进行迭代，直到残差小于预定值。

根据本发明的井间及井地联合地震波ct成像方法及装置，其优点在于：

同时设置井间检波器和地面检波器，能够提高数据的观测角度，增加反演区域内射线覆盖密度和穿过异常体垂向及横向的射线，提高反演精度和改善层析成像的分辨率，实现空间三维探测；

井间及井地联合地震波ct成像方法，激发点都位于井中，地震波能量传播距离短，与探测目标体接近，同时降低了低速带对波能量的吸收，地震波具有较高的频率，初至波较清晰；

利用滑动时窗能量比法进行人机交互拾取野外实际初至波旅行时，能够提高拾取精度以及工作效率。

本发明的方法和装置具有其它的特性和优点，这些特性和优点从并入本文中的附图和随后的具体实施例中将是显而易见的，或者将在并入本文中的附图和随后的具体实施例中进行详细陈述，这些附图和具体实施例共同用于解释本发明的特定原理。

附图说明

通过结合附图对本发明示例性实施例进行更详细的描述，本发明的上述以及其它目的、特征和优势将变得更加明显，其中，在本发明示例性实施例中，相同的参考标号通常代表相同部件。

图1示出了根据本发明的井间及井地联合地震波ct成像方法的流程图。

图2示出了根据本发明的示例性实施例井间及井地联合地震波ct成像方法施工布置示意图。

图3a、图3b和图3c示出了根据本发明的示例性实施例的反演得到的目标地区井间的三个速度剖面图。

图4示出了根据本发明的示例性实施例的反演得到的目标地区井地的速度剖面图。

具体实施方式

下面将参照附图更详细地描述本发明。虽然附图中显示了本发明的优选实施例，然而应该理解，可以以各种形式实现本发明而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了使本发明更加透彻和完整，并且能够将本发明的范围完整地传达给本领域的技术人员。

根据示例性实施例的井间及井地联合地震波ct成像方法，包括以下步骤，如图1所示：

步骤1：建立井间及井地联合观测系统，井间及井地联合观测系统包括设置于井间的井间检波器和炮点，以及设置于地面的地面检波器，从而提高射线的空间覆盖密度，增加横向和纵向分辨率。

井间及井地联合地震波ct成像方法采用的观测系统，能够提高数据的观测角度，增加反演区域内射线覆盖密度和穿过异常体垂向及横向的射线，提高反演精度和改善层析成像的分辨率，实现空间三维探测。

步骤2：通过井间检波器和地面检波器拾取实际初至波旅行时。

井间及井地联合地震波ct成像方法，激发点都位于井中，地震波能量传播距离短，与探测目标体接近，同时降低了低速带对波能量的吸收，地震波具有较高的频率，初至波较清晰。

作为优选方案，利用滑动时窗能量比法拾取实际初至波旅行时，能够提高拾取精度以及工作效率。

步骤3：根据测井资料建立初始速度模型，根据初始速度模型计算理论旅行时和射线路径。可以根据实际需要选择现有方法，根据测井资料建立初始速度模型，在此不再赘述。

步骤4：以实际初至波旅行时与理论旅行时的残差作为目标函数，基于sirt算法迭代反演速度修正值，对所述初始速度模型进行修正。

作为优选方案，基于sirt算法迭代反演速度修正值，对所述初始速度模型进行修正，包括以下步骤：

步骤401：以sj表示初始模型中的第j个网格的慢度，第k个炮点对应nk个接收点，其中，j＝1,2,…,m，m为总的网格数量，k＝1,2,…,num，num为总的炮点数量；

步骤402：根据射线追踪方法计算第k个炮点对应的接收点的理论旅行时：

其中，表示在第i次迭代中第k个炮点对应的接收点的理论旅行时，i表示迭代次数，i的初始取值为1，anj表示第n条射线在第j个网格内的路径长度。

步骤403：计算第k个炮点对应的接收点的实际初至波旅行时tn和理论旅行时之间的残差：

其中，表示在第i次迭代中第k个炮点对应的接收点的实际初至波旅行时tn和理论旅行时之间的残差；

步骤404：计算每个网格的慢度修正值

在示例性实施例中，在一次迭代中使用所有相关数据同时对慢度sj进行修正。具体来说，通过以下公式计算第n条射线在第j个网格内的慢度修正值：

其中，表示在第i次迭代中第n条射线在第j个网格内的慢度修正值；

然后，将所有射线在第j个网格内的慢度修正值进行平均，作为第i次迭代中第j个网格的慢度修正值，如以下公式所示：

其中，表示针对第k个炮点在第i次迭代中第j个网格的慢度修正值，mj表示射线经过第j个网格的次数，即射线密度；

步骤405：循环执行步骤402至404，针对所有炮点求出第j个网格的慢度修正值，并计算所有炮点的第j个网格的慢度修正值的和

步骤406：根据以下公式对初始速度模型的每个网格的慢度值进行修正：

其中，表示在第i次迭代中第j个网格的慢度值；

步骤407：重复步骤402至406，对i进行迭代，直到残差小于预定值。

上述步骤为一次完整的迭代求取过程，sirt迭代算法不会出现解的奇异现象，同时能够平稳收敛，并且计算简便快捷。完成上述步骤后，可以输出反演结果。

实施例1

在该实施例中，将根据本发明的井间及井地联合地震波ct成像方法应用于山西某煤层气区域，以求取该地区精确的速度模型。图1示出了根据本发明的井间及井地联合地震波ct成像方法的流程图，其包括以下步骤：

步骤1：建立井间及井地联合观测系统，井间及井地联合观测系统包括设置于井间的井间检波器和炮点，以及设置于地面的地面检波器，图2示出了根据本发明的示例性实施例井间及井地联合地震波ct成像方法施工布置示意图；

步骤2：通过井间检波器和地面检波器拾取实际初至波旅行时。

通过滑动时窗能量比法拾取实际初至波旅行时，提高拾取精度以及工作效率。

步骤3：根据测井资料建立初始速度模型，根据初始速度模型计算理论旅行时和射线路径；

步骤4：以实际初至波旅行时与理论旅行时的残差作为目标函数，基于sirt算法迭代反演速度修正值，对初始速度模型进行修正。

图3a、图3b和图3c示出了根据本发明示例性实施例反演得到的目标地区井间的三个速度剖面图，图中色标表示速度值。分析图中的三个速度剖面图可以看出，纵向速度分布由浅入深速度递增符合地层分布规律，与相关地质资料吻合。其中部分区域出现高速异常(图中圈出)，是由于小窑开采形成的采空区。

图4示出了根据本发明示例性实施的反演得到的目标地区井地的速度剖面图。图中色标表示速度值，黑色横线表示煤层，反演速度模型符合地层分布规律，同样出现采空区导致的高速异常区。

单纯的井地ct成像地震射线能够实现浅部大范围速度重建，且快速经济。但在深部其地震波难以横向传播，因此形成倒三角速度模型，无法对深部地层做出很好的解释。因此结合井间ct成像，形成井间及井地联合地震波ct成像方法能够很好的弥补横向信息不足。

本领域技术人员应理解，上面对本发明的实施例的描述的目的仅为了示例性地说明本发明的实施例的有益效果，并不意在将本发明的实施例限制于所给出的任何示例。

以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。本文中所用术语的选择，旨在最好地解释各实施例的原理、实际应用或对市场中的技术的改进，或者使本技术领域的其它普通技术人员能理解本文披露的各实施例。