CN113011006B - 一种基于互相关函数脉冲波形匹配的目标深度估计方法 - Google Patents

Abstract

本发明公布了一种基于互相关函数脉冲波形匹配的目标深度估计方法，所述方法包括：步骤1)利用两条近似平行的声学阵列，从被测目标辐射噪声数据中提取时域伪格林函数，获得其脉冲包络；步骤2)沿深度和距离方向划分网格，然后根据声学阵列布放海域的水文环境，计算各网格点上对应的时域伪格林函数，获取每个网格点的脉冲包络；步骤3)将步骤1)提取的脉冲包络与各网格点上脉冲包络进行匹配处理，获得水中目标的深度估计结果。本发明的方法将目标深度估计和目标距离估计分离，可以在未知目标距离的情况下获得目标深度估计结果,克服了传统匹配场处理方法需要同时估计目标距离和目标深度或者先估计出目标距离然后再估计目标深度的限制。

Description

一种基于互相关函数脉冲波形匹配的目标深度估计方法

技术领域

本发明属于海洋声学领域，具体涉及一种基于互相关函数脉冲波形匹配的目标深度估计方法。

背景技术

目标的深度信息是非常重要的目标参数。根据声学传感器接收信号给出可靠的目标深度估计结果，是解决水中目标深度分辨的有效途径。匹配场处理方法是对目标进行深度估计的方法之一。这种方法充分利用了声场的干涉特性，采用声传播模型来计算拷贝场向量，然后将拷贝场与测量场进行匹配处理，从而获得未知变量的匹配处理结果。自从1976年匹配场(MFP)定位概念提出以来，国内外已发表了数百篇相关论文，进行了大量的海上实验，取得了不少进展。但是，海洋环境条件十分复杂，环境参数往往是时变、空变的，匹配场处理方法面临海洋环境参数无法准确测量的难题，从而拷贝场计算不准确，导致匹配处理结果与真实结果偏差太大，甚至无法给出置信度超过设定阈值的处理结果。也有学者提出利用声能流有功分量和无功分量正负号的变化来判别水上水下目标的方法，该方法要求将接收水听器布设在特定深度上，此深度需要根据声场分析频率和水深、底质、声速剖面等海洋环境计算得到，而实际应用中很难精确计算此深度，限制了该方法的使用。另外，可以利用脉冲信号的频散特征对目标声源进行深度估计，这种方法要求声源发出的是时域上很短的脉冲信号，而实际应用中所接收到的往往是连续的辐射噪声信号，各类船只很少会辐射这类持续时间很短但强度很大的瞬态脉冲信号。

发明内容

针对声学目标深度分辨这一需求，本发明根据浅海低频宽带声场的水平纵向互相关函数的物理特征，提出了一种基于互相关函数脉冲波形匹配的目标深度估计方法，适用于浅海海域的低频宽带被动目标深度估计，为提高浅海海域的目标深度分辨能力提供技术手段。该方法根据浅海低频宽带声场的水平纵向互相关函数的物理特征，利用两条近似平行的声学阵列提取时域伪格林函数，通过将实际提取的伪格林函数脉冲包络与在不同深度上仿真计算的伪格林函数脉冲包络进行匹配处理，给出水中目标的深度估计结果。

为了实现上述目的，本发明提供了一种基于互相关函数脉冲波形匹配的目标深度估计方法，所述方法包括：

步骤1)利用两条近似平行的声学阵列，从被测目标辐射噪声数据中提取时域伪格林函数，获得其脉冲包络；

步骤2)沿深度和距离方向划分网格，然后根据声学阵列布放海域的水文环境，计算各网格点上对应的时域伪格林函数，获取每个网格点的脉冲包络；

步骤3)将步骤1)提取的脉冲包络与各网格点上脉冲包络进行匹配处理，获得水中目标的深度估计结果。

作为上述方法的一种改进，所述步骤1)具体包括：

步骤1-1)对声学阵列A一段时间内的多拍数据进行处理，获得声学阵列A的宽带波束输出结果；获得某一目标相对于声学阵列A的随时间t变化的方位历程θ_A(t)；

步骤1-2)对声学阵列B一段时间内的多拍数据进行处理，获得声学阵列B的宽带波束输出结果；获得某一目标相对于声学阵列B的随时间t变化的方位历程θ_B(t)；

步骤1-3)若|θ_A(t)-θ_B(t)|≤α，α为阈值；转入步骤1-4)；否则改变时间段，转入步骤1-1)，直到在宽带波束输出结果中有满足|θ_A(t)-θ_B(t)|≤α的目标出现；

步骤1-4)对两条声学阵列在目标方位的波束形成输出结果进行预处理；对预处理后的波束数据进行互相关运算，并在一定时间内积累，求取其平均值，即为伪格林函数。

作为上述方法的一种改进，所述步骤1-1)具体包括：

步骤1-1-1)将声学阵列A各阵元接收到的声信号分成多拍，每拍时间长度为ΔT，前后相邻两拍重叠P％，并对每一拍信号做傅里叶变换得到频域信号；声学阵列A的频域信号为s_A(f；t)＝[s_A1(f；t),s_A2(f；t),...,s_AM(f；t)]^T；f是信号频率，t是当前拍数据的起始时刻；M是声学阵列A的阵元个数；(·)^T表示转置；

步骤1-1-2)以声学阵列A的几何中心为参考点，计算声学阵列A各频点的在波束形成扫描方位θ的加权向量w_A(f,θ)＝[w_A0(f,θ),w_A1(f,θ),...,w_AM-1(f,θ)]^T，其中

其中，τ_Am(θ)表示声学阵列A的第m+1个阵元相对于参考点的时延：

v(θ)＝-[cosθ,sinθ]^T

p_Am＝[p_xm,p_ym]^T

v(θ)是入射信号的单位向量，p_Am为第m+1个阵元相对于参考点的二维坐标；

步骤1-1-3)计算声学阵列A指向方位θ的波束形成结果b_A(f,θ；t)：

b_A(f,θ；t)＝w_A ^H(f,θ)s_A(f；t)

其中，(·)^H表示共轭转置；θ的取值范围为[θ₀-45°，θ₀+45°]，θ₀表示声学阵列A和声学阵列B的参考阵元连线的角度；

步骤1-1-4)对声学阵列A的一段时间[t₁，t₂]内的多拍数据进行步骤1-1-1)至步骤1-1-3)的处理，获得声学阵列A的宽带波束输出结果H_A(θ,t)：

其中，符号“Σ”表示求和；

步骤1-5)根据声学阵列A的宽带波束输出结果，获得某一目标相对于声学阵列A的随时间t变化的方位历程θ_A(t)。

作为上述方法的一种改进，所述步骤1-2)具体包括：

步骤1-2-1)将声学阵列B各阵元接收到的声信号分成多拍，每拍时间长度为ΔT，前后相邻两拍重叠P％，并对每一拍信号做傅里叶变换得到频域信号；声学阵列B的频域信号为s_B(f；t)＝[s_B1(f；t),s_B2(f；t),...,s_BN(f；t)]^T；f是信号频率，选取的范围为[f_L,f_H]；t是当前拍数据的起始时刻；N是声学阵列B的阵元个数；

步骤1-2-2)以声学阵列B的几何中心为参考点，计算声学阵列B各频点的在波束形成扫描方位θ的加权向量w_B(f,θ)＝[w_B0(f,θ),w_B1(f,θ),...,w_B,N-1(f,θ)]^T，其中

其中τ_Bn(θ)表示声学阵列B的第n+1个阵元相对于参考点的时延：

v(θ)＝-[cosθ,sinθ]^T

p_Bn＝[p_xn,p_yn]^T

v(θ)是入射信号的单位向量，p_Bn为第n+1个阵元相对于参考点的二维坐标；

步骤1-2-3)计算声学阵列B指向方位θ的波束形成结果b_B(f,θ；t)：

b_B(f,θ；t)＝w_B ^H(f,θ)s_B(f；t)

其中，(·)^H表示共轭转置；

步骤1-2-4)对声学阵列B的一段时间[t₁，t₂]内的多拍数据进行步骤1-2-1)至步骤1-2-3)的处理，计算声学阵列B的宽带波束输出结果H_B(θ,t)：

步骤1-2-5)根据声学阵列B的宽带波束输出结果获得同一目标相对于声学阵列B的随时间t变化的方位历程θ_B(t)。

作为上述方法的一种改进，所述步骤1-4)包括：

步骤1-4-1)对两条声学阵列在目标方位的波束形成输出结果进行预处理，预处理后声学阵列A、B在目标方位的波束形成输出

分别为：

其中，k为波数，k＝2πf/c₀；L为声学阵列A、B几何中心水平间距；

为分析时间段内起始时刻所选目标相对于声学阵列B的水平方位；

I_A(f)为所选取目标相对于声学阵列A的平均声强谱：

I_B(f)为所选取目标相对于声学阵列B的平均声强谱：

其中，N_t为总拍数；N_t的取值范围为[600/(ΔT(1-P％))，1800/(ΔT(1-P％))]；ΔT为单拍数据的时间长度；P％为前后相邻两拍重叠率；

步骤1-4-2)对预处理后的波束数据进行互相关运算，并在一定时间内积累，求取其平均值，获得频域伪格林函数

(·)^*表示共轭；

步骤1-4-3)对频域伪格林函数

做反傅里叶变换，得到时域伪格林函数

其中，ifft表示反傅里叶变换，Re表示取实部；

步骤1-4-4)时域伪格林函数

是一个脉冲函数，其包络为v(τ)：

其中，|·|表示取绝对值，H(·)表示希尔伯特变换。

作为上述方法的一种改进，所述步骤2)具体包括：

步骤2-1)按照[0,4H/5]间隔，将声学阵列布放海域在深度上划分为Q_z个网格，H为水深；按照[R_min,R_max]间隔，将声学阵列布放海域在距离上划分为Q_R个网格，R_min为满足声学阵列远场要求所需的最小距离，R_max为声学阵列对目标的最大理论探测距离；

步骤2-2)根据水深、声速剖面和海底参数这些海洋环境参数，采用声场计算方法计算简正波信息：第l号简正波的水平波数k_l＝μ_l+iη_l，μ_l、η_l分别为k_l的实部、虚部；第l号简正波的特征函数Ψ_l(z)；其中，l＝1,2，…，I；I为简正波总号数；

步骤2-3)根据简正波信息计算各网格点的时域伪格林函数；

对于第q个网格点，其时域伪格林函数为：

其中，z_q、r_q分别为第q个网格点对应的声源深度和至两条声学阵列几何中心连线中点的水平距离；τ表示时延，ifft表示反傅里叶变换，Re表示取实部；

为第q个网格点的频域伪格林函数：

其中，B(z₁,z₂,z_q)为归一化因子：

z₁,z₂表示两个水听器的深度，其水平间距为L，两个水听器与声源的水平间距分别为r₁,r₂；

步骤2-4)时域伪格林函数

是一个脉冲函数，其包络为v_r(τ；z_q,r_q)：

其中，|·|表示取绝对值，H(·)表示希尔伯特变换。

作为上述方法的一种改进，深度网格数Q_z的取值范围为[8H/5λ,40H/5λ]，λ为分析带宽内中心频率所对应的声波波长；距离网格数Q_R的取值范围为[(R_max-R_min)/10,(R_max-R_min)]。

作为上述方法的一种改进，所述步骤3)具体包括：

步骤3-1)计算第q个网格点的时域伪格林函数包络v(τ)与实际数据提取的时域伪格林函数包络v_r(τ；z_q,r_q)的相关系数W(z_q,r_q)：

其中，V(τ)和V_r(τ；z_q,r_q)分别是对v(τ)和v_r(τ；z_q,r_q)在时延维度上取滑动平均后的结果：

其中，Δτ为滑动平均窗口宽度；Δτ的取值范围为[5/(f_H-f_L)，20/(f_H-f_L)]；f_L、f_H分别为分析频段的下限频率、上限频率；

步骤3-2)对于所有的W(z_q,r_q)，其最大值对应网格点的深度即为被测目标的深度估计结果。

作为上述方法的一种改进，如果目标距离是已知的，记为r_s，那么仅对目标距离附近范围内的网格点进行局部最大值搜索即可；距离搜索的典型范围为[0.8r_s,1.2r_s]。

本发明的优势在于：

1、本发明的方法将目标深度估计和目标距离估计分离，可以在未知目标距离的情况下获得目标深度估计结果,克服了传统匹配场处理方法需要同时估计目标距离和目标深度或者先估计出目标距离然后再估计目标深度的限制；

2、本发明的方法通过匹配时域伪格林函数脉冲包络来获得目标深度估计结果，降低了对海洋环境参数精确度的要求，提高了宽容性。

附图说明

图1为本发明中水听器及声源的相对位置示意图；

图2为本发明的基于互相关函数脉冲波形匹配的目标深度估计方法的流程图；

图3为本发明一个实施例中海水平均声速剖面；水深200m，海面声速1544m/s,声速随深度逐渐减小,200m深度处声速约为1502m/s；海底声速1608m/s；图中实线(图中标记为ssp1)为计算测量场时所用声速剖面，虚线(图中标记为ssp2)为计算拷贝场时所用声速剖面；

图4(a)是本发明一个实施例中测量场的时域伪格林函数包络,分析带宽为200Hz至300Hz；

图4(b)为图4(a)对应拷贝场的时域伪格林函数包络随时延和目标深度的二维图像。

图5(a)为所对应目标深度估计结果为8m，实际目标深度为10m的示意图；

图5(b)为所对应目标深度估计结果为101m，实际目标深度为100m的示意图。

具体实施方式

针对声学目标深度分辨这一需求，本发明根据浅海低频宽带声场的水平纵向互相关函数的物理特征，提出了一种基于互相关函数脉冲波形匹配的目标深度估计方法，适用于浅海海域的低频宽带被动目标深度估计，为提高浅海海域的目标深度分辨能力提供技术手段。

本发明所采用的技术方案是：

根据简正波理论，深度为z_s的点声源在水平距离r、深度z处的产生的声场(格林函数)可表示为

其中，k_m为第m号简正波的水平波数，k_m＝μ_m(f)+jη_m(f)，μ_m、η_m分别为k_m的实部、虚部，均为频率f的函数；Ψ_m(z)为第m号简正波的特征函数。j为虚数符号；e为自然对数的基底，e＝2.718……。

如图1所示，两个圆点表示水听器，其水平间距为L，深度分别为z₁、z₂，与声源(五角星所示)水平间距分别为r₁、r₂；一个圆点位于两个水听器的中心位置，与声源水平间距为R；声源相对于两个水听器连线的方位为θ。

对于方位为θ的声源，利用图1所示两个水听器，可获得时域伪格林函数

其中，τ表示时延，ifft表示反傅里叶变换，Re表示取实部。

其中B(z₁,z₂,z_s)为归一化因子

符号“Σ”为求和算子。

公式(2)所表示的时域伪格林函数是一个脉冲函数，其波形与目标深度、水听器纵向间距及深度、海洋环境等参数有关。已知水听器纵向间距及深度、海洋环境等参数的情况下，可根据脉冲波形估计被测目标深度。

需要指出的是，以上论述是针对两个水听器的情形的。但是一般情况下，为了提高处理增益，实际应用中通常用两个近似平行的声学阵列代替上述的两个水听器，以获得空间增益，提高信噪比。

如图2所示，本发明提出了一种基于互相关函数脉冲波形匹配的目标深度估计方法，包括如下步骤：

S1，处理实际数据，提取时域伪格林函数，获得其脉冲包络v(τ)。

S1-1，对布设在水中的两条近似平行的声学阵列(声学阵列A、B)所接收的声学数据分别进行波束形成，其波束形成结果记为b_A(f,θ；t)和b_B(f,θ；t)。其中，f是信号频率，选取的范围为[f_L,f_H]；t是当前拍数据的起始时刻；θ为预成波束方位。

优选的，f_L、f_H的典型取值范围为[2c₀/H,20c₀/H]。c₀为参考声速，通常取为水体声速的均值。H为水深。

优选的，θ的典型取值范围为[θ₀-45°，θ₀+45°]。θ₀表示声学阵列A和声学阵列B的参考阵元连线的角度。

S1-2，根据上一步结果，选取声学阵列正横方位附近的水面船只作为辐射噪声的声源，其相对于声学阵列A、B的随时间t变化的方位历程分别记为θ_A(t)、θ_B(t)。

如果所分析的时间段内有多个目标出现，应尽量选取声学阵列正横方位附近的目标。

θ_A(t)、θ_B(t)应满足|θ_A(t)-θ_B(t)|≤2°。

S1-3，对两条声学阵列在目标方位的波束形成输出结果进行预处理。预处理后声学阵列A、B在目标方位的波束形成输出分别为

其中，k为波数，k＝2πf/c₀；L为声学阵列A、B几何中心水平间距；

为分析时间段内起始时刻所选目标相对于声学阵列B的水平方位；

N_t为总拍数。

优选的，N_t的典型取值范围为[600/(ΔT(1-P％))，1800/(ΔT(1-P％))]。ΔT为单拍数据的时间长度，单位s；P％为前后相邻两拍重叠率。

S1-4，对预处理后的波束数据进行互相关运算，并在一定时间内积累，求取其平均值，获得频域伪格林函数

(·)^*表示共轭。

S1-5，对频域伪格林函数进行反傅里叶变换，获得时域伪格林函数

其中，τ表示时延，ifft表示反傅里叶变换，Re表示取实部。

时域伪格林函数

是一个脉冲函数，其包络为v(τ)：

其中，|·|表示取绝对值，H(·)表示希尔伯特变换。

S2，根据声学阵列布放海域的水文环境，仿真计算时域伪格林函数，获得不同深度、距离上对应的时域伪格林函数脉冲包络v_r(τ；z_q,r_q)。

S2-1，在深度和距离上划分网格。将深度在[0,4H/5]范围内划分为Q_Z个网格，H为水深。将距离在[R_min,R_max]范围内划分为Q_R个网格，R_min为满足声学阵列远场要求所需的最小距离，R_max为声学阵列对目标的最大理论作用距离。

设网格总数为Q＝Q_R×Q_Z，Q_R为距离网格数，Q_Z为深度网格数。

优选的，Q_Z的典型取值范围为[8H/5λ,40H/5λ]。H为水深，λ为分析带宽内中心频率所对应的声波波长。

优选的，Q_R的典型取值范围为[(R_max-R_min)/10,(R_max-R_min)]。R_min、R_max的单位为km。

S2-2，根据水深、声速剖面、海底参数等海洋环境参数，采用声场计算程序(比如Kraken或KrakenC等)计算简正波信息：第m号简正波的水平波数k_m＝μ_m+jη_m，μ_m、η_m分别为k_m的实部、虚部；第m号简正波的特征函数Ψ_m(z)。其中m＝1,2，…，M。M为简正波总号数。

S2-3，根据简正波信息，计算各网格点的频域伪格林函数。

对于第q个网格点，其时域伪格林函数为

其中，τ表示时延，ifft表示反傅里叶变换，Re表示取实部。

其中，B(z₁,z₂,z_s)为归一化因子；符号“Σ”为求和算子；z_q、r_q分别为第q个网格点对应的声源深度和至两条声学阵列几何中心连线中点的水平距离；

为步骤S1中获得的目标方位。

时域伪格林函数

是一个脉冲函数，其包络为v_r(τ；z_q,r_q)：

其中，argmax(·)表示取最大值所对应的时延，|·|表示取绝对值，H(·)表示希尔伯特变换。

S3，将时域伪格林函数包络v(τ)与实际数据提取的时域伪格林函数包络v_r(τ；z_q,r_q)进行匹配处理，获得目标深度。

S3-1，计算每个网格点的时域伪格林函数包络v(τ)与实际数据提取的时域伪格林函数包络v_r(τ；z_q,r_q)的相关系数。

对于第q个网格点，其对应的相关系数为

其中，V(τ)和V_r(τ；z_q,r_q)分别是对v(τ)和v_r(τ；z_q,r_q)在时延维度上取滑动平均后的结果，

其中，Δτ为滑动平均窗口宽度。

优选的，Δτ的典型取值范围为[5/(f_H-f_L),20/(f_H-f_L)]。f_L、f_H分别为分析频段的下限频率、上限频率。

S3-2，对于所有的W(z_q,r_q)，q＝1,2,3，…Q，取局部最大值，其对应网格点的深度即为被测目标的深度估计结果。

由于上述算法对目标距离非常不敏感，因此局部最大值对应网格点的距离与目标真实距离可能有较大差异，因此不能作为目标距离的估计结果。

如果目标距离是已知的，记为r_s，那么仅对目标距离附近范围内的网格点进行局部最大值搜索即可。距离搜索的典型范围为[0.8r_s,1.2r_s]。

本发明根据浅海低频宽带声场的水平纵向互相关函数的物理特征，提出了一种基于互相关函数脉冲波形匹配的目标深度估计方法，适用于浅海海域的低频宽带被动目标深度估计，为提高浅海海域的目标深度分辨能力提供技术手段。

现通过附图和具体实施例来阐述本发明的技术方案。

浅海海域水深200m，海水平均声速剖面如图3所示，海面声速1544m/s,声速随深度逐渐减小,200m深度处声速约为1502m/s。海底声速1608m/s。图中实线(图中标记为ssp1)为计算测量场时所用声速剖面，虚线(图中标记为ssp2)为计算拷贝场时所用声速剖面。

两条声学阵列(声学阵列A、B)布放于海底表面，声学阵列A、B均为直线阵且相互平行，阵元数均为11，阵元间距均为5m。声学阵列A、B几何中心均为其第6阵元，几何中心连线长度为L＝700m，与两条声学阵列垂直，目标方位为θ＝0°。

目标距离在10分钟内由R＝30km均匀增加至R＝33km，目标方位为θ＝0°，目标深度为z_s＝10m或z_s＝100m，按照图3中声速剖面1计算声学阵列A、B所接收到的声场，作为测量场。声学阵列A、B各阵元的平均信噪比设为-5dB。然后，分别获得声学阵列A、B在目标方位的波束输出并进行归一化处理，得

和

在频域做互相关运算，并在一定时间段内取平均，得频域伪格林函数

然后对频域伪格林函数进行反傅里叶变换，获得时域伪格林函数

及其包络v(τ)。本实施例中，N_t取100，对应时间段长度为10分钟；频率分析范围为200Hz至300Hz。图4(a)显示了声源深度分别为z_s＝10m、z_s＝100m时的时域伪格林函数包络，横轴为时延，单位ms，纵轴为归一化强度。

将深度在[2.5m,200m]范围内划分为80个网格，步长2.5m。将距离在[10km,50km]范围内划分为41个网格，步长1km。根据图3中声速剖面2所示的海域环境参数，采用Kraken程序计算简正波信息，并据此计算各网格点的时域伪格林函数及其包络。对于第q个网格点，其时域伪格林函数包络为

其中，

为时域伪格林函数，

其中，τ表示时延，ifft表示反傅里叶变换，Re表示取实部。

其中，B(z₁,z₂,z_s)为归一化因子；z_q、r_q分别为第q个网格点对应的声源深度和至两条声学阵列几何中心连线中点的水平距离。

本实施例中，声学阵列A、B深度z₁＝z₂＝200m，声学阵列A、B几何中心的水平间距L＝700m，目标方位

当r_q取30km时，仿真计算得到的时域伪格林函数包络随时延和深度变化的二维图像如图4(b)所示，横轴为时延，单位ms，纵轴为深度，单位m。

然后，对v(τ)和v_r(τ；z_q,r_q)按照下式进行相关运行，获得每个网格点的相关系数

其中，V(τ)和V_r(τ；z_q,r_q)分别是对v(τ)和v_r(τ；z_q,r_q)在时延维度上取滑动平均后的结果，

其中，Δτ为滑动平均窗口宽度，本实施例中Δτ＝100ms。

当声源深度为z_s＝10m时，W(z_q,r_q)的二维图像如图5(a)所示，横轴为水平距离，单位km，纵轴为深度，单位m。当声源深度为z_s＝100m时，W(z_q,r_q)的二维图像如图5(a)所示，横轴为水平距离，单位km，纵轴为深度，单位m。

图5(a)所对应目标深度估计结果为8m，实际目标深度为10m；图5(b)所对应目标深度估计结果为101m，实际目标深度为100m。以上结果验证了本发明所示算法的有效性。

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

Claims (8)

1.一种基于互相关函数脉冲波形匹配的目标深度估计方法，所述方法包括：

步骤1)利用两条近似平行的声学阵列，从被测目标辐射噪声数据中提取时域伪格林函数，获得其脉冲包络；

步骤2)沿深度和距离方向划分网格，然后根据声学阵列布放海域的水文环境，计算各网格点上对应的时域伪格林函数，获取每个网格点的脉冲包络；

步骤3)将步骤1)提取的脉冲包络与各网格点上脉冲包络进行匹配处理，获得水中目标的深度估计结果；

所述步骤1)具体包括：

步骤1-1)对声学阵列A一段时间内的多拍数据进行处理，获得声学阵列A的宽带波束输出结果；获得某一目标相对于声学阵列A的随时间t变化的方位历程θ_A(t)；

步骤1-2)对声学阵列B一段时间内的多拍数据进行处理，获得声学阵列B的宽带波束输出结果；获得某一目标相对于声学阵列B的随时间t变化的方位历程θ_B(t)；

步骤1-3)若|θ_A(t)-θ_B(t)|≤α，α为阈值；转入步骤1-4)；否则改变时间段，转入步骤1-1)，直到在宽带波束输出结果中有满足|θ_A(t)-θ_B(t)|≤α的目标出现；

步骤1-4)对两条声学阵列在目标方位的波束形成输出结果进行预处理；对预处理后的波束数据进行互相关运算，并在一定时间内积累，求取其平均值，即为伪格林函数。

2.根据权利要求1所述的基于互相关函数脉冲波形匹配的目标深度估计方法，其特征在于，所述步骤1-1)具体包括：

步骤1-1-1)将声学阵列A各阵元接收到的声信号分成多拍，每拍时间长度为ΔT，前后相邻两拍重叠P％，并对每一拍信号做傅里叶变换得到频域信号；声学阵列A的频域信号为s_A(f；t)＝[s_A1(f；t),s_A2(f；t),...,s_AM(f；t)]^T；f是信号频率，t是当前拍数据的起始时刻；M是声学阵列A的阵元个数；(·)^T表示转置；

步骤1-1-2)以声学阵列A的几何中心为参考点，计算声学阵列A各频点的在波束形成扫描方位θ的加权向量w_A(f,θ)＝[w_A0(f,θ),w_A1(f,θ),...,w_AM-1(f,θ)]^T，其中

其中，τ_Am(θ)表示声学阵列A的第m+1个阵元相对于参考点的时延：

v(θ)＝-[cosθ,sinθ]^T

p_Am＝[p_xm,p_ym]^T

v(θ)是入射信号的单位向量，p_Am为第m+1个阵元相对于参考点的二维坐标，c₀为参考声速；

步骤1-1-3)计算声学阵列A指向方位θ的波束形成结果b_A(f,θ；t)：

b_A(f,θ；t)＝w_A ^H(f,θ)s_A(f；t)

其中，(·)^H表示共轭转置；θ的取值范围为[θ₀-45°，θ₀+45°]，θ₀表示声学阵列A和声学阵列B的参考阵元连线的角度；

步骤1-1-4)对声学阵列A的一段时间[t₁，t₂]内的多拍数据进行步骤1-1-1)至步骤1-1-3)的处理，获得声学阵列A的宽带波束输出结果H_A(θ,t)：

其中，符号“Σ”表示求和；

步骤1-5)根据声学阵列A的宽带波束输出结果，获得某一目标相对于声学阵列A的随时间t变化的方位历程θ_A(t)。

3.根据权利要求2所述的基于互相关函数脉冲波形匹配的目标深度估计方法，其特征在于，所述步骤1-2)具体包括：

步骤1-2-1)将声学阵列B各阵元接收到的声信号分成多拍，每拍时间长度为ΔT，前后相邻两拍重叠P％，并对每一拍信号做傅里叶变换得到频域信号；声学阵列B的频域信号为s_B(f；t)＝[s_B1(f；t),s_B2(f；t),...,s_BN(f；t)]^T；f是信号频率，选取的范围为[f_L,f_H]；t是当前拍数据的起始时刻；N是声学阵列B的阵元个数；

步骤1-2-2)以声学阵列B的几何中心为参考点，计算声学阵列B各频点的在波束形成扫描方位θ的加权向量w_B(f,θ)＝[w_B0(f,θ),w_B1(f,θ),...,w_B,N-1(f,θ)]^T，其中

其中τ_Bn(θ)表示声学阵列B的第n+1个阵元相对于参考点的时延：

v(θ)＝-[cosθ,sinθ]^T

p_Bn＝[p_xn,p_yn]^T

v(θ)是入射信号的单位向量，p_Bn为第n+1个阵元相对于参考点的二维坐标；

步骤1-2-3)计算声学阵列B指向方位θ的波束形成结果b_B(f,θ；t)：

b_B(f,θ；t)＝w_B ^H(f,θ)s_B(f；t)

其中，(·)^H表示共轭转置；

步骤1-2-4)对声学阵列B的一段时间[t₁，t₂]内的多拍数据进行步骤1-2-1)至步骤1-2-3)的处理，计算声学阵列B的宽带波束输出结果H_B(θ,t)：

步骤1-2-5)根据声学阵列B的宽带波束输出结果获得同一目标相对于声学阵列B的随时间t变化的方位历程θ_B(t)。

4.根据权利要求3所述的基于互相关函数脉冲波形匹配的目标深度估计方法，其特征在于，所述步骤1-4)包括：

步骤1-4-1)对两条声学阵列在目标方位的波束形成输出结果进行预处理，预处理后声学阵列A、B在目标方位的波束形成输出

分别为：

其中，k为波数，k＝2πf/c₀；L为声学阵列A、B几何中心水平间距；

为分析时间段内起始时刻所选目标相对于声学阵列B的水平方位；

I_A(f)为所选取目标相对于声学阵列A的平均声强谱：

I_B(f)为所选取目标相对于声学阵列B的平均声强谱：

其中，N_t为总拍数；N_t的取值范围为[600/(ΔT(1-P％))，1800/(ΔT(1-P％))]；ΔT为单拍数据的时间长度；P％为前后相邻两拍重叠率；

步骤1-4-2)对预处理后的波束数据进行互相关运算，并在一定时间内积累，求取其平均值，获得频域伪格林函数

(·)^*表示共轭；

步骤1-4-3)对频域伪格林函数

做反傅里叶变换，得到时域伪格林函数

其中，ifft表示反傅里叶变换，Re表示取实部；

步骤1-4-4)时域伪格林函数

是一个脉冲函数，其包络为v(τ)：

其中，|·|表示取绝对值，H(·)表示希尔伯特变换。

5.根据权利要求4所述的基于互相关函数脉冲波形匹配的目标深度估计方法，其特征在于，所述步骤2)具体包括：

步骤2-1)按照[0,4H/5]间隔，将声学阵列布放海域在深度上划分为Q_z个网格，H为水深；按照[R_min,R_max]间隔，将声学阵列布放海域在距离上划分为Q_R个网格，R_min为满足声学阵列远场要求所需的最小距离，R_max为声学阵列对目标的最大理论探测距离；

步骤2-2)根据水深、声速剖面和海底参数这些海洋环境参数，采用声场计算方法计算简正波信息：第l号简正波的水平波数k_l＝μ_l+iη_l，μ_l、η_l分别为k_l的实部、虚部；第l号简正波的特征函数Ψ_l(z)；其中，l＝1,2，…，I；I为简正波总号数；

步骤2-3)根据简正波信息计算各网格点的时域伪格林函数；

对于第q个网格点，其时域伪格林函数为：

其中，z_q、r_q分别为第q个网格点对应的声源深度和至两条声学阵列几何中心连线中点的水平距离；τ表示时延，ifft表示反傅里叶变换，Re表示取实部；

为第q个网格点的频域伪格林函数：

其中，B(z₁,z₂,z_q)为归一化因子：

z₁,z₂表示两个水听器的深度，其水平间距为L，两个水听器与声源的水平间距分别为r₁,r₂；

步骤2-4)时域伪格林函数

是一个脉冲函数，其包络为v_r(τ；z_q,r_q)：

其中，|·|表示取绝对值，H(·)表示希尔伯特变换。

6.根据权利要求5所述的基于互相关函数脉冲波形匹配的目标深度估计方法，其特征在于，深度网格数Q_z的取值范围为[8H/5λ,40H/5λ]，λ为分析带宽内中心频率所对应的声波波长；距离网格数Q_R的取值范围为[(R_max-R_min)/10,(R_max-R_min)]。

7.根据权利要求6所述的基于互相关函数脉冲波形匹配的目标深度估计方法，其特征在于，所述步骤3)具体包括：

步骤3-1)计算第q个网格点的时域伪格林函数包络v(τ)与实际数据提取的时域伪格林函数包络v_r(τ；z_q,r_q)的相关系数W(z_q,r_q)：

其中，V(τ)和V_r(τ；z_q,r_q)分别是对v(τ)和v_r(τ；z_q,r_q)在时延维度上取滑动平均后的结果：

其中，Δτ为滑动平均窗口宽度；Δτ的取值范围为[5/(f_H-f_L)，20/(f_H-f_L)]；f_L、f_H分别为分析频段的下限频率、上限频率；

步骤3-2)对于所有的W(z_q,r_q)，其最大值对应网格点的深度即为被测目标的深度估计结果。

8.根据权利要求7所述的基于互相关函数脉冲波形匹配的目标深度估计方法，其特征在于，如果目标距离是已知的，记为r_s，那么仅对目标距离附近范围内的网格点进行局部最大值搜索即可；距离搜索的典型范围为[0.8r_s,1.2r_s]。

Images