企业官网建设流程全解析

上个月帮朋友处理一批三维采集数据，文件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文件，第一件事不是急着读数据体，而是问自己三个问题：

1. 我的数据是二维测线还是三维工区？二维和三维的坐标字段解读方式完全不同。
2. 道头里存的是网格坐标（inline/crossline）还是真实地理坐标（米或经纬度）？
3. 坐标标定因子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。下面是我实际项目里最常用的一组字段。

这里有个非常容易踩的点：SEGY标准并没有规定X必须是东向、Y必须是北向。国内很多工区资料习惯把X当北向、Y当东向，而国外软件默认X=Easting、Y=Northing。坐标反了的案例我见过不下十次，尤其是从国外软件导出的成果直接套国内坐标系的时候。

2.3 道头里三个容易被忽略的“隐藏信息”

第一是17-20字节的道识别码。可视化前必须按它过滤，哑道和空道混进坐标点会造成散点图出现大量噪声。第二是71-72字节的坐标标定因子。第三是二进制头格式码。这三个字段不确认，后面任何坐标计算都没意义。

3. 坐标定位的真正麻烦：SP/网格/投影三道关

3.1 你手里到底有几个坐标来源

二维数据通常很简单，181-184字节存CDP坐标X，185-188存CDP坐标Y，115-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，说明坐标字段本身就是真实值，不需要换算。

举个例子：存储值是18495008，scalar = -100，真实坐标就是184950.08米；如果漏做除法，直接把这个数当米用，点位会偏到一百多公里外，直接飞出工区。代码里最简单的处理方式：

if scalar > 0: coord_real = stored_value * scalar elif scalar < 0: coord_real = stored_value / abs(scalar) else: coord_real = stored_value

3.3 网格坐标转经纬度，不能硬转

如果道头里存的是高斯平面坐标或UTM坐标，想叠加到底图上必须先做投影转换。用pyproj可以一行完成：

from pyproj import Transformer # 示例：UTM 50N 转 WGS84 经纬度 transformer = Transformer.from_crs("EPSG:32650", "EPSG:4326", always_xy=True) 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_geometry=True) as f: trcount = f.tracecount sample_count = f.bin[segyio.BinField.Samples] scalar = f.header[0][segyio.TraceField.SourceGroupScalar] xs = np.empty(trcount, dtype=np.float64) ys = np.empty(trcount, dtype=np.float64) codes = np.empty(trcount, dtype=np.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], s=1, marker=".", alpha=0.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, dtype=np.uint8, mode="r") 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

排查链路：

1. 先看71-72字节scalar是否为0或异常值。
2. 确认你读的是CDP坐标还是炮点坐标。有的数据只写了炮点坐标，181-188全为0。
3. 检查有没有扩展道头。rev1允许在3600字节文本头之后、正式道数据之前插入扩展文本头，某些系统写的扩展文本头会干扰偏移计算。
4. 如果以上都正常，把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字节过滤，只看code=1的地震道，散点图立刻干净很多。

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, dtype=np.uint8, mode="r") 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分钟内完成定位质控。这个习惯帮我躲过很多次“坐标漂移”的坑，建议你也留一份。