CN103513279B - 一种基于地震波波动方程的照明分析计算方法及计算装置 - Google Patents

Abstract

本发明提供一种基于地震波波动方程的照明分析计算方法包括：建立地质模型，读入观测系统参数；将地质模型离散化成网格模型；读取第一层的数据准备进行计算；设定一个初始波场；计算下一层相应位置处的波场；根据计算得到的下一层各点的波场值，计算下一层各点的照明能量；计算装置输出能量分布图。本发明提供的一种基于地震波波动方程的照明分析计算方法，炮点和检波点位置相同的点只计算一次，减少了计算量；计算时，将波场值存储在计算机内存中，降低了对计算机存储空间的要求并提高了计算的速度。本发明还提供了一种基于地震波波动方程的照明分析计算装置实施例，该实施例可以实现本发明基于地震波波动方程的照明分析计算。

Description

一种基于地震波波动方程的照明分析计算方法及计算装置

技术领域

本发明涉及地球物理勘探技术领域，特别涉及一种基于地震波波动方程的照明分析计算方法及计算装置。

背景技术

反射地震方法是获取地下结构的重要手段之一。反射地震方法的基本原理是：由震源发出地震波传播到地下，遇到目标后反射回地表，携带着地下结构信息的反射波被设置在地表的接收系统记录到后，再经相应的数据处理方法对地下结构进行成像，从而得到关于地下结构的描述。反射地震资料的现场采集和后期处理需要非常高的成本，如果能在进行这些工作之前对观测系统的潜在探测能力提供一些评估，将大大降低现场采集的成本。基于地震波波动方程的照明分析计算方法，正是在给定观测系统和地下背景结构的情况下对反射地震探测能力提供定量描述的有效方法。基于地震波波动方程的照明分析方法已经在地球物理勘探技术领域中被广泛应用在地震资料采集设计、处理和解释的各个阶段。例如，在进行实际数据采集前，以事先针对工区的特点利用照明分析对不同的采集方案进行对比和优化，从而有可能在提高采集质量和效率的同时减少对昂贵采集设备和野外人工的投入。

基于地震波波动方程的照明分析计算方法，具体的实施步骤主要包括：先根据地质构造建立地质模型并设定观测系统，在根据设定的观测系统将地质模型转换为网格模型；逐层计算一个炮点或检波点在各层的波场值，并将计算得到的波场值存储在计算机硬盘上；再用上述方法重复计算其他炮点或检波点在各层的波场值并存储在计算机硬盘上，直至计算完所有的炮点或检波点在各层的波场值；最后从计算机硬盘上读取存储的波场值，计算地质模型上各点的照明能量。

上述基于地震波波动方程的照明分析计算方法，需要将计算的波场值全部存储在计算机硬盘上，要求计算机有大容量的存储空间；同时，计算机从硬盘上读取数据需要的时间较多，影响照明分析计算的速度。

发明内容

本发明的目的在于提供一种基于地震波波动方程的照明分析计算方法及计算装置，以实现降低对计算机存储空间的要求并提高照明分析的计算速度。

一种基于地震波波动方程的照明分析计算方法，包括：

S1：计算装置根据工区地质构造建立地质模型，并读入地震观测系统的参数；

S2：计算装置根据地震观测系统，将建立的地质模型离散化成多层的网格模型；

S3：读取上述网格模型的第一层的数据至内存，所述第一层作为当前层，准备进行计算；

S4：计算装置以炮点或检波点位置为中心，在计算区域内设定一个初始波场；

S5：计算装置根据初始波场和传播算子，计算当前层的下一层相应位置处的波场；

S6：计算装置判断当前作为中心的炮点或检波点是否是观测系统中当前层的最后一个炮点或检波点，并根据判断的结果执行相应的操作；

S7：计算装置根据计算得到的所述当前层的下一层各点的波场值，计算所述当前层的下一层各点的照明能量；

S8：计算装置判断计算出照明能量的层是否是网格模型的最后一层，并根据判断的结果执行相应的操作；

S9：计算装置输出每一层的能量分布图。

其中：

所述S1中观测系统的参数，具体包括：地震波频率、炮点位置和检波点位置；

所述S6，具体包括：判断当前作为中心的炮点或检波点是否是观测系统中当前层的最后一个炮点或检波点，若是最后一个炮点或检波点，则丢弃内存中原来存储的当前层的波场值，将计算所得的当前层的下一层的波场值存储在内存中；若当前作为中心的作用点不是观测系统中当前层的最后一个炮点或检波点，则将上述S4中的中心向右偏移至下一个炮点或检波点，并重复S4-S6进行计算和判断，直至该层所有观测系统中的炮点和检波点都计算结束，丢弃内存中原来存储的当前层的波场值，将计算所得的当前层的下一层的波场值存放在内存中；

所述S8，具体包括：先判断S7计算出的照明能量是否为网格模型的最后一层的照明能量，如是，则计算结束；若计算出照明能量的层不是网格模型的最后一层，则将该计算出照明能量的层作为当前层，并重复S4-S8进行计算和判断，直至计算出网格模型中最后一层的照明能量，则计算结束。

所述的一种基于地震波波动方程的照明分析计算方法，S4中的计算区域是按照如下方法确定的：先将所有炮点和检波点都看作是作用点，检波点和炮点重复的位置算作一个作用点，所述计算区域包括以作用点为中心，从中心的左边Maxoffset/DX个网格点开始，至中心的右边Maxoffset/DX个网格点结束的区间。

所述的一种基于地震波波动方程的照明分析计算方法，S4中在计算区域内设定一个初始波场，包括：如果当前层为网格模型的第一层，则在波数域内将计算区域的中心位置处实部设为大于0的任意值，计算区域的其他位置处实部设为0，整个计算区域的虚部都取0；如果当前层不是第一层，则在波数域内将该区域根据上一层相同位置的计算区域计算所得的波场值设定为初始波场值。

所述的一种基于地震波波动方程的照明分析计算方法，S7中计算装置根据计算得到的下一层各点的波场值，计算下一层各点的照明能量，包括：先将S5计算得到的下一层的波数域波场值作傅立叶变换至空间域；然后在空间域根据观测系统中的炮检关系，对于每一个作用点，先计算出该点作为中心点时在下一层的相同位置的波场值，然后计算出与其相关的接收点作为中心点时，该点在各区域计算得到的下一层相同位置的波场值，将该点作为中心点时得到的下一层的波场值进行平方后分别和与其相关的接收点作为中心点时该点的下一层波场值的平方相乘并求和，计算所得的结果记作该计算区域中心点位置的照明能量DI，并存储在计算机的硬盘上。

所述的一种基于地震波波动方程的照明分析计算方法，S5中，传播算子取为exp(i×DX×kz)，其中，i为虚数单位，DX是网格模型的单元边长，为已知量，kz是垂直波数；所述垂直波数可以通过下式计算得到：其中，kx为水平波数；所述水平波数kx可以通过下式计算：kx＝x×2×3.1415926×f，其中，f是地震波的频率，x是当前点的位置，f和x都为已知量。

所述的一种基于地震波波动方程的照明分析计算方法，S2中将建立的地质模型离散化成网格模型，具体包括：先计算检波点间距的二分之一和炮点间距的二分之一，即RI/2和SI/2，其中，RI是检波点间距，SI是炮点间距；然后取DX作为离散后的网格模型的单元边长，DX的取值方法为：取RI/2和SI/2中较小的数作为DX的值；再以DX为单元将S1中建立的地质模型离散化成网格模型。

本发明还提供一种基于地震波波动方程的照明分析计算装置，包括：计算前处理模块、模型网格化模块、读取数据模块、波场值设定模块、波场值计算模块、第一判断模块、照明能量计算模块、第二判断模块和能量分布输出模块。其中：

计算前处理模块，用于根据工区地质构造建立地质模型，并读入观测系统的参数；

模型网格化模块，用于将计算前处理模块建立的地质模型离散化成网格模型；

读取数据模块，用于将网格模型的第一层的网格点数据读入内存中，准备进行计算；

波场值设定模块，用于将炮点或检波点位置作为中心，在计算区域内设定一个波数域的初始波场；

波场值计算模块，用于根据波场值设定模块中设定的初始波场，利用传播算子，计算下一层相应位置处的波场值；

第一判断模块，用于判断当前作为中心的炮点或检波点是否是观测系统中该层的最后一个炮点或检波点。若判断的结果为：当前作为中心的炮点或检波点是观测系统中该层的最后一个炮点或检波点，则丢弃内存中原来存储的那一层的波场值，将计算所得的下一层的波场值存储在内存中；若判断结果为：当前作为中心的炮点或检波点不是观测系统中该层的最后一个炮点或检波点，则将上述波场设定模块中的中心向右偏移至下一个炮点或检波点，并利用波场计算模块和第一判断模块进行计算和判断，直至计算完最后一个炮点或检波点，再丢弃内存中原来存储的波场值，将计算所得的下一层的波场值存放在内存中；

照明能量计算模块，用于根据波场值计算模块计算得到的下一层的各点所有波场值，计算下一层的照明能量；

第二判断模块，用于判断计算出照明能量的层是否为网格模型的最后一层；若判断结果为：照明能量计算模块中计算出照明能量的层的是网格模型的最后一层，则计算结束；若判断结果为：计算出照明能量的层不是网格模型的最后一层，则将该计算出照明能量的层作为当前层，并利用波场值设定模块、波场值计算模块、第一判断模块、照明能量计算模块和第二判断模块再次计算下一层的照明能量并进行判断，直至计算出网格模型中最后一层的照明能量，则计算结束；

能量分布输出模块，用于输出每一层的能量分布图。

所述的一种基于地震波波动方程的照明分析计算装置，波场值设定模块，用于设定波场的初始波场值；对于网格模型中第一层的数据，波场值设定模块采用的波场设定方法为：在波数域将作为计算区域中心的作用点实部设为大于0的任意值，其他网格点位置的实部设为0，整个计算区域的虚部都设为0；对于网格模型中其它层的数据，波场值设定模块采用的波场设定方法为：则将该计算区域在波数域根据上一层相同位置的计算区域计算所得的波场值设定为初始波场值。

所述的一种基于地震波波动方程的照明分析计算装置，照明能量计算模块，用于根据计算得到的下一层的所有波场值，计算下一层的照明能量；具体地，用于先将计算得到的下一层的波数域波场值作傅立叶变换至空间域；然后根据观测系统中的炮检关系，对于每一个炮点或检波点点，先计算出该点作为中心点时在下一层的相同位置的波场值，然后计算出与其相关的接收点作为中心点时，该点在各区域计算得到的下一层相同位置的波场值，将该点作为中心点时得到的下一层的波场值进行平方后分别和与其相关的接收点作为中心点时该点的下一层波场值的平方相乘并求和，计算所得的结果记作该计算区域中心点位置的照明能量DI，并存储在计算机的硬盘上。

所述的一种基于地震波波动方程的照明分析计算装置，波场值计算模块，用于在波数域将设定的初始波场乘以传播算子得到网格模型中下一层相应位置处的波场值；所述传播算子可以取为exp(i×DX×kz)，其中i为虚数单位，DX是网格模型的单元边长，为已知量，kz是垂直波数；垂直波数可以通过下式计算得到：其中，kx为水平波数；水平波数kx可以通过下式计算：kx＝x×2×3.1415926×f，其中，f是地震波的频率，x是当前点的位置，f和x都为已知量。

本发明提供的一种基于地震波波动方程的照明分析计算方法及计算装置，逐层计算所有炮点和检波点的波场值，并根据波场值计算各点的照明能量。计算时，将波场值存储在计算机内存中，在计算出了新的一层的波场值后，立即抛弃原来存储在内存中的波场值，可以大大降低对计算机存储空间的要求；计算过程中，计算的波场值不需要存入计算机硬盘中，在计算下一层波场值时，直接从内存中读取数据，读取速度快，计算的速度得到了提高；同时本发明提供的方法在计算过程中，炮点和检波点位置相同的点只计算一次，减少了计算量，可以提高计算速度。

附图说明

图1是本发明基于地震波波动方程的照明分析计算方法流程图；

图2是本发明基于地震波波动方程的照明分析计算装置实施例；

图3是本发明实施例实际工区地质模型；

图4是本发明对图3的地质模型进行离散后的网格模型；

图5是本发明对图4进行照明分析计算的示意图；

图6是本发明对图3进行照明分析计算所得的照明能量分布图。

具体实施方式

为了使本技术领域的人员更好地理解本申请中的技术方案，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都应当属于本发明保护的范围。

图1所示为本发明基于地震波波动方程的照明分析计算方法流程图。由图1可以看出，本发明基于地震波波动方程的照明分析计算方法，其步骤包括：

S1：计算装置根据工区地质构造建立地质模型，并读入地震观测系统的参数；

S2：计算装置根据观测系统，将建立的地质模型离散化成多层的网格模型；

S3：读取上述网格模型的第一层的数据至内存中，所述第一层作为当前层，准备进行计算；

S4：计算装置以炮点或检波点位置为中心，在计算区域内设定一个初始波场；

S5：计算装置根据初始波场和传播算子，计算当前层的下一层相应位置处的波场；

S6：计算装置判断当前作为中心的炮点或检波点是否是观测系统中当前层的最后一个炮点或检波点，并根据判断的结果执行相应的操作；

S7：计算装置根据计算得到的所述当前层的下一层各点的波场值，计算所述当前层的下一层各点的照明能量；

S8：计算装置判断计算出照明能量的层是否是网格模型的最后一层，并根据判断的结果执行相应的操作；

S9：计算装置输出每一层的能量分布图。其中，

所述S1根据工区地质构造建立地质模型是指，根据工区地表的地震波传播速度，地下的断层，地下的地层等已知信息，建立一个和实际工地一一对应的数学模型。所述工区，就是进行地震勘探的位置。地震勘探的目的就是找到特殊的地质构造，这些构造里能够蕴含油气，经过一次勘探的地区，会有一些构造是已知的，通过这些已知的构造，我们可以建立一个数学模型，这个数学模型就是工地实际构造的一个数学表达。所述数学模型的长度为x，深度为z，x和z都为已知量；所述S1中读入观测系统的参数，是指读入设定好的观测系统的地震波频率、炮点位置和检波点位置，炮点位置和检波点位置为人为设定，都为已知参量。通常情况下，设定的观测系统的长度比构造的地质模型的长度小很多。根据地质模型、炮点位置和检波点位置，可以得到炮点间距SI、检波点间距RI、炮点个数n、检波点个数k、最大偏移距Maxoffset。其中，炮点间距表示两个炮点之间的距离；检波点间距是指两个检波点之间的距离；炮点可以位于其中两个相邻检波点的中间位置；最大偏移距是指最远的检波点和炮点之间的距离。

所述S2中将地质模型离散化成网格模型，具体包括：先计算检波点间距的二分之一和炮点间距的二分之一，即RI/2和SI/2；然后取DX作为离散后的网格模型的单元边长，其中DX的取值方法为DX=min(RI/2,SI/2)，即取RI/2和SI/2中较小的数作为DX的值；再以DX为单元将S1中建立的地质模型离散化成网格模型，使得地质模型在水平方向上被分为x/DX列，在深度方向上被分成z/DX层。

所述S3中，取上述网格模型的第一层的数据是指将S2中离散后的网格模型中深度方向上的第一层作为开始计算的当前层。将当前层的所有网格点的数据读入到内存中，准备计算。

所述S4中，先将所有炮点和检波点都看作是作用点，检波点和炮点重复的位置算作一个作用点。计算区域是指以作用点为中心，从中心的左边Maxoffset/DX个网格点开始，至中心的右边Maxoffset/DX个网格点结束的区间。所述设定一个波场值包括：如果当前层为网格模型的第一层，则在波数域将计算区域内的中心位置处实部设为1，其他位置的网格点处实部设为0，整个计算区域的网格点的虚部都设为0；如果当前层不是第一层，则将该区域根据上一层相同位置的计算区域计算所得的波场值设定为初始波场值。

所述S5中根据初始波场和传播算子，计算下一层相应位置处的波场，具体包括：对于每一个炮点或检波点，将S4中设定好的初始波场值乘以传播算子，所得到的值即为该网格点在下一层对应位置处的波场值。所述传播算子可以取为exp(i×DX×kz)，其中i是虚数单位,kz是垂直波数。波数是波传播的方向上单位长度内的波周数目，水平波数是指波水平方向上传播的单位长度内的波周数目，垂直波数是指波垂直方向上传播的单位长度内的波周数目。垂直波数可以通过下式计算得到：

其中，kx为水平波数，kx可以通过计算得到，计算公式如下：

kx＝x×2×3.1415926×f

其中，f是地震波的频率，x是当前点的位置。

所述S6，具体包括：当前作为中心的作用点是否是观测系统中当前层的最后一个作用点，若是最后一个作用点，则丢弃内存中原来存储的当前层的波场值，将计算所得的下一层的波场值存储在内存中；若当前作为中心的作用点不是观测系统中当前层的最后一个作用点，则将上述S4中的中心向右偏移至下一个作用点，并重复S4-S6进行计算和判断，直至观测系统中该层所有的作用点都计算结束，再丢弃内存中原来存储的当前层的波场值，并将计算所得的当前层的下一层的波场值存放在内存中。

所述S7中根据计算得到的所述当前层的下一层的波场值，计算所述当前层的下一层的照明能量，具体包括：先将S6所得到的下一层的波数域波场值作傅立叶变换变换到空间域；然后根据观测系统中的炮检关系，对于每一个作用点，先计算出该点作为中心点时在下一层的相同位置的波场值，然后计算出与其相关的接收点作为中心点时，该点在各区域计算得到的下一层相同位置的波场值，将该点作为中心点时得到的下一层的波场值进行平方后分别和与其相关的接收点作为中心点时该点的下一层波场值的平方相乘并求和，计算所得的结果记作该计算区域中心点位置的照明能量DI，并存储在计算机的硬盘上。

所述S8，具体包括：先判断S7计算出的照明能量是否为网格模型的最后一层的照明能量，若是，则计算结束；若计算出照明能量的层不是网格模型的最后一层，则将该计算出照明能量的层作为当前层，并重复S4-S8进行计算和判断，直至计算出网格模型中最后一层的照明能量，则计算结束。

根据本发明基于地震波波动方程的照明分析计算方法，提供一个基于地震波波动方程的照明分析计算的装置实施例。图2所示为基于地震波波动方程的照明分析计算装置实施例，如图2所示，该照明分析计算装置实施例，依次包括：计算前处理模块、模型网格化模块、读取数据模块、波场值设定模块、波场值计算模块、第一判断模块、照明能量计算模块、第二判断模块和能量分布输出模块；其中，

计算前处理模块，用于根据工区地质构造建立地质模型，并读入地震观测系统的参数；

模型网格化模块，用于根据地震观测系统，将计算前处理模块建立的地质模型离散化成网格模型；

读取数据模块，用于将网格模型的第一层的网格点数据读入内存中，所述第一层作为当前层，准备进行计算；

波场值设定模块，用于将炮点或检波点位置作为中心，在计算区域内设定一个波数域的初始波场；

波场值计算模块，用于根据波场值设定模块中设定的初始波场，利用传播算子，计算所述当前层的下一层相应位置处的波场值；

第一判断模块，用于判断当前作为中心的炮点或检波点是否是所述当前层的最后一个炮点或检波点，并根据判断的结果执行相应的操作；

照明能量计算模块，用于根据计算得到的所述当前层的下一层的所有波场值，计算所述当前层的下一层的照明能量；

第二判断模块，用于判断计算出照明能量的层是否为网格模型的最后一层，并根据判断的结果执行相应的操作；

能量分布输出模块，用于输出每一层的能量分布图。

所述波场值设定模块，用于设定波场的初始波场值。对于网格模型中第一层的数据，波场值设定模块采用的波场设定方法为：在波数域，将作为计算区域中心的作用点位置的实部设为1，其他网格点位置的实部设为0，整个计算区域的虚部都设为0；对于网格模型中其它层的数据，波场值设定模块采用的波场设定方法为：则将该区域根据上一层相同位置的计算区域计算所得的波场值设定为初始波场值。

所述波场值计算模块中的传播算子可以取为exp(i×DX×kz)，其中，i为虚数单位，DX是网格模型的单元边长，为已知参量，kz是垂直波数；

垂直波数可以通过下式计算得到：其中，kx为水平波数；

水平波数kx可以通过下式计算：kx＝x×2×3.1415926×f，其中，f是地震波的频率，x是当前点的位置，f和x都为已知量。

所述第一判断模块，用于判断当前作为中心的炮点或检波点是否是所述当前层最后一个炮点或检波点。若判断的结果为当前作为中心的炮点或检波点是所述当前层的最后一个炮点或检波点，则丢弃内存中原来存储的那一层的波场值，将计算所得的所述当前层的下一层的波场值存储在内存中；若判断结果为当前作为中心的炮点或检波点不是所述当前层的最后一个炮点或检波点，则将上述波场设定模块中的中心向右偏移至下一个炮点或检波点，并利用波场计算模块和第一判断模块进行计算和判断，直至计算完所述当前层的最后一个炮点或检波点，再丢弃内存中原来存储的波场值，将计算所得的所述当前层的下一层的波场值存放在内存中。

照明能量计算模块，用于根据计算得到的所述当前层的下一层的所有波场值，计算所述当前层的下一层的照明能量。具体地，用于先将计算得到的所述当前层的下一层的波数域波场值作傅立叶变换变换到空间域；然后根据观测系统中的炮检关系，对于每一个作用点，先计算出该点作为中心点时在所述当前层的下一层的相同位置的波场值，然后计算出与其相关的接收点作为中心点时，该点在各区域计算得到的所述当前层的下一层相同位置的波场值，将该点作为中心点时得到的所述当前层的下一层的波场值进行平方后分别和与其相关的接收点作为中心点时该点的所述当前层的下一层波场值的平方相乘并求和，计算所得的结果记作该计算区域中心点位置的照明能量DI，并存储在计算机的硬盘上。

所述第二判断模块，用于判断计算出照明能量的层是否为网格模型的最后一层。若判断结果为照明能量计算模块中计算出照明能量的层的是网格模型的最后一层，则计算结束；若判断结果为计算出照明能量的层不是网格模型的最后一层，则将该计算出照明能量的层作为当前层，并利用波场值设定模块、波场值计算模块、第一判断模块、照明能量计算模块和第二判断模块再次进行计算和判断，直至计算出网格模型中最后一层的照明能量，则计算结束。

图3所示为本发明实施例实际工区地质模型。图3中，横坐标表示长度，单位是米；纵坐标表示深度，单位是米；灰度表示地震波传播的速度值。

利用图2所示的照明分析计算装置实施例按照图1所示的照明分析计算流程对实际工区地质模型进行照明分析计算。具体步骤包括：

S1：计算装置根据工区地质构造建立地质模型，并读入地震观测系统的参数。

计算装置的计算前处理模块根据工区地质建立地质模型，并读入观测系统的参数。图3所示为根据实际工区建立的地质模型，如图3所示，该地质模型的长度和深度为已知量，例如长度取20000米，深度取6000米。设定好观测系统，包括：地震波的频率，通常取值范围为15-30HZ，本实施例中例如取15HZ；观测系统中的首炮位置，例如取5000米处；炮间距SI，例如取20米；检波点的位置，例如取2000米至18000米处；检波点间距RI，例如取20米；炮点位于两个相邻检波点的中间。根据上述设定的观测系统，可以得到：每个炮点对应的检波点的个数，例如300个；炮点个数，例如251个；炮点和检波点总数，例如801个；最大偏移距，例如Maxoffset=300*20/2=3000米。

S2：计算装置根据地震观测系统，将建立的地质模型离散化成网格模型。

计算装置中的模型网格化模块将建立的地质模型离散化成网格模型。具体地，先取检波点间距的二分之一和炮点间距的二分之一，例如取RI/2=10米和SI/2=20米；则网格步长DX取RI/2和SI/2中较小的值，例如DX=min(RI/2,SI/2)=10米；然后以10米为单元长度将地质模型离散化成网格模型。图4所示为图3的地质模型进行离散化后的网格模型，如图4所示，在水平方向上地质模型被离散化为2000列，在深度方向上地质模型被分成600层。

S3：计算装置读取网格模型的第一层的数据，将所述第一层作为当前层，准备进行计算。

计算装置中的读取数据模块将网格模型中第一层的2000个网格点的数据读入内存中，准备进行分析计算。

S4：计算装置以炮点或检波点位置为中心，在计算区域内设定一个初始波场。

计算装置中的波场值设定模块以炮点或检波点位置为中心，在计算区域内设定一个初始波场。具体地，本实施例中，以第一层网格的第一个计算区域的确定为例，确定计算区域的方法是，在第一层的网格中，以观测系统中的第一个炮点所在位置（即横坐标5000米处）为中心，左边300个网格开始，左边300个网格结束，将这601个点作为计算区域。在上述计算区域内，给定初始波场值；对于第一层的数据，在波数域设定计算区域中心位置的实部为1，其他位置的实部为0，即第301个点处的实部是1，其他位置的实部是0，并设定整个计算区域内的所有网格点的虚部为0；

如果需要设定初始波场值的计算区域所在的层不是网格模型的第一层，则将该区域根据上一层相同位置的计算区域计算所得的波场值设定为初始波场值。

为方便说明，引入图5，图5所示是对图4进行照明分析计算的示意图。图5中的每个计算区域只取了5个炮点或检波点。设定初始波场时，以第一层为例，设定图5中第一计算区域的初始波场，将红色的作为计算区域的中心，称为作用点，此时不考虑其他的炮点或检波点，将计算区域中心点的波数域初始波场值的实部设为1，其余的网格点的波数域初始波场值的实部设为0，整个计算区域的网格点的虚部则都设为0。

S5：计算装置根据初始波场和传播算子，计算所述当前层的下一层相应位置处的波场。

计算装置中的波场值计算模块根据初始波场和传播算子，计算下一层相应位置处的波场。具体地，将S4设定的计算区域内每个网格点的初始波场乘以传播算子，所得到的值即为该炮点或检波点在下一层对应位置处的波场值。所述传播算子可以取为exp(i×DX×kz)，其中i是虚数单位,DX是网格模型的单元长度，为已知量，kz是垂直波数；

垂直波数可以通过下式计算得到：其中，kx为水平波数；

水平波数kx可以通过下式计算得到：kx＝x×2×3.1415926×f，其中，f是输入的频率，x是当前点的位置；f和x都为已知量。

S6：计算装置判断当前作为中心的炮点或检波点是否是观测系统中所述当前层的最后一个炮点或检波点，并根据判断的结果执行相应的操作；

具体地，计算装置中的第一判断模块先判断当前作为中心的作用点是否是观测系统的最后一个作用点，若是最后一个作用点，则丢弃内存中原来存储的当前层的波场值，将计算所得的下一层的波场值存储在内存中；若当前作为中心的作用点不是观测系统的最后一个作用点，则将上述S4中的中心向右偏移至下一个作用点，并重复S4-S6进行计算和判断，直至该层所有观测系统中的作用点都计算结束，再将计算所得的下一层的波场值存放在内存中。根据上述方法，对于图4所示的网格模型和观测系统，可以计算出所有包含炮点或检点位置的波场值，如果炮点和检波点在同一位置，则该位置只计算一次。这里由于炮点和检波点的总数为801个，因此共计算801次，得到801个位置的第二层波场，每个位置的波场可包含301个值。

S7：计算装置根据计算得到的所述当前层的下一层的波场值，计算所述当前层的下一层的照明能量。

计算装置中的照明能量计算模块根据计算得到的下一层的波场值，计算下一层的照明能量。具体地，本实施例中对于图5中的第五个点，在第一计算区域中，将所有点的波场值乘以传播算子，分别得到第二层相应位置处的波场值，将此时计算得到的第五点的波场值记为a₅₁；将中心向右偏移至下一个作用点，同理可以得到a₅₂、a₅₃、a₅₄、a₅₅，其中，a₅₃为第五个点在第三计算区域中作为中心点的计算结果，而其他得到的波场值被称为与其相关的接收点在计算下一层相同坐标位置的波场值时得到的波场值；然后将得到的波场值进行平方，将其他波场值的平方与作为中心时得到的波场值的平方逐一相乘并求和，所得结果即为第五点的照明能量值，即 ${\mathrm{DI}}_5=a_{51}^2\×a_{53}^2+a_{52}^2\×a_{53}^2+a_{54}^2\×a_{53}^2+a_{55}^2\×a_{53}^2.$

根据上述方法，将作为中心点计算的波场值的平方分别和与其相关的接收点在计算下一层相同坐标位置的波场值时得到的300个波场值平方相乘并求和，得到的值即为所述801个作用点在第二层位置处的能量值。计算结束后，将上述求得的能量值存储在计算机硬盘上。

S8：计算装置判断计算出照明能量的层是否是网格模型的最后一层，并根据判断的结果执行相应的操作。

计算装置中的第二判断模块先判断S7计算出的照明能量是否为网格模型的最后一层的照明能量，若判断的结果为计算出照明能量的层是网格模型的最后一层，则计算结束；若计算出照明能量的层不是网格模型的最后一层，则将该计算出照明能量的层作为当前层，并重复S3-S8进行计算和判断，直至计算出网格模型中最后一层的照明能量，计算结束。

S9：计算装置输出每一层的能量分布图。

计算装置中的能量分布输出模块输出计算出的能量分布图。图6是本发明对图3分析计算所得的照明能量分布图。如图6所示，横坐标表示水平网格数，每网格的单位是10米，纵坐标表示深度网格数，每网格单位是10米；图中灰度表示照明能量值；图中，照明能量分布图的长度为20000米，分为2000个网格；深度为6000米，分为600个网格；由于观测系统设定在5000米至18000米处，因此，只有该范围有能量值。

本发明提供的一种基于地震波波动方程的照明分析计算方法，逐层计算所有炮点和检波点的波场值，并根据波场值计算各点的照明能量，将计算出的波场值存储在计算机硬盘上，最后输出能量分布图。计算过程中，炮点和检波点位置相同的只计算一次，减少了计算量；计算时将波场值存储在计算机内存中，在计算出了新的一层的波场值后，立即抛弃原来存储在内存中的波场值，这样可以大大降低对计算机存储空间的要求；同时，计算的波场值不需要存入计算机硬盘中，读取时直接从内存中读取，读取速度快，计算的速度得到了提高。

需要说明的是，本发明方法实施例和装置实施例中的传播算子可以选取其他算子，只需要满足波动方程即可，本发明对此并不作出限定。还需要说明的是，本发明方法实施例和装置实施例中对第一层设定初始波场时，中心位置的作用点的取值可以替换为大于0的任意值，只要能表示作用点处有一个激发信号即可，本发明对此并不作出限定。

虽然通过实施例描绘了本发明，本领域普通技术人员知道，本发明有许多变形和变化而不脱离本发明的精神，希望所附的权利要求包括这些变形和变化而不脱离本发明的精神。

Claims (8)

1.一种基于地震波波动方程的照明分析计算方法，其特征在于，包括：

S1：计算装置根据工区地质构造建立地质模型，并读入地震观测系统的参数；

S2：计算装置根据地震观测系统，将建立的地质模型离散化成多层的网格模型；

S3：读取上述网格模型的第一层的数据至内存，所述第一层作为当前层，准备进行计算；

S4：计算装置以炮点或检波点位置为中心，在计算区域内设定一个初始波场；

S5：计算装置根据初始波场和传播算子，计算当前层的下一层相应位置处的波场；

S6：计算装置判断当前作为中心的炮点或检波点是否是观测系统中当前层的最后一个炮点或检波点，并根据判断的结果执行相应的操作；

S7：计算装置根据计算得到的所述当前层的下一层各点的波场值，计算所述当前层的下一层各点的照明能量；

S8：计算装置判断计算出照明能量的层是否是网格模型的最后一层，并根据判断的结果执行相应的操作；

S9：计算装置输出每一层的能量分布图；

其中，

所述S1中观测系统的参数，具体包括：地震波频率、炮点位置和检波点位置；

所述S6，具体包括：判断当前作为中心的炮点或检波点是否是观测系统中当前层的最后一个炮点或检波点，若是最后一个炮点或检波点，则丢弃内存中原来存储的当前层的波场值，将计算所得的当前层的下一层的波场值存储在内存中；若当前作为中心的作用点不是观测系统中当前层的最后一个炮点或检波点，则将上述S4中的中心向右偏移至下一个炮点或检波点，并重复S4-S6进行计算和判断，直至该层所有观测系统中的炮点和检波点都计算结束，丢弃内存中原来存储的当前层的波场值，将计算所得的当前层的下一层的波场值存放在内存中；

所述S8，具体包括：先判断S7计算出的照明能量是否为网格模型的最后一层的照明能量，如是，则计算结束；若计算出照明能量的层不是网格模型的最后一层，则将该计算出照明能量的层作为当前层，并重复S4-S8进行计算和判断，直至计算出网格模型中最后一层的照明能量，则计算结束；

所述S7中计算装置根据计算得到的下一层各点的波场值，计算下一层各点的照明能量，包括：先将S5计算得到的下一层的波数域波场值作傅立叶变换至空间域；然后在空间域根据观测系统中的炮检关系，对于每一个作用点，先计算出该点作为中心点时在下一层的相同位置的波场值，然后计算出与其相关的接收点作为中心点时，该点在各区域计算得到的下一层相同位置的波场值，将该点作为中心点时得到的下一层的波场值进行平方后分别和与其相关的接收点作为中心点时该点的下一层波场值的平方相乘并求和，计算所得的结果记作该计算区域中心点位置的照明能量DI，并存储在计算机的硬盘上。

2.如权利要求1所述的一种基于地震波波动方程的照明分析计算方法，其特征在于，所述S4中的计算区域是按照如下方法确定的：

先将所有炮点和检波点都看作是作用点，检波点和炮点重复的位置算作一个作用点，所述计算区域包括以作用点为中心，从中心的左边Maxoffset/DX个网格点开始，至中心的右边Maxoffset/DX个网格点结束的区间；所述Maxoffset表示最大偏移距，所述DX表示网格步长。

3.如权利要求2所述的一种基于地震波波动方程的照明分析计算方法，其特征在于，所述S4中在计算区域内设定一个初始波场，包括：如果当前层为网格模型的第一层，则在波数域内将计算区域的中心位置处实部设为大于0的任意值，计算区域的其他位置处实部设为0，整个计算区域的虚部都取0；如果当前层不是第一层，则在波数域内将该区域根据上一层相同位置的计算区域计算所得的波场值设定为初始波场值。

4.如权利要求3所述的一种基于地震波波动方程的照明分析计算方法，其特征在于，所述S5中，传播算子取为exp(i×DX×kz)，其中，i为虚数单位，DX是网格模型的单元边长，为已知量，kz是垂直波数；

所述垂直波数可以通过下式计算得到：其中，kx为水平波数；

所述水平波数kx可以通过下式计算：kx＝x×2×3.1415926×f，其中，f是地震波的频率，x是当前点的位置，f和x都为已知量。

5.如权利要求1所述的一种基于地震波波动方程的照明分析计算方法，其特征在于，所述S2中将建立的地质模型离散化成网格模型，具体包括：先计算检波点间距的二分之一和炮点间距的二分之一，即RI/2和SI/2，其中，RI是检波点间距，SI是炮点间距；然后取DX作为离散后的网格模型的单元边长，DX的取值方法为：取RI/2和SI/2中较小的数作为DX的值；再以DX为单元将S1中建立的地质模型离散化成网格模型。

6.一种基于地震波波动方程的照明分析计算装置，其特征在于，包括：计算前处理模块、模型网格化模块、读取数据模块、波场值设定模块、波场值计算模块、第一判断模块、照明能量计算模块、第二判断模块和能量分布输出模块；其中，

计算前处理模块，用于根据工区地质构造建立地质模型，并读入观测系统的参数；

模型网格化模块，用于将计算前处理模块建立的地质模型离散化成网格模型；

读取数据模块，用于将网格模型的第一层的网格点数据读入内存中，准备进行计算；

波场值设定模块，用于将炮点或检波点位置作为中心，在计算区域内设定一个波数域的初始波场；

波场值计算模块，用于根据波场值设定模块中设定的初始波场，利用传播算子，计算下一层相应位置处的波场值；

第一判断模块，用于判断当前作为中心的炮点或检波点是否是观测系统中该层的最后一个炮点或检波点；若判断的结果为：当前作为中心的炮点或检波点是观测系统中该层的最后一个炮点或检波点，则丢弃内存中原来存储的那一层的波场值，将计算所得的下一层的波场值存储在内存中；若判断结果为：当前作为中心的炮点或检波点不是观测系统中该层的最后一个炮点或检波点，则将上述波场设定模块中的中心向右偏移至下一个炮点或检波点，并利用波场计算模块和第一判断模块进行计算和判断，直至计算完最后一个炮点或检波点，再丢弃内存中原来存储的波场值，将计算所得的下一层的波场值存放在内存中；

照明能量计算模块，用于根据波场值计算模块计算得到的下一层的各点所有波场值，计算下一层的照明能量；具体地，用于先将计算得到的下一层的波数域波场值作傅立叶变换至空间域；然后根据观测系统中的炮检关系，对于每一个炮点或检波点，先计算出该点作为中心点时在下一层的相同位置的波场值，然后计算出与其相关的接收点作为中心点时，该点在各区域计算得到的下一层相同位置的波场值，将该点作为中心点时得到的下一层的波场值进行平方后分别和与其相关的接收点作为中心点时该点的下一层波场值的平方相乘并求和，计算所得的结果记作该计算区域中心点位置的照明能量DI，并存储在计算机的硬盘上；

第二判断模块，用于判断计算出照明能量的层是否为网格模型的最后一层；若判断结果为：照明能量计算模块中计算出照明能量的层的是网格模型的最后一层，则计算结束；若判断结果为：计算出照明能量的层不是网格模型的最后一层，则将该计算出照明能量的层作为当前层，并利用波场值设定模块、波场值计算模块、第一判断模块、照明能量计算模块和第二判断模块再次计算下一层的照明能量并进行判断，直至计算出网格模型中最后一层的照明能量，则计算结束；

能量分布输出模块，用于输出每一层的能量分布图。

7.如权利要求6所述的一种基于地震波波动方程的照明分析计算装置，其特征在于，所述波场值设定模块，用于设定波场的初始波场值；对于网格模型中第一层的数据，波场值设定模块采用的波场设定方法为：在波数域将作为计算区域中心的作用点实部设为大于0的任意值，其他网格点位置的实部设为0，整个计算区域的虚部都设为0；对于网格模型中其它层的数据，波场值设定模块采用的波场设定方法为：则将该计算区域在波数域根据上一层相同位置的计算区域计算所得的波场值设定为初始波场值。

8.如权利要求7所述的一种基于地震波波动方程的照明分析计算装置，其特征在于，所述波场值计算模块，用于在波数域将设定的初始波场乘以传播算子得到网格模型中下一层相应位置处的波场值；所述传播算子可以取为exp(i×DX×kz)，其中i为虚数单位，DX是网格模型的单元边长，为已知量，kz是垂直波数；

垂直波数可以通过下式计算得到：其中，kx为水平波数；

水平波数kx可以通过下式计算：kx＝x×2×3.1415926×f，其中，f是地震波的频率，x是当前点的位置，f和x都为已知量。