CN102176008A - Phased azimuth filtering method for three-dimensional stratum imaging - Google Patents

Abstract

The invention discloses a phased azimuth filtering method for three-dimensional stratum imaging. For a plurality of coherent single snapshot signals from an identical distance and different angles, optimal source signal waveform estimates on each scanning azimuth of a space are obtained by adopting a maximum posterior probability criterion, and the azimuth of a target is judged according to the energy of the signals on each azimuth. The method has the advantages that: 1, a covariance matrix is not required for direct or indirect decomposition processing, the limitations of a snapshot number are avoided, and active single snapshot data can be processed; 2, a coherent source can be estimated and distinguished by adopting the maximum posterior probability criterion, decorrelation processing is not required to be performed on a highly-coherent signal source in a water body or a stratum, and coherent interference resistance is high; and 3, the signals returned from the identical distance can be divided into the summation of source signals of sufficient azimuth grids, iterative convergence can be performed by utilizing a cost function and a constraint equation to ensure the optical source signal waveform estimates on each scanning azimuth, and relatively higher processing reliability and stable performance are ensured.

Description

A kind of phased method of bearing filtering of three-dimensional formation imaging

Technical field

The invention belongs to the orientation high resolution battle array process field of small-sized latent device, mainly is a kind of phased method of bearing filtering of three-dimensional formation imaging.

Background technology

Still the 3-D imaging system that does not have at present ripe stratum, seabed in the world, this programme adopt stratum, the seabed acoustics three-D imaging method of imitative Medical CT technology to realize burying thunder and survey identification.On three-dimensional imaging mechanism method, handle the realization level orientation high resolution different by the high-resolution battle array to echo, obtain depth direction apart from high resolution by pulse compression technique, launch the two-dimentional energy profile that obtains stratum, a transversal section each time by transmitting transducer, along with sonar platform travels forward, obtain another and tie up away synthetic aperture high resolution on the boat orientation, continuous sweep can obtain the three-dimensional detection imaging of stratum landforms.Mainly need to break through the gordian technique of each dimension imaging resolving power.It is little that this small-sized latent device belongs to physical size, the integrated level height, the active sonar system that translational speed is fast, when the level orientation dimension is handled, be subjected to physical size limitations, the bearing resolution that conventional wave beam forms is nowhere near, and traditional MVDR, high-resolution algorithms such as MUSIC are directly to handle or the aftertreatment of feature decomposition signal subspace from signal covariance matrix, handling low signal-to-noise ratio, the fast enough irrelevant targets of umber of beats all have outstanding effect, but during the strong coherent source of stratigraphic distribution formula under handling single fast beat of data, the on the one hand strong coherent signal distortion covariance matrix of crosstalking, single on the other hand fast beat of data can't repeatedly be added up signal covariance matrix is carried out smoothing processing, causes these traditional high resolution algorithms based on the estimate covariance matrix can't reach the requirement of azimuth dimension resolving accuracy.

Summary of the invention

Detection, drafting, three-dimensional imaging to the stratum, seabed shows in small-sized towed body in order to solve, the high resolution that particularly solves the level orientation dimension in coherent signal source in the stratum under the little situation of physical size is differentiated, and the present invention proposes a kind of phased method of bearing filtering of three-dimensional formation imaging.The objective of the invention is at same distance, different angles, single snap, heterogeneous dried source signal, avoid the thinking of traditional estimate covariance matrix, the source signal optimum waveform that adopts maximum posteriori criterion to obtain on each scan position of space is estimated, judges target direction according to the energy size of signal on each orientation; Can be applicable to landforms, landform, the stratum detecting of the middle-size and small-size latent device in ocean, for example in the AUV/UUV system.

The present invention solves the technical scheme that its technical matters adopts:

To achieve these goals, the invention provides a kind of stratum, seabed acoustics 3-D imaging system of imitative Medical CT technology, adopt synthetic aperture technique in the azimuth dimension of walking to navigate (y axle), adopt the wideband pulse compress technique in degree of depth azimuth dimension (z axle), phased azimuth filtering (the being called for short PAF) method that adopts the present invention to propose in level orientation dimension (x axle) obtains the high resolution technique of orientation angles.System platform is in the restriction that is subjected to travelling speed along the azimuth dimension of walking to navigate, need to satisfy 1/2 spatial sampling rule, system along walk to navigate azimuth dimension forward straight line move, produce the wideband frequency modulation signal by combination transducer lift-off technology, arrange the hydrophone array that polynary rectangular shape is arranged.

On the depth dimension the every batch of echoed signal that receives is carried out the broadband signal pulse compression technique, according to the relation of range resolution and the signal bandwidth and the velocity of sound, ζ=C/2B can obtain the resolution parameter value of seabed multiple goal on depth dimension;

Form the short linear array of one group of physical pore size by a plurality of receiving hydrophones of installing perpendicular to the platform motion direction on the level orientation dimension, in receiving every lot data, a large amount of water bodys that are carved into the different azimuth that reaches simultaneously, the multiple goal on stratum, strong relevant, single snap echo data, the field angle that adopts conventional beamforming algorithm to obtain is θ=λ/L, λ is a signal wavelength, the wavelength of generally getting centre frequency calculates, L is a linear array physical pore size length, small-sized latent applicator platform has limited linear array physical pore size length, to such an extent as to field angle can't be differentiated stratum " target complex " too greatly, handle or the disposal route of dividing solution subspace performance when handling relevant single fast beat of data descends serious based on the covariance matrix of classics, the present invention adopts phased method of bearing filtering that data are adopted pre-division module respectively, the array conversion module, MAP criterion estimation module, the iteration convergence module obtains " target complex " according to depth layer, the high resolution parameter value of different azimuth angle;

On the azimuth dimension of walking to navigate, after the many receptions battle array that is parallel to the platform motion direction receives certain pore size length echo data, the multiple submatrixes synthetic aperture technique of carrying out after the hardware and software motion compensation is handled, it doesn't matter for resolving power and target range, the array element physical size is relevant with receiving, ξ=D/2, wherein D is single receiving hydrophone physics array element length, the target high resolution parameter on the azimuth dimension that obtains to walk to navigate.

Obtain the transversal section two dimension energy profile of " target complex " degree of depth and the level orientation on stratum after every batch, launch propelling one by one according to the azimuth dimension of walking to navigate, frame number superposes successively, adds the resolution parameter value of walking to navigate, and obtains the 3-D display result on stratum.

Phased method of bearing filtering comprises following step:

One data preprocessing module according to actual needs, the initialization system parameter, comprise signal type and pulsewidth, target depth range, layering is apart from band, orientation angles scope, the distinguishable angle in minimum orientation, choose the time samples point that arrives echo data, the stratum echoed signal of every layer of arrival evenly is divided into L _θIndividual orientation grid, the set of orientation grid are not more than transmitting transducer emission directive property width.

Two array matrix conversion modules comprise equation equation step, parameter relationship formula step two parts, by transforming target empty bay matrix, feature space vector, eigenwert, the array manifold matrix after acquisition receives echo space of matrices, the transformed matrix space that receives echo, target empty bay matrix, weighting, set up equation expression formula and parameter relationship formula between the above-mentioned parameter.

Three MAP criterion estimation module are estimated echo intensity on the orientation with maximum posteriori criterion (MAP), according to maximum posteriori criterion, set up cost function and equation of constraint.

Four iteration convergence modules can be converted into no constrained minimization and find the solution problem according to cost function and equation of constraint, and iterating makes the cost function minimum, obtain the layering distance with on the optimal estimation value of orientation grid echoed signal intensity.

On the described array matrix conversion module, comprise equation equation step, parameter relationship formula step two parts.Under the far field condition, last each array element of a certain layering distance receives echo vector and array manifold matrix and echo signal source vector, noise vector equation and is in can be according to ARRAY PROCESSING on equation equation step:

The stratum echoed signal evenly is divided into L _θIndividual orientation grid (L _θ＞＞M (array element number)), s ₁(t), s ₂(t) ...,

Source signal for the different azimuth grid satisfies far field condition,

M ∈ (1,2 ..., M), l ∈ (1 ..., L _θ) be a series of source signal (1 ..., L _θSource signal) each array element (1 ... M array element) frequency-domain expression of the phase differential that produces, n _1...M(t) noise sequence is the zero-mean Gaussian process, and noise is separate between array element, and each sensor of array is an isotropy, the inconsistent interference of no mutual coupling and passage.The equation equation further is expressed as follows:

$Z\left(t\right)=\underset{l=1}{\overset{L_{\mathrm{\θ}}}{\mathrm{\Σ}}}v\left({\mathrm{\φ}}_l\right)s_l\left(t\right)=V({\mathrm{\φ}}_1,...,{\mathrm{\φ}}_{L_{\mathrm{\θ}}})\·S\left(t\right)+N\left(t\right)\⇒Z=V(S+\mathrm{\ΔS})=V\hat{S}$

Z (t)=[z wherein ₁(t) ..., z _M(t)] ^TBe the array element received signal vector of M * 1,

Be L _θ* 1 source signal vector,

Be the array manifold vector of M * 1 of l source signal on M array element,

Be M * L _θThe array manifold matrix, N (t)=[n ₁(t) ..., n _M(t)] ^TBe the noise vector of M * 1,

For sneaking into the target empty bay matrix after the incoherent noise, φ _l=2 π f _cDcos θ _l/ c, wherein f _cBe centre frequency, d is an array element distance, and c is the velocity of sound in the water.

In the described parameter relationship formula step, the reception echo vector that at first will obtain, array manifold matrix, target empty bay matrix three's equation equation is further derived:

By $Z=V\hat{S}$

$\⇒V^{H}Z=V^{H}V\hat{S}=\mathrm{P\Λ}{P}^H\hat{S},$ Wherein $\mathrm{\Λ}=\left[\begin{array}{cc}{\mathrm{\Λ}}_{\mathrm{fix}}& \\ & {\mathrm{\Λ}}_{\mathrm{free}}\end{array}\right]$

Λ _Fix, Λ _FreeRepresent to obtain behind the conjugate multiplication array manifold matrix characteristics of decomposition value of matrix respectively, obtain nonzero eigenvalue Λ according to the diagonal matrix descending sort _FixWith zero eigenvalue Λ _Free

Further the matrix equivalence is derived as follows

$\⇒P^H{V}^{H}Z={\mathrm{\ΛP}}^H\hat{S}\⇒\mathrm{QZ}={\mathrm{\ΛP}}^H\hat{S}(Q=P^H{V}^H)\⇒Z^{\′}={\mathrm{\ΛS}}^{\′}(Z^{\′}=\mathrm{QZ},S^{\′}=P^H\hat{S})$

Adopt the matrix equivalence to replace, receive the echo vector apart from each array element in the layering and introduce array manifold Matrix Conjugate value and feature space vector conjugation information, the equation equation left side obtains to receive the transformed matrix space Z ' of echo, with the product of equation the right feature space vector conjugation information and target empty bay matrix, the target empty bay matrix S after the weighting after equation equation the right obtains to transform '.

By above-mentioned conversion, utilize the mutually orthogonal characteristic in nonzero eigenvalue subspace and zero eigenvalue subspace, Z ' _FixAnd S ' _FixBe nonzero eigenvalue subspace, Z ' _FreeAnd S ' _FreeBe the zero eigenvalue subspace, can will receive the transformed matrix space Z ' of echo and the target empty bay matrix S after the weighting ' all carry out the subspace to expand into:

$\left[\begin{array}{c}{Z}_{\mathrm{fix}}^{\′}\\ Z_{\mathrm{free}}^{\′}\end{array}\right]=\left[\begin{array}{cc}{\mathrm{\Λ}}_{\mathrm{fix}}& \\ & {\mathrm{\Λ}}_{\mathrm{free}}\end{array}\right]\left[\begin{array}{c}{S}_{\mathrm{fix}}^{\′}\\ S_{\mathrm{free}}^{\′}\end{array}\right]=\left[\begin{array}{c}{\mathrm{\Λ}}_{\mathrm{fix}}{S}_{\mathrm{fix}}^{\′}\\ {\mathrm{\Λ}}_{\mathrm{free}}{S}_{\mathrm{free}}^{\′}\end{array}\right]$

The Z ' of main energy distribution in the subspace of acquisition Z ' and S ' _FixAnd S ' _FixRelational expression, each element can be written as:

$Z_{\mathrm{fix}}^{\′}={\mathrm{\Λ}}_{\mathrm{fix}}{S}_{\mathrm{fix}}^{\′}\⇒S_{\mathrm{fix}}^{\′}={{\mathrm{\Λ}}_{\mathrm{fix}}}^{-1}{Z}_{\mathrm{fix}}^{\′}$

$z_{\mathrm{fix},l}^{\′}={\mathrm{\λ}}_l{s}_{\mathrm{fix},l}^{\′}\⇒s_{\mathrm{fix},l}^{\′}=\frac{1}{{\mathrm{\λ}}_l}{z}_{\mathrm{fix},l}^{\′},$ l＝1，...，L _fix

Energy distribution S ' in the target empty bay matrix after the weighting _FixCan obtain by the transformed matrix space Z ' that receives echo, can obtain target empty bay matrix successively.

$S^{\′}=P^H\hat{S}$

$\⇒\hat{S}={\mathrm{PS}}^{\′}=\left[\begin{array}{cc}{P}_{\mathrm{fix}}& P_{\mathrm{free}}\end{array}\right]\left[\begin{array}{c}{S}_{\mathrm{fix}}^{\′}\\ S_{\mathrm{free}}^{\′}\end{array}\right]P_{\mathrm{fix}}{S}_{\mathrm{fix}}^{\′}+P_{\mathrm{free}}{S}_{\mathrm{free}}^{\′}$

Take turns acquisition series of parameters relational expression in the step at this, include target empty bay matrix

Target empty bay matrix S after the weighting ', feature space vector P, eigenwert Λ, array manifold matrix V, receive echo transformed matrix space Z ', receive echo space of matrices Z.Target empty bay matrix wherein

It is the parameter that needs estimation to obtain.

$\⇒\left\{\begin{array}{c}\hat{S}={\mathrm{PS}}^{\′}\\ S_{\mathrm{fix}}^{\′}={{\mathrm{\Λ}}_{\mathrm{fix}}}^{-1}{Z}_{\mathrm{fix}}^{\′}\\ Z^{\′}=P^H{V}^{H}Z\\ Z=V\hat{S}\end{array}\right.$

Described MAP criterion estimation module comprises that (MAP) estimates θ with maximum posteriori criterion to the reception data vector Z (t) of each moment t _lEcho intensity on the orientation

It is target empty bay matrix

The orientation output spectra has this orientation of expression of peak value that the strong scattering target is arranged, and the height of peak value is represented the estimated value to this orientation target scattering intensity.

Choose cost function and equation of constraint is as follows:

$\left\{\begin{array}{c}\mathrm{\Γ}=\underset{l=1}{\overset{L_{\mathrm{\θ}}}{\mathrm{\Σ}}}\ln \left({\left|{\hat{s}}_l\right|}^2\right)\\ Z=V\hat{S}\end{array}\right.$

Described iteration convergence module is for obtaining

The optimization estimated result, need to guarantee cost function Γ minimum.As can be seen, for the target empty bay matrix S after the weighting ' any S ' _Free, utilize The target empty bay matrix that obtains

Need satisfy equation of constraint, thereby, cost function and equation of constraint separate and can regard one as about S ' _FreeNo constrained minimization find the solution problem.Promptly ask for S ' _Free, satisfy

Have: $\frac{\∂\mathrm{\Γ}}{\∂{\hat{s}}_l}=\frac{2{\hat{s}}_l}{{\left|{\hat{s}}_l\right|}^2},$ Being write as matrix form is: $\frac{\∂\mathrm{\Γ}}{\∂\hat{S}}={2T}^{-1}\hat{S}$

$\&DoubleRightArrow;\frac{\&PartialD;\mathrm{\&Gamma;}}{{\&PartialD;S}_{\mathrm{free}}^{\&prime;}}=P_{\mathrm{free}}^H\frac{\&PartialD;\mathrm{\&Gamma;}}{\&PartialD;\hat{S}}{2P}_{\mathrm{free}}^H{T}^{-1}\hat{S}=0,$ Wherein $T=\mathrm{diag}({\left|{\mathrm{<figure-callout\; id="1"\; label="s"\; filenames="HSA00000406146500041.png,HSA00000406146500042.png"\; state="\{\{state\}\}">s</figure-callout>}}_1\right|}^2,{\left|{\mathrm{<figure-callout\; id="2"\; label="s"\; filenames="HSA00000406146500041.png,HSA00000406146500042.png"\; state="\{\{state\}\}">s</figure-callout>}}_2\right|}^2,\&CenterDot;\&CenterDot;\&CenterDot;,{\left|s_{L_{\mathrm{\&theta;}}}\right|}^2)$

Introduce $\hat{S}={\mathrm{PS}}^{\&prime;};$ $\&DoubleRightArrow;S_{\mathrm{free}}^{\&prime;}=-{\left(P_{\mathrm{free}}^H{T}^{-1}{P}_{\mathrm{free}}\right)}^{-1}\left(P_{\mathrm{free}}^H{T}^{-1}{P}_{\mathrm{fix}}\right)S_{\mathrm{fix}}^{\&prime;}$

$=\left(P_{\mathrm{free}}^H{\mathrm{TP}}_{\mathrm{fix}}\right){\left(P_{\mathrm{fix}}^H{\mathrm{TP}}_{\mathrm{fix}}\right)}^{-1}{S}_{\mathrm{fix}}^{\&prime;}$

Following formula obtains S ' _FreeSeparate, by iterating to S ' _FreeConvergence makes cost function Γ minimum, thereby obtains different azimuth grid place, target empty bay matrix

The effect that the present invention is useful is:

(1) this method need not to adopt covariance matrix to carry out direct or indirect resolution process, is not subject to the snap restricted number, so can handle the fast beat of data of active list.

(2) this method adopts maximum posteriori criterion to estimate to differentiate coherent source, need not relevant the processing separated in strong coherent signal source in water body or the stratum, and anti-coherent interference ability is strong.

(3) this method and CBF, MVDR, MUSIC method thinking are different, do not belong to object detection method and DOA target Bearing Estimation algorithm, the division of signal that same distance can be returned is the source signal sum of abundant orientation grid, utilize cost function and equation of constraint to carry out iteration convergence and guarantee respectively to scan source signal optimum waveform estimation on each orientation, it is stronger to handle reliability, stable performance, 1 ° azimuthal resolution can be reached in theory, thereby the orientation high-resolution imaging of small-sized latent device can be realized seabed stratigraphic distribution formula coherent source.

Description of drawings

Fig. 1: the stratum three-dimensional acoustic imaging technical scheme synoptic diagram that expression the present invention is suitable for

Fig. 2: represent acoustic imaging system platform synoptic diagram of the invention process

Fig. 3: the expression the present invention be suitable for apart from layering and orientation grid division synoptic diagram

Fig. 4: represent operational flowchart of the present invention

Fig. 5: bead CBF azimuth discrimination figure in the seawater

Fig. 6: expression is differentiated bead PAF azimuth discrimination figure in the seawater with the present invention

Fig. 7: expression is handled bead target main lobe width comparison diagram with CBF and the present invention

Fig. 8: the two arrowband coherent source spatial spectrums (SNR=10dB) that expression is estimated with the present invention

Fig. 9: the two arrowband coherent source spatial spectrums (SNR=-5dB) that expression is estimated with the present invention

Figure 10: the iteration error convergence curve of coherent source spatial spectrum is estimated in expression with the present invention.

Embodiment

The invention will be further described below in conjunction with drawings and Examples:

One of invention embodiment checking: invention application implementation platform is a three-dimensional imaging acoustic detection pull-type system in this example, connects by towing cable between the underwater platform (towed body) that sonar transducer array promptly is installed and the lash ship.Its stratum three-dimensional imaging principle design as shown in Figure 1, the present invention needs receiving hydrophone and platform motion direction homeotropic alignment to distribute, as shown in Figure 2.When towed body towed motion, the many array element nautical receiving set collection data that are parallel to direction of motion were done synthetic aperture processing, did phased azimuth filtering perpendicular to the hydrophone array of direction of motion and handled, and part array element data can be used as degree of depth wideband pulse compression and handle.Adopt stratum three-dimensional acoustic imaging system to obtain a bead target and stratum detecting data in the seawater are handled, stratum, seabed echo belongs to strong relevant reverberation.

Platform is in motion process, and each nautical receiving set obtains echoed signal, through gathering, transmits, and filtering is amplified, and UDP network transmission protocol data sharing and real time signal processing data processing are adopted in the decomposition of frame number envelope.

The inventive method is at first carried out the data pre-service in end for process, according to actual needs, the initialization system parameter, comprise signal type and pulsewidth, target depth range, layering is apart from band, the orientation angles scope, the distinguishable angle in minimum orientation is chosen the time samples point that arrives echo data, and the stratum echoed signal of every layer of arrival evenly is divided into L _θThe source signal sum that individual orientation grid arrives, the set of orientation grid are not more than transmitting transducer emission directive property width.Next carries out the array matrix conversion, comprise equation equation step, parameter relationship formula step two parts, by transforming target empty bay matrix, feature space vector, eigenwert, the array manifold matrix after acquisition receives echo space of matrices, the transformed matrix space that receives echo, target empty bay matrix, weighting, set up equation expression formula and parameter relationship formula between the above-mentioned parameter.Estimate echo intensity on the orientation with maximum posteriori criterion (MAP) then,, set up cost function and equation of constraint according to maximum posteriori criterion.According to cost function and equation of constraint, can be converted into no constrained minimization and find the solution problem at last, iterating makes the cost function minimum, obtain the layering distance with on the optimal estimation value of orientation grid echoed signal intensity.

Detailed process step of the present invention is as follows:

(1) the stratum echoed signal evenly is divided into L _θIndividual orientation grid (L _θ＞＞M (array element number)), as shown in Figure 3, s ₁(t), s ₂(t) ...,

Source signal for the different azimuth grid satisfies far field condition,

M ∈ (1,2 ..., M), l ∈ (1 ..., L _θ) be a series of source signal (1 ..., L _θSource signal) each array element (1 ... M array element) frequency-domain expression of the phase differential that produces, n _1...M(t) noise sequence is the zero-mean Gaussian process, and noise is separate between array element, and each sensor of array is an isotropy, the inconsistent interference of no mutual coupling and passage.The equation equation further is expressed as follows:

$Z\left(t\right)=\underset{l=1}{\overset{L_{\mathrm{\&theta;}}}{\mathrm{\&Sigma;}}}v\left({\mathrm{\&phi;}}_l\right)s_l\left(t\right)=V({\mathrm{\&phi;}}_1,...,{\mathrm{\&phi;}}_{L_{\mathrm{\&theta;}}})\&CenterDot;S\left(t\right)+N\left(t\right)\&DoubleRightArrow;Z=V(S+\mathrm{\&Delta;S})=V\hat{S}$

Z (t)=[z wherein ₁(t) ..., z _M(t)] ^TBe the array element received signal vector of M * 1,

Be L _θ* 1 source signal vector,

Be the array manifold vector of M * 1 of l source signal on M array element,

Be M * L _θThe array manifold matrix, N (t)=[n ₁(t) ..., n _M(t)] ^TBe the noise vector of M * 1, For sneaking into the target empty bay matrix after the incoherent noise, φ _l=2 π f _cDcos θ _l/ c, wherein f _cBe centre frequency, d is an array element distance, and c is the velocity of sound in the water.The present invention estimates θ with maximum posteriori criterion (MAP) _lEcho intensity on the orientation

Have this orientation of expression of peak value that the strong scattering target is arranged with the orientation output spectra, and the height of peak value is represented the estimated value to this orientation target scattering intensity.Choose cost function and equation of constraint

(2) the reception echo vector of Huo Deing, array manifold matrix, target empty bay matrix three's equation equation With V ^HV makes feature decomposition, V ^HV=P Λ P ^H, wherein λ _iBe i eigenwert, e _iBe i normalization proper vector.This process can realize by V is carried out svd.Because V is M * L _θMatrix, and L _θ＞＞M, L _θ* L _θMatrix V ^HThe order of V maximum is M, therefore has L at least _θ-M zero eigenvalue.Suppose that eigenwert is with descending sort:

Eigenwert is divided into two groups: nonzero eigenvalue (λ ₁, λ ₂..., λ _Fix) and zero eigenvalue

The number of two stack features values is respectively L _FixAnd L _Free=L _θ-L _Fix

Proper vector P is divided into two submatrixs, P _Fix(L _θ* L _Fix) and P _Free(L _θ* L _Free), correspond respectively to nonzero eigenvalue and zero eigenvalue.

(3) adopt the matrix equivalence to replace, receive the echo vector apart from each array element in the layering and introduce array manifold Matrix Conjugate value and feature space vector conjugation information, the equation equation left side obtains to receive the transformed matrix space Z ' of echo, with the product of equation the right feature space vector conjugation information and target empty bay matrix, the target empty bay matrix S after the weighting after equation equation the right obtains to transform '.Further the matrix equivalence is derived as follows:

$\&DoubleRightArrow;P^H{V}^{H}Z={\mathrm{\&Lambda;P}}^H\hat{S}\&DoubleRightArrow;\mathrm{QZ}={\mathrm{\&Lambda;P}}^H\hat{S}(Q=P^H{V}^H)\&DoubleRightArrow;Z^{\&prime;}={\mathrm{\&Lambda;S}}^{\&prime;}(Z^{\&prime;}=\mathrm{QZ},S^{\&prime;}=P^H\hat{S})$

(4) acquisition is calculated Z ' to the reception data vector Z (t) of each moment t.By above-mentioned conversion, utilize the mutually orthogonal characteristic in nonzero eigenvalue subspace and zero eigenvalue subspace, Z ' _FixAnd S ' _FixBe nonzero eigenvalue subspace, Z ' _FreeAnd S ' _FreeBe the zero eigenvalue subspace, can will receive the transformed matrix space Z ' of echo and the target empty bay matrix S after the weighting ' all carry out the subspace to expand into:

$\left[\begin{array}{c}{Z}_{\mathrm{fix}}^{\&prime;}\\ Z_{\mathrm{free}}^{\&prime;}\end{array}\right]=\left[\begin{array}{cc}{\mathrm{\&Lambda;}}_{\mathrm{fix}}& \\ & {\mathrm{\&Lambda;}}_{\mathrm{free}}\end{array}\right]\left[\begin{array}{c}{S}_{\mathrm{fix}}^{\&prime;}\\ S_{\mathrm{free}}^{\&prime;}\end{array}\right]=\left[\begin{array}{c}{\mathrm{\&Lambda;}}_{\mathrm{fix}}{S}_{\mathrm{fix}}^{\&prime;}\\ {\mathrm{\&Lambda;}}_{\mathrm{free}}{S}_{\mathrm{free}}^{\&prime;}\end{array}\right]$

(5) Z ' of main energy distribution in the subspace of acquisition Z ' and S ' _FixAnd S ' _FixRelational expression, each element can be written as:

$Z_{\mathrm{fix}}^{\&prime;}={\mathrm{\&Lambda;}}_{\mathrm{fix}}{S}_{\mathrm{fix}}^{\&prime;}\&DoubleRightArrow;S_{\mathrm{fix}}^{\&prime;}={{\mathrm{\&Lambda;}}_{\mathrm{fix}}}^{-1}{Z}_{\mathrm{fix}}^{\&prime;}$

$z_{\mathrm{fix},l}^{\&prime;}={\mathrm{\&lambda;}}_l{s}_{\mathrm{fix},l}^{\&prime;}\&DoubleRightArrow;s_{\mathrm{fix},l}^{\&prime;}=\frac{1}{{\mathrm{\&lambda;}}_l}{z}_{\mathrm{fix},l}^{\&prime;},$ l＝1，...，L _fix

Calculate s ' _Fix=Λ ^-1 _FixZ ' _Fix, Λ wherein ^-1 _Fix=diag (λ ₁ ^-1, λ ₂ ^-1..., λ _Fix ^-1).

(6) energy distribution S ' in the target empty bay matrix after the weighting _FixCan obtain by the transformed matrix space Z ' that receives echo, can obtain target empty bay matrix successively.

$S^{\&prime;}=P^H\hat{S}\&DoubleRightArrow;\hat{S}={\mathrm{PS}}^{\&prime;}=\left[\begin{array}{cc}{P}_{\mathrm{fix}}& P_{\mathrm{free}}\end{array}\right]\left[\begin{array}{c}{S}_{\mathrm{fix}}^{\&prime;}\\ S_{\mathrm{free}}^{\&prime;}\end{array}\right]P_{\mathrm{fix}}{S}_{\mathrm{fix}}^{\&prime;}+P_{\mathrm{free}}{S}_{\mathrm{free}}^{\&prime;}$

Initialization vector S ' _Free, all elements zero setting.Calculate

Calculate diagonal matrix

(7) to each reception data vector Z (t) of t constantly, for obtaining

The optimization estimated result, need to guarantee cost function Γ minimum.As can be seen, for the target empty bay matrix S after the weighting ' any S ' _Free, utilize

The target empty bay matrix that obtains

Need satisfy equation of constraint, thereby, cost function and equation of constraint separate and can regard one as about S ' _FreeNo constrained minimization find the solution problem.Promptly ask for S ' _Free, satisfy

Have: $\frac{\&PartialD;\mathrm{\&Gamma;}}{\&PartialD;{\hat{s}}_l}=\frac{2{\hat{s}}_l}{{\left|{\hat{s}}_l\right|}^2},$ Being write as matrix form is: $\frac{\&PartialD;\mathrm{\&Gamma;}}{\&PartialD;\hat{S}}={2T}^{-1}\hat{S}$

$\&DoubleRightArrow;\frac{\&PartialD;\mathrm{\&Gamma;}}{{\&PartialD;S}_{\mathrm{free}}^{\&prime;}}=P_{\mathrm{free}}^H\frac{\&PartialD;\mathrm{\&Gamma;}}{\&PartialD;\hat{S}}{2P}_{\mathrm{free}}^H{T}^{-1}\hat{S}=0,$ Wherein $T=\mathrm{diag}({\left|s_1\right|}^2,{\left|s_2\right|}^2,\&CenterDot;\&CenterDot;\&CenterDot;,{\left|s_{L_{\mathrm{\&theta;}}}\right|}^2)$

$\&DoubleRightArrow;S_{\mathrm{free}}^{\&prime;}=-{\left(P_{\mathrm{free}}^H{T}^{-1}{P}_{\mathrm{free}}\right)}^{-1}\left(P_{\mathrm{free}}^H{T}^{-1}{P}_{\mathrm{fix}}\right)S_{\mathrm{fix}}^{\&prime;}$

$=\left(P_{\mathrm{free}}^H{\mathrm{TP}}_{\mathrm{fix}}\right){\left(P_{\mathrm{fix}}^H{\mathrm{TP}}_{\mathrm{fix}}\right)}^{-1}{S}_{\mathrm{fix}}^{\&prime;}$

Calculate $S_{\mathrm{free}}^{\&prime;}=\left(P_{\mathrm{free}}^H{\mathrm{TP}}_{\mathrm{fix}}\right){\left(P_{\mathrm{fix}}^H{\mathrm{TP}}_{\mathrm{fix}}\right)}^{-1}{S}_{\mathrm{fix}}^{\&prime;}.$

(8) repeat (6) to (7) step, up to S ' _FreeConvergence.Can get at last

Estimator.

Wherein (1) to (3) step is only done one time, and (4) to (8) will distinguish computing apart from layering to each that choose.

Sum up (1)～(8), algorithm flow chart as shown in Figure 4.

Adopt conventional wave beam to form CBF result such as Fig. 5, adopt the inventive method result as shown in Figure 6, horizontal ordinate is an orientation angles, ordinate is the degree of depth, and through target azimuth, the bead focusing effect stratum resolving power of algorithm process of the present invention, picture contrast is all significantly better than the CBF result, with the main lobe 3dB width after handling relatively, as shown in Figure 7,5.1 ° of CBF beam angles, 1.03 ° of PAF beam angles.

Two of invention embodiment checking: employing reception battle array is 9 yuan a uniform straight line array, and array element distance is 0.167 meter, and single snap arrowband coherent signal of two frequency 10kHz is arranged, and pulsewidth is 10ms, and the orientation is respectively 90 ° and 92 °.Fig. 8, shown in Figure 9 be signal to noise ratio (S/N ratio) be 10dB and-during 5dB (signal to noise ratio (S/N ratio) is defined as the ratio of array element received signal power and noise power herein), the normalization spatial spectrum of CBF, MVDR, MUSIC and PAF algorithm relatively, all algorithms have all been done the narrow-band filtering processing, and wherein MVDR and MUSIC do and separate relevant the processing.From figure, verify and draw, CBF can't differentiate two targets, and MVDR and MUSIC are by separating the relevant distinguishable two relevant targets of handling, comparatively speaking, adopt phased azimuth filtering PAF algorithm to have high resolving power and low secondary lobe, and do not need to separate relevant processing and can differentiate coherent source.When signal to noise ratio (S/N ratio) was low, MVDR and MUSIC can not differentiate two targets, but the PAF algorithm still can be differentiated preferably.70 °～110 ° spaces evenly are divided into 201 grids, the array received data of each moment point are carried out interative computation, be averaged behind the independent experiment 512 times and obtain convergence curve as shown in figure 10.

Should illustrate that at last above example is only in order to illustrating technical scheme of the present invention and not limit therewith, but on using, can extend to other modification, change, use, and think the modification that all are such, change, use, embodiment all within the spirit and scope of the present invention.

Claims (2)

1. the phased method of bearing filtering of a three-dimensional formation imaging, it is characterized in that: at same distance, different angles, single snap, heterogeneous dried source signal, the source signal optimum waveform that adopts maximum posteriori criterion to obtain on each scan position of space is estimated, judges target direction according to the energy size of signal on each orientation; Comprise that step is as follows:

1) data preprocessing module: according to actual needs, the initialization system parameter, comprise signal type and pulsewidth, target depth range, layering is apart from band, orientation angles scope, the distinguishable angle in minimum orientation, choose the time samples point that arrives echo data, the stratum echoed signal of every layer of arrival evenly is divided into L _θIndividual orientation grid, the set of orientation grid are not more than transmitting transducer emission directive property width;

2) array matrix conversion module: comprise equation equation step, parameter relationship formula step two parts, by transforming target empty bay matrix, feature space vector, eigenwert, the array manifold matrix after acquisition receives echo space of matrices, the transformed matrix space that receives echo, target empty bay matrix, weighting, set up equation expression formula and parameter relationship formula between the above-mentioned parameter;

3) MAP criterion estimation module: estimate echo intensity on the orientation with maximum posteriori criterion,, set up cost function and equation of constraint according to maximum posteriori criterion;

4) iteration convergence module: according to cost function and equation of constraint, be converted into no constrained minimization and find the solution problem, iterating makes the cost function minimum, obtain the layering distance with on the optimal estimation value of orientation grid echoed signal intensity.

2. the phased method of bearing filtering of three-dimensional formation imaging according to claim 1 is characterized in that:

(1) the stratum echoed signal evenly is divided into L _θIndividual orientation grid L _θ＞＞M array element number, s ₁(t), s ₂(t) ...,

Source signal for the different azimuth grid satisfies far field condition,

M ∈ (1,2 ..., M), l ∈ (1 ..., L _θ) be a series of source signal, 1 ..., L _θSource signal, in each array element, 1 ... M array element, the frequency-domain expression of the phase differential of generation, n _1...M(t) noise sequence is the zero-mean Gaussian process, and noise is separate between array element, and each sensor of array is an isotropy, the inconsistent interference of no mutual coupling and passage; The equation equation further is expressed as follows:

$Z\left(t\right)=\underset{l=1}{\overset{L_{\mathrm{\&theta;}}}{\mathrm{\&Sigma;}}}v\left({\mathrm{\&phi;}}_l\right)s_l\left(t\right)=V({\mathrm{\&phi;}}_1,...,{\mathrm{\&phi;}}_{L_{\mathrm{\&theta;}}})\&CenterDot;S\left(t\right)+N\left(t\right)\&DoubleRightArrow;Z=V(S+\mathrm{\&Delta;S})=V\hat{S}$

Z (t)=[z wherein ₁(t) ..., z _M(t)] ^TBe the array element received signal vector of M * 1,

Be L _θ* 1 source signal vector,

Be the array manifold vector of M * 1 of l source signal on M array element,

Be M * L _θThe array manifold matrix, N (t)=[n ₁(t) ..., n _M(t)] ^TBe the noise vector of M * 1,

For sneaking into the target empty bay matrix after the incoherent noise, φ _l=2 π f _cDcos θ _l/ c, wherein f _cBe centre frequency, d is an array element distance, and c is the velocity of sound in the water, and (MAP) estimates θ with maximum posteriori criterion _lEcho intensity on the orientation

Have this orientation of expression of peak value that the strong scattering target is arranged with the orientation output spectra, and the height of peak value is represented the estimated value to this orientation target scattering intensity, chooses cost function and equation of constraint

(2) the reception echo vector of Huo Deing, array manifold matrix, target empty bay matrix three's equation equation

With V ^HV makes feature decomposition, V ^HV=P Λ P ^H, wherein

λ _iBe i eigenwert, e _iBe i normalization proper vector, this process realizes by V is carried out svd, supposes that eigenwert is with descending sort: Eigenwert is divided into two groups: nonzero eigenvalue (λ ₁, λ ₂..., λ _Fix) and zero eigenvalue

The number of two stack features values is respectively L _FixAnd L _Free=L _θ-L _Fix,

Proper vector P is divided into two submatrixs, P _Fix(L _θ* L _Fix) and P _Free(L _θ* L _Free), correspond respectively to nonzero eigenvalue and zero eigenvalue;

(3) adopt the matrix equivalence to replace, receive the echo vector apart from each array element in the layering and introduce array manifold Matrix Conjugate value and feature space vector conjugation information, the equation equation left side obtains to receive the transformed matrix space Z ' of echo, product with equation the right feature space vector conjugation information and target empty bay matrix, target empty bay matrix S after the weighting after equation equation the right obtains to transform ', further the matrix equivalence is derived as follows:

$\&DoubleRightArrow;P^H{V}^{H}Z={\mathrm{\&Lambda;P}}^H\hat{S}\&DoubleRightArrow;\mathrm{QZ}={\mathrm{\&Lambda;P}}^H\hat{S}(Q=P^H{V}^H)\&DoubleRightArrow;Z^{\&prime;}={\mathrm{\&Lambda;S}}^{\&prime;}(Z^{\&prime;}=\mathrm{QZ},S^{\&prime;}=P^H\hat{S})$

(4) acquisition is calculated Z ' to the reception data vector Z (t) of each moment t, by above-mentioned conversion, utilizes the mutually orthogonal characteristic in nonzero eigenvalue subspace and zero eigenvalue subspace, Z ' _FixAnd S ' _FixBe nonzero eigenvalue subspace, Z ' _FreeAnd S ' _FreeBe the zero eigenvalue subspace, will receive the transformed matrix space Z ' of echo and the target empty bay matrix S after the weighting ' all carry out the subspace to expand into:

$\left[\begin{array}{c}{Z}_{\mathrm{fix}}^{\&prime;}\\ Z_{\mathrm{free}}^{\&prime;}\end{array}\right]=\left[\begin{array}{cc}{\mathrm{\&Lambda;}}_{\mathrm{fix}}& \\ & {\mathrm{\&Lambda;}}_{\mathrm{free}}\end{array}\right]\left[\begin{array}{c}{S}_{\mathrm{fix}}^{\&prime;}\\ S_{\mathrm{free}}^{\&prime;}\end{array}\right]=\left[\begin{array}{c}{\mathrm{\&Lambda;}}_{\mathrm{fix}}{S}_{\mathrm{fix}}^{\&prime;}\\ {\mathrm{\&Lambda;}}_{\mathrm{free}}{S}_{\mathrm{free}}^{\&prime;}\end{array}\right]$

(5) Z ' of main energy distribution in the subspace of acquisition Z ' and S ' _FixAnd S ' _FixRelational expression, each element can be written as:

$Z_{\mathrm{fix}}^{\&prime;}={\mathrm{\&Lambda;}}_{\mathrm{fix}}{S}_{\mathrm{fix}}^{\&prime;}\&DoubleRightArrow;S_{\mathrm{fix}}^{\&prime;}={{\mathrm{\&Lambda;}}_{\mathrm{fix}}}^{-1}{Z}_{\mathrm{fix}}^{\&prime;}$

$z_{\mathrm{fix},l}^{\&prime;}={\mathrm{\&lambda;}}_l{s}_{\mathrm{fix},l}^{\&prime;}\&DoubleRightArrow;s_{\mathrm{fix},l}^{\&prime;}=\frac{1}{{\mathrm{\&lambda;}}_l}{z}_{\mathrm{fix},l}^{\&prime;},$ l＝1，...，L _fix

Calculate s ' _Fix=Λ ^-1 _FixZ ' _Fix, Λ wherein ^-1 _Fix=diag (λ ₁ ^-1, λ ₂ ^-1..., λ _Fix ^-1)

(6) energy distribution S ' in the target empty bay matrix after the weighting _FixObtain by the transformed matrix space Z ' that receives echo, obtain target empty bay matrix successively;

$S^{\&prime;}=P^H\hat{S}\&DoubleRightArrow;\hat{S}={\mathrm{PS}}^{\&prime;}=\left[\begin{array}{cc}{P}_{\mathrm{fix}}& P_{\mathrm{free}}\end{array}\right]\left[\begin{array}{c}{S}_{\mathrm{fix}}^{\&prime;}\\ S_{\mathrm{free}}^{\&prime;}\end{array}\right]P_{\mathrm{fix}}{S}_{\mathrm{fix}}^{\&prime;}+P_{\mathrm{free}}{S}_{\mathrm{free}}^{\&prime;}$

Initialization vector S ' _Free, all elements zero setting is calculated Calculate diagonal matrix

(7) to each reception data vector Z (t) of t constantly, for obtaining

The optimization estimated result, need to guarantee cost function Γ minimum, with cost function and equation of constraint separate and can regard one as about S ' _FreeNo constrained minimization find the solution problem, promptly ask for S ' _Free, satisfy

Have: $\frac{\&PartialD;\mathrm{\&Gamma;}}{\&PartialD;{\hat{s}}_l}=\frac{2{\hat{s}}_l}{{\left|{\hat{s}}_l\right|}^2},$ Being write as matrix form is: $\frac{\&PartialD;\mathrm{\&Gamma;}}{\&PartialD;\hat{S}}={2T}^{-1}\hat{S}$

$\&DoubleRightArrow;\frac{\&PartialD;\mathrm{\&Gamma;}}{{\&PartialD;S}_{\mathrm{free}}^{\&prime;}}=P_{\mathrm{free}}^H\frac{\&PartialD;\mathrm{\&Gamma;}}{\&PartialD;\hat{S}}{2P}_{\mathrm{free}}^H{T}^{-1}\hat{S}=0,$ Wherein $T=\mathrm{diag}({\left|s_1\right|}^2,{\left|s_2\right|}^2,\&CenterDot;\&CenterDot;\&CenterDot;,{\left|s_{L_{\mathrm{\&theta;}}}\right|}^2)$

$S_{\mathrm{free}}^{\&prime;}=-{\left(P_{\mathrm{free}}^H{T}^{-1}{P}_{\mathrm{free}}\right)}^{-1}\left(P_{\mathrm{free}}^H{T}^{-1}{P}_{\mathrm{fix}}\right)S_{\mathrm{fix}}^{\&prime;}$

$\&DoubleRightArrow;=\left(P_{\mathrm{free}}^H{\mathrm{TP}}_{\mathrm{fix}}\right){\left(P_{\mathrm{fix}}^H{\mathrm{TP}}_{\mathrm{fix}}\right)}^{-1}{S}_{\mathrm{fix}}^{\&prime;}$

Calculate $S_{\mathrm{free}}^{\&prime;}=\left(P_{\mathrm{free}}^H{\mathrm{TP}}_{\mathrm{fix}}\right){\left(P_{\mathrm{fix}}^H{\mathrm{TP}}_{\mathrm{fix}}\right)}^{-1}{S}_{\mathrm{fix}}^{\&prime;};$

(8) repeat (6) to (7) step, up to S ' _FreeConvergence, at last

Estimator; Wherein (1) to (3) step is only done one time, and (4) to (8) will distinguish computing apart from layering to each that choose.

Images