Construction of Equivalent Functions in Anisotropic Radon Tomography

Abstract

We consider a real-valued function on a plane of the form

m(x,y,θ)=A(x,y)+Bc(x,y)cos(2θ)+Bs(x,y)sin(2θ)+Cc(x,y)cos(4θ)Cs(x,y)sin(4θ)

that models anisotropic acoustic slowness (reciprocal velocity) perturbations. This “slowness function” depends on Cartesian coordinates and polar angle θ. The five anisotropic “component functions” A (x,y), Bc(x,y), Bs(x,y), Cc(x,y) and Cs(x,y) are assumed to be real-valued Schwartz functions. The “travel time” function d(u, θ) models the travel time perturbations on an indefinitely long straight-line observation path, where the line is parameterized by perpendicular distance u from the origin and polar angle θ; it is the Radon transform of m ( x, y, θ). We show that: 1) an A can always be found with the same d(u, θ) as an arbitrary (Bc,Bs) and/or an arbitrary (Cc,Cs) ; 2) a (Bc,Bs) can always be found with the same d(u, θ) as an arbitrary A, and furthermore, infinite families of them exist; 3) a (Cc,Cs) can always be found with the same d(u, θ) as an arbitrary A, and furthermore, infinite families of them exist; 4) a (Bc,Bs) can always be found with the same d(u, θ) as an arbitrary (Cc,Cs) , and vice versa; and furthermore, infinite families of them exist; and 5) given an arbitrary isotropic reference slowness function m0(x,y), “null coefficients” (Bc,Bs) can be constructed for which d(u, θ) is identically zero (and similarly for Cc,Cs ). We provide explicit methods of constructing each of these “equivalent functions”.

Share and Cite:

Menke, W. (2019) Construction of Equivalent Functions in Anisotropic Radon Tomography. Applied Mathematics, 10, 1-10. doi: 10.4236/am.2019.101001.

1. Introduction

This paper addresses the non-uniqueness of the two-dimensional inverse Radon transform, when the real-valued function m ( x , y , θ ) being transformed is presumed to be anisotropic; that is, varying with polar angle θ as well as with position ( x , y ) . This is a common special case in geotomography, where m ( x , y , θ ) represents perturbations in anisotropic acoustic slowness (reciprocal velocity) on a plane in the Earth with respect to a homogeneous and isotropic background and where its Radon transform represents the corresponding travel time perturbations on indefinitely long straight-line observation paths [1] [2] [3] [4] [5] .

As we will discuss further below, the presence of weak anisotropy results in slowness with three modes of angular variation: isotropic, 2θ and 4θ (with polar angle θ), which are described by a total of five spatially-varying “component functions” [6] [7] [8] . Previously, a modal analysis was used to prove that insufficient information is contained in travel time measurements to uniquely determine all five component functions [9] . Also previously, analytic formula were derived for the spatially-distributed 2θ components equivalent to (in the sense of having the same travel time as) an impulsive isotropic component, and vice versa [10] . While these results indicate that every isotropic mode has a 2θ equivalent, and vice versa, it does not provide a simple method for constructing equivalent modes, and it leaves unresolved the issues regarding 4θ non-uniqueness. We address these issues here.

2. Isotropic Travel Time Tomography with the Radon Transform

Assuming straight line observation paths of indefinite length, the travel time perturbation d ( u , θ ) associated with an isotropic slowness perturbation m ( x , y ) is computed via the Radon transform:

d ( u , θ ) = L ( u , θ ) m ( x , y ) d l R m (1)

The line integral is taken over the straight line L with arc length l , parameterized by its perpendicular distance u to the origin and the counter-clockwise angle θ that the perpendicular makes with the x-axis. We limit our discussions here to two-dimensional functions drawn from Schwatz space [11] in which the Fourier transform is a linear isomorphism [12] [13] . This restriction is usually acceptable in geotomography, where slowness functions can be assumed to rapidly decrease towards zero outside of a restricted area of interest. As usual, we will refer to the travel time d ( u , θ 0 ) with θ 0 held constant as a “projection”. The Projection Slice Theorem [14] shows that the one-dimensional Fourier transform F u of a projection, which takes u into k u , is the two-dimensional Fourier Transform of the slowness function, evaluated on a line of angle θ 0 in the wavenumber plane:

d ¯ ( k u , θ 0 ) = F u d ( u , θ 0 ) = m ¯ ( k u cos ( θ 0 ) , k u sin ( θ 0 ) ) with m ¯ ( k x , k y ) F x F y m (2)

Here the two-dimensional Fourier transform F x F y takes ( x , y ) into ( k x , k y ) . Consequently, a Radon transform has no null space; that is, there is no non-zero isotopic slowness function m ( x , y ) for which R m = 0 . (As we shall show below, the same is not true for anisotropic slowness functions).

Multiplying the Radon transform by a smooth real-valued function F ( θ ) of polar angle θ to yield F ( θ ) d ( r , θ ) only changes the overall scaling of a projection and, because of the Fourier Projection Theorem, corresponds to a fan filter in the ( k x , k y ) domain:

F ( θ ) m ¯ ¯ ( k x , k y ) with tan ( θ ) = k x / k y (3)

A fan filter satisfying | F ( θ ) | 1 does not increase the overall energy of the Fourier transform, so by Plancherel’s theorem [15] the energy in the ( x , y ) -domain function cannot increase. Furthermore, if m ¯ ¯ ( k x , k y ) is a Schwartz function, then so is F ( θ ) m ¯ ¯ ( k x , k y ) , implying that the fan filtered flowness function, say m ( x , y ) F y 1 F x 1 { F ( θ ) m ¯ ¯ ( k x , k y ) } , exists. The fan filter must also obey the symmetry condition F ( θ ) = F ( θ + π ) to ensure that m ( x , y ) is real-valued.

3. Anisotropic Travel Time Tomography with the Radon Transform

The slowness function associated with weak anisotropy [6] [7] [8] [16] :

m ( x , y , θ ) = A ( x , y ) + B C ( x , y ) cos ( 2 θ ) + B S ( x , y ) sin ( 2 θ ) + C C ( x , y ) cos ( 4 θ ) + C S ( x , y ) sin ( 4 θ ) (4)

contains three “modes” of angular variation specified by a total of five spatially-varying “coefficient functions”. The coefficient function A specifies the isotropic mode; the pair of coefficient functions ( B C , B S ) specify the 2θ mode and the pair of coefficient functions ( C C , C S ) specify the 4θ mode. The travel time is the Radon transform of this slowness function:

d ( u , θ ) = R A + R [ B C cos ( 2 θ ) ] + R [ B S sin ( 2 θ ) ] + R [ C C cos ( 4 θ ) ] + R [ C S sin ( 4 θ ) ] (5)

The trigonometric functions can be moved outside the Radon transforms, since they do not vary along the transform’s straight line integration path:

d ( u , θ ) = R A + cos ( 2 θ ) R B C + sin ( 2 θ ) R B S + cos ( 4 θ ) R C C + sin ( 4 θ ) R C S (6)

Note that the term cos ( 2 θ ) R B C implies fan filtering of B ¯ ¯ C ( k x , k y ) (and similarly for the other terms containing trigonometric functions). Because of their 2θ and 4θ dependence, these fan filters obey the symmetry conditions for Fourier transforms of real functions.

4. Isotropic Mode Equivalent to an Anisotropic Mode

Consider a case where only the isotropic component function A ( x , y ) is nonzero, and another case where only the 2θ anisotropic component functions, B c ( x , y ) and B s ( x , y ) , are nonzero. The two cases can be made to imply the same travel times with the choice:

A = R 1 d B with d B cos ( 2 θ ) R B C + sin ( 2 θ ) R B S (7)

Since d B is just the sum of fan-filtered versions of functions whose Radon transforms are presumed to exist, we are assured that its inverse Radon transform exists, too. Hence, an A can always be found that is equivalent to an arbitrary ( B C , B S ) (in the sense of the two having the same travel time perturbation).

Similarly, for a case where only C c ( x , y ) and C s ( x , y ) are nonzero, the equivalent A is:

A = R 1 d C with d C cos ( 4 θ ) R C C + sin ( 4 θ ) R C S (8)

Hence, an A can always be found that mimics a set of ( C C , C S ) . We demonstrate Equation (7) by proposing a ( B C , B S ) , computing the equivalent A, and noting that their Radom transforms are identical (Figure 1).

5. Isotropic Mode Equivalent to an Anisotropic Mode

First, note the identity:

A ( x , y ) = A ( x , y ) cos 2 ( 2 θ ) + A ( x , y ) sin 2 ( 2 θ ) (9)

Now consider one case where only A ( x , y ) is nonzero:

m ( x , y , θ ) = A ( x , y ) cos 2 ( 2 θ ) + A ( x , y ) sin 2 ( 2 θ ) (10)

and another case where only ( B c , B s ) and are nonzero:

m ( x , y , θ ) = B C ( x , y ) cos ( 2 θ ) + B S ( x , y ) sin ( 2 θ ) (11)

By matching terms, two cases can be made to have the same travel time:

R B C = ( cos 2 ( 2 θ ) cos ( 2 θ ) R A ) or B C = R 1 [ cos ( 2 θ ) R A ] (12)

and

R B S = ( sin 2 ( 2 θ ) sin ( 2 θ ) R A ) or B S = R 1 [ sin ( 2 θ ) R A ] (13)

Note that this approach relies upon the zeros in the cos 2 ( 2 θ ) and sin 2 ( 2 θ ) in the numerators cancelling the zeros cos ( 2 θ ) and sin ( 2 θ ) in the denominators, so that the fractions are finite. If we could find a pair of functions c ( 2 θ ) and s ( 2 θ ) with the property c ( 2 θ ) + s ( 2 θ ) = 1 and with zeros in the appropriate locations, we could write A ( x , y ) = A ( x , y ) c ( 2 θ ) + A ( x , y ) s ( 2 θ ) , and:

B C = R 1 ( c ( 2 θ ) cos ( 2 θ ) R A ) and B S = R 1 ( s ( 2 θ ) sin ( 2 θ ) R A ) (14)

would follow. Evidentially, many such pairs of function exist, including:

c ( 2 θ ) = cos 2 n ( 2 θ ) and s ( 2 θ ) = 1 cos 2 n ( 2 θ ) (15)

Figure 1. The A equivalent to an exemplary ( B C , B S ) . ((a) and (b)) The exemplary ( B C , B S ) ; (c) The A equivalent to this ( B C , B S ) , computed by direct application of Equation (7); (d) The travel times associated with ( B C , B S ) ; (e) The travel times associated with A; (f) The travel time difference. The Radon transform is computed via quadrature along lines and the inverse Radon transform via a Fourier-domain method based on the Projection Slice Theorem.

and

c ( 2 θ ) = 1 sin 2 n ( 2 θ ) and s ( 2 θ ) = sin 2 n ( 2 θ ) (16)

with n 1 . The s ( 2 θ ) function in Equation (15) has sin 2 ( 2 θ ) as a factor, as can be demonstrated by writing the unity term as [ cos 2 ( 2 θ ) + sin 2 ( 2 θ ) ] n , expanding it with the Binomial Theorem [17] and subtracting off the second term:

s ( 2 θ ) = ( cos 2 ( 2 θ ) + sin 2 ( 2 θ ) ) n cos 2 n ( 2 θ ) = k = 0 n ( k n ) cos 2 ( n k ) ( 2 θ ) sin 2 k ( 2 θ ) cos 2 n ( 2 θ ) = k = 1 n ( k n ) cos 2 ( n k ) ( 2 θ ) sin 2 k ( 2 θ ) = sin 2 ( 2 θ ) k = 1 n ( k n ) cos 2 ( n k ) ( 2 θ ) sin 2 ( k 1 ) ( 2 θ ) (17)

Similarly, the c ( 2 θ ) in Equation (16) has cos 2 ( 2 θ ) as a factor. We conclude that a ( B C , B S ) can always be found that is equivalent to an arbitrary A, and furthermore, that infinite families of such pairs exist. By replacing 2θ with 4θ in the above argument, we conclude that a ( C C , C S ) can always be found that is equivalent to an arbitrary A, and furthermore, that infinite families of such pairs exist. We demonstrate Equations (14) and (16) by proposing a ( B C , B S ) , computing the equivalent A, and showing that their Radon transforms are identical (Figure 2).

Figure 2. Two ( B C , B S ) pairs equivalent to the same exemplary A. (a) An exemplary A; ((b) and (c)) One possible ( B C , B S ) with a travel time equivalent to A, computed according to Equation (15) with n = 1 ; (d) The same exemplary A as in (a); ((e) and (f)) Another possible ( B C , B S ) pair with a travel time equivalent to A, computed according to Equation (14) and (16) with n = 2. The equivalence of the travel times has been verified by numerical calculation (not shown).

6. Two-Theta Mode Equivalent to a Four-Theta Mode, and Vice Versa

We can find the ( C C , C S ) equivalent to a ( B C , B S ) by using the method of Section 4 to find the A equivalent to ( B C , B S ) and then using the method of Section 5 to find the ( C C , C S ) equivalent to that A. Similarly, we can find the ( B C , B S ) equivalent to an ( C C , C S ) by using the method of Section 4 to find the A equivalent to ( C C , C S ) and then using the method of Section 5 to find the ( B C , B S ) equivalent to that A. Recalling that the method of Section 5 is non-unique, we conclude that a ( B C , B S ) can always be found that mimics a ( C C , C S ) , and vice versa; and furthermore, infinite families of them exist.

As an example, we find a ( C C , C S ) equivalent to an exemplary ( B C , B S ) and a ( B C , B S ) equivalent to an exemplary ( C C , C S ) (Figure 3).

7. Anisotropic Null Slowness Functions

We have demonstrated in Section 5 that two slowness functions, say m 1 ( x , y , θ ) and m 2 ( x , y , θ ) , each with only a 2θ mode, can have the same travel time; that is R m 1 = R m 2 . Because the Radom transform is a linear operator, the “null slowness function” m n u l l = m 1 m 2 must have zero Radon transform; that is, R m n u l l = 0 [18] . (This finding does not violate the principle that no non-zero isotropic slowness function can have an identically zero Radon transform,

Figure 3. A ( B C , B S ) equivalent to a ( C C , C S ) , and vice versa. ((a) and (b)) An exemplary ( B C , B S ) ; ((c) and (d)) An exemplary ( C C , C S ) ; ((e) (f)) A ( C C , C S ) equivalent to the exemplary ( B C , B S ) ; ((g) (h)) A ( B C , B S ) equivalent to the exemplary ( C C , C S ) . The equivalence of the travel times has been verified by numerical calculation (not shown).

because the travel time perturbations of a 2θ image is the sum of two fan-filtered Radon transforms). The ( B C , B S ) of an exclusively 2θ null slowness function can be constructed from an isotropic reference slowness function m 0 ( x , y ) in the following way:

B C = R 1 ( sin ( 2 θ ) d 0 ) and B S = R 1 ( cos ( 2 θ ) d 0 ) with d 0 = R m 0 (18)

The travel time d B associated with this ( B C , B S ) is identically zero:

d B = R [ cos ( 2 θ ) B C + sin ( 2 θ ) B s ] = cos ( 2 θ ) R B C + sin ( 2 θ ) R B s = cos ( 2 θ ) sin ( 2 θ ) d 0 sin ( 2 θ ) cos ( 2 θ ) d 0 = 0 (19)

This argument can be extended to the 4θ mode simply by replacing 2θ with 4θ. Hence, given an isotropic reference slowness function m 0 ( x , y ) , a 2θ null slowness function, with coefficients ( B C , B S ) , can be constructed for which the travel time perturbation is identically zero (and similarly for 4θ null slowness function). We demonstrate Equation (18) by building ( B C , B S ) and ( C C , C S ) that correspond to null slowness functions from an exemplary isotropic reference slowness function (Figure 4).

8. Conclusions

A weakly-anisotropic acoustic slowness function on a plane has three modes of angular variability: isotropic, 2θ with 4θ (where θ is polar angle) that are described by a total of five spatially-varying coefficients. The Radon transform of the slowness determines travel time perturbations with respect to a homogeneous and isotropic background, observed in an idealized experiment with indefinitely long straight line observation paths. In this paper, we show that an arbitrary travel time perturbation can be matched by a slowness function with

Figure 4. ( B C , B S ) and ( C C , C S ) corresponding to null slowness functions. (a) An exemplary isotropic reference function m 0 ( x , y ) ; ((b) and (c)) A ( B C , B S ) built from the reference function with zero travel time using Equation (18); (d) Another exemplary isotropic reference slowness function m 0 ( x , y ) ; ((e) and (f)) A ( C C , C S ) built from the reference function with zero travel time. The zero travel times have been verified by numerical calculation (not shown).

only one of any of the three mode of angular variability. Consequently, anisotropic slowness functions cannot be uniquely determined from these travel times. Many slowness functions, equivalent in the sense of having the same travel times, exist and can be computed from one another by applying straightforward procedures. Furthermore, slowness functions with exclusively two-theta variability are non-unique; infinite families of equivalent coefficients exist. The same is true of exclusively four-theta slowness functions. Consequently, some exclusively 2θ and 4θ slowness functions have identically zero travel time and, consequentially, are inherently undetectable.

In geotomography, non-uniqueness is usually handled through the addition of prior information that regularizes the problem by singling out one “best” slowness function from among all slowness functions consistent with the travel time data [18] [19] [20] . The choice of regularizing constraint is informed both by knowledge of the physical system that is being imaged and nature of the non-uniqueness. Our results are relevant to the second issue. When the true slowness function is believed on prior grounds to contain compact equidimensional heterogeneities, the equivalent slowness functions (now considered artifacts) tend to be more spatially-distributed and oscillatory, owing to the effect of the fan filtering. Smoothness regularization will tend to suppress these artifacts, as long as the smoothing response is not, itself, oscillatory. Thus, first-derivative regularization is preferred over second-derivative regularization, because the latter has an oscillatory response [21] . Smallness regularization [22] of slowness functions believed to contain strong two-theta and four-theta modes should be avoided, because it will tend to roughen them by suppressing the contribution of null components that while not observable, contribute to the smoothness of the slowness function that is being imaged.

Acknowledgements

I thank the students in the Columbia University 2018 Anisotropy Seminar that I taught for helpful discussion. This research was supported by the US National Science Foundation under awards OCE-0426369 and EAR 11-47742.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1] Hearn, T.M. (1996) Anisotropic Pn Tomography in the Western United States. Journal of Geophysical Research, 101, 8403-8414.
https://doi.org/10.1029/96JB00114
[2] Wu, H. and Lees, J.M. (1999) Cartesian Parameterization of Anisotropic Traveltime Tomography. Geophysical Journal International, 137, 64-80.
https://doi.org/10.1046/j.1365-246x.1999.00778.x
[3] Pei, S., Zhao, J., Sun, Y., Xu, Z., Wang, S., Liu, H., Rowe, C.A., ToksÖz, M.N. and Gao, X. (2007) Upper Mantle Seismic Velocities And anisotropy in China Determined through Pn and Sn Tomography. Journal of Geophysical Research, 112, Article ID: B05312.
[4] Menke, W., Zha, Y., Webb, S.C. and Blackman, D.K. (2015) Seismic Anisotropy Indicates Ridge-Parallel Asthenospheric Flow beneath the Eastern Lau Spreading Center. Journal of Geophysical Research, 120, 976-992.
https://doi.org/10.1002/2014JB011154
[5] Eddy, C.L., EkstrÖm, G., Nettles, M. and Gaherty, J.B. (2018) Age Dependence and Anisotropy of Surface-Wave Phase Velocities in the Pacific. Geophysical Journal International, 216, 640-658.
https://doi.org/10.1093/gji/ggy438
[6] Backus, G.E. (1965) Possible Form of Seismic Anisotropy of the Upper Mantle under the Oceans. Journal of Geophysical Research, 70, 3429-3439.
https://doi.org/10.1029/JZ070i014p03429
[7] Thomsen, L. (1986) Weak Elastic Anisotropy. Geophysics, 51, 1954-1966.
https://doi.org/10.1190/1.1442051
[8] Aki, K. and Richards, P.G. (2009) Quantitative Seismology. 2nd Edition, University Science Books, Herndon, Virginia.
[9] Mochizuki, E. (1997) Nonuniqueness of Two-Dimensional Anisotropic Tomography. Seismological Society of America Bulletin, 87, 261-264.
[10] Menke, W. (2015) Equivalent Heterogeneity Analysis as a Tool for Understanding the Resolving Power of Anisotropic Travel Time Tomography. Seismological Society of America Bulletin, 105, 719-733.
https://doi.org/10.1785/0120140150
[11] Reed, M. and Simon. B. (1980) Methods of Modern Mathematical Physics: Functional Analysis I. Revised and Enlarged Edition, Academic Press, San Diego.
[12] HÖrmander, L. (1990) The Analysis of Linear Partial Differential Operators I. Distribution Theory and Fourier Analysis. 2nd Edition, Springer-Verlag, Berlin.
[13] Stein, E.M. and Shakarchi, R. (2003) Fourier Analysis: An Introduction. Princeton Lectures in Analysis I. Princeton University Press, Princeton New Jersey.
[14] Deans, S.R. (1993) The Radon Transform and Some of Its Applications. Revised Edition, Dover, Mineola, New York.
[15] Bracewell, R. (1999) The Fourier Transform and Its Applications. 3rd Edition, McGraw Hill, New York.
[16] Smith, M.L. and Dahlen, F.A. (1973) Azimuthal Dependence of Love and Rayleigh-Wave Propagation in a Slightly Anisotropic Medium. Journal of Geophysical Research, 78, 3321-3333.
https://doi.org/10.1029/JB078i017p03321
[17] Seely, R.T. (1973) Calculus of One and Several Variables. 2nd Edition, Scott Foresman, Glenview.
[18] Menke, W. (2018) Geophysical Data Analysis: Discrete Inverse Theory. Fourth Edition, Elsevier, Amsterdam.
[19] Tarantola, A. and Valette, B. (1982) Inverse Problems = Quest for Information. Journal of Geophysics, 50, 159-170.
[20] Menke, W. and Menke, J. (2016) Environmental Data Analysis with MATLAB. Second Edition, Elsevier, Amsterdam.
[21] Menke, W. and Eilon, Z. (2015) Relationship between Data Smoothing and the Regularization of Inverse Problems. Pure and Applied Geophysics, 172, 2711-2726.
https://doi.org/10.1007/s00024-015-1059-0
[22] Levenberg, K. (1944) A Method for the Solution of Certain Non-Linear Problems in Least Squares. Quarterly of Applied Mathematics, 2, 164-168.
https://doi.org/10.1090/qam/10666

Copyright © 2024 by authors and Scientific Research Publishing Inc.

Creative Commons License

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.