简介SGY是地震勘探领域常用的数据存储格式其二进制结构包含文件头、道头与样本数据面向这一格式的读取脚本专注于解析这些部分适合物探工程师、地学研究人员以及需要处理相关数据的学生和技术人员使用。压缩包为RAR格式其中仅有一个脚本文件整体体积约1KB代码虽小却覆盖了从打开文件、读取文件头与道头信息、按采样间隔提取地震道数据到转化为数组并安全关闭文件的全过程。脚本还允许使用者根据实际需求调整参数例如指定路径、选择特定道或时间区间并可修改数据解释方式以匹配浮点数或整数样本这些细节充分体现了SGY标准中元数据管理的特点。目前已有七百九十三人学习。对于希望掌握SGY格式要点或快速搭建数据预处理工具的读者这份代码提供了一个直观模板能帮助把二进制地震数据转换为易于分析的矩阵进而支持石油勘探、地质灾害评估等实际应用。1. readsegy是谁把SEG-Y二进制格式变成可用的地震道数据拿到一块野外采集的SEG-Y地震数据几百兆到几十GB的二进制文件文本编辑器打开全是乱码Excel更是直接罢工——readsegy就是专为读这种SGY格式数据而生的函数。地震勘探行业里SEG-Y就是事实标准从野外仪器出来的原始炮集到处理中心交出去的偏移剖面几乎全是这个格式。readsegy干的事情很纯粹把文件里按字节排布的道头信息和地震道振幅值读出来组织成Python能直接用的数组和字典。它的价值在于你不用再去手动翻SEG-Y标准文档、对着字节偏移表算地址也不用纠结IBM浮点转IEEE这种细节一个函数调用就能拿到带坐标、带采样率的地震道数据。适合谁地震数据处理工程师、勘探开发的研究人员、还有刚接触地震数据的学生——凡是拿到SGY文件需要做振幅分析、画剖面或者喂给深度学习模型的人都绕不开这一步。2. 先看懂SEG-Y结构卷头、二进制头、道头与数据体的字节分布想用好readsegy不能只把它当黑匣子。SEG-Y格式的读取逻辑完全建立在字节偏移上不搞清楚这些偏移遇到异常数据你连报错信息都读不懂。这一章把格式拆开揉碎后面调参数才有根据。2.1 卷头与二进制头3600字节里藏着什么一个标准的SEG-Y文件开头3600字节是卷头。前3200字节是文本卷头Textual File Header用来存ASCII或EBCDIC编码的文本信息比如工区名、采集日期、处理流程说明。老资料用EBCDIC编码居多新采集的数据一般用ASCII。readsegy在读取时会把这两类编码都解出来解不出来的字节用空格兜底不会因为个别字节非法就让整个读取崩溃。紧跟着的400字节是二进制卷头Binary File Header里面全是定长的二进制字段每个字段占据的字节位置由SEG-Y标准规定死。读这个卷头有个常见误区好多人以为二进制卷头就是一堆没用的配置信息跳过就行。其实不是。字节3225到3226存的是采样间隔单位微秒字节3221到3222存的是每道采样点数字节3217到3218存的是数据格式码——这个格式码决定了后面每个采样值用几字节、什么浮点表示法。readsegy内部就是先读这400个字节把采样间隔、采样点数、数据格式解出来再决定后续按什么规则解释数据体。代码里会写import struct def read_binary_header(f): f.seek(3200) raw f.read(400) # 大端字节序SEG-Y标准规定高位在前 sample_interval struct.unpack(H, raw[5:7])[0] # 字节3225-3226 samples_per_trace struct.unpack(H, raw[9:11])[0] # 字节3221-3222 format_code struct.unpack(H, raw[15:17])[0] # 字节3217-3218 return { sample_interval: sample_interval, samples_per_trace: samples_per_trace, format_code: format_code }这段代码的逻辑很直接先用seek定位到二进制卷头起始位置读400字节再用struct.unpack按大端序抽取关键字段。为什么用H而不是H因为SEG-Y标准规定文件里的多字节整数用大端序高位字节在前这和x86机器的内存序相反。不过现实中有少数软件写出来的文件用的是小端序这就是后话readsegy里会做一个字节序探测来解决。2.2 道头每道数据的第一张身份证卷头往后文件主体就是一长串地震道。每一道由240字节的道头Trace Header和紧随其后的数据体组成。道头里记录的字段分布在240字节的不同偏移上readsegy最关心的几个是字节偏移长度字段含义1-44字节道序列号5-84字节道号从1开始9-124字节原始道号17-204字节采样点数此道实际采样数21-244字节采样间隔73-764字节道所属测线号LINE_NUMBER77-804字节道所属道集号CDP_NUMBER81-844字节道集内道号85-884字节炮点号SOURCE_NUMBER89-924字节道集内炮点号115-1162字节坐标单位1米2英尺117-1204字节炮点X坐标121-1244字节炮点Y坐标181-1844字节接收点X坐标185-1884字节接收点Y坐标这里的字节偏移是SEG-Y标准的“1-based”描述法和代码里的索引还差一位写代码时务必注意。readsegy实现时会从这个240字节的道头里把坐标、测线号、道集号提取出来放进返回的headers列表里。很多用户拿readsegy读数据真正要用的就是这个道头里的坐标信息——没有坐标地震道数据就是一堆波形数字没法定位到地下位置。2.3 数据体格式IBM浮点与IEEE浮点的坑道头后面是数据体。这里的格式不能想当然。SEG-Y标准里定义了多种数据格式码最常见的是格式码1IBM 32位浮点、格式码232位有符号整数、格式码5IEEE 32位浮点和格式码864位IEEE浮点。国内老资料大多是格式码1新处理成果多是格式码5。readsegy读取时会对格式码做判断格式码1先按4字节IBM浮点解释再把IBM转成IEEE格式码5直接按IEEE解。IBM浮点转IEEE是这里最烧脑的一步。IBM浮点的内部结构是1位符号位、7位指数位基数为16、24位尾数。转换公式是value mantissa * 16^(exponent-64)。网上一搜能搜到一堆转换代码但抄错的不少。readsegy的转换逻辑里有一个关键处理尾数要按有符号数理解否则负数的解析会翻车。另外IBM浮点有“非规格化数”存在不加判断直接套公式会得到垃圾值。这个细节后面避坑章展开。def ibm_to_ieee(raw_ibm): # raw_ibm 是 4 字节的 int 表示 sign (raw_ibm 31) 0x01 exponent (raw_ibm 24) 0x7F mantissa raw_ibm 0x00FFFFFF if mantissa 0: return 0.0 # 尾数最高位是符号位需要转成有符号数 if mantissa 0x800000: mantissa - 0x1000000 value mantissa * (16.0 ** (exponent - 64)) return -value if sign else value这段转换代码注意两点一是mantissa要按24位有符号数处理不然负尾数会被当成很大的正数二是指数基数是16不是216.0 ** (exponent - 64)这步错了结果就是天文数字。当年我第一次写转换脚本把基数误写成2读出来的振幅值比正常大了几十个数量级画出来的剖面满屏都是蓝红色块排查了整整一个下午。3. 用readsegy在本地跑通最小读取从文件句柄到地震道数组格式结构清楚了代码就可以上桌。这一章把readsegy的完整实现写出来从打开文件到输出三维数组每行代码都能直接复制运行。这里给出的是一个自包含版本不依赖任何第三方库只用了Python标准库的struct和os。3.1 最小读取代码readsegy的核心骨架import struct import os def readsegy(filename, endian): 读取 SEG-Y 文件返回 text, bin_headers, data data: numpy 数组形状 (道数, 每道采样点数) f open(filename, rb) # 1. 读文本卷头 text_raw f.read(3200) try: text text_raw.decode(ascii, errorsreplace) except UnicodeDecodeError: text text_raw.decode(cp037, errorsreplace) # EBCDIC # 2. 读二进制卷头 bin_raw f.read(400) sample_interval struct.unpack(endian H, bin_raw[5:7])[0] samples_per_trace struct.unpack(endian H, bin_raw[9:11])[0] format_code struct.unpack(endian H, bin_raw[15:17])[0] # 3. 循环读道 traces [] headers [] while True: head_raw f.read(240) if len(head_raw) 240: break # 道头关键字段 trace_seq struct.unpack(endian I, head_raw[0:4])[0] line_num struct.unpack(endian i, head_raw[72:76])[0] cdp_num struct.unpack(endian i, head_raw[76:80])[0] sx struct.unpack(endian i, head_raw[116:120])[0] sy struct.unpack(endian i, head_raw[120:124])[0] rx struct.unpack(endian i, head_raw[180:184])[0] ry struct.unpack(endian i, head_raw[184:188])[0] coord_unit struct.unpack(endian H, head_raw[114:116])[0] headers.append({ trace_seq: trace_seq, line_num: line_num, cdp_num: cdp_num, sx: sx, sy: sy, rx: rx, ry: ry, coord_unit: coord_unit }) # 读一道数据 trace_raw f.read(4 * samples_per_trace) if len(trace_raw) 4 * samples_per_trace: break # 根据格式码解析 if format_code 1: vals struct.unpack(endian I * samples_per_trace, trace_raw) trace [ibm_to_ieee(v) for v in vals] elif format_code 5: trace list(struct.unpack(endian f * samples_per_trace, trace_raw)) elif format_code 2: vals struct.unpack(endian i * samples_per_trace, trace_raw) trace [v / 32768.0 for v in vals] # 归一化到 -1~1 else: raise ValueError(fUnsupported format code: {format_code}) traces.append(trace) f.close() return text, headers, traces逻辑说明第一步读3200字节文本头先按ASCII解码失败就按EBCDIC代码页037兜底。第二步用上一章讲的关键偏移解析采样间隔、采样点数和格式码。第三步进入主循环每道先读240字节道头按偏移列表抽坐标字段再根据format_code决定数据体解译方式。返回的三元组text, headers, traces中text是工区文字信息headers是道头字典列表traces是二维列表每一行是一道采样值序列。3.2 参数怎么设endian、坐标单位与格式码的取舍readsegy在真实使用中有四个参数几乎必调不调或者调错结果就是数据读出来但坐标对不上、振幅对不上。第一个是字节序endian。标准SEG-Y是大端但1990年代有些国产软件写文件时直接用了本机的小端序这个只有在读取时才发现。判断方法简单粗暴先按大端读读出来的采样间隔如果是一个几十万的不正常数字或者道头坐标全是负数就换成小端再试。readsegy做成参数就是给这个留后路。第二个是坐标单位。道头里115-116字节有坐标单位标记1表示米2表示英尺。但很多老数据这个字段根本没填填0。遇到这种情况不能默默认作米要结合工区常识判断——一条测线跨度几千个单位的大概率是英尺。readsegy不帮你换算只是把这个原始值原样透传换算逻辑留给使用者。第三个是format_code。前面代码里默认支持格式码1、2、5但还有一些数据用格式码316位定点数和格式码864位浮点。如果公司内部的老软件处理出了格式码3你需要自己扩展解析分支——格式码3按h解出来再除以32。实际处理时遇到不支持的格式码先别急着改代码用十六进制编辑器看一眼文件头确认格式码再决定怎么扩。第四个是采样点数不一致的问题。道头里有“本道采样点数”字段字节17-20和二进制卷头里的“每道采样点数”可能不一致。最后几道有时会短一截。readsegy的实现策略是优先信道头里的采样点数如果某道读出来的字节数不够按实际长度截断。这个细节在后面的避坑章细说。3.3 从readsegy到地震道数组数据如何组织实际处理时极少有人直接拿readsegy返回的list套list去干活因为list的索引访问慢切片也麻烦。常见做法是把它转成numpy数组import numpy as np text, headers, traces readsegy(survey.sgy) data np.array(traces, dtypenp.float32) print(data.shape) # (道数, 每道采样点数)转成np.float32数组后后续的FFT、滤波、增益、画剖面都可以直接用向量化操作。这时候readsegy返回的headers列表要和data的每一行一一对应——第i道的坐标就是headers[i]。这个对应关系看似天然但实际操作中如果中间做了道编辑比如删除了坏道headers和数据行就会错位所以我在工程里会把headers直接转成结构化数组和数据一起存储import numpy as np header_dtype np.dtype([ (trace_seq, i4), (line_num, i4), (cdp_num, i4), (sx, i4), (sy, i4), (rx, i4), (ry, i4), (coord_unit, i2) ]) header_arr np.array([tuple(h.values()) for h in headers], dtypeheader_dtype)这样headers_arr和data在行维度上严格对应之后做CDP分选或者inline/xline重排时挪数据行必须同步挪headers行一个切片操作就能同时搞定。4. 常见做法对比segyio、obspy、segpy与readsegy的取舍readsegy不是市面上唯一的解法。地震数据处理领域还有segyio、obspy、segpy这些常用工具各有各的适用场景。这一章做一个横向对比帮你在实际项目中选型。4.1 四个库的能力边界对比库底层实现读取速度支持格式码适用场景readsegy自实现纯Python struct中慢取决于自定义教学、快速探索、格式码特殊的私有数据segyioC扩展快1、2、5、8大规模生产、需要随机访问道数据obspyPython numpy中1、2、5地震学波形处理、与地震目录联动segpy纯Python构造和解析慢1、2、5写SEG-Y文件、轻量解析速度差距在文件尺寸上体现非常明显。一个三维工区可能有几千条测线、每条测线上百万个地震道单个文件几十GB。纯Python写循环逐道读取一小时能读完就算不错segyio用C扩展加内存映射十分钟内存就能把道头索引建完。所以生产环境我几乎不用自实现的readsegy去读大文件它更多是用来快速确认文件格式、看道头信息、或者调试格式码解析逻辑。4.2 大文件场景readsegy的块读取与内存控制文件一大内存就成了瓶颈。一个地震道采样点数假设是2000按4字节浮点算一道占8KB一百万个地震道就是8GB。一次性全读进内存32GB的机器也能被撑到边缘。实际处理时我习惯按“块”来组织读取逻辑def readsegy_block(filename, start_trace, num_traces, endian): f open(filename, rb) f.seek(3600) f.seek(start_trace * (240 4 * n_samples), 1) block_head [] block_data [] for _ in range(num_traces): head_raw f.read(240) if len(head_raw) 240: break block_head.append(parse_trace_header(head_raw)) trace_raw f.read(4 * n_samples) block_data.append(decode_trace(trace_raw)) f.close() return block_head, block_data这里的关键是计算偏移位置每一道在文件里占据的字节数是 240道头 采样点数 × 每采样点字节数。要在文件里找第 start_trace 道直接用乘法跳到该位置就行不用从头顺序读。这在处理大型三维数据做inline切片时特别有用——先读道头索引定位到某条inline的所有道再只读那一段数据体内存占用可以压缩到原来的十分之一。4.3 什么时候别用readsegy坐标校正与道头重写场景readsegy适合读但遇到需要重写SGY文件的场景就力不从心了。比如野外采集的数据坐标写错了要做批量坐标偏移校正这时候不仅要读还要改道头、写文件、保持字节对齐。segyio在这方面有原生优势它提供按道号随机读写的接口改一个道的坐标字段不需要重写整个文件。我一般这样分工格式探索和快速绘图用readsegy正式的生产管线道头校正、重排、输出重新打包用segyio。另外obsPy是个地震学通用框架适合处理天然地震波形数据它读SEG-Y的能力是顺带支持的如果数据里要同时处理事件目录、波形拾取这些地震学概念才值得引入。5. 读取SEG-Y的五个真实踩坑字节序、IBM浮点、延长道与坐标单位这一章写的是我在实际项目中踩过的五个坑每一条都对应一个真实遇到过的数据集。写出来一是让你少花排查时间二是让你知道遇到这些异常时该去哪里看。5.1 IBM浮点转换导致振幅全乱现象读出的地震道数值范围在±1e20附近甚至出现NaN画出来的剖面全是噪点。 原因格式码为1的数据用的是IBM浮点但转换函数写错了。最常见的错误是把IBM的指数当成基2而不是基16或者尾数没按有符号数处理。还有个隐蔽情况某些老软件写IBM浮点时高低字节互换I解出来的整数整体对不上。 解决先取一道数据用struct.unpack(I, raw)得到整数再按上一章的ibm_to_ieee手动转换验证最大最小值是否在物理合理范围速度、振幅通常不超过几千到几万量级。如果转换后数值还是异常用十六进制编辑器直接打开那一道的原始字节肉眼比对IBM浮点的标准编码。5.2 大端小端之争采样间隔异常大现象读出来的采样间隔不是常见的2000微秒2毫秒或4000微秒4毫秒而是几十万微秒。 原因这是典型的字节序问题。文件实际是小端序但readsegy默认按大端解。struct.unpack(H, bin_raw[5:7])把字节序反过来之后2000微秒的旧数据变成了一个天文数字。 解决把readsegy的endian参数改成重新读一遍。如果恢复正常说明这个文件确实是用小端序写的。这时要记录下这个参数每一道解析都要用相同字节序。5.3 延长道Extended Trace导致读取中断现象文件末尾多出一大段空白道readsegy读到一半报错或者读出来的数据最后几行全是0。 原因SEG-Y标准规定每道数据体长度等于 二进制卷头里的采样点数 × 每采样点字节数。但有些处理软件为了对齐到磁盘块在文件末尾追加了“延长道”Extended Trace长度不定里面填0。readsegy不知道文件里多了这玩意按正常道解析结果最后一条道长度对不上。 解决主循环里加一条条件判断如果剩余文件长度小于“道头240 道数据长度”就停止读取。另外查看二进制卷头的字节3217-3218里有个“扩展道数”字段有些软件会写这里可以用它来判断末尾是否有延长道。我的实现里选择直接忽略延长道——它们在处理中本来就没有物理意义。5.4 坐标单位不是1也不是2CDP坐标对不上位置现象读出来的道头坐标能和工区底图对上但缩放比例完全不对比如几百米的道间距变成了几千。 原因二进制卷头里坐标单位字段写的是0或空白。很多老软件写道头时忘了填坐标单位readsegy读出来就是0使用者默认当米处理实际上原数据可能是英尺。坐标数值往往大了3.28倍。 解决读取后检查coord_unit字段如果为0需要结合工区已知参数反推。一个土办法取一条测线连续10道的X坐标算平均道间距如果间距在20-50之间可能是米在60-170之间大概是英尺。确定后自己做一次换算再传给下游。5.5 文本头解码乱码EBCDIC和ASCII的混合体现象readsegy返回的text字段里一半是英文正常字符一半是奇怪的符号和空格。 原因文本卷头标准定义可以用ASCII或EBCDIC。有些写盘软件只在第1字节写了个标识符其余还是EBCDIC另一些则干脆混合使用。readsegy只做了一轮尝试性解码碰到混合编码就翻车。 解决先检查字节3200-3200文本卷头第一个字节看它是C1EBCDIC的A还是41ASCII的A。如果是C1整段用cp037解码如果是41用ascii。混合编码的文件解码后人工看一遍把乱码段手工替换。坦白说到了这一步已经是玄学范畴了我的经验是别在文本头上纠缠工区信息最后都会写进数据库文本卷头过一遍眼确认工区名没错就行。6. 验证readsegy读得对不对与segyio交叉核对和可视化抽查读完SGY数据之后最怕的不是报错而是数据安安静静地读出来了但结果不可信。我现在的固定动作是做交叉验证用readsegy和segyio分别读同一个文件对比道头坐标和振幅统计量。import segyio import numpy as np # 用 segyio 读同一文件做核对 with segyio.open(survey.sgy, ignore_geometryTrue) as segy_f: segy_trace segy_f.trace[0] segy_header segy_f.header[0] # 用 readsegy 读同样的一道 text, headers, traces readsegy(survey.sgy) my_trace np.array(traces[0]) my_header headers[0] # 核对道头坐标 print(segyio sx:, segy_header[segyio.su.sx], readsegy sx:, my_header[sx]) # 核对振幅数组 print(max abs diff:, np.max(np.abs(segy_trace - my_trace)))如果最大振幅差接近零说明两道解析路径对上了如果差得很大把格式码和字节序检查一遍。另一个验证方法是可视化抽查随机抽几道画成地震波形图和商业软件如Vista、Omega显示的剖面做肉眼对比。波形形状一致、振幅量级一致才敢继续往下走。一次在甲方项目里我用readsegy读了他们的老磁带转储数据故事都进行到振幅补偿了最后交叉验证时发现有一半文件的坐标单位是英尺。那次之后我学了个规矩每个SGY文件进处理流程之前先跑这个双读核对脚本不通过就不进下一步。那个核对脚本至今还在我工具目录里躺着。希望帮到你。本文还有配套的精品资源点击获取
简介SGY是地震勘探领域常用的数据存储格式其二进制结构包含文件头、道头与样本数据面向这一格式的读取脚本专注于解析这些部分适合物探工程师、地学研究人员以及需要处理相关数据的学生和技术人员使用。压缩包为RAR格式其中仅有一个脚本文件整体体积约1KB代码虽小却覆盖了从打开文件、读取文件头与道头信息、按采样间隔提取地震道数据到转化为数组并安全关闭文件的全过程。脚本还允许使用者根据实际需求调整参数例如指定路径、选择特定道或时间区间并可修改数据解释方式以匹配浮点数或整数样本这些细节充分体现了SGY标准中元数据管理的特点。目前已有七百九十三人学习。对于希望掌握SGY格式要点或快速搭建数据预处理工具的读者这份代码提供了一个直观模板能帮助把二进制地震数据转换为易于分析的矩阵进而支持石油勘探、地质灾害评估等实际应用。1. readsegy是谁把SEG-Y二进制格式变成可用的地震道数据拿到一块野外采集的SEG-Y地震数据几百兆到几十GB的二进制文件文本编辑器打开全是乱码Excel更是直接罢工——readsegy就是专为读这种SGY格式数据而生的函数。地震勘探行业里SEG-Y就是事实标准从野外仪器出来的原始炮集到处理中心交出去的偏移剖面几乎全是这个格式。readsegy干的事情很纯粹把文件里按字节排布的道头信息和地震道振幅值读出来组织成Python能直接用的数组和字典。它的价值在于你不用再去手动翻SEG-Y标准文档、对着字节偏移表算地址也不用纠结IBM浮点转IEEE这种细节一个函数调用就能拿到带坐标、带采样率的地震道数据。适合谁地震数据处理工程师、勘探开发的研究人员、还有刚接触地震数据的学生——凡是拿到SGY文件需要做振幅分析、画剖面或者喂给深度学习模型的人都绕不开这一步。2. 先看懂SEG-Y结构卷头、二进制头、道头与数据体的字节分布想用好readsegy不能只把它当黑匣子。SEG-Y格式的读取逻辑完全建立在字节偏移上不搞清楚这些偏移遇到异常数据你连报错信息都读不懂。这一章把格式拆开揉碎后面调参数才有根据。2.1 卷头与二进制头3600字节里藏着什么一个标准的SEG-Y文件开头3600字节是卷头。前3200字节是文本卷头Textual File Header用来存ASCII或EBCDIC编码的文本信息比如工区名、采集日期、处理流程说明。老资料用EBCDIC编码居多新采集的数据一般用ASCII。readsegy在读取时会把这两类编码都解出来解不出来的字节用空格兜底不会因为个别字节非法就让整个读取崩溃。紧跟着的400字节是二进制卷头Binary File Header里面全是定长的二进制字段每个字段占据的字节位置由SEG-Y标准规定死。读这个卷头有个常见误区好多人以为二进制卷头就是一堆没用的配置信息跳过就行。其实不是。字节3225到3226存的是采样间隔单位微秒字节3221到3222存的是每道采样点数字节3217到3218存的是数据格式码——这个格式码决定了后面每个采样值用几字节、什么浮点表示法。readsegy内部就是先读这400个字节把采样间隔、采样点数、数据格式解出来再决定后续按什么规则解释数据体。代码里会写import struct def read_binary_header(f): f.seek(3200) raw f.read(400) # 大端字节序SEG-Y标准规定高位在前 sample_interval struct.unpack(H, raw[5:7])[0] # 字节3225-3226 samples_per_trace struct.unpack(H, raw[9:11])[0] # 字节3221-3222 format_code struct.unpack(H, raw[15:17])[0] # 字节3217-3218 return { sample_interval: sample_interval, samples_per_trace: samples_per_trace, format_code: format_code }这段代码的逻辑很直接先用seek定位到二进制卷头起始位置读400字节再用struct.unpack按大端序抽取关键字段。为什么用H而不是H因为SEG-Y标准规定文件里的多字节整数用大端序高位字节在前这和x86机器的内存序相反。不过现实中有少数软件写出来的文件用的是小端序这就是后话readsegy里会做一个字节序探测来解决。2.2 道头每道数据的第一张身份证卷头往后文件主体就是一长串地震道。每一道由240字节的道头Trace Header和紧随其后的数据体组成。道头里记录的字段分布在240字节的不同偏移上readsegy最关心的几个是字节偏移长度字段含义1-44字节道序列号5-84字节道号从1开始9-124字节原始道号17-204字节采样点数此道实际采样数21-244字节采样间隔73-764字节道所属测线号LINE_NUMBER77-804字节道所属道集号CDP_NUMBER81-844字节道集内道号85-884字节炮点号SOURCE_NUMBER89-924字节道集内炮点号115-1162字节坐标单位1米2英尺117-1204字节炮点X坐标121-1244字节炮点Y坐标181-1844字节接收点X坐标185-1884字节接收点Y坐标这里的字节偏移是SEG-Y标准的“1-based”描述法和代码里的索引还差一位写代码时务必注意。readsegy实现时会从这个240字节的道头里把坐标、测线号、道集号提取出来放进返回的headers列表里。很多用户拿readsegy读数据真正要用的就是这个道头里的坐标信息——没有坐标地震道数据就是一堆波形数字没法定位到地下位置。2.3 数据体格式IBM浮点与IEEE浮点的坑道头后面是数据体。这里的格式不能想当然。SEG-Y标准里定义了多种数据格式码最常见的是格式码1IBM 32位浮点、格式码232位有符号整数、格式码5IEEE 32位浮点和格式码864位IEEE浮点。国内老资料大多是格式码1新处理成果多是格式码5。readsegy读取时会对格式码做判断格式码1先按4字节IBM浮点解释再把IBM转成IEEE格式码5直接按IEEE解。IBM浮点转IEEE是这里最烧脑的一步。IBM浮点的内部结构是1位符号位、7位指数位基数为16、24位尾数。转换公式是value mantissa * 16^(exponent-64)。网上一搜能搜到一堆转换代码但抄错的不少。readsegy的转换逻辑里有一个关键处理尾数要按有符号数理解否则负数的解析会翻车。另外IBM浮点有“非规格化数”存在不加判断直接套公式会得到垃圾值。这个细节后面避坑章展开。def ibm_to_ieee(raw_ibm): # raw_ibm 是 4 字节的 int 表示 sign (raw_ibm 31) 0x01 exponent (raw_ibm 24) 0x7F mantissa raw_ibm 0x00FFFFFF if mantissa 0: return 0.0 # 尾数最高位是符号位需要转成有符号数 if mantissa 0x800000: mantissa - 0x1000000 value mantissa * (16.0 ** (exponent - 64)) return -value if sign else value这段转换代码注意两点一是mantissa要按24位有符号数处理不然负尾数会被当成很大的正数二是指数基数是16不是216.0 ** (exponent - 64)这步错了结果就是天文数字。当年我第一次写转换脚本把基数误写成2读出来的振幅值比正常大了几十个数量级画出来的剖面满屏都是蓝红色块排查了整整一个下午。3. 用readsegy在本地跑通最小读取从文件句柄到地震道数组格式结构清楚了代码就可以上桌。这一章把readsegy的完整实现写出来从打开文件到输出三维数组每行代码都能直接复制运行。这里给出的是一个自包含版本不依赖任何第三方库只用了Python标准库的struct和os。3.1 最小读取代码readsegy的核心骨架import struct import os def readsegy(filename, endian): 读取 SEG-Y 文件返回 text, bin_headers, data data: numpy 数组形状 (道数, 每道采样点数) f open(filename, rb) # 1. 读文本卷头 text_raw f.read(3200) try: text text_raw.decode(ascii, errorsreplace) except UnicodeDecodeError: text text_raw.decode(cp037, errorsreplace) # EBCDIC # 2. 读二进制卷头 bin_raw f.read(400) sample_interval struct.unpack(endian H, bin_raw[5:7])[0] samples_per_trace struct.unpack(endian H, bin_raw[9:11])[0] format_code struct.unpack(endian H, bin_raw[15:17])[0] # 3. 循环读道 traces [] headers [] while True: head_raw f.read(240) if len(head_raw) 240: break # 道头关键字段 trace_seq struct.unpack(endian I, head_raw[0:4])[0] line_num struct.unpack(endian i, head_raw[72:76])[0] cdp_num struct.unpack(endian i, head_raw[76:80])[0] sx struct.unpack(endian i, head_raw[116:120])[0] sy struct.unpack(endian i, head_raw[120:124])[0] rx struct.unpack(endian i, head_raw[180:184])[0] ry struct.unpack(endian i, head_raw[184:188])[0] coord_unit struct.unpack(endian H, head_raw[114:116])[0] headers.append({ trace_seq: trace_seq, line_num: line_num, cdp_num: cdp_num, sx: sx, sy: sy, rx: rx, ry: ry, coord_unit: coord_unit }) # 读一道数据 trace_raw f.read(4 * samples_per_trace) if len(trace_raw) 4 * samples_per_trace: break # 根据格式码解析 if format_code 1: vals struct.unpack(endian I * samples_per_trace, trace_raw) trace [ibm_to_ieee(v) for v in vals] elif format_code 5: trace list(struct.unpack(endian f * samples_per_trace, trace_raw)) elif format_code 2: vals struct.unpack(endian i * samples_per_trace, trace_raw) trace [v / 32768.0 for v in vals] # 归一化到 -1~1 else: raise ValueError(fUnsupported format code: {format_code}) traces.append(trace) f.close() return text, headers, traces逻辑说明第一步读3200字节文本头先按ASCII解码失败就按EBCDIC代码页037兜底。第二步用上一章讲的关键偏移解析采样间隔、采样点数和格式码。第三步进入主循环每道先读240字节道头按偏移列表抽坐标字段再根据format_code决定数据体解译方式。返回的三元组text, headers, traces中text是工区文字信息headers是道头字典列表traces是二维列表每一行是一道采样值序列。3.2 参数怎么设endian、坐标单位与格式码的取舍readsegy在真实使用中有四个参数几乎必调不调或者调错结果就是数据读出来但坐标对不上、振幅对不上。第一个是字节序endian。标准SEG-Y是大端但1990年代有些国产软件写文件时直接用了本机的小端序这个只有在读取时才发现。判断方法简单粗暴先按大端读读出来的采样间隔如果是一个几十万的不正常数字或者道头坐标全是负数就换成小端再试。readsegy做成参数就是给这个留后路。第二个是坐标单位。道头里115-116字节有坐标单位标记1表示米2表示英尺。但很多老数据这个字段根本没填填0。遇到这种情况不能默默认作米要结合工区常识判断——一条测线跨度几千个单位的大概率是英尺。readsegy不帮你换算只是把这个原始值原样透传换算逻辑留给使用者。第三个是format_code。前面代码里默认支持格式码1、2、5但还有一些数据用格式码316位定点数和格式码864位浮点。如果公司内部的老软件处理出了格式码3你需要自己扩展解析分支——格式码3按h解出来再除以32。实际处理时遇到不支持的格式码先别急着改代码用十六进制编辑器看一眼文件头确认格式码再决定怎么扩。第四个是采样点数不一致的问题。道头里有“本道采样点数”字段字节17-20和二进制卷头里的“每道采样点数”可能不一致。最后几道有时会短一截。readsegy的实现策略是优先信道头里的采样点数如果某道读出来的字节数不够按实际长度截断。这个细节在后面的避坑章细说。3.3 从readsegy到地震道数组数据如何组织实际处理时极少有人直接拿readsegy返回的list套list去干活因为list的索引访问慢切片也麻烦。常见做法是把它转成numpy数组import numpy as np text, headers, traces readsegy(survey.sgy) data np.array(traces, dtypenp.float32) print(data.shape) # (道数, 每道采样点数)转成np.float32数组后后续的FFT、滤波、增益、画剖面都可以直接用向量化操作。这时候readsegy返回的headers列表要和data的每一行一一对应——第i道的坐标就是headers[i]。这个对应关系看似天然但实际操作中如果中间做了道编辑比如删除了坏道headers和数据行就会错位所以我在工程里会把headers直接转成结构化数组和数据一起存储import numpy as np header_dtype np.dtype([ (trace_seq, i4), (line_num, i4), (cdp_num, i4), (sx, i4), (sy, i4), (rx, i4), (ry, i4), (coord_unit, i2) ]) header_arr np.array([tuple(h.values()) for h in headers], dtypeheader_dtype)这样headers_arr和data在行维度上严格对应之后做CDP分选或者inline/xline重排时挪数据行必须同步挪headers行一个切片操作就能同时搞定。4. 常见做法对比segyio、obspy、segpy与readsegy的取舍readsegy不是市面上唯一的解法。地震数据处理领域还有segyio、obspy、segpy这些常用工具各有各的适用场景。这一章做一个横向对比帮你在实际项目中选型。4.1 四个库的能力边界对比库底层实现读取速度支持格式码适用场景readsegy自实现纯Python struct中慢取决于自定义教学、快速探索、格式码特殊的私有数据segyioC扩展快1、2、5、8大规模生产、需要随机访问道数据obspyPython numpy中1、2、5地震学波形处理、与地震目录联动segpy纯Python构造和解析慢1、2、5写SEG-Y文件、轻量解析速度差距在文件尺寸上体现非常明显。一个三维工区可能有几千条测线、每条测线上百万个地震道单个文件几十GB。纯Python写循环逐道读取一小时能读完就算不错segyio用C扩展加内存映射十分钟内存就能把道头索引建完。所以生产环境我几乎不用自实现的readsegy去读大文件它更多是用来快速确认文件格式、看道头信息、或者调试格式码解析逻辑。4.2 大文件场景readsegy的块读取与内存控制文件一大内存就成了瓶颈。一个地震道采样点数假设是2000按4字节浮点算一道占8KB一百万个地震道就是8GB。一次性全读进内存32GB的机器也能被撑到边缘。实际处理时我习惯按“块”来组织读取逻辑def readsegy_block(filename, start_trace, num_traces, endian): f open(filename, rb) f.seek(3600) f.seek(start_trace * (240 4 * n_samples), 1) block_head [] block_data [] for _ in range(num_traces): head_raw f.read(240) if len(head_raw) 240: break block_head.append(parse_trace_header(head_raw)) trace_raw f.read(4 * n_samples) block_data.append(decode_trace(trace_raw)) f.close() return block_head, block_data这里的关键是计算偏移位置每一道在文件里占据的字节数是 240道头 采样点数 × 每采样点字节数。要在文件里找第 start_trace 道直接用乘法跳到该位置就行不用从头顺序读。这在处理大型三维数据做inline切片时特别有用——先读道头索引定位到某条inline的所有道再只读那一段数据体内存占用可以压缩到原来的十分之一。4.3 什么时候别用readsegy坐标校正与道头重写场景readsegy适合读但遇到需要重写SGY文件的场景就力不从心了。比如野外采集的数据坐标写错了要做批量坐标偏移校正这时候不仅要读还要改道头、写文件、保持字节对齐。segyio在这方面有原生优势它提供按道号随机读写的接口改一个道的坐标字段不需要重写整个文件。我一般这样分工格式探索和快速绘图用readsegy正式的生产管线道头校正、重排、输出重新打包用segyio。另外obsPy是个地震学通用框架适合处理天然地震波形数据它读SEG-Y的能力是顺带支持的如果数据里要同时处理事件目录、波形拾取这些地震学概念才值得引入。5. 读取SEG-Y的五个真实踩坑字节序、IBM浮点、延长道与坐标单位这一章写的是我在实际项目中踩过的五个坑每一条都对应一个真实遇到过的数据集。写出来一是让你少花排查时间二是让你知道遇到这些异常时该去哪里看。5.1 IBM浮点转换导致振幅全乱现象读出的地震道数值范围在±1e20附近甚至出现NaN画出来的剖面全是噪点。 原因格式码为1的数据用的是IBM浮点但转换函数写错了。最常见的错误是把IBM的指数当成基2而不是基16或者尾数没按有符号数处理。还有个隐蔽情况某些老软件写IBM浮点时高低字节互换I解出来的整数整体对不上。 解决先取一道数据用struct.unpack(I, raw)得到整数再按上一章的ibm_to_ieee手动转换验证最大最小值是否在物理合理范围速度、振幅通常不超过几千到几万量级。如果转换后数值还是异常用十六进制编辑器直接打开那一道的原始字节肉眼比对IBM浮点的标准编码。5.2 大端小端之争采样间隔异常大现象读出来的采样间隔不是常见的2000微秒2毫秒或4000微秒4毫秒而是几十万微秒。 原因这是典型的字节序问题。文件实际是小端序但readsegy默认按大端解。struct.unpack(H, bin_raw[5:7])把字节序反过来之后2000微秒的旧数据变成了一个天文数字。 解决把readsegy的endian参数改成重新读一遍。如果恢复正常说明这个文件确实是用小端序写的。这时要记录下这个参数每一道解析都要用相同字节序。5.3 延长道Extended Trace导致读取中断现象文件末尾多出一大段空白道readsegy读到一半报错或者读出来的数据最后几行全是0。 原因SEG-Y标准规定每道数据体长度等于 二进制卷头里的采样点数 × 每采样点字节数。但有些处理软件为了对齐到磁盘块在文件末尾追加了“延长道”Extended Trace长度不定里面填0。readsegy不知道文件里多了这玩意按正常道解析结果最后一条道长度对不上。 解决主循环里加一条条件判断如果剩余文件长度小于“道头240 道数据长度”就停止读取。另外查看二进制卷头的字节3217-3218里有个“扩展道数”字段有些软件会写这里可以用它来判断末尾是否有延长道。我的实现里选择直接忽略延长道——它们在处理中本来就没有物理意义。5.4 坐标单位不是1也不是2CDP坐标对不上位置现象读出来的道头坐标能和工区底图对上但缩放比例完全不对比如几百米的道间距变成了几千。 原因二进制卷头里坐标单位字段写的是0或空白。很多老软件写道头时忘了填坐标单位readsegy读出来就是0使用者默认当米处理实际上原数据可能是英尺。坐标数值往往大了3.28倍。 解决读取后检查coord_unit字段如果为0需要结合工区已知参数反推。一个土办法取一条测线连续10道的X坐标算平均道间距如果间距在20-50之间可能是米在60-170之间大概是英尺。确定后自己做一次换算再传给下游。5.5 文本头解码乱码EBCDIC和ASCII的混合体现象readsegy返回的text字段里一半是英文正常字符一半是奇怪的符号和空格。 原因文本卷头标准定义可以用ASCII或EBCDIC。有些写盘软件只在第1字节写了个标识符其余还是EBCDIC另一些则干脆混合使用。readsegy只做了一轮尝试性解码碰到混合编码就翻车。 解决先检查字节3200-3200文本卷头第一个字节看它是C1EBCDIC的A还是41ASCII的A。如果是C1整段用cp037解码如果是41用ascii。混合编码的文件解码后人工看一遍把乱码段手工替换。坦白说到了这一步已经是玄学范畴了我的经验是别在文本头上纠缠工区信息最后都会写进数据库文本卷头过一遍眼确认工区名没错就行。6. 验证readsegy读得对不对与segyio交叉核对和可视化抽查读完SGY数据之后最怕的不是报错而是数据安安静静地读出来了但结果不可信。我现在的固定动作是做交叉验证用readsegy和segyio分别读同一个文件对比道头坐标和振幅统计量。import segyio import numpy as np # 用 segyio 读同一文件做核对 with segyio.open(survey.sgy, ignore_geometryTrue) as segy_f: segy_trace segy_f.trace[0] segy_header segy_f.header[0] # 用 readsegy 读同样的一道 text, headers, traces readsegy(survey.sgy) my_trace np.array(traces[0]) my_header headers[0] # 核对道头坐标 print(segyio sx:, segy_header[segyio.su.sx], readsegy sx:, my_header[sx]) # 核对振幅数组 print(max abs diff:, np.max(np.abs(segy_trace - my_trace)))如果最大振幅差接近零说明两道解析路径对上了如果差得很大把格式码和字节序检查一遍。另一个验证方法是可视化抽查随机抽几道画成地震波形图和商业软件如Vista、Omega显示的剖面做肉眼对比。波形形状一致、振幅量级一致才敢继续往下走。一次在甲方项目里我用readsegy读了他们的老磁带转储数据故事都进行到振幅补偿了最后交叉验证时发现有一半文件的坐标单位是英尺。那次之后我学了个规矩每个SGY文件进处理流程之前先跑这个双读核对脚本不通过就不进下一步。那个核对脚本至今还在我工具目录里躺着。希望帮到你。本文还有配套的精品资源点击获取
