筛选指定特征的曲线段一——计算曲线特征

算法背景

很多仪表设备都能定时产生一些数据,称之为时间序列。将时间序列画在图上就是一条曲线,如下图:

..

图中横轴是时间(dd hh:mm),纵轴是数值。

时间序列走势在一定程度上反应实际生产状况,称为工况。如曲线走势平稳表示正常,快速上升或快速下降表示可能有异常情况发生。将曲线走势与工况对应,有助于分析生产活动中存在的问题,提高生产效率。工况与走势对应情况如下图:

..

准确识别曲线的走势非常重要,本算法根据历史数据识别曲线走势,作为进一步工况分析的基础。

算法思路

曲线的趋势特征有很多,如果穷举所有的走势是不可能的。使用一些原子化的特征来描述这些特征就把不可能变为了可能,就像积木,拿着几类特定形状的积木可以搭成千千万万种形状。

我们设计下面三类原子特征,并为每个特征设计了两种指数用于表征程度:

1.        升降特征:描述曲线升降及升降快慢的特征

(1)      升降指数:描述曲线升降状态的特征

(2)      升降速度指数:描述曲线升降快慢的特征

2.        振幅特征:描述曲线振动幅度变化情况的特征

(1)      振幅指数:描述曲线振动幅度的特征

(2)      振幅升降指数:描述振动幅度升降的特征

3.        振频特征:描述曲线振动频率变化情况的特征

(1)      振频指数:描述曲线振动频率的特征

(2)      振频升降指数:描述振动频率升降的特征

 

有了这三类特征中的六个指数,就可以组合出各种走势的曲线段,如:

曲线上升段=升降指数大

振荡发散=振幅升降指数大+振频升降指数下降不多

先上升后平稳=上升+平稳=升降指数大+升降指数不变

之前介绍了描述曲线走势的方法——k值拟合主线法,详情见链接:

利用主线来量化六个指数,最后经过指数组合就可以筛选出指定特征的曲线段。

篇幅原因,本文先介绍六个指数特征的计算方法。

算法过程

时间序列X

X=[x1,x2,…,xn]

其中xi是序列X的第i个元素。

主线用M表示:

M=fit_main(X,k)

其中fit_main()是主线拟合函数。

假设拟合的主线M

M=[m1,m2,…,mn]

其中mi是序列M的第i个元素。

1.      升降特征

(1)      升降指数

主线一段区间内的差称为升降指数,记为Lift。这段区间称为升降区间,记为lift_itv

lifti=mi+lift_itv\2-mi-(lift_itv-lift_itv\2)

其中mi+lift_itv\2是第i个点之后的第lift_itvl\2个点,mi-(lift_itv-lift_itv\2)是第i个点之前的第lift_itv-lift_itv\2个点,两者的索引相差lift_itv

(2)      升降速度指数

主线一段区间内的斜率称为升降速度指数,即为Lift_s。这段区间称为升降速度区间,记为lift_s_itv

计算lift_si的序列为Mi

Mi=[mi-(lift_itv-lift_itv\2), mi-(lift_itv-lift_itv\2)+1,…, mi+lift_itv\2]

Mi1lift_itv长度的序列。

Mi中元素的索引序列Ii

Ii=[1,2,…,lift_itv+1]

lift_si=linefit(Ii,Mi)

其中linefit()是最小二乘法线性拟合函数(这里只取斜率,不关系截距)。

2.      振幅特征

(1)      振幅指数

一段区间内原值偏离主线的多少称为振幅因子,记为Range_f,振幅因子中较大的一半因子的平均值称为振幅,记为Range。这段区间称为振幅区间,记为range_itv

Range_f=|X-M|

计算rangei的振幅因子序列为Range_fi

Range_fi=[range_fi-(range_itv-range_itv\2), range_fi-(range_itv-range_itv\2)+1,…, range_fi+range_itv\2]

Range_fi1range_itv长度的序列。

..

其中range_f’jRange_fi中较大的前一半因子的成员。

(2)      振幅升降指数

振幅序列比较杂乱无章,有点类似于开始的原值序列X,所以先利用拟合主线的方法平滑一下振幅序列,得到振幅主线,记为RangeM

RangM=fit_main(Range,k)

振幅主线一段区间内的差称为振幅升降指数,记为Range_Lift

振幅升降指数的计算方法和升降指数的方法相同,只是输入的序列从主线M变成振幅Range,升降区间lift_itv变成振幅升降区间range_lift_itv

Range_Lift=FLift(Range, range_lift_itv)

其中Flift(…)是升降指数的计算方法。

3.      振频特征

(1)      振频指数

一段区间内原值穿越主线的次数称为振频,记为Fre,这段区间称为振频区间,记为fre_itv。原值穿越主线的次数等于振幅因子相邻两个元素之间异号的次数。

frei=count(range_fi*range_fi-1<0)

其中count(…)是计数函数。

(2)      振频升降指数

振频升降指数的计算方法和振幅升降指数相同,只是输入的序列从振幅Range变为振频Fre,振幅升降区间range_lift_itv变成振频升降区间fre_lift_itv

Fre_Lift=FRange_lift(Fre, fre_lift_itv)

其中FRange_lift (…)是振幅升降指数的计算方法。

写出代码:

clac_feature.dfx


A

B

C

D

1

=X

/时间序列X



2

=itv




3

=ft_name




4

=func(A7,X,itv,ft_name)



5

return A4




6

/计算指数,参数:序列,区间,指数名


7

func




8


=A7.len()



9


if B8<=1

return 0


10


=B7\2



11


=to(if(B8<=B7+1,B8,B7+1))

12


=B7-B10



13


=to(if(B8>B7,B8-B7,1),B8)

14


if   C7=="lift_speed"

=A7.((s=if(#<=B7-B10,A7(B11),if(#>B8-B10,A7(B13),~[B10-B7:B10])),func(A23,s)))

/升降速度指数

15



return C14

16


else if   C7=="range"

=A7.((s=if(#<=B7-B10,A7(B11),if(#>B8-B10,A7(B13),~[B10-B7:B10])),func(A30,s)))

/幅度指数

17



return C16

18


else if   C7=="fre"

=A7.((s=if(#<=B7-B10,A7(B11),if(#>B8-B10,A7(B13),~[B10-B7:B10])),func(A35,s)))

/频率指数

19



return C18

20


else if   C7=="lift"

=A7.(if(#<=B7-B10,A7(B11.m(-1))-A7(B11.~),if(#>B8-B10,A7(B13.m(-1))-A7(B13(1)),~[B10]-~[-B12])))

/升降指数

21



return C20

22

/升降速度,参数:序列



23

func




24


=A23.len()


25


if B24<=1

return 0


26


=to(B24)



27


if   A23.id().len()==1

return 0


28


=(B24*B26.sum(~*A23(#))-B26.sum()*A23.sum())/(B24*B26.sum(~*~)-power(B26.sum(),2))

/升降速度指数(斜率)

29

/计算幅度,参数:序列



30

func




31


=A30.len()\2


32


if B31<=1

return A30.max()

33


=A30.top(-B31,abs(~)).avg()

34

/振频,参数:振幅序列



35

func




36


1e-5



37


=A35.len()


38


if B37<=1

return 0


39


=A35.count(if(#==1,,~*~[-1]<0))







计算各种特征

call_calc_feature.dfx


A

B

1

=X

/时间序列X

2

=k

/主线参数k

3

=itv

/特征区间

4

=fit_main(A1,A2)

/主线

5

=clac_feature(A4,A3,"lift")

/升降指数

6

=clac_feature(A4,A3,"lift_speed")

/升降速度指数

7

=A1--A4

/振幅因子

8

=clac_feature(A7,A3,"range")

/振幅

9

=fit_main(A8,A2)

/振幅主线

10

=clac_feature(A9,A3,"lift")

/振幅升降指数

11

=clac_feature(A7,A3,"fre")

/振频

12

=fit_main(A11,A2)

/振频主线

13

=clac_feature(A12,A3,"lift")

/振频升降指数

参数说明:

clac_feature.dfx是计算特征的脚本,需要计算特征时,调用该脚本即可。

X:输入时间序列X

itv:计算特征的区间,如升降区间,振幅区间等;

ft_name:特征名字,如lift,lift_speed,range,range_lift,fre,fre_lift

call_calc_feature.dfx是调用clac_feature.dfx的脚本,目的是展示如何计算各种特征。

X:输入时间序列X

k:主线拟合参数k

itv:特征区间(这里只设置了一个特征区间,所以所有的特征都用这个区间来算,但实际上计算各个特征时可以使用不同的特征区间)。

应用举例

时间序列X如下图:

..

图中横轴是时间(dd hh:mm),纵轴是序列中的值。

升降指数图:

..

图中横轴是时间(dd hh:mm),左纵轴是时间序列的值,右纵轴是升降指数的值。升降指数大于0时,主线上升,时间序列呈上升趋势;升降指数小于0时,主线下降,时间序列呈下降趋势;

图例中X是时间序列,Main是主线,Lift是升降指数。

升降速度指数图:

..

图例中Lift_speed是升降指数。

升降速度指数大于0时,主线上升,时间序列呈上升趋势,且升降速度指数越大,上升越快;升降速度指数小于0时,主线下降,时间序列呈下降趋势,且升降速度指数越小,下降越快;振幅指数图:

..

图例中Range是振幅指数。

振幅指数越大,曲线围绕主线振动的幅度越大。

振幅升降指数图:

..

图例中Range是振幅指数,Range_main是振幅主线,Range_lift是振幅升降指数。

振幅升降指数大于0时,振幅上升,曲线成震荡发散趋势,振幅升降指数小于0时,振幅下降,曲线成震荡收缩趋势。

振频指数图:

..

图例中Fre是振幅指数。

振频指数越大,原值穿越主线的次数越多。

振幅升降指数图:

..

图例中Fre是振幅指数,Fre_main是振幅主线,Fre_lift是振幅升降指数。

振频升降指数大于0时,振频上升,振频升降指数小于0时,振频下降。


以下是广告时间

对润乾产品感兴趣的小伙伴,一定要知道软件还能这样卖哟性价比还不过瘾? 欢迎加入好多乾计划。
这里可以低价购买软件产品,让已经亲民的价格更加便宜!
这里可以销售产品获取佣金,赚满钱包成为土豪不再是梦!
这里还可以推荐分享抢红包,每次都是好几块钱的巨款哟!
来吧,现在就加入,拿起手机扫码,开始乾包之旅



嗯,还不太了解好多乾?
猛戳这里
玩转好多乾