未知发射机位置的闭式解椭圆定位方法∗

2024-01-26郑瑞超

传感技术学报 2023年12期

白璐，郑瑞超，王刚

(宁波大学信息科学与工程学院，浙江宁波 315211)

目标定位作为一个经典的研究课题，近几十年以来一直备受关注。传统的目标定位在军事领域主要采用雷达和声纳来定位敌对目标，同时为己方目标提供协助。然而，随着半导体芯片技术、时钟同步技术、移动通信技术和计算机技术等现代技术的发展和成熟，目标定位逐渐用于民用领域，如海洋监测、医疗救援、交通控制和灾害预防等[1－3]，并日渐成为人类生产生活中不可或缺的一部分。

目标定位通常需要布置一些位置已知的传感器，通过各个传感器与目标之间的通信收集有效信息，然后运行定位算法，以此完成定位。这些信息通常包括到达时间(Time of Arrival，TOA)[4－6]、到达时间差(Time Difference of Arrival，TDOA)[7－9]以及功率信息(如接收信号强度和声音能量)[10－13]等。其中，基于TOA 的定位方法是通过寻找若干圆(二维)或球(三维)的交点来估计目标位置，所以它又称为圆(球)形定位，而基于TDOA 定位的方法则是通过寻找若干条双曲线的交点来估计目标位置，因而又被称为双曲线定位。在多输入多输出雷达及无源相干定位系统中，目标定位通常由接收机观测发射机发出并经过目标反射的信号来完成，充分利用接收机接收到的信号可以完成测距，这些测距信息可以形成以接收机和发射机位置为焦点的椭圆(椭球)，这种定位方式通常称为椭圆(椭球)定位。本文研究的正是椭圆(椭球)定位，它需要发射机首先发出信号，然后通过几个发射机－接收机对的测量传播时间分别形成的椭圆的交点来确定目标位置。现有研究表明[14－15]，在接收机和发射机之间时钟同步的情况下，椭圆定位具有比TDOA 定位更高的精度。具体来说，文献[14]通过克拉美－罗下界(Cramer-Rao Lower Bound，CRLB)分析证实了在接收机和发射机具备时钟同步的条件下，椭圆定位具有更低的CRLB。同时，文献[15]通过CRLB 分析表明，当发射机位置未知但发射机与接收机同步时，若使用直接路径信息对发射机和目标位置进行联合估计，也能够取得更高的目标定位精度。因此，椭圆定位被广泛应用于对空中无人机及水下目标等的定位，具有重要的应用价值。

接收机接收到的信号一般有两种传播途径，一种是发射机的发射信号经过目标反射后再被接收机接收到，称为间接路径传播；另一种则是不经目标反射直接被接收机接收，称为直接路径传播，在这种情况下，信号从发射机到接收机的传播过程不会携带任何关于目标位置的信息，因此，当发射机的位置完全已知时，只需获得间接路径的测量信息就可以完成目标定位。然而，当发射机的位置随着时间变化或未知时，比如发射机位于导航系统无法通信的地区，或某些刻意隐藏发射机位置以减轻硬件要求的场景，在这些情况下，目标位置的估计就会存在延时，从而使得估计出现偏差，甚至无法估计出真实的目标方位。为解决这一问题，文献[15]提出了两步加权最小二乘方法(Two-Step Weighted Least Squares，TSWLS)来联合估计目标和发射机的位置。该方法在大噪声下能够达到克拉美－罗下界，但同时会出现门限效应，使得性能迅速恶化。文献[16]提出了一种半正定规划(Semidefinite Programming，SDP)方法来联合估计目标和发射机的位置，该方法在大噪声环境下具有稳定的性能。然而，半正定松弛在小噪声下可能使半正定规划的解无法达到原问题的最优解，导致其结果偏离CRLB。

为克服现有方法中存在的问题，本文中提出了一种三步估计法，其思想是将发射机和目标位置分别估计。首先，通过对直接路径模型处理，将该问题描述为广义信赖域子问题(Generalized Trust Region Subproblem，GTRS)，以此得到发射机的初始位置估计。然后，将第一步估计的发射机位置代入间接路径模型并把该问题表述为关于目标位置的GTRS 问题，进而可以完成对于目标位置的初始估计。最后，本文通过构造线性加权最小二乘问题来联合估计目标和发射机的位置误差项，以进一步改善估计性能，从而得到精确的目标和发射机位置估计。仿真结果表明，所提出方法的性能能够在使用较少数量的接收机的情况下，即可达到克拉美－罗下界，并且与现有方法相比复杂度更低。

此处对本文所使用的符号进行统一说明。粗体的小写字母和大写字母分别表示列向量和矩阵。a(1:k)表示向量a的第1 到第k个元素组成的向量。tr‖A‖是A的迹，diag(a)和diag(A)为a的元素构成的对角矩阵和由A的对角元素构成的向量。blkdiag(A，B)表示A和B在对角线上的块对角矩阵。0k是k×1 维全零向量，Ik表示维度为k的单位矩阵。‖A‖和E[A]分别表示A的2 阶范数和数学期望。ao是变量a的真实值。

1 测量模型

考虑k维空间中的椭圆定位问题(图1)，发射机和目标的位置分别用∈Rk和∈Rk表示，接收机的位置分布如图所示，用，j＝1，…，N表示。假设发射机和接收机在时间上是同步的，且每个接收机将已知的发射机信号波形与两条路径的叠加信号相关联。考虑到发射信号有足够的带宽和持续时间，两个最强的峰值可通过互相关运算分离，其位置对应于直接和间接路径的传播时间。由于直接路径信号的传播距离较短，较小的一个是直接路径传播时间，另一个则是间接路径传播时间。

图1 椭圆定位模型

由于在信号传播速度已知的情况下，时间延迟和距离是成比的，本文在不引起歧义的情况下将它们互换使用。间接路径的真实距离可表示为:

式中:表示发射机到目标和目标到第j个接收机之间的距离之和。直接路径的真实距离为:

式中:表示发射机到第j个接收机之间的距离。

在实际中只能利用其带有误差的测量值，间接路径和直接路径的测量距离值可表示为:

式中:rj表示间接路径的测量值，εr，j表示间接路径测量噪声，dj表示直接路径的测量值，εd，j表示直接路径测量噪声。

2 三步闭式解方法

本文提出了一种三步估计方法来估计目标和发射机的位置。首先根据直接路径测量模型式(4)构造关于发射机位置的GTRS 问题，求解该GTRS 得到发射机的位置估计，然后将第一步估计的发射机位置代入间接路径测量模型式(3)中构造关于目标位置的GTRS，从而可以得到目标的位置估计。最后，通过构造线性加权最小二乘问题来联合估计目标和发射机位置的估计误差以得到更精确的发射机和目标位置估计。

2.1 第一步:发射机位置估计

对于直接路径测量式(4)两边进行平方并且忽略二阶噪声项得:

定义优化变量＝[toT，‖to‖2]T，将式(5)表示成伪线性方程的形式:

根据式(6)，可以构造如下约束加权最小二乘问题:

式中:Ik表示维度为k的单位矩阵。注意到问题(10)中的约束明确表示出了向量yd中元素之间的关系。

问题(10)是一个典型的GTRS，其虽然为一个非凸问题，但该问题可以求得全局最优解，求解方法参考文献[17－18]。根据文献[17－18]，式(10)的全局最优解具有如下形式:

式中:λ1表示拉格朗日乘子。由式(12)可看出，若获得λ1的值，GTRS 的最优解即可获得。根据文献[15－16]，λ1可通过求解如下非线性方程获得:

式中:H表示λ1的区间范围，具体表示如下:

式中:λd(∗)表示矩阵∗的最大特征值。通常，可采用二分法求解非线性方程式(13)。

由式(12)可以得到发射机的位置估计:

值得注意的是，权重矩阵Wd＝(BdQd)－1依赖于未知的发射机真实位置to。在算法实现时，本文首先设置矩阵Bd为单位矩阵即Wd＝(Qd)－1求得发射机位置的初始估计，然后用初始估计值更新权重矩阵Wd，再一次求解式(12)得到最终估计。

2.2 第二步:目标位置估计

在第二步中，本文将利用间接路径模型估计目标位置。将第一步中发射机位置估计写为＝to＋Δt，其中Δt是发射机位置估计误差，并将＝to＋Δt代入模型(3)得:

对式(16)两边平方并忽略二阶噪声项得:

为了便于运算，定义变量＝[uoT，‖uo－‖]T，由此可将式(17)表示成伪线性方程的形式:

式中:向量br为:

矩阵Ar、Br和Er表达式如下

根据式(18)可构造如下约束加权最小二乘问题:

式(23)为一个GTRS。因此，通过类似于2.1 节中的求解方法即可得到目标位置估计。令表示优化变量yr的估计，其可表示为:

式中:λ2表示拉格朗日乘子，且其是下列方程的唯一解:

由式(25)可以得到目标位置估计:

权重矩阵Wr依赖于协方差矩阵和真实目标位置。一方面，根据文献[16]，在大噪声时，第一步中的发射机估计为近似无偏且趋近于CRLB。因此，协方差矩阵Qt可使用CRLB 近似，其中，CRLB 表达式中的发射机真实位置可用发射机位置估计代替，即:

另一方面，类似2.1 节，本文首先通过设置矩阵Bd为单位矩阵即Wr＝(Qr)－1，求得目标位置的初始估计，然后用初始估计值更新一次权重矩阵Wr，再一次求解式(25)得到最终估计。

2.3 第三步:发射机和目标位置估计的进一步精确

在前两步估计中，本文依次估计发射机和目标的位置。当直接路径和间接路径测量噪声存在相关性时，分步估计将导致一定的性能损失而无法达到联合估计的CRLB。本节通过构造线性加权最小二乘问题以进一步提高发射机和目标位置的估计精度，其思想是联合估计发射机和目标的位置误差，并利用估计误差校正前两步的估计，使其能够趋近于CRLB。

本文的思想是构造关于Δu和Δt的线性加权最小二乘问题。为此，分别将式(28)和式(29)进行一阶泰勒展开，得:

令Δφ＝[ΔuT，ΔtT]T，可将式(30)和式(31)表示成如下线性方程:

式中:矩阵A＝[，]T，向量b的表达式为:

由于式(32)为关于Δφ的线性方程，其线性最小二乘估计为:

式中:W＝Q－1表示权重矩阵。

3 分析

3.1 估计性能分析

在本节中，证明所提方法的估计均方误差(Mean Square Error，MSE)在大噪声条件下能够趋近于CRLB。

将式(32)代入式(36)，误差项δ＾φ＝[δuT，δtT]T可通过以下计算得到:

根据式(38)，易知δ＾φ的MSE 为:

显然，A中包含误差项，可令A＝Ao＋ΔA。当噪声足够小时，估计误差ΔA二次项可忽略。由此，A＝Ao＋ΔA代入上式并仅保留一阶误差项，有:

根据Neumam[19]公式可得:

由于E[ΔA]＝0，所以式(41)第二项取期望后为0，即:

根据文献[15]，不难验证，CRLB 的表达式为:

因此，在足够小的噪声下，最终估计的均方误差能够趋近于CRLB。

3.2 复杂度分析

本文所提出方法需要求解两次GTRS 及一次线性加权最小二乘问题，主要计算复杂度为O(2N2m＋T1m3＋T2m3)[17]，其中N是接收机的数量，m＝k＋1，T1、T2分别为两步GTRS 问题的迭代次数。对于二维定位问题，GTRS 中的λ1和λ2为一元四次方程的根，有闭式表达式，因此不需要迭代，即T1＝T2＝1，复杂度为O(2N2m＋2m3)；对于三维定位问题，λ1和λ2为一元六次方程的根，需要进行迭代求解，复杂度为O(2N2m＋T1m3＋T2m3)，T1、T2与搜索区间相关，以T1为例，设其搜索区间为[p，q]，则T1满足T1>log2[(qp)/γ]的最小整数，其中γ为求解精度。文献[16]中SDP 方法的计算复杂度为O{4N2(m＋2)＋(2m＋5)0.5×[5(2m＋3)3＋25(2m＋3)2＋10m2＋125]}。相比之下，本文方法的计算复杂度低于SDP 方法。

4 仿真结果

本节将通过数值仿真来评估所提出的方法的性能。简易起见，本节只考虑三维情况，因为二维场景的观察结果是类似的。同时，给出TSWLS[15]、SDP[16]以及最大似然估计(Maximum Likelihood Estimator，MLE) 的性能作为比较，其中MLE 是由MATLAB 函数“fminunc”以真实的发射机和目标位置作为初始值迭代求解。通过MSE 来评估所提出方法的性能，其定义如下:

式中:tr{Qo}表示Qo的迹，σ2表示测量噪声功率。τ的定义为，其中，γ表示信号经过目标反射后的噪声损失，本文设置γ＝3 dB。

4.1 使用四个接收机

考虑在三维空间中定位一个目标，其位置坐标uo＝[－1 000，500，1 500]Tm，发射机和接收机的位置坐标是从(－4 000，4 000)×(－4 000，4 000)×(1 000，3 000)m3区域中随机产生。使用一个发射机和四个接收机来评估所提出方法的性能，噪声设置为相关噪声。图2 和图3 分别给出了随噪声水平变化时目标和发射机位置估计仿真结果。由图中可以看出，在小噪声下，本文所提出的方法能够趋近于CRLB。而由于松弛导致SDP 即使在小噪声下也无法保证全部产生秩为1 的解，其性能低于本文所提出的闭式解方法。而在大噪声下，虽然所提方法的性能略低于MLE 和SDP，但其复杂度低于这两种方法。而SDP 方法虽然性能稳定，但在接收机数量较少的情况下，由于松弛导致SDP 获得秩为1 的解较少，因此在大噪声下也无法趋近于CRLB。注意到在大噪声下MLE 的均方误差低于CRLB，这是由于MLE 问题的非凸性，无法保证取得全局最优解而可能收敛到真实位置附近的局部点，从而导致最大似然估计是有偏估计，进而使MLE 的性能甚至低于CRLB。最后，由于接收机的数量较少，无论是目标位置还是发射机位置，无论在小噪声还是大噪声下，TSWLS 方法均无法精确估计。

图2 使用4 个接收机时不同方法对于目标位置估计MSE 比较

图3 使用4 个接收机时不同方法对于发射机位置估计MSE 比较

4.2 使用五个接收机

由4.1 节知，当使用四个接收机时，TSWLS 方法无法正常工作，因此，在本节增加一个接收机，即使用五个接收机去评估所有方法的性能。增加的一个接收机位置从相同的空间中随机产生，其他参数设置均与4.1 节相同。图4 和图5 分别给出了目标和发射机位置的估计结果。由图可得，当使用五个接收机时，TSWLS 性能有了很大的提升，但仍不理想。尽管SDP 在大噪声时的性能略优于本文所提出的方法，但在小噪声时的性能仍然无法趋近于CRLB。在性能接近的情况下，本文所提出的方法可应用于计算能力不足时的场景。由于为有偏估计，MLE 在大噪声时仍然略低于CRLB。