International Journal of Rock Mechanics & Mining Sciences 38 (2001) 807–813

High resolution seismic refraction method using surface and borehole

data for site characterization of rocks
K. Hayashi, T. Takahashi*
OYO Corporation, 43 Miyukigaoka, Tsukuba, Ibaraki 305-0841, Japan
Accepted 17 July 2001

Abstract

The authors have developed a new method for the analysis of seismic refraction data, which can handle both surface and borehole
data simultaneously. The main features of the method we developed are as follows: (1) The algorithm of analysis is based on the
tomographic reconstruction technique, and refraction data can be analyzed almost automatically after traveltime picking. (2) The
method is applicable to both layered velocity structures and those showing a velocity gradient with depth. (3) In order to apply
the seismic refraction technique to a site that has a more complex structure, the method can use the data obtained from sources
and receivers on the surface, as well as those within boreholes simultaneously.
Numerical and actual field examples have shown that the method can greatly improve the accuracy and the reliability of the
seismic velocity models obtained by the seismic refraction method as compared to the conventional analysis methods. It can also
shorten the data processing time drastically and reduce the survey cost. r 2001 Elsevier Science Ltd. All rights reserved.

1. Introduction have limitations in the application to sites where surface

topography and subsurface geological structure are
In Japan, the seismic refraction method has been complex, because they make some assumptions for
widely used in many civil engineering projects in these surface and subsurface features. For example,
mountainous areas, which include railway and highway these methods analyze refraction data supposing that
tunnels, dam constructions and landslide protections. the earth is a layered medium.
The seismic refraction method provides subsurface To expand the applicability of the seismic refraction
seismic velocity models, which are very important for method to such complicated sites, we have developed a
the design and the construction phase of these projects. new method for analyzing seismic refraction data. The
For example, in tunneling sites, the seismic velocity method uses tomographic data analysis for enabling us
model has been used for the excavation plan and the to handle lateral velocity changes in a layer. The method
design of supporting steel intervals. can not only use surface data, but also borehole data for
In the typical seismic refraction method in Japan, imaging vertical structures and thin high velocity layers.
geophones are uses for receivers and dynamite is used In this paper, we describe the details of the new
for the sources with the receiver spacing of 5 or 10 m and analysis method and demonstrate its performance
the source spacing of 30–100 m. Fig. 1 shows the through some numerical examples. In addition, we
example of measured refraction profile and traveltime show two actual examples of the application of the new
curves. method to civil engineering investigations.
There are various methods for analyzing seismic
refraction data, for example, Hagiwara’s method [1],
the Plus–minus method [2], and the generalized recipro- 2. Analysis method
cal method [3]. However, these conventional methods
The method proposed uses non-linear traveltime
*Corresponding author. Tel.: +81-298-51-6621; fax: +81-298-51- tomography consisting of ray tracing for forward
5450. modeling and simultaneous iterative reconstruction
E-mail address: [email protected] (T. Takahashi). technique (SIRT) for inversion, which are commonly

1365-1609/01/$ - see front matter r 2001 Elsevier Science Ltd. All rights reserved.
PII: S 1 3 6 5 - 1 6 0 9 ( 0 1 ) 0 0 0 4 5 - 4
808 K. Hayashi, T. Takahashi / International Journal of Rock Mechanics & Mining Sciences 38 (2001) 807–813

Fig. 3. Principle of the raytracing. Some nodes are arranged on the cell
boundary. A ray is expressed as line connecting these nodes. A
traveltime between a source and a receiver is defined as the fastest
traveltime of all ray paths.

width of each cell is chosen as the receiver interval.

Every cell has an identical velocity value. The initial
velocity model is constructed so that the velocity is
simply increasing with depth.
2. First arrival traveltimes and ray paths are calculated
by the ray tracing method based on Huygens’
principle [4,5]. The ray tracing is done with the
following steps:some nodes are arranged on the cell
boundary, and a ray is expressed as line connecting
these nodes; a traveltime between a source and a
receiver is defined as the fastest traveltime of all ray
paths (Fig. 3).
3. A model is updated by SIRT. The correction value of
the slowness (DSj ) for the jth cell is given by the
following equation:
X 
Fig. 1. (a,b) The example of measured refraction profile traveltime
0c c
 X
curves. DSj ¼ T i lij Sj =Ti lij ;
i i

where Sj is the slowness of the jth cell, lij the length of

the ith ray passing through the jth cell, Ti0c the
traveltime residual of the ith ray, Tic the calculated
traveltime of the ith ray. During the iterations, the
velocity of each cell is updated under the constraint
that velocity increases with depth. In this algorithm,
we can use traveltimes of seismic waves generated or
received within boreholes together with those for the
surface sources and receivers. If the sources or
receivers within boreholes are included, we do not
have to use the constraint that velocity is increasing
with depth.

In civil engineering applications, the subsurface

velocity model is often described by a layered model,
Fig. 2. Processing flow of the new analysis method for seismic in which each layer has a constant velocity. Because, it is
refraction data.
useful for civil engineering design works to classify the
underground area of interest as simply as possible, at
used in cross-hole traveltime tomography. Fig. 2 shows least in the beginning stage of the site investigation.
the processing flow used in the method. Main features of There are also often some sites which have clear
the algorithm are as follows: boundaries of rock quality. The subsurface velocity
model for these sites should be represented by a layered
1. The velocity model is represented by quadrangle cells. one. Since the velocity tomogram obtained by the
The model is divided into many thin layers, and the method described above represents the velocity model
 K. Hayashi, T. Takahashi / International Journal of Rock Mechanics & Mining Sciences 38 (2001) 807–813 809

as a smoothed one, it is not adequate to model rock sites

that have sharp velocity changes. In that case, we apply
the following procedures to extract clear layer bound-
aries and true velocities.
1. Approximate velocities and the number of layers are
decided from the traveltime curves.
2. Some of the velocity contours are extracted as initial
boundaries of layers.
3. This layered velocity model is represented by quad-
rangle cells. The width of each cell is chosen as the
receiver interval, and the height is the thickness of Fig. 4. Velocity model used in the first numerical example.
each layer at each receiving point. A little lateral
velocity variation can be handled in the layered
velocity model.
4. First arrival traveltimes and ray paths are calculated
by ray tracing.
5. The model is updated by SIRT. Here, the correction
is not made only for the slowness of the layer, but
also for the thickness of the layer at each receiving
point [6]. Correction of the thickness for each layer is
done by modifying the height of each cell. The
correction value of the height (DHj ) for the jth cell is
given by the following equation:
P  0c 
Hj0 i Ti lij Sj =Tic Fig. 5. Initial velocity model. The initial velocity model is constructed
DHj ¼ P ;
i lij Sj so that the velocity is simply increasing with depth.

where Hj0 is the height of the jth cell, Sj the slowness

of the jth cell, lij the length of the ith ray passing
through the jth cell, Ti0c the traveltime residual of can see that there is large difference between observed
the ith ray, Tic the calculated traveltime of the ith ray. and calculated traveltimes for the initial velocity model.
Calculation of the slowness correction value for each The initial model is updated iteratively by applying ray
layer is done with the following procedures. Firstly, tracing and SIRT. Fig. 7 shows the result of the
the traveltimes of the refraction waves from the upper tomographic reconstruction. The root mean square
boundary of the target layer are considered. Then, error (RMSE) of the observed and calculated travel-
the traveltime residuals are plotted as a function of times for reconstructed model is 1.27 ms. We can see
the distance between a source and a receiver. A that the tomogram provides almost correct velocity
straight line is fitted to those points by the least model, but the boundary of each layer is not so clear.
squares method. Finally, we assume that the slope of We tried to extract true layer boundaries from it. Given
the straight line thus obtained is the slowness the approximate velocities and the number of layers, the
correction value for the layer. tomogram was converted into a layered model by
6. The iterative corrections of the velocity model are extracting some of the velocity contours. In the actual
conducted by using ray tracing and SIRT. data, the velocities and the number of layers can be
estimated by observed traveltime curves. Fig. 8 shows
the initial layered velocity model. The initial layered
3. Numerical examples velocity model is updated by repeating ray tracing and
SIRT, in which both the boundary depth and the
To examine the validity of the method, we carried out velocity of each layer are iteratively corrected. The final
numerical experiments. result is shown in Fig. 9. The calculated traveltimes for
Fig. 4 shows the velocity model used in the first the result and the observed ones are shown in Fig. 10.
numerical experiment. The traveltimes calculated for the The RMSE of the observed and calculated traveltimes
model are used as observed data. The initial velocity for reconstructed model is 0.46 ms. It can be seen that an
model is shown in Fig. 5. The model is represented by accurate velocity model is obtained and that the
21  15 (horizontal  vertical) rectangular cells. The residuals become sufficiently small.
comparison of the calculated traveltimes for the initial The more complex the underground structure, the
model and the observed ones are shown in Fig. 6. We more difficult it becomes to obtain true velocity
810 K. Hayashi, T. Takahashi / International Journal of Rock Mechanics & Mining Sciences 38 (2001) 807–813

Fig. 6. Calculated traveltimes for the initial model and observed traveltimes. There is large difference between observed and calculated traveltimes
for the initial velocity model.

Fig. 7. Result of the tomographic reconstruction. The RMSE of the Fig. 9. Final result. The initial layered velocity model is updated by
observed and calculated traveltimes for reconstructed model is 1.27 ms. repeating ray tracing and SIRT, in which both the boundary depth and
the velocity of each layer are iteratively corrected.

surface and within boreholes. The method can handle

both sources and receivers in boreholes. Here, we
present a numerical experiment using both surface and
borehole data. Fig. 11 shows the velocity model used in
the numerical experiment. There is a vertical zone of low
velocity in the deepest layer, which may be considered as
a fault or a fractured zone in the rock in case of actual
field applications. Several receivers are placed within a
borehole located close to the low velocity zone. The
traveltimes calculated for the model are used as
Fig. 8. Initial layered velocity model. The smooth velocity model
observed data. The initial velocity model is shown in
(tomogram) is converted into the initial layered model by extracting
some of the velocity contours. Fig. 12. The initial velocity model is constructed so that
the velocity is simply increasing with depth. The
resultant velocity model by automatic correction with-
structure from a surface seismic refraction method. This out the borehole data is shown in Fig. 13, and the result
difficulty can be overcome, if we place sources and using the borehole data together with the surface data is
receivers within boreholes. We have developed a joint shown in Fig. 14. It is difficult to identify the low
analysis method using data obtained both on the ground velocity zone from Fig. 13. On the other hand, we can
 K. Hayashi, T. Takahashi / International Journal of Rock Mechanics & Mining Sciences 38 (2001) 807–813 811

Fig. 10. Calculated traveltimes for the final result and observed traveltimes. The RMSE of the observed and calculated traveltimes for reconstructed
model is 0.46 ms.

Fig. 11. Velocity model used in the second numerical experiments. Fig. 13. Resultant velocity model by automatic correction without the
borehole data.

Fig. 12. Initial velocity model. The initial velocity model is constructed Fig. 14. Resultant velocity model by automatic correction using the
so that the velocity is simply increasing with depth. borehole data together with the surface data.

4. Applications to civil engineering investigations

easily see the low velocity zone from Fig. 14. The
numerical experiment shows that joint analysis using We show two examples of applications of the high
surface and borehole data is very useful for the resolution seismic refraction method to civil engineering
delineation of complex geological features. investigations.
812 K. Hayashi, T. Takahashi / International Journal of Rock Mechanics & Mining Sciences 38 (2001) 807–813

Fig. 15. Resultant velocity model obtained by conventional analysis method.

Fig. 16. Resultant velocity model obtained by the new method and the execution record.

Fig. 17. Resultant velocity model obtained by the conventional analysis method without the borehole data.

In the first example, the method was applied to a velocity lava inserted in the low velocity tuff breccia
tunneling site in Tertiary sedimentary rocks. Fig. 15 layers. The objective of this survey was to estimate the
shows the velocity model obtained by the conventional features of the lava. As long as sources and receivers are
analysis method (reciprocal method), and Fig. 16 shows placed on the ground surface, it is impossible to image
the velocity model obtained by the method presented the low velocity zone under the high velocity layer.
here, as well as the execution record. In this example, the Therefore, in this survey we placed hydrophones within
sources and receivers only on the surface are used. a borehole drilled at the center of the survey line
During the execution, the tunnel face collapsed at two together with geophones on the ground surface to
places where the seismic velocity is smaller than that of acquire waveform data both on the surface and within
its vicinity. However, velocity differences are not so the borehole simultaneously. Twenty-four hydrophones
large. This means that high resolution and high accurate with 2 m depth intervals and 101 geophones with 5 m
velocity distribution must be obtained to predict the distance intervals were used for the data acquisition.
tunnel face collapse before excavation. Therefore, the The sources were only on the ground surface with
high resolution seismic refraction method can play an approximately 30 m distance intervals and dynamite was
important role for tunnel route investigations. used for the sources. Fig. 17 shows the resultant velocity
The next example is the application of the method to a model obtained by the conventional analysis method
dam construction site which has a thin layer of high (reciprocal method) without the borehole data. Fig. 18
 K. Hayashi, T. Takahashi / International Journal of Rock Mechanics & Mining Sciences 38 (2001) 807–813 813

Fig. 18. Resultant velocity model obtained by the new method using the borehole data together with the surface data.

shows the result obtained by the method presented here, 4. Numerical and actual field examples prove that the
using the borehole data together with the surface data. method can provide the reliable velocity model, even
In the inversion, geological information obtained from for sites with very complex geological features, for
borehole is used, and the assumption that the velocity is example, vertical structures and high velocity thin
increasing with depth is not used in this example. We layers inserted in a low velocity medium.
can see the very detailed feature of high velocity lava in 5. Considering that the demand for investigations for
between the low velocity tuff breccia from Fig. 18. It complex geological feature is increasing, the method
must be noted that the resolution and reliability of the presented here can play a very important role
result is not same through the velocity model. The increasingly for such investigations.
resolution and reliability in the near surface and the
region around the borehole may be relatively high.
References

[1] Hagiwara T, Omoto S. Land creep at Mt. Tyausu-yama

5. Conclusions (Determination of slip plane by seismic prospecting). Bull Earth-
quake Res Inst 1939;17:118–37.
In this paper, we present a new analysis method for [2] Hagedoorn JG. The plus–minus method of interpreting seismic
seismic refraction data. In conclusion refraction sections. Geopys Prosp 1959;7(2):158–82.
[3] Palmer D. An introduction to the generalaized reciprocal method
1. The method enables us to analyze seismic refraction of seismic refraction interpretation. Geophysics 1981;46:1508–18.
[4] Moser TJ. Shortest path calculation of seismic rays. Geophysics
data almost automatically after picking first arrival 1991;56:59–67.
times, with high resolution and high accuracy. [5] Saito H. Traveltimes, ray paths of first arrival seismic waves:
2. The method is applicable to both layered velocity computation method based on Huygens’ principle. 59th Annual
structures and those having a velocity gradient with International Meeting, Society of Experimental Geophysics,
Expanded Abstracts, 1989, p. 244–7.
depth.
[6] Hayashi K, Saito H. Automatic correction of velocity structures in
3. The method can handle the data obtained from seismic refraction prospecting. 66th Annual International Meeting,
sources and receivers on the surface of the ground, as Society of Experimental Geophysics, Expanded Abstracts, 1996,
well as those within boreholes simultaneously. p. 1658–61.