引言

国家地震烈度速报与预警工程项目（简称“预警项目”）是一项具有防震减灾实效的民生工程和社会工程，建成后将有效提高减灾和社会服务能力。河南省预警项目中一般站烈度计设备于2021年上半年安装完毕，相关记录数据可通过网络远程访问下载，为初步计算一般站的背景噪声创造了条件。预警项目针对台站数据包格式做出了新定义，各厂家生产的硬件采集设备已按照新的格式规程生产数据，而现有计算软件^①附带的数据分道软件无法直接打开预警数据采集器产出的数据，因此，本文通过对预警项目台站数据包格式规程的解读，采用C语言开发了预警项目台站波形数据解压软件，产出相应的分道文本数据文件，并通过具体实例对预警项目连续波形数据包格式和相应的程序实现给出详细说明。

1.   预警项目台站波形数据包组成

1.1   波形数据包格式

预警项目台站数据包格式规程中定义了多种数据包格式，其标识有“cc”、“si”、“ti”、“wc”、“wt”和“ws”，本文着重研究“wc”数据包，即连续波形数据包。“wc”数据包以miniSEED数据卷为基础，做了必要的修订，主要变化为：① miniSEED数据卷固定头段区的第0至第1个字节定义为“wc”数据包标识位置；② miniSEED数据卷固定头段区的第2个字节内高5位bit到第5个字节定义为数据包顺序编号位置，用于数据传输时断点续传，而第2个字节内低3位bit定义为数据包长度索引；③ miniSEED数据卷中“纯数据子块 $[1000] $ ”后增加了8个字节的数据标识区，用于标识通道信息顺序、数据量纲和灵敏度信息；④其他少许的新定义，如第36个字节第7 bit定义为入网标记等。由此发现，预警项目对miniSEED数据卷仅做了少量修订。为提取“wc”类连续波形数据具体内容，需对该数据卷格式信息做出解读，在此仅对提取波形使用到的信息做出说明。miniSEED被定义为SEED格式标准中的数据记录格式，一个完整的预警项目miniSEED数据卷包括固定头段区、可变头段区和数据区（中国地震局，2003）。

固定头段区包含使用数据所需的最低限度信息，主要信息按字节从前向后依次为波形类别标识、数据包序号、质量指示符、台站标识、位置标识、通道标识、台网标识、记录起始时间、样本数目、采样率、子块等。前20个字节除第2～5个字节外，均采用ASCII码表示，后28个字节采用二进制表示（张旸等，2007），固定头段区共48个字节。

可变头段区由多个可选的子块构成，在预警项目中选择的子块类型仅为“纯数据子块 $[1000] $ ”。该子块包含了对数据区解码的必要信息，主要信息有编码格式代码，一般为0 B，代表Steim2压缩格式；字序代码一般为1，代表大端字节序，高位字节在前，低位字节在后；标识数据记录长度的指数代码如果为8，则数据记录长度为2⁸字节。纯数据子块 $ [1000] $ 为8个字节，均采用二进制表示。

纯数据子块 $[1000] $ 区域后的8个字节是最具有预警项目特色定义的区域，位于可变头段区内，因此在此区域增加自定义信息最适合。本块信息主要涉及有关数据标识的内容，包括数据区数据的通道顺序，如“ENZ”表示数据按照“**E”通道数据、“**N”通道数据、“**Z”通道数据的顺序排列；数据量纲取值范围为0～3，0表示无量纲，1表示位移，2表示速度，3表示加速度；灵敏度因子记录仪器记录的灵敏度因子取值范围为0～0x3FFFFFFF。预警项目定义的数据标识内容多数为二进制，实际上这8个字节内容可视为1个变相的子块放在可变头段区内，以增加对数据区的说明。

固定头段区和可变头段区数据包含了台网、台站、通道、采样率、样本数目、时间起点、编码格式、字节序、数据记录长度等数据描述信息（王晓磊等，2016），占据了每条miniSEED数据卷的前64个字节区域，这些信息是波形数据区的辅助信息（杨周胜等，2019），也是编程提取和使用的主要参数。

数据区包含实际的时间序列数据（中国地震局，2003），数据区位置从固定头段区第44至第45个字节中内容指定的字节序数处开始，这部分数据块一般采用Steim2算法压缩（何加勇等，2009），由若干数据帧组成，每个数据帧为64个字节。波形数据区具体大小及数据通道、采样率、样本数目、量纲、字节序等信息均已在上述固定头段区和可变头段区中定义。

1.2   Steim2压缩算法

预警项目中的波形数据区数据为Steim2格式，获取观测数据后，首先需获取打开Steim2格式数据压缩的钥匙才能使用。Steim2压缩算法是借用语音压缩中最简单的差分压缩算法发展而成的（王翠芳等，2011），其特点是无损压缩，目前国内外地震数据传输和本地存储普遍采用了压缩算法。地震数据采用无损压缩算法的基本原理为：在正常地震背景情况下，地震计输出的时间序列数据通常具有较强的关联性，即前一个采样点相比后一个采样点变化小，差值小，该差值占用计算机内存小，可能远远小于32 bit，而保存1个原始值却需32 bit，因此若多个差值按规则填充满32 bit，可极大地节省存储空间，这种方式可理解为数据压缩。后期为恢复原始值，须在特定位置保存一个前向积分常数。笔者认为Steim2数据格式在压缩时具有高保真、高效的特征，符合地震动波形记录“长期平静”和“长期无震”的特点（王翠芳等，2011），但此算法在地震事件发生时，压缩效率显著降低，这是无损压缩算法的特点。

根据Steim2格式数据压缩基本思想，合成压缩数据时，首先利用直接观测数据计算差值序列，然后按照规则利用差值序列整合重建32 bit数据序列。获取数据后，用户更关心数据解压过程，其步骤与压缩相反，具体如下：

（1）按规则将32 bit数据序列拆分为差值序列

数据拆分的规则钥匙随着压缩数据进行传输，压缩数据帧的结构如图1所示。

图 1  数据帧结构

Figure 1.  Data frame structure

下载: 全尺寸图片 幻灯片

图1中每个字母长度均为32位，解压的钥匙保存在每行行首W₀和其后W₁～W₁₅中最高的两位中，拆分对象是W₁～W₁₅中后30位二进制数据。钥匙（ C_k）须与锁（ W_k）一一对应，才能顺利解锁，如图2所示。

图 2  C_k与W_k的关系及dnib位置

Figure 2.  Relationship between C_k and W_k and dnib location

下载: 全尺寸图片 幻灯片

C_k为2 bit，取值为00、01、10、11，若压缩数据为Steim1格式，这4类钥匙可拆分W_k数据，但预警项目中规定使用Steim2格式压缩数据，需另外4类钥匙才可打开第2道门，从而拆分W_k数据，而打开第2道门的钥匙隐藏在W_k的最高2位中。将C_k称为“广义字节编码”，W_k的高端2位（dnib）称为“细解码”，其组合使用如下：

C_k=00：同Steim1，W_k包含非数据信息

C_k=01：同Steim1，W_k中包含4个8位差值

C_k=10：查看dnib，即W_k的高端2位

　　dnib=01：W_k中包含1个30位差值

　　dnib=10：W_k中包含2个15位差值

　　dnib=11：W_k中包含3个10位差值

C_k=11：查看dnib，即W_k的高端2位

　　dnib=00：W_k中包含5个6位差值

　　dnib=01：W_k中包含6个5位差值

　　dnib=10：W_k中包含7个4位差值

根据上述规则，按照C_k和dnib组合定义的拆分方法，可将W_k从最低位开始按照4、5、6、8、10、15、30位的方式进行切割，分离出的数据即为差值。

（2）将差值序列还原为原始数据

获取差值后，通过前向积分常数还原数据，前向积分常数保存在数据帧0中的X₀内。依次算出原始数据，计算步骤如下：

X₁=X₀+d₁

X₂=X₁+d₂

X₃=X₂+d₃

…

X_n=X_n−1+d_n

在数据帧0中的X_n内还保存有这组数据的反向积分常数，利用该反向积分常数与最后计算出的样本原值进行比较。如果二者不同，说明数据出现了错误。

2.   数据包解析实例

以预警项目一般站MI3000-C烈度计数据采集器产出数据为例，对连续波形包数据进行解析，给出的实际二进制（用十六进制表示）数据如图3所示。

图 3  MI3000-C头段和数据帧0说明

Figure 3.  MI3000-C Header and Data frame 0 description

下载: 全尺寸图片 幻灯片

从固定头段区和可变头段区可提取诸多参数信息，如可变头段区的“C0”字节低2 bit中包含数字0，说明数据单位为无量纲，此处解压后的数据单位为count。如果“C0”字节低2 bit中包含数字为1或2，则必须给出确定的灵敏度因子，可直接输出物理量值，但目前获取的数据为count值，所以需要在后端服务器软件进行实际物理量转换。在本实例中，可获取前向积分常数X₀（FFFFFFA9）为十进制−87，反向积分常数X_n（FFFFFF7A）为十进制−134。随机抽取W₆为样例，W₀和W₆数据解析及新合成的W₆数据如图4所示。

图 4  W₀和W₆数据解析及新合成的W₆数据表

Figure 4.  W₀ and W₆ data analysis and newly synthesized W₆ data

下载: 全尺寸图片 幻灯片

由图4可知，C₀、C₁、C₂的内容均为00，代表了非压缩数据信息，同时从W₀中获取W₆的广义字节编码C₆=10，如果编码为10，还需进一步从W₆中获取其自身的细解码dnib=11，综合可知W₆的bit2到bit32包含3个10位差值，按每10位拆分即可获得3个差值，即1111110011、0010100110、1110100011，换算成十进制分别为−13、166、−93，然后根据前面累加得到的原始值推算出本段压缩数据的原始值，其他的W_k压缩数据按照此方法依次拆分即可。

利用上述分析结果编写代码后才能实用化，对数据进行逐个解析，首先确定各数据的定位，需设置一个参考起点，根据第44至第45个字节的内容获取数据区的开始偏移量并进行定位。笔者通过第15至第17个字节的通道标识符进行定位，其目的是按通道解压数据，并形成相应通道的文本文件，以配套CAL79_噪声功率谱计算软件使用。

程序流程如图5所示，分为数据读取、字节序转换、解压、校验、文本文件形成等部分，此处着重对数据拆分的实现过程进行说明。

图 5  程序流程

Figure 5.  Program flow chart

下载: 全尺寸图片 幻灯片

int unzip（int * data, FILE *ftxt, int i）//解压子程序。{……//获取样本数目，MI3000C一般为100个。　　samples = （data[i + 3]） & 0x0000FFFF;//获取仪器设备的采样率，一般为每秒100次。　　rate = （（data[i + 4]） >> 16） & 0x000000FF;//获取本记录的长度。　　byte_long = （（data[i + 9]） >> 8） & 0x000000FF;//获取数据帧个数用于循环　　frame_num = （pow（2, （double）byte_long）） / 64−1;i = i + 12; //i指向第一个数据帧的开始即W0处　　begin_data = data[i + 1];//提取X0　　last_data = data[i + 2];//提取Xn　　num = 0;　　for （size_t frame = 0; frame < frame_num;frame++）//以数据帧为单位的大循环　　{　　//从W0提取C0到C15广义字节编码，分别对应W0到W15。　　　　　　for （int j = 0; j < 16; j++）　　　　　　　{C[j] = （data[i] >> （32−2 * （j + 1））） & 0x03;}　　　　i++;　　// W1到W15进行循环，Ck和Wk一一对应。　　　　for （int j = 1; j < 16; j++, i++）　　　　{　　　　　　//获取W1到W15之一的最左侧两位“细节码”。　　　　　　dnib = （data[i] >> 30） & 0x03;　　　　　　if （C[j] == 0x02）//旗下包含1个30位差值或2个15位差值或3个10位差值　　　　　　{　　　　　　　　if （dnib == 0x01）//1个30位差值标记　　　　　　　　{　　　　　　　　　　unzip[num] = data[i] & 0x3FFFFFFF;//按位提取　　　　　　　　　　if （（unzip[num] & 0x20000000） == 0x20000000）　　　　　　　　　　　　{ unzip[num] = unzip[num] |0xC0000000; }//负值处理　　　　　　　　　　num++;　　　　　　　　}　　　　　　　　else if （dnib == 0x02）//2个15位差值标记　　　　　　　　{……}　　　　　　　　else if （dnib == 0x03）//3个10位差值标记　　　　　　　　{……}　　　　　　}　　　　　　　……//其它广义字节编码判断　　　　}　　}}　　……　　//差值转原始值、数据校验、按某通道保存文本文件　　……

由程序代码可知较关键的步骤为数据定位、以数据帧为单位执行循环、采用“按位与”的方法提取信息、对出现负值的情况给予处理。笔者通过大量二进制的位操作提高执行效率。

经实际测试，该软件可对河南省一般站MI3000-C、REMOS-SIT4烈度计产出的小时波形文件进行分道解码，且可对基准站和基本站使用的EDAS-24 GN（EEW）、TDE-324 FI/CI数据采集器产出的小时波形文件进行分道解码，并在程序解码过程中设置了数据验证功能，即前向积分常数X₀加各项差值等于反向积分常数X_n，程序判断为真才能输出结果，增加了数据解码的准确性。在获取解压后的分道文本数据后，可顺利使用噪声功率谱计算软件进行台站背景噪声环境评估。

3.   结论

本文对预警项目台站波形数据包格式进行了研究，得出以下结论：

（1）预警项目对数据采集器的产出波形数据格式做出了统一规定，对数据的部分字节内容进行了重新定义，准确把握各字节位的定义和关系是数据解析与使用的前提。

（2）连续波形数据包结构中固定头段区、可变头段区和数据区是有机整体，解压算法需相关区域的参数参与计算，解压的关键在于二进制数据按位拆分。

（3）列出了实用的程序代码，利用C语言的位操作方法提取有效信息，并拆分数据。笔者查阅了部分文献，发现针对地震数据包结构和压缩的原理描述较多，且较抽象，可操作性不强，具体实现过程多用一个函数包一笔带过，对实现细节的解读较少，本文给出了具体实现方法，稍加修改即可运行。