Computation of Lightning Current from Electric Field Based on Laplace Transform and Deconvolution Method

: A method for computation of the lightning channel base current from the corresponding vertical component of lightning electric ﬁeld was presented. The algorithm was developed by applying Laplace transform. The lightning current was estimated from its deconvolution with a special transfer function. The transfer function includes information about geometry and physical properties of entire lightning impulse generation system. The method was veriﬁed for a Heidler-type base current and a MTLL model of its propagation within the lightning channel. Research was done for close, middle, and far distance to the lightning strike point. Optimum performance was obtained for the middle distance of several kilometers where the electrostatic, induction, and radiation components of the transfer function were of the same range. An analysis was done for input electric ﬁeld with and without noise superimposed on its time domain waveform. Relative uncertainties for the electric ﬁeld and calculated lightning channel base current were similar each other. The presented approach can substantially increase a number of lightning current parameters which can be identiﬁed on the basis of its electric ﬁeld signature. This method can be applied by the lightning location systems using preprocessing which increases the timing efﬁciency of the transfer function estimation.


Introduction
The estimation of the lightning current is one of the fundamental tasks of overvoltage protection [1]. Nowadays, information about the current peak value is reported by numerous databases of lightning locating systems (LLS) [2]. However, there is a lack of information about other parameters of the current at the lightning channel base. Information about steepness of the current is necessary for the evaluation of induced effects on nearby structures, especially in case of electrical power distribution systems [3,4]. The shape of the lightning channel base current waveform in time domain allows for the detection of a continuous and continuing current [5] being responsible for thermal destruction and fires [6]. Moreover, the computation of the lightning current from its electric field signatures recorded for natural lighting can be applied for investigation of physical processes within the lightning channel [7].
The lightning channel base current signature can be obtained directly and remotely. Direct measurements are based on electrical current sensors usually installed at the elevated objects. One of the most reliable studies was done by Berger and coworkers in Switzerland [8,9]. This database, consisting of 249 strokes, is a reference for numerous lightning protection standards [10]. Other tower-based measurements were done by Dellera in Italy [11], Geldenhuys in South Africa [12], Takami and Okabe in Japan [13], Visacro in Brazil [14], and Diendorfer in Austria [15]. Lightning channel base current parameters for triggered lightning were obtained at the NASA Kennedy Space Center and Alabama by Fisher et al. [16]. The continuation of rocked experiments was done at Camp Blanding, Florida, by Rakov et al. [17]. Schoeone et al. presented a statistical analysis of 206 current waveforms obtained in 46 events recorded from 1999 to 2004 [18]. The remote estimation of the lightning channel base current can be done from the lightning electromagnetic fields [19]. This group of methods, modified for the radiated field component, is commonly used by most of the lighting location systems. One of the first attempts to evaluate the current peak value from the corresponding electric field peak value was done theoretically by Uman and McLain. Rakov derived a similar empirical formula from triggered experiments [20]. Another approach to the estimation of the triggered lightning current peak value and the current steepness was done in 1998 by Willet et al. [21] using a semiempirical equation applied to the electric field data. Mallick et al. corrected Rakov formula on the basis of data obtained by a rocket-and-wire technique of lightning triggering [22]. Independently, Rachidi et al. estimated the base current peaks from LLS measurements [23]. Uman et al. [24] and Rachidi et al. [25] derived equations for the lightning current computation from far electromagnetic fields. In 2004, Shao et al. [26] presented a derivation of the theoretical formula relating radiated components of lightning electric and magnetic fields to the corresponding current in the lightning channel taking into account the nonconstant retarded time along the channel segment. Next, Thottappillil and Rakov [27] made a comment on Shao et al. [26], finding a sign error in the formula and concluding that the given method could not be applied to models with current discontinuity at the return stroke front. In response [28], Shao et al. agreed with the addressed remarks and presented an extension of their formula for the modified transmission line (MTLE) model. Delfino et al. found the time domain method of the base current estimation from magnetic field using a Volterra like integral equation [29]. Andreotti derived equations operating on the electric field in the frequency domain [30]. Lightning channel base current identification procedures for close and intermediate distances from the lightning channel were proposed by Izadi [31]. Yang et al. invented a time-domain algorithm for the lightning current estimation from magnetic field based on the deconvolution method [32].
In this paper, the author proposes a method for computation of the lightning channel base current from the corresponding vertical component of lightning electric field. The algorithm applies Laplace transform which leads to the deconvolution equation. This equation is further discretized in order to obtain iterative formula for the current waveform in time domain. The method is designed to work for most of the base current models and attenuation functions found in the literature [33]. The algorithm was verified for Heidler current source [34] and MTLL model of its propagation within the lightning channel [35]. Lightning current waveforms were derived from close, middle, and far electric fields. The performance of this method was tested for different numerical parameters of simulation.

Lightning Channel Base Current Computation Algorithm
The geometry of the lightning channel is defined in Figure 1. The lightning channel was assumed to be a vertical line of length H. An electromagnetic field was propagating over a perfect soil. Thus, the total field at observation point was a superposition of influences from a real lighting channel and its image. The electric field observed near ground surface was considered according to [36]: where E z is a vertical component of electric field, P(z') is the lightning current attenuation function, ε 0 is the vacuum permittivity, v is the current wave propagation speed, and c is the speed of light, ∆t = |z'|/v + R/c is the retardation time. The remaining parameters are defined in Figure 1. P(z') is dependent only on the height of a particular elementary current source. Therefore, the algorithm can be applied to all LCS and almost all DSC current source models [35]. The lightning channel was assumed to be a vertical line of length H. An electromagnetic field was propagating over a perfect soil. Thus, the total field at observation point was a superposition of influences from a real lighting channel and its image. The electric field observed near ground surface was considered according to [36]: where Ez is a vertical component of electric field, P(z') is the lightning current attenuation function, ε0 is the vacuum permittivity, v is the current wave propagation speed, and c is the speed of light, Δt = |z'|/v + R/c is the retardation time. The remaining parameters are defined in Figure 1. P(z') is dependent only on the height of a particular elementary current source. Therefore, the algorithm can be applied to all LCS and almost all DSC current source models [35]. Three auxiliary Functions (2)-(4) were defined to simplify further conversions:  Three auxiliary Functions (2)-(4) were defined to simplify further conversions: Functions (2)-(4) refer to the electrostatic, induction and radiation effects, respectively. The inclusion of 1/(4πε 0 ) into Functions (2)-(4) improved numerical stability by increasing the range of those functions over ten times.
Applying Laplace transform to Equation (1) taking into account Functions (2)-(4) gave:   Then, the inverse Laplace transform was applied to obtain time domain signals: An additional integration of electric field allowed for avoiding an indeterminate form during the inverse Laplace transform of its radiated component.
Next, the property (8) of Dirac delta was used to integrate through the channel length: The application of relationship (8) to Equation (7) allowed the computation of the inner integral of Formula (7): where z t ' is value of z' computed from ∆t = |z'|/v + R/c for the particular ∆t = t. This approach significantly reduced the complexity of numerical computations because the inner integral was expressed only by the product of auxiliary functions (2)-(4) with attenuation function P(z t ').
For simplicity, Equation (9) was further presented as: where w(t) is a transformation function which includes entire information about the geometry of the model, attenuation and retardation time. It is important for the efficiency of numerical computations that the transformation function w(t) can be obtained during preprocessing procedures. Further computation of the base current was done by the application of the convolution algorithm: Discretization of Equations (10) and (11) gave: where k is the number of k-th sample of electric field and ∆T is the sampling interval. The rearrangement of Equation (12) allowed to write: In the recurrent Formula (13), there appeared a division by the first element of vector w[r,z,t]. This element should differ from zero to avoid an indeterminate form. Unfortunately, it can be seen from relationship (9) that w[r,z,t] is zero for time varying from 0 up to retardation time R/c. To resolve this problem the electric field E z (r,z,t) given in Equation (13) and transformation function w(r,z,t) should be shifted by the retardation time R/c which was done in Equation (14). Therefore, time of the first sample of signals given in Formula (14) is Finally, the equation for the lightning channel base current was: There is an alternative matrix equation I B = W −1 E for the lightning channel base current which can be derived from Equation (12), where I B is the vector of the lightning channel base current, W −1 is the matrix based on transfer function, and E is the given vector of electric field samples. However, the iterative form of Equation (14) is more efficient for the computation of longer time intervals and for lower sampling times where the matrix equation can produce out of memory errors due to the large size of W −1 matrix.

Results and Discussion
At the preliminary stage of the analysis, the correctness of derivations given in Section 2 was verified by numerical simulation of lightning electric field. The vertical component of electric field was computed on the basis of Equation (9) and compared with results from similar analyses based on Equation (1) found in [19,33]. It was not necessary to place all obtained electric field waveforms in this paper because any significant differences were observed according to those well described and already available in the literature. Moreover, the scope of the paper was not the computation of electric field but the numerical estimation of current waveform at the lightning channel base. Several electric field signatures for MTLL lightning current attenuation model at different distances from the strike point are given in Figure 2. The verification procedure was done for several different lightning current attenuation models [33,35] and for the same lightning channel base current given by the Heidler formula [34]: The simulation parameters were I 01 = 9.9 kA, η = 0.845, τ 1 = 0.072 µs, τ 2 = 5 µs, I 02 = 7.5 kA, τ 3 = 100 µs, τ 4 = 6 µs, and the channel height H = 7 km, the current wave speed v = 1.3 × 10 8 m/s, the sampling time interval was ∆T = 10 ns, and the initial channel segment length was ∆z' = 2 m. The lightning channel was initially discretized to obtain the transformation function w(t). Then, w(t) was resampled to fit to the lighting channel base current sampling time. This sampling manner was necessary for the application of the recurrent Formula (12). The transformation function could also be computed directly for uniformly spaced time vector, as well. However, this method involves the additional time-consuming computation of z' vector from t = |z |/v + (z − z ) 2 + r 2 /c.  At the second step of the simulation the lightning channel base current derived from Equation (14) was compared with the assumed Heidler current given by Formula (15) (see the corresponding analysis done in Section 3.1).

Lightning Channel Distance Influence
Having initially computed electric field waveforms (see Figure 2) and transformation functions (see Figure 3), the lightning channel base current for different distances from the strike point r = 50 m, 5 km and 100 km was computed using Equation (14). The distance distribution was done according to [17,29]. For range r = 50 m the electrostatic component dominates the total electric field signature, for r = 5 km all three components are comparable and for r = 100 km the radiation component dominates the lightning electric field. Therefore, the decision to use such distance range will cover wide change of the lightning electric field peak value and its corresponding waveform shape variation which is representative for selected distances. The observation point was located at the ground surface (z = 0 km). In this paper, the propagation over perfect ground was assumed. Therefore, the electric field wave was not attenuated during its propagation. This is a reasonable approximation for lightning events observed within distance of several kilometers from the lightning channel. However, the effect of the electric field attenuation could be included using the Rubinstein formula [37] through the correction of the observed lightning electric field. The calculated currents (see Figure 4) were in good agreement with the reference model given in Formula (15). The relative uncertainty of the computation was estimated from: where iBc is the calculated lightning channel base current, iBs is the assumed Heidler-type reference current, N is the number of probes. N = 10000 for the 100 μs length of current vectors and sampling time ΔT = 10 ns. The average uncertainty of entire data set for 50 m, 5 km and 100 km was ΔiB = 3.43 × 10 −4 . Therefore, the lightning channel base current waveform was not influenced by the distance parameter. Some variation of the uncertainty was observed close to the lightning channel for ΔT > 100 ns. The slope was disturbed by the low quality of the electric field onset. In practice, the numerical simulation showed that the sampling time should be at least ten times lower than the initial slope time of the electric field. The calculated currents (see Figure 4) were in good agreement with the reference model given in Formula (15). The relative uncertainty of the computation was estimated from: where i Bc is the calculated lightning channel base current, i Bs is the assumed Heidler-type reference current, N is the number of probes. N = 10,000 for the 100 µs length of current vectors and sampling time ∆T = 10 ns. The average uncertainty of entire data set for 50 m, 5 km and 100 km was ∆i B = 3.43 × 10 −4 . Therefore, the lightning channel base current waveform was not influenced by the distance parameter. Some variation of the uncertainty was observed close to the lightning channel for ∆T 100 ns. The slope was disturbed by the low quality of the electric field onset. In practice, the numerical simulation showed that the sampling time should be at least ten times lower than the initial slope time of the electric field.

Sensitivity to Electric Field Noise
In practice, during measurements of lightning signatures, there is always some portion of additional noise superimposed on electric field waveforms. Therefore, the influence of the initial electric field noise on resulted noise to signal (N/S) level for the computed lighting channel base current was verified. The electric field noise was added in Matlab using the randn() function. The normal distribution of input noise with 95% confidence level was assumed.
An example result for the 2 V/m noise level superimposed on electric field at r = 5 km for ΔT = 10 ns is given in Figure 5. This particular simulation was selected being as the most representative for entire data set. In this specific case, the relative noise levels for electric field and base current were ΔEz = 33.3 × 10 −3 and ΔiB = 25.8 × 10 −3 , respectively. The assumed noise level of lighting electric field was based on typical values observed during registration of natural lighting electric fields. The corresponding 2.5% accuracy of the lightning channel base current estimation is good enough for most practical applications.

Sensitivity to Electric Field Noise
In practice, during measurements of lightning signatures, there is always some portion of additional noise superimposed on electric field waveforms. Therefore, the influence of the initial electric field noise on resulted noise to signal (N/S) level for the computed lighting channel base current was verified. The electric field noise was added in Matlab using the randn() function. The normal distribution of input noise with 95% confidence level was assumed.
An example result for the 2 V/m noise level superimposed on electric field at r = 5 km for ∆T = 10 ns is given in Figure 5. This particular simulation was selected being as the most representative for entire data set. In this specific case, the relative noise levels for electric field and base current were ∆E z = 33.3 × 10 −3 and ∆i B = 25.8 × 10 −3 , respectively. The assumed noise level of lighting electric field was based on typical values observed during registration of natural lighting electric fields. The corresponding 2.5% accuracy of the lightning channel base current estimation is good enough for most practical applications.

Sensitivity to Electric Field Noise
In practice, during measurements of lightning signatures, there is always some portion of additional noise superimposed on electric field waveforms. Therefore, the influence of the initial electric field noise on resulted noise to signal (N/S) level for the computed lighting channel base current was verified. The electric field noise was added in Matlab using the randn() function. The normal distribution of input noise with 95% confidence level was assumed.
An example result for the 2 V/m noise level superimposed on electric field at r = 5 km for ΔT = 10 ns is given in Figure 5. This particular simulation was selected being as the most representative for entire data set. In this specific case, the relative noise levels for electric field and base current were ΔEz = 33.3 × 10 −3 and ΔiB = 25.8 × 10 −3 , respectively. The assumed noise level of lighting electric field was based on typical values observed during registration of natural lighting electric fields. The corresponding 2.5% accuracy of the lightning channel base current estimation is good enough for most practical applications.  More detailed and completed statistics on the lightning channel base current distortion are presented in Table 1. The table shows noise levels of the lightning channel base current i B at different distances r to the lightning channel. As expected, the noise level increased with distance to the strike point. For distances exceeding several tens of kilometers, the distortion was too high to estimate the current reasonably. An accuracy of approximately 1% was obtained within 5 km to the lightning channel for the initial electric field noise not greater than 1 V/m. The obtained good performance of the algorithm in the presence of electric field noise can be explained using Equation (10) where integration was applied to the vertical component of lighting electric field. The integration acts as a lowpass filter. Assuming the high frequency character of electric field noise the lighting current noise reduction tendency can be expected. The base current and electric field uncertainties were linearly dependent each other until ∆i B 3.43 × 10 −4 . This limitation was noticed for close distances to the lightning channel (see Table 1). It should be pointed out that from the analysis of Equation (14) the linear dependence of the lightning current noise level to the electric field noise can be expected. This feature is observed for most of data given in Table 1, especially for r = 100 km. Further improvement of accuracy could be obtained by applying an auxiliary lowpass filtering on the resulted lighting channel base current waveform, e.g., with the smooth() Matlab function. Lightning channel base current noise was estimated for ∆T = 10 ns. a Electric field noise at sampling time ∆T = 100 ns, b Electric field noise at sampling time ∆T = 1 ns.
The sampling time influence was verified in two last rows of Table 1. The variation of ∆i B was the highest when the lightning channel base current was calculated from the close electric field. In this case, the best precision of ∆i B = 1.25 × 10 −4 was obtained for ∆T = 1 ns. Decreasing the sampling time is another option of algorithm accuracy improvement, but only for electric field measurements done close to the lightning channel.
The obtained results were compared to other studies done in the literature. It was difficult to find a reasonable comparison reference because most research aimed to calculate the lighting current peak value only. A similar procedure of lighting waveform current computation was found in [32]. Yang et al. presented an algorithm of lighting current calculation form the corresponding magnetic field. They used the same Heidler type of current at the lightning channel base. For distance r = 30 km and 20 kA current peak value they reported ∆i B = 1.2 × 10 −3 which was close to the noise levels presented in this paper (Table 1). It should be pointed that the magnetic field noise level was not directly reported in [32]. From graphical analysis of the figures in [32] the noise level related to the lightning field close to this is seen in Figure 5a of this paper.

Restrictions and Limitations of Algorithm
In this subsection, the conditions which should be met to use the lightning channel base current computation algorithm in practical applications were analyzed. Some of the limitations were depicted in Sections 3.1 and 3.2. The accuracy of the algorithm, estimated for a given level of electric field noise, was better when close to the lightning channel. However, its sensitivity on the sampling time variation increased as well. Therefore, the optimal condition is an intermediate distance from lightning strike point where contributions of the electrostatic, induction, and radiation components are comparable. Moreover, for the middle distance, the transformation function is well conditioned because w e (r,z,t), w i (r,z,t), w r (r,z,t) are of the same range. The better performance of the developed method for close and middle lightning strokes can be an advantageous feature for its implementation in lightning location systems. Most of the LLS are normally configured to record far fields of the lightning discharges because the lightning current peak value is estimated from the radiation component of electromagnetic field. The registering stations are usually idle during close and middle distance stroke events and a lot of their computational resources are available [2]. Therefore, the lightning current waveform computations according to presented method (14) can be done giving additional data which can be further processed and saved to enrich LLS databases. One of the important parameters of lightning current which is currently not reported in LLS databases and can be estimated using the lightning current shape is the lightning current steepness. This parameter allows for estimating the induced effects of lightning and therefore can be used for, i.e., overvoltage protection purposes.
Another noticed limitation was the duration of lightning channel base current computation. The simulation was done on the one core of 2.9 GHz i7-2675QM processor. The sampling time varying from 1 ns up to 100 ns (see Figure 6). The remaining simulation parameters were the lightning channel base current duration of 100 µs, r = 5 km, z = 0 km, ∆z' = 2 m. The simulation time consisted of two major parts. The first part was the time needed for computation of the transformation function w(r,z,t). This time was about 0.5 s, and it was relatively static. The second part was the time of deconvolution procedure. It was strongly dependent on the number of samples of the output base current vector. The author estimated the computational complexity O(N) = N(N + 1)/2, where N is the number of samples. The second order exponential character of timing performance function was obtained in Figure 6 for ∆T 10 ns. For practical applications in real-time lightning location systems, the total duration of entire computation should be lower than the average lightning interstroke interval of about 60 ms [3,10]. This limit can be obtained by a reduction of the w(t) function estimation time, e.g., using the lookup-table method [38].

Restrictions and Limitations of Algorithm
In this subsection, the conditions which should be met to use the lightning channel base current computation algorithm in practical applications were analyzed. Some of the limitations were depicted in Sections 3.1 and 3.2. The accuracy of the algorithm, estimated for a given level of electric field noise, was better when close to the lightning channel. However, its sensitivity on the sampling time variation increased as well. Therefore, the optimal condition is an intermediate distance from lightning strike point where contributions of the electrostatic, induction, and radiation components are comparable. Moreover, for the middle distance, the transformation function is well conditioned because we(r,z,t), wi(r,z,t), wr(r,z,t) are of the same range. The better performance of the developed method for close and middle lightning strokes can be an advantageous feature for its implementation in lightning location systems. Most of the LLS are normally configured to record far fields of the lightning discharges because the lightning current peak value is estimated from the radiation component of electromagnetic field. The registering stations are usually idle during close and middle distance stroke events and a lot of their computational resources are available [2]. Therefore, the lightning current waveform computations according to presented method (14) can be done giving additional data which can be further processed and saved to enrich LLS databases. One of the important parameters of lightning current which is currently not reported in LLS databases and can be estimated using the lightning current shape is the lightning current steepness. This parameter allows for estimating the induced effects of lightning and therefore can be used for, i.e., overvoltage protection purposes.
Another noticed limitation was the duration of lightning channel base current computation. The simulation was done on the one core of 2.9 GHz i7-2675QM processor. The sampling time varying from 1 ns up to 100 ns (see Figure 6). The remaining simulation parameters were the lightning channel base current duration of 100 μs, r = 5 km, z = 0 km, Δz' = 2 m. The simulation time consisted of two major parts. The first part was the time needed for computation of the transformation function w(r,z,t). This time was about 0.5 s, and it was relatively static. The second part was the time of deconvolution procedure. It was strongly dependent on the number of samples of the output base current vector. The author estimated the computational complexity O(N) = N(N + 1)/2, where N is the number of samples. The second order exponential character of timing performance function was obtained in Figure 6 for ΔT < 10 ns. For practical applications in real-time lightning location systems, the total duration of entire computation should be lower than the average lightning interstroke interval of about 60 ms [3,10]. This limit can be obtained by a reduction of the w(t) function estimation time, e.g., using the lookup-table method [38].

Conclusions
In this paper, a method for the calculation of a current at the bottom of the lightning channel was presented. The lightning channel base current was computed from the cor-

Conclusions
In this paper, a method for the calculation of a current at the bottom of the lightning channel was presented. The lightning channel base current was computed from the corresponding vertical component of a lightning electric field and a transformation function. The function included information about the geometry and physical properties of the lightning electromagnetic impulse propagation system. The lightning channel base current was computed for different distances from the lightning channel. The algorithm was verified for a Heidler current source and a MTLL model. However, an application of other base current models and attenuation functions is still possible. The results showed an acceptable performance, especially within several kilometers to the strike point. The best operation was obtained for an intermediate distance where the electrostatic, induction, and radiation components of electric field were of the same range. The sensitivity of the lightning channel base current to the electric field noise was verified. The computational complexity was estimated on the basis of deconvolution method. Further reduction of the simulation time is still possible, particularly regarding the preprocessing of the transfer function. For applications where the simulation time is not a critical parameter, the presented algorithm can be applied in its generic form.