Bayesian doa estimation method for array distortion passive synthetic aperture sonar

Abstract

The application relates to a Bayesian DOA estimation method of an array distortion passive synthetic aperture sonar and belongs to the technical field of signal processing. The method comprises the following steps: performing layered probability modeling on array observation data, the model is used for applying a binary prior distribution to a signal vector to contain sparse induction characteristics; a variational Bayesian method is used to iteratively maximize the lower bound of the marginal likelihood function of array distortion parameters and hidden variables; unknown parameters are iteratively updated according to the estimated posterior distribution; and a spatial spectrum diagram is drawn according to the optimal estimation result, and each DOA is determined according to a peak value. The method solves the problems that a traditional passive synthetic aperture direction finding algorithm cannot estimate the direction of a random source and is sensitive to array distortion.

Description

Technical Field

[0001] This invention belongs to the field of signal processing technology, and specifically relates to a novel Bayesian passive synthetic aperture method for estimating the direction of arrival (DOA) of underwater sources. Background Technology

[0002] Since the Rayleigh resolution of a linear array is proportional to the ratio of the incident signal wavelength to the array aperture, long arrays with large apertures are required to detect weak signals and resolve nearby targets in the spatial domain. Passive synthetic aperture sonar achieves this by dragging a small-scale array, converting temporal gain into spatial gain to synthesize a long, expensive array.

[0003] Classical passive synthetic aperture algorithms can be broadly classified into two categories: beam-domain algorithms and element-domain algorithms. Typical examples of beam-domain algorithms include the Yen-Carey algorithm (NCYen and W. Carey, "Application of synthetic-aperture processing to towed-array data," J. Acoust. Soc. Amer., vol. 86, no. 2, pp. 754-765, 1989) and the FFTSA algorithm (S. Stergiopoulos and H. Urban, "A new passive synthetic aperture technique for towed arrays," IEEE J. Oceanic Eng., vol. 17, no. 1, pp. 16-25, 1992). This type of algorithm achieves the goal of expanding the physical array aperture by coherently superimposing the beam outputs of sub-apertures.Typical examples of array element domain algorithms include the ETAM algorithm (S. Stergiopoulos and E.J. Sullivan, “Extended towed array processing by overlapped correlator,” J. Acoust. Soc. Amer., vol. 86, no. 1, pp. 158-171, 1989), the METAM algorithm (R. Rajagopal and P. Ramakrishna Rao, “A modified extended towed array method (METAM) for passive synthetic aperture beamforming,” International Symposium on Signal Processing and its Applications, ISSPA, Gold Coast, Australia, 1996), the TD-ETAM algorithm (S. Kim, DH Youn, and C. Lee, “Temporal domain processing for a synthetic aperture array,” IEEE J. Oceanic Eng., vol. 27, no. 2, pp. 322-327, 2002), and the Jin-Li algorithm (S. Jin, Y. Li, and...). H. Huang, “An improved passive synthetic aperture algorithm based on curvilinear maneuverability of autonomous underwater vehicles,” J. Electron. Inf. Technol., vol. 40, no. 9, pp. 2265-2272, 2018. This type of algorithm directly cross-correlates the received data of overlapping array elements in adjacent measurement samples to achieve phase compensation and synthetic aperture. However, the aforementioned passive synthetic aperture direction finding algorithms are all designed based on a linear array receiving signal model, thus their performance degrades significantly when processing received signals from distorted arrays (usually caused by array sway, non-uniform tugboat speed, or uneven mass density distribution of towed arrays). Furthermore, introducing traditional array calibration algorithms fails to utilize array motion information to expand the array aperture.Furthermore, traditional passive synthetic aperture direction finding algorithms cannot be applied to scenarios with multiple random sources. In other words, if at least one of the incident sound sources is a random source, traditional passive synthetic aperture direction finding algorithms will be unable to estimate the location of these sources. This is because the signal coherence between adjacent measurement samples cannot be guaranteed, and therefore the phase correction factor of the synthetic aperture cannot be estimated. Although the maximum likelihood passive synthetic aperture direction finding algorithm proposed by Nuttall et al. (AHNuttall, "The maximum likelihood estimator for acoustic synthetic aperture processing," IEEE J. Oceanic Eng., vol. 17, no. 1, pp. 26-29, 1992) can effectively solve the above problems, this algorithm requires multi-dimensional search, has a huge computational load, and cannot be applied to practical engineering. Summary of the Invention

[0004] The technical problem to be solved by this invention is:

[0005] To avoid the shortcomings of existing technologies, this invention provides a Bayesian DOA estimation method for array distortion passive synthetic aperture sonar, which solves the problems that traditional passive synthetic aperture direction finding algorithms cannot estimate the azimuth of random sources and are sensitive to array distortion.

[0006] To solve the above-mentioned technical problems, the technical solution adopted by the present invention is as follows:

[0007] A Bayesian DOA estimation method for array-distorted passive synthetic aperture sonar, characterized by employing an N-element towed linear array to receive K incident sound sources, including:

[0008] S1: Based on the azimuth of the incident signal, the complex amplitude vector of the incident signal, and the angle between the line connecting the array elements and the baseline, the observation data of the towed array is obtained.

[0009] S2: Define the spatial overcomplete array manifold matrix, and convert the observation data into a sparse signal model based on the spatial overcomplete array manifold matrix;

[0010] S3: Perform hierarchical probability modeling for each latent variable in the sparse signal model, transforming the sparse reconstruction problem into... Norm optimization problem;

[0011] S4: The variational EM algorithm is used to iteratively update the posterior distribution of each latent variable and the maximum likelihood estimate of each model parameter in the hierarchical probability model;

[0012] S5: Using the observed spatial grid points as the abscissa and the mean vector of the iterated sparse signal as the ordinate, plot the spatial amplitude spectrum. Obtain the first K peaks from the amplitude spectrum in descending order of amplitude. The abscissa angle value corresponding to the peak is the incident signal DOA.

[0013] A further technical solution of the present invention: the expression for the observation data in S1:

[0014] z i =A i (θ,ρ)b i +ε i

[0015] in, Let a represent the array manifold matrix at the i-th sampling time. i (θ k ,ρ) represents the azimuth θ corresponding to the k-th incident source. k And the steering vector at the i-th sampling time, ρ=[ρ1,…,ρ N-2 [ ] is the angle between the line connecting each hydrophone segment and the baseline; This represents the complex amplitude vector of the incident signal. and Let ε represent the amplitude and phase of the k-th source at the i-th sampling time, respectively. i =[ε1(t) i ),…,ε N (t i )] T This represents a Gaussian white noise vector.

[0016] A further technical solution of the present invention: The expression of the sparse signal model in S2 is:

[0017]

[0018] Among them, A i (υ,ρ)=[a i (υ0,ρ),…,a i (υ M-1 [,ρ)] is the overcomplete array manifold matrix, a i (υ m ,ρ) represents the angle υ at the i-th sampling time. m The array guide vector, υ=[υ0,υ1,…υ M-1 [] represents the observation space grid points, and M represents the preset number of spatial discrete grid points; It is a sparse vector, constructed as follows: assuming a predefined spatial discrete grid of points υ m-1 If it is not included in the actual DOA set, then The m-th element The other positions correspond The vector element values ​​represent the true amplitude of the signal incident from that direction.

[0019] A further technical solution of the present invention: S3 specifically includes:

[0020] For signal vectors Applying a sparse prior distribution-binary distribution, the likelihood function can then be expressed as:

[0021]

[0022] Among them, z ij Indicate z i The j-th element; binary indicator variable γ m It follows the Bernoulli distribution as follows:

[0023]

[0024] The m-th element of the signal vector It follows the following pattern: mean 0, variance α. m The complex Gaussian distribution:

[0025]

[0026] Integrating the two prior distributions above, the resulting marginal distribution is... Norm;

[0027] Apply the following conditions to the disturbance vector ρ: mean 0, variance ρ Gaussian priors:

[0028]

[0029] noise variance Signal power α m and the variance of array distortion All are given an inverse gamma prior distribution, which is expressed as follows: Therefore:

[0030]

[0031] p(α m )=IG(α m |a2,b2)

[0032]

[0033] Where a1, a2, and a3 represent non-negative shape parameters, and b1, b2, and b3 represent non-negative scaling parameters.

[0034] A further technical solution of the present invention: when S4 The iteration stops when the change in the mean vector of the posterior distribution is less than a preset threshold.

[0035] A computer system is characterized by comprising: one or more processors, and 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 cause the one or more processors to implement the method described above.

[0036] A computer-readable storage medium is characterized by storing computer-executable instructions, which, when executed, are used to implement the above-described method.

[0037] The beneficial effects of this invention are as follows:

[0038] The present invention provides a Bayesian DOA estimation method for array-distorted passive synthetic aperture sonar, which has the following advantages compared with the prior art:

[0039] (1) Traditional passive synthetic aperture direction finding methods rely on the coherence of the received signal to estimate the phase correction factor, thus failing to effectively synthesize the aperture for single-target signals with random or even discontinuous phases. However, the non-cooperative nature of underwater reconnaissance and the time-varying characteristics of the underwater acoustic environment cause the phase of the array received signal to change rapidly in the time domain. In such signal environments, the direction finding performance of traditional PSAS algorithms such as time-domain ETAM, FFTSA, and Yen-Carey is severely degraded. The method designed in this invention uses the maximum likelihood criterion to synthesize the virtual array, instead of obtaining the virtual aperture by estimating the correlation factors of overlapping array elements to perform phase compensation on the continuous measurement output. At the same time, the variational EM method is used to solve the problem of high computational complexity of the maximum likelihood estimation method. Therefore, when the waveform of the incident signal exhibits random spatiotemporal characteristics, the method designed in this invention can use nonlinear spectral estimation to eliminate the adverse effects of phase fluctuations caused by overlapping correlation factors, and has high computational efficiency.

[0040] (2) To improve detection capabilities, the trend in towed array development is towards thinner and longer arrays. This trend exacerbates the problem of maintaining an ideal straight-line state during towing. Among existing array shape estimation algorithms, the element position estimation method based on hydrodynamic models is sensitive to prior system parameters, and determining these parameters is a technical challenge. Most algorithms based on hydrophone data are computationally intensive, have poor real-time performance, and lack practicality. Methods using non-acoustic sensors to measure array shape are highly dependent on high-performance hardware and are costly. Furthermore, the application requirements for large-aperture sonar arrays urgently require the use of small-aperture arrays to obtain the array gain and azimuth resolution of long arrays, thereby reducing engineering implementation costs. However, existing towed array shape estimation methods do not utilize small-scale moving arrays to convert time accumulation into spatial gain, thus failing to use synthetic aperture technology to solve the problem of insufficient aperture in passive azimuth estimation of towed array sonar targets. How to combine array aperture expansion technology and a small amount of prior model knowledge to invert and obtain the entire array structure is the main problem to be solved in towed array shape estimation technology. This invention addresses the sensitivity of traditional passive synthetic aperture direction finding algorithms to array distortion in flexible arrays. It proposes leveraging the advantages of Bayesian methods for joint estimation of multiple parameters, jointly modeling array distortion parameters and signal parameters, and employing a variational Bayesian expectation-maximization algorithm to iteratively update the model parameter values. By continuously approximating the assumed model with the observed data, the optimization of model parameters is guided, ultimately enabling simultaneous optimization estimation of array parameters and signal parameters. This reduces the computational load of the direction finding process and improves the accuracy of the direction finding results. Attached Figure Description

[0041] The accompanying drawings are for illustrative purposes only and are not intended to limit the invention. Throughout the drawings, the same reference numerals denote the same parts.

[0042] Figure 1 Wavefront diagram of acoustic signals received by a flexible, thin-wire towed array;

[0043] Figure 2 Spatial spectrum comparison diagram of the method designed in this invention with ETAM and FFTSA algorithms in a random signal environment;

[0044] Figure 3 The diagram shows the variation of the DOA estimation RMSE value of the method designed in this invention and the traditional passive synthetic aperture direction finding algorithm with SNR (the incident source is a random source with 2 sources).

[0045] Figure 4 The diagram shows the variation of the DOA estimation RMSE value of the method designed in this invention and the traditional passive synthetic aperture direction finding algorithm with the number of snapshots (the incident source is a random source, and the number of sources is 2).

[0046] Figure 5 This is a flowchart of the method of the present invention. Detailed Implementation

[0047] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention. Furthermore, the technical features involved in the various embodiments of this invention described below can be combined with each other as long as they do not conflict with each other.

[0048] This invention provides a Bayesian method for estimating the DoA (DoA) of a passive synthetic aperture sonar with array distortion. This method requires only a small amount of prior array information and can obtain the maximum likelihood estimates of the DoA and array distortion parameters by iteratively maximizing the marginal likelihood function, exhibiting statistical performance asymptotically equivalent to the method proposed by Nuttall. The specific design process is as follows: First, hierarchical probability modeling is performed on the array observation data. This model applies a binary prior distribution to the signal vector to include sparsity-induced characteristics. Second, a variational Bayesian method is used to iteratively maximize the lower bound of the marginal likelihood function of the array distortion parameters and latent variables. In this way, the estimated posterior distribution can be used to iteratively update each unknown parameter and assess their uncertainty. Furthermore, the method designed in this invention can also effectively utilize array motion information to synthesize a virtual aperture.

[0049] like Figure 5 As shown, the present invention includes the following steps:

[0050] Step 1: Use an N-element towed linear array, with an element spacing of d and the number of incident sound sources set to K. Because all elements are sealed within a flexible cable, the array cannot maintain an ideal straight-line shape during towing. A piecewise linear model is used to characterize the distorted array, that is, using the line connecting the towing point and its adjacent hydrophone as the reference, and the vector ρ = [ρ1,…,ρ...]. N-2 Describe the angles between the remaining hydrophone connections and the baseline. If the array's drag speed is v, then the received data z at the i-th sampling time... i It can be represented as:

[0051] z i =A i (θ,ρ)b i +ε i

[0052] in, Let a represent the array manifold matrix at the i-th sampling time. i (θ k ,ρ) represents the azimuth θ corresponding to the k-th incident source. k and the steering vector at the i-th sampling time, This represents the complex amplitude vector of the incident signal. and Let ε represent the amplitude and phase of the k-th source at the i-th sampling time, respectively. i =[ε1(t) i ),…,ε N (t i )] T Represents a Gaussian white noise vector;

[0053] Step 2: Define a spatially overcomplete array manifold matrix A i (υ,ρ)=[a i (υ0,ρ),…,a i (υ M-1 ,ρ)], where M represents the preset number of discrete grid points in the spatial domain, then z i It can be represented as the following sparse model:

[0054]

[0055] in, It is a sparse vector, constructed as follows: assuming a predefined spatial discrete grid of points υ m-1 If it is not included in the actual DOA set, then The m-th element The other positions correspond The vector element values ​​represent the true amplitude of the signal incident from that direction.

[0056] Step 3: Perform hierarchical probability modeling for each unknown variable in the sparse signal model, so that the sparse reconstruction problem is transformed into an l0 norm optimization problem;

[0057] Step 4: Iteratively update the posterior distribution or maximum likelihood estimate of each unknown variable using the variational EM algorithm. When The iteration stops when the change in the mean vector of the posterior distribution is less than a preset threshold.

[0058] Step 5: Observe the spatial grid points υ = [υ0,υ1,…υ M-1 [] as the x-axis, with Using the mean vector of the convergent solution as the ordinate, plot the spatial amplitude spectrum. From the amplitude spectrum, obtain the first K peaks in descending order of amplitude. The horizontal coordinate angle values ​​corresponding to these peaks are the incident signal DOA.

[0059] To enable those skilled in the art to better understand the present invention, the present invention will be described in detail below with reference to specific embodiments.

[0060] The implementation steps of this invention are as follows:

[0061] Step 1: Obtain the observation data z of the towed array at the i-th sampling time. i .

[0062] Assuming K sound sources are incident on an N-element linear hydrophone array with element spacing d, then the nth hydrophone ∈ {0,1,…,N-1} will be at position t. i The signal received at any given time can be represented by the following formula:

[0063]

[0064] Among them, t i =iΔt, i = 0, 1, ..., L-1 represents the i-th sampling time, and Δt represents the time sampling interval; s k (t i ) and ε n (t i The waveforms of the k-th signal source and the ambient noise are represented respectively. ω and ω0 represent the k-th source at sampling time t, respectively. i The amplitude, phase, and angular frequency; τ nk It is the time delay from the k-th signal source, with reference to the array head, to the n-th hydrophone.

[0065] A piecewise linear model is used to model the distorted array during towing, that is, taking the line connecting the towing point and its adjacent hydrophones as the reference, and using the vector ρ=[ρ1,…,ρ N-2 Describe the angles between the remaining hydrophone connection lines and the baseline, such as... Figure 1 As shown. If the position of the first element of the array at the starting point of the timing is ((N-1)d,0), then the Cartesian coordinates of the nth hydrophone at this moment are:

[0066]

[0067]

[0068] According to the array distortion model above, ρ0 = 0. At this point, the time delay of the nth hydrophone relative to the first array element is:

[0069] τ nk =-d nk / c

[0070] Where, d nk The nth hydrophone and the reference hydrophone are located at the kth incident source azimuth θ, calculated based on the above array distortion model. k The spatial distance is expressed as: c represents the speed of sound.

[0071] When the hydrophone array is dragged at a constant speed v, each hydrophone at time ti The spatial coordinates are:

[0072] x n (t i )=x n (0)+vt i ;y n (t i )=y n (0)

[0073] At this point, the waveform of the k-th signal source incident on the n-th hydrophone can be represented as:

[0074]

[0075] Where, ω k1 and Let represent the time angular frequency and the spatial angular frequency of the signal, respectively, and their expressions are given by the following two equations:

[0076]

[0077]

[0078] The array manifold matrix is ​​represented as Where the guiding vector a i (θ k The expression for the nth element of ,ρ) is After the signal modeling described above, the array observation data z at the i-th sampling time is... i It can be represented as:

[0079] z i =A i (θ,ρ)b i +ε i

[0080] in, Let ε be a vector composed of the complex amplitudes of each incident signal. i =[ε1(t) i ),…,ε N (t i )] T Let represent the received noise at the i-th sampling time, which follows a mean of 0 and a variance of . The Gaussian distribution.

[0081] Step 2: Construct the likelihood function p(Z|υ,ρ) for the array receiving data.

[0082] Define an overcomplete array manifold matrix A i (υ,ρ)=[a i (υ0,ρ),…,a i (υM-1 ,ρ)], where M represents the preset number of spatial discrete grid points, and z i This can be represented as a sparse model as follows:

[0083]

[0084] in, It is a sparse vector, constructed as follows: assuming a predefined spatial discrete grid of points υ m-1 If it is not included in the actual DOA set, then The remaining positions The vector element values ​​represent the true amplitude of the signal incident from that direction. After the above sparse modeling, the DOA estimation problem is equivalent to determining the DOA of a given signal. The position coordinates of the non-zero elements. Furthermore, the prior probability model of the array perturbation vector ρ is expressed as having a mean of... The covariance matrix is The Gaussian distribution, without loss of generality, assumes but

[0085]

[0086] in, L is a triangular matrix. Under this sparse signal model, the likelihood function of the array observation data can be expressed as:

[0087]

[0088] Step 3: Construct a hierarchical Bayesian probability model.

[0089] The statistical characteristics of each latent variable in the hierarchical probability model are modeled. Specifically, this involves modeling the signal vector... Applying a sparse prior distribution-binary distribution, the likelihood function can then be expressed as:

[0090]

[0091] Among them, z ij Indicate z i The j-th element; binary indicator variable γ m It follows the Bernoulli distribution as follows:

[0092]

[0093] The m-th element of the signal vector It follows the following pattern: mean 0, variance α. m The complex Gaussian distribution:

[0094]

[0095] Integrating the two prior distributions above, the resulting marginal distribution is the l0 norm.

[0096] Next, corresponding prior distributions are designed for the remaining latent variables to construct a fully Bayesian probability model. Specifically, the perturbation vector ρ is subjected to the following conditions: mean 0, variance... Gaussian priors:

[0097]

[0098] noise variance Signal power α m and the variance of array distortion All are given an inverse gamma prior distribution, which is expressed as follows: Therefore:

[0099]

[0100] p(α m )=IG(α m |a2,b2)

[0101]

[0102] Where a1, a2, and a3 represent non-negative shape parameters, and b1, b2, and b3 represent non-negative scaling parameters.

[0103] Step 4: Use the variational EM algorithm to obtain the posterior distribution of each latent variable in the hierarchical probability model and the maximum likelihood estimate of each model parameter.

[0104] 4a) The initialization methods for each latent variable and model parameter are as follows:

[0105]

[0106]

[0107]

[0108]

[0109] Vectors γ and ρ are initialized to vectors containing all 1s and vectors containing all 0s, respectively.

[0110] 4b) Update latent variables and model parameters, divided into E-step and M-step:

[0111] E-Step:

[0112] posterior distribution It follows a complex Gaussian distribution:

[0113]

[0114] in,

[0115]

[0116]

[0117] For the sake of convenience in formula representation, the mathematical expectation operation will be represented by <·>.

[0118] The posterior distribution q(ρ|Z;β) of ρ is a complex Gaussian distribution as follows:

[0119]

[0120] in,

[0121]

[0122]

[0123] γ m The posterior distribution q(γ) m |Z;β) is a complex Gaussian distribution as follows:

[0124]

[0125] in,

[0126]

[0127] M-step:

[0128] Lower bound of the likelihood function For model parameter vectors Taking the partial derivatives of each element in the equation and setting the results to 0, we obtain the maximum likelihood estimates for each parameter, as shown below:

[0129]

[0130]

[0131]

[0132]

[0133] 4c) Determine the updated Does it meet the convergence condition? Where ε represents the decision threshold, the value of which is determined based on the accuracy requirements in the actual application. If the convergence condition is met, the iteration stops; otherwise, return to step 4b) and continue iterating until the above convergence condition is met.

[0134] Step 5: According to The optimal estimation results are used to plot a spatial spectrum, and each DOA is determined based on its peak value.

[0135] The result obtained in step 4 The convergence result is a sparse vector, with most elements close to 0 and containing only K non-zero values. The spatial discrete grid points corresponding to these K non-zero values ​​are the DOA of the incident signal. Let the observation space grid points be υ=[υ0,υ1,…υ M-1 [] is the x-axis (unit: degrees), with The amplitude of the mean vector of the convergent solution is taken as the base-10 logarithm (in dB), and a spatial spectrum is plotted. From this spectrum, the first K peaks are obtained in descending order, and the angle values ​​of the horizontal axis corresponding to these peaks are the required incident signal DOA.

[0136] The effects of this invention can be illustrated by the following simulations:

[0137] 1. Simulation conditions:

[0138] The uniform linear hydrophone array has 16 elements, the narrowband incident sound source has a frequency of 100Hz, the element spacing is half the signal wavelength, and the drag speed is 12 knots. When the signal sampling frequency follows the Nyquist criterion, the array moves half a wavelength every 250 time sampling intervals.

[0139] 2. Simulation content and results:

[0140] Simulation 1: Figure 2 The image shows a spatial spectrum comparison between the method designed in this invention and the ETAM and FFTSA algorithms. Two equal-power random sound sources are incident on the array from directions θ1 = 20° and θ2 = 30°, respectively, with a signal-to-noise ratio of 0 dB. During the towing process, each hydrophone deviates from its ideal position, and the position error of the array elements follows a Gaussian distribution with a mean of 0 and a standard deviation of 3°. The number of sampling snapshots is set to 8000.

[0141] from Figure 2 Simulation results show that the ETAM and FFTSA algorithms can only estimate the location of one source, while the algorithm designed in this invention can accurately distinguish between two sources. Simulation results demonstrate the disadvantage of traditional passive synthetic aperture direction-finding algorithms based on phase correction factor estimation in distinguishing random sources.

[0142] Simulation 2: The direction-finding performance of the algorithm designed in this invention is compared with other traditional algorithms through 100 Monte Carlo trials. The evaluation criterion is the root mean square error (RMSE), which is defined as follows: in and Let $\mathbf{i}$ be the true azimuth and the estimated azimuth of source $k$ in the $i$-th Monte Carlo experiment, respectively. The spatial angle discrete interval is 1°, the incident azimuths of the two random sources are -15° and 10°, respectively, the number of sampling snapshots is set to 8000, and other simulation parameters are consistent with those in Simulation 1.

[0143] Figure 3 The figure shows the variation of the DOA estimation RMSE curves with SNR for each comparison algorithm. Additionally, the Cramer-Rao Lower Bound (CRLB) for the angle estimation results is plotted in the figure. Figure 3 The results show that the ETAM, METAM, TD-ETAM, Yen-Carey, and ETAM+SPICE algorithms (X. Shang, J. Li, and P. Stoica, “Weighted SPICE algorithms for range-doppler imaging using one-bit automotive radar,” IEEE J. Sel. Topics Signal Process., vol. 15, no. 4, pp. 1041-1054, 2021) all fail to achieve high DOA estimation accuracy. This is because the received signals on overlapping array elements in adjacent measurement samples lack coherence, making it impossible to synthesize a virtual aperture. The FFTSA and Jin-Li algorithms also exhibit low DOA estimation accuracy, indicating that these two algorithms are not robust to array element position errors. Figure 3 The results shown further confirm the superior direction-finding performance of the method designed in this invention, indicating that this method is an efficient and practical passive synthetic aperture direction-finding algorithm under array distortion and random source environments.

[0144] Simulation 3: With the signal-to-noise ratio of Simulation 2 fixed at 0dB and other simulation parameters unchanged, plot the RMSE curves of each comparison algorithm as a function of the number of snapshots.

[0145] Figure 4 The simulation results show that even with a sufficient number of snapshots, the ETAM, METAM, TD-ETAM, Yen-Carey, and ETAM+SPICE algorithms cannot achieve high-precision direction-finding results. This again demonstrates the limitations of traditional passive synthetic aperture direction-finding algorithms in estimating the phase correction factor. The simulation results also show that the direction-finding accuracy of the method designed in this invention is consistently higher than other comparative algorithms and approaches the optimal estimation result CRLB, demonstrating the superior direction-finding performance of this method in distorted array configurations and random source environments.

[0146] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any person skilled in the art can easily conceive of various equivalent modifications or substitutions within the scope of the technology disclosed in the present invention, and such modifications or substitutions should all be covered within the scope of protection of the present invention.

Claims

1. A Bayesian DOA estimation method for array-distorted passive synthetic aperture sonar, characterized in that, A towed linear array with N elements is used to receive K incident sound sources, including: S1: Based on the azimuth of the incident signal, the complex amplitude vector of the incident signal, and the angle between the line connecting the array elements and the baseline, the observation data of the towed array are obtained. S2: Define the spatial overcomplete array manifold matrix, and convert the observation data into a sparse signal model based on the spatial overcomplete array manifold matrix; S3: Perform hierarchical probability modeling for each latent variable in the sparse signal model, so that the sparse reconstruction problem is transformed into an l0 norm optimization problem; S4: The variational EM algorithm is used to iteratively update the posterior distribution of each latent variable and the maximum likelihood estimate of each model parameter in the hierarchical probability model; S5: Using the observed spatial grid points as the abscissa and the mean vector of the iterated sparse signal as the ordinate, plot the spatial amplitude spectrum. Obtain the first K peaks from the amplitude spectrum in descending order of amplitude. The abscissa angle value corresponding to the peak is the incident signal DOA.

2. The Bayesian DOA estimation method for array-distorted passive synthetic aperture sonar according to claim 1, characterized in that, The expression for the observed data in S1: z i =A i (θ,ρ)b i +e i in, Let a represent the array manifold matrix at the i-th sampling time. i (θ k ,ρ) represents the azimuth θ corresponding to the k-th incident source. k And the steering vector at the i-th sampling time, ρ=[ρ1,…,ρ N-2 [ ] is the angle between the line connecting each hydrophone segment and the baseline; This represents the complex amplitude vector of the incident signal. and Let ε represent the amplitude and phase of the k-th source at the i-th sampling time, respectively. i =[ε1(t) i ),…,ε N (t i )] T This represents a Gaussian white noise vector.

3. The Bayesian DOA estimation method for array-distorted passive synthetic aperture sonar according to claim 2, characterized in that, The expression for the sparse signal model described in S2 is: Among them, A i (υ,ρ)=[a i (υ0,ρ),…,a i (υ M-1 [,ρ)] is the overcomplete array manifold matrix, a i (υ m ,ρ) represents the angle υ at the i-th sampling time. m The array guide vector, υ=[υ0,υ1,…υ M-1 [] represents the observation space grid points, and M represents the preset number of spatial discrete grid points; It is a sparse vector, constructed as follows: assuming a predefined spatial discrete grid of points υ m-1 If it is not included in the actual DOA set, then The m-th element The other positions correspond to The vector element values ​​represent the true amplitude of the signal incident from that direction.

4. The Bayesian DOA estimation method for array-distorted passive synthetic aperture sonar according to claim 3, characterized in that, S3 specifically refers to: For signal vectors Applying a sparse prior distribution-binary distribution, the likelihood function can then be expressed as: Among them, z ij Indicate z i The j-th element; binary indicator variable γ m It follows the Bernoulli distribution as follows: The m-th element of the signal vector It follows the following: mean 0, variance α m The complex Gaussian distribution: Integrating the two prior distributions above, the resulting marginal distribution is the l0 norm. Apply the following conditions to the disturbance vector ρ: mean 0, variance ? Gaussian priors: noise variance Signal power α m and the variance of array distortion All are given an inverse gamma prior distribution, which is expressed as follows: Therefore: Where a1, a2, and a3 represent non-negative shape parameters, and b1, b2, and b3 represent non-negative scaling parameters.

5. The Bayesian DOA estimation method for array-distorted passive synthetic aperture sonar according to claim 4, characterized in that, S4 was fooled The iteration stops when the change in the mean vector of the posterior distribution is less than a preset threshold.

6. A computer system, characterized in that... include: 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 cause the one or more processors to perform the method of any one of claims 1-5.

7. A computer-readable storage medium, characterized in that... The device stores computer-executable instructions, which, when executed, are used to implement the method described in any one of claims 1-5.

Images