Introduction

The infrastructure which underpins modern society is reliant on accurate time synchronisation. Such deep dependence is now prevalent in computer networks, electricity transmission grids, transportation, telecommunication and emergency services. The loss of an accurate time reference represents a high risk event and a known critical threat to societal security1. Radio based time dissemination methods are particularly vulnerable to spoofing and jamming attack vectors. A diversity of GNSS alternatives are required to mitigate the risk of an outage.

On electricity networks, the need for accurate and low latency time synchronisation is well established. A comprehensive study of timestamping requirements for distribution networks emphasises the importance of both accuracy and latency2. Applications which demand hardware timestamping, such as Fast Distributed Generation (DG) disconnection, state estimation and Islanding in Microgrids require time accuracies in the order of 1 μs and latencies of 1–50 ms. Although these applications can be served by GNSS and Precision Time Protocol (PTP) algorithms such as IEEE-1588, the vulnerability of the grid to acute losses of these technologies is reported as a major risk to electrical supplies.

Alternatives for widescale time dissemination over large areas already exist. Notably, the enhanced Long Range Navigation (eLoran) system is recognised as a potentially viable alternative to GNSS for both ___location and timing3,4,5. The system uses a number of beacons to transmit a UTC-synchronised signal which propagates largely as a ground wave, hugging the surface. The receiver calculates its position based on the well known time difference of arrival approach in combination with corrections to account for slower than speed of light propagation speeds over land and sea. Studies have demonstrated UTC-traceable 20 m ___location and ±100 ns timing accuracy in its differential form6, which meets stratum-1 requirements. However, as a wireless system, it remains vulnerable to spoofing and jamming, though due to the higher signal strength, much less so than for GNSS. Since eLoran relies on ground wave propagation, its velocity depends on meteorological conditions7 and a constantly evolving correction for such variations is required. Furthermore, many parts of the world, particularly Europe and North America, require significant new investment in transmission towers, radio equipment and underpinning infrastructure to bring eLoran to service at national scales.

In this paper, the existing National Grid Transmission Network is proposed as a Precision Timing Network (PTN), forming the basis of a non-wireless candidate for a GNSS alternative. It is shown that the propagation characteristics of typical extra-high Voltage (EHV) transmission lines present a low-latency, low-attenuation transmission medium for low energy signals coupled into their aerial modes, a technology already established in Power Line Carrier (PLC) communications. A potential advantage of this approach is the relatively constant speed of propagation which eases the difficulty of calculating the TOF between transmitter and receiver. Furthermore, since the vast majority of signal energy remains closely coupled into the power line itself, this method of time dissemination can propagate over extremely large areas whilst remaining largely immune to wireless spoofing or jamming. Therefore, the use of the existing transmission network infrastructure could provide a viable means of disseminating an accurate timing signal nationally should GNSS fail.

The use of Chirp Spread Spectrum (CSS) has recently been proposed for time dissemination on distribution networks, with reported timing accuracy of 7.5 μs8, making it suitable for most smart grid applications. However, the system is tested on relatively short distribution lines (15–45m) under laboratory conditions. Extending the technique to national scales through transmission networks requires new innovations. For example, the Time of Flight (TOF) between transmitter and receiver - which is in the order of hundreds of μs on national transmission networks, must be corrected. Propagation over longer distances also means higher attenuation and the necessity to operate at lower Signal-to-Noise Ratio (SNR). Finally, new computationally efficient methods are required to improve the resolution of the system with the constraint of limitations of computational power on standard devices like FPGAs or microcontrollers. Recently, a CSS time of arrival estimator based on LoRa modulation was proposed9, which provides a computationally efficient way to perform a huge number of correlations per unit time. This was further refined into a new PLC modulation scheme which is resilient to low SNR regimes and extreme multipath using the principle of averaging10. In the present paper, the use of averaging over repeated chirp symbols will be proposed as an effective way to improve the performance of the system in the low SNR regime.

Extending the application of CSS based time dissemination to transmission networks could extend its scope to national scales. However, the challenge of propagation of a time dissemination signal over such large distances is a problem which needs to be resolved. A key tenet of the proposed scheme is the established practice of Power Line Carrier (PLC), which has been used commercially on Transmission Lines for several decades11,12,13. Modal theory - developed largely as a means to explain the propagation of power line carrier signals - reveals that several low attenuation aerial modes support long range communication on double circuit transmission lines14. Analysis and experimental verification has reported signal attenuation over typical transmission lines of around 0.1 dB/km. Here, these channels are investigated as mediums for time dissemination signals.

In this work, we utilise the idea of chirp based time dissemination as part of a wider methodology for obtaining widescale time synchronisation across arbitrarily large transmission networks. The scheme requires TOF estimation and averaging during GNSS availability to allow a continuation of time synchronisation during GNSS outage, based only on the transmission of a repeated chirp beacon. It is shown through ATP-EMTP simulation and experimentation that sub-μs accuracy can be achieved.

Electrical transmission grids as precision timing networks

Seminal experimental work in the 1960’s and 1970’s verified the existence of low latency and low attenuation aerial modes of propagation on EHV transmission lines15,16,17,18,19. These aerial modes are responsible for long distance communication in PLC systems, which has subsequently been used commercially for several decades in protective relaying applications on transmission networks. The central hypothesis of the present paper is that the aerial modes of transmission lines can be used as the basis of a PTN to disseminate a GPS (or more generally a GNSS) timing signal. In other words, signals can be sent from a GNSS aligned transmitter to receivers on the network, and the arrival time of these signals, coupled with information on the TOF, can be used to reconstruct a precise timing estimate. Later, the use of CSS will be proposed to provide the basic requirements of the system, but first, this section will investigate the propagation modes of a typical transmission line structure.

Figure 1
figure 1

Velocity and attenuation of the L6 transmission line as a function of frequency and mode of propagation.

Full size image

Figure 1 shows details of the L6 transmission tower - typical to the Great Britain Transmission System - to be used in this study including its modal propagation parameters, expressed as the current transformation matrix, \([T_i]\), which has been calculated by the EMTP Lines and Cables Constants (LCC) routine. \([T_i]\) relates the phase ___domain and the modal ___domain according to \([I_{phase}] = [T_i] [I_{mode}]\). It should be noted that the LCC uses a reduction technique to eliminate the earth wire, but this has no bearing on the subsequent analysis. In the modal ___domain, 6 independent modes exist. Figure 1b and c plots their respective velocities and attenuation as a function of frequency. The first major observation is the distinction between the singular ‘ground mode’ and the 5 ‘aerial modes’. The ground mode (mode 1) uses the earth as the return conductor which leads to a significantly higher attenuation and slower velocity than the aerial modes, which use various arrangements of phase conductors for the forward and reverse path of the current. The difference in attenuation between the ground and aerial modes is approximately one order of magnitude which prevents the slower signal from causing ambiguity on the time of arrival estimate. It is shown later that small variations in the velocities of the aerial modes also have negligible impact on the timing estimator for distances typical of a national scale transmission system.

Observe that the modal propagation velocities of the aerial modes asymptotically approach close to the speed of light with little variation for frequencies above few 10’s of kHz. Small variations in the conductor positioning, earth impedance or conductor resistivity have little impact on the velocity as it approaches its asymptote. Attenuation in the aerial modes (modes 2 to 6) are less than 0.1 dB/km at frequencies lower than 1 MHz, which confirms the early experimental work and practical experience of long range PLC.

The presence of low-latency and low-attenuation modes of propagation makes the transmission network a good candidate for a PTN, even over large geographical areas. The advantage of such a scheme is the possibility of widescale time dissemination in the event of a GNSS outage, and the invulnerability of the system to conventional jamming and spoofing techniques which are likely to disrupt wireless GNSS alternatives. In this work, we propose using CSS operating within similar bandwidths to traditional PLC schemes. The use of CSS allows accurate timestamping, advantageous propagation properties and the possibility of averaging over time to improve timestamping accuracy.

We propose a method of constructing a TOF estimate whilst GNSS is operational and available to both the central node and the receivers (the calibration stage). The TOF - the delay between central node and receivers - can be continuously measured by each receiver if (1) The transmitter sends the beacon at exactly the top of the second, according to the rising edge of a 1PPS signal, and (2) Each receiver records the arrival of the beacon relative to the same 1PPS rising edge. Later, we propose the use of CSS and the use of statistical estimation to perform this calculation with continuous refinement over time, and the possibility to adjust for changes in propagation speed. The second step in the process commences the moment GNSS fails (the implementation stage). Now, each receiver loses its 1PPS time source, but the transmitter is able to switch to an accurate atomic clock and therefore continue to transmit a beacon which is aligned to the top of the second. Each receiver, with knowledge of the precise TOF between itself and the transmitter, can retain time synchronisation despite the loss of GNSS.

CSS based approach

The proposed method requires a means to transmit a synchronised beacon from the central node to receivers. In this section, we propose the use of time-averaged CSS as a computationally efficient means of providing accurate timestamping, TOF measurement. First, the CSS approach utilised by the LoRa physical later is described.

Mathematical description of the timestamping algorithm

In recent years, the use of chirp waveforms for time dissemination on power networks has been proposed and investigated. Specifically, the method proposed in8 and further explored in9 is based on cross-correlating a received chirp, transmitted from a time-synchronised source through a distribution network, with a local, clean copy of the chirp. The chirp’s excellent autocorrelation properties facilitate precise and repeatable timestamping. Recent advances have focused on how to make the process more computationally efficient to facilitate a practical implementation with the highest possible resolution. This development is important to the current method and will be addressed in this section.

The autocorrelation function of a chirp is a sinc function with a well-defined peak corresponding to the moment of alignment. A raw implementation of this approach requires a sliding window correlation of the received signal with a local copy of the chirp. The idea to correlate a received chirp signal with a clean, local copy of itself is the core idea underpinning chirp based timestamping. However, this approach requires a full correlation to be performed for each new sample, which is computationally intensive. Since each new sample requires a full correlation, the resolution of the system is set by the sample time, and the computational complexity is, therefore, proportionate to the inverse of the resolution. In this work, the windowed approach used by LoRa to modulate information is used as the starting point. The windowed approach does not yield a new correlation result per sample, but instead yields a set of correlation results per window. The correlation peak, and in particular its offset, can be used to determine the precise time of arrival of the chirp. First, a brief description of the LoRa physical layer will be given. Note that a more extensive description of the LoRa modulation scheme is available in literature20.

Consider a chirp with \(2^{SF}\) samples, where SF is the Spreading Factor. The symbol energy is \(E_s\), T is the LoRa sampling period, \(B=\frac{1}{T}\) is the LoRa bandwidth and \(k\in 0,1,2\dots 2^{SF-1}\) is an index representing the cyclic offset of the chirp. The \(n=0,1,2\dots 2^{SF-1}\) sample digital symbol is given as:

$$\begin{aligned} w_k(nT) = \sqrt{\frac{E_s}{2^{SF}}} e^{j2\pi \cdot (k+n) \mod 2^{SF} \cdot \frac{n}{2^{SF}}} \end{aligned}$$
(1)

\(w_k\) can be used to generate the so called base chirp (\(i=k=0\)), defined as the chirp with a cyclic offset of zero. Alongside the other chirps, it completes the set of basis functions across all possible values of k.

After transmission, \(w_k\) is corrupted by the channel and noise, becoming \(r_k\). Equation (2) shows the correlation process, expressed (in its raw mathematical form) as the inner product of the received symbol, \(r_k\) and the complex conjugate of symbol i, \(w^*_i\). In the presence of Additive White Gaussian Noise (AWGN), \(\phi _i\), the correlator output will yield a maximum of \(\sqrt{E_s} + \phi _i\) when \(i=k\), (assuming k is the transmitted symbol) and \(\phi _i\) elsewhere.

$$\begin{aligned} y(i)= \sum _{n=0}^{2^{SF}-1} r_k(nT) \cdot w^*_i(nT) = {\left\{ \begin{array}{ll} \sqrt{E_s} + \phi _i &{} \quad i=k\\ \phi _i &{} \quad i \ne k\\ \end{array}\right. } \end{aligned}$$
(2)

where \(w^*_i\) is the conjugate of all possible basis functions. When \(i=k\) a correlation peak will occur and an arg max routine can be used to determine the index of the maximum of the correlator output, y. In LoRa, the emergence of a correlation peak reveals the transmitter communication symbol, which is chosen by the transmitter as one out of k cyclically shifted chirps. The method has recently been adapted for Power Line Communications (LoRa-PLC) on networks exhibiting the dual problem of severe multipath interference and low SNR21. Multipath effects are particularly pronounced on power networks and result in multiple, delayed copies of the incident chirp at the receiver. Interestingly, LoRa-PLC demodulates each of the delayed copies separately because they correlate strongly with the progressively time shifted bank of chirps making up the LoRa basis functions. Equation (3) shows this process mathematically. Note that the correlation results for a multipath channel follow the mirror of impulse response of the channel, where index \(\alpha _{1}\dots \alpha _x\) are the magnitudes of the consecutive impulses comprising the impulse response of the channel.

$$\begin{aligned} y(i)= \sum _{n=0}^{2^{SF}-1} r_k(nT) \cdot w^*_i(nT) = {\left\{ \begin{array}{ll} \sqrt{\alpha _1 E_s} + \phi _i &{} \quad i=k\\ \sqrt{\alpha _2 E_s} + \phi _i &{} \quad i=k-1\\ \vdots &{} \quad \vdots \\ \sqrt{\alpha _x E_s} + \phi _i &{} \quad i=k-x\\ \phi _i &{} \quad \text {elsewhere}\\ \end{array}\right. } \end{aligned}$$
(3)

Plotting the correlation results in reverse order mimics the impulse response of the channel. In TOF calculations, the incident wave is of principal importance because it carries the TOF information, and later impulses are disregarded provided they are far enough away. However, the impact of mode separation and multipath effects, where the signal reaches the receiver on a secondary route which is only marginally longer than the incident, might result in competing correlation results and a deterioration in performance.

Improving the computational efficiency using the LoRa approach

The approach described in Eq. (2) requires \(2^{SF}\) complex multiplication and adds for each of the \(2^{SF}\) possible basis functions. A far more computationally efficient method is available based on the following. Consider the initial inner product between the received signal and the \(2^{SF}\) basis functions:

$$\begin{aligned} y(i)= \sum _{n=0}^{2^{SF}-1} r_k(nT) \cdot \sqrt{\frac{E_s}{2^{SF}}} e^{-j2\pi \cdot (i+n) \mod 2^{SF} \cdot \frac{n}{2^{SF}}} = \sum _{n=0}^{2^{SF}-1} r_k(nT) e^{-j2\pi \cdot \frac{n^2}{2^{SF}}}\cdot \sqrt{\frac{E_s}{2^{SF}}} e^{-j2\pi \cdot (i+n) \mod 2^{SF} \cdot \frac{n}{2^{SF}}} \end{aligned}$$
(4)

It is shown in the mathematical analysis of LoRa20 that the second exponential term in Eq. (4) is equivalent to \(e^{-j2\pi q n \frac{1}{2^{SF}}}\) for all combinations of n and k meaning:

$$\begin{aligned} y(i)= \sum _{n=0}^{2^{SF}-1} r_k(nT) e^{-j2\pi \cdot \frac{n^2}{2^{SF}}}\cdot \sqrt{\frac{E_s}{2^{SF}}}e^{-j2\pi i n \frac{1}{2^{SF}}} \end{aligned}$$
(5)

Eq. (5) can be computed in two parts. First, the \(r_k(nT) e^{-j2\pi \cdot \frac{n^2}{2^{SF}}}\) term is equivalent to the multiplication of the received signal, \(r_k\) with a downchirp. This process of ‘dechirping’ is described as:

$$\begin{aligned} d(nT)= r_k(nT) \overbrace{e^{-j2\pi \cdot \frac{n^2}{2^{SF}}}}^{\text {downchirp}} \end{aligned}$$
(6)

Leaving:

$$\begin{aligned} y(i)= \sum _{n=0}^{2^{SF}-1} d(nT) \cdot \sqrt{\frac{E_s}{2^{SF}}}e^{-j2\pi i n \frac{1}{2^{SF}}} \end{aligned}$$
(7)

Equation (7) is equivalent to the Discrete Fourier Transform (DFT) operation on d(nT), which is computed most efficiently as the FFT. It is well known from the analysis of the LoRa physical layer and applications in Frequency-Modulated Continuous-Wave (FMCW) Radar that a single dechirp operation followed by a Fast Fourier Transform (FFT) is mathematically equivalent to performing the full set of correlations, but at a much lower computational complexity. This opens up the possibility of a huge number of correlations per window - and an associated increase of resolution - on standard devices like microcontrollers and FPGAs. Next, three further innovations in the existing LoRa approach are proposed: (1) the use of a GNSS aligned downchirp to automatically reference all correlation results to a synchronised reference (2) the use of a fine offset to improve the resolution of the system and (3) the use of averaging over n historic TOF estimates.

Improving the resolution using a fine offset

In its basic form, the accuracy of this approach is limited to a resolution set by the distance between the bins in the correlation result, which in turn is set by the LoRa sampling period T and, by extension, the LoRa bandwidth \(B = \frac{1}{T}\). The resolution can be improved by decreasing T, which in turn increases B. The resolution can also be improved without the need to increase the bandwidth by repeating the demodulation process s times with a progressive fine timing offset of \(\delta = \frac{T}{s}\). For example, if the LoRa bandwidth is 100 kHz (\(T= 10~ \upmu s\)), setting \(s=100\) increases the resolution to 100 ns without requiring a bandwidth of 10 MHz. The process can be thought of as increasing the number of basis functions by a factor of s with the distinction that the closely spaced chirps are not orthogonal, yet do have the property of revealing the sinc function and the moment of precise alignment. In its non-efficient form, the process is a correlation of the received signal, \(r_k\), with \(s \cdot 2^{SF}\) chirps which are spaced by \(\frac{T}{s}\). \({\textbf{i}} \in 0, 1, 2 \dots s \cdot (2^{{SF}}-1)\) is now the index of the correlation result with s times as many entries and \(\frac{T}{s}\) greater resolution than i.

$$\begin{aligned} c({\textbf{i}})= \sum _{n=0}^{s \cdot \left( 2^{SF}-1\right) } r_k\left( n\frac{T}{s}\right) \cdot w^*_i\left( n\frac{T}{s}\right) \end{aligned}$$
(8)

Since it is not practical to compute Eq. (8) in hardware, the computationally efficient approach splits the operation into s smaller dechirp \(\rightarrow\) FFT processes, each carried out with a time shift, as shown in Eq. (9).

$$\begin{aligned} c(i+h)= \sum \limits _{n=0}^{2^{SF}-1} d(nT + h\delta ) \cdot \sqrt{\frac{E_s}{2^{SF}}}e^{-j2\pi i n \frac{1}{2^{SF}}} \end{aligned}$$
(9)

where h is an integer between 0 and s, representing the fine offset. The fine shift operation is applied to the received waveform, \(r_k(nT + h \delta )\) (see Eq. 6).

GNSS alignment of the downchirp

Another key innovation in the proposed algorithm is the use of a GNSS aligned downchirp. Equation 10 shows that the d(nT) term is GNSS aligned, which means the downchirp in Eq. (6) is aligned to GNSS:

$$\begin{aligned} c(i+h)= \sum \limits _{n=0}^{2^{SF}-1} \overbrace{d(nT + h\delta )}^{\text {GNSS Aligned}} \cdot \sqrt{\frac{E_s}{2^{SF}}}e^{-j2\pi i n \frac{1}{2^{SF}}} \end{aligned}$$
(10)

An important aspect of this approach is ensure that the downchirp is as precisely aligned to GNSS as possible, meaning the first sample is concurrent with the rising edge of the 1PPS timing signal. The generation of the chirp at the central node (transmitter) is also GNSS aligned. This means the correlation result, c, will maximise at a distance from the rightmost index which is proportionate to the TOF of the chirp from the central node to the receiver. An overview of the proposed method to obtain the TOF during the active stage of GNSS (calibration) is shown graphically in Fig. 2. An FPGA implementation of this approach is shown in the "Experimental methodology" section.

Figure 2
figure 2

Obtaining the TOF measurement during the calibration stage \({\textcircled {1}}\) The Central Node has access to a GNSS device, and a backup atomic clock, \({\textcircled {2}}\) A series of chirps are generated and aligned to 1PPS. \({\textcircled {3}}\) The chirp signal is coupled to the transmission line using established PLCC capacitive coupling. \({\textcircled {4}}\) The signal propagates through the network, predominantly on the aerial modes. \({\textcircled {5}}\) Each substation receives the chirps depending on the propagation delay. \({\textcircled {6}}\) In the calibration stage, the receiver also has access to GNSS 1PPS, and the received chirp has an offset proportionate to its propagation delay. \({\textcircled {7}}\) To guarantee this, the downchirp used in the dechirp process is aligned to \({\textcircled {8}}\) The FFT is performed, which reveals, \({\textcircled {9}}\), c correlation results (after s fine tune operations). \({\textcircled {10}}\) The distance of the maximum correlation result from the reference line reveals the TOF. \({\textcircled {11}}\) a moving average further improves the accuracy.

Full size image

Use of averaging to further improve accuracy for low SNRs

Another important aspect of the proposed algorithm is the use of averaging to improve the timing estimate in the low SNR regime. Equation (11) shows how averaging over n results in \({\overline{TOF}}\), which has a significantly narrower distribution than each individual TOF estimation.

$$\begin{aligned} {\overline{TOF}} = \frac{1}{n} \sum \limits _{j=1}^{n} \arg \max _{c} \end{aligned}$$
(11)

Application of CSS to TOF measurement and time dissemination

Treated as a black box, the timestamping algorithm summarised in Fig. 2 will give the highest spike (\(\arg \max _{c}\)) in the correlation results window (c) at a distance to the left of the ‘reference line’ which is directly proportionate to the TOF, assuming that both transmitter and receiver use GNSS aligned chirps (which is guaranteed to be the case in the calibration stage). In other words, time delayed versions of the same chirp manifest in the demodulated correlation result at progressively lower positions. Delayed multipath energy reveals itself to the left of the reference line. This property is exploited in two ways to implement the proposed method. First, in the calibration stage, both the central node (transmitter) and receiver use local dechirp operations which are synchronised to a shared clock (e.g. GPS). The position of the right most spike in the correlation result will reveal the TOF based on its distance from the base chirp position (at index 0). Later spikes (further towards the left) can be disregarded for TOF calculations. Given that c is the demodulation index, the maximum value of this index is defined as \({\mathcal {D}} = \arg \max _{c}\).

Figure 3
figure 3

Flowchart describing the proposed method.

Full size image

A flowchart summarising the proposed method is shown in Fig. 3. During the calibration stage, this equation yields the TOF as \({\mathcal {D}} = TOF = \arg \max _{c} \mid \text {GNSS active}\), where TOF is expressed in terms of the number of fine shift separations from the reference line, and, for simplicity, the reference line is designated the zero index, and the index becomes increasingly positive as it moves to the left. Conversion to time requires a multiplication by the time per fine shift, e.g. the resolution. Each demodulated chirp results in a new TOF observation, and many such observations can be made within the calibration stage. To minimise the impact of noise, an estimator can be formed based on a moving average. Second, when GNSS is lost, the central node switches to an atomic clock and continues to send a time aligned chirp. The receivers, without a local source of time synchronisation, use the time of the received chirp, and previous knowledge of the TOF, to reconstruct a timing estimate over n history terms, as shown in Eq. (11).

Figure 4
figure 4

Evolution of the demodulation index during the transition between the calibration and implementation stages.

Full size image

Figure 4 shows the evolution of \({\mathcal {D}}\) in the transition between the calibration and implementation stages. In the calibration stage, \({\mathcal {D}} = TOF = \arg \max _{c}\) and is relatively constant (disregarding the error terms), and the moving average is used to build up a continuously improved TOF estimate over time. The offset shows the difference between the 1PPS time sources at the central node and the receiver, which is close to zero during the calibration stage when both locations have access to GNSS. When GNSS is lost, the implementation stage begins and the offset becomes nonzero. The slope of the offset curve will depend on the drift of the local receiver clock relative to the central node’s atomic clock. As the receiver clock drifts away from the atomic clock of the central node, the offset increases. However, the demodulation index, \({\mathcal {D}}\), shifts as the exact inverse of this offset. Therefore, a timing offset, \({\mathcal {T}}\), can be formed as \({\mathcal {T}} = {\mathcal {D}} - {\overline{TOF}}\). \({\mathcal {T}}\) in the implementation stage contains a correction for the TOF, and therefore provides full knowledge of the offset between the time of the receiver downchirp (used in the dechirp operation) and the time of the transmitted chirp. The disciplined 1PPS timepulse is therefore obtained by counting forward the remaining fraction of the second: \(\text {1PPS}_{corrected} = \text {1PPS}_{local} + (1-{\mathcal {T}})\). Note that this assumes a chirp of 1 second duration. In practice, each chirp duration will be 1-100 ms. A reconstructed 1PPS time pulse can still be reconstructed if an exact integer number of chirps fits within 1 second.

Simulation and experimental methods

Noise modelling

Modelling the noise on an EHV or UHV line is challenging. Much theoretical and empirical work has been carried out to characterise the noise environment on low and medium voltage networks, revealing a mix of AWGN, coloured and impulsive noise22. At transmission level, corona noise must be added to this list23,24. In this work, we propose the use of a general noise model based on the \(\alpha\)-stable distribution, S(\(\alpha\),\(\beta\),\(\gamma\),\(\delta\))25. The Gaussian distribution is a special case of the \(\alpha\)-stable distribution with the first shape parameter \(\alpha =2\), the second shape parameter \(\beta =0\), the scale parameter \(\gamma =\frac{ \sigma }{\sqrt{2}}\) and the ___location parameter \(\delta = \mu\), where \(\sigma\) and \(\mu\) are the standard deviation and mean of the resulting Gaussian distribution. However, by setting \(\alpha <2\), the distribution exhibits heavier tails which model the impulsiveness of the power line channel.

Any communication scheme operating on the power line should be able to cope with high levels of impulsive noise against the general background of AWGN noise. A widely used mitigation technique is known as thresholding22. We propose the use of a thresholding technique based on clipping, similar to that recently utilised to improve the performance of Chirp-PLC in the presence of impulsive noise10,26. Below, it is shown that the \(N_{th}\) noise sample, after thresholding, is either the raw sample from the \(\alpha\)-stable distribution (if the value is less than the threshold, T) or T itself if the same exceeds T.

$$\begin{aligned} N_{th}= {\left\{ \begin{array}{ll} N_{\alpha }(n), &{} \quad \text {if}~ |N_{\alpha }(n)|<T \\ T \cdot sgn(N_{\alpha }(n)), &{} \quad \text {otherwise} \end{array}\right. } \end{aligned}$$
(12)

The threshold, T, is set to a user-defined multiple of the 90\(^{th}\) percentile of the absolute value of the measured noise samples, \(N_{90}\). Therefore, \(T = N_{90}\cdot M\), where M is the multiple. Since the \(\alpha\)-stable distribution does not exhibit finite variance the standard definition of SNR is invalid. Instead, the signal-to-dispersion ratio is used in this work, as expressed as \(SNR_D = 10\log _{10}\frac{P}{2\gamma ^2}\), where P is the signal power.

Simulation methodology

The simulation methodology, outlined in Fig. 5, utilises JMarti frequency dependent line models within the ATP-EMTP27, modelling the L6 transmission lines discussed earlier. The transmit chirp is generated within a parameterised C++ foreign model which is precompiled as a native model into the EMTP executable, with input arguments to determine the start frequency, end frequency and number of samples of the base chirp. A consecutive series of chirps from this model is coupled via TACS into the transmission network via a coupling capacitive with properties similar to that used in power line carrier systems. The ATP-EMTP model is called via DOS command from within the main Matlab script, meaning the entire methodology can be automated. Simulation results, particularly the signals from the substations are saved in the native .PL4 format and converted to the .MAT format using the PL42MAT software, which is called from within the Matlab script. The demodulation code is organised into objects, where each object executes the demodulation process with its own configurable set of parameters (bandwidth, spreading factor etc). Each object has methods for (1) Downconversion of the chirp to baseband, (2) Dechirping, (3) Running the FFT, and (4) Finding the index of the max result from the FFT. An inner loop introduces a progressive fine shift and repetition of the demodulation process s times per window, as outlined earlier. Subsequent code pieces together the \(s \cdot 2^{SF}\) correlation results into a single plot, and a final find max routine of this waveform is used to determine TOF estimate based on its distance from the reference chirp (e.g. the rightmost part of the waveform). The automated Matlab script allows Monte-Carlo simulations to be performed for arbitrary noise scenarios (varying SNR and impulsiveness, \(\alpha\)) across a user-defined number of runs. Note that a summary of the simulation and experimental parameters used in this study are shown in Table 1.

Figure 5
figure 5

Simulation methodology procedure (left). Transmission network under test (right).

Full size image
Table 1 Configuration Parameters of Experimental and Simulation Tests.
Full size table

Experimental methodology

Figure 6 shows the architecture of the developed prototype. The system incorporates a u-blox ZED-F9T-00B timing module, which is a multi-band GNSS receiver with a 1-sigma timing accuracy of 5 ns. The module provides two GNSS disciplined timing signals, which is important to the developed technique (further explained later). The CSS based algorithm is implemented in FPGA hardware. The transmitter consists of a Phase Locked Loop which divides one of the GNSS disciplined timing signals by a factor of two. In the configuration shown, the fast timing signal from the timing module, Fclk, is set to 10.485760 MHz. The choice of frequency is to ensure that an integer number of 32 chirp waveforms fit within one 1 second cycle of the slower GNSS disciplined timing signal. A sequence of 32 progressively shifted chirps represents one full cycle from which the receiver can determine the accurate time of arrival to within a resolution of Fclk. A 15-bit counter operating at a rate of Fclk/reads samples from a pre-allocated 32,768 point RAM, which contains a full upconverted base chirp. The shift logic introduces a fine shift after every chirp. As noted previously, this is a computationally efficient way to increase the resolution of the system by a factor of s, where s in this case is 32. A CIC filter provides interpolation prior to digital to analog conversion at a rate of 65 MSPS. The role of the receiver is to estimate the time of arrival of the chirp relative to the 1 Hz GNSS timing signal. Since each of the 32 chirps in a full cycle are progressively shifted, the correlation results will vary, reaching its peak when the shift happens to be closest in alignment to the receivers own chirp. Therefore, it is important that the receiver chirp is precisely aligned with the transmitter’s base chirp. During the calibration stage, each receiver will have access to a GNSS clock, which guarantees this alignment. We define the Sclk - the inverse of which is the LoRa bandwidth - with a clock which is exactly 32 times slower than Fclk/2. This slower clock is used to read from a pair of 10-bit RAMs to produce the real and imaginary parts of the baseband base chirp.

Figure 6
figure 6

Prototype system.

Full size image

An alternative method of demodulation can be used if the incoming signal (after ADC) is immediately sent to Matlab from a Signal Tap logic analyser MEX function which is triggered from the local receiver’s 1PPS clock (which will drift following the loss of GNSS). The full demodulation routine can be performed in software. This method is slower than the hardware based approach but does provide satisfactory performance in the experimental conditions outlined in this study.

Tests are performed using a 700 m length of RG58 coaxial cable with a known propagation delay of \(\approx 3.7 \upmu s\). In the calibration stage, both the central node and the receiver have access to a separate GPS timing module with good access to the sky. Half way through the experiment, the receiver’s antenna is disconnected and the 1PPS timing accuracy starts to drift. To assess the accuracy of the timing estimate in the implementation stage, the offset between the transmitter’s GPS 1PPS timepulse (the true time) and the unsynchronised receiver’s 1PPS is compared. The index of the demodulation result should drift in proportion to this offset.

Results

Effect of mode separation

The slight difference in velocities of the propagation modes could lead to ambiguity in the timing estimate. To investigate this, the received signal on Phase A, \([I_{phase}]\), is decomposed into its modal components, \([I_{mode}]\), using \([I_{mode}] = [T_i]^{-1} [I_{phase}]\), where \([T_i]\) is the current transformation matrix as calculated by the EMTP LCC routine, as outlined in Fig. 1. Figure 7 shows the correlation results for all 6 modes at distances of 1 km, 20 km and 200 km, respectively. At 1 km (Fig. 7a), aerial modes 5 and 6 dominate but there are moderate contributions from the other 3 aerial modes (2, 3 and 4). The ground mode (1) makes a small contribution. At 10 km (Fig. 7b), the ground mode has largely dissipated. At 200 km (Fig. 7c), mode 2 is largely attenuated. As the distance increases, the signal is increasingly dominated by modes 5 and 6. It is important to note that the exact composition of modes will depend initially on the coupling arrangement, but as the signal propagates it is increasingly dominated by its two strongest aerial modes. This is an important result because the propagation speeds of modes 5 and 6 are similar, and the potential of dispersion due to the slower aerial modes is mitigated by their higher attenuations.

Figure 7
figure 7

Normalised correlation results (c) after modal decomposition of Phase A. The ground mode (1) quickly dissipates. Mode 2 is largely irrelevant at 100 km. At 200 km, aerial modes 3 to 6 remain with modes 5 and 6 dominating.

Full size image

Simulated performance of the proposed method

Figure 8 shows histograms of the TOF estimate at all substations on the test network, giving representative performance during the calibration stage of the algorithm. At all substations, the histograms broaden with worsening SNR, but crucially remain centred on the true TOF, offering improved accuracy with estimation of the mean (\({\overline{TOF}}\)). There is no evidence of significant mode mixing, which would result in delayed peaks in the histograms due to the later arrivals of slower modes. We also see no evidence of multipath interference due to signal energy taking different paths to the destination substation.

Figure 8
figure 8

Simulated performance of the TOF Estimate at the 8 substations on the test network (calibration stage). Each histogram is for \(\hbox {n}=1000\,\hbox {runs}\).

Full size image

Figure 9 shows the performance of the timing estimator during the implementation stage . We vary the SNR, degree of impulsiveness (\(\alpha\)) and the clipping threshold, M. 1000 Monte Carlo simulations are performed per condition, allocating 900 runs to the calibration stage and 100 runs to the implementation stage, where the latter experiences a progressive offset between the transmit and receiver clock phase. The 900 runs in the calibration stage are used to generate an estimate of \({\overline{TOF}}\). A timing estimate is recovered in the implementation stage according to the method outlined previously. We observe sub-μs timing accuracy for SNRs of -10dB and better. At \(\hbox {SNR}=0 \,\hbox {dB}\), the standard deviation of the timing estimator is \(\sigma = 250~ ns\). We also observe no deterioration in performance under progressively impulsive conditions (\(\alpha =1.8, 1.6\)).

Figure 9
figure 9

Simulated performance of the Timing Estimate for a calibration moving average of 900 chirps. The bars show n=100 chirp observations in the implementation stage. The boxplots’ whiskers represent \(\pm 1\) standard deviation from the mean.

Full size image

Experimental performance of the proposed method

Figure 10 shows the performance of the timing estimator for the experimental testing. Note that only AWGN is applied in this case, with SNRs ranging from 0 to − 20 dB. Here, we perform calibration for 2000 chirps and implementation for an additional 2000 chirps, sampling a chirp once per second (total time of test of 33.33 minutes). We observe sub-μs timing accuracy for all cases, with a generally worsening performance with increasing SNR. AT SNR = 0 dB, timing accuracy of 200 ns is observed, which approximately concurs with the results from the simulation study. Similarly, at the worst reported case of SNR = − 20 dB, the timing accuracy of 600 ns is only slightly better than the simulated performance on a large transmission network.

Figure 10
figure 10

Experimental performance of the timing estimator for a calibration and implementation stage of 2000 chirps, respectively. The boxplots’ whiskers represent ± 1 standard deviation from the mean.

Full size image

Discussion

The need for viable alternatives to GNSS is widely reported. The value of a national scale alternative which does not rely on wireless signals could provide a compelling backup to conventional satellite and other radio based timing services. This work has demonstrated the possibility of using a national scale transmission grid for time dissemination using a combination of established PLC technology and emerging CSS modulation techniques. We report sub-μs timing accuracy for performance on a simulated large transmission network and experimentally on a point to point coaxial cable of length 700 m. The use of CSS allows the method to operate at low SNRs - down to − 20 dB, which opens up the possibility of widescale implementation over large geographical areas. The ability to use a single Stratum 1 atomic clock source to synchronise the national grid through the existing infrastructure offers an independent means of satisfying increasingly stringent timing requirements independently of other infrastructure.

Although these results are promising, more work is required to address key limitations. An important assumption of the method is the existence of a relatively constant speed of propagation through the network. It is shown that the velocities of the aerial modes of propagation - responsible for long distance communication - asymptotically converge to the speed of light at frequencies above 100 kHz. Indeed, our simulation results on typical transmission networks show a fast separation of the slower modes and the emergence of a dominant aerial mode travelling at a constant velocity approaching the speed of light. However, real life conditions like weather variation, changes in soil conditions and corona noise present a temporal set of variables which might change the propagation delays - a phenomenon well known in the ground wave propagation of eLoran signals7 . Therefore, a key question is: how predictable and stable is the propagation velocity on transmission networks? In much the same way GPS corrects for atmospheric delay terms, the proposed approach could potentially address this issue with correction terms and predictive modelling, but further investigation is required.

Further work is required to establish the performance of the proposed method on a real transmission network. The use of existing power line carrier infrastructure, in particular the coupling apparatus, could provide a practical means of conducting such a trial with minimal cost.

Conclusion

  • This paper proposes the use of electrical transmission networks as PTNs capable of disseminating a sub-μs timing signal on national scales.

  • A computationally efficient method of time-of-arrival estimation based on CSS is proposed. In the calibration stage - where GNSS is available - the TOF between a central node and all other substations can be calculated, and continuously improved based on averaging. When GNSS is lost, it is shown that the knowledge of TOF can be used to derive a timing estimate even if the local clock loses synchronisation to UTC.

  • The LoRa physical layer is used as the basis for the proposed method, with the addition of three key additions - first, the use of a fine timing offset to improve the resolution of the timestamping algorithm beyond the LoRa bandwidth. Second, the use of a GNSS aligned downchirp at the receiver ensures that all correlation results are referenced to a synchronised time (e.g. the GNSS 1 PPS). Third, the use of averaging over n historic observations has been shown to significantly improve the TOF estimate at low SNRs.

  • Simulated performance on a large transmission network demonstrate sub-μs performance at SNRs of − 20 dB or better.

  • The algorithm, incorporating the fine offset and GNSS aligned downchirp, are realised in hardware using an FPGA. Experimental results using a 700 m coaxial cable show broadly similar performance to simulation, with sub-μs timing accuracies reported at SNRs of − 20 dB.