Performance Evaluation of Group Sparse Reconstruction and Total Variation Minimization for Target Imaging in Stratified Subsurface Media

Abstract

Sparse reconstruction methods have been successfully applied for efficient radar imaging of targets embedded in stratified dielectric subsurface media. Recently, a total variation minimization (TVM) based approach was shown to provide superior image reconstruction performance over standard L1-norm minimization-based method, especially in case of non-point-like targets. Alternatively, group sparse reconstruction (GSR) schemes can also be employed to account for embedded target extent. In this paper, we provide qualitative and quantitative performance evaluations of TVM and GSR schemes for efficient and reliable target imaging in stratified subsurface media. Using numerical electromagnetic data of targets buried in the ground, we demonstrate that GSR and TVM provide comparable reconstruction performance qualitatively, with GSR exhibiting a slight superiority over TVM quantitatively, albeit at the expense of less flexibility in regularization parameters.

1. Introduction

Effective and reliable imaging of targets embedded in stratified subsurface media is highly desirable in ground penetrating radar (GPR) applications [1,2,3,4,5,6,7,8,9,10,11,12,13]. Regularization methods based on $L^{p }\left(1<p<2\right)$ spaces framework have shown promise for microwave imaging of buried targets [14,15,16]. The sparse nature of targets has also been successfully exploited in GPR image recovery through sparse reconstruction approaches [4,17,18,19,20,21]. Most sparsity-based methods model the background environment in GPR as a two-layered medium, where radar operates in the upper air layer and the targets are buried within the lower ground layer. However, natural soil typically comprises multiple layers, which, if not properly addressed, can severely degrade the imaging performance.

Recently, a generalized sparse image reconstruction approach with total variation minimization (TVM) was proposed for efficient and reliable radar imaging through multilayered background media [4]. More specifically, the multilayered subsurface Green’s function was incorporated in the imaging algorithm to model the wave propagation effects in the multilayered environment and was efficiently evaluated using the saddle point method. As compared to standard $l_1$-norm minimization-based techniques [20,22], which are based on point target model, the TVM-based approach minimizes the gradient of the image, thereby leading to better edge preservation and, in turn, reconstruction of non-point-like and extended targets [4,21].

An alternative to TVM based approach is group sparse reconstruction (GSR), which can also account for the target extent [21,23]. In high-resolution images, each extended or non-point-like target generally occupies a contiguous group of pixels rather than a single pixel. As such, the point-target-based sparse signal model can be refined to exclude from the solution space any image whose support contains isolated pixel indices [22]. The group sparsity approach incorporates the clustering of the non-zero image pixels into a small number of contiguous groups as a constraint in the sparse reconstruction problem. The specific clustering pattern can be structured to match the desired target shape.

In this paper, we provide a performance evaluation of TVM and GSR schemes for non-point-like target imaging in stratified subsurface media. To this end, we consider a multiple-input multiple output (MIMO) radar system and use numerical electromagnetic data of targets buried in a four-layered background environment. We show that the two methods provide comparable performance qualitatively under noisy measurements. Quantitatively, GSR exhibits a slight superiority in terms of image reconstruction, especially at low signal-to-noise ratio (SNR) values. However, this GSR quantitative performance is achieved at the expense of less flexibility in setting of regularization parameters.

The remainder of the paper is organized as follows. Section 2 provides a review of the Green’s function formulation for modeling the wave propagation effects in the multilayered background media and details the TVM and GSR based imaging algorithms. Section 3 describes the considered metrics for quantitative assessment and provides performance comparison using two sets of image reconstruction results of targets embedded in stratified subsurface media. Concluding remarks are provided in Section 4.

2. Sparsity-Based Image Formation through Stratified Subsurface Media

In this section, we first describe the signal model for imaging through stratified subsurface media with a co-located MIMO radar. Then, we briefly review the TVM and GSR techniques for image reconstruction.

2.1. Signal Model in Matrix Form

We consider an N-element transmit array and an M-element receive array, with ${\mathit{r}}_{tn}=\left(x_{tn},z_{tn}\right)$ and ${\mathit{r}}_{rm}=\left(x_{rm},z_{rm}\right)$ denoting the respective position vectors of the $n$th transmitter and $m$th receiver. A stepped-frequency signal, with $P$ frequencies uniformly covering the frequency band $\left[f_{\min }, f_{\max }\right],$ is used for imaging. The transmitters are assumed to be activated sequentially, while simultaneous reception at all receivers is assumed. We focus on the four-layered background media, shown in Figure 1, where the first layer is air and the remaining three are subsurface layers. The dielectric constant and conductivity of the subsurface layers are assumed to be $\left({\epsilon }_{r2},{\sigma }_2\right), \left({\epsilon }_{r3},{\sigma }_3\right),$ and $\left({\epsilon }_{r4},{\sigma }_4\right),$ while the second and third layers have a thickness of $d_2$ and $d_3$, respectively. Although the presented formulation considers only three subsurface layers, it can be readily generalized to an arbitrary number of subsurface layers.

Assuming a point target embedded in the fourth layer, the received scattered field, $E_s\left({\mathit{r}}_{rm},{\mathit{r}}_{tn},k^p\right)$, at the mth receiver with the nth transmitter active can be expressed as

$$E_s\left({\mathit{r}}_{rm},{\mathit{r}}_{tn},k^p\right)=\int G\left({\mathit{r}}_{rm},\mathit{r},k^p\right)G\left(\mathit{r},{\mathit{r}}_{tn},k^p\right)\sigma \left(\mathit{r}\right)d\mathit{r},$$

(1)

where $\sigma \left(\mathit{r}\right)$ is the scene reflectivity at position $\mathit{r}=\left(x,z\right)$, $k^p$ is the free-space wave number of the pth frequency, and $G\left({\mathit{r}}_{rm},\mathit{r},k^p\right)$ and $G\left(\mathit{r},{\mathit{r}}_{tn},k^p\right)$ are the layered media Green’s functions characterizing wave propagation from the transmitter to the target and from the target to the receiver, respectively. We note that Equation (1) essentially uses the first-order Born approximation, which ignores the multiple scattering effects. We discretize the region being imaged in the xz-plane into $K\times L$ pixels, and represent the corresponding scene reflectivity by the $K\times L$ matrix $\tilde{s}.$ Then, the $m$th received signal in Equation (1) can be expressed in matrix form as

$${\tilde{\mathit{y}}}_{m,n}={\Psi }_{m,n}\mathit{s}, \mathit{s}=\mathrm{vec}\left(\tilde{\mathit{s}}\right),$$

(2)

where $\mathrm{vec}(·)$ returns the column-wise vectorization of its matrix argument, the pth element of ${\tilde{\mathit{y}}}_{m,n}$ is ${[{\tilde{\mathit{y}}}_{m,n}]}_p=E_s\left({\mathit{r}}_{rm},{\mathit{r}}_{tn},k^p\right),$ and ${\Psi }_{m,n}$ is a $P\times KL$ dictionary matrix encompassing the stratified media effects, with its (p, q)th element given by

$${[{\Psi }_{m,n}]}_{p,q}=G\left({\mathit{r}}_{rm},{\mathit{r}}_q,k^p\right)G\left({\mathit{r}}_q,{\mathit{r}}_{tn},k^p\right), p=1,\dots ,P, q=1,\dots ,KL.$$

(3)

The layered media Green’s function for the GPR imaging configuration in Figure 1 can be expressed in closed-form using the Saddle Point Method (SPM) as [4]

$$G\left({\mathit{r}}_R,\mathit{r},k^p\right)\approx \frac{j{k}_1^2}{4}F\left(k_1{\alpha }_0\right)e^{j\left(k_1\mathsf{\Phi }\left({\alpha }_0\right)+\pi /4\right)}\sqrt{\frac{2}{\pi k_1\left|{\mathsf{\Phi }}^{″}\left({\alpha }_0\right)\right|}},$$

(4)

where ${\mathit{r}}_R=\left(x_R,z_R\right)$, $k_1$ is the wavenumber of the air layer, α is a real-valued scaling variable, and

$$\mathsf{\Phi }\left(\alpha \right)=\left[\alpha \left|x_R-x\right|+d_2\sqrt{{\epsilon }_{r2}-{\alpha }^2}+d_3\sqrt{{\epsilon }_{r3}-{\alpha }^2}+z_R\sqrt{{\epsilon }_{r1}-{\alpha }^2}+\left(z-z_3\right)\sqrt{{\epsilon }_{r4}-{\alpha }^2}\right],$$

(5)

$$F\left(k_1\alpha \right)=\frac{\left(1+R_{12}\right)\left(1+R_{23}\right)\left(1+R_{34}\right)}{R_a{R}_b{k}_{1z}},R_a=1+R_{12}{\tilde{R}}_{23}\exp \left(j2k_{2z}{d}_2\right), R_b=1+R_{23}{R}_{34}\exp \left(j2k_{3z}{d}_3\right),$$

(6)

$${\tilde{R}}_{23}=\frac{R_{23}+R_{34}\exp \left(j2k_{3z}{d}_3\right)}{1+R_{23}{R}_{34}\exp \left(j2k_{3z}{d}_3\right)};R_{ij}=\frac{k_{iz}-k_{jz}}{k_{iz}+k_{jz}}, i=1,2,3, j=i+1,$$

(7)

with $z_3=-d_2-d_3$ and $k_{iz}=\sqrt{k_i^2-k_1^2{\alpha }^2}, i=1,2,3,4.$ Additionally, $k_i, i=2,3,4,$ is the complex wavenumber of the ith subsurface layer, and the stationary phase point ${\alpha }_0$ is determined from $\mathrm{Re}\left[{\mathsf{\Phi }}^{\prime }({\alpha }_0\right)]=0$. Stacking the measurements from all $N$ transmitters and $M$ receivers yields the signal model

$$\tilde{\mathit{y}}=\Psi \mathit{s},$$

(8)

where $\tilde{\mathit{y}}={[{\tilde{\mathit{y}}}_{1,1}^T,{\tilde{\mathit{y}}}_{1,2}^T,\dots ,{\tilde{\mathit{y}}}_{M,N}^T]}^T$, $\Psi ={[{\Psi }_{1,1}^T,{\Psi }_{1,2}^T,\dots ,{\Psi }_{M,N}^T]}^T,$ and ‘$T$’ denotes matrix transpose.

The signal model in Equation (8) corresponds to the full data measurements, comprising all $P$ frequencies from all $N$ transmitters and $M$ receivers. In many practical operational scenarios, there are often cost constraints, which may limit the number of transmitters and receivers available for deployment. Note that reducing the number of frequencies at which measurements are made over the desired bandwidth may not translate into cost reduction. This is because the antennas and radio frequency (RF) front end would still be required to operate over the entire frequency band. As such, we retain the use of all $P$ frequencies and assume that $N_t<N$ transmitters and $M_r<M$ receivers are available for data collection. Under these constraints, the model in Equation (8) takes the form

$$\mathit{y}=\Lambda \Psi s=\Theta \mathit{s}.$$

(9)

Here, $\mathit{\Lambda }$ is an $M_r{N}_{t}P\times NMP$ measurement matrix given by [24]

$$\Lambda =(\Phi ⨂I_{N_{t}P})⨂(\mathit{\vartheta }⨂I_{MP}),$$

(10)

where ‘$⨂$’ denotes the Kronecker product, $I_{(·)}$ is an identity matrix with the subscript indicating its dimensions, $\mathit{\Phi }$ is an $M_r\times M$ matrix constructed by randomly selecting $M_r$ rows of ${\mathit{I}}_{N_{t}P}$, and $\mathit{\vartheta }$ is an $N_t\times N$ matrix consisting of $N_t$ randomly selected rows of ${\mathit{I}}_{MP}$.

2.2. Total Variation Minimization

Using the reduced measurements in Equation (9), the unknown scene reflectivity vector $\mathit{s}$ can be recovered by solving the TVM problem [4],

$$\widehat{\mathit{s}}=\mathrm{argmin} {\Vert \mathit{s}\Vert }_{\mathrm{TV}}, \mathrm{subject}\mathrm{to}{\Vert \mathit{y}-\Theta \mathit{s}\Vert }_2<\delta ,$$

(11)

where $\delta$ represents a small tolerance error,

$$\begin{gathered} {\Vert s\Vert }_{\mathrm{TV}}={\displaystyle \sum }_{i=0}^{K-1}{\displaystyle \sum }_{j=0}^{L-1}\Vert {(\nabla \tilde{s})}_{i,j}\Vert {}_1={\displaystyle \sum }_{i=0}^{K-1}{\displaystyle \sum }_{j=0}^{L-1}{\Vert \left[\begin{array}{c}\tilde{s}\left(i+1,j\right)-\tilde{s}\left(i,j\right)\\ \tilde{s}\left(i,j+1\right)-\tilde{s}\left(i,j\right)\end{array}\right]\Vert }_1 \\ ={\displaystyle \sum }_{i=0}^{K-1}{\displaystyle \sum }_{j=0}^{L-1}{\Vert \left[\begin{array}{c}s\left(\left(j-1\right)K+i+1\right)-s\left(\left(j-1\right)K+i\right)\\ s\left(jK+i\right)-s\left(\left(j-1\right)K+i\right)\end{array}\right]\Vert }_1. \end{gathered}$$

(12)

In this work, we use the Nesterov algorithm in the NESTA package to solve the TVM problem in Equation (11) [25]. This algorithm utilizes a regularization scheme together with a smoothed version of $l_1$-norm to achieve the solution of the underlying convex optimization problem.

2.3. Group Sparse Reconstruction

In high-resolution imaging, targets generally occupy a group of neighboring pixels whose extent depends not only on the target dimension, but also on the system resolution. This prior pixel neighborhood information can be incorporated in the image reconstruction problem using group sparsity constraints. More specifically, the scene reflectivity vector $\mathit{s}$ can be obtained by solving the convex optimization problem [23,26,27]

$$\widehat{\mathit{s}}=\mathrm{argmin} \frac{1}{2}{\Vert \mathit{y}-\Theta \mathit{s}\Vert }_2^2+\lambda {\displaystyle \sum }_{q=0}^{KL-1}\Vert {\mathit{W}}^{\left(q\right)}{\mathit{s}}_{g_q}\Vert {}_2,$$

(13)

where $g_q\subseteq \left\{0,1,\dots ,KL-1\right\}$ is an index set corresponding to the group of pixels forming a neighborhood around the $q$th pixel and the diagonal weighting matrix ${\mathit{W}}^{\left(q\right)}$ ensures that the weighting within a group is according to the desired pixel neighborhood relation. Figure 2 provides an example of how the grouping of the image pixels works for a $10\times 10$ image [23]. The number in the top left corner of each square indicates the pixel index, whereas the pixel weight of the depicted group is represented by the number in the center of each pixel. The weights are chosen such that their sum equals unity to avoid unintentional scaling of the reconstruction result. As shown in Figure 2, the index set for the group corresponding to the 12th pixel is $g_{12}=\left\{2, 11, 12, 13, 22\right\},$ with the corresponding weighting matrix ${\mathit{W}}^{\left(12\right)}=\mathrm{diag}\left(\frac{1}{8},\frac{1}{8},\frac{1}{2},\frac{1}{8},\frac{1}{8}\right).$ In this paper, we solve the reconstruction problem in Equation (13) using Primal YALL1 group [27] and utilize overlapping groups.

3. Performance Evaluation Results and Discussion

In this section, we first describe the metrics considered for quantitative performance evaluation and then present the image reconstruction assessment results of targets embedded in a stratified subsurface.

3.1. Quantitative Metrics

We consider two different metrics, namely, relative clutter power (RCP) and Earth mover’s distance (EMD), for quantitative performance evaluation of the TVM and GSR schemes.

3.1.1. Relative Clutter Peak

We define the target region, ${\mathcal{R}}_t$, as the union of rectangular or circular regions at known target positions, whereas the remainder of the image comprises the clutter region, ${\mathcal{R}}_c$,. The size of each individual target region is determined based on the ground truth and the system resolution. With ${\mathcal{A}}_{{\mathcal{R}}_t,\widehat{s}}$ and ${\mathcal{A}}_{{\mathcal{R}}_c,\widehat{s}}$ denoting the respective maximum amplitude of target and clutter regions in the reconstructed image $\widehat{s},$ the RCP is defined as [23]

$$\mathrm{RCP}\mathrm{in}\mathrm{dB}=20{\log }_{10}\left(\frac{{\mathcal{A}}_{{\mathcal{R}}_t,\widehat{s}}}{{\mathcal{A}}_{{\mathcal{R}}_c,\widehat{s}}}\right) ,{\mathcal{A}}_{{\mathcal{R}}_t,\widehat{s}}=\underset{q\in {\mathcal{R}}_t}{\max }\left|{\widehat{s}}_q\right|,{\mathcal{A}}_{{\mathcal{R}}_c,\widehat{s}}=\underset{q\in {\mathcal{R}}_c}{\max }\left|{\widehat{s}}_q\right|,$$

(14)

where ${\widehat{s}}_q$ is the $q$th element of $\widehat{s}$. The RCP metric penalizes strong clutter and favors clean images with low noise and clutter power and high target amplitudes.

3.1.2. Earth Mover’s Distance

Earth mover’s distance is defined as the minimal amount of image intensity that has to be moved to transform one image into another [28,29]. For the underlying application, we measure the EMD between the reconstructed image and the ground truth image. This metric incorporates perceptual differences between the reconstructed and ground truth images, as it measures error in terms of not only the differences in pixel values, but also physical distance away from the actual target locations. It, therefore, is a preferred metric over mean-squared error in sparse reconstruction literature [30]. In this work, we use a fast implementation of EMD [31].

3.2. Performance Comparison

We consider a MIMO radar system with 17 uniformly spaced transmitters from −0.96 to 0.96 m and 16 receivers equally spaced from −0.9 to 0.9 m. Both arrays are at a height of 0.2 m above the ground. The stepped-frequency signal covers the 0.8 to 2 GHz bandwidth with $P=49$ frequency steps. A time-domain full wave electromagnetic solver based on Finite-Difference Time-Domain (FDTD) method is used for generating the received signals from two different scenes. Fast Fourier Transform (FFT) is applied to transform the time-domain received signals to frequency domain. The radar measurement configuration and signal parameters are the same in both scenarios. White Gaussian noise is added to the frequency-domain data. For image reconstruction, the thickness and complex permittivity of each layer of the stratified subsurface media are assumed to be known a priori. In practice, however, these parameters can be estimated using an inversion scheme. The inversion of multilayered medium parameters has been well developed within the framework of one-dimensional inverse scattering in the past two decades [32,33,34,35]. For the conventional single layer subsurface, analytical methods for estimation of the dielectric slab parameters have been provided in [32,35]. For multilayered subsurface media, the number and parameters of each layer can be efficiently retrieved using a layer-stripping algorithm or global optimization inverse scattering techniques [33,34].

3.2.1. Example 1

In this example, we consider three metallic targets (two rectangular and one cylindrical) embedded in a three-layered background environment, as shown in Figure 3. The dielectric constant, conductivity, and thickness of the second layer are ${\epsilon }_{r2}=6,{\sigma }_2=0.01$ S/m, and $d_2=0.2$ m, respectively. The third layer with a dielectric constant ${\epsilon }_{r3}=3$ and conductivity ${\sigma }_3=0.005$S/m contains the three targets. The target dimensions are specified in

Table 1. Target characteristics for Example 1.

Type	Dimensions	Center Position
Rectangular	0.2 m $\times$ 0.1 m	(−0.6 m, −0.75 m)
Cylindrical	Radius = 0.1 m	(0, −0.8 m)
Rectangular	0.2 m $\times$ 0.1 m	(0.6 m, −0.75 m)

. We randomly select two transmitters (12% of the available quantity) and six receivers (38% of the available number). For each chosen transmitter-receiver pair, we utilize all 49 frequency measurements to reconstruct the image. Figure 4 and Figure 5 depict the images obtained using TVM, GSR, and standard $l_1$-norm sparse reconstruction [20] for −10 and −5 dB SNR values, respectively. We observe from Figure 4 and Figure 5 that although there are a few false reconstructions in both TVM and GSR results, the two approaches provide cleaner images of the targets as compared to the standard $l_1$-norm sparse reconstruction result.

Next, we quantitatively evaluate the performance of the TVR, GSR, and standard $l_1$-norm reconstruction schemes. We consider SNR values in the [−10 10] dB range with 5 dB increments. We perform a total of 100 Monte Carlo trials for each SNR value with different realization of noise and different randomly chosen sets of two transmitters and six receivers each time. For every trial, we reconstruct the image using TVM, GSR, and standard $l_1$-norm schemes and compute the corresponding values of the metrics. Figure 6 plots the RCP and the EMD, each averaged over 100 trials, versus SNR. The variance of the EMD for each SNR is also indicated in Figure 5. We observe that GSR and TVM provide almost identical RCP performance, while significantly outperforming standard $l_1$-norm reconstruction for all SNR values. In terms of EMD at low SNR, GSR provides the best performance as manifested by the lowest EMD average and variance values, while standard $l_1$-norm reconstruction has the worst performance with TVM in the middle. At higher SNR values, the average EMD values for all three reconstruction methods are comparable. However, both GSR and TVM yield a smaller variance for the EMD as compared to the standard $l_1$-norm reconstruction. Similar trends were observed for reconstruction performance with a random selection of four transmitters and four receivers and a random selection of four transmitters and eight receivers.

3.2.2. Example 2

In this example, we consider a large composite metallic target consisting of two semi-cylinders on top of a rectangle cylinder, which is embedded in the third layer of a three-layer background as shown in Figure 7. The target dimensions are also specified in Figure 7. The physical and electrical properties of the three background layers are the same as in Example 1. Figure 8 and Figure 9 depict the reconstruction results obtained with two randomly chosen transmitters and six randomly chosen receivers using TVM, GSR, and standard $l_1$-norm sparse reconstruction for SNR values of −5 and 0 dB, respectively. Similar to Example 1, we observe that both TVM and GSR approaches yield superior quality images as compared to the standard $l_1$-norm sparse reconstruction.

The two quantitative performance metrics, averaged over 100 Monte Carlo trials, are plotted vs. SNR in Figure 10 for TVR, GSR, and standard $l_1$-norm reconstruction schemes. For EMD, the variance is also indicated in Figure 10. Again, we observe that the performance of GSR and TVM is comparable in terms of RCP, while that of standard $l_1$-norm reconstruction is significantly lower. Unlike the three-target scene in Example 1 wherein TVM performance in terms of average EMD was approximately half-way between that of GSR and standard $l_1$-norm reconstruction, the EMD curve for TVM in the larger composite target case closely follows the corresponding curve for GSR. TVM’s capability of edge preservation better manifests itself in case of the composite target, leading to smaller performance difference with GSR as compared to the three-targets scene in Example 1.

Reconstructions with a random selection of four transmitters and four receivers, and a set of four transmitters and eight receivers, both sets chosen at random, yielded similar quantitative performance trends.

3.3. Discussion

The results provided in Section 3.2 quantify and validate the superior performance of the GSR and TVM approaches over the standard $l_1$-norm reconstruction for non-point-like targets, especially at low SNR values. A comment is in order on the choice of regularization/penalty parameters for the employed sparse reconstruction methods. The Primal YALL1 group solver requires setting of a length-2 penalty parameter vector [27], whereas NESTA, employed for TVM and standard $l_1$-norm reconstructions, requires the specification of smoothing and stopping parameters. Both the smoothing and stopping parameters in NESTA should be set to small values for higher accuracy or large values for faster convergence [36]. In general, choosing a small value for the smoothing parameter warrants a small value of the stopping parameter. For large smoothing parameter value, the stopping parameter can also be larger. Note that setting the smoothing parameter equal to zero results in use of the standard “non-smoothed” version of the $l_1$-norm in the optimization problem. For the Primal YALL1 group solver, the elements of the length-2 penalty vector are set as inversely proportional to ${\Vert {\Theta }^H\mathit{y}\Vert }_{\infty },$ where ${\Vert ·\Vert }_{\infty }$ denotes the infinity norm of the argument [23].

For the considered numerical experiments, we set these parameters for NESTA and YALL1 empirically for a nominal number of transmitters and receivers under noise-free conditions following the aforementioned guidelines, opting for higher accuracy over faster convergence for NESTA. For NESTA, no adjustments were made when the total number of transmitters and receivers were increased or decreased compared to the nominal case. However, the proportionality constants for the penalty parameters in YALL1 group had to be adjusted to account for the change in ${\Vert {\Theta }^H\mathit{y}\Vert }_{\infty }$ with any increase or decrease in the number of transmitters and receivers employed. Thus, compared to YALL1 group, NESTA was found to be more robust to changes in the amount of data employed for the sparse reconstructions.

4. Conclusions

In this paper, we conducted qualitative and quantitative performance evaluations of group sparse reconstruction and total variation minimization approaches for radar imaging through stratified subsurface media. The TVM approach minimizes the gradient of the image resulting in good edge preservation, while group sparse approach exploits prior pixel neighborhood information about extended targets for reliable imaging. Under reduced number of transmitters and receivers, numerical EM measurements with varying SNR levels were considered. Both TVM and GSR approaches demonstrated comparable qualitative performance for different subsurface scenarios. The quantitative evaluation revealed a slight performance advantage of GSR over TVM. However, this advantage came at the cost of reduced flexibility in setting regularization/penalty parameters for GSR vs. TVM.