A Meshless Coherent Signal Direction-of-Arrival Estimation Method Based on Alternating Projection

Abstract

A gridless coherent signal direction-of-arrival estimation method based on alternating projection according to the present invention is characterized in that: the patented method uses a hydrophone array installed on an underwater vehicle to receive signal data to construct a positive semi-definite matrix S0, and uses the projection method to alternately project S0 onto a rank-constrained set and a Toeplitz matrix with averaged main diagonal elements. Secondly, the convergence degree of S0 is judged, and then T i (u) is subjected to reduced-dimensional space smoothing-singular value decomposition. Finally, z is obtained by using the polynomial root-finding method, and the incident angle of the coherent signal is solved by using the corresponding mapping relationship. The patented method effectively realizes the angle estimation of coherent signals. Secondly, the polynomial root-finding method is used instead of spectral peak search, which effectively reduces the computational amount, improves the robustness and real-time performance of the algorithm, and has high estimation accuracy, having important engineering application value.

Description

Technical Field

[0001] The present invention relates to a method for estimating the direction of arrival of a hydrophone array. Specifically, it is a method for estimating the direction of arrival of coherent signals without a grid based on alternating projection. Background Art

[0002] The estimation of the direction of arrival (DOA) has always been a research hotspot in array signal processing and has wide applications in fields such as radar and sonar. Traditional DOA estimation algorithms mainly include the Multiple Signal Classification (MUSIC) algorithm and the Estimation of Signal Parameters via Rotational Invariance Techniques (ESPRIT) algorithm. The above algorithms are mainly based on eigenvalue decomposition and have a large computational complexity. In addition, there are subspace fitting algorithms, and the representative one is the Maximum Likelihood algorithm. The estimation performance of this type of algorithm is improved compared with the eigenvalue decomposition algorithms, but its solution process is realized through multi-dimensional search and the computational amount is huge, so its practicability is low.

[0003] In recent years, the theory of compressive sensing has provided a new idea for signal reconstruction. By using the sparsity of signals in the spatial domain, DOA estimation can be effectively carried out. However, traditional compressive sensing DOA estimation algorithms still have defects. That is, the estimation accuracy of its algorithm is related to the sparse division of the signal space. If the sparse grid division density is low, for actual signals, they may not exactly fall on the grid points but between two grid points, resulting in grid mismatch. If the sparse grid division density is increased, the equal-distance constraint condition between basis atoms will not be satisfied, leading to a decline in DOA estimation performance.

[0004] In response to this, Tang proposed the atomic norm minimization method to accurately recover signals and redefined the problem as a semidefinite programming problem. In the case of uncorrelated signals, this method was proven to be a multi-snapshot implementation of maximum likelihood estimation, but no specific implementation process was given for coherent signals. Summary of the Invention

[0005] Therefore, in view of the above problems, this patent improves the algorithm and proposes a method for estimating the direction of arrival of coherent signals without a grid based on alternating projection. Its characteristics are as follows: The method of this patent uses the hydrophone array of an underwater vehicle to receive complex numbers to construct a semidefinite matrix S0, and uses the projection method to alternately project S0 onto the rank constraint set P R and the Toeplitz matrix P with the average of the main diagonal elementsT First, judge the convergence degree of S0, and then perform T i (u) for dimensionality reduction space smoothing - singular value decomposition. Finally, use the method of polynomial root finding to find z, and use the corresponding mapping relationship to solve the incident angle of the coherent signal. The method of this patent effectively realizes the angle estimation of coherent signals. Secondly, using the method of polynomial root finding instead of spectral peak search effectively reduces the computational complexity and improves the real-time performance of algorithm utilization, which has important application value.

[0006] A gridless coherent signal DOA estimation method based on alternating projection in this patent has high robustness and azimuth estimation accuracy. The signal processing of the method of this patent includes the following steps:

[0007] Step 1: The number of array elements of the hydrophone is M, the array element spacing is d, and the received coherent signal data is L is the number of sampling data points of the hydrophone array elements. According to the AIC decision criterion, judge the number of targets K, set the convergence limit ε of the matrix, and form a positive semi-definite matrix with the received coherent signal data Y and the data covariance matrix T0(u) where T0(u) = YY H , T0(u) is a conjugate symmetric Toeplitz matrix generated by the vector u, Z is a free conjugate symmetric matrix, (·) H represents the conjugate transpose operator. Initialize Z0 = zeros(L, L), and zeros(L, L) is an all-zero matrix of L*L;

[0008] Step 2: Project the positive semi-definite matrix onto the rank constraint set to obtain The specific process of its projection processing is as follows: First, perform eigenvalue decomposition on the positive semi-definite matrix eig(·) is the eigenvalue decomposition operator. Secondly, make the eigenvalues real Dz = real(Dz), where real(·) is the real part taking operator. Finally, reconstruct the projection matrix as S according to the eigenvalues and their corresponding eigenvectors 0-i =(Vz*Dz*Vz H +(Vz*Dz*Vz H ) H ) / 2;

[0009] Step 3: Project T0(u) onto the Toeplitz matrix with the main diagonal elements averaged to obtain T i (u). The specific process of its projection processing is as follows: First, according to the formula u(i) = mean([diag(T0(u),(i - 1));(diag(T0(u),-(i - 1))) *), for \(i = 1, 2, \ldots, M\), the updated vector \(u=\left[u(1), u(2), \ldots, u(M)\right]\) is obtained, where \(\text{diag}(\cdot, \text{num})\) represents taking the elements on the diagonal to form a vector. When \(\text{num} = 0\), it means taking the elements on the main diagonal of the matrix to form a vector. When \(\text{num}>0\), it means taking the diagonal parallel to the main diagonal to form a vector. When \(\text{num}<0\), it means taking the diagonal parallel to the sub - diagonal to form a vector. \(\text{mean}(\cdot)\) represents taking the mean of each column vector of the matrix to form a row vector, \((\cdot)^*\) represents the conjugate symbol, and \([\cdot;\cdot]\) represents concatenating matrices by columns. Then, the Toeplitz matrix is generated using the vector \(u(i)\).

[0010] Step Four: Substitute \(T i (u)\) into the corresponding position of \(S 0-i to obtain \(S i ;\)

[0011] Step Five: Calculate \(\lambda=\left\|S i - S_0\right\|\) F where \(\|\cdot\|\) F is the Frobenius norm. If \(\lambda\leq\varepsilon\), then take out \(T i (u)\) in the positive semi - definite matrix \(S i \) and enter Step Six. Otherwise, set \(S_0 = S i \) and return to Step Two;

[0012] Step Six: Perform reduced - dimensional space smoothing - singular value decomposition on \(T i (u)\). The process of performing reduced - dimensional space smoothing on \(T i (u)\) to obtain \(T i-ss (u)\) is as follows: First, calculate the smoothing times \(N = M-(M sub - 1)\), where \(M sub \) is the number of array elements in the sub - array. The covariance matrix corresponding to each sub - array is \(\text{sub}(T i (u, i)), i = 1, 2, \ldots, N\), and \(T i-ss (u)=\sum(\text{sub}(T i (u, i))) / N\), where \(\text{sub}(T i (u, i))=T i (u)(i:i + M sub - 1, i:i + M sub - 1)\) and \(\sum(\cdot)\) is the summation symbol;

[0013] Step Seven: Perform singular value decomposition on \(T i-ss (u)\) as \([U,\sum,V]=\text{svd}(T i-ss (u))\), where \(U N= U(:, K + 1:M - 1), where U is the matrix corresponding to the left singular vectors. where

[0014] Step 8: Convert the solved polynomial roots z to the estimated values of the signal's direction of arrival angles. The conversion formula is θ = -arcsin(angle(2×z)) * 180 / π, where angle(·) represents taking the phase angle corresponding to the complex number in [-π / 2, π / 2]. Description of the Drawings

[0015] Figure 1 is the relationship curve between the root mean square error and the signal-to-noise ratio of the direction of arrival estimation method of this patent;

[0016] Figure 2 is the relationship curve between the root mean square error and the number of data sampling points of the direction of arrival estimation method of this patent. Detailed Embodiments

[0017] Now, the present invention will be further described in conjunction with embodiments and drawings:

[0018] The 1st embodiment: The simulation conditions are as follows: The receiving hydrophone array is a uniform linear array, the number of array elements M = 20, the array element spacing d = 0.06 m is half the wavelength of the signal, the incident angles of the coherent signals are randomly generated in (-90°, 90°], and the coherent signals are incident on the uniform linear array from different angles. The number of data points received by the hydrophone array elements is L = 200. With the above simulation conditions, the specific implementation process is as follows:

[0019] Step 1: Combine the data received by the array elements and determine the number of targets K = 3 according to the AIC decision criterion. Set the convergence limit of the matrix ε = 1×10 -2 , and for the data of the received coherent signals where is the number of single data points received by the array elements, and the data covariance matrix forms a positive semi - definite matrix where T0(u) = YY H , the first column is the conjugate - symmetric Toeplitz matrix of u, Z is a free conjugate - symmetric matrix, and initialize Z0 = zeros(200, 200);

[0020] Step 2: Project the positive semi - definite matrix onto the rank - constrained set to obtain First, perform eigenvalue decomposition on the positive semi - definite matrix Secondly, make the eigenvalues real - valued Dz = real(Dz), where real(·) is the real - part operator. Finally, reconstruct the matrix S according to the eigenvalues and their corresponding eigenvectors 0-i =(Vz * Dz * VzH +(Vz * Dz * Vz H ) H ) / 2;

[0021] Step 3: Project T0(u) onto the Toeplitz matrix with averaged elements along the main diagonal to obtain T i (u). First, according to the formula u(1) = mean([diag(T0(u), 0); (diag(T0(u), 0)) * ) u(2) = mean([diag(T0(u), 1); (diag(T0(u), -1)) * ) and update u(3), u(4) successively until u(20) = mean([diag(T0(u), 19); (diag(T0(u), -19)) * ) to obtain the vector update value u = [u(1), u(2),..., u(20)]. Then, use the vector u to generate the Toeplitz matrix

[0022] Step 4: Substitute T i (u) into the corresponding position of S 0-i to obtain

[0023] Step 5: Calculate λ = ||S i - S0|| F . If λ ≤ 1 × 10 -2 , then extract T i from the positive semi - definite matrix S i (u) and enter Step 6. Otherwise, set S0 = S i and return to Step 2;

[0024] Step 6: Perform dimensionality - reduced space smoothing - singular value decomposition on T i (u). Process T i (u) by dimensionality - reduced space smoothing to obtain T i-ss (u). The specific process of its dimensionality - reduced space smoothing is as follows: First, calculate the smoothing times 2 = N = M - (M sub - 1) = 20 - (19 - 1), where M sub = 19 is the number of sub - array elements. The covariance matrices corresponding to each sub - array are sub(T i (u, 1)), sub(T i (u, 2)) T i - ss (u) = (sub(T i (u, 1)) + sub(T i (u, 2))) / 2

[0025] Step Seven: For T i-ss (u), perform singular value decomposition [U, ∑, V]=svd(T i-ss (u)), where U N =U(:, K + 1:M - 1)=U(:, 4:19). Here, U is the matrix corresponding to the left singular vectors. where Use the built-in fminsearch function in the matlab software to solve for the value of z;

[0026] Step Eight: Convert the solved value of z into the estimated value of the direction of arrival angle of the signal. The conversion formula is θ=-arcsin(angle(2×z))*180 / π. Here, angle(·) represents taking the phase angle corresponding to the complex number in [-π / 2, π / 2].

[0027] The second embodiment: Calculate the curve of the relationship between the estimated performance of the direction of arrival estimation method of this patent and the signal-to-noise ratio, and obtain the effect diagram as Figure 1 shown. The computer simulation conditions in this embodiment are as follows:

[0028] Adopt a uniform linear array with 20 array elements. The element spacing d = 0.06m is half a wavelength. The signal incident angle is randomly generated in (-90°, 90°]. The number of snapshots L = 200. Start changing the signal-to-noise ratio from -1dB and increase it in steps of 2dB to 15dB. Conduct 300 Monte Carlo experiments, simulate through the matlab software system, and observe the simulation results. The simulation results are as Figure 1 shown.

[0029] From Figure 1 it can be seen that as the signal-to-noise ratio increases, the estimation error of the method proposed in this patent is lower than that of the SS-MUSIC method. When the signal-to-noise ratio is greater than 6dB, the estimation error of the SS-MUSIC method basically remains unchanged, while for the method proposed in this patent, as the signal-to-noise ratio increases, the estimation error continues to decrease, and the estimation accuracy is higher.

[0030] The third embodiment: Calculate the relationship between the estimated performance of the direction of arrival estimation method of this patent and the number of sampled data points, and obtain the effect diagram as Figure 2 shown. The computer simulation conditions in this embodiment are as follows:

[0031] The receiving array is a uniform linear array with 20 array elements. The element spacing d = 0.06m is half a wavelength. The signal incident angle is randomly generated in (-90°, 90°]. The signal-to-noise ratio is 5dB. Change the number of sampled data points L from 30 and increase it in steps of 30 to 300. Conduct 200 Monte Carlo experiments, simulate through the matlab software system, and observe the simulation results. The simulation results are as Figure 2 shown.

[0032] From Figure 2 It can be seen that as the number of data sampling points increases, the estimation performance of the method proposed in this patent is significantly better than that of the SS-MUSIC method. As the number of sampled data points increases, the estimation errors of the two methods basically remain unchanged, and the estimation performance of the method in this patent is always better than that of SS-MUSIC.

[0033] The specific examples described in this article are only illustrative of the present invention. Those skilled in the art to which the present invention pertains can make various modifications, supplements, or use similar means of substitution to the described specific examples, but will not deviate from the present invention or exceed the scope defined by the appended claims.

Claims

1. A meshless coherent signal direction-of-arrival estimation method based on alternating projection, characterized in that: The direction-of-arrival estimation method includes the following steps: Step 1: The number of array elements of the hydrophone is M, the array element spacing is d, and the received coherent signal data is L is the number of sampling data points of the hydrophone array elements. According to the AIC decision criterion, the number of targets K is discriminated, and the convergence limit ε of the matrix is set. The received coherent signal data Y and the data covariance matrix T0(u) are combined to form a positive semi - definite matrix where T0(u)=YY H , T0(u) is a conjugate - symmetric Toeplitz matrix generated using the vector u, Z is a free conjugate - symmetric matrix, (·) H represents the conjugate transpose operator. Initialize Z0 = zeros(L,L), and zeros(L,L) is an L*L all - zero matrix; Step 2: Project the positive semi - definite matrix onto the rank - constrained set to obtain The specific process of the projection is as follows: First, perform eigenvalue decomposition on the positive semi - definite matrix eig(·) is the eigenvalue decomposition operator. Second, make the eigenvalues real Dz = real(Dz), where real(·) is the real - part extraction operator. Finally, reconstruct the projection matrix S according to the eigenvalues and their corresponding eigenvectors 0-i =(Vz * Dz * Vz H +(Vz * Dz * Vz H ) H ) / 2; Step 3: Project \(T_0(u)\) onto the Toeplitz matrix averaged along the main diagonal elements to obtain \(T(u)\), and the specific process of the projection is as follows: First, according to the formula \(u(i)=\text{mean}([\text{diag}(T_0(u),(i - 1));(\text{diag}(T_0(u),-(i - 1))])\), \(i = 1,2,\cdots,M\) to obtain the vector update value \(u=[u(1),u(2),\cdots,u(M)]\), where \(\text{diag}(\cdot,\text{num})\) represents taking the elements of the diagonal to form a vector. When \(\text{num}=0\), it means taking the elements of the main diagonal of the matrix to form a vector. When \(\text{num}>0\), it means taking the diagonal parallel to the diagonal to form a vector. When \(\text{num}<0\), it means taking the diagonal parallel to the lower diagonal to form a vector. \(\text{mean}(\cdot)\) represents taking the mean of each column vector of the matrix to form a row vector, \((\cdot)\) i represents the conjugate symbol, and \([\cdot;\cdot]\) represents concatenating the matrices by column. Then, use the vector \(u(i)\) to generate the Toeplitz matrix * * * ​ Step 4: Substitute T i (u) into the corresponding position of S 0-i to obtain S i ; Step Five: Calculate λ = ||S i - S0|| F , where ||·|| F is the Frobenius norm. If λ ≤ ε, then take out T i in the positive semi - definite matrix S i (u), and enter Step Six. Otherwise, set S0 = S i , and return to Step Two; Step 6: Take T i (u) and perform dimensionality reduction space smoothing - singular value decomposition on T i (u), and perform dimensionality reduction space smoothing processing on T i-ss (u). The specific process of its dimensionality reduction space smoothing is as follows: First, calculate the smoothing times N = M - (M sub - 1), where M sub is the number of sub - array elements. The covariance matrix corresponding to each sub - array is sub(T i (u, i)), i = 1, 2,..., N. T i-ss (u) = sum(sub(T i (u, i))) / N, where sub(T i (u, i)) = T i (u)(i:i + M sub - 1, i:i + M sub - 1), and sum(·) is the summation symbol; Step Seven: For T i-ss (u), perform singular value decomposition [U,∑,V]=svd(T i-ss (u)), where U N =U(:,K + 1:M - 1), where U is the matrix corresponding to the left singular vectors where Step 8: Convert the polynomial root z obtained by solving into the estimated value of the direction of arrival of the signal. The conversion formula is θ = -arcsin(angle(2×z))*180 / π, where angle(·) represents the phase angle corresponding to the complex number in [-π / 2, π / 2].

Images