王志峰，溫宏波，楊立學(xué)

（1.中國(guó)電子科技集團(tuán)公司第三研究所，北京 100015；2.陸軍裝備部駐北京地區(qū)軍事代表局，北京 100015）

本文采用基于概率的方法，將次聲傳播速度視為未知參數(shù)，同時(shí)將聲源位置坐標(biāo)、時(shí)間發(fā)生時(shí)間視為未知參數(shù)，通過(guò)計(jì)算各次聲分別測(cè)得目標(biāo)方位和到達(dá)時(shí)間的條件下這些未知參數(shù)的最大后驗(yàn)概率，可得到目標(biāo)聲源的準(zhǔn)確位置，本質(zhì)上屬于基于概率的定位方法［1］。

1 基本原理

在貝葉斯統(tǒng)計(jì)中，p(x|ω)表示在未知參數(shù)ω給定的條件下總體的概率密度，其中未知參數(shù)ω有先驗(yàn)概率密度π(ω)。在假設(shè)的參數(shù)ω條件下，可得到樣本的條件似然函數(shù)［2］：

式中，p(xi|ωi)為在參數(shù)ωi的條件下，產(chǎn)生樣本xi的概率密度。

樣本與參數(shù)的聯(lián)合概率密度：

式中，m(x)是樣本x的概率密度，由式（3）邊緣積分確定，且不包含參數(shù)w的信息。

參數(shù)w的后驗(yàn)概率密度如下：

式（4）即貝葉斯參數(shù)估計(jì)的概率密度函數(shù)形式，排除了和參數(shù)w無(wú)關(guān)的信息，使得后驗(yàn)概率達(dá)到最大的參數(shù)即所求參數(shù)，屬于貝葉斯點(diǎn)估計(jì)問題。

2 次聲波定位方法

2.1 建立后驗(yàn)概率模型

設(shè)各次聲臺(tái)陣水平方位角向量θ[θ1,…,θn]和到達(dá)時(shí)間向量t[t1,…,tn]，其中n為臺(tái)陣個(gè)數(shù)。每個(gè)臺(tái)陣的位置可由(xi,yi)表示，此為笛卡爾坐標(biāo)系，可由大地坐標(biāo)系轉(zhuǎn)化得到。

樣本為各臺(tái)陣的探測(cè)結(jié)果信息d=[t,θ]，待估計(jì)參數(shù)m={t0,x0,y0,v}為次聲源發(fā)生時(shí)刻、位置、傳播速度，則代估參數(shù)的后驗(yàn)概率密度分布函數(shù)為：

式中：P(m)為模型參數(shù)的先驗(yàn)概率密度分布函數(shù)；P(m|d)為觀測(cè)數(shù)據(jù)d在給定參數(shù)m下的條件似然函數(shù)；c(d)為歸一化因子，用于保證P(m|d)的積分等于1。

2.2 先驗(yàn)概率密度分布函數(shù)

由于參數(shù)向量m={t0,x0,y0,v}各變量間的獨(dú)立性，其先驗(yàn)概率可表示為：

由于傳播度受到傳播方向、距離以及大氣高度等范圍內(nèi)的大氣參數(shù)影響，很難得到準(zhǔn)確的聲速分布。為約束從聲源到次聲臺(tái)陣未知的傳播速度，本文將傳播速度約束在一定的范圍內(nèi)，定義v的先驗(yàn)概率密度分布為［3］：

式中，0.28～0.34 km/s 表示傳播速度可能的分布范圍，常數(shù)16.67 用于保證p(v)積分等于1。

根據(jù)貝葉斯估計(jì)理論［4］，在無(wú)任何先驗(yàn)知識(shí)的前提下，可取p(t0)p(x0,y0)=1。

2.3 條件似然函數(shù)

由于各次聲臺(tái)陣探測(cè)結(jié)果相互獨(dú)立，各臺(tái)陣計(jì)算得到的目標(biāo)方位角和到達(dá)時(shí)間也相互獨(dú)立，樣本的條件似然函數(shù)可定義為各臺(tái)陣水平方位角和到達(dá)時(shí)間分量的乘積：

對(duì)于臺(tái)陣i，iΘ度量了觀測(cè)的水平方位角與所選模型參數(shù)的符合度，Φi度量了觀測(cè)到達(dá)時(shí)間與給定參數(shù)下預(yù)測(cè)到達(dá)時(shí)間的誤差。

兩個(gè)分量的測(cè)量誤差均符合高斯分布，其中水平方位角的似然函數(shù)分量為：

到達(dá)時(shí)間似然函數(shù)分量為：

令di=di(x0,y0,xi,yi)表示候選聲源位置到第i個(gè)臺(tái)陣的距離，則水平方位角和到達(dá)時(shí)間偏差為：

需要說(shuō)明的是，水平方位角和到達(dá)時(shí)間的方差均包含測(cè)量誤差和模型誤差。

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

令：

則有：

由于式（15）是單調(diào)遞增函數(shù)，取最大值時(shí)等價(jià)于函數(shù)F(x,y,t,v)取得最大值。

2.4 測(cè)量誤差分析

假設(shè)水平方位角θi和到達(dá)時(shí)間ti觀測(cè)值與預(yù)測(cè)值的誤差來(lái)源于兩個(gè)不相關(guān)的誤差項(xiàng)之和，即測(cè)量誤差和模型誤差：

2.5 次聲源定位

根據(jù)以上分析，在進(jìn)行次聲源定位時(shí)一般采取偏導(dǎo)數(shù)求解方法和網(wǎng)格搜索方法兩種方法。

（1）偏導(dǎo)數(shù)求解方法。對(duì)函數(shù)F求偏導(dǎo)數(shù)，并另偏導(dǎo)數(shù)等于零，求解方程組可得到目標(biāo)解：

利用式（18）進(jìn)行求解時(shí)，需要解算復(fù)雜的非線性方程組。要利用優(yōu)化算法進(jìn)行迭代求解，迭代求解的初始值可利用各臺(tái)陣探測(cè)到的目標(biāo)方位進(jìn)行交叉定位獲得。

（2）網(wǎng)格搜索方法。一般做法是在代估計(jì)參數(shù)m={t0,x0,y0,v}的取值空間進(jìn)行網(wǎng)格搜索，找到使得后驗(yàn)概率達(dá)到最大的參數(shù)點(diǎn)即為所求。

對(duì)比這兩種方法，網(wǎng)格搜索方法計(jì)算簡(jiǎn)單、運(yùn)算快速，實(shí)際中可劃定重點(diǎn)關(guān)注區(qū)域，從而在該區(qū)域內(nèi)進(jìn)行搜索快速獲得目標(biāo)位置。

2.6 仿真分析

以兩個(gè)臺(tái)陣為例進(jìn)行分析，采用網(wǎng)格搜索方法，設(shè)兩臺(tái)陣間距分別為800 km，臺(tái)陣地址誤差為2 m，通過(guò)500 次Monte-Carlo 仿真可以得到信號(hào)頻率為0.1 Hz 和1 Hz 時(shí)，目標(biāo)距離為400 km 的平均定位誤差均小于10 km，目標(biāo)距離為1 000 km 的平均定位誤差均小于30 km，目標(biāo)距離為2 000 km 的平均定位誤差均小于80 km，如圖1 所示。

3 結(jié)論

通過(guò)分析可見，在有眾多不確定因素的條件下，利用基于概率的定位方法是解決次聲遠(yuǎn)距離定位的有效方法，具有一定的實(shí)際應(yīng)用價(jià)值。