用时序数据库 DolphinDB 实现地震波形的分析预警

 loadplugin("./plugins/mseed/pluginmseed.txt");

  //创建储实时数据的分布式数据库和分区表
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)

  //删除流表及流数据引擎
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)

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)
 */

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)

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])
}	

def streamsimulate(meta,data,taginfo){
	/*
	 * description：
	 * 	此函数用于接收 miniseed 实时流文件，并将其解析为结构化数据，存入流数据表中   
	 * input：
	 * 	meta：stream table   存放 meta 数据的流表   
	 * 	data：stream table  存放采样数据的流表  
	 * 	taginfo：dict   tagid 和 id 映射的字典     
	 */
	ret = mseed::parsestream(file("../streamminiseed/zj_a0001_40_e_i_e.20000101.mseed").readbytes(30000000))
	data_num_row = ret[`data].size()
	meta_num_row = ret[`metadata].size()
	do{
		ret = mseed::parsestream(file("../streamminiseed/zj_a0001_40_e_i_e.20000101.mseed").readbytes(30000000))
		if(ret[`data].size() > data_num_row){
			now_data_num_row = ret[`data].size()
			now_meta_num_row = ret[`metadata].size()
			qdata = ret[`data][data_num_row:now_data_num_row]
			qmeta= ret[`metadata][meta_num_row:now_meta_num_row]
			data_num_row = now_data_num_row
			meta_num_row = now_meta_num_row
			tmp = taginfo[qmeta.id]
			qdata.replacecolumn!(`id,taginfo[qdata.id])
			qmeta = select tmp as ids ,* from qmeta
			delete from qdata where id = null
			delete from qmeta where ids = null
			objbyname(data).append!(qdata)
			objbyname(meta).append!(qmeta)
		}
	}while(true)
}

meta,data= "metastream","datastream"
d = exec tagid,id from loadtable("dfs://real","taginfo")
taginfo=dict(d[`tagid],d[`id])
submitjob("streamsimulate","streamsimulate",streamsimulate,meta,data,taginfo)

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=””])

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)

createanomalydetectionengine(name,metrics,dummytable,outputtable,timecolumn,[keycolumn],[windowsize], 
[step],[garbagesize],[roundtime=true],[snapshotdir],[snapshotintervalinmsgcount],[raftgroup])

[max(data)>= 300,min(data)<=-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)	

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并行加速
}

//波形实时流数据可视化展示
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

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)