Как открыть SEGY-файл с помощью Python

Формат SEG-Y – один из стандартных форматов для хранения геофизических (как правило сейсморазведочных) данных, разработанный Society of Exploration Geophysicists (SEG) в 1973 году. Данные в формате сегвай записаны в бинарном виде, и поэтому недоступны для их просмотра без специального программного обеспечения.

Один из альтернативных способов просмотра и обработки данных в формате seg-y – считать их с помощью Python с предварительно установленным на него геофизическим пакетом пакетом obspy.

Используя стандартные функции, имеющиеся в пакете obspy можно считывать и визуализировать данные сейсморазведочных SEGY-файлов.

Для открытия sgy-файла нужно прописать следующий код:

from obspy.io.segy.core import _read_segy
from obspy.core.util import get_example_file
import numpy as np

filename = get_example_file("SEGY-file.sgy")

segyfile= _read_segy(filename, unpack_trace_headers=True)

С помощью print выводим количество трасс, содержащихся в окрытом файле .segy, а также считываем заголовок первой трассы и записываем его параметры в список:

attr= list(segyfile[0].stats.segy.trace_header)

vals=[]
for obj in segyfile[0].stats.segy.trace_header:
    vals.append(segyfile[0].stats.segy.trace_header[obj])

for i in range (0,len(attr)):
    print (attr[i], vals[i])

В результате исполнения данного кода будут выведен список всех параметров, записанных в шапке первой сейсмической трассы. В данном примере, в шапке трассы содержались следующие данные:

endian <
unpacked_header None
trace_sequence_number_within_line 1
trace_sequence_number_within_segy_file 272
original_field_record_number 272
trace_number_within_the_original_field_record 1
energy_source_point_number 0
ensemble_number 558
trace_number_within_the_ensemble 1
trace_identification_code 1
number_of_vertically_summed_traces_yielding_this_trace 0
number_of_horizontally_stacked_traces_yielding_this_trace 0
data_use 1
distance_from_center_of_the_source_point_to_the_center_of_the_receiver_group 0
receiver_group_elevation 0
surface_elevation_at_source 0
source_depth_below_surface 0
datum_elevation_at_receiver_group 0
datum_elevation_at_source 0
water_depth_at_source 0
water_depth_at_group 0
scalar_to_be_applied_to_all_elevations_and_depths 1
scalar_to_be_applied_to_all_coordinates -100
source_coordinate_x 94920
source_coordinate_y 73580
group_coordinate_x 0
group_coordinate_y 0
coordinate_units 1
weathering_velocity 0
subweathering_velocity 0
uphole_time_at_source_in_ms 0
uphole_time_at_group_in_ms 0
source_static_correction_in_ms 0
group_static_correction_in_ms 0
total_static_applied_in_ms 0
lag_time_A 0
lag_time_B 0
delay_recording_time 0
mute_time_start_time_in_ms 0
mute_time_end_time_in_ms 0
number_of_samples_in_this_trace 1001
sample_interval_in_ms_for_this_trace 4576
gain_type_of_field_instruments 0
instrument_gain_constant 0
instrument_early_or_initial_gain 0
correlated 1
sweep_frequency_at_start 0
sweep_frequency_at_end 0
sweep_length_in_ms 0
sweep_type 1
sweep_trace_taper_length_at_start_in_ms 0
sweep_trace_taper_length_at_end_in_ms 0
taper_type 1
alias_filter_frequency 0
alias_filter_slope 0
notch_filter_frequency 0
notch_filter_slope 0
low_cut_frequency 0
high_cut_frequency 0
low_cut_slope 0
high_cut_slope 0
year_data_recorded 2019
day_of_year 340
hour_of_day 11
minute_of_hour 52
second_of_minute 10
time_basis_code 1
trace_weighting_factor 0
geophone_group_number_of_roll_switch_position_one 0
geophone_group_number_of_trace_number_one 0
geophone_group_number_of_last_trace 0
gap_size 0
over_travel_associated_with_taper 0
x_coordinate_of_ensemble_position_of_this_trace 94920
y_coordinate_of_ensemble_position_of_this_trace 73580
for_3d_poststack_data_this_field_is_for_in_line_number 272
for_3d_poststack_data_this_field_is_for_cross_line_number 558
shotpoint_number 0
scalar_to_be_applied_to_the_shotpoint_number 0
trace_value_measurement_unit 0
transduction_constant_mantissa 0
transduction_constant_exponent 0
transduction_units 0
device_trace_identifier 0
scalar_to_be_applied_to_times 0
source_type_orientation 0
source_energy_direction_mantissa 0
source_energy_direction_exponent 0
source_measurement_mantissa 0
source_measurement_exponent 0
source_measurement_unit 0
unassigned b'\x00\x00\x00\x00\x00\x00\x00\x00'

При необходимости эти данные можно использовать в дальнейшем для последующей обработки сегвай-файлов, например, выбора необходимой трассы по номеру или ее координатам и т.д.

Для того, чтобы создать массивы данных, содержащие координаты (X и Y) и сейсмические сигналы для каждой сейсмотрассы можно воспользоваться следующим кодом:

x_coord=[]
y_coord=[]
Traces=[]

for tr in range (0, len(segyfile)):
    x_coord.append(segyfile[tr].stats.segy.trace_header.source_coordinate_x)
    y_coord.append(segyfile[tr].stats.segy.trace_header.source_coordinate_y)
    Traces.append (segyfile[tr].data)

Если необходимо визуализировать сигнал для определенной сейсмотрассы можно воспользоваться данным кодом:

import numpy as np
import matplotlib.pyplot as plt

dt=segyfile[0].stats.segy.trace_header.sample_interval_in_ms_for_this_trace

ntraces = len(segyfile.traces)
nsamples = len(segyfile.traces[0].data)

times= np.linspace(0, dt*nsamples, nsamples)

fig= plt.figure(1, figsize=(12, 4))
plt.plot(times, Traces[5000], label='Amplitude')

plt.grid(True)
plt.legend(loc='upper right')
plt.show()

В результате мы получим изображение сейсмической записи (в данном случае по оси Y выводится амплитуда сигнала, а по оси X время в миллисекундах) для выбранной нами сейсмотрассы под порядковым номером 5000:

Используя и другие встроенные функции пакета obspy можно производить различные операции над сейсмическими данными, например производить широкополосную фильтрацию сигналов по частоте, строить спектры сигналов и т.д.

Рекомендуем также прочитать нашу статью о том, как открывать las-файлы с помощью библиотеки lasio.

Понравилась статья? Сделай репост: