SEGY文件读取与IBM浮点转IEEE的完整实现指南

简介SEGY是地震数据交换的通用格式在油气勘探与地质工程中应用广泛。这份资源专注于IBM32浮点编码的SEGY文件读写面向需要处理地震道数据、进行格式转换或二次开发的地质与物探技术人员。压缩包共3个文件包含2个cpp源文件和1个头文件其中头文件负责类接口声明两个源文件分别承担命令行交互与SEGY解析、IBM32解码等核心实现整体结构紧凑便于学习或集成。资源包仅5KB代码量不大适合有一定C基础、希望快速理解SEGY底层结构或提取可用工具的开发者。已有337人学习使用者可从中获得完整的最小实现示例、类封装思路以及针对IBM32格式的转换算法参考便于后续的数据预处理、分析或与其他系统交换对于正在学习地学数据处理的学生或工程师这份代码也能提供清晰的入门参考。1. 拿到 SEGY 文件却读不出来时问题多半不在文件而在字节序SEGY 是地震数据行业里最顽固的格式之一不是因为它复杂而是因为它太老了。1960 年代定下的布局后来在地球物理软件栈里层层叠加导致今天你手里一个.sgy或.segy文件可能是大端写的可能是小端写的可能是 IBM 浮点可能是 IEEE 浮点甚至采样点数与道头里写的都对不上。read_segy-master_readsegyibm32_SEGY文件读写_segy数据读取这类命名通常在勘探地球物理的脚本仓库里出现核心拆开就是两件事按 SEG-Y rev1 的二进制布局去读文件以及把 IBM 32 位浮点格式正确转换成 IEEE 754。对做地震数据处理、测井数据整理或搞 AWR2243 这类雷达数据落地的工程师来说这一块的坑不在读文件而在“读对了但数值全错”。这篇文章会从字节布局讲到可复现的读取代码再落到能验证数据正确性的方法上给出一套可以直接抄走的最小实现。2. SEG-Y 格式的二进制地图3600 字节文本头、400 字节二进制头、240 字节道头2.1 为什么读 SEGY 要先读二进制头而不是文本头SEG-Y rev1 的文件结构是严格的顺序布局文件开头是 3200 字节的 EBCDIC 文本头接着是 400 字节的二进制头之后才是连续排列的数据道。很多初学者上来就去解析文本头里的字符想从中拿到采样率、道数、格式码这些信息这是个常见的误用。文本头是给人看的里面的内容各家公司写得不一致有的写 ASCII 有的写 EBCDIC有的干脆全部填空格。真正可靠的信息全在 3201 到 3600 字节这 400 字节的二进制头里而且全部是以大端序存储的 2 字节或 4 字节整数。二进制头里最关键的三个字段是道数3201 字节处的 2 字节整数、每道采样点数3221 字节处的 2 字节整数和采样间隔3217 字节处的 2 字节整数单位微秒。这三个值决定了整个文件的几何结构知道了每道采样点数和采样间隔才能算出每道的字节长度知道了道数才能确定数据体的结束位置。第三个关键字段是数据格式码位于 3217 字节偏移处从 1 到 8 分别代表 IBM 浮点、4 字节整数、2 字节整数、定点数、IEEE 浮点等。read_segy这类工具通常默认处理格式码 1也就是标题里readsegyibm32所指的 IBM 32 位浮点。2.2 EBCDIC 文本头的读取必须做转码而不能硬解码3200 字节的文本头在标准定义里是 EBCDIC 编码的行文本每行 80 字节共 40 行。用 Python 直接按 ASCII 读出来是一堆乱码但这不代表文件损坏。正确处理方式是用codecs模块把这段字节从 EBCDIC 转成 ASCII 再打印import codecs def read_text_header(segy_bytes): raw segy_bytes[0:3200] if bC in raw or b in raw and all(b 128 for b in raw): return raw.decode(ascii, errorsreplace) return codecs.decode(raw, cp037, errorsreplace)这段代码的逻辑是先判断文本头里是否含有C字节字符这是很多国产采集软件写 ASCII 文本头时留下的痕迹。如果整段都能按 ASCII 解释就直接解码否则回退到cp037也就是 EBCDIC 的常见代码页。参数errorsreplace很关键因为某些厂商会把行尾补白字节写成二进制非字形字符严格解码会直接抛异常中断整个读取流程。实际项目里这个函数返回的文本基本不会用于后续计算它只用来查野外记录班报所以解码失败时宁可替换字符也不要抛出异常。2.3 道头的 240 字节里值得读的只有几个偏移每道数据由 240 字节道头加采样数据组成。道头里存了大量冗余信息和厂商私有字段但真正在数据处理链路里会用到的就集中在几个偏移位置。第一道的道序号在字节偏移 1 到 4CDP 道号在字节偏移 5 到 8每道采样点数在偏移 17 到 20采样间隔在偏移 49 到 52CMP 点的 X 坐标和 Y 坐标分别在偏移 73 到 76 和 77 到 80。如果是三维数据体inline 号和 crossline 号常见放在偏移 189 到 192 和 193 到 196但这不是标准强制规定很多处理系统会把这个位置拿来放其他内容。读道头时最需要注意的坑是道头本身的字节序和大端小端以及“道头里写的采样点数”和“二进制头里写的采样点数”不一致。我一般以二进制头里的 3221 字节处采样点数为准来分配数组长度但拿到每道之后会校验道头 17 字节处的值两者不一致时打印警告而不是直接崩溃。import struct def read_trace_header(segy_bytes, trace_offset): hdr segy_bytes[trace_offset:trace_offset 240] fields struct.unpack(IIIIII, hdr[0:24]) trace_seq fields[0] cdp fields[1] samples fields[2] 0xFFFF interval fields[3] 0xFFFF x struct.unpack(i, hdr[72:76])[0] y struct.unpack(i, hdr[76:80])[0] return { trace_seq: trace_seq, cdp: cdp, samples: samples, interval_us: interval, x: x, y: y, }注意这是按 SEG-Y rev1 的典型偏移来读的不同采集系统改过道头布局的并不少见生产环境里最好先打印一道道头全字段人工核对后再信任这个映射。这段代码里IIIIII的意义是六个大端序无符号 4 字节整数一次取出道头前 24 字节。采样点数和采样间隔在标准文档里是 2 字节字段但因为很多采集系统以 2 字节对齐写进了 4 字节槽位直接按I解包再与0xFFFF取与能同时兼容两种写法。X 和 Y 坐标用i是因为有的坐标系里存在负坐标而无符号解包会把负数读成 42 亿级别的大数直接污染后续几何计算。3. IBM 32 位浮点的底细与 IEEE 754 的换算readsegyibm32 的核心3.1 IBM 浮点不是“大端浮点”这么简单IBM System/360 时代的 32 位浮点结构是最高 1 位符号位接着 7 位二进制指数和 24 位尾数但指数基数是 16。同样的位模式IBM 浮点能表示的数值范围与 IEEE 754 有很大差别。最直观的表现是一个在 IEEE 下是0x3F800000的 1.0放到 IBM 浮点解码方式下会变成一个天文数字或者接近零的数。很多人在读 SEGY 时发现“读出来全是 1e-30 量级的小数”或者“全是 NaN”就是因为格式码写了 1但代码里用了 IEEE 解码。IBM 浮点数值的数学表达式为value (-1)^sign * (mantissa / 2^24) * 16^(exponent - 64)其中指数是 7 位无符号数偏移量为 64。与 IEEE 的不同点在于IEEE 的尾数是规格化的1.xxx形式而 IBM 浮点尾数是纯小数而且归一化的间隔是 2 的 4 次方也就是每 4 个二进制位为一组。这意味着 IBM 浮点没有 IEEE 那样的隐式前导 1直接位运算转换时不能用 IEEE 的那套(1 23) | mantissa操作。3.2 从 IBM 位模式到 IEEE 的向量化转换实现最直接的转换方式是逐字节提取位段后用浮点运算计算数值这在单道上没问题但如果一个三维数据体有几十万道每道 2000 个采样点那就是上亿次浮点运算Python 逐位循环几乎跑不出结果。所以生产代码里的常见做法是先把整个数据体的字节流一次性按无符号 4 字节整数读入 NumPy 数组然后对整个数组做向量化的位运算和浮点转换import numpy as np def ibm2ieee(raw_u32): raw_u32 raw_u32.astype(np.uint32, copyFalse) sign (raw_u32 31) 0x01 exp ((raw_u32 24) 0x7F) - 64 mant raw_u32 0x00FFFFFF out np.zeros_like(mant, dtypenp.float64) nz mant ! 0 out[nz] mant[nz].astype(np.float64) * (16.0 ** exp[nz].astype(np.float64)) out[sign 1] * -1.0 return out.astype(np.float32)这段代码里最关键的设计是将尾数为零的采样点单独屏蔽掉原因是0 * 16**exp在 exp 为负数时会出现 0 乘无穷大的情况NumPy 会给出 NaN 而不是 0。IBM 浮点里全零位模式表示数值 0这在野外采集的死道数据里非常常见必须保留为 0。astype(np.float32)放在最后一步是为了让返回值与其他处理模块兼容但在这一步之前用float64做中间计算能避免大尾数乘大指数时损失精度。3.3 反向转换从 IEEE 写回 IBM 浮点读写是对称的输出 SEGY 文件时同样要把 IEEE 浮点转回 IBM 格式。直接按数学公式回推是可行的但要注意规格化的处理。IBM 浮点的尾数左移应该以 4 位为一档因为指数基数是 16。反向转换的正确思路是先把值的绝对值取出来不断将尾数左移直到最高 4 位非零同时减少指数这样能保证尾数的有效位数尽量多def ieee2ibm(value): fv np.float64(value) sign 0 if fv 0 else 1 fv abs(fv) if fv 0.0: return np.uint32(0) exp 64 mant fv while mant 1.0 / 16.0: mant * 16.0 exp - 1 while mant 1.0: mant / 16.0 exp 1 mant_int int(round(mant * (1 24))) if mant_int (1 24): mant_int 4 exp 1 raw (sign 31) | ((exp 0x7F) 24) | mant_int return np.uint32(raw)这段代码里的第一个while解决的是绝对值小于 1/16 时的归一化问题第二个while解决的是 1 以上的情况。mant_int (1 24)这个保护分支很重要因为四舍五入可能让尾数溢出 24 位此时必须右移 4 位并加 1 指数否则写出来的位模式会被解码成完全错误的数值。位段长度含义说明bit 311符号位0 为正1 为负bit 30-247指数无符号偏移 64基数为 16bit 23-024尾数小数部分无隐式前导位4. 用 NumPy 把整套 SEGY 数据体读进内存并批量解析道头4.1 大文件读取的 IO 策略mmap 还是整读SEGY 文件动辄几个 GB整文件读入内存不是好选择但按照道逐条seek读取又太慢。一个折中的常见做法是先用os.path.getsize读出文件总字节数再用二进制头里的道数和采样点数算出理论字节数两者对比后决定走哪条路径。如果理论字节数和实际文件大小一致说明文件里没有尾部填充直接用numpy.fromfile把整个数据体读出来然后一次性 reshape如果两者不一致就要逐道读取并跳过可能的填充字节。numpy.fromfile在大文件上的性能远好于struct.unpack逐道循环因为前者把磁盘到内存的数据搬运完全交给了 C 层实现。配合np.memmap做内存映射可以避免占用过多物理内存但对大部分工作站来说一架次的三维地震数据通常不到 8GB直接读进内存做向量化处理的速度反而比 mmap 的缺页中断开销更优。4.2 最小可用的 read_segy 函数骨架下面这个函数可以作为一个独立模块的基础。它不做任何第三方 SEGY 库的依赖只依赖 NumPy输入是文件路径输出是道头字典列表和形状为(道数, 采样点数)的浮点数组。格式码为 2 或 8 的整数型 SEGY 文件也可以在这个骨架基础上加解析分支import os import numpy as np FORMAT_CODE_TO_DTYPE { 1: ibm, 2: i4, 3: i2, 5: f4, 8: i1, } def read_segy(path): with open(path, rb) as f: f.seek(3200) bh f.read(400) nalines int.from_bytes(bh[1:3], big) ns int.from_bytes(bh[21:23], big) dt_us int.from_bytes(bh[17:19], big) fmt int.from_bytes(bh[25:27], big) trace_size 240 ns * 4 if fmt in (1, 2, 5) else 240 ns * 2 file_size os.path.getsize(path) body_size file_size - 3600 ntraces body_size // trace_size if fmt 1: raw np.fromfile(path, dtypenp.uint32, offset3600) body ibm2ieee(raw[:ntraces * ns]).reshape(ntraces, ns) else: dt np.dtype(FORMAT_CODE_TO_DTYPE[fmt]) body np.fromfile(path, dtypedt, offset3600, countntraces * ns) body body.reshape(ntraces, ns) headers [] with open(path, rb) as f: for i in range(ntraces): f.seek(3600 i * trace_size) headers.append(read_trace_header(f.read(240), 0)) return headers, body, { ntraces: ntraces, ns: ns, dt_us: dt_us, format_code: fmt, }这里有个细节值得展开说明int.from_bytes(bh[21:23], big)读取的是 3221 字节处的两个字节对应每道采样点数但有些文件在这两个字节前还有 4 字节的扩展字段导致读取位置正确但数值看起来异常。如果打印出来的采样点数是 0 或超过 100000问题通常出在文件头不是标准 3600 字节而是厂商加了自定义数据段。这种情况的处理方式是先扫描文件里从 3600 字节开始的第一个非零值区域判断实际数据体的起始位置。4.3 参数对照格式码与每道字节数的关系格式码决定了数据体的每道字节数也决定了np.fromfile的 dtype 参数。上表可以写进代码注释里备用实际上手时最常遇到的是 1IBM 浮点和 5IEEE 浮点。格式码 5 的文件读取路径最简单因为f4就是大端 IEEE 单精度NumPy 原生支持不需要任何转换函数。格式码 2 是 4 字节整数常见于早期系统读取后要乘以一个标度值才能得到真实物理量标度值存在道头偏移 69 到 70 字节处。格式码含义每采样点字节数读取方式1IBM 32 位浮点4按 uint32 读入后调用 ibm2ieee24 字节整数4大端有符号整数需结合道头标度值32 字节整数2大端有符号短整型5IEEE 32 位浮点4直接f4读取81 字节整数1大端有符号字节注意以上字节数不包含 240 字节道头计算总文件大小时必须加上ntraces * 240。5. 验证读出来的数据是否可信对拍、往返转换和统计量检查数据读出来后最忌讳的事情是“读出来是一堆数就觉得成功了”。SEGY 数据读取的错误往往不是抛异常而是静默的数值错误。我常用的验证手段有三种按成本从低到高排列能快速定位问题出在格式码判断、字节序转换还是道头偏移解析上。第一是往返转换验证。把读出来的 IEEE 数组逐块转成 IBM 位模式再调用ibm2ieee转回对比原始字节流是否一致。这个过程不依赖任何外部工具纯粹验证转换函数的自洽性。需要注意浮点转换的舍入误差所以比较时用数值容差而不是严格相等def verify_roundtrip(original_bytes, ieee_array): back ieee2ibm(original_bytes[:len(ieee_array)]) roundtrip ibm2ieee(back.astype(np.uint32)) max_err np.abs(roundtrip - ieee_array).max() rel_err max_err / np.abs(ieee_array).max() return rel_err相对误差在 1e-6 以内说明转换函数内部没有符号位或指数偏移的错误。如果误差达到 0.1 量级最常见的原因是转换成 IEEE 后没有做astype(np.float32)导致后续比较是在 float64 和 float32 之间进行。第二种验证是与工具箱对拍。如果系统里有 Seismic Unix 或者 ObsPy 的segy模块可以用它们读同一个文件对比道头和振幅统计量。第三方库同样会踩到字节序和 IBM 浮点的坑所以两边对不上时未必是自己的代码错了还要看对方的版本是否修补了非标准头的兼容问题。对比时重点看三个统计量最小值、最大值和能量np.sqrt(np.mean(x**2))。第三种验证更直接是针对未爆数据的物理合理性检查。地震数据体通常是零均值的而且每道的绝对能量在一个可预期的范围内。如果读出来的数据全部是整段的常量或者振幅超过 1e20那几乎可以断定格式码判断错了比如实际文件是 IEEE 但头里写了 1、又或者数据是按道交错存储而非连续存储。这种情况下回头检查二进制头里的格式码字节同时用xxd打印一段数据体的十六进制观察是否为可识别的浮点位模式。最后一个实用技巧是把每个文件的读取参数做成缓存字典以文件的采样点数、道数和格式码为键。因为处理同一工区的数据时这些参数通常完全一致第一次读取后缓存元数据后续文件就跳过头部解析直接读数据体对批处理几百个 SEGY 文件的场景能省下不少 IO 时间。本文还有配套的精品资源点击获取