背景
我已经工作了一段时间,尝试解决3-dimensions
和使用4
节点中的(众所周知的痛苦)到达时间差(TDoA)多面问题。如果您不熟悉此问题,则应确定给定(X,Y,Z)
节点的坐标,信号到达每个节点的时间以及信号n
的速度,从而确定某些信号源v
的坐标。
我的解决方案如下:
对于每个节点,我们编写(X-x_i)**2 + (Y-y_i)**2 + (Z-z_i)**2 = (v(t_i - T)**2
其中(x_i, y_i, z_i)
是ith
节点的坐标,而T
是发射时间。
现在,我们在4
未知数中具有4
方程。四个节点显然不足。我们可以尝试直接解决该系统,但是考虑到问题的高度非线性性质,这似乎几乎是不可能的(而且,实际上,我已经尝试了许多直接技术...并且失败了)。相反,我们通过考虑所有i/j
可能性(从等式i
中减去等式j
),将其简化为线性问题。我们获得以下形式的(n(n-1))/2
=6
方程:2*(x_j - x_i)*X + 2*(y_j - y_i)*Y + 2*(z_j - z_i)*Z + 2 * v**2 * (t_i - t_j) = v**2 ( t_i**2 - t_j**2) + (x_j**2 + y_j**2 + z_j**2) - (x_i**2 + y_i**2 + z_i**2)
看起来像Xv_1 + Y_v2 + Z_v3 + T_v4 = b
。现在,我们尝试应用标准线性最小二乘法,其中的解决方案是x
中的矩阵向量A^T Ax = A^T b
。不幸的是,如果您尝试将其输入任何标准的线性最小二乘算法中,它将变得很困难。那我们现在怎么办?
...
信号到达节点i
的时间(当然)由以下公式给出:sqrt( (X-x_i)**2 + (Y-y_i)**2 + (Z-z_i)**2 ) / v
该方程式表明到达时间T
是0
。如果我们有那个T = 0
,我们可以将T
列放在矩阵A
中,这样就大大简化了问题。的确,NumPy's
linalg.lstsq()
提供了令人惊讶的准确结果。
...
因此,我要做的是通过从每个方程式中减去最早的时间来归一化输入时间。然后,我要做的就是确定可以添加到每次的dt
,以使线性最小二乘找到的点的平方误差总和最小化。
我将某些dt
的误差定义为通过将input times
+ dt
输入最小二乘算法预测的点的到达时间之间的平方差,减去输入时间(归一化),并加总所有4
节点。
for node, time in nodes, times:
error += ( (sqrt( (X-x_i)**2 + (Y-y_i)**2 + (Z-z_i)**2 ) / v) - time) ** 2
我的问题:通过使用蛮力,我能够令人满意地做到这一点。我从
dt = 0
开始,然后向上移动了一些步骤,直到达到最大的迭代次数,或者直到达到最小的RSS错误为止,这就是我添加到归一化时间中的dt
以获得解决方案。最终的解决方案非常精确,但速度很慢。实际上,我希望能够实时解决此问题,因此将需要一种速度更快的 解决方案。我首先假设误差函数(即上面定义的
dt
与error
)将是高度非线性的-即兴,这对我来说很有意义。由于我没有实际的数学函数,因此可以自动排除需要微分的方法(例如
Newton-Raphson
)。错误函数将始终为正,因此我可以排除bisection
等。相反,我尝试简单的近似搜索。不幸的是,那失败了。然后,我尝试了禁忌搜索,然后尝试了遗传算法,并进行了其他几种搜索。他们都惨败。因此,我决定进行一些调查。事实证明,误差函数与dt的关系图看起来有点像平方根,仅根据信号源与节点的距离来右移:
其中
dt
在水平轴上,错误在垂直轴上而且,事后看来,当然可以!我将误差函数定义为包含平方根,因此,至少在我看来,这是合理的。
怎么办?
因此,我现在的问题是,如何确定与误差函数最小值对应的
dt
?我的第一次尝试(非常粗略)是在误差图上获得一些点(如上所述),使用
numpy.polyfit
对其进行拟合,然后将结果提供给numpy.root
。该根对应于dt
。不幸的是,这也失败了。我尝试使用各种degrees
以及各种要点,直至达到荒谬的点数,这样我也可以只使用蛮力。如何确定与此错误函数的最小值相对应的
dt
? 由于我们处理的是高速信号( radio 信号),因此结果的准确性和准确性非常重要,因为
dt
中的微小差异会偏离结果点。我确定我在这里所做的工作中埋藏着一些无限简单的方法,但是忽略了其他所有内容,如何找到
dt
?我的要求:
Python
和NumPy
编辑:
这是我的代码。诚然,有点困惑。在这里,我正在使用
polyfit
技术。它将为您“模拟”一个源,并比较结果:from numpy import poly1d, linspace, set_printoptions, array, linalg, triu_indices, roots, polyfit
from dataclasses import dataclass
from random import randrange
import math
@dataclass
class Vertexer:
receivers: list
# Defaults
c = 299792
# Receivers:
# [x_1, y_1, z_1]
# [x_2, y_2, z_2]
# [x_3, y_3, z_3]
# Solved:
# [x, y, z]
def error(self, dt, times):
solved = self.linear([time + dt for time in times])
error = 0
for time, receiver in zip(times, self.receivers):
error += ((math.sqrt( (solved[0] - receiver[0])**2 +
(solved[1] - receiver[1])**2 +
(solved[2] - receiver[2])**2 ) / c ) - time)**2
return error
def linear(self, times):
X = array(self.receivers)
t = array(times)
x, y, z = X.T
i, j = triu_indices(len(x), 1)
A = 2 * (X[i] - X[j])
b = self.c**2 * (t[j]**2 - t[i]**2) + (X[i]**2).sum(1) - (X[j]**2).sum(1)
solved, residuals, rank, s = linalg.lstsq(A, b, rcond=None)
return(solved)
def find(self, times):
# Normalize times
times = [time - min(times) for time in times]
# Fit the error function
y = []
x = []
dt = 1E-10
for i in range(50000):
x.append(self.error(dt * i, times))
y.append(dt * i)
p = polyfit(array(x), array(y), 2)
r = roots(p)
return(self.linear([time + r for time in times]))
# SIMPLE CODE FOR SIMULATING A SIGNAL
# Pick nodes to be at random locations
x_1 = randrange(10); y_1 = randrange(10); z_1 = randrange(10)
x_2 = randrange(10); y_2 = randrange(10); z_2 = randrange(10)
x_3 = randrange(10); y_3 = randrange(10); z_3 = randrange(10)
x_4 = randrange(10); y_4 = randrange(10); z_4 = randrange(10)
# Pick source to be at random location
x = randrange(1000); y = randrange(1000); z = randrange(1000)
# Set velocity
c = 299792 # km/ns
# Generate simulated source
t_1 = math.sqrt( (x - x_1)**2 + (y - y_1)**2 + (z - z_1)**2 ) / c
t_2 = math.sqrt( (x - x_2)**2 + (y - y_2)**2 + (z - z_2)**2 ) / c
t_3 = math.sqrt( (x - x_3)**2 + (y - y_3)**2 + (z - z_3)**2 ) / c
t_4 = math.sqrt( (x - x_4)**2 + (y - y_4)**2 + (z - z_4)**2 ) / c
print('Actual:', x, y, z)
myVertexer = Vertexer([[x_1, y_1, z_1],[x_2, y_2, z_2],[x_3, y_3, z_3],[x_4, y_4, z_4]])
solution = myVertexer.find([t_1, t_2, t_3, t_4])
print(solution)
最佳答案
似乎Bancroft方法适用于此问题?这是一个纯NumPy实现。
# Implementation of the Bancroft method, following
# https://gssc.esa.int/navipedia/index.php/Bancroft_Method
M = np.diag([1, 1, 1, -1])
def lorentz_inner(v, w):
return np.sum(v * (w @ M), axis=-1)
B = np.array(
[
[x_1, y_1, z_1, c * t_1],
[x_2, y_2, z_2, c * t_2],
[x_3, y_3, z_3, c * t_3],
[x_4, y_4, z_4, c * t_4],
]
)
one = np.ones(4)
a = 0.5 * lorentz_inner(B, B)
B_inv_one = np.linalg.solve(B, one)
B_inv_a = np.linalg.solve(B, a)
for Lambda in np.roots(
[
lorentz_inner(B_inv_one, B_inv_one),
2 * (lorentz_inner(B_inv_one, B_inv_a) - 1),
lorentz_inner(B_inv_a, B_inv_a),
]
):
x, y, z, c_t = M @ np.linalg.solve(B, Lambda * one + a)
print("Candidate:", x, y, z, c_t / c)
关于python - 使用NumPy最小化此错误功能,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/65436017/