CN109541572A - A kind of subspace direction estimation method based on linear environmental noise model - Google Patents

Abstract

The present invention relates to a kind of subspace direction estimation methods based on linear environmental noise model, for actual noise circumstance, due to the influence of eolian noise and boat noise, its spatial noise intensity distribution is heterogeneous, it is modeled using line noise model, spatial noise power distribution function is subjected to Fourier series expansion, utilizes finite term Fourier series approximate fits environmental noise model.The signal subspace under line noise Model Condition is extracted, is estimated using the orientation that the orthogonality of signal subspace and noise subspace obtains target.

Description

A kind of subspace direction estimation method based on linear environmental noise model

Technical field

The invention belongs to the fields such as array signal processing, Signal processing of sonar, are related to a kind of based on linear environment noise mode The subspace direction estimation method of type.

Background technique

Orientation estimation (Direction of Arrival, DOA) based on sensor array obtains in fields such as sonar, radars To being widely applied.Carrying out high-resolution orientation estimation to underwater and water surface multiple target using horizontal positioned hydrophone array is The important content of underwater arrays signal processing, and carry out the key technology of target acquisition and positioning in water.

Sonar array is placed in marine environment, and there are also ambient sea noises other than echo signal for array received signal. Ambient sea noise is the ambient field of ocean acoustic field, can reduce or inhibit ambient noise pair to the research mesh of ambient sea noise The interference of underwater equipment improves the performance and operating distance of underwater equipment.There are many factor of ambient sea noise, usually wrap Include pressure fluctuation caused by pressure wave caused by damp sunset, wave and turbulent flow and seismic activity, pneumatic sea, rainfall, molecule The activity etc. of biocenose in warm-up movement, ocean.Current most of direction estimation methods assume that ambient noise is Gauss white noise Sound, i.e. noise covariance matrix can be write as the product form of noise variance and unit matrix.However, practical each array element is received Ambient sea noise may be relevant, and due to the influence of the factors such as eolian noise and distant place boat noise, ambient noise The spatial distribution of intensity may also have directionality, be the slowly varying function with angle change in the horizontal direction.Therefore, exist Under complicated marine environmental conditions, the noise model that white noise is assumed may cause serious orientation estimated bias.

Current main direction estimation method has conventional beamformer algorithm and the undistorted beamforming algorithm of minimum variance Deng.Conventional beamformer algorithm is the scanning direction-finding method based on conventional beamformer, with the small spy of structure simple computation amount Point, and environmental noise model error has certain robustness, however limited by " Rayleigh criterion ", dimensional orientation is differentiated Ability is poor；Minimum variance is undistorted, and beamforming algorithm utilizes CAPON adaptive beam former realization spacescan wave beam Output power Power estimation, dimensional orientation resolution capability is better than conventional beamformer algorithm, however algorithm hypothesis noise model is White Noise Model.S.D.Somasundaram has studied the steady minimum variance undistorted wave beamformer design method in broadband (Somasundaram S D.Wideband Robust Capon Beamforming for Passive Sonar.IEEE Journal of Oceanic Engineering, 2013,38 (2): 308-322.), and it is applied to the more of passive sonar Target acquisition improves the spatial resolving power of target, but its noise model is it will again be assumed that be White Noise Model, for complex environment The orientation estimation performance of noise conditions still needs to examine.

Summary of the invention

Technical problems to be solved

In order to avoid the shortcomings of the prior art, the present invention proposes a kind of subspace based on linear environmental noise model Direction estimation method can realize height for the horizontal fixed hydrophone array laid under complicated ocean ambient noise conditions Differentiate the estimation of submarine target orientation.

Technical solution

A kind of subspace direction estimation method based on linear environmental noise model, it is characterised in that: using horizontal positioned In underwater rectilinear transducer basic matrix, array number M, array element spacing is d, and K number of targets is distributed in the far field of the horizontal array In range, estimating step is as follows:

Step 1:x (n), n=1 ..., N are the output vector of rectilinear transducer basic matrix column, and computing array output vector is adopted Sample covariance matrix

Wherein: N is data points, and subscript " H " is conjugate transposition operator；

Step 2, estimation model order: first in no undersea detection target K-0, conventional beamformer azimuth spectrum is calculated Figure:

Wherein: a (θ) is array manifold vector a (θ)=[1, exp (- j2 π f τ) ..., exp {-j2 π f (M-1) τ }]^T, f is Signal frequency, τ=dsin θ/c, d are array element spacing, and c is the velocity of sound, and θ is dimensional orientation angle；

Fourier series approaching is carried out to the azimuth spectrum CBF (θ) of conventional beamformer output, fitting order is ring The modeling order of border noise, is denoted as J；

Step 3: computation model parameter matrix

Γ=[vec (Σ₀),vec(Σ₁),…,vec(Σ_J)]

Wherein:

J=2L+1,

Step 4: the initial estimation of noise covariance matrix Σ isWherein in subscript bracket Digital representation the number of iterations；

Step 5: utilizing the noise covariance matrix Σ of (i-1)-th iteration^(i-1)At the noise prewhitening for calculating i-th iteration Covariance matrix after reason

Step 6: rightFeature decomposition is carried out, whereinIt is the dominant eigenvalue matrix that dimension is K × K,Based on it is special Levy the corresponding feature vector of vector, that is, dominant eigenvalue；The dominant eigenvalue matrix is that diagonal entry is numerical value biggish preceding K Characteristic value；

Step 7: reconstruction signal subspace matricesWherein I is unit diagonal matrix；

Step 8: vector quantities operation is carried out for signal covariance matrixSampling is assisted simultaneously Variance matrixVector quantities operation is carried out, i.e.,Matrix-vector operation vec { } refer to matrix column vector by Column combination, forms a long column vector

Step 9: the coefficient vector η of line noise model when calculating i-th iteration⁽ⁱ⁾, Pinv { } indicates pseudo-inverse operation operator；

Step 10: calculating noise covariance matrix Σ when i-th iteration⁽ⁱ⁾=vec^-1{Γη⁽ⁱ⁾, wherein vec^-1{ } table Show the inverse operation of matrix-vector；

Step 11: judging stopping criterion for iteration for abs (L⁽ⁱ⁾-L^(i-1))≤10^-3Whether true, establishment then carries out in next step Suddenly, invalid then return step 4, and iteration serial number i=i+1；

Wherein cost function| | representing matrix The value of determinant, tr { } representing matrix ask mark operation, and signed magnitude arithmetic(al) is asked in abs () expression；

Step 12: obtaining final reconstruction signal subspace matricesSubscript I indicates final the number of iterations；

Step 13: using the orthogonality of signal subspace and noise subspace, the azimuth spectrum for obtaining target is P (θ)=1/ {a^H(θ) Π a (θ) }, wherein

The maximum value of azimuth spectrum is the orientation of target.

Beneficial effect

A kind of subspace direction estimation method based on linear environmental noise model proposed by the present invention, makes an uproar for actual Acoustic environment, due to the influence of eolian noise and boat noise, spatial noise intensity distribution be it is heterogeneous, use line noise Model is modeled, and spatial noise power distribution function is carried out Fourier series expansion, utilizes finite term Fourier series Approximate fits environmental noise model.Extract the signal subspace under line noise Model Condition, using signal subspace with make an uproar The orthogonality in phonon space obtains the orientation estimation of target.

Detailed description of the invention

Fig. 1 is a kind of overall procedure block diagram of subspace direction estimation method based on linear environmental noise model；

Fig. 2 is orientation spectrogram.

Specific embodiment

Now in conjunction with embodiment, attached drawing, the invention will be further described:

The technical solution adopted by the present invention to solve the technical problems includes the following aspects:

1: one array number for being placed horizontally at underwater rectilinear transducer basic matrix is M, and array element spacing is d, at this time array The total length d (M-1) that is；

Underwater or water surface multiple target orientation is carried out using the array to estimate, if number of targets is K, and is distributed in the level In the far-field range of array, incident orientation is θ=(θ₁,…,θ_K), wherein θ_kIndicate that the incident orientation of k-th of signal, array are defeated Outgoing vector is expressed as x (n)=As (n)+e (n), n=1 ..., N, wherein N is data points, s (n)=[s₁(n),…,s_K (n)]^TFor incoming signal waveform, wherein subscript " T " indicates the operation of vector transposition, and e (n) is that array element receives noise, and A is array stream Shape matrix (steering matrix).For horizontal positioned sensor array, array manifold vector a (θ)=[1, exp (- j2 πfτ),…,exp{-j2πf(M-1)τ}]^T, wherein f is signal frequency, and τ=dsin θ/c, d are array element spacing, and c is the velocity of sound.Number It is R=E { x (n) x according to covariance matrix^H(n) }=APA^H+ Σ, wherein E { } indicates mathematic expectaion operator, and P is diagonal matrix, right Diagonal element is signal power, and subscript " H " is conjugate transposition operator, and Σ is noise covariance matrix, i.e. Σ=E { e (n) e^H (n) } sample covariance matrix, is used in practical applicationIndicate approximate data covariance matrix R, i.e.,

2: the foundation of line noise model

In given sampling time n, noise intensity is a stochastic variable v (θ, n).Therefore, array element receives noise waveform ForAssuming that noise intensity v (θ, n) is spatially uncorrelated, it is zero-mean gaussian white noise on the time Sound, then noise covariance matrix is expressed asWherein ε (θ) is spatial noise power density, empty Between Carrier To Noise Power Density ε (θ) be incident orientation angle θ periodic function, 2 π be a cycle, by spatial noise power density ε (θ) It is obtained according to Fourier series expansionWherein Fourier coefficientL is in formula Fourier series expansion order, c_lAnd s_lFourier coefficient under corresponding different rank；

Using the Fourier series of infinite order to spatial noise power density ε (θ) approximate representation, if the order is L, Then data covariance matrixWherein noise covariance square Battle array beWherein J=2L+1, WithAs L=0And work as battle array When first spacing is half-wavelength, Σ=c₀I, wherein I is unit diagonal matrix, and ambient noise is equal to zero-mean gaussian white noise at this time Sound.To noise covariance matrixCarrying out matrix-vector operation, (matrix-vector operation vec { } refers to square The column vector of battle array forms a long column vector by column combination) vec { Σ }=Γ η, matrix Γ=[vec in formula can be obtained (Σ₀),vec(Σ₁),…,vec(Σ_J)], vector η=[η₁,…,η_J]。

3: signal subspace extracts

Signal subspace is extracted by iterative manner.The initial estimation for providing noise covariance matrix Σ first isThe wherein digital representation the number of iterations in subscript bracket.It is assisted using the noise of (i-1)-th iteration Variance matrix Σ^(i-1)It calculates the relevant parameter of i-th iteration: 1. calculating the covariance matrix after noise pre -whitening processing2. rightFeature decomposition is carried out, whereinIt is the dominant eigenvalue square that dimension is K × K Battle array (diagonal entry is the biggish preceding K characteristic value of numerical value),For main feature vector (i.e. the corresponding feature of dominant eigenvalue to Amount).3. reconstruction signal subspace matrices4. carrying out vector for signal covariance matrix Change operationSimultaneously to sample covariance matrixVector quantities operation is carried out, i.e.,⑤ The coefficient vector of line noise model when calculating i-th iterationWhereinFor line noise model system Number, J is model order,Pinv { } indicates pseudo-inverse operation operator.6. when calculating i-th iteration Noise covariance matrix Σ⁽ⁱ⁾=vec^-1{Γη⁽ⁱ⁾, wherein vec^-1The inverse operation of { } representing matrix vectorization.7. judging iteration Termination condition is abs (L⁽ⁱ⁾-L^(i-1))≤10^-3Whether true, establishment then carries out next step, invalid then to return to the 1. step behaviour Make, and iteration serial number i=i+1.Wherein cost function:

| | the value of the determinant of representing matrix, tr { } representing matrix ask mark operation, and signed magnitude arithmetic(al) is asked in abs () expression.⑧ Final reconstruction signal subspace matrices areSubscript I indicates final the number of iterations.

4: target Bearing Estimation

Using the orthogonality of signal subspace and noise subspace, the azimuth spectrum of target is represented by P (θ)=1/ { a^H (θ) Π a (θ) }, whereinThe maximum value for searching for azimuth spectrum, obtains the side of target Position.

Specific steps are as follows:

(1) arrange that one is placed horizontally at underwater rectilinear transducer basic matrix, array number M, array element spacing is d.It utilizes The array received and record underwater sound signal, x (n) are array output vector；

(2) computing array sample covariance matrix

(3) estimate model order.First in no undersea detection target, conventional beamformer orientation spectrogram is calculated, i.e.,Wherein a (θ) be array manifold vector a (θ)=[1, exp (- j2 π f τ) ..., exp {-j2 π f (M-1)τ}]^T, f is signal frequency, and τ=dsin θ/c, d are array element spacing, and c is the velocity of sound, and θ is dimensional orientation angle.To conventional wave The azimuth spectrum that beam forms output carries out Fourier series approaching, and fitting order is the modeling order of ambient noise, is denoted as J；

(4) computation model parameter matrix Γ=[vec (Σ₀),vec(Σ₁),…,vec(Σ_J)], wherein J=2L+1,

(5) initial estimation of noise covariance matrix Σ isThe wherein number in subscript bracket Word indicates the number of iterations；

(6) the noise covariance matrix Σ of (i-1)-th iteration is utilized^(i-1)Calculate the noise pre -whitening processing of i-th iteration Covariance matrix afterwards

(7) rightFeature decomposition is carried out, whereinIt is that (diagonal entry is number for dominant eigenvalue matrix that dimension is K × K It is worth biggish preceding K characteristic value),For main feature vector (i.e. the corresponding feature vector of dominant eigenvalue)；

(8) reconstruction signal subspace matrices

(9) vector quantities operation is carried out for signal covariance matrixSimultaneously to sampling association side Poor matrixVector quantities operation is carried out, i.e.,

(10) coefficient vector of line noise model when calculating i-th iterationWhereinFor line Property noise model coefficient, J is model order,Pinv { } indicates pseudo-inverse operation operator；

(11) noise covariance matrix Σ when i-th iteration is calculated⁽ⁱ⁾=vec^-1{Γη⁽ⁱ⁾, wherein vec^-1{ } indicates The inverse operation of matrix-vector；

(12) judge stopping criterion for iteration for abs (L⁽ⁱ⁾-L^(i-1))≤10^-3Whether true, establishment then carries out next step, Invalid then return step (5), and iteration serial number i=i+1.Wherein cost function| | the value of the determinant of representing matrix, tr { } are indicated Signed magnitude arithmetic(al) is asked in Matrix Calculating mark operation, abs () expression.

(13) final reconstruction signal subspace matrices are exportedSubscript I indicates final the number of iterations；

(14) orthogonality of signal subspace and noise subspace is utilized, the azimuth spectrum of target may be expressed as:

P (θ)=1/ { a^H(θ) Π a (θ) }, wherein

The maximum value for searching for azimuth spectrum, obtains the orientation of target.

Specific embodiment:

Three, far field source signal is received using 10 yuan of hydrophone array, aspect is 55 °, 80 ° and 90 °, battle array First spacing is half-wavelength, and ambient noise is spatial non-uniform noise environment, uses conventional beamformer method, the side MUSIC respectively The method that method and the present invention provide calculates the orientation spectrogram of target, and as shown in Fig. 2, wherein MUSIC method is to noise circumstance More sensitive, no longer valid in such circumstances, conventional beamformer method resolution ratio is lower, and the method that the present invention provides can be with The orientation of three targets is accurately provided, and the influence for spatial non-uniform noise environment is most steady.

Claims (1)

1. a kind of subspace direction estimation method based on linear environmental noise model, it is characterised in that: use and be placed horizontally at Underwater rectilinear transducer basic matrix, array number M, array element spacing are d, and K number of targets is distributed in the far field model of the horizontal array In enclosing, estimating step is as follows:

Step 1:x (n), n=1 ..., N are the output vector of rectilinear transducer basic matrix column, the sampling association of computing array output vector Variance matrix

Wherein: N is data points, and subscript " H " is conjugate transposition operator；

Step 2, estimation model order: first in no undersea detection target K-0, conventional beamformer orientation spectrogram is calculated:

Wherein: a (θ) is array manifold vector a (θ)=[1, exp (- j2 π f τ) ..., exp {-j2 π f (M-1) τ }]^T, f is signal Frequency, τ=dsin θ/c, d are array element spacing, and c is the velocity of sound, and θ is dimensional orientation angle；

Fourier series approaching is carried out to the azimuth spectrum CBF (θ) of conventional beamformer output, fitting order is that environment is made an uproar The modeling order of sound, is denoted as J；

Step 3: computation model parameter matrix

Γ=[vec (Σ₀),vec(Σ₁),…,vec(Σ_J)]

Wherein:

J=2L+1,

Step 4: the initial estimation of noise covariance matrix Σ isThe wherein number in subscript bracket Word indicates the number of iterations；

Step 5: utilizing the noise covariance matrix Σ of (i-1)-th iteration^(i-1)After the noise pre -whitening processing for calculating i-th iteration Covariance matrix

Step 6: rightFeature decomposition is carried out, whereinIt is the dominant eigenvalue matrix that dimension is K × K,For main feature to Amount is the corresponding feature vector of dominant eigenvalue；The dominant eigenvalue matrix is that diagonal entry is the biggish preceding K feature of numerical value Value；

Step 7: reconstruction signal subspace matricesWherein I is unit diagonal matrix；

Step 8: vector quantities operation is carried out for signal covariance matrixSimultaneously to sampling covariance MatrixVector quantities operation is carried out, i.e.,Matrix-vector operation vec { } refers to matrix column vector by column group It closes, forms a long column vector

Step 9: the coefficient vector η of line noise model when calculating i-th iteration⁽ⁱ⁾,pinv{} Indicate pseudo-inverse operation operator；

Step 10: calculating noise covariance matrix Σ when i-th iteration⁽ⁱ⁾=vec^-1{Γη⁽ⁱ⁾, wherein vec^-1{ } indicates square The inverse operation of battle array vectorization；

Step 11: judging stopping criterion for iteration for abs (L⁽ⁱ⁾-L^(i-1))≤10^-3Whether true, establishment then carries out next step, no Set up then return step 4, and iteration serial number i=i+1；

Wherein cost function| | the ranks of representing matrix The value of formula, tr { } representing matrix ask mark operation, and signed magnitude arithmetic(al) is asked in abs () expression；

Step 12: obtaining final reconstruction signal subspace matricesSubscript I indicates final the number of iterations；

Step 13: using the orthogonality of signal subspace and noise subspace, the azimuth spectrum for obtaining target is P (θ)=1/ { a^H(θ) Π a (θ) }, wherein

The maximum value of azimuth spectrum is the orientation of target.