阅读说明：本技术 一种基于空间时频分布的水下多径信号参数估计方法 (Underwater multipath signal parameter estimation method based on space time-frequency distribution ) 是由 刘妹琴 韩学艳 张森林 郑荣濠 吴争光 于 2020-06-11 设计创作，主要内容包括：本发明公开了一种基于空间时频分布的水下多径信号参数估计方法,该方法将传统的阵列信号处理技术与时频分布相结合以提高系统的信号处理能力。首先,采用二次时频分布来减少交叉项和提高时频分辨率。其次,通过构造扩展的二维空间时频分布矩阵实现在被动场景下未知多径信号数目的精确估计。再次,根据信号在时频域的能量分布,引入向前向后平滑的思想以保证脊线的局部平滑性,利用属于同一多径信号的自源时频点的主特征向量相同的特性,获得每条多径信号对应脊线的时频点集合。最后,根据脊线检测结果可在欠定条件下实现多径信号瞬时频率和方位的准确估计。(The invention discloses an underwater multipath signal parameter estimation method based on space time-frequency distribution, which combines the traditional array signal processing technology with the time-frequency distribution to improve the signal processing capability of a system. First, quadratic time-frequency distribution is used to reduce cross terms and improve time-frequency resolution. Secondly, the number of unknown multipath signals is accurately estimated in a passive scene by constructing an expanded two-dimensional space time-frequency distribution matrix. And thirdly, according to the energy distribution of the signals in the time-frequency domain, introducing the idea of forward and backward smoothing to ensure the local smoothness of the ridge line, and obtaining the time-frequency point set of the ridge line corresponding to each multi-path signal by utilizing the characteristic that the main characteristic vectors of the self-source time-frequency points belonging to the same multi-path signal are the same. And finally, according to the ridge line detection result, the accurate estimation of the instantaneous frequency and the direction of the multipath signal can be realized under an underdetermined condition.)

1. An underwater multipath signal parameter estimation method based on space time-frequency distribution is characterized by comprising the following steps:

step (1): at the time t, an underwater multipath signal is received by using a uniform linear array consisting of M omnidirectional array elements, and an array output vector Y (t) at the time t is obtained₁(t) y₂(t) … y_M(t)]^TWherein, y_m(t) represents the output of the m-th array element;

step (2): solving a spatial time-frequency distribution matrix D of an array output vector Y (t) by adopting smooth pseudo-wigner distribution_YY(t,f)；

And (3): for D_YYEach sampling time in (t, f)Setting a noise threshold value for the inter-slice, filtering all noise time frequency points on a time frequency plane, and storing the time frequency points selected by the noise threshold value to a set omega;

and (4): for D_YY(t, f) carrying out pre-whitening treatment to obtain a whitening space time-frequency distribution matrix

And (5): will set omega_AThe spatial time-frequency distribution matrixes corresponding to all self-source time-frequency points in the three-order tensor form a third-order tensorL represents the set omega_AThe number of internal self-sourced time frequency points; then converting the third order tensor P into an expanded two-dimensional space time-frequency distribution matrix

and (6): ridge line detection: extracting time-frequency ridges corresponding to each multipath signal by adopting forward search and backward search to obtain self-source time-frequency points, the number of the self-source time-frequency points and the corresponding instantaneous frequency of the self-source time-frequency points contained in each ridge;

and (7): implementing multipath signal instantaneous frequency according to ridge line detection result

2. The method for estimating the underwater multipath signal parameters based on the space time-frequency distribution according to claim 1, wherein the step (1) is specifically as follows:

at time t, the output y of the m-th array element_m(t) is expressed as:

wherein, a_k,nRepresenting the attenuation factor corresponding to the nth multipath signal of the kth source; Δ f_k,nRepresents the doppler shift; Δ t_k,nRepresents the propagation delay; tau is_k,nIndicating the reception delay of the mth array element relative to the reference array element; c is the propagation velocity of the sound wave in water; n is_m(t) is the stationary zero mean additive white gaussian noise received by the mth array element; s_m(t) is the multipath signal received by the mth array element; k is the number of far field narrow band sources;is the center frequency corresponding to the kth source;

the array output vector Y (t) is represented as:

Y(t)＝[y₁(t) y₂(t) … y_M(t)]^T

＝A(Θ)BS(t)+N(t)

wherein, a (Θ) represents a direction matrix of the array, B represents a diagonal matrix, s (T) is a signal source vector at time T, T represents a transposition of the matrix, and n (T) is a noise vector of the array.

3. The method for estimating the underwater multipath signal parameters based on the space time-frequency distribution according to claim 1, wherein the step (2) is specifically as follows:

carrying out smooth pseudo-Wegener transformation on the array output vector Y (t) to obtain a space time-frequency distribution matrix D_YY(t, f), expressed as:

D_YY(t,f)＝ΨD_SS(t,f)Ψ^H+²I_M

wherein, t and f respectively represent the index of time and frequency, and (t, f) represent time frequency points, each time frequency point corresponds to a spatial time frequency distribution matrix D with M-X-M dimensions_YY(t, f); Ψ — a (Θ) B, a (Θ) representing the orientation matrix of the array, and B representing the diagonal matrix; (.)^HRepresents the complex conjugate transpose of the vector,²is a white noise variance, I_MRepresenting an identity matrix; d_SS(t, f) is a time-frequency distribution matrix of the multipath signal, expressed as:

diagonal line element

4. The method for estimating the underwater multipath signal parameters based on the space time-frequency distribution according to claim 1, wherein the step (4) is specifically as follows:

to space time frequency distribution matrix D_YY(t, f) carrying out pre-whitening treatment to obtain a whitening space time-frequency distribution matrix

where W is a whitening matrix, U ═ W Ψ denotes a unitary matrix, Ψ ═ a (Θ) B, a (Θ) denotes a direction matrix of the array, and B denotes a diagonal matrix, (W Ψ)^H＝UU^HI denotes an identity matrix;

according to whiten nullTime-frequency distribution matrixSelecting self-source time-frequency points from the set omega according to the following criteria:

wherein the content of the first and second substances,₂for the threshold value, tr {. is } represents the trace of the matrix, and the selected self-source time frequency point (t, f) is stored to the set omega_A。

5. The method for estimating the underwater multipath signal parameters based on the space time-frequency distribution according to claim 1, wherein the step (6) is specifically as follows:

step (6.1): setting parameters:

constructing a zero matrix N _ tf ═ zeros (1, N) for storing the number of self-source time frequency points contained in each ridge line;

construction of N empty sets omega_nN is more than or equal to 1 and less than or equal to N, and all self-source time-frequency points contained in each ridge line are stored;

construction of N empty collections F_nN is more than or equal to 1 and less than or equal to N, and the instantaneous frequency corresponding to the self-source time-frequency point contained in each ridge line is stored;

initializing time slice T ═ T/2, T represents the total number of samples of the signal;

step (6.2): solving the spatial time-frequency distribution matrix D_YY(t, f) N at time t_tInstantaneous frequency f corresponding to each extremum point_t,nWherein (0. ltoreq. F. ltoreq.F_max)：

Storing the corresponding time-frequency and instantaneous frequency, i.e. (t, f)_t,n)∈Ω_nAnd f_t,n∈F_nAnd N _ tf (N) ═ N _ tf (N) + 1; n is a radical of_tIndicating the number of multipath signals involved at time t, initializing N_t＝N；

Are respectively paired with N_tSpatial time-frequency distribution matrix D of self-source time-frequency points corresponding to extreme points_YY(t,f_t,n)，(1≤n≤N_t) Performing characteristic decomposition:

D_YY(t,f_t,n)＝V(t,f_t,n)Λ((t,f_t,n)V(t,f_t,n)^H

wherein, V (t, f)_t,n) Is a unitary matrix of M × M whose column vectors represent eigenvectors, Λ (t, f)_t,n) A diagonal matrix of M × M, wherein the diagonal elements are the eigenvalues corresponding to all eigenvectors, and the eigenvector corresponding to the largest eigenvalue is represented as v (t, f)_t,n) Is a reaction of N_tThe feature vector corresponding to the largest feature value is unitized to obtain the unitized feature vector, i.e. the feature vector

Step (6.3): forward search:

for each time h, wherein h is t-1, t-2, …,1, firstly, it is determined whether the self-source time-frequency points corresponding to all the extreme points at the time h +1 are edge time-frequency points, and if so, f ∈ (f ∈) is selected (f)_h+1,n-ΔF,f_h+1,n+ Delta F), wherein N is more than or equal to 1 and less than or equal to N_h+1And Δ F is a frequency range allowing searching, and the corresponding frequency points (h +1, F) all satisfy:

D_YY(h+1,f)＝0

then the self-sourced time-frequency point (h +1, f)_h+1,n) The target number N of the multipath signals at the time point of the left edge and the time h_h＝N_h+1-1；

Solving for D at h time_YY(h,f)，(1≤f≤F_max) N of (A)_hFrequency f corresponding to each extremum point_h,j，(1≤j≤N_h) According to the j (j is more than or equal to 1 and less than or equal to N) at the h moment_h) Spatial time-frequency distribution matrix D of self-source time-frequency points corresponding to extreme points_YY(h,f_h,j) Solving the corresponding unitized feature vector asAnd each of the h +1 time instantsFeature vectorBy comparison, where (1. ltoreq. N. ltoreq.N_h+1) If the following conditions are met:

the self-source time-frequency point corresponding to the extreme point belongs to the nth ridge line, i.e., (h, f)_h,j)∈Ω_nAnd f_h,j∈F_nAnd the number of self-source time-frequency points contained in the nth ridge line is N _ tf (N) ═ N _ tf (N) + 1;

step (6.4): and (3) backward searching:

for each time h, wherein h is T +1, T +2, …, T, judging whether the time frequency points corresponding to all the extreme points at the time h-1 are edge time frequency points, and if so, f ∈ (f ∈)_h-1,n-ΔF,f_h-1,n+ Delta F), wherein N is more than or equal to 1 and less than or equal to N_h-1And the corresponding time frequency points (h-1, f) all meet the following conditions:

D_YY(h-1,f)＝0

the self-sourced time-frequency point (h-1, f)_h-1,n) Is the right edge time frequency point, the target number N of the h moment_h＝N_h-1-1；

Solving for D at h time_YY(h,f)，(1≤f≤F_max) N of (A)_hFrequency f corresponding to each extremum point_h,j，(1≤j≤N_h) According to the j (j is more than or equal to 1 and less than or equal to N) at the h moment_h) Spatial time-frequency distribution matrix D of self-source time-frequency points corresponding to extreme points_YY(h,f_h,j) Solving the corresponding unitized feature vector as

the self-source time-frequency point corresponding to the extreme point belongs to the nth ridge line, i.e., (h, f)_h,j)∈Ω_nAnd f_h,j∈F_nAnd the number of self-source time-frequency points contained in the nth ridge line is N _ tf (N) ═ N _ tf (N) + 1;

step (6.5): output omega_nN _ tf (N), and F_nAnd (N is more than or equal to 1 and less than or equal to N), and obtaining self-source time-frequency points, the number of the self-source time-frequency points and the instantaneous frequency corresponding to the self-source time-frequency points contained in each ridge line.

6. The method for estimating the underwater multipath signal parameters based on the spatial time-frequency distribution according to claim 1, wherein the step (7) is specifically as follows:

step (7.1): instantaneous frequencyEstimating:

according to the instantaneous frequency set F corresponding to all self-source time-frequency points contained in each ridge line_n(N is more than or equal to 1 and less than or equal to N), and obtaining an estimated value of the instantaneous frequency corresponding to the nth multipath signalComprises the following steps:

wherein, N _ tf (N) is the number of self-source time-frequency points contained in the nth ridge line, f_iIs the instantaneous frequency, F, corresponding to the ith self-source time-frequency point contained in the nth ridge line_nThe instantaneous frequency set corresponding to all self-source time-frequency points contained in the nth ridge line;

step (7.2): orientation

calculating the average value of the space time-frequency distribution matrix corresponding to all the self-source time-frequency points contained in each ridge line, wherein the formula is as follows:

wherein the content of the first and second substances,is the average space time-frequency distribution matrix, omega, corresponding to the nth ridge line_nAll self-source time-frequency point sets contained in the nth ridge line;

to pair

wherein, ∑_SDiagonal matrix of large eigenvalues, ∑_NDiagonal matrix, U, composed of small eigenvalues_SIs a signal subspace, U_NIs a noise subspace; (.)^HRepresents a complex conjugate transpose of the vector;

and (3) carrying out orientation estimation by adopting a TF-music algorithm, and defining a spatial spectrum function:

wherein a (θ) is the direction vector of the array;

estimation value of incident angle of nth multipath signal

where argmax (·) is the argmax function.

Technical Field

The invention relates to the technical field of underwater multipath signal parameter estimation, in particular to an underwater multipath signal parameter estimation method based on space time-frequency distribution.

Background

The research of the underwater direction estimation technology plays an important role in the aspects of marine information acquisition, environment monitoring, resource development, underwater disaster early warning, marine rights and interests maintenance, safety defense and the like.

The ocean underwater acoustic channel is a time-varying-space-varying-frequency-varying random multipath transmission channel, which has the serious problems of strong background noise, multiple interference sources, narrow available bandwidth, large transmission delay, large propagation loss, multipath effect, Doppler effect and the like, and is the currently known wireless communication channel with the greatest development difficulty. The underwater acoustic signals are subjected to a great amount of reflection and refraction after being transmitted through the ocean underwater acoustic channel, and a receiving end can receive signals of direct paths transmitted by an information source and also can receive other reflected and refracted multipath signals, so that great challenges are brought to underwater target direction estimation. Therefore, researching how to realize the azimuth estimation of the multipath signals under the passive scene only according to the array received signals is always a research hotspot and difficulty in the field, and has very important significance to the application fields of marine information such as underwater communication, remote measurement and remote control, sonar and the like.

The azimuth estimation is one of the research hotspots in the field of array signal processing for a long time, and a great deal of excellent research results such as a relatively typical MUSIC algorithm, an ESPRIT algorithm, a maximum likelihood estimation, a subspace algorithm and various breakthrough improvement algorithms thereof have emerged at home and abroad at present, and a very good effect is obtained, and a set of complete and comprehensive theoretical system is gradually established. Although the current research on the direction estimation of stationary narrowband signals has achieved great results, the direction estimation of non-stationary signals has a breakthrough, and the results are relatively few. The traditional azimuth estimation method has the following defects in the azimuth estimation for processing multipath signals in a passive scene:

(1) the traditional orientation estimation method assumes that signals are a stable Gaussian process with the statistical property not changing along with time, namely linearity, stability and Gaussian. However, underwater multipath signals are essentially complex multi-modal, non-stationary signals formed by superimposing a plurality of single-modal signals of different frequencies, different amplitudes, and different time delays. The time delay of all multipath signals to reach the receiving end is different, and the frequency components of the multipath signals are changed along with the time. Therefore, the conventional azimuth estimation method is not suitable for azimuth estimation of multipath signals.

(2) The traditional azimuth estimation method utilizes space-time statistical information of array signals, and does not fully mine and utilize time-frequency information of the signals, so that a lot of useful information is inevitably lost, and the traditional azimuth estimation method has the defect of difficult avoidance on the azimuth estimation problem of processing non-stationary signals.

(3) The conventional azimuth estimation method usually requires the number of array elements to be larger than the number of source signals, i.e. the signal environment should satisfy an over-determined condition. However, due to the fact that the actual underwater environment is complex and changeable, under the condition that an information source radiation signal and an underwater acoustic propagation channel are unknown, after the information source radiation signal reaches a receiving end through multipath propagation, the number of multipath signals received by an array is often unknown and time-varying, and the number of receiving array elements is limited, so that the situation that the number of multipath signals received by the array is larger than the number of array elements is often caused, namely, azimuth estimation under an underdetermined condition is caused, and the traditional azimuth estimation method is not suitable for use.

In addition, accurate estimation of the multipath signal parameters has strong theoretical and practical significance for inverting an underwater sound propagation channel and assisting a target to complete positioning and tracking. Therefore, it is necessary to solve the problem of multipath signal parameter estimation under the implementation of the underdetermined condition in the passive scenario.

Disclosure of Invention

The invention provides an underwater multipath signal parameter estimation method based on space time-frequency distribution, aiming at the defects of the existing multipath signal azimuth estimation technology under the underdetermined condition in the passive scene. The time-frequency analysis method adopts a time-frequency analysis technology, is different from the traditional Fourier transform, not only can effectively extract frequency domain information of the signal, but also can obtain the change characteristic of the frequency along with the time, and further performs combined analysis on the time domain and the frequency domain of the signal. The invention combines the traditional array signal processing technology based on space-time processing with time-frequency analysis to obtain the space time-frequency distribution matrix of the signals. The method expands a time domain research object to a time-frequency domain, namely, a signal in the time-frequency domain selects an effective time-frequency point set through time-frequency filtering, the time-frequency point set after the time-frequency filtering can be simplified into a stable narrow-band signal, a space time-frequency distribution matrix formed by the time-frequency points is used for replacing a traditional received signal covariance matrix, and then the orientation estimation is completed by using a space spectrum estimation method. The method can separate signals which are difficult to separate in a time domain and a frequency domain, has good spatial resolution and system signal processing capacity, has good robustness to noise and interference, can realize accurate estimation of the azimuth under an underdetermined condition under the condition that the number of multipath signals is unknown and time-varying in a passive scene, can accurately estimate the number and instantaneous frequency of the multipath signals, and is easy to popularize to multipath signal parameter estimation under a simpler active scene.

In order to solve the technical problems, the invention adopts the following technical scheme:

an underwater multipath signal parameter estimation method based on space time-frequency distribution comprises the following steps:

step (1): at the time t, an underwater multipath signal is received by using a uniform linear array consisting of M omnidirectional array elements, and an array output vector Y (t) at the time t is obtained₁(t) y₂(t) … y_M(t)]^TWherein, y_m(t) represents the output of the m-th array element;

step (2): solving a spatial time-frequency distribution matrix D of an array output vector Y (t) by adopting smooth pseudo-wigner distribution_YY(t,f)；

And (3): for D_YYSetting a noise threshold value for each sampling time slice in the (t, f), filtering all noise time frequency points on a time frequency plane, and storing the time frequency points selected by the noise threshold value to a set omega;

and (4): for D_YY(t, f) carrying out pre-whitening treatment to obtain a whitening space time-frequency distribution matrixAccording toDesigning a screening criterion, screening self-source time-frequency points from the set omega, and storing the self-source time-frequency points to the set omega_A；

And (5): will set omega_AThe spatial time-frequency distribution matrixes corresponding to all self-source time-frequency points in the three-order tensor form a third-order tensorL represents the set omega_AThe number of internal self-sourced time frequency points; then converting the third order tensor P into an expanded two-dimensional space time-frequency distribution matrixSolving the rank of the matrix R by performing singular value decomposition on the matrix R to obtain the number N of signal sources;

and (6): ridge line detection: extracting time-frequency ridges corresponding to each multipath signal by adopting forward search and backward search to obtain self-source time-frequency points, the number of the self-source time-frequency points and the corresponding instantaneous frequency of the self-source time-frequency points contained in each ridge;

and (7): implementing multipath signal instantaneous frequency according to ridge line detection resultAnd orientationIs estimated.

Compared with the prior art, the invention has the beneficial effects that:

the invention combines the traditional array signal processing technology with time-frequency distribution to further improve the signal processing capability of the system.

1) The invention combines the traditional array signal processing technology based on space-time processing with time-frequency analysis to obtain the space-time-frequency distribution matrix of the signals. The time-frequency analysis technology is different from the traditional Fourier transform, and the time-frequency analysis method not only can effectively extract the frequency domain information of the signal, but also can obtain the change characteristic of the frequency along with the time, so as to carry out the joint analysis on the time domain and the frequency domain of the signal. The method expands a time domain research object to a time-frequency domain, can separate signals which are difficult to separate in the time domain and the frequency domain, has very good spatial resolution and signal processing capacity of a system, and has very good robustness to noise and interference.

2) The method adopts smooth pseudo-wigner distribution in secondary time frequency distribution to reduce cross terms, obtains higher time frequency resolution, and reduces the influence of noise and cross terms on estimation performance by setting a noise threshold and selecting self-source time frequency points. The number of unknown multipath signals is estimated only according to array receiving signals under a passive scene by constructing an expanded two-dimensional space time-frequency distribution matrix and singular value decomposition.

3) The invention fully utilizes the energy distribution of the signals provided by the signal time-frequency representation in the time-frequency domain, introduces the thought of forward and backward smoothing to ensure the local smoothness of the ridge line, and separates different multipath signal components by utilizing the characteristics that all self-source time-frequency points belonging to the same multipath signal have the same principal eigenvector to obtain the time-frequency point set of the ridge line corresponding to each multipath signal.

4) The method realizes accurate estimation of the instantaneous frequency and the direction of the multipath signal according to the ridge line detection result. The method replaces the traditional covariance matrix with the time-frequency distribution matrix, processes the multipath signal components with different time-frequency characteristics one by one, can realize the azimuth estimation under the underdetermined condition under the condition that the number of the multipath signals is unknown and time-varying in the passive scene, effectively improves the accuracy of the azimuth estimation, accurately estimates the number and instantaneous frequency of the multipath signals, and is easy to popularize to the multipath signal parameter estimation under the simpler active scene.

Drawings

In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the embodiments or the prior art will be briefly described below, and it is obvious that the drawings in the following description are only some embodiments of the present invention, and for those skilled in the art, other drawings can be obtained according to these drawings without creative efforts.

FIG. 1 is a flow chart of an underwater multipath signal parameter estimation method based on space time-frequency distribution according to the present invention;

FIG. 2 is a schematic diagram of an underwater multipath signal received by a reference array element in an embodiment of the present invention;

fig. 3 is a schematic time-frequency distribution diagram of underwater multipath signals in the embodiment of the invention.

Fig. 4 is a ridge line schematic diagram of an underwater multipath signal detected in an embodiment of the invention.

FIG. 5 is a comparison of the performance of underwater multipath signal instantaneous frequency estimation based on spatial time-frequency distribution in the embodiment of the present invention.

FIG. 6 is a spatial spectrum of underwater multipath signal orientation estimation based on spatial time-frequency distribution in an embodiment of the present invention.

FIG. 7 is a comparison of underwater multipath signal orientation estimation performance based on spatial time-frequency distribution in the embodiment of the present invention.

Detailed Description

The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.

An underwater multipath signal parameter estimation method based on space time-frequency distribution is disclosed, the basic flow is shown in figure 1, and the method mainly comprises the following steps:

s1: at the time t, an underwater multipath signal is received by using a uniform linear array consisting of M omnidirectional array elements, and an array output vector Y (t) at the time t is obtained₁(t) y₂(t) … y_M(t)]^TWherein, y_m(t) represents the output of the m-th array element;

s2: solving array outputs using smooth pseudo-wigner distributionsA space time-frequency distribution matrix D of the output vector Y (t)_YY(t,f)；

S3: for D_YYSetting a noise threshold value for each sampling time slice in the (t, f), filtering all noise time frequency points on a time frequency plane, and storing the time frequency points selected by the noise threshold value to a set omega;

s4: for D_YY(t, f) carrying out pre-whitening treatment to obtain a whitening space time-frequency distribution matrixAccording toDesigning a screening criterion, screening self-source time-frequency points from the set omega, and storing the self-source time-frequency points to the set omega_A；

S5: will set omega_AThe spatial time-frequency distribution matrixes corresponding to all self-source time-frequency points in the three-order tensor form a third-order tensor

L represents the set omega_AThe number of internal self-sourced time frequency points; then converting the third order tensor P into an expanded two-dimensional space time-frequency distribution matrixSolving the rank of the matrix R by performing singular value decomposition on the matrix R to obtain the number N of signal sources;

s6: ridge line detection: extracting time-frequency ridges corresponding to each multipath signal by adopting forward search and backward search to obtain self-source time-frequency points, the number of the self-source time-frequency points and the corresponding instantaneous frequency of the self-source time-frequency points contained in each ridge;

s7: implementing multipath signal instantaneous frequency according to ridge line detection result

And orientation

Is estimated.

In one embodiment of the present invention, a process for receiving a multipath signal is described.

At time t, using uniform linear array composed of M omnidirectional array elements to receive multipath signals to obtain array output vector Y (t) [ < y > ]₁(t) y₂(t) … y_M(t)]^TWherein the output of the m-th array element is y_m(t) of (d). The method comprises the following specific steps:

s11, K uncorrelated far-field narrow-band sources S in assumed space_k(t), (K ═ 1,2, …, K), where the kth source corresponds to a center frequency f_k ^c. Due to the influence of multipath propagation conditions in an underwater environment, it is assumed that a signal transmitted by a kth source passes through N_kThe strip path is incident on a uniform linear array consisting of M omnidirectional array elements, and N corresponding to the kth information source_kThe angle of incidence of the strip multipath signals is respectivelyWherein theta is_k,n∈[-π,π]，(1≤n≤N_k). It is assumed that all multipath signals of different sources are disjoint, and each multipath signal can be considered as a signal generated by a virtual far-field narrowband source. The number N of total multipath signals corresponding to K sources is expressed as:

s12, setting the distance of array elements in the array as d, taking array element 1 as a reference array element, considering the influence of additive white Gaussian noise, and outputting y of the mth array element at t moment_m(t) is expressed as:

wherein, a_k,nRepresenting the attenuation factor corresponding to the nth multipath signal of the kth source; Δ f_k,nRepresents the doppler shift; Δ t_k,nRepresents the propagation delay;indicating the reception delay of the mth array element relative to the reference array element; c is the propagation velocity of the sound wave in water; n is_m(t) is the stationary zero mean additive white Gaussian noise received by the mth array element, i.e., satisfies E [ n ]_m(t)]0, there is no correlation between the signal and the noise.

S13, at time t, the M × 1 dimensional array output vector y (t) can be expressed as:

wherein a (Θ) ═ a₁(Θ) A₂(Θ) … A_M(Θ)]^TA matrix of directions representing the array is shown,

wherein A is_m(Θ)＝[α(θ₁) α(θ₂) … α(θ_K))，

Is the array direction vector corresponding to the kth source,representing the center frequency corresponding to the nth multipath signal of the kth information source;

＝diag[_{1 2}…_K]wherein, in the step (A),

B＝diag[B₁B₂… B_K]wherein, in the step (A),

S(t)＝[S₁(t) S₂(t) … S_K(t)]^Tis a signal source vector, wherein the multipath signals corresponding to the k-th source are respectively expressed asWherein the content of the first and second substances,

t represents the transpose of the matrix;

N(t)＝[n₁(t) n₂(t) … n_M(t)]^Tis the noise vector of the array.

In one embodiment of the invention, the spatial time-frequency distribution matrix D is_YYThe solving process of (t, f) is introduced.

Solving a spatial time-frequency distribution matrix D of an array output vector Y (t) by adopting smooth pseudo-wigner distribution_YY(t, f), the concrete steps are as follows:

the spatial time-frequency distribution matrix of the array output vector y (t) is represented as:

wherein, t and f respectively represent time and frequency indexes, and (t, f) represent time frequency points, each time frequency point corresponds to a space time frequency distribution matrix D_YY(t, f), l and tau respectively represent time shift and frequency shift, phi (l, tau) is a kernel function of time-frequency distribution, and different phi (l, tau) is adopted to obtain different time-frequency distributions. (.)^HRepresenting the complex conjugate transpose of the vector.

The spatial time-frequency distribution matrix of the array output vector y (t) obtained by bringing formula (3) into formula (4) is represented as:

D_YY(t,f)＝ΨD_SS(t,f)Ψ^H+ΨD_SN(t,f)+D_NS(t,f)Ψ^H+D_NN(t,f) (5)

wherein there is D for each time-frequency point (t, f)_YY(t,f)∈C^M×MWherein

Psi ═ a (Θ) B contains the spatial information of the signal and the information of propagation delay, amplitude, attenuation, etc., and maps the self-source time-frequency distribution matrix and the mutual-source time-frequency distribution matrix in the time-frequency distribution matrix of the signal to the spatial time-frequency distribution matrix of the array output vector,the space time frequency distribution matrix reflects the space and time frequency information of the signal at the same time; d_SS(t, f) is the time frequency distribution matrix of the multipath signal, the invention assumes that N multipath signals are incident into the uniform linear array, therefore, for each time frequency point (t, f), D_SS(t, f) is a matrix of N × N, expressed as:

wherein the content of the first and second substances,diagonal line element thereofSelf-sourced time-frequency distribution matrix, off-diagonal elements, representing multipath signals

Expressed as a mutual source time-frequency distribution matrix; d_SN(t, f) and D_NS(t, f) respectively representing a mutual source time-frequency distribution matrix between the source and the noise, D_NN(t, f) represents a time-frequency distribution matrix of the noise. Assuming that there is no correlation between the signal and the noise component and the noise mean is zero, it can be derived that the mutual source time-frequency distribution matrix between the signal and the noise vector is zero, and therefore, the spatial time-frequency distribution matrix output by the array is represented as:

wherein the content of the first and second substances,²is a white noise variance, I_MRepresenting an identity matrix.

In a specific implementation of the invention, a denoising process of time-frequency points and a screening process of self-sourced time-frequency points are introduced.

And S3, setting a noise threshold.

In order to reduce the influence of noise, the embodiment filters D by setting a noise threshold_YYTime-frequency point corresponding to noise in (t, f), and comparing D_YYTime-frequency points with enough residual signal energy in (t, f) are stored in the set omega. The method comprises the following specific steps:

the energy of the signal in the time-frequency plane is relatively concentrated, while the energy of the noise is low and is generally uniformly distributed over the entire time-frequency plane. Thus, can be represented by D_YYEach sampled time slice t in (t, f)_iAnd setting a noise threshold to judge whether each time frequency point is a noise point, further filtering all the noise time frequency points (t, f) on the time frequency plane, and only keeping the time frequency points (t, f) with enough signal energy.

Wherein | · | purple sweet_FFor the Frobenius norm, i and j are time and frequency indices, respectively,₁the threshold value is related to the noise level (generally 0.05 when the signal-to-noise ratio SNR is 10 dB), the set Ω represents the set of time frequency points (t, f) selected by the noise threshold value, and the robustness of the azimuth estimation algorithm to noise can be effectively enhanced by removing the noise points.

And S4, selecting self-source time frequency points.

The performance of non-stationary signal orientation estimation depends on how accurately to select signal time frequency points to a great extent, the spatial time frequency distribution of the array received signals comprises self-source time frequency points generated by the self-source time frequency distribution of the signals and mutual source time frequency points generated by mutual source time frequency distribution among the signals, and the mutual source time frequency points can seriously influence the performance of signal orientation estimation. Therefore, the invention needs to optimize the time-frequency distribution of the signals to remove the mutual source time-frequency points between the signals and select the self-source time-frequency points corresponding to the signals.

First, this embodiment defines a spatial time-frequency distribution matrix D of the whitening matrix W to the array output vector y (t)_YY(t, f) performing pre-whitening processing, and then selecting self-source time frequency points (t, f) from the set omega by setting a specific criterion_aRemoving mutual source time frequency points (t, f)_cAnd storing the selected self-source time frequency point to a set omega_A. The method comprises the following specific steps:

s41, the covariance matrix of the array received signal y (t) is expressed as:

wherein R is_SCovariance matrix representing signal source vector an N × M dimensional whitening matrix W ∑ is defined^-1/2V^HWherein ∑ and V each represent the pair R_YFor any whitening matrix W, there is a unitary matrix U of dimension N × N such that U equals W Ψ, i.e., (W Ψ)^H＝UU^HI, where I denotes an identity matrix.

S42 is represented by the following formula (10) in the formula_YY(t, f) the left side and the right side are respectively multiplied by the whitening matrix W to obtain a whitening space time-frequency distribution matrix of the array output vector Y (t)

S43 random mutual source time frequency point (t, f)_cAnd the corresponding whitening space time-frequency distribution matrix meets the following conditions.

Wherein, the subscript c represents that the time frequency point is a mutual source time frequency point, and tr {. DEG } represents the trace of the matrix.

S44, selecting self-source time frequency points (t, f) from the set omega according to the following criteria_aAnd the subscript a indicates that the time frequency point is a self-source time frequency point.

Wherein the threshold value₂Is less than 1 empirical threshold (typically when SNR is 10dB ═，₂0.85), the selected self-source time-frequency point is stored to the set omega_A。

In one embodiment of the present invention, a process for estimating the number of multipath signals is described.

S5, estimation of the number of multipath signals.

This embodiment will aggregate_AThe spatial time-frequency distribution matrixes corresponding to all self-source time-frequency points in the three-order tensor form a third-order tensorThen converting the third order tensor P into an expanded two-dimensional space time-frequency distribution matrixAnd estimating the number N of signal sources by performing singular value decomposition on the matrix R. The method comprises the following specific steps:

s51, assuming self-source time frequency point set omega_AThe total number of the time frequency points is L, namely:

Ω_A＝{(t,f)_l|1≤l≤L} (14)

wherein, subscript L represents time frequency point index, and L spatial time frequency distribution matrixes corresponding to L self-source time frequency points respectively represent:

D_YY(t,f)_l＝ΨD_SS(t,f)_lΨ^H,(1≤l≤L) (15)

s52, constructing a third order tensor by L space time frequency distribution matrixes corresponding to L self-source time frequency pointsWherein the (i, j, l) th element P of P_ijlExpressed as:

P_ijl＝[D_YY(k,f)_l]_ij,(1≤i,j≤M,1≤l≤L) (16)

s53, constructing an expanded two-dimensional space time-frequency distribution matrix according to the third-order tensor PAnd satisfy R_(i-1)M+j,l＝P_ijl。

S54, defining matrix

C_l,n＝[D_SS(k,f)_l]_nnThen the matrix R can be converted into:

R＝(Ψ⊙Ψ^*)C^T(17)

in the formula (I), the compound is shown in the specification,wherein, (. cndot.,. ⊙ andrespectively representing the complex conjugate, the Khatri-Rao product and the Kronecker product.

S55, solving the rank of the matrix R to estimate the total number N of the multipath signals in the received signal Y (t).

In one embodiment of the present invention, a ridge detection process is described.

And S6, ridge line detection.

In fact, in the time-frequency plane where the time-frequency aggregation is high, the energy of the signal is mainly concentrated in the region called the ridge line. Therefore, the ridge line of the time-frequency distribution not only can represent the change situation of the signal instantaneous frequency, but also contains key information representing the signal characteristics. And, the present invention assumes that all multipath signals are disjoint in the time-frequency domain, i.e., each multipath signal component corresponds to a ridge. Therefore, as long as the ridge line belonging to each multipath signal is detected, the time-frequency point set belonging to each multipath signal can be obtained, and the accurate estimation of the instantaneous frequency and the direction of the multipath signal is further realized.

Firstly, the invention fully considers the energy distribution of the signal provided by the signal time-frequency representation in the time-frequency domain, introduces the idea of forward and backward smoothing, and searches all local extreme points on the basis of ensuring the local smoothness of the detected ridge line at the adjacent moment. Secondly, the time frequency point set of the ridge line corresponding to each multipath signal is obtained by utilizing the property that all self-source time frequency points belonging to the same multipath signal have the same principal eigenvector. The method comprises the following specific steps:

and S61, setting parameters. The invention assumes a sampling frequency f of the system_sOn the time axis, the sampling interval is Δ t equal to 1/f_sThe total number of samples of the signal is T, and on the frequency axis, the spectral line interval is delta f-f_s/N, maximum search frequency F_max＝f_sAssume that the number of multipath signals contained at time t is N_tWherein T is more than or equal to 0 and less than or equal to T.

Constructing a zero matrix N _ tf ═ zeros (1, N) for storing the number of self-source time frequency points contained in each ridge line;

construction of N empty sets omega_n＝[]N is more than or equal to 1 and less than or equal to N, and all self-source time-frequency points contained in each ridge line are stored;

construction of N empty collections F_n＝[]N is more than or equal to 1 and less than or equal to N, and is used for storing the instantaneous frequency corresponding to all self-source time-frequency points contained in each ridge line;

construction of N empty collections E_n＝[](1. ltoreq. N. ltoreq.N) for saving the set F_nAll the unit characteristic vectors corresponding to the self-source time frequency points (t, f) contained in the table

And S62, initializing.

S621, the present invention initializes time slice T to T/2, and assumes that the multipath signal can obtain the maximum source number N at the intermediate time T, that is, N_t＝N；

S622, solving D_YY(t, f) N at time t_tInstantaneous frequency f corresponding to each extremum point_t,nWherein (0. ltoreq. F. ltoreq.F_max)：

Storing the corresponding time-frequency and instantaneous frequency, i.e. (t, f)_t,n)∈Ω_nAnd f_t,n∈F_nAnd N _ tf (N) ═ N _ tf (N) + 1;

s623, respectively aligning N_tSpatial time-frequency distribution matrix D of self-source time-frequency points (t, f) corresponding to each extreme point_YY(t,f_t,n)，(1≤n≤N_t) Performing characteristic decomposition:

D_YY(t,f_t,n)＝V(t,f_t,n)Λ(t,f_t,n)V(t,f_t,n)^H(19)

wherein, V (t, f)_t,n) Is a unitary matrix of M × M whose column vectors represent eigenvectors, Λ (t, f)_t,n) The M × M diagonal matrix has its diagonal elements as eigenvalues corresponding to all eigenvectors_t,n) Each feature vector is unitized.

And saving the feature vectors after unitization, i.e.

And S63, forward searching.

S631, for each time h, where h is t-1, t-2, …,1, first, determining whether the self-source time frequency points corresponding to all the extreme points at the time h +1 are edge time frequency points, and further estimating the number N of multipath signals corresponding to the time h_h。

If any of f ∈ (f)_h+1,n-ΔF,f_h+1,n+ Delta F), wherein N is more than or equal to 1 and less than or equal to N_h+1And Δ F is a frequency range allowed to be searched, and the corresponding frequency points (h +1, F) all satisfy:

D_YY(h+1,f)＝0 (21)

then the self-sourced time-frequency point (h +1, f) is considered_h+1,n) Is the left edge time frequency point, the target number N of the h moment_h＝N_h+1-1。

S632, according to equation (18), solving for D at time h_YY(h,f)，(1≤f≤F_max) N of (A)_hFrequency f corresponding to each extremum point_h,j，(1≤j≤N_h)。

S633, obtaining the j (1 ≦ j ≦ N) at the h time point according to the equations (19) and (20)_h) Spatial time-frequency distribution matrix D of self-source time-frequency points corresponding to extreme points_YY(h,f_h,j) The corresponding unitized feature vector isEach feature vector associated with time h +1By comparison, where (1. ltoreq. N. ltoreq.N_h+1) If the following conditions are met:

the self-source time-frequency point corresponding to the extreme point belongs to the nth ridge line, i.e., (h, f)_h,j)∈Ω_nAnd f_h,j∈F_nSaving the unitized feature vectorsAnd the number of self-source time-frequency points included in the nth ridge line N _ tf (N) ═ N _ tf (N) + 1.

And S64, backward searching.

S641 determines, for each time h, where h is T +1, T +2, …, and T, whether the time-frequency points corresponding to all the extremum points at the time h-1 are edge time-frequency points according to step S631, and estimates the number N of multipath signals corresponding to the time h_h。

If any of f ∈ (f)_h-1,n-ΔF,f_h-1,n+ Delta F), wherein N is more than or equal to 1 and less than or equal to N_h-1And the corresponding time frequency points (h-1, f) all meet the following conditions:

D_YY(h-1,f)＝0 (23)

then the self-sourced time-frequency point (h-1, f) is considered_h-1,n) Is the right edge time frequency point, the target number N of the h moment_h＝N_h-1-1。

S642, according to the formula (18), solving for D at the h moment_YY(h,f)，(1≤f≤F_max) N of (A)_hFrequency f corresponding to each extremum point_h,j，(1≤j≤N_h)。

S643, the j (1 ≦ j ≦ N) at the h time is obtained according to the expressions (19) and (20)_h) Spatial time-frequency distribution matrix D of self-source time-frequency points corresponding to extreme points_YY(h,f_h,j) The corresponding unitized feature vector isEach feature vector associated with time h-1By comparison, where (1. ltoreq. N. ltoreq.N_h-1) If the following conditions are met:

the self-source time-frequency point corresponding to the extreme point belongs to the nth ridge line, i.e., (h, f)_h,j)∈Ω_nAnd f_h,j∈F_nSaving the unitized feature vectorsAnd the number of self-source time-frequency points included in the nth ridge line N _ tf (N) ═ N _ tf (N) + 1.

And S65, outputting parameters. Outputs N _ tf, Ω_nN is not less than 1 and not more than N and F_n,(1≤n≤N)。

In one embodiment of the present invention, a multipath signal parameter estimation process is described.

And S7, estimating the parameters of the multipath signal.

The invention carries out the parameter estimation of the multipath signal according to the result of ridge line detection, mainly comprises the instantaneous frequency estimation of the signalAnd orientation estimationThe method comprises the following specific steps:

and S71, instantaneous frequency estimation. According to each ridge line containingInstantaneous frequency set F corresponding to all self-source time-frequency points_n(N is more than or equal to 1 and less than or equal to N), the instantaneous frequency corresponding to the nth multipath signal can be estimated as:

and S72, estimating the azimuth. Calculating the time-frequency point set omega corresponding to each ridge line_nAnd (N is more than or equal to 1 and less than or equal to N), averaging the spatial time-frequency distribution matrixes corresponding to all the time frequency points, and then realizing the azimuth estimation of the multipath signals by using a TF-music algorithm. The method comprises the following specific steps:

and S721, constructing a space time-frequency distribution matrix corresponding to each ridge line. Calculating the time-frequency point set omega corresponding to the nth ridge line_nAnd (N is more than or equal to 1 and less than or equal to N) the average value of the space time-frequency distribution matrix corresponding to all the time-frequency points.

Wherein, N _ tf (N) represents the time-frequency point set omega corresponding to the nth ridge line_nThe number of the time frequency points in (1).

S722, average space time-frequency distribution matrix corresponding to the nth ridge lineAnd (5) performing characteristic decomposition.

Wherein, ∑_SDiagonal matrix of large eigenvalues, ∑_NDiagonal matrix, U, composed of small eigenvalues_SIs a signal subspace, U_NIs the noise subspace.

S723, defining a spatial spectrum function:

wherein a (theta) is a direction vector of the array, is related to the frequency and the direction of arrival of the signal, and is only related to the direction of arrival after the frequency of the signal is determined;

the estimated value of the incident angle of the nth multipath signal is:

where argmax (·) is the argmax function.