Target orientation estimation algorithm based on covariance matrix imaginary part reconstruction

Abstract

The invention provides a target orientation estimation algorithm based on covariance matrix imaginary part reconstruction for solving the problems that the underwater environment is complex and weak signals are difficult to detect. The algorithm includes calculating a covariance matrix from array received data; extracting an imaginary part of the covariance matrix, and reconstructing a new covariance matrix according to the characteristic value and the corresponding characteristic vector of the imaginary part; initializing a variable factor and a noise power estimation value; calculating initial values of a spatial orientation spectrum and sparseness thereof; on the basis of a sparse approximate minimum variance principle, iteratively calculating the spatial orientation spectrum until convergence; normalization processing is carried out on the iteration space azimuth spectrum, angles corresponding to the first K maximum values are target azimuth estimation values, and K is the real number of target signals. According to the method, the covariance matrix is reconstructed by extracting the covariance matrix imaginary part, so that the influence of noise on weak signals is weakened; and on the other hand, the estimation performance of the weak signal and the sparseness of the azimuth spectrum are considered by adaptively adjusting the variable factor, and the resolution capability of the weak target is enhanced.

Description

Technical Field

[0001] The invention relates to a signal processing method, in particular to a target direction estimation algorithm based on covariance matrix imaginary part reconstruction. Background Art

[0002] Signal detection technology in complex background environments is an important research topic in the field of signal processing. In engineering applications, under the combined effects of internal waves, currents, and swells, the underwater background noise environment is very complex, and the signal-to-noise ratio of the signal received by the array is generally not very large. When the interference energy is large, the accuracy of the azimuth estimation of weak targets will be affected, and may even not be reflected in the spatial azimuth spectrum at all.

[0003] In actual working scenarios, the noise received by different array elements is uncorrelated or has a small correlation coefficient. From the spatial distribution characteristics, it can be seen that the underwater environmental noise source approximately presents the characteristics of symmetrical distribution of azimuth and power statistics. Therefore, most of the noise components exist in the real part of the diagonal elements of the covariance matrix of the array receiving data. Based on this characteristic of the data covariance matrix, some people have proposed an algorithm to suppress noise by taking the diagonal elements of the data covariance matrix to zero, but this type of algorithm will lead to the loss of signal components in the diagonal elements, which will have a certain impact on the target's azimuth estimation. The basic idea of ​​the diagonal unloading method is to add a unloading amount or multiply a unloading coefficient to the diagonal elements of the array's data covariance matrix to reduce the diagonal elements, thereby weakening the impact of noise on the signal. However, it is difficult to accurately give the unloading amount or unloading coefficient in this method.

[0004] Since noise can be divided into symmetric noise and asymmetric noise, and symmetric noise only exists in the real part of the diagonal elements of the covariance matrix, the influence of symmetric noise on the array receiving data can be theoretically eliminated by extracting the imaginary part of the covariance matrix. Based on this idea, some researchers proposed to estimate the target direction by extracting the imaginary part of the covariance matrix and reconstructing its real part through the particle swarm algorithm. However, the particle swarm algorithm is sensitive to initial values ​​and is prone to fall into local optimal solutions during the calculation process. In addition, the computational complexity of this method increases rapidly with the increase in the number of targets, so this method still has room for improvement.

[0005] Therefore, it is urgent to propose corresponding solutions to the problem that weak targets are difficult to detect in complex underwater noise environments. Summary of the invention

[0006] The present invention aims to solve at least one of the technical problems existing in the prior art. To this end, the present invention proposes a target azimuth estimation algorithm based on reconstruction of the imaginary part of the covariance matrix. By reconstructing a new covariance matrix, the interference of symmetric noise on the signal is weakened. For the problem of weak target detection, the present invention compromises the sparsity of the algorithm and the estimation performance of weak signals. By adaptively adjusting the variable factors, the cancellation of weak targets is avoided while ensuring sparsity, which solves the problem of difficulty in detecting weak targets in complex underwater noise environments to a certain extent.

[0007] To achieve the above object, the present invention proposes a target direction estimation algorithm based on covariance matrix imaginary part reconstruction, comprising the following steps:

[0008] S1 calculates the covariance matrix based on the array received data;

[0009] S2 takes the imaginary part of the covariance matrix obtained in step S1, performs eigendecomposition, and extracts its non-zero eigenvalues ​​and corresponding eigenvectors;

[0010] S3 reconstructs a new covariance matrix according to the eigenvalues ​​and corresponding eigenvectors obtained in step S2;

[0011] S4 initializes variable factors and noise power estimation values;

[0012] S5 calculates the initial value of the spatial orientation spectrum and its sparsity;

[0013] S6 iteratively calculates the spatial orientation spectrum according to the sparse approximate minimum variance principle until convergence;

[0014] S7 normalizes the iterative spatial orientation spectrum to obtain the angles corresponding to the first K maximum values ​​thereof, where the angles are target orientation estimates, wherein K is the actual number of target signals.

[0015] Furthermore, the calculation of the covariance matrix is:

[0016] Assume that the time domain data received by the array is Y and the number of snapshots is N snap , then the covariance matrix R = Y*Y H / N snap , (·) H is the conjugate transpose operation.

[0017] Furthermore, the eigenvalues ​​of the imaginary part of the covariance matrix and the corresponding eigenvectors are extracted, namely:

[0018] Take the imaginary part of the covariance matrix R and perform eigendecomposition to obtain 2K non-zero eigenvalues ​​λ k and the corresponding eigenvector u k .

[0019] Furthermore, the reconstruction of the new covariance matrix includes the following steps:

[0020] S301 uses the feature vector u k Construct the matrix U = [u1,u2,K,u 2K ], using the eigenvalue λ k Construct a diagonal matrix Λ=diag{[λ1,λ2,K,λ 2K ]}, diag{·} returns a diagonal matrix containing the elements in brackets on the main diagonal or extracts a vector containing the diagonal elements of the matrix in brackets;

[0021] S303 screens the diagonal elements of the diagonal matrix Λ, sets the elements less than zero to zero, and obtains a new diagonal matrix in max{·} means extracting the maximum value of the elements in the brackets;

[0022] S305 constructs a new covariance matrix, and obtains

[0023] Furthermore, the initialization of the variable factor and the noise power estimation value comprises the following steps:

[0024] S401 The variable factor is α0, and the initial value of the variable factor α0 is set to 1;

[0025] S403 Calculate the initial value of the noise power estimate Where trace(·) is the trace of the matrix in the brackets, and M is the total number of array elements.

[0026] Furthermore, the calculation of the initial value of the spatial orientation spectrum and its sparsity comprises the following steps:

[0027] S501 determines a scanning angle set θ, where the length of θ is N. For any angle θ i ∈Θ, the expression of the array manifold vector is:

[0028] a(θ i )=[1,e j2πp ,Ke j2π(M-1)p ] T ;

[0029] in d is the array element spacing, λ is the signal wavelength, j is the imaginary unit, (·) T is the transpose operation;

[0030] The initial value P0 of the S503 spatial orientation spectrum is expressed as:

[0031]

[0032] The array manifold matrix A(θ) = [a(θ1), a(θ2), K, a(θ N )],θ i ∈Θ;

[0033] The initial value β0 of S505 sparsity is expressed as:

[0034] β0=||P0 / max(P0)||1;

[0035] Where ||·||1 is the 1-norm of the calculated vector.

[0036] Furthermore, the iterative calculation of the spatial orientation spectrum based on the sparse approximate minimum variance principle includes the following steps:

[0037] S601 updates the covariance matrix. The covariance matrix obtained in the jth iteration is The calculation formula is:

[0038]

[0039] Among them I M is an M×M identity matrix, P j-1 and p n,j-1 is the spatial orientation spectrum and noise power estimation value obtained in the j-1th iteration;

[0040] S603 Update the spatial orientation spectrum, the spatial orientation spectrum P obtained in the jth iteration j The calculation formula is:

[0041]

[0042] in(·) -1 To calculate the inverse of the matrix in brackets;

[0043] S605 updates the noise power estimation value, and the noise power estimation value p obtained in the jth iteration is n,j The calculation formula is:

[0044]

[0045] S607 updates the sparsity, and the sparsity β obtained in the jth iteration j The calculation formula is:

[0046] β j =||P j / max(P j )||1;

[0047] S609 updates the variable factor, and the variable factor α obtained in the jth iteration j The calculation formula is:

[0048]

[0049] Where abs(·) is the absolute value operation;

[0050] S611 Determine the spatial orientation spectrum P j Whether it converges, calculate the judgment factor Q = abs(β j -β j-1 ), when Q≤0.005, the iterative calculation ends and the iterative spatial orientation spectrum is obtained

[0051] Furthermore, the normalization process of the iterative spatial orientation spectrum is:

[0052] Normalized spatial orientation spectrum The calculation formula is:

[0053]

[0054] The angles corresponding to the first K maximum values ​​of are the estimated target orientations.

[0055] Compared with the prior art, the present invention has the following beneficial effects:

[0056] The present invention eliminates the real terms of the noise covariance matrix contained in the diagonal elements of the data covariance matrix by extracting the imaginary part of the covariance matrix, thereby weakening the influence of noise on weak signals and improving the performance of weak target azimuth estimation in complex noise environments; on the other hand, by adaptively adjusting variable factors, the estimation performance of weak signals and the sparsity of azimuth spectra are taken into account, thereby enhancing the resolution capability of weak targets. BRIEF DESCRIPTION OF THE DRAWINGS

[0057] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0058] Figure 1 It is a schematic diagram of the process of the present invention;

[0059] Figure 2 It is the spatial orientation spectrogram obtained by respectively processing the covariance matrix and its imaginary part by different algorithms in the specific embodiment. DETAILED DESCRIPTION

[0060] The technical solution of the present invention will be clearly and completely described below in conjunction with the embodiments. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention. Specific embodiment 1

[0062] like Figure 1 As shown, according to an embodiment of the present invention, the present invention proposes a target direction estimation algorithm based on reconstruction of the imaginary part of the covariance matrix, and the technical solution adopted includes the following steps:

[0063] Step 1: Use an M-element uniform linear array located in the xoy plane as the receiving array, with the element spacing equal to half the signal wavelength. Calculate the covariance matrix based on the array received data. Let the time domain data received by the array be Y, and the number of snapshots be N. snap , then the covariance matrix R = Y*Y H / N snap , (·) H is the conjugate transpose operation.

[0064] Step 2: Take the imaginary part of the covariance matrix R in step 1 and perform eigendecomposition to obtain 2K non-zero eigenvalues ​​λ k and the corresponding eigenvector u k .

[0065] Step 3: Use the eigenvalue λ obtained in step 2 k and the corresponding eigenvector u k , reconstruct the new covariance matrix, the specific implementation steps are as follows:

[0066] 1) Step 3-1: Using the feature vector u k Construct the matrix U = [u1,u2,K,u 2K ]; using the eigenvalue λ k Construct a diagonal matrix Λ=diag{[λ1,λ2,K,λ 2K ]}, diag{·} returns a diagonal matrix containing the elements in brackets on the main diagonal or extracts a vector containing the diagonal elements of the matrix in brackets;

[0067] 2) Step 3-2: Filter the diagonal elements of the diagonal matrix Λ, set the elements less than zero to zero, and obtain a new diagonal matrix in max{·} means extracting the maximum value of the elements in the brackets;

[0068] 3) Step 3-3: Construct a new covariance matrix

[0069] Step 4: Initialize the variable factors and noise power estimation values. The specific implementation steps are as follows:

[0070] 1) Step 4-1: Set the initial value of the variable factor α0 to 1;

[0071] 2) Step 4-2: Calculate the initial value of the noise power estimate Where trace(·) is the trace of the matrix in the brackets, and M is the total number of array elements.

[0072] Step 5: Calculate the initial value of the spatial orientation spectrum and its sparsity. The specific implementation steps are as follows:

[0073] 1) Step 5-1: Determine the scanning angle set θ, where the length of θ is N; for any angle θ i ∈Θ, the expression of the array manifold vector is:

[0074] a(θ i )=[1,e j2πp ,Ke j2π(M-1)p ] T ;

[0075] in d is the array element spacing, λ is the signal wavelength, j is the imaginary unit, (·) T is the transpose operation;

[0076] 2) Step 5-2: The initial value P0 of the spatial orientation spectrum is expressed as:

[0077]

[0078] The array manifold matrix A(θ) = [a(θ1), a(θ2), K, a(θ N )],θ i ∈Θ;

[0079] 3) Step 5-3: The initial value β0 of sparsity is expressed as:

[0080] β0=||P0 / max(P0)||1;

[0081] Where ||·||1 is the 1-norm of the calculated vector.

[0082] Step 6: Based on the sparse approximate minimum variance principle, iteratively calculate the spatial orientation spectrum until convergence. The specific implementation steps are as follows:

[0083] 1) Step 6-1: Update the covariance matrix. The covariance matrix obtained in the jth iteration is The calculation formula is:

[0084]

[0085] Among them I M is an M×M identity matrix, P j-1 and p n,j-1 is the spatial orientation spectrum and noise power estimation value obtained in the j-1th iteration;

[0086] 2) Step 6-2: Update the spatial orientation spectrum. The spatial orientation spectrum P obtained in the jth iteration is j The calculation formula is:

[0087]

[0088] in(·) -1 To calculate the inverse of the matrix in brackets;

[0089] 3) Step 6-3: Update the noise power estimate. The noise power estimate p obtained in the jth iteration is n,j The calculation formula is:

[0090]

[0091] 4) Step 6-4: Update the sparsity, the sparsity β obtained in the jth iteration j The calculation formula is:

[0092] β j =||P j / max(P j )||1;

[0093] 5) Step 6-5: Update the variable factor, the variable factor α obtained in the jth iteration j The calculation formula is:

[0094]

[0095] Where abs(·) is the absolute value operation;

[0096] 6) Step 6-6: Determine the spatial orientation spectrum P j Whether it converges, calculate the judgment factor Q = abs(β j -β j-1 ), when Q≤0.005, the iterative calculation ends and the iterative spatial orientation spectrum is obtained

[0097] Step 7: Normalize the iterative spatial orientation spectrum; Normalize the spatial orientation spectrum value The calculation formula is:

[0098]

[0099] The angles corresponding to the first K maximum values ​​of are the estimated target orientations.

[0100] In order to further verify the actual effect of the technical solution of this application, this application also provides a simulation experiment case.

[0101] This simulation experiment uses a 16-element uniform linear array as the receiving array, and the array element spacing is half the signal wavelength. Assume that K = 2 far-field narrowband signals are incident on the array, with incident angles of 50° and 120°; the signal frequency is 500Hz, the sampling frequency is five times the signal frequency, and the sampling time is 0.2 seconds; the signal-to-noise ratio is -15dB and 5dB; using the conventional beamforming algorithm and the algorithm proposed in this invention, the normalized spatial orientation spectrum is as follows: Figure 2 shown.

[0102] from Figure 2 It can be seen that when there are multiple signals of different strengths in space, the incident azimuth of weak signals in the obtained spatial orientation spectrum is difficult to identify using conventional algorithms to process the covariance matrix or to process the imaginary part of the covariance matrix alone. In contrast, the present invention successfully estimates the azimuth of all target signals in space and demonstrates strong angular resolution.

[0103] It should be noted that in the above embodiments, the description of each embodiment has its own emphasis. For parts that are not described in detail in a certain embodiment, please refer to the relevant description of other embodiments. It should be understood by those skilled in the art that the embodiments of the present invention can be provided as methods, systems or computer program products. Therefore, the present invention can take the form of a complete hardware embodiment, a complete software embodiment or an embodiment combining software and hardware. Moreover, the present invention can take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0104] The present invention is described with reference to flowcharts and / or block diagrams of methods and computer program products according to embodiments of the present invention. It should be understood that each process and / or block in the flowchart and / or block diagram and the combination of processes and / or blocks in the flowchart and / or block diagram can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded computer or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the process in the flowchart and / or block diagram. Figure 1 A process or multiple processes and / or boxes Figure 1These computer program instructions can also be stored in a computer-readable memory that can guide a computer or other programmable data processing device to work in a specific way, so that the instructions stored in the computer-readable memory produce a product including an instruction device, which implements the functions specified in the process. Figure 1 A process or multiple processes and / or boxes Figure 1 These computer program instructions can also be loaded onto a computer or other programmable data processing device, so that a series of operation steps are executed on the computer or other programmable device to produce a computer-implemented process, so that the instructions executed on the computer or other programmable device provide for implementing the process in the process. Figure 1 A process or multiple processes and / or boxes Figure 1 The steps for the functions specified in one or more boxes.

[0105] Although preferred embodiments of the present invention have been described, additional changes and modifications may occur to these embodiments once those skilled in the art understand the basic inventive concepts.

[0106] Obviously, those skilled in the art can make various changes and modifications to the present invention without departing from the spirit and scope of the present invention. Thus, if these modifications and variations of the present invention belong to the scope of the present invention and its equivalent technologies, the present invention is also intended to include these changes and variations.

Claims

1. A target direction estimation algorithm based on reconstruction of the imaginary part of the covariance matrix, characterized in that: The following steps are involved: S1 calculates the covariance matrix based on the array received data; S2 takes the imaginary part of the covariance matrix obtained in step S1, performs eigendecomposition, and extracts its non-zero eigenvalues ​​and corresponding eigenvectors; S3 reconstructs a new covariance matrix according to the eigenvalues ​​and corresponding eigenvectors obtained in step S2; S4 initializes variable factors and noise power estimation values; S5 calculates the initial value of the spatial orientation spectrum and its sparsity; S6 iteratively calculates the spatial orientation spectrum according to the sparse approximate minimum variance principle until convergence; S7 normalizes the iterative spatial orientation spectrum to obtain the angles corresponding to the first K maximum values ​​thereof, where the angles are target orientation estimates, wherein K is the actual number of target signals.

2. A target direction estimation algorithm based on covariance matrix imaginary part reconstruction according to claim 1, characterized in that, The calculation of the covariance matrix described in step S1 assumes that the time domain data received by the array is Y and the number of snapshots is N. snap , then the covariance matrix R = Y*Y H / N snap , (·) H is the conjugate transpose operation.

3. A target direction estimation algorithm based on covariance matrix imaginary part reconstruction according to claim 2, characterized in that, Extract the eigenvalues ​​of the imaginary part of the covariance matrix and the corresponding eigenvectors, take the imaginary part of the covariance matrix R, perform eigendecomposition, and obtain 2K non-zero eigenvalues ​​λ k and the corresponding eigenvector u k .

4. A target direction estimation algorithm based on covariance matrix imaginary part reconstruction according to claim 1, characterized in that, The reconstruction of the new covariance matrix described in step S3 comprises the following steps: S301 uses the feature vector u k Construct the matrix U = [u1,u2,K,u 2K ], using the eigenvalue λ k Construct a diagonal matrix Λ=diag{[λ1,λ2,K,λ 2K ]}, diag{·} returns a diagonal matrix containing the elements in brackets on the main diagonal or extracts a vector containing the diagonal elements of the matrix in brackets; S303 screens the diagonal elements of the diagonal matrix Λ, sets the elements less than zero to zero, and obtains a new diagonal matrix in max{·} means extracting the maximum value of the elements in the brackets; S305 constructs a new covariance matrix, and obtains 5. A target direction estimation algorithm based on covariance matrix imaginary part reconstruction according to claim 4, characterized in that, The initialization of the variable factor and the noise power estimation value comprises the following steps: S401 The variable factor is α0, and the initial value of the variable factor α0 is set to 1; S403 Calculate the initial value of the noise power estimate Where trace(·) is the trace of the matrix in the brackets, and M is the total number of array elements.

6. A target direction estimation algorithm based on covariance matrix imaginary part reconstruction according to claim 1, characterized in that: The calculation of the initial value of the spatial orientation spectrum and its sparsity described in step S5 includes the following steps: S501 determines a scanning angle set θ, where the length of θ is N. For any angle θ i ∈Θ, the expression of the array manifold vector is: a(θ i )[1,e j2πp ,Is j2π(M-1)p ] T 100. in d is the array element spacing, λ is the signal wavelength, j is the imaginary unit, (·) T is the transpose operation; The initial value P0 of the S503 spatial orientation spectrum is expressed as: The array manifold matrix A(θ) = [a(θ1), a(θ2), K, a(θ N )],θ i ∈Θ; The initial value β0 of S505 sparsity is expressed as: β0=||P0 / max(P0)||1; in ||·||1 is the 1-norm of the calculated vector.

7. A target direction estimation algorithm based on covariance matrix imaginary part reconstruction according to claim 6, characterized in that: The iterative calculation of the spatial orientation spectrum based on the sparse approximate minimum variance principle in step S6 includes the following steps: S601 updates the covariance matrix. The covariance matrix R obtained in the jth iteration is j The calculation formula is: Among them I M is an M×M identity matrix, P j-1 and p n,j-1 is the spatial orientation spectrum and noise power estimation value obtained in the j-1th iteration; S603 Update the spatial orientation spectrum, the spatial orientation spectrum P obtained in the jth iteration j The calculation formula is: in(·) -1 To calculate the inverse of the matrix in brackets; S605 updates the noise power estimation value, and the noise power estimation value p obtained in the jth iteration is n,j The calculation formula is: S607 updates the sparsity, and the sparsity β obtained in the jth iteration j The calculation formula is: b j =||P j / max(P j )||1; S609 updates the variable factor, and the variable factor α obtained in the jth iteration j The calculation formula is: Where abs(·) is the absolute value operation; S611 Determine the spatial orientation spectrum P j Whether it converges, calculate the judgment factor Q = abs(β j -β j-1 ), when Q≤0.005, the iterative calculation ends and the iterative spatial orientation spectrum is obtained 8. A target direction estimation algorithm based on covariance matrix imaginary part reconstruction according to claim 7, characterized in that: Step S7 is to iterate the normalization process of the spatial orientation spectrum to normalize the spatial orientation spectrum value. The calculation formula is: The angles corresponding to the first K maximum values ​​of are the estimated target orientations.

Images