GIC Influence on UHV Power Grids Based on Kalman Filter and WAMS Data

Geomagnetically induced currents (GICs) flowing through transformers affect their normal operation and can lead to accidents in power grids. The reactive power loss caused by GICs is called GIC-Q. At present, geomagnetic data are usually used to calculate GICs in research. These methods do not have real-time characteristics and cannot be used for real-time research on GICs. Response mechanisms for GICs in 1000 kV ultra-high voltage (UHV) transformers and for the effect of GIC-Q on large power grids are urgently required to prevent damage. Based on the principle of Kalman filter, this paper presents a model and method of calculating and analyzing GIC disturbance using the reactive power measured by wide area measurement systems (WAMS). The proposed method exhibits superior characteristics for real-time analysis and solves the problem of real-time shortage in GIC research. The model provides a theoretical calculation of GIC-Q disturbance in a 1000 kV autotransformer. Further, this method was verified according to the measured reactive power in three 1000 kV substations in Shandong Province and the measured geoelectric field at the Anqiu geomagnetic observatory station on September 8, 2017. The calculation results show that the reactive power measured by phasor measurement units (PMUs) can be effectively used to calculate and analyze GIC-Q in the transformers, which is significant to prevent disasters in the power grid caused by magnetic storms through dispatch and operation.


I. INTRODUCTION
Geomagnetically induced currents (GICs) are generated in a power system due to magnetic storms. GICs flowing through the transformer cause half-circle saturation resulting in reactive loss, vibration, harmonics, and increases in temperature, which further lead to successive faults in the power system, affecting the safety and stability of the system. In extreme cases, GICs can cause voltage collapse and even large-area blackouts [1]. Research has confirmed that magnetic storm disasters must be seriously considered by power system companies not only at high latitudes, but also in low latitudes and even equatorial regions [2]- [5]. To ensure the safe operation of a power system, it is necessary to conduct real-time monitoring and assess the security risk of the power system during magnetic storms. Accordingly, an intelligent dispatching and security defense scheme for power systems is formulated, and a detailed prevention strategy is developed to prevent risks The associate editor coordinating the review of this manuscript and approving it for publication was Lin Zhang . associated with magnetic storm disasters and avoid further major losses.
The North American Electricity Reliability Committee (NERC) listed magnetic storm events as a low-probability and high-risk type of natural disaster threatening the North American power grid [6]. In 2012, we performed a special reliability assessment for the power grid in North America with regard to disturbances from magnetic storms and noted that the evaluation model is still not accurate [7]. The risks to power systems caused by magnetic storms primarily include important power transformer potential fault sets and voltage stability in power systems [8]. South Africa [9], New Zealand [10], Australia, and other low-and middle-latitude countries have also suffered these problems [11]. Thus, magnetic storm disasters have attracted increasing attention from electric power operators. To study the influence of geomagnetic storms on power grids, the traditional method is to use geomagnetic data from geomagnetic observatory stations to determine the change in the geoelectric field. Then, the Earth surface potential (ESP) can be calculated; the accuracy of 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/ this calculation is key to calculating GICs in power grids [12]- [14], [19]. After years of research and exploration, some scholars have gradually developed several types of ESP analyses and calculation methods, such as plane wave theory, the CIM (complex image method), the SECS (spherical element current system) method, and DSTL (distributed source transmission line) theory. The plane wave theory and CIM are the most widely used. The plane wave method proposed by R. Pirjola et al. has been proven to have sufficient accuracy for application in mid-to low-latitude countries and regions. It is more suitable than CIM for high-latitude regions, such as Finland, Canada, and the United States [15]- [19]. In addition to the theoretical calculation of GICs, monitoring the actual GICs in power grids is important to assess possible effects and potential damage. Measuring GICs is necessary to evaluate interference, such as harmonics, temperature increases, reactive power losses, vibration, and noise. In China, the authors first proposed the installation of monitoring devices to study the impact of geomagnetic storms on China's power grid and developed a device to measure GICs at the neutral point of a transformer [20]. The authors researched and developed a device capable of monitoring harmonics, reactive power, vibration, and noise, including a background analysis system [21]. The data obtained from these devices play an important role in understanding and studying the adverse effects of geomagnetic storms on China's power grid. However, a monitoring system based on GIC effects cannot meet the demand for disaster prevention and rapidly developed defense strategies for the power grid. Thus, new research ideas and methods are necessary. At present, to study GIC-Q (the reactive power loss caused by GICs), the traditional method can be divided into four steps: 1) obtain the geomagnetic data at the geomagnetic observatory stations, 2) calculate the geoelectric field using the geomagnetic data and the Earth conductivity, 3) calculate GICs combined with the grid topology and parameters, and 4) calculate GIC-Q with the operation mode. This process involves complex calculations, specifically in the fourth step. Thus, it is necessary to simplify this problem. With the ongoing development of science and technology, the monitoring level of operation data for modern power systems has significantly improved. The wide area measurement system (WAMS) assists substantially in the real-time synchronization of power systems [22]. In particular, PMUs have been widely used in power systems because of their synchronization, rapidity, and accuracy. These advancements have initiated a technical revolution in the field of dynamic security monitoring. However, due to the large noise in PMU data, GICs cannot be measured directly; thus, this study introduces the Kalman filter to the research. Compared with the traditional filtering method, Kalman filtering has stronger filtering ability and better filtering results because of the probability method [23]. Therefore, based on summarizing traditional methods, this study proposes a method to calculate and analyze GIC disturbances using the PMU based on the Kalman filter principle in a power grid. The feasibility of this approach is verified by analyzing the PMU data at three 1000 kV substations of Shandong Province during the geomagnetic storm on September 8, 2017 with Kp = 8.

II. GIC DISTURBANCES ON REACTIVE POWER IN TRANSFORMERS
The quasi-DC GICs with a frequency of 0.0001-0.01 Hz primarily flow into or out of the power grid through the neutral points of transformers, resulting in undesirable half-circle saturation, which increases the reactive power loss. The purpose of this study is to explore the direct relationship between the characteristics of the geoelectric field and reactive power loss in a transformer measured by PMU. A 1000 kV single-phase autotransformer usually regards the non-excitation voltage regulation at the neutral point, which means that the voltage at the low-voltage side is unchanged during voltage regulation. Power transmission is conducted between the highvoltage and medium-voltage side. Reactive power compensation devices are usually installed at the low-voltage side. The low-voltage winding is connected in series with the compensation transformer and in parallel with the voltage-regulated transformer. Based on the wiring principle of the autotransformer, as shown in Fig. 1 and 2, the PMU measuring device is installed at the high-voltage and medium-voltage sides in the 1000 kV UHV transformer. PMU measures electrical parameters including fundamental reactive power, voltage, and current amplitude at both the high and medium voltage sides, while there is no PMU measurement at the low voltage side. GICs cause DC bias in the transformer. If the fundamental component of excitation current increases, the reactive power loss will increase accordingly. To explore the reactive power caused by magnetic storms during practical operation and dispatch, the GICs, GIC-Q, and PMU measured reactive power are combined to analyze the correlation between the variation of excitation currents and GIC-Q.
Transformer reactive power (Q b ) generally includes excitation reactive power and magnetic flux leakage reactive power. As shown in Fig. 1, U 1a and U 2a are the voltages at the high-voltage and medium-voltage sides, respectively. In Fig. 2, it is shown that the excitation power of a regulated-voltage transformer comes from the parallel winding between the low-voltage winding and regulated-voltage winding. They are then connected in series with the compensation winding to stabilize the voltage at the low-voltage side. As the PMU is generally located at the terminal of the transmission line, the measured values can be considered as the reactive power of the load. Thus, when the voltage at the low-voltage side is stable and reactive power compensation is fixed, we neglect the variation of reactive power in the low-voltage winding and establish the GIC-Q calculation model: where Q b represents the reactive power of the transformer, Q 1 indicates the reactive power flowing into the high-voltage winding, and Q 2L indicates the reactive power of the load at the medium-voltage side. A single-phase autotransformer is used in the 1000 kV substation. The wiring diagram is demonstrated in Fig. 1. We can deduce the current relationship: where I 1a and I 2 a are respectively currents at the high-voltage and medium-voltage side, I 0 represents the excitation current, and k a is the transformer ratio. During normal operation, the core flux presents a sinusoidal wave, as shown in Fig. 3(a) indicated by a solid line. The linear region in Fig. 3(b) denotes the operation point of the transformer. In this case, the exciting current is depicted by the solid sine wave shown in Fig. 3(c). However, when DC bias occurs, the core flux is shifted upward, as indicated in Fig. 3(a) with a dashed line, which causes the work point to enter the saturation area corresponding to the A-B section shown in Fig. 3(b). This leads to the distortion of the excitation current, resulting in the spike wave illustrated in Fig. 3(c).
Nonlinear inductance is used to express the saturation characteristics of the core. Regardless of core loss and eddy current loss, the voltage of the transformer is assumed to have sinusoidal variation; therefore, the relationship can be determined as follows: where U 0 refers to the effective value of voltage, N is the number of turns, 0 is the effective value of flux without load in the transformer, and ω is the angular frequency. The leakage reactance of the air gap in the transformer L a is far less than the reactance of the excitation branch L m . Hence, neglecting the effect of L m , the equivalent inductance in the case of DC bias can be represented by the operating point on the characteristic curve: where kn represents the saturation flux corresponding to point A in Fig. 3(b), and dc indicates the DC flux due to GICs. When the effective voltage of the excitation winding of the transformer is U 0 , the equivalent reactance is expressed as follows: According to the relationship between reactance and terminal voltage, the reactive power loss with DC bias Q b can be described as follows: The DC flux and GICs in the transformer have the following relationship: In the above formula, L m and I GIC denote the reactance of the excitation branch and size of GICs, respectively. Therefore, the reactive power can be represented by GICs: Combined with (3), the reactive power can be expressed by the geoelectric field: The coefficients a and b in the formula are only related to the distribution, topology, and parameters of the power network. E x and E y represent the eastward and northward components, respectively, of the geoelectric field.
The deduction above reveals an approximately linear relationship between the increased reactive power and GICs and a linear relationship between the geoelectric field and GICs. Thus, the reactive power variation is also approximately linear with the variation of the geoelectric field.

III. KALMAN FILTERING METHOD
Through the input and output observation data of a system, the Kalman filtering algorithm can optimally estimate the state of the system using the state equation of the system. The Kalman filtering method combines the use of statistics with a filtering algorithm. Because the observation data includes noise and interference in the system, the optimal estimation can also be regarded as a filtering process. Data filtering is a type of data processing technology used to remove noise and restore real data. Kalman filtering can estimate the state of a dynamic system from a series of data with measurement noise when the measurement variance is known. Because it is convenient to program and update and process the data collected from a site in real time, the Kalman filtering method is currently the most widely used method in many fields, such as communication, navigation, guidance, and control.
First, we introduce a state transition equation and state observation equation.
where θ (k) is the real state variable at the time k; A and C represent the state transition matrix and observation matrix, respectively; B represents the system parameter; u(k) is the control variable of the system, which is set to zero in the absence of a control variable; w (k), z (k), and v (k) denote the process noise, observation value of the state variable, and observation noise at time k. Then, the filtering equations and gain equations are represented as follows.
Filtering equations: whereθ (k) denotes the estimated value of the state variable, and e (k) is the residual. Gain equations: where K (k) represents the Kalman gain, M 1 (k) is the function of M cov (k-1); R(k), N (k), and M cov (k) are the covariance of v(k), w(k), and e(k), respectively; and e(k) is the estimation error.
Through the calculation of estimated values and errors, we can obtain the optimal values which are closest to the real values. The optimal values are the filtering results we need.

IV. MODELING OF OPTIMAL GIC AND GIC-Q
The Kalman filter uses the real-time observation data obtained by the system to continuously modify the state of the system, a process which can dynamically improve the accuracy of the estimated values describing the system state. The current GIC calculation method primarily depends on the observation of geomagnetic data, but there are not enough geomagnetic observatory stations, which are usually away from the substation. This method cannot reflect the influence of geomagnetic storm on power grid in real time, which is not conducive to the safe dispatching and control of the power grid. However, the current data obtained from PMU contain significant disturbances; thus, the accuracy is extremely poor, making it difficult to use the data directly. Thus, we introduce the principle of the Kalman filter to the research.
In this study, the Kalman filter method and PMU data are combined to develop a GIC and GIC-Q data recognition method. The details are described in following steps.
Step 1: By handling the PMU data, we can obtain the GICs and GIC-Q during the magnetic storms, denoted as I PMU and Q PMU ; these are the measured values. From previous studies [3], it is well-known that geoelectric fields can be determined from geomagnetic data and the Earth's conductivity. Then, based on the topology and parameters of the power system, the GICs and reactive power can be calculated, denoted as I E and Q E , respectively; these are used as the initial values.
Step 2: At the time (k-1), the known values Q E (k-1) and I E (k-1) are considered as estimation values. A is the state transition matrix determined from the state of the power system. Then, we can obtain the predictive values at the next moment.
Step 3: The predictive values are calculated by the following equations, in which C represents the relationship between the observation variable and state variable.
The residual is then obtained.
Step 4: The gain is updated as K (k)e (k), and the variables R and N in (17)(18) are obtained from several groups of initial observation information [24], [25]. K (k) changes with time. The optimal filtering values can then be obtained: The optimal results obtained from the Kalman filter are regarded as the GICs and GIC-Q in the power grid in which the disturbance has been removed. Using the GIC-Q as an example, a flowchart is presented in Fig. 4. To verify the accuracy of the results, we can substitute the measured optimal GICs and GIC-Q into (9). If the deviation is less than 5%, it can be proven that the filtering results are sufficiently accurate.
Accordingly, we obtain a real-time method of monitoring the influence of GIC based on the Kalman filter principle, which uses PMU real-time data as the data source with realtime synchronization. Thus, this method is of great significance to estimate the operational risks caused by geomagnetic storms.

V. CASE STUDY
We used the North China UHV power grid as an example and performed a simulation. The grid structure is illustrated in Fig. 5. In this paper, the Hai-he, Chang-le, and Lin-yi stations are selected as the research objects. The PMU data in these three stations are used to analyze the disturbance caused by GICs. The east component (y is eastward) of geoelectric fields obtained from the An-qiu geomagnetic observatory near the stations (using China Standard Time, which is 8 h ahead of UTC) during the magnetic storms on September 8, 2017 is shown in Fig. 6. The data from 5:00 to 11:00 are provided in comparison with the reactive power measured by PMU in Fig. 7. It can be observed that this was a sudden commencement magnetic storm which began at around 7:02. In the initial phase, the geoelectric field doubled to 284 mV/km, and the disturbance lasted for several hours. According to the K P index data provided by   the Space Environment Prediction Center (SEPC), there was a large geomagnetic storm that lasted for 9 h (K p = 8), a small geomagnetic storm of 6 h long (K p = 5), and geomagnetic activity of 3 h long on September 8, 2017.
Noting that the sampling intervals of geoelectric field data and PMU data are respectively 1 min and 20 s, we collect data with an interval of 1 min. We then perform a simple filtering based on the aforementioned models. In Fig. 7-9, the black line represents the calculated GIC-Q with the measured GICs   at the neutral point of the transformer, the blue line represents the measured reactive power, and the green line represents the GIC-Q obtained through the Kalman filter. Notably, the geoelectric field curve in Fig. 6 remains horizontal for several minutes. The analysis made in [14] proposed that this may be caused by temporary faults, such as poor contact due to aging of electrode measurement. According to engineering experience, because the acquisition method using PMU is more direct, the weight value is larger, which is set as 0.6. The Kalman gain and system noise of operation data from the three UHV substations are listed in Table 1.
From the data above, we can observe that the geomagnetic field and reactive power varied almost with the same trend at the beginning of the magnetic storm-that is, 7:02 Beijing time. Moreover, they suddenly changed at almost the same time. Although the waveforms at the front and rear ends were relatively inconsistent, the reactive power and excitation current changed correspondingly at the time when the geoelectric field changed. The consistency of the waveform verifies the model. We initially analyze that the difference  between the collected reactive power and calculated reactive power at 8:15 was caused by reactive power compensation at the low-voltage side. The voltage at the low-voltage side is not stable.
Consider Chang-le station as an example. At the beginning of the magnetic storm, the reactive power suddenly increased by 30 Mvar. The geoelectric field varied almost without fluctuation at 7:10, which is not consistent with the variation of collected reactive power and calculated reactive power. This may be caused by acquisition errors of the devices during the observation, a factor which does not affect the discussion in this section. After 7:15, the geoelectric field decreased significantly corresponding to an approximate trend with the excitation current. However, there was no significant change in the reactive power at this time due to other sources, i.e., reactive power compensation. However, it began to decrease significantly after 7:30, similar to the geoelectric field.
To ensure the accuracy of GICs and GIC-Q obtained by filtering, the data must be verified. We compare the GICs obtained by filtering with the GICs calculated from the geoelectric field. The results for the three substations are demonstrated in the following figures.
As shown in Fig. 10-12, in the three UHV substations, the calculated GICs with the geoelectric field are highly consistent with currents at the neutral point of the transformer obtained by the Kalman filter. However, there is a significant difference in the second half of the waveform. We later discussed this result with the dispatching personnel in the North China power grid. Our analysis determined that this may be caused by the action of reactive power compensation devices. This proves the validity of the method in this paper.

VI. DISCUSSION AND CONCLUSION
For the dispatching and operation management of power grids, it is critically important for security to monitor the real-time impact of geomagnetic storms on the power grid. This study introduces the Kalman filter into the research on GICs. Using the WAMS with the function of data collection, we propose a GIC analysis method by combining the Kalman filter with PMU data. Considering the geomagnetic storm on September 8, 2017 as a case study, the PMU-measured data obtained from three 1000 kV substations in Shandong power grid are analyzed. The data consistency verifies the accuracy of the proposed method. Compared with the traditional calculation method [14], this method substantially simplifies the analysis process and provides a foundation and possibility for real-time analysis and identification of GIC-Q disturbances. The GIC monitoring method proposed in this study can effectively improve the real-time monitoring capacity of GICs in power grids, a result which is significant for real-time measurement to prevent and control the impact of geomagnetic storms.
The GIC-Q calculation method based on the Kalman filter and PMU data is simple and effective, and it has outstanding advantages listed as follows: (a) Compared with other monitoring methods, there is no need to install more relevant equipment. The investment cost is low, and the data are accurate and reliable. This method exhibits global synchronization with the time scale of GPS. (b) Compared with traditional GIC-Q disturbance analysis, the proposed method is simple and fast and is suitable for real-time computation and simultaneous comparative analysis of large-scale grids.
SHUMING ZHANG was born in Jilin, China, in January 1987. He received the master's degree in electrical engineering from Northeast Electric Power University, Jilin, in 2013. He is currently pursing the Ph.D. degree in electrical engineering with North China Electric Power University, Beijing.
His research interests include electromagnetic analysis in power systems, power system monitoring and mitigation, modeling geomagnetically induced currents in the power grid, and assessing the effect on power systems. He has been a Professor and a Ph.D. Supervisor with the School of Electrical and Electronic Engineering, NCEPU. His research interests include safe operation and hazard prevention of power systems and power quality. Prof. Liu is a Senior Member of the Chinese Society for Electrical Engineering. He is also a Commissioner of the Chinese Space Weather Committee and the National Space Weather Monitoring and Pre-warning Technology Standard Committee.
HONGPENG LI was born in Jilin, China. He received the master's degree in electrical engineering from Northeast Electric Power University, Jilin, in 2013. His research interests include electromagnetic analysis in power systems, and projects of power transmission and transformation.