CN106019290A - Weighted broadband time reversal operator decomposition multi-target acoustic imaging method - Google Patents

Abstract

A weighted broadband time reversal operator decomposition multi-target acoustic imaging method has the core that: carrying out Fourier decomposition on signals received by the receiving and transmitting combined array, dividing the signals into sub-bands through a narrow-band filter, decomposing time reversal operator singular values corresponding to the sub-bands, judging the number of signal sub-spaces according to the change trend of the characteristic values, correlating characteristic vectors corresponding to all the signal sub-spaces with transfer vectors based on an acoustic propagation model to obtain a ambiguity function, extracting the maximum value of each sub-band ambiguity function and a corresponding two-dimensional sound field space, designing a weighting coefficient according to the characteristic values, weighting the maximum value of the ambiguity function, coherently accumulating the weighting results of all the sub-bands, coherently accumulating the results of all the signal sub-spaces, and presenting the results as a space image.

Description

During weighting broadband, inverse operator decomposes multiple target acoustic imaging method

Technical field

The invention belongs to the acoustic detection field of marine resources development and utilization, specifically one uses array to multiple The method that point target or autgmentability target carry out acoustic imaging.

Background technology

Owing to sound wave is better than light wave and electromagnetic wave in the performance of water transmission, so, it is applied to the detection base of marine resources It is based on acoustic means on Ben.Current practice remains in the spatial processing technique of sonar system with last century five, six The ten's rise and grow up with the plane wave Wave beam forming under white noise background and free field propagation conditions.But, Under neritic environment, traditional processing method based on plane wave model is the most applicable, due to acoustical signal propagation in ocean in Existing multipath/many Normal mode analysis state, it is easy to cause the pseudo-target during detection or diplopia phenomenon, affect the inspection of sonar system Survey performance.Time inverse processing start to be applied to underwater sound field from last century end, be quite to have the active target of application prospect to detect at present Means, time inverse operator decompose be by analyze target scattering matrix reach focus emission, location and the one side of imageable target Method.

But, time current, inverse operator decomposition method is applied to the multi-target imaging in marine exploration field, the side generally used Method, one assumes that the corresponding different characteristic vector of different targets, by relevant to transfer matrix for target characteristic of correspondence vector, two It is the MUSIC algorithm utilizing noise subspace, there is following limitation: the process of (1) broadband signal only considers mid frequency Or the mode of the result noncoherent accumulation of whole frequency range, does not effectively utilize the information of broadband signal；(2) multiple mesh Mark needs multiple image to show, it is impossible to represent the relative position relation of multiple target intuitively；(3) scattering of target itself is special Property, the comprehensive function such as relation of waveguide propagation effect and transducer and frequency, be eventually exhibited as target scattering characteristics and frequency Relation function, and at present processing procedure do not consider and utilize this scattering properties.

Summary of the invention

The disadvantages mentioned above of present invention prior art to be overcome, it is provided that during weighting broadband, inverse operator decomposes multiple target acoustic imaging side Method, the method on the basis of inverse operator decomposition method, effectively utilizes target scattering characteristics and wide-band-message when existing, from And reach eliminate pseudo-target, effectively display autgmentability spatial characters and show the purpose of multiple target simultaneously.

During the weighting broadband that the present invention proposes, inverse operator decomposes multiple target acoustic imaging method, according to the variation tendency of eigenvalue Determine signal subspace, extract each signal subspace, the maximum of each subband ambiguity function and corresponding two dimension sound Space, field, according to the functional relationship design weight coefficient of eigenvalue of maximum with frequency, the corresponding signal subspace of weighting and subband Ambiguity function, and by the ambiguity function noncoherent accumulation after weighting, final result is shown as sound field images, determines simultaneously The locus of multiple targets.

Below the inventive method is described further.

The parameters,acoustic such as depth of water of environment, Sound speed profile, density, sediment parameters etc. are had in the case of enough understanding, Complete three-dimensional acoustics image and the identification of Target space position of target as follows:

(1) coordinate system is set up；If the transmitting-receiving of P transducer composition is closed puts linear array, it is vertically disposed in water, with linear array as z Axle, horizontal direction is that r axle sets up coordinate system, and the water surface is zero, and first transducer degree of depth away from the water surface is z₁, P The transducer degree of depth away from the water surface is z_P；

(2) transducer array is received data snap and be arranged in column vector, according to predetermined sampling frequency to each data Snap is sampled, y_i(n)=[y_1i(n)…y_pi(n)…y_Pi(n)]^H, p=1 ..., P, i=1 ..., M, wherein n represents the moment, H represents that transposition, M represent the fast umber of beats of data；

(3) signal after sampling being carried out Short Time Fourier Transform, the n moment is represented by

$y_i\left(\mathrm{\ω}\right)=\left[\begin{array}{c}{\mathrm{\Σ}}_{l=0}^L{y}_{1i}(n-l)w\left(l\right)e^{-j\mathrm{\ω}(n-l)}\\ \·\\ \·\\ \·\\ {\mathrm{\Σ}}_{l=1}^L{y}_{Pi}(n-l)w\left(l\right)e^{-j\mathrm{\ω}(n-l)}\end{array}\right]---\left(1\right)$

Wherein w (l) is window function sequence, and L is the length of window function；

(4) all of M data snap is arranged in matrix Time then, inverse operator is K (ω)=Y (ω) Y^H(ω)；

(5) by time inverse operator be divided into Q subband, carry out singular value decomposition:

K(ω_q)=U (ω_q)Λ(ω_q)V*(ω_q), q=1 ..., Q (3)

(6) function curve of the eigenvalue after each subband noncoherent accumulation and horizontal range r is exported, according to eigenvalue of maximum Curve judge to comprise time inverse operator K corresponding to time window of target_o(ω_q), corresponding eigenvalue Λ_o(ω_q) and special Levy vector U_o(ω_q) and V_o*(ω_q)；

(7) threshold value is set, according to eigenvalue Λ_o(ω_q), q=1 ..., the variation tendency of Q judges signal subspace, it is assumed that The number of signal subspace is T；

(8) according to eigenvalue Λ_o(ω) design weight coefficient, finds out the eigenvalue λ of signal subspace_t(ω_q), q= 1 ..., the maximum λ that Q ties up in frequency^t _max, t=1 ... T, the eigenvalue of remaining each subband and maximum λ^t _maxRatio be and add Weight coefficient:

η_t(ω_q)=λ_t(ω_q)/λ^t _max, q=1 ..., Q, t=1 ... T (4)

(9) be grid by Spacial domain decomposition interested, the intersection point of grid be suppose target place position (r, Z), wherein r represents that supposition target closes the horizontal range putting battle array away from transmitting-receiving, and z represents the depth of water supposing target；

(10) Underwater Acoustic Environment faced according to formation method, determines be suitable for propagation model, respectively obtains supposition mesh Mark and transmitting-receiving are closed and are put transmission vector g (r, z, the ω between battle array_q)=[g₁(z₁,r,z,ω_q)…g_j(z_j,r,z,ω_q)…g_P(z_P, r,z,ω_q)]^T, wherein g_p(z_p,r,z,ω_q) represent the transmission function supposed between target and pth transducer.

(11) each signal subspace, the ambiguity function of each subband is:

I_t(r,z,ω_q)=| g^H(r,z,ω_q)u_t(ω_q)|², q=1 ..., Q, t=1 ..., T (5)

Wherein u_t(ω_q) represent the t signal subspace q-th subband feature value λ_t(ω_q) characteristic of correspondence vector, H table Show conjugate transpose；

(12) all signal subspaces are extracted, the maximum I of each subband ambiguity function^t _max(r,z,ω_q) and corresponding Sound field locus (r^t _max(ω_q), z^t _max(ω_q))；

(13) by the sound field locus (r of all subbands^t _max(ω_q), z^t _max(ω_q)) arrange be a set omega (r^t _max (ω_q), z^t _max(ω_q)), and ambiguity in definition degree function again:

$\begin{array}{c}{I}_t(r,z,{\mathrm{\ω}}_q)=\\ \left\{\begin{array}{cc}0,& (r,z)\∉\mathrm{\Ω}\left(r_{\max }^t\right({\mathrm{\ω}}_q),z_{\max }^t({\mathrm{\ω}}_q\left)\right)\\ I_t(r,z,{\mathrm{\ω}}_q),& (r,z)\∈\mathrm{\Ω}\left(r_{\max }^t\right({\mathrm{\ω}}_q),z_{\max }^t({\mathrm{\ω}}_q\left)\right)\end{array}\right.\end{array}---\left(6\right)$

(14) to all subbands, the ambiguity function of all signal subspaces is weighted adding up:

$I_{end}(r,z)={\mathrm{\Σ}}_{t=1}^T{\mathrm{\Σ}}_{q=1}^Q{\mathrm{\η}}_t\left({\mathrm{\ω}}_q\right)I_t(r,z,{\mathrm{\ω}}_q)---\left(7\right)$

The ambiguity function finally obtained is shown as the 3-D view of distance r and degree of depth z, multiple target can be determined simultaneously Corresponding two-dimensional space region.

The invention have the advantage that without judging the corresponding relation between each target and eigenvalue and characteristic vector, fully Utilize signal subspace information, multiple targets are presented in same piece image, clearly demonstrate the phase para-position between target Put.

Accompanying drawing explanation

Experimental arrangement and the coordinate of Fig. 1 present invention arrange figure.

Fig. 2 inventive feature value and the function relation figure of distance.

The eigenvalue of Fig. 3 target of the present invention place time window and the function relation figure of frequency.

Fig. 4 is the image of two Extended target of the present invention.

Specific embodiments

Below in conjunction with the accompanying drawings, the present invention is further described by being embodied as example.It is 1.44m a degree of depth Waveguide laboratory water tank in, three of pond patch sound eliminating tiles, the bottom in pond has spread the sand that one layer of 0.22m is thick.The ring in pond Border parameter is as follows: in water, the velocity of sound is constant, learns as c by measuring water temperature calculating₁=1493m/s, water body density p₁= 1000kg/m3, sediment parameters density p₂=1800kg/m3, the velocity of sound is c₂=1650m/s, attenuation quotient α₂=0.67dB/ λ, The density of substrate is ρ₃=1800kg/m3, velocity of sound c₃=1580m/s, attenuation quotient α₃=0.8dB/ λ.

(1), in the present embodiment, the foundation of experimental arrangement and coordinate system is as it is shown in figure 1, transmitting-receiving conjunction puts 32 array elements of battle array, entirely Structuring the formation in field, first array element, away from water surface 0.04m, array element distance 0.04m, is vertically disposed in pond；Two a diameter of 0.21m, The column type target of a length of 0.51m, a basin bottom being placed on matrix 6m, one is placed on matrix 8m, and the degree of depth is At 0.83m.Launch signal be pulse width be the mid frequency of 0.5ms be 18kHz, with the linear FM signal of a width of 5kHz.

(2) transducer array is received data snap and be arranged in column vector, be that 50kHz is to every number according to sample frequency Sample according to snap, y_i(n)=[y_1i(n)…y_pi(n)…y_Pi(n)]^H, p=1 ..., 32, i=1 ..., 32；

(3) the window length of the short time-window selected is 8 times of transmitted pulse width, i.e. count be 200 rectangular window, to sampling after letter Number carry out Short Time Fourier Transform

(4) all of 32 data snaps are arranged in matrixThen Time inverse operator be K (ω)=Y (ω) Y^H(ω)；

(5) the time inverse operator of the bandwidth of 5kHz is divided into Q=501 sub-band and carries out singular value decomposition:

K (ω)=U (ω) Λ (ω) V* (ω) (3)

(6) function curve of the eigenvalue after each subband noncoherent accumulation and horizontal range r is exported, it is judged that comprise target Time inverse operator K corresponding to time window_o(ω), corresponding eigenvalue Λ_o(ω) and characteristic vector U_o(ω) and V_o*(ω)； Fig. 2 is the function relation figure of front 5 eigenvalues and distance.As shown in Figure 2, eigenvalue You Liangge local peaking, lay respectively at two The distance at individual target place, may thereby determine that the time window containing two targets.

(7) in Fig. 3, left side figure is Λ_o(ω_q) with the variation tendency of frequency, it can be determined that signal subspace T=4；

(8) design weight coefficient, Fig. 3 top right plot is the weight coefficient of subspace corresponding to eigenvalue of maximum, and bottom-right graph is The weight coefficient of subspace corresponding to Second Eigenvalue；

η_t(ω_q)=λ_t(ω_q)/λ^t _max, q=1 ..., Q, t=1 ... T (4)

(9) being grid by Spacial domain decomposition interested, level initiates distance for 0.1m, and step-size in search is 0.05m, cuts Only distance is 15m, and vertical initial depth is 0m, and step-size in search is 0.0075m, and the cut-off degree of depth is 1.5m.The intersection point of grid is (r, z), wherein r represents and supposes that target closes the horizontal range putting battle array away from transmitting-receiving, and z represents supposition target assuming that the position at target place The depth of water；

(10) Underwater Acoustic Environment faced according to formation method, determines be suitable for propagation model, respectively obtains supposition mesh Mark and transmitting-receiving are closed and are put transmission vector g (r, z, the ω between battle array_q)=[g₁(z₁,r,z,ω_q)…g_p(z_p,r,z,ω_q)…g_P(z_P, r,z,ω_q)]^T, wherein g_p(z_p,r,z,ω_q) represent the transmission function supposed between target and pth transducer, this embodiment party Case uses Normal mode analysis propagation model, and its transmission vector is:

$g_p(z_p,r,z,{\mathrm{\ω}}_q)=-{\mathrm{j\πe}}^{{\mathrm{j\ω}}_{q}t}{\mathrm{\Σ}}_l{Z}_l\left(z_p\right)Z_l\left(z\right)H_0^{\left(2\right)}\left({\mathrm{\κ}}_{l}r\right)---\left(5\right)$

Wherein Z_lRepresent l propagating mode characteristic of correspondence function, H₀ ⁽²⁾() is Hankel function, κ_lRepresent wave number；

(11) each signal subspace, the ambiguity function of each subband is:

I_t(r,z,ω_q)=| g^H(r,z,ω_q)u_t(ω_q)|², q=1 ..., Q, t=1 ..., T (6)

Wherein u_t(ω_q) represent the t signal subspace q-th subband feature value λ_t(ω_q) characteristic of correspondence vector；

(12) all signal subspaces are extracted, the maximum I of each subband ambiguity function^t _max(r,z,ω_q) and corresponding Sound field locus (r^t _max(ω_q), z^t _max(ω_q))；

(13) by the sound field locus (r of all subbands^t _max(ω_q), z^t _max(ω_q)) arrange be a set omega (r^t _max (ω_q), z^t _max(ω_q)), and ambiguity in definition degree function again:

$\begin{array}{c}{I}_t(r,z,{\mathrm{\ω}}_q)=\\ \left\{\begin{array}{cc}0,& (r,z)\∉\mathrm{\Ω}\left(r_{\max }^t\right({\mathrm{\ω}}_q),z_{\max }^t({\mathrm{\ω}}_q\left)\right)\\ I_t(r,z,{\mathrm{\ω}}_q),& (r,z)\∈\mathrm{\Ω}\left(r_{\max }^t\right({\mathrm{\ω}}_q),z_{\max }^t({\mathrm{\ω}}_q\left)\right)\end{array}\right.\end{array}---\left(7\right)$

(14) first ambiguity function to all subbands is weighted adding up, then the result of all signal subspaces is tired out Add:

$I_{end}(r,z)={\mathrm{\Σ}}_{t=1}^T{\mathrm{\Σ}}_{q=1}^Q{\mathrm{\η}}_t\left({\mathrm{\ω}}_q\right)I_t(r,z,{\mathrm{\ω}}_q)---\left(8\right)$

And it is rendered as spatial image.

Fig. 4 is the image of two Extended target of last output, illustrates two targets at the sky residing for whole pond Between position, and its relative position.

Claims (1)

1. during weighting broadband, inverse operator decomposes multiple target acoustic imaging method, and the method comprises the steps:

(1) coordinate system is set up；If the transmitting-receiving of P transducer composition is closed puts linear array, it is vertically disposed in water, with linear array as z-axis, water Square setting up coordinate system to for r axle, the water surface is zero, and first transducer degree of depth away from the water surface is z₁, the P transducer The degree of depth away from the water surface is z_P；

(2) transducer array is received data snap and is arranged in column vector, and fast to each data according to predetermined sampling frequency Bat carries out the y that samples_i(n)=[y_1i(n)…y_pi(n)…y_Pi(n)]^H, p=1 ..., P, i=1 ..., M, n represent that moment, H represent and turn Putting, M represents the fast umber of beats of data；

(3) signal after sampling being carried out Short Time Fourier Transform, the n moment is represented by

$y_i\left(\mathrm{\ω}\right)=\left[\begin{array}{c}{\mathrm{\Σ}}_{l=0}^{L-1}{y}_{1i}(n-l)w\left(l\right)e^{-j\mathrm{\ω}(n-l)}\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{\Σ}}_{l=0}^{L-1}{y}_{Pi}(n-l)w\left(l\right)e^{-j\mathrm{\ω}(n-l)}\end{array}\right]---\left(1\right)$

Wherein ω represents that angular frequency, w (l) are window function sequence, and L is the length of window function；

(4) all of M data snap is arranged in matrix

Time then, inverse operator is

K (ω)=Y (ω) Y^H(ω)；

(5) by time inverse operator be divided into Q subband, carry out singular value decomposition:

K(ω_q)=U (ω_q)Λ(ω_q)V*(ω_q), q=1 ..., Q (3)

Wherein ω_qRepresenting the angular frequency that q-th subband is corresponding, Q is the number of subband.

(6) function curve of the eigenvalue after each subband noncoherent accumulation and horizontal range r is exported, according to the song of eigenvalue of maximum Line judges time inverse operator K corresponding to time window comprising target_o(ω_q), corresponding eigenvalue Λ_o(ω_q) and feature to Amount U_o(ω_q) and V_o*(ω_q)；

(7) threshold value is set, according to eigenvalue Λ_o(ω_q), q=1 ..., the variation tendency of Q judges signal subspace, it is assumed that signal The number of subspace is T；

(8) according to eigenvalue Λ_o(ω) design weight coefficient, finds out the eigenvalue λ of signal subspace_t(ω_q), q=1 ..., Q exists The maximum of frequency dimensionT=1 ... T, the eigenvalue of remaining each subband and maximumRatio be weight coefficient:

${\mathrm{\η}}_t\left({\mathrm{\ω}}_q\right)={\mathrm{\λ}}_t\left({\mathrm{\ω}}_q\right)/{\mathrm{\λ}}_{\max }^t,q=1.....Q,t=1,...T---\left(4\right)$

(9) be grid by Spacial domain decomposition interested, the intersection point of grid be suppose target place position (r, z), its Middle r represents that supposition target closes the horizontal range putting battle array away from transmitting-receiving, and z represents the depth of water supposing target；

(10) Underwater Acoustic Environment faced according to formation method, determines be suitable for propagation model, respectively obtain supposition target with Transmitting-receiving is closed and is put transmission vector g (r, z, the ω between battle array_q)=[g₁(z₁,r,z,ω_q)…g_p(z_p,r,z,ω_q)…g_P(z_P,r,z, ω_q)]^T, wherein g_p(z_p, r, z, ω_q) represent the transmission function supposed between target and pth transducer.

(11) each signal subspace, the ambiguity function of each subband is:

I_t(r, z, ω_q)=| g^H(r, z, ω_q)u_t(ω_q)|², q=1 ..., Q, t=1 ..., T (5)

Wherein u_t(ω_q) represent the t signal subspace q-th subband feature value λ_t(ω_q) characteristic of correspondence vector, H represents altogether Yoke transposition；

(12) all signal subspaces are extracted, the maximum I of each subband ambiguity function^t _max(r,z,ω_q) and the sound of correspondence Field locus (r^t _max(ω_q),z^t _max(ω_q))；

(13) by the sound field locus (r of all subbands^t _max(ω_q),z^t _max(ω_q)) arrange be a set omega (r^t _max (ω_q),z^t _max(ω_q)), and ambiguity in definition degree function again:

$\begin{array}{c}{I}_t\left(r,z,{\mathrm{\ω}}_q\right)=\\ \left\{\begin{array}{cc}0,& \left(r,z\right)\∉\mathrm{\Ω}\left(r_{\max }^t\left({\mathrm{\ω}}_q\right),z_{\max }^t\left({\mathrm{\ω}}_q\right)\right)\\ I_t\left(r,z,{\mathrm{\ω}}_q\right),& \left(r,z\right)\∈\mathrm{\Ω}\left(r_{\max }^t\left({\mathrm{\ω}}_q\right),z_{\max }^t\left({\mathrm{\ω}}_q\right)\right)\end{array}\right.\end{array}---\left(6\right)$

(14) first it is weighted adding up by the ambiguity function of all subbands, then the result of all signal subspaces is tired out Add:

Final result is rendered as space diagram Picture.