一种最大似然网格搜索卫星FDOA辐射源定位算法

张 威，吴 佟，张更新

（1.解放军理工大学通信工程学院，江苏南京210007;2.总参信息化部驻714军代室，江苏南京210016）

0 引言

随着卫星通信业务的飞速发展，卫星通信面临的电磁环境日益恶化，难以避免受到各种辐射源有意或无意的干扰，因此对辐射源进行准确地无源定位以采取有效的针对措施有着重要的意义。现有卫星辐射源定位体制中，采用多颗卫星接收地面辐射源信号到达时差（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.