一種基于觀測站數(shù)目最小化的TDOA/FDOA無源定位算法

孫霆，董春曦，董陽陽，劉明明

西安電子科技大學(xué) 電子工程學(xué)院，西安 710071

無源定位算法是指定位系統(tǒng)不向被定位目標發(fā)射電磁波信號，而是僅利用目標輻射源的輻射信息，確定出目標位置的一種技術(shù)[1]，被廣泛應(yīng)用在無線傳感器網(wǎng)絡(luò)[2-4]、雷達[5-6]以及聲吶系統(tǒng)[7-8]等領(lǐng)域中。對于靜止目標，無源定位技術(shù)主要是依靠不同類型的觀測量，包括到達時間(Time of Arrival，TOA)、到達時間差(Time Difference of Arrival，TDOA)、到達角(Angle of Arrival，AOA)、以及上述參數(shù)的相互結(jié)合。而當觀測站與目標存在相互運動時，可以聯(lián)合TDOA以及到達頻率差(Fequency Difference of Arrival，F(xiàn)DOA)觀測量對目標的參數(shù)進行估計。

本文主要針對TDOA/FDOA定位算法進行研究，該定位體制中由于目標參數(shù)與觀測量存在非線性關(guān)系，因此求解變得比較困難。比較經(jīng)典的TDOA/FDOA定位算法包括一些迭代算法[9-12]、兩步加權(quán)最小二乘(Two-Stage Weighted Least Squares，TSWLS)算法[13-14]、半定松弛(SeminDefine Relaxation，SDR)算法[15-16]、多維標度(MultiDimensional Scaling, MDS)算法[17-18]。文獻[9-12]均為迭代算法并且都具有良好的估計性能，但是這些迭代算法的收斂性受初值影響較大，當初值選取較差時，定位精度無法保證，并且迭代結(jié)構(gòu)要求這些算法需要繁瑣的計算。

為了克服迭代算法的缺點，文獻[14]提出了經(jīng)典的兩步加權(quán)最小二乘算法。該算法首先通過引入中間變量建立一組關(guān)于TDOA與FDOA的偽線性方程，通過加權(quán)最小二乘(Weighted Least Squares，WLS)技術(shù)得到了目標參數(shù)與中間變量的初始估計值，隨后利用中間變量與目標之間的約束關(guān)系進一步提高精度。這種TSWLS算法是一種閉式解(Closed-Form Solution)，不存在收斂問題，并且在低噪聲水平下，其定位精度能夠達到克拉美羅界(Cramér-Rao Lower Bound，CRLB)，但是TSWLS并不能很好地適應(yīng)大的測量噪聲。隨后文獻[16]將多維標度(MDS)的算法引入到TDOA/FDOA的定位體制中，文獻[17-18]使用半定松弛(SDR)算法對目標參數(shù)進行估計，在噪聲水平較大時，這些算法性能均優(yōu)于傳統(tǒng)的TSWLS算法。

MDS和SDR算法雖然有良好的性能，但是都需要大量的運算。為了能夠降低計算復(fù)雜度，文獻[19-20]分別提出了兩種改進的TDOA/FDOA運動目標兩步定位算法。這兩種算法的第1步均與傳統(tǒng)TSWLS的第1步相同，主要差異體現(xiàn)在第2步的處理上：文獻[19]中算法的第2步是對第1步定位誤差進行估計并且修正第1步得到的估計結(jié)果；文獻[20]中算法的第2步利用中間變量與目標參數(shù)之間的關(guān)系建立了新的等式并且利用WLS求解。與傳統(tǒng)的TSWLS相比，這兩種算法不僅計算簡便，測量噪聲較大時的定位精度也進一步提高。

傳統(tǒng)的TSWLS及其改進的閉式算法[19-20]雖然都能夠有效提升TDOA/FDOA定位性能。但是通常而言，對于一個3維空間的無模糊定位最少需要4個觀測站，而這些算法至少需要5個觀測站才能達到理論精度。針對這一問題，本文提出了一種基于觀測站數(shù)目最小化的TDOA/FDOA定位算法，這種算法對目標定位只需要4個 觀測站。該算法第1步首先對TSWLS中的偽線性方程重新表示，分離未知參數(shù)，根據(jù)中間變量與目標參數(shù)的關(guān)系，建立了只含有中間變量的多項式并且求解，隨后給出了目標位置以及速度的估計值；算法第2步使用泰勒級數(shù)對中間變量一階展開，將其表示為目標位置與速度估計誤差的線性函數(shù)，得到第1步定位誤差的最佳線性無偏估計量。隨后，基于第2步的誤差估計值，對第1步得到的目標參數(shù)估計值進行修正。本文所提出的算法不存在收斂問題，是一種閉式解；此外本文還推導(dǎo)了當測量噪聲較小時新算法可以得到CRLB精度。仿真實驗不僅證明了使用4個觀測站估計遠場或者近場目標，測量噪聲較小時新算法便能夠?qū)崿F(xiàn)CRLB，同時還表明在額外增加1個 觀測站的條件下新算法比改進的TSWLS算法具有更好的魯棒性。

本文中*°表示測量值*的真實值。O和0分別代表元素均為0的矩陣以及向量。

1 定位場景

目標與第i個觀測站的真實距離為

(1)

(2)

2(si-s1)Tuoi=2,3,…,M

(3)

式中：M-1組曲面的交點給出了目標位置的估計值，但是無法給出目標速度的估計值。由于目標與觀測站之間存在相對運動，因此通過觀測站可以獲得關(guān)于到達頻率差(FDOA)的信息，聯(lián)合TDOA與FDOA信息不僅可以進一步提高定位精度，同時也可以給出目標速度的估計值。

對式(1)關(guān)于時間求導(dǎo)，距離變化率與目標參數(shù)的關(guān)系為

(4)

為了利用FDOA的信息，對式(3)關(guān)于時間求導(dǎo)可以得到一組關(guān)于FDOA的方程

i=2,3,…,M

(5)

(6)

(7)

J=(m-mo)TQ-1(m-mo)T

(8)

2 降維定位算法

2.1 算法第1步

將測量向量r=ro+cΔt代入式(3)并且只保留噪聲的線性項，則式(3)變?yōu)?/p>

(9)

(10)

(11)

i=2,3,…,M

(12)

(13)

Γ1、Γ2和B1、Bt、Bf的表達式分別為

(14)

(15)

(16)

式中：θ的前3個元素為目標的位置估計值，后3個 元素為目標的速度估計值，而W1、α、β和γ的表達式分別為

(17)

式(16)給出了目標位置-速度由中間變量表示的方程，同時表明了求解目標參數(shù)需要已知中間變量。因此，接下來將給出中間變量的求解過程。

(18)

式中：G1,t、h1,t以及Bt在式(12)、式(13)和式(15)中已經(jīng)給出了定義，矩陣Wt為

(19)

α1與β1的表達式分別為

(20)

(21)

(22)

降維算法第1步計算過程需要注意以下幾點。

注意1式(17)加權(quán)矩陣W1的計算與未知參數(shù)有關(guān)。為了解決這個問題，首先令W1=Q-1，通過式(16)得到未知參數(shù)的初始估計值；隨后將初始值代入W1=(B1QB1T)-1得到更新后的加權(quán)矩陣，利用新的加權(quán)矩陣進一步給出更精確的估計值。

注意2二次多項式(21)和式(22)可能會得到多個中間變量解，因此有多組位置-速度估計值，為了解決這個問題，選擇使式(8)的ML代價函數(shù)最小化的解作為最終的估計值。

注意3在三維定位場景下，所有的觀測站不能處于同一個平面，這樣才能夠保證G1列滿秩，進而通過WLS得到唯一解。

2.2 算法第2步

第1步的估計值可以表示為

θ=θo+Δθ

(23)

國際工程教育認證以學(xué)生為中心、產(chǎn)出導(dǎo)向、持續(xù)改進作為三大基本倡導(dǎo)理念。IEET所提出的專業(yè)認證標準可以歸納為專業(yè)培養(yǎng)目標、學(xué)生、教學(xué)成效與評量、課程組成、教師、設(shè)備及空間、行政支持與經(jīng)費、領(lǐng)域認證規(guī)范、持續(xù)改善成效等九大規(guī)范，這些規(guī)范又可劃分為教學(xué)、保障和持續(xù)改進三個方面。其中，培養(yǎng)目標內(nèi)容包括培養(yǎng)目標要求、內(nèi)容和修訂機制[4]，在整個認證過程中起到綱領(lǐng)性作用，反映學(xué)生畢業(yè)5年后所能夠達到的職業(yè)層次和專業(yè)成就。

(24)

(25)

a=(u-s1)/r1

(26)

(27)

式(23)～式(25)代入式(11)，僅保留線性誤差項可以得到新的等式方程為

h2-G2Δθ=cB1n

(28)

(29)

G2,t和G2,f均為(M-1)×3維的矩陣，并且他們的每一行可以分別表示為

(30)

式(28)的解通過加權(quán)最小二乘法給出為

(31)

(32)

2.3 使用接收觀測站個數(shù)的分析

對式(11)進行簡單的變形，可以得到傳統(tǒng)TSWLS的第1步偽線性方程為

h1-G′φo=cB1n

(33)

傳統(tǒng)TSWLS的第1步對中間變量與目標位置同時求解，因此對于傳統(tǒng)的TSWLS來講，要保證式(33)是非欠定方程，在三維定位場景下必須滿足2(M-1)≥2×(3+1)，也就是說至少需要5個觀測站；而本文提出的降維算法先對參數(shù)φo中的中間變量進行求解，隨后在對θo求解，這樣避免了同時求解中間變量與目標位置，達到了降維的目的，進而保證式(11)和式(28)為非欠定方程只需滿足2(M-1)≥2×3，因此本文所提出的算法僅需要4個觀測站。

3 性能分析

本節(jié)給出降維算法的性能分析，首先證明了降維算法給出的是無偏估計值，隨后推導(dǎo)了該算法的協(xié)方差矩陣，并且與CLRB進行了比較，理論證明了在低噪聲水平下，該算法可以達到CLRB。

降維算法中第2步Δθ的估計值可以表示為

(34)

式中：δθ為估計誤差。

將式(23)和式(34)代入式(32)，則最終估計值的期望偏差與協(xié)方差矩陣可以分別表示為-E(δθ)和cov(δθ)。式(31)代入式(34)，得到估計誤差δθ的表達式為

(35)

式(35)中，由于矩陣G2以及n均包含噪聲項，因此偏差期望-E(δθ)與協(xié)方差cov(δθ)并不容易計算。但是在噪聲較小的情況下，G2中測量值的噪聲項可以忽略不計，此時式(35)是誤差項δθ關(guān)于噪聲n的一組線性等式，由于n服從零均值高斯噪聲分布，從而-E(δθ)=0。因此本文所提出的算法在噪聲較小的假設(shè)下可以給出目標參數(shù)的無偏估計值。更進一步，目標參數(shù)估計的協(xié)方差矩陣為

(36)

接下來證明當噪聲較小情況下，本文所提出的算法能夠?qū)崿F(xiàn)運動目標的克拉美羅下界。文獻[21] 已經(jīng)給出運動目標的克拉美羅下界為

(37)

4 仿真實驗

本節(jié)給出了仿真實驗進一步驗證本文所提出算法的定位性能。假設(shè)TDOA與FDOA測量誤差之間相互獨立，則距離差及其變化率測量誤差的協(xié)方差矩陣為Q=blkdiag(Qt,Qf)，其中Qt=σ2R以及Qf=0.01σ2R[20]。這里，σ2為測量噪聲的方差，R是對角線元素均為1，其他元素為0.5的矩陣。仿真通過計算機進行多次蒙特卡羅實驗獲取目標位置和速度估計的均方根誤差指標來評估算法的性能，蒙特卡羅仿真次數(shù)設(shè)置為l=10 000次，均方根誤差(RMSE)指標定義為

(38)

本文考慮兩種仿真情形：場景1在三維場景中僅使用4個觀測站進行定位，場景2在同樣的條件下增加一個觀測站進行定位。其中，在每種場景仿真中均考慮了近場條件以及遠場條件。在場景1中對比算法使用了文獻[9]中的泰勒級數(shù)迭代算法，為了保證其較好的收斂性，初值選取了目標的真實位置與速度；場景2中對比算法包括文獻[13]傳統(tǒng)的TSWLS、文獻[19-20]中改進的TSWLS算法。CRLB作為檢驗估計性能的標準在兩種場景中均給出。

如圖1(a)所示，對于近場目標位置估計，當測量噪聲低于5 dB時，降維算法與泰勒級數(shù)展開算法均能夠達到CRLB。當噪聲達到10 dB時，本文所提出的算法與泰勒級數(shù)展開的算法均開始偏離CRLB，但是泰勒級數(shù)展開的算法受噪聲影響較大，偏離程度較大，產(chǎn)生了“門限效應(yīng)(Threshold Effect)[13]”，無法給出比較準確的估計值。值得一提的是，低噪聲時泰勒級數(shù)展開的算法之所以具有如此好的估計性能，主要是因為初值選取了目標的真實位置，迭代至少進行3次收斂性及估計精度便能夠得到保證。

如圖1(b)所示，近場目標速度的估計與位置估計情況大致相同，低噪聲時兩種算法的估計性能均能達到CRLB。當噪聲達到10 dB時，兩種算法均開始偏離CRLB，但是本文提出的算法能夠更好地適應(yīng)大的測量噪聲，具有更好的魯棒性。

表1 場景1中觀測站的位置與速度Table 1 Positions and velocities of stations in Scenario 1

而對于遠場目標位置與速度的估計，如圖2所示，在測量噪聲較小的時候，兩種算法均能夠達到CRLB，隨著測量噪聲的增加，當噪聲到達-5 dB 時，泰勒級數(shù)展開的算法同樣會產(chǎn)生門限效應(yīng)，而本文提出的算法依舊十分接近CRLB。

值得注意的是，場景1對運動目標參數(shù)的估計僅僅使用4個接收觀測站便可以達到CRLB定位精度，而對于傳統(tǒng)的TSWLS或者是其改進算法還需要至少額外增加一個觀測站才能進行定位求解，本文所提出的算法使用了最少的觀測站數(shù)量。

圖1 使用4個觀測站時近場目標位置和速度估計均方根誤差對比Fig.1 RMSE comparisons of near field position and velocity estimation of in the four stations condition

圖2 使用4個觀測站時遠場目標位置和速度估計均方根誤差對比Fig.2 RMSE comparisons of far field position and velocity estimation in four stations condition

從圖3(a)中可以看出，定位系統(tǒng)存在有5個觀測站時，對于近場目標位置估計，當測量噪聲誤差不超過-10 dB時，4種算法均能夠達到CRLB。當噪聲高于-5 dB時，傳統(tǒng)的TSWLS開始偏離CRLB曲線，這主要是由于其第2步認為噪聲較小時可以忽略二次誤差項，但是當噪聲較大時忽略二次項誤差便會導(dǎo)致定位性能下降。另外，當噪聲高于10 dB時兩種改進的TSWLS算法逐漸偏離CRLB；當測量噪聲到達15 dB時，使用本文算法得到的位置估計均方根誤差比其他兩種改進的TSLWS算法提高了約1 dB。

從圖3(b)可以看出，對于近場目標速度的估計，在低噪聲水平除了傳統(tǒng)TSWLS算法，其他算法均可以達到CRLB。當測量噪聲水平較高時，對比算法性能逐漸下降。其中，文獻[19]中的算法在噪聲高于0 dB時明顯開始偏離CRLB，而當噪聲高于5 dB時，文獻[20]中的算法同樣發(fā)生偏離。本文提出的降維算法即使在較大噪聲水平下，也十分接近CRLB。當測量噪聲到達10 dB時，降維算法得到的速度估計均方根誤差比文獻[19] 的算法降低了4.5 dB，比文獻[20]的算法減少了1.5 dB。仿真實驗證明了本文所提出的定位算法具有更好的魯棒性。

而對于遠場目標位置與速度的估計，如圖4所示；在低噪聲水平下，傳統(tǒng)的TSWLS、兩種改進的TSWLS算法以及本文提出的算法均能夠?qū)崿F(xiàn)CRLB，給出良好的定位性能。而隨著噪聲的增大，傳統(tǒng)的TSWLS和文獻[19-20]中兩種改進TSWLS算法均會相繼偏離CRLB,無法給出比較精確的目標位置以及速度估計值,而本文提出的算法依舊十分接近CRLB,仿真實驗表明了本文所提出的算法能夠更好地適應(yīng)大的測量噪聲。

圖3 使用5個觀測站時近場目標位置和速度估計均方根誤差對比Fig.3 RMSE comparisons of near field position and velocity estimation in five stations condition

圖4 使用5個觀測站時遠場目標位置和速度估計均方根誤差對比Fig.4 RMSE comparisons of far field position and velocity estimation in five stations condition

5 算法計算量分析

本節(jié)主要對所提出降維算法的運算量進行了分析，為了分析簡便，主要以實數(shù)乘法進行分析；并且分別與傳統(tǒng)的TSWLS[13]和文獻[19-20]中兩種改進TSWLS計算量進行了對比。

首先分析本文算法第1步的計算量：

接下來給出第2步計算量的分析過程：

此外，由于需要對加權(quán)矩陣初始化，第1步中2)和4)各需要計算2次，因此本文提出的算法一共需要27M3+51M2+363M+889次的實數(shù)乘法運算，文獻[19]給出了改進TSWLS的計算量為48M3-48M2+768M+7 212，與傳統(tǒng)TSWLS需要48M3-48M2+768M+8 010次實數(shù)乘法相比，運算量相當；通過計算，文獻[20]改進的TSWLS算法需要進行48M3-48M2+768M+4 012次實數(shù)乘法,本文所需要的實數(shù)乘法次數(shù)小于文獻[20]中所需的實數(shù)乘法，因此本文算法運算量更小。

6 結(jié) 論

1) 在目標處于運動狀態(tài)的三維無源定位場景下，針對傳統(tǒng)TSWLS及其改進的閉式算法精確定位至少需要5個觀測站的問題，本文聯(lián)合TDOA以及FDOA，提出一種最小化觀測站數(shù)目的無源定位算法。該閉式算法僅需要4個觀測站即可進行定位，理論分析表明該算法在低噪聲水平下能夠達到CRLB。

2) 仿真實驗證明了算法的有效性：僅使用4個 觀測站，測量噪聲較小時該算法便可以達到CRLB精度；若額外增加一個觀測站，相比與傳統(tǒng)的TSWLS以及改進算法，在不增加算法計算量的前提下，該算法在測量噪聲較大時能提供更精確的目標參數(shù)估計值。

附錄A:

首先當測量噪聲較小時有如下條件成立

(A1)

(A2)

(A3)

(A4)

(A5)

根據(jù)式(15)和式(30)分別對B1以及G2的定義，-L可以進一步表示為

(A6)

i=2,3,…,M

(A7)

利用條件1和2，分別有如下近似

(A8)

(A9)

利用式(A8)和式(A9)，式(A7)可以近似為

i=2,3,…,M

(A10)

i=2,3,…,M

(A11)

同樣利用條件1和2，可以分別得到如下近似等式

(A12)

(A13)

利用(A12)與式(A13)，式(A11)可以近似表示為