CN113009418B - Target depth estimation method based on pseudo-Green function pulse delay - Google Patents

Target depth estimation method based on pseudo-Green function pulse delay Download PDF

Info

Publication number
CN113009418B
CN113009418B CN202110210711.1A CN202110210711A CN113009418B CN 113009418 B CN113009418 B CN 113009418B CN 202110210711 A CN202110210711 A CN 202110210711A CN 113009418 B CN113009418 B CN 113009418B
Authority
CN
China
Prior art keywords
acoustic array
time
depth
target
acoustic
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.)
Active
Application number
CN202110210711.1A
Other languages
Chinese (zh)
Other versions
CN113009418A (en
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.)
Institute of Acoustics CAS
Original Assignee
Institute of Acoustics CAS
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 Institute of Acoustics CAS filed Critical Institute of Acoustics CAS
Priority to CN202110210711.1A priority Critical patent/CN113009418B/en
Publication of CN113009418A publication Critical patent/CN113009418A/en
Application granted granted Critical
Publication of CN113009418B publication Critical patent/CN113009418B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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
    • G01S5/00—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations
    • G01S5/18—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations using ultrasonic, sonic, or infrasonic waves
    • G01S5/28—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations using ultrasonic, sonic, or infrasonic waves by co-ordinating position lines of different shape, e.g. hyperbolic, circular, elliptical or radial

Landscapes

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

Abstract

The invention discloses a target depth estimation method based on pseudo-Green function pulse delay, which comprises the following steps: step 1) extracting a time domain pseudo-Green function from radiation noise data of a measured target by using two approximately parallel acoustic arrays to obtain pulse time delay of the time domain pseudo-Green function; step 2) dividing the depth into a plurality of grids, calculating corresponding time domain pseudo-Green functions and corresponding pulse time delays on grid points of each depth according to the hydrological environment of the acoustic array distribution sea area, and obtaining a variation curve of the pulse time delays along with the depth; and 3) searching a grid point closest to the pulse delay of the step 1) on a change curve of the pulse delay along with the depth, wherein the corresponding depth is a depth estimation result of the measured target. The method of the invention separates the target depth estimation and the target distance estimation, and can obtain the target depth estimation result under the condition of unknown target distance and no estimation of the target distance.

Description

Target depth estimation method based on pseudo-Green function pulse delay
Technical Field
The invention belongs to the field of ocean acoustics, and particularly relates to a target depth estimation method based on pseudo-Green function pulse delay.
Background
Depth information of the target is a very important target parameter. A reliable target depth estimation result is given according to signals received by the acoustic sensor, and the method is an effective way for solving the depth resolution of the targets in water. The matching field processing method is one of methods for depth estimation of an object. The method fully utilizes the interference characteristic of the sound field, calculates the copy field vector by adopting the sound propagation model, and then carries out matching processing on the copy field and the measurement field, thereby obtaining the matching processing result of the unknown variable. Since 1976 the proposal of the concept of Matching Field (MFP) positioning, hundreds of related papers have been published at home and abroad, and a great deal of offshore experiments are carried out, thus making many progress. However, the marine environment conditions are very complex, the environmental parameters are often time-varying and space-varying, and the matching field processing method faces the difficult problem that the marine environment parameters cannot be accurately measured, so that the copy field calculation is inaccurate, the deviation between the matching processing result and the real result is too large, and even the processing result with the confidence coefficient exceeding the set threshold value cannot be given. The method requires that receiving hydrophones are arranged at a specific depth, the depth needs to be calculated according to sound field analysis frequency and marine environments such as water depth, bottom material, sound velocity profile and the like, and the depth is difficult to accurately calculate in practical application, so that the use of the method is limited. In addition, the method can utilize the dispersion characteristics of the pulse signals to carry out depth estimation on the target sound source, and the method requires that the sound source sends out the pulse signals with short time domain, and the received signals in practical application are continuous radiation noise signals, and various ships rarely radiate the transient pulse signals with short duration and high intensity.
Disclosure of Invention
Aiming at the requirement of acoustic target depth resolution, the invention provides a target depth estimation method based on pseudo-Green function pulse time delay according to the physical characteristics of a horizontal longitudinal cross-correlation function of a shallow sea low-frequency broadband sound field, which is suitable for low-frequency broadband passive target depth estimation in a shallow sea area and provides a technical means for improving the target depth resolution capability in the shallow sea area. According to the method, time domain pseudo-gray function pulse time delay is extracted by two approximately parallel acoustic arrays according to physical characteristics of a horizontal longitudinal cross-correlation function of a shallow sea low-frequency broadband sound field, and depth estimation results of targets in water are obtained by comparing the actually extracted pseudo-gray function pulse time delay with pseudo-gray function pulse time delays calculated in a simulation mode at different depths.
In order to achieve the above object, the present invention provides a target depth estimation method based on pseudo-green function pulse delay, the method comprising:
step 1) extracting a time domain pseudo-Green function from radiation noise data of a measured target by using two approximately parallel acoustic arrays to obtain pulse time delay of the time domain pseudo-Green function;
step 2) dividing the depth into a plurality of grids, calculating corresponding time domain pseudo-Green functions and corresponding pulse time delays on grid points of each depth according to the hydrological environment of the acoustic array distribution sea area, and obtaining a variation curve of the pulse time delays along with the depth;
and 3) searching a grid point closest to the pulse delay of the step 1) on a change curve of the pulse delay along with the depth, wherein the corresponding depth is a depth estimation result of the measured target.
As an improvement of the above method, the step 1) specifically includes:
step 1-1) processing multi-beat data of the acoustic array A within a period of time to obtain a broadband beam output result of the acoustic array A; obtaining the azimuth course theta of a certain target relative to the acoustic array A along the time tA(t)；
Step 1-2) processing multi-beat data of the acoustic array B within a period of time to obtain a broadband beam output result of the acoustic array B; obtaining an azimuth history theta of a target relative to the acoustic array B as a function of time tB(t)；
Step 1-3) if | θA(t)-θB(t) alpha is less than or equal to l, and alpha is a threshold value; turning to the step 1-4); otherwise, changing the time period, and turning to the step 1-1) until the output result of the broadband wave beam meets the thetaA(t)-θB(t) the appearance of a target with the value less than or equal to alpha;
step 1-4) preprocessing the output result of the two acoustic arrays in the beam forming of the target position; and performing cross-correlation operation on the preprocessed beam data, accumulating the beam data within a certain time, and calculating an average value of the beam data, namely a pseudo-Green function.
As an improvement of the above method, the step 1-1) specifically comprises:
step 1-1-1) dividing acoustic signals received by each array element of the acoustic array A into multiple beats, wherein the time length of each beat is delta T, two adjacent beats are overlapped by P%, and each beat signal is subjected to Fourier transform to obtain a frequency domain signal; the frequency domain signal of the acoustic array A is sA(f；t)＝[sA1(f；t),sA2(f；t),...,sAM(f；t)]T(ii) a f is the signal frequency, t is the starting time of the current beat data; m is the number of array elements of the acoustic array A; (.)TRepresenting a transpose;
step 1-1-2) taking the geometric center of the acoustic array A as a reference point, calculating a weighting vector w of each frequency point of the acoustic array A in a beam forming scanning direction thetaA(f,θ)＝[wA0(f,θ),wA1(f,θ),...,wAM-1(f,θ)]TWherein
Figure BDA0002952165950000021
Wherein, tauAm(θ) represents the time delay of the m +1 th array element of the acoustic array a relative to the reference point:
Figure BDA0002952165950000031
v(θ)＝-[cosθ,sinθ]T，
pAm＝[pxm,pym]T，
v (θ) is the unit vector of the incident signal, pAmThe two-dimensional coordinate of the m +1 array element relative to the reference point;
step 1-1-3) calculating a beam forming result b of the pointing direction theta of the acoustic array AA(f,θ；t)：
bA(f,θ；t)＝wA H(f,θ)sA(f；t)
Wherein, (.)HRepresents a conjugate transpose; the value range of theta is [ theta ]0-45°，θ0+45°]，θ0Representing the angle of the connecting line of the reference array elements of the acoustic array A and the acoustic array B;
steps 1-1-4) time period [ t ] for Acoustic array A1，t2]Processing the multi-beat data in the step 1-1-1) to the step 1-1-3) to obtain a broadband beam output result H of the acoustic array AA(θ,t)：
Figure BDA0002952165950000032
Wherein the symbol "Σ" denotes summation;
step 1-1-5) obtaining the azimuth course theta of a certain target relative to the acoustic array A along with the change of time t according to the broadband wave beam output result of the acoustic array AA(t)。
As an improvement of the above method, the step 1-2) specifically includes:
step 1-2-1) dividing the acoustic signals received by each array element of the acoustic array B into multiple beats, wherein the time length of each beat is delta T, two adjacent beats are overlapped by P%, and each beat is subjected to Fourier transform to obtain frequency domain signals; the frequency domain signal of the acoustic array B is sB(f；t)＝[sB1(f；t),sB2(f；t),...,sBN(f；t)]T(ii) a f is the signal frequency, selected in the range of [ fL,fH](ii) a t is the starting time of the current beat data; n is the number of array elements of the acoustic array B;
step 1-2-2) taking the geometric center of the acoustic array B as a reference point, calculating a weighting vector w of each frequency point of the acoustic array B in the beam forming scanning direction thetaB(f,θ)＝[wB0(f,θ),wB1(f,θ),...,wB,N-1(f,θ)]TWherein
Figure BDA0002952165950000033
Wherein tau isBn(θ) represents the time delay of the n +1 th array element of the acoustic array B with respect to the reference point:
Figure BDA0002952165950000041
v(θ)＝-[cosθ,sinθ]T，
pBn＝[pxn,pyn]T，
v (θ) is the unit vector of the incident signal, pBnThe two-dimensional coordinates of the (n + 1) th array element relative to the reference point;
step 1-2-3) calculating the beam forming result of the B pointing direction azimuth theta of the acoustic array
bB(f,θ；t)＝wB H(f,θ)sB(f；t)
Wherein, (.)HRepresents a conjugate transpose;
step 1-2-4) time period [ t ] for acoustic array B1，t2]Processing the multi-beat data in the step 1-2-1) to the step 1-2-3), and calculating a broadband beam output result H of the acoustic array BB(θ,t)：
Figure BDA0002952165950000042
Step 1-2-5) obtaining the azimuth course theta of the same target relative to the acoustic array B along with the change of time t according to the broadband wave beam output result of the acoustic array BB(t)。
As an improvement of the above method, the step 1-4) specifically includes:
step 1-4-1) preprocessing the beam forming output results of the two acoustic arrays in the target direction, and outputting the preprocessed beam forming output of the acoustic array A, B in the target direction
Figure BDA0002952165950000043
Respectively as follows:
Figure BDA0002952165950000044
Figure BDA0002952165950000045
wherein k is wave number, k is 2 pi f/c0(ii) a L is the horizontal spacing of the geometric centers of acoustic array A, B;
Figure BDA0002952165950000046
a horizontal orientation of the selected target relative to the acoustic array B for a starting time within the analysis time period;
IA(f) average intensity spectrum for the chosen target relative to acoustic array a:
Figure BDA0002952165950000047
IB(f) average intensity spectrum for the chosen target relative to acoustic array B:
Figure BDA0002952165950000051
wherein N istThe total number of beats; n is a radical oftThe value range of (1) is [ 600/(delta T (1-P%)), 1800/(delta T (1-P%))](ii) a Δ T is the time length of a single beat of data; p% is the overlapping rate of two adjacent beats;
step 1-4-2) performing cross-correlation operation on the preprocessed beam data, accumulating the beam data within a certain time, calculating the average value of the beam data, and obtaining a frequency domain pseudo-Green function
Figure BDA0002952165950000052
Figure BDA0002952165950000053
(·)*Represents a conjugation;
step 1-4-3) to the frequency domain pseudo-Green function
Figure BDA0002952165950000054
Performing inverse Fourier transform to obtain time domain pseudo-Green function
Figure BDA0002952165950000055
Figure BDA0002952165950000056
Wherein ifft represents inverse Fourier transform, and Re represents a real part;
step 1-4-4) time domain pseudo Green function
Figure BDA0002952165950000057
Is a pulse function, the maximum value of which envelope corresponds to the time delay point tauM：
Figure BDA0002952165950000058
Wherein argmax (·) represents the time delay corresponding to the maximum value, | · | represents the absolute value, and H (·) represents hilbert transform.
As an improvement of the above method, the step 2) specifically includes:
step 2-1) if the target distance rsIs known, then take the distance R ═ RsOtherwise, in [ R ]min，Rmax]Within the distance range, 1 distance R is selected: r ═ Rmax+Rmin)/2；RminMinimum distance, R, required to meet acoustic array far field requirementsmaxMaximum theoretical detection range of the target for the acoustic array;
step 2-2) according to [0,4H/5 ] in depth]Interval division into QzEach grid, wherein H is water depth;
step 2-3) calculating the normal wave information by adopting a sound field calculation method according to the marine environment parameters such as water depth, sound velocity profile and seabed parameters: horizontal wave number k of No. l normal wavel＝μl+iηl，μl、ηlAre each klThe real part and the imaginary part of (c); characteristic function Ψ of the ith normal wavel(z); wherein, l is 1,2, …, I; i is the total number of the normal waves;
step 2-4) calculating a time domain pseudo Green function of the q grid point on the distance R according to the normal wave information
Figure BDA0002952165950000061
Figure BDA0002952165950000062
Wherein, tau represents time delay, ifft represents inverse Fourier transform, and Re represents a real part;
Figure BDA0002952165950000063
for the frequency domain pseudo-Green function for the qth grid point on distance R:
Figure BDA0002952165950000064
wherein, B (z)1,z2,zq) Is a normalization factor:
Figure BDA0002952165950000065
z1,z2the depth of the two hydrophones is shown, the horizontal distance between the two hydrophones is L, and the horizontal distances between the two hydrophones and the sound source are r respectively1,r2；zqThe depth of a sound source corresponding to the q-th grid point on the distance R;
step 2-5) time domain pseudo Green function
Figure BDA0002952165950000066
Is a pulse function, and the time delay point corresponding to the maximum value of the envelope is marked as taur(zq)：
Figure BDA0002952165950000067
Wherein argmax (·) represents the time delay corresponding to the maximum value, | · | represents the absolute value, and H (·) represents hilbert transform.
As an improvement of the above method, the step 3) is specifically:
Figure BDA0002952165950000068
wherein argmin (·) represents the depth corresponding to the minimum value, | · | represents the absolute value, zsAnd obtaining the depth estimation result of the measured target.
The invention has the advantages that:
1. the method of the invention separates the target depth estimation from the target distance estimation, can obtain the target depth estimation result under the condition of unknown target distance and no estimation of the target distance, and overcomes the limitation that the traditional matching field processing method needs to estimate the target distance and the target depth at the same time or estimate the target distance first and then estimate the target depth;
2. the method of the invention gives the target depth estimation result according to the time domain pseudo-Green function pulse time delay, reduces the requirement on the accuracy of marine environment parameters and improves the tolerance.
Drawings
FIG. 1 is a schematic diagram of the relative positions of hydrophones and acoustic sources in accordance with the present invention;
FIG. 2 is a flow chart of a pseudo-Green function pulse delay based target depth estimation method of the present invention;
FIG. 3 is a seawater average sound velocity profile in accordance with an embodiment of the present invention; the water depth is 95m, the sea water average sound velocity profile is shown in figure 2, the sound velocity of the sea surface is 1544m/s, the sound velocity gradually decreases to 1522m/s along with the depth between the sea surface and 65m, and the depth between 65m and 95m is approximately equal to an isothermal layer; the seafloor acoustic velocity of 1608 m/s;
FIG. 4 is a graph of the relative positions of two acoustic arrays in accordance with an embodiment of the present invention; the upper right dotted line represents acoustic array a and the lower left dotted line represents acoustic array B; acoustic arrays A, B each contain 41 array elements; acoustic array A, B has a geometric center horizontal distance of 720 m; two acoustic arrays (acoustic array A, B) are deployed on the sea floor surface;
FIG. 5 is a target azimuth history map of acoustic array A in accordance with an embodiment of the present invention; the grey scale image is an output result of the acoustic array A wave beam forming, and the analysis bandwidth is 100Hz to 110 Hz; the dotted line is a change curve of the target orientation with time obtained according to the gray level map; the horizontal axis is horizontal azimuth angle, unit degree, and the vertical axis is time;
fig. 6(a) is a two-dimensional image of time domain pseudo-green function with time delay and time variation obtained in an embodiment of the present invention, where the horizontal axis is relative time delay, unit ms, and the vertical axis is absolute time;
fig. 6(b) is a time domain pseudo green function pulse delay variation curve with time, where the horizontal axis is time and the vertical axis is delay, and the unit is ms;
fig. 7(a) is a two-dimensional image of a time domain pseudo-green function changing with time delay and depth obtained by simulation calculation when the target distance is 40km, where the horizontal axis is time delay and unit ms, and the vertical axis is depth and unit m;
fig. 7(b) is a time domain pseudo-green function pulse delay variation curve with depth, where the target distance is 40km, the horizontal axis is delay, the unit ms, and the vertical axis is depth;
FIG. 8 is a graph of target depth estimation over time; the vertical axis is depth, in m; the horizontal axis is time; the maximum estimated target depth in the time period is 2.5m, and the depth measurement result conforms to the information that the target is a water surface merchant ship according to the automatic ship identification system;
FIG. 9 is a time-dependent curve of the horizontal distance of the measured target relative to the acoustic array obtained from the information of the automatic ship identification system; according to the information of the automatic identification system of the ship, the measured target can be judged to be a water surface target, and the horizontal distance of the water surface target to the acoustic array is uniformly changed from 70km to 10 km;
FIG. 10 is a time domain pseudo-Green function extracted according to the radiation noise of the underwater target, wherein the horizontal axis is time delay, unit ms, and the vertical axis is relative intensity; the time delay of the pseudo Green function in the time domain is 478.2ms according to taur(zq) Judging that the depth estimation value of the measured target is 47.5 m; the underwater target is a cooperative target with a nominal depth of 50 m.
Detailed Description
Aiming at the requirement of underwater target depth estimation, the invention provides a target depth estimation method based on pseudo-Green function pulse time delay according to the physical characteristics of the horizontal longitudinal cross-correlation function of the shallow sea low-frequency broadband sound field, which is suitable for the passive target depth estimation of the shallow sea low-frequency broadband and provides a technical means for improving the target depth resolution capability of the shallow sea.
The technical scheme adopted by the invention is as follows:
according to the theory of normal waves, the depth is zsCan be expressed as a green's function (green's function) of the generated sound field of the point sound source at horizontal distance r, depth z
Figure BDA0002952165950000081
Wherein k ismIs the horizontal wave number, k, of the No. m normal wavem＝μm(f)+jηm(f)，μm、ηmAre each kmThe real part and the imaginary part of (a) are both functions of frequency f; ΨmAnd (z) is a characteristic function of the No. m normal wave. j is an imaginary symbol; e is the base of the natural logarithm, e 2.718 … ….
As shown in FIG. 1, two dots represent hydrophones with a horizontal separation L and a depth z1、z2And the horizontal distances from the sound source (indicated by five-pointed star) are respectively r1、r2(ii) a A circular point is positioned at the central position of the two hydrophones, and the horizontal distance between the circular point and a sound source is R; the orientation of the sound source relative to the line connecting the two hydrophones is theta。
For a sound source with an orientation θ, a time-domain pseudo-green function can be obtained using the two hydrophones shown in fig. 1:
Figure BDA0002952165950000082
wherein, tau represents time delay, ift represents inverse Fourier transform, and Re represents a real part.
Figure BDA0002952165950000083
Figure BDA0002952165950000084
Wherein, B (z)1,z2,zq) Is a normalization factor; the symbol "Σ" is the summation operator.
The time domain pseudo-Green function represented by the formula (2) is an impulse function, and the time delay of the impulse function is related to parameters such as target depth, longitudinal distance and depth of hydrophones, marine environment and the like. Under the condition of known parameters such as longitudinal distance and depth of hydrophones, marine environment and the like, the depth of a measured target can be estimated according to pulse delay.
It is noted that the above discussion is for the case of two hydrophones. However, in general, in order to improve the processing gain, two approximately parallel acoustic arrays are usually used to replace the two hydrophones in practical applications, so as to obtain a spatial gain and improve the signal-to-noise ratio.
As shown in fig. 2, the present invention provides a target depth estimation method based on pseudo-green function pulse delay, which includes the following steps:
s1, processing the actual data, extracting the time domain pseudo-Green function, and obtaining the pulse time delay tauM。
S1-1, respectively forming beams of acoustic data received by two approximately parallel acoustic arrays (acoustic array A, B) arranged in water, and recording the beam forming results asbA(f, theta; t) and bB(f, θ; t). Wherein f is the signal frequency, and is selected to be in the range of [ fL,fH](ii) a t is the starting time of the current beat data; θ is the pre-beam azimuth.
Preferably, fL、fHHas a typical value range of [2c ]0/H,20c0/H]。c0Taking the reference sound velocity as the average value of the sound velocities of the water body; h is the water depth.
Preferably, the typical value range of theta is [ theta ]0-45°，θ0+45°]。θ0The angle at which the reference array elements of acoustic array a and acoustic array B are connected is shown.
S1-2, according to the result of the previous step, selecting the surface ship near the positive transverse position of the acoustic array as the sound source of the radiation noise, and the azimuth histories of the surface ship relative to the acoustic array A, B along with the change of time t are respectively marked as thetaA(t)、θB(t)。
If multiple targets are present within the analyzed time period, the targets near the lateral position of the acoustic array should be selected as much as possible.
θA(t)、θB(t) should satisfy | θA(t)-θB(t)|≤2°。
And S1-3, preprocessing the beam forming output result of the two acoustic arrays at the target position. The beamformed outputs of the pre-processed acoustic array A, B at the target azimuth are each
Figure BDA0002952165950000091
Figure BDA0002952165950000092
Wherein k is wave number, k is 2 pi f/c0(ii) a L is the horizontal spacing of the geometric centers of acoustic array A, B;
Figure BDA0002952165950000093
target-to-sound selection for start time within analysis periodHorizontal orientation of the chemical array B;
Figure BDA0002952165950000101
Figure BDA0002952165950000102
Ntthe total number of beats.
Preferably, NtTypical values of (A) are [ 600/(Delta T (1-P%)), 1800/(Delta T (1-P%))]. Δ T is the time length of a single beat of data, in units of s; p% is the overlapping rate of two adjacent beats.
S1-4, performing cross-correlation operation on the preprocessed beam data, accumulating the beam data within a certain time, calculating the average value of the beam data, and obtaining a frequency domain pseudo-Green function
Figure BDA0002952165950000103
Figure BDA0002952165950000104
(·)*Representing conjugation.
S1-5, performing inverse Fourier transform on the frequency domain pseudo Green function to obtain a time domain pseudo Green function
Figure BDA0002952165950000105
Figure BDA0002952165950000106
Wherein, tau represents time delay, ift represents inverse Fourier transform, and Re represents a real part.
Time domain pseudo green function
Figure BDA0002952165950000107
Is a pulse function whose envelope maximum corresponds to the time delayThe point is denoted as τM
Figure BDA0002952165950000108
Wherein argmax (·) represents the time delay corresponding to the maximum value, | · | represents the absolute value, and H (·) represents hilbert transform.
S2, according to the hydrological environment of the sea area where the acoustic array is distributed, simulating and calculating a time domain pseudo-Green function to obtain a variation curve tau of pulse time delay along with depthr(z)。
S2-1, in [ R ]min,Rmax]Within the distance range, 1 distance point is selected and marked as R. RminMinimum distance, R, required to meet acoustic array far field requirementsmaxThe maximum theoretical detection range of the acoustic array to the target.
Preferably, when the target distance is unknown, a typical value of R is R ═ R (R)max+Rmin) 2; if the target distance is known, then take R-RsWherein r issIs the target distance.
S2-2, dividing the grid in depth. The depth is set at [0,4H/5 ]]In the range divided into QZGrid, H is water depth.
Preferably, QZTypical values of (A) are in the range of [8H/5 lambda, 40H/5 lambda]. λ is the acoustic wavelength corresponding to the center frequency within the analysis bandwidth.
S2-3, calculating the normal wave information by adopting a sound field calculation program (such as Kraken or Krakenc) according to the water depth, the sound velocity profile, the sea bottom parameters and other marine environment parameters: horizontal wave number k of No. m normal wavem＝μm+jηm，μm、ηmAre each kmThe real part and the imaginary part of (c); characteristic function Ψ of mth normal wavem(z). Wherein M is 1,2, …, M. M is the total number of the normal waves.
And S2-4, calculating a time domain pseudo-Green function of each depth grid point on the distance R according to the normal wave information. For the qth grid point, the time-domain pseudo-Green function is
Figure BDA0002952165950000111
Wherein, tau represents time delay, ift represents inverse Fourier transform, and Re represents a real part.
Figure BDA0002952165950000112
Figure BDA0002952165950000113
Wherein, B (z)1,z2,zq) Is a normalization factor; the symbol sigma is a summation operator, and pi is a circumference ratio; z is a radical ofqThe depth of a sound source corresponding to the q-th grid point on the distance R;
Figure BDA0002952165950000114
is the target bearing obtained in step S1.
Time domain pseudo green function
Figure BDA0002952165950000115
Is a pulse function, and the time delay point corresponding to the maximum value of the envelope is marked as taur(zq)：
Figure BDA0002952165950000116
Wherein argmax (·) represents the time delay corresponding to the maximum value, | · | represents the absolute value, and H (·) represents hilbert transform.
S3, contrast τr(zq) And τMObtaining a target depth estimation result zs：
Figure BDA0002952165950000117
Wherein, argmin (·) represents the depth corresponding to the minimum value, and | · | represents the absolute value.
The invention provides a target depth estimation method based on pseudo-Green function pulse time delay according to the physical characteristics of a horizontal longitudinal cross-correlation function of a shallow sea low-frequency broadband sound field, which is suitable for low-frequency broadband passive target depth estimation in a shallow sea area and provides a technical means for improving the target depth resolution capability in the shallow sea area.
The technical solution of the present invention will now be explained by means of the accompanying drawings and specific embodiments.
A marine test was conducted in the north sea area of the south sea in 5 months of 2020. The water depth is 95m, the sea water average sound velocity profile is shown in figure 3, the sound velocity at the sea surface is 1544m/s, the sound velocity between the sea surface and the depth of 65m is gradually reduced to 1522m/s along with the depth, and the depth of 65m to 95m is approximately equal to an isothermal layer. The seafloor acoustic velocity is 1608 m/s. Two acoustic arrays (acoustic array A, B) are deployed on the sea floor surface in a matrix arrangement as shown in fig. 4.
Data were analyzed over a period of 0 to 5 points, 5/19/2020. During this time period, the data is divided into multiple beats of data of 10s per beat, with adjacent beats overlapping by 40%. Beamforming output H of acoustic array AA(θ, t) As shown in FIG. 5, there is a target at 47 ° azimuth, and the time-dependent change of azimuth is shown by the dotted line in FIG. 5. Respectively obtaining the beam output of the acoustic array A, B in the target direction and carrying out normalization processing to obtain
Figure BDA0002952165950000121
And
Figure BDA0002952165950000122
performing cross-correlation operation in the frequency domain, and averaging in a certain time period to obtain a frequency domain pseudo-Green function
Figure BDA0002952165950000123
Then, inverse Fourier transform is carried out on the frequency domain pseudo Green function to obtain a time domain pseudo Green function
Figure BDA0002952165950000124
In this example, NtTake 100 corresponding to a time period length of10 minutes; the frequency analysis range is 100Hz to 150 Hz. FIG. 6(a) shows a two-dimensional image of a time domain pseudo-Green function as a function of time delay and time, with time delay on the horizontal axis, in ms, and time on the vertical axis; fig. 6(b) is a time-domain pseudo-green function pulse delay variation curve with time, where the horizontal axis is time and the vertical axis is delay in ms.
The target distance is taken as R-40 km, at which the depth is divided into 36 grids in the range of [2.5m,80m ], step size 2.5 m. According to the sea area environment parameters shown in fig. 3, the Kraken program is adopted to calculate the normal wave information, and the time domain pseudo-green function of each grid point is calculated according to the normal wave information. For the qth grid point, its time-domain pseudo-green function is:
Figure BDA0002952165950000125
wherein, tau represents time delay, ift represents inverse Fourier transform, and Re represents a real part.
Figure BDA0002952165950000126
Figure BDA0002952165950000127
Wherein B (z)1,z2,zq) Is a normalization factor; z is a radical ofqThe sound source depth corresponding to the qth grid point;
Figure BDA0002952165950000128
is the target orientation.
In this embodiment, acoustic array A, B depth z1＝z295m, 720m horizontal spacing L of the geometric centers of the acoustic array A, B, target orientation
Figure BDA0002952165950000129
The two-dimensional image of the time domain pseudo-Green function obtained by simulation calculation along with the time delay and depth change is shown in FIG. 7(a), and the horizontal axis isDelay, in ms, and the vertical axis, in depth, in m. Fig. 7(b) is a time domain pseudo-green function pulse delay variation curve with depth, where the horizontal axis is delay, unit ms, and the vertical axis is depth.
Contrast τr(zq) And τMObtaining a target depth estimation result of zs：
Figure BDA0002952165950000131
Wherein, argmin (·) represents the depth corresponding to the minimum value, and | · | represents the absolute value.
According to the above equation, a time-dependent variation curve of the target depth estimation result is obtained, as shown in fig. 8. The vertical axis is depth, in m; the horizontal axis is time. The estimated maximum target depth in the time period is 2.5m, and the depth measurement result accords with the information that the target is the water surface merchant ship according to the automatic ship identification system. According to the information of the automatic identification system of the ship, the change curve of the horizontal distance of the water surface target relative to the acoustic array along with the time is shown in fig. 9, and the horizontal distance of the water surface target relative to the acoustic array is uniformly changed from 70km to 10km in the time period corresponding to the embodiment.
In addition, by adopting the method, wide acoustic signals transmitted by the artificial hanging sound source received by the acoustic array are processed, and the data processing frequency band is 100Hz-150 Hz. Fig. 10 is a time domain pseudo-green function extracted from an acoustic signal thereof, with time delay on the horizontal axis, unit ms, and relative intensity on the vertical axis. The time delay of the pseudo Green function in the time domain is 478.2ms according to taur(zq) The estimated depth value of the target to be measured is determined to be 47.5m, which is consistent with the actual hanging depth (50m) of the sound source.
The above results verify the validity of the method.
Finally, it should be noted that the above embodiments are only used for illustrating the technical solutions of the present invention and are not limited. Although the present invention has been described in detail with reference to the embodiments, it will be understood by those skilled in the art that various changes may be made and equivalents may be substituted without departing from the spirit and scope of the invention as defined in the appended claims.

Claims (6)

1. A pseudo-Green function pulse delay-based target depth estimation method, the method comprising:
step 1) extracting a time domain pseudo-Green function from radiation noise data of a measured target by using two approximately parallel acoustic arrays to obtain pulse time delay of the time domain pseudo-Green function;
step 2) dividing the depth into a plurality of grids, calculating corresponding time domain pseudo-Green functions and corresponding pulse time delays on grid points of each depth according to the hydrological environment of the acoustic array distribution sea area, and obtaining a variation curve of the pulse time delays along with the depth;
step 3) searching a grid point closest to the pulse delay of the step 1) on a change curve of the pulse delay along with the depth, wherein the corresponding depth is a depth estimation result of the measured target;
the step 1) specifically comprises the following steps:
step 1-1) processing multi-beat data of the acoustic array A within a period of time to obtain a broadband beam output result of the acoustic array A; obtaining the azimuth course theta of a certain target relative to the acoustic array A along the time tA(t)；
Step 1-2) processing multi-beat data of the acoustic array B within a period of time to obtain a broadband beam output result of the acoustic array B; obtaining an azimuth history theta of a target relative to the acoustic array B as a function of time tB(t)；
Step 1-3) if | θA(t)-θB(t) alpha is less than or equal to l, and alpha is a threshold value; turning to the step 1-4); otherwise, changing the time period, and turning to the step 1-1) until the output result of the broadband wave beam meets the thetaA(t)-θB(t) the appearance of a target with the value less than or equal to alpha;
step 1-4) preprocessing the output result of the two acoustic arrays in the beam forming of the target position; and performing cross-correlation operation on the preprocessed beam data, accumulating the beam data within a certain time, and calculating an average value of the beam data, namely a pseudo-Green function.
2. The pseudo-green function pulse delay-based target depth estimation method according to claim 1, wherein the step 1-1) specifically comprises:
step 1-1-1) dividing acoustic signals received by each array element of the acoustic array A into multiple beats, wherein the time length of each beat is delta T, two adjacent beats are overlapped by P%, and each beat signal is subjected to Fourier transform to obtain a frequency domain signal; the frequency domain signal of the acoustic array A is sA(f；t)＝[sA1(f；t),sA2(f；t),...,sAM(f；t)]T(ii) a f is the signal frequency, t is the starting time of the current beat data; m is the number of array elements of the acoustic array A; (.)TRepresenting a transpose;
step 1-1-2) taking the geometric center of the acoustic array A as a reference point, calculating a weighting vector w of each frequency point of the acoustic array A in a beam forming scanning direction thetaA(f,θ)＝[wA0(f,θ),wA1(f,θ),...,wAM-1(f,θ)]TWherein
Figure FDA0003209182080000021
Wherein, tauAm(θ) represents the time delay of the m +1 th array element of the acoustic array a relative to the reference point:
Figure FDA0003209182080000022
v(θ)＝-[cosθ,sinθ]T，
pAm＝[pxm,pym]T，
v (θ) is the unit vector of the incident signal, pAmThe two-dimensional coordinate of the m +1 th array element relative to a reference point, and c0 is a reference sound speed;
step 1-1-3) calculating a beam forming result b of the pointing direction theta of the acoustic array AA(f,θ；t)：
bA(f,θ；t)＝wA H(f,θ)sA(f；t)
Wherein, (.)HRepresents a conjugate transpose; the value range of theta is [ theta ]0-45°，θ0+45°]，θ0Representing the angle of the connecting line of the reference array elements of the acoustic array A and the acoustic array B;
steps 1-1-4) time period [ t ] for Acoustic array A1，t2]Processing the multi-beat data in the step 1-1-1) to the step 1-1-3) to obtain a broadband beam output result H of the acoustic array AA(θ,t)：
Figure FDA0003209182080000023
Wherein the symbol "Σ" denotes summation;
step 1-1-5) obtaining the azimuth course theta of a certain target relative to the acoustic array A along with the change of time t according to the broadband wave beam output result of the acoustic array AA(t)。
3. The pseudo-green function pulse delay-based target depth estimation method according to claim 2, wherein the step 1-2) specifically comprises:
step 1-2-1) dividing the acoustic signals received by each array element of the acoustic array B into multiple beats, wherein the time length of each beat is delta T, two adjacent beats are overlapped by P%, and each beat is subjected to Fourier transform to obtain frequency domain signals; the frequency domain signal of the acoustic array B is sB(f；t)＝[sB1(f；t),sB2(f；t),...,sBN(f；t)]T(ii) a f is the signal frequency, selected in the range of [ fL,fH](ii) a t is the starting time of the current beat data; n is the number of array elements of the acoustic array B;
step 1-2-2) taking the geometric center of the acoustic array B as a reference point, calculating a weighting vector w of each frequency point of the acoustic array B in the beam forming scanning direction thetaB(f,θ)＝[wB0(f,θ),wB1(f,θ),...,wB,N-1(f,θ)]TWherein
Figure FDA0003209182080000031
Wherein tau isBn(θ) represents the time delay of the n +1 th array element of the acoustic array B with respect to the reference point:
Figure FDA0003209182080000032
v(θ)＝-[cosθ,sinθ]T，
pBn＝[pxn,pyn]T，
v (θ) is the unit vector of the incident signal, pBnThe two-dimensional coordinates of the (n + 1) th array element relative to the reference point;
step 1-2-3) calculating the beam forming result of the B pointing direction azimuth theta of the acoustic array
bB(f,θ；t)＝wB H(f,θ)sB(f；t)
Wherein, (.)HRepresents a conjugate transpose;
step 1-2-4) time period [ t ] for acoustic array B1，t2]Processing the multi-beat data in the step 1-2-1) to the step 1-2-3), and calculating a broadband beam output result H of the acoustic array BB(θ,t)：
Figure FDA0003209182080000033
Step 1-2-5) obtaining the azimuth course theta of the same target relative to the acoustic array B along with the change of time t according to the broadband wave beam output result of the acoustic array BB(t)。
4. The pseudo-green function pulse delay-based target depth estimation method according to claim 3, wherein the steps 1-4) specifically comprise:
step 1-4-1) preprocessing the beam forming output results of the two acoustic arrays in the target direction, wherein the acoustic arrays A, B are in the target after preprocessingAzimuthal beamforming output
Figure FDA0003209182080000034
Respectively as follows:
Figure FDA0003209182080000035
Figure FDA0003209182080000036
wherein k is wave number, k is 2 pi f/c0(ii) a L is the horizontal spacing of the geometric centers of acoustic array A, B;
Figure FDA0003209182080000037
a horizontal orientation of the selected target relative to the acoustic array B for a starting time within the analysis time period;
IA(f) average intensity spectrum for the chosen target relative to acoustic array a:
Figure FDA0003209182080000041
IB(f) average intensity spectrum for the chosen target relative to acoustic array B:
Figure FDA0003209182080000042
wherein N istThe total number of beats; n is a radical oftThe value range of (1) is [ 600/(delta T (1-P%)), 1800/(delta T (1-P%))](ii) a Δ T is the time length of a single beat of data; p% is the overlapping rate of two adjacent beats;
step 1-4-2) performing cross-correlation operation on the preprocessed beam data, accumulating the beam data within a certain time, calculating the average value of the beam data, and obtaining a frequency domain pseudo-Green function
Figure FDA0003209182080000043
Figure FDA0003209182080000044
(·)*Represents a conjugation;
step 1-4-3) to the frequency domain pseudo-Green function
Figure FDA0003209182080000045
Performing inverse Fourier transform to obtain time domain pseudo-Green function
Figure FDA0003209182080000046
Figure FDA0003209182080000047
Wherein ifft represents inverse Fourier transform, and Re represents a real part;
step 1-4-4) time domain pseudo Green function
Figure FDA0003209182080000048
Is a pulse function, the maximum value of which envelope corresponds to the time delay point tauM：
Figure FDA0003209182080000049
Wherein argmax (·) represents the time delay corresponding to the maximum value, | · | represents the absolute value, and H (·) represents hilbert transform.
5. The pseudo-green function pulse delay-based target depth estimation method according to claim 4, wherein the step 2) specifically comprises:
step 2-1) if the target distance rsIs known, then take the distanceR＝rsOtherwise, in [ R ]min，Rmax]Within the distance range, 1 distance R is selected: r ═ Rmax+Rmin)/2；RminMinimum distance, R, required to meet acoustic array far field requirementsmaxMaximum theoretical detection range of the target for the acoustic array;
step 2-2) according to [0,4H/5 ] in depth]Interval division into QzEach grid, wherein H is water depth;
step 2-3) calculating the normal wave information by adopting a sound field calculation method according to the marine environment parameters such as water depth, sound velocity profile and seabed parameters: horizontal wave number k of No. l normal wavel＝μl+iηl，μl、ηlAre each klThe real part and the imaginary part of (c); characteristic function Ψ of the ith normal wavel(z); wherein, l is 1,2, …, I; i is the total number of the normal waves;
step 2-4) calculating a time domain pseudo Green function of the q grid point on the distance R according to the normal wave information
Figure FDA0003209182080000051
Figure FDA0003209182080000052
Wherein, tau represents time delay, ifft represents inverse Fourier transform, and Re represents a real part;
Figure FDA0003209182080000053
for the frequency domain pseudo-Green function for the qth grid point on distance R:
Figure FDA0003209182080000054
wherein, B (z)1,z2,zq) Is a normalization factor:
Figure FDA0003209182080000055
z1,z2the depth of the two hydrophones is shown, the horizontal distance between the two hydrophones is L, and the horizontal distances between the two hydrophones and the sound source are r respectively1,r2；zqThe depth of a sound source corresponding to the q-th grid point on the distance R;
step 2-5) time domain pseudo Green function
Figure FDA0003209182080000056
Is a pulse function, and the time delay point corresponding to the maximum value of the envelope is marked as taur(zq)：
Figure FDA0003209182080000057
Wherein argmax (·) represents the time delay corresponding to the maximum value, | · | represents the absolute value, and H (·) represents hilbert transform.
6. The pseudo-green function pulse delay-based target depth estimation method according to claim 5, wherein the step 3) is specifically:
Figure FDA0003209182080000061
wherein argmin (·) represents the depth corresponding to the minimum value, | · | represents the absolute value, zsAnd obtaining the depth estimation result of the measured target.
CN202110210711.1A 2021-02-25 2021-02-25 Target depth estimation method based on pseudo-Green function pulse delay Active CN113009418B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110210711.1A CN113009418B (en) 2021-02-25 2021-02-25 Target depth estimation method based on pseudo-Green function pulse delay

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110210711.1A CN113009418B (en) 2021-02-25 2021-02-25 Target depth estimation method based on pseudo-Green function pulse delay

Publications (2)

Publication Number Publication Date
CN113009418A CN113009418A (en) 2021-06-22
CN113009418B true CN113009418B (en) 2021-11-09

Family

ID=76386162

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110210711.1A Active CN113009418B (en) 2021-02-25 2021-02-25 Target depth estimation method based on pseudo-Green function pulse delay

Country Status (1)

Country Link
CN (1) CN113009418B (en)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103076590A (en) * 2012-12-31 2013-05-01 东南大学 Method for positioning underwater sound pulse signal on basis of frequency estimation
CN103076594A (en) * 2012-12-31 2013-05-01 东南大学 Method for positioning underwater sound pulse signal by double array elements on basis of cross-correlation
CN103176166A (en) * 2013-03-14 2013-06-26 东南大学 Tracking algorithm for time difference of arrival of signals for acoustic passive positioning
CN103197318A (en) * 2013-03-18 2013-07-10 中国科学院声学研究所 Time delay estimation method based on the Pattern delay coding underwater acoustic positioning

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106595834B (en) * 2016-11-10 2019-01-04 西北工业大学 A method of obtaining the horizontal longitudinal correlation of the big depth sound field in deep-sea

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103076590A (en) * 2012-12-31 2013-05-01 东南大学 Method for positioning underwater sound pulse signal on basis of frequency estimation
CN103076594A (en) * 2012-12-31 2013-05-01 东南大学 Method for positioning underwater sound pulse signal by double array elements on basis of cross-correlation
CN103176166A (en) * 2013-03-14 2013-06-26 东南大学 Tracking algorithm for time difference of arrival of signals for acoustic passive positioning
CN103197318A (en) * 2013-03-18 2013-07-10 中国科学院声学研究所 Time delay estimation method based on the Pattern delay coding underwater acoustic positioning

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
利用矢量海洋环境噪声提取声场格林函数;周建波 等;《声学学报》;20190531;第44卷(第3期);第337-344页 *
基于快速广义互相关时延估计算法的深度测量技术;黄清泉;《四川兵工学报》;20091231;第30卷(第12期);第89-91页 *

Also Published As

Publication number Publication date
CN113009418A (en) 2021-06-22

Similar Documents

Publication Publication Date Title
CN108226933B (en) Deep sea broadband target depth estimation method based on fringe interference structure
CN113011006B (en) Target depth estimation method based on cross-correlation function pulse waveform matching
CN108828522B (en) Underwater target radiation noise measurement method formed by utilizing vertical array LCMV wave beams
EP2263097B1 (en) Autonomous sonar system and method
CN103777177B (en) A kind of ultra-short baseline submarine target localization method based on broadband signal time delay detection
CN107272005B (en) Active positioning method based on target echo arrival time delay and arrival angle under reliable acoustic path
US7307914B1 (en) Hypothesized range and depth sonar processing method
CN107179535A (en) A kind of fidelity based on distortion towed array strengthens the method for Wave beam forming
CN103076594B (en) Method for positioning underwater sound pulse signal by double array elements on basis of cross-correlation
CN109100711B (en) Single-base active sonar low-computation-quantity three-dimensional positioning method in deep sea environment
CN103438987A (en) Ship radiation noise source distinguishing method based on super directivity small-bore cylindrical array
CN108845325A (en) Towed linear-array sonar submatrix error misfits estimation method
CN109061654B (en) Single-circular-ring-array active three-dimensional positioning method in deep sea environment
CN112269164A (en) Weak target positioning method based on interference structure matching processing under deep sea reliable acoustic path
CN101915922A (en) Towed linear array passive ranging method
CN105005026A (en) Near-field target sound source three-dimensional passive positioning method
CN108107436A (en) A kind of submarine target based on reliable acoustic path is actively classified and localization method
CN104793212A (en) Method for active-sonar remote detection by means of sound wave sub-bottom reflection
CN112098938B (en) Six-element cone vector array-based underwater acoustic target dimension reduction matching sound field positioning method
CN104714235A (en) Ranging method and system for double low-frequency vector hydrophone arrays
CN108398690B (en) Submarine backscattering intensity measuring method
Ma et al. Spatiotemporal two-dimensional deconvolution beam imaging technology
CN113009418B (en) Target depth estimation method based on pseudo-Green function pulse delay
CN115656994B (en) Real-time calibration method for bistatic active detection towing array shape
CN113009419B (en) Target depth estimation method based on frequency domain cross-correlation matching

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
GR01 Patent grant
GR01 Patent grant