A Sequential Optimization Calibration Algorithm for Near-Field Source Localization

Abstract

This paper considers the near-field source location problem for a nonuniform linear array (non-ULA) in the presence of sensor gain and phase errors. A sequential optimization calibration method is proposed to simultaneously estimate the gain and phase errors as well as the locations of calibration sources involving the ranges and the azimuths by exploiting some imprecise a-priori knowledge of calibration sources. At each iteration of the proposed method, the source locations, and the gain and phase errors are obtained iteratively. Finally, at the analysis stage, we evaluate the effectiveness of the proposed technique using some numerical simulations. Results show that the proposed algorithm shares the capability to jointly estimate the source locations and the errors.

1. Introduction

Source localization with sensor arrays has broad applications in aerospace, navigation and wireless acoustic sensor network societies [1,2,3,4,5,6]. Many effective algorithms for source locations, such as the maximum-likelihood (ML) [7], the multiple signal classification (MUSIC) [8] and the minimum-variance distortionless response (MVDR) algorithms [9], have been developed over the years. Recently, some source location algorithms based on cumulative sum and steered-response power [4,5,10,11,12] are also presented. It is worth pointing out that most of these algorithms suppose the array manifold is perfectly available. However, in practice, the array manifold is often affected by unknown array characteristics such as the sensor gain and phase errors as well as the unknown mutual coupling [13,14,15], thus resulting in the performance deterioration of these algorithms. For example, in [14], the precision of the direction-of-arrival (DOA) estimation would decrease due to the existence of the unknown mutual coupling. In [15], the accuracy of measuring acoustic intensity would degrade in the presence of the sensor gain and phase errors.

Some calibration algorithms for the gain and phase errors have been investigated. Specifically, in [16,17], the eigendecomposition method was proposed to derive the DOAs with unknown gain and phase errors. In [18], using the null characteristic of the MUSIC spectrum, a calibration technique was proposed for the sensor gain and phase uncertainties as well as location errors. In [19], the maximum a posteriori (MAP) method was presented to estimate the DOAs and the perturbations simultaneously. In [20], an eigenstructure method with the aim of simultaneously estimating the DOA and gain-phase errors without joint iteration, was proposed. In [21], based on different data models, two new estimation algorithms were presented for uniform linear array (ULA) to estimate sensor gain and phase errors. In [22], exploiting the subspace principle, the estimations of sensor gain and phase errors were addressed in subarrays-based linear sparse arrays. In [23], a new method based on the eigendecomposition of the Hadamard product of the covariance matrix and its conjugate was investigated requiring no a priori knowledge of calibration sources. In [24], the authors presented an estimation of signal parameters via rotational invariance technique (ESPRIT)-like method that can simultaneously estimate DOA, as well as the gain and phase errors in the uncalibrated portion of the ULA in close form.

Summarizing, the aforementioned works focus on the far-field scenario, where the sources are located far enough from the sensor arrays and only the DOAs are of interest. In the case of the near-field scenario, the steering vector involves the knowledge of both the azimuths and the ranges. Some works, such as [25,26,27,28,29,30,31], have been presented to obtain the ranges and DOAs of the sources but with perfect knowledge of the array manifold. However, in the presence of the sensor gain and phase errors, the estimation performance would be degraded significantly. There is very little work with investigating the near-field source location problem for a nonuniform linear array (non-ULA) in the presence of sensor gain and phase errors. In [32], the passive localization of near-field sources with partly calibrated subarray-based arrays was studied. However, the proposed method cannot apply in arbitrary arrays. We can employ the traditional active-calibration methods to estimate the gain and phase errors no matter if they are far-field targets or near-field targets requiring the exact information (azimuth or range) of the calibration source. However, from a practical point of view, the exact information cannot be available.

In this paper, we consider the near-field source location problem for a non-ULA in the presence of sensor gain and phase errors. We present a sequential optimization calibration method to simultaneously estimate the error parameters and the locations of calibration sources under minimum variance estimation criterion based on some imprecise a priori knowledge of calibration sources. Specifically, we first obtain the estimates of source azimuths, the gain and phase errors using given ranges. Second, the source ranges and the error parameters are derived according to the known azimuths. Finally, we continue the iterative procedure until convergence. Simulation results highlight that the proposed algorithm shares the capability of the joint estimates of the source locations and the error parameters. In particular, it is worth pointing out that the proposed method exploits the imprecise knowledge of calibration sources to estimate gain and phase errors. As a consequence, it is more effective than the active-calibration methods owing to combing the self-calibration and active-calibration techniques.

The rest of the paper is organized as follows. In Section 2, we formulate the signal model. Section 3 describes a self-calibrating method for near-field sources with unknown gain and phase errors. In Section 4, we provide some simulations to illustrate the effectiveness of the proposed method. Finally, in Section 5, some conclusions are derived.

2. Signal Model

Assume that a non-ULA contains M sensors, which are placed along the x-axis at $d_0,d_1,\cdots ,d_{M-1}(d_0=0)$ with unequal spacing, respectively, as shown in Figure 1. There are P uncorrelated narrow-band near-field sources impinging on the SLA with azimuth and range pairs $({\theta }_p,r_{0p})$, $p=1,2,\dots ,P$, where ${\theta }_p$ is the azimuth of the p-th source deviating normal direction of array reference point and $r_{0p}$ is the distance between the p-th source and array reference point. Let $s_p\left(n\right)$, $n=1,2,\dots ,L$, be the source waveforms. Similar to [33], the received signal vector of the non-ULA at the nth snapshot can then be expressed as,

$$\begin{array}{cc}\hfill \mathit{x}\left(n\right)& =\sum _{p=1}^P\mathbf{\Gamma }\mathit{a}({\theta }_p,r_{0p})s_p\left(n\right)+v\left(n\right)\hfill \\ & =\mathbf{\Gamma }\mathit{A}\mathit{s}\left(n\right)+v\left(n\right)n=1,2,\cdots ,L,\hfill \end{array}$$

(1)

where

• $\mathit{x}\left(n\right)={[x_0\left(n\right),x_1\left(n\right),\cdots ,x_{M-1}\left(n\right)]}^T$ is the measurement signal vector;
• $v\left(n\right)={[v_0\left(n\right),v_1\left(n\right),\cdots ,v_{M-1}\left(n\right)]}^T$ is an independent and identically distributed complex circular zero-mean Gaussian random vector with covariance matrix ${\sigma }^2{\mathbf{I}}_M$, while ${\mathbf{I}}_M$ is the M-dimensional identity matrix;
• $\mathit{A}=[\mathit{a}({\theta }_1,r_{01}),\mathit{a}({\theta }_2,r_{02}),\cdots ,\mathit{a}({\theta }_P,r_{0P})]$ is the nominal $M\times P$ steering matrix, the p-th column is

  $$\begin{array}{c}\mathit{a}({\theta }_p,r_{0p})=\left[1,\frac{r_{0p}}{r_{1p}}{e}^{-j\frac{2\pi \Delta r_{1p}}{{\lambda }_p}},\cdots ,\right.\hfill \\ {\left.\frac{r_{0p}}{r_{\left(M-1\right)p}}{e}^{-j\frac{2\pi \Delta r_{(M-1)p}}{{\lambda }_p}}\right]}^T,\hfill \end{array}$$

  (2)

  where ${\lambda }_p$ is the wavelength of the pth source and $r_{mp}$ denotes distance of the pth signal source to mth sensor, $p=1,2,\cdots ,P$, $m=0,1,\cdots ,M-1$. $\Delta r_{mp}$ is the relative distance between $r_{0p}$ and $r_{mp}$, which can be derived by geometrical relationship,

  $$\begin{array}{cc}\hfill \Delta r_{mp}& =r_{0p}-r_{mp}\hfill \\ & =r_{0p}-\sqrt{{(r_{0p}sin{\theta }_p-d_m)}^2+{(r_{0p}cos{\theta }_p)}^2}.\hfill \end{array}$$

  (3)
  To simplify the notation, we write $r_{0p}$ into $r_p$.
• $\mathbf{\Gamma }=\mathbf{diag}\{{\rho }_0{e}^{j{\phi }_0},{\rho }_1{e}^{j{\phi }_1},\cdots ,{\rho }_{M-1}{e}^{j{\phi }_{M-1}}\}$ is the error matrix of the array gain and phase, where parameters ${\rho }_m$ and ${\phi }_m$ are the gain and the phase errors associated with the m-th sensor, respectively.

3. Near-Field Calibration Method

Unlike [17], which investigated the calibration algorithm for the far-field sources where only the DOAs are of interest, this paper focuses on the case of the near-field sources involving these considerations of range and azimuth. To this end, we present a sequential optimization calibration technique to estimate the gain and phase errors, and the locations of radiating sources, simultaneously. Specifically, fixed the ranges of radiating sources, we first estimate the azimuths, the gain and phase errors, and then obtain ranges and the error parameters exploiting the estimated azimuths, and continue the procedure until convergence.

3.1. Joint Estimations of Error Matrix and Azimuthes with Known Ranges

In this subsection, we extend the method [16] to simultaneously estimate ${\mathbf{\Gamma }}_{\theta }$ and ${\left\{{\theta }_p\right\}}_{p=1}^P$ using estimated ranges ${\left\{{\widehat{r}}_p\right\}}_{p=1}^P$, where ${\mathbf{\Gamma }}_{\theta }$ is the array gain and phase error matrix under known ranges and ${\widehat{r}}_p$ denotes the estimated value of $r_p$. Precisely, given a ${\mathbf{\Gamma }}_{\theta }$, we look for P peak values corresponding to ${\left\{{\theta }_p\right\}}_{p=1}^P$, with respect to one-dimensional MVDR spectrum [9], given by,

$$\begin{array}{c}P\left(\theta \right|{\widehat{r}}_p,{\mathbf{\Gamma }}_{\theta })=\frac{1}{{\mathit{a}}^H(\theta ,{\widehat{r}}_p){\mathbf{\Gamma }}_{\theta }^H{\mathit{R}}_x^{-1}{\mathbf{\Gamma }}_{\theta }\mathit{a}(\theta ,{\widehat{r}}_p)},\hfill \\ p=1,2,\cdots ,P.\hfill \end{array}$$

(4)

Using the ${\left\{{\theta }_p\right\}}_{p=1}^P$, we estimate ${\mathbf{\Gamma }}_{\theta }$ based on minimum variance estimation criterion [9] and proceed the iteration procedure until convergence. Specifically, given ${\left\{{\theta }_p\right\}}_{p=1}^P$, we extend the proposed self-calibration algorithm [16] and develop a cost function computed as,

$$J\left({\mathbf{\Gamma }}_{\theta }\right)=\sum _{p=1}^P{\mathit{a}}^H({\theta }_p,{\widehat{r}}_p){{\mathbf{\Gamma }}_{\theta }}^H{\mathit{R}}_x^{-1}{\mathbf{\Gamma }}_{\theta }\mathit{a}({\theta }_p,{\widehat{r}}_p),$$

(5)

where $\mathit{a}({\theta }_p,{\widehat{r}}_p)$ denotes the steering manifold in Label (2) with known range ${\widehat{r}}_p$ and ${\mathit{R}}_x$ is covariance matrix computed as,

$${\mathit{R}}_x=\mathbb{E}\left[\mathit{x}{\mathit{x}}^H\right]={\mathbf{\Gamma }}_{\theta }\mathit{A}{\mathit{R}}_s{\mathit{A}}^H{{\mathbf{\Gamma }}_{\theta }}^H+{\sigma }^2{\mathit{I}}_M,$$

(6)

while ${\mathit{R}}_s=\mathbb{E}\left[\mathit{s}{\mathit{s}}^H\right]$, $\mathbb{E}[·]$ and ${(·)}^H$ represent the expectation and the Hermitian transpose operation, respectively. Note that the number P of signal source can be estimated based on Schwartz and Rissanen (MDL) criteria [34] and ${\mathit{R}}_x$ can be estimated by the sample covariance matrix (i.e., ${\mathit{R}}_x={\displaystyle \sum _{n=1}^N}\mathit{x}\left(n\right){\mathit{x}}^H\left(n\right)$, where N is the number of snapshots).

Further, let

$${\mathbf{\Gamma }}_{\theta }\mathit{a}({\theta }_p,{\widehat{r}}_p)=\tilde{\mathit{a}}({\theta }_p,{\widehat{r}}_p)\mathit{\delta },$$

(7)

in which

$${\mathbf{\Gamma }}_{\theta }=\mathbf{diag}\left(\mathit{\delta }\right),$$

(8)

and $\tilde{\mathit{a}}({\theta }_p,r_p)$ is a diagonal matrix, given by,

$$\tilde{\mathit{a}}({\theta }_p,{\widehat{r}}_p)=\mathbf{diag}\left\{\mathit{a}({\theta }_p,{\widehat{r}}_p)\right\},$$

(9)

while $\mathbf{diag}(·)$ denotes diagonal matrix formed by the entries of the vector.

Submitting Label (7) into Label (5), we have

$$J\left({\mathbf{\Gamma }}_{\theta }\right)={\mathbf{\delta }}^H\left\{\sum _{p=1}^P\tilde{\mathit{a}}{({\theta }_p,{\widehat{r}}_p)}^H{\mathit{R}}_x^{-1}\tilde{\mathit{a}}({\theta }_p,{\widehat{r}}_p)\right\}\mathbf{\delta }.$$

(10)

Hence, assuming that ${\rho }_0=1$,${\phi }_0=0$, the optimization problem ${\mathcal{P}}_1$ accounting for the constraint ${\mathbf{\delta }}^H\mathit{w}=1$ where $\mathit{w}={[1,0,\dots ,0]}^T$, can be written as,

$${\mathcal{P}}_1\left\{\begin{array}{cc}{\displaystyle \underset{{\mathbf{\Gamma }}_{\theta }}{min}}\hfill & J\left({\mathbf{\Gamma }}_{\theta }\right)\hfill \\ \mathrm{s}.\mathrm{t}.\hfill & {\mathbf{\delta }}^H\mathit{w}=1,\hfill \\ & {\mathbf{\Gamma }}_{\theta }=\mathbf{diag}\left(\mathbf{\delta }\right).\hfill \end{array}\right.$$

(11)

Employing the Lagrangian multiplier method, the optimal solution to ${\mathcal{P}}_1$ is derived as:

$$\mathbf{\delta }=\frac{{\mathit{Q}}^{-1}\mathit{w}}{{\mathit{w}}^T{\mathit{Q}}^{-1}\mathit{w}},$$

(12)

where

$$\mathit{Q}=\sum _{p=1}^P\tilde{\mathit{a}}{({\theta }_p,{\widehat{r}}_p)}^H{\mathit{R}}_x^{-1}\tilde{\mathit{a}}({\theta }_p,{\widehat{r}}_p).$$

(13)

Summarizing, the procedure of joint estimating ${\mathbf{\Gamma }}_{\theta }$ and source azimuths ${\left\{{\theta }_p\right\}}_{p=1}^P$ is summarized in Algorithm 1:

Algorithm 1 Algorithm for the joint estimations of ${\mathbf{\Gamma }}_{\theta }$ and $\theta$.

Input: $(d_0,d_1,\dots ,d_{M-1})$, ${\left\{r_p\right\}}_{p=1}^P$, w, P, ${\mathit{R}}_x$. Output: ${\left\{{\theta }_p\right\}}_{p=1}^P$ and ${\mathbf{\Gamma }}_{\theta }$. 1: For $k=0$ and ${\mathbf{\Gamma }}^{\left(0\right)}={\mathbf{I}}_M$; 2: Find the azimuths ${\left\{{\theta }_p^{\left(0\right)}\right\}}_{p=1}^P$ by searching for P highest peak using Label (4); 3: Compute $J_0$ by Label (10); 4: $k:=k+1$; 5: Compute ${\mathit{Q}}^{\left(k\right)}$ by Label (13); 6: Compute ${\mathbf{\delta }}^{\left(k\right)}$ by Label (12); 7: Construct ${\mathbf{\Gamma }}^{\left(k\right)}=\mathbf{diag}\left\{{\mathbf{\delta }}^{\left(k\right)}\right\}$; 8: Find the azimuths ${\left\{{\theta }_p^{\left(k\right)}\right\}}_{p=1}^P$ by searching for P highest peak using Label (4); 9: Compute $J_k$ by Label (10); 10: If $|J_k-J_{k-1}|>\kappa$, where $\kappa$ is a user selected parameter to control convergence, then $k=k+1$, return to step 4. Otherwise, stop and output ${\mathbf{\Gamma }}_{\theta }={\mathbf{\Gamma }}^{\left(k\right)}$, ${\left\{{\theta }_p\right\}}_{p=1}^P={\left\{{\theta }_p^{\left(k\right)}\right\}}_{p=1}^P$.

3.2. Joint Estimations of Error Matrix and Ranges with Known Azimuths

In this subsection, we focus on jointly estimating ${\mathbf{\Gamma }}_r$ and ${\left\{r_p\right\}}_{p=1}^P$ with known azimuths ${\left\{{\widehat{\theta }}_p\right\}}_{p=1}^P$, where ${\mathbf{\Gamma }}_r$ is the array gain and phase error matrix under known azimuths and ${\widehat{\theta }}_p$ denote the estimated value of ${\theta }_p$. Specifically, given a ${\Gamma }_r$, we first find ${\left\{r_p\right\}}_{p=1}^P$ by searching for one-dimension MVDR spectrum, given by

$$\begin{array}{c}P\left(r\right|{\widehat{\theta }}_p,{\mathbf{\Gamma }}_r)=\frac{1}{{\mathit{a}}^H({\widehat{\theta }}_p,r){\mathbf{\Gamma }}_r^H{\mathit{R}}_x^{-1}{\mathbf{\Gamma }}_r{\mathit{a}}^H({\widehat{\theta }}_p,r)},\hfill \\ p=1,2,\cdots ,P.\hfill \end{array}$$

(14)

Then, according to the knowledge ${\left\{r_p\right\}}_{p=1}^P$, we employ the same above procedure to derive ${\mathbf{\Gamma }}_r$ and continue the procedure until convergence.

Next, we focus on the derivation of ${\mathbf{\Gamma }}_r$ with obtained ${\left\{r_p\right\}}_{p=1}^P$. Similarly, the optimal problem ${\mathcal{P}}_2$ can be denoted as

$${\mathcal{P}}_2\left\{\begin{array}{cc}{\displaystyle \underset{{\mathbf{\Gamma }}_r}{min}}\hfill & J\left({\mathbf{\Gamma }}_r\right),\hfill \\ \mathrm{s}.\mathrm{t}.\hfill & {\mathbf{\delta }}^H\mathit{w}=1,\hfill \\ & {\mathbf{\Gamma }}_r=\mathbf{diag}\left(\mathbf{\delta }\right),\hfill \end{array}\right.$$

(15)

where

$$J\left({\mathbf{\Gamma }}_r\right)={\mathbf{\delta }}^H\left\{\sum _{p=1}^P\tilde{\mathit{a}}{({\widehat{\theta }}_p,r_p)}^H{\mathit{R}}_x^{-1}\tilde{\mathit{a}}({\widehat{\theta }}_p,r_p)\right\}\mathbf{\delta },$$

(16)

and

$$\tilde{\mathit{a}}({\widehat{\theta }}_p,r_p)=\mathbf{diag}\left\{\mathit{a}({\widehat{\theta }}_p,r_p)\right\}.$$

(17)

We can derive the solution of ${\mathcal{P}}_2$ as

$$\mathbf{\delta }=\frac{{\mathit{Q}}^{-1}\mathit{w}}{{\mathit{w}}^T{\mathit{Q}}^{-1}\mathit{w}},$$

(18)

where

$$\mathit{Q}=\sum _{p=1}^P\tilde{\mathit{a}}{({\widehat{\theta }}_p,r_p)}^H{\mathit{R}}_x^{-1}\tilde{\mathit{a}}({\widehat{\theta }}_p,r_p).$$

(19)

Finally, Algorithm 2 summarizes the procedure of jointly estimating ${\mathbf{\Gamma }}_r$ and source ranges ${\left\{r_p\right\}}_{p=1}^P$.

Algorithm 2 Algorithm for the joint estimations of ${\mathbf{\Gamma }}_r$ and $r_p$

Input: $(d_0,d_1,\dots ,d_{M-1})$, ${\left\{{\theta }_p\right\}}_{p=1}^P$, w, P, ${\mathit{R}}_x$; Output: ${\left\{r_p\right\}}_{p=1}^P$ and ${\mathbf{\Gamma }}_r$; 1: For $k=0$ and ${\mathbf{\Gamma }}^{\left(0\right)}={\mathbf{I}}_M$; 2: Find the ranges ${\left\{r_p^{\left(0\right)}\right\}}_{p=1}^P$ by searching for P highest peak using Label (14); 3: Compute $J_0$ by Label (16); 4: $k:=k+1$; 5: Compute ${\mathit{Q}}^{\left(k\right)}$ by Label (19); 6: Compute ${\mathbf{\delta }}^{\left(k\right)}$ by Label (18); 7: Construct ${\mathbf{\Gamma }}^{\left(k\right)}=\mathbf{diag}\left\{{\mathbf{\delta }}^{\left(k\right)}\right\}$; 8: Find the ranges ${\left\{r_p^{\left(k\right)}\right\}}_{p=1}^P$ by searching for P highest peak using (14); 9: Compute $J_k$ by Label (16); 10: If $|J_k-J_{k-1}|>\kappa$, then $k=k+1$, return to step 4. Otherwise, stop and output ${\mathbf{\Gamma }}_r={\mathbf{\Gamma }}^{\left(k\right)}$, ${\left\{r_p\right\}}_{p=1}^P={\left\{r_p^{\left(k\right)}\right\}}_{p=1}^P$.

3.3. Joint Estimations of Error Matrix, Azimuths and Ranges

In this subsection, the proposed iteration procedure for jointly estimating the source range and azimuth $({\theta }_p,r_p)$, $p=1,2,\cdots ,P$, as well as the gain and phase error matrix $\mathbf{\Gamma }$ is summarized in Algorithm 3. It is worth mentioning that the calibration algorithm based on eigenstructure methods in [17] can also be used to calibrate the near-field sources. However, each iteration of the calibration algorithm requires to search a two-dimensional pseudo-specturm resulting in a large computational complexity. In particular, each iteration of the proposed algorithm provides a computationally efficient calibration method through converting a two-dimensional spectrum calibration problem into two one-dimensional spectrum calibration problems, which significantly decreases the computational burden.

Finally, we point out that the proposed algorithm may not converge to an optimal solution since both optimization problems (4) and (14) are not convex. The proposed iterative algorithm ensures obtaining a quality suboptimal solution as the approximate solution of optimal solution (please see

Table 1. Theoretical and average estimated values of the gain and phase errors for 12 sensors.

Sensor		1	2	3	4	5	6	7	8	9	10	11	12
Gain error	$\rho$	1.000	0.892	1.129	1.012	1.151	1.171	1.175	0.982	1.155	0.908	1.141	1.145
	$\widehat{\rho }$	1.000	0.878	1.115	0.999	1.134	1.155	1.161	0.966	1.136	0.894	1.122	1.128
Gain error	$\phi \left(\mathrm{deg}\right)$	0	−4.067	−12.385	1.270	8.214	17.042	−1.542	11.344	14.826	−4.152	4.769	−4.553
	$\widehat{\phi }\left(\mathrm{deg}\right)$	0	−4.093	−12.249	1.453	8.285	17.232	−1.312	11.868	15.387	−3.443	5.580	−3.674

).

Algorithm 3 Algorithm for the joint estimation of $\mathbf{\Gamma }$ and ${\left\{{\theta }_p\right\}}_{p=1}^P,{\left\{r_p\right\}}_{p=1}^P$

Input: $(d_0,d_1,\dots ,d_{M-1})$, w, ${\mathbf{\Gamma }}_0$, ${\mathit{R}}_x$; Output: $\mathbf{\Gamma }$ and $({\theta }_p,r_p)$, $p=1,2,\dots ,P$; 1: For $m=0$ and estimate the number P of signals; 2: Find the P peaks corresponding to locations $({\theta }_p^{\left(0\right)},r_p^{\left(0\right)})$, $p=1,2,\dots ,P$ by searching for two-dimension MVDR spectrum $P_c(\theta ,r|{\mathbf{\Gamma }}^{\left(0\right)})=1/{\mathit{a}}^H(\theta ,r){\left({\mathbf{\Gamma }}_0\right)}^H{\mathit{R}}_x^{-1}{\mathbf{\Gamma }}_0\mathit{a}(\theta ,r)$; 3: $m:=m+1$; 4: Estimate ${\{{\theta }_p^{(m+1)}\}}_{p=1}^P$ and ${\mathbf{\Gamma }}_{\theta }$ by Algorithm 1 using ${\{r_p^{\left(m\right)}\}}_{p=1}^P$; 5: Estimate ${\{r_p^{(m+1)}\}}_{p=1}^P$ and ${\mathbf{\Gamma }}_r$ by Algorithm 2 using ${\{{\theta }_p^{(m+1)}\}}_{p=1}^P$; 6: If $\parallel {\mathbf{\Gamma }}_{\theta }-{\mathbf{\Gamma }}_r{\parallel }_F^2>\kappa$, where ${\parallel ·\parallel }_F$ denotes matrix 2-norm, back to step 3; Otherwise, stop and output $\mathbf{\Gamma }={\mathbf{\Gamma }}_{\theta }$, $(r_p,{\theta }_p)=(r_p^{(m+1)},{\theta }_p^{(m+1)})$, $p=1,2,\dots ,P$.

4. Numerical Results

In this section, we evaluate the performance of the proposed algorithm via numerical simulations. We suppose that the non-ULA composes of M = 12 isotropous elements randomly placed in the array aperture $D=100$ m with the working wavelength $\lambda =0.15$ m and $N=1000$. Additionally, we model the array gain errors and phase errors as random variables obeying uniform distribution, which are generated by [22]:

$${\rho }_m=1+\sqrt{12}{\sigma }_p{\zeta }_m,$$

$${\phi }_m=\sqrt{12}{\sigma }_{\phi }{\eta }_m,$$

where ${\zeta }_m$ and ${\eta }_m$ are independent and identically distributed random variables distributed uniformly over $[-0.5,0.5]$, ${\sigma }_p$ and ${\sigma }_{\phi }$ are the standard deviations of ${\rho }_m$ and ${\phi }_m$, respectively. Finally, the exit condition for Algorithms 1–3 is $\kappa ={10}^{-4}$.

4.1. The Joint Estimations of the Array Gain and Phase Errors and Source Locations

In this subsection, we focus on jointly estimating the array gain and phase errors and source locations using the imprecise location knowledge of calibration sources. Specifically, without loss of generality, we consider a scenario involving three calibration sources (note that we can exploit a calibration source or two calibration sources) located at $(1400\mathrm{m},{10}^{\circ }),(1500\mathrm{m},22.5^{\circ }),(1600\mathrm{m},-{15}^{\circ })$, respectively assuming that all signal-to noise ratio (SNRs) are 20 dB. In particular, we suppose the imprecise distances (this is reasonable due to the measurement error) of the three sources are $1385\mathrm{m},1485\mathrm{m},1615\mathrm{m}$, respectively, which are chosen as initializations of Algorithm 3. Additionally, we randomly generate $L_1=50$ experiments for gain and phase errors under ${\sigma }_p=0.1$ and ${\sigma }_{\phi }=\pi /18$. For each experiment corresponding to a set of fixed gain and phase errors, we conduct $L_2=50$ Monte Carlo trials for eliminating the impact on the noise. To this end, the RMSE (Root Mean Square Error) of the gain and phase errors are defined as, respectively,

$${\mathrm{RMSE}}_1=\sqrt{\frac{{\displaystyle \sum _{l=1}^{L_2}}{\displaystyle \sum _{m=1}^M}{\left|{\rho }_m-{\rho }_m^{\left(l\right)}\right|}^2}{M{L}_2}},$$

(20)

$${\mathrm{RMSE}}_2\left(\mathrm{in}\mathrm{deg}\right)=\sqrt{\frac{{\displaystyle \sum _{l=1}^{L_2}}{\displaystyle \sum _{m=1}^M}{\left|{\phi }_m-{\phi }_m^{\left(l\right)}\right|}^2}{M{L}_2}},$$

(21)

where ${\rho }_m^{\left(l\right)}$ and ${\phi }_m^{\left(l\right)}$, respectively, are estimation values of the gain and phase errors of m-th sensor for l-th Monte Carlo trial.

Figure 2 depicts the average estimation values of the range and the azimuth of calibration source versus different experiments. Note that the obtained results are average over 50 Monte Carlo trials. We observe that all range estimation values of three calibration sources hover around their true ranges in Figure 2a, respectively. In particular, the maximum estimation error is about 9 m. This is reasonable since the obtained solution is suboptimal in Algorithm 3. Interestingly, in Figure 2b, all azimuth estimation values overlap perfectly with the theoretical values. These performance behaviors indicate that the proposed algorithm can accurately estimate the azimuth of calibration source but sharing a slight error in range.

Figure 3 shows RMSEs of the gain and phase errors versus different experiments. Results reveal that the different gain and phase errors would result in different estimation errors due to the existence of range estimation errors. In particular, it can be observed that the mean RMSEs of the gain and phase errors are about 0.015 and ${3.5}^{\circ }$, respectively. Finally, it is worth pointing out that the proposed technique shares the capability of the robustness to estimate the different gain and phase errors in correspondences of the analyzed parameters.

Next, we randomly select a set of values of the gain and phase errors among 50 experiments. In Table 1, we report the theoretical and average estimated values of the gain and phase errors for 12 sensors, where $\rho$ and $\widehat{\rho }$ are respectively theoretical and estimated values of the gain error, and $\phi$ and $\widehat{\phi }$ are, respectively, theoretical and estimated values of the phase error. Note that the estimated values for each sensor are average results of 50 Monte Carlo trials. Interestingly, it can be seen that the estimated values are close to the true values showing that the proposed technique can well estimate the gain and phase errors.

Next, we analyze the MVDR spectrum given by

$$\begin{array}{c}P(r,\theta )=\frac{1}{{\mathit{a}}^H(\theta ,r){\mathbf{\Gamma }}^H{\mathit{R}}_x^{-1}\mathbf{\Gamma }\mathit{a}(\theta ,r)}.\hfill \end{array}$$

(22)

In Figure 4, we plot the two-dimensional MVDR spatial spectrums for three calibration sources considering two cases of before correction (Figure 4a) and after correction (Figure 4b). It can be observed that the high sidelobe level emerges around the locations of true targets before calibration due to the nominal steering vector imperfectly matching the real one when existing the gain and phase errors. After calibration, as expected, the sidelobe levels significantly decrease and three notable peaks can be observed.

In Figure 5, we show the estimation values of range and azimuth of calibration source for 50 Monte Carlo trials. Again, we observe that the obtained ranges show a slight oscillatory in comparison with the true ranges, whereas the estimated azimuths perfectly overlap with true values. These performance behaviors reveal that the proposed technique can better estimate the locations of calibration sources.

4.2. Array Gain and Phase Error Compensation for Near-Field Source Localization

In this subsection, we focus on estimating the locations of near-field radiating sources exploiting the array gain and phase error matrix $\mathbf{\Gamma }$ obtained by Algorithm 3. Specifically, we consider a scenario where four near-fieldsources are located at $(1250\mathrm{m},-{35}^{\circ }),(1450\mathrm{m},{15.5}^{\circ }),(1650\mathrm{m},{20.1}^{\circ }),(1840\mathrm{m},-5^{\circ })$, respectively. In particular, for the following simulations, we consider the theoretical and estimated values of the gain and phase errors in Table 1. In addition, we search the range and azimuth of interest with steps $1\mathrm{m}$ and ${0.1}^{\circ }$, respectively.

In Figure 6, we plot the two-dimensional MVDR spatial spectrums for four sources with the same SNR = 20 dB considering two cases of before compensation and after compensation. As expected, in Figure 6a, many high sidelobe levels emerge around the true targets locations due to the effect of array gain and phase errors. However, after exploiting $\mathbf{\Gamma }$ to compensate the nominal steering vector, in Figure 6b, we can observe that the spectrum peaks can be readily found, even though the estimation may be precise. Figure 7a,b depict one-dimensional MVDR spatial spectrums at $\theta =-{35}^{\circ }$ and $r=1251$ m for a source with the estimated location $(1251\mathrm{m},-{35}^{\circ })$, respectively, considering two cases of before compensation and after compensation. The lower sidelobe levels are obtained compared with before compensation. For example, in Figure 7a, the peak sidelobe level (PSL) is about −13 dB before compensation, whereas it becomes about −27 dB after compensation.

Next, we assess the estimation results of four sources $(1250\mathrm{m},-{35}^{\circ }),(1450\mathrm{m},15.5^{\circ })$,$(1650\mathrm{m},20.1^{\circ })$,

$(1840\mathrm{m},-5^{\circ })$ with SNRs= 10 dB, 25 dB, 30 dB, 5 dB, respectively. In Figure 8, we provide the two-dimensional MVDR spatial spectrums for two cases of before compensation (Figure 8a) and after compensation (Figure 8b). We can observe that the peak of the weak target location (i.e., $(1840\mathrm{m},-5^{\circ })$) drowns by high

sidelobe levels. In particular, the four peaks are more easy to obtain after compensation compared with Figure 8a. Furthermore, we show the one-dimensional MVDR spatial spectrums at $-{35}^{\circ }$ and 1251 m in Figure 9a,b, respectively, considering two cases of before compensation and after compensation. Particularly, the low sidelobe levels can be observed after compensation while showing higher sidelobe levels in comparison with Figure 7b since the low power. Finally, it is worth pointing out that these performance behaviors reflect that error parameters estimated by the devised algorithm can compensate manifold very well and significantly decrease the sidelobe level of MVDR spectrum.

In

Table 2. Estimated values of the locations for four sources.

Target	1	2	3	4
Range (m)	1651	1451	1251	1841
Azimuth (deg)	20.1	15.5	−35	−5

, we summarize the estimated values of the locations for four sources. Results exhibit that all estimated ranges share the error values of 1 m with the true values and there is no error for the obtained azimuths. Consequently, these performance behaviors again show the effectiveness of the proposed algorithm.

5. Conclusions

In this paper, we have addressed the near-field source location problem for a non-ULA embedded in gain and phase errors. We have presented a sequential optimization calibration algorithm to simultaneously estimate error parameters and locations of calibration sources under minimum variance estimation criterion exploiting some imprecise information. The source locations and the gain and phase errors are obtained iteratively with known ranges or azimuths at each iteration. Numerical results have shown that the array gain and phase errors obtained by the devised algorithm can compensate manifold very well and reduce MVDR spectrum sidelobe levels significantly. It is worth pointing out that the calibration method framework can be extended to other DOA algorithms, i.e., MUSIC, for arbitrary array geometries. Additionally, future work would be possible to consider the study of the convergence rate of the proposed sequential optimization algorithm. Last but not least, since the considered problem represents a case of optimization under uncertainty, it would be interesting to evaluate the integration of the proposed method with robust optimization and stochastic programming techniques [35,36,37].