基于MCMC的無(wú)線傳感器網(wǎng)絡(luò)目標(biāo)定位算法

張沛朋，李冬鋒

(1.濟(jì)源職業(yè)技術(shù)學(xué)院 藝術(shù)設(shè)計(jì)系，河南 濟(jì)源 459000； 2.華北水利水電大學(xué) 信息工程學(xué)院，河南 鄭州 450045)

0 引 言

由于在智能交通、目標(biāo)跟蹤等軍事和民用領(lǐng)域的廣泛應(yīng)用，基于時(shí)延-多普勒頻率的無(wú)線傳感器網(wǎng)絡(luò)(wireless sensor networks，WSN)目標(biāo)定位問(wèn)題近年來(lái)備受關(guān)注[1-4]。針對(duì)該問(wèn)題，現(xiàn)有算法可分成兩類：兩步定位和直接定位(direct position determination，DPD)。兩步定位的基本框架是：第一步從接收信號(hào)中提取時(shí)延[5,6]和多普勒[7,8]參數(shù)；第二步基于時(shí)延-多普勒，估計(jì)目標(biāo)位置參數(shù)[7,8]。兩步定位會(huì)損失一部分有用信息，導(dǎo)致最終目標(biāo)定位誤差增大，特別是在信噪比或信號(hào)采樣點(diǎn)數(shù)較低時(shí)。近年來(lái)，隨著通信帶寬和芯片計(jì)算能力的提升，直接定位開始得到關(guān)注，并發(fā)展出一系列有效算法[9-15]，這類直接定位算法的思路是從接收信號(hào)中直接估計(jì)出目標(biāo)位置參數(shù)，無(wú)需提取中間定位參數(shù)，避免了額外的信息損失，因此比傳統(tǒng)兩步定位具有更高的精度和穩(wěn)健性。然而，上述直接定位算法[9-15]存在計(jì)算復(fù)雜度過(guò)高、初值依賴和局部收斂問(wèn)題，特別是對(duì)于目標(biāo)位置-速度估計(jì)這樣的高維問(wèn)題求解。

為提高基于時(shí)延-多普勒的WSN目標(biāo)定位的精度、穩(wěn)健性和計(jì)算效率，本文提出一種基于馬爾科夫鏈-蒙特卡羅(Markov-chain Monte-Carlo，MCMC)[16-18]的直接定位算法。首先從傳感器接收信號(hào)中推導(dǎo)目標(biāo)位置和速度估計(jì)的優(yōu)化函數(shù)，然后將其轉(zhuǎn)化為目標(biāo)位置和速度的分布函數(shù)，利用MCMC方法對(duì)分布函數(shù)抽樣，通過(guò)統(tǒng)計(jì)樣本均值得到目標(biāo)位置和速度估計(jì)值。最后通過(guò)計(jì)算機(jī)仿真實(shí)驗(yàn)將估計(jì)結(jié)果與現(xiàn)有算法對(duì)比，突出本文算法的性能優(yōu)勢(shì)。

1 定位場(chǎng)景與信號(hào)模型

本文考慮的WSN運(yùn)動(dòng)目標(biāo)定位場(chǎng)景如圖1所示。

圖1 WSN目標(biāo)定位場(chǎng)景

根據(jù)上述場(chǎng)景假設(shè)，目標(biāo)與傳感器節(jié)點(diǎn)m在第k個(gè)時(shí)隙內(nèi)的距離及其變化率可以表示為

(1)

(2)

rm,k(t)=am,kxk(t-τm,k)eι2πfm,kt+wm,k(t)

(3)

假設(shè)以頻率fs對(duì)各傳感器接收信號(hào)進(jìn)行采樣，得到式(3)的離散形式如下

rm,k(n)=am,kxk(n-τm,kfs)eι2πfm,kn/fs+wm,k(n)

(4)

式中：n=0,1,…,N-1，N為每個(gè)時(shí)隙內(nèi)的采樣點(diǎn)數(shù)。

定義如下向量

n=[0,1,…,N-1]T

(5)

rm,k=[rm,k(0),rm,k(1),…,rm,k(N-1)]T

(6)

xk=[xk(0),xk(1),…,xk(N-1)]T

(7)

Dm,k=diag{eι2πfm,kn/fs}

(8)

Γm,k=diag{e-ι2πfsnτm,k/N}

(9)

F=[e-ι2π(0)n/N,e-ι2π(1)n/N,…,e-ι2π(N-1)n/N]T

(10)

wm,k=[wm,k(0),wm,k(1),…,wm,k(N-1)]T

(11)

則傳感器節(jié)點(diǎn)m在第k個(gè)時(shí)隙內(nèi)接收到的目標(biāo)基帶信號(hào)可以表示為如下的向量形式

rm,k=am,kHm,kxk+wm,k

(12)

其中，Hm,k=Dm,kFHΓm,kF，m=1,2,…,M，k=1,2,…,K。

2 基于MCMC的直接定位算法

為了以較低的計(jì)算復(fù)雜度實(shí)現(xiàn)對(duì)目標(biāo)位置和速度的估計(jì)，本節(jié)提出一種基于MCMC的直接定位算法。算法具體過(guò)程如下。

2.1 直接定位優(yōu)化函數(shù)構(gòu)建

(13)

(14)

式(14)中的未知參量中，含有冗余參數(shù)am,k， 其最大似然估計(jì)為

(15)

(16)

(17)

(18)

那么，目標(biāo)位置和速度的估計(jì)可以通過(guò)式(18)中目標(biāo)函數(shù)的極大化得到。

2.2 目標(biāo)位置參數(shù)求解

2.2.1 優(yōu)化函數(shù)全局最大化

根據(jù)Pincus提出的優(yōu)化理論[19]，式(18)中θ的估計(jì)問(wèn)題可以轉(zhuǎn)化為如下積分形式

(19)

(20)

(21)

根據(jù)伯努利大數(shù)定律，當(dāng)隨機(jī)樣本數(shù)量Ns足夠大時(shí)，式(21)的數(shù)值運(yùn)算結(jié)果即可收斂至式(19)積分結(jié)果。至此，積分問(wèn)題被轉(zhuǎn)化為了從概率密度函數(shù)pρ(θ) 中產(chǎn)生隨機(jī)樣本，也即pρ(θ) 的抽樣問(wèn)題。然而，由于pρ(θ) 的形式比較復(fù)雜，因此要從pρ(θ) 抽取樣本并不容易，而MCMC方法就是為此目的而誕生的。

2.2.2 MCMC方法

MCMC方法由兩個(gè)MC組成：第一個(gè)“MC”為“馬爾科夫鏈(Markov Chain)”縮寫，意為利用馬爾科夫鏈從目標(biāo)分布中抽樣；第二個(gè)“MC”為“蒙特卡羅(Monte Carlo)”縮寫，意為根據(jù)抽取的樣本，利用蒙特卡羅方法對(duì)積分進(jìn)行計(jì)算。具體來(lái)說(shuō)，在MCMC方法中，首先建立一個(gè)馬爾科夫鏈，使得目標(biāo)分布pρ(θ) 為其平穩(wěn)分布p(θ)， 即p(θ)=pρ(θ)， 然后運(yùn)行此馬爾科夫鏈充分長(zhǎng)時(shí)間直至其收斂至平穩(wěn)分布p(θ)， 那么從達(dá)到平穩(wěn)分布的馬爾科夫鏈中產(chǎn)生樣本，即可實(shí)現(xiàn)對(duì)目標(biāo)分布pρ(θ) 的抽樣。Metropolis-Hastings抽樣是一類常用的構(gòu)造馬爾科夫鏈的方法，包括了Metropolis抽樣、Gibbs抽樣等。對(duì)于一維隨機(jī)變量θ， Metropolis抽樣的基本流程如下：

(1)構(gòu)造提議分布q(·|θ(l))；

(2)在待估參數(shù)空間內(nèi)初始化馬爾科夫鏈的狀態(tài)θ(0)；

(3)對(duì)于l=0,1,2,…,重復(fù)以下過(guò)程進(jìn)行抽樣：

1)從提議分布q(·|θ(l)) 中產(chǎn)生候選狀態(tài)θ*；

2)生成均勻分布隨機(jī)數(shù)u～(0,1)；

3)若u<α(θ(l),θ*)， 其中α(θ(l),θ*) 為接受概率

(22)

則將馬爾科夫鏈狀態(tài)轉(zhuǎn)移至候選狀態(tài)，即令θ(l+1)=θ*； 否則馬爾科夫鏈維持當(dāng)前狀態(tài)，即令θ(l+1)=θ(l)。

4)增加l， 返回到步驟1)。

上述過(guò)程中，提議分布q(·|θ(l)) 應(yīng)為一個(gè)程序可抽樣的分布，通常可用隨機(jī)游走法構(gòu)造。隨機(jī)游走法假設(shè)候選狀態(tài)θ*從一個(gè)對(duì)稱的提議分布q(θ*|θ(l)) 中產(chǎn)生，即在每一次迭代中，從提議分布中產(chǎn)生增量Δθ，然后令θ*=θ(l)+Δθ。 典型地，可用高斯分布產(chǎn)生增量Δθ， 此時(shí)候選狀態(tài)服從高斯分布θ*～增量方差對(duì)生成馬爾科夫鏈的效果有重要影響。一般來(lái)說(shuō)，增量方差過(guò)大，會(huì)導(dǎo)致大部分候選樣本被拒絕，此時(shí)算法效率很低。如果增量方差過(guò)小，則候選樣本幾乎都被接受或拒絕，鏈幾乎不再游走，同樣效率較低。一般來(lái)說(shuō)，接受率在0.5～0.85才能保證得到的馬爾科夫鏈具有較好的性質(zhì)[20]。為此，可以設(shè)置增量方差的調(diào)節(jié)參數(shù)，并在游走過(guò)程中監(jiān)視馬爾科夫鏈的接受率，以使算法達(dá)到預(yù)設(shè)的接受率。

(23)

(24)

生成均勻分布隨機(jī)數(shù)u～(0,1)： 若令否則令

2.2.3 基于MCMC的目標(biāo)位置參數(shù)求解

基于上述MCMC方法，對(duì)于本文WSN目標(biāo)定位問(wèn)題，待估參數(shù)為目標(biāo)位置和速度參數(shù)構(gòu)成的2dim×1維列向量，因此本文采用Gibbs抽樣生成馬爾科夫鏈，采用隨機(jī)游走法構(gòu)造提議函數(shù)。算法具體實(shí)現(xiàn)過(guò)程如下：

(1)設(shè)置狀態(tài)轉(zhuǎn)移次數(shù)閾值Nb(保證馬爾科夫鏈?zhǔn)諗?，需要的樣本數(shù)量Ns；

(3)令l=0，j=1， 重復(fù)如下抽樣過(guò)程：

1)利用隨機(jī)游走法構(gòu)造提議分布

(25)

4)生成隨機(jī)數(shù)u～(0,1)：

5)調(diào)整隨機(jī)游走步長(zhǎng)：

6)如果j

7)如果l

(4)根據(jù)抽取的樣本θ(0),θ(1),…,θ(Nb+Ns)， 估計(jì)目標(biāo)位置和速度參數(shù)的估計(jì)值

(26)

3 克拉美羅界分析

(27)

其中， Re(·) 和Im(·) 分別表示復(fù)數(shù)的實(shí)部和虛部。那么，由式(13)中向量r的概率密度函數(shù)，可得待估參數(shù)的Fisher信息矩陣為

(28)

式中：H表示矩陣共軛轉(zhuǎn)置運(yùn)算。根據(jù)待估參量φ的構(gòu)成，偏導(dǎo)矩陣?r/?φT可以表示為

(29)

根據(jù)復(fù)變函數(shù)求導(dǎo)規(guī)則，式(29)右側(cè)子矩陣可以表示為

(30)

(31)

J2=diag{J2,1,J2,2,…,J2,K}

(32)

J2,k=diag{H1,kxk,H2,kxk,…,HM,kxk}

(33)

(34)

(35)

Λ1=ι2πdiag{n}

(36)

(37)

(38)

(39)

(40)

(41)

(42)

(43)

根據(jù)克拉美羅界的定義，待估參量φ的CRLB為Fisher信息矩陣的逆，即

CRLB(φ)=[FIM(φ)]-1

(44)

4 仿真實(shí)驗(yàn)

(45)

載頻fc=10 MHz。 目標(biāo)信號(hào)能夠同時(shí)被WSN的M=4個(gè)傳感器節(jié)點(diǎn)所接收。每個(gè)傳感器節(jié)點(diǎn)在K=50個(gè)時(shí)隙段對(duì)目標(biāo)輻射的電磁信號(hào)進(jìn)行接收和采樣，采樣頻率為fs=1 MHz， 每個(gè)時(shí)隙段持續(xù)時(shí)間為T=1 ms， 接收信號(hào)的幅度am,k=1。 WSN各傳感器節(jié)點(diǎn)的位置和速度見(jiàn)表1。為了模擬實(shí)際的信號(hào)數(shù)據(jù)，我們?yōu)楦鱾鞲衅鞴?jié)點(diǎn)接收信號(hào)添加零均值的高斯白噪聲，噪聲方差根據(jù)仿真實(shí)驗(yàn)場(chǎng)景調(diào)整。

算法的估計(jì)性能通過(guò)500次獨(dú)立計(jì)算機(jī)仿真實(shí)驗(yàn)統(tǒng)計(jì)得到的目標(biāo)位置和速度均方根誤差(root mean square error，RMSE)來(lái)衡量，定義如下

(46)

(47)

表1 傳感器節(jié)點(diǎn)位置和速度

首先，為了直觀展示本文所提算法的有效性，我們?cè)O(shè)置噪聲方差σ2=1， 然后利用本文所提算法對(duì)目標(biāo)進(jìn)行定位，算法的目標(biāo)函數(shù)如圖2所示。

圖2 目標(biāo)函數(shù)

從圖2(a)可以看出，目標(biāo)信號(hào)在其真實(shí)位置參數(shù)處得到峰值，并且得益于多個(gè)時(shí)隙的長(zhǎng)時(shí)間信號(hào)能量積累，原本淹沒(méi)在噪聲中的目標(biāo)信號(hào)數(shù)據(jù)積累后的峰值遠(yuǎn)高于周圍噪聲。從圖2(b)中可以看到，通過(guò)調(diào)整參數(shù)ρ的值，目標(biāo)函數(shù)的峰值被進(jìn)一步強(qiáng)化，噪聲被一定程度抑制，這將有利于后續(xù)MCMC方法抽樣。

本文算法利用MCMC方法抽取樣本的效果如圖3所示。

圖3 MCMC方法抽樣(ρ=5)

從圖3中生成的馬爾科夫鏈可以看到，抽樣過(guò)程開始后，隨機(jī)游走的增量方差較大，體現(xiàn)在馬爾科夫鏈的路徑步長(zhǎng)較大。較大的步長(zhǎng)有利于馬爾科夫鏈迅速遍歷參數(shù)空間，并快速收斂至平穩(wěn)分布。收斂至平穩(wěn)分布后，馬爾科夫鏈的樣本在目標(biāo)位置參數(shù)附近游走，步長(zhǎng)較小。對(duì)生成樣本統(tǒng)計(jì)發(fā)現(xiàn)，MCMC抽樣的接受率在0.7左右，說(shuō)明了生成馬爾科夫鏈具有較好的性質(zhì)。

如上文所示，本文所提算法涉及到的兩個(gè)重要參數(shù)——樣本數(shù)量和參數(shù)ρ。接下來(lái)通過(guò)仿真實(shí)驗(yàn)分析樣本數(shù)量和參數(shù)ρ對(duì)定位精度的影響。仿真實(shí)驗(yàn)結(jié)果如圖4所示。

圖4 MCMC方法參數(shù)對(duì)定位精度的影響

首先，為定量分析樣本數(shù)量對(duì)目標(biāo)定位精度的影響，并確定樣本數(shù)量的最優(yōu)選取，我們統(tǒng)計(jì)了不同樣本數(shù)量下的目標(biāo)位置和速度估計(jì)RMSE，結(jié)果如圖4(a)所示。從圖4(a)可以看到，與直觀預(yù)期一致，整體上隨著樣本數(shù)量的增加，目標(biāo)位置和速度估計(jì)的RMSE呈下降趨勢(shì)。然而，當(dāng)樣本數(shù)量增加至3000以后，目標(biāo)位置和速度估計(jì)的RMSE下降速度趨近于0，此時(shí)再增加樣本數(shù)量將無(wú)法顯著提高目標(biāo)位置和速度的估計(jì)精度，反而會(huì)增加計(jì)算復(fù)雜度。因此，對(duì)于本文所設(shè)置的定位場(chǎng)景，樣本數(shù)量設(shè)置為3000為宜，后續(xù)仿真實(shí)驗(yàn)亦如此設(shè)置。

然后，為了確定參數(shù)ρ的合理取值，我們?cè)诓煌膮?shù)ρ下進(jìn)行定位仿真，并統(tǒng)計(jì)目標(biāo)位置和速度估計(jì)RMSE，結(jié)果如圖4(b)所示。從圖4(b)可以看到，在參數(shù)ρ的取值較小時(shí)，與前文所述一致，增大參數(shù)ρ的值可以顯著提高目標(biāo)位置和速度的估計(jì)精度，這是因?yàn)橐粋€(gè)較大的ρ可以抑制噪聲影響，使得目標(biāo)分布峰值更加尖銳，繼而使得生成樣本更加集中于目標(biāo)位置處；然而，當(dāng)參數(shù)ρ取值過(guò)大時(shí)，也會(huì)導(dǎo)致目標(biāo)位置和速度的估計(jì)精度降低，這是因?yàn)閰?shù)ρ取值過(guò)大會(huì)導(dǎo)致目標(biāo)分布峰值過(guò)于尖銳，從而導(dǎo)致MCMC抽樣困難，并且參數(shù)ρ取值越大意味著指數(shù)運(yùn)算次數(shù)越高，計(jì)算復(fù)雜度越高。因此，對(duì)于本文定位場(chǎng)景，參數(shù)ρ取值為10左右為宜。后續(xù)仿真實(shí)驗(yàn)中設(shè)置為ρ=10。

接下來(lái)，為了評(píng)估算法在不同信噪比條件下的定位性能，我們將各傳感器節(jié)點(diǎn)接收信號(hào)添加不同方差的噪聲，從而得到不同信噪比的接收信號(hào)，并在不同的信噪比條件下利用本文算法進(jìn)行仿真定位實(shí)驗(yàn)。為了突出本文算法的定位性能，我們將現(xiàn)有5種算法，包括文獻(xiàn)[7]中傳統(tǒng)的兩步定位方法、文獻(xiàn)[9]中基于網(wǎng)格搜索的直接定位算法(DPD-GS)、文獻(xiàn)[13]中基于期望最大化的直接定位算法(DPD-EM)、文獻(xiàn)[14]中基于粒子濾波的直接定位算法(DPD-PF)、文獻(xiàn)[15]中基于改進(jìn)粒子濾波的直接定位算法(DPD-MPF)，作為對(duì)比算法。對(duì)比結(jié)果如圖5所示。

圖5 不同信噪比條件下算法的定位精度比較

圖5中結(jié)果表明，隨著接收信號(hào)信噪比的提升，幾種算法的定位精度總體上也都隨之提升，但幾種算法的定位精度卻不盡相同。傳統(tǒng)的兩步定位算法由于其兩步處理體制會(huì)損失一部分信息，導(dǎo)致最終目標(biāo)位置和速度估計(jì)RMSE始終明顯高于其余幾種直接定位算法以及CRLB。DPD-GS算法通過(guò)網(wǎng)格搜索確定目標(biāo)位置和速度參數(shù)，定位精度受到搜索網(wǎng)格大小的限制，在高信噪比條件下定位精度無(wú)法進(jìn)一步提高。DPD-EM算法通過(guò)期望最大化迭代算法確定目標(biāo)位置和速度參數(shù)，在高信噪比條件下定位精度可以接近CRLB，但是當(dāng)信噪比較低時(shí)，目標(biāo)位置和速度估計(jì)RMSE高于DPD-PF、DPD-MPF以及本文DPD-MCMC算法，這說(shuō)明了期望最大化算法作為一種局部最優(yōu)算法，在低信噪比條件下的穩(wěn)健性較差。DPD-PF算法和DPD-MPF算法的定位性能相當(dāng)，總體上定位誤差略高于本文DPD-MCMC算法。本文DPD-MCMC算法的定位誤差要低于現(xiàn)有幾種算法，并且在一般的信噪比條件下定位精度可以達(dá)到CRLB。

算法的平均運(yùn)行時(shí)間也是衡量算法性能的重要指標(biāo)。因此，我們?cè)谙嗤挠?jì)算機(jī)仿真環(huán)境下比較了現(xiàn)有算法的平均運(yùn)行時(shí)間。計(jì)算機(jī)仿真環(huán)境配置如下：電腦Lenovo ThinkCentre M920t；內(nèi)存32 GB；處理器Intel I9-9900(8核，3.1 GHz)；操作系統(tǒng)Win10 64位專業(yè)版；軟件：MATLAB R2019b。仿真結(jié)果見(jiàn)表2。

表2 算法計(jì)算復(fù)雜度比較

表2中6種算法的平均運(yùn)行時(shí)間比較結(jié)果表明，直接定位算法的計(jì)算復(fù)雜度要高于傳統(tǒng)的兩步定位算法，但是直接定位算法可以得到更高的定位精度。在幾種直接定位算法的比較中，本文DPD-MCMC算法的計(jì)算復(fù)雜度要顯著低于DPD-GS、DPD-PF、DPD-MPF算法，與DPD-EM算法相當(dāng)。這是由于本文DPD-MCMC算法通過(guò)MCMC方法獲取參數(shù)樣本，繼而利用樣本統(tǒng)計(jì)結(jié)果實(shí)現(xiàn)目標(biāo)位置和速度參數(shù)估計(jì)的方式，無(wú)需像DPD-GS算法那樣在多維空間內(nèi)進(jìn)行網(wǎng)格搜索，也無(wú)需像DPD-PF和DPD-MPF算法那樣產(chǎn)生大量粒子，有效降低了計(jì)算復(fù)雜度。需要指出的是，雖然本文DPD-MCMC算法在本文仿真環(huán)境下平均需要2.68 s時(shí)間才能估計(jì)出目標(biāo)位置和速度，但是利用更為專業(yè)和強(qiáng)大的計(jì)算平臺(tái)，運(yùn)行時(shí)間有望進(jìn)一步縮短至可對(duì)目標(biāo)實(shí)時(shí)定位的程度。

5 結(jié)束語(yǔ)

本文提出了一種針對(duì)無(wú)線傳感器網(wǎng)絡(luò)的目標(biāo)直接定位方法。首先依據(jù)最大似然準(zhǔn)則推導(dǎo)了基于時(shí)延和多普勒頻率的運(yùn)動(dòng)目標(biāo)直接定位數(shù)學(xué)優(yōu)化模型。針對(duì)該模型以矩陣特征值的形式給出而難以解析求解的問(wèn)題，提出了一種基于馬爾科夫鏈-蒙特卡羅的數(shù)值求解方法。此外，本文推導(dǎo)了基于時(shí)延和多普勒頻率的無(wú)線傳感器網(wǎng)絡(luò)目標(biāo)定位克拉美羅界，從理論上分析了算法定位誤差的理論極限。計(jì)算機(jī)仿真實(shí)驗(yàn)結(jié)果表明：相比于傳統(tǒng)的兩步定位方法，本文算法具有更高的定位精度；相比于現(xiàn)有的直接定位算法，本文算法的計(jì)算復(fù)雜度更低，并且在一般信噪比條件下，定位性能能夠逼近克拉美羅界。