基于正负区间统计分析的信号趋势项自适应消除方法

An adaptive trend elimination method based on statistical analysis of positive and negative intervals of a signal

摘要

趋势项的存在会严重影响信号异常检测方法的准确性，因此需要在信号预处理过程中进行消除。高通滤波法是一种简单且有效的趋势项去除方法，应用最为广泛，但其参数选择基于经验，无法实现实时最优，因此方法的鲁棒性和泛化性差。以高通滤波法为基础，首先提出了信号正负区间统计分布曲线的概念，并基于这一概念引出了趋势项判断准则。在此基础上，设计了高通滤波器截止频率的实时优化选择方法以实现信号趋势的自适应消除，并对所提方法的影响要素进行分析。结合实际应用问题改进了算法，提高了迭代收敛速度，减少了计算量。最终利用石油管道泄漏监测问题中的压力信号数据证明了所提方法可以获得实时最优的滤波效果，实现趋势项的消除，且预设参数简单，算法计算量小。

Abstract

The signal trend term can seriously affect the accuracy of signal anomaly detection methods. Therefore， it is very important to eliminate the trend term in the process of signal preprocessing. High⁃pass filtering has been proved to be a simple and efficient trend term elimination method. In this work， the positive and negative interval statistical distribution of the signal is first obtained. Then， a judgment criterion for the trend item is proposed. On this basis， a real-time optimal selection method of the cut-off frequency of a high⁃pass filter is designed to realize the adaptive elimination of the signal trend， and the factors influencing the efficacy of the proposed method are analyzed. The algorithm was then improved to give the high convergence speed and reduced calculation amount required for practical application. Finally， the proposed method was employed using the pressure signal data in a pipeline leakage monitoring problem. The experimental results show that this method can realize real-time optimal selection of the cut-off frequency， and the calculation quantity of the algorithm is small.

Graphical abstract

关键词

趋势项消除 / 信号预处理 / 高通滤波 / 管道泄漏监测 / 截止频率

Key words

trend item elimination / signal preprocessing / high-pass filter / pipeline leakage monitor / cut-off frequency

引言

信号趋势项是存在于信号中的缓慢变化、周期长且频率低的成分^［1］。趋势项的产生有多种原因，如工作环境，设备故障或传感器、变送器零位漂移^［2］等。趋势项广泛存在于各类信号中，如电化学噪声数据^［3］、爆破振动信号^［4］、抽油机光杆位移加速度信号^［5］、运载火箭和导弹等的遥测信号^［6］、管道泄漏检测的压力或声波信号^［7］等。检测信号中含有趋势项将改变原始有效信号的时域和频域统计特性，进而对最终信号分析、异常检测结果产生不良影响。因此，趋势项消除是信号预处理的重要步骤。

Brockwell等^［8］用数学公式的形式描述了趋势项在信号中的存在方式

x(t)=m_t+s_t+y(t)（1）

式中，m_t为趋势分量，是一个缓慢变化的函数，s_t为一个表示周期变化分量的周期性函数，y(t)为信号有效成分。去趋势项的过程就是提取确定性分量m_t和s_t，得到有效分量—平稳时序数据y(t)的过程。

由于信号趋势项存在的广泛性及其对信号异常检测结果的不利影响，研究者们提出了多种去除信号趋势项的方法，如多项式拟合法^［9］、模型法^［10］、小波法^［11］、经验模态分解法^［12］和数字滤波法^［13］等。 上述研究已经被证明在趋势项的消除方面可以取得良好的结果，但仍然存在方法的预设参数无法根据滤波效果实时自适应调整的不足。以高通滤波法为例，滤波器通带频率对趋势项消除效果有重要影响，但只能预先设定该参数而无法实时自适应调整。不仅是在高通滤波法中，预设参数无法自适应调整的问题在小波法、经验模态分解法中同样存在。

针对趋势项消除方法存在的上述问题，本文设计了基于信号正负区间统计分析的趋势项判断准则和高通滤波器的通带频率实时优化选择方法，可以实现趋势项的实时自适应消除。首先提出信号正负区间统计分布这一概念，该分布曲线与信号直方图类似，但与信号直方图相比参数设定简单且更加平滑。再以正负区间统计分布为基础，提出趋势项判断准则。在此基础上，采用迭代法基于趋势判断准则实时确定最优通带频率，实现趋势项消除。所提方法预设参数少且选择简单，对信号先验知识要求少，具有较好的泛化性。在管道泄漏检测问题上的实验结果表明本文方法在消除趋势项的同时可以保留异常信号的有效成分，提高了分析和检测结果的准确性，是一种有效的信号预处理方法。

1 信号正负区间统计分布及基于分布曲线的趋势项判定准则

直方图是一种表示信号统计分布特性的图示方法。本文在直方图思想的基础上提出了信号正负区间分布曲线，该曲线横坐标为分组次序，纵坐标为以分组端点为零点对信号进行分割所获得的正区间和负区间的个数，其计算流程如下。

算法1

步骤1 采集信号数据p=\{p(\mathrm{1}),p(\mathrm{2}),\dots ,p(N)\}。

步骤2 对信号p进行去噪。

步骤3 计算该帧信号的最大值p_{\mathrm{m}\mathrm{a}\mathrm{x}}和最小值p_{\mathrm{m}\mathrm{i}\mathrm{n}}。

步骤4 确定分组数M。

步骤5 按照分组数确定步长

step=\frac{p_{\mathrm{m}\mathrm{a}\mathrm{x}}-p_{\mathrm{m}\mathrm{i}\mathrm{n}}}{M}（2）

获得等分线序列

\{lin{e}_k=p_{\mathrm{m}\mathrm{i}\mathrm{n}}+k\times step,k=\mathrm{0,1},\mathrm{2},\dots ,M\}

步骤6 以lin{e}_k的位置为坐标轴零点对p进行处理，得到双极性信号p_k=p-lin{e}_k，对p_k作正负区间划分，统计区间个数N{C}_k，其中k=\mathrm{0,1},\mathrm{2},\dots ,M。

步骤7 经过M次横坐标轴平移后，以序列号k为横坐标，以归一化后的N{C}_k为纵坐标，绘制信号正负区间分布曲线NCs=\{N{C}_k,k=\mathrm{0,1},\mathrm{2},\dots ,M\}。

以管道泄漏检测中的压力信号为例，对信号正负区间统计分布曲线的获取方法进行介绍。图1（a）、（b）展示了管道泄漏检测中实际测得的压力信号A和压力信号B，采样时长为10 min，采样频率为100 Hz，信号长度N=60 000，该段时间内压力信号的最大值记为p_{\mathrm{m}\mathrm{a}\mathrm{x}}，最小值记为p_{\mathrm{m}\mathrm{i}\mathrm{n}}。

通过观察波形可知，压力信号A中不包含趋势项，但是在第50 000到第60 000个采样点间（红色圈内）包含一个异常信号，压力信号B包含趋势项和异常信号（红色圈内）。

根据算法1计算压力信号A和压力信号B的step和lin{e}_k（当M取99时，lin{e}_k已用蓝色点划线在图1中标出，为避免线条展示过于密集，图中只展示lin{e}_k的k=0，1，9，18，27，…，90，99时的位置；步长step也已在图1中示意），进而得出N{C}_k，最终正负区间分布曲线（NCs）如图2所示。图3绘制了图1中信号A和信号B的直方图，组数为99。

对比图3和图1可见，直方图也可以反映信号中趋势项的存在与否。但是，从图2、图3可以看出，图3所示信号（尤其信号B）直方图曲线振荡更为剧烈，图2所示信号正负区间分布曲线较为平滑，因此相比直方图，将正负区间分布曲线作为趋势项判断准则，能提高迭代计算的收敛速度，在消除趋势项的同时，更好地保留信号的异常成分。

结合图1、图2和图3可以得到如下趋势项判断准则：在分组数M取合适值的条件下，当信号的正负区间分布曲线为单峰函数或近似单峰函数时，信号中不包含趋势项；当信号的正负区间分布曲线具有多个局部极值时，信号中包含趋势项。

2 最优通带频率实时选择方法

由于趋势项主要为低频信号，因此利用高通滤波器进行趋势项消除是最为简便且有效的手段。滤波器的截止频率即高通滤波器所能通过信号的最低频率，若信号的频率低于截止频率则不能通过。

2.1 基于趋势项判断准则的最优通带频率实时选择方法

在得到信号的正负区间分布曲线后，利用分布曲线的差分序列来判断正负区间统计分布曲线的单调性。基于趋势项判断准则来自适应地调节滤波器的截止频率直至趋势项被消除，具体步骤如算法2所示。算法2的总体流程如图4所示。

算法2

步骤1 令迭代次数n=\mathrm{0}。

步骤2 采集信号数据并进行去噪得到信号p，并赋值给p_n。

步骤3 计算信号p_n的正负区间分布曲线NCs， 并记录曲线的峰值点位置peakpo{s}_n。

步骤4 对NCs进行差分处理，得到差分序列D\_NCs： \{D\_N{C}_{\mathrm{0}},D\_N{C}_{\mathrm{1}},\dots ,D\_N{C}_{M-\mathrm{1}}\}。

步骤5 如果峰值点peakpo{s}_n之前的D\_NCs值都为负，峰值点peakpo{s}_n之后的D\_NCs值都为正，说明信号趋势项已被消除，停止迭代计算，得到去除趋势项后的信号p_n；否则，执行步骤6。

步骤6n=n+\mathrm{1}，对信号p作截止频率为n\times f_{\mathrm{c}}的高通滤波得到信号p_n；返回执行步骤3。

2.2 算法鲁棒性和泛化性分析

2.2.1 分组数M的优化选择及影响因素

由图4可知，所提算法有两个预设参数—分组数M和初始截止频率f_{\mathrm{c}}。最终截止频率需通过多次迭代得到，因此初始截止频率对算法结果的影响可忽略。M取值越大，算法计算量越高。对不同M值下的最终结果和计算时间进行分析，对比迭代终止时的下限截止频率。

图5为一对实测的管道上下游压力信号，可以看出两个信号均存在趋势项。选择巴特沃斯二阶高通滤波器，初始下限截止频率f_{\mathrm{c}}设为0.001 7 Hz。M的取值决定了正负区间分布曲线横坐标轴的范围，需保证可以呈现出正负区间分布曲线的形状和趋势，满足可以用差分序列进行趋势判断的条件，因此M不宜过小；此外，还应保证正负区间分布曲线的计算量不宜过大，以满足滤波器截止频率的实时优化选择的要求。显然，M越大正负区间分布曲线的计算量越大，同时结合实验信号的压力值波动范围，本文实验中的M取值不超过500。所以，在0~500范围内，取等分线lin{e}_k条数为整数50、100、200、300、400和500，对应分组数M=49，99，199，299，399，499（因｛line_k =p_min+k*step，k=0，1，2，…，M｝中k从0变化到M，共计M+\mathrm{1}个数，所以等分线lin{e}_k条数为M+\mathrm{1}）进行分析，比较不同的M取值条件下信号C和信号D的趋势消除效果和计算时间。其中，M=\mathrm{99}和M=\mathrm{499}时的信号趋势消除结果如图6、7所示，计算时间和迭代计算终止时的高通滤波器的下限截止频率如表1所示。

由以上结果可知，采用巴特沃斯二阶高通滤波器，在M取不同值的条件下信号的趋势最终都可被消除。从趋势消除效果的对比角度可得，不同的M值对趋势消除与否影响不大。从异常信号检测的角度，对比图6和图7中红圈标出的部分，可以看出该部分存在异常信号，当M取值过大（M=\mathrm{499}）时，异常信号信噪比变低，不利于信号异常检测。由表1可知，从计算时间的角度，随着M的增大，计算时间呈现逐渐增加的趋势，其中当M=49时，出现了一个例外的情况。因此，综合信号异常检测效果和计算时间，M值不宜过大或过小，本文取为99。

2.2.2 应用于实际工程对象的改进

在2.1节中算法2以差分序列在峰值点前后的正负性作为判据，迭代终止条件较为苛刻，在实际应用中受异常成分或其他干扰因素的影响会导致迭代次数过多、计算时间过长。此外，在管道泄漏检测应用中，严格按照该终止条件消除趋势得到的异常信号将消除部分负压波信号的暂态信息，不利于管道泄漏的准确检测。因此，为保证算法在应用于实际对象时的有效性，本文对算法进行了改进，设计了趋势消除比阈值r，以降低趋势项判据要求，简化计算。改进后的算法计算过程如下。

算法3

步骤1 确定阈值r，令迭代次数n=\mathrm{0}。

步骤2 采集信号数据并进行去噪得到信号p，并赋值给p_n。

步骤3 计算信号p_n的正负区间分布曲线NCs，并记录曲线的峰值点位置peakpo{s}_n。

步骤4 对NCs进行差分处理，得到差分序列\{D\_N{C}_{\mathrm{0}},D\_N{C}_{\mathrm{1}},\dots ,D\_N{C}_{M-\mathrm{1}}\}。

步骤5 在峰值点位置peakpo{s}_n之前（即）满足D\_N{C}_k\le \mathrm{0}的点的个数记为L\_cn{t}_n，在峰值点位置peakpo{s}_n之后（即peakpo{s}_n\le k\le M-\mathrm{1}）满足D\_N{C}_k\ge \mathrm{0}的点的个数记为R\_cn{t}_n，可得到rati{o}_n的计算公式

rati{o}_n=\mathrm{1}-\frac{L\_cn{t}_n+R\_cn{t}_n}{M}（3）

如果说明信号趋势项的消除已满足异常检测要求，停止迭代计算，得到去除趋势项后的信号p_n；否则，执行步骤6。

步骤6n=n+\mathrm{1}，对信号p作截止频率为n\times f_{\mathrm{c}}的高通滤波得到信号p_n；返回执行步骤3。

以实测管道上游压力信号C为对象，采用巴特沃斯二阶高通滤波器，初始下限截止频率f_{\mathrm{c}}=\mathrm{0.001}\mathrm{7}\mathrm{H}\mathrm{z}，观察随着迭代次数增加信号C的趋势消除效果和趋势消除比rati{o}_n之间的关系。图8展示了第1、第3、第13和第34次迭代计算时的信号C波形。分析可知，第1次高通滤波（下限截止频率为f_{\mathrm{c}}）后得到的信号波形还存在明显的趋势，负压波信号的暂态特性（异常成分）也没有丢失，如图8（a）所示。第3次高通滤波（下限截止频率为\mathrm{3}{f}_{\mathrm{c}}）后的信号波形已经较好地消除了趋势，但仍然保留了完整的负压波信号的暂态特性，如图8（b）所示。随着高通滤波器下限截止频率的增加，负压波信号的暂态特性逐渐丢失，具体表现为负压波信号的宽度逐渐变小，信号逐渐变尖，并且出现了负、正双极性异常信号（分别对应负压波信号的下降沿和上升沿），如图8（c）、（d）所示。信号C每次迭代的趋势消除比rati{o}_n和对应的高通滤波器下限截止频率如图9所示。表2为信号C第1、第3、第13和第34次迭代计算时的趋势消除比rati{o}_n，可以看出随着4次迭代信号C趋势消除效果的改变，其趋势消除比rati{o}_n越来越小。

由上述结果可见，当信号趋势消除比小于某一阈值时，不仅信号的趋势被消除，信号的异常（负压波信号的暂态特性）也得到消除，因此合理设置趋势消除比阈值r，当趋势消除比rati{o}_n小于阈值r时便停止迭代计算，可在保证趋势项消除效果的同时降低计算量，且较完整地保留信号异常成分，特别是信号的暂态特性，这对信号的异常检测非常有利，从而可实现截止频率的最优选择。该预设参数的选择对最终趋势项消除结果无影响，只改变计算时间，因此对阈值r的设置可利用几个含有异常成分的信号离线完成。结合表2和图8、图9可知，r可设在0.09~0.15之间，本文在后续实验中取r=0.1。

3 算法有效性验证

3.1 在管道泄漏检测上的工程应用

本文以管道泄漏检测问题为研究对象，对所提趋势项消除算法的有效性进行验证，并与其他方法进行对比。实验所采用的信号来自某化工企业石油管道，通过负压波法进行实时管道泄漏检测。采样频率为100 Hz，10 min为1帧，每帧信号共计60 000个采样数据点。利用巴特沃斯二阶高通滤波器，初始下限截止频率f_{\mathrm{c}}为0.001 7Hz，分组数M=\mathrm{99}，选择20个带有趋势项的异常实验信号，按照2.1节所述基本方法并结合2.2.2节面向工业对象的改进算法进行趋势项判断与消除，以验证方法的有效性。实验采用的电脑处理器为AMD FX-7500 Radeon R7，内存为4 GB。

在以往研究中，直方图可以用于判断信号中是否存在趋势项^{［3，14-15］}。表3为本文所提判断方法及其改进方法与以直方图为判断准则的方法的对比结果。从表3的20个实验信号结果来看，基于直方图的判断方法每次迭代均达到了迭代上限，且不能保证信号有效异常成分的保留。相比之下，本文改进后的方法在消除信号趋势的同时能更好地保留利于异常检测的信号的异常成分；且本文所提改进后方法平均计算时间为2.58 s，改进前方法平均计算时间为32.46 s，而基于直方图的方法的平均计算时间为65.24 s。可以看出改进后方法的计算时间明显缩短，可满足实时检测要求。

3.2 在流化床结块检测上的工程应用

为进一步验证方法的实用性，将改进后的方法应用于聚乙烯流化床结块故障检测过程中产生的带有趋势的声波压力信号。声波压力信号的气固流化床结块声波采集系统主要由聚乙烯流化床中试装置和信号采集传输系统两部分组成。流化床中试装置包括流化床反应器和四路传感器，信号采集传输系统包括放置于装置现场的监测终端和位于监控室的监控机。系统采用多路非介入式传感器，完全不会影响中试装置的生产运行，同时将监测终端放置于电气防爆箱内，使装置与危险源隔离，保证系统可以运行在安全可靠的环境中。流化床结块声波采集系统的具体设计可参考文献［16］。图10为传感器B和传感器D在2017年4月25日同一时间段内采集的带有趋势和异常的声波压力信号，信号采样频率为250 Hz，利用巴特沃斯二阶高通滤波器，初始下限截止频率f_{\mathrm{c}}为0.001 7 Hz，分组数M=\mathrm{99}，利用本文改进方法进行趋势消除后的信号波形如图11所示。

从图10和图11可以看出，本文所提方法在很好地消除了声波压力信号趋势的同时还保留了异常成分。经计算消除趋势的平均所需时间为2.83 s，满足实时检测要求。

4 结束语

本文提出了一种信号趋势项自适应消除方法，克服了传统高通滤波法最优截止频率难以实时选择的缺点。方法首先定义了信号正负区间统计分布曲线，并据此提出了趋势项判断准则；在此基础上，利用迭代法实现了最优截止频率的实时选择。应用基于高通滤波的趋势项消除算法时，所提方法的预设参数少、最优参数选择范围大、易于设置、鲁棒性强。将改进后方法应用在管道泄漏检测中的压力信号时，在消除信号趋势的同时还较完整地保留了信号的异常成分，计算时间相较于基于直方图的滤波方法和改进前方法有明显减少，满足实时检测的要求。趋势项广泛存在于各类信号中，本文方法可以很容易地推广到其他带有趋势项的信号预处理运算中，具有广泛的应用领域。

参考文献

[1] ZHANG P， CHANG J， QU B Y， et al. Denoising and trend terms elimination algorithm of accelerometer signals［J］. Mathematical Problems in Engineering， 2016， 2016： 2759092.

[2] HASHEMIAN H M， JIANG J. Pressure transmitter accuracy［J］. ISA Transactions， 2009， 48（4）： 383-388.

[3] LENTKA L， SMULKO J. Methods of trend removal in electrochemical noise data⁃overview［J］. Measurement， 2019， 131： 569-581.

[4] 张胜， 凌同华， 曹峰， 等. 模式自适应连续小波去除趋势项方法在爆破振动信号分析中的应用［J］. 爆炸与冲击，2017， 37（2）： 255-261.

[5] ZHANG S， LING T H， CAO F， et al. Application of removal trend method of pattern adapted continuous wavelet to blast vibration signal analysis［J］. Explosion and Shock Waves， 2017， 37（2）： 255-261. （in Chinese）

[6] 何鹏举， 冯亮. 加速度信号随机噪声及趋势项实时消除方法研究［J］. 电子设计工程， 2013， 21（14）： 18-22.

[7] HE P J， FENG L. Study on the real⁃time elimination method of random noise and trend terms in acceleration signal［J］. Electronic Design Engineering， 2013， 21（14）： 18-22. （in Chinese）

[8] DENG G Q， HAO G Y， LV Y J， et al. Research on intermittent failure mechanism of the electrical connector of the missile launch vehicle［M］//LONG S， DHILLON B. Man⁃Machine⁃Environment System Engineering.MMESE 2018. Lecture Notes in Electrical Engineering，vol 527. Singapore： Springer， 2019.

[9] YU Z R， SHAO S， LIU N， et al. Cable tension identification based on near field radiated acoustic pressure signal［J］. Measurement， 2021， 178： 109354.

[10] BROCKWELL P J， DAVIS R A. Introduction to time series and forecasting［M］. 3rd ed. New York： Springer， 2016.

[11] CAZALS F， POUGET M. Estimating differential quantities using polynomial fitting of osculating jets［J］. Computer Aided Geometric Design， 2005， 22（2）： 121-146.

[12] WANG H W， FENG P L. Fuzzy modeling of multirate sampled nonlinear systems based on multi⁃model method［J］. Journal of Systems Engineering and Electronics， 2020， 31（4）： 761-769.

[13] COTRONEI M， ROSSINI M， SAUER T， et al. Filters for anisotropic wavelet decompositions［J］. Journal of Computational and Applied Mathematics， 2019， 349： 316-330.

[14] LANG X， REHMAN N U， ZHANG Y F， et al. Median ensemble empirical mode decomposition［J］. Signal Processing， 2020， 176： 107686.

[15] PARTHIPAN C G， KOKIL P. Stability of digital filters with state⁃delay and external interference［J］. Circuits， Systems， and Signal Processing， 2021， 40（8）： 3866-3883.

[16] BREMER J， ALFONSI L， BENCZE P， et al. Long⁃term trends in the ionosphere and upper atmosphere parameters［J］. Annals of Geophysics， 2004， 47（2/3）： 1009-1029.

[17] BUNDE A， HAVLIN S， KANTELHARDT J W， et al. Correlated and uncorrelated regions in heart⁃rate fluctuations during sleep［J］. Physical Review Letters， 2000， 85（17）： 3736-3739.

[18] 王喆. 基于声信号的聚乙烯流化床结块故障监测方法研究［D］. 北京： 北京化工大学，2018.

[19] WANG Z. Agglomeration fault monitoring of polyethylene fluidized bed based on acoustic signal［D］. Beijing： Beijing University of Chemical Technology， 2018. （in Chinese）

基金资助

国家自然科学基金(61973023)