上个月帮朋友处理一批三维采集数据文件7.4GB诉求很简单“把工区里的炮点画到图上我要看覆盖次数。”我当时按老规矩先拆道头结果第一炮点的坐标直接落在海里。排查了半天才发现不是坐标写错了而是道头里scalar标定因子是负的我没做除法把存储值当成真实米数用了。从那之后我养成了一个习惯任何SEGY文件到手先读道头、先算标定、先看坐标分布再谈画图。这篇文章就围绕SEGY道头字段展开重点解决三个问题240字节的道头到底存了什么、坐标字段怎么正确解读、如何用Python快速把几百万道数据变成可视化图件。适合刚接触地震数据可视化的开发、处理员、研究生也适合那些被“坐标漂移”折磨过但一直没找到原因的从业者。5分钟内你至少能知道该读哪些字节、该踩哪些坑。1. SEGY不是玄学先搞清楚“谁在描述谁”1.1 一个从盒式磁带时代走来的格式为什么今天还在用SEG Y是1975年SEG协会发布的勘探地震数据记录格式rev0、rev1、rev2一路演进核心结构四十多年没变过。整个文件可以理解成一个快递盒文件头是快递单道头是每个包裹上的说明贴纸数据体才是包裹里的实物。3600字节文本头人可读的作业信息包含工区名、施工日期、处理参数等。3200字节二进制头记录了采样点数、采样率、数据格式码等机器必须知道的信息。每道240字节道头描述这一道地震数据的属性道号、炮点、CDP、坐标、偏移距都在这里。之后是数据体每道由若干采样点组成浮点数或整型视格式码而定。野外采集回来的原始炮集、处理系统导出的成果剖面、叠前时间偏移后的道集绝大多数都以SEGY形式交付。可视化之前不读道头你根本不知道这几十万道在地面上是什么位置画出来等于盲画。1.2 可视化前必须先回答的三个问题拿到一个SEGY文件第一件事不是急着读数据体而是问自己三个问题我的数据是二维测线还是三维工区二维和三维的坐标字段解读方式完全不同。道头里存的是网格坐标inline/crossline还是真实地理坐标米或经纬度坐标标定因子scalar是几是正数还是负数这三个问题的答案全部藏在道头字段里。很多人用绘图软件直接load二进制当灰度图或者调用现成库读SEGY但完全不看道头最后画出来的剖面翻转、坐标漂移还以为是数据坏了其实只是没读字段。2. 240字节的道头里面到底装了什么2.1 读道头之前先看文件级信息SEGY文件开头的3600字节是文本头按40行乘80字符组织。规范上推荐用EBCDIC编码不过现在多数处理系统也写ASCII。用Python读的时候要先判断内容样式不然出来全是乱码。紧接着是3200字节的二进制头采样点数存在3221-3222字节采样间隔存在3217-3218字节数据格式码存在3225-3226字节。这三个参数是后续解析数据体的关键只有知道采样点数和格式码才能算出“一道到底占多少字节”。举个例子格式码为1代表IBM浮点4字节为5代表IEEE浮点4字节为3代表16位整型。一道数据体长度 采样点数 × 单样点字节数。这个数直接决定你能否用内存映射方式快速跳过数据体、只读道头。2.2 道头字段核心表按字节号对号入座240字节道头从字节1开始编号。注意SEGY规范从1开始而Python切片从0开始读取时字段起始字节需要减1。下面是我实际项目里最常用的一组字段。字节起始字节数字段名说明常见取值/单位14道序号文件内从1递增整数54原始道号野外记录时的道号整数94道集号CDP号或炮集号整数134道集内道号道集内从1开始整数174道识别码1地震道2哑道3空道4时间信号整数214偏移距炮点到检波点的距离米或英尺292道距/增量相邻道的距离0.01米或英尺712坐标标定因子所有坐标都要用它换算正负整数732采样点数该道样点数整数772采样间隔微秒为单位整数812增益类型1fixed等整数1034道源号震源点号整数1154炮点X坐标震源X存储值1194炮点Y坐标震源Y存储值1334接收点X坐标检波器X存储值1374接收点Y坐标检波器Y存储值1814道集X坐标CDP或道集X存储值1854道集Y坐标CDP或道集Y存储值1892道集坐标标定rev1新增覆盖71-72正负整数这里有个非常容易踩的点SEGY标准并没有规定X必须是东向、Y必须是北向。国内很多工区资料习惯把X当北向、Y当东向而国外软件默认XEasting、YNorthing。坐标反了的案例我见过不下十次尤其是从国外软件导出的成果直接套国内坐标系的时候。2.3 道头里三个容易被忽略的“隐藏信息”第一是17-20字节的道识别码。可视化前必须按它过滤哑道和空道混进坐标点会造成散点图出现大量噪声。第二是71-72字节的坐标标定因子。第三是二进制头格式码。这三个字段不确认后面任何坐标计算都没意义。3. 坐标定位的真正麻烦SP/网格/投影三道关3.1 你手里到底有几个坐标来源二维数据通常很简单181-184字节存CDP坐标X185-188存CDP坐标Y115-118和119-122存炮点坐标。画测线位置图用CDP坐标即可。三维规则工区不同很多叠后成果三维数据道头里存的是CDP的平面坐标但部分采集系统导出的SEGY道头只写了网格线号inline/crossline编号真实坐标必须由工区起算点和线距推算出来。我见过最坑的一种情况是道头115-118字节存的是类似1001、1002、1003的小整数看着像坐标实际是crossline编号。如果你照着这个值直接散点画图得到的是一条直线而不是工区平面。3.2 scalar标定因子一半坐标乌龙案的元凶坐标字段在SEGY里通常存的是整数而不是浮点数真实坐标需要靠71-72字节的scalar来还原如果scalar 0真实坐标 存储值 × scalar。如果scalar 0真实坐标 存储值 ÷ |scalar|。如果scalar 0说明坐标字段本身就是真实值不需要换算。举个例子存储值是18495008scalar -100真实坐标就是184950.08米如果漏做除法直接把这个数当米用点位会偏到一百多公里外直接飞出工区。代码里最简单的处理方式if scalar 0: coord_real stored_value * scalar elif scalar 0: coord_real stored_value / abs(scalar) else: coord_real stored_value3.3 网格坐标转经纬度不能硬转如果道头里存的是高斯平面坐标或UTM坐标想叠加到底图上必须先做投影转换。用pyproj可以一行完成from pyproj import Transformer # 示例UTM 50N 转 WGS84 经纬度 transformer Transformer.from_crs(EPSG:32650, EPSG:4326, always_xyTrue) lon, lat transformer.transform(x_m, y_m)但前提是你必须知道数据用的什么投影、哪个带号。不知道的话别硬猜老老实实翻原始资料或问数据提供方。判断带号有个土办法看东向坐标的前两位比如Y值开头是50多半是UTM 50带如果是带中央经线的500000左右要结合当地经度反推带号。3.4 三维数据里常见的“伪坐标”某些处理系统导出的SEGY道头坐标字段里写的是inline/crossline编号而不是真实地面坐标。怎么快速判断取连续几道道头看181-184字节的值是不是在递增的小整数比如1001、1002、1003。如果是基本可以确定是网格号。此时要拿到工区起点坐标、inline方向和crossline方向的线距自己算真实平面坐标real_x origin_x (inline - inline_start) * inline_spacing real_y origin_y (crossline - crossline_start) * crossline_spacing注意inline和crossline哪个对应行、哪个对应列不同系统定义不一样这也是一个大坑。4. 手把手Python解析SEGY道头并画出工区图4.1 选segyio还是手写struct正式项目我首选segyio它是社区标准库底层用C实现遍历几十万道道头速度很快。但segyio的字段映射固定遇到非标厂商把坐标写在备用字节时就得手写struct或numpy解析。所以两条路都得会。4.2 完整示例代码读道头画散点下面这段代码读取SEGY文件的CDP坐标应用scalar标定过滤哑道最终用matplotlib画出工区炮点/CDP点分布图。import segyio import numpy as np import matplotlib.pyplot as plt path survey.sgy with segyio.open(path, r, ignore_geometryTrue) as f: trcount f.tracecount sample_count f.bin[segyio.BinField.Samples] scalar f.header[0][segyio.TraceField.SourceGroupScalar] xs np.empty(trcount, dtypenp.float64) ys np.empty(trcount, dtypenp.float64) codes np.empty(trcount, dtypenp.int32) for i in range(trcount): h f.header[i] codes[i] h[segyio.TraceField.TraceIdentificationCode] xs[i] h[segyio.TraceField.CDP_X] ys[i] h[segyio.TraceField.CDP_Y] # 应用标定因子 if scalar 0: xs xs * scalar ys ys * scalar elif scalar 0: xs xs / abs(scalar) ys ys / abs(scalar) # 过滤地震道和非地震道保留1地震数据 mask codes 1 # 坐标全为0或明显异常的点也过滤掉 mask (xs ! 0) (ys ! 0) fig, ax plt.subplots(figsize(10, 8)) ax.scatter(xs[mask], ys[mask], s1, marker., alpha0.5) ax.set_aspect(equal) ax.set_xlabel(X (m)) ax.set_ylabel(Y (m)) ax.set_title(SEGY CDP Location Map) plt.show()这段代码对于二维测线和三维工区都适用。二维数据画出来是一条测线三维数据画出来是一个面状分布。如果画出来是横七竖八的线段说明道头坐标可能没写或者写的是网格号。4.3 手写numpy版复杂场景的兜底方案遇到segyio不支持的字段布局时用numpy结构化数组自己解。思路是先读二进制头的采样点数和格式码算出一道总字节数再批量解析道头import numpy as np def read_trace_coords(path, trace_bytes, trace_count): dtype np.dtype([ (seq, i4, (1,)), # 1-4 (trace, i4, (1,)), # 5-8 (cdp, i4, (1,)), # 9-12 (code, i4, (1,)), # 17-20 (scalar, i2, (1,)), # 71-72 (sx, i4, (1,)), # 115-118 (sy, i4, (1,)), # 119-122 (cdp_x, i4, (1,)), # 181-184 (cdp_y, i4, (1,)), # 185-188 ]) # 仅读取每个道的道头部分跳过数据体 # 使用mmap避免把整个大文件读进内存 arr np.memmap(path, dtypenp.uint8, moder) trace_start 3600 3200 raw arr[trace_start: trace_start trace_count * trace_bytes] raw raw.reshape(trace_count, trace_bytes) headers raw[:, :240] # 按需复制出字段 cdp_x_bytes headers[:, 180:184].copy() cdp_y_bytes headers[:, 184:188].copy() cdp_x cdp_x_bytes.view(np.int32).ravel() cdp_y cdp_y_bytes.view(np.int32).ravel() return cdp_x, cdp_y注意字节切片要按0基索引字段起始字节减1。这个方案只copy需要的部分不会一次性把数据体全部导入内存适合几个GB甚至上百GB的文件。5. 坐标错误的三类真实翻车现场5.1 现场一坐标全在同一位置或者全是0排查链路先看71-72字节scalar是否为0或异常值。确认你读的是CDP坐标还是炮点坐标。有的数据只写了炮点坐标181-188全为0。检查有没有扩展道头。rev1允许在3600字节文本头之后、正式道数据之前插入扩展文本头某些系统写的扩展文本头会干扰偏移计算。如果以上都正常把240字节道头的前64个字节用十六进制dump出来肉眼对比几道是否一致判断是否整道头都没写。这个翻车现场最常见的结论是数据提供方压根没在道头里写坐标后续需要靠SPS文件或者单独的位置表来关联。5.2 现场二坐标和图件对不上差了一点点或直接翻转坐标偏移几个公里、几十公里常见原因有三个scalar标定漏做、X/Y弄反、投影带搞错。判断方法很实用取数据中两个相邻CDP点算平面距离看是否约等于理论道距。比如理论道距是25米算出来是25000米那基本是标定因子漏了算出来距离对但方位角反了优先怀疑X/Y互换。X/Y互换还有个特征画出来的测线走向和实际工区的长方形边界成90度旋转关系。很多系统里“X北、Y东”和“X东、Y北”两种习惯并存读取后一定要和工区范围描述核对。5.3 现场三三维工区画出来像被撕碎的网格三维数据可视化时出现“锯齿形”和“撕裂感”通常是因为道号顺序和文件存储顺序不一致或者采集时炮点按非规则顺序排列。解决办法不是改绘图代码而是用CDP_X和CDP_Y重新排序或者先用inline/crossline编号排序再画。还有一种情况是dummy道占了位置但道识别码没写对。处理前先按17-20字节过滤只看code1的地震道散点图立刻干净很多。6. 大文件实操怎么从上百GB的SEGY里几秒只读道头6.1 用np.memmap配合步长抽取上百GB的SEGY文件不可能整体读进内存。最实用的是内存映射加只读道头切片。二进制头里的采样点数和格式码先拿到的然后算出每道总字节数用下面这套逻辑import numpy as np def memmap_trace_headers(path, trace_start, trace_bytes, trace_count): arr np.memmap(path, dtypenp.uint8, moder) data arr[trace_start: trace_start trace_count * trace_bytes] data data.reshape(trace_count, trace_bytes) # 取出每道前240字节作为道头 headers data[:, :240] return headers配合5%抽稀做质控图几百GB的数据几秒就能出轮廓。实际使用中我一般先map后只取坐标字段的字节段减少复制量。6.2 两步抽稀法先轮廓后细节直接渲染几百万个散点会让matplotlib卡死。我自己习惯分两步第一次只取前1万道或均匀抽稀5%看整体工区轮廓和坐标范围第二步再全量读取所有道头坐标输出成CSV或GeoJSON交给前端Leaflet或Cesium显示。工区轮廓没问题后再谈逐道渲染。6.3 输出GeoJSON给Web端处理完的坐标最好直接落成GeoJSON后续无论是Web可视化还是GIS叠加都很方便import json features [] for x, y in zip(xs[mask], ys[mask]): features.append({ type: Feature, geometry: { type: Point, coordinates: [float(x), float(y)] }, properties: {} }) geojson {type: FeatureCollection, features: features} with open(cdp_points.geojson, w) as fp: json.dump(geojson, fp)如果后面还要和地形数据、DEM叠合坐标系必须统一否则经纬度对了、高程差几十米的情况也会出现。最后再分享一个我自己的习惯拿到任何SEGY数据我先不画图先做一张质控图——把道头里的道集号按顺序着色看颜色是否连续。跳跃太大说明文件边界或道序有问题这时候画的任何“漂亮图”都是不可信的。道头解析最好做成一个脚本入参是文件路径、坐标字段选择、是否做投影转换输出GeoJSON加PNG这样以后每来一份新数据都能在5分钟内完成定位质控。这个习惯帮我躲过很多次“坐标漂移”的坑建议你也留一份。
SEGY道头解析与坐标标定：Python读取道头并绘制工区坐标图
上个月帮朋友处理一批三维采集数据文件7.4GB诉求很简单“把工区里的炮点画到图上我要看覆盖次数。”我当时按老规矩先拆道头结果第一炮点的坐标直接落在海里。排查了半天才发现不是坐标写错了而是道头里scalar标定因子是负的我没做除法把存储值当成真实米数用了。从那之后我养成了一个习惯任何SEGY文件到手先读道头、先算标定、先看坐标分布再谈画图。这篇文章就围绕SEGY道头字段展开重点解决三个问题240字节的道头到底存了什么、坐标字段怎么正确解读、如何用Python快速把几百万道数据变成可视化图件。适合刚接触地震数据可视化的开发、处理员、研究生也适合那些被“坐标漂移”折磨过但一直没找到原因的从业者。5分钟内你至少能知道该读哪些字节、该踩哪些坑。1. SEGY不是玄学先搞清楚“谁在描述谁”1.1 一个从盒式磁带时代走来的格式为什么今天还在用SEG Y是1975年SEG协会发布的勘探地震数据记录格式rev0、rev1、rev2一路演进核心结构四十多年没变过。整个文件可以理解成一个快递盒文件头是快递单道头是每个包裹上的说明贴纸数据体才是包裹里的实物。3600字节文本头人可读的作业信息包含工区名、施工日期、处理参数等。3200字节二进制头记录了采样点数、采样率、数据格式码等机器必须知道的信息。每道240字节道头描述这一道地震数据的属性道号、炮点、CDP、坐标、偏移距都在这里。之后是数据体每道由若干采样点组成浮点数或整型视格式码而定。野外采集回来的原始炮集、处理系统导出的成果剖面、叠前时间偏移后的道集绝大多数都以SEGY形式交付。可视化之前不读道头你根本不知道这几十万道在地面上是什么位置画出来等于盲画。1.2 可视化前必须先回答的三个问题拿到一个SEGY文件第一件事不是急着读数据体而是问自己三个问题我的数据是二维测线还是三维工区二维和三维的坐标字段解读方式完全不同。道头里存的是网格坐标inline/crossline还是真实地理坐标米或经纬度坐标标定因子scalar是几是正数还是负数这三个问题的答案全部藏在道头字段里。很多人用绘图软件直接load二进制当灰度图或者调用现成库读SEGY但完全不看道头最后画出来的剖面翻转、坐标漂移还以为是数据坏了其实只是没读字段。2. 240字节的道头里面到底装了什么2.1 读道头之前先看文件级信息SEGY文件开头的3600字节是文本头按40行乘80字符组织。规范上推荐用EBCDIC编码不过现在多数处理系统也写ASCII。用Python读的时候要先判断内容样式不然出来全是乱码。紧接着是3200字节的二进制头采样点数存在3221-3222字节采样间隔存在3217-3218字节数据格式码存在3225-3226字节。这三个参数是后续解析数据体的关键只有知道采样点数和格式码才能算出“一道到底占多少字节”。举个例子格式码为1代表IBM浮点4字节为5代表IEEE浮点4字节为3代表16位整型。一道数据体长度 采样点数 × 单样点字节数。这个数直接决定你能否用内存映射方式快速跳过数据体、只读道头。2.2 道头字段核心表按字节号对号入座240字节道头从字节1开始编号。注意SEGY规范从1开始而Python切片从0开始读取时字段起始字节需要减1。下面是我实际项目里最常用的一组字段。字节起始字节数字段名说明常见取值/单位14道序号文件内从1递增整数54原始道号野外记录时的道号整数94道集号CDP号或炮集号整数134道集内道号道集内从1开始整数174道识别码1地震道2哑道3空道4时间信号整数214偏移距炮点到检波点的距离米或英尺292道距/增量相邻道的距离0.01米或英尺712坐标标定因子所有坐标都要用它换算正负整数732采样点数该道样点数整数772采样间隔微秒为单位整数812增益类型1fixed等整数1034道源号震源点号整数1154炮点X坐标震源X存储值1194炮点Y坐标震源Y存储值1334接收点X坐标检波器X存储值1374接收点Y坐标检波器Y存储值1814道集X坐标CDP或道集X存储值1854道集Y坐标CDP或道集Y存储值1892道集坐标标定rev1新增覆盖71-72正负整数这里有个非常容易踩的点SEGY标准并没有规定X必须是东向、Y必须是北向。国内很多工区资料习惯把X当北向、Y当东向而国外软件默认XEasting、YNorthing。坐标反了的案例我见过不下十次尤其是从国外软件导出的成果直接套国内坐标系的时候。2.3 道头里三个容易被忽略的“隐藏信息”第一是17-20字节的道识别码。可视化前必须按它过滤哑道和空道混进坐标点会造成散点图出现大量噪声。第二是71-72字节的坐标标定因子。第三是二进制头格式码。这三个字段不确认后面任何坐标计算都没意义。3. 坐标定位的真正麻烦SP/网格/投影三道关3.1 你手里到底有几个坐标来源二维数据通常很简单181-184字节存CDP坐标X185-188存CDP坐标Y115-118和119-122存炮点坐标。画测线位置图用CDP坐标即可。三维规则工区不同很多叠后成果三维数据道头里存的是CDP的平面坐标但部分采集系统导出的SEGY道头只写了网格线号inline/crossline编号真实坐标必须由工区起算点和线距推算出来。我见过最坑的一种情况是道头115-118字节存的是类似1001、1002、1003的小整数看着像坐标实际是crossline编号。如果你照着这个值直接散点画图得到的是一条直线而不是工区平面。3.2 scalar标定因子一半坐标乌龙案的元凶坐标字段在SEGY里通常存的是整数而不是浮点数真实坐标需要靠71-72字节的scalar来还原如果scalar 0真实坐标 存储值 × scalar。如果scalar 0真实坐标 存储值 ÷ |scalar|。如果scalar 0说明坐标字段本身就是真实值不需要换算。举个例子存储值是18495008scalar -100真实坐标就是184950.08米如果漏做除法直接把这个数当米用点位会偏到一百多公里外直接飞出工区。代码里最简单的处理方式if scalar 0: coord_real stored_value * scalar elif scalar 0: coord_real stored_value / abs(scalar) else: coord_real stored_value3.3 网格坐标转经纬度不能硬转如果道头里存的是高斯平面坐标或UTM坐标想叠加到底图上必须先做投影转换。用pyproj可以一行完成from pyproj import Transformer # 示例UTM 50N 转 WGS84 经纬度 transformer Transformer.from_crs(EPSG:32650, EPSG:4326, always_xyTrue) lon, lat transformer.transform(x_m, y_m)但前提是你必须知道数据用的什么投影、哪个带号。不知道的话别硬猜老老实实翻原始资料或问数据提供方。判断带号有个土办法看东向坐标的前两位比如Y值开头是50多半是UTM 50带如果是带中央经线的500000左右要结合当地经度反推带号。3.4 三维数据里常见的“伪坐标”某些处理系统导出的SEGY道头坐标字段里写的是inline/crossline编号而不是真实地面坐标。怎么快速判断取连续几道道头看181-184字节的值是不是在递增的小整数比如1001、1002、1003。如果是基本可以确定是网格号。此时要拿到工区起点坐标、inline方向和crossline方向的线距自己算真实平面坐标real_x origin_x (inline - inline_start) * inline_spacing real_y origin_y (crossline - crossline_start) * crossline_spacing注意inline和crossline哪个对应行、哪个对应列不同系统定义不一样这也是一个大坑。4. 手把手Python解析SEGY道头并画出工区图4.1 选segyio还是手写struct正式项目我首选segyio它是社区标准库底层用C实现遍历几十万道道头速度很快。但segyio的字段映射固定遇到非标厂商把坐标写在备用字节时就得手写struct或numpy解析。所以两条路都得会。4.2 完整示例代码读道头画散点下面这段代码读取SEGY文件的CDP坐标应用scalar标定过滤哑道最终用matplotlib画出工区炮点/CDP点分布图。import segyio import numpy as np import matplotlib.pyplot as plt path survey.sgy with segyio.open(path, r, ignore_geometryTrue) as f: trcount f.tracecount sample_count f.bin[segyio.BinField.Samples] scalar f.header[0][segyio.TraceField.SourceGroupScalar] xs np.empty(trcount, dtypenp.float64) ys np.empty(trcount, dtypenp.float64) codes np.empty(trcount, dtypenp.int32) for i in range(trcount): h f.header[i] codes[i] h[segyio.TraceField.TraceIdentificationCode] xs[i] h[segyio.TraceField.CDP_X] ys[i] h[segyio.TraceField.CDP_Y] # 应用标定因子 if scalar 0: xs xs * scalar ys ys * scalar elif scalar 0: xs xs / abs(scalar) ys ys / abs(scalar) # 过滤地震道和非地震道保留1地震数据 mask codes 1 # 坐标全为0或明显异常的点也过滤掉 mask (xs ! 0) (ys ! 0) fig, ax plt.subplots(figsize(10, 8)) ax.scatter(xs[mask], ys[mask], s1, marker., alpha0.5) ax.set_aspect(equal) ax.set_xlabel(X (m)) ax.set_ylabel(Y (m)) ax.set_title(SEGY CDP Location Map) plt.show()这段代码对于二维测线和三维工区都适用。二维数据画出来是一条测线三维数据画出来是一个面状分布。如果画出来是横七竖八的线段说明道头坐标可能没写或者写的是网格号。4.3 手写numpy版复杂场景的兜底方案遇到segyio不支持的字段布局时用numpy结构化数组自己解。思路是先读二进制头的采样点数和格式码算出一道总字节数再批量解析道头import numpy as np def read_trace_coords(path, trace_bytes, trace_count): dtype np.dtype([ (seq, i4, (1,)), # 1-4 (trace, i4, (1,)), # 5-8 (cdp, i4, (1,)), # 9-12 (code, i4, (1,)), # 17-20 (scalar, i2, (1,)), # 71-72 (sx, i4, (1,)), # 115-118 (sy, i4, (1,)), # 119-122 (cdp_x, i4, (1,)), # 181-184 (cdp_y, i4, (1,)), # 185-188 ]) # 仅读取每个道的道头部分跳过数据体 # 使用mmap避免把整个大文件读进内存 arr np.memmap(path, dtypenp.uint8, moder) trace_start 3600 3200 raw arr[trace_start: trace_start trace_count * trace_bytes] raw raw.reshape(trace_count, trace_bytes) headers raw[:, :240] # 按需复制出字段 cdp_x_bytes headers[:, 180:184].copy() cdp_y_bytes headers[:, 184:188].copy() cdp_x cdp_x_bytes.view(np.int32).ravel() cdp_y cdp_y_bytes.view(np.int32).ravel() return cdp_x, cdp_y注意字节切片要按0基索引字段起始字节减1。这个方案只copy需要的部分不会一次性把数据体全部导入内存适合几个GB甚至上百GB的文件。5. 坐标错误的三类真实翻车现场5.1 现场一坐标全在同一位置或者全是0排查链路先看71-72字节scalar是否为0或异常值。确认你读的是CDP坐标还是炮点坐标。有的数据只写了炮点坐标181-188全为0。检查有没有扩展道头。rev1允许在3600字节文本头之后、正式道数据之前插入扩展文本头某些系统写的扩展文本头会干扰偏移计算。如果以上都正常把240字节道头的前64个字节用十六进制dump出来肉眼对比几道是否一致判断是否整道头都没写。这个翻车现场最常见的结论是数据提供方压根没在道头里写坐标后续需要靠SPS文件或者单独的位置表来关联。5.2 现场二坐标和图件对不上差了一点点或直接翻转坐标偏移几个公里、几十公里常见原因有三个scalar标定漏做、X/Y弄反、投影带搞错。判断方法很实用取数据中两个相邻CDP点算平面距离看是否约等于理论道距。比如理论道距是25米算出来是25000米那基本是标定因子漏了算出来距离对但方位角反了优先怀疑X/Y互换。X/Y互换还有个特征画出来的测线走向和实际工区的长方形边界成90度旋转关系。很多系统里“X北、Y东”和“X东、Y北”两种习惯并存读取后一定要和工区范围描述核对。5.3 现场三三维工区画出来像被撕碎的网格三维数据可视化时出现“锯齿形”和“撕裂感”通常是因为道号顺序和文件存储顺序不一致或者采集时炮点按非规则顺序排列。解决办法不是改绘图代码而是用CDP_X和CDP_Y重新排序或者先用inline/crossline编号排序再画。还有一种情况是dummy道占了位置但道识别码没写对。处理前先按17-20字节过滤只看code1的地震道散点图立刻干净很多。6. 大文件实操怎么从上百GB的SEGY里几秒只读道头6.1 用np.memmap配合步长抽取上百GB的SEGY文件不可能整体读进内存。最实用的是内存映射加只读道头切片。二进制头里的采样点数和格式码先拿到的然后算出每道总字节数用下面这套逻辑import numpy as np def memmap_trace_headers(path, trace_start, trace_bytes, trace_count): arr np.memmap(path, dtypenp.uint8, moder) data arr[trace_start: trace_start trace_count * trace_bytes] data data.reshape(trace_count, trace_bytes) # 取出每道前240字节作为道头 headers data[:, :240] return headers配合5%抽稀做质控图几百GB的数据几秒就能出轮廓。实际使用中我一般先map后只取坐标字段的字节段减少复制量。6.2 两步抽稀法先轮廓后细节直接渲染几百万个散点会让matplotlib卡死。我自己习惯分两步第一次只取前1万道或均匀抽稀5%看整体工区轮廓和坐标范围第二步再全量读取所有道头坐标输出成CSV或GeoJSON交给前端Leaflet或Cesium显示。工区轮廓没问题后再谈逐道渲染。6.3 输出GeoJSON给Web端处理完的坐标最好直接落成GeoJSON后续无论是Web可视化还是GIS叠加都很方便import json features [] for x, y in zip(xs[mask], ys[mask]): features.append({ type: Feature, geometry: { type: Point, coordinates: [float(x), float(y)] }, properties: {} }) geojson {type: FeatureCollection, features: features} with open(cdp_points.geojson, w) as fp: json.dump(geojson, fp)如果后面还要和地形数据、DEM叠合坐标系必须统一否则经纬度对了、高程差几十米的情况也会出现。最后再分享一个我自己的习惯拿到任何SEGY数据我先不画图先做一张质控图——把道头里的道集号按顺序着色看颜色是否连续。跳跃太大说明文件边界或道序有问题这时候画的任何“漂亮图”都是不可信的。道头解析最好做成一个脚本入参是文件路径、坐标字段选择、是否做投影转换输出GeoJSON加PNG这样以后每来一份新数据都能在5分钟内完成定位质控。这个习惯帮我躲过很多次“坐标漂移”的坑建议你也留一份。
