一種4站情況下基于TDOA/FDOA的無源定位方法

國強，李文韜

哈爾濱工程大學 信息與通信工程學院，哈爾濱 150001

無源定位是一種在自身不發射電磁波的情況下利用被探測目標自身所發出的信號進行定位的目標探測技術，目前基于多站的無源定位已廣泛應用于軍事和民用領域[1-3]，包括傳感器網絡[4]、無線通信[2]、雷達[5]、導航等[6-8]。多站無源定位可以利用的參數包括到達角(Angel of Arrival，AOA)、到達時間(Time of Arrival，TOA)、到達時間差(Time Difference of Arrival，TDOA)、到達頻率差(Frequency Difference of Arrival，FDOA)等。而聯合多種測量參數的定位體制可以融合不同參數的優勢，提升對目標輻射源信號類型的適應能力，并在一定程度上提高定位精度。本文主要研究在無約束條件下的三維空間中，聯合TDOA以及FDOA信息對運動輻射源進行定位的解算方法。在無約束條件下的三維空間中，傳統的針對基于TDOA/FDOA的多站無源定位解算方法的研究大多都是在5個及以上接收站的條件下進行的，對在4個接收站的情況下基于TDOA/FDOA的定位解算方法的研究較少，因此本文將重點對4站情況下基于TDOA/FDOA的定位解算方法進行研究。

基于TDOA/FDOA的定位方程是非線性的，因此求解過程復雜且困難，對于如何求解這一問題已經進行了大量的研究，目前對TDOA/FDOA方程進行解算的方法主要有迭代法、解析法、搜索法。

泰勒級數法[9]的主要思想是在初始值處進行泰勒級數展開并忽略高階項，然后通過迭代進行求解。文獻[10-11]通過采用半正定松弛方法將非凸的時頻差定位問題松弛為一個凸的SDP問題, 然后再利用內點法等方法進行求解。這些方法都是典型的迭代法，迭代法的抗噪性能較好，但迭代法對初始值的選擇較為敏感，初始值誤差較大時容易出現發散。

文獻[12]提出了一種經典的解析法，通過引入輔助變量將非線性方程轉化為偽線性方程，之后利用兩步加權最小二乘(Two-Stage Weighted Least Squares，TSWLS)法求解方程，該算法在噪聲較低的時候能夠快速對目標進行精確定位。在此基礎上文獻[13]考慮了接收站的站址誤差，對TSWLS方法進行了改進，在站址誤差較大時該方法的定位精度要高于原TSWLS方法。文獻[14]認為在噪聲較大時，TSWLS第2步中等式誤差所包含的二階誤差項并不能夠被忽略，因此舍棄了TSWLS方法的第2步，直接利用泰勒級數展開法對TSWLS第1步得到的結果進行修正。作為文獻[12]的擴展，文獻[15]提出了約束加權最小二乘(Constrained Weighted Least Squares，CWLS)方法，與TSWLS方法相比該方法具有更好的抗噪性能，但它是一種迭代法，需要初始值且算法復雜度遠遠高于TSWLS方法。與迭代法相比，解析法運算速度快且不需要初值，但是以上基于偽線性方程的解析法都至少需要5個接收站才能對目標進行定位，而且定位精度受噪聲影響較大。目前對4站情況下的解析法的研究很少，文獻[16]提出了一種可以在4站情況下實現定位的解析法，該方法通過聯合干擾變量的定義式與時頻差方程的加權最小二乘(Weighted Least Squares，WLS)解構造出一組新的方程來對TSWLS方法的第1步進行改進，改進后的方法在 4站情況下也可以得到一個解，之后利用文獻[14]中的泰勒級數修正法對這個初始的解進行修正，但是這種方法在第2步進行修正的過程中忽略了二階以上的誤差項，因此定位精度受噪聲影響較大，在噪聲較高時定位精度較低。

隨著智能尋優算法的發展，搜索法也逐漸地被應用到TDOA/FDOA定位當中，文獻[17]利用改進的布谷鳥算法(Modified Cuckoo Search，MCS)來對目標進行定位，在噪聲較高的情況下取得了優于TSWLS方法的定位效果。文獻[18]則先是在4站情況下利用改進的粒子群算法得到一個初始值，之后利用TSWLS方法的第2步對初始值進行修正，其定位效果較直接使用粒子群算法進行搜索有明顯的提升。傳統的搜索法均可以在只有4個接收站的情況下對目標進行定位，不需要由解析法提供初值，而且其抗噪性能普遍要優于解析法，但在搜索過程中收斂速度較慢導致搜索法實時性很低[12]，因此如何在4站的情況下得到一種高效的搜索法是一個值得研究的問題。

針對在只有4個接收站的情況下搜索法所存在的收斂速度慢導致實時性很低的問題，本文提出了一種基于改進的加權最小二乘(Modified Weighted Least Squares，MWLS)法與螢火蟲算法(Firefly Algorithm，FA)相結合的無源定位方法(MWLS-FA)。首先利用文獻[16]第1步中的改進的加權最小二乘法在4站情況下得到一個初始估計值，之后利用這個初始估計值為FA的搜索區域進行限制，同時本文也在約束條件的添加和參數選擇兩個方面針對性地對FA方法做出了調整和改進，最后利用FA方法得到最終的估計值。仿真實驗不僅證明了本文提出的MWLS-FA方法在4站情況下對目標的定位精度可以達到克拉美羅下限(Cramer-Rao Lower Bound，CRLB)，而且相比于傳統的搜索法以及文獻[16]中提出的方法，MWLS-FA方法具有更好的抗噪性能，而且在實時性上較傳統的搜索法有明顯的提升。同時也證明了在5個接收站的情況下MWLS-FA方法的定位性能要優于TSWLS和CWLS方法。

1 TDOA/FDOA定位模型

(1)

(2)

(3)

TDOA定位方程只能夠估計出目標的位置而不能估計其速度，而且僅基于TDOA信息的定位方程可能不足以對移動目標的位置進行精確的定位。而當輻射源和接收站之間存在相對的運動時，基于TDOA/FDOA的定位算法可以同時估計移動目標的位置和速度。

對式(1)取時間導數可以得到距離變化率的表達式為

(4)

對式(3)取時間導數可得FDOA方程為

(5)

r=ro+cΔt

(6)

(7)

(8)

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

(9)

2 改進的加權最小二乘法

將包含噪聲的TDOA測量值代入式(3)中并忽略二階誤差項可以得到

(10)

(11)

(12)

式中：

(13)

利用加權最小二乘法對θ進行估計可得

(14)

式中：W為權重矩陣，文獻[12]給出了W的表達式為

(15)

(16)

W=Q-1

(17)

α=(GTWG)-1GTWh1

(18)

β=(GTWG)-1GTWh2

(19)

γ=(GTWG)-1GTWh3

(20)

式中：

(21)

僅利用式(9)得到的加權最小二乘解為

(22)

αt=(GtTWtGt)-1GtTWtht1

(23)

(24)

式中：

(25)

設Qt為所測TDOA數據所含噪聲的協方差矩陣，則

Wt=B-TQt-1B-1

(26)

(27)

(28)

針對式(27)和式(28)可能存在多個解導致最終得到多個位置-速度估計值的情況，選擇使式(9) 的代價函數最小化的解作為最終的估計值。

在MWLS-FA方法中FA方法需要一個目標函數。因此本文根據式(12)得到基于MWLS方法的代價函數，并將其作為FA方法的目標函數,即

f(θ)=(h-Gθ)TW(h-Gθ)

(29)

式中的θ受約束于：

(30)

(31)

此時基于TDOA/FDOA的定位問題可以轉化為約束優化問題，接下來本文將利用FA方法對該約束優化問題進行求解。

3 螢火蟲算法

螢火蟲算法是Yang提出的一種仿生智能優化算法[20]，其仿生原理為：利用搜索空間中的點來模擬自然界中的螢火蟲個體，將搜索和尋優的過程模擬成螢火蟲的吸引和移動過程，同時利用待解問題的目標函數來衡量螢火蟲個體所處位置的優劣。該算法具有參數少、概念簡單等優點，已經在優化問題[21]、圖像壓縮處理[22]、聚類[23]等多種領域得到應用，并在這些領域中取得了良好的效果。

在螢火蟲算法中，有3個理想化的假設為：

1) 螢火蟲不分性別，只會被更亮的螢火蟲所吸引。

2) 對于任意兩只螢火蟲，亮度較暗的那只螢火蟲會朝著較亮的那只移動。吸引力決定了移動的距離，而且吸引力會隨著距離的增加而減小。最亮的螢火蟲會隨機選擇方向進行移動。

3) 螢火蟲的亮度由其所在位置所對應的目標函數決定。

如上所述，螢火蟲算法主要包含兩個要素，螢火蟲的亮度和吸引度，亮度決定螢火蟲的移動方向，吸引度決定螢火蟲的移動距離，螢火蟲不斷通過向較亮的螢火蟲移動來完成優化，直到達到預定條件為止。假設第i和第j只螢火蟲的位置分別為xi=[xi1,xi2,…,xiD]、xj=[xj1,xj2,…,xjD]，rij為第i只螢火蟲到第j只螢火蟲的距離，則rij的表達式為

(32)

螢火蟲i相對螢火蟲j的亮度定義為

I(rij)=I0e-γrij

(33)

式中:I0為螢火蟲的初始亮度，由式(12)得到的目標函數值決定，目標函數值越優初始亮度越高。γ光強吸收系數，為一個常數。rij為螢火蟲i和j的距離。

螢火蟲i相對螢火蟲j的吸引力定義為

β(rij)=βoe-γrij2

(34)

式中：β0表示最大吸引度，一般取值為1，表示光源處(r=0)螢火蟲的吸引度。由式(34)可知，吸引度會隨著距離的增加和吸收光強系數的增大而減小。

螢火蟲i被螢火蟲j吸引向其移動的位置更新公式為

xi(t+1)=xi(t)+β(rij)(xj(t)-xi(t))+ε·αi

(35)

式中：αi為步長因子;ε是在[-0.5,0.5]上服從均勻分布的隨機因子。不同亮度的螢火蟲隨機分布在D維空間中。螢火蟲的亮度和吸引力分別由式(32)和式(33)計算。較暗的螢火蟲會向較亮的螢火蟲移動。為了避免陷入局部最優,將擾動項ε·αi添加到位置更新的過程。最后，螢火蟲將聚集在最高亮度的螢火蟲周圍，從而得到最優結果。

4 MWLS-FA方法

FA方法雖然在很多領域內得到了應用，但其在基于TDOA/FDOA的無源定位解算問題中的應用仍然存在著諸多需要解決的問題。首先目標的位置和速度是未知的，因此難以為FA方法提供一個合適的搜索區域；其次利用式(29)得到的目標函數受到約束條件的限制，對約束條件的處理也是一個需要研究的問題；此外，FA方法中的參數的選擇也是一個亟待解決的問題。為了解決這些問題，本文在搜索區域、約束條件處理、參數選擇3個方面針對性地對FA方法做出了改進和調整。

4.1 螢火蟲算法的搜索區域

搜索法的定位精度較高，而且在4站情況下可以實現定位，但是搜索法需要在全局范圍內進行搜索，因此相對于解析法和迭代法要花費更多的時間，如果對搜索區域進行合理的限制，那么包括FA方法在內的搜索法可以在很大程度上減少搜索所需的時間。本節利用MWLS方法所獲得的初始結果為中心確定一個搜索范圍，為FA方法提供一個有限的搜索區域，這可以在很大程度上提高FA方法的收斂速度，減少搜索時間。

當FA方法的搜索區域過小時，目標真實位置和速度值可能不在搜索區域內從而影響算法的精度，而當搜索范圍過大時，又會對FA方法的收斂速度產生影響，甚至可能會導致算法不收斂。因此，搜索范圍的選擇應由第1步得到的估計值與目標位置和速度的真實值之間的差值所決定，當搜索范圍的上限在各個維度上都大于初始值在各個維度上加上差值，而下限在各個維度上都小于初始值在各個維度上減去差值時，可以在理論上將真實值包含在搜索范圍內。而第1步得到的目標估計位置和速度與目標真實位置和速度之間的差值是不確定的，所以一個固定的搜索范圍無法同時保證FA方法在不同情況下的實時性和收斂性。

影響第1步所得到的估計值與真實值的差值的因素有很多，比如噪聲功率、目標距離各個接收站的平均距離、布站方式、運動目標相對于各個接收站的速度等，而根據實際的仿真情況可以得出其中兩個影響最大的因素分別為噪聲功率的大小和目標距離各個接收站的平均距離。而在實際的定位過程中，無法得知真實目標到各個接收站的平均距離，因此采用由第1步 得到的目標位置到各個接收站的距離取平均值作為代替，第1步得到目標位置雖然跟真實位置相比有一定的誤差，但可以大致反映出目標到各接收站的遠近情況。綜上所述，本文在計算搜索范圍時，將考慮噪聲功率的大小和第1步 得到的目標位置到各個接收站的平均距離對搜索范圍的影響。

假設va∈[-20,20] dB、s∈[0,500] m，分別對va和s在各自的區間上取20個值，并利用MWLS方法得到相對應的d值，這樣就得到了400組數據。通過這些數據對二元函數d=h(va,s)進行擬合可以得到：

h(va,s)=(a+exp(b×va+c))×(s3+q)

(36)

式中:a=2.2×10-7、b=0.36、c=-19.5、q=-4 000。 對函數v=g(va,s)進行擬合可以得到：

g(va,s)=(k+exp (l×va+m))×(s3+n)

(37)

式中:a=1.22×10-7、b=0.36、c=-19.4、q=9 500，擬合函數d=h(va,s)、v=g(va,s)的圖像如圖1所示。

圖1 MWLS方法得到的估計值與真實值的差值與 s和va之間的函數關系Fig.1 Functional relationship among difference between estimated value obtained by MWLS method and real value, va and s

(38)

(39)

(40)

(41)

利用式(38)～式(41)得到的上下界在將真實值包含在搜索范圍內的同時也保證了算法的實時性。此外，這種求解搜索范圍上下限的方法不光適用于FA方法，對其他搜索法也同樣適用。

4.2 約束條件的處理

(42)

(43)

在處理約束優化問題時，等式約束通常需要轉換為不等式約束，于是θ在以上2個約束條件上的約束違反程度可表示為

Gi(x)=max{|gi(x)|-δ,0}

(44)

式中：i=1,2；δ為等式約束的容忍參數。目前在約束優化問題中應用最為廣泛的為罰函數法，相較于傳統的罰函數法，自適應罰函數法是一種通用性更強的約束處理方法，因為它可以根據從螢火蟲種群進化過程中獲取的信息動態的調整懲罰系數，本文采用了文獻[24]提出的一種自適應的罰函數法，其適應度函數的形式為

(45)

式中：λ(t)在每一代中按如下方式更新：

(46)

式中：t為當前迭代次數；β1=β2=2；初始λ值為1;case 1表示在過去g(由用戶自己定義)代中找到的最好個體均為可行解;case 2表示在過去g代中找到的最好個體均為不可行解，其原理可以理解為：若在之前的更新過程中找到的最好個體均為可行解，則說明懲罰系數已足夠大，可適當降低懲罰系數；若此前找到的最好個體均為不可行解，則說明懲罰系數過小，需增大懲罰系數。

4.3 螢火蟲算法的參數選擇

FA方法的性能優劣主要取決于全局搜索和局部搜索之間的平衡，在搜索過程初期，算法應具有較強的全局搜索能力以跳出局部最優，而在搜索過程后期，算法應具有更強的局部搜索能力以在小范圍內尋找更為精確的最優解。步長因子α對這種平衡起著至關重要的作用，由于FA在初始階段需要盡量對更多的區域進行搜索，此時α越大越好，而在算法的后期，為了保證收斂性，需要步長越小越好，因此理論上α應隨著迭代次數的增加而減少。本節在FA算法中加入了一種計算方法來控制步長，即

αt+1=αt×(1-Δαt)

(47)

式中：t表示迭代次數；Δαt=1-(10-4/0.9)1/t。可以看出式(47)中步長因子隨著迭代次數的增長而非線性遞減。

值得注意的是，本文提出的MWLS-FA方法對搜索范圍進行了合理的限制，因此無需再對除步長因子外的參數進行額外的動態處理。從文獻[25] 可知，對于幾何結構較為簡單的目標函數，步長因子α的最優取值范圍為α∈(0,0.5]，光吸收系數γ的最優取值范圍為γ∈(0,1]，而在α和γ設置合理的情況下，當迭代次數大于500次時，幾何結構較為簡單的目標函數求解精度基本穩定。因此本文的MWLS-FA方法初始步長因子α設置為0.5，光吸收系數γ設置為1，最大迭代次數設置為500次。

綜上所述，MWLS-FA方法的步驟為

MWLS-FA方法第1步:1.令W=Q-1、Wt=Qt-1。2.將式(27)和式(28)得到的ro1和r·o1代入到式(14)中計算θ^。3.利用得到的估計值θ^計算B1,并根據式(15)和式(26)重新計算W和Wt。4.重復以上兩個步驟2～3次,即可得到一個初始估計值θ^。第2步:1.令目標函數為式(29)。2.在由式(38)～式(41)所限定的區域內隨機產生一組初始的螢火蟲xi=[xi1,xi2,…,xi6],其中i=1,2,…,N。3.設置最大迭代次數Tmax=500,當前迭代次數t=0。4.whilet

5 仿真實驗

本節通過進行實驗仿真對本文所提出的定位算法的性能進行驗證。MWLS-FA方法中螢火蟲個數為40個，初始步長因子為0.5，光強吸收系數為1，最大迭代次數為500次。蒙特卡羅循環次數設置為10 000次，為TDOA和FDOA測量值添加均值為零的高斯白噪聲，其中TDOA的噪聲的方差為δd2，協方差矩陣為Rt=δd2σ，其中σ為對角線上元素為1、其他元素為0.5的矩陣。FDOA的噪聲方差為TDOA噪聲方差的0.1倍[12]，其協方差矩陣為Rf=0.1δd2σ，利用均方根誤差RMSE進行性能分析，表達式為

(48)

本文考慮2種仿真場景，場景1為在三維環境下使用4個接收站對目標進行定位，而場景2則在同樣的條件下使用5個接收站對目標進行定位， 在每個場景中又分別在近場和遠場條件下對目標進行定位。在4站情況下使用文獻[17]中的MCS算法和文獻[16]中提出的降維算法進行對比，在5站情況下使用文獻[12]中的TSWLS方法、文獻[15]中的CWLS方法進行對比，同時在兩種場景中均使用CRLB作為檢驗估計性能的標準。計算機仿真條件為Windows10，64位操作系統，Core i5-9300H處理器，16 GB內存，MATLAB9.0版本。

圖2給出了在僅有4個接收站的近場情況下， 估計精度隨著噪聲誤差的增加的變化情況。如圖2所示，在只有4個接收站的近場情況下，3種 方法在噪聲較低時均能達到CRLB。MCS算法和降維算法在噪聲高于0 dB時開始偏離CRLB， 而MWLS-FA方法則是在測量噪聲高于4 dB的時候開始偏離CRLB，且偏離程度要小于MCS算法和降維算法。因此在只有4個接收站的近場情況下，3種方法在低噪聲的情況下都可以達到CRLB，而相對于其他2種方法，MWLS-FA方法在測量噪聲較高時具有更好的魯棒性。

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

在4站情況下對遠場目標的位置和速度進行估計的結果如圖3所示，從圖3中可以看出，3種方法在遠場情況下的定位精度均不如在近場情況下的定位精度，當噪聲小于-12 dB時MCS算法和降維算法均可達到CRLB，而當噪聲大于-12 dB時，MCS算法和降維算法開始偏離CRLB，MWLS-FA方法在噪聲大于-8 dB時開始偏離CRLB，而且偏離程度要小于其他2種方法，可以看出在4站遠場情況下，相較于其他2種方法，MWLS-FA方法也具有更好的抗噪性能。

圖2 在近場情況下使用4個接收站對目標的位置及速度進行估計的均方根誤差Fig.2 RMSE required for estimation of position and velocity of near target in 4-receiver case

圖3 在遠場情況下使用4個接收站對目標的位置及速度進行估計的均方根誤差Fig.3 RMSE required for estimation of position and velocity of far target in 4-receiver case

圖4是在有5個接收站的情況下對近場目標的定位性能分析。從圖4中可以看出，3種方法在低噪聲的情況下均能達到CRLB，而TSWLS方法受噪聲影響最大，在噪聲到達0 dB后便開始偏離CRLB，CWLS方法在4 dB的時候開始偏離CRLB，MWLS-FA方法則在10 dB時開始偏離CRLB，而且在噪聲較大的時候也十分貼近CRLB，抗噪聲性能要優于其他2種方法。

圖4 在近場情況下使用5個接收站對目標的位置及速度進行估計的均方根誤差Fig.4 RMSE required for estimation of position and velocity of near target in 5-receiver case

圖5給出了在5個接收站情況下對遠場目標定位的性能分析。如圖5所示， 3種算法在遠場情況下的估計精度均不如近場情況，TSWLS方法和CWLS方法分別在噪聲高于-10 dB和-6 dB時開始偏離CRLB，MWLS-FA方法在-4 dB時開始偏離CRLB，不過偏離程度相較于其他兩種方法較小，這表明了MWLS-FA方法在該場景下相對于TSWLS和CWLS方法也具有更好的抗噪聲能力。

圖5 在遠場情況下使用5個接收站對目標的位置及速度進行估計的均方根誤差Fig.5 RMSE required for estimation of position and velocity of far target in 5-receiver case

接下來對上述方法的估計時間進行分析，由于各種定位方法對近場目標和遠場目標進行估計所需的時間均非常接近，因此本節接下來將只在近場情況下對各種方法的估計時間進行分析。

如表2所示，在只有4個接收站的情況下，降維方法的平均估計時間為0.002 s，MCS算法的平均估計時間約為9.46 s，而MWLS-FA方法的平均估計時間則為2.24 s，可以看出因為MWLS-FA方法是一種搜索法，有著一個搜索過程，而降維方法是一種閉式解析方法，沒有搜索或者迭代的過程，因此MWLS-FA方法的估計時間要高于降維方法。而相對于同樣采用搜索策略的MCS方法，MWLS-FA由于對搜索區域進行了有效的限制，因此在很大程度上加快了算法的收斂速度，節省了估計時間，大幅度提高了搜索法的實時性。

如表3所示，在5個接收站的情況下，TSWLS方法是一種閉式解析方法，沒有搜索和迭代過程，因此估計時間很短，平均估計時間約為0.002 s。CWLS方法的平均估計時間約為3.44 s。

表2 估計時間在4站情況下隨噪聲的變化

表3 估計時間在5站情況下隨噪聲的變化

而MWLS-FA方法的平均估計時間約為2.24 s，可以看出，由于MWLS-FA方法是一種搜索法，有一個搜索的過程， 因此實時性不如TSWLS方法，但平均估計時間要少于CWLS方法。

6 結 論

1) 在基于時頻差的運動目標無源定位系統中，針對傳統方法在4個接收站的情況下難以對目標進行實時及精確定位的問題，提出了一種基于改進的加權最小二乘方法與螢火蟲算法相結合的時頻差無源定位方法，該方法在4個接收站的情況下對目標的定位精度可以達到CRLB，而且在實時性和抗噪性方面要優于傳統的搜索法。

2) 仿真實驗證明了MWLS-FA方法的有效性和實時性，該方法在4站情況下對目標的定位精度可以達到CRLB，而且在實時性和抗噪性方面較傳統的搜索法有較大的提升。結果還表明，在5站情況下，該算法在抗噪性方面要優于TSWLS方法和CWLS方法。

[21] GANDOMI A H, YANG X S, ALAVI A H. Mixed variable structural optimization using firefly algorithm[J]. Computers & Structures, 2011, 89(23-24): 2325-2336.