Faulty Feeder Detection Method Based on VMD–FFT and Pearson Correlation Coe ﬃ cient of Non-Power Frequency Component in Resonant Grounded Systems

: Through analyzing the transient components and transient characteristics in transient zero-sequence current (TZSC), a novel fault feeder detection method based on the transient correlation of non-power frequency components (NPFCs) for the resonant grounded system is proposed. Firstly, using variational mode decomposition combined with fast Fourier transformation (VMD–FFT) to decompose the TZSC, by removing the power frequency components and noise signals, the transient NPFCs can be obtained. Secondly, to reﬂect the overall changing trend between faulty and healthy currents, the moving average ﬁlter is introduced to smooth the NPFCs; in this way, the fault transient features can be accurately revealed. Finally, the faulty feeder can be detected by comparing the threshold with the maximum di ﬀ erence value of comprehensive correlation coe ﬃ cient of NPFCs. The detection results show that the proposed fault detection method can accurately select the faulty feeder; it is una ﬀ ected by fault resistances, fault phase angles, etc. Moreover, the detection method can resist noise interference.


Introduction
Among the existing operation modes of power distribution networks, the resonance grounding system is most widely used. Single-phase ground (SPG) faults occur most frequently in this type of network, accounting for more than 70% of all fault types. When SPG faults occur in resonant grounding systems, existing line selection techniques have relatively low accuracy due to factors such as unstable arcs and weak fault characteristics [1][2][3][4].
The current signal at the moment of fault encompasses rich transient characteristic information. The transient method [11][12][13][14][15] is favored by many scholars. Liu et al. [16] used the Prony algorithm to extract the transient main frequency components of individual lines. This method requires no Energies 2020, 13 settings but has low anti-noise capacity. Shu et al. [17] used the peak-valley features of zero-sequence current morphology for faulty line selection; this approach does not work well under small-angle fault conditions. Xue et al. [18] proposed a line selection method based on the difference in projection coefficients of the zero-sequence current on the zero-sequence voltage. The selection ability for low-resistance faults showed room for improvement. Wang et al. [19] used wavelet transformation to decompose a transient zero-sequence current (TZSC), but the choice of wavelet basis function showed undue influence on the decomposition results. Wang et al. [20] introduced a novel cut-in point for fault detection issues; the power component can be calculated based on the fast Fourier transformation (FFT) method. Firstly, the power frequency component is eliminated by FFT, and then the non-power frequency component (NPFC) is obtained; the final results show that the fault detection method has high accuracy. However, when the TZSC contains harmonics, the detection method needs to be further verified. The Empirical Mode Decomposition (EMD) algorithm has been proposed for faulty line selection to resolve the above problems [21][22][23]. EMD does not require a priori selection of the basis function and can directly decompose the input signal. Shu et al. [22] first used EMD to decompose a TZSC and obtain the high-frequency components. The faulty line was determined based on the magnitude of the correlation coefficient of the waveform. Though effective to some extent, this EMD decomposition process is prone to modal aliasing.
A faulty line selection method based on variational mode decomposition combined with fast Fourier transformation (VMD-FFT) and the Pearson correlation coefficient was developed in this study to improve the accuracy of faulty line selection in resonant grounded systems. The VMD first adaptively calculates the center frequency of the signal to decompose the zero-sequence current of each outgoing line into intrinsic mode functions (IMFs), intuitively reflecting their compositions. The low-frequency IMFs of each outgoing line are then subjected to FFT analysis to obtain the power frequency component of the fault signal. Then, the NPFCs of each outgoing line are obtained, smoothed, and subjected to a Pearson correlation analysis. Finally, the faulty line is determined. This prevents the inaccurate separation of power frequency components otherwise caused by modal aliasing and boundary effects in traditional algorithms.

Composition and Characteristics of TZSC for Faults
When a SPG fault occurs in the resonant ground system, the equivalent circuit for the transient analysis is as shown in Figure 1. In this circuit, U 0 f represents the zero-sequence voltage source, R 0 and L 0 are the equivalent grounding resistance and inductance of the zero-sequence circuit, R P and L P are the equivalent resistance and equivalent inductance of the arc suppression coil, and C is the total ground capacitance of the system. The current signal at the moment of fault encompasses rich transient characteristic information. The transient method [11][12][13][14][15] is favored by many scholars. Liu et al. [16] used the Prony algorithm to extract the transient main frequency components of individual lines. This method requires no settings but has low anti-noise capacity. Shu et al. [17] used the peak-valley features of zero-sequence current morphology for faulty line selection; this approach does not work well under small-angle fault conditions. Xue et al. [18] proposed a line selection method based on the difference in projection coefficients of the zero-sequence current on the zero-sequence voltage. The selection ability for lowresistance faults showed room for improvement. Wang et al. [19] used wavelet transformation to decompose a transient zero-sequence current (TZSC), but the choice of wavelet basis function showed undue influence on the decomposition results. Wang et al. [20] introduced a novel cut-in point for fault detection issues; the power component can be calculated based on the fast Fourier transformation (FFT) method. Firstly, the power frequency component is eliminated by FFT, and then the non-power frequency component (NPFC) is obtained; the final results show that the fault detection method has high accuracy. However, when the TZSC contains harmonics, the detection method needs to be further verified.
The Empirical Mode Decomposition (EMD) algorithm has been proposed for faulty line selection to resolve the above problems [21][22][23]. EMD does not require a priori selection of the basis function and can directly decompose the input signal. Shu et al. [22] first used EMD to decompose a TZSC and obtain the high-frequency components. The faulty line was determined based on the magnitude of the correlation coefficient of the waveform. Though effective to some extent, this EMD decomposition process is prone to modal aliasing.
A faulty line selection method based on variational mode decomposition combined with fast Fourier transformation (VMD-FFT) and the Pearson correlation coefficient was developed in this study to improve the accuracy of faulty line selection in resonant grounded systems. The VMD first adaptively calculates the center frequency of the signal to decompose the zero-sequence current of each outgoing line into intrinsic mode functions (IMFs), intuitively reflecting their compositions. The low-frequency IMFs of each outgoing line are then subjected to FFT analysis to obtain the power frequency component of the fault signal. Then, the NPFCs of each outgoing line are obtained, smoothed, and subjected to a Pearson correlation analysis. Finally, the faulty line is determined. This prevents the inaccurate separation of power frequency components otherwise caused by modal aliasing and boundary effects in traditional algorithms.

Composition and Characteristics of TZSC for Faults
When a SPG fault occurs in the resonant ground system, the equivalent circuit for the transient analysis is as shown in Figure 1. In this circuit, 0 f U represents the zero-sequence voltage source, 0 R and 0 L are the equivalent grounding resistance and inductance of the zero-sequence circuit, P R and P L are the equivalent resistance and equivalent inductance of the arc suppression coil, and C is the total ground capacitance of the system. As shown in Figure 1, the circuit equation is: As shown in Figure 1, the circuit equation is: where i C and i L are the transient capacitance current and inductance current; U m is the relative ground voltage amplitude of the system; ω b is the reference frequency; and ϕ is the initial phase angle of the power supply voltage at the moment of grounding. According to Equations (1) and (2), i 0 , the TZSC flowing through the fault point can be determined as: where I L and I C are the amplitudes of the inductor current and capacitor current, ω f and δ are the transient free oscillation frequency and attenuation coefficient, δ = 1/τ c , τ c is the attenuation time constant of the capacitor current, and τ L is the attenuation time constant of the inductor current. Equation

Principles of Variational Mode Decomposition
VMD is a widely used, non-linear, pre-scaled signal processing algorithm that was first proposed in 2014. In essence, the VMD algorithm transforms the signal decomposition problem into a constrained optimization problem. The optimal solution obtained is the decomposed single-component amplitude-modulated-frequency-modulated (AM-FM) signal. In other words, VMD decomposes the TZSC signal i 0 j of line j into K intrinsic modal functions u k (t) with a center angle frequency ω k , where K is the preset number of modal components. Unlike EMD decomposition, VMD defines each IMF as an AM-FM function u k (t): where instantaneous amplitude value A k (t) ≥ 0; ω k (t) = φ k (t).
The step-wise operation of the VMD algorithm is as follows [24].
(1) The TZSC signal i 0 j is subjected to Hilbert transformation, resulting in the analytical signals of the K modal component u k (t) and the u k (t) single-sided spectrum: where δ(t) is the impulse function. (2) The frequency spectrum of each mode is modulated to the base frequency band: (3) The square L 2 norm of the gradient of Formula (6) is calculated. The bandwidth of each modal signal was estimated. The variational problem is established as follows: where ∂ t is the partial derivative of t; {u k } := {u 1 , u 2 , . . . , u K }; {ω k } := {ω 1 , ω 2 , . . . , ω K }. (4) In order to solve the above constrained variational problem, the quadratic penalty factor α and the Lagrangian multiplication operator λ(t) can be introduced into Equation (7) to turn it into a non-constrained variational problem. The extended Lagrangian expression is: The alternate direction method of multipliers (ADMM) is applied to solve the extended Lagrangian expression. The specific steps are: (a) Initialize variables û 1 k , ω 1 k ,λ 1 , n; (b) Execute cycles: n = n + 1; (c) Update u k , ω k and λ;û where τ is the noise margin parameter. (d) Steps b-d are repeated until the iteration termination conditions are met: where e is the condition for iteration termination.

of 18
In order to discuss the decomposition effect of VMD on TZSC, a previously described ideal TZSC signal [25] is taken as an example for verification: The ideal signal can be generated from the characteristics of the TZSC when a SPG fault occurs in a resonant grounding system. VMD was used here to analyze x(t) while the number of decomposition levels K = 3 and the secondary penalty factor α = 2000. The resulting decomposed IMF1, IMF2, and IMF3 are shown in Figure 2.
In order to discuss the decomposition effect of VMD on TZSC, a previously described ideal TZSC signal [25] is taken as an example for verification:  As shown in Figure 2, IMF1 is a low-frequency component containing power frequency components and attenuated DC components. IMF2 is a transient high-frequency attenuation component and IMF3 is noise signal. VMD effectively separated low-frequency components, highfrequency components, and white noise components in an ideal TZSC in this case. As shown in Figure 2, IMF1 is a low-frequency component containing power frequency components and attenuated DC components. IMF2 is a transient high-frequency attenuation component and IMF3 is noise signal. VMD effectively separated low-frequency components, high-frequency components, and white noise components in an ideal TZSC in this case.

Separation of Power Frequency Components by FFT
As shown in Figure 2, VMD cannot directly separate the power frequency components of the zero-sequence current. Further processing is needed. Steady-state harmonic signal analysis is usually conducted via FFT. In order to accurately obtain NPFCs and avoid any interference from transient signals, the second power frequency period of IMF1 was subjected here to FFT analysis. First, an FFT analysis of the second power frequency cycle of IMF1 was performed to obtain power frequency components; IMF1 (minus the obtained power frequency component) and IMF2 were then added to obtain the NPFC of the TZSC.
The actual NPFC are shown in the blue solid line in Figure 3 with the ideal zero sequence current signal as an example. An FFT analysis was performed on the IMF1 shown in Figure 2 to obtain the NPFC of the ideal zero-sequence current, as marked with a red dashed line in Figure 3.
analysis of the second power frequency cycle of IMF1 was performed to obtain power frequency components; IMF1 (minus the obtained power frequency component) and IMF2 were then added to obtain the NPFC of the TZSC.
The actual NPFC are shown in the blue solid line in Figure 3 with the ideal zero sequence current signal as an example. An FFT analysis was performed on the IMF1 shown in Figure 2 to obtain the NPFC of the ideal zero-sequence current, as marked with a red dashed line in Figure 3. In Figure 3, the actual NPFC coincide with the NPFC obtained by the proposed algorithm. This method can accurately separate the power frequency components and noise signals, ultimately revealing the NPFC of the transient zero sequence current.

Data Smoothing and Pearson Correlation Coefficient
To highlight the overall trend characteristics of NPFCs of faulty lines and non-faulty lines, the paper selects a moving average filter to smooth the NPFCs, and the results obtained can better describe the trend characteristics of NPFCs.
The NPFCs of TZSCs 0i i and 0 j i were selected for faulty and non-faulty lines, as shown in Figure 4a,c, to smooth the element point at a certain moment as an example. This point was placed as the center and g elements (odd number) were selected as an interval. (g = 199 was found by trial and error to meet the line selection requirements of this study.) The elements in the interval were fitted by the second-order polynomial to obtain a weighted regression curve; this step was repeated for the remaining element points to obtain their respective weighted regression curves. Finally, the centers of all regression curves were connected together to obtain a smoothed signal curve. In Figure 3, the actual NPFC coincide with the NPFC obtained by the proposed algorithm. This method can accurately separate the power frequency components and noise signals, ultimately revealing the NPFC of the transient zero sequence current.

Data Smoothing and Pearson Correlation Coefficient
To highlight the overall trend characteristics of NPFCs of faulty lines and non-faulty lines, the paper selects a moving average filter to smooth the NPFCs, and the results obtained can better describe the trend characteristics of NPFCs.
The NPFCs of TZSCs i 0i and i 0 j were selected for faulty and non-faulty lines, as shown in Figure 4a,c, to smooth the element point at a certain moment as an example. This point was placed as the center and g elements (odd number) were selected as an interval. (g = 199 was found by trial and error to meet the line selection requirements of this study.) The elements in the interval were fitted by the second-order polynomial to obtain a weighted regression curve; this step was repeated for the remaining element points to obtain their respective weighted regression curves. Finally, the centers of all regression curves were connected together to obtain a smoothed signal curve.  The results before and after smoothing the NPFCs are shown in Figure 4. This method appears to perform well without altering the overall trend of 0i i and 0 j i , which is conducive to subsequent correlation analysis.
The Pearson correlation coefficient is based on a concept proposed by Karl Pearson. It is typically used to measure the degree of correlation between two variables resulting in a value between −1 and The results before and after smoothing the NPFCs are shown in Figure 4. This method appears to perform well without altering the overall trend of i 0i and i 0 j , which is conducive to subsequent correlation analysis. The Pearson correlation coefficient is based on a concept proposed by Karl Pearson. It is typically used to measure the degree of correlation between two variables resulting in a value between −1 and 1 [26,27]. It is calculated as follows: where N is the number of sampling points; i 0i and i oj are the average values of N elements of i 0i and i 0 j , respectively. Pearson correlation has good analytical performance for the linear correlation. However, when signals are interfered by the noises, the strong correlation between actual signals may be broken. As a result, the Pearson correlation coefficients may be close to 0 [28]. The smoothing method can filter the high-frequency noises and restore the original characteristics of the actual signals. Therefore, the smoothing method can help to accurately reflect the true Pearson correlations between faulty NPFC and non-faulty NPFCs. The Pearson correlation coefficient of the two lines i 0i and i 0 j is −0.2901 before processing ( Figure 4) and that of ρ ij after smoothing is −0.9980. This suggests that smoothing highlights the overall trend characteristics of NPFCs in faulty and non-faulty lines, which yields an accurate correlation line selection.

Line Selection Method and Process
(1) The zero-sequence voltage of the system busbar, U 0 , is first continuously monitored on-line. When (3) This matrix M is normalized to obtain the comprehensive correlation coefficient P of each line. The comprehensive correlation coefficient P i of the i-th circuit relative to other m−1 circuits is calculated as follows: (4) The absolute value S of the difference between the maximum value P max and the minimum value P min of the comprehensive correlation coefficient is calculated. If S is greater than the set threshold P set (0.3 in this case), the corresponding line of P min is judged as faulty; otherwise, it is judged as a bus fault.
A flow chart of this faulty line selection process is given in Figure 5.    The line impedance parameters [29] are shown in Table 1. The transformers in this setup are 110/10.5 kV. The resistance of the single-phase neutral point coil on the high-voltage side is 0.40 Ω. The inductance is 12.2 Ω. The resistance of the single-phase coil on the low-voltage side is 0.006 Ω and the inductance is 0.183 Ω. The excitation current is 0.672 Ω, the excitation flux is 202.2 Wb, and the magnetic circuit resistance is 400 kΩ. The delta connection method was used in all cases, ZL = 400 + j20 Ω. Overcompensation was simulated in the arc suppression coil Energies 2020, 13, 4724 9 of 18 grounding system with a 10% overcompensation degree and 1.2819 H arc suppression coil inductance. The resistance value of the arc suppression coil is 10% of the reactance value, 40.2517 Ω. The model sampling frequency fs = 10 5 Hz. The period T = 0.02 s. The SPG fault occurred at 0.02 s in the simulation, which was 0.06 s long in total.

Simulation Modeling
A cable hybrid system simulation model was established in the electromagnetic transient simulation software ATP (3.8, SINTEF Energy Research, Norway), as shown in Figure 6, which includes four lines:  The line impedance parameters [29] are shown in Table 1. The transformers in this setup are 110/10.5 kV. The resistance of the single-phase neutral point coil on the high-voltage side is 0.40 Ω. The inductance is 12.2 Ω. The resistance of the single-phase coil on the low-voltage side is 0.006 Ω and the inductance is 0.183 Ω. The excitation current is 0.672 Ω, the excitation flux is 202.2 Wb, and the magnetic circuit resistance is 400 kΩ. The delta connection method was used in all cases, ZL = 400 + j20 Ω. Overcompensation was simulated in the arc suppression coil grounding system with a 10% overcompensation degree and 1.2819 H arc suppression coil inductance. The resistance value of the arc suppression coil is 10% of the reactance value, 40.2517 Ω. The model sampling frequency fs = 10 5 Hz. The period T = 0.02 s. The SPG fault occurred at 0.02 s in the simulation, which was 0.06 s long in total.

Overhead Line Fault
The A-phase grounding fault on the line L 1 at a distance of 5 km from the bus line is discussed here as an example. The initial phase angle of the fault is 90 • and the grounding resistance is 100 Ω. After Gaussian white noise with a signal to noise ratio (SNR) of 10 dB was injected into the collected signal, the zero-sequence current waveform of the line was derived as shown in Figure 7.

Overhead Line Fault
The A-phase grounding fault on the line L1 at a distance of 5 km from the bus line is discussed here as an example. The initial phase angle of the fault is 90° and the grounding resistance is 100 Ω. After Gaussian white noise with a signal to noise ratio (SNR) of 10 dB was injected into the collected signal, the zero-sequence current waveform of the line was derived as shown in Figure 7.   As shown in Figure 7, the zero-sequence current of each line contains obvious power frequency components and transient high-frequency components. The maximum amplitude of the faulty line reaches about 150 A, which is much larger than that of the non-faulty line. The lines have opposing polarities. The TZSC of each outgoing line was decomposed in three layers using VMD, resulting in each outgoing line IMF1, as shown in Figure 8. As shown in Figure 7, the zero-sequence current of each line contains obvious power frequency components and transient high-frequency components. The maximum amplitude of the faulty line reaches about 150 A, which is much larger than that of the non-faulty line. The lines have opposing polarities. The TZSC of each outgoing line was decomposed in three layers using VMD, resulting in each outgoing line IMF1, as shown in Figure 8. As shown in Figure 8, the IMF1 of line L1 has the largest amplitude and the IMF1 of each outgoing line is composed mainly of low-frequency components, with power frequency components as the main components. FFT analysis can be used to separate the power frequency components. As shown in Figure 8, the IMF1 of line L 1 has the largest amplitude and the IMF1 of each outgoing line is composed mainly of low-frequency components, with power frequency components as the main components. FFT analysis can be used to separate the power frequency components.
As shown in Figure 9, IMF2 is dominated by transient high-frequency components without power frequency components. This indicates that VMD can decompose low-frequency components and high-frequency components with favorable decomposition effects.
Energies 2020, 13, x FOR PEER REVIEW 11 of 19 As shown in Figure 9, IMF2 is dominated by transient high-frequency components without power frequency components. This indicates that VMD can decompose low-frequency components and high-frequency components with favorable decomposition effects. As shown in Figure 10, IMF3 is the noise signal of each line. In effect, VMD can accurately decompose the interference noise of each line and has a good denoising ability. After VMD decomposition was completed, Hj was obtained by the line selection method discussed above, as shown in Figure 11. As shown in Figure 10, IMF3 is the noise signal of each line. In effect, VMD can accurately decompose the interference noise of each line and has a good denoising ability. After VMD decomposition was completed, H j was obtained by the line selection method discussed above, as shown in Figure 11. As shown in Figure 10, IMF3 is the noise signal of each line. In effect, VMD can accurately decompose the interference noise of each line and has a good denoising ability. After VMD decomposition was completed, Hj was obtained by the line selection method discussed above, as shown in Figure 11.

Cable Line Fault
The A-phase grounding fault with a grounding resistance of 500 Ω at the cable line L4, 3 km away from the bus, is presented here as an example, where the initial phase angle of the fault is 90°. Gaussian white noise with an SNR of 30 dB was injected into the collected zero sequence current of each outgoing line as shown in Figure 12.  Figure 11 is the NPFC H j obtained here after smoothing and filtering the NPFC of each line. The Pearson correlation coefficient of each outgoing line H j was calculated to form the correlation coefficient matrix M: The comprehensive correlation coefficient P = [−0.9982, 0.3327, 0.3327, 0.3318] was calculated after normalization. S is much larger than P set in this case. As a result, the line corresponding to the minimum value, line L 1 can be judged as a faulty line. The line selection result is accurate.

Cable Line Fault
The A-phase grounding fault with a grounding resistance of 500 Ω at the cable line L 4 , 3 km away from the bus, is presented here as an example, where the initial phase angle of the fault is 90 • .
Gaussian white noise with an SNR of 30 dB was injected into the collected zero sequence current of each outgoing line as shown in Figure 12. minimum value, line L1 can be judged as a faulty line. The line selection result is accurate.

Cable Line Fault
The A-phase grounding fault with a grounding resistance of 500 Ω at the cable line L4, 3 km away from the bus, is presented here as an example, where the initial phase angle of the fault is 90°. Gaussian white noise with an SNR of 30 dB was injected into the collected zero sequence current of each outgoing line as shown in Figure 12.  The NPFC H j after the smoothing and filtering of each outgoing line is shown in Figure 13. The NPFC Hj after the smoothing and filtering of each outgoing line is shown in Figure 13.

Bus Fault
A SPG fault with a grounding resistance of 100 Ω was next imposed at the bus bar, where the initial phase angle of the fault is 0°. Gaussian white noise with an SNR of 30 dB was injected into the collected zero-sequence current of each outgoing line as shown in Figure 14. The comprehensive correlation coefficient P = [0.4486, 0.4512, 0.4516, −0.6473] was obtained after the Pearson correlation coefficient was calculated. S was greater than P set in this case, so the line corresponding to the minimum value (L 4 ) can be judged as faulty. The line selection result is accurate.

Bus Fault
A SPG fault with a grounding resistance of 100 Ω was next imposed at the bus bar, where the initial phase angle of the fault is 0 • . Gaussian white noise with an SNR of 30 dB was injected into the collected zero-sequence current of each outgoing line as shown in Figure 14.
corresponding to the minimum value (L4) can be judged as faulty. The line selection result is accurate.

Bus Fault
A SPG fault with a grounding resistance of 100 Ω was next imposed at the bus bar, where the initial phase angle of the fault is 0°. Gaussian white noise with an SNR of 30 dB was injected into the collected zero-sequence current of each outgoing line as shown in Figure 14. As shown in Figure 14, when the ground fault occurred on the bus, the initial polarities of the zero-sequence current of each line were the same. The NPFC Hj of each smoothed and filtered outgoing line is shown in Figure 15. As shown in Figure 14, when the ground fault occurred on the bus, the initial polarities of the zero-sequence current of each line were the same. The NPFC H j of each smoothed and filtered outgoing line is shown in Figure 15. As shown in Figure 15, the polarities of all outgoing lines Hj are the same. The comprehensive correlation coefficient P = [0.9937 0.9974 0.9975 0.9974] can be calculated accordingly. S is less than set P here, so a bus fault is identified. The line selection result is accurate.

Applicability Analysis
The fault resistance, initial phase angle of the fault, intensity of the injected noise, and fault location were changed in succession for the purposes of further analysis. The proposed algorithm was used to obtain the the absolute value S of the difference between the comprehensive correlation coefficient P of the NPFCs of each line and the maximum and minimum values. The results, as presented below, confirm the applicability of the faulty line selection method based on VMD-FFT NPFC correlation analysis.
The SPG fault of the overhead line L1 is presented here as an example. At this time, the initial phase angle of the fault is 90°, the fault is 5 km away from the bus, the SNR is 30 dB, and the fault resistance increases from 1 Ω to 1000 Ω. The changes in the NPFCs P and S of each feeder are shown in Figure 16. As shown in Figure 15, the polarities of all outgoing lines H j are the same. The comprehensive correlation coefficient P = [0.9937 0.9974 0.9975 0.9974] can be calculated accordingly. S is less than P set here, so a bus fault is identified. The line selection result is accurate.

Applicability Analysis
The fault resistance, initial phase angle of the fault, intensity of the injected noise, and fault location were changed in succession for the purposes of further analysis. The proposed algorithm was used to obtain the the absolute value S of the difference between the comprehensive correlation coefficient P of the NPFCs of each line and the maximum and minimum values. The results, as presented below, confirm the applicability of the faulty line selection method based on VMD-FFT NPFC correlation analysis.
The SPG fault of the overhead line L 1 is presented here as an example. At this time, the initial phase angle of the fault is 90 • , the fault is 5 km away from the bus, the SNR is 30 dB, and the fault Energies 2020, 13, 4724 14 of 18 resistance increases from 1 Ω to 1000 Ω. The changes in the NPFCs P and S of each feeder are shown in Figure 16.
was used to obtain the the absolute value S of the difference between the comprehensive correlation coefficient P of the NPFCs of each line and the maximum and minimum values. The results, as presented below, confirm the applicability of the faulty line selection method based on VMD-FFT NPFC correlation analysis.
The SPG fault of the overhead line L1 is presented here as an example. At this time, the initial phase angle of the fault is 90°, the fault is 5 km away from the bus, the SNR is 30 dB, and the fault resistance increases from 1 Ω to 1000 Ω. The changes in the NPFCs P and S of each feeder are shown in Figure 16.  As the resistance increases, the comprehensive correlation coefficients of lines L 2, L 3 , and L 4 become basically approximated. Although the comprehensive correlation coefficient of line L 1 increases to some extent, it is still smaller than the other three-circuit lines. The S value indicates that line L 1 is the faulty line. These results show that the proposed algorithm can be used to accurately select the faulty line as the grounding resistance changes.

Changing Initial Fault Phase Angle
When a SPG fault occurs on the overhead line L 1 , the fault resistance is 5 Ω, and the fault is 5 km away from the bus. As the initial phase angle of the fault changes from 0 • to 90 • , the changes in P and S of the NPFCs of each outgoing line are as shown in Figure 17. As the resistance increases, the comprehensive correlation coefficients of lines L2, L3, and L4 become basically approximated. Although the comprehensive correlation coefficient of line L1 increases to some extent, it is still smaller than the other three-circuit lines. The S value indicates that line L1 is the faulty line. These results show that the proposed algorithm can be used to accurately select the faulty line as the grounding resistance changes.

Changing Initial Fault Phase Angle
When a SPG fault occurs on the overhead line L1, the fault resistance is 5 Ω, and the fault is 5 km away from the bus. As the initial phase angle of the fault changes from 0° to 90°, the changes in P and S of the NPFCs of each outgoing line are as shown in Figure 17. As shown in Figure 17, when the initial phase angle of the fault changes, the faulty line can be accurately selected based on the P and S of each feeder.

Changing Fault Location
The SPG fault of the cable line L4 was next given a grounding resistance of 100 Ω, initial phase angle of the fault of 0°, and faults 2 km, 3 km, 4 km, and 5 km from the bus, respectively. The changes in P and S of the NPFCs of all outgoing lines are shown in Figure 18. As shown in Figure 17, when the initial phase angle of the fault changes, the faulty line can be accurately selected based on the P and S of each feeder.

Changing Fault Location
The SPG fault of the cable line L 4 was next given a grounding resistance of 100 Ω, initial phase angle of the fault of 0 • , and faults 2 km, 3 km, 4 km, and 5 km from the bus, respectively. The changes in P and S of the NPFCs of all outgoing lines are shown in Figure 18. As shown in Figure 17, when the initial phase angle of the fault changes, the faulty line can be accurately selected based on the P and S of each feeder.

Changing Fault Location
The SPG fault of the cable line L4 was next given a grounding resistance of 100 Ω, initial phase angle of the fault of 0°, and faults 2 km, 3 km, 4 km, and 5 km from the bus, respectively. The changes in P and S of the NPFCs of all outgoing lines are shown in Figure 18.  As shown in Figure 18, when the fault location changes, the faulty line can still be quickly and accurately selected based on the P and S values of each feeder.

Changing Injected Noise Intensity
When a SPG fault occurs on the cable hybrid line L 3 , the fault resistance is 100 Ω, the initial phase angle of the fault is 90 • , the fault is 3 km, and the SNRs are −10 dB, 0 dB, 10 dB, 20 dB, and 30 dB, respectively, the changes in P and S of the NPFCs of each feeder are as shown in Figure 19.
Energies 2020, 13, x FOR PEER REVIEW 16 of 19 As shown in Figure 18, when the fault location changes, the faulty line can still be quickly and accurately selected based on the P and S values of each feeder.

Changing Injected Noise Intensity
When a SPG fault occurs on the cable hybrid line L3, the fault resistance is 100 Ω, the initial phase angle of the fault is 90°, the fault is 3 km, and the SNRs are −10 dB, 0 dB, 10 dB, 20 dB, and 30 dB, respectively, the changes in P and S of the NPFCs of each feeder are as shown in Figure 19. As shown in Figure 19, when the intensity of the noise changes, the P and S of each feeder change very slightly. This suggests that the proposed method can be used to accurately select faulty lines with high noise tolerance.

Comparison Analysis
In order to verify the superiority and effectiveness of the proposed method, another two commonly used fault detection methods are compared, the one is based on the EMD algorithm and As shown in Figure 19, when the intensity of the noise changes, the P and S of each feeder change very slightly. This suggests that the proposed method can be used to accurately select faulty lines with high noise tolerance.

Comparison Analysis
In order to verify the superiority and effectiveness of the proposed method, another two commonly used fault detection methods are compared, the one is based on the EMD algorithm and correlation analysis [22], the second one is based on the Prony algorithm and correlation analysis [16].
When the SPG fault occurs on the hybrid line L 2 , and the fault initial phase angle is 90 • , the fault resistance is 1 Ω, 500 Ω, and 1000 Ω, respectively. The final fault detection results and the comparisons are listed in Table 2. As shown in Table 2, when the fault resistance is equal to 1 Ω, the three methods can detect the fault line accurately. When the fault resistance is equal to 500 Ω, both the proposed method and EMD-based method achieve accurate fault detection results, while the Prony-based method misjudges the faulty line. When the resistance increases to 1000 Ω, only the proposed method can accurately detect the faulty line, while the other two give the wrong results. Therefore, when a SPG fault occurs with a large grounding resistance, and the fault detection results based on the EMD algorithm and the Prony algorithm are no longer reliable, the superiority and effectiveness of the proposed method can be verified.
VMD and FFT were combined in this study to accurately separate power frequency components and noise signals in a TZSC. The proposed algorithm was found to be more accurate than EMD-based method and Prony-based method.

Conclusions
By analyzing the composition and characteristics of the TZSC of a resonant grounding system, a novel faulty line selection method was established in this study based on the correlation analysis of NPFCs via VMD-FFT. The main conclusions of this work can be summarized as follows.

1.
VMD and FFT are combined to prevent any inaccurate separation of power frequency components caused by modal aliasing and boundary effects, and to provide a sound basis for subsequent line selection.

2.
The moving average filter smoothes and filters the NPFCs to reveal the overall trend of the NPFCs. This highlights the correlations between faulty lines and non-faulty lines. 3.
The proposed line selection algorithm is relatively unaffected by factors such as fault resistance and fault initial phase angle and shows strong anti-noise characteristics.