面向IEEE COMTRADE格式的海量录波数据并行压缩／解压算法

桂 勋

（电子科技大学 航空航天学院，四川 成都 610054）

0 引言

电力系统电磁暂态过程是电力系统中短暂的、但非常重要的物理过程。所有的故障都伴随着相应的电磁暂态过程，现代电力系统故障诊断、广域保护和稳定分析愈加依赖于电网动态监测系统采集的含丰富故障信息的暂态过程信号。因此，对电磁暂态的监测、记录和分析一直是电力系统中一个重要的研究方向。近年来随着电子技术的快速发展，电力系统内暂态录波逐步向高采样率、连续稳态记录和海量存储方向发展[1]。为了提高故障录波数据的传输效率和降低存储空间，众多文献提出了各种压缩算法[2-15]，其中可分为直接离散小波变换压缩[2-7]、直接熵编码[8-9]压缩以及这 2 种算法的综合[10-13]共 3 类。其中在离散小波压缩研究领域，由于第1代离散小波变换Mallat算法[14]相对于第2代基于提升格式的离散小波变换算法[15]无法实现无损压缩，且计算量大、需要额外内存[10-12，15]，第 2 代基于提升格式的离散小波变换和熵编码的综合压缩算法[10-12]成为当前此领域内的研究热点。

然而到目前为止这些压缩算法在实际应用中却存在以下问题。

a.以往压缩算法是为了解决电话线传输速率过低的问题，但随着电力通信网络的改善，百兆、千兆以太网大规模安装，通信效率已经不再是难以克服的瓶颈，因此，目前众多厂家已经放弃了原有的复杂压缩、解压方案，而直接在高速以太网上传输录波文件。而另一方面，保信系统厂家不可能让负责上传录波文件的子站和主站软件系统为每个录波器厂家实现复杂的压缩／解压通信协议。正因为这两方面的原因，限制了原有压缩算法的实际应用范围。

b.以往压缩算法是为了解决海量录波数据的存储问题。但在电力系统内，目前所有的暂态记录格式逐渐被统一到了 IEEE COMTRADE 1999[16]（下文简称COMTRADE）标准上，而几乎所有的第三方暂态分析软件都是以COMTRADE标准为输入格式的。因此，即使在保信系统上采用了压缩／解压通信协议，在主站上还是要解压后以COMTRADE标准保存，否则第三方分析软件无法使用。因此压缩算法也没有解决实际中面临的海量COMTRADE数据存储问题。

c.随着多核和嵌入式多核的大规模普及，带来了一场计算技术革命：传统的串行编程技术被完全颠覆，逐渐被基于多线程模式的并行编程技术所取代。软件要提高运行速度和效率已经不能依赖CPU主频，而要使软件系统本身在设计上适应新的计算机体系结构。而此领域内以往提出的算法都是串行算法，已经无法适应计算机体系结构的发展。

笔者经过现场调研后认为，只有直接面向COMTRADE标准文件的压缩算法才能彻底解决问题a、b，假如实现了通用的COMTRADE标准压缩／解压算法及其工具，就可利用所有保信系统内都具有的文件传输协议来直接传输，在主站内就可免去解压，直接进行存储，需要分析的时候，只需调用相应的解压软件。这种应用模式可把以往算法和传输的紧密关系完全解耦，更加简单可靠，而各种第三方分析软件系统也可携带此工具，为各种电力用户节约存储空间。基于这种思路，笔者在已有研究的基础上，结合多年的电力系统故障分析软件研发经验，提出了面向COMTRADE标准的并行压缩／解压算法，并构建了相应的压缩／解压工具。算法根据COMTRADE中不同种类记录信息的特点，采用不同的压缩算法或混合压缩算法，获得极大的压缩比，且具有线性加速比，可较好地解决以往算法在实际应用中出现的问题。

1 IEEE COMTRADE格式标准

IEEE Std C37.111—1999规定和记录信息相关的文件有4个：头标文件（文件扩展名为HDR）、配置文件（文件扩展名为CFG）、数据文件（文件扩展名为DAT）以及1个可选的信息文件（文件扩展名为INF）。其中对于第三方暂态数据分析软件而言最重要的是配置文件和数据文件，其他文件是可选的。对于压缩算法而言，除了数据文件外，其他文件都可视为纯文本文件，以文本流的方式进行压缩。

1.1 配置文件

要求计算机程序能够读取配置文件，获取内部的属性信息，并采用这些属性信息正确地读取数据文件（*.DAT）中的记录数据。配置文件为ASCII文本文件，包含着计算机程序为了正确解读数据文件而需要的信息。配置文件由若干行组成，每个字段都有对应的明确物理意义。其具体内容为：站名，记录装置的特征，COMTRADE标准的修改年份；模拟量和状态量通道的数量；模拟量通道采样信息（包括模拟通道索引号、通道识别名称等配置信息）；状态量通道采样信息（包括状态通道索引号、通道识别名称等配置信息）；被采样线路的频率；采样频率个数，实际采样频率及其对应的结束采样下标；第1个数据点的日期和时间；触发点的日期和时间；数据文件类型，表示数据文件是ASCII文本文件还是二进制文件；时间标记倍乘系数。

1.2 数据文件

数据文件记录着每个采样通道中的每个采样数值。数据文件可以是ASCII或二进制格式。二进制数据文件和ASCII数据文件格式类似，每行应分为T+2列，其中T是配置文件中模拟量通道和状态量通道的总和，另外2列是采样序号和时间标记。二进制文件以BIT形式集中存放状态量通道数据。各个数据之间没有分隔符，每组采样值之间没有回车换行符隔开。各列的具体内容如下：第1列为采样序号；第2列为采样数据的时间标记；第3组的列为表示模拟量通道信息的采样数据值；第4组的列为表示状态量通道信息的采样数据值。其数据文件格式的例子如图1所示。

图1 数据文件格式Fig.1 Format of data file

2 IEEE COMTRADE数据文件压缩 ／解压算法

从图1可见数据文件分为采样序号、时间标记数据、模拟量通道数据、状态量通道数据4个部分。其中采样序号在压缩时可以完全丢弃，下面分别讨论其他3个部分的压缩算法。

2.1 时间标记数据压缩算法

由于采样时间标记是一个递增的量，可通过配置文件中的采样频率间接计算得到，因此压缩此项数据时只需记录数据文件中记录的第1个时间标记数据、配置文件中的采样率个数和实际采样频率及其对应的结束采样下标、配置文件中的时间标记乘数3种信息，即可在解压时进行无损恢复。解压算法公式为：

其中，ta为累加偏移时间，初值为0；f为当前采样频率；n为当前采样频率下从0开始的序号；t1为数据文件中的第1个采样时间；tm为配置文件中的时间标记乘数；ts为采样时间标记。算法首先计算累加偏移时间，并将采样间隔时间单位转换为COMTRADE标准规定的微秒。将累加偏移时间除以时间标记乘数后加上t1，再取整就可无损恢复出所有的时间标记。

2.2 状态量数据压缩算法

在COMTRADE标准中，BIANRY格式数据文件中的状态量数据是以BIT形式保存的，而在ASCII码中则是以1个字符（“0”或“1”）来表示。在BINARY格式中，状态量通道所占存储空间的计算公式为：

其中，mod为取整运算，%为取余数运算。

由式（2）可以知道，即使是只有1个状态量通道，状态量数据也要占据2 Byte的存储空间，即对于COMTRADE标准而言状态量数据保存是有冗余量的。并且通过电力系统内大量COMTRADE标准历史数据的验证，在所有状态量通道中，95%以上的通道记录根本没有状态变化，而只有极少量的状态量通道记录了开关量变位信息。

对于这种只有2种变化的量，结合COMTRADE状态量的特点，算法设计了最简游程编码RLE（Run Length Encoding），即采用游程L的正表示状态1，L的负表示状态0，而L的绝对值表示游程长度，例如对于 111111110000000，就可用（+8，-7）来表示，而对于没有状态量变换的状态量通道，将其分别归为全1状态量通道和全0状态量通道类，其中每一类只是记录对应的通道编号，通道编号的整型字节长度依据通道数量的多少而定，通常情况下数字量通道数都在256以内，所以可以将通道编号整型长度的初值设为1Byte。通过这种压缩编码就可对具有大量状态量信息的COMTRADE标准记录文件进行大幅压缩。

2.3 高频采样模拟量数据压缩算法

高频采样模拟量数据在IEEE COMTRADE数据文件中是以2 Byte的整型数保存的。因此，基于浮点数运算的Mallat算法将无法实现无损压缩，且计算量也比第2代提升格式的离散小波变换至少高出一倍。因此本算法中针对模拟量数据的压缩首先采用提升格式的离散小波变换，根据用户选择是否进行有损压缩，进一步采用硬阈值方法进行量化[7]，最后采用熵编码算法进行压缩的混合压缩算法。

2.3.1 提升格式小波变换

提升格式小波变换与Mallat算法相比的区别在于不依赖傅里叶变换，在时域内就可实现小波构造。对信号进行一次提升变换就相当于进行一次小波分解，一个典型的提升格式离散小波变换包括分裂、预测和更新3个步骤。而重构则是此过程的逆，即更新、预测和合并。整个提升小波的分解和重构如图2所示。

图2 小波提升方法的分解与重构Fig.2 Decomposition and reconstruction of lifting wavelet transform

从图2可知，提升格式小波变换就是不断对偶数序列（把原始序列也看作偶数序列）减半后进行分解的过程，而重构过程则是分解过程的逆。对应于计算内部实现，对其分解和重构的原理如图3所示。

图3 提升方法的计算机实现Fig.3 Computer implementation of lifting wavelet transform

从图3可知，小波分解的第1步分裂，就是将原始序列Sj的偶数部分移动到前半部分，奇数部分移动到后半部分，而后进行预测和更新，一直持续该过程，而重构是此过程的逆，不需要任何格外的内存就可完成全部变换。对于一个有2H个点的输入原始数据，经过h（h＜H）次小波变换后，尺度系数的点为序列左边的2H-h个点，其余点为小波系数，由此可见小波分解层数越深，可供量化的小波系数也就越多，故压缩效率也就更高。本文算法中提供了多种整数小波变换，可供用户进行选择，它们在尺度j上的变换公式如下所示。

其中，Sj、Dj为经过j-1层分解后得到的偶序列和奇序列，l为序列数组下标。

2.3.2 熵编码压缩算法

熵编码压缩主要是利用数据或数据序列出现概率的分布特性来寻求概率与码字长度间的最优匹配。 常用的熵 编码有 Huffman、LZ77、LZ78、LZW、RLE、算术编码等[17]。众多文献已经证明了单纯采用一种熵编码很难取得高压缩比，文献[12]采用了RLE和LZW的混合编码，取得了较好的压缩效果。本文中采用了著名的混合熵编码Deflate算法，其采用了经过优化后的LZ77和Huffman的混合熵编码算法，压缩比和速度更快。由于Deflate算法很好地平衡了压缩效率和速度，所以很多压缩软件（PKZIP、WinZip、gzip）均采用了Deflate算法或由其衍生出来的其他算法，Deflate算法的解压算法名称为Inflate。

LZ77编码算法是由Lempel和Ziv于1977年提出的一种字典压缩算法，其核心思想是把已输入的数据流的一部分作为字典，编码器为输入流开一个滑动窗口，随着字符串编码的输入把窗口中的数据从右向左移动，并在其中寻找数据中相同的部分，即寻找最长字符串匹配，由此在编码过程中形成一种内容形式为（距离，匹配长度）的动态字典。LZ77在压缩时需要进行大量的字符串匹配工作，这也是LZ77效率上最低的地方，为此Deflate算法采用哈希表来提高匹配速度。在哈希表匹配过程中，采用当前字符串中头3个字符来计算1个哈希值，并将具有相同哈希值的匹配字符串利用链表连接起来，而寻找最长匹配字符串的过程，即通过哈希值找到链表头后，遍历链表，最后确认最长匹配字符串。

在经过LZ77编码输出后，Deflate算法采用Huffman编码进一步压缩数据。Deflate算法在对一块数据进行Huffman编码前，会同时建立静态Huffman树和动态Huffman树；然后根据要输出的内容和生成的Huffman树，计算静态和动态Huffman树的最终生成块大小，而后进行比较，采用生成块较小的算法进行编码。采用这种方法，Deflate算法避免了由于文件较小时动态Huffman树编码比静态Huffman编码生成块大的情况，保证最终输出块始终小于编码输入块，这一点是其他算法很难做到的。

2.4 低频采样模拟量数据压缩算法

在电力系统内部，很多COMTRADE标准记录文件中的最后一个采样段都是低频采样的（小于100Hz）模拟量数据，这些低频数据通常是当前模拟量通道的有效值。这些数据对于故障分析没有太大意义，并且通常也是没有任何变化的量，对于这类低频数据本算法缺省采用RLE算法，将数据压缩成（数值，长度）对。

3 最佳小波分解尺度下的最佳数据区间划分算法

由于COMTRADE标准记录文件内通常包含多个采样频率下的数据，压缩算法不可能一次性加载全部模拟量通道数据，并且待压缩数据点数也不太可能是2的幂级数，因此可能需要在采样频段尾部补零，压缩算法就面临一次读取多少采样数据进行小波分解、压缩的问题。为此本文提出了一种采用配置文件中提供的采样频率信息的压缩数据区间最佳划分算法。

为防止小波分解层数过多而导致基频分量在小波变换过程中被分解，影响压缩效率，本文算法采用式（7）[7]计算最佳小波分解层数，并由此倒推出要获取最佳压缩效果的最小数据点数公式见式（8）：

式（8）的意义是在采样频率fs下，经过n层小波分解，最后刚好有2个点为尺度系数，即实际中要求的待分析数据量必须为 B 的 2k（k=0，1，2，…）倍，才能获取最佳压缩效果。根据式（8）可得到在采样频率为fs、采样数据总数为K时的最小补零个数为：

假设处理数据缓冲区大小为SB=2m（m＞n+1），于是可得以SB划分的数据区间个数M为：

而对余下的数据则需要计算出其相对于B的倍数L。

对式（11）中的L以二进制方式从最高BIT位起依次判断P位是0还是1，如果是1则利用式（12）获取一个数据划分大小Y：

通过以上算法就可得到一组可最少补零的数据划分，例如在fs=9600 Hz时，共采集了15543点数据，当前计算缓冲区为4096个整型数据，则采用以上算法可自动将数据划分为3个大小为4096的区间、1个大小为2048的区间、1个大小为1024的区间和1个大小为256的区间。最小的一个区间需要补零73个。由此可见采用此算法所划分的区间，满足了最佳小波分解层数的要求，且补零个数最少，在压缩效果和压缩效率间找到了一个较好的平衡点。

4 并行压缩算法

如图4所示，并行压缩算法首先解析COMTRADE配置文件，再利用其内信息创建压缩文件头并写入压缩文件后，将COMTRADE配套文本文件（*.cfg等）依次采用Deflate算法进行压缩，并写入压缩文件中。随后通过第3节提供的划分算法，依据采样率信息自动创建数据压缩划分区间，并根据此区间依次将原始数据读入并解析到模拟量通道数据压缩共享管理器内预先分配好的通道数据缓冲区中。此时通知激活已经创建并等待的模拟量压缩线程组和状态量压缩线程组（并行压缩算法采用的线程CPU核心映射原则是：每个核心映射1个模拟量压缩线程和1个状态量压缩线程）。而后这些线程将向模拟量通道数据压缩管理器请求待压缩通道序号。

对于模拟量数据而言，压缩线程通过判断采样频率是否高于100Hz选用不同的压缩策略。对于有损压缩，当采用提升格式小波变换后，采用文献[7]中介绍的阈值处理方法，该方法在阈值选择过程中考虑了噪声的小波变换系数在小波系数空间的传播特性，能够很好地抑制噪声对压缩性能的影响。压缩后的数据直接写入压缩文件内。

图4 COMTRADE并行压缩算法Fig.4 COMTRADE parallel compression algorithm

其中，tj为在尺度j下的硬阈值；σ为噪声强度；N为尺度j下的小波系数个数；γ为常数，选择范围为1～3。其中噪声强度σ可根据第i次小波变换后的小波系数计算得到，计算公式为：

其中，di为第i次小波变换后的小波系数，为均值，N为小波系数个数。本算法中，γ取值1、2、3分别被映射为“一般”、“较好”、“最好”3种压缩选项。

对于状态量数据而言，压缩线程每次更新在状态量通道数据缓冲区的压缩数据，只有当压缩线程遍历状态量通道的全部状态量数据后，才提炼得出全1通道序号数组、全0通道序号数组和状态量发生变化通道压缩数据，最后将这3类压缩数据写入文件。

这种基于申请序号的线程竞争计算方式，可使压缩算法效率随着CPU核心数量的增加而提高。

5 压缩文件格式

经过第4节的并行压缩算法，最终压缩文件形成了如图5所示的压缩存储格式，其中每个部分都是一个可根据压缩数据变化的量，不能直接采用固定大小结构体设计，为此设计了如下所示的通用变长结构体。

此结构的大小计算程序公式为：

其中，Count为实际存储数据大小，采用此大小分配内存后，对于type类型数据Data，可直接用数组形式访问，这样就形成了一个变长的结构体。

图5 压缩文件格式Fig.5 Format of compressed file

压缩文件的读取就是利用nStructSize来确定读写大小，采用nStructType来判断数据类型，进而准确读取所有的压缩数据块。在此压缩存储格式中，模拟量数据以压缩时数据划分的形式顺序保存，而状态量数据则集中保存在文件最后，在压缩文件头中有指向状态量压缩数据的具体偏移数据。

6 并行解压算法

图6 并行解压算法Fig6 Parallel decompression algorithm

对于第5节阐述的压缩文件格式，采用如图6所示的并行解压算法进行解压。在解压主线程内，首先采用Inflate算法对COMTRADE配套文本文件压缩数据进行解压，并保存到文件。之后解压主线程通过压缩文件头内提供的状态量压缩数据偏移量信息，直接定位并一次性将压缩数据读入到状态量压缩数据共享存储器内。随后依据压缩文件头内提供的数据划分信息，将一次数据划分的全部数据读入模拟量通道压缩数据共享管理器内，之后解压主线程将激活模拟量通道数据解压线程，并等待解压完毕（并行压缩算法采用的线程CPU核心映射原则是每个核心映射1个模拟量解压线程）。激活后的各个解压线程向共享管理器申请解压通道序号，根据原始采样率判断是否为高频模拟量压缩数据，之后采用不同的解压策略。当一次数据划分并行解压完毕后，解压主线程等待结束，随后采用状态量压缩数据解压出当前数据划分下的状态量信息，再将时间信息还原后加上采样编号，最后写入解压文件。持续此过程直到所有数据划分区间被读取并被并行解压完毕。同第4节所述并行压缩算法类似，此并行解压算法采用基于申请序号的线程竞争处理方式，可使解压效率随着CPU核心数量的增加而提高。

7 试验

将本文算法应用于4核CPU计算机，配置为：单核主频 2.4GHz，内存 1GByte，硬盘转速 7 200 r／min。

从表1可见Deflate算法的速度很快，而其解压算法Inflate的速度是Deflate算法是5～9倍。从表2可见，对于同样大小的原始计算数据，不同插值点的提升格式小波变换计算时间的差别并不明显，以6.4×104个计算点为例，CDF（2，2）和 CDF（8，2）计算时间相差208 ms。

表1 Deflate／Inflate算法效率试验结果Tab.1 Results of Deflate／Inflate algorithm efficiency test

表2 提升格式小波变换效率试验结果Tab.2 Results of lifting wavelet transform efficiency test

为了体现算法在实际应用中的压缩效果，选用了在电力系统内部具有典型意义的4种COMTRADE文件（文件详细信息如表3所示），其具体分类如下所表述。

文件1：典型的单相接地故障记录数据文件，记录末尾有少量低频采样数据。

文件2：典型的保信系统上传文件，记录大量意义不大的稳态数据及其开关量变位信息，记录末尾有较多低频采样数据。

文件3：故障持续时间较长，采样频率高，记录中没有低频数据。

文件4：系统发生较长时间扰动，录波器连续进行长时间录波，完整记录整个过程，此类连续长时间采样也是未来电力系统自动化录波的发展趋势。

为和以往文献试验结果对比，试验采用以往众多文献推荐的 CDF（4，2）小波，从表 4、5 所示的压缩比试验结果上看，本文算法均可获得非常大的压缩比，当COMTRADE文件逐渐变大后，压缩比也逐渐变大，特别是文件4，在有损压缩时取得了高达96.6和287.0的压缩比。而对于无损压缩，此算法也取得了很好的效果，文件4取得了5.7和16.9的高无损压缩比，而文件2由于尾部有较多低频数据，取得了更高的压缩比。

表3 压缩／解压试验COMTRADE数据Tab.3 COMTRADE data of compression／decompression test

表4 有损压缩试验结果Tab.4 Results of loss compression test

表5 无损压缩试验结果Tab.5 Results of lossless compression test

表6中针对文件4类型的COMTRADE文件，采用4种小波分别进行有损压缩，可见随着参与插值的点数增加，获得了更小的压缩文件。笔者还通过大量试验证明，在诸多文献中提到的CDF（4，2）的压缩效果最好的结论，实际上过于片面，试验中CDF（6，2）、CDF（8，2）小波对于COMTRADE中大量存在的小扰动系统录波数据具有更高的压缩比。并且由表2可见插值点数的增加对于提升小波变换速度的影响有限，所以笔者认为对于实际中的COMTRADE文件而言，推荐采用 CDF（6，2）和 CDF（8，2）这类插值点更多的小波，使得变换后得到的小波系数更加平滑，有利于压缩。

表6 不同小波选择的压缩结果Tab.6 Results of compression by different wavelets

表7显示了对4种不同文件的内故障记录通道有损压缩后和原始通道数据计算均方误差，采用的小波为CDF（8，2）。从表中可见有损压缩后的解压数据和原始数据误差都在小数点后3位，满足有损压缩后的可用要求。

表7 故障通道数据有损压缩最大均方误差Tab.7 Maximum deviation of loss compression for different faulty channels

表8中采用文件4类型COMTRADE数据在2核、3核、4核计算机上进行了加速比试验。可见，此算法可随着未来COMTRADE数据文件的加大、CPU核心数的增加，获取到线性加速比。

表8 并行压缩／解压加速比试验结果Tab.8 Results of parallel compression／decompression acceleration ratio test

8 结语

以往压缩算法文献都只是针对模拟量数据进行讨论的，而完全忽视了电力系统内绝大部分录波数据都是以IEEE COMTRADE格式进行保存的情况，并且由于电力系统通信网的快速发展，原有算法的应用背景也发生较大变化，原有研究已经不能适应电力系统的发展，因此本文提出了直接面向IEEE COMTRADE格式的并行压缩／解压算法，并构建了相应的通用压缩／解压工具，大量试验证明对于绝大部分现场运行，COMTRADE文件都可获得很大的压缩比，并且算法适应了多核并行计算的发展趋势，具备线性加速比。