引言
在“碳达峰”和“碳中和”以及新型电力系统建设的大背景下，我国能源需求不断增加。随着轧钢、电弧炼钢等非线性负荷的大量应用[1]，负荷的电流呈现出复杂宽频带谐波扰动的特点。电流互感器作为负荷电流变换采集的设备，其传变特性是影响电能计量综合误差的主要因素。因此，围绕电流互感器宽频带传变特性的快速准确辨识方法进行研究，对于提高电能计量的准确性具有重要意义。
为解决CT传变特性时频域的快速准确系统辨识的理论问题，本文首先推导了CT的Jiles-Atherton (J-A) 模型传递函数，其次提出了CT压缩感知系统辨识架构，构造了CT系统传变函数的多频正弦信号稀疏字典以产生压缩感知辨识激励信号，随后构建了Toeplitz降维测量矩阵，通过压缩采样匹配追踪算法辨识CT传变特性。最后通过仿真实验对压缩感知系统辨识方法的性能进行了验证。
1 J-A模型电流互感器传递函数的推导
为了进一步推导CT的传递函数，结合图 2 ，分别从一次侧、二次侧频域部分进行传递函数的推导。根据图 2(a) 可以得到CT一次侧传递函数为
式中，Z1为一次绕组阻抗且Z1=R1+sL1。当C1很小时，G1(s)=1，表示电流信号在一次侧幅值没有衰减，相位没有偏移。
在CT为理想变量器条件下，根据CT链式守恒原理N1×I1=N2×I2，一次侧、二次侧等效电路间的关系如式(2)所示。
式中，K=N2/N1。根据图 2(b) ，得到二次侧部分CT传递函数为
式中，$Z_2=\frac{\frac{1}{s C_2} \times R_3}{\frac{1}{s C_2}+R_3}+\left(R_2+s L_2\right), Z_{\mathrm{e}}=R_{\mathrm{e}}+s L_{\mathrm{e}} $ ，将其代入式(3)中可得
因此，结合CT一次侧和二次侧部分可以得到CT的传递函数为
根据式(1)、(4)、(5)，CT的传递函数可以表示为
式中，$R_{\mathrm{e}}+R_2 \gg R_3, R_{\mathrm{e}} \gg R_2, L_{\mathrm{e}} \gg L_2, \omega_0=R_{\mathrm{e}} / L_{\mathrm{e}} $ , $P=L_{\mathrm{e}} / R_3 C_2\left(L_2+L_{\mathrm{e}}\right), \omega_1=1 / R_3 C_2, \omega_2=\left(R_2+R_{\mathrm{e}}\right) /$ $\left(L_2+L_{\mathrm{e}}\right) $ 。式(6)表明，传统电流互感器的传递函数为具有1个单零点和两个极点的线性系统。
2 压缩感知重构
式中，s∈RN×1为x在变换域Ψ中的稀疏信号。对于信号s，其压缩采样过程如式(8)所示。
式中，y∈RM×1(M≪N)为观测采样值，Φ为测量矩阵，θ为传感矩阵。根据式(8)可将压缩感知重构过程表示成如图 3 所示。
同时相应的测量矩阵Φ需满足有限等距特性(restricted isometry property, RIP)[17]。对于K稀疏的信号s∈RN×1和常数δk∈(0, 1)，满足
在满足RIP准则的基础上，压缩感知重构过程可转化为求解欠定方程组的问题，即求解最小l0范数的问题，结合式(8)可表示为
最终，通过采样信号y求出稀疏信号 $\hat{s} $ ，再通过$\hat{x}=\boldsymbol{\varPsi}^{-1} \hat{s} $ 求出最终的重构信号$ \hat{x}$ 。
为进行CT传变特性系统辨识，结合压缩感知理论，可将CT看作一个参数客观存在的待辨识系统，而传感矩阵作为该系统的输入信号，由此提出CT频域传变特性压缩感知系统辨识架构。
3 电流互感器频域传变特性的压缩感知辨识架构与辨识方法
3.1 构建电流互感器传变函数的稀疏字典
在CT系统辨识过程中，由于无法直接采用压缩感知测量矩阵作为输入信号进行传变特性的辨识，故需要将系统的输入输出调整为电流信号。为解决CT传变特性的系统辨识问题，提出压缩感知系统辨识架构，如图 4 所示。
根据图 4 ，构造一种多频率正弦信号稀疏字典，产生压缩感知辨识信号，输入CT系统，再根据系统输出信号，基于所提出的架构完成CT系统传变特性的压缩感知辨识。
首先构造多频率正弦信号，设定长度为M点的截断离散单频正弦信号如式(11)所示。
式中，fs为信号采样频率，各频率分量下的点数m=1, 2, …, M，fl=l×f1且当l=1, 2, …, N时，得到N个不同单一频率且长度为M的离散正弦信号。定义bml=s(tm, l)，由式(11)进一步推导可得
式中，ω1=2π(f1/fs)表示初始角频率。
其次，采用N个M点截断离散单频正弦信号构建矩阵如式(13)所示。
式中，矩阵LM×N中第l列表示角频率为ωl的离散正弦信号。
然后，构建CT的传变函数向量如式(14)所示。
式(14)表示待辨识的CT系统。由式(13)、(14)可知所构建的矩阵LM×N包含了CT系统的各个频率分量，可以保证系统解的稀疏性。因此，矩阵LM×N为CT传变函数的稀疏字典，满足压缩感知稀疏条件。
3.2 电流互感器频域传变特性压缩感知辨识过程
为了实现电流互感器的精确辨识，同时降低压缩感知硬件的复杂度，根据图 4 ，考虑实际部署的便捷性，采用Toeplitz矩阵Φ作为测量矩阵，且当信号稀疏度K≤c1M3/ln (N/M)(N、M分别为测量矩阵的行、列数，c1为常数)时，Toeplitz矩阵满足RIP准则[18]。
根据式(13)、(14)，CT的压缩感知时频域系统输出信号为
根据式(15)，可将CT传变函数HN×1的辨识问题拟合为压缩感知重构问题，根据式(8)，可将CT二次侧的输出信号进一步写为
其中，ΦM×M为Toeplitz测量矩阵。由于式(16)为一个NP难题，无法直接求出稀疏向量HN×1的值，根据式(10)，可将式(16)写为求解最小l0范数的问题
则可以通过求解式(17)非零元素的值，再根据构造的压缩感知稀疏字典矩阵LM×N确定HN×1中非零元素所对应的频率分量，进而拟合出CT的传变函数和二次侧电流有效值。所采用的构造压缩感知矩阵来表示辨识激励信号的方法，有效地避免了常规压缩感知理论中使用压缩感知采样电路的难题。
3.3 电流互感器频域传变特性的压缩感知重构算法
(1) 输入传感矩阵θ=ΦM×MLM×N，采样信号y(y=YM×1′)以及稀疏度K；输出系数σK∈HN×1。
(2) 初始化设定：t=1，残差ro=y，索引集Λ0=Ø，候选原子集合ψ0=Ø。
(3) 结合传感矩阵θ与残差r内积，选择最大的2K个原子对应的索引集合，即$ \varLambda_k=\arg \max\limits_i\left|\boldsymbol{a}_i^{\mathrm{T}} \boldsymbol{r}_{k-1}\right|$ , 2K。
(4) 更新候选原子集合ψk，并在ψk与rk－1的内积中选择绝对最大的K个原子的索引集合，即$ \varLambda_k=\arg \max\limits_i\left|\boldsymbol{\psi}_i^{\mathrm{T}} \boldsymbol{r}_{k-1}\right|$ , K。
(5) 利用最小二乘法求出系数$\sigma^k=\boldsymbol{\theta}_{\varLambda_{k t}}^{\dagger} \boldsymbol{r}_{k-1} $ 。
(6) 更新残差并令t=t+1，如果t>K或残差rk=0则停止，否则返回步骤(3)。
4 电流互感器频域传变特性压缩感知辨识性能分析
4.1 扫频仿真实验与结果分析
由式(6)可知，CT为包含单零点和两个极点的线性系统。基于PSCAD软件平台中电流互感器的J-A模型，结合式(6)电流互感器的传递函数，根据电流互感器的实际参数情况，设定PSCAD软件平台中电流互感器的额定参数，具体如表 1 所示。建立J-A模型电流互感器作为辨识用电流互感器。
表 1 电流互感器额定参数Table 1 Rated parameters of the current transformer |
| 参数 | 设定值 |
| 输入电流/A | 100 |
| 线圈匝数比 | 20/200 |
| 铁心长度/m | 0.637 7 |
| 铁心横截面积/cm2 | 30 |
| 二次绕组电阻/Ω | 0.5 |
| 二次侧漏感/H | 0.8×10-3 |
| 二次负载电阻/Ω | 20 |
扫频系统辨识过程主要通过调整CT一次侧的电流信号频率依次得到不同频率下的二次侧电流信号，再通过Matlab计算得到CT的传变特性。CT的传变特性包括比差传变特性和角差传变特性，计算方法如下。
1) 比差
式中，K为CT一次侧和二次侧电流信号的额定变比；I1为CT一次侧实际输入的电流；I2为对应测量得到的二次侧电流。
2) 角差
式中，ϕ1为CT一次侧实际输入的电流信号相位；ϕ2为对应测量得到的二次侧电流信号相位。
结合式(18)、(19)，扫频系统辨识流程如图 5 所示。
在仿真辨识中，保持CT一次侧电流幅值I1n为100 A(I1n表示不同单一频率下的一次侧电流信号，n=1, 2, …, N)，根据图 5 ，通过调整CT一次侧电流频率(10~3 500 Hz)，利用PSCAD中的CT系统得到二次侧电流，经Matlab计算后，得到CT的比差、角差传变特性如图 6 所示。
从图 6 可以看出，J-A模型电流互感器的通带频率范围为25~2 500 Hz，低于或超出此频率范围，电流互感器将逐步产生大量误差。
4.2 压缩感知仿真实验与结果分析
为评价CT传变特性压缩感知系统辨识的性能及准确度，引入以下3个参数。
1) 各频率分量下CT频域传变特性重构后的相对误差，计算公式如式(20)所示。
式中，J表示原系统参数，J*表示辨识得到的系统参数。
2) 为评价CT压缩感知系统辨识方法的准确度，采用压缩感知评价方法，计算CT压缩感知系统辨识方法与扫频系统辨识方法二次侧输出电流有效值间的均方误差，计算公式如式(21)所示。
式中，$\hat{I}_{2 n} $ 表示不同频率分量下CT压缩感知辨识方法二次侧输出电流信号的有效值，N表示频率分量数。由于CT的系统辨识不涉及压缩系数和信噪比，因此本文不使用其作为评判指标。
3) CT压缩感知系统辨识方法和扫频系统辨识方法所需时间对比。
压缩采样匹配追踪(CoSaMP)算法下CT传变特性的系统辨识结果如图 7 所示。
上述辨识结果表明，压缩感知辨识误差与扫频辨识误差处于同一量级，扫频辨识传变特性曲线与压缩感知时频域传变特性辨识曲线也具有非常好的一致性。
在谐波电能测量频率范围内，整数次谐波频率点下的CT扫频系统辨识方法和压缩感知系统辨识方法的二次侧输出电流幅值如表 2 所示。
表 2 CT系统二次侧输出电流幅值Table 2 Output amplitude of the secondary side of the CT system |
| 频率/Hz | 电流幅值/A | 频率/Hz | 电流幅值/A | |||
| 扫频辨识 | 压缩感知辨识 | 扫频辨识 | 压缩感知辨识 | |||
| 25 | 9.994 7 | 9.999 7 | 550 | 9.997 7 | 9.996 3 | |
| 50 | 9.997 7 | 9.998 7 | 650 | 9.997 5 | 9.997 2 | |
| 100 | 9.997 6 | 9.996 8 | 750 | 9.997 3 | 9.995 3 | |
| 150 | 9.997 5 | 9.998 3 | 850 | 9.997 7 | 9.995 7 | |
| 200 | 9.997 7 | 9.996 0 | 950 | 9.997 5 | 9.996 6 | |
| 250 | 9.997 4 | 9.999 0 | 1 100 | 9.997 7 | 9.998 2 | |
| 300 | 9.997 7 | 9.996 0 | 1 500 | 9.995 2 | 9.997 3 | |
| 350 | 9.997 6 | 9.996 7 | 2 000 | 9.989 2 | 9.990 7 | |
| 400 | 9.997 4 | 9.998 7 | 2 500 | 9.896 2 | 9.906 6 | |
| 450 | 9.997 7 | 9.998 0 | 2 550 | 9.884 7 | 9.886 7 | |
根据式(21)和表 2 中CT系统二次侧输出幅值，计算得到压缩感知系统辨识方法二次侧输出电流有效值的均方误差为0.018%，能够很好地满足0.5级CT系统辨识的误差要求，表明CT传变特性的压缩感知系统辨识方法具有较高的准确度。
在谐波电能测量频率范围内，以50 Hz间隔扫频，分别计算CT扫频系统辨识方法和压缩感知系统辨识方法的运行时间。其中，扫频系统辨识方法所需时间为1 560 s，压缩感知系统辨识方法所需时间为108.948 s，是扫频系统辨识方法的1/14，说明CT传变特性的辨识效率得到提高。
5 结论
(1) 构造了一种多频率正弦信号稀疏字典，产生压缩感知辨识激励信号，在此基础上提出了CT压缩感知系统辨识架构。该架构通过多频率正弦信号稀疏字典产生CT系统的输入信号，满足了CT传变特性的辨识要求和压缩感知重构的先验条件。
(2) 提出了CT传变特性的时频域压缩感知准确辨识方法。在谐波电能测量的频率范围内，比差传变特性压缩感知辨识方法各频率点辨识相对误差平均值为0.06%，角差传变特性压缩感知辨识方法各频率点辨识相对误差平均值为0.49′，表明CT传变特性的压缩感知系统辨识方法具有较高的准确度。
(3) 仿真实验结果表明，所提的压缩感知系统辨识方法所需时间是扫频系统辨识方法的1/14，提高了CT传变特性的辨识效率。