专利名称：地震数据的压缩方法和装置的制作方法
技术领域：
本发明涉及地震数据的压缩方法和装置。
数据压缩(或缩减)是为了降低被处理的数据量而在处理中不损失主要信息的数字信号处理技术。这主要是通过除去数据中的冗余而实现的，并可能涉及舍去数据的无关紧要的部分。这种压缩的结果可能造成数据精确性的某些损失。允许精确重构原始数据的数据压缩在文献中常称为无损的数据压缩。涉及到精确性的某种降低的数据压缩称为有损的数据压缩。“舍入法”及“下采样法”(down sampling)是数据压缩普通的例子；这两个方法通常是有损的。
地震数据的探测需要进行大量的地震实验以获得可靠的地表下图象。每一实验要使用适当的声源产生声波并由大量的接收器测量地球的响应。大规模的地震观测这样就产生了通常为数字格式的大量的数据，这些数据需要被传输，存储和处理。为了便于对这种大量的数据进行处理，可以采用数据压缩。
地震数据探测中通常采用的有损数据压缩技术是分组形成法。这涉及在固定大小的组内大量的相邻接收器的保持和传输及处理，而不是个别的测量。
对数据压缩主要不使用分组形成法。分组形成法抑制随机的环境噪声，并抑制诸如陆地地震中的地滚等低表观速度的波。这种分组形成法减弱了数据的高的空间频率成分。然而，这一减弱是以粗糙的方式进行的，即仅部分地抑制了表面上呈缓慢传播的波并改变了其余的数据。因而，有充分理由从探测阶段就省去分组形成法，并个别记录每一接收器的输出。这时允许采用更为精致的方法以降低随机及相干噪声。然而，在探测阶段抛弃分组形成法则大大增加了随后的各阶段要处理的数据量。
在IEEE Int．SYM．Circuits & Systym，New Orleans，LA，1-3 May1990，Vol．2，1573-6中，A．Spanias，S．Johnson等人阐述了几种用于地震数据压缩的基于变换的方法。这些方法包括离散傅立叶变换(DFT)，离散余弦变换(DCT)，Walsh-Hardamard变换(WHT)，及Karhunen-Loeve变换(KLT)。然而，在该刊物中所描述的并应用于N个数据点的滑动帧的DCT能够用于在几个不同的变换之间进行相对比较。当用作数据压缩方法时，滑动帧在变换域产生大量的冗余数据。
因而本发明的目的是提供用于压缩地震数据的一种方法。本发明的另一目的是要提供一种不使用分组形成法的压缩地震数据的方法。
本发明提供了其中把Ⅳ型的局部空间或时间的离散三角(即，正弦或余弦)变换用于地震数据信号的第一级压缩。Ⅳ型离散的正弦／余弦变换本身是公知的。例如由H．S．Malvar在“Lapped transforms for efficienttransform／subband coding”，IEEE ASSP，vol．38，no．6，June 1990中给出了一般的阐述。局部空间或时间离散余弦／正弦变换得到的结果是得到比原始数据更加紧凑和较少相关性的变换系数。在以下的数据处理步骤中最好采用这两个性质。
在一个处理步骤中采用了变换系数的紧凑性，这步骤可被描述为重新量化(requantization)或舍入步骤。这一步骤的目的是要以高精度保留所选择的系数，并以较低的精度保留其它系数，以便降低描述系数所需的数据量并从而实现进一步的数据压缩。
降低的变换系数的相关性提供了采用编码方案的机会，以至进一步降低要被存储或传输的数据量。可采用的编码方案本身是公知的例如Huffmann编码或幅度编码。
采用该方法的地震数据信号一般是从多个接收器获得的跟踪图形，例如从地下测音器或水下测音器获得的。这些接收器可组合安置，所有这些都是现有技术中已经知道的。例如，其中之一是排布在数个平行线中的地下测音器的线性阵列的这种传统的3维陆地地震学布局。方法中使用局部变换允许对包含在每一条线上的一定数目接收器的数据进行压缩。如本专业内所熟知，局部变换是这样一种变换，其中变换施加到所定义的跟踪图形的窗口。这样，在每一相继的变换阶段施加了局部变换的跟踪图形的数目称为空间窗口，而这一窗口可能随着所施加的变换类型而改变。
局部变换窗口由窗口的功能所定义，而窗口功能是这样选择的，使得变换是正交的并且是可逆的。窗口功能的选择使得变换施加到与邻接的窗口有重叠的的中心窗口，最好是重叠这些窗口的一半。
变换可以分两个步骤进行，第一步骤是一折叠步骤，其中的中心窗口与邻接的半个窗口组合而产生一折叠信号，且第二步骤是压缩对折叠信号所进行的余弦变换。
除了局部空间变换之外，最好还向数据施加局部时间变换。两种局部变换的组合可达到较好的压缩比。局部时间变换最好是Ⅳ型局部时间离散正弦／余弦变换。然而，其它的信号变换和分解也可使用，诸如普通的局部离散余弦变换，及局部傅立叶变换。可以按任何顺序施加局部空间三角变换和局部时间三角变换。
表示变换之后原始数据的变换系数形成了一个数据集，不同的压缩方法可施加到该数据集。这些压缩方法可被共同称为(重新)量化和编码。当用于压缩数据时量化过程通常包括定标步骤和舍入。量化过程是为了降低高频成分或系数，同时以较高精度保留低频成分。
定标最好通过以表示低频的第一标量系数来除并以表示高频的第二标量系数来除而实现。第一标量选择得小于第二标量，这是由于标量越大将实现的压缩越大。这样，在地震分析中特别有意义的表示低频的系数将不会象表示高频的系数那样被压缩，从而保持前者的精度。
通过使用最接近整数函数的均匀量化，或统计舍入，或非均匀重新量化能够实现定标。
定标或量化参数可随时间、空间、或空间或时间的频率而变化。
在本发明的进一步的优选实施例中，压缩比是由地震信号中的噪声水平自动确定的。最好使用不包含由地震源产生的信号的一部分信号或多个部分信号来测量噪声水平。这样，在“第一到达”之前记录的跟踪图形的一部分或所谓噪声记录，即在没有地震源时的记录的跟踪图形，能够用来确定噪声水平。甚至更好的是，在确定噪声水平的步骤之前对信号滤波，以避免对噪声水平的过分估计。
最好选择压缩比，即量化误差，使其等于或低于噪声水平。
由数据压缩引起的精度的降低增加了数据的冗余度。因而根据本发明进一步的一个特点，数据的冗余度可用来进一步缩减被压缩的数据，这最好是通过施加可变长度数据编码，诸如幅度编码或Huffman编码。例如，可对定标系数进行幅度编码，使得每一系数与最大系数的绝对值成正比。
根据本发明的地震数据压缩方法可用于所有类型的地震数据，包括陆地，过渡带，海洋或海底探测的2维和3维观测数据。可能的数据还包括预记录数据或预处理数据，诸如冲击收集(shot gathers)，共用中点(CMP)收集，堆(stacks)，迁移段，或单传感器记录。还适用于二维地震探测几何结构，诸如包括一条线接收器的陆地布局，或对于使用浮筒与水下测音器的二维或三维海洋地震观测，及对于实际区域接收器布局(线内／交叉线应用)的两个垂直方向的应用。
使用专用芯片组，能够做到在向野外工作室(field boxes)或探测车或船传输之前在单个接收器内压缩数据。还能够在数据传输中在数据处理的更远的“下游”使用，例如对于数据向处理中心传输或对于中间阶段使用。在传输和／或存储压缩数据之后，能够使用相同的反向步骤并以反向顺序重构或解压缩出原始数据。
本发明还在于用来执行这里所述的方法的装置。
本发明还在于数据受到根据上述一种或多种方法压缩的地震观测方式。
从以下的详细说明和附图，本发明的这些或其它的特点，优选实施例及其变种，可能的用途和优点能够为本领域普通专业人员所接收和理解。
图1表示根据优选的窗口功能设计的典型的窗口。
图2表示描述当使用根据本发明的方法时在空间频域中预计的定标或量化误差的曲线。
图3表示描述当使用根据本发明的方法时在空间频域中实际的定标或量化误差的曲线。
一般来说，在传统的3维陆地地震学布局，接收器安排成数个平行线的直线接收器阵列中。在这一实施例中，数据压缩方法是沿着每一接收器行施加的，且每一接收器行可以被单独处理。这种应用称为按行的应用。
传统的分组形成法过程，是对各组接收器求和，在数学上表示为[1]g→(i)=1NgΣn=0Ng-1S→(n+iNg);]]>i=O,...,Ig-1]]>其中Ng是每组接收器的个数，Ig是每接收器行的组数，向量g包含在第i组中成组形成的数据，向量s(n)包含接收器n中测量的数据。
在方程式[1]中可见，分组形成数据中的样本数只是原始数据量的一个因子1／Ng。
根据本发明的地震数据压缩方法不依赖于这一分组形成法，这种方法是只保留接收信号的空间低频成分的粗糙方法。另一方面，该方法避免了保留来自各个接收器的所有数据。新方法以高的精确性保留了更重要的空间低频成分，并以降低的(但仍然是重要的)精确性保留了空间高频成分。精确性的降低意味着每样本数据需要较少的位，因而实现了数据的压缩。
根据本发明的方法涉及以下将说明的数个阶段阶段1空间变换Ⅳ型局部余弦变换，即向数据施加有限数目接收器的窗口内的余弦变换，C→km[2]=Σn=m/23M/2-1s→(n+nM)h(n)2Mcos(πM(k+12)(n+12))]]>对于k=0，…，M-1，且对于m=0，…，P-1。
在等式[2]中，向量Ckm是局部空间DCT-IV系数，向量S(n)是接收器n中所测量的数据，M是局部余弦变换中每窗口内的接收器数目，h(n)是窗口函数，而P是每接收器行的窗口数。
在这一具体的方程式中，假设M是偶数，当然对变换作适当的修改M可选择为奇数。如果窗口函数h(n)满足以下条件，则以上变换是正交的且可逆的0≤h(n)≤1h(n)=0for n≤-M2and n≥3M2,]]>[3]h(n)=h(M-n-1),h(n)2+h(n+M)2=1for-M2≤n<M2·]]>
公式[2]中的变换称为DCT-TV(第四型离散余弦变换)，并具有与快速傅立叶变换(FFT)可比拟的有效实现性。公式[2]中的变换将称为局部DCT-IV。对于局部DCT-IV，计算的复杂性与窗口数目的乘积成正比，且执行DCT-IV所需的工作量与(N／M)*M log(M)∝N成正比，其中N是信号长度，即每接收器行的接收器数目，M是窗口长度。这与FFT比较更有利，后者是对信号的整个长度进行的，将需要与N log(N)成正比的多次操作。这一变换与普通的局部DCT之间的主要差别在于它能够用于重叠的窗口。普通的局部DCT被限制在不相交的矩形窗口中。
在等式[2]中可见，每一组系数{Ckm；k=0，…，M-1｝的计算需要2M个接收器的数据贡献，其中的M来自本身窗口内，M／2来自每一邻接的窗口。但是，全部的变换系数的个数等于原始数据样本数，即p*M=N。数据的开始和结束可通过假设周期性或使用分开的开始和结束窗口函数处理。
在多个窗口内可同时计算变换系数。此外，存在一种分两步骤执行等式[2]中的变换的有效实现方法1．确定每一窗口m中的折叠信号向量fm的折叠步骤对于0≤n＜M／2；以及[4′]s→(n+mM)h(n)-s→(2M-n-1+mM)h(2M-n-1)]]>对于M／2≤n＜M。
2．折叠信号fm的余弦变换是这样的[5]c→km=Σn=0M-1f→m(n)2Mcos(πM(k+12)(n+12))]]>
对于k=0，…，M-1，以及m=0，…，P-1。
在IEEE ASSP，vol．38，No．6．，June 1990，Lapped transform forefficient tranform／subband coding中，H．S．Maler讨论了局部DST-IV(第四型局部正弦变换)的性质和实施方法，这一变换除了余弦由正弦代替之外与局部DCT-IV相等。
阶段2时间变换局部DCT-IV系数仍然是记录时间的函数。在所提出的方法中，也向它们施加局部时间DCT-IV。窗口的长度和窗口函数是从局部空间DCT-IV中所使用的窗口的长度和窗口函数中各自独立地选择的。
虽然其它纵坐标不再简单地指时间，但是向量Ckm的表示法不变。
阶段3重新量化正交可逆变换的方便性是，它们是保持能量的，即它们满足Parseval定理。这意味着，变换域中的平方的量化误差等于原始域中求得的平方误差。对于原始域中的平方误差同样是这样。对于所有的数据能量的平方误差也是这样。然而，这种保持性质对于数据的最大绝对振幅并不存在。
量化涉及到划分诸如系数等的数据振幅的范围，并根据所使用的具体的量化方法通过舍入振幅以便指定另一个振幅值而降低出现的数据量。
最普通的量化方法是均匀量化方法振幅范围被划分为相等的级差，并对中点舍入振幅。这样作得到的结果是数据样本的固定点(整数)表示。如果级差大小(比如△)足够小，则量化误差是以方差或能量△2／12均匀分布的白噪声。如果在变换域中应用均匀量化方法，则原始域中的量化误差也将表现为不变能量的白噪声。原始域中量化误差不一定是均匀分布。如果使用非均匀量化器，则这些说法不成立。可以设想对于增加的振幅具有下降精度的量化器，诸如浮点表示中使用的量化器。
如果达到大的压缩比，则量化器对于(在变换域中)数据的某些部分变为粗糙的。这时误差的随机分析开始失效，并出现对变换域中的数据的滤波。然而如果所使用的变换给出地震数据内容的良好的压缩，诸如象使用DCT-IV那样，则能够避免数据的最重要的内容受到这种滤波的作用。
局部空间DCT-IV的系数表示地震数据的局部空间频率内容，即在m号窗口中系数向量Ckm，K=0，…，M-1表示空间频率内容。指标k决定了所考虑的空间频率。低的k对应于低的局部空间频率，高的k对应于高的局部空间频率。
该方法涉及低和高的局部空间频率内容的分开的重新量化。这可通过对Ckm’规则舍入(或均匀量化)进行以便给出[6]c→km=NINT(c→km/δL)δL]]>对于0≤k≤km-1，以及[7]c→km=NINT(c→km/δH)δH]]>对于km≤k≤M-1，其中方程式[6]用于低局部空间频率，而方程式[7]用于高局部空间频率。这两个方程式中，NINT是最近整数函数，km是高精度的局部DCT-IV系数的数目，而δL，δH是用于对空间频率舍入的标量。在以上的表达式中，在δL，δH变得较大时量化误差变得较大。为了保持低局部空间频率在精度上高于高局部空间频率内容，选择这些标量使得δL＜δH。标量或量化误差越大，则需要的每样本的位数越少，从而能够达到更大的压缩。
因而，调节这些标量的值提供了自动选择压缩比的一种方法。首先，调节标量的一个方法是要确定对被记录的地震数据中噪声的一个估计值。这可通过比较已知没有地震信号的记录信号部分来实现，即比较最好是在第一次到达之前所记录的数据部分，或在所谓噪声冲击期间记录的数据。从这一“无信号”数据，可通过通常的统计方法推导出噪声的一个估计值。如果给出这一估计值，就能够按和它的预定的关系对压缩的标量进行调节。而且，精度降低的系数越多(即km越小)，则能够达到的压缩越大。例如，若给定-10dB的数据信号电平和-50dB的噪声电平(使用来自无信号的数据部分的数据测量)，则可设置这些标量，以便至少达到18∶1的压缩比。
在局部的和开窗口的方法中，在低的(传统的)和高的(添加的)波数带之间所需要的分离不可能是完美的。从高到低的波数带的量化噪声小量的泄漏是不可避免的。为了使噪声泄漏最小，要仔细选择对于局部余弦系数中量化噪声的窗口设计。还要通过牺牲某些压缩性能来达到满意的低水平泄漏。在以下的阶段0讨论窗口的设计。
阶段4振幅编码为了进一步压缩数据，可以采用降低精度，因而降低了每样本的位数，这增加了数据中的冗余度。这由对变换系数(Ckm’)的振幅编码所利用。为此目的，形成少量(通常为8个)系数的运算。在每一运算中，最大的绝对值决定了系数所要使用的位数。代码中，必须把系数的每一运算放在所使用的位数的后面，或者更具体来说，放在每次运算所需的位数的Huffman编码获得的代码的后面。
阶段0参数设置和窗口设计这一阶段需要选择窗口函数和重新量化参数。
该方法保留数据的低的空间频率内容，同时降低高的空间频率内容的精度。正如前面所述，这是通过改变表示地震数据的局部空间频率内容的局部DCT-IV系数的精度做到的。在以上按行应用中，理论上能够预测重新量化局部空间频率(即局部空间DCT-IV系数)对沿整个接收器行的数据的空间频率内容的作用效果。
在空间频带中，在被认为低的和被认为高的之间作出区别。在通过整个接收器行观察时，数据的低的空间频率内容的精度使高的局部空间频率的精度降低，这是不可避免的。在低空间频率引起的精度的损失由窗口函数h(还涉及其长度2M)、高精度km的局部DCT-IV系数的数目、以及由δL和δH决定的所选择的精度确定。
该方法涉及由以下步骤组成的窗口设计过程-规定沿接收器行的低空间频带；-最好与数据中的噪声预定估计值相关地设置低空间频带中最大可接受精度损失的阈值；-选择M，km，δL和δH；-对所有允许的窗口函数h(从方程式[3])使低空间频带的精度损失最小；-如果所得的精度损失低于阈值，则终止该过程，否则按高精度(km)增加局部DCT-IV系数的数目并重复前面的步骤。
影响量化噪声电平(△)的泄漏量的量为局部截止指标(km)，用于DCT-IV的长度M、及最后是窗口长度(≤2M)和形状。为了使泄漏最小，要求△尽可能地小，且km和M尽可能地大。然而为了使压缩最大，要求△尽可能地大且km尽可能地小。在窗口设计中M和△是固定的。使用对km的初始选择，针对窗口h把噪声泄漏减小到最小。如果噪声超过泄漏阈值，则增加截止指标km并计算新的窗口。
现在提出使用以下各项小规模实际数据的例子-128个接收器，因而空间频率数为64；-空间窗口长度M=16，因而空间窗口函数长度2M=32；-空间窗口数P=8；-在4ms时间每接收器有1024个样本；-时间窗口长度为64个样本；-24位定点样本值。
选择低空间频带使其由4个最低空间频率(8个实值傅立叶系数)组成。这意味着，由于窗口数为8，故每窗口高精度(即km)局部DCT-IV系数(即，Km)的个数至少为1。
设置低空间频带量化误差的阈值为-115dB。设置低量化误差为-119dB且高量化误差为-68dB。在使用这一例子初步实验期间，如果高精度(km)局部DCT-IV系数的数目从1增加到4，则误差只落到阈值以下。这意味着以高精度保留了系数的25％(每16个窗口4个)。图1中展示了所设计的窗口。预测的量化误差与表示低和高空间频率之间分割的矩形曲线一同示于图2中。
从实际数据所获得的量化误差在图3中示出。在预测的和实际的量化误差之间可观察到良好的吻合。
原始数据为24位固定点。用来获得所需精度的每局部DCT-IV系数的位数对于低DCT-IV系数(25％)为26个，对于高DCT-IV系数(75％)为17个。这是一个平均为每样本19．25位的初始压缩。然而，增加的冗余度能使得在阶段4的振幅编码把这降低到每样本4．4位。
在另一实施例中，主要目的不是使数据压缩方法与用传统方法获得的数据相抵触。这意味着大量的分组接收器要以高精度保留，而其余的地震数据以降低的精度表示。这一实施例的实施方案比第一实施例简单。涉及FFT和线性阵列的该方法如下阶段1空间变换[8]f→km=Σn=0M-1s→m(n)e exp(-2πinkM)]]>对于k=0，…，M-1，并对所有的m。
窗口大小(M)简单地等于分组大小(Ng)。这里分组数目作为数据向量Sm(n)中的下标m给出以表示所有的分组被分开处理，并且它们可以与接收器行成一角度。
阶段2时间变换对于时间变换，如同第一实施例中那样使用时间局部DCT-IV(3．1节)。
阶段3重新量化类似于第一实施例中的方程式[6]和[7]那样进行重新量化，例如[9]f→km=NINT(f→Om/δL)δL]]>(和)以及[10]f→km=NINT(f→km/δH)δH]]>对于1≤k≤M-1(高局部空间频率)。
阶段4振幅编码如同第一实施例中所述进行振幅编码。
由于这一实施例中，只以高精度保留大量的接收器分组，不限制使用线性阵列的接收器行的地震探测几何结构。例如可以应用到区域阵列的探测几何结构。
其它实施例如同以上讨论的实施例，但是有以下特点a)阶段1和2互换；b)阶段3中重新量化方程式[6]和[7]中的系数的任何其它方法。例如，统计舍入或非均匀重新量化；c)方程式[6]和[7]中重新量化参数随时间、位置或频率而变化；d)奇窗口长度M；e)在局部空间或时间DCT-IV中窗口长度不变；f)DST-IV代替DCT-IV；g)阶段1中的局部时间DCT-IV由任何其它信号变换或分解代替，诸如用(局部)DCT、局部DST、(局部／短时间)FFT、子波变换或子频带分解代替；h)局部空间DCT-IV由任何其它信号变换或分解代替，诸如用(局部)DCT，局部DST，(局部／短时间)FFT，子波变换或子频带分解；i)阶段4中的振幅编码由利用数据冗余度的任何其它方法代替。
权利要求
1．基于离散三角变换的地震数据压缩方法，所述方法的特征在于以下步骤-选择窗口函数，使得变换施加到中心窗口并与邻接窗口有重叠；-向所述数据施加Ⅳ型(DCT-IV，DST-IV)局部空间和／或时间离散三角变换以便在变换域中产生数据；以及-压缩所述变换域中的数据。
2．根据权利要求1的方法，其中按降低的精度通过保留其所选择的部分而压缩变换域中的数据。
3．根据权利要求2的方法，其中按降低的精度保留变换域中数据的高频部分。
4．根据权利要求1的方法，其中通过对变换域中数据进行较粗糙的重新量化而压缩变换域中的数据。
5．根据权利要求1的数据压缩方法，其中确定非压缩数据中的噪声估计值，并根据所述估计值自动选择对于该压缩的压缩比。
6．根据以上任何权利要求的数据压缩方法，其中变换域中的数据的降低的冗余度用来进一步缩减被压缩的数据。
7．根据权利要求6的数据压缩方法，其中通过可变长度编码压缩变换域中的数据。
8．权利要求1，4或6的方法，还包括以下步骤-传输和／或存储压缩的数据；以及-解压缩并反向变换所述压缩的数据。
9．用于压缩地震数据的装置，所述装置包括-产生窗口函数使得变换施加到中心窗口并与邻接窗口有重叠的装置；-用于向所述数据施加Ⅳ型(DCT-IV，DST-IV)局部空间和／或时间离散三角变换以便在变换域中产生数据的装置；以及-用于压缩所述变换域中的数据的装置。
全文摘要
使用施加于数据集的Ⅳ型局部离散余弦/正弦变换,描述了用于地震数据压缩的方法和装置,数据集由不变的窗口函数限制,该窗口函数降低了与邻接窗口的重叠。
文档编号G01V1/22GK1206470SQ9619933
公开日1999年1月27日 申请日期1996年11月22日 优先权日1995年12月1日
发明者皮特·伦纳德·弗米尔 申请人:格库公司