Next Article in Journal
A Learning—Based Particle Swarm Optimizer for Solving Mathematical Combinatorial Problems
Next Article in Special Issue
2-D Elastodynamic Time-Reversal Analysis for Surface Defects on Thin Plate Using Topological Sensitivity
Previous Article in Journal
Randomly Stopped Sums with Generalized Subexponential Distribution
Previous Article in Special Issue
Reconstructing Loads in Nanoplates from Dynamic Data
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Using Alternating Minimization and Convexified Carleman Weighted Objective Functional for a Time-Domain Inverse Scattering Problem

by
Nguyen Trung Thành
Department of Mathematics, Rowan University, 201 Mullica Hill Rd., Glassboro, NJ 08028, USA
Axioms 2023, 12(7), 642; https://doi.org/10.3390/axioms12070642
Submission received: 30 April 2023 / Revised: 20 June 2023 / Accepted: 26 June 2023 / Published: 28 June 2023

Abstract

:
This paper considers a 1D time-domain inverse scattering problem for the Helmholtz equation in which penetrable scatterers are to be determined from boundary measurements of the scattering data. It is formulated as a coefficient identification problem for a wave equation. Using the Laplace transform, the inverse problem is converted into an overdetermined nonlinear system of partial differential equations. To solve this system, a Carleman weighted objective functional, which is proved to be strictly convex in an arbitrary set in a Hilbert space, is constructed. An alternating minimization algorithm is used to minimize the Carleman weighted objective functional. Numerical results are presented to illustrate the performance of the proposed algorithm.

1. Introduction

Scattering theory of waves is concerned with the effect scatterers have on incident waves [1]. There are two kinds of problems in the scattering theory. Forward scattering problems aim to determine the scattered waves, which are generated when incident waves interact with known scatterers. Inverse scattering problems (ISPs) aim to detect and identify scatterers from knowledge of the scattered wave(s) generated by one or multiple incident waves. The identification of scatterers is based on their geometrical or/and physical properties. Scatterers can be categorized as penetrable or impenetrable. Penetrable scatterers allow waves to penetrate and are characterized by both physical properties (such as the refractive index or permittivity) and geometrical properties (such as location and shape). By contrast, impenetrable scatterers (which are also called obstacles) do not allow waves to penetrate and are usually characterized by geometrical properties. ISPs are also divided into two types: inverse medium scattering problems and inverse obstacle scattering problems. The former aim to determine physical properties (which in turn provide information about geometrical properties) of penetrable scatterers, whereas the latter aim to determine the geometrical properties of impenetrable scatterers.
ISPs play a central role in many technologies such as ground-penetrating radar (GPR), microwave imaging, ultrasound imaging, and seismic imaging; see, e.g., [2,3,4,5,6] and the references therein. These technologies are used in various applications. GPR is an effective tool for detecting buried objects such as landmines, underground structures, or archaeological sites [5]. Ultrasound is widely used as a medical imaging tool. Recent research has shown the potential of microwaves in detecting breast tumors or brain strokes [7,8,9,10,11]. Seismic imaging is used for geophysical explorations [12].
In this paper, we are interested in inverse medium scattering problems for the Helmholtz equation in the time domain, which can be stated as follows. Let Ω be a convex bounded domain in R d , d 1 , with a piecewise-smooth boundary Ω . We denote by S T the lateral boundary of Ω × ( 0 , T ) , i.e., S T : = Ω × 0 , T . Consider a coefficient function c x which satisfies the following conditions:
c C 1 + α R d , c ̲ c x c ¯ , x R d , c x = 1 , x R d Ω ,
for some given constants c ̲ and c ¯ such that 0 < c ̲ 1 c ¯ . Here, C m + α denotes Hölder spaces, where m 0 is an integer and α 0 , 1 . Consider the following Cauchy problem:
c x u t t = Δ u x , t + f ( x , t ) , x , t R d × 0 , ,
u x , 0 = 0 , u t x , 0 = 0 .
The coefficient c x represents a physical property of the medium (including the scatterers) in which the waves propagate. For electromagnetic waves, it represents the dielectric constant, whereas for acoustic waves, it represents the refractive index. Here, c ( x ) is normalized so that its value in the background medium, i.e., in R d Ω , equals 1. The state function u x , t represents the wave (or one component of the wave) propagating in the medium. The function f x , t 0 represents the source of the wave. In this work, we assume that the source is located outside of the domain Ω , i.e.,
f ( x , t ) = 0 , x Ω .
The inverse medium scattering problem is stated as a coefficient identification problem as follows:
CIP: Suppose that conditions (1)–(4) are satisfied. Determine the coefficient c ( x ) for x Ω , from the boundary data:
g 0 ( x , t ) : = u ( x , t ) , g 1 x , t : = u ( x , t ) ν , ( x , t ) S ,
where ν = ν ( x ) is the outward unit normal vector of Ω.
Remark 1.
1. In practice, the Neumann data, i.e., g 1 , are usually not measured. However, if the Dirichlet data, i.e., g 0 , are measured on the whole boundary Ω , we can compute g 1 by solving the initial value problem (2) and (3) in R d Ω .
2. In computation, the infinite time interval [ 0 , ) in (5) can be replaced by a finite interval [ 0 , T ] . The reason is that the incident wave is usually excited for a short period of time. Therefore, the wave function u is usually too small to be useful after some finite time. Moreover, in the numerical method we discuss in this paper, data at large time instants are essentially eliminated when a Laplace transform is used. The infinite interval is used for the convenience of theoretical analysis only.
In this work, we assume that the CIP with noiseless data has a unique solution. Our investigation is devoted to numerical methods for solving it. The most common numerical methods for solving this inverse problem use the least-squares minimization approach. These methods find an approximation of the coefficient by minimizing a least-squares objective functional. Due to the nonlinearity of the CIP, the objective functional is nonconvex. Therefore, a good initial guess is usually needed for the convergence of gradient-based iterative optimization methods to the global minimizer of the objective functional.
We consider in this paper a numerical method for solving the above CIP without requiring a good initial guess. The idea of this method is to transform the CIP into a system of nonlinear partial differential equations with overdetermined boundary conditions. To solve the resulting nonlinear system, we minimize a Carleman weighted objective functional which can be proved to be strictly convex in an arbitrary domain if the Carleman weight chosen is sufficiently large. The approach of convexification using Carleman weighted objective functional was introduced by Klibanov in the 1990s and then expanded to several inverse problems by himself and his collaborators; see, e.g., [13,14,15,16,17,18,19]. This paper extends the method we developed in our previous work [14]. The novel contributions of the current paper include the following: (1) A new first-order PDE system is obtained instead of a second-order system as in [14]; (2) the convexity of the discrete Carleman weighted objective functional is proved; and (3) we propose an alternating minimization method to solve the Carleman weighted objective functional instead of the steepest gradient algorithm used in [14]. The advantages of the proposed method include the following:
  • Unlike the least-squares approach, the proposed method does not require a good first guess. Indeed, it was demonstrated through numerical examples in our previous work [14] that the least-squares method failed to converge at initial guesses at which the convexification method was able to provide a reasonably good approximation of the coefficient.
  • Compared to our previous work [14] and the above references in the convexification approach, each minimization sub-problem in our proposed method has a much smaller number of variables than the original nonlinear system. Consequently, the sub-problems can be solved more quickly. We observed in our numerical tests that the alternating minimization method converges much faster than the steepest descent method.
As another interesting approach for solving ISPs that has received a lot of attention in the last year, we refer to various deep learning-based methods; see, e.g., [20,21,22,23,24] for an incomplete list of recent works focusing on this direction. Some results for microwave imaging data for medical applications have been reported in [25,26,27]. Numerical results of this approach appear promising, although the mathematical foundation of deep learning is still to be investigated.
The rest of the paper is organized as follows. In Section 2 we describe the transformation from the CIP into a nonlinear system of coupled partial differential equations. In Section 3 we describe the Carleman weighted objective functional and the alternating minimization algorithm. Section 4 is devoted to the discretized problems and its convexity. Numerical examples are presented in Section 5. Finally, we make some concluding remarks in Section 6.

2. Transformation of the CIP into a Nonlinear System

We first transform the CIP into a nonlinear system of coupled partial differential equations using the method described in our previous work [14]. For that purpose, we use the Laplace transform:
u ˜ ( x , s ) : = L u x , s = 0 u x , t e s t d t , s > 0 ,
where s is referred to as the pseudofrequency. Let s ̲ be a positive constant depending only on c ¯ such that the Laplace transforms of u and its derivatives D β u , β = 1 , 2 , converge absolutely for all s s ̲ . It follows from (2) and (4) that u ˜ satisfies the equation
Δ u ˜ ( x , s ) s 2 c ( x ) u ˜ ( x , s ) = 0 , x Ω , s s ̲ .
For appropriate choices of the source function f, we can prove that u ˜ ( x , s ) > 0 for s s ̲ ; see Theorem 3.1 of [28]. In this work, we assume that the source function is chosen such that u ˜ ( x , s ) > 0 for s s ̲ . We define a new vector-valued function v ( x , s ) as
v ( x , s ) : = u ˜ ( x , s ) s 2 u ˜ ( x , s ) .
From (7) we obtain
Δ u ˜ ( x , s ) = · ( s 2 v ( x , s ) u ˜ ( x , s ) ) = s 2 ( · v ( x , s ) ) u ˜ ( x , s ) + s 4 | v ( x , s ) | 2 u ˜ ( x , s ) .
Substituting (8) into (6), we obtain the following equation:
· v ( x , s ) + s 2 | v ( x , s ) | 2 = c ( x ) , x Ω .
Equation (9) provides an explicit way to compute the coefficient c ( x ) if v ( x , s ) is known. Unfortunately, both functions c and v are unknown in this equation. However, since the right-hand side of (9) does not depend on s, we can eliminate the unknown function c ( x ) by differentiating (9) with respect to s. Denoting by q ( x , s ) : = v ( x , s ) / s , we obtain the following integral differential equation for q:
· q ( x , s ) 2 s 2 q ( x , s ) · s q ( x , τ ) d τ + 2 s | s q x , τ d τ | 2 = 0 , x Ω , s s ̲ .
In addition, q satisfies the following boundary condition:
q ( x , s ) · ν = φ ( x , s ) : = s g ˜ 1 ( x , s ) s 2 g ˜ 0 ( x , s ) , x Ω ,
with g ˜ 0 and g ˜ 1 being the Laplace transforms of g 0 and g 1 , respectively.
We note that Equation (10) is not easy to solve. Using the same approach as in our previous work [14], we represent the function q ( x , s ) as a series with respect to an orthonormal basis of L 2 ( s ̲ , ) . For this purpose, we consider Laguerre polynomials [29]:
L n s = e s / 2 k = 0 n 1 k C n k s k k ! , s 0 , , C n k = n ! n k ! k ! .
It can be verified that the set f n s n = 0 , where f n s : = L n s s ̲ , s s ̲ , , is an orthonormal basis in L 2 s ̲ , . We then approximate the function q x , s as follows.
q x , s n = 0 N 1 q n x f n s , s s ̲ ,
where N is a sufficiently large integer. The choice of N will be discussed in our numerical examples presented in Section 5. Substituting the truncated series (12) into (10), we approximate the coefficients q n ( x ) , n = 0 , , N 1 , by the solution of the following equation:
n = 0 N 1 · q n ( x ) f n ( s ) 2 s 2 m = 0 N 1 n = 0 N 1 q m ( x ) · q n ( x ) f m ( s ) s f n ( τ ) d τ + 2 s m = 0 N 1 n = 0 N 1 q m ( x ) · q n ( x ) s f m ( τ ) d τ s f n ( τ ) d τ = 0 .
Multiplying both sides of (13) by f k ( s ) , integrating over ( s ̲ , ) , and using the orthonormality of f n ( s ) n = 0 in L 2 ( s ̲ , ) , we obtain the following system of coupled nonlinear equations:
· q k ( x ) + m = 0 N 1 n = 0 N 1 F k m n q m ( x ) · q n ( x ) = 0 , k = 0 , , N 1 , x Ω ,
where the numbers F k m n , k , m , n { 0 , , N 1 } , are given by
F k m n = s ̲ 2 s f k ( s ) ( s f m ( τ ) d τ s f n ( τ ) d τ ) d s s ̲ 2 s 2 f k ( s ) f m ( s ) ( s f n ( τ ) d τ ) d s .
The boundary conditions for q n are obtained by substituting again the truncated series (12) into (11). More precisely, we have
q k ( x ) · ν = φ k x , x Ω ,
with
φ k x = s ̲ φ ( x , s ) f k ( s ) d s .

3. Carleman Weighted Objective Functional and Alternating Minimization Algorithm

We note that the unknown functions q k ,   k = 0 , , N 1 , in system (14) and (15) are coupled. Moreover, this system is overdetermined. Therefore, direct methods for solving this system are not applicable. To overcome this difficulty, we convert the above system into a minimization problem. For simplicity of notation, we denote by Q x : = q 0 ( x ) , . . . , q N 1 ( x ) . We approximate the solution of (14) and (15) by minimizing the following objective functional:
J ( Q ) = 1 2 k = 0 N 1 Ω · q k ( x ) + m = 0 N 1 n = 0 N 1 F k m n q m ( x ) · q n ( x ) 2 W 2 ( x ) d x ,
in the set
G : = { Q H 1 ( Ω ) : q k H 1 ( Ω ) R , q k ( x ) · ν = φ k ( x ) , x Ω } ,
where R is a given positive constant. In (16), W ( x ) is a Carleman weight function. The key idea of adding this weight function to the objective functional J is to make it convex in the set G. The choice of the weight function W ( x ) depends on the geometry of the domain Ω . We will discuss the choice of this function in our numerical examples.
To minimize the objective functional J in (16), we use a proximal block alternating minimization method. This method updates the unknown functions q k , k = 0 , , N 1 , alternatively; see, e.g., [30,31,32,33]. The algorithm is described as follows.
The sequence { Q ( i ) } obtained by the alternating minimization satisfies the following square summable property.
Lemma 1.
Let Q ( i ) : = ( q 0 ( i ) , , q N 1 ( i ) ) , i = 0 , 1 , , be obtained by Algorithm 1. Then:
i = 1 Q ( i ) Q ( i 1 ) H 1 ( Ω ) 2 : = i = 1 k = 0 N 1 q k ( i ) q k ( i 1 ) H 1 ( Ω ) 2 < .
Algorithm 1: Proximal block alternating minimization algorithm for minimizing the objective functional J ( Q ) .
  • Given an initial guess Q ( 0 ) : = ( q 0 ( 0 ) , q 1 ( 0 ) , , q N 1 ( 0 ) ) and a positive constant L.
  • Repeat the following steps until a stopping criterion is met:
    • For k = 0 , 1 , , N 1 :
           Find q k ( i ) as the minimizer of the objective functional:
J k ( q k ) : = 1 2 J ( q 0 ( i ) , , q k 1 ( i ) , q k , q k + 1 ( i 1 ) , , q N 1 ( i 1 ) ) + L 2 q k q k ( i 1 ) H 1 ( Ω ) 2 .
Proof. 
Since q k ( i ) is the minimizer of J k , it follows from (18) that J k ( q k ( i ) ) J k ( q k ( i 1 ) ) . This inequality implies that
J ( q 0 ( i ) , , q k 1 ( i ) , q k ( i ) , q k + 1 ( i 1 ) , , q N 1 ( i 1 ) ) + L 2 q k ( i ) q k ( i 1 ) H 1 ( Ω ) 2 J ( q 0 ( i ) , , q k 1 ( i ) , q k ( i 1 ) , q k + 1 ( i 1 ) , , q N 1 ( i 1 ) ) .
Hence,
L 2 q k ( i ) q k ( i 1 ) H 1 ( Ω ) 2 J ( q 0 ( i ) , , q k 1 ( i ) , q k ( i 1 ) , q k + 1 ( i 1 ) , , q N 1 ( i 1 ) ) J ( q 0 ( i ) , , q k 1 ( i ) , q k ( i ) , q k + 1 ( i 1 ) , , q N 1 ( i 1 ) ) .
Taking the sum of inequality (19) for k = 0 , 1 , , N 1 , we obtain
L 2 Q ( i ) Q ( i 1 ) H 1 ( Ω ) 2 J ( Q ( i 1 ) ) J ( Q ( i ) ) .
Thus,
L 2 i = 1 Q ( i ) Q ( i 1 ) H 1 ( Ω ) 2 J ( Q ( 0 ) ) lim n J ( Q ( n ) ) .
The proof is complete. □
Remark 2.
The global convergence of the alternating minimization algorithm has been proved for objective functionals which satisfy the so-called Kurdyka–Lojasiewicz (KL) property. One case in which this KL property is satisfied is when the objective functional is locally strongly convex; see, e.g., [34,35]. In the case d = 1 , as we will prove in Theorem 1, the objective functional J ( Q ) in (16) is strongly convex in an arbitrary given set. Therefore, the global convergence of the alternating minimization follows.

4. Discretized Objective Functional and Convexity in 1D

In this section, we consider the discretized version of the objective functional J ( Q ) in (16) in 1D. Our focus in this section is to prove the convexity of the discretized version of J in an arbitrary domain in an appropriate Hilbert space. Although results concerning the convexity of continuous Carleman weighted objective functionals associated with second-order elliptic systems have been reported [14], this appears to be the first time that the convexity for the discretized objective functional for the first-order system (14) is proved. The multi-dimensional case is still open.
In the following, let Ω = ( 0 , b ) . In this 1D problem, we denote the spatial variable by x instead of x . The objective functional J can be rewritten as
J ( Q ) = 1 2 k = 0 N 1 0 b q k ( x ) + m = 0 N 1 n = 0 N 1 F k m n q m ( x ) q n ( x ) 2 W 2 ( x ) d x ,
Consider a partition of the interval ( 0 , b ) into M sub-intervals by a uniform set of grid points 0 = x 0 < x 1 < < x M = b with grid size h : = b / M . Denote by q k i an approximate value of q k ( x i ) , i = 0 , , M . In the following we use Q to represent the matrix Q = ( q k i ) . The objective functional J is replaced by the following discrete version:
J h ( Q ) : = h 2 k = 0 N 1 i = 0 M 1 ( R k i ( Q ) 2 W i 2 ,
where W i = W ( x i ) and R k i ( Q ) is given by
R k i ( Q ) = q k i + 1 q k i h + m = 0 N 1 n = 0 N 1 F k m n q m i q n i .
To simplify the notation, we define the following bilinear function:
G k i ( Q , Q ˜ ) : = m = 0 N 1 n = 0 N 1 F k m n q m i q ˜ n i ,
for two matrices Q = ( q k i ) and Q ˜ = ( q ˜ k i ) .
The boundary condition (15) can be rewritten as
q k 0 = φ k ( 0 ) , q k M = φ k ( b ) .
The minimization problem now becomes a finite dimensional problem in which the unknown matrix Q belongs to the set M of real matrices of sizes N × ( M + 1 ) . In the following, we consider the standard Frobenius norm of matrices in M defined by
Q F : = k = 0 N 1 i = 0 M ( q k i ) 2 1 / 2 .
We also consider the H 1 norm in M :
Q H 1 2 : = k = 0 N 1 i = 0 M 1 q k i + 1 q k i h 2 + Q F 2 .
We also use , F to denote the Frobenius inner product of two matrices in M .
For each positive constant R, we define the set B ( Q ^ , R ) in M as
B ( Q ^ , R ) : = { Q M : Q Q ^ H 1 < R , q k 0 = φ k ( 0 ) , q k M = φ k ( b ) } .
where Q ^ : = ( q ^ k i ) , with q ^ k i : = φ k ( 0 ) + ( φ k ( b ) + φ k ( 0 ) ) i h / b , i = 0 , , M .
As we mentioned in Section 2, the choice of the Carleman weight function W ( x ) for the convexity of the objective functional J (or J h ) depends on the domain Ω . In this 1D case, we can choose W ( x ) = e λ x , where λ is a positive constant referred to as the Carleman weight coefficient. With this choice, we can prove that the objective functional J h is convex in the domain B ( Q ^ , R ) for any fixed radius R given that the Carleman weight coefficient λ chosen is large enough. To prove this property, we need the following estimate, which can be considered as a discrete Carleman estimate. Its proof is presented in Appendix A.
Lemma 2.
Let ξ i , i = 0 , , M , be scalar values in R with ξ 0 = 0 . Then, for W i : = e 2 i λ h , i = 0 , , M , λ > 0 , h > 0 , the following estimate holds:
i = 0 M 1 W i 2 ξ i + 1 ξ i h 2 e λ h 1 h 2 i = 0 M ( ξ i ) 2 W i 2 .
Now we are ready to state and prove the convexity of the objective functional J h .
Theorem 1.
Let R be an arbitrary positive number. Then, there exists a sufficiently large positive number λ 0 depending only on the parameters R , h , F k m n , Q ^ , N such that the objective functional J h ( Q ) defined by (23) is strongly convex on B ¯ ( Q ^ , R ) for all λ λ 0 . More precisely, there exists a constant C * = C * ( R , h , F k m n , Q ^ , N ) > 0 such that for arbitrary matrices Q ˜ , Q B ¯ ( Q ^ , R ) , the following estimate holds:
J h ( Q + δ Q ) J h ( Q ) J h ( Q ) ( δ Q ) C * δ Q H 1 2 , for all λ λ 0 ,
where J h is the gradient of J h .
Proof. 
It follows from (24) that
R k i ( Q + δ Q ) = q k i + 1 q k i h + δ q k i + 1 δ q k i h + G k i ( Q + δ Q , Q + δ Q ) .
Since G k i is bilinear, it follows that
G k i ( Q + δ Q , Q + δ Q ) = G k i ( Q , Q ) + G k i ( Q , δ Q ) + G k i ( δ Q , Q ) + G k i ( δ Q , δ Q ) .
Hence,
R k i ( Q + δ Q ) = R k i ( Q ) + R k i ( δ Q ) + L k i ( Q , δ Q ) ,
where L k i ( Q , δ Q ) = G k i ( Q , δ Q ) + G k i ( δ Q , Q ) . Thus,
[ R k i ( Q + δ Q ) ] 2 [ R k i ( Q ) ] 2 = [ R k i ( Q + δ Q ) + R k i ( Q ) ] [ R k i ( Q + δ Q ) R k i ( Q ) ] = [ 2 R k i ( Q ) + R k i ( δ Q ) + L k i ( Q , δ Q ) ] [ R k i ( δ Q ) + L k i ( Q , δ Q ) ] = 2 R k i ( Q ) [ δ q k i + 1 δ q k i h + L k i ( Q , δ Q ) ] + 2 R k i ( Q ) G k i ( δ Q , δ Q ) + [ R k i ( δ Q ) + L k i ( Q , δ Q ) ] 2 .
Substituting this equality into (23), we obtain
J h ( Q + δ Q ) J h ( Q ) = h 2 k = 0 N 1 i = 0 M 1 [ R k i ( Q + δ Q ) ] 2 [ R k i ( Q ) ] 2 W i 2 = h 2 k = 0 N 1 i = 0 M 1 2 W i 2 R k i ( Q ) [ δ q k i + 1 δ q k i h + L k i ( Q , δ Q ) ] + h 2 k = 0 N 1 i = 0 M 1 2 W i 2 R k i ( Q ) G k i ( δ Q , δ Q ) + h 2 k = 0 N 1 i = 0 M 1 W i 2 [ R k i ( δ Q ) + L k i ( Q , δ Q ) ] 2 .
Since the last two terms of (30) are quadratic functions of δ Q and the first term on the right-hand side is a linear function of δ Q , we have
J h ( Q ) δ Q = h k = 0 N 1 i = 0 M 1 W i 2 R k i ( Q ) [ δ q k i + 1 δ q k i h + L k i ( Q , δ Q ) ] .
Hence,
J h ( Q + δ Q ) J h ( Q ) J h ( Q ) δ Q = h k = 0 N 1 i = 0 M 1 W i 2 R k i ( Q ) G k i ( δ Q , δ Q ) + h 2 k = 0 N 1 i = 0 M 1 W i 2 [ R k i ( δ Q ) + L k i ( Q , δ Q ) ] 2 .
Now we estimate the right-hand side of (31). First, we have
[ R k i ( δ Q ) + L k i ( Q , δ Q ) ] 2 = δ q k i + 1 δ q k i h + G k i ( δ Q , δ Q ) + L k i ( Q , δ Q ) 2 = δ q k i + 1 δ q k i h 2 + 2 δ q k i + 1 δ q k i h G k i ( δ Q , δ Q ) + L k i ( Q , δ Q ) + G k i ( δ Q , δ Q ) + L k i ( Q , δ Q ) 2 .
Using the Cauchy–Schwarz inequality, we have
2 δ q k i + 1 δ q k i h G k i ( δ Q , δ Q ) + L k i ( Q , δ Q ) 1 2 δ q k i + 1 δ q k i h 2 2 G k i ( δ Q , δ Q ) + L k i ( Q , δ Q ) 2 .
Substituting this inequality into (32), we obtain
[ R k i ( δ Q ) + L k i ( Q , δ Q ) ] 2 1 2 δ q k i + 1 δ q k i h 2 G k i ( δ Q , δ Q ) + L k i ( Q , δ Q ) 2 .
From (31) and (33) it follows that
J h ( Q + δ Q ) J h ( Q ) J h ( Q ) δ Q h k = 0 N 1 i = 0 M 1 W i 2 R k i ( Q ) G k i ( δ Q , δ Q ) h 2 k = 0 N 1 i = 0 M 1 W i 2 G k i ( δ Q , δ Q ) + L k i ( Q , δ Q ) 2 + h 4 k = 0 N 1 i = 0 M 1 W i 2 δ q k i + 1 δ q k i h 2 .
Since Q, Q + δ Q are in B ( Q ^ , R ) , there exists a constant C > 0 depending only on R, Q ^ , h, and F f m n such that
h k = 0 N 1 i = 0 M 1 W i 2 R k i ( Q ) G k i ( δ Q , δ Q ) h 2 k = 0 N 1 i = 0 M 1 W i 2 G k i ( δ Q , δ Q ) + L k i ( Q , δ Q ) 2 C h k = 0 N 1 i = 0 M 1 W i 2 ( δ q k i ) 2 .
It follows from Lemma 2 that
h 4 k = 0 N 1 i = 0 M 1 W i 2 δ q k i + 1 δ q k i h 2 h 8 k = 0 N 1 i = 0 M 1 W i 2 δ q k i + 1 δ q k i h 2 + h 8 e λ h 1 h 2 k = 0 N 1 i = 0 M W i 2 ( δ q k i ) 2 .
Note that W i 2 e 2 λ b . From (34)–(36) we obtain
J h ( Q + δ Q ) J h ( Q ) J h ( Q ) δ Q h 8 e 2 λ b k = 0 N 1 i = 0 M 1 δ q k i + 1 δ q k i h 2 + h 8 e λ h 1 h 2 C h e 2 λ b δ Q F 2 .
It is clear that there exists λ 0 > 0 such that h 8 e λ h 1 h 2 C h > 0 for all λ λ 0 . Hence, (29) follows from (37) with
C * = h e 2 λ b min 1 8 , 1 8 e λ h 1 h 2 C .
The proof is complete. □

5. Numerical Implementation and Examples

In this section we describe details of our numerical implementation of the proposed algorithm.

5.1. Generating Simulated Data

In this work, we tested the proposed algorithm using simulated data. To generate these data, we solved the forward problem (2) and (3) in a bounded interval ( l 1 , l 2 ) such that l 1 < 0 < b < l 2 . The source function was chosen as f ( x , t ) = δ ( x x 0 ) f ( t ) , where x 0 l 1 . The total wave function u ( x , t ) was rewritten as u ( x , t ) = u i ( x , t ) + u s ( x , t ) , where u i is the incident wave and u s is the scattered wave. The incident wave u i was given by the following formula:
u i ( x , t ) = f ( t | x x 0 | ) , t | x x 0 | , 0 , 0 t < | x x 0 | .
As in [14], we also used the absorbing boundary condition to approximate the boundary conditions at the endpoints x = l 1 and x = l 2 . The forward problem was then rewritten as follows:
c ( x ) u t t s ( x , t ) u x x s ( x , t ) = [ 1 c ( x ) ] u t t i ( x , t ) , ( x , t ) ( l 1 , l 2 ) × ( 0 , T ) ,
u s ( x , 0 ) = u t s ( x , 0 ) = 0 , x ( l 1 , l 2 ) ,
  u x s ( l 1 , t ) = u t s ( l 1 , t ) , u x s ( l 2 , t ) = u t s ( l 2 , t ) , t ( 0 , T ) .
In the numerical examples presented below, we chose l 1 = x 0 = 0.2 , b = 1 , l 2 = 1.5 , and T = 2 . The waveform f ( t ) of the incident wave was chosen to be
f ( t ) = A ( t 0.2 ) e ω 2 ( t 0.2 ) 2 ,
where ω = 30 and A = 2 ω e 1 / 2 . The constant A was used as the normalization factor. Problem (38)–(40) was solved by an explicit finite difference scheme with a uniform grid in both the x and t directions with step sizes of Δ x = 0.005 and Δ t = 0.001 . This results in 141 grid points in space and 2001 points in time.
In order to simulate noisy measurements, we added additive noise of 5 % (in the L 2 norm) to the simulated data.

5.2. Numerical Examples of the CIP

In the following examples, the parameters were chosen as follows. The pseudofrequencies were s 4 , 15 with step size Δ s = 0.05 . The number of Laguerre’s functions was N = 11 . The Carleman weight coefficient λ was chosen to be λ = 3 . The proximal coefficient L was L = 10 4 . These parameters were chosen numerically for the best reconstruction of one simulated dataset, then they were fixed for all other examples. We have observed in numerical tests that a larger number of Laguerre’s functions do not necessarily result in substantially more accurate reconstruction results. The possible reason is due to the fact that the objective functional is quite insensitive to changes in the high-frequency modes near the minimizer.
We also remark that the data associated with large pseudofrequencies are dominated by that at low pseudofrequencies. This is due to the exponential decaying nature of the Laplace transform. Therefore, the interval of the pseudofrequencies chosen should not be too large.
Example 1.
In the first example, we considered a smooth coefficient c ( x ) in the following form:
c ( x ) = 1 + e 100 ( x 0.4 ) 2 , 0 x 1 .
The reconstructed coefficient is depicted in Figure 1 together with the exact coefficient. We can see that the algorithm was able to reconstruct the coefficient quite accurately. The largest error is near x = 0.6 , which is near the “back” of the scattered. We note that the source illuminates the scatterer from the left. As a result, the measured signal at x = 0 is stronger than that at x = 1 . That made the reconstruction on the left side of the scatterer easier to reconstruct.
Example 2.
In the second example, we tested the algorithm for a more challenging smooth function c ( x ) . More precisely, c ( x ) is a combination of two Gaussian peaks. Its formula is given by
c ( x ) = 1 + e 400 ( x 0.3 ) 2 + 2 e 400 ( x 0.7 ) 2 , 0 x 1 .
The reconstructed coefficient is depicted in Figure 2. Although the reconstruction is not as accurate as in Example 1, it still captured the main peaks of the coefficient function. The largest peak was reconstructed with a relative error of approximately 17% (2.5 vs. 3) and the second peak was reconstructed with a relative error of approximately 15% (1.7 vs. 2).
Example 3.
Finally, we tested the algorithm for a piecewise constant coefficient with large and small peaks. The coefficient was given as
c ( x ) = 1 + 10 χ [ 0.2 , 0.4 ] + 3 χ [ 0.7 , 0.8 ] ,
where χ [ a , b ] is the characteristic function of interval [ a , b ] .
The reconstructed coefficient is depicted in Figure 3. As we can see, the algorithm was still able to reconstruct the peaks of the coefficient, with a relative error of approximately 10%, even though the quantitative values are not so accurate. This is a very challenging example which mimics the case when a “weak” object and a “strong” object are in the same medium. In this case, it is usually difficult to reconstruct the weaker target. To see the effect of the number of Laguerre’s functions on the reconstruction accuracy in this challenging example, we also tested using 20 Laguerre’s functions but the result (not shown in the figure) was not really improved.
Concerning computational cost, the algorithm was set to run for 2000 iterations, which took approximately 5 min on a Lenovo Thinkpad X1 Yoga laptop computer with an Intel i7 1.8 GHz processor with four cores (but only one was used in the computation).

6. Conclusions

In this paper we combined the Carleman weighted convexification method with an alternating minimization method for solving a 1D inverse medium scattering problem. The convexity of the Carleman weighted objective functional was proved in the discrete setting, which was actually used in the optimization algorithm. The global convergence of the alternating minimization was shown. Numerical examples indicated that the proposed algorithm was able to provide a good solution without requiring a good first guess. If more accurate results are expected, we can use locally convergent gradient-based methods with the result of the proposed algorithm as an initial guess to refine the reconstruction.

Funding

This research received no external funding.

Data Availability Statement

No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Acknowledgments

The author is grateful to the constructive comments and suggestions from the referees.

Conflicts of Interest

The author declares no conflict of interest.

Appendix A

Proof of Lemma 2.
The proof follows the same approach as in [36], Lemma 3.1. Let η i = ξ i e i λ h . Then ξ i = η i e i λ h . Hence,
W i 2 ξ i + 1 ξ i h 2 = e 2 i λ h η i + 1 e ( i + 1 ) λ h η i e i λ h h 2 = η i + 1 e λ h η i h 2 = 1 h 2 η i + 1 ( e λ h 1 ) + ( η i + 1 η i ) 2 = e λ h 1 h 2 ( η i + 1 ) 2 + ( e λ h 1 ) h 2 2 η i + 1 ( η i + 1 η i ) + 1 h 2 ( η i + 1 η i ) 2 e λ h 1 h 2 ( η i + 1 ) 2 + ( e λ h 1 ) h 2 2 η i + 1 ( η i + 1 η i ) .
Taking the sum for i = 0 , , M 1 , we obtain
i = 0 M 1 W i 2 ξ i + 1 ξ i h 2 e λ h 1 h 2 i = 0 M 1 ( η i + 1 ) 2 + ( e λ h 1 ) h 2 i = 0 M 1 2 η i + 1 ( η i + 1 η i ) .
Since η 0 = ξ 0 = 0 , it follows that
i = 0 M 1 2 η i + 1 ( η i + 1 η i ) = i = 0 M 1 ( η i η i + 1 ) 2 + ( η M ) 2 0 .
Hence,
i = 0 M 1 W i 2 ξ i + 1 ξ i h 2 e λ h 1 h 2 i = 0 M 1 ( η i + 1 ) 2 = e λ h 1 h 2 i = 0 M ( ξ i ) 2 W i 2 .
The proof is complete. □

References

  1. Colton, D.; Kress, R. Inverse Acoustic and Electromagnetic Scattering Theory, 3rd ed.; Springer: New York, NY, USA, 2013. [Google Scholar]
  2. Ammari, H. An Introduction to Mathematics of Emerging Biomedical Imaging; Springer: Berlin/Heidelberg, Germany, 2008. [Google Scholar]
  3. Baum, C. Detection and Identification of Visually Obscured Targets; Taylor and Francis: Philadelphia, PA, USA, 1998. [Google Scholar]
  4. Buchanan, J.L.; Gilbert, R.P.; Wirgin, A.; Xu, Y.S. Marine Acoustics: Direct and Inverse Problems; Society for Industrial and Applied Mathematics: Philadelphia, PA, USA, 2004. [Google Scholar]
  5. Daniels, D.J. Ground Penetrating Radar; The Institute of Electrical Engineers: London, UK, 2004. [Google Scholar]
  6. Pastorino, M. Microwave Imaging; John Wiley & Sons Inc.: Hoboken, NJ, USA, 2010. [Google Scholar]
  7. Ireland, D.; Bialkowski, M. Microwave heat imaging for stroke detection. Prog. Electromagn. Res. M 2011, 21, 163–175. [Google Scholar] [CrossRef] [Green Version]
  8. Nikolova, N.K. Microwave imaging for breast cancer. IEEE Microw. Mag. 2011, 12, 78–94. [Google Scholar] [CrossRef]
  9. Persson, M.; Fhager, A.; Trefna, H.D.; Yu, Y.; McKelvey, T.; Pegenius, G.; Karlsson, J.E.; Elam, M. Microwave-Based Stroke Diagnosis Making Global Prehospital Thrombolytic Treatment Possible. IEEE Trans. Biomed. Eng. 2014, 61, 2806–2817. [Google Scholar] [CrossRef] [Green Version]
  10. Scapaticci, R.; Tobon, J.; Bellizzi, G.; Vipiana, F.; Crocco, L. Design and Numerical Characterization of a Low-Complexity Microwave Device for Brain Stroke Monitoring. IEEE Trans. Antennas Propag. 2018, 66, 7328–7338. [Google Scholar] [CrossRef]
  11. Tournier, P.H.; Aliferis, I.; Bonazzoli, M.; de Buhan, M.; Darbas, M.; Dolean, V.; Hecht, F.; Jolivet, P.; Kanfoud, I.E.; Migliaccio, C.; et al. Microwave tomographic imaging of cerebrovascular accidents by using high-performance computing. Parallel Comput. 2019, 85, 88–97. [Google Scholar] [CrossRef]
  12. Yilmaz, O. Seismic Data Imaging; Society of Exploration Geophysicists: Tulsa, OK, USA, 1987. [Google Scholar]
  13. Klibanov, M.V.; Ioussoupova, O.V. Uniform strict convexity of a cost functional for three-dimensional inverse scattering problem. SIAM J. Math. Anal. 1995, 26, 147–179. [Google Scholar] [CrossRef]
  14. Klibanov, M.V.; Thành, N.T. Recovering dielectric constants of explosives via a globally strictly convex cost functional. SIAM J. Appl. Math. 2015, 75, 518–537. [Google Scholar] [CrossRef] [Green Version]
  15. Klibanov, M.V.; Kolesov, A.E. Convexification of a 3-D coefficient inverse scattering problem. Comput. Math. Appl. 2019, 77, 1681–1702. [Google Scholar] [CrossRef] [Green Version]
  16. Klibanov, M.V.; Kolesov, A.E.; Nguyen, L.; Sullivan, A. Globally strictly convex cost functional for a 1-D inverse medium scattering problem with experimental data. SIAM J. Appl. Math. 2017, 77, 1733–1755. [Google Scholar] [CrossRef] [Green Version]
  17. Klibanov, M.V.; Nguyen, D.L.; Nguyen, L.H.; Liu, H. A globally convergent numerical method for a 3D coefficient inverse problem with a single measurement of multi-frequency data. Inverse Probl. Imaging 2018, 12, 493–523. [Google Scholar] [CrossRef] [Green Version]
  18. Kolesov, A.E.; Klibanov, M.V.; Nguyen, L.H.; Nguyen, D.L.; Thành, N.T. Single measurement experimental data for an inverse medium problem inverted by a multi-frequency globally convergent numerical method. Appl. Numer. Math. 2017, 120, 176–196. [Google Scholar] [CrossRef] [Green Version]
  19. Nguyen, D.L.; Klibanov, M.V.; Nguyen, L.H.; Fiddy, M.A. Imaging of buried objects from multi-frequency experimental data using a globally convergent inversion method. J. Inverse Ill-Posed Probl. 2018, 26, 501–522. [Google Scholar] [CrossRef] [Green Version]
  20. Chen, X.; Wei, Z.; Li, M.; Rocca, P. A review of deep learning approaches for inverse scattering problems. Prog. Electromagn. Res. 2020, 167, 67–81. [Google Scholar] [CrossRef]
  21. Fan, Y.; Ying, L. Solving Inverse Wave Scattering with Deep Learning. arXiv 2019, arXiv:1911.13202. [Google Scholar] [CrossRef]
  22. Khoo, Y.; Ying, L. SwitchNet: A Neural Network Model for Forward and Inverse Scattering Problems. SIAM J. Sci. Comput. 2019, 41, A3182–A3201. [Google Scholar] [CrossRef] [Green Version]
  23. Xu, K.; Darve, E. Physics Constrained Learning for Data-driven Inverse Modeling from Sparse Observations. J. Comput. Phys. 2022, 453, 110938. [Google Scholar] [CrossRef]
  24. Wei, Z.; Chen, X. Deep-Learning Schemes for Full-Wave Nonlinear Inverse Scattering Problems. IEEE Trans. Geosci. Remote Sens. 2019, 57, 1849–1860. [Google Scholar] [CrossRef]
  25. Ambrosanio, M.; Franceschini, S.; Pascazio, V.; Baselice, F. An End-to-End Deep Learning Approach for Quantitative Microwave Breast Imaging in Real-Time Applications. Bioengineering 2022, 9, 651. [Google Scholar] [CrossRef]
  26. Dachena, C.; Fedeli, A.; Fanti, A.; Lodi, M.B.; Fumera, G.; Pastorino, M.; Randazzo, A. Microwave Medical Imaging of the Human Neck using a Neural-Networks-Based Inversion Procedure: A Phantom Study. In Proceedings of the 2023 17th European Conference on Antennas and Propagation (EuCAP), Florence, Italy, 26–31 March 2023; pp. 1–5. [Google Scholar]
  27. Salucci, M.; Polo, A.; Vrba, J. Multi-Step Learning-by-Examples Strategy for Real-Time Brain Stroke Microwave Scattering Data Inversion. Electronics 2021, 10, 95. [Google Scholar] [CrossRef]
  28. Thành, N.T.; Beilina, L.; Klibanov, M.V.; Fiddy, M.A. Imaging of Buried Objects from Experimental Backscattering Time-Dependent Measurements Using a Globally Convergent Inverse Algorithm. SIAM J. Imaging Sci. 2015, 8, 757–786. [Google Scholar] [CrossRef] [Green Version]
  29. Abramowitz, M.; Stegun, I.A. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables; National Bureau of Standards Applied Mathematics Series; U.S. Government Printing Office: Washington, DC, USA, 1964; Volume 55. [Google Scholar]
  30. Beck, A. On the convergence of alternating minimization for convex programming with applications to iteratively reweighted least squares and decomposition schemes. SIAM J. Optim. 2015, 25, 185–209. [Google Scholar] [CrossRef] [Green Version]
  31. Grippo, L.; Sciandrone, M. Globally convergent block-coordinate techniques for unconstrained optimization. Optim. Methods Softw. 1999, 10, 587–637. [Google Scholar] [CrossRef]
  32. Tupitsa, N.; Dvurechensky, P.; Gasnikov, A.; Guminov, S. Alternating minimization methods for strongly convex optimization. J. Inverse Ill-Posed Probl. 2021, 29, 721–739. [Google Scholar] [CrossRef]
  33. Wright, S.J. Coordinate descent algorithms. Math. Program. 2015, 151, 3–34. [Google Scholar] [CrossRef]
  34. Attouch, H.; Bolte, J.; Redont, P.; Soubeyran, A. Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Łojasiewicz inequality. Math. Oper. Res. 2010, 35, 438–457. [Google Scholar] [CrossRef] [Green Version]
  35. Xu, Y.; Yin, W. A Block Coordinate Descent Method for Regularized Multiconvex Optimization with Applications to Nonnegative Tensor Factorization and Completion. SIAM J. Imaging Sci. 2013, 6, 1758–1789. [Google Scholar] [CrossRef] [Green Version]
  36. Thành, N.T. Convexity of a discrete Carleman weighted objective functional in an inverse medium scattering problem. J. Inv. Ill-Posed Probl. 2020, 30, 485–493. [Google Scholar] [CrossRef]
Figure 1. Reconstruction of the coefficient c ( x ) in Example 1. The data were measured at x = 0 and x = 1 .
Figure 1. Reconstruction of the coefficient c ( x ) in Example 1. The data were measured at x = 0 and x = 1 .
Axioms 12 00642 g001
Figure 2. Reconstruction of the coefficient in Example 2. The data were measured at x = 0 and x = 1 .
Figure 2. Reconstruction of the coefficient in Example 2. The data were measured at x = 0 and x = 1 .
Axioms 12 00642 g002
Figure 3. Reconstruction of the coefficient in Example 3. The data were measured at x = 0 and x = 1 .
Figure 3. Reconstruction of the coefficient in Example 3. The data were measured at x = 0 and x = 1 .
Axioms 12 00642 g003
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Thành, N.T. Using Alternating Minimization and Convexified Carleman Weighted Objective Functional for a Time-Domain Inverse Scattering Problem. Axioms 2023, 12, 642. https://doi.org/10.3390/axioms12070642

AMA Style

Thành NT. Using Alternating Minimization and Convexified Carleman Weighted Objective Functional for a Time-Domain Inverse Scattering Problem. Axioms. 2023; 12(7):642. https://doi.org/10.3390/axioms12070642

Chicago/Turabian Style

Thành, Nguyen Trung. 2023. "Using Alternating Minimization and Convexified Carleman Weighted Objective Functional for a Time-Domain Inverse Scattering Problem" Axioms 12, no. 7: 642. https://doi.org/10.3390/axioms12070642

APA Style

Thành, N. T. (2023). Using Alternating Minimization and Convexified Carleman Weighted Objective Functional for a Time-Domain Inverse Scattering Problem. Axioms, 12(7), 642. https://doi.org/10.3390/axioms12070642

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop