A generic formulation for microtremor exploration methods using three-component records from a circular array

We present a generic formulation for analysis methods that estimate phase velocities of Rayleigh and Love waves, by way of intermediary quantities called ‘spectral ratios’, using three-component records of microtremors from a circular array of sensors. At each time in-stant, the set of records are expanded in a Fourier series with respect to azimuth, so that we obtain a set of Fourier coefﬁcients that are represented in the form of complex time histories. We then estimate power- and cross-spectral densities of those Fourier coefﬁcients. The spectral densities, thus obtained, generally contain information on the phase velocities, powers and arrival directions of individual modes of Rayleigh and Love waves. By taking the quotient of two different sorts of such spectral densities, we can cancel out information on their powers and arrival directions, and extract information on their phase velocities alone. The spectral ratios have to be estimated, in practice, on the basis of records from a ﬁnite number of seismic sensors that are either evenly or unevenly spaced around a circumference. We describe a general procedure for their estimation, and discuss the effects of directional aliasing that the ﬁnite number of sensors and their conﬁguration have on the estimates of spectral ratios. We also discuss biases in the estimates of spectral ratios caused by the presence of incoherent noise. By using our method, it is also possible to estimate the central arrival direction of the microtremors, the ellipticity of the Rayleigh waves, and whether the Rayleigh waves are prograde or retrograde. full view of the entire theoretical perspective derived from Aki’s (1957) maximal use of the wonderful potential inherent in his algorithm. In the present paper, we further extend Henstridge’s (1979) generalization of Aki’s (1957) method and build up a consistent theory that integrates the whole family of relevant approaches. The aim of our study lies in presenting, in the most general formulation possible, the method to estimate phase velocities of Rayleigh and Love using three-component of microtremors from a circular array of sensors.

recast Aki's (1957) lengthy formulation in a simpler form using a spectral representation of random wavefields, and he also devised a more general theory that included Aki's (1957) method as a special case, except that it was applicable to the vertical motion alone. Although Henstridge's (1979) theory is founded on the strong assumption that a single mode (the fundamental mode) dominates the wavefield, Cho et al. (2004) have recently been able to lift this assumption, and they modified and reformulated Henstridge's (1979) theory to allow for cases where more than one wave modes are present. Aki's (1965) efforts for publicity notwithstanding, Aki's (1957) theory was slow to see its practical adaptation to microtremor exploration. Starting with the pioneering work by Okada & Sakajiri (1983), field implementation of Aki's method has been vigorously investigated in Japan, but most of the earliest reports were published in Japanese. It was first shown by Okada & Matsushima (1989) how one can deal with cases where both Rayleigh and Love waves are present in the horizontal motion. However, most of those studies were based on the assumption that a single mode (the fundamental mode) dominates the field of surface waves, and so have not fully exploited all the latent possibilities inherent in Aki's (1957) theory.
Different research articles have cited Aki's (1957) theory in slightly different ways with partial modifications, and each of them has its own advantages and disadvantages. This makes it difficult to take a full view of the entire theoretical perspective derived from Aki's (1957) approach and to make maximal use of the wonderful potential inherent in his algorithm. In the present paper, we further extend Henstridge's (1979) generalization of Aki's (1957) method and build up a consistent theory that integrates the whole family of relevant approaches. The aim of our study lies in presenting, in the most general formulation possible, the method to estimate phase velocities of Rayleigh and Love waves using three-component records of microtremors from a circular array of sensors.
We first represent the microtremor wavefield as a random field of surface waves that is stationary in time and space. The basic equations are derived in such a form that allows for the coexistence of Rayleigh and Love waves in the horizontal motion (e.g. Okada & Matsushima 1989) and also allows for the coexistence of multiple wave modes (e.g. Cho et al. 2004). In the first half of the present paper, we first restate Aki's (1957) original theory, and then enlarge and generalize it by introducing Henstridge's (1979) ideas, before finally illustrating some specific ways it can be adapted to microtremor exploration. It is noteworthy that our theory of analysis is formulated in such a way that the microtremor records, obtained with sensors placed around the circumference of a circle, are first integrated into intermediary quantities called 'spectral ratios' before information on the phase velocities is extracted from them. Our method has the potential to be able to analyse, with high accuracy and stability, very long wavelengths relative to the array radius in some cases, because there is no need to resolve individual plane wave components (Cho et al. 2004).
Our theory is built on the basic assumption that seismograms are available everywhere around a circumference of a certain radius, and also at its centre if necessary. In practice, however, seismograms are available only at a finite number of discrete locations where sensors are placed. It is most desirable to place sensors at even intervals around the circumference if at all possible, but it is likely that the constraint of the circular shape obliges one, in many practical situations, to place them at uneven intervals. In the latter half of our paper, we propose a general method to estimate Fourier coefficients of different azimuthal orders from seismograms obtained with sensors that are either evenly or unevenly spaced around a circumference, and we theoretically evaluate the effects of the finite number of sensors and their configuration on the estimates of spectral ratios. We also present a method to theoretically evaluate biases in the estimates of spectral ratios that are caused by the presence of incoherent noise specific to each component record from each sensor.

M I C RO T R E M O R S A S A S TAT I O N A RY R A N D O M F I E L D
We define the Fourier series expansion of a regular periodic function f (φ) of period 2π as where i is the imaginary unit. We also define the Fourier transform F of f (t) and the inverse transform F −1 as We denote the ensemble mean of a random function f (·) as where f n (·) indicates the nth sample of f (·). Suppose that we measure the field of microtremors with a circular array of sensors with radius r, and assume that three-component seismograms are available at any arbitrary location (r, θ ) on the circumference and at the origin O of coordinates which is the centre of the circle. As shown in Fig. 1, we label the vertical, radial and transverse axes as Z, R and T respectively, with the upward, centrifugal and counter-clockwise directions taken to be positive.
Elastic wave propagation theory predicts that, when vibration sources of microtremors lie at large enough distances and near the ground surface, the field of microtremors is dominated by surface waves, namely Rayleigh waves and/or Love waves that arrive at the seismic array in the form of plane waves. In the present section, we derive basic equations representing the wavefield in such situations.
We first consider, as incident Rayleigh and Love waves, single, deterministic, harmonic plane wave components which arrive at the seismic array from azimuth φ with angular frequency ω and absolute wavenumber k. If we denote by ζ R (ω, k, φ) and ζ L (ω, k, φ) the complex Fourier spectra that describe both their amplitude and phase at the origin of coordinates, the contribution of the harmonic waves to the Z-, Rand T-component seismograms at location (r, θ ) is obtained by correcting these spectra for the phase delay and for the rotation of the axes: where the superscripts R and L indicate that the quantity in question pertains to the Rayleigh and Love waves, respectively. The complex function h(ω, k), describing the amplitude ratio (reciprocal of the ellipticity) and phase delay between the vertical and horizontal components of the Rayleigh waves, takes purely imaginary values, because the phase delay is always either ±π /2. The ground motion is retrograde and prograde when Arg[h(ω, k)] = −π/2 and π/2, respectively. By triply integrating the above harmonic wave components over all frequencies, wavenumbers and arrival directions, we obtain the following expressions for each component seismogram: Full Rand T-component seismograms are obtained by summing up the contributions of Rayleigh and Love waves: We shall next consider microtremors, or a field of endless signals generated by a large number of unspecified vibration sources. The basic concept remains the same, except that the formalities have to be modified into more appropriate ones that are based on the theory of random processes; see textbooks by Yaglom (1962), Koopmans (1974) and Priestley (1981) for details. To be more specific, we have only to replace the part ζ R (ω, k, φ) dωk dk dφ in eq. (11) with a set function ζ R (dω, dk, dφ) which is called an integrated spectrum or a random spectral measure (e.g. Yaglom 1962). This gives the general representation for the vertical component {Z} of the random field of microtremors: The same argument holds for the horizontal components {R l } and {T l } (l = R or L): We regard microtremors as a random field stationary in time and space. This is equivalent to assuming the following orthogonality relation for the integrated spectra ζ R (dω, dk, dφ) and ζ L (dω, dk, dφ) with respect to frequency ω, wavenumber k, and arrival direction φ: where an asterisk denotes the complex conjugate and δ(·) stands for the Dirac delta function. The functions F R (ω, k, φ) and F L (ω, k, φ), which we call the frequency-wavenumber-direction (FWD) spectral densities, represent the intensity of plane wave components of Rayleigh and Love waves, respectively, arriving from direction φ with frequency ω and wavenumber k. The orthogonality represented by δ(ω − ω ) in eq. (23) corresponds to the stationarity in time, while the orthogonality represented by δ(k − k )δ(φ − φ ) corresponds to the stationarity in space.
If the power of Rayleigh and Love waves are concentrated on their discrete modes, the FWD spectral densities can be expressed in the following way, since the wavenumber then becomes a multivalued function of frequency: where f R(q) (ω, φ) and f L(q) (ω, φ) are frequency-direction spectral densities representing the intensity of the (q − 1)th mode of Rayleigh and Love waves, respectively, arriving from direction φ with frequency ω. We assume the Rayleigh and Love wave signals to be mutually uncorrelated, under the presumption that their vibration sources are mutually uncorrelated in the statistical sense: We have heretofore set up the following assumptions regarding the field of microtremors: (1) The microtremor field is represented as an ensemble of plane waves arriving from different directions with different intensities and different phase velocities, and is dominated by Rayleigh and Love waves (eqs 18-22).
(2) The field of Rayleigh waves and that of Love waves are each represented as a random field that is stationary in time and space (eq. 23). This corresponds to the assumption that waves arriving from different directions with different phase velocities are mutually uncorrelated.
(3) The wavenumber (or phase velocity) of Rayleigh waves and that of Love waves are each a multivalued function of frequency (eq. 24).

A K I ' S M E T H O D
The originality of Aki's (1957) approach lies in the idea of integrating the whole information on a given wavefield, represented by eqs (16)-(22), into a single quantity which is called the azimuthal average of the spatial autocorrelation function. We would like to start our argument by reformulating Aki's (1957) original theory in as plain a form as possible. We define the spatial autocorrelation function between the vertical-motion record obtained at a point (r, θ ) on the circumference and that obtained at its centre O as Substituting eqs (18), (23) and (24) into eq. (26) yields where the function h(ω, k R(q) (ω)) has been rewritten as h (q) (ω) for the sake of simplicity. We define the azimuthally averaged spatial autocorrelation coefficientρ Z 0 (ω, r ) as the azimuthal averageρ Z (ω, r ) of the spatial autocorrelation function ρ Z (ω, r, θ ) normalized by the valueρ Z (ω, 0) it takes at the centre. This denominator, incidentally, equals the power-spectral density of vertical motion at the centre: By substituting eqs (27) and (29) into eq. (28) and using the relationship (A1), we havē where J m (·) denotes the Bessel function of the first kind of order m, and is the power partition ratio of the (q − 1)th mode to the total power of the Rayleigh waves in vertical motion, and satisfies the following relationship: It is possible to estimate the value of the argument rk R(q) (ω) by way of eq. (30), because the azimuthally averaged spatial autocorrelation coefficientρ Z 0 (ω, r ) can be estimated from circular array measurement records. Of special interest here is the case where only the fundamental mode of Rayleigh waves dominates (N R = 1). In this case, the value of rk R (ω) (we henceforth omit the superscript (1), representing quantities pertaining to the fundamental mode, wherever we assume that it dominates) can be estimated by simply inverting the function J 0 (·) using the observed value ofρ Z 0 (ω, r ). One can subsequently estimate the wavenumber k R (ω) since the array radius r is known, and finally the phase velocity by c R (ω) = ω/k R (ω).
In the general case, where the dominance of a single mode cannot be assumed, we have two unknowns for each mode, namely k R(q) (ω) and α (q) (ω)(q = 1, . . . , N R ). Aki (1957) proposed deploying arrays of (2N R − 1) different radii, setting up the observation eq. (30) for each radius and solving them simultaneously with the compatibility condition (32).
Similarly, let us define the spatial autocorrelation functions for the horizontal motion: By substituting eqs (16) and (17) into eqs (33) and (34) and making use of eq. (25), we have Substituting eqs (19)-(22) into the above equations and using eqs (23) and (24), we obtain, after some algebra, We define the azimuthal averagesρ R (ω, r ) andρ T (ω, r ) of the spatial autocorrelation functions, and the azimuthally averaged spatial autocorrelation coefficientsρ R0 (ω, r ) andρ T 0 (ω, r ), in the same way as we did so for the vertical motion: Substituting eqs (37), (38), (40) and (42) into eqs (39) and (41) and using eqs (A6) and (A7), we obtain where is the power partition ratio of the (q − 1)th mode of Rayleigh or Love waves to the total power of both Rayleigh and Love waves in horizontal motion, and satisfies the following condition: Eqs (43) and (44) are the observation equations for the horizontal motion, which correspond to eq. (30) in the case of the vertical motion. As we have stated in Section 1, Aki (1957) accounted only for special situations where either Rayleigh waves alone (N L = 0) or Love waves alone (N R = 0) are present, and this constituted a weakness of his method. To overcome this, Okada & Matsushima (1989), Ferrazzini et al. (1991), Chouet et al. (1998) and Saccorotti et al. (2003), under the assumption that a single mode dominates (N R = N L = 1), first estimated k R (ω) from the observed value ofρ Z 0 (ω, r ) by inverting the observation eq. (30), and then solved the two observation eqs (43) and (44) simultaneously to search for a set of k L (ω) and γ L (ω) that gave the best fit to the observed values ofρ R0 (ω, r ) andρ T 0 (ω, r ) (note the compatibility condition γ R (ω) + γ L (ω) = 1). Métaxian et al. (1997) solved the two eqs (43) and (44) for the three unknowns k R (ω), k L (ω) and γ L (ω), this being made possible by their constraining assumption that k R (ω) and k L (ω) are functions of ω of a particular form.
In the general case where the dominance of a single mode cannot be assumed, we have three unknowns, namely k R(q) (ω), α (q) (ω) and γ R(q) (ω) (q = 1, . . . , N R ), for each mode of Rayleigh waves and two unknowns, namely k L(q) (ω) and γ L(q) (ω) (q = 1, . . . , N L ), for each mode of Love waves. The number of independent unknowns is 3N R + 2N L − 2, since the above unknowns are constrained by the two compatibility conditions (32) and (46), and so they have to be solved for with at least as many observation equations. Note, however, that there have to be at least (N R − 1) observation equations for vertical motion because α (q) (ω) can be determined from vertical-motion records alone; and that there have to be at least (N R + 2N L − 1) observation equations for horizontal motion because k L(q) (ω), γ R(q) and γ L(q) can be determined from horizontal-motion records alone.
On the other hand, if we have at least (2N R + 2N L − 1) observation equations for horizontal motion, we can determine k R(q) (ω), k L(q) (ω), γ R(q) and γ L(q) , even though we cannot determine α (q) (ω). Similarly, if we have at least (2N R − 1) observation equations for vertical motion, we can determine α (q) (ω) and k R(q) (ω), even though we cannot determine k L(q) (ω), γ R(q) and γ L(q) , a situation that we have already described above.

G E N E R A L I Z I N G A K I ' S M E T H O D
In this section we enlarge and generalize, by introducing Henstridge's (1979) ideas, Aki's (1957) theory that we have described in the foregoing section. Henstridge (1979) reinterpreted the azimuthal average of the spatial autocorrelation function, appearing in Aki's (1957) formulation, as a quantity related to the zeroth-order coefficient in the Fourier series expansion of the microtremor field with respect to azimuth. Henstridge's (1979) ingenuity lies in the discovery of the utility of information that is contained in the Fourier coefficients of higher orders.
Let us expand the time histories Z (t, r, θ ), R(t, r, θ) and T(t, r, θ ) in Fourier series with respect to θ : where X stands for either Z, R or T. Substituting into eq. (48) either eq. (18), eqs (16), (19) and (21), or eqs (17), (20) and (22), and making use of the formulae (A1), (A4) and (A5), we obtain the following expressions for Z m (t, r), R m (t, r) and T m (t, r): where ζ l m (dω, dk) is the mth order coefficient in the Fourier series expansion of the integrated spectrum ζ l (dω, dk, dφ)(l = R or L) with respect to azimuth φ: The above Fourier coefficients are related to the Fourier coefficients for the FWD spectral densities by the relationship where δ ll denotes the Kronecker delta. If the power of the Rayleigh and Love waves are concentrated on their discrete modes, the Fourier coefficients for the FWD spectral densities can be expressed, in the same way as in eq. (24), by The meaning of the Fourier coefficients Z m (t, 0), R m (t, 0) and T m (t, 0) at the array centre shall be explained in a later section.

Spectral densities
The cross-or power-spectral density G XmYn (r 1 , r 2 ; ω) of the Fourier coefficients X m (t, r 1 ) and Y n (t, r 2 ), represented by eqs (49)-(51) (X and Y stand for either Z, R or T), is defined as: where Wiener-Khintchine's theorem states that the spectral density G XmYn (r 1 , r 2 ; ω) is represented as the Fourier transform of the cross-or autocorrelation function: Substituting eqs (49)-(51) into eq. (57) and making use of eqs (53) and (54), we obtain concrete expressions for different kinds of spectral densities in the following way using FWD spectral densities: Henstridge (1979) derived a representation for the spectral density of vertical motion, corresponding to G ZmZn (r 1 , r 2 ; ω) in our notation, under the limitative assumption that a single mode dominates. Our equations shown above are an enlargement of his theory to account for cases where multiple modes are present, and are also a generalization that includes the case of horizontal motion.

Spectral ratios
The spectral densities represented by eqs (58)-(60) contain, via f l(q) m−n (ω), information on the intensities and arrival directions of different wave modes. By taking the quotient of two appropriately selected spectral densities, we can cancel that information out and derive different formulae that are useful for the analysis of phase velocities, e.g.
Incidentally, the denominator G Z 0Z 0 (0, 0; ω) appearing in some of the above equations equals 4π 2 times the power-spectral density of vertical motion at the centre of the array, while the denominator G T 1T 1 (0, 0; ω) (= G R1R1 (0, 0; ω)) equals π 2 times the power-spectral density of horizontal motion at the centre of the array. This means that the horizontal-to-vertical (H/V) spectrum of microtremors, which has often been used for the purpose of evaluating amplification characteristics of shallow subsurface structures (e.g. Nakamura 1989; Konno & Ohmachi 1998;Arai & Tokimatsu 2004), is represented in the following way using our notation: As we have mentioned above, we have three unknowns for each mode of Rayleigh waves, namely k R(q) (ω), α (q) (ω) and γ R(q) (ω)(q = 1, . . . , N R ), and two unknowns for each mode of Love waves, namely k L(q) (ω) and γ L(q) (ω)(q = 1, . . . , N L ). However, we do not necessarily have to estimate all these unknowns, since our main interest in phase velocity analysis lies in estimating k R(q) (ω) and k L(q) (ω). All we have to do is to pick up necessary ones out of eqs (62)-(72) and to solve them simultaneously with the compatibility conditions (32) and (46) concerning the power partition ratios α (q) (ω) and γ l(q) (ω). If we are to estimate the phase velocities of Rayleigh waves alone, it is sufficient to pick up one or more out of the spectral ratios (62)-(64) and to solve them simultaneously with the compatibility condition (32).
As Tokimatsu et al. (1992c), Tokimatsu (1997) pointed out, it is in many cases safe to assume that the fundamental mode dominates the field of surface waves when, in the medium in question, the seismic wave velocities increase monotonically with depth. A number of data analysis results (e.g. Cho et al. 2004) suggest that this assumption of the dominance of a single mode is not far off the reality. In such cases, the procedure to estimate phase velocities becomes very simple. Any one out of eqs (62)-(64), (69) and (70) is sufficient for the estimation of the phase velocities of Rayleigh waves, whereas either one of eqs (71) and (72) is sufficient for the estimation of the phase velocities of Love waves. In fact, by setting Also, by setting N L = 1 we have Once the values of the spectral ratios on the left-hand side are known from measurements, the wavenumbers k R (ω) and k L (ω) can be estimated by inverting the functions on the right-hand side, since the array radii r, r 1 and r 2 are known. One finally obtains the phase velocities by c R (ω) = ω/k R (ω) and c L (ω) = ω/k L (ω). The functions appearing on the right-hand side of eqs (62)-(72) and (74) The method we have described in the present paper is related to existing theories in the following manner: (1) Eq. (62) is equivalent to eq. (30) which was derived by Aki (1957). The quantityρ Z 0 (ω, r ) on the left-hand side of eq. (30) is obtained by first taking the cross-spectral density ρ Z (ω, r , θ ) between the record Z(t, r, θ ) on the circumference and the record Z(t, 0, θ ) at the centre, integrating it with respect to azimuth θ , and finally normalizing it by the power-spectral density at the centre. On the other hand, the quotient G Z 0Z 0 (0, r ; ω)/G Z 0Z 0 (0, 0; ω) on the left-hand side of eq. (62) is obtained by first integrating the record Z(t, r, θ ) on the circumference with respect to azimuth θ, taking the cross-spectral density G Z 0Z 0 (0, r ; ω) between this and the record Z(t, 0, θ) at the centre, and finally normalizing it by the power-spectral density at the centre. The two quantities are in fact identical to each other, although, in their derivation, the acts of integrating with regard to azimuth and of taking the cross-spectral density come in opposite orders. Similarly, eqs (65) and (66) are equivalent, respectively, to eqs (43) and (44) which were derived by Aki (1957). See Appendix B for details.
(2) Aki (1965) proposed using records from around a semicircle of radius r, instead of records from around a full circle, in averaging, with respect to azimuth θ , the cross-spectral density ρ Z (ω, r, θ ) between the record at centre and the records on the circumference. With our theoretical formulation, however, we cannot substitute data from around a circle with data from around a semicircle: see Appendix C for details.
(3) By assuming the dominance of a single mode and setting N R = 1 in eq. (63), we obtain a quantity which Henstridge (1979) called the 'scatter function.' (4) Eq. (64) is the same spectral ratio that Cho et al. (2004) defined and used in their analysis of phase velocities.
In practice, spectral ratios pertaining to a field of microtremors, which is a stationary random process, have to be estimated on the basis of a finite number of sample records, and so their estimates are inherently subject to uncertainties. This translates to uncertainties in the estimates of the phase velocities of surface waves when we invert eqs (74) and (75) under the assumption that a single mode dominates. Cho et al. (2004) pointed out that the variance in the estimate of the spectral ratio represented by eq. (64) is proportional to the magnitude squared of the spectral ratio itself: On the basis of this idea, they theoretically evaluated, under the assumption that a single mode dominates, the variance in the phase velocity estimated by the third equation in (74). Exactly parallel arguments are possible on the variance pertaining to any other spectral ratio which is represented as a quotient of two power-spectral densities (eqs 63, 67, 68, 70 and 72).
In the following, we discuss methods to estimate the values of different spectral ratios using field records from a finite number of sensors around a circumference, and we also present a method to theoretically evaluate biases contained in the estimates thus obtained.

P RO C E S S I N G D ATA F RO M A F I N I T E N U M B E R O F S E N S O R S
Let us denote the three-component seismograms, obtained with a sensor at location (r, θ ), in a vector form as where u 3 (t, r, θ ), u 1 (t, r, θ ) and u 2 (t, r, θ) represent the vertical, east-west and north-south components, with the upward, eastward and northward motion taken to be positive respectively (Fig. 1). The Z-, Rand T-component seismograms at location (r, θ) ( Fig. 1) are obtained by the following coordinate transformation: The parameter θ becomes a dummy when r = 0, so in that case we use the notation u(t, 0) = (u 3 (t, 0), u 1 (t, 0), u 2 (t, 0)) T . The Fourier coefficients for the record at the centre is obtained in the following manner by setting r = 0 in eq. (78) and substituting it into eq. (48): Suppose that we are to estimate the Fourier coefficients (eq. 48) from records of N sensors placed at locations (r , θ j ) ( j = 1, . . . , N ) around a circumference of radius r (>0). First, we approximate the infinite sum (eq. 47), appearing in the Fourier series expansion of X (t, r, θ) with respect to θ , by a finite sum of terms of orders up to ± K: By setting θ = θ 1 , . . . , θ N in eq. (82), we obtain a system of simultaneous equations concerning X m (t, r): where The relation N ≥ 2K + 1 must hold so that the system (83) of equations may not be underdetermined. Let us minimize, in the least squares sense, the approximation errors in eq. (83). We define the squared norm L 2 as where the superscript H indicates the Hermitian transpose (complex conjugate transpose). In order to evaluate the estimateĉ of c that minimizes L 2 , we partially differentiate L 2 with respect to the real and imaginary parts of the jth element of c, respectively, and equate both classes of derivatives to zero: By combining the above equations, we obtain namely, This is the well-known normal equation. If Φ H Φ has an inverse matrix, we obtain the estimatesX m (t, r ) for the Fourier coefficients of different orders bŷ or, if we denote the pq-th component of matrix ( When Φ is a square matrix (N = 2K + 1) which has an inverse matrix, (Φ H Φ) −1 Φ H = Φ −1 , and eq. (91) reduces to a simple inverse-matrix operationĉ = Φ −1 X. When N(≥ 2K + 1) sensors are placed equidistantly around a circumference, we can put, without loss of generality, In this case, it can be shown that so that we obtain, by substituting this into eq. (91), In the above formulation, we have defined vectorĉ as a vertical row of estimatesX m (t, r ) (m = −K , . . . , K ) for the (2K + 1) complex Fourier coefficients. Among them, however, only (K + 1) are mutually independent, because they are constrained by the conjugate relationŝ X m (t, r ) =X * −m (t, r ). Alternatively, it is also possible to defineĉ in the following way as a real vector of dimension (2K + 1), by regarding the real and imaginary parts of the Fourier coefficient of each order as separate unknowns (e.g. Cho et al. 2004). In this case, all (2K + 1) elements of vectorĉ are mutually independent.
When we use the above-mentioned array of equidistantly spaced sensors, the real and imaginary parts of the Fourier coefficient of each order are estimated by the following equations derived by transforming eq. (94): In general, the operation to transform the three-component microtremor records from M sensors, placed around a circle or more than one circles of different radii and also at their centre if necessary, into Z-, Rand T-components via the coordinate transform defined by eq. (78), can be written, formally, in the following way in terms of an appropriate transform matrix B:

V(t) ≡ (Z(t), R(t), T(t)) T ,
where u 3 (t), u 1 (t) and u 2 (t) are vectors of dimension M in which the vertical, east-west and north-south records from each sensor are aligned in a vertical row, respectively, with the jth element corresponding to the record from the jth sensor. Similarly, Z(t), R(t) and T(t) are vectors of dimension M in which the Z-, Rand T-component records at each sensor location, respectively, are aligned in a vertical row. Since no coordinate transform can be defined at the centre r = 0, we simply rename u 3 (t, 0), u 1 (t, 0) and u 2 (t, 0) to Z(t, 0), R(t, 0) and T(t, 0), respectively, with no modification. The procedure to estimate the Fourier coefficients for the Z-, Rand T-components by eq. (91) using the Z-, Rand T-component time-history seismograms at each sensor location around a circle or more than one circles of different radii (including r = 0) can be written, formally, in the following way in terms of an appropriate transform matrix A: whereĉ Z (t),ĉ R (t) andĉ T (t) are vectors in which the estimatesX m (t, r ) (X = Z , R, T ) for the Fourier coefficients for the Z-, Rand T-components, respectively, obtained for one or more radii (including r = 0), are aligned in a vertical row. The estimates of the cross-and power-spectral densities between all combinations of Fourier coefficients, can be written in a collective form as a cross-spectral density matrix with respect to vectorĈ: In the present article, we have represented the cross-spectral densitiesĜ XmY n (r 1 , r 2 ; ω) as the Fourier transform of cross-correlation functions between two time-series Fourier coefficients, but in practice, the use of the direct method (Bendat & Piersol 1971), which Fourier transforms the original time-series Fourier coefficients by FFT, is preferable from the viewpoint of computational efficiency. As long as the microtremor signals are stationary (homogeneous) in time and space, the cross-spectral densities have the following mathematical properties, even when the signals do not arrive homogeneously from all directions, which could be exploited to increase the flexibility of array design for measurements: (1) We do not necessarily have to make measurements simultaneously at all points around the circle to estimate the cross-spectral densities G XmYn (r, 0; ω) between the Fourier coefficients for the record at the centre and for the records on the circumference. For instance, G XY may be estimated with the following algorithm, derived from eq. (91), by repeating, with changing θ j , the procedure to get records simultaneously at a point A j (r , θ j ) on the circumference and at its centre O and then to estimate the cross-spectral densities F X (s, r, θ j )Y * n (s − t, 0) between them (Fig. 3a): (2) In this case, the set of records from point A j on the circumference and from its centre O can be substituted by a set of records from two other locations, translated from their original locations without changing their relative geometry (Fig. 3b).
(3) The above consideration leads us to note that, if we place an even number 2N of sensors equidistantly around a circumference, the relative geometry between point A j and the centre O is merely a translation of the relative geometry between the centre O and point A j+N . From the viewpoint of efficiency, it is therefore preferable to place an odd number of sensors around the circle if we are to deploy sensors at even intervals around a circle plus another at its centre.
(4) As long as the wavefield is stationary in time and space, we do not necessarily have to use spectral densities, estimated with records from the same time and locations, as the numerator and the denominator in calculating spectral ratios. For example, we may use the power-spectral density at a point on the circumference, instead of that at its centre, as the denominator when we estimate the spectral ratio G Z 0Z 0 (0, r ; ω)/G Z 0Z 0 (0, 0; ω).

B I A S E S I N T H E E S T I M AT E S O F S P E C T R A L D E N S I T I E S
In this section, we present a method to theoretically evaluate the biases in the estimates of the spectral densities defined in the foregoing section. We first decompose the record u (eq. 77) into the signal part u s that is modelled by eqs (16)-(25) and the noise part u n that is uncorrelated with the signal: u s (s, r, θ)u n H (s − t, r, θ) = 0.
Substituting eq. (105) into eq. (100) yieldŝ where the superscripts s and n on U(t) reflect those on their component vectors u(t, r, θ ). Substituting this into (103) and using the uncorrelatedness between signal and noise, we obtain This means that the estimatesĜ s (ω) of the spectral densities for the signal part of the records, on the one hand, and the biasesĜ n (ω) in the estimates of the spectral densities due to the presence of noise, on the other hand, can be evaluated independently of each other. Note that the estimatesĜ s (ω) contain biases coming from the finite number of sensors. We hereafter denote elements of the cross-spectral matricesĜ s (ω) andĜ n (ω) asĜ s XmY n (r 1 , r 2 ; ω) andĜ n XmY n (r 1 , r 2 ; ω), respectively.

Biases in the estimates of the spectral densities for the signal part: directional aliasing
The estimatesĜ s XmY n (r 1 , r 2 ; ω) of the cross-or power-spectral densities for the signal part of the records are given by the Fourier transform of the cross-or autocorrelation functions between the estimatesX s m (t, r ) of the Fourier coefficients of the signal part: G s XmY n (r 1 , r 2 ; ω) = F X s m (s, r 1 )Ŷ s * n (s − t, r 2 ) , when r > 0, andŶ s n (t, r ) is defined in the same manner. The estimatesX s m (t, r ) of the Fourier coefficients generally do not agree with their true values X s m (t, r) because of directional aliasing, due to the fact that the Fourier coefficients are estimated using records from discrete locations alone (e.g. Henstridge 1979). Note, however, that directional aliasing brings no bias to the estimates of the Fourier coefficients X s m (t, 0) at the centre.
In the case of an array of evenly spaced sensors, we have, by substituting eq. (92) into eq. (112) and referring to eq. (93), In other words, the estimateX s m (t, r ) of a Fourier coefficient, in this case, equals the sum of its true value X s m (t, r) and Fourier coefficients X s m+jN (t, r) of higher orders ( j = ±1, ±2, . . .). In general, increasing the number N of sensors suppresses the effects of directional aliasing. Eq. (112) corresponds to a generalization of eq. (113) which accounts for arrays of unevenly spaced sensors, and demonstrates that the effects of directional aliasing depends not only on the number N of sensors but also on their configuration.
The effects of directional aliasing also control the deviation ofĜ s XmY n (r 1 , r 2 ; ω), given by eq. (111), from their true values. In the case of an array of evenly spaced sensors, for instance, we can obtain expressions for the estimates of different types of spectral densities, including the effects of directional aliasing, by substituting eqs (49)-(51) into eq. (113), evaluating the estimates for the Fourier coefficients, and by substituting those results into eq. (111).
Carrying this procedure out withĜ s Zm Zn (r 1 , r 2 ; ω), for example, we have, by virtue of the relationship (53), All terms in eqs (114) and (115) with suffixes other than j = l = 0 represent errors due to directional aliasing. The fact that all directional aliasing terms include, as multipliers, the Fourier coefficients f R(q) m (ω) of the frequency-direction spectral densities for each wave mode, indicates that the effects of directional aliasing depend on the directionality of the wavefield. The appearance of high-order Bessel functions in the multipliers indicates that the effects of directional aliasing tend to be larger for larger values of rk, namely in smaller wavelength ranges.

Biases in the estimates of spectral densities due to the presence of noise
We denote byX n m (t, r ) (X = Z , R, T ) the biases in the estimates of Fourier coefficients, brought about by the presence of noise u n (t, r, θ ) which is a stationary random process uncorrelated with the signal. The corresponding biases in the estimates of the spectral densities due to the presence of noise are written, by virtue of eq. (110), aŝ G n XmY n (r 1 , r 2 ; ω) = F X n m (s, r 1 )Ŷ n * n (s − t, r 2 ) .
To exemplify some ways to utilize the above equation, let us assume, for simplicity, that noise appearing in records of different sensors, or in different component records from an identical sensor, are mutually uncorrelated. We can think of stationary electric noise, generated inside a recording system, as an example of real noise that falls into this category. The cross-spectral densities of such incoherent noise components satisfy the following orthogonality relation: where f n l (r , θ ; ω) is the power-spectral density of the noise that contaminates the u l -component record of the sensor at location (r, θ ). In the present article, we do not consider noise of the sort that shows correlation between different components or different sensors.
If, in addition, the power-spectral densities f n l (r , θ; ω) of the noise are identical for all sensors and if the noise has identical intensities for the two horizontal components ( f n 1 (r , θ; ω) = f n 2 (r , θ ; ω)), we can reduce eq. (122), by virtue of eq. (121), to: This means that, out of all estimates of the spectral densities from records of an array of evenly spaced sensors, only the power-spectral densities of Fourier coefficients of identical component records and of identical orders are subject to the effects of noise under the above-mentioned conditions.

Examples
As an example to illustrate how to evaluate biases in the estimates of spectral densities coming from directional aliasing and the presence of incoherent noise, let us assume that we are going to estimate the spectral ratio represented by the left-hand side of eq. (74) from verticalcomponent records of microtremors obtained with an array of three evenly spaced sensors. We postulate, for simplicity, that only the fundamental-mode Rayleigh waves dominate the vertical-component field of microtremors, and that waves arrive with equal intensities from all directions. This means that the Fourier coefficient f R m (ω) of the frequency-direction spectral density of the waves is non-zero if and only if m = 0. We also assume that noise in the records have identical power-spectral densities for all sensors.
By setting N R = 1 and N = 3 in eq. (114) and using f R m (ω) = 0 (m = 0), we obtain the following expressions for the estimates of different spectral densities with the effects of directional aliasing included: The spectral ratios represented by eq. (74) can be estimated by taking mutual quotients of the above quantities. Fig. 4 shows their estimates, thus calculated, along with their true values (corresponding to the j = 0 term in eq. 124) that contains no effect of directional aliasing. The directional aliasing effects manifest themselves in different ways for different types of spectral ratios, but we can see from the figure that the effects tend to be small in small rk ranges, as is so expected from theoretical considerations. As is expected from eq. (124), directional aliasing has no effect on the estimate of the spectral ratio G Z 0Z 0 (0, r ; ω)/G Z 0Z 0 (0, 0; ω) over all frequency ranges. In fact, directional aliasing has no effect, even when N = 1, on this spectral ratio which Aki (1957) used in his microtremor exploration, as long as the waves arrive isotropically. This is exactly what Aki (1957) meant when he stated that there is no need to carry out integration with respect to azimuth when waves arrive isotropically.   The biases in the estimates of spectral densities due to the presence of incoherent noise are given, by virtue of eqs (118), (119) and (123), By using eqs (124) and (125), we obtain expressions for the estimates of spectral ratios where the effects of directional aliasing and of the presence of noise are both taken into account:   Figure 5. Same as Fig. 4 except that the effects of incoherent noise have been taken into account. The signal-to-noise ratio is set at either 10 or 100. Fig. 5 plots eq. (126) as a function of rk, with the SN ratio set at 10 and 100, respectively. The presence of incoherent noise we have postulated in this example simply diminishes the spectral ratioĜ Z 0Z 0 (0, r ; ω)/Ĝ Z 0Z 0 (0, 0; ω) by a fixed ratio 1/(1 + ε). For the spectral ratioŝ G Z 0Z 0 (r, r ; ω)/Ĝ Z 0Z 0 (0, 0; ω) andĜ Z 1Z 1 (r, r ; ω)/Ĝ Z 0Z 0 (0, 0; ω), however, the presence of noise has the more complicated effect of adding a small term to their value diminished by a fixed ratio, and so the effects of noise tend to be large in frequency ranges where their true value approaches zero. For the spectral ratioĜ Z 0Z 0 (r, r ; ω)/Ĝ Z 1Z 1 (r, r ; ω), noise has the effect of adding small terms to both spectral densities in the numerator and the denominator, so the effects of noise grow large in frequency ranges where its true value either approaches zero or tends to infinity.

D I S C U S S I O N S A N D C O N C L U S I O N
We have presented a generic formulation for methods to analyse phase velocities of Rayleigh and Love waves, by way of intermediary quantities called 'spectral ratios,' using three-component records of microtremors from a circular array of sensors. Our theory, which is an enlargement of the approaches devised by Aki (1957) and Henstridge (1979), involves no procedure to resolve individual plane wave components contained in the field of microtremors, in contrast to the frequency-wavenumber spectral method that has been popularly used. Because of this property, our method is expected to demonstrate a higher stability of analysis than the frequency-wavenumber spectral method in long wavelength ranges. In fact, Cho et al. (2004) were able to stably analyse the phase velocities of surface waves in extremely long wavelength ranges, up to 30 times the array radius, by using our eq. (64).
Our theory is built on the assumption that microtremor records are available everywhere around the circumference, but in practice, the spectral ratios have to be estimated on the basis of records obtained with a finite number of sensors placed around the circumference, and the records are subject to the influence of noise. We have thus (i) presented a general method of data processing for estimating spectral ratios from seismograms obtained with a finite number of sensors that are spaced either evenly or unevenly around the circumference, (ii) theoretically evaluated the biases in the estimates of spectral ratios due to the finite number of sensors and their configuration, and also (iii) theoretically evaluated the biases in the estimates of spectral ratios due to the presence of incoherent noise. How the errors in the estimates of spectral ratios translate to errors in the estimates of the phase velocities of surface waves, and how one can theoretically deal with coherent noise that shows correlation between different sensors or different components, have not been discussed in the present article and thus are subjects for future investigations.
The present article has been intended as a complete theoretical description of our generic method of microtremor exploration, and so we did not demonstrate examples of implementation to field data. Until we show them in later publications to come, the reader is referred to Cho et al. (2004) for more information on such practicalities as data quality checks, data processing and the evaluation of analysis results obtained.
As we have stated earlier, it has been customary to assume the dominance of a single wave mode when adapting Aki's method to field data (e.g. Ferrazzini et al. 1991;Malagnini et al. 1993;Métaxian et al. 1997;Chouet et al. 1998;Morikawa et al. 2004). In the present article, however, we have demonstrated how the presence of multi mode Rayleigh and Love waves can be dealt with theoretically. This is not just for the sake of theoretical completeness but also out of our conviction that multi mode analysis is a promising subject for future researches. Although it is difficult to envisage, at the current stage, the whole variety of ways this could be implemented to real data analysis, Tokimatsu et al.'s (1992bTokimatsu et al.'s ( ,c, 1997 and Ohori et al.'s (2002) approach deserves attention for its broad potential applicability. They determined subsurface velocity structures by a method that is roughly tantamount to optimizing the fit between the right-hand side of eq. (62), theoretically expected when you consider the presence of multiple modes, and the left-hand side of eq. (62) that is estimated from microtremor records, although unlike ours, their method of analysis is 1-D in nature.
Though we have not mentioned it so far, it is also possible to derive, from our formulation, methods to estimate other attributes of microtremors, such as the central arrival direction of waves, the ellipticity of the Rayleigh waves and whether the Rayleigh waves are prograde or retrograde. This presents another interesting subject for future investigations, and the interested reader is referred to Appendices D and E for details.

A C K N O W L E D G M E N T S
The present work was partially supported by a Grant-in-Aid for Scientific Research (no. 16360281) of the Japan Society for the Promotion of Science (JSPS), and a Grant-in-Aid for Scientific Research (no. 17651101) of the Japanese Ministry of Education, Culture, Sports, Science and Technology (MEXT). so that At first sight, this equation might give the impression that Aki's (1957) spectral ratioρ Z 0 (ω, r ) is not equivalent to the spectral ratio we have defined by eq. (62). However, Z m (t, 0) is identically zero when m = 0 and, as can be verified from eq. (58), G ZmZm (r, 0; ω) is non-zero when and only when m = 0, so that the two spectral ratios are in fact mutually equivalent: For the R-component of horizontal motion, similarly, we have, from eq. (33), Averaging this over all azimuths yields so that Just as in the case of vertical motion, Aki's (1957) spectral ratioρ R0 (ω, r ) might at first sight seem inequivalent to the spectral ratio we have defined by eq. (65). Recalling that R m (t, 0) is identically zero when m = ±1, however, we find out that the two spectral ratios are in fact closely related to each other: We can see from eq. (59) that G R(−1)R(−1) (r , 0; ω) = G R1R1 (r , 0; ω), so that or that the two spectral ratios are in fact mutually equivalent. The same argument holds forρ T 0 (ω, r ): = G T 1T 1 (r, 0; ω)/G T 1T 1 (0, 0; ω) (B11)

A P P E N D I X C : P O S S I B I L I T Y O F U S I N G D ATA F RO M A RO U N D A S E M I C I RC L E I N TA K I N G T H E A Z I M U T H A L AV E R A G E
Aki (1965) suggested using records from around a semicircle of radius r, instead of a full circle, when you calculateρ Z (ω, r ) by eq. (29), or when you take the azimuthal average of the cross-spectral density between records at the centre and records on the circle. In fact, the symmetry relation can be derived from eq. (27), as well as from the spatio-temporal translatability that is derived from the stationarity of the microtremor field. By using this symmetry, eq. (29) is transformed as For the horizontal components, similarly, it follows from the symmetry that In our theoretical formulation, however, the records on the circle are integrated with respect to azimuth before the cross-spectral density between the records at the centre and the records on the circle is calculated. Unfortunately, the data from around a full circle cannot be substituted with data from around a semicircle by consideration of symmetry when the calculation is done in this order. It is still possible, however, that you happen to use records from sensors placed only around a semicircle when you attempt to estimate azimuthal integrals with records from a finite number of sensors placed at discrete locations around a full circle.

A P P E N D I X D : A R R I VA L D I R E C T I O N O F WAV E S
The method of microtremor exploration we have described in the present article is basically aimed at analysing phase velocities of surface waves, but it is also capable of providing information on the arrival direction of waves under situations where the wavefield can be regarded as dominated by a single mode. We assume, in this section, that N R = 1 for Rayleigh waves and that N L = 1 for Love waves. Henstridge (1979) pointed out that, under situations where waves arriving from a certain direction φ R 0 overpower those from all other directions and the FWD spectral density f R (ω, φ) in eq. (24) is given by it is possible to determine φ R 0 by circular array analysis, because the Fourier coefficients for the FWD spectral densities with all orders other than zero contain information on the arrival direction of waves. In fact, the complex argument of the mth order Fourier coefficient is given, in that case, by Arg f R m (ω) = −mφ R 0 (ω).
In the long-wavelength range of rk R (ω) ≤ 3.8 (wavelength λ(ω) ≥ 1.6r ), where J 1 (rk R (ω)) takes positive values, we have In fact, the same class of analysis is possible even when we relax Henstridge's (1979) conditions and assume that waves arrive from a finite range of directions. Suppose, for example, that the wave power, peaked at azimuth φ R 0 , linearly decreases until it falls to a constant value for all azimuths deviated more than φ R d (≤ π ) from the peak direction. The frequency-direction spectral density f R (ω, φ), defined by eq. (24), can then be written as When p 0 (ω) p(ω) and φ R d = π , this represents the situation where waves arrive from all directions with almost identical intensities, except that they arrive from azimuth φ R 0 with a slightly larger intensity. Even in this case, the complex argument of the Fourier coefficient f R m (ω), obtained by the Fourier series expansion of eq. (D5), is still represented by eq. (D2), and the formula (D4) remains valid.

A P P E N D I X E : E L L I P T I C I T Y O F T H E R AY L E I G H WAV E S
The amplitude ratio |h (q) (ω)| (reciprocal of the ellipticity) and phase delay (±π /2) between the vertical and horizontal components of Rayleigh waves can, in principle, also be estimated using our method of microtremor exploration. If we are able to estimate these quantities along with the phase velocities of Rayleigh and Love waves, they are expected to impose much greater constraints on the model of the underground velocity structure which they are functions of (e.g. Tokimatsu et al. 1992c;Tokimatsu 1997).
In terms of the H/V spectrum of microtremors (eq. 73), this can be written in the following way: We need the estimates of α (q) (ω) and γ R(q) (ω) if we are to evaluate the reciprocal ellipticity |h (q) (ω)| by way of eq. (E1) or (E2). It is, in principle, possible to regard α (q) (ω) and γ R(q) (ω) as unknowns and estimate them simultaneously with the phase velocities using microtremor records from a circular array. From a practical point of view, however, it is in many cases sufficient to assume appropriate values of α (q) (ω) and γ R(q) (ω) that are either calculated by medium response analysis or are empirically known to be reasonable, and evaluate |h (q) (ω)| using the value of the H/V spectrum obtained at a single point, as was done so by Konno & Ohmachi (1998) and Arai & Tokimatsu (2004).
The phase difference Arg[h (q) (ω)] between the vertical and horizontal components of Rayleigh waves can be estimated when and only when the dominance of a single mode can be assumed. In that case, we obtain, from eq. (61), where φ R 0 is the central arrival direction of the Rayleigh waves as described in the foregoing section. We can estimate φ R 0 from eqs (D4), (D7)-(D9), for example, while the spectral density G Z0R1 (0, 0; ω) can be obtained from measurement records from a single point. Theoretically, Arg[h(ω)] is expected to be either ±π/2.