0 引言

通过建立自动高效且标准通用的途径，将长期、大量的高质量地震波形数据用于地震活动性和危险性评估、地震破裂过程以及地球内部结构等研究，是地震工作者和地震数据中心的共同目标。自中国地震局“十五”项目以来，全国建成的数字地震台站已达1100个左右(刘瑞丰等，2008)。地震台站产出的波形数据经实时数据流传输至省监测台网中心和中国地震台网中心，用于各省和全国的地震速报和编目工作。中国地震台网中心的波形数据以miniSEED标准格式(中国地震局，2003)的文件进行永久存储，通过国家地震科学数据中心和国家数字测震台网数据备份中心(郑秀芬等，2009)共享给地震科技人员研究使用。

中国地震台网中心每天存储的波形数据量约为43GB，2008年以来累计存储的地震波形数据总量已接近150TB(图 1)，数据文件近1400万个。

随着人工智能深度学习等依赖数据驱动的新技术的广泛应用，越来越多的地震学研究需要大量波形数据的支撑(Ross et al，2018)。面对多年积累的大量波形数据，国家地震科学数据中心以往主要通过手动截取和人工拷贝方式提供服务，仅能服务于数据需求量较少的研究工作，并不能满足较短时间内获取大量波形定制数据的用户需求，因此迫切需要为其存储的波形数据集建设一个可提供自动、7×24h不间断服务的新平台，以提高波形截取和数据服务能力，满足用户需求。

国际数字地震台网联盟(International Federation of Digital Seismograph Networks，FDSN)制定了地震数据交换的Web服务接口标准(International Federation of Digital Seismograph Networks，2019)，规范了服务的名称、查询参数和返回结果。目前，美国IRIS数据中心和德国GEOFON台网中心等约20个国家或地区的地震数据中心已提供FDSN标准的Web服务接口。地震数据服务接口全球标准化的实现大大降低了数据使用者获取数据的难度，研究人员使用Obspy(Krischer et al，2015)等开源地震学工具，可以从全球任何实现FDSN标准Web服务接口的地震数据中心获取数据。

2015年，中国地震台网中心“地震科学数据共享”项目“地震台网产品研发及共享等应用研究”子课题建成了地震数据服务云平台，实现了2013—2015年中国地震台网的波形数据、2009—2015年中国地震统一编目中地震目录数据以及台站参数等元数据的FDSN标准Web服务接口，其以私有云的方式(MacCarthy et al，2020)在地震行业内网布设，可根据数据请求量的多少弹性调整云平台上Web服务进程的数目。

我们对中国地震台网中心存储的2008—2021年波形数据进行了归档和整理，同时对整理后数据的可用率和断记缺失情况进行调查统计。在对全部数据文件建立索引数据库的基础上，搭建符合FDSN国际标准Web服务接口的波形数据服务平台，并通过具体应用实例对平台的稳定性和服务效率进行了初步测试。

1 地震波形数据的归档和整理

全国各省监测台网的实时波形数据包每天不间断地传输至中国地震台网中心的JOPENS流服务器上，并以简单追加的方式存储在miniSEED格式的磁盘文件中作永久保存。虽然这种保存方式可以确保数据不会丢失，但并不能改变数据传输过程中出现的重复和乱序情况，因此初始存储下来的波形数据需要进行归档和整理。

归档地震波形数据的文件系统包含了一系列有规律组织起来的文件夹和数据文件(图 2)，其中数据文件依照“台网代码. 台站代码. 位置号. 通道代码. 数据类型. 年. 儒略日”方式来命名，根据文件名称可确定其所在的位置路径为“年/台网代码/台站代码/通道代码. 数据类型/”。将初始存储文件中的miniSEED数据包逐一取出，根据数据包头段中的台网、台站、位置和通道类型的信息和采样点的时间范围，将其分发到归档系统中名称匹配的数据文件中。由于归档过程是按照初始存储文件中数据包下载时的排列顺序进行分发，因此初始存储文件中出现重复和乱序的数据包依然存在于归档系统中的数据文件中，需要进一步对归档系统中的文件进行去重和排序整理。

以贵州省地震台网毕节台站垂直方向(GZ.BJT.00.BHZ) 2019年6月17日14：56：00—14：57：30之间的原始存储数据波形记录(图 3(a))和全部92个miniSEED数据包的波形(图 3(b))为例，数据包波形从上到下按存储下载时的顺序排列，可以看到有25个迟到的冗余数据包以及1.52s的数据缺失。对2008—2021年间归档完成后的波形数据文件逐一进行如下整理：以文件中数据包为单位，按数据包内最后采样点的时间对数据包进行排序，剔除重复数据包。数据去重时优先保留采样点数较多的数据包，且整理过程并未对波形数据进行重新打包或重新采样等操作。图 3(c)为图 3(b)中数据包经过整理后的结果，波形从上到下按照数据包起始时间顺序排列，其中25个迟到的冗余数据包已经被剔除，仅剩缺失的数据待补足。

2 地震波形数据的可用率和缺失断记情况统计

对归档并整理后的连续波形数据文件中的采样点数、断记数目和断记时长进行统计。图 4展示了每日采样点数目的变化情况，可以看出，2011年前波形数据存储的缺失情况较多，可用率偏低，2012年以后数据可用率均保持在93%以上。

2008—2021年归档台站通道(由台网名、台站名、位置号和通道代码四项联合指定)总数为4382个，其中运行时间超过1年的台站通道数3944个。将台站通道数据可用率定义为某一“台网. 台站. 位置. 通道”指定数据流在其运行时间内的实际存有采样点总数与应有采样点总数的比值。运行时间超过1年的台站通道数据可用率的分布情况如图 5所示，其中可用率大于等于0.9的台站通道数为2761个，占比70%；可用率大于等于0.8的台站通道数为3326个，占比84.3%；可用率低于0.6的台站通道数为219个，占比6.6%。

我们计算了某一台站所有通道数据可用率的平均值，记为台站数据可用率。运行时间超过1年的台站数据可用率的空间分布情况，如图 6所示，从图中可见西藏(XZ)台网和贵州(GZ)台网中台站可用率偏低的台站数较多，此外，北京和天津部分配备短周期地震仪的台站可用率在70%~80%之间。

对每年波形数据的缺失断记情况按6个时长范围进行统计，结果如图 7所示。可以看出，无论是波形的缺失断记数目，还是断记的时间长度，均有年度的变化。在2014年之前，1~10s长度的数据缺失数目最多，2014年之后1s时长以内的缺失数目逐渐增多，特别是2018年和2020年的波形数据中出现大量1s内的断记，2021年数据整体缺失断记数目最少。

对断记个数和时长的关系进行进一步研究，若定义N_Gap(T)为2008—2021年间数据缺失断记时间长度大于等于T的数据缺失断记数目，N_Gap(T)随着T的减小而增大，两者的关系在双对数坐标系下呈现明显的线性分布特征(图 8)，可表示为

$ \lg N_{\text {Gap }}=-0.71 \times \lg T+8.47 $

(1)

该特征符合统计分形的碎形分布(陈颙等，2005；栾锡武等，2004)，表明了导致数据断计的主要原因为随机原因。

3 波形数据服务新平台的搭建与测试 3.1 波形数据服务新平台的建设

归档并整理后的波形数据文件全部保存在国家地震科学数据中心服务器的磁盘阵列中，该服务器配置有2个8核CPU(主频为2.1GHz)、64G内存和120个容量为8TB的SAS硬盘，组装了RAID5磁盘阵列。波形数据服务新平台是一个布设在上述服务器上的软件系统，7×24h不间断地运行，自动回应用户的数据请求(图 9)。

平台采用美国IRIS数据中心开源的portable-fdsnws-dataselect程序作为数据与用户中间层的服务软件。为了从大量连续波形数据中快速找到用户所需的数据，需要在数据库中对归档系统中全部连续波形数据建立两个索引表。图 10中展示了第一个索引表中的一行内容，为贵州毕节台2019年6月17日的波形数据文件GZ.BJT.00.BHZ.D.2019.168内数据片段的索引情况，包括该数据片段ID、时间范围、采样率、所属文件及其路径、该片段起始点在所属文件中的字节偏移量及其所占字节数等。另一个索引表为Summary表，记录了每个数据ID对应波形数据的总起止时间。

最后，选择在平台的服务端运行portable-fdsnws-dataselect，其提供了国际标准的FDSN Web服务接口，该软件通过调取索引数据库，返回符合用户请求的某台站通道在某时间段内的miniSEED波形数据，通过HTTP协议直接传输至客户端的应用程序。

3.2 中国地震台网波形数据服务平台的测试

为验证中国地震台网波形数据服务平台的可靠性、高效性和应用场景的多样性，搭建了conda环境，安装地震学的开源python工具包Obspy(Krischer et al，2015)，利用Obspy中的FDSN Web服务客户端对事件波形和震相波形的批量下载、震相智能检测、连续波形下载及其日能量比计算进行了测试。

(1) 测试用例1：下载2008—2021年836个震级≥4.5、震中距≤200km台站的三分量波形数据，起止时间为发震时刻前1min至震后5min，使用的输入文件包括events.txt和stations.txt。其中，events.txt为事件列表，包含5列数据，分别是发震日期YYYY-MM-DD、发震时间HH：MM：SS.SS、震中纬度、震中经度、震中深度和震级；stations.txt为台站列表，包含三列数据，分别为台站ID(台网名. 台站名)、台站纬度和台站经度。测试耗时36min，测试python脚本如下：

#!/usr/bin/env python3

from obspy.clients.fdsn import Client

from obspy import read, UTCDateTime, Stream

import numpy as np

import glob, os

from obspy.geodetics import degrees2 kilometers

from obspy.taup.taup_geo import calc_dist_azi

client=Client("http://服务器IP地址：端口号, service_mappings={′event′：None, ′station′：None})

#读取事件参数列表，格式为5列，分别为：发震日期YYYY-MM-DD发震时刻HH：MM：SS.SS事件纬度事件经度事件深度震级

orid, orit, evla, evlo, evdp, m=np.loadtxt(′events.txt′, dtype=bytes, usecols=(0, 1, 2, 3, 4, 5), unpack=True).astype(str)

evla=evla.astype(float)

evlo=evlo.astype(float)

m=m.astype(float)

#读取台站参数列表，格式为3列，分别为：台网名. 台站名. 位置号. 通道号台站纬度台站经度

stid, stla, stlo=np.loadtxt(“stations.txt”, dtype=bytes, usecols=(0, 1, 2), unpack=True).astype(str)

stla=stla.astype(float)

stlo=stlo.astype(float)

for i in range(len(orid))：

ori_part1=orid[i].replace(′-′, ″)

ori_part2=orit[i].replace(′：′, ″)

evid=ori_part1+ori_part2+′_′+str(m[i])

os.makedirs(′./′+evid, exist_ok=True)

ori=UTCDateTime(orid[i]+′T′+orit[i])-8*3600

t1=ori-1*60

t2=ori+5*60

for j in range(len(stid))：

net=stid[J].split(′.′)[0]

stnm=stid[J].split(′.′)[1]

gcarc=calc_dist_azi(evla[i], evlo[i], stla[J].stlo[J].6371.0, 0)[0]

dist=degrees2 kilometers(gcarc)

dist=round(dist, 2)

if dist＜=200：

for ch in [′*E′, ′*N′, ′*Z′]：

try：

st=client.get_waveforms(network=net, station=stnm, location=′00′, channel=ch, starttime=t1, endtime=t2)

except：

print(“No data available for”+net+′.′+stnm+′ between′+str(t1)+′ and′+str(t2))

else：

stream_name=st[0].id+′_′+str(dist)

st.write(′./′+evid+′/′+stream_name+′.mseed′, format=′MSEED′)

测试用例1所下载的部分波形数据，如图 11所示。

(2) 测试用例2：根据中国地震统一正式目录，下载2017年震级≥3.0、震中距≤500km的初至P震相波形6501条(占应下载波形数据的96.5%)，波形长度为P震相到时前2min至震后2min，保存至本地磁盘，共耗时11.2h，数据大小为16GB。应用人工智能震相检测模型EQTransformer(Mousavi et al，2020)对下载波形进行检测，主要参数设置分别为：detection_threshold=0.1，P_threshold=0.1，S_threshold=0.1。P波初至震相检出率为94.3%，检测到时与目录中人工拾取到时之差的统计分布如图 12所示，其中到时差在0.1s内的震相占比68.4%，到时差在1s内的震相占比89.5%。

(3) 测试用例3：下载2019年6月25日云南台网56个台站一天的三分量连续波形，并在内存中直接计算每5min的三分量波形数据能量均值之比(E/Z，N/Z，E/N)(Pedersen et al，2020)，取每天的中位数作为当天的能量比，用于评估该台网在当天信号质量。下载和处理数据共耗时6min，统计结果见图 13。

其中，云南(YN)台网鹿泉台(LUQ)N/Z比值为72.2，显示出异常，为此下载该台站6月全部的三分量波形数据，并计算了每日能量比(图 14)。由图 14可见该台站的能量比仅在6月25日一天内出现异常，图 15展示的BHN分量波形数据也证实了该异常。

4 结语

在对2008—2021年近150TB的中国地震台网连续波形数据归档和整理的基础上，建设波形数据服务平台，该平台采用了国际标准的Web服务接口，实现了用户通过应用程序从中国地震台网中心自动获取波形数据的功能。

此外，对波形数据可用率和断记情况进行统计，结果为平台用户更有效地获取数据提供了参考，同时有助于跟进数据的查漏补缺工作。未来我们将继续开展更详细的波形数据质量评估，并搭建震源参数和地震台站元数据的自动服务接口。

本文采取将波形数据以文件形式存储至一个磁盘阵列的方式，该方式在进行大量数据处理时存在瓶颈。因此，利用大数据技术提供分布式的数据查询下载或在集群服务器上使用多个处理器进行数据的分布式处理(陈通等，2022)，是中国地震台网波形数据服务未来的发展方向。

致谢: 在整理波形数据和搭建服务平台的过程中，使用了IRIS数据中心(https://github.com/iris-edu)的多个开源程序，包括dataselect、mseedindex、msi和portable-fdsnws-dataselect等，在此表示感谢。