CN103513238B - 一种规整化最小二乘子空间相交的目标方位测向方法 - Google Patents

Abstract

本发明涉及一种规整化最小二乘子空间相交的目标方位测向方法，该方法将吉洪诺夫规整化方法与最小二乘子空间相交的算法结合起来，处理浅海声纳装置的接收信号，在利用子空间相交的原理构造成最小二乘问题的基础上，再通过规整化方法确定最优的规整化因子，最终，实现对浅海目标的检测和精确定位。本发明将Tikhonov规整化方法与最小二乘子空间相交的算法结合起来，形成新的测向算法，该方法比常规波束形成算法具有更高的方位分辨力和准确度，并且数值计算稳定。

Description

一种规整化最小二乘子空间相交的目标方位测向方法

技术领域

本发明属于声纳数字信号处理领域，特别涉及一种适用于浅海环境中的规整化最小二乘子空间相交测向方法。

背景技术

在浅海环境中，声传播受海洋波导界面的影响显著，特别是海底的影响，在这种情况下，传统的以平面波模型为基础的目标方位估计算法已不再适用，简正波模型可以更准确地描述声场，根据简正波理论，声源在海洋波导中会激发出若干号简正波，阵列的各阵元接收到的信号是各号简正波迭加的结果，各号简正波的相速度是不同的，特别是对于低频信号，相速度的差别尤为明显，这种差别在时域上反映为多途效应，采用平面波假设必然会带来目标方位估计偏差。匹配场处理（MFP）测向方法对声场用简正波模型表示，消除了测向算法中多途效应的影响，与常规平面波波束形成法相比，有明显的优点。

MFP测向的思想出现较早，因巨大的建模和信号处理计算量，未能在工程中使用。早在上世纪90年代，为了减小浅海低频测向的误差，就有人研究了相干特性对大尺度阵列的匹配场测向的性能影响，如文献1“P.S.Naidu,“Onsubspacemethodforsourcelocalization,”JASA,90(5),1991:2489-2491”；后来子空间的思想被用于匹配测向，如文献2“J.V.Candy，“Oceanacousticsignalprocessing:amodel-basedapproach”,JASA,1992”；同一频率的声波，各模态相速度不同，参考声速的选取不当导致阵列在端射附近出现明显的测向误差，文献3“宫在晓，林京，郭良浩，“浅海声传播相速度对测向精度的影响”，声学学报，27(6),2002:492-496”提出用声场匹配的方法做测向，同时指出匹配场测向方法需要知道海洋环境参数，而且需要对三维空间进行极大值搜索，计算量很大，计算量和环境失配是要面临的问题。

为解决MFP测向计算量大的问题，文献4“S.Lakshmipathi,andG.V.Anand,“Subspaceintersectionmethodofhigh-resolutionbearingestimationinshallowocean”,SignalProcessing,84,2004:1367-1384”提出了一种子空间相交（Subspaceintersection，SI）的算法，在已知各号简正波的波数的情况下，只需要对角度进行搜索，可以以较小的计算量实现匹配测向，较好地解决了建模计算量的问题。SI算法使用QR分解，在要分解的矩阵条件数很大时，数值计算不稳定；文献5“张爱民，基于子空间的目标方位估计的若干问题研究，中科院声学所博士学位论文”将子空间相交问题描述为最小二乘问题，提出了最小二乘SI的算法，并提出了总体最小二乘的算法，使算法在存在系数的扰动时稳健性能提高。

然而，正如所有的最小二乘问题一样，上述最小二乘SI算法中的方程的求解面临着病态问题，通常情况下，不同号数简正波的波数相差很小，在有限阵长条件下，很难满足在每个扫描向量线性无关的条件，矩阵的条件数很大，即该最小二乘问题的方程是病态的，数值计算的稳定性很差。

总体而言，最小二乘SI算法比传统波束形成算法具有无可比拟的性能优势，但数值计算稳定性差，测向精度低，可行性不高。实际应用中需要一种能够在模型参数扰动的情形下依然可用的处理方法。

发明内容

本发明目的在于，为克服现有最小二乘SI测向方法的病态问题，提出了一种规整化最小二乘子空间相交测向方法，从而使在浅海环境中的测向精度和数值计算稳定度大大提高。

为实现上述发明目的，本发明的规整化最小二乘子空间相交的测向方法，该方法将吉洪诺夫（Tikhonov）规整化方法与最小二乘子空间相交的算法结合起来，形成新的测向算法，其中Tikhonov规整化方法在文献6“Chne.H.Golub,etal.Tikhnovregularizationandtotleleastsquare”中有详细论述。该方法比常规波束形成算法具有更高的方位分辨力和准确度，并且数值计算稳定。

本发明的一种规整化最小二乘子空间相交的目标方位测向方法，该方法将吉洪诺夫规整化方法与最小二乘子空间相交的算法结合起来，处理浅海声纳装置的接收信号，在利用子空间相交的原理构造成最小二乘问题的基础上，再通过规整化方法确定最优的规整化因子，最终，实现对浅海目标的检测和精确定位。

所述的浅海声纳装置由多个水听器组成，是一条海底直线阵或是拖曳阵。

该方法的步骤具体包括：

1）所述的浅海声纳装置由多个水听器组成，数目为N，阵元间距为d，目标入射方向θ，快拍长度为L；用线阵接收空间信号，得到N个阵元的时域信号x(t)；

2）对L个时域快拍的数据进行快速傅里叶变换：

$x_{k,m}\left(t\right)\⇒X_{k,m}\left(f_i\right)$

3）根据傅里叶变换后的快拍数据估计阵列协方差矩阵R：

R=E[XX^H]

上式中，X是阵列接收数据矩阵，其中元素为步骤2）中的X_k,m(f_i)，H表示共轭转置，E表示对不同的快拍求统计平均；

4）对阵列协方差矩阵R进行特征值分解：

R=U∑U

式中，U为特征矢量矩阵；

其中，由特征值组成的对角阵Σ为：

如果有J个声源，则上式中的特征值满足如下关系：

λ₁≥λ₂≥…λ_J≥λ_J+1=…λ_N=σ²

上式中，λ表示阵列协方差的特征值，σ表示噪声的标准差；

由此，定义两个对角阵：

其中，Σ_S为大特征值组成的对角阵，Σ_N为小特征值组成的对角阵；

再将特征矢量矩阵分为与特征值对应的两部分：

大特征值对应的信号子空间：

U_S=[u₁,u₂,…,u_J]

和小特征值对应的噪声子空间：

U_N=[u_J+1,u_J+2,…,u_N]

5）根据海洋环境建模得到简正波响应矩阵：

A(θ)=[a(θ，k₁),...,a(θ，k_M)]∈C^N×M

上式中，A表示简正波响应矩阵，k_i表示第i号简正波的波数，M表示简正波号数，a(θ，k_i)表示当目标方位为θ时，接收阵列相对于第i号简正波的导向矢量；

6）利用子空间相交的原理，若某一方位存在目标，利用信号子空间和简正波响应矩阵构造成最小二乘问题：

A(θ)x=U_sy

上式中，A表示简正波响应矩阵，U_S表示大特征值对应的特征向量，x，y代表最小二乘问题的未知向量；

7）通过规整化方法，遍历所有可行的规整化因子，将规整化解向量的范数与对应的余差作出问题的L曲线，确定最优的规整化因子λ；

8）最后，得到问题的规整化最小二乘解：

x=(A^TA+λI)^-1ATU_sy

继而得到目标方位的判别量：

$P\left(\mathrm{\θ}\right)=\frac{1}{|A\left(\mathrm{\θ}\right)x-U_{s}y|};$

上式中，P(θ)表示θ方位的判别量，将θ从1~360遍历，可得到全方位的判别量，其中最大值对应的方位即为目标方位。

本发明的优点在于：本发明将Tikhonov规整化方法与最小二乘子空间相交的算法结合起来，形成新的测向算法，该方法比常规波束形成算法具有更高的方位分辨力和准确度，并且数值计算稳定。

附图说明

图1是本发明的算法流程图；

图2是本发明实施例中的解向量示意图；

图3是本发明实施例中的海试数据处理得到的L曲线示意图；

图4是采用本发明、未采用本发明处理得到噪声目标波束输出对比示意图；

图5是采用本发明、未采用本发明处理得到宽带声源目标波束输出对比示意图。

具体实施方式

下面结合附图对本发明进行进一步说明。

本发明的基本构思是：实际应用系统中，所有的最小二乘问题的求解都面临着病态问题，最小二乘SI算法也不例外。通常情况下，不同号数简正波的波数相差很小，在有限阵长条件下，很难满足在每个扫描向量线性无关的条件，矩阵的条件数很大，即该最小二乘问题的方程是病态的，数值计算的稳定性很差。

为实现精确测向的目的，设实际阵元数目为N，阵元间距d，目标入射方向θ，阵元接收信号表示为x(t)；快拍长度为L。本发明提出一种规整化最小二乘子空间相交的测向方法包含如下步骤：

1）用线阵接收空间信号，得到N个阵元的时域信号；

2）对L个时域快拍的数据做快速傅里叶变换；

3）使用傅里叶变换后的快拍数据估计阵列协方差矩阵R

R=E[XX^H]

4）对R进行特征值分解如下

R=U∑U

式中，U为特征矢量矩阵，其中由特征值组成的对角阵Σ如下：

上式中的特征值满足如下关系：

λ₁≥λ₂≥…λ_J≥λ_J+1=…λ_N=σ²

定义如下两个对角阵：

前者为大特征值组成的对角阵，后者为小特征值组成的对角阵。

将特征矢量矩阵分为与特征值对应的两部分：一是与大特征值对应的信号子空间

U_S=[u₁,u₂,…,u_J]

二是小特征值对应的噪声子空间

U_N=[u_J+1,u_J+2,…,u_N]

5）根据海洋环境建模得到简正波响应矩阵

A(θ)=[a(θ，k₁),...,a(θ，k_M)]∈C^N×M

6）利用子空间相交的原理，若某一方位存在目标，则信号子空间式及简正波响应矩阵式可以构造成最小二乘问题

A(θ)x=U_sy

7）作出问题的L曲线，确定最优的规整化因子λ；

L曲线是遍历所有可行的规整化因子，将规整化解向量的范数与对应的余差画成一条曲线，将该曲线画在对数坐标轴上，它往往显示为L型，这是它名字的由来。是规整化方法中的常规手段。

8）得到问题的规整化最小二乘解：

x=(A^TA+λI)^-1A^TU_sy

继而得到目标方位的判别量：

$P\left(\mathrm{\θ}\right)=\frac{1}{|A\left(\mathrm{\θ}\right)x-U_{s}y|}$

上述步骤5）中的建模操作，在文献“M.Porter，Anumericalmethodforcomputingoceanacousticnormalmodes.PHD1984”中有详细的阐述。对本领域的技术人员来说，理解和实现是没问题的。

下面结合某次海试数据和附图对本发明的具体实施方式做进一步的详细描述。

试验参数：海深30m，46个阵元，2m、4m不等间隔水平直线阵，目标位于声阵的20°方位左右，信号采样率fs=2048Hz。试验中采用了两种声源信号，分别是试验船的辐射噪声信号和由信号发生器发出频率700-800Hz的宽带信号，其中，试验船辐射噪声的处理频带为150-300Hz，声速1512m/s，快拍长度60。

本发明的规整化最小二乘子空间相交测向方法，如图1所示，包括以下步骤：

步骤1：对应图1中的101，用线阵接收空间信号，得到46个阵元的时域信号；

步骤2：对应图1中的102和103，对60个时域快拍的数据做快速傅里叶变换，如下式所示，行表示时间采样，列表示阵元。

步骤3：对应图1中的104，使用傅里叶变换后的快拍数据估计阵列协方差矩阵R

R=E[XX^H]

步骤4：对应图1中的105，对R进行特征值分解如下

R=U∑U

大特征值对应的空间为信号子空间，在本例中，只有1个声源，故J=1，得到

U_S=[u₁]

步骤5：对应图1中的106，根据海洋环境建模得到简正波响应矩阵，根据本例中的海洋环境，使用Kraken建模工具，得到波导所能承载的简正波，共有10号简正波，得到阵列的简正波响应矩阵：

A(θ)=[a(θ，k₁),...,a(θ，k₁₀)]∈C^46×10

步骤6：对应图1中的107，利用信号子空间式及简正波响应矩阵式构造最小二乘问题，在本例中，方程简化为：

A(θ)x=u₁

图2是本方程采用原算法的解向量示意图，可以看到，该方程呈现明显的病态，系数矩阵的大小特征值相差很大（图中的三角），约差5个数量级；而系数并不随特征值减小（图中的星），这样造成的结果是在最终解当中小特征值所对应的分量占主导地位（图中的圈），所以该方程的解受噪声和计算取舍误差的影响很大。

步骤7：对应图1中的108，作出问题的L曲线，如图3所示，根据此图确定最优的规整化因子λ=0.31168；

步骤8：对应图1中的109，得到问题最终的的规整化最小二乘解：

x=(A^TA+λI)^-1A^Tu₁

继而得到目标方位的判别量：

$P\left(\mathrm{\θ}\right)=\frac{1}{|A\left(\mathrm{\θ}\right)x-u_1|}.$

图4是应用本发明前后的船噪声目标的测向结果比较，菱形连线是原方法的处理结果，三角连线是本发明的处理结果，可以看到规整化之后波束主瓣宽度明显改善，平台效应消除，另外旁瓣也较规整化之前低，表明本发明对于原方法是由明显优势的。

图5是700-800Hz的宽带声源数据的处理结果，菱形连线是原方法的处理结果，三角连线是本发明的处理结果。同样，应用本方法后测向性能得到改善。

总之，本发明能够准确稳定地得到浅海中目标方位估计值。

最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (1)

1.一种规整化最小二乘子空间相交的目标方位测向方法，该方法将吉洪诺夫规整化方法与最小二乘子空间相交的算法结合起来，处理浅海声纳装置的接收信号，在利用子空间相交的原理构造成最小二乘问题的基础上，再通过规整化方法确定最优的规整化因子，最终，实现对浅海目标的检测和精确定位；

所述的浅海声纳装置由多个水听器组成，是一条海底直线阵或是拖曳阵；

所述方法的步骤包括：

1)所述的浅海声纳装置由多个水听器组成，数目为N，阵元间距为d，目标入射方向θ，快拍长度为L；用线阵接收空间信号，得到N个阵元的时域信号x(t)；

2)对L个时域快拍的数据进行快速F傅F里T叶变换：

$x_{k,m}\left(t\right)\⇒X_{k,m}\left(f_i\right)$

3)根据傅里叶变换后的快拍数据估计阵列协方差矩阵R：

R＝E[XX^H]

上式中，X是阵列接收数据矩阵，其中元素为步骤2)中的X_k,m(f_i)，H表示共轭转置，E表示对不同的快拍求统计平均；

4)对阵列协方差矩阵R进行特征值分解：

R＝UΣU

式中，U为特征矢量矩阵；

其中，由特征值组成的对角阵Σ为：

如果有J个声源，则上式中的特征值满足如下关系：

λ₁≥λ₂≥…λ_J≥λ_J+1＝…λ_N＝σ²

上式中，σ表示噪声的标准差；

由此，定义两个对角阵：

其中，Σ_S为大特征值组成的对角阵，Σ_N为小特征值组成的对角阵；

再将特征矢量矩阵分为与特征值对应的两部分：

大特征值对应的信号子空间：

U_S＝[u₁,u₂,…,u_J]

和小特征值对应的噪声子空间：

U_N＝[u_J+1,u_J+2,…,u_N]

5)根据海洋环境建模得到简正波响应矩阵：

A(θ)＝[a(θ,k₁),...,a(θ,k_M)]∈C^N×M

上式中，A表示简正波响应矩阵，k_i表示第i号简正波的波数，M表示简正波号数，a(θ,k_i)表示当目标方位为θ时，接收阵列相对于第i号简正波的导向矢量；

6)利用子空间相交的原理，若某一方位存在目标，利用信号子空间和简正波响应矩阵构造成最小二乘问题：

A(θ)x＝U_sy

上式中，A表示简正波响应矩阵，U_S表示大特征值对应的特征向量，x，y代表最小二乘问题的未知向量；

7)通过规整化方法，遍历所有可行的规整化因子，将规整化解向量的范数与对应的余差作出问题的L曲线，确定最优的规整化因子λ；

8)最后，得到问题的规整化最小二乘解：

x＝(A^TA+λI)^-1A^TU_sy

继而得到目标方位的判别量：

$P\left(\mathrm{\θ}\right)=\frac{1}{|A\left(\mathrm{\θ}\right)x-U_{s}y|};$

上式中，P(θ)表示θ方位的判别量，将θ从1～360遍历，可得到全方位的判别量，其中最大值对应的方位即为目标方位。