A New Vibration Analysis Approach for Monitoring the Working Condition of a High-Voltage Shunt Reactor

With the rapid development of power grids, the voltage level and capacity of high-voltage shunt reactors (HVSRs) are increasing year by year, and HVSR faults are increasing, especially HVSRs core winding faults. The HVRSs core winding faults seriously threaten the safe and stable operation of power grids. To solve this problem, we propose a feature extraction method for HVSR core winding faults identification to improve the accuracy of faults identification. This method relies on the combined application of qualitative analysis of phase space reconstruction (PSR) and quantitative calculation of an improved K-means clustering method. First, we employ PSR to qualitatively analyze the HVSR vibration signal and extract the phase trajectory feature quantity of the vibration signal. Then, this paper uses an improved K-means clustering optimized by an improved grasshopper optimization algorithm (GOA) to cluster the phase trajectory features, so as to obtain the cluster center coordinates. Furthermore, we employ the distance from the cluster center displacement vector to the origin and the angle change to perform quantitative calculations to realize fault identification. Finally, the fault simulation experimental data sets of 10kV and 20kV HVSRs verify the effectiveness of the proposed method. The experiment results show that the proposed method has better accuracy and can truly reflect the state characteristics of the HVSR core winding. The proposed high-accuracy method helps to improve the efficiency of on-site HVSRs condition assessment and maintenance.


I. INTRODUCTION
A high-voltage shunt reactor (HVSR) is a core component of power systems. An HVSR has variable functions, such as compensating reactive power, reducing line power loss, suppressing the increase in power frequency voltage, and preventing resonant voltage [1], [2]. The healthy operation of an HVSR is of great significance to sustain the stability of a power system and ensure efficient energy transmission [3]. However, most of the core winding mechanical failures and The associate editor coordinating the review of this manuscript and approving it for publication was Baoping Cai . electrical failures that cause abnormal HVSR vibrations, such as winding looseness, iron core looseness, internal components falling off, and interturn short circuits, may lead to a reduction in HVSR short-circuit resistance, local power failures or even large-scale power outages [4]. For example, a 220 kV overhead line trip caused by internal winding faults in an HVSR resulted in serious power failure in a large area and led to enormous economic losses [5]. In view of this threat, it is vital to devise an effective core winding failure feature extraction approach for core winding failure of HVSRs with abnormal vibrations. Such an approach is necessary to diagnose the core winding failure of HVSRs as soon as VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ possible. As a direct result, an early-warning mechanism is significant for guaranteeing the safe and reliable operation of HVSRs.
In recent years, scholars have performed many theoretical and practical studies on HVSR core winding failure. The main analysis methods for HVSR core winding failure include frequency response analysis [4], diagnosis gas analysis [6], sweep frequency measurement [7], and vibration signal analysis [8]. Among them, the vibration signal analysis method has been widely studied. Because it can sensitively reflect the operation state of a reactor, it is easy to realize live detection without direct electrical connection with the system. It should be noted that the reactor state indicates the condition of the reactor itself. Additionally, since the vibration signal during the operation of the reactor contains sufficient state information (it can effectively reflect the condition of the reactor), the vibration analysis method can accurately and effectively diagnose core winding faults in a reactor. It is worth mentioning that it is vital to diagnose an HVSR by vibration signals to extract considerable feature information.
To date, researchers have made contributions and achieved some progress in vibration signal fault feature extraction [9]- [11]. Vibration signal processing techniques for condition monitoring, fault diagnosis and prediction are classified into three dominant types: time domain, frequency domain and time frequency analysis. All of them are suitable for linear systems with stable signals, which are called traditional signal processing methods in this paper. However, once the above traditional methods are applied to nonlinear systems, there will be problems such as frequency leakage, endpoint effects and poor adaptability [12], [13]. It should be noted that the vibrational signal of an HVSR during operation presents strong time-varying and nonlinear characteristics. Therefore, traditional signal processing methods may not fully or veraciously capture the rich state characteristics contained in a vibration signal, where the state characteristics refer to the physical quantity that can effectively reflect the state information of the reactor. For confused nonlinear and nonstationary vibration signals, nonlinear time domain analysis methods (e.g., chaos and phase space reconstruction (PSR) technology) can completely extract the state information hidden in the vibration signal and effectively describe the essential characteristics and internal laws of the vibration signal. In [14], ergodic theory held that all condition information of signals can be extracted or recovered from the one-dimensional time series of signals. Takens et al. proposed an embedding theorem [15]. They proved that the phase space constructed by the time series retains the dynamic characteristics of the original time series. From this result, we can see that the change in the phase trajectory can reflect the relevant information of the original system. Yang et al. [16], [18] and Ruan et al. [17] effectively extracted the vibration signal characteristics at the moment of circuit breaker closing based on the chaos theory. Therefore, chaos technology can effectively extract the hidden state information in the vibration signal.
Although significant research on the analysis of vibration signals by chaos theory has been carried out, it has mainly focused on the qualitative analysis of phase trajectories. Meng et al. [19] employed phase trajectory features as training samples to predict the failure of rolling bearing vibration signals, but ignored the vector features and quantitative calculations of reconstructed signals, which may lead to misjudgment. Zhou et al. [20] used PSR to qualitatively analyze the vibration signal of the transformer, and judged whether the transformer winding is faulty or not based on the phase trajectory. However, the influence of the embedding dimension on the phase trajectory is ignored, and its diagnostic accuracy needs to be improved. Zhao et al. [21] and Li et al. [22] proposed the fault diagnosis model of vibration signal power spectrum and Markov chain, and used the phase space coefficient to identify the state characteristics of the vibration signal. But other sensitive features of the reconstructed signal in the high-dimensional space are not considered, and the recognition ability is limited. Huang and Gan [23] employed the percentage of the 90-degree angle and the non-90-degree angle of the phase trajectory point as the feature extraction standard, and extracted the vibration signal feature, which achieved a good fault classification. However, the time-delay and embedding dimension have a great influence on the determination of the angle of the phase trajectory point, and the reliability of the feature extraction result is insufficient. Fei [24] used PSR as an intermediate bridge to solve the problem of difficult extraction of weak feature signals, and constructed a vibration signal extraction method based on wavelet change-phase space reconstruction-singular value decomposition to accurately extract features. But the vector characteristic contained in the chaotic attractor is lost. In [25], a feature extraction method based on the combination of PSR and singular value decomposition was proposed, which extracted the state features of the vibration signal in different states of the circuit breaker. However, due to the large amount of calculation of the characteristic matrix formed by the phase trajectory, its calculation efficiency is low. In [26], a detection method of transformer short-circuit turns based on chaos theory-fractal dimension index was proposed, which combined with data mining method to realize the function of automatic diagnosis. The disadvantage is that the fractal dimension index has limited ability to classify fault types, and it is difficult to accurately identify the degree of fault.
Summarizing the related research on fault feature extraction based on chaos theory and vibration signal in recent years [19]- [26], we found that the current research mainly uses qualitative analysis of phase trajectory and chaotic index analysis for feature extraction and fault diagnosis. Research on the vector characteristics and quantitative calculations of reconstructed signals is limited. Regarding quantitative calculations, some scholars have proposed a performance indicator detection toolbox [27], [28], but the practicality needs further research. Therefore, in order to improve the pattern recognition ability and accurately warn the HVSR state change, this paper uses the clustering algorithm to quantitatively calculate and analyze the phase trajectory under the premise of fully considering the vector characteristics. We try to identify the state change of the HVSR via the distance from the cluster center displacement vector to the origin and the angle change. Cluster analysis classifies a series of unclassified objects based on the attributes of the objects and classifies a group of individuals into several categories according to similarity. It is worth noting that K-means clustering is widely used in processing numerical data. Here, considering that the characteristic quantities of phase trajectories change with the state of the reactor core and winding, we use cluster analysis to calculate the distribution of characteristic quantities of phase trajectories. This not only can effectively realize the state detection of the reactor core and winding, but also can improve the pattern recognition ability and quantitative analysis ability.
In light of the above considerations, we propose a fault diagnosis method combining PSR and an improved grasshopper optimization algorithm with K-means (IGOA-K-means) clustering to further enhance the accuracy of reactor fault diagnosis. This paper makes the following contributions: 1) To improve the accuracy of HVSR core winding faults diagnosis, we use improved K-means clustering to cluster the reconstructed phase trajectories. According to the cluster center space coordinates obtained by clustering, we calculate the distance from the cluster center to the origin as the first feature and the angle between the cluster center displacement vectors in the healthy state as the second feature value. The core winding state of HVSR is judged by the first feature value and the second feature value to realize fault warning and diagnosis.
2) We propose a method to optimize the original cluster centers of K-means clustering to overcome the defect of random selection of the K-means initial cluster center and improve the accuracy of the clustering results. We take the minimum distance from the cluster center as the fitness function and use the IGOA to optimize the selection of the initial cluster center of K-means clustering. Taking into account the problem of the adaptive difference of the linear parameters of the GOA algorithm, we use adaptive weight coefficients to optimize and improve the GOA algorithm.
3) We employ 10 kV and 20 kV HVSR fault simulation experiment cases to verify the effectiveness of the method proposed in this paper, summarize the clustering characteristics and change rules of different categories of faults, and discuss the precautions for field applications. Note that the HVSR faults simulated in this work are mainly core winding faults that cause abnormal vibration signals.
The remainder of this paper is organized as follows. Section II presents the PSR method of HVSR vibration signals and the proposed IGOA-K-means clustering. Section III carries out two case studies for the evaluation of the proposed method using experimental datasets. A discussion is provided in Section IV. Finally, Section V concludes the present work.

II. METHODS OF ANALYSIS
This section first uses the small data set method to judge the chaotic characteristics of the reactor vibration signal. Then the embedding dimension and time delay of the reconstructed phase space of the vibration signal are determined. Subsequently, a method for clustering analysis of phase trajectories using the optimized K-means algorithm is proposed. Finally, according to the obtained cluster center coordinates, we calculate the distance from the cluster center displacement vector to the origin and the change of the included angle.

A. RECONSTRUCTION OF THE PHASE PORTRAIT FROM THE MEASURED VIBRATIONS
An HVSR is a very complex nonlinear dynamic system. In the process of operation, the vibration signal generated by an HVSR not only includes the dynamic characteristics of real-time output but also covers the information of any variable and the state evolution of the HVSR system.

1) CHAOTIC CHARACTERISTICS OF HVSR VIBRATION SIGNALS
Before using PSR theory to analyze HVSR vibration signals, it is necessary to judge whether the HVSR vibration signal has chaotic characteristics. The basic characteristic of chaotic motion is that the motion is very sensitive to the initial conditions. The orbits generated by two extremely close initial values are separated exponentially over time. This phenomenon can be quantitatively described by the Lyapunov exponent [29]. The Lyapunov exponent (LE) quantitatively describes the divergence rate of the nearby orbit. In the course of movement, the distance between two adjacent phase tracks is δx(0) at the initial time. After time t, the maximum component of the distance is as follows: where λ is the Lyapunov exponent.
In chaos research and practical engineering applications, it is not necessary to compute all LEs of the time series, but it is sufficient to calculate the largest Lyapunov exponent (LLE) that determines the orbital divergence rate. A value of LLE greater than zero means that the vibration signals have chaotic characteristics [29]. The methods currently used to compute the LLE mainly include the physical definition method, Wolf algorithm, Jocobian method, P-norm method and small data set method. Considering the characteristics of a small data group, high reliability, small amount of calculation and easy operation [29], this paper uses the small data set method to calculate the LLE for HVSR vibration signals.

2) RECONSTRUCTION OF HVSR VIBRATION SIGNAL PHASE SPACE
Regarding the dynamic factor analysis of nonlinear time series, the delayed coordinate state space reconstruction method is widely used at present, which is PSR technology. The phase trajectory diagram obtained by PSR technology can effectively predict and judge the dynamic characteristics and state characteristics of a dynamic system. Packard et al. suggested reconstructing the phase space with the delay coordinates of a certain variable in the initial system [30]. Takens [15] proved that a favorable embedding dimension could be found, that is, if the dimension of the delay coordinate m ≥ 2d + 1 (d is the dimension of the dynamic system), the regular trajectory (attractor) can be recovered in this embedding dimension space. On the trajectory in the reconstructed m-dimensional space, the motive force system maintains diffeomorphism. Suppose the measured time series of the one-dimensional vibration signal of the HVSR is: where t 0 is the starting time point of reactor vibration signal monitoring; t is the time interval; and N is the length of the time sequence. According to the one-dimensional time series, the phase distribution of an m-dimensional phase space can be reconstructed as follows: where ζ = N − (m − 1)τ , m is the embedding dimension, and τ is the time delay. These ζ vectors constitute the phase space of the reconstructed signal of the HVSR system. The suitable embedding dimension and time delay are significant to PSR. According to Takens theorem [15], for an ideal infinite and noise-free one-dimensional time series, the embedding dimension m and the time delay τ can take arbitrary values. However, the actual time series are usually limited in length and there is noise, so appropriate m and τ should be calculated to achieve accurate reconstruction of the phase space.

3) DETERMINATION OF EMBEDDING DIMENSION
The embedding dimension refers to the minimum spatial dimension that completely covers all the state information of the reconstructed signal of the reactor. When the dimension of the reconstruction space reaches or exceeds the embedding dimension, the reconstructed signal can be fully expanded to eliminate the self-intersection phenomenon [31]. Therefore, proper selection of the embedding dimension may thoroughly reveal the dynamic parameters of the HVSR system. The commonly used methods include the G-P algorithm proposed by Grassberger and Procaccia and the FNN (false nearest neighbor) algorithm [31], [32]. Here, we employ the G-P algorithm to determine the optimal embedding dimension. The G-P algorithm mainly includes the following points: (1) set a relatively small reconstruction embedding dimension and reconstruct the reactor vibration signals of formula (2) to obtain a set of reconstructed phase spaces as shown in formula (3); (2) calculate the Euclidean distance r ij (m) of any pair of phase points in the phase space obtained from the reconstruction of the reactor vibration signals: (3) compute the correlation function C(r, m): where θ(x) is the Heaviside function, expressed as formula (6); C(r, m) denotes a cumulative distribution function, representing the probability that the distance between two phase points in the reconstructed phase space is less than r; within an appropriate range of r, the correlation dimension d(m) of the HVSR reconstruction signal and the cumulative distribution function C(r, m) should satisfy a linear relationship, namely, d(m) = lnC(r, m)/ln (r), so that the correlation dimension estimate d(m 0 ) corresponding to m 0 can be obtained by fitting; (4) increase the embedding dimension so that m 1 > m 0 , and repeat (3) until the corresponding dimension estimate d(m) no longer increases with the increase of m and remains unchanged within a certain error range. At this time, the d(m) obtained is the correlation dimension of the number m, and the corresponding m is the optimal embedding dimension required.

4) CALCULATION TIME DELAY
The methods used to calculate the time delay currently include the autocorrelation function method and mutual information method (MIM) [16], [33]. However, the autocorrelation function is based on the correlation calculation of linear theory, which is suitable for the time delay calculation of linear signals. According to Shannon information entropy, the MIM uses the mutual information value (MIV) to measure the degree of random correlation between two random variables, which is suitable for nonlinear signals. Taking into account the nonlinear characteristics of the HVSR vibration signals, we employ the MIM to obtain the best time-delay of the HVSR power system. Suppose the original vibration signal of the HVSR is x ={x(i),1,2, . . . , N }, and the vibration signal after the timedelay is rewritten as x τ ={x(i + τ ),1,2,. . . , N }. Then, relying on the definition of Shannon information entropy and probability distribution [33], the mutual information of the vibration signal is as follows: 46490 VOLUME 9, 2021 where H (x) and H (x τ ) denote, respectively, the information entropy of x and x τ , H (x, x τ ) is the joint information entropy of x and x τ , P(x i ) and P(x i+τ ) show, respectively, the probability of occurrence of x i and x i+τ , and P(x i , x i+τ ) denotes the joint probability distribution of x i and x i+τ .

B. VIBRATION SIGNAL ANALYSIS BASED ON IGOA-K-MEANS
After the HVSR vibration signals are reconstructed, a mass of phase points distributed in the high-dimensional phase space are formed as a phase trajectory diagram. Different mechanical states of the HVSR correspond to different phase trajectories. At present, it is difficult to construct a mathematical model of the change in the mechanical characteristics of the reactor. Here, we start from the phase trajectory diagram of the HVSR vibration signals to analyze the dynamic characteristics of the HVSR under different mechanical states and different fault conditions. IGOA-K-means clustering is used to quantitatively describe and identify the distribution characteristics of the reconstructed high-dimensional phase trajectory map.

1) K-MEANS
The basic concept of K-means clustering is to divide a large number of high-dimensional data points into k clusters according to the characteristics of the data and extract the most representative points of each cluster as the prototype of the data. Then, the prototype points of the k clusters can represent the repetitive all-status information of the structure signals [34], [35]. We use K-means to calculate that the HVSR vibration signals are reconstructed into a high-dimensional phase point, as shown in formula (3). The specific process contains the following points: 1) Determine the initial number of cluster centers k and the initial distribution of cluster centers; 2) According to the principle of closest distance, the reconstructed phase points are allocated to the corresponding cluster centers to form new k clusters; 3) Compute the cluster centers of the newly generated k clusters; 4) Repeat 2) and 3) until the position of the cluster center is stable and classify the phase points; 5) Compute the distance J (C k ) between the center of each cluster and the phase point in the cluster, and sum the distance between each cluster center to obtain the overall distance B, which is: where σ k denotes the cluster center and K is the total number of cluster centers.

2) GRASSHOPPER OPTIMIZATION ALGORITHM (GOA)
We employ the GOA to optimize the selection of K-means initial cluster centers to solve the problem of random selection of K-means initial cluster centers and improve the accuracy of the clustering results. The GOA was proposed by Shahrzad Saremi in 2017, and its basic concept was inspired by the foraging behavior of grasshopper groups [36]. The GOA, like most intelligent algorithms, has two procedures of exploration and development, which guarantees that the GOA has a robust global search capability and can avoid local optimization. Here, we regard each phase point of the reconstructed vibration signals of the HVSR as a grasshopper. The implementation steps of employing GOA to select the initial cluster center are as follows: (i) Initialize the population position P i , the number of grasshoppers M , the parameters c min , c max and the maximum number of iterations L; where S i defines the social interaction between individual grasshoppers; G i is the gravity received by the i-th grasshopper; A i denotes the wind convection received by the i-th grasshopper; p i shows the i-th grasshopper; s is the social force; when s > 0, the grasshoppers attract each other; when s < 0, the grasshoppers repel each other; g denotes the gravity constant;ê g is the unit vector pointing to the center of the earth; d ij = | p j − p i | shows the distance between the i-th grasshopper and the j-th grasshopper; u denotes the wind constant; andê w defines the unit vector pointing to the wind direction.
(ii) Compute the fitness value of each grasshopper in the population and save the grasshopper positionT with the best fitness value; (iiii) Employ the following formula (11) to update c: where c max denotes the maximum value, c min defines the minimum value, l is the current iteration, and L shows the maximum number of iterations. We assign 1, 0.0001 and 100 to c max , c min and L, respectively. (iv) Update the grasshopper position P d i according to formula (12), compute the fitness value of each grasshopper, and update the optimal grasshopper positionT d in the d-dimensional space; where c is the decreasing coefficient, which affects the range of comfort zone, repulsion zone and attraction zone; ub d VOLUME 9, 2021 and lb d indicate, respectively, the upper and lower bounds of grasshopper in the d-dimensional space.
(v) If the fitness value of the grasshopper obtained in (iv) is better than the previous fitness value, the optimal position of the grasshopper is updated; (vi) Repeat the processes (iii)-(v) to judge whether the maximum number of iterations is reached; if the maximum number of iterations is reached, the algorithm ends, and the optimal position valueT d is output.

3) THE IMPROVED GOA (IGOA)
The parameter settings in the standard GOA algorithm directly affect the performance of the GOA, and is very important to the function of the algorithm [37], [38]. We employ the adaptive inertial weight w to improve formula (12) to better balance the global exploration capability and local development capability of the GOA. Expression (12), after introducing w, is shown in the following formula (13): In the early stage of GOA, the grasshoppers are scattered to a relatively high degree. We employ a larger w to make the grasshoppers have a strong global exploration ability and quickly find a better position in a large area. In the late stage of the GOA, the grasshoppers gather at a higher level, and a smaller w is helpful for the grasshoppers to search in a local area. The specific improvement strategies are as follows.
Assuming that the fitness value of the individual grasshopper in the t-th iteration is f pi(t) , we calculate the individual fitness value of the grasshopper, arrange the values in ascending order and divide them into two parts, and find the average fitness values f avg1 and f avg2 of the two parts (where f avg1 < f avg2 ). We compare the fitness value f pi(t) of each grasshopper with f avg1 and f avg2 and divide the grasshopper population into three levels, each with different inertia weights.
(a) f pi(t) < f avg1 : the grasshoppers at this level are farther from the optimal grasshopper positionT d in the population. At this time, w should take a larger inertia weight value to improve the global exploration capability, so we take w = 0.9 in this paper. (b) f pi(t) < f avg2 : the position of the grasshoppers in the hierarchy is closer to the optimal grasshopper positionT d , and w should take a smaller inertia weight value, which is conducive to the local development and convergence of the algorithm. In this case, we take w = 0.2. (c) f avg1 < f pi(t) < f avg2 : the grasshoppers belong to ordinary grasshoppers at this level. At this time, the random inertia weighting strategy is adopted, and the value of w is randomly selected in [0.4, 0.6]. This random selection method helps to improve the optimization ability of the algorithm.

4) IGOA-K-MEANS CLUSTER ANALYSIS OF THE PHASE TRAJECTORY OF HVSR VIBRATION SIGNALS
We employ the IGOA to optimize the initial K-means clustering centers and construct the objective function, namely, the fitness function, as shown in formula (14).
The flow chart of clustering and identifying phase trajectory features using the IGOA-K-means algorithm is shown in Fig. 1. According to Fig. 1, we calculate the overall distance J (C) when K is 1-15. In general, as the number of cluster centers K increases, the overall distance J (C) will continue to decrease [35]. This paper takes the K value when the overall distance J (C) decreases and slows down as the number of cluster centers. At the same time, it is worth noting that, relying on the calculation process in Fig. 1, we can also obtain the cluster center coordinates. In this paper, the change in the internal state of the HVSR is quantitatively identified based on the distance from the cluster center displacement vector to the origin and the change in the included angle. Here, to promote the understanding of the abovementioned HVSR chaotic analysis phase trajectory feature extraction and IGOA-K-means clustering identification, we summarize the detailed process as follows: Step 1) The collected raw vibration signal is preprocessed for noise reduction and detrending.
Step 2) Choose the best measurement point according to the vibration intensity of the measurement point.
Step 3) Compute the value of LLE to determine whether the vibration signal has chaotic characteristics.
Step 4) Employ the G-P algorithm to calculate the best embedding dimension m and the mutual information method to select the appropriate delay time τ .
Step 5) The one-dimensional vibration signal of HVSR is reconstructed into the multidimensional phase space, and the phase trajectory characteristic quantity of the vibration signal is recovered.
Step 6) Select the number of cluster centers and use the IGOA-K-means clustering algorithm proposed in this paper to cluster the phase trajectory features and calculate the cluster center coordinates.
Step 7) The first feature value and the second feature value are solved by using the cluster center coordinates. According to the first feature value and the second feature value, the core winding fault of the HVSR is classified and identified.

III. CASE STUDY
This section introduces two case studies to test the sensitivity and stability of the proposed diagnosis method for HVSR internal faults. Case I uses the data collected in a 10 kV HVSR fault simulation experiment system to evaluate the performance of the proposed method. These datasets include three kinds of typical internal HVSR defects. Case II performs a fault simulation on a 20 kV HVSR, whose structure is different from the 10 kV HVSR, and uses the collected data to further assess the performance of the method proposed in this paper. Case II simulates three different levels of internal faults as well as healthy states.
A. CASE I: 10-kV, JSRT-ACL11H-30 HVSR This case employs a single-phase oil-immersed HVSR tank surface vibration signal to verify the proposed method. The rated voltage of the HVSR is 10 kV, the rated capacity is 30 kvar, and the rated frequency is 50 Hz. It should be noted that a customized HVSR is used in this case. The difference from an ordinary reactor is that the customized reactor in this case adds three insulated terminals, which are connected to the winding taps from different layers of the internal windings. The five terminals of the reactor are numbered 1, 2, 3, 4, and 5, which are connected to the first turn of the winding, which is the beginning of the winding, the 1017th turn, the 1026th turn, the 1035th turn, and the 1044th turn, which is the end of the winding.
This case established a 10 kV reactor fault simulation experiment platform, as shown in Fig. 2. This highperformance data acquisition platform was mainly composed of seven parts: 1) an AC power supply, which provides power for the HVSR and vibration signal acquisition system; 2) a voltage regulator, which increases the laboratory voltage to the rated voltage of the low-voltage side of the transformer; 3) a current limiting resistor, used to limit the inrush current at startup; 4) a transformer (low-voltage side: 380 V, high-voltage side: 10 kV), used for boosting, which increases the 380 V voltage to 10 kV to reach the rated working voltage of the HVSR; 5) an HVSR, for fault simulation; 6) a capacitor, which provides capacitive reactive power to offset the inductive reactive power of the HVSR and keeps the power supply system stable; 7) a vibration signal acquisition system, used for data acquisition of vibration signals, including highsensitivity vibration acceleration sensors, a DH5922D vibration signal acquisition instrument, data transmission lines and a computer.
During the operation of the HVSR, the vibration acceleration sensor was used to obtain the surface vibration signal of the HVSR tank with a sensitivity of 500 mV/g, a range of −10 g-10 g, an installation resonance frequency of >17 kHz, an operating temperature of −50-80 • , and a maximum lateral sensitivity of 5%. Taking into account the distribution range of the vibration frequency components of the HVSR, which are mainly distributed from 1-2000 Hz, we set the sampling frequency of the acquisition system to 20 kHz. The sampling time of each record was kept at 5-10 s. Aiming to collect the vibration signal of the HVSR more efficiently and accurately, we employed the array layout method to layout the acceleration measurement points, as shown in Fig. 3. On the four side surfaces, we set up a total of forty-eight measurement points in three rows and sixteen columns, numbered 1-48; on the top surface, we set up five measurement points, numbered T1-T5.
Three categories of representative fault defects were artificially simulated in the HVSR: 1) A core and winding axial loosening failure (CWALF) is generated by adjusting the axial tightening bolt, as shown in Fig. 4(a). During the transportation and installation of the HVSR, the iron core and winding may be loose. We simulated two kinds of potential looseness faults:  a light-axial-loosening fault of the core and winding (LALFCW), which applies an 80% axial preload, and a serious-axial-loosening fault of the core and winding (SALFCW), which applies a 50% axial preload under rated voltage. 2) The internal component drop malfunction (ICDM) is simulated by breaking the fastening strap of the winding, as shown in Fig. 4(b). The internal components of the reactor were aging for long-term operation, and when the short-circuit current was impacted, the aging parts might fall off (such as winding spacers). Without affecting the normal operation of the reactor, we simulated the failure of the HVSR internal components under the rated voltage. 3) A short-circuit fault between winding turns (SCFWT), short-circuit insulated terminal 2 and insulated terminal 3, that is, short-circuit 9 turns, to simulate short-circuit fault between turns, as shown in Fig. 4(c). The aging of winding insulation may lead to short-circuit faults between turns. Taking into account the safety of the experiment and the service life of the HVSR, the vibration signal acquisition time of the interturn short circuit fault condition was controlled within 1.5-2.5 s.
The above three types of fault defects were built to simulate actual situations.
Aiming to obtain vibration signals with a better signalto-noise ratio and more accurately obtain the features of the vibration signal of the reactor in different states, we select the best measurement point on the side. Taking into account the different strengths of vibration signals collected at different positions on the surface of the reactor box, we select the maximum vibration data of each measurement point under the healthy condition of the reactor for comparative analysis. We take the measurement point with the highest vibration intensity among all measurement points as the best measurement point. Fig. 5 exhibits the distribution of the maximum vibration intensity at different measurement points. It can be seen in Fig. 5 that the measurement point with the largest vibration amplitude is the third point in the 14th column, that is, measurement point 42, and its corresponding HVSR surface position is shown in Fig. 3(c). In this case, we employ the HVSR vibration signal data collected at measurement point 42 and top surface T1-T5 measurement points for analysis. The experiments were carried out in the spring from April to May at a room temperature of approximately 20 • C. After denoising and preprocessing, the vibration signals of the HVSR at the measurement points T1 and 42 under the three kinds of simulated defects are shown in Fig. 6.

1) RECONSTRUCTION OF THE HVSR VIBRATION SIGNALS
We employ the small data amount method to judge the chaotic characteristics of the HVSR vibration signals. This paper computes the LLE of the HVSR healthy/faulty vibration signals. Fig. 7 and Fig. 8 illustrate the LE spectra of healthy/faulty states of measurement points T1 and 42, respectively. We employ the least square method to fit the LE spectrum and obtain the red imaginary line in the graph, where the slope of the fitted line is LLE.  Considering the readability and length of the paper, this paper lists the LLE results of the vibration signals at the top measurement points T2-T4 in Table 1. Because the value of LLE is positive in both healthy and faulty conditions, we can determine the vibration signals of the HVSR chaotic characteristics. Now, according to the theory of Section II-A-2, we reconstruct the HVSR vibration signals into a multidimensional  phase space. Two significant parameters, namely, the time delay and embedding dimension, should be determined before PSR. Fig. 9 shows the mutual information value curve of the vibration signals at the measurement points T1 and 42 in the healthy/faulty states. According to the MIV curve, the first minimum value of the MIV is selected as the time delay value. It can be seen in Fig. 9 that the time delay values of the same measurement point in different states are different. The time delay calculation results of different measurement points are listed in Table 2. Relying on the vibration signals collected from measurement points T1 and 42, this paper calculates the curve of ln C(r, m) with ln r, as shown in Fig. 10, where the embedding dimension m corresponding to the actual vibration signals increases from 1 to 9. In this article, the best embedding VOLUME 9, 2021 dimension m is the value corresponding to when the linear part of the curve no longer changes. Fig. 10 shows that when m is 3, the linear part of the curve basically no longer changes, so 3 is chosen as the best embedding dimension. It should be noted that the best embedding dimension of the measurement points T1 and 42 in different states is 3, and the best embedding dimensions of measurement points T2-T4 is also 3 in this case. Therefore, the embedding dimensions of different measurement points are all taken as 3. After determining the time delay and the embedding dimension, the HVSR vibration signals are reconstructed in the high-dimensional phase space according to II-A-2) reconstruction theory. Fig. 11 indicates the reconstructed phase space results of the vibration signal in the healthy/faulty states. It is worth noting that the phase trajectories of different faulty states and healthy states have changed significantly, as shown in Fig. 11. We take the vibration signals of the measurement point T1 as an example for simple analysis. In a healthy state, the vibration signal phase trajectory condenses into an irregular cluster. In the faulty conditions of the axial loosening of the iron core and winding, the phase trajectory of the vibration signals opens in an irregular ring structure, and the more severe the looseness of the iron core and winding, the greater the degree of opening of the irregular ring structure of the phase trajectory. In the state of internal component drop failure and interturn short-circuit failure, the phase trajectories are also irregular ring structures, but the phase trajectories are obviously different. Certainly, this phenomenon also exists in the phase trajectory of the measurement point 42 in different fault states and the phase trajectories of different states are different. It can be seen that the change in the phase trajectory can effectively and qualitatively analyze the state change of the HVSR. Thus, we employ the phase trajectory of the vibration signal to extract the characteristic quantities of the HVSR different states.

2) IGOA-K-MEANS CLUSTER ANALYSIS
Aiming to identify the faults more accurately, we employ the improved IGOA-K-means clustering algorithm in this paper to quantitatively analyze the phase trajectory characteristics of the reconstructed signals in different states of the HVSR. Before cluster analysis, the number K of cluster centers needs to be determined first. This paper chooses the number of cluster centers when the J (C) drop rate is below 20% and remains basically stable to ensure the accuracy of clustering. As shown in Fig. 12, we compute the J (C) change curve of the cluster center of the reconstructed signal phase trajectories of the measurement points T1 and 42 as the value of K increases according to the overall distance formula (9) of the cluster center.   12 depicts the variation curve of the total distance between cluster centers as the number of cluster centers increases. It can be seen in Fig. 12(a) that the total distance J (C) of cluster centers decreases as the number of cluster centers increases. When K was 7, the decrease rate of J (C) at the T1 measurement point decreased from the initial 35.66%-71.25% to 11.60%-18.81% and remained stable. As shown in Fig. 12(b), when K was 6, the drop rate of J (C) at the measurement point 42 decreased from the initial 43.39%-71.39% to 14.81%-19.56% and remained relatively stable. Thus, the number K of the cluster centers at the measurement point T1 is 7 and the K value of the measurement point 42 is 6. In this case, considering that the vibration signal of the same surface has the same change rule, we take the K value of the top surface T2-T4 to be 7.
Here, we employ the IGOA-K-means clustering method proposed in this paper to calculate the position of the cluster center in the phase space under the healthy/faulty states. The calculation results are shown in Fig. 13. Fig. 13 illustrates the position distribution of the cluster centers of the measurement points T1-T5 and 42 in the phase space in the healthy/faulty states. It can be seen in Fig. 13 that after a fault occurs inside the HVSR, the cluster center position is significantly shifted from the healthy state. At the same time, the center offsets of different fault types are different, and this phenomenon exists at each measurement point. The reason for this phenomenon may be that the internal magnetic field distribution of the HVSR and the force of the internal structure of the HVSR change after the internal failure, which affects the transmission capacity of the HVSR vibration signal. For example, the latent failure of the HVSR core and winding axial looseness will cause the internal electromagnetic field distribution and natural frequency of the HVSR to change during operation [39], [40]. The short-circuit fault between the windings of the HVSR will also cause the internal electromagnetic field to change, which will affect the vibration signal transmission [41]. It can be seen that the cluster center distribution of the reconstructed signal of the HVSR could effectively reflect the change of the internal structure of the HVSR.
This paper proposes to employ the distance between the cluster center displacement vector and the origin and the angle between the cluster center displacement vector in the fault states and the health state cluster center displacement vector to calculate the cluster center offset in the phase space to quantitatively judge the condition change of the HVSR. According to the cluster center coordinates of each measurement point in the healthy/faulty states in Fig. 13, the vector summation operation is performed on the cluster center coordinates to extract the coordinates of the displacement vector. Depending on the coordinates of the displacement vector of the cluster centers, we can calculate the distance from the origin and the angle of the displacement vector between the faulty states and the healthy state.
Note that for the sake of comparison, we call the distance from the cluster center displacement vector to the origin the first feature value in this paper. Taking the cluster center displacement vector in the healthy state as the benchmark, we call the angle formed by the cluster center displacement vector in the faulty states and the cluster center displacement vector in the healthy state the second feature value. In particular, the second feature value represents the degree of difference between the vibration mode in the current state and the vibration mode in the normal state. The first feature value and the second feature value of the cluster in Fig. 13 are shown in Fig. 14 and Fig. 15, respectively. The first feature values of different measurement points in the healthy/faulty state are presented in Fig. 14. It is worth mentioning that the first feature value of each measurement point not only can quantitatively determine the change in the internal state of the HVSR but also can identify different fault types. It can be seen in Fig. 14 that there is a corresponding relationship between the change of the first feature value of the IGOA-K-means cluster and the state change of the HVSR.   15 shows the second feature value curve of the faulty states at different measurement points. In particular, the larger the second feature value is, the greater the difference between the vibration mode in the faulty state and the vibration mode in the healthy state. It can be seen in Fig. 15 that the second feature value of different faulty states has changed greatly, and the change of the internal state of the HVSR may be quantitatively determined.
Specifically, in this case, the summary and analysis of the change law of the first feature value and the second feature value of IGOA-K-means clustering in different states of the HVSR are as follows: (1) When he core and winding axial loosening fault occurs, the first feature value of different measurement points obviously increases, and as the loosening fault strengthens, the first characteristic value gradually increases. In addition, the second feature values of different measurement points have a positive correlation with the degree of looseness of the iron core and windings. The reasons for the above phenomenon may be as follows: 1) when the iron core becomes loose, the leakage magnetic field between the gap joints of the silicon steel sheets and the layers increases, the electromagnetic force between the silicon steel sheets increases, and the axial vibration generated increases; 2) after the winding is loosened, the compression force on the coil is reduced, and the electromotive force causes the acceleration of the coil to increase, causing the vibration to be further strengthened; and 3) there is a certain degree of coupling between the vibration caused by the loosening of the iron core and the winding. (2) When an internal component fails, the first feature value of different measurement points is between the first feature value of LALFCW and the first feature value of SALFCW. The second feature values of different measurement points have increased. The main reason may be that the component simulated in this case falls off, which leads to an increase in the electric power received by the winding in the radial direction of the HVSR during operation, which causes an increase in the total vibration amplitude. (3) When a short-circuit fault exists between winding turns, the first feature values of the measurement points on the top surface are close to the first feature value when it is a SALFCW. The reason may be that after a short-circuit fault occurs between the winding turns, the internal magnetic field strength of the HVSR increases, and the electromagnetic force generated by the silicon steel sheet increases, which causes the axial vibration to also increase. However, the first feature value of the measurement point 42 is greater than the first feature value in the SALFCW fault state. This is mainly because the current flowing through the winding increases rapidly after the winding is short-circuited between turns, and the resulting radial electromotive force is sharply increased, resulting in a substantial increase in the vibration intensity received by the tank wall on the side of the HVSR. Compared with the healthy state, the second feature values of different measurement points in the interturn short circuit are greatly increased. The analysis results show that the subtle changes of the HVSR vibration signal may be clearly expressed by the first feature value and the second feature value proposed in this paper. Therefore, this also indicates that the proposed IGOA-K-means clustering method is effective in identifying the phase trajectory characteristic quantities of the HVSR vibration signals.
B. CASE II: 20 kV, BKD-6700/20 HVSR This case further verifies the effectiveness and versatility of the proposed HVSR vibration feature extraction and fault identification method. In this case, a 20 kV, 6700 kvar HVSR is used for vibration signal collection, as shown in Fig. 16(a) for the field collection diagram. Taking into account the economy and reusability of the HVSR, we avoided destructive and unrecoverable failure experiments. We simulated three typical latent failures of the core and winding loosening: a 60% looseness of the core and winding (60% LCW), 100% looseness of the core and winding (100% LCW), and a complete loosening of winding (CLW). Of course, vibration signal acquisition of the healthy state (HS) was also carried out. It should be noted that the LCW was set by an axial fastening bolt, and the CLW was realized by complete dislocation of the top pressing pad under a marked preload, as shown in Fig. 16(b). The vibration acceleration sensors are arranged on the surface of the HVSR tank, and the specific layout of each side and top surface is shown in Fig. 16 (c)-(f). The vibration signal data acquisition system is the same as Case I, and the sampling frequency was set to 20 kHz. The sampling time of each record was 5-10 s.
Here, we first perform noise reduction and detrend preprocessing on the collected original vibration signals, according to the steps described in Section II B.4., and then select the best measurement point. The maximum intensity distribution of the vibration signals at each measurement point on the side is illustrated in Fig. 17. It can be seen in Fig. 17 that the best measurement point for the vibration signals is the measurement point in the third row and the fifth column, that is, the measurement point 23, as shown in Fig. 16(e). Fig. 18 shows the time-domain diagram of the vibration signals under different loosening faulty states at the measurement point 23. Here, comparing the waveform of the vibration signals, we can conclude that the shape of the healthy/loose state is similar, except that the amplitude of the iron core and winding is completely loosened, which makes the fault information difficult to find.  Furthermore, we employ the steps introduced in Case I to identify whether the vibration signal has chaotic characteristics. Fig. 19 shows the LE evolution curve of the vibration signals at measurement point 23, and the LLE value is obtained by fitting. Table 3 shows the calculation results of LLE at D1-D5 on the top surface. Here, from the calculation results in Fig. 19 and Table 3, we can see that the LLE of the vibration signal at different measurement points are all positive. Therefore, the vibration signals of HVSR in the healthy/faulty states have chaotic characteristics.      The reconstructed trajectories of vibration signals in the healthy/faulty states of the 20 kV HVSR at the measurement point 23 are shown in Fig. 22. It is not difficult to find that the phase trajectories of different states have changed significantly, and each phase trajectory has its own characteristics. Under healthy conditions, the phase trajectory presents an irregular aggregate ring. When the iron core and the winding are loosened, the phase trajectory opens in an irregular ring shape, and as the looseness of the iron core and the winding increases, the ring opening degree gradually increases. As described in Case I, IGOA-K-means clustering is used to identify phase trajectory features. Following the steps of Case I, here, we first select the appropriate cluster centers K . Fig. 23 shows the curve of the total distance J (C) from the cluster center of the reconstructed signal phase trajectory of the measurement points 23 and D1. It can be seen in Fig. 23 that the overall distance J (C) of cluster centers reduces with the increase of the number of cluster centers.
According to the selection method of the K value in Case I: for measurement point 23, when K ≥ 6, the drop rate remains below 20% and remains relatively stable.
For measurement point D1, when K ≥ 7, the drop rate remains below 20% and remains relatively stable. Therefore, the K value of measurement point 23 is taken as 6, and the K value of D1 is taken as 7. At the same time, the K value of other measurement points on the top surface of the HVSR is also taken as 7.  It is not difficult to find that the iron core and windings are loosened, the position of the cluster center in the phase space is significantly shifted from the healthy state, and the more severe the loosening, the greater the shift. Based on the cluster center coordinates in the above four states, the distance from the cluster center displacement vector to the origin and the included angle with the cluster center displacement vector in the healthy state are calculated. It is worth noting that the first feature value and the second feature value can reflect the deviation of the total position of the cluster center and effectively identify the change of the internal state of the HVSR. Fig. 25 shows the change in the first feature values and the second feature values under healthy/faulty conditions. Here, we find that the first feature value is positively correlated with the degree of core and winding looseness failure. Moreover, in the core and winding loose faults, the second feature value is also positively correlated with the degree of looseness of the core and windings. Therefore, the changes in the first feature value and the second feature value proposed in this paper can clearly reflect the internal state of the HVSR, and the classification performance is quite good. It is noteworthy that this method can accurately classify and recognize HVSR internal faults without machine learning or a large number of training samples and has the advantage of high sensitivity to identify faults.

A. COMPARISON OF CLUSTERING METHODS
Here, to verify the superiority of the IGOA-K-means clustering proposed in this paper, we employ the phase trajectory data extracted in Case I with rich vibration mode information to compare the performance of the traditional K-means, GOA-K-means and IGOA-K-means clustering methods. We mainly compare and analyze the first feature value and the second feature value under the healthy/faulty states obtained from the measurement points T1-T5 and 42 in Case I. We first calculate the coordinates of the cluster center and then solve the first feature value and the second feature value.
The histogram and curve shown in Fig. 26 represent the first feature value and the second feature value calculated by K-means clustering, GOA-K-means clustering and IGOA-K-means clustering, respectively. Comparing the histogram amplitude of the first feature value obtained by K-means clustering in Fig. 26(a), we find that the first feature has changed significantly under the healthy/faulty states of the same measurement point. In our work, these shortcomings are reflected in the fact that the first feature value is hard to discriminate between the healthy state and the faulty states, and the second feature value is also difficult to judge. For example, the first feature values of the measurement point T3 in LALFCW and SALFCW states are smaller than the first feature value in the healthy state, which is inconsistent with the actual situation, and there are serious clustering errors. In addition, the second feature values in the LALFCW state, ICDM state and SCFWT state of the measurement point T3 are all small, which is difficult to distinguish from the second feature value in the healthy state. As such, we conclude that the clustering results of traditional K-means have poor accuracy.
The first feature value and the second feature value of GOA-K-means clustering are shown in Fig. 26(b). It should be noted that GOA-K-means clustering greatly improves the problem of random selection of cluster centers in traditional K-means clustering. However, there are still problems of low sensitivity and low clustering accuracy, which affect the effective extraction of fault characteristics of HVSR vibration signals. For example, in Fig. 26(b), the first feature values of measurement points T1, T2, and 42 in the LALFCW state have a small difference from the first feature value in the HC, which may cause the leakage or loss of latent faults. The second feature value of GOA-K-means clustering has changed significantly and the clustering effect is better than that of K-means clustering. However, the disadvantage is that it cannot accurately reflect the change rule of aggravating failure degree in LALFCW and SALFCW faults. For example, the second feature value in the SALFCW state of the T1 and T2 measurement points is smaller than the second feature value in the LALFCW state.
Based on the above analysis of traditional K-means and GOA-K-means clustering, it is undeniable that there is a certain fault recognition effect, but there are many shortcomings. The IGOA-K-means clustering proposed in this paper may effectively eliminate the shortcomings of randomness, low sensitivity and poor clustering accuracy of the initial cluster center selection, and can accurately identify different types of HVSR faults. In Fig. 26(c), the first feature value of the IGOA-K-means cluster has a positive correlation with the amplitude of the HVSR vibration signal, which effectively reflects the cumulative effect of changes in the state of the iron core and winding. More specifically, the first feature values of the LALFCW state at different measurement points are 1.26-1.85 times the first feature value in the HC, and the first feature value of the SALFCW state is 2.02-2.67 times the first feature value in the HC. It is worth noting that the second feature value gradually increases as the degree of looseness increases, showing a positive correlation, which can effectively distinguish the degree of looseness of the iron core and the winding. For the ICDM fault, its first feature value is between that of the LALFCW states and the SALFCW states, which is 1.41-2.26 times the first feature value in the HC. In addition, in the SCFWT fault, the first feature value of the measurement point 42 shows a large change, which is approximately 4 times that in the HC, which may be used to distinguish other fault types.
In summary, the clustering effect of IGOA-K-means is better than that of GOA-K-means and traditional K-means. It has strong classification and recognition capabilities and can effectively identify different faults. The clustering methods proposed in this paper are summarized in Table 4.

B. COMPARISON WITH CURRENT SIGNAL PROCESSING METHODS
To verify the superiority of the proposed method, we employ the 20kV HVSR vibration data sets in Case II to study the wavelet packet energy (WPE) feature extraction and analysis method widely used in vibration signal detection [24] and the signal processing method based on fast Fourier transform (FFT). At the same time, we use K-nearest neighbor algorithm (KNN) [42], support vector machine (SVM) [24] and decision tree (DT) [43] to classify and test the processed signals respectively. We combine two signal processing methods and three classification learning methods to form a total of six combinations: WPE+KNN, WPE+SVM, WPE+DT, FFT+KNN, FFT+SVM and FFT+DT.
We use 20 measuring points (No.16-35 measuring point) in Case II, totaling 1600 data sets in four different states for calculation and comparison. We employ WPE to extract features from the original vibration signal to obtain a 1600 × 64 feature matrix, and use FFT to extract frequency domain features to obtain a 1600 × 10 feature matrix.  The feature matrices obtained by the two methods were input into KNN, SVM, and DT for training and testing (85% for training and 15% for testing). The test results of the accuracy of the above methods and the proposed method are shown in Table 5.
It can be seen in Table 5 that the accuracy rates of WPE+KNN, WPE+SVM, WPE+DT, FFT+KNN, FFT+SVM and FFT+DT are 73.9%, 78.9%, 88.7%, 78.6%, 83.8% and 91.9%. Compared with the above six methods, the results obtained by the proposed method are satisfactory, and its accuracy rate is at least 8% higher than the above methods. This implies that the proposed method has a better fault diagnosis capability than the above method. It is worth mentioning that the proposed method can accurately classify HVSR core winding faults without using a large amount of data training of machine learning. Moreover, the proposed method has the advantages of unsupervised learning and sensitivity to weak fault characteristics, so it can be tried on-site practical applications.

C. THRESHOLD SETTING AND POPULATION
In the feature space, the threshold of the clustering result of the HVSR vibration signal phase trajectory should be based on the changes of the two features. As mentioned above, the two features can successfully report the changes in the internal state of the HVSR. Depending on the experimental results of Case I, we analysis the two feature values of the T5 measurement point. The first feature value of the LALFCW state is 1.42 times that of the HC, and the corresponding second feature value is 91.2 • . The first feature value of the SALFCW state is 2.01 times that of the HC, and the second feature value is 94.3 • . The first feature value of the ICDM is 1.65 times that of the HC, and the corresponding second feature value is 61.7 • . The first feature value of the SCFWT state is 2.37 times that of the HC, and the corresponding second feature value is 126.9 • .
Through the above analysis, the threshold setting of the first feature value and the second feature value in different states of the T5 measurement point is now performed. Considering that the fluctuation of the HVSR voltage in actual operation may cause slight changes in the vibration signal, we set the threshold to guarantee the accuracy of identifying the fault type. The results are shown in Table 6. Note that the threshold ranges of the first feature value and the second feature value of different measurement points are different, and need to be set according to the specific measurement point. Table 6 introduces the threshold range of the first feature value and the second feature value of the four faulty states of the measurement point T5. It should be noted that if both the first feature value and the second feature value meet the threshold range, they are corresponding fault types; otherwise, they are unknown fault kinds. Of course, the premise is that it is different from the first feature value and the second feature value of HC.
The threshold range setting for unknown fault types will be supplemented in the following two ways: A) Through laboratory fault simulation, the threshold range of unknown faults is supplemented. When the first feature value and the second feature value of the simulated fault are not within the range of the known fault feature value, it is considered an unknown fault. The state of the HVSR is indicated by the first feature value and the second feature value of the same measurement point.
B) The field data is combined with laboratory fault simulation data to fill in fault types together.

D. FIELD APPLICATION CHALLENGES
The challenges that may be faced in the field application of the methods proposed in this paper include the following: 1) Limited by site conditions, how can sensors be arranged on the top of an HVSR?
2) Due to the complex operating environment of an HVSR, how can the electromagnetic interference problem be effectively solved?
For the first challenge, our recommended solutions are as follows: a) When installing the reactor, we fixedly install the vibration acceleration sensor on the top of the reactor in synchronization with the reactor.
b) The reactor is routinely repaired after a power failure every year, and we can install the sensor during the power failure.
In view of the problem of electromagnetic interference in the HVSR field test, we recommend the following measures: a) Select a sensor with strong anti-interference ability, and use shielding technology to install the signal acquisition system in a cabinet made of metal materials to prevent on-site electromagnetic interference. b) Avoid long-distance data transmission lines. If it is unavoidable that long-distance transmission lines must be used, transmission lines with strong anti-electromagnetic interference capabilities (such as shielded twisted pair lines) can be used for transmission. c) Add filter circuits (e.g., RC filter and AC power filter) to the detection device, and use filter technology to suppress electromagnetic interference. d) Determine the noise mode in the reactor vibration signal by collecting a large number of vibration signals during steady-state operation of the reactor, and use noise reduction algorithms (such as median filtering, wavelet decomposition and reconstruction) to preprocess the data before data analysis.

V. CONCLUSION
This paper proposes and describes a vibration analysis method for monitoring the abnormal state of HVSRs core winding. We use PSR and IGOA-K-means clustering to select and identify HVSR core winding fault features. Here, we use the improved GOA to optimize the traditional K-means clustering, which solves the randomness problem of the initial cluster center selection of K-means clustering and improves the accuracy of the clustering results. By employing the two characteristics of the distance between the cluster center displacement vector and the origin and the angle between the faulty state and the healthy state cluster center displacement vector, we effectively identify different fault types and different fault levels of the same fault type. The application of the method proposed in this paper realizes the fault diagnosis of unsupervised learning, which is very important for the online state estimation of an HVSR. Our work may help identify core winding faults in an actual HVSR. We have explored the challenges of applying the methods presented in this paper in the field, which is the focus of ongoing and future work.