地震波形数据存储解决方案

转到 MiniSeed 教程 dolphindb/DolphinDBPlugin，在 bin 目录下，找到需要的插件版本并下载。通过以下方式加载插件：

 loadPlugin("./plugins/mseed/PluginMseed.txt");

MiniSeed 插件提供了 read、write、parse、parseStream 等接口，为 MiniSeed 文件的读、写、解析提供了强大的帮助，用户只需加载 MiniSeed 插件，就能快速地调用这些接口。

存储采样历史数据的数据库设计方案如下：

• 选择 TSDB 存储引擎，采用组合分区，一级分区按日期进行值分区，二级分区按设备 Id 进行值分区，每个分区存储一个 id 一天的历史数据
• 分区列为 ts 和 id，索引列为 id 和 ts。ts, id, value 列均采用 Delta（delta-of-delta encoding） 压缩算法，每个分区内 sortKey 相同的数据全部保留

通过以下方式创建分布式数据库和分区表：

  //创建储实时数据的分布式数据库和分区表
if( existsDatabase("dfs://real") ){ dropDatabase("dfs://real") }
create database "dfs://real" partitioned by VALUE(2023.03.01..2023.03.10), VALUE(1..3900), engine='TSDB'
create table "dfs://real"."realData"(
	id INT[compress="delta"],
	ts TIMESTAMP[compress="delta"],
	value INT[compress="delta"]
)
partitioned by ts, id,
sortColumns=[`id, `ts],
keepDuplicates=ALL

  //创建存储时延计算结果的分布式数据库和分区表
if(existsDatabase("dfs://delay")){dropDatabase("dfs://delay")}
create database "dfs://delay" partitioned by VALUE(2023.03.01..2023.03.10), HASH([INT, 10]), engine='TSDB'	
create table "dfs://delay"."delayData"(
	id INT[compress="delta"],
	tagid SYMBOL,
	startTime TIMESTAMP[compress="delta"],
	receivedTime TIMESTAMP[compress="delta"],
	delay INT[compress="delta"]
)
partitioned by startTime, id,
sortColumns=[`id, `startTime],
keepDuplicates=ALL			

  //创建存储台网、台站、位置、通道等基础信息的维度表 
if(existsTable("dfs://real","tagInfo")){ dropTable(database("dfs://real"),"tagInfo") }
create table "dfs://real"."tagInfo"(
	id INT[compress="delta"],
	net SYMBOL,
	sta SYMBOL,
	loc SYMBOL,
	chn SYMBOL,
	tagid SYMBOL
)
sortColumns=[`id]

库表创建后，可以先行模拟台网、台站、位置、通道等基础信息的维度表数据，示例代码如下：

// 导入台网、台站、位置、通道等基础信息的维度表信息数据
net = ["ZJ","YN","XZ","XJ","TJ"]
sta = ["A0001","A0002","A0003","A0004","A0005","A0006","B0001","B0002","B0003","C0001"]
tmp = `EIE`EIN`EIZ
netList = stretch(net,150)
staList = take(stretch(sta,30),150)
locList = take(`40,150)
chn = take(tmp,150)
colt =   array(STRING)
for(i in 0..(chn.size()-1)){
	colt.append!( chn[i].split()[0] + "_" + chn[i].split()[1] + "_" +chn[i].split()[2] )
}
tagid = "XFDSN:"+netList+"_"+staList+"_"+locList+"_"+colt
t = table(1..150 as id,netList as net,staList as sta,locList as loc,chn,tagid)
pt = loadTable("dfs://real","tagInfo")
pt.append!(t)

流数据表用于接收实时流数据以及进行时延计算，定义四张流数据表，分别用于接收采样信息、MiniSeed 文件信息、时延计算结果以及异常告警信息：

• dataStream：接收采样数据，其结果有三个走向，分别是存入分布式数据库、接受外部 API 的订阅、实时展示
• metaStream：保存 MiniSeed 插件解析文件的元数据（如文件的数据行数、解析的时间、文件的时间等），对元数据进行时延计算后将结果注入到 delay 结果表中
• delayStream：时延计算结果表
• abnormalStream：异常告警信息结果表

上述流数据表创建及持久化方式如下：

  //删除流数据表及流数据引擎
unsubscribeTable(tableName = `dataStream,actionName = `append_data_into_dfs)
unsubscribeTable(tableName = `dataStream,actionName = `abnormalDetect)
unsubscribeTable(tableName = `metaStream,actionName = `calculate_delay)
unsubscribeTable(tableName = `delayStream,actionName = `append_delay_into_dfs)	
try{ dropStreamTable(`metaStream) }catch(ex){ print(ex) }
try{ dropStreamTable(`dataStream) }catch(ex){ print(ex) }
try{ dropStreamTable(`delayStream) }catch(ex){ print(ex) }
try{ dropStreamEngine(`engine) }catch(ex){ print(ex) }
try{ dropStreamTable(`abnormalStream) }catch(ex){ print(ex) }	

  //创建建立接收实时流数据的流数据表
st1 = streamTable(1000000:0,`id`tagid`startTime`receivedTime`actualCount`expectedCount`sampleRate,[INT,SYMBOL,TIMESTAMP,TIMESTAMP,INT,INT,DOUBLE])
st2 = streamTable(1000000:0,`id`ts`data,[INT,TIMESTAMP,INT])
st3 = streamTable(1000000:0,`id`tagid`startTime`receivedTime`delay,[INT,SYMBOL,TIMESTAMP,TIMESTAMP,INT])
st4 = streamTable(10000:0,`Time`id`anomalyType`anomalyString,[TIMESTAMP,STRING,INT,STRING])
enableTableShareAndPersistence(table=st1, tableName=`metaStream, asynWrite=true, compress=true, cacheSize=2000000, preCache = 100000)
enableTableShareAndPersistence(table=st2, tableName=`dataStream, asynWrite=true, compress=true, cacheSize=2000000000, preCache = 100000)
enableTableShareAndPersistence(table=st3, tableName=`delayStream, asynWrite=true, compress=true, cacheSize=2000000, preCache = 100000)
enableTableShareAndPersistence(table=st4, tableName=`abnormalStream, asynWrite=true, compress=true, cacheSize=2000000, preCache = 100000)

dataStream 表结构见表 2，metaStream、delayStream、abnormalStream 表结构分别见表 4、表 5、表 6。

调用 MiniSeed 插件的 write 函数向 MiniSeed 文件写入采样数据，模拟地震计源源不断地生成 MiniSeed 实时流过程。以下代码以 submitJob 提交任务的形式模拟通道“ZJ_A0001_40_E_I_E”一直往 Miniseed 文件写入实时流数据的场景。

def writeMseed(){
	/*
	 * Description：
	 * 	此函数用于模拟实际生产环境中，每3s向台网中心发送 MiniSeed 数据包的过程，即每隔3s向 aimPath 路径上的 MiniSeed 文件写数据
	 */
	sid = "XFDSN:ZJ_A0001_40_E_I_E"
	aimPath = "../streamMiniSeed/ZJ_A0001_40_E_I_E.20000101.mseed"
	sampleRate = 100
	cnt = 0
	do{
		startTime = now()
		mseed::write(aimPath, sid, startTime, sampleRate, rand(-3000..500,300))
		sleep(3000)
		cnt += 1
	}while(cnt < 40)	
}
	

mseedStreamJobid = submitJob("wrtieMseedStream","wrtieMseedStream",witeMseed) //提交后台任务
/*
 * 使用以下方式取消实时写入任务
 * cancelJob(mseedStreamJobid)
 */

调用 MiniSeed 插件的 write 函数模拟历史数据。以下代码以 submitJob 的方式提交了模拟 2023.03.02 所有通道 1 天内波形数据的任务。

def writeMseed(filePath,sid,startTime,sampleRate,valueList){
	/*
	 * Description：
	 * 	此函数用于将一段连续的采样值写入到 MiniSeed 文件
	 * Input：
	 * 	filePath：STRING MiniSeed 文件存储路径
	 * 	sid：STRING 文件块id 
	 * 	startTime：TIMESTAMP 开始时间
	 * 	sampleRate：DOUBLE 采样频率
	 * 	valueList：INT VECTOR 采样值向量
	 */
	mseed::write(filePath, sid, startTime, sampleRate, take(valueList,8639999))
}

def parallelWrite(){
	/*
	 * Description：
	 * 	此函数用于多线程 写 MiniSeed 文件
	 */
	realDbName,dtName = "dfs://real","tagInfo"
	filePathList = exec "../miniSeed/"+net+"."+sta+"."+loc+"."+chn+"."+"20230302.mseed" from loadTable(realDbName,dtName)
	sidList = exec tagid from loadTable(realDbName,dtName)
	startTime = 2023.03.02T00:00:00.000
	sampleRate = 100
	ploop(writeMseed{,,startTime,sampleRate, rand(-3000..500,300)},filePathList,sidList)
}

submitJob("writeMseed","writeMseed",parallelWrite)

调用 parseStream 函数多线程解析 MiniSeed 文件，该接口返回一个包含三个键值对的字典 ret：ret[`data]、 ret[`size] 和 ret[`metaData]，其中，ret[`data]为采样数据，ret[`metaData]为 MiniSeed 文件块信息。以下代码以 submitJob 的方式提交了解析所有通道 2023.03.02 1 天内波形数据的任务。

def parseAndImportIntoDfs(realDbName,realTbNname,dtName,tagInfo,fileParse){
	/*
	 * Description：
	 * 	此函数用于解析 miniSeed 文件，并将结构化数据存入分布式数据库
	 * Input：
	 * 	realDbName,realTbNname,dtNam：STRING 
	 * 	tagInfo：dict 键为台网、台站、位置通道拼接起的字符串，值为对应的 id
	 * 	fileParse：STRING 解析的miniSeed文件存储路径
	 */
	ret = mseed::parseStream(file(fileParse).readBytes(30000000))
	data = ret[`data]
	data.replaceColumn!(`id,tagInfo[data.id])
	delete from data where id = NULL
	pt = loadTable(realDbName,realTbNname)
	pt.append!(data)	 	
}

def parallelInToDfs(realDbName,realTbNname,dtName,tagInfo,filePathList){
	/*
	 * Description：
	 * 	此函数用于并行解析 MiniSeed 文件，并存入分布式数据库
	 * Input：
	 * 	realDbName,realTbNname,dtName：均为 STRING 类型常量，分别代表分布式数据库名、分区表名、维度表名
	 * 	tagInfo：dict 键为由台网、台站、位置、通道组成的字符串，值为对应的id
	 * 	filePathList：STRING VECTOR 需要解析的MiniSeed文件路径向量
	 */
	ploop(parseAndImportIntoDfs{realDbName,realTbNname,dtName,tagInfo,},filePathList)
}

realDbName,realTbNname,dtName = "dfs://real","realData","tagInfo"
filePathList = exec "../miniSeed/"+net+"."+sta+"."+loc+"."+chn+"."+"20230302.mseed" from loadTable(realDbName,dtName)
filePathListGroup = cut(filePathList,48)
d = exec tagid,id from loadTable(realDbName,dtName)
tagInfo=dict(d[`tagid],d[`id])
for(i in 0..(filePathListGroup.size()-1)){
	submitJobEx("parse_mseed_into_dfs"+string(i),"parse_mseed_into_dfs"+string(i),i%10,48,parallelInToDfs,realDbName,realTbNname,dtName,tagInfo,filePathListGroup[i])
}	

将历史数据解析入库并强制刷盘后，观察磁盘空间占用，其磁盘空间占用大小为 1.23GB，观察原始 MiniSeed 文件大小为 3.12GB，由此可见，MiniSeed 文件解压后存储为结构化数据的体积大大减少了。

通过流数据表的发布 - 订阅功能可实现对流数据的处理。对于采样数据，通过内置函数 subscribeTable 订阅 dataStream，实现入库存储。subscribeTable 函数调用方式如下：

subscribeTable([server],tableName,[actionName],[offset=-1],handler,[msgAsTable=false],[batchSize=0],
[throttle=1],[hash=-1],[reconnect=false],[filter],[persistOffset=false],[timeTrigger=false],
[handlerNeedMsgId=false],[raftGroup],[userId=””],[password=””])

对于时延计算，首先，定义时延计算函数 calDelay

def calDelay(msg){
	/*
	 * msg：数据流
	 * delayTableName：时延结果流数据表
	 */
	t = select first(id) as id,first(tagid) as tagid,last(startTime) as startTime,receivedTime,last((receivedTime+8*60*60*1000+3000)-(startTime+(actualCount-1)*10)) as delay from msg group by receivedTime
	objByName(`delayStream).append!(t) 	
}

其计算规则如下：

receivedTime-(startTime+(actualCount)*10)

举例说明：receivedTime 为"2023.03.09T12:33:03.995"，startTime 为"2023.03.09T12:33:00.984"，actualCount 为 300，则该 MiniSeed 文件块的时延为 2023.03.09T12:33:03.995-(2023.03.09T12:33:00.984+300*10) = 11ms。

然后，再通过 subscribeTable 函数订阅 metaStream，将计算结果注入 delayStream，实现时延计算。时延计算结果示例如表 9 所示：

DolphinDB 内置异常告警引擎 createAnomalyDetectionEngine，其调用方式如下：

createAnomalyDetectionEngine(name,metrics,dummyTable,outputTable,timeColumn,[keyColumn],[windowSize], 
[step],[garbageSize],[roundTime=true],[snapshotDir],[snapshotIntervalInMsgCount],[raftGroup])

对于异常告警，首先根据 createAnomalyDetectionEngine 函数定义异常告警引擎 engine，其告警规则为：

[max(data)>= 300,min(data)<=-2500]

即，3 秒长度的时间窗口内，若采样数据最大值大于 300 或采样数据最小值小于 -2500，则告警。

然后通过 subscribeTable 函数订阅 dataStream，将告警结果注入 abnormalStream，实现异常告警。异常告警结果示例如表 10 所示：

在 abnormalStream 中，若同一个 id 在同一个时间戳内出现了两条数据，则说明，在 3s 内，同时出现了波形数据的最大值大于 300、最小值小于 -2500 的现象。

实时流数据存储、时延计算、异常告警代码如下所示：

//实时流数据入库存储
subscribeTable(tableName=`dataStream, actionName="append_data_into_dfs", offset=0, handler=loadTable("dfs://real","realData"), msgAsTable=true,batchSize=100000, throttle=1, reconnect=true)
//时延计算
subscribeTable(tableName=`metaStream, actionName="calculate_delay", offset=0, handler=calDelay, msgAsTable=true,batchSize=100000, throttle=1, reconnect=true)
//时延计算结果入库存储
subscribeTable(tableName=`delayStream , actionName="append_delay_into_dfs", offset=0, handler=loadTable("dfs://delay","delayData"), msgAsTable=true,batchSize=100000, throttle=1, reconnect=true)
//建立异常告警引擎，设置告警规则
engine = createAnomalyDetectionEngine(name="engine", metrics=<[max(data)>= 300,min(data)<=-2500]>, dummyTable=objByName(`dataStream), outputTable=objByName(`abnormalStream), timeColumn=`ts, keyColumn=`id, windowSize =3000 , step = 3000)
//订阅实时流数据，进行异常告警
subscribeTable(tableName="dataStream", actionName="abnormalDetect", offset=0, handler=append!{engine}, msgAsTable=true)	

历史数据查询与导出测试 case 如下：

• case1：某个台网、某个台站一天所有通道的历史观测数据降采样查询，精确到分钟
• case2：某个事件发生时刻前 2 分钟、后 11 分钟历史观测数据查询，关联台站为 20-30 个
• case3：某个台网、某个台站一天内所有通道数据导出为 MiniSeed 文件
• case4：Java API 查询，某个台网、某个台站一天所有通道的历史观测数据

性能测试结果如表 11 所示：

查询 case1-3 代码如下所示：

netAim,staAim,locAim = `ZJ,`A0003, `40  //查询case的查询条件
// case1
timer(10){
	aim  = exec id from loadTable("dfs://real","tagInfo") where net = netAim and  sta = staAim and loc = locAim
	t = select avg(value)  as sample_value from loadTable("dfs://real","realData") where ts >= 2023.03.02T00:00:00.000 and ts < 2023.03.03T00:00:00.000 and id in  aim   group by bar(ts,60000),id
}

// case2
timer(10){
	aim = exec id from loadTable("dfs://real","tagInfo") where net = netAim
	t = select *  from loadTable("dfs://real","realData") where ts >= 2023.03.02T00:08:00.000 and ts <= 2023.03.02T00:21:00.000 and id in aim		 	
}

// case3
timer {
	aim = exec id from loadTable("dfs://real","tagInfo") where net = netAim and sta = staAim
	t = select * from loadTable("dfs://real","realData") where ts >= 2023.03.02T00:00:00.000 and ts < 2023.03.03T00:00:00.000 and id in aim
	t = select ts,loc,net,sta,chn,value,tagid from ej(loadTable("dfs://real","tagInfo"),t,`id) // 导出的结果集数据
	
	sidList = exec distinct tagid from t // 涉及的所有sid
	startTimeTable = select min(ts) as startTime from t group by tagid order by tagid desc  // 文件开始采样时间
	dataList = [] // 涉及的所有采样值
	for(i in 0..(sidList.size()-1)){
		tmpValue = exec value from t where tagid = sidList[i]
		dataList.append!(tmpValue)
	}
	
	aimFilePath= "<YOURDIR>/"+sidList+'.20230302.mseed' // 目标文件
	sampleRate = 100.00 // 采样频率，表示1秒钟有100个采集点，即10ms采集一次
	ploop(mseed::write{,,,sampleRate,},aimFilePath,sidListid,startTimeTable[`startTime],dataList) // 使用mseed::write函数写入文件，使用ploop并行加速
}

结构化实时流数据可视化展示通过 grafana 实现。在 grafana 中，可通过编写简单的查询语句，便能获取需要展示的数据。grafana 中执行以下查询语句对采样数据实时流进行可视化展示：

//波形实时流数据可视化展示
real_db_name,dt_name = "dfs://real","tagInfo"
t = select avg(data)  as value from dataStream group by bar(ts,6000) as ts, id
select t.value from t left join loadTable(realDbName,dtName) rt on t.id = rt.id  pivot by t.ts,rt.tagid

执行以下查询语句对时延数据进行可视化展示：

//时延流数据可视化展示
select delay from delayStream pivot by startTime,tagid

• 《大数据架构下的地震波形数据分析应用浅析》
• 《数字地震台网综合处理系统通用架构研究与设计》
• 《地震台网智慧服务平台的开发与应用》
• DolphinDB 流数据教程
• DolphinDB 分布式存储教程
• DolphinDB-grafana-datasource 插件教程

用户可按照行文顺序运行文章中所列出的脚本复现解决方案，此外为方便用户复现，本节对文章中所涉及的脚本整理成文件，用户只需下载并执行附件压缩包中 1~8 的 txt 文件，然后按照如下所示脚本调用文件中的函数即可复现本文的解决方案：

login(`admin,`123456)
undef(all);
clearAllCache();
go;

/***
 * ==============================
 * 在运行此脚本此之前，需要做以下事情：
 * 	1.建立与server同级的目录  miniSeed
 * 	2.建立与server同级的目录  streamMiniSeed
 * 	3.建立与server同级的目录  outputMiniSeed
 * 	4.在./server/plugins/目录下建立mseed目录，并将libPluginMseed.so文件和PluginMseed.txt文件放置在改目录下
 * 	5. 运行脚本1-8
 * ==============================
 */



netAim,staAim,locAim = `ZJ,`A0003, `40  //查询case的查询条件


  //step1：加载插件、创建分布式库表、流数据表持久化
mainPrepare()

  //step2：MiniSeed实时流数据模拟与解析
mainStreamSimulate()

  //step3：实时流入库、时延计算及异常告警
mainBuidSubscribe()

  //step4：MiniSeed历史数据模拟( 耗时1min左右，该步骤后台作业完成后再执行step5)
mainHistoryMiniSeedSimulate()


  //step5：MiniSeed历史数据入库( 耗时5min左右，该步骤后台作业完成后再执行step6和step7)
mainHistoryMiniSeedParse()


  //step6：磁盘空间占用统计(单位：GB)
mainDiskUsage()

  //step7：查询case1-case4
queryRecord =  mainQueryTest(netAim,staAim,locAim)

  //step8：查询case5
outputRecord = mainOutputMiniSeed(netAim,staAim)

测试脚本压缩包：test_scripts.zip