4. 环境准备
本节主要内容为加载插件、创建流数据表等,是后续波形实时数据模拟、异常检测的前备工作。
4.1 加载插件
为实现波形数据异常检测,DolphinDB 开发了对应的插件,见附录(目前只提供离线版,插件正式发布后会提供下载链接)。下载插件,在 plugins 目录下建立 filterpicker、rtseis、tf 三个文件夹,将插件解压到对应的文件夹里。
注意:添加环境变量前需关闭 DolphinDB Server。
Linux 添加环境变量:
export LD_LIBRARY_PATH=<YOUR DolphinDB Dir>/server/plugins/rtseis/:$LD_LIBRARY_PATH
export LD_LIBRARY_PATH=<YOUR DolphinDB Dir>/server/plugins/tf/:$LD_LIBRARY_PATH
启动 DolphinDB Server,通过 loadPlugin()
函数加载插件
try{ loadPlugin("./plugins/filterpicker/PluginFilterPicker.txt") }catch(ex){print(ex) }
try{ loadPlugin("./plugins/rtseis/PluginRTSeis.txt") }catch(ex){print(ex) }
try{ loadPlugin("./plugins/tf/PluginTf.txt") }catch(ex){print(ex) }
4.2 流数据表环境清理
清理流数据表之前需要取消订阅,通过 unsubscribeTable()
函数取消订阅,然后再通过 dropStream()
函数删除流数据表:
unsubscribeTable(tableName = "dataStream", actionName="filterPickerPredict");go
unsubscribeTable(tableName = "pickerStream",actionName="tensorFlowPredict");go
try{dropStreamTable("dataStream")}catch(ex){print(ex)}
try{dropStreamTable("pickerStream")}catch(ex){print(ex)}
try{dropStreamTable("tensorStream")}catch(ex){print(ex)}
以 enableTableShareAndPersistence()
函数创建流数据共享表并持久化:
//dataStream,接收波形实时流数据
enableTableShareAndPersistence(table=streamTable(1000000:0, `tagid`ts`data, [INT,TIMESTAMP, INT]), tableName=`dataStream, cacheSize=1000000);
//pickerStream
enableTableShareAndPersistence(table=streamTable(100000:0, `ts`id, [TIMESTAMP,INT]), tableName=`pickerStream, cacheSize=1000000);
//tensorStream
enableTableShareAndPersistence(table=streamTable(100000:0, `ts`id, [TIMESTAMP, INT]), tableName=`tensorStream, cacheSize=1000000);
创建分布式数据库和维度表:
if(existsDatabase("dfs://seis01")) dropDatabase("dfs://seis01");
//创建分布式数据库
create database "dfs://seis01" partitioned by VALUE(1..10), engine='TSDB'
//创建维度表
create table "dfs://seis01"."tagInfo"(
xfdsn SYMBOL,
net SYMBOL,
sta SYMBOL,
loc SYMBOL,
chn SYMBOL,
id INT[compress="delta"]
)
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] )
}
xfdsn = "XFDSN:"+netList+"_"+staList+"_"+locList+"_"+colt
t = table(xfdsn,netList as net,staList as sta,locList as loc,chn,1..150 as id)
loadTable( "dfs://seis01","tagInfo").append!(t);
5. 方案实现
本章包括以下内容
- 波形实时数据模拟
- 一级异常检测
- 二级异常检测
5.1 波形实时数据模拟
实际生产环境中每个通道每10毫秒采集一次数据,每个台站有三个通道。以下代码模拟了三个台站的波形实时数据:
/*
* 模拟实时数据
*/
def insertIntoDataStream(){
do{
tagidList = 1 2 3 4 5 6 7 8 9
ts = now()+(0..400)*10
t = table(stretch(tagidList,3600)as tagid,take(ts,3600) as ts,randNormal(-2500, 1000, 3600) as data)
objByName(`dataStream).append!(t)
sleep(3500)
}while(true)
}
jobId = submitJob("simulate","simulate",insertIntoDataStream);
//通过以下方式取消Job
//cancelJob(jobId)
5.2 一级异常检测
通过 subscribe()
函数订阅实时数据,进行一级异常检测,处理逻辑及代码如下所示:
**step1:**将实时数据按照 tagid 分组
**step2:**对每一 tagid,调用 filterPicker::createFilterPicker()
函数创建 filter handler
**step3:**调用 filterPicker::filerPicker(pickerPtr, timeCol, dataCol, [fixedSize])
函数进行模板匹配,得到突峭点。其中,各个参数含义如下所示:
- pickerPtr:每一 tagid 对应的 filter handler;
- timeCol:每一 tagid 对应的时间列;
- dataCol:每一 tagid 对应的采样值;
- fixedSize:批计算数据量大小,若数据不足 fixedSize,会累积到下次凑足 fixedSize 再进行计算
其中,突峭点算法参考 ALomax - FilterPicker - a general purpose, broad-band, phase detector and picker
/*
*一级异常检测处理函数
*计算波形实时数据中的突峭时间点
*返回值为突峭时间点和通道id
*/
def pickerHandler(mutable pickerSt, mutable res, vz, msg){
tags=groups(msg[`tagid]) //按tagid分组
for (tag in tags.keys()){
re=res[tag]
if(re==NULL){
//创建filter handler
re= filterPicker::createFilterPicker()
res[tag]=re
}
//将数据传入filterPicker进行计算
vt=filterPicker::filterPicker(re,msg[tags[tag]][`ts], msg[tags[tag]][`data].float(), 1000)
if(vt.size()>0){
vid=take(tag,vt.size())
pickerSt.append!(table(vt as ts,vid as id))
}
}
}
vz=exec id from loadTable("dfs://seis01","tagInfo") where chn ilike("%z") //只需要计算Z分量的数据
res=dict(INT, ANY)
//订阅dataStream,进行一级异常检测,异常检测结果输出到filterStream中
subscribeTable(tableName = "dataStream", actionName="filterPickerPredict", offset=-1,handler=pickerHandler{objByName("pickerStream"), res, vz, }, msgAsTable=true, batchSize = 2000, throttle=1,reconnect=true)
一级异常检测结果样例如下所示:
ts | id |
---|---|
2023.04.20T14:37:57.838 | 1 |
上述异常时间点前后几秒内的采样数据如下:
tagid | ts | data |
---|---|---|
… | … | … |
1 | 2023.04.20T14:37:57.797 | -3120 |
1 | 2023.04.20T14:37:57.807 | -3005 |
1 | 2023.04.20T14:37:57.818 | -521 |
1 | 2023.04.20T14:37:57.828 | -3886 |
1 | 2023.04.20T14:37:57.838 | 1441 |
1 | 2023.04.20T14:37:57.848 | 1220 |
1 | 2023.04.20T14:37:57.858 | -974 |
1 | 2023.04.20T14:37:57.868 | -396 |
1 | 2023.04.20T14:37:57.878 | -943 |
… | … | … |
从采样数据来看,“2023.04.20T14:37:57.838”时间点前后振幅明显变大,确实是一个突峭点。
5.3 二级异常检测
通过 subscribe()
函数订阅 pickerStream,进行去趋势、滤波、归一化以及二级异常检测,处理逻辑及代码如下所示:
**step1:**调用 tf::load(modelFileName, inputName, outputName, [shape])
函数导入训练好的 TensorFlow 模型(本文使用的深度学习模型不方便对外提供,读者需自行训练模型并调用)。其中各参数含义如下:
- modelFileName:模型路径。如果是 keras 生成的 h5 格式的模型文件,可用 Keras to TensorFlow 将 h5 格式的模型文件转成 TensorFlow 的 pb 文件后再进行加载。
- inputName:模型中输入节点名称,如果不知道,可用附录的 printTensorName.py 文件进行打印,输出的第一个即为输入节点名称。
- outputName:模型中输出节点名称,如果不知道,可用附录的 printTensorName.py 文件进行打印,输出的最后一个即为输出节点名称。
- shape:N 维 (N>=3) 数组的 shape,当 N>=3 时才需要提供,类型为 int 类型的 vector,如[3, 5, 8]表示 shape 为3×5×8。如果提供了 shape,后面用于预测的 data 请使用
flatten(x)
函数转成的一维向量。
**step2:**提取突峭点前后两秒内三分量数据
**step3:**剔除趋势成分(本文将趋势视为线性的,采用最小二乘法估计线性函数的参数)
**step4:**调用 rtseis::sosBandPass(x, nTaps, f_min, f_max, windowType, fsamp, ripple)
函数对三通道数据进行 sosBandPass 变换。其中,各个参数含义如下所示:
- x:需要处理的数据
- nTaps:筛选器参数
- f_min:最低频率
- f_max:最高频率
- windowType:窗口类型,0-3分别代表 BUTTERWORTH,BESSEL,CHEBYSHEV1,CHEBYSHEV2
- fsamp:采样频率
- ripple:windowType 为2或3时的波纹大小
**step5:**先对滤波数据进行归一化处理,然后调用 tf::predict(model, data)
函数进行预测,当预测概率 prob_P 或者 prob_s 的概率大于0.9时,认为该段信号很可能是地震波 P 波或者 S 波信号,然后输出对应的 filterpicker 的突峭点。其中,tf::predict(model, data)
各参数含义如下所示:
- model:tensorflow 的 model,为
tf::load()
函数的返回值 - data:需要预测的数据,仅支持同种数据类型的 matrix、vector、table 或者 dictionary。如果在
tf::load()
时提供了 shape 参数,则 data 必须为使用flatten(x)
函数转成的一维向量。
def tensorHandler(mutable tensorSt, mutable infoTable, mutable model, mutable cache, msg){
if(msg.type()!=0) cache.append!(msg)
if(cache.size()==0) return
currTime = now() - 2000
outed = select * from cache where timestamp < currTime
for(s in outed){
//获取同个台站id信息
info = select net,sta,loc from infoTable where id == s[`id]
netInfo = info[`net][0]
staInfo = info[`sta][0]
bhzInfo = info[`loc][0]
bheId = exec id from infoTable where net == netInfo and sta == staInfo and loc == bhzInfo and ilike(chn,"%e")
bhnId = exec id from infoTable where net == netInfo and sta == staInfo and loc == bhzInfo and ilike(chn,"%n")
bhzId = exec id from infoTable where net == netInfo and sta == staInfo and loc == bhzInfo and ilike(chn,"%z")
//获取三个分量数据
begT = timestamp(s[`timestamp] - 2000)
endT = timestamp(s[`timestamp] + 1999)
ve = exec data from dataStream where tagid ==bheId[0] and ts between(begT:endT)
vn = exec data from dataStream where tagid ==bhnId[0] and ts between(begT:endT)
vz = exec data from dataStream where tagid ==bhzId[0] and ts between(begT:endT)
if(ve.size()== 400 && vn.size() == 400 && vz.size() == 400) {
//去趋势
tem = ols(vn, 0..399)
bhn = vn - (tem[0] + tem[1]*0..399)
tem = ols(ve, 0..399)
bhe = ve - (tem[0] + tem[1]*0..399)
tem = ols(vz, 0..399)
bhz = vz - (tem[0] + tem[1]*0..399)
//滤波
bhn = rtseis::sosBandPass(bhn,4,3.0,20.0,0,100.0,0.0);
bhe = rtseis::sosBandPass(bhe,4,3.0,20.0,0,100.0,0.0);
bhz = rtseis::sosBandPass(bhz,4,3.0,20.0,0,100.0,0.0);
m = matrix(bhn,bhe,bhz)
//矩阵数据归一化
m1 = m / max(max(abs(m)))
//输入模型进行预测, ts返回值为预测结果
ts = tf::predict(model, float(m1))
//prob_p 是P波的概率,prob_S是S波概率, prob_N是噪声的概率。
prob_P = ts[0]
prob_S = ts[1]
prob_N = ts[2]
if(prob_P > 0.90 || prob_S > 0.90){
tensorSt.append!(table(s[`timestamp] as timestamp, s[`id] as id))
}
}
}
delete from cache where timestamp < currTime
}
//导入模型
model = tf::load("./model_pol.pb", "conv1d_1_input", "activation_7_1/concat")
infoTable = select * from loadTable("dfs://seis01","tagInfo")
vz=exec id from loadTable("dfs://seis01","tagInfo") where chn ilike("%z") //只需要计算Z分量的数据
cache = table(10:0, `timestamp`id, [TIMESTAMP, INT])
//订阅pickerStream,进行二级异常检测
subscribeTable(tableName = "pickerStream",actionName="tensorFlowPredict", offset=-1, handler=tensorHandler{objByName("tensorStream"),infoTable, model,cache,}, msgAsTable=true, batchSize = 1000, throttle=0.4, timeTrigger=true,reconnect=true)
对一级异常检测的结果进行二级异常检测后,结果样例如下表所示:
ts | id |
---|---|
2023.04.20T14:37:57.838 | 1 |
6. 附录
插件及脚本文件
自我介绍一下,小编13年上海交大毕业,曾经在小公司待过,也去过华为、OPPO等大厂,18年进入阿里一直到现在。
深知大多数嵌入式工程师,想要提升技能,往往是自己摸索成长或者是报班学习,但对于培训机构动则几千的学费,着实压力不小。自己不成体系的自学效果低效又漫长,而且极易碰到天花板技术停滞不前!
因此收集整理了一份《2024年嵌入式&物联网开发全套学习资料》,初衷也很简单,就是希望能够帮助到想自学提升又不知道该从何学起的朋友,同时减轻大家的负担。
既有适合小白学习的零基础资料,也有适合3年以上经验的小伙伴深入学习提升的进阶课程,基本涵盖了95%以上嵌入式&物联网开发知识点,真正体系化!
由于文件比较大,这里只是将部分目录大纲截图出来,每个节点里面都包含大厂面经、学习笔记、源码讲义、实战项目、讲解视频,并且后续会持续更新
如果你觉得这些内容对你有帮助,可以+V:Vip1104z获取!!! (备注:嵌入式)
![img](https://img-blog.csdnimg.cn/img_convert/334139e6fb15c74b4754b69204f3a966.jpeg)
最后
资料整理不易,觉得有帮助的朋友可以帮忙点赞分享支持一下小编~
你的支持,我的动力;祝各位前程似锦,offer不断,步步高升!!!
.(img-orNAz2K9-1712238797809)]
由于文件比较大,这里只是将部分目录大纲截图出来,每个节点里面都包含大厂面经、学习笔记、源码讲义、实战项目、讲解视频,并且后续会持续更新
如果你觉得这些内容对你有帮助,可以+V:Vip1104z获取!!! (备注:嵌入式)
![img](https://img-blog.csdnimg.cn/img_convert/334139e6fb15c74b4754b69204f3a966.jpeg)
最后
资料整理不易,觉得有帮助的朋友可以帮忙点赞分享支持一下小编~
你的支持,我的动力;祝各位前程似锦,offer不断,步步高升!!!