step three: model parameter broadening treatment: and (3) taking equal weight to the regular term and the minimum module term in the formula (10), neglecting the constant term, and writing into a suboptimal regular cost function form:

it can be seen that the optimization solution process of equation (16) does not depend on the model parameters of Cauchy distribution, thereby improving the robustness of the algorithm. The Q matrix of equation (12) becomes:

\tilde{Q} = diag {{YY}^{H}} - - - (17)

while preserving the lambda parameter in equation (11) to satisfy the positive nature of the matrix.

Step four: iterative solution of high-resolution spatial spectra:

the iterative solution steps are as follows:

2) The ith iteration of

Y⁽ⁱ⁾＝Q^(i-1)F^H(λI_M+FQ^(i-1)F^H)^-1X(18)

The Q matrix is updated using the equation:

Q⁽ⁱ⁾＝diag{Y⁽ⁱ⁾Y^(i)H}；

3) judging a convergence condition:

|J⁽ⁱ⁾-J^(i-1)|/|J^(i-1)|＜ε(19)

ε is a small amount. If the condition is not met, making i equal to i +1, and returning to the step 2); if the condition is met, outputting the space spectrum Y of the current iteration cycle⁽ⁱ⁾。

The method comprises the following specific steps:

the method comprises the following steps: considering an M-element equidistant linear array with an array element spacing of d, the frequency domain signal vector of the k-th frequency point can be expressed as:

X(k)＝[X₀(k)，...，X_M-1(k)]^T(7)

wherein,

the signal of the mth array element at the kth frequency point, M is 0, 1, …, M-1, wm (k) is the noise signal of the corresponding kth frequency point, c is the sound velocity, D is the signal source number, and θ i is the azimuth of the ith signal relative to the positive horizontal direction.

The nth beam output signal of the kth frequency point can be expressed as

The above formula is a space-domain zero-filling DFT form, N is the number of space-domain DFT points, and the nth discrete space direction satisfies

θ_n＝asin{cn/(Nf_kd)}(10)

Equations (3) and (4) show that frequency domain beamforming can be rapidly achieved by zero-padding spatial DFT. However, similar to the time-domain to frequency-domain transform, the azimuthal resolution of conventional spectral estimation cannot be improved by padding zeros in the spatial dimension, which depends on the effective physical aperture size of the array.

According to equation (3), the inverse transform (IDFT) is

Expressed by a system of linear equations as

X(k)＝FY(k)(12)

Wherein, F is a Fourier transform matrix, and the m and n elements are as follows: [ f ] of]_m，n(1/M) exp (j2 pi mn/N). Equation (6), the underdetermined linear equation set problem, can be solved by various linear inversion calculation methods. I.e. without using conventional supplementsAnd zero processing, namely inverting the sparse power distribution of the signals in the spatial direction by using limited physical array element data information. By adopting proper constraint optimization processing, the azimuth resolution can be improved, and the relative strength of signals can be reflected.

Step two: based on a Cauchy-Gaussian signal model, solving a constraint optimization regular solution of space domain DFT linear inversion:

by utilizing the sparse distribution characteristics of the space incident plane wave signals in the azimuth set, the space domain DFT inversion can be regularized, the resolution of a space azimuth spectrum can be improved, and the space complex amplitude spectrum is approximately regarded as obeying Cauchy distribution:

<math><mrow><mi>p</mi><mo>{</mo><mi>Y</mi><mrow><mo>(</mo><mi>n</mi><mo>,</mo><mi>k</mi><mo>)</mo></mrow><mo>}</mo><mo>&Proportional;</mo><mfrac><mn>1</mn><mrow><mn>1</mn><mo>+</mo><msup><mrow><mo>|</mo><mi>Y</mi><mrow><mo>(</mo><mi>n</mi><mo>,</mo><mi>k</mi><mo>)</mo></mrow><mo>|</mo></mrow><mn>2</mn></msup><mo>/</mo><msubsup><mi>σ</mi><mi>Y</mi><mn>2</mn></msubsup></mrow></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>13</mn><mo>)</mo></mrow></mrow></math>

accordingly, the distribution of the vector Y in which the elements are independent of each other is p { Y (k) } ═ ii_np{Y(n，k)}。

Representing the variance of the variable Y.

Assuming that the noise W obeys a complex Gaussian distributionThe conditional distribution of the matrix frequency domain signals is as follows:

according to the Bayes (Bayes) criterion, the posterior distribution of the vector Y can be expressed as:

using equations (7) and (8), the maximum a posteriori probability solution (MAP) with equation (9) maximized, i.e., the cost function of the following equation, is minimized:

the first term is a regular term constructed according to Cauchy distribution of angular spectrum amplitudes and serves as a measure for describing angular spectrum sparsity. According to

Adjusting the sparsity of DFT inversion; the second term is the minimum mode error, i.e. the constraint term of the array signal itself, by noise power

The weights between the regularization term and the min-modulo term are adjusted.

The derivation of equation (10) and making it zero gives the optimal solution of the cost function as

Y＝(λQ^-1+F^HF)^-1F^HX(17)

Wherein,

q is the diagonal matrix:

equation (11) is formally a damped minimum modulus solution. Since Q is related to Y, an estimator satisfying equation (11) must be obtained by iterative solutionDue to the fact that

F^H(λI_M+FQF^H)＝(λQ^-1+F^HF)QF^H(19)

By positive characterization of the matrix in brackets on both sides, there are

(λQ^-1+F^HF)^-1F^H＝QF^H(λI_M+FQF^H)^-1(20)

Whereby formula (11) can be written as

Y＝QF^H(λI_M+FQF^H)^-1X(21)

The advantage of using equation (15) instead of equation (11) is that only the inverse of the M × M dimensional matrix needs to be solved, reducing the amount of computation.

Step three: model parameter broadening treatment:

and step two, a high-resolution spatial spectrum estimation algorithm based on a Cauchy-Gaussian model has the advantages that the covariance matrix does not need to be smoothed by utilizing multi-data snapshots, and high-resolution orientation estimation is realized in single-frequency-domain data snapshots. The above algorithm requires estimating parametersAnd

the estimation accuracy of which will directly affect the performance of the azimuth spectrum estimation. When in useFar larger than the size of the angular spectrum, the spectral output of equation (15) will degrade to the conventional spatial spectrum, whereas, suitably, it

The values can recover the sparsely distributed characteristics of the signal in space, so that the multi-target directions which are closely spaced can be more finely distinguished.

According to the above analysis, when

Sometimes, resolution can be improved, and thus, the use ofAnd (3) taking equal weight for the regular term and the minimum mode term in the formula (10), and neglecting the constant term to obtain a suboptimal tolerant orientation spectrum estimation based on the Cauchy-Gaussian model:

\tilde{Q} = diag {{YY}^{H}} - - - (23)

and the lambda in the formula (11) is reserved to meet the positive nature of the matrix, and by utilizing the probability distribution theory, the lambda can approximately take a fixed constant, and the empirical value can take 1/M.

Step four: iterative solution of high-resolution spatial spectra:

the iterative solution steps are as follows:

1) with Y⁽⁰⁾DFT { X } (conventional beamforming) as an initial value, and substituting equation (17) to obtain Q⁽⁰⁾；

2) The ith iteration of

Y⁽ⁱ⁾＝Q^(i-1)F^H(λI_M+FQ^(i-1)F^H)^-1X(24)

The Q matrix is updated using the equation:

Q⁽ⁱ⁾＝diag{Y⁽ⁱ⁾Y^(i)H}；

3) judging a convergence condition:

|J⁽ⁱ⁾-J^(i-1)|/|J^(i-1)|＜ε(25)

The specific implementation mode of the invention is as follows:

(1) transforming the matrix receiving time domain signal to a frequency domain:

X(k)＝[X₀(k)，...，X_M-1(k)]^T(26)

θ_n＝asin{cn/(Nf_kd)}(28)

X(k)＝FY(k)(6)

wherein, F is Fourier transform matrix, and Y (k) beam output signal vector of the k frequency point.

(3) Constructing a regular cost function:

(4) taking a fixed empirical value (such as 1/M) for the adjustment parameter lambda in the formula (11), and then iteratively solving the high-resolution spatial spectrum in the following manner:

2) The ith iteration of

Y⁽ⁱ⁾＝Q^(i-1)F^H(λI_M+FQ^(i-1)F^H)^-1X(30)

The Q matrix is updated using the equation:

Q⁽ⁱ⁾＝diag{Y⁽ⁱ⁾Y^(i)H}；

3) judging a convergence condition:

|J⁽ⁱ⁾-J^(i-1)|/|J^(i-1)|＜ε(31)

Simulation analysis example:

a simulation model: the signal is 3 unit amplitude sinusoidal plane waves, and the spatial normalized wave number (fdsin theta/c) is 0.2, 0.25 and-0.25 respectively; the normalized frequencies were 0.2, 0.2 and 0.35. The signal time length is 100 samples, the standard deviation of Gaussian noise is 0.2, and a Hamming window is added to each channel signal to transform the channel signal into a frequency domain. The frequency-wavenumber spectra of both the CBF and the forgiving CG-SPEC algorithms (noted RCG-SPEC) were compared.

Fig. 1 and 2 show normalized frequency-wavenumber spectra of two algorithms. It can be seen that CBF cannot resolve two signals at frequency 0.2 and the side lobes are high. As can be seen from FIG. 2, the RCG-SPEC algorithm proposed by the present invention has high robustness despite the small amount of clutter and the broadening of the mainlobe.

Analyzing sea test data:

and (3) checking the performance of the algorithm by using underwater high-resolution three-dimensional imaging sea test data. The transmitting signal is linear frequency modulation, and a narrow-band echo azimuth spectrum of the 9-element linear array perpendicular to the navigation direction at a frequency point of 12kHz is observed. The CBF and RCG-SPEC algorithm is used for single-snapshot depth-azimuth slice images at a certain navigation position, and bright points in the slice images are small balls suspended in water. Therefore, compared with the conventional beam forming method, the method provided by the invention effectively improves the azimuth resolution.

Because the aperture of the array is limited by the size of the underwater towed body, and the towed body is in motion, the conventional high-resolution algorithm based on covariance matrix estimation is difficult to be effectively applied, the method provided by the invention provides a feasible solution for the occasions where the physical aperture is limited and higher resolution is needed.

In addition to the above embodiments, any technical solutions formed by equivalent substitutions or equivalent transformations fall within the protection scope of the claims of the present invention.

Claims

1. A novel high-resolution direction estimation method based on a Cauchy-Gaussian model is characterized by comprising the following steps: the method comprises the following steps:

(1) transforming the matrix receiving time domain signal to a frequency domain:

X(k)＝[X₀(k)，...，X_M-1(k)]^T (1)

θ_n＝asin{cn/(Nf_kd)} (3)

X(k)＝FY(k) (6)

(3) constructing a regular cost function:

(4)Y＝(λQ^-1+F^HF)^-1F^HX (11)

wherein,

q is the diagonal matrix:

\tilde{Q} = diag {{YY}^{H}} - - - (17)

2) The ith iteration of

Y⁽ⁱ⁾＝Q^(i-1)F^H(λI_M+FQ^(i-1)F^H)^-1X (5)

The Q matrix is updated using the equation:

Q⁽ⁱ⁾＝diag{Y⁽ⁱ⁾Y^(i)H}；

3) judging a convergence condition:

|J⁽ⁱ⁾-J^(i-1)|/|J^(i-1)|＜ε (6)