A Location Method for Partial Discharge Using Time Reversal and Improved Whale Optimization Algorithm

Time difference of arrival method is commonly used in the location of partial discharges. In this method, the main problem are as follows: (1) time delay error is inevitable. (2) location algorithm is sensitive to time delay error. To solve the problems above, an ultra-high frequency localization method for transformer partial discharge based on time reversal and improved whale optimization algorithm is proposed in this article. The method is consisted of three parts, that including time delay extraction, time reversal imaging and optimal solution. Firstly, the time delay is extracted by energy accumulation method. And then the approximate position of partial discharge source is determined by time reversal imaging. Combined with the time reversal imaging, the strategy of initial population selection and weight setting in the whale optimization algorithm is improved. Finally, the improved whale optimization algorithm is used to solve the optimal solution. The results are demonstrated through simulation and experiment. Compared with other localization algorithms, it is shown that the improved whale optimization algorithm has a better control performance in stability and accuracy.


I. INTRODUCTION
Deterioration degree of power equipment's insulation performance can be expressed by partial discharge [1]. At the same time, PD (partial discharge, PD) will exacerbate the deterioration of insulation [2]. Thus, to ensure the safety and reliability of power equipment, it is essential to conduct partial discharge testing on power equipment regularly, find partial discharges in time, and locate the discharge source accurately [3]. Therefore, the location of PD source in electrical equipment has always been a hot spot in the field of PD research.
On the other hand, time difference of arrival algorithm is the basic method of UHF (ultra-high frequency, HUF) PD location [4]. The principle of TDOA (time difference of arrival, TDOA) method is simple, in which the accuracy of The associate editor coordinating the review of this manuscript and approving it for publication was Dazhong Ma . time delay extraction determines the accuracy of PD location [5]. At present, many methods of time delay extraction have been proposed, such as initial peak method [6], cross-correlation function method [7] and energy accumulation method [8]. The determination of the initial peak is affected by many factors, such as the sampling rate of the oscilloscope, on-site interference and the experience of the testing personnel [9]. The cross-correlation function method has a relatively accurate delay estimation value at low SNR (signal-to-noise ratio, SNR), but a high SNR may cause a large error in the time delay estimation [10]. Owing to the influence of interference and other conditions, the inflection point may appear in energy accumulation curve before the arrival of the pulse signal. This leads to the misjudgment of the signal starting time [11]. The methods above are mainly to distinguish and extract the waveform features, which inevitably have some errors. Therefore, it can only default the existence of the error to continue to solve the positioning VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ algorithm. In reference [12], PSO (particle swarm optimization, PSO) algorithm is used to solve PD positioning equation and the effectiveness of the algorithm is verified by an example. In reference [11], the particle swarm optimization algorithm is improved by combining Tabu search algorithm with PSO, which avoids the problems of non-convergence, non-uniqueness of solutions and high initial value requirements of the least square method. In reference [13], the genetic algorithm is used to find the appropriate initial point for Taylor expansion, in which the Taylor algorithm is used to solve the optimal solution. Hence, the error caused by the randomness of the initial point is avoided. However, the positioning algorithm is extremely sensitive to the time delay error.
In order to reduce the influence of interference signal on PD location, BSS (blind source separation, BSS) is used to extract PD signal in [14]. In [15], BSS and kurtosis are used to identify the PD sources fast. In [16], BSS is applied to TDOA location method. It not only improves the accuracy of time delay extraction, but also reduces the number of iterations of the positioning algorithm.
TR (time reversal, TR) technique is adaptive, which can achieve the focus of the signal source when the transmission characteristics of the transmission medium and sensor are unknown [17]. TR imaging technology can ignore the influence of waveform mode, and does not need to calculate the time delay. The position information of signal can be obtained by TR imaging. The steps of TR imaging are as follows: (1) Record the signal from the signal source received by different sensors.
(2) The received signal is transformed by TR. Then the signal is transmitted back to signal source from the corresponding sensor. This is an effective method to reconstruct the signal source [18].
Motivated by the above discussion, this article has considered the application of TR imaging method to PD localization. It is combined with IWOA (the improved whale optimization algorithm, IWOA) to further improve the accuracy of TDOA positioning algorithm. The present work seeks to improve the positioning accuracy of PD by making the following contributions.
(1) TR imaging is performed on PD source. The results of TR imaging before and after BSS are observed. The TR imaging is compared with the time delay extraction method, such as initial peak value method, cross correlation method and accumulated energy method.
(2) The initial population and weight strategy of WOA are set by TR imaging. The fitness function of TDOA localization algorithm is improved and applied to WOA. The TR_IWOA algorithm is applied to the TDOA method of UHF PD location.
(3) Through the simulation and experiment, the TR_IWOA is compared with particle swarm optimization and least square method in TDOA. The results have verified the feasibility and accuracy of TR_IWOA.

II. PRINCIPLE OF PARTIAL DISCHARGE LOCATION
At present, the common positioning method in engineering field is TDOA. In this article, the time reversal imaging is also applied to PD location.

A. PRINCIPLE OF PARTIAL DISCHARGE LOCATION BY TIME DIFFERENCE OF ARRIVAL
When PD occurs in the transformer, the electromagnetic wave generated by the PD signal will propagate in the form of spherical waves to the surroundings. At last, it will be received by UHF sensors in different positions through different propagation paths. Fig. 1 is the location principle diagram of TDOA method. Assume that there are N sensors in total. Considering the interference of noise during propagation, the time for the PD signal to reach sensor i is: where t i is the theoretical time for the PD signal to reach the sensor, n t i is the error due to noise interference subject to a Gaussian distribution with mean 0 and variance δ 2 d . Suppose d i is the distance from PD signal to S i and d i,j is the measurement value of the distance difference between PD signal to S i and S j , one has where i, j = 1, 2, 3, 4 and i = j, n d i is the noise distance error, subject to a Gaussian distribution with mean 0 and variance δ 2 d . According to the distance difference between the PD signal and two different sensors, a hyperbolic equation can be established. The location of the PD source can be obtained by solving the hyperbolic equations established by multiple sensors through an optimization algorithm.

B. PRINCIPLE OF TIME REVERSAL FOCUS IMAGING
Transmission media and sensors form a signal transmission system [19]. If the sensor receives the signal from the signal source A(ω), the received signal S(ω) can be expressed as follows: where H (r, ω) is the transfer function of the transmission medium, r is the distance from the sensor to the signal source and ω is the angular frequency.
In the time domain, the signal S(t) is transformed by TR, that is, conjugate in the frequency domain of the signal S(ω). the signal S(−t) after TR transformation can be expressed as follows: where * is the convolution symbol. According to the principle of reciprocity, when the propagation path is unchanged, the position of the sensor and the signal source can be exchanged and the frequency response is unchanged. Therefore, the signal S * (ω) is sent by the corresponding sensor and the signal C(ω) received at the signal source can be expressed as follows: where H * (r, ω) · H (r, ω) is a real, even and positive function. Inverse Fourier transforms of the same phase will superimpose each other at time zero and generate the main correlation peak. When signals received by several sensors are all operated by this method, multiple signals are accumulated simultaneously. The signal amplitude where the signal occurs is significantly enhanced due to focusing effect. This is the basic principle of TR focusing. When PD occurs at the point P, the UHF signal generated will propagate along the path of Fig. 1 and finally be received by the sensor. Assume that the PD signal is sent out at time t 0 . Sensor i last received the signal. Take t N +1 > t i . The UHF signal in t N +1 − t 0 time period is intercepted for TR transformation, and the resulting signal is f n . Nf n are sent out on the corresponding sensor to focus the signal at the PD source P at the same time. Then the focusing time t f can be expressed as: If one point Q(x i , y j , z k ) is selected in the whole detection space, then the energy amplitude of the point at the focusing time can be obtained as follows: where A ijk is the energy amplitude of point Q by superposition, K is the amplification factor of the sensor, f n (t nijk ) is the amplitude of the signal received by the sensor n after TR transformation to point Q, C n (x n , y n , z n ) is the coordinate of the sensor N and v is the signal propagation speed. Fig. 2 is a flow chart of amplitude focusing of PD UHF signal based on TR focusing imaging.

III. IMPROVED WHALE OPTIMPROVED ALGORITHM BASED ON TIME REVERSAL
In TR imaging, the point with a larger focus amplitude and a darker color can be considered to be the more likely the location of PD source. Therefore, TR imaging can be regarded as a probability distribution. The larger the amplitude, the greater the probability of PD source. In view of the situation that the delay error causes the optimization algorithm to fall into a local optimal or no solution, combined with the probability distribution, the initial population and weight strategy of TDOA location algorithm are improved. On the other hand, the spiral optimization path in WOA is similar to TR imaging, as shown in the Fig. 3. Therefore, an improved whale optimization algorithm suitable for TODA PD location is proposed by combining TR imaging with WOA. Detailed improvements are as follows.

A. WHALE OPTIMIZATION ALGORITHM
Inspired by the method of whale predation [20]. WOA mainly includes three steps: surrounding predation, bubble attack and random food search

1) STAGE OF SURROUNDING PREDATION
Humpback whales can identify the location of prey and circle it. Since the location of the optimal design in the search space is not a priori known, WOA assumes that the current optimal candidate solution is the target prey or close to the optimal solution. After the best whale position is defined, other whales will try to update their position to the best whale position. This behavior is expressed by the following equations: where t is the current number of iterations, is the best position in the generation VOLUME 8, 2020 and − → X (t) is the current position of the individual in the t generation.
− → A and − → C are determined by the following calculation: where a ∈ [0, 2] and gradually decrease as the number of iterations increases and r is a random number and r ∈ [0, 1].

2) STAGE OF BUBBLE ATTACK
In order to establish the mathematical model of humpback whale bubble attack, the following two methods are designed: The contractive encirclement is realized by reducing the value of a in equation (13). It can be seen from equation (13) 1], any position between the whale and its prey may become the next position where the whale swims.

b: SPIRAL POSITION UPDATING MECHANISM
First, this article calculates the distance between the whale's current position and its prey. Then, the spiral motion of humpback whales is simulated and a spiral equation between the two points is established. The humpback whale approaches its prey along the trajectory of spiral equation. The mathematical model is described as follows: where − → D is the distance between the position of the current individual and the position of the current optimal individual, b is the constant that defines the logarithmic spiral shape and l is a random number and l ∈ [−1, 1].
In addition, humpback whales will shrink around their prey along a spiral path. In order to simulate this behavior. It is assumed that the choice between the contraction encircling mechanism and the spiral mode is used to update the position of the whale and the probability of each behavior is 50%. The mathematical model is as follows: where p is a random number and p ∈ [0, 1].

3) STAGE OF FINDING FOOD
Similarly, variables are used to search for prey. In fact, humpback whales search randomly based on each other's location. It is assumed that, when |A| > 1, the humpback whale is forced to abandon the current best individual and update its position. This can enhance the global search ability of the algorithm. The mathematical model is expressed as follows: where −−→ X rand is a random individual selected from the current population.

B. IMPROVED FITNESS FUNCTION
At present, in PD positioning of TDOA method, one sensor will be selected as the reference sensor. For example, literatures [21] chose sensor S 1 as the reference sensor to solve the positioning equations. Assume that sensor 1 is used as a reference sensor. Since d i,1 are independently distributed and follow a Gaussian distribution with mean µ = d i − d 1 and variance δ 2 = v 2 δ 2 d , the position of the PD source P(x, y, z) can be estimated by the maximum likelihood estimation method. The likelihood function is described as follows: Calculating the maximum coordinate value of equation (20) is equivalent to the following: When the value of (21) is the smallest, the estimated coordinates of the position of the PD source are obtained. Since it is a complex nonlinear function, it is difficult to be solved directly. WOA is used to solve the problem. According to equation (21), the fitness function, with S 1 as the reference sensor, can be obtained as follows: where X is the position vector of the individual whale. (22) is sensitive to t 1 . If only S 1 is used as the reference sensor, the value of (22) will increase obviously when the error of t 1 is large. This will distort the fitness function and cannot reflect the quality of the solution exactly. Therefore, different sensors S i can be set as reference sensor to calculate the fitness value and reduce the occurrence of accidental situations. Assume that the reference sensor is S i , the fitness function f i (X ) is described as: where i = 1, 2, 3, · · · , N and i = j. According to equation (23), N fitness values can be calculated. Then the minimum value is substituted as the improved fitness function f (X ) into equation (21) to solve the optimal solution. The improved fitness function and the objective function are as follows:

C. IMPROVED INITIAL POPULATION BASED ON TIME REVERSAL IMADING
The distribution of the initial population has a great effect on convergence of the algorithm. Since there is no prior knowledge, most of the initial positions of intelligent optimization algorithms are generated randomly. This article employs the ideas of [22] and introduces results of the TR imaging into the solution process of WOA. The initial population of the WOA is generated based on the TR imaging. Fig. 3 is an imaging diagram in a two-dimensional plane using the TR method. It can also be viewed as a probability distribution map. The darker the color, the more likely it is to become the PD source. Then the roulette method is used to select the initial population.

D. WEIGHT SETTING BASED ON TIME REVERSAL IMAGING
The weighting factor can balance the global search ability and local optimization ability of the algorithm. Referring to the weight setting [23], combined with the TR imaging map, a weight setting strategy based on TR imaging results is proposed. The specific settings are as follows: Step 1: First, the fitness value f of each whale individual is calculated and arranged in order from small to large, then all fitness values are divided into two groups. The average fitness values are f avg1 and f avg2 and f avg1 < f avg2 . The fitness value of the current whale individual is compared with a and b respectively, and different inertia weights are set for each different whale individual.

1) f (X ) ≤ f avg1
The fitness of the current whale individual is better than the average fitness of the better group, indicating that the position of the individual may be relatively better in the population. Meanwhile the weight is taken ω 1 ∈ (0.8, 1.2). It is useful for excellent individuals searching in a small range to find the local optimal value.
The current individual fitness of whales is in the middle of the population, which belongs to general individuals. So, it can be made close to the optimal position according to the original path. At this case, assume ω 1 = 1.

3) f (X ) > f avg2
The current fitness of individual whales is at a poor level among the population, which belongs to unsatisfactory individuals. In this case, the weight ω 1 ∈ (0.3, 0.6) or ω 1 ∈ (1.3, 1.6) is taken and each of their probability is 50%. It will help poor individuals to search for other prey and enhance the global search capability.
Step 2: In order to solve the problem of large positioning error caused by time delay estimation error, the weight is set according to the positioning result combined with TR method. Referring to the probability diagram in Fig. 3, assume that the probability of the individual whale's position in the probability diagram is p, the weight can be set as ω 2 = p or ω 2 = 2 − p, each accounting for 50%.
Finally, let the weighting factors ω 1 and ω 2 account for half of the weight in the overall optimization process. The final weighting factor is shown as follows: Whale individuals can modify the optimization path according to the TR imaging results in the process of optimization, and reduce the impact of time delay error on the positioning results. By introducing the weighting factor, formula (17) can be rewritten as follows:

E. STEPS OF THE IMPROVED WHALE OPTIMIZATION ALGORITHM
The main steps of PD TDOA location using IWOA are as follows: (1) Initialization of population. Use the method in Section 2.3 to initialize the position of each individual whale. At the same time, the parameters are initialized, such as population number N , weight factor ω, logarithmic spiral shape constant b, random number l, initial iterations t and maximum iterations M .
(2) Calculate the initial fitness value. To calculate the fitness value of the initial whale individual and find the position X * of the best individual in the population from (24). (3) The first iteration. The first iteration is carried out according to the traditional WOA.
(4) Calculate fitness value and weight. After the iteration, the fitness value of the individual whale is calculated according to (24), and the weight is set according to the weight strategy in Section 2.4.  (4) and (5) until the maximum number of iterations is reached or the optimal position is found. The optimal individual position is output as the optimal solution of PD localization.
Flow chart of TR_IWOA for TODA positioning method is as follows:

IV. THE EXPERIMENT OF TIME REVERSAL IMAGING
In order to observe the TR imaging effect before and after BSS. the following experiments were designed.
The GIS model of the laboratory is used to carry out the double channel positioning experiment. Since PD location in GIS belongs to linear dual channel positioning model, it can complete positioning without positioning algorithm. The TR imaging is compared with various time delay extraction methods. The needle plate discharge model [24] made in the laboratory is used as the simulated PD source, as shown in Fig. 5. The model of the oscilloscope is lecroy8254, the bandwidth is 2.5GHz, and the maximum sampling rate is 20GSa/s. To eliminate the time delay error in the signal transmission process, the UHF sensor is directly connected to the oscilloscope through coaxial cables of the same length. The specific test platform is shown in Figure 6. The coordinate   BSS is performed using the method in [25]. The PD signal before and after BSS are shown in Fig. 7 and Fig. 8. TR imaging is shown in Fig 9. As shown in Fig. 9, the PD source is most likely to be in the red area of the TR imaging. It takes a period of time from the starting point of PD signal to the maximum amplitude point. Therefore, the maximum amplitude point in TR imaging is not the position of PD source. But the PD source must be in the red range with large amplitude. Table 1 shows the position of the center point and the range of red area in TR imaging. It can be seen from Table 1 that all the red areas of TR imaging cover the true position of PD source. It is found that blind source separation improves the focus of the TR imaging.    Table 2 shows the results of time delay extraction methods. In the PD position of GIS, the time delay is extracted first. Then the position of PD source is obtained by the following formula: where x is the distance between the PD source and sensor 1, L is the distance between the two sensors and v is the speed of the UHF signal with v = 3 × 10 8 m/s in GIS.
Next, compare the positioning results in Table 2 with the center point position of the TR imaging. The distance between the center point and PD source are less than the positioning errors in Table 2. Therefore, the following results are obtained: when PD source is located in GIS, the position of center point can be regarded as the location of PD source.

V. SIMULATION ANALYSIS
In order to verify the feasibility and effectiveness of IWOA algorithm in PD location, simulation is carried out in twodimensional plane. IWOA algorithm is compared with least square algorithm and particle swarm optimization algorithm.
First of all, a plane model of 400cm ×400cm is established and four sensors S 1 -S 4 are arranged. The coordinates of the sensors and PD are shown in Table 3. The theoretical time delay τ 12 , τ 13 and τ 14 can be achieved through formula (2) and (3). When the theoretical time delay is added to the  time delay errors τ 12 , τ 13 and τ 14 of different sizes respectively, the time delays τ 12 , τ 13 and τ 14 of the simulated actual measurement can be obtained. Take the UHF wave velocity in oil as v = 2.031 × 10 8 m/s [26]. The position of the center point of TR imaging is set randomly. Assume distance between the center point and PD source is r.

A. INFLUENCE OF DIFFERENT TIME DELAY ERRORS ON POSITIONING RESULTS
Assume r ≤ 15cm. Under the premise of known time delay τ 1i , sensor coordinate S i (x i , y i , z i ) and TR imaging results, the PD location is solved by least square method, PSO algorithm and IWOA algorithm respectively. The results are shown in Fig. 10. It can be seen that due to the influence and limitation of TR imaging, IWOA positioning results are closely around the simulated PD source. They are more concentrated than PSO positioning results. But the results of least square method are more scattered. The positioning errors under different time delay errors are shown in Fig. 11. It can be seen from the figure that:  (1) When τ 1i ≤ 0.5ns, the three methods can locate the PD source effectively.
(2) When 0.5ns < τ 1i ≤ 1ns, the positioning error of the least square method increases, so it can't locate effectively.
(3) When τ 1i > 1ns, the positioning error of PSO also increases gradually. The positioning error of IWOA fluctuates slightly with the change of time delay error. But its positioning error is always less than 10cm.
It is shown that the IWOA optimization algorithm can effectively locate the PD source. It has good stability and accuracy.

B. INFLUENCE OF DIFFERENT TIME REVERSAL IMAGINGS ON POSITIONING RESULTS
To analyze the influence of different TR imaging center point on PD location accuracy, different TR imaging are simulated by changing r. The distance between the imaging center point and the simulated PD source can be changed. Assume the distance between the TR imaging center and the simulated PD source( r) are 0 cm, 10 cm, 20 cm, 30 cm, 60 cm and 100 cm. Assume τ 1i are 0ns, 0.2ns, 0.4ns, 0.6ns, 0.8ns, and 1ns. r and τ 1i are matched one by one to analyze the positioning error. The specific positioning errors d is shown in Table 4.
It can be seen from the table: (1) When r < 30cm, the positioning results are accurate. It is shown that TR imaging plays a positive role in positioning when the time delay error is too large.
(2) When τ 1i < 0.6ns, even if r is too large, the positioning error is still less than 10cm.
(3) When r > 30cm and τ 1i > 0.6ns, with the increase of time delay error, the positioning error also increases gradually. And the positioning accuracy is already slightly worse than that of PSO.
In summary, when the accuracy of the TR imaging is high, the TR imaging map can play a positive guiding role in the process of solving the IWOA optimal solution. When the accuracy of the TR imaging is low, the advantages of the IWOA disappear.

VI. ANALYSIS OF LABORATORY RESULTS
In order to test the practical effect of IWOA optimization algorithm in PD positioning, PD positioning experiment is carried out in the 35kV transformer in the laboratory. The size of the transformer is 260cm ×220cm ×270cm. The needle plate discharge model and detection device in Section 3 are used. The specific test platform is shown in Fig. 12. The coordinates of four sensors were S1(220cm, 130cm, 10cm), S2(110cm, 260cm, 200cm), S3(0cm, 130cm, 10cm), and S4(110cm, 0cm, 150cm). Multiple measurements were obtained by changing the position of the PD source. The specific position of the PD source and the time delay measured by the energy accumulation method with S1 as the reference sensor are shown in Table 5. At the same time, TR imaging   of PD source at each location were obtained by TR method. Fig. 13 is a TR imaging of PD1.
Three methods of IWOA algorithm, PSO algorithm and least square method are used to solve the equation system. The positioning results are shown in Table 6.
It can be seen from the table 6: (1) when the delay error is small, the three methods can obtain the ideal positioning result. However, small changes in time delay will greatly change the positioning result of the least square method. It is shown that the positioning algorithms are very sensitive to delay errors. A slight increase in delay error may make the positioning algorithm unsolvable. For example, the least square method has no solution in the positioning process of pd3-pd5.
(2) The average positioning error of the PSO algorithm is 14.856cm. The overall situation is good. But the positioning error would increase with the increase of the delay error.
(3) Due to the restriction of TR imaging, the IWOA had low sensitivity to time delay error relatively. The average positioning error is 8.51cm, which meets the positioning accuracy requirements.
It is seen that the IWOA algorithm can reduce the sensitivity of the positioning algorithm to time delay errors. The reliability and accuracy of positioning are improved.

VII. CONCLUSION
To solve the problem that the TDOA method used in PD positioning is sensitive to delay errors, some researches are given as follows: (1) The TR method can ignore the influence of the waveform mode and does not need to estimate the time delay. It can show the energy distribution of PD signal arrival time. The general position of the PD source can be judged by TR imaging. It has a guiding effect on precise positioning of PD. It is found through experiments that BSS can not only improve the accuracy of time delay extraction, but also improve the focus of the TR imaging. VOLUME 8, 2020 (2) An improved whale optimization algorithm for PD location is proposed. The sensitivity of location algorithm to delay error is reduced. The average error of laboratory positioning results is 8.51cm. TR_IWOA has good positioning accuracy and stability.
(3) Since the PD signal is weak and the field interference is serious, the TR focusing is realized by data processing. In practical application, we can make TR imaging of PD signal obtained by oscilloscope in simulation scene. Or the PD signal is enhanced and transmitted from the position of the sensor receiving the signal, and TR imaging is performed in the actual scene. His research interest includes application of artificial intelligence in power systems. VOLUME 8, 2020