Robust Localization for Near- and Far-Field Signals with an Unknown Number of Sources

Abstract

Source location is a constant issue of importance of both theoretical study and practical engineering. Many pioneers have put out the corresponding solutions for near- or far-field signals, and preferred contributions are suggested. To our best knowledge, there are currently few focused approaches to the complicated situation where both near- and far-field signals exist with an unknown number of sources. Additionally, the robustness of the method must be taken into account when the additive background noise does not follow Gaussian or super-Gaussian distribution. To solve these problems, a novel method based on phased fractional lower-order moment (PFLOM) is proposed to simultaneously better preserve the signal and suppress the noise. Secondly, the whole procedure of the method containing direction of arrival (DOA) estimation, range estimation, separation of near-and far-field sources, and crucial parameter settings are studied in detail. Finally, comprehensive Monte Carlo experiments are carried out in the simulation to demonstrate the superiority of the proposed method compared to the existing competitive methods. Due to the novel method’s effectiveness with an unknown number of sources and robustness against various noises, it is believed that it could be fully utilized in more fields.

1. Introduction

Source location is an essential and crucial area of study in modern science, with a wide range of applications including sound acoustics [1,2,3,4,5,6], monitoring system [7,8], mobile communication [9,10,11,12], satellite communication [13], and so on. Based on their ranges, the sources are classified into two categories: near-field sources and far-field sources. Scholars have thoroughly investigated feasible solutions to the source location for near- or far-field ones [14,15,16,17,18,19]. In some complicated cases, the coexistence of near- and far-field sources leads to the failure of the above-mentioned solutions.

To address this tough problem, a large number of methods are provided one after another [20,21,22,23,24,25]. Among these methods, the two technologies of subspace decomposition and sparse reconstruction are introduced, and satisfactory algorithm performance is achieved [26,27,28,29]. Nonetheless, the two technologies have certain limitations. On the one hand, subspace decomposition requires the existence of the second or high-order statistics of signals with various noises [30,31]. On the other hand, sparse reconstruction is sensitive to the influence of noises [32,33,34]. It means that these targeted methods are commonly assumed and applied in the additive noise obeying Gaussian or super-Gaussian distribution. When the noise does not meet the above conditions, these methods will show significant performance degradation or even unavailability [35,36].

In addition, the estimation of the number of sources is another important problem that needs to be explored in depth. Usually, the number of łsources closely relating to the sparsity is given as known prior knowledge in research, which greatly reduces the difficulty of sparse reconstruction [37,38]. In our design, the estimation of the number is divided into two steps, total number estimation and then separation of near and far-field sources. The high total accuracy is based on the high accuracy of each step over the whole procedure [39].

To our best knowledge, there has not yet been any research focusing on source location for both near and far-field signals with an unknown number of sources in the presence of non-Gaussian noise. Therefore, more reliable and robust methodologies need to be developed. In this paper, we first employ alpha-stable distribution instead of traditional Gaussian distribution to model the additive background noise. Alpha-stable distribution equipped with four parameters has a strong generalization and can generate a wide range of distributions, including Gaussian distribution, Cauchy distribution, Levy distribution, and so on. Secondly, we provide an appropriate method to modify the covariance and suppress the noise based on PFLOM [33]. Both amplitude information and phase information are seriously analyzed to obtain high computational accuracy. Finally, we step by step derive the processes of the novel method, containing DOA estimation, range estimation, and separation of near and far-field sources.

The remainder of this paper is organized as follows: In Section 2, the signal model and noise model are briefly introduced. In Section 3, the entire procedure and crucial parameter setting in the proposed method are theoretically studied. In Section 4, Monte Carlo experiments are carried out to demonstrate the algorithm performance of the proposed method compared to the existing competitive methods. In Section 5, some remarkable conclusions are finally drawn.

2. Related Work

2.1. Signal Model

The signal model is given by the following scenario. On the sending side, there are K narrow band sources containing $K_1$ near-field ones and $K-K_1$ far-field ones. On the receiving side, there are $2M+1$ sensors forming a symmetric uniform linear array. The phase reference point is set in the center. Accordingly, the receiving signal can be expressed as

$$x(m,t)=\sum _{k=1}^{K}s(k,t){\mathrm{e}}^{\mathrm{j}\tau (m,k)}+n(m,t),-M\le m\le M$$

(1)

where $s(k,t)$ denotes the kth source, $n(m,t)$ denotes the additive noise, and $\tau (m,k)$ denotes the phase shift caused by the kth source’s propagation delay between the center and mth sensors. It is worth noting that near- and far-field sources share the same carrier frequency (wavelength). If not, the sources with different carrier frequencies and Baud rates can be separated by cyclostationary signal processing. In this study, we focus on the hard issue of sources with the same carrier frequency. Furthermore, $\tau (m,k)$ can be expanded by

$$\tau (m,k)=-\frac{2\pi dsin{\theta }_k}{\lambda }m+\frac{\pi d^2{cos}^2{\theta }_k}{\lambda r_k}{m}^2,$$

(2)

where d denotes the distance between adjacent sensors, $\lambda$ denotes the wavelength of both near and far-field sources. ${\theta }_k$ and $r_k$ denote the kth DOA and range of source, respectively. For the far-field sources (i.e., $r_k$ is very large), the second term in (2) is approximately equal to zero. From (1) and (2), the matrix expression of the signal model is given by

$$\mathbf{X}=A_1{\mathbf{S}}_1+A_2{\mathbf{S}}_2+\mathbf{N}$$

(3)

$$A_1\in {\mathbb{C}}^{(2M+1)\times K_1},A_2\in {\mathbb{C}}^{(2M+1)\times (K-K_1)},{\mathbf{S}}_1\in {\mathbb{R}}^{K_1\times L},{\mathbf{S}}_2\in {\mathbb{R}}^{(K-K_1)\times L}$$

(4)

where ${\mathbf{S}}_1$ and $A_1$ denote near-field sources and the corresponding steering matrix. As such, ${\mathbf{S}}_2$ and $A_2$ denote far-field sources and the corresponding steering matrix. ${\mathbf{S}}_1$, ${\mathbf{S}}_2$ are represented by $s(k,t)$, and $A_1$, $A_2$ are calculated from $\tau (m,k)$, and L denotes the signal length of sources. $\mathbf{N}$ is the matrix form of $n(m,t)$, and $\mathbf{X}$ is the matrix form of $x(m,t)$.

2.2. Noise Model

In traditional scientific research and engineering practice, additional background noise is commonly modeled based on Gaussian distribution derived from the central limit theorem (CLT), so-called Gaussian noise [40,41]. Nonetheless, Gaussian noise is an ideal model due to the fixed characteristic exponent, and it is not suitable to describe noise showing intensive impulsiveness. To overcome the weakness, alpha-stable distribution derived from the generalized central limit theorem (GCLT) is introduced to model the additional background noise, so-called alpha-stable noise [42,43,44]. As alpha-stable distribution does not have an explicit formula for the general probability density function, its characteristic function is usually used instead. The expression is given by

$${\phi }_X\left(t\right)=\mathbf{E}\left[{\mathrm{e}}^{\mathrm{j}tX}\right]=exp\left\{\mathrm{j}\delta t-{\gamma }^{\alpha }{\left|t\right|}^{\alpha }\left[1-\mathrm{j}\beta \mathrm{sgn}\left(t\right)\varphi (t,\alpha )\right]\right\},$$

(5)

where

$$\varphi (t,\alpha )=\left\{\begin{array}{cccc}\hfill & tan\left(\frac{\pi \alpha }{2}\right)\hfill & \hfill & \alpha \ne 1\hfill \\ \hfill & -\frac{2}{\pi }ln\left|t\right|\hfill & \hfill & \alpha =1.\hfill \end{array}\right.$$

(6)

${\phi }_X\left(t\right)$ denotes the characteristic function of a random variable X obeying alpha-stable distribution. $\alpha$, $\beta$, $\gamma$, and $\delta$ and four important parameters to model alpha-stable distribution $\mathcal{S}(\alpha ,\beta ,\gamma ,\delta )$. Specifically, $\alpha \in (0,2]$ denotes the characteristic exponent describing the two tails of the distribution. The larger $\alpha$ is, the heavier the tails are. $\beta \in [-1,1]$ denotes the skewness describing the slope of the distribution. When $\beta =0$, the distribution is completely symmetric about $x=\delta$ and is called symmetric alpha-stable (S$\alpha$S) distribution. $\gamma \in [0,\infty )$, and $\delta \in (-\infty ,\infty )$ denote skewness, scale, and location, respectively. $\gamma$ and $\delta$ are very much like $\sigma$ and $\mu$ used in Gaussian distribution $\mathcal{N}(\mu ,{\sigma }^2)$. Alpha-stable distribution has strong generalization, and Gaussian distribution as a sub-case of it has the following correspondence:

$$\mathcal{N}(\mu ,{\sigma }^2)\leftrightharpoons \mathcal{S}(2,\beta ,\gamma ,\delta ).$$

(7)

When $\alpha =2$, we have $\varphi (t,\alpha )=0$, and it means the value of $\beta$ has nothing to do with ${\phi }_X\left(t\right)$. In (1), the additive background noise $n(m,t)$ can be further expressed as

$$n(m,t)=\sum _{k=1}^{K}n(k,t),$$

(8)

where $n(k,t)$ denotes the additive alpha-stable noise of each source $s(k,t)$. If $n(k,t)$ with different values of k are independent of each other, we can draw an important conclusion that $n(m,t)$ is also alpha-stable noise. Particularly, when all $n(k,t)$ have the same characteristic exponent, the four parameters in $n(m,t)$ can be calculated according to the following mathematical property [45]:

Property 1. 

If$X_1\sim \mathcal{S}(\alpha ,{\beta }_1,{\gamma }_1,{\delta }_1)$and$X_2\sim \mathcal{S}(\alpha ,{\beta }_2,{\gamma }_2,{\delta }_2)$are independent, then

$$X=X_1+X_2\sim \mathcal{S}(\alpha ,\beta ,\gamma ,\delta )$$

(9)

$$\beta =\frac{{\beta }_1{\gamma }_1^{\alpha }+{\beta }_2{\gamma }_2^{\alpha }}{{\gamma }_1^{\alpha }+{\gamma }_1^{\alpha }},\gamma =\sqrt[\alpha ]{{\gamma }_1^{\alpha }+{\gamma }_2^{\alpha }},\delta ={\delta }_1+{\delta }_2$$

(10)

3. Proposed Method

3.1. DOA Estimation

The covariance of the $m_1$th and $m_2$th sensors is given by

$$C_x(m_1,m_2)=\mathbf{E}\left[x(m_1,t)x^{*}(m_2,t)\right],$$

(11)

where $\mathbf{E}$ denotes the time averaging. It is different from the same symbol $\mathbf{E}$ denoting the statistical averaging used in (5). From the perspective of the matrix given in (3), the operation is computing the inner product between the $m_1$th row of $\mathbf{X}$ and the $m_2$th row of ${\mathbf{X}}^{*}$. From (1) and (11), it can be seen that, when the additive background noise obeys alpha-stable distribution, $C_x(m_1,m_2)$ does not exist in theory due to the second order statistics. To solve this problem, we employ PFLOM to replace the traditional covariance to guarantee the existence. The PFLOM of the $m_1$th and $m_2$th sensors is defined by

$$C_x^{\langle p\rangle }(m_1,m_2)=\mathbf{E}\left[x^{\langle p\rangle }(m_1,t)x^{-\langle p\rangle }(m_2,t)\right]$$

(12)

where

$$\left\{\begin{array}{l}{x}^{\langle p\rangle }={\left|x\right|}^{p-1}x\hfill \\ x^{-\langle p\rangle }={\left[x^{\langle p\rangle }\right]}^{*}={\left|x\right|}^{p-1}{x}^{*}\hfill \end{array}\right.$$

(13)

In (13), p denotes the fractional lower order and satisfies the domain of $0\le p<\alpha /2\le 1$. Furthermore, $C_x^{\langle p\rangle }(m_1,m_2)$ is expanded as follows:

$$\begin{array}{cc}\hfill & \mathbf{E}\left[x^{\langle p\rangle }(m_1,t)x^{-\langle p\rangle }(m_2,t)\right]\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill & =\mathbf{E}\left\{|x(m_1,t){|}^{p-1}{\left|x(m_2,t)\right|}^{p-1}x(m_1,t)x^{*}(m_2,t)\right\}\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill & =\mathbf{E}\left\{A_0\left[{\displaystyle \sum _{k=1}^K}s(k,t){\mathrm{e}}^{\mathrm{j}\tau (m_1,k)}+n(m_1,t)\right]\left[{\displaystyle \sum _{k=1}^K}{s}^{*}(k,t){\mathrm{e}}^{-\mathrm{j}\tau (m_2,k)}+n^{*}(m_2,t)\right]\right\}\hfill \\ \hfill & =\mathbf{E}\left\{A_0{\displaystyle \sum _{k_1=1}^K}{\displaystyle \sum _{k_2=1}^K}s(k_1,t)s^{*}(k_2,t){\mathrm{e}}^{\mathrm{j}[\tau (m_1,k_1)-\tau (m_2,k_2)]}\right\}+\mathbf{E}\left[A_{0}n(m_1,t)n^{*}(m_2,t)\right]\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill & +\mathbf{E}\left\{A_{0}n(m_1,t){\displaystyle \sum _{k=1}^K}{s}^{*}(k,t){\mathrm{e}}^{-\mathrm{j}\tau (m_1,k)}\right\}+\mathbf{E}\left\{A_0{n}^{*}(m_2,t){\displaystyle \sum _{k=1}^K}s(k,t){\mathrm{e}}^{\mathrm{j}\tau (m_2,k)}\right\}\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill & =A\mathbf{E}\left\{{\displaystyle \sum _{k=1}^K}|s(k,t){|}^2{\mathrm{e}}^{\mathrm{j}[\tau (m_1,k)-\tau (m_2,k)]}\right\}=A{\displaystyle \sum _{k=1}^K}{P}_k{\mathrm{e}}^{\mathrm{j}[\tau (m_1,k)-\tau (m_2,k)]}\hfill \end{array}$$

(18)

In (14), $A_0$ denotes the decay term of the amplitude information, and has nothing to do with the phase information. In (18), $P_k$ denotes the power of $s(k,t)$, and, from (16) to (18), is based on the condition that all sources and noises are independent from each other. When $m_1=-m_2=-m$, we can further simplify to obtain the definition of $T_x\left(m\right)$:

$$T_x\left(m\right)\triangleq C_x^{\langle p\rangle }(-m,m)=A\sum _{k=1}^K{P}_{k}exp\left(\frac{\mathrm{j}4\pi dsin{\theta }_k}{\lambda }m\right)=A\sum _{k=1}^K\varphi (m,k)P_k$$

(19)

Accordingly, $\mathbf{T}$ is the vector form of $T_x\left(m\right)$ that can be obtained by traversing all the values of m taken:

$$\left[\begin{array}{c}{T}_x(-M)\\ T_x(-M+1)\\ ⋮\\ T_x\left(M\right)\end{array}\right]=A\left[\begin{array}{cccc}\varphi (-M,1)& \varphi (-M,2)& \cdots & \varphi (-M,K)\\ \varphi (-M+1,1)& \varphi (-M+1,2)& \cdots & \varphi (-M+1,K)\\ ⋮& ⋮& \ddots & ⋮\\ \varphi (M,1)& \varphi (M,2)& \cdots & \varphi (M,K)\end{array}\right]\left[\begin{array}{c}{P}_1\\ P_2\\ ⋮\\ P_K\end{array}\right]$$

(20)

$$\mathbf{T}={\mathbf{\Phi }}_0{\mathbf{P}}_0$$

(21)

Regretfully, $P_k$ can not be directly obtained because of the unknown number of sources K. When the source and the DOA are one-to-one, we can extend the matrix by traversing the angle from $-{90}^{\circ }$ to ${90}^{\circ }$ and selecting $1^{\circ }$ as the interval, leading to

$${\theta }_k\in {\theta }_l,{\theta }_l=l-91,1\le l\le 181.$$

(22)

Correspondingly, the expansion of $\mathbf{T}$ is revised as follows:

$$\left[\begin{array}{c}{T}_x(-M)\\ T_x(-M+1)\\ ⋮\\ T_x\left(M\right)\end{array}\right]=\left[\begin{array}{cccc}\varphi (-M,1)& \varphi (-M,2)& \cdots & \varphi (-M,181)\\ \varphi (-M+1,1)& \varphi (-M+1,2)& \cdots & \varphi (-M+1,181)\\ ⋮& ⋮& \ddots & ⋮\\ \varphi (M,1)& \varphi (M,2)& \cdots & \varphi (M,181)\end{array}\right]\left[\begin{array}{c}{P}_1\\ P_2\\ ⋮\\ P_{181}\end{array}\right]$$

(23)

$$\mathbf{T}={\mathbf{\Phi }}_1{\mathbf{P}}_1$$

(24)

It is easy to see that (24) is an under-determined equation, and the nonzero term $P_k$ can be estimated by sparse reconstruction.

The approaches for sparse reconstruction are mainly divided into two categories: greedy algorithm and convex optimization. Generally speaking, the greedy algorithm has less computational complexity, while convex optimization has higher computational accuracy. Typical greedy algorithms include: OMP $\mathcal{O}\left(M_0{N}_0{K}_0\right)$[46], SP $\mathcal{O}\left(M_0{N}_0{L}_0\right)$[47], and IHT $\mathcal{O}\left(M_0{N}_0{L}_0\right)$[48]. Typical convex optimizations include: BPDN $\mathcal{O}\left(M_0^2{N}_0^3\right)$[38] and Reweighted ${\ell }_1$ minimization $\mathcal{O}\left(M_0^2{N}_0^3{L}_0\right)$[49]. $M_0$, $N_0$, $K_0$, and $L_0$ denote the dimension of the observation vector, the dimension of the original vector, the sparsity of the original vector, and the number of iterations, respectively. Of course, there are also some examples of exceptions. Some greedy algorithms with tree search and backtracking technology have computational accuracy equal to or even higher than convex optimization at the cost of computational complexity [50,51,52]. These greedy algorithms are suitable for use in the case of large sparsity. In the specific scenario of DOA estimation, we have

$$M_0=2M+1,N_0=181,K_0=K.$$

(25)

At the same time, the sparse reconstruction in the complex domain will increase the computational complexity to a certain extent. After full consideration of all aspects, we can accept the computational complexity of convex optimization and think that it is a suitable choice for the whole procedure. The estimation of ${\mathbf{P}}_1$ can be expressed as

$${\widehat{\mathbf{P}}}_1=arg\underset{{\mathbf{P}}_1}{min}\mathcal{L}\left({\mathbf{P}}_1\right)=arg\underset{{\mathbf{P}}_1}{min}\parallel \mathbf{T}-{\mathbf{\Phi }}_1{\mathbf{P}}_1{\parallel }_2+{\lambda }_1{\parallel {\mathbf{P}}_1\parallel }_1,$$

(26)

where $\mathcal{L}(·)$ denotes loss function. If the kth term in ${\widehat{\mathbf{P}}}_1$ is a peak, there is a near or far-field source in this direction such that

$$\left\{{\widehat{P}}_k\mathrm{is}\mathrm{a}\mathrm{peak}\mathrm{in}{\widehat{\mathbf{P}}}_1\right\}\Rightarrow {\widehat{\theta }}_k=k-91.$$

(27)

3.2. Range Estimation

According to the DOAs achieved in the first process, we secondly estimate the range of both near and far-field sources by a subspace-based method. Substituting the estimated DOA ${\theta }_k$ into (2), we have

$$\begin{array}{cc}\hfill \tau (m,{\widehat{\theta }}_k,r_k)& =-\frac{2\pi dsin{\widehat{\theta }}_k}{\lambda }m+\frac{\pi d^2{cos}^2{\widehat{\theta }}_k}{\lambda r_k}{m}^2\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill {\mathrm{e}}^{\mathrm{j}\tau (m,{\widehat{\theta }}_k,r_k)}& =exp\left(-\frac{\mathrm{j}2\pi dsin{\widehat{\theta }}_k}{\lambda }m\right)exp\left(\frac{\mathrm{j}\pi d^2{cos}^2{\widehat{\theta }}_k}{\lambda r_k}{m}^2\right).\hfill \end{array}$$

(29)

Furthermore, the matrix form of ${\mathrm{e}}^{\tau (m,{\widehat{\theta }}_k,r_k)}$ with different values of m is given by

$$\mathbf{H}={[{\mathrm{e}}^{\mathrm{j}\tau (-M,{\widehat{\theta }}_k,r_k)},{\mathrm{e}}^{\mathrm{j}\tau (-M+1,{\widehat{\theta }}_k,r_k)},\dots ,{\mathrm{e}}^{\mathrm{j}\tau (M-1,{\widehat{\theta }}_k,r_k)},{\mathrm{e}}^{\mathrm{j}\tau (M,{\widehat{\theta }}_k,r_k)}]}^{\mathsf{T}}.$$

(30)

Then, the range $r_k$ can be estimated by searching the peak of the following spectrum:

$$Q({\widehat{\theta }}_k,r_k)=\frac{1}{{\mathbf{H}}^{\mathsf{H}}\mathbf{U}{\mathbf{U}}^{\mathsf{H}}\mathbf{H}}$$

(31)

where $\mathbf{U}$ denotes the noise subspace of the eigenvectors. These eigenvectors are corresponding to the $2M+1-K$ eigenvalues with the smallest absolute value. In order to reduce the computational complexity, we employ an effective Root-MUSIC like method for reference. According to (28), $\mathbf{H}$ can be further decomposed into the product of two exponential terms such that

$$\left[\begin{array}{c}{\mathrm{e}}^{\mathrm{j}\tau (-M,{\widehat{\theta }}_k,r_k)}\\ ⋮\\ {\mathrm{e}}^{\mathrm{j}\tau (M,{\widehat{\theta }}_k,r_k)}\end{array}\right]=\left[\begin{array}{ccc}exp(\frac{\mathrm{j}2\pi dsin{\widehat{\theta }}_k}{\lambda }M)& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & exp(-\frac{\mathrm{j}2\pi dsin{\widehat{\theta }}_k}{\lambda }M)\end{array}\right]\left[\begin{array}{c}g{({\widehat{\theta }}_k,r_k)}^{{(-M)}^2}\\ ⋮\\ g{({\widehat{\theta }}_k,r_k)}^{M^2}\end{array}\right]$$

(32)

$$\mathbf{H}=\Lambda \mathbf{G},$$

(33)

and $g{({\widehat{\theta }}_k,r_k)}^{m^2}$ is given by

$$g{({\widehat{\theta }}_k,r_k)}^{m^2}=exp\left(\frac{\mathrm{j}\pi d^2{cos}^2{\widehat{\theta }}_k}{\lambda r_k}{m}^2\right).$$

(34)

We define a new polynomial about $g({\widehat{\theta }}_k,r_k)$ and make it equal to zero to obtain $g({\widehat{\theta }}_k,r_k)$.

$$f[{\widehat{\theta }}_k,g({\widehat{\theta }}_k,r_k)]\triangleq g{({\widehat{\theta }}_k,r_k)}^{M^2}{\mathbf{H}}^{\mathsf{H}}\mathbf{U}{\mathbf{U}}^{\mathsf{H}}\mathbf{H}=g{({\widehat{\theta }}_k,r_k)}^{M^2}{\mathbf{G}}^{\mathsf{H}}\Lambda \mathbf{U}{\mathbf{U}}^{\mathsf{H}}\Lambda \mathbf{G}=0$$

(35)

Among a large number of roots, the root that is closest to the unit circle is chosen since $|g({\widehat{\theta }}_k,r_k)|=1$. Finally, based on the result of $\widehat{g}({\widehat{\theta }}_k,r_k)$, we have

$${\widehat{r}}_k=\frac{\pi d^2{cos}^2{\widehat{\theta }}_k}{\lambda \mathbf{angle}[\widehat{g}({\widehat{\theta }}_k,r_k)]}.$$

(36)

3.3. Separation of Near and Far-Field Sources

According to the ranges achieved in the second process, we finally separate the near and far-field sources by the Fresnel region [16]:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \left\{\begin{array}{cc}{\widehat{r}}_k\in [0.62\sqrt{D^3/\lambda },2D^2/\lambda ],\hfill & \mathrm{near}-\mathrm{field}\mathrm{source}\hfill \\ {\widehat{r}}_k\notin [0.62\sqrt{D^3/\lambda },2D^2/\lambda ],\hfill & \mathrm{far}-\mathrm{field}\mathrm{source}\hfill \end{array}\right.\end{array}\end{array}$$

(37)

where $D=2Md$ denotes the array aperture. If a source is determined to be far-field, the estimated range previously is meaningless. This is because there is no $r_k$ in $\tau (m,k)$ of a far-field source according to (2). In other words, the estimated range is only effective for near-field sources.

3.4. Crucial Parameter Setting

3.4.1. Regularization Coefficient ${\lambda }_1$

Calculating the partial derivative of ${\mathbf{P}}_1$ from (26), we have

$$\frac{\partial \mathcal{L}\left({\mathbf{P}}_1\right)}{\partial {\mathbf{P}}_1}=\frac{{\mathbf{\Phi }}_1^{\mathsf{H}}({\mathbf{\Phi }}_1{\mathbf{P}}_1-\mathbf{T})}{\parallel \mathbf{T}-{\mathbf{\Phi }}_1{\mathbf{P}}_1{\parallel }_2}+{\lambda }_1\mathbf{C},$$

(38)

where $\mathbf{C}=\mathrm{ones}(181,1)={[1,1,\cdots ,1]}^{\mathsf{T}}$ because $P_k\ge 0$. Let the partial derivative be zero, and the analytic expression of $\lambda$ is given by

$${\lambda }_1=\frac{{\mathbf{C}}^{\mathsf{T}}{\mathbf{\Phi }}_1^{\mathsf{H}}(\mathbf{T}-{\mathbf{\Phi }}_1{\mathbf{P}}_1)}{181\parallel \mathbf{T}-{\mathbf{\Phi }}_1{\mathbf{P}}_1{\parallel }_2}.$$

(39)

In fact, the solution of ${\mathbf{P}}_1$ and ${\lambda }_1$ is a cross-iterative process. Most researchers use the toolboxes and a large number of training sets to estimate ${\lambda }_1$ in various applications. In addition, if the source signals are non-stationary, $\mathbf{T}$ is a high dimensional random variable, which makes it difficult to find the optimal solution of ${\lambda }_1$ under different conditions.

In our scenario, we carry out some experiments to observe the impact of ${\lambda }_1$ on DOA estimation. From Figure 1, ${\lambda }_1=1.5$ is a proper value for the sparse reconstruction as the four inside peaks with high amplitudes and peak-to-average powers (PAPRs).

3.4.2. Evaluation Standard

Regarding the traditional DOA estimation evaluation standard, we consider the following equation using the absolute deviation and the negative exponential function to evaluate the DOA estimation:

$$f_1=\frac{1}{K}\sum _{k=1}^{K}exp(-|{\widehat{\theta }}_k-{\theta }_k|)$$

(40)

when the number of peaks is smaller than K (e.g., Figure 1c), and the DOA estimation is incorrect and returns $f_1=0$. When the number of peaks is larger than K (e.g., Figure 1a,d), the result is calculated from the K highest peaks. In the same way, the range estimation is evaluated as follows:

$$f_2=\frac{1}{K_1}\sum _{k=1}^{K_1}exp(-|{\widehat{r}}_k-r_k|)$$

(41)

The range estimation must be predicated on the successful DOA estimation and successful source separation as a foundation. Specifically, if the number of peaks is smaller than K or the estimated number of near-field sources is not $K_1$, return $f_2=0$.

4. Simulation

Both near and far-field sources are quadrature phase-shift keying (QPSK) signals with identical mapping mode and power. Because the second-order statistics related to the power of alpha-stable noise do not exist, the generalized power is introduced to replace the traditional power. In order to quickly generate a large number of noises with the same generalized power for Monte Carlo experiments, we set skewness $\beta$ and location $\delta$ to zero to make alpha-stable noise without losing the generality. Consequently, the generalized signal-to-noise (GSNR) is defined by

$$\mathrm{GSNR}=10lg\left(\frac{P_x}{{\gamma }^{\alpha }}\right),\beta =0,\delta =0.$$

(42)

where ${\gamma }^{\alpha }$ denotes the generalized power, which is very similar to ${\sigma }^2$ indicating the power of $\mathcal{N}(0,\sigma )$. With the preset characteristic exponent $\alpha$ and GSNR, scale $\gamma$ can be directly calculated by

$$\gamma =\sqrt[\alpha ]{\frac{P_x}{{10}^{\mathrm{GSNR}/10}}},\beta =0,\delta =0.$$

(43)

Correntropy is another effective method to deal with alpha-stable noise with a broad consensus in research [53]. It is necessary to take this existing competitive method into consideration. Nevertheless, in this application, correntropy can not work since its result is a real number without any phase information [54,55]. Therefore, correntropy is not added to the simulation. Instead, the MS-subspace method and its modified version employing PFLOM are added [56]. As there is no sparse reconstruction step in the MS-subspace method and its modified version, the regularization coefficient ${\lambda }_1$ has nothing to do with them. The parameters used in the simulation are listed in

Table 1. Parameter setting in the simulation.

Parameter	Value
DOAs of near- and far-field sources	$[-6^{\circ },{13}^{\circ }]$, $[-{40}^{\circ },{31}^{\circ }]$
Ranges of near-field sources	$[10\lambda$, $15\lambda ]$
Baud rates of near- and far-field sources	$5\times {10}^6$ Baud, ${10}^7$ Baud
Carrier frequency and sampling frequency	20 Mhz, 100 Mhz
Snapshots	5000
GSNR	$10:5:20$ dB
Characteristic exponent, $\alpha$	$0.6:0.1:2.0$
Scale parameter, $\gamma$	calculated by GSNR and $\alpha$

.

Figure 2, Figure 3 and Figure 4 show the performance comparison of different methods according to the evaluation standards of $f_1$ and $f_2$. Each point in the figures is calculated from the statistical average of one hundred Monte Carlo experiments. From the figures, it can be seen that the fractional order p has a great impact on both $f_1$ and $f_2$, which is consistent with the principle of PFLOM. The conventional covariance-based method and MS-subspace method do very well in Gaussian noise ($\alpha =2$). However, with the decreasing of character exponent $\alpha$ from 2 to $0.6$, only the PFLOM-based method ($p=0.2$) and the modified MS-subspace method ($p=0.2$) demonstrate strong robustness, but the other methods are severely degraded due to the growing frequency and intensity of the impulses in alpha-stable noise. In comparison, the PFLOM-based method ($p=0.2$) is better than the modified MS-subspace method ($p=0.2$) in accuracy.

In Figure 2a, when $\alpha =2$, the covariance-based method and MS-subspace method perform best in line with the conclusion that the second- and fourth- order statistics perform well with Gaussian and super-Gaussian noises. When $\alpha =0.6$, the PFLOM-based method ($p=0.2$) performs best because of its algorithm characteristic and the smallest value of p. The performance of the PFLOM-based methods ($p=0.4:0.2:0.8$) is between the PFLOM-based method ($p=0.2$) and covariance-based method. Similarly, the performance of the modified MS-subspace methods ($p=0.4:0.2:0.8$) is between the modified Ms-subspace method ($p=0.2$) and Ms-subspace method. In addition, the PFLOM-based methods with different values of p are monotonic, while the modified MS-subspace methods are not.

The only difference between Figure 2a,b is that the value of ${\lambda }_1$ is taken differently. Overall, the five curves of the covariance-based method and PFLOM-based method in Figure 2a łhave trends similar to the corresponding curves in Figure 2b, i.e., the MS-subspace method and its modified version are not related to ${\lambda }_1$. The PFLOM-based methods ($p=0.2:0.2:0.8$) in Figure 2b perform a bit worse than those in Figure 2a, which means ${\lambda }_1$ is another important factor to effect the methods containing the step of convex optimization. Unlike p, ${\lambda }_1$ is not monotonic in the PFLOM-based methods, and $[1.6,1.9]$ is an acceptable definition domain. It is easy to see that the effect of p on the algorithm performance is higher than that of ${\lambda }_1$.

Comparing Figure 2a–d, we see that $f_1>f_2$ under all conditions including $\alpha$, ${\lambda }_1$, and p. This means that the accuracy of the range estimation is lower than that of the DOA estimation, which is in line with the objective facts that the range estimation must be predicated on the correct DOA estimation. In other words, the range estimation as a further process is more difficult than the DOA estimation. In Figure 2c,d, the covariance-based method and MS-subspace method perform best when $\alpha =2$ but has no robustness against the decreasing of $\alpha$. The PFLOM-based method ($p=0.2$) and the modified MS-subspace method ($p=0.2$) have the strongest robustness at the cost of slightly worse performance at $\alpha =2$. In comparison, the PFLOM-based method ($p=0.2$) is better than the modified MS-subspace method ($p=0.2$) in accuracy.

Figure 3 and Figure 4 show the performance comparisons of different methods under the conditions of $\mathrm{GSNR}=15\mathrm{dB}$ and $\mathrm{GSNR}=20\mathrm{dB}$. From the two figures, it can be seen that all methods perform better than those in Figure 2 because of the higher GSNRs, and maintain similar trends because of their own algorithm characteristics.

5. Conclusions

In this paper, a novel method for robust localization for near and far-field signals with an unknown number of sources is provided. The proposed method employs PFLOM to maintain the signals of interest and suppress the various noises, and further achieves high accuracy DOA estimation, range estimations, and separation of near- and far-field sources. From the numerical results in Monte Carlo experiments, the impacts of fractional order p and regularization coefficient ${\lambda }_1$ are clearly shown, which is of most importance to describing the special characteristics of the methodology. The proposed method with a suitable fractional order demonstrates good effectiveness and strong robustness under different conditions. It is thought that it might be used to its full potential in additional areas.