CC BY
13
RUSSIAN JOURNAL OF EARTH SCIENCES, VOL. 11, ES1007, doi:10.2205/2009ES000415, 2009
PROCEEDINGS OF THE INTERNATIONAL CONFERENCE
Electronic Geophysical Year: State of the Art and Results
3-6 June 2009 • Pereslavl-Zalessky, Russia
Use of distributed computing systems in seismic wave form inversion
I. M. Aleshin,1 M. N. Zhizhin,2 V. N. Koryagin,1 D. P. Medvedev,2 D. Yu. Mishin,2 D. V. Peregoudov,1 and K. I. Kholodkov1
Received 15 October 2009; accepted 16 November 2009; published 24 November 2009.
Seismic anisotropy presents a unique possibility to study tectonic processes at depths inaccessible for direct observations. In our previous study to determine the mantle anisotropic parameters we performed a joint inversion of SKS and receiver functions waveforms, based on approximate methods because of time consuming synthetic seismograms calculation. Using parallel calculation and GRID technology allows us to get the exact solution of the problem: we can perform direct calculation of cost function on uniform grid within model parameter space. Calculations were performed for both synthetic models and real data. It is shown that the application of the joint inversion of SKS and receiver function from the one hand improves resolution for the determination of base anisotropic parameters, from the other hand requires careful analysis of the consistence of different groups of data. Ignoring the possible disagreement of different groups of data can lead to significant errors in the estimation of anisotropies parameters. KEYWORDS: seismic anisisotropy, SKS spliting, receiver function, waveform inversion, distributed calculation, grid.
Citation: Aleshin, I. M., M. N. Zhizhin, V. N. Koryagin, D. P. Medvedev, D. Yu. Mishin, D. V. Peregoudov, and K. I. Kholodkov (2009), Use of distributed computing systems in seismic wave form inversion, Russ. J. Earth. Sci., 11, ES1007, doi:10.2205/2009ES000415.
Introduction
The investigation of seismic anisotropy gives us the opportunity for direct study of the mantle structure [Nikolas and Christensen, 1987]. Currently seismic waves of different types are used for studying the anisotropic mantle parameters [Dziewonski and Anderson, 1981; Babuska et al., 1984; Kosarev et al., 1984; Vinnik et al., 1984]. The most known method based on the analysis of SKS waveforms and related phases was first used by Vinnik et al., [1984] and reviewed by Savage, [1999]. If there are anisotropic rocks on the path from the core-mantle boundary to the receiver, the transverse wave splits inside the rocks into two quasi-transverse waves, traveling with different speeds. The time delay of one wave relative to another is formed inside the anisotropic layer, and outside of it both phases travel with the same velocity, that makes impossible to determine the absolute
1Institute of Physics of the Earth RAS, Moscow, Russia
2Geophysical Center RAS, Moscow, Russia
Copyright 2009 by the Russian Journal of Earth Sciences. http://elpub.wdcb.ru/journals/rjes/doi/2009ES000415.html
depth of the anisotropic layers [Menke and Levin, 2003].
Another method of mantle anisotropy investigation is using the P ^ S exchange on the environment heterogeneities. The direct and converted waves are traveling with different speed after leaving the layer, that allows us to measure the depth of the anisotropic layer. The first usage of converted waves for the investigation of mantle anisotropic parameters was described by Kosarev [1984]. Further progress of this method was made by Girardin and Farra, [1998], Vinnik et al., [2002, 2007]. The method's drawback is associated with the low amplitude of the converted waves being weakened by both the small anisotropic parameters and the small contrast of the heterogeneities velocities (only the first factor weakens the SKS).
The combined use of SKS waveforms and receiver functions of converted waves was suggested for the investigation of the Earth anisotropic parameters [Vinnik et al., [2002]. It is based on the assumption that the different observed effects in both data sets have the same nature - regions at 250-300 kilometers depth with anisotropic rocks (this depth limitation is caused by the converted waves allocation method). In this case both the SKS and converted waveforms can be explained by the same anisotropic model.
ES1007
1 of 9
Figure 1. Tangential component of converted waves (panels a, c) for different azimuths and its azimuthal filtration result (panels b, d). The top panels (a, b) - illustrate the azimuthal filtration procedure for data got as a result of model calculation (Table 1), the bottom panels (c, d) - results for the same data but with gaussian noise addition.
In a simple case the model under the station can be presented as a stack of homogeneous layers on a homogeneous isotropic half-space. To investigate the model parameters we need to perform waveform inversion. However even for the small amount of layers the time needed to browse the whole models space is quite large. One of the ways to browse the whole parameters space is to use a distributed computing systems. We used the GRID - geographically distributed infrastructure that combines a large amount of different resources (processors, long-term and operating memory, storages and databases, networks) that allows access for users from any place and any location. Computational GRID is oriented at the applications running and controlling the big amounts of parallel computations on distributed computational clusters. The computations performed using the EGEE GRID-infrastructure [Renard et al., 2009; Foster and Kesselman, 1999]. Initially the EGEE project provided the computational resources for the analysis of data from the Large Hadron Collider (LHC) in CERN, Geneva. The research we presented is the first attempt of geophysical computations using the EGEE infrastructure supported by the eEarth virtual organization, the Russian virtual organization for geophysical computations in GRID.
Another aspect is associated with the simultaneous use of
the two data types. In this case two objective functions appear. So it's considered to be a vector function. The problem of vector function minimization is a generalization of a scalar function (e.g. [Sakawa, 1993]). In general the multipurpose optimization goal is to determine the Pareto parameters, when every component of a multipurpose function can't be improved without worsening the another components. However in our case we should not use this method. A vector multipurpose function can be used when there are different contradictory goals. If two data groups can be described by one model (what is meant), then the both components of the objective function would reach the minimum simultaneously. The differences can be caused by the measurement errors. In our case we can speak about compatibility of different data types rather than choosing the compromise model.
Problem Statement
The procedure of the data processing is described elsewhere [Vinnik et al., 2007]. We assume that we have the normalized records of converted waves with good azimuthal
coverage (15-20 degrees) and SKS phases records from two-three different azimuths. This is typical for real observations.
Suppose the environment under the station is modeled by a stack of plane homogeneous layers, partly anisotropic, on the isotropic half-space Kosarev et al., 1979. Our goal is to choose the layer parameters to explain the observations made. The environment with hexagonal symmetry with a horizontal axis is used as a model of the anisotropic layer. To describe the hexagonal media we have use five elastic constants that can be presented by following physical parameters: average by direction velocities of transverse and longitudinal waves vS vp, anisotropy coefficients of transverse and longitudinal waves velocities as and ap, defined by correlations
v(p's) ^max
v(p,s)/ (l + 1/2a(p,s))
>'s) (1 + a(p,s)) ,
(1)
viiaX - maximum and minimum wave velocities. The
last anisotropic parameter n is considered to be equal 1.03. Having analyzed the data published by BenIsmail and Main-price, [1998], we assume that the anisotropy coefficients are linearly dependent: ap = 1.5as. We consider the maximum velocity to conform the hexagonal axis (a(p,s) > 0). Thereby the variable parameters of the model are layer thickness d, anisotropy coefficient of S-waves as and anisotropy axis azimuth 0.
The important feature of the chosen model is the n-periodicity of the waveforms by the azimuth. It can be used to separate the weak anisotropic effects from noise and lateral heterogeneity effects. For that we used az-imuthal filtration of records. Assume the measured seis-mograms Sobs(t, ifi) for azimuths set ipi. Then the filtering procedure can be put down as
ffots (t,i>) = ^ Sobs (t, (pi)g('pi - i>), (2)
g(x) — n-periodic function, N = ^ ■ g2(ifii — — normalization coefficient. The g(^) function choice is arbitrary. It is shown by Girardin and Farra, [1998] that we can choose g(^) just as a trigonometric function of double angle.
As we mentioned in the introduction, the anisotropic part of the converted waves P ^ S is very weak. It makes the azimuthal filtering necessary. As an example there are synthetic SH-components of the converted waves for different azimuths with noise addition and the result of azimuthal filtering using the formula above in Figure 1.
In the first works making use of SKS waves the azimuthal waveforms filtering was also performed [Kosarev et al., 1984]. However the difficulty and often impossibility of the observation of SKS waves with good azimuthal coverage forces to use only two-three records for the environment anisotropic parameters determination [Savage, 1999]. It prevents us from the ability to make azimuthal filtering for excluding the noise and possible inclined boundaries affection.
Let us assume the objective functions as an average quadratic deviation of converted waves synthetic seismograms and SKS waves from the corresponding observed data. The
Figure 2. GRID cluster tasks running scheme.
problem is to find the relation of the objective functions to the model parameters. If both function minimums have a shape of an extended "valley" then the common usage of exchange and SKS waves is meaningful if the "valleys" cross at a considerable angle.
We used multiprocessor systems for objective function values calculation. This allowed minimizing the processing time. The calculations on a loosely coupled cluster assume that the task is split into blocks calculated without any communication between each other and sequence independent. In our case every block included the objective function values calculation for a number of points from the models space. The amount of points in a block affects the block startup time on the cluster. In our case the calculations on each model were split into 67 blocks containing 100,000 tasks each (Figure 2). This scheme allowed the balanced loading of all cluster nodes and made the startup time of each block considerably less that the block processing time.The blocks amount allowed loading all the available nodes of the computing cluster.
To run the test tasks we used the Geophysical Center cluster. It has one computing element and two working nodes with two processors each. The final calculations were performed on SINP MSU cluster. Currently it contains three computing elements: lcg02.sinp.msu.ru - 48 processors (AMD), lcg06.sinp.msu.ru - 44 processors (Intel XEON) and lcg38.sinp.msu.ru - 76 processors (Intel XEON).
One model calculation time in GC RAS was 16 hours, SINP - 7 hours.
Synthetic Examples
To solve numerically the minimization problem we need to calculate the objective function values on a grid. We chose the following discretization step values: anisotropy axis azimuth — 5 degrees, anisotropy coefficient 0.01, layer thickness — 10 km. Considering the geophysical task applications, we should mention the different importance of the parameters. The most important is the anisotropy axis direction in the layer allowing making the important tectonic conclusions [Nikolas and Christensen, 1987]. The anisotropy strength or the anisotropy layer thickness are not so important and their variations much less affect the objective function.
Table 1. Model parameters for the first synthetic example. The model consists of three layers - one isotropic and two anisotropic - in the isotropic half-space.
Layer number Vp, km s 1 Vs, km s 1 P, g cm Thick, H, km Azimuth, 0, deg . Anisotropy coefficient as
1 6.4 3.7 2.9 40 — 0.00
2 8.1 4.5 3.3 60 30 0.04
3 8.1 4.5 3.3 100 100 0.02
4 8.5 4.72 3.4 TO — 0.00
To analyze the data we can construct the function that reflects the models distribution by anisotropy axes directions independently from the other parameter values. We consider the base model with two anisotropic layers. The distribution function of anisotropy axes azimuths is:
qa(0i,02)= ^ e (Ca - Ca(m)). (3)
m6MOi >2)
e(x) - Heaviside function, m - models space vector, M(01, 02) - models subspace with fixed anisotropic axes azimuths values, Ca - objective function value for data group a = SKS, RF. Ca values were chosen empirically.
Let's start the research from the synthetic example with the parameters in Table 1 (this model was used in [Vinnik et al., 2007]).
The synthetic seismograms of converted waves (18 azimuth values) and SKS waves (3 azimuth values) were calculated for this model. Using it as "observations data" we have got the vector objective function values (CRF,Csks). It is important that we used the three-level model with two anisotropic layers to calculate these functions. The results analysis shows that in this (ideal) case the global extremum of each component of the objective function match and correspond to the original model (Figure 3).
The picture shows the objective function of converted waves is more sensitive to the first layer anisotropy axis direction change. It could be caused by the bigger (compared to the second layer) contrast on its top boundary. The data objective function shape improves the localization os summary extremum when used with converted waves data though.
Figure 3. Azimuth distribution function derived from waveforms inversion results for the model in Table 1. The functions maximums of distribution functions by converted and SKS waves azimuths are marked with "+" and "x" symbols accordingly.
Figure 4. The same as on Figure 3 with addition of the gaussian noise to the source synthetic traces.
In the next experiment we added the gaussian noise convolved from the longitudinal component spectrum to the converted waves signal (Figure 1). Figure 4 shows the inversion procedure also gives us the answer almost coincident with the truthful.
It becomes worse if we have a high depth anisotropic layer. As we mentioned in the introduction even if we consider the absence of the inclined boundaries on the SKS wave path, it summarizes the information about the anisotropy along the path from the kernel-mantle border to the Earth surface. To illustrate the effects of ignoring this fact let's examine the model in Table 2.
The upper part is the same as the previous model, but two layers are added to the lower part: one is isotropic, another is anisotropic. This addition can't seriously affect the converted wave waveforms in the time interval from 0 to 20 seconds. But change of the SKS signal is significant. As
in the previous example, we used the synthetic seismograms built for this model as a source data. It is important that the inversion was made using the two-layer model, as before.
Figure 5 shows the results of the calculations. As we expected, the minimum of the converted waves objective function still matches the "old" model (or the upper part of the new model), and the matching extremum for SKS moved significantly. If we minimize the scalar objective function (wa are weighting coefficients)
^ WaCa(m), (4)
a
like it was made by Vinnik et al., [2007], we will get the model marked with the sign "*" in Figure 5. The quality of both groups selection is quite good (Figure 6) that shows the potential danger of data misinterpretation.
Table 2. Model parameters for the second synthetic example. The upper part of the model matches the previous one. Two layers are added to the lower part, one is anisotropic.
Layer number Vp, km s 1 Vs, km s 1 P, g cm Thick, H, km Azimuth, ф, deg . Anisotropy coefficient as
1 6.4 3.7 2.9 40 — 0.00
2 8.1 4.5 3.3 60 30 0.04
3 8.1 4.5 3.3 100 100 0.02
4 8.1 4.5 3.3 20 — 0.00
5 8.1 4.5 3.3 30 60 0.03
6 8.5 4.72 3.4 TO — 0.00
Figure 5. Azimuths distribution function built using the waveforms inversion results for the model in Table 2. Symbols "+" and "x" matches the distribution functions maximums of the converted and SKS waves. The "*" symbol shows the minimum of the scalarized object function.
Figure 6. Fit quality illustration in minimizing the scalarized objective function (formula 4). The input data was calculated using the Table 2 model, but the model Table 1 was used for the inversion. The black lines are the input data and grey lines match the model marked in Figure 5 with the "*" symbol.
Figure 7. The azimuth distribution function derived from a result of the waveform inversion for the CHM station. The notifications are the same as in Figure 5. The minimum value of the objective function marked with "*" symbol was obtained in [Vinnik et al., 2002; Vinnik et al., 2007].
Real Data
Let us explore the data received in Tian-Shan anisotropy studies [Vinnik et al., 2002, 2007]. We chose two stations as an example: CHM (coordinates: 74.8o E, 43.0o N) and AKSU (80.1o E, 41.1o N). To explain the corresponding waveforms the most simple models were used in [Vinnik et al., 2007]: three layers (the isotropic first one and two anisotropic lower ones) on the homogeneous isotropic halfspace.
Using this model, the values of the objective function (Crf,Csks) were calculated on the grid described in the previous chapter. The calculation results for CHM and AKSU stations are shown in Figure 7 and Figure 8 accordingly. It's easy to see that for both stations converted and SKS waves data do not coincide each other. We can see some data groups matching for CHM station, at least for the second layer anisotropy direction, but for the AKSU station the objective function minimums of converted and SKS waves are reached in absolutely different models. Using the scalar multiplication (4) as an objective function gives a model not being the "average" of the parameters of the first and second data groups models. The relief features of the objective functions like narrow "valleys" make the optimal in scalar
multiplication model to be remote from both minimums. It is most noticeable for the AKSU station (Figure 8).
To conclude the above synthetic example, the result has an obvious explanation: the SKS wave was affected by the additional (compared to the converted waves) influence. It could be either a deep (deeper than 250-300 km) anisotropic layer or an inclined boundary, or another lateral heteroge-neousness.
Considering the analysis made we can conclude the impossibility of both data groups consistent interpretation for AKSU and CHM stations in the selected model boundaries.
Conclusion
As we can see from the first two synthetic examples, ideally the combined inversion of converted and SKS wave waveforms allows us improving the upper mantle anisotropic parameters determination. However, the third example shows, and the real analysis data convinces us that the lack of SKS data and the following impossibility of azimuthal filtering results in data interpretation errors. Therefore before making the common interpretation of such data we need to make
0 20 40 60 80 100 120 140 160 180
Azimuth, deg
Figure 8. AKSU station. The same with the on Figure 7.
sure the objective functions minimums of converted and SKS waves are reached in close models. Otherwise either the combined inversion should be refused or the enough amount of SKS waves records with a good azimuthal coverage should be available. Also the azimuthal filtering can be done to reduce the inclined boundaries and the affection of deep anisotropic layers.
The use of GRID technologies allowed us to make the calculations in acceptable time. The advantage of the technology used is also the use of the GRID monitoring tools to control the processing including the calculations finishing time. The method allows analyzing the more complex models by extending the amount of GRID nodes used. But further use of the method would require a user interface development that would automate the task start process and, more importantly, further analysis of the results.
References
Babuska, V., J. Plomerova, J. Sfleny (1984), Spatial variations of P residuals and deep structure of the European lithosphere, Geophys. J. R. Astron. Soc., 79, 363.
Benismail, W., D. Mainprice (1998), An olivine fabric database: an overview of upper mantle fabrics and seismic anisotropy, Tectonophysics, 296, 145. doi:10.1016/S0040-1951(98)00141-3
Dziewonski, A. M., D. L. Anderson (1981), Preliminary Reference Earth Model (PREM), Phys. Earth Planet. Inter., 25, 297. doi:10.1016/0031-9201(81)90046-7 Foster, I., C. Kesselman (1999), The Grid: Blueprint for a New Computing Infrastructure, Morgan Kaufmann Publishers, ISBN 1-55860-475-8. Girardin, N., V. Farra (1998), Azimuthal anisotropy in the upper mantle from observations of P-to-S converted phases: application to the southern Australia, Geophys. J. Int., 133, 615. doi:10.1046/j.1365-246X.1998.00525.x Kosarev, G. L., L. I. Makeeva, E. F. Savarenskiy,
E. M. Chesnokov (1979), Anisotropy affection under a seis-mostation at the volume waves, Earth Physics, 26. Kosarev, G. L., L. I. Makeeva, L. P. Vinnik (1984), Anisotropy of the mantle inferred from observations of P to S converted waves, Geophys. J. R. Astr. Soc., 76, 209. Menke, W., V. Levin (2003), The cross-convolution method for interpreting SKS spliting observations, with application to one and two-layer anisotropic earth models, Geophys. J. Int., 154, 379. doi:10.1046/j.1365-246X.2003.01937.x Nicolas, A., N. I. Christensen (1987), Formation of anisotropy in upper mantle peridotites — a review, Composition, Structure and Dynamics of the Lithosphere Asthenosphere System, Fuchs, K., Froideveaux, C. (Eds.), 111, Am. Geophys. Union, Geodyn. Ser., 16, Wash.. Renard, P., V. Badoux, M. Petitdidier, R. Cossu (2009), Grid Computing for Earth Science, Eos, 90(14), 117. doi:10.1029/ 2009E0140002
Sakawa, M. (1993), Fuzzy sets and interactive multiobjective
optimization, Plentum Press, NewYork. Savage, M. K. (1999), Seismic anisotropy and mantle deformation: what have we learned from shear wave splitting?, Rew. Geophys., 37(1), 65. doi:10.1029/98RG02075
Turban, E., J. Meredith (1994), Fundamentals Management Science, McGraw-Hill, Boston, USA.
Vinnik, L. P., G. L. Kosarev, L. I. Makeeva (1984), Litosphere anisotropy by SKS and SKKS observations, DAN USSR, 278, 1335.
Vinnik, L. P., D. Peregoudov, L. I. Makeeva, S. Oreshin (2002), Towards 3D fabric in the continental lithosphere and asthenosphere: the Tien Shan, Geoph. Res. Lett., 29, 1795. doi:10.1029/2001GL014588
Vinnik, L. P., I. M. Aleshin, S. G. Kiselev, G. L. Kosarev,
L. I. Makeeva (2007), Depth localized azimuthal anisotropy from SKS and P receiver functions: The Tien Shan, Geophys. J. Int., 169, 1289. doi:10.1111/j.1365-246X.2007.03394.x
I. M. Aleshin, K. I. Kholodkov, V. N. Koryagin, D. V. Pere-goudov, Institute of Physics of the Earth RAS, Moscow, Russia (ima@ifz.ru)
D. P. Medvedev, D. Yu. Mishin, M. N. Zhizhin, Geophysical Center RAS, Moscow, Russia (jjn@wdcb.ru)
