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

Abstract

When overhead lines are impacted by lightning, the traveling wave of the fault contains a wealth of fault information. The accurate extraction of feature quantities from transient components and their classification are fundamental to the identification of lightning faults. The extraction process may involve modal aliasing, optimal wavelet base issues, and inconsistencies between the lightning strike distance and the fault point. These factors have the potential to impact the effectiveness of recognition. This paper presents a method for identifying lightning strike faults by utilizing Kullback–Leibler (KL) divergence enhanced Variational Mode Decomposition (VMD) and Symmetric Geometry Mode Decomposition (SGMD) improved with Permutation Entropy (PE) to address the aforementioned issues. A model of a 220 kV overhead line is constructed using real faults to replicate scenarios of winding strike, counterstrike, and short circuit. The three-phase voltage is chosen and then subjected to Karenbaren decoupling in order to transform it into zero mode, line mode 1, and line mode 2. The zero-mode voltage is decomposed using KL-VMD and PE-SGMD methods, and the lightning identification criteria are developed based on various transient energy ratios. The research findings demonstrate that the criteria effectively differentiate between winding strike, counterstrike, and short-circuit faults, thus confirming the accuracy and efficacy of the lightning fault identification criteria utilizing KL-VMD and PE-SGMD.

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]. The accurate detection of signal anomalies and the extraction of feature from complex transient components are vital for the prompt and precise identification of the fault type in power systems. This is essential to maintain their stable operation.

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] suggest 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 the 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] propose a method for distinguishing between non-fault lightning strike currents and fault traveling waves. This method involves comparing the maximum initial modulus and the maximum second modulus of transient traveling waves. This comparison helps in identifying whether the waves are truncated by the fault point (i.e., a fault traveling wave) or not truncated (i.e., a non-fault lightning strike current). This approach aids in the identification of non-fault and fault lightning strikes. Sima et al. [7] utilize 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 the wavelet modulus maxima and discharge time to detect winding strike and counterstrike faults. Guo et al. [8] differentiate 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] utilize 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] employ 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] employ wavelet transform for multi-scale decomposition and utilize 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] employ 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, the EMD decomposition process is subject to phenomena such as over-enveloping, under-enveloping, and modal mixing. Zhao et al. [13] integrate 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] differentiate 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 the manual configuration of the modal components k and penalty factor α [15] in the 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. In Figure 1, ‘C’ represents the left system busbar, and ‘D’ represents the right system busbar. 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. The power supply voltage is 1.05 times the rated voltage, the value is 230 kV, the system sampling frequency is 1 MHz, and the simulation step is 0.02 μs. The system parameters are shown in

Table A1. The system parameters.

System Parameters	Numeric Value
power supply voltage	230 kV
sampling frequency	1 MHZ
simulation time step	0.02 μs
rise time	1.2 μs
half-value time	50 μs
lightning channel impedance	300 Ω
total length of the circuit	156 km
ground resistance	10 Ω

in Appendix A.

2.2. Lightning Current Model

The empirical data indicate 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 the lightning current, while its waveform is illustrated in Figure 2.

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

(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.

2.4. Insulator Flashover Model

To minimize the disparity between the 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 the insulator flashover. A capacitor paralleled with a voltage-controlled switch is used to simulate the flashover process [19], as depicted in Equation (2):

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

(2)

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

2.5. Simulation Model Verification

Regulatory standards are used to determine the withstand capability against winding strokes and counterstrike, with the formulas given as Equations (4) and (5).

$$U_{50\%}=533L_X+132$$

(3)

$$I_1=U_{50\%}/100$$

(4)

$$I_2=\frac{U_{50\%}}{(1-k)\beta R_g+\left(\frac{h_a}{h_t}-k\right)\beta \frac{L_t}{2.6}+\left(1-\frac{h_g}{h_c}{k}_0\right)\frac{h_c}{2.6}}$$

(5)

In Equation (5), $k_0$ is the coupling coefficient.

According to the regulatory standards, the withstand capability against winding strokes is determined to be 12 kA, while the withstand capability against back flashover is found to be 75 kA. The lightning current is conducted at intervals of 0.5 kA, observing whether flashovers occur on the insulator. The simulation results indicate that the withstand capability against winding strokes is 12.5 kA, while the withstand capability against back flashover is 77.5 kA. It complies with the regulatory standards, and the simulated three-phase voltage waveform characteristics, including polarity, amplitude, and shape, align with the descriptions of the fault waveforms in references [20,21]. Therefore, this validates the accuracy of the simulations conducted in this study.

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 the VMD algorithm. In accordance with modulation standards, IMFs are characterized as amplitude-frequency modulated signals, as depicted in Equation (6):

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

(6)

In Equation (6), φ_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 [22], and the constraint problem is expressed as follows:

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

(7)

In Equation (7), {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 [23] as follows:

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

(8)

In Equation (8), * is the convolution operator. The 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:

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

(9)

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

(10)

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

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

(11)

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 (12):

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

(12)

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

(13)

In Equation (12), P(x) is the probability density function and h is the window parameter, also referred to as the smoothing parameter. Equation (13) 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. [25] 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 [26], 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\right(1),x(1+\lambda ),\cdots ,x(1+(m-1)\lambda \left)\right\}\hfill \\ X\left(k\right)=\left\{x\right(i),x(i+\lambda ),\cdots ,x(i+(m-1)\lambda \left)\right\}\hfill \\ ⋮\hfill \\ X(N-(m-1\left)\lambda \right)=\left\{x\right(N-(m-1)\lambda ),x(N-(m-2)\lambda ),\cdots ,x(N\left)\right\}\hfill \end{array}\right\}$$

(14)

In Equation (14), 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)=\left\{x\right(i+(j_1-1)\lambda )\le x(i+(j_2-1)\lambda )\le \cdots \le x(i+(j_m-1)\lambda \}$$

(15)

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)=[j_1,j_2,\cdots ,j_m]$$

(16)

In Equation (16), g = 1, 2, ..., k, k ≤ m!. The m symbols [j₁, j₂, ..., j_m] amount to m factorial permutations. According to various permutations [27], 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,{\sum }_{g=1}^k{P}_g=1$.

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

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

(17)

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(m!)$$

(18)

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) [29] can be topologically equivalent to the construction of a multidimensional time series matrix, denoted as the trajectory matrix X in Equation (19):

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

(19)

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 analytical results.

4.3.2. Symplectic Orthogonal Matrix QR Decomposition

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

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

(20)

Let N = M², then both M and N be 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]$$

(21)

In Equation (21), 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 \sum _{p=1}^k}{y}_{p,k-p+1}^{*},& 1\le k\le d \ast \\ \frac{1}{d}{\displaystyle \sum _{p=1}^{d*}}{y}_{p,k-p+1}^{*},& d \ast \le k\le m \ast \\ \frac{1}{n-k+1}n{\displaystyle \sum _{p=k-m*+1}^{n-m*+1}}{y}_{p,k-p+1}^{*},& m \ast <k\le n\end{array}\right.$$

(22)

Equation (22) 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 the 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 the PE of every row in the SGC matrix, and the outcomes are saved in array z. 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 the SGC matrix. As z values are iterated through, the corresponding rows of SGC are added to the respective s arrays based on the specified conditions. The arrays are stacked vertically and transposed to derive the 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 v5.0.2 software for simulating the model depicted in Figure 1, single-phase short-circuit, winding strike, and counterstrike faults are configured to take place at a distance of 50 km from the C end of the line. The voltage traveling waves in three phases, measured at the 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 6a, 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.

(1)

In Figure 6b,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 [31], as demonstrated in Equation (23).

$$\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)$$

(23)

In Equation (23), 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 the C end, there is a high detected transient energy ratio of the high-frequency component. Conversely, when the lightning strike is near the 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, calculating the combined energy ratios from both ends is necessary to mitigate the influence of the 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. In Figure 7, IMF1 represents the first intrinsic mode component of zero-mode voltage decomposed by KL-VMD; IMF2 represents the second intrinsic mode component of zero-mode voltage decomposed by KL-VMD; IMF3, IMF4, and so on.

Through optimization, the m intrinsic mode components are decomposed at the D end, while the n intrinsic mode components are decomposed at the 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}{{\sum }_{i=1}^m{I}_i}+\frac{F_n}{{\sum }_{i=1}^n{F}_i}$$

(24)

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

Table 1. Summation of the energy ratio of high-frequency components of relative entropy in zero mode-voltage at both ends of the line under different conditions E_K1.

Fault Type	PC/km	PD/km	EK1
Short-Circuit Fault	30	126	1.3326 × 10−5
Short-Circuit Fault	40	116	1.3516 × 10−5
Short-Circuit Fault	50	106	1.4201 × 10−5
Short-Circuit Fault	60	96	1.5051 × 10−5
Winding Strike	30	126	0.0159
Winding Strike	40	116	0.0356
Winding Strike	50	106	0.5490
Winding Strike	60	96	0.2999
Counterstrike	30	126	0.1838
Counterstrike	40	116	0.1629
Counterstrike	50	106	0.1717
Counterstrike	60	96	0.1789

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 the C end. The data presented in Table 1 indicate that the 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, the 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_{\mathrm{p}}}{{\displaystyle \sum _{i=1}^{\mathrm{p}}{M}_i}}+\frac{W_{\mathrm{q}}}{{\displaystyle \sum _{i=1}^{\mathrm{q}}{W}_i}}$$

(25)

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

The data presented in

Table 2. Summation of the energy ratio of the high-frequency component in PE symplectic geometric modal component of zero-mode voltage at both ends of the line under different conditions E_k2.

Fault Type	PC/km	PD/km	Ek2
Winding Strike	30	126	0.0026
Winding Strike	40	116	0.0043
Winding Strike	50	106	0.0033
Winding Strike	60	96	0.0029
Winding Strike	70	86	0.0068
Counterstrike	30	126	0.8985
Counterstrike	40	116	0.9310
Counterstrike	50	106	0.6551
Counterstrike	60	96	0.8733
Counterstrike	70	86	0.4693

indicate disparities in the E_K2 values between winding strike faults and counterstrike faults. The E_K2 values for winding strike faults are consistently below 0.08, whereas the 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 the 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 and 70 kA, while the counterstrike lightning current amplitude range is set between 50 and 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 the 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 are imported into MATLAB for identification, and certain calculation results are presented in

Table 3. Partial calculation results.

Distance/km	30 km			50 km			70 km
	EK1	EK2	Result	EK1	EK2	Result	EK1	EK2	Result
Phase A ground short circuit, Rg = 30 Ω, θ = 90°	1.3326 × 10−5	—	Short Circuit	1.4201 × 10−5	—	Short Circuit	2.0681 × 10−5	—	Short Circuit
Phase A ground short circuit, Rg = 50 Ω, θ = 0°	1.6745 × 10−5	—	Short Circuit	1.8363 × 10−5	—	Short Circuit	2.7329 × 10−5	—	Short Circuit
Phase AB ground short circuit, Rg = 30 Ω, θ = 90°	2.9690 × 10−5	—	Short Circuit	2.6474 × 10−5	—	Short Circuit	3.0331 × 10−5	—	Short Circuit
Phase AB ground short circuit, Rg = 50 Ω, θ = 90°	3.1702 × 10−5	—	Short Circuit	2.8562 × 10−5	—	Short Circuit	3.0805 × 10−5	—	Short Circuit
Winding Strike, Imax = 30 kA, θ = 0°, 1.2/50 μs	0.3145	0.0035	Winding Strike	0.5705	0.0088	Winding Strike	0.5166	0.0088	Winding Strike
Winding Strike, Imax = 40 kA, θ = 0°, 2.6/50 μs	0.3168	0.0029	Winding Strike	0.5781	0.0087	Winding Strike	0.5196	0.0091	Winding Strike
Winding Strike, Imax = 50 kA, θ = 90°, 1.2/50 μs	0.0159	0.0026	Winding Strike	0.5492	0.0033	Winding Strike	0.2793	0.0068	Winding Strike
Winding Strike, Imax = 60 kA, θ = 90°, 2.6/50 μs	0.0165	0.0026	Winding Strike	0.5526	0.0039	Winding Strike	0.2817	0.0063	Winding Strike
Counterstrike, Imax = 90 kA, θ = 0°, 1.2/50 μs	0.1951	1.1181	Counterstrike	0.1906	0.8240	Counterstrike	0.1645	0.4694	Counterstrike
Counterstrike, Imax = 100 kA, θ = 90°, 2.6/50 μs	0.1726	0.5824	Counterstrike	0.1566	0.5809	Counterstrike	0.1455	0.3876	Counterstrike
Counterstrike, Imax = 110 kA, θ = 90°, 1.2/50 μs	0.1838	0.8985	Counterstrike	0.1717	0.6551	Counterstrike	0.1513	0.4693	Counterstrike
Counterstrike, Imax = 120 kA, θ = 0°, 2.6/50 μs	0.1788	0.5922	Counterstrike	0.1820	0.5773	Counterstrike	0.1607	0.3888	Counterstrike

.

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.

Selecting different lightning current amplitudes, distances, and initial phase angles, 30 sets each for winding strike faults and counterstrike faults were chosen. The results were compared with the methods in references [32,33], as shown in

Table 4. Recognition performance of different methods.

Identification Method	Number of Test Samples	Fault Detection Accuracy
		Winding Strike Fault	Counterstrike Fault
identification methods in the literature 32	60	90.0%	88.3%
identification methods in the literature 33	60	95.0%	95.0%
KL-VMD+PE-SGMD	60	96.6%	98.3%

. it can be observed that the recognition accuracy of the method proposed in this paper is higher than the other two methods in all cases. This also validates the effectiveness of the criteria used.

7. Conclusions

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, along with the 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.