High-Precision Time Difference of Arrival Estimation Method Based on Phase Measurement

Abstract

In unmanned aerial vehicle (UAV)-based time difference of arrival (TDOA) positioning technique, baselines are limited due to communication constraints. In this case, the accuracy is highly sensitive to the TDOA measurements’ error. This article primarily addresses the problem of short-baseline high-precision time synchronization and TDOA measurement. We conducted a detailed analysis of error models in TDOA systems, considering both the time and phase measurement. We utilize the frequency division wireless phase synchronization technique in TDOA systems. Building upon this synchronization scheme, we propose a novel time delay estimation method that relies on phase measurements based on the integer least squares method. The performance of this method is demonstrated through Monte Carlo simulations and outdoor experiments. The standard deviations of synchronization and TDOA measurements in experiments are 1.12 ps and 1.66 ps, respectively. Furthermore, the circular error probable (CEP) accuracy is improved from 0.33%R to 0.02%R, offering support for the practical application of distributed short-baseline high-precision passive location techniques.

1. Introduction

Unmanned aerial vehicle (UAV)-based passive location has always been a research hotspot in the academic and civil fields. Compared with the active location, it has the advantages of low power consumption and longer positioning distance. With the increasing maturity of low-cost commercial-off-the-shelf (COTS) software-defined radio (SDR), radio frequency system-on-chip (RFSoC), and other products [1], the cost of distributed passive location systems has reduced, making it highly valuable in search and rescue [2], intelligent logistics [3], urban positioning [4], and other fields [5,6].

Time difference of arrival (TDOA) localization offers a higher positioning accuracy compared to other localization techniques [7,8,9,10,11]. Passive localization systems capture the same pulse signal and estimate the TDOA. When the baselines between nodes are long in comparison to the distance to the emitter, current global navigation satellite system (GNSS)-based time synchronization and position techniques meet the system requirements. However, when the operational range of localization system is limited, the baselines are constrained by the emitter beamwidth or UAV communication range. Under such conditions, TDOA positioning accuracy is highly susceptible to TDOA measurements. To maintain the same localization performance compared with long-baseline TDOA, the requirements for TDOA measurements are more stringent [12,13]. To achieve high-precision TDOA positioning, three problems must be taken into consideration: high-precision relative positioning, high-precision time synchronization, and high-precision time delay estimation (TDE). In this article, we focus on high-precision time delay estimation and high-precision time synchronization techniques.

The issue of achieving high-precision wireless time synchronization has always been a significant challenge in various fields. In wireless sensor networks, time synchronization is composed of centralized methods and distributed methods [14,15,16], where the synchronization convergence speed of a large-scale network is the main concern. On the other hand, in multi/bistatic radar, where the network scale is smaller, synchronization accuracy takes precedence. To enhance accuracy, the bistatic SAR satellites [17,18,19] adopt a time division and co-frequency method to exchange pulse signals, achieving time and phase synchronization between two satellites. Building upon this, some studies in the literature have designed multi-node wireless synchronization protocol [11,12,20], improving the efficiency of a two-way time transfer for wireless nodes via time division broadcasting. In active distributed radar, since the radar waveform and the transmission time are known in advance, the radar mode and synchronous mode are avoided by the time division design or orthogonal waveform design. However, in passive localization applications, the waveform and launch time of the emitter are unknown; if the synchronization link frequency band overlaps with the emitter frequency band, it may introduce interference to the localization system.

Besides the time synchronization error, the accuracy of TDOA estimation also determines the localization performance. When the signal-to-noise (SNR) ratio and the bandwidth are fixed, the time delay estimation (TDE) method plays a crucial role in approaching the performance bounds in TDOA estimation. Most traditional methods use the cross-correlation function (CCF) or general cross-correlation (GCC) function of peak position estimation methods, such as the sinc interpolation method [21], the quadratic least squares (QLS) method [22,23], the sinc nonlinear least squares (sinc-LS) method [24], and the matched filter least squares (MFLS) method [25]. The pursuit of high-precision TDE is essential not only in TDOA localization, but also in many other fields. For instance, in real-time kinematic (RTK) technique, precise carrier phase measurements are introduced to realize high positioning [26,27]. However, most current methods in TDOA only utilize the peak position information, thereby underutilizing the full performance potential of the localization system.

The contributions in the paper include the following:

(1)

We conducted short-baseline TDOA experiments using a precise time and phase wireless synchronization method based on frequency division.

(2)

We propose a TDOA estimation using phase measurements based on the integer least squares method.

(3)

We verify the effect of frequency instability on peak position-based estimation through model derivation and experiments.

This article is organized as follows. In Section 2, the error sources of the distributed TDOA system link are modeled in detail from the two dimensions of time observation and phase observation. In Section 3, error analysis is conducted on time observation and phase observation for the synchronization system link, and a frequency division synchronization scheme is designed for distributed passive localization systems. In Section 4, a high-precision delay estimation method based on phase measurement is introduced. In Section 5, we outline the principal tests and detail a three-sensor TDOA localization experiment. The experimental results align closely with the signal models outlined in Section 2 and Section 3. The standard deviations of synchronization and TDOA measurements in the experiments are 1.12 ps and 1.66 ps, respectively. Furthermore, the CEP accuracy is improved from 0.33%R to 0.02%R. Finally, Section 6 concludes the article.

2. System Models

In the scenario described in Figure 1, after the disaster, N unmanned aerial vehicles (UAVs) are settled to search survivors by locating their portable station. The distance between UAVs is limited by communication and synchronization link as well as the portable station antenna beamwidth. Synchronization errors between UAVs are compensated through wireless synchronization links. In this section, we will analyze the error model of the TDOA system to provide a foundation for high-precision TDOA estimation.

2.1. Hardware Block Model

For simplicity, we consider the hardware model for two sensors as shown in Figure 2. For each sensor, the transmitted signal is received by CH1. Since the TDOA system receives the same pulse from the emitter, we ignore the effects of non-idealities and the initial phase of the radiation signal. Then, sensor i receives the radiation signal after time delay, ${\tau }_{tof}^i$, at the antenna phase center as $s_{ant}^i(t)$:

$$s_{ant}^i(t)=s_{BD}\left(t-{\tau }_{tof}^i\right)\exp \left(j2\pi f_c\left(t-{\tau }_{tof}^i\right)\right)$$

(1)

where $s_{BD}\left(t\right)$ is the baseband signal and $f_c$ is the carrier frequency.

2.2. Signal Model with Different Oscillators

Considering the short-term instability of the crystal frequency of the TDOA system, we assume that the instantaneous frequency remains constant within the pulse width of the emitter, but varies in different pulses [19]. Define $f_{LO}^i\left(t_0\right)$ as the instantaneous local frequency of sensor i at time $t_0$. We have

$$f_{LO}^i\left(t_0\right)=f_{LO}+\delta f_{LO}^i\left(t_0\right)$$

(2)

where $f_{LO}$ is the nominal LO frequency and $\delta f_{LO}^i$ is the time-variant frequency instability. The phase noise of the LO signal ${\theta }_{pn}^i\left(t_0\right)$ at time $t_0$ can be calculated as follows [28]:

$${\theta }_{pn}^i\left(t_0\right)={\displaystyle {\int }_0^{t_0}2\pi \delta f_{LO}^i\left(\eta \right)d\eta }.$$

(3)

So, when the radiation signal travels to the mixer at time $t_{mix}^i$, the LO signal $s_{LO}^i\left(t\right)$ is

$$s_{LO}^i\left(t\right)=\exp \left[j2\pi \delta f_{LO}^i\left(t_{mix}^i\right)\cdot t+j2\pi f_{LO}t\right]\cdot \exp \left[j{\theta }_{pn}^i\left(t_{mix}^i\right)+j{\theta }_{ini}^i\right]$$

(4)

where ${\theta }_{ini}^i$ is the initial phase of the LO signal. The instrument delay and phase shift are defined as ${\tau }_{inst}^i$ and ${\theta }_{inst}^i$, respectively [29,30]. The intermediate frequency signal $s_{IF}^i(t)$ is down-converted with $s_{LO}^i(t)$ as follows:

$$\begin{array}{c}{s}_{IF}^i(t)=s_{BD}\left(t-{\tau }_{tof}^i-{\tau }_{inst}^i\right)\exp \left[-j2\pi f_c{\tau }_{tof}^i+j{\theta }_{inst}^i\right]\\ ·\exp \left[j2\pi f_{IF}^i\left(t_{mix}^i\right)t-j{\theta }_{pn}^i\left(t_{mix}^i\right)-j{\theta }_{ini}^i\right]\end{array}$$

(5)

where the intermediate frequency $f_{IF}^i$ (omit $t_{mix}^i$ for simplicity) is

$$f_{IF}^i=f_c-f_{LO}^i\left(t_{mix}^i\right).$$

(6)

After $s_{IF}^i(t)$ is sampled by the realistic but unknown clock, we have to reconstruct the continuous signal with the nominal clock. Under the assumption from Equation (2), we ignore the rapid time jitter in a pulse as well as the quantization error. The sample point maps to the ideal time stamps [31]. The time synchronization error is defined as positive when the realistic clock lags the ideal clock. In this case, the reconstructed signal will be ahead, which is equivalent to the translation and expansion transformation of the actual sampled signal [32], as graphically depicted in Figure 3a. In contrast, for the DAC signal, as we have no knowledge about the realistic clock, we have to calculate the signal sample point value under the ideal clock. In this case, the generated signal will lag the wanted signal, as shown in Figure 3b. Note that the difference between nominal frequency and realistic frequency is exaggerated to shown the effect.

We assume that the frequency bias, when averaged over seconds and compared to the nominal frequency, is negligible due to locking it to the GNSS one-pulse-per-second (1 PPS) signal [33]. Consequently, we only concentrate on the time synchronization error. Define the initial time synchronization error $\Delta {\tau }_{clk}^i$ between the actual sampling clock and the nominal sampling clock during ADC sampling. After ADC sampling and digital down-conversion, we then obtain the baseband signal $s_i(t)$:

$$\begin{array}{cc}\hfill s_i(t)& =s_{IF}^i\left(t+\Delta {\tau }_{clk}^i\right)\cdot \exp \left(-j2\pi f_{IF}t\right)\hfill \\ & =s_{BD}\left(t+\Delta {\tau }_{clk}^i-{\tau }_{tof}^i-{\tau }_{inst}^i\right)·\exp \left[-j2\pi f_c{\tau }_{tof}^i\right]\hfill \\ & ·\exp \left[j2\pi f_{IF}\Delta {\tau }_{clk}^i\right]\cdot \exp \left[-j2\pi \delta f_{LO}^i\left(t+\Delta {\tau }_{clk}^i-{\tau }_{tof}^i-{\tau }_{inst}^i\right)\right]\hfill \\ & ·\exp \left[j{\theta }_{inst}^i-j{\theta }_{pn}^i\left(t_{mix}^i\right)-j{\theta }_{ini}^i\right]\hfill \end{array}.$$

(7)

If the clock generator and the frequency synthesizer are driven by the same oscillator, the relationship between the time synchronization error and phase error can be expressed as follows [28]:

$$\frac{d\Delta {\tau }_{clk}^i}{dt}=-\frac{1}{2\pi f_{LO}}·\frac{d{f}_{pn}^i}{dt}=-\frac{\delta f_{LO}^i}{f_{LO}},$$

(8)

$$\Delta {\tau }_{clk}^i=-\frac{1}{f_{LO}}{\displaystyle {\int }_0^t\delta f_{LO}^i\left(\eta \right)d\eta +C_0}=-\frac{1}{2\pi f_{LO}}{\theta }_{pn}^i+C_0$$

(9)

where $C_0$ is an unknown constant. Since we are interested in the time-variant element, the phase error can be simplified as

$$2\pi f_{IF}\cdot \Delta {\tau }_{clk}^i\left(t_0\right)-{\theta }_{pn}^i\left(t_0\right)=-\frac{f_c}{f_{LO}}{\theta }_{pn}^i\left(t_0\right).$$

(10)

Substitute Equation (10) with Equation (7), to obtain $s_i(t)$:

$$\begin{array}{cc}\hfill s_i(t)& =s_{BD}\left(t+\Delta {\tau }_{clk}^i-{\tau }_{tof}^i-{\tau }_{inst}^i\right)\cdot \exp \left(-j2\pi f_c{\tau }_{tof}^i\right)\hfill \\ & \cdot \exp \left[-j\cdot f_c/f_{LO}\cdot {\theta }_{pn}^i\left(t_{mix}^i\right)\right]\cdot \exp \left[-j2\pi \delta f_{LO}^i\left(t+\Delta {\tau }_{clk}^i-{\tau }_{tof}^i-{\tau }_{inst}^i\right)\right]\hfill \\ & \cdot \exp \left[j{\theta }_{inst}^i-j{\theta }_{ini}^i\right]\hfill \end{array}.$$

(11)

Similarly, the data $s_j\left(t\right)$ received by sensor j are

$$\begin{array}{cc}\hfill s_j(t)& =s_{BD}\left(t+\Delta {\tau }_{clk}^j-{\tau }_{tof}^i-{\tau }_{inst}^i\right)\cdot \exp \left(-j2\pi f_c{\tau }_{tof}^j\right)\hfill \\ & \cdot \exp \left[-j\cdot f_c/f_{LO}\cdot {\theta }_{pn}^j\left(t_{mix}^j\right)\right]\cdot \exp \left[-j2\pi \delta f_{LO}^j\left(t+\Delta {\tau }_{clk}^j-{\tau }_{tof}^i-{\tau }_{inst}^i\right)\right]\hfill \\ & \cdot \exp \left[j{\theta }_{inst}^j-j{\theta }_{ini}^j\right]\hfill \end{array}.$$

(12)

The received signal experiences additional time delays and phase shifts caused by the instruments and time synchronization error. Furthermore, there is an instantaneous frequency bias over the pulse width due to the instability of oscillators.

2.3. The Effect of Different Oscillators on TDOA Estimation

In general, we estimate TDOA by measuring the peak position of the CCF. $R_{ij}\left(\tau \right)$ is defined as the CCF [34]:

$$R_{ij}\left(\tau \right)={\displaystyle {\int }_{-\infty }^{\infty }{s}_i\left(t\right)\cdot s_j^{*}\left(t-\tau \right)dt}.$$

(13)

The detailed derivation is shown in Appendix A. Here, we present the results directly. The amplitude of $R_{ij}\left(\tau \right)$ is

$$\left|R_{ij}\left(\tau \right)\right|=rect\left(\frac{\tau -\gamma }{2T}\right)\cdot \left(T-\left|\tau -\gamma \right|\right)\cdot \sin \mathrm{c}\left(\pi \left[k\left(\tau -\gamma \right)-\delta f_{LO}^{ij}\left(t_{mix}^i,t_{mix}^j\right)\right]\left(T-\left|\tau -\gamma \right|\right)\right)$$

(14)

where $\gamma ={\tau }_{tof}^i-{\tau }_{tof}^j+{\tau }_{inst}^i-{\tau }_{inst}^j-\Delta {\tau }_{clk}^i+\Delta {\tau }_{clk}^j$ is the time delay difference of the received signals between sensor i and sensor j and $\delta f_{LO}^{ij}\left(t_{mix}^i,t_{mix}^j\right)$ is the difference of the instantaneous frequency at $t_{mix}^i$ and $t_{mix}^j$, which will be simplified as $\delta f_{LO}^{ij}$ in the following sections. Denote ${(\cdot )}^{ij}={(\cdot )}^i-{(\cdot )}^i$ and ${(\cdot )}^{i+j}={(\cdot )}^i+{(\cdot )}^i$. The position of the CCF peak is the maximum likelihood estimation of the delay. The CCF peak position ${\tau }_{pk}^{ij}$ is

$${\tau }_{pk}^{ij}\approx \gamma +\frac{\delta f_{LO}^{ij}}{k}.$$

(15)

The CCF peak angle ${\theta }_{pk}^{ij}$ is

$${\theta }_{pk}^{ij}=Ang\left\{R_{ij}\left({\tau }_{pk}^{ij}\right)\right\}=-2\pi f_c{t}_{mix}^{ij}+\pi \delta f_{LO}^{i+j}\left(-\frac{\delta f_{LO}^{ij}}{k}\right)-\frac{f_c}{f_{LO}}{\theta }_{pn}^{ij}\left(t_0\right)+\left({\theta }_{inst}^{ij}-{\theta }_{ini}^{ij}\right).$$

(16)

Denote the time synchronization compensation as $\Delta {\tau }_{clk,c}^{ij}$. After the delay calibration and time synchronization of the instruments are performed, the TDOA estimation, $TDO{A}_{\tau }^{ij}$, derived from the CCF peak position is

$$TDO{A}_{\tau }^{ij}={\tau }_{tof}^{ij}+\frac{\delta f_{LO}^{ij}}{k}-\left(\Delta {\tau }_{clk}^{ij}-\Delta {\tau }_{clk,c}^{ij}\right).$$

(17)

Denote the phase synchronization compensation as ${\theta }_{pn,c}^{ij}\left(t_0\right)$. After conducting the phase shift calibration and time (phase) synchronization of the instruments, the TDOA estimation, $TDO{A}_{\theta }^{ij}$, derived from the CCF peak phase is

$$\begin{array}{cc}\hfill TDO{A}_{\theta }^{ij}& =\frac{1}{2\pi f_c}\left\{-{\theta }_{pk}^{ij}+2N_{ij}\pi -\frac{f_c}{f_{LO}}{\theta }_{pn}^{ij}\left(t_0\right)+{\theta }_{pn,c}^{ij}\left(t_0\right)\right\}-\frac{\delta f_{LO}^{i+j}}{2f_c}\cdot \frac{\delta f_{LO}^{ij}}{k}\hfill \\ & \approx \frac{1}{2\pi f_c}\left\{-{\theta }_{pk}^{ij}+2N_{ij}\pi -\frac{f_c}{f_{LO}}{\theta }_{pn}^{ij}\left(t_0\right)+{\theta }_{pn,c}^{ij}\left(t_0\right)\right\}\hfill \end{array}$$

(18)

where $N_{ij}\in Z$ is the integer ambiguity. From Equation (17), the term $\delta f_{LO}^{ij}/k$ degrades the performance of peak position-based methods because of the instability of oscillators or even VCOs. Since the order of signal pulse widths is in μs, $\delta f_{LO}^{ij}$ is corelated with the Allan Deviation (ADEV) at the average time of μs levels. According to the relationship between the ADEV ${\sigma }_y\left(\tau \right)$ and the phase noise spectrum $L\left(f\right)$ [35],

$${\sigma }_y^2\left(\tau \right)=\frac{4}{{\pi }^2{f}_c^2}{\displaystyle {\int }_{f_l}^{f_h}L\left(f\right)H\left(f\right)df}$$

(19)

where $H\left(f\right)={\sin }^4\left(\pi f\tau \right)/{\tau }^2$, $\tau$ is the average time to calculate the instantaneous frequency, and $f_l$ and $f_h$ are the lower and upper cutoff frequencies respectively. Figure 4 shows $H\left(f\right)$ when $\tau =10$ μs, and the main error of $\delta f_{LO}^{ij}$ is located in $L\left(f\right)$ at about 50 kHz, which is a fast noise-like error. It means that $\delta f_{LO}^{ij}$ changes rapidly and only a few coherences remain among pulses, while in Equation (18), both the elements $\delta f_{LO}^{i+j}/\left(2f_c\right)$ and $\delta f_{LO}^{ij}/k$ are small and can be ignored. In Section 5.3, we will show the experiment results.

3. Time Synchronization Method

In this section, we take the advantage of the time-variant phase error measurement from the synchronization link and compensate for the TDOA phase measurement.

3.1. Synchronization Scheme

The synchronous link is shown in Figure 2, and each sensor utilizes CH2 to exchange the synchronization signal. To ensure that CH1 can continuously detect the emitter signal, the synchronization link uses another frequency band to observe the synchronization error between nodes. The timing diagram is shown in Figure 5. Each sensor periodically transmits the synchronization pulse signal, and the launch time of each sensor is staggered to ensure that only one sensor transmits the signal each time. At time $t$, sensor 1 transmits a synchronization pulse with the pulse width of $T_p$, and sensor 2 and 3 receives it. After waiting for a fixed delay ${\tau }_{sys}$, sensor 2 transmits a synchronization signal with the pulse width of $T_p$, and sensor 1 and 3 receives it. After another fixed delay ${\tau }_{sys}$, sensor 3 transmits the synchronization signal and sensor 1 and 2 receives it. The above process is repeated with pulse repetition time (PRT). The synchronization process continues throughout the entire working time.

3.2. Compensation Time and Phase

The process of the DAC transmitting the synchronization signal is similar to that of ADC sampling. After transmitting, the actual wave is equivalent to a translation and scaling transformation from the designed one, as shown in Figure 3b. The signal transmitted by the DAC is

$$s_{DA}^i\left(t\right)=s_{BD}\left(t-\Delta {\tau }_{clk}^i\right)·\exp \left[j2\pi f_{IF}\left(t-\Delta {\tau }_{clk}^i\right)\right].$$

(20)

The synchronization process is similar to the previous section. So, we present the result directly. After receiving the signal and pulse compression processing, the CCF is

$$\begin{array}{cc}\hfill R^{i\to j}\left(\tau \right)& =rect\left(\frac{\tau }{2T}\right)\left(T-\left|\tau \right|\right)\sin \mathrm{c}\left[\pi \left(k\tau +\delta f_{LO2}^{ij}\right)\left(T-\left|\tau \right|\right)\right]\hfill \\ & \cdot \exp \left[-j2\pi f_{LO2}\cdot {\tau }_{tof}^{i\to j}\right]\hfill \\ & \cdot \exp \left(j\pi \delta f_{LO2}^{ij}\tau \right)\hfill \\ & \cdot \exp \left[j\frac{f_c}{f_{LO2}}\left({\theta }_{pn}^{i,ch2}-{\theta }_{pn}^{j,ch2}\right)\right]\hfill \\ & \cdot \exp \left(j{\theta }_{ini}^{i,ch1}-j{\theta }_{ini}^{j,ch2}+j{\theta }_{inst,Tx}^{i,ch2}+j{\theta }_{inst,Rx}^{j,ch2}\right)\hfill \end{array}.$$

(21)

The peak position estimation and the peak phase of sensor i are

$${\tau }_{pk}^{i\to j}={\tau }_{tof}^{i\to j}+\Delta {\tau }_{clk}^{ij}-\frac{\delta f_{LO2}^{ij}}{k}+{\tau }_{inst,Tx}^{i,ch2}+{\tau }_{inst,Rx}^{j,ch2},$$

(22)

$$\begin{array}{cc}\hfill {\theta }_{pk}^{i\to j}& =-j2\pi f_{LO2}\cdot {\tau }_{tof}^{i\to j}+\pi \delta f_{LO2}^{ij}\frac{\delta f_{LO2}^{ij}}{k}\hfill \\ & +\frac{f_c}{f_{LO2}}\left({\theta }_{pn}^{i,ch2}-{\theta }_{pn}^{j,ch2}\right)+{\theta }_{inst,Tx}^{i,ch2}+{\theta }_{inst,Rx}^{j,ch2}+{\theta }_{ini}^{i,ch1}-{\theta }_{ini}^{j,ch2}\hfill \end{array}.$$

(23)

Similarly, when sensor j transmits a signal to sensor i, we obtain the peak position estimation ${\tau }_{pk}^{j\to i}$ and the peak phase ${\theta }_{pk}^{j\to i}$. Assuming the instrument time delay (phase shift) as well as initial phase have been calibrated, finally, we obtain the time synchronization compensation $\Delta {\tau }_{clk,c}^{ij}$ and phase synchronization compensation ${\theta }_{pn,c}^{ij}\left(t_0\right)$:

$$\Delta {\tau }_{clk,c}^{ij}=\frac{{\tau }_{pk}^{i\to j}-{\tau }_{pk}^{j\to i}}{2}=\Delta {\tau }_{clk}^{ij}+\frac{1}{2k}\Delta \delta f_{LO2}^{ij}\left(t_0,{\tau }_{sys}\right),$$

(24)

$${\theta }_{pn,c}^{ij}=\frac{{\theta }_{pk}^{i\to j}-{\theta }_{pk}^{j\to i}}{2}=\frac{f_c}{f_{LO}}{\theta }_{pn}^{ij}\left(t_0\right)+\Delta {\theta }_{sycn}^{ij}\left(t_0,{\tau }_{sys}\right)$$

(25)

where $\Delta \delta f_{LO2}^{ij}\left(t_0,{\tau }_{sys}\right)=\left[\delta f_{LO2}^{ij}\left(t_0\right)-\delta f_{LO2}^{ji}\left(t_0+{\tau }_{sys}\right)\right]$ represents the instantaneous frequency difference changes during the interval ${\tau }_{sys}$, and $\Delta {\theta }_{sycn}^{ij}\left(t_0,{\tau }_{sys}\right)=\pi f_c/f_{LO2}\cdot {\displaystyle {\int }_{t_0}^{t_0+{\tau }_{sys}}\delta f_{LO2}^{ij}\left(\eta \right)d\eta }$ is the phase noise change during the interval ${\tau }_{sys}$ of the time division duplex at the time instant $t_0$.

Substitute Equation (24) with Equation (17), and Equation (25) with Equation (18); we then have the TDOA estimation after compensation by the synchronization link:

$$TDO{A}_{\tau }^{ij}={\tau }_{tof}^{ij}+\frac{\delta f_{LO}^{ij}}{k}+\frac{1}{2k}\Delta \delta f_{LO2}^{ij}\left(t_0,{\tau }_{sys}\right),$$

(26)

$$TDO{A}_{\theta }^{ij}=\frac{1}{2\pi f_c}\left\{-{\theta }_{pk}^{ij}+2N_{ij}\pi +\Delta {\theta }_{sycn}^{ij}\left(t_0,{\tau }_{sys}\right)\right\}.$$

(27)

Comparing Equation (26) with Equation (27), Equation (26) has an extra error due to the unequal instantaneous frequencies $\delta f_{LO2}^{ij}$ and $\delta f_{LO2}^{ji}$. According to Figure 5, the exchange is a time division duplex and the residual error is affected by the coherence of oscillators during ${\tau }_{sys}$ once the oscillator is determined.

4. High-Precision Time Delay Estimation

In this section, we concentrate on the time delay estimation using phase measurements. Before we clarify the problem of TDOA measurements, we add some constraints to simplify the problem. We assume that the phase can be correctly unwrapped, which means that the carrier phase change between consecutive pulses is less than π. Furthermore, we assume that the carrier frequency of the portable station is known. So, when the TDOA system can capture the signal continuously, the integer ambiguity should be a constant value. Once the integer ambiguity is solved, high-precision TDOA estimation can be calculated according to Equation (27).

In the following, the proposed integer ambiguity solution as well as the Cramer–Rao lower bound (CRLB) are present. In Section 4.3, we will analyze the performance and the constraints through simulations.

4.1. Ambiguity Integer Solution

Denote $M$ as the captured pulse number. Given $M$ coarse estimations $TDO{A}_{\tau }^{ij}$ and $M$ ambiguous phase measurements ${\theta }_{pk}^{ij}$, the integer ambiguity problem can be defined as

$$\underset{N}{\min }{‖TDO{A}_{\tau }^{ij}-\frac{{\theta }_{pk}^{ij}}{2\pi f_c}-\frac{N}{f_c}‖}^2,\hspace{0.33em}N\in \mathbb{Z}.$$

(28)

We use the integer least squares (LS) estimation to solve the problem; the estimation consists of two steps:

(1)

Obtain the float solution $\widehat{N}$ with a standard LS problem:

$$\underset{N}{\min }{‖TDO{A}_{\tau }^{ij}-\frac{{\theta }_{pk}^{ij}}{2\pi f_c}-\frac{N}{f_c}‖}^2,\hspace{0.33em}N\in \mathbb{R}.$$

(29)

By taking the first the derivative of the above equation and by setting it to zero, we have the float LS solution:

$$\widehat{N}=\frac{f_c}{M}\cdot {\displaystyle \sum _{i=1}^M\left[TDO{A}_{\tau }^{ij}\left(i\right)-\frac{{\theta }_{pk}^{ij}\left(i\right)}{2\pi f_c}\right]}.$$

(30)

(2)

Integer ambiguity $\tilde{N}$ is calculated as

$$\tilde{N}=\left[\widehat{N}\right]$$

(31)

where $\left[ \right]$ means that the value is rounded to the nearest integer.

The solution can only be worked out after all $M$ pulses are captured. We also present a sequential LS solution:

$$\left\{\begin{array}{cc}\widehat{N}\left(1\right)=f_c\cdot \left[TDO{A}_{\tau }^{ij}\left(1\right)-\frac{{\theta }_{pk}^{ij}\left(1\right)}{2\pi f_c}\right],& m=1\hfill \\ \widehat{N}\left(m\right)=f_c\cdot \frac{m-1}{m}\widehat{N}\left(m-1\right)+\frac{f_c}{m}\left[TDO{A}_{\tau }^{ij}\left(m\right)-\frac{{\theta }_{pk}^{ij}\left(m\right)}{2\pi f_c}\right],& m\ge 2\hfill \end{array}\right.,$$

(32)

$$N\left(m\right)=\left[\widehat{N}\left(m\right)\right].$$

(33)

Note that when the mth pulse arrives, $N\left(m\right)$ might change, and all the first m TDOA estimations ${\tau }_{\phi }\left(m\right)$ will change.

In (28), the process is similar to the integer ambiguity resolution of double-difference carrier phase in RTK technique [26,27]. We utilize the peak position-based measurement to replace pseudorange measurements in RTK. However, there are a few differences between the TDOA model and the double-difference model. The phase shift of satellites and receivers can be eliminated through double-difference from the same epoch, leaving only the integer phase ambiguity [30]. It is difficult for passive localization to have a collaborative emitter all the time. Therefore, the system should be consistent in design and calibrated before use.

4.2. Cramer–Rao Lower Bound

The peak position and phase of CCF are extracted from Equation (46). To simplify the problem, we set the instantaneous frequency $\delta f_{LO}^{ij}=0$. If the TDE is derived from the peak position, the CRLB [19] is

$${\sigma }_{pos}^2\ge \frac{3}{2{\pi }^2{B}^2\cdot f_s\cdot T_p\cdot SNR}.$$

(34)

If the TDE is derived from the peak phase, the CCF can be simplified as

$$R\left[n\right]=\left|R\left[n\right]\right|\cdot \exp \left(j\phi \right).$$

(35)

When $R\left[n\right]$ is superimposed with Gaussian white noise with a variance of ${\sigma }_N^2$, the CRLB of the phase measurement [36] is

$${\sigma }_{\phi }^2\ge \frac{{\sigma }_N^2/2}{{\displaystyle {\sum }_{n=0}^{N-1}{\left|\frac{\partial s}{\partial \phi }\right|}^2}}.$$

(36)

Assuming that the sampling interval is sufficiently small, Equation (36) is approximated as

$${\sigma }_{\phi }^2\ge \frac{{\sigma }_N^2/2}{{\displaystyle {\sum }_{n=0}^{N-1}{\left|\frac{\partial s}{\partial \phi }\right|}^2}}\approx \frac{{\sigma }_N^2/2}{\frac{1}{\Delta t}·{\displaystyle {\int }_0^{T_p}{\left|\frac{\partial s}{\partial \phi }\right|}^{2}dt}}.$$

(37)

Denote ${\displaystyle {\int }_0^{T_p}{\left|\frac{\partial s}{\partial \phi }\right|}^{2}dt}={\displaystyle {\int }_0^{T_p}{\left|s\right|}^{2}dt}=P_{av}\cdot T_p$, $f_s=1/\Delta t$ and signal-to-noise ratio as $SNR=P_{av}/{\sigma }_N^2$, we then have

$${\sigma }_{\phi }^2\ge \frac{{\sigma }_N^2}{2f_s{P}_{av}\cdot T_p}=\frac{1}{2f_s\cdot T_p\cdot SNR}.$$

(38)

If the integer ambiguity is solved correctly, then the CRLB of the time delay estimation is

$${\sigma }_{pha}^2={\left(\frac{1}{2\pi f_c}\right)}^2{\sigma }_{\phi }^2\ge {\left(\frac{1}{2\pi f_c}\right)}^2\cdot \frac{1}{2f_s\cdot T_p\cdot SNR}.$$

(39)

Assume the signal bandwidth is 10 MHz and the carrier frequency is 3.6 GHz, the ratio of the upper bound of the performance of the two methods is

$$10{\log }_{10}\left(\frac{{\sigma }_{pha}}{{\sigma }_{pos}}\right)=10{\log }_{10}\left(\frac{\sqrt{3}}{6}\cdot \frac{B}{f_c}\right)\approx -30.96\hspace{0.33em}\mathrm{d}\mathrm{B}.$$

(40)

4.3. Anaylsis and Simulations

The accuracy of the proposed method highly relies on the biased error of peak position-based TDE $TDO{A}_{\tau }^{ij}$. Apart from noises, the bias error is inevitable due to a finite number of coefficients and discretization [21]. When the CCF peak locates between two adjacent sample points, the bias error can be reduced by interpolation methods, such as the sinc interpolation method [21], QLS [22,23], and sinc-LS [24]. QLS uses a parabola to fit the CCF main lobe, for instance, the CCF peak and the two adjacent points, and calculates the peak of the parabola as the TDE. Sinc-LS uses the sinc function to fit the CCF, but it requires iteration. Details of those algorithms are described in the cited references.

Different from those peak position-based methods, Hatch filter [37,38] is a kind of smooth algorithm combining the pseudorange and carrier phase in GNSS position, which is

$$\left\{\begin{array}{cc}TDO{A}_{hatch}^{ij}\left(n\right)=\frac{1}{n}TDO{A}_{\tau }^{ij}\left(n\right)+\frac{n-1}{n}\left[TDO{A}_{hatch}^{ij}\left(n-1\right)+\frac{{\theta }_{pk}^{ij}\left(m-1\right)-{\theta }_{pk}^{ij}\left(m\right)}{2\pi f_c}\right],& n<M\\ TDO{A}_{hatch}^{ij}\left(n\right)=\frac{1}{M}TDO{A}_{\tau }^{ij}\left(n\right)+\frac{M-1}{M}\left[TDO{A}_{hatch}^{ij}\left(n-1\right)+\frac{{\theta }_{pk}^{ij}\left(m-1\right)-{\theta }_{pk}^{ij}\left(m\right)}{2\pi f_c}\right],& n\ge M\end{array}\right..$$

(41)

The Monte Carlo simulations are carried out to evaluate the performance of the TDOA estimation. To this end, a radio signal is received by a two-node TDOA system with different SNRs. The simulation parameters are shown in

Table 1. Monte Carlo simulation parameters.

Parameter	Value
Waveform	chirp
Signal Bandwidth	10 MHz
Center Frequency	3.6 GHz
Pulse Repetition Time	1 ms
Signal Pulse Width	3 μs
Receiving Pulse	5 μs
Sample Rate	2400 MHz
Pulse Num per Trial	100
TDOA Change Rate	−27.8 ns/s

. To ensure that the adjacent pulse phase difference is less than $\pi$, we set the TDOA change rate to be −27.8 ns/s, which results in an adjacent pulse phase difference of −0.2$\pi$ radians. The simulation has 100 trials and each trial has 100 pulses. We choose the root mean square (RMS) of the difference between TDOA estimations and TDOA true values, which is

$$RMS=\sqrt{\frac{1}{M}{\displaystyle \sum _{n=1}^M{\left|\widehat{TDO{A}_n^{ij}}-TDO{A}_n^{ij}\right|}^2}}.$$

(42)

According to Equations (34) and (39), the root mean square (RMS) of the phase method is about two orders of magnitude smaller than that of the peak position method. Besides the proposed method, we choose the peak position-based methods and Hatch filter mentioned above. In the Hatch filter method, we use the sinc-LS results as the unambiguous TDE. The general value of SNR in short-baseline TDOA systems [13,39] is within the range from −9 dB to 25 dB. So, we set the SNR to begin from −10 dB. The relationship between the performance of TDE algorithms and the SNR is depicted in Figure 6. When the SNR is below 10 dB, the proposed algorithm performs similarly to QLS, sinc-LS, and Hatch filter. The error of the proposed algorithm is primarily limited by the performance of disambiguation. When the SNR is within the range of 10 dB to 14 dB, the RMS TDOA error of the proposed method is much higher than peak phase CRLB, but the performance is still better than the Hatch filter and peak position-based methods. As the SNR still increases, the performance of peak position-based methods gradually plateaus and does not further improve, but the proposed algorithm approaches the CRLB.

After performing the digital down-conversion mentioned in Equation (7), the baseband CCF has an advantage in extracting the peak phase. Using the same simulation parameters in Table 1, we plot a CCF without noise in Figure 7a, where the phase in the main lobe remains constant. This means that there is no need for CCF peak interpolation to extract the peak phase. However, when the carrier frequency is unknown and needs to be estimated, there is an estimation error $\Delta f_{esti}$. The phase varies in the main lobe, as shown in Figure 7b. The CCF peak phase accuracy is related to the CCF peak position accuracy (or TDE accuracy).

Once the emitter is noncooperative, it is meaningful to analyze the effect of $\Delta f_{esti}$. Denote the peak position estimation error as $\Delta {\tau }_{pk}$. The formulation of phase error $\Delta {\theta }_{pk}$ is given in Equation (A2), which is the third exponential item, as

$$\Delta {\theta }_{pk}=2\pi \Delta f_{esti}\Delta {\tau }_{pk}.$$

(43)

Denote the true peak phase as ${\theta }_{pk}$. When using the proposed method, the TDOA error $\Delta \epsilon$ is

$$\Delta \epsilon =\frac{{\theta }_{pk}+\Delta {\theta }_{pk}}{2\pi \left(f_c+\Delta f_{esti}\right)}-\frac{{\theta }_{pk}}{2\pi f_c}.$$

(44)

By substituting Equation (43) with Equation (44), we have

$$\Delta \epsilon =\frac{\Delta f_{esti}}{f_c+\Delta f_{esti}}\left(\Delta {\tau }_{pk}+{\epsilon }_0\right)$$

(45)

where ${\epsilon }_0=-{\theta }_{pk}/\left(2\pi f_c\right)$ is the true TDOA value. From Equation (45), the TDOA error $\Delta \epsilon$ is related to the carrier frequency $f_c$, frequency error $\Delta f_{esti}$, peak position estimation error $\Delta {\tau }_{pk}$, and the true TDOA value ${\epsilon }_0$.

Notice that ${\epsilon }_0$ is related to the geometric configuration. In Figure 8, we label the emitter as green points A and B and TDOA sensors as red points. According to Equation (45), the relationship between $\Delta \epsilon$ and $\Delta f_{esti}$ is shown in Figure 9. When the emitter is at point A, $\Delta \epsilon$ is sensitive to $\Delta f_{esti}$. From Figure 6, the performance can reach sub-picosecond levels when the SNR reaches 20 dB. However, when considering unknown carrier frequency estimation, the frequency estimation accuracy would degrade the final TDOA measurements. When the emitter is at point B, the TDOA error is less sensitive compared to point A. In general, the proposed method is sensitive to the geometry.

5. Experiments

In this section, we present the results from selected experiments that demonstrate the performance of synchronization and TDOA estimation.

5.1. Hardware Implementation

We utilize the system prototype to carry out the demonstration experiments, as shown in Figure 2, including receiver, GNSS receiver, and antennas. The receiver has two down-conversion receiving channels, with CH1 operating at $f_1$ Hz and CH2 operating at $f_2$ Hz. The GNSS receiver provides 1 PPS signal and 100 MHz frequency reference signal, achieving coarse time synchronization on the order of nanoseconds. Subsequently, fine synchronization is achieved through the synchronization link consisting of CH2 and omnidirectional antennas. The synchronization error is observed by periodically transmitting a chirp signal with a bandwidth of 50 MHz at a center frequency of 2.5 GHz.

The signal source, acting as the “portable station” in Figure 2, generates a repeated pulsed chirp signal modulated at a center frequency of 3.6 GHz every 10 ms. The receivers can detect the signal from CH1 and be able to collect the complete pulse every PRT. Detailed parameters are presented in

Table 2. System parameters of passive localization validation systems.

Parameter	Value
Power	35 dBm
Center Frequency	3.6 GHz
Pulse Width	10 μs
Analog Bandwidth	1000 MHz
Sample Rate	4800 MHz

.

5.2. Performance of Synchronization

In this part, we concentrate on the performance of the synchronization link using the proposed method. Figure 10a,b show the time error estimated by the peak position-based method and peak phase-based method between sensor i and j. Both of the two kinds of measurements share the same trend, but the peak phase-based method has less fluctuations. To evaluate the performance of the synchronization approach, we first remove the slow varying trend by calculating differences between adjacent data of time error estimation. Then, we calculate the standard deviation (STD) of the residual error, which we assume is independent with time. So, the performance should be divided by $\surd 2$. The STD of the proposed method is 1.12 ps, as opposed to 8.52 ps by the QLS method.

5.3. TDOA Measurements with Different Chirp Rates

In this experiment, we demonstrate the lower bound performance of the TDOA measured by the GNSS receiver under different chirp rates. The sensor arrangement is the same as in Section 5.2. The results are shown in Figure 11. For the QLS method [22,23], we observed a decrease in STD with an increase in bandwidth. In contrast, for the peak phase measurements, the STD remained stable. Using the QLS method and the proposed method from Section 4, we perform Monte Carlo simulations under the same SNR in the experiment. We find that the STD of the simulation results are much lower than the experiment. Since we only consider the white Gaussian noise in simulation, there must be another main noise source in the system. We hypothesize that the incoherence of oscillators might be the main error source.

To verify the hypothesis, we calculate the STD of instantaneous frequency. According to Equation (15), the STD of instantaneous frequency at each signal pulse should be

$${\sigma }_{LO}^{ij}=k·{\sigma }_{\tau }=\frac{Bw}{T_p}{\sigma }_{\tau }.$$

(46)

It is worth noting that these results are achieved using the same oscillators, indicating that they should be consistent across all trials. Therefore, ${\sigma }_{LO}^{ij}$ should be consistent across the same sensors.

Table 3. Estimated instantaneous frequency STD with different bandwidths in the same pulse width.

Result	Bandwidth
	20 MHz	40 MHz	60 MHz	80 MHz	100 MHz
Instantaneous frequency STD	141.49 Hz	133.232 Hz	147.48 Hz	137.294 Hz	144.88 Hz

displays the TDOA precision for different bandwidths with the same pulse width. We observed a consistent trend in the product of the TDOA STD and chirp rates.

5.4. Instantaneous Frequency Measurement

In this experiment, we measure the instantaneous frequency to further verify the main error source assumption in Section 5.3. We configure the signal source to emit a sine wave with a pulse width of 20 μs at 3.6 GHz. After performing the Hilbert Transform and digital down-conversion, the phase angle difference between the two sensors is identified and is shown in Figure 12a. We utilize the least squares method to compute the frequency based on the angle data of each 10 μs interval with a step of 0.1 μs in every PRT, as shown in Figure 12b,c. Although the GNSS receiver can provide a coarse time synchronization, there remains non-coherence. From Figure 12b, even in the same pulse, the instantaneous frequency can change rapidly, which will cause extra errors in time division synchronization using the peak position-based method. But from Figure 12a, the angle changes less due to the integration operation, which remains coherent in the same pulse. The result is shown in Figure 12d; the STD of the instantaneous frequency is about 140 Hz, which is consistent with that presented in Section 5.3.

5.5. Three-Sensor TDOA Localization

In this experiment, we conducted a three-sensor TDOA localization in Hebei, China. The three sensors were arranged along a straight line at equal intervals of 30 m, as shown in Figure 13. We set up the emitter on a tower with a height of 45 m, which was 500 m away from sensor a. We used a differential GNSS technique to measure the accurate positions of the tower and the sensors.

We performed TDOA localization five times, and used the differential global positioning system (DGPS) positions technique to calibrate the first-trial TDOA measurement. Then, we used the same calibration data to correct the mean error of the other three trials. The corrected TDOA measurements are shown in Figure 14, and we found that there remained a linear trend only in the GNSS synchronization. The STD is shown in

Table 4. The STDs of 5 TDOA measurement trials.

Test	Result	1st Trial	2nd Trial	3rd Trial	4th Trial	5th Trial
sensor a and b	Pk-pos	18.05 ps	23.29 ps	18.19 ps	23.82 ps	17.05 ps
	Pk-pha	1.94 ps	2.05 ps	1.75 ps	1.75 ps	1.66 ps
sensor a and c	Pk-pos	43.56 ps	26.74 ps	34.64 ps	26.15 ps	37.90 ps
	Pk-pha	1.34 ps	1.37 ps	1.22 ps	1.26 ps	1.20 ps

. The localization result is shown in Figure 15, where the accuracy is 0.02%R (circular error probable 50%, CEP50%, proposed method), in comparison with 0.33%R (CEP50%, peak position-based method) and 4%R (CEP50%m, GNSS synchronization only).

6. Conclusions

This paper introduces a short-baseline TDOA location experiment using a wireless phase synchronization method. We propose an integer LS-based time delay estimation method that utilizes carrier phase measurements. The validity of the method is verified by simulation and experiment. The experiment shows that compared with the traditional correlation peak position-based method, the STD of synchronization error is reduced from 8.52 ps to 1.12 ps, and the STD of TDOA estimation accuracy is improved from 17.05 ps to 1.66 ps. In such a case, the CEP accuracy is improved from 0.33%R to 0.02%R.

Ambiguity resolution is the core problem of high-precision delay estimation. On moving platforms, such as UAVs, the difference of TDOA values between adjacent pulses may exceed π. This can lead to phase unwrapping failure and requires thorough experimental verification. To support the application of UAV-based distributed high-precision passive positioning, the research work of high-precision positioning on moving platforms will be further carried out by us and verified by experiments.