这篇文章发表在 André Chr. Andersen 的博客上,http://blog.andersen.im/2012/07/signal-emitter-positioning-using-multilateration
很遗憾,他的博客现在打不开了,可能他已经没有维护了。
我们第一个版本的定位引擎使用的TDOA算法,就来自他的这篇文章。André Chr. Andersen在文章后面还提供了 MATLAB 代码,Paul Hayes把Andersen的MATLAB代码翻译为 Python 代码,可以在 github 上找到 https://github.com/paulhayes/MultilaterationExample
这个算法是针对超声波定位的,但是本质上跟 UWB 定位是一样的,重要的是这是一个简单易实现的 Multilateration 算法,改一改就可以用在 uwb 的坐标计算上。
非常感谢 André Chr. Andersen 和 Paul Hayes 两位,如果没有他们的算法和代码,也许我们这个项目就提前中止了。
Andersen 的算法可以用,但是也有一些问题。主要体现在两个方面:
(1)、计算出来的坐标不太准。标签在区域中心的时候比较准,但是当标准在区域边缘时,误差会变得很大。
(2)、算法中迭代的过程中可能会出错,导致有时算不出坐标。因为各种因素影响,我们得到的TDOA值肯定是不准确的、有误差的,实际上我们期望得到的是一个近似值。但是误差太大,有时会得出很离谱的坐标,这个就比较尴尬了。
我们后来自己研发了两个算法,最终采用的是最小二乘法,在速度和精度方面得到较好的平衡。
以下是 Paul Hayes 重写的 Python 代码。
###########################
#
# Python rewrite of multilateration technique by André Andersen in his [blog post](http://blog.andersen.im/2012/07/signal-emitter-positioning-using-multilateration).
#
############################

from numpy import *
from numpy.linalg import *
import json
#speed of sound in medium
v = 3450
numOfDimensions = 3
nSensors = 5
region = 3
sensorRegion = 2

#choose a random sensor location
emitterLocation = region * ( random.random_sample(numOfDimensions) - 0.5 )
sensorLocations = [ sensorRegion * ( random.random_sample(numOfDimensions)-0.5 ) for n in range(nSensors) ]
p = matrix( sensorLocations ).T

#Time from emitter to each sensor
sensorTimes = [ sqrt( dot(location-emitterLocation,location-emitterLocation) ) / v for location in sensorLocations ]

c = argmin(sensorTimes)
cTime = sensorTimes[c]

#sensors delta time relative to sensor c
t = sensorDeltaTimes = [ sensorTime - cTime for sensorTime in sensorTimes ]

ijs = range(nSensors)
del ijs[c]

A = zeros([nSensors-1,numOfDimensions])
b = zeros([nSensors-1,1])
iRow = 0
rankA = 0

for i in ijs:
	for j in ijs:
		A[iRow,:] = 2*( v*(t[j])*(p[:,i]-p[:,c]).T - v*(t[i])*(p[:,j]-p[:,c]).T )
		b[iRow] = v*(t[i])*(v*v*(t[j])**2-p[:,j].T*p[:,j]) + \
		(v*(t[i])-v*(t[j]))*p[:,c].T*p[:,c] + \
		v*(t[j])*(p[:,i].T*p[:,i]-v*v*(t[i])**2)
		rankA = matrix_rank(A)
		if rankA >= numOfDimensions :
			break
		iRow += 1
	if rankA >= numOfDimensions:
		break

calculatedLocation = asarray( lstsq(A,b)[0] )[:,0]

print "Emitter location: %s " % emitterLocation
print "Calculated position of emitter: %s " % calculatedLocation

版权所有

本站文章均为原创,版权属 Zhang Xiaolong 所有。

所有人可以转载,但必须注明出处以及作者信息。