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

算法背景

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

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

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

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

算法思路

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

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

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

（1） 升降指数：描述曲线升降状态的特征

（2） 升降速度指数：描述曲线升降快慢的特征

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

（1） 振幅指数：描述曲线振动幅度的特征

（2） 振幅升降指数：描述振动幅度升降的特征

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

（1） 振频指数：描述曲线振动频率的特征

（2） 振频升降指数：描述振动频率升降的特征

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

曲线上升段=升降指数大

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

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

之前介绍了描述曲线走势的方法——《拟合曲线趋势》

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

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

算法过程

时间序列X：

X=[x₁,x₂,…,x_n]

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

主线用M表示：

M=fit_main(X,k)

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

假设拟合的主线M

M=[m₁,m₂,…,m_n]

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

1. 升降特征

（1） 升降指数

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

lift_i=m_{i+lift_itv\2}-m_{i-(lift_itv-lift_itv\2)}

其中m_{i+lift_itv\2}是第i个点之后的第lift_itvl\2个点，m_{i-(lift_itv-lift_itv\2)}是第i个点之前的第lift_itv-lift_itv\2个点，两者的索引相差lift_itv。

（2） 升降速度指数

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

计算lift_s_i的序列为M_i

M_i=[m_{i-(lift_itv-lift_itv\2)}, m_{i-(lift_itv-lift_itv\2)}+1,…, m_{i+lift_itv\2}]

M_i是1个lift_itv长度的序列。

M_i中元素的索引序列I_i

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

lift_s_i=linefit(I_i,M_i)

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

2. 振幅特征

（1） 振幅指数

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

Range_f=|X-M|

计算range_i的振幅因子序列为Range_f_i

Range_f_i=[range_f_{i-(range_itv-range_itv\2)}, range_f_{i-(range_itv-range_itv\2)}+1,…, range_f_{i+range_itv\2}]

Range_f_i是1个range_itv长度的序列。

其中range_f’_j是Range_f_i中较大的前一半因子的成员。

（2） 振幅升降指数

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

RangM=fit_main(Range,k)

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

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

Range_Lift=F_Lift(Range, range_lift_itv)

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

3. 振频特征

（1） 振频指数

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

fre_i=count(range_f_i*range_f_i-1<0)

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

（2） 振频升降指数

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

Fre_Lift=F_{Range_lift}(Fre, fre_lift_itv)

其中F_{Range_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时，振频下降。