Soft decision line spectrum detection method based on measurement information assistance

Abstract

The invention discloses a soft decision line spectrum detection method based on measurement information assistance, belongs to the field of underwater acoustic signal processing, and particularly relates to the soft decision line spectrum detection method based on the measurement information assistance. The objective of the invention is to solve the problems that an active sound intensity likelihood function is difficult to obtain and the traditional TBD technology has model mismatch so that the line spectrum detection accuracy is low in the prior art. The invention discloses a soft decision line spectrum detection method based on measurement information assistance. The method comprises the following steps: 1, synthesizing multichannel data received by a vector hydrophone to obtain active sound intensity; 2, analyzing statistical characteristics of active sound intensity based on the step 1; 3, constructing a line spectrum track-before-detection model based on the first and second steps; 4, estimating noise parameters of the line spectrum track-before-detect model constructed in the step 3 in real time based on a maximum likelihood criterion; and 5, obtaining a line spectrum existence probability based on the noise parameters estimated in real time in the step 4, wherein the line spectrum existence probability is a soft decision detection result.

Description

Technical Field

[0001] The present invention belongs to the field of underwater acoustic signal processing and relates to a line spectrum detection technology that can be applied to fields such as underwater target detection and ship radiation line spectrum detection; specifically, it relates to a soft decision line spectrum detection method based on the assistance of measurement information. Background Art

[0002] The narrowband line spectrum in ship radiated noise has the characteristics of high intensity and good stability, and is often used as an important feature for target detection in passive sonar. However, with the development of vibration reduction and noise reduction technology, the signal-to-noise ratio of underwater target radiation line spectrum has decreased. The currently popular line spectrum detector based on the Neyman-Pearson (NP) criterion ([1]Wang Y, Ma S, Fan Z, et al. Robust DFT-based generalized likelihood ratio test for underwater tone detection [J]. IET Radar, Sonar & Navigation, 2017, 11 (12): 1845-1853. [2]Wan CR, Goh J T, Chee H T. Optimal tonal detectors based on the power spectrum [J]. IEEE Journal of Oceanic Engineering, 2000, 25 (4): 540-552.) is difficult to effectively detect the presence of targets.

[0003] To avoid missed targets caused by traditional line spectrum detectors, Tracking Before Detection (TBD) and vector hydrophone technologies have been applied to weak line spectrum detection. TBD technology eliminates the threshold judgment process and simultaneously detects and tracks the data collected by the hydrophone, thus avoiding the problem of missed weak line spectrum detection. Vector hydrophone technology synthesizes the received multi-channel data to obtain active sound intensity, improving the signal-to-noise ratio of the received signal and further enhancing line spectrum detection capabilities. However, when applying the traditional TBD technology ([3]Zhang D,Gao L,Sun D,et al.Soft-decision detection of weak tonalsfor passive sonar using track-before-detect method[J].AppliedAcoustics,2022,188:108549.) to perform line spectrum detection on active sound intensity, due to the existence of the cross term of the noise and signal product, the statistical characteristics of the detection statistics no longer meet the model assumptions. Using the existing TBD model to perform line spectrum detection on active sound intensity will lead to model mismatch, which will seriously degrade the detection performance.

[0004] Considering the limitations of the above algorithms, it is of great significance to propose a TBD line spectrum detection method that uses active sound intensity as the detection statistic. Summary of the Invention

[0005] The purpose of the present invention is to solve the problem that the active sound intensity likelihood function is difficult to obtain in existing methods and the traditional TBD technology has model mismatch, resulting in low accuracy of line spectrum detection, and to propose a soft decision line spectrum detection method based on measurement information assistance.

[0006] The specific process of a soft decision line spectrum detection method based on measurement information assistance is as follows:

[0007] Step 1: synthesize the multi-channel data received by the vector hydrophone to obtain the active sound intensity;

[0008] Step 2: Analyze the statistical characteristics of active sound intensity based on step 1;

[0009] Step 3: construct a line spectrum pre-detection tracking model based on steps 1 and 2;

[0010] Step 4: Estimate the noise parameters of the line spectrum pre-detection tracking model constructed in step 3 in real time based on the maximum likelihood criterion;

[0011] Step 5: Based on the noise parameters estimated in real time in step 4, the line spectrum existence probability is obtained, which is the soft decision detection result.

[0012] The beneficial effects of the present invention are:

[0013] To improve the detection of weak line spectra, this paper applies TBD technology to vector hydrophone line spectrum detection. Given the difficulty in obtaining the likelihood function for active sound intensity and the model mismatch issues inherent in traditional TBD technology, this paper introduces additional measurement information to aid in obtaining the likelihood function. This results in a Track Before Detect based on Measurement Information Assisted (MIA-TBD) technique.

[0014] The present invention synthesizes the multi-channel signals collected by the vector hydrophone to obtain active sound intensity and improve the signal-to-noise ratio of the received signal. The MIA-TBD technology uses active sound intensity as a detection statistic to increase the signal-to-noise ratio gain brought by vector synthesis, effectively improving the detection accuracy of weak line spectra. Because the detection results of the MIA-TBD technology are characterized by probability, the present invention refers to this detection method as soft decision detection. The present invention belongs to a line spectrum detection technology that can be applied to fields such as underwater target detection and ship radiation line spectrum detection.

[0015] The present invention can achieve a detection probability of 0.6 for line spectrum when the false alarm probability is 0.001 and the P channel spectrum level signal-to-noise ratio is 0dB.

[0016] When the false alarm probability is 0.001 and the detection probability is 0.5, the present invention has a signal-to-noise ratio gain of 1 to 7 dB compared with traditional TBD and NP methods.

[0017] The present invention performs well in processing experimental data. Under conditions of low signal-to-noise ratio and low false alarm probability, compared with traditional TBD and NP methods, the present invention can increase the line spectrum detection probability by 20% to 60%. BRIEF DESCRIPTION OF THE DRAWINGS

[0018] Figure 1 This is a processing flow chart of the soft decision detector of the present invention;

[0019] Figure 2 The flowchart for obtaining active sound intensity using the complex sound intensity method;

[0020] Figure 3 When SNR=0dB, the ROC curves of the method of the present invention, the traditional TBD method and the NP method are shown;

[0021] Figure 4 When the false alarm probability is 0.001, the ROC curves of the method of the present invention, the traditional TBD method and the NP method are shown;

[0022] Figure 5 When processing experimental data, the ROC curve diagram of the method of the present invention, traditional TBD method and NP method is shown. DETAILED DESCRIPTION

[0023] Specific implementation method 1: The specific process of the soft decision line spectrum detection method based on measurement information assistance in this implementation method is as follows:

[0024] Step 1: synthesize the multi-channel data received by the vector hydrophone to obtain the active sound intensity;

[0025] Step 2: Analyze the statistical characteristics of active sound intensity based on step 1 to provide a theoretical basis for subsequent model construction;

[0026] Step 3: construct a line spectrum pre-detection tracking model based on steps 1 and 2;

[0027] Step 4: Estimate the noise parameters of the line spectrum pre-detection tracking model constructed in step 3 in real time based on the maximum likelihood criterion;

[0028] Step 5: Based on the noise parameters estimated in real time in step 4, the line spectrum existence probability is obtained, which is the soft decision detection result.

[0029] Specific embodiment 2: The difference between this embodiment and specific embodiment 1 is that: in step 1, the multi-channel data (p(t), v x (t), v y (t)) is synthesized to obtain the active sound intensity; the specific process is:

[0030] In order to improve the signal-to-noise ratio of the received signal, this method constructs a complex sound intensifier to synthesize the sound pressure and vibration velocity signals to obtain a synthetic sound signal.

[0031] Perform discrete Fourier transform (DFT) on the sound pressure channel signal to obtain the sound pressure spectrum P(f);

[0032] Perform discrete Fourier transform (DFT) on the vibration velocity channel signal to obtain the vibration velocity spectrum V(f);

[0033] The active sound intensity is obtained by conjugating and multiplying the sound pressure spectrum P(f) and the vibration velocity spectrum V(f) and taking the real part.

[0034] 1) When using a vector hydrophone for line spectrum detection, not only sound pressure information but also vibration velocity information can be obtained. To simplify the description, assume that the normalized seawater acoustic impedance constant ρc = 1, where ρ is the seawater density and c is the seawater sound velocity.

[0035] Based on the normalized seawater acoustic impedance constant ρc=1, assuming that the target detected by passive sonar is located in the far field of the receiver, the pitch angle of the far-field radiation signal is 0, and the received signal is expressed as

[0036]

[0037] Where t represents time, s(t) represents the signal radiated by the target, and θ represents the horizontal azimuth angle of the target;

[0038] p(t), v x (t), v y (t) represents the signals received by the sound pressure channel, the x-direction vibration velocity channel, and the y-direction vibration velocity channel, respectively; the x-direction and y-direction are the x-direction and y-direction in the two-dimensional Cartesian coordinate system, where the z-axis of the Cartesian coordinate system points from the sea surface to the seabed, and the xy plane is perpendicular to the z-axis;

[0039] g p (t), g x (t), g y (t) represents the noise received by the sound pressure channel, x-direction vibration velocity channel, and y-direction vibration velocity channel respectively;

[0040] Assume g p The probability density function of (t) has a mean of 0 and a variance of Gaussian distribution;

[0041] Assume g x The probability density function of (t) has a mean of 0 and a variance of Gaussian distribution;

[0042] Assume g y The probability density function of (t) has a mean of 0 and a variance of Gaussian distribution;

[0043] In actual working scenarios, the noise variances of the two vibration velocity channels are generally equal, that is,

[0044] 2) When processing the velocity channel data of the vector hydrophone, the signal v received by the velocity channel in the x direction is x (t) and the signal v received by the y-direction vibration channel y (t) is vector synthesized to obtain the vibration velocity Expressed as:

[0045]

[0046] Where o is the target direction, which is known a priori in this paper;

[0047] represents the unit vector in the x direction in the two-dimensional Cartesian coordinate system, Represents the unit vector in the y direction in the two-dimensional Cartesian coordinate system;

[0048] The vibration speed The modulus value is expressed as

[0049]

[0050] Where g v (t) is the noise, g v (t) = g x (t)cos(o)+g y (t)sin(o),g v (t) has a mean of 0 and a variance of Gaussian distribution; 3)

[0052] Perform Fourier transform on the sound pressure channel p(t) to obtain the spectrum P(f) of the received signal of the sound pressure channel at frequency f;

[0053] Vibration speed The modulus value Perform Fourier transform to obtain the frequency spectrum V of the received signal of the velocity channel at frequency f c (f);

[0054] 4) In order to improve the signal-to-noise ratio of the received signal, this method constructs a complex sound intensity analyzer to measure the spectrum P(f) of the sound pressure and the spectrum V(f) of the vibration velocity. c (f) Perform synthesis to obtain the active sound intensity I c (f); This method will I c (f) As the detection statistic of the line spectrum detection process.

[0055] Other steps and parameters are the same as those in the first embodiment.

[0056] Specific embodiment three: This embodiment differs from specific embodiment one or two in that: in 3) the sound pressure channel p(t) is subjected to Fourier transform to obtain the spectrum P(f) of the sound pressure channel received signal at frequency f;

[0057] Vibration speed The modulus value Perform Fourier transform to obtain the frequency spectrum V of the received signal of the velocity channel at frequency f c (f);

[0058] Expressed as:

[0059]

[0060] Where,

[0061] S(f) represents the spectrum of the signal s(t), where f is the frequency;

[0062] G P (f) is complex Gaussian noise

[13] ,and The mean is 0 and the variance is The complex Gaussian distribution, σ P is the standard deviation of the noise component in P(f), and N is the number of sampling points of discrete Fourier transform;

[0063] G V (f) is complex Gaussian noise, and The mean is 0 and the variance is The complex Gaussian distribution, σ V is the standard deviation of the noise component in V(f).

[0064] Other steps and parameters are the same as those in the first or second embodiment.

[0065] Specific embodiment 4: This embodiment differs from any one of the specific embodiments 1 to 3 in that: in the above 4), in order to improve the signal-to-noise ratio of the received signal, this method constructs a complex sound intensity device to measure the spectrum P(f) of the sound pressure and the spectrum V(f) of the vibration velocity. c (f) Perform synthesis to obtain the active sound intensity I c (f); This method will I c (f) is the detection statistic of the line spectrum detection process; it is expressed as:

[0066]

[0067] Where, V c The conjugate of (f), Re{·}, represents the real part.

[0068] The other steps and parameters are the same as those in the first to third embodiments.

[0069] Specific embodiment 5: This embodiment differs from any one of specific embodiments 1 to 4 in that: in step 2, the statistical characteristics of the active sound intensity are analyzed based on step 1 to provide a theoretical basis for subsequent model construction; the specific process is:

[0070] The likelihood function of active sound intensity follows the Laplace distribution when only noise exists;

[0071] When the signal is superimposed with noise, the active sound intensity likelihood function is a complex multiple integral. The present invention decomposes the active sound intensity into two related random variables U1(f) and U2(f);

[0072] U1(f) obeys Gaussian distribution, and U2(f) obeys Laplace distribution;

[0073] Regardless of whether there is a signal, the likelihood function of the active sound intensity is related to the noise parameters σ of the P channel and V channel. P σ V related.

[0074] 1) Due to the existence of multiple cross-terms of signal and noise products, I c The probability density function of (f) is a complex integral form, and it is difficult to derive a simple function expression through mathematical derivation. c (f) The statistical characteristics of this method are c (f) structure analysis, first of all, S(f), G P (f), G V (f) Expand to plural form:

[0075]

[0076] Where y(f) represents the real part of S(f);

[0077] j represents the imaginary unit, j 2 =-1;

[0078] represents the imaginary part of S(f);

[0079] n P (f) represents G P The real part of (f), Indicates G P the imaginary part of (f);

[0080] n V (f) represents G V The real part of (f), Indicates G V the imaginary part of (f);

[0081] 2) Substitute formula (5) into formula (3) and formula (4) to obtain

[0082] I c (f)=|S(f)| 2 +U1(f)+U2(f) (6)

[0083]

[0084]

[0085] Where U1(f) represents the cross term of the product of signal and noise, and U2(f) represents the noise component (random variable) independent of the signal;

[0086] Obviously, U1(f) and U2(f) are two related random variables, U1(f) is the cross term of the product of signal and noise, and U2(f) is the noise component independent of the signal;

[0087] According to the properties of Gaussian random variables, U1(f) obeys Gaussian distribution. The mean is 0 and the variance is Gaussian distribution;

[0088] U2(f) obeys the parameter σ P σ V The probability density function f of the random variable U2(f) U (u) is:

[0089]

[0090] Where x1 represents the first integral variable of equation (9), x2 represents the second integral variable of equation (9), and u represents the specific value of the noise component U2(f) independent of the signal;

[0091] 3) Order Substituting variables in formula (9), we can get

[0092]

[0093] Where, f U (u) represents the probability density function of the random variable U2(f), ζ represents the first integral variable of formula (10), and r represents the second integral variable of formula (10);

[0094] Depend on get

[0095]

[0096] In the formula, a represents a constant, b represents a constant, and l represents an integral variable;

[0097] That is, the probability density function of U2(f) is subject to the location parameter 0 and the scale parameter σ P σ V Laplace distribution, denoted as La(0,σ P σ V )distributed.

[0098] The other steps and parameters are the same as those in the first to fourth embodiments.

[0099] Specific embodiment 6: This embodiment differs from any one of specific embodiments 1 to 5 in that: in step 3, a tracking model before line spectrum detection is constructed based on steps 1 and 2; the specific process is:

[0100] When using the pre-detection tracking technology for line spectrum detection, it is necessary to build a corresponding line spectrum detection and tracking model. From formula (6), we can know that I c The statistical characteristics of I do not satisfy the generalized Rayleigh distribution (GRD). c When performing line spectrum detection, the performance is seriously degraded. In addition, when the target exists, I c The likelihood function is a complex integral form, which is difficult to simplify into a simple function, and is not conducive to the construction of the pre-detection tracking model. Based on this, this method introduces the detection statistic I c Other measurement information to assist in obtaining I c Likelihood function of this measurement information is P(f)+V c The real and imaginary parts of (f) are obtained, and a line spectrum tracking model is constructed based on this. This method is called the MIA-TBD method.

[0101] (1) Equation of state:

[0102] 1) Assume that the line spectrum state at time k consists of three parts: the line spectrum frequency f at time k k , the spectrum amplitude A at time k after Fourier transform k and the phase at time k

[0103] Since the line spectrum is a relatively stable quantity in actual scenarios, and the amplitude and frequency generally do not change significantly in a short period of time, this method models the spectrum amplitude and line spectrum frequency as a 0-speed model;

[0104] Based on the 0-speed model, assuming that the line spectrum state at time k is Then the state equation is:

[0105] X k =FX k-1 +Q k (12)

[0106]

[0107] Where F represents the state transfer matrix, is the state noise, where q f ,q A , They are frequency noise, amplitude noise and phase noise respectively, and the superscript T means to take the transpose;

[0108] X k-1 Indicates the line spectrum state at time k-1;

[0109] ΔT represents the time interval between two adjacent frames of data (the discrete Fourier transform is piecewise discrete, and one piece of data is 1 frame of data);

[0110] To describe the detection problem in tracking before detection, we set an existence variable E k Indicates whether the line spectrum exists. When the line spectrum exists at time k, E k =1, otherwise E k =0;

[0111] 2) Assume that the transfer process of the existence variable is a discrete-time Markov chain, and the transfer probability of the existence variable is expressed as

[0112] p birth =p(E k =1|E k-1 =0) (14)

[0113] p death =p(E k =0|E k-1 =1) (15)

[0114] Where p birth represents the probability of particle creation, p(E k =1|E k-1 =0) indicates that variable E exists k-1 =0 transfer to E k =1 probability, E k-1 =0 means that the line spectrum does not exist at time k-1;

[0115] p death represents the probability of particle death, p(E k =0|E k-1 =1) indicates that variable E exists k-1 =1 transfer to E k =0 probability, E k-1 =1 means that the line spectrum exists at time k-1;

[0116] Based on p birth and p death Construct a two-step Markov transfer matrix π, which is expressed as:

[0117]

[0118] Therefore, the Markov chain is constructed by specifying p birth and p death to confirm.

[0119] (2) Measurement equation:

[0120] Compared with traditional line spectrum detection methods, the biggest feature of pre-detection tracking technology is that it does not pre-set thresholds for detection, and uses the data received at the current moment as the measurement value for pre-detection tracking. The measurement value is related to the line spectrum state. At time k, the measurement equation is written as

[0121] Z k =h(E k ,X k ,R k ) (17)

[0122] Where Z k Represents the measurement value generated by the state quantity, h(E k ,X k ,R k ) represents the mapping relationship between state quantity and measurement value, represents the measurement noise, which is a four-dimensional Gaussian random variable. is random noise, and The mean is 0 and the variance is Gaussian distribution, The mean is 0 and the variance is Gaussian distribution;

[0123] In this method, according to I c The relationship between (f) and S(f) defines h(·) as:

[0124]

[0125]

[0126]

[0127] Where u 1,k represents the random variable corresponding to U1(f) in equation (6) at time k, u 2,k represents the random variable corresponding to U2(f) in equation (6) at time k;

[0128] u 2,k ~La(0,σ P σ V );

[0129] The mean is 0 and the variance is A random variable, represents the channel noise variance of P(f) V(f) channel noise variance of and state A at time k k squared The product of

[0130] La(0,σ P σ V ) means that the position parameter is 0 and the scale parameter is σ P σ V Laplace distribution, σ P σ V is the standard deviation of the P(f) channel noise σ P and V(f) channel noise standard deviation σ V The product of

[0131] (3) Likelihood function:

[0132] 1) According to formula (18), when E k =0, Z k (f)~La(0,σ P σ V ),at this time

[0133]

[0134] Where Z k (f) represents the active sound intensity I at frequency f at time k c (f), p(Z k (f)|E k =0) indicates E k =0 when Z k (f) likelihood function;

[0135] 2) According to formula (18), when E k =1, Z k The likelihood function of (f) is a complex integral form, let μ(f) be

[0136]

[0137] Where μ(f) represents the random variable corresponding to U1(f) in equation (6), which is another form of U1(f);

[0138] a(f) and b(f) represent additionally introduced measurement information;

[0139] Represents the phase information of the line spectrum;

[0140] a(f) and b(f) are expressed as:

[0141] a(f)=Re(P(f)+V c (f))=2y(f)+n P (f)+n V (f) (23)

[0142]

[0143] Where,

[0144] P(f) represents the spectrum of the received signal of the sound pressure channel at frequency f;

[0145] V c (f) represents the spectrum of the received signal of the vibration velocity channel at frequency f;

[0146] Re() represents the real part operator;

[0147] Im() represents the imaginary part operator;

[0148] Then the active sound intensity I in formula (6) is c (f) is rewritten as

[0149]

[0150] At this time I c (f)-μ(f)~La(0,σ p σ V );

[0151] Z k (f) and I c (f) corresponds to E k =1, Z k (f)-μ k (f)~La(0,σ p σ V )Right now

[0152]

[0153] Where,

[0154] p(Z k (f)|E k =1) indicates E k =1 when Z k (f) likelihood function;

[0155] μ k (f) represents the frequency at time k, which is obtained by equation (22) according to the line spectrum state A. k , The generated variables are μ(f) in Eq. (22); A k represents the signal spectrum amplitude at time k, i.e. |S(f)| in Equation (22); represents the signal spectrum phase at time k, that is,

[0156] 3) Based on formula (21) and formula (26), the likelihood ratio function L(Z k (f),μk (f),E k ).

[0157] The other steps and parameters are the same as those in the first to fifth embodiments.

[0158] Specific embodiment 7: The difference between this embodiment and any one of the specific embodiments 1 to 6 is that in the above 3), the likelihood ratio function L(Z) is obtained based on formula (21) and formula (26). k (f),μ k (f),E k ); expressed as:

[0159]

[0160] Where,

[0161] L(Z k (f),μ k (f),E k ) represents the likelihood ratio function;

[0162] p(Z k (f)|E k =0) indicates E k =0 when Z k (f) Likelihood function.

[0163] The other steps and parameters are the same as those in the first to sixth embodiments.

[0164] Specific embodiment eight: This embodiment differs from any one of specific embodiments one to seven in that: in step four, the noise parameter σ of the tracking model before line spectrum detection is estimated in real time based on the maximum likelihood criterion. P σ V ; The specific process is:

[0165] In the actual line spectrum detection process, the noise standard deviation is time-varying and unknown, so it is necessary to estimate σ in real time. P σ V ;

[0166] The maximum likelihood estimation uses known sample data to estimate the sample parameters that maximize the likelihood function. The present invention uses the maximum likelihood criterion to estimate the noise parameters σ of the P channel and the V channel. P σ V Make an estimate.

[0167] 1) Discretize f into M frequency points, and the observation set is Z k =[Z k (f1),Z k (f2),…,Z k (f m ),…,Z k (fM )] T , the line spectrum is located at the mth frequency point. Except for the mth frequency point, the rest are all noise points, so σ P σ V The value of is estimated in real time by the remaining noise points except the line spectrum;

[0168] At this time E k =0, Z k The likelihood function is

[0169]

[0170] Where,

[0171] p(Z k |E k =0) indicates E k =0 when Z k The joint probability density of

[0172] Z k (f m ) represents f m Active sound intensity at frequency point I c (f) the measured value;

[0173] p(Z k (f b )|E k =0) indicates E k =0 when Z k (f b )’s likelihood function;

[0174] m represents the mth frequency point after f is discretized, and M represents the total number of frequency points after f is discretized;

[0175] b represents the sum index;

[0176] The superscript T means to find the transpose;

[0177] Taking the logarithm of equation (28) and taking the partial derivative we get

[0178]

[0179] Let equation (29) be equal to 0, and the noise parameter σ under maximum likelihood estimation can be obtained P σ V An estimate of ; expressed as:

[0180]

[0181] Where, represents the noise parameter σ at time k P σ V estimated value.

[0182] Other steps and parameters are the same as those in Specific Embodiments 1 to 7-1.

[0183] Specific embodiment 9: This embodiment differs from any one of specific embodiments 1 to 8 in that: the noise parameter σ at time k P σ V Estimated value of Expressed as:

[0184] In practical applications, considering the simultaneous presence of multiple line spectra, if all frequency bins except the mth frequency bin are considered noise, then the frequencies containing other line spectra will also be treated as noise, affecting the estimation results. To address this issue, a reference window of a certain length should be used for noise bin selection. Furthermore, due to the limited number of DFT points, line spectra often have sidelobes of a certain width, so a guard window of a certain length must also be selected in addition to the reference window.

[0185]

[0186] Where W1 is the reference window length, and W2 is the protection window length.

[0187] The other steps and parameters are the same as those in the specific implementation modes 1 to 8-1.

[0188] Specific embodiment 10: This embodiment differs from any one of specific embodiments 1 to 9 in that: in step 5, the line spectrum existence probability is obtained based on step 4, and the line spectrum existence probability is the soft decision detection result; the specific process is:

[0189] The pre-detection tracking algorithm requires an approximate method to estimate the line spectrum state when implemented. The present invention uses particle filtering technology to estimate the line spectrum state; the soft decision detection result is the ratio of the number of particles indicating the existence of the line spectrum to the total number of particles, and its physical meaning is the estimation of the posterior probability of the existence of the line spectrum.

[0190] In actual situations, a frequency band often has multiple line spectra. To achieve multi-line spectrum detection, this method uses the sparse characteristics of line spectra in the frequency domain to divide the frequency band to be detected into M frequency points, and performs soft decision detection of line spectra on these frequency points simultaneously. Assume that the center frequency of the frequency point is {f m ,m=1,2,…M}, the frequency range is (f m -Δf / 2, f m +Δf / 2), Δf represents the frequency interval between adjacent frequency points.

[0191] The implementation of the tracking before detection technology requires an approximate method to complete the estimation of the line spectrum state. This method uses particle filtering technology to estimate the line spectrum state. The soft decision detection result is the existence of the variable E k=1 to the total number of particles. Its physical meaning is the estimation of the posterior probability of the existence of the line spectrum. The recursive process of the line spectrum soft decision detector based on the MIA-TBD algorithm is as follows:

[0192] Step 1. Particle initialization (k=0); the process is:

[0193] First generate N p particles, at the mth frequency point, the initial state and weight of the i-th particle are expressed as

[0194] in represents the prior probability distribution, Indicates the presence of a line spectrum, Indicates that there is no line spectrum. represents the initial state of the i-th particle at the m-th frequency point at the initial moment, represents the weight of the i-th particle at the m-th frequency point at the initial moment;

[0195] Step 2. State prediction (k>0); the process is:

[0196] Before state prediction, we must first calculate the particle existence variables According to formula (16), there are variables There are four types of transfer:

[0197] The demise of staves:

[0198] Line spectrum persists:

[0199] New line music:

[0200] The line spectrum persists:

[0201] Based on existing variables Four transfer situations, status The forecast is divided into three parts:

[0202] For the particle whose line spectrum does not exist at time k (line spectrum disappears, line spectrum does not exist continuously), let the state of particle i be is a zero vector; the particle whose line spectrum does not exist at time k corresponds to the disappearance of the line spectrum and the continuous non-existence of the line spectrum;

[0203] For the particle with a new line spectrum at time k (new line spectrum), let the state of the i-th particle satisfy the initial distribution, that is,

[0204] For the particle with continuous line spectrum at time k (line spectrum continues to exist), let the state of particle i be Transfer according to formula (12);

[0205] Step 3. Calculate particle weights; the process is:

[0206] This method calculates particle weights based on the values of existing variables and likelihood ratio functions;

[0207] When particles have variables When the particle weight Set to 1, when When , the particle weight is calculated by formula (27); it is expressed as:

[0208]

[0209] Finally, normalize all the calculated particle weights so that the sum of all particle weights is 1;

[0210] The normalized weight of each particle is expressed as:

[0211] is the normalized particle weight of the i-th particle;

[0212] Step 4. Resampling; the process is:

[0213] When performing particle filtering, particle degradation may occur. In this case, resampling is needed to suppress particle degradation. In the above steps, the weights corresponding to different particle states are the predicted posterior probability density functions. During the resampling process, the weight distribution of particles is improved by copying as many particles with large weights as possible and deleting particles with small weights.

[0214] The present invention adopts a system resampling method to resample particles. First, the particle weights of all normalized particles are sorted to obtain a weight sequence, and the weight sequence is divided into N according to the particle weight value. p intervals;

[0215] For example, there are 5 particles in total. The normalized particle weight of the first particle is 0.2, the particle weight of the second particle is 0.1, the particle weight of the third particle is 0.3, the particle weight of the fourth particle is 0.3, and the particle weight of the fifth particle is 0.1. The normalized particle weights of all particles are sorted to obtain a weight sequence (0.2, 0.1, 0.3, 0.3, 0.1). The weight sequence is divided into 5 intervals according to the particle weight value: 0-0.2, 0.2-0.3, 0.3-0.6, 0.6-0.9, 0.9-1.

[0216] Generate a random number uniformly distributed between [0,1] and determine if the random number falls within N p If this random number falls in which interval interval, then the state of the vth particle is used and existential variables Replace the state of the i-th particle and existential variables v is the particle index, i is the traversal 1~N p The particle index of

[0217] After resampling, the weights of all particles are expressed as:

[0218] The new particle is represented by As the state at the previous moment in the state transition process when the next frame of data arrives;

[0219] It represents the state of the i-th particle after resampling at the m-th frequency point at time k, It represents the existence variable of the i-th particle after resampling at the m-th frequency point at time k;

[0220] Step 5. State estimation; the process is:

[0221] Three state quantities of the line spectrum: posterior probability density function State quantity Probability of line spectrum existence Make estimates;

[0222] The existence variable of the i-th particle after resampling at the m-th frequency point at time k The state of the i-th particle after resampling at the m-th frequency point at time k The posterior probability density function is estimated to be

[0223]

[0224] The state quantity is estimated to be

[0225]

[0226] Probability of line spectrum existence That is, the soft decision detection result is obtained after resampling. The number of particles is determined by

[0227]

[0228] in, Represents the state of the i-th particle after resampling at the m-th frequency point at time k;

[0229] represents the Dirac function, when This value is 1 when This value is 0;

[0230] Represents the existence variable of the i-th particle after resampling at the m-th frequency point at time k.

[0231] At time k, the MIA-TBD-based line spectrum soft-decision detector can provide the probability of the presence of a line spectrum at the mth frequency point, rather than a binary decision of 0 or 1. Even at low signal-to-noise ratios, the soft-decision detector can produce a larger detection result, making it easier to detect the presence of a line spectrum.

[0232] The other steps and parameters are the same as those in Specific Embodiments 1 to 9.

[0233] The present invention may have many other embodiments. Without departing from the spirit and essence of the present invention, those skilled in the art may make various corresponding changes and modifications based on the present invention, but these corresponding changes and modifications should all fall within the scope of protection of the claims attached to the present invention.

Claims

1. A soft decision line spectrum detection method based on measurement information assistance, characterized by: The specific process of the method is: Step 1: synthesize the multi-channel data received by the vector hydrophone to obtain the active sound intensity; Step 2: Analyze the statistical characteristics of active sound intensity based on step 1; Step 3: construct a line spectrum pre-detection tracking model based on steps 1 and 2; Step 4: Estimate the noise parameters of the line spectrum pre-detection tracking model constructed in step 3 in real time based on the maximum likelihood criterion; Step 5: Based on the noise parameters estimated in real time in step 4, the line spectrum existence probability is obtained, which is the soft decision detection result.

2. The method for detecting soft decision line spectrum based on measurement information according to claim 1, characterized in that: In step 1, the multi-channel data received by the vector hydrophone is synthesized to obtain the active sound intensity; the specific process is: 1) Assume the normalized seawater acoustic impedance constant ρc = 1, where ρ is the seawater density and c is the seawater sound velocity; Based on the normalized seawater acoustic impedance constant ρc=1, assuming that the target detected by passive sonar is located in the far field of the receiver and the pitch angle of the far-field radiation signal is 0, the received signal is expressed as: Where t represents time, s(t) represents the signal radiated by the target, and θ represents the horizontal azimuth angle of the target; p(t), v x (t), v y (t) represents the signals received by the sound pressure channel, the x-direction vibration velocity channel, and the y-direction vibration velocity channel, respectively; the x-direction and y-direction are the x-direction and y-direction in the two-dimensional Cartesian coordinate system, where the z-axis of the Cartesian coordinate system points from the sea surface to the seabed, and the xy plane is perpendicular to the z-axis; g p (t), g x (t), g y (t) represents the noise received by the sound pressure channel, x-direction vibration velocity channel, and y-direction vibration velocity channel respectively; Assume g p The probability density function of (t) has a mean of 0 and a variance of Gaussian distribution; Assume g x The probability density function of (t) has a mean of 0 and a variance of Gaussian distribution; Assume g y The probability density function of (t) has a mean of 0 and a variance of Gaussian distribution; 2) The signal v received by the x-direction vibration channel x (t) and the signal v received by the y-direction vibration channel y (t) is vector synthesized to obtain the vibration velocity Expressed as: Where o is the target direction; represents the unit vector in the x direction in the two-dimensional Cartesian coordinate system, Represents the unit vector in the y direction in the two-dimensional Cartesian coordinate system; The vibration speed The modulus value is expressed as Where g v (t) is the noise, g v (t) = g x (t)cos(o)+g y (t)sin(o),g v (t) has a mean of 0 and a variance of Gaussian distribution; 3) Perform Fourier transform on the sound pressure channel p(t) to obtain the spectrum P(f) of the received signal of the sound pressure channel at frequency f; Vibration speed The modulus value Perform Fourier transform to obtain the frequency spectrum V of the received signal of the velocity channel at frequency f c (f); 4) The spectrum of sound pressure P(f) and the spectrum of vibration velocity V c (f) Perform synthesis to obtain the active sound intensity I c (f); I c (f) As the detection statistic of the line spectrum detection process.

3. The method for detecting soft decision line spectrum based on measurement information according to claim 2, characterized in that: In the above 3), the sound pressure channel p(t) is subjected to Fourier transform to obtain the spectrum P(f) of the received signal of the sound pressure channel at the frequency f; Vibration speed The modulus value Perform Fourier transform to obtain the frequency spectrum V of the received signal of the velocity channel at frequency f c (f); Expressed as: Where S(f) represents the spectrum of the signal s(t), and f is the frequency; G P (f) is complex Gaussian noise [13] ,and The mean is 0 and the variance is The complex Gaussian distribution, σ P is the standard deviation of the noise component in P(f), and N is the number of sampling points of discrete Fourier transform; G V (f) is complex Gaussian noise, and The mean is 0 and the variance is The complex Gaussian distribution, σ V is the standard deviation of the noise component in V(f).

4. The method for detecting soft decision line spectrum based on measurement information according to claim 3, characterized in that: The spectrum of sound pressure P(f) and the spectrum of vibration velocity V in 4) are c (f) Perform synthesis to obtain the active sound intensity I c (f); I c (f) is the detection statistic of the line spectrum detection process; it is expressed as: Where, V c The conjugate of (f), Re{·}, represents the real part.

5. The method for detecting soft decision line spectrum based on measurement information as claimed in claim 4, characterized in that: The statistical characteristics of the active sound intensity are analyzed based on step 1 in step 2; the specific process is: 1) First, replace S(f) and G in formula (3) P (f), G V (f) Expand to plural form: Where, y(f) represents the real part of S(f); j represents the imaginary unit, j 2 =-1; represents the imaginary part of S(f); n P (f) represents G P The real part of (f), Represents G P the imaginary part of (f); n V (f) represents G V The real part of (f), Represents G V the imaginary part of (f); 2) Substitute formula (5) into formula (3) and formula (4) to obtain I c (f)=|S(f)| 2 +U1(f)+U2(f) (6) Where U1(f) represents the cross term of the product of signal and noise, U2(f) represents the noise component independent of the signal; U1(f) obeys Gaussian distribution, The mean is 0 and the variance is Gaussian distribution; U2(f) obeys the parameter σ P σ V The probability density function f of the random variable U2(f) U (u) is: Where x1 represents the first integral variable of equation (9), x2 represents the second integral variable of equation (9), and u represents the specific value of the noise component U2(f) independent of the signal; 3) Order Substituting variables in formula (9), we can get Where, f U (u) represents the probability density function of the random variable U2(f), represents the first integral variable of formula (10), r represents the second integral variable of formula (10); Depend on get In the formula, a represents a constant, b represents a constant, and l represents an integral variable; That is, the probability density function of U2(f) is subject to the location parameter 0 and the scale parameter σ P σ V Laplace distribution, denoted as La(0,σ P σ V )distributed.

6. The method for detecting soft decision line spectrum based on measurement information as claimed in claim 5, characterized in that: In step 3, a tracking model before line spectrum detection is constructed based on steps 1 and 2. The specific process is as follows: (1) Equation of state: 1) Assume that the line spectrum state at time k consists of three parts: the line spectrum frequency f at time k k , the spectrum amplitude A at time k after Fourier transform k and the phase at time k Model the spectrum amplitude and line spectrum frequency as a 0-speed model; Based on the 0-speed model, assuming that the line spectrum state at time k is Then the state equation is: X k =FX k-1 +Q k (12) Where F represents the state transfer matrix, is the state noise, where q f ,q A , They are frequency noise, amplitude noise and phase noise respectively, and the superscript T means to take the transpose; X k-1 Indicates the line spectrum state at time k-1; ΔT represents the time interval between two adjacent frames of data; Set an existing variable E k Indicates whether the line spectrum exists. When the line spectrum exists at time k, E k =1, otherwise Ek=0; 2) Assume that the transfer process of the existence variable is a discrete-time Markov chain, and the transfer probability of the existence variable is expressed as p birth =p(E k =1|E k -1=0) (14) p death =p(E k =0|E k-1 =1) (15) Where p birth represents the probability of particle creation, p(E k =1|E k-1 =0) indicates that variable E exists k-1 =0 transfer to E k =1 probability, E k-1 =0 means that the line spectrum does not exist at time k-1; p death represents the probability of particle death, p(E k =0|E k-1 =1) indicates that variable E exists k-1 =1 transfer to E k =0 probability, E k-1 =1 means that the line spectrum exists at time k-1; Based on p birth and p death Construct a two-step Markov transfer matrix π, which is expressed as: (2) Measurement equation: At time k, the measurement equation is written as Z k =h(E k ,X k ,R k ) (17) Where Z k Represents the measurement value generated by the state quantity, h(E k ,X k ,R k ) represents the mapping relationship between state quantity and measurement value, represents the measurement noise, is random noise, and The mean is 0 and the variance is Gaussian distribution, The mean is 0 and the variance is Gaussian distribution; According to I c The relationship between (f) and S(f) defines h(·) as: Where u 1,k represents the random variable corresponding to U1(f) in equation (6) at time k, u 2,k represents the random variable corresponding to U2(f) in equation (6) at time k; The mean is 0 and the variance is A random variable, represents the channel noise variance of P(f) V(f) channel noise variance of and state A at time k k squared The product of La(0,σ P σ V ) means that the position parameter is 0 and the scale parameter is σ P σ V Laplace distribution, σ P σ V is the standard deviation of the P(f) channel noise σ P and V(f) channel noise standard deviation σ V The product of (3) Likelihood function: 1) According to formula (18), when E k =0, Z k (f)~La(0,σ P σ V ),at this time Where Z k (f) represents the active sound intensity I at frequency f at time k c (f), p(Z k (f)|E k =0) indicates E k =0 when Z k (f) likelihood function; 2) According to formula (18), when E k =1, let μ(f) be: Where, μ(f) represents the random variable corresponding to U1(f) in Equation (6); a(f) and b(f) represent additionally introduced measurement information; Represents the phase information of the line spectrum; a(f) and b(f) are expressed as: a(f)=Re(P(f)+V c (f))=2y(f)+n P (f)+n V (f) (23) Where P(f) represents the spectrum of the received signal of the sound pressure channel at frequency f, V c (f) represents the spectrum of the received signal of the vibration velocity channel at frequency f; Re() represents the real part operator; Im() represents the imaginary part operator; Then the active sound intensity I in formula (6) is c (f) is rewritten as at this time I c (f)-μ(f)~La(0,σ p s V ); Z k (f) and I c (f) corresponds to E k =1, Z k (f)-μ k (f)~La(0,σ p σ V )Right now Where, p(Z k (f)|E k =1) indicates E k =1 when Z k The likelihood function of (f), μ k (f) represents the frequency at time k, which is obtained by equation (22) according to the line spectrum state A. k , Generated variables; 3) Based on formula (21) and formula (26), the likelihood ratio function L(Z k (f),μ k (f),E k ).

7. The method for detecting soft decision line spectrum based on measurement information according to claim 6, characterized in that: In 3) above, the likelihood ratio function L(Z k (f),μ k (f),E k ); expressed as: Where, L(Z k (f),μ k (f),E k ) represents the likelihood ratio function; p(Z k (f)|E k =0) indicates E k =0, Z k (f) Likelihood function.

8. The method for detecting soft decision line spectrum based on measurement information as claimed in claim 7, characterized in that: In the fourth step, the noise parameter σ of the tracking model before line spectrum detection is estimated in real time based on the maximum likelihood criterion. P σ V ; The specific process is: Discretize f into M frequency points, and the observation set is Z k =[Z k (f1),Z k (f2),…,Z k (f m ),…,Z k (f M )] T , the line spectrum is located at the mth frequency point. Except for the mth frequency point, the rest are all noise points, so σ P σ V The value of is estimated in real time by the remaining noise points except the line spectrum; At this time E k =0, Z k The likelihood function is Where, p(Z k |E k =0) indicates E k =0 when Z k The joint probability density of Z k (f m ) represents f m Active sound intensity at frequency point I c (f) the measured value; p(Z k (f b )|E k =0) indicates E k =0 when Z k (f b )’s likelihood function; m represents the mth frequency point after f is discretized, and M represents the total number of frequency points after f is discretized; b represents the sum index; The superscript T means to find the transpose; Taking the logarithm of equation (28) and taking the partial derivative we get Let equation (29) be equal to 0, and the noise parameter σ under maximum likelihood estimation can be obtained P σ V An estimate of ; expressed as: Where, represents the noise parameter σ at time k P σ V estimated value.

9. The method for detecting soft decision line spectrum based on measurement information as claimed in claim 8, characterized in that: The k-time noise parameter σ P σ V Estimated value of Expressed as: Where W1 is the reference window length, and W2 is the protection window length.

10. The method for detecting soft decision line spectrum based on measurement information as claimed in claim 9, characterized in that: In step 5, the line spectrum existence probability is obtained based on step 4, and the line spectrum existence probability is the soft decision detection result; the specific process is: Step 1. Particle initialization (k=0); the process is: First generate N p particles, at the mth frequency point, the initial state and weight of the i-th particle are expressed as in represents the prior probability distribution, Indicates the presence of a line spectrum, Indicates that there is no line spectrum. represents the initial state of the i-th particle at the m-th frequency point at the initial moment, represents the weight of the i-th particle at the m-th frequency point at the initial moment; Step 2. State prediction (k>0); the process is: Before state prediction, we must first calculate the particle existence variables According to formula (16), there are variables There are four types of transfer: The demise of staves: Line spectrum persists: New line music: The line spectrum persists: Based on existing variables Four transfer situations, status The forecast is divided into three parts: For a particle whose line spectrum does not exist at time k, let the state of the i-th particle be is a zero vector; the particle whose line spectrum does not exist at time k corresponds to the disappearance of the line spectrum and the continuous non-existence of the line spectrum; For the particle with the new line spectrum at time k, let the state of the i-th particle satisfy the initial distribution, that is, For the particle whose line spectrum persists at time k, let the state of particle i be Transfer according to formula (12); Step 3. Calculate particle weights; the process is: When particles have variables When the particle weight Set to 1, when When , the particle weight is calculated by formula (27); it is expressed as: Finally, normalize all the calculated particle weights so that the sum of all particle weights is 1; The normalized weight of each particle is expressed as: is the normalized particle weight of the i-th particle; Step 4. Resampling; the process is: First, sort the normalized particle weights of all particles to obtain a weight sequence, and divide the weight sequence into N p intervals; Generate a random number uniformly distributed between [0,1] and determine if the random number falls within N p If this random number falls in which interval interval, then the state of the vth particle is used and existential variables Replace the state of the i-th particle and existential variables v is the particle index, i is the traversal 1~N p The particle index of After resampling, the weights of all particles are expressed as: The new particle is represented by As the state at the previous moment in the state transition process when the next frame of data arrives; It represents the state of the i-th particle after resampling at the m-th frequency point at time k, It represents the existence variable of the i-th particle after resampling at the m-th frequency point at time k; Step 5. State estimation; the process is: Three state quantities of the line spectrum: posterior probability density function State quantity Probability of line spectrum existence Make estimates; The existence variable of the i-th particle after resampling at the m-th frequency point at time k The state of the i-th particle after resampling at the m-th frequency point at time k The posterior probability density function is estimated to be The state quantity is estimated to be Probability of line spectrum existence That is, the soft decision detection result is obtained after resampling. The number of particles is determined by in, Represents the state of the i-th particle after resampling at the m-th frequency point at time k; represents the Dirac function, when This value is 1, when This value is 0; Represents the existence variable of the i-th particle after resampling at the m-th frequency point at time k.

Images