CN113221059B

CN113221059B - Fast conjugate gradient direction finding algorithm without constructing covariance matrix

Info

Publication number: CN113221059B
Application number: CN202010721858.2A
Authority: CN
Inventors: 闫锋刚; 刘帅; 王军; 赫靓; 李享; 孟祥天
Original assignee: Harbin Institute of Technology Weihai
Current assignee: Harbin Institute of Technology Weihai
Priority date: 2020-07-24
Filing date: 2020-07-24
Publication date: 2023-01-17
Anticipated expiration: 2040-07-24
Also published as: CN113221059A

Abstract

The invention belongs to the technical field of array signal processing, and particularly relates to a fast conjugate gradient direction finding algorithm which can remarkably reduce the calculation complexity without constructing a covariance matrix under the condition of not losing the precision.

Description

Fast conjugate gradient direction finding algorithm without constructing covariance matrix

The technical field is as follows:

the invention belongs to the technical field of array signal processing, and particularly relates to a fast conjugate gradient direction finding algorithm which can remarkably reduce the computational complexity without constructing a covariance matrix under the condition of not losing the precision.

Background art:

estimation of the direction of arrival of a signal is an important research topic often encountered in applications such as radar, sonar, wireless communication, and passive positioning. The subspace algorithm represented by the multiple signal classification and rotation invariant subspace is provided, the leap of the traditional spatial spectrum estimation to super-resolution angle measurement is realized, but the huge calculation amount of the MUSIC algorithm and the lower estimation precision of the ESPRIT algorithm hinder the engineering progress of the super-resolution algorithm. To address this problem, the emergence of a Krylov subspace-based Conjugate Gradient (CG) algorithm was motivated.

The CG algorithm is one of the most classical Direction of Arrival (DOA) estimation methods based on Krylov subspace. Different from the traditional characteristic subspace algorithm, the CG algorithm utilizes the Krylov subspace to replace the signal subspace, and the effective algorithm of the rapid DOA estimation is realized by gradient iteration. Since the CG algorithm is at the beginning of iterationBoth the initial value setting and the iteration process need to calculate an Array Covariance Matrix (ACM), and the calculation complexity of the ACM is about O (M) ² N), if the calculation of ACM is avoided during the gradient iteration, the computational complexity will be greatly reduced, especially at large fast beat numbers. More importantly, most existing CG algorithms maintain good estimation performance while reducing computational complexity compared to the mucci algorithm. Therefore, a Krylov subspace-based CG technique should be chosen over the conventional approach.

The algorithm enriches the theoretical connotation of CG, but needs to estimate a signal subspace by means of iteration in practical engineering application. As is well known, the setting of the initial value of the iteration affects the efficiency and accuracy of the algorithm, and improper handling may result in the iteration failing to converge or converging to a wrong solution.

The invention content is as follows:

aiming at the problem of high complexity of a CG algorithm caused by the conventional calculation of the ACM, the invention provides a fast conjugate gradient direction finding algorithm which reduces the calculation complexity of the CG algorithm by replacing the ACM with a received data matrix through improving a wienerhoff equation and realizes the purpose of remarkably reducing the calculation complexity without constructing a covariance matrix under the condition of not losing the precision.

The invention can be achieved by the following measures:

a fast conjugate gradient direction finding algorithm without constructing a covariance matrix is characterized in that after a radiation source signal is received, a wiener Hough equation is improved by using a least square method to obtain a least square wiener Hough equation, the least square wiener Hough equation is a normal equation, a reference signal only related to received data is constructed by using the normal equation, an iteration initial value is set, conjugate gradient iteration is carried out, a signal subspace is estimated according to an information source number, and a signal wave arrival direction is obtained.

The method for improving the wiener Hough equation by using the least square method specifically comprises the following steps of:

step 1-1: wienerhoff equation RW = r _xd Simplifying the wienerhoff equation, R and R _xd Respectively expressed as:

wherein X is a matrix of received data signals, and has X = [ X (0), X (1), \ 8230;, X (M-1)]D is reference signal, d = [ d (0), d (1), \8230;, d (M-1)]；

The above formula is substituted into the wienerhoff equation

After simplification, X is obtained ^H W＝d ^H ；

Step 1-2: defining a least square wienerhoff equation as r + X ^H W＝d ^H Wherein r is a least square residual vector, and the least square wiener Hough equation is a normal equation;

step 1-3: to avoid computing the array covariance matrix R and to obtain an improved transition vector from the least squares wiener Hoff equation of v = X ^H p；

Step 1-4: defining the intermediate variable normal equation residual as q = Xr, where q is related only to the array received data and the reference signal.

The construction of the reference signal related to the received data only specifically comprises the following steps:

step 2-1: defining the reference signal as

Wherein a (theta) is a search steering vector, and theta e [ -90 DEG, 90 DEG]D is a reference signal;

step 2-2: the data matrix of the received signal of the k array element is expressed as

When theta = theta _i When i belongs to {1, \8230;, K }, the initial value q of the residual error of the normal equation _cg,0 Is shown as

It can be known that the normal equation residual error initial value q _cg,0 For linear combination of signal steering vectors, then q _cg,0 The K-dimensional Krylov subspace is a subset of the signal subspace, denoted as

span{Q _cg,K }≡span{A(θ)}，

Wherein

Q _cg,K ＝{q _cg,0 ,q _cg,1 ,…,q _cg,K-1 }；

Step 2-3: when theta is not equal to theta _i I belongs to {1, \8230;, K }, and the initial value q of the residual error of the normal equation _cg,0 Expressed as:

at this time, the normal equation residual vector is also a linear combination of the guide vectors, is combined by K signal vectors and a search guide vector, namely, a K + 1-dimensional Krylov subspace spanned by the normal equation residual vector at this time is also a subset of a subspace spanned by the signal guide vectors, and is expressed as

span{Q _cg,K+1 }≡span{A(θ),a(θ)}，

Wherein

Q _cg,K+1 ＝{q _cg,0 ,q _cg,1 ,…,q _cg,K-1 ,q _cg,K }。

The method specifically comprises the following steps of setting an iteration initial value and carrying out conjugate gradient iteration:

step 3-1: in the setting of the initial value of the conjugate gradient, the initial weight vector is: w is a ₀ =0, search vector initial value: p is a radical of ₁ ＝-r ₁ = -d, where d is a reference signal;

step 3-2: the least square residual vector initial value in the initial value setting is as follows: r is ₀ ＝d ^H -X ^H W, search vector initial value: p is a radical of ₁ ＝q ₀ ，γ ₀ ＝||q ₀ || ² The initial value of the normal equation residual vector is:q ₀ ＝Xr ₀ ；

step 3-3: when the iteration number is k =1,2, \8230, M, an intermediate transition variable is obtained: v. of _k ＝X ^H p _k And determining the step size: alpha is alpha _k ＝γ _i-1 /||v _k || ² To obtain a weight vector: w is a _k ＝w _k-1 +α _k p _k Using the least squares residual: r is _k ＝r _k-1 -α _k v _k To obtain the normal equation residual: q. q.s _k ＝Xr _k Determining a conjugate direction vector:

the method for estimating the signal subspace according to the information source number and obtaining the signal direction of arrival specifically comprises the following steps:

step 4-1: given the number of sources K, only K iterations are needed to obtain the Krylov subspace equivalent to the signal subspace: w is a group of _cg,K ＝{w _cg,0 ,w _cg,1 ,…,w _cg,K-1 The spectrum peak search function of the MUSIC algorithm is as follows:

where a (theta) is search guide vector, theta belongs to (-90 DEG, 90 DEG), and U _N Is a noise subspace;

step 4-2: if the source number is known a priori, the iteration number is selected to be K =1,2, \8230;, K, the new algorithm spectral peak search function is:

wherein I _M Is an identity matrix, W _cg,K Is a Krylov subspace equivalent to the signal subspace.

The invention receives the radiation source signal and expresses as follows: assuming that L mutually independent array elements form a Uniform Linear Array (ULA) at equal intervals d, K far-field narrow-band signals in space are considered to be incident to the array. Wherein, assuming that K is known a priori, d satisfies d ≦ λ 2 to avoid phase ambiguity, λ is the wavelength of the narrowband signal, and the array receiving radiation source signal is:

where a (θ) is an array flow pattern matrix of dimension L × K, S (t) is an incident signal vector of dimension K × 1, N (t) is an additive white gaussian noise vector of dimension L × 1, α (θ) is a column vector of a (θ), which can be expressed as:

the method realizes that the problem of high complexity caused by setting an iteration initial value and calculating the ACM in the iteration process is avoided by improving the wiener Hough equation and replacing the ACM with the received data matrix, and provides technical support for the engineering realization of the direction of arrival.

Description of the drawings:

FIG. 1 is a flow chart of the present invention.

FIG. 2 is a comparison graph of the spectral peak search of the present invention and the MUSIC algorithm, CG algorithm, where M =10, SNR =10dB, N =100, K =2, θ = ₁ ＝20°，θ ₂ ＝40°。

FIG. 3 is a graph of RMSE as a function of input SNR, where M =10, the signal-to-noise ratio is varied from-20 dB to 0dB, N =200, K =2, and θ, in 5dB steps ₁ ＝10°,θ ₂ ＝30°。

FIG. 4 is a graph comparing the computational efficiency of the present invention and the MUSIC algorithm, CG algorithm, RV-CG algorithm, where M is from 10 to 70, the step size is 10,SNR =10dB, N =200, K =2, θ = ₁ ＝10°,θ ₂ ＝30°。

The specific implementation mode is as follows:

aiming at the problem that initial values and an iteration process both need to calculate ACM to cause high algorithm complexity, the invention provides a method which does not need to construct a covariance matrix DOA estimator, replaces the ACM by a received data matrix through improving a wiener Hough equation to reduce the calculation complexity of a CG algorithm, and theoretical analysis and test results show that the method can remarkably reduce the calculation complexity without losing precision and can effectively avoid wrong solutions in the iteration process, thereby providing theoretical reference for actual engineering of CG.

As shown in the attached figure 1, the invention relates to a new conjugate gradient method based on a Krylov subspace without constructing a covariance matrix, which comprises the following specific steps:

a first step of receiving a radiation source signal using an antenna array, the first step comprising the steps of: (1) Assuming that L mutually independent array elements form a Uniform Linear Array (ULA) at equal intervals d, K far-field narrow-band signals in space are considered to be incident to the array. Wherein, assuming that K is known a priori, d is satisfied to be ≦ λ/2 to avoid phase ambiguity, λ being the wavelength of the narrowband signal. The array receives the radiation source signals as follows:

and a second step of improving the wienerhoff equation by using a least square method, wherein the second step comprises the following steps:

(1) Wienerhoff equation RW = r _xd The classical conjugate gradient solution process is:

first, setting an initial weight vector as w ₀ =0, initial value of search vector is p ₁ ＝-r ₁ = d, when k =1,2 \ 8230, M, the step size is determined as:

the transition vector is: v = Rp _k And updating the weight vector: w is a _k ＝w _k-1 +α _k p _k And obtaining a gradient vector: r is _k+1 ＝r _k +α _k Rd _k Determining a conjugate direction vector:

searching for a direction vector p _k Is R-conjugated, i.e. for all i ≠ j, there is

Memory matrix P _k ＝[p ₁ ,p ₂ ,…,p _k ]Is an M x k dimensional matrix composed of search direction vectors, then P _k Is column full rank. The conjugate gradient algorithm obtained by the analysis can be converged to a real solution in the k +1 gradient direction and P step at most by M steps of iteration _k Spanned k-dimensional space orthogonality:

after k iterations, the negative gradient direction vector r ₁ ，r ₂ ，…，r _k Are orthogonal to each other. Search direction p ₁ ,p ₂ ,…,p _k Spanned space, spanned space in gradient direction, matrix R and initial gradient vector R ₁ The Krylov subspace formed is equal:

thereby obtaining a signal subspace or a noise subspace, and thus performing DOA estimation.

(2) Reduction of the wienerhoff equation, R and R _xd Can be respectively expressed as:

wherein X is a matrix of received data signals, and has X = [ X (0), X (1), \ 8230;, X (M-1)]D is reference signal, and has d = [ d (0), d (1), \8230;, d (M-1)]。

Therefore, the above formula is introduced into the wienerhoff equation

After simplification, X is obtained ^H W＝d ^H 。

(3) At the moment, X in the wienerhoff equation does not meet the structure of the Hermite conjugate matrix, so in order to meet the requirement on conjugate gradient iteration in the Krylov subspace, the least square wienerhoff equation is defined to be r + X ^H W＝d ^H Wherein r is a least square residual vector, and the least square wiener Hough equation is a normal equation.

(4) In the conventional conjugate gradient algorithm, the transition vector is defined as v = Rp,

wherein p is a search vector, to avoid calculating the array covariance matrix R, and to obtain an improved transition vector according to the least squares wiener hough equation, v = X ^H p。

(5) Defining the residual of the normal equation of the intermediate variable as q = Xr, wherein q is only related to the array received data and the reference signal.

A third step of constructing a reference signal related to only the received data, the third step comprising the steps of:

(1) Defining the reference signal as

Where a (theta) is a search steering vector, theta ∈ [ -90 DEG, 90 DEG]D is a reference signal;

(2) The data matrix of the received signal of the kth array element is expressed as

From the above formula, it can be seen that the residual error initial value q of the normal equation _cg,0 For linear combination of signal steering vectors, then q _cg,0 The K-dimensional Krylov subspace is a signal subspaceIs represented as

span{Q _cg,K Is ≡ span { A (θ) }, where Q _cg,K ＝{q _cg,0 ,q _cg,1 ,…,q _cg,K-1 }。

(3) When theta is not equal to theta _i I ∈ {1, \8230;, K }, the initial value q of the residual error of the normal equation _cg,0 Expressed as:

the above equation shows that the normal equation residual vector is also a linear combination of the guide vectors, and is a combination of K signal vectors and the search guide vector, that is, the K + 1-dimensional Krylov subspace spanned by the normal equation residual vector is also a subset of the subspace spanned by the signal guide vectors, and is expressed as

span{Q _cg,K+1 Is { identical to ] span { A (theta), a (theta) }, wherein Q _cg,K+1 ＝{q _cg,0 ,q _cg,1 ,…,q _cg,K-1 ,q _cg,K }。

A fourth step of setting an iteration initial value and performing conjugate gradient iteration, the fourth step including the steps of:

(1) The initial weight vector in the classical conjugate gradient initial value setting is:

w ₀ =0, search vector initial value: p is a radical of formula ₁ ＝-r ₁ = -d, where d is a reference signal.

(2) The initial value of the least square residual vector in the initial value setting of the new algorithm is as follows: r is a radical of hydrogen ₀ ＝d ^H -X ^H W, search vector initial value: p is a radical of formula ₁ ＝q ₀ ，γ ₀ ＝||q ₀ || ² The initial value of the normal equation residual vector is: q. q.s ₀ ＝Xr ₀ 。

(3) When the number of iterations is k =1,2, \8230, M, an intermediate transition variable is obtained: v. of _k ＝X ^H p _k And determining the step size: alpha is alpha _k ＝γ _i-1 /||v _k || ² To obtain a weight vector: w is a _k ＝w _k-1 +α _k p _k Using the least squares residual: r is _k ＝r _k-1 -α _k v _k To obtain the normal equation residual: q. q.s _k ＝Xr _k Determining the conjugate direction vector:

and a fifth step of estimating a signal subspace according to the number of the information sources and obtaining DOA, wherein the fifth step comprises the following steps:

(1) Given that the number of the sources is K, only K iterations are needed to obtain a Krylov subspace equivalent to the signal subspace: w _cg,K ＝{w _cg,0 ,w _cg,1 ,…,w _cg,K-1 }。

(2) The spectrum peak search function of the MUSIC algorithm is as follows:

where a (theta) is search guide vector, theta belongs to (-90 DEG, 90 DEG), and U _N Is the noise subspace.

(3) If the source number is known a priori, the iteration times are selected to be K =1,2, \8230, K, and the new algorithm spectral peak search function is:

Example (b):

the performance of the invention is illustrated by example simulations as follows:

simulation conditions are as follows: assuming a ULA array type with 10 array elements and an array element spacing of d = λ 2, the direction of two incident signals is θ ₁ =10 ° and θ ₂ =30 °. To further evaluate the performance of the present invention, the number of monte carlo experiments was set to 500, and Root Mean Square Error (RMSE) was used as an evaluation index.

Simulation content and results:

simulation 1, the spectral peak search graphs of the present invention, the MUSIC algorithm and the CG algorithm are compared, and the results are shown in fig. 2.

Setting the initial reference signal to d in fig. 2 ₀ = Xa (θ)/| | Xa (θ) |, where θ ∈ (-90 °,90 °). The graph shows that the space spectrum peak of the least square CG algorithm is sharper than the space spectrum peaks of the classic CG algorithm and the MUSIC algorithm, the estimation performance is good, and the calculation complexity of the WCC-CG algorithm is lower than that of the classic CG algorithm, so that the WCC-CG algorithm can be better applied to practice.

Simulation 2, comparing the RMSE of the present invention with different algorithms with the input (Signal-to-noise ratio, SNR), shows the results in fig. 3.

It can be seen from fig. 3 that the estimation accuracy of the RV-CG algorithm and the estimation accuracy of the WCC-CG algorithm are both almost close to the estimation accuracy of the MUSIC algorithm, while the RV-CG algorithm and the WCC-CG algorithm are real-valued operations and do not need to calculate the covariance matrix of the received data and the eigenvalue decomposition, so that the computation amount of the RV-CG algorithm and the WCC-CG algorithm is relatively small under the condition that the accuracy is almost the same, and is between-20 dB and-7 dB, the accuracy of the WCC-CG algorithm is better than that of the other two algorithms, and the performance of the WCC-CG algorithm is poorer between-5 dB and 0dB, but is within an acceptable range. Through comparison of the three algorithms, the RV-CG algorithm and the WCC-CG algorithm have certain advantages, conditions required by actual engineering are further met, and the new algorithm is suitable for the conditions of low signal-to-noise ratio and small fast beat number.

And 3, simulating by using different array elements to compare the calculation efficiency of the method with different algorithms. By running MATLAB code in the same environment, the simulation results are given by the PC side with AMD Ryzen 5 3500U with Radon Vega Mobile Gfx 2.10GHz and 8.00 GBRAM. The equivalent evaluation of the calculation efficiency from the viewpoint of CPU time is performed, and the result is shown in fig. 4.

In fig. 4, the calculation time of the algorithm for searching spectral peaks is removed, and only the calculation of the reference signal and the process of obtaining the signal subspace through gradient iteration are recorded. It can be seen from the figure that the computation speed of the CG algorithm and the improved CG algorithm is faster than that of the MUSIC algorithm, and the greater the number of array elements is, the more obvious the advantages of the algorithm are. The WCC-CG algorithm has the fastest calculation speed, and the RV-CG algorithm and the CG algorithm are arranged in the second place.

The invention realizes that the reception data matrix is used for replacing ACM by improving the wienerhoff equation, avoids the problem of high complexity caused by setting an iteration initial value and calculating the ACM in the iteration process, and provides technical support for the engineering realization of the direction of arrival.

Claims

1. A fast conjugate gradient direction-finding algorithm without constructing a covariance matrix is characterized in that after a radiation source signal is received, a wienerhoff equation is improved by using a least square method to obtain a least square wienerhoff equation, the least square wienerhoff equation is a normal equation, a reference signal only related to received data is constructed by using the normal equation, an iteration initial value is set, conjugate gradient iteration is carried out, a signal subspace is estimated according to the number of information sources, and the direction of arrival of the signal is obtained; the method for improving the wienerhoff equation by using the least square method specifically comprises the following steps:

wherein X is a matrix of received data signals, and has X = [ X (0), X (1), \ 8230;, X (M-1)]D is reference signal, and has d = [ d (0), d (1), \8230;, d (M-1)]；

The above formula is introduced into the wienerhoff equation, one

After simplification, X is obtained ^H W＝d ^H ；

Step 1-2: defining a least square wienerhoff equation as r + X ^H W＝d ^H Wherein r is a least square residual vector, and the least square wienerhoff equation is a normal equation;

Step 1-4: defining the residual error of the intermediate variable normal equation as q = Xr, wherein q is only related to array received data and a reference signal;

the construction of the reference signal relating only to the received data comprises in particular the steps of:

step 2-1: defining the reference signal as

step 2-2: the data matrix of the received signal of the kth array element is expressed as

span{Q _cg,K }≡span{A(θ)}，

Wherein

Q _cg,K ＝{q _cg,0 ,q _cg,1 ,…,q _cg,K-1 }；

span{Q _cg,K+1 }≡span{A(θ),a(θ)}，

Wherein

Q _cg,K+1 ＝{q _cg,0 ,q _cg,1 ,…,q _cg,K-1 ,q _cg,K }。

2. The fast conjugate gradient direction-finding algorithm without constructing covariance matrix as claimed in claim 1, wherein the setting of iteration initial value and conjugate gradient iteration specifically comprises the following steps:

step 3-2: the least square residual vector initial value in the initial value setting is as follows: r is a radical of hydrogen ₀ ＝d ^H -X ^H W, search vector initial value: p is a radical of ₁ ＝q ₀ ，γ ₀ ＝||q ₀ || ² The initial value of the normal equation residual vector is: q. q of ₀ ＝Xr ₀ ；

Step 3-3: when the number of iterations is k =1,2, \8230, M, an intermediate transition variable is obtained: v. of _k ＝X ^H p _k And determining the step size: alpha (alpha) ("alpha") _k ＝γ _i-1 /||v _k || ² And obtaining a weight vector: w is a _k ＝w _k-1 +α _k p _k Using the least squares residual: r is a radical of hydrogen _k ＝r _k-1 -α _k v _k Obtaining a normal equation residual: q. q.s _k ＝Xr _k Determining the conjugate direction vector:

3. the fast conjugate gradient direction-finding algorithm without constructing covariance matrix as claimed in claim 1, wherein said estimating signal subspace according to source number and obtaining signal direction of arrival specifically comprises the following steps:

step 4-1: given that the number of the sources is K, only K iterations are needed to obtain a Krylov subspace equivalent to the signal subspace: w is a group of _cg,K ＝{w _cg,0 ,w _cg,1 ,…,w _cg,K-1 The spectrum peak search function of the MUSIC algorithm is as follows:

step 4-2: if the source number is known a priori, the iteration times are selected to be K =1,2, \8230, K, and the new algorithm spectral peak search function is:

4. A fast conjugate gradient direction finding algorithm without constructing covariance matrix as claimed in claim 1 wherein the received radiation source signal is represented as follows: assuming that L mutually independent array elements form a Uniform Linear Array (ULA) by equal spacing d, considering that K far-field narrow-band signals exist in a space and enter the array, wherein K is known a priori, d satisfies d ≦ λ 2 to avoid phase ambiguity, λ is the wavelength of the narrow-band signals, and the array receiving radiation source signals are as follows:

where a (θ) is an array flow pattern matrix of dimension L × K, S (t) is an incident signal vector of dimension K × 1, N (t) is an additive white gaussian noise vector of dimension L × 1, α (θ) is a column vector of a (θ), and is expressed as: