基于改进Prony 算法的新能源配网谐波和振荡复合扰动参数提取

李 磊，都成刚，王维兵，齐少春，杨学林，杨瑞卿

（1.内蒙古电力（集团）有限责任公司，呼和浩特 010020；2.南京南瑞继保工程技术有限公司，南京 211122；3.内蒙古电力（集团）有限责任公司阿拉善供电分公司，阿拉善盟 750300；4.华北电力大学电力工程系，保定 071003）

目前，以新能源发电并网逆变器和变流器为代表的电力电子设备的大规模应用[1]，使得配电网中的背景谐波问题愈发严重。同时，新型电力电子设备及其控制系统与传统电磁设备存在强耦合关系，会导致电力系统的动态行为更加复杂[2]。在此背景下，电网更易出现宽频振荡等安全问题[3]。

配电网中成分复杂的谐波会增加网络损耗，并使得发生网络谐波谐振的概率增加，从而引起谐波电流放大[4]，严重时甚至影响设备的正常工作。而当配电网发生宽频振荡时，如果没有得到合理的处置，电压电流的大幅波动很容易导致新能源发电设备大量脱网，严重影响电网的安全稳定运行，对与主网电气距离间隔较大的配网的危害尤其明显[5]。各类电能质量问题对电力用户的安全和电网的安全稳定运行构成了新的威胁。因此，准确提取电网中背景谐波和宽频振荡参数，以此指导电能质量问题的治理也显得迫在眉睫。

电能质量参数提取方法主要有离散傅里叶变换DFT（discrete Fourier transform）、小波变换、希尔伯特-黄变换HHT（Hilbert-Huang transform）、S 变换ST（Stockwell transform）和Prony 算法等。离散傅里叶变换可以获得扰动信号的视频信息，但难以准确提取暂态过程信号参数；小波变换适合分析暂态信号，但对背景谐波这类小扰动信号参数提取能力较弱[6]；希尔伯特-黄变换存在端点效应、模态混叠等问题[7]，不适宜分析复合扰动；S 变换有良好的时频域分析精度，但难以直观区分扰动信号类型[8]，需额外添加辅助算法对扰动类型进行分类[9]；Prony算法能辨识同一时刻内多种扰动的参数信息，但算法对采样信号中的噪声特别敏感，抗噪声能力较差[10]。同时，Prony 算法信号辨识结果的精度很大程度上取决于模型有效秩阶数的选取是否合理，合理的阶数能大大降低计算量，而不合理的阶数会导致部分参数辨识结果遗漏。

针对Prony 算法的抗噪声能力差，使用时需搭配降噪算法，常用的有经验模态分解EMD（empirical mode decomposition）、模糊滤波法和奇异值分解SVD（singular value decomposition）法。EMD 和模糊分解法存在模态混叠和端点效应，分析复合扰动时降噪精度较差，且容易滤除有效信息[11]；而SVD 法在分析复合扰动时不易遗漏有效成分，降噪效果也较为理想[12]。对于有效秩阶数的计算，主要有经验小波变换EWT（empirical wavelet transform）法和二阶导数理论SDM（second derivative method）法。其中EWT 法需要对信号进行预处理，会降低参数辨识速度，多适用于信号成分存在耦合的情况[13]；而SDM法通过计算有效秩阶数增长率为0的点，能够快速地确定扰动对应的有效秩阶数。

为此，本文提出了一种基于SVD 和SDM-Prony算法的谐波和振荡复合扰动参数提取方法。先应用SVD 法对采样信号进行降噪，再引入SDM 法，针对复合扰动信号自适应地设定Prony算法的有效秩阶数，最大程度地提高参数辨识精度。仿真分析表明，在进行复合扰动特征分析时，引入的2 种优化方法减少了Prony算法参数提取的相对误差。

1 算法基本原理

Prony算法的实质是一种基于多项式线性拟合的参数辨识方法。设采样信号数据为x(n)，则拟合x(0)，x(1)，…，x(N-1)，可以记为(n)，表示为

式中：p为模型的阶数；N为信号采用点的个数；bi和zi为复数，其中A为复数模值，θi为复数相角，αi为复数衰减因子，fi为复数频率。并且，为确保拟合信号尽可能接近实际采样信号，这里需要采用误差平方和ε最小作为优化目标，对应其目标函数，表示为

将采样信号函数x(n)处理为

式中：r(i,j)为矩阵R中的元素；pe为矩阵R的阶数。根据式（5），构造采样样本矩阵R，矩阵中元素和结构表示为

将式（6）中矩阵R进行奇异值分解或最小二乘估计，可以求取有效秩pe，其大小取决于信号扰动量的个数。当信号没有噪声时，pe的物理含义为将样本矩阵R奇异值分解后奇异值矩阵中非零元素的个数。

保留有效秩阶数pe所对应的矩阵R部分元素，然后联立求解，有

式中：ap为特征向量；up为扰动变量。利用式（7）中的ai求解多项式（8），得到其对应的根zi。然后将zi代入式（1），并进行变换，可得到矩阵方程为

式（9）最左侧矩阵是一个N×p维的范德蒙矩阵，矩阵满秩且可逆，记为矩阵V，因此可进行变换为

接着，就可以计算出混合扰动信号的特征参数，主要包括：宽频振荡的幅值Ai、衰减因子αi和相位θi、谐波信号的频率fi，分别表示为

式中，i=0，1，…，N-1，表示扰动量的数目。在参数计算时，可以将谐波信号看作是衰减因子αi=0 的振荡参数。

2 基于SVD 的SDM-Prony 算法原理

本文提出的基于SVD 的SDM-Prony 算法提取谐波和振荡复合扰动参数的整体流程如图1所示。

图1 参数提取方法流程Fig.1 Flow chart of parameter extraction method

首先，对采集到的含噪声和扰动的离散电信号进行SVD 降噪，滤除噪声信号后，重新生成采样信号序列。接着，应用SDM 法，计算样本函数二阶导函数为0 时所对应的序列号，从而确定样本矩阵R的有效秩p，并依据p设定Prony 算法的模型精度和特征方程的阶数。最后应用Prony算法提取复合扰动的参数信息。

2.1 SVD 降噪算法原理

由于采集到的并网点处电信号包含着大量的高斯噪声，同时Prony 算法参数的提取结果受噪声影响较大。因此，使用Prony 算法时必须搭配使用降噪算法，先对采集到的电信号进行降噪处理，然后使用Prony算法提取扰动参数特征。

基于SVD 的降噪算法原理是将原本的采样信号向量空间分为噪声信号空间和有效信号空间，且噪声信号和有效信号空间彼此正交。针对采样点处采集到的离散信号：x=[x1,x2,…,xN]，可构造Hankel矩阵X为

其中：如果采样点个数N为奇数，则n=（N+1）/2，m=（N+1）/2；如果N为偶数，则n=N/2，m=N/2+1。

接着对矩阵X进行奇异值分解为

式中：q为非奇异值元素阶数；U为m×m维的正交矩阵；Σ为n×n维矩阵，其中Δ是对角矩阵，对角线上元素是矩阵X的奇异值；V为n×n维可逆的正交矩阵。

依据奇异值分解理论，若存在k维有效信号空间，则对角矩阵Δ中一定存在与之对应的k个奇异值，余下的q-k个奇异值就对应噪声信号空间的维数。k就是有效信号空间所对应的阶数。

SVD 算法即是利用矩阵Σ奇异值序列的分布特点，甄别有效信号空间所对应的阶数k。对角矩阵Δ对角线上的元素即为矩阵X对应的奇异值序列，将该序列按照奇异值元素绝对值大小进行排序，可以写为：σ11,σ22,…,σkk,σ(k+1)(k+1), …,σnn。其中奇异值序列中σ11,σ22,…,σkk元素数值较大，且数值变化幅度较小。而第k+1个奇异值σ(k+1)(k+1)数值较第k个奇异值σkk会出现大幅下降，之后，序列σ(k+1)(k+1),σ(k+2)(k+2),…,σnn数值变化会再次趋于平稳，没有明显的波动，其原因是，采样信号中噪声分量相较于有效信号分量要小得多。

对含有20 dB高斯白噪声和复合扰动的电信号进行采样，采样频率设置为2 000 Hz。由采样信号生成的Hankel 矩阵奇异值序列如图2 所示。若将奇异值序列视为一个按照大小排布的数列，那么数列数值变化的突变点所在位置就是有效阶数k。由图2可以直观地看出，其有效阶数为4。

图2 含噪信号Hankel 矩阵奇异值序列Fig.2 Sequence of singular values of Hankel matrix with noisy signals

在实际计算过程中，通过设定阈值函数v（k），可求取矩阵的有效阶数，则有

选取阈值函数v（k）大于阈值时所对应的最小k值即为有效阶数。实际SVD算法使用时，通常阈值选为0.90～0.95时效果最好。

如果k值选取过小，虽然滤波性能良好，但部分有效信号分量也被滤除，会导致重构的采样型号特征分量缺失和提取结果不准确。而如果k值选取过大，虽然能很好保留有效信号分量，但滤波性能就会下降。

之后，对矩阵Σ中的分块矩阵Δ进行处理，将q-k阶对角元素置0，修改后的矩阵记为Σk。然后对其进行变换，有

2.2 SDM 法基本原理

式（6）所示的矩阵R，其有效秩p的选取结果会影响Prony算法提取结果的精度。由于待检测信号复合扰动数目是未知的，传统Prony 算法往往事先设定模型阶数。如果有效秩p选取过小，容易在参数提取时遗漏部分参数信息，使得参数提取结果不完备。而如果有效秩p选取过大，会造成辨识结果参数信息冗余[14]。

基于SDM 的矩阵有效秩提取方法是通过将式（7）中矩阵R的奇异值按照大小顺序进行排序，进而得到序列，即

由于实际电信号中噪声的存在，σp+1及其后的奇异值并不为0，是一组数值趋于0 的序列。通过拟合离散的奇异值数据点生成函数，再求取函数二阶导数，如满足

即存在二阶导数为0的临界点。其中变量p对应正整数序列N+，表明在该临界点前后矩阵奇异值发生了突变。该点就可以视为有效秩的分界点，此临界点对应的序号即表示有效秩阶数p。

求取的有效秩阶数p可代入式（7）～式（9），极大减少了计算量，同时也使得式（11）提取出的参数数目更接近于真实扰动量的数目，极大减少了无效数据的数量。

3 算例仿真分析

3.1 算例仿真

电网中由于非线性负载的存在，往往含有一定比例的背景谐波，同时，伴随逆变器工作状态的改变，会产生宽频振荡这一过渡过程。仿真算例所研究的信号模型是同时包含了谐波和振荡的复合扰动，其数学表达式为

式中：Acos(100πt) 为电网中基波信号；Aicos(2πfit+φi)e-αit为复合扰动信号分量。背景谐波可以看作是衰减因子为0 的宽频振荡，信号分析宽频振荡和背景谐波可以通过衰减因子进行区别。

复合扰动参数的设定如表1 所示，同时在信号中加入20 dB的高斯白噪声。

表1 信号复合扰动参数Tab.1 Signal composite disturbance parameters

工矿企业中大量使用的整流器是5次和7次背景谐波的主要来源，而非线性的照明设备如节能灯，是3次背景谐波的主要来源。

信号采用频率设定为2 000 Hz，仿真时间设定为1 s。采样序列共包含2 000个元素，用采样序列拟合原始信号，其波形图3所示。

图3 含噪声原始信号波形Fig.3 Waveform of raw signal with noise

应用SVD 降噪算法后对原始信号进行处理，重新生成采样序列，其拟合的波形如图4所示。

图4 SVD 降噪后信号波形Fig.4 Waveform of signal after SVD noise reduction

对比图3 和图4 可以看出，经过SVD 法重构的信号序列，其噪声分量得到了有效的滤除。对降噪重构后的信号波形应用Prony算法进行复合扰动参数提取，可以提高提取精度。其中，暂态振荡参数提取结果和相对误差如表2所示。

表2 宽频振荡参数提取结果Tab.2 Extraction results of broadband oscillation parameters

从表2 可以看出，基于SDM 的SVD-Prony 算法能够有效提取信号中的振荡参数，对于宽频振荡参数，由于其振荡幅度一般相对较大，扰动特征明显，因此提取也相对容易，结果也更为准确。谐波参数提取结果如表3所示。

表3 谐波参数提取结果Tab.3 Harmonic parameter extraction results

从表3 可以看出，即使是扰动特征不很明显的背景谐波，其参数提取结果的相对误差也在可接受误差范围内。

误差主要来源于SVD噪声滤除环节，由于SVD降噪算法中需要设定有效信号空间阶数k。阶数k如果设置过大，虽然滤波效果会更好，但会使得信号部分有效成分丢失，反而影响后续Prony算法提取结果的准确性，严重时甚至会遗漏部分参数信息。

3.2 算法比较分析

仍采用式（12）所示的复合信号模型，并分别采用原Prony 算法、SVD-Prony 算法和基于SDM 法的SVD-Prony 算法进行扰动信号的参数提取。3 种方法对扰动信号频率fi、参数幅值Ai和衰减因子αi提取结果对比如表4所示。

表4 Prony 算法及其改进法参数提取结果对比Tab.4 Comparisonofparameterextractionresult between Prony algorithm and its improved versions

可以看出，如果不加SVD 降噪算法，提取结果误差很大，且频率越高的扰动受噪声信号影响越大，提取结果的精度越低。宽频振荡由于其振荡范围明显，特征较为突出，反而受噪声影响较低。基于SVD-Prony 算法由于对噪声进行了滤除，其精度有了很大程度的提升，但相对误差仍然较大。基于SDM 的SVD-Prony 算法更进一步，从有效秩的角度入手，在滤波操作过后，运用二阶导数理论，进一步计算扰动模型的阶数，使得相对误差能够做到3种方法中最小。但与之对应，参数提取时间相较于原算法更长。

同时，希尔伯特-黄变换和S变换也能提取复合扰动参数信息。而离散傅里叶变换仅能提取谐波参数，无法提取暂态振荡参数。同时对比其他算法对同一信号扰动参数的提取结果，如表5所示。

表5 改进Prony 算法与其他谐波检测算法提取结果对比Tab.5 Comparison of extraction result between improved Prony algorithm and other harmonic detection algorithms

从表5 可以看出，对于暂态振荡扰动，S 变换、希尔伯特-黄变换和改进Prony 算法提取结果的相对误差接近。但对于背景谐波扰动，改进Prony算法的相对误差小得多，原因在于改进Prony 算法添加了噪声滤除环节和模型参数选取环节。此外，3 种不同类型算法的参数提取时间对比如表6所示。

表6 算法参数提取时间对比Tab.6 Comparison of parameter extraction time among different algorithms

从表6 可以看出，S 变换和希尔伯特-黄变换参数提取所需时间更长。经仿真算例验证，改进Prony算法能准确提取谐波和振荡复合扰动参数信息，且较其他类型算法所需时间更短。

3.3 现场数据分析

为了检验本文提出的改进Prony 算法的有效性，现选取现场试验数据进行分析，数据来源于厦门某变电站内。当电力电子设备，如逆变器和变流器在运行过程中进行工作状态切换时，会产生电压、电流振荡和谐波问题，如图5所示[15]。

图5 现场电压、振荡和含背景谐波数据Fig.5 Data of field voltage,oscillations and background harmonics

对变电内采集到的现场数据，运用基于SVD和SDM-Prony算法进行参数提取，其结果如表7所示。

表7 现场数据参数提取结果Tab.7 Field data parameter extraction results

可以看出，背景谐波多为特征次谐波，复合扰动包含3、5、7、11 和13 次背景谐波和一个频率为232.64 Hz 的振荡。将该参数提取结果与现场电能质量分析仪器的提取结果进行对比，结果如表8所示。

表8 现场数据参数提取结果对比Tab.8 Comparison of filed data parameter extraction result

电能质量分析仪基于微机系统，针对谐波扰动采用DFT算法提取扰动参数信息，针对振荡扰动采用改进S变换法提取扰动信息，由于经过多次迭代计算，其提取结果相对准确，但耗时长得多。通过对比现场数据分析结果，本文提出的改进Prony 算法的扰动参数提取结果的相对误差也处于一个较小的区间内。

4 结 论

本文提出的基于SDM 和SVD-Prony 算法具有以下的特点和优势：

（1）针对原有的Prony 算法在进行扰动参数提取时易受噪声影响且提取结果很大程度受模型阶数影响，本文提出的改进方法首先采用SVD 降噪法，对电网中含噪声的信号进行滤波并重构采样信号。接着，运用SDM 法，对重构后的信号序列进行分析，计算拟合函数后二阶导函数为0 的序列点，进而确定Prony 算法模型的阶数。最后应用Prony算法提取扰动信号的频率、幅值、衰减因子等参数信息。

（2）仿真实验和现场数据分析证明了方法的有效性。同时对比分析了3 种算法对同一信号波形的参数提取结果，对比发现基于SVD 的SDM-Prony算法具有最优的提取精度。