Windowed 2D deconvolved multi-beam power spectrum estimation algorithm

Abstract

The invention belongs to the field of sonar signal processing, and discloses a windowed two-dimensional deconvolution multi-beam power spectrum estimation algorithm, the steps of which include precalculating and storing a full set of guidance vectors of each scanning direction (θ, F) according to an array structure and a processing frequency band, namely, a point spread function; performing a windowed Fourier transform using k-point array element time domain data to obtain an n-th batch of array element frequency domain data; performing a windowed beamformer processing using a pre-stored guidance vector set to obtain a multi-beam power spectrum estimation result; performing a two-dimensional deconvolution processing on a multi-beam power spectrum windowed in a space-time domain using an R-L iterative algorithm in a discrete two-dimensional integral form to estimate the target's real energy distribution estimation PB. out ; After the unwinding process is completed, PB out A two-dimensional spectrum peak search is performed, the minimum value among the extreme points is used as the background estimate, and all non-extreme point backgrounds are smoothed to obtain a new multi-beam power spectrum estimate. The method of the invention effectively improves the performance of passive sonar narrowband warning and has strong engineering practice significance.

Description

Technical Field

[0001] The invention belongs to the field of sonar signal processing, and in particular relates to a narrow-band line spectrum detection technology for passive sonar. Background Art

[0002] Affected by the leakage of strong interference energy and features in the air and frequency domains, the passive sonar's detection performance for weak target features is reduced. In order to avoid the error in extracting the line spectrum features of the target of interest caused by the leakage of the air and frequency domain side lobes of strong line spectrum interference, the two-dimensional integral expression of the multi-beam power spectrum is derived, and it is proposed to unwrap the output results of the multi-beam power spectrum in the azimuth-frequency two-dimensional process, so as to reduce the influence of the leakage of the air and frequency domain side lobes. At the same time, in order to avoid the problem that the weak target under low signal-to-noise ratio conditions cannot be correctly detected due to the strong interference side lobes being too strong, it is proposed to first use the beamforming technology of the air and time two-dimensional windowing to obtain the multi-beam power spectrum after the windowing process, and then unwrap the processing results in two dimensions to obtain the multi-beam power spectrum estimation result of the method of the present invention. The proposed algorithm inherits the low sidelobe advantage of the windowing process and has better weak target detection capabilities. At the same time, due to the use of two-dimensional unwrapping processing, the main lobe broadening problem caused by windowing is weakened, the beam resolution and frequency resolution are improved, and it has better resolution and weak target detection capabilities. Finally, the output result background is reasonably smoothed to extract the line spectrum features of the target beam of interest, and the target features with purified line spectrum leakage are obtained. Simulation verification shows that the proposed algorithm can effectively suppress the spatial and frequency domain leakage of strong interference, improve the detection ability and extraction accuracy of the target features of interest, and effectively improve the actual performance of multi-beam power spectrum estimation.

[0003] Passive sonar often uses frequency domain space processing to detect the energy and characteristics of targets, that is, the time domain data received by the array element is Fourier transformed into the frequency domain, and then several sub-bands that meet the narrowband assumption are divided. Each narrow sub-band is subjected to frequency domain beamforming to obtain a multi-beam power spectrum result. Under the condition of setting the detection domain reasonably, the azimuth and line spectrum of the target are detected. In the actual use environment of sonar, there are interferences such as fishing boats and surface ships, and their energy and characteristics are seriously leaked in the air and frequency domains, which affects the effective detection, tracking and feature extraction of passive sonar for the target of interest. In order to solve the problem of weak target detection under strong interference conditions, Capon proposed an MVDR processor, which can adaptively form a null in the interference direction, thereby suppressing the sidelobe leakage of the interference. However, it is sensitive to various mismatches, and its performance is seriously degraded in the actual use environment of sonar with serious mismatches. In order to improve the robustness of adaptive beamformers, researchers have conducted a series of robust adaptive algorithm studies, including diagonal loading, feature space, and covariance matrix reconstruction methods. However, due to the particularity of the sonar system, in actual use environments, array element errors and amplitude phase errors are more serious than those in array signal processing fields such as radar and wireless communications, which reduces the performance of robust adaptive algorithms. TC Yang first pointed out in 2016 that in the cosθ spatial domain, the output of a conventional beamformer is the convolution result of the beam pattern (point spread function) and the spatial delta function [Deconvolved conventional beamforming for a horizontal line array, IEEE Journal of oceanic engineering, 2016]. In this way, the RL iterative algorithm can be used to deconvolute the real spatial energy distribution of the target. Since the spatial distribution of target energy is similar to the delta function, the spatial energy leakage is weakened, and the sidelobe level is lowered while the main lobe is reduced. Since the output of the conventional beamformer is post-processed, the deconvolution process inherits its robustness under various mismatch conditions and obtains good weak target detection and multi-target resolution capabilities.

[0004] The literature [High-resolution power spectral estimation method using deconvolution, IEEE Journal of oceanic engineering, 2018] discussed the impact of spatial windowed deconvolution on the main lobe width. As the number of iterations increases, the main lobe widths of different windowed processes gradually tend to be consistent. However, the improvement of weak target detection performance under strong interference by spatial windowing is ignored. Under strong interference conditions, the energy of its spatial sidelobe leakage can be much greater than the main lobe output of the weak target. Due to the errors caused by noise, guidance vector mismatch and other influences in the iterative processing of deconvolution, simply increasing the number of iterations under this condition cannot effectively detect weak targets and weak line spectra.

[0005] The following is a simulation example: a 32-element uniform linear array is located on the positive half axis of the x-axis, the reference array element is located at the origin of the coordinate system, the array is arranged with half a wavelength, the negative half axis of the x-axis is 0°, and it scans 180° clockwise. The entire airspace is scanned in an equal cosine division method, and the scanning beam is 128. There are 4 independent narrowband signals incident on the array in space, one of which is interference, with an incident direction of 30° and an interference-to-noise ratio of 35dB. The other three are weak targets to be detected, with incident directions of 45°, 93°, and 98°, respectively. The signal-to-noise ratios are -5dB, 0dB, and 0dB, respectively, and the number of snapshots is set to 100. CBF, CBF-chebwin, dcv, and dcv-chebwin represent conventional beamforming, spatial Chebyshev weighted beamforming, traditional deconvolution, and Chebyshev weighted post-deconvolution processing, respectively. Figure 7 The normalized spatial energy spectrum results of different processing methods are given. The number of deconvolution iterations is 500.

[0006] It can be seen that due to the influence of the sidelobe leakage of the strong interference at 30°, CBF cannot effectively detect and correctly distinguish the three weak targets at 45°, 93°, and 98°. The Chebyshev windowing process lowers the sidelobe level, but sacrifices the spatial resolution and widens the mainlobe width, making it impossible to distinguish the 30° interference from the 45° target, and at the same time, it is impossible to distinguish the 93° from the 98° target. The Dcv process uses the spatial spectrum of CBF for iterative deconvolution. Since the energy of the interference leakage in the 45° direction is much stronger than the output energy of the 45° target beam, the iterative operation cannot effectively detect the 45° target. Dcv-chebwin performs iterative deconvolution based on the CBF-chebwin process, improves the beam resolution, sharpens the mainlobe width, and can effectively detect and distinguish the 45° target near the 30° interference, and can also correctly distinguish the weak targets at 93° and 98°.

[0007] In view of the above problems, the present invention proposes to use spatial windowing to alleviate the spatial energy leakage of strong interference, and to use time domain windowing to alleviate the frequency domain energy leakage of strong line spectrum. In view of the main lobe widening caused by windowing, it is proposed to perform two-dimensional deconvolution of the multi-beam power spectrum obtained by spatial-temporal windowing beamforming, so as to sharpen the main lobe width of the spatial beam and the main lobe width of the frequency domain beam, reduce the sidelobe level, and suppress the interference energy and feature leakage. The simulation test results show that the proposed algorithm can effectively suppress the energy leakage of interference in the spatial and frequency domains, and improve the detection performance of weak target energy and features under the influence of strong interference. Summary of the invention

[0008] The windowed two-dimensional deconvolved multi-beam power spectrum estimation algorithm comprises the following steps:

[0009] Step 1: Pre-calculate and store each scan direction according to the array structure and processing frequency band The complete set of guiding vectors, i.e., point spread function, defines the negative semi-axis of x as the 0° direction, and scans 180° clockwise to obtain This is a four-dimensional matrix, where represents vector dot product, w s (k) represents the spatial window function, w t (k) represents the time domain window function, It represents the frequency domain data obtained by Fourier transforming a single frequency signal with a frequency of F after time domain windowing;

[0010] Step 2: Use the k-point array element time domain data to perform a windowed Fourier transform to obtain the n-th batch of array element frequency domain data

[0011]

[0012] Where X(k) is the time domain data of the array element, so the cross-spectral density matrix is ​​R(f)=E[Y t (f,n)Y t H (f,n)];

[0013] Step 3: Use the pre-stored steering vector set to perform windowed beamforming processing to obtain the multi-beam power spectrum estimation result

[0014]

[0015] Step 4: Use the discrete two-dimensional integral form of the RL iterative algorithm to perform two-dimensional deconvolution of the spatial-temporal windowed multi-beam power spectrum to estimate the true energy distribution of the target.

[0016]

[0017] in, and

[0018] is the energy-normalized two-dimensional point spread function, Assign the initial value of After i+1 iterations, the final spatial energy distribution estimate is obtained

[0019] Step 5: After the decompression process is completed, out Perform a two-dimensional spectrum peak search and find all the maximum points recorded as Select the minimum value among the extreme points and record it as About PB out All non-extreme points are background-smoothed to obtain a new multi-beam power spectrum estimate

[0020]

[0021] The technical effect of the present invention is that the basic principle and implementation scheme of the present invention have been verified by computer numerical simulation, and the results show that the method of the present invention effectively improves the performance of passive sonar narrowband warning, and has strong engineering practice significance. The present invention derives a two-dimensional discrete integral expression of the multi-beam power spectrum, and proposes a two-dimensional deconvolution of the multi-beam power spectrum obtained by spatial-temporal windowed beamforming, so as to sharpen the main lobe width of the spatial-frequency domain beam, reduce the sidelobe level, and suppress interference energy and feature leakage. The proposed algorithm is not affected by the signal propagation model and array structure, is suitable for the development of various types of passive sonar, and is universal. After verification by computer numerical simulation, the results show that the proposed algorithm improves the spatial resolution and frequency resolution of narrowband warning, can improve the system's detection performance of weak target line spectrum characteristics under the influence of strong interference, and effectively improves the actual performance of passive sonar. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] Figure 1 It is a flow chart of the present invention.

[0023] Figure 2 is pointing to The two-dimensional point spread function of .

[0024] Figure 3 (a) is the multi-beam power spectrum obtained by the CBF method.

[0025] Figure 3 (b) is the multi-beam power spectrum obtained by the CBF-chebwin method.

[0026] Figure 3 (c) is the multi-beam power spectrum obtained by the dcv-chebwin (20 iterations) method.

[0027] Figure 3(d) is the multi-beam power spectrum obtained by the DCV-Chebwin (50 iterations) method.

[0028] Figure 4 (a) is the one-dimensional deconvolved multi-beam power spectrum.

[0029] Figure 4 (b) is the two-dimensional deconvolved multi-beam power spectrum.

[0030] Figure 5 (a) is the spatial spectrum slice corresponding to the normalized frequency 0.036.

[0031] Figure 5 (b) is the spatial spectrum slice corresponding to the normalized frequency 0.0426.

[0032] Figure 6 (a) is the power spectrum slice diagram corresponding to the target beam with an azimuth of 80°.

[0033] Figure 6 (b) is the power spectrum slice diagram corresponding to the target beam with an azimuth of 120°.

[0034] Figure 7 It is the spatial energy spectrum normalized by CBF, CBF-chebwin, DCV, and DCV-chebwin processing methods. DETAILED DESCRIPTION

[0035] The technical solution of the present invention is further described below in conjunction with the accompanying drawings and specific implementation examples:

[0036] First, let’s introduce some preliminary knowledge:

[0037] Without loss of generality, consider the far-field target model of a uniform linear array. A uniform linear array with M array elements is located on the positive half axis of the x-axis. The first array element is the reference array element located at the origin of the coordinate system, and the array element spacing is d. Define the negative half axis of the x-axis as 0° and scan 180° clockwise. Assume that J independent signals in space are sent from θ j The model of the k-th snapshot of the signal received by the array element can be expressed as

[0038]

[0039] Among them, a(θ j )=[1,exp(i2πfdcosθ j / c),…,exp(i2πf(M-1)dcosθ j / c)] T and j (k) respectively represent the steering vector and signal component of the j-th signal, N(k) represents the noise received by the array, f represents the processing frequency, and c represents the sound speed.

[0040] Express the above formula in vector form

[0041] X(k)=A(θ)S(k)+N(k) (2)

[0042] Among them, A(θ)=[a(θ1),a(θ2),…,a(θ J )], S(k)=[s1(k); s2(k);…;s J (k)].

[0043] When the noises between array elements are uncorrelated and of equal power, the theoretical covariance matrix is:

[0044]

[0045] Among them, the symbol [] H and E[] represent conjugate transpose and expectation respectively, represents the jth signal power, represents the array element noise power, and I represents the diagonal unit matrix.

[0046] The output of a conventional beamformer can be expressed as

[0047] P CBF (θ) = a H (θ)Ra(θ) (4)

[0048] In practice, it is difficult to obtain the theoretical covariance matrix, so it is estimated using finite snapshots. The covariance matrix estimate of K snapshots can be expressed as

[0049]

[0050] When the entire space is scanned in an equal cosine division mode, without considering the noise, the output of the conventional beamformer can be expressed as a convolution form

[0051]

[0052] in, Main lobe pointing for conventional beamformer Under this condition, the beam response in each cosθ direction is the point spread function of the conventional beamformer. From the spatial properties of the beam pattern of the conventional beamformer, it can be seen that in the cos(θ) domain Satisfies translation invariance, that is This integral operation can be converted into a convolution operation Among them, * represents convolution. Represents the actual energy distribution of the target in space, as shown below

[0053]

[0054] in, is the signal power, δ() is the delta function, that is, the unit impulse function.

[0055] When the guidance vector scans in the θ domain, the beam pattern does not satisfy the translation invariance, and equation (8) should be expressed in the form of discrete integral:

[0056]

[0057] in, The beam pattern under each main lobe pointing condition after energy normalization. It can be seen that the beam pattern after energy normalization is different under different main lobe pointing conditions. All conditions are met Secondly, the relevant theories and specific contents of each step are explained as follows:

[0058] 1) The relevant theories and specific contents of step 1 are as follows:

[0059] For the problem of multi-beam power spectrum estimation with fixed array structure and pre-defined processing frequency band, each scanning direction can be The two-dimensional point spread function is pre-calculated, stored and retrieved using

[0060]

[0061] It can be seen that this is a four-dimensional matrix, in which represents vector dot product, w s (k) represents the spatial window function, w t (k) represents the time domain window function, It represents the frequency domain data obtained by Fourier transforming the single-frequency signal with frequency F after time-domain windowing. Without loss of generality, this paper adopts Chebyshev weighting suitable for linear arrays. Taking the window function of point K as an example, the expression of point k is

[0062]

[0063] in, SL is the specified sidelobe level.

[0064] 2) The relevant theories and specific contents of step 2 are as follows:

[0065] In order to overcome the problem of weak target energy and feature detection errors increasing and performance degradation caused by strong interference leakage in the spatial and frequency domains, the following two-dimensional windowing in the spatial and time domains is used to obtain the multi-beam power spectrum estimation result.

[0066] The k-point array element time domain data is used for windowed Fourier transform to obtain the n-th batch of array element frequency domain data.

[0067]

[0068] Among them, FFT[] is fast Fourier transform, is the vector dot product, w t (k) is the time domain window function. Then, the cross-spectral density matrix R(f)=E[Y t (f,n)Y t H (f,n)].

[0069] 3) The relevant theories and specific contents of step three are as follows:

[0070] The multi-beam power spectrum estimation result is obtained by using a spatial windowed beamformer

[0071]

[0072] Since the azimuth-beam response of each frequency point is different, the two-dimensional point spread function in the (θ, f) domain does not have the translation invariance property. Formula (13) can be expressed as a two-dimensional integral:

[0073]

[0074] Among them, f l and f h Respectively represent the upper and lower frequency limits of processing, Indicates line spectrum F and the main lobe points to Under these conditions, the two-dimensional directivity response function of each θ direction and line spectrum f is the point spread function. Indicates direction The true energy distribution of frequency F, Indicates power. 4) The relevant theory and specific content of step 4 are as follows:

[0075] The spatial-temporal windowed multi-beam power spectrum is unwrapped in two dimensions to estimate the true energy distribution of the target. Since the beam patterns of different frequencies are different, the point spread function of the multi-beam power spectrum does not satisfy the translation invariance. Therefore, the RL iterative algorithm in the form of discrete two-dimensional integral is used to unwrap the target energy distribution estimation result. The expression is as follows:

[0076]

[0077] in, and is the energy-normalized two-dimensional point spread function. The initial value of After i+1 iterations, the final spatial energy distribution is obtained Similar to the one-dimensional deconvolution process, compared with the two-dimensional deconvolution process in the form of convolution, the iterative deconvolution process in the form of two-dimensional integration does not have the boundary blur problem, so there is no need to expand the data of PB(θ,f).

[0078] 5) The relevant theories and specific contents of step five are as follows:

[0079] After the deconvolution process is completed, some scanning points (θ, f) produce nulls. In order to facilitate the subsequent single-beam feature extraction, the deconvolution multi-beam power background can be smoothed. out Perform a two-dimensional spectrum peak search and find all the maximum points recorded as Select the minimum value among the extreme points and record it as About PB out All non-extreme points are background-smoothed to obtain a new multi-beam power spectrum estimate

[0080]

[0081] For the near-field spherical wave expansion model, the sound field propagation model and other array layout models, the pre-stored two-dimensional point spread function can be used to perform iterative unconvolution processing in integral form to obtain the multi-beam power spectrum estimation of the proposed algorithm. This shows that the proposed algorithm is not affected by the signal expansion model or array structure and is universal.

[0082] Embodiment 1

[0083] According to the process of step 1 to step 5, the performance of the present invention is investigated by simulation. A uniform linear array with 32 array elements is used to verify the performance of the proposed algorithm. The array element spacing is half the wavelength corresponding to the highest processing frequency, and the processing frequency band is normalized to the sampling frequency and ranges from [0.03, 0.05]. Two independent signals in space are incident on the array, and the noise is independent Gaussian white noise between array elements. The target spectrum structure includes broadband continuous spectrum and line spectrum characteristics covering the entire processing frequency band, wherein the interference incident direction is 80°, the interference-noise ratio is 10dB, and it includes a line spectrum component 15dB higher than the continuous spectrum background, and the normalized line spectrum frequency is 0.036. The target incident direction is 120°, the signal-to-noise ratio is 0dB, and it includes a line spectrum component 15dB higher than the continuous spectrum background, and the normalized line spectrum frequency is 0.0426. The spatial scanning beam is divided in an equal cosine manner, and 0° to 180° is 4 times the number of array elements. The frequency domain data of the array element is obtained by Fourier transform using the data of 0.4096s. It can be seen that the frequency resolution is 1 / 0.4096Hz. The number of frequency domain snapshots n=16, and the number of deconvolution iterations is 20. The sidelobe level of Chebyshev weighting is set to -40dB. First, as an example, The two-dimensional point spread function is as follows Figure 2 shown.

[0084] It can be seen that since the beam responses of each frequency point are not exactly the same, the point spread function does not have translation invariance, so it is necessary to use iterative unconvolution processing in the form of integration to obtain multi-beam power estimation. Figure 3 The energy-normalized multi-beam power spectra of CBF, CBF-chebwin, dcv-chebwin (20 iterations), and dcv-chebwin (50 iterations) are given, respectively, showing that the dynamic range is unified at 40 dB.

[0085] It can be seen that the interference features of traditional CBF processing are seriously leaked in the spatial and frequency domains. The feature leakage in the frequency domain will cause errors in the extraction results of the weak line spectrum features of the beam, and the spatial feature leakage will cause errors in the extraction results of the target features of the adjacent beams, that is, the interference features are misjudged as target features, which will affect the target recognition results. Since the target features are stronger than the continuous spectrum background, there are also problems of energy leakage in the spatial and frequency domains, which will cause the adjacent targetless beams to extract the same line spectrum, resulting in false alarms. The CBF-chebwin processing partially suppresses the sidelobe leakage of the interference, but it sacrifices the spatial and frequency resolution, which is manifested as widening the main lobe in the spatial and frequency domains, and reducing the feature resolution results of adjacent beams and adjacent frequencies. The dcv-chebwin processing sharpens the main lobe width in the spatial and frequency domains, improves the resolution, and maintains low sidelobe performance, effectively purifies the spatial and frequency domain leakage of interference energy and features, and improves the detection and feature extraction performance of weak targets. By comparing the multi-beam power spectrum results of dcv-chebwin with different iteration numbers, it can be seen that as the number of iterations increases, the main lobe of the beam becomes closer and narrower, and the side lobe becomes further lower. However, the corresponding amount of calculation will increase, so a compromise can be sought between the main lobe performance and the amount of calculation.

[0086] The following is a three-dimensional diagram of the multi-beam power spectrum obtained by one-dimensional deconvolution of each frequency point in the full spatial domain. Figure 4 As shown in (a), Figure 4 (b) is the power spectrum obtained after the proposed two-dimensional unconvolution.

[0087] It can be seen that the independent one-dimensional unwrapping of each frequency point will cause different spatial noise background heights in different frequency bands, thus forming a spatial depression in the strong line spectrum processing. When the maximum value of each beam is extracted as the line spectrum estimate after the unwrapping process, a false line spectrum caused by the unwrapping process will be generated, causing misjudgment. At the same time, since each frequency point is processed independently, the leakage of the strong line spectrum in the frequency domain cannot be suppressed. Under the condition of insufficient time-frequency analysis points, the influence of frequency domain leakage on the detection of weak line spectrum features cannot be improved. The proposed two-dimensional processing algorithm considers the entire multi-beam power spectrum uniformly, reduces the background level of the beam-frequency dimension as a whole, and performs a two-dimensional extreme point search on the multi-beam power spectrum to obtain the line spectrum azimuth and frequency estimation of the target.

[0088] Figure 5 and Figure 6 The two-dimensional slice diagrams of multi-beam power spectra of different methods are given respectively. Figure 5 is the spatial spectrum slice corresponding to a single frequency point, Figure 6 It is the power spectrum slice corresponding to a single beam. It can be seen that although the windowing process suppresses the sidelobe leakage of the conventional beamformer, it widens the main lobe of the beam and sacrifices some spatial resolution and frequency resolution. Observing the power spectrum results of the 120° beam, compared with the conventional beamformer and windowing process, the method of the present invention effectively suppresses the spatial leakage of the interference line spectrum characteristics. At the same time, due to the use of deconvolution processing, it maintains good spatial-frequency resolution and has the best multi-target, multi-line spectrum resolution, as well as weak target and weak line spectrum detection capabilities.

[0089] According to the implementation examples, the method of the present invention can improve the spatial resolution and frequency resolution of the system in the actual use environment of passive sonar with serious leakage of strong line spectrum interference, effectively improve the line spectrum detection and feature extraction results of weak targets, and avoid the problems of difficult detection of weak targets and feature misjudgment caused by leakage of spatial and frequency domain features.

Claims

1. A windowed two-dimensional deconvolved multi-beam power spectrum estimation algorithm, characterized by: The following steps are involved: Step 1: Pre-calculate and store the full set of steering vectors for each scanning direction (θ, F), i.e., point spread function, according to the array structure and processing frequency band; Step 2: Perform a windowed Fourier transform on the k-point array element time domain data to obtain the n-th batch of array element frequency domain data, and then obtain the cross-spectral density matrix R; Step 3: Use the pre-stored steering vector set to perform windowed beamforming processing to obtain a multi-beam power spectrum estimation result; Step 4: Use the discrete two-dimensional integral form of the RL iterative algorithm to perform two-dimensional deconvolution of the spatial-temporal windowed multi-beam power spectrum to estimate the final spatial energy distribution of the target PB out ; Step 5: After the decompression process is completed, out Perform a two-dimensional spectrum peak search, use the minimum value among the extreme points as the background estimate, smooth the background of all non-extreme points, and obtain a new multi-beam power spectrum estimate; The two-dimensional deconvolution processing formula in step 4 is in, and is the energy-normalized two-dimensional point spread function.

2. According to the windowed two-dimensional deconvolved multi-beam power spectrum estimation algorithm described in claim 1, it is characterized in that: the point spread function in step one defines the negative semi-axis of x as the 0° direction, scans 180° clockwise, and obtains a four-dimensional matrix.

3. The windowed two-dimensional deconvolved multi-beam power spectrum estimation algorithm according to claim 2, characterized in that: The four-dimensional matrix is in Represents vector dot product, w s (k) represents the spatial window function, w t (k) represents the time domain window function, It represents the frequency domain data obtained by Fourier transforming a single-frequency signal with a frequency of F in the time domain after windowing it.

4. The windowed two-dimensional deconvolved multi-beam power spectrum estimation algorithm according to claim 1, characterized in that: The windowed Fourier transform formula in step 2 is Where X(k) is the time domain data of the array element, and the cross-spectral density matrix is ​​R(f)=E[Y t (f,n)Y t H (f, n)].

5. The windowed two-dimensional deconvolved multi-beam power spectrum estimation algorithm according to claim 1, characterized in that: The processing formula in step 3 is 6. The windowed two-dimensional deconvolved multi-beam power spectrum estimation algorithm according to claim 1, characterized in that: for Assign the initial value of After i+1 iterations, the final spatial energy distribution estimate is obtained 7. The windowed two-dimensional deconvolved multi-beam power spectrum estimation algorithm according to claim 1, characterized in that: Search out all the maximum points and record them as Select the minimum value among the extreme points and record it as About PB out All non-extreme point backgrounds are smoothed to obtain a new multi-beam power spectrum estimate.

8. The windowed two-dimensional deconvolved multi-beam power spectrum estimation algorithm according to claim 7, characterized in that: The new multi-beam power spectrum is estimated as

Images