地震学习笔记 | Python读写segy数据 | segyio示例

如果你读过Anaconda + VScode | Python环境安装与管理, 相信你已经配置好了python环境，本文介绍利用Python读写segy数据。

segyio是一个可以读写segy文件的python库。其详细说明可以参考官方文档或Github主页。

这里仅介绍一个segyio使用示例。

1. 准备数据

本文以Marmousi模型为例，读取segy文件。

首先，下载marmousi模型示例数据。

2. 安装segyio

首先, 打开Anaconda Prompt, 创建一个环境seismic (也可以在之前创建好的环境/base环境中直接安装) ：

conda create -n seismic python=3.11

激活该环境：

conda activate seismic 

安装segyio：

pip install segyio

安装画图用的matplotlib库:

pip install matplotlib #或conda install matplotlib

3. 读取segy文件

打开vscode，创建.ipynb文件，选择已安装好segyio的环境seismic。

import segyio
import matplotlib.pyplot as plt
import numpy as np

def read_segy(data_dir,shotnum=0):
    with segyio.open(data_dir,'r',ignore_geometry=True) as f:
        sourceX = f.attributes(segyio.TraceField.SourceX)[:]
        trace_num = len(sourceX) #number of all trace
        if shotnum:
            shot_num = shotnum 
        else:
            shot_num = len(set(sourceX)) #shot number 
        len_shot = trace_num//shot_num   #The length of the data in each shot data
        time = f.trace[0].shape[0]
        print('start read segy data')
        data = np.zeros((shot_num,time,len_shot))
        for j in range(0,shot_num):
            data[j,:,:] = np.asarray([np.copy(x) for x in f.trace[j*len_shot:(j+1)*len_shot]]).T
        return data

读入segy文件（前面下载好的marmousi模型）

segyfile = '../elastic-marmousi-model/model/MODEL_P-WAVE_VELOCITY_1.25m.segy'
marmousi = read_segy(segyfile,shotnum=1)
print(marmousi.shape)

画图：

plt.figure(figsize=(15,3))
plt.imshow(marmousi[0])
plt.colorbar()

试一下另一个文件：

segyfile1 = '../elastic-marmousi-model/processed_data/SEGY-Depth/SYNTHETIC.segy'
marmousi_depth = read_segy(segyfile1,shotnum=1)
print(marmousi_depth.shape)
plt.figure(figsize=(10,5))
plt.imshow(marmousi_depth[0],vmin=-0.01,vmax=0.01,aspect=1)

4. 写入segy文件

创建一个简单的一维数组（可以当成速度），并扩展到二维，画图：

velo = np.arange(1000,5000,100)

def one_to_2d(v_1d, nx =1000):
    return np.ones((v_1d.shape[0],nx))*v_1d.reshape(v_1d.shape[0],1)

velo2d = one_to_2d(velo)
plt.figure(figsize=(8,3))
plt.imshow(velo2d,aspect=10)
plt.colorbar()

定义写入segy的函数：

def creat_velocity_segy(vel, dir, dx = 10, dt = 10):
    seis = vel
    spec = segyio.spec()
    spec.samples = list(np.int32(np.arange(vel.shape[0])* dt))
    spec.format = 5
    spec.tracecount =vel.shape[1]
    with segyio.create(dir, spec) as f:
    ## fill the file with data
        for i in range(vel.shape[1]):
            t = i 
            f.trace[t] = seis[:,i]
            f.header[t][segyio.TraceField.FieldRecord] = 0
            f.header[t][segyio.TraceField.ShotPoint] = 0
            f.header[t][segyio.TraceField.TraceNumber] = i+1
            f.header[t][segyio.TraceField.ReceiverGroupElevation] = 10
            f.header[t][segyio.TraceField.GroupX] = i*dx
            f.header[t][segyio.TraceField.SourceX] = i*dx
            f.header[t][segyio.TraceField.TRACE_SAMPLE_INTERVAL] = 10

写入segy文件：

savepath = './savepath/velomodel.segy'
creat_velocity_segy(velo2d,dir=savepath,dx=10,dt=10)

读取刚刚写的segy文件，画图检查是否正确：

modelcheck = read_segy(savepath,shotnum=1)
plt.figure(figsize=(8,3))
plt.imshow(modelcheck[0],aspect=10)
plt.colorbar()

Tips 将炮集数据抽成共检波点道集

paths= 'common_source.segy'
with segyio.open(paths,'r',ignore_geometry=True) as f:
    sourceX = f.attributes(segyio.TraceField.SourceX)[:]
    trace_num = len(sourceX) #number of trace, The sourceX under the same shot is the same character.
    print('trace_num = ',trace_num)

    shot_num = len(set(sourceX)) #shot number 
    print('source_num = ',shot_num)

    len_shot = trace_num//shot_num   #The length of the data in each shot data
    print('len_shot(trace num of per source) = ',len_shot)

    time = f.trace[0].shape[0]
    print('sampling time = ', time)

    print('start read segy data')
    data = np.zeros((shot_num,time,len_shot))
    for j in range(0,shot_num):
        data[j,:,:] = np.asarray([np.copy(x) for x in f.trace[j*len_shot:(j+1)*len_shot]]).T
    print(data.shape)
    
    crp_data = np.zeros((len_shot,time,shot_num))
    for i in range(len_shot):
        for t in range(time):
            crp_data[i, t, :] = data[:, t, i]
    print(crp_data.shape)