GeophyAI

Python与地球物理数据处理

0%

SEGY格式数据读取之SEGYIO

发表于 分类于 Valine：
简要介绍了如何在Python中使用SEGYIO读取SEGY格式地震数据。

SEGYIO包使用简介

本文简要介绍python包SEGYIO(安装方法：pip install segyio==1.5.1)的使用。
SEGIOGit链接：https://github.com/equinor/segyio

在使用马达加斯加处理地震数据时，遇到了需要读取SEGY格式数据的情况(book/rsf/su/rsflab16)。为了更好地利用python对国际通用格式的地震数据进行处理，结合book/rsf/su/rsflab16中的需求，我们简要介绍如何在python中使用SEGYIO包快速读取此类格式的数据。

我们利用SEGYIO主要解决马达加斯加源码示例book/rsf/su/rsflab16中的以下需求：

  1. 读取SEGY格式文件中的剖面数据；
  2. 读取SEGY格式文件中偏移距和CDP号数据；
  3. 根据CDP道集号从$1$中抽取数据

注：在book/rsf/su/rsflab16文件夹下执行scons即可下载本次实验用到的SEGY文件seismic.segy。

读取SEGY道集信息

1
2
3
4
5
6
7
8
import segyio
root_path = r'your_file_path'
with segyio.open(os.path.join(root_path, 'seismic.segy')) as segyfile:
    # Memory map file for faster reading (especially if file is big...)
    segyfile.mmap()
 	index = 0
    # 打印第index道的道头信息
    print(segyfile.header[index])

输出(只列出了部分信息)：

1
2
3
4
5
6
7
{TRACE_SEQUENCE_LINE: 1, TRACE_SEQUENCE_FILE: 1, FieldRecord: 3, TraceNumber: 1, 
EnergySourcePoint: 101, CDP: 1, CDP_TRACE: 1, offset: -3237, 
ReceiverGroupElevation: -10, SourceSurfaceElevation: -6, SourceDepth: 0, ReceiverDatumElevation: 0,
SourceDatumElevation: 0, SourceWaterDepth: 0, GroupWaterDepth: 0, ElevationScalar: 1, 
SourceGroupScalar: 1, SourceX: 3237, SourceY: 0, GroupX: 0,
GroupY: 0, CoordinateUnits: 3, MuteTimeEND: 48, TRACE_SAMPLE_COUNT: 1500, 
TRACE_SAMPLE_INTERVAL: 4000

读取剖面数据

如果只读取其中的数据，则可以利用以下方法：

1
2
3
data = segyio.tools.cube(os.path.join(root_path, 'seismic.segy'))
data = np.squeeze(data)
# data.shape=(120120, 1500)

输出data的shape时可以发现，这里将所有道的数据都读取到了一个数组中，知道了每炮对应的道数后，我们可以手动对这些数据按炮进行划分：
1
2
3
4
5
DATA = []
for i in tqdm.trange(120120//120):
	DATA.append(data[i*120:(i+1)*120,:])
DATA = np.array(DATA)
# DATA.shape = (1001, 1500, 50)

偏移距和CDP号读取

偏移距和CDP号存储在每一道的道头中，我们可以使用如下方法获得二者信息：

1
2
3
4
5
with segyio.open(os.path.join(root_path, 'seismic.segy')) as segyfile:
    # Memory map file for faster reading (especially if file is big...)
    segyfile.mmap()
    CDP = segyfile.attributes(segyio.TraceField.CDP)[:]
    OFFSET = segyfile.attributes(segyio.TraceField.offset)[:]

以上代码中CDP和OFFSET中存储的偏移距和CDP号数据与SConstruct源码中的如下代码所生成的cdp.rsf和offset.rsf文件中数据是相对应且一致的。
1
2
3
4
# SConstruct 源码读取偏移距和CDP号
for key in ('cdp','offset'):
    Flow(key,'tseismic',
         'dd type=float | headermath output=%s' % key)

根据CDP道集号从$1$中抽取数据

首先来看SConstruct源码：

1
2
3
4
5
6
# Capture a single CMP
Flow('cdpmask','cdp','mask min=265 max=265')
Flow('cdp265','seismic cdpmask',
     'headerwindow mask=${SOURCES[1]} | pow pow1=2')
Flow('offset265','offset cdpmask',
     'headerwindow mask=${SOURCES[1]}')

以上源码中，第一个Flow生成了一个与cdp中数据等长的全零元素数组，令与cdp中数字265对应的索引为1，其余位置为0。另外，第三个Flow还利用cdpmask.rsf文件从偏移距数据中抽取了CDP号为265的位置对应的数值。例如：
1
2
3
4
5
6
# 若
cdp = [2, 33, 256, 265, 8, 265]
offset = [1, 8, 16, 16, 2, 32]
# 则
cdpmask = [0, 0, 0, 1, 0, 0]
offset265= [16, 32]

然后从原始数据中抽取了CDP号为265的道形成了新的数据，在利用sfpow对新数据进行增益处理后生成了cdp265.rsf文件。

这相当于利用numpy数组进行如下操作：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
cdp_data = data[np.where(cdp==265)[0]]
def power(data, p = 2, o = 0., d = 1., axis = -1):
	# p：增益指数项
	# o：起始位置
	# d：时间采样间隔
	# axis：沿哪个轴进行增益
    nx = data.shape[axis]
    reps = data.shape[0]
    scale = np.zeros((1, nx))
    x_range = np.arange(0, nx) + 1
    scale = np.power(o+x_range*d, p)
    return data*scale, scale
cdp_data_power, scale = power(cdp_data, d=0.004)
offset265 = OFFSET[np.where(cdp==265)[0]]

代码中np.where(cdp==265)输出了数组cdp中元素值为265的坐标索引，例如：cdp=[2, 33, 256, 265, 8, 265],np.where(cdp==265)=[3, 5]，由于原cdp数据为二维数组，因此这里我们只输出其第一个维度索引值。上述代码中numpy数组cdp_data_power中的数据与cdp265.rsf文件的数据是一致的，并且offset265与文件offset265.rsf文件中数据也是一致的。