SEGY是地震勘探存储数据的标准格式,最早在1975年发布revision 0。由于那时候地震数据的处理一般是在服务器上进行的,所以SEGY格式打上了时代印记,例如SEGY文件采用大字节序,数据点是IBM 32位浮点数。2002年发布的版本revision 1增加了IEEE 32位浮点数,但还是采用大字节序。
PC的处理器都是小字节序的,32位浮点数是IEEE浮点数,因此在PC上读取SEGY格式的文件时就比较麻烦,主要体现在:
- 读取SEGY文件里的整数或IEEE浮点数之后,必须进行字节序交换,转换为小字节序,才能正确地表示为PC上的数。
32位的IBM浮点数与32位IEEE浮点数格式不一样,不是简单的交换字节序的问题,需要经过计算才能将读取的IBM浮点数转换为PC上的浮点数。
本人在编写的一个数据处理软件中,需要读取其他采集设备转换出的标准SEGY文件,有的SEGY文件采用IBM浮点数,就涉及到将IBM浮点数转换为IEEE浮点数值。
关于32位IBM和IEEE浮点数结构及其转换方法可参考下面的一个文献:
符茂松《32位IEEE和IBM浮点数结构及其转换方法》,工程地球物理学报,2011年12月。
图A. 32位IBM浮点数结构
32位IBM浮点数的结构如图A所示,其具体解释见参考文献。
下面是读取一个4字节IBM浮点数,然后将其转换为IEEE浮点数的示意代码(只是从文件读取IBM浮点数、字节序交换、转换为IEEE浮点数数值的示意代码,不是读取一个文件内容的完整代码)。
QFile aFile(aFileName); //以文件方式读出
if (!(aFile.open(QIODevice::ReadOnly)))
return false;
unsigned long IBMFloatBytes; //即quint32
aFile.read((char *)&IBMFloatBytes, 4);//用无符号4字节整数读出IBM 浮点数的4个字节
// 交换字节,变成如图A的字节序
IBMFloatBytes=(((IBMFloatBytes>>24)&0xff) | ((IBMFloatBytes&0xff)<<24) |
((IBMFloatBytes>>8)&0xff00) | ((IBMFloatBytes&0xff00)<<8));
float aSampValue=IBM2IEEE_float(IBMFloatBytes);//4字节IBM浮点数的实际数值
- 从文件中读取4字节IBM浮点数时,采用是读取一个4字节无符号整数的方法,这样32位数据原始不动的保存到变量IBMFloatBytes里。
变量IBMFloatBytes的4个字节实际是大字节序的IBM浮点数,将其进行字节序交换,变成小字节序,也就是图A所示的顺序。
- 调用函数IBM2IEEE_float()将变量IBMFloatBytes存储的小字节序的4个字节转换为IEEE浮点数。PC上C语言里的float就是4字节IEEE浮点数。
下面是函数IBM2IEEE_float( ) 的完整代码,它的功能是返回如图A存储的4个字节表示的浮点数的数值。
float QsgyFileReader::IBM2IEEE_float(const unsigned long IBMFloatBytes)
{ //IBM 浮点数转换为IEEE浮点数, IBMFloatBytes是已经转换为小字节序的 IEEE浮点数32位编码
unsigned long ibmCode=IBMFloatBytes;
float ieeeFloat=0;
//IBM浮点数: SEEEEEEE MMMMMMMM MMMMMMMM MMMMMMMM
//浮点数值 Value = (-1)^s * M * 16^(E-64)
unsigned long signCode=(ibmCode>>31);//获取符号位, S=0或1, sign=00 00 00 0000000S
int sign=1; //正数
if (signCode==1)
sign=-1;//负数
ibmCode= (ibmCode<<1); // 左移移出符号位, ibmCode= EEEEEEEM MMMMMMMM MMMMMMMM MMMMMMM0
int exponent = (int)(ibmCode >> 25); // 获取阶数, exponent=00 00 00 0EEEEEEE
unsigned long fraction= (ibmCode << 7); //移出符号位和阶数剩余的部分:尾数部分fraction=MMMMMMMM MMMMMMMM MMMMMMM 00000000
fraction= (fraction >> 8); //fraction=00000000 MMMMMMMM MMMMMMMM MMMMMMM
static double ratio=qPow(2,24); //使用静态局部变量,避免重复计算
if ((exponent==0) && (fraction==0)) //00000000 00 00 00 或10000000 00 00 00 都表示0
ieeeFloat=0;
else //IBMfloat= (sign)*Fraction*16^(exponent-64);
{
double P2=qPow(16, exponent-64);
double P1=fraction/ratio;
ieeeFloat=sign*P1*P2;
}
return ieeeFloat;
}
IBM2IEEE_float( )函数的代码通过多次的移位操作,得到浮点数的符号、指数部分和尾数部分,然后计算出其表示的浮点数。代码不详细解释了,看注释即懂。
这里只实现了读取IBM浮点数,转换为的IEEE浮点数是什么样的存储结构由CPU自动处理。在导出标准格式SEGY文件时,使用IEEE浮点数格式,避免进行IEEE到IBM格式的转换。
图B是本人编写的数据处理软件中的数据波形显示,读取的数据文件是数据点格式为32位IBM浮点数的标准SEGY文件。
图B 读取数据点格式为IBM浮点数的标准SEGY文件并显示其波形