Evaluation of lightning location using measured induced voltage obtained from distribution power networks

Abstract

Lightning location is a significant issue in the protection of transmission lines, renewable energy sources, and electrical equipment. In this article, a new technique for the determination of lightning striking points is been proposed. This method is depending on measured values of lightning-induced voltage obtained from distribution power lines in the vicinity of the lightning channel. The proposed method considers lightning location as well as lightning current wave shape at the channel base. It was validated using a number of measured return stroke currents, and the outcomes are debated accordingly. All electromagnetic field components were taken into consideration in the proposed method, and lightning current wave shape was determined in contrast with the other, widely used methods. Evaluated lightning current wave shapes, based on local measured data, can be used to examine and improve existing protection schemes in renewable energy generation farms, transmission lines, and other power sites. Moreover, lightning maps and ground flash density can be obtained from the evaluated lighting location data to be used in lightning risk assessment studies.

KEYWORDS:

• Lightning-induced voltage
• lightning location
• lightning
• electromagnetic fields
• lightning current

Abbreviations

BG	=	Bruce and Golde model
ECM	=	Engineering current models
LC	=	Lightning current
LEF	=	Lightning electric fields
LIV	=	Lightning-induced voltage
LL	=	Lightning location
MDF	=	Magnetic direction finder
MTLE	=	Modified transmission line with exponential decay model
MTLL	=	Modified transmission line with linear decay model
PSO	=	Particle swarm optimization algorithm
TCS	=	Travelling current source model
TDOA	=	The time difference of arrival
TL	=	Transmission line model
TOA	=	Time of arrival technique
TOGA	=	Time of group arrival
VLF	=	Very low frequency

Introduction

Protection of power units and power farms against lightning is an important issue. Site selection for power networks and renewable farms needs to consider lightning risk assessment based on lightning maps and ground flash density at the project site [1,2]. Therefore, evaluation of LL and LC wave shape, and development of related data banks can be helpful toward increasing the protection level of existing power sites and designing proper protection schemes for new power stations [3]. Several methods have been proposed and used to evaluate lightning location. These have normally been based on the radiation part of LEF [4–6]. However, they contain inherent errors due to their neglect of other field components and ground effects on propagation fields. The most widely used method can only evaluate current peak or current values within a certain range of frequency when the full shape of the current is a significant issue in both the design of surge protection devices and system performance [7,8]. One of the important challenges for electric field-based methods is measuring electric fields and extracting lightning electric fields from other field pollutants, as they have been used in algorithms(as input data).

In this study, a new method is proposed to evaluate LL as well as the full shape of LC based on measured values of LIV on the distribution lines in the vicinity of the lightning channel. In this method, the measured induced voltage data are used as input data as all field components were taken into account indirectly through induced voltage values. Therefore, lightning at different ranges of radial distance can be considered using the proposed method, in contrast with the widely used electric field-based methods. Likewise, the presented method can evaluate the full shape of LC to prepare the required data for protection and evaluation of system performance in power stations. The proposed method is a voltage-bases algorithm. Measurement of voltage on the distribution line has fewer operational obstacles compared to field-based methods. Moreover, the measuring points can be developed easily on different distribution lines. In this article, the performance of the proposed procedure was evaluated using different recorded return stroke current samples. Results will be analysed and considered accordingly. The fundamental presumptions in this paper are as follows:

1. Return stroke velocity as a constant value

2. Lightning as a vertical channel

3. Lightning branching and its effect are ignored

4. Ground conductivity is set to the perfect condition

Lightning location systems and inverse procedure algorithms

Several studies that have evaluated lightning location are summarized in Table 1.

Table 1. Overview of previous studies on LL systems.

Method	Description
TOA [9,10]	This technique is based on the arrival time of radiated fields as they have been measured at different stations - By using a number of sensors through spherical geometry - Sense electromagnetic field radiated by lightning flash in the frequency range 1 kHz to 1MHz
TDOA [4,11]	In this method, the LIV on a cable with two terminals at both ends have been considered and the time differences between travelling times of wave between main source and two terminals and also the main component and first reflected wave have been taken into account.
TOGA [12,13]	This method has developed TOA method in a range of VLF by considering on the phase slope with respect to frequency (6–22 kHz) that can be obtained from first microsecond of wave
MDF [12,14]	In this method by using two current loops perpendicular to each other at each measuring station, the induced current due to lightning at different directions have been evaluated and by comparison and cross-section between field vectors and considering their angles (at different stations), lightning location has be determined.

Moreover, different inverse procedure algorithms have been proposed in previous research that considered lightning current wave shape along the channel as listed in Table 2.

Table 2. Overview on the previous researches on the lightning inverse procedure algorithms.

Method	Description
Evaluation of current peak depending on radiation term of LEF [12]	In these methods, based of radiation part of vertical electric field, the current peak based of known current model has been evaluated as the radiation part is more effective only at far distances with respect to lightning. In these models, the TL or TCS current models (see Table 3) have been set as a basic assumption.
Evaluation of current at certain frequency range using electric field sensors with a certain distance with lightning channel [8,15]	In these methods, the field sensors should be installed in the range 5–10 km with respect to lightning channel while before lightning event, the striking point is not clear and setting of sensor location is a significant challenge. Moreover, the current has been evaluated in a certain frequency range (50Hz–1MHz) and the full shape of current cannot be evaluated. In these methods, lightning location and the radial distance have been set as known information.
Evaluation of current based on both electric and magnetic fields in a station [8,16]	In these methods, the lightning location and the radial distance have been set as known information and the lightning location has been set at closed distances (static component). Moreover, they evaluate the current value at certain frequency range.

Proposed method

In order to evaluate LL as well as LC wave shape at the channel base, lightning-induced voltage signals were used in the proposed method as input data. Figure 1 shows the geometry of the problem.

Figure 1. Geometry of the study problem.

Figure 2 illustrates a block diagram of the study, where in the induced voltage signals were obtained from overhead distribution lines using high voltage probes. They were used in the proposed algorithm that will be presented in the next part. Output of the proposed method was the LL as well as the LC wave shape. Noted that the LC at the channel base (striking point with ground), expressed using various current functions as the sum of two Heilder functions, has been used in this study as presented by Equation (1): (1) $I\left(0,t\right)=\frac{I_{01}}{{\eta }_{01}}\frac{{\left(\frac{t}{{\mathrm{\Gamma }}_1}\right)}^n}{1+{\left(\frac{t}{{\mathrm{\Gamma }}_1}\right)}^n}\exp \left(\frac{-t}{{\mathrm{\Gamma }}_2}\right)+\frac{I_{02}}{{\eta }_{02}}\frac{{\left(\frac{t}{{\mathrm{\Gamma }}_3}\right)}^n}{1+{\left(\frac{t}{{\mathrm{\Gamma }}_3}\right)}^n}\exp \left(\frac{-t}{{\mathrm{\Gamma }}_4}\right)$(1) where $I_{01},I_{02}$ are the peaks of the LC, ${\mathrm{\Gamma }}_1,{\mathrm{\Gamma }}_2$ are the front time parameters, ${\mathrm{\Gamma }}_3,{\mathrm{\Gamma }}_4$ are the decay-time parameters, n is the exponent (the rate is between 2 to10), $\begin{array}{rl}{\eta }_{01}& =\mathrm{e}\mathrm{x}\mathrm{p}\left[-\left({\mathrm{\Gamma }}_1/\phantom{{\mathrm{\Gamma }}_1{\mathrm{\Gamma }}_2}{\mathrm{\Gamma }}_2\right){\left(n\times \frac{{\mathrm{\Gamma }}_2}{{\mathrm{\Gamma }}_1}\right)}^{\frac{1}{n}}\right]\end{array}$ $\begin{array}{rl}{\eta }_{02}& =\mathrm{e}\mathrm{x}\mathrm{p}\left[-\left({\mathrm{\Gamma }}_3/\phantom{{\mathrm{\Gamma }}_3{\mathrm{\Gamma }}_4}{\mathrm{\Gamma }}_4\right){\left(n\times \frac{{\mathrm{\Gamma }}_4}{{\mathrm{\Gamma }}_3}\right)}^{\frac{1}{n}}\right]\end{array}$

Figure 2. Block diagram of the study.

Lightning current at various altitudes along its channel can be expressed using disparate models, e.g. dynamic current models [17,18] electromagnetic current models, distributed circuit models, ECM. In this article, the generic analytical form of the ECM was used as declared in Equation (2): (2) $I\left(z^{\mathrm{\prime }},t\right)=P\left(z^{\mathrm{\prime }}\right)\times I\left(0,t-\frac{z^{\mathrm{\prime }}}{v}\right)\times u\left(t-\frac{z^{\mathrm{\prime }}}{v_f}\right)$(2) where I(z’,t) is the current function along channel, I(0,t) is the current at channel base, P(z’) is a reduction factor (depending on charge height), z’ is the impermanent charge altitude along the channel, v is the current wave propagation velocity, $v_f$ is the skyward propagation front velocity, u is the Heaviside function defined as $u\left(t-\frac{z^{\mathrm{\prime }}}{v_f}\right)=\left\{\begin{array}{c}1\text{for}t\ge \frac{z^{\mathrm{\prime }}}{v_f}\\ 0\text{for}t<\frac{z^{\mathrm{\prime }}}{v_f}\end{array}\right.$

Table 3 widely used ECM based on Equation (3), where λ is a decay constant ($1\sim 2$ km), $H_c$ is the cloud altitude, c is$3\times {10}^8\frac{\text{m}}{\text{s}}$, and $v_f$ is a constant (between c/2 to 2c/3) [19].

Table 3. Parameters of widely used engineering current models [19].

Model	P(z’)	v
BG	1	$\mathrm{\infty }$
TCS	1	–c
TL	1	$v_f$
MTLL	$\left(1-{\displaystyle \frac{z^{\mathrm{\prime }}}{H_c}}\right)$	$v_f$
MTLE	$\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{z^{\mathrm{\prime }}}{\lambda }}\right)$	$v_f$

Lightning-induced voltage for each term of Equation (1), based on ECM, can be evaluated from this equation as the observation point is set on an overhead line (Equation 3). It should be mentioned that the sum of two Heidler functions has been used in this study. The induced voltage due to the first and the second terms in Equation (1) can be evaluated by using $(I_{01},{\mathrm{\Gamma }}_1,{\mathrm{\Gamma }}_2,{\eta }_{01})$and ($I_{02},{\mathrm{\Gamma }}_3,{\mathrm{\Gamma }}_4,{\eta }_{02}$) in Equation (3), respectively, and the total induced voltage is the sum of two obtained induced voltages from two terms of Equation (1) [20]. Noted that the expressed induced voltage expression is based on Agrawal coupling model. (3) $\begin{array}{rl}{V}_{ind}\left(\begin{array}{c}x\\ y=d\\ z=h\\ t=t_n\end{array}\right)& =\left[A_1^{\mathrm{\prime }}\left(\begin{array}{c}x\\ y=d\\ z=h\\ t=t_n\end{array}\right)+A_2^{\mathrm{\prime }}\left(\begin{array}{c}x\\ y=d\\ z=h\\ t=t_n\end{array}\right)+A_3^{\mathrm{\prime }}\left(\begin{array}{c}x\\ y=d\\ z=h\\ t=t_n\end{array}\right)\right]\\ & \times U\left(t_n-\frac{\sqrt{x^2+d^2+h^2}}{c}\right)\end{array}$(3) where $V_{ind}$ is the value of LIV $\begin{array}{rl}{A}_1^{\mathrm{\prime }}\left(\begin{array}{c}x\\ y\\ z\\ t\end{array}\right)& =\frac{-1}{2\mathrm{\Delta }h}\sum _{m^{\mathrm{\prime }}=1}^{k^{\mathrm{\prime }}+1}{b}_m\times \left[{\overrightarrow{E}}_z\left(\begin{array}{c}x\\ y\\ z=\left(m^{\mathrm{\prime }}-1\right)\times \mathrm{\Delta }h\\ t=t_{n-1}\end{array}\right)+\mathrm{\Delta }t\right.\\ & \times \sum _{i=1}^{n_t}\sum _{m=1}^{k+1}\left\{a_m{F}_2\left(\begin{array}{c}x\\ y\\ z=\left(m^{\mathrm{\prime }}-1\right)\times \mathrm{\Delta }h\\ t=t_n\\ z^{\mathrm{\prime }}=h_{m,i}\end{array}\right)-{a^{\mathrm{\prime }}}_m{F}_2\left(\begin{array}{c}x\\ y\\ z=\left(m^{\mathrm{\prime }}-1\right)\times \mathrm{\Delta }h\\ t=t_n\\ z^{\mathrm{\prime }}={h^{\mathrm{\prime }}}_{m,i}\end{array}\right)\right\}\\ & \times \left.U\left(t-\frac{\sqrt{r^2+z^2}}{c}\right)\phantom{\begin{array}{c}x\\ y\\ z=\left(m^{\mathrm{\prime }}-1\right)\times \mathrm{\Delta }h\\ t=t_{n-1}\end{array}}\right]\end{array}$ $\begin{array}{rl}{A}_2^{\mathrm{\prime }}\left(\begin{array}{c}x\\ y\\ z\\ t\end{array}\right)& =-\frac{1}{4\mathrm{\Delta }x}\times \sum _{m^{\mathrm{\prime }}=1}^{q+1}{c}_m\times \left[{\overrightarrow{E}}_x\left(\begin{array}{c}x=x_{m^{\mathrm{\prime }}}\\ y=d\\ z=h\\ t=t_{n-1}-\frac{x_{m^{\mathrm{\prime }}}-x}{c}\end{array}\right)+\mathrm{\Delta }t\right.\\ & \times \sum _{i=1}^{n_t}\sum _{m=1}^{k+1}\left\{a_m{F}_1\left(\begin{array}{c}x=x_{m^{\mathrm{\prime }}}\\ y=d\\ z=h\\ t=t_n-\frac{x_{m^{\mathrm{\prime }}}-x}{c}\\ z^{\mathrm{\prime }}=h_{m,i}\end{array}\right)-{a^{\mathrm{\prime }}}_m{F}_1\left(\begin{array}{c}x=x_{m^{\mathrm{\prime }}}\\ y=d\\ z=h\\ t=t_n-\frac{x_{m^{\mathrm{\prime }}}-x}{c}\\ z^{\mathrm{\prime }}={h^{\mathrm{\prime }}}_{m,i}\end{array}\right)\right\}\\ & \times \left.U\left(t_n-\frac{x_{m^{\mathrm{\prime }}}-x}{c}-\frac{\sqrt{{x_{m^{\mathrm{\prime }}}}^2+d^2+h^2}}{c}\right)\phantom{\begin{array}{c}x\\ y\\ z=\left(m^{\mathrm{\prime }}-1\right)\times \mathrm{\Delta }h\\ t=t_{n-1}\end{array}}\right]\end{array}$ $\begin{array}{rl}{A}_3^{\mathrm{\prime }}\left(\begin{array}{c}x\\ y\\ z\\ t\end{array}\right)& =\frac{1}{4\mathrm{\Delta }{x}^{\mathrm{\prime }}}\times \sum _{m^{\mathrm{\prime }}=1}^{q+1}{c}_m\times \left[{\overrightarrow{E}}_x\left(\begin{array}{c}x=x{}_{m^{\mathrm{\prime }}}^{\mathrm{\prime }}\\ y=d\\ z=h\\ t=t_{n-1}+\frac{x{}_{m^{\mathrm{\prime }}}^{\mathrm{\prime }}-x}{c}\end{array}\right)+\mathrm{\Delta }t\right.\\ & \times \sum _{i=1}^n\sum _{m=1}^{k+1}\left\{a_m{F}_1\left(\begin{array}{c}x=x{}_{m^{\mathrm{\prime }}}^{\mathrm{\prime }}\\ y=d\\ z=h\\ t=t_n+\frac{x{}_{m^{\mathrm{\prime }}}^{\mathrm{\prime }}-x}{c}\end{array}\right)-a_m^{\mathrm{\prime }}{F}_1\left(\begin{array}{c}x=x{}_{m^{\mathrm{\prime }}}^{\mathrm{\prime }}\\ y=d\\ z=h\\ t=t_n+\frac{x{}_{m^{\mathrm{\prime }}}^{\mathrm{\prime }}-x}{c}\end{array}\right)\right\}\\ & \times \left.U\left(t_n+\frac{x{}_{m^{\mathrm{\prime }}}^{\mathrm{\prime }}-x}{c}-\frac{\sqrt{x{{}_{m^{\mathrm{\prime }}}^{\mathrm{\prime }}}^2+d^2+h^2}}{c}\right)\phantom{\begin{array}{c}x\\ y\\ z=\left(m^{\mathrm{\prime }}-1\right)\times \mathrm{\Delta }h\\ t=t_{n-1}\end{array}}\right]\end{array}$ $\begin{array}{rl}\mathrm{\Delta }{x}^{\mathrm{\prime }}& =\frac{x+\frac{{\left(ct-x\right)}^2-h^2-d^2}{2\left(c{t}_n-x\right)}}{q},\end{array}$ $\begin{array}{rl}\mathrm{\Delta }x& =-\mathrm{\Delta }{x}^{\mathrm{\prime }},x{}_{m^{\mathrm{\prime }}}^{\mathrm{\prime }}=-\frac{{\left(c{t}_n-x\right)}^2-h^2-d^2}{2\left(c{t}_n-x\right)}+\left(m^{\mathrm{\prime }}-1\right)\times \mathrm{\Delta }{x}^{\mathrm{\prime }},\end{array}$ $q$ is division factor (≥2), $\begin{array}{rl}{b}_m& =c_m=\left\{\begin{array}{c}2\text{for}{m}^{\mathrm{\prime }}=1\text{and}{m}^{\mathrm{\prime }}=q+1\\ 1\text{for}\text{ others}\end{array}\right.\end{array}$ $\begin{array}{rl}{F}_1& =\frac{I_{01}P\left(z^{\mathrm{\prime }}\right)x\left(z-z^{\mathrm{\prime }}\right)A_2}{4\pi {ϵ}_0{\eta }_{01}\left[{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^n+1\right]}\\ & \times \left\{\frac{{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^n\left(3c^2{\mathrm{\Gamma }}_2^2-3c{\mathrm{\Gamma }}_{2}R+R^2\right)}{c^2{\mathrm{\Gamma }}_2^2{R}^5}+\frac{\begin{array}{c}3n([{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^n+1]{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^{n-1}\\ -{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^{2n-1})\end{array}}{c{\mathrm{\Gamma }}_1{R}^4\left[{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^n+1\right]}\right.\\ & +\frac{2n\left({\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^{2n-1}-\left[{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^n+1\right]{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^{n-1}\right)}{c^2{\mathrm{\Gamma }}_1{\mathrm{\Gamma }}_2{R}^3\left[{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^n+1\right]}\\ & +\left.\frac{\begin{array}{c}2n^2{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^{3n-2}+n\left(n-1\right){\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^{n-2}{\left[{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^n+1\right]}^2\\ +\left(-3n^2+n\right){\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^{2n-2}\left[{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^n+1\right]\end{array}}{c^2{\mathrm{\Gamma }}_1^2{R}^3{\left[{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^n+1\right]}^2}\phantom{\frac{{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^n\left(3c^2{\mathrm{\Gamma }}_2^2-3c{\mathrm{\Gamma }}_{2}R+R^2\right)}{c^2{\mathrm{\Gamma }}_2^2{R}^5}}\right\}\end{array}$ $\begin{array}{rl}{F}_2& =\frac{I_{01}P\left(z^{\mathrm{\prime }}\right)A_2}{4\pi {ϵ}_0{\eta }_{01}\left[{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^n+1\right]}\\ & \times \left\{\frac{2{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^n\left(c{\mathrm{\Gamma }}_1{\mathrm{\Gamma }}_2-{\mathrm{\Gamma }}_{1}R\right)+2n{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^{n-1}{\mathrm{\Gamma }}_{2}R}{c{\mathrm{\Gamma }}_1{\mathrm{\Gamma }}_2{R}^3}-\frac{n\left(n-1\right)\left(x^2+y^2\right){\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^{n-2}}{c^2{\mathrm{\Gamma }}_1^2{R}^3}\right.\\ & +\frac{\begin{array}{c}\left(x^2+y^2\right)\left[{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^n\left(-3c^2{\mathrm{\Gamma }}_1{\mathrm{\Gamma }}_2^2+3c{\mathrm{\Gamma }}_1{\mathrm{\Gamma }}_{2}R-{\mathrm{\Gamma }}_1{R}^2\right)\right.\\ +\left.n{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^{n-1}\left(-3c{\mathrm{\Gamma }}_2^{2}R+2{\mathrm{\Gamma }}_2{R}^2\right)\right]\end{array}}{c^2{t}_1{t}_2^2{R}^5}\\ & +\frac{\begin{array}{c}\left(x^2+y^2\right)\left[2n^2{\mathrm{\Gamma }}_{2}R{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^{2n-2}+3nc{\mathrm{\Gamma }}_1{\mathrm{\Gamma }}_2{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^{2n-1}\right.\\ +\left.n\left(n-1\right){\mathrm{\Gamma }}_{2}R{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^{2n-2}-2n{\mathrm{\Gamma }}_{1}R{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^{2n-1}\right]\end{array}}{c^2{\mathrm{\Gamma }}_1^2{\mathrm{\Gamma }}_2{R}^4\left[{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^n+1\right]}\\ & -\left.\frac{2nc{\mathrm{\Gamma }}_1{\mathrm{\Gamma }}_2{R}^2{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^{2n-1}}{c^2{\mathrm{\Gamma }}_1^2{\mathrm{\Gamma }}_2{R}^4\left[{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^n+1\right]}-\frac{2n^2\left(x^2+y^2\right){\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^{3n-2}}{c^2{\mathrm{\Gamma }}_1^2{R}^3{\left[{\left(\frac{A_1}{{\mathrm{\Gamma }}_1}\right)}^n+1\right]}^2}\right\}\end{array}$ $\begin{array}{rl}{\overrightarrow{E}}_x\left(\begin{array}{c}x\\ y\\ z\\ t_n\end{array}\right)& ={\overrightarrow{E}}_x\left(\begin{array}{c}x\\ y\\ z\\ t_{n-1}\end{array}\right)+\mathrm{\Delta }t\times \sum _{i=1}^{n_t}\sum _{m=1}^{k+1}\left\{a_m{F}_1\left(\begin{array}{c}x\\ y\\ z\\ t=t_n\\ z^{\mathrm{\prime }}=h_{m,i}\end{array}\right)-a_m^{\mathrm{\prime }}{F}_1\left(\begin{array}{c}x\\ y\\ z\\ t=t_n\\ z^{\mathrm{\prime }}={h^{\mathrm{\prime }}}_{m,i}\end{array}\right)\right\}\\ & \times U\left(t_n-\frac{\sqrt{r^2+z^2}}{c}\right)\end{array}$ $\begin{array}{rl}{\overrightarrow{E}}_z\left(\begin{array}{c}x\\ y\\ z\\ t_n\end{array}\right)& ={\overrightarrow{E}}_z\left(\begin{array}{c}x\\ y\\ z\\ t_{n-1}\end{array}\right)+\mathrm{\Delta }t\times \sum _{i=1}^{n_t}\sum _{m=1}^{k+1}\left\{a_m{F}_2\left(\begin{array}{c}x\\ y\\ z\\ t=t_n\\ z^{\mathrm{\prime }}=h_{m,i}\end{array}\right)-a_m^{\mathrm{\prime }}{F}_2\left(\begin{array}{c}x\\ y\\ z\\ t=t_n\\ z^{\mathrm{\prime }}={h^{\mathrm{\prime }}}_{m,i}\end{array}\right)\right\}\\ & \times U\left(t_n-\frac{\sqrt{r^2+z^2}}{c}\right)\end{array}$ $\begin{array}{rl}R& =\sqrt{x^2+y^2+{\left(z-z^{\mathrm{\prime }}\right)}^2}\end{array}$ $\begin{array}{rl}{A}_1& =t-\frac{R}{c}-\frac{|z^{\mathrm{\prime }}|}{v}\end{array}$ $\begin{array}{rl}{A}_2& =exp\left(\frac{\frac{R}{c}-t+\frac{|z^{\mathrm{\prime }}|}{v}}{{\mathrm{\Gamma }}_2}\right)\end{array}$ ${\overrightarrow{E}}_x$ is the LEF at the x direction,

${\overrightarrow{E}}_z$ is the LEF at the z direction,

$\mathrm{\Delta }t$ is the time step,

$n_t$ is the quantity of time steps, $\begin{array}{rl}{t}_n& =\left(n_t-1\right)\mathrm{\Delta }t{n}_t=1,2,\dots ,n_{max}\end{array}$ $\begin{array}{rl}{a}_m& =\left\{\begin{array}{c}\frac{\mathrm{\Delta }{h}_i}{2\times k}\text{for}m=1\text{and}m=k+1\\ \frac{\mathrm{\Delta }{h}_i}{k}\text{for}\text{others}\end{array}\right.\end{array}$ $\begin{array}{rl}a{}_m^{\mathrm{\prime }}& =\left\{\begin{array}{c}\frac{\mathrm{\Delta }h{}_i^{\mathrm{\prime }}}{2\times k}\text{for}m=1\text{and}m=k+1\\ \frac{\mathrm{\Delta }h{}_i^{\mathrm{\prime }}}{k}\text{for}\text{others}\end{array}\right.\end{array}$ k is division factor($\ge 2$), $\begin{array}{rl}\mathrm{\Delta }{h}_i& =\left\{\begin{array}{c}\sqrt{{\left(\frac{r}{\chi }\right)}^2+{\left(\text{β}c{t}_{i-1}-z\right)}^2}-\sqrt{{\left(\frac{r}{\chi }\right)}^2+{\left(\text{β}c{t}_i-z\right)}^2}+\text{β}{\chi }^2\left\{\right(c{t}_i-c{t}_{i-1}\left)\right\}\\ \left\{-\left(\text{β}z-c{t}_i\right)-\sqrt{{\left(\text{β}c{t}_i-z\right)}^2+{\left(\frac{r}{\chi }\right)}^2}\right\}\times \text{β}{\chi }^2\text{for}i=1\end{array}\right.\end{array}$ $\begin{array}{rl}\mathrm{\Delta }h{}_i^{\mathrm{\prime }}& =\left\{\begin{array}{c}-\sqrt{{\left(\frac{r}{\chi }\right)}^2+{\left(\text{β}c{t}_{i-1}+z\right)}^2}+\sqrt{{\left(\frac{r}{\chi }\right)}^2+{\left(\text{β}c{t}_i+z\right)}^2}+\text{β}{\chi }^2\left\{\right(c{t}_{i-1}-c{t}_i\left)\right\}\\ \left\{\sqrt{{\left(\text{β}c{t}_i+z\right)}^2+{\left(\frac{r}{\chi }\right)}^2}-\left(\text{β}z+c{t}_i\right)\right\}\times \text{β}{\chi }^2\text{for}i=1\end{array}\right.\end{array}$ $\begin{array}{rl}{h}_{m,i}& =\left\{\begin{array}{c}{h}_{m=k+1,i-1}+\frac{(m-1)\times \mathrm{\Delta }{h}_i}{k}\\ \frac{(m-1)\times \mathrm{\Delta }{h}_i}{k}fori=1\end{array}\right.\end{array}$ $\begin{array}{rl}h{}_{m,i}^{\mathrm{\prime }}& =\left\{\begin{array}{c}h{}_{m=k+1,i-1}^{\mathrm{\prime }}+\frac{(m-1)\times \mathrm{\Delta }h{}_i^{\mathrm{\prime }}}{k}\\ \frac{(m-1)\times \mathrm{\Delta }h{}_i^{\mathrm{\prime }}}{k}fori=1\end{array}\right.\end{array}$

In order to evaluate LL and LC wave shape using two measured signals of LIV in first and second distribution lines (based on Figure 1), the relation between measured signals at different time steps of each measuring point was equal to the corresponding analytical expression as mentioned in Equation (3). $V_m^1$ and $V_m^2$ were defined as the measured, induced voltages in the first and second distribution lines(based on Figure 1), respectively. The ratio between measured induced voltage values and the corresponding analytical expressions in the first and second distribution lines were expressed as in Equation (4). Noted that the channel base current function with unbeknown parameters were set as the sum of two Heidler functions, and the analytical expression of LIV consisted of both terms of equation (1). Moreover, the MTLE model, with unbeknown parameters, was set as the current model along the channel. Noted that the MTLE model had been frequently used in previous studies as it could prepare the simulated induced voltage with a proper range of accuracy compared to the corresponding measured one however the current model in this paper can be replaced with other engineering current models(as tabulated in Table 3) [6,21–22]. The coordinates of lightning location were set as unknown parameters, and the ground level effect was ignored in this paper. Equation (4) can be converted to a nonlinear equation system for both measuring stations as expressed in Equation (5): (4) $\begin{array}{r}\left\{\begin{array}{c}{V}_m^1\left(\begin{array}{c}{x}_1\\ y_1\\ z_1\\ t=t_1\end{array}\right)=\left[A_1^{\mathrm{\prime }}\left(\begin{array}{c}{x}_1-x_L\\ y_1-y_L\\ z_1\\ t=t_1\end{array}\right)+A_2^{\mathrm{\prime }}\left(\begin{array}{c}{x}_1-x_L\\ y_1-y_L\\ z_1\\ t=t_1\end{array}\right)+A_3^{\mathrm{\prime }}\left(\begin{array}{c}{x}_1-x_L\\ y_1-y_L\\ z_1\\ t=t_1\end{array}\right)\right]U\left(t_1-{\xi }_1\right)\\ V_m^1\left(\begin{array}{c}{x}_1\\ y_1\\ z_1\\ t=t_2\end{array}\right)=\left[A_1^{\mathrm{\prime }}\left(\begin{array}{c}{x}_1-x_L\\ y_1-y_L\\ z_1\\ t=t_2\end{array}\right)+A_2^{\mathrm{\prime }}\left(\begin{array}{c}{x}_1-x_L\\ y_1-y_L\\ z_1\\ t=t_1\end{array}\right)+A_3^{\mathrm{\prime }}\left(\begin{array}{c}{x}_1-x_L\\ y_1-y_L\\ z_1\\ t=t_1\end{array}\right)\right]U\left(t_2-{\xi }_1\right)\\ ⋮\\ V_m^1\left(\begin{array}{c}{x}_1\\ y_1\\ z_1\\ t=t_n\end{array}\right)=\left[A_1^{\mathrm{\prime }}\left(\begin{array}{c}{x}_1-x_L\\ y_1-y_L\\ z_1\\ t=t_n\end{array}\right)+A_2^{\mathrm{\prime }}\left(\begin{array}{c}{x}_1-x_L\\ y_1-y_L\\ z_1\\ t=t_n\end{array}\right)+A_3^{\mathrm{\prime }}\left(\begin{array}{c}{x}_1-x_L\\ y_1-y_L\\ z_1\\ t=t_n\end{array}\right)\right]U\left(t_n-{\xi }_1\right)\\ V_m^2\left(\begin{array}{c}{x}_2\\ y_2\\ z_2\\ t=t_1\end{array}\right)=\left[A_1^{\mathrm{\prime }}\left(\begin{array}{c}{x}_2-x_L\\ y_2-y_L\\ z_2\\ t=t_1\end{array}\right)+A_2^{\mathrm{\prime }}\left(\begin{array}{c}{x}_2-x_L\\ y_2-y_L\\ z_2\\ t=t_1\end{array}\right)+A_3^{\mathrm{\prime }}\left(\begin{array}{c}{x}_2-x_L\\ y_2-y_L\\ z_2\\ t=t_1\end{array}\right)\right]U\left(t_1-{\xi }_2\right)\\ V_m^2\left(\begin{array}{c}{x}_2\\ y_2\\ z_2\\ t=t_2\end{array}\right)=\left[A_1^{\mathrm{\prime }}\left(\begin{array}{c}{x}_2-x_L\\ y_2-y_L\\ z_2\\ t=t_2\end{array}\right)+A_2^{\mathrm{\prime }}\left(\begin{array}{c}{x}_2-x_L\\ y_2-y_L\\ z_2\\ t=t_2\end{array}\right)+A_3^{\mathrm{\prime }}\left(\begin{array}{c}{x}_2-x_L\\ y_2-y_L\\ z_2\\ t=t_2\end{array}\right)\right]U\left(t_2-{\xi }_2\right)\\ ⋮\\ V_m^2\left(\begin{array}{c}{x}_2\\ y_2\\ z_2\\ t=t_n\end{array}\right)=\left[A_1^{\mathrm{\prime }}\left(\begin{array}{c}{x}_2-x_L\\ y_2-y_L\\ z_2\\ t=t_n\end{array}\right)+A_2^{\mathrm{\prime }}\left(\begin{array}{c}{x}_2-x_L\\ y_2-y_L\\ z_2\\ t=t_n\end{array}\right)+A_3^{\mathrm{\prime }}\left(\begin{array}{c}{x}_2-x_L\\ y_2-y_L\\ z_2\\ t=t_n\end{array}\right)\right]U\left(t_n-{\xi }_2\right)\end{array}\right.\end{array}$(4) where $\begin{array}{rl}{\xi }_1& =\frac{\sqrt{{\left(x_1-x_L\right)}^2+{\left(y_1-y_L\right)}^2+{z_1}^2}}{c}\end{array}$ $\begin{array}{rl}{\xi }_2& =\frac{\sqrt{{\left(x_2-x_L\right)}^2+{\left(y_2-y_L\right)}^2+{z_2}^2}}{c}\end{array}$ (5) $\begin{array}{r}\left\{\begin{array}{c}{V}_m^1\left(\begin{array}{c}{x}_1\\ y_1\\ z_1\\ t=t_1\end{array}\right)-\left[A_1^{\mathrm{\prime }}\left(\begin{array}{c}{x}_1-x_L\\ y_1-y_L\\ z_1\\ t=t_1\end{array}\right)+A_2^{\mathrm{\prime }}\left(\begin{array}{c}{x}_1-x_L\\ y_1-y_L\\ z_1\\ t=t_1\end{array}\right)+A_3^{\mathrm{\prime }}\left(\begin{array}{c}{x}_1-x_L\\ y_1-y_L\\ z_1\\ t=t_1\end{array}\right)\right]U\left(t_1-{\xi }_1\right)=0\\ V_m^1\left(\begin{array}{c}{x}_1\\ y_1\\ z_1\\ t=t_2\end{array}\right)-\left[A_1^{\mathrm{\prime }}\left(\begin{array}{c}{x}_1-x_L\\ y_1-y_L\\ z_1\\ t=t_2\end{array}\right)+A_2^{\mathrm{\prime }}\left(\begin{array}{c}{x}_1-x_L\\ y_1-y_L\\ z_1\\ t=t_2\end{array}\right)+A_3^{\mathrm{\prime }}\left(\begin{array}{c}{x}_1-x_L\\ y_1-y_L\\ z_1\\ t=t_2\end{array}\right)\right]U\left(t_2-{\xi }_1\right)=0\\ ⋮\\ V_m^1\left(\begin{array}{c}{x}_1\\ y_1\\ z_1\\ t=t_n\end{array}\right)-\left[A_1^{\mathrm{\prime }}\left(\begin{array}{c}{x}_1-x_L\\ y_1-y_L\\ z_1\\ t=t_n\end{array}\right)+A_2^{\mathrm{\prime }}\left(\begin{array}{c}{x}_1-x_L\\ y_1-y_L\\ z_1\\ t=t_n\end{array}\right)+A_3^{\mathrm{\prime }}\left(\begin{array}{c}{x}_1-x_L\\ y_1-y_L\\ z_1\\ t=t_n\end{array}\right)\right]U\left(t_n-{\xi }_1\right)=0\\ V_m^2\left(\begin{array}{c}{x}_2\\ y_2\\ z_2\\ t=t_1\end{array}\right)-\left[A_1^{\mathrm{\prime }}\left(\begin{array}{c}{x}_2-x_L\\ y_2-y_L\\ z_2\\ t=t_1\end{array}\right)+A_2^{\mathrm{\prime }}\left(\begin{array}{c}{x}_2-x_L\\ y_2-y_L\\ z_2\\ t=t_1\end{array}\right)+A_3^{\mathrm{\prime }}\left(\begin{array}{c}{x}_2-x_L\\ y_2-y_L\\ z_2\\ t=t_1\end{array}\right)\right]U\left(t_1-{\xi }_2\right)=0\\ V_m^2\left(\begin{array}{c}{x}_2\\ y_2\\ z_2\\ t=t_2\end{array}\right)-\left[A_1^{\mathrm{\prime }}\left(\begin{array}{c}{x}_2-x_L\\ \begin{array}{c}{y}_2-y_L\\ z_2\\ t=t_2\end{array}\end{array}\right)+A_2^{\mathrm{\prime }}\left(\begin{array}{c}{x}_2-x_L\\ y_2-y_L\\ z_2\\ t=t_2\end{array}\right)+A_3^{\mathrm{\prime }}\left(\begin{array}{c}{x}_2-x_L\\ y_2-y_L\\ z_2\\ t=t_2\end{array}\right)\right]U\left(t_2-{\xi }_2\right)=0\\ ⋮\\ V_m^2\left(\begin{array}{c}{x}_2\\ y_2\\ z_2\\ t=t_n\end{array}\right)-\left[A_1^{\mathrm{\prime }}\left(\begin{array}{c}{x}_2-x_L\\ y_2-y_L\\ z_2\\ t=t_n\end{array}\right)+A_2^{\mathrm{\prime }}\left(\begin{array}{c}{x}_2-x_L\\ y_2-y_L\\ z_2\\ t=t_n\end{array}\right)+A_3^{\mathrm{\prime }}\left(\begin{array}{c}{x}_2-x_L\\ y_2-y_L\\ z_2\\ t=t_n\end{array}\right)\right]U\left(t_n-{\xi }_2\right)=0\end{array}\right.\end{array}$(5)

In Equation (5), the known parameters are ${\text{V}}_{\text{m}}^1$ and ${\text{V}}_{\text{m}}^2$ at different time steps and coordinates in both distribution lines. Moreover, the unknown parameters are those of the current function and model as well as coordinates of the lightning striking points. In order to solve the above-mentioned nonlinear equation system (equation 5), various numerical techniques can be applied. In this research, the PSO was used to solve this system. In this method, the absolute values of the left-hand equations in Equation (5) were set as objectives. The roots of the equation were determined by the PSO algorithm [23,24] to minimize all equations by a multi-objective approach; the main goal of minimization is zero on the right side of Equation (5). Noted that the PSO searching algorithm had been proposed based on the social behaviour of a flock of flying birds during their search for food while the position and velocity of the swarm particles are strangely depend on the cooperative communication among all the particles and each individual’s own experience at the same time. This algorithm had been used in different engineering problems to minimize the objective functions. In order to minimize the function, each PSO particle represents a candidate potential solution with a velocity vector and a position vector. The advantages of this algorithm are insensitive to scaling of design variables, simple implementation, easily parallelized for concurrent processing, very few algorithm parameters and very efficient global search algorithm however it has slow convergence in the local search. The unknown parameters of equation 5 are $I_{01},I_{02},{\mathrm{\Gamma }}_1,{\mathrm{\Gamma }}_2,{\mathrm{\Gamma }}_3,{\mathrm{\Gamma }}_4,n,v_f,\lambda ,X_L\text{and}{Y}_L.$ Other optimization algorithms of the same approach can be applied to whatever convergence rate of PSO is an acceptable level.

Results and discussion

In order to evaluate and validate the presented method, the following steps had been considered;

First step – Two samples of lightning current (based on measured currents from previous researches) were utilized and the values of lightning-induced voltage at two different distances with respect to striking point were estimated.

Second step – The evaluated induced voltage signals on two different overhead lines were used as input data to the proposed algorithm.

Third step – Lightning location and lightning current were determined by the presented method

Forth step – Both evaluated lightning location and lightning current were compared with the corresponding utilized ones (in first step).

Results are herein discussed accordingly. In the first step, Based on the geometry of problem for first case as shown in Figure 3, two signals of lightning-induced voltage (from two overhead lines with a 10 m height above the ground surface)were used to feed the proposed method as illustrated in Figures 4 and 5.

Figure 3. Geometry of problem for first case.

Figure 4. First LIV on the distribution line with a 10 m height (based on first sample of LC).

Figure 5. Second LIV on the distribution line with a 10 m height (based on first sample of LC).

With regard to the coordinates of the measuring stations and the measured values of induced voltage at both stations(as illustrated in Figures 4 and 5), the evaluated lightning return stroke current was compared with the corresponding reference one (Figure 6) [22]. Noted that actual current (reference current) was utilized in the first step of validation. The detailed parameters of reference current based on the sum of two Heidler functions (Equation 1) are $I_{01}=19.5\text{kA},I_{02}=12\text{kA},{\mathrm{\Gamma }}_1=1\mu \mathrm{s},{\mathrm{\Gamma }}_2=2\mu \mathrm{s},{\mathrm{\Gamma }}_3=8\mu \mathrm{s},{\mathrm{\Gamma }}_4=30\mu \mathrm{s},n=2$. Moreover, the lightning location was determined with radial distances of 51.66m and 1502.2m with respect to the first and second stations, respectively (based on Figure 3). The error percentage was <1% with respect to actual location. It should be noted that the actual and evaluated lightning locations were (200,800) and (198,799), respectively. In addition, real radial distance with respect to the first measuring station was 50m (248,786) while the actual value for the second station was 1500 m (1640, 1220). In this case, the non linear equations system obtained from Equation (5) (at different time steps) had been solved using PSO method and the current and lightning location parameters had been evaluated. The cost function value at last stage was 0.00002.

Figure 6. Comparison between the evaluated return stroke current and the corresponding reference one based on the first sample of LC.

In order to consider the performance of the presented algorithm based on geometry of problem (see Figure 7), two other induced voltage signals from the two overhead lines with 7.5 and 10m heights with respect to ground surface were used as input data as demonstrated in Figures 8 and 9, respectively. Evaluated LC at the channel base and the corresponding reference one are shown in Figure 10 [22]. They were in a good agreement. Moreover, the determined radial distances with respect to the first and second stations were presented as 502.2 and 2499 m, respectively (based on Figure 7) The maximum error percentage of <1% was in an acceptable range with respect to previous research. Note that the actual and evaluated lightning locations were (1100, 1500) and (1102, 1499). Likewise, the real radial distances with respect to the first measuring station were 500 m (700, 1800)while the actual value for the second station was 2500 m (3100, 3000). In this case, the cost function value of PSO at last stage was 0.000142.

Figure 7. Geometry of problem for second case.

Figure 8. First LIV on the distribution line with a 7.5 m height (based on the second sample of LC).

Figure 9. Second LIV on the distribution line with a 10 m height (based on the second sample of LC).

Figure 10. Comparison between the evaluated return stroke current and the corresponding reference one (based on the second sample of LC).

Results showed that the proposed method can evaluate LC wave shape as well as lightning location in an acceptable range of accuracy. The proposed method is based on measured induced voltage from the phase conductor of distribution overhead lines. All field components (static, induction, and radiation) were taken into account indirectly through induced voltage values in tandem with a number of inverse procedure algorithms. These are algorithms that only consider the radiation field component that is effective at far distances to the striking point. In addition, the considered method can evaluate the wave shape of LC straight forwardly in the time domain, at variance with a number of previous studies where the current peak of current values in a certain frequency range were estimated. It should be mentioned noted here that the evaluated LC wave shapes can be helpful in developing lightning protection standards depending on local data. On the other hand, the proposed method is depending on measured LIV from existing distribution lines; given the vastness of distribution networks, the system operation is easily accessible compared with field base methods. The considered method can be applied on a greater number of measuring stations and can be placed in a network throughout the considered district. In this research, Equation (1) was utilized to express the generalized form of LC function. The MTLE current model was applied to present current behaviour at various altitudes along the channel as it maintained a good range of accuracy in previous studies. However, the proposed method is very flexible and other current functions and models can be updated with the same approach.

Conclusion

In this research, a new algorithm is proposed to assign LL and LC at the channel base based on measured LIV on overhead distribution networks (in the proximity to a channel base). The proposed method was evaluated and validated using two samples of return stroke currents. Based on the proposed algorithm, all different components of the field were taken into consideration indirectly through LIV and were used for various spaces with respect to the striking point. Furthermore, it can evaluate the wave shape of LC straightforwardly in the time domain as opposed to the findings in a number of previous studies evaluating current peak or current values in a certain frequency range. The considered method is depending on measured LIV on a distribution line. The system operation is easily accessible compared with field-based methods. This study can be applied to different current models and functions and the obtained results can be used for preparing data banks on lightning flash maps and lightning flash density maps for risk assessment studies. They can also enhance lightning current database preparation to improve existing lightning protection schemes and standards based on local information.

Disclosure statement

No potential conflict of interest was reported by the author(s).