Precise Lightning Strike Detection in Overhead Lines Using KL-VMD and PE-SGMD Innovations

1. Introduction

High-voltage overhead lines serve as the primary conduits for the transmission of electric power and constitute a crucial component of the power grid [1,2]. These entities are situated in expansive regions with diverse natural and geographical features, rendering them particularly vulnerable to lightning strikes. More than 50% of power system accidents in Japan result from lightning strikes on transmission lines. In nations such as the former Soviet Union and the United States, 60% of incidents involving 275-500 kV overhead lines are attributed to lightning strikes. In China, lightning strikes account for 40% to 70% of all trip accidents [3]. The transient stage of the fault contains abundant fault information [4]. Accurately detecting signal singularities and extracting feature quantities from complex transient components, and promptly and precisely identifying the type of fault, are crucial for ensuring the stable operation of the power system.

Concurrently, a significant body of research has been dedicated to the development of feature extraction and identification methods for lightning strikes on transmission lines. Identification methods primarily rely on two aspects: the high-frequency transient components and the amplitude of lightning. Depending on the approach to feature extraction, the primary methods for identifying lightning strikes include the amplitude method, time-domain waveform method, modal extremum method, wavelet energy spectrum method, and high-frequency attenuation ratio method.

Zou et al. [5] suggests that the mutation points of transient traveling waves correspond one-to-one with the wavelet transform modulus maxima points. The strength of the traveling wave changes is characterized by the magnitude of the wavelet modulus, whereas the direction of signal change is indicated by the polarity of the wavelet changes. Consequently, the wavelet transform modulus maxima can be utilized for the identification of the type of lightning strike. Shu et al. [6] employs the ratio of the maximum initial modulus to the maximum second modulus of transient traveling waves to distinguish between non-fault lightning strike current waves that are not truncated by the fault point and fault traveling waves that are truncated by the grounding point. This approach aids in the identification of non-fault and fault lightning strikes. Sima et al. [7] utilizes wavelet transform modulus maxima for the identification of the polarity of initial transient voltage surges and the duration of surge arrester discharge. This is achieved by analyzing the polarity of wavelet modulus maxima and discharge time to detect winding strike and counterstrike faults. Guo et al. [8] differentiates between short-circuit and lightning faults by analyzing their high-frequency distribution. This is achieved by using the modulus maxima of the detail coefficients obtained through wavelet transform as criteria. These identification methods [5,6,7,8] exhibit high reliability and remain unaffected by the initial phase angle of the fault, the type of fault, and the transient fault resistance. However, the identification of modulus maxima presents a challenge and lacks a standardized approach.

Zhong et al. [9] utilizes mathematical morphology for the analysis of variations in fault waveforms. It compares the clustering effects of single-scale and multi-scale fractal curves to classify and identify lightning faults. Si et al. [10] employs a linear regression to model the disturbed traveling wave by analyzing the variations in the waveform of the transient traveling wave over time. The fit is satisfactory for cases of faults but inconsistent for non-fault lightning strikes. The waveform's consistency coefficient is utilized for assessing the occurrence of a fault on the line; however, the establishment of its threshold value lacks theoretical justification. Gao et al. [11] employs wavelet transform for multi-scale decomposition and utilizes the proportion of the transient signal energy frequency band characteristic values as identification criteria. The wavelet transform encounters challenges in selecting the optimal wavelet basis and is vulnerable to noise, which can result in potentially inaccurate outcomes. Chen et al. [12] employs the Hilbert-Huang Transform method to classify lightning types by analyzing the energy distribution across various frequency bands. This approach comprises two essential stages: Empirical Mode Decomposition (EMD) and Hilbert Transform. However, EMD decomposition process is subject to phenomena such as over-enveloping, under-enveloping, and modal mixing. Zhao et al. [13] integrates Fluorescence-Imaged MicroDeformation (FIMD) and Hilbert Transform for the detection of line traveling waves, however, the signal decomposition process is susceptible to modal mixing phenomena. Han et al. [14] differentiates non-fault lightning strikes from short-circuit faults by assessing whether the product of the time difference between the first two wavefronts at the reference side bus and the wave speed exceeds the total length of the line. However, this approach does not account for situations in which the short-circuit point or lightning strike point is located at the midpoint of the line.

In response to the identified limitations and the requirement for manual configuration of the modal components k and penalty factor α [15] in Variational Mode Decomposition (VMD) process, VMD has been enhanced through the utilization of Kullback-Leibler (KL) divergence. The optimized VMD algorithm not only effectively mitigates modal mixing phenomena but also exhibits improved resilience to noise.

This study examines the zero-mode voltage of 220 kV overhead lines under the influence of winding strike, counterstrike, and short-circuit faults following phase-mode transformation decoupling. The analysis employs KL-VMD and Permutation Entropy-Symmetric Geometry Mode Decomposition (PE-SGMD). In light of this, a proposed criterion is capable of distinguishing between winding strike, counterstrike, and short-circuit faults. The criterion has been validated using MATLAB to extensively identify Power Systems Computer Aided Design (PSCAD) simulation data, demonstrating its effectiveness in distinguishing different types of line faults.

2. Lightning Simulation Model

2.1. 220 kV Overhead Line Model

A specific 220 kV overhead line encounters a lightning fault, spanning a total length of 156 km and comprising 480 towers, with an average spacing of 325 m between them. The simulation model is depicted in Figure 1. The system employs a Phase model, encountering a lightning fault at a distance C from the 50 km end busbar, with dual lightning lines installed along the line.

2.2. Lightning Current Model

Empirical data indicates that a substantial proportion, ranging from 70% to 90%, of natural lightning carries a negative charge. Therefore, a 1.2/50 μs negative polarity double exponential lightning current model is employed for analysis [16]. Equation (1) represents the calculation for lightning current, while its waveform is illustrated in Figure 2.

$$i\left(t\right)=A{I}_0\left(e^{-\beta t}-e^{-\delta t}\right)$$

(1)

In Equation (1): A is the lightning current amplitude correction factor; I₀ is the maximum value of the lightning current; and β and δ are the corresponding waveform decay coefficients.

Figure 3 illustrates overhead lines experiencing winding strikes and counterstrikes. In this representation, I_L denotes a controlled current source, while Z₀ denotes the lightning channel's equivalent wave impedance, which is specified as 300 Ω [17].

2.3. Tower Model

Various simulation models for transmission towers are commonly used, including the standard model, inductive model, and wave impedance model. The standard model does not account for the phenomenon of overvoltage change, while the inductive model overlooks the influence of the tower on the lightning current wave [18]. This study employs a multi-wave impedance model derived from Bergeron to replicate the behavior of 220 kV cup-type iron towers. The model accurately reproduces the refraction and reflection of lightning current waves on the tower, as depicted in Figure 4. In this context, Z_H denotes the transverse carrier wave impedance, Z_M denotes the main brace wave impedance of the tower, and R_g represents the grounding resistance, which is set at 10 Ω.

2.4. Insulator Flashover Model

To minimize the disparity between simulation outcomes and real-world scenarios, the intersection method (involving the intersection of the overvoltage waveform curve with the volt-second characteristic curve) is selected for the characterization of insulator flashover. A capacitor paralleled with a voltage-controlled switch is used to simulate the flashover process [19], as depicted in Equation (2):

$$u\left(t\right)=400L_X+\frac{710L_X}{t^{0.75}}$$

(2)

In Equation (2): L_X represents the length of the insulator string.

3. VMD Algorithm

The parameters for VMD are frequently established using empirical values or the center frequency observation method. KL divergence is capable of indicating the similarity between two signals, and is therefore employed as the objective function for optimizing VMD algorithm. In accordance with modulation standards, IMFs are characterized as amplitude-frequency modulated signals, as depicted in Equation (3):

$$u_k\left(t\right)=A_k\left(t\right)cos\left[{\phi }_k\left(t\right)\right]$$

(3)

In Equation (3): φ_k(t) is the phase function and A_k(t) is the instantaneous amplitude.

The original signal is decomposed into k intrinsic mode components known as Intrinsic Mode Functions (IMFs), along with the amplitude u_k(t) and center frequency ω_k(t) of each IMF. The bandwidth of each modal signal is estimated through the application of Gaussian smoothness [20], and the constraint problem is expressed as follows:

$$\left\{\begin{array}{c}\underset{\left\{u_k\right\},\left\{{\omega }_k\right\}}{min}\left\{{\displaystyle {\displaystyle \sum }_k}|\left|{\partial }_t\left[\left(\delta \left(t\right)+j/\pi t\right)\times u_k\left(t\right)\right]e^{-j{\omega }_{k}t}\right|{|}_2^2\right\}\hfill \\ s.t.{\displaystyle {\displaystyle \sum }_k}{u}_k=f\end{array}\right.$$

(4)

In Equation (4): {u_k} are the k decomposed modal components of the signal and {ω_k} are the central frequencies of each component.

To determine the optimal solution for constrained variational mode mingling, a quadratic penalty factor α and Lagrange multipliers λ are introduced, leading to the equation of an augmented Lagrangian function [21] as follows:

$$\begin{array}{c}L\left(\left\{u_k\right\},\left\{{\omega }_k\right\},\lambda \right)=|\left|f\left(t\right)-{\displaystyle {\displaystyle \sum }_k}{u}_k\left(t\right)\right|{|}_2^2+〈\lambda \left(t\right),f\left(t\right)-{\displaystyle {\displaystyle \sum }_k}{u}_k\left(t\right)〉+\hfill \\ \alpha {\displaystyle {\displaystyle \sum }_k}|\left|{\partial }_t\left[\left(\delta \left(t\right)+j/\pi t\right)*u_k\left(t\right)\right]e^{-j{\omega }_{k}t}\right|{|}_2^2\hfill \end{array}$$

(5)

In Equation (5): * is the convolution operator. Alternation Direction Method of Multipliers (ADMM) is a technique that involves iteratively updating {u_k} and {ω_k} in order to seek the Lagrange saddle point, until the iteration termination criteria are satisfied. This process is expressed as follows:

$$u_k^{n+1}={\left|\left|f\left(t\right)-{\displaystyle \sum }_k{u}_k\left(t\right)+\lambda \left(t\right)/2\right|\right|}_2^2+\underset{u_k\in X}{argmin}\left\{\alpha {\left|\left|{\partial }_t\left[\left(\delta \left(t\right)+j/\pi t\right)\times u_k\left(t\right)\right]e^{-j{\omega }_{k}t}\right|\right|}_2^2\right\}$$

(6)

$${\omega }_k^{n+1}=\underset{{\omega }_k\in X}{argmin}\left\{{\left|\left|{\partial }_t\left[\left(\delta \left(t\right)+j/\pi t\right)\times u_k\left(t\right)\right]e^{-j{\omega }_{k}t}\right|\right|}_2^2\right\}$$

(7)

KL divergence, also referred to as relative entropy, quantifies the difference between two probability distributions P and Q [22]. It is defined in Equation (8):

$$D_{KL}(P||Q)={\displaystyle {\displaystyle \sum }_i}P\left(i\right)\log \frac{P\left(i\right)}{Q\left(i\right)}$$

(8)

The probability density distribution of a signal is determined by applying the kernel density estimation method to the probability distribution and probability density distribution functions of two signals, as depicted in Equation (9):

$$P\left(x\right)=\frac{1}{nh}{\displaystyle {\displaystyle \sum }_{i=1}^n}k\left[\frac{x_i-x}{h}\right],x\in R$$

(9)

$$k\left(u\right)=\frac{1}{\sqrt{2\pi }}{e}^{-u^2/2}$$

(10)

In Equation (9): P(x) is the probability density function and h is the window parameter, also referred to as the smoothing parameter. Equation (10) defines the Gaussian kernel function. The process by which KL divergence enhances the performance of VMD is depicted in Figure 5.

4. PE Algorithm

4.1. Introduction to PE

PE was introduced by Bandt et al. [23] as a technique for identifying randomness and dynamical variations in time series data. Principal Component Analysis (PCA) is recognized for its simplicity, rapid computational speed, and robust resistance to interference, rendering it especially well-suited for nonlinear data. The computational method can be outlined in Equation (11).

Upon conducting phase space reconstruction [24], the time series {x(i), i = 1, 2, ..., N} of length N transforms into the following time series:

$$\left.\begin{array}{c}X\left(1\right)=\left\{x\left(1\right),x\left(1+\lambda \right),\cdots ,x\left(1+\left(m-1\right)\lambda \right)\right\}\hfill \\ X\left(k\right)=\left\{x\left(i\right),x\left(i+\lambda \right),\cdots ,x\left(i+\left(m-1\right)\lambda \right)\right\}\hfill \\ ⋮\hfill \\ X\left(N-\left(\mathrm{m}-1\right)\lambda \right)=\left\{x\left(N-\left(m-1\right)\lambda \right),x\left(N-\left(m-2\right)\lambda \right),\cdots ,x\left(N\right)\right\}\hfill \end{array}\right\}$$

(11)

In Equation (11): m is the embedding dimension and λ is the time delay.

The m vectors of X(i), denoted as X(i) = {x(i), x(i+λ), ..., x(i+(m-1)λ)}, are rearranged in ascending order as follows:

$$X\left(i\right)=\{x\left(i+\left(j_1-1\right)\lambda \right)\le x\left(i+\left(j_2-1\right)\lambda \right)\le \cdots \le x(i+\left(j_m-1\right)\lambda \}$$

(12)

If x(i+(j_i1-1)λ) = x(i+(j_i2-1)λ) exists, the values are arranged based on the magnitude of j, such that when j_k1 < j_k2, it implies: x(i+(j_a-1)λ ≤ x(i+(j_a-1)λ).

Therefore, any vector X(i) has the potential to generate a set of symbolic sequences:

$$S\left(g\right)=\left[j_1,j_2,\cdots ,j_m\right]$$

(13)

In Equation (13): g = 1, 2, ..., k, k ≤ m!. The m symbols [j₁, j₂, ..., j_m] amount to m factorial permutations. According to various permutations [25], there are m! corresponding permutations. There are various symbolic sequences, with S(g) representing just one of m!. The probability of occurrence for each type of symbolic sequence is calculated through $P_1,P_2,\cdots ,P_k,{{\displaystyle \sum }}_{g=1}^k{P}_g=1$.

PE of the time series {x(i), i = 1, 2, ..., N} can be expressed in terms of Shannon entropy as follows [26]:

$$H_p\left(m\right)=-{{\displaystyle \sum }}_{g=1}^k{P}_g\ln P_g$$

(14)

When P_g=1/m!, the function H_p(m) reaches its maximum value at In(m!). PE H_p(m) can be normalized using the In(m!) transformation.The expression is presented as follows:

$$H_p=H_p\left(m\right)/In\left(m!\right)$$

(15)

The range of H_p values falls within the interval 0 ≤ H_p ≤ 1. The magnitude of H_p serves as an indicator of the level of randomness present in the time series. A higher H_p value indicates a greater degree of randomness in the time series, while a lower H_p value suggests a more predictable pattern in the time series.

4.2. Selection of PE Parameters

In the calculation of PE, it is necessary to determine three parameters: the length of the time series N, the embedding dimension m, and the time delay λ. The selection of the embedding dimension, denoted as m, typically falls within the range of 3 to 7. If m is too small, the resulting reconstructed vector will contain an insufficient number of states, thereby undermining the meaningfulness of the algorithm. Conversely, if m is too large, the reconstruction of the phase space will lead to a homogenization of the time series, resulting in increased computational demands and an inability to capture subtle changes in the series. As a result, this paper chooses m = 5. The time delay parameter λ exerts a negligible influence on the calculation of the time series. In this study, a value of λ = 1 is selected.

4.3. Symplectic Geometry Mode Decomposition

4.3.1. Phase Space Reconstruction

According to Takens' embedding theorem, a one-dimensional discrete original signal s = (x₁, x₂, ..., x_n) [27] can be topologically equivalent to the construction of a multidimensional time series matrix, denoted as the trajectory matrix X in Equation (16):

$$X=\left[\begin{array}{cccc}{x}_1& x_{1+\tau }& \cdots & x_{1+\left(d-1\right)\tau }\\ x_2& x_{2+\tau }& \cdots & x_{2+\left(d-1\right)\tau }\\ ⋮& ⋮& & ⋮\\ x_m& x_{m+\tau }& \cdots & x_{m+\left(d-1\right)\tau }\end{array}\right]$$

(16)

In Equation (16): d is the difference dimension and τ is the delay time. The trajectory matrix X is characterized by two primary parameters, d and τ. Varied values of d and τ can yield distinct matrices X, thereby exerting a substantial influence on the the analytical results.

4.3.2. Symplectic Orthogonal Matrix QR Decomposition

Let A = X^T X construct M as a Hamilton matrix [28], namely:

$$M=\left[\begin{array}{cc}{A}^T& 0\\ 0& -A\end{array}\right]$$

(17)

Let N = M², then both M and N are Hamilton matrices. A symplectic orthogonal matrix Q is to be constructed such that:

$$Q^{T}NQ=\left[\begin{array}{cc}B& R\\ 0& B^T\end{array}\right]$$

(18)

In Equation (18): B is an upper triangular matrix, b_ij = 0(i > j+1). The transformation of N using the defined symplectic orthogonal type for 2 yields B, with eigenvalues λ₁, λ₂, ..., λ_d. According to the given properties, ${\sigma }_i=\sqrt{{\mathsf{\lambda }}_i}$, where Q_i represents the eigenvectors of matrix A corresponding to σ_i. By constructing transformation coefficients $S_i=Q_i^T{X}^T$ using the matrix Q_i and the trajectory matrix X, and defining Y = Q_iS_i, the matrix X_i = Y^T can be obtained. The original matrix X for phase space reconstruction comprises d components, which can be expressed as X = X₁+X₂+...+X_d.

4.3.3. Diagonal Averaging

This step involves transforming each component X_k (1 ≤ k ≤ d) into a time series of length n, resulting in d sets of time series, each with a length equivalent to the original time series n.

Let X_m×d = (x_ij)_m×d, d* = min(m,d), m* = max(m,d), n = m+(d-1)τ.

$$y_k=\left\{\begin{array}{cc}\frac{1}{k}{\displaystyle {\displaystyle \sum }_{p=1}^k}{y}_{p,k-p+1}^{*},& 1\le k\le d*\\ \frac{1}{d}{\displaystyle {\displaystyle \sum }_{p=1}^{d*}}{y}_{p,k-p+1}^{*},& d*\le k\le m*\\ \frac{1}{n-k+1}n{\displaystyle {\displaystyle \sum }_{p=k-m*+1}^{n-m*+1}}{y}_{p,k-p+1}^{*},& m*<k\le n\end{array}\right.$$

(19)

Equation (19) yields a time series Y_k = (y₁, y₂, ..., y_n) based on X_i, leading to d groups of independent superimposed components in the time domain, specifically d of Symplectic Geometry Component (SGC) components.

4.4. PE Improved SGMD

Initially, an empty array z is initialized. The multiscale PE function is employed for computing PE of every row in SGC matrix, and the outcomes are saved in z array. Subsequently, the values in array z are sorted, and the resulting sorted indices are retained. The minimum and maximum values in array z are calculated and then assigned to d₁ and d₂, respectively. The value of d is determined by the equation d = (d₂ - d₁)/4. Four arrays, namely s₁, s₂, s₃, and s₄, are initialized and populated based on the conditions associated with the value of z. These arrays store the sums of rows of SGC matrix. As z values are iterated through, the corresponding rows of SGC are added to the respective arrays s based on the specified conditions. The arrays are stacked vertically and transposed to derive x array, which is subsequently assigned to the output variable.

5. Method for Identifying Lightning Faults

5.1. Identification of Short Circuit and Direct Strike Faults

The study primarily focuses on direct strikes as induced lightning poses minimal threat to 220 kV overhead lines. Utilizing PSCAD software for simulating the model depicted in Figure 1, single-phase short circuits, winding strikes, and counterstrike faults are configured to take place at a distance of 50 km from C end of the line. The voltage traveling waves in three phases, measured at C end, are depicted in Figure 6:

The simulated waveforms were observed, and the following differences were illustrated as follows:

(1) As depicted in Figure 6(a), in the event of a single-phase ground short circuit experienced by phase A, its voltage amplitude is lower than that of phases B and C, with a difference of approximately 100 kV. The voltage value of the short-circuit process fluctuates in accordance with a sine function.

(2) As depicted in Figure 6(b) and 6(c), the voltage amplitude of winding strike faults exceeds that of counterstrike faults for the lightning current in Figure 2. Specifically, the first wavefront voltage amplitudes reach approximately 2800 kV and 500 kV, respectively, and both exhibit steep shapes. When lightning bypasses lightning rods and poles and directly strikes the phase A conductor, the polarity of its voltage traveling wave is opposite to that of the other two phases. When the lightning rod or pole is directly struck by lightning, the transient voltage traveling wave polarities of all three phases are identical.

The conclusions drawn are based on the comparison of time-domain waveforms, and it is important to acknowledge that there are inherent limitations in practical applications. For instance, variations in fault distances may potentially affect the accuracy of identification results. This paper examines the transient energy ratio of high-frequency components of zero-mode voltage from a theoretical standpoint, and subsequently proposes a more comprehensive method for identification.

The electromagnetic coupling among the three-phase transmission lines necessitates the use of Karhunen-Loeve phase-mode transformation to mitigate its impact [29], as demonstrated in Equation (20).

$$\left(\begin{array}{c}{U}_0\\ U_1\\ U_2\end{array}\right)=\frac{1}{3}\left(\begin{array}{ccc}1& 1& 1\\ 1& -1& 0\\ 1& 0& -1\end{array}\right)\left(\begin{array}{c}{U}_a\\ U_b\\ U_c\end{array}\right)$$

(20)

In Equation (20): U₀, U₁, and U₂ denote the zero-mode component, line mode 1 component, and line mode 2 component of the phase voltages U_a, U_b, and U_c, respectively, as per Karhunen-Loeve transformation.

When the lightning strike occurs near C end, there is a high detected transient energy ratio of the high-frequency component, whereas near D end, the detected transient energy ratio of the high-frequency component is low, and vice versa. In light of the unpredictable characteristics of lightning strikes on overhead lines, it has been determined that calculating the combined energy ratios from both ends is necessary to mitigate the influence of strike distance on the transient energy ratio of high-frequency components. The zero-mode voltage at both ends of the line is extracted and then subjected to KL-VMD decomposition, as depicted in Figure 7.

Through optimization, m intrinsic mode components are decomposed at D end, while n intrinsic mode components are decomposed at C end. The absolute values are computed to derive the energy distribution of the IMF (I₁, I₂, ..., I_m and F₁, F₂, ..., F_n), extracting the high-frequency components of the highest layer I_m and F_n.

The calculation involves determining the sum of the transient energy ratios of the high-frequency components of the zero-mode voltage relative entropy at both ends of the line, denoted as E_K1. The definition is as follows:

$$E_{K1}=\frac{I_m}{{{\displaystyle \sum }}_{i=1}^m{I}_i}+\frac{F_n}{{{\displaystyle \sum }}_{i=1}^n{F}_i}$$

(21)

In Equation (21): m and n are the number of KL-VMD decomposed components IMF at C and D ends, respectively; and i is the starting index.

Table 1 presents the summation of the energy ratios of the high-frequency components of the relative entropy of the zero-mode voltage at both ends of the line, denoted as E_K1, for short circuit, winding strike, and counterstrike faults at varying distances from C end. The data presented in Table 1 indicates that E_K1 value for lightning faults significantly exceeds that for single-phase ground short circuits. Consequently, a threshold of E₁ = 0.001 is established. If E₁ > E_K1, it corresponds to a short circuit fault; if E₁ is less than or equal to E_K1, it corresponds to a direct strike fault. The lightning current functions as a unipolar high-frequency impulse signal, whereas a short circuit fault can be likened to the application of a reverse voltage source. Consequently, E_K1 value for short circuit faults is lower. Extensive simulation data calculations demonstrate that the method is capable of accurately identifying short circuit faults, but it does not differentiate between winding strike and counterstrike faults.

5.2. Identification of Winding Strike and Counterstrike Faults

To accurately differentiate between winding and counterstrike faults, it is essential to employ a suitable method that effectively identifies the transient differences between the two. The SGMD method demonstrates proficiency in signal decomposition by addressing the challenge of empirical parameter selection in EEMD methods and wavelet transforms, as well as circumventing the issue of forced component in LCD methods. Consequently, it achieves enhanced accuracy in fault identification.

At the outset, the zero-mode voltage at both terminals of the line is decomposed using PE-SGMD, with the objective of optimizing the decomposition into p intrinsic mode components at the C terminal and q intrinsic mode components at the D terminal. The absolute values are computed to derive the energy distribution of the IMF (M₁, M₂, ..., M_p and W₁, W₂, ..., W_q), by extracting the high-frequency components M_p and W_q from the highest layer. Finally, the calculation involves determining the sample entropy, symplectic geometry mode decomposition, and the high-frequency component energy ratio E_K2 of the zero-mode voltage at both ends of the line, which is defined as follows:

$$E_{K2}=\frac{M_p}{{\displaystyle \sum _{i=1}^p{M}_i}}+\frac{W_q}{{\displaystyle \sum _{i=1}^q{W}_i}}$$

(22)

In Equation (22): p and q are the number of PE-SGMD decomposed components IMF at C and D ends, respectively; and i is the starting index.

The data presented in Table 2 indicates disparities in E_K2 values between winding strike faults and counterstrike faults. E_K2 values for winding strike faults are consistently below 0.08, whereas E_K2 values for counterstrike faults consistently exceed 0.08. Consequently, a threshold of E_K2 = 0.08 has been established for differentiation. The flowchart illustrating the lightning fault identification criterion can be observed in Figure 8.

6. Simulation Validation

Utilizing the overhead line model depicted in Figure 1, simulations are performed for various operational conditions, such as short circuit, winding strike, and counterstrike faults. The double exponential model for lightning current is standardized with a 1.2/50 μs and 2.6/50 μs waveform, with initial phase angles set at 0° and 90°. The amplitude range of the winding lightning current is set between 30-70 kA, while the counterstrike lightning current amplitude range is set between 50-90 kA. The symmetrical line configuration places the lightning fault point P at a distance of 2080 km from the busbar. The fault range for single-phase ground short circuit and two-phase ground short circuit is set between 2070 km. Transition resistances are established at 10 Ω, 30 Ω, 50 Ω, 70 Ω, and 100 Ω, resulting in a total of 300 sample data sets. The data is imported into MATLAB for identification, and certain calculation results are presented in Table 3.

The analysis of the data suggests that the method yields high recognition results for short circuit, winding strike, and counterstrike faults across various operating conditions. The use of KL-VMD and SE-SGMD for extracting high-frequency components, in combination with specific criteria, has been shown to facilitate the efficient and accurate identification of overhead line fault types in diverse scenarios.

7. Conclusion

This study has developed a model for identifying overhead line lightning faults by employing KL-VMD and SE-SGMD to decompose the zero-mode voltage. The study also established lightning recognition criteria through extensive data calculations, resulting in the following conclusions:

• (1) In situations where signal decomposition is hindered by challenges such as modal mingling, the utilization of KL-VMD can automatically optimize the decomposition layers and penalty factors. This approach effectively extracts transient characteristic quantities, demonstrating its strong adaptability in fault signal decomposition;
• (2) A criterion is proposed for identifying winding strike, counterstrike, and short circuit faults by analyzing the fault stage traveling wave amplitude, wavefront polarity, rate of change, and modal energy distribution using KL-VMD and PE-SGMD. Following thorough data calculations, the validity and accuracy of this criterion have been confirmed;
• (3) The criterion demonstrates high reliability in accurately distinguishing between short circuit faults and lightning conditions under various lightning current amplitudes, distances, and initial phase angles. It also provides a reference for line fault identification.