2010, Vol. 53
2. 长江大学地球物理与石油资源学院, 荆州 434023
2. School of Geophysics and Oil Resources, Yangtze University, Jingzhou 434023, China
地震勘探是最重要的石油勘探方法之一，随着石油勘探开发工作的不断深入，地震勘探正向多维、多分量、高精度方向发展[1, 2]，使得勘探数据呈爆炸性趋势增长.例如，2008年9月，在胜利油田罗家地区实施的高密度三维三分量（3D3C）地震先导试验采集中，仪器记录道数33600道，单炮数据量达到940.8 M，形成了海量数据[2].不仅如此，近几年在地震勘探数据处理的方法研究方面所使用的实际地震资料的数据量也非常大.2008年，常旭等[3]在研究深水崎岖海底地震数据成像方法时，所选用的实际地震数据是一条近40km长的二维测线，共1521炮，每炮198道，记录长度6s，采样间隔4ms，每道1500样点.若按每样点4字节计算，这条二维测线的地震数据量就高达约1.8G字节.
海量数据给地震勘探数据的传输、存储和处理提出了严峻挑战.具体表现为，给地震数据的存储所带来的主要问题是要占用极大的存储空间，从而增加存储成本；而给传输所带来的主要问题是极大地增加了传输带宽需求，同时也增加了传输时间.更严重的是因数据量太大、传输时间过长往往会导致网络断线而必须重传，甚至可能完成不了完整的数据传输，以致所接收到的数据无法使用，这给跨国或跨探区作业的地震勘探质量监测带来了困难.地震数据压缩是解决海量地震数据传输和存储问题的关键.虽说地震数据的压缩和图像数据的压缩有很多相似之处，但又有许多不同点，如地震数据的动态范围大、中低频信息丰富等[4].因此，地震勘探数据压缩方法的研究受到国内外学者的广泛关注[5~10].地震勘探数据压缩方法按照有无信息丢失可分为无损压缩[11~13]和有损压缩[6, 14, 15]两大类.无损压缩又被称为无失真压缩，其基本原理是通过去除原始数据中的冗余，用尽可能少的码字来描述原始数据，从而减少数据的量[5].顾名思义，这类方法的最大优点就是能够保证压缩重建后的数据没有任何的信息丢失，因此特别适合于对传输或存储后还需要进一步处理的地震数据进行压缩.这类方法的弱点是压缩效率低，一般压缩倍数在2左右[5].有损压缩是一种以牺牲部分信息获得较高倍压缩的方法，这类方法的基本原理是利用某种正交变换（如离散余弦变换DCT，小波变换等）将原始数据变换到一种更易于压缩的域（如频率域），并对变换域的数据进行量化，然后对量化后的变换域数据进行无失真的熵编码.信息的损失主要是由量化所引起的.这类方法适合于对传输或存储后不需要再处理的地震数据进行压缩，例如用于地震勘探质量监测.
地震勘探数据一般都是以SEG（Societyof Exploration Geophysicists）的标准格式（如SEG D，SEG Y等）记录和存储的，按照这些格式所记录的数据文件中除了包含地震样点数据之外，还包含头段数据，以对诸如文件中所记录的道数、每道样点数、样点数据格式、覆盖次数、采样间隔、仪器因素等很多重要信息进行必要的说明.为了保证对地震数据的正确处理和使用，不允许头段数据出现任何差错.因此，无论是有损还是无损的地震数据压缩，对头段数据的压缩都必须是无失真的.头段数据在整个地震数据文件中占有一定比例，能否有效压缩头段数据将对整个地震数据文件压缩倍数产生重要影响.本文将定义头段数据所占比例，研究头段数据压缩倍数与地震数据文件压缩倍数之间的关系，并据此阐明对头段数据进行高倍压缩的必要性.在此基础上，本文将重点研究SEG Y格式头段数据的特点，并根据这些特点提出适合于对SEG Y头段数据进行高倍无失真压缩的方法，对提高SEG Y数据文件的压缩倍数具有极其重要的实际意义.
2 SEG Y格式简介地震勘探数据交换格式在地震数据采集和处理中发挥着非常重要的作用.为了满足勘探技术发展的需要，SEG分别于不同时期发布了多种不同的地震数据交换标准，包括SEG A，SEGB，SEGX[16]，SEGC[17]，SEG Y[18]和SEG D[19].这些标准被分别简称为A格式，B格式，X格式，C格式，Y格式和D格式.其中D格式分别于1994年、1997年、2006年和2009年做过4次修改，其最新版本是2009年的SEG-Drev3.0；而Y格式于2002年升级为1.0版[20].同时，1975年版的Y格式[18]被改称为SEG Yrev0.目前在地震数据采集和处理中最常用的数据格式分别是D格式和Y格式，其中SEG Yrev0仍然在普遍使用.本文主要介绍1975年版的Y格式（即SEG Yrev0），为方便起见，以下将其简称为SEG Y或Y格式，并将Y格式地震数据集称为SEG Y数据文件或SEG Y记录.
如图 1所示，SEG Y数据文件主要包括卷头（Reel identification header）数据块（灰色背景）和地震道数据两大部分.其中，卷头数据块共3600字节，又分为3200字节的EBCDIC[18]卡和400字节的二进制编码头两部分；而地震道数据包括m个地震道数据块（Trace data block），每个地震道数据块由240字节的道头数据和多个样点数据组成.根据样点数据的类型不同，每个样点数据由4字节（浮点型和定点4字节型）或2字节（定点2字节型）构成.本文将卷头数据块和所有道头数据统称为头段数据.下面对头段数据作进一步介绍：（1）EBCDIC卡从第1字节开始到第3200字节结束，包含40行×80列的文本字符，可转换为ASCII后显示；（2）从第3201字节开始到第3600字节结束共400字节是二进制编码数据，其中后340个字节（3261~3600）作为备用信息未被定义，而在前60个字节（3201~3260）中定义了对整个SEG Y数据文件有效的信息，包括：工作识别号（Job identification number）、测线号（Line number）、卷号（Reel number）、地震道数、辅助道数、采样率、每道样点数、样点数据类型、覆盖次数等27个数据项.除工作识别号、测线号和卷号3个数据项各为4字节之外，其余24项皆为每项2字节；（3）每个道头的240字节中，前180字节定义了71个数据项（52个2字节项和19个4字节项），后60个字节备用未被定义，本文将这60个字节划分为15个4字节项，因此道头共有86个数据项；在71个SEG定义的数据项中，有6项是SEG建议必须记录的，它们是：测线中的道序号（第1~第4字节），原始野外记录号（第9~第12字节），原始野外记录中的道号（第13~第16字节），道识别号（第29~第30字节），样点数/道（第115~第116字节），本道采样间隔（第117~第118字节）.关于卷头和道头的更详细介绍请参阅SEG Y标准[18].
|
图 1 SEG Y格式数据文件构成 (ni, 为第i道中的样点数，一般情况下n1=n2=…=nm=n) Fig. 1 SEG Y format (ni is the number of samples in trace i, usually, n1=n2=…=nm=n) |
用H，D和V分别表示SEG Y格式原始头段数据量（字节数）、原始地震道样点数据量（字节数）和原始SEG Y数据文件的数据量（字节数），并用m，n和k分别表示SEG-Y数据文件中的地震道数、样点数/道和字节数/样点，则
|
(1) |
|
(2) |
|
(3) |
其中，3600为卷头（Reel Header）字节数，240为每道的道头字节数，k视样点数据的类型可以为2或4.
用α表示H在V中所占比例（称为头段数据所占比例，简称头段比例），则
|
(4) |
当240m>>3600，且nk>>240时，
|
(5) |
240m>>3600的物理意义是，SEG Y文件中所有地震道的道头数据量之和远大于卷头数据量；而nk>>240的物理意义是，地震道中各样点数据量之和远大于道头数据量.一般SEG Y格式地震数据文件都满足这两个条件.实际上，240m>>3600等价于m>>15，而且对任何k值的情况，只要n>>120则nk>>240成立.所以（5）式说明，只要地震道数和每道样点数足够大，则头段数据所占比例α只与样点数/道和字节数/样点有关，而与地震道数无关.特别地，对于样点数据为浮点型或定点4字节的情况，k=4，此时α≈60/n；而对于样点数据为定点2字节的情况，k=2，则α≈120/n.例如，对于道数为2000，每道3000个样点的浮点型SEG Y文件（即m=2000，n=3000，k=4），利用（5）式计算得到α=0.02，与利用（4）式计算得到的α=0.01975199726非常接近.
3.2 压缩倍数公式本文将对SEG Y数据文件的压缩倍数称为文件压缩倍数（或总压缩倍数），它等于SEG Y数据文件压缩前的字节数与压缩后的字节数之比，用Cr表示；将对头段数据的无失真压缩倍数称为头段压缩倍数，它等于SEG Y数据文件中的头段数据压缩前的字节数与压缩后的字节数之比，用CrH表示；而将对所有样点数据进行压缩的倍数称为样点压缩倍数，它等于SEG Y数据文件中所有的样点数据压缩前的字节数与压缩后的字节数之比，用CrD表示.若用Hc，Dc和Vc分别表示对原始头段数据、原始地震样点数据和原始SEG Y数据文件进行压缩后的数据量，则
|
(6) |
|
(7) |
|
(8) |
|
(9) |
由（4）式和（6）~（7）式可得
|
(10) |
其中，
|
(11) |
当α<<1时，（10）式也可近似为
|
(12) |
（10）式表明，SEG Y文件压缩倍数与头段比例（α）、头段压缩倍数（CrH）和样点压缩倍数（CrD）有关，并以CrLim（CrH与α之比）为理论极限.也就是说，对于一定的CrH和α，当CrD趋于∞时Cr的极限为CrLim.因此，寻求适合于SEG Y格式头段数据的高倍无失真压缩方法，对于提高SEG Y数据文件的总压缩倍数具有极其重要的意义.例如，对于α=0.02的情况，若不对头段数据进行压缩（即CrH=1），则总压缩倍数的极限为CrLim=50，而当样点数据的压缩倍数CrD高达50时，总压缩倍数Cr仅为25.3；即使CrD高达100，Cr也只有33.6；另一方面，要是能对头段数据实现30倍压缩，还是对于α=0.02，就有CrLim=1500，此时，若CrD=50，Cr就能达到49.3.
4 SEG Y格式地震数据头段部分的特点分析通过对江苏油田和江汉油田的10个SEG Y数据文件进行分析，发现SEG Y格式卷头数据具有如下特点：（1）在3200字节的EBCDIC部分存在大量连续的空格字符（EBCDIC码为0X40，所对应的ASCII码为0X20），尤其是在第C23~第C40行之间，绝大部分都为空格字符；（2）若将400字节的二进制编码部分看作400个单字节的整数，其中存在大量的零值，主要表现在除了已定义的60个字节中存在部分零值外，其余未被定义的340个备用字节基本上全都为零；（3）若将整个卷头看作3600个单字节整数，记为R1，R2，…，R3600，另外令R0=0，定义Ri的一阶差分为
|
(13) |
并对ΔRi（i=1，2，…，3600）进行统计，其中零值比例约占67%，而且出现连续零的几率较大，在ΔRi（i=1，2，…，3600）中出现高达343个连零的现象基本上在每个文件中都存在.
每一道头内的86个数据项之间存在着相关性，例如，测线中的道序号（第1~第4字节）与卷中的道序号（第5~第8字节）一般是相同的.但是，道头数据的冗余主要表现为相邻地震道的道头数据的对应数据项之间的强相关性.设SEG Y数据文件中共记录m个地震道，用xi，1，xi，2，…，xi，86表示第i道（i=1，2，…，m）的道头数据中的86个数据项，令x0，1=x0，2=…=x0，86=0和x-1，1=x-1，2=…=x-1，86=0，定义道头数据项的二阶差分为
|
(14) |
通过对江苏油田和江汉油田各5个SEG Y数据文件（共10个文件）的道头数据进行统计分析，发现绝大多数Δ2xi，j（i=1，2，…，m；j=1，2，…，86）为零，并且从第3道至第m道的Δ2xi，j（i=3，4，…，m）全为零的数据项很多，在86个数据项中平均只有7.9个数据项的Δ2xi，j（i=3，4，…，m）不全为零.即二阶差分不全为零的数据项仅占约9.2%，这一统计特性为有效地压缩道头数据提供了依据.
5 头段数据无失真预测编码本文的头段数据无失真压缩分为两个部分：卷头编码器和道头编码器.卷头编码器是根据前节所述的卷头数据的特点而设计的，具体编码过程如下：
（1）将整个卷头看作3600个单字节整数数据项，并根据（13）式计算每个数据项所对应的一阶差分；
（2）从第1项开始，依次对各差分值进行编码：
（a）当差分值为零时使用游程编码；
（b）当差分值非零时使用变长编码（VLC）.
卷头解码是卷头编码的逆过程.
同样，道头编码器是根据前节所发现的道头数据的统计特性而设计的，设共有m个地震道，具体编码过程如下：
（1）设置86个标志位，记为Bf（j）（j=1，2，…，86），分别与86个数据项相对应；
（2）根据（14）式计算所有地震道的与86个数据项相对应的二阶差分；
（3）对第j个数据项（j=1，2，…，86），判断Δ2x3，j=Δ2x4，j=…=Δ2xm，j=0是否成立？若是，则置Bf（j）=0；否则，置Bf（j）=1；
（4）按位输出（或存储）Bf（j）（j=1，2，…，86）；
（5）对Δ2x1，j和Δ2x2，j（j=1，2，…，86）进行变长编码；
（6）对每个j（j=1，2，…，86），若Bf（j）为1，则依次对二阶差分值Δ2x3，j，Δ2x4，j，…，Δ2xm，j进行编码：
（a）当二阶差分值为零时使用游程编码；
（b）当二阶差分值非零时使用变长编码（VLC）.
道头解码是道头编码的逆过程.实际上无论是对一阶差分值还是对二阶差分值的编码都等价于对某种预测器的预测误差的编码，于是卷头编码器和道头编码器都是基于预测方法的编码器，所以本文将这种编码方法称为头段数据无失真预测编码（Header Identification Data Lossless Prediction Coding，简记为HIDLPC）方法.
6 实验结果及分析为了验证本文所提出的HIDLPC方法，作者不仅开发了HIDPLC软件用于对SEG Y头段数据进行压缩和解压缩，而且还另外开发了两个辅助程序：（1）SEG Y头段数据提取程序；（2）SEG Y头段数据文件比较程序.其中，头段数据提取程序的作用是从SEG Y数据文件中将卷头数据和道头数据提取出来并另存为头段数据文件，以便用HIDPLC软件、Winzip软件和WinRAR软件对其进行压缩；而SEG Y头段数据文件比较程序的作用是对压缩前和经过压缩重建后的SEG Y头段数据文件进行比较，看是否存在差别.用头段数据提取程序所另存的头段数据文件如图 2所示.
|
图 2 头段数据文件构成 Fig. 2 File structure of the header data |
利用本文所设计的HIDLPC软件对江苏油田和江汉油田的16个SEG Y数据文件（如表 1所示）的头段数据进行了压缩和解压缩，经比较确认皆为无失真压缩.在表 1中，m、n、k、V和H分别表示道数、样点数/道、字节数/样点、原始SEG Y文件大小（字节数）和原始头段数据文件大小（字节数）；以JS开头的SEGY文件来自江苏油田，而以JH开头的SEG Y文件来自江汉油田，其中JS1~JS5和JH1~JH5曾在本文第4节中用于对SEG Y头段数据进行统计分析.表 2中列出了分别用HIDLPC软件、Winzip和WinRAR对这16个文件进行无失真压缩的压缩倍数，可以看出，HIDLPC对16个文件的压缩倍数有12个超过了100倍，最低压缩倍数为38.64，最高压缩倍数达到1077.72，并且全部都远高于Winzip和WinRAR的对应压缩倍数；根据表中结果还可计算出HIDLPC、Winzip和WinRAR对16个文件的平均压缩倍数分别为约200、22和36.由此可见，本文所提出的HIDLPC方法是一种能够对SEG Y头段数据实现高倍压缩的方法.与Winzip和WinRAR相比，HIDLPC最适合于SEG Y头段数据的无失真压缩，而WinRAR比Winzip更适合于SEG Y头段数据压缩.
|
表 1 用于压缩实验的16个SEG Y数据文件 Table 1 The 16 SEG Y files used for compresion experiments |
|
|
表 2 SEG Y头段压缩倍数比较 Table 2 Comparison of the compression ratios of SEG Y header data |
本文对SEG Y头段数据进行了研究和分析，发现了卷头数据和道头数据各自的特点，并在此基础上提出了一种头段数据无失真预测编码（HIDLPC）方法，实验结果表明，该方法的平均压缩倍数达到约200，最高压缩倍数超过1000，无论是平均压缩倍数还是对具体SEG Y头段数据文件的压缩倍数都远高于Winzip和WinRAR的相应压缩倍数.这使得对SEG Y数据文件的总压缩倍数的理论极限大大提高，对研究高效的SEG Y文件压缩方法具有重要的实际意义.
| [1] | 熊翥. 高精度三维地震(Ⅰ):数据采集. 勘探地球物理进展 , 2009, 32(1): 1–11. Xiong Z. High precision 3-D seismic:Part I Data acquisition. Progress in Exploration Geophysics (in Chinese) , 2009, 32(1): 1-11. |
| [2] | 赵殿栋. 高精度地震勘探技术发展回顾与展望. 石油物探 , 2009, 48(5): 425–436. Zhao D D. Review and prospect on high-precision seismic exploration technique. Geophysical Prospecting for Petroleum (in Chinese) , 2009, 48(5): 425-436. |
| [3] | 常旭, 刘伊克, 杜向东, 等. 深水崎岖海底地震数据成像方法与应用. 地球物理学报 , 2008, 51(1): 228–234. Chang X, Liu Y K, Du X D, et al. Seismic imaging under the irregular deep water bottom. Chinese J. Geophys. (in Chinese) , 2008, 51(1): 228-234. |
| [4] | 张学工, 李衍达, 盛硕. 几种现代信息处理技术和计算机网络在勘探地球物理中的应用. 地球物理学报 , 1997, 40(Suppl.): 275–291. Zhang X G, Li Y D, Sheng S. Application of some modern information processing and computer network techniques in exploration geophysics. Chinese J. Geophys. (Acta Geophysica Sinica) (in Chinese) , 1997, 40(Suppl.): 275-291. |
| [5] | 武文波, 胡福祥, 徐凌, 等. 地震数据压缩方法与应用. 北京: 石油工业出版社, 2007 . Wu W B, Hu F X, Xu L, et al. Seismic Data Compression Methods and Applications (in Chinese). Beijing: Petroleum Industry Press, 2007 . |
| [6] | 余平, 马小虎, 陈恒金. 基于提升小波的地震数据压缩编码算法. 苏州大学学报(工科版) , 2009, 29(1): 7–11. Yu P, Ma X H, Chen H J. The encoding algorithm of seismic data compression based on lifting wavelet transform. Journal of Suzhou University (Engineering Science Edition) (in Chinese) , 2009, 29(1): 7-11. |
| [7] | 刘财, 王培茂, 王世煜, 等. 联合数据压缩技术在地震勘探中的应用. 世界地质 , 2006, 25(4): 434–439. Liu C, Wang P M, Wang S Y, et al. Application of united data compression technique in seismic prospecting. Global Geology (in Chinese) , 2006, 25(4): 434-439. |
| [8] | Wood L C. Seismic data compression methods. Geophysics , 1974, 39(4): 499-525. DOI:10.1190/1.1440443 |
| [9] | Jonsson S B, Spanias A S. Seismic data compression. IEEE Conference Proceedings-Annual Phoenix Conference, 1990. 276~279 |
| [10] | RØsten T, Ramstad T A, Amundsen L. Optimization of sub-band coding method for seismic data compression. Geophysical Prospecting , 2004, 52(5): 359-378. DOI:10.1111/gpr.2004.52.issue-5 |
| [11] | Steams S D, Tan L, Magotra N. A technique for lossless compression of seismic data. International Geoscience and Remote Sensing Symposium, 1992(IGARSS'92), 1992. 681~683 |
| [12] | Mandyam G, Magotra N, McCoy W. Lossless seismic data compression using adaptive linear prediction. International Geoscience and Remote Sensing Symposium (IGARSS), vol. 2, 1996. 1029~1031 |
| [13] | Abanmi A O, Alshebeili S A, Alamri T H. Lossless compression of seismic data. Journal of the Franklin Institute , 2006, 104(6): 340-351. |
| [14] | 刘财, 王培茂, 薛亚丽, 等. 离散余弦变换(DCT)编码在地震勘探数据压缩中的应用. 吉林大学学报(地球科学版) , 2004, 34(2): 277–282. Liu C, Wang P M, Xue Y L, et al. The application of discrete cosine transform coding to seismic data compression. Journal of Jilin University (Earth Science Edition) (in Chinese) , 2004, 34(2): 277-282. |
| [15] | Spanias S, Jonsson S B, Stearns S D. Transform methods for seismic data compression. IEEE Trans. Geoscience and Remote Sensing , 1991, 29(3): 407-416. DOI:10.1109/36.79431 |
| [16] | Northwood E J, Weisinger R C, Bradley J J. Recommended standards for digital tape formats. Geophysics , 1967, 32(6): 1073-1084. DOI:10.1190/1.32060004.1 |
| [17] | Meiners E P, Lenz L L, Dalby A E, et al. Recommended standards for digital tape formats. Geophysics , 1972, 37(1): 36-44. DOI:10.1190/1.1440250 |
| [18] | Barry K M, Cavers D A, Kneale C W. Report on recommended standards for digital tape formats. Geophysics , 1975, 40(2): 344-352. DOI:10.1190/1.1440530 |
| [19] | SEG Tech. Std. Comm.. Digital tape standards-SEG-A, SEG-B, SEG-C, SEG-Y, and SEG-D formats, plus SEG-D revisions 1 and 2. Soc. Expl. Geophys., 1977 |
| [20] | SEG Tech. Std. Comm.. SEG Y rev 1 Data Exchange format. Soc. Expl. Geophys., release 1.0, May 2002 |