CN113325365A - Quaternion-based coherent signal two-dimensional DOA estimation method - Google Patents

Abstract

The invention provides a quaternion-based coherent signal two-dimensional DOA estimation method, which comprises the following steps: establishing an observation signal complex model, establishing an observation signal quaternion model, constructing a correlation matrix, simplifying quaternion cross-correlation vectors, realizing information source decorrelation through correlation matrix smoothing, estimating a signal pitch angle by using a propagation operator method, estimating a propagation operator, estimating an incident signal horizontal azimuth angle, and reconstructing an array response matrix. The present invention solves the problem of how to make full use of the array data while eliminating the additional additive noise introduced.

Description

Quaternion-based coherent signal two-dimensional DOA estimation method

Technical Field

The invention relates to accurate estimation of the direction of arrival of a plurality of coherent signal sources by using a vector sensor array, in particular to a quaternion-based coherent signal two-dimensional DOA estimation method, belonging to the field of array signal processing.

Background

Compared with the traditional scalar sensor, the acoustic vector sensor can measure the sound pressure of a certain point in a sound field and simultaneously provide the three-dimensional particle vibration velocity information of the point, so that the sound field information can be more comprehensively and fully acquired, and therefore the acoustic vector sensor plays an important role in the field of underwater acoustic array signal processing in recent years. In the traditional method, output signals of the vector hydrophone are modeled in a complex domain, and all components of the signals are spliced into a long vector form, so that the orthogonal relation among all components of the output signals of the vector sensor is destroyed, and the structural information of the array is not completely utilized. In recent years, therefore, many scholars have begun to attempt to model vector signals using quaternions and have achieved numerous achievements.

Accurate estimation (DOA) of directions of arrival of a plurality of signal sources is one of the main contents studied in the field of array signal processing, and the fields of sonar, radar, seismic detection, communication system and the like have the shadow of the related technology, and one of the representatives of the technology is a subspace method. The traditional subspace method obtains a signal subspace or a noise subspace by performing characteristic decomposition on a data covariance matrix, and then obtains DOA estimation of a signal according to orthogonality (MUSIC) of the noise subspace and the signal subspace or rotation invariance (ESPRIT) of the signal subspace in a certain array layout form, wherein the characteristic decomposition of a high-dimensional data covariance matrix is time-consuming and labor-consuming, and a propagation operator method is generated to solve the problem. According to the method, firstly, a covariance matrix is subjected to subarray division, then, a noise subspace or a signal subspace can be obtained through linear operation, then, the arrival angle of a signal is obtained in a mode similar to an MUSIC or ESPRIT algorithm, but the performance of the algorithm under low signal-to-noise ratio is not ideal due to the fact that a propagation operator method does not fully utilize the spatial domain characteristic of noise. Therefore, many algorithms (e.g., sumpe, CODE) based on array cross-correlation and applying propagation operator techniques have been proposed, which successfully suppress additive noise by array cross-correlation, thereby providing guarantees for good performance of the propagation operator method.

However, the disadvantage of removing additive noise by array cross-correlation is that: (1) the step of sub-array partitioning the array loses array aperture; (2) efficient use of array data cannot be achieved by only considering the cross-correlation between each subarray and ignoring the autocorrelation information within the subarray.

Therefore, how to make full use of the array data (i.e., taking into account both the cross-correlation of the array and the autocorrelation information of the array elements) and eliminate the extra additive noise (introduced by the array element autocorrelation) is a key breakthrough to further gain the performance of such algorithms.

Disclosure of Invention

In order to overcome the defects of the prior art, the invention provides a brand new thought for modeling signals in a quaternion domain, and aims to solve the problem of how to eliminate the introduced extra additive noise while fully utilizing array data. It should be noted that the method has a certain universality, and the excellent effect obtained by combining the method with the SUMWE algorithm (AQ-SUMWE) further proves the feasibility and effectiveness of the method.

The purpose of the invention is realized as follows: the method comprises the following steps:

the method comprises the following steps: establishing an observation signal complex model;

step two: establishing an observation signal quaternion model; two quaternion signal models were constructed as follows:

model one:

z(t)＝p(t)+x(t)j+y(t)j

＝A(s(t)+Γ_xs(t)j+Γ_ys(t)j)+n(t)

model two:

step three: constructing a correlation matrix: considering the cross-correlation between the quaternion model of the m-th array element output signal and the quaternion model of the entire array output signal, there are:

will r is_mDeployed, in the following form:

step four: quaternion cross-correlation vectorial reduction, r_mThe final form of (a) is:

wherein the content of the first and second substances,

R_s＝E{s(t)s^H(t) } is a signal covariance matrix;

step five: the information source decoherence is realized through the correlation matrix smoothing;

step six: estimating the signal pitch angle by using a propagation operator method:

step seven: estimating a propagation operator, and dividing R as follows:

wherein R is₁，R₂Are respectively the front K line and the back M-2K +1 line of R₂＝P₀R₁

Thus, the propagation operator can be expressed as:

estimating propagation factors from sample data

Pitch angle

Can be selected from

K characteristic values of

Obtaining:

step eight: the azimuth of the incident signal level is estimated,

concatenating x (t) and y (t) given in (a) to form the following array output vector z' (t):

the following 2M × M cross-correlation matrix R' is introduced:

wherein

And n_p(t) are spatially independent of each other, i.e.

Thus statistically, the noise in R' has been eliminated; from R'_(1:M,:)、R′_(M+1:2M,:)Respectively represent the front M rows and the rear M rows of R', will

Substituting, there are:

w_x,k,w_y,kk is W, 1,2, K is W, respectively_x,W_yThe kth line of (1); from the relationship: u. of_k tanθ_k＝v_kIt is known that

tanθ_kw_x,k＝w_y,k

Step nine: reconstructing an array response matrix, with reference to a (phi) ═ a (phi)₁),a(φ₂),...,a(φ_K)]And use

K1, 2, …, K to reconstruct the matrix a, i.e.

W_x,W_yThe estimated values of (c) are:

the azimuth angle of the kth signal can be estimated as:

therefore, the pitch angle and the azimuth angle of the incident signal are obtained, and the pitch angle and the azimuth angle are automatically matched.

The invention also includes such structural features:

1. the first step is specifically as follows: assuming that the array is a linear array of vector sensors arranged at equal intervals of half wavelength along the z-axis, and a single vector sensor comprises a sound pressure sensor and two particle vibration velocity sensors respectively pointing to the x-axis and the y-axis, the manifold vector of the vector sensor can be written as follows:

wherein, theta_kAnd phi_kHorizontal and pitch angles, u, of the kth signal, respectively_kAnd v_kThe direction cosine corresponds to a particle vibration velocity sensor pointing to an x axis and a y axis respectively; assuming that K far-field narrow-band signals are incident on the sensor array, receiving signal models of the sound pressure sensor array, the x-axis-oriented particle velocity sensor array and the y-axis-oriented particle velocity sensor array can be respectively established in a complex number domain, and the following steps are performed in sequence:

wherein A (phi) is [ a (phi) ]₁),a(φ₂),...,a(φ_K)]Is a manifold of an array of acoustic pressure sensors,

for the steering vector corresponding to the k-th signal,

is space domain phase factor between sensors, d is array element spacing, lambda is signal wavelength, and signal vector s (t) is [ s ]₁(t),...,s_K(t)]^TIs a complex gaussian random process, beta ═ beta₁,…,β_K]^TIs a complex attenuation coefficient of which beta ₁1 and beta_k≠0,k＝2,3,…K，n_lWhere l is p, x, y are noise vectors that are uncorrelated with the signal and independent of each otherQuantity, Γ_x＝diag{cosθ₁sinφ₁,…,cosθ_K sinφ_K}，Γ_y＝diag{sinθ₁sinφ₁,...,sinθ_K sinφ_K}。

2. The fifth step is specifically as follows: considering K coherent signals, we pass through the correlation matrix r_mSmoothing to construct a matrix

Can be expressed as:

wherein A is₀The first M-K +1 line containing A,

is the selection of the matrix or matrices,

and

is two zero matrices, Π_m＝[ψ_m,Φ¹ψ_m,...,Φ^K-1ψ_m],

It can be found that,

is K, so

Can be used to estimate the direction of the K coherent signals;

the following correlation matrix R is defined for subsequent processing:

wherein pi ═ pi₁,Π₂,...,Π_M]。

3. The sixth step is specifically as follows: propagation operator method based on pair matrix A₀Is divided into sub-arrays, i.e.

Wherein, because A₁Is a full rank matrix of K rows, A₂Is a matrix of M-2K +1 rows; a. the₀Is linearly independent, then A₀The remaining rows of (a) may be represented as a₀Linear combination of the first K rows; propagation matrix P₀Is defined as an AND matrix A₁、A₂The relevant linear operator:

P₀A₀＝A₁

then, the following steps are obtained:

wherein the content of the first and second substances,

by P_a,P_bRespectively representing the front M-K line and the rear M-K line of P, using A_a,A_bIs represented by A₀Front M-K lines and rear M-K lines; comprises the following steps:

wherein, P_a,P_bThe relationship of Φ is:

wherein the content of the first and second substances,

Φ＝diag{ξ₁,ξ₂,...,ξ_Kthe kth diagonal element in Φ corresponds to the sensor space phase factor of the kth signal;

and Φ is a similarity matrix, pair

Is equal to the diagonal element xi_kK is 1,2,. K; pitch angle estimation result phi_kObtained from the formula:

compared with the prior art, the invention has the beneficial effects that: 1. modeling the vector sensor output signal in the quaternion domain preserves the vector characteristics of the signal and also makes full use of the orthogonal relationship between the components of the signal. Compared with a long vector modeling mode in a complex domain, the quaternion model provides a more compact and convenient processing mode for multi-component complex signals. 2. The spatial domain phase information between adjacent sensors in the acoustic vector sensor array and the inherent direction cosine information in the vibration velocity sensor can be fully utilized, and the gain of the algorithm performance is realized by more fully utilizing the information. 3. The method not only considers the cross-correlation information among the array elements, but also considers the self-correlation information. In the traditional method, additive noise introduced by autocorrelation cannot be removed, which seriously limits the performance of a propagation operator method, but in the method, the additive noise is statistically eliminated by processing data in a quaternion domain, so that the propagation operator method obtains a more accurate DOA estimation result.

Drawings

FIG. 1 is a block flow diagram of the present invention;

FIG. 2 is a schematic view of array layout;

FIG. 3 is a comparison of algorithm performance under different SNR (in the figure, AQ-SUMWE is the method proposed in this patent, V-SUMWE and FBSS-ESPRIT are the existing algorithms, the same below);

FIG. 4 is a comparison of algorithm performance for different fast beat numbers;

FIG. 5 is a comparison of algorithm performance for different array element numbers.

Detailed Description

The present invention will be described in further detail with reference to the accompanying drawings and specific embodiments.

With reference to fig. 1 to 5, the steps of the present invention are as follows:

the method comprises the following steps: establishing an observation signal complex model:

assuming that the array is a linear array of vector sensors arranged at equal intervals of half wavelength along the z-axis, and a single vector sensor comprises a sound pressure sensor and two particle vibration velocity sensors respectively pointing to the x-axis and the y-axis, the manifold vector of the vector sensor can be written as follows:

wherein, theta_kAnd phi_kHorizontal and pitch angles, u, of the kth signal, respectively_kAnd v_kDirection cosines, corresponding to particle velocity sensors pointing in the x-axis and y-axis, respectively. Assuming that K far-field narrow-band signals are incident on the sensor array, receiving signal models of the sound pressure sensor array, the x-axis-oriented particle velocity sensor array and the y-axis-oriented particle velocity sensor array can be respectively established in a complex number domain, and the following steps are performed in sequence:

wherein A (phi) is [ a (phi) ]₁),a(φ₂),...,a(φ_K)]Is a manifold of an array of acoustic pressure sensors,

for the steering vector corresponding to the k-th signal,

is space domain phase factor between sensors, d is array element spacing, lambda is signal wavelength, and signal vector s (t) is [ s ]₁(t),...,s_K(t)]^TIs a complex gaussian random process, beta ═ beta₁,...,β_K]^TIs a complex attenuation coefficient of which beta ₁1 and beta_k≠0,k＝2,3,…K，n_lP, x, y are noise vectors which are uncorrelated with the signal and independent of one another, Γ_x＝diag{cosθ₁sinφ₁,...,cosθ_K sinφ_K}，Γ_y＝diag{sinθ₁sinφ₁,...,sinθ_K sinφ_K}。

Step two: establishing an observation signal quaternion model:

now, the following two quaternion signal models are constructed according to the above signal model:

model one:

model two:

step three: constructing a correlation matrix:

taking into account the cross-correlation between the quaternion model of the m-th array element output signal and the quaternion model of the entire array output signal, there

Will r is_mDeployed, in the following form:

step four: the quaternion cross-correlation vector is simplified,

according to the quaternion operation rule: for the

Satisfy jx ═ x^*j, considering the signal as a circular signal, we have

E{s(t)js^H(t)}＝E{βs₁(t)s₁(t)β^Tj}＝0 (4.1)

Therefore, r_mCan be simplified into the following form:

wherein the content of the first and second substances,

when l ≠ m, there are

When l is m and the noise is a round signal, there are

In the above derivation, we used the following two assumptions: 1) the noise being spatially independent, i.e.

2) Noise vector n_p(t),n_x(t),n_y(t) are all round signals, i.e., E { n }_l,m(t)n_l,m(t) } 0, l ═ p, x, y. By this approach, we successfully eliminated r at a statistical level using quaternion algebra theory_mAdditive noise in (1).

According to the above derivation, r_mThe final form of (a) is:

wherein the content of the first and second substances,

R_s＝E{s(t)s^H(t) } is the signal covariance matrix.

Step five: source decoherence is achieved by correlation matrix smoothing:

considering K coherent signals, we pass through the correlation matrix r_mSmoothing to construct a matrix

Can be expressed as

Wherein A is₀The first M-K +1 line containing A,

is the selection of the matrix or matrices,

and

are two matrices of zero which are then used,

it can be found that,

is K, becauseThis is achieved by

Can be used to estimate the direction of the K coherent signals.

Finally, we define the following correlation matrix R for subsequent processing:

wherein pi ═ pi₁,Π₂,...,Π_M]。

Step six: estimating the signal pitch angle by using a propagation operator method:

propagation operator method based on pair matrix A₀Is divided into sub-arrays, i.e.

Wherein, because A₁Has a Vandermonde structure, and is therefore a full-rank (nonsingular) matrix of K rows, A₂Is a matrix of M-2K +1 rows. A. the₀Is linearly independent, then A₀The remaining rows of (a) may be represented as a₀Linear combinations of the first K rows. Propagation matrix P₀Can be uniquely defined as a and matrix A₁、A₂The relevant linear operator:

P₀A₀＝A₁ (6.2)

from, get

Wherein the content of the first and second substances,

by P_a,P_bRespectively represent the front M-K line and the rear M-K line of P, and similarly, use A_a,A_bIs represented by A₀The front M-K row and the back M-K row. We have found thatIs provided with

Wherein, P_a,P_bThe relationship of phi

Here, the first and second liquid crystal display panels are,

Φ＝diag{ξ₁,ξ₂,...,ξ_Kthe kth diagonal element in Φ corresponds to the sensor space phase factor of the kth signal. According to the formula, the compound can be obtained,

and Φ is a similarity matrix, pair

Is equal to the diagonal element xi_kK is 1, 2. Therefore, the pitch angle estimation result φ_kCan be obtained from the following formula

Step seven: estimating a propagation operator:

to estimate the propagation operator, we partition R as follows

Wherein R is₁，R₂Are respectively the front K line and the rear M-2K +1 line of R, and are respectively provided with

R₂＝P₀R₁ (7.2)

Thus, the propagation operator can be expressed as

Estimating propagation factors from sample data

Pitch angle

Can be selected from

K characteristic values of

Obtaining:

step eight: the azimuth of the incident signal level is estimated,

here we obtain an azimuth estimate by using directional information embedded in the acoustic particle velocity component. First, by concatenating x (t) and y (t) given in (a), the following array output vector z' (t) can be formed:

then, the following 2M × M cross-correlation matrix R' is introduced:

it is noted that,

and n_p(t) spatially from each otherThis is independent, i.e.

Thus, statistically, the noise in R' has been eliminated. We are with R'_(1:M,:)、R′_(M+1:2M,:)Respectively represent the front M rows and the rear M rows of R', will

Into, there is

w_x,k,w_y,kK is W, 1,2, K is W, respectively_x,W_yThe k-th row of (1). From the relationship: u. of_k tanθ_k＝v_kIt is known that

tanθ_kw_x,k＝w_y,k (8.4)

Step nine: the array response matrix is reconstructed and the array response matrix,

to estimate the array response matrix a, we can refer to a (Φ) ═ a (Φ) for₁),a(φ₂),…,a(φ_K)]And use

K1, 2, …, K to reconstruct the matrix a, i.e.

Thus, W_x,W_yIs estimated as

Then, from the sum, the azimuth of the k-th signal can be estimated as

Therefore, the pitch angle and the azimuth angle of the incident signal are obtained, and the pitch angle and the azimuth angle are automatically matched without additional operation. In fig. 3, fig. 4, and fig. 5, we respectively compare the performance of the proposed algorithm with that of the existing algorithms (V-sum, FBSS-ESPRIT) under different parameters (signal-to-noise ratio, fast beat number, array element number), and simulation results show that the proposed algorithms all obtain more accurate estimation results.

The specific implementation steps of the algorithm are shown in fig. 1, and it should be noted that the "method for removing additive noise by using quaternion" proposed in this patent has certain universality, and can be regarded as a preprocessing on data, so that the method can be combined with multiple algorithms (for example, sumpe in this patent) to further improve the performance of the algorithm. Without departing from the technical principle of the invention, a person skilled in the art can also make several modifications to the method (for example, combine the method with other existing algorithms), and these modifications should also fall within the protection scope of the invention.

Claims (4)

1. A quaternion-based coherent signal two-dimensional DOA estimation method is characterized by comprising the following steps: the method comprises the following steps:

the method comprises the following steps: establishing an observation signal complex model;

step two: establishing an observation signal quaternion model; two quaternion signal models were constructed as follows:

model one:

z(t)＝p(t)+x(t)j+y(t)j

＝A(s(t)+Γ_xs(t)j+Γ_ys(t)j)+n(t)

model two:

step three: constructing a correlation matrix: considering the cross-correlation between the quaternion model of the m-th array element output signal and the quaternion model of the entire array output signal, there are:

will r is_mDeployed, in the following form:

step four: quaternion cross-correlation vectorial reduction, r_mThe final form of (a) is:

wherein the content of the first and second substances,

R_s＝E{s(t)s^H(t) } is a signal covariance matrix;

step five: the information source decoherence is realized through the correlation matrix smoothing;

step six: estimating the signal pitch angle by using a propagation operator method:

step seven: estimating a propagation operator, and dividing R as follows:

wherein R is₁，R₂Are respectively the front K line and the back M-2K +1 line of R₂＝P₀R₁

Thus, the propagation operator can be expressed as:

estimating propagation factors from sample data

Pitch angle

Can be selected from

K characteristic values of

Obtaining:

step eight: the azimuth of the incident signal level is estimated,

concatenating x (t) and y (t) given in (a) to form the following array output vector z' (t):

the following 2M × M cross-correlation matrix R' is introduced:

wherein

And n_p(t) are spatially independent of each other, i.e.

Thus statistically, the noise in R' has been eliminated; from R'_(1:M,:)、R′_(M+1:2M,:)Respectively represent the front M rows and the rear M rows of R', will

Substituting, there are:

w_x,k,w_y,kk is W, 1,2, K is W, respectively_x,W_yThe kth line of (1); from the relationship: u. of_k tanθ_k＝v_kIt is known that

tanθ_kw_x,k＝w_y,k

Step nine: reconstructing an array response matrix, with reference to a (phi) ═ a (phi)₁),a(φ₂),...,a(φ_K)]And use

To reconstruct the matrix A, i.e.

W_x,W_yThe estimated values of (c) are:

the azimuth angle of the kth signal can be estimated as:

therefore, the pitch angle and the azimuth angle of the incident signal are obtained, and the pitch angle and the azimuth angle are automatically matched.

2. The quaternion-based coherent signal two-dimensional DOA estimation method of claim 1, wherein: the first step is specifically as follows: assuming that the array is a linear array of vector sensors arranged at equal intervals of half wavelength along the z-axis, and a single vector sensor comprises a sound pressure sensor and two particle vibration velocity sensors respectively pointing to the x-axis and the y-axis, the manifold vector of the vector sensor can be written as follows:

wherein, theta_kAnd phi_kHorizontal and pitch angles, u, of the kth signal, respectively_kAnd v_kThe direction cosine corresponds to a particle vibration velocity sensor pointing to an x axis and a y axis respectively; assuming that K far-field narrow-band signals are incident on the sensor array, receiving signal models of the sound pressure sensor array, the x-axis-oriented particle velocity sensor array and the y-axis-oriented particle velocity sensor array can be respectively established in a complex number domain, and the following steps are performed in sequence:

wherein A (phi) is [ a (phi) ]₁),a(φ₂),...,a(φ_K)]Is a manifold of an array of acoustic pressure sensors,

for the steering vector corresponding to the k-th signal,

is space domain phase factor between sensors, d is array element spacing, lambda is signal wavelength, and signal vector s (t) is [ s ]₁(t),…,s_K(t)]^TIs a complex gaussian random process, beta ═ beta₁,…,β_K]^TIs a complex attenuation coefficient of which beta₁1 and beta_k≠0,k＝2,3,…K，n_lP, x, y are noise vectors which are uncorrelated with the signal and independent of one another, Γ_x＝diag{cosθ₁sinφ₁,…,cosθ_Ksinφ_K}，Γ_y＝diag{sinθ₁sinφ₁,…,sinθ_Ksinφ_K}。

3. A quaternion-based coherent signal two-dimensional DOA estimation method according to claim 1 or 2 characterized by: the fifth step is specifically as follows: considering K coherent signals, we pass through the correlation matrix r_mSmoothing to construct a matrix

Can be expressed as:

wherein A is₀The first M-K +1 line containing A,

is the selection of the matrix or matrices,

and

is two zero matrices, Π_m＝[ψ_m,Φ¹ψ_m,…,Φ^K-1ψ_m],

It can be found that,

is K, so

Can be used to estimate the direction of the K coherent signals;

the following correlation matrix R is defined for subsequent processing:

wherein pi ═ pi₁,Π₂,...,Π_M]。

4. The quaternion-based coherent signal two-dimensional DOA estimation method of claim 3, wherein: the sixth step is specifically as follows: propagation operator method based on pair matrix A₀Is divided into sub-arrays, i.e.

Wherein, because A₁Is a full rank matrix of K rows, A₂Is a matrix of M-2K +1 rows; a. the₀Is linearly independent, then A₀The remaining rows of (a) may be represented as a₀Linear combination of the first K rows; propagation matrix P₀Is defined as an AND matrix A₁、A₂The relevant linear operator:

P₀A₀＝A₁

then, the following steps are obtained:

wherein the content of the first and second substances,

by P_a,P_bRespectively representing the front M-K line and the rear M-K line of P, using A_a,A_bIs represented by A₀Front M-K lines and rear M-K lines; comprises the following steps:

wherein, P_a,P_bThe relationship of Φ is:

wherein the content of the first and second substances,

Φ＝diag{ξ₁,ξ₂,...,ξ_Kthe kth diagonal element in Φ corresponds to the sensor space phase factor of the kth signal;

and Φ is a similarity matrix, pair

Is equal to the diagonal element xi_kK is 1,2,. K; pitch angle estimation result phi_kObtained from the formula:

Images