周榮艷, 陳建峰, 李曉強(qiáng), 譚偉杰

(1. 西北工業(yè)大學(xué)航海學(xué)院, 陜西 西安 710072; 2. 南陽(yáng)理工學(xué)院, 河南 南陽(yáng) 473004;3. 貴州大學(xué)公共大數(shù)據(jù)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 貴州 貴陽(yáng) 550025)

0 引 言

分布式定位系統(tǒng)的定位精度不僅依賴于傳感器節(jié)點(diǎn)的定位算法、節(jié)點(diǎn)個(gè)數(shù)、采樣頻率等因素,節(jié)點(diǎn)的部署位置也是影響系統(tǒng)定位性能的重要因素之一[1-4]。節(jié)點(diǎn)部署算法能夠在不改變定位算法和不增加節(jié)點(diǎn)個(gè)數(shù)的情況下,有效提高系統(tǒng)的定位精度,得到國(guó)內(nèi)外學(xué)者的廣泛研究,并逐步成為分布式定位系統(tǒng)研究中的一個(gè)重要分支[5-6]。

在目標(biāo)探測(cè)與定位領(lǐng)域,克拉美羅界(Cramer-Rao bound, CRB)提供了任何無(wú)偏估計(jì)量方差所能達(dá)到的一個(gè)理論下界,并且CRB的表達(dá)式是關(guān)于目標(biāo)與節(jié)點(diǎn)相對(duì)幾何位置的非線性函數(shù),因此常把CRB作為優(yōu)化節(jié)點(diǎn)與目標(biāo)相對(duì)位置的評(píng)價(jià)函數(shù)[7-10]。近年來(lái),關(guān)于CRB的研究主要集中在兩個(gè)方面:一是根據(jù)相關(guān)理論的推導(dǎo)得到針對(duì)某些特殊目標(biāo)位置或接收特殊陣列時(shí)距離與角度估計(jì)的CRB[11-13];二是根據(jù)某一設(shè)計(jì)指標(biāo)要求,按照相應(yīng)的優(yōu)化原則,將求解CRB的問(wèn)題轉(zhuǎn)化為最優(yōu)化問(wèn)題,采用智能優(yōu)化算法進(jìn)行迭代,尋優(yōu)從而逼近CRB[14-15]。這些方法都是分析如何減小定位估計(jì)誤差的CRB,雖然可以較好地提高目標(biāo)在某一特定位置的定位精度,但并沒(méi)有考慮目標(biāo)服從高斯分布時(shí)定位誤差的CRB。

文獻(xiàn)[16]基于到達(dá)時(shí)間(time of arrival, TOA)定位算法,以最小化平均CRB為優(yōu)化準(zhǔn)則,得出基于節(jié)點(diǎn)觀測(cè)性能的最優(yōu)布局與聲源出現(xiàn)概率分布直接相關(guān)的結(jié)論。文獻(xiàn)[17]研究了目標(biāo)服從高斯先驗(yàn)分布的條件下?tīng)顟B(tài)估計(jì)的協(xié)方差矩陣和信息矩陣,為傳感器的選擇和最優(yōu)布局提供了衡量標(biāo)準(zhǔn),并得出純方位估計(jì)中視距的方向向量垂直于協(xié)方差矩陣最大特征向量的結(jié)論。此結(jié)論應(yīng)用于文獻(xiàn)[18]中,并對(duì)目標(biāo)高斯先驗(yàn)分布下異構(gòu)傳感器網(wǎng)絡(luò)中的節(jié)點(diǎn)進(jìn)行優(yōu)化配置,但需要將初始協(xié)方差矩陣進(jìn)行若干次迭代更新,計(jì)算量大且較為復(fù)雜。

近年來(lái)有些學(xué)者通過(guò)分析坐標(biāo)旋轉(zhuǎn)在目標(biāo)定位和探測(cè)方面的應(yīng)用,為相關(guān)研究提供了新的思路。文獻(xiàn)[19]證明了二維笛卡爾坐標(biāo)系統(tǒng)中幾何坐標(biāo)平移與旋轉(zhuǎn)對(duì)純方位估計(jì)算法的偽線性估計(jì)(pseudolinear estimation, PLE)與最大似然估計(jì)(maximum likelihood estimation, MLE)沒(méi)有影響,而總體最小二乘(total least squares, TLS)估計(jì)通過(guò)坐標(biāo)旋轉(zhuǎn)發(fā)生變化,并依此為基礎(chǔ)提出了一種減小TLS估計(jì)偏差的算法,但并未考慮在三維空間中的相關(guān)問(wèn)題。文獻(xiàn)[20]通過(guò)對(duì)方位角和俯仰角的測(cè)量方程進(jìn)行線性化,將未知節(jié)點(diǎn)估計(jì)轉(zhuǎn)化為PLE模型求解,并證明了此三維偽線性估計(jì)算法通過(guò)笛卡爾坐標(biāo)系旋轉(zhuǎn)可以顯著改變算法性能,從而有效地減小系統(tǒng)偏差。但此算法分別基于xy坐標(biāo)軸與z坐標(biāo)軸,并未對(duì)三維空間中的坐標(biāo)系整體進(jìn)行研究。

本文在其他學(xué)者研究的基礎(chǔ)上采用到達(dá)角(angle of arrival, AOA)定位算法,首先基于目標(biāo)高斯分布的FIM信息矩陣,研究滿足CRB的條件;之后推導(dǎo)基于三維坐標(biāo)旋轉(zhuǎn)的最大后驗(yàn)概率估計(jì)(maximum a posteriori estimation, MAP)定位算法,將目標(biāo)先驗(yàn)概率分布的非對(duì)角協(xié)方差矩陣轉(zhuǎn)化為對(duì)角陣;得到滿足CRB的最小跡時(shí)節(jié)點(diǎn)的最優(yōu)部署位置。

1 問(wèn)題描述

在三維空間中,假設(shè)一個(gè)目標(biāo)聲源的位置s=[x,y,z]T,N個(gè)具有三維測(cè)向能力的節(jié)點(diǎn)位置坐標(biāo)pk=[px k,py k,pz k]T,k=1,2,…,N。如圖1所示,任一節(jié)點(diǎn)對(duì)目標(biāo)的觀測(cè)都有方位角φk和俯仰角θk兩個(gè)估計(jì)量,表示為

圖1 三維AOA定位算法中方位角和俯仰角示意圖Fig.1 Illustration of azimuth angle and elevation angle in three dimensional AOA positioning algorithm

(1)

(2)

在噪聲環(huán)境中,角度測(cè)量值為真實(shí)的角度測(cè)量值和測(cè)量誤差的疊加[21],即

(3)

(4)

Φ=M0+JTΣ-1J=

M0+Φ1+Φ2

(5)

(6)

CRB為任何無(wú)偏估計(jì)量方差的確定下界,CRB為費(fèi)希爾信息矩陣的逆矩陣。由文獻(xiàn)[22]可知,三維空間中AOA定位算法估計(jì)誤差的FIM可寫為

(7)

(8)

2 基于三維坐標(biāo)旋轉(zhuǎn)的MAP方法

2.1 MAP方法

p(s)=

(9)

由式(3)可得,角度測(cè)量值的似然函數(shù)為

(10)

(11)

(12)

最大后驗(yàn)成本函數(shù)JMAP(s) 可表示為

(13)

則殘差為

Γ(s)=[e(s);r(s)]

(14)

其中e(t)和r(t)為

r(s)=s0-s

(15)

(16)

J1i=

(17)

式中:

(18)

設(shè)3×3矩陣J2i表示r(s)的雅克比矩陣,可得

(19)

將式(17)與式(19)結(jié)合,即可得式(13)雅克比矩陣

Ji=-[J1i;J2i]

(20)

采用高斯-牛頓算法對(duì)MAP進(jìn)行計(jì)算,則可得迭代關(guān)系式[25]

(21)

2.2 基于三維坐標(biāo)旋轉(zhuǎn)的MAP方法

由式(8)可知,當(dāng)目標(biāo)位置誤差的協(xié)方差矩陣P0=diag([a,b,c])時(shí),則M0=diag([a-1,b-1,c-1]),由式(5)可知只需滿足等號(hào)右邊后兩項(xiàng)Φ1與Φ2最小時(shí)，得到tr(CRB)的最小值。當(dāng)P0為對(duì)稱矩陣時(shí),文獻(xiàn)[17]通過(guò)將初始協(xié)方差矩陣P0進(jìn)行若干次迭代更新，將其轉(zhuǎn)化為對(duì)角陣，尋找純方位定位算法中節(jié)點(diǎn)的最優(yōu)部署方向,但計(jì)算量大且較為復(fù)雜。

由于目標(biāo)的先驗(yàn)協(xié)方差矩陣P0在物理上相當(dāng)于一個(gè)橢球,其內(nèi)任一點(diǎn)都是目標(biāo)可能出現(xiàn)的位置。因?yàn)槿S坐標(biāo)旋轉(zhuǎn)并不影響橢球體的大小,且協(xié)方差矩陣P0對(duì)于任何相似變換UPUT具有不變性(其中U為酉矩陣)?？刹捎萌S坐標(biāo)旋轉(zhuǎn)的方法將P0對(duì)角化,通過(guò)合適的繞x軸、y軸、z軸旋轉(zhuǎn)角度α、β、γ組成的三維旋轉(zhuǎn)矩陣R,將初始協(xié)方差矩陣P0對(duì)角化,即可得到tr(CRB)的最小值。

在本節(jié)中首先定義三維空間中的旋轉(zhuǎn)矩陣,再對(duì)三維旋轉(zhuǎn)后的MAP方法進(jìn)行推導(dǎo)。定義三維空間中的坐標(biāo)旋轉(zhuǎn)矩陣如下所示。

(22)

其中α、β、γ分別為繞x軸、y軸、z軸逆時(shí)針的旋轉(zhuǎn)角,如圖2所示。

圖2 繞x軸、y軸、z軸逆時(shí)針的旋轉(zhuǎn)角 Fig 2 Rotation angle counterclockwise about the x -,y -,z - axis

R=RxRyRz

(23)

當(dāng)三維空間中進(jìn)行坐標(biāo)旋轉(zhuǎn)時(shí),可得旋轉(zhuǎn)后定位系統(tǒng)中各參量為

(24)

(25)

因此,旋轉(zhuǎn)后的俯仰角和方位角表示為

(26)

為了計(jì)算旋轉(zhuǎn)后的協(xié)方差矩陣可采用一階泰勒級(jí)數(shù)近似來(lái)計(jì)算旋轉(zhuǎn)后的角度噪聲。將式(3)代入到式(26)中,可得

(27)

由誤差傳播定律可知,在三維旋轉(zhuǎn)后,第k個(gè)節(jié)點(diǎn)的噪聲協(xié)方差矩陣可以寫為

(28)

旋轉(zhuǎn)后包含噪聲分量的角度估計(jì)為

(29)

旋轉(zhuǎn)后的tr和qr計(jì)算得到無(wú)噪聲的真實(shí)目標(biāo)位置可寫為

(30)

因此,旋轉(zhuǎn)后的MAP 可表示為

(31)

采用式(24)中旋轉(zhuǎn)后目標(biāo)先驗(yàn)分布中的協(xié)方差矩陣,結(jié)合上式中的Kr,則MAP的協(xié)方差矩陣為

(32)

則式(31)可寫為

(33)

殘差為

Γ(sr)=[e(sr);r(sr)]

(34)

式中:

(35)

(36)

采用高斯牛頓法計(jì)算MAP,迭代關(guān)系式如下所示

(37)

2.3 最優(yōu)節(jié)點(diǎn)部署方法

本節(jié)主要研究已知目標(biāo)先驗(yàn)分布的情況下純方位定位系統(tǒng)中節(jié)點(diǎn)的最優(yōu)部署。首先,由第2.2節(jié)的結(jié)論可知,在目標(biāo)聲源服從高斯分布的假設(shè)下,三維坐標(biāo)旋轉(zhuǎn)能夠解決目標(biāo)先驗(yàn)協(xié)方差矩陣為非對(duì)角陣的問(wèn)題,可通過(guò)三維坐標(biāo)旋轉(zhuǎn)的方法將目標(biāo)先驗(yàn)分布中的非對(duì)角協(xié)方差矩陣轉(zhuǎn)化為對(duì)角陣,所以在后續(xù)研究中,只考慮P0為對(duì)角陣時(shí)的CRB。由式(5)～式(8)可知CRB的跡為

tr(CRB)=tr(Φ-1)≥

(38)

當(dāng)以下等式同時(shí)成立時(shí),式(38)的等號(hào)成立。

(39)

由于通過(guò)三維旋轉(zhuǎn)后P0為對(duì)角陣,則使CRB對(duì)角化的條件為[21]

(40)

滿足式(40)的方位角的集合為

C={{θ1,θ2,…,θN}|θk∈{0,±π/2},k=1,2,…,N}

{{0,0,…,0},{±π/2,±π/2,…,±π/2}}

(41)

俯仰角的集合為

Z={{φ1,φ2,…,φN}|φk∈{0,±π/2},k=1,2,…,N}

{±π/2,±π/2,…,±π/2}

(42)

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

本節(jié)首先采用梯度下降算法,驗(yàn)證本文基于目標(biāo)高斯分布的定位系統(tǒng)中節(jié)點(diǎn)最優(yōu)部署方法。假設(shè)已知目標(biāo)的先驗(yàn)概率分布,系統(tǒng)中有3個(gè)節(jié)點(diǎn)進(jìn)行定位,并將節(jié)點(diǎn)位置隨意放置,文獻(xiàn)[21]中未經(jīng)三維旋轉(zhuǎn)的最優(yōu)部署位置和本文方法得到的最優(yōu)部署進(jìn)行對(duì)比,證明本文方法的優(yōu)越性。

3.1 最優(yōu)節(jié)點(diǎn)部署方法實(shí)驗(yàn)驗(yàn)證

采用牛頓梯度下降算法對(duì)定位誤差CRB的跡進(jìn)行最小化迭代尋優(yōu)[26],本實(shí)驗(yàn)仿真中對(duì)牛頓梯度下降算法執(zhí)行10 000次。

實(shí)驗(yàn) 1目標(biāo)聲源和傳感器之間的最小距離用r0表示,r01=r02=r03=100 m,s0=[0,0,0]T,協(xié)方差矩陣

σφ1=σφ2=0.5°,σφ3=1°。節(jié)點(diǎn)初始迭代坐標(biāo)位置為p1=[100,-200,0],p2=[300,200,200],p3=[-200,100,-200]。圖3(a)為節(jié)點(diǎn)從初始迭代位置移動(dòng)到最優(yōu)部署位置的迭代軌跡。目標(biāo)在最終位置時(shí)的角度為θ1=-91.07°,θ2=89.07°,θ3=-178.07°,φ1=-0.06°,φ2=-0.013°,φ3=-88.76°。

圖3 節(jié)點(diǎn)迭代軌跡圖Fig.3 Iterative trajectory of the sensors

實(shí)驗(yàn) 2在其他條件都不改變的條件下,設(shè)先驗(yàn)分布中的協(xié)方差矩陣為非對(duì)稱矩陣

通過(guò)三維坐標(biāo)旋轉(zhuǎn)的方法將其轉(zhuǎn)化為對(duì)稱陣,旋轉(zhuǎn)角度為α=14.5°,β=342.1°,γ=15.4°,可得到對(duì)角矩陣

圖3(b)為節(jié)點(diǎn)從最初迭代位置到最優(yōu)位置的迭代軌跡,最終角度為θ1=-92.12°,θ2=88.27°,θ3=-176.72°,φ1=-0.02°,φ2=-0.047°,φ3=-87.82°。

由此可知在協(xié)方差矩陣為對(duì)角陣和非對(duì)角陣的情況下,通過(guò)最小化CRB的跡得到節(jié)點(diǎn)最優(yōu)部署,與理論推導(dǎo)結(jié)論一致。

3.2 節(jié)點(diǎn)部署對(duì)比實(shí)驗(yàn)

本節(jié)將3種不同部署方法所得CRB的跡進(jìn)行對(duì)比,驗(yàn)證該算法的優(yōu)越性如表1所示。

表1 3種節(jié)點(diǎn)部署tr(CRB)的對(duì)比

部署方法1:假設(shè)節(jié)點(diǎn)部署位置即為第3.1節(jié)中節(jié)點(diǎn)的初始迭代位置,其余條件相同。

部署方法2:文獻(xiàn)[21]中3個(gè)節(jié)點(diǎn)的最優(yōu)部署位置,即θ1=0°,θ2=90°,θ3=90°,φ1=90°,φ2=0°,φ3=0°。

部署方法3:第3.1節(jié)中節(jié)點(diǎn)部署位置。

由表1可知,隨意位置部署時(shí)的tr(CRB)值較大,而文獻(xiàn)[21]的結(jié)果由于受高斯分布中協(xié)方差矩陣的影響,也不再適用于此情況,通過(guò)三維旋轉(zhuǎn)后得到的節(jié)點(diǎn)部署可獲得tr(CRB)的最小值。

4 結(jié) 論

為了解決定位系統(tǒng)中目標(biāo)位置服從高斯先驗(yàn)分布假設(shè)下的最優(yōu)節(jié)點(diǎn)部署問(wèn)題,本文探討了基于三維坐標(biāo)旋轉(zhuǎn)的MAP方法,提出了一種基于目標(biāo)高斯分布的求解克拉美羅下界的方法,并依此給出了最優(yōu)的節(jié)點(diǎn)部署配置,最后通過(guò)梯度下降方法驗(yàn)證了理論推導(dǎo)的正確性。由于本文中最優(yōu)部署方法是基于AOA定位算法,這也將在后續(xù)研究中為AOA節(jié)點(diǎn)定位算法的最優(yōu)布局和節(jié)點(diǎn)選擇等問(wèn)題提供了理論指導(dǎo)和評(píng)價(jià)標(biāo)準(zhǔn),以便在實(shí)際工程應(yīng)用中獲得最佳目標(biāo)定位與探測(cè)性能。