地震三维观测系统的共中心点面元属性的分析方法及装置

摘要

本发明提供了一种地震三维观测系统的CMP面元属性的分析方法及装置，该方法包括：根据计算资源所允许容纳的数据量，将地震三维观测系统的待勘探区域的面元格网划分为多个覆盖次数计算分区；采用多个计算资源并行分别计算多个覆盖次数计算分区中面元的覆盖次数；根据多个覆盖次数计算分区中面元的覆盖次数和计算资源所允许容纳的中点数量，将待勘探区域的面元格网划分为多个中点属性计算分区；采用多个计算资源并行分别计算多个中点属性计算分区的中点属性值，中点属性值包括：炮检距、炮检方位角。避免了串行依次对每个面元的面元属性进行计算，从而提高了分析CMP面元属性的效率、资源利用率。

说明书

技术领域

本发明涉及石油勘探技术领域，特别涉及一种地震三维观测系统的CMP （Common Mid-point，共中心点）面元属性的分析方法及装置。

背景技术

为了合理设计野外观测系统方案，需要对不同观测系统的CMP面元属性进行对 比分析。常规CMP面元属性分析有：覆盖次数、炮检距和炮检方位角。对覆盖次数、 炮检距和炮检方位角的直观和精确描述，是优化观测系统设计的关键因素之一，而 CMP面元属性计算则是实现观测系统面元属性分析的基本技术手段。随着高密度地 震勘探的推广应用，野外生产的接收道数不断增多，百万乃至千万级炮检点的布设越 来越普遍，由于一个面元的面元属性值是对该面元上覆盖的多个关系片计算后的统计 和累加结果，CMP面元属性计算的工作量大、耗时长，成为生产应用的瓶颈。在对 CMP面元属性分析的过程中，由于传统的面元属性计算方法是串行执行的，即依次 计算每个面元的面元属性，因此，在多核硬件环境下，不能充分利用计算资源；由于 串行执行CMP面元属性计算的效率低，而不能满足高密度地震勘探带来海量CMP 面元属性计算的需求，从而导致了分析CMP面元属性时，效率低、资源利用率低的 问题。

发明内容

本发明实施例提供了一种地震三维观测系统的CMP面元属性的分析方法及装 置，解决了现有技术中分析CMP面元属性效率低、资源利用率低的技术问题。

本发明实施例提供了一种地震三维观测系统的CMP面元属性的分析方法，该方 法包括：根据计算资源所允许容纳的数据量，将地震三维观测系统的待勘探区域的面 元格网划分为多个覆盖次数计算分区；采用多个计算资源并行分别计算多个覆盖次数 计算分区中面元的覆盖次数；根据上述多个覆盖次数计算分区中面元的覆盖次数和计 算资源所允许容纳的中点数量，将上述待勘探区域的上述面元格网划分为多个中点属 性计算分区；采用多个计算资源并行分别计算多个中点属性计算分区的中点属性值， 其中，上述中点属性值包括：炮检距、炮检方位角。

在一个实施例中，根据计算资源所允许容纳的数据量，将所述地震三维观测系统 的待勘探区域的面元格网划分为多个覆盖次数计算分区，包括：在大地坐标系下，对 所述待勘探区域的面元格网中的面元进行编号，其中，以上述大地坐标系原点处的面 元为起点，先按照以RMin所在的方向从左到右的顺序，再按照RMax所在的方向从 下到上的顺序对上述面元格网中的面元进行递增编号，上述RMin是上述待勘探区域 横向长度的最小者，上述RMax是上述待勘探区域纵向长度的最大者；

根据以下公式确定覆盖次数计算分区的数量： 其中，上述M是覆盖次数计算分区的数量， $\lambda =\frac{(\mathrm{MaxBinRowInl}+1)\times (\mathrm{MaxBinColCrl}+1)\times \mathrm{FoldDataPerBin}}{\mathrm{DataPerKernel}},$为向上取整符号， NumOfKernel参与计算的计算资源的数量，DataPerKernel为每个计算资源所允许容 纳的数据量，FoldDataPreBin为每个面元的覆盖次数属性项数据，MaxBinRowInl为 上述待勘探区域的纵向最大面元行号，MaxBinColCrl为上述待勘探区域的横向最大 面元列号。

在一个实施例中，根据计算资源所允许容纳的数据量，将所述地震三维观测系统 的待勘探区域的面元格网划分为多个覆盖次数计算分区，包括：通过以下公式确定每 个覆盖次数计算分区包括的面元行数：BinLinePerZone=min(L_max,η)，其中， $\eta =\frac{\mathrm{DataPerKernel}}{\mathrm{FoldDataPerBin}\times (\mathrm{MinBinIDOfLine}+1)},$为向上 取整符号，BinLinePerZone为每个覆盖次数计算分区包括的面元行数，min为取二个 数值中的较小者，DataPerKernel为每个计算资源所允许容纳的数据量，FoldDataPreBin 为每个面元的覆盖次数属性项数据，MaxBinIDOfLine是每个覆盖次数计算分区的行 数和列数中的最大者；MinBinIDOfLine是每个覆盖次数计算分区的行数和列数中的 最小者，NumOfKernel是参与计算的计算资源的数量；通过以下公式确定每个覆盖次 数计算分区的起始面元的编号： BinIDBegin_i=BinLinePerZone×(MinBinIDOfLine+1)×i，其中，上述i是上述覆盖次数 计算分区的编号，i∈[0,M-1]；通过以下公式确定每个覆盖次数计算分区的终止面元 的编号：BinIDEnd_i=BinLinePerZone×(MinBinIDOfLine+1)×(i+1)-1。

在一个实施例中，采用上述多个计算资源并行分别计算上述多个覆盖次数计算分 区中面元的覆盖次数包括：根据每个覆盖次数计算分区的起始面元的编号和终止面元 的编号，确定上述每个覆盖次数计算分区的起始面元和终止面元的坐标值；根据上述 每个覆盖次数计算分区的起始面元和终止面元的坐标值，确定每个覆盖次数计算分区 在相对坐标系下的包络范围，其中，该相对坐标系是面元线方位角为零的大地坐标系； 将上述地震三维观测系统的各个关系片内的物理点的坐标值与上述相对坐标系的坐 标值进行比较，确定上述各个关系片在上述相对坐标系下的包络范围，其中，上述各 个关系片在上述相对坐标系下的包络范围是包含物理点的最小矩形；判断上述覆盖次 数计算分区的包络范围与上述各个关系片的包络范围是否有交集，若有，则计算该关 系片内每对炮点和检波点的中心点的坐标值；判断计算出的中心点的坐标值是否包含 在上述覆盖次数计算分区的面元内，若是，则该面元的覆盖次数加一，否则，不累加。

在一个实施例中，根据上述多个覆盖次数计算分区中面元的覆盖次数和上述计算 资源所允许容纳的中点数量，将上述待勘探区域的上述面元格网划分为上述多个中点 属性计算分区包括：根据以下公式确定中点属性计算分区的数量：其中， 上述N是中点属性计算分区的数量；T是累加每个覆盖次数计算分区中每个面元的覆 盖次数得到上述待勘探区域的覆盖次数；C是所述计算资源所允许容纳的中点数量， 通过以下公式确定上述计算资源所允许容纳的中点数量：其中， $\mu =\frac{T}{\mathrm{NumOfKernel}},\delta =\frac{\mathrm{DataPerKernel}}{\mathrm{DataPerMidPnt}},$DataPreMidPnt为每次覆盖对应的中点属性 项的数据量。

在一个实施例中，根据上述多个覆盖次数计算分区中面元的覆盖次数和上述计算 资源所允许容纳的中点数量，将上述待勘探区域的面元格网划分为上述多个中点属性 计算分区，包括：根据以下公式确定每个中点属性计算分区的起始面元的编号： BinIDBegin_k+1=BinIDEnd_k+1，其中，上述K是上述中点属性计算分区的编号， k∈[0,N-1]，BinIDBegin_k+1为第k+1个中点属性计算分区的起始面元的编号， BinIDEnd_k为第k个中点属性计算分区的终止面元的编号；根据以下公式确定每个中 点属性计算分区的终止面元的编号：BinIDEnd_k+1=BinIDBegin_k+1+R-1，其中， BinIDEnd_k+1为第k+1个中点属性计算分区的终止面元的编号，上述R为上述中点属 性计算分区包括的面元数。

在一个实施例中，采用上述多个计算资源并行分别计算上述多个中点属性计算分 区的中点属性值之后，上述方法还包括：按照面元的编号的顺序，对每个中点属性计 算分区中面元的中点、炮检距和炮检方位角进行排序，其中，面元编号相同的中点、 炮检距和炮检方位角相邻且连续排列。

本发明实施例还提供了一种地震三维观测系统的CMP面元属性的分析装置，该 装置包括：第一划分模块，用于根据计算资源所允许容纳的数据量，将地震三维观测 系统的待勘探区域的面元格网划分为多个覆盖次数计算分区；第一计算模块，用于采 用多个计算资源并行分别计算上述多个覆盖次数计算分区中面元的覆盖次数；第二划 分模块，用于根据上述多个覆盖次数计算分区中面元的覆盖次数和上述计算资源所允 许容纳的中点数量，将上述待勘探区域的上述面元格网划分为多个中点属性计算分 区；第二计算模块，用于采用多个计算资源并行分别计算上述多个中点属性计算分区 的中点属性值，其中，上述中点属性值包括：炮检距、炮检方位角。

在一个实施例中，上述第一划分模块包括：编号单元，用于在大地坐标系下，对 上述待勘探区域的面元格网中的面元进行编号，其中，以上述大地坐标系原点处的面 元为起点，先按照以RMin所在的方向从左到右的顺序，再按照RMax所在的方向从 下到上的顺序对上述面元格网中的面元进行递增编号，上述RMin是上述待勘探区域 横向长度的最小者，上述RMax是上述待勘探区域纵向长度的最大者；覆盖次数计算 分区确定单元，用于根据以下公式确定覆盖次数计算分区的数量： 其中，上述M是覆盖次数计算分区的数量， $\lambda =\frac{(\mathrm{MaxBinRowInl}+1)\times (\mathrm{MaxBinColCrl}+1)\times \mathrm{FoldDataPerBin}}{\mathrm{DataPerKernel}},$为向上取整符号， NumOfKernel参与计算的计算资源的数量，DataPerKernel为每个计算资源所允许容 纳的数据量，FoldDataPreBin为每个面元的覆盖次数属性项数据，MaxBinRowInl为 上述待勘探区域的纵向最大面元行号，MaxBinColCrl为上述待勘探区域的横向最大 面元列号。

在一个实施例中，上述第一划分模块还包括：面元行数确定单元，用于根据以下 公式确定每个覆盖次数计算分区包括的面元行数：BinLinePerZone=min(L_max,η)，其 中，$\eta =\frac{\mathrm{DataPerKernel}}{\mathrm{FoldDataPerBin}\times (\mathrm{MinBinIDOfLine}+1)},$为 向上取整符号，BinLinePerZone为每个覆盖次数计算分区包括的面元行数，min为取 二个数值中的较小者，DataPerKernel为每个计算资源所允许容纳的数据量， FoldDataPreBin为每个面元的覆盖次数属性项数据，MaxBinIDOfLine是每个覆盖次 数计算分区的行数和列数中的最大者，MinBinIDOfLine是每个覆盖次数计算分区的 行数和列数中的最小者，NumOfKernel是参与计算的计算资源的数量；第一起始面元 确定单元，用于通过以下公式确定每个覆盖次数计算分区的起始面元的编号： BinIDBegin_i=BinLinePerZone×(MinBinIDOfLine+1)×i，其中，上述i是上述覆盖次数 计算分区的编号，i∈[0,M-1]；第一终止面元确定单元，用于通过以下公式确定每个 覆盖次数计算分区的终止面元的编号： BinIDEnd_i=BinLinePerZone×(MinBinIDOfLine+1)×(i+1)-1。

在一个实施例中，上述第一计算模块包括：坐标确定单元，用于根据每个覆盖次 数计算分区的起始面元的编号和终止面元的编号，确定上述每个覆盖次数计算分区的 起始面元和终止面元的坐标值；第一包络范围确定单元，用于根据上述每个覆盖次数 计算分区的起始面元和终止面元的坐标值，确定每个覆盖次数计算分区在相对坐标系 下的包络范围，其中，上述相对坐标系是面元线方位角为零的大地坐标系；第二包络 范围确定单元，用于将上述地震三维观测系统的各个关系片内的物理点的坐标值与上 述相对坐标系的坐标值进行比较，确定上述各个关系片在上述相对坐标系下的包络范 围，其中，上述各个关系片在上述相对坐标系下的包络范围是包含物理点的最小矩形； 计算单元，用于判断上述覆盖次数计算分区的包络范围与上述各个关系片的包络范围 是否有交集，若有，则计算该关系片内每对炮点和检波点的中心点的坐标值；计数单 元，用于判断计算出的中心点的坐标值是否包含在上述覆盖次数计算分区的面元内， 若是，则该面元的覆盖次数加一，否则，不累加。

在一个实施例中，第二划分模块包括：中点属性计算分区确定单元，用于根据以 下公式确定中点属性计算分区的数量：其中，上述N是中点属性计算分区 的数量；累加单元，用于累加每个覆盖次数计算分区中每个面元的覆盖次数得到上述 待勘探区域的覆盖次数T；中点数量确定单元，用于通过以下公式确定上述计算资源 所允许容纳的中点数量C：其中，DataPreMidPnt为每次覆盖对应的中点属性项的数据量。

在一个实施例中，上述第二划分模块还包括：第二起始面元确定单元，用于根据 以下公式确定每个中点属性计算分区的起始面元的编号： BinIDBegin_k+1=BinIDEnd_k+1，其中，上述K是上述中点属性计算分区的编号， k∈[0,N-1]，BinIDBegin_k+1为第k+1个中点属性计算分区的起始面元的编号， BinIDEnd_k为第k个中点属性计算分区的终止面元的编号；第二终止面元确定单元， 用于根据以下公式确定每个中点属性计算分区的终止面元的编号： BinIDEnd_k+1=BinIDBegin_k+1+R-1，其中，BinIDEnd_k+1为第k+1个中点属性计算分 区的终止面元的编号，上述R为上述中点属性计算分区包括的面元数。

在一个实施例中，上述装置还包括：排序模块，用于按照面元的编号的顺序，对 每个中点属性计算分区中面元的中点、炮检距和炮检方位角进行排序，其中，面元编 号相同的中点、炮检距和炮检方位角相邻且连续排列。

在本发明实施例中，通过将待勘探区域的面元格网划分成若干个独立的、粒度适 中的计算分区，并采用若干个计算资源并行分别计算若干个计算分区的面元属性，避 免了串行依次对每个面元的面元属性进行计算，从而可以提高计算资源的利用率，尤 其有利于多核和多机环境下的并行计算；还可以提高面元属性的计算效率，以满足大 型探区、高密度地震勘探的CMP面元属性计算的需求，进而实现了分析CMP面元 属性时，提高效率、提高资源利用率的技术效果。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，并不 构成对本发明的限定。在附图中：

图1是本发明实施例提供的一种地震三维观测系统的CMP面元属性的分析方法 的流程图；

图2是本发明实施例提供的一种对面元格网中面元进行编号的示意图；

图3是本发明施例提供的一种将面元格网划分为若干个覆盖次数计算分区的示 意图；

图4是本发明实施例提供的一种对中点属性计算分区中面元的属性值进行排序 的示意图；

图5是本发明实施例提供的一种地震三维观测系统的CMP面元属性的分析装置 的结构框图；

图6是本发明实施例提供的一种第一划分模块的结构框图；

图7是本发明实施例提供的一种第一计算模块的结构框图；

图8是本发明实施例提供的一种第二划分模块的结构框图；

图9是本发明实施例提供的另一种地震三维观测系统的CMP面元属性的分析装 置的结构框图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，下面结合实施方式和附图， 对本发明做进一步详细说明。在此，本发明的示意性实施方式及其说明用于解释本发 明，但并不作为对本发明的限定。

在本发明实施例中，提供了一种计算地震三维观测系统的CMP面元属性的方法， 如图1所示，该方法包括：

步骤101：根据计算资源所允许容纳的数据量，将地震三维观测系统的待勘探区 域的面元格网划分为多个覆盖次数计算分区；

步骤102：采用多个计算资源并行分别计算上述多个覆盖次数计算分区中面元的 覆盖次数；

步骤103：根据上述多个覆盖次数计算分区中面元的覆盖次数和上述计算资源所 允许容纳的中点数量，将上述待勘探区域的上述面元格网划分为多个中点属性计算分 区；

步骤104：采用多个计算资源并行分别计算上述多个中点属性计算分区的中点属 性值，其中，上述中点属性值包括：炮检距、炮检方位角。

由图1所示流程可以得知，本发明实施例中，首先，根据计算资源所允许容纳的 数据量，将地震三维观测系统的待勘探区域的面元格网划分为多个覆盖次数计算分 区，以便采用多个计算资源并行分别计算上述多个覆盖次数计算分区中面元的覆盖次 数，即，每个计算资源负责一个覆盖次数计算分区中面元的覆盖次数的计算，多个计 算资源并行计算；然后，根据上述多个覆盖次数计算分区中面元的覆盖次数和上述计 算资源所允许容纳的中点数量，将上述待勘探区域的面元格网划分为多个中点属性计 算分区，以便采用多个计算资源并行分别计算上述多个中点属性计算分区的中点属性 值，其中，该中点属性值包括：炮检距、炮检方位角，即，每个计算资源负责一个中 点属性计算分区的中点属性值的计算，多个计算资源并行计算。

由此可见，本发明实施例通过将待勘探区域的面元格网划分成若干个独立的、粒 度适中的计算分区，并采用若干个计算资源并行分别计算若干个计算分区的面元属 性，各计算分区间并行计算，计算分区内顺序计算，结果面元属性值既不重复，也不 会遗漏，避免了串行依次对每个面元的面元属性进行计算，从而可以提高计算资源的 利用率，尤其有利于多核和多机环境下的并行计算；还可以提高面元属性的计算效率， 以满足大型探区、高密度地震勘探的CMP面元属性计算的需求，进而实现了分析CMP 面元属性时，提高效率、提高资源利用率的技术效果。

具体实施时，先将设计好的地震三维观测系统方案布设到待勘探区域，得到炮点、 检波点的位置及排列片关系，再将上述地震三维观测系统的待勘探区域的上述面元格 网划分为多个覆盖次数计算分区。具体的，在将上述地震三维观测系统的待勘探区域 的上述面元格网划分为多个覆盖次数计算分区之前，可以通过以下方式将待勘探区域 划分成面元格网，例如，如图2所示，在以Ox、Oy为观测系统原点的大地坐标下， 可以将待勘探区域划分成由宽为w，高为h的矩形面元组成的面元格网，总行数为 row，总列数为col，本发明实施例对矩形面元的大小不做限定，可以根据具体需求 情况，确定矩形面元的宽和高。

在将待勘探区域划分成面元格网后，为了确定合理的、适中的覆盖次数计算分区 数，在将面元格网划分成多个覆盖次数计算分区时，可以通过以下方式实现：首先， 对待勘探区域的面元格网中的面元进行编号，然后，再确定覆盖次数计算分区的数量。

具体的，可以通过以下方式对待勘探区域的面元格网中的面元进行编号：在大地 坐标系下，对待勘探区域的面元格网中的面元进行编号，其中，以上述大地坐标系原 点处的面元为起点，先按照以RMin所在的方向从左到右的顺序，再按照RMax所在 的方向从下到上的顺序对上述面元格网中的面元进行递增编号，即按照如图2中虚线 箭头所示的方向进行编号，其中，上述RMin是上述待勘探区域横向长度的最小者， 上述RMax是上述待勘探区域纵向长度的最大者。

具体的，可以根据以下公式确定覆盖次数计算分区的数量：

$\lambda =\frac{(\mathrm{MaxBinRowInl}+1)\times (\mathrm{MaxBinColCrl}+1)\times \mathrm{FoldDataPerBin}}{\mathrm{DataPerKernel}}---\left(2\right)$

其中，M是覆盖次数计算分区的数量，为向上取整符号，NumOfKernel参与 计算的计算资源的数量，DataPerKernel为每个计算资源所允许容纳的数据量， FoldDataPreBin为每个面元的覆盖次数属性项数据，MaxBinRowInl为上述待勘探区 域的纵向最大面元行号，MaxBinColCrl为上述待勘探区域的横向最大面元列号。

在确定出覆盖次数计算分区的数量之后，为了快速地、准确地确定每个覆盖次数 计算分区包含的面元，进而确定每个覆盖次数计算分区的范围对应的坐标值，在本实 施例中，分别确定了每个覆盖次数计算分区的面元行数、起始面元的编号、终止面元 的编号。

具体的，可以通过以下公式确定每个覆盖次数计算分区包括的面元行数：

BinLinePerZone=min(L_max,η)（3）

$\eta =\frac{\mathrm{DataPerKernel}}{\mathrm{FoldDataPerBin}\times (\mathrm{MinBinIDOfLine}+1)}---\left(4\right)$

其中，为向上取整符号，BinLinePerZone为每个覆盖次数计算分区包括的面 元行数，min为取二个数值中的较小者，DataPerKernel为每个计算资源所允许容纳的 数据量，FoldDataPreBin为每个面元的覆盖次数属性项数据，MaxBinIDOfLine是每 个覆盖次数计算分区的行数和列数中的最大者，MinBinIDOfLine是每个覆盖次数计 算分区的行数和列数中的最小者，NumOfKernel是参与计算的计算资源的数量。

具体的，可以通过以下公式确定每个覆盖次数计算分区的起始面元的编号：

BinIDBegin_i=BinLinePerZone×(MinBinIDOfLine+1)×i     （6）

其中，上述i是上述覆盖次数计算分区的编号，i∈[0,M-1]。

通过以下公式确定每个覆盖次数计算分区的终止面元的编号：

BinIDEnd_i=BinLinePerZone×(MinBinIDOfLine+1)×(i+1)-1     （7）

确定出每个覆盖次数计算分区的面元行数、起止面元的编号后，如图3所示，可 以快速地、准确地确定覆盖次数计算分区0、1、2、3包含的面元，根据面元在大地 坐标系下的坐标值，可以进而确定每个覆盖次数计算分区的包络范围对应的坐标值。

确定出每个覆盖次数计算分区的面元行数、起止面元的编号之后，为了准确地、 高效地计算出每个计算分区中面元的覆盖次数，在本实施例中，采用多个计算资源并 行分别计算上述多个覆盖次数计算分区的覆盖次数，具体实施时，首先，确定每个覆 盖次数计算分区的起始面元和终止面元的坐标值，再确定每个覆盖次数计算分区在相 对坐标系下的包络范围，以及确定各个关系片在上述相对坐标系下的包络范围，最后， 判断上述覆盖次数计算分区的包络范围与上述各个关系片的包络范围是否有交集，若 有，则计算该关系片内每对炮点和检波点的中点的坐标值，并计算中点的坐标值所在 的覆盖次数计算分区的面元的覆盖次数。

具体的，确定每个覆盖次数计算分区的起始面元和终止面元的坐标值的过程包 括：根据每个覆盖次数计算分区的起始面元的编号和终止面元的编号，确定每个覆盖 次数计算分区的起始面元和终止面元的坐标值。

具体的，确定每个覆盖次数计算分区在相对坐标系下的包络范围的过程包括：根 据每个覆盖次数计算分区的起始面元和终止面元的坐标值，确定每个覆盖次数计算分 区在相对坐标系下的包络范围，其中，该相对坐标系是面元线方位角为零的大地坐标 系；即为了后续的计算便捷，将覆盖次数计算分区在大地坐标系下的坐标值转换为在 相对坐标系下的坐标值，例如，大地坐标系下任意点M(E,N)到相对坐标系下M(x,y) 的转换关系为：

$\left(\begin{array}{c}x=(N-\mathrm{Ox})\times \sin \theta +(E-\mathrm{Oy})\times \cos \theta \\ y=-(N-\mathrm{Ox})\times \cos \theta +(E-\mathrm{Oy})\times \sin \theta \end{array}\right)---\left(8\right)$

上式中Ox、Oy为大地坐标系的原点，θ为面元线方位角。

具体的，确定各个关系片在上述相对坐标系下的包络范围的过程包括：将上述地 震三维观测系统的各个关系片内的物理点的坐标值与上述相对坐标系的坐标值进行 比较，确定上述各个关系片在上述相对坐标系下的包络范围，其中，上述各个关系片 在上述相对坐标系下的包络范围是包含物理点的最小矩形；即，将各个关系片内的物 理点在大地坐标系下的坐标值转换到相对坐标系下，并将物理点在相对坐标系下的坐 标值与相对坐标系的坐标值进行比较，确定包含物理点的最小矩形的包络范围。

具体的，计算关系片内每对炮点和检波点的中点的坐标值的过程包括：判断上述 覆盖次数计算分区的包络范围与上述各个关系片的包络范围是否有交集，若有，则计 算该关系片内每对炮点和检波点的中点的坐标值，例如，可以采用以下公式计算每对 炮点和检波点的中点的坐标值：

$\left(\begin{array}{c}x=[(\mathrm{Sx}+\mathrm{Rx})\times 100/2]\\ y=[(\mathrm{Sy}+\mathrm{Ry})\times 100/2]\end{array}---\left(9\right)\right)$

其中，为取整符号，(Sx,Sy)为炮点在大地坐标系下的坐标值，(Rx,Ry)为检波 点在大地坐标系下的坐标值。

具体的，计算覆盖次数计算分区中面元的覆盖次数的过程包括：判断计算出的中 点的坐标值是否包含在上述覆盖次数计算分区的面元内，若是，则该面元的覆盖次数 加一，否则，不累加。

在计算出覆盖次数计算分区中面元的覆盖次数之后，为了快速、高效地计算上述 待勘探区域的中点属性值，在本实施例中，提供了一种有效确定中点属性计算分区数 的方法，即，根据上述多个覆盖次数计算分区中面元的覆盖次数和上述计算资源所允 许容纳的中点数量，将上述待勘探区域的上述面元格网划分为多个中点属性计算分 区。

具体实施时，可以根据以下公式确定中点属性计算分区的数量：

其中，N是中点属性计算分区的数量。

具体的，T是待勘探区域的覆盖次数，可以累加每个覆盖次数计算分区中每个面 元的覆盖次数得到上述待勘探区域的覆盖次数T；例如，可以采用以下公式来计算待 勘探区域的覆盖次数T，

$T=\underset{i=0}{\overset{n}{\Sigma }}\underset{j=0}{\overset{m}{\Sigma }}{F}_{\mathrm{ij}}---\left(11\right)$

其中，i为面元号，i∈[0,n]，n为最大的面元编号，j为面元的第j次覆盖，m为 第i个面元的最高覆盖次数，F_ij为第i个面元第j次覆盖。

C是计算资源所允许容纳的中点数量，具体的，可以通过以下公式确定上述计算 资源所允许容纳的中点数量：

$\mu =\frac{T}{\mathrm{NumOfKernel}}---\left(13\right)$

$\delta =\frac{\mathrm{DataPerKernel}}{\mathrm{DataPerMidPnt}}---\left(14\right)$

其中，DataPreMidPnt为每次覆盖对应的中点属性项的数据量。

确定出中点属性计算分区的数量之后，为了快速、准确地确定每个中点属性计算 分区包含的面元，在本实施例中，可以确定每个中点属性计算分区包括的面元。具体 实施时，分别确定每个中点属性计算分区的起始面元的编号和终止面元的编号。

具体的，可以根据以下公式确定每个中点属性计算分区的起始面元的编号：

BinIDBegin_k+1=BinIDEnd_k+1     （15）

其中，上述K是上述中点属性计算分区的编号，k∈[0,N-1]，BinIDBegin_k+1为 第k+1个中点属性计算分区的起始面元的编号，BinIDEnd_k为第k个中点属性计算分 区的终止面元的编号；

具体的，可以根据以下公式确定每个中点属性计算分区的终止面元的编号：

BinIDEnd_k+1=BinIDBegin_k+1+R-1     （16）

其中，BinIDEnd_k+1为第k+1个中点属性计算分区的终止面元的编号，上述R为 上述中点属性计算分区包括的面元数。另外，还可以计算每个中点属性计算分区的覆 盖次数，例如，可以根据以下公式来计算：其中，U为相邻且连续的R个 面元的覆盖次数之和，F_i为其中第i个面元的覆盖次数，U≤S且U取最接近S的整数。

确定出每个中点属性计算分区包含的面元后，在计算每个中点属性计算分区的中 点属性值的过程中，可以采用以下公式计算得到炮检距和炮检方位角：

$D=\sqrt{{(\mathrm{Sx}-\mathrm{Rx})}^2+{(\mathrm{Sy}-\mathrm{Ry})}^2}---\left(17\right)$

β=atan2((Sx-Rx),(Sy-Ry))     （18）

其中，D为炮检距，β为炮检方位角，(Sx,Sy)为炮点在大地坐标系下的坐标值， (Rx,Ry)为检波点在大地坐标系下的坐标值。

采用上述多个计算资源并行分别计算上述多个中点属性计算分区的中点属性值 之后，为了得到上述待勘探区域的整体的面元属性，在本实施例中，上述方法还包括： 按照面元的编号的顺序，对每个中点属性计算分区中面元的中点、炮检距和炮检方位 角进行排序，其中，面元编号相同的中点、炮检距和炮检方位角相邻且连续排列。例 如，如图4所示，图中的方格部分显示有数字9的6个方格表示为中点数量为9的6 个面元，显示有数字10的9个方格表示为中点数量为10的9个面元，以图中相邻的 中点数量分别为9和10的两个面元为例，对其的中点进行排序，则如图4中序列号 部分所示，中点数量为9的面元的中点相邻且连续排列，即9个中点编号为1至9 依次排列，对相邻的中点数量为10的面元的中点排序则是10个中点编号为10至19 相邻且连续排列，并整体与前一面元的中点相邻排列。将上述各个中点属性计算分区 中面元的中点、炮检距和炮检方位角排序后，输出到存储设备中进行存储，即可得到 整个观测系统的中点、炮检距和炮检方位角。

具体实施时，计算资源可以是计算机、多核计算机等具备计算功能的设备。

以下结合具体实施例对上述地震三维观测系统的CMP面元属性的分析方法进行 详细的描述，该方法包括如下步骤：

步骤1：将设计好的地震三维观测系统布设到待勘探区域后，得到炮点、检波点 位置及炮检关系（排列关系片）。例如，某待勘探区域设计观测系统为：待勘探区域 范围纵向（inline）长29540m，横向（crossline）长26240m，面元大小20m×20m， 面元线方位角为336度，总激发点数25438，总接收点数92773，关系片个数16218。

步骤2：根据面元线方位角参数将大地坐标系下的观测系统旋转到面元线方位角 为零的相对坐标系下。

步骤3：将待勘探区域划分成由矩形面元组成的面元格网，并对面元格网中的面 元进行编号；

以宽w=20m，高h=20m的矩形单元将待勘探区域划分为总行数row=1477，总列 数col=1312的面元格网；从大地坐标系原点处的面元开始，先以RMin所在方向按 从左到右的顺序，再以RMax方向按从下到上的顺序对面元递增编号，即，按照如图 2中虚线箭头所示的方向进行编号，其中，RMin代表待勘探区域inline长度中的最小 者所在方向，RMax代表待勘探区域crossline长度中的最大者所在方向，面元号取值 范围是[0,1937823]，即起始面元编号0，最大面元编号为1937823。

步骤4：利用上述观测系统的关系片集合P，先将各个关系片中物理点在大地坐 标系下的坐标值转换为在相对坐标系下的坐标值，再将物理点在相对坐标系下的坐标 值与相对坐标系下的坐标值进行比较，确定各个关系片Patch的包络范围，其中，各 个关系片的包络范围为包含物理点的最小矩形。

步骤5：按每个计算资源所容纳的数据量，将面元格网划分为多个覆盖次数计算 分区；

如，当NumOfKernel=4、DataPerKernel=200兆字节、FoldDataPreBin=8字节、 MaxBinRowInl=1476、MaxBinColCrl=1311，代入公式（2）计算得：

$\lambda =\frac{(1476+1)\times (1311+1)\times 8}{209715200}=0.07$

综合覆盖次数分区数的取值条件公式（1），M取值为4。

当NumOfKernel=4、MaxBinIDOfLine=1476、MinBinIDOfLine=1311时，代入公 式（5）和（4）得：综合计算公式（3）得每个分区的面元行数BinLinePerZone=370行。

所述覆盖次数计算分区的范围由起始面元号和终止面元号唯一决定，4个分区的 起始面元号和终止面元号分别由公式（6）和（7）计算，分别得到[0，485439]、[485440， 970879]、[970880，1456319]、[1456320，1937823]，覆盖次数计算分区如图3所示。

步骤6：利用多个计算资源和观测系统关系片集合得到炮点和检波点的中点，并 根据中点得到各面元的覆盖次数；

步骤6所述利用多个计算资源计算炮点和检波点的中点，即将相邻的4个覆盖次 数计算分区分别分配给4个计算资源同时并行进行计算，每个计算资源负责计算一个 覆盖次数计算分区，相同时间段内有且只有4个相邻覆盖次数计算分区参与计算；

所述的利用关系片得到炮点和检波点的中点是，将覆盖次数计算分区的起始、终 止面元在大地坐标系下的坐标值转换为在相对坐标系下的坐标值，并根据相对坐标系 下起始、终止面元的坐标值得到覆盖次数计算分区包络范围Rect；遍历关系片集合P 内每个关系片Patch，利用步骤4所得关系片Patch包络矩形与Rect求交，如果相交， 则计算该关系片中每对炮点、检波对的中点的坐标值，否则放弃该关系片。

所述的根据中点得到各面元覆盖次数是：如图2所示，如果炮点S和检波点R 的中点M的坐标值落在由起始面元号和终止面元号决定的覆盖次数计算分区中某面 元Bin内，则中点位置有效，面元Bin的覆盖次数累加一次；否则中点位置无效，面 元Bin的覆盖次数不予累加。面元Bin最终得到的中点总数即该面元的覆盖次数。

步骤7：用步骤6所得面元覆盖次数结果统计出观测系统全部面元覆盖次数总和 T，按照每个计算资源所能容纳的中点数量和各面元覆盖次数将面元格网划分为多个 中点属性计算分区，并求取各中点属性计算分区所对应的起止面元号；

步骤7所述的统计观测系统范围覆盖次数总和是，按公式（11）将所有面元的覆 盖次数相加得到T=92105054次。

如，当NumOfKerne=4、DataPerKernel=200兆字节、DataPerMidPnt=20字节时， 代入公式（13）和（14）得综合取值条件公式（12）得，每个计算资源所 能容纳中点数量C取值为10485760个。根据求得的中心点数量C，计算中点属性计 算分区个数

所述中点属性计算分区的范围由起始、终止面元号唯一决定，9个分区的起始、 终止面元号由公式（15）、（16）计算，分别得到[0,257890]、[257891,454270]、 [454271,650526]、[650527,858516]、[858517,1070146]、[1070147,1268183]、 [1268184,1483795]、[1483796,1701043]、[1701044,1937823]。

步骤8：利用多个计算资源和中点属性计算分区得到中点属性值，并对每个中点 属性计算分区中面元的中点、炮检距和炮检方位角排序，以输出到介质设备；

采用公式（17）（18）计算得到炮检距和炮检方位角；

步骤8所述的对每个中点属性计算分区中面元的中点、炮检距和炮检方位角排 序，是以面元号为关键字，对同一中点属性计算分区中面元的中点、炮检距和炮检方 位角进行排序。如图4所示，图中的方格部分显示有数字9的6个方格表示为中点数 量为9的6个面元，显示有数字10的9个方格表示为中点数量为10的9个面元，以 图中相邻的中点数量分别为9和10的两个面元为例，对其的中点进行排序，则如图 4中序列号部分所示，中点数量为9的面元的中点相邻且连续排列，即9个中点编号 为1至9依次排列，对相邻的中点数量为10的面元的中点排序则是10个中点编号为 10至19相邻且连续排列，并整体与前一面元的中点相邻排列。并把所有中点属性计 算分区中面元的中点、炮检距和炮检方位角排序后，输出到存储设备，即可得到整个 观测系统的中点、炮检距和炮检方位角。

基于同一发明构思，本发明实施例中还提供了一种计算地震三维观测系统的 CMP面元属性的装置装置，如下面的实施例所述。由于计算地震三维观测系统的CMP 面元属性的装置解决问题的原理与一种计算地震三维观测系统的CMP面元属性的方 法相似，因此一种计算地震三维观测系统的CMP面元属性的装置的实施可以参见一 种计算地震三维观测系统的CMP面元属性的方法的实施，重复之处不再赘述。以下 所使用的，术语“单元”或者“模块”可以实现预定功能的软件和/或硬件的组合。尽管以 下实施例所描述的装置较佳地以软件来实现，但是硬件，或者软件和硬件的组合的实 现也是可能并被构想的。

图5是本发明实施例的一种计算地震三维观测系统的CMP面元属性的装置的一 种结构框图，如图5所示，该装置包括：第一划分模块501、第一计算模块502、第 二划分模块503、第二计算模块504，下面对该结构进行说明。

第一划分模块501，用于根据计算资源所允许容纳的数据量，将地震三维观测系 统的待勘探区域的面元格网划分为多个覆盖次数计算分区；第一计算模块502，与上 述第一划分模块501连接，用于采用多个计算资源并行分别计算上述多个覆盖次数计 算分区中面元的覆盖次数；

第二划分模块503，与上述第一计算模块502连接，用于根据上述多个覆盖次数 计算分区中面元的覆盖次数和上述计算资源所允许容纳的中点数量，将上述待勘探区 域的上述面元格网划分为多个中点属性计算分区；

第二计算模块504，与上述第二划分模块503连接，用于采用多个计算资源并行 分别计算上述多个中点属性计算分区的中点属性值，其中，上述中点属性值包括：炮 检距、炮检方位角。

为了确定合理的、适中的覆盖次数计算分区数，在将面元格网划分成多个覆盖次 数计算分区的过程中，如图6所示，上述第一划分模块501包括编号单元5011和覆 盖次数计算分区确定单元501，以下对该结构进行详细说明。

编号单元5011，用于在大地坐标系下，对上述待勘探区域的面元格网中的面元 进行编号，其中，以上述大地坐标系原点处的面元为起点，先按照以RMin所在的方 向从左到右的顺序，再按照RMax所在的方向从下到上的顺序对上述面元格网中的面 元进行递增编号，上述RMin是上述待勘探区域横向长度的最小者，上述RMax是上 述待勘探区域纵向长度的最大者。

覆盖次数计算分区确定单元5012，与上述编号单元5011连接，用于根据以下公 式确定覆盖次数计算分区的数量：其中，M 是覆盖次数计算分区的数量， $\lambda =\frac{(\mathrm{MaxBinRowInl}+1)\times (\mathrm{MaxBinColCrl}+1)\times \mathrm{FoldDataPerBin}}{\mathrm{DataPerKernel}},$为向上取整符号， NumOfKernel参与计算的计算资源的数量，DataPerKernel为每个计算资源所允许容 纳的数据量，FoldDataPreBin为每个面元的覆盖次数属性项数据，MaxBinRowInl为 上述待勘探区域的纵向最大面元行号，MaxBinColCrl为上述待勘探区域的横向最大 面元列号。

为了快速地、准确地确定每个计算分区包含的面元，进而确定每个计算分区的范 围对应的坐标值，在本实施例中，上述第一划分模块501还包括面元行数确定单元 5013、第一起始面元确定单元5014和第一终止面元确定单元5015，以下对该结构进 行详细说明。

面元行数确定单元5013，用于根据以下公式确定每个覆盖次数计算分区包括的 面元行数：BinLinePerZone=min(L_max,η)，其中， $\eta =\frac{\mathrm{DataPerKernel}}{\mathrm{FoldDataPerBin}\times (\mathrm{MinBinIDOfLine}+1)},$为向上 取整符号，BinLinePerZone为每个覆盖次数计算分区包括的面元行数，min为取二个 数值中的较小者，DataPerKernel为每个计算资源所允许容纳的数据量，FoldDataPreBin 为每个面元的覆盖次数属性项数据，MaxBinIDOfLine是每个覆盖次数计算分区的行 数和列数中的最大者，MinBinIDOfLine是每个覆盖次数计算分区的行数和列数中的 最小者，NumOfKernel是参与计算的计算资源的数量。

第一起始面元确定单元5014，与上述面元行数确定单元5013连接，用于通过以 下公式确定每个覆盖次数计算分区的起始面元的编号： BinIDBegin_i=BinLinePerZone×(MinBinIDOfLine+1)×i，其中，上述i是上述覆盖次数 计算分区的编号，i∈[0,M-1]。

第一终止面元确定单元5015，与上述面元行数确定单元5013连接，用于通过以 下公式确定每个覆盖次数计算分区的终止面元的编号： BinIDEnd_i=BinLinePerZone×(MinBinIDOfLine+1)×(i+1)-1。

为了准确地、高效地计算出每个计算分区的覆盖次数，在本实施例中，如图7 所示，上述第一计算模块502包括坐标确定单元5021、第一包络范围确定单元5022、 第二包络范围确定单元5023、计算单元50244和计数单元5025，以下对该结构进行 详细说明。

坐标确定单元5021，用于根据每个覆盖次数计算分区的起始面元的编号和终止 面元的编号，确定每个覆盖次数计算分区的起始面元和终止面元的坐标值。

第一包络范围确定单元5022，与上述坐标确定单元5021连接，用于根据每个覆 盖次数计算分区的起始面元和终止面元的坐标值，确定每个覆盖次数计算分区在相对 坐标系下的包络范围，其中，上述相对坐标系是面元线方位角为零的大地坐标系。

第二包络范围确定单元5023，与上述第五确定单元5022连接，用于将上述地震 三维观测系统的各个关系片内的物理点的坐标值与上述相对坐标系的坐标值进行比 较，确定上述各个关系片在上述相对坐标系下的包络范围，其中，上述各个关系片在 上述相对坐标系下的包络范围是包含物理点的最小矩形。

计算单元5024，与上述第二包络范围确定单元5023连接，用于判断上述覆盖次 数计算分区的包络范围与上述各个关系片的包络范围是否有交集，若有，则计算该关 系片内每对炮点和检波点的中心点的坐标值。

计数单元5025，与上述计算单元5024连接，用于判断计算出的中心点的坐标值 是否包含在上述覆盖次数计算分区的面元内，若是，则该面元的覆盖次数加一，否则， 不累加。

为了快速、高效地计算上述待勘探区域的中点属性值，在本实施例中，如图8 所示，第二划分模块503包括中点属性计算分区确定单元5031、累加单元5032和中 点数量确定单元5033，以下对该结构进行详细说明。

中点属性计算分区确定单元5031，用于根据以下公式确定中点属性计算分区的 数量：其中，N是中点属性计算分区的数量。

累加单元5032，用于累加每个覆盖次数计算分区中每个面元的覆盖次数得到上 述待勘探区域的覆盖次数T。

中点数量确定单元5033，与上述累加单元5032连接，用于通过以下公式确定上 述计算资源所允许容纳的中点数量C：其中，DataPreMidPnt为每次覆盖对应的中点属性项的数据量。

为了快速、准确地确定每个中点属性计算分区包含的面元，在本实施例中，上述 第二划分模块503还包括第二起始面元确定单元5034和第二终止面元确定单元 5035，以下对该结构进行详细说明。

第二起始面元确定单元5034，用于根据以下公式确定每个中点属性计算分区的 起始面元的编号：BinIDBegin_k+1=BinIDEnd_k+1，其中，上述K是上述中点属性计算 分区的编号，k∈[0,N-1]，BinIDBegin_k+1为第k+1个中点属性计算分区的起始面元 的编号，BinIDEnd_k为第k个中点属性计算分区的终止面元的编号；

第二终止面元确定单元5035，与上述第二起始面元确定单元5034连接，用于根 据以下公式确定每个中点属性计算分区的终止面元的编号： BinIDEnd_k+1=BinIDBegin_k+1+R-1，其中，BinIDEnd_k+1为第k+1个中点属性计算分 区的终止面元的编号，上述R为上述中点属性计算分区包括的面元数。

为了得到上述待勘探区域的整体的面元属性，在本实施例中，如图9所示，上述 装置还包括：排序模块505，用于按照面元的编号的顺序，对每个中点属性计算分区 中面元的中点、炮检距和炮检方位角进行排序，其中，面元编号相同的中点、炮检距 和炮检方位角相邻且连续排列。

在本发明实施例中，通过将待勘探区域的面元格网划分成若干个独立的、粒度适 中的计算分区，并采用若干个计算资源并行分别计算若干个计算分区的面元属性，避 免了串行依次对每个面元的面元属性进行计算，从而可以提高计算资源的利用率，尤 其有利于多核和多机环境下的并行计算；还可以提高面元属性的计算效率，以满足大 型探区、高密度地震勘探的CMP面元属性计算的需求，进而实现了分析CMP面元 属性时，提高效率、提高资源利用率的技术效果。

显然，本领域的技术人员应该明白，上述的本发明实施例的各模块或各步骤可以 用通用的计算装置来实现，它们可以集中在单个的计算装置上，或者分布在多个计算 装置所组成的网络上，可选地，它们可以用计算装置可执行的程序代码来实现，从而， 可以将它们存储在存储装置中由计算装置来执行，并且在某些情况下，可以以不同于 此处的顺序执行所示出或描述的步骤，或者将它们分别制作成各个集成电路模块，或 者将它们中的多个模块或步骤制作成单个集成电路模块来实现。这样，本发明实施例 不限制于任何特定的硬件和软件结合。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技 术人员来说，本发明实施例可以有各种更改和变化。凡在本发明的精神和原则之内， 所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。