A Review of Receiver Function Techniques for Estimation of One-Dintensional Velocity Structures

Two receiver function techniques both broadly used by seismologists to estimate one-dimensional shear wave velocity structures of the crust and upper mantle beneath seismic stations have been evaluated. One employs a deconvolution filtering,. which is directly accomplish�d in the time domain (Vinnik, 1977; Kosarev et aL, 1987). The other is completed through a source. . equalization, which is perf or1ned by spectra division in the frequency domain (Langston, 1977; Owens et al., 1984). In this study, the performance of these techniques is examined from two sets of synthetic seismograms that are computed from two one-dimensional models by using Thomson and Haskell's method (Haskell, 1962). The results suggest that both techniques can mostly recover the ass11med models very well when the energy of multiples (or crustal reverberations) in the direction of P-wave propagation is minor. If the model is more complicated or the multiple energy becomes stronger, however, the results from Vinnik-Kosarev technique appear to be better thap the other in both modeling the converted P-SV phases and inverting the structures. Combining the results of the synthetic tests and theoretic comparisons, it is concluded that the differences resulting from both techniques are primarily caused by the processes through difTerent domains. Besides, both techniques consistently indicate that the inversion results are dependent upon the incident angles. For the nearly vertical incidence, neither technique could resolve the models very well. (


INTRODUCTION
The modeling of converted phases of P-SV type from teleseismic wavefonns, known as receiver function analysis, has often been used to detennine one-dimensional shear wave velocity structures of the crust and upper mantle beneath a seismic station.One main process 182 TAO, Vol.6, No.2, June 1995 of these analyses is to recover the converted phases by removing the source and deep man tle effects from tele-seismograms.These converted phases are generated at the boundaries between homogeneous layers of the crust and upper mantle which the primary waves pass through.To carry out this process, two different kinds of techniques have been indepen dently developed.One is the source-equalization that was fi rst proposed by Langston (1977) to deal with long-period teleseismic P-waves collected from several Wo rld-Wide Standard Seismograph Network (WWSSN) stations.Then Owens et al. (1984) and Owens (1987) incorporated this technique with an inversion scheme in the time domain to recover detailed infonnation present in the broadband teleseismic P-wavefonns recorded at several stations of the Department of E nergy Regional Seismic Te st Network (RSTN).The other technique, initially reported by Vinnik (1977) and later modified by Kosarev et al. (1987;1993), is the procedure of using the, time domain deconvolution fi ltering of Berkhout ( 1977).
This study is undertaken in order to evaluate the two techniques mentioned above, and it is based on synthetic tests and theoretic comparisons.At the beginning of this study, both technical procedures of obtaining receiver functions are briefly described.Then both tech niques are examined by modeling synthetic seismograms that were created by the Thomson Haskell method (Haskell, -1962).Finally, the results of the synthetic tests along with some theoretic comparisons are discussed.

METHODOLOGY
2.1 Langston and Owens Technique Langston (1979) and Owens et al. (1984) used a source-equalization procedure to determine the velocity structure of the crust and upper mantle beneath a seismic station (Figure la).They stated that a teleseismic P-wave, D(t), recorded at a station can be theoretically expressed by : Dv (t) ==: I(t) * S(t) * E v (t), DR(t) I(t) * S(t) * ER(t), Dr(t) I(t) * S(t) * Er(t), (1) where the subscripts V, R, and T are the vertical, radial, and tangential components, respec tively; I( t) is the instrument response of the recording system; S( t) contains the effects of the seismic source and deep mantle, .and E( t) is the effects of the crust and upper mantle structure.Asterisks denote the convolution operator.
The main purpose of the source-equalization process is to isolate the response of crust and upper mantle structures from the other factors that interact with it to form the observed seismograms recorded at teleseismic distances.The most important assumption made in this technique is that E v( t) = 6( t ), which means crustal multiples (reverberations) and converted phases on the vertical component of steeply incident P-waves are considered negligible.Thus, components, the source and deep mantle effect in the radial component can be removed by deconvolving the vertical components, Dv (t), from the radial component, DR(t).
As a matter of fact, the procedure for performing the deconvolution is to divide the spectrum of the radial component by that of the vertical component in the frequency-domain\, The process is given by:

Dv(w)
(2) To reduce the• effects of noise generated by dividing by very small numbers due to troughs in the spectrum, a substitute water-level process (Dey-Sarkar and Wiggins, 1976) is given by: where and D\,(w) is the complex conjugate of Dv(w).The function of F88(w) can be thought of as the auto-correlation of Dv( w) with any spectral troughs filled to a level as deterrnined by the water-level value ' 'c''.The parameter ' 'a'' is used to control the width of the Gaussian function.This substitute Equation (3) introduces the minimum allowable amplitude level of the amplitude spectrum of the vertical component.
After spectral division, the resulting spectrum is transforrned back to the time domain• • to obtain the observed receiver function of the radial component, ER(t).Meanwhile, the calculated receiver function is generated by applying the source-equalization procedure de scribed above to the synthetic seismograms that are calculated by a fast ray-tracing scheme (Langston, 1977).Finally, one-dimensional shear wave velocity structures are inverted by minimizing the difference between the observed receiver function and the calculated one.

Vinnik-Kosarev Technique
Vinnik (1977) and Kosarev et al. (1993) considered seismograms in an L, H, and T coordinate system, where L is along the P polarization direction, H is in the plane of propagation and nonnal to L, and Tis normal to both H and L (Figure lb).This coordinate system •is chosen to optimize the detection of SV.These authors also assumed the instrument response to be the same in the three components.Thus, the response of the teleseismic P-wave, D(t), can be presented as: . DL(t) -S(t) * EL(t) DH(t) S(t) * EH(t), (4) where S( t) is the source of the P-wave in the half-space, E( t) is the response of the layer structure, and asterisks are the convolution operator., Once the seismograms are rotated, three steps are followed.The first step is to generate a deconvolution filter, F( t), in the time domain by minimizing the least-square difference between the output of the filtering operating on the L component and a spike-like function of the normalized amplitude.This deconvolution filter i • s essentially equivalent to the Wiener filter (BerKfiou• t;--•1 CJ7 7)," • --anaits:-c-purpose.IS. -to-recovef--tfie• -corifrtoutions of the-P-SV.• phase.
conversions to the H .. component.The second step is the application of tll is deconvolution filter to the H component,. which results in the response of the medium in • the H direction to.a no11nalized spike in the L direction.This step is valid on • ly if the L component is not contaminated �significantly by multiples from near-surface structures, which means the structure response in the L component, EL ( t ), is .considered as a spike-like function.The final step is to invert the structure by minimizing the difference between the observed receiver functions and _ the c�mputed responses of :the medium.,

SYNTHETIC TESTS
To evaluate thy perfonnance of both techniques mentioned above, two groups of syn thetic seismograms have been created using the Thomson-Haskell technique (Figure 2).The synthetic seismograms are constructed from two assumed one-dimensional models (Table I).
The Ml model is composed of three layers with a half space, and velocities gradually increase with _ q�pths.In contrast, the more complicated M2 model is composed of seven layers with a half Wace, including a low velocity layer at depths between 6 and 10 km.These synthetic 1seismograms are the structure responses to a spike-like function of P input with a sampling interval of 0.2 seconds� The width of the spike is 2 • seconds.Each group includes four pairs of synthetic seismograms corresponding to four representative incidences (I=l 0, 10°, 20° and 30°), respectively. -•-.

•
The synthetic seismograms clearly show that the energy of multiples (or reverberations) and converted P-S phases is strongly dependent upon the incident angles of the input seismic waves as well as the models (Figure 2).The energy of the converted P-SV phases signifi cantly increases with the incident angles.For the 30° incidence, the maximum amplitude of the P-S converted phases reaches about l 0% of the first P-wave.For the nearly vertical inci dence (1=1°), on the other hand, there is • almost no energy of the P-SV converted phases in the H-component.Moreover, the energy o .f multiples decreases as the inc • ident angle becomes Cheng-Horng Lin   horizontal.To illustrate, the energy of multiples in the L components from I 0 incidence is about twice that fr om 30° incidence for the case of model Ml.
Figure 2 also shows that the multiples and converted P-SV phases computed from the two models are quite different in not only shape but also amplitude.For example, the average energy as shown in Figure 2a is only about half of that in Figure 2b.Besides, it should be noted that the maximum energy of the multiples in the L-component is about 10-15% of the first P-wave arrivals, depending on the model and incidence angle.These results are not suiprising because the reflected or transmitted energy from each boundary of the model is dete11nined by the velocity and density contrast between the layers.

Tests of the Langston-Owens Technique
To test both techniques, the synthetic responses as shown in Figure 2 are taken to be seismograms observed at surface.At first, the observed _ receiver functlons in the radial component were extracted from these seismograms using the source-equalization procedure (Langston, 1977;Owens et al., 1984).Then these observed receiver functions are expected to be exactly the same as the synthetic responses on surface because the input function at the bottom of the models is equivalent to a spike function.
The result of the source-equalization procedure shows that the Langston-Owens tech nique can mostly recover observed receiver functions in the radial component, except at some multiples or converted P-SV phases (Figure 3).In the case of the nearly vertical incidence, (l=l 0) for example, main differences are found at the the arrival of the first multiple in which a reverse amplitude between the observed receiver functions and synthetic seismogramshas been seen at 5 sec in Figure 3a.Besides, several slight differences are seen at the later ar rivals of the multiple or converted P-SV phase.These differences become even more visible when the model is more complicated.For instance, there are some reverse amplitude at the about 7 sec in Figure 3b.Such results suggest that the observed receiver functions have been contaminated by the energy of multiples that have not been properly considered in the source-equalization process.
The inverted models are also compared with the true models that are used to create the synthetic seismograms (Figure 4).The results of the inverted models depend on the complexity of the models and the incident angle of the seismic waves.For the simple model (M 1 ), consisting of three layers over a half space, the inverted models are almost the same as the true models, and the receiver functions match well (Figures 4a and 5a).When the models are more complicated, however, the receiver functions do not match very well (Figure 5b).Furthermore, the differences between the inverted and true models (Figure 4b) become remarkable even though the general trend can be recovered.At the upper part of the models, the low velocity layer at depths from 6 to 10 km is hardly resolved.At the lower part of the models, the true models are not even limited within the error bars of the inverted models.It is believed that these mis-fits appear to have been primarily caused by the contamination of multiple energy in_ the observed receiver functions.187 The inversion res�lts vary significantly depending on the incident angle (Figure 4).In general, the fit of the inversion models is improved with the increase in the incident angles.One possible explanation for this fe ature is that when the incidence ray is close to the vertical direction, the energy of the converted P-SV phases is small, and then it is easily contaminated by the neglected multiple energy in the observed receiver functions.

Tes ts of the Vinnik-Kosarev Technique
The synthetic seismograms (Figure 2), used for testing the Langston-Owens technique above, are also employed for testing the Vinnik-Kosarev te chnique in this section.In the same way, it is expected that the receiver functions in the H-component are equivalent to
Model M1 The deconvolution results in the H-component of the Vinnik-Kosarev technique can mostly recover most of the receiver functions (Figure 6).Although there are also some visible differences between the synthetic seismograms and observed receiver functions, one can not find at least any of reverse amplitudes between both seismograms.Therefore, it is believed that these results are better than those computed by the source-equalization procedure in the Langston-Owens technique.
The inversion results show that the models are generally resolved quite well even though the fit is slightly affected by the complexity of the model.For the simple model (Ml ), the inverted models are similar to the true ones (Figure 7a).For the complicated model (M2), the inverted results are only slightly different from those of the true models (Figure 7b ).

•
The calculated seismogram in the L-component is almost equivalent to that in the deconvolved one which is taken as the wavefo1n1 of the incoming P-wave (Figure 8).In other words, the structure response in the L-component can be considered a delta function or Ei(t)=8(t).-Besides; the receiver functions'in the H�component match very• well (Figure 9).These results imply.most energy of mu, J: tiples, :which has signifi cantly contaminated the receiver functions in the • Langston-Owens .technique,.may have already been removed or suppressed by the deconvolution processing in the Vinnik-Kosarev technique.
---------r---- •2 0 2 4 6 8 10 12 14 16 18 20 22 24 . .,_ .A brief summary that can be drawn from the te�ts is that the Vinnik-Kosarev technique can mostly recover the structures from all of the test models.On the other hand, the Langston Owens technique can successfully deal with simple structures.When the structures become complicated, care must be taken in using this technique because the calculated receiver functions for inverting the structures may be contarninated by multiples.

DISCUSSION
Although the results of the synthetic tests show that both the Langston-Owens and Vi nnik-Kosarev techniques can mostly recover most of the true models, some visible differ ences have been found in their receiver functions and inverted models.In order to understand   (1) The selection of coordinate systems: Langston-Owens chose a traditional coordinate system (R and Z) to easily deal with most vertical angles of incidence.In fact, the incident angle of the teleseismic P-wave may not, however, be vertical.On the other hand, the coordinate system (L and H) chosen by Vinnik-Kosarev is more optimal in distinguishing the particle motions between the Pand SV-waves regardless of the incidence angle of the teleseismic P-waves.However, it is felt that the noticeable difference in the inversion results of the synthetic tests could hardly be caused -3 -1 Incidence= 20 .. that this process is generally quite sensitive to noise in the input data.In the analysis of receiver functions, unfortunately, one of the most remarkable noises is due to multiples which are always neglected because their energy on the vertical component of steep incident P-waves is minor.In other words, it may be expected that the structural response in the vertical direction is a spike-like function or that Ev(t)=8(t).From the synthetic seismograms, however, it has been found that the energy of crustal multiples and converted phases strongly depends on the incidence angles and the complexity of structures (Figure 2).The energy of maximum multiples can reach up to about 10-15% of the energy of the first P-wave arrival.Thus the assumption that Ev(t)=b(t) is too simple to fit the particular cases where multiple energy is not considered as minor .

193
Therefore, this simple assumption of minor multiples probably causes some problems in the receiver function analysis under the Langston-Owens technique (Figure 3).In addition to the multiples ignored in the receiver function analyses, some high fre quency noise may be introduced due to the process of discrete Fourier Tran sf onn (Brigham, 1974).For example, the sampling rate produces the aliasing problem, and the truncated length in the time domain causes rippling in the fr equency domain.In the Langston-Owens technique, these noises may be filtered out by using the adjustable water-level process (or the Gaussian smoother).In the meantime, however, some origi nal convolution signals may have truely lost all information at these fr equencies where the spectra are below the water-level value.
Just like the Langston-Owens technique, the Vinnik-Kosarev technique assumes that the energy of crustal multiples is minor (Kosarev et al., 1993).In the synthetic tests above (Figure 6), however, the inverted structures seemingly have not been contaminated by multiple energy because the energy of • multiples has been removed or suppressed in the calculation of the deconvolution filter (a Wiener filtering) in the time domain.As a result, the inverted structures are usually better than those obtained by using the source-equalization process through the frequency domain.
(3) The forward calculation: the Vi nnik-Kosarev technique uses the Thomson-Haskell method to calculate the exact receiver functions to invert the structures.To save com puter time, on the other hand, the Langston-Owens technique employs a substitute method, a Gaussian smoother, to find the receiver functions.Although this substitute method does not generate the complete response of the structure, the appropriateness of this choice has been verified by comparing it with the alternative in the Thomson Haskell method (Owens et al., 1984).

S. CONCLUSIONS
The results of synthetic tests and theoretical comparisons reveal that the inversion of teleseismic wavef orrns to estimate the one-dimensional shear wave velocity structures beneath a seismic station can be performed equally well by both the Langston-Owens and Vinnik Kosarev techniques if the energy of multiples on the direction of P-wave propagation is minor.However, when the structures are more complicated or multiple energy becomes stronger, the Vinnik-Kosarev technique appears to be better than the other.One of the reasons for this is the process in the different domains.The process of a source-equalization, which is performed by spectra division in the frequency domain, is very sensitive to the multiples that have been ignored in the analysis of the receiver functions.On the other hand, the process of a deconvolution fi ltering in the time domain can successfully handle multiples.In any case, both techniques consistently indicate that the inversion results are strongly dependent upon the incident angle of a seismic wave.In fact, for the nearly vertical incidence, neither technique can resolve the models very well because the energy of the converted P-SV phases is very small.

Fig. 1 .
Fig. 1.Simplified models to indicate the crust and upper mantle structures and the •incidence of teleseismic waves: (a) coordinate system (V, R, T) used in the I ,angston-Owens.technique, and (b) coordinate system (L, H, T) used in the Vinnik-Kosarev technique.

Fig. 3 .
Fig. 3. Synthetic seismograms in the radial component (solid lines) generated from (a) the Ml model, and (b) the M2 model.Dotted lines show the observed receiver functions in the radial component obtained from the source-equalization procedure.
functions in the radial component (solid lines) and . .: calculated receiver function� (dotted .lines) obtained from the sourcef the Langston-Ow�ns techn i que: (a) the Ml model, • . .and (b) the M2 model.
the models• can be recovered not orily at the low velocity layer but also at the lower part of•the models .•This has not been successfully : resolved by • the Langston-Owens te chnique . .: Therefore, • it is • .concluded, .that the results computed .from the Vinnik-Kosarev te chnique are clearly better than those from the Langston� • ow ens technique.

Fig. 6 .
Fig. 6.Synthetic seismograms in the ff-component (solid lines) generated from .(a) the Ml model, and (b) the M2 model.Dotted lines show the observed receiver functions in the H-component from the deconvolution procedure of Vinnik-Kosarev technique.
• by the selection of different coordinate systems.(2)The processing domain: The source-equalization of the Langston ... Owens technique is perfo1•1ned by a process of spectra division in the frequency domain.It is well known ------� ��--�.

Fig. 8 .
Fig. 8. Deconvolved seismograms on the L-component (solid lines) and calcu-, lated responses (dotted lines) by using the Thomson-Haskell technique: (a) the Ml model, and (b) the M2 model.

Fig. 9 .
Fig. 9. Observed receiver functions on the H-component (solid lines) and calculated receiver functions (dotted lines) obtained from the Vinnik-Kosarev technique: (a) the Ml model, and (b) the M2 model.

Table I .
Tw o one-dimensional models for generating synthetic seismograms. .