基于NILM的污染源监测方法应用初探

刘春蕾,庞鹏飞,石纹赫,胡伟俊,戚 军

(1. 国网河北省电力有限公司保定供电分公司，河北 保定 071000;2. 浙江工业大学信息工程学院，浙江 杭州 310023)

环境污染降低了自然资源的利用率，影响人的身体健康[1]，环境监测是抑制环境污染事件发生的重要一环。我国现阶段污染监测技术主要分为地表水系统监测、空气质量监测以及污染源监测三种形式，目前都存在一些问题，例如：地表水系统监测对于监测位置选择具有较高要求，城市中空气质量监测站数量尚无法满足监测准确性和完整性要求，安装于每个排污设备端和治污设备端的污染源监测仪不仅实施代价高而且存在作弊问题[2]。考虑到很多排污设备和治污设备都由电力驱动，以及设备运行时电压电流数据通常无法造假的特点，本文尝试将NILM技术应用于污染源监测中，通过实时监测排污企业入户电力数据来获取排污设备和治污设备的启停信息，从而在确保监测实时性和准确性的同时降低污染源监测的难度。

NILM由Hart博士于20世纪80年代提出，随着测量技术的发展和算法研究的深入[3-4]，NILM技术出现在广阔的应用领域，包括：需求侧管理中为电价与供需联动提供基础数据[5]，供电监测中实现异常负荷监测以及故障定位[6]，用电行为分析中提升用户体验和资源利用率[7]。NILM中设备投切事件检测的准确性和快速性对于后续负荷特征提取以及类型识别具有显著影响。早期研究中使用较多的启发式方法及匹配滤波方法原理都较为简单，在判断事件是否发生过程中，启发式方法关注负荷的变化率是否超过阈值，匹配滤波方法则将采集信号与已知模板进行匹配，两者仅适用于设备种类较少且干扰较小的应用场景。基于概率模型的变点检测方法是当前的研究热点，Rehman等[8]提出用事件通用似然比(Generalized likelihood ratio，GLR)来检测事件，其通过调整滑动窗口的长度、功率方差判定阈值等参数来提高算法的精度；Anderson等[9]在GLR法的基础上引入了投票机制，提高GLR法的准确性；史帅彬等[10]提出了一种基于复合滑动窗的累积和(Cumulative Sum，CUSUM)暂态事件检测算法，能够准确检测到容易被忽略的具有缓慢爬坡特性的暂态事件。由于需要采集足够长的数据才能作出较为准确的判断，变点检测方法很难及时定位事件发生时刻，从而使测量和存储高频负荷数据变得困难。

本文以污染源监测需求为出发点，提出并设计了基于NILM技术的污染源监测方法及实现流程，主要包括基于小波滤波算法的实测数据去噪环节、基于最小二乘预测法(Least-square prediction method，LSPM)的事件检测环节以及基于相位和频率矫正的负荷特征提取环节。最后在BLUED数据库和实际污染源设备监测工程中分别进行仿真验证和实际测试。

1 基于NILM的污染源监测系统

基于NILM的污染源监测系统主要由工业用户端、数据采集卡、负荷分析仪、电网公司以及相关环保部门组成，具体如图1所示。电网公司在工业用户的入口端安装数据采集卡采集电力信号；数据采集卡中的信息被传输到负荷分析仪；负荷分析仪完成负荷分析计算并实现污染源监测，将监测结果通过网络上传到云端；电网公司以及环保部门能够读取云端数据进行数据可视化展示，最后由环保部门对工业用户的排污和治污行为作进一步监管。

图1 NILM污染源监测系统结构

基于NILM的污染源监测流程主要包括数据采集、小波去噪、事件检测以及特征提取等环节，详见图2。其中，数据采集环节需在工厂入户端加装高频(建议10kHz及以上)数据采集卡，实时测量入户电力线路上的电压和电流数据。由于传感、测量、转换和传输等过程存在各种干扰，采集到的电压电流波形难免带有毛刺或畸变，直接影响后续事件检测及特征提取环节的准确性[11]，因此在小波去噪环节将对采集到的高频电压和电流信号进行在线滤波以减少噪声干扰影响。在事件检测环节根据滤波后电流的变化情况，采用LSPM法来判断是否有事件发生。一旦检测到事件，将存储事件前后总线电压和电流波形数据。在负荷特征提取环节将基于相位和频率矫正技术从总线波形中获取投切设备电压和电流波形，并进一步获取波形所包含的负荷特征信息。根据详细的负荷特征信息，比对历史负荷特征库，即可判断出当前启停事件对应的设备类型，从而实现对污染源的实时在线监测。

1.1 小波去噪环节

去噪环节中采用的是小波滤波法，相较于中值滤波与均值滤波法无法保留信号中突变信息的缺点，小波滤波法具有突变信息保留更完整、波形更稳定的优点。此外，小波滤波法效率高、原理简单，基本能满足工程使用中的实时性要求。小波滤波法具体步骤如下：

步骤1) 小波分解：选择Symlet系小波基，假设分解层数总数为m层，i表示第i层分解层数，预先设置i=1，第1层输入信号x(t)由含噪声的总线电压或电流测量信号构成。

图2 污染源监测流程

步骤1-1) 对信号x(t)进行分解，分解公式如式(1)所示。

(1)

其中k表示函数位移量，x(t)为尺度函数族φi,k(t)和Symlet系小波函数族ψi,k(t)的线性组合，低频分量ci[k]与高频分量di[k]的计算公式如式(2)和式(3)所示[12]。

ci[k]=〈x(t),φi,k(t)〉

(2)

di[k]=〈x(t),ψi,k(t)〉

(3)

其中<>表示内积计算。对x(t)完成分解后，令i=i+1，并转入步骤1-2)。

步骤1-2) 若i≤m，进入步骤1-3)继续下一层小波分解；否则，结束小波分解，进入步骤1-4)。

步骤1-3) 更新第i层的分解信号x(t)为ci-1，然后返回步骤1-1)。

步骤1-4) 分解结束得到由m层高频分量与1层低频分量所组成的小波系数{d1,…,dm,cm}，第i层系数对应的长度为n/2i，n为原始信号长度。

步骤2) 阈值选择：阈值选取直接影响到重构信号的还原程度，本文采用通用阈值规则进行阈值选取，阈值T的计算公式如下[13]：

(4)

其中σ为噪声标准方差，计算公式如式(5)所示。

(5)

其中Median函数为中值求解函数。

步骤3) 阈值去噪：由于信号中大部分噪声存在于高频信号中，因此需要通过阈值函数对每一层高频分量di去噪。由于常用的软、硬阈值函数对信号去噪后存在重构信号偏差较大的缺点[14]，本文选择软硬阈值折衷法，其表达式如下[15]：

(6)

其中sgn函数在输入值大于0时返回1，等于0时返回0，小于0时返回-1。软硬阈值折衷法通过在0与1之间适当调整α的大小，获得较好去噪效果，文中取α=0.5，对高频分量中小于阈值T的部分进行过滤，大于阈值T的部分保留并进行收缩，得到新的高频分量d′i。

步骤4) 小波重构：将去噪后的高频分量d′i与低频分量cm进行离散逆变换就可以将原始信号重构，得到去噪后的负荷信号f′(t)，信号重构公式如式(7)所示。

(7)

1.2 LSPM事件检测环节

具体实现步骤如下：

步骤1) 基于滤波后的电流波形提取工频周期最大值，其中W1的电流最大值序列为{Imax(k-m),Imax(k-m+1),…,Imax(k-1)}。

步骤2) 在检测窗口中构造一条待拟合曲线，函数表达式如下所示：

Y=Xθ

(8)

其中，待拟合直线对应的系数θ∈R2×1，输入值X∈Rm×2，实测值Y∈Rm×1，矩阵的表达式如式(9)和式(10)所示。

(9)

(10)

图3 LSPM检测算法示意图

根据构造的曲线，能够得到其误差函数如式(11)所示。

J(θ)=(Xθ-Y)T(Xθ-Y)=

θTXTXθ-θTXTY-YTXθ+YTY

(11)

步骤3) 为求得使误差函数最小的系数θ，对误差函数中的J(θ)求偏导，如式(12)所示。

(12)

求导等式进行整理之后能够得到待拟合曲线系数的最优解θ，如式(13)所示。

θ=(XTX)-1XTY

(13)

(14)

在检测到C(k)超出阈值设定范围时，判定采样点k处发生设备投切事件，转入步骤5)，否则返回步骤1)，将W1的起始点更新为k-m+1，通过新的电流最大值序列继续进行事件检测。

步骤5) 根据W2中数据确定投切设备由暂态进入稳态的时间点，判断依据如式(15)所示。

(15)

步骤6) 假设设备在Imax(x)序列中第k+l个采样点处结束暂态，l为设备暂态过程周期个数。从投切前稳态起始采样点k-l，暂态起始采样点k以及投切后新稳态的起始采样点k+l开始分别保存三个时间段的电压与电流波形，并分别命名为upre-event，ipre-event，utran，itran，upost-event，ipost-event，每段波形的长度都为Nl，N为工频下每个周期数据长度。

1.3 负荷特征提取环节

由于用电设备都是并联连接，且供电线路上的压降通常较小，因此用电设备的端电压可以近似于总线测量节点上的电压。总线上的电流信号由多个设备的电流信号叠加而成，负荷m投切后总线上电流信号itran(t)由投切前总线电流信号ipre-event(t)与用电负荷m的电流信号iload(t)叠加而成，总线电流和每个用电负荷电流信号的表达式都如式(16)所示。

(16)

其中，f为电网频率，ψ为电流信号相位角，A为电流信号幅值，x为信号中谐波次数。

理论上从itran(t)与ipost-event(t)中剔除ipre-event(t)即可提取投切设备的暂-稳态电流波形iload(t)，但是实际离散测量信号面临电网频率波动带来的相位偏差以及周期采样点数量的变化，将电流信号简单地相减会使电流波形出现误差。因此本文提出一种相位和频率矫正方法来获取较为准确的电流信号。特征提取环节具体实现步骤如下：

步骤1) 相位矫正：由于电流信号相角由电压的起始相位决定，因此需要保证三个时段的电流信号对应的电压波形在同一起始相角下。当upre-event第一个周期中的采样点upre-event(m)满足式(17)时，则认为该点为电压过零点，并将ipre-event的稳态起始点更新为m，ipre-event最后一个周期的结束点同样更新为对应的电压过零点。

(17)

同理，通过式(17)将itran与ipost-event的起始与结束采样点更新为utran与upost-event第一和最后一个周期中的电压过零点，与此同时将upre-event，utran，upost-event的起始与结束点也更新为电压过零点。

步骤2) 频率判断：通过实测信号数据计算出三个时段对应的电网频率，计算公式如下所示：

(18)

其中fs为NILM装置采样频率，本文设置为10kHz,l′为相位矫正后波形中的周期个数，n为矫正后波形长度。若实测频率与额定频率偏差不超过0.1%，直接转入步骤4)提取投切设备的电压与电流波形，否则转入步骤3)进行频率矫正。

步骤3) 频率矫正：由于不同时段的电压和电流信号频率矫正过程相同，因此将待矫正信号统一记为Ypre-cor(k)，k=1,2,…,n，将采样时刻根据实测频率f与额定频率f0的比值进行缩放，第k个采样点对应时刻t(k)的计算公式如下：

(19)

将t(0)与t(k)之间时段按照1/fs等间隔采样，记为t′(j)，j=1,2,…,Nl′，与t′(j)对应的Ycor(j)即为频率矫正后的信号，计算公式如式(20)所示。

Ycor(j)=

(20)

步骤4) 波形提取：相位及频率矫正后的itran，ipost-event与ipre-event分别相减能够得到投切设备的暂态电流波形iload -tran以及稳态电流波形iload-steady，负荷电压波形uload-tran和uload-steady与总线电压utran和upost-event相同，如式(21)～式(24)所示。

iload-tran=itran-ipre-event

(21)

iload-steady=ipost-event-ipre-event

(22)

uload-tran=uload

(23)

uload-steady=upost-event

(24)

步骤5) 特征计算：负荷特征包括暂态时长Ttran、稳态有功功率P和无功功率Q、暂态电流绝对平均值|I|mean、有效值Irms，最大峰值Ipp以及电流谐波的幅值Ax与相角φx，其中x为信号谐波次数[16]。

步骤6) 负荷类型判断：通过与历史负荷特征数据库进行对比，判断当前检测到的事件所属设备类别。

2 BLUED数据库验证

2.1 去噪方法仿真分析

本文选择使用BLUED数据库[17]中电力数据来验证该监测技术的准确性及有效性。首先选择BLUED数据库中A相电流信号2个周期的采样数据共计400个采样点，分别采用小波滤波、均值滤波以及中值滤波进行去噪处理。去噪前后的波形如图4所示，与中值滤波及均值滤波法相比，小波滤波法去掉了大部分尖峰，使电流的波形变得较为平滑，并且较好地保留了波形原始信息。

图4 去噪前后波形对比

2.2 事件检测仿真结果分析

本文采用召回率Precall、精确率Pprecision、F分数Fm以及事件检测时间Tavr作为事件检测的评判指标[18]，其定义如式(25)～式(27)所示。

(25)

(26)

(27)

其中TP(true positive)为真阳性，TN(true negative)为真阴性，FP(false positive)为假阳性，FN(false negative)为假阴性。Precall为正确检测到阳性事件占所有实际阳性事件的比例，Pprecision为正确检测到阳性事件占全部检测结果为阳性事件的比例，Fm为Precall和Pprecision的调和平均值，用于综合考虑精确率和召回率，Tavr用来衡量完成一次事件检测所需要的平均时间。

选取BLUED数据库中2011年10月20～21日A相24小时的实际测量数据进行事件检测仿真对比，投切事件实际发生的总次数为132次。为验证去噪方法对事件检测效果的影响，本文加入数据滤波前的检测结果作为对比，为考察本文所使用的LSPM法的性能，增加GLR法以及滑动窗法[19]作为对比实验，三种检测方法的结果如表1所示。

表1 事件检测结果对比

由表1可知，由于滤波需要花费一定的计算时间，滤波前事件检测的Tavr都小于滤波后的Tavr；由于滤波前数据附带的许多尖峰与毛刺给事件检测带来误差，造成LSPM法和滑动窗法检测的准确率显著下降。其中GLR法运用有功功率曲线进行事件检测，而有功功率计算过程相当于进行了均值滤波，因此导致滤波前后其事件检测结果并无明显差别。

从小波去噪之后的LSPM法、GLR法以及滑动窗法的性能对比可以看出，本文所提出的LSPM法的Precall、Pprecision及Fm基本与GLR法相当，但显著优于滑动窗法；从检测时间对比来看，LSPM法的Tavr明显短于GLR法，该检测方法在事件发生后能够快速反映事件检测结果。因此下文的实际污染源设备监测工程中，选用准确性与GLR相当，但检测速度明显优于GLR的LSPM法来检测排污设备和治污设备的启停。

2.3 特征提取仿真分析

选取BLUED数据库中A相某一时刻冰箱开启前后的电流波形进行特征提取。如图5所示，分别为相位与频率矫正前与矫正后提取的设备电流和电压波形。由矫正前后波形对比图可知，矫正后的电流波形更为平滑且幅值更为准确，由此可见相位及频率偏差对于所提取的设备电流波形影响非常大。根据提取出的冰箱设备的电压电流波形，计算得到如表2所示的用电设备特征，由表2可知，矫正后用电设备的P和Q与准确值(即设备投切前后总线稳态功率之差)更为接近。由此可见提取设备电流波形前进行相位与频率矫正的必要性。

3 实际污染源设备监测工程验证

对某工厂电压和电流信号采用小波滤波进行去噪处理，其中5个周波电流波形去噪前后的对比如图6所示。小波去噪之后过滤掉大部分噪声，总线上电流波形更为平滑，并且较好地保留了信号的真实特征，有助于提高后续事件检测及特征提取的精度，准确定位发生事件对应的时刻。

图5 矫正前后冰箱电流波形

设定事件检测阈值为200A，对2021年5月2～9日这一周的事件检测结果进行分析。这一周中发生的投切事件次数为84次，其中三相事件检测结果如表3所示。可以看出三相电流波形检测的召回率能够达到88%以上，精确率达到了90%以上，说明该方法运用在实际设备的事件检测中漏检率较低。考虑到事件发生频次可能较高，NILM监测系统需要在较短时间内完成事件检测，如表1中Tavr对比所示，LSPM法的检测速度能够满足实时分析需求，实际运行结果也证实了这一点。

表2 相位和频率矫正前后的负荷特征

图6 去噪前后工厂电流波形对比

表3 工厂事件检测结果

经相位与频率矫正之后提取某设备电流波形如图7所示。在设备启动瞬间，设备的电流峰值接近900A，随后电流开始缓慢下降，最终在20s左右达到稳态，稳态时设备电流峰值稳定在240A左右。根据所提取的设备电流波形计算得到负荷特征如表4所示。对比历史负荷特征库，即可判断出当前投切设备所对应的设备类型，从而判断该工厂实时的排污和治污状况，为环保监管提供及时准确的现场信息。

图7 投切设备的电流波形

表4 设备负荷特征

4 结 论

本文设计了基于NILM的污染源监测系统，其采用小波滤波法对所采集的电力数据进行去噪，采用LSPM法对总线电流数据进行事件检测，采用相位和频率矫正法对总线上电压和电流波形进行矫正以获取用电设备特征。最后在BLUED数据库和实际工程案例中进行了详细的测试和验证，结果表明：本文设计的污染源监测系统中使用的小波滤波法具有效率高、原理简单的优势；提出的LSPM法具有较高的召回率、准确性以及事件检测时效性；经过相位和频率矫正后的电流波形更为准确，提取的特征能较好满足污染源监测的需求。本文的研究为污染源实时监测提供了新的技术视角和实现途径，具有较好的应用前景。

本文主要着眼于污染源设备投切事件检测和负荷特征提取，但并未对设备类型识别作深入探讨，后续可在本文基础上应用机器学习算法对此作进一步研究。