CN110031898A - 数据优化方法及积分法叠前深度偏移方法 - Google Patents

数据优化方法及积分法叠前深度偏移方法 Download PDF

Info

Publication number
CN110031898A
CN110031898A CN201910397720.9A CN201910397720A CN110031898A CN 110031898 A CN110031898 A CN 110031898A CN 201910397720 A CN201910397720 A CN 201910397720A CN 110031898 A CN110031898 A CN 110031898A
Authority
CN
China
Prior art keywords
imaging
equivalent
value
point
optimization
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.)
Granted
Application number
CN201910397720.9A
Other languages
English (en)
Other versions
CN110031898B (zh
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.)
Institute of Geology and Geophysics of CAS
Original Assignee
Institute of Geology and Geophysics of CAS
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 Institute of Geology and Geophysics of CAS filed Critical Institute of Geology and Geophysics of CAS
Priority to CN201910397720.9A priority Critical patent/CN110031898B/zh
Publication of CN110031898A publication Critical patent/CN110031898A/zh
Priority to US16/556,430 priority patent/US11209563B2/en
Application granted granted Critical
Publication of CN110031898B publication Critical patent/CN110031898B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/362—Effecting static or dynamic corrections; Stacking
    • 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/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T1/00—General purpose image data processing
    • G06T1/20—Processor architectures; Processor configuration, e.g. pipelining
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/50—Corrections or adjustments related to wave propagation
    • G01V2210/51—Migration
    • G01V2210/512—Pre-stack

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Remote Sensing (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geophysics (AREA)
  • Acoustics & Sound (AREA)
  • Geology (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Computing Systems (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Image Processing (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明提供了一种数据优化方法及积分法叠前深度偏移方法，涉及地震勘探技术领域，该优化方法包括获取待优化的目标矩阵；根据目标矩阵生成第一序列；根据预设的网格密度抽稀第一序列得到第二序列各个元素的取值位置，并基于最小二乘原理求取第二序列各个元素的值；对第二序列进行插值得到第三序列；计算第三序列对应的目标矩阵；计算待优化的目标矩阵与第三序列对应的目标矩阵之间的误差；当该误差小于第一误差阈值时，记录上述第二序列对应的目标矩阵为待优化的目标矩阵的优化目标矩阵。本发明实施例提供的一种数据优化方法，可以使得利用GPU进行补偿介质吸收的Kirchhoff积分法叠前深度偏移技术并行计算的效率大幅提高。

Description

数据优化方法及积分法叠前深度偏移方法
技术领域
本发明涉及地震勘探技术领域，尤其是涉及一种数据优化方法及积分法叠前深度偏移方法。
背景技术
现行地震叠前偏移方法的实践中，以Kirchhoff方法为代表的积分法叠前偏移技术占据着主体地位。Kirchhoff方法地震叠前偏移主要分为两类，Kirchhoff叠前时间偏移和Kirchhoff叠前深度偏移。
为实现大地对地震波的吸收补偿，我们需要在Kirchhoff叠前深度偏移实现流程中引入与地震波传播路径有关的等效Q值，这一等效Q值与地震波自炮点至成像点，再至检波点的传播路径上的所有层Q值有关。因此，这一补偿介质吸收Kirchhoff叠前深度偏移技术中，要得到正确的偏移幅值，除需事先计算并存储走时表外，还需额外存储一个等效Q值表。然而，对于具有高覆盖次数的高空间、时间采样密度地震数据，其走时表与等效Q值表所需存储量巨大，并且与其覆盖次数是呈正相关的。与之形成鲜明对比的是，GPU卡的主要特点是擅长计算，其存储量非常有限。
目前，采用CPU读取走时表和等效Q值表，再通过CPU-GPU之间的通讯，以解决利用GPU进行偏移幅值计算时的走时表、等效Q值表需求的策略，会大幅降低GPU并行计算的计算效率，使得利用GPU进行补偿介质吸收Kirchhoff叠前深度偏移计算时的计算效率低下。
发明内容
有鉴于此，本发明的目的在于提供一种数据优化方法及积分法叠前深度偏移方法，可以实现以较小的“走时表和等效Q值表”逼近原先较大的“走时表和等效Q值表”，在满足计算精度的条件下，大大减少对GPU的存储需求，使得利用GPU进行补偿介质吸收的Kirchhoff积分法叠前深度偏移技术并行计算的效率大幅提高。
第一方面，本发明实施例提供了一种数据优化方法，包括：获取待优化的目标矩阵，该目标矩阵为走时表或等效Q值表；根据该目标矩阵生成第一序列；根据预设的网格密度抽稀第一序列，得到第二序列各个元素的取值位置，并基于最小二乘原理求取第二序列各个元素的值；第二序列的元素个数远少于第一序列的元素个数；对第二序列进行插值得到第三序列；该第三序列和第一序列的元素个数相等；计算第三序列对应的目标矩阵；计算待优化的目标矩阵与第三序列对应的目标矩阵之间的误差；比较该误差与预设第一误差阈值的大小，当该误差大于预设第一误差阈值时，按预设规则增大网格密度，并继续执行根据预设的网格密度抽稀第一序列得到第二序列各个元素的取值位置的步骤；当该误差小于第一误差阈值时，记录上述第二序列对应的目标矩阵为该待优化的目标矩阵的优化目标矩阵。
结合第一方面，本发明实施例提供了第一方面的第一种可能的实施方式，其中，上述根据该目标矩阵生成第一序列的步骤，包括：当该目标矩阵为走时表时，先将走时表转化为慢度表，再根据该慢度表的所有元素生成第一序列；当该目标矩阵为等效Q值表时，根据等效Q值表的所有元素生成第一序列。
结合第一方面，本发明实施例提供了第一方面的第二种可能的实施方式，其中，上述基于最小二乘原理求取第二序列各个元素的值的步骤，包括：对第二序列中各个元素求偏导数，并令各个元素的偏导数为0，得到关于各个元素的方程组；求解该方程组得到各个元素的值。
结合第一方面，本发明实施例提供了第一方面的第三种可能的实施方式，其中，上述计算待优化的目标矩阵与第三序列对应的目标矩阵之间的误差的步骤，包括：计算第三序列对应的目标矩阵与该待优化的目标矩阵中对应各元素差值的绝对值；统计该差值的绝对值大于预设偏差阈值的元素数量占该待优化的目标矩阵的元素总数的百分比；将该百分比作为第三序列对应的目标矩阵与该待优化的目标矩阵的误差。
结合第一方面，本发明实施例提供了第一方面的第四种可能的实施方式，其中，在上述当该误差小于第一误差阈值时，记录第二序列对应的目标矩阵为该待优化的目标矩阵的优化目标矩阵的步骤之后，该方法还包括：记录该网格密度，将该网格密度对应的采样个数及间距作为该目标矩阵的共享优化参数。
第二方面，本发明实施例还提供了一种积分法叠前深度偏移方法，应用于CPU，该方法包括：获取目标工区的叠前地震数据、偏移参数，以及该目标工区成像空间的成像线走时表和成像线等效Q值表；运用第一方面提供的数据优化方法优化该成像线走时表和该成像线等效Q值表得到成像线优化走时表和成像线优化等效Q值表；以成像线为单位存储该成像线优化走时表和该成像线优化等效Q值表；将叠前地震数据、偏移参数，以及成像线优化走时表和成像线优化等效Q值表以成像线为单位发送给GPU，以根据该成像线优化走时表、该成像线优化等效Q值表、该叠前地震数据及该偏移参数，计算目标工区成像空间内各个成像点对应不同偏移距的偏移幅值。
结合第二方面，本发明实施例提供了第二方面的第一种可能的实施方式，其中，上述以成像线为单位存储该成像线优化走时表和该成像线优化等效Q值表的步骤，包括：以成像线为单位构建成像线优化走时表的第一索引文件和成像线优化等效Q值表的第二索引文件；该第一索引文件包括成像线的成像块划分参数，每个成像块内的共享走时表优化参数，以及每个成像块内第一个成像点的优化走时表在成像线优化走时表的数据文件中的存储位置；该第二索引文件包括成像块划分参数，每个成像块内的共享等效Q值表优化参数，以及每个成像块内第一个成像点的优化等效Q值表在成像线优化等效Q值表的数据文件中的存储位置；存储第一索引文件和第二索引文件，并以无符号短整型变量类型存储成像线优化走时表的数据文件和成像线优化等效Q值表的数据文件。
结合第二方面的第一种可能的实施方式，本发明实施例提供了第二方面的第二种可能的实施方式，其中，上述以无符号短整型变量类型存储成像线优化走时表的数据文件和成像线优化等效Q值表的数据文件的步骤，包括：根据第一关系式、第二关系式分别逐成像块存储该成像线优化走时表的数据文件、该成像线优化等效Q值表的数据文件；第一关系式为：
第二关系式为：
式中，ixi为成像点在成像块内位置x方向索引，izi为成像点在成像块内位置z方向索引，为该成像点优化走时表中炮点在x方向索引，为该成像点优化走时表中炮点在y方向索引；为该成像点优化等效Q值表中炮点在x方向索引，为该成像点优化等效Q值表中炮点在y方向索引；k为maxτ为该目标工区的成像空间优化走时表中走时的最大值，b为0.0；表示该成像块优化走时表数据文件的存储顺序；表示该成像块优化走时表的数据格式；表示该成像块优化等效Q值表的数据格式；表示该成像块优化等效Q值表数据文件的存储顺序。
第三方面，本发明实施例还提供了一种积分法叠前深度偏移方法，应用于GPU，该方法包括：接收CPU发送的目标工区的叠前地震数据、偏移参数，以及该CPU发送的以成像线为单位的，目标工区成像空间的成像线优化走时表和成像线优化等效Q值表；该成像空间被划分为多条成像线，该成像线被划分为多个成像块，该成像块包含多个成像点；该成像线优化走时表和该成像线优化等效Q值表，为上述CPU运用第一方面提供的数据优化方法优化该目标工区的成像线走时表和成像线等效Q值表得到；根据成像线优化走时表、成像线优化等效Q值表、叠前地震数据及偏移参数，计算目标工区成像空间内各个成像点对应不同偏移距的偏移幅值。
结合第三方面，本发明实施例提供了第三方面的第一种可能的实施方式，其中，上述根据该成像线优化走时表、成像线优化等效Q值表、叠前地震数据及偏移参数，计算该目标工区成像空间内各个成像点对应不同偏移距的偏移幅值的步骤，包括：获取叠前地震道的炮点坐标和检波点坐标；确定GPU线程块的大小和GPU线程块网格的大小；根据该GPU线程块的大小和该GPU线程块网格的大小计算该成像线中每个成像点的线程索引，该线程索引包括深度位置索引和横向位置索引；根据该线程索引线程读取成像线的第一索引文件和第二索引文件，以确定每个成像点所在的成像块的走时表优化参数和等效Q值表优化参数，以及该成像点所在成像块内第一个成像点的优化走时表在成像线优化走时表的数据文件中的存储位置，和该成像点所在成像块内第一个成像点的优化等效Q值表在成像线优化等效Q值表的数据文件中的存储位置；根据该炮点坐标、该检波点坐标，以及该成像线对应的第一索引文件、第二索引文件、成像线优化走时表的数据文件和成像线优化等效Q值表的数据文件，线程获取叠前地震道对应的炮点到该成像点的走时、等效Q值，叠前地震道对应的检波点到该成像点的走时、等效Q值；根据成像公式线程计算各个成像点的偏移幅值；该成像公式为：
式中，I(h，xp，yp，zp)为该成像点的偏移幅值，h为该叠前地震道偏移距，j为虚部单位，ω为频率，ω0为地震数据主频，F(ω)为该叠前地震道变换至频率域的序列，rs，rg分别为该叠前地震道对应的炮点、检波点至该成像点的距离，τs、Qs分别为该叠前地震道对应的炮点到该成像点的走时、等效Q值，τg、Qg分别为该叠前地震道对应的检波点到该成像点的走时、等效Q值。
结合第三方面的第一种可能的实施方式，本发明实施例提供了第三方面的第二种可能的实施方式，其中，上述线程获取该叠前地震道对应的炮点到该成像点的走时，该叠前地震道对应的检波点到该成像点的走时的步骤，包括：根据该炮点坐标、检波点坐标、走时表优化参数，计算该炮点坐标对应的优化走时表的第一网格及该第一网格四个顶点的第一索引值，以及该检波点坐标对应的优化走时表的第二网格及该第二网格四个顶点的第二索引值；根据该第一索引值、第二索引值，分别计算该第一网格四个顶点、该第二网格四个顶点的走时在成像线优化走时表的数据文件中的第一位置、第二位置，并根据该第一位置、第二位置读取第一网格四个顶点的走时、第二网格四个顶点的走时；根据该第一网格四个顶点的走时、第二网格四个顶点的走时分别计算第一网格四个顶点的慢度值、第二网格四个顶点的慢度值；根据该第一网格四个顶点的慢度值、第二网格四个顶点的慢度值分别插值得到该成像点的第一慢度值、第二慢度值；根据该第一慢度值、第二慢度值分别计算该叠前地震道对应的炮点到该成像点的走时、该叠前地震道对应的检波点到该成像点的走时。
结合第三方面的第一种可能的实施方式，本发明实施例提供了第三方面的第三种可能的实施方式，其中，上述线程获取该叠前地震道对应的炮点到该成像点的等效Q值，该叠前地震道对应的检波点到该成像点的等效Q值的步骤，包括：根据该炮点坐标、检波点坐标、等效Q值优化参数，计算该炮点坐标对应的优化等效Q值的第三网格及该第三网格四个顶点的第三索引值，以及该检波点坐标对应的优化等效Q值的第四网格及该第四网格四个顶点的第四索引值；根据该第三索引值、第四索引值，分别计算该第三网格四个顶点、第四网格四个顶点的等效Q值在成像线优化等效Q值表的数据文件中的第三位置、第四位置，并根据该第三位置、第四位置读取该第三网格四个顶点的等效Q值、第四网格四个顶点的等效Q值；根据该第三网格四个顶点的等效Q值、第四网格四个顶点的等效Q值，分别插值得到该叠前地震道对应的炮点到该成像点的等效Q值、该叠前地震道对应的检波点到该成像点的等效Q值。
第四方面，本发明实施例还提供了一种积分法叠前深度偏移方法，包括：CPU获取目标工区的叠前地震数据及偏移参数；该CPU根据该叠前地震数据及偏移参数计算目标工区的激发点走时表和激发点等效Q值表；该CPU根据该激发点走时表和激发点等效Q值表，计算目标工区成像空间的成像空间走时表和成像空间等效Q值表；该成像空间走时表被划分为多个成像线走时表，该成像空间等效Q值表被划分为多个成像线等效Q值表；该CPU运用上述第一方面提供的数据优化方法优化该成像线走时表和成像线等效Q值表，得到并存储成像线优化走时表和成像线优化等效Q值表；该CPU以成像线为单位将成像线优化走时表、成像线优化等效Q值表发送给GPU，并将该叠前地震数据及偏移参数发送给GPU；该GPU线程读取成像线优化走时表和成像线优化等效Q值表，以根据该成像线优化走时表、成像线优化等效Q值表、叠前地震数据及偏移参数计算目标工区成像空间内各个成像点对应不同偏移距的偏移幅值。
本发明实施例带来了以下有益效果：
本发明实施例提供的一种数据优化方法及积分法叠前深度偏移方法，该数据优化方法包括：获取待优化的目标矩阵，该目标矩阵为走时表或等效Q值表；根据该目标矩阵生成第一序列；根据预设的网格密度抽稀第一序列得到第二序列各个元素的取值位置，并基于最小二乘原理求取第二序列各个元素的值；第二序列的元素个数远少于第一序列的元素个数；对第二序列进行插值得到第三序列；该第三序列和第一序列的元素个数相等；计算第三序列对应的目标矩阵；计算待优化的目标矩阵与第三序列对应的目标矩阵之间的误差；当该误差大于预设第一误差阈值时，按预设规则增大网格密度，并继续执行根据预设的网格密度抽稀第一序列得到第二序列各个元素的取值位置的步骤；当该误差小于第一误差阈值时，记录上述第二序列对应的目标矩阵为该待优化的目标矩阵的优化目标矩阵。本发明实施例提供的一种数据优化方法，基于最小二乘原理引入优化的走时表和等效Q值表逼近补偿介质吸收Kirchhoff叠前深度偏移实施流程中的走时表和等效Q值表，在满足计算精度的条件下，大大减少对GPU硬件的存储需求，提高了基于GPU进行补偿介质吸收Kirchhoff叠前深度偏移并行计算的计算效率。
本公开的其他特征和优点将在随后的说明书中阐述，或者，部分特征和优点可以从说明书推知或毫无疑义地确定，或者通过实施本公开的上述技术即可得知。
为使本公开的上述目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附附图，作详细说明如下。
附图说明
为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例提供的一种数据优化方法的流程示意图；
图2为本发明实施例提供的一种成像点走时表或等效Q值表优化示意图；
图3a、图3b、图3c分别为本发明实施例提供的深度3km的成像点走时表优化前慢度、优化后再插值慢度，以及优化前后走时差值的示意图；
图4a、图4b、图4c分别为本发明实施例提供的深度7km的成像点走时表优化前慢度、优化后再插值慢度，以及优化前后走时差值的示意图；
图5为本发明实施例提供的某成像线的走时表与等效Q值表优化前后存储空间对比示意图；
图6为本发明实施例提供的一种积分法叠前深度偏移方法的流程示意图；
图7为本发明实施例提供的另一种积分法叠前深度偏移方法的流程示意图；
图8为本发明实施例提供的一种GPU线程块内多线程读取优化的走时表或等效Q值表数据文件示意图；
图9为本发明实施例提供的某计算节点利用GPU卡运行补偿介质吸收Kirchhoff叠前深度偏移作业未应用优化的走时表和等效Q值表与应用后计算用时对比示意图；
图10a、图10b、图10c、图10d分别为本发明实施例提供的某工区单点高密度采集数据补偿介质吸收Kirchhoff叠前深度偏移处理X方向切片、Y方向切片、Z＝3km切片和Z＝5km切片；
图11为本发明实施例提供的另一种积分法叠前深度偏移方法的流程示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
随着勘探目标由埋深浅、构造简单的油气藏转向埋深大、构造复杂的油气藏或者岩性油气藏，地震勘探也向宽方位、宽频带、高密度三维地震勘探发展。其特点是，在野外地震采集中，采用较宽的方位、高覆盖次数、小面元、单点或小组合接收等系列精细措施。所获取的高空间、时间采样密度资料给地震资料处理技术带来了极大挑战。
作为地震资料处理技术流程中最重要的环节，地震偏移成像技术对寻找复杂构造与岩性圈闭以及油气、流体检测等有重要作用。高空间、时间采样密度资料给地震叠前成像带来的挑战主要存在两个方面：
其一，海量地震数据对计算机硬盘的高存储量需求、高覆盖次数对偏移中成像点偏移幅值计算时的大内存、高计算效率需求；
其二，理论上讲，单点或小组合接收的地震数据具有宽频带的优势，但由于大地在地震波传播过程中的吸收作用，尤其是对埋深大的目标，其地震同相轴的高频成分损失很大，从而造成单点高密度采集地震数据宽频带优势不能体现、地震成像结果分辨率不高的后果。
现行地震叠前偏移方法的实践中，以Kirchhoff方法为代表的积分法叠前偏移技术占据着主体地位。Kirchhoff方法地震叠前偏移主要分为两类，Kirchhoff叠前时间偏移和Kirchhoff叠前深度偏移。
从具体实现来看，Kirchhoff叠前时间偏移对速度模型的处理方式体现了利用双曲近似来描述波的传播过程：对于某一成像点来说，通过假设在偏移孔径范围内存在一个横向均匀的速度模型来对走时进行解析计算，因此其走时的计算简单且高效。而在Kirchhoff叠前深度偏移中，对速度模型的处理方式则体现了一个物理传播的过程：其在走时的计算中所用到的速度信息是偏移孔径范围所确定的层速度模型，这一模型的速度是横向变化的，相应的走时计算必须采用射线追踪这样复杂的方法进行，且计算量大。
人们为了克服Kirchhoff叠前深度偏移计算量大的缺点，通常采用偏移前计算并存储走时表，偏移中自走时表读取相应的走时的策略提高其计算效率。实际应用中，这两类偏移方法的选择与所处理问题的地质特点也是紧密相关的：对反射构造简单或者构造复杂但速度的横向变化不是很剧烈的地质问题，叠前时间偏移是可以较好成像的，但针对速度横向较强变化的地质目标，如埋深大的复杂构造油气藏或者岩性油气藏，则必须采用叠前深度偏移方法。
能够发挥高空间、时间采样密度资料宽频带优势的Kirchhoff叠前时间偏移技术已经有了较充分的研究。通过在常规Kirchhoff叠前时间偏移流程中新引入等效Q值参数，对大地的吸收效应沿地震波的实际传播路径进行分频补偿，即可发挥单点高密度采集数据的宽频带优势，实现对勘探目标的高分辨率成像。这一补偿介质吸收Kirchhoff叠前时间偏移技术中，与对速度模型的处理方法类似，成像点的吸收补偿也仅与该点等效Q值有关。因此，虽然等效Q值的引入带来的分频率偏移幅值计算带来了巨大计算量，我们可以利用GPU并行计算大幅提高其计算效率，使其计算周期控制在可接受的范围。然而如前所述，这一技术在埋深大的复杂构造油气藏或者岩性油气藏勘探中并不适用。
能够发挥高空间、时间采样密度资料宽频带优势的Kirchhoff叠前深度偏移技术实现则存在计算效率低的巨大障碍。这主要与该技术的实现方式和GPU卡的特点有关。为实现大地对地震波的吸收补偿，我们需要在Kirchhoff叠前深度偏移实现流程中引入与地震波传播路径有关的等效Q值，这一等效Q值与地震波自炮点至成像点，再至检波点的传播路径上的所有层Q值有关。因此，这一补偿介质吸收Kirchhoff叠前深度偏移技术中，要得到正确的偏移幅值，除需事先计算并存储走时表外，还需额外存储一个等效Q值表。
然而，对于具有高覆盖次数的高空间、时间采样密度地震数据，其走时表与等效Q值表所需存储量巨大，并且与其覆盖次数是呈正相关的。与之形成鲜明对比的是，GPU卡的主要特点是擅长计算，其存储量非常有限。
目前，采用CPU读取走时表和等效Q值表，然后再通过CPU-GPU之间的通讯解决利用GPU进行偏移幅值计算时的走时表、等效Q值表需求的策略，会大幅降低GPU并行计算的计算效率。因此，走时表和等效Q值表的高效存取是制约补偿介质吸收Kirchhoff叠前深度偏移这一复杂构造油气藏或者岩性油气藏勘探关键技术工业应用的瓶颈。
基于此，本发明实施例提供的一种数据优化方法及积分法叠前深度偏移方法，可以实现以较小的“走时表和等效Q值表”逼近原先较大的“走时表和等效Q值表”，在满足计算精度的条件下，大大减少对GPU的存储需求，使得利用GPU进行补偿介质吸收的Kirchhoff积分法叠前深度偏移技术并行计算的效率大幅提高。并且，通过上述数据优化方法及积分法叠前深度偏移方法，能够高效实现补偿介质吸收Kirchhoff叠前深度偏移技术的工业化应用，充分、快速挖掘高空间、时间采样密度地震资料的宽频带优势，提高这类成本高昂的地震反射数据对地下复杂构造以及岩性油气藏的分辨能力，获得关于地下深层-超深层勘探目标更为精细、准确的构造和流体信息，从而对深层-超深层复杂构造、岩性油气藏的勘探具有重要应用价值。
为便于对本实施例进行理解，首先对本发明实施例所公开的一种数据优化方法进行详细介绍。
实施例一：
如图1所示，为本发明实施例提供的一种数据优化方法的流程示意图，由图1可见，该优化方法包括以下步骤：
步骤S102：获取待优化的目标矩阵，该目标矩阵为走时表或等效Q值表。
这里，待优化的对象称为目标矩阵，该目标矩阵可以是走时表，也可以是等效Q值表，其中，走时表是关于走时的矩阵，等效Q值表是关于等效Q值的矩阵。
以某目标工区的地震数据叠前深度偏移处理为例，具体以成像空间中的成像点的走时表和等效Q值表的优化进行说明，其中，任一成像点的走时表中的每个元素对应为地震波自一个炮点传播至该成像点的走时，并且，该成像点的等效Q值表中的每个元素对应地震波自一个炮点传播至该成像点的等效Q值。
步骤S104：根据该目标矩阵生成第一序列。
当该目标矩阵为走时表时，先将走时表转化为慢度表，再根据该慢度表的所有元素生成第一序列；
当该目标矩阵为等效Q值表时，根据等效Q值表的所有元素生成第一序列。
步骤S106：根据预设的网格密度抽稀第一序列得到第二序列各个元素的取值位置，并基于最小二乘原理求取第二序列各个元素的值；第二序列的元素个数远少于第一序列的元素个数。
这里，利用预设的网格密度对第一序列抽稀，可以得到第二序列各个元素在该第一序列中的取值位置，但是对于第二序列的各个元素的值是未知的。在抽稀第一序列时，通常会设置一个较稀的网格密度，使得抽稀得到的第二序列的元素个数远少于第一序列的元素个数。
对于第二序列中的各个元素，可以根据最小二乘原理求取各个元素的值。具体地，可以以第二序列的所有元素为离散点集，应用插值算法，逼近第一序列各个元素的值，以此确定目标函数，这里，插值算法可以是线性插值、样条插值、克里金插值等；然后，对第二序列中各个元素求偏导数，并令各个元素的偏导数为0，得到关于各个元素的方程组；最后，通过求解该方程组得到各个元素的值。这样，即得到了经过抽稀的第二序列。
步骤S108：对第二序列进行插值得到第三序列；该第三序列和第一序列的元素个数相等。
在得到第二序列之后，利用插值算法扩充序列的元素个数得到第三序列，使得第三序列的元素个数和第一序列的元素个数相等。这里，该插值算法与上述步骤S106中应用的插值算法相同。
为了更好理解上述抽稀和插值的运算，参见图2，为本发明实施例提供的一种成像点走时表或等效Q值表优化示意图，在图2中，中央黑色矩形为成像点在地表投影，圆点为优化后的走时表或等效Q值表所取的点位置，黑色“×”号为优化前走时表或等效Q值表所取的点位置，其中，走时表或等效Q值表优化就是基于最小二乘原理应用圆点处的慢度或等效Q值，应用插值算法逼近黑色“×”号处的慢度或等效Q值。
步骤S110：计算第三序列对应的目标矩阵。
经过上述抽稀得到元素位置和解方程组求取元素值的运算，得到了第二序列。第二序列与第一序列相比较，占用空间很小，但同时也可以通过简单的插值运算逼近第一序列，有利于以该第二序列代替第一序列应用于一些空间复杂度要求严格的算法，例如应用于基于GPU的叠前偏移算法，从这个意义上讲，由目标矩阵的所有元素构成的第一序列得到了“优化”，形成第二序列。
这里，根据第三序列计算其对应的走时表或等效Q值表。
步骤S112：计算待优化的目标矩阵与第三序列对应的目标矩阵之间的误差。
为了评估“优化”后得到的第二序列是否符合精度要求，需要计算基于第二序列应用插值算法得到的第三序列对应的目标矩阵与优化前的目标矩阵之间的误差。
在其中一种实施方式中，上述计算待优化的目标矩阵与第三序列对应的目标矩阵之间的误差的步骤包括：
步骤11：计算第三序列对应的目标矩阵与该待优化的目标矩阵中对应各元素差值的绝对值；
步骤12：统计该差值的绝对值大于预设偏差阈值的元素数量占该待优化的目标矩阵的元素总数的百分比；
步骤13：将该百分比作为第三序列对应的目标矩阵与该待优化的目标矩阵的误差。
步骤S114：比较该误差与预设第一误差阈值的大小。
将计算得到的上述误差与预设的第一误差阈值进行比较，其中，第一误差阈值用于衡量优化后的目标矩阵与优化前目标矩阵之间的差异，为了满足计算精度，需要将优化前后目标矩阵的误差控制在设计的范围内。
步骤S116：当该误差大于预设第一误差阈值时，按预设规则增大网格密度，并继续执行根据预设的网格密度抽稀第一序列得到第二序列各个元素的取值位置的步骤。
这里，若待优化的目标矩阵与第三序列对应的目标矩阵之间的误差大于预设第一误差阈值，表明优化结果尚不满足要求，需要将用于抽稀第一序列的网格密度增大，以增加参与插值的原始数据个数，从而减小拟合插值后的第三序列与第一序列之间的误差。此时，可以按照预设规则增大网格密度，并继续执行根据预设的网格密度抽稀第一序列得到第二序列各个元素的取值位置的步骤，也即，只要误差还大于预设的第一误差阈值，则循环进行上述步骤S106至步骤S114，通过不断增大网格密度，使得上述误差不断逼近预设的第一误差阈值。
步骤S118：当该误差小于第一误差阈值时，记录上述第二序列对应的目标矩阵为该待优化的目标矩阵的优化目标矩阵。
若待优化的目标矩阵与第三序列对应的目标矩阵之间的误差小于预设第一误差阈值，表明优化结果满足要求，此时，记录该第二序列对应的目标矩阵为待优化的目标矩阵的优化目标矩阵。
这样，通过上述方法即可优化走时表和等效Q值表得到相应的优化走时表和优化等效Q值表。
在实际操作中，由于一个工区的地震数据成像空间对应大量的成像点，在对成像点的走时表和等效Q值表进行优化时，为了简化计算过程，在至少一种可能的实施方式中，可以将成像空间划分为多条成像线，将成像线划分成多个成像块，每个成像块里包含多个成像点，其中，对于任一成像块，可以选取成像块内的典型成像点优先进行上述优化运算，在得到典型成像点的优化目标矩阵时，相应记录此次循环用于抽稀序列的网格密度对应的采样个数及间距作为该目标矩阵的共享优化参数。这样，可以将该共享优化参数共享给同一成像块中的其他成像点进行优化运算，其他成像点在优化走时表和等效Q值表时，直接以该网格密度进行序列的抽稀，默认优化结果满足误差要求，从而节省了各个成像点进行循环计算以得到满足误差要求的网格密度的时间。
参见图3a、图3b、图3c，分别为本发明实施例提供的深度3km的成像点走时表优化前慢度、优化后再插值慢度，以及优化前后走时差值的示意图。其中，图3a为优化前根据走时表得到的慢度分布图，其存储占用61×61＝3721浮点类型数值；图3b是在由存储的优化走时表得到的慢度基础上，应用插值算法得到的慢度分布图，走时表优化后仅需存储11×9＝99个短整形数据，走时表优化后存储量低至原来的1.33％；图3c是图3a与图3b所示慢度转换为走时后的差值。
参见图4a、图4b、图4c，分别为本发明实施例提供的深度7km的成像点走时表优化前慢度、优化后再插值慢度，以及优化前后走时差值的示意图。其中，图4a为优化前根据走时表得到的慢度分布图，其存储占用141×141＝19981浮点类型数值；图4b是在由存储的优化走时表得到的慢度基础上应用插值算法得到的慢度分布图，走时表优化后仅需存储39×39＝1521个短整形数据，走时表优化后存储量低至原来的3.83％；图4c是图4a与图图4b所示慢度转换为走时后的差值。
可见，运用本实施例提供的数据优化方法，可以实现以较小的“走时表和等效Q值表”逼近原先较大的“走时表和等效Q值表”，在满足计算精度的条件下，大大减少对GPU的存储需求，这样，可以充分发挥GPU计算的优势，使得利用GPU进行补偿介质吸收的Kirchhoff积分法叠前深度偏移技术并行计算的效率大幅提高。
本发明实施例提供的一种数据优化方法，该方法包括获取待优化的目标矩阵，该目标矩阵为走时表或等效Q值表；根据该目标矩阵生成第一序列；根据预设的网格密度抽稀第一序列得到第二序列各个元素的取值位置，并基于最小二乘原理求取第二序列各个元素的值；第二序列的元素个数远少于第一序列的元素个数；对第二序列进行插值得到第三序列；该第三序列和第一序列的元素个数相等；计算第三序列对应的目标矩阵；计算待优化的目标矩阵与第三序列对应的目标矩阵之间的误差；当该误差大于预设第一误差阈值时，按预设规则增大网格密度，并继续执行根据预设的网格密度抽稀第一序列得到第二序列各个元素的取值位置的步骤；当该误差小于第一误差阈值时，记录上述第二序列对应的目标矩阵为该待优化的目标矩阵的优化目标矩阵。该方法基于最小二乘原理引入优化的走时表和等效Q值表逼近补偿介质吸收Kirchhoff叠前深度偏移实施流程中的走时表和等效Q值表，在满足计算精度的条件下，大大减少对GPU硬件的存储需求，提高了基于GPU进行补偿介质吸收Kirchhoff叠前深度偏移并行计算的计算效率。
实施例二：
为了更好理解上述数据优化方法的过程，以及优化后的效果，这里，以某目标工区地震数据叠前深度偏移处理中，其成像空间的成像点的走时表和等效Q值表的优化进行说明。
对于该目标工区，其激发点走时表为目标工区平面范围内所有炮点走时表的集合τ[nystotal，nxstotal，npy，npx，npz]，其中，nystotal为目标工区炮线数目，nxstotal为目标工区每条炮线包含炮数目，目标工区共计nystotal×nxstotal炮，npy，npx，npz为每一炮覆盖的成像空间在Y，X，Z方向的采样个数。
这里，炮点走时表是指地震波自地面某一炮点(ys，xs，0)激发，传播到该炮地震记录覆盖范围内所有成像点的旅行时集合，记为τ[npy，npx，npz]，其每一个元素定义为：
式中，n为地震波传播到坐标为(yp，xp，zp)的成像点经过的地层数目，Δzi为地层厚度，vi为地层层速度。
另外，该目标工区的激发点等效Q值表指目标工区平面范围内所有炮点等效Q值表的集合Q[nystotal，nxstotal，npy，npx，npz]。
炮点等效Q值表是指地震波自地面某一炮点(ys，xs，0)激发，传播到该炮地震记录覆盖范围内所有成像点的等效Q值的集合，记为Q[npy，npx，npz]，其每一个元素定义为：
式中，n为地震波传播到坐标为(yp，xp，zp)的成像点经过的地层数目，Δτi为地震波在地层内传播旅行时，qi为地层层Q值。
定义该目标工区的成像空间ISPACE[ny，nx，nz]，其中ny，nx，nz，分别为成像空间在y，x，z方向成像点的个数。通常，nz＝npz，ny＞＞npy，nx＞＞npx。
对目标工区的某一炮点，根据炮点横向坐标(ys，xs)、成像点横向坐标(yp，xp)以及该成像点的有效成像孔径范围识别该炮点是否在成像点的有效成像孔径范围内，若是，则根据炮点横向坐标(ys，xs)由目标工区的激发点走时表τ[nystotal，nxstotal，npy，npx，npz]中抽取该炮点至该成像点的走时并由目标工区的激发点等效Q值表Q[nystotal，nxstotal，npy，npx，npz]中抽取地震波沿该炮点至该成像点的传播路径上的等效Q值遍历目标工区的nystotal×nxstotal炮，形成该成像点的走时表和等效Q值表其中，nysP和nxsP分别为该成像点的有效成像孔径范围包含的y方向和x方向炮点个数。
成像空间内所有成像点的走时表集合组成目标工区的成像空间走时表τ[ny，nx，nz，nysP，nxsP]。成像空间内所有成像点的等效Q值表集合组成目标工区的成像空间等效Q值表Q[ny，nx，nz，nysP，nxsP]。
对目标工区的成像空间ISPACE[ny，nx，nz]内的任一成像点(yp，xp，zp)，该成像点的走时表中的每个元素对应地震波自一个炮点(iysP，ixsP)传播至该成像点的走时。
(1)将该成像点的走时表的元素按照下述公式计算：
式中，为该成像点到炮点(iysP，ixsP)的直线距离，形成该成像点的慢度表并由该慢度表的所有元素得到一个包含nysP×nxsP个数的第一慢度序列slowi，i＝1，2，3，…nysP×nxsP。
(2)根据预设的网格密度抽稀上述第一慢度序列，并基于最小二乘原理，求取上述mxsP×mysP个炮点网格上的慢度值组成第二慢度序列i＝1，2，3，…mysP×mxsP；
令目标函数1
最小。式中aki为由该成像点坐标、炮点坐标、以及所应用的插值方法确定的权系数。可得到关于k＝1，2…mxsP×mysP的方程组：
式中，
这里，上述令目标函数1最小也即：对 求偏导数并令之为0。
求解该方程组即可得到mxsP×mysP个炮点网格上的慢度值组成第二慢度序列i＝1，2，3，…mysP×mxsP。
(3)基于第二慢度序列应用插值算法得到原始nysP×nxsP个炮点网格上的第三慢度序列 式中aki与步骤(2)中的权系数相同。
(4)基于步骤(3)得到的原始nysP×nxsP个炮点网格上的第三慢度序列
由该第三慢度序列组成新的慢度矩阵并对其每个元素应用公式得到新的走时矩阵
(5)求取上述新的走时矩阵和原始矩阵的每个元素的差并进行统计分析。
当差的绝对值大于1毫秒的炮(元素)个数大于有效成像孔径范围总炮数(矩阵元素总个数)nysP×nxsP的1％时，按预设规则增加网格密度mxsP、mysP，这里，在其中一种实施方式中，对于上述预设的网格密度mxsP×mysP，其中mxsP、mysP均为奇数，且每次增加网格密度时，x、y方向各增加两个点。
增加该网格密度后，重新回到上述第(2)步，循环计算直至上述差的绝对值大于1毫秒的炮个数小于有效成像孔径范围内总炮数nysP×nxsP的1％，此时，记录达到条件的网格密度对应的采样个数及间距为走时表共享优化参数，具体地，记录采样个数mxsP为记录mysP为同时记录此时走时表炮点网格的间距这里，即为走时表共享优化参数。
相似地，对于上述成像点(yp，xp，zp)的等效Q值表的优化过程如下。
其中，该成像点的等效Q值表中的每个元素对应地震波自一个炮点(iysP，ixsP)传播至该成像点的等效Q值。
(6)将该成像点的等效Q值表的所有元素组成一个包含nysP×nxsP个数的第一等效Q值序列Qi，i＝1，2，3，…nysP×nxsP。
(7)根据预设的网格密度抽稀上述第一等效Q值序列，并基于最小二乘原理，求取上述mxsP×mysP个炮点网格上的等效Q值组成第二等效Q值序列i＝1，2，3，…mysP×mxsP；
令目标函数2：
最小。式中aki为由该成像点坐标、炮点坐标、以及所应用的插值方法确定的权系数。可得到关于k＝1，2…mxsP×mysP的方程组：
式中，
这里，令上述目标函数2最小也即：对k＝1，2…mxsP×mysP求偏导数，并令之为0。
求解该方程组即可得到坐标为(yp，xp，zp)的成像点处mxsP×mysP个炮点网格上的等效Q值组成的第二等效Q值序列i＝1，2，3，…mysP×mxsP。
(8)基于该第二等效Q值序列应用插值算法得到原始nysP×nxsP个炮点网格上的第三等效Q值序列i＝1，2，3，…nysP×nxsP。式中aki与步骤(7)中的权系数相同。
(9)基于步骤(8)得到的原始nysP×nxsP个炮点网格上的等效Q值序列i＝1，2，3，…nysP×nxsP组成新的等效Q值矩阵
(10)求取该新的等效Q值矩阵和原始矩阵的每个元素的差，并进行统计分析；
当差的绝对值大于5的炮(元素)个数，大于有效成像孔径范围总炮数(矩阵元素总个数)nysP×nxsP的1％时，按预设规则增加mxsP×mysP中的mxsP、mysP的值，并重新回到第(7)步，循环计算直至上述差的绝对值大于5的炮个数小于有效成像孔径范围内总炮数nysP×nxsP的1％，记录达到条件的网格密度对应的采样个数及间距为等效Q值表共享优化参数，具体地，记录采样个数mxsP为记录mysP为同时记录此时走时表炮点网格的间距这里，即为等效Q值表的共享优化参数。
这样，通过上述步骤即可得到目标工区的成像空间ISPACE[ny，nx，nz]内的任一成像点的优化走时表和优化等效Q值表。
在实际操作中，为了简化计算过程，可以将成像空间划分为多条成像线，将成像线划分成多个成像块，每个成像块里包含多个成像点，其中，对于任一成像块，可以选取成像块内的典型成像点优先进行上述优化运算，在得到典型成像点的优化走时表和优化等效Q值表时，相应记录该典型成像点的走时表共享优化参数和等效Q值表共享优化参数。这样，可以将共享优化参数共享给同一成像块中的其他成像点进行优化运算，以节省成像块中其他成像点进行循环计算以得到满足误差要求的网格密度的时间。具体地，操作如下：
首先，根据目标工区的成像空间ISPACE[ny，nx，nz]内包含的成像点y坐标，将该目标工区的成像空间走时表τ[ny，nx，nz，nysP，nxsP]划分为不同成像线的走时表集合τiy[nx，nz，nysP，nxsP]，iy＝1，2，3…ny。将该目标工区的成像空间等效Q值表Q[ny，nx，nz，nysP，nxsP]划分为不同成像线的等效Q值表集合Qiy[nx，nz，nysP，nxsP]，iy＝1，2，3…ny。每条成像线包含nx×nz个成像点。
其次，对每条成像线，根据其所包含成像点上的速度在x方向和z方向上的横向和纵向变化情况，将其划分为mxb×mzb个成像块，每个成像块在x方向采样个数为nxb，在z方向采样个数为nzb。
然后，对目标工区的成像空间ISPACE[ny，nx，nz]内某一成像线的某一成像块，依据条件1选取成像块内的典型成像点，其中，条件1为：成像点的z坐标最大且速度横向梯度最大。
对于任一成像块，选取其典型成像点后，通过上述(1)至(10)的步骤获得该典型成像点的走时表共享优化参数和以及其等效Q值表共享优化参数和
这里，典型成像点所在的成像块内共享相同的走时表优化参数和令对该典型成像点所在的成像块内的任意一个成像点应用步骤(1)-(2)得到该成像块的优化慢度序列式中，(yp，xp，zp)为成像块内成像点的空间坐标。再应用 得到各个成像点的优化走时表式中，ri(yp，xp，zp)为第i个炮点到成像点(yp，xp，zp)的直线距离。
遍历上述典型成像点所在的成像块内的所有成像点，得到该成像块的优化走时表
同理，对上述典型成像点所在的成像块内共享相同的等效Q值表优化参数和令对该典型成像点所在的成像块内的任意一个成像点应用步骤(6)-(7)得到该成像点的优化等效Q值序列并形成该成像点的优化等效Q值表式中，(yp，xp，zp)为该成像点的空间坐标。
遍历该典型成像点所在的成像块内的所有成像点，得到该成像块的优化等效Q值表
同理，遍历上述成像线的所有成像块，得到该成像线的优化走时表和等效Q值表。
并且，遍历目标工区的成像空间ISPACE[ny，nx，nz]内所有的成像线，即可得到该目标工区的成像空间ISPACE[ny，nx，nz]的优化走时表和等效Q值表。
参见图5，为本发明实施例提供的某成像线的走时表与等效Q值表优化前后存储空间对比示意图，其中，该成像线在X方向有1124个采样点，Z方向有2000个采样点，共包含2248000个成像点，优化前走时表的大小为72.03G，优化后走时表的大小缩减至2.2G，优化前等效Q值表的大小为72.03G，优化后缩减至0.79G。
由图5示出的效果可见，通过本发明实施例提供的走时表与等效Q值表优化方法，可以大大降低成像线的走时表与等效Q值表对存储空间的要求。
实施例三：
本发明实施例还提供了一种积分法叠前深度偏移方法，该方法应用于CPU，参见图6，为该方法的流程示意图，由图6可见，该方法包括以下步骤：
步骤S602：获取目标工区的叠前地震数据、偏移参数，以及该目标工区成像空间的成像线走时表和成像线等效Q值表。
这里，目标工区的成像空间被划分为多条成像线，该成像线被划分为多个成像块，该成像块包含多个成像点。
步骤S604：运用实施例一提供的数据优化方法优化该成像线走时表和该成像线等效Q值表得到成像线优化走时表和成像线优化等效Q值表。
在CPU中以成像线为单位，逐条优化成像线对应的成像线走时表和该成像线等效Q值表，得的各条成像线对应的成像线优化走时表和成像线优化等效Q值表。这里，对于成像线走时表和成像线等效Q值表的优化步骤和过程，可以参见实施例二中的相应内容，在此不再赘述。
步骤S606：以成像线为单位存储该成像线优化走时表和该成像线优化等效Q值表。
在CPU中将优化后得到的各条成像线对应的成像线优化走时表和该成像线优化等效Q值表，以成像线为单位进行存储。
具体地，上述以成像线为单位存储该成像线优化走时表和该成像线优化等效Q值表的步骤，包括：
首先，以成像线为单位构建成像线优化走时表的第一索引文件和成像线优化等效Q值表的第二索引文件。其中，该第一索引文件包括成像线的成像块划分参数，每个成像块内的共享走时表优化参数，以及每个成像块内第一个成像点的优化走时表在成像线优化走时表的数据文件中的存储位置。并且，该第二索引文件包括成像块划分参数，每个成像块内的共享等效Q值表优化参数，以及每个成像块内第一个成像点的优化等效Q值在成像线优化等效Q值表的数据文件中的存储位置。
其次，存储第一索引文件和第二索引文件，并以无符号短整型变量类型存储成像线优化走时表的数据文件和成像线优化等效Q值表的数据文件。
这里，根据第一关系式、第二关系式分别逐成像块存储该成像线优化走时表的数据文件、该成像线优化等效Q值表的数据文件。第一关系式为：
第二关系式为：
式中，ixi为成像点在成像块内位置x方向索引，izi为成像点在成像块内位置z方向索引，为该成像点优化走时表中炮点在x方向索引，为该成像点优化走时表中炮点在y方向索引；为所述成像点优化等效Q值表中炮点在x方向索引，为所述成像点优化等效Q值表中炮点在y方向索引；k为maxτ为该目标工区的成像空间优化走时表中走时的最大值，b为0.0；表示该成像块优化走时表数据文件的存储顺序；表示该成像块优化走时表的数据格式；表示该成像块优化等效Q值表的数据格式；表示该成像块优化等效Q值表数据文件的存储顺序。
这里，在对数据文件进行存储时，针对某一成像块，存在一个数组顺序重排的操作，是为了提高GPU线程对走时表和等效Q值表数据文件的读取速度，而依据第一、第二关系式以无符号短整形数据格式存储优化走时表和优化等效Q值表的存储策略可在满足计算精度要求的情况下将GPU卡的内存需求减少一半。
当GPU进行积分法叠前深度偏移计算时，GPU以成像线为单位读取上述成像线优化走时表和成像线优化等效Q值表。
步骤S608:将叠前地震数据、偏移参数，以及成像线优化走时表和成像线优化等效Q值表以成像线为单位发送给GPU，以根据该成像线优化走时表、该成像线优化等效Q值表、该叠前地震数据及该偏移参数，计算目标工区成像空间内各个成像点对应不同偏移距的偏移幅值。
在进行叠前偏移计算时，CPU首先将叠前地震数据、偏移参数发送给GPU，然后，以成像线为单位将成像线优化走时表和成像线优化等效Q值表发送给GPU，在其中一种实施方式中，GPU中一次只接受一条成像线的成像线优化走时表和成像线优化等效Q值表，并根据该成像线优化等效Q值表、该叠前地震数据及该偏移参数计算该成像线中各个成像点对应不同偏移距的偏移幅值，在计算完该条成像线后，再接收下一条成像线对应的成像线优化走时表和成像线优化等效Q值表以进行新一轮的计算，如此往复，直到计算完目标工区成像空间内各个成像点对应不同偏移距的偏移幅值。
参见图7，为本发明实施例提供的另一种积分法叠前深度偏移方法，该方法应用于GPU，由图7可见，该方法包括以下步骤：
步骤S702：接收CPU发送的目标工区的叠前地震数据、偏移参数，以及该CPU发送的以成像线为单位的，目标工区成像空间的成像线优化走时表和成像线优化等效Q值表；该成像空间被划分为多条成像线，该成像线被划分为多个成像块，该成像块包含多个成像点；该成像线优化走时表和该成像线优化等效Q值表，为上述CPU运用第一方面提供的数据优化方法优化该目标工区的成像线走时表和成像线等效Q值表得到。
在利用本方法进行叠前偏移计算时，在其中一种可能的实施方式中，由CPU和GPU协同并行计算，此时，GPU可以从CPU获取目标工区的叠前地震数据及偏移参数，并且，该GPU一次只接受一条成像线的成像线优化走时表和成像线优化等效Q值表，并根据该成像线优化等效Q值表、该叠前地震数据及该偏移参数计算该成像线中各个成像点对应不同偏移距的偏移幅值，在计算完该条成像线后，再接收下一条成像线对应的成像线优化走时表和成像线优化等效Q值表以进行新一轮的计算。
对于上述目标工区的成像空间，该成像空间被划分为多条成像线，该成像线被划分为多个成像块，该成像块包含多个成像点。这里，对于成像线走时表和成像线等效Q值表的优化步骤和过程，可以参见实施例二中的相应内容，在此不再赘述。
步骤S704：根据成像线优化走时表、成像线优化等效Q值表、叠前地震数据及偏移参数，计算目标工区成像空间内各个成像点对应不同偏移距的偏移幅值。
在GPU以成像线为单位读取了待计算的成像线优化走时表和成像线优化等效Q值表之后，结合叠前地震数据及偏移参数即可计算各个成像点对应不同偏移距的偏移幅值，在至少一种可能的实施方式中，上述计算偏移幅值的步骤包括：
步骤31：获取叠前地震道的炮点坐标和检波点坐标。
步骤32：确定GPU线程块的大小和GPU线程块网格的大小。
步骤33：根据该GPU线程块的大小和该GPU线程块网格的大小计算该成像线中每个成像点的线程索引，该线程索引包括深度位置索引和横向位置索引。
步骤34：根据该线程索引线程读取成像线的第一索引文件和第二索引文件，以确定每个成像点所在的成像块的走时表优化参数和等效Q值表优化参数，以及该成像点所在成像块内第一个成像点的优化走时表在成像线优化走时表的数据文件中的存储位置，和该成像点所在成像块内第一个成像点的优化等效Q值表在成像线优化等效Q值表的数据文件中的存储位置。
步骤35：根据该炮点坐标、该检波点坐标，以及该成像线对应的第一索引文件、第二索引文件、成像线优化走时表的数据文件和成像线优化等效Q值表的数据文件，线程获取叠前地震道对应的炮点到该成像点的走时、等效Q值，叠前地震道对应的检波点到该成像点的走时、等效Q值。
这里，线程获取该叠前地震道对应的炮点到该成像点的走时和该叠前地震道对应的检波点到该成像点的走时的步骤，包括：
根据该炮点坐标、检波点坐标、走时表优化参数及走时表炮点网格间距，计算该炮点坐标对应的优化走时表的第一网格及该第一网格四个顶点的第一索引值，以及该检波点坐标对应的优化走时表的第二网格及该第二网格四个顶点的第二索引值；
根据该第一索引值、第二索引值，分别计算该第一网格四个顶点、该第二网格四个顶点的走时在成像线优化走时表的数据文件中的第一位置、第二位置，并根据该第一位置、第二位置读取第一网格四个顶点的走时、第二网格四个顶点的走时；
根据该第一网格四个顶点的走时、第二网格四个顶点的走时分别计算第一网格四个顶点的慢度值、第二网格四个顶点的慢度值；根据该第一网格四个顶点的慢度值、第二网格四个顶点的慢度值分别插值得到该成像点的第一慢度值、第二慢度值；
根据该第一慢度值、第二慢度值分别计算该叠前地震道对应的炮点到该成像点的走时、该叠前地震道对应的检波点到该成像点的走时。
另外，线程获取该叠前地震道对应的炮点到该成像点的等效Q值和该叠前地震道对应的检波点到该成像点的等效Q值的步骤，包括：
根据该炮点坐标、检波点坐标、等效Q值优化参数及等效Q值表炮点网格间距，计算该炮点坐标对应的优化等效Q值的第三网格及该第三网格四个顶点的第三索引值，以及该检波点坐标对应的优化等效Q值的第四网格及该第四网格四个顶点的第四索引值；
根据该第三索引值、第四索引值，分别计算该第三网格四个顶点、第四网格四个顶点的等效Q值在成像线优化等效Q值表的数据文件中的第三位置、第四位置，并根据该第三位置、第四位置读取该第三网格四个顶点的等效Q值、第四网格四个顶点的等效Q值；
根据该第三网格四个顶点的等效Q值、第四网格四个顶点的等效Q值，分别插值得到该叠前地震道对应的炮点到该成像点的等效Q值、该叠前地震道对应的检波点到该成像点的等效Q值。
步骤36：根据成像公式线程计算各个成像点的偏移幅值。该成像公式为：
式中，I(h，xp，yp，zp)为该成像点的偏移幅值，h为该叠前地震道偏移距，j为虚部单位，ω为频率，ω0为地震数据主频，F(ω)为该叠前地震道变换至频率域的序列，rs，rg分别为该叠前地震道对应的炮点、检波点至该成像点的距离，τs、Qs分别为该叠前地震道对应的炮点到该成像点的走时、等效Q值，τg、Qg分别为该叠前地震道对应的检波点到该成像点的走时、等效Q值。
本实施例提供的积分法叠前深度偏移方法，将成像空间进行划分，同一成像线内划分为多个成像块，不同成像块的优化参数可以不同，即减少了成像线走时表和等效Q值表的存储空间，又可适应成像线区域内构造非均匀变化。并且，通过引入短整形数据存储优化的走时表，可在满足走时存储精度要求的前提下进一步减少GPU存储需求；通过引入索引文件及数据文件相结合的走时表和等效Q值表存储方式，以及优化数据文件中走时表和等效Q值表的存储顺序，有利于GPU线程实现对走时及等效Q值的快速读取。
实施例四：
为了更好理解上述实施例三提供的积分法叠前深度偏移方法中GPU线程块内多线程读取优化的走时表和等效Q值表以及偏移幅值计算的过程，本实施例仍以上述实施例二中的目标工区成像空间为例进行说明。
对于任一个计算机节点，该节点包含CPU和GPU，其中，在该CPU中，目标工区的成像空间优化走时表和等效Q值表采用分成像线存储的方式。每条成像线的优化走时表和等效Q值表，均由1个索引文件和1个数据文件构成。
成像线的优化走时表索引文件描述了成像线成像区域的分块情况，包括成像线区域成像块划分参数mxb、mzb、nxb、nzb；每个成像块内的共享走时表优化参数以及每个成像块内第一个成像点的优化走时表在成像线优化走时表的数据文件中的存储位置
成像线的优化等效Q值表索引文件也描述了成像线成像区域的分块情况，包括成像线区域成像块划分参数mxb、mzb、nxb、nzb；每个成像块内的共享走时表优化参数以及每个成像块内第一个成像点的优化等效Q值表在成像线优化等效Q值表的数据文件中的存储位置
并且，成像线的优化走时表和优化等效Q值表的数据文件以无符号短整型变量类型存储。
优化走时表数据文件按照上述实施例三中的第一关系式存储地震波自优化后得到的个炮点网格上的每一炮点至成像点的走时。
优化等效Q值表数据文件按照上述实施例三中的第二关系式直接存储地震波自优化后得到的个炮点网格上的每一炮点至成像点的等效Q值。
另外，在GPU计算成像点对应不同偏移距的偏移幅值时，其步骤如下：
步骤41：CPU将该目标工区成像空间内各个成像点按照其y坐标划分为不同的成像线。
步骤42：CPU根据该成像线的y方向坐标自该目标工区的成像空间优化走时表和等效Q值表中确定该成像线对应的成像线的优化走时表和优化等效Q值表存储文件。
步骤43：CPU读取叠前地震道，获取其炮点坐标sx，sy、检波点坐标gx，gy以及地震数据并发送至GPU。
步骤44：CPU读取成像线的优化走时表和优化等效Q值表存储文件并发送至GPU。
步骤45：确定GPU线程块的大小：blockDim.x＝32blockDim.y＝16
根据成像线在x方向和z方向的采样个数nx，nz确定GPU线程块网格的大小：
步骤46：GPU确定成像线中每个成像点的线程索引：
深度位置索引izi＝blockIdx.x×blockDim.x+threadIdx.x
横向位置索引ixi＝blockIdx.y×blockDim.y+threadIdx.y
步骤47：GPU每个线程读取成像线的优化走时表和优化等效Q值表索引文件，根据该成像点(yp，xp，zp)的x方向坐标、z方向坐标确定该成像点所在的成像块、走时表优化参数和该成像点所在成像块内第一个成像点的优化走时表在成像线的优化走时表数据文件中的存储位置等效Q值表优化参数和该成像点所在成像块内第一个成像点的优化等效Q值表在成像线的优化等效Q值表数据文件中的存储位置这里，可参见图8，为GPU线程块内多线程读取优化的走时表或等效Q值表数据文件示意图。
步骤48：GPU每个线程根据该叠前地震道炮点坐标sx，sy，检波点坐标gx，gy以及该成像线对应的成像线的优化走时表和优化等效Q值表存储文件，分别获取该叠前地震道对应的炮点到该成像点的走时τs、等效Q值Qs，该叠前地震道对应的检波点到该成像点到的走时τg、等效Q值Qg。
具体地，计算叠前地震道对应的炮点到该成像点的走时τs、检波点到该成像点的走时τg的步骤包括：
步骤481：由该叠前地震道炮点坐标sx，sy及优化走时表参数间距得到该叠前地震道炮点坐标sx，sy对应的优化走时表网格及该网格的四个顶点的索引值：
(ixs1，iys1)，(ixs2，iys2)，(ixs3，iys3)，(ixs4，iys4)
对于(ixsj，iysj)，j＝1，2，3，4计算其走时在成像线的优化走时表数据存储文件中的位置：
并读取相应走时j＝1，2，3，4。
步骤482：利用公式j＝1，2，3，4，k为maxτ为所述目标工区的成像空间优化走时表中走时的最大值，b为0.0，计算该网格的四个顶点的慢度值并应用上述实施例一中插值方法得到该成像点的慢度值Ss，对其求倒数得到该叠前地震道对应的炮点到该成像点的走时τs。
同理，将上述步骤481～482中的sx更换为gx，sy更换为gy，即可得到该叠前地震道对应的检波点到该成像点的走时τg。
另外，计算该叠前地震道对应的炮点到该成像点的等效Q值Qs、检波点到该成像点到等效Q值Qg的步骤包括：
步骤483：由该叠前地震道炮点坐标sx，sy及优化等效Q值表参数 间距得到该叠前地震道炮点坐标sx，sy对应的优化等效Q值表网格及该网格的四个顶点的索引值：
对于j＝1，2，3，4计算其等效Q值在成像线的优化等效Q值表数据存储文件中的位置：
并读取相应等效Q值j＝1，2，3，4。
步骤484：应用上述实施例一中插值方法得到该叠前地震道对应的炮点到该成像点传播路径的等效Q值Qs。
同理，将上述步骤483～484中的sx更换为gx，sy更换为gy，即可得到该叠前地震道对应的检波点到该成像点传播路径的等效Q值Qg。
步骤49：GPU每个线程按上述实施例三中的成像公式计算成像点的偏移幅值。
这里，参见图9为本发明实施例提供的某计算节点利用GPU卡运行补偿介质吸收Kirchhoff叠前深度偏移作业未应用优化的走时表和等效Q值表与应用后计算用时对比示意图，其中，该计算节点分配的输入地震叠前数据大小为194G。在图9示出的结果中，走时表和等效Q值表优化前计算用时为14.58小时，走时表和等效Q值表优化后计算用时为2.32小时，应用走时表和等效Q值表优化后，计算效率提高了6.29倍。
另外，参见图10a、图10b、图10c和图10d，分别为本发明实施例提供的某工区单点高密度采集数据补偿介质吸收Kirchhoff叠前深度偏移处理X方向切片、Y方向切片、Z＝3km切片和Z＝5km切片。
可见，通过本实施例提供的积分法叠前深度偏移方法，大幅降低了对GPU硬件的存储需求，在保证计算精度的前提下，有效提高了基于GPU进行补偿介质吸收Kirchhoff叠前深度偏移并行计算的计算效率。
实施例五：
参见图11，为本发明实施例还提供的另一种积分法叠前深度偏移方法，由图10可见，该方法包括以下步骤：
步骤1002：CPU获取目标工区的叠前地震数据及偏移参数。
对于单个计算节点，该计算节点包含CPU和GPU，在其中一种实施方式中，计算节点可以从集群计算机主节点获取分配到该计算节点的叠前地震数据及偏移参数。
步骤1004：该CPU根据该叠前地震数据及偏移参数计算目标工区的激发点走时表和激发点等效Q值表。
在至少一种可能的实施方式中，集群计算机主节点根据目标工区的平面范围分配激发点走时表和等效Q值表计算任务至各计算节点，集群计算机各计算节点利用CPU，应用快速波前法计算本计算节点分配的平面范围内的激发点走时表和等效Q值表。并且，最后由该主节点回收并形成该目标工区的激发点走时表和等效Q值表。
步骤1006：该CPU根据该激发点走时表和激发点等效Q值表，计算目标工区成像空间的成像空间走时表和成像空间等效Q值表。该成像空间走时表被划分为多个成像线走时表，该成像空间等效Q值表被划分为多个成像线等效Q值表。
这里，目标工区成像空间的成像空间走时表和成像空间等效Q值表的计算过程可参见上述实施例二中的有关内容，在此不再赘述。
步骤1008：该CPU运用上述实施例一提供的数据优化方法优化该成像线走时表和成像线等效Q值表，得到成像线优化走时表和成像线优化等效Q值表。
步骤1010：该CPU以成像线为单位将成像线优化走时表、成像线优化等效Q值表发送给GPU，并将该叠前地震数据及偏移参数发送给GPU。
步骤1012：GPU接收叠前地震数据及偏移参数，且以成像线为单位接收成像线优化走时表和优化等效Q值表，并将其存储在全局内存中。
步骤1014：该GPU线程读取成像线优化走时表和成像线优化等效Q值表，以根据该成像线优化走时表、成像线优化等效Q值表、叠前地震数据及偏移参数计算目标工区成像空间内各个成像点对应不同偏移距的偏移幅值。
在另一种实施方式中，在上述步骤S1014之后，该积分法叠前深度偏移方法还包括：
GPU以成像线为单位将目标工区成像空间内各个成像点对应不同偏移距的幅值发送给CPU；
CPU以成像线为单位接收来自GPU的目标工区成像空间内各个成像点对应不同偏移距的幅值，并将之输出到持久化存储介质。
本实施例五介绍了CPU和GPU协同并行计算实现积分法叠前深度偏移计算的过程，其中，CPU和GPU的详细工作过程可参看上述实施例一至实施例四的相应内容，在此不再赘述。
本发明实施例提供的积分法叠前深度偏移方法，与上述实施例一、实施例二提供的数据优化方法具有相同的技术特征，所以也能解决相同的技术问题，达到相同的技术效果。
除非另外具体说明，否则在这些实施例中阐述的部件和步骤的相对步骤、数字表达式和数值并不限制本发明的范围。
在这里示出和描述的所有示例中，任何具体值应被解释为仅仅是示例性的，而不是作为限制，因此，示例性实施例的其他示例可以具有不同的值。
附图中的流程图和框图显示了根据本发明的多个实施例的***、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段或代码的一部分，所述模块、程序段或代码的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。也应当注意，在有些作为替换的实现中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个连续的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或动作的专用的基于硬件的***来实现，或者可以用专用硬件与计算机指令的组合来实现。
另外，在本发明实施例的描述中，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。
在本发明的描述中，需要说明的是，术语“中心”、“上”、“下”、“左”、“右”、“竖直”、“水平”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”、“第三”仅用于描述目的，而不能理解为指示或暗示相对重要性。
本发明实施例所提供的进行数据优化方法的计算机程序产品，包括存储了处理器可执行的非易失的程序代码的计算机可读存储介质，所述程序代码包括的指令可用于执行前面方法实施例中所述的方法，具体实现可参见方法实施例，在此不再赘述。
所述功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个处理器可执行的非易失的计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。
最后应说明的是：以上所述实施例，仅为本发明的具体实施方式，用以说明本发明的技术方案，而非对其限制，本发明的保护范围并不局限于此，尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，其依然可以对前述实施例所记载的技术方案进行修改或可轻易想到变化，或者对其中部分技术特征进行等同替换；而这些修改、变化或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。

Claims (13)

1.一种数据优化方法，其特征在于，包括：
获取待优化的目标矩阵，所述目标矩阵为走时表或等效Q值表；
根据所述目标矩阵生成第一序列；
根据预设的网格密度抽稀所述第一序列，得到第二序列各个元素的取值位置，并基于最小二乘原理求取所述第二序列各个元素的值；所述第二序列的元素个数远少于所述第一序列的元素个数；
对所述第二序列进行插值得到第三序列；所述第三序列和所述第一序列的元素个数相等；
计算所述第三序列对应的目标矩阵；
计算所述待优化的目标矩阵与所述第三序列对应的目标矩阵之间的误差；
比较所述误差与预设第一误差阈值的大小，当所述误差大于所述预设第一误差阈值时，按预设规则增大所述网格密度，并继续执行根据预设的网格密度抽稀所述第一序列得到第二序列各个元素的取值位置的步骤；
当所述误差小于所述第一误差阈值时，记录所述第二序列对应的目标矩阵为所述待优化的目标矩阵的优化目标矩阵。
2.根据权利要求1所述的数据优化方法，其特征在于，所述根据所述目标矩阵生成第一序列的步骤，包括：
当所述目标矩阵为走时表时，先将所述走时表转化为慢度表，再根据所述慢度表的所有元素生成第一序列；
当所述目标矩阵为等效Q值表时，根据所述等效Q值表的所有元素生成第一序列。
3.根据权利要求1所述的数据优化方法，其特征在于，所述基于最小二乘原理求取所述第二序列各个元素的值的步骤，包括：
对所述第二序列中各个元素求偏导数，并令所述各个元素的偏导数为0，得到关于所述各个元素的方程组；
求解所述方程组得到所述各个元素的值。
4.根据权利要求1所述的数据优化方法，其特征在于，所述计算所述待优化的目标矩阵与所述第三序列对应的目标矩阵之间的误差的步骤，包括：
计算所述第三序列对应的目标矩阵与所述待优化的目标矩阵中对应各元素差值的绝对值；
统计所述差值的绝对值大于预设偏差阈值的元素数量占所述待优化的目标矩阵的元素总数的百分比；
将所述百分比作为所述第三序列对应的目标矩阵与所述待优化的目标矩阵的误差。
5.根据权利要求1所述的数据优化方法，其特征在于，在所述当所述误差小于所述第一误差阈值时，记录所述第二序列对应的目标矩阵为所述待优化的目标矩阵的优化目标矩阵的步骤之后，所述方法还包括：
记录所述网格密度，将所述网格密度对应的采样个数及间距作为所述目标矩阵的共享优化参数。
6.一种积分法叠前深度偏移方法，其特征在于，应用于CPU，所述方法包括：
获取目标工区的叠前地震数据、偏移参数，以及所述目标工区成像空间的成像线走时表和成像线等效Q值表；
运用权利要求1所述的数据优化方法优化所述成像线走时表和所述成像线等效Q值表得到成像线优化走时表和成像线优化等效Q值表；
以成像线为单位存储所述成像线优化走时表和所述成像线优化等效Q值表；
将所述叠前地震数据、所述偏移参数，以及所述成像线优化走时表和所述成像线优化等效Q值表以成像线为单位发送给GPU，以根据所述成像线优化走时表、所述成像线优化等效Q值表、所述叠前地震数据及所述偏移参数，计算所述目标工区成像空间内各个成像点对应不同偏移距的偏移幅值。
7.根据权利要求6所述的积分法叠前深度偏移方法，其特征在于，所述以成像线为单位存储所述成像线优化走时表和所述成像线优化等效Q值表的步骤，包括：
以成像线为单位构建成像线优化走时表的第一索引文件和成像线优化等效Q值表的第二索引文件；所述第一索引文件包括成像线的成像块划分参数，每个所述成像块内的共享走时表优化参数，以及每个所述成像块内第一个成像点的优化走时表在成像线优化走时表的数据文件中的存储位置；所述第二索引文件包括所述成像块划分参数，每个所述成像块内的共享等效Q值表优化参数，以及每个所述成像块内第一个成像点的优化等效Q值在成像线优化等效Q值表的数据文件中的存储位置；
存储所述第一索引文件和所述第二索引文件，并以无符号短整型变量类型存储所述成像线优化走时表的数据文件和所述成像线优化等效Q值表的数据文件。
8.根据权利要求7所述的积分法叠前深度偏移方法，其特征在于，所述以无符号短整型变量类型存储所述成像线优化走时表的数据文件和所述成像线优化等效Q值表的数据文件的步骤，包括：
根据第一关系式、第二关系式分别逐成像块存储所述成像线优化走时表的数据文件、所述成像线优化等效Q值表的数据文件；
所述第一关系式为：
所述第二关系式为：
式中，ixi为成像点在成像块内位置x方向索引，izi为成像点在成像块内位置z方向索引，为所述成像点优化走时表中炮点在x方向索引，为所述成像点优化走时表中炮点在y方向索引；为所述成像点优化等效Q值表中炮点在x方向索引，为所述成像点优化等效Q值表中炮点在y方向索引；k为maxτ为所述目标工区的成像空间优化走时表中走时的最大值，b为0.0；表示所述成像块优化走时表数据文件的存储顺序；表示所述成像块优化走时表的数据格式；表示所述成像块优化等效Q值表的数据格式；表示所述成像块优化等效Q值表数据文件的存储顺序。
9.一种积分法叠前深度偏移方法，其特征在于，应用于GPU，所述方法包括：
接收CPU发送的目标工区的叠前地震数据、偏移参数，以及所述CPU发送的以成像线为单位的，所述目标工区成像空间的成像线优化走时表和成像线优化等效Q值表；所述成像空间被划分为多条成像线，所述成像线被划分为多个成像块，所述成像块包含多个成像点；所述成像线优化走时表和所述成像线优化等效Q值表，为所述CPU运用权利要求1所述的数据优化方法优化所述目标工区的成像线走时表和成像线等效Q值表得到；
根据所述成像线优化走时表、所述成像线优化等效Q值表、所述叠前地震数据及所述偏移参数，计算所述目标工区成像空间内各个成像点对应不同偏移距的偏移幅值。
10.根据权利要求9所述的积分法叠前深度偏移方法，其特征在于，所述根据所述成像线优化走时表、所述成像线优化等效Q值表、所述叠前地震数据及所述偏移参数，计算所述目标工区成像空间内各个成像点对应不同偏移距的偏移幅值的步骤，包括：
获取叠前地震道的炮点坐标和检波点坐标；
确定GPU线程块的大小和GPU线程块网格的大小；
根据所述GPU线程块的大小和所述GPU线程块网格的大小计算所述成像线中每个成像点的线程索引，所述线程索引包括深度位置索引和横向位置索引；
根据所述线程索引线程读取成像线的第一索引文件和第二索引文件，以确定每个所述成像点所在的成像块的走时表优化参数和等效Q值表优化参数，以及所述成像点所在成像块内第一个成像点的优化走时表在成像线优化走时表的数据文件中的存储位置，和所述成像点所在成像块内第一个成像点的优化等效Q值表在成像线优化等效Q值表的数据文件中的存储位置；
根据所述炮点坐标、所述检波点坐标，以及所述成像线对应的所述第一索引文件、所述第二索引文件、所述成像线优化走时表的数据文件和所述成像线优化等效Q值表的数据文件，线程获取所述叠前地震道对应的炮点到所述成像点的走时、等效Q值，所述叠前地震道对应的检波点到所述成像点的走时、等效Q值；
根据成像公式线程计算各个成像点的偏移幅值；所述成像公式为：
式中，I(h，xp，yp，2p)为所述成像点的偏移幅值，h为所述叠前地震道偏移距，j为虚部单位，ω为频率，ω0为地震数据主频，F(ω)为所述叠前地震道变换至频率域的序列，rs，rg分别为所述叠前地震道对应的炮点、检波点至所述成像点的距离，τs、Qs分别为所述叠前地震道对应的炮点到所述成像点的走时、等效Q值，τg、Qg分别为所述叠前地震道对应的检波点到所述成像点的走时、等效Q值。
11.根据权利要求10所述的积分法叠前深度偏移方法，其特征在于，所述线程获取所述叠前地震道对应的炮点到所述成像点的走时，所述叠前地震道对应的检波点到所述成像点的走时的步骤，包括：
根据所述炮点坐标、所述检波点坐标和走时表优化参数，计算所述炮点坐标对应的优化走时表的第一网格及所述第一网格四个顶点的第一索引值，以及所述检波点坐标对应的优化走时表的第二网格及所述第二网格四个顶点的第二索引值；
根据所述第一索引值、所述第二索引值，分别计算所述第一网格四个顶点、所述第二网格四个顶点的走时在成像线优化走时表的数据文件中的第一位置、第二位置，并根据所述第一位置、所述第二位置读取所述第一网格四个顶点的走时、所述第二网格四个顶点的走时；
根据所述第一网格四个顶点的走时、所述第二网格四个顶点的走时分别计算所述第一网格四个顶点的慢度值、所述第二网格四个顶点的慢度值；
根据所述第一网格四个顶点的慢度值、所述第二网格四个顶点的慢度值分别插值得到所述成像点的第一慢度值、第二慢度值；
根据所述第一慢度值、所述第二慢度值分别计算所述叠前地震道对应的炮点到所述成像点的走时、所述叠前地震道对应的检波点到所述成像点的走时。
12.根据权利要求10所述的积分法叠前深度偏移方法，其特征在于，所述线程获取所述叠前地震道对应的炮点到所述成像点的等效Q值，所述叠前地震道对应的检波点到所述成像点的等效Q值的步骤，包括：
根据所述炮点坐标、所述检波点坐标和等效Q值优化参数，计算所述炮点坐标对应的优化等效Q值的第三网格及所述第三网格四个顶点的第三索引值，以及所述检波点坐标对应的优化等效Q值的第四网格及所述第四网格四个顶点的第四索引值；
根据所述第三索引值、所述第四索引值，分别计算所述第三网格四个顶点、所述第四网格四个顶点的等效Q值在成像线优化等效Q值表的数据文件中的第三位置、第四位置，并根据所述第三位置、所述第四位置读取所述第三网格四个顶点的等效Q值、所述第四网格四个顶点的等效Q值；
根据所述第三网格四个顶点的等效Q值、所述第四网格四个顶点的等效Q值，分别插值得到所述叠前地震道对应的炮点到所述成像点的等效Q值、所述叠前地震道对应的检波点到所述成像点的等效Q值。
13.一种积分法叠前深度偏移方法，其特征在于，包括：
CPU获取目标工区的叠前地震数据及偏移参数；
所述CPU根据所述叠前地震数据及所述偏移参数计算所述目标工区的激发点走时表和激发点等效Q值表；
所述CPU根据所述激发点走时表和所述激发点等效Q值表，计算所述目标工区成像空间的成像空间走时表和成像空间等效Q值表；所述成像空间走时表被划分为多个成像线走时表，所述成像空间等效Q值表被划分为多个成像线等效Q值表；
所述CPU运用权利要求1所述的数据优化方法优化所述成像线走时表和所述成像线等效Q值表，得到并存储成像线优化走时表和成像线优化等效Q值表；
所述CPU以成像线为单位将所述成像线优化走时表、所述成像线优化等效Q值表发送给GPU，并将所述叠前地震数据及所述偏移参数发送给所述GPU；
所述GPU线程读取所述成像线优化走时表和所述成像线优化等效Q值表，以根据所述成像线优化走时表、所述成像线优化等效Q值表、所述叠前地震数据及所述偏移参数计算所述目标工区成像空间内各个成像点对应不同偏移距的偏移幅值。
CN201910397720.9A 2019-05-14 2019-05-14 数据优化方法及积分法叠前深度偏移方法 Active CN110031898B (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
CN201910397720.9A CN110031898B (zh) 2019-05-14 2019-05-14 数据优化方法及积分法叠前深度偏移方法
US16/556,430 US11209563B2 (en) 2019-05-14 2019-08-30 Data optimization method and integral prestack depth migration method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910397720.9A CN110031898B (zh) 2019-05-14 2019-05-14 数据优化方法及积分法叠前深度偏移方法

Publications (2)

Publication Number Publication Date
CN110031898A true CN110031898A (zh) 2019-07-19
CN110031898B CN110031898B (zh) 2019-12-13

Family

ID=67242069

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910397720.9A Active CN110031898B (zh) 2019-05-14 2019-05-14 数据优化方法及积分法叠前深度偏移方法

Country Status (2)

Country Link
US (1) US11209563B2 (zh)
CN (1) CN110031898B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112305590A (zh) * 2020-09-23 2021-02-02 中国石油天然气集团有限公司 粘声介质叠前时间偏移计算方法及装置
US11774615B2 (en) 2021-12-09 2023-10-03 Saudi Arabian Oil Company Method and systems for computational efficiency 3D prestack Kirchhoff depth migration

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US11836490B2 (en) * 2019-11-14 2023-12-05 Nvidia Corporation Kernel coefficient quantization
CN116774293B (zh) * 2023-08-25 2023-10-27 浙江大学海南研究院 一种同相轴自动拾取方法、***、电子设备及介质

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101545986A (zh) * 2009-05-06 2009-09-30 匡斌 基于最大能量旅行时计算的三维积分叠前深度偏移方法
CN103020945A (zh) * 2011-09-21 2013-04-03 中国科学院电子学研究所 一种多源传感器的遥感图像配准方法
CN106291691A (zh) * 2016-08-22 2017-01-04 中国石油天然气集团公司 一种地震偏移成像方法及装置
CN106548462A (zh) * 2016-11-02 2017-03-29 西安电子科技大学 基于薄板样条插值的非线性sar图像几何校正方法
CN108710148A (zh) * 2018-05-29 2018-10-26 中国科学院地质与地球物理研究所 三维倾角域稳相叠前深度偏移方法和装置
US20190154857A1 (en) * 2017-11-17 2019-05-23 Cgg Services Sas Seismic exploration using image-based reflection full waveform inversion to update low wavenumber velocity model

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20160013773A1 (en) * 2012-11-06 2016-01-14 Pavel Dourbal Method and apparatus for fast digital filtering and signal processing
US11112516B2 (en) * 2018-04-30 2021-09-07 Schlumberger Technology Corporation Data fusion technique to compute reservoir quality and completion quality by combining various log measurements

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101545986A (zh) * 2009-05-06 2009-09-30 匡斌 基于最大能量旅行时计算的三维积分叠前深度偏移方法
CN103020945A (zh) * 2011-09-21 2013-04-03 中国科学院电子学研究所 一种多源传感器的遥感图像配准方法
CN106291691A (zh) * 2016-08-22 2017-01-04 中国石油天然气集团公司 一种地震偏移成像方法及装置
CN106548462A (zh) * 2016-11-02 2017-03-29 西安电子科技大学 基于薄板样条插值的非线性sar图像几何校正方法
US20190154857A1 (en) * 2017-11-17 2019-05-23 Cgg Services Sas Seismic exploration using image-based reflection full waveform inversion to update low wavenumber velocity model
CN108710148A (zh) * 2018-05-29 2018-10-26 中国科学院地质与地球物理研究所 三维倾角域稳相叠前深度偏移方法和装置

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
王华忠等: "大规模三维地震数据Kirchhoff叠前深度偏移及其并行实现", 《石油地球物理勘探》 *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112305590A (zh) * 2020-09-23 2021-02-02 中国石油天然气集团有限公司 粘声介质叠前时间偏移计算方法及装置
CN112305590B (zh) * 2020-09-23 2023-12-22 中国石油天然气集团有限公司 粘声介质叠前时间偏移计算方法及装置
US11774615B2 (en) 2021-12-09 2023-10-03 Saudi Arabian Oil Company Method and systems for computational efficiency 3D prestack Kirchhoff depth migration

Also Published As

Publication number Publication date
US20200363549A1 (en) 2020-11-19
CN110031898B (zh) 2019-12-13
US11209563B2 (en) 2021-12-28

Similar Documents

Publication Publication Date Title
CN110031898A (zh) 数据优化方法及积分法叠前深度偏移方法
CN108445538B (zh) 基于反射地震资料建立深度域层q模型的方法和***
US9759826B2 (en) System and method for generating an implicit model of geological horizons
CN110031896A (zh) 基于多点地质统计学先验信息的地震随机反演方法及装置
CN109085648B (zh) 叠前深度偏移方法和装置
CN102066980A (zh) 地震层位骨架化
CN108845351A (zh) 一种vsp地震资料转换波全波形反演方法
CN107462924B (zh) 一种不依赖于测井资料的绝对波阻抗反演方法
CN106094029A (zh) 利用偏移距矢量片地震数据预测储层裂缝的方法
CN103926619B (zh) 一种三维vsp数据的逆时偏移方法
CN109061728B (zh) 一种滩坝砂体精细预测方法
AU2007329168A1 (en) Method of building a sub surface velocity model
CN102590862A (zh) 补偿吸收衰减的叠前时间偏移方法
CN109001813A (zh) 一种压制多次波的方法、装置及***
CN109188506A (zh) 一种适用于高铁隧底地震ct的纯地表三维观测***
CN105372705A (zh) 一种基于多波资料的地层切片方法
CN112508399B (zh) 一种地下空间开发适宜性评价方法及***
WO2021127382A1 (en) Full waveform inversion in the midpoint-offset domain
CN110927793A (zh) 一种基于序贯随机模糊模拟的储层预测方法及***
CN105137479A (zh) 一种面元覆盖次数的计算方法及装置
Hole et al. Interface inversion using broadside seismic refraction data and three‐dimensional travel time calculations
CN116859478B (zh) 一种基于瞬变电磁法成像的地下水模拟方法及***
CN103116183B (zh) 一种石油地震采集面元覆盖次数属性体切片成图方法
CN108363097A (zh) 一种地震资料偏移成像方法
CN103543478A (zh) 地质形态插值的km方法

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
GR01 Patent grant
GR01 Patent grant