阅读说明：本技术 一种长时波束相位统计特性的线谱检测方法 (Line spectrum detection method for long-time beam phase statistical characteristics ) 是由 吴艳群 张兵兵 郭微 朱家华 王俊 张文 徐国军 胡正良 于 2021-07-13 设计创作，主要内容包括：本发明属于水声工程、海洋工程和声呐技术等领域,具体涉及一种长时波束相位统计特性的线谱检测方法,先通过水平阵列的波束形成获得水平各个方位各个频点输出的复波束序列,在对准目标方位的输出序列信噪比得到了空间增强；进而以第一快拍的复波束序列为参考,计算其它时刻的复波束序列与参考复波束序列的互谱,获得互谱序列。本发明的特点在于首先对水平阵列进行波束形成,在增强线谱信号的信噪比的同时一定程度上抑制了其它方向上的强干扰；同时,基于长时互谱序列的相位统计特性来区分线谱与各类宽带噪声,其检测门限易于确定且相对稳定。(The invention belongs to the fields of underwater acoustic engineering, ocean engineering, sonar technology and the like, and particularly relates to a line spectrum detection method for long-time beam phase statistical characteristics, wherein a complex beam sequence output by each frequency point in each horizontal direction is obtained through beam forming of a horizontal array, and the signal-to-noise ratio of the output sequence aligned to a target direction is spatially enhanced; and further taking the complex beam sequence of the first snapshot as a reference, calculating the cross spectrum of the complex beam sequence at other moments and the reference complex beam sequence, and obtaining a cross spectrum sequence. The method is characterized in that firstly, the horizontal array is subjected to beam forming, so that the signal-to-noise ratio of a line spectrum signal is enhanced, and strong interference in other directions is inhibited to a certain extent; meanwhile, line spectrum and various broadband noises are distinguished based on the phase statistical characteristics of the long-term cross-spectrum sequence, and the detection threshold is easy to determine and relatively stable.)

1. A line spectrum detection method for long-time beam phase statistical characteristics is characterized by comprising the following steps: the method utilizes a hydrophone array composed of M hydrophone array elements to perform frequency domain beam forming on received acoustic signals, solves the long-time cross-spectrum phase variance on each frequency point of a target beam direction as a threshold value, and improves the performance of line spectrum detection, and specifically comprises the following steps:

s1, carrying out Fourier transform on the acoustic signal matrix received by M array elements of the hydrophone array along each column windowing to obtain a frequency spectrum matrix P (f)_i) (ii) a The method comprises the following specific steps:

s1.1 hydrophoneThe sound signal received by M array elements is represented as N_TMatrix p in x M dimension:

wherein each row p of the matrix p₁(t) p₂(t) … p_M(t) denotes the acoustic signal received by M elements of the hydrophone array at time t, t being 1,2, …, N_T(ii) a Let the time sampling rate be f_sAnd the total observation time is T, the total sampling point number of each array element is N_T＝floor[Tf_s]Where floor denotes a round-down operation;

s1.2 arbitrarily one is not more than N_TPositive integer of (1)_winAs the window length, and one of them is less than N_winPositive integer of (1)_overlapAs the window overlap length, each column of data of the matrix p is divided into K floor [ (N)_T-N_win)/(N_win-N_overlap)+1]A window;

s1.3 choosing the kth window of each column of matrix p, i.e., (k-1) (N)_win-N_overlap) Line +1 to line (k-1) (N)_win-N_overlap)+N_winData of a line is processed by N_winDiscrete Fourier transform of the points to obtain a frequency spectrum matrix P (f)_i)，k＝1,2,...,K：

Is provided withThe minimum working frequency f of the array is the frequency axis vector after Fourier transform_minAnd the highest operating frequency f_maxAt f_FFTThe corresponding numbers in the vector are i_minAnd i_max，i_min＝floor[f_minN_win/f_s]+1，i_max＝floor[f_maxN_win/f_s]+ 1; recording the ith frequency point f after the Fourier transform of the kth window in the mth column_iHas a value of P_m(f_iK) in which f_i＝(i-1)f_s/N_win，m＝1,2,...,M，i＝i_min,i_min+1,...,i_max(ii) a Sequentially traversing K windows of the M-th column, and then traversing all the columns to obtain a K multiplied by M-dimensional frequency spectrum matrix P (f)_i)：

Wherein the matrix P (f)_i) Data P (f) of k-th line_iK) data called kth snapshot: p (f)_i,k)＝[P₁(f_i,k),...,P_M(f_i,k)]K is called fast beat number;

s2, for the spectrum matrix P (f)_i) Kth line data P (f)_iK) performing beamforming to obtain a beamformed complex beam sequence b (f)_iTheta, k), and then calculating a broadband wave beam power output result B (theta, k) in the working frequency band of the hydrophone array and a possible target azimuth angle corresponding to the maximum value of the broadband wave beam power output result B (theta, k); the method comprises the following specific steps:

s2.1 pairs of P (f)_iK) performing beamforming to obtain a complex beam sequence

b(f_i,θ,k)＝w^H(f_i,θ)P^T(f_i,k)

In the formula, a superscript H represents the complex conjugate transpose of a matrix, and a superscript T represents the transpose of the matrix; w (f)_iAnd theta) is a guide vector of the hydrophone array under the far-field plane wave model; theta is horizontal azimuth search angle, angle search interval is delta theta, and minimum value and maximum value are theta respectively_minAnd theta_max；

S2.2 Complex Beam sequence b (f) obtained according to S2.1_iθ, k), solving the broadband wave beam power output result in the working frequency band of the hydrophone array:

s2.3, recording all azimuth angles theta corresponding to the maximum B (theta, k)_skS, S represents the number of B (θ, k) local maximum values, i.e., the kth snapshotThe next possible target number;

s3, calculating the azimuth angle theta of the S-th maximum value_skAt a frequency point f_iTime, complex beam sequence b (f)_iθ, k) complex cross-spectrum value C (f) of sequence value of kth snapshot and sequence value of 1 st snapshot_iS, k) and its phase value phi (f)_i,s,k)：

C(f_i,s,k)＝b(f_i,θ_sk,k)b^*(f_i,θ_s1,1)

φ(f_i,s,k)＝arg{C(f_i,s,k)}

Wherein the superscript denotes the complex conjugate of the vector, arg {. denotes the operation of taking the phase of the complex number;

s4, repeating S3 until all K fast beats are traversed to obtain a long-term statistical complex cross-spectrum phase sequence phi (f) with the length of K_iS) and then calculates the variance δ φ (f) of the sequence_i,s)；

The long-term statistical complex cross-spectral phase sequence can be expressed as

φ(f_i,s)＝arg{[C(f_i,s,1),C(f_i,s,2),...,C(f_i,s,K)]^T}

The variance is

WhereinThe average phase of the ith frequency point is;

s5, mixing the ith_minTo i_maxTotal i_max-i_min+1 frequency points all phase variances delta phi (f)_iAverage of reciprocal of s)As a detection threshold value, delta phi (f) at different frequency points is used_iS) are each separately compared withAnd (3) comparison: the value exceeding the threshold value is regarded as a line spectrum, and the value lower than the threshold value is regarded as random noise; wherein the threshold valueCan be expressed as

S6, at the S-th maximum azimuth angle theta_skHere, S3-S5 are repeated until all i are traversed_max-i_min+1 frequency points, and obtaining a line spectrum detection result of the s-th possible target;

and S7, repeating S6 until all S possible targets are traversed, and finishing the line spectrum detection.

2. A line spectrum detection method for long-term beam phase statistics according to claim 1, characterized in that: s2.1, forming a guide vector w (f) of a corresponding hydrophone array by using plane wave beams under a far-field model_iθ) can be expressed as

3. A line spectrum detection method for long-term beam phase statistics according to claim 1, characterized in that: w (f)_iAnd theta) can also be a steering vector of a plane wave beam forming corresponding hydrophone array under a far-field model by a non-equidistant array or an array with other geometrical configurations.

Technical Field

The invention belongs to the fields of underwater acoustic engineering, ocean engineering, sonar technology and the like, and particularly relates to a line spectrum detection method for long-time wave beam phase statistical characteristics.

Background

The ship radiation noise is generally formed by superposition of a continuous spectrum and a line spectrum, wherein the line spectrum noise is caused by the reciprocating motion of an auxiliary machinery, the frequency of the line spectrum noise is low, the line spectrum noise can be remotely transmitted, and the spectrum level is 10-25dB higher than that of the continuous spectrum, so that the line spectrum is an important characteristic for remotely detecting and identifying underwater targets. The line spectrum detection mainly adopts methods such as power spectrum estimation, modern spectrum estimation and the like, and both adopt a hypothesis test line spectrum detection method based on Discrete Fourier Transform (DFT) and construct a detector according to the statistical characteristics of signals. The Low-Frequency line spectrum Analysis history map (Low Frequency Analysis And Recording, LOFAR) is obtained by continuously performing Short-Time Fourier Transform (STFT) calculation on a segment of received data by using the characteristic that a line spectrum signal has Time continuity. The main principle of line spectrum detection in the LOFAR diagram is to utilize the visual accumulation effect of human eyes to automatically track the frequency of continuous line spectrum signals within a certain time.

In order to improve The line spectrum detection capability, Wagstaff et al of The American naval laboratories invented a Wagstaff's Integration Site Processor (WISPR) and a series of derivative processing techniques, which can perform phase compensation on a line spectrum signal with slowly changing phase by using a phase difference method so as to achieve a long-time integration of The line spectrum signal to obtain sufficient processing gain (see: Exploiting phase fluctuation temporal coherence, IEEE J.Ocean. Eng. to, 2004,29(2):498 and 510. and The AWSUMfilter: A20-dB gain fluctuation-based processor, IEEE J.Ocean. Eng. 1997,22(1):110 and 118). Caesar et al, combined with the output gain of the spatial domain array, effectively improve the resolution of the SNR gain for detecting the target line spectral components and the estimation capability of the DOA azimuth by differentially aligning the fluctuation phase of the output signal of the beam domain (see: research on signal fluctuation detection algorithm based on the radiation noise of the target, J. am. electronics & informatics, 2013, 35(4): 844-851).

However, the existing detection method has the following problems: on one hand, the signal-to-noise ratio of the line spectrum fluctuates due to fluctuation of background noise of the marine environment or interference of other ships, and when some part of the line spectrum is submerged in the interference due to low signal-to-noise ratio due to fluctuation of the signal-to-noise ratio, the part of the line spectrum cannot be detected. Although the use of the array can suppress interference signals from other directions in the horizontal plane through the beamforming technology, and the signal-to-noise ratio of the beam domain signal in the desired target direction is improved to some extent, the desired target line spectrum can still be submerged in noise in the presence of strong interference; on the other hand, due to the interference effect of the ocean waveguide, the line spectrum power received by the receiver fluctuates when the target moves, and the fluctuation range can reach dozens of decibels or even twenty decibels, so that the line spectrum can be detected when the energy is high, and cannot be detected when the energy is low. Considering the random jitter of the noise signal phase, the phase variance is large during long-time observation, the line spectrum signal phase is relatively stable, the phase variance is very small during long-time observation, and even if Doppler frequency shift caused by relative motion of a target or a receiver exists, the phase variance of the line spectrum signal can still keep a small value within a certain time.

Disclosure of Invention

The invention aims to solve the technical problem of providing a line spectrum detection method with long-term beam phase statistical characteristics aiming at the defects of insufficient anti-interference capability and difficult determination of detection threshold of the existing line spectrum detection. According to the method, long-time phase quadratic statistical characteristics (namely phase variance) of the beam domain output frequency domain signals in the target direction are used as the detection quantity of the detector, the array gain is obtained, the output signal-to-noise ratio is improved, meanwhile, line spectrums and various broadband noises (including broadband interference energy leaked to the line spectrum direction) can be distinguished according to the phase statistical characteristics of the line spectrums and the noises, and the method has the advantages that the detection threshold is easy to determine and the anti-interference capability is high.

The technical scheme adopted by the invention is as follows:

a line spectrum detection method for long-time beam phase statistical characteristics utilizes a hydrophone array formed by M hydrophone array elements to perform frequency domain beam forming on received acoustic signals, and solves long-time cross-spectral phase variances on each frequency point of a target beam direction to serve as threshold values, so that the line spectrum detection performance is improved. The method comprises the following steps:

s1, carrying out Fourier transform on the acoustic signal matrix received by M array elements of the hydrophone array along each column windowing to obtain a frequency spectrum matrix P (f)_i) (ii) a The method comprises the following specific steps:

s1.1 representing acoustic signals received by M array elements of a hydrophone array as N_TMatrix p in x M dimension:

wherein each row p of the matrix p₁(t) p₂(t) … p_M(t) denotes the acoustic signal received by M elements of the hydrophone array at time t, t being 1,2, …, N_T(ii) a Let the time sampling rate be f_sAnd the total observation time is T, the total sampling point number of each array element is N_T＝floor[Tf_s]Where floor denotes a round-down operation.

S1.2 arbitrarily one is not more than N_TPositive integer of (1)_winAs the window length, and one of them is less than N_winPositive integer of (1)_overlapAs the window overlap length, each column of data of the matrix p is divided into K floor [ (N)_T-N_win)/(N_win-N_overlap)+1]A window;

s1.3 choosing the kth window of each column of matrix p, i.e., (k-1) (N)_win-N_overlap) Line +1 to line (k-1) (N)_win-N_overlap)+N_winData of a line is processed by N_winDiscrete Fourier transform of the points to obtain a frequency spectrum matrix P (f)_i)，k＝1,2,...,K：

Is provided withThe minimum working frequency f of the array is the frequency axis vector after Fourier transform_minAnd the highest operating frequency f_maxAt f_FFTThe corresponding numbers in the vector are i_minAnd i_max，i_min＝floor[f_minN_win/f_s]+1，i_max＝floor[f_maxN_win/f_s]+ 1; recording the ith frequency point f after the Fourier transform of the kth window in the mth column_iHas a value of P_m(f_iK) in which f_i＝(i-1)f_s/N_win，m＝1,2,...,M，i＝i_min,i_min+1,...,i_max(ii) a Sequentially traversing K windows of the M-th column, and then traversing all the columns to obtain a K multiplied by M-dimensional frequency spectrum matrix P (f)_i)：

Wherein the matrix P (f)_i) Data P (f) of k-th line_iK) data called kth snapshot: p (f)_i,k)＝[P₁(f_i,k),...,P_M(f_i,k)]And k is called the fast beat number.

S2, for the spectrum matrix P (f)_i) Kth line data P (f)_iK) performing beamforming to obtain a beamformed complex beam sequence b (f)_iTheta, k), and then calculating a broadband wave beam power output result B (theta, k) in the working frequency band of the hydrophone array and a possible target azimuth angle corresponding to the maximum value of the broadband wave beam power output result B (theta, k); the method comprises the following specific steps:

s2.1 pairs of P (f)_iK) performing beamforming to obtain a complex beam sequence

b(f_i,θ,k)＝w^H(f_i,θ)P^T(f_i,k)

In the formula, a superscript H represents the complex conjugate transpose of a matrix, and a superscript T represents the transpose of the matrix; w (f)_iTheta) is far field flatA guide vector of the hydrophone array under the surface wave model; theta is horizontal azimuth search angle, angle search interval is delta theta, and minimum value and maximum value are theta respectively_minAnd theta_max。

S2.2 Complex Beam sequence b (f) obtained according to S2.1_iθ, k), solving the broadband wave beam power output result in the working frequency band of the hydrophone array:

s2.3, recording all azimuth angles theta corresponding to the maximum B (theta, k)_skS, S represents the number of B (θ, k) maxima, i.e., the number of possible targets at the kth snapshot.

S3, calculating the azimuth angle theta of the S-th maximum value_skAt a frequency point f_iTime, complex beam sequence b (f)_iθ, k) complex cross-spectrum value C (f) of sequence value of kth snapshot and sequence value of 1 st snapshot_iS, k) and its phase value phi (f)_i,s,k)：

C(f_i,s,k)＝b(f_i,θ_sk,k)b^*(f_i,θ_s1,1)

φ(f_i,s,k)＝arg{C(f_i,s,k)}

Where superscript denotes the complex conjugate of the vector and arg {. denotes the operation of phasing the complex number.

S4, repeating S3 until all K fast beats are traversed to obtain a long-term statistical complex cross-spectrum phase sequence phi (f) with the length of K_iS) and then calculates the variance δ φ (f) of the sequence_i,s)。

The long-term statistical complex cross-spectral phase sequence can be expressed as

φ(f_i,s)＝arg{[C(f_i,s,1),C(f_i,s,2),...,C(f_i,s,K)]^T}

The variance is

WhereinIs the average phase of the ith frequency point.

S5, mixing the ith_minTo i_maxTotal i_max-i_min+1 frequency points all phase variances delta phi (f)_iAverage of reciprocal of s)As a detection threshold value, delta phi (f) at different frequency points is used_iS) are each separately compared withAnd (3) comparison: values above the threshold are considered as line spectra and values below the threshold are considered as random noise. Wherein the threshold valueCan be expressed as

S6, at the S-th maximum azimuth angle theta_skHere, S3-S5 are repeated until all i are traversed_max-i_minAnd +1 frequency points, and obtaining a line spectrum detection result of the s-th possible target.

And S7, repeating S6 until all S possible targets are traversed, and finishing the line spectrum detection.

The invention is based on the following principle: the invention obtains the complex wave beam sequence output by each frequency point in each horizontal direction through the wave beam formation of the horizontal array, and the signal-to-noise ratio of the output sequence aligned to the target direction is enhanced spatially. And further taking the complex beam sequence of the first snapshot as a reference, calculating the cross spectrum of the complex beam sequence at other moments and the reference complex beam sequence, and obtaining a cross spectrum sequence. For line spectrum signals, the phase of the cross-spectrum sequence reflects the long-term stability of the phase of the line spectrum signals in the direction of the target beam. For noise or other broadband signals, the phase long-time statistics of the cross-spectral sequence are random. The method is characterized in that firstly, the horizontal array is subjected to beam forming, so that the signal-to-noise ratio of a line spectrum signal is enhanced, and strong interference in other directions is inhibited to a certain extent; meanwhile, the line spectrum is distinguished from various types of broadband noise (including broadband interference energy leaked to the line spectrum direction) based on the phase statistical characteristics of the long-term cross-spectrum sequence, and the detection threshold is easy to determine and is relatively stable.

Compared with the prior art, the invention has the beneficial effects that:

the invention obtains the threshold value of line spectrum detection by performing beam forming on the horizontal array and calculating the cross-spectrum phase variance. Compared with the traditional method, the threshold acquisition method has the advantages of self-adaptability, small calculated amount, good robustness and convenience for engineering realization; compared with the traditional power spectrum method, the method can obtain higher line spectrum signal-to-noise ratio.

Drawings

FIG. 1 is an overall flow diagram of an embodiment of the present invention;

FIG. 2 is a graph of array beamforming azimuth history results;

FIG. 3 is a LOFAR results plot of beam output directed at a desired target direction;

fig. 4 is a comparison graph of the detection result (a) of the LOFAR output in the target beam direction and the detection result (b) of the phase statistical property of the beam cross-spectrum output of the embodiment of the present invention, in which the five-pointed star represents the line spectrum signal emitted by the sound source.

Detailed Description

The invention is described in detail below with reference to the drawings and a specific embodiment of the sea trial data processing. The line spectrum detection method for the long-term beam phase statistical characteristic provided by the invention is characterized in that the line spectrum is judged as the abnormal long-term phase through statistical analysis of the beam output cross-spectrum phase. FIG. 1 shows a general flow chart of an embodiment of the present invention, which will be described in detail with reference to data obtained from a certain sea test.

At sea this timeIn the test, an unequal-pitch array with the number M of hydrophones being 28 is adopted, and the array shape is close to a linear array. Establishing a xoy rectangular coordinate system by taking the 1 st array element of the array as an origin, wherein the coordinate of the m array element is x_m＝[x_m,y_m]^T. The horizontal azimuth search angle theta is defined as the angle with the positive x-axis, the angle search range is from 0 degrees to 359 degrees, and the angle search interval delta theta is 1 degree.

S1, carrying out Fourier transform on the acoustic signal matrix received by M array elements of the hydrophone array along each column windowing to obtain a frequency spectrum matrix P (f)_i) (ii) a The method comprises the following specific steps:

s1.1 representing acoustic signals received by M array elements of a hydrophone array as N_TMatrix p in x M dimension:

wherein each row p of the matrix p₁(t) p₂(t) … p_M(t) denotes the acoustic signal received by M elements of the hydrophone array at time t, t being 1,2, …, N_T(ii) a Let the time sampling rate be f_sAnd the total observation time is T, the total sampling point number of each array element is N_T＝floor[Tf_s]Where floor denotes a round-down operation.

S1.2 arbitrarily one is not more than N_TPositive integer of (1)_winAs the window length, and one of them is less than N_winPositive integer of (1)_overlapAs the window overlap length, each column of data of the matrix p is divided into K floor [ (N)_T-N_win)/(N_win-N_overlap)+1]A window;

s1.3 selects the kth (K ═ 1, 2.., K) window of each column of the matrix p, i.e., (K-1) (N)_win-N_overlap) +1 line to

(k-1) (N)_win-N_overlap)+N_winData of a line is processed by N_winDiscrete fourier transform of the points, obtaining the spectral matrix p (fi):

is provided withThe minimum working frequency f of the array is the frequency axis vector after Fourier transform_minAnd the highest operating frequency f_maxAt f_FFTThe corresponding numbers in the vector are i_minAnd i_maxThen i_min＝floor[f_minN_win/f_s]+1，i_max＝floor[f_maxN_win/f_s]+ 1; recording the ith frequency point f after the k window Fourier transform of the mth column (M is 1,2_iHas a value of P_m(f_iK) in which f_i＝(i-1)f_s/N_win(i＝i_min,i_min+1,...,i_max) (ii) a Sequentially traversing K windows of the M-th column and then traversing all the columns to obtain a K multiplied by M-dimensional frequency spectrum matrix P (f)_i) I.e. by

Wherein the matrix P (f)_i) The data of the k-th line is referred to as data of the k-th snapshot, and k is referred to as a snapshot number.

In this embodiment, the time sampling rate f_s3276.8Hz, the total observation time is T-5 min, and the total sampling point number of each hydrophone is N_T983040. The number of time window points used in the Fourier transform is N_win65536 points, i.e. the observation duration per time window is 20 s. The number of overlapping points with the last data every time the data is selected is N_overlap45875, a total of 47 windows K can be divided into the total sampling time T. Thus, the matrix P (f)_i) Has dimensions of 47 x 28. The time window point number is 65536 points, so the frequency point number after Fourier transformation is 65536 points. The working frequency range of the array is 40Hz to 410Hz, and the vector f is arranged on the frequency axis_FFTThe corresponding frequency point numbers are i respectively_min802 and i_max8202, the total number of frequency points is therefore i_max-i_minAnd +1 ═ 7401.

S2, for the spectrum matrix P (f)_i) To (1) aPerforming beam forming on the k rows of data to obtain a complex beam sequence b (f) after beam forming_iTheta, k) and a broadband beam power output result B (theta, k) within the working frequency band thereof and a possible target azimuth angle corresponding to the maximum value thereof; the method comprises the following specific steps:

s2.1 treatment of P (f)_i,k)＝[P₁(f_i,k),...,P_M(f_i,k)]Is expressed as a spectrum matrix P (f)_i) The kth row of data (namely sound signals of the kth snap of M array elements of the hydrophone array) are subjected to Fourier transform to obtain frequency spectrum signals of the ith frequency point. To P (f)_iK) performing beamforming to obtain a complex beam sequence

b(f_i,θ,k)＝w^H(f_i,θ)P^T(f_i,k)

In the formula, the superscript H represents the complex conjugate transpose of the matrix, and T represents the transpose of the matrix; theta is horizontal azimuth search angle, angle search interval is delta theta, and minimum value and maximum value are theta respectively_minAnd theta_maxThus θ ═ θ_min,θ_min+Δθ,...,θ_max]^T(ii) a Without loss of generality, let θ_min＝-180°，θ_maxIn practical operation, the range of the horizontal azimuth search angle can be narrowed according to the prior information of the effective detection azimuth angle range of the array to the target, for example, the array can only carry out [ -90 DEG, 90 DEG ] due to the existence of left and right ambiguities]Beamforming of the range of angles of incidence, θ_minAnd theta_maxSet at-90 ° and 90 ° respectively; w (f)_iAnd theta) is a hydrophone array steering vector corresponding to the plane wave under the far-field model. In general, plane wave beamforming under far field model forms the steering vector w (f) of the corresponding hydrophone array_iθ) can be expressed as

For the convenience of analysis, the following contents will take the above-mentioned guide vector as an example. Similarly, w (f)_iTheta) or non-equidistant array or array with other geometric configurations forms a corresponding hydrophone array by plane wave beam under a far-field modelColumn steering vectors, expressions for different types of hydrophone array steering vectors are found in the literature: nielsen Sonar Signal processing Artech House, London, 1991.

S2.2 further solving the broadband power output result in the array working frequency band according to the complex beam sequence obtained in S2.1

S2.3 recording all azimuth angles theta corresponding to the maximum B (theta, k)_sk(S ═ 1, 2., S), where S denotes the number of B (θ, k) maxima, i.e., the number of possible targets at the kth snapshot. Preferably, a person skilled in the art may reduce the number of possible targets (for example, assuming that there are 3 maxima recorded in this step, each corresponding to an azimuth angle θ) according to the prior information of the desired target position (for example, the prior information of the desired target position is obtained by an Automatic Identification System (AIS) of the vessel)_1k、θ_2k、θ_3k79 degrees, 80 degrees and 100 degrees respectively, and the prior information according to the experimental process shows that the possible targets only exist near 80 degrees, and when the azimuth angle exclusion threshold is set to be 5 degrees, the maximum value corresponding to 100 degrees of azimuth angles can be excluded, so that the number of the possible targets is reduced to 2).

In this embodiment, the horizontal azimuth spatial scanning angle is 1 to 180 degrees, and the angle interval is 1 degree. Thus, a different frequency point f is obtained_iAnd in the case of different snapshots k, the complex beam sequence in 180-degree direction outputs b (f)_iθ, k). The target azimuth history shown in fig. 2 is obtained in terms of the power output B (θ, k) of the beam sequence. FIG. 2 shows the output of B (θ, k) in 5 minutes. From a priori information in the experimental conditions, the target is expected to be located around 80 degrees. As can be seen from FIG. 2, within 5 minutes, the strong interference azimuth is around 8 degrees, and the desired target angle is substantially constant at 78 degrees, so θ₁₁＝...＝θ_1K78 degrees. According to fig. 1, the interference energy is higher than the target. As shown in FIG. 3, the output beam signal characteristics in the desired target direction, i.e., LOFAR analysis result | b (f)_i,θ_1k,k)|². Comparing fig. 2 and fig. 3, it can be seen that although the interference direction and the target direction are greatly different, the energy of the strong interference still leaks to the target direction through the side lobe, so that the low-frequency part of the LOFAR map in the target direction has strong ship interference energy, and many line spectrum energies are unstable, and cannot be observed from the LOFAR map.

S3, calculating the azimuth angle theta of the S-th maximum value_skAt a frequency point f_iThen, the complex beam sequence has a complex cross-spectrum value C (f) between the sequence value of the kth snapshot and the sequence value of the 1 st snapshot_iS, k) and its phase value phi (f)_i,s,k)：

C(f_i,s,k)＝b(f_i,θ_sk,k)b^*(f_i,θ_s1,1)

φ(f_i,s,k)＝arg{C(f_i,s,k)}

Where superscript denotes the complex conjugate of the vector and arg denotes the operation of taking the phase of the complex number.

In this embodiment, since the number of data points is divided into 47 windows, the total number of frequency points is i_max-i_min+ 1-7401, so that 47 complex cross-spectral phase values phi (f) are obtained in total_iK), i.e. every 47 snapshots the subsequent line spectrum detection is performed. Therefore, for 7401 frequency points in the range of the array working frequency band, the acquired cross-spectrum phase phi (f)_iAnd k) are 47 in number.

S4, repeating S3 until all K fast beats are traversed, obtaining a long-term statistical complex cross-spectrum phase sequence with the length of K, and then calculating the variance delta phi (f) of the sequence_i,s,k)。

The long-term statistical complex cross-spectral phase sequence can be expressed as

φ(f_i,s,k)＝arg{[C(f_i,s,1),C(f_i,s,2),...,C(f_i,s,K)]^T}

The variance is

WhereinIs the average phase of the ith frequency point.

S5, mixing the ith_minTo i_maxTotal i_max-i_min+1 frequency points all phase variances delta phi (f)_iAverage of reciprocal of s)As a detection threshold value, delta phi (f) at different frequency points is used_iS) are each separately compared withAnd (3) comparison: values above the threshold are considered as line spectra and values below the threshold are considered as random noise. Wherein the threshold valueCan be expressed as

S6, at the S-th maximum azimuth angle theta_skHere, S3-S5 are repeated until all i are traversed_max-i_minAnd +1 frequency points, and obtaining a line spectrum detection result of the s-th possible target.

And S7, repeating S6 until all S possible targets are traversed, and finishing the line spectrum detection.

In this embodiment, for each frequency point within the array operating band, the cross-spectrum phase phi (f) of the desired target is obtained_i1, k) mean valueSum variance δ φ (f)_i,1). FIG. 4(b) shows the phase variance δ φ (f) for all frequency bins_i,1). To verify the effectiveness of the present invention, fig. 4 shows a comparison of the conventional power spectroscopy (fig. 4(a)) and the results of the present invention (fig. 4 (b)). Conventional power spectroscopy is performed by taking the desired objective as shown in FIG. 3Averaging of characteristics of output beam signals in a target direction, i.e. taking multiple snapshotsIn the present example, there is a wide-band strong interference target and the signal-to-interference ratio is low, so that it is difficult for the conventional power spectroscopy to detect the line spectrum by the signal energy. Even if the trend term of the power spectrum is smoothed by detrending in fig. 4(a), the signal-to-noise ratio of the line spectrum cannot be improved, and the signal is still difficult to detect. The invention separates the line spectrum signal from the noise by utilizing the statistical law of the phase stability of the line spectrum signal and the statistical law of the broadband noise random phase, rather than adopting the signal energy characteristic, and has better detection effect. Meanwhile, the phase variances of different noise frequency points are very close, namely the noise phase variance spectrum is very flat, which provides convenience for threshold determination and is convenient for system implementation.

Finally, it should be noted that: although the present invention has been described in detail with reference to the foregoing embodiments, it will be apparent to those skilled in the art that modifications may be made to the embodiments described above, or equivalents may be substituted for elements thereof. Any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention should be included in the protection scope of the present invention.