基于HCM 的光伏反窃电算法研究

2020-03-20陈海峰应国德潘成峰

浙江电力 2020年2期

陈海峰，应国德，曹杰，林超，潘成峰，金潮

（1.国网浙江温岭市供电有限公司，浙江温岭 317500；2.国网浙江省电力有限公司台州供电公司，浙江台州 318020；3.温岭市非普电气有限公司，浙江温岭 317500)

0 引言

由于分布式光伏发电享受的补贴主要取决于自身发电量，在利益的驱使下，某些用户通过技术手段使分布式光伏上网电表多计量发电量，进而获取高额补贴，这种骗取补贴的行为称为光伏窃电行为。光伏窃电行为严重影响了中国新能源扶植政策的落实和发电市场的公正性，用户因窃电私自改接线路给供配电带来了巨大的安全隐患，影响了光伏发电行业的正常发展。

光伏窃电是一个较新的课题，目前研究较少。但对传统反窃电手段进行分析，可以为本课题的研究提供思路。目前国内外的反窃电手段主要包括仪表自身数据排查、窃电识别技术、大数据分析技术[1-2]。

仪表自身数据排查主要是从仪器数据自身的层面对数据进行简单分析，排查可能的窃电行为，这是目前研究和应用都比较多的策略。詹乔松对常见的两种窃电接线方式进行了分析，针对这两种窃电接线方式制订了分析策略。这种策略虽然是针对传统窃电行为的，但依据其逆向思路也可以发现光伏窃电[3]。

窃电识别技术是通过人工神经网络、支持向量机、逻辑识别等人工智能算法来发现窃电行为。胡林等对高压采集终端、低压采集终端和无线数字公用通信网的数据进行协同分析，利用逻辑识别的手段检测窃电行为。这种算法在常德供电局已得到应用，运行效果良好[4]。

大数据分析是利用数据挖掘的相关算法，结合高性能计算机，实现对窃电行为的分析。王颖琛等利用高维随机矩阵分析了窃电行为的统计学特性，结合大数据分析技术，实现了对窃电发生的判别、窃电发生时间的确定、窃电地点的精确定位、窃电类型的判别[5]。杨兰从配电环节进行数据分析，计算配电管理系统的特征值，从而对窃电行为进行判断[6]。

由于巨大的商业利益，窃电领域是一个攻防升级都很快的领域。近几年，窃电方式已经由传统的线路改造、电表改造，升级为智能电表改造[7]、虚假数据智能生成[8-9]等多种手段并举的模式。本文所涉及的光伏领域，窃电案例逐年上涨，如果不尽快遏制此类行为的苗头，将给电网行业、社会信用带来难以挽回的损失。

聚类算法是机器学习中涉及对数据进行分组的一种算法，它包含HCM（硬聚类）和FCM（模糊聚类）。HCM 是在基于目标函数的聚类算法理论中比较完善、应用比较广泛的方法。本文从智能电表入手，以电力数据特征量的提取为手段，综合HCM 算法，对窃电行为进行挖掘，并通过在浙江台州的实际应用，对该算法的有效性进行验证。

1 光伏出力曲线特征的提取

1.1 出力曲线实例分析

近几年，随着通信技术的提升，智能电表分时段甚至按照小时、分钟上送发电量已经常态化。细粒度的数据为精细化的窃电检测提供了可行的基础，同时也抬高了样本的维数。考虑到用户数量多，这种高维数为数据的分析带来了极大的挑战。提取数据中的关键信息，降低数据复杂度，是解决此类问题的必由之路。

图1 以浙江台州某光伏用户的数据为依据，选取了1 条冬季出力曲线和1 条夏季出力曲线。2条曲线均每隔0.5 h 进行一次采样，每天有48 个量测点。夏季光照强，随着太阳的升起，光伏曲线迅速抬升，达到出力极限以后将不再飙升；冬季光照弱，上午的出力曲线缓慢爬升，但始终不会达到出力极限。

江浙地区极少出现雷暴天气[10]，因此光照强度很少出现断崖式的变化，而是呈现出比较平缓的出力曲线[11]，这一特点为本文的特征提取算法提供了依据。

图1 不同季节的出力曲线示例

1.2 出力曲线的理论分析

光伏出力遵从下述公式[12]：

式中：Ep为单位时间段内光伏设备的发电量（以图1 为例，Ep为0.5 h 内光伏设备的发电量），由于时间间隔已经取定，因此该量与发电功率是一致的；H 为单位时间、单位面积内的太阳能总辐射量（以图1 为例，就是0.5 h 内单位面积的太阳能总辐射量）；S 为光伏设备中光伏组件的总有效面积；K1为光伏组件的转化效率；K2为系统综合发电效率。从式（1）可以看出，对于一台固定的设备，光伏设备的发电量只与光照强度有关。

在江浙地区，由于雷暴天气极少，一天内的光照强度变化一般比较平缓。因此，可以直接根据几何学原理，推算太阳光照强度的变化规律[13]，即：式中：Hmax，Hmin分别为一天中光照强度的最大、最小值；t0为日出时间，在式（2）中换算为相对于0:00 的秒数；t3为日落时间。公式的第一行是针对白天的，近似服从三角函数的变化（越接近太阳直射点的地区，上述公式越准）；公式的第二行是针对黑夜的，圆月和残月光照强度不同，黄昏后、深夜和黎明前光照强度也不同，但是这样强度的光照远远弱于白天，因此其波动可以忽略不计。

综合式（1）和式（2）可以看出，当白天光照强度没有达到光伏设备的发电极限时，发电功率近似可以看成三角函数。

为了得到更精确的结果，可以采用小波分析和支持向量机相结合的方式进行预测。

采用拉格朗日乘子，可以得出回归函数：

式中：k（xi·x）为核函数；为径向参数；c 为偏移。

式中：g 为核函数参数；z 为中心点。通过设置g和z，可以得到更好的预测函数。

1.3 出力曲线的特征提取

本文中曲线采用一天48 点的采样频率，在实际应用中，一天24 点、96 点、288 点的采样频率均可能出现。如果直接对量测点的数据进行分析，则难以评估曲线自身的变化规律，而且会造成“维数灾难”。不同用户的采样频率可能不同，这也为统一分析带来了挑战。因此，需要从曲线中提取特征量，作为后续分析的依据，以此降低数据维数。从上一节的分析可知，光伏设备的出力曲线包含时间、负荷量两大特征。

首先看光伏出力曲线时间特征，如图2 所示。从图中可以看出，时间特征主要包含t0，t1，t2，t3这4 个时间点。图中每0.5 h 进行一次量测，因此横轴数字1 代表0:00，2 代表0:30，3 代表1:00，依此类推。t0，t1，t2，t3分别代表出力开始明显上升、出力临近极大值、出力开始下降、出力临近极小值的时间点。可以对比前一个采样点和后一个采样点的出力：两者之差的绝对值第一次大于阈值时，意味着t0或t2的到来；两者之差的绝对值第一次小于阈值时，意味着t1或t3的到来。事实上，上述4 个时间点中，t0，t3与日出、日落时间是吻合的，可依据所在地区的纬度和日期得出。光伏设备如果没有达到出力极限，那么t1，t2取值相同，对应日中时间。由于江浙地区的实际时间略晚于东八区时间，所以t1，t2取值相同时（未达到发电极限），该时间略晚于正午12:00。

负荷量的特征值主要包括最高单位时间发电量EPmax和最低发电量EPmin。此外，需要拟合发电功率爬升和下降阶段的曲线。

图2 光伏出力曲线的时间特征

一般而言，拟合曲线可以借助小波分析、傅里叶级数[14]等多种方式。从前文分析已经可以看出，爬升和下降阶段的曲线与三角函数相似。因此，可以将曲线拟合为三角函数，拟合算法借助了傅里叶级数的思想。假设t0，t1，t2，t3这4 个时间点对应的量测点编号为m0，m1，m2，m3则：

式中：a 为起始点参量，相当于傅里叶级数中的第一个参量（平均值）；b 为一阶分量，与曲线的波动程度直接相关；Ep（n）为n 点的单位时间发电量；Eps（n）为利用公式拟合后所得到出力曲线。

光伏出力曲线的拟合如图3 所示。可以看出，利用式（5）所得到的拟合效果与曲线的实际情况相似度较高。从经验来看：如果被拟合日的气象状况比较平稳，那么拟合效果较好；如果被拟合日的气象发生剧烈变化，那么拟合效果较差。但是，如果用户存在光伏窃电行为，这种拟合效果极可能不佳，除非窃电人刻意以三角函数为模型上送发电数据。

图3 光伏出力曲线的拟合

基于上述分析，需要对曲线特征进行进一步的提取。式中：am0，m1（Ep）为上升段的平均负荷；σ 为上升段实际曲线与拟合曲线的方差，通过σ 表示2 条曲线的差异。

此外，判断是否发生窃电的另一个重要方法是：后一日发电量是否相对于之前一段时间（如30天）发生了显著的增加。

式中：Ep，sum，分别为光伏用户当日发电量、此前30 天的日均发电量；δ 为两日发电差值率。如果δ 数值较大，那么说明天气由阴转晴、用户扩容，或者存在窃电行为。

至此，可以得出光伏发电曲线的特征向量：

式中：σ′为下降段曲线方差；Epmax，Epmin分别为Ep的最大值、最小值。

可以看出：式（8）的前7 个分量与负荷的具体数值无关，表示了负荷曲线本身的形状，后文的聚类算法主要针对前7 个分量进行分析；后3个分量涉及了具体的数值，与辅助的反窃电手段有关。

2 HCM 聚类算法

2.1 算法选择依据

HCM是一种典型的动态聚类算法，它有明确的聚类中心，自适应能力差，但速度快；FCM 则没有固定的聚类中心，需要反复迭代，但自适应能力很强[15-17]。

在光伏发电的场景下，正常的发电用户一般遵循较为相似的发电曲线[18-23]。如果能够找到正常出力的代表曲线，以之作为聚类中心；针对常见的窃电手段，考察其曲线特征，以之作为窃电用户的聚类中心。基于这种相对固定的特征，可以采用HCM 算法进行聚类。

2.2 聚类中心

聚类中心主要有三类：第一类是正常的出力曲线，第二类是设备故障等原因导致的非正常出力曲线，第三类是存在窃电行为的出力曲线。

对于正常的光伏出力曲线，首先指定可信的光伏设备。这种设备一般为电力公司自营的光伏设备，可以随时检查其运行工况。如果反窃电系统未接入此类设备，那么以日期、气象因数进行计算，或者任取若干个设备的出力曲线，计算其平均特征。如果系统内包含多种差异较大的光伏设备，可以从每一类设备中选取一个聚类中心。

对于故障场景，其突出特点往往是比正常的光伏设备发电量要低很多。其特征为b，Epmax的数值都很小，发生故障的第一天，δ 的数值为负数。基于该特征，可以设置两个故障设备的聚类中心，具体如下：

对于存在窃电行为的光伏设备，其窃电目标是在数据上造假，比正常的数值要大很多。基于这种考虑，其特征包括：

（1）t0，t1，t2，t3这4 个时间点明显与昼夜变化背离。作为上升段，发生窃电时t0，t1一般比正常值偏早；作为下降段，发生窃电时t2，t3一般比正常值偏晚。

（2）σ 和σ′数值很大，少数情况下很小。如果正常设备的σ 较大（当日天气发生了剧烈变化），那么对于σ 较小的设备，可以判定存在窃电行为。

（3）δ 数值远高于正常值，可以判定存在窃电行为，但该判据仅适用于窃电行为发生的第一天。

（4）Epmax，Epmin数值很大，这是一部分以改接线路为手段进行窃电的典型特征。

对于窃电行为，可能同时存在上述4 种特征，也可能只存在部分特征。以正常曲线的聚类特征为基础，将式（8）中的特征分量分别替换为异常值，就可以得到多个代表窃电的聚类中心。窃电情形的聚类中心与正常情形的聚类中心越近，窃电检查的标准就越严格，误报率越高；反之，则漏报率越高。

2.3 聚类算法

对于所有样本，正常光伏用户的集合记为Sn，疑似窃电的用户集合记为Ss，故障设备的集合记为Sb，未确定是否正常的用户集合记为Su。在确认聚类中心以后，按照下述步骤挖掘窃电行为。

（1）考察样本的δ。设关于δ 的阈值为δΔ，聚类中心记为cΔ。如果δ-cΔ≥δΔ，表明当日发电量远高于前一日，而且超出了正常范围，该样本疑似窃电，记入Ss；如果cΔ-δ≥δΔ，表明与其他用户相比，当日发电量低于前一日，该样本疑似故障，计入Sb；剩余样本计入Su继续考察。

（2）对保留在Su中的样本考察t0和t1，定义基于这两个量的距离如下：

以上述距离为基础，利用正常、故障、窃电3 个聚类中心以及HCM 算法，将样本分别计入Su，Sb和Ss。

针对t2和t3也进行相似的分析。

（3）对保留在Su中的样本考察σ。定义σ 的距离如下：

在σ 值较小时，只在比σ 数值大的一侧设置窃电的聚类中心；在σ 值比较大时，在该值的两侧都设置窃电的聚类中心。由于方差的影响因素较多，在实际应用中可以适当调整聚类中心的远近，以决定方差的敏感性。应用HCM 算法，将样本分别计入Su和Ss。

对于σ′进行类似操作。

（4）对Epmax，Epmin，b 设置阈值，当参数值高于阈值时，则判定存在窃电行为。Su中剩余的样本判断为正常光伏设备，计入Sn。

至此，完成样本分析。对于Ss中的光伏发电设备，可以通过上门检查等手段进行最终确认。

3 辅助发现手段

通过负荷曲线以外的数据信息，可以辅助发现窃电行为，进一步提高判定的准确率。这主要包括以下几种信息：

（1）电压、电流量。功率是通过电压和电流计算出来的，如果改变接线，那么可以提高电压或电流值，以此获得更高的光伏补贴。设置电压、电流的阈值，当这些数据处于非正常水平时，可以判定窃电。

（2）用户接入信息。用户接入当日，不进行第2 节的分析，这样避免窃电误报。次日开始，接入线路本身的功率阈值，可以作为判断Epmax是否超过限值的依据。

（3）光伏设备的信息。包括光伏设备采用单晶硅还是多晶硅、峰值功率、峰值电压等，以及逆变器的额定输出容量、整机效率。将这些信息与Epmax对比，可以作为判定窃电的依据。

（4）气象信息。主要指光照变化情况，利用光照信息可以拟合出正常样本的σ，σ′，δ 等数据。如果与光伏设备本身的参数结合，这个数据会更准确。

（5）线路信息。结合配电网数据，考察各线路的线损，线损较高的线路存在窃电的可能性较大。针对该线路，将相应的窃电情形的聚类中心参数设置为与正常情形较接近的值，提高窃电检测的严格度。

4 应用实例

自2019 年3 月以来，在浙江台州选择了948家光伏用户进行了实例验证。这948 家光伏用户以居民为主，峰值功率都在3～5 kW，用户申请光伏并网时注册的设备信息均可查。光伏设备的数据每隔0.5 h 上送至电力公司，数据只包含每小时的发电量。由于不能接入气象、配电网数据，因此在查找窃电行为时未考虑此类数据。

这些用户中包括台州一家电力营业厅，其峰值发电功率为5 kW。基于其可信任性，将该用户设置为正常用户的聚类中心。以5 月8 日首次运行本文算法为例：t0为330（以分钟为单位，对应5:30），t1为510，t2为960，t3为1 140；σ，σ′分别为0.253 9 和0.103 9；δ 为0.012；Epmax为245 kWh，Epmin为1 kWh，b 为240.48 kWh。对于故障场景的聚类中心，将Epmax设置为0；对于窃电场景的聚类中心，将δ 设置为1。

经过本文算法分析，得到如图4 所示的散点图。