CN113283055A - Multi-signal separation and direction finding combined processing method based on parallel factor model - Google Patents

Abstract

The invention discloses a multi-signal separation and direction finding joint processing method based on a parallel factor model. A uniform area array is used to collect a mixed electromagnetic signal containing multiple signals as a received sampling signal; the received sampling signal is preprocessed, and a subspace is used Rotation algorithm performs angle estimation to obtain a preliminary angle estimation value; constructs an initial estimated steering matrix according to the obtained preliminary angle estimation value; re-models the received sampling signal as a PARAFAC model; uses the constructed steering matrix as a PARAFAC decomposition For the initial value, use the TALS algorithm to fit the PARAFAC model until the convergence conditions are met, and obtain the source matrix and steering matrix; extract the separation signal from the estimated source matrix; estimate the angle of arrival corresponding to the separation signal from the estimated steering matrix. The invention has low computational complexity, can effectively separate source signals and estimate corresponding angle parameters.

Description

A Combined Processing Method of Multiple Signal Separation and Direction Finding Based on Parallel Factor Model

technical field

The invention belongs to the technical field of array signal processing, and in particular relates to a multi-signal separation and direction finding joint processing method based on a parallel factor model applied to a uniform area array.

Background technique

Multi-signal separation and direction finding are key issues in radio spectrum monitoring, electronic reconnaissance and wireless communication. The main implementation method is to use an array antenna to receive electromagnetic signals propagated in space, and to estimate the arrival directions of different signals through digital signal processing. extracted from the mixed signal.

At present, beamforming and subspace estimation techniques are mainly used in multi-signal separation and direction finding, and widely used subspace rotation methods such as ESPRIT, PM and other algorithms, which use the signal subspace rotation invariance of the signal vector for parameter estimation , but these two algorithms can only estimate the angle, cannot separate the source signal, and the angle estimation performance is limited. The beamforming method can separate the signal in a certain direction, but this method has low angular resolution.

Tensor decomposition is an emerging technology in signal processing and data analysis, which has been widely used in biomedicine, wireless communication, and machine learning. Among them, parallel factorization (PARAFAC) decomposition is one of the most commonly used tensor decomposition methods. The parallel factor model is applied to the array signal processing, and the trilinear alternating least squares (TALS) can be used to realize source separation and direction finding, and can achieve better performance. But this standard PARAFAC decomposition method has the problems of slow convergence speed and high computational complexity.

SUMMARY OF THE INVENTION

Purpose of the invention: The present invention provides a multi-signal separation and direction finding joint processing method based on a parallel factor model, which can improve the performance of signal separation and angle estimation under uniform area arrays, and reduce the computational complexity of standard PARAFAC decomposition.

Technical solution: A method for multi-signal separation and direction finding joint processing based on a parallel factor model described in the present invention includes the following steps:

(1) Collecting the mixed electromagnetic signal containing multiple signals by the uniform area array as the received sampling signal;

(2) Preprocessing the received sampled signal, and using the subspace rotation algorithm to estimate the angle to obtain a preliminary angle estimate;

(3) Construct the steering matrix of the initial estimate;

(4) remodeling the received sampling signal in step (1) as a PARAFAC model;

(5) Take the steering matrix as the initial value of the PARAFAC decomposition, use the TALS algorithm to fit the PARAFAC model until the convergence conditions are met, and obtain the source matrix and the steering matrix;

(6) extracting the separation signal from the source matrix estimated in step (5);

(7) Estimate the angle of arrival corresponding to the separated signal from the steering matrix estimated in step (5).

Further, described step (2) comprises the following steps:

(21) Estimate the covariance matrix using the received sampled signal

的分块子矩阵计算传播算子

的最小二乘估计，构造信号子空间矩阵E；(22) Using the covariance matrix

The block submatrix computes the propagation operator

The least squares estimation of , constructs the signal subspace matrix E;

(23) Use the block submatrix of the signal subspace matrix E to calculate the rotation matrix

进行特征分解，求得其特征值，根据特征值计算出角频率

其中K为信号数目；(24) Pair rotation matrix

Perform eigendecomposition, obtain its eigenvalues, and calculate the angular frequency according to the eigenvalues

where K is the number of signals;

(25) Reconstruct the signal subspace matrix E to obtain a new signal subspace matrix E', repeat steps (23) (24), and calculate another angular frequency

和方位角

(26) Calculate the signal pitch angle according to the formula u _k = sinθ _k cosφ _k and v _k =sinθ _k sinφ _k

and azimuth

Further, the initial estimated steering matrices described in step (3) are respectively:

A _x =[a _x (u ₁ ),a _x (u ₂ ),...,a _x (u _K )]

A _y =[a _y (v ₁ ),a _y (v ₂ ),...,a _y (v _K )]

where a _x (u _k ) and a _y (v _k ) are the steering vectors on the x and y axes, respectively.

Further, described step (4) realization process is as follows:

表示,沿着三个不同的维度切分和拼接得到三个数据矩阵：

和

其中N和M分别表示均匀面阵沿x轴和y轴的阵元数目，L表示时域采样点数。The received sampled signal is reused as a third-order tensor according to the parallel factor model

Representation, splitting and splicing along three different dimensions to obtain three data matrices:

and

Among them, N and M represent the number of array elements along the x-axis and y-axis of the uniform surface array, respectively, and L represents the number of sampling points in the time domain.

Further, described step (5) comprises the following steps:

和

(51) Initialize with the initial steering matrix

and

的最小二乘估计，

其中，符号⊙表示Khatri–Rao积，上标

表示伪逆；(52) Calculation

The least squares estimate of ,

where the symbol ⊙ represents the Khatri–Rao product, and the superscript

represents a pseudo-inverse;

的最小二乘估计，

(53) Calculation

The least squares estimate of ,

的最小二乘估计，

(54) Calculation

The least squares estimate of ,

(55) Determine whether the set convergence condition is reached, and the algorithm stops when it is reached; otherwise, go back to step (52) to continue calculating a new estimated value.

Further, described step (7) comprises the following steps:

和

的列向量进行归一化，使首项等于1；(71) will

and

The column vector of is normalized so that the first term is equal to 1;

和

的第k列向量，计算r_x＝-angle(a_xk)，r_y＝-angle(a_yk)，angle()表示计算相位角；(72) represented by a _xk and a _yk

and

The k-th column vector of , calculate r _x =-angle(a _xk ), r _y =-angle(a _yk ), angle() means to calculate the phase angle;

(73) According to the relationship between phase and angular frequency, calculate the least squares estimate of angular frequency:

in,

B _x =[0,2πd _x /λ,...,2π(N-1)d _x /λ] ^T

B _y =[0,2πd _y /λ,...,2π(M-1)d _y /λ] ^T

Among them, d _x = _dy =d is the spacing of the array elements, and λ is the wavelength;

和方位角

(74) Calculate the signal pitch angle according to the formula u _k = sinθ _k cosφ _k and v _k =sinθ _k sinφ _k

and azimuth

Beneficial effects: Compared with the prior art, the beneficial effects of the present invention: 1. The multi-signal separation and direction finding joint processing method provided by the present invention applies the parallel factor analysis model to the field of array signal processing, and measures the The direction accuracy is better than that of traditional PM, ESPRIT and other algorithms, and no angle pairing is required, and separate signals can be obtained while measuring the angle;

2. The method of the present invention uses the PM algorithm to perform the initial estimation of the angle, so that the TALS algorithm can reach the convergence stagnation point faster, and the calculation amount of the standard TALS algorithm using random initialization is reduced.

Description of drawings

Fig. 1 is the flow chart of the present invention;

2 is a schematic diagram of a uniform surface array structure involved in the method of the present invention;

3 is a schematic diagram of the separated signal waveform obtained by the processing method of the present invention when the signal-to-noise ratio is 10dB;

4 is a scatter plot of angle estimation obtained by the processing method of the present invention when the signal-to-noise ratio is 10dB;

5 is a comparison diagram of the convergence speed of the method of the present invention and the standard PARAFAC algorithm under the same array structure and the same number of snapshots;

6 is a comparison diagram of the angle estimation performance of the method of the present invention and the four methods of PM, ESPRIT and standard PARAFAC under the same array structure and the same number of snapshots;

7 is a comparison diagram of the signal separation performance of the present invention and the standard PARAFAC two methods under different signal-to-noise ratio conditions;

Detailed ways

The present invention will be further described in detail below in conjunction with the accompanying drawings.

The present invention provides a multi-signal separation and direction finding joint processing method based on a parallel factor model, as shown in Figure 1, which specifically includes the following steps:

Step 1: Use a uniform area array to collect a mixed electromagnetic signal containing multiple signals as a received sampling signal.

FIG. 2 is a schematic diagram of a uniform area array structure involved in the method of the present invention. The surface array has a total of N×M array elements, which are uniformly distributed, and the spacing between adjacent array elements is d=d _{x =} dy , and d≤λ/2 (λ is the wavelength). Suppose that there are K incoherent far-field signals incident on this uniform surface array, and their arrival directions are (θ _k , φ _k ), k=1,2,…,K, where θ _k and φ _k represent the first The pitch and azimuth of the k signals. The sampled signal received by the array can be written as:

为含有噪声的数据矩阵；L表示时域采样点数；

为信源矩阵；N为高斯白噪声矩阵；

为阵列流形矩阵，可以表示：in,

is a data matrix containing noise; L represents the number of sampling points in the time domain;

is the source matrix; N is the Gaussian white noise matrix;

is an array manifold matrix, which can be expressed as:

表示Kronecker积，u_k＝sinθ_kcosφ_k，v_k＝sinθ_ksinφ_k；a_x(u_k)和a_y(v_k)分别为x轴上和y轴上的导向矢量，可以表示为：Among them, the symbol

Representing the Kronecker product, u _k = sinθ _k cosφ _k , v _k = sinθ _k sinφ _k ; a _x (u _k ) and a _y (v _k ) are the steering vectors on the x-axis and y-axis, respectively, and can be expressed as:

a _x (u _k )=[1,exp(-j2πd _x u _k /λ),...,exp(-j2π(N-1)d _x u _k /λ)] ^T (3)

a _y (v _k )=[1,exp(-j2πd _y v _k /λ),...,exp(-j2π(M-1)d _y v _k /λ)] ^T (4)

Formula (2) can be further expressed as:

A=[A _y ⊙A _x ] (5)

where ⊙ represents the Khatri–Rao product, A _x =[a _x (u ₁ ),a _x (u ₂ ),...,a _x (u _K )], A _y =[a _y (v ₁ ), a _y (v ₂ ),...,a _y (v _K )].

Step 2: Preprocess the received sampling signal obtained in Step 1, and use the subspace rotation algorithm to perform angle estimation to obtain a preliminary angle estimation value.

After getting the received data matrix, use the PM algorithm to estimate the initial angle, including the following steps:

估计协方差矩阵

(1) Using the received sampled signal matrix

Estimating the covariance matrix

根据

进行分块，其中

是

的前K个列向量组成的子矩阵，

是剩余列向量组成的子矩阵；(2) Convert the covariance matrix

according to

block, where

Yes

The submatrix consisting of the first K column vectors of ,

is a submatrix composed of the remaining column vectors;

计算传播算子的最小二乘解

(3) According to the formula

Compute the least squares solution of the propagation operator

其中I_K是K×K的单位矩阵，构造矩阵E_x＝E_1:N(M-1),1:K和E_y＝E_N+1:MN,1:K，E_x为矩阵E的第1到第N(M-1)行组成的子矩阵，E_y为矩阵E的第N+1到第MN行组成的子矩阵；(4) Define the matrix

Wherein I _K is the identity matrix of K×K, the construction matrix E _x =E _1:N(M-1),1:K and E _y =E _N+1:MN,1:K , E _x is the matrix E A sub-matrix consisting of rows 1 to N (M-1), E _y is a sub-matrix consisting of rows N+1 to MN of matrix E;

(5) Calculate the matrix according to the formula

进行特征值分解，求得其特征值

计算出角频率：

其中λ为信号波长，angle(·)表示计算相位角；(6) Right

Perform eigenvalue decomposition to find its eigenvalues

Calculate the angular frequency:

where λ is the signal wavelength, and angle( ) represents the calculated phase angle;

(7) Reconstruct the matrix E to obtain the matrix E', and construct the matrix E _x '=E' _{1:M(N-1), 1:K} and E _y '=E' _{M+1:MN,1: K} , E _x ' is a submatrix composed of rows 1 to M(N-1) of matrix E', and E _y ' is a submatrix composed of rows M+1 to MN of matrix E';

(8) Calculate the matrix according to the formula

进行特征值分解，求得其特征值

计算出另一角频率：

(9) pair

Perform eigenvalue decomposition to find its eigenvalues

Calculate the other corner frequency:

其中

是第k个信号俯仰角的估计，

是第k个信号方位角的估计。(10) Calculated according to the formula

in

is the estimation of the pitch angle of the kth signal,

is an estimate of the kth signal azimuth.

Step 3: Construct an initial estimated steering matrix according to the initial estimated angle value obtained in Step 2.

After obtaining the initial estimated value of the angle calculated by the PM algorithm, the steering vector can be constructed according to formula (3) and formula (4), and then the initial estimated steering matrix can be composed:

A _x =[a _x (u ₁ ),a _x (u ₂ ),...,a _x (u _K )]

A _y =[a _y (v ₁ ),a _y (v ₂ ),...,a _y (v _K )]

where a _x (u _k ) and a _y (v _k ) are the steering vectors on the x and y axes, respectively.

Step 4: Remodel the received sampled signal in Step 1 as a PARAFAC model.

According to the parallel factor (PARAFAC) model, the array received signal can be expressed in the form of a trilinear model:

的第(n,l,m)个元素。Among them, A _x (n, k), A _y (m, k), S (l, k) are the (n, k)th element of the x-axis direction matrix A _x and the y-axis direction matrix A _y . (m,k) elements and the (l,k)th element of the source matrix S, x _n,l,m is a third-order tensor

The (n,l,m)th element of .

沿着三个不同的维度切分和拼接得到三个数据矩阵

和

其中N和M分别表示步骤1中所使用的均匀面阵沿x轴和y轴的阵元数目，L表示时域采样点数。Will

Three data matrices are obtained by splitting and splicing along three different dimensions

and

Among them, N and M respectively represent the number of array elements along the x-axis and y-axis of the uniform surface array used in step 1, and L represents the number of sampling points in the time domain.

Step 5: Use the steering matrix constructed in Step 3 as the initial value of the PARAFAC decomposition, and use the TALS algorithm to fit the PARAFAC model until the convergence conditions are met.

和

(1) Initialize the initial steering matrix obtained after the initial estimation of PM

and

的最小二乘估计，

(2) Calculation

The least squares estimate of ,

的最小二乘估计，

(3) Calculation

The least squares estimate of ,

的最小二乘估计，

(4) Calculation

The least squares estimate of ,

(5) Calculate the convergence rate of the residual sum of squares, and the algorithm stops when it is less than a certain small value set; otherwise, go back to (2) to continue to calculate a new estimated value.

Step 6: Extract the separated signal from the source matrix estimated in Step 5.

按行提取，得到K个分离信号向量{s_k|1≤k≤K}。After the TALS algorithm is completed, the final source matrix is

Extract by row, and obtain K separate signal vectors {s _k |1≤k≤K}.

Step 7: Estimate the angle of arrival corresponding to the separated signal from the steering matrix estimated in Step 5.

和

的列向量进行归一化，使首项等于1；(1) The last obtained in step 5

and

The column vector of is normalized so that the first term is equal to 1;

和

的第k列向量，计算r_x＝-angle(a_xk)，r_y＝-angle(a_yk)，angle(·)表示计算相位角；(2) Represented by a _xk and a _yk

and

The kth column vector of , calculate r _x =-angle(a _xk ), r _y =-angle(a _yk ), angle( ) represents the calculated phase angle;

其中B_x＝[0,2πd_x/λ,...,2π(N-1)d_x/λ]^T，B_y＝[0,2πd_y/λ,...,2π(M-1)d_y/λ]^T，d_x＝d_y＝d为阵元间距；(3) Calculate the least squares estimate of the angular frequency according to the relationship between the phase and the angular frequency

where B _x =[0,2πd _x /λ,...,2π(N-1)d _x /λ] ^T , B _y =[0,2πd _y /λ,...,2π(M-1) dy /λ] ^T , d _{x =dy =} d is the array element spacing;

其中符号|·|表示计算复数的模值。(4) Calculate the final angle estimate according to the formula

where the symbol |·| denotes the calculation of the modulo value of a complex number.

Taking an 8×8 uniform area array as an example, it is assumed that there are three different types of typical modulation signals in the space, which are the single carrier frequency signal s ₁ (t)=cos(2π×5×10 ⁶ t), the linear frequency modulation Signal s ₂ (t)=cos(π×10 ¹² t ² +2π×2×10 ⁶ t), amplitude modulation signal s ₃ (t)=cos(2π×3×10 ⁵ t)sin(2π×5× 10 ⁶ t), the arrival angles are (θ ₁ , φ ₁ )=(10°, 15°), (θ ₂ , φ ₂ )=(20°, 25°) and (θ ₃ , φ ₃ )=( 30°, 35°), the sampling frequency is 100MHz.

Fig. 3 is a schematic diagram of the separated signal waveform obtained by the processing method of the present invention when the signal-to-noise ratio is 10dB, Fig. 4 is a scatter diagram of the angle estimation obtained by the processing method of the present invention when the signal-to-noise ratio is 10dB, the number of sampling points in the time domain L= 800. It can be seen from FIG. 3 and FIG. 4 that the method of the present invention can effectively separate the source signal and estimate the corresponding angle parameter.

5 is a comparison diagram of the computational complexity of the method of the present invention and the standard PARAFAC algorithm under the conditions of the same array structure and the same number of sampling points, the number of sampling points in the time domain is L=800, and the number of simulation statistics is 1000. It can be seen from FIG. 5 that the computational complexity of the method of the present invention is lower than that of the standard PARAFAC decomposition method, and the convergence speed is increased by more than ten times in this embodiment.

Fig. 6a and Fig. 6b are the comparison charts of the angle estimation performance of the method of the present invention and the four methods of PM, ESPRIT and standard PARAFAC under the same array structure and the same number of snapshots, the number of sampling points in the time domain is L=800, and the simulation statistics The number of times is 1000. RMSE represents the root mean square error of the angle. It can be seen from the figure that the angle estimation performance of the method of the present invention is better than that of the 2D-PM and 2D-ESPRIT algorithms, and is close to the estimation performance of the standard PARAFAC algorithm.

7 is a comparison diagram of the signal separation performance of the method of the present invention and the standard PARAFAC method under different signal-to-noise ratio conditions, the number of sampling points in the time domain is L=800, and the number of simulation statistics is 1000. ρ represents the average similarity coefficient between the separated signal and the original signal. As can be seen from the figure, the signal separation performance of the method of the present invention is close to the estimation performance of the standard PARAFAC algorithm, and with the improvement of the signal-to-noise ratio, the separation performance is also improved.

As mentioned above, the above embodiments are only used to illustrate the technical solutions of the present invention, but not to limit them. Although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand: The technical solutions described in the embodiments are modified, or some technical features thereof are equivalently replaced; and these modifications or replacements, without departing from the principles of the present invention, should be regarded as the protection scope of the present invention.

Claims (6)

1. a multi-signal separation and direction finding joint processing method based on a parallel factor model, is characterized in that, comprises the following steps: (1) Collecting the mixed electromagnetic signal containing multiple signals by the uniform area array as the received sampling signal; (2) Preprocessing the received sampled signal, and using the subspace rotation algorithm to estimate the angle to obtain a preliminary angle estimate; (3) Construct the steering matrix of the initial estimate; (4) remodeling the received sampling signal in step (1) as a PARAFAC model; (5) Take the steering matrix as the initial value of the PARAFAC decomposition, use the TALS algorithm to fit the PARAFAC model until the convergence conditions are met, and obtain the source matrix and the steering matrix; (6) extracting the separation signal from the source matrix estimated in step (5); (7) Estimate the angle of arrival corresponding to the separated signal from the steering matrix estimated in step (5). 2. the multi-signal separation and direction finding joint processing method based on parallel factor model according to claim 1, is characterized in that, described step (2) comprises the following steps: (21) Estimate the covariance matrix using the received sampled signal

(22) Using the covariance matrix

The block submatrix computes the propagation operator

The least squares estimation of , constructs the signal subspace matrix E; (23) Use the block submatrix of the signal subspace matrix E to calculate the rotation matrix

(24) Pair rotation matrix

Perform eigendecomposition, obtain its eigenvalues, and calculate the angular frequency according to the eigenvalues

k=1,2,...,K; where K is the number of signals; (25) Reconstruct the signal subspace matrix E to obtain a new signal subspace matrix E', repeat steps (23) (24), and calculate another angular frequency

k=1,2,...,K; (26) Calculate the signal pitch angle according to the formula u _k = sinθ _k cosφ _k and v _k =sinθ _k sinφ _k

and azimuth

3. the multi-signal separation and direction finding joint processing method based on the parallel factor model according to claim 1, is characterized in that, the described initial estimation steering matrix of step (3) is respectively: A _x =[a _x (u ₁ ),a _x (u ₂ ),...,a _x (u _K )] A _y =[a _y (v ₁ ),a _y (v ₂ ),...,a _y (v _K )] where a _x (u _k ) and a _y (v _k ) are the steering vectors on the x and y axes, respectively. 4. the multi-signal separation and direction finding joint processing method based on parallel factor model according to claim 1, is characterized in that, described step (4) realization process is as follows: The received sampled signal is reused as a third-order tensor according to the parallel factor model

Representation, splitting and splicing along three different dimensions to obtain three data matrices:

and

Among them, N and M represent the number of array elements along the x-axis and y-axis of the uniform surface array, respectively, and L represents the number of sampling points in the time domain. 5. multi-signal separation and direction finding joint processing method based on parallel factor model according to claim 1, is characterized in that, described step (5) comprises the following steps: (51) Initialize with the initial steering matrix

and

(52) Calculation

The least squares estimate of ,

where the symbol ⊙ represents the Khatri–Rao product, and the superscript

represents a pseudo-inverse; (53) Calculation

The least squares estimate of ,

(54) Calculation

The least squares estimate of ,

(55) Determine whether the set convergence condition is reached, and the algorithm stops when it is reached; otherwise, go back to step (52) to continue calculating a new estimated value. 6. the multi-signal separation and direction finding joint processing method based on parallel factor model according to claim 1, is characterized in that, described step (7) comprises the following steps: (71) will

and

The column vector of is normalized so that the first term is equal to 1; (72) represented by a _xk and a _yk

and

The k-th column vector of , calculate r _x =-angle(a _xk ), r _y =-angle(a _yk ), angle() means to calculate the phase angle; (73) According to the relationship between phase and angular frequency, calculate the least squares estimate of angular frequency:

in, B _x =[0,2πd _x /λ,...,2π(N-1)d _x /λ] ^T B _y =[0,2πd _y /λ,...,2π(M-1)d _y /λ] ^T Among them, d _x = _dy =d is the spacing of the array elements, and λ is the wavelength; (74) Calculate the signal pitch angle according to the formula u _k = sinθ _k cosφ _k and v _k =sinθ _k sinφ _k

and azimuth

Images