Space-time near field source three-dimensional parameter estimation method based on cumulant domain step-by-step decoupling

Abstract

The invention discloses a space-time near field source three-dimensional parameter estimation method based on cumulant domain step-by-step decoupling, and the method comprises the steps: receiving a signal through a cross array, and carrying out the phase shift through Fresnel approximation; the method comprises the following steps: constructing a four-order cumulant matrix by considering signal amplitude attenuation and combining with signals received by a cross array, and decomposing a three-dimensional parameter estimation problem into a one-dimensional parameter estimation problem through step-by-step decoupling: firstly, decomposing a pitch angle and an azimuth angle of an incident signal into included angles with an x axis and a y axis; then angle spectrum peak searchers constructed by a reduced-rank MUSIC algorithm are adopted for estimation, then distance spectrum peak searchers are constructed for distance estimation, and finally the pitch angle and the azimuth angle are calculated according to the estimated values of the included angles between the incident signals and the x axis and the y axis, so that three-dimensional parameter estimation is achieved. The method has the advantages that the influence of signal amplitude attenuation is considered, the high-order cumulant is used, three-dimensional parameter estimation is achieved, and the estimation precision is higher.

Description

Technical Field

[0001] The present invention relates to a near-field signal source positioning technology, and in particular to a space-time near-field source three-dimensional parameter estimation method based on cumulative quantity domain step-by-step decoupling. Background Art

[0002] The DOA (Direction Of Arrival) positioning problem of the signal source is an important topic in the field of array signal processing. This positioning problem is widely used in many fields such as wireless communication, radar, sonar, etc. DOA estimation can be divided into two cases: near-field and far-field direction finding according to the distance between the signal source and the array. In the far-field direction finding case, the distance between the signal source and the reference array element satisfies r>>2D 2 / λ, where r represents the distance between the signal source and the reference array element, D represents the array aperture, and λ is the signal wavelength. In this case, signal propagation can be regarded as a plane wave, and the positioning problem is related to the incident angle. However, as the array aperture increases and the signal frequency increases, the Rayleigh distance 2D 2 / λ increases, and when the distance between the signal source and the reference array element satisfies In this case, the signal source will be in the near field of the array. If the far-field direction finding method is still used, the performance of parameter estimation will be reduced. Therefore, it is necessary to remodel the array signal model in the near field and design a new algorithm. At this time, the signal propagation can no longer be approximated as a plane wave, but a spherical wave, and the positioning problem is related to both angle and distance.

[0003] Most of the existing near-field direction finding algorithms do not consider the influence of signal amplitude attenuation. When the signal source is close to the sensors in the array, the signal amplitudes received by different sensors are often quite different. If the influence of the difference in signal amplitude attenuation when the array element receives the signal is ignored, it may cause mismatch of the signal model, thereby reducing the performance of parameter estimation. Many algorithms for locating near-field signal sources have been proposed, such as the maximum likelihood method and the second-order statistics method. However, most of the above methods only consider the two-dimensional parameters of the near-field source, namely azimuth and distance. However, for the DOA positioning problem of near-field signal sources in three-dimensional scenes, it is necessary to jointly estimate three-dimensional parameters, namely azimuth, elevation and distance. In addition, most of the existing DOA estimation algorithms for near-field signal sources use the second-order statistical characteristics of the signal. However, using high-order cumulants can obtain more information about the received signal, and can well suppress Gaussian noise and obtain better parameter estimation performance. Therefore, considering the influence of signal amplitude attenuation, it is crucial to study the method of three-dimensional parameter estimation of near-field signal sources in the cumulant domain. Summary of the invention

[0004] The technical problem to be solved by the present invention is to provide a three-dimensional parameter estimation method for space-time near-field sources based on step-by-step decoupling of the cumulant domain, which takes into account the influence of signal amplitude attenuation and uses high-order cumulants to achieve three-dimensional parameter estimation with higher estimation accuracy.

[0005] The technical solution adopted by the present invention to solve the above technical problems is: a method for estimating three-dimensional parameters of space-time near-field sources based on step-by-step decoupling of cumulant domain, characterized by comprising the following steps:

[0006] Step 1: Establish a cross array in a three-dimensional coordinate system: The cross array is used as a signal receiving array, which includes a uniform linear array located on the x-axis of the three-dimensional coordinate system and a uniform linear array located on the y-axis of the three-dimensional coordinate system. The uniform linear array located on the x-axis is composed of M x The uniform linear array on the y-axis consists of M array elements evenly distributed. y The two uniform linear arrays share the central element at the origin, and the central element is used as the reference element of the cross array. Each uniform linear array is symmetrical with the reference element as the center, and the distance between two adjacent elements in each uniform linear array is d; where M x Represents the total number of uniform linear array elements located on the x-axis, M x =2N x +1, N x represents the number of single-sided array elements of the uniform linear array located on the x-axis, M y Represents the total number of uniform linear array elements located on the y-axis, M y =2N y +1, N y It represents the number of single-sided array elements of the uniform linear array on the y-axis. The total number of array elements of the cross array is M, M = M x +M y -1, d = λ / 8, λ represents the signal wavelength, and the array element is the sensor;

[0007] Step 2: Assume that there are K space-time independent narrowband near-field sources that simultaneously incident signals on each element in the cross array, and each element receives all the signals that arrive at it; the three-dimensional parameters of the K space-time independent narrowband near-field sources are expressed as (θ1, ψ1, r 0，1 ),…,(θ k , ψ k , r 0，k ),…,(θ K , ψ K , r 0，K ), the signal received by the array element at position (m, 0) in the uniform linear array on the x-axis is recorded as x m (t), The signal received by the array element at position (0, n) in the uniform linear array on the y-axis is denoted as y n (t), Among them, N x ≥K and N y ≥K,K≥1,1≤k≤K,θ k represents the elevation angle of the signal incident from the kth narrowband near-field source, ψ k represents the azimuth of the signal incident from the kth narrowband near-field source, r 0，k represents the distance between the kth narrowband near-field source and the center of the reference array element, m∈{-N x , -N x +1,…,N x -1, N x}, n∈{-N y , -N y +1,…,N y -1, N y}, t is the time variable, r m，k represents the distance between the kth narrowband near-field source and the center of the array element at position (m, 0), α k represents the angle between the signal incident from the kth narrowband near-field source and the x-axis, α k With θ k , k The relationship is cos(α k )=sin(θ k )cos(ψ k ), r n，k represents the distance between the kth narrowband near-field source and the center of the array element at position (0, n), β k represents the angle between the signal incident from the kth narrowband near-field source and the y-axis, β k With θ k , k The relationship is sin(β k )=sin(θ k )sin(ψ k ), p represents the a priori known signal amplitude attenuation exponent value, e represents the natural constant, j represents the imaginary unit, s k (t) represents the baseband signal received and demodulated by any array element from the signal incident by the kth narrowband near-field source, w m (t) represents the additive white Gaussian noise on the array element at position (m, 0), w n (t) represents the additive Gaussian white noise on the array element at position (0, n), γ m，k and γ n，k are the intermediate quantities introduced,

[0008] Then the phase shift r in formula (1) m，k -r 0，k Using the Fresnel approximation, formula (1) can be rewritten as: Similarly, for the phase shift r in formula (2), n，k -r 0，k Using the Fresnel approximation, formula (2) can be rewritten as: in,

[0009] Based on formula (3), the received signal vector of all array elements in the uniform linear array on the x-axis is recorded as x(t), x(t) = A x s(t)+w x (t); Similarly, based on formula (4), the received signal vectors of all array elements in the uniform linear array located on the y-axis are recorded as y(t), y(t) = A y s(t)+w y (t); where A x represents the array manifold matrix of the uniform linear array located on the x-axis, A x =[a x (α1, r 0，1 ), ..., a x (α k , r 0，k ), ..., a x (α K , r 0，K )],a x (α k , r 0，k ) represents the steering vector formed by the incident signal of the kth narrowband near-field source and the x-axis, [] T Represents the transpose of a matrix or vector, A y represents the array manifold matrix of the uniform linear array located on the y-axis, A y =[a y (β1, r 0，1 ), ..., a y (β k , r 0，k ), ..., a y (β K , r 0，K )],a y (β k , r 0，k) represents the steering vector formed by the incident signal of the kth narrowband near-field source and the y-axis, s(t) represents the baseband signal vector, s(t) = [s1(t), ..., s k (t), ..., s K (t)] T , w x (t) represents the noise vector on the uniform linear array located on the x-axis, w y (t) represents the noise vector on the uniform linear array located on the y-axis,

[0010] Step 3: Using the symmetry of the uniform linear array, construct the fourth-order cumulant of the signal received by the array element at position (-m, 0) and the array element at position (0, -n) after a delay of τ, which is related to m and n and is recorded as C(m, n, τ). The fourth-order cumulant matrix of the received signal after the delay τ is further obtained, denoted as C(τ), And C(τ) is expressed as Then, under the premise that the lengths of x(t) and y(t) are both J=N+L, the maximum likelihood estimation is used for x(t) and y(t), and the delay τ=lT is obtained. s The virtual delayed receiving data after is recorded as r(τ). r(τ) is the maximum likelihood approximation of C(τ); after L delays, τ = T s , 2T s , ..., LT s The sampling of the virtual pseudo-snapshot receiving data matrix based on the delay is obtained, which is recorded as R, R = [vec (r (T s )),vec(r(2T s )), ..., vec(r(LT s ))]; where τ is the delay variable, Indicates finding x m (t) The cumulative amount of x is , the superscript "*" indicates taking the conjugate, -m (t+τ) represents the signal x received by the element at position (-m, 0) in the uniform linear array on the x-axis. -m (t) The signal obtained after a delay of τ, y -n (t+τ) represents the signal y received by the array element at position (0, -n) in the uniform linear array on the y-axis. -n (t) is the signal obtained after a delay of τ, E[] represents the mathematical expectation, Indicates k (t) is the fourth-order cumulant after a delay of τ, Indicates seeking k (t), The cumulative amount, s k (t+τ) represents s k (t) The signal obtained after a delay of τ, r -m，k represents the distance between the kth narrowband near-field source and the center of the array element at position (-m, 0), r -n，k represents the distance between the kth narrowband near-field source and the center of the array element at position (0, -n), and are the intermediate quantities introduced, C s (τ) represents the fourth-order cumulant matrix of the baseband signal after delay τ, “⊙” represents Khatri-Rao product, N represents the number of snapshots, L represents the number of pseudo snapshots, l=1, 2, ..., L, T s represents the sampling time, Indicates M x The identity matrix of order, Indicates M y The unit matrix of order, the superscript "H" indicates the conjugate transpose, x(t+lT s ) represents the delay of x(t) lT s The signal vector obtained after that is y(t+lT s ) represents the delay of y(t) lT s The signal vector obtained after that, vec(·) represents the vectorization operation of the matrix;

[0011] Step 4: Calculate the covariance matrix Z of R, Z = E(RR H )=RR H / L; then perform eigendecomposition on Z to obtain Among them, U signal represents the signal subspace of Z. All eigenvalues ​​obtained by eigendecomposition of Z are sorted from large to small. signal It is composed of the eigenvectors corresponding to the first K largest eigenvalues, ∑ signal is the signal diagonal matrix, U noise represents the noise subspace of Z, U noise It consists of the eigenvectors corresponding to all eigenvalues ​​except the first K largest eigenvalues, ∑ noise is a noise diagonal matrix;

[0012] Step 5: According to γ -m，k =γ m，kThe symmetry of the k-th narrowband near-field source and the steering vector a formed by the x-axis x (α k , r 0，k ) is expressed as: Similarly, according to γ -n，k =γ n，k The symmetry of the signal incident from the kth narrowband near-field source and the steering vector a formed by the y-axis y (β k , r 0，k ) is expressed as:

[0013]

[0014] Step 6: Let α, β, and r represent the vector variable composed of the angles between the signals incident from the K narrow-band near-field sources and the x-axis, the vector variable composed of the angles between the signals incident from the K narrow-band near-field sources and the y-axis, and the vector variable composed of the distances between the centers of the K narrow-band near-field sources and the reference array element; according to the reduced-rank MUSIC algorithm, construct the spectrum peak searchers corresponding to β, α, and r respectively, and perform a one-dimensional search on the spectrum peak searchers to obtain the estimated values ​​of β, α, and r respectively. The specific process is as follows:

[0015] Step 6.1: Obtain an estimate of β: Let Then Decompose it and get:

[0016] Then, according to the reduced-rank MUSIC algorithm, the spectrum peak searcher corresponding to β is constructed, which is expressed as: Then perform a one-dimensional search on the spectrum peak searcher corresponding to β, and select the angles corresponding to the largest K spectrum peaks as the estimated value of β, recorded as in, represents the Kronecker product, a y (β k , r 0，k ) is replaced by β, r to obtain a y (β, r), a x (α k , r 0，k ) is replaced by α, r to get a x (α, r), ξ y (β k ) is replaced by β to obtain ξ y (β),υ y (β k , r 0，k) is replaced by β, r to obtain υ y (β, r), ξ x (α k ) is replaced by α to obtain ξ x (α),υ x (α k , r 0，k ) is replaced by α, r to obtain υ x (α, r), P(β) represents the spectrum peak searcher corresponding to β, and det() represents taking the determinant;

[0017] Step 6.2: Get an estimate of α: Further decomposition gives: ; Then according to the reduced rank MUSIC algorithm, combined with the obtained step 6.1 Construct α k , k = 1, 2, ..., K corresponding to the spectrum peak searcher, expressed as: Then for α k , k=1, 2, ..., K corresponding to the spectrum peak searcher to perform a one-dimensional search, select the angle corresponding to the largest spectrum peak as α k , k = 1, 2, ..., the estimated value of K is denoted as That is, we get the estimated value of α, which can be expressed as in, P(α k ) represents α k The corresponding peak searcher, the parameter in G(α) is replaced by α k Get G(α k ), the parameters in B(β) are replaced by get

[0018] Step 6.3: Obtain an estimated value of r: Based on the reduced-rank MUSIC algorithm, combined with the estimated value obtained in step 6.1 and obtained in step 6.2 Construct 0,k , k = 1, 2, ..., K corresponding to the spectrum peak searcher, expressed as: ; then to r 0,k , k=1, 2, ..., K corresponding to the spectrum peak searcher to perform a one-dimensional search, select the distance value corresponding to the largest spectrum peak as r 0,k , k = 1, 2, ..., the estimated value of K is denoted as That is, we get the estimated value of r, which is recorded as Among them, P(r 0,k ) represents r 0，k The corresponding peak searcher, The parameters in are replaced by r 0，k get

[0019] Step 7: According to cos(α k )=sin(θ k )cos(ψ k ) and sin(β k )=sin(θ k )sin(ψ k ),as well as and Calculate θ1, ..., θ k , ..., θ K Their respective estimated values ​​are denoted as At the same time, ψ1, ..., ψ k , ..., ψ K Their respective estimated values ​​are denoted as Then we can get the estimated values ​​of the three-dimensional parameters of the K narrow-band near-field sources, which are:

[0020] Compared with the prior art, the advantages of the present invention are:

[0021] 1) The method of the present invention uses a cross array to receive signals and processes phase shifts through Fresnel approximation; considering the attenuation of signal amplitude, the signal received by the cross array is combined to construct a fourth-order cumulant matrix, and the three-dimensional parameter estimation problem is decomposed into a one-dimensional parameter estimation problem through step-by-step decoupling. First, the pitch angle and azimuth angle of the incident signal are decomposed into angles with the x-axis and the y-axis, and then the angle spectrum peak searcher constructed by the reduced-rank MUSIC algorithm is used for estimation; then, a distance spectrum peak searcher is constructed to estimate the distance; finally, the pitch angle and azimuth angle are calculated according to the estimated angle between the incident signal and the x-axis and the y-axis to achieve three-dimensional parameter estimation. The method of the present invention only requires multiple one-dimensional searches, has low computational complexity, and has good parameter estimation performance.

[0022] 2) Through a large number of Monte Carlo simulations, it is verified that the method of the present invention can effectively estimate the pitch angle and azimuth angle of the signal incident from the near-field source, as well as the distance between the near-field source and the center of the reference array element, and that the two-dimensional angle and one-dimensional distance parameters can be automatically matched.

[0023] 3) Compared with the existing near-field source parameter estimation methods, the fourth-order cumulant construction adopted by the method of the present invention can effectively suppress the influence of Gaussian white noise and has better estimation performance. BRIEF DESCRIPTION OF THE DRAWINGS

[0024] Figure 1 A schematic diagram of a cross array established in the method of the present invention and a narrow-band near-field source incident on the cross array;

[0025] Figure 2 It is a scatter point matching result diagram of the two-dimensional angle parameters and one-dimensional distance parameters of three narrow-band near-field sources estimated by the method of the present invention in the simulation experiment. DETAILED DESCRIPTION

[0026] The present invention is further described in detail below with reference to the accompanying drawings.

[0027] The present invention relates to a method for estimating three-dimensional parameters of space-time near-field sources based on step-by-step decoupling of the cumulant domain. The method receives a signal incident from a near-field source on a cross array, and jointly constructs a fourth-order cumulant matrix for the signals received by the array elements in the uniform linear array located on the x-axis and the signals received by the array elements in the uniform linear array located on the y-axis, and considers the influence of signal amplitude attenuation; subsequently, the joint steering vector formed by the steering vector formed by the signal incident from the narrowband near-field source and the x-axis and the steering vector formed by the signal incident from the narrowband near-field source and the y-axis are decomposed, and a spectrum peak searcher is constructed using a reduced-rank MUSIC algorithm to obtain corresponding angle and distance estimates. When performing DOA estimation, the method of the present invention first decomposes the joint steering vector into a matrix only related to β and a matrix simultaneously related to α, β, and r, and constructs a spectrum peak searcher only related to β through the reduced-rank MUSIC algorithm, performs a spectrum peak search, and obtains an estimated value of β, which is recorded as Then the joint steering vector is decomposed into a matrix related only to α and β and a matrix related to α, β, and r at the same time. Through the reduced-rank MUSIC algorithm, a spectrum peak searcher related to α and β is constructed and brought into Perform a peak search and obtain an estimated value of α, denoted as Then, according to the MUSIC algorithm, a spectrum peak searcher related to α, β, and r is constructed, and the estimated and Substitute them in turn to get the estimated value of the distance, recorded as Finally, and Decomposed into estimates of the elevation angle θ and the azimuth angle ψ.

[0028] The present invention proposes a method for estimating three-dimensional parameters of space-time near-field sources based on step-by-step decoupling of cumulant domain, which comprises the following steps:

[0029] Step 1: Create a cross array in a three-dimensional coordinate system: Figure 1 As shown, the cross array is used as a signal receiving array, which includes a uniform linear array located on the x-axis of the three-dimensional coordinate system and a uniform linear array located on the y-axis of the three-dimensional coordinate system. The uniform linear array located on the x-axis is composed of M x The uniform linear array on the y-axis consists of M array elements evenly distributed. yThe two uniform linear arrays share the central element at the origin, and the central element is used as the reference element of the cross array. Each uniform linear array is symmetrical with the reference element as the center, and the distance between two adjacent elements in each uniform linear array is d; where M x Represents the total number of uniform linear array elements located on the x-axis, M x =2N x +1, N x represents the number of single-sided array elements of the uniform linear array located on the x-axis, M y Represents the total number of uniform linear array elements located on the y-axis, M y =2N y +1, N y represents the number of single-sided array elements of the uniform linear array located on the y-axis. In this embodiment, N x =N y =5, the total number of elements in the cross array is M, M = M x +M y -1, d=λ / 8, λ represents the signal wavelength, in this embodiment, λ=1, and the array element is a sensor.

[0030] Step 2: Assume that there are K space-time independent narrowband near-field sources that simultaneously incident signals on each element in the cross array, and each element receives all the signals that arrive at it; the three-dimensional parameters of the K space-time independent narrowband near-field sources are expressed as (θ1, ψ1, r 0，1 ),…,(θ k , ψ k , r 0，k ),…,(θ K , ψ K , r 0，K ), the signal received by the array element at position (m, 0) in the uniform linear array on the x-axis is recorded as x m (t), The signal received by the array element at position (0, n) in the uniform linear array on the y-axis is denoted as y n (t), Among them, N x ≥K and N y ≥K,K≥1,1≤k≤K,θ k represents the elevation angle of the signal incident from the kth narrowband near-field source, ψ k represents the azimuth of the signal incident from the kth narrowband near-field source, r 0，k represents the distance between the kth narrowband near-field source and the center of the reference array element, m∈{-N x , -N x +1,…,Nx -1, N x}, n∈{-N y , -N y +1,…,N y -1, N y}, t is the time variable, r m，k represents the distance between the kth narrowband near-field source and the center of the array element at position (m, 0), α k represents the angle between the signal incident from the kth narrowband near-field source and the x-axis, α k With θ k , k The relationship is cos(α k )=sin(θ k )cos(ψ k ), r n，k represents the distance between the kth narrowband near-field source and the center of the array element at position (0, n), β k represents the angle between the signal incident from the kth narrowband near-field source and the y-axis, β k With θ k , k The relationship is sin(β k )=sin(θ k )sin(ψ k ), p represents the a priori known signal amplitude attenuation exponent value, e represents the natural constant, e = 2.71828, j represents the imaginary unit, s k (t) represents the baseband signal received and demodulated by any array element from the signal incident by the kth narrowband near-field source, w m (t) represents the additive white Gaussian noise on the array element at position (m, 0), w n (t) represents the additive Gaussian white noise on the array element at position (0, n), γ m，k and γ n，k are the intermediate quantities introduced,

[0031] Then the phase shift r in formula (1) m，k -r 0，k Using the Fresnel approximation, formula (1) can be rewritten as: Similarly, for the phase shift r in formula (2), n，k -r 0，k Using the Fresnel approximation, formula (2) can be rewritten as: in,

[0032]

[0033] Based on formula (3), the received signal vector of all array elements in the uniform linear array on the x-axis is recorded as x(t), x(t) = A x s(t)+w x (t); Similarly, based on formula (4), the received signal vectors of all array elements in the uniform linear array located on the y-axis are recorded as y(t), y(t) = A y s(t)+w y (t); where A x represents the array manifold matrix of the uniform linear array located on the x-axis, A x =[a x (α1, r 0，1 ), ..., a x (α k , r 0，k ), ..., a x (α K , r 0，K )],a x (α k , r 0，k ) represents the steering vector formed by the incident signal of the kth narrowband near-field source and the x-axis, []T represents the transpose of a matrix or vector, A y represents the array manifold matrix of the uniform linear array located on the y-axis, A y =[a y (β1, r 0，1 ), ..., a y (β k , r 0，k ), ..., a y (β K , r 0，K )],a y (β k , r 0，k ) represents the steering vector formed by the incident signal of the kth narrowband near-field source and the y-axis, s(t) represents the baseband signal vector, s(t) = [s1(t), ..., s k (t), ..., s K (t)] T , w x (t) represents the noise vector on the uniform linear array located on the x-axis, w y (t) represents the noise vector on the uniform linear array located on the y-axis,

[0034] Step 3: Using the symmetry of the uniform linear array, construct the fourth-order cumulant of the signal received by the array element at position (-m, 0) and the array element at position (o, -n) related to m and n after a delay of τ, denoted as C(m, n, τ). The fourth-order cumulant matrix of the received signal after the delay τ is further obtained, denoted as C(τ), And C(τ) is expressed as Then, under the premise that the lengths of x(t) and y(t) are both J=N+L, the maximum likelihood estimation is used for x(t) and y(t), and the delay τ=lT is obtained. s The virtual delayed receiving data after is recorded as r(τ).

[0035] r(τ) is the maximum likelihood approximation of C(τ); after L delays, τ = T s , 2T s , ..., LT s The sampling of the virtual pseudo-snapshot receiving data matrix based on the delay is obtained, which is recorded as R, R = [vec (r (T s )),vec(r(2T s )), ..., vec(r(LT s ))]; where τ is the delay variable, Indicates finding x m (t), The cumulative amount of x is , the superscript "*" indicates taking the conjugate, -m (t+τ) represents the signal x received by the element at position (-m, 0) in the uniform linear array on the x-axis. -m (t) The signal obtained after a delay of τ, y -n (t+τ) represents the signal y received by the array element at position (0, -n) in the uniform linear array on the y-axis. -n (t) is the signal obtained after a delay of τ, E[] represents the mathematical expectation, Indicates k (t) is the fourth-order cumulant after a delay of τ, Indicates seeking k (t), The cumulative amount, s k (t+τ) represents s k (t) The signal obtained after a delay of τ, r -m，k represents the distance between the kth narrowband near-field source and the center of the array element at position (-m, 0), r -n，krepresents the distance between the kth narrowband near-field source and the center of the array element at position (0, -n), and are the intermediate quantities introduced, C s (τ) represents the fourth-order cumulant matrix of the baseband signal after delay τ, “⊙” represents Khatri-Rao product, N represents the number of snapshots, L represents the number of pseudo snapshots, l=1, 2, ..., L, T s represents the sampling time, Indicates M x The identity matrix of order, Indicates M y The unit matrix of order, the superscript "H" indicates the conjugate transpose, x(t+lT s ) represents the delay of x(t) lT s The signal vector obtained after that is y(t+lT s ) represents the delay of y(t) lT s The signal vector obtained after , vec(·) represents the vectorization operation of the matrix.

[0036] Step 4: Calculate the covariance matrix Z of R, Z = E(RR H )=RR H / L; then perform eigendecomposition on Z to obtain Among them, U signal represents the signal subspace of Z. All eigenvalues ​​obtained by eigendecomposition of Z are sorted from large to small. signal It is composed of the eigenvectors corresponding to the first K largest eigenvalues, ∑ signal is the signal diagonal matrix, U noise represents the noise subspace of Z, U noise It consists of the eigenvectors corresponding to all eigenvalues ​​except the first K largest eigenvalues, ∑ noise is a diagonal noise matrix.

[0037] Step 5: According to γ -m，k =γ m,k The symmetry of the k-th narrowband near-field source and the steering vector a formed by the x-axis x (α k , r 0,k ) is expressed as: Similarly, according to γ -n，k =γ n，k The symmetry of the signal incident from the kth narrowband near-field source and the steering vector a formed by the y-axis y (β k, r 0，k ) is expressed as: Among them, γ -m，k according to Get, γ -n，k according to get.

[0038] Step 6: Let α, β, and r represent the vector variable composed of the angles between the signals incident from the K narrow-band near-field sources and the x-axis, the vector variable composed of the angles between the signals incident from the K narrow-band near-field sources and the y-axis, and the vector variable composed of the distances between the centers of the K narrow-band near-field sources and the reference array element; according to the reduced-rank MUSIC algorithm, construct the spectrum peak searchers corresponding to β, α, and r respectively, and perform a one-dimensional search on the spectrum peak searchers to obtain the estimated values ​​of β, α, and r respectively. The specific process is as follows:

[0039] Step 6.1: Obtain an estimate of β: Let Then Decompose it and get:

[0040] Then, according to the reduced-rank MUSIC algorithm, the spectrum peak searcher corresponding to β is constructed, which is expressed as: Then perform a one-dimensional search on the spectrum peak searcher corresponding to β, and select the angles corresponding to the largest K spectrum peaks as the estimated value of β, recorded as in, represents the Kronecker product, a y (β k , r 0，k ) is replaced by β, r to obtain a y (β, r), a x (α k , r 0，k ) is replaced by α, r to get a x (α, r), ξ y (β k ) is replaced by β to obtain ξ y (β),υ y (β k , r 0，k ) is replaced by β, r to obtain υ y (β, r), ξ x (α k ) is replaced by α to obtain ξ x (α),υ x (α k , r0，k ) is replaced by α, r to obtain υ x (α, r), P(β) represents the spectrum peak searcher corresponding to β, and det() represents taking the determinant.

[0041] Step 6.2: Get an estimate of α: Further decomposition gives: Then, according to the reduced-rank MUSIC algorithm, combined with the obtained in step 6.1 Construct α k , k = 1, 2, ..., K corresponding to the spectrum peak searcher, expressed as: ; Then for α k , k=1, 2, ..., K corresponding to the spectrum peak searcher to perform a one-dimensional search, select the angle corresponding to the largest spectrum peak as α k , k = 1, 2, ..., the estimated value of K is denoted as That is, we get the estimated value of α, which can be expressed as in, P(α k ) represents α k The corresponding peak searcher, the parameter in G(α) is replaced by α k Get G(α k ), the parameters in B(β) are replaced by get

[0042] Here, it is equivalent to using the reduced-rank MUSIC algorithm and combining it with the one obtained in step 6.1 Construct the peak searcher corresponding to α1 Perform a one-dimensional search on the spectrum peak searcher corresponding to α1, and select the angle corresponding to the largest spectrum peak as the estimated value of α1, which is recorded as Similarly, according to the reduced-rank MUSIC algorithm, combined with the obtained in step 6.1 Construct the spectrum peak searcher corresponding to α2 Perform a one-dimensional search on the spectrum peak searcher corresponding to α2, and select the angle corresponding to the largest spectrum peak as the estimated value of α2, denoted as By analogy, after constructing K peak searchers and performing one-dimensional searches, a total of K angles are obtained as estimated values ​​of α, denoted as

[0043] Step 6.3: Obtain an estimated value of r: Based on the reduced-rank MUSIC algorithm, combined with the estimated value obtained in step 6.1 and obtained in step 6.2 structure The corresponding peak searcher is expressed as: ; then to r 0,k, k=1, 2, ..., K corresponding to the spectrum peak searcher to perform a one-dimensional search, select the distance value corresponding to the largest spectrum peak as r 0,k , k = 1, 2, ..., the estimated value of K is denoted as That is, we get the estimated value of r, which is recorded as Among them, P(r 0，k ) represents r 0，k The corresponding peak searcher, The parameters in are replaced by r 0，k get

[0044] Here, it is equivalent to using the reduced-rank MUSIC algorithm and combining it with the one obtained in step 6.1 and obtained in step 6.2 Construct 0，1 Corresponding peak searcher R 0，1 The corresponding spectrum peak searcher performs a one-dimensional search and selects the distance value corresponding to the largest spectrum peak as r 0，1 The estimated value of Similarly, according to the reduced-rank MUSIC algorithm, combined with the obtained in step 6.1 and obtained in step 6.2 Construct 0，2 Corresponding peak searcher R 0，2 The corresponding spectrum peak searcher performs a one-dimensional search and selects the distance value corresponding to the largest spectrum peak as r 0，2 The estimated value of By analogy, after constructing K peak searchers and performing one-dimensional searches, a total of K distance values ​​are obtained as the estimated values ​​of r, denoted as

[0045] Step 7: According to cos(α k )=sin(θ k )cos(ψ k ) and sin(β k )=sin(θ k )sin(ψ k ),as well as and Calculate θ1, ..., θ k , ..., θ K Their respective estimated values ​​are denoted as At the same time, ψ1, ..., ψ k , ..., ψ K Their respective estimated values ​​are denoted as Then we can get the estimated values ​​of the three-dimensional parameters of the K narrow-band near-field sources, which are:

[0046] In order to verify the effectiveness and accuracy of the method of the present invention, a simulation test is carried out on the method of the present invention, as follows:

[0047] In the simulation, the spacing between two adjacent array elements in each uniform linear array is set to d = λ / 8, and λ is taken as 1. It is assumed that there are three spatially and temporally uncorrelated narrowband near-field sources in the spatial domain, and the signals they emit are incident on the cross array. The actual values ​​of the three-dimensional parameters of the three narrowband near-field sources are (40°, 155°, λ), (120°, 45°, 1.5λ) and (135°, 75°, 2λ), respectively. It is assumed that the signal amplitude attenuation of the signal incident from the narrowband near-field source to the array element, i.e., the sensor, follows the physical inverse square law, and the signal amplitude attenuation exponent value p is taken as 2. The noise on the array element is additive white Gaussian noise. N x =5,N y =5, the Monte Carlo number is 1000. Set the signal-to-noise ratio SNR to 20dB, the number of snapshots N to 800, and the number of pseudo snapshots L to 100.

[0048] The method of the present invention is used to estimate the three-dimensional parameters of each of the three narrowband near-field sources, namely the pitch angle of the signal incident from each narrowband near-field source, the azimuth angle of the signal incident from each narrowband near-field source, and the distance between each narrowband near-field source and the center of the reference array element. Figure 2 The scatter point matching result diagram of the two-dimensional angle parameters and one-dimensional distance parameters of the three narrow-band near-field sources estimated by the method of the present invention is given. Figure 2 It can be seen from the figure that the method of the present invention can realize the automatic matching of angle parameters (the pitch angle and azimuth angle of the signal incident by the narrowband near-field source) and distance parameters (the distance between the narrowband near-field source and the center of the reference array element), and ensure a certain accuracy, with the angle error not exceeding 1° and the distance error not exceeding 0.1λ. Therefore, it is fully demonstrated that the method of the present invention can realize the three-dimensional parameter estimation of the near-field source, and can obtain good parameter estimation performance under a certain number of snapshots and signal-to-noise ratio.

Claims

1. A three-dimensional parameter estimation method for space-time near-field sources based on step-by-step decoupling in the cumulant domain. Features The steps include: Step 1: Establish a cross array in a three-dimensional coordinate system: The cross array is used as a signal receiving array, which includes a uniform linear array located on the x-axis of the three-dimensional coordinate system and a uniform linear array located on the y-axis of the three-dimensional coordinate system. The uniform linear array located on the x-axis is composed of M x The uniform linear array on the y-axis consists of M array elements evenly distributed. y The two uniform linear arrays share the central element at the origin, and the central element is used as the reference element of the cross array. Each uniform linear array is symmetrical with the reference element as the center, and the distance between two adjacent elements in each uniform linear array is d; where M x represents the total number of uniform linear array elements located on the x-axis, M x =2N x +1, N x represents the number of single-sided array elements of the uniform linear array located on the x-axis, M y Represents the total number of uniform linear array elements located on the y-axis, M y =2N y +1, N y It represents the number of single-sided array elements of the uniform linear array on the y-axis. The total number of array elements of the cross array is M, M = M x +M y -1, d = λ / 8, λ represents the signal wavelength, and the array element is the sensor; Step 2: Assume that there are K space-time independent narrowband near-field sources that simultaneously incident signals on each element in the cross array, and each element receives all the signals that arrive at it; the three-dimensional parameters of each of the K space-time independent narrowband near-field sources are expressed as (θ 1 , ψ 1 , r 0，1 ),…,(θ k , ψ k , r 0，k ),…,(θ K , ψ K , r 0，K ), the signal received by the element at position (m, 0) in the uniform linear array on the x-axis is recorded as x m (t), The signal received by the array element at position (0, n) in the uniform linear array on the y-axis is denoted as y n (t), Among them, N x ≥K and N y ≥K,K≥1,1≤k≤K,θ k represents the elevation angle of the signal incident from the kth narrowband near-field source, ψ k represents the azimuth of the signal incident from the kth narrowband near-field source, r 0，k represents the distance between the kth narrowband near-field source and the center of the reference array element, m∈{-N x , -N x +1,…,N x -1, N x }, n∈{-N y , -N y +1,…,N y -1, N y }, t is the time variable, r m，k represents the distance between the kth narrowband near-field source and the center of the array element at position (m, 0), α k represents the angle between the signal incident from the kth narrowband near-field source and the x-axis, α k With θ k , k The relationship is cos(α k )=sin(θ k )cos(ψ k ), r n，k represents the distance between the kth narrowband near-field source and the center of the array element at position (0, n), βk represents the angle between the signal incident from the kth narrowband near-field source and the y-axis, β k With θ k , k The relationship is sin(β k )=sin(θ k )sin(ψ k ), p represents the a priori known signal amplitude attenuation exponent value, e represents the natural constant, j represents the imaginary unit, s k (t) represents the baseband signal received and demodulated by any array element from the signal incident by the kth narrowband near-field source, w m (t) represents the additive white Gaussian noise on the array element at position (m, 0), w n (t) represents the additive Gaussian white noise on the array element at position (0, n), γ m，k and γ n，k are the intermediate quantities introduced, Then the phase shift r in formula (1) m，k -r 0，k Using the Fresnel approximation, formula (1) can be rewritten as: Similarly, for the phase shift r in formula (2), n，k -r 0，k Using the Fresnel approximation, formula (2) can be rewritten as: in, Based on formula (3), the received signal vector of all array elements in the uniform linear array on the x-axis is recorded as x(t), x(t) = A x s(t)+w x (t); Similarly, based on formula (4), the received signal vector of all array elements in the uniform linear array located on the y-axis is recorded as y(t), y(t) = A y s(t)+w y (t); where A x represents the array manifold matrix of the uniform linear array located on the x-axis, A x =[a x (α 1 , r 0，1 ), ..., a x (α k , r 0，k ), ..., a x (α K , r 0，K )],a x (α k , r 0，k ) represents the steering vector formed by the incident signal of the kth narrowband near-field source and the x-axis, [] T Represents the transpose of a matrix or vector, A y represents the array manifold matrix of the uniform linear array located on the y-axis, A y =[a y (β 1 , r 0，1 ), ..., a y (β k , r 0，k ), ..., a y (β K , r 0，K )],a y (β k , r 0，k ) represents the steering vector formed by the incident signal of the kth narrowband near-field source and the y-axis, s(t) represents the baseband signal vector, s(t) = [s 1 (t), ..., s k (t), ..., s K (t)] T , w x (t) represents the noise vector on the uniform linear array on the x-axis, w y (t) represents the noise vector on the uniform linear array located on the y-axis, Step 3: Using the symmetry of the uniform linear array, construct the fourth-order cumulant of the signal received by the array element at position (-m, 0) and the array element at position (0, -n) after a delay of τ, which is related to m and n and is recorded as C(m, n, τ). The fourth-order cumulant matrix of the received signal after the delay τ is further obtained, denoted as C(τ), And C(τ) is expressed as Then, under the premise that the lengths of x(t) and y(t) are both J=N+L, the maximum likelihood estimation is used for x(t) and y(t), and the delay τ=lT is obtained. s The virtual delayed receiving data after is recorded as r(τ). r(τ) is the maximum likelihood approximation of C(τ); after L delays, τ = T s , 2T s , ..., LT s The sampling of the virtual pseudo-snapshot receiving data matrix based on the delay is obtained, which is recorded as R, R = [vec (r (T s )),vec(r(2T s )), ..., vec(r(LT s ))]; where τ is the delay variable, Indicates finding x m (t), Cumulative amount, the superscript "*" indicates taking the conjugate, x -m (t+τ) represents the signal x received by the element at position (-m, 0) in the uniform linear array on the x-axis. -m (t) The signal obtained after a delay of τ, y -n (t+τ) represents the signal y received by the array element at position (0, -n) in the uniform linear array on the y-axis. -n (t) is the signal obtained after a delay of τ, E[] represents the mathematical expectation, Indicates k (t) is the fourth-order cumulant after a delay of τ, , Indicates seeking k (t), The cumulative amount, s k (t+τ) represents s k (t) The signal obtained after a delay of τ, r -m，k represents the distance between the kth narrowband near-field source and the center of the array element at position (-m, 0), r -n，k represents the distance between the kth narrowband near-field source and the center of the array element at position (0, -n), and are the intermediate quantities introduced, C s (τ) represents the fourth-order cumulant matrix of the baseband signal after delay τ, "⊙" represents the Khatri-Rao product, N represents the number of snapshots, L represents the number of pseudo snapshots, l = 1, 2, ..., L, T s represents the sampling time, Indicates M x The identity matrix of order, Indicates M y The unit matrix of order, the superscript "H" indicates the conjugate transpose, x(t+lT s ) represents the delay of x(t) lT s The signal vector obtained after that is y(t+lT s ) represents the delay of y(t) lT s The signal vector obtained after that, vec(·) represents the vectorization operation of the matrix; Step 4: Calculate the covariance matrix Z of R, Z = E(RR H )=RR H / L; then perform eigendecomposition on Z to obtain Among them, U signal represents the signal subspace of Z. All eigenvalues ​​obtained by eigendecomposition of Z are sorted from large to small. signal It is composed of the eigenvectors corresponding to the first K largest eigenvalues, ∑ signal is the signal diagonal matrix, U noise represents the noise subspace of Z, U noise It consists of the eigenvectors corresponding to all eigenvalues ​​except the first K largest eigenvalues, ∑ noise is a noise diagonal matrix; Step 5: According to γ -m，k =γ m,k The symmetry of the k-th narrowband near-field source and the steering vector a formed by the x-axis x (α k , r 0，k ) is expressed as: Similarly, according to γ -n，k =γ n，k The symmetry of the signal incident from the kth narrowband near-field source and the steering vector a formed by the y-axis y (β k , r 0，k ) is expressed as: Step 6: Let α, β, and r represent the vector variable composed of the angles between the signals incident from the K narrow-band near-field sources and the x-axis, the vector variable composed of the angles between the signals incident from the K narrow-band near-field sources and the y-axis, and the vector variable composed of the distances between the centers of the K narrow-band near-field sources and the reference array element; according to the reduced-rank MUSIC algorithm, construct the spectrum peak searchers corresponding to β, α, and r respectively, and perform a one-dimensional search on the spectrum peak searchers to obtain the estimated values ​​of β, α, and r respectively. The specific process is as follows: Step 6.1: Obtain an estimate of β: Let Then Decompose it and get: Then, according to the reduced-rank MUSIC algorithm, the spectrum peak searcher corresponding to β is constructed, which is expressed as: Then perform a one-dimensional search on the spectrum peak searcher corresponding to β, and select the angles corresponding to the largest K spectrum peaks as the estimated value of β, recorded as in, represents the Kronecker product, a y (β k , r 0，k ) is replaced by β, r to obtain a y (β, r), a x (α k , r 0，k ) is replaced by α, r to get a x (α, r), ξ y (β k ) is replaced by β to obtain ξ y (β),υ y (β k , r 0，k ) is replaced by β, r to obtain υ y (β, r), ξ x (α k ) is replaced by α to obtain ξ x (α),υ x (α k , r 0，k ) is replaced by α, r to obtain υ x (α, r), P(β) represents the spectrum peak searcher corresponding to β, and det() represents taking the determinant; Step 6.2: Get an estimate of α: Further decomposition gives: Then, according to the reduced-rank MUSIC algorithm, combined with the obtained in step 6.1 Construct α k , k = 1, 2, ..., K corresponding to the spectrum peak searcher, expressed as: ; Then for α k , k=1, 2, ..., K corresponding to the spectrum peak searcher to perform a one-dimensional search, select the angle corresponding to the largest spectrum peak as α k , k = 1, 2, ..., the estimated value of K is denoted as That is, we get the estimated value of α, which can be expressed as in, P(α k ) represents α k The corresponding peak searcher, the parameter in G(α) is replaced by α k Get G(α k ), the parameters in B(β) are replaced by get Step 6.3: Obtain an estimated value of r: Based on the reduced-rank MUSIC algorithm, combined with the estimated value obtained in step 6.1 and obtained in step 6.2 structure The corresponding peak searcher is expressed as: ; then to r 0,k , k=1, 2, ..., K corresponding to the spectrum peak searcher to perform a one-dimensional search, select the distance value corresponding to the largest spectrum peak as r 0,k , k = 1, 2, ..., the estimated value of K is denoted as That is, we get the estimated value of r, which is recorded as Among them, P(r 0,k ) represents r 0，k The corresponding peak searcher, The parameters in are replaced by get Step 7: According to cos(α k )=sin(θ k )cos(ψ k ) and sin(β k )=sin(θ k )sin(ψ k ),as well as and Calculate θ 1 , ..., θ k , ..., θ K Their respective estimated values ​​are denoted as At the same time, we calculate ψ 1 , ..., ψ k , ..., ψ K Their respective estimated values ​​are denoted as Then we can get the estimated values ​​of the three-dimensional parameters of the K narrow-band near-field sources, which are:

Images