Regularization of systems of nonlinear ill-posed equations: II. Applications

In part I we introduced modified Landweber-Kaczmarz methods and have established a convergence analysis. In the present work we investigate three applications: an inverse problem related to thermoacoustic tomography, a nonlinear inverse problem for semiconductor equations, and a nonlinear problem in Schlieren tomography. Each application is considered in the framework established in the previous part. The novel algorithms show robustness, stability, computational efficiency and high accuracy.


Introduction
In [8] we analyzed novel iterative Landweber-Kaczmarz methods for solving systems of ill-posed equations F i (x) = y δ,i , i = 0, . . . , N − 1 . (1) Here F i : D i ⊆ X → Y are operators between Hilbert spaces X and Y and y δ,i are approximations of noise free data y i with y δ,i − y i ≤ δ i for i ∈ {0, . . . , N − 1}.
In [8] we have proven stability and convergence for both the eLK and lLK method. In this article we apply the methods to three different problems: a linear inverse problem related to thermoacoustic tomography, an inverse problem for semiconductors and Schlieren tomography. For validation, the results are compared to the classical Landweber-Kaczmarz (LK) method [18,11] x n+1 = x n − F [n] (x n ) * (F [n] (x n ) − y δ, [n] ) .
The outline of this article is as follows. In Section 2 we apply both methods, the lLK method and the eLK method, to an inverse problem related to thermoacoustic computed tomography, which mathematically can be reduced to the inversion of the circular mean transform. In Section 3 we consider the exponentially ill-posed nonlinear inverse problem of estimating the doping profile in a semiconductor from voltage-current measurements. In Section 4 we consider Schlieren tomography, an imaging method for reconstructing three dimensional pressure fields.

Thermoacoustic computed tomography
Thermoacoustic computed tomography is a promising new imaging modality that has the potential to become a mayor non-invasive imaging method. In thermoacoustic imaging an object of interest is illuminated by short electromagnetic pulses, such as optical illumination or radio waves, which results in an excitation of acoustic waves (pressure variations). The spatially varying initial pressure distribution inside the object caries valuable structural and functional information and is reconstructed from acoustical data which is recorded with detectors outside of the object. Among the several publications on thermoacoustic imaging we mention [2,12,22,23,24].

Mathematical model
Assume that the initial pressure distribution x(ξ 1 , ξ 2 ) is independent of the third spatial coordinate, let ξ := (ξ 1 , ξ 2 ) and let B R denote the open disc (in R 2 ) with radius R.
with noisy data y δ,i , where is the scaled mean value of x ∈ C ∞ 0 (B R ) over the circle with center ξ i ∈ ∂B R and radius t ∈ [0, 2R]. Here S 1 denotes the surface of the unit disc in R 2 and dΩ the line measure on S 1 and (ξ, 0) corresponds to the locations of an acoustic receiver, see Figure 1. This section is devoted to the stable solution of (7).
Recently, we proposed a detection technique [5] where an array of parallel line detectors is rotated around the around the object and pressure signals along the line detectors are recorded. Then the solution of (7) allows for reconstructing of a fully three dimensional initial pressure distribution. In experiments line detectors are realized by thin laser beams that are either part of a Fabry-Perot or a Mach-Zehnder interferometer [5,19].

Abstract formulation in Hilbert-spaces
In the following we show that the problem of recovering of a function x from its circular means can be put in the abstract framework of [8].
In the remainder of this section, let L 2 (B R ) be the Hilbert space of square integrable functions on B R with x 2 := B R x(ξ) 2 dξ. Moreover, we denote by L 2 ([0, 2R], tdt) the Hilbert space of all functions y : [0, ∞) → R (observable quantities) with support in [0, 2R] with Finally we denote by the associated inner product on L 2 ([0, 2R], tdt).
The scaled circular mean operators are defined by where M i (x) is defined by (8).
Proof. Assume that x ∈ C ∞ 0 (B R ). From the definition of M i (x), the Cauchy Schwarz inequality and Fubini's theorem it follows that Hence M i is bounded from C ∞ 0 (Ω) into L 2 ([0, 2R], tdt) and therefore can be extended to a bounded linear operator mapping on L 2 (B R ) with M i ≤ 1. In particular M i has a bounded adjoint.
Let x ∈ C ∞ 0 (B R ) and y ∈ C ∞ 0 ([0, 2R]). From Fubini's theorem and the substitution ξ = ξ i + tσ it follows that This shows (9). We note that, for linear bounded operators, the tangential cone condition [8, Eq. (15)] is satisfied with η = 0. Since M i is linear and M i ≤ 1, the lLK method and the eLK method provide convergent regularization methods for solving (7).

Numerical reconstruction
The lLK method applied to thermoacoustic computed tomography reads as with In the numerical implementation the spaces L 2 (B R ) and L 2 ([0, 2R], tdt) are approximated by the linear spans of piecewise linear splines. Each spline is represented by vectors in R M ×M and R M , respectively. For the numerical evaluation of M i the integration over S 1 in (8)  . This is the same complexity that is needed to perform one step in the Landweber iteration.
In the following numerical examples we consider N = 80 measurements, where the centers ξ i = R(sin(πi/N ) cos(πi/N )) are uniformly distributed on the semicircle S + := {(ξ 1 , ξ 2 ) ∈ ∂B R : ξ 1 ≥ 0}. The phantom x, shown in the right picture in Figure 2, consists of a superposition of characteristic functions and one Gaussian kernel. The data M i (x) was calculated via numerical integration and 5% uniformly distributed noise was added. Micro-local analysis predicts, that in such limited angle situation certain details of x outside the convex hull of S + cannot be recovered in a stable way [16,25]. The numerical reconstruction  (10), (11). The right picture shows the reconstruction with the LK method using the same number of n * δ of iterations.  with the lLK method, implemented with τ = 2.0, is depicted in Figure 3. The stopping rule becomes active after 8 cycles. For comparison purposes, the result of the LK method is plotted after 8 cycles. Figure 5 shows results obtained with the eLK method using the same phantom. We remind that the iterates of the eLK method are vectors For the results depicted in Figure  5 we used the average x n := n over all components.

Discussion
The left image in Figure 4 shows the decrease of error for both the LK method and the lLK method. The number of actually computed Landweber steps in the lLK method (number of iterations with ω n = 1 within one cycle) of the lLK method is shown in the right image in Figure 4. The lLK method provides a better approximation of the exact solution than the LK method. Moreover, since more than half of the iterates are loped (see right picture in Figure 4), the numerical effort is remarkably smaller. Since in the eLK method an averaging over all components of the solution vector x n δ is performed, the quality of the reconstruction is slightly better than for the lLK method (and other tested regularization methods). Shortcoming of the eLK method are the larger amount of memory and the larger number of iteration cycles needed.

An inverse problem for semiconductors
In this section we investigate the solution of an inverse problem for nondestructive testing of semiconductors. More precisely, we aim to recover the doping profile on the basis of a simplified drift diffusion model from measurements given by the voltage-current (VC) map [15,14]. The precise implantation of the doping profile is crucial for the desired performance of semiconductor devices in practice. In order to minimize manufacturing costs of semiconductors as well as for quality control, there is substantial interest in replacing expensive laboratory testing by numerical simulation for non-destructive evaluation.

Mathematical modelling
Let Ω ⊂ R d , d ∈ {1, 2, 3}, be a domain representing the semiconductor device and let ν denote the unit normal vector to its boundary ∂Ω. The boundary of Ω is divided into two nonempty disjoint parts ∂Ω N and ∂Ω D . The Dirichlet part of the boundary ∂Ω D models the Ohmic contacts, where an electrostatic potential is induced. The Neumann part of the boundary ∂Ω N corresponds to insulating surfaces.
The doping profile C : Ω → R (unknown parameter) models the preconcentration of ions in the crystal, which is produced by diffusion of different materials into the silicon crystal and by implantation with an ion beam. In particular, C = C + − C − , where C + and C − are concentrations of positive and negative ions, respectively. In those subregions of Ω in which the preconcentration of negative ions predominate (P-region), we have C < 0. Analogously, we define the N-region, where C > 0. The boundaries between the P-regions and N-regions (where C changes sign) are called the PN-junctions, see Figure 6.
We review the linearized stationary unipolar mathematical model, which is described by a drift diffusion equation (see [15,14]): in Ω (12) and Here, the function V : Ω → R denotes the electrostatic potential (−∇V is the electric field) and u : Ω → R corresponds, up to an exponential transformation, to the concentration of free carriers of negative charge (the concentration of free positive charge is assumed to vanish; the reason why the model is called unipolar). The positive constant µ n is related to the mobility of electrons and λ 2 is the Debye length.
At the Dirichlet part of the boundary ∂Ω D the potential V and the concentration u are prescribed: V bi is a given logarithmic function [14], and the function U denotes the applied potential. At the the Neumann part of the boundary ∂Ω N zero current flow (17) and a zero electric field in the normal direction (14) are prescribed.
In this section we are concerned with determining the inverse doping problem for the VC map, which consists in the determination of the doping profile C in system (12)(13)(14)(15)(16)(17) from measurements of the VC map which maps an applied potential U at ∂Ω D to the current flow Σ C (U ) through the contact Γ 1 . Here (V, u) solve (12)(13)(14)(15)(16)(17) and ∂u/∂ν is the normal derivative of u. For results on existence and uniqueness of H 1 -solutions of system (12)(13)(14)(15)(16)(17), we refer to [15]. Since systems (12)(13)(14) and (15)(16)(17) are decoupled, we can split the inverse problem in two parts: First, the inverse problem of identifying the parameter x(ξ) := e V (ξ) in (15)(16)(17) from measurements of the Dirichlet to Neumann (DN) map is solved. The second step consists in the determination of the doping profile from (12), namely C = x − λ 2 ∆(ln x). Notice that the second step is related to twice numerical differentiation which is mildly ill-posed (see [7]). On the other hand, the first step corresponds to the inverse problem of impedance tomography (EIT) with partial data, which is known to be severely ill-posed. (A review of EIT can be found in [3].) For the sake of clarity of presentation we will focus on the problem of identifying x in the first step above.

Abstract formulation in Hilbert space
Due to the nature of the practical experiments, some restrictions on the data have to be taken into account: 1. The voltage profiles U ∈ H 3/2 (∂Ω D ) must satisfy U | Γ1 = 0.
2. x has to be determined from a finite number of measurements, i.e. from the given data where U i ∈ H 3/2 (∂Ω D ) are prescribed voltage profiles satisfying U i | Γ1 = 0.
In order to apply the abstract framework of [8] we write the inverse doping problem as a system of operator equations F i (x) = y δ,i , i = 0, . . . , N − 1 . Here x ∈ L 2 (Ω) =: X is the unknown parameter, y δ,i ∈ R denotes measurement data and F i : X → Y are the parameter to output maps, with domain of definition Although the operators F i are Fréchet differentiable, they do not satisfy the tangential cone condition [8,Eq. (13)]. Therefore, the convergence results derived in [8, Section 2] cannot be applied.

Numerical reconstruction
In the following numerical examples we assume that N = 9 Dirichlet-Neumann pairs (U i , F i (x)) of measurement data are available. The domain Ω ⊂ R 2 is the unit square, and the boundary parts are defined as follows The fixed inputs U j , are chosen to be piecewise constant functions supported in Γ 0 where the points s j are uniformly distributed on Γ 0 and h = 1/32. The doping profile to be reconstructed is shown in Figure 7 (a). In Figure 7 (b) a typical voltage profile U j (applied at Γ 0 ) is shown as well as the corresponding solution u of (15)(16)(17). In these pictures, as well as in the forthcoming ones, Γ 1 appears on the lower left edge and Γ 0 on the top right edge (the origin corresponds to the upper right corner). In Figure 8 we show the evolution of the iteration error for both the LK method (a1-a3) and the lLK method (b1, b2). The number of actually computed Landweber steps within each cycle of the lLK method is shown in Figure 8 (c).  The stopping rule for the lLK method with τ = 2.0 is satisfied after 52 cycles. In Figures 8 (b1, b2) one can see the iteration error after 10 and 52 cycles. For comparison purposes, the iteration error for the LK method is plotted in Figures 8 (a1-a3) after 10, 50 and 100 cycles.

Discussion
The LK method requires almost twice as much cycles as the lLK method in order to obtain a similar accuracy. The efficiency of the lLK method becomes even more evident when we compare the total number of actually performed iterations. Each cycle of the LK method requires the computation of 9 iterations, while in the lLK method the number of actually performed iterations is very small after a few iteration cycles. As one can see in Figure 8 (c), no more than 3 Landweber steps are computed after the 20-th cycle of the lLK method. In total for the computation of the approximation in Figure 8

Reconstruction of transducer pressure fields from Schlieren images
This section is concerned with the problem of reconstructing three dimensional pressure fields generated by a medical ultrasound transducer from Schlieren data.
The data are collected with a Schlieren optical system where in the experiment an acoustic pressure is emitted into a water tank (see Figure 9). The Schlieren optical system outputs the intensity pattern of light passing through the tank which is proportional to the square of the line integral of the pressure along the light path.
For practical aspects on the realization of Schlieren systems and more background information on such measurement devices we refer to [4,6,10,13,20,21,26].   Figure 10. In the following the height parameter z is fixed. We aim to reconstruct the pressure function x : D → R defined by x(ξ) := p(ξ, z) from data F i (x)(s) := P i (s, z). Now Schlieren tomography can be formulated as the problem of solving the system of operator equations F i (x) = y δ,i , i = 0, . . . , N − 1 , for given noisy data y δ,i .

Abstract formulation in Hilbert-space
denotes the Radon transform, and F i = R 2 i denotes the Schlieren transform. Moreover, F i is Fréchet differentiable and the Fréchet derivative of F i at x is given by Proof. Let x ∈ C ∞ 0 (D). By applying the Cauchy Schwarz inequality two times, Using Sobolev's embedding theorem (see [1]), it follows that for some positive constant C. Therefore, F i extends, by continuity, to an operator H 1 0 (D) → L 2 (D). Now let F i [x] be defined by (21). Obviously F i [x] is a bounded linear operator. Moreover, from the Cauchy Schwarz inequality it follows that Therefore, we have We denote by R i the adjoint of R i , considered as an operator from L 2 (D) into L 2 (I) [17], which is given by With this operator we can define the adjoint of F i [x]: in H 1 0 (D). Here ∆ denotes the Laplace operator on H 1 0 (D).
Proof. Let x ∈ H 1 0 (D). We remark that F i [x] * is defined by for all y ∈ L 2 (D) and x ∈ H 1 0 (D). (24), the left hand side of (26) can be rewritten as follows

Using Fubini's theorem and
2. The left hand side of (26) reads as follows since h vanishes at the boundary of D.

Numerical reconstruction
According to (25) the lLK method for Schlieren tomography reads as follows In the numerical implementation we approximated functions defined on D and I by piecewise linear splines. In the numerical experiment we used τ = 2.2 and N = 250. The synthetic data set y δ,i , i ∈ 0, . . . , N − 1, were generated by adding normal distributed random noise with 0.01% noise level to the exact data. The phantom to be reconstructed is shown in Figure 11. It is a superposition of several characteristic functions. The reconstructions with the lLK method and the LK method were performed using the constant function x 0 (ξ) = 0.01 as initial guess.
We note that in all numerical experiments the lLK method applied to Schlieren tomography converged, even if the tangential cone condition was not satisfied.

Discussion
In contrast to the FPB algorithm the lLK method and the LK method are able to reconstruct both the positive and negative part of u. The lLK method reduces the artifact which is present in the reconstruction via the LK method.  The left image in Figure 12 shows that the number of actually computed Landweber steps per cycle in the lLK method is rapidly decreasing. Moreover, the norm of the error in the reconstruction with the lLK method is below the error of the LK method. Actually, the regularized solution x n δ * of the lLK method is a better approximation of the true solution than all iterates of the LK method.
Since the used measuring data are squared numbers, there are two solutions, one with a positive sign and one with a negative sign. Our numerical simulations showed that a strictly positive initial guess x 0 leads to a numerical reconstruction x n δ * with positive mean value.

Conclusion
We applied the abstract theory of the first part of this article to thermoacoustic computed tomography, to an inverse problem for semiconductors and to Schlieren tomography. In all applications the lLK method turned out to be an efficient iterative regularization method.