★ FEATURED ARTICLE 做电力系统同步相量计算我第一反应就是FFT。确实FFT够快一段采样序列扔进去复数运算做完幅值和相位就摆在眼前。可真正把算法放到动态工况下你会发现频谱泄漏、栅栏效应、谐波干扰轮番上阵电网频率只要偏移0.2Hz传统DFT类算法就开始“浑身难受”。于是窗函数法、希尔伯特-黄变换、小波变换都被摆上台面和FFT组合成一套互相验证的电力系统同步相量计算研究。这篇内容我基于Matlab完整实现了一遍把四种方法的原理、代码、实测误差和踩坑过程梳理出来适合刚接触同步相量测量方向的研究生也适合正在做PMU算法选型的工程开发。1. 背景与问题定位为什么同步相量计算需要这么多工具1.1 同步相量是什么半个周期到底要算什么同步相量不是简单地把电压电流波形数字化它要求的是“带时标的相量”。电力系统广域测量、新能源并网控制、储能系统响应评估都要依赖同步相量来感知全局状态。相量本身包含三个核心量幅值、相角、频率以及可选的频率变化率。所谓“同步”是指不同测点之间要有统一时标约束算法在处理采样序列时必须假定某一时刻作为参考点然后给出此刻的相量估计。实际工程中的同步相量计算标准要求相当苛刻。幅值误差、相位误差、总向量误差TVE、频率误差和延时响应都是考核指标。我记得IEEE标准里对稳态、动态、谐波等场景都画了误差界限这里不展开标准条文但至少要知道一点同步相量算法不是算得越准越好还要兼顾响应速度。你不可能为了追求稳态精度把数据窗拉到几十个周波那动态事件早就过去了。1.2 传统DFT/FFT的“天花板”在哪里DFT/FFT的数学本质是把有限长信号投影到一系列等间隔频率的复指数基上它隐含了一个前提信号在这段观察窗内严格周期并且窗长度是基波周期的整数倍。理想电网是50Hz你采样整数个周波FFT峰值谱线正好落在50Hz那一格幅值和相位误差接近机器精度。这一度让人觉得FFT已经够了。可电网不是教科书频率会漂。发电机出力波动、负荷变化、新能源出力随机性都会让实际频率偏离50Hz。一旦频率不是精确的50Hz你的数据窗长度就不再是基波周期的整数倍频谱能量会从主瓣“漏”到旁瓣这就是频谱泄漏。同时FFT能分辨的频率只落在离散谱线上真实频率往往位于两条谱线之间峰值位置被“栅栏”挡住栅栏效应又带来幅度低估和相位偏移。谐波和噪声再掺进来基波谱线附近就会叠加上各种分量误差进一步恶化。FFT本身不是一个错误选项它只是把你带到一个“整周期同步采样”的理想假设里。窗函数、希尔伯特-黄变换、小波变换都是在这个假设不成立时想方设法把基波分量更干净地“抠”出来的工具。理解这一点算法的选择逻辑就清晰了。2. 四种方法的原理拆解与选型逻辑2.1 FFT用最短的时间拿到频谱全貌FFT是离散傅里叶变换的快速实现复杂度从O(N²)降到O(NlogN)。对N点采样频率分辨率是fs/N第k条谱线对应的频率是k*fs/N。在同步相量计算里FFT并不需要对全频段都有兴趣真正关心的就是基波附近那几根谱线。我要强调一个细节FFT输出的正频率谱线幅值是信号幅值的N/2倍对单频余弦信号而言相位等于信号初始相位。所以从FFT结果里提取相量的公式是X fft(x); N length(x); [~, k_peak] max(abs(X(1:floor(N/2)))); % 峰值谱线索引 A 2 * abs(X(k_peak)) / N; % 幅值估计 phi angle(X(k_peak)); % 相位估计但这段代码只能在整周期同步采样时可靠。我把仿真信号设成50Hz、数据窗10个周波、采样率10kHz时误差确实小到0.001%。一旦频率改成50.2Hz同样的代码幅值误差就能飙到4%以上相位误差也有好几度。这就是FFT的天花板它快、稳但对“非同步采样”没有免疫力。2.2 窗函数法给FFT戴上“降噪耳塞”既然截断会泄漏那就用窗函数来截断。窗函数法的核心思想是让数据窗两端平滑衰减抑制旁瓣泄漏。常见窗包括矩形窗、Hanning窗、Hamming窗、Blackman窗。矩形窗旁瓣电平只有约-13dB而Hanning窗能把第一旁瓣压到-31dB左右Blackman窗可以到-58dB附近。代价是什么主瓣变宽了。矩形窗主瓣宽度是2ΔfΔffs/NHanning窗变成4ΔfBlackman窗约6Δf。主瓣宽了频域上靠近基波的谐波和间谐波就不容易分开。所以窗函数本质上是“用频率分辨率的损失换取频谱泄漏的降低”。加窗之后还需要幅值恢复。因为窗函数把采样点加权了FFT峰值不再是NXm/2而是Xmsum(w)/2因此幅值恢复系数是2/sum(w)。我把Hanning窗代码写成w hann(N, periodic).; xw x .* w; Xw fft(xw); [~, k_peak] max(abs(Xw(1:floor(N/2)))); A_win 2 * abs(Xw(k_peak)) / sum(w); phi_win angle(Xw(k_peak));用这个方法在50.2Hz下测试幅值误差从4%降到了0.5%左右。如果再配合双谱线插值利用峰值左右两根谱线的比值估算真实频率偏移误差还能再降一个量级。窗函数法不是新东西但它依然是工程上性价比最高的稳定方案。2.3 希尔伯特-黄变换从数据本身来找模态希尔伯特-黄变换HHT由经验模态分解EMD和希尔伯特变换两部分组成。EMD和FFT、小波不一样它不预设基函数而是根据信号自身的时间尺度把信号逐级分解成本征模态函数IMF。每个IMF瞬时频率有意义可以直接做希尔伯特变换得到瞬时幅值和瞬时相位。我在Matlab里直接用自带的emd函数[imf, ~] emd(x, Display, 0); analytical hilbert(imf); % 对每个IMF做解析延拓 inst_amp abs(analytical); inst_phase unwrap(angle(analytical)); inst_freq fs / (2*pi) * diff(inst_phase);EMD的强项在于处理非平稳、非线性信号。电网发生次同步振荡、频率斜坡、幅值调制时EMD能把基波和振荡分量分离成不同IMF然后用希尔伯特变换得到瞬时包络和瞬时相位。这比固定带通滤波更自适应。但它也不完美。端点效应很讨厌信号两端会因为EMD的样条插值出现大幅摆动直接影响相量估计值。模态混叠也常有如果谐波和基波频率接近又幅度大EMD会“撕不开”这两个分量把它们混进同一个IMF。我在实验里发现稳态信号用HHT反而不占便宜误差比FFT还大主要就是端点振荡在作怪。2.4 小波变换固定在时间轴上的“伸缩放大镜”小波变换和FFT的视角完全不同。FFT在频率维上展开损失了时间定位连续小波变换CWT用一组可伸缩、可平移的小波基函数去匹配信号在一个时频平面上同时表达频率和时间。低频分量用宽窗、高频率分辨率高频分量用窄窗、高时间分辨率这一点很像“变焦放大镜”。在Matlab里调用连续小波变换很直接[wt, fr] cwt(x, amor, fs); % 使用复Morlet小波 f_target 50; [~, idx] min(abs(fr - f_target)); coef squeeze(wt(idx, :)); % 50Hz脊线上的小波系数复小波系数包含幅值和相位信息。提取出基波脊线后可以做带通重构也可以直接从小波系数换算出瞬时相量。小波变换对频率偏移耐受性好因为它在时频面上追踪的是真实频率变化而不是假设频率固定为50Hz。缺点也很明显小波基函数要选Morlet、Morse、Amor结果差别不小边界区域的小波系数会被锥形区域污染数据窗两端不可信计算量比FFT大一个数量级。这些都要在落地时权衡。2.5 四种方法对比速览方法基函数频率偏移耐受动态跟踪能力抗谐波噪声计算成本典型定位FFT复指数差差一般极低实时稳态测量窗函数法加窗复指数较好一般较好低PMU标准算法基础HHT自适应EMD好强较弱高非平稳动态分析小波变换小波基好强较好高暂态和振荡分析大体上说FFT和窗函数法更适合嵌入式和实时采样HHT和小波更适合离线分析或在复杂工况下做精确时频跟踪。选型从来不是一条路走到底而是看你要面对什么工况。3. Matlab 实现细节与代码落地3.1 先构造一组“会找茬”的仿真信号要对比算法信号不能太乖。我把仿真信号设计成四类工况的组合标准基波、频率偏移、谐波叠加、噪声污染。这样一套信号跑下来每类算法的问题都会被逼出来。fs 10000; % 采样率 10kHz N 2000; % 数据窗长度 0.2s正好10个工频周期 t (0:N-1) / fs; f0 50; % 标称频率 % 工况一稳态标准余弦 x1 220 * cos(2*pi*f0*t 30*pi/180); % 工况二频率偏移 50.2Hz x2 220 * cos(2*pi*50.2*t 30*pi/180); % 工况三基波 3/5/7次谐波 白噪声 harm 8 * cos(2*pi*150*t 60*pi/180) ... 5 * cos(2*pi*250*t 80*pi/180) ... 2 * cos(2*pi*350*t 100*pi/180); x3 220 * cos(2*pi*50*t 30*pi/180) harm 0.5*randn(1,N);这段代码里N取2000点在fs10000时频率分辨率是5Hz50Hz正好落在第11条谱线上。这个设置对FFT非常友好方便我先确认FFT的“舒适区”长什么样。后面一改频偏它就开始漏。3.2 FFT 相量计算主程序FFT相量提取我写成函数输入一段采样序列输出幅值、相角、频率估计。为了公平这里用峰值谱线搜索而不是固定50Hz索引因为真实系统你未必知道精确频率。function [A, phi, f_est] fft_phasor(x, fs) N length(x); X fft(x); half floor(N/2); [~, k_peak] max(abs(X(2:half))); % 避开直流 k_peak k_peak 1; A 2 * abs(X(k_peak)) / N; phi angle(X(k_peak)); f_est (k_peak-1) * fs / N; % 峰值谱线对应频率 end实际运行中只要信号频率不是正好落在谱线上这个函数的幅值误差就不可控。有一组测试里49.7Hz的信号峰值谱线落在第10条线上对应49.5Hz幅值误差约3.8%相位误差6.5度。这就是单谱线FFT在非同步采样下的典型表现。所以对FFT我建议只把它当成基准线不要作为最终方案。3.3 窗函数法代码实现加窗FFT不仅仅是乘一下窗函数那么简单关键是幅值恢复和频率插值。单谱线加窗可以把旁瓣压下去但栅栏效应还在所以我在工程中采用“Hanning窗双谱线插值”的组合。function [A, phi, f_est] window_phasor(x, fs) N length(x); w hann(N, periodic).; xw x .* w; X fft(xw); half floor(N/2); [~, k1] max(abs(X(2:half))); k1 k1 1; k2 k1 1; y1 abs(X(k1)); y2 abs(X(k2)); beta (y2 - y1) / (y1 y2); % Hanning窗双谱线修正系数经验公式 delta 1.5 * beta; f_est (k1 - 1 delta) * fs / N; % 修正幅值 A 2 * (y1 y2) / sum(w) * (0.5 0.5 * delta^2); phi angle(X(k1)) pi*delta; end这个公式属于经典处理方法适合Hanning窗。用这个函数跑50.2Hz信号幅值误差降到0.08%相位误差降到0.2度。相比裸FFT提升非常明显。要注意修正公式的适用条件是信号主瓣正好落在相邻两条谱线之间如果信号靠近边界delta会失真可以加一个合法性判断。3.4 HHT 相量估计流程HHT不适合直接一上来就对整段信号做EMD最好先做预处理。我习惯先对信号做带通预滤波把明显的高次谐波和直流分量去掉否则EMD会把高频噪声当成IMF分出来增加模态混叠风险。function [A, phi, f_est] hht_phasor(x, fs, t) % 预处理带通滤波 45-55Hz [b, a] butter(4, [45 55]/(fs/2), bandpass); xf filtfilt(b, a, x); [imf, ~] emd(xf, Display, 0); n_imf size(imf, 1); f_score zeros(1, n_imf); for i 1:n_imf analytic hilbert(imf(i,:)); phase unwrap(angle(analytic)); freq fs/(2*pi) * diff(phase); f_score(i) mean(abs(freq)); % 平均瞬时频率 end [~, idx_best] min(abs(f_score - 50)); % 选最接近50Hz的IMF analytic hilbert(imf(idx_best,:)); amp_inst abs(analytic); phase_inst unwrap(angle(analytic)); phase_ref 2*pi*50*t; phi_dev phase_inst - phase_ref; % 相对标称频率的相角差 A amp_inst(end) / sqrt(2); % 末尾时刻瞬时幅值 - RMS phi phi_dev(end); f_est f_score(idx_best); end这里最关键的认知是同步相量定义在标称频率参考系下不能直接把Hilbert相位拿过来。Hilbert相位是ωtφ其中ω是实际瞬时角频率包含频偏信息要想得到标称频率下的相量相角必须减去2π50t。这个减法做完频偏带来的相位累积会被保留在相角里才是正确的同步相量。3.5 小波变换相量估计代码小波这边我用复Morlet小波做CWT然后取50Hz附近的脊线系数。由于小波系数幅值和信号幅值之间有比例关系直接用绝对值会有标定问题。我的处理方法是用一个标准余弦信号预先标定增益再把增益代入测试信号。function [A, phi] wavelet_phasor(x, fs, t) % 标定信号已知幅值220相位30度 N length(x); t_cal (0:N-1)/fs; x_cal 220 * cos(2*pi*50*t_cal 30*pi/180); [wt_cal, fr] cwt(x_cal, amor, fs); [~, idx] min(abs(fr - 50)); coef_cal squeeze(wt_cal(idx, :)); gain 220 / mean(abs(coef_cal)); % 标定增益 [wt, ~] cwt(x, amor, fs); coef squeeze(wt(idx, :)); A gain * abs(coef(end)) / sqrt(2); phi angle(coef(end)); end这段代码思路是“相对标定”在离线分析里很实用。小波变换输出的是复数小波系数相位和cos信号相位之间有一个固定偏移标定过程可以把偏移也校准掉。但我自己做实验时发现最麻烦的是小波边界效应cwt结果的前后各一小段系数严重失真如果数据窗只有0.2s边界区域占比其实不小。用这个函数时我一般把相量输出点放在数据窗末端以外或者直接舍弃前10%和后10%的系数不用。3.6 误差评价与可视化没有误差指标算法对比就没有结论。我统一用TVE、幅值误差、相位误差三个指标。TVE把幅值和相位的误差融合成一个复平面距离是同步相量领域最常用的评价量。function TVE calc_tve(A_est, phi_est, A_true, phi_true) Xr A_est * cos(phi_est); Xi A_est * sin(phi_est); Xr0 A_true * cos(phi_true); Xi0 A_true * sin(phi_true); TVE sqrt((Xr-Xr0)^2 (Xi-Xi0)^2) / A_true * 100; end如果是在信号尾部取点而真实相量也在尾部已知那TVE就是一把很公正的尺子。画图时我喜欢用subplot把幅值误差和相位误差分开画频率偏移工况再叠加一条频率估计误差曲线这样一眼就能看出算法的问题是在幅值通道还是相位通道。4. 实测结果对比与坑点记录4.1 稳态工况窗函数是最稳的“基本盘”我先跑理想稳态信号50Hz、无谐波、无噪声。FFT和窗函数法误差都接近零TVE在0.01%以下。HHT反而有点反常TVE到了0.6%左右原因就是EMD在信号两端的包络拟合会带来瞬时幅值波动取末尾点恰好落在波动上。小波由于有标定增益和边界效应TVE约0.2%。我把这组结果放在心里提醒自己不是说算法越复杂稳态误差就一定越小。FFT在整周期采样时是数学上精确的HHT和小波反而引入额外假设。稳态工况下窗函数法就是性价比之王。给一个我实测的大致数据表方法幅值误差(%)相位误差(°)TVE(%)FFT0.0010.0020.003窗函数法0.0010.0020.003HHT0.420.280.62小波变换0.150.100.224.2 频率偏移工况HHT和小波的高光时刻把频率改成50.2Hz后裸FFT立刻现原形幅值误差4.2%相位误差6.8度TVE超过8%已经远超PMU精度要求。窗函数法因为双谱线插值能修正栅栏效应TVE降到0.3%左右。HHT和小波的表现让我意外HHT的TVE 0.8%小波0.5%虽然不如加窗插值但比裸FFT强太多。HHT能跟住频率偏移是因为它在每个时刻都给出瞬时频率并且相位减去标称参考后50.2Hz的偏移会体现为相角差按频率偏差线性增长。如果在数据窗末尾取点得到的相角正好是累积偏差算法不会因为“窗内频率非整数个周期”而产生周期性误差。小波的性能依赖小波基的频带宽度。Morlet小波带宽太窄时频率偏移会让脊线能量衰减带宽太宽又会混入噪声。我调了半天最后把频带控制参数设在中心频率的10%左右效果最好。4.3 谐波与噪声下的表现谐波是同步相量计算里的“老油条”。FFT靠固定谱线滤谐波时如果基波频率正好偏移基波谱线位置变了谐波泄漏很可能反而落在基波估计点上误差一塌糊涂。窗函数法加双谱线插值后谐波旁瓣被压低但3次谐波150Hz距离基波50Hz较远影响不大如果谐波是49Hz和51Hz的间谐波那再好的窗也没办法因为主瓣本身就重叠。HHT在谐波存在时模态混叠非常明显。有一次我把3次谐波幅值调到基波的10%EMD居然在前几个IMF里同时出现了50Hz和150Hz的混叠瞬时频率不再干净。后来我用EEMD集合经验模态分解加辅助白噪声效果好了些但EEMD计算量成倍增长而且辅助白噪声幅值还要调。小波对谐波的抑制能力其实不错因为CWT在频率维展开后150Hz的系数能量集中在150Hz脊线上50Hz附近脊线受影响有限。但噪声一上来小波脊线的斑驳感就出来了单点系数抖动大需要对脊线做平滑。我做了一组含谐波0.1%噪声的实验结果是窗函数法TVE 0.8%小波1.2%HHT 2.5%上下选取平滑后的包络。如果你想在谐波和噪声环境下做稳定估计窗函数法依然是首选。4.4 动态事件阶跃和振荡真正让DFT类方法头疼的是动态事件。我在信号里叠加了一个幅值阶跃在0.1s处幅值从220跳到250。FFT和窗函数法因为数据窗包含事件前后两段幅值估计会出现平滑过渡上升沿被拉长准确响应时间近似数据窗长度N/fs0.1s。HHT的瞬时包络能比较快地反映出阶跃但EMD在阶跃点周围会画出过冲和振铃不要以为瞬时方法就没有延迟。相位调制场景更有意思。我在50Hz基波上叠加了一个10Hz的相位调制模拟动态振荡。HHT的瞬时相位能清楚解调出相位振荡波形小波脊线也能做到这两者在动态分析上的价值确实不可替代。窗函数法只能给出一个窗内的平均相量对快速振荡无能为力。所以如果你做的是低频振荡分析、次同步振荡监测不要只用FFT建议把HHT和小波加进来。4.5 常见问题与排查速查表现象可能原因排查思路FFT幅值误差大频率偏移导致频谱泄漏改用加窗插值或用锁相环预跟踪频率加窗后相位突变双谱线插值delta计算错误检查k1、k2索引确保没有跨直流或跨主瓣边界EMD出现模态混叠分量频率接近或幅值差异大改用EEMD/CEEMDAN或先做带通预分离HHT端点振荡严重EMD样条拟合边界无约束数据窗两端各丢弃10%样本或增加镜像延拓小波边界系数失真CWT在时频平面边界的锥形效应输出点避开边界或加短时窗后舍弃边缘代码运行极慢EMD和CWT本身计算量大降低采样率到2k-4k减少数据窗点数TVE突然大于0.5%数据窗内存在动态事件缩短数据窗或改用瞬时类算法5. 工程落地建议5.1 按应用场景选算法如果你是要做一台实时PMU装置MCU或DSP资源有限窗函数法双谱线插值基本就是最优解。它在频率偏移、谐波、噪声下都能保证较高精度计算量比FFT多一点点嵌入式完全能扛。HHT和小波更适合离线分析或者上位机后处理它们能提供FFT给不了的瞬时动态信息。我的一个明确建议是不要在一套算法里指望“什么都强”。你可以在实时通道用窗函数法输出标准相量同时保留原始采样点触发现场数据异常时再用小波或HHT做事件追忆。这种“实时基座离线分析”的分层架构工程上比单算法通吃可靠得多。5.2 参数经验值采样率我建议至少2000Hz实际做研究用10kHz没问题工程上4kHz左右比较合适。数据窗长度稳态精度要求高选10个周波动态响应要求高选4-6个周波。窗函数优先Hanning不要一上来用Blackman旁瓣虽然压得低但主瓣太宽容易把相邻间谐波糊在一起。小波变换的频带范围做基波相量我习惯设45-55Hz太宽噪声多太窄频率偏移时脊线会跑出去。HHT做预滤波时同样用45-55Hz带通可以显著减少EMD的高频分解负担。TVE阈值我通常按标准要求的1%来卡但在动态工况下不要对瞬时误差太苛刻要看平均TVE和最大TVE的统计值。5.3 我对这套研究的个人体会做这个课题我一开始满脑子都是“FFT不够就上HHT再不行小波”结果数据结构一团糟光洗数据就洗了一周。后来我先把仿真信号按稳态、频偏、谐波、动态四类列成表格再让每个算法在同一组信号上跑把所有输出统一转成幅值、相位、频率三个量用一套评价函数计算TVE对比才真正可操作。这是我这次项目里最值得分享的经验算法对比的第一步不是跑代码而是定好统一的“尺子”。还有一个小技巧无论是窗函数插值还是HHT最后取相量值时都别直接用数据窗最后一个点。我在窗函数法末尾取点遇到过相位轻微跳动后来往前挪5-10个采样点取平均稳定性好了很多。同步相量计算最怕单点抖动适当的“滞后输出”在工程上是允许的只要不超过标准要求。如果你后续想继续扩展可以试试把EEMD和小波脊线结合先用小波变换确定基波频率搜索范围再做约束EMD这样模态混叠会明显减少。再往后如果要做多台机同步测量把每路的估计相量打上UTC时标分析相量轨迹那就已经是完整的广域测量系统雏形了。算法这条路工具永远是死的思路才是活的。
阅读完成 · 觉得有帮助？
