A direct location method for coherent distributed sources based on parallel factors

Abstract

The present invention discloses a method for direct positioning of coherent distributed information sources based on parallel factors, comprising the following steps: S1, constructing a coherent distributed information source positioning model using angle and delay information to obtain the received signal of each array in the frequency domain; S2, extracting useful frequency points to obtain a known inverse Fourier transform matrix, and constructing a parallel factor model; S3, obtaining a position information matrix using a trilinear alternating least squares method; S4, constructing a direct positioning cost function related to the position information matrix, performing a three-dimensional search on the direct positioning cost function, and obtaining the grid point position with the smallest Euclidean distance to the position information matrix as the estimated position of the coherent information source. The present invention is applicable to antenna arrays with arbitrary spacing, and has application prospects in multipath environments such as dense urban areas.

Description

Technical Field

[0001] The present invention relates to the field of wireless positioning technology, and in particular to a coherent distributed source direct positioning method based on parallel factors. Background Art

[0002] In some multipath environments, such as high-precision positioning of mobile communication equipment and "low, slow and small" drone targets in dense urban areas, the signals emitted by these sources will be scattered by a large number of surrounding scatterers, resulting in the signal received by the array antenna being the superposition of these multipath components. In such cases, the source cannot be simply regarded as a point source, but should be established as a distributed source model. Traditional direct positioning technology based on point sources will experience a serious decline in positioning performance when solving positioning problems in this scenario.

[0003] According to different actual situations, the correlation between multipath components of the same distributed source can be different, and it can be divided into partially coherently distributed (PCD) source, coherently distributed (CD) source and incoherently distributed (ICD) source. In 2015, Ofer first proposed a direct positioning method based on CD source, providing a maximum likelihood estimation operator using joint angle and delay information and other suboptimal operators such as least squares estimation. The positioning performance of maximum likelihood estimation can reach the Cromer-Rao lower bound. Later, Lu Zhiyu improved the direct positioning method based on CD source, and used the conjugate characteristics of non-circular signals in the data domain to enhance the positioning performance and increase its degree of freedom. Although Ofer's method has the best performance and is applicable to any array structure, the computational complexity it faces is extremely high, which is impractical for some low-performance hardware processing scenarios and cannot bring real-time positioning effects. As for Lu Zhiyu's method, it only uses angle information and does not deeply explore potential delay information, so the positioning accuracy is not good enough. Summary of the invention

[0004] Purpose of the invention: The purpose of the present invention is to provide a coherent distributed source direct positioning method based on parallel factors, which can reduce the computational complexity while enhancing the positioning performance by utilizing angle and delay information.

[0005] Technical solution: A method for direct location of coherent distributed sources based on parallel factors, comprising the following steps:

[0006] S1, build a coherent distributed source positioning model using angle and delay information to obtain the received signal x of each array in the frequency domain l(k, j);

[0007] S2, extract useful frequency points, obtain a known inverse Fourier transform matrix, and construct a parallel factor model;

[0008] S3, using the trilinear alternating least squares method to obtain the position information matrix

[0009] S4, constructing a direct positioning cost function related to the position information matrix, performing a three-dimensional search on the direct positioning cost function, and obtaining the grid point position with the smallest Euclidean distance to the position information matrix as the estimated position of the coherent signal source.

[0010] Further, in step S1, the obtained l-th array received signal in the coherent distributed signal source scenario is:

[0011]

[0012] Where K is the number of snapshots, J is the number of sampling points; f j =jf s / J is the center frequency of the jth filter channel, f s is the sampling frequency; s(k, j) represents the output signal of the jth filter channel after discretization; n l (k, j) represents the noise received by the jth filter channel; τ l (p) is the delay information; b l is the generalized array steering vector, whose position is x l，m The received signal on the mth array element is expressed as:

[0013]

[0014] Among them, θ l (p) is the central arrival angle of the distributed source measured by the array, λ is the wavelength, σ l (p) is the angular extension and e is a natural constant.

[0015] Further, in step S2, the constructed parallel factor model is as follows:

[0016]

[0017] in, is an ordered set, L is the number of base stations, M is the number of array elements; P is the number of useful frequency points extracted in the frequency domain, is a corresponding known inverse Fourier transform matrix, S is the matrix of all signals; H is the matrix containing angle and delay information, and its specific expression is:

[0018]

[0019]

[0020] Among them, h L (M) is the information vector received by the Mth array element in the Lth base station, Represents the complex number field.

[0021] Further, in step S3, by optimizing the following cost function:

[0022]

[0023] The estimated signal matrix is:

[0024]

[0025] in, is the estimated information matrix of the previous iteration;

[0026] Then, continue to optimize another cost function:

[0027]

[0028] Among them, Z j yes The Gaussian white noise matrix;

[0029] The new information matrix estimate is:

[0030]

[0031] Repeat the above steps until the set convergence threshold is reached, and finally perform descaling and permutation blurring to obtain the final information matrix

[0032] Furthermore, in step S4, according to the structure of the position information matrix, a direct positioning cost function is constructed to measure the similarity with the information matrix:

[0033]

[0034] in, f cp is the frequency corresponding to the pth useful frequency point, τ L1 is the relative delay measured by the Lth base station; The matrix The column vector consisting of the p-th column of ;

[0035] The direct positioning cost function is searched and solved in three dimensions, and the maximum value position is obtained as the source target position.

[0036] Compared with the prior art, the present invention has the following significant effects:

[0037] 1. The present invention uses angle and delay information for direct position estimation, and its positioning performance is better than that of the method of direct positioning using only angle information;

[0038] 2. The present invention introduces a parallel factor model, under which noise can be effectively suppressed; and it greatly reduces the amount of computation required to solve the direct positioning cost function; the present invention significantly reduces the computational complexity while maintaining excellent positioning performance, and has broader prospects in practical applications. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] Figure 1 is a flow chart of the present invention;

[0040] Figure 2 Positioning description graph for distributed information sources;

[0041] Figure 3 It is a logarithmic curve diagram of the computational complexity of the present invention and the comparative algorithm under different numbers of base stations;

[0042] Figure 4 This is a spectrum peak diagram under one experiment of the present invention;

[0043] Figure 5 Root mean square error performance curves of the present invention and the comparison algorithm under different signal-to-noise ratios;

[0044] Figure 6 It is the error cumulative distribution function diagram of the present invention and the comparison algorithm. DETAILED DESCRIPTION

[0045] The present invention is further described in detail below in conjunction with the accompanying drawings and specific implementation methods.

[0046] In order to achieve the best balance between positioning performance and required computational complexity, the present invention proposes a coherent distributed source direct positioning method based on parallel factor (PARAllel FACtor, PARAFAC), the flow chart is as follows: Figure 1 As shown, the steps include:

[0047] Step 1: construct a coherent distributed source positioning model using angle and delay information to obtain the combined received signals of each array in the frequency domain.

[0048] First, refer to Figure 2 The target positioning diagram in the distributed information source scenario shown in FIG. 1 is a schematic diagram of the target positioning in the distributed information source scenario, assuming that the information source position is p = (x p ,y p ) T , the Euclidean distance between the source and the lth base station is dl (p). The source is expanded to a radius of γ p Surrounded by scatterers with circular symmetric random distribution. Assuming that the total observation time of each base station is T0, the noise-free signal received by the lth base station at time t is:

[0049]

[0050] Where Ω is the angular field of view, φ is an angle belonging to the angular field of view Ω; ζ(φ, t-τ l (p); η l ) is called the angular signal density. l (φ) is the steering vector, and its array element position is x l，m The steering vector expression is Where λ is the wavelength. η l =[θ l (p),σ l (p)] is the position vector, θ l (p) is the central arrival angle of the distributed source measured by the array, σ l (p) is the angle spread; τ l (p) represents the delay.

[0051] Assuming that the channel coherence time is greater than the observation interval length, so that the array response remains unchanged, the source should be modeled as a coherent distributed source. In this case, the noise-free signal received by the lth base station at time t is expressed as follows:

[0052]

[0053] Among them, s(t-τ l (p)) is the time delay τ l (p) signal; in addition, it indicates the position x l，m The received signal expression on the mth array element is:

[0054]

[0055] Here, e represents a natural constant.

[0056] Considering the additive white Gaussian noise n l (t), we have:

[0057] x l (t) = b l s(t-τ l (p))+n l (t), 0≤t≤T0 (3)

[0058] Among them, b lis called the generalized array steering vector. In addition, for convenience, [b l (θ l (p))] m Abbreviated as b lm .

[0059] In order to extract the time delay information, the present invention uses a filter bank to convert the time delay information into an independent phase offset term through discretization processing, as shown below:

[0060]

[0061] Where, j=1, ...J; k=1, ...K, K represents the number of snapshots, and J represents the number of sampling points in each snapshot. represents the center frequency of the jth filter channel, f s is the sampling frequency; s(k, j) represents the output signal of the jth filter channel after discretization; n l (k, j) represents the noise received by the jth filter channel.

[0062] Step 2: Collect the frequency domain signal data received from all base stations, extract the useful frequency points, and obtain a known inverse Fourier transform matrix to construct a parallel factor model.

[0063] In the frequency domain, for each array element, J frequency points are collected, and the frequency is f = [f1, ..., f J ] T , the array element output is:

[0064]

[0065] in, represents the complex field; x lm，k =[x lm (k, 1), …, x lm (k,J)] T , x lm (k, J) is the Jth frequency point, x l The mth element of (k, J), that is, the signal received by the mth array element; and n lm，k =[n lm (k, 1), …, n lm (k,J)] T , n lm (k, J) is the noise received by the Jth frequency point and the mth array element; diag() means taking the diagonal elements of the matrix to form a column vector.

[0066] Using an iDFT matrix (inverse Fourier transform matrix) F H, the array element output is expressed in the time domain as x lm (t) = F H x lm，k According to the Nyquist sampling theorem, in the signal spectrum s k There are a large number of zero-value points or points close to zero in the network, so useful frequency points can be extracted. Assume that the frequency corresponding to the useful frequency point is represents the field of real numbers, f cp represents the frequency of the pth useful frequency point, p = 1, ..., P; the number of useful frequency points is P, and obviously P < J. Then:

[0067]

[0068] in, Represent the Fourier coefficients corresponding to the useful frequency points and F H The corresponding columns are reorganized into iDFT matrices; is an information vector containing angle and delay information.

[0069] After collecting K snapshots, we have:

[0070]

[0071] in, and N lm =[n lm (1),…,n lm (K)].

[0072] After combining all array elements of all base stations, the total frequency point is expressed as:

[0073]

[0074] in, is an ordered set, L is the number of base stations, M is the number of array elements; D u The function (·) is to extract the u-th row vector of the matrix and form a diagonal matrix; N u is the noise vector after receiving K snapshots and J frequency points at each array element of each base station. H is an information matrix containing angle and delay components, which is defined as:

[0075]

[0076] Among them, h L (M) is the information vector received by the Mth array element in the Lth base station.

[0077] Under noise-free conditions,

[0078]

[0079] Formula (10) can be regarded as a three-way tensor The u-th forward slice of .

[0080] Similarly, transverse slices can be obtained and horizontal slices They are:

[0081]

[0082]

[0083] Thus, the PARAFAC model (parallel factor model) can be expressed as the sum of three matrix products:

[0084]

[0085] Among them, x j，k，u express The (j, k)th element of i，p represents the element in the i-th row and p-th column of matrix A, where A is S.H.

[0086] Step 3: Use the trilinear alternating least squares method to perform trilinear decomposition on the PARAFAC model (parallel factor model) to obtain a position information matrix containing position information.

[0087] First, optimize the cost function of the total frequency point:

[0088]

[0089] The estimated signal matrix S is:

[0090]

[0091] in, is the estimated information matrix of the previous iteration.

[0092] Then, continue to optimize another cost function:

[0093]

[0094] Among them, Z j yes The Gaussian white noise matrix;

[0095] The new information matrix estimate can be obtained as:

[0096]

[0097] because It is known that no estimation is needed. Repeat the above steps until the set convergence threshold is reached, and finally perform descaling and permutation blurring to obtain the final information matrix In the following form:

[0098]

[0099] in, Represents the delay relative to the reference base station, l=1,…,L.

[0100] Step 4, construct a direct positioning cost function related to the position information matrix, perform a three-dimensional search on the direct positioning cost function, and obtain the grid point position with the smallest Euclidean distance to the position information matrix, which is the estimated position of the coherent source (target estimated position).

[0101] From the information matrix in step 3 It can be seen that there is a one-to-one mapping relationship between the source position p and it, so the following direct positioning cost function can be constructed:

[0102]

[0103] in, Representation Information Matrix The pth column of .

[0104] When the grid point is closer to the source, e(f cp ,p,γ p )and The greater the similarity, the smaller the Euclidean distance between the two. Therefore, it is only necessary to perform a three-dimensional search on the X-axis, Y-axis, and expansion radius to find the highest spectrum peak point in the power spectrum generated by the cost function. The corresponding position of the highest spectrum peak point is the estimated position of the source.

[0105] Figure 3 This is a curve chart obtained after performing computational complexity analysis on the method of the present invention (CDP-DPD) and other comparative algorithms. The algorithms compared in the coherent distributed source scenario are: 1) GS-DPD: a direct positioning algorithm using generalized subspace data fusion; 2) Angle-only AMLE: a direct positioning algorithm using approximate maximum likelihood estimation using only angle information; 3) ESSI-DPD: a direct positioning method using angle and time delay information based on symmetric rotation invariance; 4) JTDA-MUSIC: a direct positioning method using subspace orthogonality of angle and time delay information. Figure 3 It can be seen that the computational complexity of the method of the present invention is much smaller than that of JTDA-MUSIC, and is substantially the same as that of a low-complexity algorithm such as ESSI-DPD that only requires two-dimensional search, which demonstrates the superiority of the present invention in computational efficiency.

[0106] Figure 4 This is the spectrum peak diagram of the present invention when the signal-to-noise ratio is 0 dB. This is the two-dimensional spectrum peak diagram determined only by the position obtained by the present invention by eliminating the extension radius dimension through simple summation. It can be seen that at the position (110, 120) T The maximum power is at point m, which can intuitively illustrate that the present invention can effectively locate the position of the distributed information source.

[0107] Figure 5 This is a performance comparison curve of various algorithms on the signal-to-noise ratio under the root mean square error (RMSE). First, the definition of RMSE is given:

[0108]

[0109] Among them, Mon is the total number of Monte Carlo, is the estimated value of the source position in the month experiment.

[0110] from Figure 5 It can be seen that the performance of JTDA-MUSIC is closest to the Cromerow lower bound, but the complexity of this algorithm is the highest in the complexity analysis; and although the complexity of ESSI-DPD is low, it sacrifices the positioning performance. The method CDP-DPD of the present invention can achieve a better balance between computational complexity and positioning performance. While having the same high computational efficiency as ESSI-DPD, it also has superior positioning performance close to or even exceeding JTDA-MUSIC under high signal-to-noise ratio.

[0111] In addition to the RMSE performance indicator, another indicator that can be used to evaluate positioning performance is the error cumulative distribution function (CDF), which is defined as:

[0112]

[0113] Among them, error represents the set error value, Num(x) is the Boolean function of x, when x is true, the function count is 1, otherwise it is 0.

[0114] Figure 6 It is a comparison curve diagram of the error cumulative distribution function of each algorithm. It can be clearly classified that the present invention and JTDA-MUSIC algorithm have the smallest error under the same CDF. Although the error of the present invention is larger than that of JTDA-MUSIC in most cases, the CDF of CDP-DPD can reach 100% when the error is 0.07m, and its error variance is smaller, and its stability is also the best in actual positioning.

[0115] In summary, from Figure 3 to Figure 6 Analysis shows that the coherent distributed source direct positioning method based on PARAFAC analysis proposed in the present invention can effectively estimate the position of the coherent distributed source, and can achieve the excellent positioning performance of the JTDA-DPD algorithm while maintaining the same low complexity as the ESSI-DPD algorithm, and has a broader application prospect in actual positioning.

[0116] The embodiments of the present invention are described in detail above with reference to the accompanying drawings, but the present invention is not limited to the above embodiments, and various changes can be made within the knowledge scope of ordinary technicians in this field without departing from the purpose of the present invention.

Claims

1. A method for direct location of coherent distributed sources based on parallel factors, characterized in that: The following steps are involved: S1, build a coherent distributed source positioning model using angle and delay information to obtain the received signal of each array in the frequency domain S2, extract useful frequency points, obtain a known inverse Fourier transform matrix, and construct a parallel factor model; S3, using the trilinear alternating least squares method to obtain the position information matrix S4, constructing a direct positioning cost function related to the position information matrix, performing a three-dimensional search on the direct positioning cost function, and obtaining the grid point position with the smallest Euclidean distance to the position information matrix as the estimated position of the coherent distribution source; In step S1, the coherent distributed information source scenario is obtained. The signal received by each base station is: Where K is the number of snapshots, J is the number of sampling points; f j =jf s / J is the center frequency of the jth filter channel, f s is the sampling frequency; s(k, j) represents the output signal of the jth filter channel after discretization processing; represents the noise received by the jth filter channel; It is the time delay information; is the generalized array steering vector, and its position is The received signal on the mth array element is expressed as: in, is the central arrival angle of the distributed source measured by the array, λ is the wavelength, is the angular extension and e is a natural constant.

2. The method for direct location of coherent distributed sources based on parallel factors according to claim 1, characterized in that: In step S2, the constructed parallel factor model is as follows: in, is an ordered set, L is the number of base stations, M is the number of array elements; P is the number of useful frequency points extracted in the frequency domain, is a corresponding known inverse Fourier transform matrix, S is the matrix of all signals; H is a matrix containing angle and delay information, and its specific expression is: Among them, h L (M) is the information vector received by the Mth array element in the Lth base station, represents a complex domain; Indicates the frequency corresponding to the useful frequency point, Represents the real number domain, and the number of useful frequency points is P.

3. The method for direct location of coherent distributed sources based on parallel factors according to claim 2, characterized in that: In step S3, the following cost function is optimized: The estimated signal matrix is: in, is the estimated information matrix of the previous iteration; Then, continue to optimize another cost function: Among them, Z j yes The Gaussian white noise matrix; The new information matrix estimate is: Repeat the above steps until the set convergence threshold is reached, and finally perform descaling and permutation blurring to obtain the final information matrix 4. The method for direct location of coherent distributed sources based on parallel factors according to claim 2, characterized in that: In step S4, according to the structure of the position information matrix, a direct positioning cost function is constructed to measure the similarity with the information matrix: in, f cp is the frequency corresponding to the pth useful frequency point, is the relative delay measured by the Lth base station; The matrix The column vector composed of the pth column of γ p is the expansion radius of the source target; The direct positioning cost function is searched and solved in three dimensions, and the maximum value position is obtained as the source target position.

Images