張 威，吳 佟，張更新

（1.解放軍理工大學(xué)通信工程學(xué)院，江蘇南京210007;2.總參信息化部駐714軍代室，江蘇南京210016）

0 引言

隨著衛(wèi)星通信業(yè)務(wù)的飛速發(fā)展，衛(wèi)星通信面臨的電磁環(huán)境日益惡化，難以避免受到各種輻射源有意或無意的干擾，因此對輻射源進(jìn)行準(zhǔn)確地?zé)o源定位以采取有效的針對措施有著重要的意義。現(xiàn)有衛(wèi)星輻射源定位體制中，采用多顆衛(wèi)星接收地面輻射源信號到達(dá)時(shí)差（TDOA，Time Difference of Arrival）進(jìn)行定位的技術(shù)已經(jīng)相對成熟，主要的算法有平面相交[1－2]、球面相交[3]、球面內(nèi)插[4]、泰勒級數(shù)[5]、最小二乘[6－7]和粒子濾波[8－9]，相應(yīng)的算法已經(jīng)投入到實(shí)際應(yīng)用;而基于信號到達(dá)頻差（Frequency Difference of Arrival，F(xiàn)DOA）的定位方程較為復(fù)雜，但隨著定位參數(shù)測量技術(shù)的發(fā)展，F(xiàn)DOA定位技術(shù)成為目前衛(wèi)星無源定位技術(shù)研究的熱點(diǎn)之一。

FDOA定位方程組一般具有非線性的特點(diǎn)，難以給出其解析解，而采用迭代算法對其求解則需要設(shè)置初始位置?；诰W(wǎng)格的最大似然算法并不直接求解FDOA的定位方程，而是將輻射源可能存在的區(qū)域進(jìn)行網(wǎng)格劃分，將劃分所得的網(wǎng)格點(diǎn)位置的FDOA參數(shù)與測量所得參數(shù)進(jìn)行匹配，并將匹配誤差最小的網(wǎng)格點(diǎn)作為下次搜索的區(qū)域中心，縮小搜索區(qū)域，再次進(jìn)行網(wǎng)格劃分搜索和匹配，直到得到滿足誤差門限的網(wǎng)格點(diǎn)為止，并將此網(wǎng)格點(diǎn)作為最終的輻射源估計(jì)位置。網(wǎng)格搜索算法是針對實(shí)際問題求解非線性方程的一種有效方法，具有算法過程簡單、精度高和頑健性強(qiáng)等特性。

1 FDOA定位方程

設(shè)輻射源在ECEF（Earth Centered Earth Fixed）坐標(biāo)系下的位置矢量為，由于衛(wèi)星的星歷已知，可以計(jì)算出各衛(wèi)星在ECEF坐標(biāo)系下的位置矢量為，速度矢量為，i=1，…，M，M為觀測衛(wèi)星的數(shù)目。則輻射源信號到達(dá)衛(wèi)星的多普勒頻移如式（1）所示:

式中，fc為輻射源信號載波頻率，c為信號傳播速度，假設(shè)地球半徑為R，將地球面的約束條件引入FDOA定位方程組，則定位方程如式（2）所示，ΔF1i為FDOA參數(shù)，ΔF1i=Δfi－Δf1。

2 算法描述

由式（2）可知頻差的方程如式（3）所示:

將所有的FDOA測量值ΔF1i組成測量向量，如式（4）所示:

則地球面約束的FDOA最大似然網(wǎng)格位置所搜方程如式（5）所示:

式中，ML表示最大似然搜索解，D為搜索的范圍，為將搜索的位置向量代入式（3）所得的FDOA向量，arg表示取變量，這里變量為。式（5）的基本含義是搜索范圍D內(nèi)的所有點(diǎn)，搜索使得最小的位置向量作為最終目標(biāo)位置。為了使對范圍D的搜索可實(shí)現(xiàn)，必須在范圍D內(nèi)取有限的點(diǎn)進(jìn)行搜索，這里采用網(wǎng)格化的取點(diǎn)方法，具體方法如下。

為了優(yōu)化搜索范圍，減少搜索點(diǎn)數(shù)，加快搜索收斂的速度，本文采用一種分步搜索的方法。

第一步，確定初始搜索范圍D0，如式（6）所示:

式中，λ與φ分別為地球表面的經(jīng)度與緯度，λ1與λ2分別為觀測衛(wèi)星共視區(qū)在地球面的經(jīng)度的最小與最大值，φ1與φ2分別為共視區(qū)緯度的最小與最大值。對初始搜索范圍D0進(jìn)行基于經(jīng)度和緯度的N等分割網(wǎng)格化，形成N2網(wǎng)格，則共有N+（）12經(jīng)緯度網(wǎng)格點(diǎn)。當(dāng)取N=10時(shí)，搜索范圍D0的網(wǎng)格劃分結(jié)果如圖1所示。將此網(wǎng)格點(diǎn)代入式（5），搜索與測量FDOA向量誤差值最小的網(wǎng)格點(diǎn)，設(shè)此點(diǎn)為K1。

圖1 將初始搜索范圍N等分割網(wǎng)格化示意圖（N=10）

第二步，以K1為基準(zhǔn)點(diǎn)，建立新的搜索范圍D1，這里D1滿足式（7）所示關(guān)系:

式中，λd與φd分別取D0中每個(gè)網(wǎng)格的經(jīng)度與維度跨度的2倍。然后搜索范圍D1進(jìn)行基于經(jīng)度和緯度的N等分割網(wǎng)格化，這里取N=10，對分割后的網(wǎng)格點(diǎn)進(jìn)行再次搜索，得到與測量FDOA向量誤差值最小的網(wǎng)格點(diǎn)K2。如此進(jìn)行L步，直到搜索結(jié)果處FDOA向量與實(shí)際測量FDOA向量的誤差值小于規(guī)定門限時(shí)，將此時(shí)的搜索結(jié)果作為輻射源的最終估計(jì)位置。

3 均方根定位誤差與克拉美羅下限

為了能更好的評價(jià)算法性能，引入2個(gè)概念，分別為均方根定位誤差（Root Mean Square Error，RMSE）與FDOA定位算法的克拉美羅下限[10－12]（Cramér－Rao Lower Bound，CRLB），表達(dá)式如式（8）和式（9）所示:

式（10）中，JFDOA為FDOA定位算法的Fisher信息矩陣，如式（11）所示:

式中，

式中，F(xiàn)是與約束有關(guān)的未知參數(shù)的梯度矩陣，當(dāng)F為零向量時(shí)，式（14）退化為式（10），對于有地球約束的FDOA定位，，CRLB是任何無偏參數(shù)估計(jì)均方根誤差的下限。

5 算法性能仿真分析

為便于仿真，選取3顆高軌觀測衛(wèi)星的星歷如表1所示，輻射源位于廣州（東經(jīng)113.3°、北緯23.1°和高程0km），地面接收站位于北京（東經(jīng)116.4°、北緯39.9°和高程0km），設(shè)FDOA測量時(shí)刻為1Jul 2011 12:00:00.000，則經(jīng)過計(jì)算可得，衛(wèi)星、輻射源及地面接收站該時(shí)刻在ECEF坐標(biāo)系中的位置矢量如表2所示，衛(wèi)星的速度矢量如表3所示。

表1 4顆觀測衛(wèi)星的星歷（1 Jul 2011 12:00:00.000）

表2 FDOA測量時(shí)刻衛(wèi)星、輻射源及地面接收站在ECEF坐標(biāo)系下的位置矢量

表3 FDOA測量時(shí)刻衛(wèi)星在ECEF坐標(biāo)系下的速度矢量

表4列出了FDOA測量誤差的標(biāo)準(zhǔn)差σ分別為10－4Hz、10－3Hz、10－2Hz、10－1Hz及1Hz時(shí)，基于網(wǎng)格搜索的最大似然算法的RMSE。其中，RMSE均為5000次獨(dú)立實(shí)驗(yàn)的統(tǒng)計(jì)結(jié)果，網(wǎng)格搜索算法每次進(jìn)行10等分割網(wǎng)格化，搜索結(jié)果處FDOA向量值與實(shí)際測量FDOA向量值的誤差值小于10－7Hz時(shí)停止搜索。圖2給出了基于網(wǎng)格搜索的最大似然算法與CRLB曲線的比較。

表4σ不同取值情況下各算法的RMSE/km

圖2 基于網(wǎng)格搜索的最大似然算法與CRLB曲線比較圖

由表4及圖2可見，基于網(wǎng)格搜索的最大似然算法在所給的FDOA測量誤差情況下都能較好的接近克拉美羅下限。仿真結(jié)果中FDOA測量誤差的標(biāo)準(zhǔn)差σ較低時(shí)網(wǎng)格搜索算法的RMSE與CRLB之間的偏離較大，總結(jié)其主要原因?yàn)榉抡嬷幸恍┕潭▍?shù)的選取誤差，如地球半徑、地球橢圓偏心率、坐標(biāo)系轉(zhuǎn)換中的春分點(diǎn)角等參數(shù)的選取誤差，隨著測量誤差的增大，這些參數(shù)選取誤差的影響越來越小。

6 結(jié)束語

針對FDOA定位方程的非線性，求解過程較為復(fù)雜的特點(diǎn)，首先針對FDOA定位方程組進(jìn)行了分析，然后詳細(xì)研究了其基于網(wǎng)格的最大似然算法的求解過程，并推導(dǎo)了其克拉美羅下線。經(jīng)過理論和仿真分析，獲得了如下結(jié)論:該算法求解過程較為簡便，收斂速度較快，能夠有效地逼近CRLB，是最優(yōu)的定位估計(jì)器，可以用于實(shí)際工程，為研究FDOA定位的相關(guān)人員提供一定的參考。

[1]SCHMIDT R.A New Approach to Geometry of Range Difference Location[J].IEEE Trans Aerospace Electron，1972，8（11）:821－835.

[2]SCHMIDT R.Least Square Range Difference Location[J].IEEE Trans Aerospace and Electronic Systems，1996，32（1）:234－242.

[3]SHAU H C，ROBINSON A Z.Passive Source Localization Employing Intersecting Spherical Surfaces From Time－of－arrival Differences[J].IEEE Trans on ASSP，1987，35（8）:1123－1225.

[4]SMITH J O，ABEL J S.Closed－form Least－squares Source Location Estimation from Range－difference Measurements[J].IEEE Trans on ASSP，1987，35（12）:1661－1669.

[5]FOY W K.Position－location Solutions by Taylor－series Estimation[J].IEEE Trans Aerospace and Electronic Systems 1976，12（2）:187－194.

[6]CHAN Y T，HO K C.A Simple and Efficient Estimator for Hyperbolic Location[J].IEEE Trans Signal Processing，1994，42（8）:1905－1915.

[7]HO K C，CHAN Y T.Geolocation of a Known Altitude Object from TDOA and FDOA Measurements[J].IEEE Trans Aerospace and Electronic Systems，1997，33（3）:770－783.

[8]GUSTAFSSON F.Positioning Using Time-difference of Arrival Measurements[C]//Proc.IEEE Conf Acoustics，Speech and Signal Processing（ICASSP），Hong Kong，China，2003:553－556.

[9]GUSTAFSSON F，BERGMAN N，F(xiàn)ORSSELL，et al.Particle Filter for Positioning，Navigation，and Tracking[J].IEEE Trans Signal Processing，2002，50（2）:425－437.

[10]GORMAN J D，HERO A O.Lower Bounds on Parametric Estimators with Constraints[C]//4th ASPP Workshop Spectrum Estimation Modeling，Minneapolis，MN，1988:223－228.

[11]GORMAN J D，HERO A O.Lower Bounds for Parametric Estimation with Constraints[J].IEEE Transactions on Information Theory，1990，26（6）:1285－1301.

[12]MARZETTA T L.A Simple Derivation of the Constrained Multiple Parameter Cramér－Rao Bound[J].IEEE Transactions on Acoustics，1993，41（6）:2247－2249.