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

摘要

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

说明书

技术领域

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

背景技术

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

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

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

发明内容

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

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

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\times a_{53}^2+a_{52}^2\times a_{53}^2+a_{54}^2\times a_{53}^2+a_{55}^2\times a_{53}^2.$

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

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

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

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

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

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

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

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