Extraction of Anisotropic Characteristics of Scattering Centers and Feature Enhancement in Wide-Angle SAR Imagery Based on the Iterative Re-Weighted Tikhonov Regularization

Abstract

The anisotropic characteristic reflects discriminating features of the geometry of a scattering center. In this study, we propose a novel method based on the iterative re-weighted Tikhonov regularization (IRWTR) to achieve the extraction of anisotropic characteristics of scattering centers in the wide-angle SAR synthetic aperture radar (SAR) imaging. Moreover, based on the extracted anisotropic scattering behaviors, the incomplete imaging results of the distributed scattering centers are restored. In this paper, we first discussed the scattering property in SAR imagery from the perspective of attributed scattering center model (ASCM). The reason for the incomplete imaging results of the distributed scattering centers was also discussed based on the ASCM. Subsequently, we modeled the aspect-dependent amplitude responses of a scattering center as a linear combination of a set of orthogonal basis. Consequently, the extraction of anisotropic characteristics can be transformed into an inverse problem, which was solved by the proposed IRWTR with high efficiency and accuracy. After the extraction, we attempted to restore the complete image of a distributed scattering center, which consisted of only two points. The enhanced SAR image can reveal the actual shape of a target. Processing results of electromagnetic computation data validated that the proposed method is effective and efficient.

1. Introduction

Wide-angle SAR synthetic aperture radar (SAR) imagery can provide SAR images with higher azimuth resolution in comparison with the traditional small-angle SAR. In addition to the finer cross-range resolution, wide-angle SAR also benefits the detection and extraction of anisotropic characteristics of scatterers. Anisotropic characteristics, which can be identified in wide-angle SAR imaging, provide us with more interesting and important information about targets. Anisotropic characteristic implies the amplitude response of a scattering center can vary with the look angle [1,2,3]. We can infer geometric features of a scatterer from its anisotropic characteristic. Therefore, it is well worth studying the method to extract anisotropy characteristics from radar measurements to broaden the understanding of targets in wide-angle SAR image.

Because of its importance, some publications have addressed the methods of extracting anisotropic characteristics. In [4,5,6,7], the extraction of anisotropic characteristics was achieved via the parameter estimation of the parametric scattering models. In these publications, the attributed scattering center model (ASCM) [8,9,10] was used to represent the scattering behavior. The parameters of the ASCM that can make the model fit the radar measurements were searched for. Then, based on the parameters and the ASCM, the anisotropic characteristics could be retrieved from the radar measurements. In [8,9], the ASCM was introduced and used for automatic target recognition (ATR). In [4], the authors used the modified RELAX algorithm to extract ASCM parameters of scattering centers. The high-frequency multipeak model was derived from the ASCM and was used to describe the scattering behavior in the SAR imagery in [5]. This model was used for joint frequency–polarization scatter classification. In [6], the algorithm based on the sparse representation with dictionary-refinement was proposed to extract attributed scattering centers and estimate ASCM parameters. A statistical signal processing algorithm was proposed to extract ASCs attributed scattering centers (ASCs) from interferometric circular SAR measurements in [7]. Except for the methods based on the parametric models, [2] proposed an effective method to directly extract anisotropic characteristics based on the azimuth-dependent scattering feature representation. In [2], an over-complete dictionary containing the possible anisotropic scattering behaviors was adopted and an efficient method to search for the elements in the dictionary matching the scatterers’ scattering behaviors was also designed.

Following the idea in [2], we adopted a set of the orthogonal basis to represent the amplitude responses of a scattering center and model the radar measurements as the sum of responses from individual scattering centers. Then, the inversion of the amplitude responses of scattering centers at different aspect angles from the radar measurements could be achieved by solving the inverse problem. To solve this inverse problem with high accuracy, we added the Tikhonov regularization with multiple regularization parameters (MRP) to ‘regularize’ the solution. The regularization parameters are relevant to the anisotropic characteristics of scattering centers. As the scattering attributes of the scatterers cannot be known in advance, we used the Expectation Maximization (EM) algorithm to solve the MRP Tikhonov regularization problem, which resulted in an iteratively re-weighting algorithm. Therefore, we refer to the proposed method as iterative re-weighted Tikhonov regularization (IRWTR). The IRWTR can select proper regularization parameter for each scattering center based on its own scattering property adaptively. By using this method, an accurate estimation of the anisotropic characteristics can be achieved with high efficiency.

After the extraction of anisotropic characteristics, we tried to restore the incomplete image of the distributed scattering center, which consisted of only two points based on the extracted anisotropic characteristics. According to the ASCM and the phenomenon in the practice, we found that the imaging result of a distributed scattering center in a wide-angle SAR image may have only two separated points, which can impair the SAR-image understanding. With the anisotropic characteristics, we found that some points in the SAR image were actually the imaging results of the ends of distributed scattering centers. Furthermore, the defective images can be restored based on the anisotropic characteristics.

The content of this paper has been organized as follows: in Section 2, the signal model is presented and we analyze the reason why the imaging result of a distributed scattering center may only contain two isolated points; the proposed algorithm is introduced in detail in Section 3; the results obtained via the proposed method are displayed in Section 4. The conclusion is drawn in Section 5.

2. Analysis of the Signal Model in Wide-Angle SAR Imaging

2.1. Problem Formulation

The simplified two-dimensional (2D) wide-angle SAR imaging geometric model is displayed in Figure 1 based on the contents in [1,2,3]. During the acquisition time, the radar platform moves along the flight path in Figure 1. $O$ is the center of the antenna footprint and $p$ is a scattering center in the imaging scene. The radar transmitted signals, which illuminated the imaging scene. The signals reflected from the imaging scene were collected by the radar, which were determined by the scattering characteristics and the positions of the targets in the imaging scene. The radar transmitted and collected the signals at different positions on the flight path. The radar antenna always illuminated the imaging scene during the procedure. With the assumption of far-field and plane-wave propagation, $\phi$ denotes the look angle of the radar that is identical for all the scattering centers in the imaging scene.

In wide-angle SAR imaging, the dependence of the amplitude response on the look angle $\phi$ should be taken into consideration. In wide-angle SAR imaging, the amplitude response of the $p$th scattering center can be expressed as $E_p^{}\left(f,\phi ,x_p^{},y_p^{}\right)$, where $f$ is the frequency of the transmitted signal. $\left(x_p^{},y_p^{}\right)$ denotes the position of the scattering center. From the perspective of ASCM, $E_p^{}\left(f,\phi ,x_p^{},y_p^{}\right)$ can be expressed as:

$$E_p^{}\left(f,\phi ,x_p^{},y_p^{}\right)=A_p^{}{\left(j\frac{f}{f_c^{}}\right)}_{}^{{\alpha }_p}\cdot \exp \left(-2\pi {\gamma }_{p}f\sin \phi \right)\exp \left(\left(-j\frac{4\pi f}{c}\right)\cdot \left(x_p\cos \phi +y_p\sin \phi \right)\right)$$

(1)

for a local scattering center or:

$$\begin{array}{c}{E}_p\left(f,\phi ,x_p,y_p\right)=A_p{\left(j\frac{f}{f_c}\right)}^{{\alpha }_p}\mathrm{sinc}\left(\left(\frac{2\pi f{L}_p}{c}\right)\cdot \sin \left(\phi -{\overline{\phi }}_{0p}\right)\right)\\ \cdot \exp \left(\left(-j\frac{4\pi f}{c}\right)\cdot \left(x_p\cos \phi +y_p\sin \phi \right)\right)\end{array}$$

(2)

for a distributed scattering center. Trihedral, top hat, and sphere are canonical local scattering centers. Distributed scattering centers include dihedral, flat plate, cylinder, and broadside edge diffraction. Generally, the local scattering center shows isotropic characteristic. In contrast, distributed scattering center presents anisotropic characteristic. In (1) and (2), $\sin \mathrm{c}\left(x\right)=\frac{\sin x}{x}$ and $j=\sqrt{-1}$. $A_p$ is the relative amplitude of the $p$th scatterer; ${\alpha }_p$ reflects the frequency dependence; $\left\{{\gamma }_p,L_p,{\overline{\phi }}_{0p}\right\}$ reflects the aspect dependence. ${\gamma }_p$ denotes the tiny aspect dependence of a local scattering center, whose influence can be omitted in practical. $\left\{L_p,{\overline{\phi }}_{0p}\right\}$ are the length and the orientation angle of a distributed scattering center. The orientation angle is the angle between the flank of a distributed scattering center and X-axis of the imaging scene.

The SAR image of a local scattering center is a point. The images of a distributed scattering center with different orientations are quite different. For a distributed scattering center, if $\phi ={\overline{\phi }}_{0p}$ at a look angle, the image will be a line whose length equals the length of the distributed scattering center. However, if $\phi \ne {\overline{\phi }}_{0p}$ in the acquiring procedure, the imaging result of the distributed scattering center will be composed of two isolated points. If $\phi \ne {\overline{\phi }}_{0p}$ in the acquiring procedure, $\mathrm{sinc}\left(\left(\frac{2\pi f{L}_p}{c}\right)\cdot \sin \left(\phi -{\overline{\phi }}_{0p}\right)\right)$ in (1) can be written as:

$$\begin{array}{cc}\mathrm{sinc}\left[\left(2\pi f{L}_p/c\right)\cdot \sin \left(\phi -{\overline{\phi }}_{0p}\right)\right]& =\frac{\sin \left(\left(2\pi f{L}_p/c\right)\cdot \sin \left(\phi -{\overline{\phi }}_{0p}\right)\right)}{\left(2\pi f{L}_p/c\right)\cdot \sin \left(\phi -{\overline{\phi }}_{0p}\right)}\hfill \\ & =\frac{e^{j\left(2\pi f{L}_p/c\right)\cdot \sin \left(\phi -{\overline{\phi }}_{0p}\right)}-e^{-j\left(2\pi f{L}_p/c\right)\cdot \sin \left(\phi -{\overline{\phi }}_{0p}\right)}}{2j\left(2\pi f{L}_p/c\right)\cdot \sin \left(\phi -{\overline{\phi }}_{0p}\right)}\hfill \end{array}$$

(3)

From (3), a distributed scattering center can be regarded as two isolated local scattering centers and their responses can be expressed as:

$$E_{p1}\left(f,\phi \right)=\frac{A_p{\left(j\cdot f/f_c\right)}^{{\alpha }_p}}{2j\left(2\pi f{L}_p/c\right)\cdot \sin \left(\phi -{\overline{\phi }}_{0p}\right)}\exp \left(\frac{-j4\pi f}{c}\cdot \left(x_1\cos \phi +y_1\sin \phi \right)\right)$$

(4)

and:

$$E_{p2}\left(f,\phi \right)=\frac{-A_p{\left(j\cdot f/f_c\right)}^{{\alpha }_p}}{2j\left(2\pi f{L}_p/c\right)\cdot \sin \left(\phi -{\overline{\phi }}_{0p}\right)}\exp \left(\frac{-j4\pi f}{c}\cdot \left(x_2\cos \phi +y_2\sin \phi \right)\right)$$

(5)

According to (4) and (5), a distributed scattering center becomes two local scattering centers, and they are at the ends of the distributed scattering center. Consequently, the imaging result of the distributed scattering center only comprises of two points. To illustrate this phenomenon, the backscattered data of the distributed scattering centers with different orientation angles are simulated by FEKO [11] and the corresponding images are shown in Figure 2 and Figure 3. The relative information of this simulation is given in Section 4.2.

The blue lines in the Figure 2a and Figure 3a are the look angles. The straight line connecting a point on the blue line and the origin of the coordinate defines the line of sight. In Figure 2a, during the acquisition procedure, the incident direction is vertical to the flank of the distributed scattering center at a look angle. The corresponding imaging result shown in Figure 2b is complete. In comparison, in Figure 3b, there are only two points. From Figure 3b, we cannot know the actual shape of the distributed scattering center. Furthermore, the distributed scattering center with an incomplete image can be easily mistaken for two local scattering centers.

2.2. Signal Model of Return Signals

The frequency-dependent factor $g\left(f\right)={\left(j\cdot f/f_c\right)}^{{\alpha }_p}$ and aspect-dependent factor $s_p\left(f,\phi \right)=\mathrm{sinc}\left[\left(2\pi f{L}_p/c\right)\cdot \sin \left(\phi -{\overline{\phi }}_{0p}\right)\right]$ (for the distributed scattering centers) or $s_p\left(f,\phi \right)=\exp \left(-2\pi {\gamma }_{p}f\sin \phi \right)$ (for the local scattering centers) can be simplified as $g\left(f\right)\approx 1$ and $s_p\left(f,\phi \right)\approx s_p\left(f_c,\phi \right)$, where $f_c^{}$ is the center frequency. The approximations are based on the fact that the bandwidth of a normal radar system is quite small in comparison with the center frequency [1,4]. Then, the simplified ASCM can be expressed as:

$$E_p\left(f,\phi \right)=A_p{s}_p\left(f_c,\phi \right)\exp \left(\left(-j\frac{4\pi f}{c}\right)\cdot \left(x_p\cos \phi +y_p\sin \phi \right)\right)$$

(6)

The extraction of anisotropic characteristics can be achieved by estimating $s_p\left(f_c,\phi \right)$.

Assuming there are $P$ scattering centers in the scene and the radar measurements are collected at $K$ discrete frequencies and $N$ discrete look angles. The received signal $y\left(f_k,{\phi }_n\right)$, which is collected by the radar at the frequency $f_k$ and the look angle ${\phi }_n$, can be expressed as:

$$\begin{array}{cc}y\left(f_k,{\phi }_n\right)& ={\displaystyle \sum _{p=1}^P{E}_p\left(f_k,{\phi }_n,x_p,y_p\right)}+n\left(f_k,{\phi }_n\right)\hfill \\ & ={\displaystyle \sum _{p=1}^P{A}_p\cdot s_p\left(f_c,{\phi }_n\right)\cdot \exp \left(\left(-j4\pi f_k/c\right)\cdot \left(x_p\cos {\phi }_n+y_p\sin {\phi }_n\right)\right)}+n\left(f_k,{\phi }_n\right)\hfill \end{array}$$

(7)

where $n\left(f_k,{\phi }_n\right)$ is the zero mean additive complex white Gaussian noise and $A_1^{},A_2^{},\cdots ,A_P^{}$ are the relative amplitudes of the $P$ scattering centers, respectively. We used a vector set $B=\left[b_1,b_2,\cdots ,b_m,\cdots ,b_M\right]$ including $M$ column vectors to represent $s_p={\left[s_p\left(f_c,{\phi }_1\right),\cdots ,s_p\left(f_c,{\phi }_n\right),\cdots s_p\left(f_c,{\phi }_N\right)\right]}^{\mathrm{T}}$. $b_m$ is a column vector containing N elements. Then, $s_p$ can be represented as:

$$s_p={\displaystyle \sum _{m=1}^M{a}_{p,m}{b}_m}$$

(8)

The matrix $B$ plays an important role in designing the processing algorithm of anisotropic characteristic extraction. We will discuss $B$ in Section 3. Subsequently, we can use (8) to express (7) in matrix form:

$$y=\Phi a+n$$

(9)

The variables in (9) are listed below:

$$\begin{array}{l}y={\left[y_1,y_2,\cdots y_K\right]}^T\\ y_k=\left[y\left(f_k,{\phi }_1\right),y\left(f_k,{\phi }_2\right),\cdots ,y\left(f_k,{\phi }_N\right)\right]\\ \Phi =\left[\begin{array}{cccc}{A}_1·\left(B\odot \left(r_{1,1}{1}_M^T\right)\right)& A_2·\left(B\odot \left(r_{2,1}{1}_M^T\right)\right)& \cdots & A_P·\left(B\odot \left(r_{P,1}{1}_M^T\right)\right)\\ A_1·\left(B\odot \left(r_{1,2}{1}_M^T\right)\right)& A_2·\left(B\odot \left(r_{2,2}{1}_M^T\right)\right)& \cdots & A_P·\left(B\odot \left(r_{P,2}{1}_M^T\right)\right)\\ ⋮& ⋮& \ddots & ⋮\\ A_1·\left(B\odot \left(r_{1,K}{1}_M^T\right)\right)& A_2·\left(B\odot \left(r_{2,K}{1}_M^T\right)\right)& \cdots & A_P·\left(B\odot \left(r_{P,K}{1}_M^T\right)\right)\end{array}\right]\\ a={\left[a_1,a_2,\cdots a_P\right]}^T\\ a_p=\left[a_{p,1},a_{p,2},\cdots a_{p,M}\right]\\ r_{p,k}={\left[e^{-j\frac{4\pi f_k}{c}\left(x_p\cos {\phi }_1+y_p\sin {\phi }_1\right)},e^{-j\frac{4\pi f_k}{c}\left(x_p\cos {\phi }_2+y_p\sin {\phi }_2\right)},\cdots ,e^{-j\frac{4\pi f_k}{c}\left(x_p\cos {\phi }_N+y_p\sin {\phi }_N\right)}\right]}^T\\ n={\left[n_1,n_2,\cdots ,n_K\right]}^T\\ n_k=\left[n\left(f_k,{\phi }_1\right),n\left(f_k,{\phi }_2\right),\cdots n\left(f_k,{\phi }_N\right)\right]\end{array}$$

where $\odot$ represents Hadamard product and $1_M$ is a $M\times 1$ vector of ones. According to (9), we can get $a$ by solving an inverse problem. After obtaining $a$, the azimuth-dependent amplitudes of scattering centers can be retrieved by $B$ and (8).

3. The IRWTR Algorithm and Feature Enhancement of Wide-Angle SAR Images

In this section, we firstly present the IRWTR algorithm for the extraction of anisotropic characteristics. Secondly, based on the anisotropic characteristics extracted by the IRWTR and the minimization of the least squares error criterion, we propose a method to restore the incomplete imaging results of certain distributed scattering centers.

A diagram to illustrate the proposed algorithm is shown in Figure 4. Figure 4a shows the imaging scene, which includes a metal rod and metal ball. The corresponding imaging result is shown in Figure 4b. There are three points in the image, which cannot reflect the actual imaging scene. We used the IRWTR algorithm to extract the aspect-dependent amplitudes of the scatterers that are shown in Figure 4c. From the aspect-dependent amplitudes, we restored the incomplete imaging result of the metal rod, which is shown in Figure 4d.

In the following subsections, the main steps of the proposed method are introduced.

3.1. Selection of Scattering Centers and Estimation of the Relative Amplitudes and Positions

In the first part of the proposed scheme, the positions and relative amplitudes of scatterers were estimated. Firstly, we used the Backprojection (BP) algorithm [12] to achieve the SAR image. Based on the imaging result, we could perform the estimation of the positions of scatterers. After the position estimation, the RELAX method [4] was utilized to accurately estimate the relative amplitudes of the scatterers. The positions and relative amplitudes were required in forming $\Phi$.

3.2. Anisotropic Characteristic Estimation

Apart from the positions and the relative amplitudes of the scattering centers, we still needed an appropriate matrix $B$ to form matrix $\Phi$. In this paper, we proposed a parsimony $B$ which is an identity matrix with the dimensions $N\times N$ and can be expressed as:

$$B=\left[\begin{array}{cccc}1& 0& \cdots & 0\\ 0& 1& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & 1\end{array}\right]$$

(10)

The physical meaning of $B$ is clear and sample. Its $n$th column is the basis that only represents the amplitude of a scattering center at the $n$th look angle. Then, the estimation of the scattering centers’ aspect-dependent amplitudes can be achieved by solving $a$ from $y$, which can be expressed as the following problem:

$$\widehat{a}=\underset{a}{\arg \min }{\Vert y-\Phi a\Vert }_2^2$$

(11)

Problem (11) is a least squares problem, and its analytic solution is:

$$\widehat{a}=\left({\left({\Phi }^H\Phi \right)}^{-1}\Phi \right)y$$

(12)

Equation (12) only provides a simple and reasonable solution. Under the condition of noise and other perturbations, the performance of (12) may be poor.

Because of the poor performance of LS method, it was necessary to impose a regularization term on the problem in (11), which could make the solution more accurate and suppress noise’s influence. We used the Tikhonov regularization term to make the solution more accurate. The rationality of using Tikhonov regularization was based on the fact that the aspect-dependent amplitude function of look angle of a scatterer was normally continuous and smooth. The objective function with the Tikhonov regularization can be expressed as:

$$\widehat{a}=\underset{a}{\arg \min }{\Vert y-\Phi a\Vert }_2^2+\lambda {\Vert Fa\Vert }_2^2$$

(13)

where $F$ is a gradient operator and can be expressed as:

$$F=\left[\begin{array}{cccc}D& & & \\ & D& & \\ & & \ddots & \\ & & & D\end{array}\right]$$

(14)

$D$ is the matrix of first-order derivatives with the dimensions $\left(M-1\right)\times M$ and can be expressed as:

$$D=\left[\begin{array}{ccccc}-1& 1& 0& \cdots & 0\\ 0& -1& 1& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& -1& 1\end{array}\right]$$

(15)

$\lambda$ is the regularization parameter to control the tradeoff between ${\Vert y-\Phi a\Vert }_2^2$ and ${\Vert Fa\Vert }_2^2$. With a relative smaller regularization parameter, the solution will vary rapidly. Consequently, a smaller regularization parameter was suitable to retrieve the anisotropic scattering. In contrast, a larger regularization parameter was appropriate for the inversion of the isotropic scattering. To use a uniform regularization parameter was not appropriate. It was necessary to choose specific $\lambda$ for each scattering center based on its own scattering property. By taking the individual scattering attributes into consideration, we proposed the MRP Tikhonov regularization. The problem with the MRP Tikhonov regularization can be expressed as:

$$\widehat{a}=\arg \min \left[{\Vert y-\Phi a\Vert }_2^2+{\Vert F_{w}a\Vert }_2^2\right]$$

(16)

where $F_w$ is expressed as:

$$F_w=\left[\begin{array}{cccc}{\lambda }_1^{}D& & & \\ & {\lambda }_2^{}D& & \\ & & \ddots & \\ & & & {\lambda }_P^{}D\end{array}\right]$$

(17)

In (17), ${\lambda }_p^{}$ is the regularization parameter for the $p$th scattering center. Since $\lambda ={\left[{\lambda }_1^{},{\lambda }_2^{},\cdots ,{\lambda }_P^{}\right]}^T$ has a significant influence on the solution, the selection of it is crucial. In order to accurately select a suitable $\lambda$, we discussed the exact physical meaning of $\lambda$, based on which selection method was introduced.

Similar to [13], we also revealed the regularization parameter’s meaning from the perspective of the maximum a posteriori (MAP) estimation. If we use the MAP theory to understand (16), the procedure of estimating $a$ is:

$$\begin{array}{cc}\widehat{a}& =\arg \max P\left(a|y\right)\hfill \\ & =\arg \max \left[P\left(y-\Phi a\right)·P\left(a\right)\right]\hfill \end{array}$$

(18)

where $P\left(y-\Phi a\right)$ can be expressed as:

$$\begin{array}{cc}P\left(y-\Phi a\right)& =P\left(n\right)\hfill \\ & =\frac{1}{{\left(\pi {\sigma }_n^2\right)}^{NK}}\exp \left(-\frac{{\Vert y-\Phi a\Vert }_2^2}{{\sigma }_n^2}\right)\hfill \end{array}$$

(19)

where $P\left(n\right)$ is the probability density function (PDF) of the noise which is a complex Gaussian PDF and ${\sigma }_n^2$ is the noise variance. ${\sigma }_n^2$ can be estimated from the SAR imaging result [13]. To estimate ${\sigma }_n^2$ from the data, we selected an area of the imaging result of the data which contained no targets. Then, the estimation of ${\sigma }_n^2$ can achieved by calculating the mean power of this area. If $I$th pixels of the image are used, the mean square error of the estimator is $\frac{{\sigma }_n^2}{I}$. P(a) is the PDF of the amplitudes of scatterers and can be expressed as:

$$P\left(a\right)=P\left(a_1^{}\right)P\left(a_2^{}\right)P\left(a_3^{}\right)\cdots P\left(a_p^{}\right)\cdots P\left(a_P^{}\right)$$

(20)

since the distributions of the amplitudes of different scattering centers were independent. In $a_p^{}$, $a_{p,n}$ is the amplitude response of the $p$th scattering center at the $n$th look angle, which indicates that adjacent elements in $a_p^{}$ show correlation. Consequently, it is reasonable to model $a_p^{}$ as the first-order Gaussian-Markov process [14,15]:

$$a_{p,n}=a_{p,n-1}+w_{p,n}$$

(21)

where $w_{p,n}$ is a random variable and $w_{p,n}\sim CN\left(0,{\sigma }_p^2\right)$. For a different $n$, $w_{p,n}$ is independent and identically distributed (i.i.d). The real part and imaginary part of $a_{p,1}$ are two random variables with the independent uniform distributions whose PDF are a constant $C_p$. With a larger ${\sigma }_p^2$, the variation between two adjacent elements of $a_p^{}$ would be more bigger, which indicates that the anisotropic feature of $p$th scattering center is more significant. In Figure 5, we show an example to illustrate the influence of ${\sigma }_p^2$. In Figure 5, it is clear that the elements in $a_p^{}$ varies more significantly as ${\sigma }_p^2$ increases.

Based on (21) and the PDF of $a_{p,1}$, $P\left(a_p\right)$ could be written as:

$$\begin{array}{cc}P\left(a_p\right)& =P\left(a_{p,1},a_{p,2},\cdots a_{p,N}\right)\hfill \\ & =P\left(a_{p,1}\right)P\left(w_{p,2}\right)P\left(w_{p,3}\right)\cdots P\left(w_{p,N}\right)\hfill \\ & =C_p^2{\displaystyle \prod _{n=2}^N\frac{1}{\pi {\sigma }_p^2}\exp \left(-\frac{{\left(w_{p,n}\right)}^2}{{\sigma }_p^2}\right)}\hfill \end{array}$$

(22)

According to (18), (19) and (22), one has:

$$\begin{array}{cc}\widehat{a}& =\arg \max \left[P\left(y-\Phi a\right)·P\left(a\right)\right]\hfill \\ & =\arg \max \left[\ln P\left(y-\Phi a\right)+\ln P\left(a\right)\right]\hfill \\ & =\arg \min \left[{\Vert y-\Phi a\Vert }_2^2+{\sigma }_n^2{\displaystyle \sum _{p=1}^P\frac{1}{{\sigma }_p^2}{\displaystyle \sum _{n=2}^N{\left(w_{p,n}\right)}^2}}\right]\hfill \\ & =\arg \min \left[{\Vert y-\Phi a\Vert }_2^2+{\sigma }_n^2{\displaystyle \sum _{p=1}^P\frac{1}{{\sigma }_p^2}{\displaystyle \sum _{n=2}^N{\left(a_{p,n}-a_{p,n-1}\right)}^2}}\right]\hfill \\ & =\arg \min \left[{\Vert y-\Phi a\Vert }_2^2+{\displaystyle \sum _{p=1}^P\frac{{\sigma }_n^2}{{\sigma }_p^2}{\Vert D{a}_p\Vert }_2^2}\right]\hfill \end{array}$$

(23)

From (16), (17) and (23), $\lambda$ can be expressed as:

$$\lambda ={\left[{\lambda }_1^{},{\lambda }_2^{},\cdots ,{\lambda }_P^{}\right]}^T={\left[\frac{{\sigma }_n^2}{{\sigma }_1^2},\frac{{\sigma }_n^2}{{\sigma }_2^2},\cdots \frac{{\sigma }_n^2}{{\sigma }_P^2}\right]}^T$$

(24)

From (24), $\lambda$ is related to the anisotropic characteristics of scattering centers, which we cannot know before the problem (16) is solved. Thus, we needed to jointly estimate $a$ and $\lambda$. Based on [16], Gamma distribution was used as the prior of ${\lambda }_p$ and its PDF can be expressed as:

$$P\left({\lambda }_p^{}\right)=\frac{{\lambda }_p^{\alpha -1}}{{\beta }^{\alpha }\Gamma \left(\alpha \right)}\exp \left(-\frac{{\lambda }_p^2}{\beta }\right)$$

(25)

where $\mathsf{\Gamma }(\cdot )$ is the Gamma function and $\alpha$ and $\beta$ are the hyper-parameters. The MAP estimation of $a$ and $\lambda$ can be achieved by solving the following problem:

$$\left\{\widehat{a},\widehat{\lambda }\right\}=\underset{a,\lambda }{\arg \min }\left({\Vert y-\Phi a\Vert }_2^2+{\sigma }_n^2{\displaystyle \sum _{p=1}^P\frac{1}{{\sigma }_p^2}{\Vert D{a}_p\Vert }_2^2}-\log \left({\displaystyle \sum _p{\lambda }_p^{}}\right)+{\beta }^{-1}{\displaystyle \sum _p{\lambda }_p^{}}\right)$$

(26)

As ${\sigma }_n^2$ can be estimated from the SAR image, we mainly focused on the estimation of $\sigma ={\left[{\sigma }_1^2,\cdots ,{\sigma }_p^2,\cdots {\sigma }_P^2\right]}^{\mathrm{T}}$. Based on (25) and (26), the MAP estimation of $a$ and $\sigma$ can be expressed as:

$$\left\{\widehat{a},\widehat{\sigma }\right\}=\underset{a,\sigma }{\arg \min }\left({\Vert y-\Phi a\Vert }_2^2+{\sigma }_n^2{\displaystyle \sum _{p=1}^P\frac{1}{{\sigma }_p^2}{\Vert D{a}_p\Vert }_2^2}-\log \left({\displaystyle \sum _p\frac{1}{{\sigma }_p^2}}\right)+{\beta }^{-1}{\sigma }_n^2{\displaystyle \sum _p\frac{1}{{\sigma }_p^2}}\right)$$

(27)

It is difficult to solve the problem (27) directly. In this study, the EM algorithm was employed to solve (27), which yielded an iterative algorithm. We use ${\sigma }_{t-1}$ ($t$ is the iteration number) to calculate ${\widehat{a}}_t$ by the following equation:

$${\widehat{a}}_t=\arg \min \left[{\Vert y-\Phi a\Vert }_2^2+{\Vert F_{w,t-1}a\Vert }_2^2\right]$$

(28)

where $F_{w,t-1}$ is expressed as

$$F_{w,t-1}=\left[\begin{array}{cccc}\frac{{\sigma }_n^2}{{\sigma }_{1,t-1}^2}D& & & \\ & \frac{{\sigma }_n^2}{{\sigma }_{2,t-1}^2}D& & \\ & & \ddots & \\ & & & \frac{{\sigma }_n^2}{{\sigma }_{P,t-1}^2}D\end{array}\right]$$

(29)

Then, ${\sigma }_t={\left[{\sigma }_{1,t}^2,\cdots ,{\sigma }_{p,t}^2,\cdots ,{\sigma }_{P,t}^2\right]}^T$ is calculated by the following equation:

$${\sigma }_{p,t}^2=\mathrm{var}\left(D{a}_{p,t}\right)$$

(30)

where $\mathrm{var}\left(D{a}_{p,t}\right)$ denotes the variance of $D{a}_{p,t}$. After ${\sigma }_t$ is obtained, the matrix $F_{w,t-1}$ is updated for the next iteration.

Problem (28) has an analytic solution. In the following, we derived the solution with the closed form. The conjugate gradient of $J\left(a_t\right)={\Vert y-\Phi a_t\Vert }_2^2+{\Vert F_{w,t-1}{a}_t\Vert }_2^2$ can be expressed as

$$\frac{\partial J\left(a_t\right)}{\partial a_t^{*}}=-{\Phi }^{H}y+{\Phi }^H\Phi a_t+F_{w,t-1}^T{F}_{w,t-1}{a}_t$$

(31)

By setting $\frac{\partial J\left(a_t\right)}{\partial a_t^{*}}=0$, the minimizer of (28) can be expressed as

$$\left({\Phi }^H\Phi +F_{w,t-1}^T{F}_{w,t-1}\right){\widehat{a}}_t={\Phi }^{H}y$$

(32)

Equation (32) is the closed form of the solution to (28), so (28) can be solved directly. Consequently, the proposed method can be implemented with high efficiency.

To sum up, the IRWTR algorithm consists of two steps. In the first step, we calculate $a_t$ by (28) and ${\sigma }_{t-1}$. In the second step, a new ${\sigma }_t$ is calculated via (30) to update the regularization parameters. This procedure will be repeated until the difference between ${\sigma }_{t-1}$ and ${\sigma }_t$ in two successive iterations is smaller than a small threshold.

3.3. The Restoration of the Incomplete Images of Certain Distributed Scattering Centers

In this subsection, we discuss the method to restore the incomplete imaging results of the distributed scattering centers. We have mentioned that some points in the SAR image are actually the ends of distributed scattering centers. With the aspect-dependent amplitudes extracted by the IRWTR, we can distinguish these ends from the local scattering centers. Furthermore, we need to know which two points belong to the same distributed scattering center and restore the incomplete imaging result.

Firstly, we distinguished the ends of the distributed scattering centers from the local scattering centers based on the variance of the aspect-dependent amplitudes. A threshold was set and the scattering centers’ aspect-dependent-amplitude variances below the threshold will be identified to be local scattering centers, while the rest were identified to be the ends of distributed scattering centers.

Secondly, these ends were arranged in pairs, of which the procedure is referred to as the matching procedure. The matching procedure exploits the ASCM and the minimization of least squares error criterion. With the aspect-dependent amplitudes extracted by the IRWTR, the position and the relative amplitude, the estimated backscatter of the $p$th scatterer at $f_c$ can be expressed as:

$$E_p^{\mathrm{est}}\left(f_c^{}\right)={\left[E_p^{\mathrm{est}}\left(f_c^{},{\phi }_1^{}\right),\cdots ,E_p^{\mathrm{est}}\left(f_c^{},{\phi }_n^{}\right),\cdots ,E_p^{\mathrm{est}}\left(f_c^{},{\phi }_N^{}\right)\right]}^{\mathrm{T}}$$

(33)

where $E_p^{\mathrm{est}}\left(f_c^{},{\phi }_n^{}\right)$ is calculated by (6). If the $i$ and $j$th scattering center are the two ends of the same distributed scattering center, the least squares error $e_{i,j}^{\mathrm{LS}}$, expressed as the following equation, will be minimized:

$$e_{i,j}^{\mathrm{LS}}={\Vert E_i^{\mathrm{est}}\left(f_c\right)+E_j^{\mathrm{est}}\left(f_c\right)-E_{i,j}^{\mathrm{est}}\left(f_c\right)\Vert }_2^2$$

(34)

where $E_{i,j}^{\mathrm{est}}\left(f_c\right)$ represents the backscatter of the distributed scattering center whose ends are at the $i$th and $j$th scatterers’ positions. The length $L_{i,j}^{\mathrm{est}}$, orientation angle ${\overline{\phi }}_{0,i,j}^{\mathrm{est}}$, position $\left(x_{i,j}^{\mathrm{est}},y_{i,j}^{\mathrm{est}}\right)$ and relative amplitude $A_{i,j}^{\mathrm{est}}$ of this distributed scattering center are calculated by the following equations:

$$L_{i,j}^{\mathrm{est}}=\sqrt{{\left(x_i^{\mathrm{est}}-x_j^{\mathrm{est}}\right)}^2+{\left(y_i^{\mathrm{est}}-y_j^{\mathrm{est}}\right)}^2}$$

(35)

$${\overline{\phi }}_{0,i,j}^{\mathrm{est}}={\tan }^{-1}\left(\frac{y_i^{\mathrm{est}}-y_j^{\mathrm{est}}}{x_i^{\mathrm{est}}-x_j^{\mathrm{est}}}\right)$$

(36)

$$x_{i,j}^{\mathrm{est}}=\frac{x_i^{\mathrm{est}}+x_j^{\mathrm{est}}}{2}$$

(37)

$$y_{i,j}^{\mathrm{est}}=\frac{y_i^{\mathrm{est}}+y_j^{\mathrm{est}}}{2}$$

(38)

$$A_{i,j}^{\mathrm{est}}=\frac{A_i^{\mathrm{est}}+A_j^{\mathrm{est}}}{2}$$

(39)

where $\left(x_i^{\mathrm{est}},y_i^{\mathrm{est}}\right)$, $\left(x_j^{\mathrm{est}},y_j^{\mathrm{est}}\right)$ are the estimated positions of the $i$th and $j$th scatterer. $A_i^{\mathrm{est}}$ and $A_j^{\mathrm{est}}$ are the estimated relative amplitides of the $i$th and $j$th scatterer. The backscatter $E_{i,j}^{\mathrm{est}}\left(f_c\right)$ can be calculated by (1) and the parameters in (35)–(39).

At the beginning of the matching procedure, we chose one end and found the other end that was paired with it by minimizing (34). Then, this procedure was repeated among the residual ends. The matching procedure was stopped until all the ends were arranged in pairs.

Thirdly, we restored the incomplete imaging results of the distributed scattering centers based on the matching result. We added a signal to the original data in the frequency-aspect domain. The imaging result of this signal was the flank of a distributed scattering center. For the distributed scattering center whose ends are the $i$th and $j$th scatterer, we added the signal expressed as the following equation to the original data:

$$E_{i,j}^{}\left(f_k^{},{\phi }_n^{},x_{i,j}^{\mathrm{est}},y_{i,j}^{\mathrm{est}}\right)=A_{i,j}^{\mathrm{est}}\sin \mathrm{c}\left(\left(\frac{2\pi L_{i,j}^{\mathrm{est}}}{c}\right)\left(f_k^{}\sin \left({\phi }_n^{}-{\overline{\phi }}_{0,i,j}^{\mathrm{est}}\right)+f_c\sin {\overline{\phi }}_{0,i,j}^{\mathrm{est}}\right)\right)$$

(40)

where the parameters are calculated by (35)–(39). After we fixed all the incomplete imaging results of the distributed scattering centers by using (40), a BP algorithm was performed to yield the improved SAR image.

Based on the aforementioned contents, we present the flowchart of the whole algorithm in Figure 6. In the first part, the extraction of the scattering centers’ positions and relative amplitudes is performed. The relative amplitudes and positions are used to form the matrix $\Phi$ in (9). In the second part, the IRWTR is utilized to estimate the aspect-dependent amplitudes of the scattering centers. In the third part, the incomplete imaging results of the distributed scattering centers are restored based on the extracted anisotropic characteristics.

In the next section, the processing results by the proposed method will be shown.

4. Experimental Results

In this section, we apply the proposed algorithm to the simulated data and the corresponding results are shown. In Section 4.1, we show the IRWTR processing result of the data simulated by MATLAB.

4.1. Results of Matlab Simulated Data and Analysis

The processing result of Matlab simulated data will be shown in this section. There are three scattering centers whose positions are $\left(0m,1m\right)$ $\left(0m,0m\right)$ and $\left(1.5m,0m\right)$, respectively, in the imaging scene. Their amplitude functions of azimuth were a constant, Sinc function and hamming window function. The radar transmitted step sweep frequency signal with a start frequency of 9.5 GHz and bandwidth of 1 GHz. The number of frequency signals was 32. The look angles ranged from −22.5° to 22.5°. The number of the look angles was 64.

Two methods were used to process the simulated data: solving (13) and the IRWTR algorithm. The processing was undertaken under different noise conditions (10 dB and 0 dB) and the results are shown in Figure 7 and Figure 8, where the blue dashed line shows the real value and the red solid line shows the estimation.

In Figure 7, the performance of the estimation of the Sinc-like amplitude function was the best. In contrast, the other estimations were inaccurate, for the regularization parameter was not large enough. This phenomenon is more significant in Figure 8. According to estimation of the Sinc-like amplitude function in Figure 8b, the regularization parameter was larger than the requirement, for the estimation was too smooth and deviated from the real values. In contrast, the regularization parameters for the rest amplitude functions were smaller than the requirement, which made the estimations oscillate. Compared with the regularization with a single regularization parameter, the MRP regularization provided better estimations. Both in Figure 7 and Figure 8, the estimations obtained via IRWTR were more close to the real values.

To evaluate the performance of different methods under different conditions of signal-to-noise ratio (SNR) SNR, a Monte Carlo simulation was used, in which 50 trials were performed for each SNR. In each trial, the relative mean squared error (RMSE) [17] was calculated. The statistical analysis is shown in Figure 9. In Figure 9, the performance of the IRWTR was more superior as SNR decreased.

4.2. The Results of Electromagnetic Simulation Data and Analysis

In this section, we apply the proposed processing procedure to the data simulated by FEKO software and the processing results will be shown.

The imaging scene is shown in Figure 10, which included two metal rods and three metal balls. The radar transmitted a step sweep frequency signal with a start frequency of 8 GHz, bandwidth 4 GHz. The azimuth start and end angle were −22.5° and 22.5°. We applied the BP imaging algorithm to the simulation data and the result is shown in Figure 11.

In Figure 11, there are only seven points. In fact, we knew that there were two distributed scattering centers in the scene from Figure 10. However, because their specific orientations, we could not recognize these two distributed scattering centers from Figure 11. We applied the proposed method in Section 3.2 to the data and extract the aspect-dependent amplitudes of all points. The estimated aspect-dependent amplitudes are presented in Figure 12. According to Figure 12, we can find that the amplitudes of certain scatterers were almost independent of the look angle. These scatterers with slowly varying amplitudes were the imaging results of the metal balls, which were the actual local scattering centers. The variations in the amplitude responses of the rest were more remarkable. The degrees of variation in amplitude responses could help us distinguish the ends of distributed scattering centers from the local scattering centers. Based on the azimuth-dependent scattering behaviors shown in Figure 12, the ends of distributed scattering centers are indicated by red circles in Figure 11.

The estimated aspect-dependent amplitudes of the scatterer 1, 2, 4, and 7 were compared with the ones calculated by FEKO. The comparisons are shown in Figure 13.

From Figure 13, it can be found that the estimated amplitudes corresponded closely with the ones calculated by FEKO. Moreover, we arranged the ends in pairs based on the extracted amplitudes. The signals to form the images of the flanks were constructed based on the arrangement and added to the original return signals. The imaging result of the reproduced data is shown in Figure 14. In Figure 14, the distributed scattering centers were easier to be distinguished from the local scattering centers.

4.3. The Results of the Backhoe Data and Discussion

In this subsection, we apply the proposed method to the backhoe data. The backhoe model is shown in Figure 15.

The corresponding processing results by the proposed method are presented in Figure 16b,e. The original imaging result of the backhoe with the look angles ranging from 25°to 70° is presented in Figure 16a. The portion within the white box in Figure 16a is shown in Figure 16d. The imaging result with the look angles ranging from −10°to 35° is shown in Figure 16c and the corresponding enlargement is shown in Figure 16f. In Figure 16a,d, some distributed scattering centers cannot be seen, which can be found in Figure 16c,f. The proposed algorithm was applied to the data with look angles ranging from 25° to 70° to restore the incomplete images of distributed scattering centers. The restored image is shown in Figure 16b,e. According to Figure 16b,e, some incomplete images of distributed scattering centers were restored.

However, in the restored image, we also found that some images of the distributed scattering centers were not restored correctly. According to the inaccurate restorations, we discussed some problems and limitations of the proposed method. In the discussion part, we mainly focus on the processing result of the bucket of the backhoe, whose images are within the white boxes in Figure 16.

The black and red lines in Figure 15 mark the distributed scatterers whose images are not complete in Figure 16d. The red lines in Figure 15 denote the distributed scatterers with the correctly restored images in Figure 16e, corresponding to the areas within the red dashed boxes. The distributed scatterers that were not restored correctly were marked by the black lines on the backhoe model in Figure 15. These scatterers correspond to the areas within the black dashed boxes. The two distributed scatterers marked by the black lines on the model were inside the bucket under the condition that the look angles ranged from 25° to 70°, since we can only find the images of their ends on one side in Figure 16e. Consequently, although the algorithm detected the two ends, it cannot perform the correct pairing, since the ends that can be correctly paired with them cannot be illuminated by the signals. However, the proposed algorithm has to find the ends to be paired with them. As a result, the wrong pairing occured. Another phenomenon we did not consider in this algorithm is the interactions [10] among the scatterers. Interactions among the scatterers are very common in a complex scene, especially for the bucket of the backhoe which can be regarded as a cavity. The interactions may give rise to the images of scatterers with wrong positions. This may influence the anisotropic characteristics extraction and the pairing.

5. Conclusions

A novel method for extracting anisotropic characteristics of scattering centers and the feature enhancement in the wide-angle SAR imagery was proposed in this paper. The main target of this study was to restore the incomplete images of the distributed scattering centers in a wide-angle SAR image based on the anisotropic characteristics. In the proposed scheme, firstly, the anisotropic characteristics were extracted by the proposed IRWTR algorithm. Secondly, the anisotropic characteristics were used to determine if a point in the SAR image was actually the image of an end of a distributed scattering center. Thirdly, these point-like images of the ends of distributed scattering centers will be lined up in pairs based on their anisotropic characteristics extracted by IRWTR.