专利名称：一种地震勘探相关/叠加数据处理方法
技术领域：
本发明涉及一种大型地震仪器野外地震勘探数据采集快速相关/叠加数据处理 方法。
背景技术：
利用可控震源作为地震采集信号激发源，进行野外地震勘探数据采集施工已经非 常普遍。对采集的地震数据进行相关/叠加运算与处理的最常使用的方法是通过硬件的手 段制作相关运算器与叠加运算器，进行实时的运算处理。以国外某种大型地震仪器为例，叠 加运算与相关运算分别由插在计算机总线上的噪声编辑处理板卡(即，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混合算法)地震勘探数据采集中的相关运算在从时域到频域，频域到时域的变换过程中展 开，正变换后对频域中的数据进行如下构造 其反变换过程则有 以FFT正变换的形式实现了从频域向时域的反变换过程。(2) FFT-Hartly混合算法中的分段与求和将一组实序列分成几段，分别进行FFT运算并保存中间结果及kx.y. 0 )，下 标j为正整数，代表第几段，当这组实序列所有分段全部运算结束后，通过将保存的中间结 果进行累加求出及hxy Ck )再进行后续的运算；设需要输出的相关rxy(n)样点长度为L，扫描记录y (n)样点长度为M，采集记录 x (n)样点长度为N，将y (n)分为每段长为p的J段，则有M = JXp ；对x (n)以R为长度分 段，且满足K = 2kl彡L+p，kl代表某一满足条件的正整数，当& = 2kl > L+p时，要通过补 零的方法使得补零后的L+p段与&相等才能进行运算；如果x(n)某一分段，通常是最后一 段，由于样点数的不足，没有达到的要求，也应该通过补零的手段填满剩余的部分再进行运 算； 求出i hxy {k )可以根据FFT-Hartly混合算法中讨论过的结论求出rxy(n)；(3)多核、多CPU并行处理的线程安排将一个地震道采集数据多个分段中的一个数据分段作为一个单独的处理线程，以 多核CPU中的一个物理内核作为一个独立处理单元，由这个处理单元完成对这个数据分段
的运算与处理，并以中间结果的形式存入临时缓冲区，(即，灭^。(k))。随后会有其它数据分段处理线程分配到这个处理单进行运算，周而复始直至所有 数据分段的处理线程全部计算完毕，再到缓冲区中取出中间结果，以数据分段的应对关系
进行相加求和 这样每一个地震道采集数据都会求得一个及Hxy (k )，然后再将每一个地震道采
集数据的及t (A;)作为一个处理线程，分配到多核cpu中的一个物理内核中进行最后 的运算，求出最后的计算结果(lf，rxy(n))。这样软件相关/叠加运算器在多核、多CPU的条件下，可以同一时刻完成对多个地 震道采集数据进行并行计算与处理。综上所述，通过对算法结构的有机调整、计算过程中最大程度的时间复用以及大 量采用多线程技术的软件工程手段，使得相关/叠加运算器的计算效率、处理能力得到了 巨大的提高。
发明的效果通过在IBM 3755服务器(4核)上进行测试可以得到如下指标在不加载噪声控 制情况下，2毫秒采样率，24秒扫描信号，32秒地震采集，相关/叠加运算后得到8秒记录长 度的情况下，处理4000道时间不大于3秒，处理8000道不大于4秒水平。对于加载噪声控 制部分时(分段加权叠加)运算8000道(每道记录长度与前面相同)不大于7秒。
图IFFT-Hartly分段工作处理流程2分段逆功率加权叠加算法及处理流程图
具体实施例方式以下以分为三段相关处理为例对单个线程而言，对每一道的数据先进行分段，然 后对最后一段不够部分进行补零；再对其各段数据进行傅立叶变换，变换后的各段数据与 事先已经分段傅立叶变换后的相应的各段数据进行Hartly点乘；将Hartly点乘后的各段 数据相累加；再进行傅立叶变换，得到相关处理结果。分段逆功率加权叠加算法的计算过程是将地震道信号分成若干个窗口，分别计算 各个窗口的功率，先按窗口功率的反比例对该窗口信号加权，然后逐段叠加，最后再进行归 一化处理。这样可有效地降低受噪声污染的信号的影响，提高地震记录的信噪比。本发明所述的地震勘探相关/叠加数据处理方法是由(1)快速傅里叶变换算法与哈特利变换性质的结合(FFT-Hartly混合算法)哈特 利变换DHT在对实序列进行运算的特性适合于从时域到频域，频域到时域的变换过程。地 震勘探数据采集中的相关运算，均在两实序列中展开，因此特别哈特利变换特性非常符合 在地震勘探数据采集过程中的应用。哈特利变换DHT也有快速算法FHT，但FHT的运算效率远远低于快速傅里叶变换变 换算法FFT。FFT算法是复序列运算，所以我们可以将FHT运算中的偶部与奇部，以FFT算 法中的实部与虚部的形式进行构造。运用FFT算法实现FHT的运算，提高其计算效率。设x(n)为长度为N的有限长实序列，离散哈特利变换DHT则有
N-iInnk 由于X(n)X(n)为实序列，所以Xh(k)仍然为实数，所以离散哈特利变换DHT由性质 可知Xh(-k) =Xh(N-k)。在以上情况下离散哈特利变换的逆变换IDHT，除了常数因子N以 外，正反变换的形式完全相同，则有x(n) = IDHT[Xh(k)] 如果设X(k)是x(n)的离散傅立叶变换，XK(k)和Xjk)分别代表X(k)的实部与虚部，Xh(k)是x(n)的离散哈特利变换，和^C(A)分别代表Xh(k)的偶部与奇部，则X(k)
和Xh(k)有下列相互表示关系Xh(k) =XE(k)-XI(k)XR (/t) = Xhe (k) = (k) + Xh (N-；t)]Xi (k) = -Xho (k) =(N-k)-Xh ⑷]根据DHT循环相关定理，设冷(n)是x (n)和y (n)的循环相关，则有 0/ N_l 令为冷(…的DHT,在满足DHT的长度N为2的正数幂，且N>M+L(M为相关运算中移动序列的长度，L代表相关运算后输出序列的长度)则 有Rhxy (k) = Xh (k)rf (k)-Xh{N-k、Yoh {k)根据上述情况的分解与推理，FFT-Hartly混合算法则可以这样描述，正变换后对 频域中的数据进行如下构造 其反变换过程则有rxy{n) = IDHT[Rhxy{k)\^^DHT[Rhxy{k)}=去⑷])-lm{DFT[Rhxy㈨])}综上所述，FFT-Hartly混合算法将实序列从时域转换到频域后，进行一系列的重 新构建，以FFT正变换的形式实现了从频域向时域的反变换过程。通过这种设计，即充分利 用了哈特利变换的特性，又发挥了快速傅里叶变换算法的计算高效性。由于整个算法中实 际的计算过程均由快速傅里叶算法FFT完成，因此使得软件相关/叠加运算器在算法实现 的结构与构成上，变得简洁明了，容易设计。(2) FFT-Hartly混合算法中的分段与求和由于FFT-Hartly混合算法在进行相关运算时，其主要的计算过程均由快速傅里 叶算法FFT完成，因此在某种程度上可以将FFT-Hartly混合算法的运算过程视为“FFT”直 接运算。FFT运算的长度必须满足N为2的正数幂，且N彡M+L(见上节文字的介绍)。当 需要FFT运算的实序列相对较长，为了满足FFT运算长度的条件(通过补零等手段)，可能 会使得FFT的运算程度变得非常的长，从而大大的增加FFT的运算时间，这极不利于数据实 时处理的要求。
8
解决这一问题的方法是将一组实序列分成几段，分别进行FFT运算并保存中间结 果凡(下标j为正整数，代表第几段)，当这组实序列所有分段全部运算结束后，通过 将保存的中间结果进行累加求出再进行后续的运算。设需要输出的相关rxy(n)样点长度为L，扫描记录y (n)样点长度为M，采集记录 x(n)样点长度为N。将y(n)分为每段长为p的J段，则有M = JXp;对x(n)以R为长度 分段，且满足K = 2kl彡L+p (kl代表某一满足条件的正整数)。当& = 2kl > L+p时，要通 过补零的方法使得补零后的L+p段与&相等才能进行运算；如果x(n)某一分段(通常是 最后一段)由于样点数的不足，没有达到的要求，也应该通过补零的手段填满剩余的部分 再进行运算；设Xj (n)、Yj (n)分别为满足上述条件的x (n)、y (n)样点个数达到R要求的分段， rXJyJ(n)则有 x,(n) = x (n+jp)
0彡n彡N,-l
j = 0， J-i 注若序列X(n)最后一段的长度不够，可在该分段的后端补零，以达到长度要求.
令k = l+jp，且0彡n彡l时有 所以可以得出 求出i hxy (k )可以根据FFT-Hartly混合算法中讨论过的结论求出rxy(n)。分段求和法的应用，使得利用软件相关/叠加运算器实时处理长实序列成为了可 能。这样的设计为软件相关/叠加运算器带来了两大好处一是利用算法结构上的变换，从 根本上解决了快速傅里叶变换FFT算法对长序列进行运算时，由于算法本身特性的原因造
9成的计算耗时过程的问题；二是由于计算是分段进行的，所以软件相关/叠加运算器不需 要等到地震采集数据全部收集完整再进行计算，只需要一个分段中的数据收集完整就可以 开始计算，从而充分复用了地震数据采集回收的等待时间，在人的主观上大大缩短了软件 相关/叠加运算器的计算处理时间。(3)多核、多CPU并行处理的线程安排将一个地震道采集数据多个分段中的一个数据分段作为一个单独的处理线程，以 多核CPU中的一个物理内核作为一个独立处理单元，由这个处理单元完成对这个数据分段
的运算与处理，并以中间结果的形式存入临时缓冲区(即，及(k))。随后会有其它数据分段处理线程分配到这个处理单进行运算，周而复始直至所有 数据分段的处理线程全部计算完毕，再到缓冲区中取出中间结果，以数据分段的应对关系
进行相加求和(
)。这样每一个地震道采集数据都会求得一个i O )，然后再将每一个地震道采 集数据的^0 )作为一个处理线程，分配到多核CPU中的一个物理内核中进行最后的运 算，求出最后的计算结果(即，rxy(n))。这样软件相关/叠加运算器在多核、多CPU的条件下，可以同一时刻完成对多个地 震道采集数据进行并行计算与处理。综上所述，通过对算法结构的有机调整、计算过程中最大程度的时间复用以及大 量采用多线程技术的软件工程手段，使得相关/叠加运算器的计算效率、处理能力得到了 巨大的提高。
权利要求
一种地震勘探相关/叠加数据处理方法，其特征在于是由(1)快速傅里叶变换算法与哈特利变换性质的结合混合算法地震勘探数据采集中的相关运算在从时域到频域，频域到时域的变换过程中展开，正变换后对频域中的数据进行如下构造 <mrow><msubsup> <mi>R</mi> <mi>xy</mi> <mi>h</mi></msubsup><mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo></mrow><mo>=</mo><msub> <mi>X</mi> <mi>R</mi></msub><mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo></mrow><mo>[</mo><msub> <mi>Y</mi> <mi>R</mi></msub><mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo></mrow><mo>+</mo><msub> <mi>Y</mi> <mi>I</mi></msub><mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo></mrow><mo>]</mo><mo>-</mo><msub> <mi>X</mi> <mi>I</mi></msub><mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo></mrow><mo>[</mo><msub> <mi>Y</mi> <mi>R</mi></msub><mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo></mrow><mo>-</mo><msub> <mi>Y</mi> <mi>I</mi></msub><mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo></mrow><mo>]</mo> </mrow> <mrow><mi>k</mi><mo>=</mo><mn>0</mn><mo>,</mo><mo>.</mo><mo>.</mo><mo>.</mo><mo>.</mo><mo>.</mo><mo>.</mo><mo>,</mo><mfrac> <mi>N</mi> <mn>2</mn></mfrac> </mrow> <mrow><msubsup> <mi>R</mi> <mi>xy</mi> <mi>h</mi></msubsup><mrow> <mo>(</mo> <mi>N</mi> <mo>-</mo> <mi>k</mi> <mo>)</mo></mrow><mo>=</mo><msub> <mi>X</mi> <mi>R</mi></msub><mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo></mrow><mo>[</mo><msub> <mi>Y</mi> <mi>R</mi></msub><mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo></mrow><mo>-</mo><msub> <mi>Y</mi> <mi>I</mi></msub><mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo></mrow><mo>]</mo><mo>+</mo><msub> <mi>X</mi> <mi>I</mi></msub><mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo></mrow><mo>[</mo><msub> <mi>Y</mi> <mi>R</mi></msub><mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo></mrow><mo>+</mo><msub> <mi>Y</mi> <mi>I</mi></msub><mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo></mrow><mo>]</mo> </mrow>其反变换过程则有 <mrow><msub> <mi>r</mi> <mi>xy</mi></msub><mrow> <mo>(</mo> <mi>n</mi> <mo>)</mo></mrow><mo>=</mo><mi>IDHT</mi><mo>[</mo><msubsup> <mi>R</mi> <mi>xy</mi> <mi>h</mi></msubsup><mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo></mrow><mo>]</mo> </mrow> <mrow><mo>=</mo><mfrac> <mn>1</mn> <mi>N</mi></mfrac><mi>DHT</mi><mo>[</mo><msubsup> <mi>R</mi> <mi>xy</mi> <mi>h</mi></msubsup><mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo></mrow><mo>]</mo> </mrow> <mrow><mo>=</mo><mfrac> <mn>1</mn> <mi>N</mi></mfrac><mo>{</mo><mi>Re</mi><mrow> <mo>(</mo> <mi>DFT</mi> <mo>[</mo> <msubsup><mi>R</mi><mi>xy</mi><mi>h</mi> </msubsup> <mrow><mo>(</mo><mi>k</mi><mo>)</mo> </mrow> <mo>]</mo> <mo>)</mo></mrow><mo>-</mo><mi>Im</mi><mrow> <mo>(</mo> <mi>DFT</mi> <mo>[</mo> <msubsup><mi>R</mi><mi>xy</mi><mi>h</mi> </msubsup> <mrow><mo>(</mo><mi>k</mi><mo>)</mo> </mrow> <mo>]</mo> <mo>)</mo></mrow><mo>}</mo> </mrow>以FFT正变换的形式实现了从频域向时域的反变换过程；(2)FFT Hartly混合算法中分段与求和将一组实序列分成几段，分别进行FFT运算并保存中间结果下标j为正整数，代表第几段，当这组实序列所有分段全部运算结束后，通过将保存的中间结果进行累加求出再进行后续的运算；设需要输出的相关rxy(n)样点长度为L，扫描记录y(n)样点长度为M，采集记录x(n)样点长度为N，将y(n)分为每段长为p的J段，则有M＝J×p；对x(n)以N1为长度分段，且满足N1＝2k1≥L+p，k1代表某一满足条件的正整数，当N1＝2k1＞L+p时，要通过补零的方法使得补零后的L+p段与N1相等才能进行运算；如果x(n)某一分段，通常是最后一段，由于样点数的不足，没有达到的要求，也应该通过补零的手段填满剩余的部分再进行运算； <mrow><msubsup> <mi>R</mi> <mi>xy</mi> <mi>h</mi></msubsup><mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo></mrow><mo>=</mo><munderover> <mi>Σ</mi> <mrow><mi>j</mi><mo>=</mo><mn>0</mn> </mrow> <mrow><mi>J</mi><mo>-</mo><mn>1</mn> </mrow></munderover><msubsup> <mi>R</mi> <mrow><msub> <mi>x</mi> <mi>j</mi></msub><msub> <mi>y</mi> <mi>j</mi></msub> </mrow> <mi>h</mi></msubsup><mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo></mrow> </mrow>0≤k≤N 10≤j≤J 1求出可以根据FFT Hartly混合算法中讨论过的结论求出rxy(n)；(3)多核、多CPU并行处理的线程安排将一个地震道采集数据多个分段中的一个数据分段作为一个单独的处理线程，以多核CPU中的一个物理内核作为一个独立处理单元，由这个处理单元完成对这个数据分段的运算与处理，并以中间结果的形式存入临时缓冲区，随后会有其它数据分段处理线程分配到这个处理单进行运算，周而复始直至所有数据分段的处理线程全部计算完毕，再到缓冲区中取出中间结果，以数据分段的应对关系进行相加求和，这样每一个地震道采集数据都会求得一个然后再将每一个地震道采集数据的作为一个处理线程，分配到多核CPU中的一个物理内核中进行最后的运算，求出最后的计算结果rxy(n)。F2009100867487C0000017.tif,F2009100867487C0000021.tif,F2009100867487C0000023.tif,F2009100867487C0000024.tif,F2009100867487C0000025.tif,F2009100867487C0000026.tif,F2009100867487C0000031.tif
全文摘要
本发明涉及一种地震勘探相关/叠加数据处理方法，(1)快速傅里叶变换算法与哈特利变换性质的结合混合算法；(2)FFT-Hartly混合算法中分段与求和；(3)多核、多CPU并行处理的线程安排；通过在IBM 3755服务器(4核)上进行测试可以得到如下指标在不加载噪声控制情况下，2毫秒采样率，24秒扫描信号，32秒地震采集，相关/叠加运算后得到8秒记录长度的情况下，处理4000道时间不大于3秒，处理8000道不大于4秒水平，对于加载噪声控制部分时(分段加权叠加)运算8000道(每道记录长度与前面相同)不大于7秒。
文档编号G01V1/28GK101930079SQ200910086748
公开日2010年12月29日 申请日期2009年6月26日 优先权日2009年6月26日
发明者刘益成, 张洁, 张红震, 李永权, 沈孝科, 王浩, 穆群英, 赵培根, 魏启 申请人:西安英诺瓦物探装备有限公司