Webshell动态沙箱日志信号处理（signal process）和时序分析（time series analysis）问题初探
1. 时序问题简介
时序数据（time series data）广泛存在于现实生活中，是指同一种现象在不同时间上的相继观察值排列而成的一组数字序列，其时间轴上的采样值通常又被称为特征。
由于时序数据与时间相关联，因而其数据量一般都是非常庞大的，这就对时序数据挖掘技术提出了更高的要求。
2. 信号和时间序列的来源
信号并但是计算机领域的专有名词，地球这颗行星以及她周围的空间都是信号的来源，例如：
- 测量太阳黑子的数量
- 不同地区温度的变化
- 风速
- 小行星的速度
- 股票价格及相关衍生品价格走势
- 大小企业的销售额
- 制造业，互联网相关活动
- 能源生产变化趋势
- 政治和社会事件发生趋势
- 大脑活动（EGG）生物信号
- 心脏活动（ECG）
- 肌肉张力（EMG）
- 从可穿戴设备上收集的数据比如脉搏、基于加速度计测量的活动、睡眠、压力指标
Relevant Link:
http://www.python88.com/topic/27573
3. 经典时间序列和信号分析方法
学界和工业界针对时间序列和信号分析领域创造出一系列数学模型和方法。本章对一些最重要模型进行概要介绍。
0x1：时域分析
时域分析，主要是关于“看”时间序列在一段时间内如何演变的分析方法。它包括时序的长度宽度、统计特征以及其他“可见的”特征。
时域信号处理中，时域信号压缩是一类比较常见的应用。
0x2：频域分析
很多信号很难用时间变化来描述，但是可以用幅度及其变化来表示。傅里叶分析和小波分析是这类方法的代表。
def fourier_transform(): import numpy as np import matplotlib.pyplot as plt def signal(sinuses): ''' This function creates a function f(x) = sin(a*pi*2)+sin(b*p*2)+... plots the graph of f(x) and saves it Args: sinuses(Float) : The list of constants a,b,c,... Returns: t (numpy array): X axis s (numpy array): Y axis ''' # Create the singal as a sum of different sinuses t = np.linspace(0, 0.5, 800) s = 0 for i in range(len(sinuses)): s += np.sin(sinuses[i] * 2 * np.pi * t) # Plot the signal plt.style.use('seaborn') fig = plt.figure(figsize=(8, 4)) ax = fig.add_subplot(1, 1, 1) ax.plot(t, s, label=r'$y=f(x)$') ax.set_title(" Signal ", fontsize=20) ax.set_ylabel("Amplitude") ax.set_xlabel("Time [s]") ax.legend(loc='best') ax.grid(True) plt.show() fig.savefig('signal.png') return s, t def Fourier(s, t, alg="False"): '''This function performs the FFT in the function f(x) defined in signal function, plots and saves the figure Args: t (numpy array): X axis s (numpy array): Y axis alg (string): Variable to determine whether to use my FFT or numpy's ''' # Perform the Fourier Transform if alg == "True": fft = FFT(s) else: fft = np.fft.fft(s) T = t[1] - t[0] # sample rate N = s.size # 1/T = frequency f = np.linspace(0, 1 / T, N) # Plot the signal Decomposition plt.style.use('seaborn') fig = plt.figure(figsize=(8, 4)) ax = fig.add_subplot(1, 1, 1) ax.set_title(" Decomposed Signal ", fontsize=20) ax.set_ylabel("Amplitude") ax.set_xlabel("Frequency [Hz]") print "x: ", f[:N // 2] print "y: ", np.abs(fft)[:N // 2] * 1 / N ax.bar(f[:N // 2], np.abs(fft)[:N // 2] * 1 / N, width=1.5) # 1 / N is a normalization factor ax.grid(False) plt.show() fig.savefig("Decomposed_signal.png") def FFT(x): """Compute the discrete Fourier Transform of the signal x Args: x (numpy array): The function f(x) Returns: M (numpy array): The FFT """ x = np.asarray(x, dtype=float) N = x.shape[0] n = np.arange(N) k = n.reshape((N, 1)) M = np.exp(-2j * np.pi * k * n / N) return np.dot(M, x) sinuses = [50, 20] alg = True s, t = signal(sinuses) # Create, plot and save the signal Fourier(s, t, alg) # Decompose it, plot it and save it if __name__ == '__main__': fourier_transform()
Relevant Link:
https://github.com/fotisk07/Fourier-Transform
0x3：近邻分析
有时候我们仅仅需要对比两个信号或者测量这两个信号之间的距离，我们不能使用常规的距离矩阵比如欧几里得矩阵，因为信号的长度和相似度会同时改变。
动态时间扭曲（dynamic time warping）是运用在时间序列方面的典型矩阵技术。
0x4：(S)AR(I)MA(X)模型
这是一类非常流行的数学模型，它根据时间序列中线性自相关来解释将来的波动。
The equation for a SARMA(p,q)(P,Q) model
0x5：时间信号分解
一个重要的预测方法是将时间序列分解成几个逻辑部分：趋势部分、季节部分、以及剩余部分。这些被被分解的部分可以通过加总和乘法还原成原来的序列。
STL（Seasonal and Trend decomposition using Loess）就是这方面的一个典型算法。
0x6：非线性动力学
非线性动力学，是指将信号或者时间序列当成动态系统，用微分方程（常规微分方程，偏微分方程，随机微分方程以及其他微分方程）作为研究工具来研究他们。
0x7：机器学习/深度学习
以RNN、LSTM为代表的深度学习技术，作为一种专为序列分析开发的神经网络，它能保存潜在的模式并学习时间依赖关系，是完全的图灵机并能对付任何长度的序列。
0x8：异常检测
另一个与序列数据（常常是流数据）相关的重要的和流行的任务是异常检测，即发现数据中“非在此刻希望看到的”情况。
通常我们用设置阈值方法和测量距离（有时是随机距离）来处理这类问题。
Relevant Link:
https://www.zhihu.com/question/304796104/answer/547001050 https://zhuanlan.zhihu.com/p/52089047
4. 时域数据分类
0x1：时序分类问题背景
在时序数据挖掘的研究与应用领域，时序数据分类是重要任务之一。
例如，
- 依据语音信号的波形识别出说话人的性别和年龄
- 依据心电图的时序波形识别出病者所患的病症
- 依据地震波的历史数据，去识别地震的类型
- 依据在机器运转过程中进行故障检测和识别故障类型
- 在客户关系管理中根据某段时间的客户购买信息，识别不同的消费群体等等.
衡量分类技术优劣的核心指标是分类准确率，而提高分类准确率途径有两种：
- 一是改进分类器
- 二是采用特征提取技术（feature extraction）
特征提取是在分类前对数据时间采样值上进行适量的归约，以达到减少数据量同时提高分类准确率（底线是不牺牲分类准确率）的目的。
时间序列除了具有的趋势性、季节性、周期性等一般特征之外，不同的时序数据又存在不同的个别特征。
- 金融数据，普遍具有“高峰厚尾”和“平方序列微弱而持续的自相关”的特点
- 地震波则具有强度随时序延伸而减弱的特点
- 语音信号幅度具有一定的范围，并以零幅和近零幅的概率高，而且长时间的语音信号会有相当多的无信号区间，即所谓的语音寂静区间
- 心电信号则具有很强的周期性，它的主要特征是低电压(0.8~1mV)，小电流(12uA)，重复频率低，每个波段具有各自的频率
针对时序数据的这些特征，所选择的特征提取方法应该能提取出时间序列中具有较好分类能力特征，进行特征提取后的特征矢量能够很好地代表原有的时间序列数据，这样才能取得良好的分类效果。
0x2：针对时序分类问题的经典特征提取方法
本小节针对分类中常用的时序特征提取方法归纳和总结，现阶段特征提取方法主要有四类，分别是：
- 基于基本统计方法的特征提取
- 基于模型的特征提取
- 基于变换的特征提取
- 基于分形理论的特征提取
1、 基于基本统计方法的特征提取
基本统计方法的特征提取，就是提取数据波形统计特征特征，来代表原有的时序数据作为特征矢量。
时域的常见基本统计征有：
- 均值
- 方差
- 极值
- 过零点
- 边界点
- 波段的长短峰值等
频域的基本统计方法有：
- 功率谱
- 功率密度比
- 中值频率
- 平均功率频率等
基于统计特征提取后的特征矢量，可采用线性判别式构造分类器，也可采用神经网络进行分类。它适用于信号波形的统计特征比较明显的时间序列数据，如 EEG 信号、ECG 信号等医学数据。
2、基于模型的特征提取
基于模型的特征提取，是指用模型去刻画时间序列数据，然后提取模型的系数作为特征矢量。
对于平稳时间序列，可用通用 ARMA 模型（自回归滑动平均模型）及其特殊情况的自回归模型、滑动平均模型（MA）或组合-ARMA 模型等来进行拟合。
对于非平稳时间序列，则要先将观测到的时间序列进行差分运算，转化为平稳时间序列，再用适当模型去拟合这个差分序列。
当然，不同的时序数据都会有比较适合它的模型去进行特征提取，例如：
- 对于 ECG 数据，通常采用 AR 模型进行特征提取，而相对应的 MAR 模型则用于多通道的 ECG 数据中。AR 模型是一个线性的，二阶矩平稳模型，比较适合短数据分析，不仅在 ECG 数据中具有优势
- 对于金融时序数据，常用的多元线性回归和 ARMA 模型等都不再适合，针对金融市场的价格波动聚集现象，建议采用 ARCH 模型（自回归条件异方差族计量模型）。而采用 SVAR（结构向量自回归）模型更加适合刻画金融数据中常出现的“高峰厚尾”和“弱自相关”现象
选择合适的模型进行特征提取，产生的模型系数作为特征矢量。
针对这类特征矢量，往往采用简单的判别函数即可达到理想的分类效果。
3、基于变换的特征提取
通过变换的手段, 使适合分类的特性突显出来，也是经常用到的特征提取方法。
变换包括：
- 时频变换：时频变换中具有代表性的方法是：
- 快速傅立叶变换
- 短时傅立叶变换
- 倒谱系数等
- 线性变换：线性变换主要有
- PCA
- ICA
- SVD
- 线性判别式分析
- 要素分析
- 映射等等
1）基于时频变换的特征提取
时频变换是将信号从时域变换到频域的一种手段.
时间序列数据在特征提取中常采用傅立叶变换、倒谱系数等时频变换方法。其中最经典的是傅里叶变换。
傅里叶变换是将时域的信号变换成频域的信号，它是把时域的信号变换成由频率、幅值、和相位的正弦波的组合。
将时序数据进行傅立叶变换，然后选择它的系数作为特征矢量，
- 若选择前面的系数，则代表了信号的低频特性
- 选择较大的系数，则代表了信号的能量特征
在心电信号 QRS 波识别中用 DFT 有效进行降维和提取特征，然后用神经网络进行分类。在股市中应用傅立叶变换与反变换，可以有效去除噪声和进行数据约简。
2）基于线性变换的特征提取
线性变换中有很多特征提取方法颇为经典，特别是 PCA 和小波变换。
PCA 为主成分分析，它与独立成分分析 ICA、K －L 变换、奇异值分解极为相似，有殊途同归的效果，在数据降维中都经常用到。研究表明，利用 PCA 变换可以在信息损失最小的前提下，用较少的分量代替原来的高维数据，达到降维的效果，从而使得处理数据的时问和费用大大降低。
小波变换是将信号分解到不同尺度的线性变换。
基于小波变换的特征提取有多种，如：
- 提取模极大值特征
- 能量特征
- 熵特征
- 以及适应性小波网络的特征提取等
4、基于分形理论的特征提取
分形是指具有: 无限精细、非常不规则、无穷自相似结构和非整数的点集。
在大自然中，海岸线、雪花，云雾这些不规则形体都属于分形，即部分与整体有自相似性。
分形学与混沌学息息相关，成为非线性科学的两大重要组成部分，分形理论真正发展只有十余年，但应用于特征提取已越来越广泛，特别是在时间序列中。
分形理论进行特征提取时，并不是所有信号都适合，而是要看信号在某个尺度下是否具有可分形特征，即不同状态下的分形维数是不同的，这样才具有可分性。
如故障诊断中，将信号分为 N 个状态，每个状态可提取一个分形维数作为特征参数，把这个分维数与学习好的故障分类结果比较，可判断是否有故障发生。
分形理论用于特征提取时，主要是针对非线性信号，是用它的定量分析指标分维数作来特征矢量。
Relevant Link:
http://www.c-s-a.org.cn/ch/reader/create_pdf.aspx?file_no=20121053
5. 时域数据轨迹聚类
轨迹聚类在使用中有着很多挑战和局限性。
- 噪声问题：时间序列需要尽可能高信噪比，我们不希望存在大量的噪声，这会引入较大的误差，因此之前需要做去噪。
- 信息留存问题：降噪过程伴随着很多降维操作，不可避免地导致了信息丢失，但是我们要避免过度丢失信息。
- 错误匹配问题：时间序列很多时候存在错位匹配的情况，这个需要相似性度量算法来归并解决，实际中也要根据场景，额外做些处理
- 调参问题：聚类的方法和参数选择也是一个问题，需要借助大量的领域知识来弥补漏洞。
0x1：时间序列的表示
时间序列的表示并没有什么硬性限制条件，抛开与相似性度量的配合，目的只有一个：
在尽可能保留完整的信息量的前提下，尽量压缩原始信号的长度。
而相似性度量一般都会有一些规范需要遵循，否则定义出来的相似性就失去了物理含义，也无法服务后续的聚类等分析方法。
为了定义的方便，我们把相似性视为距离的倒数，这样相似性的定义就转为距离的定义。
在距离的定义中其中最常见的、也是最基本的就是以下三个条件：
- 两个时间序列的距离是非负的。
- 满足对称性，也即：d(a,b)=d(b,a)。
- 三角不等式，也即：d(a,c) ≤ d(a,b)+d(b,c)
1、欧氏距离度量
我们直接把时间序列的数值作为时间序列的表示，用对应时间点之间的欧式距离之和作为距离（假设两个时间序列的长度一致），那么我们就得到了最简单的定义。看上去一些似乎都很顺利，我们拿到了一个结果，然后就可以去做后面的聚类了。但是在实际的应用中，会面临很多问题。
欧式距离最大的问题就是会被噪声或是离群点所影响。比如以下两对时间序列：
- 第一组是十个时间点、均值为0方差为1的时间序列，包含一个离群点
- 第二组是十个时间点、均值为0方差为0.6的时间序列
我们可以调整离群点的值使得两对时间序列的欧式距离接近。
但是如果在物理意义上，我们期望这两组时间序列的距离是不一致的，这就说明我们的定义是不合理的，或者说这不是我们期望的定义。
因此我们希望我们的距离定义能够具有普适性，以适应不同时间序列特性的区分，比如Minkowski距离。
Minkowski距离可以视为一组距离的定义，可以通过参数p调整距离对时间序列特性的适应，特别的，当p=2的时候就是欧式距离。
- 如果我们希望突出两个时间点之间存在差异，而非差异度，我们可以让p值调小
- 反之，我们希望突出两个时间点之间的差异度，那我们可以让p值调大
- 极端的情况，当p趋近于0，结果等价于统计有几对时间点直接存在差异，没对差异计数1
- 当p趋于无穷大，结果等价于时间点对之间距离对最大值
因此，如果我们要剔除离群点的影响，可以把p值调小，要剔除噪声的影响可以把p值调大。
2、离散化抽样
针对离群点和噪声，我们还可以采用离散化的思路来解决，把一定范围内的距离差视为一类，把大于一定阈值的距离差视为一类，这样就弱化了它们带来的影响。
其实这就是一种符号匹配的思想，将时间序列语义化映射到符号或者字符串上，不仅仅是离散化，也可以直接来表示时间点之间的差值或者变化率等信息。
基于这样的表达方式，我们就不能采用数值型的距离度量来进行相似度的计算了，但是我们可以引入文本相似度的计算逻辑，比如说Hamming距离或者Jaccard距离。
3、信号降维
在进行时序信号分析的时候，我们也会遇到机器学习中最常见的问题：维度灾难。
也即，当时间序列过长时，会导致距离的差异会逐渐接近，信号信息过于分散，从而无法区分。这个时候需要做的和机器学习中一样，对时间序列进行降维。
时间序列的降维思路很多，常见的有两大类：
- 全局拟合或分段拟合，一般是用线性函数或是多项式函数，也可以根据物理公式或者领域经验模型（例如某种编程语言中经常会出现的某种时序信号子串）去拟合
- 频域变换，通过频谱特征来表示时间序列。通常可以去做64、128或256点的FFT，也可以使用小波变换等方法
Relevant Link:
https://zhuanlan.zhihu.com/p/52202274
6. 基于领域经验对时序信号进行降维
尽管有很多经典的时序信号分析方法，但是具体到具体的问题场景，我们还是可以将我们的领域经验融合进来，创造出一些新的时序信号降维方法。
0x1：基于CFG（code flow graph）进行降维
在编程语言文本检测中，我们常常会通过沙箱等手段，从样本中提取出opcode执行序列，这种序列本质上也可以看成一种时序信号。
例如下面这个php源代码：
<?php $map = [ 1 => ')', 2 => '(', 3 => ';', 4 => '$_GET[1]', 5 => "eval", ]; $a = ""; if ($_GET[2] == 5 ){ $a .= $map[5]; } if ($_GET[3] == 2 ){ $a .= $map[2]; } if ($_GET[4] == 4 ){ $a .= $map[4]; } if ($_GET[5] == 1 ){ $a .= $map[1]; } if ($_GET[6] == 3 ){ $a .= $map[3]; }else{ $a=""; } eval($a); ?>
经过zend编译后，得到的op_array序列为：
474848484826265051112B511E2A5051112B511E2A5051112B511E2A5051112B511E2A5051112B511E2A26493E
翻译为zend op_name：
END_INIT_ARRAY,
ZEND_ADD_ARRAY_ELEMENT,
ZEND_ADD_ARRAY_ELEMENT,
ZEND_ADD_ARRAY_ELEMENT,
ZEND_ADD_ARRAY_ELEMENT,
ZEND_ASSIGN,ZEND_
ASSIGN,ZEND_FETCH_R,
ZEND_FETCH_DIM_R,
ZEND_IS_EQUAL,
ZEND_JMPZ,
ZEND_FETCH_DIM_R,
ZEND_ASSIGN_CONCAT,
ZEND_JMP,
ZEND_FETCH_R,
ZEND_FETCH_DIM_R,
ZEND_IS_EQUAL,
ZEND_JMPZ,
ZEND_FETCH_DIM_R,
ZEND_ASSIGN_CONCAT,
ZEND_JMP,
ZEND_FETCH_R,
ZEND_FETCH_DIM_R,
ZEND_IS_EQUAL,
ZEND_JMPZ,
ZEND_FETCH_DIM_R,
ZEND_ASSIGN_CONCAT,
ZEND_JMP,
ZEND_FETCH_R,
ZEND_FETCH_DIM_R,
ZEND_IS_EQUAL,
ZEND_JMPZ,
ZEND_FETCH_DIM_R,
ZEND_ASSIGN_CONCAT,
ZEND_JMP,
ZEND_FETCH_R,
ZEND_FETCH_DIM_R,
ZEND_IS_EQUAL,
ZEND_JMPZ,
ZEND_FETCH_DIM_R,
ZEND_ASSIGN_CONCAT,
ZEND_JMP,
ZEND_ASSIGN,
ZEND_INCLUDE_OR_EVAL,
ZEND_RETURN
从op_array中可以看出明显的多个code block的结构。例如下图：
0x2：连续重复信号归并
php中，因为循环、数组操作等原因，vmc序列中常常会有很多连续出现的信号。
为了降低信号维度，可以采用抽样降维的方式，提取信号中的变化信息，过滤后的vmc序列如下：
ZEND_INIT_ARRAY,
ZEND_ADD_ARRAY_ELEMENT,
ZEND_FETCH_DIM_R,
ZEND_FETCH_DIM_R,
ZEND_JMP,
ZEND_FETCH_DIM_R,
ZEND_FETCH_DIM_R,
ZEND_JMP,
ZEND_FETCH_DIM_R,
ZEND_FETCH_DIM_R,
ZEND_JMP,
ZEND_FETCH_DIM_R,
ZEND_FETCH_DIM_R,
ZEND_JMP,
ZEND_FETCH_DIM_R,
ZEND_FETCH_DIM_R,
ZEND_JMP,
ZEND_INCLUDE_OR_EVAL,
ZEND_RETURN
