筛选指定特征的曲线段一——计算曲线特征
筛选指定特征的曲线段一——计算曲线特征
算法背景
很多仪表设备都能定时产生一些数据，称之为时间序列。将时间序列画在图上就是一条曲线，如下图：
图中横轴是时间(dd hh:mm)，纵轴是数值。
时间序列走势在一定程度上反映实际生产状况，称为工况。如曲线走势平稳表示正常，快速上升或快速下降表示可能有异常情况发生。将曲线走势与工况对应，有助于分析生产活动中存在的问题，提高生产效率。工况与走势对应情况如下图：
准确识别曲线的走势非常重要，本算法根据历史数据识别曲线走势，作为进一步工况分析的基础。
算法思路
曲线的趋势特征有很多，如果穷举所有的走势是不可能的。使用一些原子化的特征来描述这些特征就把不可能变为了可能，就像积木，拿着几类特定形状的积木可以搭成千千万万种形状。
我们设计下面三类原子特征，并为每个特征设计了两种指数用于表征程度：
1. 升降特征：描述曲线升降及升降快慢的特征
（1） 升降指数：描述曲线升降状态的特征
（2） 升降速度指数：描述曲线升降快慢的特征
2. 振幅特征：描述曲线振动幅度变化情况的特征
（1） 振幅指数：描述曲线振动幅度的特征
（2） 振幅升降指数：描述振动幅度升降的特征
3. 振频特征：描述曲线振动频率变化情况的特征
（1） 振频指数：描述曲线振动频率的特征
（2） 振频升降指数：描述振动频率升降的特征
有了这三类特征中的六个指数，就可以组合出各种走势的曲线段，如：
曲线上升段=升降指数大
振荡发散=振幅升降指数大+振频升降指数下降不多
先上升后平稳=上升+平稳=升降指数大+升降指数不变
之前介绍了描述曲线走势的方法——《拟合曲线趋势》
利用主线来量化六个指数，最后经过指数组合就可以筛选出指定特征的曲线段。
篇幅原因，本文先介绍六个指数特征的计算方法。
算法过程
时间序列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]
Mi是1个lift_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_fi是1个range_itv长度的序列。
其中range_f’j是Range_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时，振频下降。
