An improved natural frequency based transmission line fault location method with full utilization of frequency spectrum information

Funding information National Natural Science Foundation of China, Grant/Award Number: 51807119 Abstract An improved natural frequency based method is proposed to accurately locate faults in transmission lines. The method only requires three phase instantaneous current measurements at the local terminal of the line and is compatible with protective relays using highspeed tripping techniques. First, the theoretical peak frequencies along the entire frequency spectrum of the measured currents are derived, for different fault types, locations and resistances. Next, the fault location is determined by finding the minimum average distance between the theoretical and the measured peak frequencies. Compared to the existing natural frequency based fault location methods, the proposed method solves the mode mixing problem during single phase to ground faults and systematically considers the effect of fault resistances. In addition, the proposed method fully utilizes the frequency spectrum information instead of only the dominant frequency, and therefore overcomes the difficulty to extract the dominant frequency. Numerical experiments have shown that, compared to the existing method, the proposed method presents much higher fault location accuracy during single phase to ground faults, and slightly higher (comparable) fault location accuracy for other types of faults, regardless of fault locations and fault resistances.


INTRODUCTION
Accurate fault location method is beneficial for reducing the power outage time, as well as, operational costs of power systems. Fault location methods of transmission lines are studied in many literatures. The most widely adopted fault location methods are the fundamental frequency phasor based methods, which utilize fundamental frequency phasors to calculate the impedance and afterwards the distance between the fault and one terminal of the line [1][2][3][4][5][6]. The main limitation of the fundamental frequency phasor based methods is the dependency on the accurate extraction of fundamental frequency phasors, which could be challenging especially during system transients. Therefore, for transmission lines equipped with protective relays using high-speed tripping techniques, the available data window to extract the phasors during faults could be very short (e.g. fraction of a cycle), and these fault location methods may present compromised fault location results. To overcome the above limitation, researchers proposed advanced fault location methods that are compatible with the short data window, including time domain model based methods, traveling wave based methods and natural frequency based methods. Time domain model based methods determine the location of the fault with the help of the accurate time domain transmission line model before and during the fault. Specifically, in "solving equation" methods, the fault location is obtained by solving analytical equations that describe the physical laws between the instantaneous measurements and the fault location [7][8][9]. However, to ensure analytical expressions, these transmission line models are usually limited to lumped parameter models or multi-section π models, which may increase the fault location errors. To further utilize more accurate transmission line models, the voltage methods calculate the voltage distributions along the transmission line at either terminal and then find the intersection point of the two voltage distributions to obtain the fault location. These methods do not need to formulate the analytical equations and therefore more accurate distributed parameter transmission line models can be applied such as the Bergeron model (with assumptions of lumped resistances) [10,11] or even models described via the matrix form partial differential equations (fully distributed parameter models) [12]. Nevertheless, these methods usually require dual ended voltage and current measurements, and synchronization among the measurements are typically required.
Traveling wave based methods determine the fault location by detecting the arrival time of wavefronts at line terminals after the occurrence of the fault, including single ended and dual ended algorithms. Single ended methods can be further classified into types A, C, E and F methods (according to the generation of the traveling waves). The effectiveness of single ended algorithms [18][19][20] relies on the successful detection of multiple reflections of traveling waves between line terminal and fault location. Dual ended algorithms [13][14][15][16][17] observe the arrival time of first wavefronts at both line terminals and then obtain the fault location by the time difference, which can be further classified into types B and D methods (according to synchronization techniques). However, the reliable detection of wavefronts is challenging with small fault inception angles or high fault resistances, which increases difficulty to accurately capture the fault location.
To avoid reliability issues in detecting wavefronts of traveling waves, researchers proposed the natural frequency based methods. Instead of accurately detecting the arrival time of wavefronts of traveling waves, these methods adopt the frequency characteristics of the voltages/current measurements to locate faults. The preliminary relationship between the fault location and the natural frequency was studied by Swift in 1979 [21], where only special system conditions are considered: the system impedance is assumed to either infinite or zero. In recent years, researchers largely improve the performance of the natural frequency based fault location method by considering more general system boundary conditions. These natural frequency based fault location approaches have been applied in both HVAC [24,25] and HVDC [22,23] systems. The main idea of these methods is to first extract the dominant frequency from the single ended instantaneous current or voltage measurements of a certain transmission line mode, and afterwards use the wave speed and the reflection coefficient of a certain transmission line mode at the local terminal to locate faults. In fact, for three phase AC systems, the natural frequency based fault location methods work well for phase to phase faults, double phase to ground faults and three phase faults. Nevertheless, when a single phase to ground fault occurs, there exists the "mode mixing phenomenon" [23,24], where the frequency spectrum of a certain mode is mixed with that of other modes. Therefore, since the derivations of the existing methods are usually within one transmission line mode, the mode mixing phenomenon may cause difficulty to extract the dominant frequency and also compromised fault location results. In addition, the existing methods typically assume zero fault resistance, which may further increase fault location errors especially during high resistance faults.
In this paper, a novel natural frequency based method is proposed. It only needs three phase instantaneous current measurements at the local terminal of the line, and no communication channel or time synchronization is required. Unlike the existing method which only uses the dominant frequency, the proposed method fully exploits the frequency spectrum information and utilizes a series of peak frequencies of the frequency spectrum. First, the frequency spectrum of the measured currents for different fault types, locations and resistances is derived in detail. Afterwards, the relationship between the fault location and the peak frequencies is studied, with full consideration of the mode mixing phenomenon and the fault resistance. Finally, the fault location is determined by solving an optimization problem that finds the best match between the theoretical and the measured peak frequencies. The original contributions of this paper are highlighted as follows: • The method only needs single ended current measurements and a relatively short time window; compared to the traveling wave based methods, the method does not require reliable wavefront detection and works with zero fault inception angle. • Instead of only using the dominant frequency, the proposed method strictly and systematically derives all the peak frequencies along the frequency spectrum of the measured currents for different fault types, locations and resistances; the proposed method fully utilizes the information embedded inside the frequency spectrum, and therefore avoids the difficulty to extract the dominant frequency. • During single phase to ground faults, the proposed method systematically considers the coupling among different line modes as well as the fault resistance compared to the existing method, and therefore fully solves the mode mixing problem. • During other types of faults (phase to phase, double phase to ground, and three phase faults), the proposed method considers the fault resistance compared to the existing method; mathematical proofs illustrate that the existing method is a special case of the proposed method with zero fault resistance.
• Numerical experiments prove that the proposed method presents much higher fault location accuracy for single phase to ground faults, and slightly higher (comparable) fault location accuracy for other types of faults, compared to existing natural frequency based method.
The rest of the paper is organized as follows. Section 2 introduces the fundamental principles and reviews in detail the existing natural frequency based fault location method as well as the limitations. Sections 3 and 4 provide the strictly mathematical derivation to obtain the peak frequencies of the frequency spectrum for different types of faults. Section 5 introduces proposed the fault location algorithm. Section 6 compares the proposed method and the existing method via numerical experiments to show the advantages of the proposed method. Section 7 discusses the effects of different factors on the fault location accuracy. Section 8 draws a conclusion.

Review of the healthy three phase lossless line model
The partial differential equations representing a healthy three phase lossless transmission line model are [26], where matrices L and C are the inductance and capacitance matrices per unit length, u m (x, t ) and i m (x, t )(m = a, b, c) are per phase voltages and currents on the transmission line at location x and time t.
To solve Equation (1), the modal decomposition is usually adopted to transform phase components abc into mode components 0αβ (here the Clarke transformation is taken as an example), where the Clarke transformation matrix T is defined as T = [1∕3, 1∕3, 1∕3; 2∕3, −1∕3, −1∕3; 0, 1∕ (1) and (4), where and diag([⋅]) is a diagonal matrix with the diagonal vector [⋅]. From Equation (5), the mode voltages and currents are fully decoupled for healthy lines. The solution to Equation (5) is [26], where j = 0, α, β represents different modes. For each mode j, z j = √ L j ∕C j is the surge impedance, f j + and f j − are the forward and backward waves, with speed v j = 1∕ √ L j C j .

Review of the existing natural frequency based fault location method
The frequency spectrum of the fault voltage/current measurements is called the natural frequency. After modal decomposition, the natural frequency corresponding to each mode can be represented by roots of Equation (9) [22], where P (s) = e −sT , Γ 1 (s) and Γ f (s) are reflection coefficients at the line end and fault location in Laplace domain respectively, which can be calculated by the parameters of the system. The spectrum of the signal can be extracted using MUSIC algorithm [22]. There are peaks in the frequency spectrum, and the frequency corresponding to the first peak is selected as the dominant frequency to calculate the fault location. The dominant frequency has the largest amplitude other than the fundamental frequency, and can be identified easily. The existing method firstly extracts the frequency spectrum of the terminal current of a certain mode, and then finds the dominant frequency in the spectrum. With the calculated wave speed and the reflection coefficient at the line end of a certain mode, the fault location is determined as, where f d is the dominant frequency of the current of the certain mode measured at the line end after the occurrence of the fault, v is the wave speed of the certain mode (or average wave speed of several modes if mode mixing phenomenon occurs), and θ 1 is the angle of reflection coefficient of the certain mode at the line terminal. Note that the fault resistances are assumed to be zero during the derivation. The limitations of the existing method are described in detail as follows. First, as mentioned in refs. [23,24], when a single phase to ground fault occurs on the transmission line, there exists the mode mixing phenomenon, which increases the difficulty of accurate fault location. For example, for a three phase transmission line with Clarke transformation for modal decomposition, mode α and mode 0 are coupled (mixed) together during single phase to ground faults, that is, the spectrum extracted by the mode α current will consist of the frequency information of both mode α and mode 0, vice versa. Consequently, the spectrum of the mode current is rather complicated, resulting in challenge to find the correct dominant frequency. Moreover, even with the correct dominant frequency, it is also challenging to determine the accurate wave speed for the mixed modes. In refs. [23,24], the existing methods approximate the two mixed modes as one mode, with the wave speed as the average speed 1. The existing method only utilizes the dominant frequency; other information of the frequency spectrum is not fully utilized; 2. Due to mode mixing phenomenon especially during single phase to ground faults, accurate extraction of dominant frequency and the wave speed of the mixed mode is challenging; 3. Derivations are within one transmission line mode, and the fault resistances are assumed to be zero.

PROPOSED NATURAL FREQUENCY BASED FAULT LOCATION METHOD FOR SINGLE PHASE TO GROUND FAULTS
To overcome the limitations of the existing natural frequency based method, a novel natural frequency based fault location method is proposed in this section. The proposed method fully utilizes the information of the frequency spectrum instead of only one dominant frequency. The mathematical relationship among the frequency spectrum, the fault location and the fault resistance is rigorously derived for different types of faults. Afterwards, the fault location is determined with full utilization of the information in the frequency spectrum. Specifically, the mode mixing problem for single phase to ground faults as well as the effect of fault resistances are systematically considered using the proposed fault location methodology.

Frequency spectrum during single phase to ground faults
Here phase A to ground faults are taken as examples to show the frequency spectrum as functions of fault location and fault resistances. Also, the relationship between the existing method and the proposed method during single phase to ground faults is explained. The example power system of interest is shown in Figure 1. Some definitions of variables are as follows. The entire transmission line during faults is divided into two healthy transmission line sections as Line 1 (left side of the fault) and Line 2 (right side of the fault). For mode j (j = 0, α, β) or phase j (j = a, b, c), the voltages and currents on Line 1 and Line 2 are defined as u j (1) (x 1 , t ), i j (1) (x 1 , t ), u j (2) (x 2 , t ) and i j (2) (x 2 , t ), respectively, where x 1 and x 2 denote the distance to the left and the right terminal, respectively; the voltage of the equivalent sources at the left and the right side of the line is u j s1 (t ) and u j s2 (t ), respectively; the voltage at the location of the fault and the current flowing into the fault are defined as u j f (t ) and i j f (t ). The voltage and current at the corresponding voltage and current vectors are defined as where w can be replaced by u (1) , i (1) , u (2) , i (2) , u f , i f , u s1 (t ) and u s2 (t ). The fault is at distance d from the left terminal. The length of the entire transmission line is l .

Boundary conditions at the fault location and line terminal
For a phase A to ground fault, at the location of the fault, where with fault resistance R f . For the mode network, from Equation (4), at the location of the fault, where The boundary conditions at the fault location can be obtained by observing that Line 1 and Line 2 share the same node. First, the voltage at the right terminal of Line 1 and the left terminal of Line 2 are the same (Equation (15)). Second, the Kirchhoff 's Current Laws at the fault location should be obeyed (Equation (16)). Therefore, The boundary conditions at the line ends can be obtained by observing Kirchhoff's Voltage Laws from the equivalent sources to the transmission line terminals. where and are the left side source resistance and inductance matrices for the mode network, R abc s1 and L abc s1 are the corresponding matrices for the phase networks.
The definitions are similar for the right side source (with superscript 's2').

Solving the peak frequencies of the frequency spectrum
In this section, since the line with fault consists of two healthy lines (Line 1 and Line 2) as well as, the fault, the frequency spectrum can be obtained using the solution in Section 2.1 and boundary conditions in Section 3.2. Here the α mode and 0 mode are taken as examples. Through the derivation, the wave speeds of each mode are defined in Equation (8).
Substitute Equation (8) into Equation (15) and take the first two rows, Substitute Equations (8) and (15) into Equation (16) and take the first two rows, Substitute Equation (8) into Equations (17) and (18) and take the first two rows, To derive the frequency spectrum, we transform Equations (22)-(25) into Laplace domain in Equations (26)-(29), formulated according to the solution of these 8 unknowns. Due to space limitations, here we only present the α mode current at the left terminal i (1) (0, s) (other voltages and currents can be similarly obtained). After converting Equation (8) into Laplace domain and substituting the solutions, (32) In Equation (30), one can observe that A(s) includes the source voltages U s1 (s) and U s2 (s), denoting that A(s) forms the frequency peak at the fundamental frequency (source frequency: 50 or 60 Hz). Therefore, if we let s = j = j 2 f , where ω and f are the angular frequency and the frequency, the peaks of the frequency spectrum (peak frequencies) except the peak of fundamental frequency can be represented as the local minimum of B( j 2 f ) or local maximum of where the function B(⋅) is defined in the expression of the mode current. It can be observed that there could be many local maxima of 1∕B( j 2 f ) (corresponds to the peak frequencies of the mode current). These peak frequencies are functions of the pre-defined coefficients, and therefore are functions of the fault location d and the fault resistance R f .

Relationship between the existing and the proposed method during single phase to ground faults
For the proposed method, from Equations (30) and (34), it can be seen that the peak frequencies through the frequency spectrum are related to both mode α and mode 0 networks (for example in Equation (30), variables a, p(s), m(s), c 1 (s) and c 2 (s) are related to mode α; variables b, q(s), n(s), d 1 (s) and d 2 (s) are related to mode 0), proving the correctness of the mode mixing phenomenon (mode α and mode 0 are mixed together and cannot be simply decoupled). One can also observe that the frequency spectrum is a function of the fault resistance (variables a and b). On the other hand, for the existing method, from Equation (10), the dominant frequency (corresponding to one of the peak frequencies in the proposed method) is derived accordingly to network of only one certain mode and the fault resistance is not considered. Even if the wave speed in Equation (10) is approximated as the average speed in mode α and mode 0 for single phase to ground fault [23,24], the existing method will still present errors. Therefore, during single phase to ground faults, the proposed method systematically considers the mode mixing phenomenon and the fault resistance compared to the existing method.

PROPOSED NATURAL FREQUENCY BASED FAULT LOCATION METHOD FOR OTHER TYPES OF FAULTS
Next, the frequency spectrum is derived during other types of faults (including phase to phase, double phase to ground and three phase faults). For these faults, at least one certain mode is not coupled with the rest of the modes, and therefore the derivation of the frequency spectrum is much simplified. Also, the relationship between the existing method and the proposed method during these faults is explained. The example power system of interest is shown in Figure 1. Note that the definitions of variables are consistent with those in Section 3.

Phase to phase faults
Here phase B to phase C faults are taken as examples. At the location of the fault, the expressions are the same as Equations (11) and (13) except that Y abc f and Y 0 f are updated as and The boundary conditions are the same as Equations (15)- (18). From the diagonal structure of Y 0 f , it can be seen that the mode 0αβ are fully decoupled. Here take mode β as an example. Similarly, substitute Equation (8) into Equation (15), substitute Equations (8) and (15) into Equation (16), substitute Equation (8) into Equations (17) and (18), and take the third row (corresponding to mode β). Afterwards, transform them into Laplace domain,  Table 2.

Double phase to ground faults
Here phase BC to ground faults are taken as examples. At the location of the fault, the expressions are the same as Equations (11) and (13) except that Y abc f and Y 0 f are updated as and The boundary conditions are the same as Equations (15)- (18). Similarly, one can observe from the structure of Y 0 f that mode β is not coupled with the rest of the modes and therefore mode β is selected. Similarly, the expression of i (1) (0, s) can be derived and the result is the same as Equation (41) except that the definition of a is updated as a = R f ∕z . Finally, the peak frequencies can be found through Equation (34).

Three phase faults
At the location of the fault, the expressions are the same as Equations (11) and (13) except that Y abc f and Y 0 f are updated as and The boundary conditions are the same as Equations (15)- (18). Similarly, the expression of i (1) (0, s) and the peak frequencies can be derived. The results are exactly the same as those of double phase to ground faults in Section 4.2.

Relationship between the existing and the proposed method during other types of faults
For the proposed method, when dealing with the faults other than single phase to ground faults, the peak frequencies through the frequency spectrum can be found from Equations (41) and (34), and are only related to mode β, proving that there is no mode mixing phenomenon for these types of faults. Nevertheless, one can observe that the peak frequencies are still functions of the fault resistance (variable a). On the other hand, for the existing method, from Equation (10), the expression of the dominant frequency (corresponding to one of the peak frequencies in the proposed method) does not consider the fault resistance. Therefore, during faults other than single phase to ground faults, the proposed method considers the fault resistance compared to the existing method. In fact, the existing method is a special case of the proposed method when the fault resistance is zero. Next, a brief proof is provided.
Proof. If R f = 0, from the definition of the variable a, a = 0. Therefore, Equation (41) can be simplified as, Therefore, the peak frequencies can be found when, Substitute s = j = j 2 f into Equation (49), we have, where the reflection coefficient 1 is the angle of Γ s , and The phase angles at both sides of Equation (50) should be equal or with differences 2kπ, Therefore, the first peak frequency corresponds to k = 0. We can conclude that the expression is equivalent to Equation (10). Q.E.D.

PROPOSED FAULT LOCATION ALGORITHM
The peak frequencies of the measurements during faults are extracted. Here the three phase instantaneous current measurements at the local terminal are utilized. The three phase currents are first transformed into mode currents. Afterwards, the spectrum of one certain mode current can be calculated by MUSIC algorithm [22]. Within a certain frequency range, in the spectrum of the current, there exist many peak frequencies denoted as f i measured (i = 1, 2, …, n) other than the fundamental frequency of the system (50 or 60 Hz).
Next, the extracted peak frequencies of the measurements are compared to the theoretical peak frequencies to obtain the fault location. These theoretical peak frequencies can be easily extracted offline (even before the installation of the fault location device) by first substituting different frequencies into the objective function in Equation (34) and afterwards applying a numerical peak-searching algorithm. The derivations in Sections 3 and 4 demonstrate that the theoretical peak frequencies f i theoretical (i = 1, 2, …, n) are functions of the fault location d and the fault resistance R f . The validity of the proposed method is based on the following fact: The correct fault location and fault resistance should correspond to the best match between f i measured and f i theoretical . Therefore, this is equivalently to solve an optimization problem, where is the objective function that describes the average difference between the theoretical and measured peak frequencies.
In practice, the number of theoretical and measured peak frequencies may not be exactly the same within a certain frequency band (due to measurement errors, parameter errors etc.). Assume that the number of theoretical peak frequencies is m while the number of measured peak frequencies is still n. To improve the robustness of the proposed fault location method, the proposed method first finds the closest theoretical and measured peak frequency pairs and calculates their absolute frequency differences |f j theoretical −f k measured |. Afterwards, these two frequencies are removed from the two peak frequency datasets, and the same procedure is executed. In order to further improve the robustness of this fault location algorithm and reject bad data/outliers, the procedure is executed until a certain percentage (e.g. 80%) of either the theoretical or measured peak frequency dataset is covered. The procedure of calculating g(d, R f ) with given d and R f is summarized as follows.
There are many optimization tools to solve Equation (54). Since here the theoretical frequencies with different d and R f can be calculated offline (even prior to the fault), the main calculation burden is to find the minimum value of the matrix A as well as some summations and multiplications after extracting the measured peak frequencies, which are not quite high. Therefore, a very effective and straightforward way is to draw a mesh with different d and R f , calculate all the values of g(d, R f ), and directly find the minimum. Here we utilize a two-step algorithm to reduce the computational burden.

NUMERICAL EXPERIMENTS
The proposed method is validated via a 50 Hz system with a 500 kV, 200 km three phase ideally transposed transmission line built in PSCAD/EMTDC. The example test system is the same as  The theoretical peak frequencies are calculated offline prior to the fault. Here we consider the frequency spectrum with the maximum frequency of 10 kHz. The corresponding variables are selected as follows: Δl f,step1 = 0.5 km, ΔR f,step1 = 5 ohm, Δl f,step2 = 0.1 km, ΔR f,step2 = 1 ohm. Also, for single phase to ground faults, the theoretical peak frequencies can be calculated by Equations (30) and (34), with R f,max = 600 ohm; for other types of faults, the theoretical peak frequencies can be calculated by Equations (41) and (34), with R f,max = 50 ohm (these selection of maximum fault resistances can cover extreme situations). For the computational burden of the algorithm, here we take the single phase to ground faults (with more computational burden since R f,max is larger) as examples. The maximum calculation time without parallel computing is less than 19 seconds using Matlab R2018b on a personal computer with Intel i7-7700 CPU (with 48,100 values of g(d, R f ); on average for each g(d, R f ) the calculation time is less than 400 μs). This calculation time for fault location is acceptable in practice. In addition, the calculation time can be further minimized if the values of g(d, R f ) with different d and R f are calculated using parallel computing techniques, since the values of g(d, R f ) with different d and R f are independent and can be computed in parallel.
Next, the performances of the proposed method are compared to those of the existing natural frequency based fault location method (as described in Section 2.2 in detail) via different test cases. Specifically, for single phase to ground faults, low and high resistance faults (0.01 ohm, 1 ohm, 10 ohm; 100 ohm, 300 ohm, 500 ohm) are studied to demonstrate the effect of mode mixing phenomenon as well as the fault resistance on both methods. For other types of faults, a specific fault resistance of 10 ohm is provided to demonstrate the effect of fault resistance on both methods.

6.1
Low resistance single phase to ground faults A 0.01 ohm A-G fault at 50 km from the local terminal is taken as an example. For this fault event, the three phase current instantaneous measurements at the local terminal are shown in Figure 3 (the current measurements are not shown for other fault events due to space limitations). The values of g(d, R f ) in step 1 and step 2 of the proposed algorithm are shown in Figure 4 Figure 5 depicts the comparison between the theoretical and the measured frequency spectrums in detail. The theoretical frequency spectrum is with the optimal fault location d = 49.6 km and fault resistance R f = 1 ohm. One can observe that the peak frequencies of the theoretical frequency spectrum are very consistent with those of the measured frequency spectrum, proving the effectiveness of the proposed calculation method of the theoretical peak frequencies as well as the proposed fault location method. Next, the performance of the existing method is also studied. From Figure 5, there are many peak frequencies in the measured frequency spectrum and it is not trivial to extract the dominant natural frequency. Here the first peak frequency (other than the fundamental frequency of 50 Hz) of the measured frequency spectrum is not the dominant frequency: The first peak frequency is 542.9 Hz, resulting in 204.34 km fault location according to Equation (10) (with wave speed approximated as the average speed of mode α and mode 0 [23,24]). In fact, if we determine the actual dominant frequency by finding the peak frequency that results in the most accurate fault location, the dominant frequency is 1822.9 Hz, as marked in Figure 5. It can be observed that it is very hard to determine this dominant frequency. This fact is caused by the mode mixing phenomenon. Second, even we assume that the dominant frequency can somehow be found, the fault location result of the existing method according to Equation (10) is 48.38 km (the absolute fault location error is 0.81%). Therefore, the fault location accuracy of the proposed method is improved compared to the existing method.
The effectiveness of proposed method is further validated using a group of low fault resistance (0.01 ohm, 1 ohm, 10 ohm) A-G faults at different fault locations through the line (21 fault locations, 2 km, 198 km and every 10 km from ref. [10], 190 km from local terminal). The results are shown in Figure 6(a) and Table 3. Here we assume that the dominant frequency can somehow be accurately extracted for the existing method, to show the best performance of the existing method for comparison. The results demonstrate that compared to the existing method, fault location accuracy of the proposed method is also significantly improved.

High resistance single phase to ground faults
A 500 ohm A-G fault at 90 km from the local terminal is taken as an example. The values of g(d, R f ) in step 1 and step 2 of the proposed algorithm are shown in Figure 7(a),(b). The fault location and fault resistance that minimize g(d, R f ) are 90.0 km and 470 ohm, respectively, which are close to actual fault location and fault resistance (the absolute fault location error is 0%). Figure 8 depicts the comparison between the theoretical and the measured frequency spectrums in detail. Similarly, one can observe that the peak frequencies of the theoretical frequency spectrum are very consistent with those of the measured frequency spectrum.
Next, the performance of the existing method is also studied. Similarly, there are many peak frequencies in the measured frequency spectrum and it is not trivial to extract the dominant natural frequency (similarly, if we determine the actual dominant frequency by finding the peak frequency that results in the most accurate fault location, the dominant frequency is 1009.8 Hz, as marked in Figure 8). Even we assume that the dominant frequency can somehow be found, the fault location result of the existing method according to Equation (10) Table 4. Similarly, here we assume that the dominant frequency can somehow be accurately extracted for the existing

Other types of faults
Next, phase to phase faults, double phase to ground faults and three phase faults are considered. Due to space limitations, here we only take a 10 ohm phase B to phase C fault occurs at 100 km as an example (for other types of faults the results are similar). The values of g(d, R f ) in step 1 and step 2 of the proposed algorithm are shown in Figure 9(a),(b). The fault location and fault resistance that minimize g(d, R f ) are 100.1 km and 9 ohm, respectively, which are close to actual fault location and fault resistance (the absolute fault location error is 0.05%). Figure 10 depicts the comparison between the theoretical and the measured frequency spectrums in detail. Similarly, one can observe that the peak frequencies of the theoretical frequency spectrum are very consistent with those of the measured frequency spectrum. The performance of the existing method is also studied. In this case, it is straightforward to determine the dominant natural frequency: it is the first peak frequency 1123.1 Hz, and the fault location result of the existing method according to Equation (10) is 100.16 km (the absolute fault location error is 0.08%). The existing method has similar fault location accuracy compared to the proposed method. This is because for other types of faults (other than single phase to ground faults) there does not exist the mode mixing phenomenon and therefore the fault location result is rather accurate for the existing method.
The effectiveness of proposed method is further validated using phase B to phase C faults, phase BC to ground faults, and three phase faults with 10 ohm fault resistance and different fault locations through the line (21 fault locations, 2 km, 198 km and every 10 km from ref. [10], 190 km from local terminal). The results are shown in Figure 11 and Table 5. The results demonstrate that the proposed method presents slightly higher (comparable) fault location accuracy than the existing method. This is because the proposed method considers the effect of fault resistance. Nevertheless, for these types of faults, both the proposed method and the existing method present accurate fault location results, and the improvements while considering the fault resistance using the proposed method are not significant.

DISCUSSION
In this section, the effectiveness of proposed method is further discussed with different factors. All the conditions of example test system are the same as those in Section 6. For all the following studies (if not specially mentioned), 0.01 ohm A-G faults with different fault locations (21 fault locations, 2 km, 198 km and every 10 km from [10,190] km from the local terminal) are taken as examples.

Effect of loading conditions
The effect of loading conditions (phase angle of source 1 leading that of source 2 for 10 • , 20 • , 30 • , 40 • and 50 • ) is studied. Figure 12(a) and Table 6 summarize the results of fault location. The fault location results of both the proposed and the existing methods do not change much with different loading conditions.

Effect of measurement noises
The effect of measurement noises (Gaussian distribution, with 0.2%, 0.5% and 1% standard deviations) is studied. Figure 12(b)  Table 7 summarize the results of fault location. The fault location errors generally increase with higher measurement noises. The fault location accuracy of the proposed method is still higher compared to the existing method.

Effect of parameter errors
The effect of parameter errors (all parameter matrices of the line with parameter errors of 0.2%, 0.5% and 1%) is studied. Figure 12(c) and Table 8 summarize the results of fault location. Similarly, the fault location errors generally increase with higher parameter errors. The fault location accuracy of proposed method is still higher compared to the existing method.

Effect of sampling rates
The effect of sampling rate (20 kHz, 50 kHz, 100 kHz and 200 kHz) is studied. Here a group of A-G faults with 0.01 ohm fault resistance is selected as examples. The fault location results are shown in Figure 13 and Table 9. From the figure, one can observe that when the sampling rate is relatively low (20 kHz or 50 kHz), the existing method does not have solution for faults close to the local terminal (in this case 2 km, since the corresponding dominant frequency is too high to be captured). In

Effect of line lengths
The effect of different line length (300 km, 400 km and 500 km) is studied. The fault location results are summarized in Figure 14 and Table 10. One can observe that the fault location accuracy of the proposed method does not change much, and the proposed method still presents high fault location accuracy compared to the existing method.

Effect of inception angles and advantages towards traveling wave based fault location methods
First, the effect of fault inception angles (including 0 • , 10 • , 30 • ) is studied and the fault location results are shown in Figure 15(a) and Table 11. One can observe that the fault location accuracy does not change much for both the existing and the proposed method, and the proposed method still presents high fault location accuracy compared to the existing method.
To further evaluate the performance of the proposed method during zero inception angle, a group of A-G faults with different fault resistances (1 ohm, 100 ohm, 500 ohm) is studied. The results are shown in Figure 15(b) and Table 12. One can observe that the proposed method can accurately locate faults even during extreme conditions, with zero inception angle and very high fault resistances. In fact, the validity of the traveling wave based methods is highly dependent on the reliable detection of the wavefronts. During faults with small inception angles or high fault resistances, the intensity of fault-created traveling wave is weak and the detection of the wavefronts could be problematic. For example, if the fault inception angle is close to zero, the wavefronts cannot be detected and the traveling wave based method will fail. On the contrary, the proposed method can still accurately locate faults even with zero fault inception angle. This is because the derivation procedure of the proposed method is based on frequency spectrum of the measurements, and is independent of the intensity of wavefronts. These results demonstrate clear advantages of the proposed method towards traveling wave based methods.

7.7
Advantages towards conventional single-ended impedance based method Comparing to the conventional single-ended impedance based fault location (IBFL) methods [27], the advantages of the proposed method are summarized as follows. (a) Short data window: Since the IBFL methods are based on fundamental frequency phasors, accurate extraction of phasors is mandatory for accurate IBFL results. However, if the transmission line is tripped with high speed after the fault, the available data window is short (e.g. 10 ms) and in this case the phasor estimation accuracy will be limited by the system transients and the decaying DC offset. In comparison, the proposed method is based on instantaneous measurements, does not require phasor extraction and is immune to decaying DC offset. Therefore, the proposed method can accurately locate faults with the short time window of 10 ms after the fault and is compatible with transmission lines using high-speed tripping techniques. (b) Accurate transmission line model: IBFL methods are usually based on lumped parameters of the transmission line, and therefore result in compromised fault location accuracy. In comparison, the derivation of the proposed method is based on distributed parameter model and the modelling accuracy of the transmission line is much improved. As a result, the proposed method presents higher fault location accuracy compared to the conventional IBFL methods especially for relatively long transmission lines.

7.8
Other discussions and future work In section 7.2, the measurement noises are assumed to be Gaussian distributed. In practice, if the system operates with conventional instrumentation current transformers (CTs) instead of optical CTs, the current measurements during faults may experience systematic errors due to CT saturations etc. In this case, the measurement errors could possibly compromise the fault location accuracy. Therefore, detailed modelling of CTs might be required to enable accurate recovery of primary currents during the fault [28] before applying the proposed natural frequency based fault location method. Besides, the proposed fault location methodology could potentially be applied to other transmission lines. Detail verifications of the proposed methodology for other transmission lines will be covered in future publications.

CONCLUSION
This paper proposes an improved natural frequency based fault location method for transmission lines. Only three phase instantaneous current measurements at the local terminal are required. The relationship between the theoretical peak frequencies and the fault location is first strictly derived. The fault location can be afterwards obtained by minimizing the average differences between the theoretical peak frequencies and measured peak frequencies. Compared to the existing natural frequency based fault location method, the proposed method systematically considers the mode mixing phenomenon and fault resistances, and fully utilizes the information embedded in the frequency spectrum of traveling waves instead of only the dominant frequency, resulting in improved fault location accuracy. Numerical experiments prove that compared to the existing method, the proposed method has much higher fault location accuracy during single phase to ground faults, and slightly higher (comparable) fault location accuracy during other types of faults, with different fault types, fault locations and fault resistances. The method also works with zero fault inception angle, where the system is without any fault-induced transient traveling wave.