使用Numpy linalg.lstsq求解线性系统时获得(明显)不准确的值

Question

我试图使用Numpy的linalg.lstsq()来遵循本文中的逻辑。我在这里重复一遍。

---

假设我在三维空间中有四个麦克风(或光接收器等(。我想用每个麦克风记录枪击事件之间的时间差，或者用到达每个麦克风的时间来确定枪手的坐标。因此，我们知道了四个麦克风的坐标x_i, y_i和声音的速度v。我们不知道枪手的坐标，也不知道开枪的时间。

对于每个麦克风，我们写：

(X-x_i)**2 + (Y-y_i)**2 + (Z-z_i)**2 = (v(t_i - T))**2

其中，X,Y,Z是射击手的坐标，T是射击的时间。

如果我们考虑所有i/j的可能性，从方程j中减去方程i，这个问题就可以简化。我们最终得到了形式为的6方程(通常，给定n麦克风，n(n-1)/2方程(

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) = 2 * 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，其中v_i是矢量。然后我们(应该(能够使用线性最小二乘拟合来"；解决"；(或者，至少可以粗略估计(X,Y,Z和作为副产品的T。

---

我已经尝试了很多不同的方法来解决这个问题，所有迹象都指向使用链接文章中描述的方法(即，通过线性最小二乘法解决问题，然后使用这个"解决方案"作为更"精确"算法的初始估计(。因此，我编写了以下(诚然，有些粗糙(Python代码来开始。它解决了问题的线性最小二乘部分：

from dataclasses import dataclass
from numpy import *
from random import randrange
@dataclass
class Vertexer:
@staticmethod
def find():
# Pick microphones 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 shooter to be at random location
x = randrange(100); y = randrange(100); z = randrange(100)
# Set velocity (ok, it's a ray gun...)
c = 299792 # km/s
# 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

A = array([[2 * (x_1 - x_2), 2 * (y_1 - y_2), 2 * (z_1 - z_2), 2 * c * c * (t_2 - t_1)],
[2 * (x_1 - x_3), 2 * (y_1 - y_3), 2 * (z_1 - z_3), 2 * c * c * (t_3 - t_1)],
[2 * (x_1 - x_4), 2 * (y_1 - y_4), 2 * (z_1 - z_4), 2 * c * c * (t_4 - t_1)],
[2 * (x_2 - x_3), 2 * (y_2 - y_3), 2 * (z_2 - z_3), 2 * c * c * (t_3 - t_2)],
[2 * (x_2 - x_4), 2 * (y_2 - y_4), 2 * (z_2 - z_4), 2 * c * c * (t_4 - t_2)],
[2 * (x_3 - x_4), 2 * (y_3 - y_4), 2 * (z_3 - z_4), 2 * c * c * (t_4 - t_3)]
])
b = array([2 * c * c * (t_2 * t_2 - t_1 * t_1) + (x_1 * x_1 + y_1 * y_1 + z_1 * z_1) - (x_2 * x_2 + y_2 * y_2 + z_2 * z_2),
2 * c * c * (t_3 * t_3 - t_1 * t_1) + (x_1 * x_1 + y_1 * y_1 + z_1 * z_1) - (x_3 * x_3 + y_3 * y_3 + z_3 * z_3),
2 * c * c * (t_4 * t_4 - t_1 * t_1) + (x_1 * x_1 + y_1 * y_1 + z_1 * z_1) - (x_4 * x_4 + y_4 * y_4 + z_4 * z_4),
2 * c * c * (t_3 * t_3 - t_2 * t_2) + (x_2 * x_2 + y_2 * y_2 + z_2 * z_2) - (x_3 * x_3 + y_3 * y_3 + z_3 * z_3),
2 * c * c * (t_4 * t_4 - t_2 * t_2) + (x_2 * x_2 + y_2 * y_2 + z_2 * z_2) - (x_4 * x_4 + y_4 * y_4 + z_4 * z_4),
2 * c * c * (t_4 * t_4 - t_3 * t_3) + (x_3 * x_3 + y_3 * y_3 + z_3 * z_3) - (x_4 * x_4 + y_4 * y_4 + z_4 * z_4)])
solved, residuals, rank, s = linalg.lstsq(A, b)
print(solved)

myVertexer = Vertexer()
myVertexer.find()

然而，不幸的是，它似乎没有产生非常合理的估计。如果你愿意，你可以自己运行代码，或者这里有一组运行，将找到的值与"的坐标进行比较；模拟的"；算法产生的源。

Predicted: [ 2.26739519e+00  4.80191502e+00 -5.07181020e+00 ]
Actual:    78 35 57
Predicted: [ 7.72173135e-01 -2.61250803e+00  4.10306750e+00 ]
Actual:    30 75 48
Predicted: [-1.65368110e+00  8.82919123e-01  1.43648336e+00 ]
Actual:    67 56 44
Predicted: [ 4.08698715e+00 -3.19377148e-01  8.81706364e-01 ]
Actual:    9 82 13

诚然，我不确定初始最小二乘解的精度是多少。然而，这些价值观对我来说似乎不太好。我哪里错了？

请让我知道我是否可以提供任何进一步的澄清。

请注意，我知道还有其他方法可以解决这个问题(值得注意的是，只需插入带有非线性求解器的多个软件包中的一个，就可以获得相当好的解决方案(。我感兴趣的是确定为什么这个解决方案(我知道链接的海报已经成功地使用了它，还有其他几个(不起作用，这既是出于好奇，也是因为我想更进一步，其中"；自定义"；代码将非常有用。

Answer 1

有两个独立的问题：

1. 您在原始帖子中提到的减法方程中有一个错误：在RHS上，它应该是c * c * (t_i * t_i - t_j * t_j)而不是2 * c * c * (t_i * t_i - t_j * t_j)
2. 第二个问题是，在矩阵A中，您有一列tau(第四列(，而在计算t_1到t_4时，您隐含地假设tau=0，因此您不能将其作为矩阵A中的变量；因此需要移除第四列

修订后的代码如下：

from random import randrange, seed
from numpy import array, linalg, set_printoptions
from dataclasses import dataclass
set_printoptions(suppress=True, linewidth=1000, precision=6)
seed(5)
@dataclass
class Vertexer:
@staticmethod
def find():
# Pick microphones 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 shooter to be at random location
x = randrange(100); y = randrange(100); z = randrange(100)
# Set velocity (ok, it's a ray gun...)
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
A = array([[2 * (x_1 - x_2), 2 * (y_1 - y_2), 2 * (z_1 - z_2)],
[2 * (x_1 - x_3), 2 * (y_1 - y_3), 2 * (z_1 - z_3)],
[2 * (x_1 - x_4), 2 * (y_1 - y_4), 2 * (z_1 - z_4)],
[2 * (x_2 - x_3), 2 * (y_2 - y_3), 2 * (z_2 - z_3)],
[2 * (x_2 - x_4), 2 * (y_2 - y_4), 2 * (z_2 - z_4)],
[2 * (x_3 - x_4), 2 * (y_3 - y_4), 2 * (z_3 - z_4)]
])
b = array([1 * c * c * (t_2 * t_2 - t_1 * t_1) + (x_1 * x_1 + y_1 * y_1 + z_1 * z_1) - (x_2 * x_2 + y_2 * y_2 + z_2 * z_2),
1 * c * c * (t_3 * t_3 - t_1 * t_1) + (x_1 * x_1 + y_1 * y_1 + z_1 * z_1) - (x_3 * x_3 + y_3 * y_3 + z_3 * z_3),
1 * c * c * (t_4 * t_4 - t_1 * t_1) + (x_1 * x_1 + y_1 * y_1 + z_1 * z_1) - (x_4 * x_4 + y_4 * y_4 + z_4 * z_4),
1 * c * c * (t_3 * t_3 - t_2 * t_2) + (x_2 * x_2 + y_2 * y_2 + z_2 * z_2) - (x_3 * x_3 + y_3 * y_3 + z_3 * z_3),
1 * c * c * (t_4 * t_4 - t_2 * t_2) + (x_2 * x_2 + y_2 * y_2 + z_2 * z_2) - (x_4 * x_4 + y_4 * y_4 + z_4 * z_4),
1 * c * c * (t_4 * t_4 - t_3 * t_3) + (x_3 * x_3 + y_3 * y_3 + z_3 * z_3) - (x_4 * x_4 + y_4 * y_4 + z_4 * z_4)])
solved, residuals, rank, s = linalg.lstsq(A, b, rcond=None)
print('Actual:', [x, y, z])
print('Predicted:', solved)
myVertexer = Vertexer()
myVertexer.find()

输出为：

Actual: [31, 48, 69]
Predicted: [31. 48. 69.]