# Soft Shrinkage Thresholding Algorithm for Nonlinear Microwave Imaging

Hidayet Zaimaga, Marc Lambert

► To cite this version:

Hidayet Zaimaga, Marc Lambert. Soft Shrinkage Thresholding Algorithm for Nonlinear Mi-crowave Imaging. Journal of Physics: Conference Series, 2016, 756 (1), pp.012011. 10.1088/1742-6596/756/1/012011. [https://doi.org/10.1088/1742-6596/756/1/012011](https://doi.org/10.1088/1742-6596/756/1/012011). [hal-01390779](https://HAL.inhal-01390779)

HAL Id: hal-01390779

https://centralesupelec.hal.science/hal-01390779

Submitted on 16 Jul 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Soft Shrinkage Thresholding Algorithm for Nonlinear Microwave Imaging

To cite this article: Hidayet Zaimaga and Marc Lambert 2016 *J. Phys.: Conf. Ser.* **756** 012011

View the [article online](#) for updates and enhancements.

## Related content

- Development of a thresholding algorithm for calcium classification at multiple CT energies

LY. Ng, M. Alssabbagh, A. A. Tajuddin et al.

- Sparsity regularization in inverse problems

Bangti Jin, Peter Maaß and Otmar Scherzer

- Note

Bumseok Namgung, Peng Kai Ong, Yun Hui Wong et al.

## Recent citations

- Target Tracking via Particle Filter and Convolutional Network

Hongxia Chu et al

<figure><img src="image_2.png" /></figure>

IOP ebooks™

Bringing together innovative digital publishing with leading authors from the global scientific community.

Start exploring the collection-download the first chapter of every title for free.
# Soft Shrinkage Thresholding Algorithm for Nonlinear Microwave Imaging

Hidayet Zaimaga¹, Marc Lambert²

¹ Laboratoire des signaux et systèmes-CNRS UMR8506, CentraleSupélec-Univ.Paris-Sud, Université Paris-Saclay

² GeePs — Group of electrical engineering - Paris, UMR CNRS 8507, CentraleSupélec, Univ. Paris-Sud, Université Paris-Saclay, Sorbonne Universités, UPMC Univ Paris 06

E-mail: hidayet.zaimaga@12s.centralesupelec.fr

**Abstract.** In this paper, we analyze a sparse nonlinear inverse scattering problem arising in microwave imaging and numerically solved it for retrieving dielectric contrast from measured fields. In sparsity reconstruction, contrast profiles are a priori assumed to be sparse with respect to a certain base. We proposed an approach which is motivated by a Tikhonov functional incorporating a sparsity promoting $l_1$-penalty term. The proposed iterative algorithm of soft shrinkage type enforces the sparsity constraint at each nonlinear iteration. The scheme produces sharp and good reconstruction of dielectric profiles in sparse domains by adapting Barzilai and Borwein (BB) step size selection criteria and positivity by maintaining its convergence during the reconstruction.

## 1. Introduction

Development of efficient reconstruction methods and techniques that exploit sparseness regularized formulations have been widely emerged for solving inverse electromagnetic scattering problems in recent years. High demand of such methods in various applications such as material characterization, subsurface prospecting, remote sensing, and non-destructive testing and evaluation [1,2] enforces the importance and the need of effective and accurate methods. Inverse electromagnetic scattering problems reconstruct material properties such as permittivity and conductivity in an unknown region from measured electromagnetic fields. However, implementation of such stable, reliable, and efficient reconstruction algorithms is challenging because of the nonlinearity of the scattering equations and ill-posedness of the problem [2–4].

Several approaches have been proposed in order to overcome these issues. Global optimization tools, multi-step information retrieval techniques, and qualitative methods have been introduced to alleviate the non-linearity or its effects. Moreover, first order approximations such as diffraction tomography, Born and Rytov approximations have been proven in order to linearize the problem [2, 3, 5]. On the other hand, innovative sparseness-regularized formulations have recently emerged as an effective recipe to overcome the non-uniqueness and/or numerical instability of the inversion process [6, 7]. The reason behind this is that many images have sparse representations with respect to their expansion basis in the wavelet domain and this yields new developing approaches that minimize the cost functions with zeroth/first norm penalty terms

<figure><img src="image_3.png" /></figure>
using highly effective iterative shrinkage thresholding algorithms. Such an increased interest is proven by number of publications in wide domains [4, 6–10].

The aim of the present work is to develop a novel sparse reconstruction algorithm and demonstrate its performance in microwave imaging. In the proposed work sparse reconstruction algorithm of iterative soft shrinkage type is studied. Sparsity constraint is directly applied to the problem of reconstructing the complex internal dielectric properties of an object based on knowledge of the external scattered field which is generated by the interaction between the object and a known incident field. The nonlinear optimization problem is solved by iterative algorithm of soft shrinkage in order to enforce the sparsity constraint by a soft thresholding function. We exploit an adaptive step size selection by BB rule and a priori constraint called as positivity in order to have improved reconstructions. The proposed method provides necessary and sufficient conditions to yield well-posedness and convergence [11–13].

## 2. Sparse Problem Statement and Discretization

Consider a 2D configuration with transverse magnetic polarization (TM) case where the object under investigation illuminated by a given source numbered as $i, i = 1, ..., Ns$ as in fig. 1. Let $E^{\text{inc}}$ be an incident electric field which is generated by a source and polarized along the z-axis with an implicit time factor $\exp(-i\omega t)$. The object is considered in an investigation domain $D$ and the different media are characterized by their propagation constant $k(\mathbf{r})$ such that $k(\mathbf{r})^2 = \omega^2\varepsilon_0\varepsilon_r(\mathbf{r})\mu_0 + i\omega\mu_0\sigma(\mathbf{r})$, where $\omega$ is the angular frequency, $\varepsilon_0$ and $\mu_0$ are the permittivity and the permeability of the air respectively, $\varepsilon_r(\mathbf{r})$ and $\sigma(\mathbf{r})$ are the relative permittivity and conductivity of the medium as $\mathbf{r} \in D$ is an observation point. The dielectric properties of $D$ are described by the inhomogeneous contrast function, $\chi(\mathbf{r})$, which is defined for non-magnetic area such as $\chi(\mathbf{r}) = (k(\mathbf{r})^2 - k_B^2)$ where $k_B$ is the propagation constant of the embedding medium.

A known incident field interacts with the scatterer yielding to a total field which is the sum of the incident and scattered fields. Assuming no magnetic media and considering the boundary and radiation conditions we can obtain two coupled contrast source integral equations by applying Green’s theorem to Helmholtz wave equations [1, 5] as following

$$E^{\text{tot}}(\mathbf{r}) = E^{\text{inc}}(\mathbf{r}) + \int_{D} G(\mathbf{r}, \mathbf{r}')J(\mathbf{r}')d\mathbf{r}' \quad \forall \mathbf{r} \in D \quad (1)$$

$$E^{\text{scatt}}(\mathbf{r}) = \int_{D} G(\mathbf{r}, \mathbf{r}')J(\mathbf{r}')d\mathbf{r}' \quad \forall \mathbf{r} \in L \quad (2)$$

where $G(\mathbf{r}, \mathbf{r}') = \frac{i}{4}H_0^{(1)}(k_B\|\mathbf{r} - \mathbf{r}'\|)$, $H_0^{(1)}$ is the 0-th order Hankel function of the first kind and $E^{\text{scatt}}$ is the scattered field. We also defined the contrast source induced within the object by the incident wave such as $J(\mathbf{r}') = \chi(\mathbf{r})E^{\text{tot}}(\mathbf{r}')$, where $E^{\text{tot}}$ being the total field in the object. Solving the equations (1) and (2) the unknown contrast function $\chi(\mathbf{r})$ can be determined. Following this, the scattering equations are discretized with the point-matching Method of Moments as in [5] by considering pixel basis functions

$$\psi_k(\mathbf{r}) = \begin{cases} 1 & \mathbf{r} \in D_k, \\ 0 & \text{otherwise,} \end{cases} \quad (3)$$

$D_k$ being the $k$-th pixel and the unknown contrast defined as $\chi(\mathbf{r}) = \sum_{k=1}^{N} \chi_k \psi_k(\mathbf{r})$, $N$ being the number of pixels.
### 3. Nonlinear Optimization

Soft shrinkage is an approach which minimizes a nonlinear Tikhonov functional with sparsity promoting penalty term. The algorithm is based on the iterated soft shrinkage approach originated for linear operators in the work [6]. A generalization to nonlinear inverse problems has been studied in [4, 14].

The algorithm performs a gradient descent step which involves the adjoint gradient of the cost function with a step size $\kappa$, and then a shrinkage step. The latter enforces the sparsity of the reconstruction by setting the small coefficients to zero. The solution of the inverse problem can be obtained by minimizing the cost function which is of the form

$$F(\chi) = \frac{1}{2} \|\xi(\chi) - E^{\text{scatt}}\|^2 + \alpha \|\chi\|_{l_1} \quad (4)$$

whereas $\xi(\chi) = G(\mathbf{r}, \mathbf{r}')\chi(\mathbf{r}')[I - G(\mathbf{r}, \mathbf{r}')\chi(\mathbf{r})]^{-1}E^{\text{inc}}$. The $l_1$ penalty can promote a-priori knowledge of the sparse representation.

Iterative soft shrinkage has the form as (4) where $\xi: X \to Y$ is counted to be a bounded and nonlinear operator with respect to unknown contrast. The algorithm is started by choosing an initial guess $\chi^1$, and continues the iteration as

$$\chi^{j+1} = S_{\alpha} \left( \chi^j - \kappa \xi^{*} (\chi^j) \left[ \xi(\chi^j) - E^{\text{scatt}} \right] \right) \quad (5)$$

where $\xi'(\chi)$ is the gradient of nonlinear function $\xi(\chi)$, and $\xi^{*}(\chi)$ is the adjoint operator. $S_{\alpha}$ is the soft shrinkage operator defined componentwise by [6]

$$(S_{\alpha}(\chi))_i = \begin{cases} (|\chi_i| - \alpha)\text{sign}(\chi_i), & \text{if } |\chi_i| > \alpha \\ 0, & \text{otherwise.} \end{cases} \quad (6)$$

It truncates small values to zero and shrinks large values. The term $\xi^{*}(\chi^j)[\xi(\chi^j) - E^{\text{scatt}}]$ is the gradient of the discrepancy $\frac{1}{2}\|\xi(\chi) - E^{\text{scatt}}\|^2$. The algorithm which has been computed, is as following:

<caption>**Algorithm 1** Steepest descent reconstruction algorithm with sparsity constraint</caption>

1: Initialize $\chi^1$ and $\alpha$
2: **for** $j = 1, \dots, J$ **do**
3: Solve the direct problem $E^{\text{scatt}}(\chi^j)$
4: Compute the gradient $K'(\chi^j) = \nabla_{\chi} K(\chi)|_{\chi^j}$
5: Determine the step size $\kappa_j$
6: Update inhomogeneity by $\chi^{j+1} = \chi^j - \kappa_j K'(\chi^j)$
7: Threshold $\chi^{j+1}$ by $S_{\alpha}(\chi^{j+1})$
8: check stopping criterion
9: **end for**
10: **output** approximate the minimizer of (4).

### 3.1. Step Size Selection

Selecting the proper step size adaptively is important as it can increase the convergence speed. Therefore, the step size $\kappa$ can be determined in order to accelerate the algorithm. The iterative
<figure><img src="image_4.png" /><figcaption>Figure 1: Measured configuration of actual permittivity profile and source-receiver locations on $(x(m), y(m))$ axis.</figcaption></figure>

soft shrinkage algorithm with a fixed step size favors the classical Landweber method. Thus, the motivation for increasing the rate of convergence is the comparison with the classical Landweber iteration whose slow convergence results from using a constant step size which is very small [4, 10]. Hence, we select the step size in a way to increase the convergence speed where we only consider the steepest descent operation $\chi^j - \kappa K'(\chi^j)$ of the algorithm. The selection is done by the two-point rule of Barzilai and Borwein which calculates the step size as [15]

$$ \kappa_j = \frac{\langle \chi^j - \chi^{j-1}, K'(\chi^j) - K'(\chi^{j-1}) \rangle}{\langle K'(\chi^j) - K'(\chi^{j-1}), K'(\chi^j) - K'(\chi^{j-1}) \rangle} \quad (7) $$

### 3.2. Priori Constraints

Imposing a-priori constraints can sometimes improve the quality of solutions to the inverse problems in a great portion [16]. Towards this end, non-negativity is important in applications like imaging [6, 16]. We know that in order to have a physical solution the unknowns we are dealing with should have constraints, the latter being $\epsilon_r(\mathbf{r}) \ge 1$ and $\sigma(\mathbf{r}) \ge 0$. However, in the general case (as embedded obstacle in half-space) those constraints do not impose non-negativity to the real and imaginary parts of the contrast function. We can consider the constraint by “Projection” where at each iteration the following projection $\epsilon_r(\mathbf{r}) = \max(\epsilon_r(\mathbf{r}), 1)$ and $\sigma(\mathbf{r}) = \max(\sigma(\mathbf{r}), 0)$ is applied.

## 4. Numerical Results

In the following, the basis functions $\psi_k$ have been chosen as pixel-based. In our example (see fig. 1), we have five objects of physical characteristics described in Tab. 1. Investigation area $D$ is a $l = 6 \times \lambda$-sided square. The discretization size is $n \times n = 80 \times 80$ for the forward problem and $n \times n = 30 \times 30$ for the inverse problem. The number of transmitters and receivers located around the investigation area are 36 evenly distributed on a circle of radius $r = 7$ m. The frequency of the transmitters is 200 MHz. The measured field samples are generated by adding 10 dB white Gaussian noise.

Let us defined the relative error norm as a comparison criteria expressed as

$$ \epsilon_r^{\text{rec}} = \frac{\|\epsilon_r^{\text{rec}} - \epsilon_r\|_2}{\|\epsilon_r\|_2}, \quad \sigma^{\text{err}} = \frac{\|\sigma^{\text{rec}} - \sigma\|_2}{\|\sigma\|_2}, \quad (8) $$

where $\epsilon_r^{\text{rec}}$ and $\sigma^{\text{err}}$ are the reconstructed permittivity and conductivity respectively and $\epsilon_r$ and $\sigma$ the exact ones. The minimum of cost function being the one reached at the end of the process for each $\alpha$.

One of the key point of such an inversion is the choice of the regularization parameter $\alpha$ in (6). Different tests have been performed in order to evaluate the sensitivity of the choice of
<table><caption>Table 1: Description of the five obstacles, $x$, $y$ being the coordinate of the center of the obstacle (in m), $L_x$, $L_y$ its lengths (in m) and $\varepsilon_r$ and $\sigma$ its relative permittivity and conductivity (the latter in S m$^{-1}$)</caption><thead><tr><th>#</th><th>x</th><th>y</th><th>L<sub>x</sub></th><th>L<sub>y</sub></th><th>ε<sub>r</sub></th><th>σ</th></tr></thead><tbody><tr><td>1</td><td>-1.5</td><td>-1.5</td><td>0.5</td><td>0.5</td><td>1.5</td><td>0.0022</td></tr><tr><td>2</td><td>-0.33</td><td>0.44</td><td>0.5</td><td>0.5</td><td>2.25</td><td>0</td></tr><tr><td>3</td><td>1.5</td><td>-1.5</td><td>0.5</td><td>0.33</td><td>2</td><td>0.00445</td></tr><tr><td>4</td><td>0</td><td>1.5</td><td>0.5</td><td>0.33</td><td>2</td><td>0.02225</td></tr><tr><td>5</td><td>1</td><td>1.33</td><td>0.5</td><td>0.5</td><td>1</td><td>0.02225</td></tr></tbody></table>

regularization on the solution. In our case it can be observed that the higher the value we choose for $\alpha$ the sharper the reconstructions are. However, the choice of $\alpha$ is not arbitrary as can be seen in Fig.2. We can get good reconstruction by choosing $\alpha$ within the range of $\alpha = 1.5 \times 10^{-3}$ and $\alpha = 1 \times 10^{-2}$ for which $\varepsilon_r^{\text{err}}$ and $\sigma^{\text{err}}$ are the smallest as shown fig. 2b and fig. 2c respectively whereas those errors get bigger for larger $\alpha$.

Moreover, addition of projection to the sparsity reconstruction method improve the quality of the reconstruction even more. Fig. 2 shows us that we can choose the regularization parameter in a wider range when we add projection to our algorithm.

As a complementary validation Figs. 3 and 4 present a comparison of the map of the permittivity and conductivity respectively for various inversion parameters. The influence of the $\alpha$ parameter without any projection constraint onto the permittivity can be shown by comparing fig. 3b and fig. 3c whereas the influence of the projection constraint is presented fig. 3d. Clearer image is obtained by addition of projection constraint with an appropriate $\alpha$.

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Figure 2: Error on the cost function (as defined in eq. 4), $\varepsilon_r^{\text{err}}$, and $\sigma^{\text{err}}$ (as defined in eq. 8) as a function of $\alpha$. •-with positivity via projection and ×-without positivity</figcaption></figure>

## 5. Conclusion

Soft iterative thresholding is used to solve the 2-D electromagnetic inverse scattering problem based on $l_1$ norm penalty term. Retrieval of arbitrary complex shaped targets from simulated data shows that this approach is effective and gives us sharper reconstructions by adopting wavelet basis function. Taking advantage of BB method for step size selection and adding projection enhance the effectiveness of proposed method. Even if choosing the regularization parameter is challenging, addition of projection can lead us to choose it in a wider range possibilities.

Even though the results are promising and advantageous, investigations in experimental settings need to be done to show the feasibility of the algorithm whereas this would be our next step. The development of nonlinear inversion with sparsity regularization for higher permittivity
<figure><img src="image_8.png" /><figcaption>Figure 3: Permittivity $\varepsilon_r$ obtained by using sparsity and positivity with 10 dB noise data.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 4: Conductivity $\sigma$ obtained by using sparsity and positivity with 10 dB noise data.</figcaption></figure>

values will be studied. Smoothness constraint for the gradient and projection on gradient (as mentioned in section 3.2) will be added to the proposed algorithm in order to analyze the improvement in the reconstruction quality.

## Acknowledgment

This work is supported by DIGITEO France through the SIRENA project (2014–2017) under the “Call for Chairs 2014”.

## References

[1] Pastorino M 2010 *Microwave imaging* vol 208 (John Wiley & Sons)

[2] Poli L, Oliveri G and Massa A 2012 *IEEE Trans. Antennas Propagat.* **60** 2865–2879

[3] Di Benedetto F, Estatico C, Nagy J G and Pastorino M 2009 *Electron. Trans. Numer. Anal.* **33** 105–125

[4] Ramlau R and Teschke G 2006 *Numer. Math.* **104** 177–203

[5] Richmond J H 1965 *IEEE Trans. Antennas Propagat.* **13** 334–341

[6] Daubechies I, Defrise M and De Mol C 2003 *arXiv preprint math/0307152*

[7] Jin B and Maass P 2012 *Inverse Prob.* **28** 123001

[8] Desmal A and Bagci H 2015 *Prog. Electromagn. Res.*

[9] Grasmair M, Haltmeier M and Scherzer O 2008 *Inverse Prob.* **24** 055020

[10] Hanke M, Neubauer A and Scherzer O 1995 *Numer. Math.* **72** 21–37

[11] Bagci H and Desmal A 2014 *IEEE Trans. Antennas Propagat.* **62** 3878–3884

[12] Gehre M, Kluth T, Lipponen A, Jin B, Seppänen A, Kaipio J P and Maass P 2012 *J. Comput. Appl. Math.* **236** 2126–2136

[13] Jin B, Khan T and Maass P 2012 *Int. J. Numer. Methods Eng.* **89** 337–353

[14] Bonesky T, Bredies K, Lorenz D A and Maass P 2007 *Inverse Prob.* **23** 2041

[15] Barzilai J and Borwein J M 1988 *IMA J. Numer. Anal.* **8** 141–148

[16] Vogel C R 2002 *Computational methods for inverse problems* vol 23 (SIAM)