CC BY
30
2017
Геология
Том 16, № 3
ГЕОФИЗИКА, ГЕОФИЗИЧЕСКИЕ МЕТОДЫ ПОИСКОВ ПОЛЕЗНЫХ ИСКОПАЕМЫХ
УДК550.344.094+550.34.01
Seismic Refraction Modeling Using Finite Difference Method and its Implications in the Understanding of the First Arrivals
Y.P. Goyesa, S. Khuramab, G. Reinab, G. Jimenezc
aPerm State University. 15 Bukireva Str., Perm 614990, Russia. E-mail: goyes.yesid@gmail.com
bFaculty of Geology, Santander Industrial University, Cra 27 Calle 9 Ciudad Universitaria, Bucaramanga, Santander680002 Colombia E-mail: skhurama@uis.edu.co
cResearch Group in Basic and Applied Geology (GIGBA), Santander Industrial University, Cra 27 Calle 9 Ciudad Universitaria, Bucaramanga, Santander 680002, Colombia. E-mail: gjimenezd@gmail.com
(Статья поступила в редакцию 2 сентября 2016 г.)
This paper presents the results of the seismic refraction modeling using finite differences method (FDM), implemented in the program from the open source package Seismic Unix (SUFDMOD2), and its comparison with the modeling by time-term method (TTM) from the SEISIMAGER© software. The applied velocity model corresponds to the typical measurements situation in a seismic refraction survey. The depth of refraction interfaces varies approximately from 2 to 6 meters, thus allowing examining the propagation of head waves with the variable dip angles. The synthetic seis-mograms allow us identifying the first arrivals of head waves, which are subsequently the travel-time curves. The analysis of results obtained with FDM and TTM algorithms shows a high correlation with the refraction arrival times, but a low correlation with the arrival times of the direct wave. Finally, the obtained results allow concluding that the seismic modeling of the propagation of head waves using the method of finite differences makes it possible to accurately determine the first arrival time within the complex geological conditions and the velocity dispersion between the layers. The goal of this work is to show how these findings can be applied to seismic modeling by the method of finite differences for calculation of the first arrival time and clarify the results of measurements conducted in real conditions.
Keywords: seismic modeling, seismic refraction, time-term inversion, finite differences, travel time curve, first arrival time, Seismic Unix. DOI: 10.17072/psu.geol.16.3.256
Introduction
As described in Dobrin (1976), the seismic modeling has been an important tool for
analyzing various problems mainly in the exploration of hydrocarbons. Their algorithms have been improved in terms of the numerical approximation and propagation in
© Goyes Y.P., Khurama S., Reina G., Jimenez G., 2017
256
the elastic mediums (Sayers, 2009). According to Bridle (2016), and due to the inclusion of geophysical methods in the determination of exploration geotechnical data, the seismic modeling was taken to study the behavior of waves on different geometries and acoustic contrast between layers (Budhu, 2010). For example, such as the contact between soil and residual colluvial deposits, where the acoustic impedance is low, and the contact between the colluvial deposits and bedrock, where the impedances high (Redpath, 1973).
In this work, we study the propagation of head waves in an elastic medium over a contact between unconsolidated sediments and bedrock, which enables to see the changes in the first arrivals of the P wave on the surface, showing significant changes in the traveltime curves in the areas with tilt refractors (Lankston, 1988). Traditionally, the modeling of the seismic refraction focuses in calculating the first arrival time (FAT), which is performed using mainly the time-term method (TTM). An example of this method is shown in (Iwasaki, 2002) and (Berry et al, 1966). In this work, to study the first arrivals picking, the method of finite differences (FDM) was used, because it allows one to analyze the times (FAT) and the record of the wave propagation for each geophone. Finally, a comparison using linear fitting, coefficient of determination, and residual mean square between the results obtained by FDM and TTM modeling was made.
Methodology
The modeling by time-term method (TTM) was conducted using the module of the Plotrefra (SEISIMAGER), which includes the capability of creating a custom velocity model for forward modeling purposes. The first arrival time was calculated by time-term method, according to (Bath, 1978). This technique is a linear least-squares approach to determine the best discrete-layer solution to the data. The math behind this technique is comparatively simple. A two layers single receiver-source model is showed in Fig. 1 with the next configuration;
Vi and V2 are the velocities for the layer 1 and 2 respectively; ic is the critical angle; Zi and Z2 are the perpendicular depth measure from source and receiver position respectively. In this case, the time (t) is calculated as shown in the following equation:
t = czx + cz2 + xS2
(1)
where Sr = 1/v1? S2 = 1/v2 ,and c = 2Si cos(z'c), Si, S2 -slowness in layer 1 and 2 respectively, c - substitution constrain slowness-dependent, x -distance from source to receiver.
Fig. 1. Two layers scheme to explain the calculation oof the trajectory and time by TTM method. Modified from (Iwasaki, 2002)
Generalizing the equation (1) for a model with multiple layers,
t =llCJkZk + XjS 2
(2)
k=1
we get in a matrix form:
v C ,
V m1
'12
'22
'32
• c
1n
• c
2n
X,
X,
\f -7 \
•c
3n
42
••
•c
X
' m2
•c
S,
m JV 2 J
ft \
, t , V m J
where number of travel times j = i..m, and number of receivers k = i..n (depths to be calculated). In this case, it is possible to solve
c
1
t
c
z
21
2
2
c
X
z
3
3
3
t
c
z
n
X
the matrix for Z1...Zn and S2.
The first arrival time calculation using the finite difference method (FDM) was performed using the free software package SeismicUnix (program SUFDMOD2), which works on the Linux operating system (in this work, Fedora 20). This program uses finite-difference method (Stockwell, 1999) to solve the 2D acoustic wave equation:
d 2((x, z, t)
dt2
= v 2( x, z )V (( x, z, t ), (3)
where v(x,z) is the velocity in the acoustic media.
The Laplacian operator can be approximated with central difference operators. Where the velocity, spatial sampling rate, and grid spacing are in consistent units.
V 2(
j - 2( + j j - 2( + j
(4)
Ax2 Az2
The time derivative is calculated by a second order finite difference scheme:
^ 2((t )
dt2
((t + At) - 2((t) + ((t - At) At2
(5)
Vmax At s 1
Ax V 2
(6)
Vmax - maximum velocity in the model (m/s);
At - time interval in seconds;
Ax - spatial interval in the x-axis (meters).
Table 1. Grid spacing parameters
Model nx Ax nz Az
M1 24 2m 100 0.14m
M2 48 1m 100 0.14m
M3 500 0.0922m 500 0.0281m
Finally, each finite difference scheme has a stability condition (Slawinski et al, 1999). The stability condition for the second order scheme used by SUFMOD2 is shown in the equation 6. This program does not require the grid spacing to be equal in the horizontal and vertical directions. Figure 2 describes the general flow diagram with the list of processes required in SUFDMOD2 for the generation of synthetic seismic records. To create models with complex geometry, we recommend using the package TRIMODEL, which calculates a triangulated velocity model of subsurface. Subsequently, this model is made uniform by using the program TRI2UNI. In this work, three variants of the same model, which differ only in the grid spacing, were studied. In the Table 1 Az and Ax are the spatial interval in the z-axis and x-axis respectively, and nx and nz are the sample numbers in the z-axis and x-axis respectively.
Fig. 2. A general processing flow diagram of using SUFDMOD2program in SeismicUnix
We used for both modeling methods the general model shown in Fig. 3 with the area of thickness of 2 m simulating a geometry of the completely flat refractor. The zone B (1250 m/s) has a variable thickness 0 - 4 m that corresponds to a refractor with complex geometry simulating the advanced processes of weathering on the bedrock indicated in the general model by the area C. We used a sequence of five shot-points with different location. Table 2 shows the position of each shot.
Table 2. Shot points location
Shot point Location
SP1 1 m
SP2 12 m
SP3 23 m
SP4 34 m
SP5 45 m
-2 -4
■B -6 0 -8 -10 -12
0 10 20 30 40 46
Distance (m)
Fig. 3. General model. (A) Colluvial zone Vpi=300m/s; (B) residual soil zone Vp2=i250m/s with flat refractor at z=-6m and complex refractor at x= (20, 42) m; (C) bedrock zone Vp3=2500m/s
Distance (m)
Fig. 4. Uniformly sampled velocity model of layered medium (M3) in binary format used for SUFDMOD2 program
The uniform models (Fig. 4) were built using linear interpolation of the data obtained by seismic tomography. SeismicUnix creates different models using the program UNI2, which generates a 2-D uniformly sampled velocity profile from a layered model. In each layer, velocity is a linear function of position.
Results
The modeling with finite differences for the velocity model with 3 different grid spacing (Table 1) produced a variation in the optimization of stable At (6) and the number of iterations (Table 3).
A seismogram with the traces corresponding to the registration of the P wave to geophone spread was obtained for each model settings, which differ in values nx (number of samples in x-axis) and nz (number of depth-samples).
It was found that the variation in the number of geophones used for seismic acqui-
Table 3. Computing parameters during the model generation by the SUFDMOD2 program for different grid spacing
Model No. ofitera- CPU Stable
tions time At(s)
M1 3450 1.101 s 2.898e-05
M2 3751 1.745 s 2.665e-05
M3 18761 436.17s 5.338e-06
sition has a significant impact on the quality of registration of seismic events. In the seismogram with nx=24 (Fig. 5), the first break of refracted wave is reliably identified, but the accurate arrival time of the direct wave close to position of the source is difficult to define. Using nx=48 (Fig.6), we are able to improve the sampling of the waves between 0.04 and 0.05 seconds, and only seismogram obtained with nx=500 (Fig. 7) shows clearly the head waves 1 and 2, and direct wave.
The results obtained using TTM and FDM methods are illustrated in Fig. 8. The time of arrival of the P wave on the same model using nx = 24 was calculated with each method. In Fig. 8A, the comparison between the travel-time curves shows a high dispersion of the data for time less than 0.015 s. For times exceeding 0.015 s (first refraction time in Fig. 8C), the evidenced changes in the slopes of the curve corresponding to the different layers of the subsoil are clearly recognized. The highest difference (FD -TTM) about 0.008 s (Fig. 8B) corresponds to the direct wave. To analyze the dispersion of the results, we performed a linear fit between them shown in Fig. 8C.
This curve was divided into two zones (high correlation zone and low correlation zone). The value of calculated coefficient of determination was determined as 0.956866 and the straight fit equation was defined as timeRT=1.07871478*timeFD - 0.00284832, where timeRT and timeFD are the times calculated from TTM and FDM data respectively.
Analysis of the seismic survey data
A seismic refraction survey in Barranca-bermeja (Colombia) was developed with the
aim to identify changes in the subsurface and define the depths and types of soils.
The seismic acquisition was performed using a seismograph ES-3000 with 24 vertical 14 Hz geophones separated 2 meters. The total length was 46 meters. The result of the seismic processing is shown in Fig. 9. The results of seismic tomography processing correlate with the data of the geotechnical
survey.
Fig. 7. Shot gather with 500 geophones, point x=1m
shot
Fig. 5. Shot gather with 24 geophones, shot point x=1m
Fig. 6. Shot gather with 48 geophones, shot point x=1m
For optimization the temporal parameters in future acquisition surveys over the same area, we performed a seismic modeling using FDM. Additionally, the modeling allowed us to know the validity of the obtained results by comparing the synthesized traces with the survey traces (Fig. 10).
The events associated with the refraction waves were found between 0 and 0.05 seconds. The seismic survey should decrease the maximum time of the acquisition but increasing, on other hand, the detail of the head waves to improve the process of first arrivals picking.
The arrival times obtained in the seismic record x=0m using FDM and real data, procurement of field (RD) are shown in Fig. 11.
The slope preserves the general trend, but the data differ in terms of lap times. In blue curve shows the difference between two times. The maximum difference is 0.011 seconds.
Conclusions
Modeling using the finite difference method showed precision in the calculation of the first arrival time of P-wave, which was verified using statistical correlations among the times calculated by the TTM method and the times obtained from the seismogram (first arrival time picking using the FDM method).
Fig. 8. First arrival times obtained using FDM and TTM. (A) travel time curves, (B) time difference between methods, (C) first time arrivals scattering between methods and linear fitting to compare the correlation
Red Clay with sand
Brown Clay with veins
Yellow Clay high cohesion
0 5 10 15 20 25 30 35 40 45 "" Distance (m)
Fig. 9. Seismic tomography in Barrancabermeja refraction survey
Fig. 11. Travel time for FDM and field , in x=0m source.
set
Fig. 10.Shot gathers in x=0m. (a) FDM modeling and (b) survey field data set with first time arrival (red line)
The modeling allowed an analysis and study of seismic survey on models with different geometrical configurations and variation of the seismic velocity with depth. With respect to the number of geophones necessary to record the seismic events involved in seismic refractions, 24 and/or 48 geophones
are sufficient, and this is a configuration similar to the standard equipment used in real seismic survey.
This type of seismic modeling allows studying the arrival times and synthetic seis-mograms in layered and uniform subsurface models. The application varies from verification of data inversion algorithms to the optimal seismic survey geometry determination for various structure of the substrate in complex areas.
Acknowledgement
The authors express their acknowledgment to the company of geotechnical exploration and engineering (TORRES ING S. A. S.) and Associate Professor Oleg Kovin for their recommendations to improve this work.
References
Bath M. 1978. An analysis of the time term method in refraction seismology. Tecto-nophysics. 51(3—4): 155-169. doi:10. 1016/0040- 1951(78)90238-X Berry M.J., West G.F. 1966. An interpretation of the first arrival data of the Lake Superior experiment by the time-term method, Bull. Seism. Soc. Am., 56:141171.
Bridle R. 2016.Introduction to this special section: Near-surface modeling and imaging. The Leading Edge. 35(11):938-939. doi: 10.1190/tle35110938.1
Budhu M. 2010. Soil mechanics and foundation. John Wiley & Sons. p. 781. Dobrin M.B., 1976. Introduction to Geophysical Prospecting, p. 629. Iwasaki T. 2002. Extended time-term method for identifying lateral structural variations from seismic refraction data. Earth Planets Space. 54:663-677. doi: 10.1186/BF03351718 Lankston R. W.1988. High Resolution Refraction Data Acquisition and Interpretation. Symposium on the Application of Geophysics to Engineering and Environmental Problems. pp. 349-408.doi: 10.4133/1.2921806 Redpath B. 1973. Seismic refraction exploration for engineering site investigations. Technical Report E-73-4 U.S. Army Corps of Engineers Explosive Excavation Research Laboratory, Livermore, California. p. 55. Sayers C., Chopra S. 2009. Introduction to special section: Seismic modeling. The Leading Edge, 28:528-529. doi: 10.1190/1.3124926 Slawinski L.R., Bording R.P., 1999. A recipe for stability of finite-difference wave-equation computations. Geophysics. 64:967-969. Stockwell Jr.J.W. 1999. The CWP/SU: Seismic Un*x Package, Computers and Geo-sciences, 25(4): 415 - 419. doi: 10.1016/S0098-3004(98)00145-9
Моделирование преломленных волн с использованием метода конечных разностей и его применение для исследования времён первовступлений
Е.П. Гойес a, С.Курама^ Г. Рейна^ Г. Химене^
^Пермский государственный национальный исследовательский университет, 614990, Пермь, ул. Букирева, 15. E-mail: goyes.yesid@gmail.com
^антандерский индустриальный университет, 680002, Букараманга, Колумбия. E-mail: sait.khurama@uis.edu.co
Исследовательская группа в области фундаментальной и прикладной геологии, Сантандерский индустриальный университет, 680002, Букараманга, Колумбия. E-mail: gjimenezd@gmail.com
Представлены результаты моделирования преломленных волн методом конечных разностей (МКР), полученные с помощью программы с открытым кодом SUFDMOD2 (пакет Seismic Unix), и сравнения их с моделированием по методу годографа (ТТМ), проведенным с использованием программ пакета SEISIMAGER. Использовалась скоростная модель, соответствующая типичным измерениям методом преломленных волн (МПВ). Преломляющая граница расположена наклонно с минимальной и максимальной глубинами 2 и 6 м соответственно, что позволяет изучать распространения головных волн на границах с переменным углом падения. Синтетические сейсмограммы позволили выявить первые вступления головных волн, которые изучаются годографами. Анализ результатов, полученных с помощью МКР и ТТМ, показывают высокую корреляцию на участке годографа прихода преломленных волн и низкую корреляцию на участке годографа прямой волны. Наконец, полученные результаты позволяют сделать вывод, что сейсмическое моделирование распространения головных волн с использованием метода конечных разностей даёт возможность четко определить время первого вступления сложной геометрии и разброс скоростей между слоями. Цель этой работы - показать, как можно применить сейсмическое моделирование по методу конечных разностей для расчета времени первого вступления и уточнения результатов проведения измерений в реальных условиях.
Ключевые слова: сейсмомоделирование, МПВ, годограф, времена первовступле-ний, Seismic Unix.
