A Complete Framework for Acousto-Electric Tomography With Numerical Examples

Acousto-electric tomography (AET) involves three steps to retrieve the distribution of conductivity in a domain of interest (DOI): measure the potential on the DOI boundary (usually with a limited number of electrodes set there), compute the power density from this potential, use it to retrieve the distribution sought after. Almost all developed algorithms for AET assume that the power density is known, so their focus is mostly on the last step. A complete framework for AET is proposed herein to connect the three steps, the complete electrode model (CEM) being used to simulate the voltages measured on electrodes. The potential on the whole DOI boundary is reconstructed from such voltages. Then, the power density is computed, and the conductivity distribution in the DOI retrieved. A method based on singular value decomposition (SVD) is proposed. This method and the iterative Levenberg-Marquardt method are used for numerical illustration. The SVD-based method yields the potential on the whole DOI boundary, and a gross map of the conductivity distribution is also obtained, to serve as initial guess of the Levenberg-Marquardt method to yield the conductivity contribution with higher accuracy.


I. INTRODUCTION
Electrical impedance tomography (EIT) has attractive applications in medical imaging [1]- [3]. In the measurement, a number of electrodes are placed at the boundary of the domain of interest (DOI) and a given electrical current is injected between a pair of electrodes. The potentials on the other electrodes are then measured. It goes on likewise for every pair, and the corresponding measurements then yield a whole set of data which can be used to produce an image of the electrical conductivity distribution within the DOI, two-dimensional at least if proper symmetries or simplifications, three-dimensional at best. However, the total number of measurements is still much smaller than the number of points to reconstruct, thus EIT is an ill-posed inverse problem, and the spatial resolution remains low compared to other medical imaging methods.
The associate editor coordinating the review of this manuscript and approving it for publication was Ge Wang . Some hybrid imaging methods [2] have been developed to overcome ill-posedness of electrical impedance tomography (EIT) and to correspondingly improve the image quality [3]. Other methods, e.g., learning structural sparsity via a Bayesian approach [4], [5] and some based on multi-physics models [6], [7] exist as well. 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 locally changes the electrical conductivity, hence 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 larger the change. This deformation also influences the resolution of the imaging [1]. The resulting changes of the boundary potential can be recorded with EIT measurements [8], which are then used to compute the power density [1], [9] for conductivity reconstruction.
There are three steps in AET to retrieve the conductivity in the DOI as well-known: measure the potential on the DOI boundary; compute the power density from this potential; use it to reconstruct the distribution of conductivity in the DOI. This is sketched in Fig. 1. In the second step, the reconstruction of the power density can be achieved from potentials on the whole DOI boundary, following [1]. Yet a limited number of electrodes is set in practice on the DOI boundary [9], and only their voltages are measured. That is, the second step cannot be achieved directly, since one needs to retrieve the potential on the whole DOI boundary before computing the power density.
In the existing literature, almost all applications of algorithms for AET consider that the power density is known, so focus is mostly on the third step of AET as investigated earlier in [10]- [12]. In the present contribution, all three steps are considered, CEM being used to simulate the measured voltages on electrodes since it appears as the model that provides simulated data best matching measured ones [9].
In a previous work [10], it was shown that the boundary potential converges much faster compared to the convergence rate of the conductivity in an iterative reconstruction process for EIT or AET -mainly because the boundary voltages measured on the DOI boundary are not sensitive to the interior changes of conductivity, e.g., [8]. This enabled to accelerate the convergence of the optimization procedure for the AET based on CEM and the Levenberg-Marquardt algorithm, a good performance being achieved.
In the present contribution, this observation has led to a method to reconstruct the boundary potential from an EIT measurement (voltages on a limited number of electrodes), which provides the necessary data to compute the power density for the retrieval of the conductivity in the DOI. In doing so, the aforementioned second step is completed, and hence a complete framework for AET follows.
To illustrate the feasibility of the proposed framework for AET, the CEM is used to simulate the voltages on the electrodes. A method based on singular value decomposition (SVD) is proposed and applied to obtain the boundary potentials from the boundary voltages measured, a lowaccuracy conductivity distribution being also obtained during this first step. The latter is taken as initial guess in the next step, a Levenberg-Marquardt method as discussed in [10], [13] yielding the conductivity distribution with high accuracy from the power density.
The contribution is organized as follows. In section II, the complete electrode model is introduced and a singular value decomposition method for EIT is investigated. The Levenberg-Marquardt method based on a continuum forward model is also presented. In Section III, numerical experiments are carried out. The proposed strategy is illustrated and validated with the methods introduced in the previous section. Conclusions and perspectives are in section IV.

II. RECONSTRUCTION ALGORITHMS
A DOI imaged by AET is modeled as a bounded Lipschitz domain ⊂ R n , n ≥ 2. The changes caused by the acoustic wave can be recorded with EIT measurements on the boundary ∂ [14], [15]. The power density in is defined as Here, σ is the conductivity distribution in domain , and u(σ ) is the electrical potential produced by applying either voltage or electric current on ∂ . Given noisy measurements E δ (σ ) of the true power densities E(σ ), the problem is to find the conductivity map σ by minimizing the functional E(σ ) can be retrieved from the boundary measurements [1], [16]. Yet the non-linear relationship between E(σ ) and σ renders this problem nonlinear [17]. The proposed strategy for AET is sketched in Strategy 1. Data U l, and U l are the measured voltages on the electrodes e l , l = 1, 2, . . . , N , with and without the perturbation of the focused ultrasonic wave, respectively. In this contribution, those are simulated with the CEM model, as the most practical model for EIT, and able to simulate the electrical potential with very good accuracy [18]. It is given with equations (3): N electrodes are attached at the boundary ∂ . A known total current I l is injected through the l-th electrode. Here, ν is the outward unit normal vector to boundary ∂ , e l the area occupied by the l-th electrode, and ∂/∂ν indicates the directional derivative of u along ν, z l being the contact impedance corresponding to the l-th electrode. To ensure existence and uniqueness of the solution, this model must include the law of VOLUME 8, 2020 charge conservation: N l=1 I l = 0, and N l=1 U l = 0 is also required to uniquely determine the solution of (3). Replacing the boundary conditions (3b)-(3d) in equations (3) yields the continuum model.
It is assumed that the conductivity map σ ∈ H s ( ) for s > n/2, where H s ( ) is the Sobolev space. It implies that the Sobolev space H s ( ) is a Banach algebra w.r.t point-wise multiplication [19]. The ''discontinuous'' Robin-type boundary conditions (3d) are such that the regularity of the CEM solution u(σ ) is limited in H 2− ( ) [9] for > 0, so E(σ ) : H s ( ) → H 1− ( ) with singularities at the edges of the electrodes -those can be eliminated with a contact impedance which exhibits smoothed value changes set at edges of the electrodes.
As indicated already, this step is performed with SVD. Potentials u b and u b, on the entire boundary of the DOI are retrieved from U l and U l, , respectively, then the power density proceeds from u b and u b, by using the method introduced in [1].

Strategy 1:
Proposed strategy for the retrieval of the conductivity distribution in the DOI from boundary measurements with a limited number of electrodes.
Data: Measured U l and U l, on electrodes e l .

Result:
The conductivity distribution σ r in DOI. 1 determine u b and σ 0 from U l based on EIT; 2 determine u b, and σ 0, from U ,l based on EIT; 3 determine E(σ ) from u b and u b, ; 4 determine σ r from E(σ ) with σ 0 as initial guess.

A. RETRIEVAL OF u b AND E(σ )
With the electrodes configuration as in (3), the space is introduced, letting l = 1, . . . , N . Here χ l is the indicator function of the l-th electrode and L 2 (∂ ) = {f ∈ L 2 (∂ ) : ∂ f = 0}. The weak formulation of the CEM reads [18] where It is assumed U = H 1 ( ) herein. As CEM only provides N (N − 1)/2 independent measurements, conductivities can only be recovered to a degree-of-freedom less than the number of independent measurements. The forward CEM can be represented with current-to-voltage map ) being a small perturbation, and σ 0 ∈ S + as background conductivity in the DOI. Linearizing around σ 0 yields which enables to predict the influence of h. Here, is the Fréchet derivative of T (σ, I ) along the direction h. Due to the symmetry of T (σ 0 , I ), easily validated by (5), the map T (σ, I ) is represented with a symmetric matrix in the order of dimension N − 1, written as where I i and I j are the ith and jth column vectors of an is thus represented in matrix form [20] as u(I i ) as unique solution of (5) w.r.t. the electrode current I i .
In numerical settings, can be described via a triangulation T , so as (10) is approximated by (12) With (9) and (12), (8) can be written into a linear system with h = {h 1 , h2, . . . , h T }, and b a vector obtained via reshaping the matrix The conductivity distribution can be obtained with an iterative algorithm [21] or singular value decomposition, as herein. In vector form, the reconstructed perturbation h r is Here, s i are singular values of T , u i and v i the corresponding singular vectors, N the number of singular values used. With h r , the boundary potential u b can be obtained by solving (5) with σ r = σ 0 + h r . With the perturbation introduced by the focused ultrasonic wave, the boundary potential u b, can be retrieved as the one for u b . According to [1], the power density can in principle be obtained from u b and u b, with where ν(x) is the ratio between perturbed and original volumes, and where φ I is the current distribution on the boundary of the DOI, both of them being known functions.
An iterative method based on the Levenberg-Marquardt algorithm (LMA) [22]- [24], as a well-known solution method for non-linear least-square problems, can yield the distribution of σ from E(σ ). 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 letting α k be the Tikhonov regularization parameter, and δσ = σ − σ k . This minimization problem regarding to (16) is equivalent to the linearized one letting F (σ k ) * be the adjoint of F (σ k ). In [13] and [24], It has been proven that, with properly chosen parameter α k and an initial guess σ 0 sufficiently close to the desired solution, the LMA converges to a solution σ δ of F(σ δ ) = E(σ ) δ . The LMA based on (17) can be seen 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 it 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 it behaves like a Gauss-Newton method, which is faster than the steepest descent. This method is used here for our fourth step in Strategy 1. The Fréchet derivative and adjoint operator needed in (17) are computed with a method proposed in [13], as follows. It is assumed that σ is positive, and σ ∈ H l ( ) for l > n/2. For f ∈ H l+ 1 2 ( ), then (3) has a unique solution u(σ ) ∈ H l+1 ( ), which indicates that the operators E(σ ), E (σ ) and E (σ ) * are all maps from H s ( ) to H s ( ) [13].
Here, it is assumed that τ vanishes in a domain = {x ∈ |d(x, ∂ ) < } for a sufficiently small . So, all terms associated to τ in the boundary conditions disappear. It is well-known that the solution of equations (20) uniquely exists in the weak sense. u (σ ) as defined in (20) is the Fréchet derivative of u(σ ).
β is a regularization parameter. Upon integration by parts, B * ζ is the solution of the Neumann problem with ς ∈ L 2 ( ). Refer to [13] 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 (17), 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 (30), Here, M τ is easily obtained with (21), (23). After computing τ k from (30), the conductivity map σ k obtained at the k-th iteration is updated by σ k+1 = σ k + τ k for a new iteration.
The iterative construction method based on the above system is named as LM-CM, refer to Algorithm 2 for a single measurement. Since ν(x) is not known for the current numerical study, the data of E δ are simulated by solving (3a) with 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. Since η cannot be determined in practice, a relatively large value can be given to α 0 , and α k slowly decreased to ensure convergence. The iteration is stopped if τ k L 2 is smaller than the expected error or if the maximum number of iterations is reached.

Algorithm 2: LM-SCEM Algorithm for Retrieving the Conductivity Map of a Domain From a Single Measurement of Power Density
Data: Measured power density E δ and initial guess σ 0 Result: Retrieved conductivity map σ r with relative error η 1 N : maximum number of iterations; 2 σ k ← σ 0 ; 3 α k ← α 0 ; 4 num ← 1; 5 norm ← 1; 6 while norm > δ and num < N do 7 Update u k from σ k with CEM; (24)  The numerical experiments are carried out with a lung-heart model depicted in Fig. 2a. Three tissues are considered: heart (red, σ = 0.7 S/m), lung (cyan, σ = 0.26 S/m), and soft-tissues (blue, σ = 0.33 S/m). The model is in a circular region with background material (white, σ = 0.20 S/m) and radius r = 23 cm. So, the conductivity distribution of the phantom is piecewise. The boundaries between two regions with different conductivities are smoothed via a convolution between conductivity distribution and mollification equation which reduces interpolation errors at the boundaries. The mollified conductivity distribution of the model is displayed in Fig. 2b,, and it is referred to as σ 0 . All electrodes are uniformly distributed on the boundary of the DOI, indicated as the red squares in Fig. 2a.  A perturbation h, refer to Fig. 3(a), is introduced in the model. The perturbation value in the red circular domain is given as 0.2 S/m, and the values in other parts are all zero. The conductivity distribution of the perturbed model is shown in Fig. 3(b), which is simple σ 0 + h. The relative error of the reconstruction is defined ccording to the L 2 norm, δ = where M is the (maximum) current amplitude. The vector I n is a current pattern, which is a discrete approximation to M cos(nθ ) for 1 < n < L/2 or M sin(nθ ) for L/2 − 1 ≤ n ≤ L − 1. The vectors I n are orthonormalized to form a complete set of bases for the space V . After computing the matrices T σ 0 +h , T σ 0 and T from the data obtained with the complete electrode model, the vector h can be computed from (14).  The singular values of the matrix T are ordered from large to small and shown in Fig. 4. Generally, the singular vectors associated to the large ones give more information about the region close to the DOI boundary, and those assocoated to the smaller ones give more information about the interior region. However, if the singular value s i is too small, a very large influence will be put on h r since s i is at the denominator of (14). Thus, only a limited number of singular values is used to reconstruct h r .
The variation of the relative error δ with the number of singular values employed in the reconstruction is shown in Fig.6. The best results are obtained by using the first 175 largest singular values with relative error δ = 0.6141. The retrieved h in the DOI is shown in Fig. 5a. Obviously, position and profile of the perturbation can be suitably retrieved with the SVD-based method, but the relative error is not good enough for a good reconstruction of h r . This is mainly caused by the Robin-type boundary condition applied in the CEM. This leads to the computed potential existing in H 2− when σ ∈ H 2 , thus the computed power density E(σ k ) is in H 1− for any > 0. In particular they are not in H with > n/2 even in one dimension [27], [28]. Yet, it will be shown below that the value obained is good enough for determination of the boundary potential u b .
To verify that the measured boundary potential in EIT is not very sensitive to the changes of the interior conductivity, the boundary potentials for the perturbed conductivity map σ p,t and the retrieved one σ p,r are computed using  the complete electrode model with different current patterns applied on the electrodes. u b,t and u b,r are the boundary potentials computed with σ p,t and σ p,r , respectively. The relative error is defined as The variation of u b,t with central angle of the observation point on the boundary is shown in Fig. 7, and the difference |u b,t − u b,r | is in Fig. 8. It is seen that the difference between boundary potentials with σ p,t and σ p,r is very small. Specifically, the relative error δ b is shown in Fig.9. It is close to zero save in the range where u b,t is approximately zero. Indeed, |u b,t | is the denominator in the definition of δ b ; when |u b,t | is about zero, a very small error in |u b,t − u b,r | leads to a large value in δ b .
According to (15), i.e., the relation between E(σ ) and u b, − u b , the noise level in E(σ ) is decided by the noise in u b, − u b . So, as long as the error in u b, − u b is reasonable,  the power density E(σ ) can be computed to retrieve the conductivity distribution with methods well developed in the literature.
The Levenberg-Marquardt algorithm (LMA) can be used to get the conductivity distribution σ from the obtained power density E(σ ), which is the 4th step in the strategy 1. To compute the power density from equation (15) calls for the knowledge of ν(x), which cannot be obtained here. Therefore, E(σ ) is computed directly from its definition σ |∇u| 2 , where u is obtained with the complete electrode model introduced in section II. The DOI is represented by 15690 triangles for the forward computation and 14301 triangles for the inversion. The complete electrode model and current patterns I 1 and I 16 are used to compute potentials u 1 and u 16 from the σ r . One stops when δ < 10 −3 or the number of iterations is larger than 30. Parameter α k for the k-th iteration is α k = α 0 /a k with α 0 = 50.0 and a = 1.5. The value of β is set as 10 −2 , but it could be greater or smaller whether needed. A large value of β can provide a better noise tolerance of the Levenberg-Marquardt method. The result shown in Fig. 5b obtained via the SVD-based method is taken as initial guess in the iteration. Noise is appraised with signal-to-noise ratio (SNR) SNR = 20 log 10 where N is a Gaussian white noise distribution. Numerical examples involve different levels of noise in the data E(σ ). Fig. 11a gives the results with SNR = ∞, i.e., no noise in the data E(σ ). Results are much better that those with the SVD-based method. But, in practice, noise is faced with. So, the noise tolerance of the Levenberg-Marquardt approach is a key factor. It is investigated by adding to the data E(σ ) different levels of noise. SNR = 40, 20, and 10 dB are considered. The variation of the relative error δ is given in Fig. 10. With 40 dB, LMA provides a good convergence, and the relative error decreases as fast as without noise, and quickly achieved at δ = 0.028 at iteration step 12, see Fig. 11b. Afterwards, a divergence is observed, and its reason is given in the previous section.
But when the noise is increased to 10 dB, the performance of LMA worsens. However, it is observed from Fig. 11c for 20 dB that the result is still much better than the one obtained with the SVD-based method.
When SNR = 10 dB, in effect the method does not produce convergent results anymore, see Fig. 10, noticing that more than 30% of noise is added to the data E(σ ). In this case, a better noise tolerance is required. The value of β is increased from 1 × 10 −2 to 0.1, and a good convergence is observed from the dash-dotted curve in Fig. 10, the best result obtained at step 12 being given in Fig. 11d, which is much better than with the SVD-based method.
In brief, the large value of β ensures proper noise tolerance, yet this value cannot be too large to make the algorithm converge well or to avoid over-smoothed results.
These numerical simulations show that the Levenberg-Marquardt approach provides the conductivity distribution of a closed domain with high accuracy, and has good noise tolerance. With high noise level, the conductivity can still be well mapped if proper value of β.

IV. CONCLUSION
An hybrid strategy has been proposed for acousto-electric tomography in order to retrieve the conductivity distribution of the DOI from voltages measured on electrodes attached at its boundary. How to implement this strategy has been detailed, and illustrated by comprehensive numerical simulations.
An SVD-based method is used as boundary potential reconstruction method to compute the electrical potential on the whole boundary, and it is shown that, with a low quality reconstruction of the interior perturbation of the conductivity in the DOI, the boundary potential can be computed with good accuracy. Meanwhile, good noise tolerance for achieving high-quality conductivity distribution from the power density is critical for the final result.
The iterative Levenberg-Marquardt method is used to retrieve the conductivity from the computed power density. Strong noise tolerance is observed and the quality of the final conductivity map much improves, e.g, a good retrieval with relative error less than 3% is achieved in a dozen of iterations for low-level noise, and good results are still obtained with higher-level noise.
The work summarized herein can be seen as within a continuum of investigations of AET as illustrated by the companion contribution [10]. Next, one should on the one hand focus onto higher accuracy potential reconstruction methods, and on the other hand onto comprehensive validations from experimental data in laboratory-controlled environments.