An Approach for Estimating Lightning Current Parameters Using the Empirical Mode Decomposition Method

Lightning parameters are needed in different engineering applications. For the prediction of the severity of transient voltages in power systems, an accurate knowledge of the parameters of lightning currents is essential. All relevant standards and technical brochures recommend that lightning characteristics should be classified according to geographical regions instead of assuming that these characteristics are globally uniform. Many engineers and scientists suggest that better methods for lightning current measurements and analyses need to be developed. A system for direct lightning current measurements installed on Mount Lovćen is described in this paper. Observed data were analyzed, and statistical data on parameters that are of interest for engineering applications were obtained, as well as correlations between various lightning parameters. Furthermore, a novel approach for classifying and analyzing lightning data from direct measurements based on empirical mode decomposition (EMD) is proposed. Matlab was used as a tool for signal processing and statistical analysis. The methodology implemented in this work opens possibilities for automated analysis of large data sets and expressing lightning parameters in probabilistic terms from the data measured on site.


Introduction
With the continuous increase in the complexity of the electrical power system, growing exposure to the adverse effects of different environmental factors is further aggravating problems associated with reliability and safe operation of power system. Extreme weather phenomena and lightning in particular is one of the most common causes of faults and power supply interruptions. Systematic reviews of numerous observations of lightning events all around the world and information thus obtained provide valuable tools to reduce vulnerability and improve overall performance of power system. Damage to the power system is caused by both direct and indirect lightning strokes. In order to assess lightning effects and to design effective protection systems accurate lightning current parameters must be used. Lightning current parameters are of great importance in insulation coordination procedure, e.g., if this procedure is not properly determined the energy of lightning discharge can exceed energy handling capability of power system components [1][2][3].
Three approaches can be used to obtain lightning data: direct measurements using instrumented towers, direct measurements using the technique of the artificial initiation of lightning and lightning location systems [4,5]. Formatted lightning data from modern lightning location systems include: time and date of lightning stroke, GPS coordinates (2D), lightning current amplitude, lightning type, height (for inter-cloud lightning) and 2D statistical error [6,7]. Lightning location systems do not provide data about front and tail time of lightning current waveform. It is well known that front time has a high influence on insulation in power systems, while, for example, energy stresses of surge arresters strongly depend on tail time of the overvoltage wave [8][9][10][11]. Therefore, when conducting simulations of power system transients it is necessary to predict these data, as it is recommended in the CIGRE Brochure 549 (2013) and last IEEE review of parameters of lightning strokes (2005) [12,13].
In the field of lightning research for engineering applications, the most important data are obtained by the analysis of directly measured lightning current waveform. When designing such measuring systems, the first important step is the selection of a location with high lightning activity. Modern lightning location systems should be used for selection of regions with high lightning activity. Furthermore, it is important to correctly select the components of the measuring system, which will be in accordance with the specific characteristics of the lightning current that should be measured. In addition, it is necessary to develop tools for adequate and precise processing and analysis of measured data.
If a record of lightning current waveform is available, it is then possible, using appropriate numerical technique, to determine different parameters associated with that specific waveform. In every measurement, however, the measuring device is affected by environmental disturbances, referred to as noise, that alter characteristics of the output signal. Presence of noise can cause serious errors in measurement signal processing. In the case of lightning current measurement, in general, parameter determination may be difficult task due to the fact that all measured lightning current data are contaminated by considerable levels of noise, so additional processing steps must be undertaken in order to minimise effects of noise. The classification of recorded lightning current waveforms based on polarity and multiplicity is another important consideration in lightning studies. When dealing with a large amount of data from lightning monitoring systems, it is impractical to classify and analyze data manually. For such studies it becomes necessary to develop methodologies for automated classification and extraction of waveform parameters from the recorded data.
In the recent literature on lightning research, processing techniques are insufficiently considered. Several established methods are frequently reported: the Fourier transform, the short Fourier transform (STFT) [14], and the wavelet transform (WT) [15,16]. In [17], a time domain digital processing system for lightning current waveform parameters extraction is described. Using this approach procedure for parameters extraction from negative lightning flash with only one stroke was developed. In [18][19][20], the authors use Empirical Mode Decomposition (EMD) for discharge electric field pulse analyses, but in the recent literature, this method was not used for lightning current waveshape analyses. Therefore, this paper proposes a novel approach for analyzing lightning current waveform parameters and it is based on EMD. Continuing the work in [17], the authors expanded capabilities of previously developed signal processing procedure in terms of introducing new algorithm, increasing number of analysed features and including all typical types of discharges.
In this study, data from direct lightning current measurement system is analyzed using novel signal processing and parameter estimation technique and detailed statistics for one year observation period is presented. Correlations between various lightning parameters are established. This paper includes following contributions: • a novel approach to lightning current waveform processing based on EMD for more accurate automatic lightning classification and lightning parameter extraction is introduced; • statistical properties of lightning current parameters that are of great importance for engineering applications in region of mountain Lovćen (Montenegro) is presented; • empirical expressions for cumulative peak current distributions of first and subsequent strokes are determined.
The rest of the paper is organized as follows. Section 2 introduces the types of lightning discharges and lightning current parameters. Description of observation site and lightning monitoring system is given in Section 3. Section 4 describes novel approach for lightning data analyses. Statistics of lightning current parameters and correlations between parameters are reported in Section 5 and the paper is concluded in Section 6.

Lightning Discharges and Lightning Parameters for Engineering Application
According to [4,21], cloud-to-ground lightning discharges are divided into four main types: upward or downward (by the direction of the motion of the initial leader) and positive or negative (by the sign of the charge deposited along the channel). Classification in [4] includes only "unipolar flashes" that transport charge of one polarity to ground. Lightning flashes that transport both negative and positive charges to ground called "bipolar flashes" are not included in this classification. More than one lightning stroke can hit the same place on the ground in short time interval. To identify number of strokes in a single flash the term multiplicity is introduced. Usually first strokes have larger currents than subsequent strokes that occur both in new and in previously formed channel. Further details on lightning phenomenon can be found in [4].
Different types of lightning discharges have different impacts on power systems. Therefore, it is very important to identify the parameters of lightning current. According to CIGRE publication [22], typical lightning current waveshape, shown in Figure 1    Based on the above parameters, the time duration of current front, t f , is defined as time interval from t 0 to t p and is determined as shown in Figure 1. The time to half value, t h , represents time interval from t 0 to the 50% value point of the first peak (t 50 ). The energy in a lightning flash is assessed generally by its charge, Q, defined as: and specific energy E is defined as:

Location
Mountain Lovćen, with peak altitude of 1749 m is located in southwestern Montenegro, near the Adriatic Sea. Geographical location of mountain Lovćen and tower on which measuring equipment is installed are shown in Figure 2. Lightning current measurement equipment was installed on the 88 m high broadcasting tower, one of the most important communications hub in the region. The decision to install measurement equipment on this site was made based on previous reports and on data available from lightning location system (LLS).
LLS data have revealed that Lovćen, with 1063 strokes per square kilometer per year, has more than 100 times above median value in this region. Another contributing factor when choosing this site was the 500 kA maximum lightning current amplitude recorded by LLS reported in [23].

Measuring System
The lightning measurement system was constructed from a sensor, recording unit, power supply unit, central processing unit and user interface. A installed hardware is presented in Figure 3, while detailed block diagram of the system is shown in Figure 4.
The sensor unit containing current transformer, electric field sensor and IP camera, was installed on tower top.
The lightning current sensor used was current transformer with 500 kA input range. Changes in electric field are registered using electric field sensor BOLTEK-EFM 100 Atmospheric Electric Field Monitor. IP camera (UFG1122 HD IP Camera) with 120 fps (frames per second), equipped with SD card and infrared cut filter for day/night operations. The ecording unit is based on an industrial computer that records data from sensor unit. The acquisition unit is a four-channel card with an acquisition sampling rate of 8 MSa/s per channel and 15 bits vertical resolution. Accurate timing is provided by an integrated GPS receiver. Local ethernet connection is used for communication with the remote server. Recording and processing units were installed inside the broadcasting tower and are supplied from AC mains. A low loss cable (RG218) was used to connect output from the current transformer to the input of acquisition unit. A voltage attenuator with the ratio 10/1 was installed at the acquisition card input. Data are transferred in real-time via internet to the central server and stored into the integral information system. Detailed information on system architecture is provided in [24].

Signal Processing and Parameter Estimation
Signals obtained directly from the measuring system contain considerable amount of noise. Important lightning current parameters can be distinguished without filtering directly from the measured current shape, but it is not possible to determine the exact values. Therefore, in order to extract the values of important parameters, it is necessary to apply the appropriate signal processing technique. To improve the parameter extraction process, empirical mode decomposition was introduced for lightning current waveform denoising.

EMD Algorithm and Parameters Determination
Empirical mode decomposition represents a method of breaking down a signal without leaving the time domain. This method is a powerful tool for analyzing natural signals, which are mostly non-linear and non-stationary. It serves as a signal processing technique based on an empirical and algorithm defined method. EMD can adaptively decompose a complex signal into a set of complete, almost orthogonal components. These components are known as Intrinsic Mode Functions (IMFs).
EMD filters out IMFs without requiring any preliminary understanding of the nature and quantity of the IMF components in the data. The main advantage of EMD compared with the widely used wavelet-based technique is that EMD can be used to decompose a signal without specifying the basics functions in advance, and the degree of decomposition is adaptively determined in accordance with the nature of the signal to be decomposed [20].
Due to its performance, EMD has been widely used in many disciplines. EMD was first proposed by Huang et al. in [25] and this approach is used in computational neuroscience, biomedical signal processing, climate signal analysis, audio signal processing, image processing, and seismic signal and discharge electric field pulse analyses [26]. Therefore, details of the EMD algorithm and denoising principles can be found in the literature [25,[27][28][29].
This paper introduces the EMD algorithm into analyses of lighting current waveforms for parameter determination. This study presents the concept of EMD and its application to lightning current signal processing. Figure 5 shows the proposed EMD-based adaptive thresholding lightning current enhancement concept.
The basic steps of proposed method are as follows: • Step 1 Applying EMD algorithm to the raw data (noisy lightning current waveshape), which decomposes input signal in to IMFs. • Step 2 IMFs segmentation into frames. • Step 3 Frame classification into noise dominant and signal dominant frames. • Step 4 Adaptive thresholding. • Step 5 Combining of denoised IMFs. • Step 6 Parameter determination from enhanced signal.
Proposed novel lightning current signal processing and parameter determination method was implemented in MATLAB.

Evaluation of Used Methodology
In order to evaluate accuracy of proposed lightning current parameters estimation procedure several experiments were conducted. The proposed processing method was applied to a set of three types of synthetic signals. Three types of standard CIGRE concave lightning current waveforms, with parameters given in Table 1, are generated with sampling rate of 8 Msamples/s. It was assumed that the measured lightning current can be represented at most in accordance to CIGRE concave lightning current model. In addition, it is assumed that the noise observed in measured signals is additive white Gaussian noise. These assumptions are reasonable due to fact that most of recorded lightning strokes are very similar to the assumed CIGRE model [12,13].  Large number of randomly generated synthetic signals were generated in Monte Carlo simulations that were used for evaluation of performance of proposed method. For each type of standard lightning current waveshapes from Table 1, 1000 synthetic signals with a signal-to-noise ratio (SNR) in the range of 0 to 25 dB were generated. These signals were then processed, using the proposed method described in Section 4.1; enhanced signals were obtained and subjected to classification and parameter estimation algorithms.
The difference between estimated parameters obtained from enhanced signals relative to original signal parameters (in Table 2) are listed in Table 2 and were used as criterion for performance evaluation. From Table 2, it is clear that accuracy of some estimated parameters (peak current, tail time, total charge and specific energy) is almost independent of noise, while parameters such as front time and steepens are very sensitive to noise level. Estimated peak current values are almost constant in entire range of simulated SNR.
For peak currents greater than 10 kA, the average relative error was ±2.44% (from 0.81% to 6.65% with relative standard deviation below 7%). For lower peak currents (lower than 10 kA) average estimation error is slightly higher (±8.53%, with greatest error at 0 dB with value of 25.63%), and for SNR greater than 5 dB average estimation error was ±5.11%. These results suggests that proposed procedure can estimate with high accuracy peak current values in wide range of SNR (from 5 dB to 25 dB). As expected, for low current amplitudes and for SNR below 5 dB accuracy is decreasing.
As it can be seen from Equations (1) and (2), total charge and specific energy are functions of lightning current and due to this fact these parameters are also estimated with high accuracy within entire region of simulated SNR with average relative error of ±2.84% and ±0.77%, respectively. Noise level does not significantly affect these parameters since integration, in principle, represents a low pass filter. Time duration and steepness parameters estimation, however, are in general more variable and more sensitive to SNR. Tail time duration is estimated with the average relative error of ±3.55%. Waveform parameters front time and steepness in the investigated range of values and noise levels are estimated with higher average relative errors of 34.99% and 16.71%, respectively. The estimation of these parameters is significantly affected by the noise level. For fast rising currents (t f around 1 µs) in extreme case (at SNR = 0 dB), estimation errors may be up to 200%, and in this case, estimation is not reliable. However, in the more common range of SNR values (>10 dB), the average relative error for the front time is ±15.04%, while for the steepness, it is ±7.32%. Considering the standard tolerances given in [30,31] (for front time ±30% and for tail time ±20%), the obtained results are within the acceptable range.

Results of Observation
During the observation period of one year, 163 lightning events were recorded. Using the developed approach for automated data analysis, different types of lightning discharges were identified. The total number of lightning flashes was 64. The analyzed events were classified as given in Table 3. Detailed statistics were performed only for negative strokes due to the representative sample size.

Statistical Distribution of Lightning Parameters
It is generally agreed that the statistical distribution of lightning parameters follows the log-normal distribution, where the statistical variation of the logarithm of a random variable, x, follows the Gaussian distribution. The log-normal probability density function, p(x), is defined as in [13]: where σ lnx is standard deviation of lnx, and x m is median value of x. Therefore, x m and σ lnx need to be known to estimate the statistical distribution of a lightning parameter. The cumulative probability, P c (x), that the parameter will exceed x, is given by integrating Equation (3) between u 0 and ∞, resulting in: For approximating the log-normal distribution P c of lightning current peak, a simplified equation given by Anderson in [22,32] is also used: where µ and ρ are calculated from empirical data. Various correlations among lightning parameters have been found [13]. Assuming log-normal distributed random variables x and y, relationship between x and y can be expressed as:

Negative Flashes
As can be noted from Table 3, 59 first negative strokes and 86 subsequent negative strokes were analyzed. For the purpose of such analysis novel proposed processing method was used. The statistical distribution of multicomponent lightning flashes recorded in this study and compared with Anderson and Ericson (in [22]) is given in Figure 6. The frequency of the occurrence of multicomponent flashes in this region is very similar to that of Anderson and Ericson which is widely accepted.
Classification, analysis and parameter determination are more challenging tasks for lightning flashes that consists of more stokes than for lightning flashes with single stroke. Therefore, as an example, a multicomponent lightning flash that consists of four negative strokes is presented in Figure 7. Figure 7 shows the originally measured signal and enhanced signal. Determined parameters that are important for engineering applications are presented in Table 4.
Statistical parameters resulting from the measurements from this study during observation period of one year are given in Tables 5 and 6. First and subsequent negative currents were considered. After the log transformation, the Lilliefors test for normality, considering the 95% level of significance, was applied for the complete data set. It proved to be significant for most parameters of first negative strokes, similarly to the results presented by Anderson and Ericson in [22,33].     Cumulative statistical distributions of various parameters for the first and subsequent strokes are presented in the figures below (from Figure 8), as well as probability plots for lognormal distribution. From the figures, it can be seen that the measured data for the first and subsequent negative stokes are in good agreement with the theoretical cumulative distribution function. As indicated by the p-value from Table 5 at a significance level of 95%, it can be concluded that most of the analyzed parameters of the first negative strokes are distributed according to log-normal distribution. The total charge and maximum steepness for first negative strokes, as well as most parameters of the subsequent strokes, are similarly distributed, but failed the Lillifors test and, therefore, the log-normal distribution of these parameters cannot be confirmed. Very important formulas for cumulative probability as a function of the peak current for the first and subsequent stokes are performed from these results. The approximate expressions for cumulative probability of first and subsequent negative strokes current are given by Equations (7) and (8) (8) It is well known that these expressions have direct application in assessment of the lightning performance of electrical systems especially in insulation coordination studies. Therefore, it is of great importance to develop such formulas for different regions worldwide.

Correlations between Parameters for First Stoke in Negative Flashes
Correlations between various parameters of recorded lightning current waveshapes are considered by using fitting curves given by Equation (6) or it can be also represented with equation: ln(y) = ln(a) + dln(x) Correlation coefficients among parameters are given in Table 7. Results indicated that correlations are observed between almost all parameters. Few correlations that are not significant are marked in tables. Based on p-values from Table 8 statistically significance of correlations are confirmed. Strongest correlations were observed between peak current and all other parameters except tail time. Table 9 presents correlation expressions (described with function given in (9)) along with correlation coefficients r for such parameters, following logarithmic linear regression.   Similarly, as it is published in the literature, in this study, a correlation between peaks current and front time was observed. Additionally, a similar correlation exists between the peak current and maximum steepness (see Figure 9), with a similar correlation coefficients in available literature. As expected, a very strong correlation was observed between the peak current and total charge and specific energy. An interesting observation is that both the front and tail time are negatively correlated with steepness (see Figure 10).

Positive Flashes
Four positive lightning flashes were observed, each with single stroke. The parameters for positive flashes are given in Table 10. As an example, the original and enhanced positive stroke are shown in Figure 11. Original signal Enhanced signal Figure 11. Positive lighting flash original and enhanced signal.

Bipolar Flashes
During the observation period of one year, one bipolar lightning flash was recorded on the 26th of February 2016, and it is presented in Figure 12. The parameters for the bipolar flash are given in Table 11. All parameters were determined for the positive and negative parts of the lightning current waveshape.

Conclusions
This paper presents the statistics of lightning current parameters obtained by processing the data collected by direct measurement at the broadcasting tower on Mount Lovćen. The analyzed data were collected during the one-year observation period. A new EMD-based approach for more accurate lightning classification and lightning parameter extractionwas applied to the data analysis. The introduction of the EMD algorithm significantly improved the accuracy of determining lightning current parameters compared to the methods previously used by the authors. Unlike conventional filters, using this algorithm in the proposed scheme, the phase shift of the signal is almost eliminated.
Based on the cumulative distributions of the peak current of the first and subsequent strokes, the formulas for determining the probability of the occurrence of the peak current and the expression relating the peak current and other important parameters were generated. These expressions can be used directly in power system analyses. The statistical data for this region showed that most of the parameters for the first negative stroke are distributed according to the lognormal distribution and are very similar to their representation in contemporary literature. This study also confirmed that most of recorded events are negative lightning strokes, while positive and bipolar lightning strokes are rare. However, positive and bipolar lightning flashes are very dangerous, especially bipolar flashes (which transmit a large amount of energy to the ground), and should be taken into account when designing lightning protection.
Many uncertainties regarding lightning events presently exist, and therefore better methods for lightning current measurements and waveshape analyses should be developed. It should be continued with efforts to collect data for formulation of lightning parameters according to geographical regions and for developing important formulas for power system lightning protection as well as developing correlation expressions among lighting parameters.
For the signal processing methodologies used, it was shown that the accuracy of determining the front time and steepness should be improved. For this improvement, better acquisition units should be installed in the measurement system, for example, with a much higher sampling rate (50 MSa/s or 100 MSa/s) in order to be able to better record front time that has a very short duration (from 1 to several microseconds). In the future, the dynamic characteristics of the measurement system should be taken into account when analyzing the signal-to-noise ratio in order to further improve the signal processing.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.