CN105607033A - 基于正交均匀线阵的水下波达方向估计方法及*** - Google Patents

Abstract

本发明公开了一种基于正交均匀线阵的水下波达方向估计方法及***，所述估计***包括数据采集处理与控制模块、发射模块、接收模块、输出模块和电源模块，其中，所述接收模块中超声波接收探头阵列为正交的均匀直线阵列，其中每个阵列包括多个超声接收探头，并且数量相同。本发明实现声波在水中传播速度不确定的情况下对目标进行精确定位，对ESPRIT算法进行了改进，在算法中消去了声速，使得最后的估计结果和声速无关，从而排除了声速的不确定对波达方向估计的影响，使得其实用性更高，并且可在不同水声环境下进行测量，水下目标定位的精确性也更高。

Description

基于正交均匀线阵的水下波达方向估计方法及***

技术领域

本发明涉及目标定位的技术领域，特别涉及一种基于正交均匀线阵的水下波达方向估计方法及***。

背景技术

空间信号的波达方向(简称：DOA)估计问题是阵列信号处理所研究的一类重要问题。所谓阵列信号处理就是指在空间的不同位置放置多个传感器组成阵列，利用阵列对空间信号进行接收和处理，目的是抑制无用的干扰和噪声，提取有用的信号特征或信息。DOA估计便是用这种处理方法来估计信号的入射角。

现有的DOA估计方法主要有波束形成法、高分辨率谱估计法、子空间拟合类算法和子空间分解类算法。第一，波束形成法将传统时域傅立叶谱估计扩展成相应的空域形式，即通过空域各阵列的接收数据代替传统时域处理中的时域数据。但扩展至空域后，阵列的角度分辨力受到空域“傅立叶限”的限制。空域“傅立叶限”就是阵列的物理孔径限，即对位于一个波束宽度内的空间目标不可分辨。所以，提高空域处理精度的有效方法就是增大天线孔径(等效于减小波束宽度)，但对于许多实际应用环境而言，增大天线孔径是不现实的。第二，高分辨率谱估计法是一种空域谱估计方法，它由许多时域非线性谱估计方法推广而成。它假定信号源在空间中是连续分布的，信号是空间平稳的随机过程。但这样的假设在大多数空间谱估计中是不成立的，因而这些方法具有局限性；并且这种方法不能有效利用加性噪声的统计特性，因而其分辨性能较差。第三，子空间拟合类算法将各种不同的方位估计方法用统一的算法结构联系起来，使其协方差矩阵的估计误差达到最小。这种算法精确度高，分辨率强，但是一种迭代算法，运算量非常大。第四，子空间分解类算法通过对阵列接收数据的数学分解(如奇异值分解、特征分解及QR分解等)，将接收数据划分为两个相互正交的子空间：一个与信号源的阵列流型空间一致的信号子空间和一个与信号子空间正交的噪声子空间。它利用两个子空间的正交性可以大大提高算法的分辨力。典型代表算法为ESPRIT算法。ESPRIT算法建立在子空间旋转不变技术的基础上，不需要全空间搜索，减少了运算量。因此，ESPRIT算法具有现实可行、分辨率高和运算量较小的特点，从而克服了前三种方法存在的不足，更具优越性。但是ESPRIT算法假定信号在介质中的传播速度是确定不变的，没有考虑速度对算法结果的影响，在速度变化较大、不确定的情况下容易产生较大的误差，大大降低了目标定位的精确度。

发明内容

本发明的目的在于克服现有技术的缺点与不足，提供一种基于正交均匀线阵的水下波达方向估计方法及***，实现声波在水中传播速度不确定的情况下对目标进行精确定位，从而提高水下目标定位的精确性，且该方法分辨力高，运算量适中。

根据公开的实施例，本发明第一方面提出了一种基于正交均匀线阵的水下波达方向估计方法，所述估计方法包括下列步骤：

S1、建立具有正交均匀线阵的超声波接收探头阵列，所述正交均匀线阵为水平均匀线阵和竖直均匀线阵；

S2、推导所述水平均匀线阵的接收信号模型；

S3、计算所述水水平均匀线阵和所述水竖直均匀线阵的接收信号表达式；

S4、分别计算所述水平均匀线阵的两个协方差矩阵R_hxx和R_hxy以及所述竖直均匀线阵的两个协方差矩阵R_vxx和R_vxy，对上述协方差矩阵进行特征值分解，利用求出的最小特征值得出C_hxx和C_hxy、C_vxx和C_vxy，分别计算矩阵束{C_hxx,C_hxy}和{C_vxx,C_vxy}的广义特征值分解，得到特征值λ_x1,λ_x2,…,λ_xK和λ_y1,λ_y2,…,λ_yK；

S5、建立声波在所述水平均匀线阵以及所述竖直均匀线阵上的阵列方向角θ_xi和θ_yi之间的关系式；

S6、根据阵列方向角θ_xi和θ_yi之间的关系式，对所述水平均匀线阵的第一旋转算子以及所述竖直均匀线阵的第二旋转算子中的对角元素进行配对；

S7、根据配对结果，并联系对角元素与θ_xi之间的关系式，计算出θ_xi的表达式。

进一步地，所述步骤S1具体为：

在水中放置两个正交的所述水平均匀线阵和所述竖直均匀线阵，其中，每个均匀线阵都有M个接收阵元，且阵元间间距为d，K个窄带目标声源分别为S₁,S₂,…,S_K，中心频率为f，上述目标声源对应于所述水平均匀线阵的方向角分别为θ_x1,θ_x2,…,θ_xK，对应于所述竖直均匀线阵的方向角分别为θ_y1,θ_y2,…,θ_yK。

进一步地，所述步骤S2具体为：

S21、以第一个阵元为参考点，则第一个阵元接收的信号为：

$x_1\left(t\right)=\underset{i=1}{\overset{K}{\mathrm{\Σ}}}{s}_i\left(t\right)+n_1\left(t\right)$

其中s_i(t)表示第i个源信号，n₁(t)表示第一个阵元上的噪声；

S22、第m个阵元在同一时刻接收到的信号为：

$x_m\left(t\right)=\underset{i=1}{\overset{K}{\mathrm{\Σ}}}{s}_i\left(t\right)e^{-j(m-1)\frac{2\mathrm{\π}}{{\mathrm{\λ}}_i}d{\mathrm{sin\θ}}_{xi}}+n_m\left(t\right),m=1,2,...,M$

其中λ_i表示第i个目标源反射回来的声波波长，n_m(t)表示第m个阵元上的噪声；

S23、将各阵元的接收信号排列成列向量形式，则所述水平均匀线阵接收的信号可用以下矢量式子表示：

X(t)＝AS(t)+N(t)(公式1)

其中，为M×K的导向矢量矩阵，

X(t)＝[x₁(t),x₂(t),…,x_M(t)]^T为M×1的接收信号矩阵，

S(t)＝[s₁(t),s₂(t),…,s_K(t)]^T为K×1的源信号矩阵，

N(t)＝[n₁(t),n₂(t),…,n_M(t)]^T为M×1的噪声矩阵。

进一步地，所述步骤S3中计算所述水平均匀线阵的接收信号表达式过程如下：

S301、将所述水平均匀线阵中的M个阵元分为两个平移矢量为d的子阵列Z_hx和Z_hy，

所述子阵列Z_hx由所述水平均匀线阵的第一到第M-1个阵元组成，则有：

x_h1(t)＝x₁(t),x_h2(t)＝x₂(t),…,x_h(M-1)(t)＝x_M-1(t)

其中，x_h1(t),x_h2(t),…,x_h(M-1)(t)分别是所述子阵列Z_hx上第一个阵元到第M-1个阵元接收到的信号；x₁(t),x₂(t),…,x_M(t)分别是水平阵列上第一个阵元到第M个阵元。

所述子阵列Z_hy由所述水平均匀线阵的第二到第M个阵元组成，则有：

y_h1(t)＝x₂(t),y_h2(t)＝x₃(t),…,y_h(M-1)(t)＝x_M(t)

其中，y_h1(t),y_h2(t),…,y_h(M-1)(t)分别是所述子阵列Z_hy上第一个阵元到第M-1个阵元接收到的信号；

S302、将所述子阵列Z_hx和Z_hy中第m个阵元的接收信号分别表示为：

$x_{hm}\left(t\right)=\underset{i=1}{\overset{K}{\mathrm{\Σ}}}{s}_i\left(t\right)a_{xm}\left({\mathrm{\θ}}_i\right)+n_{hxm}\left(t\right),m=1,2,...,M$

$y_{hm}\left(t\right)=\underset{i=1}{\overset{K}{\Σ}}{s}_i\left(t\right)e^{j\frac{2\mathrm{\π}d}{{\mathrm{\λ}}_i}{\mathrm{sin\θ}}_{xi}}{a}_{xm}\left({\mathrm{\θ}}_i\right)+n_{hym}\left(t\right),m=1,2,...,M$

其中n_hxm(t)和n_hym(t)分别为子阵Z_hx和Z_hy上第m个阵元的加性噪声；

S303、将上式改写成矢量形式：

X_h(t)＝A_xS(t)+N_hx(t)

Y_h(t)＝A_xΦ_xS(t)+N_hy(t)

其中

矩阵Φ_x为K×K的对角矩阵，又称第一旋转算子，其对角元素包含K个信号的波前在任意一个阵元偶之间的相位延迟信息，表示为：

${\mathrm{\Φ}}_x=diag\left\{e^{j2\mathrm{\π}\frac{d}{{\mathrm{\λ}}_1}{\mathrm{sin\θ}}_{x1}},e^{j2\mathrm{\π}\frac{d}{{\mathrm{\λ}}_2}{\mathrm{sin\θ}}_{x2}},\mathrm{...},e^{j2\mathrm{\π}\frac{d}{{\mathrm{\λ}}_k}{\mathrm{sin\θ}}_{xk}}\right\}$ (公式2)；

所述步骤S3中计算所述竖直均匀线阵的接收信号表达式过程如下：

S3X1、将所述竖直均匀线阵中的M个阵元分为两个平移矢量为d的子阵列Z_vx和Z_vy，

所述子阵列Z_vx由所述竖直均匀线阵的第一到第M-1个阵元组成，则有：

x_v1(t)＝y₁(t),x_v2(t)＝y₂(t),…,x_v(M-1)(t)＝y_M-1(t)

其中，x_v1(t),x_v2(t),…,x_v(M-1)(t)分别是所述子阵列Z_vx上第一个到第M-1个阵元接收到的信号；y₁(t),y₂(t),…,y_M(t)分别是竖直阵列上第一个阵元到第M个阵元。

所述子阵列Z_vy由所述竖直均匀线阵的第二到第M个阵元组成，则有：

y_v1(t)＝y₂(t),y_v2(t)＝y₃(t),…,y_v(M-1)(t)＝y_M(t)

其中，y_v1(t),y_v2(t),…,y_v(M-1)(t)分别是所述子阵列Z_vy上第一个阵元到第M-1个阵元接收到的信号；

S3X2、将所述子阵列Z_vx和Z_vy中第m个阵元的接收信号分别为：

$x_{vm}\left(t\right)=\underset{i=1}{\overset{K}{\mathrm{\Σ}}}{s}_i\left(t\right)a_{ym}\left({\mathrm{\θ}}_i\right)+n_{vxm}\left(t\right),m=1,2,...,M$

$y_{vm}\left(t\right)=\underset{i=1}{\overset{K}{\Σ}}{s}_i\left(t\right)e^{j\frac{2\mathrm{\π}d}{{\mathrm{\λ}}_i}{\mathrm{sin\θ}}_{yi}}{a}_{ym}\left({\mathrm{\θ}}_i\right)+n_{vym}\left(t\right),m=1,2,...,M$

其中n_vxm(t)和n_vym(t)分别为子阵Z_vx和Z_vy上第m个阵元的加性噪声；

S3X3、将上式写成矢量形式：

X_v(t)＝A_yS(t)+N_vx(t)

Y_v(t)＝A_yΦ_yS(t)+N_vy(t)

其中

矩阵Φ_y为K×K的对角矩阵，又称第二旋转算子，其对角元素包含了K个信号的波前在任意一个阵元偶之间的相位延迟信息，表示为：

${\mathrm{\Φ}}_y=diag\left\{e^{j2\mathrm{\π}\frac{d}{{\mathrm{\λ}}_1}{\mathrm{sin\θ}}_{y1}},e^{j2\mathrm{\π}\frac{d}{{\mathrm{\λ}}_2}{\mathrm{sin\θ}}_{y2}},\mathrm{...},e^{j2\mathrm{\π}\frac{d}{{\mathrm{\λ}}_K}{\mathrm{sin\θ}}_{yk}}\right\}$ (公式3)。

进一步地，所述步骤S4中计算矩阵束{C_hxx,C_hxy}的广义特征值分解过程如下：

S401、X_h(t)的协方差矩阵表示为：

R_hxx＝E[X_h(t)X_h ^H(t)]＝A_xR_ssA^H+σ_x ²I

其中R_ss＝E{S(t)S^H(t)}，为信源部分协方差矩阵；

S402、X_h(t)和Y_h(t)的互协方差矩阵表示为：

R_hxy＝E{X_h(t)Y_h ^H(t)}＝A_xR_ssΦ_x ^HA_x ^H+σ_x ²Z；

S403、对上述协方差矩阵以及互协方差矩阵进行特征值分解得到最小特征值为σ²，利用σ²可以得到矩阵束其中C_hxy＝R_hxy-σ_x ²Z＝A_xR_ssΦ_x ^HA_x ^H；

S404、计算所述矩阵束{C_hxx,C_hxy}的广义特征值分解，得到非零特征值λ_x1,λ_x2,…,λ_xK,因此由公式(2)可记：

$e^{j2\mathrm{\π}\frac{d}{{\mathrm{\λ}}_i}{\mathrm{sin\θ}}_{xi}}={\mathrm{\φ}}_{xi},i=1,2,...,K$ (公式4)

其中φ_xi为矩阵Φ_x上的对角元素，且φ_xi∈{λ_x1,λ_x2,…,λ_xK},i＝1,2,…,K；

所述步骤S4中计算矩阵束{C_vxx,C_vxy}的广义特征值分解过程如下：

S4X1、X_v(t)的协方差矩阵表示为：

R_vxx＝E[X_v(t)X_v ^H(t)]＝A_yR_ssA_y ^H+σ_y ²I

其中R_ss＝E{S(t)S^H(t)}，为信源部分协方差矩阵；

S4X2、X_v(t)和Y_v(t)的互协方差矩阵表示为：

R_vxy＝E{X_v(t)Y_v ^H(t)}＝A_yR_ssΦ_y ^HA_y ^H+σ_y ²Z；

S4X3、对上述协方差矩阵以及互协方差矩阵进行特征值分解得到最小特征值σ_y ²，利用σ_y ²得到矩阵束y其中C_vxy＝R_vxy-σ_y ²Z＝A_yR_ssΦ_y ^HA_y ^H；

S4X4、计算所述矩阵束{C_vxx,C_vxy}的广义特征值分解，得到非零特征值λ_y1,λ_y2,…,λ_yK,因此由公式(3)可记：

$e^{j2\mathrm{\π}\frac{d}{{\mathrm{\λ}}_i}{\mathrm{sin\θ}}_{yi}}={\mathrm{\φ}}_{yi},i=1,2,...,K$ (公式5)

其中φ_yi为矩阵Φ_y上的对角元素，且φ_yi∈{λ_y1,λ_y2,…,λ_yK},i＝1,2,…,K。

进一步地，所述步骤S5具体为

步骤S5、建立声波在所述水平均匀线阵以及所述竖直均匀线阵上的阵列方向角θ_xi和θ_yi之间的关系式；

S51、分别以所述水平均匀线阵以及所述竖直均匀线阵所在的方向建立一个直角坐标系，其中，θ_v1和θ_h1分别为第一象限中声波与两个均匀线阵法线的夹角，θ_v2和θ_h2分别为第二象限中声波与两个均匀线阵法线的夹角；

S52、当声波从第一象限中入射时，此时有θ_xi＝θ_h1，θ_yi＝θ_v1，则θ_xi和θ_yi的关系为：

${\mathrm{\θ}}_{xi}+{\mathrm{\θ}}_{yi}=\frac{\mathrm{\π}}{2}$ (公式6)；

S53、当声波从第二象限中入射时，此时有θ_xi＝-θ_h2，θ_yi＝θ_v2，则θ_xi和θ_yi的关系为：

$-{\mathrm{\θ}}_{xi}+{\mathrm{\θ}}_{yi}=\frac{\mathrm{\π}}{2}$ (公式7)；

S54、根据公式(6)和公式(7)得到：

${\mathrm{sin\θ}}_{yi}=sin(\frac{\mathrm{\π}}{2}\±{\mathrm{\θ}}_{xi})={\mathrm{cos\θ}}_{xi}$ (公式8)。

进一步地，所述步骤S6具体为

S61、将公式(8)带入公式(5)，则有：

$e^{j2\mathrm{\π}\frac{d}{{\mathrm{\λ}}_i}{\mathrm{cos\θ}}_{xi}}={\mathrm{\φ}}_{yi},i=1,2,...,K$ (公式9)；

S62、根据公式(4)和公式(9)可知，若配对成功，则有以下式子成立：

${\arg }^2\left({\mathrm{\φ}}_{xi}\right)+{\arg }^2\left({\mathrm{\φ}}_{yi}\right)=\frac{4{\mathrm{\π}}^2{d}^2}{{\mathrm{\λ}}_i^2},i=1,2,...,K$ (公式10)

其中，λ_i是第i个目标源反射回来的声波波长，且有λ_i＝v_i/f；

S63、将arg(λ_x1),arg(λ_x2),…,arg(λ_xK)按照各自的平方大小顺序从大到小排列得到序列Η；将arg(λ_y1),arg(λ_y2),…,arg(λ_yK)按照各自的平方大小顺序从小到大排列得到序列V，于是有：

$\left\{\begin{array}{c}\arg \left({\mathrm{\φ}}_{xi}\right)=h_i\\ \arg \left({\mathrm{\φ}}_{yi}\right)=v_i\end{array}\right.,i=1,2,...,K$ (公式11)

其中h_i为序列Η中的第i个元素；v_i为序列V中的第i个元素。

进一步地，所述θ_xi的表达式；

${\mathrm{\θ}}_{xi}=arcsin\left(\frac{\arg \left({\mathrm{\φ}}_{xi}\right)}{\sqrt{{\arg }^2\left({\mathrm{\φ}}_{xi}\right)+{\arg }^2\left({\mathrm{\φ}}_{yi}\right)}}\right),i=1,2...,K$ (公式12)

当arg(φ_xi)大于0时，θ_xi取正值，即目标声源在x轴负方向；当arg(φ)_xi小于0时，θ_xi取负值，即目标声源在x轴正方向。

根据公开的实施例，本发明第二方面提出了一种基于正交均匀线阵的水下波达方向估计***，所述估计***包括数据采集处理与控制模块、发射模块、接收模块、输出模块和电源模块，

其中，所述数据采集处理与控制模块包括处理器、D/A转换器和一个A/D转换器，并且分别与所述发射模块、所述接收模块、所述输出模块和所述电源模块相连，能够通过设定参数来控制所述发射模块，使所述发射模块发射指定的声波信号；还能够对所述接收模块传送回来的数据进行采集与处理，计算出精确的波达方向角；该模块还能够将处理完毕的数据传送给所述输出模块，然后通过USB接口输出到外部装置或者在显示器上显示出来。

所述发射模块包括至少一个超声波发送探头，并且与所述数据采集处理与控制模块和所述电源模块相连，能够根据所述数据采集处理与控制模块的指令，获取指定的信号并对其进行超声发射；

所述接收模块包括两个以相互垂直方式放置的超声波接收探头阵列，并与所述数据采集处理与控制模块和所述电源模块相连，能够接收目标声源发射回来的信号，并将其传送给所述数据采集处理与控制模块；

所述输出模块包括至少一个USB接口和至少一个显示器组成，并且与所述数据采集处理与控制模块模块和所述电源模块相连，能够提供人机交互，将所述数据采集处理与控制模块模块中处理好的数据通过USB接口输出到外部装置或者在显示器上显示出来；

所述电源模块与上述各模块相连，为上述各模块供电。

进一步地，所述数据采集处理与控制模块中处理器为TMS320VC5509A型号的DSP芯片，所述接收模块中超声波接收探头阵列为正交的均匀直线阵列，其中每个阵列包括多个超声接收探头，并且数量相同。

本发明相对于现有技术具有如下的优点及效果：

1、本发明相对于利用传统ESPRIT算法来进行水下目标定位的方法更具有实用性，并且精确度也更高。在复杂的水声环境中，声波的速度并不确定，即每个目标声源反射回来的声波速度大小不相同，而传统ESPRIT算法假定声速是不变且相同的，没有考虑声速的不确定性对估计算法的影响，这使得在水声环境下传统ESPRIT算法不具有实用性。本发明采用了两个正交的均匀线阵，对ESPRIT算法进行了改进，在算法中消去了声速，使得最后的估计结果和声速无关，从而排除了声速的不确定对波达方向估计的影响，使得其实用性更高，并且可在不同水声环境下进行测量，目标定位的精度也更高。

2、本发明对传统ESPRIT算法进行了改进，从而保留了ESPRIT算法分辨率高的优点，并且改进后算法的运算量和复杂度也没有增加，保证了算法的可行性。

3、本发明的估计***可行性强，安装简单。除此之外，现代处理器计算处理能力的不断提高，这使得本发明所使用的处理器等芯片的集成度高，并且计算能力强，从而保证了本发明的可行性。

附图说明

图1是本发明估计***的硬件结构模块图；

图2是本发明估计***所用的正交均匀线阵模型；

图3是本发明中水平均匀线阵的接收信号模型；

图4是本发明公开的基于正交均匀线阵的水下波达方向估计方法的流程图。

具体实施方式

为使本发明的目的、技术方案及优点更加清楚、明确，以下参照附图并举实施例对本发明进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

实施例一

本实施例提供的基于正交均匀线阵的水下波达方向估计方法，该估计方法的流程步骤图参照附图4所示，包含以下几个具体步骤：

步骤S1、建立正交均匀线阵的场景模型；

在水中放置两个正交的均匀线阵，其模型场景如图2所示。每个均匀线阵都有M个接收阵元，且阵元间间距为d，K个窄带目标声源分别为S₁,S₂,…,S_K，中心频率为f，这些目标声源对应于水平线阵的方向角分别为θ_x1,θ_x2,…,θ_xK，对应于竖直线阵的方向角分别为θ_y1,θ_y2,…,θ_yK。

步骤S2、推导水平均匀线阵的接收信号模型；

水平均匀线阵的模型场景如图3所示。以第一个阵元为参考点，则第一个阵元接收的信号为：

$x_1\left(t\right)=\underset{j=1}{\overset{K}{\mathrm{\Σ}}}{s}_i\left(t\right)+n_1\left(t\right)$

其中s_i(t)表示第i个源信号，n₁(t)表示第一个阵元上的噪声。

接收信号满足窄带条件，即当信号延迟远小于带宽倒数时，延迟作用相当于使基带信号产生一个相移。那么第m个阵元在同一时刻接收到的信号为：

$x_m\left(t\right)=\underset{i=1}{\overset{K}{\mathrm{\Σ}}}{s}_i\left(t\right)e^{-j(m-1)\frac{2\mathrm{\π}}{{\mathrm{\λ}}_i}d{\mathrm{sin\θ}}_{xi}}+n_m\left(t\right),m=1,2,...,M$

其中λ_i表示第i个目标源反射回来的声波波长，n_m(t)表示第m个阵元上的噪声。将各阵元的接收信号排列成列向量形式，则整个水平线阵接收的信号可用以下矢量式子表示：

X(t)＝AS(t)+N(t)(公式1)

其中，为M×K的导向矢量矩阵，

X(t)＝[x₁(t),x₂(t),…,x_M(t)]^T为M×1的接收信号矩阵，

S(t)＝[s₁(t),s₂(t),…,s_K(t)]^T为K×1的源信号矩阵，

N(t)＝[n₁(t),n₂(t),…,n_M(t)]^T为M×1的噪声矩阵。

步骤S3、计算水平阵列和竖直阵列的接收信号表达式；

将水平均匀线阵中的M个阵元分为两个平移矢量为d的子阵列Z_hx和Z_hy。子阵列Z_hx由水平均匀线阵的第一到第M-1个阵元组成，则有：

x_h1(t)＝x₁(t),x_h2(t)＝x₂(t),…,x_h(M-1)(t)＝x_M-1(t)

其中，x_h1(t),x_h2(t),…,x_h(M-1)(t)分别是子阵列Z_hx上第一个阵元到第M-1个阵元接收到的信号；x₁(t),x₂(t),…,x_M(t)分别是水平阵列上第一个阵元到第M个阵元。

子阵列Z_hy由水平阵列的第二到第M个阵元组成，则有：

y_h1(t)＝x₂(t),y_h2(t)＝x₃(t),…,y_h(M-1)(t)＝x_M(t)

其中，y_h1(t),y_h2(t),…,y_h(M-1)(t)分别是子阵列Z_hy上第一个阵元到第M-1个阵元接收到的信号。

那么两个子阵列中第m个阵元的接收信号分别为：

$x_{hm}\left(t\right)=\underset{i=1}{\overset{K}{\mathrm{\Σ}}}{s}_i\left(t\right)a_{xm}\left({\mathrm{\θ}}_i\right)+n_{hxm}\left(t\right),m=1,2,...,M$

$y_{hm}\left(t\right)=\underset{i=1}{\overset{K}{\Σ}}{s}_i\left(t\right)e^{j\frac{2\mathrm{\π}d}{{\mathrm{\λ}}_i}{\mathrm{sin\θ}}_{xi}}{a}_{xm}\left({\mathrm{\θ}}_i\right)+n_{hym}\left(t\right),m=1,2,...,M$

其中n_hxm(t)和n_hym(t)分别为子阵Z_hx和Z_hy上第m个阵元的加性噪声。将上式写成矢量形式：

X_h(t)＝A_xS(t)+N_hx(t)

Y_h(t)＝A_xΦ_xS(t)+N_hy(t)

其中

矩阵Φ_x为K×K的对角矩阵，它是把子阵Z_hx和Z_hy的输出联系起来的酉阵，也称第一旋转算子，其对角元素包含了K个信号的波前在任意一个阵元偶之间的相位延迟信息，表示为：

${\mathrm{\Φ}}_x=diag\left\{e^{j2\mathrm{\π}\frac{d}{{\mathrm{\λ}}_1}{\mathrm{sin\θ}}_{x1}},e^{j2\mathrm{\π}\frac{d}{{\mathrm{\λ}}_2}{\mathrm{sin\θ}}_{x2}},\mathrm{...},e^{j2\mathrm{\π}\frac{d}{{\mathrm{\λ}}_k}{\mathrm{sin\θ}}_{xk}}\right\}$ (公式2)

同样，将竖直均匀线阵中的M个阵元也分为两个平移矢量为d的子阵列Z_vx和Z_vy。子阵列Z_vx由竖直均匀线阵的第一到第M-1个阵元组成，则有：

x_v1(t)＝y₁(t),x_v2(t)＝y₂(t),…,x_v(M-1)(t)＝y_M-1(t)

其中，x_v1(t),x_v2(t),…,x_v(M-1)(t)分别是子阵列Z_vx上第一个阵元到第M-1个阵元接收到的信号。

子阵列Z_vy由竖直阵列的第二到第M个阵元组成，则有：

y_v1(t)＝y₂(t),y_v2(t)＝y₃(t),…,y_v(M-1)(t)＝y_M(t)

其中，y_v1(t),y_v2(t),…,y_v(M-1)(t)分别是子阵列Z_vy上第一个阵元到第M-1个阵元接收到的信号；y₁(t),y₂(t),…,y_M(t)分别是竖直阵列上第一个阵元到第M个阵元。

那么两个子阵列中第m个阵元的接收信号分别为：

$x_{vm}\left(t\right)=\underset{i=1}{\overset{K}{\mathrm{\Σ}}}{s}_i\left(t\right)a_{ym}\left({\mathrm{\θ}}_i\right)+n_{vxm}\left(t\right),m=1,2,...,M$

$y_{vm}\left(t\right)=\underset{i=1}{\overset{K}{\Σ}}{s}_i\left(t\right)e^{j\frac{2\mathrm{\π}d}{{\mathrm{\λ}}_i}{\mathrm{sin\θ}}_{yi}}{a}_{ym}\left({\mathrm{\θ}}_i\right)+n_{vym}\left(t\right),m=1,2,...,M$

其中n_vxm(t)和n_vym(t)分别为子阵Z_vx和Z_vy上第m个阵元的加性噪声。将上式写成矢量形式：

X_v(t)＝A_yS(t)+N_vx(t)

Y_v(t)＝A_yΦ_yS(t)+N_vy(t)

其中矩阵Φ_y为K×K的对角矩阵，它是把子阵Z_vx和Z_vy的输出联系起来的酉阵，也称第二旋转算子，其对角元素包含了K个信号的波前在任意一个阵元偶之间的相位延迟信息，表示为：

${\mathrm{\Φ}}_y=diag\left\{e^{j2\mathrm{\π}\frac{d}{{\mathrm{\λ}}_1}{\mathrm{sin\θ}}_{y1}},e^{j2\mathrm{\π}\frac{d}{{\mathrm{\λ}}_2}{\mathrm{sin\θ}}_{y2}},\mathrm{...},e^{j2\mathrm{\π}\frac{d}{{\mathrm{\λ}}_K}{\mathrm{sin\θ}}_{yk}}\right\}$ (公式3)

步骤S4、分别计算水平均匀线阵的两个协方差矩阵R_hxx和R_hxy以及竖直均匀线阵的两个协方差矩阵R_vxx和R_vxy，对上述矩阵协方差矩阵进行特征值分解，利用求出的最小特征值得出C_hxx和C_hxy、C_vxx和C_vxy，分别计算矩阵束{C_hxx,C_hxy}和{C_vxx,C_vxy}的广义特征值分解，得到特征值。

求解第一旋转算子Φ_x和第二旋转算子Φ_y。X_h(t)的协方差矩阵可以表示为：

R_hxx＝E[X_h(t)X_h ^H(t)]＝A_xR_ssA_x ^H+σ_x ²I

其中R_ss＝E{S(t)S^H(t)}，为信源部分协方差矩阵。

X_h(t)和Y_h(t)的互协方差矩阵为：

R_hxy＝E{X_h(t)Y_h ^H(t)}＝A_xR_ssΦ_x ^HA_x ^H+σ_x ²Z

对矩阵协方差矩阵进行特征值分解得到最小特征值为σ²，利用σ²可以得到矩阵束{C_hxx,C_hxy}，其中C_hxx＝R_hxx-σ_x ²I＝A_xR_ssA_x ^H，C_hxy＝R_hxy-σ_x ²Z＝A_xR_ssΦ_x ^HA_x ^H。计算矩阵束{C_hxx,C_hxy}的广义特征值分解，得到非零特征值λ_x1,λ_x2,…,λ_xK,它们一一对应着矩阵Φ_x对角线上的元素，但对应关系并不确定，因此由公式(2)可记：

$e^{j2\mathrm{\π}\frac{d}{{\mathrm{\λ}}_i}{\mathrm{sin\θ}}_{xi}}={\mathrm{\φ}}_{xi},i=1,2,...,K$ (公式4)

其中φ_xi为矩阵Φ_x上的对角元素，且φ_xi∈{λ_x1,λ_x2,…,λ_xK},i＝1,2,…,K。

根据以上步骤，同理可以求得竖直均匀线阵的两个协方差矩阵R_vxx和R_vxy，X_v(t)的协方差矩阵可以表示为：

R_vxx＝E[X_v(t)X_v ^H(t)]＝A_yR_ssA_y ^H+σ_y ²I

其中R_ss＝E{S(t)S^H(t)}，为信源部分协方差矩阵。

X_v(t)和Y_v(t)的互协方差矩阵为：

R_vxy＝E{X_v(t)Y_v ^H(t)}＝A_yR_ssΦ_y ^HA_y ^H+σ_y ²Z

其中R_ss＝E{S(t)S^H(t)}，为信源部分协方差矩阵。X_v(t)和Y_v(t)的互协方差矩阵为：对矩阵协方差矩阵进行特征值分解得到最小特征值σ_y ²，利用σ_y ²得到矩阵束{C_vxx,C_vxy}，其中C_vxx＝R_vxx-σ_y ²I＝A_yR_ssA_y ^H，C_vxy＝R_vxy-σ_y ²Z＝A_yR_ssΦ_y ^HA_y ^H。

计算矩阵束{C_vxx,C_vxy}的广义特征值分解，得到非零特征值λ_y1,λ_y2,…,λ_yK,它们一一对应着矩阵Φ_y对角线上的元素，但对应关系并不确定，因此由公式(3)可记：

$e^{j2\mathrm{\π}\frac{d}{{\mathrm{\λ}}_i}{\mathrm{sin\θ}}_{yi}}={\mathrm{\φ}}_{yi},i=1,2,...,K$ (公式5)

其中φ_yi为矩阵Φ_y上的对角元素，且φ_yi∈{λ_y1,λ_y2,…,λ_yK},i＝1,2,…,K。

步骤S5、建立声波在水平均匀线阵以及竖直均匀线阵的阵列方向角θ_xi和θ_yi之间的关系式；

在图2中分别以这两个正交均匀线阵所在的方向可以建立一个直角坐标系。因为只是在水面的位置对水面下方的目标进行定位，因此可只考虑两个象限的声波入射情况。其中，θ_v1和θ_h1分别为第一象限①中声波与两个均匀线阵法线的夹角；θ_v2和θ_h2分别为第二象限②中声波与两个均匀线阵法线的夹角。

当声波从第一象限①中入射时，此时有θ_xi＝θ_h1，θ_yi＝θ_v1，则θ_xi和θ_yi的关系为：

${\mathrm{\θ}}_{xi}+{\mathrm{\θ}}_{yi}=\frac{\mathrm{\π}}{2}$ (公式6)

当声波从第二象限②中入射时，由于处在x轴上的阵列信号是以处在x轴最负方向的阵元为参考阵元的，并且子阵Z_hx也在子阵Z_hy的负x轴方向。因此当声波从象限②中入射时，参考阵元是最晚接收到信号的，子阵Z_hx中的阵元也比子阵Z_hy中对应的阵元晚接收到信号，从而可以得到时延参数τ小于0，又因为所以此时有θ_xi＝-θ_h2，θ_yi＝θ_v2，则θ_xi和θ_yi的关系为：

$-{\mathrm{\θ}}_{xi}+{\mathrm{\θ}}_{yi}=\frac{\mathrm{\π}}{2}$ (公式7)

根据公式(6)和公式(7)可以得到：

${\mathrm{sin\θ}}_{yi}=sin(\frac{\mathrm{\π}}{2}\±{\mathrm{\θ}}_{xi})={\mathrm{cos\θ}}_{xi}$ (公式8)

将公式(8)带入公式(5)，则有：

$e^{j2\mathrm{\π}\frac{d}{{\mathrm{\λ}}_i}{\mathrm{cos\θ}}_{xi}}={\mathrm{\φ}}_{yi},i=1,2,...,K$ (公式9)

步骤S6、根据阵列方向角θ_xi和θ_yi之间的关系式，对两个阵列的旋转算子中的对角元素进行配对；

对矩阵Φ_x和矩阵Φ_y上的对角元素φ_xi与φ_yi进行配对。

将公式(8)带入公式(5)，则有：

$e^{j2\mathrm{\π}\frac{d}{{\mathrm{\λ}}_i}{\mathrm{cos\θ}}_{xi}}={\mathrm{\φ}}_{yi},i=1,2,...,K$ (公式9)

根据公式(4)和公式(9)可知，若配对成功，则有以下式子成立：

${\arg }^2\left({\mathrm{\φ}}_{xi}\right)+{\arg }^2\left({\mathrm{\φ}}_{yi}\right)=\frac{4{\mathrm{\π}}^2{d}^2}{{\mathrm{\λ}}_i^2},i=1,2,...,K$ (公式10)

其中，λ_i是第i个目标源反射回来的声波波长，且有由于每个目标源反射回来的声波速度v_i差异不大，故对于每个目标源来说，arg²(φ_xi)与arg²(φ_yi)的和都是一定的。

将arg(λ_x1),arg(λ_x2),…,arg(λ_xK)按照各自的平方大小顺序从大到小排列得到序列Η；将arg(λ_y1),arg(λ_y2),…,arg(λ_yK)按照各自的平方大小顺序从小到大排列得到序列V。于是有：

$\left\{\begin{array}{c}\arg \left({\mathrm{\φ}}_{xi}\right)=h_i\\ \arg \left({\mathrm{\φ}}_{yi}\right)=v_i\end{array}\right.,i=1,2,...,K$ (公式11)

其中h_i为序列Η中的第i个元素；v_i为序列V中的第i个元素。

步骤S7、根据配对结果，并联系对角元素与θ_xi之间的关系式，计算出θ_xi的表达式；

根据配对结果求出θ_xi的大小。

根据公式(4)和公式(9)可以约去λ_i得：

${\mathrm{\θ}}_{xi}=arcsin\left(\frac{\arg \left({\mathrm{\φ}}_{xi}\right)}{\sqrt{{\arg }^2\left({\mathrm{\φ}}_{xi}\right)+{\arg }^2\left({\mathrm{\φ}}_{yi}\right)}}\right),i=1,2...,K$ (公式12)

当arg(φ_xi)大于0时，θ_xi取正值，即目标声源在x轴负方向；当arg(φ)_xi小于0时，θ_xi取负值，即目标声源在x轴正方向。

根据以上算法流程可知，本发明提出的改进算法不需要知道声速的大小就可以对θ_xi进行精确的估计，即可以在声速不确定的情况下估计出波达方向角θ_xi的值，克服了传统ESPRIT算法的缺点。

实施例二

本实施例提供的基于正交均匀线阵的水下波达方向估计***，其硬件结构模块组成如图1所示，包括数据采集处理与控制模块、发射模块、接收模块、输出模块和电源模块。

数据采集处理与控制模块由一个处理器、一个D/A转换器和一个A/D转换器组成，并且与发射模块、接收模块、输出模块和电源模块相连。它能够通过设定参数来控制发射模块，使发射模块发射指定的声波信号；还能够对接收模块传送回来的数据进行采集与处理，运用本发明算法计算出精确的波达方向角。除此之外，它还能够将处理好的数据传送给输出模块，使它们通过USB接口输出到外部装置或者在显示器上显示出来。

发射模块由一个超声波发送探头组成，并且与数据采集处理与控制模块和电源模块相连。它能够根据数据采集处理与控制模块的指令，获取指定的信号并对其进行超声发射；

接收模块由两个以相互垂直方式放置的超声波接收探头阵列组成，并与数据采集处理与控制模块和电源模块相连。它能够接收目标声源发射回来的信号，并将其传送给数据采集处理与控制模块。

输出模块由一个USB接口和一个显示器组成，并且与数据采集处理与控制模块模块和电源模块相连。它能够提供人机交互，将数据采集处理与控制模块模块中处理好的数据通过USB接口输出到外部装置或者在显示器上显示出来。

电源模块由一个电源组成，并且与数据采集处理与控制模块、发射模块、接收模块和输出模块相连。它能够为这些模块供电。

本发明装置的主要工作流程如下：操作者通过数据采集处理与控制模块输入对应的参数，使该模块产生所需要的数字信号，该模块中的D/A转换器将其转换成模拟信号后，传送给发射模块，超声波发射探头将对其进行发射。接收模块中的垂直阵列接收到信号后将其传送给数据采集处理与控制模块，信号经过该模块中的A/D转换器数字信号后变为数字信号，然后模块中的处理器利用本发明算法对其进行计算。最后数据采集处理与控制模块将计算得出的结果传给输出模块，使结果通过USB接口输出给外部装置或者在显示屏上显示出来。工作的全程所有模块都由电源模块供电。

实施例三

本实施例提供的基于正交均匀线阵的水下波达方向估计***，包括数据采集处理与控制模块、发射模块、接收模块、输出模块和电源模块。数据采集处理与控制模块采用DSP芯片实现(如：TI公司的TMS320VC5509A型号的DSP芯片)，此DSP芯片可实现A/D转换和D/A转换的功能，并能够实现均匀线阵的旋转算子和最终波达方向的计算；发射模块使用一个超声波发射探头；接收模块使用两个正交的均匀直线阵列，其中每个阵列包括多个超声接收探头，并且数量相同；输出模块使用一个USB接口和一个LCD显示屏。图1即为本实施例所述装置的硬件结构模块图。

本发明的主要工作步骤具体如下：

步骤1：在水下放置4个目标声源，与水平阵列的方向角分别为35°，60°，-30°和-70°。设定阵列接收模块和发射模块中的各项参数。将接收模块中的每个均匀线阵中的阵元个数M定为9，则每个均匀线阵中前8个阵元为子阵Z_x，后8个阵元为子阵Z_y。两个子阵列之间的距离d设为5mm。通过数据采集处理与控制模块控制发射模块使待发射的超声频率为100kHz，脉冲长度5ms。

步骤2：对超声接收探头接收到的目标声源信号进行采样；水平方向均匀线阵接收到的信号分别为x_x1(t),x_x2(t),…,x_x8(t)和y_x1(t),y_x2(t),…,y_x8(t)，竖直方向均匀线阵接收的信号分别为x_y1(t),x_y2(t),…,x_y8(t)和y_y1(t),y_y2(t),…,y_y8(t)。共采样接收200次，并将接收到的信号传递给数据采集处理与控制模块模块做运算处理。

步骤3：信号在数据采集处理与控制模块中的处理步骤具体如下：

31)将处在水平方向上的均匀线阵接收到的信号排成矢量形式X_h(t)和Y_h(t)，计算X_h(t)的协方差矩阵R_hxx＝E[X_h(t)X_h ^H(t)]，X_h(t)和Y_h(t)之间的互协方差矩阵R_hxy＝E{X_h(t)Y_h ^H(t)}。同时对竖直方向上的均匀线阵接收到的信号也进行相同处理，得到R_vxx＝E[X_v(t)X_v ^H(t)]和R_vxy＝E{X_v(t)Y_v ^H(t)}。

32)对水平阵列中的两个协方差矩阵R_hxx和R_hxy进行特征值分解，得到最小的特征值σ²，从而有C_hxx＝R_hxx-σ²I＝AR_ssA^H和C_hxy＝R_hxy-σ²Z＝AR_ssΦ^HA^H。同时对竖直阵列中的两个协方差矩阵进行相同的处理，得到C_vxx和C_vxy。

33)分别计算矩阵束{C_hxx,C_hxy}和{C_vxx,C_vxy}的广义特征值分解，得到λ_x1,λ_x2,…,λ_xK和λ_y1,λ_y2,…,λ_yK。

34)将arg(λ_x1),arg(λ_x2),…,arg(λ_xK)按照各自的平方大小顺序从大到小排列得到序列Η，将arg(λ_y1),arg(λ_y2),…,arg(λ_yK)按照各自的平方大小顺序从小到大排列得到序列V。然后把H中第i个元素h_i的值赋给arg(φ_xi)，V中第i个元素v_i的值赋给arg(φ_yi)。

34)根据匹配得出的arg(φ_xi)和arg(φ_yi)，联立式和式对θ_xi进行求解，最终求得

${\mathrm{\θ}}_{xi}=arcsin\left(\frac{\arg \left({\mathrm{\φ}}_{xi}\right)}{\sqrt{{\arg }^2\left({\mathrm{\φ}}_{xi}\right)+{\arg }^2\left({\mathrm{\φ}}_{yi}\right)}}\right).$

步骤4：将计算出的方向角信息存储下来，并传送给输出模块，使其通过USB接口输出给外部装置或者显示在LCD显示屏上。根据本发明算法估计出的方向角分别为35°，60°，-30°和-70°，与实际角度相同，说明估计结果正确，本发明方法及***可行。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (10)

1.一种基于正交均匀线阵的水下波达方向估计方法，其特征在于，所述估计方法包括下列步骤：

S1、建立具有正交均匀线阵的超声波接收探头阵列，所述正交均匀线阵为水平均匀线阵和竖直均匀线阵；

S2、推导所述水平均匀线阵的接收信号模型；

S3、计算所述水平均匀线阵和所述竖直均匀线阵的接收信号表达式；

S4、分别计算所述水平均匀线阵的两个协方差矩阵R_hxx和R_hxy以及所述竖直均匀线阵的两个协方差矩阵R_vxx和R_vxy，对上述协方差矩阵进行特征值分解，利用求出的最小特征值得出C_hxx和C_hxy、C_vxx和C_vxy，分别计算矩阵束{C_hxx,C_hxy}和{C_vxx,C_vxy}的广义特征值分解，得到特征值λ_x1,λ_x2,…,λ_xK和λ_y1,λ_y2,…,λ_yK；

S5、建立声波在所述水平均匀线阵以及所述竖直均匀线阵上的阵列方向角θ_xi和θ_yi之间的关系式；

S6、根据阵列方向角θ_xi和θ_yi之间的关系式，对所述水平均匀线阵的第一旋转算子以及所述竖直均匀线阵的第二旋转算子中的对角元素进行配对；

S7、根据配对结果，并联系对角元素与θ_xi之间的关系式，计算出θ_xi的表达式。

2.根据权利要求1所述的一种基于正交均匀线阵的水下波达方向估计方法，其特征在于，所述步骤S1具体为：

在水中放置两个正交的所述水平均匀线阵和所述竖直均匀线阵，其中，每个均匀线阵都有M个接收阵元，且阵元间间距为d，K个窄带目标声源分别为S₁,S₂,…,S_K，中心频率为f，上述目标声源对应于所述水平均匀线阵的方向角分别为θ_x1,θ_x2,…,θ_xK，对应于所述竖直均匀线阵的方向角分别为θ_y1,θ_y2,…,θ_yK。

3.根据权利要求2所述的一种基于正交均匀线阵的水下波达方向估计方法，其特征在于，所述步骤S2具体为：

S21、以第一个阵元为参考点，则第一个阵元接收的信号为：

$x_1\left(t\right)=\underset{i=1}{\overset{K}{\mathrm{\Σ}}}{s}_i\left(t\right)+n_1\left(t\right)$

其中s_i(t)表示第i个源信号，n₁(t)表示第一个阵元上的噪声；

S22、第m个阵元在同一时刻接收到的信号为：

$x_m\left(t\right)=\underset{i=1}{\overset{K}{\Σ}}{s}_i\left(t\right)e^{-j(m-1)\frac{2\mathrm{\π}}{{\mathrm{\λ}}_i}d{\mathrm{sin\θ}}_{xi}}+n_m\left(t\right),m=1,2,...,M$

其中λ_i表示第i个目标源反射回来的声波波长，n_m(t)表示第m个阵元上的噪声；

S23、将各阵元的接收信号排列成列向量形式，则所述水平均匀线阵接收的信号可用以下矢量式子表示：

X(t)＝AS(t)+N(t)(公式1)

其中，为M×K的导向矢量矩阵，

X(t)＝[x₁(t),x₂(t),…,x_M(t)]^T为M×1的接收信号矩阵，

S(t)＝[s₁(t),s₂(t),…,s_K(t)]^T为K×1的源信号矩阵，

N(t)＝[n₁(t),n₂(t),…,n_M(t)]^T为M×1的噪声矩阵。

4.根据权利要求3所述的一种基于正交均匀线阵的水下波达方向估计方法，其特征在于，所述步骤S3中计算所述水平均匀线阵的接收信号表达式过程如下：

S301、将所述水平均匀线阵中的M个阵元分为两个平移矢量为d的子阵列Z_hx和Z_hy，

所述子阵列Z_hx由所述水平均匀线阵的第一到第M-1个阵元组成，则有：

x_h1(t)＝x₁(t),x_h2(t)＝x₂(t),…,x_h(M-1)(t)＝x_M-1(t)

其中，x_h1(t),x_h2(t),…,x_h(M-1)(t)分别是所述子阵列Z_hx上第一个阵元到第M-1个阵元接收到的信号；x₁(t),x₂(t),…,x_M(t)分别是水平阵列上第一个阵元到第M个阵元；

所述子阵列Z_hy由所述水平均匀线阵的第二到第M个阵元组成，则有：

y_h1(t)＝x₂(t),y_h2(t)＝x₃(t),…,y_h(M-1)(t)＝x_M(t)

其中，y_h1(t),y_h2(t),…,y_h(M-1)(t)分别是所述子阵列Z_hy上第一个阵元到第M-1个阵元接收到的信号；

S302、将所述子阵列Z_hx和Z_hy中第m个阵元的接收信号分别表示为：

$x_{hm}\left(t\right)=\underset{i=1}{\overset{K}{\mathrm{\Σ}}}{s}_i\left(t\right)a_{xm}\left({\mathrm{\θ}}_i\right)+n_{hxm}\left(t\right),m=1,2,...,M$

$y_{hm}\left(t\right)=\underset{i=1}{\overset{K}{\Σ}}{s}_i\left(t\right)e^{j\frac{2\mathrm{\π}d}{{\mathrm{\λ}}_i}{\mathrm{sin\θ}}_{xi}}{a}_{xm}\left({\mathrm{\θ}}_i\right)+n_{hym}\left(t\right),m=1,2,...,M$

其中n_hxm(t)和n_hym(t)分别为子阵Z_hx和Z_hy上第m个阵元的加性噪声；

S303、将上式改写成矢量形式：

X_h(t)＝A_xS(t)+N_hx(t)

Y_h(t)＝A_xΦ_xS(t)+N_hy(t)

其中

矩阵Φ_x为K×K的对角矩阵，又称第一旋转算子，其对角元素包含K个信号的波前在任意一个阵元偶之间的相位延迟信息，表示为：

所述步骤S3中计算所述竖直均匀线阵的接收信号表达式过程如下：

S3X1、将所述竖直均匀线阵中的M个阵元分为两个平移矢量为d的子阵列Z_vx和Z_vy，

所述子阵列Z_vx由所述竖直均匀线阵的第一到第M-1个阵元组成，则有：

x_v1(t)＝y₁(t),x_v2(t)＝y₂(t),…,x_v(M-1)(t)＝y_M-1(t)

其中，x_v1(t),x_v2(t),…,x_v(M-1)(t)分别是所述子阵列Z_vx上第一个到第M-1个阵元接收到的信号；y₁(t),y₂(t),…,y_M(t)分别是竖直阵列上第一个阵元到第M个阵元；

所述子阵列Z_vy由所述竖直均匀线阵的第二到第M个阵元组成，则有：

y_v1(t)＝y₂(t),y_v2(t)＝y₃(t),…,y_v(M-1)(t)＝y_M(t)

其中，y_v1(t),y_v2(t),…,y_v(M-1)(t)分别是所述子阵列Z_vy上第一个阵元到第M-1个阵元接收到的信号；

S3X2、将所述子阵列Z_vx和Z_vy中第m个阵元的接收信号分别为：

$x_{vm}\left(t\right)=\underset{i=1}{\overset{K}{\mathrm{\Σ}}}{s}_i\left(t\right)a_{ym}\left({\mathrm{\θ}}_i\right)+n_{vxm}\left(t\right),m=1,2,...,M$

$y_{vm}\left(t\right)=\underset{i=1}{\overset{K}{\Σ}}{s}_i\left(t\right)e^{j\frac{2\mathrm{\π}d}{{\mathrm{\λ}}_i}{\mathrm{sin\θ}}_{yi}}{a}_{ym}\left({\mathrm{\θ}}_i\right)+n_{vym}\left(t\right),m=1,2,...,M$

其中n_vxm(t)和n_vym(t)分别为子阵Z_vx和Z_vy上第m个阵元的加性噪声；

S3X3、将上式写成矢量形式：

X_v(t)＝A_yS(t)+N_vx(t)

Y_v(t)＝A_yΦ_yS(t)+N_vy(t)

其中

矩阵Φ_y为K×K的对角矩阵，又称第二旋转算子，其对角元素包含了K个信号的波前在任意一个阵元偶之间的相位延迟信息，表示为：

5.根据权利要求4所述的一种基于正交均匀线阵的水下波达方向估计方法，其特征在于，所述步骤S4中计算矩阵束{C_hxx,C_hxy}的广义特征值分解过程如下：

S401、X_h(t)的协方差矩阵表示为：

R_hxx＝E[X_h(t)X_h ^H(t)]＝A_xR_ssA_x ^H+σ_x ²I

其中R_ss＝E{S(t)S^H(t)}，为信源部分协方差矩阵；

S402、X_h(t)和Y_h(t)的互协方差矩阵表示为：

R_hxy＝E{X_h(t)Y_h ^H(t)}＝A_xR_ssΦ_x ^HA_x ^H+σ_x ²Z；

S403、对上述协方差矩阵以及互协方差矩阵进行特征值分解得到最小特征值为σ_x ²，利用σ_x ²可以得到矩阵束{C_hxx,C_hxy}，其中C_hxy＝R_hxy-σ_x ²Z＝A_xR_ssΦ_x ^HA_x ^H；

S404、计算所述矩阵束{C_hxx,C_hxy}的广义特征值分解，得到非零特征值λ_x1,λ_x2,…,λ_xK,因此由公式(2)可记：

其中φ_xi为矩阵Φ_x上的对角元素，且φ_xi∈{λ_x1,λ_x2,…,λ_xK},i＝1,2,…,K；

所述步骤S4中计算矩阵束{C_vxx,C_vxy}的广义特征值分解过程如下：

S4X1、X_v(t)的协方差矩阵表示为：

R_vxx＝E[X_v(t)X_v ^H(t)]＝A_yR_ssA_y ^H+σ_y ²I

其中R_ss＝E{S(t)S^H(t)}，为信源部分协方差矩阵；

S4X2、X_v(t)和Y_v(t)的互协方差矩阵表示为：

R_vxy＝E{X_v(t)Y_v ^H(t)}＝A_yR_ssΦ_y ^HA_y ^H+σ_y ²Z；

S4X3、对上述协方差矩阵以及互协方差矩阵进行特征值分解得到最小特征值σ_y ²，利用σ_y ²得到矩阵束其中C_vxy＝R_vxy-σ_y ²Z＝A_yR_ssΦ_y ^HA_y ^H；

S4X4、计算所述矩阵束{C_vxx,C_vxy}的广义特征值分解，得到非零特征值λ_y1,λ_y2,…,λ_yK,因此由公式(3)可记：

其中φ_yi为矩阵Φ_y上的对角元素，且φ_yi∈{λ_y1,λ_y2,…,λ_yK},i＝1,2,…,K。

6.根据权利要求5所述的一种基于正交均匀线阵的水下波达方向估计方法，其特征在于，所述步骤S5具体为：

S51、分别以所述水平均匀线阵以及所述竖直均匀线阵所在的方向建立一个直角坐标系，其中，θ_v1和θ_h1分别为第一象限中声波与两个均匀线阵法线的夹角，θ_v2和θ_h2分别为第二象限中声波与两个均匀线阵法线的夹角；

S52、当声波从第一象限中入射时，此时有θ_xi＝θ_h1，θ_yi＝θ_v1，则θ_xi和θ_yi的关系为：

S53、当声波从第二象限中入射时，此时有θ_xi＝-θ_h2，θ_yi＝θ_v2，则θ_xi和θ_yi的关系为：

S54、根据公式(6)和公式(7)得到：

7.根据权利要求6所述的一种基于正交均匀线阵的水下波达方向估计方法，其特征在于，所述步骤S6具体为：

S61、将公式(8)带入公式(5)，则有：

S62、根据公式(4)和公式(9)可知，若配对成功，则有以下式子成立：

其中，λ_i是第i个目标源反射回来的声波波长，且有λ_i＝v_i/f；

S63、将arg(λ_x1),arg(λ_x2),…,arg(λ_xK)按照各自的平方大小顺序从大到小排列得到序列Η；将arg(λ_y1),arg(λ_y2),…,arg(λ_yK)按照各自的平方大小顺序从小到大排列得到序列V，于是有：

其中h_i为序列Η中的第i个元素；v_i为序列V中的第i个元素。

8.根据权利要求7所述的一种基于正交均匀线阵的水下波达方向估计方法，其特征在于，所述θ_xi的表达式；

当arg(φ_xi)大于0时，θ_xi取正值，即目标声源在x轴负方向；当arg(φ_x)_i小于0时，θ_xi取负值，即目标声源在x轴正方向。

9.一种基于正交均匀线阵的水下波达方向估计***，其特征在于，所述估计***包括数据采集处理与控制模块、发射模块、接收模块、输出模块和电源模块，

其中，所述数据采集处理与控制模块包括处理器、D/A转换器和一个A/D转换器，并且分别与所述发射模块、所述接收模块、所述输出模块和所述电源模块相连，能够通过设定参数来控制所述发射模块，使所述发射模块发射指定的声波信号；还能够对所述接收模块传送回来的数据进行采集与处理，计算出精确的波达方向角；该模块还能够将处理完毕的数据传送给所述输出模块，然后通过USB接口输出到外部装置或者在显示器上显示出来；

所述发射模块包括至少一个超声波发送探头，并且与所述数据采集处理与控制模块和所述电源模块相连，能够根据所述数据采集处理与控制模块的指令，获取指定的信号并对其进行超声发射；

所述接收模块包括两个以相互垂直方式放置的超声波接收探头阵列，并与所述数据采集处理与控制模块和所述电源模块相连，能够接收目标声源发射回来的信号，并将其传送给所述数据采集处理与控制模块；

所述输出模块包括至少一个USB接口和至少一个显示器组成，并且与所述数据采集处理与控制模块模块和所述电源模块相连，能够提供人机交互，将所述数据采集处理与控制模块模块中处理好的数据通过USB接口输出到外部装置或者在显示器上显示出来；

所述电源模块与上述各模块相连，为上述各模块供电。

10.根据权利要求9所述的一种基于正交均匀线阵的水下波达方向估计***，其特征在于，所述数据采集处理与控制模块中处理器为TMS320VC5509A型号的DSP芯片，所述接收模块中超声波接收探头阵列为正交的均匀直线阵列，其中每个阵列包括多个超声接收探头，并且数量相同。