CN101930079A - 一种地震勘探相关/叠加数据处理方法 - Google Patents

Abstract

本发明涉及一种地震勘探相关/叠加数据处理方法，(1)快速傅里叶变换算法与哈特利变换性质的结合混合算法；(2)FFT-Hartly混合算法中分段与求和；(3)多核、多CPU并行处理的线程安排；通过在IBM 3755服务器(4核)上进行测试可以得到如下指标：在不加载噪声控制情况下，2毫秒采样率，24秒扫描信号，32秒地震采集，相关/叠加运算后得到8秒记录长度的情况下，处理4000道时间不大于3秒，处理8000道不大于4秒水平，对于加载噪声控制部分时(分段加权叠加)运算8000道(每道记录长度与前面相同)不大于7秒。

Description

一种地震勘探相关/叠加数据处理方法

技术领域

本发明涉及一种大型地震仪器野外地震勘探数据采集快速相关/叠加数据处理方法。

背景技术

利用可控震源作为地震采集信号激发源，进行野外地震勘探数据采集施工已经非常普遍。对采集的地震数据进行相关/叠加运算与处理的最常使用的方法是通过硬件的手段制作相关运算器与叠加运算器，进行实时的运算处理。以国外某种大型地震仪器为例，叠加运算与相关运算分别由插在计算机总线上的噪声编辑处理板卡(即，NEP板)和快速傅立叶处理板卡(即，FTP板)来完成。其中FTP板拥有2片摩托罗拉96002DSP处理器，利用快速傅里叶变换算法(FFT)在频域执行相关运算。这两块处理板配合使用在两毫秒采样率二十六秒采集长度的条件下，最大处理能力为1800道。

使用硬件进行相关/叠加运算的优点是非常明显的，那就是实时性好、反应迅速。例如：摩托罗拉96002DSP处理器为精简指令集，具有每秒一千三百五十万次的处理能力。

但使用硬件进行相关/叠加运算也存在一些难以避免的问题：

(1)随着地震勘探施工方法的不断发展与丰富，使得同一时刻参与数据采集的地震道数不断扩大，同一时刻需要进行相关/叠加运算的数据海量增加，使得使用硬件相关/叠加运算器的数据吞吐能力出现了瓶颈。

(2)由于受到硬件设计目标局限性的影响，以及选用的电子器件技术指标的影响，这种“瓶颈”很难在后期使用阶段通过某种扩展手段进行根本性的解决。

(3)相关/叠加运算器的硬件电路设计复杂，研发周期长，难以在短时间保持研发与施工需求变化同步。

(4)某些行业进行地震采集施工时(如：煤田地质行业)对同一时刻进行的地震采集道数要求比较小，如果配置处理能力过强的相关/叠加运算器硬件设备，从资源成本的角度上讲，显得过于浪费。但为此专门开发一种相应的相关/叠加运算器的硬件设备又得不偿失。

(5)研发相关/叠加运算器硬件设备，资金投入大，风险高。

随着通用计算机CPU运算能力的不断提高，硬件相关/叠加运算器有逐步转为全软件方式实现。

使用软件进行相关运算常用的方法为利用傅立叶变换性质和快速傅立叶算法，对一组数字信号利用两次正变换和一次逆变换，进行运算求解。同时域算法相比，提高了计算效率。

由于在地震勘探数据采集过程中进行相关运算的信号均属于实序列，而傅立叶变换是复序列，在运算过程中存在很大的冗余，限制了存储效率与运算速度的进一步提高。利用实序列傅立叶变换的对称性进行一步简化运算方法，通过将两组实序列构造成复序列结构进行运算，可通过一次变换得到两组实序列的傅立叶变换。这比直接采用傅里叶变换减少了一次傅立叶变换的处理过程，从而节省了近三分之一的运算量，大大缩短了运算时间。

以国外某种大型地震仪器为例，使用SUN Blade 2500工作站(8G内存)在2毫秒采样率，26秒采集长度的条件下，最大处理能力为1800道。

使用软件进行相关/叠加运算的优点是：

(1)可扩展性好，用户可以根据施工处理需要，对运行相关/叠加运算的工作站配置指标进行适当的调整。

(2)相对于研发硬件相关/叠加运算器而言，开发软件相关/叠加运算器研发周期短。

(3)利用软件手段实现通用的相关/叠加算法较为简单，因此与研发硬件相关/叠加运算器相比风险小。

(4)研发所需资金投入较小。

现今使用的软件相关/叠加运算急需解决的问题：

(1)同硬件相关/叠加运算器相比，实时性差，处理同等数量级的数据时软件相关/叠加运算器需要更长的时间。

(2)难以适应地震勘探施工中大道数地震数据采集的发展趋势。

发明内容

本发明的目的是提供一种对FFT-Hartly混合算法的计算流程和运算步骤进行重新的编排与设计，增加相关运算器的运算效率，提高运算速度，减少运算时间的地震勘探相关/叠加数据处理方法。从地震采集施工数据流回收工作流程出发，结合多核、多CPU工作站(服务器)运算与并行处理机制，设计、编排相关/叠加运算的工作处理顺序。在最大限度上复用数据流回收工作流程中，步骤间工作状态相互转换的空白时间，予以实现相关/叠加运算器拥有快速处理大道数地震采集施工中地震采集数据的处理能力。

相关/叠加运算器处理能力的预期目标：

(1)在同一时刻进行8000道数据采集，地震采集施工技术指标为，2毫秒采样率，可控震源扫描长度为24秒，进行32秒地震数据采集，相关/叠加运算后输出8秒记录长度时，相关/叠加运算的处理时间应在8秒内完成。

(2)在同一时刻进行4000道数据采集，地震采集施工技术指标为，2毫秒采样率，可控震源扫描长度为24秒，进行32秒地震数据采集，相关/叠加运算后输出8秒记录长度时，相关/叠加运算的处理时间应在4秒内完成。

本发明所述的地震勘探相关/叠加数据处理方法是由

(1)快速傅里叶变换算法与哈特利变换性质的结合(FFT-Hartly混合算法)

地震勘探数据采集中的相关运算在从时域到频域，频域到时域的变换过程中展开，正变换后对频域中的数据进行如下构造：

$R_{\mathrm{xy}}^h\left(k\right)=X_R\left(k\right)[Y_R\left(k\right)+Y_I\left(k\right)]-X_I\left(k\right)[Y_R\left(k\right)-Y_I\left(k\right)]$

$k=0,......,\frac{N}{2}$

$R_{\mathrm{xy}}^h(N-k)=X_R\left(k\right)[Y_R\left(k\right)-Y_I\left(k\right)]+X_I\left(k\right)[Y_R\left(k\right)+Y_I\left(k\right)]$

其反变换过程则有：

$r_{\mathrm{xy}}\left(n\right)=\mathrm{IDHT}\left[R_{\mathrm{xy}}^h\left(k\right)\right]$

$=\frac{1}{N}\mathrm{DHT}\left[R_{\mathrm{xy}}^h\left(k\right)\right]$

$=\frac{1}{N}\{\mathrm{Re}\left(\mathrm{DFT}\right[R_{\mathrm{xy}}^h\left(k\right)\left]\right)-\mathrm{Im}\left(\mathrm{DFT}\right[R_{\mathrm{xy}}^h\left(k\right)\left]\right)\}$

以FFT正变换的形式实现了从频域向时域的反变换过程。

(2)FFT-Hartly混合算法中的分段与求和

将一组实序列分成几段，分别进行FFT运算并保存中间结果

下标j为正整数，代表第几段，当这组实序列所有分段全部运算结束后，通过将保存的中间结果进行累加求出

再进行后续的运算；

设需要输出的相关r_xy(n)样点长度为L，扫描记录y(n)样点长度为M，采集记录x(n)样点长度为N，将y(n)分为每段长为p的J段，则有M＝J×p；对x(n)以N₁为长度分段，且满足N₁＝2^k1≥L+p，k1代表某一满足条件的正整数，当N₁＝2^k1＞L+p时，要通过补零的方法使得补零后的L+p段与N₁相等才能进行运算；如果x(n)某一分段，通常是最后一段，由于样点数的不足，没有达到的要求，也应该通过补零的手段填满剩余的部分再进行运算；

$R_{\mathrm{xy}}^h\left(k\right)=\underset{j=0}{\overset{J-1}{\mathrm{\Σ}}}{R}_{x_j{y}_j}^h\left(k\right)$ 0≤k≤N-1         0≤j≤J-1

求出可以根据FFT-Hartly混合算法中讨论过的结论求出r_xy(n)；

(3)多核、多CPU并行处理的线程安排

将一个地震道采集数据多个分段中的一个数据分段作为一个单独的处理线程，以多核CPU中的一个物理内核作为一个独立处理单元，由这个处理单元完成对这个数据分段的运算与处理，并以中间结果的形式存入临时缓冲区，(即，

随后会有其它数据分段处理线程分配到这个处理单进行运算，周而复始直至所有数据分段的处理线程全部计算完毕，再到缓冲区中取出中间结果，以数据分段的应对关系进行相加求和(即，

这样每一个地震道采集数据都会求得一个

然后再将每一个地震道采集数据的作为一个处理线程，分配到多核CPU中的一个物理内核中进行最后的运算，求出最后的计算结果(即，r_xy(n))。

这样软件相关/叠加运算器在多核、多CPU的条件下，可以同一时刻完成对多个地震道采集数据进行并行计算与处理。

综上所述，通过对算法结构的有机调整、计算过程中最大程度的时间复用以及大量采用多线程技术的软件工程手段，使得相关/叠加运算器的计算效率、处理能力得到了巨大的提高。

发明的效果

通过在IBM 3755服务器(4核)上进行测试可以得到如下指标：在不加载噪声控制情况下，2毫秒采样率，24秒扫描信号，32秒地震采集，相关/叠加运算后得到8秒记录长度的情况下，处理4000道时间不大于3秒，处理8000道不大于4秒水平。对于加载噪声控制部分时(分段加权叠加)运算8000道(每道记录长度与前面相同)不大于7秒。

附图说明

图1FFT-Hartly分段工作处理流程图

图2分段逆功率加权叠加算法及处理流程图

具体实施方式

以下以分为三段相关处理为例：对单个线程而言，对每一道的数据先进行分段，然后对最后一段不够部分进行补零；再对其各段数据进行傅立叶变换，变换后的各段数据与事先已经分段傅立叶变换后的相应的各段数据进行Hartly点乘；将Hartly点乘后的各段数据相累加；再进行傅立叶变换，得到相关处理结果。

分段逆功率加权叠加算法的计算过程是将地震道信号分成若干个窗口，分别计算各个窗口的功率，先按窗口功率的反比例对该窗口信号加权，然后逐段叠加，最后再进行归一化处理。这样可有效地降低受噪声污染的信号的影响，提高地震记录的信噪比。

本发明所述的地震勘探相关/叠加数据处理方法是由

(1)快速傅里叶变换算法与哈特利变换性质的结合(FFT-Hartly混合算法)哈特利变换DHT在对实序列进行运算的特性适合于从时域到频域，频域到时域的变换过程。地震勘探数据采集中的相关运算，均在两实序列中展开，因此特别哈特利变换特性非常符合在地震勘探数据采集过程中的应用。

哈特利变换DHT也有快速算法FHT，但FHT的运算效率远远低于快速傅里叶变换变换算法FFT。FFT算法是复序列运算，所以我们可以将FHT运算中的偶部与奇部，以FFT算法中的实部与虚部的形式进行构造。运用FFT算法实现FHT的运算，提高其计算效率。

设x(n)为长度为N的有限长实序列，离散哈特利变换DHT则有：

$X^h\left(k\right)=\mathrm{DHT}\left[x\left(n\right)\right]=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}x\left(n\right)\mathrm{CAS}\frac{2\mathrm{\πnk}}{N}$ (k＝0，......，N-1)

注： $\mathrm{CAS}\frac{2\mathrm{\πnk}}{N}=\cos \frac{2\mathrm{\πnk}}{N}+\sin \frac{2\mathrm{\πnk}}{N}$

由于x(n)x(n)为实序列，所以X^h(k)仍然为实数，所以离散哈特利变换DHT由性质可知X^h(-k)＝X^h(N-k)。在以上情况下离散哈特利变换的逆变换IDHT，除了常数因子N以外，正反变换的形式完全相同，则有：

x(n)＝IDHT[X^h(k)]

$=\frac{1}{N}\mathrm{DHT}\left[X^h\left(k\right)\right]$

$=\frac{1}{N}\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{X}^h\left(k\right)\mathrm{CAS}\frac{2\mathrm{\πnk}}{N}$ n＝0，......，N-1

如果设X(k)是x(n)的离散傅立叶变换，X_R(k)和X_I(k)分别代表X(k)的实部与虚部，X^h(k)是x(n)的离散哈特利变换，

和

分别代表X^h(k)的偶部与奇部，则X(k)和X^h(k)有下列相互表示关系：

X^h(k)＝X_R(k)-X_I(k)

$X_R\left(k\right)=X_e^h\left(k\right)=\frac{1}{2}[X^h\left(k\right)+X^h(N-k)]$

$X_I\left(k\right)=-X_o^h\left(k\right)=\frac{1}{2}[X^h(N-k)-X^{h}k)]$

根据DHT循环相关定理，设是x(n)和y(n)的循环相关，则有：

n＝0，......，N-1

令

为

的DHT，在满足DHT的长度N为2的正数幂，且

N≥M+L(M为相关运算中移动序列的长度，L代表相关运算后输出序列的长度)则有：

$R_{\mathrm{xy}}^h\left(k\right)=X^h\left(k\right)Y_e^h\left(k\right)-X^h(N-k)Y_o^h\left(k\right)$

根据上述情况的分解与推理，FFT-Hartly混合算法则可以这样描述，正变换后对频域中的数据进行如下构造：

$R_{\mathrm{xy}}^h\left(k\right)=X_R\left(k\right)[Y_R\left(k\right)+Y_I\left(k\right)]-X_I\left(k\right)[Y_R\left(k\right)-Y_I\left(k\right)]$

$k=0,......,\frac{N}{2}$

$R_{\mathrm{xy}}^h(N-k)=X_R\left(k\right)[Y_R\left(k\right)-Y_I\left(k\right)]+X_I\left(k\right)[Y_R\left(k\right)+Y_I\left(k\right)]$

其反变换过程则有：

$r_{\mathrm{xy}}\left(n\right)=\mathrm{IDHT}\left[R_{\mathrm{xy}}^h\left(k\right)\right]$

$=\frac{1}{N}\mathrm{DHT}\left[R_{\mathrm{xy}}^h\left(k\right)\right]$

$=\frac{1}{N}\{\mathrm{Re}\left(\mathrm{DFT}\right[R_{\mathrm{xy}}^h\left(k\right)\left]\right)-\mathrm{Im}\left(\mathrm{DFT}\right[R_{\mathrm{xy}}^h\left(k\right)\left]\right)\}$

综上所述，FFT-Hartly混合算法将实序列从时域转换到频域后，进行一系列的重新构建，以FFT正变换的形式实现了从频域向时域的反变换过程。通过这种设计，即充分利用了哈特利变换的特性，又发挥了快速傅里叶变换算法的计算高效性。由于整个算法中实际的计算过程均由快速傅里叶算法FFT完成，因此使得软件相关/叠加运算器在算法实现的结构与构成上，变得简洁明了，容易设计。

(2)FFT-Hartly混合算法中的分段与求和

由于FFT-Hartly混合算法在进行相关运算时，其主要的计算过程均由快速傅里叶算法FFT完成，因此在某种程度上可以将FFT-Hartly混合算法的运算过程视为“FFT”直接运算。

FFT运算的长度必须满足N为2的正数幂，且N≥M+L(见上节文字的介绍)。当需要FFT运算的实序列相对较长，为了满足FFT运算长度的条件(通过补零等手段)，可能会使得FFT的运算程度变得非常的长，从而大大的增加FFT的运算时间，这极不利于数据实时处理的要求。

解决这一问题的方法是将一组实序列分成几段，分别进行FFT运算并保存中间结果

(下标j为正整数，代表第几段)，当这组实序列所有分段全部运算结束后，通过将保存的中间结果进行累加求出

再进行后续的运算。

设需要输出的相关r_xy(n)样点长度为L，扫描记录y(n)样点长度为M，采集记录x(n)样点长度为N。将y(n)分为每段长为p的J段，则有M＝J×p；对x(n)以N₁为长度分段，且满足N₁＝2^k1≥L+p(k1代表某一满足条件的正整数)。当N₁＝2^k1＞L+p时，要通过补零的方法使得补零后的L+p段与N₁相等才能进行运算；如果x(n)某一分段(通常是最后一段)由于样点数的不足，没有达到的要求，也应该通过补零的手段填满剩余的部分再进行运算；

设x_j(n)、y_j(n)分别为满足上述条件的x(n)、y(n)样点个数达到N₁要求的分段，r_xjyj(n)则有：

x_j(n)＝x(n+jp)       0≤n≤N₁-1       j＝0，...J-1

注：若序列x(n)最后一段的长度不够，可在该分段的后端补零，以达到长度要求。

$y_j\left(n\right)=\left\{\begin{array}{cc}y(n+\mathrm{jp})& 0\≤n\≤p-1\\ 0& p\≤n\≤N_1\end{array}\right.,j=0,...,J-1$

$r_{x_j{y}_j}\left(n\right)=\underset{l=0}{\overset{p-1}{\mathrm{\Σ}}}{x}_j(n+l)y_j\left(l\right)$

$=\underset{l=0}{\overset{p-1}{\mathrm{\Σ}}}x(n+l+\mathrm{jp})y(l+\mathrm{jp})$

令k＝l+jp，且0≤n≤l时有：

$\underset{k=\mathrm{jp}}{\overset{(j+1)p-1}{\mathrm{\Σ}}}x(n+k)y\left(k\right)$

$\underset{j=0}{\overset{J-1}{\mathrm{\Σ}}}{r}_{x_j{y}_j}\left(n\right)=\underset{j=0}{\overset{J-1}{\mathrm{\Σ}}}\left(\underset{k=\mathrm{jp}}{\overset{(j+1)p-1}{\mathrm{\Σ}}}x(n+k)y\left(k\right)\right)$

$=\underset{k=0}{\overset{\mathrm{Jp}-1}{\mathrm{\Σ}}}x(n+k)y\left(k\right)$

$=\underset{k=0}{\overset{M-1}{\mathrm{\Σ}}}x(n+k)y\left(k\right)$

$=r_{\mathrm{xy}}\left(n\right)$

所以可以得出：

$R_{\mathrm{xy}}^h\left(k\right)=\underset{j=0}{\overset{J-1}{\mathrm{\Σ}}}{R}_{x_j{y}_j}^h\left(k\right)$ 0≤k≤N-1     0≤j≤J-1

求出

可以根据FFT-Hartly混合算法中讨论过的结论求出r_xy(n)。

分段求和法的应用，使得利用软件相关/叠加运算器实时处理长实序列成为了可能。这样的设计为软件相关/叠加运算器带来了两大好处：一是利用算法结构上的变换，从根本上解决了快速傅里叶变换FFT算法对长序列进行运算时，由于算法本身特性的原因造成的计算耗时过程的问题；二是由于计算是分段进行的，所以软件相关/叠加运算器不需要等到地震采集数据全部收集完整再进行计算，只需要一个分段中的数据收集完整就可以开始计算，从而充分复用了地震数据采集回收的等待时间，在人的主观上大大缩短了软件相关/叠加运算器的计算处理时间。

(3)多核、多CPU并行处理的线程安排

将一个地震道采集数据多个分段中的一个数据分段作为一个单独的处理线程，以多核CPU中的一个物理内核作为一个独立处理单元，由这个处理单元完成对这个数据分段的运算与处理，并以中间结果的形式存入临时缓冲区(即，

随后会有其它数据分段处理线程分配到这个处理单进行运算，周而复始直至所有数据分段的处理线程全部计算完毕，再到缓冲区中取出中间结果，以数据分段的应对关系进行相加求和(即，

这样每一个地震道采集数据都会求得一个

然后再将每一个地震道采集数据的

作为一个处理线程，分配到多核CPU中的一个物理内核中进行最后的运算，求出最后的计算结果(即，r_xy(n))。

这样软件相关/叠加运算器在多核、多CPU的条件下，可以同一时刻完成对多个地震道采集数据进行并行计算与处理。

综上所述，通过对算法结构的有机调整、计算过程中最大程度的时间复用以及大量采用多线程技术的软件工程手段，使得相关/叠加运算器的计算效率、处理能力得到了巨大的提高。

Claims (1)

1.一种地震勘探相关/叠加数据处理方法，其特征在于：是由

(1)快速傅里叶变换算法与哈特利变换性质的结合混合算法：

地震勘探数据采集中的相关运算在从时域到频域，频域到时域的变换过程中展开，正变换后对频域中的数据进行如下构造：

$R_{\mathrm{xy}}^h\left(k\right)=X_R\left(k\right)[Y_R\left(k\right)+Y_I\left(k\right)]-X_I\left(k\right)[Y_R\left(k\right)-Y_I\left(k\right)]$

$k=0,......,\frac{N}{2}$

$R_{\mathrm{xy}}^h(N-k)=X_R\left(k\right)[Y_R\left(k\right)-Y_I\left(k\right)]+X_I\left(k\right)[Y_R\left(k\right)+Y_I\left(k\right)]$

其反变换过程则有：

$r_{\mathrm{xy}}\left(n\right)=\mathrm{IDHT}\left[R_{\mathrm{xy}}^h\left(k\right)\right]$

$=\frac{1}{N}\mathrm{DHT}\left[R_{\mathrm{xy}}^h\left(k\right)\right]$

$=\frac{1}{N}\{\mathrm{Re}\left(\mathrm{DFT}\right[R_{\mathrm{xy}}^h\left(k\right)\left]\right)-\mathrm{Im}\left(\mathrm{DFT}\right[R_{\mathrm{xy}}^h\left(k\right)\left]\right)\}$

以FFT正变换的形式实现了从频域向时域的反变换过程；

(2)FFT-Hartly混合算法中分段与求和

将一组实序列分成几段，分别进行FFT运算并保存中间结果

下标j为正整数，代表第几段，当这组实序列所有分段全部运算结束后，通过将保存的中间结果进行累加求出

再进行后续的运算；

设需要输出的相关r_xy(n)样点长度为L，扫描记录y(n)样点长度为M，采集记录x(n)样点长度为N，将y(n)分为每段长为p的J段，则有M＝J×p；对x(n)以N₁为长度分段，且满足N₁＝2^k1≥L+p，k1代表某一满足条件的正整数，当N₁＝2^k1＞L+p时，要通过补零的方法使得补零后的L+p段与N₁相等才能进行运算；如果x(n)某一分段，通常是最后一段，由于样点数的不足，没有达到的要求，也应该通过补零的手段填满剩余的部分再进行运算；

$R_{\mathrm{xy}}^h\left(k\right)=\underset{j=0}{\overset{J-1}{\mathrm{\Σ}}}{R}_{x_j{y}_j}^h\left(k\right)$ 0≤k≤N-1    0≤j≤J-1

求出

可以根据FFT-Hartly混合算法中讨论过的结论求出r_xy(n)；

(3)多核、多CPU并行处理的线程安排

将一个地震道采集数据多个分段中的一个数据分段作为一个单独的处理线程，以多核CPU中的一个物理内核作为一个独立处理单元，由这个处理单元完成对这个数据分段的运算与处理，并以中间结果的形式存入临时缓冲区，

随后会有其它数据分段处理线程分配到这个处理单进行运算，周而复始直至所有数据分段的处理线程全部计算完毕，再到缓冲区中取出中间结果，以数据分段的应对关系进行相加求和，

这样每一个地震道采集数据都会求得一个

然后再将每一个地震道采集数据的作为一个处理线程，分配到多核CPU中的一个物理内核中进行最后的运算，求出最后的计算结果r_xy(n)。

Images