An adaptive approach for the detection of contrast targets for the through-wall imaging

Abstract

Through-wall imaging is capable of detecting various living and non-living things behind the wall. The characteristics of the wall under the investigation, amount of clutter and noise govern the quality and reliability of the image as well as the detection ability of the targets using through the wall imaging system. The characteristics of the wall are not known prior, in the literature only the intensity profile is investigated for the unknown wall characteristics using a single dielectric target and the effect of the wall characteristics on the contrast imaging and impact on time or frequency domain features are not investigated. The target with less dielectric is having less reflectivity; hence its detection in the presence of a high reflective target and a noisy environment becomes difficult. In this paper, to enhance the detection ability of the imaging system attenuation constant (α) of the wall is estimated with the proposed wall parameter estimation methods and used as a normalizing factor. To achieve effective beamforming different weighting strategies are developed and the obtained images are compared with the traditional beamforming. Furthermore, a novel approach to finding the effective rank in the low-rank estimation using a statistical model and multi-objective genetic algorithm is proposed for de-noising.

Introduction

In civilian and military applications, through-wall imaging (TWI) is the popular technology in rescue missions and target detection behind the wall. In locating and detecting the target, characteristics of the wall such as thickness and dielectric constant play an important role [1]. Incorrect estimation of the wall characteristics can smear and blur the image quality and also shift the target position. There are several studies available in the literature where the characteristics of the wall are known in advance [2–4], but in real-time scenarios characteristics of the wall are not known and hence autofocusing techniques become important to improve the image quality. Autofocusing for TWI application using various methods is proposed in [2–9], but these techniques do not estimate the wall parameters and hence, do not compensate for the attenuation factor for the wall. To our knowledge many studies are performed with the same dielectric constant target and the effect of the presence of contrast targets in microwave imaging is missed out. To improve the intensity profile and to reduce the spreading of the target in the image, standardized moments are used in [1], unlike this paper here standardized moments are used to determine the parameters of the wall such as thickness and permittivity. Once the parameters of the wall are estimated then attenuation due to the wall can be calculated and compensated. The image pixel value depends upon the complex composite signal received from the target and neighboring scene and it consists of focusing delay, intensity profile, and weights assigned to each pixel. To our knowledge very little work is available in the literature on the contrast target imaging and decisions making of these weighing strategies, hence weighing strategies based on the swath of the antenna is proposed in this paper.

TWI data are generally corrupted by the clutter and noise; to remove system inherent noise external calibration and some pre-processing steps are required [3]. In this paper, to improve the contrast target detection hybrid technique is proposed based on opt-shrink and wavelet thresholding. Peak signal to noise ratio (PSNR) can be increased by the proper separation between signal and noise, to recover signal matrix (X) from the noisy data matrix (Y) low-rank approximation is a popular technique among signal and image processing researchers [10]. In singular value decomposition (SVD), received signal (Y) is represented by separate right and left singular matrixes as, $Y = USV^T = \sum\nolimits_{i = 1}^{\min ( m, n) } {\lambda _iu_iv_i^T }$, and λ_i the numbers of singular values (SVs) are truncated by hard thresholding or soft thresholding. In a hard thresholding number of SVs can be chosen by cross-validation but this approach is not reliable [11, 12]. In [13] it is mentioned that in the case of TWI first eigenvalue from SVD corresponds to the wall and only the second eigenvalue corresponds to the target. In [14] authors show that apart from the second eigenvalue the target subspace is spread over all the eigenvalues in the noise subspace; hence another efficient approach is required to exploit the signal from the noise. Opt-Shrink algorithm is proposed in [15] to approximate the SVs optimally by weighing the singular vector from the measurement matrix (S), this algorithm estimates the weights optimally in the form of shrinkage and thresholding. Two threshold methods are discussed in the literature namely soft threshold also called shrinkage function T _s(x) = sgn(x) × max(|x| − T,  0) and hard threshold i.e. T _h(x) = x × 1{|x| > T} otherwise T _h(x) = 0. Soft threshold shrinks SVs using the function ${\mathop{\lambda}\limits^{\frown}}_i = \lambda _i( 1-{\tau / {\lambda _i}}) +$, where i is the variable for each threshold value. Candès et al. [14] have proposed a soft threshold τ > 0 by Stein's unbiased risk estimate (SURE). Dong et al. [16] extended the work done by Candès et al. using the principle of wavelet Bayes-Shrink, which estimates SVs based on local signal and overall noise variance. Xie et al. [17] proposed that noise energy is distributed over all eigenvalues; hence the relative significance of the SV in SVD is needed to be estimated for effective de-noising.

In our proposed approach, the relative significance of SVs with wall parameters is established to recover images better in terms of PSNR using a multi-objective genetic algorithm (GA). To improve the intensity profile proper threshold for the microwave image is required. In [18] parameter “n” is defined in the function ${\mathop{\lambda}\limits^{\frown}}_i = \lambda _i( 1-{{\tau ^n} / {\lambda _i^n }}) +$ to choose between the hard and soft threshold, if “n” is ∞ then hard thresholding is selected and if it is 1 then soft thresholding is selected. A data-driven soft threshold for image de-noising is proposed in [19], the framework for the proposed threshold is Bayesian and it works on wavelet coefficients which are normally distributed (μ = 0,  σ = 1). In the above framework, the threshold is in closed form and also adaptive since it depends upon the data-driven estimates of the parameter.

The objectives of this work are to determine the characteristics of the wall under investigation and to detect the contrast target behind the wall by removing the noise in the low-contrast part of the image. The novelty of the work is (a) the parameters for the wall under investigation are estimated using higher-order standardized moments and estimation is improved further using GA. The different weights based on the swath of the antenna and the distance of the target from the transceiver are proposed to improve the quality and reliability of the image. (b) The hybrid technique based on Opt-Shrink and Bayes-Shrink algorithms is developed to reduce weak noise-like patterns in low-contrast parts of the image. The limitation of this work is that the study for the characterization of the wall is carried out for the homogenous wall if the wall is inhomogenous computational complexity can be more.

This paper is organized as follows: Section “TWI experimental setup and pre-processing” illustrates the TWI imaging setup used in the experimental work and reviews de-noising related work available in the literature. Section “Proposed wall characterization and focusing weights assignment” demonstrates the wall parameters estimation and the effect of the different weights applied in TWI beamforming. Section “Proposed technique for wall parameter estimation and image de-noising” describes the proposed hybrid technique for wall removal and de-noising using data-driven Bayes-Shrink in the wavelet domain. Section “Conclusion” concludes the work.

TWI experimental setup and pre-processing

In our experimental work, raw data are acquired by scanning the whole wall, TWI setup for which is shown in Fig. 1. It consists of a single antenna and transceiver; the technique used is stepped frequency continuous wave (SFCW). The raw data are collected for the targets of different dielectrics such as metal (ɛ = ∞), wood (ɛ = 2.9), and Teflon (ɛ = 2.1) behind the wall at different distances. The transceiver is used to generate 201 numbers of points in the frequency range 1–3 GHz and the reflection coefficient S ₁₁ is measured.

Fig. 1. TWI setup and geometry.

The combinations and arrangement of the targets behind the wall are described in Table 1 and the transceiver system parameters are given in Table 2.

Table 1. Different combinations and arrangements for the targets behind the wall

Table 2. Transceiver system parameters

The antenna system is moved in the horizontal direction to scan the whole wall. “L” locations x_ak (k = 1…L) are chosen at equal distances of S = 5 cm; the imaging geometry for TWI is shown in Fig. 2. Consider an antenna is placed at L positions {x_ak = (x_ak, (y_off + d))}, where k varies from 1 to L and y_off is the standoff distance between the antenna and a uniform wall of the thickness (d) and dielectric constant (ɛ_r) located in the x–y plane.

Fig. 2. 2-D geometry for TWI scenario.

In our experimental work, the bandwidth used for imaging is 2 GHz; at this bandwidth, the pixel size in the cross-range is x = 2.5 cm and down-range is y = 3.75 cm, interested readers may refer [20] for more details about cross-range and down-range resolution. The imaging region lies along the positive y-axis and can be divided into M*N pixels behind the wall as shown in Fig. 3, along with cross-range and down-range. The wideband signal s(t) is assumed and wideband delay and sum beamforming (DSBF) [21] is used for generating the image from the raw data; this involves electronically focusing the beams across all pixels of the scene. Each transmitter transmits a pulse s(t) with relative time delays so that all pulses transmitted reach a particular location simultaneously; for example, the qth pixel located at x_q = (x_i, y_j) with i varying from 1 to M and j varying from 1 to N. These focusing delays are applied on transmission and reception and adjusted to sweep the beams across all voxels in the image. To control the shape and side-lobe structure of the beams, additional weights further described in the Section “Proposed wall characterization and focusing weights assignment” can be applied.

Fig. 3. Formation of imaging grid.

A single point target located at position x_q = (x_i, y_j), the output of the delay-and-sum beam-former corresponding to the qth voxel at x_q is given in [21–23]:

(1)$$z_q( t ) = \mathop \sum \limits_{k = 1}^L w_ka( {x_p} ) e^{-\alpha ( {2^\ast l_{kq, wall}} ) }s( {t-\tau_{k, p} + \tau_{k, q}} ) $$

where a(x_p) is the target reflectivity and w_k is the weights applied on transmit and receive antenna (both are same). α is the attenuation constant of the wall, l_kq, _wall is the distance traveled by the wave inside the wall on transmission and reception. τ_k,p is the propagation delay encountered by the signal as it propagates from the kth transmitter to the target and back to the kth receiver, and τ_k,q is the aggregate focusing delay applied to the kth transmitter and the output of the kth receiver. The propagation delay and focusing delay is calculated for i = p and for i = q respectively:

(2)$$\tau _{k, i} = \displaystyle{{2^\ast l_{ki, air, 1}} \over c} + \displaystyle{{2^\ast l_{ki, wall}} \over v} + \displaystyle{{2^\ast l_{ki, air, 2}} \over c}$$

where v = c/$\sqrt {\varepsilon _r}$, i = p for the target and i = q for the qth pixel, c is the velocity of light in free space, and the subscripts “air,1” “wall,” “air,2” denote the distances before, through, and beyond the wall respectively to x_q for i = q and x_p for i = p. The complex amplitude image value I(x_q) corresponding to the qth pixel is then computed by applying a matched filter to z_q(t), matched to s(t), and sampling the filtered data, which is given by

(3)$$I( {x_q} ) = z_q( t ) \ast h( t ) \vert _{t = 0}$$

(4)$$I( {x_q} ) = \sum\limits_{k = 1}^L {w_k\,} a( {x_p} ) e^{-\alpha ( {2^\ast l_{kq, wall}} ) }s( {t-\tau_{k, p} + \tau_{k, q}} )^ \ast h( t ) \vert _{t = 0}$$

where h(t) = s*(−t) is the impulse response of the matched filter, the superscript * denotes complex conjugate, and “*” denotes convolution operator. To generate the complete image of the scene, the process, described in equations (1)–(4), is performed for all P pixels in the region of interest:

(5)$$I( {x_q} ) = \sum\limits_{k = 1}^L {\sum\limits_{\,p = 1}^p {w_k\,} } a( {x_p} ) e^{-\alpha ( {2^\ast l_{kq, wall}} ) }s( {t-\tau_{k, p} + \tau_{k, q}} )^ \ast h( t ) \vert _{t = 0}$$

For a scene consisting of P point targets, the complex amplitude pixel value I(x_q) can be obtained by the superposition of the target reflections. The equivalent frequency-domain representation of the qth voxel value can be expressed as

(6)$$\eqalign{I( {x_q} ) & = \sum\limits_{k = 1}^L {\sum\limits_{\,p = 1}^p {w_k} } \cdot a( {x_p} ) e^{-\alpha ( {2^\ast l_{kq, wall}} ) }\displaystyle{1 \over {2\pi }}\cr& \int_{-\infty }^\infty {\vert S( \omega ) } \vert ^2\exp ( {\,j\omega ( {t-\tau_{k, p} + \tau_{k, q}} ) } )} $$

where, S(ω) is the Fourier transform of the transmitted signal s(t). The exact knowledge about the dielectric constant and wall thickness allows us to calculate focusing delays. Further pre- and post-processing to remove system inherent noise and conversion from the frequency domain to time domain and from the time domain to spatial domain is carried out using different steps given in [3].

Calculation for wave propagation distances

First, we have to consider that the wave travels from the kth transmit position x_ak to the qth voxel positioned at x_q. The distances l_kq, _{air, 1}, l_kq, _wall, and l_kq, _{air, 2} are expressed as derived in [21, 22]:

(7)$$l_{kq, air, 1} = \displaystyle{{y_{off}} \over {{\rm cos}( {\theta_{kq}} ) }}$$

(8)$$l_{kq, wall} = \displaystyle{d \over {{\rm cos}( {\varphi_{kq}} ) }}$$

(9)$$l_{kq, air, 2} = \displaystyle{{y_j} \over {{\rm cos}( {\theta_{kq}} ) }}$$

If φ_kq is the refraction angle and θ_kq is the incidence angle, they are evaluated by solving

(10)$$n_1\sin ( {\theta_{mq}} ) = n_2\sin ( {\varphi_{mq}} ) $$

If the exact focusing distance is not used, then there is a shift in the actual position of the target, which causes huge errors in estimating the exact position and size of the targets. This happens because refraction phenomena occur through the wall. Although exact calculation requires some computational burden, it cannot be avoided if accurate results are required.

Autofocusing is a popular technique to eliminate the image distortion caused by errors in wall parameter estimation. In [1] different contrast measures are proposed to investigate the quality and reliability of the TWI image. These contrast measures generally do not give satisfactory results so higher-order standardized moments are proposed, both conventional and higher-order standardized moments are summarized here.

1. (1) Normalized sum of image intensity:

   (11)$$C_1 = \displaystyle{{\mathop \sum \nolimits_{q = 1}^Q {\vert {I( {x_q} ) } \vert }^2} \over {{\left({\mathop \sum \nolimits_{q = 1}^Q \vert {I( {x_q} ) } \vert } \right)}^2}}$$

where I(x _q) = intensity of qth pixel and Q = total number of pixels in the image.

1. (2) Normalized sum of squared intensity:

   (12)$$C_2 = \displaystyle{{\mathop \sum \nolimits_{q = 1}^Q {\vert {I( {x_q} ) } \vert }^4} \over {{\left({\mathop \sum \nolimits_{q = 1}^Q \vert {I( {x_q} ) } \vert } \right)}^4}}$$

2. (3) Negative of image entropy:

   (13)$$C_3 = \mathop \sum \limits_{q = 1}^Q ss( {x_q} ) {\rm ln}( {ss( {x_q} ) } ) $$

where $ss( {x_q} ) = {{{\vert {I( {x_q} ) } \vert }^2} \over {\mathop \sum \nolimits_{q = 1}^Q \vert {I( {x_q} ) } \vert }}$.

1. (4) Ratio of standard deviation to mean amplitude:

   (14)$$C_4 = \displaystyle{{\sqrt {\mathop \sum \nolimits_{q = 1}^Q {\left[{\vert {I( {x_q} ) } \vert -( 1/Q) \mathop \sum \nolimits_{q = 1}^Q \vert {I( {x_q} ) } \vert } \right]}^2} } \over {\mathop \sum \nolimits_{q = 1}^Q \vert {I( {x_q} ) } \vert }}$$

2. (5) Higher-order metrics:

   (15)$$\gamma _n = \displaystyle{{\mathop \sum \nolimits_{q = 1}^Q {( {P( {x_q} ) -\hat{\mu }} ) }^n} \over {( {Q-1} ) {\hat{\sigma }}^n}}$$
   (16)$$Skew[ P ] = \gamma _3 = \gamma _n = \displaystyle{{\mathop \sum \nolimits_{q = 1}^Q {( {P( {x_q} ) -\hat{\mu }} ) }^3} \over {( {Q-1} ) {\hat{\sigma }}^3}}$$
   (17)$$Kurt[ P ] = \gamma _4-3 = \displaystyle{{\mathop \sum \nolimits_{q = 1}^Q {( {P( {x_q} ) -\hat{\mu }} ) }^4} \over {( {Q-1} ) {\hat{\sigma }}^4}}-3$$

where γ denotes the higher order metrics, σ is the standard deviation, and μ is the mean, while $\hat{\mu }$ is the estimated mean. In this paper, unlike [1] we have used the above measures to estimate the thickness and dielectric constant of the wall when the either-or parameter is considered [24]. In the category of higher-order metrics, Kurtosis is the fourth-order moment and gives greater details about the probability distribution of the data hence it is considered in the estimation for wall parameters.

The raw image is obtained by DSBF; generally, this image is corrupted by clutter and noise, due to this clutter and noise detection of weak targets such as Teflon or wood in the presence of strong target such as metal becomes difficult. The major clutter contribution in the TWI image is from the wall, clutter removal techniques based on eigenvalue decomposition such as SVD is the popular technique. SVD truncates the SVs by hard thresholding; authors in [20] show that SVD cannot detect weak targets when eigenvalue truncates with hard thresholding; hence in this paper a hybrid technique based on Opt-Shrink and data-driven adaptive thresholding is proposed.

Opt-Shrink shrinkage and thresholding algorithm

Opt-Shrink algorithm has been introduced recently in [15] to truncate eigenvalues from the noisy measurement matrix by soft thresholding; the principle of the random matrix theory is used in this paper for characterization of a large matrix. The Opt-Shrink algorithm optimally weights the coefficients obtained for the large class noise model such as i.i.d. Gaussian noise case. D-transform is the natural choice in the Opt-Shrink algorithm; it is an analog form of the log-Fourier transforms which gives a relation between the distribution of independent and individual singular matrices. If the pre-processed raw image matrix is R; different steps for the implementation of Opt-Shrink algorithm are as follows:

1. (1) Estimate the effective rank for low-rank signal matrix

2. (2) Compute SVD for R by $\hat{R} = \sum\nolimits_{i = 1}^q {{\hat{\lambda }}_i{\hat{u}}_i{\hat{v}}_i^T }$

3. (3) Compute $\sum\nolimits_{\hat{r}} { = diag( {\hat{\lambda }}_{r + 1} \ldots {\hat{\lambda }}_q) \in R^{( {n-\hat{r}} ) \times ( {m-\hat{r}} ) }}$

4. (4) Compute D-transform for $\hat{D}( \hat{\lambda }_i, \;\sum\nolimits_{\hat{r}} {} )$ and $\widehat{{{D}^{\prime}}}( \hat{\lambda }_i, \;\sum\nolimits_{\hat{r}} {} )$

   $$\hat{D}( {z, \;x} ) = 1/nTr( {z{( {z^2I-xx^H} ) }^{{-}1}} ) \cdot 1/mTr( {z{( {z^2I-x^Hx} ) }^{{-}1}} ) $$
   $$\eqalign{{{\hat D}^{\prime}}( {z, \;x} ) & = 1/nTr((z{(z^{2}I-xx^{H})}^{-1}).1/mTr({-2z{(z^{2}I-x^{H}x)}^{-2}}) \cr & + {(z^{2}I-x^{H}x)}^{-1}) + 1/mTr\;( z{(z^{2}I-x^{H}x)}^{-1}) \cr & \cdot 1/nTr(-2z^{2}{(z^{2}I-xx^{H})}^{-2} + {(z^{2}I-x^{H}x)}^{-1}}$$

5. (5) Compute $\omega _{i, \hat{r}}^{opt} = -2\,\hat{D}( \hat{\lambda }_i, \;\sum\nolimits_{\hat{r}} {} ) /\widehat{{{D}^{\prime}}}( \hat{\lambda }_i, \;\sum\nolimits_{\hat{r}} {} )$

6. (6) Evaluate $\widehat{{S_{opt}}} = \sum\nolimits_{i = 1}^{\hat{r}} {\omega _{i, \hat{r}}^{opt} {\hat{u}}_i{\hat{v}}_i^T }$

To estimate the effective rank, a novel approach based on multi-objective GA optimization is proposed in this paper. Opt-Shrink algorithm is principally based on SVD in terms of least squares. As stated in [25] mean square error (MSE) for these methods is larger than the Bayesian estimator, also these methods give higher outlier values and generate weak noise-like patterns in low-contrast areas of the image [10]. To remove these weak noise patterns effective thresholding is required.

Wavelet-based Bayesian data adaptive thresholding

Consider the received signal at the antenna is $\{ {S_{ij}, \;\,i, \;\,j = 1 \hbox{ to }N} \} \exists N = i^2$ where i is the integer. If this received signal is corrupted by i.i.d. Gaussian noise (0,  σ ²), our goal is to remove the noise and estimate $\widehat{{S_{ij}}}$ which minimizes MSE:

(18)$$MSE( {\hat{S}} ) = 1/N^2\mathop \sum \limits_{i, j = 1}^N ( {\widehat{{S_{ij}}}-S_{ij}} ) $$

A data-driven Bayesian technique is proposed in [19] to remove the weak noise from low-contrast areas of the image. In this paper, wavelet de-noising is carried out for $\widehat{{S_{opt}}}$ received by Opt-Shrink after modifying for TWI application. Let G is the received signal corrupted by noise, S is the original signal, and € is the additive noise the boldfaced letters represent matrix under consideration:

(19)$$G = S + \unicode{x20AC}$$

Let Z = WG, where Z is the wavelet transform output, G is the wavelet coefficients derived from 2D-orthogonal dyadic wavelet W; reference [26] may be referred by interested readers for more details about the dyadic wavelet transform. Wavelet filters give estimate $\hat{Z}$ for each coefficient Z, this estimate for the de-noised wavelet transform output is obtained by $\hat{S} = W^{{-}1}{\boldsymbol X}$ where W ⁻¹ is the inverse wavelet transform and X is the input matrix. Wavelet thresholding keeps low-resolution coefficients intact while removing noise from the detail coefficients. This is the main reason to remove weak noise from low-contrast areas wavelet de-noising is selected, which can be helpful in the detection of low-dielectric and contrast targets in TWI successfully.

Proposed wall characterization and focusing weights assignment

The raw data are collected for materials having different dielectric constants such as Teflon, wood, and metal at different distances from the wall having 2ⁿ − 1 = 7 combinations where n is the number of targets, using our TWI system. Raw images developed using the DSBF algorithm for Teflon, metal, and contrast targets such as wood and metal are shown in Figs 4(a), 4(b), and 4(c), respectively, for illustrative purposes; it can be seen from the raw image that the edge details are completely lost and actual length of the target cannot be identified hence post-processing is necessary to improve the image quality and is described further.

Fig. 4. Raw images using DSBF: (a) target ID 03: Teflon, (b) target ID 01: metal, and (c) target ID 7: metal and wood.

Wall effect minimization

Accurate wall attenuation is always difficult to predict without effective modeling due to a broad variety of building materials. Properties for the wall depend upon the architectural need, the electromagnetic properties such as complex propagation constant (γ), the transmission coefficient (T), and reflection coefficient (Γ) show large variation with different building materials. Electrical characteristics (μ, ɛ, σ) show varied and random effects with the variation of building material hence becomes very difficult to predict. Complex propagation constant $\gamma = \sqrt {\mu _0\mu _r( \varepsilon -j\sigma /\omega } ) = \alpha + j\beta$ where α is the attenuation constant (Np/m), β = phase shift (rad/m), σ is the conductivity; σ = ∞ for a perfect conductor and σ = 0 for a perfect dielectric. Most of the time building materials are lossy and having finite ɛ and σ. The ratio of two (ɛ/σ) called relaxation time (τ′) and if σ ≫ ωɛ or if ωτ ≪ 1 then the material is a good conductor. In the case of lossless material, there will only be phase shift while for lossy material amplitude accompanies by the phase shift. In this paper, we concentrate on the attenuation constant (α) as it impacts intensity profile of the image majorly and velocity correction is done using estimated permittivity value for the wall.

In [24], the attenuation of different types of concrete and brick walls is estimated in the frequency range 1–5 GHz. In equation (5) there is a term $e^{-\alpha ( {2^\ast l_{kq, wall}} ) }$ related to attenuation due to wall and to nullify the effect of this factor we have to multiply I(x _q) with $e^{\alpha ( {2^\ast l_{kq, wall}} ) }$. The attenuation corrected image amplitude is given by final imaging equation:

(20)$$\eqalign{I( {x_q} ) & = I{\rm ^{\prime}}( {x_q} ) .\;e^{\alpha ( {2{^\ast} l_{kq, wall}} ) } {I}^{\prime}( {x_q} ) \cr& = \sum\limits_{k = 1}^L {\sum\limits_{\,p = 1}^p {w_k} } \cdot a( {x_p} ) \displaystyle{1 \over {2\pi }}\int_{-\infty }^\infty {\vert S( \omega ) } \vert ^2\exp ( {\,j\omega ( {t-\tau_{k, p} + \tau_{k, q}} ) } ) } $$

Correct estimation of the wall characteristic is important to determine the attenuation constant (α) and to develop a quality image using the TWI system. When various contrast measures mentioned in the Section “TWI experimental setup and pre-processing” are calculated for metal target and Teflon target different relations are obtained between the actual value and estimated value defined as estimation error (E_r) for the dielectric constant and thickness of the wall. First, we consider the thickness of the wall is known i.e. (d = 14 cm) which is a normal brick wall used in the experimental work. The relations obtained for the dielectric constant of the wall are shown in Figs 5 and 6, the contrast measures attained minimum or maximum value within some range depending on the type and number of targets.

Fig. 5. Variation of contrast measures with a dielectric constant for metal target.

Fig. 6. Variation of contrast measures with a dielectric constant for Teflon target.

The corresponding value of the dielectric constant further considered for estimating the thickness of the wall. In this case, the contrast measures attained either minimum or maximum value in the range 12–15 cm and this can be observed in Figs 7 and 8.

Fig. 7. Variation of contrast measures with wall thickness for metal target.

Fig. 8. Variation of contrast measures with wall thickness for Teflon target.

It can be observed in Figs 5–8 that the estimation error (E_r) for kurtosis at different values of wall thickness and dielectric constant attained some minimum or maximum value at different wall parameters depending upon the reflectivity of the target placed behind the wall.

The intensity profile is responsible for the variation in the error values as shown in Fig. 9. The attenuation factor for the wall is estimated after estimating the wall parameters by observing Figs 5–8. The estimation error governs the PSNR for the developed microwave image. To improve the wall parameter prediction, we further investigate the effect of the above estimation error on the received PSNR values with different combinations for the dielectric constant and thickness of the wall. PSNR is defined as

(21)$$\eqalign{PSNR = & 10\,{\rm lo}{\rm g}_{10}( 1/MSE) \cr MSE = & 1/M^\ast N( im-\widehat{{im) }}^2} $$

where im = raw image, $\widehat{{im}}$ = estimated image, M is the number of horizontal scanning points, and N is the number of vertical scanning points.

Fig. 9. Intensity profile for different targets behind the wall.

It can be observed in Table 3 that PSNR = 15.5694 is optimum for the contrast target at wall dielectric (ɛ _r = 5.3) and thickness (d = 13.5 cm). The attenuation constant (α) for the wall parameter values is 0.54695 from [24]. This value of α can be put in equation (20) to compensate for the attenuation due to the wall, processed images after amplitude attenuation correction are shown in Fig. 10.

Fig. 10. Processed TWI raw images for (a) target ID 3: Teflon and (b) target ID 1: metal after wall amplitude attenuation correction.

Table 3. PSNR values when the thickness of the wall (d) or dielectric constant (ɛ_r) of the wall is known

Application of weights (w_k)

Earlier in the section, Fig. 4 shows raw images developed by scanning the whole wall using M number of horizontal and N number of vertical points hence 2D data of dimension MN is obtained. Raw images developed using all MN data produce smooth edges, image saturation, and consume a lot of time for the data collection. To reduce the computation time as well as to produce an image with a minimum error weighing system is proposed. To overcome the effect of unnecessary rounding effect instead of assigning uniform weights w_k = 1, we assign weights based on the swath of the antenna which is calculated by the following equation and performance is compared using parameters like target size, shape, and PSNR:

(22)$$D_{swath} = 2{^\ast} [ {y_{off}\cdot {\rm tan}( {\theta_{kq}} ) + D_{wall}\cdot {\rm tan}( {\varphi_{kq}} ) + y_j{\rm tan}( {\theta_{kq}} ) } ] $$

where y _off = distance of antenna from wall, D _wall = the thickness of wall, y _j = distance of pixel for which weight needs to be estimated, θ _kq = incidence angle = beam-width/2, φ _kq = refraction angle through wall where, θ _kq and φ _kq are related by Snell's law:

$$\sqrt {\varepsilon _{air}}{\hskip2pt \ast} \sin ( {\theta_{kq}} ) = \sqrt {\varepsilon _{wall}} {\hskip2pt \ast} \sin ( {\varphi_{kq}} ) $$

Figure 11 shows the assignment of the weights for each pixel coming within the swath of the antenna; weights for each pixel can be governed by the following relations.

1. (i) Linear weights: In this scheme, when we are calculating the intensity corresponding to an image pixel, we give maximum weight to the nearest antenna present in front of the pixel and linearly decrease weights to the antennas further away from the pixel in all directions. It means we keep on decreasing the weights as we take readings from the antenna above and below the pixel as well as toward the right or left of the pixel.

Fig. 11. Weight assignment strategy.

For the computation of these weights, we assign linear values to only those antenna positions which are within the swath. So now, if there are N_t antennas with w _k = 1, then

(23)$$wt( i ) = 1-\displaystyle{{\vert {mid-i} \vert \cdot p} \over {mid-1}}$$

where mid = (N_t + 1)/2; p = slope parameter to decide the extent of variation, it is calculated using the Pythagorean theorem, where the vertical distance is the d (target distance from the antenna) and the horizontal distance is S = 5 cm (antenna movement for scanning) hence the run: d ² + S ² = slope length².

• (ii) Gaussian weights: In this scheme, the distribution of the weights is done like the above-described procedure but instead of linearly decreasing the weights, they are decreased according to the Gaussian distribution i.e. in addition to finding the swath and the exact number of antennas lying within the swath, we give more weight to the antenna closest to the desired pixel and less weight to the farthest (according to Gaussian distribution). A simple Gaussian distribution is given by

  (24)$$\eqalign{\,f( {x{\rm \vert }\mu , \;\sigma^2} ) & = \displaystyle{1 \over {\sqrt {2\pi \sigma ^2} }}e^{-{( {x-\mu } ) }^2/2\sigma ^2} \cr& wt( {x_i} ) = f( x ) \vert _{x = x_i}} $$

where i varies from 1 to N_t scanning antenna locations.

1. (iii) Euclidean weights: In this scheme, the weight distribution is done by considering a straight-line distance between the two points in Euclidean space. If the antenna is placed at location x_k(x, y) = x_k(x_i, y_j) and desired pixel is at location q = x_q(x_j, y_j) in Fig. 12 then the Euclidean distance between location p and q is given by the shortest distance between the two points and maximum weight is assigned to the minimum Euclidean distance among the number of points present in the pixel:

   (25)$$d( {\,p, \;q} ) = d( {q, \;p} ) = \sqrt {\mathop \sum \limits_{i = 1}^{N_t} {( {q_i-p_i} ) }^2} $$

Fig. 12. Geometry for different weighting strategies: (a) linear, (b) Gaussian, and (c) Euclidean.

Table 4 gives a comparison between the different PSNR values obtained after applying focusing delay, corrected attenuation coefficient, and different weights. We can infer that the images by applying linear weights are better in terms of PSNR.

Table 4. PSNR values in different scenarios

It can be seen in Figs 13(a) and 13(c) that even though the intensity and resolution of the targets are improved compared to raw images apart from the target, weak noise, and strong reflections due to the wall are also present. To remove the wall clutter and weak noise-like pattern novel hybrid technique is proposed in the next section.

Fig. 13. TWI images for the target ID 1: metal target with different weights: (a) linear, (b) Gaussian, and (c) Euclidean.

Proposed technique for wall parameter estimation and image de-noising

In this section, hybrid technique using the Opt-Shrink algorithm modified for the TWI application and data-driven Bayes shrink algorithm is proposed. Opt-Shrink algorithm gives optimal weights for the coefficients from the large class of noise models so Opt-Shrink algorithm is a good choice for de-noising of the TWI images. Implementing the steps for the conventional Opt-Shrink algorithm is described in the Section “TWI experimental setup and pre-processing” and the processed TWI images with the conventional algorithm are shown in Fig. 14.

Fig. 14. TWI images for (a) target ID 3: Teflon target and (b) target ID 1: metal target using conventional Opt-Shrink algorithm.

It can be seen by inspection that the noise is reduced compared to images in Fig. 13 still strong reflections from the wall are present as well as weak noise is present in the low-contrast areas of the images. In our earlier work [20] we mentioned that the first eigenvalue from low-rank approximation corresponds to the homogenous wall so Opt-Shrink algorithm for TWI can be modified as follows.

Steps to remove wall clutter and finding the optimum rank

The different steps for the implementation for Opt-shrink algorithm is given in the Section “Opt-shrink shrinkage and thresholding algorithm,” where we can modify steps 2 and 6 to remove wall clutter by ignoring the first eigenvalue corresponding to the wall i.e. $\hat{R} = \sum\nolimits_{i = 2}^q {{\hat{\sigma }}_i{\hat{u}}_i{\hat{v}}_i^T }$ and evaluate $\widehat{{S_{opt}}} = \sum\nolimits_{i = 2}^{\hat{r}} {\omega _{i, \hat{r}}^{opt} {\hat{u}}_i{\hat{v}}_i^T }$ where $\hat{r}$ is the effective rank, as mentioned in [14, 27] target sub-space is not spanned by 1D but by multi-dimensional sub-space and these sub-spaces split into signal and noise depending upon the factors like target electrical property, target location, target size, and the number of targets in the scene. In [15] it is recommended that multi-dimensional sub-space can be estimated from eyeballing the “knee” of the SV plot or from a random matrix theory-based test. This approach is not adequate for estimating the effective rank; hence in the next section we propose a novel approach using a GA to estimate the effective rank.

Estimation for effective optimum rank for low-rank signal matrix and improved wall parameter estimation

In this section, to estimate the effective rank for low-rank estimation we use the fourth moment statistical measure known as Kurtosis. The effective rank ($\hat{r}$) is to be determined at maximum PSNR value for different wall parameters such as thickness and permittivity. As discussed in the Section “Proposed wall characterization and focusing weights assignment,” these parameters are decided at minima or maxima for the estimation error.

A GA is a popular choice to find the local minima or maxima for any function. We will investigate the use of GA to find the effective rank at maximum PSNR value for the image by solving the multi-objective function in terms of thickness and permittivity. To find the optimum effective rank, evaluate $\widehat{{S_{opt}}} = \sum\nolimits_{i = 1}^{\hat{r}} {\omega _{i, \hat{r}}^{opt} {\hat{u}}_i{\hat{v}}_i^T }$ for the corresponding eigenvalue using SVD. To achieve this following fitness function is defined:

(26)$$[ {a, \;b} ] = \matrix{ {max} \cr {\hat{r}\;\ge 1} \cr } \mathop \sum \limits_{i = 1}^{\hat{r}} \omega _{i, \hat{r}}\hat{u}_i\hat{v}_i^T $$

(27)$$\eqalignno{& \hskip1pc{\rm i}.{\rm e}.\,\,F( {\hat{r}} ) = [ {\,f_1( {\hat{r}} ) } ] , \;\,1 \le ( {\hat{r}} ) \le length( {\omega_{i, \hat{r}}} ) \cr & \hskip1pc{\rm such\ that\ }f_1( {\hat{r}} ) = 1\quad {\rm i}.{\rm e}.\,{\rm maximum\ normalized\ PSNR\ for} \cr & \hskip1pc{\rm the\ target\ image}} $$

where $\omega _{i, \hat{r}}$ is the diagonal eigenvalue matrix using SVD, this eigenspace is split into signal and noise; hence effective rank estimation is necessary to improve the PSNR. GA is implemented for fixed value of b (either thickness or permittivity) and a is obtained within the bound i.e. for the wall thickness (13–15 cm) and the permittivity (4–6) as obtained by higher-order standardized moments in the Section “Proposed wall characterization and focusing weights assignment.” To verify the proposed method, we measured the wall parameters without placing any target and found that the simulated value for the wall permittivity is 5.3419 and thickness is 13.70 which are nearer to the actual values of the wall used in the experiment.

Table 5 shows the result obtained using GA for the estimation of low rank when permittivity is considered as a constant in the case of metal target.

Table 5. PSNR values for different values of rank ($\hat{r}$) at ɛ_r = 5.3

It can be observed in Table 5 that maximum PSNR is obtained at wall thickness d = 14.27 cm which is approximately the same as obtained in the Section “Wall effect minimization.” Furthermore, we verify permittivity value for the wall by considering the wall thickness as constant and found it also the same as shown in Table 6.

Table 6. PSNR values for different values of rank ($\hat{r}$) at d = 14 cm

At these parameter values, the effective rank can be determined and the same procedure can be repeated for Teflon and contrast targets. The estimated effective rank bound for the targets under investigation is shown in Fig. 15(a). It can be observed that it is in the range 3–5 and individual evaluated value is put for $\widehat{{S_{opt}}}{\rm \;}$ in step 6 of the Section “Steps to remove wall clutter and finding the optimum rank.” The plot for the different rank estimation methods for detecting the contrast targets in our experimental work at maximum PSNR is shown in Fig. 15(b). It is observed that PSNR is decreasing after the 10th eigenvalue in all cases. The comparison between the different methods for rank estimation available in the literature and our method is illustrated in Table 7.

Fig. 15 (a) Relationship between rank and wall characteristics. (b) Plot for estimated rank versus PSNR.

Table 7. Comparison between different rank estimation methods

Removal of the weak noise pattern and artifacts

After computing the $\widehat{{S_{opt}}}$ by using the above approach mentioned in the Sections “Steps to remove wall clutter and finding the optimum rank” and “Estimation for effective optimum rank for low-rank signal matrix and improved wall parameter estimation,” the next task is to remove the artifacts due to large MSE. After critically analyzing available filtering techniques in literature [28], we choose Daubechies wavelet filter bank which gives the orthogonal multi-resolution analysis in terms of low-pass and high-pass filter coefficients. The detail and approximate coefficients for the optimum eigenvalues then can find out using multilevel 2D wavelet decomposition. As mentioned in the Section “TWI experimental setup and pre-processing” wavelet thresholding keep low-resolution coefficient intact while filtering detail wavelet coefficients. We use wavelet Bayesian thresholding [19] modified as data-driven Bayesian thresholding for estimating the low-resolution coefficients. The data-driven parameters such as standard deviation and noise variance can be estimated as follows.

Estimating the data parameters

In this subsection, we will describe the method for effective thresholding of the image. If we consider the noise is i.i.d. Gaussian then observation model is represented as Y = X + V, with X and V are independent. To estimate the σ_x i.e. standard deviation, noise variance σ ² determined first from the detail wavelet coefficient of the high-pass filter at level 1 of the filter bank using robust median estimator [29, 30]:

(28)$${\mathop{\sigma}\limits^{\frown}} = \displaystyle{{median( Y_{ij}) } \over {0.6745}}Y_{ij}\in subband\ HH_1$$

(29)$$\sigma ^2_Y = \sigma ^2_X + \sigma ^2$$

$\sigma ^2_Y$ is the variance of Y having zero mean hence:

(30)$$\sigma ^2_Y = \displaystyle{1 \over {n^2}}\sum\limits_{i, j = 1}^n {Y_{ij}^2 } $$

n × n is the size of sub-band then

(31)$$Th( \sigma _X) = {{\sigma ^2} / {\sigma ^2_Y }}$$

where

(32)$$\sigma _X = \sqrt {\max ( \sigma ^2_Y -\sigma ^2, \;0) } $$

If σ ² > $\sigma ^2_Y$ then Th(σ _X) = 0, otherwise max(|Y_ij|).

The flow chart for the overall methodology along with the proposed hybrid technique is shown in Fig. 16. After implementing the procedure described for wall characterization and weighting in the Section “TWI experimental setup and pre-processing” and de-noising procedure in the Sections “Steps to remove wall clutter and finding the optimum rank,” “Estimation for effective optimum rank for low-rank signal matrix and improved wall parameter estimation,” and “Removal of the weak noise pattern and artifacts,” the final images obtained using the TWI system are shown in Fig. 17.

Fig. 16. Flow chart for proposed methodology including hybrid technique.

Fig. 17. Final TWI images: (a) target ID 3: Teflon target, (b) target ID 1: metal target, (c) target ID 7: metal and wood target, and (d) target ID 2: wood target using the proposed hybrid technique.

Conclusion

Through the wall imaging is becoming a popular technology in military and civilian applications. In real-time scenario characteristics of the wall are not known “a priori.” The characteristics of the wall used in the experimental work are estimated and verified with the thickness of the wall having characteristics d = 14 cm and ɛ_r = 5.3. The thickness of the wall obtained in the simulated result is d = 13.5 cm hence the error is less than 4% using higher-order standardized moments. The method we proposed to estimate different characteristics of the wall can be generalized for any type of wall. We observed an improvement in the image quality compared to a raw image when the estimated attenuation coefficient and linear weights are incorporated in the imaging process.

To remove the wall reflection from the image conventional, Opt-Shrink algorithm is modified for TWI application and wall clutter is removed successfully. To improve the PSNR of the developed microwave image noise space in the eigenvalue decomposition must be exploited effectively, this effective rank in Opt-Shrink algorithm is estimated using multi-objective GA. The hybrid technique based on Opt-Shrink and data-driven Bayesian thresholding is implemented which will be able to detect low-dielectric (low contrast) and high-dielectric (high contrast) targets successfully.