An Accurate and Efficient Time Delay Estimation Method of Ultra-High Frequency Signals for Partial Discharge Localization in Substations

Partial discharge (PD) localization in substations based on the ultra-high frequency (UHF) method can be used to efficiently assess insulation conditions. Localization accuracy is affected by the accuracy of the time delay (TD) estimation, which is critical for PD localization in substations. A review of existing TD estimation methods indicates that there is a need to develop methods that are both accurate and computationally efficient. In this paper, a novel TD estimation method is proposed to improve both accuracy and efficiency. The TD is calculated using an improved cross-correlation algorithm based on full-wavefronts of array UHF signals, which are extracted using the minimum cumulative energy method and zero-crossing points searching methods. The cross-correlation algorithm effectively suppresses the TD error caused by differences between full-wavefronts. To verify the method, a simulated PD source test in a laboratory and a field test in a 220 kV substation were carried out. The results show that the proposed method is accurate even in case of low signal-to-noise ratio, but with greatly improved computational efficiency.


Introduction
The early faults of high-voltage equipment in substations can be detected effectively by locating the partial discharge (PD) source based on the ultra-high frequency (UHF) method [1]. The time difference of arrival (TDOA) method is commonly used as the localization algorithm [2]. The key of the TDOA method is to calculate the time delay (TD) between any two UHF signals captured by the sensor array. With three or four TDs and the coordinates of the sensor array, the localization equations can be established, from which the coordinate of the PD source can be obtained. The accuracy of the TDs has a direct impact on positioning accuracy [3]. Improving TD accuracy has received substantial attention in the field of PD localization in substations.
Conventional TD estimation methods include the threshold method [4], cumulative energy method [5], TD estimation method based on the Akaike information criterion (AIC) [6], and the cross-correlation algorithm [7]. The first three methods have a common feature, that is, the onset time is extracted from the UHF signal with respective methods, then TD is obtained by calculating the difference between the two onset times. The principle of the cross-correlation method is the construction of a cross-correlation function of the two UHF signals, and the offset of the maximum of the function is the TD of the two signals [8]. Theoretically, the TD can be accurately and efficiently calculated using these methods. However, in real substations, errors in TD estimation may occur due to two factors. First, the signal-to-noise ratio (SNR) of the UHF signal is reduced because of periodic narrowband noise and Gaussian white noise, and the onset time is changed with the SNR [9]. Therefore, the accuracy of the first three methods must be affected. Second, the consistency of the UHF array signals is severely reduced by the multipath effect [10,11], which affects the accuracy of the cross-correlation algorithm [12]; nevertheless, this algorithm has certain advantages when the SNR is low.
To obtain an accurate TD, many improved TD estimation methods have been developed. In [13], an improved threshold method was proposed, in which the high energy component of UHF signals was extracted using wavelet decomposition method, and the threshold value was set based on the wavelet coefficient. Subsequently, the onset time of the UHF signal could be obtained. However, the high energy component of the captured UHF signal is dominated by the field noise when the SNR is low, and hence the onset time cannot be obtained with reasonable accuracy. In [14], the continuous wavelet transform (CWT)-based binary map method was proposed, in which the onset time is determined by the earliest non-zero pixel of the binary map obtained by Otsu's method based on the CWT. The problem of the method is same to the threshold method, and the earliest non-zero pixel is affected by noise. In [6], part of the UHF signal was extracted using a time-domain window, which consists of the wavefront, and subsequently, TD is calculated using the cross-correlation algorithm based on the two extracted signals. When the effective wavefronts in the two parts are different, TD estimation accuracy is inevitably affected. In [15], a TD estimation method based on the interpolation cross-correlation algorithm was proposed to reduce the error generated by limiting the system sample rate, which enhanced the accuracy of TD estimation to a certain degree. However, the error caused by inconsistency in the array UHF signals has not been effectively solved. In [16], a TD estimation method using the cross-correlation algorithm based on a half-wavefront was proposed, which was used to locate the PD source in air-insulated substations. Two problems may affect this method: first, the half-wavefront extraction approach is complex, and this step takes most of the TD estimation time. Second, the half-wavefront extraction result is affected by the SNR. To resolve the influence of non-Gaussian noise on the cross-correlation algorithm, a method based on high-order statistics was proposed in [12]. This method effectively mitigates the effect of field noise. The UHF signal is divided into several segments, whose higher-order cumulants is calculated. The computational burden is large, which severely reduces efficiency. Designing a method that is both accurate and computationally efficient is desirable and challenging.
To improve the accuracy and efficiency of TD estimation and to develop a method applicable to realistic substation conditions, a novel TD estimation method is proposed in this paper. First, the initial time range of the effective UHF signal is obtained using the minimum cumulative energy method. Second, the full-wavefront of the UHF signal is extracted using the proposed zero-crossing point searching method. Third, TD is calculated using an improved cross-correlation algorithm. To verify the method, a simulated PD source was established in the laboratory, and array UHF signals with different SNRs were used as the samples for testing the method. A field test in a 220 kV substation with a real PD source was subsequently carried out. The accuracy and efficiency of the method were verified by comparing it with selected published methods.

Model of Array UHF Signals
To locate the PD source in the substation, a UHF sensor array was used to capture the signal emitted by the PD source. UHF signals typically comprise three components: the PD electromagnetic signal, periodic narrowband noise and the white noise, as shown in Equation (1). The y i (t) is the UHF signal received by the i-th sensor, v i (t) is the electromagnetic signal radiated from the PD source, g i (t) is the periodic narrowband noise, and w i (t) is the white noise. The frequency of UHF signals is within the range of 300 MHz-3 GHz [17]. Considering the frequencies of wireless communication signals used in China, the frequencies of the periodic narrowband noise are 470 MHz, 900 MHz and 1800 MHz [18]. White noise is generated by the heating effect of high-voltage apparatus: Because of the substation's complex environment, the propagation paths of UHF signals received by different sensors differ, as shown in Figure 1, where Ω represents the reflection surface, κ xx represents different propagation paths, and v i (t) is the superposition of the direct wave and several reflected waves. v i (t) can be expressed as Equation (2), where e i0 (t) is the direct wave; θ and ϕ are the azimuth angle and pitch angle, respectively, between the sensor and any point in the space; and h(θ 0 , ϕ 0 , t) is the transfer function of the sensor at bearing (θ 0 , ϕ 0 ). * is the convolution operator, ζ k is the attenuation coefficient of the electromagnetic signal, e(t + ∆t κ ) is the expression of the reflected wave in the time domain, ∆t κ is the TD between the direct wave and reflected wave, and h(θ κ , ϕ κ , t) is the transfer function of the sensor at bearing (θ κ , ϕ κ ). is the periodic narrowband noise, and wi(t) is the white noise. The frequency of UHF signals is within the range of 300 MHz-3 GHz [17]. Considering the frequencies of wireless communication signals used in China, the frequencies of the periodic narrowband noise are 470 MHz, 900 MHz and 1800 MHz [18]. White noise is generated by the heating effect of high-voltage apparatus: Because of the substation's complex environment, the propagation paths of UHF signals received by different sensors differ, as shown in Figure 1, where Ω represents the reflection surface, κxx represents different propagation paths, and vi(t) is the superposition of the direct wave and several reflected waves. vi(t) can be expressed as Equation (2), where ei0(t) is the direct wave; θ and φ are the azimuth angle and pitch angle, respectively, between the sensor and any point in the space; and h(θ0, φ0, t) is the transfer function of the sensor at bearing (θ0, φ0). * is the convolution operator, ζk is the attenuation coefficient of the electromagnetic signal, e(t + Δtκ) is the expression of the reflected wave in the time domain, Δtκ is the TD between the direct wave and reflected wave, and h(θκ, φκ, t) is the transfer function of the sensor at bearing (θκ, φκ).
The positions of the sensors in the sensor array differ, thus, the propagation paths of the UHF signals received by different sensors differ, and the bearings (θ, φ) of each signal differ, which would lead the inconsistencies of the array UHF signals. According to Figure 1, the propagation path of the reflected wave is longer than that of the direct wave, so the reflect wave is at the end of the UHF signal. According to Equation (2), if the reflection surfaces are different, the bearings (θκ, φκ) will be different. The gain of the sensor is related to directionality [19]. Therefore, the inconsistencies of the array UHF signals caused by the multipath effect exist at the end of the signals, and the wavefronts are unaffected or only slightly affected. Figure 2 depicts the array UHF signals radiated from a PD defect model in the laboratory. The signals were captured using two Vivaldi antennas [20] and were recorded by a LeCry WR640Zi oscilloscope with a bandwidth of 4 GHz and a sampling rate of 20 GS/s. The distance between the two antennas was 1.5 m, and the distances between the PD defect model and the two antennas were 4.24 and 3.71 m, respectively. The propagation speed of the electromagnetic wave was 3.00 × 10 8 m/s, and thus the TD of the two UHF signals was 1.767 ns, or 35-36 samples. To verify the effectiveness of

Measured Array UHF Signals
The positions of the sensors in the sensor array differ, thus, the propagation paths of the UHF signals received by different sensors differ, and the bearings (θ, ϕ) of each signal differ, which would lead the inconsistencies of the array UHF signals. According to Figure 1, the propagation path of the reflected wave is longer than that of the direct wave, so the reflect wave is at the end of the UHF signal. According to Equation (2), if the reflection surfaces are different, the bearings (θ κ , ϕ κ ) will be different. The gain of the sensor is related to directionality [19]. Therefore, the inconsistencies of the array UHF signals caused by the multipath effect exist at the end of the signals, and the wavefronts are unaffected or only slightly affected. Figure 2 depicts the array UHF signals radiated from a PD defect model in the laboratory. The signals were captured using two Vivaldi antennas [20] and were recorded by a LeCry WR640Zi oscilloscope with a bandwidth of 4 GHz and a sampling rate of 20 GS/s. The distance between the two antennas was 1.5 m, and the distances between the PD defect model and the two antennas were 4.24 and 3.71 m, respectively. The propagation speed of the electromagnetic wave was 3.00 × 10 8 m/s, and thus the TD of the two UHF signals was 1.767 ns, or 35-36 samples. To verify the effectiveness of the TD estimation method at a low SNR, the field noise was measured at the substation, and the noise was added to the array UHF signals. The noisy UHF signals are presented in Figure 2b, and the SNR was 3.86 dB.

Measured Array UHF Signals
Sensors 2018, 18, x 4 of 15 the TD estimation method at a low SNR, the field noise was measured at the substation, and the noise was added to the array UHF signals. The noisy UHF signals are presented in Figure 2b, and the SNR was 3.86 dB.
(a) (b) Figure 2. Array UHF signals radiated from the PD defect model and the noisy array UHF signals: (a) array UHF signals; (b) noisy array UHF signals.

TD Estimation Method Using the Improved Cross-Correlation Algorithm Based on the Full-Wavefront
According to the analysis in Section 2.1.1, the wavefront of the UHF signal has the least distortion; furthermore, it contains the largest proportion of the direct wave energy. If the fullwavefront is used for TD estimation instead of the complete UHF signal, then the error caused by the multipath effect can be minimized. The cross-correlation algorithm can reduce the influence of uncorrelated noise, leading to more accurate TD estimate. However, the effect of signal consistency on the result must be considered. The methods for accurately and effectively extracting the wavefront from the noisy UHF signal and improving the cross-correlation algorithm are key to the TD estimation method.

Full-Wavefront Extraction Method
To obtain an accurate full-wavefront, periodic narrowband noise and white noise should be considered. The full-wavefront extraction procedure is as follows: Step 1: Suppressing the white noise w(t). The original UHF signal yi(n) is filtered using a neighborhood smoothing filter, and the signal is recorded as yi'(n). The neighborhood smoothing filter is defined in Equation (3), where M is the width of the neighborhood: Step 2: Suppressing the periodic narrowband noise g(t).
Step 2.1: The cumulative energy curve of yi'(n) is calculated using Equation (4); subsequently, the minimum cumulative energy curve of yi'(n) is calculated using Equation (5), where EN is the maximum of E(n). The sampling number of the minimum of Emin(n) is the apparent onset time of yi(n); the sampling number is subsequently extracted and recorded as Nmin. Note that the apparent onset time is not the actual onset time, which is affected by the noise and the procedure in Step 1.
Step 2.2: The noise frame of the UHF signal w'(n) is extracted from the first sampling point to Nmin. The zero-crossing points and crossing direction are obtained using the Horizontal axis10ns/div Vertical axis 25mV/div CH1 CH2 CH1

TD Estimation Method Using the Improved Cross-Correlation Algorithm Based on the Full-Wavefront
According to the analysis in Section 2.1.1, the wavefront of the UHF signal has the least distortion; furthermore, it contains the largest proportion of the direct wave energy. If the full-wavefront is used for TD estimation instead of the complete UHF signal, then the error caused by the multipath effect can be minimized. The cross-correlation algorithm can reduce the influence of uncorrelated noise, leading to more accurate TD estimate. However, the effect of signal consistency on the result must be considered. The methods for accurately and effectively extracting the wavefront from the noisy UHF signal and improving the cross-correlation algorithm are key to the TD estimation method.

Full-Wavefront Extraction Method
To obtain an accurate full-wavefront, periodic narrowband noise and white noise should be considered. The full-wavefront extraction procedure is as follows: Step 1: Suppressing the white noise w(t). The original UHF signal y i (n) is filtered using a neighborhood smoothing filter, and the signal is recorded as y i (n). The neighborhood smoothing filter is defined in Equation (3), where M is the width of the neighborhood: Step 2: Suppressing the periodic narrowband noise g(t).
Step 2.1: The cumulative energy curve of y i (n) is calculated using Equation (4); subsequently, the minimum cumulative energy curve of y i (n) is calculated using Equation (5), where E N is the maximum of E(n). The sampling number of the minimum of E min (n) is the apparent onset time of y i (n); the sampling number is subsequently extracted and recorded as N min . Note that the apparent onset time is not the actual onset time, which is affected by the noise and the procedure in Step 1.
Step 2.2: The noise frame of the UHF signal w (n) is extracted from the first sampling point to Nmin. The zero-crossing points and crossing direction are obtained using the method shown in Figure 3, where L(j) is the sampling number of the zero-crossing point, F(j) is the crossing direction of the zero-crossing point, F(j) = 1 indicates the rising edge, and F(j) = −1 indicates the falling edge.

End
Calculating the product pro Figure 3. Flow chart of the zero-crossing point searching method.
Step 2.3: Depending on the results acquired in Step 2.2, the signal w"(n) is extracted from the noise frame w'(n) between the first zero-crossing point L(0) and the last zerocrossing point, which has the same crossing direction as L(0). The signal w"(n) is extended r times, and the length of the extended signal is checked to ensure it is equal to the length of yi(n) from point L(0) to the end of the signal, which is recorded as w"'(n). The signal from yi'(n) between L(0) and the end is extracted; subsequently, the difference between the extracted signal and w"'(n) is calculated and recorded as yi-dif'(n). Subsequently, L(0) zeros are padded at the front of yi-dif'(n), which is recorded as yi"(n), and the signal yi"(n) is the UHF signal after the periodic narrowband noise is suppressed.
Step 3: Extracting the full-wavefront Step 3.1: The minimum cumulative energy curve of yi"(n) is calculated using Equations (4) and (5), and the sampling number N"min is extracted, which corresponds to the minimum of the curve. Two signals (duration: 3 ns) are extracted before and after N"min, and the extracted signal is the effective signal frame, which is recorded as yiwf"(n).
Step 3.2: The signal segment of yi"(n) is extracted from the start of yi"(n) to the sample N"min, and the mean value of the absolute value of all crests and troughs is calculated, which is recorded as Aave. The zero-crossing points of yi"(n) are obtained using the Step 2.3: Depending on the results acquired in Step 2.2, the signal w (n) is extracted from the noise frame w (n) between the first zero-crossing point L(0) and the last zero-crossing point, which has the same crossing direction as L(0). The signal w (n) is extended r times, and the length of the extended signal is checked to ensure it is equal to the length of y i (n) from point L(0) to the end of the signal, which is recorded as w (n). The signal from y i (n) between L(0) and the end is extracted; subsequently, the difference between the extracted signal and w (n) is calculated and recorded as y i-dif (n). Subsequently, L(0) zeros are padded at the front of y i-dif (n), which is recorded as y i (n), and the signal y i (n) is the UHF signal after the periodic narrowband noise is suppressed.
Step 3: Extracting the full-wavefront Step 3.1: The minimum cumulative energy curve of y i (n) is calculated using Equations (4) and (5), and the sampling number N min is extracted, which corresponds to the minimum of the curve. Two signals (duration: 3 ns) are extracted before and after N min , and the extracted signal is the effective signal frame, which is recorded as y iwf (n). Step 3.2: The signal segment of y i (n) is extracted from the start of y i (n) to the sample N min , and the mean value of the absolute value of all crests and troughs is calculated, which is recorded as A ave . The zero-crossing points of y i (n) are obtained using the method presented in Figure 3, and the absolute value of the crest and trough between every two zero-crossing points is counted. With δ·A ave as the threshold, to find the first absolute value that is larger than δ·A ave , the sampling number of the previous zero-crossing point before the absolute value is the start of the full-wavefront. The third zero-crossing point is the end of the full-wavefront. With the start and the end samples, the full-wavefront is extracted from the original signal y(n). This indicates that the full-wavefront contains a crest and a trough.
Taking the first signal in Figure 2b as an example, the full-wavefront extraction procedure is illustrated in Figure 4.  Figure 3, and the absolute value of the crest and trough between every two zero-crossing points is counted. With δ·Aave as the threshold, to find the first absolute value that is larger than δ·Aave, the sampling number of the previous zero-crossing point before the absolute value is the start of the fullwavefront. The third zero-crossing point is the end of the full-wavefront. With the start and the end samples, the full-wavefront is extracted from the original signal y(n). This indicates that the full-wavefront contains a crest and a trough.
Taking the first signal in Figure 2b as an example, the full-wavefront extraction procedure is illustrated in Figure 4.   Figure 4a presents the result of suppressing the white noise, and the neighborhood smoothing filter parameter M is set as 5. Figure 4b presents the result of suppressing the periodic narrowband noise and the full-wavefront extraction process. In Step 2, the value of δ is set as 1.1. Figure 4c shows the extracted full-wavefront.

TD Calculation Using the Improved Cross-Correlation Algorithm
As explained in Section 2.1, the values of the wavefronts of the array UHF signals differ because of the difference in the bearings and directionality of the UHF sensor. During full-wavefront extraction, only the parameter δ needs to be set manually. Because of the difference among the array UHF signals, if the parameter δ is unreasonable, then the array full-wavefronts will deviate. Assume that the value of δ·A ave is greater than the amplitude of the first half wavefront; consequently, the full-wavefront extraction result will deviate, as described in Figure 5. The wave in Figure 5 can be divided into three half-waves: A(t), B(t), and C(t). The corect full-wavefront is A(t) + B(t), and the full-wavefront with deviation is A(t) + C(t). Only if the extracted full-wavefronts are all correct or have the same deviations can the correct TD be obtained using the cross-correlation algorithm. However, the parameter δ·A ave cannot satisfy all the UHF signals.  Figure 4a presents the result of suppressing the white noise, and the neighborhood smoothing filter parameter M is set as 5. Figure 4b presents the result of suppressing the periodic narrowband noise and the full-wavefront extraction process. In Step 2, the value of δ is set as 1.1. Figure 4c shows the extracted full-wavefront.

TD Calculation Using the Improved Cross-Correlation Algorithm
As explained in Section 2.1, the values of the wavefronts of the array UHF signals differ because of the difference in the bearings and directionality of the UHF sensor. During full-wavefront extraction, only the parameter δ needs to be set manually. Because of the difference among the array UHF signals, if the parameter δ is unreasonable, then the array full-wavefronts will deviate. Assume that the value of δ·Aave is greater than the amplitude of the first half wavefront; consequently, the fullwavefront extraction result will deviate, as described in Figure 5. The wave in Figure 5 can be divided into three half-waves: A(t), B(t), and C(t). The corect full-wavefront is A(t) + B(t), and the fullwavefront with deviation is A(t) + C(t). Only if the extracted full-wavefronts are all correct or have the same deviations can the correct TD be obtained using the cross-correlation algorithm. However, the parameter δ·Aave cannot satisfy all the UHF signals.
To obtain the correct TD when the two wavefronts are inconsistent, the cross-correlation algorithm should be improved. The half wavefront is recorded as vwf(t); assume that the center times of A(t), B(t), and C(t) are t, t − T1, and t + T2, respectively. The expressions of the three half-waves are presented in Equation (6) In the case of deviation, the expression of the two full-wavefronts s1(t) and s2(t) is as Equation (7), where ζ is the ratio of the amplitudes of the two wavefronts: The cross-correlation function of s1(t) and s2(t) is shown as Equation (8), where E is the expected operator: To obtain the correct TD when the two wavefronts are inconsistent, the cross-correlation algorithm should be improved. The half wavefront is recorded as v wf (t); assume that the center times of A(t), B(t), and C(t) are t, t − T 1 , and t + T 2 , respectively. The expressions of the three half-waves are presented in Equation (6), where a, b, and c are the amplitudes of A(t), B(t), and C(t), respectively.
In the case of deviation, the expression of the two full-wavefronts s 1 (t) and s 2 (t) is as Equation (7), where ζ is the ratio of the amplitudes of the two wavefronts: The cross-correlation function of s 1 (t) and s 2 (t) is shown as Equation (8), where E is the expected operator: The result of Equation (8) is the sum of the cross-correlation functions, A(t) and A(t + τ), A(t) and C(t + τ), A(t + τ) and B(t), and B(t) and C(t + τ). According to the principle of the cross-correlation algorithm, when t = −τ, the cross-correlation function of A(t) and A(t + τ) reaches the extremum value of ζ·a 2 , which is the correct TD of s 1 (t) and s 2 (t). The time and values of the extremum of the four cross-correlation functions are listed in Table 1. The data in the table indicate that ζ·a 2 and ζ·bc are both positive numbers, and that ζ·ab and ζ·ac are both negative numbers. Thus, theoretically, there are two crests and two troughs on the cross-correlation curve. Because the values of T 1 and T 2 are close, the two troughs overlap. The amplitudes of the two crests are related to ζ·a 2 and ζ·bc, respectively. The correct TD corresponds to the crest whose amplitude is ζ·a 2 , and the correct TD can be obtained by judging the magnitude of ζ·a 2 and ζ·bc.

Cross-Correlation Function Time of the Extremum Value of the Extremum
How to determine parameter a is key to the improved cross-correlation algorithm. Per the results of Step 3.1 in Section 2.2.1, the sampling number N min is close to the actual onset time of the UHF signal. The difference between N min and the sample point of the first extremum was calculated, and the wavefront with a smaller difference was the correct one to be used for determining parameter a.
Consider the signals in Figure 2b as an example. The full-wavefront and cross-correlation curve are presented in Figure 6a,b, respectively. The peak of the cross-correlation curve occurs at sample point 965. The numbers of samples for the two UHF signals are both 1000. The estimated number of samples for the TD result is 35, which is the correct value. The full-wavefronts of another two UHF signals with deviation and minimum cumulative energy curve are listed in Figure 6c, and the full-wavefronts are recorded as s 1wf (t) and s 2wf (t). The cross-correlation curve is illustrated in Figure 6d. According to the aforementioned estimation method, s 1wf (t) is the correct full-wavefront. The amplitudes of the crests and troughs of the two full-wavefronts are determined: max(s 1wf (t)) = 29.40 mV, min(s 1wf (t)) = −48.76 mV, max(s 2wf (t)) = 26.83 mV, and min(s 1wf (t)) = −30.38 mV. Three extremums are shown on the curve. This indicates that a 2 = min(s 1wf (t)) · min(s 1wf (t)) and bc = max(s 1wf (t)) · max(s 2wf (t)). Because a 2 > bc, the sampling point of the second peak is selected, which is 964. The TD estimation result is 36, which is consistent with the real TD.

B(t), C(t + τ)
-T1-T2-τ ζ·bc How to determine parameter a is key to the improved cross-correlation algorithm. Per the results of Step 3.1 in Section 2.2.1, the sampling number N''min is close to the actual onset time of the UHF signal. The difference between N''min and the sample point of the first extremum was calculated, and the wavefront with a smaller difference was the correct one to be used for determining parameter a.

Comparative Analysis of Several TD Estimation Methods
To verify the accuracy and efficiency of the proposed method in this paper, the performance of the proposed method was compared with that of four published methods: the minimum cumulative energy method [5], cross-correlation algorithm, cross-correlation algorithm based on half-wavefront [16], and TD estimation method based on high-order statistics [12]. The four methods are described in Table 2.
With the array UHF signals in Figure 2b as an example, the TD estimation results provided by the four methods are listed in Figure 7. In Figure 7a, the TD is 1.55 ns, calculated using Method A, and the error is 5 samples. In Figure 7b, the number of samples at the peak of the cross-correlation function curve is 1160. The TD is −160 samples, and the error is because of the inadequate consistency of the two UHF signals. The TD calculated using Method C is 73 sampling points, as shown in Figure 7c, and the error is 38 samples, which is because of the deviation between the two half-wavefronts. In Figure 7d, the TD calculated using Method D is 36, which is consistent with the actual TD.

TD Estimation Method Method Name Principle
Method A Minimum cumulative energy method (1) The minimum cumulative energy curve is calculated using Equations (4) and (5).
(2) The sampling number correspond to the minimum of the curve is determined, which is the start of the UHF signal frame.
(3) The difference in the sampling numbers of the two signals is calculated, which is the TD.

Method B
Cross-correlation algorithm (1) The cross-correlation function of two signals is calculated.
(2) The offset of the maximum value of the cross-correlation function is determined, which is the TD.

TD Estimation Method Method Name Principle
Method C TD estimation method proposed in [16] (1) The time-frequency transform of the UHF signal is carried out.
(2) The signal in the range 0.5-0.8 GHz is extracted, and the signal is reconstructed.
(3) The signal containing the wavefront is extracted on the basis of the reconstructed signal, and the number and positions of the crest and trough are found.
(4) The half-wavefront is extracted, and the TD is calculated using the cross-correlation algorithm.

Method D TD estimation method based on high-order statistics
The principle of this method is described in [12].
Sensors 2018, 18, x 10 of 15 high-order statistics Fifty array UHF signals were measured in the laboratory, and field noise of various amplitudes was added. These signals were then used to verify the accuracy and efficiency of the five methods. The SNRs were set to 25, 20, 15, 10 and 5 dB. The means and standard deviations of the 50 TD estimation results were used to evaluate accuracy, and the total computational time was used to evaluate efficiency. The computer configuration was as follows: an Intel Core i5-4460 system with 3.2 GHz of clock speed, 8 GB of RAM, and a 64-bit Windows 7 operating system. The means, standard deviations, and time consumptions for the 50 array UHF signals are listed in Table 3.   Fifty array UHF signals were measured in the laboratory, and field noise of various amplitudes was added. These signals were then used to verify the accuracy and efficiency of the five methods. The SNRs were set to 25, 20, 15, 10 and 5 dB. The means and standard deviations of the 50 TD estimation results were used to evaluate accuracy, and the total computational time was used to evaluate efficiency. The computer configuration was as follows: an Intel Core i5-4460 system with 3.2 GHz of clock speed, 8 GB of RAM, and a 64-bit Windows 7 operating system. The means, standard deviations, and time consumptions for the 50 array UHF signals are listed in Table 3.
As the SNR decreased, the TD estimation accuracy of the five methods also decreased. For SNR = 25 dB, the mean value derived using Method B deviated substantially from the actual TD, whereas other methods could obtain accurate results. This indicates that the array signal consistency has stronger effect on the cross-correlation algorithm. When the SNR ≤ 10 dB, the standard deviations of Methods A and C increased sharply, indicating that field noise has a greater effect on these two methods than on other methods. For Method A, the position of the trough in the minimum energy accumulation curve changed depending on field noise. For Method C, the amplitude of the noise exceeded the amplitude of the half-wavefront, which caused the deviation of the half-wavefront extraction. In this case, the conventional cross-correlation algorithm will fail to obtain the correct time delay. The proposed method in this paper and Method D are more accurate than the other three methods. This indicates that the two methods suppress noise and array signal consistency effectively. The time consumptions of the proposed method and methods A, B, and C have the same order of magnitude, and the time consumptions are all less than 200 ms. The time consumption of Method C is greater than that of the proposed method because the time-frequency analysis is used to extract the half-wavefront, which is a time-intensive process. The time required for Method D is 3860.5 s, which is 40,213 times as much as that for the proposed method. Because numerous matrix operations are performed during the TD calculation, the efficiency of the algorithm is severely reduced. By contrast, only the basic array operations are used in the proposed method, thus greatly improving the computational efficiency. In summary, the proposed TD estimation method is highly accurate and efficient for array UHF signals with a low SNR.

Test Arrangement
To verify the effectiveness of proposed method, a field test in a 220 kV substation was carried out. The field test setup is shown in Figure 8a. The measurement system comprised an antenna array, a rotating platform, and a data acquisition system. The parameters of the antenna array and oscilloscope were identical to those described in Section 2.1.2. The antenna array was placed on the rotating platform, whose center coincides with that of the antenna array. A PD source was placed in a disc insulator with an air-gap defect, as marked by the red ellipse in Figure 8a. The result of the online monitoring system is presented in Figure 8b.
In the substation, the coordinate system was established, and the plane where the array sensors are placed is set as the xoy plane. The location range of the disc insulator was {x|x ∈ (3.3, 3.4)}, {y|y ∈ (8.55, 9)}, and {z|z ∈ (2.6, 3)} in meters. The origin point was selected as the measurement point, o (0, 0, 0). One of the measured array UHF signals is listed in Figure 8d, and the rotation angle of the antenna array is −15 • . The TD range of the two UHF signals is 0.5-0.6 ns, and the corresponding samples are theoretically 10-12. The amplitude of the field noise is close to that of the signal, and the SNR is low.
out. The field test setup is shown in Figure 8a. The measurement system comprised an antenna array, a rotating platform, and a data acquisition system. The parameters of the antenna array and oscilloscope were identical to those described in Section 2.1.2. The antenna array was placed on the rotating platform, whose center coincides with that of the antenna array. A PD source was placed in a disc insulator with an air-gap defect, as marked by the red ellipse in Figure 8a. The result of the online monitoring system is presented in Figure 8b. In the substation, the coordinate system was established, and the plane where the array sensors are placed is set as the xoy plane. The location range of the disc insulator was {x|x  (3.3, 3.4)}, {y|y  (8.55, 9)}, and {z|z  (2.6, 3)} in meters. The origin point was selected as the measurement point, o (0, 0, 0). One of the measured array UHF signals is listed in Figure 8d, and the rotation angle of the antenna array is −15°. The TD range of the two UHF signals is 0.5-0.6 ns, and the corresponding samples are theoretically 10-12. The amplitude of the field noise is close to that of the signal, and the SNR is low.

Result Analysis of the TD Estimation
The TD is calculated as follows: First, the full-wavefronts are extracted, as shown in Figure 9a, and a deviation separates the two full wavefronts. Figure 9b illustrates the cross-correlation function curve. According to Section 2.2.2, the sampling point of the maximum is 4011; thus, the TD is 11 sampling points (0.55 ns). Figure 10 presents a comparison of the TD estimation result of the four published methods. The calculated TD was 1.05, 3.4, 0.5, and 0.55 ns. The results calculated using Methods C and D are accurate. At this rotating angle, TD estimation results obtained considering 100 array UHF signals as samples are shown in Table 4. Method B has the largest error. The SNR of some UHF signals is excessively low, which leads to the deviation in the inflection point and the halfwavefront. Thus, the means of TDs calculated using Methods A and B deviate from the actual TD.

Result Analysis of the TD Estimation
The TD is calculated as follows: First, the full-wavefronts are extracted, as shown in Figure 9a, and a deviation separates the two full wavefronts. Figure 9b illustrates the cross-correlation function curve. According to Section 2.2.2, the sampling point of the maximum is 4011; thus, the TD is 11 sampling points (0.55 ns). Figure 10 presents a comparison of the TD estimation result of the four published methods. The calculated TD was 1.05, 3.4, 0.5, and 0.55 ns. The results calculated using Methods C and D are accurate. At this rotating angle, TD estimation results obtained considering 100 array UHF signals as samples are shown in Table 4. Method B has the largest error. The SNR of some UHF signals is excessively low, which leads to the deviation in the inflection point and the half-wavefront. Thus, the means of TDs calculated using Methods A and B deviate from the actual TD. In the substation, the coordinate system was established, and the plane where the array sensors are placed is set as the xoy plane. The location range of the disc insulator was {x|x  (3.3, 3.4)}, {y|y  (8.55, 9)}, and {z|z  (2.6, 3)} in meters. The origin point was selected as the measurement point, o (0, 0, 0). One of the measured array UHF signals is listed in Figure 8d, and the rotation angle of the antenna array is −15°. The TD range of the two UHF signals is 0.5-0.6 ns, and the corresponding samples are theoretically 10-12. The amplitude of the field noise is close to that of the signal, and the SNR is low.

Result Analysis of the TD Estimation
The TD is calculated as follows: First, the full-wavefronts are extracted, as shown in Figure 9a, and a deviation separates the two full wavefronts. Figure 9b illustrates the cross-correlation function curve. According to Section 2.2.2, the sampling point of the maximum is 4011; thus, the TD is 11 sampling points (0.55 ns). Figure 10 presents a comparison of the TD estimation result of the four published methods. The calculated TD was 1.05, 3.4, 0.5, and 0.55 ns. The results calculated using Methods C and D are accurate. At this rotating angle, TD estimation results obtained considering 100 array UHF signals as samples are shown in Table 4. Method B has the largest error. The SNR of some UHF signals is excessively low, which leads to the deviation in the inflection point and the halfwavefront. Thus, the means of TDs calculated using Methods A and B deviate from the actual TD.  The TD obtained using Method C is close to the actual result, while the standard deviation is large. This indicates that some of the calculation results have large deviations. Both Method D and the proposed method have obtained TD with higher accuracy. In terms of computational efficiency, time consumption is 41,270 times as much as that of the proposed method.

Conclusions
In this paper, a TD estimation method for locating PD sources in substations is proposed, in which the full-wavefront of the UHF signal is extracted to replace the original UHF signal. The crosscorrelation algorithm was improved to enhance TD accuracy. A simulated model was established, and a field test was conducted to verify the method. The conclusions are summarized as follows: 1) The full-wavefront was extracted for calculating the TD, which reduces the influence of the field noise and the inconsistency of the array UHF signals and thus effectively improves the accuracy of TD estimation. 2) TD estimation efficiency was enhanced using the proposed full-wavefront extraction method and improved cross-correlation algorithm, making it more suitable for field applications. 3) An experimental platform with a PD defect model was established in the laboratory. Three hundred array UHF signals with different SNRs were used as the samples to verify the proposed method, and the results of the proposed method were compared with that of four published methods. The proposed method yielded accurate TD estimation based on high-  The TD obtained using Method C is close to the actual result, while the standard deviation is large. This indicates that some of the calculation results have large deviations. Both Method D and the proposed method have obtained TD with higher accuracy. In terms of computational efficiency, time consumption is 41,270 times as much as that of the proposed method.

Conclusions
In this paper, a TD estimation method for locating PD sources in substations is proposed, in which the full-wavefront of the UHF signal is extracted to replace the original UHF signal. The cross-correlation algorithm was improved to enhance TD accuracy. A simulated model was established, and a field test was conducted to verify the method. The conclusions are summarized as follows: (1) The full-wavefront was extracted for calculating the TD, which reduces the influence of the field noise and the inconsistency of the array UHF signals and thus effectively improves the accuracy of TD estimation. (2) TD estimation efficiency was enhanced using the proposed full-wavefront extraction method and improved cross-correlation algorithm, making it more suitable for field applications. (3) An experimental platform with a PD defect model was established in the laboratory. Three hundred array UHF signals with different SNRs were used as the samples to verify the proposed method, and the results of the proposed method were compared with that of four published methods. The proposed method yielded accurate TD estimation based on high-order statistics.
Regarding efficiency, computational time was reduced from 3,860,500 ms to 96 ms, and the efficiency increased more than 40,000 times. (4) A field test was carried in a 220 kV substation. An air-gap discharge in a disc insulator was used as the testing target. Compared with the four previously published methods, the proposed method yielded the most accurate TD with lowest time requirement, indicating the high accuracy and efficiency of the proposed method.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: PD partial discharge UHF ultra-high frequency TD time delay SNR signal noise ratio TDOA time difference of arrival CWT continuous wavelet transform