一种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.