基于到达时差(time difference of arrival, TDOA)的无源定位，利用目标辐射的电磁波到达不同观测站的到达时间差，再乘以电磁波传播速度，可以得到距离差，从而绘制出多组以观测站为焦点、距离差为长轴的双曲线，通过多个双曲线交点实现定位。该方法又被称为双曲线定位法。

定位方程解算方法的研究是一个十分经典的问题, 根据TDOA定位方程高度非线性的特点, 解决该问题大致可以通过3类解算方法。第1种方法是迭代法, 该方法可以得到数值解, 其典型代表是泰勒级数法^[1]。该方法需要一个迭代初始值, 在该值处进行一阶泰勒展开并进行循环迭代, 将每次迭代的误差与预先设置的阈值比较, 当小于阈值时，即得到最终结果。文献[2]基于约束总体最小二乘的思想, 利用牛顿法迭代得到数值解。文献[3]设立多个不同的参考站, 用泰勒展开法求解TDOA定位的最大似然问题。第2种方法是解析法, 该方法可以得到闭式解, 球面插值法^[4]、球面相交法^[5]以及平面相交法^[6]是该方法早期的典型代表。不同于迭代法, 其最大特点是计算量小, 无需迭代初值, 但由于没有修正过程, 其精度较差。在此基础上, 两步加权最小二乘^[7]法被提出, 该方法在第一步引入冗余变量线性化TDOA方程, 估计出位置的粗略值; 该方法在第二步利用辐射源位置和冗余变量的关系修正估计值, 称该算法为Chan算法。文献[8]针对加权非线性最小二乘法需要进行复杂矩阵运算的问题, 提出一种非线性期望最大化算法, 在取得与原始算法接近性能的同时大大减少了计算量。对于运动目标定位, 针对Ho等人提出的经典两步加权最小二乘^[9]存在阈值效应的问题, 文献[10]、文献[11]对最小二乘式进行不同的改进, 定位性能在较高噪声水平下有了一定提高。第3种方法可以描述为基于最大似然目标函数的优化算法。粒子群算法^[12-13]利用粒子种群对位置的多次迭代得到每个个体的局部最优解, 进而通过搜索获得全局最优解, 精度较高。另外，遗传算法^[14]、蚁群算法^[15]等也逐渐被引入, 解决了该最大似然优化问题。文献[16]利用电磁波到达方向(direction of arrival, DOA)和TDOA的信息融合来更好地估计辐射源位置, 该方法先利用DOA得到更加准确的时差, 再利用Chan算法进行位置解算。

以上算法不考虑站址误差, 但无人机集群是高机动、高灵活性平台, 往往存在较大站址误差。在利用TDOA建模求解的过程中, 需要考虑误差因素。在文献[9]的基础上, Ho研究了考虑观测站位置误差的前提下, 基于TDOA^[17]以及TDOA/到达频差(frequency difference of arrival, FDOA)^[18]的算法改进, 将位置误差项加入到加权最小二乘式中, 使定位精度得到提高。此后, 文献[19]、文献[20]对该闭式算法进行了改进, 提高了噪声鲁棒性。文献[21]在联合时频差对移动辐射源进行位置速度估计基础上, 引入了多普勒频移, 通过仿真证明小噪声时该算法可以接近克拉美罗下界(Cramer-Rao lower bound, CRLB)。针对传感器位置随机且难以精确估计的问题, Ho提出了通过放置校准发射器来提高传感器的位置估计精度的方法^[22], 推导出校准器的最佳位置准则。文献[23]基于时频差测量, 利用半定规划法求解, 并将定位场景扩展到了非视距传输和多径衰落等场景。文献[24]模拟设置了多种不同的校准发射器, 通过测量TDOA并开发了一种半定规划算法，有效地降低了站址误差对定位精度的影响。在此基础上, 文献[25]进一步将时钟同步偏差考虑进来, 扩大了算法的应用场景。文献[26]提出一种基于加权最小二乘的闭式高精度定位方法, 该方法无需迭代运算, 计算量较小。文献[27]提出基于到达时间(time of arrival, TOA)变化率的高精度定位方法, 忽略径向加速度，进行有偏估计。

本文首先基于无人机集群对单目标定位进行建模, 研究TDOA单目标定位方程的解算方法, 在考虑站址误差的前提下改进Chan算法和泰勒算法, 改进后的算法分别称为解析法和迭代法。将粒子群算法的解作为泰勒迭代的初值, 提出粒子群泰勒协同算法, 仿真实验证明该方法具有较高的定位精度。

在TDOA定位中, 辐射源到达每个辅站与主站的时间会存在一组时间差, 将时间差乘以电磁波传播速度即得到距离差, 每组距离差可确定一组双曲线, 多组双曲线的交点即为辐射源的位置。

二维定位场景原理如图 1所示。要实现单目标定位，至少需要两组双曲线, 即至少需要3架无人机。设3架无人机位于s₁=(x₁, y₁)^T、s₂=(x₂, y₂)^T、s₃=(x₃, y₃)^T处。把s₁设为主机, 则另外两架飞机分别可与主机确定一组双曲线, 辐射源就位于双曲线交点上。为了去除虚假定位点并提高定位精度, 往往需要增加无人机的数量。

三维空间单目标定位模型如图 2所示, 每架辅机都可与主机确定一组三维空间中的双曲面, 在不存在TDOA观测误差的情况下, 至少需要4架无人机才能确定单目标的位置, 要想去除虚假点, 需要进一步增加无人机的数量。

在图 2中，假设共有M架无人机, 目标位置坐标为u =(x, y, z)^T, 主机真实位置为s₁^o=(x₁^o, y₁^o, z₁^o)^T, 辅机的真实坐标分别为s₂^o=(x₂^o, y₂^o, z₂^o)^T, s₃^o=(x₃^o, y₃^o, z₃^o)^T, s₄^o=(x₄^o, y₄^o, z₄^o)^T, …, s_M^o=(x_M^o, y_M^o, z_M^o)^T。设第i架无人机到辐射源的距离为r_i^o(i=1, 2, …, M), 则有

(1)$\begin{gathered}r_i^o=\left|\boldsymbol{u}-\boldsymbol{s}_i^o\right|= \\\sqrt{\left(x-x_i^o\right)^2+\left(y-y_i^o\right)^2+\left(z-z_i^o\right)^2}\end{gathered}$

辅机i与主机接收到信号的实际时间差为

(2)$t_{i 1}^o=\frac{r_{i 1}^o}{c}=\frac{r_i^o-r_1^o}{c}$

测量到的时间差带有随机噪声, 即t_i1=t_i1^o+Δt_i1。其中，Δt_i1为时差测量噪声, t_i1为带有误差的时间差测量值, 且r_i1=t_i1×c, 是存在误差的距离差测量值。其中，c为光速，c=3×10⁸ m/s。令真实距离差向量为r^o=[r₂₁^o, r₃₁^o, …, r_M1^o]^T, 则

(3)$\begin{gathered}\boldsymbol{r}=\boldsymbol{r}^o+\boldsymbol{v}= \\{\left[r_{21}^o, r_{31}^o, \cdots, r_{M 1}^o\right]^{\mathrm{T}}+\left[n_{21}, n_{31}, \cdots, n_{M 1}\right]^{\mathrm{T}}}\end{gathered}$

式中: n_i1=c·Δt_i1; v为由到达时间差测量误差导致的距离差误差向量。误差的协方差矩阵Q =E[vv^T]。

另外, 作为高机动平台, 无人机自身位置一般是由全球定位系统测量的, 准确性和实时性难以保证, 自身位置会存在误差, 则定位过程中所用的带有误差的站址为

(4)$\boldsymbol{s}_i=\boldsymbol{s}_i^o+\Delta \boldsymbol{s}_i$

式中: Δs_i为无人机位置的误差向量。

将不考虑站址误差的Chan算法和Taylor算法直接应用到具有较大站址误差的无人机集群定位场景中时, 会出现定位精度严重偏离CRLB的情况。根据文献[28], 当无人机高低角为30°, 飞行高度为5 km时, 卫星对其位置估计的误差约为30 m。根据文献[29], 由于近年来卫星定位精度的提高, 无人机在低空慢速飞行时, 误差可控制在10 m以内。本节将站址误差考虑进来, 将改进后的Chan算法称为解析法^[29], 改进后的Taylor算法称为迭代法^[29], 并针对Taylor算法需要较为精确的初始值的问题, 提出了一种粒子群泰勒协同的算法。

解析法定位模型为第1节阐述的三维场景, 该算法可得出无人机具有位置误差时的闭式解。对距离关系式r_i1^o+r₁^o=r_i^o两侧进行平方，得到(r_i1^o+r₁^o)²=r_i^o2, 再将式(1)代入该式, 有

(5)$r_{i 1}^{o 2}+2 r_{i 1}^o r_1^o=R_i^{o 2}-R_1^{o 2}-2\left(s_i^o-s_1^o\right)^{\mathrm{T}} \boldsymbol{u}$

式中: $R_i^o=\sqrt{x_i^{o 2}+y_i^{o 2}+z_i^{o 2}} $。无人机真实三维坐标是无法得到的; 由于存在TDOA测量误差, 真实距离差r_i1^o也是无法得到的。因此，式(5)中的R_i^o与r_i1^o实际上无法得到, 求解时只能使用带有误差的量R_i和r_i1。通过引入额外变量r₁^o，使得等式可以看作辐射源位置u的线性方程, 设第一次加权最小二乘的待求向量为u₁=[x, y, z, r₁^o]^T, 在TDOA测量值和无人机自身站址存在误差的情况下, 将式(5)中的真实值替换为含有噪声的测量值, 则式(5)的误差向量为

(6)$\boldsymbol{\varepsilon}_1=\boldsymbol{h}_1-\boldsymbol{G}_1 \boldsymbol{u}_1$

式中:

(7)$\left\{\begin{array}{l}\boldsymbol{h}_1=\left[\begin{array}{c}r_{21}^2-R_2^2+R_1^2 \\r_{31}^2-R_3^2+R_1^2 \\\vdots \\r_{M 1}^2-R_M^2+R_1^2\end{array}\right] \\\boldsymbol{u}_1=\left[\begin{array}{llll}x & y & z & r_1^o\end{array}\right]^{\mathrm{T}} \\\boldsymbol{G}_1=-2\left[\begin{array}{cc}\left(\boldsymbol{s}_2-\boldsymbol{s}_1\right)^{\mathrm{T}} & r_{21} \\\left(\boldsymbol{s}_3-\boldsymbol{s}_1\right)^{\mathrm{T}} & r_{21} \\\vdots & \vdots \\\left(\boldsymbol{s}_M-\boldsymbol{s}_1\right)^{\mathrm{T}} & r_{21}\end{array}\right]\end{array}\right.$

式中: $\boldsymbol{s}_i=\boldsymbol{s}_i^o+\Delta \boldsymbol{s}_i ; R_i^2=x_i^2+y_i^2+z_i^2 $。式(6)中u₁的解为

(8)$\boldsymbol{u}_1=\left(\boldsymbol{G}_1^{\mathrm{T}} \boldsymbol{W}_1 \boldsymbol{G}_1\right)^{-1} \boldsymbol{G}_1^{\mathrm{T}} \boldsymbol{W}_1 \boldsymbol{h}_1$

式中: $ \boldsymbol{W}_1=\mathrm{E}\left[\boldsymbol{\varepsilon}_1 \boldsymbol{\varepsilon}_1^{\mathrm{T}}\right]^{-1}=\left[\boldsymbol{B}_1 \boldsymbol{Q}_t \boldsymbol{B}_1^{\mathrm{T}}+\boldsymbol{D}_1 \boldsymbol{Q}_p \boldsymbol{D}_1^{\mathrm{T}}\right]^{-1}$。其中, $\boldsymbol Q_t=c^2 \mathrm{E}\left[\Delta \boldsymbol{t} \Delta \boldsymbol{t}^{\mathrm{T}}\right], \boldsymbol{Q}_p=\mathrm{E}\left[\Delta \boldsymbol{s} \Delta \boldsymbol{s}^{\mathrm{T}}\right], \boldsymbol{B}_1=2 \operatorname{diag}\left\{r_2^o, r_3^o, \cdots, r_M^o\right\} $, 且

(9)$\boldsymbol{D}_1=2\left[\begin{array}{cccc}-\left(\boldsymbol{u}-\boldsymbol{s}_1^{\boldsymbol{o}}\right)^{\mathrm{T}} & \left(\boldsymbol{u}-\boldsymbol{s}_2^{\boldsymbol{o}}\right)^{\mathrm{T}} & \bf{0}^{\mathrm{T}} & \bf{0}^{\mathrm{T}} \\\vdots & \bf{0}^{\mathrm{T}} & \ddots & \vdots \\-\left(\boldsymbol{u}-\boldsymbol{s}_1^{\boldsymbol{o}}\right)^{\mathrm{T}} & \bf{0}^{\mathrm{T}} & \cdots & \left(\boldsymbol{u}-\boldsymbol{s}_{\mathrm{M}}^{\boldsymbol{o}}\right)^{\mathrm{T}}\end{array}\right]$

式中：0是一个3×1的列向量。然而，B₁和D₁中均包含无误差的位置和距离, 实际中无法得到, 只能取近似值, 近似方法在文献[17]、文献[18]的W₁初始值获取中已有说明。为了获得W₁初始值, 在目标位置远离无人机集群时, 有r₁≈ r₂≈…≈r_M, 则B₁可近似为2r₁[I], D₁可近似为2r₁[－1P], 1是大小为(M－1)×3的矩阵, P是大小为(M－1)×3(M－1)的对角矩阵, 定义为

(10)$\boldsymbol{P}=\left[\begin{array}{c}(1, 1, 1) \\(1, 1, 1) \\\ddots \\(1, 1, 1)\end{array}\right]$

加权矩阵W₁近似写为

(11)$\boldsymbol{W}_1=\frac{1}{4}\left(\boldsymbol{Q}_t+\left[\begin{array}{ll}-\bf{1} & \boldsymbol{P}\end{array}\right] \boldsymbol{Q}_p\left[\begin{array}{ll}-\bf{1} & \boldsymbol{P}\end{array}\right]^{\mathrm{T}}\right)^{-1}$

将式(11)代入式(8)中，可得到初始目标位置u₁, 再通过计算B₁和D₁得到性能更优的W₁, 进行第一次加权最小二乘。大量仿真结果表明^[18], 定位精度对加权矩阵中的噪声不敏感, 所以对初始加权矩阵W₁近似是可行的。u₁的协方差矩阵由下式给出：

(12)$\operatorname{cov}\left(\boldsymbol{u}_1\right)=\left(\boldsymbol{G}_1^{\mathrm{T}} \boldsymbol{W}_1 \boldsymbol{G}_1\right)^{-1}$

接下来，进行第二次加权最小二乘

(13)$\boldsymbol{\varepsilon}_2=\boldsymbol{h}_2-\boldsymbol{G}_2 \boldsymbol{u}_2$

式中:

(14)$\left\{\begin{array}{l}\boldsymbol{h}_2=\left[\begin{array}{l}\left(u_1(1)-x_1\right)^2 \\\left(u_1(2)-y_1\right)^2 \\\left(u_1(3)-z_1\right)^2 \\u_1(4)^2\end{array}\right], \\\boldsymbol{G}_2=\left[\begin{array}{lll}1 & 0 & 0 \\0 & 1 & 0 \\0 & 0 & 1 \\1 & 1 & 1\end{array}\right] \\\boldsymbol{u}_2=\left[\left(x-x_1\right)^2, \left(y-y_1\right)^2, \left(z-z_1\right)^2\right]^{\mathrm{T}}\end{array}\right.$

则式(13)的解为

(15)$\boldsymbol{u}_2=\left(\boldsymbol{G}_2^{\mathrm{T}} \boldsymbol{W}_2 \boldsymbol{G}_2\right)^{-1} \boldsymbol{G}_2^{\mathrm{T}} \boldsymbol{W}_2 \boldsymbol{h}_2$

式中:

(16)$\boldsymbol{W}_2=\left[\boldsymbol{B}_2 \operatorname{cov}\left(\boldsymbol{u}_1\right) \boldsymbol{B}_2^{\mathrm{T}}+\boldsymbol{D}_2 \boldsymbol{Q}_p \boldsymbol{D}_2\right]^{-1}$

(17)$\left\{\begin{array}{c}\boldsymbol{B}_2=2\left[\begin{array}{cccc}x-x_1 & 0 & 0 & 0 \\0 & y-y_1 & 0 & 0 \\0 & 0 & z-z_1 & 0 \\0 & 0 & 0 & r_1^o\end{array}\right] \\\boldsymbol{D}_2=2\left[\begin{array}{cccc}\bf{0}^{\mathrm{T}} & \bf{0}^{\mathrm{T}} & \cdots & \bf{0}^{\mathrm{T}} \\\bf{0}^{\mathrm{T}} & \bf{0}^{\mathrm{T}} & \cdots & \bf{0}^{\mathrm{T}} \\\left(\boldsymbol{u}-\boldsymbol{s}_1^o\right)^{\mathrm{T}} & \bf{0}^{\mathrm{T}} & \boldsymbol{\cdots} & \bf{0}^{\mathrm{T}}\end{array}\right]\end{array}\right.$

B₂和D₂中的无人机位置和辐射源位置可以将第一步加权最小二乘解代入。根据u₂与待求辐射源位置关系, 并在消除符号模糊后, 求得最终的目标位置为

(18)$\boldsymbol{u}=\boldsymbol{P}\left[\sqrt{u_2(1)}, \sqrt{u_2(2)}, \sqrt{u_2(3)}\right]^{\mathrm{T}}+\boldsymbol{s}_1$

式中: $ \boldsymbol{P}=\operatorname{diag}\left\{\operatorname{sgn}\left(\boldsymbol{u}_1(1: 3)-\boldsymbol{s}_1\right)\right\}$。

解析法在Chan算法的基础上将站址误差考虑进来, 得到的定位精度有一定提升, 但与Chan算法类似, 解析法在求解过程中做了两次近似, 因此精度的提升有限。

在Taylor法的基础上, 把无人机的站址误差考虑进来, 所得到的改进算法称为迭代法^[29]。令向量θ=[u^T, s₁^oT, …, s_M^oT]^T为真实目标位置和无人机位置的向量, 且

(19)$\boldsymbol{f}(\boldsymbol{\theta})=\boldsymbol{T}=\boldsymbol{M}-\boldsymbol{E}$

式中: T=[r^oT, s^oT]^T, r^o=[r₂₁^o, r₃₁^o, …, r_M1^o]^T为真实的距离差向量, 且s^o=[s₁^oT, s₂^oT, …, s_M^oT]^T为无人机实际位置。M=[r₂₁, r₃₁, …, r_M1, s₁^T, s₂^T, …, s_M^T]^T为测量向量, E=[cΔt₂₁, cΔt₃₁, …, cΔt_M1, Δs₁^T, Δs₂^T, …, Δs_M^T]^T为误差向量。假设E服从正态分布, 则E的协方差矩阵为

(20)$\boldsymbol{Q}=\left[\begin{array}{ll}\boldsymbol{Q}_t & \boldsymbol{O}^{\mathrm{T}} \\\boldsymbol{O} & \boldsymbol{Q}_p\end{array}\right]$

式中: O是一个3M×(M-1)的零矩阵。Q_t=c²E[ΔtΔt^T], Q_p=E[ΔsΔs^T], 并且Δs=[Δs₁^T, Δs₂^T, …, Δs_M^T]^T。

设置向量θ_g=[u_g^T, s_{1, g}^T, …, s_{M, g}^T]^T为迭代的初始值, 将其代入式中, 并忽略二阶以上的误差, 得到

(21)$\left.\boldsymbol{f}(\boldsymbol{\theta})\right|_{\boldsymbol{\theta}=\boldsymbol{\theta}_g}+\boldsymbol{A} \delta \boldsymbol{\theta} \simeq \boldsymbol{M}-\boldsymbol{E}$

式中：δθ为一个3M+3维的列向量, 其值为

(22)$\delta \boldsymbol{\theta}=\left[\delta \boldsymbol{u}^{\mathrm{T}}, \delta \bf{s}_1^{\mathrm{T}}, \cdots, \delta \bf{s}_M^{\mathrm{T}}\right]^{\mathrm{T}}=\boldsymbol{\theta}-\boldsymbol{\theta}_g$

A是一个(4M-1)×(3M+3)的矩阵, 其计算方法为

(23)$\boldsymbol{A}=\left.\left[\begin{array}{ll}\frac{\partial \boldsymbol{f}(\boldsymbol{\theta})}{\partial \boldsymbol{u}} & \frac{\partial \boldsymbol{f}(\boldsymbol{\theta})}{\partial \boldsymbol{s}}\end{array}\right]\right|_{\boldsymbol{\theta}=\boldsymbol{\theta}_{\boldsymbol{g}}}$

计算r_{i, g}^o时, 将真实辐射源位置替换为估计的辐射源位置。通过计算得到A为

(24)$\boldsymbol{A}=\left[\left.\begin{array}{ll}\frac{\partial \boldsymbol{r}^o}{\partial \boldsymbol{u}} & \frac{\partial \boldsymbol{r}^o}{\partial \boldsymbol{s}} \\\frac{\partial \boldsymbol{s}^o}{\partial \boldsymbol{u}} & \frac{\partial \boldsymbol{s}^o}{\partial \boldsymbol{s}}\end{array}\right|_{\boldsymbol{\theta}=\boldsymbol{\theta}_{\boldsymbol{g}}}\right.$

式中:

(25)$\left.\frac{\partial \boldsymbol{r}^o}{\partial \boldsymbol{u}}\right|_{\boldsymbol{\theta}=\boldsymbol{\theta}_g}=\left[\begin{array}{c}\frac{\left(\boldsymbol{u}_g-\boldsymbol{s}_{2, g}\right)^{\mathrm{T}}}{r_{2, g}^o}-\frac{\left(\boldsymbol{u}_g-\boldsymbol{s}_{1, \theta}\right)^{\mathrm{T}}}{r_{1, g}^o} \\\frac{\left(\boldsymbol{u}_g-\boldsymbol{s}_{3, g}\right)^{\mathrm{T}}}{r_{3, g}^o}-\frac{\left(\boldsymbol{u}_g-\boldsymbol{s}_{1, \theta}\right)^{\mathrm{T}}}{r_{1, g}^o} \\\vdots \\\frac{\left(\boldsymbol{u}_g-\boldsymbol{s}_{M, g}\right)^{\mathrm{T}}}{r_{M, g}^o}-\frac{\left(\boldsymbol{u}_g-\boldsymbol{s}_{1, g}\right)^{\mathrm{T}}}{r_{1, g}^o}\end{array}\right]\\\begin{aligned}& \left.\frac{\partial \boldsymbol{r}^o}{\partial \boldsymbol{s}}\right|_{\boldsymbol \theta=\boldsymbol \theta_\boldsymbol g}= \\& {\left[\begin{array}{cccc}\frac{\left(\boldsymbol{u}_g-\boldsymbol{s}_{1, \theta}\right)^{\mathrm{T}}}{r_{1, g}^o} & \frac{-\left(\boldsymbol{u}_g-\boldsymbol{s}_{2, \theta}\right)^{\mathrm{T}}}{r_{2, g}^0} & \cdots & \bf{0}^{\mathrm{T}} \\\frac{\left(\boldsymbol{u}_g-\boldsymbol{s}_{1, \theta}\right)^{\mathrm{T}}}{r_{1, g}^{\circ}} & \bf{0}^{\mathrm{T}} & \cdots & \bf{0}^{\mathrm{T}} \\\vdots & \vdots & \ddots & \vdots \\\frac{\left(\boldsymbol{u}_g-\boldsymbol{s}_{1, \theta}\right)^{\mathrm{T}}}{r_{1, g}^{\circ}} & \bf{0}^{\mathrm{T}} & \cdots & -\frac{\left(\boldsymbol{u}_g-\boldsymbol{s}_{M, g}\right)^{\mathrm{T}}}{r_{M, g}^o}\end{array}\right]} \\& \left.\frac{\partial s^o}{\partial \boldsymbol u}\right|_{\boldsymbol \theta=\boldsymbol \theta_\boldsymbol g}=\boldsymbol{O} \\& \left.\frac{\partial \boldsymbol{s}^o}{\partial \boldsymbol{s}}\right|_{\boldsymbol \theta=\boldsymbol \theta_\boldsymbol g}=\boldsymbol{I} \\&\end{aligned}$

式中: 0是一个3×1的列向量; O是3M×3的零矩阵; I是3M的单位矩阵, 式(21)可重新写为

(26)$\boldsymbol{A} \delta \boldsymbol{\theta}=\boldsymbol{W}-\boldsymbol{E}$

式中:

(27)$\boldsymbol{W}=\boldsymbol{M}-\left.\boldsymbol{f}(\boldsymbol{\theta})\right|_{\boldsymbol \theta=\boldsymbol \theta_\boldsymbol \theta}$

当加权矩阵为Q^－1时, δθ的加权最小二乘估计为

(28)$\delta \boldsymbol{\theta}=\left[\boldsymbol{A}^{\mathrm{T}} \boldsymbol{Q}^{-1} \boldsymbol{A}\right]^{-1} \boldsymbol{A}^{\mathrm{T}} \boldsymbol{Q}^{-1} \boldsymbol{W}$

得到的解等于

(29)$\boldsymbol{\theta}^{(1)}=\boldsymbol{\theta}_g+\delta \boldsymbol{\theta}$

接下来，代入泰勒迭代的算法中, 即有当$ \left\|\delta \boldsymbol{\theta}^{(i)}\right\|>\varepsilon$

(30)$\left\{\begin{array}{l}\delta \boldsymbol{\theta}^{(i+1)}=\left[\boldsymbol{A}^{(i) T} \boldsymbol{Q}^{-1} \boldsymbol{A}^{(i)}\right]^{-1} \boldsymbol{A}^{(i) \mathrm{T}} \boldsymbol{Q}^{-1} \boldsymbol{W}^{(i)} \\\boldsymbol{\theta}^{(i+1)}=\boldsymbol{\theta}^{(i)}+\delta \boldsymbol{\theta}^{(i+1)} \\i=i+1\end{array}\right.$

式中: ε为控制迭代的预设阈值。

根据式(30)进行循环迭代, 直到误差小于预设的终止值, 停止迭代, 目标的估计位置即为向量θ的前3个值, 且目标位置后的每3个值就是一架飞机的位置估计。

由文献[30]、利用两步定位法，将在传感器获取的目标距离信息中获得的目标位置估计作为初值启发, 粒子群算法属于先根据最大似然原理设定目标函数、再利用粒子群搜索目标函数最优值的方法。将其用在时差定位方程的求解时, 关键是找到非线性最大似然目标函数。

粒子群算法的计算过程如下。在D维搜索空间中, 限定搜索范围, 并生成N个粒子, 每个粒子所在位置都可能是问题的解, 粒子通过自身惯性、个体经验和群体经验三大因素得到每次迭代的速度更新值, 进而得到本次迭代的位置更新公式。设第i个粒子当前位置为x_i=[x_{i, 1}, x_{i, 2}, …, x_{i, D}]^T(i=1, 2, …, N)，飞行速度是[v_{i, 1}, v_{i, 2}, …, v_{i, D}]^T(i=1, 2, …, N), 该粒子通过不断学习更新目前搜索到的个体的最优解为p_i=[p_{i, 1}, p_{i, 2}, …, p_{i, D}]^T(i=1, 2, …, N)。经过一定次数的学习迭代后, 每个粒子都到达一个自身最优位置, 称为个体最优解。由于群体的社会性, 粒子之间可进行信息共享, N个粒子的种群就可以从中选出群体最优解, 使得目标函数最优。设全局最优解为p_g, p_g=[p_{g, 1}, p_{g, 2}, …, p_{g, D}]^T(i=1, 2, …, N)。当迭代终止时, p_g即为粒子群算法求出的最优解。粒子的更新公式为

(31)$\left\{\begin{array}{l}\boldsymbol{v}_i^{t+1}=\omega \boldsymbol{v}_i^t+c_1 \boldsymbol{r}_1\left(\boldsymbol{p}_i^t-\boldsymbol{x}_i^t\right)+c_2 \boldsymbol{r}_2\left(\boldsymbol{p}_g^t-\boldsymbol{x}_i^t\right) \\\boldsymbol{x}_i^{t+1}=\boldsymbol{x}_i^t+\boldsymbol{v}_i^t \\i=1, 2, \cdots, N\end{array}\right.$

式中: ω称为惯性权重, 表征上代速度对当代速度的影响程度；c₁和c₂为学习因子, 分别表征自身经验和社会经验对下一步动作的影响程度；r₁和r₂是对角线元素，服从[0, 1]均匀分布的对角阵。式(31)中的两个式子分别是速度更新公式和位置更新公式。令算法经过t次迭代后停止。

在时差定位中, 若不存在测量误差, 则每个观测值确定的双曲线最终交于一点, 可得到精确定位结果, 但存在观测误差时, 双曲线不能交于一点, 可以利用误差概率密度函数得出最大似然式, 并将其作为目标函数, 对目标函数利用粒子群算法求解。假设TDOA的测量噪声为高斯噪声, 服从n_{i, 1}~N(0, δ²)分布:

(32)$\begin{gathered}\prod\limits_{i=2}^N \frac{1}{\sqrt{2 \pi} \delta} \exp \left\{-\frac{\boldsymbol{r}-\left(\boldsymbol{r}_i-\boldsymbol{r}_1\right)}{2 \delta^2}\right\}= \\\left(\frac{1}{\sqrt{2 \pi} \delta}\right)^{N-1} \exp \left\{-\frac{\left[\boldsymbol{r}-\left(\boldsymbol{r}_i-\boldsymbol{r}_1\right)\right]^{\mathrm{T}}\left[\boldsymbol{r}-\left(\boldsymbol{r}_i-\boldsymbol{r}_1\right)\right]}{2 \delta^2}\right\}\end{gathered}$

式中: r的定义如式(3)所示, r_i=[r₂, r₃, …, r_M]^T, r₁=[r₁, r₁, …, r₁]^T。观察式(32), 若使该函数最大, 则等价于求下式的最小值

(33)$(\hat{x}, \hat{y}, \hat{z})=\arg \min\limits _{(x, y, z)}\left\{\left[\boldsymbol{r}-\left(\boldsymbol{r}_i-\boldsymbol{r}_1\right)\right]^{\mathrm{T}}\left[\boldsymbol{r}-\left(\boldsymbol{r}_i-\boldsymbol{r}_1\right)\right]\right\}$

式(32)进一步可等价为

(34)$(\hat{x}, \hat{y}, \hat{z})=\min \left\{\sum\limits_{i=2}^m\left[\boldsymbol{r}-\left(\boldsymbol{r}_i-\boldsymbol{r}_1\right)\right]^2\right\}$

考虑站址误差时, 解析法依旧存在较大误差, 迭代法在精度上有较大优势, 但需要较为精确的初始值, 否则在噪声较大时存在不收敛的问题。这时，解析法的结果不能再作为迭代法的初值代入, 应考虑将一种精度较高算法的结果作为迭代初值, 粒子群算法符合这一要求, 其代价函数需要在存在站址误差时进行修改。原本的代价函数为

(35)$(\hat{x}, \hat{y}, \hat{z})=\min \left\{\sum\limits_{i=2}^m\left[\boldsymbol{r}-\left(\boldsymbol{r}_i-\boldsymbol{r}_1\right)\right]^2\right\}$

式中: r=r^o+v=[r₂₁^o, r₃₁^o, …, r_M1^o]^T+[n₂₁, n₃₁, …, n_M1]^T; r_i= [r₂, r₃, …, r_M]^T; r₁=[r₁, r₁, …, r₁]^T。根据t_i1^o=r_i1^o/c=(r_i^o－r₁^o)/c, r为带有测量误差的距离差向量, 与时差测量误差有关, 向量r_i需要知道无人机的确切位置。因此, 需要对其进行修改。测量得到带有站址误差的无人机位置s_i=s_i^o+Δs_i=[x_i, y_i, z_i]^T, 则令

(36)$\begin{gathered}\tilde{r}_i=\left|\boldsymbol{u}-\boldsymbol{s}_i\right|= \\\sqrt{\left(x-x_i\right)^2+\left(y-y_i\right)^2+\left(z-z_i\right)^2}\end{gathered}$

并且$ \tilde{\boldsymbol{r}}_i=\left[\tilde{r}_2, \tilde{r}_3, \cdots, \tilde{r}_M\right]^{\mathrm{T}}, \tilde{\boldsymbol{r}}_1=\left[\tilde{r}_1, \tilde{r}_1, \cdots, \tilde{r}_1\right]^{\mathrm{T}}$, 则代价函数可修改为

(37)$(\hat{x}, \hat{y}, \hat{z})=\min \left\{\sum\limits_{i=2}^m\left[\boldsymbol{r}-\left(\tilde{\boldsymbol{r}}_i-\tilde{\boldsymbol{r}}_1\right)\right]^2\right\}$

存在站址误差时, 改进的Chan算法存在较大误差, 粒子群算法虽然相对于CRLB有一定偏离, 但其结果可以作为初始值代入改进的泰勒算法, 并能收敛到较为精确的结果, 称该算法为粒子群泰勒(Taylor-PSO-Mod)算法。

实验1   解析法定位解算性能研究

设集群中包含5架无人机, 按矩形布站, 辐射源位置为(5, 10, 0)km, 飞机位置为(0, 0, 3.3)km、(-3, 0, 3)km、(3, 0, 3)km、(0, 3, 3)km、(0, -3, 3)km。

为了研究存在无人机位置误差时解析法对Chan算法的改进效果, 设置站址误差范围为1~40 m, 进行1 000次蒙特卡罗实验, 以均方根误差(root mean square error, RMSE)作为定位结果的评价标准, 图 3显示了定位结果的对比。

图例中带有Mod标识的为改进方法, 下同。由仿真结果和上述理论推导, 由于在求解过程中做了两次近似, 所以解析法对Chan算法的改进并不明显, 在站址误差较低时就出现了阈值效应, 抗噪声性能并不理想。

实验2    迭代法定位解算性能研究

设置实验条件与实验1相同, 图 4显示了迭代法对Taylor算法的改进效果。

由仿真结果可知, 在站址误差较大时, 迭代法相对于泰勒算法有较好的抗噪声性能, RMSE明显降低, 但在站址误差较小时偏离了CRLB。且迭代法需要一个较为准确的初始估计值, 否则会无法收敛, 实验1的结果表明定位误差较大的解析法已不能满足要求。

实验3    Taylor-PSO-Mod算法定位解算性能研究

设目标、无人机位置与实验1中的位置相同, 图 5给出了解析法、迭代法、改进粒子群算法以及Taylor-PSO-Mod算法随站址误差变化的定位精度。

图 5中的紫色线代表本文提出的Taylor-PSO-Mod算法, 观察比较该算法与迭代法、粒子群算法的性能。在低站址误差时, Taylor-PSO-Mod算法表现最好, 在站址误差较大情况下, Taylor-PSO-Mod算法的定位精度高于粒子群算法, 与迭代法有相近的性能, 但它解决了迭代法的初值问题。其在较宽的无人机位置误差范围内, 定位结果与RMSE均最小。

图 6中经纬度误差指的是目标的经纬度, 从曲线可以看出, 不论是经度误差还是纬度误差, 解析法的表现均较差, 粒子群算法的纬度误差较大, 而Taylor-PSO-Mod算法在经纬度误差上的表现均良好, 这就使得接下来的工作有了一定的仿真基础。

为了进一步探究当无人机与目标相距不同距离时各种算法的定位性能变化, 在以上初始位置的基础上, 设置站址误差为5 m, 时差测量误差为20 ns, 令无人机以v=(40, 70, 0)m/s的速度向目标飞行, 则随飞行时间变化位置估计的RMSE曲线如图 7所示。

根据仿真结果, 在时差测量误差和站址误差确定时, 随着无人机飞行时间增大，Taylor-PSO-Mod算法对目标定位的RMSE始终保持最小, 说明了所提算法的稳定性。且无人机与辐射源相距越近, 精度越高。

进一步地, 为观察密集布站情况下所提方法的定位性能和稳定性, 参考文献[31]中的布站方式, 新增一架无人机, 无人机集群分别位于(0.3, 0.1, 0.15)km, (0.4, 0.15, 0.1)km, (0.3, 0.5, 0.2)km, (0.35, 0.2, 0.1)km, (-0.1, -0.1, -0.1)km, (-0.2, -0.3, -0.2)km。辐射源位置为(285, 325, 275)m。观察密集布站条件下所提算法的性能。该实验对曲线作了如下处理, 因为所设站址误差较小, 所以对站址误差取了对数米, 得到更为平滑的曲线，以更好地观察实验现象, 通过实验获得的仿真结果如图 8所示。

根据图 8的仿真结果, 当设置较小的站址误差且密集布站的条件下, 所提的Taylor-PSO-Mod算法定位RMSE仍然最小, 基本符合CRLB, 说明在多种仿真条件下所提算法性能具备稳定性。值得注意的是, 图 8中的粒子群算法在站址误差较大时定位误差小于CRLB, 这是由于粒子群算法是一种有偏估计, 而CRLB是无偏估计的最小方差, 在某些布站条件下, 当目标位置估计偏差较大时, 可能出现定位的RMSE小于CRLB的情况。

本文主要对基于无人机集群的时差测量单目标定位方法进行研究, 建立了单目标定位模型。作为高机动平台, 无人机无法得到自身精确位置, 需要在经典算法的基础上将站址误差考虑进来。本文重点研究存在站址误差情况下提高单目标定位精度的方法。首先，对Chan算法和Taylor算法进行改进, 得到解析法和迭代法, 在存在站址误差的情况下算法性能有不同程度的提高；然后，提出一种Taylor-PSO-Mod解算方法, 仿真实验证明了这种方法相对于解析法和迭代法的优点, 其在多种仿真场景下均能达到较高的定位精度, 定位性能稳定, 同时还解决了迭代法的初值问题。从图 5各个算法的表现来看, 站址误差较大时，其RMSE误差比性能较好的改进的粒子群算法减小了300 m左右。