二维TDOA自校准新解法

传感器网络TDOA自校准：二维复杂性分析与解决方案

1. 引言

无线传感器网络已被广泛研究[1, 2,3, 4]，并已成功应用于多个领域，如定位[5],建图[6],麦克风阵列校准[7]和波束成形[8]。为了正确使用，必须首先对网络进行校准，即确定其节点的位置。这可以通过例如在已知位置设置发射机并对网络节点进行三边测量来实现。

然而，一些应用需要同时对接收器和发射器进行定位[9, 10]。这被称为网络自校准，也是本文的主要主题。

考虑一个包含 m个接收器和 n个发射器的网络，简称为 mr/ns。根据测量结果，将有mn个方程可用于求解自校准问题。现在，我们考虑两种不同的场景：

• 同步的接收机和发射机：在这种情况下，信号的发送和测量时刻是已知的，因此可以测量接收机与发射机之间的距离。对于每一对接收机‐发射机，得到以下多项式方程
f2 ij= d2ij= ‖ri − sj‖ 2, (1)
其中 fij为测得的距离， dij为位于ri的接收机与位于sj的发射机之间的距离。总共我们将有 K(m+ n) − G自由度（DoF），其中 K表示空间维度（K= 2,3）， G为规范自由（2D中为 G= 3，3D中为 G= 6）。由于我们仅测量距离，位置只能在欧几里得变换意义下恢复。实际上，这意味着求解坐标时所采用的坐标系可以自由选择，从而减少自由度。该模型称为到达时间（TOA）[11]。

• 同步的接收机与不同步的发射机：在这种情况下，所有接收器将在同一时钟帧中测量到达时间。对于每一对发射机‐接收器，可得到以下多项式方程
(fij − o j) 2= d2ij = ‖ri − s j‖ 2. (2)
如方程(2)所示，每个发射机引入了一个额外的未知量 o j，即发射机本地时钟与接收器时钟之间的时间偏移。因此，自由度的总数为 K(m+ n)+ n− G。这种表述方式称为Time-Diff erence-Of-Arrival（到达时间差）(TDOA)[12]。

上述方程可通过数值迭代方法[13, 14]求解。然而，这些问题可能存在若干问题，例如陷入局部最小值、收敛缓慢以及对异常值敏感。研究表明，代数非迭代方法能够实现更高的精度[11, 12, 15, 16, 17]。此外，将代数解作为迭代方法的初始值，可实现更快更精确的收敛。

本文重点研究二维到达时间差(2D‐TDOA)（K= 2, G= 3），即我们同时确定接收器和发射机的位置以及发射机的时间偏移。总的自由度数量为 2m+ 3n − 3，当自由度数量等于方程数量时的配置情况

这些配置被称为最小构型，而使用最少数量的接收器和发射器进行网络自校准的问题被称为最小问题。二维到达时间差的最小构型是6r/3s和4r/5s。以往的研究已经为非最小情况开发了代数求解器，例如7r/6s、5r/6s[12]和8r/4s [18],，然而最小问题目前尚无法通过代数方法解决。

在本文中，我们填补了最小问题方面的空白，提出了一种新方法，能够解决此前未解决的两种配置： 9r/3s和 6r/4s，如表1所示。与以往先求解偏移量再处理剩余到达时间（TOA）问题的到达时间差（ TDOA）方法不同，我们的方法结合了TOA和TDOA的思想，联合求解偏移量和位置，从而降低了整个流程的计算复杂度。此外，我们还对若干未解决问题的配置提供了计算复杂度的定量估计。

本文结构如下：在第二节中，回顾了到达时间（ TOA）和到达时间差（TDOA）求解技术的当前研究现状。在第三节中，阐述了提出的方法，并在第四节中对我们的求解器进行了基准测试。最后，在第五节中得出了结论。

2. 相关工作

在这些章节中，我们回顾了迄今为止用于代数求解到达时间和到达时间差的最先进的技术。

2.1. 到达时间

根据距离 di j ，我们可以定义压缩矩阵[11] ˜D ∈ R (m−1)×(n−1)，使得
[ ˜ D]ij = d 2 i+1, j +1 − d 2 1, j +1 − d 2 i+1, j + d 2 11 . (3)
通过代数运算可以证明，以下分解成立
D˜= −2R T S, (4)
其中Ri =[ri+1 −r 1]对于 i= 1… m − 1，类似地S j =[sj +1 − s 1]对于 j= 1… n− 1。

这种分解不是唯一的，假设我们有a fac‐˜ ˜ ˜torization D= RS。显然，对于每个满秩矩阵L，都有 ˜D= ˜R−1L˜L = ˜S T ˜R ˜S。在通过例如奇异值分解计算出 ˜R 和 S 之后，问题就简化为确定矩阵L，使得R −T ˜L˜ R且 S LS。现在可以如下对接收器和发射器进行参数化
r1= 0 s1= Lb, ri= L−T R˜i−1, i= 2… m,
sj= L(−1 2 S˜j−1+ b), j= 2… n,
(5)
其中b是一个待确定的向量。最后，定义H=(LTL)−1，可推导出以下方程
(A) d211= bT H−1b, (B) d21j − d211= 1 4 S˜jT−1H −1 S˜j−1 − bT H−1 S˜j−1, (C) d2i1 − d211= R˜T i−1H R˜i−1 − 2bT R˜i−1, (6)
其中 i= 2… m和 j= 2… n。由于矩阵H是对称的，因此可以用3个未知数进行参数化，因此H和b总共依赖于5个参数。给定 m个接收机和 n个发射机，可得到 1个3次A类方程、 n − 1个2次B类方程以及 m − 1个1次 C类方程。最终可通过这些方程利用Gr¨obner基[19, 17, 20]推导出快速多项式求解器。该方法已在[11]中用于求解到达时间最小问题。

2.2. 到达时间差

对于到达时间差，D取决于偏移量，因此无法进行数值分解。因此，前一节的参数化方法无法使用。然而，从R和S的维度可以看出，rank D= 2。如果 m> 3且 n> 3，这意味着D是秩亏的，因此所有 3 × 3子行列式必须等于零。尽管通常情况下有(m −1 3) ·(n −1 3)个子行列式，但以下定理成立。

定理1 给定一个秩为2的矩阵 A ∈ R(m×n), m, n> 2，则可获得 (m − 2)(n − 2)个独立的秩约束。[12]

由于压缩矩阵为 ˜D ∈ R(m−1)×(n−1)，给定m个接收器和 n个发射器，可获得(m − 3)(n − 3)个独立约束。该方法曾用于求解[12],中的某些到达时间差配置，其中利用秩约束确定所有偏移量，从而将问题简化为到达时间问题。但需要注意的是，仅靠秩约束无法求解到达时间差的最小情况。对于最小问题6r/3s，压缩矩阵已具有两列，因此无法推导出任何秩约束。对于另一种最小情况4接收器/5发射器，

只能获得两个独立的约束，这不足以求解五个未知的偏移量。本文中，我们提出了一种不同的压缩矩阵符号分解方法，该方法能够利用参数化表示(5)-(6)，从而为之前未解决的情况提供新的求解器。

3. 提出的方法

在本节中，我们描述了解决到达时间差（TDOA）问题所采用的数值技术。核心思想是利用式（4）中的分解，生成同时依赖于偏移量和坐标的新的方程。与以往的 TDOA方法不同，我们的目标是一步求解所有未知量。

首先，我们展示当仅有三个发射器存在时，如何获得一个简单的分解。接着，我们说明如何将其推广到更多发射器的情形。在我们所有的公式中，通过施加r1= 0和r2=[r2x, 0]T来固定规范自由。

3.1. 三个发射器

如果只有三个发射器可用，则无法施加秩约束。然而，可以注意到以下情况成立
D˜=(D ˜T)T I, (7)
其中I是单位矩阵。因此，我们可以通过施加 ˜R= ˜ DT和 ˜S= I来建立式(6)中的方程，得到m+ 2个关于8个未知数的方程（其中5个来自H和b，3个来自偏移量）。如前所述，最小情况为6r/3s。需要注意的是，现在C类方程尽管在H和b上线性，但由于R依赖于未知的偏移量，整体上仍为3次方程。对于次最小情况(m> 6)，我们的约束多于未知数。这就提出了一个问题：应如何从 m+2可用的方程中选取八个？首先，引入记号 abc表示使用 a个A类方程、 b个B类方程和 c个C类方程的组合形式。接下来，对每种组合形式，我们计算其对应的多项式理想的单项式基 B[21]。先前的研究 [19, 17],表明，单项式基的大小可作为该问题的复杂度度量：其值越小，最终求解器的速度越快且精度越高。仿真的结果如表2所示。

从表2可以看出，C类方程导致最低的计算复杂度，而A类方程导致最高的计算复杂度。特别是对于9r/3s情况，该问题可以仅使用C类方程求解。表中的最后一行对应于最小问题6r/3s。观察表1，这个新求解器使得仅用3个发射机即可校准网络，而之前的最先进的求解器无法处理仅有3个发射机的网络。

总结前面的讨论，我们得到以下新 9r/3s求解器的最终数值步骤：
1. 使用因式分解（7）和来自（6）的C类方程建立多项式方程。
2. 利用[19]和[17],中描述的技术，为该多项式系统生成求解器。
3. 使用（5）计算接收器和发射机的位置。
4. 使用非线性优化（例如 Levenberg‐Marquardt算法）优化上一步的解。

3.2. 多于三个发射器

如果我们有超过三个发射器，则S将不是方阵，因此无法直接替换为单位矩阵。因此，为了使用(7)中的分解，我们需要舍弃一些发射器。特别是，让我们考虑之前未解决的情况 6r/4s，利用(6)中的参数化方法，我们总共得到9个未知数，其中5个来自H和b，另外4个来自偏移量。现在我们的压缩矩阵D的尺寸为 5 × 3。为了解决该构型，我们考虑新的压缩矩阵ˆ˜D，即舍弃D的最后一列后得到的矩阵。采用前一节中描述的相同方法，我们得到1个A类方程、2个B类方程和5个C类方程，共计8个方程。此外，由于原始压缩矩阵D必须为秩2，通过令其所有子式为0，我们可得到关于偏移量的10个方程，根据定理1，其中仅有3个是独立的。我们再次针对上述方程的不同组合计算了标准单项式基，并确定最高效率的求解器应包含全部5个C类方程、1个B类方程以及全部10个秩约束。该形式化方法导致标准单项式基的大小为 |B| = 22。尽管仅需3个秩约束就已足够，但数值实验表明，使用全部10个方程能得到更稳定的求解器。特别是，如果仅使用3个秩方程，则会得到|B| = 66。一旦偏移量以及参数H和b被确定，接收器和前3个发射器便可由(5)确定。最后一个发射器s 4 最终通过三边测量法计算得出。

4. 结果

本节介绍了用于定量评估我们求解器的实验。这些求解器在合成数据上进行了基准测试。接收器和发射器的位置是从标准差为10的零均值正态分布中采样的。时间偏移量由标准正态分布生成。我们生成的Matlab求解器能够得到精确的解，并在Intel i7‐8565U处理器上以 ∼ 200 ms的速度运行，因此适用于近实时应用。

除了提出的两种求解器外，我们还使用本文描述的方法计算了介于先前已解决和最小配置之间的标准单项式基。表3所示的结果可用于粗略评估未解决情况的计算复杂度，并对其可行性提供提示。

4.1. 干净数据

为了评估我们的求解器，我们随机生成输入数据，并利用这些数据求解自校准问题。在此步骤中，我们考虑无噪声数据。我们还将我们的求解器与迭代算法在随机位置和偏移量初始化下的表现进行比较，这些随机值来自与真实值相同的分布。通过1000次不同随机数据运行求解器所得到的相对误差分布如图1所示，中位相对误差见表4。可以看出，即使对于干净数据，随机初始化也无法收敛。而我们的求解器则能产生精确的解，其相对误差实际上等于计算机的机器精度。此外，以我们的求解器结果作为起点，迭代算法在极少的迭代次数后即可收敛。

4.2. 噪声数据

我们还研究了我们的求解器在噪声数据下的表现，通过向测量值 fij 添加均值为零、标准差为 σ的高斯噪声。仿真的结果如图2所示。正如之前讨论的，到达时间差问题对初始值的选择很敏感，较差的初始化可能导致收敛问题。虽然随机初始化会失败，但使用我们的求解器可以消除这一问题，即使在较高的噪声水平下也能实现快速且准确的收敛。

5. 结论

本文中，我们研究了基于二维TDOA测量的传感器网络自校准问题。我们提出了一种新算法，由此开发出新的鲁棒且高效的多项式求解器，能够解决此前未解决的两种配置。我们证明了所提出方法在干净和含噪声数据下均具有稳定的解。