第一次拿到SEGY文件时我盯着十六进制编辑器里的字节心里想的只有一个问题这堆数字到底哪一截是坐标后来师傅甩给我一份道头字段表我才明白SEGY能通吃地震行业几十年靠的就是把每道地震数据的“元信息”固定地塞进了那240个字节里。你不需要背下所有字段但你必须知道和坐标定位、可视化的那几个位置长在哪。这篇文章就是给准备做地震数据可视化的人看的。无论你是处理叠前道集还是画叠后剖面的测线位置只要想把SEGY里的地震道放到它该在的空间位置上道头字段就是你的第一道关口。我会把最核心的字段、坐标换算的野路子、以及我踩过的坑一次说清。1. SEGY到底是什么先说清这个“地震数据快递箱”1.1 从磁带时代走来的格式标准SEGY全称是SEG-Y由Society of Exploration GeophysicistsSEG提出最初是为了把地震数据从采集车上卸下来装到磁带上方便不同厂商的机器互相读取。这个标准从七十年代用到现在虽然中间出过Rev0、Rev1这些版本但整体骨架没变一个SEGY文件就像一车快递最前面是总运单后面是一箱一箱的货物每一箱上贴着一张标签。落实到文件结构上就是三段3200字节的文本卷头Textual Header存放工区名、施工日期、采样率等说明性文字本质上是ASCII文本。400字节的二进制卷头Binary Header用固定字节记录采样点数、采样间隔、数据格式代码等机器要用的参数。地震道数据Trace Data一条一条排下来每条道又分成240字节的道头Trace Header和紧随其后的振幅采样数据。道头就是那张“货物标签”它告诉你这一箱东西是从哪里来的、往哪运的。SEG-Y Rev1标准的道头固定是240字节所有关键信息都按字节偏移写死。你做可视化时等于是先读标签再决定把这箱货放哪。1.2 道头字段是可视化绕不开的“户籍档案”做地震数据可视化最基础的需求就是把每一道放到平面图上。要放对位置至少需要知道这一道属于哪条测线、哪个CDP、炮点坐标是多少、检波点坐标是多少。这些信息全都写在道头里。如果你对道头字段一知半解看到一堆坐标数据也不知道是什么单位、需不需要标量换算。更麻烦的是不同处理系统在写SEGY时可能会用不同扩展槽记录GPS坐标这时候不懂道头就完全没法下手。换句话说道头字段就是每道地震数据的户籍档案不读它你手里只有一堆顺序排列的振幅值没有空间属性可视化就无从谈起。2. 道头字段全解析5分钟掌握关键字节的“密码本”2.1 道头核心字段与字节位置速查表下面这张表是从SEG-Y Rev1里摘出来的常用字段字节位置按1开始计数和大多数SEGY文档保持一致。做可视化和坐标定位重点看标了加粗的坐标相关字段。字节位置字段名称说明1-4Trace sequence number within line测线内道序号5-8Trace sequence number within file文件内道序号9-12Original field record number原始记录号通常是炮号13-16Trace number within field record原始记录内道号17-20Energy source point number能量源点号桩号21-24CDP ensemble numberCDP号/共深度点号25-28Trace number within CDP ensembleCDP内道号29-30Trace identification code道识别码1地震道2死道3哑道31-32Vertically summed traces垂直叠加道数33-34Horizontally stacked traces水平叠加道数35-36Data use数据用途37-40Distance from source to receiver group偏移距41-44Receiver group elevation检波点高程45-48Surface elevation at source炮点地表高程49-52Source depth below surface炮点深度53-56Datum elevation at receiver group检波点基准面高程57-60Datum elevation at source炮点基准面高程61-64Water depth at source炮点水深65-68Water depth at group检波点水深69-70Scalar for elevations and depths高程/深度标量71-72Scalar for coordinates坐标标量73-76Source X coordinate炮点X坐标77-80Source Y coordinate炮点Y坐标81-84Group X coordinate检波点X坐标85-88Group Y coordinate检波点Y坐标89-90Coordinate units坐标单位1米2英尺91-92Weathering velocity风化层速度93-94Subweathering velocity风化层下速度95-96Uphole time at source炮点井口时间97-98Uphole time at group检波点井口时间你可能会发现道头里各种高程字段也带标量但坐标标量是独立在71-72字节的千万别拿高程标量去换算坐标。坐标单位字段在89-90字节很多老数据会漏填默认按米处理问题不大但严谨起见还是先看一眼。2.2 字节序与整型陷阱大端序下的常见翻车现场道头字段绝大多数是4字节或2字节的带符号整数。SEGY沿用了磁带时代大端序big-endian的习惯也就是高字节在前。如果你用Python的struct去读格式符要用i而不是i否则读出来的数字会完全不对。举个例子坐标字段在字节73-76假设原始数据是十六进制00 00 B0 74。按大端规则读取得到45300按小端规则读变成1953005568完全不是一个量级。当你发现读出来的坐标是几亿甚至上百亿时第一反应应该是字节序错了而不是去怀疑坐标精度。还有一个和道头无关、但处理SEGY时一定会遇到的坑振幅数据很多是IBM浮点格式不是平时CPU直接支持的IEEE float。IBM浮点的二进制位分配和IEEE完全不同如果你直接把数据块按float32解析波形就是一团噪点。本文重点在道头不展开振幅格式但你心里要有数后面做波形可视化时一定会碰到。3. 坐标定位实操指南把道头里的坐标变成地图上的点3.1 坐标标量换算正负号决定乘还是除道头里的坐标值并不总是真实坐标。为了节省存储空间早期磁带会把坐标放大或缩小后存成整数等到使用时再靠“坐标标量”还原。标量在道头71-72字节是一个2字节带符号整数。规则是这样的标量大于0真实值 道头值 × 标量标量小于0真实值 道头值 ÷ |标量|标量等于0不做变换这个规则确实和直觉反过来很多人第一次用会搞错。比如道头里读到的坐标是4529400标量是-100真实坐标就是45294.00米如果标量是100则真实坐标是452940000那就非常离谱了。所以实际数据里标量通常都是负数表示要除以10、100或1000来还原小数位。写代码时标量处理应该封装成一个函数不要到处手写乘除否则很容易在某个分支上漏掉符号判断。同样的逻辑也适用于高程标量但高程标量和坐标标量是不同的字段要分开读。3.2 炮点坐标 vs 检波点坐标 vs CDP坐标道头里同时写了炮点坐标73-80字节和检波点坐标81-88字节。这两个坐标是叠前资料的“元数据核心”。做可视化之前先想清楚自己要画什么画炮点在平面上的分布用Source X/Y能看出震源走线是否均匀。画检波点分布用Group X/Y能看出接收排列的覆盖范围。要画CDP点位置更常见的做法是用炮点和检波点的中点MX (SX GX) / 2MY (SY GY) / 2。这个中点其实是道集的中心坐标也是叠后道数据最常用的空间位置。如果你手里是叠后数据道头里通常不会填炮点检波点坐标而是在CDP号附近放CDP坐标。不同处理中心写的位置不一样有的把CDP坐标放在73-76和77-80有的放在道头扩展槽。遇到坐标全是0的情况别急着怀疑文件坏了先检查是不是叠后数据再查处理报告。3.3 坐标系与投影UTM、经纬度不要混着画道头坐标最常见的是UTM投影坐标和经纬度坐标两种坐标混着画图出来的测线形状会被严重拉伸。我的经验是先看数值量级判断坐标系UTM东向通常在六到七位数比如500000左右北向也在六位数经纬度则是几十到一百多带小数位。如果是经纬度坐标X对应经度Y对应纬度画图时必须设置等纵横比。用Matplotlib画散点加一行ax.set_aspect(equal)就行不然纬度一度和经度一度的屏幕长度不同测线会像被压扁一样。如果是UTM坐标X就是东向EastingY就是北向Northing同样要设置等纵横比保证测线方位角不变形。如果在同一工区里有多个SEGY文件有些用UTM有些用经纬度最好统一投影到同一个坐标系再画。最省事的办法是用pyproj把经纬度转成UTM或者把UTM转成经纬度统一之后再合并画图。4. 用Python快速实现SEGY道头解析与可视化4.1 环境准备安装segyio和绘图库解析SEGY的Python库不少但最省心的是segyio它由Statoil现Equinor维护对道头字段的封装很完整底层直接用C实现读大文件也快。另一个库是obspy偏地震学研究处理反射地震勘探数据不如segyio顺手。安装命令很简单pip install segyio matplotlib numpy如果你在Windows上安装segyio遇到问题多半是缺少编译环境。建议优先从PyPI下载预编译的wheel或者先装Anaconda再用conda install -c conda-forge segyio。4.2 读取道头坐标并绘制测线图的完整代码示例下面这段代码是我目前在用的“SEGY坐标快检”脚本专门用来在拿到新工区数据后5分钟之内画出一张空间分布图。它读取所有道的炮点、检波点坐标换算标量然后计算中点坐标并绘制散点图。import segyio import numpy as np import matplotlib.pyplot as plt def load_coordinates(segy_path): with segyio.open(segy_path, r, ignore_geometryTrue) as f: sx np.asarray(f.attributes(segyio.TraceField.SourceX)[:], dtypenp.float64) sy np.asarray(f.attributes(segyio.TraceField.SourceY)[:], dtypenp.float64) gx np.asarray(f.attributes(segyio.TraceField.GroupX)[:], dtypenp.float64) gy np.asarray(f.attributes(segyio.TraceField.GroupY)[:], dtypenp.float64) scalar np.asarray(f.attributes(segyio.TraceField.SourceGroupScalar)[:], dtypenp.float64) units np.asarray(f.attributes(segyio.TraceField.CoordinateUnits)[:], dtypenp.int64) return sx, sy, gx, gy, scalar, units def apply_scalar(values, scalar): out values.copy() pos scalar 0 neg scalar 0 out[pos] values[pos] * scalar[pos] out[neg] values[neg] / np.abs(scalar[neg]) return out def plot_coordinates(sx, sy, gx, gy): # 中点坐标也就是叠后CDP点的空间位置 mx (sx gx) / 2.0 my (sy gy) / 2.0 plt.figure(figsize(10, 9)) plt.scatter(mx, my, s4, cblue, alpha0.6, labelMidpoint) plt.scatter(sx, sy, s1, cred, alpha0.3, labelSource) ax plt.gca() ax.set_aspect(equal) ax.set_xlabel(X coordinate) ax.set_ylabel(Y coordinate) ax.set_title(SEGY Trace Position Overview) ax.legend() ax.grid(True, linestyle--, alpha0.5) plt.show()调用方式sx, sy, gx, gy, scalar, units load_coordinates(line.sgy) sx apply_scalar(sx, scalar) sy apply_scalar(sy, scalar) gx apply_scalar(gx, scalar) gy apply_scalar(gy, scalar) # 如果坐标单位是英尺转成米 # if np.all(units 2): # sx * 0.3048 # sy * 0.3048 # gx * 0.3048 # gy * 0.3048 plot_coordinates(sx, sy, gx, gy)ignore_geometryTrue这个参数很关键。很多文件没有写完整的测线几何信息如果不加这一项segyio打开时会提示找不到几何结构甚至报错。加上之后它会忽略测线几何只按道顺序读取对我们解析道头坐标完全够用。f.attributes(字段)[:]会一次性返回所有道的该字段值内部是C实现比在Python里逐道循环快上百倍。数据量几万道时几乎没有区别但几十万道时就能明显感到差别了。4.3 如果segyio不好用手写一个极简道头坐标读取器有时候环境受限装不了segyio或者你只是想快速验证一个文件的字节布局手写一个极简读取器也能应急。下面这个函数只读前几道的道头坐标用struct按大端序解包import struct def peek_coords_from_segy(filename, num_traces3, nsamp1000, sample_format1): 简易SEGY道头坐标读取器。 nsamp每道采样点数需要从二进制卷头取得这里先按已知值填。 sample_format1为IBM浮点每样本4字节2为4字节整数5为IEEE浮点。 sample_bytes {1: 4, 2: 4, 3: 2, 5: 4, 6: 8} data_len nsamp * sample_bytes.get(sample_format, 4) with open(filename, rb) as f: f.seek(3600) # 跳过3200字节文本头 400字节二进制头 for t in range(num_traces): header f.read(240) sx struct.unpack(i, header[72:76])[0] sy struct.unpack(i, header[76:80])[0] gx struct.unpack(i, header[80:84])[0] gy struct.unpack(i, header[84:88])[0] scalar struct.unpack(h, header[70:72])[0] print(fTrace {t}: SX{sx}, SY{sy}, GX{gx}, GY{gy}, Scalar{scalar}) f.seek(data_len, 1) # 跳过振幅数据注意这里header[72:76]对应的是0-based偏移72也就是文档里的字节73-76header[70:72]对应坐标标量。如果你不知道每道的采样点数和采样格式需要先读二进制卷头里的相应位置这里只是展示原理。生产环境请老老实实用segyio手写方案只适合应急或教学。5. 常见问题与排查技巧实录5.1 坐标全是0或者“天文数字”怎么办坐标全是0最常见的原因是文件是叠后数据处理系统没往标准坐标槽位写坐标坐标可能在CDP号扩展区或辅助道里。另一个原因是你读的字段不对比如读成了高程。坐标是“天文数字”九成是标量没换算。先用一个简单的办法验证把前三个道头的SX、SY、GX、GY原始值和标量打出来手算一下看看除以标量后的数值是不是落在工区合理范围。如果除以标量后仍然不合理再怀疑字节序和坐标单位问题。5.2 道头读取出来的值像负数、乱码如果你是用手写方式解析读出来的坐标不定期出现负数先检查是不是把有符号整型当成了无符号整型或者字节序反了。struct.unpack(i)读4字节有符号大端整数是标准做法。道头里有些字段本来就是有符号的比如高程可以有负值但坐标一般为正。如果你用segyio读出来的还是乱那就不是代码问题而是文件真的被写得不规范。可以用十六进制编辑器直接看道头字节比如SX字段在文件第3600t×(240data_len)72字节的位置看一眼是不是像正常的整数。5.3 画出来的测线方向不对或整体翻转坐标颠倒是很常见的坑尤其是野外采集的SEGY里Source和Group的X/Y定义可能被某些采集系统写反。解决办法很简单把炮点坐标和检波点坐标交换重新计算中点并画图对比哪一版和已知的测线位置吻合。另外有些数据源的Y方向是南而不是北导致整个测线镜像。这种情况多半是坐标系定义差异比如用了南半球UTM或自定义投影。拿到数据时如果能问一句处理方“坐标是哪个UTM带基准面是WGS84还是本地基准面”就能省下大量排查时间。5.4 扩展道头里藏着GPS信息SEGY Rev1允许在道头之后插入扩展文本头但这和坐标关系不大。真正麻烦的是部分采集设备会把高精度GPS信息放在道头的“用户槽”里而不是标准73-88的位置。如果你发现自己熟悉的所有标准字段都读不出合理坐标可以打一个前几道的十六进制dump找有没有类似X:、Y:这种ASCII标记或者和采集队确认道头字节分配表。我有一年在处理两套衔接工区数据时就遇到过一个文件用标准坐标字段另一个用自定义槽位。最后是靠两者的叠加覆盖图才互相印证出来位置。这件事给我的教训是永远不要默认SEGY文件是“标准”的只有打印出来的道头才是真相。6. 最后再分享一个实用习惯每次拿到新的SEGY文件我的固定流程是先打印前三个道头的关键坐标字段然后画炮点、检波点、中点的散点分布图确认空间形态没问题再做后续的波形显示或属性分析。这个过程通常不超过5分钟但能帮我避开至少一半的“数据看起来是好的、画出来全错”的尴尬场景。道头解析这门手艺说白了就是“按图索骥”。SEGY格式耐用了这么多年核心就是它把所有元数据固定在了一套可预期的方式里。你只要愿意花几分钟把道头字段表过一遍再配合几个小函数后面不管哪一年、哪个工区来的数据都能快速定位到正确位置。如果后续有时间我再写一篇关于SEGY二进制卷头里采样率、数据格式和坐标网格之间关系的文章那个问题也很有意思。