CN113093195A - 基于时频分布相关的水下瞬态声源深度距离估计方法 - Google Patents

Abstract

本发明属于信号与信息处理技术领域，公开了基于时频分布相关的水下瞬态声源深度距离估计方法。利用简正波在距离上的频散特性和深度上模态振幅的分布，通过对接收信号的时频表示提取声源的深度和距离信息，最后对实际接收信号的时频分布与模拟接收信号的时频分布图作相关实现对声源的深度距离估计。对时频分布图作相关时，利用向量相关方法间接相关时频分布图对目标定位。利用向量化时频矩阵向量的信号子空间代替向量化时频矩阵向量与接收信号向量化时频矩阵的噪声子空间作相关对目标定位。直接求取实际接收信号和模拟接收信号原始时频矩阵的协方差矩阵，分别计算其噪声子空间和信号子空间，对以上两种子空间作相关实现声源深度距离估计。

Description

基于时频分布相关的水下瞬态声源深度距离估计方法

技术领域

本发明属于信号与信息处理技术领域，特别涉及到一种利用水听器接收瞬态信号时频分布估计水下目标深度和距离的方法，适用于声纳信号处理。

背景技术

水下目标深度距离估计，一直是水声学领域的一个关键问题。根据水声传播的简正波理论，声源在水下会激发多阶模态，接收点处的声压由各阶模态叠加而成，而声源深度信息和距离信息就包含于接收信号之中。同时，由于简正波为频散波导，各阶简正波在传播过程中会产生模间频散和模内频散，将接收信号表示在时频分布图上可观察到这种现象。频散程度与声源和接收端的相对距离有关，各阶模态在接收端的振幅与声源激发的模态幅度直接相关。

关于时频分布技术可参考文献《Source depth estimation using modaldecomposition and time-frequency representations》，其中有关于简正波理论与时频分析的相关理论。

假定声源位于水下某距离和深度处，则此位置激发的模态到达接收端，所呈现出的时频分布结果，与实际接收信号的时频分布结果相关性强。而其它距离和深度的位置，如果有声源，到达接收端所呈现出的时频分布结果，与实际接收信号的时频分布结果相关性弱。因此，可以以此作为依据，选择合适的算法，以相关性强弱实现水下目标深度距离估计。

发明内容

本发明为了能同时确定水下声源的深度和距离，提供了一种基于时频分布相关的水下瞬态声源深度距离估计方法。利用简正波在距离上的频散特性和深度上模态振幅的分布，通过对接收信号的时频表示提取声源的深度和距离信息，最后对实际接收信号的时频分布与模拟接收信号的时频分布图作相关实现对声源的深度距离估计。

水下瞬态声源在水下激发出不同阶次的模态，各阶模态的振幅与声源的深度有关，声波以各阶模态所对应的波向远处传播。同时，经声场传播之后，各阶模态具有不同的速度，导致在接收端会出现频散现象，从低到高的各阶模态传输速度逐渐变慢，导致第一阶模态最先到达，其它各阶模态传输速度相对变慢，依次到达。随着收发两端相对距离的增大，模态之间的频散现象越来越明显。

利用接收瞬态信号的时频分布，可以观察到各阶模态的分布情况，时频分布图中，模态的振幅和模态的频散特性对于不同深度和距离的水下目标，呈现的时频分布结果将有差异。可以以时频分布的相关性，作为水下目标的定位标准。

本发明的技术方案是：

基于时频分布相关的水下瞬态声源深度距离估计方法，包括步骤如下：

步骤1：根据海洋环境模型对声源目标区域进行网格化处理；

步骤2：计算接收信号；

根据简正波理论，在深度z_s处宽带声源发射一个低频宽带信号，频率带宽为f₁～f₂，经过海洋波导传播至距离r和深度z处对应的接收信号x(r,z,t)为：

其中，j为虚数，r为接收点与声源的水平距离，z为接收点深度，z_s为声源深度，ξ_n(ω)为频率为ω的声源波数k₀的水平分量。Z_n(z,ω)和Z_n(z_s,ω)是接收点和发射点处所对应的第n阶模态振幅，ω_min为最小声源频率，ω_max为最大声源频率，t为时间序列，S(ω)为发生信号频谱。

观察式(1)可发现声源深度信息包含于模态振幅Z_n(z_s,ω)中。具体表现为：对于不同声源深度，接收端处接收到的各阶模态之间比例倍数不同，或者其对应的能量比例关系不同。

又因为简正波为频散波导，各阶模态传播速度不同。第n阶模态信号传播至距离r处的到达时间t_n(ω)满足以下关系：

t_n(ω)＝r/c_gn(ω) (2)

其中c_gn(ω)为频率为ω的第n阶模态的群速度。

由式(2)可得：

其中ΔT_mn(ω)表示第n阶模态和第m阶模态到达接收点的时间差，上式中c_gm(ω)是第m阶模态的群速度。

式(3)决定了各阶模态到达接收端的时间差，也决定了时频分布上的频散程度。对于不同距离的声源，模态之间的频散程度不同。

对式(1)的时域信号作短时傅里叶变换，进行时频分析：

其中，h(t)为短时傅立叶变换窗函数，τ表示在原信号t时刻附近τ时段的信号。

步骤3：对接收信号进行时频分析；

将步骤2得到的模拟接收时域信号和实际接收时域信号分别进行时频分析，得到实际接收信号归一化时频矩阵X_r和网格点处假定声源的模拟接收信号归一化时频矩阵X_c；

在实际过程中由于发射信号频谱未知，为了方便在时频矩阵上作相关处理，因此需要对式(4)中每个频率ω处对应的信号分别进行归一化处理：

其中X(r,z,ω_i,t)表示第i个ω对应的信号向量，X_j(r,z,ω_i,t)表示第i个ω对应的信号向量中的第j个元素，

表示第i个ω对应的信号向量中的第j个归一化后的元素。

将搜索区域进行网格划分，对不同距离和不同深度上的每个网格点，计算到达接收端的时频分布。同时，对接收水听器实际接收信号x_r(r,z,t)进行短时傅立叶变换，实际接收信号的时频分布为X_r(r,z,ω,t)。

步骤4：声源深度距离估计；

利用步骤3得到的实际接收信号的时频矩阵X_r与每个网格点处假定声源信号的时频矩阵X_c，同时计算出实际接收信号和模拟接收信号关于声源深度的各阶模态能量比以及关于声源距离的各阶模态时间差；

在信号的时频分布图中包含声源深度和距离信息，从理论上讲，深度信息包含在模态振幅中，不同声源深度处的各阶模态的振幅所占比重不同，图1和图2给出了两个声源深度处的模态振幅和接收信号时频分布图；而距离信息是根据简正波的频散特性，随着距离的增大，各阶模态都会产生频散现象，通过计算实际接收信号的各阶模态之间到达接收点的时间差，然后再将其与模拟接收信号的各阶模态之间达到接收点时间差作相关就可估计声源距离，图3给出了在频率为300Hz情况下不同距离的各阶模态的时间差示意图。

利用实际接收信号的时频分布矩阵，与距离r和深度z网格点处假定声源信号的时频分布矩阵相匹配进行声源定位。可采用的方法如下：

方法一：

式中x_rnm为实际接收信号时频矩阵X_r中的第n行、第m列的元素，x_cnm(r,z)为在实测海洋环境条件下，以距离r和深度z仿真计算所得的数据时频矩阵X_c中的第n行、第m列的元素，N和M分别为时-频矩阵的行数和列数。通过式(6)的代价函数，在网格化区域内进行峰值搜索确定声源的位置。

方法二：

将N行、M列的时频矩阵X_N×M进行列向量化：

K＝[k₁,k₂,…,k_M] (7)

其中k₁、k₂、k_M分别为时-频矩阵X_N×M的第1列、第2列和第M列元素，它们是N×1维的向量。

然后对矩阵K的M个N×1维的向量求转置再进行矩阵合并：

其中矩阵U中的元素u₁,u₂,u_N×M分别对应时频矩阵X_N×M中的元素x₁₁,x₂₁,x_NM。

因为对接收信号进行时频分析得到的矩阵元素为复数，我们需要将其取模然后进行之后的处理：

U＝[|u₁|,|u₂|,…,|u_N×M|] (9)

对向量化的矩阵作处理，计算实际接收信号向量化时频矩阵的协方差矩阵：

其中U_r表示实际接收信号时频矩阵向量化后的1×(N×M)维的时频矩阵。

以距离r和深度z仿真计算所得的时频矩阵数据向量化后向量为U_c(r,z)，通过下式可估计水下目标。

方法三：

该方法对噪声有一定程度上的抑制能力。

方法四：

基于矩阵特征分解的相关理论，利用噪声子空间与信号子空间的相互正交这一关系将噪声子空间与网格点处的模拟接收信号作相关，在网格化区域内进行峰值搜索，确定出声源的位置。

其中N_r＝[q_k,q_k+1,…,q_M×N]为噪声子空间，q_k为对协方差矩阵R_r进行特征分解得到的第k个特征值对应的特征向量。

方法五：

式(13)中的U_c(r,z)是通过(N×M)维的时频矩阵向量化得到的1×(N×M)维时频矩阵，求模拟接收信号向量化时频矩阵的协方差矩阵：

对式(14)进行矩阵特征分解得到信号子空间：

S_c(r,z)＝[q₁,q₂,…,q_k-1] (15)

其中S_c(r,z)中的元素为对协方差矩阵R_c(r,z)进行特征分解得到的第1到第k-1个特征值对应的特征向量。

使用式(15)中的S_c(r,z)代替方法四中U_c(r,z)与式(13)中的N_r作相关得到目标函数：

方法六：

首先求(N×M)维的实际接收信号时频矩阵的协方差矩阵：

利用矩阵特征分解得到噪声子空间：

N_xr＝[q_k,q_k+1,…,q_M] (18)

其中q_k为对协方差矩阵R_xr进行特征分解得到的第k个特征值对应的特征向量。

对模拟接收信号的(N×M)时频矩阵采用相同的处理方式得到信号子空间：

S_xc＝[q₁,q₂,…,q_k-1] (19)

其中S_xc中的元素为对协方差矩阵R_xr进行特征分解得到的第1到第k-1个特征值对应的特征向量。

使用式(19)计算网格点处假定声源的接收信号的(N×M)维时频矩阵的信号子空间，然后与式(18)作相关在网格化区域内峰值搜索确定声源位置。

本发明的有益效果为，利用包含声源深度距离信息的时频分布图作相关，时频分布图可以清楚展现出假定声源的深度和距离信息。对时频分布图作相关时的处理手段，对表示时频分布图的时频矩阵进行列向量化，利用向量相关方法一～四间接相关时频分布图对目标定位。利用向量化时频矩阵向量的信号子空间代替向量化时频矩阵向量与接收信号向量化时频矩阵的噪声子空间作相关，对应于声源深度距离估计方法五。直接求取实际接收信号和模拟接收信号原始时频矩阵的协方差矩阵，分别计算其噪声子空间和信号子空间，对以上两种子空间作相关实现声源深度距离估计；对应于声源深度距离估计方法六。

附图说明

图1模态振幅随声源深度分布的示意图。

图2是不同深度处的声源在接收点处的时域信号时频分布图。(a)声源深度为14.5m处的信号时频图，(b)声源深度为7.5m处的信号时频图。

图3是不同距离处的声源在接收点处的时域信号时频分布图。(a)声源距离为10km处的信号时频图，(b)声源距离为20km处的信号时频图。

图4本发明使用的浅海Pekeris波导环境模型。

图5是方法一声源的定位模糊度平面。

图6是方法二声源的定位模糊度平面。

图7是方法三声源的定位模糊度平面。

图8是方法四声源的定位模糊度平面。

图9是方法五声源的定位模糊度平面。

图10是方法六声源的定位模糊度平面。

具体实施方式

以下结合方案和附图详细叙述本发明的具体实施例子。

步骤1：计算模态振幅和声源波数水平分量。

根据图4所示的海洋环境模型对声源可能存在的区域：深度0～40m，18000～19000m进行网格化处理，深度搜索步长为0.5m，距离搜索步长为10m，然后对每一个网格点处的假定声源和深度20m、距离18500m的实际声源使用Kraken声学工具包计算其模态振幅Z_n和声源波数的水平分量ξ_n。

步骤2：计算接收信号。

在某一点处的声源发射一个低频宽带信号，频率带宽为200～300Hz,中心频率为250Hz，经过海洋波导传播在接收点表示为：

其中，j为虚数，r为接收点与声源的水平距离，z为接收点深度，z_s为声源深度，ξ_n(ω)为频率为ω的声源波数k₀的水平分量。Z_n(z,ω)和Z_n(z_s,ω)是接收点和发射点处所对应的第n阶模态振幅，ω_min为最小声源频率，ω_max为最大声源频率，t为时间序列，S(ω)为发生信号频谱。

将步骤1得到的模态振幅和水平分量带入上式，得到网格点处假定声源和实际声源在接收点处的时域接收信号。

步骤3：对接收信号进行时频分析。

将步骤2得到的模拟接收时域信号和实际接收时域信号分别进行时频分析得到实际接收信号归一化时频矩阵X_r和网格点处假定声源的模拟接收信号归一化时频矩阵X_c。

图2是距离为18000m，深度14.5m和7.5m处的接收信号的时频图。相比于图2(b)，观察发现图2(a)的第三条曲线没有表现出来，而图2(b)的第三条曲线的亮度却最高。这个现象正好对应图1中的模态振幅的分布情况，在图1中当深度为14.5m时，第三阶模态振幅接近于0，当深度为7.5m时，图1中的第三阶模态振幅达到最大值，对应于图2(b)中对应的第三条曲线亮度最高，仔细观察图1的模态振幅大小和图2中曲线的明亮程度，可发现它们两者之间是一一对应的。

图3是深度为20m，距离分别为10km和20km处的接收信号的时频图。观察10km和20km的各阶模态时间差明显不同，随着距离的增大，时间差变大，因此可用不同距离处各阶模态时间差的差异唯一确定声源距离。

步骤4：声源深度距离估计。

利用步骤3得到的实际接收信号的时频矩阵与每个网格点处假定声源信号的时频矩阵，同时也就计算出了实际接收信号和模拟接收信号关于声源深度的各阶模态能量比以及关于声源距离的各阶模态时间差。

对步骤3的实际接收信号时频矩阵X_r进行列向量化得：

K_r＝[k_r1,k_r2,…,k_rM]

然后对K_r求转置再进行矩阵合并：

对上式中的各个元素取模：

U_r＝[|u_r1|,|u_r2|,…,|u_r(N×M)|]

对步骤3的模拟接收信号时频矩阵X_c采用上述相同的处理得到U_c：

U_c＝[|u_c1|,|u_c2|,…,|u_c(N×M)|]

计算得到U_r的协方差矩阵R_r：

对R_r进行矩阵特征分解得到噪声子空间：

N_r＝[q_k,q_k+1,…,q_M×N]

其中，q_k为对协方差矩阵R_r进行特征分解得到的第K个小特征值对应的特征向量。

计算得到U_c的协方差矩阵R_c：

对R_c进行特征分解得到信号子空间S_c：

S_c＝[q₁,q₂,…,q_k-1]

计算得到步骤3中(N×M)维的时频矩阵X_r的协方差矩阵R_xr：

对R_xr进行特征分解求得噪声子空间N_xr：

N_xr＝[q_k,q_k+1,…,q_M]

计算得到步骤3中(N×M)维的时频矩阵X_c的协方差矩阵R_xc：

对R_xc进行特征分解得到信号子空间S_xc：

S_xc＝[q₁,q₂,…,q_k-1]

其中S_xc中的元素为对协方差矩阵R_xr进行特征分解得到的K-1个大特征值对应的特征向量。

然后使用下列六种方法估计水下目标位置。

方法一：

将步骤3计算得到X_r和X_c代入下式实现对声源的定位。

方法二：

将上述计算得到的U_c和R_r代入下式实现声源定位。

方法三：

将上述计算得到的U_c和R_r代入下式实现声源定位。

方法四：

将上述计算得到的U_c和N_r代入下式实现声源定位。

方法五：

将上述计算得到的S_c和N_r代入下式实现声源定位。

方法六：

将上述计算得到的S_xc和N_xr代入下式实现声源定位。

图5、6、7、8、9、10分别是以上六种方法的定位结果。所提出的方法能够定位声源位置。

Claims (6)

1.基于时频分布相关的水下瞬态声源深度距离估计方法，其特征在于，包括步骤如下：

步骤1：根据海洋环境模型对声源目标区域进行网格化处理；

步骤2：计算接收信号；

声源发射一个低频宽带信号，频率带宽为f₁～f₂，经过海洋波导传播至距离r和深度z处对应的接收信号x(r,z,t)为：

其中，j为虚数，r为接收点与声源的水平距离，z为接收点深度，z_s为声源深度，ξ_n(ω)为频率为ω的声源波数k₀的水平分量；Z_n(z,ω)和Z_n(z_s,ω)是接收点和发射点处所对应的第n阶模态振幅，ω_min为最小声源频率，ω_max为最大声源频率，t为时间序列，S(ω)为发生信号频谱；

分别计算步骤1得到的网格中假定声源和实际声源在接收点处的模态振幅和声源波数的水平分量；带入公式(1)，分别得到网格点处假定声源和实际声源在接收点处的模拟接收时域信号和实际接收时域信号；

步骤3：对接收信号进行时频分析；

将步骤2得到的模拟接收时域信号和实际接收时域信号分别进行时频分析，得到实际接收信号归一化时频矩阵X_r和网格点处假定声源的模拟接收信号归一化时频矩阵X_c；

步骤4：声源深度距离估计；

利用步骤3得到的实际接收信号的时频矩阵X_r与每个网格点处假定声源信号的时频矩阵X_c，同时计算出实际接收信号和模拟接收信号关于声源深度的各阶模态能量比以及关于声源距离的各阶模态时间差；

将步骤3计算得到X_r和X_c代入下式实现对声源的定位；

式中x_rnm为实际接收信号时频矩阵X_r中的第n行、第m列的元素，x_cnm(r,z)为在实测海洋环境条件下，以距离r和深度z仿真计算所得的数据时频矩阵X_c中的第n行、第m列的元素，N和M分别为时-频矩阵的行数和列数。

2.根据权利要求1所述的基于时频分布相关的水下瞬态声源深度距离估计方法，其特征在于，步骤4：声源深度距离估计；

对步骤3的实际接收信号时频矩阵X_r进行列向量化得：

K_r＝[k_r1,k_r2,…,k_rM]

然后对K_r求转置再进行矩阵合并：

对上式中的各个元素取模：

U_r＝[|u_r1|,|u_r2|,…,|u_r(N×M)|]

对步骤3的模拟接收信号时频矩阵X_c采用上述相同的处理得到U_c：

U_c＝[|u_c1|,|u_c2|,…,|u_c(N×M)|]

计算得到U_r的协方差矩阵R_r：

使用下述方法估计水下目标位置；

将计算得到的U_c和R_r代入下式实现声源定位：

3.根据权利要求2所述的基于时频分布相关的水下瞬态声源深度距离估计方法，其特征在于，使用下述方法估计水下目标位置；

将计算得到的U_c和R_r代入下式实现声源定位：

4.根据权利要求2所述的基于时频分布相关的水下瞬态声源深度距离估计方法，其特征在于，对R_r进行矩阵特征分解得到噪声子空间：

N_r＝[q_k,q_k+1,…,q_M×N]

其中，q_k为对协方差矩阵R_r进行特征分解得到的第k个特征值对应的特征向量；

将上述计算得到的U_c和N_r代入下式实现声源定位；

5.根据权利要求2所述的基于时频分布相关的水下瞬态声源深度距离估计方法，其特征在于，计算得到U_c的协方差矩阵R_c：

对R_c进行特征分解得到信号子空间S_c：

S_c＝[q₁,q₂,…,q_k-1]；

将上述计算得到的S_c和N_r代入下式实现声源定位；

6.根据权利要求1所述的基于时频分布相关的水下瞬态声源深度距离估计方法，其特征在于，计算步骤3中(N×M)维的时频矩阵X_r的协方差矩阵R_xr：

对R_xr进行特征分解求得噪声子空间N_xr：

N_xr＝[q_k,q_k+1,…,q_M]

计算步骤3中(N×M)维的时频矩阵X_c的协方差矩阵R_xc：

对R_xc进行特征分解得到信号子空间S_xc：

S_xc＝[q₁,q₂,…,q_k-1]

其中S_xc中的元素为对协方差矩阵R_xr进行特征分解得到的第1到第k-1个特征值对应的特征向量；

将上述计算得到的S_xc和N_xr代入下式实现声源定位；

Images