基于首达波与次达波到达时差的深海浅层移动声源定位

高飞， 孙磊,2

(1.海军研究院 海洋环境研究所, 天津 300061； 2.哈尔滨工程大学 水声工程学院, 黑龙江 哈尔滨 150001)

0 引言

水下移动目标定位一直以来都是水声学研究的热点问题，广泛应用于海洋渔业、水下搜救、水声制导等领域，极具军事和民用价值。传统的水声定位通常使用多个水听器阵元组合成一定阵型，例如垂线阵、水平阵、L形阵、球面阵等，基于阵列信号采用Bayes匹配场算法处理[1-3]，可达到良好的目标定位效果。

单水听器声源定位主要利用宽带声信号的频散特征及多途声波的声到达时差，可有效节约观测成本，并提高水下探测效率。对于水下宽带声信号，基于声波能量传输与群速度的相关性，可利用简正波理论模型计算得到宽带信号的模内离散、模间离散信息，并依据各频率分量时域上的差异构建代价函数，进行声源定位及声传播环境参数的反演[4-6]。

多途波是指声源单次辐射信号沿不同传播路径到达同一接收点的声波，由于传播路径总长度的差异，通常可从时域对其进行区分。Aubauer等[7]通过声线几何关系推导，在不考虑声速垂向变化的前提下，利用直达波与海底反射波的声到达时差对浅海鱼群进行定位；Nosal等[8]则通过坐底水听器数据得到直达波、海面反射波声到达时差，融合实测声速剖面反演浅海鱼群位置；Gebbie等[9]和Duan等[10]则采用连续声源信号，基于直达波、海面反射波声到达时差反演移动声源位置，并验证了定位结果的准确性。

事实上，海洋环境的不确定性、信号接收的偶然性及数据处理的方法差异都会影响多途声到达时差的计算结果，直接关系到声源定位的精度。本文采用扩展卡尔曼滤波(EKF)算法进行移动声源定位，以首达波与次达波声到达时差为观测值进行连续滤波，可有效克服观测值与实际值间的偏差对定位结果的影响，得到稳定、高精度的声源定位参数。

1 定位算法理论

构建水平不变深海多途声传播环境如图1所示。移动声源(深度zs，水平距离rs)和水听器(深度zr，水平距离0)位于浅层水域，其深度远小于水深H. 图1中GPS、SOFAR分别为差分全球定位系统、深海声道。

分别设直达波、海面反射波、海底反射波、海面海底反射波传播路径长度分别为mD、mSR、mBR、mSBR，传播时间分别为τD、τSR、τBR、τSBR，则有

(1)

(2)

rs通常取104m以上量级，对于深海浅层收发分置的水声设备，zs、zr深度值约为102m量级，故可得mD≈mSR≈rs. 当目标信号脉冲长度较大时(约500 ms)，直达波与海面反射波重合，难以从时域对其进行分离，故下文将直达波和海面反射波统称为首达波。

(3)

(4)

值得一提的是，在实际海洋环境中，受声速垂向变化的影响，声波并非沿直线传播，上文声线路径长度是一种近似推导过程。

1.1 多途波识别

受海洋环境噪声影响，通常难以直接识别原始水听器时间序列数据(见图2(a))中的目标信号。对于频率已知的拖曳式移动声源，利用功率谱计算或带通滤波处理原始信号，可分别从频域和时域识别目标信号。

选取水听器与移动声源相对距离为10.5 km的一段水听器原始数据(见图2(a))，移动声源间隔10 s发射一次中心频率为500 Hz、脉冲长度为500 ms的正弦波激励信号，分段计算其功率谱，其中频率为400～600 Hz的功率谱计算结果如图2(b)所示。从中不难发现，声波多路径到达特征明显，图2(b)中方框、圆圈、虚线方框分别标注不同路径到达波。

为确认多途到达波的具体传播路径及相对能量大小，通过Butterworth带通滤波处理时域信号(见图2(c))。对比电压幅度可知，到达声波C、A、B能量依次增加。采用Bellhop射线模型(模型代码下载地址：http:∥oalib.saic.com)模拟试验数据海洋环境中的声线分布特征(见图2(d))，前3个到达波束分别为首达波、次达波、第2次海底反射波，其相对到达时间与声波A、B、C吻合较好，故可判定声波A、B、C分别为首达波、次达波、第2次海底反射波，且次达波能量最大，首达波次之，第2次海底反射波最小，故下文选择直达波与第1次海底反射波的到达时差进行声源定位。详尽多途分类及各到达波的能量相对大小可参照文献[11]。

1.2 多途波时延算法

相关性函数用于描述随机过程中不同时段或空间内的相互关系，对于连续水声信号时间序列，时域相关性函数[12-13]可表示为

(5)

式中：相关系数Γ为时延Δτ的函数；τU、τL分别为分段数据的时间上限、下限(s)；pm、pe为匹配、实测声压时域波形数据。通过pm=sin (2πfτ)计算匹配声压，其中f=500 Hz为频率，τ=τL∶1/fs∶τU，fs=20 000 Hz为采样频率。

本文取匹配声压长度为0.5 s(与移动声源激励信号长度相同)，即τU=τL+0.5，沿原始采样时域数据逐个采样点平移，并获取各组最大相关系数。以计算得到的与首达波、次达波的最大相关系数分别确定其相对到达时间τF、τS，进而计算声到达时差ΔτF-S(见图3)。

1.3 EKF算法

EKF算法通过建立非线性系统模型，主要为描述过程变化动态的状态方程和系统参量与测量值间关系的测量方程，以实现其最小方差条件下的递推，从而达到多维平稳、非平稳随机过程的最优估计。对于速度为v(m/s)的移动声源，其状态方程可表述为

xk|k-1=Φk-1|k-1+qk-1，

(6)

式中：状态矢量x=[rs,vr,zs,vz]T，vr、vz分别为水平、垂直方向的声源移动速度；xk|k-1为k-1时刻得到的k时刻状态预测结果；k-1|k-1为k-1时刻的状态最优估计值；qk-1～N(0，αk-1)为过程噪声，αk-1为过程噪声均方差；状态转移矩阵Φ可表示为

(7)

Δt为迭代时间步长(s)。在水平距离不变的深海环境中，k时刻的首达波与次达波声到达时差ΔτF-S可表示为声速剖面c(z)、水听器位置以及H、rs、zs的函数，假设c(z)、H、水听器位置已知，则测量方程可表示为

ΔτF-S,k=h(rs,k|k-1,zs,k|k-1)+pk,

(8)

式中：pk～N(0，βk)为测量噪声，βk为测量噪声均方差；函数h(·)为声线传播模型，无具体表达形式。于是测量方程可用Taylor级数展开，并取1阶近似：

ΔτF-S,k≈h(s,k-1|k-1,s,k-1|k-1)+s,k-1|k-1)+s,k-1|k-1),

(9)

式中：上标“^”表示已得到的参数估计值。利用Bellhop模型计算得到的首达波与次达波声到达时差如图4所示。

图4中水听器深度位于30 m，声速剖面如图5所示，海表声速混合层厚度约45 m，45～200 m水层为强跃层。从中不难发现ΔτF-S随声源深度、水平距离呈负梯度变化特征，即(9)式中的∂h/∂zs、∂h/∂rs均小于0，说明随着移动声源深度、距离的增加，ΔτF-S逐渐减小。

依据(9)式，测量方程的雅克比矩阵为

(10)

(10)式中的∂h/∂rs、∂h/∂zs虽无具体表达式，但可依据Bellhop射线模型计算各空间点处的声到达时差(见图4)，进而得到声到达时差随声源距离、深度变化的梯度值。

按照非线性离散系统下的EKF进行状态更新迭代：

Pk|k-1=Φk-1|k-1ΦT+αk-1，

(11)

(12)

k|k=xk|k-1+Kk(ΔτF-S,k-Hkxk|k-1)，

(13)

k|k=Pk|k-1-KkHkPk|k-1，

(14)

式中：Pk|k-1为基于k-1时刻预测的k时刻协方差矩阵；k|k为k时刻的协方差最优估计值矩阵；Kk为k时刻的EKF增益矩阵。(11)式、(12)式、(13)式、(14)式分别为预测协方差方程、EKF增益方程、状态更新方程、协方差更新方程。EKF声源定位的算法具体步骤如下：

1)k=0，状态矢量0|0、协方差0|0初始化；

2)基于(6)式，利用k-1时刻状态k-1|k-1预测k时刻的状态xk|k-1；

3)基于(11)式，计算预测协方差矩阵Pk|k-1；

4)依据预测状态xk|k-1和Bellhop射线模型计算声到达时差梯度Hk，进而基于(12)式，计算EKF增益Kk；

5)基于(13)式，输入声到达时差ΔτF-S，得到更新状态k|k；

7)单次迭代计算结束，返回第2步进行k+1时刻的EKF.

2 声源定位与结果验证

试验海域水深H≈4 150 m，移动声源间隔10 s发射中心频率f=1 000 Hz的正弦波，声速剖面如图5所示，利用位于30 m水深的水听器数据计算首达波与次达波声到达时差(见图6)。

利用1.2节中的时域相关性算法计算得到首达波与次达波声到达时差，所用数据时长为30 min，ΔτF-S的测量值分布于1.4～1.8 s之间，且随时间呈减小的趋势。根据图4，可初步推测移动声源的水平距离逐渐增加，或者所处深度逐步增大。

本文中深海水声调查区域，海表声速混合层约45 m，受海洋声速剖面影响，混合层以下直达波迅速减弱消失；海表至约20 m水层受散射衰减影响，直达波幅度较弱，故选择30 m水层数据。

采用压力传感器采集记录移动声源实时深度，移动声源与水听器间的水平距离通过GPS信号计算得到，水平移动速度通过测量船航速共享可得。由于无移动声源垂向速度实测数据，这里不对vz进行分析。水声试验过程中，连接移动声源的拖缆放缆长度保持不变，声源深度约49.1 m，且起伏较小。移动声源以3.40 m/s沿直线近似匀速离开水听器垂线阵。

分析可知，EKF得到的声源位置与实测数据总体吻合较好。受海洋环境及试验过程中的不确定性影响，试验获取的首达波与次达波声到达时差随距离和深度呈现严格单调的变化趋势(见图6)，利用(5)式计算的ΔτF-S具有一定误差，这也是导致滤波结果起伏的重要因素。

计算得到EKF迭代稳定后的声源深度(见图7(b))和速度(见图7(c))的平均值分别为48.8 m、3.42 m/s，实测数据的均值分别为49.1 m、3.38 m/s，二者差别较小，验证了本文声源定位结果的准确性。

3 结论

1)深海浅层水域是海洋鱼群、水下航行器等活动频繁的水层，也是声纳等水下探测设备经常布放的深度。本文利用深海水声调查数据，并依据其建立了声源定位模型，采用时域相关性分析计算多途声到达时差，以声源深度、水平距离、速度为变量参数构建了EKF状态方程((6)式)，并结合Bellhop射线模型得到EKF测量方程((8)式)，考虑声速剖面，进而进行了移动声源定位。

2)通过将EKF定位结果与压力传感器、GPS及母船航速对比，结果表明迭代稳定后的定位结果与实测数据吻合较好，可验证本文算法的可靠性。

3)准确的识别多途信号中的直达波与第1次海底反射波并计算其声到达时差，是影响EKF算法精度的重要因素，需要考虑深海声速剖面、波导截止频率对浅表层直达波的影响。