# Estimation Method for Mixed Near-Field and Far-Field Sources Based on Four-Order Cumulant and Polarization Sensitive Array

Authors: Liguo Wang; Binbin Yu
Journal: Security and Communication Networks
 (2022)

---

## Abstract

A new multidimensional parameters joint estimation method of mixed near-field and far-field sources based on polarization sensitive array, fourth-order cumulant, joint diagonalization technology, and propagator method is presented, which can realize the joint estimation of DOA (Direction of Arrival), range, frequency, polarization auxiliary angle, and polarization phase difference without multidimensional spectral peak searching and parameter pairing, and it is suitable for any additive Gaussian noise environment and effective for reducing the loss of array aperture. The algorithm skillfully constructs the fourth-order cumulant matrix by using the output on the label of specific dipole pairs of the received array, which effectively avoids the matrix rank reduction caused by the far-field sources coexistence situation. In addition, the presented algorithm utilizes the orthogonal propagation algorithm for subspace decomposition and uses the total least square solution to replace the orthogonal solution of singular value decomposition, which effectively reduces the computational complexity. The experiment proved the effectiveness of the proposed algorithm.


---

## Body

## 1. Introduction

Source parameter estimation is widely used in radar, sonar, indoor positioning, military confrontation, and many other fields. It can be separated into two types, namely, far-field source and near-field source. For the far-field source, the signal is far away from the array, and it is in the Fraunhofer region, which meets the requirements ofr≥2D2/λ. This moment, the wavefront of the signal is described as the plane wave, and the position information of the signal is determined only by DOA. However, for near-field source, the signal is close to the receiving array in Fresnel region, which meets the requirements of r<2D2/λ, and the wavefront is described as the spherical wave, so the range parameter of the signal cannot be ignored. At this moment, the position information of the signal is described by DOA and range. The scholars at home and abroad have conducted in-depth study on the parameter estimation problem of far-field source and near-field source and presented many good algorithms [1–6]. Nevertheless, in many practical applications, the received signals contain both far-field source and near-field source simultaneously; this moment, the traditional algorithm will not work correctly.Many scholars have studied the parameter estimation problem of mixed far-field and near-field sources in recent years. Xu et al. presented a multiparameters estimation method of mixed sources based on fourth-order cumulants. By selecting the output of a specific sequence number array element to construct the fourth-order cumulant matrix, the algorithm effectively avoided the problem of matrix rank reduction due to the existence of far-field sources at the same time and suits for any Gaussian noise environment [7]. Liu et al. proposed a parameter estimation strategy of mixed sources based on spatial difference method. According to the characteristic structure differences of covariance matrix of far-field sources and near-field sources, the algorithm uses oblique projection technology and spatial difference technology to restrain the influence of far-field sources [8]. Zheng et al. presented a mixed source parameter estimation algorithm based on mixed order statistics and MUSIC. The second-order statistics is used to estimate the DOA of far-field sources, and the four-order cumulant is used to estimate the DOA and range parameters of near-field [9]. Liang et al. presented a two-step MUSIC algorithm for mixed sources parameters estimation. The algorithm needs to not only construct high-order cumulant matrices, but also conduct MUSIC spectral peak searches, so the computation is huge [10]. Wang et al. designed an enhanced symmetric nested array, constructed a special cumulant matrix to achieve more conservative lags and more unique lags, and estimated the DOA of the mixed source, and then used one-dimensional spectral peaks to search the distance parameters of the near-field source. The algorithm proves that the designed enhanced symmetric nested array has excellent performance in the accuracy and resolution of DOA and distance [11]. Huang et al. proposed a low complexity far-field and near-field hybrid source location algorithm based on discrete Fourier transform (DFT) and orthogonal matching pursuit (OMP). By DFT on the antiangle elements of the covariance matrix of the array received signal, the DOAs of the mixed sources are estimated, and the far-field and near-field mixed sources are classified. Then, the sparse reconstruction problem related to the distance parameters is constructed by using the estimated angle parameters and covariance and solved by OMP algorithm. Experimental results show the superiority of the algorithm [12]. Li et al. introduced the large aperture coprime array into the DOA estimation of exponential coherent distributed sources, combined the sparsity of array space domain with the fourth-order cumulant characteristics of signals, expanded the array aperture, and improved the estimation accuracy and degree of freedom. Experiments show that when the array elements are the same, the accuracy and degree of freedom of the algorithm are higher than those of distributed signal parameter estimator (DSPE) algorithm and rotation invariant least squares algorithm [13]. Qian et al. studied the high-dimensional search of noncircular sources by using multinested arrays and proposed a reduced dimensional subspace data fusion algorithm to reduce the complexity and improve the estimation accuracy. The algorithm improves the localization positioning accuracy when the spatial degree of freedom is higher [14]. Su et al. constructed a convolution neural network for near-field and far-field mixed source classification and location by using the geometric structure of symmetrical nested array. Simulation results show that the proposed learning based method can improve the accuracy of mixed source localization [15]. He et al. developed a far-field and near-field mixed source localization algorithm based on accurate spatial geometry. The algorithm does not need to simplify the time-delay model of the received signal but establishes a unified nonapproximate model framework to solve it with a very simple mathematical method. The characteristic of the algorithm is that it does not need signal separation or isolation to classify the source type and is suitable for linear arrays with arbitrary spacing and can adapt to arbitrary and possibly unknown propagation loss [16]. Tian et al. studied the mixed sources localization problem when the number of array elements is large and has the same order of magnitude as the number of snapshots. The algorithm uses the phase compensation result of peak covariance matrix to construct a modified spectral function to estimate the DOA estimation of far-field source. Then, the oblique projection operation is used to extract the near-field source, and the DOA and range estimation of the near-field source are realized through two one-dimensional (1-D) spectral searches [17]. Zheng et al. proposed a mixed source localization method based on symmetric permutation coprime array. The array is composed of two permutation coprime arrays and a shared subarray. For a given number of sensors, it has a closed form expression of sensor position and continuous virtual sensors. The array configuration further reduces the mutual coupling between sensors, so it has better estimation performance [18]. Amir et al. presented an algorithm for multiple mixed sources based on component separation. The algorithm constructed a special cumulant matrix and extracted noncoherent DOA by MUSIC method. Another cumulant matrix is constructed to estimate the range and signal classify the incoherent sources [19]. Wu et al. presented a meshless mixed sources parameter estimation algorithm by low rank matrix reconstruction. Two special matrices are constructed by fourth-order cumulants, which are only related to DOA or range. Then, two reconstruction models of low rank matrix are established to estimate the DOA and range [20].In this paper, a new low complexity parameter estimation algorithm for far-field and near-field mixed sources is presented based on fourth-order cumulant, joint diagonalization technology, and orthogonal propagator method (OPM) [21, 22]. The proposed algorithm adopts the polarization sensitive array composed of 2P+1 dipole pairs, which can receive the two-dimensional electric field vector. Compared with scalar sensor, polarization sensitive array is characterized by strong anti-interference ability, robust detection ability, and high resolution. The algorithm is implemented in two steps: in the first step, four cumulant matrices were constructed by the output of dipole pairs of specific labels to estimate the DOA, frequency parameters, and polarization parameters of mixed source by the joint diagonalization method [23]. In the second step, we constructed a second-order statistics matrix by the output of polarization array, then decomposed it utilizing the OMP to calculate the noise subspace, finally, adopted the MUSIC algorithm to estimate the range of the near-field source and, at the same time, separated the near-field source and far-field source. The algorithm’s merit is low complexity, which is manifested in two points. In the first aspect, the algorithm adopted the joint diagonalization technique, which can reduce the dimension of the high-dimensional cumulant matrices and then cut down the computational burden. In the second aspect, orthogonal propagator method is used to reduce the computational complexity of subspace decomposition. The algorithm is only applicable to the case where the DOAs of near-field source and far-field source are different.


## 2. Signal Model

Suppose there areL mutually independent mixed far-field and near-field sources incident on the uniform polarization sensitive array, composed of 2P+1 dipole pairs, and spacing of the adjacent dipole pairs satisfies d≤λ/4. Figure 1 is the structure diagram of receiving array, θl is the DOA of the mixed sources, and rl is the range parameter of the near-field source. Selecting the dipole pairs where the coordinate origin is located as the phase reference point, slt is the l-th narrowband stationary source, and nmt is the noise on the m-th dipole pairs, τml represents the phase difference in which the l-th signal source is incident on between the reference dipole pair (0,0) and the m-th dipole pair, and then the received signal of the m-th dipole pair in x direction and y direction can be written as follows:where γl is named polarization auxiliary angle, ηl is named polarization phase difference, and vl can be expressed as(1)zmxt=∑l=1Lslt+τmlθl,rlvlxγl,ηl+nmxtm=−P,…,0,…,P,(2)zmyt=∑l=1Lslt+τmlθl,rlvlyγl,ηl+nmytm=−P,…,0,…,P,where γl is named polarization auxiliary angle, ηl is named polarization phase difference, and vl can be expressed as(3)vl=vlxvly=sinγlcosθlcosφlejηl−cosγlsinφlsinγlcosθlsinφlejηl+cosγlcosφl.Figure 1
Structure diagram of receiving array.Since the receiving array is a uniform linear array, and the signal is located in thexoy plane, in this case φl=π/2; then the above formula can be reduced to(4)vl=vlxvly=−cosγlsinγlcosθlejηl.For the near-field source,τml, can be written as(5)τml=2πλlrl2+md2−2rlmdsinθl−rl,where λl, θl, and rl represent the wavelength, DOA, and range parameters of the mixed sources, respectively.By the Fresnel approximation, the above formula will be simplified as(6)τm,l≈2πrlλm2d22rl2−mdsinθlrl−m2d2sin2θl2rl2≈ωlm+ϕlm2+Od2rl2≈ωlm+ϕlm2,where(7)ωl=−2πλldsinθl,(8)ϕl=πd2λlrlcos2θ.In this case, formulas (1) and (2) can be rewritten as(9)zmxt=∑l=1Lsltej2πfltejωlm+ϕlm2vlxγl,ηl+nmxt,zmyt=∑l=1Lsltej2πfltejωlm+ϕlm2vlyγl,ηl+nmyt.However, for the far-field sources, the source is far away from the receiving array, and it can be regarded as a near-field source with infinite range. That is,rl=∞, and ϕl=0. So, formula (6) can be rewritten as(10)τm,l=ωlm.Then, formulas (1) and (2) can be rewritten as(11)zmxt=∑l=1Lsltej2πfltejωlmvlxγl,ηl+nmxt,zmyt=∑l=1Lsltej2πfltejωlmvlyγl,ηl+nmyt.According to the above analysis, the matrix form of the mixed source model can be written as(12)Zx,yt=Astvlx,y+Nx,yt=AFsFtvlx,y+ANsNtvlx,y+Nx,yt,where A=AFAN, and s=sFTsNT. AF and AN are the direction matrix, and their expressions are as follows:(13)AN=aNω1,ϕ1,...,aNωK,ϕK,AF=aFωK+1,...,aFωL.In the above formulas,K is the number of near-field sources, and L-K is the number of far-field sources. aFωl and aNωl,ϕl represent the l-th column of the far-field source and near-field sources’ direction matrix, respectively, which can be express as(14)aNωl,ϕl=ej−ωlM+ϕlM2,...,1,...,ejωlM+ϕlM2T,aFωl=e−jMωl,...,1,...,ejMωlT.The superscriptT represents the transpose of the matrix. sNt and sFt are the signal vectors.(15)sNt=s1t,...,sKtT,sFt=sK+1t,...,sLtT.For the convenience of theoretical derivation, three lemmas are given below.Lemma 1 (see [23]).
The matrixB and the matrix E are essentially equal; if there is a matrix F, such that B=EF, where there is only one nonzero element in each row and column and its module equals 1, then B=˙E.Lemma 2 (see [23]).
The essential and unique condition of joint diagonalization is thatG=G1,G2 is the set of two matrices. For l=1,2,⋯,L, every matrix Gl can be expressed in the form of UΦlUH, and U is a unitary matrix. Then any joint diagonalizer of G˜ is essentially equal to the matrix U, only if(16)∀≤p≠q≤L∃l=1,2,…,Lφpl≠φql.This paper makes the following assumptions:(1)
The number of mixed sources is known, and it is a zero mean, statistically independent stochastic narrowband stationary process(2)
The noise is additive white Gaussian noise with a mean value of zero and variance of one, and it is uncorrelated with the mixed source(3)
θi≠θj, ωi≠ωj, andϕi≠ϕj, if i≠j(4)
The range between the receiving array dipole pairsd and the number of sources meet the requirements of L<2P+1


## 3. Proposed Algorithm

According to formula (4) and formula (10), it can be concluded that the parameter estimation problem of far-field source is a special case of the near-field source; it is the near-field source parameter estimation problem when the range approaches infinity.Their phase difference has essentially the same form. For the far-field source, the value of the parameterϕl is zero, and ωl is the common parameter of them, which contains DOA information. Therefore, the idea of the algorithm is to construct four fourth-order cumulant matrices by the output on different polarization array dipole pairs to realize the separation of angle and range. The joint diagonalization technology is used to estimate the DOAs, carrier frequency, and polarization parameters of mixed sources directly and then utilizes the orthogonal propagator method to estimate the ranges of near-field and classify the near-field and far-field.


### 3.1. DOA Estimation of Mixed Sources

The traditional array signal parameter estimation algorithms mostly use the second-order statistical characteristics of the signal. Four-order cumulants are insensitive to Gaussian process and have many good properties in mathematical form. Therefore, in dealing with Gaussian noise, four-order cumulants have better estimation performance than second-order statistics and can solve the problem that second-order statistics are difficult to solve, such as Gaussian colored noise.According to the property of fourth-order cumulant [6], define the following matrix C1:(17)C1m,n=cumzmxt,z−mxt∗,znyt∗,z−nyt=∑w=1L∑g=1L∑k=1L∑l=1Lcumswtej2πfwtejmωw+m2ϕwvwx,sg∗te−j2πfgtej−mωg+m2ϕgvgx∗,sk∗te−j2πfkte−jnωk+n2ϕkvkx∗,sltej2πfltej−nωl+n2ϕlvlx=∑l=1Lvlx2vly2c4slej2m−nωl.Since the fourth-order cumulant of Gaussian noise is zero, there is no noise term in the above formula.It can be expressed in the form of matrix(18)C1m,n=BωlC4sBHωl,where(19)C4s=diagv1x2v1y2c4s1,…,vLx2vLy2c4sL,B=bω1,…,bωL2P+1×L,c4sl=cumsl∗t,slt,sl∗t,slt,bωl=e−j2Pωl,e−j2P−1ωl,…,1,…,ej2PωlT.It can be seen from assumption (1) that since the signal source is a narrowband stationary signal source andzt is a slowly varying amplitude modulation function (also known as the real envelope function), there is zt+1≈zt. Using the time lag of the received signal and rewriting the fourth-order cumulant matrix C1, the following matrix can be obtained:(20)C2m,n=cumzmxt+1,z−mxt∗,znyt∗,z−nyt=∑l=1Lvlx2vly2c4,slej2πflej2m−nωl.The matrix form is(21)C2m,n=BωlC4sΦ1BHωl,where(22)Φ1=diagej2πf1,...,ej2πfL.The superscript H represents the conjugate transpose of the matrix.Similarly, we constructed the following four-order cumulant matrices:(23)C3m,n=cumzm−1xt,z1−mxt∗,znyt∗,z−nyt=∑l=1Lvlx2vly2c4,sle−j2ωlej2m−nωl,C4m,n=cumzmxt∗,z−myt,znyt∗,z−nyt=−∑l=1Lvlx2vly2c4,slcosθltanγle−jηlej2m−nωl.The matrices form is(24)C3m,n=BωlC4sΦ2BHωl,C4m,n=BωlC4sΦ3BHωl,where(25)Φ2=diage−j2ω1,...,e−j2ωL,Φ3=diag−cosθ1tanγ1e−jη1,...,−cosθLtanγLe−jηL.The joint diagonalization of multiple matrices was first proposed to consider the common Principal Component Analysis ofK covariance matrices [24]. Later, Dr. Cardoso and Dr. Belochrani proposed the joint diagonalization of multiple cumulant matrices and covariance matrices from the scale of blind signal separation [23, 25]. The mathematical problem of joint diagonalization is as follows: given K symmetric matrices A1,A2,...,AK of m×m dimensional, find a full rank matrix U, so that K matrices diagonalize at the same time, that is,(26)Ak=UΛkUH,k=1,2,...,K,where U is named the Joint diagonalizer.It should be pointed out that the (exact) joint diagonalization of twon×n Hermitian matrices A and B is equivalent to the generalized eigenvalue decomposition of Hermitian matrix pencil A,B.Joint diagonalizationAk=UΛkUH is exact joint diagonalization. However, the actual joint diagonalization is approximate joint diagonalization; given the matrix set Ak=A1,…,Ak, we want to find a joint diagonalizer U∈ℂm×n and K corresponding n×n diagonal matrices, to minimize the objective function [25]:(27)minJ1U,A1,,AK=min∑k=1KωkUAHkU−ΛkF2,or [26](28)minJ2U,A1,,AK=min∑k=1KωkAk−UΛkUHF2,where ω1,ω2,⋯,ωK is a positive weight coefficient.The joint diagonalization algorithm can be implemented in two steps. Firstly, the whitening matrix is constructed to whiten other matrices, and then the joint diagonalization is realized.For simplifying derivation process, the output of the array is normalized:(29)yt=Ast+nt=∑l=1Lalclclslt+nt.That is to say, the amplitude of thel-th source is induced to the direction vector al.We carried out the eigenvalue decomposition for the matrixC1 and constructed the whitened matrix W with the maximum eigenvalue, the eigenvectors corresponding to the maximum eigenvalue, and the average value of the minimum eigenvalue(30)W=ε1−σ2−1/2h1,...,εL−σ2−1/2hLH.In the above formula,ε1,...,εL is the L maximum eigenvalues of matrix C1, h1,...,hL is the corresponding eigenvector, and σ is the average value of small eigenvalues and satisfies the following conditions:(31)WC1WH=WBBHWH=I.It can be inferred from the above formula that if matrixW is whitening matrix, then WA is unitary matrix with L×L dimension. That is to say, for any whitening matrix W, there is an L×L dimensional unitary matrix satisfying U=WB, then the matrix U can be expressed as B=W†U, and † denote pseudo inverse. This also makes the original problem from the determination of 2N+1×L dimension matrix B to the L×L dimension matrix U, so as to achieve the purpose of reducing the dimension of the matrix and the computational complexity.Next, perform the following transformation:(32)H1=WC2WH=UΦ1UH,H2=WC3WH=UΦ2UH,H3=WC4WH=UΦ3UH.IfHi=UΦiUH, 1,2,3, then it can only be diagonalized by a union diagonalizer essentially equal to U, and there exists V=˙U such that(33)LG˜,V=∑i=13offVHGiV,where offM means the squares sum of absolute values of nonprincipal diagonal elements of M.The joint diagonalization of matricesH1, H2, and H3 can be realized through a series of given rotations. The specific process can refer to [23].LetV1, V2, and V3 be the eigenvalue matrix obtained by joint diagonalization, that is, the estimated value of matrix Ri. Then the phase parameters of the mixed sources can be obtained according to formula (22) and formula (25).(34)ω∧l=−0.5argV1,f∧l=argV2.According to formula (7) and formula (25), the DOA parameters, polarization auxiliary angle, and polarization phase difference will be calculated.(35)θ∧l=arcsinγ^lλl2πd,γ∧l=tan−1V∧3llcosθ∧l,η∧l=−arg−V∧3ll.


### 3.2. Mixed Sources Classification and Range Estimation of Near Field

Due to the DOA, frequency and polarization parameters of the mixed source have been estimated; two tasks need to be completed in the next step: the first task is to classify the estimated DOA, frequency, and polarization parameters, that is, to determine which are the near-field sources and far-field sources; the second task is to estimate the range parameters of the near-field sources. For the second task, the traditional MUSIC algorithm can be used; that is, the MUSIC spatial spectrum function is constructed by noise subspace of the data received by the array, and the DOA and range parameters of near-field sources are estimated through multidimensional spectral peak search.However, we have estimated the DOA, frequency, and polarization parameters in the first stage, so substituting these parameters into the MUSIC spatial spectral function, the multidimensional spectral peak search is transformed into unidimensional search. In this process,θ∧l and r∧l are automatically paired without additional pairing program, because the MUSIC space spectrum function corresponding to the DOA of a certain source searches the range satisfying r∧l∈0.62D3/λ1/2,2D2/λ; then the DOA is a near-field source; if satisfying r∧l∈2D2/λ,+∞, then the source corresponding to the DOA is a far-field source, thus realizing the source classification.Usually, MUSIC algorithm mainly uses the array to receive the signal, construct the covariance matrix or high-order cumulant matrix, perform singular value decomposition, estimate the noise subspace, construct the MUSIC spectral function, and then estimate the source parameters through spectral peak search. In this course, whether constructing high-dimensional covariance matrix or high-order cumulant matrix, or constructing noise subspace by singular value decomposition, the amount of calculation is quite large. For the sake of solving the problem, we use the propagator method as an effective mean.The propagator is a linear operator based upon a partition of the steering vectors. The propagator method is a fast subspace decomposition method, which does not need high-dimensional singular value decomposition but utilizes the total least squares method to estimate the noise subspace, thus greatly reducing the computational complexity. The basis of the noise subspace estimated by the propagator method is not orthogonal solution. Compared with the MUSIC algorithm based on singular value decomposition, the performance of this algorithm is relatively poor.In order to solve this problem, Sylvie Marcos et al. improved the propagator method and proposed the orthogonal propagator method [21, 22], where the estimated solution of noise subspace is orthogonal solution. In the medium or high SNR, the estimation performance is similar to MUSIC algorithm based on singular value decomposition, but it still has advantages in the computational complexity. Therefore, the orthogonal propagator method is used to estimate the noise subspace in this paper.According to the signals received by the array, the following matrix is constructed:(36)M=EzxtzxtH=Aθl,rl,fl,γl,ηlRsAHθl,rl,fl,γl,ηl+σ2I.Divided the matrixM into two parts:(37)M=M1M2,where(38)M1=M1:L,:,M2=ML+1:2P+1,:,where L is the number of mixed sources.Because the matrixM1 is a row full rank matrix, there exists a unique propagator matrix Q which satisfies the below formula:(39)QHM1=M2.Estimate the matrixQ by the total least square method:(40)Q=M1HM1−1M1HM2.Define the matrixU:(41)U=QT,−I2P+1−KT,and then U is noise subspace and satisfies the following formula:(42)UHAf,θ,r,γ,η=O.Now, introduce a projection operator onto the noise subspace; the matrixU can be replaced by its orthonormalized version:(43)U0=UUHU−1/2.Then, the MUSIC spectral function can be expressed:(44)FOPM−MUSICf,θ,r,γ,η=1aHf,θ,r,γ,ηU0U0Haf,θ,r,γ,η.We substituted frequencies, DOAs, polarization auxiliary angle, and polarization phase difference into formula (44), and the above formula can be rewritten as follows:(45)r∧=argmax1aHf∧,θ∧,r,γ∧,η∧U0U0Haf∧,θ∧,r,γ∧,η∧.


### 3.3. Complexity Analysis

The algorithm adopts mixed order statistics. The algorithm is divided into two stages. In the first stage, the fourth-order cumulant is used to estimate the azimuth, frequency, and polarization parameters, and in the second stage, the second-order statistics is used to estimate the range and classify the near-field source and far-field. The main computational complexity of the algorithm is analyzed as follows: in the first stage, the computational complexity of constructing four fourth-order cumulant matrices isO(36(2N + 1)2K), the computational complexity of constructing one whitening matrix is O(2N + 1)3/3, the computational complexity of whitening reduction dimension for three matrices is O(12L(N+1)2)+O(6 L2(N+1)), and the computational complexity of jointly diagonalizing three matrices is O(3 N3). In the second stage, the computational complexity of constructing a second-order statistics matrix is O((2N + 1)2K), and the computational complexity of estimating the noise subspace by using the orthogonal propagation operator method is O(L(2N + 1)2), where K is the number of snapshots. Therefore, the main computational complexity is O(36(2N + 1)2K+(2N + 1)3/3 + 12L(N+1)2+6 L2(N+1)+3 N3+(2N + 1)2K + L(2N + 1)2).


## 3.1. DOA Estimation of Mixed Sources

The traditional array signal parameter estimation algorithms mostly use the second-order statistical characteristics of the signal. Four-order cumulants are insensitive to Gaussian process and have many good properties in mathematical form. Therefore, in dealing with Gaussian noise, four-order cumulants have better estimation performance than second-order statistics and can solve the problem that second-order statistics are difficult to solve, such as Gaussian colored noise.According to the property of fourth-order cumulant [6], define the following matrix C1:(17)C1m,n=cumzmxt,z−mxt∗,znyt∗,z−nyt=∑w=1L∑g=1L∑k=1L∑l=1Lcumswtej2πfwtejmωw+m2ϕwvwx,sg∗te−j2πfgtej−mωg+m2ϕgvgx∗,sk∗te−j2πfkte−jnωk+n2ϕkvkx∗,sltej2πfltej−nωl+n2ϕlvlx=∑l=1Lvlx2vly2c4slej2m−nωl.Since the fourth-order cumulant of Gaussian noise is zero, there is no noise term in the above formula.It can be expressed in the form of matrix(18)C1m,n=BωlC4sBHωl,where(19)C4s=diagv1x2v1y2c4s1,…,vLx2vLy2c4sL,B=bω1,…,bωL2P+1×L,c4sl=cumsl∗t,slt,sl∗t,slt,bωl=e−j2Pωl,e−j2P−1ωl,…,1,…,ej2PωlT.It can be seen from assumption (1) that since the signal source is a narrowband stationary signal source andzt is a slowly varying amplitude modulation function (also known as the real envelope function), there is zt+1≈zt. Using the time lag of the received signal and rewriting the fourth-order cumulant matrix C1, the following matrix can be obtained:(20)C2m,n=cumzmxt+1,z−mxt∗,znyt∗,z−nyt=∑l=1Lvlx2vly2c4,slej2πflej2m−nωl.The matrix form is(21)C2m,n=BωlC4sΦ1BHωl,where(22)Φ1=diagej2πf1,...,ej2πfL.The superscript H represents the conjugate transpose of the matrix.Similarly, we constructed the following four-order cumulant matrices:(23)C3m,n=cumzm−1xt,z1−mxt∗,znyt∗,z−nyt=∑l=1Lvlx2vly2c4,sle−j2ωlej2m−nωl,C4m,n=cumzmxt∗,z−myt,znyt∗,z−nyt=−∑l=1Lvlx2vly2c4,slcosθltanγle−jηlej2m−nωl.The matrices form is(24)C3m,n=BωlC4sΦ2BHωl,C4m,n=BωlC4sΦ3BHωl,where(25)Φ2=diage−j2ω1,...,e−j2ωL,Φ3=diag−cosθ1tanγ1e−jη1,...,−cosθLtanγLe−jηL.The joint diagonalization of multiple matrices was first proposed to consider the common Principal Component Analysis ofK covariance matrices [24]. Later, Dr. Cardoso and Dr. Belochrani proposed the joint diagonalization of multiple cumulant matrices and covariance matrices from the scale of blind signal separation [23, 25]. The mathematical problem of joint diagonalization is as follows: given K symmetric matrices A1,A2,...,AK of m×m dimensional, find a full rank matrix U, so that K matrices diagonalize at the same time, that is,(26)Ak=UΛkUH,k=1,2,...,K,where U is named the Joint diagonalizer.It should be pointed out that the (exact) joint diagonalization of twon×n Hermitian matrices A and B is equivalent to the generalized eigenvalue decomposition of Hermitian matrix pencil A,B.Joint diagonalizationAk=UΛkUH is exact joint diagonalization. However, the actual joint diagonalization is approximate joint diagonalization; given the matrix set Ak=A1,…,Ak, we want to find a joint diagonalizer U∈ℂm×n and K corresponding n×n diagonal matrices, to minimize the objective function [25]:(27)minJ1U,A1,,AK=min∑k=1KωkUAHkU−ΛkF2,or [26](28)minJ2U,A1,,AK=min∑k=1KωkAk−UΛkUHF2,where ω1,ω2,⋯,ωK is a positive weight coefficient.The joint diagonalization algorithm can be implemented in two steps. Firstly, the whitening matrix is constructed to whiten other matrices, and then the joint diagonalization is realized.For simplifying derivation process, the output of the array is normalized:(29)yt=Ast+nt=∑l=1Lalclclslt+nt.That is to say, the amplitude of thel-th source is induced to the direction vector al.We carried out the eigenvalue decomposition for the matrixC1 and constructed the whitened matrix W with the maximum eigenvalue, the eigenvectors corresponding to the maximum eigenvalue, and the average value of the minimum eigenvalue(30)W=ε1−σ2−1/2h1,...,εL−σ2−1/2hLH.In the above formula,ε1,...,εL is the L maximum eigenvalues of matrix C1, h1,...,hL is the corresponding eigenvector, and σ is the average value of small eigenvalues and satisfies the following conditions:(31)WC1WH=WBBHWH=I.It can be inferred from the above formula that if matrixW is whitening matrix, then WA is unitary matrix with L×L dimension. That is to say, for any whitening matrix W, there is an L×L dimensional unitary matrix satisfying U=WB, then the matrix U can be expressed as B=W†U, and † denote pseudo inverse. This also makes the original problem from the determination of 2N+1×L dimension matrix B to the L×L dimension matrix U, so as to achieve the purpose of reducing the dimension of the matrix and the computational complexity.Next, perform the following transformation:(32)H1=WC2WH=UΦ1UH,H2=WC3WH=UΦ2UH,H3=WC4WH=UΦ3UH.IfHi=UΦiUH, 1,2,3, then it can only be diagonalized by a union diagonalizer essentially equal to U, and there exists V=˙U such that(33)LG˜,V=∑i=13offVHGiV,where offM means the squares sum of absolute values of nonprincipal diagonal elements of M.The joint diagonalization of matricesH1, H2, and H3 can be realized through a series of given rotations. The specific process can refer to [23].LetV1, V2, and V3 be the eigenvalue matrix obtained by joint diagonalization, that is, the estimated value of matrix Ri. Then the phase parameters of the mixed sources can be obtained according to formula (22) and formula (25).(34)ω∧l=−0.5argV1,f∧l=argV2.According to formula (7) and formula (25), the DOA parameters, polarization auxiliary angle, and polarization phase difference will be calculated.(35)θ∧l=arcsinγ^lλl2πd,γ∧l=tan−1V∧3llcosθ∧l,η∧l=−arg−V∧3ll.


## 3.2. Mixed Sources Classification and Range Estimation of Near Field

Due to the DOA, frequency and polarization parameters of the mixed source have been estimated; two tasks need to be completed in the next step: the first task is to classify the estimated DOA, frequency, and polarization parameters, that is, to determine which are the near-field sources and far-field sources; the second task is to estimate the range parameters of the near-field sources. For the second task, the traditional MUSIC algorithm can be used; that is, the MUSIC spatial spectrum function is constructed by noise subspace of the data received by the array, and the DOA and range parameters of near-field sources are estimated through multidimensional spectral peak search.However, we have estimated the DOA, frequency, and polarization parameters in the first stage, so substituting these parameters into the MUSIC spatial spectral function, the multidimensional spectral peak search is transformed into unidimensional search. In this process,θ∧l and r∧l are automatically paired without additional pairing program, because the MUSIC space spectrum function corresponding to the DOA of a certain source searches the range satisfying r∧l∈0.62D3/λ1/2,2D2/λ; then the DOA is a near-field source; if satisfying r∧l∈2D2/λ,+∞, then the source corresponding to the DOA is a far-field source, thus realizing the source classification.Usually, MUSIC algorithm mainly uses the array to receive the signal, construct the covariance matrix or high-order cumulant matrix, perform singular value decomposition, estimate the noise subspace, construct the MUSIC spectral function, and then estimate the source parameters through spectral peak search. In this course, whether constructing high-dimensional covariance matrix or high-order cumulant matrix, or constructing noise subspace by singular value decomposition, the amount of calculation is quite large. For the sake of solving the problem, we use the propagator method as an effective mean.The propagator is a linear operator based upon a partition of the steering vectors. The propagator method is a fast subspace decomposition method, which does not need high-dimensional singular value decomposition but utilizes the total least squares method to estimate the noise subspace, thus greatly reducing the computational complexity. The basis of the noise subspace estimated by the propagator method is not orthogonal solution. Compared with the MUSIC algorithm based on singular value decomposition, the performance of this algorithm is relatively poor.In order to solve this problem, Sylvie Marcos et al. improved the propagator method and proposed the orthogonal propagator method [21, 22], where the estimated solution of noise subspace is orthogonal solution. In the medium or high SNR, the estimation performance is similar to MUSIC algorithm based on singular value decomposition, but it still has advantages in the computational complexity. Therefore, the orthogonal propagator method is used to estimate the noise subspace in this paper.According to the signals received by the array, the following matrix is constructed:(36)M=EzxtzxtH=Aθl,rl,fl,γl,ηlRsAHθl,rl,fl,γl,ηl+σ2I.Divided the matrixM into two parts:(37)M=M1M2,where(38)M1=M1:L,:,M2=ML+1:2P+1,:,where L is the number of mixed sources.Because the matrixM1 is a row full rank matrix, there exists a unique propagator matrix Q which satisfies the below formula:(39)QHM1=M2.Estimate the matrixQ by the total least square method:(40)Q=M1HM1−1M1HM2.Define the matrixU:(41)U=QT,−I2P+1−KT,and then U is noise subspace and satisfies the following formula:(42)UHAf,θ,r,γ,η=O.Now, introduce a projection operator onto the noise subspace; the matrixU can be replaced by its orthonormalized version:(43)U0=UUHU−1/2.Then, the MUSIC spectral function can be expressed:(44)FOPM−MUSICf,θ,r,γ,η=1aHf,θ,r,γ,ηU0U0Haf,θ,r,γ,η.We substituted frequencies, DOAs, polarization auxiliary angle, and polarization phase difference into formula (44), and the above formula can be rewritten as follows:(45)r∧=argmax1aHf∧,θ∧,r,γ∧,η∧U0U0Haf∧,θ∧,r,γ∧,η∧.


## 3.3. Complexity Analysis

The algorithm adopts mixed order statistics. The algorithm is divided into two stages. In the first stage, the fourth-order cumulant is used to estimate the azimuth, frequency, and polarization parameters, and in the second stage, the second-order statistics is used to estimate the range and classify the near-field source and far-field. The main computational complexity of the algorithm is analyzed as follows: in the first stage, the computational complexity of constructing four fourth-order cumulant matrices isO(36(2N + 1)2K), the computational complexity of constructing one whitening matrix is O(2N + 1)3/3, the computational complexity of whitening reduction dimension for three matrices is O(12L(N+1)2)+O(6 L2(N+1)), and the computational complexity of jointly diagonalizing three matrices is O(3 N3). In the second stage, the computational complexity of constructing a second-order statistics matrix is O((2N + 1)2K), and the computational complexity of estimating the noise subspace by using the orthogonal propagation operator method is O(L(2N + 1)2), where K is the number of snapshots. Therefore, the main computational complexity is O(36(2N + 1)2K+(2N + 1)3/3 + 12L(N+1)2+6 L2(N+1)+3 N3+(2N + 1)2K + L(2N + 1)2).


## 4. Simulation and Results

In order to verify the effectiveness, the proposed algorithm is compared with [7] and [9]. The simulation conditions are as follows: assuming that two equal power and statistically independent narrowband sources are incident on the receiving array, the array is a uniform linear polarization sensitive array which is composed of 2P+1 = 15 omnidirectional dipole pairs; the carrier frequencies are {f1 = 2 MHz, f2 = 4 MHz}, respectively. The polarization parameter of the first source is γ1=35∘,η1=55∘, the polarization parameter of the second source is γ2=45∘,η2=65∘, the noise received from the dipole pair is additive Gaussian white noise, and the mean equals zero and the variance equals one. The signal-to-noise ratio is expressed as 10log10σl2/σn2, where σl2 is the power of the signal and σn2 is the power of the noise.The units of frequency, DOA, and range are MHz, degree, and wavelength, respectively. Take the mean value of 500 independent experiments for each estimated value, and use the RMSE (Root Mean Square Error) as the standard of performance evaluation, as shown in formula (46).(46)RMSE=1Nmc∑r=1Nmcx∧r,l−xk2,where Nmc is the number of Monte Carlo experiments, x∧r,l represents the estimated value, and xl represents the true value.Suppose that the first source is a near-field source with the positionθ1=30∘,r1=1.5λ, and the second source is a far-field source at position θ2=45∘,r2=+∞. If the number of snapshots is set at 256 and the signal-to-noise ratio changes from 0 dB to 25 dB, the change curves of DOA, range, frequency, and polarization parameters are shown in Figures 2–6. When the signal-to-noise ratio is set at 10 dB and the number of snapshots is changed from 256 to 1536, the change curves of DOA, range, frequency, and polarization parameters are shown in Figures 7–11.Figure 2
The RMSE curve of the DOA in SNR.Figure 3
The RMSE curve of the range in SNR.Figure 4
The RMSE curve of the frequency in SNR.Figure 5
The RMSE curve of the polarization auxiliary angle in SNR.Figure 6
The RMSE curve of the polarization phase difference in SNR.Figure 7
The RMSE curve of the DOA in snapshots.Figure 8
The RMSE curve of the range in snapshots.Figure 9
The RMSE curve of the frequency in snapshots.Figure 10
The RMSE curve of the polarization auxiliary angle in snapshots.Figure 11
The RMSE curve of the polarization phase difference in snapshots.From Figures2 to 11, it can be concluded that the proposed algorithm can effectively classify mixed sources and can accurately estimate multidimensional parameters, such as DOAs, range, frequency, and polarization parameters, and the calculation performance is better than that of [7] and [9]. The reasons are as follows: firstly, the algorithm adopts the fourth-order cumulant, which can effectively remove the influence of Gaussian noise. Secondly, it utilizes the joint diagonalization method to directly estimate the DOA, range, frequency, and polarization parameters, without additional parameter pairing. Finally, the MUSIC spectral function of range parameter near-field is constructed by the orthogonal propagation operator method, and the one-dimensional spatial search method is utilized to classify the far-field source and near-field source and estimate the range parameter of near-field sources.


## 5. Conclusion

The paper proposed a new multidimensional parameter joint estimation algorithm for mixed far-field and near-field sources based on polarization sensitive array, fourth-order cumulant, and orthogonal and propagator method. The algorithm was divided into two stages. In the first stage, four fourth-order cumulant matrices were constructed according to the specific array elements, and the joint diagonalization technology was utilized to directly estimate the DOA, frequency, and polarization parameters without parameter pairing. Meanwhile at the same time, it can effectively avoid the phenomenon of rank reduction due to the coexistence of far-field sources. In the second stage, firstly, we constructed a correlation matrix by the second-order statistics and implemented the orthogonal propagator method to estimate the noise subspace of correlation matrix and used the total least square solution of noise subspace to replace the orthogonal solution, so as to reduce the computational complexity. Secondly, we constructed the MUSIC spatial spectral function and then substituted the DOA, frequency, and polarization parameters estimated in the first step into the spatial spectral function to obtain the one-dimensional MUSIC spatial spectral function only about the range. Finally, we performed one-dimensional spectral peak search to estimate the range and classified the near-field and far-field. The algorithm did not need parameter pairing and multidimensional spectral peak search and has low computational complexity. Simulation results show that the proposed algorithm has good estimation performance.



---

*Source: 2994794-2022-05-16.xml*