一種最大似然網格搜索衛星FDOA輻射源定位算法

張 威，吳 佟，張更新

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

0 引言

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

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

1 FDOA定位方程

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

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

2 算法描述

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

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

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

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

為了優化搜索范圍，減少搜索點數，加快搜索收斂的速度，本文采用一種分步搜索的方法。

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

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

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

第二步，以K1為基準點，建立新的搜索范圍D1，這里D1滿足式（7）所示關系:

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

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

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

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

式中，

式中，F是與約束有關的未知參數的梯度矩陣，當F為零向量時，式（14）退化為式（10），對于有地球約束的FDOA定位，，CRLB是任何無偏參數估計均方根誤差的下限。

5 算法性能仿真分析

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

表1 4顆觀測衛星的星歷（1 Jul 2011 12:00:00.000）

表2 FDOA測量時刻衛星、輻射源及地面接收站在ECEF坐標系下的位置矢量

表3 FDOA測量時刻衛星在ECEF坐標系下的速度矢量

表4列出了FDOA測量誤差的標準差σ分別為10－4Hz、10－3Hz、10－2Hz、10－1Hz及1Hz時，基于網格搜索的最大似然算法的RMSE。其中，RMSE均為5000次獨立實驗的統計結果，網格搜索算法每次進行10等分割網格化，搜索結果處FDOA向量值與實際測量FDOA向量值的誤差值小于10－7Hz時停止搜索。圖2給出了基于網格搜索的最大似然算法與CRLB曲線的比較。

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

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

由表4及圖2可見，基于網格搜索的最大似然算法在所給的FDOA測量誤差情況下都能較好的接近克拉美羅下限。仿真結果中FDOA測量誤差的標準差σ較低時網格搜索算法的RMSE與CRLB之間的偏離較大，總結其主要原因為仿真中一些固定參數的選取誤差，如地球半徑、地球橢圓偏心率、坐標系轉換中的春分點角等參數的選取誤差，隨著測量誤差的增大，這些參數選取誤差的影響越來越小。

6 結束語

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

[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，FORSSELL，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.