Single‐Station Lightning Location Using Azimuth and Time of Arrival of Sferics

Lightning strikes produce electromagnetic waves, now referred to as sferics, in the very low frequency (3 kHz–30 kHz) and the extremely low frequency (3 Hz–3 kHz) bands. Within these frequency bands, the Earth and ionosphere form a waveguide in which sferics propagate long distances with low attenuation. The structure of the received sferic waveform is mainly a function of propagation distance and the waveguide's parameters. This suggests that each observed sferic waveform contains information about the distance that this sferic has propagated which can be used to geolocate lightning. There are various approaches for analyzing received sferics, which mostly rely on measurements from multiple stations. However, in these methods, each station imposes an additional cost for building, maintenance, and synchronization. Here we present a novel method to estimate both the emission time and location of lightning, which works by measuring sferics recorded at a single station. We first process the sferic waveforms to obtain the arrival times of the very low frequency and extremely low frequency radiation components which propagate with different speeds. Once these two separate arrival times are determined, we use them to approximate the distance the sferic propagated in the Earth‐ionosphere waveguide. We have used this novel method in combination with a method to find sferic direction to geolocate a significant number of lightning strikes for 4 July 2013. Using this proposed method, the distance of propagation estimates are accurate to within 6.7% of the National Lightning Detection Network‐determined propagation distance, and the direction of propagation estimates are accurate to within ∼ 1.3% of the National Lightning Detection Network‐determined direction.


Introduction
Electromagnetic waves originating from lightning are called radio atmospherics, or sferics. The very low frequency (VLF) and extremely low frequency (ELF) components of the sferics propagate global distances through the Earth-ionosphere waveguide. A number of studies have used these components individually or together to geolocate lightning from recorded sferic data (Price et al., 2002;Said et al., 2010). Figure 1 shows the VLF and ELF components of a typical sferic waveform.
The VLF component, highlighted in pink on the left, is the initial burst, which is followed by the ELF component, highlighted on the right. In Figure 1, the magnitudes of both components are normalized (the plotted magnitude is M−mean (M) max (M) , where M is the measured time domain signal amplitude). The VLF and ELF components of different sferics vary in time duration, but in general the VLF component is considerably shorter than the ELF component. The data set used in this paper was recorded at two stations, one located at Arrival Height, Antarctica (78 • S, 167 • E), and the other at Sondrestrom, Greenland (66.99 • N,50.95 • W), on 6 June and 4 July 2013, respectively. In each of the stations, sferics were measured by two orthogonal magnetic loop antennas in north-south (NS) and east-west (EW) geomagnetic directions. The data are recorded with a sampling frequency of 100 kHz, and antennas are capable of measuring horizontal magnetic field in frequencies ranging from 300 Hz to 40 kHz.
Using the above data and our proposed method, we first estimated the distances of the lightning strikes and then compared the accuracy of our estimations with the locations estimated by the National Lightning Detection Network's (NLDN) multistation method.

Current Lightning Location Methods
Several prior studies have attempted to utilize observed sferics at multiple ground stations to locate lightning. Some of them are based on multistation detections and use a central processor to employ time of arrival and magnetic direction finding (MDF) techniques (Nag et al., 2015;Said, 2017). NLDN and the Global Lightning Detection Network (GLD360) are examples of multistation networks used to locate lightning. NLDN, operated by Vaisala, Inc., employs a network of low-frequency (LF) and VLF sensors to provide lightning location within the continental United States. Using NLDN, the arrival time and azimuth are measured with accuracies of 1.5 μs and 1 • , respectively (Cummins & Murphy, 2009). The detection efficiency is estimated to be ∼ 60-76% for cloud to ground strokes . The median location error, as an important evaluation metric (Murphy, 2018), was reported around 308 m (Nag et al., 2011). In GLD360, each station correlates observed sferic data with a waveform bank. Multiple stations then send their results to a central processor to approximate the lightning's time and locations. This method achieves 57% flash detection efficiency and ∼2-5 km accuracy (Said et al., 2010(Said et al., , 2013. Alternative methods exist to locate lightning using mostly VLF measurements from a single station. For example, in Byerley et al. (1992) the ratio between radiation field component (Er) and an electrostatic field component (Es) is used for close-range lightning detection. The ratio of electric and magnetic fields or wave impedance also have been used to estimate propagation distance using a single station (Burke & Jones, 1995). However, these methods require the availability of both electric and magnetic field antennas, introducing additional cost and complexity in each station, which also is a concern in multistation locating methods. (Ramachandran et al., 2007) described another method which works based on lightning-generated VLF sferics received in a single station. The reports 8.8% average distance estimation error, but it does not benefit from using ELF components of the sferic to improve distance estimations.
On the other hand, Wait developed a theory proposing that lightning can be located by examining the structure of recorded sferic data from a single location. Wait (1960a) also proposed a theory to estimate propagation distance by modeling the ELF component of the sferics.
However, Wait's method suffers from limited accuracy. An average error of 45.1% for daytime conditions with a standard deviation of 13.8% is reported for his method (Mackay & Fraser-Smith, 2010). Wait assumed a certain analytical form for the lightning current moment, which nonideally can be used to model the ELF component of the sferic. Later more accurate approximation for source current moment was introduced by Cummer and Inan (2000) and used by Mackay and Fraser-Smith (2010) which shows an improvement in lightning distance estimation. However, this improvement is also limited by assumptions made regarding the source and current of the lightning discharge. On the other note, Ogawa has found that by examining the background noise in 1 Hz to 11 kHz, it is possible to observe a secondary waveform caused by the same lightning strikes. The secondary (antipodal) waveform results from the sferic propagating along the opposite direction around the world. With data from the direct and secondary waveforms, lightning location can be estimated (Ogawa & Komatsu, 2007), though in many cases it is hard identify the secondary waveforms . To overcome these problems, we introduce a lightning detection method based on observations from a single station that do not suffer from these limitations.

Evolution of Sferics Time Domain Structure With Distance
Before introducing our method, we begin by inspecting some recorded sferics shown in Figure 2. By evaluating the data in this figure, we can make several key observations that will form the basis of our proposed algorithm. Each row in the figure corresponds to a lightning discharge reported by the NLDN data set. We also used this data set to calculate the propagation distance and to determine whether the sferic propagated during the daytime or nighttime. All of the sferics shown in Figure 2 propagated during the night. The zero point for each row is an NLDN time estimate for a lightning strike. Although the sferic magnitude is generally attenuated over propagation distance, here the magnitudes have been normalized. Therefore, the magnitude attenuation is not shown in Figure 2. This normalization will help us to see the patterns and variations in the waveform structure for different propagation distances, which is a key for our method. There are two key observations to make: As the propagation distances increase, (1) the time interval between the zero time and the sferic start point increases, and (2) the separation time between the VLF component (initial burst) and the ELF component (subsequent component) also increases. In the following sections, we will show that these two time values are both proportional to propagation distance. This, along with phase velocity derivations described in sections 4 and 5, will be the basis of our method to estimate the propagation distances.

Galejs's Model for the Earth-Ionosphere Wave Guide
In this work we use phase velocities derived from the Earth-ionosphere model developed by Galejs. Galejs formulated a model to study wave propagation in the spherical shell between the Earth and ionosphere (Galejs, 1972). The model makes a number of simplifying assumptions. It assumes the Earth and ionosphere boundaries are concentric spheres forming a spherical waveguide. The ionosphere is assumed to be a sharply bounded and homogeneous ionized medium. Galejs further assumes the height and conductivity of the ionosphere remain constant through the propagation path. Figure 3 illustrates this model with two boundaries located at r = a and r = a + h, where a is Earth's radius and h is the ionosphere height. The source is a vertical dipole is located at = 0, r = a + Z s , and the receiving antenna (shown as orthogonal loops) is located at = 0 , r = a + Z o .
Since the ionization level of the ionosphere is different in day and night, the height and conductivity of the ionosphere vary in daytime and nighttime (Wait, 1960a). Considering this fact, in this paper and in the Galejs model, different values are assumed for the height and conductivity of the ionosphere depending on whether sferics propagate during day or night. During the day h is set to 70 km and the conductivity of the ionosphere is assumed to be 10 −6 S/m. At night h is set to 85 km and the conductivity of the ionosphere is assumed to be 10 −5 S/m. The ground's conductivity is assumed 10 −3 S/m for day and night. The solution of radial fields observed at r = a+Z o is expressed in terms of spherical Bessel functions as

10.1029/2018RS006627
and H m ( +0.5) (u) is the Hankel function of kind m and order + 0.5, u = k 0 r and k 0 is the wave number and is equal to (3) 0 and 0 are the permeability and permittivity of free space, respectively. Due to lack of simple representations for the spherical Bessel functions of large order and u, a solution for approximation for thin shell was developed, where h r → 0. The factor r is also approximated by the average value in the shell (r m = a + 0.5h). With this assumption, equation (1) reduces to equation (4) as below: and S is the propagation parameter, which can be interpreted as the ratio between the wave number along the surface of the earth k and the wave number of the free space k 0 . In this paper, S is used to find the velocity of both the VLF and ELF components of sferics. Galejs defined S using modal equation of Transverse Magnetic (TM) modes for the VLF range in equation (6): where n is the T.M. mode number and can be 1 or an arbitrary integer. Galejs shows equation (6) is strictly valid for a conducting ground surface and a perfectly reflective ionosphere. However, these assumptions may idealize the boundaries in the VLF range and the results may slightly diverge from observations. In the ELF range, where only the TEM mode propagates, the waveguide boundaries appear as nearly perfect conductors and equation (6) can be reduced to equation (7): where and Δ e and Δ g represent the normalized impedances for the ionosphere and ground, respectively. i , g , and 0 are the permittivity of the ionosphere, ground, and free space, respectively. i and g are the electrical conductivity of the ionosphere and ground, respectively. Having defined S, Galejs went on to derive the phase velocities. In our method, the phase velocities are essential parameters for estimating lightning distances.
In the next section we present the phase velocities for the VLF and ELF components of sferics using Galejs's model.

Modeling the Waves's Phase Velocities
The phase velocity is inversely related to the real part of the propagation parameter as follows: is the free-space velocity of electromagnetic waves. Since the Galejs model is most accurate at lower frequencies, where the height of the waveguide is smaller than the wavelength (Galejs, 1972), we can confidently apply this model to obtain the propagation parameter for low frequency bands. Figures 4 and 5 show S in ELF and VLF frequencies, respectively. Later we use these graphs along with equation (10) to find the phase velocity for these components of the sferic. Using equation (7), we calculated the real part of propagation parameter in the ELF frequency band for both daytime and nighttime, which are shown in Figure 4.
It follows from equations (7) and (3) that for an increase in frequency, S will decrease. Figure 4 shows this as a downward trend in S over frequency. The first cutoff frequency for the Earth-ionosphere waveguide is 1,580 Hz (Inan, et al., 2000). As shown in Figure 4, at frequencies less 1,580 Hz, S is greater than 1. Since S is greater than 1, according to equation (10), the phase velocity in the ELF range is less than c. However, as the waveguide boundaries become perfect conductors the phase velocity of the ELF waves becomes nearly equal to c. This is because, in perfect conductors, normalized impedances (Δ e and Δ g ) approach 0, causing S to approach 1. Figure 5 shows the propagation parameter of the VLF waves, calculated using equation (6),  with different values for daytime and nighttime.The value reported in this graph is consistent with phase values reported by Wait (1970). In contrast with the ELF range, the propagation parameter in the VLF range shows different behavior. As shown in equation (6), an increase in frequency will result in an increase in S. This accounts for the upward trend observed in Figure 5.
One interesting observation from equation (6): S is imaginary for lower frequencies and, as a result, the waves in that frequency range are severely attenuated. But as the frequency approaches the cutoff frequencies, k 0 = k c = (n−0.5) h and S approaches 0. Thus, according to equation (10), the phase velocity will be infinite at the cutoff frequencies. In the frequencies higher than the cutoff frequencies, S will gradually becomes greater than 1, and when that occurs, the phase velocity becomes less than c.

Methodology
Section 5 described a model how to obtain the phase velocities for the VLF and ELF components of the sferic. These values will be used as inputs in our method to estimate lightning distance. Our methodology can be divided into three main steps: (1) Sferic identifications in the recorded data, (2) distance and time estimation of lightning using ELF and VLF components of the identified sferic, and (3) direction estimation. Finally, we employ all three of these steps to geolocate the lightning. Each of these steps are described in the following sections.

Sferic Identification in Recorded Data
As mentioned before, the method presented here to find the lightning emission time and location is based on characteristics of the associated ELF and VLF components of sferics. Prior to using sferics, the ELF and VLF components must first be identified and isolated in the recorded data. An example of this isolation process is shown in Figure 6. Figure 6a illustrates a typical sferic waveform generated by lightning in approximately 300 Hz to 40 kHz frequency range. Figures 6b and 6c show the isolated VLF and ELF components of the same sferic, after passing it through high-and low-pass filters, respectively. We also eliminated the phase shifts introduced by filters by compensating the group delays.
The waveform shown above was received by a NS directed antenna . However, the closer the sferic's arrival direction is to the EW direction, the less the NS antenna will be able to capture the signal. As a result, to capture a sferic coming from any possible direction, we need to consider recorded data from the antenna directed in EW, the orthogonal direction, as well. Figure 7a shows the sferic displayed in Figure 6 but now includes both NS and EW directions, shown in red and blue, respectively.
To identify a sferic's waveform coming from any direction, we calculate the envelope of the signal (NS 2 + EW 2 ) 0.5 from the recorded data, shown in Figure 7b. To determine the arrival time of the VLF component (t VLF ), we passed the envelope signal through a high pass filter with a cutoff frequency of 1,500 Hz. The high-pass envelope of the corresponding sferic is shown in Figure 7c. Using the high-pass envelope, we determine t VLF to be the global maximum point in that envelope. To find the arrival time of the ELF component (t ELF ), we passed the envelope signal through a low pass filter using the same cutoff frequency as the VLF component. The low pass envelope of the corresponding sferic is shown in Figure 7d. We also defined the first global peak as the arrival time of the sferic.
It should be noted that the detection efficiency of proposed distance estimation method is limited to the percentage of sferics with clear ELF components. There exists a number of factors that limits the percentage of clear sferics. One of these factors is the interference with other recorded sferics which reduces the percentage of clear and usable sferics. Mackay (2012) shows that on average, 70% of the recorded sferics are interference-free. Another limiting factor is noise threshold. For a sferic to be detectable in a particular station, the sferic power must be higher than the noise level threshold in that station. It means the magnitude of peak current of the sferic should be greater than the station noise cutoff and the attenuation of propagation path.

Distance Estimation Algorithm
All of the previous sections have provided us with required inputs for a novel algorithm that is described in the following sections. To start, assume that the lightning discharge happened at t initial at d distance from our receiver. And also assume t VLF and t ELF are the arrival times of the VLF and ELF components, respectively. Also assume that the waves in ELF and VLF bands propagate along the same path to reach the antennas. Given these parameters and according to the time-velocity relationship: Solving equations (11) and (12) leads to where a is a constant and In equation (14), t ELF − t VLF is the separation time between the arrival of waves in the ELF and the VLF bands and can be calculated from the waveform's characteristics as described in section 6.1. As deduced from equations (13) and (14), we can see that only the velocity and arrival times for the ELF and VLF components from a single station are necessary to find emission time and the propagation distance of the lightning strikes. This frees us from the traditional locating methods, which required recording data from multiple stations. Since the profile of the ionosphere differs during daytime and night time, characteristics of a sferic also varies depending on whether it propagated during the day, at night, or across the day-night terminator. The proposed distance estimation method assumes the sferic propagated during the day or at night. When the sferic crosses the day-night terminator, signal propagation should be modeled over nonuniform waveguide which leads to a much more complex problem. As a possible solution, propagation parameters can be averaged to determine attenuation rates and phase velocities of sferics crossing the day-night terminator. A more complex solution involves modeling attenuation as a dependent of propagation distance and time of day (Mackay & Fraser-Smith, 2010).

Additional Distance Estimation Method
In certain circumstances, there is an additional method that can be used along with our method. When these two methods are used together, they yield more accurate distance approximations from recorded data. Ogawa found that with a strong lightning strike it is possible to observe a secondary waveform caused by the same lightning strike (Ogawa & Komatsu, 2007). These antipodal waves propagate in the opposite direction across the globe (Figure 8). The difference in arrival times between the antipodal wave and direct wave can also be used to estimate the distance between antennas and individual lighting occurrences. Figure 8 shows an ELF sferic generated by a single lightning return stroke, which has three noticeable components. The first burst is correlated to the shortest path from the lightning to the receiver; antipodal path generate the second burst and the third burst is the shortest path and a delay from one around-the-world path. Since the antipodal wave needs to travel at least half of Earth's circumference (∼ 20 Mm), the VLF component of the lightning will be significantly attenuated. However, the ELF component will not display the same level of attenuation and, because of that it should generally be visible in the spectrum data radiated by antipodal waves. Considering this fact, the time-velocity relationship for the direct and antipodal waves are shown in equations (16) and (17): where V ELF is the phase velocity of waves in ELF range d is the shortest distance between return stroke and receiver antenna. t Direct and t Antipodal are correlated to direct and antipodal waves, respectively.

Radio Science
10.1029/2018RS006627 Solving equations (16) and (17) for d leads to equation (18): In conclusion, since both equations (14) and (18) estimate propagation distance separately, using both equations and averaging their distance estimations can potentially lead to a more accurate result. While the Ogawa method is also a single-station method for estimating the propagation distances, it has some restrictions that limit it uses when compared to our method. In particular it is limited to cases where lightning strikes are strong enough that attenuated ELF waves are still visible in the recorded data. Another limitation is that this method is restricted to cases where the lightning occurred not too close to the receiver, because in this case the antipodal wave will be highly attenuated and not visible in the received data. The other disadvantage of Ogawa method is that in his equation, the velocity of the waves stays the same for the whole propagation path. However, the antipodal waves must cross the day-night terminator, which requires considering different wave velocities for day and night.
In comparison to a single-station method presented in Ogawa and Komatsu (2007), the proposed method described in section 6.2 does not require observing antipodal peak on the sferic's structure, making it more widely applicable. Also compared to Mackay (2012), another single-station method, this method is less resource intensive and requires less processing cycles and memory usage.

Direction Estimation
Given the estimated distance between lightning strike and receiver, to obtain the location of the lightning we only need to find the sferic's incoming direction. Horizontal magnetic field measurements could be used in finding the direction of lightning discharges (Füllekrug, 2017). To do this, we use the Magnetic Direction Finding (MDF) method for the single station which was introduced by Said et al. (2010). This method also uses data from two orthogonal loop antennas that are positioned in NS and EW directions. The incoming sferic excites the NS and EW antenna, and as shown in Figure 9, inducing signal B 1 (t) and B 2 (t) in NS and EW antennas, respectively.
Since the ratio of B 1 (t) and B 2 (t) is proportional to the tangent of the azimuth, they can be sampled and used to find the incoming direction of the sferic. Figure 10 illustrates an implementation of Said et al. (2010) method using our data. As the method is developed for the VLF component of the sferic, we first isolated the VLF component by passing it through a high pass filter (1,500 Hz). After isolation, we plotted the NS recording against the EW recording for the first 500 μs. The tangent of the azimuth is the slope of the least squares fitting line of the plot. However, if the data recorded from one of the antennas is noisy, the accuracy of this method will dramatically drop. We found out that an improvement can be achieved in direction estimation if we instead consider sferic waveform in full bandwidth and not limit ourselves to only the VLF part of the sferics. This would add more data points to the model and potentially helps fitting a more accurate line. Figure 11 illustrates an our improvement on Said's method on a noisy recorded data. One potential limitation to this method could be the cases where the ELF component has a lower SNR than the VLF component, where considering noisy ELF data, could drop the accuracy of fitting line.
Also since the crossed antennas at the site are never perfectly aligned to the geomagnetic north-south and east-west, we used equation (19) to compute a corrected arrival azimuth, corr , using the calculated azimuth, calc , and the correction factors specific to the recording site.

Results
To estimate the distance of propagation, we employed our proposed method (described by equation (14)) along with the Ogawa method (described by equation (17)) when possible. By using both equations and averaging their distance estimations, we can potentially gain a more accurate result. However, the Ogawa method could only be applied to a limited number of sferics where special conditions mentioned in section 6.3 are satisfied. Method correctness requires that the velocity computations and arrival-time measurements be computed in the same frequency, since both are frequency-dependent. In this work, we have arbitrarily selected 500 Hz and 30 kHz as operating frequencies for the ELF and VLF computations, respectively. To use our proposed method we need to find the constant a (equation (14)). This constant is a function of the ELF and VLF velocities, which are themselves a function of their respective propagation parameters. For the ELF component, we applied = 500 Hz, i = 10 −5 S/m, and we obtained S = 1.05 for the nighttime profile. Substituting the S value in equation (10) gives a phase velocity of ∼ 2.85 × 10 8 m/s. For the VLF component, we used equation (6) where = 30 kHz, i = 10 −5 S/m, and we obtained S = 1.005 for the night. Substituting the S value in equation (6), gives a phase velocity around 2.98 × 10 8 m/s for the night.
Substituting the above phase velocity values in equation (14) gives distance estimations described below in equation (20).
Note that equation (20) only applies for sferics in which height of the ionosphere is constant in their propagation path, that is, those sferics that do not cross the day-night terminator. Figure 12 shows a histogram of the estimated error for lightning distance for a set of 100 sferics that occurred at night. As seen from Figure 12, distance estimation has an average error of ∼ 6.7% with ∼ 5% standard deviation. We further compared our time emission estimations with those of the NLDN for the same strikes. It is worth noting that this average error is in percentage of the actual distance between station and the lightning stroke and would increase as station-lightning stroke distance increases.
The histogram of the time emission approximation error is shown in Figure 13. The result shows ∼ 2×10 −4 % average error and ∼ 3 × 10 −4 % standard deviation.
The coverage range of the method can be estimated using distribution of the lightning's peak currents and also minimum detectable peak current around the globe. Sondrestorm station mostly covers lightnings occurring in the Northern American and Western European regions (Mackay, 2012). The error estimation results for the direction finding method described in section 6.4, are also shown in Figure 14a. This method yields average error ∼ 6.7% with ∼ 17.6% standard deviation. Using proposed improvements in section 6.4, we also computed a histogram of the error for direction approximation for same set of sferics . As seen in Figure 14b the average error dropped by ∼ 5.4% in comparison with the method proposed by Said et al. (2010).

Conclusion
We have introduced a novel method to estimate the propagation distance and the emission time of sferics using recorded data from a single station. Our method is more accurate and more widely applicable than the Ogawa method while being less resource intensive than other single-station methods (Mackay & Fraser-Smith, 2010) as well as multistation methods, like the NLDN.
These improvements were achieved first by making key observations into the time relationship between ELF and VLF components of sferics. With these observations and the Galejs method for approximating the phase velocities we estimated the propagation distance of the sferics. Additional improvements were achieved by estimating the azimuth of the sferic by considering the whole sferic waveform and not just the VLF part.
In short, while our method may be less accurate than multiple station methods, it is more accurate and less resource intensive than previous single-station methods. Given the advantages of using single-station measurements relative to using multiple-station methods such as reduced cost for building, maintenance, and synchronization of the stations, it is a great alternative to multiple-station methods in many applications where the precise location of individual lightning discharges is not the significant factor. Identifying lightning storm regions is the important factor which can be achieved using the proposed method by clustering the measurements of each single lightning within the storm.