CN106019214A - DOA estimation method for broadband coherent signal source - Google Patents

Abstract

The invention belongs to the field of signal processing to solve the DOA estimation problem of a distributed type wideband coherent signal source, and can estimate DOA of a signal source in a distributed type point source model and the expansion angle generated by factors such as scattering and multi-path. The technical schemes are related to a DOA estimation method for broadband coherent signal source, which comprises the following steps: firstly, constructing a signal model of a sensor array; secondly, performing fractional Fourier Transformation (FrFT) to the obtained signals; thirdly, simplifying the signal model in the frequency domain; and fourthly, using the TLS-ESPRIT algorithm for DOA estimation; when the signal models of two sensor sub-arrays are simplified, seeking for their respective covariance matrices; using singular value decomposition to separate the signal subspace from the noise subspace, and finally obtaining the DOA estimates of all the signal sources. The invention can be used on various occasions for signal processing.

Description

DOA estimation method of broadband coherent signal source

Technical Field

The invention belongs to the field of signal processing, and relates to the accurate estimation of the direction of arrival of a plurality of signal sources by using a sensor array. And more particularly, to a broadband coherent signal source DOA estimation algorithm.

Background

The DOA estimation of a signal source is one of main contents of research in the field of array signal processing, and has wide application in sonar, radar, seismic detection and communication systems.

The DOA estimation has been studied for a long time, forming a series of classical DOA estimation algorithms. The conventional beamforming algorithm is one of the earliest DOA estimation algorithms, and its implementation is simple, but its resolution is low and its noise immunity is weak. Then, it is proposed to apply a subspace decomposition technique to DOA estimation, and the two algorithms, namely a multiple signal classification algorithm (MUSIC) and a rotation invariant technique-based signal parameter estimation algorithm (ESPRIT), are the most classical and commonly used, and both of them obtain a signal subspace and a noise subspace through covariance matrix decomposition, and construct a power spectrum function to estimate the arrival angle of a signal source, and the algorithm has high accuracy and strong environment adaptability. The ESPRIT algorithm is more convenient to use because the sensor information is not required to be known in advance, and a series of improved algorithms such as LS-ESPRIT, TLS-ESPRIT and the like are derived. The algorithms described above are generally suitable for processing only narrow-band signals.

In wideband array signal processing, spatial parameter estimation is usually performed based on a point source signal model, and the energy of a point source signal is always transmitted along the DOA direction of the signal. During signal transmission, the resulting angular spread is in many cases not negligible due to spatial scattering and multipath effects. Therefore, a method is needed to simultaneously and accurately estimate the angle of arrival and the spread angle in the distributed signal model.

Usually, a Bayesian method and a high-order statistic method can be adopted to estimate spatial parameters of distributed broadband signal sources, but the methods are not suitable for complex scenes with more signal sources and serious point source scattering and the like. In fact, the MUSIC algorithm can extend signal processing from narrow band to wide band by improvement, and can solve the influence caused by signal scattering, but the algorithm is not applicable in the case of time-varying position vectors, and in the actual environment, the position vectors are time-varying in most cases. According to the spatial characteristics of the time-varying position vector, a traditional distributed signal source parameter estimation algorithm (DSPE) can be adopted for processing, but the algorithm requires that the number of source signals cannot exceed the number of sensors, and signal sources are independent, so that the application of the method has great limitation.

Disclosure of Invention

In order to overcome the defects of the prior art, the invention aims to solve the DOA estimation problem of a distributed broadband coherent signal source, and can estimate the DOA of the signal source in a distributed point source model and the extension angle generated by scattering, multipath and other factors. The technical scheme adopted by the invention is that the DOA estimation method of the broadband coherent signal source comprises the following steps:

the method comprises the following steps: firstly, constructing a signal model of a sensor array, supposing that Q signal sources in a space are incident to a pair of sensor arrays which are uniformly and linearly arranged, wherein the number of array elements in each group is P, P is less than Q, and the spacing between the array elements is d₁,d₂Then, the model of the received signal at the ith sensor in each group of subarrays at time t is expressed as:

in the formula,which represents a position vector at the sensor or sensors,is the angular density of the ith signal source, where_iFor unknown position parameters including central angle and spread angle,mean value of 0, variance of o²White gaussian noise, which represents the signal x_iThe noise of (a) the noise of (b),is the angle of incidence;

for coherent signal sources, its angular densityBy usingPerforming substitution, wherein s (t) is a random variable reflecting the signal source characteristics,an angular signal density, which is a complex equation with respect to the angle θ, representing the spatial distribution of the signal;

step two: performing fractional Fourier transform (FrFT) on the obtained signal to obtain

In the formula,is a position vector matrix of the sensor array, S_i(α, u) is the ith signal source, N_x(α, u) is the noise at the sensor array;

according to the rotation characteristic of the fractional Fourier transform, the frequency modulation domain and the energy concentration domain are perpendicular to each other, and the rotation angles of the signals in the two domains are respectively expressed as

α_d＝arctanμ₀

α_e＝-arccotμ₀

In the formula, mu₀The frequency modulation of the signal is such that the two rotation angles satisfy α_d＝α_e+ π/2, thus obtaining

X(α_d,u)＝F[X(α_e,u)]

In the formula, F (-) represents fractional Fourier transform, the above formula represents that the fractional Fourier transform of the signal in the modulation frequency domain is equal to the fractional Fourier transform of the signal in the energy concentration domain, and in the actual measurement, the signal is extracted by a peak search method according to the energy value of the signal, namely the energy concentration domain signal X (α)_eU) converting it into the FM domain by means of fractional Fourier transform to obtain a signal X (α)_dU) to determine the rotation angle α of the signal source_d；

Step three: simplifying signal models in the FM domain

Define parameter b (α)_d,ψ_i(ii) a u) is

Thus, the sensor array signal model is represented in the frequency-modulated domain as

Further obtain

X(α_d,u)＝B(α_d,ψ；u)·S(α_d,u)+N_x(α_d,u)

In the formula,

B(α_d,ψ；u)＝[b(α_d,ψ₁；u),...,b(α_d,ψ_Q；u)]

S(α_d,u)＝[s₁(α_d,u),...,s_Q(α_d,u)]^T

and the other set of sub-arrays is expressed in the frequency modulation domain as

Y(α_d,u)＝C(α_d,ψ；u)S(α_d,u)+N_y(α_d,u)

In the formula,

C(α_d,ψ；u)＝[c(α_d,ψ₁；u),...,c(α_d,ψ_Q；u)]

for a uniform linear array, the position vector matrix of the sensor array is

And wherein the position vector at each sensor is represented in the frequency modulated domain as

Wherein,representing the time delay of the signal arriving at the kth sensor, and c is the speed of light;

due to the estimated angle of the ith signalFrom its actual angle theta_0iClose, i.e. areThus making use of TaylorThe series expands the above formula at theta

Is further shown as

Wherein

Due to the fact that in the formulaAbout an angle theta_0iNormalized symmetry parameter of, andthus M_n,iThe medium odd term is 0, resulting in:

$b(\mathrm{\α},{\mathrm{\ψ}}_i)=\underset{k=0}{\overset{\mathrm{\∞}}{\mathrm{\Σ}}}\frac{A^{\left(2k\right)}(\mathrm{\α},{\mathrm{\θ}}_{0i})}{\left(2k\right)!}{M}_{2k,i}$

same c (α)_d,ψ_i) After simplification become

$c({\mathrm{\α}}_d,{\mathrm{\ψ}}_i)\≈e^{jf\left(\mathrm{\θ}\right)}\underset{k=0}{\overset{\mathrm{\∞}}{\mathrm{\Σ}}}\frac{A^{\left(2k\right)}({\mathrm{\α}}_d,{\mathrm{\θ}}_{0i})e^{jf\left({\mathrm{\θ}}_{0i}\right)}}{\left(2k\right)!}{M}_{2k,i}$

Thus, c (α, psi) can be obtained_i) And b (α, psi)_i) Satisfy the following relationship

Step four: and carrying out DOA estimation by using a TLS-ESPRIT algorithm, respectively calculating respective covariance matrixes of the two sensor sub-arrays after simplifying signal models of the two sensor sub-arrays, separating a signal sub-space from a noise sub-space by using singular value decomposition, and finally obtaining DOA estimation values of all signal sources.

Processing the signal model in the step one, which can be divided into the following steps:

1) converting the obtained signal into a frequency modulation domain for processing by utilizing fractional Fourier transform, wherein the expression is

X(α_d,u)＝B(α_d,ψ；u)S(α_d,u)+N_x(α_d,u)

The fractional order correlation matrix of the signal is

R_x(α_d,u)＝E[X(α_d,u)2X(α_d,u)^H]S(α_d,u)+N_x(α_d,u)

＝B(α_d,ψ；u)E[S(α_d,u)·S(α_d,u)^H]B(α_d,ψ；u)^H+E[N_X(α_d,u)·N_X(α_d,u)^H]

+B(α_d,ψ；u)E[S(α_d,u)·N(α_d,u)^H]+E[N_X(α_d,u)·S_X(α_d,u)^H]B(α_d,ψ；u)^H

Since the signal is uncorrelated with noise and the array element noises are uncorrelated with each other, then

R_x(α_d,u)＝B(α_d,ψ；u)E[S(α_d,u)·S(α_d,u)^H]B(α_d,ψ；u)^H+E[N_X(α_d,u)·N_X(α_d,u)^H]

＝B(α_d,ψ；u)E_SB(α_d,ψ；u)^H+E_N

The signal subspace E can be obtained by decomposition according to the formula_SAnd noise subspace E_NRotating the obtained received signals in a frequency modulation domain, wherein each signal obtains different rotation angles and can obtain an energy spectrum of the signal in a (α, u) domain;

2) searching a peak value in the energy spectrum to obtain the number of all sources;

3) because source signals with different frequency modulation rates have different rotation angles, the source signals have different spatial representations in an energy concentration domain, but in the domain, only one signal energy value reaches the maximum, and other signals and noise energy are relatively dispersed;

4) according to X (α)_d,u)＝F[X(α_e,u)]The filtered signal is subjected to fractional Fourier transform and convertedInto the frequency-modulated domain, the fractional order correlation matrix becomesThe signal subspaces E can be obtained respectively by eigenvalue decomposition_SAnd noise subspace E_N。

Utilizing TLS-ESPRIT algorithm to carry out DOA estimation, and the method comprises the following specific steps:

an improved algorithm TLS-ESPRIT based on ESPRIT algorithm is adopted to estimate DOA and an expansion angle, and a matrix Z (alpha, u) is defined firstly

$Z(\mathrm{\α},u)=\left[\begin{array}{c}X(\mathrm{\α},u)\\ Y(\mathrm{\α},u)\end{array}\right]$

Signal subspace matrix E for obtaining covariance matrix thereof_s(α, u), which may be denoted as E according to the model of reception of the array signal_s(α,u)＝[E₁(α,u)E₂(α,u)]^TThus, it can be known that there is a column full rank matrix T, so that

E₁(α,u)＝A₁(α,u)T(α,u)

E₂(α,u)＝A₂(α,u)T(α,u)＝A₁(α,u)Φ(α,u)T(α,u)

In the formula,due to E₁,E₂Are all full rank, so there is only one non-singular matrix Ψ (α, u) ═ T^-1(α, u) Φ (α, u) T (α, u) is such that

E₂(α,u)＝E₁(α,u)Ψ(α,u)

In actual environment, due to the influence of noise and other environmental factorsAndto replace E₁(α, u) and E₂(α, u) processing Ψ (α, u) with the TLS estimate

${\widehat{\mathrm{\Ψ}}}_{TLS}(\mathrm{\α},u)=-V_{1,2}(\mathrm{\α},u)V_{2,2}(\mathrm{\α},u)$

In the formula, V_1,2(α,u)，V_2,2(α, u) is determined by feature decomposition

$\left[\begin{array}{c}{\hat{E}}_1(\mathrm{\α},u)\\ {\hat{E}}_2(\mathrm{\α},u)\end{array}\right]\[\begin{array}{cc}{\hat{E}}_1(\mathrm{\α},u)& {\hat{E}}_2(\mathrm{\α},u)\end{array}\]=\left[\begin{array}{cc}{V}_{1,1}(\mathrm{\α},u)& V_{1,2}(\mathrm{\α},u)\\ V_{2,1}(\mathrm{\α},u)& V_{2,2}(\mathrm{\α},u)\end{array}\right]L(\mathrm{\α},u)\left[\begin{array}{cc}{V}_{1,1}^H(\mathrm{\α},u)& V_{1,2}^H(\mathrm{\α},u)\\ V_{2,1}^H(\mathrm{\α},u)& V_{2,2}^H(\mathrm{\α},u)\end{array}\right]$

And then will contain thereinAndestimating to obtain an estimated position vector

$\begin{array}{c}{\hat{A}}_1(\mathrm{\α},u)={\hat{E}}_1(\mathrm{\α},u){\hat{T}}^{-1}(\mathrm{\α},u)={\hat{E}}_2(\mathrm{\α},u){\hat{T}}^{-1}(\mathrm{\α},u){\widehat{\mathrm{\Φ}}}^{-1}(\mathrm{\α},u)\\ =\frac{1}{2}\{{\hat{E}}_1(\mathrm{\α},u){\hat{T}}^{-1}(\mathrm{\α},u)+{\hat{E}}_2(\mathrm{\α},u){\hat{T}}^{-1}(\mathrm{\α},u){\widehat{\mathrm{\Φ}}}^{-1}(\mathrm{\α},u)\}\end{array}$

According toThus estimating the angular information contained therein; and then, combining the peak value information of the signal source in the energy concentration domain to solve the extension angle information contained in the corresponding DOA estimation value. The invention has the characteristics and beneficial effects that:

1. the method can accurately estimate DOA values of all signal sources under the condition that the number of the signal sources is unknown;

2. the invention simplifies the signal model and reduces the operation amount;

3. the method solves the problem that the traditional DOA estimation algorithm cannot carry out accurate DOA estimation when the number of the sensors is less than the number of the sources;

4. the invention utilizes the improved TLS-ESPRIT algorithm to carry out DOA estimation of the signal source, thereby improving the estimation precision.

Description of the drawings:

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

FIG. 2 is a schematic diagram of a sensor array distribution;

FIG. 3 is a schematic diagram of the transformation of different signals into the frequency modulated domain;

FIG. 4 is a spectral diagram of a signal transformed to a modulation domain using the algorithm of the present invention;

fig. 5 is a schematic diagram of signal separation using a band pass filter.

Detailed Description

The realization process of the invention is as follows:

the method comprises the following steps: firstly, a signal model of a sensor array is constructed, and if Q signal sources in a space are incident to a pair of sensor arrays which are uniformly and linearly arranged, the number of array elements in each group is P, a received signal model at the ith sensor in each group of subarrays at the time t can be represented as

In the formula,which represents a position vector at the sensor or sensors,is the angular density of the ith source, where ψ is an unknown positional parameter containing the center angle and spread angle,mean value of 0, variance of o²White gaussian noise, which represents the signal x_iThe noise of (a) the noise of (b),is the angle of incidence.

For coherent signal sources, its angular densityCan usePerforming substitution, wherein s (t) is a random variable reflecting the signal source characteristics,to represent the angular signal density of the signal spatial distribution, it is a complex equation with respect to the angle θ.

Step two: performing fractional Fourier transform (FrFT) on the obtained signal

According to the rotation characteristic of the fractional Fourier transform, a frequency modulation domain (demodulation domain) and an Energy-concentrated domain (Energy-concentrated domain) are perpendicular to each other. The rotation angles of the signals in the two domains satisfy the following relation

α_d＝arctanμ₀

α_e＝-arccotμ₀

In the formula, mu₀The two rotation angles are α_d＝α_e+ π/2. Thus can obtain

X(α_d,u)＝F[X(α_e,u)]Wherein F (-) represents a fractional Fourier transformIn practical measurement, we can extract the signal by a peak search method according to the energy value of the signal, namely the signal X in the energy concentration domain (α)_eU) converting it into the FM domain by means of fractional Fourier transform to obtain a signal X (α)_dU) to determine the rotation angle α of the signal source_d。

Step three: the signal model in the frequency-modulated domain is simplified.

Defining the parameter b (α, psi)_i(ii) a u) is

Thus, the sensor array signal model may be represented in the frequency-modulated domain as

$X({\mathrm{\α}}_d,u)=\underset{i=1}{\overset{Q}{\mathrm{\Σ}}}b({\mathrm{\α}}_d,{\mathrm{\ψ}}_i;u)S_i({\mathrm{\α}}_d,u)+N_x({\mathrm{\α}}_d,u)$

And then can obtain

X(α_d,u)＝B(α_d,ψ；u)·S(α_d,u)+N_x(α_d,u)

In the formula,

B(α_d,ψ；u)＝[b(α_d,ψ₁；u),...,b(α_d,ψ_Q；u)]

S(α_d,u)＝[s₁(α_d,u),...,s_Q(α_d,u)]^T

and the other set of sub-arrays is expressed in the frequency modulation domain as

Y(α_d,u)＝C(α_d,ψ；u)·S(α_d,u)+N_y(α_d,u)

In the formula,

C(α_d,ψ；u)＝[c(α_d,ψ₁；u),...,c(α_d,ψ_Q；u)]

for a uniform linear array, the position vector matrix of the sensor array is

And wherein the position vector at each sensor is represented in the frequency modulated domain as

Wherein,representing the time delay of arrival of the signal at the kth sensor,is the angle of incidence, c is the speed of light.

Due to the estimated angle of the ith signalFrom its actual angle theta_0iClose, i.e. areThe Taylor series can therefore be used to expand the above equation at θ

Is further shown as

$b(\mathrm{\α},{\mathrm{\ψ}}_i)=\underset{n=0}{\overset{\mathrm{\∞}}{\mathrm{\Σ}}}\frac{A^{\left(n\right)}(\mathrm{\α},{\mathrm{\θ}}_{0i})}{n!}{M}_{n,i}$

Wherein

Due to the fact that in the formulaAbout an angle theta_0iNormalized symmetry parameter of, andthus M_n,iThe odd term is 0, can be obtained

$b(\mathrm{\α},{\mathrm{\ψ}}_i)=\underset{k=0}{\overset{\mathrm{\∞}}{\mathrm{\Σ}}}\frac{A^{\left(2k\right)}(\mathrm{\α},{\mathrm{\θ}}_{0i})}{\left(2k\right)!}{M}_{2k,i}$

Same c (α)_d,ψ_i) After simplification become

$c({\mathrm{\α}}_d,{\mathrm{\ψ}}_i)\≈e^{jf\left(\mathrm{\θ}\right)}\underset{k=0}{\overset{\mathrm{\∞}}{\mathrm{\Σ}}}\frac{A^{\left(2k\right)}({\mathrm{\α}}_d,{\mathrm{\θ}}_{0i})e^{jf\left({\mathrm{\θ}}_{0i}\right)}}{\left(2k\right)!}{M}_{2k,i}$

Thus, c (α, psi) can be obtained_i) And b (α, psi)_i) Satisfy the following relationship

And simplifying the signal model, then calculating a correlation function of the signal model, decomposing the signal model by using a singular value to obtain a signal subspace and a noise subspace, searching a peak value and filtering the signal, converting the filtered signal into a modulation frequency domain, and then decomposing the subspace.

Step four: DOA estimation is performed using the TLS-ESPRIT algorithm. And according to the signal subspace obtained after filtering in the third step, estimating the DOA and the expansion angle of the signal source by utilizing a TLS-ESPRIT algorithm.

The invention provides a broadband coherent signal source DOA estimation algorithm based on a TLS-ESPRIT algorithm, and a specific implementation process can be divided into three parts. The present invention will now be described more fully hereinafter with reference to the accompanying drawings.

The first part is the construction of a signal model.

Assuming that Q coherent signal sources are incident to a pair of sensor arrays which are uniformly and linearly arranged in space, the number of each group of array elements is P, the array distribution of the sensors is as shown in fig. 2, and then a received signal model at the ith sensor in each group of subarrays at time t can be represented as

The traditional DOA estimation algorithm is not applicable under the condition that the number of sensor array elements is less than that of signal sources, and in order to solve the problem, P is set to be less than Q in parameter design.

The second section mainly describes the signal model processing, which can be divided into the following steps.

1. Converting the obtained signal into a frequency modulation domain by utilizing fractional Fourier transform, and simplifying a signal model by utilizing the method in the third step to obtain an expression

X(α_d,u)＝B(α_d,ψ；u)S(α_d,u)+N_x(α_d,u)

The fractional order correlation matrix of the signal is

R_x(α_d,u)＝E[X(α_d,u)·X(α_d,u)^H]S(α_d,u)+N_x(α_d,u)

＝B(α_d,ψ；u)E[S(α_d,u)·S(α_d,u)^H]B(α_d,ψ；u)^H+E[N_X(α_d,u)·N_X(α_d,u)^H]

+B(α_d,ψ；u)E[S(α_d,u)·N(α_d,u)H]+E[N_X(α_d,u)·S_X(α_d,u)^H]B(α_d,ψ；u)^H

Since the signal is uncorrelated with noise and the array element noises are uncorrelated with each other, then

R_x(α_d,u)＝B(α_d,ψ；u)E[S(α_d,u)·S(α_d,u)^H]B(α_d,ψ；u)^H+E[N_X(α_d,u)·N_X(α_d,u)^H]

＝B(α_d,ψ；u)E_SB(α_d,ψ；u)^H+E_N

The signals obtained after the above simplified processing are rotated in the frequency modulation domain, and each signal obtains a different rotation angle, as shown in fig. 3, and simultaneously, the energy spectrum in the (α, u) domain can be obtained.

2. As shown in fig. 4, a peak search is performed in the energy spectrum, and the number of all sources can be obtained.

3. Since source signals with different frequencies have different rotation angles, they have different spatial representations in the energy concentration domain. But in this domain only one signal energy value is maximized, while the other signal and noise energy are relatively dispersed. Based on this property, it can be separated by a band-pass filter. As shown in fig. 5, the bandpass filter uses a rectangular window function whose bandwidth is the maximum of the distance between adjacent signal peaks, which can be used to eliminate the influence of other signal sources.

4. According to X (α)_d,u)＝F[X(α_e,u)]The fractional Fourier transform is carried out on the filtered signal and the signal is converted into a frequency modulation domain, and then the fractional correlation matrix becomesThe signal subspaces E can be obtained respectively by eigenvalue decomposition_SAnd noise subspace E_N。

The third part is DOA estimation using TLS-ESPRIT algorithm.

The conventional distributed signal source parameter estimation method (DSPE) performs parameter estimation based on noise feature vector minimization, i.e. for the spatial distribution parameter ψ to be estimated, we can find its estimated value by the following formula

$\widehat{\mathrm{\ψ}}=\arg \underset{\mathrm{\ψ}}{max}\frac{1}{\left|\right|b^H\left(\mathrm{\ψ}\right)E_n\left|\right|}$

In the formula,

although the algorithm can estimate the signal space parameters, the resolution is low, and the number of the sources has a large influence on the estimation accuracy, so that the improved algorithm TLS-ESPRIT based on the ESPRIT algorithm is adopted in the invention to estimate the DOA and the extension angle.

Here, a matrix Z (α, u) is first defined

$Z(\mathrm{\α},u)=\left[\begin{array}{c}X(\mathrm{\α},u)\\ Y(\mathrm{\α},u)\end{array}\right]$

From the second part, a signal subspace matrix E is obtained whose covariance matrix is_s(α, u), which may be represented as E_s(α,u)＝[E₁(α,u)E₂(α,u)]^T。

It can thus be seen that there is a column full rank matrix T, such that

E₁(α,u)＝A₁(α,u)T(α,u)

E₂(α,u)＝A₂(α,u)T(α,u)＝A₁(α,u)Φ(α,u)T(α,u)

In the formula,due to E₁,E₂Are all full rank, so there is only one non-singular matrix Ψ (α, u) ═ T^-1(α, u) Φ (α, u) T (α, u) is such that

E₂(α,u)＝E₁(α,u)Ψ(α,u)

In a real environment, usingAndto replace E₁(α, u) and E₂(α,u)。

Processing Ψ (α, u) using TLS estimates

${\widehat{\mathrm{\Ψ}}}_{TLS}(\mathrm{\α},u)=-V_{1,2}(\mathrm{\α},u)V_{2,2}(\mathrm{\α},u)$

In the formula, V_1,2(α,u)，V_2,2(α, u) can be found by feature decomposition

$\left[\begin{array}{c}{\hat{E}}_1(\mathrm{\α},u)\\ {\hat{E}}_2(\mathrm{\α},u)\end{array}\right]\[\begin{array}{cc}{\hat{E}}_1(\mathrm{\α},u)& {\hat{E}}_2(\mathrm{\α},u)\end{array}\]=\left[\begin{array}{cc}{V}_{1,1}(\mathrm{\α},u)& V_{1,2}(\mathrm{\α},u)\\ V_{2,1}(\mathrm{\α},u)& V_{2,2}(\mathrm{\α},u)\end{array}\right]L(\mathrm{\α},u)\left[\begin{array}{cc}{V}_{1,1}^H(\mathrm{\α},u)& V_{1,2}^H(\mathrm{\α},u)\\ V_{2,1}^H(\mathrm{\α},u)& V_{2,2}^H(\mathrm{\α},u)\end{array}\right]$

And can further include thereinAndand (4) estimating.

An estimated position vector may be derived

$\begin{array}{c}{\hat{A}}_1(\mathrm{\α},u)={\hat{E}}_1(\mathrm{\α},u){\hat{T}}^{-1}(\mathrm{\α},u)={\hat{E}}_2(\mathrm{\α},u){\hat{T}}^{-1}(\mathrm{\α},u){\widehat{\mathrm{\Φ}}}^{-1}(\mathrm{\α},u)\\ =\frac{1}{2}\{{\hat{E}}_1(\mathrm{\α},u){\hat{T}}^{-1}(\mathrm{\α},u)+{\hat{E}}_2(\mathrm{\α},u){\hat{T}}^{-1}(\mathrm{\α},u){\widehat{\mathrm{\Φ}}}^{-1}(\mathrm{\α},u)\}\end{array}$

According to the aboveThe angle information contained therein can be estimated. And then, combining the peak value information of the signal source in the energy concentration domain to solve the extension angle information contained in the corresponding DOA estimation value.

Claims (3)

1. A DOA estimation method of a broadband coherent signal source is characterized by comprising the following steps:

the method comprises the following steps: firstly, constructing a signal model of a sensor array, supposing that Q signal sources in a space are incident to a pair of sensor arrays which are uniformly and linearly arranged, wherein the number of array elements in each group is P, P is less than Q, and the spacing between the array elements is d₁,d₂Then, the model of the received signal at the ith sensor in each group of subarrays at time t is expressed as:

in the formula,which represents a position vector at the sensor or sensors,is the angular density of the ith signal source, where_iFor unknown position parameters including central angle and spread angle,is a mean of 0 and a variance of o²White gaussian noise, which represents the signal x_iThe noise of (a) the noise of (b),is the angle of incidence;

for coherent signal sources, its angular densityBy usingPerforming substitution, wherein s (t) is a random variable reflecting the signal source characteristics,an angular signal density, which is a complex equation with respect to the angle θ, representing the spatial distribution of the signal;

step two: performing fractional Fourier transform (FrFT) on the obtained signal to obtain

In the formula,is a position vector matrix of the sensor array, S_i(α, u) is the ith signal source, N_x(α, u) is the noise at the sensor array;

according to the rotation characteristic of the fractional Fourier transform, the frequency modulation domain and the energy concentration domain are perpendicular to each other, and the rotation angles of the signals in the two domains are respectively expressed as

α_d＝arctanμ₀

α_e＝-arccotμ₀

In the formula, mu₀The frequency modulation of the signal is such that the two rotation angles satisfy α_d＝α_e+ π/2, thus obtaining

X(α_d,u)＝F[X(α_e,u)]

In the formula, F (-) represents fractional Fourier transform, the above formula represents that the fractional Fourier transform of the signal in the modulation frequency domain is equal to the fractional Fourier transform of the signal in the energy concentration domain, and in the actual measurement, the signal is extracted by a peak search method according to the energy value of the signal, namely the energy concentration domain signal X (α)_eU) converting it into the FM domain by means of fractional Fourier transform to obtain a signal X (α)_dU) to determine the rotation angle α of the signal source_d；

Step three: simplifying signal models in the FM domain

Define parameter b (α)_d，ψ_i(ii) a u) is

Thus, the sensor array signal model is represented in the frequency-modulated domain as

$X({\mathrm{\α}}_d,u)=\underset{i=1}{\overset{Q}{\Σ}}b({\mathrm{\α}}_d,{\mathrm{\ψ}}_i;u)S_i({\mathrm{\α}}_d,u)+N_x({\mathrm{\α}}_d,u)$

Further obtain

X(α_d,u)＝B(α_d,ψ；u)·S(α_d,u)+N_x(α_d,u)

In the formula,

B(α_d,ψ；u)＝[b(α_d,ψ₁；u),...,b(α_d,ψ_Q；u)]

S(α_d,u)＝[s₁(α_d,u),...,s_Q(α_d,u)]^T

and the other set of sub-arrays is expressed in the frequency modulation domain as

Y(α_d,u)＝C(α_d,ψ；u)·S(α_d,u)+N_y(α_d,u)

In the formula,

C(α_d,ψ；u)＝[c(α_d,ψ₁；u),...,c(α_d,ψ_Q；u)]

for a uniform linear array, the position vector matrix of the sensor array is

And wherein the position vector at each sensor is represented in the frequency modulated domain as

Wherein,representing the time delay of the signal arriving at the kth sensor, and c is the speed of light;

due to the estimated angle of the ith signalFrom its actual angle theta_0iClose, i.e. areTherefore, the Taylor series is used to expand the above formula at theta

Is further shown as

$b(\mathrm{\α},{\mathrm{\ψ}}_i)=\underset{n=0}{\overset{\mathrm{\∞}}{\Σ}}\frac{A^{\left(n\right)}(\mathrm{\α},{\mathrm{\θ}}_{0i})}{n!}{M}_{n,i}$

Wherein

Due to the fact that in the formulaAbout an angle theta_0iNormalized symmetry parameter of, andthus M_n,iThe medium odd term is 0, resulting in:

$b(\mathrm{\α},{\mathrm{\ψ}}_i)=\underset{k=0}{\overset{\mathrm{\∞}}{\Σ}}\frac{A^{\left(2k\right)}(\mathrm{\α},{\mathrm{\θ}}_{0i})}{\left(2k\right)!}{M}_{2k,i}$

same c (α)_d，ψ_i) After simplification become

$c({\mathrm{\α}}_d,{\mathrm{\ψ}}_i)\≈e^{jf\left(\mathrm{\θ}\right)}\underset{k=0}{\overset{\mathrm{\∞}}{\Σ}}\frac{A^{\left(2k\right)}({\mathrm{\α}}_d,{\mathrm{\θ}}_{0i})e^{jf\left({\mathrm{\θ}}_{0i}\right)}}{\left(2k\right)!}{M}_{2k,i}$

Thus, c (α, psi) can be obtained_i) And b (α, psi)_i) Satisfy the following relationship

Step four: and carrying out DOA estimation by using a TLS-ESPRIT algorithm, respectively calculating respective covariance matrixes of the two sensor sub-arrays after simplifying signal models of the two sensor sub-arrays, separating a signal sub-space from a noise sub-space by using singular value decomposition, and finally obtaining DOA estimation values of all signal sources.

2. A DOA estimation algorithm for a wide band coherent signal source as claimed in claim 1, characterized in that the signal model in step one is processed, which is subdivided into the following steps:

1) converting the obtained signal into a frequency modulation domain for processing by utilizing fractional Fourier transform, wherein the expression is

X(α_d,u)＝B(α_d,ψ；u)·S(α_d,u)+N_x(α_d,u)

The fractional order correlation matrix of the signal is

R_x(α_d,u)＝E[X(α_d,u)·X(α_d,u)^H]S(α_d,u)+N_x(α_d,u)

＝B(α_d,ψ；u)E[S(α_d,u)·S(α_d,u)^H]B(α_d,ψ；u)^H+E[N_X(α_d,u)·N_X(α_d,u)^H]

+B(α_d,ψ；u)E[S(α_d,u)·N(α_d,u)^H]+E[N_X(α_d,u)·S_X(α_d,u)^H]B(α_d,ψ；u)^H

Since the signal is uncorrelated with noise and the array element noises are uncorrelated with each other, then

R_x(α_d,u)＝B(α_d,ψ；u)E[S(α_d,u)·S(α_d,u)^H]B(α_d,ψ；u)^H+E[N_X(α_d,u)·N_X(α_d,u)^H]

＝B(α_d,ψ；u)E_SB(α_d,ψ；u)^H+E_N

The signal subspace E can be obtained by decomposition according to the formula_SAnd noise subspace E_NRotating the obtained received signals in a frequency modulation domain, wherein each signal obtains different rotation angles and can obtain an energy spectrum of the signal in a (α, u) domain;

2) searching a peak value in the energy spectrum to obtain the number of all sources;

3) because source signals with different frequency modulation rates have different rotation angles, the source signals have different spatial representations in an energy concentration domain, but in the domain, only one signal energy value reaches the maximum, and other signals and noise energy are relatively dispersed;

4) according to X (α)_d,u)＝F[X(α_e,u)]The fractional Fourier transform is carried out on the filtered signal and the signal is converted into a frequency modulation domain, and then the fractional correlation matrix becomesThe signal subspaces E can be obtained respectively by eigenvalue decomposition_SAnd noise subspace E_N。

3. The DOA estimation method of a broadband coherent signal source according to claim 1, wherein the DOA estimation is performed by using TLS-ESPRIT algorithm, and the specific steps are as follows:

an improved algorithm TLS-ESPRIT based on ESPRIT algorithm is adopted to estimate DOA and an expansion angle, and a matrix Z (alpha, u) is defined firstly

$Z(\mathrm{\α},u)=\left[\begin{array}{c}X(\mathrm{\α},u)\\ Y(\mathrm{\α},u)\end{array}\right]$

Signal subspace matrix E for obtaining covariance matrix thereof_s(α, u), which may be denoted as E according to the model of reception of the array signal_s(α,u)＝[E₁(α,u)E₂(α,u)]^TThus, it can be known that there is a column full rank matrix T, so that

E₁(α,u)＝A₁(α,u)T(α,u)

E₂(α,u)＝A₂(α,u)T(α,u)＝A₁(α,u)Φ(α,u)T(α,u)

In the formula,due to E₁，E₂Are all full rank, so there is only one non-singular matrix Ψ (α, u) ═ T^-1(α, u) Φ (α, u) T (α, u) is such that

E₂(α,u)＝E₁(α,u)Ψ(α,u)

In actual environment, due to the influence of noise and other environmental factorsAndto replace E₁(α, u) and E₂(α, u) processing Ψ (α, u) with the TLS estimate

${\widehat{\mathrm{\Ψ}}}_{TLS}(\mathrm{\α},u)=-V_{1,2}(\mathrm{\α},u)V_{2,2}(\mathrm{\α},u)$

In the formula, V_1,2(α,u)，V_2,2(α, u) is determined by feature decomposition

$\left[\begin{array}{c}{\hat{E}}_1(\mathrm{\α},u)\\ {\hat{E}}_2(\mathrm{\α},u)\end{array}\right]\left[\begin{array}{cc}{\hat{E}}_1(\mathrm{\α},u)& {\hat{E}}_2(\mathrm{\α},u)\end{array}\right]=\left[\begin{array}{cc}{V}_{1,1}(\mathrm{\α},u)& V_{1,2}(\mathrm{\α},u)\\ V_{2,1}(\mathrm{\α},u)& V_{2,2}(\mathrm{\α},u)\end{array}\right]L(\mathrm{\α},u)\left[\begin{array}{cc}{V}_{1,1}^H(\mathrm{\α},u)& V_{1,2}^H(\mathrm{\α},u)\\ V_{2,1}^H(\mathrm{\α},u)& V_{2,2}^H(\mathrm{\α},u)\end{array}\right]$

And then will contain thereinAndestimating to obtain an estimated position vector

$\begin{array}{c}{\hat{A}}_1(\mathrm{\α},u)={\hat{E}}_1(\mathrm{\α},u){\hat{T}}^{-1}(\mathrm{\α},u)={\hat{E}}_2(\mathrm{\α},u){\hat{T}}^{-1}(\mathrm{\α},u){\widehat{\mathrm{\Φ}}}^{-1}(\mathrm{\α},u)\\ =\frac{1}{2}\{{\hat{E}}_1(\mathrm{\α},u){\hat{T}}^{-1}(\mathrm{\α},u)+{\hat{E}}_2(\mathrm{\α},u){\hat{T}}^{-1}(\mathrm{\α},u){\widehat{\mathrm{\Φ}}}^{-1}(\mathrm{\α},u)\}\end{array}$

According toThus estimating the angular information contained therein; and then, combining the peak value information of the signal source in the energy concentration domain to solve the extension angle information contained in the corresponding DOA estimation value.