给定三个接收器的位置以及它们接收信号的时间(到达时间延迟),如何定位信号?

How to localize a signal given the location of three receivers and the times at which when they receive the signal (Time Delay of Arrival)?

我有 3 个接收器(A、B 和 C)和一些位置未知的信号产生源(比如声音或光)。给定 A、B、C 的位置,以及每个接收器“听到”信号的时间,我想确定信号源的方向。

我知道可以通过 TDoA multilateration/trilateration 实现这一点,但是我在执行计算时遇到了问题。对于那些完全不熟悉该主题的人来说,没有很多关于此的清晰,详细的信息。外面的东西对我来说是模糊的,更理论化的,或者有点太深奥了。

关于 SO 的一些类似帖子(但不是我想要的): TDOA multilateration to locate a sound source Trilateration of a signal using Time Difference(TDOA)

这也很有趣,但假设我们有一些界限: Multiliteration implementation with inaccurate distance data

@Dave 还评论了一个优秀且相当容易访问的资源 https://sites.tufts.edu/eeseniordesignhandbook/files/2017/05/FireBrick_OKeefe_F1.pdf,但它不够深入,以至于人们可能能够在代码中实际实现它(至少,对于没有深入了解的人来说)回归,找到所得双曲线的交点等)。

[编辑]:我应该补充一点,我可以假设 3 个传感器 源在地球表面,地球曲率的影响是可以忽略不计(即我们可以在二维中工作)。

最简单(但不是最快)的方法是用 gradient descent 求解方程。

我假设我们知道

  • 接收器A、B、C不在同一直线上的位置;
  • 未知源X到A、B、C的伪距

直观上,我们模拟了一个物理系统,其中三个ideal springs配置成这样,其中每个spring的平衡长度是相应的伪距。

  A
  |
  X
 / \
B   C

spring 距离太小时推,距离太大时拉。 X 的近似静止位置应该是一个合理的估计(尽管根据您的应用程序,您可能需要进行额外的验证)。

下面是一些示例 Python 代码,使用复数作为 2D 向量,应该很容易音译。

import random


def distance(p, q):
    return abs(p - q)


# Force exerted by an ideal spring between variable y and fixed q of equilibrium
# length dxq.
def force(y, q, dxq):
    dyq = distance(y, q)
    return (dxq - dyq) * (y - q) / dyq


# Trilateration via gradient descent.
def trilaterate(
    a, dxa, b, dxb, c, dxc, *, max_iterations=1000000, gamma=0.001, precision=1e-12
):
    # Use the centroid of the receivers as the initial estimate.
    y = (a + b + c) / 3
    for i in range(max_iterations):
        f = force(y, a, dxa) + force(y, b, dxb) + force(y, c, dxc)
        y += gamma * f
        if abs(f) <= precision:
            return y
    return None


def random_point():
    return complex(random.random(), random.random())


def test_error():
    a = random_point()
    b = random_point()
    c = random_point()
    x = random_point()
    y = trilaterate(a, distance(x, a), b, distance(x, b), c, distance(x, c))
    return distance(x, y)


if __name__ == "__main__":
    print(test_error())

有趣的问题。我懒得推导出代数解的方程。相反,为什么不适合结果?

因此,使用任何能够找到局部解决方案的拟合方法(使用某些误差值的 optimization/minimization)简单地拟合 2D(或更高)位置。当我使用我的简单 来适应位置时,结果看起来很不错。

算法是:

  1. 遍历范围内的“所有”位置

    of coarse 并不是所有的试探法都会大大减少问题。

  2. 在每个测试位置上计算将要测量的增量时间

    从测试位置到接收站的简单行进时间。

  3. 标准化所有增量时间,以便从零开始

    所以从所有接收者时间中减去最小的到达时间。实际测量时间也是如此。这确保时间不涉及相对偏移。

  4. 计算实际测量时间和计算时间之间的差异

    简单的腹肌差异就足够了。使用此值作为拟合参数(优化)。

这里使用我上面 link 中的近似值 class 执行此操作的小 C++ 示例:

//---------------------------------------------------------------------------
// TDoA Time Difference of Arrival
//---------------------------------------------------------------------------
const int n=3;
double recv[n][3];  // (x,y) [m] receiver position,[s] time of arrival of signal
double pos0[2];     // (x,y) [m] object's real position
double pos [2];     // (x,y) [m] object's estimated position
double v=340.0;     // [m/s] speed of signal
double err=0.0;     // [m] error between real and estimated position
//---------------------------------------------------------------------------
void compute()
    {
    int i;
    double x,y,a,da,t0;
    //---------------------------------------------------------
    // init positions
    da=2.0*M_PI/(n);
    for (a=0.0,i=0;i<n;i++,a+=da)
        {
        recv[i][0]=256.0+(220.0*cos(a));
        recv[i][1]=256.0+(220.0*sin(a));
        }
    pos0[0]=300.0;
    pos0[1]=220.0;
    // simulate measurement
    t0=123.5;                   // some start time
    for (i=0;i<n;i++)
        {
        x=recv[i][0]-pos0[0];
        y=recv[i][1]-pos0[1];
        a=sqrt((x*x)+(y*y));    // distance to receiver
        recv[i][2]=t0+(a/v);    // start time + time of travel
        }
    //---------------------------------------------------------
    // normalize times into deltas from zero
    a=recv[0][2]; for (i=1;i<n;i++) if (a>recv[i][2]) a=recv[i][2];
    for (i=0;i<n;i++) recv[i][2]-=a;
    // fit position
    int N=6;
    approx ax,ay;
    double e,dt[n];
              // min,  max,step,recursions,&error
    for (ax.init( 0.0,512.0, 32.0        ,N,   &e);!ax.done;ax.step())
     for (ay.init(  0.0,512.0, 32.0       ,N,   &e);!ay.done;ay.step())
        {
        // simulate measurement -> dt[]
        for (i=0;i<n;i++)
            {
            x=recv[i][0]-ax.a;
            y=recv[i][1]-ay.a;
            a=sqrt((x*x)+(y*y));    // distance to receiver
            dt[i]=a/v;              // time of travel
            }
        // normalize times dt[] into deltas from zero
        a=dt[0]; for (i=1;i<n;i++) if (a>dt[i]) a=dt[i];
        for (i=0;i<n;i++) dt[i]-=a;
        // error
        e=0.0; for (i=0;i<n;i++) e+=fabs(recv[i][2]-dt[i]);
        }
    pos[0]=ax.aa;
    pos[1]=ay.aa;
    //---------------------------------------------------------
    // compute error
    x=pos[0]-pos0[0];
    y=pos[1]-pos0[1];
    err=sqrt((x*x)+(y*y));  // [m]
    }
//---------------------------------------------------------------------------

此处预览:

蓝点是接收器,红点是物体的真实位置,黄色十字是它的估计位置。该区域是 512x512 m,我用初始步骤 32 m6 递归拟合它,导致错误 ~0.005 m

我对结果非常满意...您可以更改接收器的数量 n 而无需对源或算法进行任何实际更改。我将接收器位置启动为均匀分布在圆上,但位置可能是任何其他位置(并非全部在粗略的单行上)