CN105652264A - High-order cumulant-based method for multi-path propagation and separation of acoustic signals - Google Patents

Abstract

The invention discloses a high-order cumulant-based method for the multi-path propagation and separation of acoustic signals. According to the technical scheme of the invention, an original second-order cumulant-based active broadband signal separation method is extended to a fourth-order cumulant. According to the fourth-order cumulant, a corresponding duplicate vector and noise sub-space is constructed. Meanwhile, an additive noise is assumed to be a colored noise during the acoustic signal propagation process. Finally, an acoustic ray path is found out based on the orthogonality of the noise sub-space. The invention also discloses an acoustic ray propagation time tomography-based ocean acoustic tomography method and a sound source localization method. Compared with the prior art, the above method is higher in separation accuracy and can correctly separate acoustic ray paths of smaller intervals. Meanwhile, the method can be applied to conditions wherein the number of sensors is smaller than the number of ray paths. At the same time, the influence of both Gaussian noises and non-Gaussian noises can be inhibited.

Description

Multipath propagation acoustic signal based on Higher Order Cumulants

Technical field

The present invention relates to Underwater Acoustic channels technical field, particularly relate to a kind of multipath propagation acoustic signal based on Higher Order Cumulants.

Background technology

Underwater Acoustic channels technology has been widely used at present in many fields such as Underwater Detection, communication and earthquake, biomedical engineering etc. Owing to acoustical signal propagating in media as well exists with the form of multipath often, what sensor received is the aliased data of the ripple of these multiple propagated, therefore, the premise that acoustical signal to carry out subsequent treatment is to be separated in every sound ray path from the multipath propagation acoustical signal that sensor receives, and this Underwater Acoustic channels method is multipath propagation acoustic signal.

For marine acoustic tomography, Ocean Acoustic Tomography is one of most important technology by velocity of sound method detection ocean interior, and it utilizes sound change of spread speed in ocean to come change or even the marine climate change of inverting ocean temperature. According to different inverse models, can be divided into: chromatography (2) normal-mode propagation time in (1) sound ray propagation time chromatography; (3) peak value coupling chromatography; (4) normal mode phase place chromatography; (5) normal mode horizontal refraction chromatography; (6) Matched Field chromatography.

The forward problem of Shallow Sea Acoustic tomography needs identify different ray paths, be used for solving inversion problem and estimating sonic velocity change by the time of advent of ray path afterwards. High-quality refutation process is based upon (1) discernible thread path number; (2) the overlayable Oceanic waveguide uniform space of thread path. The important properties that Shallow Sea Acoustic tomography utilizes is acoustical signal multipath propagation character in water, it may be assumed that due to refraction or the reflection in sea or seabed, and acoustical signal is the duplicate that form is propagated and each thread path is transmitting signal with multipath. Fig. 1 is an example of multipath propagation. Owing to each thread path covers different sea areas, multipath propagation character can provide more information for the refutation process of Shallow Sea Acoustic tomography. But multipath propagation produces interference region simultaneously, therefore, first needing with array-processing techniques thread path to be separated and estimate the time of advent of different thread path in Shallow Sea Acoustic tomography, this is to obtain the matter of utmost importance that accurate inversion result to solve.

In Shallow Sea Acoustic chromatography is applied, existing array-processing techniques is mainly Wave beam forming class algorithm. Wave beam forming is algorithm simple and classical in array-processing techniques. Originally, it is used in the test structure of point-to-point, it may be assumed that a signal source and a sensor. In order to improve separation accuracy, Wave beam forming is extended to again a little in the structure of array, namely with a vertical sensor array as signal receiving end, the direction of arrival of signal to sensor array has been successfully separated, as discriminant parameter, the situation that part cannot separate in point-to-point topology. While it is true, when thread path arrives array with only small interval, Wave beam forming still can not successfully be isolated identification at point in array structure. Recently, researcher is had to propose a kind of new dualbeam formation algorithm for Signal separator under sound tomography background, this algorithm foundation sound is in the reciprocity of water transmission, signal source array is adopted at signal transmitting terminal, set up the array two-dimentional test structure to array, thus introducing new discriminant parameter and launch angle, test shows to the method increase separation accuracy than beamforming algorithm. For improving algorithm resolving power further, Jiang et al. is in " Raypathseparationwithhighresolutionprocessing " literary composition, propose in point to array structure, in conjunction with spatial-frequency domain smoothing method and actively broadband multi-signals separation algorithm, it is proposed that smooth active broadband multi-signals separation algorithm (smoothingMultipleSignalClassificationActiveLarge-band:sm oothing-MUSICAL). Mainly comprising the following steps of the method: take turns doing spatial-frequency domain smoothing processing after the data Fourier transformation that sensor array is received, seek the covariance matrix of smooth rear data, covariance matrix is done feature decomposition, structure copy vector, structure estimator, finally search out the time of advent and the angle in sound ray path. The advantage that the method has following tripartite face: (1) improves resolving power than Beamforming Method: when the thread path interval time of advent is very little, smoothing-MUSICAL can be successfully separated part signal; (2) the high-resolution separation problem of signal when thread path is correlated with or is relevant is solved; (3) generated data and actual tests result show that smoothing-MUSICAL algorithm has good noise immunity.

But the algorithm in smoothing-MUSICAL still suffers from problems with: (1) resolving power still can not meet the requirement of practical application; (2) algorithm is based on second-order statistic, need to assume that signal is gaussian signal; (3) array number must more than signal number to be separated; (4) pre-treatment step need to be increased and separate relevant or coherent signal.

Summary of the invention

The technical problem to be solved is in that to overcome the deficiency of existing smooth active broadband multi-signals separation algorithm, a kind of multipath propagation acoustic signal based on Higher Order Cumulants is provided, it has higher separation accuracy, the sound ray path that interval is less can be properly separated out, and can be suitably used for the number of sensors situation less than ray path, may also suppress the impact of Gaussian noise and non-Gaussian noise simultaneously.

The present invention specifically solves above-mentioned technical problem by the following technical solutions:

Based on the multipath propagation acoustic signal of Higher Order Cumulants, for isolating every sound ray path from the multipath propagation acoustical signal received by acoustic sensor array, comprise the following steps:

Step 1, the laggard line space territory of sensor received acoustical signal Fourier transformation-frequency domain smoothing is processed, it is thus achieved that K=(2K_s+1)(2K_f+ 1) composition matrix estimated by individual arrowbandx _{Ks, kf}, K_s��K_fRespectively spatial domain smooth order, frequency domain smoothing order;And K is more than sound source number P;

Step 2, utilize following formula calculate matrixx _{Ks, kf}Fourth order cumulant

$\begin{array}{c}\hat{C}=E\left\{({\underset{\‾}{x}}_{ks,kf}\⊗{{\underset{\‾}{x}}_{ks,kf}}^{*}){({\underset{\‾}{x}}_{ks,kf}\⊗{{\underset{\‾}{x}}_{ks,kf}}^{*})}^H\right\}-E\{{\underset{\‾}{x}}_{ks,kf}\⊗{{\underset{\‾}{x}}_{ks,kf}}^{*}\}E\left\{\right({\underset{\‾}{x}}_{ks,kf}\⊗\\ {{\underset{\‾}{x}}_{ks,kf}}^{*}{)}^H\}-E\{\left({\underset{\‾}{x}}_{ks,kf}{{\underset{\‾}{x}}_{ks,kf}}^H\right\}\⊗E\left\{{\left({\underset{\‾}{x}}_{ks,kf}{{\underset{\‾}{x}}_{ks,kf}}^H\right)}^{*}\right\}\end{array}$

Wherein, E represent ask expectation, * to represent asking conjugation, H to represent seeks conjugate transpose,Represent Kronecker product;

Step 3, to fourth order cumulantCarry out EVD feature decomposition, and by obtained (MF)²Individual eigenvalue is arranged as from big to small:Wherein, the frequency number that F chooses when being frequency domain smoothing; M is the number of sensors in acoustic sensor array;

Step 4, rightObtain (MF) after carrying out feature decomposition²Individual eigenvalueWith wherein P²Individual bigger eigenvalueCharacteristic of correspondence vector structure signal subspace, with remaining (MF)²-P²Individual less eigenvalueCharacteristic of correspondence vector structure noise subspace;

Step 5, for copy vector a (��, T) below sound ray path configuration:

$a(\mathrm{\θ},T)=b(\mathrm{\θ},T)\⊗b{(\mathrm{\θ},T)}^{*}$

$b(\mathrm{\θ},T)=\left[\begin{array}{c}e\left(v_1\right)e^{-2{\mathrm{i\πv}}_{1}T}d(v_1,\mathrm{\θ})\\ \mathrm{...}\\ e\left(v_F\right)e^{-2{\mathrm{i\πv}}_{F}T}d(v_F,\mathrm{\θ})\end{array}\right]$

$d(v_i,\mathrm{\θ})=\[1,e^{-2{\mathrm{i\πv}}_i{\mathrm{\τ}}_{1,2}\left(\mathrm{\θ}\right)},...,e^{-2{\mathrm{i\πv}}_i{\mathrm{\τ}}_{1,M-1}\left(\mathrm{\θ}\right)}\]$

Wherein, �� represents the angle of arrival in sound ray path, and T represents the time of advent in sound ray path, e (v_i) it is that acoustical signal is at frequency v_iThe amplitude at place, i=1,2 ..., F,Represent Kronecker product, ��_1,j(��) represent that sound ray path arrives jth sensor relative to arriving the time delay of the 1st sensor as reference sensor, j=2 ..., M-1;

Step 6, the search copy vector orthogonal with described noise subspace, these copy vectors are the sound ray path being finally recovered out, and its �� and T is the corresponding arrival direction in sound ray path and the time of advent.

Techniques below scheme can also be obtained according to identical invention thinking:

A kind of Ocean Acoustic Tomography method based on sound ray propagation time chromatography, sound change of spread speed in ocean is utilized to carry out inverting ocean environment parameter, it is characterized in that, first with the above-mentioned multipath propagation acoustic signal based on Higher Order Cumulants, from the multipath propagation acoustical signal received by acoustic sensor array, isolate each sound ray path; Then it is finally inversed by ocean environment parameter the time of advent according to sound ray path.

A kind of sound localization method, first with the above-mentioned multipath propagation acoustic signal based on Higher Order Cumulants, isolates each sound ray path from the multipath propagation acoustical signal received by acoustic sensor array; Then according to the angle of arrival in sound ray path and the orientation determining sound source the time of advent in the two dimensional surface spatially and temporally constituted.

Compared to existing technology, the method have the advantages that

1, the present invention has higher sound ray path separation accuracy;

2, the present invention can properly separate out the sound ray path that interval is less;

3, the present invention can properly separate, when number of sensors is less than ray path, thread path of speaking;

4, for the modeling of acoustic signal propagation, existing smooth active broadband multi-signals separation algorithm assumes that additive noise is Gaussian noise, but practical situation is often really not so, therefore can cause signal model mismatch; And present invention assumes that additive noise is chromatic noise, assume than traditional Gaussian noise the marine environment situation more meeting reality can suppress all kinds of environment noise (including Gaussian noise and non-Gaussian noise) better;

5, the artifact that the present invention produces is less.

Accompanying drawing explanation

Fig. 1 is an example of acoustical signal multipath propagation;

Fig. 2 is an example of spatial domain smoothing method;

Fig. 3 is an example of frequency domain smoothing method;

Fig. 4 is the Contrast on effect of the inventive method and existing smoothing-MUSICAL method.

Detailed description of the invention

Below in conjunction with accompanying drawing, technical scheme is described in detail:

At point in array structure, according to the propagation in media as well of sound ray path, by parameters such as direction of arrival, set up the mathematical model of the physical process that can rationally describe acoustical signal multipath propagation.The number assuming signal source is P, and the number of sensor is M, then the signal received on m-th sensor is represented by:

$x_m\left(t\right)=\underset{p=1}{\overset{P}{\Σ}}{a}_{p}e(t-{\mathrm{\τ}}_{m,p})+n_m\left(t\right)---\left(1\right)$

��_{M, p}=T_p+t_m(��_p)(2)

(1) formula is done Fourier transformation and is represented by:

$x_m\left(v\right)=\underset{p=1}{\overset{P}{\Σ}}{a}_{p}e\left(v\right)-\exp (-j2{\mathrm{\π\τ}}_{m,p})+n_m\left(v\right)---\left(3\right)$

On above signal model basis, multipath propagation acoustic signal of the present invention specifically includes following steps: step 1, the received acoustical signal of sensor carries out spatial domain-frequency domain smoothing process, it is thus achieved that K=(2K_s+1)(2K_f+ 1) composition matrix estimated by individual arrowbandx _{Ks, kf}, K_s��K_fRespectively spatial domain smooth order, frequency domain smoothing order; And K is more than sound source number P.

Sound ray path is by the refraction through extended media of the signal source or the multipath propagation ray that reflects to form, it is concerned with completely in these sound ray paths, the Cross-spectral Moment rank of matrix of signal is 1, but only when Cross-spectral Moment rank of matrix is at least sound source number P, just sound ray path can be efficiently separated, therefore, in order to increase the quantity of rank of matrix, it is necessary to use spatial domain-frequency domain smoothing technology. Spatial domain-frequency domain smoothing is prior art, below it is briefly described.

Spatial domain smoothly through on average realizing spatial domain subarray, its principle such as Fig. 2. M sensor is divided into some equivalently-sized, partly overlapping subarray, it is assumed that subarray is linearly consistent, then on subarray sensor, acute variation will not occur the intensity of signal. Group matrix number is non-single order matrix more than or equal to ray path number time space spectrum matrix. The smooth order of hypothesis space is K_s, then each subarray is of a size of M-2K_s, the number of subarray is 2K_s+1��

Frequency domain smoothing method according to operation to as if time domain data or frequency domain data can be divided into the following two kinds: (1) weighted correlation matrix; (2) frequency domain sub-band is average. The present invention adopts frequency domain sub-band averaging method, and its concrete principle is shown in Fig. 3. The smooth order assuming frequency domain is K_f, then 2K can be obtained after smoothing_f+ 1 is of a size of M-2K_fSubband.

The excess smoothness of spatial domain or frequency domain all can be caused that cross-spectrum matrix produces serious error, in order to reduce estimation difference, it is necessary to combined in spatial domain and the smooth of frequency domain.

From one group of observing matrix x, it is possible to produce 2K_sThe observation data that+1 spatial translation occursx _ks. This 2K_sThe matrix of+1 reproduction can obtain again K=(2K respectively through the smooth of frequency domain_s+1)(2K_f+ 1) composition matrix estimated by individual arrowbandx _{Ks, kf}. In order to better ray path be separated from noise, it is necessary to arrange the value of K more than P.

Step 2, calculating matrixx _{Ks, kf}Fourth order cumulant

Fourth order cumulantComputing formula as follows:

$\begin{array}{c}\hat{C}=E\left\{({\underset{\‾}{x}}_{ks,kf}\⊗{{\underset{\‾}{x}}_{ks,kf}}^{*}){({\underset{\‾}{x}}_{ks,kf}\⊗{{\underset{\‾}{x}}_{ks,kf}}^{*})}^H\right\}-E\{{\underset{\‾}{x}}_{ks,kf}\⊗{{\underset{\‾}{x}}_{ks,kf}}^{*}\}E\left\{\right({\underset{\‾}{x}}_{ks,kf}\⊗\\ {{\underset{\‾}{x}}_{ks,kf}}^{*}{)}^H\}-E\{\left({\underset{\‾}{x}}_{ks,kf}{{\underset{\‾}{x}}_{ks,kf}}^H\right\}\⊗E\left\{{\left({\underset{\‾}{x}}_{ks,kf}{{\underset{\‾}{x}}_{ks,kf}}^H\right)}^{*}\right\}\end{array}---\left(4\right)$

In formula, E represent ask expectation, * to represent asking conjugation, H to represent seeks conjugate transpose,Represent Kronecker product.

Step 3, to fourth order cumulantCarry out EVD feature decomposition, and by obtained (MF)²Individual eigenvalue is arranged as from big to small:Wherein, the frequency number that F chooses when being frequency domain smoothing; M is the number of sensors in acoustic sensor array;

EVD feature decomposition is prior art, and its expression formula is as follows:

$\hat{C}={\Σ}_{k=1}^{{\left(MF\right)}^2}{\mathrm{\λ}}_k{\mathrm{\μ}}_k{{\mathrm{\μ}}_k}^{*}.---\left(5\right)$

Wherein, the frequency number that F chooses when being frequency domain smoothing; M is the number of sensors in acoustic sensor array; ��_kFor rightCarry out the kth eigenvalue that feature decomposition obtains; ��_kFor rightThe kth eigenvalue �� obtained after carrying out feature decomposition_kCorresponding unit character vector; Eigenvalue is arranged as from big to small:

Step 4, rightObtain (MF) after carrying out feature decomposition²Individual eigenvalueWith wherein P²Individual bigger eigenvalueCharacteristic of correspondence vector structure signal subspace, with remaining (MF)²-P²Individual less eigenvalueCharacteristic of correspondence vector structure noise subspace;

Then fourth order cumulantCan be expressed as again:

$\hat{C}={\Σ}_{k=1}^{{\left(MF\right)}^2}{\mathrm{\λ}}_k{\mathrm{\μ}}_k{{\mathrm{\μ}}_k}^{*}={\Σ}_{k=1}^{P^2}{\mathrm{\λ}}_k{\mathrm{\μ}}_k{{\mathrm{\μ}}_k}^{*}+{\Σ}_{k=P^2+1}^{{\left(MF\right)}^2}{\mathrm{\λ}}_k{\mathrm{\μ}}_k{{\mathrm{\μ}}_k}^{*}.---\left(6\right)$

So P²The signal subspace projection that individual bigger eigenvalue characteristic of correspondence vector characterizes is represented by:

${\hat{C}}_s={\Σ}_{k=1}^{P^2}{\mathrm{\μ}}_k{{\mathrm{\μ}}_k}^{*}.---\left(7\right)$

(MF)²-P²The noise subspace projection that individual less eigenvalue characteristic of correspondence vector characterizes is represented by:

${\hat{C}}_n={\Σ}_{k=P^2+1}^{{\left(MF\right)}^2}{\mathrm{\μ}}_k{{\mathrm{\μ}}_k}^{*}.---\left(8\right)$

Step 5, for sound ray path configuration copy vector a (��, T):

Including the information of parameter to be estimated in copy vector, rational copy vector not only can improve location separation accuracy, and if add launch signal spectrum information, it is also possible to improve the algorithm robustness to noise. The present invention utilizes the amplitude launching signal and phase information when constructing copy vector, utilizes signal wide-band-message at frequency domain simultaneously, improves separation accuracy and the noise immunity of algorithm. The specific configuration form of copy vector a (��, T) is as follows:

$a(\mathrm{\θ},T)=b(\mathrm{\θ},T)\⊗b{(\mathrm{\θ},T)}^{*}---\left(9\right)$

$b(\mathrm{\θ},T)=\left[\begin{array}{c}e\left(v_1\right)e^{-2{\mathrm{i\πv}}_{1}T}d(v_1,\mathrm{\θ})\\ \mathrm{...}\\ e\left(v_F\right)e^{-2{\mathrm{i\πv}}_{F}T}d(v_F,\mathrm{\θ})\end{array}\right]---\left(10\right)$

$d(v_i,\mathrm{\θ})=\[1,e^{-2{\mathrm{i\πv}}_i{\mathrm{\τ}}_{1,2}\left(\mathrm{\θ}\right)},...,e^{-2{\mathrm{i\πv}}_i{\mathrm{\τ}}_{1,M-1}\left(\mathrm{\θ}\right)}\]---\left(11\right)$

Wherein, �� represents the angle of arrival in sound ray path, and T represents the time of advent in sound ray path, e (v_i) it is that acoustical signal is at frequency v_iThe amplitude at place, i=1,2 ..., F,Represent Kronecker product, ��_1,j(��) represent that sound ray path arrives jth sensor relative to arriving the time delay of the 1st sensor as reference sensor, j=2 ..., M-1;

Step 6, the search copy vector orthogonal with described noise subspace, these copy vectors are the sound ray path being finally recovered out, and its �� and T is the corresponding arrival direction in sound ray path and the time of advent;

Existing document shows, signal subspace is orthogonal with noise subspace, accordingly, it is possible to structure is suitable for the estimator of higher order statistical. The form of the present invention constructed estimator P (��, T) is as follows:

$P(\mathrm{\θ},T)=\frac{1}{a{(\mathrm{\θ},T)}^{*}{\hat{C}}_{n}a(\mathrm{\θ},T)}---\left(12\right)$

When function P (��, T) obtains maximum, copy vector and noise subspace are mutually perpendicular to, and what copy vector parameters �� now represented is exactly the direction in sound ray path, and T represents that sound ray path arrives the time of sensor.

After completing the separation of sound ray path, it is possible to carry out follow-up process, for instance, the orientation of signal source can be determined in the two dimensional surface that angle and time are constituted by parameter �� and T; Or, according to parameter T, utilize sound ray propagation time chromatography method to be finally inversed by ocean environment parameter.

In order to verify the effect of the present invention, it is carried out simulation comparison experiment with existing smooth active broadband multi-signals separation algorithm (being called for short smoothing-MUSICAL). The relevant information of ray path that experiment uses is: five ray paths are respectively as follows: 0,1 ,-1,1.6 ,-1.6 the time delay between each sensor; Article five, ray path arrives the time respectively 10,12,12.5,14,14 of sensor; Choosing first sensor is reference sensor; Sensor sample data length is 129; Spatial domain and the smooth order of frequency domain that two kinds of methods adopt are 1.

Smoothing-MUSICAL algorithm is based on the second-order cumulant of sensor sample data, and assumes that noise is Gaussian noise; The inventive method (abbreviation 4-Smoothing-MUSICAL) is based on the fourth order cumulant of sensor sample data, assumes that noise is chromatic noise, more meets the marine environment of reality simultaneously. Experimental subject for same: be respectively as follows: 0,1 ,-1,1.6 ,-1.6 time delay between sensor, signal arrive sensor time respectively 10,12,12.5,14,14 five ray paths:

1) not plus noise, number of sensors is 4, and choosing first sensor is reference sensor; The result of Smoothing-MUSICAL algorithm and 4-Smoothing-MUSICAL algorithm shown in two width figure of the first row, is the result adopting Smoothing-MUSICAL algorithm to obtain on the left of wherein respectively in Fig. 4, and right side is the result of the inventive method.

2) not plus noise, number of sensors is 6, and choosing first sensor is reference sensor; The result of Smoothing-MUSICAL algorithm and 4-Smoothing-MUSICAL algorithm shown in two width figure of the second row, is the result adopting Smoothing-MUSICAL algorithm to obtain on the left of wherein respectively in Fig. 4, and right side is the result of the inventive method.

3) not plus noise, number of sensors is 7, and choosing first sensor is reference sensor; The result of Smoothing-MUSICAL algorithm and 4-Smoothing-MUSICAL algorithm shown in the two width figure of the third line, is the result adopting Smoothing-MUSICAL algorithm to obtain on the left of wherein respectively in Fig. 4, and right side is the result of the inventive method.

4) variegating noise, signal to noise ratio is 20dB, and number of sensors is 4, and choosing first sensor is reference sensor; The result of Smoothing-MUSICAL algorithm and 4-Smoothing-MUSICAL algorithm shown in two width figure of fourth line, is the result adopting Smoothing-MUSICAL algorithm to obtain on the left of wherein respectively in Fig. 4, and right side is the result of the inventive method.

5) variegating noise, signal to noise ratio is-5dB, and number of sensors is 4, and choosing first sensor is reference sensor; The result of Smoothing-MUSICAL algorithm and 4-Smoothing-MUSICAL algorithm shown in two width figure of fifth line, is the result adopting Smoothing-MUSICAL algorithm to obtain on the left of wherein respectively in Fig. 4, and right side is the result of the inventive method.

6) variegating noise, signal to noise ratio is-10dB, and number of sensors is 4, and choosing first sensor is reference sensor; The result of Smoothing-MUSICAL algorithm and 4-Smoothing-MUSICAL algorithm shown in two width figure of the 6th row, is the result adopting Smoothing-MUSICAL algorithm to obtain on the left of wherein respectively in Fig. 4, and right side is the result of the inventive method.

Visible the inventive method of contrast and experiment according to Fig. 4 has the advantage that compared to existing Smoothing-MUSICAL algorithm

1) ray path that interval is less can be properly separated out;

2) number of sensors situation less than ray path can be isolated

3) impact of non-Gaussian noise is suppressed;

4) high separation accuracy;

5) less artifact is produced.

Claims (4)

1. based on the multipath propagation acoustic signal of Higher Order Cumulants, for isolating every sound ray path from the multipath propagation acoustical signal received by acoustic sensor array, it is characterised in that comprise the following steps:

Step 1, the received acoustical signal of sensor is done the laggard line space territory of Fourier transformation-frequency domain smoothing process, it is thus achieved that K=(2K_s+1)(2K_f+ 1) composition matrix estimated by individual arrowbandx _{Ks, kf}, K_s��K_fRespectively spatial domain smooth order, frequency domain smoothing order; And K is more than sound source number P;

Step 2, utilize following formula calculate matrixx _{Ks, kf}Fourth order cumulant

$\begin{array}{c}\hat{C}=E\left\{\left({\underset{\‾}{x}}_{ks,kf}\⊗{{\underset{\‾}{x}}_{ks,kf}}^{*}\right){\left({\underset{\‾}{x}}_{ks,kf}\⊗{{\underset{\‾}{x}}_{ks,kf}}^{*}\right)}^H\right\}-E\left\{{\underset{\‾}{x}}_{ks,kf}\⊗{{\underset{\‾}{x}}_{ks,kf}}^{*}\right\}E\left\{\right({\underset{\‾}{x}}_{ks,kf}\⊗\\ {{\underset{\‾}{x}}_{ks,kf}}^{*}{)}^H\}-E\left\{{\underset{\‾}{x}}_{ks,kf}{{\underset{\‾}{x}}_{ks,kf}}^H\right\}\⊗E\left\{{\left({\underset{\‾}{x}}_{ks,kf}{{\underset{\‾}{x}}_{ks,kf}}^H\right)}^{*}\right\}\end{array}$

Wherein, E represent ask expectation, * to represent asking conjugation, H to represent seeks conjugate transpose,Represent Kronecker product;

Step 3, to fourth order cumulantCarry out EVD feature decomposition, and by obtained (MF)²Individual eigenvalue is arranged as from big to small:Wherein, the frequency number that F chooses when being frequency domain smoothing; M is the number of sensors in acoustic sensor array;

Step 4, rightObtain (MF) after carrying out feature decomposition²Individual eigenvalueWith wherein P²Individual bigger eigenvalue ��₁������Characteristic of correspondence vector structure signal subspace, with remaining (MF)²-P²Individual less eigenvalueCharacteristic of correspondence vector structure noise subspace;

Step 5, for copy vector a (��, T) below sound ray path configuration:

$a(\mathrm{\θ},T)=b(\mathrm{\θ},T)\⊗b{(\mathrm{\θ},T)}^{*}$

$b(\mathrm{\θ},T)=\left[\begin{array}{c}e\left(v_1\right)e^{-2{\mathrm{i\πv}}_{1}T}d(v_1,\mathrm{\θ})\\ ...\\ e\left(v_F\right)e^{-2{\mathrm{i\πv}}_{F}T}d(v_F,\mathrm{\θ})\end{array}\right]$

$d(v_i,\mathrm{\θ})=\[1,e^{-2{\mathrm{i\πv}}_i{\mathrm{\τ}}_{1,2}\left(\mathrm{\θ}\right)},...,e^{-2{\mathrm{i\πv}}_i{\mathrm{\τ}}_{1,M-1}\left(\mathrm{\θ}\right)}\]$

Wherein, �� represents the angle of arrival in sound ray path, and T represents the time of advent in sound ray path, e (v_i) it is that acoustical signal is at frequency v_iThe amplitude at place, i=1,2 ..., F,Represent Kronecker product, ��_1,j(��) represent that sound ray path arrives jth sensor relative to arriving the time delay of the 1st sensor as reference sensor, j=2 ..., M-1;

Step 6, the search copy vector orthogonal with described noise subspace, these copy vectors are the sound ray path being finally recovered out, and its �� and T is the corresponding arrival direction in sound ray path and the time of advent.

2. multipath propagation acoustic signal as claimed in claim 1, it is characterised in that described frequency domain smoothing uses frequency domain sub-band averaging method.

3. the Ocean Acoustic Tomography method based on sound ray propagation time chromatography, sound change of spread speed in ocean is utilized to carry out inverting ocean environment parameter, it is characterized in that, first with multipath propagation acoustic signal based on Higher Order Cumulants described in claim 1 or 2, from the multipath propagation acoustical signal received by acoustic sensor array, isolate each sound ray path; Then it is finally inversed by ocean environment parameter the time of advent according to sound ray path.

4. a sound localization method, it is characterised in that first with multipath propagation acoustic signal based on Higher Order Cumulants described in claim 1 or 2, isolate each sound ray path from the multipath propagation acoustical signal received by acoustic sensor array; Then according to the angle of arrival in sound ray path and the orientation determining sound source the time of advent in the two dimensional surface spatially and temporally constituted.