I have $3$ ground-based radio receivers (call them $i$, $j$, and $k$), with known coordinates $x_i,y_i,z_i$. A transmitter (call it $m$) with unknown coordinates $x(t)_m, y(t)_m, z(t)_m$, moving with unknown velocity $v$, emits radio signal pulses with known velocity $c$ at an unknown, varying rate.
We may assume that the transmitter is traveling in a straight line. We know nothing about its acceleration.
The ground-based receivers record the time of signal arrival at each receiver for each pulse. I'm able to sort through these timestamps, and "pair" them up. So, I have a table that looks something like this
+---------+------------+------------+------------+
| Pulse # | Receiver i | Receiver j | Receiver k |
+---------+------------+------------+------------+
| 1 | 1100.12 | 1100.27 | 1100.28 |
| 2 | 1101.24 | 1101.39 | 1101.41 |
| ... | ... | ... | ... |
+---------+------------+------------+------------+
I've just come up with some random values, for illustrative purposes
I'm then, of course, able to determine the rate at which pulses are emitted (which may or may not be useful in solving my problem… I think not).
My goal is to determine the velocity of the transmitter.
What I've tried.
I'm familiar with such approaches as the MUSIC algorithm, however these all rely on directional antennae (which my receivers are not), the frequency of the signal, etc, which I do not have.
My first thought was to use Time Delay of Arrival (TDoA) techniques to determine the transmitter location at each pulse, and then use this (coupled with the timestamps) to approximate velocity.
I've solved this problem before (finding the location of a transmitter in $3D$ space given time of signal arrival at $n$ stations) using $4$ stations, however, as far as I can tell, a single solution using only three receivers is not possible in $3D$.
I recently posted a similar question, an answer to which claims that one can solve for two points, one of which will be the actual transmitter coordinates, and one of which will be a reflection of that point across the plane of the receivers.
Unfortunately, when one attempts to actually solve this system (below), it quickly reduces to a singular matrix. I've tried various solvers (those provided by GSL, as well as a few provided by Python's SciPy and SymPy packages) to no avail. Anecdotally, in speaking with an applied maths prof., I was told that the given system cannot be solved (an answer to this question says likewise). In any case, I'll keep trying.
$$\sqrt{(x – x_i)^2 + (y – y_i)^2 + (z – z_i)^2} + c \cdot (t_j – t_i) = \sqrt{(x – x_j)^2 + (y – y_j)^2 + (z – z_j)^2}$$
$$\sqrt{(x – x_j)^2 + (y – y_j)^2 + (z – z_k)^2} + c \cdot (t_k – t_j) = \sqrt{(x – x_k)^2 + (y – y_k)^2 + (z – z_k)^2}$$
$$\sqrt{(x – x_i)^2 + (y – y_i)^2 + (z – z_i)^2} + c \cdot (t_k – t_i) = \sqrt{(x – x_k)^2 + (y – y_k)^2 + (z – z_k)^2}$$
Ideally, this should give two solutions, one of which is correct, which would be fine for my purposes
Since what I'm truly after is the velocity of the transmitter, I've created this question in hopes there's a better (perhaps simpler) way to accomplish this.
Best Answer
Let the velocity of the transmitter be $\ \mathbf{v}\ $, the times at which it emits $\ T\ $ radio pulses be $\ t_1, t_2,\dots, t_T\ $, and its position at time $\ t=t_1\ $ be $\ \mathbf{x}_0\ $. All these quantities are unknown.
Let the positions of the $\ R\ $ receivers be $\ \mathbf{p}_r, r=1,2,\dots,R\ $ and the time when the $\ i^{\,\text{th}}\ $ pulse arrives at the $\ r^{\,\text{th}}\ $ receiver be $\ a_{ir}\ $. All these quantities are known.
The position of the transmitter at time $\ t\ $ will be $\ \mathbf{x}_0+\mathbf{v}\big(t-t_1\big)\ $, so the distance between the transmitter and the $\ r^{\,\text{th}}\ $ receiver at the time when it transmits the $\ i^{\,\text{th}}\ $ pulse will be $\ \big\|\big(t_i-t_1\big)\mathbf{v}+\mathbf{x}_0-\mathbf{p}_r\big\|\ $. The time the pulse will take to travel this distance is $\ c^{-1}\big\|\big(t_i-t_1\big)\mathbf{v}+\mathbf{x}_0-\mathbf{p}_r\big\|\ $, therefore it will arrive at the $\ r^{\,\text{th}}\ $ receiver at time $\ t_i+$$c^{-1}\big\|\big(t_i-t_1\big)\mathbf{v}+\mathbf{x}_0-\mathbf{p}_r\big\|\ $, which we know to be $\ a_{ir}\ $. Therefore, \begin{align} (1)\hspace{1em}c\big(a_{ir}-t_i\big)&=\big\|\big(t_i-t_1\big)\mathbf{v}+\mathbf{x}_0-\mathbf{p}_r\big\| \end{align} for $\ i=1,2,\dots,T\ $ and $\ r=1,2,\dots,R\ $. Squaring this gives \begin{align} (2)\hspace{1em}c^2\big(a_{ir}-t_i\big)^2&=\big\|\mathbf{v}\big\|^2\big(t_i-t_1\big)^2+2\big(t_i-t_1\big)\big\langle\mathbf{v},\mathbf{x}_0-\mathbf{p}_r\big\rangle+\big\|\mathbf{x}_0-\mathbf{p}_r\big\|^2\ . \end{align} In these equations, we can take $\ t_i,\ $$i=1,2,\dots,T\ $, $\big\langle\mathbf{v},\mathbf{x}_0-\mathbf{p}_r\big\rangle,\ $$\big\|\mathbf{x}_0-\mathbf{p}_r\big\|^2,\ $$r=1,2,\dots,R\ $, and $\ \big\|\mathbf{v}\big\|^2\ $ as the unknown real numbers whose values we need to recover. There are a total of $\ T+2R+1\ $ of these, and we have $\ TR\ $ equations available to determine them, so if $\ TR\ge T+2R+1\ $, then there should only be a small finite number of solutions. If $\ R\ge2\ $, this will be true if and only if $\ T\ge\frac{2R+1}{R-1}\ $. For $\ R=2\ $, $\ T=5\ $ will be enough to achieve this, while for $\ R=3\ $, $\ T=4\ $ will suffice.
Even if one of these conditions is satisfied, there's unlikely to be a convenient formula for the solutions, but it should be possible to obtain good approximations to them by a numerical algorithm, such as the Newton-Raphson method, for instance. When the system of equations is overdetermined $\big($that is, $\ TR>T+2R+1\ $, and all the equations are $\text{used}\big)$, it will, in practice, have no exact solution anyway, because imprecision in the measurements of $\ a_{ir}\ $ means that the values we have for them are only approximations. In this case, the Newton-Raphson method will find the values of the unknown quantities that minimise the sum of squares $$ \sum_{i=1}^T\sum_{r=1}^R\Big(c^2\big(a_{ir}-t_i\big)^2-\big\|\big(t_i-t_1\big)\mathbf{v}+\mathbf{x}_0-\mathbf{p}_r\big\|^2\Big)^2\ . $$ Once we have determined the values of all the unknown quantities, we can obtain the speed of the transmitter by taking the square root of $\ \|\mathbf{v}\|^2\ $, and we can try to use the now known values of $\ \|\mathbf{v}\|,$$\,\big\|\mathbf{x}_0-\mathbf{p}_r\big\|^2$ and $\ \big\langle\mathbf{v},\mathbf{x}_0-\mathbf{p}_r\big\rangle\ $ to recover possible values for $\ \mathbf{x}_0\ $ and $\ \mathbf{v}\ $.
However, when $\ R=2\ $, and $\ U:\mathbb{R}^3\rightarrow\mathbb{R}^3\ $ is any rotation about the line through $\ \mathbf{p}_1\ $ and $\ \mathbf{p}_2\ $, then $\ U(\mathbf{x})=\mathbf{u}_0+\Theta\,\mathbf{x}\ $ for some orthogonal matrix $\ \Theta\ $ with $\ \mathbf{p}_r-\Theta\mathbf{p}_r=\mathbf{u}_0\ $ for $\ r=1,2\ $ $\big($because $\ U(\mathbf{p})=\mathbf{p}\ $ for all points $\ \mathbf{p}\ $ lying on the line through $\ \mathbf{p}_1\ $ and $\ \mathbf{p}_2\ \big)$, and \begin{align} \big\|\big(t_i-t_1\big)\mathbf{v}+\mathbf{x}_0-\mathbf{p}_r\big\|&=\big\|\Theta\big(\big(t_i-t_1\big)\mathbf{v}+\mathbf{x}_0-\mathbf{p}_r\big)\big\|\\ &=\big\|\big(t_i-t_1\big)\Theta\mathbf{v}+\Theta\mathbf{x}_0-\Theta\mathbf{p}_r\big\|\\ &=\big\|\big(t_i-t_1\big)\Theta\mathbf{v}+\Theta\mathbf{x}_0+\mathbf{u}_0-\mathbf{p}_r\big\|\ . \end{align} Thus, if $\ \mathbf{v}\ $,$\ \mathbf{x}_0\ $ and $\ t_i\ $, $\ i=1,2,\dots,T\ $, is any solution for the system of equations $(1)$, then so is $\ \mathbf{v}'=\Theta\,\mathbf{v}\ $, $\ \mathbf{x}_0'=\Theta\,\mathbf{x}_0+\mathbf{u}_0\ $ and $\ t_i\ $, $\ i=$$1,2,$$\dots,T\ $. Thus, with only two receivers it's impossible to recover $\ \mathbf{v}\ $ or $\ \mathbf{x}_0\ $ completely. Nevertheless, with a sufficient number of received pulses, the speed $\ \|\mathbf{v}\|\ $ of the receiver should still be uniquely determined.
When $\ R=3\ $, and $\ \mathbf{p}_1$, $\mathbf{p}_2$, $\mathbf{p}_3\ $ are not collinear, we will have known values $\ v=\|\mathbf{v}\|\ $, \begin{align} d_r&=\big\|\mathbf{x}_0-\mathbf{p}_r\big\|\ ,\ \text{and}\\ e_r&=\big\langle\mathbf{v},\mathbf{x}_0-\mathbf{p}_r\big\rangle \end{align} for $\ r=1,2,3\ $. Since $\ \mathbf{p}_3-\mathbf{p}_1\ $ and $\ \mathbf{p}_3-\mathbf{p}_2\ $ are linearly independent, we can solve the linear equations \begin{align} d_3^2-d_1^2-\big\|\mathbf{p}_3\big\|^2+\big\|\mathbf{p}_1\big\|^2&=2\big\langle\mathbf{p}_1-\mathbf{p}_3,\mathbf{x}_0\big\rangle\\ d_3^2-d_2^2-\big\|\mathbf{p}_3\big\|^2+\big\|\mathbf{p}_2\big\|^2&=2\big\langle\mathbf{p}_2-\mathbf{p}_3,\mathbf{x}_0\big\rangle \end{align} to obtain $\ \mathbf{x}_0=\mathbf{a}+\lambda\mathbf{u}\ $, where $\ \mathbf{u}\ $ is a unit vector perpendicular to $\ \mathbf{p}_3-\mathbf{p}_1\ $ and $\ \mathbf{p}_3-\mathbf{p}_2\ $, $\ \mathbf{a}\ $ a known point lying in the plane of $\ \mathbf{p}_1$, $\mathbf{p}_2$, and $\mathbf{p}_3\ $, and $\ \lambda\ $ a real parameter whose value must be either $\ \sqrt{d_1^2-\big\|\mathbf{a}-\mathbf{p}_1\big\|^2}\ $ or $\ {-}\sqrt{d_1^2-\big\|\mathbf{a}-\mathbf{p}_1\big\|^2}\ $. If $\ e_r\ $ is value recovered for $\ \big\langle\mathbf{v},\mathbf{x}_0-\mathbf{p}_r\big\rangle\ $ from the numerical solution of the system of equations $(2)$, then for each of the two possible values $\ \mathbf{a}\pm\Big(\sqrt{d_1^2-\big\|\mathbf{a}-\mathbf{p}_1\big\|^2}\Big)\mathbf{u}\ $ of $\ \mathbf{x}_0\ $, we can now solve the linear equations $$ e_r=\big\langle\mathbf{v},\mathbf{x}_0-\mathbf{p}_r\big\rangle $$ for $\ r=1,2,3\ $ to obtain the corresponding value of $\ \mathbf{v}\ $. The solution of either of these systems will always be the reflection of the solution of the other in the plane of $\ \mathbf{p}_1$, $\mathbf{p}_2$, and $\mathbf{p}_3\ $.