Subspace-based direction of arrival estimation method based on distributed algorithm

Abstract

This invention belongs to the field of signal processing, and particularly relates to direction-of-arrival (DOA) estimation for electromagnetic and sonar signals. Specifically, it is a subspace-based DOA estimation method based on a distributed algorithm, which can be used for passive localization and target detection. The method includes the following steps: S1, establishing a distributed array model at the receiver; S2, calculating the signal subspace of the received signal; S3, calculating the weight w(θ); S4, calculating the spatial spectrum and performing a spatial domain search to find the DOA value. The proposed method achieves distributed solution of the algorithm, thus fundamentally avoiding the shortcomings of centralized algorithms. Simultaneously, it maintains good estimation performance, similar to that of centralized algorithms, and outperforms distributed MUSIC methods at low signal-to-noise ratios.

Description

Technical Field

[0001] This invention belongs to the field of signal processing, and particularly relates to the direction of arrival estimation of electromagnetic signals and sonar signals. Specifically, it is a subspace-based direction of arrival estimation method based on a distributed algorithm, which can be used for passive localization and target detection. Background Technology

[0002] DOA (Direction of Arrival) estimation is an important research direction in array signal processing. It refers to using sensor arrays at different geographical locations to receive signals, obtain discrete observation data, and then process the data through a series of mathematical algorithms to obtain useful information from the source signal while suppressing interference. It has wide applications in military and civilian fields such as autonomous driving, wireless communication, and radar.

[0003] Traditional DOA estimation techniques are centralized algorithms, requiring all array elements to transmit received signal information to a central processor for processing. However, this approach is limited in many scenarios. For example, in a large array composed of geographically dispersed subarrays, a large amount of observation data needs to be transmitted between subarrays to a central processor; however, the capacity of the communication links limits data transmission between subarrays. Furthermore, centralized processing fails completely in many scenarios, such as in sonar environments, where signals received by widely spaced subarrays lose coherence, making centralized coherent processing impossible. Therefore, to improve the adaptability of the algorithm, a decentralized DOA estimation algorithm needs to be designed. Summary of the Invention

[0004] The purpose of this invention is to address the shortcomings of the prior art by proposing a subspace-based direction-of-arrival estimation method based on a distributed algorithm. This method decentralizes the solution process of the subspace-based method, thereby avoiding the influence of traditional centralized algorithms, and achieving estimation performance similar to that of centralized algorithms.

[0005] The specific technical solution adopted in this invention is as follows:

[0006] A subspace-based direction-of-arrival estimation method based on a distributed algorithm, comprising the following steps:

[0007] S1. Establish a distributed array model at the receiving end;

[0008] S2. Calculate the signal subspace of the received signal.

[0009] S3. Calculate the weight w(θ);

[0010] S4. Calculate the spatial spectrum and perform a spatial domain search to find the DOA value.

[0011] A further improvement of the present invention is that, in step S1, a uniform linear array, denoted as array x, is first constructed using M array elements with a spacing equal to half the wavelength of the incident narrowband signal; it is assumed that there are K far-field narrowband signals θ. k For k = 1, ..., K, after receiving L snapshots, the signal received by array x is represented as:

[0012] Y = A(θ)X + N

[0013] Where Y is the array element received data matrix, X is the incident signal data matrix, and N is the noise data matrix, with noise and signal being independent of each other; A(θ)=[a(θ1),…,a(θ) K ] is an array manifold matrix, and is the array steering vector a(θ) k The set of ), a(θ) k )=[1,e iπcosθ ,…,e iπ ( M-1 ) cosθ ] T .

[0014] These M arrays are divided into P subarrays (nodes) with equal array elements, meaning the number of array elements in each subarray is C = M / P. Each subarray has a subprocessor that can only store and process data accessible to that subarray. Furthermore, the neighboring nodes of each node are defined, with the neighboring node of the i-th node represented as... Therefore, the original signal is also divided into components belonging to each subarray. The definition of the original model is added as follows: (■) j This represents the vector (or matrix) component belonging to the j-th subarray, where j = 1, ..., P. For example, the incident signal X can be represented as X = [x1, x2, ..., xp]. P ], where x1 is the incident signal received by the first subarray.

[0015] In a further improvement of the present invention, in step S2, the signal subspace of the received data is calculated. Specifically:

[0016] S21. First, calculate the first eigenvector of the covariance matrix R, randomly initialize the vector u(0), multiply it by the covariance matrix, and iterate I. pw Next, that is

[0017]

[0018] in,(·) H For the matrix (and the conjugate of the vector), x j (t) is the incident signal vector component belonging to the j-th subarray in the t-th snapshot, u 1,j (Ipw ) for the first pw The eigenvector components belonging to the j-th subarray in the next iteration. AC(·) is the average consensus protocol operation, which means averaging the data within the brackets of all subarrays, i.e.

[0019]

[0020] The average consensus protocol is an algorithm based on the Laplace matrix of a distributed network. The Laplace matrix is ​​defined as follows: when node i can communicate with node j, l ij = -1; otherwise, it is 0; the node's degrees of freedom (how many nodes it can communicate with) is l. ii (l jj The value of ), that is

[0021]

[0022] Suppose that the Laplace matrix B has R distinct eigenvalues ​​λ1,…,λ2. R , and λ R =0, then the initial weight is

[0023]

[0024] The subsequent weights are:

[0025] W(t) = B - λ t+1 I

[0026] If we want to find the mean of x(t), i.e., AC(x j If (t)), then the iterative formula for this protocol is:

[0027] x(t+1)=W(t)x(t)

[0028] S22. Calculate the magnitude of the first eigenvector:

[0029]

[0030] S23. Each node normalizes the data it receives:

[0031]

[0032] S24, as in S21, through a sufficient number of iterations I pw Calculate the remaining K-1 eigenvectors belonging to the large eigenvalues.

[0033]

[0034] S25. After obtaining the result, normalize it as in S23 and S24.

[0035] S26. Find the signal subspace. E x =[u1,…,u K The large eigenvalues ​​of the covariance matrix constitute the structure.

[0036]

[0037] Among them, u i,j Let i be the component of the i-th eigenvector that belongs to the j-th submatrix.

[0038] A further improvement to the present invention is that in step S3, the weight w(θ) is solved, specifically in the form of:

[0039]

[0040] First, the distributed conjugate gradient method is used to calculate... This is an iterative solution method. The specific steps for each iteration are as follows:

[0041] S31. Calculate the search step size α i :

[0042]

[0043] Where, r i Let d be the residual of the i-th iteration. i This represents the search direction for the i-th iteration.

[0044]

[0045] S32, Each node calculates the update target h. i :

[0046] h i+1,j =h i,j +α i d i,j

[0047] S33. Update the residual at each node:

[0048] r i+1,j =r i,j -α i R yy d i =r i,j -α i o

[0049] S34. Update the Schmitt component β used for the search direction. i :

[0050]

[0051] in, Distributed solutions in the form of S31 Seeking the method.

[0052] S35. Each node updates the search direction:

[0053] d i+1,j =r i+1,j +β i d i,j

[0054] When the number of iterations exceeds the set maximum number of iterations or the target result updated between two consecutive iterations changes very little, Then stop iterating;

[0055] S36. After obtaining the repeating unit h, w ES The formula for calculating (θ) is:

[0056]

[0057] Furthermore, S4 calculates the spatial spectrum P(θ):

[0058]

[0059] Let AC(x) j,i w j )=f i ,but

[0060]

[0061] Because f i Each processor receives the data, so each subarray processor can process L data points belonging to its own node (summing or averaging). but

[0062]

[0063] The DOA estimate can be obtained by scanning the spatial domain.

[0064] The beneficial effects of this invention are as follows: The method proposed in this invention realizes the distributed solution of the algorithm, thereby avoiding the shortcomings of the centralized algorithm in essence; at the same time, it can maintain good estimation performance, similar to the estimation performance of the centralized algorithm, and outperforms the distributed MUSIC method at low signal-to-noise ratio. Attached Figure Description

[0065] Figure 1 This is the overall flowchart of the present invention.

[0066] Figure 2 This is a schematic diagram of the distributed structure of the present invention.

[0067] Figure 3This is a performance graph of an example of the present invention under low signal-to-noise ratio.

[0068] Figure 4 This is a performance graph of an example of the present invention when photographed on a small block.

[0069] Figure 5 This is a performance chart illustrating the examples of the present invention and the centralized method.

[0070] Figure 6 This is a performance graph of the example of the present invention and the distributed method. Detailed Implementation Plan

[0071] To enhance understanding of the present invention, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. These embodiments are only used to explain the present invention and do not constitute a limitation on the scope of protection of the present invention.

[0072] like Figure 1 and Figure 2 As shown, a subspace-based direction-of-arrival estimation method based on a distributed algorithm includes the following steps:

[0073] S1. Establish a distributed array model at the receiving end;

[0074] First, a uniform linear array, denoted as array x, is constructed using M array elements with a spacing equal to half the wavelength of the incident narrowband signal. Then, assume there are K far-field narrowband signals θ. k For k = 1, ..., K, after receiving L snapshots, the signal received by array x is represented as:

[0075] Y = A(θ)X + N

[0076] Where Y is the array element received data matrix, X is the incident signal data matrix, and N is the noise data matrix, with noise and signal being independent of each other. A(θ)=[a(θ1),…,a(θ) K ] is an array manifold matrix, and is the array steering vector a(θ) k The set of ), a(θ) k )=[1,e iπcosθ ,…,e iπ(M-1)cosθ ] T .

[0077] These M arrays are divided into P subarrays with equal array elements, meaning the number of array elements in each subarray is C = M / P. Each subarray has a subprocessor that can only store and process data accessible to that subarray. Furthermore, the neighboring nodes of each node are defined, with the neighboring nodes of the i-th node represented as... Therefore, the original signal is also divided into components belonging to each subarray. The definition of the original model is added as follows: (■) jThis represents a vector (or matrix) component of dimension C that belongs to the j-th submatrix, where j = 1, ..., P. For example, the incident signal X can be represented as X = [x1, x2, ..., x...]. P ], where x1 is the C×1 dimensional incident signal received by the first subarray.

[0078] S2. Calculate the signal subspace of the received signal using the distributed power method.

[0079] S21. First, calculate the first eigenvector of the covariance matrix R, randomly initialize the vector u(0), multiply it by the covariance matrix, and iterate I. pw Next, that is

[0080]

[0081] in,(·) H For the matrix (and the conjugate of the vector), x j (t) is the incident signal vector component belonging to the j-th subarray in the t-th snapshot, u 1,j (I pw ) for the first pw The eigenvector components belonging to the j-th subarray in the next iteration. AC(·) is the average consensus protocol operation, which means averaging the data within the brackets of all subarrays, i.e.

[0082]

[0083] The average consensus protocol is an algorithm based on the Laplace matrix of a distributed network. The Laplace matrix is ​​defined as follows: when node i can communicate with node j, l ij = -1; otherwise, it is 0; the node's degrees of freedom (how many nodes it can communicate with) is l. ii (l jj The value of ), that is

[0084]

[0085] Suppose that the Laplace matrix B has R distinct eigenvalues ​​λ1,…,λ2. R , and λ R =0, then the initial weight is

[0086]

[0087] The subsequent weights are,

[0088] W(t) = B - λ t+1 I

[0089] The iterative formula for this protocol is:

[0090] x(t+1)=W(t)x(t)

[0091] S22. Calculate the magnitude of the first eigenvector:

[0092]

[0093] S23. Each node normalizes the data it receives:

[0094]

[0095] S24, as in S21, through a sufficient number of iterations I pw Calculate the remaining K-1 eigenvectors belonging to the large eigenvalues.

[0096]

[0097] S25. After obtaining the result, normalize it as in S23 and S24.

[0098] S26. Find the signal subspace. E x =[u1,…,u K The large eigenvalues ​​of the covariance matrix constitute the structure.

[0099]

[0100] Among them, u i,j Let i be the component of the i-th eigenvector that belongs to the j-th submatrix.

[0101] S3. Calculate the weight w(θ) using the distributed conjugate gradient method, specifically in the form of:

[0102]

[0103] First, the distributed conjugate gradient method is used to calculate... This is an iterative solution method. The specific steps for each iteration are as follows:

[0104] S31. Calculate the search step size α i :

[0105]

[0106] Where, r i Let d be the residual of the i-th iteration. i This represents the search direction for the i-th iteration.

[0107]

[0108] S32, Each node calculates the update target h. i :

[0109] h i+1,j =h i,j +α i d i,j

[0110] S33. Update the residual at each node:

[0111] r i+1,j =r i,j -α i R yy d i =r i,j -α i o

[0112] S34. Update the Schmitt component β used for the search direction. i :

[0113]

[0114] in, Distributed solutions in the form of S31 Seeking the method.

[0115] S35. Each node updates the search direction:

[0116] d i+1,j =r i+1,j +β i d i,j

[0117] When the number of iterations exceeds the set maximum number of iterations or the target result updated between two consecutive iterations changes very little, Then stop iterating.

[0118] S36. After obtaining the repeating unit h, w ES The formula for calculating (θ) is:

[0119]

[0120] S4. Calculate the spatial spectrum and perform a spatial domain search to find the DOA value. Calculate the spatial spectrum P(θ):

[0121]

[0122] Let AC(x) j,i w j )=f i ,but

[0123]

[0124] Because f i Each processor receives the data, so each subarray processor can process L data points belonging to its own node (summing or averaging). but

[0125]

[0126] The DOA estimate can be obtained by scanning the spatial domain.

[0127] The effects of the present invention will be further described below with reference to simulation examples. Figure 3 , Figure 4 , Figure 6 In this context, MUSIC (Multiple Signal Classification) is an angle estimation method; the Capon method is an angle estimation method; and the distributed ES method represents an example method.

[0128] The simulation examples below all establish a matrix of 12 ULAs with an element spacing of λ / 2, and divide it into P = 6 subarrays, defining the neighboring nodes of each subarray. In the subarray division, the neighboring nodes of each subarray are defined as follows:

[0129]

[0130] Simulation Example 1: Considering a scenario with low signal-to-noise ratio and low snapshot count. Assume that two far-field uncorrelated sources are incident on the distributed array from two directions [-6°, 6°]. (1) Low signal-to-noise ratio: the number of snapshots is L = 100, and the signal-to-noise ratio SNR = -5dB; (2) the number of snapshots is L = 10, and the signal-to-noise ratio SNR = 10dB. Figure 3 Spatial spectral diagrams of the distributed ES algorithm, Capon algorithm, and MUSIC algorithm at low signal-to-noise ratio are given. Figure 4 Spatial spectrograms of the three methods are presented at low snapshot numbers. The dashed lines represent the correct source direction of arrival, i.e., the DOA estimation result. It can be seen that all three distributed methods can correctly estimate the DOA result. The side lobes of the distributed ES algorithm are significantly lower than the other two distributed algorithms, and the main lobe width remains relatively thin.

[0131] Simulation Example 2: Verifying whether the distributed Elasticsearch method can approximate the estimation performance of the centralized Elasticsearch method. RMSE (Root Mean Square Error) is used as a metric, and its calculation formula is...

[0132]

[0133] Where mt is the total number of Monte Carlo experiments, and K is the number of information sources. Let θ represent the DOA estimate of the k-th source in the t-th experiment. kLet be the k-th true source value. Assume two uncorrelated far-field signals are incident on the array from [-10° 10°]. The centralized algorithm is built on the complete ULA without partitioning. The signal-to-noise ratio ranges from -20 to 10 dB, with a step size of 4 and a snapshot number of L = 100. The total number of Monte Carlo experiments is mt = 500. Figure 5 The changes in RMSE with SNR for the two methods are presented. It can be seen that, except at -15dB, the distributed ES method is slightly worse than the centralized method, but at other times the distributed method can achieve the same DOA estimation performance as the centralized method.

[0134] Simulation Example 3: Comparing the estimation performance of the distributed Elasticsearch method with other distributed algorithms. The experimental assumptions are consistent with those of Simulation Example 2. Figure 6 The estimation performance of the three methods as a function of signal-to-noise ratio is presented. The comparison shows that at low signal-to-noise ratios, the estimation performance of distributed MUSIC is inferior to the other two methods, while the estimation performance of distributed ES algorithm is good.

[0135] In summary, the method proposed in this invention achieves distributed solution of the algorithm, thereby fundamentally avoiding the shortcomings of centralized algorithms; at the same time, it can maintain good estimation performance, similar to that of centralized algorithms, and outperforms distributed MUSIC methods at low signal-to-noise ratios.

[0136] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely illustrative of the principles of the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the present invention as claimed. The scope of protection of the present invention is defined by the appended claims and their equivalents.

Claims

1. A subspace-based direction-of-arrival estimation method based on a distributed algorithm, characterized in that, Includes the following steps: S1. Establish a distributed array model at the receiving end; S2. Calculate the signal subspace of the received signal. ; S3, Calculate weights ; S4. Calculate the spatial spectrum and perform a spatial domain search to find the DOA value; In step S1, a uniform linear array, denoted as array x, is first constructed using M array elements with a spacing equal to half the wavelength of the incident narrowband signal. K far-field narrowband signals are then set. After receiving L snapshots, the signal received by array x is represented as: ; Where Y is the array element received data matrix, X is the incident signal data matrix, and N is the noise data matrix, with noise and signal being independent of each other. It is an array manifold matrix, and is the array steering vector. The set, ; Divide these M arrays into P subarrays with equal array elements, that is, the number of array elements in each subarray is 1. Each subarray has a subprocessor, which can only store and process data accessible to that subarray. The neighboring nodes of each node are defined, and the neighboring nodes of the i-th node are represented as... The following additions are made to the definition of the variables: Let represent the vector component belonging to the j-th submatrix. The incident signal X is represented as ,in, The incident signal received by the first subarray; In step S2, the signal subspace of the received data is calculated. Specifically: S21. First, calculate the first eigenvector of the covariance matrix R, and then randomly initialize the vector. Multiply it by the covariance matrix and iterate. Next, that is , in, To take the conjugate of a matrix and a vector, It is the incident signal vector component belonging to the j-th subarray in the t-th snapshot. In the first The eigenvector components belonging to the j-th submatrix in the next iteration This is an average consensus protocol operation, which refers to averaging the data within the brackets of all subarrays, i.e. ； The Laplace matrix is ​​defined as follows: when node i can communicate with node j, ... Otherwise, it is 0; the degree of freedom of the node itself is... The value of, i.e. ； Suppose that the Laplace matrix B has R distinct eigenvalues. ,and The initial weights are ； The subsequent weights are ； Require The mean is The iterative formula for this protocol is: ; S22. Calculate the magnitude of the first eigenvector: ; S23. Each node normalizes the data it receives: ; S24 and S21 were passed through a sufficient number of iterations. Calculate the remaining Eigenvectors belonging to large eigenvalues: S25. After obtaining S23 and S24, normalize them. S26. Find the signal subspace. It is composed of the large eigenvalues ​​of the covariance matrix, i.e. ； in, Let i be the component of the i-th eigenvector that belongs to the j-th submatrix.

2. The subspace-based direction-of-arrival estimation method based on a distributed algorithm according to claim 1, characterized in that, In step S3, the weights are solved. The specific form is 。 3. The subspace-based direction-of-arrival estimation method based on a distributed algorithm according to claim 2, characterized in that, In step S3, the distributed conjugate gradient method is used to calculate... The specific steps for each iteration are as follows: S31. Calculate the search step size : ; in, Let be the residual of the i-th iteration. Let be the search direction for the i-th iteration. ； S32, Each node calculates the update target. : ； S33. Update the residual at each node: ； S34. Update the Schmidt component used for the search direction. : ； in, Distributed solutions in the form of S31 Seeking the Dharma S35. Each node updates the search direction: ; When the number of iterations exceeds the set maximum number of iterations or the target result updated between two consecutive iterations changes very little, Then stop iterating. S36, Obtain repeating units back, The calculation formula is: .

4. The subspace-based direction-of-arrival estimation method based on a distributed algorithm according to claim 3, characterized in that, Step S4 calculates the spatial spectrum. : make ,but ； because Each processor acquires the data, so each subarray processor can process data belonging to its own node. Data, record ,but ； The DOA estimate can be obtained by scanning the spatial domain.

Images