Inversion of the elliptical Radon transform arising in migration imaging using the regular Radon transform

Inversion of the elliptical Radon transform arising in migration imaging using the regular Radon transform

Sunghwan Moon111shmoon@unist.ac.kr and Joonghyeok Heo
Department of Mathematical Sciences, Ulsan National Institute of Science and Technology (UNIST), Ulsan 44919, Korea
UM Energy Institute, University of Michigan, Ann Arbor, MI 48109, U.S.
Abstract

In recent years, many types of elliptical Radon transforms that integrate functions over various sets of ellipses/ellipsoids have been considered, relating to studies in bistatic synthetic aperture radar, ultrasound reflection tomography, radio tomography, and migration imaging. In this article, we consider the transform that integrates a given function in over a set of ellipses (when ) or ellipsoids of rotation (when ) with foci restricted to a hyperplane. We show a relation between this elliptical Radon transform and the regular Radon transform, and provide the inversion formula for the elliptical Radon transform using this relation. Numerical simulations are performed to demonstrate the suggested algorithms in two dimensions, and these simulations are also provided in this article.

Keywords: Radon transform; tomography; elliptical; migration; seismic reflection

AMS Subject Classification: 44A12; 65R10; 92C55

1 Introduction

Radon-type transforms that assign to a given function its integrals over various sets of ellipses/ellipsoids arise in migration imaging under an assumption that the medium is acoustic and homogeneous. The aim of migration is to construct an image of the inside of the earth from seismic reflections recorded at its surface [9, 25]. A graphical approach called classical migration was developed systematically by Hagedoorn [11]. Classical migration had been abandoned after the wave-equation method was introduced by Claerbout [5] in 1971. Gazdag and Sguazzero pointed out that the classical migration procedures that existed at that time were not based on a completely sound theory [9]. However, a correct construction for the wave-equation method was often difficult to find because the experiment data did not fit into a single wave equation. To adapt the diffraction stack to borehole seismic experiments, a new approach to seismic migration was found. This approach gave classical migration a sound theory. After that discovery, classical migration has attracted many researchers in the field. The underlying idea is that seismic data in the far field can be regarded as if the data are coming from integrals of the earth’s acoustic scattering potential over surfaces determined by the velocity model [16]. These Radon-type transforms relate to migration imaging as well as Bistatic Synthetic Aperture Radar (BiSAR) [2, 6, 13, 14, 32], Ultrasound Reflection Tomography (URT) [1, 10, 26], and radio tomography [27, 28, 29].

Because of these applications, there have been several papers devoted to the topic of elliptical Radon transforms. The family of ellipses with one focus fixed at the origin and the other one moving along a given line was considered in [14]. In the same paper, the family of ellipses with a fixed focal distance was also studied. The authors of [1, 10] dealt with the case of circular acquisition, when the two foci of ellipses with a given focal distance are located on a given circle. A family of ellipses with two moving foci was also handled in [6]. Radio tomography, which uses a wireless network of radio transmitters and receivers to image the distribution of attenuation within the network, was discussed in [27, 28, 29]. They approximated the obtained signal by the volume integral of the attenuation over this ellipsoid. One work [17, 18] derived two inversion formulas of this volume integral of the attenuation over this ellipsoid and studied its properties. Many works found an approximate inversion for elliptical Radon transforms.

Here we consider migration imaging and introduce a new type of an elliptical Radon transform obtained by restricting the position of the source and receiver in migration imaging. We find an explicit inversion formula for this elliptical Radon transform arising in migration imaging which is the line/area integral of the function over the ellipse/ellipsoid with foci restricted to a hyperplane.

The rest of this paper is organized as follows. The problem of interest is stated precisely and the elliptical Radon transform is formulated in section 2. In section 3, we show how to reduce the elliptical Radon transform to the regular Radon transform. The numerical simulation to demonstrate the suggested 2-dimensional algorithm is presented in section 4.

2 Formulation of the problem

Let and represent 3-dimensional source and receiver positions, respectively. For fixed points , an isochron surface 222Hagedoorn called this set a surface of maximum concavity (see [11]). is a surface consisting of image points associated by the travel time function , which gives the travel time between two points and with a known velocity (see [16]).

Mathematically, can be described as the set of the image points satisfying the constraint that the total travel time from the source though the image point to the receiver is constant and equal to . The isochron surface can be represented as

Seismic experiments yield data which are functions of the source position , receiver position , and time . Assuming an object function on , the data is modelled by the integral of over , i.e.,

The two detectors have nonzero sizes and time is also passing, so it is reasonable to assume that what we measure is the “average”-concentrated near the location of the detectors and nearly zero sufficiently far away from the detectors-over a small region of space and a small time interval preceding the time . In mathematical terms, our data can be written as

where is the Dirac delta function. (Most works [1, 2, 6, 10, 13, 14] dealing with an elliptical Radon transform consider the arc length measure instead.) Miller, Oristaglio, and Beylkin suggested this model and approximately recovered the object function by appropriately weighted back projection of the data in [16]. However, we present the explicit formula for reconstructing using the specific minimum data of by changing the variables of positions of the two detectors.

The total dimension of the data is 7. To reduce the overdeterminacy, we assume that and are located on a hyperplane, say and also that the difference vector is parallel to a given line, say the -axis. (In 2 dimensions, the difference vector is automatically parallel to the -axis because the hyperplane is the line.) Hence, for any , we write and . Also, assuming a velocity is set to a constant 1, becomes the distance between and . Since satisfies , a point can be described by the ellipsoid equation

Let us choose a constant between 0 and 1. We consider the fixed value of satisfying the equation , which is a new additional requirement. Then a point can be also described by the equation

If we set , then our data become

where . Here, we restrict the positions of the source and receiver in migration imaging, say, and for . In general, reducing the space of positions for two devices is more useful and economical. If the function is odd with respect to , then is equal to zero. We thus assume the function to be even with respect to : where . We call an elliptical Radon transform, since this is the surface area integral of over the set of these ellipsoids. We generalize this -dimension setup to -dimension and define a more general form of an elliptical Radon transform.

Definition 2.1.

Let be given numbers, let denote the diagonal matrix with diagonal entries , and let be a locally integrable function on , even with respect to . The elliptical Radon transforms of is defined by

Note that we do not need the condition any more. If all , , are equal to 1, the elliptical Radon transform becomes the spherical Radon transform with the centers of the sphere of integration located on the hyperplane, a well-studied transform (see [3, 4, 7, 8, 15, 12, 19, 20, 21, 22, 23, 24, 30, 31]).

3 Inversion formula

Here we assume that the object function does not touch the detectors; that is, the support of does not intersect the hyperplane where the two detectors are located.

To obtain an inversion formula for the elliptical Radon transform, we manipulate . By changing the variables , we can write

(1)

where and . Here is the area measure of the unit sphere.

Let and . We define a map by

where is the diagonal matrix.

Proposition 3.1.
  • The map is a bijection with the inverse map defined by

  • We have

    where .

The map transforms an ellipsoid into a hyperplane​. Changing variables using this map plays a critical role in reducing the elliptical Radon transform to the regular Radon transform.

Proof.

We can easily check that for and for , so is a bijection.

Consider for with :

(2)

We define the function on by

where . This is equivalent for , to

where . By the evenness of , we have

(3)

Let the regular Radon transform be defined by

where and for

This can be represented by

(4)

for and . The following theorem shows the relation between and .

Theorem 3.2.

Let be even in and have compact support in . Then we have

(5)

Again, .

Proof.

Combining two equations (1) and (3), we have

Let us consider the variable change

The Jacobian of this transformation is , so

Since has compact support in , we have

(6)

We recognize the right hand side as the integral along the hyperplane perpendicular to

with (signed) distance

from the origin. In this case, the measure for the hyperplane becomes Setting

in equation (4) we have the desired formula from equation (6). ∎

Using the projection slice theorem for the regular Radon transform, we obtain an analog of the projection slice theorem for the elliptical Radon transform :

Theorem 3.3.

Let be even in and have compact support in . Then we have for ,

where is the -dimensional Fourier transform of .

Proof.

The projection slice theorem implies

To use this theorem, we multiply equation (5) by and integrate from to infinity:

Since has compact support in , the plane perpendicular to with distance from the origin less than does not intersect the compact support of . Hence, is equal to

where in the second line, we used Theorem 3.2. Setting and completes our proof. ∎

Taking the inverse Fourier transform of a function and using equation (3), we have the following inversion formula:

Corollary 3.4.

Let be even in and have compact support in . Then can be reconstructed as follows:

(7)

where

Proof.

We get for ,

where in the first line, we used equation (3) and in the last line, we used Theorem 3.3 and

The evenness of in gives us

Changing the variables , is equal to

By definition of , can be determined through

When , equation (7) becomes

since in the distribution sense,

4 2-dimensional numerical implementation

Here we discuss the results of 2-dimensional numerical implementations.

In 2-dimension, equation (5) becomes

(8)

Setting and

in equation (8), we have

where

(While is used in Theorem 3.2, is used to utilize the built-in function in Matlab. Hence a small change is required.) Again, since has compact support in , the plane perpendicular to with distance from the origin less than does not intersect the compact support of . Therefore, we have

(9)

First of all, and were set to be 0.8 and 1, respectively. In the experiments presented here we use the phantom shown in Figure 1 (a). The phantom, supported within the rectangle , is the sum of eight characteristic functions of disks. Notice that it has to be even with respect to the -axis and there are four characteristic functions of disks centered at , , , and with radii 0.2 0.15 0.05 and 0.05, whose values are 1, 0.5, 1.5, and 2 above the -axis. Hence we also include their reflection below the -axis. (Actually, our phantom has support in

This implies that has support in the unit ball and this makes it sufficient to consider the range in .) The images are used in Figure 1. To reconstruct the image in Figure 1 (b), we have projections for and in equation (9). All projections are computed by numerical integration. After finding the function using the usual inversion code for the regular Radon transform, we obtain the function using equation (3). (When using the inversion code for the regular Radon transform, the built-in function “iradon” in Matlab was used. The function “iradon” is the inversion of the built-in function “radon” in Matlab which considers the number of the pixels where the line passes through. Thus when computing , we also considers the number of the pixels where the ellipse passes through. We used the default version of the function in which the filter, whose aim is to deemphasize high frequencies, is set to the Ram-Lak filter and the interpolation is set to be linear.)

While Figure 1 (b) demonstrates the image reconstructed from the exact data, Figure 1 (c) shows the absolute value of the reconstruction from noisy data. The noise is modeled by normally distributed random numbers and this is scaled so that its norm was equal to 5 of the norm of the exact data. In Figure 1 (c) the noisy data is modeled by adding the noise values scaled to 5 of the norm of the exact data to the exact data. In Figure 2, the surface plots of Figure 1 (a) and (b) are provided. Another phantom and reconstruction are shown in Figure 3. The reconstructed two dimensional data sets consisting of 256256 projections using the implemented inversion formula have been computed in less than one second (around 0.3 second) on a Intel(R) Cor(TM) i5-3470 CPU @ 3.20 GHz.

Figure 1: Reconstructions in 2 dimensions: (a) the phantom, (b) the reconstruction form exact data, and (c) the reconstruction from noisy data
Figure 2: Surface plots (a) the phantom and (b) the reconstruction from exact data
Figure 3: Reconstruction in 2D: (a) the phantom and (b) the reconstruction from exact data

5 Conclusion

This paper is devoted to the study of the elliptical Radon transform arising in seismic imaging. We provide an inversion formula for the elliptical Radon transform by reducing this transform to the regular Radon transform. Also, we demonstrate our algorithm by providing numerical simulations.

Acknowledgements

The first author thanks C. Jung for fruitful discussions. The first author was supported in part by Basic Science Research Program through the National Research Foundation of Korea(NRF) funded by the Ministry of science, ICT and future planning (2015R1C1A1A02037662)).

References

  • [1] G. Ambartsoumian, J. Boman, V.P. Krishnan, and E.T. Quinto. Microlocal analysis of an ultrasound transform with circular source and receiver trajectories. Contemporary Mathematics, 598, 2013.
  • [2] G. Ambartsoumian, R. Felea, V.P. Krishnan, C. Nolan, and E.T. Quinto. A class of singular Fourier integral operators in Synthetic Aperture Radar imaging. Journal of Functional Analysis, 264(1):246 – 269, 2013.
  • [3] L. Andersson. On the determination of a function from spherical averages. SIAM Journal on Mathematical Analysis, 19(1):214–232, 1988.
  • [4] A. Beltukov. Inversion of the spherical mean transform with sources on a hyperplane. arXiv:0919.1380v1, 2009.
  • [5] J. Claerbout. Toward a unified theory of reflector mapping. GEOPHYSICS, 36(3):467–481, 1971.
  • [6] J.D. Coker and A.H. Tewfik. Multistatic SAR image reconstruction based on an elliptical-geometry Radon transform. In Waveform Diversity and Design Conference, pages 204 –208, June 2007.
  • [7] R. Courant and D. Hilbert. Methods of Mathematical Physics. Number V. 2 in Methods of Mathematical Physics. Interscience Publishers, 1962.
  • [8] J. Fawcett. Inversion of n-dimensional spherical averages. SIAM Journal on Applied Mathematics, 45(2):336–341, 1985.
  • [9] J. Gazdag and P. Sguazzero. Migration of seismic data. Proceedings of the IEEE, 72(10):1302–1315, 1984.
  • [10] R. Gouia-Zarrad and G. Ambartsoumian. Approximate inversion algorithm of the elliptical Radon transform. In Mechatronics and its Applications (ISMA), 2012 8th International Symposium on, pages 1 –4, April 2012.
  • [11] J. G. Hagedoorn. a process of seismic reflection interpretation*. Geophysical Prospecting, 2(2):85–127, 1954.
  • [12] J. Klein. Inverting the spherical Radon transform for physically meaningful functions. arXiv:math/0307348, 2003.
  • [13] V. Krishnan and E.T. Quinto. Microlocal aspects of common offset synthetic aperture radar imaging. Inverse Problems and Imaging, 5:659–674, 2011.
  • [14] V.P. Krishnan, H. Levinson, and E.T. Quinto. Microlocal analysis of elliptical Radon transforms with foci on a line. In Irene Sabadini and Daniele C Struppa, editors, The mathematical legacy of Leon Ehrenpreis, volume 16 of Springer Proceedings in Mathematics, pages 163–182. Springer Milan, 2012.
  • [15] M.M. Lavrentiev, V.G. Romanov, and V.G. Vasiliev. Multidimensional inverse problems for differential equations. Lecture notes in mathematics. Springer-Verlag, 1970.
  • [16] D. Miller, M. Oristaglio, and G. Beylkin. A new slant on seismic imaging: Migration and integral geometry. Geophysics, 52(7):943–964, 1987.
  • [17] S. Moon. Properties of some integral transforms arising in tomography. Dissertation, Texas A&M University, 2013.
  • [18] S. Moon. On the determination of a function from an elliptical Radon transform. Journal of Mathematical Analysis and Applications, 416(2):724 – 734, 2014.
  • [19] E.K. Narayanan and Rakesh. Spherical means with centers on a hyperplane in even dimensions. Inverse Problems, 26(3):035014, 2010.
  • [20] F. Natterer and F. Wübbeling. Mathematical methods in image reconstruction. SIAM Monographs on mathematical modeling and computation. SIAM, Society of industrial and applied mathematics, Philadelphia (Pa.), 2001.
  • [21] S. Nilsson. Application of fast backprojection techniques for some inverse problems of integral geometry. Linköping studies in science and technology: Dissertations. Department of Mathematics, Linköping University, 1997.
  • [22] S. J. Norton. Reconstruction of a reflectivity field from line integrals over circular paths. The Journal of the Acoustical Society of America, 67(3):853–863, 1980.
  • [23] V. P. Palamodov. Reconstructive integral geometry, volume 98. Springer, 2004.
  • [24] N.T. Redding and G.N. Newsam. Inverting the circular Radon transform. DTSO Research Report DTSO-Ru-0211, August 2001.
  • [25] E.A. Robinson. Migration of geophysical data. International Human Resources Development Corp., 1983.
  • [26] S. Roy, V. P. Krishnan, P. Chandrashekar, and AS V. Murthy. An efficient numerical algorithm for the inversion of an integral transform arising in ultrasound imaging. Journal of Mathematical Imaging and Vision, pages 1–14, 2014.
  • [27] J. Wilson and N. Patwari. Through-wall tracking using variance-based radio tomography networks. ArXiv e-prints, September 2009.
  • [28] J. Wilson and N. Patwari. Radio tomographic imaging with wireless networks. IEEE Transactions on Mobile Computing, 9(5):621 –632, May 2010.
  • [29] J. Wilson, N. Patwari, and F.G. Vasquesz. Regularization methods for radio tomographic imaging. Virginia Tech Symposium on Wireless Personal Communications, June 2009.
  • [30] C.E. Yarman and B. Yazici. Inversion of circular averages using the Funk transform. In Acoustics, Speech and Signal Processing, 2007. ICASSP 2007. IEEE International Conference on, volume 1, pages I–541–I–544, 2007.
  • [31] C.E. Yarman and B. Yazici. Inversion of the circular averages transform using the Funk transform. Inverse Problems, 27(6):065001, 2011.
  • [32] C.E. Yarman, B. Yazici, and M. Cheney. Bistatic synthetic aperture radar imaging for arbitrary flight trajectories. Image Processing, IEEE Transactions on, 17(1):84–93, Jan 2008.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
65291
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description