Data fusion and frequency focusing broadband radiation source direct positioning algorithm

Abstract

The invention discloses a broadband radiation source direct positioning algorithm based on data fusion and frequency focusing, and the algorithm comprises the steps: carrying out the J-point FFT of data received by each observation station; then, under each frequency point, carrying out fusion processing on the received data of all the observation stations to obtain J groups of received data; then, covariance matrixes of the J groups of data are calculated, eigenvalue decomposition is carried out on the covariance matrixes, and signal subspaces are separated out; then, calculating a focusing matrix corresponding to each frequency point to obtain a final covariance matrix; then, carrying out eigenvalue decomposition on the final covariance matrix, and uploading a noise subspace to a central station; thirdly, gridding the positioning area at the central station, and calculating the function value of the position estimation function in each grid coordinate; and finally, searching a plurality of grid coordinates of the radiation source with the maximum function value, and taking the position as a positioning result. Compared with the existing algorithm, the method provided by the invention has higher positioning precision in a direct positioning scene of the broadband radiation source.

Description

Technical Field

[0001] The present invention relates to the field of wireless positioning technology, and in particular to a broadband radiation source direct positioning algorithm of data fusion and frequency focusing. Background Art

[0002] As a key technology in the field of array signal processing, the radiation source localization technology has broad application prospects in the fields of human-computer interaction, traffic noise localization and drone detection. Traditional emitter localization methods usually include two main steps. First, the parameters related to the source position are estimated, including arrival angle, arrival time, arrival time difference, etc. Subsequently, the estimation from the first step is used to determine the location information of the source. Due to the lack of intrinsic relationship between the parameters obtained from independent measurements of different arrays, the performance of these algorithms is suboptimal.

[0003] Existing direct positioning algorithms are generally designed for narrowband radiation sources. However, with the rapid development of signal processing technology, the electromagnetic environment has become more and more complex, and the frequency distribution range of signals has gradually expanded. Compared with narrowband signals, broadband signals contain more information, have stronger anti-interference ability and higher resolution, and there are more and more application scenarios for broadband radiation sources. Therefore, the study of broadband source positioning methods has important theoretical significance and practical application value. In order to directly locate broadband radiation sources, the present invention proposes a broadband radiation source direct positioning algorithm with data fusion and frequency focusing. Summary of the invention

[0004] The technical problem to be solved by the present invention is to provide a broadband radiation source direct positioning algorithm with data fusion and frequency focusing in view of the defects involved in the background technology.

[0005] The present invention adopts the following technical solutions to solve the above technical problems:

[0006] A broadband radiation source direct positioning algorithm based on data fusion and frequency focusing comprises the following steps:

[0007] Step 1) Set up L observation stations with precisely known positions in the area to be located, each equipped with an M-element uniform linear array, with a spacing of d between each element, and their positions are u l =[x l ,y l ], l = 1, 2, ..., L, used to locate K far-field radiation sources;

[0008] The positions of the K far-field radiation sources are respectively p k =[x k ,y k ], k = 1, 2, ..., K;

[0009] Let the far-field radiation source emit independent broadband stationary signals, the starting frequency and cutoff frequency of the signal are f0 and f1 respectively;

[0010] The signal x received by the uniform linear array at the lth observation station at time t is l (t) is:

[0011]

[0012] Among them, ρ l,k is the unknown channel attenuation generated by the lth observation station in the process of receiving the kth radiation source, n l (t) is the signal additive noise received by the lth observation station at time t, s k (t) is the signal emitted by the kth far-field radiation source at time t, α l (p k ) represents the array steering vector when the lth observation station receives the kth far-field radiation source, τ l (p k ) is the time delay of the signal received by the lth observation station compared to the time when the kth signal is received by the first observation station;

[0013] Step 2), divide the received signal of each observation station in the Ts time period into N sub-time periods, perform a fast Fourier transform of J points on each sub-time period, where J is a preset power of 2, and obtain the frequency domain model x of the signal received by the lth observation station at the jth frequency point in the nth sub-time period l,n (f j ) is as follows, n is an integer greater than or equal to 1 and less than or equal to N, and j is an integer greater than or equal to 1 and less than or equal to J:

[0014]

[0015] In the formula, is the signal frequency of the jth frequency point;

[0016] in, c is the speed of light;

[0017] s k,n (f j ) is s k The Fourier transform result of the nth time segment of (t); l,n (f j ) is n l The Fourier transform result of the nth time segment of (t);

[0018] Step 3), for each observation station, the signal of the jth frequency point in all time periods is expressed as the following matrix:

[0019] X l,n (f j )=A l (f j )Φ l S n (f j )+N l (f j )

[0020] Among them A l (f j )=[α l (f j ,p1),α l (f j ,p2),…,α l (f j ,p K )] T is an M×K dimensional matrix, which is the steering vector of the uniform linear array of the lth observation station at the jth frequency point; Φ l =diag(ρ l,1 ,…,ρ l,K );S n (f j )=[s 1,n (f j ),…s K,n (f j )] T ; N l (f j )=[n l,1 (f j ), …, n l,N (f j )];

[0021] Step 4) Combine the signals received by all observation stations at the jth frequency point in N time periods to obtain the joint reception data model X n (f j )as follows:

[0022]

[0023] In the formula,

[0024] Step 5), calculate the sampling covariance matrix of the joint reception data model right Perform eigenvalue decomposition to obtain the matrix U(f j ); Take the matrix U(f j ) to obtain the signal subspace U at the jth frequency point s (f j );

[0025] Step 6), select the signal subspace U corresponding to the starting frequency f0 in the signal subspace under J frequency points s (f0), let the focusing error matrix of the jth frequency point be For C(f j ) performs singular value decomposition to obtain the matrix U composed of left singular vectors c (f j ) and the matrix V consisting of the right singular vectors c (f j );

[0026] Then the focusing matrix T(f j )=U c (f j )V c (f j ) H , the covariance matrix after focusing corresponding to the jth frequency point

[0027] Step 7) Average the focused covariance matrices corresponding to each frequency point to obtain the final covariance matrix

[0028] right Perform eigenvalue decomposition to obtain the matrix U composed of its eigenvectors, and take out the K+1th column to the last column of the matrix U as the final noise subspace matrix U n ;

[0029] Step 8), dividing the area to be located according to a preset grid side length threshold to form a number of grids, and setting the coordinates of the center point of each grid as the coordinates of the grid;

[0030] Step 9), construct the cost function about the grid coordinate p as follows:

[0031]

[0032] In the formula,

[0033] because γ=[ρ1,…,ρ L ] T is the unknown channel attenuation generated by the L observation stations in the process of receiving the radiation source;

[0034] The cost function is converted to:

[0035]

[0036] where e1=[1,0,...,0]T is an L×1 dimensional matrix;

[0037] Step 10), let Solving the converted cost function using the Lagrange multiplier method yields:

[0038]

[0039] Substitute the solved γ into the cost function The radiation source position estimation function is obtained

[0040] Step 11), input the coordinates of each grid into the position estimation function in turn The estimated value corresponding to each grid is obtained, and the grids corresponding to the largest K estimated values ​​are the positions of the K far-field radiation sources to be determined.

[0041] Compared with the prior art, the present invention adopts the above technical solution and has the following technical effects:

[0042] In the scenario of locating broadband radiation sources, compared with the traditional ISM-DPD method (not using focusing processing) and FSS-DPD method (not fusing all observation station data and losing some information), the method of the present invention significantly improves the positioning accuracy while ensuring that the calculation complexity is not increased. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0044] Figure 2 Positioning scene graphs for distributed arrays;

[0045] Figure 3 It is the two-dimensional spatial spectrum of the method of the present invention;

[0046] Figure 4 A scatter plot of the results of locating multiple radiation sources by the method of the present invention;

[0047] Figure 5 Schematic diagram of root mean square error performance of positioning results of the present invention and the traditional positioning method under different signal-to-noise ratios;

[0048] Figure 6 Schematic diagram of root mean square error performance of positioning results of the present invention and the traditional positioning method under different snapshot numbers. DETAILED DESCRIPTION

[0049] The technical solution of the present invention is further described in detail below in conjunction with the accompanying drawings:

[0050] The present invention can be implemented in many different forms and should not be considered to be limited to the embodiments described herein. On the contrary, these embodiments are provided to make this disclosure thorough and complete, and will fully express the scope of the present invention to those skilled in the art. In the accompanying drawings, components are enlarged for clarity.

[0051] Symbols: In this paper, bold uppercase letters, bold lowercase letters, and italic letters, such as A, a, and a, represent matrices, vectors, and scalars, respectively. T , (·) H and(·) -1 denote the matrix transpose, conjugate transpose and inverse operations, E(·), ||·|| F , tr(·), Re(·), diag(·) represent the expectation, F-norm, trace, taking the real part, and converting the vector to a diagonal matrix, respectively, with triangular superscript letters represents the estimated value of R.

[0052] The present invention provides a broadband source weighted direct positioning method based on subspace focusing, such as Figure 1 As shown, the specific steps include:

[0053] Step 1), consider the location area as follows Figure 2 Positioning scenario. Set up L observation stations with precisely known positions, each equipped with an M-element uniform linear array, with a spacing of d between each element, and their positions are u l =[x l ,y l ], l = 1, 2, ..., L, used to locate K far-field radiation sources. The positions of the K far-field radiation sources are p k =[x k ,y k ], k = 1, 2, ..., K; let the far-field radiation source emit independent broadband stationary signals, the starting frequency and cutoff frequency of the signal are f0 and f1 respectively;

[0054] The signal x received by the uniform linear array at the lth observation station at time t is l (t) is:

[0055]

[0056] Among them, ρ l,k is the unknown channel attenuation generated by the lth observation station in the process of receiving the kth radiation source, n l (t) is the signal additive noise received by the lth observation station at time t, s k (t) is the signal emitted by the kth far-field radiation source at time t, α l (p k0 represents the array steering vector when the lth observation station receives the kth far-field radiation source, τ l (p k ) is the time delay of the signal received by the lth observation station compared to the time when the kth signal is received by the first observation station;

[0057] Step 2), divide the received signal of each observation station in the Ts time period into N sub-time periods, perform a fast Fourier transform of J points on each sub-time period, where J is a preset power of 2, and obtain the frequency domain model x of the signal received by the lth observation station at the jth frequency point in the nth sub-time period l,n (f j 0 is as follows, n is an integer greater than or equal to 1 and less than or equal to N, and j is an integer greater than or equal to 1 and less than or equal to J:

[0058]

[0059] In the formula, is the signal frequency of the jth frequency point;

[0060] in, c is the speed of light;

[0061] s k,n (f j ) is s k The Fourier transform result of the nth time segment of (t); l,n (f j ) is n l The Fourier transform result of the nth time segment of (t);

[0062] Step 3), for each observation station, the signal of the jth frequency point in all time periods is expressed as the following matrix:

[0063] X l,n (f j )=A l (f j )Φ l S n (f j )+N l (f j )

[0064] Among them A l (f j )=[α l (f j ,p1),α l (f j ,p2),…,α l (f j ,p K )]T is an M×K dimensional matrix, which is the steering vector of the uniform linear array of the lth observation station at the jth frequency point; Φ l =diag(ρ l,1 ,…,ρ l,K );S n (f j )=[s 1,n (f j ),…s K,n (f j )] T ; N l (f j )=[n l,1 (f j ), …, n l,N (f j )];

[0065] Step 4), using the relevant information between the observation stations can improve the positioning accuracy of the algorithm to a certain extent. By combining the data of L receiving observation stations, we can get the joint receiving data model as follows:

[0066]

[0067] In the formula,

[0068] Step 5), calculate the sampling covariance matrix of the joint reception data model right Perform eigenvalue decomposition to obtain the matrix U(f j ); Take the matrix U(f j ) to obtain the signal subspace U at the jth frequency point s (f j );

[0069] Step 6), select the signal subspace U corresponding to the starting frequency f0 in the signal subspace under J frequency points s (f0), let the focusing error matrix of the jth frequency point be For C(f j ) performs singular value decomposition to obtain the matrix U composed of left singular vectors c (f j ) and the matrix V consisting of the right singular vectors c (f j );

[0070] Then the focusing matrix T(f j )=U c (f j )V c (f j ) H, the covariance matrix after focusing corresponding to the jth frequency point

[0071] Step 7) Average the focused covariance matrices corresponding to each frequency point to obtain the final covariance matrix

[0072] right Perform eigenvalue decomposition to obtain the matrix U composed of its eigenvectors, and take out the K+1th column to the last column of the matrix U as the final noise subspace matrix U n ;

[0073] Step 8), dividing the area to be located according to a preset grid side length threshold to form a number of grids, and setting the coordinates of the center point of each grid as the coordinates of the grid;

[0074] Step 9), construct the cost function about the grid coordinate p as follows:

[0075]

[0076] In the formula,

[0077] because γ=[ρ1,…,ρ L ] T is the unknown channel attenuation generated by the L observation stations in the process of receiving the radiation source;

[0078] The cost function is converted to:

[0079]

[0080] where e1=[1,0,...,0] T is an L×1 dimensional matrix;

[0081] Step 10), let Solving the converted cost function using the Lagrange multiplier method yields:

[0082]

[0083] Substitute the solved γ into the cost function The radiation source position estimation function is obtained

[0084] Step 11), input the coordinates of each grid into the position estimation function in turn The estimated value corresponding to each grid is obtained, and the grids corresponding to the largest K estimated values ​​are the positions of the K far-field radiation sources to be determined.

[0085] like Figure 3 The figure shows the two-dimensional spatial spectrum obtained by calculating the cost function in the possible area by the method of the present invention, where the radiation sources are located at (23, 52) m and (68, 38) m. It can be seen that the two-dimensional spatial spectrum has two maximum values ​​at the positions corresponding to the two radiation sources, which means that the radiation source can be located by finding the K positions with the largest two-dimensional spatial spectrum.

[0086] Figure 4 It is a scatter plot of the positioning results of the method of the invention, where two radiation sources are considered to be located at p1 = [23m, 52m] respectively. T and p2 = [68m, 38m] T The four observation stations are located at u1=[0m,0m] T 、u2=[0m,100m] T 、u3=[100m,0m] T and u4=[100m,100m] T . Take the number of snapshots J = 500, the number of array elements M = 8, and the signal-to-noise ratio as 10dB. Figure 4 It can be seen that the method of the present invention can locate the radiation source with high precision.

[0087] The performance estimation standard of the present invention is the root mean square error (RMSE), and the RMSE of the radiation source positioning estimation is defined as:

[0088]

[0089] in, It represents the estimated result of the kth radiation source position in the mcth test, MC represents the number of Monte Carlo simulations, and the number of Monte Carlo simulations in this section is 500.

[0090] Figure 5 The performance curve of the root mean square error of the positioning result obtained by the method of the present invention and the signal-to-noise ratio is shown in FIG. Figure 4 Same. From Figure 5 It can be seen that the present invention achieves higher positioning accuracy.

[0091] Figure 6 The performance curve of the root mean square error of the positioning result obtained by the method of the present invention and the signal-to-noise ratio is shown in FIG. Figure 4 Same. From Figure 6 It can be seen that the present invention achieves higher positioning accuracy.

[0092] It will be understood by those skilled in the art that, unless otherwise defined, all terms (including technical and scientific terms) used herein have the same meaning as those generally understood by those skilled in the art in the art to which the present invention belongs. It should also be understood that terms such as those defined in common dictionaries should be understood to have meanings consistent with the meanings in the context of the prior art, and will not be interpreted with idealized or overly formal meanings unless defined as herein.

[0093] The specific implementation methods described above further illustrate the objectives, technical solutions and beneficial effects of the present invention in detail. It should be understood that the above description is only a specific implementation method of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A broadband radiation source direct positioning algorithm with data fusion and frequency focusing, characterized in that it includes the following steps: Step 1) Set up L observation stations with precisely known positions in the area to be located, each equipped with an M-element uniform linear array, with a spacing of d between each element, and their positions are u l =[x l ,y l ], l = 1, 2, ..., L, used to locate K far-field radiation sources; The positions of the K far-field radiation sources are respectively p k =[x k ,y k ], k = 1, 2, ..., K; Let the far-field radiation source emit independent broadband stationary signals, the starting frequency and cutoff frequency of the signal are f0 and f1 respectively; The signal x received by the uniform linear array at the lth observation station at time t is l (t) is: Among them, ρ l,k is the unknown channel attenuation generated by the lth observation station in the process of receiving the kth radiation source, n l (t) is the signal additive noise received by the lth observation station at time t, s k (t) is the signal emitted by the kth far-field radiation source at time t, α l (p k 0 represents the array steering vector when the lth observation station receives the kth far-field radiation source, τ l (p k ) is the time delay of the signal received by the lth observation station compared to the time when the kth signal is received by the first observation station; Step 2), divide the received signal of each observation station in the Ts time period into N sub-time periods, perform a fast Fourier transform of J points on each sub-time period, where J is a preset power of 2, and obtain the frequency domain model x of the signal received by the lth observation station at the jth frequency point in the nth sub-time period l,n (f j 0 is as follows, n is an integer greater than or equal to 1 and less than or equal to N, and j is an integer greater than or equal to 1 and less than or equal to J: In the formula, is the signal frequency of the jth frequency point; in, c is the speed of light; s k,n (f j ) is s k The Fourier transform result of the nth time segment of (t); l,n (f j ) is n l The Fourier transform result of the nth time segment of (t); Step 3), for each observation station, the signal of the jth frequency point in all time periods is expressed as the following matrix: X l,n (f j )=A l (f j )Φ l S n (f j )+N l (f j ) Among them A l (f j )=[α l (f j ,p1),α l (f j ,p2),…,α l (f j ,p K )] T is an M×K dimensional matrix, which is the steering vector of the uniform linear array of the lth observation station at the jth frequency point; Φ l =diag(ρ l,1 ,…,ρl l,K );S n (f j )=[s 1,n (f j ),…s K,n (f j )] T ; N l (f j )=[n l,1 (f j ), …, n l,N (f j )]; Step 4) Combine the signals received by all observation stations at the jth frequency point in N time periods to obtain the joint reception data model X n (f j )as follows: In the formula, Step 5), calculate the sampling covariance matrix of the joint reception data model right Perform eigenvalue decomposition to obtain the matrix U(f j ); Take the matrix U(f j ) to obtain the signal subspace U at the jth frequency point s (f j ); Step 6), select the signal subspace U corresponding to the starting frequency f0 in the signal subspace under J frequency points s (f0), let the focusing error matrix of the jth frequency point be For C(f j ) performs singular value decomposition to obtain the matrix U composed of left singular vectors c (f j ) and the matrix V consisting of the right singular vectors c (f j ); Then the focusing matrix T(f j )=U c (f j )V c (f j ) H , the covariance matrix after focusing corresponding to the jth frequency point Step 7) Average the focused covariance matrices corresponding to each frequency point to obtain the final covariance matrix right Perform eigenvalue decomposition to obtain the matrix U composed of its eigenvectors, and take out the K+1th column to the last column of the matrix U as the final noise subspace matrix U n ; Step 8), dividing the area to be located according to a preset grid side length threshold to form a number of grids, and setting the coordinates of the center point of each grid as the coordinates of the grid; Step 9), construct the cost function about the grid coordinate p as follows: In the formula, because γ=[ρ1,…,ρ L ] T is the unknown channel attenuation generated by the L observation stations in the process of receiving the radiation source; The cost function is converted to: where e1=[1,0,...,0] T is an L×1 dimensional matrix; Step 10), let Solving the converted cost function using the Lagrange multiplier method yields: Substitute the solved γ into the cost function The radiation source position estimation function is obtained Step 11), input the coordinates of each grid into the position estimation function in turn The estimated value corresponding to each grid is obtained, and the grids corresponding to the largest K estimated values ​​are the positions of the K far-field radiation sources to be determined.

Images