Python读取SU数据

SU数据格式

以下简称 Seismic Unix 为 SU 。SU 格式是 SEGY 的简化，没有前面3600字节文件头，故需要从第一道的道头获取需要的信息。这里默认所有道的采样点数即 nt 是一样的。

SEGY 和 SU 格式的道头都有 240 字节，但 SEGY 只有前 180 字节有信息，SU 的 181-240 字节定义了画图相关信息，可通过查看 SU 源码中的 segy.h 头文件获知。

SU 数据道头中主要的信息见下表：

Python代码

from matplotlib.colors import Normalize
import matplotlib.pyplot as plt
import numpy as np
import os
import struct


def read_su(path):
	data_size = os.stat(path).st_size
	byte2ushot = lambda x: struct.unpack('H', x)
	byte2float = lambda x: struct.unpack('f', x)
	with open(path, 'rb') as f:
		f.seek(114, 0)
		nt = byte2ushot(f.read(2))[0]
		dt = byte2ushot(f.read(2))[0]
		f.seek(70, 1)
		dx = byte2float(f.read(4))[0]
		nx = int(data_size / (4 * nt + 240))
		f.seek(0, 0)
		data = np.zeros((nx, nt), dtype='float32')
		for i in range(nx):
			f.seek(240, 1)
			dat = f.read(4 * nt)
			data[i] = np.array(struct.unpack('f' * nt, dat))
	return data, dt, dx


data, dt, dx = read_su('syndata3')
nx, nt = data.shape
dt /= 1e6

# 坐标轴刻度标签，注意根据需要调整
xtick = list(range(0, nx, 10))
xticklabel = [f'{i * dx:.2f}' for i in xtick]
ytick = list(range(0, nt, 10))
yticklabel = [f'{i * dt:.2f}' for i in ytick]
# 画图缩放参数，注意根据需要调整
zoom = 1
norm = Normalize(vmin=np.min(data) * zoom, vmax=np.max(data) * zoom)

fig, ax = plt.subplots()
ax.imshow(data.T, aspect='auto', norm=norm)
ax.xaxis.set_ticks(xtick)
ax.xaxis.set_ticklabels(xticklabel)
ax.yaxis.set_ticks(ytick)
ax.yaxis.set_ticklabels(yticklabel)
ax.set_xlabel('distance (m)')
ax.set_ylabel('time (s)')
plt.show(block=1)

运行结果

 注意

因为我 SU 装在虚拟机里，故不需要大小端转换，如有需要可参考Python读取大端SEGY数据。

画图时根据自己的数据调整参数。

写 SU 数据可参考我之前发的Python读写SEGY数据的简单实现，但注意那篇文章是将信息写在前面3600字节文件头中，写 SU 数据时需将信息写入道头。