Application of Iterative Reverse Time Migration Procedure on Transcranial Thermoacoustic Tomography Imaging

high quality and

It applies the principle of the thermoacoustic effect [6], an effect that arises from the combination of the pressure oscillations of a sound wave with the accompanying adiabatic temperature oscillations [7]. Within TAT an input microwave impulse stimulates thermo-expansion in biotissue and consequently generates acoustic waves (P wave) to be recorded by transducers arranged outside the tissue. When the tissue is relatively uniform, the initial local acoustic amplitude is approximately proportional to the absorption ratio of the microwave [1], [8], [9]. Consequently, the TAT imaging problem is to retrieve the distribution of the initial acoustic amplitude based on recorded acoustic wave energy, and thus it can be characterized as a problem of "multiple sources localization".
In recent years, imaging by using thermoacoustic effect has raise increasing concerns [10], [11]. In which TAT has also been well studied for brain imaging [8], [9], [12]. The reasons include the following: first and far more important, TAT takes advantage of deep penetration of the electromagnetic impulse and the high resolution of the ultrasonic wave for deep imaging. Second, brain tissue is fundamentally uniform and isotropic. For example, the acoustic velocity of human brain is narrowly ranged within 1483-1521 m/s [13]. Acoustic wave propagation inside the brain is very close to one-way line-of-sight transmission without much aberration. The most commonly used imaging algorithm of TAT is to back-project the wave energy recorded by each transducer along the ray paths to all possible locations inside the imaging domain [1], [9]. The back-projection method [15]- [18] has been reported to be well performed on its simple scheme with high cost-effectiveness [19]. However, the approximation of oneway sight transmission used in the back-projection algorithm is no longer valid when high velocity contrast exists, for example, when the skull with an acoustic speed of 2500-2900 m/s [20] is included in the imaging domain. For compensating skull-related aberration Jin et al. [13] have developed a strategy based on the approximation of ray-tracing; however, it may still suffer from the difficulty of accurately calculating Green's function when the velocity structure becomes more irregular.
Currently reverse-time migration (RTM) has emerged as a more precise and powerful imaging tool in the exploration geophysics community [21]- [24]. RTM takes full advantage of the wave equation that includes all the dynamic features of a propagating wave field. Different from back-projection, RTM is based on the insensitivity of the wave equation's solution to the directionality of time. During RTM, by solving the wave equation with either the finitedifferences time domain (FDTD) method [25] or the pseudo-spectral time domain method [26], all transducers act as a virtual source by broadcasting their own records back to the domain in a time-reversed manner. If the velocity model is precise, the reversed-time wave field should converge and be enhanced at the origins of the to-be-imaged structures. Previous studies [19], [27], [28] have compared back-projection method and RTM by using complex velocity structure models, and have confirmed better results by using RTM. However, with the advantage of RTM outlined previously, the barriers to RTM adoption are also evident. Compared with back-projection, RTM appears to be more sensitive to accuracy of given velocity model. Once the velocity model is inaccurate, the quality of imaging results from RTM may be insignificant in comparison to back-projection [19], [27]. Consequently the solutions for creating accurate velocity model become extremely vital to guarantee the imaging quality of RTM. For estimate velocity model correctly before RTM, to the best of author's horizon there are two types of approaches. One is tomography, such as ultrasonic transmission tomography. Although it has been reported to be effectiveness for velocity correction from Jin and Wang in 2006 [29], prior of RTM detection, additional laboratory experiment has to be conducted for tomography data collection. This setup obviously increases the complexity of imaging procedure and financial cost. The other approach is an iterative imaging procedure developed by Whitmore [23]. In his approach for estimate the velocity correctly, an initial guess of the velocity model is setup by incorporating all known external velocity information as the first step. After that the TAT image result is derived by RTM, and a comparison between the reconstructed image and the velocity model is made. The differences derived by this comparison are attempted to be eliminated by generating the new velocity model. This procedure is repeated until both RTM results and velocity models are unchangeable. This iterative RTM procedure will be applied in this work.
Based on our previous studies [28], [30], [31] in this paper we closely concentrate on two topics: 1) comparison between back-projection and RTM on transcranial TAT; 2) Application of iterative RTM procedure on velocity correction. The synthetic dataset derived from a twodimensional (2D) brain model and real datasets acquired from laboratory experiments by Xu and Wang [9] (hereinafter referred as XW06) on rhesus monkey heads will be used for test case. The comparison of imaging results derived by both methods supports the finding of previous studies that RTM is superior to back-projection in imaging quality and accuracy. Meanwhile by analysis the imaging results derived by this iterative RTM procedure, conclusion shows that iterative RTM procedure should be a promising approach for achieving transcranial TAT. This paper is organized as follows. Firstly back-projection and RTM are briefly introduced in Section 2. Section 3 illustrates the validation of RTM and back-projection methods through both synthetic data and laboratory data. The iterative RTM procedure will be described in detail and demonstrated with the applications to synthetic data and real laboratory data in Section 4. Section 5 is detailed discussion and analysis on the results coming from Section 3 and Section 4. Finally, in Section 6 we restate the major findings as conclusion.

Back-projection and reverse time migration
In this section we briefly introduce back-projection (or named Kirchhoff migration in exploration seismic engineering) and RTM, two different migration approaches used for transcranial imaging. Generally speaking migration is an inversion operation involving rearrangement of seismic or acoustic information elements from time to depth domain so that reflections and diffractions are plotted at their true locations [32]. It plays a function of moving dipping events to its original location and increasing the accuracy of the image. From seismic exploration the seismic wave is generated by the source on the ground to propagate downward the earth. Once it hits the non-continuity along sub-surface, reflection wave is generated and detected by transducers located on the ground [33]. The misfits exist between observed data and subsurface's real location even sources and receivers are perfectly overlapped. The reasons include: 1) Due to the low frequency and non-infinitesimal width of the input beam, the reflection wave could propagate to multiple direction rather than following the single line; 2) Due to the possible non-horizontal structure owned by sub-surfaces, their locations could be "shifted" to wrong position. As Figure 1 shows, combinations of source and receiver are collocated (black solid dots) along the ground line AB. Suppose the velocity is uniform underground, the slope reflectors along CD is shifted incorrectly to C'D' from the results collected by transducers with wave's travel time unchanged. And therefore, the purpose of migration is to shift C'D' back to CD, and hence the correct structure can be reached.
The migration problem can be transformed as a more straight-forward case when it is applied on TAT imaging. The purpose of migration on TAT imaging is to shift the locations of different acoustic sources to their correct locations. In Figure 2, the shadowed area is the to-be-imaged tissue, in which each point is an independent acoustic source after receiving microwave impulse. The dashed circle surrounded shows the location of receiver array. Suppose the velocity distribution in all 2D space is relatively uniform, and tissue's acoustic wavelet is a single narrow pulse, the possible location of certain source S is constrained on a circle with the center as transducer T N . Its radius is the distance calculated by the wave propagation time took from S to T N . During the TAT imaging, for achieving the exact location of S, back-projection and RTM are applied by different principles. Back-projection can be summed as a kind of spatial integration method. It redistributes the received acoustic energies and makes summation in space. In Figure 2, the results observed by different transducers are redistributed along different circles in to-be-imaged domain. If there is a large number of transducers, the amplitude at the location S will be enhanced due to interference of multi circles, otherwise the amplitudes should be minimum [33]. This approach is easily to achieve when the velocity  Application of migration on TAT imaging. T 1 and T N are any two transducers of transducers' array shown as dash circle surrounded with the target tissue (gray shadow). By applying back-projection, the acoustic amplitude of S will be enhanced by multiple transducers. model is uniform; However when the velocity distribution inside the to-be-imaged domain is complex, the possible location of S will be distorted from circle to an irregular shape. Combined with multiple scattering, the assumption of one-way sight transmission in back-projection is no longer established. Different from back-projection, RTM works by running the wave equation backward for all transducers. It sets up all transducers as pseudo sources, and the records are broadcasted with a time-reversed manner-from the end of the trace to time zero. If the velocity model is presumed correctly, in TAT imaging the wave field at time 0 in space will coincide with the original source distribution [19]. Compared with back-projection, by applying full wave equation, RTM is capable of involving all wave field phenomenon including diffraction, aberration and multiple scattering. However it requires much higher usage of CPU time and memory [27].

Validation via synthetic data
To test the effectiveness of back-projection and RTM, we have built a 2D synthetic human brain model with an intact skull. In this model, the brain is made of gray matter and white matter [14], and the skull is made of three layers, namely the inner table, diploe, and outer table [20]. To mimic a real laboratory experiment similar to [8], [13], we modeled the space outside the skull as mineral oil. The reason for us using mineral oil is that it provides both back-ground with uniform acoustic speed and tiny loss coupling between acoustic transducers and human head surface. Additionally when compared with water, mineral oil has much weaker microwave absorption ratio (~0), this will guarantee the minimum errors are introduced during modeling. Instead of mineral oil, a plenty kinds of fluid with weak microwave absorption ratio and comparable acoustic velocity to bio-tissue could also been used in further realistic test. To mimic pathological changes and build a benchmark for results analysis, in the synthetic model we replaced a small area of the brain with blood. This area was located at the left cerebral hemisphere and defined as elliptically-shaped. The distribution of acoustic velocity in the synthetic model is shown in Figure 3a. From a review of several literatures, mechanical parameters of all related bio-tissues were collected and are listed in Table I, in which the velocities and densities of grey matter and white matter were measured from lamb brain using acoustic frequency of 1 MHz [14]. The skull's velocities, densities and thicknesses of different layers are applied from datasets in research [20], [34] and [35]. The loss factor is defined in [35] to depict the amplitude of the propagating acoustic wave's energy decay. Its values come from [14] for brain and [35] for skull, respectively.
After establishing the mechanical properties, we calculated the initial acoustic amplitudes in the synthetic model. We assumed that the initial acoustic pressure (dyne/cm2) in mineral oil and the skull was 0, so the acoustic wave field was entirely generated by multiple acoustic sources in the brain at time zero. Their amplitudes were closely linked to the microwave absorption ratio. As expressed in [6], the relationship between the power intensity I (W/cm2) of absorbed microwaves and the generated peak acoustic pressure P0 (dyne/cm2) is shown as Eq. 1, where c is the sound velocity, β is the volumetric thermo-expansion coefficient, and C p is heat capacity. The definition of absorbed power intensity I can be expressed as Eq. 2 [8], where E, ρ, and σ are the maximum amplitude of the radiated electromagnetic field, tissue's density, and electrical conductivity, respectively. By combining Eq. 1 and Eq. 2, Eq. 3 is derived as the theoretical relationship between tissue's electrical conductivity and initial acoustic pressure generated based on the thermo-acoustic effect.
The volumetric thermo-expansion coefficient β and heat capacity Cp of the brain are nearly uniform: (β=12.3X10 -5 /℃ [36] and C p =4160 KJ/m 3 ℃ [37]). The acoustic velocity and density vary little among gray matter, white matter, and blood (Table 1). Consequently the initial acoustic pressure can be well approximated as proportional to the electrical conductivity σ.
On the other hand, Table 2 shows that for input microwaves the electrical conductivities of white matter, gray matter, and blood are significantly different [8]. For blood in particular, σcan be more than 1.5 times higher than for the other two kinds of tissue. Using the value of electrical conductivity for microwave of 900 MHz from [8] and Eq. 3, the distribution of initial acoustic pressure is derived as shown in Figure. 3b.    Table 2. Selected electrical conductivities of human brain Once the synthetic model was established, we applied the FDTD method to forward modeling acoustic wave propagation. Using the initial acoustic pressure shown in Figure. 3b, a zerooffset Ricker wavelet with a central frequency of 0.15 MHz was applied to each point of the brain as the acoustic source. The model space was meshed as 512 x 512 grids with a spatial interval of 0.5 mm. A total of 1600 time steps with time intervals of 0.115 µs were judged long enough to allow the acoustic waves to propagate to each receiver from the most remote grid in the brain. The outgoing acoustic signal was recorded by 240 receivers located outside the skull (white circle shown in Figure. 3). The derived synthetic time traces are shown in Figure. 4, which provides the input dataset for later imaging.
The velocity model used in both migration algorithms is critical to successful imaging. In this study we applied two velocity models: one (abbreviated as V1) assumes the average acoustic velocity is uniformly 1540 m/s in the model space. Essentially it approximates a "bare brain" model with the effect of the skull excluded. The second model (abbreviated as V2) includes the effect of the high acoustic speed of the skull, which is almost double the speed of brain tissues. In our study, velocity models V1 and V2 were applied to both back-projection and RTM. Due to the velocity variance in V2, we applied V2 to two migration methods by different approaches. In back-projection, ray-tracing was applied from every transducer to all directions. This procedure was similar to the methods described in [13], but for simplicity we considered only the rays' travel time caused by velocity variance, and aberration around the skull was ignored. Different from back-projection, as a kind of full-wave migration, RTM uses the same scheme as forward modeling methods such as FDTD. Consequently V2 can be applied to RTM in a straightforward manner. Comparisons of migration imaging results are shown in Figure. 5-7.    Figure 5a andFigure 6a are results derived by back-projection and RTM using the "bare brain" model V1. When compared with the distribution of initial acoustic pressure of the original model ( Figure 5b or Figure 6b) we can clearly observe two kinds of imaging artifacts. First, the area with higher initial pressure at the left cerebral hemisphere, which is set up artificially in an elliptical-shape, is seriously enlarged by back-projection ( Figure 5a) and falsely elongated along the major axis of the ellipse by RTM (Figure 6a); Second, delicate features such as the gyrus, located in the outer part of the brain, and the gap which separates the left and right cerebral hemispheres in our model are totally blurred in both results when using velocity model V1. In contrast, from Figure 5c and Figure 6c, which are results based on velocity model V2 with the skull's velocity included, these two misfits are substantially reduced. From all of these comparisons we can see that exclusion of the skull leads to severe error and distortion in migration imaging for both back-projection and RTM.
To further examine the differences between back-projection and RTM, we reorganized the results using back-projection, RTM, along with the original velocity model V2 as shown in Figure 7. From it we can see that although both back-projection and RTM can transform most of the wave field back to its original location correctly, there are obvious differences in imaging quality between the two methods. Compared with the original model shown in Fig. 5b, the detail features are blurred in back-projection result (Figure 7a) but appear to be clean and sharp in RTM result (Figure 7c). These visual differences can be further amplified through 1-D comparison along the x-direction along the horizontal white dashed line shown in Figure 7 chosen to cross the brain's left boundary, artificially blooded area, interhemispheric fissure, and right boundary. The source amplitudes along this profile for back-projection RTM, and the original model are shown in Figure 8. It is obvious that the result of RTM is far more superior than that of back-projection. Compared with the back-projection result (the dotted line), result of RTM (the solid line) has larger variance for depicting structures such as interhemispheric fissure around 12.50 cm as well as the boundaries of the artificially blooded area around 7.5 cm and 10 cm. These differences can also be clearly observed on brain boundaries around 5 cm and 20 cm, where the amplitude from RTM decays as sharp as the original model, but back-projection's result is obviously incorrect and decay much slower outside of the brain.  Fig.5a-c, which passes brain's left boundary, artificially blooded area, interhemispheric fissure and the right boundary. The dotted-broken line, dotted line, and solid line show the original model (Fig. 5b), the results by back-projection (Fig. 5a), and RTM (Fig. 5c) using the model V2.

Validation through laboratory data
The back-projection and RTM algorithms were also tested by using the laboratory data acquired by XW06. In their experiment, the monkey's head was decapitated and fixed by a clamp and completely immersed in mineral oil. During TAT detection, this specimen was stimulated by 3-GHz microwave pulses, and the derived acoustic wave field was recorded by a transducer with a 1 MHz central frequency and about 0.8 MHz bandwidth. The transducer was positioned from 6-14 cm to the center of the monkey's head, and the sampling frequency was 20 MHz. During the experiment, the clamp fixing the monkey's head was mounted on a rotary table driven by a stepper motor with a step size of 2.25 degrees. Accordingly in this laboratory application, the outgoing acoustic wave was observed by 160 receivers surrounding the head in a 2D circle. With data processing performed through the procedure of [38] for high frequency enhancement, only the segment with a spectrum of 0.3-1 MHz of the observed data was picked up and enhanced for imaging by back-projection and RTM. We applied the estimated average acoustic velocity to both image approaches, since any velocity information on velocity distribution of earlier experiment in XW06 was unknown, which may introduce some error and reduce the image quality. Figure 9 shows the results based on a dataset collected from a one-month-old monkey head with a skull thickness of less than 1 mm. The velocity model is assumed to be uniform as 1437 m/s, as the same as acoustic velocity of mineral oil. The region shown is 53mm by 51mm along the coronal cross section. From the experiment of XW06, three steel needles with diameters of 0.9 mm were inserted in the approximate locations as shown in Figure. 9a (XW06). The results derived from back-projection ( Figure. 9b) and RTM ( Figure. 9c) are shown side-by-side for comparison. Both imaging algorithms show the three needles, and the black dot located at the center is believed to be an air bubble introduced by inserting the needles (XW06). Compared with the result derived from back-projection, the needle A in RTM result owns sharper edges. Meanwhile, back-projection provides less visibility for needle C than RTM. From the plots along the x cross-section shown in Figure. 9d, although both needle A and B can be detected by using back-projection and RTM, the back-projection image is much nosier. The existence of this noise causes seriously reduction of signal noise ratio and fussy in whole image by backprojection. However needle B is seriously distorted from results of both back-projection and RTM. This shows that due to the technical limitations of the coarsely estimated average acoustic velocity for TAT reconstruction, neither of these two methods can provide satisfying image quality.

Iterative RTM procedure
The previous result shows that the image quality of RTM is restricted by the precision of the velocity model. In some sense, the increased requirement of RTM sensitivity to indistinct velocity model brings a paradox: If the velocity model is well known by operation or any invaded techniques, RTM is not necessary; If the velocity model is not known, RTM lacks an essential input. For solving this paradox and improving the RTM result, rather than single RTM run, we developed iterative RTM procedure based on the theory of Whitmore [23]. The iterative RTM procedure works based on an assumption between velocity model and RTM result. Suppose there exists a function to map velocity from RTM result, complex velocity model could be renewed after RTM for one time. By repeating this procedure, velocity model will continue be updated after multiple RTM runs until the result becomes promising. This iterative RTM procedure is shown in Figure 10 as a schematic. An initial guess of a model is made, incorporating all known external velocity information. The image is derived through RTM. After that a new velocity model is made based on current RTM result. This procedure is repeated until both RTM results and velocity models are unchangeable.
As Figure 10 shows, the key procedure of iterative RTM is to update velocity model after each RTM. This requires a function for mapping velocity from image derived by RTM. Since for TAT imaging RTM image is the reconstruction of electrical conductivity distribution, the desired function of velocity and electrical conductivity could be built by referring to Table 1 and Table 2. Based on the given velocity and electrical conductivity values of white matter and gray matter, a function is established by using the second order polynomial fitting as shown in Figure 11. Obviously this function is coarsely approximated. It could be improved by involving more velocity-electrical conductivity pairs in further research. Meanwhile it's noteworthy that the mapping function can't derive correct velocity on blood, water and mineral oil, consequently a mask must be used to preserve the region outside sample from wrongly updating. At present the function shown in Figure 11 is used in our research to map the velocity from electrical conductivity derived by RTM; in this way the velocity model will be updated for the next RTM iteration.  Mapping function for calculating acoustic velocity through RTM derived electrical conductivity. This function is second-order polynomial fitted by velocity-electrical conductivity pairs of white matter and grey matter given by Table 1 and Table 2.
An example of applying iterative RTM procedure is shown in Figure 12. The laboratory data of this case comes from XW06, which has been introduced in Section V. By referring to XW06, the velocity model would not be known during TAT. The only knowledge of the velocity field in this case is 1437 m/s of mineral oil, which submerges the sample during TAT detection. To demonstrate the iterative RTM procedure, we make an initial guess that velocity model is uniform as 1437 m/s. By using the initial velocity model, the first RTM output is derived as shown in Figure 12b, which is as the same as Figure 9c. Velocity model will then be updated as new velocity input for RTM. This procedure is repeated for several times until RTM output is unchangeable. The TAT results derived from iterative RTM procedure on step 3 and 5 are displayed on Figure 12b and Figure 12c.   Figure 12, by using iterative RTM procedure the improvement on TAT image can be observed from several aspects: First, by comparing with Figure 12b (Figure 9c), the distortion of needle B has been well corrected after five times iteration. Second, the original RTM result shown in Figure 12b provides less visibility for needle C. After two iterations, Figure 12c shows that the boundary of needle C has been largely enhanced. Additionally, comparing among Figure 12 a-c, the cross-section of needle A has been focused gradually with iteration times increasing. All of this shows that iterative RTM procedure has great potential of correcting coarsely estimated velocity model and therefore enhance imaging capability of RTM.

Discussion
By combining Section II and III, It is clear that RTM is superior to back-projection in terms of imaging quality and higher signal to noise ratio. Compared with back-projection, which makes a ray approximation, RTM bases its entire algorithm on solving a full-wave equation, without substantial approximation, and holds the original dynamic features of the wave field intact. The handle of velocity heterogeneity is fundamentally intrinsic. Usually, adapting ray tracing in back-projection is time consuming and has a limited improvement on image quality. Unlike back-projection, the quality of RTM's results is independent of the complexity of the velocity model. This feature makes the wave propagation in the domain highly accurate compared with using ray-tracing. Figure 5-8 show that RTM is able to recover almost all features of the brain to their original potion when the skull's velocity is included.
RTM can recover a complex structure's boundary sharply. This has been reported by [19], [27] and is proved by our results in Figure 6. When both velocity models V1 and V2 are used in RTM, the boundaries of tiny features can be clearly seen. Especially, even though obvious distortions exist in the result using V1 (Figure 6a), with the exclusion of the skull, all features are still relatively un-blurred in comparison with the back-projection results (Figure 5a). By comparison, looking at the cross-section in Figure 8, the sharp edges of brain are well recovered by RTM but seriously smeared by back-projection. Further, in Figure 9 when the skullexcluded model is applied, the contour of Needle C is well recovered by RTM but not by backprojection.
Nevertheless, it is noteworthy that the image quality of RTM is still limited by the precision of the velocity model. Consequently the key to capitalizing on the benefit of RTM is to build better velocity models before applying RTM. For achieving this goal, we developed iterative RTM procedure to update velocity model iteratively. Its principle is based on an assumption that acoustic velocity could be mapped from RTM image through a certain function. The results shown in Figure 12 demonstrated that currently even by using our coarsely approximated velocity function, the quality of RTM image can be well improved after several iterations within iterative RTM procedure. It could be improved by involving more velocity-electrical conductivity pairs in our further research.

Conclusion
In this paper we have compared back-projection and RTM for transcranial TAT imaging. Compared with back-projection, RTM offers better performance with regard to velocity variance, imaging quality, and noise suppression caused by spatial aliasing. The capability of RTM can be further improved by iterative RTM procedure, which aimed to provide velocity model with higher accuracy. Generally speaking RTM owns great potential for achieving acoustic localization on transcranial TAT imaging with high quality and accuracy.

Author details
Zijian Liu * and Lanbo Liu *Address all correspondence to: lzjknight34305@gmail.com School of Engineering, University of Connecticut, Storrs, CT, USA