CN111273301A - Frequency spectrum reconstruction method for underwater sound target radiation noise linear array wave beam output signal - Google Patents

Frequency spectrum reconstruction method for underwater sound target radiation noise linear array wave beam output signal Download PDF

Info

Publication number
CN111273301A
CN111273301A CN202010100407.7A CN202010100407A CN111273301A CN 111273301 A CN111273301 A CN 111273301A CN 202010100407 A CN202010100407 A CN 202010100407A CN 111273301 A CN111273301 A CN 111273301A
Authority
CN
China
Prior art keywords
frequency
dictionary
azimuth
sparse
signal
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202010100407.7A
Other languages
Chinese (zh)
Inventor
曾向阳
陆晨翔
乔彦
杨爽
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Northwestern Polytechnical University
Original Assignee
Northwestern Polytechnical University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Northwestern Polytechnical University filed Critical Northwestern Polytechnical University
Priority to CN202010100407.7A priority Critical patent/CN111273301A/en
Publication of CN111273301A publication Critical patent/CN111273301A/en
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S15/00—Systems using the reflection or reradiation of acoustic waves, e.g. sonar systems
    • G01S15/88—Sonar systems specially adapted for specific applications
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S15/00—Systems using the reflection or reradiation of acoustic waves, e.g. sonar systems
    • G01S15/02—Systems using the reflection or reradiation of acoustic waves, e.g. sonar systems using reflection of acoustic waves
    • G01S15/04—Systems determining presence of a target
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/52—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00
    • G01S7/523—Details of pulse systems
    • G01S7/526—Receivers
    • G01S7/527—Extracting wanted echo signals
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/52—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00
    • G01S7/534—Details of non-pulse systems
    • G01S7/536—Extracting wanted echo signals
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/52—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00
    • G01S7/539—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Acoustics & Sound (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)

Abstract

The invention provides a method for reconstructing a frequency spectrum of an output signal of a linear array beam of underwater sound target radiation noise, which comprises the steps of firstly transforming a linear array multi-channel signal to a frequency domain through Fourier transform, then constructing a space azimuth dictionary, respectively constructing observation vectors for each frequency point in a selected frequency band range, adding space response consistency constraint among the frequency points, constructing a multi-frequency-point multi-task sparse decomposition model, finally deducing and estimating posterior mean values of sparse decomposition coefficients by using Bayesian variation, and obtaining frequency spectrum amplitude values of the frequency points in all azimuths including a target azimuth. Compared with the traditional conventional beam forming method, the method can better recover the target frequency spectrum of the medium-high frequency band of 500-1000Hz in the radiation noise signal.

Description

Frequency spectrum reconstruction method for underwater sound target radiation noise linear array wave beam output signal
Technical Field
The invention relates to the field of signal analysis and processing, and can be used in the field of analysis and processing of underwater artificial target radiation signals acquired by linear arrays. The invention can realize the requirement of constant beam width of the beam on different frequency points by adding the structuralized sparse constraint on the basis of the spatial response consistency of different frequency bands while obtaining the array gain by beam forming, thereby reducing the influence of the beam signal distortion caused by the change of the beam width along with the frequency. Compared with the traditional beam forming method, the method can better recover the medium-high frequency signal components of the target radiation noise.
Background
Sonar hydrophone arrays are used for collecting radiation noise signals of underwater artificial targets, and whether the underwater targets exist or not is detected and the characteristics of the underwater targets are analyzed. The underwater artificial target radiation noise signal collected based on the hydrophone array is often a broadband signal, and for the fixed linear array, the beam width becomes narrow as the frequency increases, which causes the distortion of the output broadband radiation noise signal. To solve this problem, it is necessary to ensure that the array beam has a constant beam width, and therefore, it is necessary to design a constant beam width beamformer. The calculation amount required by designing the broadband constant beam width beam by adopting the optimization method is large, and the application in practice is still not limited.
Disclosure of Invention
In order to overcome the defects of the prior art, the invention provides a method for reconstructing a frequency spectrum of an underwater sound target radiation noise linear array beam output signal. The invention realizes the effect of constant beam width through structured sparse constraint, obtains broadband target radiation noise beam signals with smaller distortion degree, and has higher operation efficiency. The method comprises the steps of firstly transforming linear array multi-channel signals to a frequency domain through Fourier transform, then constructing a space orientation dictionary, respectively constructing observation vectors for each frequency point in a selected frequency band range, adding space response consistency constraint among the frequency points, constructing a multi-frequency-point multi-task sparse decomposition model, finally deducing and estimating posterior mean values of sparse decomposition coefficients by using Bayesian variational components, and obtaining frequency spectrum amplitude values of each frequency point in all orientations including a target orientation.
The technical scheme adopted by the invention for solving the technical problem comprises the following detailed steps:
step 1: falseSetting a linear array to have N array elements, wherein time domain signals acquired by each array element channel are x respectively1(n),…,xN(n) for each channel signal x in the array1(n),…,xN(n) performing Fast Fourier Transform (FFT) to obtain frequency domain data xF1(f),xF2(f),…,xFN(f)：
Figure BDA0002386703980000021
In the formula (1), FFT is fast Fourier transform;
step 2: the user selects a signal frequency band needing to be processed for the concerned frequency band, selects a frequency domain signal of a corresponding frequency band in the step 1, the upper frequency limit nf _ h is determined by the actual requirement of the user, and 1000Hz can be selected because the artificial target radiation noise is concentrated below 1000 Hz; if the amplitude modulation spectrum above 1000Hz needs to be paid attention to, the upper limit frequency can be set to 3000-4000Hz, and the specific numerical setting needs to be combined with the sampling frequency of the underwater acoustic equipment and the prior information judgment of a user on the target radiation noise signal; setting a lower frequency limit nf _ l by combining the actual requirements of a user; recording each channel frequency domain signal in the selected formulated frequency band as x'F1(f),…,x'FN(f) Then, there are:
Figure BDA0002386703980000022
and step 3: setting an incident space azimuth vector theta with the length ND according to the linear array aperture size, the resolution requirement and the actual user requirement, and then constructing a space azimuth dictionary D based on the selected space azimuth setf，DfThe number of rows is N, the number of columns is ND, NF — nh — NF — l +1 frequency points are selected in step 2, and correspondingly, each frequency point f corresponds to one dictionary DfTotal NF dictionary, frequency point f corresponding to space azimuth dictionary DfThe expression is as follows:
Figure BDA0002386703980000023
wherein e is a natural index, i is an imaginary number symbol, pi is a circumferential rate, c is a sound velocity, p is a position vector consisting of N array elements, theta is a selected incident space azimuth angle vector, and T represents transposition;
and 4, step 4: to orientation dictionary DfAnd (3) carrying out normalization:
Figure BDA0002386703980000024
in the formula (4), sqrt represents an evolution, and diag () represents an operation of taking diagonal elements of a matrix to form a vector for a matrix operation;
and 5: modeling by adopting a structured sparse model;
for any frequency point f in the frequency band selected in the step 2, the frequency domain observation vector of each array element of the array is Yf＝[xF1(f),xF2(f),…,xFN(f)]Wherein x isFi(f) For observation of frequency point f for ith array element, YfLength is array element number N, azimuth dictionary DfIs NXND, YfThrough space orientation dictionary DfDecomposition coefficient X offThe length is the number ND of incident azimuth angles, and for each frequency point, the following sparse decomposition model is established:
Yf＝DfXf(5)
wherein the target is sparse in space, and thus
Figure BDA0002386703980000031
The medium and large coefficient elements are sparse, so that the formula (5) is YfDictionary D based on spatial orientationfDecomposing to obtain XfA sparse decomposition process of (c);
because of the spatial orientation dictionary DfIs a frequency dependent quantity and is therefore for each frequency point f within the frequency band selected in step 21～fNAll the requirements are as follows:
Figure BDA0002386703980000032
obtained by performing sparse decomposition on all frequency points selected in the formula (6)
Figure BDA0002386703980000033
In, the positions where the large coefficient elements appear should be the same or very close; in order to fully utilize the constraint information, the N sparse decomposition problems in the formula (6) are jointly solved; setting joint observation vector on all frequency points
Figure BDA0002386703980000034
X is Y channel joint space azimuth dictionary
Figure BDA0002386703980000035
D is a three-dimensional matrix, the third dimension corresponds to D at different frequency pointsfThe N sparse decomposition models of equation (6) form a multitask sparse decomposition model by joint solution:
Y＝DX (7)
jointly solving the multitask sparse decomposition model, namely adding the spatial response consistency constraint of each frequency point, wherein the spatial response consistency constraint is reflected in a sparse decomposition coefficient vector X, and a structured sparse constraint is added to the vector X large elements; obtaining a beam signal with the characteristic of constant beam width while carrying out broadband beam scanning;
since the actual signal contains noise, the following model is obtained by considering the noise n:
Y＝DX+n (8)
step 6: deducing and solving the multitask sparse decomposition model represented by the formula (8) by adopting a Bayesian variational algorithm; assuming that the noise n follows a multivariate Gaussian distribution, in equation (9)
Figure BDA0002386703980000041
Representing a multivariate gaussian distribution:
Figure BDA0002386703980000042
equation (9) assumes that the hidden variable λ of the model obeys Gamma distribution based on parameters a and b, Γ represents Gamma distribution, a and b are model adjustable parameters, and the specific value needs to be determined by continuously adjusting according to the characteristics and the actual effect of the actual signal:
p(λ)＝Γ(λ|a,b) (10)
equation (10) assumes the sparse decomposition coefficients to be solved as random vectors in the probabilistic model and assumes XfElement X in (1)kfAre subject to a mean of 0 and an accuracy of gammakAnd independently of each other, XfObeying multivariate Gaussian distribution, wherein the prior probability distribution is as follows:
Figure BDA0002386703980000043
for the accuracy parameter, it is still assumed that γkObey a Gamma prior distribution and are independent of each other, Gamma-Gamma1,γ2,...,γND]Is given by the following equation:
Figure BDA0002386703980000044
and 7: according to the step 6, deducing hidden variables lambda and gamma in the multitask sparse decomposition model, and finally taking the posterior probability mean value of the decomposition coefficient X as an estimated value of the decomposition coefficient X in the formula (7);
and 8: finding the position of large coefficient in X, and mapping the position to dictionary DfAzimuth of corresponding position as appearance azimuth theta of targettThen, the orientation theta of the decomposition coefficient X to the target appearance is selectedtAll frequency point results XθtObtaining an azimuth angle estimation result and a frequency domain beam amplitude of the one-time processing snapshot duration; and finally, processing the continuous time signals repeatedly according to the steps 1-7, and splicing the continuous time signals according to time sequence to obtain a continuous time azimuth history chart and a frequency domain wave beam output signal time frequency spectrum.
Compared with the traditional conventional beam forming method, the method has the advantage that the target frequency spectrum of the medium-high frequency band of 500-1000Hz in the radiation noise signal can be better recovered. Now, beam gain is performed on a certain underwater artificial target signal acquired by a 96-array element linear array, fig. 3 shows a time-frequency diagram obtained by a conventional beam forming method, the horizontal axis of the time-frequency diagram is frequency, and the vertical axis of the time-frequency diagram is time, so that it can be seen that the signal energy of the middle-high frequency part is remarkably reduced due to the fact that the low-frequency beam width is wide and the low-frequency component energy of the output beam signal is high. And fig. 4 is a target time-frequency diagram obtained by the method, it can be seen that the target radiation noise in the 300-1000Hz frequency band is better recovered, especially in the 300Hz-700Hz frequency band, the details of the target radiation noise are clearer, and meanwhile, the periodic pulse interference signal below 100Hz is also significantly suppressed.
Drawings
Fig. 1 is a working flow chart of the method for reconstructing a frequency spectrum of an output signal of a linear array beam of radiation noise based on a structured sparse underwater acoustic target according to the present invention.
FIG. 2 is a diagram of a multi-frequency multi-task sparse decomposition model according to the present invention.
Fig. 3 is a time-frequency diagram of the beam output signal obtained based on the conventional beam forming method.
FIG. 4 is a time-frequency diagram of a beam output signal obtained by the method for reconstructing a frequency spectrum of a linear array beam output signal of an underwater acoustic target radiation noise based on structured sparsity.
Detailed Description
The invention is further illustrated with reference to the following figures and examples.
The method is based on the principle of consistent spatial response of different frequency bands, and achieves the purpose of constraining the same beam width of different frequency beams by using a multi-frequency-point multi-task structured sparse algorithm model, so as to better recover the broadband radiation noise of the underwater target.
The specific work flow chart is shown in the attached figure 1. The structure diagram of the multi-frequency point multi-task sparse decomposition model is shown in figure 2.
The technical scheme adopted by the invention for solving the technical problem comprises the following detailed steps:
step 1: the linear array is assumed to have N array elements, and time domain signals acquired by each array element channel are x respectively1(n),...,xN(n), array pairingSignals x of individual channels in a column1(n),...,xN(n) performing Fast Fourier Transform (FFT) to obtain frequency domain data xF1(f),xF2(f),...,xFN(f)：
Figure BDA0002386703980000051
In the formula (1), FFT is fast Fourier transform;
step 2: the user selects a signal frequency band needing to be processed for the concerned frequency band, selects a frequency domain signal of a corresponding frequency band in the step 1, the upper frequency limit nf _ h is determined by the actual requirement of the user, and 1000Hz can be selected because the artificial target radiation noise is concentrated below 1000 Hz; if the amplitude modulation spectrum above 1000Hz needs to be paid attention to, the upper limit frequency can be set to 3000-4000Hz, and the specific numerical setting needs to be combined with the sampling frequency of the underwater acoustic equipment and the prior information judgment of a user on the target radiation noise signal; on one hand, the lower frequency limit nf _ l needs to refer to the frequency response design index of the hydrophones, is set to be 50Hz for most hydrophones, and selects a lower value for the hydrophone with better low-frequency response according to the design index, and meanwhile needs to be set by combining the actual requirements of users; recording each channel frequency domain signal in the selected formulated frequency band as x'F1(f),...,x'FN(f) Then, there are:
Figure BDA0002386703980000061
and step 3: setting an incident space azimuth vector theta with the length of ND according to the size of a linear array aperture, the resolution requirement and the requirement of an actual user, and if the azimuth range concerned by the user is 0-180 degrees, setting the vector length ND as 181 azimuth vectors at an interval of 1 degree; then constructing a space orientation dictionary D based on the selected space orientation setf，DfThe number of rows is N, the number of columns is ND, NF — nh — NF — l +1 frequency points are selected in step 2, and correspondingly, each frequency point f corresponds to one dictionary DfTotal NF dictionary, frequency point f corresponding to space azimuth dictionary DfThe expression is as follows:
Figure BDA0002386703980000062
wherein e is a natural index, i is an imaginary number symbol, pi is a circumferential rate, c is a sound velocity, p is a position vector consisting of N array elements, theta is a selected incident space azimuth angle vector, and T represents transposition;
and 4, step 4: to orientation dictionary DfAnd (3) carrying out normalization:
Figure BDA0002386703980000063
in the formula (4), sqrt represents an evolution, and diag () represents an operation of taking diagonal elements of a matrix to form a vector for a matrix operation;
and 5: modeling by adopting a structured sparse model;
for any frequency point f in the frequency band selected in the step 2, the frequency domain observation vector of each array element of the array is Yf＝[xF1(f),xF2(f),…,xFN(f)]Wherein x isFi(f) For observation of frequency point f for ith array element, YfLength is array element number N, azimuth dictionary DfIs NXND, YfThrough space orientation dictionary DfDecomposition coefficient X offThe length is the number ND of incident azimuth angles, and for each frequency point, the following sparse decomposition model is established:
Yf＝DfXf(5)
wherein the target is sparse in space, and thus
Figure BDA0002386703980000071
The medium and large coefficient elements are sparse, so that the formula (5) is YfDictionary D based on spatial orientationfDecomposing to obtain XfA sparse decomposition process of (c);
because of the spatial orientation dictionary DfIs a frequency dependent quantity and is therefore for each frequency point f within the frequency band selected in step 21～fNAll areSatisfies the following conditions:
Figure BDA0002386703980000072
because the direction of the target is determined and is independent of each frequency point value in the selected frequency band, all the frequency points selected in the formula (6) are obtained through sparse decomposition
Figure BDA0002386703980000073
In, the positions where the large coefficient elements appear should be the same or very close; in order to fully utilize the constraint information, the N sparse decomposition problems in the formula (6) are jointly solved; setting joint observation vector on all frequency points
Figure BDA0002386703980000074
X is Y channel joint space azimuth dictionary
Figure BDA0002386703980000075
D is a three-dimensional matrix, the third dimension corresponds to D at different frequency pointsfThe N sparse decomposition models of equation (6) form a multitask sparse decomposition model by joint solution:
Y＝DX (7)
jointly solving the multitask sparse decomposition model, namely adding the spatial response consistency constraint of each frequency point, wherein the spatial response consistency constraint is reflected in a sparse decomposition coefficient vector X, and a structured sparse constraint is added to the vector X large elements; obtaining a beam signal with the characteristic of constant beam width while carrying out broadband beam scanning;
since the actual signal contains noise, the following model is obtained by considering the noise n:
Y＝DX+n (8)
step 6: deducing and solving the multitask sparse decomposition model represented by the formula (8) by adopting a Bayesian variational algorithm; assuming that the noise n follows a multivariate Gaussian distribution, in equation (9)
Figure BDA0002386703980000081
Representing a multivariate gaussian distribution:
Figure BDA0002386703980000082
equation (9) assumes that the hidden variable λ of the model obeys Gamma distribution based on parameters a and b, Γ represents Gamma distribution, a and b are model-adjustable parameters, a takes a smaller value, for example, 0.01, b takes 10, and the specific value needs to be determined by continuously adjusting according to the actual signal characteristics and the actual effect:
p(λ)＝Γ(λ|a,b) (10)
equation (10) assumes the sparse decomposition coefficients to be solved as random vectors in the probabilistic model and assumes XfElement X in (1)kfAre subject to a mean of 0 and an accuracy of gammakAnd independently of each other, XfObeying multivariate Gaussian distribution, wherein the prior probability distribution is as follows:
Figure BDA0002386703980000083
for the accuracy parameter, it is still assumed that γkObey a Gamma prior distribution and are independent of each other, Gamma-Gamma1,γ2,…,γND]Is given by the following equation:
Figure BDA0002386703980000084
and 7: according to the step 6, deducing hidden variables lambda and gamma in the multitask sparse decomposition model, and finally taking the posterior probability mean value of the decomposition coefficient X as an estimated value of the decomposition coefficient X in the formula (7);
and 8: finding the position of large coefficient in X, and mapping the position to dictionary DfAzimuth of corresponding position as appearance azimuth theta of targettThen, the orientation theta of the decomposition coefficient X to the target appearance is selectedtAll frequency point results XθtObtaining an azimuth angle estimation result and a frequency domain beam amplitude of the one-time processing snapshot duration; finally repeating the continuous time signalProcessing according to the steps 1-7, and splicing according to time to obtain a continuous-time azimuth history chart and a frequency spectrum of the frequency domain wave beam output signal.
The examples are as follows:
step 1: underwater target radiation noise data x adopting one section of 96-array element linear array1(n),x2(n)…x96(n), duration 600s, sampling frequency 2048 Hz. Reading 96 array element data of the 1 st second, and carrying out fast Fourier transform on each channel signal to obtain data xf1(n),xf2(n)…xf96(n)。
Step 2: considering that the frequency response range of the hydrophone is 50Hz-100kHz, the data sampling frequency is 2048Hz, the lower limit nf _ l of the frequency band of the signal to be recovered is 50Hz, and the upper limit nf _ h is 1024 Hz. The extracted corresponding frequency band signal is x'f1(n),…,x'f96(n)。
And step 3: the array is a uniform linear array, p1＝0m,p2＝4m,...,p96380 m. The azimuth range is selected from 0 to 180 degrees with a resolution of 1 degree. According to the frequency range of 50-1024Hz set in the previous step, the frequency interval is 1Hz, and an orientation dictionary D is respectively constructed for each frequency point according to the formula (3)f。
And 4, step 4: all the orientation dictionaries D constructed in the previous step are processed according to the formula (4)fAnd (6) carrying out normalization.
And 5: constructing a joint observation vector Y ═ x 'of all frequency points in a 50-1024Hz frequency band'f1(n),x'f2(n),...,x'f96(n)]。
For each frequency point in the frequency band, the orientation dictionary D constructed on the basis of the frequency point is respectively usedfAnd constructing a sparse decomposition model, and then carrying out joint solution on the sparse decomposition models of the multiple frequency points.
Step 6: the noise n in the model follows a gaussian distribution with an accuracy (inverse variance) of λ. The prior distribution of λ is a Gamma distribution, whose two parameters a 0.01 and b 10 are initialized. The precision parameter Gamma of the sparse decomposition coefficient X in the model obeying Gaussian distribution obeys Gamma prior distribution, and two parameters c ═ a and d ═ b of the precision parameter Gamma are initialized.
And 7: according to the prior distribution hypothesis, the posterior probability inference is carried out on the lambda, the gamma and the X, and the mean value of the X is used as an estimation result.
And 8: and taking a target with the position of 124 degrees, and extracting the average value of X calculated by the positions in the current second to be used as the spectral amplitude of the beam output signal weighted by the linear array in the current second.
And step 9: the above operation is repeated for the remaining 599 seconds. And finally, splicing and accumulating the spectrum amplitude of each second according to a time axis to obtain a radiation noise time-frequency spectrogram of the target with the azimuth of 124 degrees, as shown in FIG. 4.

Claims (1)

1. A method for reconstructing the frequency spectrum of an underwater sound target radiation noise linear array wave beam output signal is characterized by comprising the following steps:
step 1: the linear array is assumed to have N array elements, and time domain signals acquired by each array element channel are x respectively1(n),…,xN(n) for each channel signal x in the array1(n),…,xN(n) performing Fast Fourier Transform (FFT) to obtain frequency domain data xF1(f),xF2(f),…,xFN(f)：
Figure FDA0002386703970000011
In the formula (1), FFT is fast Fourier transform;
step 2: the user selects a signal frequency band needing to be processed for the concerned frequency band, selects a frequency domain signal of a corresponding frequency band in the step 1, the upper frequency limit nf _ h is determined by the actual requirement of the user, and 1000Hz can be selected because the artificial target radiation noise is concentrated below 1000 Hz; if the amplitude modulation spectrum above 1000Hz needs to be paid attention to, the upper limit frequency can be set to 3000-4000Hz, and the specific numerical setting needs to be combined with the sampling frequency of the underwater acoustic equipment and the prior information judgment of a user on the target radiation noise signal; setting a lower frequency limit nf _ l by combining the actual requirements of a user; recording each channel frequency domain signal in the selected formulated frequency band as x'F1(f),...,x'FN(f) Then, there are:
Figure FDA0002386703970000012
and step 3: setting an incident space azimuth vector theta with the length ND according to the linear array aperture size, the resolution requirement and the actual user requirement, and then constructing a space azimuth dictionary D based on the selected space azimuth setf，DfThe number of rows is N, the number of columns is ND, NF — nh — NF — l +1 frequency points are selected in step 2, and correspondingly, each frequency point f corresponds to one dictionary DfTotal NF dictionary, frequency point f corresponding to space azimuth dictionary DfThe expression is as follows:
Figure FDA0002386703970000013
wherein e is a natural index, i is an imaginary number symbol, pi is a circumferential rate, c is a sound velocity, p is a position vector consisting of N array elements, theta is a selected incident space azimuth angle vector, and T represents transposition;
and 4, step 4: to orientation dictionary DfAnd (3) carrying out normalization:
Figure FDA0002386703970000021
in the formula (4), sqrt represents an evolution, and diag () represents an operation of taking diagonal elements of a matrix to form a vector for a matrix operation;
and 5: modeling by adopting a structured sparse model;
for any frequency point f in the frequency band selected in the step 2, the frequency domain observation vector of each array element of the array is Yf＝[xF1(f),xF2(f),…,xFN(f)]Wherein x isFi(f) For observation of frequency point f for ith array element, YfLength is array element number N, azimuth dictionary DfIs NXND, YfThrough space orientation dictionary DfDecomposition coefficient X offThe length is the number ND of incident azimuth angles, and for each frequency point, the following sparse decomposition model is established:
Yf＝DfXf(5)
wherein the target is sparse in space, and thus
Figure FDA0002386703970000022
The medium and large coefficient elements are sparse, so that the formula (5) is YfDictionary D based on spatial orientationfDecomposing to obtain XfA sparse decomposition process of (c);
because of the spatial orientation dictionary DfIs a frequency dependent quantity and is therefore for each frequency point f within the frequency band selected in step 21～fNAll the requirements are as follows:
Figure FDA0002386703970000023
obtained by performing sparse decomposition on all frequency points selected in the formula (6)
Figure FDA0002386703970000024
In, the positions where the large coefficient elements appear should be the same or very close; in order to fully utilize the constraint information, the N sparse decomposition problems in the formula (6) are jointly solved; setting joint observation vector on all frequency points
Figure FDA0002386703970000025
X is Y channel joint space azimuth dictionary
Figure FDA0002386703970000026
D is a three-dimensional matrix, the third dimension corresponds to D at different frequency pointsfThe N sparse decomposition models of equation (6) form a multitask sparse decomposition model by joint solution:
Y＝DX (7)
jointly solving the multitask sparse decomposition model, namely adding the spatial response consistency constraint of each frequency point, wherein the spatial response consistency constraint is reflected in a sparse decomposition coefficient vector X, and a structured sparse constraint is added to the vector X large elements; obtaining a beam signal with the characteristic of constant beam width while carrying out broadband beam scanning;
since the actual signal contains noise, the following model is obtained by considering the noise n:
Y＝DX+n (8)
step 6: deducing and solving the multitask sparse decomposition model represented by the formula (8) by adopting a Bayesian variational algorithm; assuming that the noise n follows a multivariate Gaussian distribution, in equation (9)
Figure FDA0002386703970000031
Representing a multivariate gaussian distribution:
Figure FDA0002386703970000032
equation (9) assumes that the hidden variable λ of the model obeys Gamma distribution based on parameters a and b, Γ represents Gamma distribution, a and b are model adjustable parameters, and the specific value needs to be determined by continuously adjusting according to the characteristics and the actual effect of the actual signal:
p(λ)＝Γ(λ|a,b) (10)
equation (10) assumes the sparse decomposition coefficients to be solved as random vectors in the probabilistic model and assumes XfElement X in (1)kfAre subject to a mean of 0 and an accuracy of gammakAnd independently of each other, XfObeying multivariate Gaussian distribution, wherein the prior probability distribution is as follows:
Figure FDA0002386703970000033
for the accuracy parameter, it is still assumed that γkObey a Gamma prior distribution and are independent of each other, Gamma-Gamma1,γ2,...,γND]Is given by the following equation:
Figure FDA0002386703970000034
and 7: according to the step 6, deducing hidden variables lambda and gamma in the multitask sparse decomposition model, and finally taking the posterior probability mean value of the decomposition coefficient X as an estimated value of the decomposition coefficient X in the formula (7);
and 8: finding the position of large coefficient in X, and mapping the position to dictionary DfAzimuth of corresponding position as appearance azimuth theta of targettThen, the orientation theta of the decomposition coefficient X to the target appearance is selectedtAll frequency point results XθtObtaining an azimuth angle estimation result and a frequency domain beam amplitude of the one-time processing snapshot duration; and finally, processing the continuous time signals repeatedly according to the steps 1-7, and splicing the continuous time signals according to time sequence to obtain a continuous time azimuth history chart and a frequency domain wave beam output signal time frequency spectrum.
CN202010100407.7A 2020-02-18 2020-02-18 Frequency spectrum reconstruction method for underwater sound target radiation noise linear array wave beam output signal Pending CN111273301A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010100407.7A CN111273301A (en) 2020-02-18 2020-02-18 Frequency spectrum reconstruction method for underwater sound target radiation noise linear array wave beam output signal

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010100407.7A CN111273301A (en) 2020-02-18 2020-02-18 Frequency spectrum reconstruction method for underwater sound target radiation noise linear array wave beam output signal

Publications (1)

Publication Number Publication Date
CN111273301A true CN111273301A (en) 2020-06-12

Family

ID=71002128

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010100407.7A Pending CN111273301A (en) 2020-02-18 2020-02-18 Frequency spectrum reconstruction method for underwater sound target radiation noise linear array wave beam output signal

Country Status (1)

Country Link
CN (1) CN111273301A (en)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112198515A (en) * 2020-10-13 2021-01-08 湖南国天电子科技有限公司 Parametric array shallow-section difference frequency conversion performance optimization method
CN113051739A (en) * 2021-03-17 2021-06-29 电子科技大学长三角研究院(衢州) Robustness self-adaptive processing method based on sparse constraint
CN113050025A (en) * 2021-02-02 2021-06-29 中国电子科技集团公司第二十九研究所 Method for improving direction-finding precision of millimeter wave signals without frequency information based on partition direction-finding
CN114036975A (en) * 2021-10-19 2022-02-11 中国科学院声学研究所 Target signal extraction method based on frequency domain-wavenumber domain deconvolution

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20160209538A1 (en) * 2013-08-05 2016-07-21 Schlumberger Technology Corporation Apparatus for Mode Extraction Using Multiple Frequencies
CN107703477A (en) * 2017-09-11 2018-02-16 电子科技大学 The steady broadband array signal Wave arrival direction estimating method of standard based on block management loading
US20180287822A1 (en) * 2017-03-31 2018-10-04 Mitsubishi Electric Research Laboratories, Inc. System and Method for Channel Estimation in mmWave Communications Exploiting Joint AoD-AoA Angular Spread
CN108919240A (en) * 2018-04-23 2018-11-30 东南大学 A kind of underwater acoustic target radiated noise modulation spectrum reconstruction method based on group sparsity structure
CN110208735A (en) * 2019-06-12 2019-09-06 西北工业大学 A kind of DOA Estimation in Coherent Signal method based on management loading

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20160209538A1 (en) * 2013-08-05 2016-07-21 Schlumberger Technology Corporation Apparatus for Mode Extraction Using Multiple Frequencies
US20180287822A1 (en) * 2017-03-31 2018-10-04 Mitsubishi Electric Research Laboratories, Inc. System and Method for Channel Estimation in mmWave Communications Exploiting Joint AoD-AoA Angular Spread
CN107703477A (en) * 2017-09-11 2018-02-16 电子科技大学 The steady broadband array signal Wave arrival direction estimating method of standard based on block management loading
CN108919240A (en) * 2018-04-23 2018-11-30 东南大学 A kind of underwater acoustic target radiated noise modulation spectrum reconstruction method based on group sparsity structure
CN110208735A (en) * 2019-06-12 2019-09-06 西北工业大学 A kind of DOA Estimation in Coherent Signal method based on management loading

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
LIFAN ZHAO ET AL.: "Computationally Efficient Wide-Band DOA Estimation Methods Based on Sparse Bayesian Framework", 《IEEE TRANSACTIONS ON VEHICULAR TECHNOLOGY》 *
QISONG WU ET AL.: "Multi-Task Bayesian Compressive Sensing Exploiting Intra-Task Dependency", 《IEEE SIGNAL PROCESSING LETTERS》 *
刘清宇等: "一种基于组稀疏结构的高分辨调制谱重构方法", 《中国科学:信息科学》 *
王彪等: "一种宽带水声目标DOA估计方法", 《江苏科技大学学报(自然科学版)》 *
王彪等: "一种快速稀疏贝叶斯学习的水声目标方位估计方法研究", 《声学学报》 *
陆晨翔等: "水下目标信号的结构化稀疏特征提取方法", 《哈尔滨工程大学学报》 *

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112198515A (en) * 2020-10-13 2021-01-08 湖南国天电子科技有限公司 Parametric array shallow-section difference frequency conversion performance optimization method
CN112198515B (en) * 2020-10-13 2021-06-29 湖南国天电子科技有限公司 Parametric array shallow-section difference frequency conversion performance optimization method
US11237258B1 (en) 2020-10-13 2022-02-01 Hunan Guotian Electronic Technology Co., Ltd. Method for optimization of a parametric array shallow profile difference frequency conversion performance
CN113050025A (en) * 2021-02-02 2021-06-29 中国电子科技集团公司第二十九研究所 Method for improving direction-finding precision of millimeter wave signals without frequency information based on partition direction-finding
CN113050025B (en) * 2021-02-02 2022-07-15 中国电子科技集团公司第二十九研究所 Method for improving direction-finding precision of millimeter wave signals without frequency information based on partition direction finding
CN113051739A (en) * 2021-03-17 2021-06-29 电子科技大学长三角研究院(衢州) Robustness self-adaptive processing method based on sparse constraint
CN113051739B (en) * 2021-03-17 2023-08-18 电子科技大学长三角研究院(衢州) Robustness self-adaptive processing method based on sparse constraint
CN114036975A (en) * 2021-10-19 2022-02-11 中国科学院声学研究所 Target signal extraction method based on frequency domain-wavenumber domain deconvolution

Similar Documents

Publication Publication Date Title
CN111273301A (en) Frequency spectrum reconstruction method for underwater sound target radiation noise linear array wave beam output signal
US5473759A (en) Sound analysis and resynthesis using correlograms
EP3414919B1 (en) Microphone probe, method, system and computer program product for audio signals processing
CN109683134B (en) High-resolution positioning method for rotary sound source
CN109212500A (en) A kind of miscellaneous covariance matrix high-precision estimation method of making an uproar of KA-STAP based on sparse reconstruct
CN109307855B (en) Grid-error-model-based non-grid sparse approximate minimum variance DOA estimation method
CN114782745B (en) Ocean sound velocity profile classification method and device based on machine learning
CN111175727B (en) Method for estimating orientation of broadband signal based on conditional wave number spectral density
JP7254938B2 (en) Combined source localization and separation method for acoustic sources
CN110632605B (en) Wide-tolerance large-aperture towed linear array time domain single-beam processing method
CN109658944B (en) Helicopter acoustic signal enhancement method and device
CN109061551B (en) Grid-free sparse spectrum estimation method based on polynomial root finding
CN115587291A (en) Denoising characterization method and system based on crack ultrasonic scattering matrix
CN114755628A (en) Method for estimating direction of arrival of acoustic vector sensor array under non-uniform noise
CN112505640B (en) Time-frequency analysis method for expanded B distribution pulse signal based on parameter self-adaption
CN114035157A (en) Sub-band delay estimation method and system based on expectation maximization algorithm
CN112363217A (en) Random noise suppression method and system for seismic data
CN112649787B (en) Target azimuth estimation method based on low-frequency circular array
Wu et al. Low-frequency Underwater Acoustic Signal Denoising Method in the Shallow Sea with a Low Signal-to-noise Ratio
Bayram A multichannel audio denoising formulation based on spectral sparsity
CN113820693B (en) Uniform linear array element failure calibration method based on generation of countermeasure network
CN115223580A (en) Voice enhancement method based on spherical microphone array and deep neural network
CN114061730B (en) Target scattering echo variable step length rapid self-adaptive estimation method
CN117572338A (en) Method for estimating DOA of space adjacent broadband signal based on mutual mass array
CN113921031A (en) Multichannel audio data processing method and device, computer equipment and storage medium

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
WD01 Invention patent application deemed withdrawn after publication

Application publication date: 20200612

WD01 Invention patent application deemed withdrawn after publication