A non-cooperative sonar signal recognition method based on multi-domain feature joint processing

Abstract

The application discloses a non-cooperative sonar signal recognition method based on multi-domain feature joint processing, which comprises the following steps: S10, acquiring underwater active sonar signal time domain sampling data; S20, calculating the power spectral density function of the underwater active sonar signal time domain sampling data; S30, recognizing the signal according to the obtained power spectral density function and the cyclic spectrum domain feature calculated and acquired. The method solves the problems of low recognition rate and low robustness of the underwater active sonar signal in a complex marine environment, has small calculation amount, can recognize multiple signal types (CW, BPSK, QPSK, MSK, 2FSK, 4FSK, LFM, NLFM), and can be suitable for low-power non-cooperative underwater active sonar signal recognition occasions in engineering.

Description

Technical Field

[0001] This invention belongs to the field of signal processing, specifically relating to a non-cooperative sonar signal recognition method based on multi-domain feature joint processing. Background Technology

[0002] The ability to stably identify non-cooperative active underwater sonar signals in complex marine environments has been a long-standing research focus for scholars in the field of underwater acoustics.

[0003] Signal recognition methods can generally be divided into two main categories: decision theory methods based on maximum likelihood hypothesis testing and statistical pattern recognition methods based on feature extraction. When using maximum likelihood to identify signals, the likelihood function of the received signal is first calculated, and then compared with a threshold value to identify the signal type. Since maximum likelihood-based methods are based on the minimum false positive cost criterion, they are optimal classification algorithms in the Bayesian sense. However, two problems exist: firstly, the maximum likelihood recognition algorithm requires a large number of signal parameters; secondly, the existence of unknown parameters makes the calculation expression of the likelihood ratio function very complex. Simplifying the likelihood ratio function can lead to the loss of classification information and reduced recognition performance; therefore, many suboptimal algorithms have been proposed.

[0004] Feature extraction-based statistical pattern recognition methods extract time-domain features or features from various transform domains of a signal, such as the frequency domain, wavelet domain, and cyclic spectrum domain, and then identify the signal by observing the differences in these features. Although feature extraction-based statistical pattern recognition methods may not be optimal, they are easy to implement in engineering, and with proper design, their performance can approach that of the best.

[0005] Kim et al. proposed a method for identifying BPSK and QPSK signals based on decision theory. However, this method has poor robustness and requires prior information about the signal, such as carrier frequency, initial phase, and symbol rate. Subsequently, Hsue, SZ et al. distinguished between FSK and PSK signals by extracting the variance feature of the zero-crossing interval. The zero-crossing interval is essentially a measure of the instantaneous frequency of the signal. In FSK signals, the zero-crossing interval is a step function, while in PSK signals, it is a constant. Simulations show that this feature can distinguish between FSK and PSK signals in low signal-to-noise ratio environments. Next, Nandi et al. proposed a feature extraction-based method to identify 2ASK, 4ASK, BPSK, QPSK, 2FSK, and 4FSK signals. The main extracted feature is the maximum value γ of the zero-center normalized instantaneous amplitude spectral density. max The standard deviation σ of the absolute value of the instantaneous phase nonlinear component of the zero-center non-weak signal segment apThe standard deviation σ of the instantaneous phase nonlinear component of the zero-center non-weak signal segment dp Spectral symmetry P, standard deviation σ of the zero-center normalized instantaneous amplitude absolute value aa The standard deviation σ of the absolute value of the instantaneous frequency of the non-weak signal segment normalized to zero center af The channel is assumed to be a Gaussian channel, and the classifier is a tree classifier. Simulations show that when the signal-to-noise ratio is greater than 15dB, the overall recognition rate of the identified signal can reach more than 90%.

[0006] In recent years, with the development of machine learning theory, many signal recognition classifiers based on machine learning theory have been proposed. B. Kim et al. used 21 signal statistical features and a three-layer neural network classifier to identify BPSK, QPSK, 8PSK, 16QAM, and 64QAM signals. When the signal-to-noise ratio (SNR) is greater than 10 dB, the method achieves an overall signal recognition rate of over 95%. Weihua Jiang et al. combined principal component analysis features with artificial neural networks to identify underwater acoustic non-cooperative communication signals. The signal types they could identify were BPSK, QPSK, and MFSK signals. Simulations showed that when the SNR is greater than 5 dB, the method achieves an overall signal recognition rate of over 80%. Afan Ali et al. extracted constellation diagram features of signals and combined them with a deep learning network to identify BPSK, 4QAM, 16QAM, and 64QAM signals. Simulations showed that when the SNR is greater than 5 dB, the method achieves an overall signal recognition rate of over 95%. These machine learning-based recognition methods still perform well even in low signal-to-noise ratio environments. However, most of these methods are computationally intensive, require a large number of training samples, and are prone to overfitting, making them difficult to use in practice.

[0007] As can be seen from the above introduction, there is currently a lack of research on unified identification methods for non-cooperative underwater acoustic active sonar signals, such as CW, BPSK, QPSK, 2FSK, 4FSK, MSK, LFM, and NLFM, and the need to study the identification of these underwater acoustic signals is quite urgent. Summary of the Invention

[0008] To address the shortcomings of existing technologies, this invention proposes a non-cooperative sonar signal recognition method based on multi-domain feature joint processing. This method is highly stable, can identify a wide variety of signal types, and has low computational complexity, making it suitable for various low-power engineering applications.

[0009] The technical solution for implementing this invention is: a non-cooperative sonar signal identification method based on multi-domain feature joint processing, comprising the following steps:

[0010] S10: Acquire time-domain sampling data of underwater active sonar signal, then proceed to S20.

[0011] S20. Obtain the power spectral density function of the time-domain sampled data of the underwater active sonar signal and normalize it.

[0012] S30. Identify the signal based on the normalized power spectral density function and the calculated cyclic spectral domain characteristics.

[0013] Compared with the prior art, the significant advantages of this invention are: it makes full use of the inherent characteristics of underwater active sonar signals, can identify a variety of signal types, has strong stability, and solves the problem of difficulty in determining the number of true line spectra in the signal power spectral density function under low signal-to-noise ratio environments. At the same time, it solves the problem of poor signal recognition effect caused by the instability of the cyclic spectrum cross-section characteristics due to unreasonable setting of the number of frequency smoothing points. It is suitable for rapid and robust identification of non-cooperative underwater active sonar signals in engineering. Attached Figure Description

[0014] Figure 1 This is a flowchart of the method of the present invention.

[0015] Figure 2 This is a signal recognition performance diagram according to an embodiment of the present invention. Detailed Implementation

[0016] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of them. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0017] Furthermore, the technical solutions of the various embodiments of the present invention can be combined with each other, but only if they are feasible to those skilled in the art. If the combination of technical solutions is contradictory or cannot be implemented, it should be considered that such combination of technical solutions does not exist and is not within the scope of protection claimed by the present invention.

[0018] The following section will further introduce the specific implementation method, as well as the technical difficulties and inventive points of this invention, using this design example as an example.

[0019] Current research on unified identification methods for CW, BPSK, QPSK, 2FSK, 4FSK, MSK, LFM, and NLFM underwater acoustic signals in non-cooperative active sonar signal identification methods is relatively lacking, highlighting the urgent need for such unified identification methods. The non-cooperative active sonar signal identification method of this invention fully utilizes the inherent properties of signals, can identify a wide range of signal types, requires minimal computation, and solves the problem of difficulty in determining the true number of line spectra in the signal power spectral density function under low signal-to-noise ratio environments. Furthermore, it addresses the issue of poor signal identification performance caused by unstable cyclic spectrum cross-sectional characteristics due to unreasonable frequency smoothing point settings. The algorithm exhibits strong stability.

[0020] Reference Figure 1 The present invention provides a non-cooperative sonar signal identification method based on multi-domain feature joint processing, comprising the following steps:

[0021] S10. Obtain time-domain sampling data of underwater active sonar signals, as follows:

[0022] S11. The underwater active sonar signals from N sampling points received in real time from the sensor or stored in the memory are used as the time-domain sampling data x(n) of the underwater active sonar signals. The number of sampling points n = 0, 1, ..., N-1, where the data sampling frequency is denoted as f. s N is usually chosen as an integer power of 2.

[0023] Transfer to S20.

[0024] S20. Obtain the power spectral density function of the time-domain sampled data of the underwater active sonar signal and normalize it, as follows:

[0025] S21. Perform a discrete Fourier transform on the time-domain sampled data x(n) of the underwater active sonar signal to obtain discrete data X(k). The calculation process is as follows:

[0026]

[0027] Where j represents the imaginary part and k represents the discrete frequency index.

[0028] S22. The power spectral density function P(k) of x(n) is obtained from X(k). The calculation process is as follows:

[0029] P(k)=|X(k)| 2 / N,k=0,1,…,N / 2-1

[0030] Here, |·| represents the modulo operation.

[0031] S23. Normalize P(k) to obtain the normalized power spectral density function P2(k). The calculation process is as follows:

[0032]

[0033] Here, max{·} represents the operation of finding the maximum value.

[0034] Switch to S30.

[0035] S30. The signal is identified based on the normalized power spectral density function and the calculated cyclic spectral domain features, as follows:

[0036] S31. Count the number of line spectra in P2(k) whose amplitude is greater than the threshold Th1, and denote it as Num1. If Num1 = 0, then x(n) is determined to be an unknown signal and the signal identification process ends. If Num1 = 1, then x(n) is determined to be a CW signal. Otherwise, go to S32.

[0037] S32. If Num1≥2, then search for the line spectrum frequency index k0 corresponding to the maximum amplitude of P2(k), that is:

[0038]

[0039] S33. Set the value of P2(k) near k0 to 0 to obtain the modified power spectrum P3(k), k = 0, 1, ..., N / 2-1, that is:

[0040]

[0041] Where B1 represents the number of frequency points that need to be set to 0, and the formula for calculating B1 is as follows: B1 = round{20 / (f s / N)} Here, round{·} represents rounding operation.

[0042] S34. Continue searching for the spectral frequency index k1 corresponding to the maximum amplitude of P3(k), i.e.:

[0043]

[0044] in This represents the discrete frequency index corresponding to the maximum value of P3(k) within the range of 0≤k≤N / 2-1.

[0045] S35. If k1 < k0, then: k1 = k0, k0 = k1.

[0046] S36. Calculate the bandwidth estimate BW1 between the two line spectra, i.e.:

[0047] BW1=k1-k0

[0048] S37. Search for the number of line spectra of P2(k) that are greater than the threshold Th2 near k0-BW1, that is, search for the number of line spectra of P2(k) whose amplitude is greater than Th2 when k is between [k0-BW1-B1, k0-BW1+B1], and denot it as T1.

[0049] S38. Search for the number of line spectra of P2(k) that are greater than the threshold Th2 near k1+BW1, that is, search for the number of line spectra of P2(k) whose amplitude is greater than Th2 when k is between [k1+BW1-B1, k1+BW1+B1], and denot it as T1_2.

[0050] S39. Search for the number of line spectra of P2(k) that are greater than the threshold Th2 in the region between k0 and k1. That is, search for the number of line spectra of P2(k) whose amplitude is greater than Th2 when k is between [(k0+k1) / 2-B1, (k0+k1) / 2+B1], and denote it as T1_3.

[0051] S310. Search for the number of line spectra with amplitude greater than Th2 when P2(k) is between [k0-2BW1-B1,k0-2BW1+B1], and denote it as T1_4.

[0052] S311. Search for the number of line spectra with amplitude greater than Th2 when P2(k) is between [k0+2BW1-B1,k0+2BW1+B1], and denote it as T1_5.

[0053] S312. Calculate half of BW1 and round it to obtain the second bandwidth estimate BW2, i.e.:

[0054] BW2 = round{BW1 / 2}

[0055] S313. Search for the number of line spectra with amplitude greater than Th2 when P2(k) is between [k0-BW2-B1, k0-BW2+B1], and denote it as T2_1.

[0056] S314. Search for the number of line spectra with amplitude greater than Th2 when P2(k) is between [k1+BW2-B1,k1+BW2+B1], and denote it as T2_2.

[0057] S315. Search for the number of line spectra with amplitude greater than Th2 when P2(k) is in the range [round{((k0+k1) / 2+k0) / 2}-B1,round{((k0+k1) / 2+k0) / 2}+B1], and denote it as T2_3.

[0058] S316. Search for the number of line spectra with amplitude greater than Th2 when P2(k) is between [round{((k0+k1) / 2+k1) / 2}-B1,round{((k0+k1) / 2+k1) / 2}+B1], and denote it as T2_4.

[0059] S317. Search for the number of line spectra with amplitude greater than Th2 when P2(k) is between [k0+round{(k1-k0) / 3}-B1,k0+round{(k1-k0) / 3}+B1], and denote it as T3.

[0060] S318. Search for the number of line spectra with amplitude greater than Th2 when P2(k) is between [k0+2*round{(k1-k0) / 3}-B1,k0+2*round{(k1-k0) / 3}+B1], and denote it as T3_1.

[0061] S319, Search P2(k)

[0062] [(2k0+3*round{(k1-k0) / 3}) / 2-B1,(2k0+3*round{(k1-k0) / 3}) / 2+

[0063] The number of line spectra with amplitudes greater than Th2 between B1] is denoted as T3_2.

[0064] S320, Search P2(k)

[0065] The number of line spectra with amplitudes greater than Th2 between [round{(2k0+round{(k1-k0) / 3}) / 2}-B1,round{(2k0+round{(k1-k0) / 3}) / 2}+B1], denoted as T3_3.

[0066] S321. Search for the number of line spectra with amplitude greater than Th2 when P2(k) is in the range [round{(k0+2*round{(k1-k0) / 3}+k1) / 2}-B1,round{(k0+2*round{(k1-k0) / 3}+k1) / 2}+B1], and denote it as T3_4.

[0067] S322. Determine whether the following conditions are true or false, i.e., determine:

[0068] ((T1>0)&(T1_3<1))|((T1_2>0)&(T1_3<1))|((T1_4>0)&(T1_3<1))|((T1_5>0) &(T1_3<1))|((T2>0)&(T2_3<1)&(T2_4<1))|((T2_1>0)&(T2_3<1)&(T2_4<1))|

[0069] ((T2_2>0)&(T2_3<1)&(T2_4<1))|((T3>0)&(T3_2<1)&(T3_3<1)&(T3_4<1))|((T3_1>0)&(T3_2<1)&(T3_3<1)&(T3_4<1))

[0070] If the condition is met, then x(n) is determined to be a 4FSK signal; otherwise, proceed to S323. Here, & represents the AND operation, and | represents the OR operation.

[0071] S323. Determine whether the following conditions are true or false, i.e., determine:

[0072] If (T1<1)&(T1_2<1)&(T1_3<1) is true, then determine if x(n) is a 2FSK signal; otherwise, go to S324.

[0073] S324. Calculate the cyclic spectrum cross section of x(n) when the cyclic spectrum frequency is zero, and take its positive cyclic frequency portion. The calculation formula is as follows:

[0074]

[0075] Here, α is the cycle frequency of the cyclic spectrum, and its value ranges from 0 < α ≤ f. s / 2, and the cycle frequency resolution is f s / N, where M is the number of frequency smoothing points, and m is an integer. This indicates the floor function.

[0076] S325, Search The number of line spectra with amplitudes greater than the threshold Th3 on the cross section is denoted as Num2.

[0077] S326. If Num2 < 10, square x(n) and remove the mean to obtain the squared signal data xSqrt(n). The calculation process is as follows:

[0078] xSqrt(n)=x 2 (n)-mean{x(n)},n=0,1,…,N-1

[0079] S327. Calculate the power spectral density function of the squared signal data xSqrt(n), take its positive cyclic frequency component, and normalize it to obtain the normalized power spectral density function P of the squared signal data. sqrt (k), k=0,1,…,N / 2-1, its calculation process is shown in S20.

[0080] S328, Search P sqrt(k) The number of line spectra with amplitude greater than the threshold Th2 is denoted as Num3. If Num3 = 1, x(n) is determined to be a BPSK signal; otherwise, x(n) is determined to be an MSK signal.

[0081] S329. If Num2 ≥ 10, then find the cyclic spectrum cross section when the cyclic spectrum frequency of xSqrt(n) is zero. The calculation formula is shown in S324. (Search) The number of line spectra with amplitudes greater than the threshold Th3 on the cross section is recorded as xSqrtNum2; if xSqrtNum2 < 5, then x(n) is determined to be a QPSK signal, otherwise go to S330.

[0082] S330. Search for points in P2(k), k = 0, 1, ..., N / 2-1 whose amplitude is greater than Th3, and denote them as P3(k), k = 0, 1, ..., Num2.

[0083] S331. Find the mean of P3(k), k = 0, 1, ..., round{Num2 / 2}, denoted as P 31 Similarly, find the mean of P3(k), k = , round{Num2 / 2}, ..., Num2, and denote it as P. 32 ;like

[0084] |P 31 -P 32 |<0.1

[0085] If x(n) is an LFM signal, then x(n) is determined to be an NLFM signal; otherwise, x(n) is determined to be an NLFM signal.

[0086] Example 1

[0087] The simulation signal parameters are set as follows: the signal sampling frequency f in the simulation. s =10kHz, signal length N=16384, noise is white noise, Figure 2 The diagram shows the signal recognition probability of the method of the present invention. As can be seen from the results of the embodiments, when the signal-to-noise ratio is greater than 0dB, the overall signal recognition rate of this method can reach 90%, demonstrating good recognition performance.

Claims

1. A non-cooperative sonar signal recognition method based on multi-domain feature joint processing, characterized in that, Includes the following steps: S10: Acquire time-domain sampling data of underwater active sonar signal, then proceed to S20; S20. Obtain the power spectral density function of the time-domain sampled data of the underwater active sonar signal and normalize it; S30. The signal is identified based on the normalized power spectral density function and the calculated cyclic spectral domain features, as follows: S31, Statistics The number of line spectra with an amplitude greater than the threshold Th1 is denoted as Num1; if Num1=0, then... If the signal is unknown, the signal recognition process ends; if Num1=1, then a decision is made. If it is a CW signal, otherwise switch to S32; S32. If Num1 ≥ 2, then search The line spectrum frequency index corresponding to the maximum amplitude ,Right now: ； S33, will exist The nearby values ​​are set to 0 to obtain the modified power spectrum. ,Right now: ； in, This indicates the number of frequency points that need to be set to 0. The calculation formula is as follows: ； here This indicates rounding operations. Where N is the data sampling frequency, and N is the number of sampling points received from the sensor in real time or stored in the memory; S34, Continue searching The line spectrum frequency index corresponding to the maximum amplitude ,Right now: ； in Indicates in Search within range The discrete frequency index corresponding to the maximum value; S35, if Then: exchange and The value; S36. Calculate the bandwidth estimate between the two line spectra. ,Right now: ； S37, Search exist The number of spectral lines in the vicinity that are greater than the threshold Th2, i.e., the search In the In The number of line spectra with amplitudes greater than Th2 between these two points is denoted as T1; S38, Search exist The number of spectral lines in the vicinity that are greater than the threshold Th2, i.e., the search In the In The number of line spectra with amplitudes greater than Th2 between these two points is denoted as T1_2; S39, Search exist and The number of spectral lines greater than the threshold Th2 near the middle region, i.e., the search In the In The number of line spectra with amplitudes greater than Th2 between these two points is denoted as T1_3; S310, Search exist The number of line spectra with amplitudes greater than Th2 between these two points is denoted as T1_4; S311, Search exist The number of line spectra with amplitudes greater than Th2 between these two points is denoted as T1_5; S312, Calculation Half of that, rounded to the nearest whole number, gives the second bandwidth estimate. ,Right now: ； S313, Search exist The number of line spectra with amplitudes greater than Th2 between these two points is denoted as T2_1. S314, Search exist The number of line spectra with amplitudes greater than Th2 between these two values ​​is denoted as T2_2. S315, Search exist The number of line spectra with amplitudes greater than Th2 between these two values ​​is denoted as T2_3. S316, Search exist The number of line spectra with amplitudes greater than Th2 between these two points is denoted as T2_4; S317, Search exist The number of line spectra with amplitudes greater than Th2 between these two points is denoted as T3; S318, Search exist The number of line spectra with amplitudes greater than Th2 between these two points is denoted as T3_1; S319, Search exist The number of line spectra with amplitudes greater than Th2 between these two values ​​is denoted as T3_2. S320, Search exist The number of line spectra with amplitudes greater than Th2 between these two points is denoted as T3_3; S321, Search exist The number of line spectra with amplitudes greater than Th2 between these two points is denoted as T3_4; S322. Determine whether the following conditions are true or false, i.e., determine: ((T1>0)&(T1_3<1))|((T1_2>0)&(T1_3<1))|((T1_4>0)&(T1_3<1))|((T1_5>0)& (T1_3<1))|((T2_1>0)&(T2_3<1)&(T2_4<1))|((T2_1>0)&(T2_3<1)&(T2_4<1))| ((T2_2>0)&(T2_3<1)&(T2_4<1))|((T3>0)&(T3_2<1)&(T3_3<1)&(T3_4<1))|((T3_1>0)&(T3_2<1)&(T3_3<1)&(T3_4<1)); Is it true? If it is true, then determine... If it is a 4FSK signal, otherwise go to S323; here & means AND operation, and | means OR operation; S323. Determine whether the following conditions are true or false, i.e., judge: If (T1<1)&(T1_2<1)&(T1_3<1) is true, then determine whether the following condition holds true. If it is a 2FSK signal, otherwise switch to S324; S324, Calculation The cross section of the cyclic spectrum when the cyclic frequency is zero, and its positive cyclic frequency portion is taken. The calculation formula is as follows: ； here The cyclic frequency of the cyclic spectrum, and its range of values. And the cycle frequency resolution is , For the number of frequency smoothing points, It is an integer. This indicates the floor function; S325, Search The number of line spectra with amplitudes greater than the threshold Th3 on the cross section, denoted as Num2; S326. If Num2 < 10, for Perform a square operation and subtract the mean to obtain the squared signal data. The calculation process is as follows: = ； S327, Calculate the squared signal data The power spectral density function is obtained by taking its positive cyclic frequency component and normalizing it, thus obtaining the power spectral density function of the normalized squared signal data. The calculation process is shown in S20; S328, Search The number of line spectra with an amplitude greater than the threshold Th2 is recorded as Num3. If Num3=1, then a judgment is made. If it is a BPSK signal, otherwise determine... This is the MSK signal; S329. If Num2 ≥ 10, then calculate... Circular spectrum cross section when the circular spectrum frequency is zero The calculation formula is shown in S324. (Search) The number of line spectra with amplitudes greater than the threshold Th3 on the cross section is denoted as xSqrtNum2; if xSqrtNum2 < 5, then... If it is a QPSK signal, otherwise switch to S330; S330, Search Points with a magnitude greater than Th3 are denoted as ; S331, Seek The mean, denoted as Similarly, to obtain The mean, denoted as ; like ； Then determine If it is an LFM signal, otherwise determine... This is an NLFM signal.

2. The non-cooperative sonar signal recognition method based on multi-domain feature joint processing according to claim 1, characterized in that, S10 includes the following steps: S11. The underwater active sonar signals from N sampling points received in real time from the sensor or stored in the memory are used as the time-domain sampling data of the underwater active sonar signals. Number of sampling points The data sampling frequency is denoted as N is an integer power of 2.

3. The non-cooperative sonar signal recognition method based on multi-domain feature joint processing according to claim 2, characterized in that, S20 includes the following steps: S21. Time-domain sampling data of underwater active sonar signals Perform a discrete Fourier transform to obtain discrete data. The calculation process is as follows: ； Where j represents the imaginary part and k represents the discrete frequency index; S22, according to Seeking power spectral density function The calculation process is as follows: ； here This indicates the modulo operation; S23, to Normalization is performed to obtain the normalized power spectral density function. The calculation process is as follows: ； here This indicates the operation of finding the maximum value.

4. The non-cooperative sonar signal recognition method based on multi-domain feature joint processing according to claim 1, characterized in that, In S30, the threshold values ​​are Th1 = 0.4, Th2 = 0.25, and Th3 = 0.6.

Images