1

I am having some trouble to find or implement an algorithm to find a signal source. The objective of my work is to find the sound emitter position.

To accomplish this I am using three vibration sensors. The technique that I am using is multilateration that is based on the time difference of arrival.

The time difference of arrival between each sensor are found using Cross Correlation of the received signals.

I already implemented the algorithm to find the time difference of arrival, but my problem is more on how multilateration works, it's unclear for me based on my reference, and I couldn't find any other good reference for this that are free/open.

I can't figure out how to do it having only the time differences of arrival ( Tab,Tac,Tbc)

Any help on this would be much appreciated

2 Answers2

2

Below is a closed-form linearized solution. It looks for intersection of hyperboloids and makes approximation.

For 3D localisation minimum of 4 anchors (receivers) are required. To solve equation, 0'th anchor is used as reference for linearization. It does not matter which anchor you select for linearization.

\begin{split} \begin{aligned} N & \text{ : Number of anchors} \\ Si(x_i,y_i,z_i) & \text{ : The coordinates of anchors (i=0,1,2,3 ... N-1)} \\ S(x,y,z) & \text{ : The coordinate of tag (transmitter)} \\ d_{i} & \text{ : Distance between i'th anchor and tag} \\ d_{0i} & \text{ : Distance difference between 0'th anchor and i'th anchor} \\ \tau _{0i} & \text{ : Time difference of arrival between i'th and 0'th anchors (i=1,2,3 ... N-1) } \\ c & \text{ : Speed of light/sound } \end{aligned} \end{split}

\begin{split} \tau _{0i} & = \frac{1}{c} |S-S_0| - \frac{1}{c} |S-S_i|\\ & = \frac{1}{c} (d_0 - d_i) \end{split}

\begin{split} d_{0i} = \tau _{0i}.c = d_{0} - d_{i} \end{split}

\begin{equation} \begin{split} \Rightarrow d_{i} = d_0 - d_{0i} \end{split} \end{equation}

\begin{split} d_{0}^2 - d_{i}^2 &= |S-S_0|^2 - |S-S_i|^2\\ &=(x-x_0)^2 + (y-y_0)^2 + (z-z_0)^2- (x-x_i)^2 - (y-y_i)^2 - (z-z_i)^2\\ &=x_0^2-x_i^2+2x(x_i-x_0)+y_0^2-y_i^2+2y(y_i-y_0) + z_0^2-z_i^2+2z(z_i-z_0)\\ \end{split}

Substitude $d_i$ for $d_0-d_{0i}$

\begin{equation} \begin{split} d_{0}^2 - (d_0 - d_{0i})^2 &= x_0^2-x_i^2 + 2x(x_i-x_0) + y_0^2-y_i^2+2y(y_i-y_0) + z_0^2-z_i^2+2z(z_i-z_0)\\ 2d_0.d_{0i} - d_{0i}^2 &= x_0^2-x_i^2 + 2x(x_i-x_0) + y_0^2-y_i^2+2y(y_i-y_0) + z_0^2-z_i^2+2z(z_i-z_0)\\ \end{split} \end{equation}

Regrouping the terms lead to

\begin{equation} \begin{split} x(x_0-x_i) + y(y_0-y_i) + z(z_0-z_i) + d_{0i}.d_0 &= \frac{1}{2}(x_0^2-x_i^2 + y_0^2-y_i^2 + z_0^2-z_i^2 + d_{0i}^2)\\ \end{split} \end{equation}

These linear system of equations are easily written in matrix form as below. (i = 1,2,3 .. N-1)

\begin{gather} \begin{bmatrix} x_0-x_1 & y_0-y_1 & z_0-z_1 & d_{01}\\ x_0-x_2 & y_0-y_2 & z_0-z_2 & d_{02}\\ x_0-x_3 & y_0-y_3 & z_0-z_3 & d_{03}\\ ... & ... & ... & ...\\ x_0-x_n & y_0-y_n & z_0-z_n & d_{0n}\\ \end{bmatrix} \begin{bmatrix} x\\ y\\ z\\ d_0\\ \end{bmatrix} = \begin{bmatrix} \frac{1}{2}(x_0^{2} - x_1^{2} + y_0^{2} - y_1^{2} + z_0^{2} - z_1^{2} + d_{01}^{2}) \\ \frac{1}{2}(x_0^{2} - x_2^{2} + y_0^{2} - y_2^{2} + z_0^{2} - z_2^{2} + d_{02}^{2}) \\ \frac{1}{2}(x_0^{2} - x_3^{2} + y_0^{2} - y_3^{2} + z_0^{2} - z_3^{2} + d_{03}^{2}) \\ ... \\ \frac{1}{2}(x_0^{2} - x_n^{2} + y_0^{2} - y_n^{2} + z_0^{2} - z_n^{2} + d_{0n}^{2}) \\ \end{bmatrix} \end{gather}

It's an ordinary $A\vec{x}=\vec{b}$ equation. Matrice A and vector $\vec{b}$ are known. The problem requires determination of $\vec{x}$ such that $A\vec{x} \approx \vec{b}$ minimizing sum of squares of residuals. It's called least square regression.

Assuming that error is $\varepsilon$

$\varepsilon(\vec{x}) = A\vec{x} - \vec{b}$

Minimize square of error function $\parallel A\vec{x} - \vec{b}\parallel^2$ by deriving it with respect to $\vec{x}$ and then equate zero.

\begin{equation} \begin{split} \frac{\partial}{\partial \vec{x}} \parallel A\vec{x}-\vec{b}\parallel^2 &=0\\ \frac{\partial}{\partial \vec{x}} [(A\vec{x}-\vec{b})^T(A\vec{x}-\vec{b})] &= 0\\ \frac{\partial}{\partial \vec{x}} [\vec{x}^T A^T A \vec{x} - 2\vec{x}^TA^T \vec{b} + \vec{b}^T \vec{b} ] &= 0\\ 2A^T A \vec{x} - 2A^T \vec{b} &= 0 \\ \Rightarrow \vec{x} &= (A^TA)^{-1} A^T\vec{b} \end{split} \end{equation}

The vector $\vec{x}$ is solution.

If anchors are not placed uniformly then matrice $A$ is not singular so it can be invertible and solution exists, otherwise QR decomposition can be used. It's another topic.

It's the simplest way for solution. At this point, mostly $\vec{x}$ is used as initial value of further complex iterations. There are lots of different approach for fine localization.

KetZoomer
  • 115
inninaro
  • 121
  • 6
0

If you look here, you will find something very similar to inninaro's answer.

However, when this preliminary step has been done (which is just a multilinear regression), I strongly suggest you use nonlinear regression and fit the real model $$t=T+\frac 1 c\sqrt{(x-X)^2+(y-Y)^2+(z-Z)^2}$$ in order to polish the parameters.

Reusing the same data as in the link $(c=300)$

$$\left( \begin{array}{cccc} x_i & y_i & z_i & t_i \\ 123 & 234 & 456 & 50.20 \\ 234 & 456 & 789 & 49.48 \\ 345 & 678 & 901 & 48.92 \\ 456 & 789 & 12 & 49.10 \\ 567 & 890 & 123 & 48.54 \\ 789 & 12 & 345 & 49.40 \\ 890 & 123 & 456 & 48.85 \end{array} \right)$$

the first step leads to $$X=1241.10 \qquad Y=995.36 \qquad Z=656.36 \qquad T=45.624$$

Using these values as estimates for the nonlinear regression leads to the following results $$\begin{array}{clclclclc} \text{} & \text{Estimate} & \text{Standard Error} & \text{Confidence Interval} \\ X & 1232.27 & 6.1325 & \{1205.88,1258.65\} \\ Y & 987.669 & 4.4095 & \{968.696,1006.64\} \\ Z & 653.595 & 2.6010 & \{642.365,664.824\} \\ T & 45.6817 & 0.0229 & \{45.5832,45.7803\} \\ \end{array}$$ The data were generated using $X=1234$, $Y=987$, $Z=654$, $T=45.678$ and the times were rounded to the second decimal place.

As you can see, the second step changes significantly the results.