Bayesian robust doa estimation method for passive synthetic aperture sonar

Abstract

The application relates to a Bayesian robust DOA estimation method of passive synthetic aperture sonar, which comprises the following steps: obtaining frequency domain towed array receiving data based on a water flow power model and a fast Fourier transform, and carrying out beam domain pretreatment on the towed array receiving data to obtain beam domain receiving data; using a hierarchical probability model to represent the beam domain receiving data; using a variational expectation maximization EM algorithm to iteratively optimize each unknown variable in the hierarchical probability model until convergence; after the variational expectation maximization EM algorithm converges, determining a spectral peak in a spatial spectrum, and sequentially carrying out one-dimensional search in the adjacent space of each angle coarse estimation value to obtain a direction fine estimation value of a signal source. The method solves the problem of performance deterioration of a traditional passive synthetic aperture direction finding method in a random signal source and array shape distortion environment.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of signal processing, and particularly relates to a Bayesian Direction-of-Arrival (DOA) estimation method which is insensitive to source amplitude / phase flicker and array shape distortion. BACKGROUND

[0002] Traditional beamforming methods can realize spatial filtering of array receiving data, and the spatial resolution achieved by the method is determined by the ratio of the wavelength of the incident signal to the array aperture. Therefore, high-resolution direction finding of low-frequency signals depends on a large array aperture, which sharply increases the implementation cost of the array direction finding system, and further limits the application range of conventional direction finding methods. To solve this problem, the synthetic aperture concept of converting "time processing gain" into "spatial processing gain" is proposed to achieve the purpose of synthesizing a virtual large-aperture array by using the motion information of a towed array.

[0003] Traditional passive synthetic aperture algorithms are divided into two categories: array element domain and beam domain. The ETAM algorithm proposed by Stergiopoulos and Sullivan belongs to the array element domain processing method, which coherently superimposes the array signals collected by the towed array at different motion times to obtain a receiving signal of a virtual large-aperture array; Rajagopal and Rao proposed an improved ETAM (METAM) algorithm, which not only solves the defect that the ETAM algorithm can only process a single source, but also further reduces the sampling time required for synthesizing a given length aperture; Kim, Youn and Lee modified the estimation method of the phase correction factor in the ETAM algorithm (i.e., replacing single-point sampling with time averaging after upsampling) to improve the direction finding performance of the algorithm in a low signal-to-noise ratio environment. Yen and Carey proposed a passive synthetic aperture direction finding algorithm in the beam domain, which first performs beamforming on each measurement sample, and then coherently adds the beam outputs of each segment, and documents have verified the superior performance of the algorithm relative to the ETAM algorithm in the array shape distortion environment; Stergiopoulos and Urban proposed a fast Fourier transform synthetic aperture (FFTSA) technology, which can be regarded as an improvement of the above-mentioned beam domain synthetic aperture processing technology to further improve the calculation efficiency; In recent years, documents have pointed out that combining the beam domain processing method with the array element domain processing method can further enhance the adaptability to non-ideal signal environment without reducing the spatial resolution, and even when the towed trajectory is a curve, multiple targets can still be distinguished.

[0004] Traditional passive synthetic aperture direction finding algorithm can only deal with phase stationary signal, that is, the signal shows strong correlation in space and time domain. However, in actual underwater environment, the space-time correlation of the received signal is difficult to guarantee, so that the coherent superposition of each measurement sample cannot be carried out, and the array aperture cannot be expanded. Although the maximum likelihood estimation algorithm (MLE) proposed by Nuttall can deal with phase random signal, the method needs to carry out multi-dimensional search when estimating the direction of multiple sources, and the operation amount is large, which does not have practical engineering application value. In addition, the algorithm does not carry out beam domain preprocessing, and is sensitive to array shape distortion. SUMMARY

[0005] The technical problem solved by the present application is:

[0006] In order to avoid the shortcomings of the prior art, the present application provides a Bayesian robust DOA estimation method of passive synthetic aperture sonar, to solve the problem of performance deterioration of traditional passive synthetic aperture direction finding method in random source and array shape distortion environment.

[0007] In order to solve the above technical problems, the technical scheme adopted by the present application is:

[0008] A Bayesian robust DOA estimation method of passive synthetic aperture sonar, characterized in that it comprises:

[0009] S1: obtaining the frequency domain towed array received data based on the flow power model and fast Fourier transform, and carrying out beam domain preprocessing on the towed array received data to obtain beam domain received data;

[0010] S2: using a hierarchical probability model to represent the beam domain received data;

[0011] S3: using variational expectation maximization (EM) algorithm to iteratively optimize each unknown variable in the hierarchical probability model until convergence;

[0012] S4: after the variational expectation maximization (EM) algorithm converges, determining the spectral peak in the spatial spectrum, and sequentially performing one-dimensional search in the adjacent space domain of each angle coarse estimate value to obtain the direction fine estimate value of the source.

[0013] The further technical scheme of the present application is that S1 comprises:

[0014] S11: generating a distorted array shape in the actual underwater towed environment by using a flow power model, and obtaining a received data signal model according to the geometric relationship of the distorted array shape;

[0015] S12: converting the received signal of the received data signal model to the frequency domain by using fast Fourier transform;

[0016] S13: performing beam domain preprocessing on the converted frequency domain signal to obtain beam domain received data.

[0017] Further technical solutions of the present application: the received data signal model in S11:

[0018] z i = A i (θ) b i + ε i , i = 0, 1,..., L-1

[0019] Wherein, A i (θ) = [a i (θ0)... a i (θ K-1 )], The steering vector of the kth source is represented, wherein represents the steering vector of the kth source, and represents the steering vector of the kth source. And respectively represent the amplitude and phase of the kth source at t i , and ε i represents the noise vector received by the array, and L is the total number of sampling points.

[0020] Further technical solutions of the present application: S12 is specifically:

[0021] Take B time intervals of τ, respectively, and do fast Fourier transform on each segment of observation data, and represent its dependence on τ, to obtain the array received data in the frequency domain:

[0022]

[0023] Wherein, z b represents the superposition of the spectrum coefficient of z i in the frequency band [f0(1-v / c), f0(1+v / c)], c represents the sound speed, f0 represents the signal frequency, is obtained by Fourier transform of the array manifold matrix, b b and respectively represent the Fourier coefficients of b i and ε i in the bth time period.

[0024] Further technical solutions of the present application: S13 is specifically:

[0025] The N×J-dimensional beam forming matrix W is used to convert the high-dimensional element domain data into low-dimensional beam domain data, to obtain:

[0026]

[0027] where T = W(W H W) -1 / 2 , and denote the array manifold matrix and the received noise in the beam domain, respectively.

[0028] The further technical solution of the present application is that S2 is specifically:

[0029] A spatially super-complete array manifold matrix of JxM dimensions is introduced The beam domain data is spatially discretized, where M » K, to obtain:

[0030]

[0031] where is a super-complete representation of b b , that is, the positions of the K non-zero elements in correspond to the directions of the K incident signals;

[0032] The likelihood function of the beam domain received data after spatial sparsification obeys the following complex Gaussian distribution:

[0033]

[0034] where

[0035] The further technical solution of the present application is that S3 is specifically:

[0036] The variational expectation-maximization algorithm is used to iteratively optimize the unknown parameters in the hierarchical probability model; in a single iteration, the partial derivative of the log-likelihood function of the observed data with respect to each unknown variable is maximized until the local maximum of the log-likelihood is converged.

[0037] A computer system, characterized in comprising: one or more processors, a computer readable storage medium for storing one or more programs, wherein when the one or more programs are executed by the one or more processors, the one or more processors implement the above method.

[0038] A computer readable storage medium, characterized in storing computer executable instructions, the instructions being executed to implement the above method.

[0039] A computer program product, characterized in comprising computer executable instructions, the instructions being executed to implement the above method.

[0040] The present application has the beneficial effects that:

[0041] ​The application provides a passive synthetic aperture sonar Bayesian robust DOA estimation method, and has the following beneficial effects:

[0042] (1) Although the passive synthetic aperture technology provides a feasible path for resolving the contradiction between the array physical aperture and the spatial azimuth resolution performance, the ability of the traditional passive synthetic aperture algorithm to realize high-precision direction finding is being challenged by the increasingly complex signal environment. The existing passive synthetic aperture direction finding algorithm estimates the phase correction factor by using the coherence of the received signal, and therefore cannot effectively synthesize the aperture for a single target signal with random phase or even discontinuous phase. However, the non-cooperation of underwater reconnaissance and the time-varying characteristics of the underwater acoustic environment will cause the phase of the array received signal to change rapidly in the time domain, and in the above signal environment, the traditional passive synthetic aperture direction finding algorithm fails. The method designed in the application is based on in-depth analysis of the basic principles of the passive synthetic aperture technology and the shortcomings of related algorithms, and by combining the advantages of the Bayesian sparse reconstruction technology in low signal-to-noise ratio, model mismatch adaptability and super-resolution capability with the actual application requirements of the passive synthetic aperture technology, the method realizes accurate estimation of the azimuth of a phase random target signal in a distorted towed array array shape environment. The method designed in the application can synthesize the virtual aperture by using the continuous measurement output of the array during the movement according to the Bayesian multi-layer probability model, and avoids estimating the phase correction factor by using the cross-correlation of the overlapping subarrays, so that the aperture synthesis effect is not limited by the signal coherence.

[0043] (2) The application requirement of large aperture sonar array urgently needs to obtain the array gain and azimuth resolution of long array by moving small aperture array, thereby reducing the engineering implementation cost, and the existing towed array array shape estimation method does not utilize the small scale motion array to convert time accumulation into spatial gain, so that the aperture shortage problem existing in the passive azimuth estimation of towed array sonar target cannot be solved by using synthetic aperture technology. How to combine the array aperture expansion technology and a small amount of model prior knowledge to inversely obtain the whole array structure is the main problem to be solved in the towed array array shape estimation technology. At present, there is no solution to this problem at home and abroad. The present application provides a robust passive synthetic aperture processing method for towed array which is tolerant to array element position error, and realizes off-grid DOA estimation by using the idea of local grid refinement. Since the actual array cannot be completely calibrated, the beam domain preprocessing method is used to reduce the sensitivity of the array processing algorithm to the array element position error, and then the separation of each correlation signal (caused by the multipath effect and other factors commonly existing in actual application) is realized by means of Bayesian factor analysis technology, so as to avoid the influence of the mismatch between the commonly used independent prior distribution assumption and the actual signal model on the direction finding accuracy. After the convergence of the spatial feature reconstruction process of the incident signal, in order to eliminate the influence of angle quantization error, the spatial search is carried out in the smaller angle interval in each signal spectrum peak interval, and the maximum spectrum peak position of the observation data likelihood function in these intervals corresponds to the signal incident direction, and this angle precision estimation process will not reduce the calculation efficiency of the algorithm. BRIEF DESCRIPTION OF DRAWINGS

[0044] The accompanying drawings are included to provide a further understanding of the application and are incorporated in and constitute a part of this specification, illustrate embodiments of the application and together with the description serve to explain the principles of the application. In the drawings:

[0045] Figure 1 The flow chart of the method of the present application.

[0046] Figure 2 In a random signal environment, the normalized spatial spectrum comparison chart of the algorithm designed by the present application and the ETAM and FFTSA algorithms (the number of signal sources is 2).

[0047] Figure 3 The DOA estimation RMSE change curve of the direction finding algorithm designed by the present application and the traditional passive synthetic aperture direction finding algorithm with SNR (the incident signal source is a random source, and the number of signal sources is 2).

[0048] Figure 4 The DOA estimation RMSE change curve of the direction finding algorithm designed by the present application and the traditional passive synthetic aperture direction finding algorithm with the number of snapshots (the incident signal source is a random source, and the number of signal sources is 2).

[0049] Figure 5The DOA estimation RMSE curve of the designed direction finding algorithm and the traditional passive synthetic aperture direction finding algorithm varies with the angle interval (the incident source is a random source, and the number of sources is 2). DETAILED DESCRIPTION

[0050] In order to make the purpose, technical scheme and advantages of the present application clearer, the present application is further described in detail below in combination with the drawings and examples. It should be understood that the specific examples described herein are only used to explain the present application and do not limit the present application. In addition, the technical features involved in each embodiment of the present application described below can be combined with each other as long as they do not conflict with each other.

[0051] The present application provides a Bayesian robust DOA estimation method for passive synthetic aperture sonar. First, the flow power model and the fast Fourier transform technology are used to obtain the frequency domain towed array receiving data, and the beam domain preprocessing is performed on the data. Second, a hierarchical probability model is used to represent the above-mentioned beam domain receiving data, and the spatial sparsity of the signal incident direction is mined. Third, the variational expectation-maximization (EM) algorithm is used to obtain the maximum a posteriori estimation value of each unknown variable in the above-mentioned probability model. Finally, one-dimensional search is performed in the adjacent spatial domain of the coarse angle estimate to obtain the fine angle estimation result. The method comprises the following steps:

[0052] Step 1: a towed line array with N elements is used, the element spacing is set as d, and it is assumed that the incident signal is K (K < N) narrowband signals with a frequency f0. The water flow power model is used to generate the distorted array shape in the actual underwater environment. If the towed speed of the array is v, the observation vector of each element at t i time can be expressed as follows:

[0053] z i =A i (θ)b i +ε i ,i=0,1,…,L-1

[0054] wherein A i (θ)=[a i (θ0)…a i (θ K-1 )], represents the steering vector of the kth source, wherein i represents and respectively represent the amplitude and phase of the kth source at t i time, and ε i represents the noise vector received by the array.

[0055] Take B time intervals of τ, respectively, for each segment of the observation data to do fast Fourier transform, and characterize its dependence on τ, get the array receiving data in the frequency domain:

[0056]

[0057] Where, z b represents z i The superposition of the frequency spectrum coefficient in the frequency band [f0(1-v / c),f0(1+v / c)], c represents the sound velocity, is the array manifold matrix obtained by Fourier transform, b b and respectively represent b i And ε i The Fourier coefficient in the bth time interval.

[0058] After obtaining the above frequency domain data model, the beam forming matrix W is used to convert the high-dimensional array element domain data into low-dimensional beam domain data, and the following is obtained:

[0059]

[0060] Where, T=W(W H W) -1 / 2 , and respectively represent the array manifold matrix and the received noise in the beam domain.

[0061] Step 2: By introducing a JxM(M>>K) dimensional spatial super complete array manifold matrix The beam domain data is discretized in space to obtain:

[0062]

[0063] Where, is the over complete representation form of b b , that is, The position of the K non-zero elements in corresponds to the direction of the K incident signals.

[0064] The likelihood function of the beam domain receiving data after spatial sparsification obeys the following complex Gaussian distribution:

[0065]

[0066] Where, On this basis, the specific meaning and distribution form of each related variable are described in detail.

[0067] Step 3: Based on the hierarchical probabilistic model, the unknown parameters in the model are iteratively optimized by using the variational EM algorithm. In a single iteration, the partial derivatives of the log-likelihood function of the observed data with respect to each unknown variable are maximized until the local maximum of the log-likelihood is converged.

[0068] Step 4: After the variational EM algorithm converges, the spatial spectrum can be determined , and a one-dimensional search is performed in the vicinity of each angle coarse estimate to obtain the fine estimate of the direction of the source.

[0069] The above steps are as follows:

[0070] Step 1: A flow dynamics model is used to generate the distorted array shape in the actual underwater towing environment, and a model of the received data signal is obtained based on this geometric relationship. Then, the received signal is converted to the frequency domain using fast Fourier transform technology, and the frequency domain signal is preprocessed in the beam domain.

[0071] The flow dynamics model uses a zero-order partial differential equation to describe the propagation characteristics of the horizontal offset of the first array element on the array, and a first-order Euler approximation is used to discretize the partial differential operator of the equation to obtain the following state equation:

[0072] l i =Fl i-1 +ω i +ε i

[0073] The above equation illustrates the relationship between the heading of the first array element at time iΔt and (i-1)Δt, where Δt is the time discretization interval. In the above equation, F represents the transition matrix, ω i represents the driving term, and ε i represents the discretization error of the partial differential operator. The solution of the above discrete equation is obtained by numerical methods, and then converted into the spatial coordinates of the nth hydrophone at time iΔt:

[0074]

[0075] where (x0, y0) is the spatial coordinates of the first array element, which can generally be obtained by measurement equipment, and d is the spatial distance between adjacent hydrophones.

[0076] Assume that K (K < N) single-frequency sound sources are incident on a uniform linear array of N elements, the sound source frequency is f0, the number of array elements is N, the signal sampling time is t i = iΔt, i = 0, 1, …, L-1, where L is the total number of sampling points. Without loss of generality, the first array element is taken as the reference array element, i.e. the first array element coordinates are (x0, y0) = (0, 0). Assuming the towing speed is v, then the nth hydrophone at t iThe sound pressure field of the kth sound source received at time t is:

[0077]

[0078] wherein, and respectively represent the amplitude and phase of the kth sound source at time t i , wherein represents the kth sound source is the time frequency of the kth sound source, ω0=2πf0 is the angular frequency, θ k is the incident direction of the kth sound source, and c represents the sound speed. k represents the spatial frequency of the kth sound source, and the expression is as follows:

[0079]

[0080] wherein λ0=c / f0 represents the signal wavelength.

[0081] The received signals of each array element at time t i are spliced according to the array element number to obtain the following N-dimensional array received signal:

[0082] z i =A i (θ)b i +ε i , i=0, 1, …, L-1

[0083] wherein ε i represents an N×1-dimensional array received noise vector, the elements are independent of each other, and subject to a complex Gaussian distribution with a mean of 0 and a variance of σ 2 , A i (θ)=[a i (θ0)…a i (θ K-1 )]. represents the steering vector of the kth sound source.

[0084] The entire observation time period (L-1)Δt is divided to obtain B observation intervals: 0, (L-1)Δt / B, …, (L-1)Δt, and the duration of each observation interval is Then, the fast Fourier transform is performed on each segment of observation data at the angular frequency ω0 to obtain the following array received data in the frequency domain:

[0085]

[0086] wherein z b represents z ithe superposition of the spectral coefficients in the frequency band [f0(1-v / c), f0(1+v / c)], is the array manifold matrix A i (θ) is obtained by Fourier transform, and the expression of the kth column is:

[0087]

[0088] where b b and are the b i and ε i Fourier coefficients in the b th time interval, respectively. The superscript τ of b and

[0089] implies the dependence of these data on τ.

[0090]

[0091] where j=0, 1, …, J-1, J represents the number of designed spatial beams, and satisfies the following relationship K max and θ min are the upper and lower boundaries of the preset beam illumination area, The beamforming matrix is orthogonalized in the following way:

[0092] T=W(W H W) -1 / 2

[0093] And T H is multiplied on the left of the array element domain data z b to obtain the following beam domain data:

[0094]

[0095] where, and represent the array manifold matrix and the received noise in the beam domain, respectively.

[0096] Step 2: Use the hierarchical probability model to represent the above beam domain received data and mine the spatial sparsity of the signal incident direction.

[0097] Discretize the beam domain data obtained in step 1 in space to obtain:

[0098]

[0099] where, Let J×M (M>>K) dimensional spatial overcomplete array manifold matrix be constructed in the same way as... Similarly, only Replace with preset spatial discrete angles That's all. It is b b The overcomplete representation form, that is, in addition to the actual location of the incident signal, All other elements are 0. In other words, The positions of the K non-zero elements correspond to the orientations of the K incident signals.

[0100] The likelihood function of the spatially sparse beam domain received data follows a complex Gaussian distribution:

[0101]

[0102] in, The prior distribution of each column is modeled as a complex Gaussian distribution:

[0103]

[0104] Where C is an M×D dimensional "factor loading" matrix, and μ is an M×1 dimensional vector containing... The non-zero components of the mean, It is an M×M dimensional diagonal covariance matrix, η b Let b∈{1,…,B} be a D×1 dimensional latent variable vector that follows a complex Gaussian distribution with mean vector μ0 and covariance matrix Σ0:

[0105]

[0106] Using the formulas for probability addition and multiplication, we obtain the following: Marginal distribution:

[0107]

[0108] To simplify the model, let μ0 = 0. Furthermore, assume that the columns in C are statistically independent and follow a complex Gaussian distribution as follows:

[0109]

[0110] Among them, c i Let α represent the i-th column of C, and let α be the hyperparameter. i The precision is represented by the following inverse Gamma distribution:

[0111]

[0112] where shape parameter δ and scale parameter γ are both non-negative, denotes the gamma function. Using p(αi i ), we can obtain the following c i distribution of the edge:

[0113]

[0114] where, is the third kind of modified Bessel function with index ω. The distribution shown in the above formula is a generalized hyperbolic distribution.

[0115] Step 3: Use the variational EM algorithm to optimize the unknown variables.

[0116] Based on the hierarchical probabilistic model, the unknown parameters and the latent variables are iteratively optimized using the variational EM algorithm. In a single iteration process, the algorithm maximizes the partial derivative of the log-likelihood function of the observed data with respect to the posterior distribution of H and Λ (corresponding to the E step and M step of the variational EM algorithm, respectively), until it converges to the local maximum of the log-likelihood.

[0117] Based on the above variational EM algorithm, the iterative update formula of the posterior distribution of the latent variables in the E step is easily derived as:

[0118]

[0119] where,

[0120] denotes the i-th row of C,

[0121] In the M step, the iterative update formula of the model parameters is easily derived as:

[0122]

[0123] where, is called the double gamma function.

[0124] The above parameter estimates can be used in the update of the posterior distribution of the latent variables in the E step in the next iteration process.

[0125] Step 4: Perform a one-dimensional search in the vicinity of the coarse angle estimate to obtain a fine angle estimate.

[0126] After the convergence of the variational EM algorithm, we can determine spectral peaks in the spatial spectrum, and the angle set contained in each spectral peak is Only the ζ=(ζ1,...,ζ M ) T In the value at the corresponding orientation, and the elements of the remaining positions are set to 0, then the log-likelihood function of the observation data with respect to the above truncated ζ vector is:

[0127]

[0128] Where, Indicates the log-likelihood function composed of the remaining terms after removing the terms related to ζ j in ζ For convenience of representation, define:

[0129]

[0130] The refined estimates of the orientation and power of the jth source can be calculated by the following formula:

[0131]

[0132] Thus, the refined estimate of ζ j is:

[0133]

[0134] And substitute this refined estimate into the above joint refined estimate formula of the orientation and power, to obtain the refined estimate formula of θ j :

[0135]

[0136] Where, Ω j is the angular refined search interval of the jth source. Apply the above method to the spectral peaks in turn, and the refined estimate of the orientation of the source can be obtained. The angular refined estimation method proposed in the present application only involves one-dimensional search, so compared with the traditional maximum likelihood estimation method involving multi-dimensional search, it greatly reduces the operation complexity while ensuring high estimation accuracy.

[0137] The effects of the present application can be illustrated by the following simulation:

[0138] 1. Simulation conditions:

[0139] The towed array is a uniform linear array, the number of array elements is 16, the element spacing is half of the wavelength of the incident signal, the incident signal frequency is 100 Hz, the observation time is 100 s, the towed speed is 12 knots, and the signal sampling frequency is 200 Hz, that is, every 250 sampling time points, the array moves half a wavelength.

[0140] 2. Simulation content and results:

[0141] Simulation 1: The incident DOAs of two equal-power random sources are -15.5° and 13.5°, respectively, and the signal SNR is 0dB. In both the FFTSA algorithm and the designed algorithm, each observation interval contains 250 signal time-domain sampling points. Furthermore, in beam domain preprocessing, the beam coverage range is [-20°, +20°], and the number of beams is 10. Figure 2 The image shows a comparison of the normalized spatial spectra of the algorithm designed in this invention, the ETAM algorithm, and the FFTSA algorithm.

[0142] according to Figure 2 The simulation results show that neither the ETAM nor FFTSA algorithms can simultaneously estimate the orientation of two sources, thus proving that these two algorithms fail in random source environments. In contrast, the algorithm designed in this invention can accurately estimate the orientation of all targets in random source environments.

[0143] Simulation 2: 100 Monte Carlo trials were conducted to obtain ETAM, METAM, TD-ETAM, FFTSA, Yen-Carey, Jin-Li, and ETAM+JLZA (MMHyder and K. Mahata, "Direction-of-arrival estimation using a mixed l 2,0 The root mean square error of the DOA estimation algorithm ("norm approximation," IEEE Trans. Signal Process., vol. 58, no. 9, pp. 4646-4655, 2010.) is defined as... In the formula and Let represent the estimated and true values ​​of the DOA of the k-th source in the i-th trial, respectively. Furthermore, the Cramerlow lower bound (CRLB) is used as a comparison benchmark to evaluate the estimation accuracy of the above algorithm. In the simulation, the incident azimuths of the two random signals are -12.5° and 10.5°, respectively, the number of snapshots is 8000, and the other parameters are the same as in Simulation 1. Figure 3 The figure shows the curve of RMSE of DOA estimation as a function of SNR for the above algorithm.

[0144] Figure 3 The results show that the algorithm designed in this invention has the smallest RMSE, while the ETAM, METAM, TD-ETAM, Yen-Carey and ETAM+JLZA algorithms have large direction-finding errors in this simulation environment. The reason for this is that, except for the algorithm designed in this invention, the successful application of the other algorithms all depend on the strong correlation of the received signals during the observation period in order to accurately estimate the phase correction factor, but this simulation environment obviously does not meet the above requirements.

[0145] Simulation 3: The SNR of the two sources in simulation 2 is fixed as 0dB, and the rest of the parameters remain unchanged, the RMSE of each algorithm is plotted with the number of snapshots, as shown in Fig. 3. Figure 4

[0146] Figure 4 The simulation results shown in Fig. 3 show that the algorithm designed in the application can accurately estimate the bearing of the two sources under different snapshots, and the direction finding result error of the rest of the algorithms is larger, which is due to the fact that these algorithms can only be applied to the direction finding environment of coherent signals.

[0147] Simulation 4: The number of snapshots in simulation 3 is fixed as 8000, and the rest of the parameters remain unchanged, the angle interval is changed from 1° to 8°, and the RMSE of each algorithm is plotted with the angle interval, as shown in Fig. 4. Figure 5

[0148] Figure 5 The simulation results shown in Fig. 4 show that in the case of spatially adjacent sources, the direction finding performance of the algorithm designed in the application is still much better than that of the other comparative algorithms.

[0149] The above is only a specific embodiment of the application, but the protection scope of the application is not limited thereto, and any person skilled in the art can easily think of various equivalent modifications or replacements within the technical range disclosed in the application, and these modifications or replacements should be covered within the protection scope of the application.​​

Claims

1. A Bayesian robust DOA estimation method for passive synthetic aperture sonar, characterized in that, The method comprises the following steps: S1: obtaining frequency domain towed array receiving data based on a water flow force model and a fast Fourier transform, and performing beam domain preprocessing on the towed array receiving data to obtain beam domain receiving data; S2: using a hierarchical probability model to represent the beam domain receiving data; S3: using a variational expectation maximization (EM) algorithm to iteratively optimize each unknown variable in the hierarchical probability model until convergence is achieved; S4: after the variational expectation maximization (EM) algorithm converges, determining a spectral peak in a spatial spectrum, and sequentially performing one-dimensional search in a neighboring space of each angle coarse estimate value to obtain an azimuth fine estimate value of the signal source.

2. The Bayesian robust DOA estimation method of a passive synthetic aperture sonar according to claim 1, characterized in that, S1 comprises the following steps: S11: generating a distorted array shape in an actual underwater towed environment using a water flow force model, and obtaining a receiving data signal model according to the geometric relationship of the distorted array shape; S12: converting the receiving signal of the receiving data signal model to the frequency domain using a fast Fourier transform; S13: performing beam domain preprocessing on the converted frequency domain signal to obtain beam domain receiving data.

3. The Bayesian robust DOA estimation method of a passive synthetic aperture sonar according to claim 2, characterized in that, The receiving data signal model in S11 comprises the following steps: z i = A i (θ)b i + ε i i = 0, 1,..., L - 1 where A i (θ) = [a i (θ0)…a i (θ K-1 )], denotes the steering vector of the kth source, where and denote the amplitude and phase of the kth source at time t i , respectively, and i denotes the noise vector received by the array, and L is the total number of samples.

4. The Bayesian robust DOA estimation method of a passive synthetic aperture sonar according to claim 3, characterized in that, S12 specifically comprises the following steps: B time intervals each of which is τ are taken, and fast Fourier transforms are performed on each segment of observation data to represent the dependence on τ, thereby obtaining array receiving data in the frequency domain: where z b represents z i The superposition of the frequency spectrum coefficients in the frequency band [f0(1-v / c),f0(1+v / c)], c represents the sound speed, f0 represents the signal frequency, is obtained by Fourier transform of the array manifold matrix, b b and respectively represent b i and ε i The Fourier coefficient in the bth time period.

5. The Bayesian robust DOA estimation method of a passive synthetic aperture sonar according to claim 4, characterized in that, S13 specifically comprises the following steps: An N×J-dimensional beam forming matrix W is used to convert high-dimensional array element domain data into low-dimensional beam domain data, thereby obtaining: where T = W(W H W) -1 / 2 , and denote the array manifold matrix and the receive noise in the beam domain, respectively.

6. The Bayesian robust DOA estimation method of a passive synthetic aperture sonar according to claim 5, characterized in that, S2 specifically comprises the following steps: Introducing a JxM dimensional spatially super-complete array manifold matrix Spatially discretizing the beam domain data, where M » K, yields: wherein, is a b b overcomplete representation, i.e. the positions of the K non-zero elements in correspond to the bearings of the K incident signals; The likelihood function of the spatially sparse beam domain receiving data is subject to the following complex Gaussian distribution: wherein 7. The Bayesian robust DOA estimation method of a passive synthetic aperture sonar according to claim 6, characterized in that, S3 specifically comprises the following steps: The variational expectation maximization algorithm is used to iteratively optimize unknown parameters in the hierarchical probability model; in a single iteration, the partial derivative of the log-likelihood function of the observation data with respect to each unknown variable is maximized until the local maximum of the log-likelihood is reached.

8. A computer system, characterized by The method comprises the following steps: One or more processors, a computer readable storage medium, and one or more programs stored in the computer readable storage medium, wherein when the one or more programs are executed by the one or more processors, the one or more processors implement the method of claim 1.

9. A computer-readable storage medium, characterized in that Computer executable instructions are stored, and the instructions are used to implement the method of claim 1 when executed.

10. A computer program product, characterised in that Computer executable instructions are stored, and the instructions are used to implement the method of claim 1 when executed.

Images