CN105929386A - 一种基于高阶累积量的波达估计方法 - Google Patents

Abstract

本发明公开了一种基于高阶累积量的波达估计方法，属于信号处理技术领域。该方法利用等间距直线传感器阵列所接收到的观测信号，估计出信号源的波达方向及波达时间；包括以下步骤：步骤1、对观测信号做傅里叶变换后进行空域‑频域平滑处理；步骤2、构造出空域‑频域平滑处理后信号的四阶累积量矩阵；步骤3、利用迭代地部分SVD方法，根据四阶累积量矩阵构建观测信号的信号子空间和噪声子空间；步骤4、根据观测信号的信号子空间和噪声子空间之间的正交性，估计出信号源的波达方向及波达时间。本发明还公开了一种基于声线传播时间层析的海洋声层析方法及一种定位方法。本发明可大幅降低算法的计算复杂度，提高算法实时性并降低硬件资源消耗。

Description

一种基于高阶累积量的波达估计方法

技术领域

本发明涉及信号处理技术领域，尤其涉及一种基于高阶累积量的波达估计方法。

背景技术

波达方向(DOA)估计指的是要确定同时处在空间某一区域内多个感兴趣信号的空间位置，即各个信号到达阵列参考阵元的方向角。波达方向技术是阵列信号处理中的重要研究课题。波达方向估计技术的关键在于利用处于空间不同位置的天线阵列，接收多个不同方向的信号源发出的信号，运用现代处理方法快速且高精度地估计出信号源的方向。波达方向估计技术目前广泛应用于雷达、声纳、地震等领域。

八十年代以后，学术界提出一类基于矩阵特征值分解的谱估计方法，以Schmidt等人提出的MUSIC(Multiple Signal Classification)算法为代表，它是基于信号子空间和噪声子空间的正交性，具有很好的角度分辨能力。但是当信号相关或者相干时，样本的协方差矩阵会不满秩，MUSIC算法则不能给出正确的分离结果。为了解决这个问题，学者们给出了许多方法。特别地，Jiang等人在“Raypath separation with high resolution processing”一文中，提出在点到阵列结构中，结合空域-频域平滑方法与主动宽带多信号分离算法，给出平滑的主动宽带多信号分离算法smoothing-MUSICAL(smoothing Multiple Signal ClassificationActive Large-band)。该算法仍然存在着一些问题：(1)分辨能力并不能满足实际应用的需求；(2)该算法基于二阶统计量，需要假设信号为高斯信号；(3)传感器阵列的阵元数目必须大于待分离信号源的数目。

为了解决这些问题，一篇中国发明专利申请(发明名称为“基于高阶累积量的多路径传播声信号分离方法”，申请号为201610006173.3，申请日为2016-1-5)提出了采用基于高阶累积量的阵列信号处理方法。该方法的具体步骤如下：

步骤1、对传感器所接收声信号做傅里叶变换后进行空域-频域平滑处理，获得K＝(2K_s+1)(2K_f+1)个窄带估计组成矩阵x _ks，kf，K_s、K_f分别为空域平滑阶次、频域平滑阶次；且K大于声源个数P；

步骤2、利用下式计算矩阵x _ks，kf的四阶累积量

$\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\{{\underset{\‾}{x}}_{ks,kf}{{\underset{\‾}{x}}_{ks,kf}}^H\}\⊗H\{{\left({\underset{\‾}{x}}_{ks,kf}{{\underset{\‾}{x}}_{ks,kf}}^H\right)}^{*}\}\end{array}$

其中，E表示求期望，*表示求共轭，H表示求共轭转置，表示克罗内克积；

步骤3、对四阶累积量进行EVD特征分解，并将所得到的(MF)²个特征值从大到小排列为：其中，F为频域平滑时选取的频率数；M为声传感器阵列中的传感器数量；

步骤4、对进行特征分解后得到(MF)²个特征值以

其中P²个较大的特征值对应的特征向量构造信号子空间，以剩下的(MF)²-P²个较小特征值对应的特征向量构造噪声子空间；

步骤5、为声线路径构造以下的副本矢量a(θ，T)：

$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)},\mathrm{...},e^{-2{\mathrm{i\πv}}_i{\mathrm{\τ}}_{1,M-1}\left(\mathrm{\θ}\right)}\]$

其中，θ表示声线路径的到达角度，T表示声线路径的到达时间，e(v_i)是声信号在频率v_i处的幅值，i＝1,2,…,F，表示克罗内克积，τ_1,j(θ)表示声线路径到达第j个传感器相对于到达作为参考传感器的第1个传感器的时间延迟，j＝2,…,M-1；

步骤6、搜索与所述噪声子空间正交的副本矢量，这些副本矢量即为最终分离出的声线路径，其θ和T即为声线路径相应的到达方向和到达时间。

由于高阶累积量对任意的高斯噪声有自然盲性，基于高阶累积量的算法使原有的DOA估计所适应的观测噪声扩展到高斯空间有色噪声或者对称的非高斯空间有色噪声和白噪声。

然而，上述基于高阶累积量的算法，由于引入了样本的高阶累积量矩阵，这将导致矩阵规模极其庞大，随之而来的问题是对更大的存储空间和内存的需求，同时基于特征值分解的这样一个计算过程的时间消耗也会指数级增长。

发明内容

本发明所要解决的技术问题在于克服现有技术不足，提供一种基于高阶累积量的波达估计方法，一方面可获得更高的信号估计精度并能适用于传感器数目小于射线路径的情况，同时还可抑制高斯噪声和非高斯噪声的影响；另一方面可大幅降低算法的计算复杂度，提高算法实时性并降低硬件资源消耗。

本发明具体采用以下技术方案解决上述技术问题：

一种基于高阶累积量的波达估计方法，利用等间距直线传感器阵列所接收到的观测信号，估计出信号源的波达方向及波达时间；包括以下步骤：

步骤1、对所述观测信号做傅里叶变换后进行空域-频域平滑处理，得到K＝(2K_s+1)(2K_f+1)个窄带估计作为列向量所组成的矩阵X，K_s、K_f分别为空域平滑阶次、频域平滑阶次；且K大于信号源个数N；

步骤2、根据下式构造出空域-频域平滑处理后信号的四阶累积量矩阵C：

$\begin{array}{c}C=E\left\{(X\⊗X^{*}){(X\⊗X^{*})}^H\right\}-E\{X\⊗X^{*}\}E\left\{{(X\⊗X^{*})}^H\right\}\\ -E\left\{{\mathrm{XX}}^H\right\}\⊗E\left\{{\mathrm{XX}}^H\right\}\end{array}$

X表示空域-频域平滑处理之后得到的矩阵，E表示期望，上标*表示共轭，上标H表示共轭转置，表示克罗内克积；

步骤3、根据所述四阶累积量矩阵构建观测信号的信号子空间和噪声子空间；

步骤4、根据观测信号的信号子空间和噪声子空间之间的正交性，估计出信号源的波达方向及波达时间；

所述步骤3具体包括以下子步骤：

步骤301、对两个规模均为(MF)²×p的正交矩阵U、V进行随机初始化，记为U⁰、V⁰，矩阵右上角的数字标号表示迭代的次数，M为所述等间距直线传感器阵列的阵元个数，F为空域-频域平滑处理过程中选取的频率通道数，p＝K²且p＜＜(MF)²；

步骤302、对所述两个正交矩阵U、V进行迭代处理，并在每一次迭代后判断预设条件是否满足，如是，则停止迭代并转至下一步骤，否则，继续下一次迭代：对于第k次迭代：

首先对矩阵做部分SVD如下：

$U^{k^H}C=V^{k+1}{S}_{temp1}{V_{temp1}}^H$

再对矩阵CV^k做部分SVD如下：

CV^k＝U^k+1S_temp2V_temp2 ^H

其中：

矩阵右上角的数字标号k表示第k次迭代的初值，每一次的迭代更新该值；

U^k，V^k表示第k次迭代的初值；

U^k+1，V^k+1表示第k次迭代的更新之后的结果，亦作为下一次迭代的初值；

U和V的规模均为：(MF)²×p；

下标temp1和temp2表示该矩阵为临时变量，不需要保存；

所述预设条件具体如下：

范数1：

范数2：

其中，表示计算矩阵A最大的奇异值，||A||_F表示矩阵A的Frobenius范数，ε为预设阈值；

步骤303、截取当前矩阵U的前N²列，用表示，即为观测信号的信号子空间，。

优选地，所述根据观测信号的信号子空间和噪声子空间之间的正交性，估计出信号源的波达方向及波达时间，具体如下：

求取使得如下估计函数取最大值时的θ、T，即为信号源的波达方向及波达时间：

$P(\mathrm{\θ},T)=a{(\mathrm{\θ},T)}^{*}{U}_{N^2}{U_{N^2}}^{*}a(\mathrm{\θ},T)$

其中：

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

$b(\mathrm{\θ},T)=\left[\begin{array}{c}e\left({\mathrm{\υ}}_1\right)e^{-2{\mathrm{j\π\υ}}_{1}T}d({\mathrm{\υ}}_1,\mathrm{\θ})\\ .\\ .\\ .\\ e\left({\mathrm{\υ}}_F\right)e^{-2{\mathrm{j\π\υ}}_{F}T}d({\mathrm{\υ}}_F,\mathrm{\θ})\end{array}\right]$

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

θ表示信号路径的波达角度；T表示信号路径的波达时间；e(υ_i)表示信号在频率υ_i处的幅值，i＝1,2,...,F；τ_1,j(θ)表示信号路径到达第j个传感器相对于到达参考传感器的时间延迟，j＝2,…,M-1。

优选地，阈值ε的取值范围为1e-3～1e-6。

优选地，所述空域-频域平滑中的频域平滑使用频域子带平均方法。

上述波达估计方法可广泛用于雷达、声纳、地震等领域，以下为两个具体应用方案：

一种基于声线传播时间层析的海洋声层析方法，利用声音在海洋中传播速度的变化来反演海洋环境参数，首先利用以上任一技术方案所述基于高阶累积量的波达估计方法对从声传感器阵列所接收到的多路径传播声信号进行波达估计，从而分离出每一条声线路径；然后根据声线路径的到达时间反演出海洋环境参数。

一种定位方法，首先利用以上任一技术方案所述基于高阶累积量的波达估计方法进行波达方向估计，然后利用估计出的波达方向确定信号源的位置。

相比现有技术，本发明具有以下有益效果：

首先，和基于二阶统计量的技术相比，本发明具有更高的信号路线分离精度，可以分离出间隔较小的路径，且本发明可以在传感器阵列阵元数目小于待分离信号源数目的情况下正确分离出信号的路径。同时本发明中假设的加性噪声是彩色噪声，而传统的技术都是假设噪声为高斯噪声，而彩色噪声会更符合实际情况，本发明可以更好地抑制各类环境噪声，包括高斯噪声和非高斯噪声。

其次，和基于特征值分解的高阶累积量DOA算法相比，本发明利用迭代地部分SVD的方法，设定收敛标准，可以极大地加快算法的计算速度，减少耗时。

附图说明

图1为空域平滑方法的一个原理示例；

图2为频域平滑方法的一个原理示例；

图3为本发明方法的流程示意图；

图4为本发明方法与Smoothing-MUSICAL算法以及四阶累积量基于特征值分解方法的对比实验结果。

具体实施方式

下面结合附图对本发明的技术方案进行详细说明：

本发明的发明思路是针对基于高阶累积量的阵列信号处理方法(CN201610006173.3)所存在的计算复杂度过高的问题，在根据四阶累积量矩阵构建观测信号的信号子空间和噪声子空间过程中，利用迭代地部分奇异值分解的方法可以大幅减少计算复杂度，从而大幅节约算法的耗时。

为了便于公众理解，下面以一个具体实施例来对本发明技术方案进行详细说明。

首先建立信号模型：

各待测信号源是相关的或相干的，设其数目为N。传感器阵列中阵元数目为M，并不需要限制M＞N，传感器阵列是等间距直线阵，各阵元的特性相同，间隔为d，且阵元间隔不大于最高频率信号的1/2波长。信号源的传播过程中有彩色加性噪声。

第m个传感器上接收到的信号可以表示为：

$x_m\left(t\right)=\underset{n=1}{\overset{N}{\Σ}}{a}_{n}e(t-{\mathrm{\τ}}_{m,n})+{\mathrm{noise}}_m\left(t\right)---\left(1\right)$

其中a_n为第n个射线路径对第m传感器的影响，e(t)为信号源发射的信号，τ_m,n为第n个射线路径到达第m个传感器的时间，noise_m(t)表示第m个传感器上的加性彩色噪声。

到达时间可以表示为：

τ_m,n＝T_n+t_m(θ_n) (2)

其中T_n表示第n个射线路径到达参考传感器的时间，t_m(θ_n)则表示第n个射线路径到达第m个传感器相比参考传感器的时间延迟，角度θ_n则表示第n个射线路径到达传感器的方向。

式(1)的频域表示为：

$x_m\left(\mathrm{\υ}\right)=\underset{n=1}{\overset{N}{\Σ}}{a}_{n}e\left(\mathrm{\υ}\right)\exp (-j2{\mathrm{\π\τ}}_{m,n})+{\mathrm{noise}}_m\left(\mathrm{\υ}\right)---\left(3\right)$

在以上信号模型的基础上，本发明方法具体包括以下步骤：

步骤1.对传感器所接收信号进行空域-频域的平滑处理，获得K＝(2K_s+1)(2K_f+1)个窄带估计作为列向量所组成的矩阵X，其中K_s和K_f分别为空域平滑阶次和频域平滑阶次，且K大于信号源个数N。

空域-频域平滑为现有技术，下面对其进行简要说明。

空域平滑通过对空域子阵列的平均来实现。如图1所示为空域平滑的原理图。M个传感器被分成若干尺寸相同，部分重叠的子阵列，假设子阵列是线性一致的，则子阵列传感器上信号的强度不会发生剧烈变化。当子阵列数目大于或等于射线路径数时空间谱矩阵为非单秩矩阵。假设空间平滑阶次为K_s，则每个子阵列的尺寸为M-2K_s，子阵列的数目为2K_s+1。

频域平滑原理图如图2所示，根据操作对象是时域数据还是频域数据可分为如下两种：(1)加权相关矩阵；(2)频域子带平均。本发明优选采用频域子带平均法。假设频域的平滑阶次为K_f，则平滑后可以得到2K_f+1个尺寸为M-2K_f的子带。

对空域或频域的过度平滑均会导致互谱矩阵产生严重的误差，为了减小估计误差，需要将空域和频域的平滑相结合。

传感器阵列的阵元个数用M表示，选取的频率通道数用F表示，则平滑后以K个窄带估计作为列向量所组成的矩阵X的规模为：(MF)×K。

步骤2.构造平滑后的矩阵X的四阶累积量矩阵如下：

$\begin{array}{c}C=E\left\{(X\⊗X^{*}){(X\⊗X^{*})}^H\right\}-E\{X\⊗X^{*}\}E\left\{{(X\⊗X^{*})}^H\right\}\\ -E\left\{{\mathrm{XX}}^H\right\}\⊗E\left\{{\mathrm{XX}}^H\right\}\end{array}---\left(4\right)$

其中：X表示平滑处理之后的观测矩阵；C表示四阶累积量矩阵；E表示期望；上标*表示共轭；上标H表示共轭转置；表示克罗内克积。构造出的四阶累积量的规模为：(MF)²×(MF)²。

步骤3.该步骤是一个迭代的过程：

(1)迭代前先初始化，生成两个随机的正交矩阵：U⁰、V⁰，作为U、V的初始值，这意味着这两个正交矩阵的规模均为：(MF)²×p。

其中：

p＝K²且p＜＜(MF)²。

I表示单位矩阵。

矩阵右上角的数字标号表示迭代的次数，下面沿用该表示方式。

(2)迭代的更新(以第k次迭代为例，迭代从第0次开始)：

首先对矩阵做部分SVD如下：

$U^{k^H}C=V^{k+1}{S}_{temp1}{V_{temp1}}^H---\left(5\right)$

再对矩阵CV^k做部分SVD如下：

CV^k＝U^k+1S_temp2V_temp2 ^H (6)

其中：

矩阵右上角的数字标号k表示第k次迭代的初值，每一次的迭代更新该值；U^k，V^k表示第k次迭代的初值；U^k+1，V^k+1表示第k次迭代的更新之后的结果，亦作为下一次迭代的初值；U和V的规模均为：(MF)²×p。下标temp1和temp2表示该矩阵我们并不关注。

(3)终止条件：设定作为收敛判定标准的阈值ε的取值，经大量实验发现，ε的量级选取在1e-3至1e-6左右时精确度以及计算复杂度都会比较好，由于越接近最优值收敛的速度越慢，如果标准设定在1e-7及更小的数量级，迭代的次数会大幅增加，相比直接计算四阶累积量矩阵的特征值其计算复杂度的优势不会那么明显，同样的，如果收敛标准的值过大，虽然会极快地地达到收敛标准停止迭代，但是精确度也会降低，因此本发明优选在该范围内取值，例如可取ε＝1e-5。当然本领域技术人员也可根据需要确定其取值，如需要突出精确度，则可取较小的值，反之，如需要快速计算，则可取较大的值。

计算两个范数：

范数1：

范数2：

${||\[\left(V^k{V}^{k^H}\right)-I\]{\mathrm{CV}}^k{\left(U^{k^H}{\mathrm{CV}}^k\right)}^H||}_F+{||\[\left(U^k{U}^{k^H}\right)-I\]C^H{U}^k\left(U^{k^H}{\mathrm{CV}}^k\right)||}_F---\left(8\right)$

其中；计算的是矩阵A最大的奇异值；||A||_F是Frobenius范数，计算的是矩阵A协方差矩阵的迹的二次方根。

每一次迭代我们都作一次判定，若：则停止迭代。

最终得到的(MF)²×p矩阵U和V则分别是对应于四阶累积量矩阵的前p个最大的奇异值的左奇异向量子空间和右奇异向量子空间。

信号源个数用N来表示，截取矩阵U的前N²列，用表示，即为我们需要的信号子空间。

现有方法是直接进行累积量矩阵C的特征值分解计算，从而来获得信号子空间和噪声子空间，然而应用了高阶累积量之后，C的规模会变得很大，四阶累积量的情况下其规模为(MF)²×(MF)²。直接对其做特征值分解的计算量耗时相当严重，而在本发明方法中，并不直接对该规模的矩阵做特征值分解，而是在有限次的迭代过程中做规模为(MF)²×p的矩阵的部分SVD，由于p＜＜(MF)²，因此计算复杂度可以极大地降低，从而节约计算时的时间消耗。

步骤4.求取如下估计函数的值：

$P(\mathrm{\θ},T)=a{(\mathrm{\θ},T)}^{*}{U}_{N^2}{U_{N^2}}^{*}a(\mathrm{\θ},T)---\left(9\right)$

其中：

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

$b(\mathrm{\θ},T)=\left[\begin{array}{c}e\left({\mathrm{\υ}}_1\right)e^{-2{\mathrm{j\π\υ}}_{1}T}d({\mathrm{\υ}}_1,\mathrm{\θ})\\ .\\ .\\ .\\ e\left({\mathrm{\υ}}_F\right)e^{-2{\mathrm{j\π\υ}}_{F}T}d({\mathrm{\υ}}_F,\mathrm{\θ})\end{array}\right]---\left(11\right)$

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

θ表示信号路径的波达角度；T表示信号路径的波达时间；e(υ_i)表示信号在频率υ_i处的幅值，i＝1,2,...,F；τ_1,j(θ)表示信号路径到达第j个传感器相对于到达作为参考传感器的时间延迟，j＝2,…,M-1。

估计函数的最大值所对应的θ和T即信号源的达波方向和达波时间。整个处理过程的流程参见图3。

为了验证本发明的效果，将其与现有的smoothing-MUSICAL算法进行仿真对比实验，同时与基于特征值分解的四阶累积量算法进行对比。实验使用的射线路径的相关信息为：五条射线路径在各个传感器间的延迟时间分别为：0、1、-1、1.6、-1.6；五条射线路径到达传感器的时间分别为10、12、12.5、14、14；选取第一个传感器为参考传感器；传感器采样数据长度为129；两种方法所采用的空域和频域的平滑阶次均为1。且该实验中收敛标准设定为1e-3，可以更突出其在计算复杂度方面的提升。算法运行的操作系统是Red Hat Enterprise LinuxServer release 5.6(Tikanga)，内存24GB，处理器是Intel X5650@2.67GHz。

对比实验的结果如图4所示，图中每行的3个结果分别是smoothing-MUSICAL(左侧的结果)和四阶累积量基于特征值分解(中间的结果)和本发明方法(右侧的结果)的方法使用同样的观测数据得出的结果。同样的实验对象是5个信号源：传感器间的延迟时间分别为：0、1、-1、1.6、-1.6，信号到达传感器的时间分别为10、12、12.5、14、14的五条射线路径，图中从最上面一行至最后一行的实验条件依次如下：

(1)不加噪声，传感器数目为4，选取第一个传感器为参考传感器；

(2)加彩色噪声，信噪比为20dB，传感器数目为4，选取第一个传感器为参考传感器；

(3)加彩色噪声，信噪比为-5dB，传感器数目为4，选取第一个传感器为参考传感器；

(4)加彩色噪声，信噪比为-10dB，传感器数目为4，选取第一个传感器为参考传感器；

(5)不加噪声，传感器数目为6，选取第一个传感器为参考传感器；

(6)加彩色噪声，信噪比为20dB，传感器数目为6，选取第一个传感器为参考传感器；

(7)加彩色噪声，信噪比为-5dB，传感器数目为6，选取第一个传感器为参考传感器；

(8)加彩色噪声，信噪比为-10dB，传感器数目为6，选取第一个传感器为参考传感器；

表1为该实验中四阶累积量DOA方法分别基于特征值分解(即CN201610006173.3中的方法)和基于部分SVD迭代算法(即本发明方法)获取信号子空间过程的时间消耗。

表1.四阶累积量DOA算法分别基于特征值分解和部分SVD迭代的时间消耗对比

根据以上对比实验可知，与基于二阶统计量的技术相比，本发明具有更高的信号路线分离精度，可以分离出间隔较小的路径，且本发明可以在传感器阵列阵元数目小于待分离信号源数目的情况下正确分离出信号的路径。同时本发明可以抑制非高斯噪声的影响。

和基于特征值分解的高阶累积量DOA算法相比，本发明基于部分SVD迭代的方法，可以极大地加快算法的计算速度，减少耗时。

本发明波达估计方法可可广泛用于雷达、声纳、地震等领域，以下为两个具体应用方案：

一种基于声线传播时间层析的海洋声层析方法，利用声音在海洋中传播速度的变化来反演海洋环境参数，首先利用以上任一技术方案所述基于高阶累积量的波达估计方法对从声传感器阵列所接收到的多路径传播声信号进行波达估计，从而分离出每一条声线路径；然后根据声线路径的到达时间反演出海洋环境参数。

一种定位方法，首先利用以上任一技术方案所述基于高阶累积量的波达估计方法进行波达方向估计，然后利用估计出的波达方向确定信号源的位置。

Claims (6)

1.一种基于高阶累积量的波达估计方法，利用等间距直线传感器阵列所接收到的观测信号，估计出信号源的波达方向及波达时间；包括以下步骤：

步骤1、对所述观测信号做傅里叶变换后进行空域-频域平滑处理，得到K＝(2K_s+1)(2K_f+1)个窄带估计作为列向量所组成的矩阵X，K_s、K_f分别为空域平滑阶次、频域平滑阶次；且K大于信号源个数N；

步骤2、根据下式构造出空域-频域平滑处理后信号的四阶累积量矩阵C：

$\begin{array}{c}C=E\left\{(X\⊗X^{*}){(X\⊗X^{*})}^H\right\}-E\{X\⊗X^{*}\}E\left\{{(X\⊗X^{*})}^H\right\}\\ -E\left\{{\mathrm{XX}}^H\right\}\⊗E\left\{{\mathrm{XX}}^H\right\}\end{array}$

X表示空域-频域平滑处理之后得到的矩阵，E表示期望，上标*表示共轭，上标H表示共轭转置，表示克罗内克积；

步骤3、根据所述四阶累积量矩阵构建观测信号的信号子空间和噪声子空间；

步骤4、根据观测信号的信号子空间和噪声子空间之间的正交性，估计出信号源的波达方向及波达时间；

其特征在于，所述步骤3具体包括以下子步骤：

步骤301、对两个规模均为(MF)²×p的正交矩阵U、V进行随机初始化，记为U⁰、V⁰，矩阵右上角的数字标号表示迭代的次数，M为所述等间距直线传感器阵列的阵元个数，F为空域-频域平滑处理过程中选取的频率通道数，p＝K²且p＜＜(MF)²；

步骤302、对所述两个正交矩阵U、V进行迭代处理，并在每一次迭代后判断预设条件是否满足，如是，则停止迭代并转至下一步骤，否则，继续下一次迭代：对于第k次迭代：

首先对矩阵C做部分SVD如下：

$U^{k^H}C=V^{k+1}{S}_{temp1}{V_{temp1}}^H$

再对矩阵CV^k做部分SVD如下：

${\mathrm{CV}}^k=U^{k+1}{S}_{temp2}{V_{temp2}}^H$

其中：

矩阵右上角的数字标号k表示第k次迭代的初值，每一次的迭代更新该值；

U^k，V^k表示第k次迭代的初值；

U^k+1，V^k+1表示第k次迭代的更新之后的结果，亦作为下一次迭代的初值；

U和V的规模均为：(MF)²×p；

下标temp1和temp2表示该矩阵为临时变量，不需要保存；

所述预设条件具体如下：

范数1：

范数2：

其中，表示计算矩阵A最大的奇异值，表示矩阵A的Frobenius范数，ε为预设阈值；

步骤303、截取当前矩阵U的前N²列，用表示，即为观测信号的信号子空间。

2.如权利要求1所述波达估计方法，其特征在于，所述根据观测信号的信号子空间和噪声子空间之间的正交性，估计出信号源的波达方向及波达时间，具体如下：

求取使得如下估计函数取最大值时的θ、T，即为信号源的波达方向及波达时间：

$P(\mathrm{\θ},T)=a{(\mathrm{\θ},T)}^{*}{U}_{N^2}{U_{N^2}}^{*}a(\mathrm{\θ},T)$

其中：

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

$b(\mathrm{\θ},T)=\left[\begin{array}{c}e\left({\mathrm{\υ}}_1\right)e^{-2{\mathrm{j\π\υ}}_{1}T}d({\mathrm{\υ}}_1,\mathrm{\θ})\\ \begin{array}{c}.\\ .\\ .\end{array}\\ e\left({\mathrm{\υ}}_F\right)e^{-2{\mathrm{j\π\υ}}_{F}T}d({\mathrm{\υ}}_F,\mathrm{\θ})\end{array}\right]$

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

θ表示信号路径的波达角度；T表示信号路径的波达时间；e(υ_i)表示信号在频率υ_i处的幅值，i＝1,2,...,F；τ_1,j(θ)表示信号路径到达第j个传感器相对于到达参考传感器的时间延迟，j＝2,…,M-1。

3.如权利要求1所述波达估计方法，其特征在于，阈值ε的取值范围为1e-3～1e-6。

4.如权利要求1所述波达估计方法，其特征在于，所述空域-频域平滑中的频域平滑使用频域子带平均方法。

5.一种基于声线传播时间层析的海洋声层析方法，利用声音在海洋中传播速度的变化来反演海洋环境参数，其特征在于，首先利用权利要求1～4任一项所述基于高阶累积量的波达估计方法对从声传感器阵列所接收到的多路径传播声信号进行波达估计，从而分离出每一条声线路径；然后根据声线路径的到达时间反演出海洋环境参数。

6.一种定位方法，其特征在于，首先利用权利要求1～4任一项所述基于高阶累积量的波达估计方法进行波达方向估计，然后利用估计出的波达方向确定信号源的位置。