# Patent

## Title
System frequency measurement method, synchrophasor measurement method and device thereof

## Abstract

Provided are a system frequency measurement method, synchrophasor measurement method and device. The synchrophasor measurement method comprises: sampling voltage and current signals based on an absolute time reference; removing out-of-band interference in the sampled voltage and current signals; first-resampling at least one of the voltage and current signals with the out-of-band interference being removed based on an estimated system frequency; calculating phasor for the signals obtained from the first resampling and updating the estimated system frequency; for a reporting timing, second-resampling the voltage and current signals with the out-of-band interference being removed based on the updated estimated system frequency; and calculating synchrophasor of the voltage and current signals obtained from the second resampling. The system frequency measurement method, synchrophasor measurement method and device are not only applicable to the stable operation of power systems but also to the dynamic operation thereof, and provide measurement result with high precision.

## Background
FIELD OF INVENTION
The present invention relates to power system automatic measurement, and more specifically to a high-precision system frequency measurement method, synchrophasor measurement method and device thereof.
BACKGROUND
Synchronized phasor (synchrophasor) measurement unit (PMU) has become an important measurement component in a power system and can provide high precision voltage and current phasors and frequency which are synchronized with an absolute time reference. The applications of PMU have been well recognized for real time dynamic monitoring, operation and control of the power system.
Traditionally, a power system may be monitored through a combination of non-synchronized measurements, for example, generator/load power and voltage magnitudes are fed into a software model, and then the complete voltage angle and magnitude information for each bus for specified load and generator real power and voltage conditions are calculated. However, these quantities can be calculated directly using PMU.
With the growing demand for synchrophasor measurement technology, a standard related to this measurement has been set up by IEEE to ensure the interoperability of devices from different vendors since 2005. There are now many PMUs in compliance with the standard in the market. However, as the standard is only designed for steady-state measurement and cannot meet further demand for dynamic-state measurement (such as during power swing, phase angle change, frequency consecutive change and even fault disturbance) with a lack of comparability and interoperability among synchrophasors measured by different products, some vendors have begun to use their enterprise standard defined by themselves. Based on the new demand, the IEEE and IEC have been working together to establish a new standard including dynamic-state measurement. This standard may be released soon.
A synchrophasor measurement method and device applicable to both the steady-state measurement and dynamic-state measurement in a power system is thus needed. Such a method and device is able to comply with the new standard that covers the dynamic-state measurement in the power system.
SUMMARY
In consideration of the above issues, the present invention is proposed. It aims at providing a high precision synchrophasor measurement method and device, which obtains synchronized voltage phasor and synchronized current phasor by removing the out-of-band interference from raw samples of voltage and current signals, then a first resampling to obtain an estimated system frequency followed by a second resampling for each reporting time. The synchrophasor measurement method selects different calculation window length for steady operation and dynamic operation of the power system respectively, and calculates the estimated system frequency for each sample, thus high-precision may still be achieved while the method is applied to the dynamic operating condition and the requirements of the new standard may be met.
According to one aspect of the present invention, it is provided a synchrophasor measurement method, comprising: sampling voltage and current signals in a power system based on an absolute time reference; removing out-of-band interference in the sampled voltage and current signals; first-resampling at least one of the voltage and current signals with the out-of-band reference being removed based on an estimated system frequency; calculating phasor for signals obtained from the first resampling, and updating the estimated system frequency; for a reporting time, second-resampling the voltage and current signals with the out-of-band reference being removed based on the updated estimated system frequency; and calculating synchrophasor for the voltage and current signals obtained from the second resampling.
Preferably, in the synchrophasor measurement method, the out-of-band interference in the sampled voltage and current signals is removed with a bandpass filter which is customized according to performance class, nominal system frequency of the power system and reporting rate.
Preferably, in the synchrophasor measurement method, calculating phasor for signals obtained from the first resampling comprises: determining the length of a first calculation window according to an operation status of the power system; and calculating phasor for signals obtained from the first resampling in the first calculation window.
Preferably, in the synchrophasor measurement method, updating the estimated system frequency comprises: for each sample obtained from the first resampling, calculating the estimated system frequency; smooth filtering the calculated estimated system frequency to update the estimated system frequency; and setting a corresponding time tag for the updated estimated system frequency, said the time tag compensating the delay introduced by the out-of-band interference removing and the delay introduced by the smooth filtering.
Preferably, in the synchrophasor measurement method, second-resampling the voltage and current signals with the out-of-band reference being removed comprises: estimating the system frequency at a reporting time based on the updated estimated system frequency and its corresponding time tag; and seond-resampling the voltage and current signals with the out-of-band interference being removed based on the estimated system frequency at the reporting time, which comprises: determining the center of a second calculation window based on the reporting time; determining the length of the second calculation window based on a reporting rate, and second-resampling the voltage and current signals with the out-of-band interference being removed based on the estimated system frequency at the reporting time in the second calculation window.
Preferably, the synchrophasor measurement method further comprises calculating rate of change of system frequency by using curve fitting based on the updated estimated system frequency.
According to another aspect of the present invention, it is provided a synchrophasor measurement device, comprising: a raw sampler, for sampling voltage and current signals in a power system based on an absolute time reference; a filter, for removing out-of-band interference in the sampled voltage and current signals; a first resampler, for first-resampling at least one of the voltage and current signals with the out-of-band reference being removed based on an estimated system frequency; a system frequency estimation means, for calculating phasor for signals obtained from the first resampling, and updating the estimated system frequency; a second resampler, for second-resampling the voltage and current signals with the out-of-band reference being removed based on the updated estimated system frequency for a reporting time; and a synchrophasor estimation means, for calculating synchrophasor for the voltage and current signals obtained from the second resampling.
Preferably, said filter is a bandpass filter being customized according to performance class, nominal system frequency of the power system and reporting rate.
Preferably, the system frequency estimation means determines the length of a first calculation window according to an operation status of the power system, and calculates phasor for signals obtained from the first resampling in the first calculation window.
Preferably, the system frequency estimation means calculates the estimated system frequency for each sample obtained from the first resampling, smooth filters the calculated estimated system frequency to update the estimated system frequency, and sets a corresponding time tag for the updated estimated system frequency, said time tag compensating the delay introduced by the out-of-band interference removing and the delay introduced by the smooth filtering.
Preferably, the second resampler estimates the system frequency at a reporting time based on the updated estimated system frequency and its corresponding time tag, determines the center of a second calculation window based on the reporting time, determines the length of the second calculation window based on a reporting rate, and secondly resamples the voltage and current signals with the out-of-band interference being removed based on the estimated system frequency at the reporting time in the second calculation window.
According to still another aspect of the present invention, it is provided a system frequency measurement method of a power system, comprising: sampling voltage and current signals in a power system based on an absolute time reference; removing out-of-band interference in the sampled voltage and current signals; first-resampling at least one of the voltage and current signals with the out-of-band reference being removed based on an estimated system frequency; and calculating phasor for signals obtained from the first resampling, and updating the estimated system frequency.
Preferably, in the system frequency measurement method, calculating phasor for signals obtained from the first resampling comprises: determining the length of a first calculation window according to an operation status of the power system; and calculating phasor for signals obtained from the first resampling in the first calculation window.
Preferably, in the system frequency measurement method, updating the estimated system frequency comprises: for each sample obtained from the first resampling, calculating the estimated system frequency; smooth filtering the calculated estimated system frequency to update the estimated system frequency; and setting a corresponding time tag for the updated estimated system frequency, said time tag compensating the delay introduced by the out-of-band interference removing and the delay introduced by the smooth filtering.
Preferably, the system frequency measurement method further comprises: for a reporting time, estimating the system frequency at the reporting time based on the updated estimated system frequency and its corresponding time tag; second-resampling the voltage and current signals with the out-of-band reference being removed based on the estimated system frequency at the reporting time; and calculating synchrophasor for the voltage and current signals obtained from the second resampling.
Preferably, in the system frequency measurement method, second resampling the voltage and current signals with the out-of-band interference being removed comprises: determining the center of a second calculation window based on the reporting time; determining the length of the second calculation window based on a reporting rate, and second-resampling the voltage and current signals with the out-of-band interference being removed based on the estimated system frequency at the reporting time in the second calculation window, wherein, calculating synchrophasor for voltage and current signals obtained from the second resampling comprises: calculating synchrophasor for the voltage and current signals obtained from the second resampling in the second calculation window.
Preferably, the system frequency measurement method further comprises: calculating rate of change of system frequency by using curve fitting based on the updated estimated system frequency.
With the synchrophasor measurement method and device according to the present invention, the raw samples of voltage and current signals are first filtered to remove the out-of-band interference, then a specific calculation window length is selected for each sample in the first resampling according to the operation status of the power system to obtain the estimated system frequency, and the synchronized voltage phasor and synchronized current phasor at the reporting time are obtained with the second resampling, thus it is guaranteed to be not only applicable to the steady state of the power system but also the dynamic state thereof, and at the same time, high-precision synchrophasor measurement result is provided.



BRIEF DESCRIPTION OF THE DRAWINGS
The above and other purposes, characteristics and advantages of the present invention will become more apparent to those skilled in the art with the detailed description to the embodiments of the present invention by referring to the accompanying drawings, in which:
FIG. 1 is a flow chart of the synchrophasor measurement method according to an embodiment of the present invention;
FIGS. 2A and 2B are examples of a band pass filter designed for a first performance class (M class);
FIGS. 3A and 3B are examples of a band pass filter designed for a second performance class (P class);
FIG. 4 is a detailed flow chart of the step for updating the estimated system frequency according to an embodiment of the present invention;
FIG. 5 is a detailed flow chart of the step for second-resampling and synchrophasor calculating according to an embodiment of the present invention;
FIG. 6 is a schematic block diagram of a synchrophasor measurement device according to an embodiment of the present invention.



DETAILED DESCRIPTION
A synchrophasor measurement method and device according to embodiments of the present invention will now be described by referring to the accompanying figures.
First, a synchrophasor measurement method according to an embodiment of the present invention will be described by referring to FIG. 1.
At Step S105, the synchrophasor measurement method 100 according to the embodiment of the present invention is initialized.
At Step S110, voltage and current signals in a power system are sampled based on an absolute time reference. As known in the art, the voltage and current signals in the power system are sampled in a synchronized sampling system using an A/D converter so as to obtain raw samples of the voltage and current signals. Typically, the triggering sampling pulse used in the synchronized sampling system is locked in phase with an absolute time reference, such as the 1PPS signal of GPS or other absolute time reference from other sources. For example, the triggering sampling pulse used in the synchronized sampling system may be a divided signal of the 1PPS signal of GPS; or a pulse signal generated from local oscillators, and in this case, the local oscillator is synchronized with the absolute time reference at every predetermined interval (for example, 1 s).
Then at Step S120, the out-of-band interference is removed from the raw samples of the voltage and current signals. The out-of-band interference signal is a signal at frequency f where: |f−f0|≧Fs/2, wherein f0 is the nominal system frequency of the power system; and Fs is the reporting rate of PMU (synchronized phasor measurement unit), i.e. the reporting times per second. For example, the reporting rate of PMU may be 10 Hz, 20 Hz, 60 Hz, 120 Hz and so on. Taking an example of 10 Hz as the reporting rate, then the reporting times per second is 10 times with the corresponding reporting time being at 0.1 s, 0.2 s, 0.3 s, . . . 0.9 s, 1 s. More generally, when the reporting interval is labeled as T0, the reporting time may be indicated correspondingly as T0, 2T0, 3T0, . . . , kT0, (k+1)T0, . . . , NT0.
At Step S130, based on an estimated system frequency, a first resampling is performed on at least one of the voltage and current signals with the out-of-band reference being removed. After the initialization of the system, the estimated system frequency is firstly set to the nominal system frequency of the power system.
Then, at Step S140, phasor is calculated for signals obtained from the first resampling, and the estimated system frequency is updated.
Next, the synchrophasor measurement method 100 according to the embodiment of the present invention returns to Step S130, and performs the first resampling with the estimated system frequency updated at Step S140.
On the other hand, after Step S140, the synchrophasor measurement method 100 according to the embodiment of the present invention advances to Step S150, wherein for each reporting time, the voltage and current signals with the out-of-band reference being removed are second-resampled based on the updated estimated system frequency. For example, at Step S150, for each reporting time, the system frequency at the reporting time is first estimated based on the updated estimated system frequency, then the voltage and current signals with the out-of-band reference being removed are second-resampled based on the estimated system frequency at the reporting time.
At last, at Step S160, synchrophasor is calculated for the voltage and current signals obtained from the second resampling.
Next, exemplary operation of Step S120 for removing the out-of-band reference signal in the synchrophasor measurement method 100 according to the embodiment of the present invention will be briefly illustrated by referring to FIGS. 2 and 3. It should be understood that the operation of Step S120 for removing the out-of-band reference signal in the synchrophasor measurement method 100 according to the embodiment of the present invention is not limited to the examples described as follows, and a person skilled in the art may design different filters upon actual requirements.
As described above, the out-of-band interference signal is a signal at frequency f where: |f−f0|≧Fs/2, wherein f0 is the nominal system frequency of the power system; and Fs is the reporting rate of PMU, i.e. the reporting times per second.
The filters may be designed with considerations to specified sampling rate, reporting rate, nominal system frequency of the power system and performance class. Regarding the performance class, there are typically two different kinds in the synchrophasor measurement, namely a first performance class (M class) intended for high precision synchrophasor measurement, and a second performance class (P class) intended for fast response synchrophasor measurement. The M class filter requires a higher measurement precision and is more tolerant of reporting latency, while the P class filter requires a faster response (i.e. a shorter reporting latency) and lower measurement precision.
For example, for a M class filter, when the sampling rate is 4800 Hz, the reporting rate is 60 Hz and the nominal system frequency is 60 Hz, a finite impulse response (FIR) filter with the amplitude-frequency response curve as shown in FIG. 2A may be used.
As an alternative example, for a M class filter, when the sampling rate is 4800 Hz, the reporting rate is 20 Hz and the nominal system frequency is 60 Hz, a FIR filter with the amplitude-frequency response curve as shown in FIG. 2B may be used.
As it is considered that P class is intended for applications requiring fast response and allowing short reporting latency. Therefore, there is no requirement for P class to remove all out-of-band interference.
For example, for a P class filter, when the sampling rate is 4800 Hz, the reporting rate is 60 Hz and the nominal system frequency is 60 Hz, a filter with the amplitude-frequency response curve as shown in FIG. 3A may be used and it meets the requirement of P class application for very small reporting latency. Compared with the amplitude-frequency response curve of FIG. 2A, a shorter reporting latency is obtained in the curve of FIG. 3A at the cost of the flatness of the curve of FIG. 3A being dramatically reduced. In the case where the reporting rate is less than 60 Hz, a filter with an amplitude-frequency response similar to that of FIG. 3A may also be used.
As an alternative example, for a P class filter, when the sampling rate is 4800 Hz, the reporting rate is 120 Hz and the nominal system frequency is 60 Hz, a filter with the amplitude-frequency response curve as shown in FIG. 3B may be used and it meets the requirement of P class application for very small reporting latency. FIG. 3B shows an extreme case of the bandpass filter i.e. a low-pass filter. Furthermore, in the case where the reporting rate is larger than 60 Hz, a filter with an amplitude-frequency response similar to that of FIG. 3B may also be used.
After the filter having been designed, the filtered signals may be obtained as follows:




  
    
      
        
          
            
              
                x
                ⁡
                
                  (
                  k
                  )
                
              
              =
              
              ⁢
              
                
                  ∑
                  
                    m
                    =
                    0
                  
                  
                    L
                    -
                    1
                  
                
                ⁢
                
                  
                    y
                    ⁡
                    
                      (
                      
                        k
                        -
                        m
                      
                      )
                    
                  
                  ·
                  
                    h
                    ⁡
                    
                      (
                      m
                      )
                    
                  
                
              
            
          
        
        
          
            
              =
              
              ⁢
              
                
                  ∑
                  
                    m
                    =
                    0
                  
                  
                    
                      
                        L
                        -
                        1
                      
                      2
                    
                    -
                    1
                  
                
                ⁢
                
                  
                    [
                    
                      
                        y
                        ⁡
                        
                          (
                          
                            k
                            -
                            m
                          
                          )
                        
                      
                      +
                      
                        y
                        ⁡
                        
                          (
                          
                            k
                            -
                            L
                            +
                            1
                            -
                            m
                          
                          )
                        
                      
                    
                    ]
                  
                  ·
                
              
            
          
        
        
          
            
              
              ⁢
              
                
                  h
                  ⁢
                  
                    (
                    m
                    )
                  
                
                +
                
                  
                    y
                    (
                    
                      k
                      -
                      
                        
                          L
                          -
                          1
                        
                        2
                      
                    
                    )
                  
                  ·
                  
                    h
                    (
                    
                      
                        L
                        -
                        1
                      
                      2
                    
                    )
                  
                
              
            
          
        
      
    
    
      
        (
        1
        )
      
    
  




Wherein x(•) is the filtered voltage and current signals; h(•) is the coefficients of the filter; y(•) is raw samples of the voltage and current signals; and L is the number of samples used in the filter.
Next, the exemplary operation of Step S140 of the synchrophasor measurement method 100 according to an embodiment of the present invention will be described. It should be understood that the operation of Step S140 of the synchrophasor measurement method 100 according to the embodiment of the present invention is not limited to the examples described as follows, and a person skilled in the art may design different calculation window.
At Step S140, different window length is selected based on the operating state of the power system, i.e. whether the power system is in a steady state or a dynamic state. For example, if the power system is in a steady state, then one cycle window length may be selected for either M class or P class; while if the power system is in a dynamic state, then 1⁄2 or 1⁄4 cycle or other proper value may be selected for use as the window length.
For the steady state, the one-phase voltage phasor is calculated by using the following formula with a current voltage obtained from the first resampling and N previous voltages obtained from the first resampling:




  
    
      
        
          
            X
            .
          
          ⁡
          
            (
            i
            )
          
        
        =
        
          
            
              2
            
            N
          
          ⁢
          
            
              ∑
              
                k
                =
                0
              
              
                N
                -
                1
              
            
            ⁢
            
              
                x
                ⁡
                
                  (
                  
                    i
                    -
                    N
                    +
                    k
                  
                  )
                
              
              ·
              
                e
                
                  
                    -
                    j
                  
                  ⁢
                  
                    
                      2
                      ⁢
                      π
                    
                    N
                  
                  ⁢
                  
                    (
                    
                      i
                      +
                      k
                      -
                      
                        1
                        2
                      
                    
                    )
                  
                
              
            
          
        
      
    
    
      
        (
        2
        )
      
    
  




Wherein {dot over (X)}(i) is the current voltage phasor of one-phase, N is the number of samples in one cycle and i is the index number of the current voltage obtained from the first resampling.
Similarly, the current obtained from the first resampling and N previous currents obtained from the first resampling may be used to calculate the phasor of the one-phase-current, by using formula (2).
For the dynamic state, when the calculation window length is selected to be 1⁄2 cycle, the one-phase voltage phasor may be calculated by using the following formula with the current voltage obtained from the first resampling and (N/2) previous voltage samples obtained from the first resampling:




  
    
      
        
          
            X
            .
          
          ⁡
          
            (
            i
            )
          
        
        =
        
          
            
              2
              ⁢
              
                2
              
            
            N
          
          ⁢
          
            
              ∑
              
                k
                =
                0
              
              
                
                  N
                  /
                  2
                
                -
                1
              
            
            ⁢
            
              
                x
                ⁡
                
                  (
                  
                    i
                    -
                    
                      N
                      /
                      2
                    
                    +
                    k
                  
                  )
                
              
              ·
              
                e
                
                  
                    -
                    j
                  
                  ⁢
                  
                    
                  
                  ⁢
                  
                    
                      2
                      ⁢
                      π
                    
                    N
                  
                  ⁢
                  
                    (
                    
                      i
                      +
                      k
                      -
                      
                        1
                        2
                      
                    
                    )
                  
                
              
            
          
        
      
    
    
      
        (
        3
        )
      
    
  




Similarly, the current obtained from the first resampling and (N/2) previous current samples obtained from the first resampling may be used to calculate the phasor of the one-phase-current, by using formula (3).
For the dynamic state, when the calculation window length is selected to be 1⁄4 cycle, the one-phase voltage phasor is calculated by using the following formula with the current voltage obtained from the first resampling and (N/4) previous voltage samples obtained from the first resampling and N being equal to 80, i.e. 80 samples per cycle:




  
    
      
        
          [
          
            
              
                
                  
                    X
                    R
                  
                  ⁡
                  
                    (
                    i
                    )
                  
                
              
            
            
              
                
                  
                    X
                    I
                  
                  ⁡
                  
                    (
                    i
                    )
                  
                
              
            
          
          ]
        
        =
        
          
            [
            
              
                
                  0.1684
                
                
                  0.1073
                
              
              
                
                  0.1073
                
                
                  0.1684
                
              
            
            ]
          
          ·
          
            [
            
              
                
                  
                    
                      ∑
                      
                        k
                        =
                        
                          i
                          -
                          
                            N
                            /
                            4
                          
                          +
                          1
                        
                      
                      i
                    
                    ⁢
                    
                      
                        cos
                        (
                        
                          
                            
                              2
                              ⁢
                              π
                            
                            N
                          
                          ·
                          
                            (
                            
                              k
                              -
                              
                                1
                                2
                              
                            
                            )
                          
                        
                        )
                      
                      ·
                      
                        x
                        ⁡
                        
                          (
                          k
                          )
                        
                      
                    
                  
                
              
              
                
                  
                    
                      ∑
                      
                        k
                        =
                        
                          i
                          -
                          
                            N
                            /
                            4
                          
                          +
                          1
                        
                      
                      i
                    
                    ⁢
                    
                      
                        cos
                        (
                        
                          
                            
                              2
                              ⁢
                              π
                            
                            N
                          
                          ·
                          
                            (
                            
                              k
                              -
                              
                                1
                                2
                              
                            
                            )
                          
                        
                        )
                      
                      ·
                      
                        x
                        ⁡
                        
                          (
                          k
                          )
                        
                      
                    
                  
                
              
            
            ]
          
        
      
    
    
      
        (
        4
        )
      
    
  




Wherein XR(i) is the real part of the one-phase current voltage and XI (i) is the imaginary part of the one-phase current voltage.
Similarly, the current obtained from the first resampling and (N/4) previous current samples obtained from the first resampling may be used to calculate the phasor of the one-phase-current, by using formula (4).
It should be understood that the above formulas (2), (3) and (4) are examples only and the present invention are not limited thereto, and alternatively formulas (2), (3) and (4) may be achieved with recursion.
Next, it will be described as how to update the estimated system frequency with the calculated phasor resulted from the above formulas by referring to FIG. 4.
At Step S140, for each sample obtained from the first resampling, the estimated system frequency is calculated. For example, the estimated system frequency is updated according to the angle difference between the current estimated phasor and the previous phasor being spaced a predetermined time interval from the current instant.
For example, after an one-phase voltage phasor (i.e. un-synchronized voltage phasor) is obtained, the positive sequence voltage phasor is given by:




  
    
      
        
          
            
              V
              .
            
            1
          
          ⁡
          
            (
            i
            )
          
        
        =
        
          
            1
            3
          
          ⁢
          
            (
            
              
                
                  
                    V
                    .
                  
                  A
                
                ⁡
                
                  (
                  i
                  )
                
              
              +
              
                α
                ⁢
                
                  
                
                ⁢
                
                  
                    
                      V
                      .
                    
                    B
                  
                  ⁡
                  
                    (
                    i
                    )
                  
                
              
              +
              
                
                  α
                  2
                
                ⁢
                
                  
                    V
                    C
                  
                  ⁡
                  
                    (
                    i
                    )
                  
                
              
            
            )
          
        
      
    
    
      
        (
        5
        )
      
    
  




Wherein {dot over (V)}A(i)=VAR(i)+jVAI(i)

    
    
        {dot over (V)}B(i)=VBR(i)+jVBI(i)
        {dot over (V)}C(i)=VCR(i)+jVCI(i)
        α=ej2π/3 
    
    


Wherein {dot over (V)}A(i) is the A-phase voltage phasor, VAR(i) is the real part of the A-phase voltage phasor, VAI(i) is the imaginary part of the A-phase voltage phasor; {dot over (V)}B(i) is the B-phase voltage phasor, VBR (i) is the real part of the B-phase voltage phasor, VBI(i) is the imaginary part of the B-phase voltage phasor; {dot over (V)}C(i) is the C-phase voltage phasor, VCR(i) is the real part of the C-phase voltage phasor, VCI(i) is the imaginary part of the C-phase voltage phasor.
Similarly, after an one-phase current phasor is obtained, the positive sequence current phasor may be calculated by using formula (5).
After the current positive sequence voltage phasor for each sample obtained from the first resampling is calculated, the estimated system frequency is calculated by using formula (6) according to the angle difference between the positive sequence voltage phasor at the current sample time and the positive sequence voltage phasor at a sample time that is spaced the calculation window length from the current sample time.




  
    
      
        
          f
          k
        
        =
        
          
            f
            
              k
              -
              1
            
          
          +
          
            
              f
              
                k
                -
                1
              
            
            ⁢
            
              
                
                  Δ
                  ⁢
                  
                    
                  
                  ⁢
                  
                    φ
                    k
                  
                
                
                  2
                  ⁢
                  π
                
              
              ·
              
                N
                m
              
            
          
        
      
    
    
      
        (
        6
        )
      
    
  




Wherein fk is the estimated system frequency at the current time (k), fk-1 is the estimated system frequency at the previous time (k−1) before the current time (k), Δφk=φk−φk-m is the angle difference of the positive sequence voltage phasor of the current time (k) and the positive sequence voltage phasor of the time (k−m) that is spaced calculation window of m from the current time, m is the calculation window length represented by the number of samples, i.e. the number of samples in the calculation window, for example taking N (one cycle)=80, for the case where the calculation window is selected to be one cycle, m=80; for the case where the calculation window is selected to be 1⁄2 cycle, m=40; and for the case where the calculation window is selected to be 1⁄4 cycle, m=20. However, the present invention is not limited thereto and the length of the calculation window may be selected according to the requirements and the corresponding m may be obtained.
Δφk may be calculated by using the following formula:




  
    
      
        
          
        
        ⁢
        
          
            Δ
            ⁢
            
              
            
            ⁢
            
              φ
              k
            
          
          =
          
            {
            
              
                
                  
                    
                      
                        a
                        ⁢
                        
                          
                        
                        ⁢
                        
                          tan
                          ⁡
                          
                            (
                            
                              
                                sin
                                ⁡
                                
                                  (
                                  
                                    Δ
                                    ⁢
                                    
                                      
                                    
                                    ⁢
                                    
                                      φ
                                      k
                                    
                                  
                                  )
                                
                              
                              
                                cos
                                ⁡
                                
                                  (
                                  
                                    Δ
                                    ⁢
                                    
                                      
                                    
                                    ⁢
                                    
                                      φ
                                      k
                                    
                                  
                                  )
                                
                              
                            
                            )
                          
                        
                      
                    
                    
                      
                        
                          if
                          ⁢
                          
                            
                          
                          ⁢
                          
                            sin
                            ⁡
                            
                              (
                              
                                Δ
                                ⁢
                                
                                  
                                
                                ⁢
                                
                                  φ
                                  k
                                
                              
                              )
                            
                          
                        
                        >
                        
                          10
                          ⁢
                          
                            cos
                            ⁡
                            
                              (
                              
                                Δ
                                ⁢
                                
                                  
                                
                                ⁢
                                
                                  φ
                                  k
                                
                              
                              )
                            
                          
                        
                      
                    
                  
                  
                    
                      
                        
                          π
                          2
                        
                        ⁢
                        
                          sign
                          ⁡
                          
                            (
                            
                              sin
                              ⁡
                              
                                (
                                
                                  Δ
                                  ⁢
                                  
                                    
                                  
                                  ⁢
                                  
                                    φ
                                    k
                                  
                                
                                )
                              
                            
                            )
                          
                        
                      
                    
                    
                      else
                    
                  
                
                ⁢
                
                  

                
                ⁢
                
                  sin
                  ⁡
                  
                    (
                    
                      Δ
                      ⁢
                      
                        
                      
                      ⁢
                      
                        φ
                        k
                      
                    
                    )
                  
                
              
              =
              
                
                  
                    
                      Im
                      (
                      
                        
                          
                            V
                            .
                          
                          1
                        
                        ⁡
                        
                          (
                          k
                          )
                        
                      
                      )
                    
                    ·
                    
                      Re
                      (
                      
                        
                          
                            V
                            .
                          
                          1
                        
                        ⁡
                        
                          (
                          
                            k
                            -
                            m
                          
                          )
                        
                      
                      )
                    
                  
                  -
                  
                    
                      
                        Re
                        (
                        
                          
                            
                              V
                              .
                            
                            1
                          
                          ⁡
                          
                            (
                            k
                            )
                          
                        
                        )
                      
                      ·
                      
                        Im
                        (
                        
                          
                            
                              V
                              .
                            
                            1
                          
                          ⁡
                          
                            (
                            
                              k
                              -
                              m
                            
                            )
                          
                        
                        )
                      
                    
                    ⁢
                    
                      

                    
                    ⁢
                    
                      cos
                      ⁡
                      
                        (
                        
                          Δ
                          ⁢
                          
                            
                          
                          ⁢
                          
                            φ
                            k
                          
                        
                        )
                      
                    
                  
                
                =
                
                  
                    
                      Im
                      (
                      
                        
                          
                            V
                            .
                          
                          1
                        
                        ⁡
                        
                          (
                          k
                          )
                        
                      
                      )
                    
                    ·
                    
                      Im
                      (
                      
                        
                          
                            V
                            .
                          
                          1
                        
                        ⁡
                        
                          (
                          
                            k
                            -
                            m
                          
                          )
                        
                      
                      )
                    
                  
                  +
                  
                    
                      Re
                      (
                      
                        
                          
                            V
                            .
                          
                          1
                        
                        ⁡
                        
                          (
                          k
                          )
                        
                      
                      )
                    
                    ·
                    
                      Re
                      (
                      
                        
                          
                            V
                            .
                          
                          1
                        
                        ⁡
                        
                          (
                          
                            k
                            -
                            m
                          
                          )
                        
                      
                      )
                    
                  
                
              
            
          
        
      
    
    
      
        (
        7
        )
      
    
  




Wherein Im(•) represents the imaginary part of “•”, and Re(•) represents the real part of “•”.
Similarly, the estimated system frequency may be calculated by using formulas (6) and (7) according to the angle difference of the positive sequence voltage phasor of the current time and that of the time spaced the calculation window length from the current time.
However, the present invention is not limited thereto, a person skilled in the art may also use other method to calculate estimated system frequency as the updated estimated system frequency where the positive sequence voltage phasor or the positive sequence current phasor is not calculated.
More advantageously, at Step S420, the estimated system frequency is updated by performing smooth filtering on the calculated estimated system frequency fk. For example, the smooth filtering may be performed with the estimated system frequency fk at the current time (k) and previous m estimated system frequencies by using the following formula:




  
    
      
        
          f
          k
          ′
        
        =
        
          
            
              ∑
              
                i
                =
                
                  k
                  -
                  m
                
              
              k
            
            ⁢
            
              f
              i
            
          
          m
        
      
    
    
      
        (
        8
        )
      
    
  




Wherein f′k is the estimated system frequency after the smooth filtering, and is used as the updated estimated system frequency.
However, the smooth filtering is not limited to the case expressed in formula (8), a person skilled in the art may perform smooth filtering upon actual requirements with the weighted average.
Furthermore, the smooth filtering may not be necessary, in other words, the estimated system frequency calculated with formula (6) may be used directly as the updated estimated system frequency.
Then at Step S430, after the updated estimated system frequency (formula (6)) is obtained, a time tag corresponding to the system frequency is calculated, wherein the time tag compensates the delay introduced by the out-of-band interference removing. Furthermore, in the case where the smooth filtering is performed (the case of formula (8)), the time tag also compensates the delay introduced by the smooth filtering.

t=tres−tfd−tad  (9)

Wherein t is the time tag of the updated estimated system frequency, tres is the time in the raw samples of voltage/current signal corresponding to the sample obtained at the current time k from the first resampling, tfd is the delay introduced by the out-of-band interference removing, and tad is the delay introduced by the smooth filtering.
In addition, after the updated estimated system frequency is obtained (formula (6) or (8)), the new sample interval (i.e. the interval between the current time and the next sampling time) can be calculated by the following formula:




  
    Δ
    ⁢
    
      
    
    ⁢
    t
  
  =
  
    1
    
      N
      ·
      
        f
        k
        ′
      
    
  




Next, an exemplary operation of the Step S150 of the synchrophasor measurement method 100 according to the embodiment of the present invention will be described by referring to FIG. 5, wherein for a specific reporting time, the system frequency at the specific reporting time is estimated based on the updated estimated system frequency at Step S140, then the voltage and current signals with the out-of-band interference being removed is second-resampled. It should be understood that the operation of the phasor calculation Step S150 in the synchrophasor measurement method 100 according to the embodiment of the present invention is not limited to the examples described as follows, a person skilled in the art may design different calculation windows upon actual requirements.
As each reporting time does not necessarily correspond to the time tag of the updated estimated system frequency after the first resampling, the system frequency of each reporting time should be estimated based on the updated estimated system frequency and its corresponding time tag.
First at Step S510, the system frequency at the reporting time is estimated based on the updated estimated system frequency and its corresponding time tag.
For example, the system frequency at the reporting time can be obtained by the following two steps: obtaining two estimated system frequencies with the time tag being closest to the reporting time, and interpolating between the two estimated system frequencies based on the time tags of the two estimated system frequencies. The system frequency of at the reporting time can be estimated by using the following formula:




  
    
      
        
          f
          out
        
        =
        
          
            f
            m
            ′
          
          +
          
            
              
                
                  f
                  
                    m
                    +
                    1
                  
                  ′
                
                -
                
                  f
                  m
                  ′
                
              
              
                
                  t
                  
                    m
                    +
                    1
                  
                
                -
                
                  t
                  m
                
              
            
            ·
            
              (
              
                
                  t
                  c
                
                -
                
                  t
                  m
                
              
              )
            
          
        
      
    
    
      
        (
        11
        )
      
    
  




Wherein tc is the reporting time, fout is the system frequency at reporting time, tm is the closest time tag before the reporting time, tm+1 is the closest time tag after the reporting time, f′m, is the estimated system frequency with the time tag tm, and f′m+1 is the estimated system frequency with the time tag tm+1.
Next, in Steps S520-S540, the voltage and current signals with the out-of-band interference being removed are second-resampled.
Specifically, at Step S520, the center of a second calculation window is determined based on the reporting time. For example, the reporting time may be the center of the second calculation window.
At Step S530, the length of the second calculation window may be determined based on the reporting rate. For example, when the reporting rate is 60 Hz, the length of the second calculation window may be selected as one cycle, i.e. having N samples. However, the present invention is not limited thereto, and with a higher reporting rate such as 120 Hz, the length of the second calculation window may be selected as 1⁄2 cycle, i.e. having N/2 samples.
Then at Step S540, in the second calculation window, the voltage and current signals with the out-of-band interference being removed are second-resampled based on the estimated system frequency at the reporting time.
For example, the sample interval of the second resampling can be calculated with the estimated system frequency at the reporting time by using the following formula:




  
    
      
        
          Δ
          ⁢
          
            
          
          ⁢
          
            t
            out
          
        
        =
        
          1
          
            N
            ·
            
              f
              out
            
          
        
      
    
    
      
        (
        12
        )
      
    
  




Wherein N is the number of samples in one cycle.
Therefore, in the case where the reporting time is the center of the second calculation window, the corresponding time in the raw samples of voltage and current signals for each sample in the second calculation window can be calculated as follows:




  
    t
    pi
  
  =
  
    
      t
      c
    
    +
    
      t
      fd
    
    -
    
      t
      wid
    
    +
    
      
        1
        2
      
      ⁢
      
        i
        ·
        Δ
      
      ⁢
      
        
      
      ⁢
      t
    
  




Wherein tpi is the corresponding time in the raw samples of voltage and current signals for the ith (1≦i≦N) sample in the second calculation window; tc is the center of the second calculation window i.e. the reporting time of the current calculation; tfd is the delay introduced by the out-of-band removing; and twid is the half-window delay, i.e. half of the total time length of the second calculation window.
At last, after all the steps in FIG. 5 are performed, the synchrophasor measurement method according to the embodiment of the present invention advances to Step S160 as shown in FIG. 1, wherein in the second calculation window, a synchrophasor is calculated for the voltage and current signals obtained from the second resampling.
In the second calculation window, a synchrophasor is calculated for the voltage and current signals obtained from the second resampling, so as to be used as the reporting synchrophasor at the reporting time corresponding to the second calculation window. When the calculation window is one cycle, the calculation formula is given by:




  
    
      
        
          
            
              X
              .
            
            out
          
          ⁡
          
            (
            i
            )
          
        
        =
        
          
            
              2
            
            N
          
          ⁢
          
            
              ∑
              
                k
                =
                
                  
                    -
                    N
                  
                  /
                  2
                
              
              
                
                  N
                  /
                  2
                
                -
                1
              
            
            ⁢
            
              
                x
                (
                
                  i
                  +
                  k
                  +
                  
                    1
                    2
                  
                
                )
              
              ·
              
                e
                
                  -
                  
                    j
                    ⁡
                    
                      (
                      
                        
                          
                            
                              2
                              ⁢
                              π
                            
                            N
                          
                          ⁢
                          
                            (
                            
                              i
                              +
                              k
                              +
                              
                                N
                                2
                              
                            
                            )
                          
                        
                        -
                        π
                        +
                        
                          π
                          N
                        
                      
                      )
                    
                  
                
              
            
          
        
      
    
    
      
        (
        13
        )
      
    
  




Wherein {dot over (X)}out (i) is the one-phase current voltage or current, N is the number of samples in one cycle, k is the index number of the current voltage or current obtained from the second resampling and i corresponds to the index number of reporting time.
Furthermore, a high-precision rate of change of frequency (ROCOF) may be obtained after Step S140 wherein the estimated system frequency is updated by the synchrophasor measurement method according to the embodiments of the present invention, wherein the updated estimated system frequency may be the calculation result from formula (8). In order to obtain high-precision, a curve fitting within a predetermined period may be used and the predetermined period may be selected as two continuous cycles. However, the present invention is not limited thereto, and the predetermined period may be selected as other length upon actual system design requirements.

f=kf·t+fini  (14)

Wherein kf is the slope of the curve fitting i.e. ROCOF; t is the specific time within the predetermined period, f is the estimated system frequency corresponding to the specific time within the predetermined period; and fini is a constant.
Furthermore, because the foregoing amplitude-frequency response cannot be completely flat, especially for filters of class P (as can be seen from FIGS. 3A and 3B), it is necessary to compensate the calculated phasor's amplitude in order to achieve high precision. Accordingly, the amplitude deviation due to filtering can be compensated after Step S160 where the synchrophasor is calculated for the voltage and current signals obtained from the second resampling by the synchrophasor measurement method according to the embodiments of the present invention. For example, such amplitude deviation may result from the unflatness of the passband and most of the edge of a bandpass filter.
Ideally, all designed filter's amplitude response can be pre-calculated off-line and saved in a buffer according to piecewise-defined frequency. Therefore, the amplitude deviation due to unflatness of the passband and most of the edge of a bandpass filter can be calculated using the following formula:

{dot over (X)}′out(fout)={dot over (X)}out(fout)·h(fout)  (15)

Where {dot over (X)}out(fout) is the sychrophasor at frequency fout that is calculated at the reporting time; {dot over (X)}′out(fout) is the synchrophasor at the frequency fout after compensation; and h(fout) is the reciprocal of amplitude coefficient of the designed filter at the frequency fout.
Next, a synchrophasor measurement device 600 according to an embodiment of the present invention will be briefly described by referring to FIG. 6.
The synchrophasor measurement device 600 comprises a raw sampler 610, a filter 620, a first resampler 630, a system frequency estimation means 640, a second resampler 650 and a synchrophasor estimation means 660.
The raw sampler 610 samples the voltage and current signals of the power system based on an absolute time reference.
The filter 620 filters the out-of-band interference of the sampled voltage and current signals. As described above, the filter 620 may be customized according to performance class, nominal system frequency of the power system and reporting rate, and the filter 620 may be a bandpass filter.
The first resampler 630 first-resamples at least one of the voltage and current signals with the out-of-band interference being removed.
The system frequency estimation means 640 calculates phasor for signals obtained from the first resampling and updates the estimated system frequency.
Preferably, the system frequency estimation means 640 determines the length of the calculation window according to the operation status of the power system, calculates phasor for signals obtained from the first resampling in the calculation window, and updates the estimated system frequency according to the calculated phasor.
Preferably, the system frequency estimation means 640 calculates the estimated system frequency for each sample obtained from the first resampling in order to update the estimated system frequency, and sets corresponding time tag for the updated estimated system frequency, the time tag compensating the delay introduced by filtering the out-of-band interference.
More preferably, after the system frequency estimation means 640 calculates the estimated system frequency for each sample obtained from the first resampling, it also performs smooth filtering on the calculated estimated system frequency so as to update the estimated system frequency. In this case, the corresponding time tag set for the updated estimated system frequency also compensates the delay introduced by the smooth filtering.
The second resampler 650 second-resamples the voltage and current signals with the out-of-band interference being removed based on the updated estimated system frequency for a reporting time.
Preferably, the second resampler 650 estimates the system frequency at the reporting time based on updated estimated system frequency and its corresponding time tag, determines the center of the second calculation window based on the reporting time, determines the length of the second calculation window based on the reporting rate, and in the second calculation window second-resamples the voltage and current signals with the out-of-band interference being removed based on the estimated system frequency at the reporting time.
The synchrophasor estimation means 660 calculates the synchrophasor for the voltage and current signals obtained from the second resampling.
The raw sampler 610, the first resampler 630 and the second resampler 650 can be implemented individually in hardware samplers, or the raw sampler 610 may be implemented in a hardware A/D sampler, while the first resampler 630 and the second resampler 650 may be implemented in special purposed processor, general purposed processor or FGPA etc.
The filter 620 may be implemented in special purposed filter circuits, or in special purposed processor, general purposed processor or FGPA etc.
The system frequency estimation means 640 and the synchrophasor estimation means 660 may be implemented individually in separate special purposed processors, or may be implemented on the same general purposed processor or FGPA etc.
With the synchrophasor measurement method and device according to the present invention, the raw samples of voltage and current signals are first filtered to remove the out-of-band interference, then specific calculation window length is selected for each sample in the first resampling according to the operation status of the power system to obtain the estimated system frequency, and the synchronized voltage phasor and synchronized current phasor at the reporting time are obtained with the second resampling, thus it is guaranteed to be not only applicable to the steady state of the power system but also the dynamic state thereof, and at the same time, high-precision synchrophasor measurement result is provided.
While exemplary embodiments have been described by referring to the accompanying drawings herein, it should be understood that the above exemplary embodiments are illustrative only and are not intended to limit the scope of the present invention thereto. An ordinary person skilled in the art may make various modifications and alterations without departing from the scope and spirit of the present invention. It is intended that all these modifications and alterations be encompassed within the appended claims.

## Claims

1. A synchrophasor measurement method, for a power system for real time dynamic monitoring, operation and control of the power system, comprising:
sampling voltage and current signals in a power system based on an absolute time reference;
removing out-of-band interference in the sampled voltage and current signals;
first-resampling at least one of the voltage and current signals with the out-of-band interference being removed based on an estimated system frequency;
calculating phasor for signals obtained from the first resampling, and updating the estimated system frequency;
for a reporting time, second-resampling the voltage and current signals with the out-of-band interference being removed based on the updated estimated system frequency;
calculating synchrophasor for the voltage and current signals obtained from the second resampling,
wherein updating the estimated system frequency comprises:
for each sample obtained from the first resampling, calculating the estimated system frequency;
smooth filtering the calculated estimated system frequency to update the estimated system frequency; and
setting a corresponding time tag for the updated estimated system frequency, the time tag compensating the delay introduced by the out-of-band interference removing and the delay introduced by the smooth filtering,

and using the synchrophasor to execute a steady operation and a dynamic operation of the power system.

2. The synchrophasor measurement method according to claim 1, wherein,
the out-of-band interference in the sampled voltage and current signals is removed with a bandpass filter which is customized according to performance class, nominal system frequency of the power system and reporting rate.

3. The synchrophasor measurement method according to claim 1, wherein calculating phasor for signals obtained from the first resampling comprises:
determining the length of a first calculation window according to an operation status of the power system; and
calculating phasor for signals obtained from the first resampling in the first calculation window.

4. The synchrophasor measurement method according to claim 1, wherein calculating the estimated system frequency for each sample obtained from the first resampling comprises:
calculating the estimated system frequency according to the angle difference of the phasor between the sample obtained from the first resampling and a previous sample obtained from the first resampling being spaced a predetermined time length from the sample.

5. The synchrophasor measurement method according to claim 1, wherein second-resampling the voltage and current signals with the out-of-band interference being removed comprises:
estimating the system frequency at a reporting time based on the updated estimated system frequency and its corresponding time tag; and
second-resampling the voltage and current signals with the out-of-band interference being removed based on the estimated system frequency at the reporting time, comprising:
determining the center of a second calculation window based on the reporting time;
determining the length of the second calculation window based on a reporting rate, and
second-resampling the voltage and current signals with the out-of-band interference being removed based on the estimated system frequency at the reporting time in the second calculation window,
wherein, calculating synchrophasor for the voltage and current signals obtained from the second resampling comprises:
calculating synchrophasor for the voltage and current signals obtained from the second resampling in the second calculation window.

6. The synchrophasor measurement method according to claim 5, wherein estimating the system frequency at the reporting time comprises:
obtaining system frequencies estimated at two times closest to the reporting time; and
interpolating between the system frequencies estimated at the two times closest to the reporting time, in order to calculate the system frequency at the reporting time.

7. The synchrophasor measurement method according to claim 1, further comprising:
calculating rate of change of system frequency by using curve fitting based on the updated estimated system frequency.

8. A synchrophasor measurement device, for a power system for real time dynamic monitoring, operation and control of the power system, comprising:
a raw sampling device, for sampling voltage and current signals in a power system based on an absolute time reference;
a filter circuit, for removing out-of-band interference in the sampled voltage and current signals;
a first resampling device, for first-resampling at least one of the voltage and current signals with the out-of-band interference being removed based on an estimated system frequency;
a system frequency estimation means, for calculating phasor for signals obtained from the first resampling, and updating the estimated system frequency;
a second resampling device, for second-resampling the voltage and current signals with the out-of-band interference being removed based on the updated estimated system frequency for a reporting time; and
a synchrophasor estimation means, for calculating synchrophasor for the voltage and current signals obtained from the second resampling,
wherein the system frequency estimation means calculates the estimated system frequency for each sample obtained from the first resampling, smooth filters the calculated estimated system frequency to update the estimated system frequency, and sets a corresponding time tag for the updated estimated system frequency, the time tag compensating the delay introduced by the out-of-band interference removing and the delay introduced by the smooth filtering, and
wherein the synchrophasor measurement device uses the synchrophasor to execute a steady operation and a dynamic operation of the power system.

9. The synchrophasor measurement device according to claim 8, wherein,
said filter circuit is a bandpass filter circuit being customized according to performance class, nominal system frequency of the power system and reporting rate.

10. The synchrophasor measurement device according to claim 8, wherein the system frequency estimation means determines the length of a first calculation window according to an operation status of the power system, and calculates phasor for signals obtained from the first resampling in the first calculation window.

11. The synchrophasor measurement device according to claim 8, wherein the second resampling device estimates the system frequency at a reporting time based on the updated estimated system frequency and its corresponding time tag, determines the center of a second calculation window based on the reporting time, determines the length of the second calculation window based on a reporting rate, and second-resamples the voltage and current signals with the out-of-band interference being removed based on the estimated system frequency at the reporting time in the second calculation window.

12. A system frequency measurement method for a power system for real time dynamic monitoring, operation and control of the power system, comprising:
sampling voltage and current signals in a power system based on an absolute time reference;
removing out-of-band interference in the sampled voltage and current signals;
first-resampling at least one of the voltage and current signals with the out-of-band interference being removed based on an estimated system frequency;
calculating phasor for signals obtained from the first resampling, and updating the estimated system frequency,
wherein updating the estimated system frequency comprises:
for each sample obtained from the first resampling, calculating the estimated system frequency;
smooth filtering the calculated estimated system frequency to update the estimated system frequency; and
setting a corresponding time tag for the updated estimated system frequency, the time tag compensating the delay introduced by the out-of-band interference removing and the delay introduced by the smooth filtering,

and using a synchrophasor to execute a steady operation and a dynamic operation of the power system.

13. The system frequency measurement method according to claim 12, wherein calculating phasor for signals obtained from the first resampling comprises:
determining the length of a first calculation window according to an operation status of the power system; and
calculating phasor for signals obtained from the first resampling in the first calculation window.

14. The system frequency measurement method according to claim 12, further comprising:
for a reporting time, estimating the system frequency at the reporting time based on the updated estimated system frequency and its corresponding time tag;
second-resampling the voltage and current signals with the out-of-band interference being removed based on the estimated system frequency at the reporting time; and
calculating synchrophasor for the voltage and current signals obtained from the second resampling.

15. The system frequency measurement method according to claim 14, wherein second-resampling the voltage and current signals with the out-of-band interference being removed comprising:
determining the center of a second calculation window based on the reporting time;
determining the length of the second calculation window based on a reporting rate, and
second-resampling the voltage and current signals with the out-of-band interference being removed based on the estimated system frequency at the reporting time in the second calculation window,
wherein, calculating synchrophasor for the voltage and current signals obtained from the second resampling comprises:
calculating synchrophasor for the voltage and current signals obtained from the second resampling in the second calculation window.

16. The system frequency measurement method according to claim 12, further comprising:
calculating rate of change of system frequency by using curve fitting based on the updated estimated system frequency.