1.   引言

随着现代战争环境的日益复杂，无源定位系统由于其自身不辐射电磁信号^[1]，战场生存能力强而备受关注。无源定位技术的主要原理是通过获取不同类型的观测量来估计目标状态，常用的观测量有到达角度(Angle Of Arrival, AOA)、到达时间(Time Of Arrival, TOA)、到达时间差(Time Difference Of Arrival, TDOA)、到达频率差(Frequency Difference Of Arrival, FDOA)、以及上述观测量的相互结合^[2]。其中，基于AOA的定位算法较简单，但由于目前能实现的角度测量误差较大，故该类算法对远距离目标定位时，定位精度低^[3]，不能满足战场需求。基于TOA定位需要满足目标与观测站时间同步，对系统硬件要求高^[4]，战场环境往往无法满足其要求。而使用TDOA进行目标定位，可消除TOA测量中引起定位误差的时钟偏差，从而解决了时间同步问题，系统构成简单。但只利用TDOA进行定位的算法仅能获取到目标的位置信息，而无法获取运动目标的速度信息。为实现对动目标的无源定位，需要进一步引入包含接收站与目标之间相对速度信息的FDOA^[5,6]。联合TDOA和FDOA定位可在大幅度提高目标定位精度的同时，获取到目标的速度信息，这对于后续的军事战略部署、敌情监控具有重大意义^[7]。

在实战应用中，接收站自定位精度往往不够高，文献[3]中已证明，对于未考虑接收站位置和速度误差的算法，即使在存在很小站址误差的情况下，算法的定位精度也会大幅度下降，故对动目标定位时需考虑站址误差。现有的考虑站址误差的动目标定位算法主要有两大类：迭代类和解析类。文献[8]中提出了一种迭代类算法—约束最小二乘法(Constrained Total Least Squares, CTLS)，解决了定位方程的非线性问题，但所有迭代类算法都需要拟定适当的初值，不当的初值拟定会导致算法定位性能的下降。为克服这个问题，Ho在文献[9]中提出了一种解析类算法—两步加权最小二乘(Two-Stage Weighted Least Squares, TSWLS)算法，该算法无需拟定初值，更适用于对非合作目标定位的场景。但在该算法中，当目标接近参考站的任一坐标轴时，会产生缺秩问题，从而导致算法在特定点处定位误差显著增大。为了避免这个问题，文献[10]提出了一种改进的TSWLS算法(Improved Two-Stage Weighted Least Squares, ITSWLS)，显著改善了算法的稳健性。尽管上述基于TSWLS的算法具有无需拟定初值以及计算量较小的优点，但由于其第2步存在平方、开方等非线性运算，可能会导致算法产生很大的估计误差，在中等噪声条件下，算法的定位精度仍不尽人意。为了同时避免TSWLS算法中的非线性运算和矩阵缺秩问题，文献[11]提出一种修正定位误差(Localization Error Refinement, LER)的算法，改进了TSWLS算法的步骤2，提高了算法的定位精度和稳健性，但是该算法只可用来定位静目标。之后，刘洋等人^[12]将其扩展应用到动目标定位上来，但其中的LER算法没有考虑站址误差，仍存在严重的精度损失，不适用于真实无源定位场景。综上所述，现有动目标无源定位算法仍存在3个缺陷：(1)在计算过程中出现非线性运算；(2)未考虑站址误差；(3)定位稳健性弱。

针对以上3个缺陷，本文在考虑站址误差的条件下，提出了一种动目标稳健无源定位的改进算法。算法在与TSWLS算法步骤1相同的基础上，提出了全新的步骤2。新的步骤2基于误差校正的思想，首先对步骤1的定位误差进行估计，然后用步骤1中的初估值减去定位误差估计值，从而得到更精确的目标定位结果。本文算法可有效避免TSWLS算法中的非线性运算，具有更高的定位精度，且不存在矩阵缺秩问题，因此具有更强的稳健性。数字仿真结果表明，相比于现有算法^[8-12]，本文所提算法具有更优的抗噪性和稳健性。

2.   定位模型

本文考虑三维定位场景，场景中存在M个移动接收站，利用测量值TDOA和FDOA对一个位置${{{u}}^{\rm{o}}} = {[ {{x^{\rm{o}}}}\ \ {{y^{\rm{o}}}}\ \ {{z^{\rm{o}}}} ]^{\rm{T}}}$，速度${{\dot{ u}}^{\rm{o}}} = {[ {{{\dot x}^{\rm{o}}}}\ \ {{{\dot y}^{\rm{o}}}}\ \ {{{\dot z}^{\rm{o}}}} ]^{\rm{T}}}$的动目标进行侦察与定位。本文用${( \cdot )^{\rm{o}}}$表示真实值，$( \cdot )$表示含噪声的测量值。接收站的真实位置和速度分别用${{s}}_i^{\rm{o}} = {[ {x_i^{\rm{o}}}\ \ {y_i^{\rm{o}}}\ \ {z_i^{\rm{o}}} ]^{\rm{T}}}$和${\dot{ s}}_i^{\rm{o}} = {[ {\dot x_i^{\rm{o}}}\ \ {\dot y_i^{\rm{o}}}\ \ {\dot z_i^{\rm{o}}} ]^{\rm{T}}}$表示，其中$i = 1, 2, ··· ,M$。不失一般性地假设第1个接收站为参考站^[8-12]，则目标到第$i$个接收站与到参考站之间的距离差及距离差变化率可表示为

其中，$r_i^{\rm{o}}$, $\dot r_i^{\rm{o}}$分别为目标到第$i$个接收站的距离及距离变化率,$i = 2,3, ··· ,M$，可通过TDOA, FDOA计算得到

考虑到实际情况中存在测量误差$\Delta {r_{i1}}$及$\Delta {\dot r_{i1}}$，真实测量值${r_{i1}}$, ${\dot r_{i1}}$分别可表示为$r_{i1}^{\rm{o}} + \Delta {r_{i1}}$, $\dot r_{i1}^{\rm{o}} + \Delta {\dot r_{i1}}$。为方便表示，记

设测量向量$\alpha \!=\! {[ {{{{r}}^{\rm{T}} }}\ \ {{{{\dot{ r}}}^{\rm{T}}}} ]^{\rm{T}}} \!=\! { \alpha ^{\rm{o}}} + \Delta \alpha$，其中，${ \alpha ^{\rm{o}}} \!= {[ {{{{r}}^{{\rm{o}}{\rm{T}} }}}\ \ {{{{\dot{ r}}}^{{\rm{o}}{\rm{T}} }}} ]^{\rm{T}}}$，假设测量误差向量$\Delta \alpha = {[ {\Delta {{{r}}^{\rm{T}}}}\ \ {\Delta {{{\dot{ r}}}^{\rm{T}} }} ]^{\rm{T}}}$服从高斯零均值分布，协方差${{{Q}}_{ \alpha} } = E[\Delta {{\alpha }}\Delta {{{\alpha}} ^{\rm{T}}}]$^[9]。

另外，由于站址误差的存在，实际情况中只能获取到含有误差的接收站位置${{{s}}_i} = {{s}}_i^{\rm{o}} + \Delta {{{s}}_i}$及速度${{\dot{ s}}_i} = {\dot{ s}}_i^{\rm{o}} + \Delta {{\dot{ s}}_i}$, $i = 1,2, ··· ,M$。用向量表示为

设站址向量${\beta } \!=\! {[ {{{{s}}^{\rm{T}}}}\ \ {{{{\dot{ s}}}^{\rm{T}}}} ]^{\rm{T}}} \!=\! {{\beta }^{\rm{o}}} + \Delta {\beta }$，其中，${{\beta }^{\rm{o}}}\! = {[ {{{{s}}^{{\rm{o}}{\rm{T}} }}}\ \ {{{{\dot{ s}}}^{{\rm{o}}{\rm{T}} }}} ]^{\rm{T}}}$，假设站址误差向量$\Delta {\beta } = {[ {\Delta {{{s}}^{\rm{T}}}}\ \ {\Delta {{{\dot{ s}}}^{\rm{T}}}} ]^{\rm{T}}}$服从高斯零均值分布，且与$\Delta \alpha$相互独立，协方差${{{Q}}_{{\beta }}} = E[\Delta {\beta }\Delta {{\beta }^{\rm{T}} }]$^[9]。本文将利用测量值TDOA, FDOA以及含有误差的站址向量${\beta }$实现对目标位置及速度的估计。

3.   定位算法

本文所提算法共有两步，其中，步骤2为创新点，是算法性能提升的关键。简便起见，本文仅简述步骤1流程并对其中必要的参量进行说明，而详细介绍步骤2。

3.1   初步定位

重构式(1)得$r_{i1}^{\rm{o}} + r_1^{\rm{o}} = r_i^{\rm{o}}$，将其两边平方，分别用含噪声值的${r_{i1}} - \Delta {r_{i1}}$和${{{s}}_i} - \Delta {{{s}}_i}$替换$r_{i1}^{\rm{o}}$和${{s}}_i^{\rm{o}}$，忽略2阶误差项得到一组TDOA等式

对式(5)关于时间求导得到一组FDOA等式

其中，${\varepsilon _{i,t}}$和${\varepsilon _{i,f}}$分别表示每个TDOA, FDOA等式的误差值。设${ \theta} _1^{\rm{o}} = {[ {{{{u}}^{{\rm{o}}{\rm{T}} }}}\ \ {r_1^{\rm{o}}}\ \ {{{{\dot{ u}}}^{{\rm{o}}{\rm{T}} }}}\ \ {\dot r_1^{\rm{o}}} ]^{\rm{T}} }$，将所有$i = 2, 3, ··· ,M$的TDOA, FDOA等式用矩阵表示为

其中，${{{\varepsilon}} _1} = {[\begin{array}{*{20}{c}} {{{\varepsilon}} _t^{\rm{T}} }&{{{\varepsilon }}_f^{\rm{T}} } \end{array}]^{\rm{T}} }$, ${{{\varepsilon}} _t} = {[ {{\varepsilon _{2,t}}}\ \ {{\varepsilon _{3,t}}}\ \ ··· \ \ {{\varepsilon _{M,t}}} ]^{\rm{T}}}$, ${{{\varepsilon}} _f} = {[ {{\varepsilon _{2,f}}}\ \ {{\varepsilon _{3,f}}}\ \ ···\ \ {{\varepsilon _{M,f}}} ]^{\rm{T}}}$, ${{{h}}_1} = {[\begin{array}{*{20}{c}} {{{h}}_t^{\rm{T}} }&{{{h}}_f^{\rm{T}} } \end{array}]^{\rm{T}}}$, ${{{G}}_1} =$${[ {{{G}}_t^{\rm{T}} }\ \ {{{G}}_f^{\rm{T}} } ]^{\rm{T}}}$

其中，${{O}}$表示全零矩阵

未知向量${ \theta} _1^{\rm{o}}$的最小二乘解为

加权矩阵${{{W}}_1}$可以表示为

定义步骤1的估计误差

由式(7)可知${{{\varepsilon}} _1} = {{{B}}_1}\Delta {{\alpha}} + {{{D}}_1}\Delta {\beta }$，其中，$\Delta {{\alpha }}$, $\Delta {\beta }$服从零均值高斯分布，因此在噪声较小时，$\Delta {{{\theta}} _1}$均值近似为零，也就是说${\hat {{\theta}} _1}$是渐进无偏的，其协方差矩阵为

3.2   误差校正

TSWLS算法中存在平方、开方等非线性运算，会进一步放大步骤1中的定位误差，影响算法的最终定位精度。为避免非线性运算带来的精度损失，本文提出全新步骤2，通过对定位误差进行估计并用初步定位结果减去估计误差，得到更准确的定位结果。

从步骤1中可以得到

将式(17)中的$r_1^{\rm{o}}$及$\dot r_1^{\rm{o}}$在${\hat{ u}} = {\hat{{ \theta}} _1}(1:3)$, ${\hat {\dot { u}}} = {\hat {{\theta}} _1}(5:7)$, ${{{s}}_1}$和${{\dot{ s}}_1}$处进行1阶泰勒展开

其中，

另外，根据文献[11,13]，可得另外一组等式

注意到，式(18)和式(20)都是关于$\Delta {{u}}$和$\Delta {\dot{ u}}$的线性方程，设${{\theta }}_2^{\rm{o}} = {[ {\Delta {{{u}}^{\rm{T}}}}\ \ {\Delta {{{\dot{ u}}}^{\rm{T}} }} ]^{\rm{T}}}$，由式(18)和式(20)可得

其中，

其中，${{{I}}_{i \times i}}$表示维度为$i \times i$的单位矩阵。利用加权最小二乘可得步骤1定位误差${{\theta }}_2^{\rm{o}}$的估计值

其中，

${\rm{cov}} ({{{\theta}} _1})$可由式(16)获得。${\hat{{ \theta }}_2}$的估计误差可表示为

其协方差为

用步骤1的初估值减去${\hat {{\theta }}_2}$得到目标的最终定位结果

需要注意的是，加权矩阵${{{W}}_1}$是关于目标位置的函数，所以在步骤1中不可直接使用。因此，先假设${{{W}}_1}$为等维度的单位阵，带入式(13)得到一个初始估计值${\hat {{\theta}} _1}$，然后用该初估值计算得到新的${{{W}}_1}$，再通过式(13)更新${\hat {{\theta }}_1}$，本文共更新两次。${{{W}}_2}$中的目标位置和速度可用步骤1的估计结果代替。

4.   性能分析

4.1   CRLB分析

本节通过量化所提算法的估计误差，并将其与克拉美罗下界(Cramer-Rao Lower Bound, CRLB)相比较，证明了本文算法在噪声较小的情况下可达CRLB。

令${{{\theta}} ^{\rm{o}}} = {[ {{{{u}}^{\rm{o}{\rm{T}} }}}\ \ {{{{\dot{ u}}}^{\rm{o}{\rm{T}} }}} ]^{\rm{T}}}$，将式(28)中的$[ {{{\hat{{ \theta}} }_1}{{(1:3)}^{\rm{T}} }} \ \ {{{\hat {{\theta}} }_1}{{(5:7)}^{\rm{T}}}} ]^{\rm{T}}$表示为${{{\theta }}^{\rm{o}}} + {[ {\Delta {{{u}}^{\rm{T}} }}\ \ {\Delta {{{\dot{ u}}}^{\rm{T}} }} ]^{\rm{T}}}$，等号两边同时减去${{{\theta}} ^{\rm{o}}}$，可得

本文在3.1节中已证明$\Delta {{{\theta}} _1}$均值近似为零，结合式(21)、式(26)、式(29)可得$\Delta {{{\theta}} _3}$均值近似为零，所提算法近似为无偏估计，其协方差可表示为

在一般无源定位场景下，可以认为目标到各接收站的距离近似相等^[9,14]，即$r_1^{\rm{o}} \simeq r_2^{\rm{o}} \simeq ··· \simeq r_M^{\rm{o}} = {r^{\rm{o}}}$。从文中式(7)—式(9)可以看出${{{\varepsilon}} _1}$正比于${r^o}$，从式(15)可以看出${{{\varepsilon}} _1}$还正比于$\Delta {{{\theta}} _1}$，即$\Delta {{{\theta}} _1}$与${r^{\rm{o}}}$成正比，故站址误差$\Delta {{\beta}}$通常远小于定位估计误差$\Delta {{{\theta}} _1}$，式(21)中${{{\varepsilon}} _2}$的大小主要由第1项决定，可表示为${{{\varepsilon}} _2} \simeq {{{B}}_2}\Delta {{{\theta}} _1}$，则有

且根据矩阵求逆引理，有

将式(32)代入式(16)，再将所得结果代入式(31)有

最后将式(33)代入式(30)可得

其中，${{{G}}_3} = {{B}}_1^{ - 1}{{{G}}_1}{{B}}_2^{ - 1}{{{G}}_2}$, ${{{G}}_4} = {{B}}_1^{ - 1}{{{D}}_1}$。

文献[9]中已给出站址误差条件下的定位CRLB

其中，${{X}} = {({{{\partial {{\alpha }^{\rm{o}}}} / {\partial {\theta} }}^{\rm{o}}})^{\rm{T}}}{{Q}}_{ \alpha} ^{ - 1}({{{\partial {{\alpha }^{\rm{o}}}} / {\partial {{\theta}} }}^{\rm{o}}})$, ${{Y}} = ({\partial {{\alpha }^{\rm{o}}}} / {\partial {{\theta }^{\rm{o}}}})^{\rm{T}} {{Q}}_{ \alpha} ^{ - 1}({{{\partial {{ \alpha }^{\rm{o}}}} / {\partial {\beta }}}^{\rm{o}}})$, ${{Z}} = {({{{\partial {{\alpha }^{\rm{o}}}} / {\partial {\beta }}}^{\rm{o}}})^{\rm{T}}}{{Q}}_{ \alpha} ^{ - 1}({\partial {{\alpha }^{\rm{o}}}} / {\partial {\beta }}^{\rm{o}}) + {{Q}}_{{\beta }}^{ - 1}$。显然式(34)和式(35)矩阵形式相同，通过比较可证得${{{G}}_3} = {{\partial {{\alpha }^{\rm{o}}}} / {\partial {{\theta }^{\rm{o}}}}}$, ${{{G}}_4} = {{{\partial {{\alpha }^{\rm{o}}}} / {\partial {\beta }}}^{\rm{o}}}$。故在噪声较小的情况下，所提算法的估计精度可达CRLB。(证明过程见附录A)

4.2   稳健性分析

出于对定位连贯性的要求，在对空间目标定位的过程中，不可出现误差过大的特殊位置点^[15]，即空间目标的地理位置应与算法的定位性能无关。

TSWLS算法中矩阵${{{B}}_2}$为

其中包含有与目标位置相关的元素，当目标位置接近参考站的任一坐标轴时，$({{u}} - {{{s}}_1})$中某个元素接近为零，${{{B}}_2}$变为缺秩矩阵，在式(25)中对${{{B}}_2}$求逆时就会导致定位误差较大，最终影响算法定位精度。而根据本文3.2节可知，所提算法步骤2中矩阵${{{B}}_2}$为

由式(37)可看出，本文${{{B}}_2}$内所有元素与目标位置无关，有效地避免了矩阵缺秩问题，使得本文算法具有更高的稳健性。

4.3   定位精度分析

TSWLS算法的步骤2利用辅助变量与目标位置信息间的关系构建方程。算法首先对辅助变量$r_1^{\rm{o}}$和$\dot r_1^{\rm{o}}$进行平方和乘积运算

并将其用以估计$({{{u}}^{\rm{o}}} - {{{s}}_1}) \odot ({{{u}}^{\rm{o}}} - {{{s}}_1})$和$({{\dot{ u}}^{\rm{o}}} - {{\dot{ s}}_1}) \odot ({{{u}}^{\rm{o}}} - {{{s}}_1})$，其中，$\odot$代表矩阵的Schur积(对应元素相乘)。在该过程中，TSWLS步骤2的方程误差可以表示为

式(39)将方程误差分为了两个部分：1阶误差项和2阶误差项。TSWLS算法在平方及乘积运算后直接忽略了2阶误差项，接着对步骤2的结果开方得到最终定位结果${{\breve{u} }}$和${{\breve{\dot u} }}$。算法中，平方及开方等非线性运算都会使定位误差增大，且在噪声较大的情况下，2阶误差会显著增大以至于在计算过程中不能被忽略，将其直接忽略将会造成严重的精度损失。

所提算法的步骤2不涉及非线性运算，虽然式(18)也只保留了1阶项，但其高阶相都与$\left\| {{{u}} - {{{s}}_1}} \right\|$成反比，例如${\hat r_1}$的2阶项

当目标距接收站较远时，2阶相非常小以至可以被忽略。故相比于TSWLS算法，所提算法具有更高的定位精度。

相比于文献[12]，所提算法将接收站的位置及速度误差考虑在内，更贴近现实情况，消除了由于无法获得接收站的准确位置和速度而产生的精度下降问题。虽然所提算法计算量大约为文献[12]的5倍，但考虑到所提算法定位精度的提升，一定计算量的牺牲是值得的。(计算量计算过程见附录B)

5.   仿真结果与分析

本节通过蒙特卡洛仿真实验将本文算法与ITSWLS^[10], TSWLS^[9], LER^[12]及CRLB进行比较，验证了本文所提算法的定位性能。

仿真中，本节将以L次独立的蒙特卡洛实验得到的均方根误差(Root Mean Squares Error, RMSE)来定量评估算法的定位性能。TDOA和FDOA的观测误差协方差矩阵设为${{{Q}}_{ \alpha} } \!=\! \sigma _t^2{\rm{diag}} \left\{ {{{{J}}_{M - 1}}} \ \ {0.1{{{J}}_{M - 1}}} \right\}$，接收站站址误差协方差矩阵设为${{{Q}}_{ \beta} } = \sigma _s^2{\rm{diag}} \left\{ {{{{J}}_{3M}}} \ \ {0.1{{{J}}_{3M}}} \right\}$。其中，${\sigma _t}$为TDOA测量误差，${\sigma _s}$为接收站位置误差，${{{J}}_i}$为$i$维方阵，其对角线元素为1，其余元素均为0.5。

5.1   仿真1

仿真1的目的是量化算法精度，以文献[9]中的CRLB为基准，对比所提算法与TSWLS, ITSWLS, LER算法的定位精度。与文献[9]相同，本文选取6个移动接收站，其位置与速度如表1所示

表  1  接收站位置(m)及速度(m/s)

Table  1.  Position(m) and velocity(m/s) of receivers

接收站	${x_i}$	${y_i}$	${z_i}$	${\dot x_i}$	${\dot y_i}$	${\dot z_i}$
1	300	100	150	30	–20	20
2	400	150	100	–30	10	20
3	300	500	200	10	–20	10
4	350	200	150	10	20	30
5	–100	–100	–100	–20	10	10
6	200	–300	–200	20	–10	10

下载: 导出CSV 

| 显示表格

对一个真实位置为${[\begin{array}{*{20}{c}} {2000}&{2500}&{3000} \end{array}]^{\rm{T}}}\;{\rm{m}}$，速度为${[ {-20}\ \ {15}\ \ {40} ]^{\rm{T}}}\;{\rm{m}} /{\rm s}$的动目标进行定位。分别设TDOA测量误差${\sigma _t} = {10^{ - 2}}\;{\rm{m}}$，接收站位置误差${\sigma _s} = {10^{ - 1}}\;{\rm{m}}$，并以$0.05\;{\rm{m}}$为步进长度不断增大${\sigma _s}$至${10^0}\;{\rm{m}}$，经过$L = 500$次蒙特卡洛实验后的仿真结果如图1所示。

图  1  不同算法性能比较

Figure  1.  Performance comparison among different methods

下载: 全尺寸图片 幻灯片

图1给出了当站址误差${\sigma _s}$从${10^{ - 1}}\;{\rm{m}}$变化到${10^0}\;{\rm{m}}$时，不同算法的定位RMSE和误差的变化曲线。从图中可以看出，当${\sigma _s}$不断增大时，所提算法的定位RMSE能更好地贴合CRLB，且拥有更小的定位误差，尽管在${\sigma _s} = 1\;{\rm{m}}$时仍存在轻微的门槛效应，但相比于其他算法在${\sigma _s} \approx 0.4\;{\rm{m}}$就已经开始偏离CRLB，且在${\sigma _s} = 1\;{\rm{m}}$时误差较大，所提算法的抗噪性仍明显优于其他对比算法，这也验证了4.2节中的理论分析。

5.2   仿真2

仿真2的目的是量化算法的稳健性，仍然以文献[9]中的CRLB为基准，对比所提算法与ITSWLS, LER, TSWLS算法的稳健性。为使仿真结果更直观，选取6个速度为$0\;{\rm{m}} /{\rm{s}}$，真实位置分别为${[ 0\ \ 0\ \ { - 150} ]^{\rm{T}}}\;{\rm{m}}$, ${[\begin{array}{*{20}{c}} 0&{300}&0 \end{array}]^{\rm{T}}}\;{\rm{m}}$, ${[ { - 300}\ \ 0\ \ 0 ]^{\rm{T}}}\;{\rm{m}}$, ${[\begin{array}{*{20}{c}} 0&{ - 300}&0 \end{array}]^{\rm{T}}}\;{\rm{m}}$, ${[ {300}\ \ 0\ \ 0 ]^{\rm{T}}}\;{\rm{m}}$, ${[ {150}\ \ 0\ \ 0 ]^{\rm{T}}}\;{\rm{m}}$的固定接收站，用以定位一个以$20\;{\rm{m}} /{\rm{s}}$的速度在$3000\;{\rm{m}}$高空做半径为$2000\;{\rm{m}}$圆周运动的动目标，目标的位置及速度表示为

其中，方位角$\varphi \in [0,2\pi ]$。

设测量误差${\sigma _t} \!=\! {10^{ - 2}}\;{\rm{m}}$，站址误差${\sigma _s} \!=\! {10^{ - 1}}\;{\rm{m}}$，目标方位角从0～$2\pi$变化^[14]，经过$L = 100$次蒙特卡洛实验后的仿真结果如图2所示。

图  2  不同算法稳健性比较

Figure  2.  Robustness comparison among different methods

下载: 全尺寸图片 幻灯片

图2对比展示了所提算法与ITSWLS, LER, TSWLS算法的稳健性。可以看出，ITSWLS, LER算法及本文算法在方位角变化时定位精度都始终能达到CRLB，表现稳定。而TSWLS算法由于存在矩阵缺秩问题，即使在噪声较小的情况下，每当目标方位角接近$k\pi /2$, $k = 0,1, ··· ,4$时，算法的RMSE都会显著增大。例如，当$\varphi = \pi /2$时，向量${({{{u}}^{\rm{o}}} - {{s}}_1^{\rm{o}})^{\rm{T}}} = {(0,1700,3000)^{\rm{T}} }{\rm{m}}$，此时式(36)中${{{B}}_2}$缺秩，对其求逆时会导致算法误差增大，定位RMSE激增，在图2中呈锯齿状。而本文所提算法则不存在矩阵缺秩问题，算法的定位性能与目标位置无关，具有更强的稳健性。

6.   结论

在存在站址误差的条件下，本文提出了一种基于误差校正的高精度稳健定位算法。所提算法中，全新的步骤2通过对步骤1中目标位置和速度估计值的误差$\Delta {{u}}$和$\Delta {\dot{ u}}$进行估计，并用步骤1的初步估计值减去误差估计值，最终获得精确的目标位置和速度。避免了传统基于TSWLS算法中的非线性运算及矩阵缺秩问题。理论分析和仿真实验都表明，相比于现有算法，所提算法具有更好的抗噪性，更小的定位误差及更强的稳健性。

附录A

设${{{C}}_1}(i,:) = [ { - {{a}}_1^{\rm{T}} }\ \ \ {{{O}}_{3(i - 1) \times 1}^{\rm{T}} }\ \ \ {{{a}}_{i + 1}^{\rm{T}} }\ \ \ {{{O}}_{3(M - i - 1) \times 1}^{\rm{T}} } ]$, ${{{C}}_2}(i,:) \!=\! [ { - {{b}}_1^{\rm{T}} }\ \ \ {{{O}}_{3(i - 1) \times 1}^{\rm{T}} }\ \ \ {{{b}}_{i + 1}^{\rm{T}} }\ \ \ {{{O}}_{3(M - i - 1) \times 1}^{\rm{T}} } ]$, ${{{D}}_1}(i,:) = [\begin{array}{*{20}{c}} { - {{c}}_1^{\rm{T}} }&{{{O}}_{3(i - 1) \times 1}^{\rm{T}} }&{{{c}}_{i + 1}^{\rm{T}} }&{{{O}}_{3(M - i - 1) \times 1}^{\rm{T}} } \end{array}]$, ${{{D}}_2}(i,:) = [\begin{array}{*{20}{c}} { - {{d}}_1^{\rm{T}} }&{{{O}}_{3(i - 1) \times 1}^{\rm{T}} }&{{{d}}_{i + 1}^{\rm{T}} }&{{{O}}_{3(M - i - 1) \times 1}^{\rm{T}} } \end{array}]$。其中，${{{a}}_i} = {{({{{u}}^{\rm{o}}} - {{s}}_i^{\rm{o}})} / {r_i^{\rm{o}}}}$, ${{{b}}_i} = {{({{{\dot{ u}}}^{\rm{o}}} - {\dot{ s}}_i^{\rm{o}})} / {r_i^{\rm{o}}}} - {{\dot r_i^{\rm{o}}({{{u}}^{\rm{o}}} - {{s}}_i^{\rm{o}})} / {r_i^{{\rm{o}}2}}}$, ${{{c}}_i} = {{({\hat{ u}} - {{{s}}_i})} / {{r_i}}}$, ${{{d}}_i} = {{({\hat{\dot { u}} }- {{{\dot{ s}}}_i})} / {{r_i}}} - {{{{\dot r}_i}({\hat{ u}} - {{{s}}_i})} / {r_i^2}}$, $i = 2,3, \cdots, M$

由文献[9]，噪声较小时${{\Delta {{u}}} / {r_i^{\rm{o}}}} \approx 0$, ${{\Delta {{{\dot{ u}}}_i}} / {r_i^{\rm{o}}}} \approx 0$, ${{\Delta {{{s}}_i}} / {r_i^{\rm{o}}}} \approx 0$, ${{\Delta {{{\dot{ s}}}_i}} / {r_i^{\rm{o}}}} \approx 0$

同理可得${{{b}}_i} = {{{d}}_i}$，故${{{G}}_4} = {{{\partial {{\alpha }^{\rm{o}}}} / {\partial {\beta }}}^{\rm{o}}}$

噪声较小时${{\Delta {r_{i1}}} / {r_i^{\rm{o}}}} \approx 0$，有

且已证得${{{c}}_i} = {{{a}}_i}$

同理，当${{\Delta {{u}}} / {r_i^{{\rm{o}}2}}} \approx 0$，${{\Delta {\dot{ u}}} / {r_i^{\rm{o}}}} \approx 0$，${{\Delta {{{s}}_i}} / {r_i^{{\rm{o}}2}}} \approx 0$，${{\Delta {{{\dot{ s}}}_i}} / {r_i^{\rm{o}}}} \approx 0$时

故${{{G}}_3} = {{{\partial {{\alpha }^{\rm{o}}}} / {\partial {\theta }}}^{\rm{o}}}$。得证${\rm{cov}} {({\hat {{\theta}} _3})^{ - 1}} \approx {\rm{CRLB}} {({{{\theta}} ^{\rm{o}}})^{ - 1}}$。

附录B

本文主要以算法涉及到的实数乘法次数为标准^[2]，对比分析了所提算法与文献[9]及文献[12]的计算量。此外，由于迭代类算法需要多次迭代才能获得较好的估计性能，计算量往往远大于解析类，故在此不再具体分析迭代类算法的计算量。

所提算法步骤1的计算量

(1) 计算${{{W}}_1}$：包括计算(a)${{{B}}_1}{{{Q}}_{ \alpha} }{{B}}_1^{\rm{T}}$; (b)${{{D}}_1}{{{Q}}_{ \beta} }{{D}}_1^{\rm{T}}$; (c)求逆运算。以上每个步骤分别需要进行$16{(M - 1)^3}$; $72(M - 1){M^2} + 24{(M - 1)^2}M$以及$8{(M - 1)^3}$次实数乘法。

(2) 计算${{\hat{ \varphi }}_1}$：包括计算(a)${{G}}_1^{\rm{T}} {{{W}}_1}{{{G}}_1}$; (b)${({{G}}_1^{\rm{T}} {{{W}}_1}{{{G}}_1})^{ - 1}}$; (c)${({{G}}_1^{\rm{T}} {{{W}}_1}{{{G}}_1})^{ - 1}}{{G}}_1^{\rm{T}} {{{W}}_1}{h_1}$。一共需要$64{(M - 1)^2} + 272(M - 1) + 512$次实数乘法。

所提算法步骤2的计算量

(1) 计算${{{W}}_2}$：包括计算(a)${{{B}}_2}{\rm{cov}}({{{\varphi }}_1}){{B}}_2^{\rm{T}}$; (b)${{{D}}_2}{{{Q}}_{ \beta} }{{D}}_2^{\rm{T}}$; (c)${{{B}}_2}{({{G}}_1^{\rm{T}}{{{W}}_1}{{{G}}_1})^{ - 1}}{{G}}_1^{\rm{T}}{{{W}}_1}{{{D}}_1}{{{Q}}_{ \beta} }{{D}}_2^{\rm{T}}$; (d)${{{D}}_2}{{{Q}}^{\rm{T}}_{ \beta} } \cdot {{D}}_1^{\rm{T}}{{{W}}^{\rm{T}}_1} {{{G}}_1} {({{G}}_1^{\rm{T}}{{{W}}_1}{{{G}}_1})^{ - {\rm{T}} }}{ B}_2^{\rm{T}}$; (e)求逆运算。一共需要$1184{M^2} + 832M + 3712$次实数乘法。

(2) 计算${{\hat{ \varphi }}_2}$：包括计算(a)${{G}}_2^{\rm{T}} {{{W}}_2}{{{G}}_2}$; (b)${({{G}}_2^{\rm{T}} {{{W}}_2}{{{G}}_2})^{ - 1}}$; (c)${({{G}}_2^{\rm{T}} {{{W}}_2}{{{G}}_2})^{ - 1}}{{G}}_2^{\rm{T}} {{{W}}_2}{{{h}}_2}$。一共需要1608次实数乘法。

需要注意的是，由于所提算法步骤1中需要先初始化加权矩阵${{{W}}_1}$，并更新两次得到估计结果，故步骤1中(a)需要计算两次(b)需要计算三次，因此所提算法共需要$(240{M^3} + 992{M^2} + 1312M + 5992)$次实数乘法。文献[9]中TSWLS算法在步骤2后仍需利用公式$\breve{ u} = {\rm{diag}} [{{{\varphi }}_1}(1:3) - {{ s}_1}]\sqrt {{{{\varphi }}_2}(1:3)} + {{ s}_1}$得到最终定位结果，故相比于TSWLS算法，所提算法计算量稍小。而相比于文献[12]算法，由于所提算法多考虑了接收站位置和速度的误差，在计算${{{W}}_1}$和${{{W}}_2}$时比文献[12]共多出$192{M^3} + 944{M^2} + 880M + 2176$次实数乘法。