简介本资源是一份面向地球物理勘探人员、地质工程专业师生及MATLAB地震数据处理初学者的实用工具脚本聚焦解决SGYSeg-Y格式地震数据在MATLAB平台上的解析难题。SGY作为国际通用的地震数据存储标准其二进制结构复杂、元数据丰富手动读取门槛高该资源提供轻量级、可直接调用的readsegy.m函数完整实现文件头解析、道头提取、样本数据读取与矩阵化转换全流程显著降低地震数据导入与预处理难度。压缩包仅含1个MATLAB源文件.m大小仅1KB精炼无冗余便于嵌入现有分析流程或二次开发。目前已有789人学习下载适用于石油勘探、地震反演、教学实验等场景读者可直接部署运行快速获取结构化地震道矩阵及关键元数据字段为后续滤波、叠加、成像等分析奠定可靠数据基础。1. 从一次数据读取失败说起SGY格式的“隐形门槛”最近在帮一个做地球物理勘探的朋友处理一批地震数据他发来一个压缩包里面全是.sgy文件。他原话是“哥们儿用Python帮我读一下画个剖面看看应该很简单吧”我心想读取一个二进制文件用struct或者numpy.fromfile不就搞定了吗结果一上手直接懵了。文件是读进来了但显示出来的剖面图全是乱码要么就是波形挤成一团要么就是振幅值完全不对根本没法看。这就是SGYSEG-Y格式给很多初入地球物理数据处理甚至是一些有经验的程序员设下的“隐形门槛”。它看起来就是一个后缀为.sgy或.segy的二进制文件网上也能搜到所谓的“标准”格式说明。但当你真正动手去读时会发现所谓的“标准”里充满了厂商自定义、历史遗留的“坑”。readsegy这个关键词背后远不是一句f.read()那么简单它涉及对地球物理数据采集、存储历史的深刻理解以及对二进制数据结构的精细操作。简单来说SGY是勘探地球物理领域的“事实标准”用于存储地震勘探采集到的地震道数据。一次勘探可能产生成千上万个这样的文件每个文件包含几十到几千条“地震道”每条道则是随时间变化的地震波振幅值。我们的核心任务就是准确无误地将这些二进制数字还原成有物理意义的矩阵供后续偏移成像、属性分析等使用。这个过程就像在破解一份没有固定密码本的加密电报你需要根据“报头”里的线索去推断数据的真实排列方式。2. SGY文件结构深度拆解不只是3200字节的文本头很多人以为SGY文件就是“一个文本头 一堆数据”。这个认知太粗略了是导致读取失败的主要原因。我们必须像外科手术一样精确地解剖它的结构。一个标准的SGY文件可以划分为四个逻辑部分每一部分都有其特定的使命和陷阱。2.1 第一部分3200字节的EBCDIC文本头这是文件开头的3200个字节。第一个坑就在这里它通常是用EBCDIC编码的文本而不是我们熟悉的ASCII或UTF-8。如果你直接用‘utf-8’去解码得到的就是一堆乱码。with open(‘survey.sgy‘, ‘rb‘) as f: ebcdic_header f.read(3200) # 错误做法直接尝试解码为utf-8 # text_header ebcdic_header.decode(‘utf-8‘, errors‘ignore‘) # 会得到乱码 # 正确做法先转换编码需安装ebcdic库或使用转换表 # 方法一使用codecs部分系统支持 import codecs try: text_header ebcdic_header.decode(‘cp500‘) # IBM EBCDIC 500 是常见编码 except: # 方法二更通用的做法是使用专用库或手动映射复杂 pass这段文本头里包含了工区名、采集公司、处理历史等描述性信息对于数据管理至关重要但对于纯粹的数据读取readsegy的核心来说它不是必须正确解析的。一个务实的做法是如果不需要这些文本信息可以跳过这3200字节或者仅将其作为二进制块保存。如果需要则必须处理EBCDIC编码问题。2.2 第二部分400字节的二进制文件头紧接在文本头之后是400个字节的二进制文件头。这是第一个关键数据结构。它包含了描述整个文件数据体的全局参数这些参数控制着如何解释后续所有的地震道数据。这里有几个生死攸关的字段错一个全盘皆错JobID,LineNumber,ReelNumber等字节位置 0-120这些是标识信息相对安全。DataTracesPerEnsemble字节 120-122每个“道集”包含的数据道数。对于简单的2D数据这个值通常就是文件里的总道数或者一个很大的数如32767。但要注意它是short类型2字节。AuxTracesPerEnsemble字节 122-124辅助道数通常为0。SampleInterval字节 116-118采样间隔单位是微秒μs。这是核心中的核心它决定了时间轴的刻度。例如值4000代表采样间隔是4毫秒ms。常见的有1ms1000、2ms2000、4ms4000。读取时必须确认其字节顺序。SamplesPerTrace字节 114-116每道采样点数。它和SampleInterval共同确定了每道数据的时间长度。例如1501个采样点4ms采样间隔则一道的时间长度是 (1501-1)*4ms 6秒。DataSampleFormatCode字节 124-126数据样本格式码。这是第二个核心关键直接决定数据体的二进制格式。以下是常见值及其含义1 IBM 32位浮点数旧格式麻烦2 32位有符号整数不常见3 16位有符号整数最常见5 IEEE 32位浮点数现代格式越来越常见8 8位有符号整数不常见绝大多数野外采集的原始数据格式码是3即每个振幅值用2字节16位的有符号整数存储。EnsembleFold字节 126-128覆盖次数处理时有用。MeasurementSystem字节 128-130度量系统1米2英尺。影响道头中坐标的解释。SegyFormatRevisionNumber字节 300-302 SEG-Y格式版本号。0100代表Rev10200代表Rev2。Rev2有较大扩展但支持的工具较少。读取二进制文件头时必须使用struct模块或numpy.dtype来精确指定每个字段的字节位置、类型和字节顺序。字节顺序大端‘或小端‘是另一个大坑通常SEG-Y标准规定为大端序但有些软件或设备产出的是小端序。import numpy as np import struct # 假设我们已经跳过3200字节的文本头文件指针位于二进制文件头开始处 with open(‘survey.sgy‘, ‘rb‘) as f: f.seek(3200) # 定位到二进制文件头开始 bin_header f.read(400) # 使用struct解包关键字段假设大端序 ‘‘ # ‘H‘: unsigned short (2字节), ‘h‘: signed short (2字节), ‘i‘: signed int (4字节) samples_per_trace, struct.unpack(‘H‘, bin_header[114:116]) # 注意标准中这些字段可能是short需查证准确位置和类型 sample_interval, struct.unpack(‘H‘, bin_header[116:118]) data_format, struct.unpack(‘h‘, bin_header[124:126]) print(f“每道采样点数: {samples_per_trace}“) print(f“采样间隔(μs): {sample_interval}“) print(f“数据格式码: {data_format}“)实操心得1不要相信任何一个现成的“偏移量表”。不同来源的SGY文件其二进制文件头内部字段的精确字节位置可能存在变体。最可靠的方法是结合一个已知正确的SGY文件和十六进制编辑器手动验证关键字段的位置。或者使用成熟的库如segyio,obspy来读取它们内部处理了这些变体。2.3 第三部分道头与数据体循环结构这是文件的主体由N个道块依次排列而成。每个道块包括道头240字节的二进制头描述这一道的属性。数据体SamplesPerTrace个数据样本每个样本的格式由DataSampleFormatCode决定。道头里包含的信息对于后续处理至关重要道序号字节 0-4这道在文件中的顺序号。CDP号字节 20-24共深度点号是核心的横向坐标。坐标字节 72-84, 84-96X和Y坐标。注意其缩放因子字节 140-144, 144-148实际坐标 头值 * 10^(-缩放因子)。缩放因子为负时是除法为正时是乘法这里极易出错。接收点高程/深度等。读取数据体的逻辑是循环读取240字节道头。根据SamplesPerTrace和DataSampleFormatCode计算出数据体的大小。读取数据体并按照指定的格式和字节顺序解码成数值数组。# 续前文假设已知关键参数 samples samples_per_trace format_code data_format # 假设是3即16位有符号整数 trace_data_size samples * 2 # 因为格式3是2字节/样本 traces [] with open(‘survey.sgy‘, ‘rb‘) as f: f.seek(3600) # 跳过3200文本头 400二进制头定位到第一个道块 while True: trace_header f.read(240) if not trace_header: break # 解析道头信息例如CDP号 cdp, struct.unpack(‘i‘, trace_header[20:24]) # 读取数据体 trace_body_bytes f.read(trace_data_size) if len(trace_body_bytes) trace_data_size: break # 文件可能意外结束 # 根据格式码解码数据体 if format_code 3: # 16位有符号整数 # 使用numpy从缓冲区直接转换指定大端序 ‘‘ 和类型 ‘i2‘ (有符号16位整数) trace_data np.frombuffer(trace_body_bytes, dtype‘i2‘) elif format_code 5: # IEEE 32位浮点数 trace_data np.frombuffer(trace_body_bytes, dtype‘f4‘) else: raise ValueError(f“不支持的数据格式码: {format_code}“) # 存储道数据和头信息 traces.append({‘cdp‘: cdp, ‘data‘: trace_data}) print(f“成功读取 {len(traces)} 道数据。“)2.4 第四部分道头映射与变体这才是最让人头疼的地方。SEG-Y标准特别是早期的Rev1只定义了道头中部分字段的含义大量字节位置是留给“用户自定义”的。不同的采集设备制造商、不同的处理软件会把这些信息如炮点号、接收点号、偏移距、方位角等放在不同的道头字节位置。例如偏移距这个关键参数标准没有规定其位置。软件A可能把它放在字节37-40软件B可能放在字节69-72。如果你用软件A的映射表去读软件B生成的文件偏移距就会读错导致后续动校正等处理完全失败。实操心得2在开始大规模读取数据前务必先进行数据验证。选择文件中间的一道用十六进制编辑器查看其道头并与已知的、正确的道头信息例如从处理软件中导出的道头列表进行比对确认关键字段CDP、偏移距、坐标的字节位置是否正确。或者使用一个能正确显示该数据的商业/开源软件读入导出其道头再与你自己的解析结果对比。3. 实战“readsegy”手动实现与成熟库的抉择理解了结构我们就可以动手实现一个readsegy函数了。但在此之前必须做一个重要的抉择是手动造轮子还是使用成熟库3.1 方案一手动实现核心读取逻辑对于学习原理或处理格式非常固定的数据手动实现是有价值的。下面是一个高度简化的示例聚焦于核心数据体的读取忽略了复杂的道头解析和变体处理。import numpy as np import struct def read_segy_basic(filepath, format_code3, endian‘‘): “““ 一个基础的SGY读取函数仅读取数据体为numpy数组。 假设文本头3200B二进制头400B道头240B数据格式和采样点数已知或通过参数传入。 此函数省略了从文件头自动检测参数的过程实际应用需补全。 “““ # 在实际应用中这些参数应从二进制文件头读取 # 此处作为参数传入或硬编码仅用于演示 samples_per_trace 1501 trace_header_size 240 # 根据格式码确定dtype if format_code 3: data_dtype np.dtype(endian ‘i2‘) # 有符号16位整型 bytes_per_sample 2 elif format_code 5: data_dtype np.dtype(endian ‘f4‘) # 32位浮点型 bytes_per_sample 4 else: raise ValueError(f“暂不支持格式码: {format_code}“) trace_data_size samples_per_trace * bytes_per_sample trace_block_size trace_header_size trace_data_size data_list [] with open(filepath, ‘rb‘) as f: # 跳过文件头 f.seek(3600) while True: # 跳过道头 f.seek(trace_header_size, 1) # 读取数据体 trace_bytes f.read(trace_data_size) if not trace_bytes or len(trace_bytes) trace_data_size: break # 转换为numpy数组并存储 trace_data np.frombuffer(trace_bytes, dtypedata_dtype).copy() # .copy()使其内存连续 data_list.append(trace_data) # 将列表转换为2D numpy数组 [道数, 采样点数] if data_list: data_matrix np.vstack(data_list) return data_matrix else: return np.array([]) # 使用示例 try: seismic_data read_segy_basic(‘your_data.sgy‘, format_code3, endian‘‘) print(f“数据形状: {seismic_data.shape}“) # 输出 (道数, 1501) except Exception as e: print(f“读取失败: {e}“)这个简易函数的局限性无法自动从文件头读取SamplesPerTrace和DataSampleFormatCode。完全忽略了道头信息丢失了所有空间属性。没有处理字节顺序自动检测。没有处理文件可能存在的额外尾随字节。3.2 方案二使用成熟开源库推荐对于生产环境或科研强烈推荐使用成熟的开源库它们已经踩过了所有的坑。segyioPython和C库性能极佳API简洁是业界许多商业软件的幕后英雄。它提供了类似numpy的切片操作是处理大型SGY文件的首选。import segyio with segyio.open(‘survey.sgy‘, ‘r‘, ignore_geometryTrue) as f: # 快速读取所有数据到一个numpy数组注意内存 data segyio.tools.cube(f) # 或者逐道读取 for trace in f.trace: process(trace) # 获取道头信息 cdp_numbers f.attributes(segyio.TraceField.CDP)[:]obspy地震学领域的瑞士军刀。它的read函数可以直接读取SGY并返回一个Stream对象内部是Trace对象每个Trace都带有丰富的头信息非常方便。from obspy import read stream read(‘survey.sgy‘) print(stream) # 显示包含多少道 for trace in stream: print(trace.stats) # 查看道头信息 print(trace.data) # 数据数组lasio虽然主要针对测井数据LAS但有时也被用来读取简单的SGY不推荐用于复杂情况。实操心得3segyio在性能和内存控制上更优适合处理海量数据。obspy在数据封装和地震学特定处理上更友好。如果你的数据能被obspy正确识别道头那么用obspy会省心很多。如果obspy读取出错通常是道头映射问题可以尝试用segyio的strictFalse模式或者直接使用segyio。4. 读取后的关键步骤数据验证与质量检查成功将二进制数据读入内存成为numpy数组只是万里长征第一步。接下来必须进行严格的数据验证否则垃圾数据进去垃圾结果出来。4.1 可视化快速检查绘制几个典型的道或者绘制整个剖面的灰度图/波形图。import matplotlib.pyplot as plt data read_segy_basic(‘data.sgy‘, ...) # 或从segyio/obspy获取 # 检查单道波形 plt.figure(figsize(12, 4)) plt.subplot(1, 2, 1) plt.plot(data[100, :]) # 第101道 plt.title(‘Trace 100‘) plt.xlabel(‘Sample Index‘) plt.ylabel(‘Amplitude‘) # 检查剖面通常需要增益调整 plt.subplot(1, 2, 2) # 对数据进行裁剪或归一化以便显示 clip_percentile 98 vmax np.percentile(np.abs(data), clip_percentile) plt.imshow(data.T, aspect‘auto‘, cmap‘seismic‘, vmin-vmax, vmaxvmax) # 转置使时间为纵轴 plt.colorbar(label‘Amplitude‘) plt.title(‘Seismic Section‘) plt.xlabel(‘Trace Number‘) plt.ylabel(‘Time Sample‘) plt.tight_layout() plt.show()观察点单道波形是否合理有无明显的直流偏移波形整体不在零线附近剖面中的同相轴是否连续有无明显的条带状噪声或异常道4.2 统计信息检查计算数据的基本统计量与预期对比。print(f“数据形状: {data.shape}“) print(f“数据范围: [{data.min():.2f}, {data.max():.2f}]“) print(f“数据均值: {data.mean():.2f} (应接近0)“) print(f“数据标准差: {data.std():.2f}“) print(f“NaN值数量: {np.isnan(data).sum()}“) print(f“Inf值数量: {np.isinf(data).sum()}“)预期原始地震数据均值应非常接近0。如果均值很大可能存在直流分量。NaN或Inf值表明读取过程或原始数据有问题。4.3 道头信息一致性检查如果你读取了道头信息如CDP号检查它们是否单调递增有无跳号或重复。cdp_numbers ... # 从道头读取的CDP号数组 plt.plot(cdp_numbers) plt.xlabel(‘Trace Index‘) plt.ylabel(‘CDP Number‘) plt.title(‘CDP Number Sequence‘) plt.grid(True) plt.show() # 检查跳变 diff np.diff(cdp_numbers) print(f“CDP号最大跳变: {diff.max()}“) print(f“CDP号减少的次数: {(diff 0).sum()}“)不连续的CDP号可能意味着数据缺失或排序问题需要在后续处理中注意。5. 高级话题与避坑指南5.1 字节顺序Endianness的自动检测最稳妥的方法是尝试两种字节顺序来读取二进制文件头中的某个已知字段。例如SamplesPerTrace假设在3200114字节处应该是一个合理的正数比如几百到几千。我们可以用大端和小端分别去解包这个字段哪个结果看起来合理就采用哪种字节顺序。def detect_endianness(filepath): with open(filepath, ‘rb‘) as f: f.seek(3314) # 3200 114 SamplesPerTrace的位置 bytes_ f.read(2) # 尝试大端序解包 samples_be, struct.unpack(‘H‘, bytes_) # 尝试小端序解包 samples_le, struct.unpack(‘H‘, bytes_) # 判断哪个值更“合理” if 100 samples_be 10000: # 合理范围 return ‘‘ # 大端序 elif 100 samples_le 10000: return ‘‘ # 小端序 else: # 如果都不合理可以尝试其他已知字段如格式码应在1-8之间 raise ValueError(“无法自动检测字节顺序请手动指定。“)5.2 处理“非标准”SGY文件你可能会遇到没有3200字节文本头直接就是二进制文件头。这时需要调整偏移量。文件末尾有多余字节有些软件会在文件末尾添加一些注释或填充。在循环读取道块时如果最后读到的数据块不足一个道块大小应优雅地跳出循环。数据体被压缩或加密极少见通常需要专门的软件或密钥。5.3 内存管理与大数据处理一个三维工区的SGY文件可能达到GB甚至TB级别。一次性读入内存segyio.tools.cube或numpy.memmap的完全加载会导致内存溢出。解决方案使用segyio的迭代器for trace in segyfile.trace:可以逐道处理。内存映射文件对于格式规整的文件可以用numpy.memmap但需要自己精确计算偏移量非常复杂。分块读取使用segyio或自定义代码每次只读取一个范围的道如traces[1000:2000]。使用Dask对于超大规模数据可以考虑使用dask.array来构建延迟计算的数组但需要编写适配器。实操心得4在处理超大文件前先用segyio或快速扫描获取总道数和每道采样数估算内存占用总道数 * 每道采样数 * 4字节对于浮点数。如果远超物理内存就必须设计流式或分块处理方案。5.4 从读取到应用数据坐标系的建立正确读取振幅数据只是拿到了“像素值”。要得到有地理意义的地震剖面必须结合道头信息建立坐标系。横向坐标通常使用CDP号作为索引。更精确的则需要使用道头中的CDP_X和CDP_Y坐标需应用缩放因子。纵向坐标时间轴根据二进制文件头中的SampleInterval微秒和SamplesPerTrace计算时间向量time np.arange(samples_per_trace) * (sample_interval / 1e6)单位转换为秒。绘制剖面使用matplotlib的imshow时通过extent参数传入横向和纵向的范围即可将图像坐标映射到物理坐标。# 假设已获取 sample_interval_us 4000 # 微秒 samples 1501 cdp_numbers np.array([...]) # 从道头读取的CDP号 # 创建时间轴单位秒 time_axis np.arange(samples) * (sample_interval_us / 1e6) # 绘制带有物理坐标的剖面 plt.imshow(data.T, aspect‘auto‘, cmap‘seismic‘, extent[cdp_numbers.min(), cdp_numbers.max(), time_axis.max(), time_axis.min()]) # 注意时间轴上下翻转 plt.xlabel(‘CDP Number‘) plt.ylabel(‘Time (s)‘) plt.colorbar() plt.show()读取SGY文件是一个连接数字世界与物理世界的过程。每一个字节背后都对应着地下某一点在某个时刻的振动强度。代码的精确性直接决定了我们看到的“地下影像”是否真实可靠。从混乱的二进制流中精准地提取出有意义的数字矩阵和空间标签是进行任何高级地震解释和处理不可动摇的基石。我个人的体会是与其在后期处理中纠结算法效果不好不如花双倍的时间在前期的数据读取和验证上确保“喂”给算法的数据是干净、正确的。很多时候问题不是出在复杂的数学公式上恰恰就出在最基础的struct.unpack那个格式字符串里。本文还有配套的精品资源点击获取 SEG-Y地震数据读取全解析：从二进制结构到Python实战
简介本资源是一份面向地球物理勘探人员、地质工程专业师生及MATLAB地震数据处理初学者的实用工具脚本聚焦解决SGYSeg-Y格式地震数据在MATLAB平台上的解析难题。SGY作为国际通用的地震数据存储标准其二进制结构复杂、元数据丰富手动读取门槛高该资源提供轻量级、可直接调用的readsegy.m函数完整实现文件头解析、道头提取、样本数据读取与矩阵化转换全流程显著降低地震数据导入与预处理难度。压缩包仅含1个MATLAB源文件.m大小仅1KB精炼无冗余便于嵌入现有分析流程或二次开发。目前已有789人学习下载适用于石油勘探、地震反演、教学实验等场景读者可直接部署运行快速获取结构化地震道矩阵及关键元数据字段为后续滤波、叠加、成像等分析奠定可靠数据基础。1. 从一次数据读取失败说起SGY格式的“隐形门槛”最近在帮一个做地球物理勘探的朋友处理一批地震数据他发来一个压缩包里面全是.sgy文件。他原话是“哥们儿用Python帮我读一下画个剖面看看应该很简单吧”我心想读取一个二进制文件用struct或者numpy.fromfile不就搞定了吗结果一上手直接懵了。文件是读进来了但显示出来的剖面图全是乱码要么就是波形挤成一团要么就是振幅值完全不对根本没法看。这就是SGYSEG-Y格式给很多初入地球物理数据处理甚至是一些有经验的程序员设下的“隐形门槛”。它看起来就是一个后缀为.sgy或.segy的二进制文件网上也能搜到所谓的“标准”格式说明。但当你真正动手去读时会发现所谓的“标准”里充满了厂商自定义、历史遗留的“坑”。readsegy这个关键词背后远不是一句f.read()那么简单它涉及对地球物理数据采集、存储历史的深刻理解以及对二进制数据结构的精细操作。简单来说SGY是勘探地球物理领域的“事实标准”用于存储地震勘探采集到的地震道数据。一次勘探可能产生成千上万个这样的文件每个文件包含几十到几千条“地震道”每条道则是随时间变化的地震波振幅值。我们的核心任务就是准确无误地将这些二进制数字还原成有物理意义的矩阵供后续偏移成像、属性分析等使用。这个过程就像在破解一份没有固定密码本的加密电报你需要根据“报头”里的线索去推断数据的真实排列方式。2. SGY文件结构深度拆解不只是3200字节的文本头很多人以为SGY文件就是“一个文本头 一堆数据”。这个认知太粗略了是导致读取失败的主要原因。我们必须像外科手术一样精确地解剖它的结构。一个标准的SGY文件可以划分为四个逻辑部分每一部分都有其特定的使命和陷阱。2.1 第一部分3200字节的EBCDIC文本头这是文件开头的3200个字节。第一个坑就在这里它通常是用EBCDIC编码的文本而不是我们熟悉的ASCII或UTF-8。如果你直接用‘utf-8’去解码得到的就是一堆乱码。with open(‘survey.sgy‘, ‘rb‘) as f: ebcdic_header f.read(3200) # 错误做法直接尝试解码为utf-8 # text_header ebcdic_header.decode(‘utf-8‘, errors‘ignore‘) # 会得到乱码 # 正确做法先转换编码需安装ebcdic库或使用转换表 # 方法一使用codecs部分系统支持 import codecs try: text_header ebcdic_header.decode(‘cp500‘) # IBM EBCDIC 500 是常见编码 except: # 方法二更通用的做法是使用专用库或手动映射复杂 pass这段文本头里包含了工区名、采集公司、处理历史等描述性信息对于数据管理至关重要但对于纯粹的数据读取readsegy的核心来说它不是必须正确解析的。一个务实的做法是如果不需要这些文本信息可以跳过这3200字节或者仅将其作为二进制块保存。如果需要则必须处理EBCDIC编码问题。2.2 第二部分400字节的二进制文件头紧接在文本头之后是400个字节的二进制文件头。这是第一个关键数据结构。它包含了描述整个文件数据体的全局参数这些参数控制着如何解释后续所有的地震道数据。这里有几个生死攸关的字段错一个全盘皆错JobID,LineNumber,ReelNumber等字节位置 0-120这些是标识信息相对安全。DataTracesPerEnsemble字节 120-122每个“道集”包含的数据道数。对于简单的2D数据这个值通常就是文件里的总道数或者一个很大的数如32767。但要注意它是short类型2字节。AuxTracesPerEnsemble字节 122-124辅助道数通常为0。SampleInterval字节 116-118采样间隔单位是微秒μs。这是核心中的核心它决定了时间轴的刻度。例如值4000代表采样间隔是4毫秒ms。常见的有1ms1000、2ms2000、4ms4000。读取时必须确认其字节顺序。SamplesPerTrace字节 114-116每道采样点数。它和SampleInterval共同确定了每道数据的时间长度。例如1501个采样点4ms采样间隔则一道的时间长度是 (1501-1)*4ms 6秒。DataSampleFormatCode字节 124-126数据样本格式码。这是第二个核心关键直接决定数据体的二进制格式。以下是常见值及其含义1 IBM 32位浮点数旧格式麻烦2 32位有符号整数不常见3 16位有符号整数最常见5 IEEE 32位浮点数现代格式越来越常见8 8位有符号整数不常见绝大多数野外采集的原始数据格式码是3即每个振幅值用2字节16位的有符号整数存储。EnsembleFold字节 126-128覆盖次数处理时有用。MeasurementSystem字节 128-130度量系统1米2英尺。影响道头中坐标的解释。SegyFormatRevisionNumber字节 300-302 SEG-Y格式版本号。0100代表Rev10200代表Rev2。Rev2有较大扩展但支持的工具较少。读取二进制文件头时必须使用struct模块或numpy.dtype来精确指定每个字段的字节位置、类型和字节顺序。字节顺序大端‘或小端‘是另一个大坑通常SEG-Y标准规定为大端序但有些软件或设备产出的是小端序。import numpy as np import struct # 假设我们已经跳过3200字节的文本头文件指针位于二进制文件头开始处 with open(‘survey.sgy‘, ‘rb‘) as f: f.seek(3200) # 定位到二进制文件头开始 bin_header f.read(400) # 使用struct解包关键字段假设大端序 ‘‘ # ‘H‘: unsigned short (2字节), ‘h‘: signed short (2字节), ‘i‘: signed int (4字节) samples_per_trace, struct.unpack(‘H‘, bin_header[114:116]) # 注意标准中这些字段可能是short需查证准确位置和类型 sample_interval, struct.unpack(‘H‘, bin_header[116:118]) data_format, struct.unpack(‘h‘, bin_header[124:126]) print(f“每道采样点数: {samples_per_trace}“) print(f“采样间隔(μs): {sample_interval}“) print(f“数据格式码: {data_format}“)实操心得1不要相信任何一个现成的“偏移量表”。不同来源的SGY文件其二进制文件头内部字段的精确字节位置可能存在变体。最可靠的方法是结合一个已知正确的SGY文件和十六进制编辑器手动验证关键字段的位置。或者使用成熟的库如segyio,obspy来读取它们内部处理了这些变体。2.3 第三部分道头与数据体循环结构这是文件的主体由N个道块依次排列而成。每个道块包括道头240字节的二进制头描述这一道的属性。数据体SamplesPerTrace个数据样本每个样本的格式由DataSampleFormatCode决定。道头里包含的信息对于后续处理至关重要道序号字节 0-4这道在文件中的顺序号。CDP号字节 20-24共深度点号是核心的横向坐标。坐标字节 72-84, 84-96X和Y坐标。注意其缩放因子字节 140-144, 144-148实际坐标 头值 * 10^(-缩放因子)。缩放因子为负时是除法为正时是乘法这里极易出错。接收点高程/深度等。读取数据体的逻辑是循环读取240字节道头。根据SamplesPerTrace和DataSampleFormatCode计算出数据体的大小。读取数据体并按照指定的格式和字节顺序解码成数值数组。# 续前文假设已知关键参数 samples samples_per_trace format_code data_format # 假设是3即16位有符号整数 trace_data_size samples * 2 # 因为格式3是2字节/样本 traces [] with open(‘survey.sgy‘, ‘rb‘) as f: f.seek(3600) # 跳过3200文本头 400二进制头定位到第一个道块 while True: trace_header f.read(240) if not trace_header: break # 解析道头信息例如CDP号 cdp, struct.unpack(‘i‘, trace_header[20:24]) # 读取数据体 trace_body_bytes f.read(trace_data_size) if len(trace_body_bytes) trace_data_size: break # 文件可能意外结束 # 根据格式码解码数据体 if format_code 3: # 16位有符号整数 # 使用numpy从缓冲区直接转换指定大端序 ‘‘ 和类型 ‘i2‘ (有符号16位整数) trace_data np.frombuffer(trace_body_bytes, dtype‘i2‘) elif format_code 5: # IEEE 32位浮点数 trace_data np.frombuffer(trace_body_bytes, dtype‘f4‘) else: raise ValueError(f“不支持的数据格式码: {format_code}“) # 存储道数据和头信息 traces.append({‘cdp‘: cdp, ‘data‘: trace_data}) print(f“成功读取 {len(traces)} 道数据。“)2.4 第四部分道头映射与变体这才是最让人头疼的地方。SEG-Y标准特别是早期的Rev1只定义了道头中部分字段的含义大量字节位置是留给“用户自定义”的。不同的采集设备制造商、不同的处理软件会把这些信息如炮点号、接收点号、偏移距、方位角等放在不同的道头字节位置。例如偏移距这个关键参数标准没有规定其位置。软件A可能把它放在字节37-40软件B可能放在字节69-72。如果你用软件A的映射表去读软件B生成的文件偏移距就会读错导致后续动校正等处理完全失败。实操心得2在开始大规模读取数据前务必先进行数据验证。选择文件中间的一道用十六进制编辑器查看其道头并与已知的、正确的道头信息例如从处理软件中导出的道头列表进行比对确认关键字段CDP、偏移距、坐标的字节位置是否正确。或者使用一个能正确显示该数据的商业/开源软件读入导出其道头再与你自己的解析结果对比。3. 实战“readsegy”手动实现与成熟库的抉择理解了结构我们就可以动手实现一个readsegy函数了。但在此之前必须做一个重要的抉择是手动造轮子还是使用成熟库3.1 方案一手动实现核心读取逻辑对于学习原理或处理格式非常固定的数据手动实现是有价值的。下面是一个高度简化的示例聚焦于核心数据体的读取忽略了复杂的道头解析和变体处理。import numpy as np import struct def read_segy_basic(filepath, format_code3, endian‘‘): “““ 一个基础的SGY读取函数仅读取数据体为numpy数组。 假设文本头3200B二进制头400B道头240B数据格式和采样点数已知或通过参数传入。 此函数省略了从文件头自动检测参数的过程实际应用需补全。 “““ # 在实际应用中这些参数应从二进制文件头读取 # 此处作为参数传入或硬编码仅用于演示 samples_per_trace 1501 trace_header_size 240 # 根据格式码确定dtype if format_code 3: data_dtype np.dtype(endian ‘i2‘) # 有符号16位整型 bytes_per_sample 2 elif format_code 5: data_dtype np.dtype(endian ‘f4‘) # 32位浮点型 bytes_per_sample 4 else: raise ValueError(f“暂不支持格式码: {format_code}“) trace_data_size samples_per_trace * bytes_per_sample trace_block_size trace_header_size trace_data_size data_list [] with open(filepath, ‘rb‘) as f: # 跳过文件头 f.seek(3600) while True: # 跳过道头 f.seek(trace_header_size, 1) # 读取数据体 trace_bytes f.read(trace_data_size) if not trace_bytes or len(trace_bytes) trace_data_size: break # 转换为numpy数组并存储 trace_data np.frombuffer(trace_bytes, dtypedata_dtype).copy() # .copy()使其内存连续 data_list.append(trace_data) # 将列表转换为2D numpy数组 [道数, 采样点数] if data_list: data_matrix np.vstack(data_list) return data_matrix else: return np.array([]) # 使用示例 try: seismic_data read_segy_basic(‘your_data.sgy‘, format_code3, endian‘‘) print(f“数据形状: {seismic_data.shape}“) # 输出 (道数, 1501) except Exception as e: print(f“读取失败: {e}“)这个简易函数的局限性无法自动从文件头读取SamplesPerTrace和DataSampleFormatCode。完全忽略了道头信息丢失了所有空间属性。没有处理字节顺序自动检测。没有处理文件可能存在的额外尾随字节。3.2 方案二使用成熟开源库推荐对于生产环境或科研强烈推荐使用成熟的开源库它们已经踩过了所有的坑。segyioPython和C库性能极佳API简洁是业界许多商业软件的幕后英雄。它提供了类似numpy的切片操作是处理大型SGY文件的首选。import segyio with segyio.open(‘survey.sgy‘, ‘r‘, ignore_geometryTrue) as f: # 快速读取所有数据到一个numpy数组注意内存 data segyio.tools.cube(f) # 或者逐道读取 for trace in f.trace: process(trace) # 获取道头信息 cdp_numbers f.attributes(segyio.TraceField.CDP)[:]obspy地震学领域的瑞士军刀。它的read函数可以直接读取SGY并返回一个Stream对象内部是Trace对象每个Trace都带有丰富的头信息非常方便。from obspy import read stream read(‘survey.sgy‘) print(stream) # 显示包含多少道 for trace in stream: print(trace.stats) # 查看道头信息 print(trace.data) # 数据数组lasio虽然主要针对测井数据LAS但有时也被用来读取简单的SGY不推荐用于复杂情况。实操心得3segyio在性能和内存控制上更优适合处理海量数据。obspy在数据封装和地震学特定处理上更友好。如果你的数据能被obspy正确识别道头那么用obspy会省心很多。如果obspy读取出错通常是道头映射问题可以尝试用segyio的strictFalse模式或者直接使用segyio。4. 读取后的关键步骤数据验证与质量检查成功将二进制数据读入内存成为numpy数组只是万里长征第一步。接下来必须进行严格的数据验证否则垃圾数据进去垃圾结果出来。4.1 可视化快速检查绘制几个典型的道或者绘制整个剖面的灰度图/波形图。import matplotlib.pyplot as plt data read_segy_basic(‘data.sgy‘, ...) # 或从segyio/obspy获取 # 检查单道波形 plt.figure(figsize(12, 4)) plt.subplot(1, 2, 1) plt.plot(data[100, :]) # 第101道 plt.title(‘Trace 100‘) plt.xlabel(‘Sample Index‘) plt.ylabel(‘Amplitude‘) # 检查剖面通常需要增益调整 plt.subplot(1, 2, 2) # 对数据进行裁剪或归一化以便显示 clip_percentile 98 vmax np.percentile(np.abs(data), clip_percentile) plt.imshow(data.T, aspect‘auto‘, cmap‘seismic‘, vmin-vmax, vmaxvmax) # 转置使时间为纵轴 plt.colorbar(label‘Amplitude‘) plt.title(‘Seismic Section‘) plt.xlabel(‘Trace Number‘) plt.ylabel(‘Time Sample‘) plt.tight_layout() plt.show()观察点单道波形是否合理有无明显的直流偏移波形整体不在零线附近剖面中的同相轴是否连续有无明显的条带状噪声或异常道4.2 统计信息检查计算数据的基本统计量与预期对比。print(f“数据形状: {data.shape}“) print(f“数据范围: [{data.min():.2f}, {data.max():.2f}]“) print(f“数据均值: {data.mean():.2f} (应接近0)“) print(f“数据标准差: {data.std():.2f}“) print(f“NaN值数量: {np.isnan(data).sum()}“) print(f“Inf值数量: {np.isinf(data).sum()}“)预期原始地震数据均值应非常接近0。如果均值很大可能存在直流分量。NaN或Inf值表明读取过程或原始数据有问题。4.3 道头信息一致性检查如果你读取了道头信息如CDP号检查它们是否单调递增有无跳号或重复。cdp_numbers ... # 从道头读取的CDP号数组 plt.plot(cdp_numbers) plt.xlabel(‘Trace Index‘) plt.ylabel(‘CDP Number‘) plt.title(‘CDP Number Sequence‘) plt.grid(True) plt.show() # 检查跳变 diff np.diff(cdp_numbers) print(f“CDP号最大跳变: {diff.max()}“) print(f“CDP号减少的次数: {(diff 0).sum()}“)不连续的CDP号可能意味着数据缺失或排序问题需要在后续处理中注意。5. 高级话题与避坑指南5.1 字节顺序Endianness的自动检测最稳妥的方法是尝试两种字节顺序来读取二进制文件头中的某个已知字段。例如SamplesPerTrace假设在3200114字节处应该是一个合理的正数比如几百到几千。我们可以用大端和小端分别去解包这个字段哪个结果看起来合理就采用哪种字节顺序。def detect_endianness(filepath): with open(filepath, ‘rb‘) as f: f.seek(3314) # 3200 114 SamplesPerTrace的位置 bytes_ f.read(2) # 尝试大端序解包 samples_be, struct.unpack(‘H‘, bytes_) # 尝试小端序解包 samples_le, struct.unpack(‘H‘, bytes_) # 判断哪个值更“合理” if 100 samples_be 10000: # 合理范围 return ‘‘ # 大端序 elif 100 samples_le 10000: return ‘‘ # 小端序 else: # 如果都不合理可以尝试其他已知字段如格式码应在1-8之间 raise ValueError(“无法自动检测字节顺序请手动指定。“)5.2 处理“非标准”SGY文件你可能会遇到没有3200字节文本头直接就是二进制文件头。这时需要调整偏移量。文件末尾有多余字节有些软件会在文件末尾添加一些注释或填充。在循环读取道块时如果最后读到的数据块不足一个道块大小应优雅地跳出循环。数据体被压缩或加密极少见通常需要专门的软件或密钥。5.3 内存管理与大数据处理一个三维工区的SGY文件可能达到GB甚至TB级别。一次性读入内存segyio.tools.cube或numpy.memmap的完全加载会导致内存溢出。解决方案使用segyio的迭代器for trace in segyfile.trace:可以逐道处理。内存映射文件对于格式规整的文件可以用numpy.memmap但需要自己精确计算偏移量非常复杂。分块读取使用segyio或自定义代码每次只读取一个范围的道如traces[1000:2000]。使用Dask对于超大规模数据可以考虑使用dask.array来构建延迟计算的数组但需要编写适配器。实操心得4在处理超大文件前先用segyio或快速扫描获取总道数和每道采样数估算内存占用总道数 * 每道采样数 * 4字节对于浮点数。如果远超物理内存就必须设计流式或分块处理方案。5.4 从读取到应用数据坐标系的建立正确读取振幅数据只是拿到了“像素值”。要得到有地理意义的地震剖面必须结合道头信息建立坐标系。横向坐标通常使用CDP号作为索引。更精确的则需要使用道头中的CDP_X和CDP_Y坐标需应用缩放因子。纵向坐标时间轴根据二进制文件头中的SampleInterval微秒和SamplesPerTrace计算时间向量time np.arange(samples_per_trace) * (sample_interval / 1e6)单位转换为秒。绘制剖面使用matplotlib的imshow时通过extent参数传入横向和纵向的范围即可将图像坐标映射到物理坐标。# 假设已获取 sample_interval_us 4000 # 微秒 samples 1501 cdp_numbers np.array([...]) # 从道头读取的CDP号 # 创建时间轴单位秒 time_axis np.arange(samples) * (sample_interval_us / 1e6) # 绘制带有物理坐标的剖面 plt.imshow(data.T, aspect‘auto‘, cmap‘seismic‘, extent[cdp_numbers.min(), cdp_numbers.max(), time_axis.max(), time_axis.min()]) # 注意时间轴上下翻转 plt.xlabel(‘CDP Number‘) plt.ylabel(‘Time (s)‘) plt.colorbar() plt.show()读取SGY文件是一个连接数字世界与物理世界的过程。每一个字节背后都对应着地下某一点在某个时刻的振动强度。代码的精确性直接决定了我们看到的“地下影像”是否真实可靠。从混乱的二进制流中精准地提取出有意义的数字矩阵和空间标签是进行任何高级地震解释和处理不可动摇的基石。我个人的体会是与其在后期处理中纠结算法效果不好不如花双倍的时间在前期的数据读取和验证上确保“喂”给算法的数据是干净、正确的。很多时候问题不是出在复杂的数学公式上恰恰就出在最基础的struct.unpack那个格式字符串里。本文还有配套的精品资源点击获取