A fast soft threshold iterative deconvolution beamforming method in complex domain based on probability mapping

Abstract

The present invention discloses a complex domain fast soft threshold iterative deconvolution beamforming method based on probability mapping, which takes the conventional beamforming result of the echo signal of the sonar array as the initial value of the first iteration; performs complex gradient descent on the initial value of the iteration, and uses fast Fourier transform to accelerate the matrix multiplication calculation in the gradient descent to obtain the intermediate result of the gradient descent; establishes a probability mapping model, adds beam clustering prior, and uses multi-task Bayesian compressed sensing to quickly solve the beam result and the corresponding Gaussian distribution parameters; performs momentum update according to the beam result to obtain the initial value of the iteration of the next iteration, and accelerates the iterative convergence. The deconvolution beamforming is extended to the complex domain, which makes full use of the beam phase information, is more suitable for practical applications, effectively reduces the main lobe width and side lobe intensity, and has a high noise resistance, thereby improving the array imaging quality.

Description

Technical Field

[0001] The invention relates to the technical field of array narrowband signal beamforming, and in particular to a complex domain fast soft threshold iterative deconvolution beamforming method based on probability mapping. Background Art

[0002] Array beamforming is widely used in acoustic imaging systems such as sonar and radar. The observation image is obtained by transmitting sound waves to penetrate the observation scene and receiving echoes to form beams in all directions. For example, patent document CN109283536A discloses a multi-beam bathymetric sonar water imaging beamforming algorithm, and patent document CN114200457A discloses a beamforming method. Traditional beamforming methods have a wide main lobe width, weak sidelobe suppression capability, and low beam resolution.

[0003] Deconvolution beamforming methods are used to improve beam quality. Existing deconvolution methods are mostly based on beam intensity values, that is, real domain beam results are deconvolved, assuming that the target is incoherent and ignoring phase information. However, for high-resolution array imaging, narrowband signals are mostly used, and the target is coherent. Real domain methods have model mismatch and poor results. Summary of the invention

[0004] In view of the above, an object of the present invention is to provide a complex domain fast soft threshold iterative deconvolution beamforming method based on probability mapping to effectively improve the beam resolution.

[0005] To achieve the above-mentioned object of the invention, the present invention provides a complex domain fast soft threshold iterative deconvolution beamforming method based on probability mapping, comprising the following steps:

[0006] Step 1, using the conventional beamforming result of the echo signal of the sonar array as the initial value of the first iteration;

[0007] Step 2: Perform complex gradient descent on the initial value of the iteration, and use fast Fourier transform to accelerate the matrix multiplication calculation in the gradient descent to obtain the intermediate result of the gradient descent;

[0008] Step 3: Use the minimum mean square error estimation to perform probability mapping between the expected beam result and the intermediate result of gradient descent, and abstract the expected beam result into a Gaussian distribution parameter representation;

[0009] Step 4: Based on the expected beam result represented by the Gaussian distribution parameters, add the beam clustering prior to reconstruct the clustered expected beam and the corresponding Gaussian distribution parameters;

[0010] Step 5, after quickly solving the clustered expected beam and its Gaussian distribution parameters according to the RVM framework in the multi-task Bayesian compressed sensing method, the expected beam result of the current iteration of the expected beam direction in the observation domain is calculated;

[0011] Step 6: Update the momentum according to the expected beam result of the current iteration to obtain the initial value of the next iteration;

[0012] Step 7, repeat iterative steps 2 to 6 until the iteration is terminated to obtain the final expected beam result.

[0013] For a two-dimensional sonar planar array with M array elements of size R×C, the observation range contains T targets, and the observation domain is divided into P×Q, with a total of N expected beam directions. For the nth expected beam direction The corresponding conventional beamforming result It is expressed as:

[0014]

[0015] Where, subscript t is the target number, x t Indicates the echo signal strength of the target number, λ is the signal wavelength, represents the direction angle and pitch angle of the target incident, sinc() represents the non-normalized Single function, u nt = sin(θ n )-sin(θ t ), PSF() is the point spread function of deconvolution. For all N expected beam directions within the observation range, the beamforming formula can be written in matrix form:

[0016] B=Φx

[0017] Among them, B represents the beamforming result of N beams, x represents the position vector of T targets, Φ is the spatial mapping from target to beam, that is, the PSF matrix, whose n rows and t columns are

[0018] Preferably, in step 2, the initial value of the iteration is rewritten into a matrix form, and based on the Wirtinger derivative formula, the gradient of the complex PSF matrix is ​​deduced as:

[0019]

[0020] Among them, y (k) is the initial value of the kth iteration, the deconvolution objective function The calculation formula for the intermediate result of gradient descent is as follows:

[0021] r (k) =y (k) +μ▽f(y (k) )=y (k) +μΦ H (Φy (k) -B)

[0022] Among them, r (k) is the intermediate result of the gradient descent of the kth iteration, μ is the gradient step size, and its value is the reciprocal of the Lipschitz constant L of the matrix Φ in the point spread function PSF. Lipschitz constant L=max[diag(Φ H Φ)], max() is the function for taking the maximum eigenvalue, and diag() is the function for constructing a diagonal matrix;

[0023] Because the array element spacing in two directions in the two-dimensional planar array is equal, it is proved that Φ is a Hermitian matrix with complex conjugate symmetry. The point multiplication after fast Fourier transform is used to accelerate the matrix multiplication operation. The intermediate result of the gradient descent accelerated by fast Fourier transform is calculated as:

[0024]

[0025] in, represents the fast Fourier transform, represents inverse fast Fourier transform.

[0026] Preferably, in step 3, the expected beam result x (k) And the intermediate result r of gradient descent (k) Make a probability mapping, expressed as:

[0027]

[0028] in, Represents a probability mapping, and the specific operations are as follows:

[0029] The minimum mean square error estimate for the kth iteration is According to the Bayesian probability formula, when the probability distributions of r and x are known a priori, to obtain the posterior probability P(x|r), we only need to determine the conditional probability P(r|x). Assume r (k) =Ex (k) +ε, E is the identity matrix, ε is the variance σ 2 The conditional probability is as follows:

[0030]

[0031] Considering that the target is sparse in the entire observation range, it is assumed that x (k) It also conforms to the zero-mean Gaussian distribution with a variance of α -1 To further ensure sparsity, add Gamma distribution constraints to α -1 , x (k) The abstract Gaussian probability distribution is expressed as follows:

[0032]

[0033] Where i represents the i-th beam direction, It means that 0 is the mean. is a Gaussian distribution with variance, Γ(α i |a,b) represents the Gamma distribution with a and b as parameters, where a and b are the Gamma distribution hyperparameters added to α;

[0034] For complex probabilities, we try to make the real and imaginary parts into two independent parts that share priors and parameters. (k) Each element in The probability of is calculated as follows:

[0035]

[0036] in, and represent the real and imaginary parts respectively.

[0037] Preferably, in step 4, clustering priors are added to the four adjacent beams in the horizontal and vertical directions to share Gaussian distribution parameters. After the clustering priors are added, the clustering expected beams are expressed as:

[0038]

[0039] in, and represents the clustering expected beam and the corresponding Gaussian distribution parameters, Represents a complex number domain.

[0040] Preferably, in step 5, the Gaussian distribution parameter α of the expected beam is separated and clustered according to the RVM framework in the multi-task Bayesian compressed sensing method. s The maximum probability point estimate, the Gaussian distribution parameter α that makes the probability zero s as follows:

[0041]

[0042] in, C=I4+I s A -1 I s , r H represents the conjugate transpose of r, d represents the array element spacing, I si Indicates I s The i-th column vector of , I4 is a 4x4 matrix of all 1s, A = diag(α s ), represents the removal of the i-th column vector after the inversion of C, a is the aforementioned Gamma distribution hyperparameter, I s E is obtained by adding clustering constraints, and the intermediate process matrix I is s ′ and E′ are Is and E vector is reshaped from N×1 to a matrix of size P×Q, then I s The elements of ′ are (p=1,...,P-1; q=1,...,Q-1):

[0043]

[0044] I s By matrix I′ s Compressed into a vector, we get;

[0045] The Gaussian distribution parameter α s The Gaussian distribution parameter α of the original beam before clustering is restored, and the calculation formula is:

[0046]

[0047] After obtaining the Gaussian distribution parameter α of the original beam, the expected beam result of the current iteration of the expected beam direction in the observation domain is calculated according to the Gaussian distribution parameter α.

[0048] Preferably, in step 6, momentum updating is performed according to the expected beam result of the current iteration, including:

[0049]

[0050]

[0051] Among them, y (k+1) is the initial value of the next iteration, t (k+1) Update the scale parameter for the momentum.

[0052] Preferably, in step 7, y is obtained (k+1) After that, the initial value y of the current round of iteration (k) After the modulus is determined, the difference is calculated and the absolute value of the difference is calculated. If the absolute value is less than the iteration termination threshold, the iteration is exited and the final beam result x is returned. (k) , otherwise continue iterating.

[0053] Compared with the prior art, the present invention has the following beneficial technical effects:

[0054] The complex domain fast soft threshold iterative deconvolution beamforming method based on probability mapping provided by the present invention extends conventional deconvolution to the complex domain, uses complex gradient descent and probability mapping with clustering prior to make it closer to the actual beamforming model, improves the beam resolution accuracy and anti-noise robustness; uses fast Fourier transform (FFT) to accelerate large-scale matrix multiplication calculations, effectively reducing computational complexity; and selects a better iterative initial value through momentum update to accelerate iterative convergence, which helps to reduce computational delay. Therefore, the present invention is suitable for narrowband deconvolution wave formation of arrays, and improves array imaging quality. BRIEF DESCRIPTION OF THE DRAWINGS

[0055] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0056] Figure 1 A schematic flow chart of a complex domain fast soft threshold iterative deconvolution beamforming method based on probability mapping provided by an embodiment;

[0057] Figure 2 A schematic diagram of a planar array far-field signal model provided in an embodiment;

[0058] Figure 3 A schematic diagram of matrix transformation after adding clustering priors provided in an embodiment;

[0059] Figure 4 A comparison chart of the noise resistance of the deconvolution method provided in the embodiment and the conventional intensity-based deconvolution method;

[0060] Figure 5 The conventional beamforming test results of the planar array provided by the embodiment;

[0061] Figure 6 A beam result ratio diagram of the deconvolution method provided in the embodiment and the conventional intensity-based deconvolution method. DETAILED DESCRIPTION

[0062] To make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific implementation methods described herein are only used to explain the present invention and do not limit the scope of protection of the present invention.

[0063] In order to solve the problem of poor performance of traditional deconvolution methods, the embodiment provides a complex domain fast soft threshold iterative deconvolution beamforming method based on probability mapping on the basis of the conventional intensity-based deconvolution method, so as to realize narrowband beamforming high-resolution deconvolution and improve the resolution of the beam.

[0064] like Figure 1 As shown, the complex domain fast soft threshold iterative deconvolution beamforming method based on probability mapping provided by the embodiment includes the following steps:

[0065] Step 1: Using the conventional beamforming result of the echo signal of the sonar array as the initial value of the first iteration.

[0066] The deconvolution beamforming method provided in the embodiment uses the conventional beamforming (CBF) result as the initial value, from which the beam is accurately estimated. Figure 2 As shown, for a two-dimensional sonar planar array containing M array elements of size R×C, the two-dimensional sonar planar array is located in the XOY plane of the three-dimensional space, the sides of the two-dimensional rectangular array are parallel to the X and Y coordinate axes respectively, the center of the array is located at the origin of the coordinates, and the coordinates of the (r, c)th transducer array element are expressed as (x r ,y c ), the array element spacing is d x ,d y , the beam focusing range is the far field area, and the observation range is determined by two directional parameters θ and Indicates that they are parallel to the X and Y axes respectively, and the value range is (-90°, 90°). The observation range contains T targets, and the observation domain is divided into P×Q, a total of N expected beam directions. For the nth expected beam direction The corresponding conventional beamforming result It is expressed as:

[0067]

[0068] Where, subscript t is the target number, x t Indicates the echo signal strength of the target number, λ is the signal wavelength, represents the direction angle and pitch angle of the target incident, sinc() represents the non-normalized Single function, u nt = sin(θ n )-sin(θ t ), PSF() is the point spread function of deconvolution. For all N expected beam directions within the observation range, the beamforming formula can be written in matrix form:

[0069] B=Φx

[0070] Among them, B represents the beamforming result of N beams, x represents the position vector of T targets, Φ is the spatial mapping from target to beam, that is, the PSF matrix, whose n rows and t columns are

[0071] Step 2: Perform complex gradient descent on the initial value of the iteration, and use fast Fourier transform to accelerate the matrix multiplication calculation in the gradient descent to obtain the intermediate result of the gradient descent.

[0072] In the embodiment, the deconvolution calculation is performed based on the complex domain of the beam, and the beam calculation is rewritten into a matrix form: B = Φx + ξ, where and noise The initial value of the first iteration is the conventional beamforming result in step 1, and the initial value of subsequent iterations is the momentum update result of the previous iteration. The initial value of the iteration is subjected to complex domain gradient descent, and the intermediate result of the gradient descent is obtained as the input of the probability mapping.

[0073] In the embodiment, the gradient of the complex PSF matrix is ​​derived based on the Wirtinger derivative formula:

[0074]

[0075] Among them, y (k) is the initial value of the kth iteration, the deconvolution objective function The calculation formula for the intermediate result of gradient descent is as follows:

[0076]

[0077] Among them, r (k) is the intermediate result of the gradient descent of the kth iteration, μ is the gradient step size, and its value is the reciprocal of the Lipschitz constant L of the matrix Φ in the point spread function PSF. Lipschitz constant L=max[diag(Φ H Φ)], max() is the function for taking the maximum eigenvalue, and diag() is the function for constructing a diagonal matrix;

[0078] Because the spacing between the array elements in two directions in a two-dimensional planar array is equal, that is, d x =d y = d, it is proved that Φ is a Hermitian matrix with complex conjugate symmetry, so it can be compressed into The point multiplication after fast Fourier transform is used to accelerate the matrix multiplication operation. The intermediate result of the gradient descent accelerated by fast Fourier transform is calculated as:

[0079]

[0080] in, represents the fast Fourier transform, represents inverse fast Fourier transform.

[0081] Step 3: Use the minimum mean square error estimation (MMSE) to perform probability mapping between the expected beam result and the intermediate result of gradient descent, and abstract the expected beam result into a Gaussian distribution parameter representation.

[0082] In the embodiment, the expected beam result x (k) And the intermediate result r of gradient descent (k) Make a probability mapping, expressed as:

[0083]

[0084] in, Represents a probability mapping, and the specific operations are as follows:

[0085] The minimum mean square error estimate for the kth iteration is According to the Bayesian probability formula, when the probability distributions of r and x are known a priori, to obtain the posterior probability P(x|r), we only need to determine the conditional probability P(r|x). Assume r (k) =Ex (k) +ε, E is the identity matrix, ε is the variance σ 2 The conditional probability is as follows:

[0086]

[0087] Considering that the target is sparse in the entire observation range, it is assumed that x (k) It also conforms to the zero-mean Gaussian distribution with a variance of α -1 To further ensure sparsity, add Gamma distribution constraints to α -1 , x (k) The abstract Gaussian probability distribution is expressed as follows:

[0088]

[0089] Where i represents the i-th beam direction, It means that 0 is the mean. is a Gaussian distribution with variance, Γ(α i |a,b) represents the Gamma distribution with a and b as parameters. a and b are the Gamma distribution hyperparameters added to α. The default values ​​are both 0 and are assigned according to the data.

[0090] For complex probabilities, we try to make the real and imaginary parts into two independent parts that share priors and parameters. (k) Each element in The probability of is calculated as follows:

[0091]

[0092] in, and represent the real and imaginary parts respectively.

[0093] Step 4: Based on the expected beam result represented by the Gaussian distribution parameters, add the beam clustering prior to reconstruct the clustered expected beam and the corresponding Gaussian distribution parameters.

[0094] Acoustic target imaging results are usually sparse in the entire observation range, but clustered near the target. Based on this feature, clustering priors are added to the beams to further make the probability model fit the actual situation. In the embodiment, clustering priors are added to the four adjacent beams in the horizontal and vertical directions, and the Gaussian distribution parameters are shared. After the clustering priors are added, the beam and parameter matrix change process is as follows: Figure 3 As shown, the clustered expected beam is expressed as:

[0095]

[0096] in, and represents the clustering expected beam and the corresponding Gaussian distribution parameters, Represents a complex number domain.

[0097] Step 5, after quickly solving the clustered expected beam and its Gaussian distribution parameters according to the RVM framework in the multi-task Bayesian compressed sensing method, the expected beam result of the current iteration of the expected beam direction in the observation domain is calculated.

[0098] In the embodiment, the Gaussian distribution parameter α of the expected beam is separated and clustered according to the RVM framework in the multi-task Bayesian compressed sensing method. s The maximum probability point estimate, the Gaussian distribution parameter α that makes the probability zero s as follows:

[0099]

[0100] in, C=I4+I s A -1 I s , r H represents the conjugate transpose of r, d represents the array element spacing, I si Indicates I s The i-th column vector of , I4 is a 4x4 matrix of all 1s, A = diag(α s ), represents the inversion of C after removing the i-th column vector, and a is the aforementioned Gamma distribution hyperparameter. s It is obtained by adding clustering constraints to E. For ease of understanding, let the intermediate process matrix I s ′ and E′ are I s and E vector is reshaped from N×1 to a matrix of size P×Q, then I s The elements of ′ are (p=1,...,P-1; q=1,...,Q-1):

[0101]

[0102] I s By matrix I′ sCompressed into a vector. By Gaussian distribution parameter α s The Gaussian distribution parameter α of the original beam before clustering is restored, and the calculation formula is:

[0103]

[0104] After obtaining the Gaussian distribution parameter α of the original beam, the expected beam result of the current iteration of the expected beam direction in the observation domain is calculated according to the Gaussian distribution parameter α. Specifically, the Gaussian distribution parameter α is substituted into the formula in step 3 to calculate the expected beam result x (k) .

[0105] Step 6: Update the momentum according to the expected beam result of the current iteration to obtain the initial value of the next iteration.

[0106] In the embodiment, momentum updating is performed according to the expected beam result of the current iteration, including:

[0107]

[0108]

[0109] Among them, y (k+1) is the initial value of the next iteration, t (k+1) Update the scale parameter for the momentum.

[0110] Step 7, repeat iterative steps 2 to 6 until the iteration is terminated to obtain the final expected beam result.

[0111] In the embodiment, when obtaining y (k+1) After that, the initial value y of the current round of iteration (k) After the modulus is determined, the difference is calculated and the absolute value of the difference is calculated. If the absolute value is less than the iteration termination threshold, the iteration is exited and the final beam result x is returned. (k) , otherwise continue to iterate and calculate the prior for the next iteration by doing momentum update to accelerate the iterative convergence.

[0112] In order to demonstrate the advantages of this method in imaging quality and noise immunity, the embodiment uses a 48×48 sonar array element, a planar array with half-wavelength equal spacing for comparison with other deconvolution methods, with an observation range of 60°×60° and an expected beam direction of 128×128. The comparison methods include FFT-NNLS, DAMAS2, CLEAN, RL, and FISTA. The complex domain fast soft threshold iterative deconvolution beamforming method is abbreviated as CFISTA.

[0113] In the noise experiment, the noise variance was gradually increased, and 100 random experiments were performed for each method. The normalized root mean square error (RMSE) between the deconvolution beam result under the current variance and the ideal situation was recorded. The experimental results are shown in Figure 4As can be seen from the figure, the RMSE of this method is significantly smaller than that of the comparison method under various noise intensities, and it has strong anti-noise ability.

[0114] In the beam imaging quality experiment, 8 target points were selected as the target direction in the expected observation direction, and the theoretical beam intensity was set to 1. The CBF beam results are shown in Figure 5 As can be seen from the figure, the main lobe of conventional beamforming is wide, and the beam directivity is poor after the target is gathered at the center. Using the CBF result as the initial value, deconvolution calculation is performed, and the results are as follows: Figure 6 As shown in the figure, a~f are the deconvolution results of FFT-NNLS, DAMAS2, CLEAN, RL, FISTA and this method respectively. Figure 6 It can be seen that several comparison methods based on intensity values ​​have more side lobes, and some main lobes are wider or have relative intensity deviations. This method completely suppresses the side lobes (the side lobe intensity is less than the cutoff intensity -40dB), and the beam result has only 8 expected directions of the beam main lobe, which is basically consistent with the ideal situation. The side lobe suppression and RMSE indicators are shown in Table 1 below.

[0115] Table 1

[0116] method Sidelobe Suppression RMSE CBF 11.93dB 2.3729 FFT-NNLS 23.82dB 0.0432 DAMAS2 17.24dB 0.1206 CLEAN 22.26dB 0.0618 RL 27.76dB 0.2379 FISTA 35.76dB 0.0276 This method ≥40dB <![CDATA[3.16×10 -6 ]]>

[0117] It can be seen from the table that the main indicators of this method are better than those of the comparison method, and the beam imaging quality is significantly improved.

[0118] The specific implementation methods described above provide a detailed description of the technical solutions and beneficial effects of the present invention. It should be understood that the above is only the most preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, supplements and equivalent substitutions made within the scope of the principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A complex domain fast soft threshold iterative deconvolution beamforming method based on probability mapping, characterized in that: The following steps are involved: Step 1, taking the conventional beamforming result of the echo signal of the sonar array as the initial value of the first iteration, including: for a two-dimensional sonar planar array including M array elements of size R×C, T targets are included in the observation range, and the observation domain is divided into P×Q total N expected beam directions, for the nth expected beam direction The corresponding conventional beamforming result It is expressed as: Where, subscript t is the target number, x t represents the echo signal strength of the target number, λ is the signal wavelength, and the array element spacing is d x ,d y , represents the direction angle and pitch angle of the target incident, sinc() represents the non-normalized Single function, u nt = sin(θ n )-sin(θ t ), PSF() is the point spread function of deconvolution. For all N expected beam directions within the observation range, the beamforming formula can be written in matrix form: B=Φx Among them, B represents the beamforming result of N beams, x represents the position vector of T targets, Φ is the spatial mapping from target to beam, that is, the PSF matrix, whose n rows and t columns are Step 2: Perform complex gradient descent on the initial value of the iteration, and use fast Fourier transform to accelerate the matrix multiplication calculation in the gradient descent to obtain the intermediate result of the gradient descent, including: rewriting the initial value of the iteration into a matrix form, and based on the Wirtinger derivative formula, deriving the gradient of the complex PSF matrix as: Among them, y (k) is the initial value of the kth iteration, the deconvolution objective function The calculation formula for the intermediate result of gradient descent is as follows: Among them, r (k) is the intermediate result of the gradient descent of the kth iteration, μ is the gradient step size, and its value is the reciprocal of the Lipschitz constant L of the matrix Φ in the point spread function PSF. Lipschitz constant L=max[diag(Φ H Φ)], max() is the function for taking the maximum eigenvalue, and diag() is the function for constructing a diagonal matrix; Because the array element spacing in two directions in the two-dimensional planar array is equal, it is proved that Φ is a Hermitian matrix with complex conjugate symmetry. The point multiplication after fast Fourier transform is used to accelerate the matrix multiplication operation. The intermediate result of the gradient descent accelerated by fast Fourier transform is calculated as: in, represents the fast Fourier transform, represents inverse fast Fourier transform; Step 3: Use the minimum mean square error estimation to perform probability mapping between the expected beam result and the intermediate result of gradient descent, and abstract the expected beam result into a Gaussian distribution parameter representation; Step 4: Based on the expected beam result represented by the Gaussian distribution parameters, add the beam clustering prior to reconstruct the clustered expected beam and the corresponding Gaussian distribution parameters; Step 5, after quickly solving the clustered expected beam and its Gaussian distribution parameters according to the RVM framework in the multi-task Bayesian compressed sensing method, the expected beam result of the current iteration of the expected beam direction in the observation domain is calculated; Step 6: Update the momentum according to the expected beam result of the current iteration to obtain the initial value of the next iteration; Step 7, repeat iterative steps 2 to 6 until the iteration is terminated to obtain the final expected beam result.

2. The complex domain fast soft threshold iterative deconvolution beamforming method based on probability mapping according to claim 1, characterized in that: In step 3, the expected beam result x (k) And the intermediate result r of gradient descent (k) Make a probability mapping, expressed as: in, Represents a probability mapping, and the specific operations are as follows: The minimum mean square error estimate for the kth iteration is According to the Bayesian probability formula, when the probability distribution of r and x is known a priori, to obtain the posterior probability P(x|r), we only need to determine the conditional probability P(r|x). Assuming r (k) =Ex (k) +ε, E is the identity matrix, ε is the variance σ 2 The conditional probability is as follows: Considering that the target is sparse in the entire observation range, it is assumed that x (k) It also conforms to the zero-mean Gaussian distribution with a variance of α -1 To further ensure sparsity, add Gamma distribution constraints to α -1 , x (k) The abstract Gaussian probability distribution is expressed as follows: Where i represents the i-th beam direction, It means that 0 is the mean. is a Gaussian distribution with variance, Γ(α i |a,b) represents the Gamma distribution with a and b as parameters, where a and b are the Gamma distribution hyperparameters added to α; For complex probabilities, we try to make the real and imaginary parts into two independent parts that share priors and parameters. (k) Each element in The probability of is calculated as follows: in, and represent the real and imaginary parts respectively.

3. The complex domain fast soft threshold iterative deconvolution beamforming method based on probability mapping according to claim 1, characterized in that: In step 4, clustering priors are added to the four adjacent beams in the horizontal and vertical directions, and the Gaussian distribution parameters are shared. After the clustering priors are added, the clustering expected beam is expressed as: in, and represents the clustering expected beam and the corresponding Gaussian distribution parameters, represents the complex field, and i represents the i-th beam direction.

4. The complex domain fast soft threshold iterative deconvolution beamforming method based on probability mapping according to claim 2, characterized in that: In step 5, the Gaussian distribution parameter α of the expected beam is separated and clustered according to the RVM framework in the multi-task Bayesian compressed sensing method. s The maximum probability point estimate, the Gaussian distribution parameter α that makes the probability zero s as follows: in, C=I4+I s A -1 I s , r H represents the conjugate transpose of r, d represents the array element spacing, I si Indicates I s The i-th column vector of , I4 is a 4x4 matrix of all 1s, A = diag(α s ), represents the removal of the i-th column vector after the inversion of C, a is the aforementioned Gamma distribution hyperparameter, I s E is obtained by adding clustering constraints, and the intermediate process matrix I is s ' and E' are I s and E vector is reshaped from N×1 to a matrix of size P×Q, then I s The elements of ' are (p=1,...,P-1;q=1,...,Q-1): I s By matrix I s 'Compress it into a vector to get; The Gaussian distribution parameter α s The Gaussian distribution parameter α of the original beam before clustering is restored, and the calculation formula is: After obtaining the Gaussian distribution parameter α of the original beam, the expected beam result of the current iteration of the expected beam direction in the observation domain is calculated according to the Gaussian distribution parameter α.

5. The complex domain fast soft threshold iterative deconvolution beamforming method based on probability mapping according to claim 2, characterized in that: In step 6, momentum is updated according to the expected beam result of the current iteration, including: Among them, y (k+1) is the initial value of the next iteration, t (k+1) Update the scale parameter for the momentum.

6. The complex domain fast soft threshold iterative deconvolution beamforming method based on probability mapping according to claim 5, characterized in that: In step 7, we get y (k+1) After that, the initial value y of the current round of iteration (k) After the modulo, make the difference and calculate the absolute value of the difference. If the absolute value is less than the iteration termination threshold, the iteration is exited and the final beam result x is returned. (k) , otherwise continue iterating.

Images