工业过程数据稳态判断

1 minute read

Published:

如化工过程等连续的工业过程往往需要处于稳态运行以满足安全性和产品质量等要求,但是由于过程上游输入的不确定性,以及内部设备和操作参数的变化,稳态并不唯一。过程在不同稳态之间的切换形成了暂态。如DCS等所记录的时序数据中往往同时包含了稳态和暂态数据,对其中稳态数据的识别有助于过程建模和优化等工作。


一、算法原理

1.1 连续和离散小波变换

不同的小波 $\Psi$ 具有不同形状,每类小波 $\it\Psi$ 的基波表达式为:

\[\it\Psi_{s,\tau}(t) = \frac{1}{\sqrt{s}}\it\Psi\left(\frac{t - \tau}{s}\right) \tag{1}\]

其中,参数 $s$ 为尺度参数(scale),其值越大,对应小波尺度越大,频率越低;$\tau$ 为时间位移参数(translation),用于控制小波在时间方向的移动。

连续小波变换(continuous wavelet transform, CWT) 中,$s$ 和 $\tau$ 连续取值,计算量较大且不同频率结果之间存在冗余(对于信号重构而言),而在 离散小波变换(discrete wavelet transform, DWT) 中,$s$ 和 $\tau$ 的取值则是离散的:

\[\it\Psi_{j, k}(t) = \frac{1}{\sqrt{s_0^j}}\it\Psi\left(\frac{t - k\tau_0 s_0^j}{s_0^j} \right) \tag{2}\]

其中,$j, k \in Z$,通常取 $\tau_0 = 1$,$s_0 = 2$。

注意,对于离散小波变换 DWT:

  • 尺度参数表征的是频率,在子小波中尺度参数以2的倍数增长(即小波的“长度”被“拉长”了2倍),那么子小波对应能检测到的频率值也会以1/2的倍数缩小。母小波所对应的频谱位于频率谱的高端,具有最大的频率谱范围,而其他的子小波的频率谱则依次向频谱图的低频端移动,同时它们所覆盖的频率谱范围也相应地递减;
  • 在理想情况下,所有的滤波器应首尾相接互相覆盖。比如在第一级分解中,尺度参数 $s_0$ 值较小,因此滤波细节信号 $D_1$ 保留了大量高频信息(高通滤波),剩余的近似信号 $A_1$ 则进入下一步采用更大尺度参数 $2\times s_0$ 的分解,依此类推。因此,小波分解所得近似信号的频率逐级递减。

1.2 单变量稳态数据识别

设待识别变量为 $X_{\rm r}$,其时序数据为 $X_{\rm r} = {x_{\rm r, 1}, \cdots, x_{\rm r, N}}$,其中 $x_{\rm r, t}$ 为第 $t$ 个时刻的数据。对 $X_{\rm r}$ 进行离散小波变换 DWT,得到 $X_{\rm r}$ 的 $J$ 级小波分解系数 $W_{\rm r} = {W_{\rm r, j, k}}$,其中 $j = 1, 2, \cdots, J$ 为分解级数,$k = 1, 2, \cdots, N$ 为分解系数的时间索引。保留 $W_{\rm r}$ 中的低频分量,将其余分量置零,然后对 $W_{\rm r}$ 进行逆小波变换,得到新的时序数据 $X$。

接下来,从滤波后的 $X = {x_1, x_2, \cdots, x_N}$ 中提取稳态数据。设在时刻 $t$ 的记录 $x_t = f(t)$,则此时 $X$ 变化的一阶和二阶导数分别为 $f’(t)$ 和 $f’‘(t)$,如果两个导数变化绝对值均低于设定阈值,则可认为此时 $X$ 处于稳态。计算过程如下:

  1. 确定稳态判断阈值参数 $T_s$、$T_w$ 和 $T_u$。其中,$T_s$ 和 $T_w$ 分别为过程稳态数据的一阶和二阶导数绝对值的90%或95%分位数。$T_u = \alpha \cdot T_s$,其中 $\alpha$ 为超参数,为 (2, 5) 之间可调的整数参数。

  2. 对于每个时刻 $t$:
    • 首先,计算系数 $\gamma(t)$:

      \[\gamma(t) = \left\{ \begin{align*} &0, \text{if } |f''(t)|\leq T_w \\ &\frac{f(t)-T_w}{2 T_w}, \text{if } T_w < |f''(t)|\leq 3T_w\\ &1, \text{else} \end{align*} \right.\]
    • 接下来,计算系数 $\theta(t)$:

      \[\theta(t) = (1 + \gamma(t)) |f''(t)|\]
    • 然后,计算变量的瞬时稳态系数 $\beta(t)$:

      \[\beta(t) = \left\{ \begin{align*} &0, \text{if } \theta(t) > T_u \\ &0.5\left(\cos\left(\frac{\theta(t) - T_s}{T_u - T_s}\right) + 1\right), \text{if } T_s < \theta(t) < T_u \\ &1, \text{else} \end{align*} \right.\]

      稳态系数 $\beta(t)$ 的取值范围为 [0, 1]。当 $\beta(t) = 0$ 时,表示时刻 $t$ 处于非稳态,$\beta(t)$ 的取值越接近1,表示过程在时刻 $t$ 越平稳。

    • 最后,使用稳态系数阈值 $T_{\beta} = 1 - \varepsilon$ 对当前时刻的稳态系数 $\beta(t)$ 进行判断,其中 $\varepsilon \leq 0.1$ 为超参数。若 $\beta(t) \geq T_{\beta}$,则认为时刻 $t$ 处于稳态,否则处于非稳态。

  3. 当对所有时刻 $t$ 都进行了稳态判断后,便获得了稳态和非稳态数据的时刻记录 $T_{ss}$ 和 $T_{ns}$。如果要求稳态数据长度不小于某一阈值 $L$,则进一步地从 $T_{ss}$ 中提取出所有满足单段连续长度不小于 $L$ 的稳态数据段,最后形成稳态时刻片段集合 $T_{ss}^{\rm seg}$,并提取对应的稳态数据片段集合 $\left{X_{ss}^{\rm seg}\right}$。

1.3 多变量稳态数据识别

在含有 $P$ 个变量的过程中,若每个变量均达到稳态,则过程整体达到稳态。因此,可通过如下的加权方式获得过程整体的瞬时稳态系数:

\[B(t) = \prod_{i=1}^{P} \beta_i(t)^{w_i / \sum_{i=1}^P w_i} \tag{3}\]

其中,$w_i$ 为每个变量的权重,可由用户指定。若 $B(t) \geq T_{\beta}$,则认为整个过程在时刻 $t$ 处于稳态。同样地,可按照上述2.2节中第3步求解满足最低稳态样本量要求 $L$ 的稳态样本片段集合。


二、算例实现

2.1 单变量稳态数据识别

从样本中,变量 $x_1$ 的原始数据记录如下:

其中红色标记区域为人工挑选的用于计算参数 $T_s$、$T_w$ 的稳态数据片段,长度为500。对该片段采用db4小波进行变换,得到滤波结果如下:

可见,滤波后结果更为平滑。最终求得参数值:$T_s = 0.28$、$T_w = 0.21$、$T_u = 0.56$。

接下来,对整个 $x_1$ 信号进行稳态识别,得到结果如下:

2.2 多变量稳态数据识别

在2.1节的基础上,对同时含有多个变量过程的稳态进行识别。过程变量 $X_1$、$X_2$ 和 $X_3$ 的原始数据记录如下:

其中,红色标记区域为人工挑选的用于计算参数 $T_s$、$T_w$ 的稳态数据片段,长度为500。计算获得各变量稳定系数和过程稳定系数变化如下:

最终获得过程稳态识别结果如下:


三、参考文献

  1. T. Jiang, B. Chen, X. He, et al. Application of Steady-State Detection Method Based on Wavelet Transform, Computers & Chemical Engineering, 2002.