CN106019290B - Inverse operator decomposes multiple target acoustic imaging method when weighting broadband - Google Patents
Inverse operator decomposes multiple target acoustic imaging method when weighting broadband Download PDFInfo
- Publication number
- CN106019290B CN106019290B CN201610356675.9A CN201610356675A CN106019290B CN 106019290 B CN106019290 B CN 106019290B CN 201610356675 A CN201610356675 A CN 201610356675A CN 106019290 B CN106019290 B CN 106019290B
- Authority
- CN
- China
- Prior art keywords
- subband
- target
- signal
- max
- function
- 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
Links
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
- G01S15/89—Sonar systems specially adapted for specific applications for mapping or imaging
Abstract
Inverse operator decomposes multiple target acoustic imaging method when a kind of weighting broadband, and core is：Transmitting-receiving is closed and sets the signal progress Fourier decomposition that battle array receives, it is divided into each subband by narrow band filter, inverse operator singular value decomposition when each subband is corresponding, the number of signal subspace is judged according to the variation tendency of characteristic value, ambiguity function will be obtained after the corresponding feature vector of all signal subspaces and the transmission vector correlation based on propagation model, extract the maximum value of each subband ambiguity function and corresponding two-dimensional acoustic field space, weighting coefficient is designed with characteristic value simultaneously, the maximum value of Weighted Fuzzy degree function, by the weighted results coherent accumulation of all subbands, again by the result coherent accumulation of all signal subspaces, it is rendered as spatial image.
Description
Technical field
The invention belongs to marine resources development and the acoustic detection field utilized, specifically a kind of to use array to multiple
The method that point target or autgmentability target carry out acoustic imaging.
Background technology
Since the performance that sound wave is propagated in water is better than light wave and electromagnetic wave, so, it is applied to the detection base of marine resources
It is based on acoustic means on this.Current practice is still with last century five, six in the spatial processing technique of sonar system
The ten's rise and grow up with the plane wave Wave beam forming under white noise background and free field propagation conditions.However,
Under neritic environment, traditional processing method based on plane wave model is no longer applicable in, since propagation of the acoustical signal in ocean is in
Existing multipath/more Normal mode analysis state, it is easy to cause the pseudo- target in detection process or diplopia phenomenon, influence the inspection of sonar system
Survey performance.When inverse processing since last century end be applied to underwater sound field, be at present quite have application prospect active target detect
Means, when inverse operator to decompose be to reach a kind of side of focus emission, positioning and imageable target by analyzing target scattering matrix
Method.
However, inverse operator decomposition method is applied to the multi-target imaging in marine exploration field, the side of generally use when current
Method, one assumes that different targets corresponds to different feature vectors, and the corresponding feature vector of target is related to transfer matrix, two
It is the MUSIC algorithms using noise subspace, there are following limitations：(1) processing of broadband signal only considers centre frequency
Or the mode of the handling result noncoherent accumulation of entire frequency range, do not effectively utilize the information of broadband signal；(2) multiple mesh
Mark needs multiple image to show, cannot intuitively indicate the relative position relation of multiple targets；(3) scattering of target itself is special
The comprehensive functions such as property, waveguide propagation effect and the relationship of energy converter and frequency, are eventually exhibited as target scattering characteristics and frequency
Relation function, and do not consider and utilize this scattering properties in processing procedure at present.
Invention content
The present invention will overcome the disadvantages mentioned above of the prior art, and inverse operator decomposes multiple target acoustic imaging side when providing weighting broadband
Method, this method on the basis of inverse operator decomposition method, effectively utilize target scattering characteristics and wide-band-message when existing, from
And achievees the purpose that eliminate pseudo- target, effectively display autgmentability spatial characters and show multiple targets simultaneously.
Inverse operator decomposes multiple target acoustic imaging method when weighting broadband proposed by the present invention, according to the variation tendency of characteristic value
It determines signal subspace, extracts each signal subspace, the maximum value of each subband ambiguity function and corresponding two-dimentional sound
Weighting coefficient is designed in space according to the functional relation of maximum eigenvalue and frequency, weights corresponding signal subspace and subband
Ambiguity function, and by the ambiguity function noncoherent accumulation after weighting, final result is shown as sound field images, it determines simultaneously
The spatial position of multiple targets.
The method of the present invention is described further below.
Have to the parameters,acoustic such as depth of water, Sound speed profile, density, the sediment parameters of environment in the case of understanding enough,
The identification of the three-dimensional acoustics image and Target space position of target is completed as follows：
(1) coordinate system is established；If linear array is set in the transmitting-receiving conjunction of P energy converter composition, it is disposed vertically in water, is z with linear array
Axis, horizontal direction are that r axis establishes coordinate system, and the water surface is coordinate origin, and first depth of the energy converter away from the water surface is z1, P
Depth of the energy converter away from the water surface is zP；
(2) transducer array is received into data snap and is arranged in column vector, according to predetermined sampling frequency to each data
Snap is sampled, yi(n)=[y1i(n)…ypi(n)…yPi(n)]H, p=1 ..., P, i=1 ..., M, wherein n the expression moment,
H indicates that transposition, M indicate data number of snapshots；
(3) Short Time Fourier Transform is carried out to the signal after sampling, the n moment is represented by
Wherein w (l) is window function sequence, and L is the length of window function；
(4) M all data snaps is arranged in matrix
Inverse operator is K (ω)=Y (ω) Y when thenH(ω)；
(5) by when 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 for exporting characteristic value and horizontal distance r after each subband noncoherent accumulation, according to maximum eigenvalue
Curve judge the when inverse operator K corresponding to the time window comprising targeto(ωq), corresponding characteristic value Λo(ωq) and it is special
Levy vector Uo(ωq) and Vo*(ωq)；
(7) threshold value is set, according to characteristic value Λo(ωq), the variation tendency of q=1 ..., Q judge signal subspace, it is assumed that
The number of signal subspace is T；
(8) according to characteristic value Λo(ω) designs weighting coefficient, finds out the eigenvalue λ of signal subspacet(ωq), q=
1 ..., the maximum value λ that Q is tieed up in frequencyt max, t=1 ... T, characteristic value and the maximum value λ of remaining each subbandt maxRatio be plus
Weight coefficient：
ηt(ωq)=λt(ωq)/λt max, q=1 ..., Q, t=1 ... T (4)
(9) be grid by interested Spacial domain decomposition, the intersection point of grid be where assuming target position (r,
Z), wherein r indicates to assume that target closes the horizontal distance for setting battle array away from transmitting-receiving, and z indicates to assume the depth of water of target；
(10) Underwater Acoustic Environment faced according to imaging method determines the propagation model being applicable in, respectively obtains hypothesis mesh
Transmission vector g (r, z, the ω set between battle array is closed in mark and transmitting-receivingq)=[g1(z1,r,z,ωq)…gj(zj,r,z,ωq)…gP(zP,
r,z,ωq)]T, wherein gp(zp,r,z,ωq) indicate to assume the transmission function between target and p-th of energy converter.
(11) ambiguity function of each signal subspace, each subband is：
It(r,z,ωq)=| gH(r,z,ωq)ut(ωq)|2, q=1 ..., Q, t=1 ..., T (5)
Wherein ut(ωq) indicate t-th of signal subspace, q-th of subband feature value λt(ωq) corresponding feature vector, H tables
Show conjugate transposition；
(12) all signal subspaces, the maximum value I of each subband ambiguity function are extractedt max(r,z,ωq) and its it is corresponding
Sound field spatial position (rt max(ωq), zt max(ωq))；
(13) by the sound field spatial position (r of all subbandst max(ωq), zt max(ωq)) arrange as a set omega (rt max
(ωq), zt max(ωq)), and ambiguity in definition degree function again：
(14) to all subbands, the ambiguity functions of all signal subspaces is weighted cumulative：
The ambiguity function finally obtained is shown as to the 3-D view of distance r and depth z, can determine multiple targets simultaneously
Corresponding two-dimensional space region.
It is an advantage of the invention that：Without judging the correspondence between each target and characteristic value and feature vector, fully
Using signal subspace information, multiple targets are presented in same piece image, clearly demonstrate the opposite position between target
It sets.
Description of the drawings
Experimental arrangement and coordinate the setting figure of Fig. 1 present invention.
The function relation figure of the characteristic value and distance of Fig. 2 present invention.
The characteristic value of time window and the function relation figure of frequency where Fig. 3 targets of the present invention.
Fig. 4 is the image of two extension targets of the present invention.
Specific embodiment
Below in conjunction with the accompanying drawings, the present invention is further described by specific embodiment.It is 1.44m in a depth
Waveguide laboratory water tank in, the three face paste sound eliminating tiles in pond, the bottom in pond has spread the sand of one layer of 0.22m thickness.The ring in pond
Border parameter is as follows：The velocity of sound is constant in water, is learnt as c by measuring water temperature calculating1=1493m/s, water body density p1=
1000kg/m3, sediment parameters density p2=1800kg/m3, velocity of sound c2=1650m/s, attenuation coefficient α2=0.67dB/ λ,
The density of substrate is ρ3=1800kg/m3, velocity of sound c3=1580m/s, attenuation coefficient α3=0.8dB/ λ.
(1) in the present embodiment, the foundation of experimental arrangement and coordinate system sets 32 array elements of battle array as shown in Figure 1, transmitting-receiving is closed, entirely
Field is structured the formation, and first array element is vertically disposed in away from water surface 0.04m, array element spacing 0.04m in pond；Two a diameter of 0.21m,
Length is the cylindrical type target of 0.51m, and a basin bottom for being placed on matrix 6m, one is placed on matrix 8m, and depth is
At 0.83m.Transmitting signal is that the centre frequency that pulse width is 0.5ms is 18kHz, and bandwidth is the linear FM signal of 5kHz.
(2) transducer array is received into data snap and is arranged in column vector, according to sample frequency be 50kHz to every number
It is sampled according to snap, yi(n)=[y1i(n)…ypi(n)…yPi(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., points for 200 rectangular window, to the letter after sampling
Number carry out Short Time Fourier Transform
(4) 32 all data snaps are arranged in matrix
Inverse operator is K (ω)=Y (ω) Y when thenH(ω)；
(5) the when inverse operator of the bandwidth of 5kHz is divided into Q=501 son band and carries out singular value decomposition：
K (ω)=U (ω) Λ (ω) V* (ω) (3)
(6) function curve for exporting characteristic value and horizontal distance r after each subband noncoherent accumulation, judges comprising target
When inverse operator K corresponding to time windowo(ω), corresponding characteristic value Λo(ω) and feature vector Uo(ω) and Vo*(ω)；
Fig. 2 is the function relation figure of preceding 5 characteristic values and distance.As shown in Figure 2, there are two local peakings for characteristic value, are located at two
At distance where a target, the time window for containing two targets may thereby determine that.
(7) left side figure is Λ in Fig. 3o(ωq) with the variation tendency of frequency, it can be determined that signal subspace T=4；
(8) weighting coefficient is designed, Fig. 3 top right plots are the weighting coefficients of subspace corresponding to maximum eigenvalue, and bottom-right graph is
The weighting coefficient of subspace corresponding to Second Eigenvalue；
ηt(ωq)=λt(ωq)/λt max, q=1 ..., Q, t=1 ... T (4)
(9) it is grid by interested Spacial domain decomposition, horizontal starting distance is 0.1m, step-size in search 0.05m, is cut
Only distance is 15m, and vertical initial depth is 0m, and step-size in search 0.0075m, cut-off depth is 1.5m.The intersection point of grid is
It is assumed that the position (r, z) where target, wherein r indicates to assume that target closes the horizontal distance for setting battle array away from transmitting-receiving, and z indicates to assume target
The depth of water；
(10) Underwater Acoustic Environment faced according to imaging method determines the propagation model being applicable in, respectively obtains hypothesis mesh
Transmission vector g (r, z, the ω set between battle array is closed in mark and transmitting-receivingq)=[g1(z1,r,z,ωq)…gp(zp,r,z,ωq)…gP(zP,
r,z,ωq)]T, wherein gp(zp,r,z,ωq) indicate to assume the transmission function between target and p-th of energy converter, this embodiment party
Case uses Normal mode analysis propagation model, transmits vector and is：
Wherein ZlIndicate the corresponding characteristic function of l propagating modes, H0 (2)() is Hankel function, κlIndicate wave number；
(11) ambiguity function of each signal subspace, each subband is：
It(r,z,ωq)=| gH(r,z,ωq)ut(ωq)|2, q=1 ..., Q, t=1 ..., T (6)
Wherein ut(ωq) indicate t-th of signal subspace, q-th of subband feature value λt(ωq) corresponding feature vector；
(12) all signal subspaces, the maximum value I of each subband ambiguity function are extractedt max(r,z,ωq) and its it is corresponding
Sound field spatial position (rt max(ωq), zt max(ωq))；
(13) by the sound field spatial position (r of all subbandst max(ωq), zt max(ωq)) arrange as a set omega (rt max
(ωq), zt max(ωq)), and ambiguity in definition degree function again：
(14) first the ambiguity function of all subbands is weighted cumulative, then the result of all signal subspaces is tired out
Add：
And it is rendered as spatial image.
Fig. 4 is the image of the two extension targets finally exported, illustrates sky of two targets residing for entire pond
Between position and its relative position.
Claims (1)
1. inverse operator decomposes multiple target acoustic imaging method when weighting broadband, this method comprises the following steps：
(1) coordinate system is established；If linear array is set in the transmitting-receiving conjunction of P energy converter composition, it is disposed vertically in water, using linear array as z-axis, water
Square establish coordinate system to for r axis, the water surface is coordinate origin, and first depth of the energy converter away from the water surface is z1, the P energy converter
Depth away from the water surface is zP；
(2) transducer array is received into data snap and is arranged in column vector, and is fast to each data according to predetermined sampling frequency
Bat carries out sampling yi(n)=[y1i(n)…ypi(n)…yPi(n)]H, p=1 ..., P, i=1 ..., M, n indicate that moment, H indicate to turn
It sets, M indicates data number of snapshots；
(3) Short Time Fourier Transform is carried out to the signal after sampling, the n moment is represented by
Wherein ω indicates that angular frequency, w (l) are window function sequence, and L is the length of window function；
(4) M all data snaps is arranged in matrix
Inverse operator is when then
K (ω)=Y (ω) YH(ω)；
(5) by when inverse operator be divided into Q subband, carry out singular value decomposition：
K(ωq)=U (ωq)Λ(ωq)V*(ωq), q=1 ..., Q (3)
Wherein ωqIndicate that the corresponding angular frequency of q-th of subband, Q are the number of subband;
(6) function curve for exporting characteristic value and horizontal distance r after each subband noncoherent accumulation, according to the song of maximum eigenvalue
Line judges the when inverse operator K corresponding to the time window comprising targeto(ωq), corresponding characteristic value Λo(ωq) and feature to
Measure Uo(ωq) and Vo*(ωq)；
(7) threshold value is set, according to characteristic value Λo(ωq), the variation tendency of q=1 ..., Q judge signal subspace, it is assumed that signal
The number of subspace is T；
(8) according to characteristic value Λo(ω) designs weighting coefficient, finds out the eigenvalue λ of signal subspacet(ωq), q=1 ..., Q exists
The maximum value of frequency dimensionT=1 ... T, the characteristic value and maximum value of remaining each subbandRatio be weighting coefficient：
(9) it is grid by interested Spacial domain decomposition, the intersection point of grid is the position (r, z) assumed where target,
Middle r indicates to assume that the horizontal distance that target sets battle array away from transmitting-receiving conjunction, z indicate to assume the depth of water of target；
(10) Underwater Acoustic Environment faced according to imaging method determines the propagation model that is applicable in, respectively obtain assume target and
Transmission vector g (r, z, the ω set between battle array is closed in transmitting-receivingq)=[g1(z1,r,z,ωq)…gp(zp,r,z,ωq)…gP(zP,r,z,
ωq)]T, wherein gp(zp, r, z, ωq) indicate to assume the transmission function between target and p-th of energy converter;
(11) ambiguity function of each signal subspace, each subband is：
It(r, z, ωq)=| gH(r, z, ωq)ut(ωq)|2, q=1 ..., Q, t=1 ..., T (5)
Wherein ut(ωq) indicate t-th of signal subspace, q-th of subband feature value λt(ωq) corresponding feature vector, H indicates altogether
Yoke transposition；
(12) all signal subspaces, the maximum value I of each subband ambiguity function are extractedt max(r,z,ωq) and its corresponding sound
Field spatial position (rt max(ωq),zt max(ωq))；
(13) by the sound field spatial position (r of all subbandst max(ωq),zt max(ωq)) arrange as a set omega (rt max
(ωq),zt max(ωq)), and ambiguity in definition degree function again：
(14) first the ambiguity function of all subbands is weighted cumulative, then the result of all signal subspaces is tired out
Add：
Final result is rendered as spatial image.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610356675.9A CN106019290B (en) | 2016-05-26 | 2016-05-26 | Inverse operator decomposes multiple target acoustic imaging method when weighting broadband |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610356675.9A CN106019290B (en) | 2016-05-26 | 2016-05-26 | Inverse operator decomposes multiple target acoustic imaging method when weighting broadband |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN106019290A CN106019290A (en) | 2016-10-12 |
| CN106019290B true CN106019290B (en) | 2018-11-13 |
Family
ID=57094384
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610356675.9A Active CN106019290B (en) | 2016-05-26 | 2016-05-26 | Inverse operator decomposes multiple target acoustic imaging method when weighting broadband |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106019290B (en) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110161546B (en) * | 2019-05-23 | 2021-04-16 | 杭州中科微电子有限公司 | Satellite orientation device and method using iterative weighted ambiguity function method |
| CN110456342B (en) * | 2019-07-11 | 2023-03-21 | 西安电子科技大学 | Far-field multi-moving-object detection method of single-transmitting-antenna radar |
| CN110455931B (en) * | 2019-09-05 | 2020-07-28 | 中国科学院声学研究所 | Target detection and positioning method in layered medium |
| CN111856474B (en) * | 2020-07-30 | 2023-07-25 | 重庆大学 | Subarray-based space-time domain conditional coherence coefficient ultrasonic imaging method |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104280737A (en) * | 2014-08-29 | 2015-01-14 | 浙江工业大学 | Weighted broadband time reversal operator resolution acoustic imaging method |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7649805B2 (en) * | 2007-09-12 | 2010-01-19 | Schlumberger Technology Corporation | Dispersion extraction for acoustic data using time frequency analysis |
| CN101387701B (en) * | 2008-10-24 | 2011-01-05 | 西北工业大学 | Passive time reversal reverberation suppression method based on forward prediction |
| CN103076590A (en) * | 2012-12-31 | 2013-05-01 | 东南大学 | Method for positioning underwater sound pulse signal on basis of frequency estimation |
| CN103197282B (en) * | 2013-03-18 | 2015-12-02 | 哈尔滨工程大学 | Anti-during MVDR based on Amplitude Compensation focus on localization method |
-
2016
- 2016-05-26 CN CN201610356675.9A patent/CN106019290B/en active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104280737A (en) * | 2014-08-29 | 2015-01-14 | 浙江工业大学 | Weighted broadband time reversal operator resolution acoustic imaging method |
Non-Patent Citations (2)
| Title |
|---|
| "Experimental investigation of selective localization by decomposition of the time reversal operator and subspace-based technique";Li. et al;《IET RADAR Sonar and Navigation》;20081231;第2卷(第6期);第428页 * |
| "干扰环境下窄带信号时反算子分解聚焦性能分析";李春晓 等;《哈尔滨工程大学学报》;20080831;第29卷(第8期);第868页左栏倒数第2行-右栏第1行 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN106019290A (en) | 2016-10-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US11733381B2 (en) | Sound velocity profile inversion method based on inverted multi-beam echo sounder | |
| CN108226933B (en) | Deep sea broadband target depth estimation method based on fringe interference structure | |
| CN112083404B (en) | Single-vector hydrophone sound source depth estimation method based on multi-path feature matching | |
| CN106019290B (en) | Inverse operator decomposes multiple target acoustic imaging method when weighting broadband | |
| Gervaise et al. | Passive geoacoustic inversion with a single hydrophone using broadband ship noise | |
| Duan et al. | A reliable acoustic path: Physical properties and a source localization method | |
| CN107272005B (en) | Active positioning method based on target echo arrival time delay and arrival angle under reliable acoustic path | |
| CN104820218B (en) | Shallow sea bottom single-parameter inversion method based on frequency domain autocorrelation | |
| RU2603724C2 (en) | Method and device to control acoustic characteristics of network of acoustic nodes located along towed acoustic linear antennae | |
| Gingras et al. | Electromagnetic matched-field processing: Basic concepts and tropospheric simulations | |
| Zhao et al. | Open-lake experimental investigation of azimuth angle estimation using a single acoustic vector sensor | |
| Jiang et al. | Estimation of geoacoustic properties of marine sediment using a hybrid differential evolution inversion method | |
| Cho et al. | Impact of array tilt on source-range estimation in shallow water using the array invariant | |
| Pan et al. | IoUT based underwater target localization in the presence of time synchronization attacks | |
| CN104280737B (en) | Weighted broadband time reversal operator resolution acoustic imaging method | |
| Lucifredi et al. | Gray whale target strength measurements and the analysis of the backscattered response | |
| Qin et al. | Sound propagation from the shelfbreak to deep water? | |
| CN113126029B (en) | Multi-sensor pulse sound source positioning method suitable for deep sea reliable acoustic path environment | |
| Zhang et al. | Simulation of ship radiated noise field in deep sea based on statistical characteristics of sound source | |
| Gerstoft et al. | Estimation of transmission loss in the presence of geoacoustic inversion uncertainty | |
| Yang et al. | Bayesian passive acoustic tracking of a cooperative moving source in shallow water | |
| Fang et al. | The echolocation transmission beam of free-ranging Indo-Pacific humpback dolphins (Sousa chinensis) | |
| Aksenov et al. | Invariance of the effective phase velocity of a hydroacoustic field in the deep ocean | |
| CN117406172A (en) | Sea depth estimation method based on tug interference characteristics in deep sea environment | |
| Talebpour | Single Hydrophone Underwater Localization Approach in Sallow Waters |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |