An algorithm for estimating the detection efficiency of a lightning location system

Abstract

This article puts forward an algorithm for estimating the detection efficiency (DE) of a lightning location system (LLS). The algorithm can be applied to lightning flash/stroke density correction and LLS performance evaluation. Fundamentally, the generalized extreme value (GEV) distribution was found to best fit the probability distribution of the signal strengths of the lightning observed by the ADvanced TOA and Direction system (ADTD) detectors in Beijing, China. With respect to this probability distribution, we estimated the single-station acceptance damped by the uneven underlying land surface conductivity, which is determined by overlapping the soil-type polygons with the propagating paths of electromagnetic lightning signals. Accounting for the multi-detector location modes supported by single-station acceptance, the iterative algorithm was applied for deducing the network DE. In this case study, the DE estimates were lower in the mountainous areas than in the plains. This can be due to the low underlying conductivity of the mountainous areas, which creates a high attenuation effect on the lightning electromagnetic signals, and the greater distances from the lightning detectors. The CG lightning flash/stroke density corrected for the DEs revealed that the CG lightning flash/stroke densities in the mountains are lower than that in the highly urbanized plains.

Keywords:

• Lighting location system
• detection efficiency
• cloud-to-ground lightning
• flash/stroke density

1. Introduction

Lightning location systems (LLSs) capture natural lightning electromagnetic signals and consequently locate lightning sources. LLS observations can be used to derive and retrieve lightning parameters that indicate lighting characteristics, e.g. the cloud-to-ground (CG) lightning flash density, the positive/negative CG flash/stroke density, the percent of positive/negative flashes, the first stroke negative and positive peak currents, and the multiplicity for negative and positive flashes (Orville et al. 2002; 2011; Rudlosky and Fuelberg 2010; Taszarek et al. 2015). These types of data have been widely used in weather forecasting, severe weather warning systems and climatological analysis, e.g. assimilation in weather forecasting models (e.g. Fierro et al. 2011; Fierro and Reisner 2011), synoptic analysis of thunderstorms (e.g. Changnon 1993; Cummins and Murphy 2009), presentations of lightning characteristics (e.g. Schulz et al. 2005; Drüe et al. 2007; Mäkelä et al. 2010), lightning climatology and climate pattern recognition (e.g. Antonescu and Burcea 2010; Williams 2005; Taszarek et al. 2015; Etherington and Perry 2017), and lightning risk assessment (Hu et al. 2014). However, a serious problem emerges when a LLS provides an inhomogeneous detection efficiency (DE) within the detectable distance or even a very low DE at the edges of the detectable distance (Naccarato and Pinto 2009), which greatly restricts the applicability of the LLS data (Hu et al. 2014; Yao et al. 2012).

DE is determined by the performance and sensitivity of the sensors, the sensor network geometry, the underlying ground conductivity (Schütte et al. 1988; Naccarato and Pinto 2009; Mäkelä et al. 2010), etc. The unevenly distributed DEs contribute to LLS data uncertainty and ruin data validity (Williams 2005; Yao et al. 2012; Hu et al. 2014). Thus, it is desirable to further develop algorithms to estimate the DE, which can be applied for correcting the lightning flash/stroke density derived from the LLS data.

Moreover, such an algorithm is useful in LLS performance evaluation and network deployment validation and would help to determine optimal detector deployment (Cummins et al. 1998; Idone et al. 1998; Schulz et al. 2005). This methodology was originally intended to test network deployment optimization and network performance and eventually became applicable to improve LLS data quality. Schütte et al. (1987, 1988) suggested that the lightning impulse signal strength (LISS) weakens linearly with propagation distance and confirmed the Weibull probability distribution of LISSs received by a single station within a detectable distance. Naccarato and Pinto (2009) estimated the DE of the Brazilian lightning detection network (BrasilDAT) based on the individual sensor DE probability functions, which were derived from a large amount of CG stroke data provided by the network.

Methodologically, Naccarato and Pinto  (2009) estimated the BrasilDAT DE based on the sensor individual DE probability functions, which are derived from a large amount of CG stroke data provided by the network considering different distances from the sensor and specific peak current ranges, and minimized the unrealistic enhancement of the DE closer to the network boundaries. Their method can be briefly described as: (1) One smaller region (450 × 450 km) region was chosen as the center quad (QUAD) area, which was used to assess the reference cumulative peak current distribution (PCD); (2) The coverage area of the LDN was divided into blocks and the number of CG lightning flashes was computed for each block. The PCD for each cell was computed individually and adjusted to the reference PCD. The ratio between the adjusted and the reference PCD at the 0 kA axis corresponded to the RFDE of the network in that cell. Schütte et al. (1988) assumed that the distribution of lightning signal strengths is that of Weibull distribution, and the acceptance of each sensor in the network could be the cumulative probability of detected signal strength range. The damping model for signal propagation, considering the acceptance/distance relationship for different ground conductivities was also taken into account (Schütte et al. 1987). Then, the acceptance of lightning detectors can be calculated knowing the lightning signal-strength distribution, the receiver thresholds, and the damping conditions for the propagation of lightning pulses.

These applicable approaches obviously possess reference values in the aftermath of DE estimation. Note that the precedent work can be cited as a referential source; however, it is inappropriate to simply copy the methods of earlier researchers if considering different data resources and circumstances. For example, the Weibull distribution of LISSs identified by Schütte et al. (1987, 1988) only reflected the lightning characteristics observed by the European LLS network at that time. However, our study revealed that the probability distribution of the LISSs detected by the ADTD sensors deployed around Beijing, China, would fit the generalized extreme value (GEV) distribution, rather than the Weibull distribution. This finding demonstrates the necessity of modifying the precedent methodology to estimate single-station acceptance.

Moreover, the DE and location accuracy of a LLS can be influenced by the lightning location mode, which is defined as the number of detectors synchronously receiving the lightning electromagnetic signal in one lightning location. The algorithm should consider the multi-detector location mode related to the number of detectors in the lightning location. Usually, at least two magnetic direction finder (MDF) detectors and three time-of-arrival (TOA) detectors are needed in each lightning location (Schütte et al. 1988). As observed by a LLS network with newly upgraded IMPACT detectors (combined MDF and TOA technology), a lightning source can be located with signals synchronously detected by two or more such detectors (Bourscheidt et al. 2012). With multi-detector location modes, here, an iterative algorithm is introduced to derive the DE from single-station acceptance.

Underlying surface conductivity plays an important role in damping lightning signals and eventually influences the DE (Honma et al. 1998; Rudlosky and Fuelberg 2010). Schütte et al. (1987) suggested that land surface damping influences both the DE and location accuracy and that the TOA network appears more sensitive to land surface damping. At this point, the algorithm should consider the damping effect on signal strength, which is induced by land surface conductivity.

The algorithm was constructed by following the above mentioned fundamentals. In this methodology, the LISSs were calculated using a linear signal propagation model suggested by Cummins et al. (1998). The damping coefficients in a grid system were estimated by addressing the land surface conductivity of specific land use and land cover (LULC) types. After correcting the signal strength magnitudes of the damping factor, we identified the GEV probability distribution of the LISSs. Based on this probability density function (PDF), we estimated the single-station acceptance of lightning source signals and ultimately calculated the grid DEs using an iterative algorithm. We corrected the CG lightning stroke density for the DEs and found an improvement in the estimates of the CG lightning stroke density in comparison with the uncorrected density.

Admittedly, the uncertainty cannot be reduced completely through DE estimation because the DE estimation fundamentals involve simulating signal attenuation propagated over certain distances and dampening due to land surface conductivity. The improvement in the LLS data quality (e.g., DE and location accuracy) basically depends on the implementation of optimal detector deployments and detector upgrading, and DE estimates are critical in promoting the applicability of the LLS data.

2. Data description

The LLS data collected from 2007–2016 by the ADTD (Advanced TOA and Direction system, where TOA denotes time-of-arrival) network deployed by the China Meteorology Administration (CMA) were used to (1) identify the probability distribution of the LISSs and (2) derive the CG flash/stroke density, which is corrected for the DE. These data include the time, location, amperage and polarity of the CG lightning strokes.

The ADTD consists of more than 301 sensors (as of March 2011) in China (Yao et al. 2012). In Beijing, 9-14 ADTD-1 sensors (improved IMPACT (combined MDF and TOA) sensors) can detect 1–450 kHz (very low-frequency (VLF)/low-frequency (LF) band) lightning sources (Figure 1). The ADTD-1 sensors use the combined MDF and time-of-arrival (TOA) method for position retrieval (Ma 2015). In this method, if a lightning source is only detected by two ADTD-1 sensors, the algorithm uses one TOA hyperbolic curve and two MDF vectors to retrieve the position. If a lightning source is detected by three sensors in a non-duplicate region, the TOA algorithm is used to retrieve the position directly. In contrast, if the TOA is the first to find a duplicate location, then the MDF method is used to identify the true location. If a lightning event is detected by four or more sensors, a TOA least squares method is used to find a more precise location. Thus, a lightning location detected by four or more sensors is more precise than that detected by fewer sensors. In the ADTD data, the percentage of lightning sources reported by four or more sensors relative to the total number of detected sources is 63.3%. Meanwhile, the ADTD-observed + CG and –CG lightning peak currents are in the ranges of 0.08 kA to 995.9 kA and 0.258 kA to 992.6 kA, respectively (Figure 1). Ma (2015) concluded that the peak current estimates between 10 kA and 100 kA could be relatively consistent with ground truth currents, while the estimates (<10 kA or >100 kA) are 23% lower than ground truth peak current measurements.

Figure 1. Histogram of the peak current probability density of (a) the –CG lightning, (b) +CG lightning and + IC lightning identified from their peak currents less than 15 kA and (c) the total CG lightning.

The ADTD sensor manufacturers claim that the DE of the ADTD sensors could be 90% at distances between 300 and 600 km, with a median location accuracy error of 1 km. However, only the flash DE can be 90%, and the stroke detection efficiency (SDE) is lower. The first stroke peak current in a multiple-stroke CG flash can be greater than twice its subsequent stroke peak current (Rakov and Uman 1990). Thus, the sensors can capture the first larger peak stroke but miss the weaker subsequent stroke (Rudlosky and Fuelberg 2010). Moreover, some weak CG strokes (including single-stroke CG flashes) cannot be detected due to the signal attenuation induced by long-distance propagation and damping effects in lower conductivity mountainous regions (Schütte et al. 1988).

We estimated the SDEs of the ADTD in a grid system (1 km × 1 km, see Figure 2) and corrected the lightning stroke density using the SDE. The SDE estimates approximate those of the U.S. NLDN (National Lightning Detection Network) in 1998, which was reported to be 62% (Idone et al. 1998). Hence, the DE level of the ADTD is equivalent to that of the NLDN, at least in 1998, suggesting that a considerable potential for improvement remains in terms of network upgrades and that DE estimation is necessary for promoting the applicability of the ADTD data.

Figure 2. Fitted probability density function curves for (a) the Weibull (red line), (b) the log normal (light yellow line), (c) the extreme value (deep yellow line) and (d) the GEV (blue line) distributions of the LISSs with the histogram of the sample data.

3. Method

The probability distribution of the LISSs is critical in estimating single-station acceptance. Being prior probability distribution, it can be deduced using samples of lightning location data. The likely probability distribution of these data would be that of either the Weibull, generalized extreme value (GEV), lognormal (LN), extreme value (EV) or gamma distribution (Schütte et al. 1988; Evans et al. 1993). Hereinafter, the GEV distribution is primarily introduced (Refer to the literature of Evans et al. (1993) to see the fundamentals of the other distributions). As we found in the latter, the GEV would provide the fittest probability distribution of the LISSs in the Beijing region.

3.1. GEV distribution

The GEV distribution has a cumulative distribution function (CDF) (1) $$F\left(x\right|\mu ,\sigma ,\xi )=\hspace{0.17em}\text{exp}\hspace{0.17em}\left\{-\left[1+\xi {\left(\frac{x-\mu }{\sigma }\right)}^{-1/\xi }\right]\right\}$$(1) for $1+\xi (x-\mu )/\sigma >0,$ where $\mu \in R$ is the location parameter, $\sigma >0$ is the scale parameter and $\xi \in R$ is the shape parameter. Thus, for $\xi >0,$ the CDF is only valid when $x>\mu -\sigma /\xi ;$ for $\xi <0,$ the CDF is valid and has a short-tailed Gumbel distribution when $x<\mu +\sigma /(-\xi ).$ For $\xi =0,$ the CDF is formally undefined and can be replaced by the result obtained by taking the limit as $\xi \to 0.$ Then, (2) $$F\left(x\right|\mu ,\sigma ,0)=\hspace{0.17em}\text{exp}\hspace{0.17em}\left[-\hspace{0.17em}\text{exp}\hspace{0.17em}\left(-\frac{x-\mu }{\sigma }\right)\right]$$(2) with no restriction on x.

The density distribution function is consequently (3) $$f\left(x\right|\mu ,\sigma ,\xi )=\frac{1}{\sigma }\left[1+\xi {\left(\frac{x-\mu }{\sigma }\right)}^{(-1/\xi )-1}\right]\hspace{0.17em}\text{exp}\hspace{0.17em}\left\{-{\left[1+\xi \left(\frac{x-\mu }{\sigma }\right)\right]}^{-1/\xi }\right\}$$(3)

When $\xi =0,$ the density is positive over the entire real line and equal to the Weibull distribution: (4) $$f\left(x\right|\mu ,\sigma ,0)=\frac{1}{\sigma }\text{exp}\hspace{0.17em}\left[-\left(\frac{x-\mu }{\sigma }\right)\right]\hspace{0.17em}\text{exp}\hspace{0.17em}\left\{-\hspace{0.17em}\text{exp}\hspace{0.17em}\left[\left(-\frac{x-\mu }{\sigma }\right)\right]\right\}$$(4)

3.2. Fundamentals of the GEV distribution applied for estimating single-station acceptance

Consider a normalized (for example, to a distance of 100 km), undamped reference signal-strength distribution with the frequency function f(s₀) and an expression for the signal strength dependence on the distance s₀ =D⁻¹(s), where s₀ is the signal strength normalized to a standard distance r_0. The simplest case for an undamped radiation field is the inverse power law s= s₀r₀/r, expressing the signal-strength attenuation with distance r, where s is the signal-strength. The acceptance, A(r), defined as the percentage of detected cloud-to-ground lightning flashes (Schütte et al. 1987), can be calculated as a function of lightning signal strength distribution f(s₀), the distance r and the lower and upper signal-strength thresholds s_min and s_max, respectively: (5) $$\begin{array}{c}A\left(r\right)={\int }_{s{D}_{\min }}^{s{D}_{\max }}f\left(s_0\right)\mathrm{d}{s}_0,s_0=D^{-1}\left(s\right),\mathrm{d}{s}_0=D^{-1^{\prime }}\left(s\right)\mathrm{d}s\\ ={\int }_{s{D}_{\min }}^{s{D}_{\max }}f\left[D^{-1}\right(s\left)\right]D^{-1^{\prime }}\mathrm{d}s=F\left[D^{-1}\right(s\left)\right]|\begin{array}{c}{s}_{\max }\\ s_{\min }\end{array}\end{array}$$(5)

For undamped propagation, we can obtain (6) $$A\left(r\right)=F(\mathit{sr}/r_0)|\begin{array}{c}{s}_{\max }\\ s_{\min }\end{array}$$(6)

In the case of the GEV distribution, combined with Eqs. (1)(1) $$F\left(x\right|\mu ,\sigma ,\xi )=\hspace{0.17em}\text{exp}\hspace{0.17em}\left\{-\left[1+\xi {\left(\frac{x-\mu }{\sigma }\right)}^{-1/\xi }\right]\right\}$$(1) and (6)(6) $$A\left(r\right)=F(\mathit{sr}/r_0)|\begin{array}{c}{s}_{\max }\\ s_{\min }\end{array}$$(6) , the acceptance is given by (7) $$A\left(r\right)=\{\begin{array}{l}0\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cccc}& & & \end{array}& & & \end{array}& & & \end{array}& & & \end{array}& r<c{r}_0/s_{\max }& & \end{array}& & & \end{array}\\ \hspace{0.17em}\text{exp}\hspace{0.17em}\left\{\left[-[1+\xi {\left(\frac{s_{\max }r/r_0-\mu }{\sigma }\right)}^{-1/\xi }\right]\right\}\begin{array}{cccc}\begin{array}{c}\end{array}& c{r}_0/s_{\max }<r\le c{r}_0/s_{\min }& & \end{array}\end{array}$$(7)

where c is the signal bias and can be replaced by the minimum value of the signal samples.

The effective radius, $\rho ,$ describing the properties of a lightning counter or direction finder can be defined as (8) $$\pi {\rho }^2=2\pi {\int }_0^{\infty }A\left(r\right)r\mathrm{d}r$$(8)

This definition means that an ideal lightning counter accepts all lightning signal s with $r\le \rho$ and none with $r>\rho .$

When the signal strength is GEV distributed and the counter has as a sufficient dynamic range (s_min ≪ s_max), the effective radius can be calculated analytically: (9) $$\begin{array}{c}{\rho }^2=2{\int }_{\mu r_0/s_{\max }}^{\mu r_0/s_{\min }}r\cdot \hspace{0.17em}\text{exp}\hspace{0.17em}\left\{-[1+\xi {\left(\frac{s_{\max }r/r_0-\mu }{\sigma }\right)}^{-1/\xi }]\right\}\mathrm{d}r\\ =2{\int }_{c{r}_0/s_{\min }}^{\infty }r\cdot \hspace{0.17em}\text{exp}\hspace{0.17em}\left\{-\left[1+\xi {\left(\frac{s_{\max }r/r_0-\mu }{\sigma }\right)}^{-1/\xi }\right]\right\}\mathrm{d}r\\ -2{\int }_{c{r}_0/s_{\max }}^{\infty }r\cdot \hspace{0.17em}\text{exp}\hspace{0.17em}\left\{-\left[1+\xi {\left(\frac{s_{\min }r/r_0-\mu }{\sigma }\right)}^{-1/\xi }\right]\right\}\mathrm{d}r.\end{array}$$(9)

By the substitution $$z=-\left[1+\xi {\left(\frac{s\cdot r/r_0-\mu }{\sigma }\right)}^{-1/\xi }\right]$$ the integral (10) $${\int }_{c{r}_0/s}^{\infty }r\cdot \hspace{0.17em}\text{exp}\hspace{0.17em}\left\{-\left[1+\xi {\left(\frac{s_{}r/r_0-\mu }{\sigma }\right)}^{-1/\xi }\right]\right\}\mathrm{d}r$$(10) becomes (11) $${\int }_0^{\infty }\left(\frac{\sigma }{\xi \cdot s}{z}^{-\frac{1}{\xi }}+\frac{\sigma +\xi \cdot \mu \cdot r_0}{\xi \cdot s}\right)e^{-z}\mathrm{d}z=\frac{1}{s}\left[\frac{\sigma }{\xi }\Gamma \left(\frac{-1}{\xi }\right)+{\int }_0^{\infty }\left(\frac{\sigma +\xi \cdot \mu \cdot r_0}{\xi }\right)e^{-z}\mathrm{d}z\right]$$(11)

Thus, the substitution $B={\int }_0^{\infty }\left(\frac{\sigma +\xi \cdot \mu \cdot r_0}{\xi }\right)e^{-z}\mathrm{d}z$ does not include variable s and is integrable. The equation to calculate the effective radius can then be formulated as (12) $${\rho }^2=\left(1/s_{\min }^2-1/s_{\max }^2\right)\cdot \left[\frac{\sigma }{\xi }\Gamma \left(\frac{-1}{\xi }\right)+B\right]$$(12)

Because S_min is far less than S_max, here, S_min can greatly influence the effective radius, while S_max has little influence on the effective radius.

3.3. Damping modes

The LISS, s, and return stroke peak current, $I_{\mathit{peak}},$ are in a linear relationship: (13) $$s=k\cdot {\left(\frac{r}{r_0}\right)}^p.I_{\mathit{peak}}$$(13) where Eq. (13)(13) $$s=k\cdot {\left(\frac{r}{r_0}\right)}^p.I_{\mathit{peak}}$$(13) is calculated using the coefficient k and the signal strength, has a non-parametric unit (a.u.). The empirical value of the coefficient k is 100/3, equivalent to a signal strength of 150 (a.u.) produced by a return stroke peak current of 45 kA over a distance of r₀=100 km. The coefficient p is also an empirical value, and here, it is set to 1. Eq. (13)(13) $$s=k\cdot {\left(\frac{r}{r_0}\right)}^p.I_{\mathit{peak}}$$(13) indicates a linear relationship between signal attenuation and the signal propagating distance. Some applications prefer a non-linear relationship; for example, Cummins et al. (1998) attempted to determine an unbiased non-linear relationship to calculate the signal strength detected by the NLDN, in which p was set to 1.13. However, our application found that the non-linear and linear relationships were similar (not given) for our results.

The lightning electromagnetic signal will be attenuated when propagated over land and/or the ocean surface, and this attenuation is determined by the underlying surface conductivity. Based on a combination of experiments and calculation, Schütte et al. (1988) drew an analytic approximation for the damping made by curve-fitting, which can be cited in estimates of signal attenuation.

Here, we introduce the parameter (14) $$\gamma =r/\theta$$(14) where r is the distance and $\theta$ is the ground conductivity. The damping, Dm, is expressed as (Schütte et al. 1988) (15) $$\mathit{Dm}(\gamma )=1.652-0.0491\hspace{0.17em}\text{ln}\hspace{0.17em}(\gamma /{\gamma }_0)$$(15) in which $\gamma >{10}^6{m}^2\Omega ,$ ${\gamma }_0=1m^2\Omega .$ The damping factor Dm is multiplied by the 1/r propagation law for the radiation fields in order to calculate the signal-strength (16) $$s\left(r\right)=\mathit{Dm}(r/\theta )s_0{r}_0/r$$(16)

Therefore, in the damping modes, the single-station acceptance can be calculated as (17) $$A\left(r\right)=\{\begin{array}{l}0\begin{array}{cc}& \end{array}\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cccc}& & & \end{array}& & & \end{array}& & & \end{array}& & & \end{array}\begin{array}{cc}\begin{array}{cc}& \end{array}& \end{array}& r<c{r}_0/s_{\max }& & \end{array}& & & \end{array}\\ \hspace{0.17em}\text{exp}\hspace{0.17em}\left\{-\left[1+\xi {\left(\frac{D(r/\sigma )s_{\max }r/r_0-\mu }{\sigma }\right)}^{-1/\xi }\right]\right\}\begin{array}{cccc}\begin{array}{ccc}& & \end{array}& c{r}_0/s_{\max }<r\le c{r}_0/s_{\min }& & \end{array}\end{array}$$(17)

3.4. Iterative algorithm for DE estimation

After the acceptance of each detector in a lightning location network is confirmed, the DE of the network can be estimated using an iterative algorithm. The numbers of detectors dedicated to locate lightning flashes/strokes are determined by lightning locating methods, e.g. MDF and TOA. The MDF method requires at least two detectors synchronously detecting one lightning return stroke (Schütte et al. 1987), while the TOA method requires four detectors since the signals received by two additional TOA detectors are needed to eliminate the falsely duplicated locations detected by the other two detectors. The IMPACT detector networks locate a lightning stroke using at least two detectors. Thus, an iterative algorithm is needed in the deduction of DEs, especially with respect to lightning locations that require three or more detectors in LLS network coverage.

The algorithm is expressed as follows.

1. This iterative algorithm is designed for DE estimation in two-detector networks: (18) $$F(2,2)=A\left(r_1\right)·A\left(r_2\right)$$(18) where F(2,2) is the DE of the two-detector location mode in a two-detector network and A(r₁) and A(r₂) are acceptances of the detectors r₁ and r₂, respectively.

2. In three-detector networks (19) $$F(3,2)=A\left(r_1\right)·A\left(r_2\right)·A\left(r_3\right)+A\left(r_1\right)·A\left(r_2\right)·R\left(r_3\right)+A\left(r_1\right)·R\left(r_2\right)·A\left(r_2\right)+R\left(r_1\right)·A\left(r_2\right)·A\left(r_2\right)$$(19) where F(3,2) is the DE of the two-detector location mode in a three-detector network, A(r₃) is the acceptance of detector r₃, and R(r₁), R(r₂) and R(r₃) are the non-acceptances of detectors r₁, r₂ and r₃, respectively.

3. In an n-detector network (n > 3), (20) $$F(n,m)=R\left(r_n\right)·F(n-1,m)+A\left(r_n\right)·F(n,m-1)$$(20) where F(n,m) is the DE of the m-detector location mode in an n-detector network.

In Eq. (20)(20) $$F(n,m)=R\left(r_n\right)·F(n-1,m)+A\left(r_n\right)·F(n,m-1)$$(20) , (21) $$R\left(r_i\right)=1-A\left(r_i\right),i=1,2,\dots ,n$$(21)

4. Test and case study

4.1. Identifying the probability distribution and estimating parameters

The maximum likelihood method allows the determination of the most effective estimators, and the sample data histogram is helpful in identifying the probability distribution (Bishop 2006; Kotowski and Kaźmierczak 2013). The sample data of 240804 lightning CG strokes detected by the ADTD network from 2007 to 2016 were used to determine the probability distribution of the LISSs in Beijing and the parameters of the probability density function (PDF). To obtain the fittest probability distribution, which can be used to calculate the probability above a certain signal strength threshold, we estimated the parameters of the Weibull, the lognormal, the extreme value, the gamma, and the generalized extreme value (GEV) distributions from the sample data. In a comparison of the fitness of these probability distribution curves and the sample data histogram, the GEV curve is identified as the best fitted (Figure 2). The parameters of these probability distribution functions are the maximum likelihood estimates (MLE) computed by the machinery learning iterative algorithm. The fitness of these PDFs can be indicated by the negative of the log-likelihood for these distributions (Table 3), which the magnitude is the smaller and the fitter. Obviously, the negative of the log-likelihood for the GEV distribution is the smallest. Therefore, the GEV probability distribution is regarded as the distribution with the best fit and is used in deducing single-station acceptance. The probability distribution of the LISSs is empirically established in certain regions, and it should be properly predefined before applying the results to deduce acceptance.

Table 3. Parameters estimated using maximum likelihood method and the negative of the log-likelihood nlogL for the listed distributions.

	The Weibull	The log normal	The extreme value	The GEV
Parameters of MLE	The scale 302.826, and the shape 0.982	The mean 5.284, and the standard deviation 0.816	The location 1831.2, and the scale 12596.1	The shape 0.508, the scale 97.450, and the location 143.282
nlogL	Infinite	862,290.272	1401,342.58	856,512.42

4.2. DE estimates

The ADTD system observes the lightning parameters of the CG lightning strokes, e.g. the peak current, current intense slope, and latitude and longitude of lightning stroke position. However, the CG lightning strokes were not grouped into lightning flashes by origin. We estimated the stroke DEs rather than the flash DEs.

Usually, if the first return stroke or one of its subsequent return strokes in a multiple lightning flash is located, that flash can be counted as detected (Naccarato and Pinto 2009). Therefore, the stroke DEs would be lower than the flash DEs. Naccarato and Pinto (2009) confirmed that the ratio of the stroke DEs to the flash DEs in the Paraíba Valley, Brazil, should be 1.6, implying a larger magnitude of stroke DE.

Using the software of Matlab 2012b, we deduced the parameters of the GEV distribution of the LISSs using the sampled lightning electromagnetic signal data (Figure 1): $\sigma$=97.449, $\xi$=0.5083, $\mu$=143.2815. Next, the single-station signal acceptance was deduced using Eq. (17)(17) $$A\left(r\right)=\{\begin{array}{l}0\begin{array}{cc}& \end{array}\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cccc}& & & \end{array}& & & \end{array}& & & \end{array}& & & \end{array}\begin{array}{cc}\begin{array}{cc}& \end{array}& \end{array}& r<c{r}_0/s_{\max }& & \end{array}& & & \end{array}\\ \hspace{0.17em}\text{exp}\hspace{0.17em}\left\{-\left[1+\xi {\left(\frac{D(r/\sigma )s_{\max }r/r_0-\mu }{\sigma }\right)}^{-1/\xi }\right]\right\}\begin{array}{cccc}\begin{array}{ccc}& & \end{array}& c{r}_0/s_{\max }<r\le c{r}_0/s_{\min }& & \end{array}\end{array}$$(17) .

When dealing with the signal attenuation induced by underlying land surface conductivity, we used the soil-type distribution in the Beijing region to estimate the damping coefficient in the grids. The Beijing Bureau of Geology and Mineral Resources (BBGMR) provided the soil-type distribution maps, which had been digitized as GIS polygons of the soil type (not displayed here). These polygons can be overlapped with the propagating paths of electromagnetic lightning signals to obtain the path segments over different soil-type polygons, which exhibit different conductivities and attenuation signals. The predefined soil-type conductivities are listed in Table 1. Additionally, we generalized the LULC (Land use and land cover) as regions of mountains and plains where we lacked the detailed soil-type distribution. Accounting for the most extensively distributed bare rocks, gravels, and sands in the mountainous areas, which exhibit lower conductivity, the mountainous conductivity magnitude is simplified to be 1500 Ω m. On the other hand, the conductivity in the plains, which are mostly composed of loamy soil, is 60 Ω m. Therefore, using Eqs. (14)(14) $$\gamma =r/\theta$$(14) to (17)(17) $$A\left(r\right)=\{\begin{array}{l}0\begin{array}{cc}& \end{array}\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cccc}& & & \end{array}& & & \end{array}& & & \end{array}& & & \end{array}\begin{array}{cc}\begin{array}{cc}& \end{array}& \end{array}& r<c{r}_0/s_{\max }& & \end{array}& & & \end{array}\\ \hspace{0.17em}\text{exp}\hspace{0.17em}\left\{-\left[1+\xi {\left(\frac{D(r/\sigma )s_{\max }r/r_0-\mu }{\sigma }\right)}^{-1/\xi }\right]\right\}\begin{array}{cccc}\begin{array}{ccc}& & \end{array}& c{r}_0/s_{\max }<r\le c{r}_0/s_{\min }& & \end{array}\end{array}$$(17) , we can estimate the signal strength influenced by the damping effects induced by the land surface conductivity.

Table 1. Lists of underlying soils conductivities.

Types	Subtypes	Conductivity (Ω m)	Conductivity ranges of different conditions (Ω m)
			Moist areas	Dry areas	Saline and alkaline underground water
Soils	Clay soils	10	5–20	10–100	3–10
	Clay loams	20	10–30	50–300	3–30
	Heavy loamy texture soils	50	30–100	50–300	10–30
	Moderate loamy texture soils	60	30–100	50–300	10–30
	Light loamy texture soils	100	30–300	80–1000	10–80
	Sandy loamy texture soils	200	100–200	250	30
	Loose sandy soils	300	100–1000	1000	30–100
	Sandy soils		300	/	/
Sands	Sands	1000	250–1000	1000–2500	/
Gravels	Gravels and pebble beaches	5000	250–1000	1000–2500	/
	Bare rocks	5000	/	/	/
Others	Urban areas	2000	/	/	/
	Mountains	1500	/	/	/
	Plains	60	30–100	50–300	10–30
	Rivers and reservoirs	5	/	/	/

Using Eq. (17)(17) $$A\left(r\right)=\{\begin{array}{l}0\begin{array}{cc}& \end{array}\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cccc}& & & \end{array}& & & \end{array}& & & \end{array}& & & \end{array}\begin{array}{cc}\begin{array}{cc}& \end{array}& \end{array}& r<c{r}_0/s_{\max }& & \end{array}& & & \end{array}\\ \hspace{0.17em}\text{exp}\hspace{0.17em}\left\{-\left[1+\xi {\left(\frac{D(r/\sigma )s_{\max }r/r_0-\mu }{\sigma }\right)}^{-1/\xi }\right]\right\}\begin{array}{cccc}\begin{array}{ccc}& & \end{array}& c{r}_0/s_{\max }<r\le c{r}_0/s_{\min }& & \end{array}\end{array}$$(17) , we calculated the single-station acceptances A(r) of approximately 11 ADTD detectors around Beijing. The distribution maps indicate that the acceptances in areas closer to the detectors are clearly higher and that the acceptances rapidly decrease toward the mountainous areas (see Figure 3).

Figure 3. Distributions of the single-station acceptances of the detectors in (a) Zhangjiakou, (b) Guangxiangtai, and (c) Tanggu.

Subsequently, we deduced the ADTD network DE (Figure 4) using the iterative algorithm based on single-station acceptances. According to the recordings of the location modes in the ADTD data detected from 2007 to 2016, these CG lightning strokes were located after their signal sources were received synchronously by 2 to 7 detectors (see Table 2 for the proportions of each location mode, indicating number of detectors dedicated to locate one lightning flash/stroke). Note that only 11 ADTD detectors in the Beijing region had been selected to estimate the DE for the case study, and thus, the results are only valid in this region. The results indicate that the DEs in the Beijing metropolitan area are relatively higher, almost all are above 60%, and the maximum is 78% (see Figure 4).

Figure 4. Deduced CG lightning stroke DE of the ADTD network around Beijing (a) not accounting for the land surface conductivity damping and (b) accounting for the land surface conductivity damping.

Table 2. Proportions of the ADTD location modes recorded in the lightning location data.

Location modes	2-sensor location	3-sensor location	4-sensor location	5-sensor location	6-sensor location	7-sensor location
Proportion	34.25%	2.45%	4.3%	19.18%	37.6%	2.22%

It is reasonable to expect that the CG lightning stroke DEs can reach 60% or higher. Based on a video observation, Idone et al. (1998) found that from 1994 to 1996, the in-site NLDN stroke DEs were most likely between 47% and 67%, while the flash DEs reached between 72% and 86%. Regarding the Vaisala reports, Orville et al. (2011) suggested that the flash DEs of the North American Lightning Detection Network (NALDN) in the contiguous U.S. and southern Canada should surpass 90%, while the overall stroke DEs would be between 60% and 80%.

According to the DE estimates of the ADTD network (Figure 4), in contrast to the low DEs in the mountainous northern and western regions, the southern plains exhibit a higher DE. Relatively exceptionally low DEs can be found in parts of the Beijing district. These areas are closer to the detector at Guangxiangtai station but far away from the other closest detectors (e.g., Zhangjiakou), meaning that the signal would propagate a long distance over mountains before being received by these detectors; as a consequence, the signal would be attenuated.

We estimated the DEs of the ADTD network by both accounting and not accounting for the damping effects induced by the underlying conductivity. In the comparison of the DEs deduced in the damping and non-damping modes, the DE deduced in the non-damping mode is clearly higher than that in the damping mode, and this discrepancy is exaggerated in mountainous areas (see Figure 4), where the lower conductivity leads to a greater attenuation of lightning electromagnetic signals.

4.3. CG lightning stroke density corrected for the deduced DEs

The CG lightning stroke densities (stroke/yr km²) were derived from the ADTD lightning location data. The derivation was performed on the 1 $\times$ 1 km resolution grid, which is 182 km N–S and 181 km E–W and covers the entire Beijing district. The GIS “contain” operator is applied for counting the lightning CG stroke numbers in the grid cells. Then, the CG lightning stroke density N_grid can be corrected for the grid DEs F_grid, and the correction formula is (22) $${N^{\prime }}_{\text{grid}}=N_{\text{grid}}/F_{\text{grid}}$$(22) where N′_grid (stroke/yr km²) is the corrected CG lightning stroke density.

The derivation indicated that the CG lightning stroke densities in the mountainous northern and western regions are clearly lower than that of the other areas (see Figure 5). This characteristic is similar to that of the DE estimates. In the mountainous regions, the elevated terrains and low conductivity of the underlying surface, as well as rough and non-flat ground configuration (Li et al. 2017), have a strong influence on signal attenuation (Cummins and Murphy 2009), which in turn reduces the DE and location accuracy of the ADTD network. However, this result does not convincingly prove that the low CG lightning stroke density in mountainous areas is determined by the low DEs.

Figure 5. Distributions of the ADTD-derived (a) non-corrected CG lightning stroke density and (b) CG lighting stroke density corrected for the DE.

After the CG lightning stroke density was corrected for the DE, the densities in the mountainous northern and western regions increased slightly but remained lower than those of the other areas (see Figure 5). Higher densities are still present in the mountainous areas of eastern and southwestern Beijing and in the plains.

A high CG stroke density zone appears in the urban area and its downwind areas (see A in Figure 5b), in agreement with the hypothesis that an increased urban land-surface roughness could intensify thunderstorm activities (Shepherd et al. 2002; Rose et al. 2008; Stallins et al. 2006; Stallins and Rose 2008; Hu et al. 2014; Hu 2015; Sun and Shu 2007).

Additionally, the aerosols serve as cloud condensation nuclei (CCN) and have a substantial effect on cloud properties and the initiation of precipitation (Steiger et al. 2002; Naccarato et al. 2003). The higher CG lighting strikes predicted in the plains agree with the hypothesis of the urban aerosol effects on lightning enhancement since the plains usually have more polluted air than the mountainous areas (Tie et al. 2015). Moreover, an urban heat island (UHI) can trigger an urban-to-rural thermal circulation and force a convection, which produces enhanced thunderstorms in the metropolitan regions (Jiang and Liu 2006; Sun and Shu 2007).

Therefore, it is justifiable that the corrected CG lightning stroke densities would be higher in the plains, almost twice (mostly above 2.0 stroke/yr km²) the density in the northern and northwestern areas. Therefore, the urbanization in the plains plays an important role in inducing thunderstorm enhancement that results in a higher CG lightning stroke density without accounting for the higher DEs in the plains.

5. Conclusion

The DE is determined by the LLS performance and sensitivity of the detectors, detector network geometry, and damping factors induced by the underlying land surface conductivity. For this reason, DEs vary due to the underlying land surface and become very low in some regions. This situation leads to a discrepancy between the derived and actual lightning flash/stroke density. This characteristic discounts the applicability of the LLS data. Therefore, we developed an algorithm to estimate DEs which can be used for lightning flash/stroke density correction. The method and its implementation in the CG lightning stroke density correction can be summarized as follows:

1. It is critical to identify a suitable probability distribution of the LISSs, which is determined for specific regions. In comparison with the lognormal, the Weibull, the extreme value (EV) and the GEV probability distributions, we found that the GEV distribution provides the best fit to the histogram of the sample data.

2. Based on formula derivations, it was verified that the GEV probability distribution fit would be suitable in estimating the DEs and that the single-station effective radius and acceptance would be determined by the minimum detectable signal strength (S_min).

3. To address the damping effects induced by the inconsistent underlying land surface conductivity, we estimated the single-station acceptance in the damping modes with respect to the GEV probability distribution of the LISSs.

4. The DEs of the ADTD network were deduced based on an iterative algorithm using single-station acceptances. According to the LLS location data in Beijing, the number of detectors required to locate the CG lightning strokes was generally between 2 and 7. Thus, the complexity of the iterative algorithm was low, and a convergent solution was found.

5. The damping effects induced by underlying land surface conductivity were considered in the DE estimates. Thus, the estimate results exhibit that the DEs in the mountains are clearly lower than that in the plains. This finding may be due to the lightning electromagnetic signals propagating over mountainous areas with low-conductivity underlying land surfaces, as well as rough and non-flat ground configuration, which considerably attenuates the signal.

6. The CG lightning stroke density corrections for the DEs reflected that although the densities in the mountainous northern and western areas increase to some degree, they are still lower than those in the highly urbanized plains. This phenomenon can be related to the increased urban land-surface roughness, the aerosols and the UHI effects on the thermodynamic activity in urban areas.

The DE estimates can be effectively used for correcting the CG lightning flash/stroke density derived from the LLS data. However, we suggest that the location error correction and uncertainty reduction in the lightning detection should be continued by adding new detectors and/or network upgrading and optimizing the network sensor deployment, which would vitally enhance the LLS performance and sensor sensitivity and achieve the goal of DE improvement.

Acknowledgements

This study has been supported by the Climate Change Research Projects of China Meteorological Administration (Project: CCSF201717) and the Research Projects of Institute of Urban Meteorology, CMA (Project: IUMKY201716).