CN115902849A - 一种基于波束输出强度重采样的深海声源深度估计方法 - Google Patents

Abstract

本发明涉及一种基于波束输出强度重采样的深海声源深度估计方法基于波束输出强度重采样的深海声源深度估计方法，基于深海可靠声路径环境下直达和海面反射路径的干涉特征，提出了一种适用于该环境下垂直线列阵的近海面单频运动声源深度估计方法。所提方法基于声速剖面信息，对波束输出强度进行重采样，随后对重采样后的波束输出强度进行修正傅立叶变换，变换输出的峰值位置即为声源深度估计结果。相比于已有的修正傅立叶变换声源深度估计方法。本发明所提方法能够有效估计近海面附近单频运动声源的深度，在所给出的典型实施例中，本发明所提方法(相比于原始方法)能够将深度估计偏差从13.01m降低至1.23m，对应相对估计偏差从16.25％降至1.54％。

Description

一种基于波束输出强度重采样的深海声源深度估计方法

技术领域

本发明属于海洋工程、水声工程、阵列信号处理和声呐技术等领域，涉及一种基于波束输出强度重采样的深海声源深度估计方法，适用于深海大深度垂直线列阵对近海面单频运动声源深度的被动估计问题。

背景技术

深海声源深度估计是声源定位中一个重要且具有挑战性的问题，一直以来受到广泛关注。在深海可靠声路径环境下(声源位于海面附近，接收阵列布放在海底附近)，阵列接收信号主要由直达路径(D)和海面反射路径(SR)组成，两条路径之间的干涉引起波束输出的强度起伏特征与声源深度密切相关。因此，可靠声路径环境下，直达-海面反射路径的干涉(D-SR干涉)特征是估计深海声源深度行之有效的方法之一。

对于宽带声源，D-SR干涉体现在接收信号在频域内的强度起伏，在声源信号垂直到达角已知的情况下，声源深度能够通过匹配波束形成输出信号的频谱幅度起伏特征与模型计算的理论频谱幅度得到。对于窄带/单频运动声源，D-SR干涉特征体现在波束输出强度随声源信号垂直到达角(或声源距离)的起伏。在窄带条件下，声源深度不仅可以通过匹配的方法得到，还可以通过McCargar和Zurk所提出的修正傅立叶变换方法进行估计[R.McCargar and L.M.Zurk,"Depth-based signal separation with vertical linearrays in the deep ocean,"J.Acoust.Soc.Am.,vol.133,no.4,pp.EL320-EL325,2013.]。该方法基于劳埃德镜理论，对波束输出强度序列进行修正傅立叶变换，变换输出的峰值直接对应声源深度。

相比于匹配类方法，基于修正傅立叶变换的声源深度估计方法具有实际应用方便，计算速度快等优势。但是，现有的基于修正傅立叶变换的声源深度估计方法是在劳埃德镜理论(即等声速假设)下推导得到的，该方法在实际随深度变化的声速剖面环境中可能会存在估计偏差。

发明内容

要解决的技术问题

为了避免现有技术的不足之处，本发明提出一种基于波束输出强度重采样的深海声源深度估计方法，解决现有基于修正傅立叶变换的深海单频声源深度估计方法在实际随深度变化声速剖面环境中估计偏差较大的问题。首先基于声速剖面信息对波束输出强度进行重采样，而后使用修正傅立叶变换进行声源深度估计。相比于现有的修正傅立叶变换声源深度估计方法，本发明所提方法在典型深海环境下具有更高的估计精度。

技术方案

一种基于波束输出强度重采样的深海声源深度估计方法，其特征在于：在海底布放垂直线列阵接收海面附近径向运动声源发出的单频信号，垂线阵由M个阵元组成，阵元间隔为d，阵列中心深度为z_r；声源深度估计步骤如下：

步骤1：初始时刻声源距离垂线阵的水平距离为r₀，声源频率为f，深度为z_s；观测时间内声源的水平径向运动孔径为r_SPAN，观测结束时声源水平距离为r₀+r_SPAN；

总观测时间为T，波束输出强度即方位谱观测间隔为t_intv，在声源运动过程中一共观测K＝ceil(T/t_intv)次波束输出强度即方位谱，其中ceil(·)表示向下取整；

步骤2：针对第k个观测时间点，k＝1,2,…,K，对阵列接收信号进行常规波束形成，获得声源信号的波束形成方位谱：

B^k(θ)＝|Y^k(θ,f)|²

其中：Y^k(θ,f)为频点f处指向角度θ的波束形成器输出；

步骤3：对所有时刻的波束输出方位谱B^k(θ)，在时间-方位平面上提取声源轨迹，流程为：以所有时刻波束输出强度即方位谱峰值的归一化强度最大的L个观测时刻，记为L₁,L₂,…,L_L；

提取L个时刻声源对应的垂直到达角，即方位谱

峰值对应的角度，记为θ_L1,θ_L2,…,θ_LL，对应声源在t_L1,t_L2,…,t_LL时刻的角度；

拟合所有K个观测时间点的声源信号接收角，得到声源信号在k＝1,2,…,K时刻的垂直到达角，记为θ₁,θ₂,…,θ_K；

K个观测时间点的声源信号接收角即为声源在时间-方位平面上提取声源轨迹；

步骤4：对k＝1,2,…,K时刻，提取常规波束形成在θ₁,θ₂,…,θ_K方向上的波束输出强度，将其与sin(θ₁),sin(θ₂),…,sin(θ_K)对应，获得接收角正弦域内的波束输出强度序列，记为x(u)＝B^k(θ_k)，u＝sin(θ₁),…,sin(θ_K)；

步骤5：对波束输出强度序列x(u)进行重采样；

步骤6：对重采样后的波束输出强度序列进行修正的傅里叶变换，重采样序列修正傅里叶变换输出X_rsp(z)的峰值所在位置即为声源深度的估计结果。

所述步骤2声源信号的波束形成方位谱的计算为：对每个时刻阵元接收信号进行快速傅立叶变换，得到频域接收信号，提取频率为f处的接收信号，记k时刻，第j个阵元的频域接收信号为

频率f处指向角度θ的常规波束形成器权值：

其中上标“T”表示转置运算，c为波束形成参考声速，i为虚数单位；第k时刻，频点f处指向角度θ的波束形成器输出为：

其中上标“H”表示共轭转置运算，“×”表示向量相乘；

则波束形成方位谱为：B^k(θ)＝|Y^k(θ,f)|²。

所述步骤5的对波束输出强度序列x(u)进行重采样：x_rsp(u)＝x[h(u)]，其中h(u)表示重采样算子：h(u)＝ξ^-1(u)，ξ^-1(u)表示的ξ(u)反函数，ξ(u)如下：

其中：c_zr为垂线阵中心点深度处的声速，z_s-orig为一个初始声源深度估计结果，由现有的修正傅立叶变换方法提供，c_w为一个声速平均值：

其中，z_s-max为最大可能声源深度，c(z)为深度z处的海水声速。

所述步骤6对重采样后的波束输出强度序列进行修正的傅里叶变换：

其中u_min-rsp＝ξ(sinθ₁)和u_max-rsp＝ξ(sinθ_K)分别表示重采样后序列x_rsp(u)中u的取值下限和上限；u_n为重采样波束强度序列中的第n个采样点，N为总的重采样点数，Δu_rsp为重采样波束强度序列的采样间隔。

所述声源深度不小于10m。

所述声源径向运动孔径r_SPAN不小于3km。

所述声源频率不小于100Hz。

有益效果

本发明提出的一种基于波束输出强度重采样的深海声源深度估计方法，基于深海可靠声路径环境下直达和海面反射路径的干涉特征，提出了一种适用于该环境下垂直线列阵的近海面单频运动声源深度估计方法。所提方法基于声速剖面信息，对波束输出强度进行重采样，随后对重采样后的波束输出强度进行修正傅立叶变换，变换输出的峰值位置即为声源深度估计结果。相比于已有的修正傅立叶变换声源深度估计方法。

本发明所提方法具有更高的估计精度，同时还保留了计算速度快，实际应用方便的优点。本发明所提方法的基本原理和实施方案经过了计算机数值仿真的验证，其结果表明，在典型深海环境下，本发明所提方法能够有效估计近海面附近单频运动声源的深度，在所给出的典型实施例中，本发明所提方法(相比于原始方法)能够将深度估计偏差从13.01m降低至1.23m，对应相对估计偏差从16.25％降至1.54％。

附图说明

图1是仿真场景声速剖面示意图。

图2是声源信号以及阵列接收信号仿真示例，子图(a)为声源信号时域波形及时频分析结果，子图(b)为声源距离为6km时的阵列接收信号。

图3是本发明所提一种基于波束输出强度重采样的深海声源深度估计方法流程图。

图4是所有时刻阵列接收信号常规波束形成方位谱(归一化结果)。图中黑色圆圈表示归一化强度大于-3dB的波束输出峰值位置，黑色虚线表示提取到的声源轨迹。

图5为由图4中提取到的接收角正弦域内的波束输出强度归一化结果。

图6为重采样后的接收角度正弦域内的波束输出强度归一化结果。

图7为原始的修正傅立叶变换声源深度估计方法结果。

图8为本发明所提方法声源深度估计结果。

具体实施方式

现结合实施例、附图对本发明作进一步描述：

本发明解决其技术问题所采用的技术方案：一种基于波束输出强度重采样的深海声源深度估计方法，其特点是包括下述步骤：

步骤一：在典型深海环境下，在海底附近布放一套垂直线列阵接收海面附近径向运动声源发出的单频信号。垂线阵由M个阵元组成，阵元间隔为d，阵列中心深度为z_r。声源深度为z_s，频率为f，初始时刻声源距离垂线阵的水平距离为r₀，观测时间内声源的水平径向运动孔径为r_SPAN，观测结束时声源水平距离为r₀+r_SPAN。总观测时间为T，波束输出强度观测间隔为t_intv，在声源运动过程中一共观测K＝ceil(T/t_intv)次波束输出强度，其中ceil(·)表示向下取整。

步骤二：针对第k个观测时间点，k＝1,2,…,K，对阵列接收信号进行常规波束形成，获得声源信号的波束输出强度。其具体流程如下：

针对每个时刻，对阵元接收信号进行快速傅立叶变换，得到频域接收信号，提取频率为f处的接收信号。记k时刻，第j个阵元的频域接收信号为

频率f处指向角度θ的常规波束形成器权值为

其中上标“T”表示转置运算，c为波束形成参考声速，i为虚数单位。第k时刻，频点f处指向角度θ的波束形成器输出为

其中上标“H”表示共轭转置运算，“×”表示向量相乘。

波束形成方位谱由下式计算：

B^k(θ)＝|Y^k(θ,f)|²  (3)

步骤三：联合所有时刻的波束输出方位谱B^k(θ)，在时间-方位平面上提取声源轨迹。具体流程如下：

由于声源信号的波束强度会随时间起伏，因此可以找到所有波束输出强度较大的(方位谱峰值的归一化强度)L个观测时刻，记为L₁,L₂,…,L_L。提取这L个时刻声源对应的垂直到达角[即方位谱

峰值对应的角度]，记为θ_L1,θ_L2,…,θ_LL，对应声源在t_L1,t_L2,…,t_LL时刻的角度。根据上述L个时刻的声源角度，拟合所有K个时刻的声源信号接收角，得到声源信号在k＝1,2,…,K时刻的垂直到达角，记为θ₁,θ₂,…,θ_K。

步骤四：对k＝1,2,…,K时刻，提取常规波束形成在θ₁,θ₂,…,θ_K方向上的波束输出强度，将其与sin(θ₁),sin(θ₂),…,sin(θ_K)对应，获得接收角正弦域内的波束输出强度序列，记为x(u)＝B^k(θ_k)，u＝sin(θ₁),…,sin(θ_K)。

步骤五：按照下式对波束输出强度序列x(u)进行重采样：

x_rsp(u)＝x[h(u)]  (4)

其中h(u)表示重采样算子，由下式确定

h(u)＝ξ^-1(u)  (5)

其中ξ^-1(u)表示的ξ(u)反函数

上式中，c_zr为垂线阵中心点深度处的声速，z_s-orig为一个初始声源深度估计结果，可以由现有的修正傅立叶变换方法提供，c_w为一个声速平均值，由下式计算：

其中，z_s-max为最大可能声源深度，c(z)为深度z处的海水声速。

步骤六：按照下式，对重采样后的波束输出强度序列进行修正的傅里叶变换：

其中u_min-rsp＝ξ(sinθ₁)和u_max-rsp＝ξ(sinθ_K)分别表示重采样后序列x_rsp(u)中u的取值下限和上限。u_n为重采样波束强度序列中的第n个采样点，N为总的重采样点数，Δu_rsp为重采样波束强度序列的采样间隔。

最后，重采样序列修正傅里叶变换输出X_rsp(z)的峰值所在位置即为声源深度的估计结果。

其适用于径向运动的窄带/单频声源，声源径向运动可以是由远及近或由近及远。

具体实施例：

1.深海波导环境、声源以及接收垂线阵配置

为验证本发明方法的有效性，利用计算机进行仿真实验。本实施例考虑一个典型深海环境，海深为3950m，海水声速剖面如附图1所示，海水密度为1.0g/cm³；海底半空间声速为1600m/s，密度为1.5g/cm³，海底底质压缩波衰减系数为0.14dB/λ。接收所使用的垂直线列阵由16个阵元组成，阵元间隔为2m，阵列中心深度为3716m。声源深度为80m，辐射频率f＝535Hz的单频信号，初始水平距离为r₀＝6km，水平径向运动孔径为r_SPAN＝5km，声源径向运动速度为5m/s。声源运动期间垂线阵对声源信号波束输出强度进行1000次观测，观测间隔为t_intv＝1s，观测水平距离间隔为5m。

2.垂直线列阵接收信号

本实施例假设声源信号为一个单频CW脉冲信号，脉冲长度为0.5s，幅度为1，其时域波形和时频分析结果如附图2(a)所示。阵列采样频率为f_s＝10kHz，阵列在声源信号发出的同时开始采集，接收机开机时间为20s。

阵列接收信号的仿真方法为：针对某一时刻的声源距离，使用Bellhop声场模型计算出声源与接收阵列(所有传感器)之间声线的到达时间和幅度。针对某一传感器，使用这些声线到达结构构成信道冲击响应，将声源信号与信道冲击响应进行时域卷积，即可得到该传感器接收信号。最后根据信噪比在仿真的接收信号中添加噪声。假设各个传感器在T_R内的接收信噪比均为SNR＝0dB。依次对所有接收时刻的所有接收阵元进行上述仿真操作，即可得到仿真的阵列接收信号。声源距离为6km时的所有阵元的接收信号如附图2(b)所示。

3.一种基于波束输出强度重采样的深海声源深度估计方法

如附图3所示，本发明所提一种基于波束输出强度重采样的深海声源深度估计方法具体实施流程如下：

步骤一：在典型深海环境下，在海底附近布放一套垂直线列阵接收海面附近径向运动声源发出的窄带单频信号。垂线阵由M＝16个阵元组成，阵元间隔为d＝2m，阵列中心深度为z_r＝3716m。声源深度为z_s＝80m，频率为f＝535Hz，初始时刻声源距离垂线阵的水平距离为r₀＝6km，观测时间内声源的水平径向运动孔径为r_SPAN＝5km，观测结束时声源水平距离为r₀+r_SPAN＝11km。总观测时间为T＝1000s，波束输出强度观测间隔为t_intv＝1s，在声源运动过程中一共观测K＝ceil(T/t_intv)＝1000次，其中ceil(·)表示向下取整。本实施例中这一步已通过Bellhop模型的声场模型仿真完成，最终得到1000组对应声源不同位置的阵列接收信号，每组阵列接收信号的采样点数均为2×10⁵。

步骤二：针对第k个观测时间点，k＝1,2,…,1000，对阵列接收信号进行常规波束形成，获得声源信号的波束输出强度。其具体流程如下：

针对每个时刻，对阵元接收信号进行快速傅立叶变换，得到频域接收信号，提取频率为f＝535Hz处的接收信号，其具体操作方法为：对长度为2×10⁵的阵元接收信号进行快速傅立叶变换，得到长度为2×10⁵的变换结果，提取第10701个点的值为535Hz处的接收信号。记k时刻，第j个阵元的频域接收信号为

频率f处指向角度θ的常规波束形成器权值为

其中上标“T”表示转置运算，c＝1521.4m/s为波束形成参考声速，i为虚数单位。第k时刻，频点f处指向角度θ的波束形成器输出为

其中上标“H”表示共轭转置运算，“×”表示向量相乘。

波束形成方位谱由下式计算：

B^k(θ)＝|Y^k(θ,f)|²   (11)

通过上述步骤，所有观测时刻的波束形成方位谱组成的时间-方位平面如附图4所示。

步骤三：联合所有时刻的波束输出方位谱B^k(θ)，在时间-方位平面上提取声源轨迹。具体流程如下：

如附图4所示，声源信号的波束强度会随时间起伏，因此可以找到所有波束输出强度较大(归一化强度大于-3dB)的212个观测时刻，记为L₁,L₂,…,L₂₁₂。提取这212个时刻声源对应的垂直到达角[即方位谱

峰值对应的角度]，记为θ_L1,θ_L2,…,θ_L212，对应声源在t_L1,t_L2,…,t_L212时刻的角度(如附图4中黑色圆圈所示)。根据上述212个时刻的声源角度，对所有1000个时刻的声源信号接收角进行拟合，得到声源信号在k＝1,2,…,1000时刻的垂直到达角，记为θ₁,θ₂,…,θ₁₀₀₀，如附图4中黑色虚线所示。

步骤四：对k＝1,2,…,1000时刻，分别提取常规波束形成在θ₁,θ₂,…,θ₁₀₀₀方向上的波束输出强度，将其与sin(θ₁),sin(θ₂),…,sin(θ₁₀₀₀)对应，获得接收角度正弦域内的波束输出强度序列，记为x(u)＝B^k(θ_k)，u＝sin(θ₁),…,sin(θ₁₀₀₀)。针对本实施例，提取到的接收角正弦域内的波束输出强度序列归一化结果如附图5所示，u的取值范围为0.2804～0.5079。

步骤五：按照下式对波束输出强度序列x(u)进行重采样：

x_rsp(u)＝x[h(u)]  (12)

其中h(u)表示重采样算子，由下式确定

h(u)＝ξ^-1(u)  (13)

其中

上式中，c_zr为垂线阵中心点深度处的声速，为1521.4m/s，z_s-orig＝93.01m为一个初始声源深度估计结果，由现有的修正傅立叶变换方法提供(其结果如附图7所示)，c_w为一个声速的平均值，由下式计算：

其中，z_s-max＝300m为最大可能声源深度，c(z)为深度z处的海水声速。

重采样后的波束输出强度序列如附图6所示，重采样后u的取值范围为0.2119～0.4732。u域的重采样间隔为3×10^-4，重采样点数为1000点。

步骤六：按照下式，对重采样后的波束输出强度序列进行修正的傅里叶变换：

其中u_min-rsp＝ξ(sinθ₁)＝0.211和u_max-rsp＝ξ(sinθ_K)＝0.4732分别表示重采样后序列x_rsp(u)中u的取值下限和上限。u_n为重采样波束强度序列中的第n个采样点，Δu_rsp＝3×10^-4为重采样波束强度序列的采样间隔。

重采样序列修正傅里叶变换输出X_rsp(z)如附图8所示，其峰值位置即为声源深度估计结果，为81.23m。相比于原始方法(估计结果为93.01m)具有更高的估计精度。可以看出，相比于原始方法，本发明所提方法能够将估计偏差从13.01m降低至1.23m，对应相对估计偏差从16.25％降至1.54％。

Claims (7)

1.一种基于波束输出强度重采样的深海声源深度估计方法，其特征在于：在海底布放垂直线列阵接收海面附近径向运动声源发出的单频信号，垂线阵由M个阵元组成，阵元间隔为d，阵列中心深度为z_r；声源深度估计步骤如下：

步骤1：初始时刻声源距离垂线阵的水平距离为r₀，声源频率为f，深度为z_s；观测时间内声源的水平径向运动孔径为r_SPAN，观测结束时声源水平距离为r₀+r_SPAN；

总观测时间为T，波束输出强度即方位谱观测间隔为t_intv，在声源运动过程中一共观测K＝ceil(T/t_intv)次波束输出强度即方位谱，其中ceil(·)表示向下取整；

步骤2：针对第k个观测时间点，k＝1,2,…,K，对阵列接收信号进行常规波束形成，获得声源信号的波束形成方位谱：

B^k(θ)＝|Y^k(θ,f)|²

其中：Y^k(θ,f)为频点f处指向角度θ的波束形成器输出；

步骤3：对所有时刻的波束输出方位谱B^k(θ)，在时间-方位平面上提取声源轨迹，流程为：以所有时刻波束输出强度即方位谱峰值的归一化强度最大的L个观测时刻，记为L₁,L₂,…,L_L；

提取L个时刻声源对应的垂直到达角，即方位谱

峰值对应的角度，记为θ_L1,θ_L2,…,θ_LL，对应声源在t_L1,t_L2,…,t_LL时刻的角度；

拟合所有K个观测时间点的声源信号接收角，得到声源信号在k＝1,2,…,K时刻的垂直到达角，记为θ₁,θ₂,…,θ_K；

K个观测时间点的声源信号接收角即为声源在时间-方位平面上提取声源轨迹；

步骤4：对k＝1,2,…,K时刻，提取常规波束形成在θ₁,θ₂,…,θ_K方向上的波束输出强度，将其与sin(θ₁),sin(θ₂),…,sin(θ_K)对应，获得接收角正弦域内的波束输出强度序列，记为x(u)＝B^k(θ_k)，u＝sin(θ₁),…,sin(θ_K)；

步骤5：对波束输出强度序列x(u)进行重采样；

步骤6：对重采样后的波束输出强度序列进行修正的傅里叶变换，重采样序列修正傅里叶变换输出X_rsp(z)的峰值所在位置即为声源深度的估计结果。

2.根据权利要求1所述基于波束输出强度重采样的深海声源深度估计方法，其特征在于：所述步骤2声源信号的波束形成方位谱的计算为：对每个时刻阵元接收信号进行快速傅立叶变换，得到频域接收信号，提取频率为f处的接收信号，记k时刻，第j个阵元的频域接收信号为

频率f处指向角度θ的常规波束形成器权值：

其中上标“T”表示转置运算，c为波束形成参考声速，i为虚数单位；第k时刻，频点f处指向角度θ的波束形成器输出为：

其中上标“H”表示共轭转置运算，“×”表示向量相乘；

则波束形成方位谱为：B^k(θ)＝|Y^k(θ,f)|²。

3.根据权利要求1所述基于波束输出强度重采样的深海声源深度估计方法，其特征在于：所述步骤5的对波束输出强度序列x(u)进行重采样：x_rsp(u)＝x[h(u)]，其中h(u)表示重采样算子：h(u)＝ξ^-1(u)，ξ^-1(u)表示的ξ(u)反函数，ξ(u)如下：

其中：c_zr为垂线阵中心点深度处的声速，z_s-orig为一个初始声源深度估计结果，由现有的修正傅立叶变换方法提供，c_w为一个声速平均值：

其中，z_s-max为最大可能声源深度，c(z)为深度z处的海水声速。

4.根据权利要求1所述基于波束输出强度重采样的深海声源深度估计方法，其特征在于：所述步骤6对重采样后的波束输出强度序列进行修正的傅里叶变换：

其中u_min-rsp＝ξ(sinθ₁)和u_max-rsp＝ξ(sinθ_K)分别表示重采样后序列x_rsp(u)中u的取值下限和上限；u_n为重采样波束强度序列中的第n个采样点，N为总的重采样点数，Δu_rsp为重采样波束强度序列的采样间隔。

5.根据权利要求1所述基于波束输出强度重采样的深海声源深度估计方法，其特征在于：所述声源深度不小于10m。

6.根据权利要求1所述基于波束输出强度重采样的深海声源深度估计方法，其特征在于：所述声源径向运动孔径r_SPAN不小于3km。

7.根据权利要求1所述基于波束输出强度重采样的深海声源深度估计方法，其特征在于：所述声源频率不小于100Hz。

Images