CN106500820A - 一种用于二维波达方向估计的声速测量方法及装置 - Google Patents

Abstract

本发明公开了一种用于二维波达方向估计的声速测量方法及装置。测量方法包括如下步骤：确定水平均匀线阵和竖直均匀线阵的接收信号模型；计算水平均匀线阵和竖直均匀线阵的旋转算子；建立声波在正交均匀线阵上两个方向角之间的关系；根据两方向角之间的关系求出探测路径上声速的大小。装置包括数据采集处理与控制模块、发射模块、接收模块、输出模块和电源模块；数据采集处理与控制模块由一个处理器、一个D/A转换器和一个A/D转换器组成，并且与发射模块、接收模块、输出模块和电源模块相连；本发明的一种用于二维DOA估计的声速测量方法，实现对二维DOA估计中的探测路径声速大小进行精确的测量，本发明装置可行性强，安装简单。

Description

一种用于二维波达方向估计的声速测量方法及装置

技术领域

本发明涉及声速测量的技术领域，尤其是指一种用于二维波达方向估计的声速测量方法及装置。

背景技术

空间信号的波达方向(简称：DOA)估计是阵列信号处理领域的一项重要技术，它利用阵列对空间信号进行接收和处理，提取有用的信号特征或信息，以此来估计信号的入射角，并常常用于水下目标定位中。二维DOA估计可以对目标的方向角和距离进行估计，但最终的估计结果受探测路径上的声速大小影响很大，如果探测路径声速测量得不精准，就会对二维DOA估计的结果产生影响。而在进行水下目标定位时，二维DOA估计通常是利用一片海水区域内的声速来代替每条探测路径上的声速，但是海水环境复杂，且受水流因素影响较大，每一条探测路径上的声速其实是不一样的，即二维DOA估计中对探测路径声速的确定是不准确的，这将导致二维DOA估计的估计误差较大。因此精准的水下探测路径声速测量对二维DOA估计至关重要。

目前，海水声速测量方法主要分为两种：一种是直接法，另外一种是间接法。直接法是用水声设备直接在海洋现场测量海水声速，它可分为时差法、共振声谱法、驻波干涉法与相位比较法等方法。而间接法则是根据海水中的深度、温度和盐度三个主要影响海水声速大小的因素，经过成百上千次的分析和实践，提出一种经验公式，利用此经验公式来测量海水声速。根据经验公式的不同，间接法可分为Del Grosso声速算法、Wilson声速算法和Chen-Millero-Li声速算法三种声速算法。直接法和间接法虽然原理简单且操作较方便，但它们都需要在目标位置已知的情况下才能对某一条探测路径上的声速进行测量，而在二维DOA估计中，目标位置都未知的且是用来估计的；同时，间接法测量的是一个声速分布较均匀的海水环境中的声速，它测量的并不是某一条路径上的声速，即点对点声速，因此直接法和间接法对二维DOA估计的水下探测路径声速大小的测量并不具有实用性。

为了克服在利用直接法和间接法测量二维DOA估计的声速过程中存在的局限性，本专利使用两个相互正交的均匀阵列作为接收阵列，利用从目标反射回来的声波在两个阵列中形成的方向角间的关系，得到二维DOA估计中探测路径声速与两个旋转算子之间的关系式，然后再计算出声速的大小，实现了在二维DOA估计中对水下探测路径声速的精准测量。因此，本专利方法对二维DOA估计中的声速测量更具有实用性。

发明内容

本发明的目的在于克服现有方法存在的各种不足，提供一种用于二维DOA估计的声速测量方法，实现对二维DOA估计中的探测路径声速大小进行精确的测量。

另外，本发明另一目的在于提供一种用于二维DOA估计的声速测量装置。

本发明方法的原理是基于从目标反射回来的声波在两个正交的均匀线阵上形成的方向角相加或相减等于的关系，得出两个均匀线阵上的旋转算子与二维DOA估计中的探测路径声速间的关系式，最后通过计算两均匀线阵上的旋转算子来求得最终探测路径声速大小的。

为了实现上述目的，本发明至少通过以下技术方案实现。

一种用于二维波达方向估计的声速测量方法，其包括如下步骤：

步骤一：确定水平均匀线阵和竖直均匀线阵的接收信号模型；

步骤二：计算水平均匀线阵和竖直均匀线阵的旋转算子；

步骤三：建立声波在正交均匀线阵上两个方向角之间的关系；

步骤四：根据两方向角之间的关系求出探测路径上声速v的大小。

进一步地，步骤一具体包括：

正交均匀线阵的水平和竖直均匀线阵都有M个接收阵元，M为正整数，且阵元间间距为d，窄带目标声源为S，中心频率为f，目标对应于水平线阵的方向角为θ_x，对应于竖直线阵的方向角为θ_y；其中水平均匀线阵的M个阵元为x₁,x₂,…,x_M，竖直均匀线阵的M个阵元为y₁,y₂,…,y_M；

将水平线阵中的M个阵元分为平移矢量为d的两个子阵列X_h和Y_h；子阵列X_h由水平线阵中的第一到第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)分别是子阵列X_h上第一个阵元到第M-1个阵元接收到的信号；

子阵列Y_h由水平线阵的第二到第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)分别是子阵列Y_h上第一个阵元到第M-1个阵元接收到的信号；

x₁(t),x₂(t),…,x_M(t)是水平线阵中第一个阵元到第M个阵元的接收信号；

子阵列X_h接收信号中，以第一个阵元x_h1为参考点，则第一个阵元接收的信号为：

x_h1(t)＝s(t)+n_hx1(t)

其中s(t)表示目标信号，n_hx1(t)表示子阵列X_h第一个阵元上的噪声；

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

x_hm(t)＝s(t)a_m(θ_x)+n_hxm(t),m＝1,2,…,M-1

其中a_m(θ_x)中v表示声波在探测路径上的速度，n_hxm(t)表示子阵列X_h上第m个阵元上的噪声；

由于子阵列Y_h和子阵列X_h的相对平移矢量为d，那么子阵列Y_h上第m个阵元在同一时刻接收到的信号为：

其中n_hym(t)表示子阵列X_h上第m个阵元上的噪声；

将子阵列X_h和子阵列Y_h上的各阵元的接收信号排列成列向量形式，则上面两式可写成矢量形式：

X_h(t)＝A_xs(t)+N_hx(t) (公式1)

Y_h(t)＝A_xΦ_xs(t)+N_hy(t) (公式2)

其中，为(M-1)×1的子阵列X_h导向矢量矩阵，

X_h(t)＝[x_h1(t),x_h2(t),…,x_h(M-1)(t)]^T为(M-1)×1的子阵列X_h接收信号矩阵，

Y_h(t)＝[y_h1(t),y_h2(t),…,y_h(M-1)(t)]^T为(M-1)×1的子阵列Y_h接收信号矩阵，

N_hx(t)＝[n_hx1(t),n_hx2(t),…,n_hx(M-1)(t)]^T为(M-1)×1的子阵列X_h噪声矩阵，

N_hy(t)＝[n_hy1(t),n_hy2(t),…,n_hy(M-1)(t)]^T为(M-1)×1的子阵列Y_h噪声矩阵；

而Φ_x是把子阵X_h和Y_h的输出联系起来的一个因子，称旋转算子，它包含了目标反射回来的信号波前在任意水平均匀线阵中一个阵元偶之间的相位延迟信息，表示为：

对于竖直均匀线阵，同样将其分为两个平移矢量为d的子阵列X_v和Y_v；子阵列X_v由竖直阵列的第一到第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)分别是子阵列X_v上第一个阵元到第M-1个阵元接收到的信号；

子阵列Y_v由竖直阵列的第二到第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)分别是子阵列Y_v上第一个阵元到第M-1个阵元接收到的信号；

y₁(t),y₂(t),…,y_M(t)是竖直线阵中第一个阵元到第M个阵元的接收信号；

以第一个阵元x_v1为参考阵元，那么子阵列X_v和子阵列Y_v中第m个阵元的接收信号分别为：

x_vm(t)＝s(t)a_m(θ_y)+n_vxm(t),m＝1,2,…,M-1

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

X_v(t)＝A_ys(t)+N_vx(t) (公式4)

Y_v(t)＝A_yΦ_ys(t)+N_vy(t) (公式5)

其中，为(M-1)×1的子阵列X_v导向矢量矩阵，

X_v(t)＝[x_v1(t),x_v2(t),…,x_v(M-1)(t)]^T为(M-1)×1的子阵列X_v接收信号矩阵，

Y_v(t)＝[y_v1(t),y_v2(t),…,y_v(M-1)(t)]^T为(M-1)×1的子阵列Y_v接收信号矩阵，

N_vx(t)＝[n_vx1(t),n_vx2(t),…,n_vx(M-1)(t)]^T为(M-1)×1的子阵列X_v噪声矩阵，

N_vy(t)＝[n_vy1(t),n_vy2(t),…,n_vy(M-1)(t)]^T为(M-1)×1的子阵列Y_v噪声矩阵；

同样，Φ_y是把子阵X_v和Y_v的输出联系起来的一个因子，也称旋转算子，它包含了目标反射回来的信号波前在任意竖直均匀线阵中一个阵元偶之间的相位延迟信息，表示为：

进一步地，步骤二具体包括：

对于水平均匀线阵，计算X_h(t)的协方差矩阵R_hxx与X_h(t)和Y_h(t)的互协方差矩阵R_hxy，然后对其进行特征值分解得到最小特征值为σ_x ²，利用σ_x ²计算C_hxx＝R_hxx-σ_x ²I和C_hxy＝R_hxy-σ_x ²Z；其中I是(M-1)×(M-1)阶的单位矩阵，Z也是(M-1)×(M-1)阶的矩阵，并且

最后计算矩阵束{C_hxx,C_hxy}的广义特征值分解，得到非零特征值λ_x,此特征值就是旋转算子的值，即：

对于竖直均匀线阵，同样计算X_v(t)的协方差矩阵R_vxx与X_v(t)和Y_v(t)的互协方差矩阵R_vxy，然后对其进行特征值分解得到最小特征值为σ_y ²，利用σ_y ²计算C_vxx＝R_vxx-σ_y ²I和C_vxy＝R_vxy-σ_y ²Z。最后计算矩阵束{C_vxx,C_vxy}的广义特征值分解，得到非零特征值λ_y,此特征值就是旋转算子的值，即：

进一步地，步骤三具体包括：

在水面位置对水下的探测路径声速进行测量，只考虑从水平均匀线阵下方反射的声波，反射回来的声波分为两种，一种是声波①在水平均匀线阵上到达的第一个阵元不是参考阵元；另一种是声波②在水平均匀线阵上到达的第一个阵元是参考阵元，，设θ_x1和θ_y1分别为声波①与两个均匀线阵法线的夹角，它们相加之和为θ_x2和θ_y2分别为声波②与两个均匀线阵法线的夹角，它们的相加之和也为

当声波在水平均匀线阵上到达的第一个阵元不是参考阵元，即属于声波①的情形时，参考阵元x₁是最晚接收到信号的，因此时延参数τ是小于0的，又因为所以此时声波在水平线阵上的方向角是负的，则有θ_x＝-θ_x1，θ_y＝θ_y1，由此可得θ_x和θ_y的关系为：

当声波在水平均匀线阵上到达的第一个阵元是参考阵元，即属于声波②的情形时，有θ_x＝θ_x2，θ_y＝θ_y2，由此可得θ_x和θ_y的关系为：

根据公式(9)和公式(10)可以得到：

进一步地，步骤四具体包括：

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

则联立公式(7)和公式(12)消去θ_x可得：

整理可得速度v的最终表达为：

因此，利用正交均匀线阵接收到的声波信号求出两个旋转算子的值即λ_x和λ_y，然后将其带入公式(14)便可计算出二维DOA估计中探测路径上声速v的大小。

进一步地，当有多个目标即多条探测路径时，两个均匀线阵上的旋转算子是一个对角矩阵，每一个对角元素都对应于一个目标，求得的特征值也有多个，它们一一对应着对角矩阵上的对角元素，最后利用这些特征值求得二维DOA估计中每条探测路径上的声速。

实现前述测量方法的装置，其包括数据采集处理与控制模块、发射模块、接收模块、输出模块和电源模块；数据采集处理与控制模块由一个处理器、一个D/A转换器和一个A/D转换器组成，并且与发射模块、接收模块、输出模块和电源模块相连；处理器能够通过设定参数来控制发射模块，使发射模块发射指定的声波信号；还能够对接收模块传送回来的数据进行采集与处理，然后运计算出二维DOA估计中精准的水下探测路径声速大小；处理器还将处理好的数据传送给输出模块，再通过USB接口输出到外部装置或者在显示器上显示出来。

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

接收模块由两个相互正交的超声波接收探头阵列和阻抗匹配电路组成，与数据采集处理与控制模块和电源模块相连，能够接收从目标发射回来的信号，并传送给数据采集处理与控制模块；

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

本发明与现有技术相比，具有如下优点和技术效果：

1、在原理上，本发明方法利用正交阵列可以对二维DOA估计的水下探测路径声速进行精准的测量，克服了直接法和间接法需要在目标位置已知情况下才能对一条路径上的声速进行测量的局限性，因此本发明算法在二维DOA估计中的水下探测路径声速测量中更具有实用性。

2、在操作复杂度上，本发明的设备较直接法的简单，并且操作方便；同时，本发明的测量方法不受海水环境的影响，无需像间接法一样需要对水的深度、温度和盐度进行测量，因此本发明算法在二维DOA估计中的水下探测路径声速测量中更具方便性。

3、本发明方法可以与传统的二维DOA估计算法相结合，如MUSIC算法和ESPIRT算法，使得这些传统的二维DOA估计算法的估计精度更佳准确。

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

附图说明

图1为实例中装置的硬件结构模块图。

图2为两个相互正交的超声波接收探头阵列模型图。

图3为实例中所用的正交均匀线阵模型。

图4为水平均匀线阵中子阵列X_h的接收信号的示意图。

图5为实例中测量方法的流程图。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

如图5所示，本实例的测量方法具体包括以下几步：

步骤一：推导水平均匀线阵和竖直均匀线阵的接收信号模型。

正交均匀线阵模型场景如图3所示。水平和竖直均匀线阵都有M个接收阵元，且阵元间间距为d，窄带目标声源为S，中心频率为f，目标对应于水平线阵的方向角为θ_x，对应于竖直线阵的方向角为θ_y。其中水平均匀线阵的M个阵元为x₁,x₂,…,x_M，竖直均匀线阵的M个阵元为y₁,y₂,…,y_M。

将水平线阵中的M个阵元分为平移矢量为d的两个子阵列X_h和Y_h。子阵列X_h由水平线阵中的第一到第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)分别是子阵列X_h上第一个阵元到第M-1个阵元接收到的信号。

子阵列Y_h由水平线阵的第二到第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)分别是子阵列Y_h上第一个阵元到第M-1个阵元接收到的信号。

x₁(t),x₂(t),…,x_M(t)是水平线阵中第一个阵元到第M个阵元的接收信号。

子阵列X_h接收信号的示意图如图4所示。以第一个阵元x_h1为参考点，则第一个阵元接收的信号为：

x_h1(t)＝s(t)+n_hx1(t)

其中s(t)表示目标信号，n_hx1(t)表示子阵列X_h第一个阵元上的噪声。

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

x_hm(t)＝s(t)a_m(θ_x)+n_hxm(t),m＝1,2,…,M-1

其中a_m(θ_x)中v表示声波在探测路径上的速度，n_hxm(t)表示子阵列X_h上第m个阵元上的噪声。

由于子阵列Y_h和子阵列X_h的相对平移矢量为d，那么子阵列Y_h上第m个阵元在同一时刻接收到的信号为：

其中n_hym(t)表示子阵列X_h上第m个阵元上的噪声。

将子阵列X_h和子阵列Y_h上的各阵元的接收信号排列成列向量形式，则上面两式可写成矢量形式：

X_h(t)＝A_xs(t)+N_hx(t) (公式1)

Y_h(t)＝A_xΦ_xs(t)+N_hy(t) (公式2)

其中，为(M-1)×1的子阵列X_h导向矢量矩阵，

X_h(t)＝[x_h1(t),x_h2(t),…,x_h(M-1)(t)]^T为(M-1)×1的子阵列X_h接收信号矩阵，

Y_h(t)＝[y_h1(t),y_h2(t),…,y_h(M-1)(t)]^T为(M-1)×1的子阵列Y_h接收信号矩阵，

N_hx(t)＝[n_hx1(t),n_hx2(t),…,n_hx(M-1)(t)]^T为(M-1)×1的子阵列X_h噪声矩阵，

N_hy(t)＝[n_hy1(t),n_hy2(t),…,n_hy(M-1)(t)]^T为(M-1)×1的子阵列Y_h噪声矩阵。

而Φ_x是把子阵X_h和Y_h的输出联系起来的一个因子，称旋转算子，它包含了目标反射回来的信号波前在任意水平均匀线阵中一个阵元偶之间的相位延迟信息，表示为：

对于竖直均匀线阵，同样将其分为两个平移矢量为d的子阵列X_v和Y_v。子阵列X_v由竖直阵列的第一到第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)分别是子阵列X_v上第一个阵元到第M-1个阵元接收到的信号。

子阵列Y_v由竖直阵列的第二到第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)分别是子阵列Y_v上第一个阵元到第M-1个阵元接收到的信号。

y₁(t),y₂(t),…,y_M(t)是竖直线阵中第一个阵元到第M个阵元的接收信号。

以第一个阵元x_v1为参考阵元，那么子阵列X_v和子阵列Y_v中第m个阵元的接收信号分别为：

x_vm(t)＝s(t)a_m(θ_y)+n_vxm(t),m＝1,2,…,M-1

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

X_v(t)＝A_ys(t)+N_vx(t) (公式4)

Y_v(t)＝A_yΦ_ys(t)+N_vy(t) (公式5)

其中，为(M-1)×1的子阵列X_v导向矢量矩阵，

X_v(t)＝[x_v1(t),x_v2(t),…,x_v(M-1)(t)]^T为(M-1)×1的子阵列X_v接收信号矩阵，

Y_v(t)＝[y_v1(t),y_v2(t),…,y_v(M-1)(t)]^T为(M-1)×1的子阵列Y_v接收信号矩阵，

N_vx(t)＝[n_vx1(t),n_vx2(t),…,n_vx(M-1)(t)]^T为(M-1)×1的子阵列X_v噪声矩阵，

N_vy(t)＝[n_vy1(t),n_vy2(t),…,n_vy(M-1)(t)]^T为(M-1)×1的子阵列Y_v噪声矩阵。

同样，Φ_y是把子阵X_v和Y_v的输出联系起来的一个因子，也称旋转算子，它包含了目标反射回来的信号波前在任意竖直均匀线阵中一个阵元偶之间的相位延迟信息，表示为：

步骤二：计算水平均匀线阵和竖直均匀线阵的旋转算子。

对于水平均匀线阵，计算X_h(t)的协方差矩阵R_hxx与X_h(t)和Y_h(t)的互协方差矩阵R_hxy，然后对其进行特征值分解得到最小特征值为σ_x ²，利用σ_x ²计算C_hxx＝R_hxx-σ_x ²I和C_hxy＝R_hxy-σ_x ²Z。其中I是(M-1)×(M-1)阶的单位矩阵，Z也是(M-1)×(M-1)阶的矩阵，并且

最后计算矩阵束{C_hxx,C_hxy}的广义特征值分解，得到非零特征值λ_x,此特征值就是旋转算子的值，即：

对于竖直均匀线阵，同样计算X_v(t)的协方差矩阵R_vxx与X_v(t)和Y_v(t)的互协方差矩阵R_vxy，然后对其进行特征值分解得到最小特征值为σ_y ²，利用σ_y ²计算C_vxx＝R_vxx-σ_y ²I和C_vxy＝R_vxy-σ_y ²Z。最后计算矩阵束{C_vxx,C_vxy}的广义特征值分解，得到非零特征值λ_y,此特征值就是旋转算子的值，即：

步骤三：建立声波在正交均匀线阵上两个方向角之间的关系。

由于是在水面位置对水下的探测路径声速进行测量，因此可只考虑从水平均匀线阵下方反射的声波。反射回来的声波分为两种，一种是声波在水平均匀线阵上到达的第一个阵元不是参考阵元，如图2中的声波①；另一种是声波在水平均匀线阵上到达的第一个阵元是参考阵元，如图2中的声波②。设θ_x1和θ_y1分别为声波①与两个均匀线阵法线的夹角，它们相加之和为θ_x2和θ_y2分别为声波②与两个均匀线阵法线的夹角，它们的相加之和也为

当声波在水平均匀线阵上到达的第一个阵元不是参考阵元，即属于声波①的情形时，参考阵元x₁是最晚接收到信号的，因此时延参数τ是小于0的，又因为所以此时声波在水平线阵上的方向角是负的，则有θ_x＝-θ_x1，θ_y＝θ_y1，由此可得θ_x和θ_y的关系为：

当声波在水平均匀线阵上到达的第一个阵元是参考阵元，即属于声波②的情形时，有θ_x＝θ_x2，θ_y＝θ_y2，由此可得θ_x和θ_y的关系为：

根据公式(9)和公式(10)可以得到：

步骤四：根据两方向角之间的关系求出探测路径上声速v的大小。

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

则联立公式(7)和公式(12)消去θ_x可得：

整理可得速度v的最终表达为：

因此，利用正交均匀线阵接收到的声波信号求出两个旋转算子的值(即λ_x和λ_y)，然后将其带入公式(14)便可计算出二维DOA估计中探测路径上声速v的大小。

根据以上算法流程可知，本发明提出的声速测量方法可以对二维DOA估计中的水下探测路径声速进行精确的测量，克服了现有直接测量法和间接测量法存在的缺陷。

同时，本发明方法可以推广到有多个目标的情况中去，这时候的两个均匀线阵上的旋转算子是一个对角矩阵，每一个对角元素都对应于一个目标，求得的特征值也有多个，它们一一对应着对角矩阵上的对角元素，最后利用这些特征值可以求得二维DOA估计中每条探测路径上的声速。

本实例用于二维DOA估计的声速测量装置如图1所示，装置包括数据采集处理与控制模块、发射模块、接收模块、输出模块和电源模块。数据采集处理与控制模块可以用DSP芯片实现(如：TI公司的TMS320VC5509A型号的DSP芯片)，此DSP芯片可实现A/D转换和D/A转换的功能，并能够实现正交均匀线阵的旋转算子和最终声波传播路径速度的计算；发射模块使用一个超声波发射探头，利用阻抗匹配电路与D/A连接；接收模块使用两个正交的均匀线阵，其中每个阵列包括多个超声接收探头，并且数量相同，，采用阻抗匹配电路与A/D连接；输出模块使用一个USB接口和一个LCD显示屏。图1即为本发明所述装置的硬件结构模块图。

本装置的主要工作流程如下：操作者通过数据采集处理与控制模块输入对应的参数，使该模块产生所需要的数字信号，该模块中的D/A转换器将其转换成模拟信号后，传送给发发射模块，超声波发射探头将对其进行发射。接收模块中的正交均匀线阵接收到从目标源发射回来的信号后将其传送给数据采集处理与控制模块，信号经过该模块中的A/D转换器后变为数字信号，然后模块中的处理器利用本发明算法对其进行计算，得出二维DOA估计中水下探测路径声速的大小。最后数据采集处理与控制模块将计算得出的速度值传给输出模块，使其通过USB接口输出给外部装置或者在显示屏上显示出来。工作的全程中，电源模块为所有模块供电。

作为举例，本实例的装置主要工作步骤具体如下：

步骤1：在一个淡水库中放置一个目标S。设定阵列接收模块和发射模块中的各项参数：将接收模块中的每个均匀线阵中的阵元个数M定为9，则每个均匀线阵中前8个阵元为子阵X_h和X_v，后8个阵元为子阵Y_h和Y_v；两个子阵列之间的距离d设为5mm；通过数据采集处理与控制模块控制发射模块使待发射的超声频率为100kHz，脉冲长度5ms。

步骤2：对超声接收探头接收到的目标声源信号进行采样；水平方向均匀线阵接收到的信号为x₁(t),x₂(t),…,x₉(t)，竖直方向均匀线阵接收的信号为y₁(t),y₂(t),…,y₉(t)。共采样接收200次，并将接收到的信号传递给数据采集处理与控制模块模块做运算处理。

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

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

2)对水平线阵中的两个协方差矩阵R_hxx和R_hxy进行特征值分解，得到最小的特征值，从而有和同时对竖直线阵中的两个协方差矩阵进行相同的处理，得到C_vxx和C_vxy。

3)分别计算矩阵束{C_hxx,C_hxy}和{C_vxx,C_vxy}的广义特征值分解，得到特征值λ_x和λ_y。

4)根据计算得出的λ_x和λ_y，然后联立式子对v进行求解，最终求得

步骤4：将计算出的水下探测路径上的声速信息存储下来，并传送给输出模块，使其通过USB接口输出给外部装置或者显示在LCD显示屏上。根据本发明算法测量出二维DOA估计中的水下探测路径声速大小为1456.4m/s。

Claims (8)

1.一种用于二维波达方向估计的声速测量方法，其特征在于包括如下步骤：

步骤一：确定水平均匀线阵和竖直均匀线阵的接收信号模型；

步骤二：计算水平均匀线阵和竖直均匀线阵的旋转算子；

步骤三：建立声波在正交均匀线阵上两个方向角之间的关系；

步骤四：根据两方向角之间的关系求出探测路径上声速v的大小。

2.根据权利要求1所述的一种用于二维波达方向估计的声速测量方法，其特征在于步骤一具体包括：

正交均匀线阵的水平和竖直均匀线阵都有M个接收阵元，M为正整数，且阵元间间距为d，窄带目标声源为S，中心频率为f，目标对应于水平线阵的方向角为θ_x，对应于竖直线阵的方向角为θ_y；其中水平均匀线阵的M个阵元为x₁,x₂,…,x_M，竖直均匀线阵的M个阵元为y₁,y₂,…,y_M；

将水平线阵中的M个阵元分为平移矢量为d的两个子阵列X_h和Y_h；子阵列X_h由水平线阵中的第一到第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)分别是子阵列X_h上第一个阵元到第M-1个阵元接收到的信号；

子阵列Y_h由水平线阵的第二到第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)分别是子阵列Y_h上第一个阵元到第M-1个阵元接收到的信号；

x₁(t),x₂(t),…,x_M(t)是水平线阵中第一个阵元到第M个阵元的接收信号；

子阵列X_h接收信号中，以第一个阵元x_h1为参考点，则第一个阵元接收的信号为：

x_h1(t)＝s(t)+n_hx1(t)

其中s(t)表示目标信号，n_hx1(t)表示子阵列X_h第一个阵元上的噪声；

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

x_hm(t)＝s(t)a_m(θ_x)+n_hxm(t),m＝1,2,…,M-1

其中a_m(θ_x)中v表示声波在探测路径上的速度，n_hxm(t)表示子阵列X_h上第m个阵元上的噪声；

由于子阵列Y_h和子阵列X_h的相对平移矢量为d，那么子阵列Y_h上第m个阵元在同一时刻接收到的信号为：

$y_{hm}\left(t\right)=s{\left(t\right)}^{j\frac{2\mathrm{\π}f}{v}{\mathrm{dsin\θ}}_x}{a}_m\left({\mathrm{\θ}}_x\right)+n_{hym}\left(t\right),m=1,2,\mathrm{...},M-1$

其中n_hym(t)表示子阵列X_h上第m个阵元上的噪声；

将子阵列X_h和子阵列Y_h上的各阵元的接收信号排列成列向量形式，则上面两式可写成矢量形式：

X_h(t)＝A_xs(t)+N_hx(t) (公式1)

Y_h(t)＝A_xΦ_xs(t)+N_hy(t) (公式2)

其中，为(M-1)×1的子阵列X_h导向矢量矩阵，

X_h(t)＝[x_h1(t),x_h2(t),…,x_h(M-1)(t)]^T为(M-1)×1的子阵列X_h接收信号矩阵，

Y_h(t)＝[y_h1(t),y_h2(t),…,y_h(M-1)(t)]^T为(M-1)×1的子阵列Y_h接收信号矩阵，

N_hx(t)＝[n_hx1(t),n_hx2(t),…,n_hx(M-1)(t)]^T为(M-1)×1的子阵列X_h噪声矩阵，

N_hy(t)＝[n_hy1(t),n_hy2(t),…,n_hy(M-1)(t)]^T为(M-1)×1的子阵列Y_h噪声矩阵；

而Φ_x是把子阵X_h和Y_h的输出联系起来的一个因子，称旋转算子，它包含了目标反射回来的信号波前在任意水平均匀线阵中一个阵元偶之间的相位延迟信息，表示为：

对于竖直均匀线阵，同样将其分为两个平移矢量为d的子阵列X_v和Y_v；子阵列X_v由竖直阵列的第一到第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)分别是子阵列X_v上第一个阵元到第M-1个阵元接收到的信号；

子阵列Y_v由竖直阵列的第二到第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)分别是子阵列Y_v上第一个阵元到第M-1个阵元接收到的信号；

y₁(t),y₂(t),…,y_M(t)是竖直线阵中第一个阵元到第M个阵元的接收信号；

以第一个阵元x_v1为参考阵元，那么子阵列X_v和子阵列Y_v中第m个阵元的接收信号分别为：

x_vm(t)＝s(t)a_m(θ_y)+n_vxm(t),m＝1,2,…,M-1

$y_{vm}\left(t\right)=s{\left(t\right)}^{j\frac{2\mathrm{\π}f}{v}{\mathrm{dsin\θ}}_y}{a}_m\left({\mathrm{\θ}}_y\right)+n_{vym}\left(t\right),m=1,2,\mathrm{...},M-1$

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

X_v(t)＝A_ys(t)+N_vx(t) (公式4)

Y_v(t)＝A_yΦ_ys(t)+N_vy(t) (公式5)

其中，为(M-1)×1的子阵列X_v导向矢量矩阵，

X_v(t)＝[x_v1(t),x_v2(t),…,x_v(M-1)(t)]^T为(M-1)×1的子阵列X_v接收信号矩阵，

Y_v(t)＝[y_v1(t),y_v2(t),…,y_v(M-1)(t)]^T为(M-1)×1的子阵列Y_v接收信号矩阵，

N_vx(t)＝[n_vx1(t),n_vx2(t),…,n_vx(M-1)(t)]^T为(M-1)×1的子阵列X_v噪声矩阵，

N_vy(t)＝[n_vy1(t),n_vy2(t),…,n_vy(M-1)(t)]^T为(M-1)×1的子阵列Y_v噪声矩阵；

同样，Φ_y是把子阵X_v和Y_v的输出联系起来的一个因子，也称旋转算子，它包含了目标反射回来的信号波前在任意竖直均匀线阵中一个阵元偶之间的相位延迟信息，表示为：

3.根据权利要求2所述的一种用于二维波达方向估计的声速测量方法，其特征在于步骤二具体包括：

对于水平均匀线阵，计算X_h(t)的协方差矩阵R_hxx与X_h(t)和Y_h(t)的互协方差矩阵R_hxy，然后对其进行特征值分解得到最小特征值为σ_x ²，利用σ_x ²计算C_hxx＝R_hxx-σ_x ²I和C_hxy＝R_hxy-σ_x ²Z；其中I是(M-1)×(M-1)阶的单位矩阵，Z也是(M-1)×(M-1)阶的矩阵，并且

最后计算矩阵束{C_hxx,C_hxy}的广义特征值分解，得到非零特征值λ_x,此特征值就是旋转算子的值，即：

对于竖直均匀线阵，同样计算X_v(t)的协方差矩阵R_vxx与X_v(t)和Y_v(t)的互协方差矩阵R_vxy，然后对其进行特征值分解得到最小特征值为σ_y ²，利用σ_y ²计算C_vxx＝R_vxx-σ_y ²I和C_vxy＝R_vxy-σ_y ²Z。最后计算矩阵束{C_vxx,C_vxy}的广义特征值分解，得到非零特征值λ_y,此特征值就是旋转算子的值，即：

4.根据权利要求3所述的一种用于二维波达方向估计的声速测量方法，其特征在于步骤三具体包括：

在水面位置对水下的探测路径声速进行测量，只考虑从水平均匀线阵下方反射的声波，反射回来的声波分为两种，一种是声波①在水平均匀线阵上到达的第一个阵元不是参考阵元；另一种是声波②在水平均匀线阵上到达的第一个阵元是参考阵元，，设θ_x1和θ_y1分别为声波①与两个均匀线阵法线的夹角，它们相加之和为θ_x2和θ_y2分别为声波②与两个均匀线阵法线的夹角，它们的相加之和也为

当声波在水平均匀线阵上到达的第一个阵元不是参考阵元，即属于声波①的情形时，参考阵元x₁是最晚接收到信号的，因此时延参数τ是小于0的，又因为所以此时声波在水平线阵上的方向角是负的，则有θ_x＝-θ_x1，θ_y＝θ_y1，由此可得θ_x和θ_y的关系为：

当声波在水平均匀线阵上到达的第一个阵元是参考阵元，即属于声波②的情形时，有θ_x＝θ_x2，θ_y＝θ_y2，由此可得θ_x和θ_y的关系为：

根据公式(9)和公式(10)可以得到：

5.根据权利要求4所述的一种用于二维波达方向估计的声速测量方法，其特征在于步骤四具体包括：

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

则联立公式(7)和公式(12)消去θ_x可得：

整理可得速度v的最终表达为：

因此，利用正交均匀线阵接收到的声波信号求出两个旋转算子的值即λ_x和λ_y，然后将其带入公式(14)便可计算出二维DOA估计中探测路径上声速v的大小。

6.根据权利要求5所述的一种用于二维波达方向估计的声速测量方法，其特征在于：当有多个目标即多条探测路径时，两个均匀线阵上的旋转算子是一个对角矩阵，每一个对角元素都对应于一个目标，求得的特征值也有多个，它们一一对应着对角矩阵上的对角元素，最后利用这些特征值求得二维DOA估计中每条探测路径上的声速。

7.实现权利要求1～6任一项所述测量方法的装置，其特征在于包括数据采集处理与控制模块、发射模块、接收模块、输出模块和电源模块；数据采集处理与控制模块由一个处理器、一个D/A转换器和一个A/D转换器组成，并且与发射模块、接收模块、输出模块和电源模块相连；处理器能够通过设定参数来控制发射模块，使发射模块发射指定的声波信号；还能够对接收模块传送回来的数据进行采集与处理，然后运计算出二维DOA估计中精准的水下探测路径声速大小；处理器还将处理好的数据传送给输出模块，再通过USB接口输出到外部装置或者在显示器上显示出来。

8.根据权利要求7所述的装置，其特征在于发射模块由一个阻抗匹配电路和一个超声波发送探头组成，并且与数据采集处理与控制模块和电源模块相连，能够根据数据采集处理与控制模块的指令，获取指定的信号并进行超声发射；

接收模块由两个相互正交的超声波接收探头阵列和阻抗匹配电路组成，与数据采集处理与控制模块和电源模块相连，能够接收从目标发射回来的信号，并传送给数据采集处理与控制模块；

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