近年来，无源定位算法由于其本身不辐射电磁信号，具有电磁隐蔽性^[1]等优点受到了国内外学者的广泛关注，它已经被大量应用在无线传感器网络、雷达、麦克风阵列以及声呐系统^[2-9]等领域中。无源定位技术可以归纳为依靠传感器网络获取不同类型的观测量进行求解，这些观测量主要包括到达时间(Time of Arrival, TOA)、到达时间差(Time Difference of Arrival, TDOA)、到达角(Angle of Arrival, AOA)、以及上述参数的相互结合。当运动目标或者传感器与目标之间存在相对运动关系时，还可以利用到达频率差(Fequency Difference of Arrival, FDOA)观测量对目标的位置以及速度进行估计。

本文主要针对运动目标无源定位算法进行研究。由于其相关的最大似然问题存在高度的非线性以及非凸性，因此寻找全局最优解变得十分困难。为此，国内外学者提出了大量的算法：线性化的最大似然估计(Maximum Likelihood Estimation, MLE)^[10-11]是解决定位问题十分有效的算法，但是这种算法的性能以及收敛性依赖于一个好的初始位置值，一旦该值选取较差其性能也会变差。为了克服这个问题，Ho和Xu提出了著名的两步加权最小二乘(Two-Stage Weighted Least Squares, TSWLS)算法^[12]，TSWLS算法第1步引入额外变量并且建立了伪线性方程，随后给出了加权最小二乘(Weighted Least Squares, WLS)解，第2步利用额外变量与目标位置的关系建立了新的等式方程以提高定位精度。虽然TSWLS算法实时性较高，但是在定位精度方面还有待进一步提高；半定松弛(SeminDefine Relaxation, SDR)算法^[13-14]首先将定位问题描述为具有二次约束条件的优化问题，随后利用合理的近似以及适当的松弛条件将其转化为了半定规划(SeminDefine Programming, SDP)问题，最后利用优化工具箱给出了具有高精度的估计值。但是这个问题的计算复杂度较高，需要在实时性以及精度两者之间进行较好的折中。文献[15]给出了一种基于定位误差修正的算法，第1步与传统的TSWLS相同，第2步对估计误差做了修正并给出了最终的定位估计值。文献[16]首先将定位问题描述为具有二次约束条件的二次规划(Quadratically Constrained Quadratic Programming, QCQP)问题，随后利用WLS的解，将二次约束转化为了线性约束，即将QCQP问题转化为了线性约束二次规划(Linearly Constrained Quadratic Programming, LCQP)问题；最后利用广义逆矩阵的性质对LCQP问题求解，并且形成一种迭代算法，这种算法具有闭式解的优点，同时在高斯噪声模型下，其定位性能优于现存的技术。

除上述提到的非凸性以及非线性问题，传感器参数存在误差时也会对无源算法的性能产生较大的影响，文献[17]详细阐述了不考虑传感器参数误差时产生的影响并且将TSWLS算法扩展到了存在传感器参数误差时的场景中；随后文献[18]中的TSWLS在其基础上进一步提高了算法精度，并且该算法适用非相交多源的定位解算。文献[13]中的SDR算法同样给出了一种考虑了传感器参数误差的情况，但是SDR算法计算复杂度远高于TSWLS算法。而近几年，文献[19-20]分别提出了存在传感器位置误差情况下改进的TSWLS算法以及多维标度(MultiDimensional Scaling, MDS)算法；但是，这些算法是基于TDOA观测量对目标定位解算，并不能应用到对运动目标参数估计的场景之中。因此针对存在传感器参数误差的情况下，研究一种运动目标无源定位的闭式算法十分有必要。

本文主要针对TDOA/FDOA运动目标定位场景，在考虑了传感器存在参数误差的情况下，提出了一种改进的两步加权最小二乘算法；该算法分为两步，其中第1步与文献[18]中经典的TSWLS算法第1步相同，首先引入了辅助向量，并且建立了伪线性方程，给出了辅助向量的WLS解；第2步中，与文献[18]有所区别，利用辅助向量中额外变量与目标参数之间的关系，构造了新的等式方程，并且利用WLS给出最终解。理论分析表明该算法在测量噪声较小时下可以实现克拉美罗下界(Cramér-Rao Lower Bound, CRLB)^[21]，并且该算法仍然是闭式解，计算复杂度与TSWLS算法相当，远小于SDR算法。此外，经过适当的维度调整，该算法同样适用于对多非相交源的定位场景。仿真实验验证了所提算法的有效性。

1 定位场景

考虑在三维空间中使用M个传感器对一个运动目标的位置u^o=[x^o, y^o, z^o]^T以及速度$\dot{\boldsymbol{u}}^{\circ}=$$\left[\dot{x}_{i}^{\mathrm{o}}, \dot{y}_{i}^{\mathrm{o}}, \dot{z}_{i}^{\mathrm{o}}\right]^{\mathrm{T}}$进行估计。第i个传感器的真实位置以及速度分别表示为$\boldsymbol{s}_{i}^{\circ}=\left[x_{i}^{0}, y_{i}^{\circ}, z_{i}^{\circ}\right]^{\top}$和$\dot{\boldsymbol{S}}_{i}^{\mathrm{o}}=$$\left[\dot{x}_{i}^{\circ}, \dot{y}_{i}^{o}, \dot{z}_{i}^{\circ}\right]^{\mathrm{T}}$，其中，i=1, 2, …, M。然而，只能获取得到不准确传感器的位置和速度，这些参数可以分解为真实参数值加一些随机噪声：

$ \begin{array}{*{20}{l}} {{\mathit{\boldsymbol{s}}_i} = \mathit{\boldsymbol{s}}_i^{\rm{o}} + \Delta {\mathit{\boldsymbol{s}}_i} = {{\left[ {{x_i},{y_i},{z_i}} \right]}^{\rm{T}}}}\\ {{{\mathit{\boldsymbol{\dot s}}}_i} = \mathit{\boldsymbol{\dot s}}_i^{\rm{o}} + \Delta {{\mathit{\boldsymbol{\dot s}}}_i} = {{\left[ {{{\dot x}_i},{{\dot y}_i},{{\dot z}_i}} \right]}^{\rm{T}}}} \end{array} $	（1）

式中：s_i和${\mathit{\boldsymbol{\dot s}}_\mathit{i}}$分别为第i个传感器可利用的位置以及速度；Δs_i和$\Delta {\boldsymbol{\dot s}_i}$为零均值高斯噪声向量。

将所有可利用的传感器位置和速度值使用向量形式表示为

$ \begin{array}{l} \mathit{\boldsymbol{s}} = {\left[ {\mathit{\boldsymbol{s}}_1^{\rm{T}},\mathit{\boldsymbol{s}}_2^{\rm{T}}, \cdots ,\mathit{\boldsymbol{s}}_M^{\rm{T}}} \right]^{\rm{T}}} = {\mathit{\boldsymbol{s}}^{\rm{o}}} + \Delta \mathit{\boldsymbol{s}}\\ \mathit{\boldsymbol{\dot s}} = {\left[ {\mathit{\boldsymbol{\dot s}}_1^{\rm{T}},\mathit{\boldsymbol{\dot s}}_2^{\rm{T}}, \cdots ,\mathit{\boldsymbol{\dot s}}_M^{\rm{T}}} \right]^{\rm{T}}} = {{\mathit{\boldsymbol{\dot s}}}^{\rm{o}}} + \Delta \mathit{\boldsymbol{\dot s}} \end{array} $

（2）

式中：${\mathit{\boldsymbol{s}}^\circ } = {\left[ {\mathit{\boldsymbol{s}}_1^{{\rm{oT}}},\mathit{\boldsymbol{s}}_2^{{\rm{oT}}}, \cdots ,\mathit{\boldsymbol{s}}_M^{{\rm{oT}}}} \right]^{\rm{T}}}{\rm{且}} \cdot {\mathit{\boldsymbol{\dot s}}^{\rm{o}}} = \left[ {\mathit{\boldsymbol{s}}_1^{{\rm{oT}}},\mathit{\boldsymbol{s}}_2^{{\rm{oT}}}, \cdots } \right.,$${\left. {\mathit{\boldsymbol{\dot s}}_M^{{\rm{oT}}}} \right]^{\rm{T}}};\Delta \mathit{\boldsymbol{s}} = {\left[ {\Delta \mathit{\boldsymbol{s}}_1^{\rm{T}},\Delta \mathit{\boldsymbol{s}}_2^{\rm{T}}, \cdots ,\Delta \mathit{\boldsymbol{s}}_M^{\rm{T}}} \right]^{\rm{T}}}{\rm{且}}\Delta \mathit{\boldsymbol{\dot s}} = \left[ {\Delta \mathit{\boldsymbol{\dot s}}_{\rm{1}}^{\rm{T}}} \right.,$${\left. {\Delta \mathit{\boldsymbol{\dot s}}_{\rm{2}}^{\rm{T}}, \cdots ,\Delta \mathit{\boldsymbol{\dot s}}_M^{\rm{T}}} \right]^{\rm{T}}}$。

记$\mathit{\boldsymbol{\beta }} = {\left[ {{\mathit{\boldsymbol{s}}^{\rm{T}}}, {{\mathit{\boldsymbol{\dot s}}}^{\rm{T}}}} \right]^{\rm{T}}}$和${\mathit{\boldsymbol{\beta }}^\circ } = {\left[ {{\mathit{\boldsymbol{s}}^{^\circ {\rm{T}}}}, {\mathit{\boldsymbol{s}}^{{\rm{oT}}}}} \right]^{\rm{T}}}$，那么总的测量噪声向量可以表示为

$ \Delta \mathit{\boldsymbol{\beta }} = \mathit{\boldsymbol{\beta }} - {\mathit{\boldsymbol{\beta }}^{\rm{o}}} = {\left[ {\Delta {\mathit{\boldsymbol{s}}^{\rm{T}}},\Delta {{\mathit{\boldsymbol{\dot s}}}^{\rm{T}}}} \right]^{\rm{T}}} $	（3）

式中：总的测量噪声向量Δβ服从零均值高斯分布，其协方差矩阵为E(Δβ^TΔβ)=Q_β。

不失一般性地，选取第1个传感器为参考传感器。因此由M个传感器可以得到(M－1)个TDOA值，乘以信号传播速度之后，得到如下表达式：

$ r_{i1}^{\rm{o}} = r_i^{\rm{o}} - r_1^{\rm{o}}\;\;\;i = 2,3, \cdots ,M $	（4）

式中：r_i1^o为TDOA乘以信号传播速度得到的距离差(Range Difference, RD)。r_i^o为目标与第i个传感器之间的距离，其表达式为

$ r_i^{\rm{o}} = \left\| {{\mathit{\boldsymbol{u}}^{\rm{o}}} - \mathit{\boldsymbol{s}}_i^{\rm{o}}} \right\|\;\;\;\;i = 2,3, \cdots ,M $	（5）

对式(4)两边关于时间求导有

$ \dot r_{i1}^{\rm{o}} = \dot r_i^{\rm{o}} - \dot r_1^{\rm{o}}\;\;\;i = 2,3, \cdots ,M $	（6）

式中：$\dot{r}_{i 1}^{\mathrm{o}}$为距离差变化率(Range Difference Rate, RDR)，并且可以通过FDOA计算得到^[12]。$\dot r_{i1}^{\rm{o}}$为距离变化率，且有

$ \dot r_i^{\rm{o}} = {\left( {{\mathit{\boldsymbol{u}}^{\rm{o}}} - \mathit{\boldsymbol{s}}_i^{\rm{o}}} \right)^{\rm{T}}}\left( {{{\mathit{\boldsymbol{\dot u}}}^{\rm{o}}} - \mathit{\boldsymbol{\dot s}}_i^{\rm{o}}} \right)/r_i^{\rm{o}} $	（7）

分别记RD以及RDR的测量向量为r和$\mathit{\boldsymbol{\dot r}}$，那么它们可以表示为

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{r}} = {\left[ {{r_{21}},{r_{31}}, \cdots ,{r_{M1}}} \right]^{\rm{T}}} = {\mathit{\boldsymbol{r}}^{\rm{o}}} + \Delta \mathit{\boldsymbol{r}}\\ \mathit{\boldsymbol{\dot r}} = {\left[ {{{\dot r}_{21}},{{\dot r}_{31}}, \cdots ,{{\dot r}_{M1}}} \right]^{\rm{T}}} = {{\mathit{\boldsymbol{\dot r}}}^{\rm{o}}} + \Delta \mathit{\boldsymbol{\dot r}} \end{array} \right. $

（8）

式中：r^o=[r₂₁^o, r₃₁^o, …, r_M1^o]^T和${\mathit{\boldsymbol{\dot r}}^{\rm{o}}} = $$ = {\left[ {\Delta {{\dot r}_{21}},\Delta {{\dot r}_{31}}, \cdots ,\Delta {{\dot r}_{M1}}} \right]^T}$为真实的RD和RDR向量；Δr=[Δr₂₁, Δr₃₁, …, Δr_M1]^T以及$\Delta \mathit{\boldsymbol{\dot r}} = \left[ {\Delta {{\dot r}_{21}},\Delta {{\dot r}_{31}}} \right.,{\left. { \cdots ,\Delta {{\dot r}_{M1}}} \right]^{\rm{T}}}$为服从零均值高斯分布的噪声向量。

同样，记α=[r^T, ${\mathit{\boldsymbol{\dot r}}^{\rm{T}}}$]^T为总的RD和RDR测量向量，那么噪声向量可以表示为

$ \Delta \mathit{\boldsymbol{\alpha }} = \mathit{\boldsymbol{\alpha }} - {\mathit{\boldsymbol{\alpha }}^{\rm{o}}} = {\left[ {\Delta {\mathit{\boldsymbol{r}}^{\rm{T}}},\Delta {{\mathit{\boldsymbol{\dot r}}}^{\rm{T}}}} \right]^{\rm{T}}} $	（9）

式中：${\mathit{\boldsymbol{\alpha }}^\circ } = {\left[ {{\mathit{\boldsymbol{r}}^{{\rm{oT}}}}, {{\mathit{\boldsymbol{\dot r}}}^{{\rm{oT}}}}} \right]^{\rm{T}}}$为真实的RD和RDR向量值；Δα服从零均值高斯分布并且与Δβ相互独立，其协方差矩阵为E(Δα^TΔα)=Q_α。

2 算法

本节具体给出所提出算法的推导过程，由于所提出的算法第1步与文献[18]的第1步相同，因此第1步仅给出简单的总结；在新的第2步中，详细给出了推导过程。

2.1 算法第1步

首先将额外变量r₁^o在s₁处进行一阶泰勒级数展开有

$ r_1^{\rm{o}} = \left\| {{\mathit{\boldsymbol{u}}^{\rm{o}}} - \mathit{\boldsymbol{s}}_1^{\rm{o}}} \right\| \approx d_1^{\rm{o}} + \mathit{\boldsymbol{\rho }}_{{\mathit{\boldsymbol{u}}^{\rm{o}}},{\mathit{\boldsymbol{s}}_1}}^{\rm{T}}\Delta {\mathit{\boldsymbol{s}}_1} $	（10）

式中：$d_1^0 = \left\| {{\mathit{\boldsymbol{u}}^\circ } - {{\boldsymbol{s}}_1}} \right\|且, {\mathit{\boldsymbol{\rho }}_{{\mathit{\boldsymbol{u}}^0}, {{\boldsymbol{s}}_1}}} = \left( {{\mathit{\boldsymbol{u}}^\circ } - {s_1}} \right)/\left\| {{{\boldsymbol{u}}^\circ } - {{\boldsymbol{s}}_1}} \right\|$。

对式(4)移项并且两边进行平方，将RD测量值r_i1、可利用的传感器位置s_i以及式(10)代入，同时仅保留线性误差项可以得到一组TDOA等式为

$ \begin{array}{l} {\varepsilon _i} = 2r_i^{\rm{o}}\Delta {n_{i1}} + 2{\left( {{\mathit{\boldsymbol{u}}^{\rm{o}}} - {\mathit{\boldsymbol{s}}_i}} \right)^{\rm{T}}}\Delta {\mathit{\boldsymbol{s}}_i} - 2\mathit{\boldsymbol{\eta }}_{1i}^{\rm{T}}\Delta {\mathit{\boldsymbol{s}}_1} = \\ \;\;\;r_{i1}^2 + {\left\| {{\mathit{\boldsymbol{s}}_1}} \right\|^2} - {\left\| {{\mathit{\boldsymbol{s}}_i}} \right\|^2} + 2{\left( {{\mathit{\boldsymbol{s}}_i} - {\mathit{\boldsymbol{s}}_1}} \right)^{\rm{T}}}{\mathit{\boldsymbol{u}}^{\rm{o}}} + 2{r_{i1}}d_1^o\\ \;\;\;\left( {i = 2,3, \cdots ,M} \right) \end{array} $

（11）

式中：${\mathit{\boldsymbol{\eta }}_{1i}} = {\mathit{\boldsymbol{u}}^\circ } - {\mathit{\boldsymbol{s}}_1} + {r_{i1}}{\mathit{\boldsymbol{\rho }}_{{u^0}, {s_1}}}$。

对式(11)两边关于时间求导可以得到一组FDOA方程为

$ \begin{array}{l} {{\dot \varepsilon }_i} = 2\dot r_i^{\rm{o}}\Delta {n_{i1}} + 2r_i^{\rm{o}}\Delta {{\dot n}_{i1}} + 2{\left( {{{\mathit{\boldsymbol{\dot u}}}^{\rm{o}}} - {{\mathit{\boldsymbol{\dot s}}}_i}} \right)^{\rm{T}}}\Delta {\mathit{\boldsymbol{s}}_i} - \\ \;\;\;\;\;\;2\mathit{\boldsymbol{\eta }}_{2i}^{\rm{T}}\Delta {\mathit{\boldsymbol{s}}_1} + 2{\left( {{\mathit{\boldsymbol{u}}^o} - {\mathit{\boldsymbol{s}}_i}} \right)^{\rm{T}}}\Delta {{\mathit{\boldsymbol{\dot s}}}_i} - 2\mathit{\boldsymbol{\eta }}_{1i}^{\rm{T}}\Delta {{\mathit{\boldsymbol{\dot s}}}_1} = \\ \;\;\;\;\;\;2{\left( {{{\mathit{\boldsymbol{\dot s}}}_i} - {{\mathit{\boldsymbol{\dot s}}}_1}} \right)^{\rm{T}}}{\mathit{\boldsymbol{u}}^{\rm{o}}} + 2{\left( {{\mathit{\boldsymbol{s}}_i} - {\mathit{\boldsymbol{s}}_1}} \right)^{\rm{T}}}{{\mathit{\boldsymbol{\dot u}}}^{\rm{o}}} + \\ \;\;\;\;\;\;2{r_{i1}}d_1^{\rm{o}} + 2{r_{i1}}\dot d_1^{\rm{o}} + 2\left( {{{\dot r}_{i1}}{r_{i1}} + \mathit{\boldsymbol{\dot s}}_1^{\rm{T}}{\mathit{\boldsymbol{s}}_1} - \mathit{\boldsymbol{\dot s}}_i^{\rm{T}}{\mathit{\boldsymbol{s}}_i}} \right)\\ \;\;\;\;\;\;\left( {i = 2,3, \cdots ,M} \right) \end{array} $

（12）

式中：$\dot{d}_{1}^{\circ}=\left(\boldsymbol{u}^{\circ}-\boldsymbol{s}_{1}\right)^{\mathrm{T}}\left(\boldsymbol{u}^{o}-\dot{\boldsymbol{s}}_{1}\right) /\left\|\boldsymbol{u}^{o}-\boldsymbol{s}_{1}\right\|$，并且

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{\eta }}_{2i}} = \left( {{{\mathit{\boldsymbol{\dot u}}}^{\rm{o}}} - {{\mathit{\boldsymbol{\dot s}}}_1} + \gamma {r_{i1}} + {{\dot r}_{i1}}{\mathit{\boldsymbol{\rho }}_{{\mathit{\boldsymbol{u}}^{\rm{o}}},{\mathit{\boldsymbol{s}}_1}}}} \right)\\ \mathit{\boldsymbol{\gamma }} = \left( {{{\mathit{\boldsymbol{\dot u}}}^{\rm{o}}} - {{\mathit{\boldsymbol{\dot s}}}_1} + \dot d_1^{\rm{o}}{\mathit{\boldsymbol{\rho }}_{{\mathit{\boldsymbol{u}}^{\rm{o}}},{\mathit{\boldsymbol{s}}_1}}}} \right)/d_1^{\rm{o}} \end{array} \right. $

（13）

定义$\boldsymbol{\varphi}_{1}^{\circ}=\left[\boldsymbol{u}^{\mathrm{oT}}, d_{1}^{\mathrm{o}}, \dot{\boldsymbol{u}}^{\mathrm{oT}}, \dot{d}_{1}^{\mathrm{o}}\right]^{\mathrm{T}}$为未知的辅助向量，其中d₁^o和$\dot d_1^o$为新的额外变量。将所有的TDOA和FDOA等式用矩阵形式表示，那么可以得到

$ {\mathit{\boldsymbol{\varepsilon }}_1} = {\mathit{\boldsymbol{B}}_1}\Delta \alpha + {\mathit{\boldsymbol{D}}_1}\Delta \beta = {\mathit{\boldsymbol{h}}_1} - {\mathit{\boldsymbol{G}}_1}\mathit{\boldsymbol{\varphi }}_1^o $	（14）

式中：h₁=[h_t^T, h_f^T]^T为一个2(M－1)维的列向量；G₁=[G_t^T, G_f^T]^T为一个2(M－1)×8的矩阵；ε₁为服从零均值高斯分布的等式误差向量；h_t和h_f的每个元素分别表示为

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{h}}_{\rm{t}}}\left( {i - 1} \right) = r_{i1}^2 - \mathit{\boldsymbol{s}}_i^{\rm{T}}{\mathit{\boldsymbol{s}}_i} + \mathit{\boldsymbol{s}}_1^{\rm{T}}{\mathit{\boldsymbol{s}}_1}\\ {\mathit{\boldsymbol{h}}_{\rm{f}}}\left( {i - 1} \right) = 2\left( {{r_{i1}}{{\dot r}_{i1}} - \mathit{\boldsymbol{\dot s}}_i^{\rm{T}}{\mathit{\boldsymbol{s}}_i} + \mathit{\boldsymbol{\dot s}}_1^{\rm{T}}{\mathit{\boldsymbol{s}}_1}} \right) \end{array} \right. $

（15）

同样，矩阵G_t和G_f的每一行元素分别为

$ \begin{array}{l} {\mathit{\boldsymbol{G}}_{\rm{t}}}\left( {i - 1,:} \right) = - 2\left[ {\begin{array}{*{20}{l}} {{{\left( {{\mathit{\boldsymbol{s}}_i} - {\mathit{\boldsymbol{s}}_1}} \right)}^{\rm{T}}}}&{{r_{i1}}}&{{\bf{0}}_{3 \times 1}^{\rm{T}}}&0 \end{array}} \right]\\ {\mathit{\boldsymbol{G}}_{\rm{f}}}\left( {i - 1,:} \right) = - 2\left[ {{{\left( {{{\mathit{\boldsymbol{\dot s}}}_i} - {{\mathit{\boldsymbol{\dot s}}}_1}} \right)}^{\rm{T}}}{{\dot r}_{i1}}\quad {{\left( {{\mathit{\boldsymbol{s}}_i} - {\mathit{\boldsymbol{s}}_1}} \right)}^{\rm{T}}}{r_{i1}}} \right] \end{array} $

（16）

另外，$\boldsymbol{B}_{1}=\left[\begin{array}{ll} \boldsymbol{B}_{1 t} & \boldsymbol{O} \\ \boldsymbol{B}_{1 f} & \boldsymbol{B}_{1 t} \end{array}\right]$且B_1t=2diag(r₂^o, r₃^o, ..., r_M^o)以及B_1f=2diag $\left(\dot{r}_{2}^{o}, \dot{r}_{3}^{o}, \cdots, \dot{r}_{M}^{o}\right)$；D₁ $=\left[\begin{array}{ll} \boldsymbol{D}_{1 t} & \boldsymbol{O} \\ \boldsymbol{D}_{1 \mathrm{f}} & \boldsymbol{D}_{1 t} \end{array}\right]$，且D_1f和D_1t，的每一行元素分别表示为

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{D}}_{1{\rm{t}}}}\left( {i - 1,:} \right) = 2\left[ { - \mathit{\boldsymbol{\eta }}_{1i}^{\rm{T}},{\bf{0}}_{3(i - 2) \times 1}^{\rm{T}},{{\left( {{\mathit{\boldsymbol{u}}^o} - {\mathit{\boldsymbol{s}}_i}} \right)}^{\rm{T}}},} \right.\\ \;\;\;\;\;\;\left. {{\bf{0}}_{3\left( {M - i} \right) \times 1}^{\rm{T}}} \right]\\ {{\mathit{\boldsymbol{\dot D}}}_{1{\rm{f}}}}\left( {i - 1,:} \right) = 2\left[ { - \mathit{\boldsymbol{\eta }}_{2{\rm{i}}}^{\rm{T}},{\bf{0}}_{3(i - 2) \times 1}^{\rm{T}},{{\left( {{{\mathit{\boldsymbol{\dot u}}}^o} - {{\mathit{\boldsymbol{\dot s}}}_i}} \right)}^{\rm{T}}},} \right.\\ \;\;\;\;\;\;\left. {{\bf{0}}_{3\left( {M - i} \right) \times 1}^{\rm{T}}} \right] \end{array} \right. $

（17）

矩阵方程(14)的加权最小二乘解为

$ {{\mathit{\boldsymbol{\hat \varphi }}}_1} = {\left( {\mathit{\boldsymbol{G}}_1^{\rm{T}}{\mathit{\boldsymbol{W}}_1}{\mathit{\boldsymbol{G}}_1}} \right)^{ - 1}}\mathit{\boldsymbol{G}}_1^{\rm{T}}{\mathit{\boldsymbol{W}}_1}{\mathit{\boldsymbol{h}}_1} $	（18）

式中：u=φ₁(1:3)和$\dot{\boldsymbol{u}}=\boldsymbol{\varphi}_{1}$(5:7)为第1步得到的位置和速度估计值；${\mathit{\boldsymbol{W}}_1} = \left( {{\mathit{\boldsymbol{B}}_1}{\mathit{\boldsymbol{Q}}_a}\mathit{\boldsymbol{B}}_1^{\rm{T}} + } \right.{\left. {{\mathit{\boldsymbol{D}}_1}{\mathit{\boldsymbol{Q}}_\beta }\mathit{\boldsymbol{D}}_1^ \top } \right)^{ - 1}}$为加权矩阵。另外，令d₁=φ₁(4)且${\dot d_1} = {\mathit{\boldsymbol{\varphi }}_1}(8)$分别为d₁^o和$\dot d_1^o$的估计值。

记$\Delta {\mathit{\boldsymbol{\varphi }}_1} = {\left[ {\Delta {\mathit{\boldsymbol{u}}^{\rm{T}}},\Delta {d_1},\Delta {{\mathit{\boldsymbol{\dot u}}}^{\rm{T}}},\Delta {{\dot d}_1}} \right]^{\rm{T}}} = {\mathit{\boldsymbol{\hat \varphi }}_1} - \mathit{\boldsymbol{\varphi }}_1^{\rm{o}}$为估计误差，则利用式(18)估计误差可以表示为

$ \Delta {\mathit{\boldsymbol{\varphi }}_1} = {\left( {\mathit{\boldsymbol{G}}_1^{\rm{T}}{\mathit{\boldsymbol{W}}_1}{\mathit{\boldsymbol{G}}_1}} \right)^{ - 1}}\mathit{\boldsymbol{G}}_1^{\rm{T}}{\mathit{\boldsymbol{W}}_1}{\mathit{\boldsymbol{\varepsilon }}_1} $	（19）

由于式(19)中ε₁均值为0，因此第1步的估计值在测量噪声较小时近似为无偏估计；经过简单的矩阵运算，式(18)估计值的协方差矩阵近似为${\mathop{\rm cov}} {\mathit{\boldsymbol{\widehat \varphi }}_1} = {\left( {\mathit{\boldsymbol{G}}_1^{\rm{T}}{\mathit{\boldsymbol{W}}_1}{\mathit{\boldsymbol{G}}_1}} \right)^{ - 1}}$。

2.2 算法第2步

第1步中引入的辅助向量是假设目标参数与额外变量d₁^o以及$\dot d_1^o$之间相互独立，但是实际上，它们之间存在非线性关系为

$ \left\{ {\begin{array}{*{20}{l}} {d_1^{{\rm{o}}2} = {{\left( {{\mathit{\boldsymbol{u}}^{\rm{o}}} - {\mathit{\boldsymbol{s}}_1}} \right)}^{\rm{T}}}\left( {{\mathit{\boldsymbol{u}}^{\rm{o}}} - {\mathit{\boldsymbol{s}}_1}} \right)}\\ {\dot d_1^{\rm{o}}d_1^{\rm{o}} = {{\left( {{\mathit{\boldsymbol{u}}^{\rm{o}}} - {\mathit{\boldsymbol{s}}_1}} \right)}^{\rm{T}}}\left( {{{\mathit{\boldsymbol{\dot u}}}^{\rm{o}}} - {{\mathit{\boldsymbol{\dot s}}}_1}} \right)} \end{array}} \right. $

（20）

文献[18]的第2步同样考虑了这种非线性关系，并且建立了等式方程，而本文第2步建立了一个不同于文献[18]的等式方程；另外，文献[18]中的TSWLS第2步得到结果之后还需要进行一些其他操作(如平方根运算等)才能得到最终结果，相比之下，算法第2步更简洁明了，可以直接得到最终的WLS估计结果；第2步的具体过程如下。

根据第1步的估计误差Δφ₁，真实向量值φ₁^o中的各个元素可以表示为：

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{u}}^{\rm{o}}} = \mathit{\boldsymbol{u}} - \Delta \mathit{\boldsymbol{u}},{\mathit{\boldsymbol{u}}^{\rm{o}}} = \mathit{\boldsymbol{u}} - \Delta \mathit{\boldsymbol{u}}}\\ {d_1^o = {d_1} - \Delta {d_1},\dot d_1^o = {{\dot d}_1} - \Delta {{\dot d}_1}} \end{array}} \right. $	（21）

将式(21)代入式(20)，并且只保留线性误差项，则有如下等式：

$ \begin{array}{*{20}{c}} {\left( {{\mathit{\boldsymbol{u}}^{\rm{T}}}\Delta \mathit{\boldsymbol{u}} - {d_1}\Delta {d_1}} \right) = 0.5\left( {{\mathit{\boldsymbol{u}}^{\rm{T}}}\mathit{\boldsymbol{u}} + \mathit{\boldsymbol{s}}_1^{\rm{T}}{\mathit{\boldsymbol{s}}_1} - d_1^2} \right) - }\\ {\mathit{\boldsymbol{s}}_1^{\rm{T}}{\mathit{\boldsymbol{u}}^{\rm{o}}}\left( {{{\mathit{\boldsymbol{\dot u}}}^{\rm{T}}}\Delta \mathit{\boldsymbol{u}} + {\mathit{\boldsymbol{u}}^{\rm{T}}}\Delta \mathit{\boldsymbol{\dot u}} - {{\dot d}_1}\Delta {d_1} - {d_1}\Delta {{\dot d}_1}} \right) = }\\ {\left( {{{\mathit{\boldsymbol{\dot u}}}^{\rm{T}}}\mathit{\boldsymbol{u}} + \mathit{\boldsymbol{\dot s}}_1^{\rm{T}}{\mathit{\boldsymbol{s}}_1} - {{\dot d}_1}{d_1}} \right) - \mathit{\boldsymbol{\dot s}}_1^{\rm{T}}{\mathit{\boldsymbol{u}}^{\rm{o}}} - \mathit{\boldsymbol{s}}_1^{\rm{T}}{{\mathit{\boldsymbol{\dot u}}}^{\rm{o}}}} \end{array} $

（22）

另外，第1步得到的位置和速度估计误差可以进一步改写为

$ \Delta \mathit{\boldsymbol{u}} = \mathit{\boldsymbol{u}} - {\mathit{\boldsymbol{u}}^{\rm{o}}},\Delta \mathit{\boldsymbol{\dot u}} = \mathit{\boldsymbol{\dot u}} - {{\mathit{\boldsymbol{\dot u}}}^{\rm{o}}} $	（23）

联合式(22)和式(23)，将它们用矩阵形式表示可以建立新的矩阵方程为

$ {\mathit{\boldsymbol{B}}_2}\Delta {\mathit{\boldsymbol{\varphi }}_1} = {\mathit{\boldsymbol{h}}_2} - {\mathit{\boldsymbol{G}}_2}\mathit{\boldsymbol{\varphi }}_2^{\rm{o}} $	（24）

式中：$\mathit{\boldsymbol{\varphi }}_2^\circ = {\left[ {{\mathit{\boldsymbol{u}}^{{\rm{oT}}}}, {{\mathit{\boldsymbol{\dot u}}}^{{\rm{oT}}}}} \right]^{\rm{T}}}$为第2步中待求的位置和速度值；并且

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{h}}_2} = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\dot u}}}\\ {0.5\left( {{{\mathit{\boldsymbol{\dot u}}}^{\rm{T}}}\mathit{\boldsymbol{\dot u}} + \mathit{\boldsymbol{s}}_1^{\rm{T}}{\mathit{\boldsymbol{s}}_1} - \hat d_1^2} \right)}\\ {\mathit{\boldsymbol{\dot u}}}\\ {\left( {{{\mathit{\boldsymbol{\hat {\dot u}}}}^{\rm{T}}}\mathit{\boldsymbol{\dot u}} + \mathit{\boldsymbol{\dot s}}_1^{\rm{T}}{\mathit{\boldsymbol{s}}_1} - \hat {\dot d}_1^2{{\hat d}_1}} \right)} \end{array}} \right],{\mathit{\boldsymbol{G}}_2} = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{I}}&\mathit{\boldsymbol{O}}\\ {\mathit{\boldsymbol{s}}_1^{\rm{T}}}&\mathit{\boldsymbol{O}}\\ \mathit{\boldsymbol{O}}&\mathit{\boldsymbol{I}}\\ {\mathit{\boldsymbol{\dot s}}_1^{\rm{T}}}&{\mathit{\boldsymbol{s}}_1^{\rm{T}}} \end{array}} \right]\\ {\mathit{\boldsymbol{B}}_2} = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{I}}&{{{\bf{0}}_{3 \times 1}}}&\mathit{\boldsymbol{O}}&{{{\bf{0}}_{3 \times 1}}}\\ {{\mathit{\boldsymbol{u}}^{\rm{T}}}}&{ - {d_1}}&{{\bf{0}}_{3 \times 1}^{\rm{T}}}&0\\ \mathit{\boldsymbol{O}}&{{{\bf{0}}_{3 \times 1}}}&\mathit{\boldsymbol{I}}&{{{\bf{0}}_{3 \times 1}}}\\ {{{\mathit{\boldsymbol{\dot u}}}^{\rm{T}}}}&{ - {{\dot d}_1}}&{{\mathit{\boldsymbol{u}}^{\rm{T}}}}&{ - {d_1}} \end{array}} \right] \end{array} \right. $

（25）

式(24)的WLS解为

$ {{\mathit{\boldsymbol{\hat \varphi }}}_2} = {\left( {\mathit{\boldsymbol{G}}_2^{\rm{T}}{\mathit{\boldsymbol{W}}_2}{\mathit{\boldsymbol{G}}_2}} \right)^{ - 1}}\mathit{\boldsymbol{G}}_2^{\rm{T}}{\mathit{\boldsymbol{W}}_2}{\mathit{\boldsymbol{h}}_2} $	（26）

式中：加权矩阵${\mathit{\boldsymbol{W}}_2} = {\left( {{\mathit{\boldsymbol{B}}_2}{\mathop{\rm cov}} \left( {{{\mathit{\boldsymbol{\widehat \varphi }}}_1}} \right)\mathit{\boldsymbol{B}}_2^T} \right)^{ - 1}}$。

记$\Delta {\mathit{\boldsymbol{\varphi }}_2} = {\mathit{\boldsymbol{\widehat \varphi }}_2} - \mathit{\boldsymbol{\varphi }}_2^\circ $为第2步的估计误差，将式(26)以及$\mathit{\boldsymbol{\varphi }}_2^\circ = {\left( {\mathit{\boldsymbol{G}}_2^{\rm{T}}{\mathit{\boldsymbol{W}}_2}{\mathit{\boldsymbol{G}}_2}} \right)^{ - 1}}\mathit{\boldsymbol{G}}_2^{\rm{T}}{\mathit{\boldsymbol{W}}_2}{\mathit{\boldsymbol{G}}_2}\mathit{\boldsymbol{\varphi }}_2^\circ $代入估计误差中，Δφ₂的可以表示为

$ \Delta {{\mathit{\boldsymbol{\hat \varphi }}}_2} = {\left( {\mathit{\boldsymbol{G}}_2^{\rm{T}}{\mathit{\boldsymbol{W}}_2}{\mathit{\boldsymbol{G}}_2}} \right)^{ - 1}}\mathit{\boldsymbol{G}}_2^{\rm{T}}{\mathit{\boldsymbol{W}}_2}{\mathit{\boldsymbol{B}}_2}\Delta {\mathit{\boldsymbol{\varphi }}_1} $	（27）

同样，由于E(Δφ₁)=0，因此当测量噪声较小时第2步的估计近似为无偏估计，其协方差矩阵近似为

$ {\rm cov}{{\mathit{\boldsymbol{\hat \varphi }}}_2} = {\left( {\mathit{\boldsymbol{G}}_2^{\rm{T}}{\mathit{\boldsymbol{W}}_2}{\mathit{\boldsymbol{G}}_2}} \right)^{ - 1}} $	（28）

至此，本文算法通过式(26)给出了最终的目标位置和速度估计值分别为$\mathit{\boldsymbol{\dot u}} = {\mathit{\boldsymbol{\hat \varphi }}_2}(1:3),\mathit{\boldsymbol{\dot u}} = {\mathit{\boldsymbol{\hat \varphi }}_\mathit{\boldsymbol{2}}}(4:6)$。

此外，还需要注意：

1) 加权矩阵W₁和W₂的获取均需要真实位置和速度，我们首先对加权矩阵近似为W₁=Q_α^－1，并利用式(18)得到初始的参数估计值；随后利用获得的参数值更新加权矩阵W₁和W₂，得到一个更精确的估计值。仿真实验表明了这种加权矩阵的近似带来的性能损失是可以忽略不计的。

2) 若使用该算法对多非相交源的参数进行估计时，只需要对相应的矩阵进行维度调整，对此本文不在赘述，仅给出大致思路：其中，第1步与文献[18]算法的第1步相同；第2步中φ₂^o扩展为全部目标的位置和速度向量，每个目标均能构造出矩阵G₂以及B₂，将这些矩阵联合为块对角矩阵即可建立新的矩阵G₂和矩阵B₂，进而类似的矩阵方程也随之建立；最后同样利用加权最小二乘方法给出估计解。

3 性能分析

本节将对所提出算法的性能进行分析，给出本文定位场景下的CRLB，并且理论推导证明在测量噪声较小时，所提算法可以达到CRLB。

上面已经提到过，本文算法给出的目标位置以及速度估计在测量噪声较小时是无偏的，因此只需要证明协方差矩阵${\mathop{\rm cov}} {\mathit{\boldsymbol{\widehat \varphi }}_2}$可以到达CRLB即可。将第1步估计值的协方差矩阵${\mathop{\rm cov}} {\mathit{\boldsymbol{\widehat \varphi }}_1}$代入${\mathop{\rm cov}} {\mathit{\boldsymbol{\widehat \varphi }}_2}$的表达式中有

$ \begin{array}{l} {\left( {{\rm cov}{{\mathit{\boldsymbol{\hat \varphi }}}_2}} \right)^{ - 1}} = \mathit{\boldsymbol{G}}_3^{\rm{T}}\mathit{\boldsymbol{Q}}_\alpha ^{ - 1}{\mathit{\boldsymbol{G}}_3} - \\ \;\;\;\;\;\;\;\mathit{\boldsymbol{G}}_3^{\rm{T}}\mathit{\boldsymbol{Q}}_\alpha ^{ - 1}{\mathit{\boldsymbol{G}}_4}{\left( {\mathit{\boldsymbol{Q}}_\beta ^{ - 1} + \mathit{\boldsymbol{G}}_4^{\rm{T}}\mathit{\boldsymbol{Q}}_\alpha ^{ - 1}{\mathit{\boldsymbol{G}}_4}} \right)^{ - 1}}\mathit{\boldsymbol{G}}_4^{\rm{T}}\mathit{\boldsymbol{Q}}_\alpha ^{ - 1}{\mathit{\boldsymbol{G}}_3} \end{array} $

（29）

式中：$\boldsymbol{G}_{3}=\boldsymbol{B}_{1}^{-1} \boldsymbol{G}_{1} \boldsymbol{B}_{2}^{-1} \boldsymbol{G}_{2}, \boldsymbol{G}_{4}=\boldsymbol{B}_{1}^{-1} \boldsymbol{D}_{1}$。

根据文献[17], CRLB(φ₂^o)的表达式为

$ {\rm{CRL}}{{\rm{B}}^{ - 1}}\left( {\mathit{\boldsymbol{\varphi }}_2^{\rm{o}}} \right) = \mathit{\boldsymbol{X}} - \mathit{\boldsymbol{Y}}{\mathit{\boldsymbol{Z}}^{ - 1}}{\mathit{\boldsymbol{Y}}^{\rm{T}}} $	（30）

式中:

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{X}} = {\left( {\frac{{\partial {\mathit{\boldsymbol{\alpha }}^{\rm{o}}}}}{{\partial \varphi _2^{\rm{o}}}}} \right)^{\rm{T}}}\mathit{\boldsymbol{Q}}_\alpha ^{ - 1}\frac{{\partial {\mathit{\boldsymbol{\alpha }}^{\rm{o}}}}}{{\partial \mathit{\boldsymbol{\varphi }}_2^{\rm{o}}}}\\ \mathit{\boldsymbol{Y}} = {\left( {\frac{{\partial {\mathit{\boldsymbol{\alpha }}^{\rm{o}}}}}{{\partial \mathit{\boldsymbol{\varphi }}_2^{\rm{o}}}}} \right)^{\rm{T}}}\mathit{\boldsymbol{Q}}_\alpha ^{ - 1}\frac{{\partial {\mathit{\boldsymbol{\alpha }}^{\rm{o}}}}}{{\partial {\mathit{\boldsymbol{\beta }}^{\rm{o}}}}}\\ \mathit{\boldsymbol{Z}} = \mathit{\boldsymbol{Q}}_\beta ^{ - 1} + {\left( {\frac{{\partial {\mathit{\boldsymbol{\alpha }}^{\rm{o}}}}}{{\partial {\mathit{\boldsymbol{\beta }}^{\rm{o}}}}}} \right)^{\rm{T}}}\mathit{\boldsymbol{Q}}_\alpha ^{ - 1}\frac{{\partial {\mathit{\boldsymbol{\alpha }}^{\rm{o}}}}}{{\partial {\mathit{\boldsymbol{\beta }}^{\rm{o}}}}} \end{array} \right. $

（31）

附录A中给出了$\partial {\mathit{\boldsymbol{\alpha }}^\circ }/\partial \mathit{\boldsymbol{\varphi }}_2^\circ $和$\partial {\mathit{\boldsymbol{\alpha }}^{\rm{o}}}/\partial {\mathit{\boldsymbol{\beta }}^{\rm{o}}}$具体的表达式，同时证明了${\mathit{\boldsymbol{G}}_3} = \partial {\mathit{\boldsymbol{\alpha }}^\circ }/\partial \mathit{\boldsymbol{\varphi }}_2^\circ $和${\mathit{\boldsymbol{G}}_4} = - \partial {{\mathit{\boldsymbol{\alpha }}}^\circ }/\partial {{\bf{ \pmb{\mathit{ β}} }}^\circ }$在噪声水平较小时成立。显然，式(29)和式(30)具有相同的矩阵形式，因此利用附录的结果可以得到

$ {\rm cov}{{\mathit{\boldsymbol{\hat \varphi }}}_2} = {\rm{CRLB}}\left( {\mathit{\boldsymbol{\varphi }}_2^{\rm{o}}} \right) $	（32）

即本文提出的算法在测量噪声较小时能够实现CRLB。

4 仿真实验

本节给出了仿真实验进一步验证本文所提出算法的定位性能。假设TDOA与FDOA测量误差之间相互独立，RD及RDD测量误差协方差矩阵设置为Q_α=σ²blkdiag(R, 0.01R)，这里，测量噪声的方差设置为σ²=10^－4，R为对角线元素均为1，其他元素为0.5的矩阵。与文献[13, 17]相同，仿真实验选取了6个传感器，所有传感器位置和速度的真实值在表 1中给出。可利用的传感器位置和速度值由真实值添加高斯噪声得到，噪声协方差矩阵的设置与文献[13]相同，为Q_β=σ_s²blkdiag(J, 0.5J)，其中σ_s²为噪声功率且J=[1, 1, 1, 2, 2, 2, 10, 10, 10, 40, 40, 40, 20, 20, 20, 3, 3, 3]。通过进行多次蒙特卡罗实验获取目标位置和速度估计的均方根误差(Root Mean Square Error, RMSE)和偏差(Bias)来评估算法的性能，蒙特卡罗仿真次数设置为l=10 000次，均方根误差以及偏差的指标定义为

$ \begin{array}{l} {\rm RMSE}\left( {{{\mathit{\boldsymbol{\hat \varphi }}}_2}} \right) = \sqrt {\sum\limits_{i = 1}^l {{{\left\| {\mathit{\boldsymbol{\hat \varphi }}_2^i - \mathit{\boldsymbol{\varphi }}_2^{\rm{o}}} \right\|}^2}} /l} \\ {\rm{bias}}\left( {{{\mathit{\boldsymbol{\hat \varphi }}}_2}} \right) = \sum\limits_{i = 1}^l {{{\left\| {\mathit{\boldsymbol{\hat \varphi }}_2^i - \mathit{\boldsymbol{\varphi }}_2^{\rm{o}}} \right\|}^2}} /l \end{array} $

（33）

式中：$\mathit{\boldsymbol{\widehat \varphi }}_2^i$为第i次蒙特卡罗实验得到的目标参数估计值。仿真实验考虑了3种不同情形下的定位场景，场景1和2分别为不同算法对近场以及远场目标参数的估计性能，对比算法选取了没有考虑传感器误差的基于定位误差修正的算法^[15]和ICWLS算法^[16]，以及考虑了传感器误差的TSWLS算法^[17]和SDR算法^[13]。场景3中将本文算法延伸到了对多个目标的定位场景中，对比算法为文献[18]中的TSWLS算法。另外，在场景3中，估计参数的均方根误差和偏差计算方法同式(33)，唯一的区别在于场景3中得到的$\mathit{\boldsymbol{\widehat \varphi }}_2^i$是第i次蒙特卡罗实验得到的包含所有目标位置和速度的估计值。

图 1和图 2给出的是不同算法对近场目标位置以及速度估计性能的比较，目标位置和速度分别选取为u^o=[305, 500, 165]^T m和${\mathit{\boldsymbol{\dot u}}^{\rm{o}}}$=[－20, 15, 40]^Tm/s。由图 1可知，不管是位置估计还是速度估计的均方根误差，没有考虑传感器噪声的ICWLS算法以及基于定位误差修正的算法即使在测量噪声较小时也不能实现CRLB，这主要是因为这些算法所忽略的传感器位置误差在估计性能中占了主导地位。当噪声较小时，本文提出的算法和对比算法^{[13, 17]}都可以到达CRLB，这也证明了第3节理论分析的正确性。而文献[17]算法的位置估计均方根误差在噪声大于-5 dB时开始逐渐偏离CRLB；它的速度估计均方根误差则在大于-15 dB开始下降，而且出现了阈值效应；这一点，也是传统TSWLS算法的缺点，即当参考传感器位置的某一维坐标与目标的位置坐标比较接近时，该算法会有较差的定位性能，这一点在文献[15]中也曾提及。而即使在测量噪声较大时，本文提出的算法和SDR算法位置和速度估计的性能也十分接近CRLB，与TSWLS相比，显然这两种算法在测量噪声较大时，具有更强的适应能力。

如图 2所示，对于近场目标估计偏差性能而言，没有考虑传感器误差的算法^[15-16]具有较大的位置和速度估计偏差。当测量噪声较小时，TSWLS算法能够提供和SDR算法以及本文提出算法比较接近的估计偏差，但是当噪声继续增大时，TSWLS的偏差性能开始下降。当测量噪声达到5 dB时，TSWLS算法的目标位置估计偏差比本文提出的算法大约高出了6 dB，同样，目标速度估计偏差比本文提出的算法大约高出了6 dB。SDR算法与本文所提出的算法在整个噪声功率范围内对于目标参数估计的偏差比较接近，均能够实现较小的位置与速度估计偏差，这两种算法具有较好的定位性能。但是需要强调的是，SDR算法求解优化问题需要大量的计算复杂度，该算法实时性难以保证；本文算法是一种闭式解法，因此具有较高的实时性(见第5节表 2)。

图 3和图 4给出的是不同算法对远场目标位置以及速度估计性能的比较，目标位置和速度选取为u^o=[2 000, 2 500, 3 000]^T m和${\mathit{\boldsymbol{\dot u}}^{\rm{o}}}$=[－20, 15, 40]^T m/s。如图 3均方根误差性能所示，对于远场目标参数估计，没有考虑传感器误差的算法^[15-16]依旧无法给出较好的定位性能。另外，当噪声功率小于－20 dB时，TSWLS算法、SDR算法以及本文算法法均能够实现CRLB，但是噪声继续增大时，3种算法出现了不同程度的偏离CRLB，而本文算法偏离程度较小，具有更好的鲁棒性。图 4给出的是远场目标参数估计的偏差性能，与近场目标情况类似，在所有算法中，SDR算法和本文算法依旧保持了最小的目标位置和速度估计偏差。但是正如前文所述，本文所提出的算法与SDR算法相比最大的优势就是其较低的计算复杂度(见第5节)。

另外，本文算法还可以扩展到对非相交多目标定位解算的场景中；图 5和图 6给出了不同算法的估计性能。目标位置分别选取为u₁^o=[305, 500, 165]^T m和u₂^o=[600, 650, 500]^T m；速度为$\mathit{\boldsymbol{\dot u}}_1^{\rm{o}}$=[－20, 15, 40]^T m/s和$\mathit{\boldsymbol{\dot u}}_2^{\rm{o}}$=[40, －20, 15]^Tm/s，存在噪声的传感器参数设置方法与前2个场景中的设置方法相同。

如图 5所示，对于位置估计均方根误差，当噪声小于-5 dB时，2种算法均能够实现CRLB。但是当噪声继续增大时，2种算法定位性能开始下降，当噪声为0 dB时，本文算法得到的均方根误差比文献[18]的TSWLS得到的均方根误差大约降低了4 dB。对于速度估计均方根误差，当噪声大于-15 dB时，文献[18]TSWLS算法的估计性能出现了较大的下降，而本文算法在噪声大于-5 dB时才开始偏离CRLB，本文算法具有更好的鲁棒性。图 6给出了2种算法对于目标参数估计的偏差性能，对于位置估计偏差性能，文献[18]的算法与本文算法十分接近；但是对于速度估计偏差，本文算法显然比文献[18]的算法具有更小的偏差，当测量噪声等于-5 dB时，本文算法比文献[18]算法的估计偏差大约减少了7 dB。

5 计算量分析

本节主要对所提出算法的计算复杂度进行分析，为了分析简便，主要以实数乘法进行分析。

首先分析本文算法第1步的计算量。

1) 计算W₁：包括计算① B₁Q_αB₁^T，② D₁Q_βD₁^T，③求逆运算，以上每个步骤分别需要进行16(M－1)³、36(M－1)M²以及8(M－1)³次实数乘法。

2) 计算${\mathit{\boldsymbol{\widehat \varphi }}_1}$：需要分别计算① G₁^TW₁，② ${\mathit{\boldsymbol{G}}}_1^{\rm{T}}{{\mathit{\boldsymbol{W}}}_1}{\mathit{\boldsymbol{G}}}$，③ ${\left( {{\mathit{\boldsymbol{G}}}_1^{\rm{T}}{{\mathit{\boldsymbol{W}}}_1}{{\mathit{\boldsymbol{G}}}_1}} \right)^{ - 1}}$，④ ${\left( {{\mathit{\boldsymbol{G}}}_1^{\rm{T}}{{\mathit{\boldsymbol{W}}}_1}{{\mathit{\boldsymbol{G}}}_1}} \right)^{ - 1}}{\mathit{\boldsymbol{G}}}_1^{\rm{T}}{{\mathit{\boldsymbol{W}}}_1}$，⑤ ${\left( {\mathit{\boldsymbol{G}}_1^{\rm{T}}{\mathit{\boldsymbol{W}}_1}{\mathit{\boldsymbol{G}}_1}} \right)^{ - 1}}\mathit{\boldsymbol{G}}_1^{\rm{T}}{\mathit{\boldsymbol{W}}_1}{\mathit{\boldsymbol{h}}_1}$，一共需要32(M－1)²+272(M－1)+512次实数乘法。

接下来给出第2步计算量的分析过程，第2步中计算与观测站数目无关，因此能够得到准确的乘法次数。

1) 计算W₂：包括计算① B₂^－T，② B₂^－Tcov^－1$\mathit{\boldsymbol{\widehat \varphi }}$₁，③ B₂^－T(cov₁ $\mathit{\boldsymbol{\widehat \varphi }}$₁)^－1B₂，一共需要进行3×8³=1 536次实数乘法。

2) 计算$\hat{\boldsymbol{\varphi}}_{2}$：包括计算① G₂^TW₂，② G₂^TW₂G₂，③ (G₂^TW₂G₂)^－1，④ (G₂^TW₂G₂)^－1G₂^T，⑤ $\left(\boldsymbol{G}_{2}^{\mathrm{T}} \boldsymbol{W}_{2} \boldsymbol{G}_{2}\right)^{-1} \boldsymbol{G}_{2}^{\mathrm{T}} \boldsymbol{W}_{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}$一共需要进行2 760次实数乘法。

此外，由于TSWLS的第1步需要对加权矩阵初始化，并且更新加权矩阵W₁一般需要重复计算两次，则第1步中1)需要计算2次，2)需要计算3次，因此本文提出的算法一共需要(120M³－120M²+688M+5 688)次的实数乘法运算。同样，文献[17-18]中的TSWLS算法第2步均与观测站数目无关，计算量主要都体现在了第1步中，因此本文算法与它们计算量相当，这从表 2中也可以看出，表 2中给出了考虑传感器参数误差存在情况下不同算法的平均CPU运行时间，显然TSWLS算法以及本文算法计算量相当。而SDR算法需要进行大量的运行时间。本文算法具有相当高的实时性，这一点在定位解算中具有重要意义。

6 结论

1) 所提出算法在测量噪声较小时能够实现CRLB。另外，经过适当的维度调整，可以应用于多非相交源的定位场景中。

2) 所提算法计算复杂度与TSWLS算法相当，同时能提供不高于SDR算法的均方根误差和偏差估计性能。

附录A

这里给出了$\partial {\mathit{\boldsymbol{\alpha }}^\circ }/\partial \mathit{\boldsymbol{\varphi }}_2^\circ 和\partial {\mathit{\boldsymbol{\alpha }}^\circ }/\partial {\mathit{\boldsymbol{\beta }}^\circ }$的具体表达式以及测量噪声较小时$\boldsymbol{G}_{3}=\partial \boldsymbol{\alpha}^{\circ} / \partial \boldsymbol{\varphi}_{2}^{\circ}$和${\mathit{\boldsymbol{G}}_4} - \partial {\mathit{\boldsymbol{\alpha }}^\circ }/\partial {\mathit{\boldsymbol{\beta }}^\circ }$的具体证明过程。首先当测量噪声较小时有如下条件成立

条件①：$\left| {\Delta {r_{i1}}} \right| \ll \left| {r_i^o} \right|与\left| {\Delta {{\dot r}_{i1}}} \right| \ll \left| {r_i^o} \right|$

条件②：$\left| {\Delta {s_i}} \right| \ll \left| {r_i^o} \right|与\left| {\Delta {{\dot s}_i}} \right| \ll \left| {r_i^o} \right|$

$\partial \boldsymbol{\alpha}^{\circ} / \partial \boldsymbol{\varphi}_{2}^{\circ}$写为分块矩阵可以表示为

$ \frac{{\partial {\mathit{\boldsymbol{\alpha }}^{\rm{o}}}}}{{\partial \mathit{\boldsymbol{\varphi }}_2^{\rm{o}}}} = \left[ {\begin{array}{*{20}{l}} {\frac{{\partial {\mathit{\boldsymbol{r}}^{\rm{o}}}}}{{\partial {\mathit{\boldsymbol{u}}^{\rm{o}}}}}}&{\frac{{\partial {\mathit{\boldsymbol{r}}^{\rm{o}}}}}{{\partial {{\mathit{\boldsymbol{\dot u}}}^{\rm{o}}}}}}\\ {\frac{{\partial {{\mathit{\boldsymbol{\dot r}}}^{\rm{o}}}}}{{\partial {\mathit{\boldsymbol{u}}^{\rm{o}}}}}}&{\frac{{\partial {{\mathit{\boldsymbol{\dot r}}}^{\rm{o}}}}}{{\partial {{\mathit{\boldsymbol{\dot u}}}^{\rm{o}}}}}} \end{array}} \right] $

（A1）

式中：$\frac{\partial \boldsymbol{r}^{\mathrm{o}}(\boldsymbol{\theta})}{\partial \dot{\boldsymbol{u}}}=\boldsymbol{O}_{(M-1) \times 3}, \frac{\partial \boldsymbol{r}^{\mathrm{o}}}{\partial \boldsymbol{u}^{\mathrm{o}}}, \frac{\partial \dot{\boldsymbol{r}}^{\mathrm{o}}}{\partial \boldsymbol{u}^{\mathrm{o}}}$与$\frac{\partial \dot{\boldsymbol{r}}^{\mathrm{o}}}{\partial \dot{\boldsymbol{u}}^{\mathrm{o}}}$的第i－1行分别表示为

$ \frac{{\partial {\mathit{\boldsymbol{r}}^{\rm{o}}}}}{{\partial {\mathit{\boldsymbol{u}}^{\rm{o}}}}}\left( {i - 1,:} \right) = \mathit{\boldsymbol{\rho }}_{{\mathit{\boldsymbol{u}}^{\rm{o}}},\mathit{\boldsymbol{s}}_i^{\rm{o}}}^{\rm{T}} - \mathit{\boldsymbol{\rho }}_{{\mathit{\boldsymbol{u}}^{\rm{o}}},\mathit{\boldsymbol{s}}_1^{\rm{o}}}^{\rm{T}} $	（A2）

$ \frac{{\partial {{\mathit{\boldsymbol{\dot r}}}^{\rm{o}}}}}{{\partial {{\mathit{\boldsymbol{\dot u}}}^{\rm{o}}}}}\left( {i - 1,:} \right) = \mathit{\boldsymbol{\rho }}_{{\mathit{\boldsymbol{u}}^{\rm{o}}},\mathit{\boldsymbol{s}}_i^{\rm{o}}}^{\rm{T}} - \mathit{\boldsymbol{\rho }}_{{\mathit{\boldsymbol{u}}^{\rm{o}}},\mathit{\boldsymbol{s}}_1^{\rm{o}}}^{\rm{T}} $	（A3）

$ \frac{{\partial {\mathit{\boldsymbol{r}}^{\rm{o}}}\left( \mathit{\boldsymbol{\theta }} \right)}}{{\partial {{\mathit{\boldsymbol{\dot u}}}^{\rm{o}}}}}\left( {i - 1,:} \right) = \mathit{\boldsymbol{\gamma }}_{{\mathit{\boldsymbol{u}}^{\rm{o}}},\mathit{\boldsymbol{s}}_i^{\rm{o}}}^{\rm{T}} - \mathit{\boldsymbol{\gamma }}_{{\mathit{\boldsymbol{u}}^{\rm{o}}},\mathit{\boldsymbol{s}}_1^{\rm{o}}}^{\rm{T}} $

（A4）

式中：$\boldsymbol{\rho}_{u^{o} \cdot s_{i}}=\frac{\boldsymbol{u}^{o}-\boldsymbol{s}_{i}^{o}}{r_{i}^{o}}, \boldsymbol{\gamma}_{u^{o} \cdot s_{i}}=\frac{\dot{\boldsymbol{u}}^{\circ}-\dot{\boldsymbol{s}}_{i}^{o}}{r_{i}^{o}}-\frac{\boldsymbol{\rho}_{u^{o}}^{o} \cdot s_{i}^{o}}{r_{i}^{o}} r_{i}^{o}$且i=2, 3, …, M－1。

同样，$\partial \boldsymbol{\alpha}^{\circ} / \partial \boldsymbol{\beta}^{\circ}$写为分块矩阵为

$ \frac{{\partial {\mathit{\boldsymbol{\alpha }}^{\rm{o}}}}}{{\partial {\mathit{\boldsymbol{\beta }}^{\rm{o}}}}} = \left[ {\begin{array}{*{20}{l}} {\frac{{\partial {\mathit{\boldsymbol{r}}^{\rm{o}}}}}{{\partial {\mathit{\boldsymbol{s}}^{\rm{o}}}}}}&{\frac{{\partial {\mathit{\boldsymbol{r}}^{\rm{o}}}}}{{\partial {{\mathit{\boldsymbol{\dot s}}}^{\rm{o}}}}}}\\ {\frac{{\partial {{\mathit{\boldsymbol{\dot r}}}^{\rm{o}}}}}{{\partial {\mathit{\boldsymbol{s}}^{\rm{o}}}}}}&{\frac{{\partial {{\mathit{\boldsymbol{\dot r}}}^{\rm{o}}}}}{{\partial {{\mathit{\boldsymbol{\dot s}}}^{\rm{o}}}}}} \end{array}} \right] $

（A5）

式中：$\frac{{\partial {\mathit{\boldsymbol{r}}^{\rm{o}}}}}{{\partial {{\mathit{\boldsymbol{\dot s}}}^{\rm{o}}}}} = {\mathit{\boldsymbol{O}}_{(M - 1) \times 3M}},\frac{{\partial {\mathit{\boldsymbol{r}}^{\rm{o}}}}}{{\partial {\mathit{\boldsymbol{s}}^{\rm{o}}}}} \cdot |\frac{{\partial {{\mathit{\boldsymbol{\dot r}}}^{\rm{o}}}}}{{\partial {\mathit{\boldsymbol{s}}^{\rm{o}}}}}{\rm{与}}\frac{{\partial {{\mathit{\boldsymbol{\dot r}}}^{\rm{o}}}}}{{\partial {{\mathit{\boldsymbol{\dot s}}}^{\rm{o}}}}}$的第i－1分别表示为

$ \begin{array}{l} \frac{{\partial {\mathit{\boldsymbol{r}}^{\rm{o}}}}}{{\partial {\mathit{\boldsymbol{s}}^{\rm{o}}}}}\left( {i - 1,:} \right) = \frac{{\partial {\mathit{\boldsymbol{r}}^{\rm{o}}}}}{{\partial {{\mathit{\boldsymbol{\dot s}}}^{\rm{o}}}}}\left( {i - 1,:} \right)\\ \;\;\;\;\;\; = \left[ {\mathit{\boldsymbol{\rho }}_{{\mathit{\boldsymbol{u}}^{\rm{o}}},\mathit{\boldsymbol{s}}_1^{\rm{o}}}^{\rm{T}},{\bf{0}}_{3\left( {i - 2} \right) \times 1}^{\rm{T}}, - \mathit{\boldsymbol{\rho }}_{{\mathit{\boldsymbol{u}}^{\rm{o}}},\mathit{\boldsymbol{s}}_i^{\rm{o}}}^{\rm{T}},{\bf{0}}_{3\left( {M - i} \right) \times 1}^{\rm{T}}} \right] \end{array} $

（A6）

$ \frac{{\partial {{\mathit{\boldsymbol{\dot r}}}^{\rm{o}}}}}{{\partial {\mathit{\boldsymbol{s}}^{\rm{o}}}}}\left( {i - 1,:} \right) = \left[ {\mathit{\boldsymbol{\gamma }}_{{\mathit{\boldsymbol{u}}^{\rm{o}}},\mathit{\boldsymbol{s}}_1^{\rm{o}}}^{\rm{T}},{\bf{0}}_{3\left( {i - 2} \right) \times 1}^{\rm{T}}, - \mathit{\boldsymbol{\gamma }}_{{\mathit{\boldsymbol{u}}^{\rm{o}}},\mathit{\boldsymbol{s}}_i^{\rm{o}}}^{\rm{T}},{\bf{0}}_{3\left( {M - i} \right) \times 1}^{\rm{T}}} \right] $

（A7）

且i=2, 3, …, M－1。

将相关矩阵代入G₃=B₁^－1G₁B₂^－1G₂以及G₄=B₁^－1D₁，经过一些简单的矩阵运算，同时利用条件①和条件②可以得到

$ {\mathit{\boldsymbol{G}}_3} = \partial {\mathit{\boldsymbol{\alpha }}^{\rm{o}}}/\partial \mathit{\boldsymbol{\varphi }}_2^{\rm{o}},{\mathit{\boldsymbol{G}}_4} = - \partial {\mathit{\boldsymbol{\alpha }}^{\rm{o}}}/\partial {\mathit{\boldsymbol{\beta }}^{\rm{o}}} $	（A8）