数据优化方法及积分法叠前深度偏移方法与流程
本发明涉及地震勘探技术领域，尤其是涉及一种数据优化方法及积分法叠前深度偏移方法。
背景技术：
现行地震叠前偏移方法的实践中，以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值表的数据文件；第一关系式为：
第二关系式为：
第三方面，本发明实施例还提供了一种积分法叠前深度偏移方法，应用于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]中抽取该炮点至该成像点的走时
成像空间内所有成像点的走时表集合组成目标工区的成像空间走时表τ[ny，nx，nz，nysp，nxsp]。成像空间内所有成像点的等效q值表集合组成目标工区的成像空间等效q值表q[ny，nx，nz，nysp，nxsp]。
对目标工区的成像空间ispace[ny，nx，nz]内的任一成像点(yp，xp，zp)，该成像点的走时表
(1)将该成像点的走时表的元素
式中，
(2)根据预设的网格密度
令目标函数1
最小。式中aki为由该成像点坐标、炮点坐标、以及所应用的插值方法确定的权系数。可得到关于
式中，
这里，上述令目标函数1最小也即：对
求解该方程组即可得到mxsp×mysp个炮点网格上的慢度值组成第二慢度序列
(3)基于第二慢度序列应用插值算法得到原始nysp×nxsp个炮点网格上的第三慢度序列
(4)基于步骤(3)得到的原始nysp×nxsp个炮点网格上的第三慢度序列
由该第三慢度序列组成新的慢度矩阵
(5)求取上述新的走时矩阵
当差的绝对值大于1毫秒的炮(元素)个数大于有效成像孔径范围总炮数(矩阵元素总个数)nysp×nxsp的1％时，按预设规则增加网格密度mxsp、mysp，这里，在其中一种实施方式中，对于上述预设的网格密度mxsp×mysp，其中mxsp、mysp均为奇数，且每次增加网格密度时，x、y方向各增加两个点。
增加该网格密度后，重新回到上述第(2)步，循环计算直至上述差的绝对值大于1毫秒的炮个数小于有效成像孔径范围内总炮数nysp×nxsp的1％，此时，记录达到条件的网格密度对应的采样个数及间距为走时表共享优化参数，具体地，记录采样个数mxsp为
相似地，对于上述成像点(yp，xp，zp)的等效q值表的优化过程如下。
其中，该成像点的等效q值表
(6)将该成像点的等效q值表的所有元素组成一个包含nysp×nxsp个数的第一等效q值序列qi，i＝1，2，3，…nysp×nxsp。
(7)根据预设的网格密度
令目标函数2：
最小。式中aki为由该成像点坐标、炮点坐标、以及所应用的插值方法确定的权系数。可得到关于
式中，
这里，令上述目标函数2最小也即：对
求解该方程组即可得到坐标为(yp，xp，zp)的成像点处mxsp×mysp个炮点网格上的等效q值组成的第二等效q值序列
(8)基于该第二等效q值序列应用插值算法得到原始nysp×nxsp个炮点网格上的第三等效q值序列
(9)基于步骤(8)得到的原始nysp×nxsp个炮点网格上的等效q值序列
(10)求取该新的等效q值矩阵
当差的绝对值大于5的炮(元素)个数，大于有效成像孔径范围总炮数(矩阵元素总个数)nysp×nxsp的1％时，按预设规则增加mxsp×mysp中的mxsp、mysp的值，并重新回到第(7)步，循环计算直至上述差的绝对值大于5的炮个数小于有效成像孔径范围内总炮数nysp×nxsp的1％，记录达到条件的网格密度对应的采样个数及间距为等效q值表共享优化参数，具体地，记录采样个数mxsp为
这样，通过上述步骤即可得到目标工区的成像空间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值表优化参数
遍历该典型成像点所在的成像块内的所有成像点，得到该成像块的优化等效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方向索引，
这里，在对数据文件进行存储时，针对某一成像块，存在一个数组顺序重排的操作，是为了提高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值表数据文件按照上述实施例三中的第二关系式直接存储地震波自优化后得到的
另外，在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方向坐标确定该成像点所在的成像块、走时表优化参数
步骤48：gpu每个线程根据该叠前地震道炮点坐标sx，sy，检波点坐标gx，gy以及该成像线对应的成像线的优化走时表和优化等效q值表存储文件，分别获取该叠前地震道对应的炮点到该成像点的走时τs、等效q值qs，该叠前地震道对应的检波点到该成像点到的走时τg、等效q值qg。
具体地，计算叠前地震道对应的炮点到该成像点的走时τs、检波点到该成像点的走时τg的步骤包括：
步骤481：由该叠前地震道炮点坐标sx，sy及优化走时表参数
(ixs1，iys1)，(ixs2，iys2)，(ixs3，iys3)，(ixs4，iys4)
对于(ixsj，iysj)，j＝1，2，3，4计算其走时在成像线的优化走时表数据存储文件中的位置：
并读取相应走时
步骤482：利用公式
同理，将上述步骤481～482中的sx更换为gx，sy更换为gy，即可得到该叠前地震道对应的检波点到该成像点的走时τg。
另外，计算该叠前地震道对应的炮点到该成像点的等效q值qs、检波点到该成像点到等效q值qg的步骤包括：
步骤483：由该叠前地震道炮点坐标sx，sy及优化等效q值表参数
对于
并读取相应等效q值
步骤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-onlymemory)、随机存取存储器(ram，randomaccessmemory)、磁碟或者光盘等各种可以存储程序代码的介质。
最后应说明的是：以上所述实施例，仅为本发明的具体实施方式，用以说明本发明的技术方案，而非对其限制，本发明的保护范围并不局限于此，尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，其依然可以对前述实施例所记载的技术方案进行修改或可轻易想到变化，或者对其中部分技术特征进行等同替换；而这些修改、变化或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。
- 还没有人留言评论。精彩留言会获得点赞！