AU749713B2 - Improvements relating to the location of lightning discharges - Google Patents

Description

-2- FIELD OF THE INVENTION This invention relates to determining the location of lightning discharges. In particular the invention relates to utilising an electromagnetic emission from a lightning discharge to determine the location where the lightning discharge occurred. Mathematical calculations are performed to determine the lightning location from a plurality of waveforms each of which represent the emission received at a respective receiver station.

BACKGROUND TO THE INVENTION At present there are systems for determining the location of a lightning discharge using a network of receivers which provide coverage of the area which needs to be monitored for lightning strikes. The receivers detect electromagnetic emissions from the lightning discharges and use the emissions to locate each strike.

Existing lightning location systems using networks of radio receivers for commercial use are of two main classes: those using magnetic direction finding from each receiving station and those based on precision timing of the arrival of the lightning emissions received at each station.

The second class can be further divided into two sub classes. In the first, one measures the time of arrival (TOA) of the very beginning of the lightning pulse at each of the receivers comprising the network. This uses only the first few microseconds of the pulse, which S.correspond to the highest frequencies received (a few MHz), to preserve the sharpness of the pulse and to avoid the sky wave (reflected off the ionosphere) which arrives slightly later.

Since the lightning pulse is dominated by the return stroke (ground to cloud), this system allows location of the ground point of the lightning strike to within a few hundred metres.

Such location is important for insurance inspectors checking claims and for power line companies locating line faults. The system is less able to locate the position of the source (cloud end of the lightning) since this can be laterally displaced from the ground point by as much as 5 kmn. The system also needs a tight network of ground stations with separations of a few hundred kilometres since it relies on ground wave propagation which has high attenuation at the high frequencies required. Such a tight network requires about -3- 100 times the ground station density of the next subclass to be described and so is not commercially feasible for large areas of low population density or over oceans.

The second subclass uses only low frequencies (long waves) which can travel several thousands of kilometres with little loss. This involves many reflections off the ionosphere and the ground so the whole wave form representing the characteristic emission received at any one station is made up of components which have made different numbers of such reflections (hops). The component which arrives first (no hops) may be much weaker than later components, and even undetectable. Thus the time of arrival (TOA) is not clearly defined.

Consequently, in the second subclass, the whole wave form representing the characteristic S. emission received at a pair of stations is compared to measure the arrival time difference (ATD) between the emission received at each station of the pair. There are two ways to measure the ATD: by cross-correlation and wide band interferometry (WBI) which is also known as broad band interferometry (BBI).

One such network using the ATD method is described in Lee, A. C An experimental study of the remote location of lightning flashes using a VLF arrival time difference technique, Q. J. R. Meteorol. Soc., 112, 203-229, (1986a), and Lee, A. C. An operational :0°0i system for the remote location of lightning flashes using a VLF arrival time difference technique, J. Atmos. Ocean. Tech., 3, 630-642, (1986b). The method described uses crosscorrelation to find the ATD. In this method, each receiver of the network is configured to form a plurality of pairs with every other receiver of the network. The waveforms 25 representing lightning emissions received at each of a pair of receivers are cross-correlated.

The resultant waveform contains a peak which is used to calculate the difference in time of arrival of the waveforms at each receiver of that pair. This method is repeated for the other of pairs of receivers in the network to produce a number of ATDs. The ATDs can then be used to determine the lightning location.

WBI involves computing the phase difference between each respective frequency component of the representative waveforms received at each station of a pair. These phase differences can then be used to determine the lightning strike location. As a technique, WBI has not been used to determine ATDs between receiver sites, but rather has only been used at VHF/UHF (very high frequencies/ultra high frequencies) for measuring angles of arrival r, W -4of lightning pulses coming from short segments of a nearby (few tens of kilometres distant) lightning strike. Two pairs of antennas at right angles are used to measure both the azimuth and elevation angles of the strike. A second two pairs, displaced some 10-20 km horizontally, are then used to measure simultaneously the azimuth and elevation angles at a second site. The 3-D position of each lightning segment is then found by trigonometry (one baseline, two azimuth and two elevation angles). From the successive positions in time of these segments, the 3-D shape of the lightning path from cloud to ground is determined such a method is described in Shao, X M, D. N. Holden and C. T. Rhodes, Broad band radio interferometryfor lightning observations, Geophysical Research Letters, 23, 1917-1920, (1996). Thus WBI is used above for angle measurement at microwave frequencies in a way similar to magnetic direction finding at VLF.

SUMMARY OF THE INVENTION It would be advantageous to have a system which utilises WBI at lower frequencies. Such a system could have lower density of receiver stations due to the long ranges allowed by the low attenuation of VLF (very low frequency) propagation in the Earth-ionosphere waveguide and also by the fact that most of the radio energy of lightning is in the VLF band. The receivers in such a network could be approximately ten times further apart (100 times lower density) than that required by some other methods. This would be particularly •suitable for covering wide areas of oceans or sparsely populated regions, where the cost and technical feasibility of a high density network may be impractical.

Therefore, it is an object of the present invention to provide a method and system for determining the location of a lightning strike using WBI at VLF using a receiver network, or at least to provide an alternative to existing methods and systems for locating lightning.

In general terms the invention involves receiving a lightning discharge emission at various locations and processing the received emission to determine arrival time differences of the emission at the various locations. The arrival time differences can then be used to calculate the source of the lightning.

In one aspect the present invention may be said to consist in a method for determining the location of a lightning discharge including: receiving a discharge emission at a plurality of receiver stations, determining a time of detection of the emission at each station, for at least three pairs of stations, comparing the emission received at one station of the pair with the emission received at another station of the pair to determine, using wide band interferometry (WBI) at VLF, the arrival time difference between the emission received at each station of the pair based on the comparison and the times of detection, and determining the location of the discharge using the arrival time differences determined for each station pair.

In another aspect the present invention may be said to consist in a method for determining the location of a lightning discharge including: receiving a discharge emission at a plurality of receiver stations, determining a time of detection of the emission at each station, for each station, comparing the emission with a reference waveform to determine an individual residual time using WBI at VLF, for at least three pairs of stations, determining a residual time difference using the individual residual times for each station of the pair, for at least three pairs of stations determining an arrival time difference for each station pair using the residual time and times of detection, and determining the lightning location using the arrival time difference for each pair.

•0 0 In another aspect the present invention may be said to consist in a system for determining the location of a lightning discharge including: a plurality of stations for receiving an emission from the discharge, the stations configured as pairs, one or more time stamping means for recording the time of detection of the emission at each station, and a signal processor associated with each station pair for implementing WBI at VLF on the emission received at each station of a station pair to determine a residual time difference for that station pair, wherein the processor is further adapted to calculate an arrival time difference .ii :of the emission at each station of an associated pair using the times of detection and the residual time difference, and is further adapted to determine the lightning location using the *arrival time differences calculated for each pair.

S-In another aspect the present invention may be said to consist in a receiver station for use in a system which determines the location of a lightening discharge from an emission, the station associated with another receiver station to form a station pair, the station including: time stamping means for recording the time of detection of the emission received at the station, and a signal processor for determining a residual time difference for the station pair from the received emission by implementing WBI at VLF, wherein the processor is adapted to calculate an arrival time difference of the emission at each station of the pair using the time of detection and the residual time difference.

-6- Preferably the representative residual time difference is calculated using the rate of change with respect to frequency of the phase difference between respective frequency components of the stored waveforms corresponding to the emission detected at the receiver station and second station. Alternatively the emission detected at each station may be compared in a similar manner with a reference waveform to produce a individual residual time difference for each station BRIEF DESCRIPTION OF FIGURES Preferred embodiments of the invention will be described with reference to the following drawings of which, Figure 1 is a schematic of a receiver network according to the invention, Figure 2 shows the network of widely separated receiver stations together with the corresponding contours of location errors when using the VLF waveforms of the sferics according to the invention, Figure 3 shows a flow diagram of the preferred method for locating lightning according to the invention, Figure 4 shows the waveforms of the emission from the same lightning discharge received at two different sites during a test of the invention, Figure 5a is a graph showing the difference in phase with respect to frequency between the two waveforms shown in figure 4, Figure 5b is a graph showing the difference in phase with respect to .frequency between the waveform received at each station and a Dirac delta function, 25 Figure 6 is a plot of 1 kHz averages of the rates of change with respect to frequency of the phase difference plotted in figure 5, to identify and reject averages degraded by spurious effects, Figure 7a shows the waveforms of a subset of 11 frequency components of a simulated pair of simultaneous and identical sferics to illustrate the precision to which the time difference (zero in this case) can be found in the absence of noise, Figure 7b shows the waveforms of the same subset of 11 frequency components as in figure 7a to illustrate the apparent degrading of time precision in the presence of noise 10dB weaker than the simulated sferics, (It l -7- Figure 8 shows a plot of phase differences versus frequency between two simulated sferics with a signal to noise ratio of 10dB, illustrating that for 1000 frequency components the slope of the straight line fit is substantially unaffected by random noise, and Figure 9 shows a block diagram of hardware associated with a receiver station according to the invention.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS Referring to the drawings it will be appreciated that a method and system according to the invention for determining the location of a lightning strike may be implemented in various forms. The following embodiments are given by way of example only.

o In a preferred embodiment the lightning strike or discharge is located by using the phase "'"information of the emitted radio waves over a wide spectrum from about 3 kHz to 25 kHz.

A network of four or more sites with radio receiver stations is used to measure the difference in arrival time of the waveforms received at each pair of sites from the rate of change of the phase difference of these waveforms with respect to frequency.

Preferably at each receiver, the wide band of radio waves (3 kHz to 25 kHz) is continuously sampled at 50 kHz and a representative waveform is stored in a buffer of about milliseconds long. It will be appreciated, however, that this is by way of example only, and that the range of VLF frequencies which are sampled may vary. For example, in an alternative embodiment the band of radio waves from 3-30 kHz could be sampled, and in a test of the system, a band with an upper frequency of 15 kHz was sampled. Preferably the 25 stored waveform represents all, or a substantial part of a VLF portion of the radio waves emitted from the discharge. Once a lightning discharge occurs a characteristic radio emission in the form of a lightning pulse train propagates from the lightning source, and is detected by the receiver. On detection of the lightning pulse train (called a "sferic"), the sampled sferic from the buffer, beginning just before detection, is kept for comparison (usually via a reference or standard waveform) with all the sampled sferics from the same lightning strike at all the other receiver sites which detected it.

Figure 1 shows a preferred embodiment of a network of four or more receiver stations used for the present invention. Each station 10, 11, 12, 13 can detect and receive the radio waves emitted 16 from a lightning discharge 15. The stations are all communicably linked, -8preferably via the Intemrnet 14, although other communication networks could be used. The stations 10, 11, 12, 13 are located to provide coverage of the area to be monitored for lightning discharges and can be spaced in the order of a few thousand kilometres apart. As the stations are all communicably linked, each station can form a pair with any of the other stations. Preferably the network of receivers will utilise every pair combination of receivers in the network to determine the lightning location. Therefore in a network comprising four receivers such as that shown in figure 1, six receiver pair combinations are possible, ie receiver pairs 10-11, 10-12, 10-13, 11-12, 11-13 and 12-13. Only three of these pairs are independent because the second three pairs can be obtained from the first three. For example, 11-12=10-12 (10-11). Note that for this pair, 11-12, the waveform of 10 serves only as an intermediary. This means that the waveform at each station could be compared to a standard waveform at the station so that only the results of the comparison need to be o transmitted to a central station (any one of the four sations).

In one possible embodiment for the present invention, six stations 20, 21, 22,23 24, are located as shown in figure 2. The likely location errors are indicated in kilometres by error contours labelled on the map, for example 3.2, 5.6, 10. This assumes that the error in the arrival time differences for all pairs are 10 is which is believed to be achievable using S .50kHz sampling.

A preferred method of locating a lightning strike using WBI and a preferred receiver network is shown in figure 3. The preferred method will also be described with reference to results obtained during an example used to test the system, in which the sampling rate was 30 kHz giving an upper frequency limit of 15 kHz. At each receiver, the wide band of radio waves is continuously sampled 30 and a representative waveform stored in a buffer. When a sferic is detected by the stations 31 the waveform in the buffer for each station is stored 32. In the test the sample length was about 35 ms (1024 samples) of which ms occurred before the trigger point. However only the first 512 samples (17 mins) were used in this test since the duration of the sferic wave train is less than 1 ms. Figure 4 shows the sampled sferics from a single lightning strike received at Tumbi Umbi in Australia and at Dunedin, New Zealand (separation distance of 2138 km) on 20 June, 1999, at 04:59:20 UT (about 5 pm local time in Dunedin, 3 pm in Tumbi Umbi).

From figure 4 it appears that the two sferics which are sampled (digitised) 30 at the sites of reception arrived at those stations almost simultaneously. This is an illusion since the -9time scale zero points are relative to the times of detection 31. These detection or trigger times are known with precision from local clocks continuously synchronised by GPS.

When a sferic is detected at a station, a time stamp is provided on the sferic to record the time of detection. For each station pair a parameter t o is determined 33 which represents the difference in the time of detection of the sferics received at each station of the pair. t o is calculated using the time stamps corresponding to the time of detection of the sferic at each station of the pair.

The time of detection at each station is typically known to within a few hundred nanoseconds. As the detection of the sferic at each station may be triggered from a different portion of the sferic waveform, this time is not the (relative) time of arrival (TOA) of the sferic, but is usually within a hundred microseconds of the TOA, and is dependent *upon which part of the sferic waveform triggered at each site.

To determine 35 the ATD for each station pair a residual time is calculated 34 to allow for the fact that the stations are triggered by a different part of the sferic waveform. The residual time (6t) is related to the ATD and t o by the equation t A TD-to. The difference residual time can be calculated directly using the rate of change with respect to frequency of the phase difference between respective frequency components of the waveforms for the emissions detected at each station of the pair.

In an alternative embodiment the residual time difference can be obtained indirectly using individual residual times, t, relating to the emission received at each station. The difference .residual time for a station pair is the difference between these times, namely &t=t/-t 2 where t, is the individual residual time in respect of the emission received at the first station and t 2 is the individual residual time in respect of the emission received at the second station.

The individual residual time for each station is calculated by performing WBI at VLF on the emission and a reference or standard waveform.

A single lightning stroke produces a sferic at each of a pair of spaced receivers. The recorded sferics will appear, as in figure 4, within about 3 sample periods (60 is for our intended sampling rate of 50kHz) of one another. This is because the sampling start is triggered locally by the reception of the sferic at each receiver. This has the virtue that the two recorded sferic waveforms overlap almost completely so that the Fourier transforms of each waveform and F 2 apply to the same relative time span. The offset, to, of the sampling spans at the two sites is known as described above.

To determine 6t directly we assume that the two sferics are identical (zero noise and propagation distortion) except for a residual displacement in time of st. This means that F (o)=F,(w).exp(qw6t) Where 6t is positive if sferic waveform #1 is a delayed version of Thus, ln(F (w))=ln(F 2 +j t=ln(F 2 +j q(w) 2 Alternatively, the individual residual times, t, can be calculated by using a standard waveform. Such a standard waveform must be precisely defined in time and cover all frequencies. The ideal waveform is a Dirac delta function which can be described as a pulse of infinite amplitude but zero width such that the product (amplitude times width) is unity. The Delta function is defined by o(to,) where to, is the detection or trigger time at station i. Its Fourier transform is exp(Ywtoa). At station #1 the Fourier transform of the sferic waveform is (wo)exp(jwt). Since F(O=exp(wt 01 then F (o)|exp(ja(t-toi), where t-to, is the individual residual time, for station #1 which is calculated on site at station Meanwhile, the same sferic is received at 25 station In the same way, the individual residual time, t 2 for station #2 is calculated on site at station Assuming, as done above in the direct comparison, that the two sferics are identical except for a residual displacement in time, then the magnitudes of their Fourier transforms are identical. Thus ln{F,(w)/F 2 (w)}=jw(t-t 2 )=-jAq o the same as obtained directly.

It should be noted that both the direct and indirect determination of 6t assumes that the sferic waveforms received at the stations are the same except for a relative displacement 11in time. This is not strictly true, since the sferics travel over different terrain and for different distances in general, so there is a small error which research so far indicates is not serious. Also, since the waveform of a sferic and that of a delta function are very different, t, and t 2 are not the respective times of arrival (TOAs) because there is an offset from the normal definition of the TOA (the time of arrival of the earliest component of the sferic) which depends on the sferic structure. However, if the sferics arriving at stations #1 and #2 have the same structure, this offset is the same for both so that 6t=t-t exactly. By direct measurement we find the direct and indirect determinations of 6t from any pairs such as those in figure 4 give exactly the same value of 6t regardless of errors since any errors or offsets through using an intermediary cancel as shown mathematically above.

The residual time difference of arrival (6t) is obtained from the rate at which the phase difference, Ad(w), changes with frequency, c: 6t =dA(o)dw the units in radians/radian/second or in cycles/Hz, both give st in seconds. Thus 6t =dAp(f)/df where 6t is the residual time difference as a function of frequency, Ap(f) is the phase i difference in cycles between respective frequency components of the waveforms detected at each station of the pair andf is the frequency in cycles/second (hertz).

25 Therefore the equations can be a function of frequency (f0 rather than angular frequency (w) to produce the following equivalent equation which can be used Ap(t)=-jln{F, 2 where Ap(f) is the phase difference in cycles between respective frequency components of the waveforms detected at each station of the pair,f is frequency in cycles/second, is the Fourier transform of the waveform representing the VLF portion of the emission detected at a first receiver station of the pair, and F 2 is the Fourier transform of the waveform representing the VLF portion of the emission detected at a second receiver station of the pair.

-12- If the sample containing both sferics is of duration T, then the frequency interval, df, in the FFT (the separation between adjacent channels in a spectrum analyser) is If the change of difference phase AO within df is a complete cycle the output is zero. This would correspond to t dAqidco T, which can never arise because the sferics are not then in the same sampling duration. Even with our present sampling duration of 35 ms as in Figure 3, the two sferics from a common lightning automatically typically overlap with a mismatch of less than 100 us.

Figure 5a shows the plot of difference phase versus frequency for the two received waveforms shown in figure 4. In practice the residual time difference cannot be calculated directly from 6t =dAp/df due to propagation delay and the presence of noise in the received signals. Therefore a representative residual time difference is determined.

Preferably to do this the full frequency range of the waveform is divided into segments and a straight line approximation of the plotted phase difference versus frequency is determined for each segment. For example to find the representative residual time difference from the values plotted in figure 5a, the graph is divided into 14 segments of lkHz, and the graph in each of the segments is fitted to least squares straight lines. The representative residual time difference for the pair is taken as the mean, or alternatively the median, of the slopes of the straight line approximation in each segment after omitting those values affected by interfering signals and resonances near mode cut-off frequencies. Therefore, using this method for finding the ATD, all or substantially all of the waveform which represents the sampled VLF portion of the lightning emission is used. It will be appreciated that while this is a preferred method for finding the representative residual time difference other methods may be envisaged by those skilled in the art.

S The average slope of each 1kHz segment of the graph of figure 5a is plotted in Figure 6 as a circle. The median (mean of middle two slopes) is shown as a horizontal line. The median value is 344ts which is the representative residual time difference for this receiver pair. The error printed (±18 4s) is half the difference between the upper and lower quartiles. Another method is to omit values which are twice this difference (±36 [is in this example) from the median and then calculate the mean of the remaining values. This would omit the two extreme values, give a mean 6t of 344s and probable error of 15 ts.

In the alternative embodiment in which an individual residual time is first calculated, an emission received at a station is compared to a reference waveform such as a Dirac delta -13function. The individual residual time is found from the rate of change with respect to frequency of the phase difference between the emission received at the station of interest and the Dirac delta function using the technique as described above. Figure 5b shows a plot of the phase difference with respect to frequency between the waveform received at each station shown in figure 4 and a Dirac delta function. The difference between the calculated individual residual times can be found to specify the residual time difference for the station pair.

If the two sferic wave forms are identical except for a displacement in time, in a further alternative method to find the ATD it would be possible to add increments of time to the first arrival until ddqYdcw O, using WBI to determine this condition. This means that all waves (or frequencies) over the whole band are exactly in phase at and only at 6t O. Thus all waves cross zero at 6t 0 making that point very sharply defined as seen in figure 7a, which uses only 11 waves (instead of 1000) for illustration purposes.

1 As seen shown in figure 7b, adding white noise of 10 dB less power than the sferic components (noise having about a third of the RMS amplitude of the sferic) appears to spoil the sharpness. However, as shown in Figure 8, by using all 1000 frequency components, *se this same level of noise (10 dB below sferic power) produces an error in t of only 0.2 ps.

This is despite random phase errors of up to 900 produced by the noise. Close inspection shows a change in the average phase but little discernible change in the phase slope.

As seen in figure 4, the sferic occupies less than 2% of the 35 ms time span. Since most .sferics last less than 0.6ms, the sample duration can be reduced to 1 ms beginning 0.2ms before the trigger time leaving 0.8ms for the rest of the sferic. This has the advantage of greatly reducing the waveform data to be transmitted to a central station or all the other stations. It also removes later sferic (like that at 15ms in Figure 4) from the sample triggered by the first sferic and allows them to trigger their own samples. Thus sferics occurring only Ims apart could all trigger samples and so their causative lightning strokes could be located. Such occurrence is rare since even world wide the average lightning rate is less than 100 strokes per second. To provide adequate frequency resolution in the FFT, the effective sample length can be increased from 50 data points (Ims duration for a sampling rate of 50 kHz) to 512 for efficient FFT calculation simply by adding zero fill.

-14- Another way of further reducing the data to be transmitted to all the other stations is to represent the whole waveform by the local (at station i) individual residual time, together with its likely error obtained from the phase regression showing errors beyond a set limit (about 10% of sferics) could be simply omitted from location measurements. The full waveforms can be archived on site for later analysis and research.

Once the ATD has been determined 35 for each receiver station pair of the network, the lightning location in latitude and longitude can be deduced 36 from these ATDs. For example, if there are 4 receiver sites, then these can form a combination of 3 independent pairs of sites, and from these 3 site pairs, 3 ATD values can be determined 35. The relative time of arrival (TOA) at each station can be calculated from the ATD values. The TOAs are relative to the TOA at the station which detected the sferic first. This can be taken as the trigger time at that station since the absolute TOA is not needed for location. This is equally true for lightning location methods which are called TOA methods.

1 Methods for determining lightning location from ATDs 36 or TOAs are well know in the art, and such methods are described in for example Lee, A. C. Ground truth confirmation and theoretical limits of an experimental VLF arrival time difference lightningflash locating system, Q. J. R. Meteorol. Soc., 115, 1147-1166, (1989). A brief description of a possible method is given here, however a detailed description will not be given as it will be appreciated that those skilled in this area of technology will be aware of the methods that can be used.

In a preferred method, the lightning location is determined 36 by first getting a rough estimate of the lightning location, for each station pair calculating a theoretical ATD value which would result if the estimate was actually the lightning location, comparing these theoretical ATDs with the actual measured ATD values to find a better estimate of the lightning location, and recalculating the theoretical ATD values from the better estimate.

This method is re-iterated until the difference between the calculated theoretical ATD values, for the lightning position and actual measured ATD values are minimised, giving more weight to those ATD values likely to be most accurate. Although the first iteration can assume a flat Earth, later ones must take into account the spheroidal Earth and nature of each path from lightning to each receiver (whether in day or night, over ground or sea, etc).

I f Lee [1989] does not state the actual method used to perform the minimisation in the iterative process. In the preferred embodiment of the invention the minimisation method titled "Powell's direction set methods in multi dimensions" disclosed in Press, W. S.

A. Teukolsky, W. T Vetterling and B. P. Flannery, Numerical Recipes in C. The Art of Scientific Computing, Cambridge University Press, Second Edition, (1992) is used.

For each iteration, ATD values for the assumed lightning position are determined by calculating the great circle path length differences on the spheroidal earth using Rudoe's reverse formula as described by Bomford, Geodesy, Clarendon Press, Oxford, (1980).

A constant phase velocity of 1.004c which is appropriate but curiously uses the phase velocity instead of the group velocity to determine the ATD. He justifies this on the empirical grounds that it reduces the value of the residual by a factor of between 2 and 4 compared to using 1.000c. However, the apparent path length depends on the assumed •altitude of the path (since the Earth is round!) in the Earth-ionosphere waveguide through which the long wave components of the lightning travel. At 50 km altitude, which is at about the electrical centre of the Earth-ionosphere waveguide, the path is 1.008 times longer than at sea level. Thus a wave propagating in the waveguide at an effective group velocity of 1.000c would have a ground speed of 0.992c.

While it is clear that the cross-correlation method should use the wave group velocity (because it matches the wave group or packet as received at two point), this is true also for the WBI method. This is because the group velocity is determined by the rate of change of the phase velocity with respect to frequency and the WBI method measures the ATD as S* -the rate of change of the phase difference with respect to frequency.

Figure 9 shows a block diagram of a preferred embodiment of hardware associated with a receiver station to be used for the present invention. It will be appreciated by those skilled in the art that this embodiment is by way of example only, and that other hardware setups and/or configurations different to that described could be used. Not all the hardware comprising a receiver station is necessarily located at the same site.

The VLF antenna 91 is a 1.5 m whip for receiving radio waves. Mounted at the base of the whip antenna is a wide band (3-50 kHz) preamplifier (PA) 92. Preferably these are placed on a tall building in a city at an Internet Service Provider's (ISP) location and can be kilometres from the remainder of the equipment, but preferably only a few tens of metres.

16- Passive filters in the preamplifier 92 reduce the effects of power line harmonics and very strong transmissions at high frequencies such AM and FM broadcasts, TV, radar, and the like on the performance of the receiver station.

Apart from the GPS 98 antenna which preferably is within line of sight of most of the sky, the rest of the equipment is preferably located together, though it need not be easily accessible. The Analogue-to-Digital Converter (ADC) 93 samples at about 50 kHz. Its frequency is not critical since it is measured by the number of counts (sample points) between successive GPS pulses (1 pulse per second) which are emitted from the GPS receiver 96 within a few hundred nanoseconds of UTC (Universal Time). The trigger time of the sferic is measured with respect to these GPS pulses by WBI to 0.Sus. The necessity for the 10 MHz crystal 94 depends on the type of ADC used. The GPS receiver also provides the time code, which can be TSIP as used in the preferred embodiment or NMEA, and provides the precision time stamp of the sferic sample °The hardware further includes a PC 97. Preferably the PC 97 is permanently on line to the Internet 14 and which enables remote access for any software modifications and operating procedures, including monitoring usage by paying clients. The PC 97 uses the Linux operating system, however it will be appreciated that any other suitable operating system could be used. The DSP 95 in conjunction with the PC 97 perform signal processing in accordance with the method described above to determine the lightning location. In a .preferred embodiment all receiving stations communicate via the Internet 14 and each of them process all the available sferic information using their respective DSPs and Pcs. This S. processing includes sampling and storing in a buffer radio waves received via the antenna 25 91, preamplifier 92 and analogue to digital converter 93, detecting when a sferic has been received, calculating the ti and t o and on receiving the counter part data from the other stations, calculating 6t for each station pair, and the lightning location from the ATDs.

Having each station process all the available sferic information avoids a "master-slave" relationship so that if any one station goes out of operation the whole network is merely degraded slightly. Such an occurrence will result in decreased lightning location accuracy in some areas but little change in others. It will be appreciated however that while this is the configuration of a preferred embodiment of a receiver network, other network configurations known in the art could be used.

-17- To minimise stations going out of operation due, for example, electric power loss, the system is preferably powered by an Uninterruptible Power Supply (UPS). As further protection, the PC 97 can be rebooted remotely by an Internet-controlled hardware switch.

The lightning location which is determined 36 by the system can be displayed to subscribers via a website 37. Preferably the website will show 37 a graphical representation of lightning strike locations. The entire process of determining the time difference, t o t, between every pair of sites 35, calculating the position of the lightning strike from these times 36 and presenting the data 37 via the Internet in a form useful to the user, is preferably done within seconds. Thus the user can see the lightning strike position on a computer screen map almost as it happens. The user also has access to historical data relating to strike positions of earlier lightning in the same storm preferably differentiated by colour on the screen, as well tabulated data covering previous months and years and/or requested areas.

o oo 6*6o 15 In such a system the lightning position (latitude and longitude) applies more to the source 6..

(in the thunder cloud) than the ground striking position. This is "close to the updraft region of a convective cell which is the source of all hazardous thunderstorm produced meteorological conditions including: severe turbulence, hail, icing (supercooled water), lightning, microbursts, windshear, gust fronts, heavy rain and tornadoes" [see Markson, R.J, and Ruhnke, L. New technologyfor mapping total lightning with groundbased and airborne sonsors, Abstract El. 2, URSIAbstracts, 265 1999]. Locations made available in •:4.6i real time (a few seconds after the lightning) would be of primary use for "now casting" advanced (by minutes or hours) warning of hazards to civil defence bodies, aircraft, boats, foresters, farmers, power transmission line operators, etc. The lower location accuracy of individual lightning strikes (a few kilometres instead of a few hundred metres) is of no consequence for this since such "now casting" is based on an advancing front of lightning strike positions.

Claims (15)

1. A method for determining the location of a lightning discharge including: receiving a discharge emission at a plurality of receiver stations, determining a time of detection of the emission at each station, for at least three pairs of stations, comparing the emission received at one station of the pair with the emission received at another station of the pair to determine, using wide band interferometry (WBI) at VLF, the arrival time difference between the emission received at each station of the pair based on the comparison and the times of detection, and determining the location of the discharge using the arrival time differences determined for each station pair.

2. A method according to claim 1 wherein using WBI to determine the arrival time difference includes: determining a residual time difference caused by the detection of the emission at each station of the pair being triggered at a different amplitude level of the emission, wherein the residual time difference for each station pair is the rate of change with respect to frequency of the phase difference between frequency components of the emission received at each station of the pair.

3. A method according to claim 2 wherein the residual time difference is determined approximately using the average rate of change with respect to frequency of the phase difference between the frequency components.

4. A method according to claim 3 wherein the average rate of change is determined by: dividing the emission received at each station of the pair into a plurality of corresponding frequency bands, calculating the average rate of change with respect to frequency of the phase difference between respective. frequency components contained in each corresponding band, and calculating the median of the average rates of change calculated for each corresponding band.

19- A method according to claim 4 wherein the arrival time difference for each station pair is calculated using the equation: atd to where atd is the arrival time difference,.t o is the difference in time of detection and a, is the residual time difference. 6. A method according to claim 5 wherein a VLF portion of the emission received at each station is used for determining the residual time differences. 7. A method according to claim 6 wherein the rate of change with respect to frequency is specified by the equation: 6 =dAp(f)/df where t is the residual time difference in seconds, Ap(f) is the phase difference in cycles between respective frequency components of the emission detected at each station of the pair, andfis frequency in cycles/second. 8. A method according to claim 7 wherein Ap(f) is specified by the equation: Ap(f)=-jln )/F 2 S. 25 where is the Fourier transform of the emission detected at a first receiver station of a pair, and F(f) is the Fourier transform of the emission detected at a second receiver station of a pair. 9. A method for determining the location of a lightning discharge including: receiving a discharge emission at a plurality of receiver stations, determining a time of detection of the emission at each station, for each station, comparing the emission with a reference waveform to determine an individual residual time using WBI at VLF, for at least three pairs of stations, determining a residual time difference using the individual residual times for each station of the pair, and for at least three pairs of stations determining an arrival time difference for each station pair using the residual time and times of detection, and determining the lightning location using the arrival time difference for each pair. 10. A method according to claim 9 wherein the individual residual time for a station is determined from the average rate of change with respect to frequency of the phase difference between frequency components of the received emission and the reference waveform. 11. A method according to claim 10 wherein the average rate of change is determined by: dividing the received emission and the reference waveform into a plurality of corresponding frequency bands, calculating the average rate of change with respect to frequency of the phase difference between respective frequency components contained in each corresponding SSo. band, and .calculating the median of the average rates of change calculated for each corresponding bands. 12. A method according to claim 11 wherein the arrival time difference for each station pair is calculated using the equation: atd to r where atd is the arrival time difference, t o is the difference in time of detection and Ar is the residual time difference. 13. A method according to claim 12 wherein the reference waveform is a Dirac delta function. 14. A system for determining the location of a lightning discharge including: a plurality of stations for receiving an emission from the discharge, the stations configured as pairs, one or more time stamping means for recording the time of detection of the emission at each station, and -21- a signal processor associated with each station pair for implementing WBI at VLF on the emission received at each station of a station pair to determine a residual time difference for that station pair, wherein the processor is further adapted to calculate an arrival time difference of the emission at each station of an associated pair using the times of detection and the residual time difference, and is further adapted to determine the lightning location using the arrival time differences calculated for each pair. A system according to claim 14 wherein the lightning location is determined using the arrival time differences calculated for each pair 16. A system according to claim 15 wherein the processor calculates an approximate residual time difference for each station pair by using the average rate of change with respect to frequency of the phase difference between frequency components of the emission received at each station of the pair. *0 0 S"17. A system according to claim 16 wherein the processor determines the average rate of change by: S•dividing the emission received at each station of the pair into a plurality of corresponding frequency bands, calculating the average rate of change with respect to frequency of the phase difference between respective frequency components contained in each corresponding band, and calculating the median of the average rates of change calculated for each 00 25 corresponding band. 18. A system according to claim 17 wherein the arrival time difference for each station o• pair is calculated using the equation: 0000 0000 atd to r where atd is the arrival time difference, t o is the difference in time of detection and &r is the residual time difference. -22- 19. A system according to claim 18 wherein the receiver stations are communicably linked via the Internet. A system according to claim 15 wherein the processor determines a residual time difference for a station pair by calculating individual residual times relating to the emission received at each station of the pair.

21. A system according to claim 20 wherein the processor calculates an residual time difference for a station by using the equation: dtr=tI- t2 where 4r is the residual time difference, t, is the individual residual time for the first station of the pair and t 2 is the individual residual time for the second station of the pair.

22. A system according to claim 21 wherein the processor calculates an individual residual time for a station by using the average rate of change with respect to frequency of the phase difference between frequency components of the emission received at the station and a reference waveform.

23. A receiver station for use in a system which determines the location of a lightning discharge from an emission, the station associated with another receiver station to form a station pair, the station including: time stamping means for recording the time of detection of the emission o 25 received at the station, and a signal processor for determining a residual time difference for the station pair from the received emission by implementing WBI at VLF, wherein the processor is adapted to calculate an arrival time difference of the emission at each station of the pair using the time of detection and the residual time difference.

24. A station according to claim 23 wherein the processor calculates an approximate residual time difference by using the average rate of change with respect to frequency of the phase difference between frequency components of the emission received at each station of the pair. -23- A station according to claim 24 wherein the processor determines the average rate of change by: dividing the emission received at each station of the pair into a plurality of corresponding frequency bands, calculating the average rate of change with respect to frequency of the phase difference between respective frequency components contained in each corresponding band, and calculating the median of the average rates of change calculated for each corresponding band.

26. A station according to claim 25 wherein the arrival time difference is calculated using the equation: atd t o +t S where atd is the arrival time difference, t o is the difference in time of detection and Ar is S•the residual time difference. S• .27. A station according to claim 23 wherein the processor determines a residual time difference for the station pair by calculating individual residual times relating to the emission received at each station of the pair.

28. A system according to claim 27 wherein the processor calculates a residual time 0:0 5 difference for a station by using the equation: dt=tl- t2 where dtr is the residual time difference, t, is the individual time difference for the first °station of the pair and t 2 is the individual time difference for the second station of the pair.

29. A system according to claim 28 wherein the processor calculates an individual residual time for a station by using the average rate of change with respect to frequency of the phase difference between frequency components of the emission received at the station and a reference waveform. -24- A receiver station according to claim 23 further including means for enabling connection to the Internet.

31. A method for determining the location of a lightning discharge substantially as hereinbefore described with reference to the accompanying drawings.

32. A system for determining the location of a lightning discharge substantially as hereinbefore described with reference to the accompanying drawings.

33. A receiver station for use in a system which determines the location of a lightning discharge from an emission substantially as hereinbefore described with reference to the accompanying drawings. *i 0^