Passive Localization of 3D Near-Field Cyclostationary Sources Using Parallel Factor Analysis

Abstract

By exploiting favorable characteristics of a uniform cross-array, a passive localization algorithm of narrowband cyclostationary sources in the spherical coordinates (azimuth, elevation, and range) is proposed. Firstly, we construct a parallel factor (PARAFAC) analysis model by computing the third-order cyclic moment matrices of the properly chosen sensor outputs. Then, we analyze the uniqueness of the constructed model and obtain three-dimensional (3D) near-field parameters via trilinear alternating least squares regression (TALS). The investigated algorithm is well suitable for the localization of the near-field cyclostationary sources. In addition, it avoids the multidimensional search and pairing parameters. Results of computer simulations are carried out to confirm the satisfactory performance of the proposed method.

1. Introduction

There has been considerable interest in bearing estimation for radar, sonar, communication, and electronic surveillance [1]. Various high-resolution algorithms, such as MUSIC [2] and ESPRIT [3], have been proposed to obtain the direction-of-arrival estimation of the far-field sources. In addition, when the sources are localized at the Fresnel region [4] of the array aperture, both the azimuth and range should be estimated. Recently, a significant amount attention has been paid to this issue and several near-field sources localization algorithms [5–7] are also available. However, all these methods as mentioned above only address the 2D problem of estimating azimuth and range and rely on the assumption of the stationary sources.

In recent years, several 3D near-field sources localization methods have been developed to obtain azimuth, elevation, and range. Meraim and Hua [8] proposed a second-order based method to cope with this issue. By translating the 1D uniform linear array of near field into a virtual rectangular array of virtual far field, Challa and Shamsunder [9] developed a fourth-order cumulants based Unitary-ESPRIT method. Moreover, [10, 11] also introduced efficient localization methods for 3-D near-field non-Gaussian stationary sources. It is obviously seen that all these methods still require that the incoming signals should be stationary ones; in addition, a multidimensional search or pairing parameters is also failed to avoid.

Cyclostationarity, which is a statistical property processed by most man-made communication signals, is related to the underlying periodicity arising from cyclic frequency of based rates. Besson et al. [12] introduced an original far-field approximation [13] based method to deal with the bearings estimation of near-field cyclostationary sources. Due to the fact that the effect caused by the mismatch between the actual spherical wavefront phase vector and assumed planar wavefront vector was alleviated, the proposed technique showed a satisfactory performance for the near-field sources far away from the sensor array. The estimations for elevation and range, however, have not been well considered.

In this paper, we consider the problem of jointly estimating elevation, azimuth, and range of the near-field cyclostationary sources; what is more, a two-stage passive localization method has been proposed. In the first stage, several third-order cyclic moment matrices of a cross-array observations data are computed, and a parallel factor analysis model in the cyclic statistic domain is constructed. In the second stage, the uniqueness of the constructed model is proved; in addition, the 3-D localizing parameters for the near-field sources have also been obtained through TALS. The algorithm developed in this paper would be well suitable for near-field cyclostationary sources, and it does not require multidimensional search or pairing parameters; in addition, it also can effectively alleviates the array aperture loss.

The rest of this paper is organized as follows. Section 2 introduces the signal model of near-field localization based on cross-array. Section 3 develops a joint estimation algorithm of three parameters in near-field. Section 4 shows simulation results. Section 5 presents the conclusion of the whole paper.

2. Near-Field Signal Model Based Cross-Array

2.1. Near-Field Signal Model

We consider a near-field scenario of 𝑀 uncorrelated narrowband signals impinging on a cross-array signed with the 𝑥 and 𝑦 axes (Figure 1), which consists of 𝐿=2𝑁+3 sensors with element spacing 𝑑.

Figure 1 

Sensor-source configuration for the near-field problem.

Let the array center be the phase reference point, and the signals received by the (𝑙,0)th and the (𝑙,0)th can be, respectively, expressed as 𝑥_𝑙,0(𝑡)=^𝑀_𝑚=1𝑠_𝑚(𝑡)exp𝑗𝛾_𝑥𝑚𝑙+𝜑_𝑥𝑚𝑙²+𝑛_𝑙,0(𝑡),𝑥_0,𝑙(𝑡)=^𝑀_𝑚=1𝑠_𝑚(𝑡)exp𝑗𝛾_𝑗𝑚𝑙+𝜑_𝑗𝑚𝑙²+𝑛_0,𝑙(𝑡), (1) where 𝑛(𝑡) donates the sensor additive noise, and 𝛾_𝑥𝑚𝑑=−2𝜋𝜆sin𝛼_𝑚cos𝜃_𝑚,𝜑_𝑥𝑚𝑑=𝜋²𝜆𝑟_𝑚1−sin²𝛼_𝑚cos²𝜃_𝑚,𝛾_𝑦𝑚𝑑=−2𝜋𝜆sin𝛼_𝑚sin𝜃_𝑚,𝜑_𝑦𝑚𝑑=𝜋²𝜆𝑟_𝑚1−sin²𝛼_𝑚sin²𝜃_𝑚, (2) where 𝛼_𝑚, 𝜃_𝑚, and 𝑟_𝑚 indicate elevation, azimuth, and range of 𝑚th signal, respectively, and 𝜆 is wavelength of source signal.

The 𝑚th source signal with amplitude 𝑧_𝑚(𝑡) can be modeled as 𝑠_𝑚(𝑡)=𝑧_𝑚(𝑡)exp𝑗𝜔𝑡𝑚=1,2,…,𝑀, (3) where 𝜔 is the center frequency.

2.2. Assumption of Signal Model

The main problem addressed in this paper is to jointly estimate the sets of parameters (𝛼_𝑚,𝜃_𝑚,and𝑟_𝑚), and then the following assumptions are assumed to hold:(1)the envelope 𝑧_𝑚(𝑡) is non-Gaussian stationary random process with zero mean and nonzero skewness;(2)the sensor noise is the additive stationary one and independent from the source signals;(3)the sensor array is a uniform linear array with element spacing 𝑑≤𝜆/4; in addition, the source number 𝑀 is not more than sensor number 𝑁.

2.3. Parallel Factor Analysis

We need to introduce the following notation that will be used in the sequel.

Definition 1 (see [14]). Let 𝑥_{𝑘,𝑛,𝑝} stand for the (𝑘,𝑛,𝑝) element of a three-dimensional tensor 𝐗, if 𝑥_{𝑘,𝑛,𝑝}=^𝑀_𝑚=1𝑎_𝑘,𝑚𝑏_𝑛,𝑚𝑐_𝑝,𝑚, (4)

where 𝑎_𝑘,𝑚 denotes the (𝑘,𝑚) element of matrix 𝐀 and similarly for the others. Equation (4) indicates 𝑥_{𝑘,𝑛.𝑝} as a sum of triple products, which is variably known as the trilinear model, trilinear decomposition, triple product decomposition, canonical decomposition, or parallel factor (PARAFAC) analysis.

Definition 2 (see [14]). For a matrix 𝐁∈𝐶^𝐼×𝐽, if all 𝐼<𝐽 columns of 𝐁 are linearly independent, but there exists a collection of 𝐼+1 linearly dependent columns of 𝐁, then it has Kruskal-rank (𝑘-rank) 𝑘_𝐁=𝐼.

Theorem 3 (see [15]). Consider a three dimensional tensor 𝐗 as defined in (4), and 𝑀 represents the common dimension, if 𝑘_𝐀+𝑘_𝐁+𝑘_𝐂^𝑇≥2𝑀+2, (5) then 𝐀, 𝐁, and 𝐂 are unique up to permutation and (complex) scaling of columns.

3. PARAFAC Based 3D Near-Field Sources Localization

3.1. Compute the Third-Order Cyclic Moments Matrices

Let 𝑀^𝑎_3,𝑥(0,−𝑙−1,−𝑙) denote the third-order cyclic moment with cyclic frequency 𝛼 defined as 𝑀^𝑎_3,𝑥(0,−𝑙−1,−𝑙)=lim_𝑇→∞1𝑇^𝑇_𝑡=1𝐸𝑥_0,0(𝑡)𝑥_{−𝑙−1,0}(𝑡)𝑥^∗_−𝑙,0(𝑡)exp−𝑗𝛼𝑡, (6) where the superscript “*” denotes the complex conjugate operation; in addition, 𝐸{𝑥_0,0(𝑡)𝑥_{−𝑙−1,0}(𝑡)𝑥^∗_−𝑙,0(𝑡)} can be given by 𝐸𝑥_0,0(𝑡)𝑥_{−𝑙−1,0}(𝑡)𝑥^∗_−𝑙,0(𝑡=_𝑀𝑚=1)𝑚_3,𝑧𝑚exp𝑗−𝛾_𝑥𝑚+𝜙_𝑥𝑚exp𝑗2𝜙_𝑥𝑚𝑙exp𝑗𝜔𝑡+𝑚^𝑀_1,𝑛𝑚=1𝑚_2,𝑧𝑚exp𝑗𝛾_𝑥𝑚(−𝑙−1)+𝜙_𝑥𝑚(−𝑙−1)²×exp𝑗2𝜔𝑡+𝑚^𝑀_1,𝑛𝑚=1𝑚_2,𝑧𝑚exp−𝑗𝛾_𝑥𝑚(−𝑙)+𝜙_𝑥𝑚(−𝑙)²+𝑚^𝑀_2,𝑛𝑚=1𝑚_1,𝑧𝑚exp𝑗𝜔𝑡+𝑚^𝑀_1,𝑛𝑚=1𝑚_2,𝑧𝑚exp𝑗−𝛾_𝑥𝑚+𝜙_𝑥𝑚exp𝑗2𝜙_𝑥𝑚𝑙+𝑚^𝑀_2,𝑛𝑚=1𝑚_1,𝑧𝑚exp𝑗𝛾_𝑥𝑚(−𝑙−1)+𝜙_𝑥𝑚(−𝑙−1)²×exp𝑗𝜔𝑡+𝑚^𝑀_2,𝑛𝑚=1𝑚_1,𝑧𝑚exp−𝑗𝛾_𝑥𝑚(−𝑙)+𝜙_𝑥𝑚(−𝑙)²×exp𝑗𝜔𝑡+𝑚_3,𝑛, (7) where 𝑚_1,𝑧𝑚, 𝑚_2,𝑧𝑚, and 𝑚_3,𝑧𝑚 denote mean, variance, and the third-order moment of 𝑧_𝑚(𝑡); in addition, 𝑚_1,𝑛, 𝑚_2,𝑛, and 𝑚_3,𝑛 represent mean, variance, and the third-order cyclic moment of 𝑛(𝑡), respectively.

With assumption being considered, we further obtain 𝑀^𝑎_3,𝑥()=0,−𝑙−1,−𝑙_𝑀𝑚=1𝑚_3,𝑧𝑚exp𝑗−𝛾_𝑥𝑚+𝜙_𝑥𝑚exp𝑗2𝜙_𝑥𝑚𝑙𝛿(𝜔−𝛼)+𝑚^𝑀_1,𝑛𝑚=1𝑚_2,𝑧𝑚exp𝑗𝛾_𝑥𝑚(−𝑙−1)+𝜙_𝑥𝑚(−𝑙−1)²×𝛿(2𝜔−𝛼)+𝑚^𝑀_1,𝑛𝑚=1𝑚_2,𝑧𝑚exp−𝑗𝛾_𝑥𝑚(−𝑙)+𝜙_𝑥𝑚(−𝑙)²𝛿(𝛼)+𝑚^𝑀_1,𝑛𝑚=1𝑚_2,𝑧𝑚exp𝑗−𝛾_𝑥𝑚+𝜙_𝑥𝑚exp𝑗2𝜙_𝑥𝑚𝑙𝛿(𝛼)+𝑚_3,𝑛𝛿(𝛼), (8) where 𝛿 is impulse function. Letting 𝑎=𝜔 and substituting (8) into (6) yield 𝑀^𝑎_3,𝑥()=0,−𝑙−1,−𝑙_𝑀𝑚=1𝑚_3,𝑧𝑚exp𝑗−𝛾_𝑥𝑚+𝜑_𝑥𝑚exp𝑗2𝑙𝜑_𝑥𝑚. (9)

Based on (9), we reconstruct the spatial third-order cyclic moment matrix 𝐌^𝑎₁, in which the (𝑘,𝑞)th element can be expressed as follows: 𝐌^𝑎₁𝑘,𝑞=𝑀^𝑎_3,𝑥=0,𝑞−𝑘−1,𝑞−𝑘_𝑀𝑚=1𝑚_3,𝑧𝑚exp𝑗−𝛾_𝑥𝑚+𝜑_𝑥𝑚×exp𝑗2𝑘−𝑞𝜑_𝑥𝑚. (10)

In a matrix form, (10) can be written as 𝐌^𝑎₁=𝐀𝐌^𝑎_𝑧𝛀^∗₁Φ₁𝐀^𝐻, (11) where the superscript “𝐻” represents conjugate transpose operation, 𝐌^𝑎_𝑧 denotes the third-order moment matrix of source signals, and 𝐀=𝐚₁,𝐚₂,…,𝐚_𝑀,𝐚_𝑚=1,exp𝑗2𝜑_𝑥𝑚,…,exp𝑗2(𝑁−1)𝜑_𝑥𝑚,𝛀₁=diagexp𝑗𝛾_𝑥1,exp𝑗𝛾_𝑥2,…,exp𝑗𝛾_𝑥𝑀,Φ₁=diagexp𝑗𝜑_𝑥1,exp𝑗𝜑_𝑥2,…,exp𝑗𝜑_𝑥𝑀. (12)

On the other hand, following the same process described above, we can easily obtain 𝑀^𝑎_3,𝑥(0,𝑙+1,𝑙)=𝑚^𝑀_3,𝑢𝑚=1𝑚_3,𝑧𝑚exp𝑗𝛾_𝑥𝑚+𝜑_𝑥𝑚exp𝑗2𝑙𝜑_𝑥𝑚,𝑀^𝑎_3,𝑥(1,−𝑙−1,−𝑙)=𝑚^𝑀_3,𝑢𝑚=1𝑚_3,𝑧𝑚exp𝑗−𝛾_𝑥𝑚+𝜑_𝑥𝑚×exp𝑗𝛾_𝑦𝑚+𝜑_𝑦𝑚exp𝑗2𝑙𝜑_𝑥𝑚,𝑀^𝑎_3,𝑥(−1,−𝑙−1,−𝑙)=𝑚^𝑀_3,𝑢𝑚=1𝑚_3,𝑧𝑚exp𝑗−𝛾_𝑥𝑚+𝜑_𝑥𝑚×exp𝑗−𝛾_𝑦𝑚+𝜑_𝑦𝑚exp𝑗2𝑙𝜑_𝑥𝑚. (13)

And the elements of the other three matrices can be expressed as 𝐌^𝑎₂𝑘,𝑞=𝑀^𝑎_3,𝑥=0,𝑘−𝑞+1,𝑘−𝑞_𝑀𝑚=1𝑚_3,𝑧𝑚exp𝑗𝛾_𝑥𝑚+𝜑_𝑥𝑚exp𝑗2𝑘−𝑞𝜑_𝑥𝑚,𝐌^𝑎₃𝑘,𝑞=𝑀^𝑎_3,𝑥=1,𝑞−𝑘−1,𝑞−𝑘_𝑀𝑚=1𝑚_3,𝑧𝑚exp𝑗−𝛾_𝑥𝑚+𝜑_𝑥𝑚×exp𝑗𝛾_𝑦𝑚+𝜑_𝑦𝑚exp𝑗2𝑘−𝑞𝜑_𝑥𝑚,𝐌^𝑎₄𝑘,𝑞=𝑀^𝑎_3,𝑥=−1,𝑞−𝑘−1,𝑞−𝑘_𝑀𝑚=1𝑚_3,𝑧𝑚exp𝑗−𝛾_𝑥𝑚+𝜑_𝑥𝑚×exp𝑗−𝛾_𝑦𝑚+𝜑_𝑦𝑚exp𝑗2𝑘−𝑞𝜑_𝑥𝑚. (14)

Finally, we have 𝐌^𝑎₂=𝐀𝐌^𝑎_𝑧𝛀₁Φ₁𝐀^𝐻,𝐌^𝑎₃=𝐀𝐌^𝑎_𝑧𝛀^∗₁Φ₁𝛀₂Φ₂𝐀^𝐻,𝐌^𝑎₄=𝐀𝐌^𝑎_𝑧𝛀^∗₁Φ₁𝛀₂Φ₂𝐀^𝐻 (15) with 𝛀₂=diagexp𝑗𝛾_𝑦1,exp𝑗𝛾_𝑦2,…,exp𝑗𝛾_𝑦𝑀,Φ₂=diagexp𝑗𝜑_𝑦1,exp𝑗𝜑_𝑦2,…,exp𝑗𝜑_𝑦𝑀. (16)

3.2. Build the Parallel Factor Analysis Model

Considering the situation of limited samples, we build a parallel factor analysis model that uses the third-order cyclic moments as 𝐌^𝑎=𝐌^𝑎()𝐌∶,∶,1_𝑎()𝐌∶,∶,2_𝑎()𝐌∶,∶,3_𝑎()=𝐌∶,∶,4^𝑎1𝐌^𝑎₂𝐌^𝑎₃𝐌^𝑎₄=𝐀𝐌^𝑎_𝑧𝛀^∗₁Φ₁𝐀^𝐻𝐀𝐌^a_z𝛀₁Φ₁𝐀^𝐻𝐀𝐌^𝑎_𝑧𝛀^∗₁Φ₁𝛀₂Φ₂𝐀^𝐻𝐀𝐌^𝑎_𝑧𝛀^∗₁Φ₁𝛀^∗₂Φ₂𝐀^𝐻+𝜈₁ (17) with 𝐂=𝐀^∗, and the Khatri-Rao product [15] for (17) shows 𝐌^𝑎=(𝐃⊗𝐀)𝐂^𝑇+𝜈₁, (18) where 𝑔𝐃=⁻¹𝐌^𝑎_𝑧𝛀^∗₁Φ₁𝑔⁻¹𝐌^𝑎_𝑧𝛀₁Φ₁𝑔⁻¹𝐌^𝑎_𝑧𝛀^∗_𝟏Φ₁𝛀₂Φ₂𝑔⁻¹𝐌^𝑎_𝑧𝛀^∗_𝟏Φ₁𝛀^∗_𝟐Φ₂, (19) with 𝑔⁻¹(𝐌^𝑎_𝑧𝛀^∗₁Φ₁) denoting a row vector consisting of diagonal matrix 𝐌^𝑎_𝑧𝛀^∗₁Φ₁.

Similarly, (18) also yields 𝐗=(𝐀⊗𝐂)𝐃^𝑇+𝜈₂,𝐘=(𝐂⊗𝐃)𝐀^𝑇+𝜈₃. (20)

3.3. Solve the Parallel Factor Analysis Model, and Estimate 3-D Parameters

As it stands, 𝐀 and 𝐂 are both Vandermonde matrices, and then they have Kruskal-rank (𝑘-rank) 𝑘_𝐀=𝑘_𝐂^𝑇=𝑀. On the other hand, the 𝑘-rank of 𝐃 will be 𝑘_𝐃=min(4,𝑀). When the condition that the number of signals being 𝑀≥2 holds, then 𝐀, 𝐂, and 𝐃 are unique up to permutation and scaling of columns. With trilinear alternating least squares regression, we obtain that 𝐂^𝑇=argmin_𝐂𝑇𝑀^𝛼−(𝐃⊗𝐀)𝐂^𝑇2_𝐹,𝐃^𝑇=argmin_𝐃𝑇𝐗−(𝐀⊗𝐂)𝐃^𝑇2_𝐹,𝐀^𝑇=argmin_𝐀𝑇𝐘−(𝐂⊗𝐃)𝐀^𝑇2_𝐹. (21)

Then using these estimates, we can associate with each pair (𝛾_𝑥𝑚,𝛾_𝑦𝑚,and𝜑_𝑦𝑚) as follows: 𝛾_𝑥𝑚=12𝐃()arg2,𝑚𝐃(),𝛾1,𝑚_𝑦𝑚=12𝐃()arg3,𝑚𝐃()𝐃()1,𝑚−arg4,𝑚𝐃(),𝜑1,𝑚_𝑦𝑚=12𝐃()arg3,𝑚𝐃()𝐃()1,𝑚+arg4,𝑚𝐃().1,𝑚 (22)

Finally, the sources parameters can be estimated as 𝛼_𝑚𝜆=asin𝛾2𝜋𝑑^2𝑥𝑚+𝛾²_𝑦𝑚^1/2𝜃_𝑚,𝛾=atan^𝑦𝑚𝛾_𝑥𝑚,𝑟_𝑚̃=𝜋𝑑²1−sin²𝛼_𝑚sin²𝜃_𝑚𝜆𝜑_𝑦𝑚. (23)(24)(25)

4. Computer Simulation Results

In this section, we explicit some simulation results to evaluate the performance of proposed method. For all examples, a uniform linear array with number of 15 sensors and element spacing 0.25𝜆 is displayed, where 𝜆 is the wavelength of the narrowband source signals. Two near-field or far-field equal power cyclostationary signals are impinging on the array with center frequency 0.25𝜋, and their envelopes can be modeled as exponentially distributed. To be compared, we simultaneously execute the method of [12], which is suitable for the 3-D near-field non-Gaussian stationary sources localization. The presented results are evaluated by the estimated root mean square error (RMSE) from the averaged results of 200 independent Monte-Carlo simulations.

In the first example, Two near-field sources are located at (35°, 46°, and 1/6𝜆) and (20°, 60°, and 2/5𝜆), respectively. The additive noise is time and spatially white Gaussian with zero mean and unit variance. In addition, the snapshot number is set equal to 1024. When SNR varies from 0 dB to 25 dB, the RMSE of the elevation, azimuth, and range estimations for two near-field cyclostationary sources can be shown in Figures 2, 3, and 4. For the comparison, the performance in the same situation for the fourth-order cumulants based method has also been displayed. From these three figures, it is obvious that the proposed method performs better in elevation, azimuth, and range estimation than the FFA based algorithm for all SNRs. What is more, the RMSE of range estimations for the first source that is closer to the array is less than the second one. This phenomenon is well agreement with the theoretical analysis that the sources closer to sensor array would hold a smaller standard deviation than the one far from the array [7].

Figure 2 

The RMSE of elevation estimations for two near-field sources using the proposed method and the cumulant based method versus SNR.

Figure 3 

RMSE of azimuth estimations for two near-field sources using the proposed method and the cumulant based method versus SNR.

Figure 4 

The RMSE of range estimations for two near-field sources using the proposed method and the cumulant based method versus SNR.

In the second example, the simulation condition is similar to the first example, except that the SNR is set at 10 dB, and the snapshot number is varied from 400 to 2000. The RMSE of elevation, azimuth, and range estimations of two near-field cyclostationary sources obtained from the proposed method as well as the cumulant based method can be displayed in Figures 5, 6, and 7. From these three figures, we can see that the proposed method still shows a more satisfactory accuracy than the cumulant based method in all available snapshots. In addition, the range estimations for the first sources are much better than the second one, while the elevation and azimuth estimations for both near-field sources are similar for the proposed method.

Figure 5 

The RMSE of elevation estimations for two near-field sources using the proposed method and the cumulant based method versus snapshot number.

Figure 6 

The RMSE of azimuth estimations for two near-field sources using the proposed method and the cumulant based method versus snapshot number.

Figure 7 

The RMSE of range estimations for two near-field sources using the proposed method and the cumulant based method versus snapshot number.

In the third example, the additive noise is time and spatially white Gaussian random process and the colored noise from a second-order AR model with parameter (−1.8,0.9), respectively. In addition, the other simulation condition is similar to the first example, except that the SNR is set at 10 dB, and snapshot number is equal to 1024. The mean and variance of the elevation, azimuth, and range estimations for two near-field cyclostationary sources can be shown in Table 1. From this table, we can easily see that the effectiveness of the proposed method is not little affected by the difference of the sensor noise.

Table 1 

RMSE of 3D parameters for near-field localization in the third example.

5. Conclusion

This paper considers the problem of the passive localization of 3-D near-field cyclostationary sources and proposes an efficient third-order cyclic moment based algorithm. The construction of the parallel factor analysis model in the cyclic domain avoids the multidimensional search and pairing parameters. Moreover, the utilization of trilinear alternating least squares regression deals with the joint estimation of elevation, azimuth, and range. From the simulation results mentioned above, the proposed method outperforms the cumulant based method in locating near-field cyclostationary sources; in addition, the stationary noise has a little influence on the estimated accuracy of the proposed method.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (61171137) and Specialized Research Fund for the Doctoral Program of Higher Education (20090061120042).