一种基于加权多维标度分析的多个非相关源TDOA/FDOA协同定位方法

王鼎，尹洁昕，*，张欣光，郑娜娥

1.中国人民解放军战略支援部队信息工程大学 信息系统工程学院，郑州 450001

2.国家数字交换系统工程技术研究中心，郑州 450002

3.北京航天长征飞行器研究所，北京 100076

4.中国人民解放军战略支援部队信息工程大学 数据与目标工程学院，郑州 450001

众所周知，辐射源定位技术已广泛应用于无线通信、目标监测、航天航空、安全管理等诸多领域［1-3］，其在工业生产和国防安全中均发挥重要作用。对于静止辐射源，常使用的定位观测量为到达时间差（Time Difference of Arrival，TDOA）［4-6］。对于运动辐射源，多普勒效应会产生到达频率差（Frequency Difference of Arrival，FDOA）［7-8］，此时联合TDOA和FDOA两种观测量进行定位可获得更高精度。近些年来，国内外学者提出若干TDOA/FDOA定位方法，主要分为迭代定位方法和闭式定位方法两大类。

迭代定位方法包括泰勒（Taylor）级数迭代方法［9］、半正定松弛方法［10］、约束总体最小二乘（Constrained Total Least Squares，CTLS）估计方法［11］、约束加权最小二乘（Constrained Weighted Least Squares，CWLS）估计方法［12］、迭代约束加权最小二乘（Iterative Constrained Weighted Least Squares， ICWLS）估计方法［13］、Bi迭代方法［14］等。虽然这些方法在一定条件下能获得渐近统计最优的估计精度，但需要复杂计算量，并且大多对迭代初始值较敏感，有些还需要设置迭代步长，易导致发散和局部收敛。

闭式定位方法能直接给出辐射源位置和速度显式表达式，虽然其需要对加权矩阵进行更新，但不用选择迭代初始值和设置迭代步长，通常对加权矩阵更新3次就能满足性能要求，并不存在发散和局部收敛问题，一般具有更低复杂度。最经典闭式定位方法是两步加权最小二乘（Two Step Weighted Least Squares，TSWLS）估计方法及其改进型方法［15-17］，该类方法通过引入中间变量建立2组伪线性观测方程，2次使用线性加权最小二乘估计器获得辐射源位置和速度估计值。经典闭式定位方法虽能有效克服迭代定位方法中的缺点，但该类方法需要推导伪线性观测方程，因此其普适性略差于迭代定位方法，更重要的是其产生“门限效应”的误差阈值相对较小，在大观测误差条件下的定位精度会出现快速下降现象。

另一类改进型闭式定位方法是加权多维标度定位方法［18-25］。多维标度分析是一种将多维空间对象简化到低维空间进行处理，同时还保留对象间原始关系的数据分析方法。近年来，国内外学者提出若干基于加权多维标度分析的定位方法。具体而言，文献［18-20］提出基于到达时间（Time of Arrival，TOA）的加权多维标度定位方法，文献［21］提出基于接收信号强度（Received Signal Strength，RSS）的加权多维标度定位方法，文献［22-23］针对无线传感网节点定位问题，提出相应加权多维标度定位方法，文献［24］提出基于TDOA的加权多维标度定位方法，文献［25］针对运动辐射源提出基于TDOA/FDOA的加权多维标度定位方法。从这些文献中不难发现，加权多维标度定位方法具备2个优势。第1个优势是可以给出辐射源位置和速度显式表达式；第2个优势是相比其他定位方法，其对大观测误差具有更强鲁棒性和稳健性，这是因为该类方法能充分利用标量积矩阵的维度信息和特征结构信息，并且标量积矩阵包含与距离及其变化率相关的全部信息。

上述定位方法均假设传感器位置和速度精确已知，当传感器安装在机载或舰载平台，又或是传感器随机布设时，传感器位置和速度精确值可能无法获知，仅能得其先验观测值，其中含有先验观测误差（也称模型误差）［26］。为了抑制模型误差的影响，国内外学者提出两类定位方法。第1类方法是将模型误差与TDOA/FDOA观测误差等同看待，通过对各种估计器重新设置加权矩阵以提高对模型误差的鲁棒性。例如，文献［27］提出基于TDOA/FDOA的鲁棒总体最小二乘（Total Least Squares，TLS）定位方法，文献［28］提出基于TDOA/FDOA的鲁棒CTLS定位方法，文献［29］提出基于TDOA/FDOA的鲁棒CWLS定位方法，文献［26］提出基于TDOA/FDOA的鲁棒TSWLS定位方法，文献［30］提出基于TDOA/FDOA的鲁棒加权多维标度定位方法，文献［31］提出基于TDOA/FDOA/DDR（Differential Doppler Rate）的鲁棒TSWLS定位方法。这些方法均可有效减弱传感器位置和速度先验观测误差对定位性能的影响。第2类方法是对辐射源与传感器位置和速度进行联合估计，该类方法不仅可以给出更准确的传感器位置和速度估计值，还能比第1类方法具有更高误差阈值。例如，文献［32-33］提出基于Taylor级数迭代的辐射源与传感器位置和速度联合估计方法。

前面提到的定位方法仅考虑对单个辐射源进行定位，在实际场景中有时需要对多个辐射源进行定位，此时应考虑多源协同定位，以提高整体定位精度。当传感器位置和速度存在先验观测误差时，即便对非相关辐射源而言，通过协同定位也可有效提升每个辐射源的定位精度。事实上，协同定位本身就可以看成是抑制模型误差的有效途径。文献［34］提出基于TSWLS的多辐射源TDOA/FDOA协同定位方法，文献［35］在文献［34］的基础上提出基于TSWLS的改进型多辐射源TDOA/FDOA协同定位方法，该方法实现了对多个非相关辐射源以及各个传感器位置和速度的联合估计，其在大观测误差条件下具有更高的定位性能。

基于上述讨论，本文在传感器位置和速度存在先验观测误差情况下，研究多个非相关辐射源TDOA/FDOA协同定位方法。为了提高在大观测误差条件下定位性能，文中首次提出一种基于加权多维标度分析的多个非相关源TDOA/FDOA协同定位方法。该方法包含2个计算阶段。阶段a基于多维标度分析原理构造两组标量积矩阵，并由此形成定位关系式，用于获得多个非相关辐射源以及各个传感器位置和速度的估计值；阶段b建立阶段a估计误差的约束优化模型，并由此得到阶段a估计误差的估计值，用于更新阶段a的估计结果。此外，文中还对新方法进行理论性能分析，从数学上证明其对辐射源以及传感器位置和速度估计精度均能渐近逼近克拉美罗界（Cramér-Rao Bound，CRB）。与现有加权多维标度定位方法相比，新方法具有以下4个优势：① 实现多辐射源协同定位；② 实现辐射源与传感器位置和速度联合估计；③ 利用2个阶段确保最终定位结果的渐近统计最优性；④ 利用矩阵正交变换技术解决定位关系式中误差协方差矩阵秩亏损问题。

这里给出文中使用的数学符号： ⊗表示矩阵Kronecker积；⊙表示矩阵点积；In表示n×n阶单位矩阵，其中第m列向量记为；On×m表示n×m阶全零矩阵；1n×m表示n×m阶全1矩阵；=[On×1In]T；diag[·]表示由向量元素构成的对角矩阵；blkdiag{·}表示由矩阵或向量作为对角元素构成的块状对角矩阵；range{·}表示矩阵列空间；Λm-n表示满足等式vec()=vec(Λm-nAnm)的0-1置换矩阵（其中Anm为任意n×m阶矩阵）；Π[·]表示矩阵列空间的正交投影矩阵，Π⊥[·] 表示矩阵列补空间的正交投影矩阵；[·]†表示矩阵Moore-Penrose逆；A≥B表示A-B为半正定矩阵；＜a＞m表示向量a中的第m个元素。

1 TDOA/FDOA定位观测模型

现有M个运动传感器利用TDOA/FDOA观测信息对N个非相关运动辐射源进行定位。第m个传感器位置和速度向量分别记为sm=并称其为第m个传感器位置-速度向量；第n个辐射源位置和速度向量分别为令并称其为第n个辐射源位置-速度向量。在辐射源定位场景中，TDOA和FDOA可分别等价为距离差和距离差变化率，下文直接利用距离差和距离差变化率进行建模和分析。

不失一般性，将第1个传感器设为参考，将第n个辐射源与第m个传感器之间的距离和距离变化率分别记为于是第n个辐射源与第m个和参考传感器之间的距离差和距离差变化率可以分别表示为

实际中获得的距离差和距离差变化率观测量都含有误差，它们可以分别表示为

式中：和表示距离差和距离差变化率观测误差。针对第n个运动辐射源，将式（2）中的2组等式写成向量形式可得

现将式（3）中的2个等式合并可得

式中：

假设观测误差向量服从均值为零、协方差矩阵为的高斯分布。

虽然辐射源互不相关，但是不同辐射源对应共同的传感器位置-速度误差，此时应考虑对多个辐射源进行协同定位，以获得协同增益。为此需要将式（5）中的N个等式合并可得

式中：

观测误差向量服从均值为0、协方差矩阵为}的高斯分布。

在实际定位中，当传感器安装在机载或舰载平台，又或是传感器随机布设时，传感器位置-速度向量可能无法精确获得，仅能得到其先验观测值，即

注释1在多源协同定位前需要首先估计辐射源个数，这属于信号处理问题。如果是利用单天线（单通道）估计信号的TDOA/FDOA参数，则可以通过信号自相关矩阵的奇异值来确定信号个数［36］，或者构造信号时频图像，并通过卷积神经网络检测信号个数［37］。如果是基于阵列天线估计信号的TDOA/FDOA参数，则可以利用信息论准则（例如AIC准则、MDL准则等）获得信号个数［38-39］，或者利用深度学习技术确定信号个数［40］。

2 新方法的基本原理与计算步骤

2.1 新方法的总体描述

为了避免发散和局部收敛等问题，本节提出一种新的定位方法，该方法是由多维标度分析衍生出，相比其他定位方法，其在大观测误差条件下具有更高定位性能。为获得渐近统计最优估计精度，新方法包含2个计算阶段（记为阶段a和阶段b），每个阶段的主要原理可描述如下：

1）阶段a首先基于多维标度分析构造2组标量积矩阵，并利用标量积矩阵性质构建定位关系式，然后通过一阶误差分析推导定位关系式中的误差协方差矩阵，并利用矩阵正交变换将其恢复为满秩矩阵，从而确定加权矩阵，最后基于定位关系式和加权矩阵获得多个非相关辐射源以及各个传感器位置和速度估计值。

2）阶段b首先基于多维标度分析中引入的中间变量，构建关于阶段a定位误差的约束优化模型，然后基于此模型对阶段a的定位误差进行估计，最后利用该估计值对阶段a的定位结果进行更新，以获得渐近统计最优的估计性能。

2.2 阶段a的计算原理

2.2.1 构造标量积矩阵

在多维标度分析中需要首先构造标量积矩阵，由于本文考虑多源协同定位，因此需要依次对每个辐射源构造标量积矩阵。首先针对第n个辐射源定义如下4维复坐标和复速度向量

式中：j表示虚数单位；dn1=rn1-rn1=0和=-=0。基于式（11）可定义如下复坐标和复速度矩阵

2.2.2 构建定位关系式

下面利用标量积矩阵构建新的定位关系式，与其它常规定位方法中的关系式不同，这里的关系式能充分利用标量积矩阵的维度信息和特征结构信息，从而提高在大观测误差条件下的定位性能。首先由第1组标量积矩阵{Wn}1≤n≤N可以得到如下结论。

命题1针对第n个辐射源，分别定义如下实向量vn和实矩阵Vn

证明首先由Moore-Penrose逆矩阵的性质可得

利用矩阵和求逆公式可知

然后将式（19）代入式（18）中可得

另一方面，将式（15）和式（16）代入式（12）中的第1式中可知

最后联合式（13）、式（20）以及式（21）可得

证毕。

若定义矩阵

则由命题1可知

将式（24）两边对时间求导可得

式中

最后将式（24）和式（25）合并可知

若忽略误差二阶项，则由式（28）可得

式中：ΔWn、、ΔTn以及分别表示、、以及中的误差一阶项。

利用一阶误差分析可以将式（29）右边第1项表示为关于观测误差向量和ε(s)的线性函数，即

式中：

利用一阶误差分析可以将式（29）右边第2项表示为关于观测误差向量ε(d)n和ε(s)的线性函数，即

利用一阶误差分析可以将式（29）右边第3项表示为关于观测误差向量ε(d)n和ε(s)的线性函数，即

式中：

其中

利用一阶误差分析可以将式（29）右边第4项表示为关于观测误差向量ε(d)n和ε(s)的线性函数，即

式中

其中

将式（30）、式（34）、式（38）以及式（42）代入式（29）中可得

式中：由式（46）可知，误差向量ξn渐近服从零均值高斯分布，并且其协方差矩阵为

式中：Qn1和Qn2均为列正交矩阵；Rn为上三角满秩矩阵。为了得到协方差矩阵为满秩的误差向量，可以将矩阵左乘以误差向量ξn，并且结合式（28）、式（46）以及式（48）可知

由式（49）可知，误差向量δn的协方差矩阵为

由于误差向量δn的维数为2(M-1)，其与TDOA/FDOA观测量个数相等，因此容易验证COV(δn)为满秩矩阵。

2.2.4 估计准则及其最优解

首先将矩阵和按列分块表示为=由式（49）可将误差向量δn重新写为

式中：

从式（49）中可以看出，N个误差向量{δn}1≤n≤N中包含共同的模型误差向量ε(s)，为了进行多辐射源协同定位，需要将式（51）中的N个等式合并，如下式所示

式中：

结合式（49）和式（53）可知，误差向量δc渐近服从零均值高斯分布，并且其协方差矩阵为

式中：

为了提高发生门限效应时的误差阈值，这里考虑对向量和ˉ进行联合估计，为此需要定义扩维参数向量和扩维误差向量然后结合式（10）和式（53）可构建阶段a中的估计准则，即

式中：

式（57）中的(COV(ζ))-1可以看成是加权矩阵，其作用在于抑制观测误差向量和模型误差向量ε(s)的影响，其最优解为

若将估计值（或者{}1≤m≤M）中的误差向量记为（或者{}1≤n≤N）以及则有

2.2.5 理论性能分析

阶段b是对阶段a的估计结果进行更新，其需要阶段a估计结果的均方误差矩阵，因此下面将对阶段a的估计结果进行理论性能分析。首先由线性加权最小二乘估计理论可知，误差向量Δρa渐近服从零均值高斯分布，并且估计值的均方误差矩阵为

式中：Pc=blkdiag{P1，P2， … ，PN}，其中

接着利用分块矩阵求逆公式可知

将式（64）代入式（62）中可得

式中：

结合式（61）和式（65）可知，误差向量渐近服从零均值高斯分布，并且估计值的均方误差矩阵为

式中需要指出的是，虽然向量是渐近无偏估计值，但尚不是渐近统计最优估计值，具体可见如下命题。

命题2若将未知参数的估计均方误差的CRB记为则在一阶误差分析理论框架下满足

命题2的证明见附录A。阶段a的估计值不具备渐近统计最优性的根本原因在于，中间变量的引入导致向量ρ中的元素间存在相关性，这种相关性会使得误差向量Δρa服从某个等式约束。阶段b将利用该等式约束进一步提高定位精度。

2.3 阶段b的计算原理

首先推导误差向量Δρa服从的等式约束。根据向量ρ第8n-4个元素和第8n个元素的定义可知

分别将式（68）和式（69）右边第1项在点处进行一阶Taylor级数展开可得

式中：

将式（70）和式（71）中的2N个等式合并，可以得到如下关于误差向量Δρa的线性等式

式中：

式（73）即为误差向量Δρa满足的等式约束。根据误差向量Δρa服从的高斯分布特性，可以建立估计误差向量Δρa的约束优化模型，即

利用拉格朗日乘子技术可知［46］，式（75）的最优解为

结合式（61）和式（76）可以得到误差向量和的估计值分别为

于是辐射源位置-速度向量和传感器位置-速度向量ˉ在阶段b的估计值分别为

2.4 新方法的计算步骤与讨论

综合前面的讨论，总结新方法的计算步骤，如图1所示。针对上述新方法给出3点注释：

图1 本文新方法的计算流程图
Fig.1 Flow chart of method proposed in this paper

3 新方法的理论性能分析

本节推导估计值的统计特性，主要推导其估计均方误差，并与相应的CRB进行比较，从而证明其渐近统计最优性，为此首先给出多个非相关源协同定位的CRB的表达式。

3.1 多个非相关源协同定位的CRB

基于文献［34-35］中的结论可知，多个非相关源协同定位的CRB矩阵可以表示为

式中：均表示函数fc(ˉ)的Jacobian矩阵，它们可以分别表示为

其中：均表示函数f(ˉ)的Jacobian矩阵。由式（79）和式（80）可知，多个非相关源协同定位的CRB与辐射源位置紧密相关。

3.2 估计均方误差

首先将式（75）中的等式约束代入式（76）中可得

式中由式（81）可知，估计值的均方误差矩阵为

结合式（77）和式（82）可知，估计值的均方误差矩阵为

3.3 渐近统计最优性分析

下面证明估计值具有渐近统计最优性，为此需要给出均方误差矩阵的另一种表达式，具体可见如下命题。

命题3均方误差矩阵可以表示为

证明：首先由式（83）可得

考虑矩阵这2个矩阵的列数之和等于8N+6M，并且满足

于是有

由此可知

将式（88）代入式（85）中可得

不难验证，将该式代入式（89）中可知式（84）成立。证毕。

命题4在一阶误差分析理论框架下，满足

证明：首先根据向量ρ的定义可得

将式（90）代入式（84）中，并且利用式（65）可知

结合式（66）和式（A8）可得

将式（92）～式（94）代入式（91）中，并且结合式（79）可知命题4成立。证毕。

4 仿真实验

本节进行仿真实验，基础实验条件为：利用7个运动传感器获得TDOA/FDOA信息（亦即距离差/距离差变化率信息）对运动辐射源进行定位，传感器位置和速度见表1，距离差/距离差变化率观测误差向量服从均值为零、协方差矩阵为的高斯分布，传感器位置-速度先验观测误差向量ε(s)服从均值为零、协方差矩阵为E(s)=σ22(IM⊗blkdiag{I3，I3/10})的高斯分布。

表1 传感器位置和速度数值
Table 1 Numerical value of sensor position and velocity

传感器序号1234567 x(s)m/m 2 400 1 700 2 300-1 600 2 500-1 100-2 200 y(s)m/m 1 800 2 100-1 900 1 700-1 600-1 400 1 500 z(s)m/m 1 300-2 200 1 600 1 500-1 200 2 400-2 500 ẋ(s)m/（m·s-1）-10 10-18-10-10 18 15 ẏ(s)m/（m·s-1）16-15-10 10 18-15-10 ż(s)m/（m·s-1）15 16 18-18 15-16 16

仿真实验1假设有2个非相关辐射源，辐射源1位置和速度向量分别为u1=[-3 300，-3 600，3 800]Tm和=[20，10，-5]Tm/s，辐射源2位置和速度向量分别为u2=[7 200，7 500，-7 700]Tm和u̇2=[-10，15，15]Tm/s，令20lg(σ1)=20lg(σ2)=-5 dBm，利用本文方法对这2个辐射源进行协同定位，并且进行5 000次蒙特卡罗实验，图2和图3给出了辐射源1和辐射源2的定位结果散点图与定位误差椭圆曲线。

图2 辐射源1的定位结果散点图与定位误差椭圆曲线（x-y平面）
Fig.2 Scatter plot of positioning results and elliptic curve of positioning error for Emitter 1 （x-yplane）

图3 辐射源2的定位结果散点图与定位误差椭圆曲线（y-z平面）
Fig.3 Scatter plot of positioning results and elliptic curve of positioning error for Emitter 2 （y-zplane）

从图2和图3中可以看出，本文方法的定位结果散点图的形状与定位误差椭圆的形状一致，并且大概率对应大椭圆，小概率对应小椭圆，从而验证了该方法的有效性。

仿真实验2假设有3个非相关辐射源，辐射源1位置和速度向量分别为u1=[4 600，5 600，5 800]Tm和=[10，-25，20]Tm/s，辐射源2位置和速度向量分别为u2=[6 200，7 200，7 400]Tm和=[-20，15，10]Tm/s，辐射源3位置和速度向量分别为u3=[8 800，9 600，9 200]Tm和=[15，-10，-10]Tm/s，下面利用本文方法对这3个辐射源进行协同定位，并将其定位精度与一些主流TDOA/FDOA定位方法进行比较，其中包括文献［27］中的TLS方法、文献［31］中的TSWLS方法、文献［34］中的TSWLS方法、文献［35］中的改进型TSWLS方法以及文献［13］中的ICWLS方法，其中前面4种是闭式定位方法，最后1种是迭代定位方法。需要指出的是，虽然文献［27］中的TLS方法、文献［31］中的TSWLS方法以及文献［13］中的ICWLS方法都是在单辐射源存在的场景下提出的，但却可以将它们推广至多辐射源协同定位的场景。此外，文献［31］中的TSWLS方法联合利用了DDR信息进行定位，但为了公平比较，这里的仿真并未利用此信息。另一方面，选择ICWLS方法作为比较对象是因为它是一种性能优越的迭代定位方法［13］，为了进行更全面的比较，这里给出其初始值为随机值和真实值2种情形下的定位性能。首先令20lg(σ2)=0 dBm，改变σ1的数值，图4给出了多辐射源位置和速度估计均方根误差随20lg(σ1)的变化曲线；图5给出了多传感器位置和速度估计均方根误差随20lg(σ1)的变化曲线。然后令20lg(σ1)=-5 dBm，改变σ2的数值，图6给出了多辐射源位置和速度估计均方根误差随20lg(σ2)的变化曲线；图7给出了多传感器位置和速度估计均方根误差随20lg(σ2)的变化曲线。

图4 多辐射源位置和速度估计均方根误差随20lg（σ1）的变化曲线
Fig.4 RMSE of multiple-source position and velocity estimation versus 20lg（σ1）

图5 多传感器位置和速度估计均方根误差随20lg（σ1）的变化曲线
Fig.5 RMSE of multiple-sensor position and velocity estimation versus 20lg（σ1）

图6 多辐射源位置和速度估计均方根误差随20lg（σ2）的变化曲线
Fig.6 RMSE of multiple-source position and velocity estimation versus 20lg（σ2）

图7 多传感器位置和速度估计均方根误差随20lg（σ2）的变化曲线
Fig.7 RMSE of multiple-sensor position and velocity estimation versus 20lg（σ2）

从图4～图7中可以看出：① 在大观测误差条件下本文方法的估计精度优于TLS方法和3种TSWLS方法，新方法发生门限效应的误差阈值更高，这验证了该方法相比其他闭式定位方法的优越性；② 在大观测误差条件下，本文方法估计精度优于初始值为随机值的ICWLS方法，并且与初始值为真实值的ICWLS方法接近，注意到将真实值作为初始值在实际计算中难以实现，而本文方法无需初始值，这验证了本文方法相比迭代定位方法的优越性；③ 本文方法对辐射源定位性能渐近逼近相应的CRB，这验证了第3节理论性能分析的有效性；④ 本文方法可以提高传感器位置和速度估计精度（相比其先验观测精度而言），并且对传感器位置和速度估计均方根误差均可以渐近逼近相应的CRB，从而再次验证了第3节理论性能分析的有效性。

仿真实验3仍然考虑3个非相关辐射源存在的场景，辐射源1位置和速度向量为u1=[5cos(γ)cos(α)，5sin(γ)cos(α)，5sin(α)]Tkm和=[10，15，-10]Tm/s，辐射源2位置和速度向量为u2=[6cos(γ)cos(α)，6sin(γ)cos(α)，6sin(α)]Tkm和=[-10，5，-15]Tm/s，辐射源3位置和速度向量为u3=[7cos(γ)cos(α)，7sin(γ)cos(α)， 7sin(α)]Tkm和=[-10，15，15]Tm/s。由于本文方法可以通过简化直接应用于单辐射源存在场景（此时退化为非协同定位方法），于是下面不妨比较多辐射源协同定位和多辐射源非协同独立定位2种方式下的定位精度，旨在验证多源协同定位能产生协同增益。此外，还将文中方法与文献［30］中的加权多维标度定位方法进行比较，这是因为该方法也考虑了传感器位置和速度先验观测误差，只是该方法是针对单辐射源的定位方法，目前尚没有进行多源TDOA/FDOA协同定位的加权多维标度定位方法。另一方面，文献［25］也是基于加权多维标度的TDOA/FDOA定位方法，但是该方法没有考虑传感器位置和速度先验观测误差，因此不适合与文中的新方法进行性能比较。令20lg(σ1)=-10 dBm，20lg(σ2)=-5 dBm以及γ=12°，通过改变α的数值来改变辐射源位置。图8给出了辐射源1的位置和速度估计均方根误差随α的变化曲线；图9给出了辐射源2的位置和速度估计均方根误差随α的变化曲线；图10给出了辐射源3的位置和速度估计均方根误差随α的变化曲线。

图8 辐射源1的位置和速度估计均方根误差随α的变化曲线
Fig.8 RMSE of position and velocity estimation of Emitter 1 versusα

图9 辐射源2的位置和速度估计均方根误差随α的变化曲线
Fig.9 RMSE of position and velocity estimation of Emitter 2 versusα

图10 辐射源3的位置和速度估计均方根误差随α的变化曲线
Fig.10 RMSE of position and velocity estimation of Emitter 3 versusα

从图8～图10中可以看出：① 多源协同定位产生的性能增益是显著的；② 本文方法（非协同定位版）的估计精度优于文献［30］中的加权多维标度定位方法，这是因为后者尚未利用中间变量包含的有用信息；③ 随着辐射源位置的改变，文献［30］中的加权多维标度定位方法的估计性能并不总能达到相应的CRB，而本文方法包含2个计算阶段，可以使其定位精度始终渐近逼近相应的CRB。

仿真实验4将本文方法与文献［27］中的TLS方法、文献［31］中的TSWLS方法、文献［34］中的TSWLS方法、文献［35］中的改进型TSWLS方法以及文献［13］中的ICWLS方法的运行时间进行比较，以间接比较6种定位方法的计算复杂度。仿真软件为MATLAB R2020a，仿真程序在安装有i7-3520CPU的PC机上运行，令20lg(σ1)=10 dBm和20lg(σ2)=10 dBm，其余条件不变，并进行5 000次蒙特卡罗实验。表2给出了6种定位方法的平均运行时间。从表2中可以看出，新方法的计算复杂度高于4种闭式定位方法，但是低于ICWLS定位方法。

表2 6种定位方法的平均运行时间
Table 2 Average running time of six localization methods

定位方法本文方法文献［27］方法文献［31］方法文献［34］方法文献［35］方法文献［13］方法平均运行时间/s 1.374 2 0.329 2 0.663 9 0.664 2 0.681 6 3.274 2

5 结论

为了提高在大观测误差条件下的定位性能，该文提出了一种基于加权多维标度分析的多个非相关辐射源TDOA/FDOA协同定位方法，主要结论包括：

1）在大观测误差条件下，本文方法的估计精度优于其他闭式定位方法，其具有更高的误差阈值。

2）在大观测误差条件下，本文方法的估计精度与迭代初始值为真实值的迭代定位方法接近（以性能较优的ICWLS定位方法为比较对象）。

3）本文方法的估计精度优于已有的加权多维标度定位方法。

4）本文方法可以进一步提高传感器位置和速度的估计精度（相比其先验观测精度而言）。

5）本文方法通过多源协同定位获得了较高的协同增益。

附录A首先定义如下向量和矩阵

式中：

于是有进一步由式（27）可得

接着将定义式代入式（A3）中的矩阵Wn、、Tn以及中可以得到关于向量和ˉ的恒等式，将该恒等式两边分别对向量和求导可知

将式（48）代入式（A4）和式（A5）中可得

分别将式（A6）和式（A7）中的N个等式合并可知

将式（A8）代入式（79）中，并且利用式（66）可得

根据矩阵J的定义可以验证JI6N+6M，由此可以将矩阵进一步表示为

最后对比式（67）和式（A10），并且利用矩阵不等式A-1≥B(BTAB)-1BT（其中A为任意正定矩阵，B为任意列满秩矩阵）可知命题2成立。

附录B

附录B中以实数乘法为参照给出了新方法的计算复杂度，具体见表B1。

表 B1 新方法的计算复杂度
Table B1 Computational complexity of proposed method

步骤步骤a-1步骤a-2步骤a-3步骤a-4步骤a-5步骤a-6步骤a-7步骤a-8步骤b-1步骤b-2步骤b-3步骤b-4各步骤计算量M(M+1)N(107M+O(53))N 15M2N 128MN3+16MN2+64N2(12M4+438M2+225M+225+O(53))N 16M(M-1)2N(16(M-1)3+24M2(M-1)+48M(M-1)2+O(8(M-1)3))N 32(M-1)2N3+128(M-1)N3+96M(M-1)N2+4M(M-1)N+(8N+6M)(6M+2(M-1)N)+O((6M+8N)3)32(M-1)2N3+128(M-1)N3+96M(M-1)N2+O((6M+8N)3)24N 12N 2(8N+6M)2N+(8N+6M)×(4N2+2N)+4N2+O(8N3)总计算量12M4N+64M2N3+192M2N2+256MN3+104M3N+56MN2+342M2N-96N3+68N2+36M2+489MN+309N+NO(8(M-1)3)+2O((6M+8N)3)+2NO(53)+O(8N3)