Python在振动信号处理中的应用（一）：振动信号消除多项式趋势项

一、概述

  在振动测试过程中，传感器或采集设备由于自身性能问题或受环境干扰（如温度、供电等影响），容易产生零漂，而偏离基线，偏离基线的大小甚至还会随时间变化。偏离基线随时间变化的过程称之为信号的趋势项。信号的趋势线将直接影响信号分析的准确性，应该将其消除。常用的去除趋势项的方法是最小二乘法，以下介绍该方法的原理。

二、最小二乘法去除趋势项原理

  令实测振动信号的采样数据为$x_k(k=1,2,3,\cdots n)$ ，由于采样数据是等时间间隔的，为简化起见，令采样间隔$\Delta t=1$ ，设置一个多项式$(2.1)$：
$${\hat{x}}_k=a_0+a_{1}k+a_2{k}^2+\cdots +a_m{k}^m(k=1,2,3,\cdots n)\text{\hspace{1em}(2.1)}$$
  其中：${\hat{x}}_k$ 为拟合数据； $a_i$为待定系数； $m$为阶数。
  对所有采样数据均进行计算，可得：
$$\hat{\mathbf{X}}=\mathbf{K}\mathbf{A}\text{\hspace{1em}(2.2)}$$  其中：$\hat{\mathbf{X}}=\left[\begin{array}{c}{\hat{x}}_1\\ {\hat{x}}_2\\ ⋮\\ {\hat{x}}_n\end{array}\right],\mathbf{K}=\left[\begin{array}{cccc}1& 1^2& \cdots & 1^m\\ 1& 2^2& \cdots & 2^m\\ ⋮& ⋮& ⋮& ⋮\\ 1& n^2& \cdots & n^m\end{array}\right],\mathbf{A}=\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_m\end{array}\right]$
  求取函数${\hat{x}}_k$的待定系数$a_i(i=0,1,\cdots m)$，使得函数${\hat{x}}_k$与采样数据$x_k$之间的误差平方和最小，即：
$$\min (E)=\min \left[(\mathbf{X}-\hat{\mathbf{X}}{)}^T(\mathbf{X}-\hat{\mathbf{X}})\right]\text{\hspace{1em}(2.3)}$$  结合式$(2.2)$，求$E$对$a_i(i=0,1,\cdots m)$的偏导，则有：
$$\frac{\partial E}{\partial A}=2\left({\mathbf{K}}^T\mathbf{K}\mathbf{A}-{\mathbf{K}}^T\mathbf{X}\right)\text{\hspace{1em}(2.4)}$$  令偏导数为零，即为误差平方和的极小值，因此有：
$$\mathbf{A}={\left({\mathbf{K}}^T\mathbf{K}\right)}^{-1}{\mathbf{K}}^T\mathbf{X}\text{\hspace{1em}(2.5)}$$  通过以上公式即可求个m+1个待定系数。当m=0时，求取的趋势项为采样数据的算术平均数；当m=1时，求取的趋势项为线性趋势项；当m≥2时，求取的趋势项为曲线趋势项。消除趋势项的公式如下：
$$y_k=x_k-{\hat{x}}_k=x_k-\left(a_0+a_{1}k+a_2{k}^2+\cdots +a_m{k}^m\right)\text{\hspace{1em}(2.6)}$$  通常去m=1~3来对振动采样数据进行去多项式趋势项分析。

三、python模块介绍

  最小二乘法求回归方程系数，是机器学习最基本的算法。本文采用python中scikit-learn模块完成振动信号去除趋势项算法的编制。在开始算法描述之前，我们先介绍一下scikit-learn及其中的线性回归LinearRegression、多项式PolynomialFeatures以及工作流Pipeline。

3.1 scikitlearn简介

  scikitlearn是一个Python第三方提供的非常强力的机器学习库，它包含了从数据预处理到训练模型的多个方面，包括：预处理、分类、回归、聚类、降维和模型选择。为了方便使用，scikitlearn函数可以归为转化器(Transformer) 和估计器(Estimator)两类。

（1）转化器(Transformer)

  转换器用于对数据的处理，例如标准化、升维、降维以及特征选择等等。基本上估计器都会有以下几个方法：
  1. fit(x, y)：该方法接受输入和标签，计算出数据变换的方式。
  2. transform(x)：根据已经计算出的变换方式，返回对输入数据x变换后的结果（不改变x）。

（2）估计器(Estimator)

  估算器就是模型，它用于对数据的预测或回归。基本上估计器都会有以下几个方法：
  1. fit(x, y)：传入数据以及标签即可训练模型，训练的时间和参数设置，数据集大小以及数据本身的特点有关
  2. score(x, y)用于对模型的正确率进行评分(范围0-1)。但由于对在不同的问题下，评判模型优劣的的标准不限于简单的正确率，可能还包括召回率、查准率等其他的指标，特别是对于类别失衡的样本，准确率并不能很好的评估模型的优劣，因此在对模型进行评估时，不要轻易的被score的得分蒙蔽。
  3. predict(x)用于对数据的预测，它接受输入，并输出预测标签，输出的格式为numpy数组。我们通常使用这个方法返回测试的结果，再将这个结果用于评估模型。
  以上仅仅是简单的概括scikitlearn的函数的一些特点。scikitlearn绝大部分的函数的基本用法大概如此。但是不同的转换器和估计器会有不同的属性。对于机器学习来说模型的好坏不仅取决于你选择的是哪种模型，很大程度上与模型的参数设置有关。因此使用scikitlearn的时候一定要去看看scikit-learn官方学习资料，以便对参数设置进行调整。

3.2 线性回归LinearRegression

  线性回归(Linear Regression)是利用数理统计中回归分析，来确定两种或两种以上变量间相互依赖的定量关系的一种统计分析方法。线性回归利用线性回归方程的误差最小二乘法，对一个或多个自变量和因变量之间关系进行建模，其计算原理见第二章节。只有一个自变量的情况称为简单回归，大于一个自变量情况的叫做多元回归。
  本文介绍的是scikit-learn中的线性回归模型，LinearRegression官方学习资料，其基本用法如下：
  sklearn.linear_model.LinearRegression(fit_intercept=True, normalize=False, copy_X=True, n_jobs=None, positive=False)
  输入参数：
  fit_intercept：bool型，默认True，是否计算该模型的截距。如果使用中心化的数据，可以考虑设置为False，不考虑截距。
  Normalize：bool型，默认False，是否对数据进行标准化处理。建议将标准化的工作放在训练模型之前，通过设置sklearn.preprocessing.StandardScaler来实现，而在此处设置为false当  fit_intercept设置为false的时候，这个参数会被自动忽略。如果为True，回归器会标准化输入参数：减去平均值，并且除以相应的二范数
  copy_X：bool型，默认True，是否对X复制。如为false，则即经过中心化，标准化后，把新数据覆盖到原数据上
  n_jobs：int or None, optional，默认None,如果选择-1则代表使用所有的CPU。计算时设置的任务个数，这一参数的对于目标个数>1（n_targets>1）且足够大规模的问题有加速作用

  输出参数：
  coef_：对于线性回归问题计算得到的feature的系数。如果输入的是多目标问题，则返回一个二维数组(n_targets, n_features)；如果是单目标问题，返回一个一维数组(n_features,)。
  rank_：矩阵X的秩，仅在X为密集矩阵时有效
  singular_：矩阵X的奇异值，仅在X为密集矩阵时有效
  intercept_：截距，线性模型中的独立项
  使用的方法：
  fit(X, y[, sample_weight]) ：训练模型，sample_weight为每个样本权重值，默认None
  predict(X)：模型预测，返回预测值
  score(X, y[, sample_weight])：模型评估，最优值为1，说明所有数据都预测正确

3.3 多项式PolynomialFeatures

  PolynomialFeatures主要用于生成多项式，PolynomialFeatures官方学习资料，其基本用法如下：
  sklearn.preprocessing.PolynomialFeatures(degree=2, interaction_only=False, include_bias=True, order=‘C’)
  degree：int型，默认值= 2，多项式阶次。
  interaction_only：bool型，默认为False，是否值包含交互项。如果为True，那么就不会有特征自己和自己结合的项（没有a ^ 2、b ^ 2）
  include_bias：bool型，默认为True，是否有全是1的一项。如果为True（默认值），那么结果中就会有 0 次幂项，即全为 1 这一列
  order：str型，{‘C’, ‘F’}，‘C’，密集情况下输出数组的顺序。F阶的计算速度更快，但可能会降低后续的估计速度。

3.4 工作流Pipeline

  Pipeline可以将许多算法模型串联起来，比如将特征提取、归一化、分类组织在一起形成一个典型的机器学习问题工作流。主要带来两点好处：
  1. 直接调用fit和predict方法来对pipeline中的所有算法模型进行训练和预测。
  2. 可以结合grid search对参数进行选择。
Scikitlearn中Pipeline为高级用法，在后面的文章中进行讲解。

四、python代码实现

  为了能更好的展示算法去趋势项的效果，本文构造了带趋势项的模拟信号，随后通过去趋势算法处理该信号，具体如下。

4.1 构造模拟信号

  构造模拟信号$signal=\sin (2\pi f)+0.1+0.1t+0.1t^2$。

# 导入模块
import numpy as np
import matplotlib.pyplot as plt

# 构造模拟信号
fre = 10
t = np.arange(0, 1, 0.001)
sig = np.sin(2*np.pi*fre*t) + 0.1 + 0.1 * t + 0.2 * t**2
sig = sig.reshape(-1, 1)
# 展示模拟信号
plt.plot(t, sig)
plt.xlabel('time(s)')
plt.ylabel('signal(m/s^2)')
plt.show()

图4.1 模拟信号

4.2 编写去除趋势项主代码

# 导入模块
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression

# 构造去除趋势项函数
def trend_remove(t, sig, degree):

    # 构造多项式回归模型
    poly_reg = Pipeline([
        ('poly',PolynomialFeatures(degree = degree)),
        ('lin_reg',LinearRegression(fit_intercept=False))])
    # 训练多项式回归模型（即求待定系数）
    poly_reg.fit(t, sig)
    # 去除趋势项
    sig_predict = poly_reg.predict(t)
    sig_res = sig - sig_predict
    # 输出趋势项待定系数
    model = poly_reg.named_steps['lin_reg']
    print('趋势项待定系数为：')
    print(model.coef_[0])
    # 返回去趋势项结果
    return sig_res

4.3 展示算法效果

sig_res = trend_remove(t, sig, 2)
# 展示去除趋势项后的信号
plt.plot(t, sig)
plt.plot(t, sig_res)
plt.xlabel('time(s)')
plt.ylabel('signal(m/s^2)')
plt.legend(['sig','sig_res'])
plt.show()

图4.2 去除趋势项效果

  从图中可以看出，虽然没有完全去除趋势项，但是效果已经很不错了，大功告成。

五、TIPS

  文中用到了sklearn模块中部分函数，建议大家先照做，后理解，再入研究，用啥学啥才是好工科生（此处省略大哭表情）。当然学有余力的同学，可以查看官方文档：scikit-learn官方学习资料，学习更多相关内容。