单水听器浅海声源被动定位方法研究

石海杰,李京华,常虹

(1．西北工业大学电子信息学院， 陕西西安 710072； 2．中北大学信息与通信工程学院， 山西太原 030051；3．西安邮电大学通信与信息工程学院， 陕西西安 710121)

自冷战结束以后,对海洋的关注更多转向近岸浅海,由于干涉结构和环境参数的非稳定性,很难直接利用接收的声信号进行目标定位,需要充分考虑声场模型和声源的先验信息。匹配场定位是典型的基于模型的方法,自诞生以来相关方法不断涌现[1-4],但该方法面临的最大挑战是模型失配问题,特别是在浅海环境下,建立确定意义上的声传播模型通常是很困难的[5]。同时,为了获取更多空间信息,匹配场定位往往需要采用阵列接收信号,而在隐蔽低耗平台条件下对运动目标定位时,需要的是一种利用单个水听器就能实现定位的方法。李晓曼等[6]利用warping变换分离宽带脉冲信号的简正波,根据群延迟理论实现单个水听器条件下的目标定位,但实际上,窄带信号具有更好的传输特性。谢亮等[7]分析典型深海环境中脉冲信号的到达结构特征,提出通过特征参数匹配搜索进行目标定位的方法,单个水听器海上实验的距离估计误差不大于8%。Zhang等[8]提出一种基于自相关函数峰值提取的无源定位方法,利用单个水听器就能实现定位,但该方法无法适应非稳定环境。Warner等[9]利用单个水听器，采用基于模型的贝叶斯估计方法实现了水声参数和海底参数估计,但并未涉及定位方法的研究。贝叶斯估计理论是一种基于概率密度模型的方法,参数估计时以后验概率的形式给出,由于可以通过时间累计效应来拓展空间信息的缺失,因此可以采用单个水听器对目标进行定位。Ogiso等[10]利用递归贝叶斯滤波方法在室内环境条件下,实现移动机器人自主定位,该方法属于高精度定位,但并不适用于水中被动定位情况。王彪等[11]将贝叶斯理论应用到多目标方向估计问题当中,利用多矢量稀疏贝叶斯学习算法重构信号空间谱,仿真结果表明该方法相对于传统的算法具有更高的空间分辨率和估计精度。Ren等[12]采用贝叶斯滤波方法,利用矢量传感器测得的船舶噪声数据跟踪海底环境参数,该方法可以有效解决环境参数非平稳变化的问题。Chu等[13]提出一种基于非同步阵列测量的贝叶斯推断方法细化点阵函数,实验结果表明该方法即使在拉普拉斯噪声干扰下也有较高的定位精度。Guo等[14]将贝叶斯方法用于冰下信道估计,利用估计的信道对脉冲噪声进行迭代估计和消除,仿真结果表明,所提出的迭代算法具有较高的估计精度。Sun等[15]将贝叶斯技术应用于速度跟踪,采用统计推断方法对未知非平稳噪声的统计特性进行动态跟踪,提出的算法可以有效地解决卡尔曼滤波对初始测量噪声方差过于依赖的问题。文献[16-18]应用贝叶斯理论对环境聚焦以期解决匹配场定位时的失配问题,其实质是将未知环境参数和声源位置进行联合优化估计,该类方法解决了环境失配问题却增加了算法复杂程度。总之,在浅海环境条件下,声信号在传播过程中受到环境的影响尤为严重,多种状态的海面和海底参数的相互作用导致了具有复杂空间干涉结构的声场。而当低频线谱声源移动时,声场的空间复杂调制特性转变成固定点接收信号的时间幅度调制特性。

本文提出的基于贝叶斯估计的目标定位方法是一种基于模型的方法,不直接利用声场环境参数进行匹配,也不聚焦估计声场环境参数,而是利用声场概率密度模型消除非稳定因素,参数估计结果采用后验概率的形式给出,能够有效克服随机环境对定位结果的影响,而且能够通过时间累计效应拓展空间信息的缺失,实现单个水听器对目标进行定位,因此基于贝叶斯估计水声目标定位方法是本文主要研究内容。

1 随机声场统计模型

为了解决非稳定性环境参数导致的模型失配问题,一种有效的方法是利用信道的统计特性,以概率密度函数描述声场传输模型。概率密度函数估计有参数化方法和非参数化方法,核密度估计(kernel density estimation,KDE)是典型的非参数化方法。KDE法不事先设定随机过程的分布规律,而是从样本本身出发研究数据分布特征,所以理论上KDE法适用于任何随机过程概率密度函数的求解。

(1)

式中：κ(x-xi;w),i=1,2,…,n为核函数,w为核函数宽度；n为核函数个数。核函数可以是均匀核函数、三角核函数、高斯核函数等。理论上任何平滑形式的峰值函数均可作为KDE的核函数来使用,只要满足函数曲线下方的积分面积等于1即可。根据(1)式可知,估计结果实际上是多个样本数据点的核函数曲线的叠加,由于邻近样本点之间会发生波形叠加,因此最终所形成的曲线形状与选择的核函数关系并不密切(类比任意波形信号都可以很好地由正弦波、方波或者三角波正交基合成)。考虑到函数在波形合成计算上的易用性以及高斯分布的普遍适用性,选择使用高斯曲线作为KDE的核函数。高斯核函数可表示为

(2)

每个核函数都是以样本点xi为均值,以w为方差的高斯分布。核函数宽度w能够控制结果曲线的整体平坦程度,也就是样本点数据在概率密度曲线形成过程中所占比重。w越宽,样本点在最终形成的曲线形状中作用范围越大,曲线越平坦。反之w越窄,样本点在最终形成的曲线形状中作用范围越小,结果曲线越陡峭。实际上,这也说明估计得到的概率密度是许多样本点分布直方图平滑的结果,这就是核密度估计方法确定概率密度函数的核心思想。

核函数宽度w按(3)式确定[19]

(3)

式中,σx为随机变量x的标准差,通过(4)式确定

(4)

(5)

2 贝叶斯估计定位方法建模与求解

2.1 浅海声场信号传输几何模型

如图1所示,声源S以径向速v[n]运动,n为离散时间,z[n],r[n]分别为声源深度和距离,Q为静止不动的水听器。声源发射的声信号为

a[n]=As[n]cos(2πf0n+φ0)

(6)

式中,As[n]为声源声压幅度。声源状态矢量s[n]包含4个分量:r[n],v[n],z[n],As[n],目标定位通过估计r[n],z[n]实现。水听器接收的信号为

y[n]=A(s[n])cos(2πfdn+φo)+w[n]

(7)

其中,A(s[n])为接收信号声压幅度。传输衰减用x[n]表示,则对于已知固定位置的水听器,x[n]由(8)式决定

A(s[n])=As[n]x[n]

(8)

图1 声场信号传输几何模型

2.2 声源状态矢量贝叶斯估计建模

基于贝叶斯估计的水声定位方法,就是在已知水听器位置的条件下,利用声场的概率密度模型和接收信号y[n],确定声源状态矢量s[n]。根据声场概率密度模型可以确定以声源状态矢量为条件的概率密度f(x[n]|s[n]),也称先验概率密度,再结合接收信号y[n]提取的信号幅度A(s[n]),利用贝叶斯估计方法求解后验概率,从而估计出声源位置参数。若以后验概率的最大化函数作为参数估计的求解条件,则声源状态矢量估计值为

(9)

根据贝叶斯公式,可得

f(s[n]|x[1],…,x[n])=

(10)

根据条件概率公式和乘法公式,(10)式变形为

(11)

(11)式中各函数含义如下:

1)f(s[n]|x[1],…,x[n])是要求解的后验概率,表示当声场中某一固定点接收到信号为x[1],…,x[n]时,声源状态矢量s[n]的概率密度。

2)f(x[n]|s[n],x[1],…,x[n-1])表示声源状态矢量为s[n],在声场中某个固定点已经接收到x[1],…,x[n-1]条件下,再次接收到x[n]的概率密度。声场中某点接收的信号由声源状态矢量和传输模型决定,与该点之前接收的信号无关,因此

f(x[n]|s[n],x[1],…,x[n-1])=f(x[n]|s[n])

(12)

3)f(s[n]|x[1],…,x[n-1])表示当声场中某一固定点接收到信号为x[1],…,x[n-1]时,声源状态矢量s[n]的概率密度,实际上就是前一时刻的后验概率。

4)f(x[n]|x[1],…,x[n-1])表示不考虑声源状态情况下,声场中某点已经接收到x[1],…,x[n-1]条件下,再次接收到x[n]的概率密度。由于该函数与声源状态无关,因此可以假设

(13)

式中,η为常数,满足后验概率对状态矢量的积分为1。

将(12)～(13)式代入(11)式,整理得

(14)

用全概率公式展开f(s[n]|x[1],…,x[n-1]),得:

(15)

式中,f(s[n]|s[n-1])=f(s[n]|s[n-1],x[1],…,x[n-1])是状态更新概率密度,状态更新过程有马尔科夫性。

至此,得到了递归形式的后验概率模型,如(15)式所示,也是基于贝叶斯估计的水声定位方法的模型,也称为贝叶斯滤波公式,由于积分运算的存在,递归过程很难求得封闭形式解,需要研究有效的求解方法。

2.3 直方图滤波法求解贝叶斯估计问题

直方图滤波是贝叶斯滤波公式的离散化实现方法,将声源状态矢量离散化成多维网格,认为单一固定网格内状态矢量的概率密度为常数。相当于利用离散的概率直方图代替连续的状态矢量概率密度函数,因此称为直方图滤波。

假设将声源状态各分量离散成Kr,Kv,Kz和KA点,声源状态矢量总网格数为K=KrKvKzKA。设第k个网格的状态矢量为sk,则贝叶斯滤波公式(15)变成离散形式

(16)

(16)式中各函数含义如下:

1)P(sk[n]|x[1],…,x[n])表示当声场中某一固定点接收到信号x[1],…,x[n]时,声源状态矢量为s[n]=sk的概率。

2)P(x[n]|sk[n])表示在n时刻,声源状态矢量为sk条件下,声场中对应点处接收到x[n]的概率。

3)P(sk[n]|sj[n-1])表示状态矢量转移概率,在整个离散状态矢量空间内,k和j的取值都可以是1～K,因此状态转移概率有K2个可能取值。

运动目标的状态更新过程可用(17)式表示

s[n]=Ans[n-1]+εn

(17)

式中,An为状态转移矩阵

(18)

对于离散的状态矢量,状态转移概率为

(19)

4)P(sj[n-1]|x[1],…,x[n-1])表示前一时刻声源状态矢量的后验概率,是前一时刻迭代的输出,也是后一时刻迭代的输入。η1是一个常数,满足后验概率对状态矢量的求和为1。

后验概率公式递推至n=1时刻,需要声源状态矢量的初始分布作为迭代公式的输入,均匀分布可以很好地表示状态矢量的初始分布,因此

(20)

以上推导过程,将贝叶斯估计问题进行了离散化,(16)式就是贝叶斯滤波问题的数值解,该方法解决了状态矢量积分难于求解的问题。

2.4 分级网格划分直方图滤波法

观察直方图滤波算法,不难发现,网格的划分方法是影响直方图滤波效果的关键,网格划分过大,估计精度必然不高,网格划分过小,算法运算量大增,针对此问题,提出一种分级网格划分的直方图滤波方法,该方法在相同网格划分精度情况下,可以大大降低算法的运算量。

图2 网格划分方法对比图

以Kr=16为例,单级网格划分方法需要16次后验概率估计运算和16次搜索运算,而采用二级网格划分方法则只需要8次后验概率估计运算和8次搜索运算,在相同的预期估计精度情况下,改进方法效率提高一倍。以Kr=100为例,单级网格划分方法需要100次后验概率估计运算和100次搜索运算,而采用二级网格划分方法则只需要20次后验概率估计运算和20次搜索运算,在相同的预期估计精度情况下,改进方法效率提高5倍。随着网格划分数的增加,效率提高的效果将会更加显著。

3 基于SWellEx-96数据的算法验证

3.1 SWellEx-96实验说明

SWellEx-96系列实验是1996年5月在加利福尼亚圣地亚哥海岸附近进行的,其中的S5实验条件为：拖曳深度为9 m的声源以2.5 m/s速度运行,采用海底布放的水平阵列接收(south horizontal line array,SHLA)。由于实验信息记录完整,GPS数据准确,实验数据质量较高,因此非常适合于浅海环境条件下目标定位方法的验证。实验海域海深约200 m,SHLA布放深度198 m,本文利用SHLA中的3个阵元接收的数据进行验证,声源含多个频率的线谱,本文利用最低频率109 Hz信号,详细的环境参数可参考http:∥swellex96.ucsd.edu/index.htm。

3.2 核密度估计法随机声场建模结果

按SWellEx-96实验提供的51组声速剖面数据,构成一组输入参数,假设海水深度服从范围180～200 m、间隔为1 m的随机均匀分布,由21个深度组成另外一组输入参数,2组参数任意组合,构成包含1 071个样本的随机环境参数,利用Kraken程序仿真1 071个对应的声场分布样本,KDE方法根据这些声场样本可以估计声场概率密度模型。图3所示为声场中一些点的概率密度函数结果展示,其中包含海面近场、靠近声源(海底近场)、海面远场和海底远场的估计结果,这些是典型的特殊位置,也包括典型的一般位置的估计结果。从样本的绿色柱状直方图对比可以看出:①不同位置的直方图是不同的,同一区域的直方图也不同,特殊位置与一般位置的直方图都是随机的,没有观察到可遵循的一般规律;②近场、特别是靠近声源位置的直方图更不规律,更容易出现多峰情况,说明这些位置的声场信号的随机性更强;③在仿真的过程中发现临近位置的直方图区别可能也很大,这是由声场干涉结构决定的,干涉结构的纹理越细,声场随机性越强,这也再次说明只有基于统计特性的信号处理方法才能适用于此处。

图3中的红色划线表示核密度方法估计的概率密度函数,与高斯模型和最大熵估计得到的概率密度模型对比可以看出:高斯模型法得出的概率密度函数始终保持高斯函数固有的钟形脉冲形式,只在样本分布为单峰情况下勉强适用;最大熵法好于高斯模型法,但对多峰分布的适应能力也不强;核函数方法估计出的概率密度函数在不同情况下都有很好的适应性,即使在靠近声源位置,直方图分布具有多个峰值的情况下,估计结果与样本分布也匹配得很好,所以,本文采用核密度方法估计声场概率密度函数f(x)。将声源状态矢量参数带入f(x),就可以求得先验概率P(x[n]|s[n])。

图3 声场概率密度函数估计结果

3.3 SWellEx-96实验数据预处理

为了验证单个水听器能够实现定位,以阵元1、阵元17和阵元32作为3个独立水听器,提取它们的数据进行前置滤波(抗混叠)、降采样、混频、滤波处理,将最终提取的经过传输调制的幅度值,作为A(s[n])的估计值,用以为后验概率估计时备用。图4所示为预处理结果图,原始数据有50 min时长,选取了5 min进行展示和预处理,可以看出1号阵元信号幅度最小,32号阵元信号幅度最大,这种规律是由阵元位置关系决定的。原始数据频域图中能够很清晰地看到有109 Hz和112 Hz信号的存在,由于声源是移动目标,因此线谱被展宽了。处理后时域、频域信号可以看出109 Hz信号被有效地提取出来。

图4 实测数据预处理结果

3.4 基于贝叶斯估计的水声定位结果

声源状态矢量按表1取值,确定各声源状态矢量的范围和网格划分数K,声源幅度为归一化范围,因此无单位。信号预处理时将实测数据采样率降为1 000 Hz,因此T=0.001 s,状态更新噪声参数按照表2取值。根据(19)和(20)式可以确定状态转移概率和初始状态分布。单次定位取20个数据进行迭代,每隔1 min定位1次。

表1 声源状态矢量参数表

表2 状态更新参数表

定位结果如图5和图6所示,图5是三阵元分别在3个时刻的定位结果,图中蓝色圆圈中心位置代表声源的真实位置,红色十字中心位置是该时刻的定位结果,灰度图表示该时刻目标在定位区域内后验概率分布。图6是定位结果与真实值随时间变化的对比情况。综合观察图5和图6的定位结果图可以得出如下结论:①表示概率分布的灰度图是以一个具有较大值区域(灰度较暗的区域)的方式跟进目标真实位置,最大概率值点与该区域有较高的一致性,而非某些突兀的概率奇点作为定位结果,这能说明本文方法是一种具有较高可靠性的方法;②开始时刻的定位精度没有后续时刻定位精度高,这是由于开始时刻声源状态矢量采用均匀分布方式表示,后续时刻声源状态矢量是前一时刻的后验概率,随着时间的推进,定位结果将得到不断纠正,越来越接近真实目标位置;③距离定位结果的稳定性好于深度定位结果,这是由于声场模型干涉结构在距离方向的尺度远大于深度方向的尺度,尺度越小干涉越复杂,定位越困难,稳定性越差。实际上,浅海环境目标的深度定位一直是一个难以解决的问题。从实测数据的定位结果来看,本文方法的可取之处在于距离方向定位的稳定性,更在于能够较准确估计目标的深度。

图5 各阵元不同时刻定位结果

图6 定位结果随时间变化情况

3.5 分级网格划分直方图滤波法效率仿真

图7所示为单级网格划分与分级网格划分实测效率对比图,仿真时以普通办公用计算机为硬件平台,以仿真工具Matlab为软件平台。网格化参数K取值:[1,25,100,400,2 500,4 900,6 400,10 000],迭代次数为50,图中的横坐标为K值取对数,纵坐标为相对时间取对数,ΔT为K值取1时单次迭代的时间开销。观察图7的效率对比图可以得出如下结论:①分级网格划分直方图滤波法对算法效率改善是非常明显的,网格化参数K取值越大,效率提高效果越明显,这一点与理论分析是一致的;②当网格化参数K取值较小时(小于25次),算法总的开销时间还会受到其他因素 (包括临时变量初始化、概率密度函数装载、后验概率空间分配等)影响,随着K值的增加,算法开销时间主要决定于后验概率的迭代求解过程,K值越大,这种决定性越强。

图7 分级网格划分法效率对比

3.6 基于贝叶斯估计目标定位方法的精度

定义(21)式为距离均方根误差,定义(22)式为深度均方根误差

(21)

(22)

可以得到表3的结果,表中均值是指3个阵元均方根误差的平均值。说明本文研究的定位方法,深度定位精度可达到35 m,距离定位精度可以达到0.69 km。考虑做均方根误差的样本包含定位开始阶段偏差较大的定位结果,因此实际定位效果要更好。r[n]范围为10 km,直观估计目标距离估计误差在7%以内。与文献[18](6%)、文献[7](8%)、文献[6](10%)方法相比,本文方法距离定位结果总体略好,并且效率高出许多。由于合适的深度定位精度对比的文献不多,所以没有进行深度定位精度的横向比较。

表3 定位结果均方根误差表

4 结 论

贝叶斯滤波原理的实质是试图用所有已知信息来构造系统状态变量的后验概率密度,即用系统模型预测状态的先验概率密度,再使用最近的量测值进行修正,得到后验概率密度。这样通过观测数据的不断更新来递推计算状态参数,由此获得状态的最优估计。本文将贝叶斯估计理论应用于水声信号的跟踪与定位算法中,根据浅海条件和单个水听器接收条件,利用声场环境先验信息和水听器接收的声场数据,通过统计算法确定声源状态矢量后验概率,充分利用环境的不确定性和水声传输模型的先验知识较好地实现了声源目标被动跟踪和定位。