CN107942372A - 井间及井地联合地震波ct成像方法及装置 - Google Patents

井间及井地联合地震波ct成像方法及装置 Download PDF

Info

Publication number
CN107942372A
CN107942372A CN201610889087.1A CN201610889087A CN107942372A CN 107942372 A CN107942372 A CN 107942372A CN 201610889087 A CN201610889087 A CN 201610889087A CN 107942372 A CN107942372 A CN 107942372A
Authority
CN
China
Prior art keywords
mrow
msubsup
msub
well
munderover
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN201610889087.1A
Other languages
English (en)
Inventor
须振华
周单
林正良
滕龙
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Petroleum and Chemical Corp
Sinopec Geophysical Research Institute
Original Assignee
China Petroleum and Chemical Corp
Sinopec Geophysical Research Institute
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China Petroleum and Chemical Corp, Sinopec Geophysical Research Institute filed Critical China Petroleum and Chemical Corp
Priority to CN201610889087.1A priority Critical patent/CN107942372A/zh
Publication of CN107942372A publication Critical patent/CN107942372A/zh
Pending legal-status Critical Current

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/282—Application of seismic models, synthetic seismograms
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/303—Analysis for determining velocity profiles or travel times
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/303—Analysis for determining velocity profiles or travel times
    • G01V1/305—Travel times
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/62—Physical property of subsurface
    • G01V2210/622—Velocity, density or impedance
    • G01V2210/6222—Velocity; travel time

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Analysing Materials By The Use Of Radiation (AREA)

Abstract

本发明公开了一种井间及井地联合地震波CT成像方法及装置。该方法包括以下步骤：步骤1：建立井间及井地联合观测系统，其包括设置于井间的井间检波器和炮点，以及设置于地面的地面检波器；步骤2：通过井间检波器和地面检波器拾取实际初至波旅行时；步骤3：根据测井资料建立初始速度模型，根据初始速度模型计算理论旅行时和射线路径；步骤4：以实际初至波旅行时与理论旅行时的残差作为目标函数，基于SIRT算法迭代反演速度修正值，对初始速度模型进行修正。该方法的优点在于能够兼具井间及井地CT成像的优点，以提高勘探分辨率，实现空间三维探测。

Description

井间及井地联合地震波CT成像方法及装置
技术领域
本发明涉及地震勘探及开发领域，更具体地，涉及一种井间及井地联合地震波CT成像方法。
背景技术
地震层析成像是一种基于模型的数据拟合方法。与医学中的计算机辅助层析成像(Computer Tomography，CT)类似并受其启发，地震层析成像通过分析台站记录到的由天然地震或者人工地震激发的地震波信号来获取地球内部结构的信息。
1976年，Aki和Lee从60个台站处提取了32个地震的走时信息，并利用这些走时信息成功地获得了美国加里福尼亚地区的地壳结构(Aki and Lee,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成像方法能够很好的弥补横向信息不足。
本领域技术人员应理解，上面对本发明的实施例的描述的目的仅为了示例性地说明本发明的实施例的有益效果，并不意在将本发明的实施例限制于所给出的任何示例。
以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。本文中所用术语的选择，旨在最好地解释各实施例的原理、实际应用或对市场中的技术的改进，或者使本技术领域的其它普通技术人员能理解本文披露的各实施例。

Claims (6)

1.一种井间及井地联合地震波CT成像方法，包括以下步骤：
步骤1：建立井间及井地联合观测系统，所述井间及井地联合观测系统包括设置于井间的井间检波器和炮点，以及设置于地面的地面检波器；
步骤2：通过所述井间检波器和地面检波器拾取实际初至波旅行时；
步骤3：根据测井资料建立初始速度模型，根据所述初始速度模型计算理论旅行时和射线路径；
步骤4：以所述实际初至波旅行时与所述理论旅行时的残差作为目标函数，基于SIRT算法迭代反演速度修正值，对所述初始速度模型进行修正。
2.根据权利要求1所述的井间及井地联合地震波CT成像方法，其中，通过滑动时窗能量比法拾取所述实际初至波旅行时。
3.根据权利要求1所述的井间及井地联合地震波CT成像方法，其中，基于SIRT算法迭代反演速度修正值，对所述初始速度模型进行修正包括：
步骤401：以Sj表示初始模型中的第j个网格的慢度，第k个炮点对应Nk个接收点，其中，j＝1,2,…,M，M为总的网格数量，k＝1,2,…,NUM，NUM为总的炮点数量；
步骤402：根据射线追踪方法计算第k个炮点对应的接收点的理论旅行时：
<mrow> <msubsup> <mi>P</mi> <mi>n</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </msubsup> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </munderover> <msub> <mi>a</mi> <mrow> <mi>n</mi> <mi>j</mi> </mrow> </msub> <msubsup> <mi>S</mi> <mi>j</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <mi>n</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>,</mo> <mn>3</mn> <mo>...</mo> <mo>...</mo> <msub> <mi>N</mi> <mi>k</mi> </msub> <mo>;</mo> </mrow>
其中，表示在第i次迭代中第k个炮点对应的接收点的理论旅行时，i表示迭代次数，i的初始取值为1，anj表示第n条射线在第j个网格内的路径长度；
步骤403：计算第k个炮点对应的接收点的实际初至波旅行时Tn和理论旅行时之间的残差：
<mrow> <msubsup> <mi>&amp;Delta;t</mi> <mi>n</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </msubsup> <mo>=</mo> <msub> <mi>T</mi> <mi>n</mi> </msub> <mo>-</mo> <msubsup> <mi>P</mi> <mi>n</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <mi>n</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>,</mo> <mn>3</mn> <mo>...</mo> <mo>...</mo> <msub> <mi>N</mi> <mi>k</mi> </msub> <mo>;</mo> </mrow>
其中，表示在第i次迭代中第k个炮点对应的接收点的实际初至波旅行时Tn和理论旅行时之间的残差；
步骤404：计算每个网格的慢度修正值：
<mrow> <msubsup> <mi>&amp;Delta;S</mi> <mi>j</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </msubsup> <mo>=</mo> <mfrac> <mn>1</mn> <msub> <mi>M</mi> <mi>j</mi> </msub> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>n</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mi>k</mi> </msub> </munderover> <mrow> <mo>(</mo> <mfrac> <mrow> <msub> <mi>a</mi> <mrow> <mi>n</mi> <mi>j</mi> </mrow> </msub> <msubsup> <mi>&amp;Delta;t</mi> <mi>n</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </msubsup> </mrow> <mrow> <msubsup> <mi>&amp;Sigma;</mi> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </msubsup> <msubsup> <mi>a</mi> <mrow> <mi>n</mi> <mi>j</mi> </mrow> <mn>2</mn> </msubsup> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>,</mo> <mrow> <mo>(</mo> <mi>j</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>...</mo> <mi>M</mi> <mo>)</mo> </mrow> </mrow>
其中，表示针对第k个炮点在第i次迭代中第j个网格的慢度修正值，Mj表示射线经过第j个网格的次数，即射线密度；
步骤405：循环执行步骤402至404，针对所有炮点求出第j个网格的慢度修正值，并计算所有炮点的第j个网格的慢度修正值的和
步骤406：对初始速度模型的每个网格的慢度值进行修正：
<mrow> <msubsup> <mi>S</mi> <mi>j</mi> <mi>i</mi> </msubsup> <mo>=</mo> <msubsup> <mi>S</mi> <mi>j</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> <mo>+</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>1</mn> </mrow> <mrow> <mi>N</mi> <mi>U</mi> <mi>M</mi> </mrow> </munderover> <msubsup> <mi>&amp;Delta;S</mi> <mi>j</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </msubsup> <mo>=</mo> <msubsup> <mi>S</mi> <mi>j</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> <mo>+</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>1</mn> </mrow> <mrow> <mi>N</mi> <mi>U</mi> <mi>M</mi> </mrow> </munderover> <mrow> <mo>(</mo> <mfrac> <mn>1</mn> <msub> <mi>M</mi> <mi>j</mi> </msub> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>n</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mi>k</mi> </msub> </munderover> <mo>(</mo> <mfrac> <mrow> <msub> <mi>a</mi> <mrow> <mi>n</mi> <mi>j</mi> </mrow> </msub> <msubsup> <mi>&amp;Delta;t</mi> <mi>n</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </msubsup> </mrow> <mrow> <msubsup> <mi>&amp;Sigma;</mi> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </msubsup> <msubsup> <mi>a</mi> <mrow> <mi>n</mi> <mi>j</mi> </mrow> <mn>2</mn> </msubsup> </mrow> </mfrac> <mo>)</mo> <mo>)</mo> </mrow> <mo>,</mo> <mrow> <mo>(</mo> <mi>j</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>...</mo> <mi>M</mi> <mo>)</mo> </mrow> </mrow>
其中，表示在第i次迭代中第j个网格的慢度值；
步骤407：重复步骤402至406，对i进行迭代，直到残差小于预定值。
4.一种井间及井地联合地震波CT成像装置，包括：
井间及井地联合观测系统，所述井间及井地联合观测系统包括设置于井间的井间检波器和炮点，以及设置于地面的地面检波器；
实际初至波旅行时拾取单元，用于基于所述井间检波器和地面检波器拾取实际初至波旅行时；
初始速度模型建立单元，用于根据测井资料建立初始速度模型，根据所述初始速度模型计算理论旅行时和射线路径；
反演单元，用于以所述实际初至波旅行时与所述理论旅行时的残差作为目标函数，基于SIRT算法迭代反演速度修正值，对所述初始速度模型进行修正。
5.根据权利要求4所述的井间及井地联合地震波CT成像装置，其中，通过滑动时窗能量比法拾取所述实际初至波旅行时单元。
6.根据权利要求4所述的井间及井地联合地震波CT成像装置，其中，基于SIRT算法迭代反演速度修正值，对所述初始速度模型进行修正包括以下步骤：
步骤401：以Sj表示初始模型中的第j个网格的慢度，第k个炮点对应Nk个接收点，其中，j＝1,2,…,M，M为总的网格数量，k＝1,2,…,NUM，NUM为总的炮点数量；
步骤402：根据射线追踪方法计算第k个炮点对应的接收点的理论旅行时：
<mrow> <msubsup> <mi>P</mi> <mi>n</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </msubsup> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </munderover> <msub> <mi>a</mi> <mrow> <mi>n</mi> <mi>j</mi> </mrow> </msub> <msubsup> <mi>S</mi> <mi>j</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <mi>n</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>,</mo> <mn>3</mn> <mo>...</mo> <mo>...</mo> <msub> <mi>N</mi> <mi>k</mi> </msub> <mo>;</mo> </mrow>
其中，表示在第i次迭代中第k个炮点对应的接收点的理论旅行时，i表示迭代次数，i的初始取值为1，anj表示第n条射线在第j个网格内的路径长度；
步骤403：计算第k个炮点对应的接收点的实际初至波旅行时Tn和理论旅行时之间的残差：
<mrow> <msubsup> <mi>&amp;Delta;t</mi> <mi>n</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </msubsup> <mo>=</mo> <msub> <mi>T</mi> <mi>n</mi> </msub> <mo>-</mo> <msubsup> <mi>P</mi> <mi>n</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <mi>n</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>,</mo> <mn>3</mn> <mo>...</mo> <mo>...</mo> <msub> <mi>N</mi> <mi>k</mi> </msub> <mo>;</mo> </mrow>
其中，表示在第i次迭代中第k个炮点对应的接收点的实际初至波旅行时Tn和理论旅行时之间的残差；
步骤404：计算每个网格的慢度修正值：
<mrow> <msubsup> <mi>&amp;Delta;S</mi> <mi>j</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </msubsup> <mo>=</mo> <mfrac> <mn>1</mn> <msub> <mi>M</mi> <mi>j</mi> </msub> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>n</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mi>k</mi> </msub> </munderover> <mrow> <mo>(</mo> <mfrac> <mrow> <msub> <mi>a</mi> <mrow> <mi>n</mi> <mi>j</mi> </mrow> </msub> <msubsup> <mi>&amp;Delta;t</mi> <mi>n</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </msubsup> </mrow> <mrow> <msubsup> <mi>&amp;Sigma;</mi> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </msubsup> <msubsup> <mi>a</mi> <mrow> <mi>n</mi> <mi>j</mi> </mrow> <mn>2</mn> </msubsup> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>,</mo> <mrow> <mo>(</mo> <mi>j</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>...</mo> <mi>M</mi> <mo>)</mo> </mrow> </mrow>
其中，表示针对第k个炮点在第i次迭代中第j个网格的慢度修正值，Mj表示射线经过第j个网格的次数，即射线密度；
步骤405：循环执行步骤402至404，针对所有炮点求出第j个网格的慢度修正值，并计算所有炮点的第j个网格的慢度修正值的和
步骤406：对初始速度模型的每个网格的慢度值进行修正：
<mrow> <msubsup> <mi>S</mi> <mi>j</mi> <mi>i</mi> </msubsup> <mo>=</mo> <msubsup> <mi>S</mi> <mi>j</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> <mo>+</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>1</mn> </mrow> <mrow> <mi>N</mi> <mi>U</mi> <mi>M</mi> </mrow> </munderover> <msubsup> <mi>&amp;Delta;S</mi> <mi>j</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </msubsup> <mo>=</mo> <msubsup> <mi>S</mi> <mi>j</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> <mo>+</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>1</mn> </mrow> <mrow> <mi>N</mi> <mi>U</mi> <mi>M</mi> </mrow> </munderover> <mrow> <mo>(</mo> <mfrac> <mn>1</mn> <msub> <mi>M</mi> <mi>j</mi> </msub> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>n</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mi>k</mi> </msub> </munderover> <mo>(</mo> <mfrac> <mrow> <msub> <mi>a</mi> <mrow> <mi>n</mi> <mi>j</mi> </mrow> </msub> <msubsup> <mi>&amp;Delta;t</mi> <mi>n</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </msubsup> </mrow> <mrow> <msubsup> <mi>&amp;Sigma;</mi> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </msubsup> <msubsup> <mi>a</mi> <mrow> <mi>n</mi> <mi>j</mi> </mrow> <mn>2</mn> </msubsup> </mrow> </mfrac> <mo>)</mo> <mo>)</mo> </mrow> <mo>,</mo> <mrow> <mo>(</mo> <mi>j</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>...</mo> <mi>M</mi> <mo>)</mo> </mrow> </mrow>
其中，表示在第i次迭代中第j个网格的慢度值；
步骤407：重复步骤402至406，对i进行迭代，直到残差小于预定值。
CN201610889087.1A 2016-10-12 2016-10-12 井间及井地联合地震波ct成像方法及装置 Pending CN107942372A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610889087.1A CN107942372A (zh) 2016-10-12 2016-10-12 井间及井地联合地震波ct成像方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610889087.1A CN107942372A (zh) 2016-10-12 2016-10-12 井间及井地联合地震波ct成像方法及装置

Publications (1)

Publication Number Publication Date
CN107942372A true CN107942372A (zh) 2018-04-20

Family

ID=61928216

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610889087.1A Pending CN107942372A (zh) 2016-10-12 2016-10-12 井间及井地联合地震波ct成像方法及装置

Country Status (1)

Country Link
CN (1) CN107942372A (zh)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109188528A (zh) * 2018-08-10 2019-01-11 武汉市工程科学技术研究院 井间弹性波层析成像系统和方法
CN110820814A (zh) * 2019-11-14 2020-02-21 山东大学 一种桩基检测装置及方法
CN112305595A (zh) * 2019-07-24 2021-02-02 中国石油化工股份有限公司 基于折射波分析地质体结构的方法及存储介质
CN112363211A (zh) * 2020-11-23 2021-02-12 同济大学 一种改进的sirt法射线走时层析成像方法
CN112485825A (zh) * 2019-09-11 2021-03-12 中国石油化工股份有限公司 一种基于初至波走时层析的微测井解释方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP0297737A2 (en) * 1987-06-29 1989-01-04 Conoco Inc. Three-dimensional iterative structural modeling of seismic data
CN102877828A (zh) * 2012-09-09 2013-01-16 山西山地物探技术有限公司 一种三维多井联合井地ct成像方法
CN103777238A (zh) * 2012-10-17 2014-05-07 中国石油化工股份有限公司 一种纯纵波各向异性波场模拟方法
CN105301639A (zh) * 2015-10-21 2016-02-03 中国石油天然气集团公司 基于vsp旅行时双加权层析反演速度场的方法及其装置

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP0297737A2 (en) * 1987-06-29 1989-01-04 Conoco Inc. Three-dimensional iterative structural modeling of seismic data
CN102877828A (zh) * 2012-09-09 2013-01-16 山西山地物探技术有限公司 一种三维多井联合井地ct成像方法
CN103777238A (zh) * 2012-10-17 2014-05-07 中国石油化工股份有限公司 一种纯纵波各向异性波场模拟方法
CN105301639A (zh) * 2015-10-21 2016-02-03 中国石油天然气集团公司 基于vsp旅行时双加权层析反演速度场的方法及其装置

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
张连伟: "井间地震CT技术及其在铁路岩溶勘察中的应用", 《中国优秀硕士学位论文全文数据库 工程科技Ⅱ辑》 *
曾来: "井中地震CT观测系统研究", 《上海国土资源》 *

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109188528A (zh) * 2018-08-10 2019-01-11 武汉市工程科学技术研究院 井间弹性波层析成像系统和方法
CN109188528B (zh) * 2018-08-10 2020-04-17 武汉市工程科学技术研究院 井间弹性波层析成像系统和方法
CN112305595A (zh) * 2019-07-24 2021-02-02 中国石油化工股份有限公司 基于折射波分析地质体结构的方法及存储介质
CN112305595B (zh) * 2019-07-24 2024-05-17 中国石油化工股份有限公司 基于折射波分析地质体结构的方法及存储介质
CN112485825A (zh) * 2019-09-11 2021-03-12 中国石油化工股份有限公司 一种基于初至波走时层析的微测井解释方法
CN112485825B (zh) * 2019-09-11 2024-04-09 中国石油化工股份有限公司 一种基于初至波走时层析的微测井解释方法
CN110820814A (zh) * 2019-11-14 2020-02-21 山东大学 一种桩基检测装置及方法
CN112363211A (zh) * 2020-11-23 2021-02-12 同济大学 一种改进的sirt法射线走时层析成像方法

Similar Documents

Publication Publication Date Title
CN107942372A (zh) 井间及井地联合地震波ct成像方法及装置
CN102426384B (zh) 一种探测地下采空区和岩溶分布的方法
CN112485823B (zh) 高效综合超前地质预报方法
Cueto et al. Karst‐induced sinkhole detection using an integrated geophysical survey: a case study along the Riyadh Metro Line 3 (Saudi Arabia)
CN112965136B (zh) 一种富水岩溶隧道的多手段超前探测方法
Esmailzadeh et al. Prediction of rock mass rating using TSP method and statistical analysis in Semnan Rooziyeh spring conveyance tunnel
Holma et al. Future prospects of muography for geological research and geotechnical and mining engineering
Takahashi et al. ISRM suggested methods for borehole geophysics in rock engineering
CN102877828A (zh) 一种三维多井联合井地ct成像方法
CN105607119B (zh) 近地表模型构建方法与静校正量求取方法
CN103592680A (zh) 一种基于正反演的测井数据和深度域地震剖面合成方法
Zhang et al. Excavation-induced structural deterioration of rock masses at different depths
CN113376695B (zh) 一种适用于煤层底板复杂陷落柱的全波形反演方法
Toney et al. Joint body‐and surface‐wave tomography of Yucca Flat, Nevada, using a novel seismic source
CN114280669A (zh) 一种基于折射波周期振幅衰减的薄煤带探测方法及系统
CN117192615A (zh) 基于透射地震尾波的采煤工作面内隐伏地质构造探测方法
Cosma et al. Seismic characterization of fracturing at the Äspö Hard Rock Laboratory, Sweden, from the kilometer scale to the meter scale
Gu et al. Investigation of fractures using seismic computerized crosshole tomography
CN109375251B (zh) 利用城市既有地下空间与地表的探测方法及系统
Gunnarsson 3D modeling in Petrel of geological CO2 storage site
CN110687591B (zh) 基于先验数据的波形匹配确定煤层及围岩物性参数的方法
Fan et al. Application of Seismic Channel Wave Technology on Small Structure Exploration in Coal Mine
Lizhi et al. The seismic CT method in measuring rock bodies
CN105242317A (zh) 一种纵波速度的确定方法及装置
CN112305594A (zh) 非均质储层的油气分布确定方法及系统

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
RJ01 Rejection of invention patent application after publication
RJ01 Rejection of invention patent application after publication

Application publication date: 20180420