基于贝叶斯点估计的次声事件定位方法

王志峰，温宏波，杨立学

（1.中国电子科技集团公司第三研究所，北京 100015；2.陆军装备部驻北京地区军事代表局，北京 100015）

本文采用基于概率的方法，将次声传播速度视为未知参数，同时将声源位置坐标、时间发生时间视为未知参数，通过计算各次声分别测得目标方位和到达时间的条件下这些未知参数的最大后验概率，可得到目标声源的准确位置，本质上属于基于概率的定位方法［1］。

1 基本原理

在贝叶斯统计中，p(x|ω)表示在未知参数ω给定的条件下总体的概率密度，其中未知参数ω有先验概率密度π(ω)。在假设的参数ω条件下，可得到样本的条件似然函数［2］：

式中，p(xi|ωi)为在参数ωi的条件下，产生样本xi的概率密度。

样本与参数的联合概率密度：

式中，m(x)是样本x的概率密度，由式（3）边缘积分确定，且不包含参数w的信息。

参数w的后验概率密度如下：

式（4）即贝叶斯参数估计的概率密度函数形式，排除了和参数w无关的信息，使得后验概率达到最大的参数即所求参数，属于贝叶斯点估计问题。

2 次声波定位方法

2.1 建立后验概率模型

设各次声台阵水平方位角向量θ[θ1,…,θn]和到达时间向量t[t1,…,tn]，其中n为台阵个数。每个台阵的位置可由(xi,yi)表示，此为笛卡尔坐标系，可由大地坐标系转化得到。

样本为各台阵的探测结果信息d=[t,θ]，待估计参数m={t0,x0,y0,v}为次声源发生时刻、位置、传播速度，则代估参数的后验概率密度分布函数为：

式中：P(m)为模型参数的先验概率密度分布函数；P(m|d)为观测数据d在给定参数m下的条件似然函数；c(d)为归一化因子，用于保证P(m|d)的积分等于1。

2.2 先验概率密度分布函数

由于参数向量m={t0,x0,y0,v}各变量间的独立性，其先验概率可表示为：

由于传播度受到传播方向、距离以及大气高度等范围内的大气参数影响，很难得到准确的声速分布。为约束从声源到次声台阵未知的传播速度，本文将传播速度约束在一定的范围内，定义v的先验概率密度分布为［3］：

式中，0.28～0.34 km/s 表示传播速度可能的分布范围，常数16.67 用于保证p(v)积分等于1。

根据贝叶斯估计理论［4］，在无任何先验知识的前提下，可取p(t0)p(x0,y0)=1。

2.3 条件似然函数

由于各次声台阵探测结果相互独立，各台阵计算得到的目标方位角和到达时间也相互独立，样本的条件似然函数可定义为各台阵水平方位角和到达时间分量的乘积：

对于台阵i，iΘ度量了观测的水平方位角与所选模型参数的符合度，Φi度量了观测到达时间与给定参数下预测到达时间的误差。

两个分量的测量误差均符合高斯分布，其中水平方位角的似然函数分量为：

到达时间似然函数分量为：

令di=di(x0,y0,xi,yi)表示候选声源位置到第i个台阵的距离，则水平方位角和到达时间偏差为：

需要说明的是，水平方位角和到达时间的方差均包含测量误差和模型误差。

将式（11）和式（12）带入式（9）和式（10），并将x0、y0、t0用x、y、t表示，可得：

令：

则有：

由于式（15）是单调递增函数，取最大值时等价于函数F(x,y,t,v)取得最大值。

2.4 测量误差分析

假设水平方位角θi和到达时间ti观测值与预测值的误差来源于两个不相关的误差项之和，即测量误差和模型误差：

2.5 次声源定位

根据以上分析，在进行次声源定位时一般采取偏导数求解方法和网格搜索方法两种方法。

（1）偏导数求解方法。对函数F求偏导数，并另偏导数等于零，求解方程组可得到目标解：

利用式（18）进行求解时，需要解算复杂的非线性方程组。要利用优化算法进行迭代求解，迭代求解的初始值可利用各台阵探测到的目标方位进行交叉定位获得。

（2）网格搜索方法。一般做法是在代估计参数m={t0,x0,y0,v}的取值空间进行网格搜索，找到使得后验概率达到最大的参数点即为所求。

对比这两种方法，网格搜索方法计算简单、运算快速，实际中可划定重点关注区域，从而在该区域内进行搜索快速获得目标位置。

2.6 仿真分析

以两个台阵为例进行分析，采用网格搜索方法，设两台阵间距分别为800 km，台阵地址误差为2 m，通过500 次Monte-Carlo 仿真可以得到信号频率为0.1 Hz 和1 Hz 时，目标距离为400 km 的平均定位误差均小于10 km，目标距离为1 000 km 的平均定位误差均小于30 km，目标距离为2 000 km 的平均定位误差均小于80 km，如图1 所示。

3 结论

通过分析可见，在有众多不确定因素的条件下，利用基于概率的定位方法是解决次声远距离定位的有效方法，具有一定的实际应用价值。