Dynamic Parameter Identification With Synchronous Spectrum Fitting Technique for Sub/Supersynchronous Oscillations in Power Systems
-
摘要: “双高”新型电力系统中的次同步振荡主要由风电、光伏的电力电子设备谐振引起，且可能伴随与之频率耦合的超同步振荡，振荡传播广且变化快，因此有必要对其进行同步动态监测。该文提出一种基于同步相量频谱拟合的电力系统次/超同步振荡动态参数辨识方法。通过构建矩阵形式的方程组还原同步相量频谱的叠加特性，可准确辨识频移基波、次/超同步分量的频率、幅值和相位。该文算法相比于现有算法的优势在于，大幅将频谱分析方法必要的时间窗长缩短至200 ms，实现利用基波同步相量的次同步振荡百毫秒量级的动态同步监测；实现了同步相量数据丢失下的参数辨识；以频谱幅值作为误差权重，提升了参数辨识结果的精度。基于模拟和仿真的同步相量数据的算例结果表明，受噪声影响且数据丢失条件下，该文所提方法仍可行有效。Abstract: Subsynchronous oscillation in power systems with high penetration of renewables and inverters is mainly caused by the resonance of the power electronic equipment of wind power and photovoltaic, which may be accompanied by frequency coupling of supersynchronous oscillation. Oscillations spread widely and change quickly. Therefore, synchronous dynamic monitoring of subsynchronous oscillation is required. In this paper, a dynamic parameter identification method of power system sub/supersynchronous oscillations based on synchronous phasor spectrum fitting is proposed. By constructing a system of equations in matrix form, the superposition characteristics of the synchronous phasor spectrum are restored. It can accurately obtain the frequency, amplitude and phase of the basic component of the frequency shift and the sub/supersynchronous components. This algorithm has great advantages over the existing algorithms. First, the algorithm significantly reduces the necessary time window for spectrum analysis methods to 200 ms. It can achieve dynamic synchronous monitoring of sub synchronous oscillations in the order of one hundred milliseconds using fundamental synchronous phasors. Then, the problem of the influence of synchronous phasor data loss on parameter identification is solved. Finally, this algorithm uses spectral amplitude as error weight to improve the accuracy of parameter identification results. The results of simulated synchronous phasor data and actual simulation data show that the proposed method is still effective even under the condition of noise and data loss.
-
0. 引言
“双碳”背景下新能源(如风力发电、光伏发电)和电力电子设备大量引入电力系统，电力电子变换器间或其与电网间的次同步振荡(subsynchronous oscillation，SSO)事件频发[1-3]。在新型电力系统中，受电力电子器件控制策略的影响，SSO可能伴随与之频率耦合的超同步振荡[4]，表现出随机时变、振荡频率范围宽的特点[2]。SSO轻则降低系统输送功率，重则导致大面积切机，危及电力系统稳定运行和设备安全，因此需要对其进行全局动态监测并及时采取抑制措施[5-6]。传统依赖于故障录波器提供的瞬时采样数据的SSO辨识方法，如Prony法[7]、变分模态分解法[8]、小波变换[9]等，均只能进行SSO的本地测量而无法满足新型电力系统的需求。
广域测量系统(wide area measurement system，WAMS)在电力系统中广泛使用，其相量测量单元(phasor measurement unit，PMU)提供的基波同步相量具有同步性且动态数据上传频率高。早在2011年，冀北电网的WAMS系统监测到SSO事件[10]，调度主站基于同步相量幅值的频谱分析有效获取了SSO频率，但同时也发现次同步振荡幅值受同步相量测量算法影响而难以反推。国外学者相关研究也得到了相似的结论[11]。文献[12]利用同步相量计算和离散傅里叶变换(discrete Fourier transform，DFT)频谱分析的线性特性，首次尝试利用同步相量幅值频谱获取SSO的频率和幅值，但其为了保证较高的频率分辨率和求解精度而不得不采用长达10 s的数据窗，算法不仅实时性差且认为振荡模态10 s不变的假设过强，严重限制该方法的实用性。此外，新疆电网的WAMS系统也监测到多次20~30 Hz的次同步振荡事件[2, 9, 13-14]。文献[14]分析次/超同步振荡分量对同步相量计算的影响。文献[9]发现同步相量复数形式可有效区分次/超同步间振荡，并利用补零离散傅里叶变换(fast Fourier transform，FFT)实现次/超同步振荡参数辨识，该方法已在WAMS主站运行并实现SSO的同步实时监测，但同时也发现SSO下PMU的相量测量误差较大，并进一步影响SSO识别精度，1 s长数据窗且受噪声影响下幅值精度仅为10%左右甚至更差。
基于基波同步相量进行SSO辨识的难点在于如何使用较短的同步相量数据窗所提供的有限信息较高精度地恢复原始SSO振荡信息，现有研究主要分为频域法和时域法两大类。频域法以基于DFT频谱分析的方法为主。前文所提方法[9-10, 12]均为频谱方法。文献[15]对同步相量的复数域序列进行频谱分析，利用插值DFT和汉宁窗对频谱泄露做近似或忽略处理，将数据窗缩短至2 s。在此基础上，文献[16]所提改进插值DFT法额外分析与次同步分量耦合的超同步分量辨识方法。文献[17]将矩形窗改进为性能更好的加阻尼Rife-Vincent M阶窗函数，虽能一定程度提高辨识精度但无法进一步缩短数据窗长。文献[18]利用频谱偏移特性辨识了超同步分量。综上，基于频谱的方法通过更好地处理频谱泄露可缩短同步相量数据窗，若继续缩短数据窗，已有方法将因无法有效处理更严重的频谱泄露而失效。然而，一方面，辨识的前提是2 s内SSO模态不变，这一假设过强而限制了算法实用性；另一方面，由于新型电力系统的SSO往往是快速变化和发展的[19]，因此，百毫秒量级的次同步振荡动态监测是十分必要的，已有的频谱方法在此场景下已不再适用。
时域法以常见的模态提取法为主。模态提取法通常利用同步相量构造特定矩阵从而实现振荡参数辨识，如矩阵束算法[20]、特征值系统实现算法[13]、变分模态分解与希尔伯特变换法[21]、人工智能算法[22-23]等。与频域法相似，这些方法也需要足够长的数据窗以包含足够多的辨识信息，目前最短的数据窗为1 s。文献[24]提出同步相量轨迹拟合算法，利用非线性拟合方程组对SSO过程中时域的同步相量椭圆轨迹进行拟合，实现基于100 ms同步相量数据窗的次/超同步振荡频率和幅值辨识。该方法在100 ms数据窗条件下对测量噪声较为敏感而可靠性欠佳，需要将数据窗延长至200 ms (详见3节)。此外，在实际电力系统中存在同步相量数据丢失问题，而已有的频域法和时域法均是利用同步相量序列的时间连续性实现参数辨识的，无法处理数据丢失这一情况。
为了利用PMU数据实现百毫秒量级的次/超同步振荡动态同步监测，本文提出基于同步相量频谱拟合的电力系统次/超同步振荡参数辨识算法。
1）不同于已有的频谱分析法，该算法通过构建方程组等效频谱叠加过程，没有对频谱泄露进行忽略或近似，由此实现了超短数据窗频谱分析的参数辨识；
2）不同于轨迹拟合法，该算法用频谱提供的先验信息即频谱幅值作为误差权重，利用频谱峰值的关键信息进行参数辨识，提高了参数的求解精度；
3）由于频谱拟合法构建的拟合方程组不依赖同步相量数据的时域连续性，数据丢失时并不影响算法结构，故频谱拟合法解决了传统方法依赖同步相量时域连续性而不能处理数据丢失的问题。
1. 次/超同步振荡中的同步相量频谱拟合
本文所使用的次同步振荡过程中电力系统瞬时信号表达式即参数辨识模型如下：
$$\begin{aligned} x(t)= & x_0 \cos \left(2 \pi f_0 t+\phi_0\right)+x_{\text {sub }} \cos \left(2 \pi f_{\text {sub }} t+\phi_{\text {sub }}\right)+ \\ & x_{\text {sup }} \cos \left(2 \pi f_{\text {sup }} t+\phi_{\text {sup }}\right) \end{aligned}$$ (1) 式中：x为幅值；f为频率；φ为相位；下标“0”表示可能频移的基波正弦分量；下标“sub”表示次同步正弦分量；下标“sup”表示超同步正弦分量。
由瞬时信号计算同步相量计算的公式如下：
$$ \dot X(k) = \frac{2}{N}\sum\limits_{n = 0}^{N - 1} {x(t){{\text{e}}^{ - {\text{j}}\frac{{2{\text{π }}}}{N}n}}} $$ (2) 基波、次/超同步分量对应的同步相量均分别由正、负频率分量组成[24]，如式(3)—(5)。
$$ \left\{ {\begin{array}{*{20}{l}} {\dot X_0^ + (k) = Q({f_0}, - 1){x_0}{{\text{e}}^{{\text{j}}{\phi _0}}}{{\text{e}}^{{\text{j}}{\omega _{\text{f}}}k}}} \\ {\dot X_0^ - (k) = Q*({f_0}, + 1){x_0}{{\text{e}}^{ - {\text{j}}{\phi _0}}}{{\text{e}}^{{\text{j}}\omega _{\text{f}}^*k}}} \end{array}} \right. $$ (3) $$ \left\{ {\begin{array}{*{20}{l}} {\dot X_{{\text{sub}}}^ + (k) = Q({f_{{\text{sub}}}}, - 1){x_{{\text{sub}}}}{{\text{e}}^{{\text{j}}{\phi _{{\text{sub}}}}}}{{\text{e}}^{{\text{j}}{\omega _{{\text{sub}}}}k}}} \\ {\dot X_{{\text{sub}}}^ - (k) = Q*({f_{{\text{sub}}}}, + 1){x_{{\text{sub}}}}{{\text{e}}^{ - {\text{j}}{\phi _{{\text{sub}}}}}}{{\text{e}}^{{\text{j}}\omega _{{\text{sub}}}^*k}}} \end{array}} \right. $$ (4) $$ \left\{ {\begin{array}{*{20}{l}} {\dot X_{\sup }^ + (k) = Q*({f_{\sup }}, + 1){x_{\sup }}{{\text{e}}^{ - {\text{j}}{\phi _{\sup }}}}{{\text{e}}^{{\text{j}}\omega _{\sup }^*k}}} \\ {\dot X_{\sup }^ - (k) = Q({f_{\sup }}, - 1){x_{\sup }}{{\text{e}}^{{\text{j}}{\phi _{\sup }}}}{{\text{e}}^{{\text{j}}{\omega _{\sup }}k}}} \end{array}} \right. $$ (5) 式中：“*”表示共轭；函数Q(f, l)表示同步相量计算。
$$ Q(f,l) = \frac{1}{N}\sum\limits_{n = 0}^{N - 1} {{{\text{e}}^{{\text{j}}\frac{{2\pi n}}{N}(\frac{f}{{{f_{\text{N}}}}} + l)}}} $$ (6) 且有：
$$ {\omega _{\text{f}}} = \frac{{2{\text{π }}{f_0}}}{{{f_{\text{s}}}}},\;\;{\omega _{{\text{sub}}}} = \frac{{2{\text{π }}{f_{{\text{sub}}}}}}{{{f_{\text{s}}}}},\;\;{\omega _{{\text{sup}}}} = \frac{{2{\text{π }}{f_{{\text{sup}}}}}}{{{f_{\text{s}}}}} $$ (7) 式中：fs为PMU中同步相量的上传频率，我国工程实际中多为100 Hz；ωf、ωsub、ωsup分别为f0、fsub和fsup对应的相对角频率。由于次/超同步振荡之间存在频率耦合关系[16, 24]，即fsub+fsup=2f0，考虑到实际电力系统中f0非常接近fN，故有${\omega _{{\text{sub}}}} = \omega _{{\text{sup}}}^*$，统一以ωs表示。因此次、超同步相量的正、负频率分量可分别合成振荡分量的正、负频率分量。
$$\left\{\begin{array}{c} \dot{X}_{\mathrm{s}}^{+}(k)=\left[Q\left(f_{\text {sub }},-1\right) x_{\text {sub }} \mathrm{e}^{\mathrm{j} \phi_{\text {sub }}}+\right. \\ \left.\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;Q^*\left(f_{\text {sup }},+1\right) x_{\text {sup }} \mathrm{e}^{-\mathrm{j} \phi_{\text {sup }}}\right] \mathrm{e}^{\mathrm{j} \omega_{\mathrm{s}} k} \\ \dot{X}_{\mathrm{s}}^{-}(k)=\left[Q^*\left(f_{\text {sub }},+1\right) x_{\text {sub }} \mathrm{e}^{-\mathrm{j} \phi_{\text {sub }}}+\right. \\ \left.\;\;\;\;\;\;\;\;\;\;\;\;\;Q\left(f_{\text {sup }},-1\right) x_{\text {sup }} \mathrm{e}^{\mathrm{j} \phi_{\text {sup }}}\right] \mathrm{e}^{\mathrm{j} \omega_{\mathrm{s}}^* k} \end{array}\right.$$ (8) 综上，同步相量$ \dot X(k) $由基波正频率$ \dot X_0^ + (k) $、基波负频率$ \dot X_0^ - (k) $、振荡正频率$ \dot X_{\text{s}}^ + (k) $、振荡负频率$ \dot X_{\text{s}}^ - (k) $4个分量组成，如式(9)，分别对应相对角频率ωf、$\omega _{\text{f}}^*$、ωs、$\omega _{\text{s}}^*$。
$$ \dot X(k) = \dot X_0^ + (k) + \dot X_0^ - (k) + \dot X_{\text{s}}^ + (k) + \dot X_{\text{s}}^ - (k) $$ (9) 以F表示对同步相量序列的DFT频谱分析，使用的同步相量序列长度为M=2K+1，频谱结果为$\dot F(m)$，m=0, 1, ···, M−1。
$$ \dot F(m) = F[\dot X(k)] = \frac{1}{M}\sum\limits_{k = - K}^K {\dot X(k){w^{mk}}} ,\;\;w = {{\text{e}}^{ - {\text{j}}\frac{{2{\text{π }}}}{M}}} $$ (10) 由于该频谱分析是线性的，则分别对应式(9)中同步相量的4个旋转分量，同步相量的频谱分析结果$ \dot F(m) $也是由基波正频谱$ \dot F_0^ + (m) $、基波负频谱$ \dot F_0^ - (m) $、振荡正频谱$ \dot F_{\text{s}}^ + (m) $、振荡负频谱$ \dot F_{\text{s}}^ - (m) $组成。
$$ \dot F(m) = \dot F_0^ + (m) + \dot F_0^ - (m) + \dot F_{\text{s}}^ + (m) + \dot F_{\text{s}}^ - (m) $$ (11) $$ \left\{ {\begin{array}{*{20}{l}} {\dot F_0^ + (m) = D(m,{\omega _{\text{f}}})\dot X_0^ + (0)} \\ {\dot F_0^ - (m) = D(m,\omega _{\text{f}}^{\text{*}})\dot X_0^ - (0)} \\ {\dot F_{\text{s}}^ + (m) = D(m,{\omega _{\text{s}}})\dot X_{\text{s}}^ + (0)} \\ {\dot F_{\text{s}}^ - (m) = D(m,\omega _{\text{s}}^*)\dot X_{\text{s}}^ - (0)} \end{array}} \right. $$ (12) 其中，引入D(ω, m)函数表示同步相量的频谱分析过程，其表达式为
$$ D(m,\omega ) = \frac{1}{M}\sum\limits_{k = - {\text{K}}}^{\text{K}} {{w^{mk}}{{\text{e}}^{\omega k}}} $$ (13) 基于上述分析和推导，在某一频点m处，可建立相应的频谱拟合方程如下：
$$\hat{\dot{F}}(m)+\dot{\varepsilon}(m)=\dot{F}(m)$$ (14) 式中：$\hat{\dot{F}}(m)$为拟合的同步相量频谱，具体表达式为式(12)的4个分量之和；$ \dot F(m) $为实际的同步相量频谱，可根据同步相量数据序列直接由式(10)计算；$ \dot \varepsilon (m) $为拟合频谱和实际频谱的误差。
上述频谱拟合方程在所有m取值上均成立，m=0, 1, ···, M−1，也就是说，全频段所有频点的频谱拟合方程可组成一个方程组，则这一全频段频谱拟合方程组可写为如式(15)所示的矩阵形式。
频谱拟合方程组中的待定量共有6个：仅包含频率信息的2个实部为零的纯虚数ωf和ωs(0 < ωs < jπ, 0 < ωs < jπ)，同时包含幅值和相位信息的4个复数$\dot X_0^ + (0)$、$\dot X_0^ - (0)$、$\dot X_{\text{s}}^ + (0)$、$\dot X_{\text{s}}^ - (0)$。一方面，从频谱的叠加特性看，频谱拟合方程组可写为式(15)的矩阵形式；另一方面，对于未知数而言，其是一个非线性超定方程组，利用数值解法对其进行求解即可求得待定量。
本文所提的频谱拟合法是频域的频谱分析方法和时域的轨迹拟合方法的优势组合。一方面，受频谱分析法[15-17]的启发，频谱拟合法仍从频谱出发，充分利用频谱分析能突出主要振荡模态便于特征提取的优势；另一方面，受轨迹拟合法[24]的启发，频谱拟合法通过构建超定非线性方程组对同步相量频谱进行全谱拟合，而并未对频谱泄露进行忽略或近似。由此，频谱拟合法在保证参数辨识精度的同时，缩短了数据窗长，并提高了算法的实时性。
$$ \left[ {\begin{array}{*{20}{c}} {D(0,{\omega _{\text{f}}})}&{D(0,\omega _{\text{f}}^{\text{*}})}&{D(0,{\omega _{\text{s}}})}&{D(0,\omega _{\text{s}}^*)} \\ {D(1,{\omega _{\text{f}}})}&{D(1,\omega _{\text{f}}^{\text{*}})}&{D(1,{\omega _{\text{s}}})}&{D(1,\omega _{\text{s}}^*)} \\ \vdots & \vdots & \vdots & \vdots \\ {D(M - 2,{\omega _{\text{f}}})}&{D(M - 2,\omega _{\text{f}}^{\text{*}})}&{D(M - 2,{\omega _{\text{s}}})}&{D(M - 2,\omega _{\text{s}}^*)} \\ {D(M - 1,{\omega _{\text{f}}})}&{D(M - 1,\omega _{\text{f}}^{\text{*}})}&{D(M - 1,{\omega _{\text{s}}})}&{D(M - 1,\omega _{\text{s}}^*)} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\dot X_0^ + (0)} \\ {\dot X_0^ - (0)} \\ {\dot X_{\text{s}}^ + (0)} \\ {\dot X_{\text{s}}^ - (0)} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {\dot \varepsilon (0)} \\ {\dot \varepsilon (1)} \\ \vdots \\ {\dot \varepsilon (M - {\text{2}})} \\ {\dot \varepsilon (M - 1)} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\dot F(0)} \\ {\dot F(1)} \\ \vdots \\ {\dot F(M - 2)} \\ {\dot F(M - 1)} \end{array}} \right] $$ (15) 2. 同步相量频谱拟合的方法和优势
2.1 超同步振荡辨识的可行性
利用上传频率为100 Hz的同步相量序列可辨识55~95 Hz的超同步振荡参数。传统采样定理中采样对象为实数域的瞬时信号，其频谱左右半幅呈共轭对称性，仅有一半频谱能提供有效信息，故可辨识的频率范围为采样频率的一半。而本文频谱分析的对象为PMU提供的复数域的基波同步相量数据，其复数域频谱的左右半幅是非共轭的，因此同步相量频谱可辨识的频率范围等于采样频率，故利用同步相量复频谱可辨识的频率范围为0~100 Hz。由于新能源发电引起的SSO往往伴随耦合的超同步振荡。因此，有必要考虑次/超同步之间的耦合关系[16, 18, 23-24, 26-27]并辨识超同步振荡参数，仅利用同步相量幅值信息则受可辨识频率范围所限而并不能辨识超同步[10, 12, 21]。
2.2 理想条件下的零误差辨识与高实时性
本文所提的频谱拟合法与传统频谱分析法不同，并未对频谱泄露进行近似或忽略处理，并由此实现了利用比传统频谱分析法更短的数据窗进行参数辨识，具体分析如下。
传统频谱分析法通过使用汉宁窗[15-16]或RV-m窗[17]对此处频谱进行处理，实现可将$ \dot F_{\text{s}}^ - (m) $、$ \dot F_0^ + (m) $、$ \dot F_0^ - (m) $忽略，均为频谱泄露而较小，进而得到此m点处有$ \dot F(m) \approx \dot F_{\text{s}}^ + (m) $(需要注意的是，此式为示意性表达，实际具体公式随方法不同而有所差异)，如图 1(a)所示，进而仅利用频谱峰值信息即可实现参数辨识。但是，当数据窗进一步缩短时，如图 1(b)所示，即M取值显著变小，D(ω, m)函数的频率筛选特性迅速变差导致频谱泄露问题凸显。这种情况下，传统频谱方法对频谱泄露做近似或忽略的处理方法不再有效，必须将频谱泄露对主导分量频谱的影响考虑在内。综上可见，传统频谱分析法的数据窗必须足够长才能包含足够多的动态信息以实现较高的频率分辨率，进而保证对频谱近似或忽略处理时仍保留足量信息以确保辨识精度。
频谱拟合法的核心思路是通过构建如式(15)的矩阵形式的方程组还原全频段的频谱叠加特性，通过拟合所求得的结果是对频谱结果的全频段最小二乘拟合。如图 2所示，频谱拟合结果所得的两个频率ωf、ωs和4个分量$\dot X_0^ + (0)$、$\dot X_0^ - (0)$、$\dot X_{\text{s}}^ + (0)$、$\dot X_{\text{s}}^ - (0)$是对$ \dot F(m) $在所有M个频点处的拟合，即m=0, 1, ···, M−1，在任意频点上的频谱不区分主导分量与否，而是以更精确的方程组描述频谱叠加特性，进而替代了传统方法在频谱峰值对频谱泄露的忽略或近似处理，并由此实现了对频谱信息更加高效的利用。综上，频谱拟合法的数据窗显著缩短至百毫秒量级，相较于传统频谱法秒级的数据窗长，可仅利用其十分之一的动态信息进行参数辨识。
![]()
图 2 频谱叠加特性及频谱拟合的构成Figure 2. Spectrum superposition characteristic and composition of the fitting spectrum在理想条件下，即瞬时信号中仅有基波分量和一对频率耦合的次、超同步振荡分量，其模型结构与式(1)一致，并且作为辨识数据的同步相量的算法与式(2)一致时，则频谱拟合法可实现零误差辨识。此时，由于从式(1)瞬时值模型至式(15)频谱拟合方程的推导过程没有任何近似环节，则式(15)所求得的拟合结果可反推出与式(1)完全一致的参数，即零误差辨识。相对地，传统频谱方法[15-17]因在频谱峰值处对频谱泄露做了近似处理，导致这些方法即便在理想条件下仍不是零误差辨识，缩短数据窗将导致参数辨识误差显著增加而算法不可用，如3节所述。在实际条件下，式(1)中的瞬时值在采样、传输过程中会受到以高斯噪声为主的各种噪声干扰，使得PMU提供的同步相量数据中将包含噪声[25]。噪声条件由于对基于基波同步相量的次同步振荡参数辨识影响较大，是目前相关研究所考虑的主要干扰因素[15-17, 20, 26]，PMU的瞬时值采样信噪比通常为45 dB[15]。为了保证噪声条件下频谱拟合算法的高鲁棒性，本文采用了数据窗长200 ms、上传频率100 Hz的同步相量数据进行频谱拟合参数辨识。
综上可见，得益于对同步相量的复数域频谱特性更精细地数学解析刻画，本文所提频谱拟合法相较于传统频谱方法大大缩短了同步相量的数据窗长，大幅提高了算法的实时性。
2.3 数据丢失条件下的频谱拟合
在实际工程应用中，PMU设备故障、通信中断等原因将导致个别同步相量数据点丢失。现有的基于同步相量参数辨识方法均使用时间连续的同步相量数据序列而无法处理数据丢失这一特殊情况。本文所提的频谱拟合算法结构并不受此影响，可实现数据丢失条件下的参数辨识。
式(10)所示的同步相量序列的频谱计算过程可写为如下矩阵形式：
$$\dot{\boldsymbol{F}}=\boldsymbol{W} \cdot \dot{\boldsymbol{X}}$$ (16) 式中：$ \dot{\boldsymbol{F}} $表示实际频谱；矩阵W表示DFT频谱的计算；w为式(10)中的计算因子；向量$ \dot{\boldsymbol{X}} $表示同步相量序列，其中，
$$ \left\{ \begin{gathered} {\boldsymbol{W}} = \left[ {\begin{array}{*{20}{c}} 1&1& \cdots &1 \\ {{w^{ - K}}}&{{w^{ - K{\text{ + 1}}}}}& \cdots &{{w^K}} \\ {{w^{ - 2K}}}&{{w^{ - 2K{\text{ + 2}}}}}& \cdots &{{w^{2K}}} \\ \vdots & \vdots & \ddots & \vdots \\ {{w^{ - K(M - 1)}}}&{{w^{( - K{\text{ + 1}})(M - 1)}}}& \cdots &{{w^{K(M - 1)}}} \end{array}} \right] \hfill \\ {\boldsymbol{X}} = {[\begin{array}{*{20}{c}} {\dot X( - K)}&{\dot X( - K + 1)}&{\dot X( - K + 2)}&{\dot X(K)} \end{array}]^{\text{T}}} \hfill \\ \end{gathered} \right. $$ (17) 相似地，根据式(13)，拟合频谱的特性系数矩阵D可写成W矩阵和R矩阵相乘，即$\boldsymbol{D}=\boldsymbol{W} \cdot \boldsymbol{R}$：
$$ {\boldsymbol{R}} = \left[ {\begin{array}{*{20}{c}} {{{\text{e}}^{ - K{\omega _{\text{f}}}}}}&{{{\text{e}}^{ - K\omega _{\text{f}}^{\text{*}}}}}&{{{\text{e}}^{ - K{\omega _{\text{s}}}}}}&{{{\text{e}}^{ - K\omega _{\text{s}}^{\text{*}}}}} \\ {{{\text{e}}^{( - K{\text{ + 1}}){\omega _{\text{f}}}}}}&{{{\text{e}}^{( - K{\text{ + 1}})\omega _{\text{f}}^{\text{*}}}}}&{{{\text{e}}^{( - {\text{K + 1}}){\omega _{\text{s}}}}}}&{{{\text{e}}^{( - K{\text{ + 1}})\omega _{\text{s}}^{\text{*}}}}} \\ \vdots & \vdots & \vdots & \vdots \\ {{{\text{e}}^{K{\omega _{\text{f}}}}}}&{{{\text{e}}^{K\omega _{\text{f}}^{\text{*}}}}}&{{{\text{e}}^{{\text{K}}{\omega _{\text{s}}}}}}&{{{\text{e}}^{K\omega _{\text{s}}^{\text{*}}}}} \end{array}} \right] $$ (18) 式中矩阵R表示由旋转分量e构成，R的各列表示同步相量4个分量分别以角频率±ωf、±ωs旋转，各行表示某一时刻同步相量4个分量旋转后的叠加。故拟合频谱可简写为
$$\hat{\dot{\boldsymbol{F}}}=\boldsymbol{W} \cdot \boldsymbol{R} \cdot \dot{\boldsymbol{X}}_{0 \mathrm{~s}}$$ (19) 当同步相量序列中某一数据点丢失时，拟合频谱和实际频谱做出相应改变即可保证频谱拟合方程的成立，进而实现丢数据条件下的参数辨识，其原理如图 3所示。以同步相量序列丢失第2个数据点为例：对于式(15)右边的实际频谱，应删除同步相量序列中的数据丢失点即$ \dot{\boldsymbol{X}} $矩阵的第2行及W中第2个数据点对应的频谱变换作用即W的第2列；对于式(15)左边的拟合频谱，同样删除W矩阵的第2列，并删除相应的旋转特性即R矩阵的第2行。由此可见，尽管丢失数据导致方程左右两边频谱因缺项变为缺项频谱，但丢数据后频谱维数不会改变，拟合方程仍然成立，仍可进行参数辨识。而因缺项频谱并不能提供所需的频谱峰值信息，传统频谱分析法[15-17]并不能处理数据丢失这一情况。
当发生多个数据丢失时，与丢失1个数据点类似，需删除所有丢失数据点及其有关特性、对应的频谱变换作用，即对丢失多个数据点的缺项频谱进行拟合，从而实现参数辨识。具体地，设L分别表示丢失数据点的集合，并以下标“L”表示数据丢失，则数据丢失时频谱拟合方程变为：
$$\boldsymbol{W}_{\mathrm{L}} \cdot \boldsymbol{R}_{\mathrm{L}} \cdot \dot{\boldsymbol{X}}_{05}+\dot{\boldsymbol{E}}_{\mathrm{L}}=\boldsymbol{W}_{\mathrm{L}} \dot{\boldsymbol{X}}_{\mathrm{L}}$$ (20) ![]()
图 3 丢失数据点时频谱拟合算法原理示意图Figure 3. Schematic diagram of the spectrum fitting algorithm when data loss occurs$$\boldsymbol{W}_{\mathrm{L}}=\left\{\begin{array}{cc} \boldsymbol{W}_{\mathrm{L}}(:, m)=\varnothing & , m \in L \\ \boldsymbol{W}_{\mathrm{L}}(:, m)=\boldsymbol{W}(:, m) & , m \notin L \end{array}\right.$$ (21) $$\boldsymbol{R}_{\mathrm{L}}=\left\{\begin{array}{cl} \boldsymbol{R}_{\mathrm{L}}(m,:)=\varnothing & , m \in L \\ \boldsymbol{R}_{\mathrm{L}}(m,:)=\boldsymbol{R}(m,:) & , m \notin L \end{array}\right.$$ (22) $$\dot{\boldsymbol{E}}_{\mathrm{L}}=\left\{\begin{array}{cc} \dot{\boldsymbol{E}}_{\mathrm{L}}(m,:)=\varnothing & , m \in L \\ \dot{\boldsymbol{E}}_{\mathrm{L}}(m,:)=\dot{\boldsymbol{E}}(m,:) & , m \notin L \end{array}\right.$$ (23) $$\dot{\boldsymbol{X}}_{\mathrm{L}}=\left\{\begin{array}{cc} \dot{\boldsymbol{X}}_{\mathrm{L}}(m)=\varnothing & , m \in L \\ \dot{\boldsymbol{X}}_{\mathrm{L}}(m)=\dot{\boldsymbol{X}}(m) & , m \notin L \end{array}\right.$$ (24) 式中：$\varnothing$为空集，(: , m)表示第m列；(m, : )表示第m行；WL为缺项的频谱变换方程；RL为缺项的旋转特性矩阵；XL为缺项的同步相量序列。
此外，若出现[1, 3, ···, 2k−1, ···, 21]$ \subseteq $L和[2, 4, ···, 2k, ···, 20]$ \subseteq $L这两种特殊情况，根据奈奎斯特采样定理，无法识别出SSO参数。但实际电力系统的数据丢失往往是随机的，出现上述情况的概率并不大。实际应用中，测量噪声会降低参数辨识结果的精度，故丢包率不能过大。工程应用中的PMU信噪比一般在45 dB左右，在40 dB噪声条件下，频谱拟合法可容许的最大丢包率为20%，具体结果见3节。
2.4 误差加权的方程求解和参数辨识
为充分利用频谱峰值提供的有效信息以进一步提升所提频谱拟合法的参数辨识精度，在式(15)的基础上，本文采用不同频点的频谱幅值作为频谱拟合方程的误差加权系数。
具体地，考虑权重的拟合目标函数如下，
$$\min \sum\limits_{m=0}^{M-1}|\dot{F}(m)| \cdot\|\dot{\varepsilon}(m)\|_2=\min \left(\operatorname{norm}_2\left(\boldsymbol{F}_{\mathrm{w}} \dot{\boldsymbol{E}}\right)\right)^2$$ (25) 式中：$\dot F(m)$为在频点m处的同步相量频谱；FW为各权重系数组成的对角矩阵，有：
$$\boldsymbol{F}_{\mathrm{w}}=\operatorname{diag}(|\dot{F}(0)|, \cdots,|\dot{F}(M-1)|)$$ (26) 则可将各权重系数写为矩阵形式并加入到式(15)中，得到加权后的频谱拟合方程为
$$\boldsymbol{F}_{\mathrm{w}} \hat{\dot{\boldsymbol{F}}}+\boldsymbol{F}_{\mathrm{w}} \dot{\boldsymbol{E}}=\boldsymbol{F}_{\mathrm{w}} \dot{\boldsymbol{F}}$$ (27) 由于加权系数矩阵FW可根据同步相量序列直接计算得到且计算量非常小，求解式(27)所示的加权频谱拟合方程的复杂度与式(15)的非加权频谱拟合方程完全一致而不会显著增加计算量。
以不同频率点处的频谱幅值作为误差权重将增加频谱尖峰、削弱非尖峰的影响，由于同步相量的频谱峰值分别对应于基波和振荡模态，以加权频谱误差最小为目标进行拟合求解，增强了基波、次/超同步振荡分量的影响，削弱了频谱中其他分量的影响。换言之，频谱拟合法充分利用频谱峰值所含的关键信息，可保证关键振荡模态的频谱误差最小，建立频谱误差与参数间的关联关系，有利于进一步提高辨识结果的精度。
相对地，轨迹拟合法[24]以轨迹误差最小为目标对时域的同步相量椭圆轨迹进行拟合。即使振荡模态不变，由于同步相量椭圆轨迹在不同时刻的相对位置不同，时域的轨迹误差并不与振荡模态误差直接相关。3节的结果表明，尤其是在相位参数辨识方面，频谱拟合法具有独特的优势。
最终，根据式(25)和(27)求得的ωf、$ \dot X_0^ + (0) $和$ \dot X_0^ - (0) $，可进一步求得基波分量的频率、幅值和相位：
$$ {f_0} = \left\{ {\begin{array}{*{20}{l}} {\frac{{{\omega _{\text{f}}}{f_{\text{s}}}}}{{{\text{j}}2{\text{π }}}}}&{,||\dot X_0^ + (0)|{|_2} < ||\dot X_0^ - (0)|{|_2}} \\ { - \frac{{{\omega _{\text{f}}}{f_{\text{s}}}}}{{{\text{j}}2{\text{π }}}}}&{,||\dot X_0^ + (0)|{|_2} > ||\dot X_0^ - (0)|{|_2}} \end{array}} \right. $$ (28) $$\left\{\begin{array}{l} Q\left(f_0,-1\right) x_0 \mathrm{e}^{\mathrm{j} \phi_0}=\dot{X}_0^{+}(0),\left\|\dot{X}_0^{+}(0)\right\|_2<\left\|\dot{X}_0^{-}(0)\right\|_2 \\ Q\left(f_0,-1\right) x_0 \mathrm{e}^{\mathrm{j} \phi_0}=\dot{X}_0^{-}(0),\left\|\dot{X}_0^{+}(0)\right\|_2 \geq\left\|\dot{X}_0^{-}(0)\right\|_2 \end{array}\right.$$ (29) 根据求得的ωs、$ \dot X_{\text{s}}^{\text{ + }}({\text{0}}) $和$ \dot X_{\text{s}}^ - ({\text{0}}) $，通过式(30)—(32)，可得求得次/超同步振荡分量的频率、幅值和相位：
$$ {f_{{\text{sub}}}} = \frac{{{\omega _{\text{s}}}{f_{\text{s}}}}}{{{\text{j}}2{\text{π }}}} $$ (30) $$ {f_{{\text{sup}}}} = 2{f_{\text{N}}} - {f_{{\text{sub}}}} $$ (31) $$\left[\begin{array}{ll} Q\left(f_{\text {sub }},-1\right) & Q^*\left(f_{\text {sup }},+1\right) \\ Q\left(f_{\text {sub }},+1\right) & Q^*\left(f_{\text {sup }},-1\right) \end{array}\right]\left[\begin{array}{c} x_{\text {sub }} \mathrm{e}^{\mathrm{j} \phi_{\text {sub }}} \\ x_{\text {sup }} \mathrm{e}^{-\mathrm{j} \phi_{\text {sup }}} \end{array}\right]=\left[\begin{array}{c} \dot{X}_{\mathrm{s}}^{+}(0) \\ \dot{X}_{\mathrm{s}}^{*-}(0) \end{array}\right]$$ (32) 综上，频谱拟合法在采用频谱幅值这一先验信息作为误差权重，建立了频谱误差与振荡参数之间的关联关系，提高了辨识结果的精度；加强含有关键信息的频谱峰值的影响，保证了噪声条件下频谱拟合算法的高鲁棒性。
3. 模拟PMU数据的验证和特性分析
为了验证频谱拟合算法的正确性和有效性，本文分别采用两种PMU同步相量数据进行验证分析。一种是由给定的瞬时值计算而来的模拟PMU数据，旨在分析在振荡参数全工况大范围内取值时的参数辨识精度，以验证算法正确性和完备性；另一种是仿真SSO场景下由仿真所得瞬时值计算而来的仿真PMU数据，旨在验证典型的接近实际SSO振荡过程中的算法有效性，两者互为补充。
本节为模拟PMU数据的验证和特性分析，在理想条件下验证频谱拟合算法的零误差辨识特性，并在此基础上叠加噪声和数据丢失的影响，验证算法鲁棒性。同时，本节采用改进插值DFT法[16]和轨迹拟合法[24]作为对比，以验证本文方法的优越性。本节使用的同步相量数据是由式(1)所示的瞬时信号模型通过式(2)的DFT同步相量算法计算得到，即满足理想条件。式(2)中瞬时值采样频率为1.6 kHz，数据窗长为20 ms即一个整周期。具体地，式(1)瞬时信号中各个参数设定如下：系统额定频率为50 Hz，基波幅值为100 A，基波相位为2.5 rad，基波频率分别取在49.5~50.5Hz范围内以0.2 Hz间隔变化，分别以次/超同步振荡的频率、幅值和相位作为变量。频谱拟合算法选用200 ms数据窗内的基波同步相量进行参数辨识，以辨识结果与真实值相对误差为指标进行分析，详细的误差统计和对比分析结果如表 1所示。由于测试结果表明次同步振荡频率的变化对辨识结果误差影响最大，本节仅以次同步振荡频率变化时的结果图为代表，见图 4—6。
表 1 同步相量频谱拟合算法与其他算法的性能比较Table 1. Performance comparison between synchronous phasor spectrum fitting algorithm and other algorithms各分量辨识结果相对误差① 噪声条件SNR/dB 数据丢失率/% 频谱拟合 文献[24] (轨迹拟合) 文献[16] (改进的插值DFT) 0.2 s 0.2 s③ 0.2 s 1 s 2 s $\frac{{|{f_0} - {{\hat f}_0}|}}{{{f_0}}}$ ∞② 0 0 10−15 10−3 10−4 10−4 40 0 2×10−3 10−15 10−3 10−4 10−4 20 2×10−3 — — — — $ \frac{{|{x_0} - {{\hat x}_0}|}}{{{x_0}}} $ ∞ 0 10−14 10−14 10−2 10−3 10−2 40 0 6×10−3 10−2 10−2 10−2 10−2 20 10−2 — — — — $ \frac{{|{\phi _0} - {{\hat \phi }_0}|}}{{{\phi _0}}} $ ∞ 0 10−14 10−14 10−2 10−2 10−2 40 0 3×10−2 10−1 10−2 10−2 10−2 20 3×10−2 — — — — $ \frac{{|{f_{{\text{sub}}}} - {{\hat f}_{{\text{sub}}}}|}}{{{f_{{\text{sub}}}}}} $ ∞ 0 10−14 10−14 10−1 10−3 10−4 40 0 5×10−3 10−2 0.5 > 1 10−3 20 10−2 — — — —/ $ \frac{{|{x_{{\text{sub}}}} - {{\hat x}_{{\text{sub}}}}|}}{{{x_{{\text{sub}}}}}} $ ∞ 0 10−12 10−12 > 1 10−1 10−2 40 0 10−2 5×10−2 > 1 1 10−1 20 2×10−2 — — — —/ $ \frac{{|{\phi _{{\text{sub}}}} - {{\hat \phi }_{{\text{sub}}}}|}}{{{\phi _{{\text{sub}}}}}} $ ∞ 0 10−12 10−12 > 1 0.15 10−1 40 0 6×10−3 10−1 > 1 > 1 10−1 20 2×10−2 — — — — ①本表所示结果为不同fsub下的辨识结果相对误差最大量级。不同算法的测试条件均为fsub在[5, 45]Hz范围内变化，fsup、xsup、ϕsup的结果误差与fsub、xsub、ϕsub的结果误差相近而未列出。
②SNR为信噪比(signal-noise ratio)，SNR=∞表示没有噪声。
③轨迹拟合在噪声条件下使用0.1 s数据窗的相位辨识结果误差比使用0.2 s数据窗的误差大，相位辨识结果误差大于10%，不满足精度要求。![]()
图 4 x0=100，xsub=30，xsup=20，ϕ0=2.5，ϕsub=3，ϕsup=−3和fsub在[5, 45]Hz范围内变化时的各分量的频率、幅值和相位估计误差Figure 4. Estimation errors with x0=100, xsub=30, xsup=20, ϕ0=2.5, ϕsub=3, ϕsup=−3 and fsub varying in the range [5, 45]Hz![]()
图 5 x0=100，xsub=30，xsup=20，ϕ0=2.5，ϕsub=3，ϕsup=−3和fsub在[5, 45]Hz范围内变化时的各分量的频率、幅值和相位估计误差(加40 dB的噪声)Figure 5. Estimation errors with x0=100, xsub=30, xsup=20, ϕ0=2.5, ϕsub=3, ϕsup=−3 and fsub varying in the range [5, 45]Hz (with 40 dB noise)![]()
图 6 40 dB噪声条件下，丢包率为20%时各分量的频率、幅值和相位估计误差Figure 6. Estimation errors when data loss rate is 20% under 40 dB noise condition1）理想条件下：即瞬时值如式(1)，有频移基波和一对次/超同步分量、无噪声和数据丢失时，表 1和图 4中，基波、次/超同步分量的频率、幅值和相位的相对误差均小于10−13，即计算机浮点数舍入误差量级，可见频谱拟合法可实现零误差辨识。不同次/超同步振荡分量的幅值和相位辨识结果均与图 4相类似，辨识结果相对误差均小于10−13，此处不再赘述。
由表 1所示理想条件下的不同算法对比结果可以看出，改进的插值DFT法由于对频谱泄露进行了近似或忽略导致即使在理想条件下辨识结果仍存在误差，其使用的最短数据窗为2 s，继续缩短数据窗至1 s甚至是0.2 s时因辨识误差过大而无法进行参数辨识。与此相比，轨迹拟合法和本文所提频谱拟合法，由于未做任何忽略或近似，在理想条件下可实现振荡参数的零误差辨识。
2）噪声的影响：由于瞬时值采样和测量噪声对基于同步相量的SSO参数辨识影响显著[15-17, 20, 26]，为了研究频谱拟合法噪声条件下的鲁棒性，在理想条件下瞬时信号中添加高斯白噪声并重复测试，结果如表 1和图 5所示，其中信噪比设定为40 dB[15]。由于文献[24]并未展示0.1 s数据窗长和噪声条件下的相位辨识结果，但实际复现结果表明，该条件下各分量相位辨识结果相对误差达50%量级而不可用，故本文以0.2 s数据窗长下轨迹拟合法结果作为对比。表 1的对比结果表明，在噪声影响下，频谱拟合法相较于其他方法鲁棒性更高。具体地，改进插值DFT法在1 s数据窗下的次同步分量幅值误差达到10%量级，在0.2 s数据窗下的辨识结果相对误差甚至大于100%，显然2 s数据窗已是作为传统频谱分析法代表的插值DFT法的极限最短窗长；轨迹拟合法尽管仅用0.2 s窗长时的频率和幅值辨识结果在1%，但唯独次同步分量相位的相对误差达到10%量级，由于扰动源定位等高级应用需要相位信息[13]，则轨迹拟合法在此场景下也几乎不可用；频谱拟合法相位辨识结果的精度明显提高，虽次同步频率辨识结果的精度有所降低，但已满足对精度的要求。
3）数据丢失的影响：理想且无噪声条件下，当丢包率小于75%时，频谱拟合法仍可实现零误差辨识。考虑实际情况下噪声的影响，在40 dB SNR时，分别取不同的丢包率进行测试，得出频谱拟合法可容许的最大丢包率为20%，其结果如图 6所示，各参数相对误差最大不超过4×10−2。表 1对比结果表明，无论是轨迹拟合法还是改进的插值DFT法，均没有就数据丢失条件下的参数辨识问题给出相应的解决方案。
4）算法计算时间：频谱拟合法相较于轨迹拟合法的一个优势在于可利用频谱峰值对应的粗略频率范围这一先验信息简化数值求解难度；而同步相量的时域轨迹无法直接提供次同步频率的粗略范围。在matlab中测试算法，计算机配置为Intel i5-12400U CPU @ 4×2.5GHz和16GB RAM，当N=21时(即取200 ms数据窗)，使用matlab自带fslove函数数值求解，频谱拟合法和轨迹拟合法每次辨识计算耗时分别为0.05和0.2 s。但频谱拟合法中拟合方程的数值求解仍有大幅优化空间，与轨迹拟合方程中的系数矩阵中个元素均为幅值为1的复数所不同[24]，频谱拟合法中的系数矩阵D因具有频率筛选特性而是一个稀疏矩阵，后续研究中可以利用其稀疏性进一步优化求解过程以缩短计算时间。
4. 仿真PMU数据的验证
新型电力系统的SSO表现出随机时变的特点，因此，现场实际情况下信号的参数可能会发生变化。为了验证该算法的工程实用性，证明其可对振荡实现百毫秒量级的动态监测，在完整数据和数据丢失两种情况下，采用仿真PMU数据分别进行验证分析。本节使用的SSO仿真模型同文献[24]，分别仿真振荡发散且只有次同步振荡、振荡不发散且次同步振荡过程中伴有超同步振荡两个不同场景。
4.1 完整数据情况下
场景一和场景二均采用A相瞬时电流信号，采样频率分别为2和20 kHz，各自截取10和6 s的振荡过程进行分析，如图 7所示，场景一的振荡是发散的，场景二的振荡是非发散的。瞬时值的频谱分析结果如图 8所示，场景一只有次同步振荡，场景二的次同步振荡过程中伴有超同步振荡。受频谱泄露影响严重，此结果仅做参考。用改进插值DFT法、轨迹拟合法和Prony算法的计算结果作为频谱拟合算法动态辨识结果的误差性能评价指标。频谱拟合法、轨迹拟合法和Prony算法均使用200 ms数据窗，改进插值DFT法使用2 s数据窗。
1）场景一下次同步振荡结果：图 8(a)中次同步频率点附近频谱拟合法辨识结果的放大及其他算法的辨识结果见图 9。其中不同颜色代表不同时间段内的辨识结果，颜色渐变表示时间的发展，以频谱拟合法为例，其对应的50组辨识结果以红色开始，黄色结束。次同步频率的辨识结果表明，本文提出的频谱拟合法、轨迹拟合法以及Prony算法所得的次同步频率在4.45~4.49 Hz范围内变化，且频率呈发散趋势，改进的插值DFT法得到的次同步振荡频率为4.481 1、4.481 1、4.480 6、4.478 4、4.462 6 Hz。综上，4种算法结果基本上完全吻合，证明在振荡发散且振荡参数快速变化的场景下，频谱拟合法能进行百毫秒量级的动态参数辨识。且相较于改进插值DFT法得到的动态过程平均化结果，频谱拟合法可准确观测到某一时刻振荡的模态。尤其是当振荡处于严重发散阶段时，频谱拟合法仍能继续追踪振荡的发展和演变过程，如图 9所示的数据窗结尾时段。
![]()
图 9 场景一下不同算法次同步振荡频率和幅值辨识结果Figure 9. Frequency and amplitude results of sub-synchronous oscillation components with different algorithms in Scene Ⅰ2）场景二下超同步振荡结果：频谱拟合算法的次同步振荡分量和超同步振荡分量的参数辨识结果具有相似的特征，因此本节仅展示完整数据情况下超同步振荡的辨识结果，如图 10所示。超同步频率的辨识结果表明，本文提出的频谱拟合法、轨迹拟合法以及Prony算法所得的超同步频率在72.5~75 Hz之间振荡变化，它们的辨识结果轨迹一致；虽然改进插值DFT法无法监测振荡的动态变化过程，但其3组结果的变化轨迹也和其他3种算法结果轨迹相吻合。证明在振荡不发散且伴有超同步振荡的场景下，频谱拟合法仍能进行百毫秒量级的动态参数辨识。
![]()
图 10 场景二下不同算法超同步振荡频率和幅值辨识结果Figure 10. Frequency and amplitude results of sup-synchronous oscillation components with different algorithms in Scene Ⅱ4.2 数据丢失情况下
考虑到两个场景在数据丢失条件下的参数辨识结果具有相似的特征，为了验证在基于实际数据丢失情况下，频谱拟合算法噪声条件下的鲁棒性，本节仅展示场景一在数据丢失时的辨识结果。当发生数据丢失时，丢包率越大，辨识结果的准确度越小。在40 dB SNR时，丢包率为20%的情况下，辨识结果如图 11(a)所示。当发生数据丢失时，次同步振荡辨识结果仍可观察到振荡发散的动态发展过程。且振荡频率和幅值辨识结果相对误差均小于2×10−2，如图 11(b)。故即使存在40 dB噪声干扰，丢包率不大于20%时，仍可利用频谱拟合法准确追踪振荡的动态发展过程。
![]()
图 11 场景一在丢包率20%和完整数据情况下次同步分量频率和幅值辨识结果Figure 11. Identification results of frequency and amplitude of next synchronization component with packet loss rate of 20% and complete data in Scene Ⅰ5. 结论
本文提出一种基于同步相量频谱拟合的电力系统次/超同步振荡动态参数辨识方法，旨在利用频谱方法实现百毫秒量级的次/超同步振荡动态监测。首先，该算法通过构建同步相量的频谱拟合方程描述频谱的叠加特性，以更精确的数学解析式全频段频谱拟合代替传统频谱方法在频谱峰值处对频谱泄露的近似或忽略处理，实现理想条件下的零误差辨识。其次，频谱拟合法利用同步相量数据丢失条件下拟合方程仍成立的特性，实现了数据丢失条件下的参数辨识，摆脱了频谱方法对同步相量数据时域连续性的依赖。此外，本文利用频谱峰值幅值作为误差权重，利用频谱误差与振荡模态之间的关联关系，提高了参数辨识精度。综上所述，本文所提的频谱拟合法实现仅利用200 ms同步相量数据窗的高实时性次/超同步振荡参数辨识。
基于模拟PMU数据和仿真PMU数据的算例验证分析结果表明，相比于传统频谱方法，本文方法可以仅用1/10长度的数据窗实现对次/超同步振荡频率、幅值和相位的参数辨识，证明了该算法的可行性和有效性。
-
![]()
图 2 频谱叠加特性及频谱拟合的构成
Figure 2. Spectrum superposition characteristic and composition of the fitting spectrum
![]()
图 3 丢失数据点时频谱拟合算法原理示意图
Figure 3. Schematic diagram of the spectrum fitting algorithm when data loss occurs
![]()
图 4 x0=100，xsub=30，xsup=20，ϕ0=2.5，ϕsub=3，ϕsup=−3和fsub在[5, 45]Hz范围内变化时的各分量的频率、幅值和相位估计误差
Figure 4. Estimation errors with x0=100, xsub=30, xsup=20, ϕ0=2.5, ϕsub=3, ϕsup=−3 and fsub varying in the range [5, 45]Hz
![]()
图 5 x0=100，xsub=30，xsup=20，ϕ0=2.5，ϕsub=3，ϕsup=−3和fsub在[5, 45]Hz范围内变化时的各分量的频率、幅值和相位估计误差(加40 dB的噪声)
Figure 5. Estimation errors with x0=100, xsub=30, xsup=20, ϕ0=2.5, ϕsub=3, ϕsup=−3 and fsub varying in the range [5, 45]Hz (with 40 dB noise)
![]()
图 6 40 dB噪声条件下，丢包率为20%时各分量的频率、幅值和相位估计误差
Figure 6. Estimation errors when data loss rate is 20% under 40 dB noise condition
![]()
图 9 场景一下不同算法次同步振荡频率和幅值辨识结果
Figure 9. Frequency and amplitude results of sub-synchronous oscillation components with different algorithms in Scene Ⅰ
![]()
图 10 场景二下不同算法超同步振荡频率和幅值辨识结果
Figure 10. Frequency and amplitude results of sup-synchronous oscillation components with different algorithms in Scene Ⅱ
![]()
图 11 场景一在丢包率20%和完整数据情况下次同步分量频率和幅值辨识结果
Figure 11. Identification results of frequency and amplitude of next synchronization component with packet loss rate of 20% and complete data in Scene Ⅰ
表 1 同步相量频谱拟合算法与其他算法的性能比较
Table 1 Performance comparison between synchronous phasor spectrum fitting algorithm and other algorithms
各分量辨识结果相对误差① 噪声条件SNR/dB 数据丢失率/% 频谱拟合 文献[24] (轨迹拟合) 文献[16] (改进的插值DFT) 0.2 s 0.2 s③ 0.2 s 1 s 2 s |f0−ˆf0|f0 ∞② 0 0 10−15 10−3 10−4 10−4 40 0 2×10−3 10−15 10−3 10−4 10−4 20 2×10−3 — — — — |x0−ˆx0|x0 ∞ 0 10−14 10−14 10−2 10−3 10−2 40 0 6×10−3 10−2 10−2 10−2 10−2 20 10−2 — — — — |ϕ0−ˆϕ0|ϕ0 ∞ 0 10−14 10−14 10−2 10−2 10−2 40 0 3×10−2 10−1 10−2 10−2 10−2 20 3×10−2 — — — — |fsub−ˆfsub|fsub ∞ 0 10−14 10−14 10−1 10−3 10−4 40 0 5×10−3 10−2 0.5 > 1 10−3 20 10−2 — — — —/ |xsub−ˆxsub|xsub ∞ 0 10−12 10−12 > 1 10−1 10−2 40 0 10−2 5×10−2 > 1 1 10−1 20 2×10−2 — — — —/ |ϕsub−ˆϕsub|ϕsub ∞ 0 10−12 10−12 > 1 0.15 10−1 40 0 6×10−3 10−1 > 1 > 1 10−1 20 2×10−2 — — — — ①本表所示结果为不同fsub下的辨识结果相对误差最大量级。不同算法的测试条件均为fsub在[5, 45]Hz范围内变化，fsup、xsup、ϕsup的结果误差与fsub、xsub、ϕsub的结果误差相近而未列出。
②SNR为信噪比(signal-noise ratio)，SNR=∞表示没有噪声。
③轨迹拟合在噪声条件下使用0.1 s数据窗的相位辨识结果误差比使用0.2 s数据窗的误差大，相位辨识结果误差大于10%，不满足精度要求。
下载: 导出CSV
-
[1] 李浩志, 谢小荣, 刘芮彤, 等. 新能源经柔直送出系统的次同步振荡分析与抑制[J]. 中国电机工程学报, 2024, 44(8): 2945-2955. DOI: 10.13334/j.0258-8013.pcsee.223006 LI Haozhi, XIE Xiaorong, LIU Ruitong, et al. Analysis and mitigation of the subsynchronous oscillation in renewable energy system connected to the MMC-HVDC [J]. Proceedings of the CSEE, 2024, 44(8): 2945-2955(in Chinese). DOI: 10.13334/j.0258-8013.pcsee.223006
[2] 谢小荣, 刘华坤, 贺静波, 等. 电力系统新型振荡问题浅析[J]. 中国电机工程学报, 2018, 38(10): 2821-2828. DOI: 10.13334/j.0258-8013.pcsee.172662 XIE Xiaorong, LIU Huakun, HE Jingbo, et al. On new oscillation issues of power systems[J]. Proceedings of the CSEE, 2018, 38(10): 2821-2828(in Chinese). DOI: 10.13334/j.0258-8013.pcsee.172662
[3] 高本锋, 刘洋洋, 王潇. 基于路径分析法的直驱风电场经柔直并网系统次同步振荡机理分析[J]. 中国电机工程学报, 2023, 43(17): 6646-6660. DOI: 10.13334/j.0258-8013.pcsee.221589 GAO Benfeng, LIU Yangyang, WANGXiao. Subsynchronous oscillation mechanism analysis of direct-drive wind farm with VSC-HVDC system based on path analysis method [J]. Proceedings of the CSEE, 2023, 43(17): 6646-6660(in Chinese). DOI: 10.13334/j.0258-8013.pcsee.221589
[4] LI Yin, FAN Lingling, MIAO Zhixin. Replicating real-world wind farm SSR events[J]. IEEE Transactions on Power Delivery, 2020, 35(1): 339-348. DOI: 10.1109/TPWRD.2019.2931838
[5] 陈曦, 吴熙, 顾文, 等. 基于次同步模态能量的双馈风机并网系统次同步振荡能量分析与溯源方法研究[J]. 中国电机工程学报, 2023, 43(15): 5863-5877. DOI: 10.13334/j.0258-8013.pcsee.220997 CHEN Xi, WU Xi, GU Wen, et al. Research on energy analysis and SSO source location method for power grids with DFIGs based on the sub-synchronous modal energy[J]. Proceedings of the CSEE, 2023, 43(15): 5863-5877(in Chinese). DOI: 10.13334/j.0258-8013.pcsee.220997
[6] WANG Liang, XIE Xiaorong, JIANG Qirong, et al. Investigation of SSR in practical DFIG-based wind farms connected to a series-compensated power system [J]. IEEE Transactions on Power Systems, 2015, 30(5): 2772-2779. DOI: 10.1109/TPWRS.2014.2365197
[7] 董青迅, 李兴源, 穆子龙, 等. 基于Prony算法的次同步谐振检测方法[J]. 电力系统及其自动化学报, 2012, 24(2): 1-4, 34. DOI: 10.3969/j.issn.1003-8930.2012.02.001 DONG QINGXUN, LI XINGYUAN, MU Zilong, et al. New method for detection of sub-synchronous resonance based on Prony algorithm[J]. Proceedings of the CSU-EPSA, 2012, 24(2): 1-4, 34(in Chinese). DOI: 10.3969/j.issn.1003-8930.2012.02.001
[8] PATERNINA M R A, TRIPATHY R K, ZAMORA- MENDEZ A, et al. Identification of electromechanical oscillatory modes based on variational mode decomposition[J]. Electric Power Systems Research, 2019, 167: 71-85. DOI: 10.1016/j.epsr.2018.10.014
[9] GAO Bo, TORQUATO R, XU W, et al. Waveform-based method for fast and accurate identification of subsynchronous resonance events[J]. IEEE Transactions on Power Systems, 2019, 34(5): 3626-3636. DOI: 10.1109/TPWRS.2019.2904914
[10] 王茂海, 高洵, 王蓓, 等. 基于广域测量系统的次同步振荡在线监测预警方法[J]. 电力系统自动化, 2011, 35(6): 98-102. WANG Maohai, GAO Xun, WANG Bei, et al. Online early-warning of sub-synchronous oscillations based on wide area measurement system[J]. Automation of Electric Power Systems, 2011, 35(6): 98-102(in Chinese).
[11] RAUHALA T, GOLE A, JÄRVENTAUSTA P. Detection of subsynchronous torsional oscillation frequencies using phasor measurement[J]. IEEE Transactions on Power Delivery, 2016, 31(1): 11-19. DOI: 10.1109/TPWRD.2015.2436814
[12] ZHANG Fang, CHENG Lin, GAO Wendong, et al. Synchrophasors-based identification for subsynchronous oscillations in power systems[J]. IEEE Transactions on Smart Grid, 2019, 10(2): 2224-2233. DOI: 10.1109/TSG.2018.2792005
[13] WANG Yang, JIANG Xiaolong, Xie Xiaorong, et al. Identifying sources of subsynchronous resonance using wide-area phasor measurements[J]. IEEE Transactions on Power Delivery, 2021, 36(5): 3242-3254. DOI: 10.1109/TPWRD.2020.3037289
[14] LIU Hao, BI Tianshu, CHANG Xiqiang, et al. Impacts of subsynchronous and supersynchronous frequency components on synchrophasor measurements [J]. Journal of Modern Power Systems and Clean Energy, 2016, 4(3): 362-369. DOI: 10.1007/s40565-016-0225-4
[15] YANG Xiaomei, ZHANG Jianing, XIE Xiaorong, et al. Interpolated DFT-based identification of sub- synchronous oscillation parameters using synchrophasor data[J]. IEEE Transactions on Smart Grid, 2020, 11(3): 2662-2675. DOI: 10.1109/TSG.2019.2959811
[16] ZHANG Fang, LI Jiaxin, LIU Jun, et al. An improved interpolated DFT-based parameter identification for sub-/super-synchronous oscillations with synchrophasors [J]. IEEE Transactions on Power Systems, 2023, 38(2): 1714-1727. DOI: 10.1109/TPWRS.2022.3175735
[17] 马钺, 蔡东升, 黄琦. 基于Rife-Vincent窗和同步相量测量数据的风电次同步振荡参数辨识[J]. 中国电机工程学报, 2021, 41(3): 790-802. DOI: 10.13334/j.0258-8013.pcsee.201400 MA Yue, CAI Dongsheng, HUANG Qi. Parameter identification of wind power sub-synchronous oscillation based on Rife-Vincent window and synchrophasor data [J]. Proceedings of the CSEE, 2021, 41(3): 790-802(in Chinese). DOI: 10.13334/j.0258-8013.pcsee.201400
[18] WANG Yuhong, WANG Hongyu, SONG Yuyan, et al. Parameter identification of Sub-synchronous/Super- synchronous oscillations based on synchrophasor rotation and spectral shift[J]. International Journal of Electrical Power & Energy Systems, 2023, 149: 109044.
[19] SHAIR J, XIE Xiaorong, YANG Jianjun, et al. Adaptive damping control of subsynchronous oscillation in DFIG-based wind farms connected to series-compensated network[J]. IEEE Transactions on Power Delivery, 2022, 37(2): 1036-1049. DOI: 10.1109/TPWRD.2021.3076053
[20] 王杨, 晁苗苗, 谢小荣, 等. 基于同步相量数据的次同步振荡参数辨识与实测验证[J]. 中国电机工程学报, 2022, 42(3): 899-908. DOI: 10.13334/j.0258-8013.pcsee.210792 WANG Yang, CHAO Miaomiao, XIE Xiaorong, et al. Identification of subsynchronous oscillation parameters and field tests based on PMU data [J]. Proceedings of the CSEE, 2022, 42(3): 899-908(in Chinese). DOI: 10.13334/j.0258-8013.pcsee.210792
[21] 王渝红, 王宏宇, 于光远, 等. 基于同步相量数据的次同步振荡检测与模态参数辨识方法[J]. 高电压技术, 2023, 49(6): 2557-2568. WANG Yuhong, WANG Hongguang, YU Guangyuan, et al. Sub-synchronous oscillation detection and modal parameter identification method based on synchrophasor [J]. High Voltage Engineering, 2023, 49(6): 2557-2568(in Chinese).
[22] QI Yuan, LI Jiaxian, LIU Hao, et al. Subsynchronous oscillation monitoring and alarm method based on phasor measurements[J]. Global Energy Interconnection, 2020, 3(5): 464-474.
[23] LIU Hao, QI Yuan, ZHAO Junbo, et al. Data-Driven subsynchronous oscillation identification using field synchrophasor measurements[J]. IEEE Transactions on Power Delivery, 2022, 37(1): 165-175.
[24] 张放, 刘军, 李佳欣, 等. 基于同步相量轨迹拟合的电力系统次同步/超同步振荡的实时参数辨识[J]. 中国电机工程学报, 2023, 43(4): 1413-1425. DOI: 10.13334/j.0258-8013.pcsee.212806 ZHANG Fang, LIU Jun, LI Jiaxin, et al. Real-time parameter identification with synchrophasor trajectory fitting technique for subsynchronous/supersynchronous oscillations in power systems[J]. Proceedings of the CSEE, 2023, 43(4): 1413-1425(in Chinese). DOI: 10.13334/j.0258-8013.pcsee.212806
[25] 顾文, 吕万, 杨宏宇, 等. 电力系统次同步振荡的检测技术综述[J]. 电测与仪表, 2022, 59(7): 9-18. GU Wen, LV Wan, YANG Hongyu, et al. Review of detection technology of sub-synchronous oscillation in power systems[J]. Electrical Measurement & Instrumentation, 2022, 59(7): 9-18(in Chinese).
[26] 刘灏, 李珏, 毕天姝, 等. 基于PMU相量的次/超同步间谐波识别方法[J]. 电网技术, 2017, 41(10): 3237-3243. LIU Hao, LI Jue, BI Tianshu, et al. Subsynchronous and supersynchronous inter-harmonic identification method based on phasor measurements[J]. Power System Technology, 2017, 41(10): 3237-3243(in Chinese).
[27] YANG Xiaomei, YANG Lin, XIAO Xianyong, et al. A novel detection method for supersynchronous resonance from synchrophasor data[J]. IEEE Transactions on Power Systems, 2023, 38(4): 3694-3706.
计量
- 文章访问数: 0
- HTML全文浏览量: 0
- PDF下载量: 0
下载:
