A source location algorithm of lightning detection networks in China

Fast and accurate retrieval of lightning sources is crucial to the early warning and quick repairs of lightning disaster. An algorithm for computing the location and onset time of cloud-to-ground lightning using the time-of-arrival (TOA) and azimuth-of-arrival (AOA) data is introduced in this paper. The algorithm can iteratively calculate the leastsquares solution of a lightning source on an oblate spheroidal Earth. It contains a set of unique formulas to compute the geodesic distance and azimuth and an explicit method to compute the initial position using TOA data of only three sensors. Since the method accounts for the effects of the oblateness of the Earth, it would provide a more accurate solution than algorithms based on planar or spherical surface models. Numerical simulations are presented to test this algorithm and evaluate the performance of a lightning detection network in the Hubei province of China. Since 1990s, the proposed algorithm has been used in many regional lightning detection networks installed by the electric power system in China. It is expected that the proposed algorithm be used in more lightning detection networks and other location systems.


Introduction
Cloud-to-Ground (CG) lightning is one of the most dangerous atmospheric phenomena, often causing injuries and deaths, power system interruptions and equipment damage.Ground-based lightning detection networks have been installed all over the world to provide useful information for CG lightning monitoring and research.Some providers are Correspondence to: Z. X. Hu (zhixianghust@gmail.com) the US National Lightning Detection Network (NLDN), the lightning detection network in Europe (LINET), Zeus in Europe and Africa, and the World Wide Lightning Location network (WWLLN) (Cummins et al., 1998a;Betz et al., 2009;Chronis and Anagnostou, 2006;Dowden et al., 2008;Rodger et al., 2005Rodger et al., , 2006)), etc.These long-range lightning detection networks primarily locate cloud-to-ground (CG) lightning which is hazardous to human society.The data provided by these networks are widely used in electric power systems and for the aviation industry and meteorology research.
Modern lightning detection systems are equipped with multiple time-sync remote sensors capable of recording the time-of-arrival (TOA) and/or azimuth-of-arrival (AOA) of the electromagnetic wave radiated by CG lightning (Krider et al., 1976;Lee, 1986).Thus the unknown lightning location and onset time can be determined by these measurements.A variety of algorithms have been proposed for the lightning location (Cummins and Murphy, 2009;Koshak et al., 2000;Koshak and Blakeslee, 2001;Cummins et al., 1998b;Dowden et al., 2002); these can be roughly divided into two categories: one works by solving the non-linear equations (Koshak et al., 2000;Koshak and Blakeslee, 2001) and the other by searching for the optimum position in solution space (Cummins et al., 1998b;Dowden et al., 2002).Under the assumption that the electromagnetic wave propagates along a flat or spherical Earth, a pedagogical effort has been made by Koshak et al. to linearize the observation equations and give a non-iterative solution.Since the ground wave propagates along the surface of an oblate spheroidal Earth, the above planar or spherical surface assumptions are not appropriate.To address this problem, an iterative oblate method taking the Earth's oblateness into account was introduced by Koshak and Blakeslee (2001), but at least four sensors are needed for providing an initial position and the iterations often degrade the initial value as well.In addition, the iterative oblate method can not incorporate the azimuth of arrival data for improving the location accuracy.In 1992, the improved accuracy from combined technology (IMPACT) was developed for lightning locations detected by NLDN (Cummins et al., 1998b).The IMPACT algorithm determines the best estimate of the lightning location by iteratively moving the stroke position along the surface of an oblate spheroidal Earth.Although the IMPACT algorithm is commercially successful, the specific software is proprietary and not widely distributed free of charge to the scientific community, and its data processing ability is constrained by computer power.In WWLLN, the lightning strokes are located using the "downhill simplex" method (Dowden et al., 2008;John and Mead, 1965), which is also a searching method similar to the IM-PACT algorithm.
In this paper, an efficient CG lightning location algorithm is presented utilizing both TOA and AOA data.This algorithm has the advantages of retrieval lightning location and of onset time directly on the oblate spheroidal Earth surface.The earlier version of this algorithm was developed by Zhao et al. (1999), and since then it has been widely used in many regional lightning detection networks installed by the electric power system in China.Because of the demand for lightning monitoring for the space launch center in Wenchang, Hainan, we plan to install several lightning sensors in this region before 2012.For this purpose, we have improved Zhao's algorithm and are evaluating it in the present study.Compared to the earlier version of this algorithm, the explicit initial position calculating method is more accurate after applying a simple unit circle model for choosing the best combination of three sensors.Furthermore, in the present version of the algorithm, a high accuracy expression of the geodesic distance calculating method is used, which enables this algorithm to be employed for broader lightning detection networks.And for the first time, the algorithm is tested using a series of numerical simulations.In the ensuing section, the mathematical theory of lightning location on an oblate spheroidal Earth surface is introduced, as well as a unique method to compute the geodesic distance and azimuth.In Sect.3, an explicit method for calculating an initial position of lightning source position using time of arrival measurements measured by three sensors is presented, and then the least-square position is computed iteratively.The performance of the proposed algorithm is tested in Sect. 4.Then, conclusions are drawn in Sect. 5.

CG lightning location principle
As shown in Fig. 1, the 1984 World Geodetic System Earth Ellipsoid (WGS-84) is used in our work: B = geodetic latitude and L = longitude.The curvature radius at the poles C 0 = 6399593.626m, the first eccentricity e 2 = 0.00669437992 and the second eccentricity e 2 = 0.006739496745.The lightning source and the sensors are assumed to be located on the Earth surface; altitude is not considered in this paper.
Consider that the position of a CG lightning is (B,L) and its onset time is t.The TOA recorded by the i-th sensor is denoted as t i , i = 1, 2, •••, n, and the AOA as α i , the arrival time and azimuth equations can be written as where c is the speed of light in air, S iP and β iP are the geodesic distance and azimuth between sensor i and the lightning source, ε T i and ε Ai are the measurement errors of arrival time and azimuth, respectively.The measurement errors are due to the effects of the nondeterministic nature of the lightning signals (Thomas et al., 2004), the propagation over complex terrain (Cummins et al., 1998a;Schulz and Diendorfer, 2000) and the clock error, and the reflexion of the signal will cause systematic error of AOA measurement (Orville, 1991).We compile all these effects into Gaussian error models, i.e., the errors are normally distributed with zero mean values and predefined variances which is denoted as σ 2 T and σ 2 A in this paper.An efficient method for computing the geodesic distance and azimuth is given below, with which one can calculate the geodesic distance and azimuth with high accuracy.These formulas were first derived by Zhao (1997) for ellipsoidal geodesy research.As shown in Fig. 2, let (B i ,L i ) denote position of the i-th sensor, then the geodesic distance S iP and azimuth β iP can be calculated by Table 1.Adjustment of the calculated value of the geodesic azimuth.
where l iP = L − L i is the longitude difference between the two points, ϕ i is a correctional angle, θ iP is radian measure distance and is the radius of curvature in prime vertical at the i-th sensor position.Equations ( 3) to ( 5) are closed formulas derived from the properties of the spherical surface; Eq. ( 6) is truncated forms from an infinite series that approaches the geodesic distance.Since the value of geodesic azimuth calculated by Eq. ( 4) lies in the range [− π 2 , π 2 ] rather than the expected range [0, π], the calculated value should be adjusted, as listed in Table 1.If the azimuth is very close to π 2 or 3π 2 which corresponds to an unstable (extremely large) tangent value, the azimuth should be carefully determined to eliminate calculation error.In these two cases, the azimuth is close to π 2 when L i < L P or close to 3π 2 when L i > L P .An exchange of the subscript of Eq. ( 3)∼(6) can give (B P + ϕ P ), β P i , θ P i and SP i .The values of SiP are found to coincide favorably with the centimeter accuracy of SP i .The geodesic distance is finally given by When S iP is shorter than 500 km, the calculation errors of geodesic distance and azimuth are below 0.05 m and 0.01 arc sec, which is accurate enough for lightning detection networks.For longer geodesic lines, we will verify this method by the results of a set of higher accuracy forms and compare it with Sodano's forms (Sodano, 1965) in Sect.4; the latter has been used by Koshak et al. for the lightning source location.Koshak et al. has proven that one can use Sodano's forms to calculate geodetic distance within centimeter accuracy.The detailed derived process of the proposed forms Eq. ( 3)∼( 6) is beyond the scope of this paper, for further details, the reader is referred to Zhao (1997).
Given multiple (n ≥ 3) sensors' recorded data of signal arrival time and azimuth, the CG lightning source can be retrieved from the equation set composed of Eqs. ( 1) and ( 2).Since the equations are nonlinear, a Newton iteration method can be applied.Good initial position and onset time are required for the algorithms' convergence.The initial position computation method will be discussed in the subsequent section, as well as the detailed procedure of the iteration method.

Explicit initial values
Time-of-arrival geolocation techniques were primarily developed for marine navigation systems like the Loran system.Razin had derived an explicit position calculation method for the Loran system (Razin, 1967).The position calculation method that we present here is similar but not identical to Razin's explicit solution.Our method is simpler than Razin's method and no osculating sphere (best-fit sphere) approach is needed for the approximation of the Earth's surface surface.This method can provide initial CG lightning position and onset time for the subsequent Newton iteration algorithm which will give an improved solution.
Consider that recorded data of sensor 1, 2 and 3 are used to calculate the lightning location and onset time, as shown in Fig. 3. S 12 , β 12 , S 13 , and S 13 can be calculated using Eq. ( 3)∼( 8).Let θ ij denote the radian measure distance from point i to j , which can be estimated by here R is the mean radius of curvature of the network composed of sensor 1, 2 and 3, According to the properties of arc lengths on a sphere, the spherical law of cosines can be written as where φ 1 = β 1P − β 12 and φ 2 = β 1P − β 13 .Let ij denote the radian measure range difference of a pair of sensors to the location of the lightning source, defined by ij = Now Substituting Eq. ( 13) into Eq.( 11) yields two equations in θ 1P and β 1P which, after transposition and division by sinθ 1P , become (cos 12 − cosθ 12 )cotθ 1P = sin 12 + sinθ 12 cosφ 1 (cos 13 − cosθ 13 )cotθ 1P = sin 13 + sinθ 13 cosφ 2 .(14) Note that φ 1 and φ 2 can be expressed by β 1P .A forward elimination of cotθ 1P gives a standard trigonometric equation below: From Eq. ( 15), an initial value of the azimuth of the geodesic from lightning source position to sensor 1 can be obtained, given by where Note that β 0 1P has two candidate values, one given by Eq. ( 16) (which should be plus 2π if it is negative), one given by subtracting the value of Eq. ( 16) from π.When four or more sensors' data are available, there are many combinations for calculating of the initial azimuths of arrival of a chosen sensor #1.Although each combination gives two candidate values of β 0 1P , the difference between all truth-values is very small.Thus the false-values can be eliminated.If only three sensors are trigged by the lightning signal, the falsevalues can be eliminated using azimuth-of-arrival measurement information.In some extreme cases, both of the calculated values of β 0 1P are used in the subsequent computation and both the two final solutions are recorded.
Hence the initial value β 0 1P is obtained; the initial value of geodesic distance from the lightning source position to the selected sensor #1 can be calculated by where θ 0 1P can be calculated from any one of the equations in Eq. ( 14).
Finally, considering the property of spherical triangles, the initial values of lightning source position and onset time can be calculated by

Least square solution
Newton iteration method is frequently used for solving nonlinear problems (Torrieri, 1984), which needs to linearize nonlinear terms in the equations and then solve them iteratively (Chan and Ho, 1994).In this work, the unknown lightning source position and onset time are determined by Eqs. ( 1) and ( 2), which can be linearized by expanding the Ann.Geophys., 28,[1981][1982][1983][1984][1985][1986][1987][1988][1989][1990][1991]2010 www.ann-geophys.net/28/1981/2010/right-hand side of Eqs. ( 1) and ( 2) in a Taylor series and keeping only terms below second order.Denote the errors of the initial values provided by the explicit method above as δB = B − B 0 P , δL = L − L 0 P and δt = t − t 0 P .The approximate set of equations can be written as where i = 1,2, •••,n, and the notations ∂ ω S iP and ∂ ω β iP are the partial derivatives of the geodesic distance and azimuth to ω, respectively, with ω = B,L.All of these derivatives are evaluated at (B 0 P ,L 0 P ).Secant-type approximations of the derivatives may be used (Kohsak and Blakeslee, 2001), but we have derived exact representations for these derivatives, where N 0 P is the radius of curvature in prime vertical at the initial lightning source location, calculated by Eq. ( 7), thus N 0 P cosB 0 P is the radius of curvature the parallel circle; and M 0 P is the radius of curvature in meridian, The equations defined by Eqs. ( 19) and ( 20) can be rewritten as an equation system for m = 2n row measurements: Denote the coefficient matrix as A, the unknown vector as x and column vector on the right-hand side as y, Eq. ( 26) can be expressed as This is an overdetermined linear system and the least square solution can be given by where W is the covariance matrix of the vector denoted by y, Therefore W −1 is regarded as a matrix of weighting coefficients.In practical applications, the matrix W −1 can be replaced by assuming that the variances of time-of-arrival measurements are σ T = 1 µs, and variances of azimuth-of-arrival measurements are σ A = π/180 rad.
Finally, an updated solution for the lightning source location and onset time can be given by This updated solution can also be used as new initial values for a second iteration, but usually the results offered by only one or two iterations are good enough.

Retrieval errors
Let P denote the covariance matrix of x.According to the statistical theory of the passive location (Torrieri, 1984), we obtain here W is the covariance matrix defined by Eq. ( 29).It is also the covariance matrix of vector [B,L,t] T , and can be rewritten as The diagonal elements of P can give the variances of the errors of lightning positions and on set time, and the error ellipses of the estimated location of lightning sources can be inferred by where M P and N P are the radius of curvature in meridian and prime vertical at the estimated lightning source location, respectively.Lengths of semi-major axis and semi-minor axis are given by the eigenvalues of the matrix P .The root mean square error (RMSE) of the source location estimator is If only arrival time measurements are used for source location, the matrix P and RMSE can be derived from Eqs. ( 21), ( 22), ( 29) and ( 32).Given the fixed network configuration   6) and ( 36) (in meter).The square is the beginning of each geodesic line.The difference is below in a large (yellow) area.
and uncertainty of time measurement, it's easy to find that the source location accuracy (retrieval error) is dependent on the azimuths of the geodesic line from the lightning source location to each sensor, i.e., β P i .We have investigated how different combinations of the azimuths β P i will influence the source location accuracy.If all the azimuths are concentrated, say max(β P i ) − min(β P i ) < 30 • , this corresponds to the situation that the lightning source is outside and far away from the network, thus the retrieval error is very large.If the azimuths are dispersive, this corresponds to the situation that the source is inside or near the network, thus the retrieval error is small.This fact can be used to explain why the necessity of "surroundedness" for adequate accuracy applies to all lightning location networks that use timing alone for location (Dowden et al., 2008).

Algorithm performance
In this section, simulation results are presented to test and verify the above algorithm.Because the Hainan lightning detection network is under installation, we chose the Hubei lightning detection network as our test object.The network was installed in 1998 and is operated by the electric power system.This network contains 9 lightning detection sensors.The longitudes of these sensors range from 110 • E to 115 • E, and latitudes range from 29 • N to 33 • N, which will be shown in later figures.First, the geodesic distance and azimuth calculation methods are verified, then the accuracy of the proposed explicit method for determining the initial position is investigated, and finally we present the estimated accuracy of a regional network (Hubei) under the assumption that all sensors participate in the calculation of lightning sources' positions.The computer generated time-of-arrival and azimuth-of-arrival measurements are given by a high accuracy expression (given later), and the retrieval algorithm uses the simplified geodesic distance and azimuth calculation forms Eq. ( 3)∼( 6).We have examined the algorithm under two retrieval scenarios: location with random measurement errors and error free.The former is used to evaluate the spatial distribution of the location accuracy of the network and the latter is used to show the best performance of the proposed algorithm.

Verification of the geodesic distance and azimuth calculations
In order to verify the accuracy of the geodesic distance and azimuth calculation, results of a high accuracy version of forms derived by Zhao (1997) are used as references.This high accuracy version of forms is obtained by replacing Eq. ( 6) with where η 2 i denotes e 2 cos 2 B i , and ψ iP denotes cosβ iP for brevity.Compared to Eq. ( 6), this expression retains many high order terms, thus the geodesic distance can be estimated more accurately.It is worth noting that the geodesic azimuths calculated by both methods are identical because only Eq. ( 6) is changed.Their differences, geodesic distances calculated by Eq. ( 36) minus those calculated by Eq. ( 6), are shown in Fig. 4. The data are obtained by assuming that the lightning source locations are defined on a grid with (113 • E, 31 • N) as its center and a grid resolution of 0.2 • and assuming a lightning detection sensor is located at the center of the grid.It can be inferred that the differences are less than 0.1 m in a large area.Given that the uncertainty of the time measurement of a lightning sensor is about 1 µs, corresponding to a ranging error cε T ≈ 300 m, the geodesic distance calculation is accurate enough for the total investigated area.But we recommend the high accuracy expression ( 36), yet noniterative, to be used in a broader lightning detection network such as WWLLN.
To the authors' knowledge, the noniterative method for calculating geodesic distance and azimuth has not been used for lightning location except Sodano's forms.Koshak et al. has verified Sodano's forms using the numerical method  36) and forms derived by Sodano (in meter).There difference is below 1 cm in a large (blue) area.(Koshak and Blakeslee, 2001, Appendix) and asserted that one can calculate the geodesic distance using Sodano's forms with centimeter accuracy.We have also compared Sodano's solution to our high accuracy expression (36), and the results are shown in Figs. 5 and 6.Note that the discrepancies between the two methods are trivial for most of the investigated area.The difference of the geodesic distance calculated by the two methods is below several centimeters all over the investigated area.And the difference of the geodesic azimuth is below several arc sec, as shown in Fig. 6.Therefore, it is plausible that the geodesic distance and azimuth calculating method are accurate enough for lightning location application.For more discussion of the geodesic distance and azimuth calculation, the reader is referred to Rapp (1993).

Error of the initial value
It is important to know how accurately the initial position and time can be calculated by the explicit method proposed in the above section.An initial position close to the true source's position is essential for the next least square solution.Because of the oblateness of the Earth, the spherical law of cosines (Eq.11) represents an approximate relationship of the variables; thus the explicit initial position must contain some error.The other factors causing the initial position to depart from the true source's position are the random timing errors.Figure 7 is the spatial distribution of the mean distances from the calculated initial position to the true source's location.The black square is the sites of three selected sen-Fig.6. Same as in Fig. 4 except for the geodesic azimuth difference (in arc sec).At the four corner of the map, the azimuth difference becomes large.Consider that the investigated map covers a large area, the azimuth calculating is accuracy enough for lightning location.
sors and the upper sensor is used as sensor #1.This small network is a part of the Hubei lightning detection network.The data are obtained by calculating 500 times the initial position of the fictitious lightning sources on a grid with resolution of 0.2 • , and then the mean distance is computed as the error of the initial position.For each of the 500 trials, an arrival time error selected from a uniform random distribution (ranging from − √ 3 µs to √ 3 µs for Fig. 7b, but error free for Fig. 7a) is added.The elliptical area that gradually becomes large (error increasing from several tens of meters to several kilometers) in Fig. 7a is due to the effect of the oblateness of the Earth surface.Thus the sensor which receives the lightning signal first (has the earliest time of arrival) should be regarded as sensor #1.It can be inferred from Fig. 7b that an initial position of the lightning source can be calculated with error below several kilometers when the source is not along the outer sensor baseline or far away from the network.And the largest error occurs in Fig. 7b when the lightning source is along the outer sensor baseline, which is caused by the timing error and amplified by the fact that the source is not surrounded by the sensors.
If more sensors are triggered by the lightning signal, many combinations of three sensors can be used (e.g., C 3 5 = 20 combinations for five sensors).The best combination can be chosen to avoid the situation that the lightning source is along or near the outer sensor baseline.The principle of selecting the best combination is that the azimuths of the geodesic than the simulated results.The coarse lines in Fig. 10 may be due to numerical error or nonrandomized characteristics of the computer generated data.All of the simulations have shown that the lightning location error is small in the periphery of the network (below 1000 m) and increases when the source is outside the network, especially along the outer sensor baseline.Thus the lightning location accuracy is affected by not only the measurements' accuracy, but also by the relative position between the source and the sensors.

Conclusions
An efficient algorithm for CG lightning source location is presented in this paper.This algorithm has been used in many regional lightning detection networks in China.We have verified the geodesic distance and azimuth calculation method, which is an important part of the proposed algorithm.And the explicit initial position calculation method is evaluated by computer generated data.A simple unit circle model is used to select the best combination of three sensors for calculating the best initial position.In addition, the least square solution is tested using the Hubei lightning detection network.All of the simulations are performed using computer generated data, with or without measurement errors.When measurement errors do not exist, the lightning source retrieval error is only affected by the error of the geodesic distance and azimuth calculation, therefore the retrieval error is favorably small.This is because the proposed algorithm has accounted for the effects of the oblateness of the Earth.When random measurement errors are added to the simulations, a spatial distribution of the lightning source location accuracy is plotted.In the future, the location accuracy of lightning events reported by this network can be evaluated by the presented simulated contour map.Also, the location accuracy may be rapidly evaluated by the covariance matrix presented in Sect.3.
The authors intend to apply the proposed algorithm in other lightning detection networks in the future.Given time, improvement of the accuracy of TOA and AOA measurements will make the proposed algorithm more attractive.Finally, we hope the proposed algorithm, along with the high accuracy expression for geodesic distance calculation, can be used for global lightning detection networks, such as Zeus and WWLLN, et al.

Fig. 1 .
Fig. 1.The oblate spheroidal Earth geometry.B is geodetic latitude, L is longitude of a point P on the oblate spheroidal Earth surface.

Fig. 2 .
Fig. 2. The geodesic path S iP and the azimuth β iP and β P i between the i-th sensor and the lightning source position.P is the lightning location.

Fig. 3 .
Fig. 3.The geometry of CG lightning source retrieval using only three sensors #1, #2 and #3.P is the unknown lightning location; β denotes azimuth; S and θ represent geodesic distance and radian measure distance, respectively.

Fig. 4 .
Fig. 4. The spatial distribution of the geodesic distance calculating difference between Eqs. (6) and (36) (in meter).The square is the beginning of each geodesic line.The difference is below in a large (yellow) area.

Fig. 5 .
Fig. 5. Same as in Fig. 4 except for the geodesic distance difference between the results obtained by Eq. (36) and forms derived by Sodano (in meter).There difference is below 1 cm in a large (blue) area.

Fig. 9 .
Fig. 9.The spatial distribution of the improved accuracy of the least square positions when the measurements are error free (in meter).The lightning location error is below 10 m in the majority area of the investigated region.

Fig. 10 .
Fig. 10.Same as Fig. 7, except for the improved accuracy of the least square positions (in kilometer): (a) lightning source location using TOA measurement only, (b) lightning source location using both TOA and AOA measurements.

Fig. 11 .
Fig. 11.Spatial distribution of the lightning source location error (RMSE) given by Eq. (35), (in kilometer): (a) only TOA measurements are used, (b) both TOA and AOA measurements are used.