The Levenberg–Marquardt Method for Acousto-Electric Tomography on Different Conductivity Contrast

: The stability and convergence performance of Levenberg–Marquardt method for acousto-electric tomography (AET) applied to different levels of conductivity contrast is studied in this paper. As a multi-physical imaging modality, acousto-electric tomography (AET) provides high spatial imaging resolution while also conserving the high contrast property of electrical impedance tomography. The Levenberg–Marquardt method is a well known iteration scheme which can be applied for the nonlinear problem of AET. However, the inﬂuence of the conductivity contrast on the stability and convergence performances of this conductivity reconstruction method is rarely discussed in the literature. In this paper, the performance of the Tikhonov regularization-based Levenberg–Marquardt method is applied to reconstruct conductivity map with different conductivity contrast between different regions of the domain of interest (DOI). Numerical investigations are carried out for phantoms with different conductivity contrast. Reconstructed results with different levels of noise are presented and discussed in detail.


Introduction
Electrical Impedance Tomography (EIT) is a noninvasive technique for reconstructing the conductivity distribution from the measured potential on the electrodes that are attached to the surface of the domain of interest (DOI), and current is injected through the same electrodes into DOI [1]. As the total number of measurements is much smaller than the total number of points to be reconstructed, the EIT problem is severely ill-posed. Acousto-electric tomography [2][3][4][5] is developed to improve the ill-posedness and spatial imaging resolution of electrical impedance tomography (EIT) while conserving the high contrast property of EIT. Other methods, e.g., learning structural sparsity via a Bayesian approach [6,7], and some based on multi-physics models [8,9], exist as well for the same purpose. AET is applied to determine the internal conductivity of a physical body with better stability and resolution compared to the ill-posed EIT. The main idea is to carry out a classical EIT measurement while a known focused ultrasonic wave travels through the DOI. The high intensity of the pressure produced by the acoustic wave creates a small local deformation within the DOI, and this deformation changes the electrical conductivity locally, and therefore the electrical current distribution, which further changes the electric voltages measured on the electrodes. The higher the electrical conductivity in the acoustic wave focal zone, the greater the change. This deformation also influences the resolution of the imaging [2]. The resulting changes of the boundary potential can be recorded with EIT measurements [10], which are then used to compute the power density [2,11] for reconstructing

The Inverse Problem and Forward Model
A physical body imaged by EIT is modeled as a bounded Lipschitz domain Ω ⊂ R n for n ≥ 2. The changes caused by the acoustic wave can be recorded with EIT measurements on boundary ∂Ω [18][19][20]. The power density in Ω is defined as where σ is the conductivity map in domain Ω. u(σ) is the electrical potential produced by applying either voltage or electric current on ∂Ω. Given noisy measurements E δ (σ) of the true power density E (σ), the problem is to reconstruct conductivity map σ via minimizing the following functional, The power density E (σ) can be reconstructed from the boundary measurements [2,21]. The nonlinear relationship between E (σ) and σ renders this problem nonlinear [22]. To develop the methods for reconstructing σ from E (σ), the continuum model with Dirichlet condition is used here f is the boundary potential, f ∈ H 1/2 (Ω). The present work assumes σ is positive. With the assumption that measurements are carried out in the low frequency domain without any interior sources or charges in Ω, Equation (3a) is obtained by simplifying time-harmonic Maxwell's equation with low-frequency approximation [17].

The Reconstruction Algorithm
The Levenberg-Marquardt algorithm (LMA) locates the minimum of a function expressed as the sum of squares of errors between function and measured data through iterative updating. For a general operator equation F(σ) = E (σ) δ , where F : X → Y and X, Y are Hilbert spaces [23]. E (σ) δ is the measured noisy data of E (σ). More practically, one could instead minimize [24,25] letting α k the Tikhonov regularization parameter, and δσ = σ − σ k . This minimization problem regarding to (4) is equivalent to the linearized one letting F (σ k ) * the adjoint of F (σ k ). It has been proven in [17,[26][27][28][29], with a properly chosen parameter α k and an initial guess σ 0 sufficiently close to the desired solution, that the LMA converges to a solution The LMA based on Equation (5) can be thought of as a combination of steepest descent and Gauss-Newton method. When the current solution is far from the correct one, a large value is given to α k , and LMA behaves like a steepest-descent method, which converges slowly but well. When the current solution is close to the correct solution, a small α k is used, and LMA behaves like a Gauss-Newton method which is known as a faster iterative method than the steepest descent. The Fréchet derivative and adjoint operator needed in (5) are computed with a method proposed in [17], as follows. It is assumed that σ is positive, and σ ∈ H l (Ω) for l > n/2. For f ∈ H l+ 1 2 (Ω), then the forward model has a unique solution u(σ) ∈ H l+1 (Ω), which indicates the operators E (σ), E (σ), and E (σ) * are all maps from H s (Ω) to H s (Ω) [17].
The boundary conditions for u (σ)τ can also be derived from Equation (6). If ξ = u (σ)τ, then ξ satisfies Here, it is assumed that τ vanishes in a domain Ω = {x ∈Ω|d(x, ∂Ω) < } for a sufficiently small . Therefore, all terms regarding to τ in the boundary conditions disappear. It is well known that the solution of Equations (8) uniquely exists in a weak sense. u (σ) as defined in (8) is the Fréchet derivative of u(σ).
The directional derivative of E (σ) is given as with terms of the order of τ 2 neglected. E (σ)τ is actually the Fréchet derivative of E (σ).
β is a regularization parameter, which enforces a second-order weak derivative. Upon integration by parts, B * ζ is the solution of the Neumann problem with ς ∈ L 2 (Ω). Refer to the work in [17] for more detail. This fourth-order partial differential equation is not practical, but it can be written into an equivalent form According to the Levenberg-Marquardt iteration in Equation (5), the formula for calculating the k-th updating step τ k for the problem at hand is given as where E δ represents the measured power density derived from u b and u b, in the last section. For the left hand side of (18), E (σ k ) * E (σ k )τ = B * Mτ with Mτ =|∇u| 2 (τ|∇u| 2 + 2σ∇u∇u (σ)τ) + 2∇u∇V(τ|∇u| 2 ) Here, Mτ is easily obtained with Equations (9) and (11). After computing τ k from (18), the conductivity map σ k obtained at the k-th iteration is updated by σ k+1 = σ k + τ k for a new iteration. To calculate B * Mτ, Equation (17) needs to be solved with ς = Mτ, which requires first solving Equations (8) and (12). As all these PDE systems are coupled to each other via τ, they need to be collected into one PDE system and solved in one variational form.
It yields the following system, Equations (20g)-(20l) must be solved for each measurement, so one additional measurement will require to solve in addition three partial differential equations. The smoothness requirement on σ can be achieved by a so-called mollification operation. To solve the above system, Thus, w can be calculated with (12), and x = B * z is the solution of (17).
The iterative construction method based on the above system is named as LM-CM, referring to Algorithm 1 for a single measurement. As ν(x) is not known for the current numerical study, the data of E δ are simulated by solving (3a) with a Dirichlet boundary condition u = u b in the present contribution. The relative error η is defined as η = σ t − σ r L 2 / σ t L 2 . Index t indicates the true value, and r means retrieved value. The parameter α k should be updated according to the value of τ k . If σ k + τ k leads to a reduction of the relative error in σ k , α k is decreased and τ k is accepted. Otherwise, τ k is discarded and α k is increased. As η cannot be determined in practice, a relatively large value can be given to α 0 , and α k is slowly decreased to ensure convergence. The iteration is stopped when τ k L 2 is smaller than the expected error or the maximum number of iterations is reached. Algorithm 1: LM-CM algorithm for retrieving the conductivity map of a domain Ω from a single measurement of power density.

Numerical Investigations
The conductivity distribution can be obtained from the power density E (σ) with Algorithm 1 based on the Levenberg-Marquardt method. In the following computation, E (σ) is computed directly from its definition σ|∇u| 2 , where u is obtained with the forward model introduced in Section 2.
Different phantoms with different contrast between the inclusion and the background are studied here. The weak forms of the equations in (20) can be easily obtained by applying integration by parts, and the equations are solved with finite element method implemented by using the open-source finite element library FEniCS [31]. Two boundary conditions with f 1 (x, y) = x and f 2 (x, y) = y are given to Equation (3a) for a reconstruction with multiple measurements. A relative error η is defined as η = σ t − σ r L 2 / σ t L 2 to indicate the convergence of the Levenberg-Marquardt iteration. Index t indicates the true value, and index r means retrieved value.
Practical measurements are always affected by various noises introduced by the devices or during the measurement process. Here the performance of the algorithm to different levels of contrast is investigated with different levels of noise. The noise level is measured with signal-to-noise ratio (SNR), which is defined as where N satisfies Gaussian white noise distribution. The noise is added to the data of E (σ) for numerical examples with SNR = ∞, 40 dB, and 20 dB. SNR = ∞ for noise free.

Numerical Study with a Simple Phantom
A phantom with a smaller circular inclusion included in a circular domain is used here for numerically studying the performance of the Levenberg-Marquardt method to various conductivity contrast. The radius of the circular DOI is 1.0 m. The circular inclusion with its center located at (−0.3, 0.3) has a radius of 0.2 m. The conductivities of the background and the inclusion are denoted as σ b and σ i , respectively. The contrast of the conductivity is defined as σ i /σ b . α k for the k-th iteration is set as α k = α 0 /a k with α 0 = 1.0, a = 2.0, and β k set as constant β 0 = 1e − 3.
The mesh for the forward model is shown in Figure 1a which has 3854 vertices and 7522 cells. When α 0 = 0.1, the results given in Figure 2a show that the method can converge fast for both low and high contrast cases. However, the relative errors for the high contrast cases, e.g., σ = 12.0 S/m, 14.0 S/m, start to diverge after achieving to the minimum relative error. This divergence becomes severe when the noise level is increased, as shown in Figure 2b. When the value of α 0 is increased, the divergence phenomenon disappear for both lower and higher noisy cases, but the speed of convergence slows down, and it takes much more iteration steps for the computations to achieve a similar level of relative error as α 0 = 0.1. A direct comparison of the results is given in Figure 2a,b. A larger value of α 0 brings higher stability even with higher noise. Therefore, a smaller α 0 will accelerate the convergence of the computation, and a larger α 0 will improve the noise tolerance and the performance of the method for high contrast conductivity distribution. To accelerate the convergence without losing the stability of the computation, α k for the k-th iteration is set as α k = α 0 /a k with α 0 = 1.0, a = 2.0.
The values of β k = β 0 are fixed in the following computations. However, a different value of β 0 is possible to provide different performance of the Levenberg-Marquardt method for high contrast cases. The model with two different levels of noise (SNR = 40 dB, 20 dB) will be considered when conductivity is 10.0 S/m, 12.0 S/m, and 14.0 S/m. The β 0 is set as 0.0001, 0.001, and 0.01. The results are shown in Figure 3a,b. When the value of β is increased from 0.0001 to 0.001, the speed of convergence is similar, this indicates the convergence is mainly dominated by the value of α k . However, the minimum relative error achieved during the iteration process is decreased. For the high contrast case, the instability is also influenced by the value of the β k , a larger β k with a better stability, and this becomes more obvious when more noise is introduced into the measured data, as shown in Figure 3b. Therefore, a smaller β 0 will make the minimum relative error smaller during computation, however it decreases the noise tolerance. Therefore, the value of β k is set to β 0 = 1e − 3 in the following computations. The values of β 0 and α 0 could be greater or smaller whether needed. Another phenomenon observed in Figure 3a,b is that when α k becomes very small, the relative error for the high contrast cases starts to increase, and it stops at certain level. This is the instability caused by the high contrast of the conductivity. Therefore, the value of α k can not be too small. In what follows, the variation of α k follows α k = α 0 /a k at the beginning of the iteration, but it is fixed to the value where the minimum relative error is achieved.

Numerical results with SNR = ∞
When SNR = ∞, the relative error corresponding to each conductivity is shown in Figure 4. Here, α k for the k-th iteration is set as α k = α 0 /a k with α 0 = 1.0 and a = 2.0, and β k is set as constant β 0 = 1e − 3. The figure shows that the algorithm works stably with good convergence when the conductivity σ i is less than 10.0 S/m, where the contrast is 5. When the contrast is more than 5, the algorithm converges to certain level, and then it starts to diverge. The curve of the relative error varying with iteration steps becomes a "V-shape". The minimum relative error becomes larger as the contrast increases, as obviously seen in Figure 4. This phenomenon is also observed in Figure 3a for high contrast cases and Figure 3b for almost all the noisy cases. Theoretically, the conductivity distribution should be in the H 2 (Ω) space, and the value of β k is closely related to the second order derivative of the conductivity distribution as shown in Equation (15). When the contrast is too high, this continuity requirement can not be satisfied, this "V-shape" curve of the relative error appears. When the conductivity of the inclusion is σ i = 4.0 S/m, the minimum value of the relative error η = 0.8%, and this value is achieved at the 12th iteration step. The reconstructed conductivity and the power density are shown in Figure 5a,b, respectively. A very good level of reconstruction is achieved. The change of the power density across the boundary between the inclusion and the background is not large. When σ i is further increased to σ i = 12.0 S/m, the minimum relative error is achieved at the 7th iteration step and η = 0.120. The reconstructed conductivity is given in Figure 6a. However, η starts to increase as the iteration proceeding. The reason of the increased error can be explored from power density, and the one at the 7th step is given in Figure 6b. Comparing to Figure 5b, there is a larger jump of the power density across the boundary between the inclusion and the background. This may cause the relative error increase and make the algorithm become unstable and poorly convergent.

Numerical Results with Noise
When SNR = 40 dB, 1% of noise is added to the computed power density for the reconstruction. Here, α k for the k-th iteration is set as α k = α 0 /a k with α 0 = 1.0 and a = 2.0, and β k is set as constant β 0 = 1 × 10 −3 . The variation of the relative errors of reconstruction corresponding to different conductivities is shown in Figure 7. The algorithm works fine with σ i < 10 S/m though the relative error increases even with such a low level noise. If the contrast is continuously increased, the instability of the algorithm becomes severe, and the minimum relative error also becomes larger.
The minimum relative error η = 0.011 is achieved at the 11th step when the conductivity of the interior circular domain is σ i = 4.0 S/m. The reconstructed conductivity distribution is shown in Figure 8a, and the corresponding power density is shown in Figure 8b. The reconstructed conductivity map is attractive although noise exists in the DOI, and the change of the power density across the boundary is not much. When the conductivity is further increased to σ i = 12.0 S/m, the computation starts to diverge after the 7th step, and the relative error begins to increase until to η = 0.142. The reconstructed conductivity map and the power density of the 7th step are given in Figure 9a However, when the noise level is further increased to SNR = 20 dB, the performance of the algorithm becomes more sensitive to the contrast. The variation of the relative errors corresponding to different conductivities is shown in Figure 10. The convergence of the computation is poor even for σ i = 8 S/m. For higher contrast, the computation converges well in the first several steps, but it starts to diverge afterwards. The minimum relative error η = 0.050 S/m is achieved at the 5th step when the conductivity of interior domain is 4.0 S/m. The reconstructed conductivity map is shown in Figure 11a, and the corresponding power density is shown in Figure 11b. The results are obviously contaminated by noise, and fluctuation is observed in the DOI. This is caused by the large noise introduced into the power density as shown in Figure 11b. To conclude, the performance of the introduced Levenberg-Marquardt method for acousto-electric tomography has a close relation with the contrast of the conductivity. It works adequately for a problem with low contrast and low noise level. This is observed when the contrast is smaller than 5. Then, the performance becomes worse with the increase of the contrast, and the minimum relative error also increases. However, higher noise will increase the sensitivity of the method to the contrast, and therefore decrease the stability of the algorithm and the quality of the reconstruction. These performances are also closely related to the smoothness of the distribution of power density in the domain of interest. However, it should not be neglected that the parameters α k and β k introduced in (20) also play an important role.

Heart Lung Model
To further investigate the performance of the algorithm, a heart-lung model is created and used here, which is shown in Figure 12a. The considered tissues are heart (red, σ = 0.7 S/m), lung (cyan, σ = 0.26 S/m), and soft tissues (blue, σ = 0.33 S/m). Refer to the work in [16] for conductivities of different tissues. An additional circular domain is added to the left lung to simulate a lung cancer, and the tissue is in purple color. The conductivity of lung cancer tissue is generally 1.6 to 3.3 times larger than the normal tissues; here, it is set to σ = 0.50 S/m, which is~1.92 times of the lung tissue. The model is placed into a circular region with a background material σ = 0.22 S/m and the radius is 0.25 m. The true conductivity distributions are shown in Figure 12b. The inverse mesh consists of 12,539 vertices and 24,730 cells.The initial guess for the interaction is set to 0.22 S/m for the whole DOI. Here, α k for the k-th iteration is set as α k = α 0 /a k with α 0 = 1.0, a = 2.0, and β k is set as constant β 0 = 1 × 10 −3 . The conductivity is reconstructed with the method introduced in this paper. The relative error and reconstructed conductivity map with SNR = ∞, 40 dB, and 20 dB are shown in Figure 13a-c, respectively. In the computation, each iteration takes about 16 s and the total 50 iterations takes about 13 min. For the noise free case and the case with very low noise level as SNR = 40 dB, the minimum relative error is~0.074. The convergence of the computation is fast, but the reconstructed results are not good, as can be seen in Figure 13a,b. When higher noise level is considered, the convergence of the computation becomes worse, and the reconstructed results are also not good. As seen in Figure 13c for SNR = 20 dB, although the general conductivity contribution is observed, and a minimum relative error 0.094 is achieved, a lot of noise-like areas appear in the whole DOI. However, different tissues can still be clearly located in the reconstructed conductivity distribution, although their shapes are not so clear. To investigate the influences of the conductivity jump on the boundaries of two subdomains in the DOI, a smooth-heart-lung model is obtained with a mollification operation. The piece-wised conductivity distribution is mollified with to produce σ ∈ H s (Ω) with s = 2 for 2D problem. The constant C > 0 is selected so that R 2 η = 1. The value of is 1 cm for heart-lung model. The true smoothed distribution of σ for the heart-lung model is shown in Figure 14. The numerical setup for the reconstruction of smoothed heart-lung model is kept. The relative error and reconstructed conductivity map of smooth-heart-lung model with SNR = ∞, 40 dB, 20 dB are shown in Figure 15. The minimum relative errors are~0.002 when SNR = ∞ and 0.007 when SNR = 40 dB, which is much smaller than the case with a piece-wised conductivity distribution. Even for the more noisy case with SNR = 20 dB, the minimum relative error is about 0.037, which is also much smaller. The convergence of the computation also performs much better. This significant decrease of relative error and increase of convergence speed are mainly caused by the smoothness of the power density which is further resulted from the smoothed conductivity distribution. The power density of corresponding to the minimum relative errors for SNR = ∞, 40 dB, 20 dB are given in Figure 16 for both piece-wised and smoothed conductivity distributions of heart-lung model, respectively. A comparison shows that the power density for the smoothed conductivity is smoother, and there is no obvious discontinuity across the boundary of different subdomains. Therefore, a rather smooth distribution of the conductivity can lead to a smoothed power density distribution, and this can well improve the stability of the algorithm and accelerate the convergence as well.

Conclusions
The Levenberg-Marquardt method for the acousto-electric tomography is applied to reconstruct the conductivity distribution in the DOI from voltages measured on its boundary. The main focus has been on the influence of the contrast on the stability and convergence performance of the Levenberg-Marquardt method. It also discussed the influences of the regularization parameters. Several numerical investigations have been carried out for the illustration. The work herein, on the one hand, focuses on the higher contrast potential reconstruction method, and on the other hand, on validations of its practicability on detecting lung cancer.
It has been discovered that the algorithm becomes unstable when the contrast is approximately beyond 5. A contrast beyond this level will bring larger minimum relative error and a divergence after several iterations. The regularization parameter α k plays an important role in the computation. A larger value of the α k can make the computation converge better, but it will also render the convergence slower. A smaller α k will bring a faster convergence as long as the solution is not too far from the true value, otherwise the algorithm will be unstable. Therefore, a strategy of decreasing α k with the progress of the iteration can accelerate the convergence and keep the computation stable.
Theoretically, the conductivity distribution should be in the H 2 (Ω) space, and the value of β k is closely related to the second order derivative of the conductivity distribution. When the contrast is too high, this continuity requirement can not be satisfied, a "V-shape" curve of the relative error appears. Therefore, the value of β k is also important for the stability and accuracy of the reconstruction. A large β k will enforce strong continuity requirement on the conductivity distribution, which will render the computation unstable. A small β k releases the continuity requirement, which will produce a poor reconstructed results. Therefore, the value of β k needs to be properly chosen. It is well known that the Tikhonov regularization can smooth the reconstructed conductivity, therefore a large α k will balance the influence of a small β k .
Besides the above parameters, a high contrast between the conductivity of different regions will cause a discontinuity of power density across the boundaries; this will not meet the continuity requirement in the reconstruction algorithm, which will also bring the instability and poor convergence, as demonstrated by the numerical results of the high noisy and high contrast cases.
The reconstructed results of a heart-lung model and a smooth-heart-lung model are compared to support the above conclusions, and a much better convergence and stability are observed with this smoothed heart-lung model. However, it has been shown that the Levenberg-Marquardt method works quite well in general, and performs much better in low contrast case. The instability caused by the conductivity contrast can be improved with a proper Tikhonov regularization parameter. A simple dynamic strategy is used in this paper for choosing the proper value of Tikhonove regularization parameter α k , which requires the ground truth of the conductivity distribution. Although it works quite well in all the numerical results presented in the paper, it is not practical, as the ground truth is generally unknown in most cases. However, one can use the Morozov's discrepancy principle for choosing this proper value, which should work better in theory.
Author Contributions: Conceptualization, C.L. and K.Z.; methodology, C.L.; software, C.L. and K.A.; formal analysis, C.L. and K.A. investigation, C.L. and K.A.; writing-original draft preparation, C.L. and K.A.; writing-review and editing, K.Z.; All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflicts of interest.