原子稀疏分解范文

原子稀疏分解范文（精选7篇）

原子稀疏分解 第1篇

风能作为可再生能源中成本较低、技术较成熟、 可靠性较高的新能源,近年来发展很快并开始在能源供应中发挥重要作用。 随着风电场规模的增大,风速的波动性[1-4]和非平稳性成为制约风电大规模、高效并网的严峻问题。 风电功率预测技术是解决风电波动[5]、风电并网[6]和电网调度[7]的关键技术之一, 这也对风电功率的预测精确提出了更高的要求。

为得到较高的预测精度,国内外很多研究集中在构造合适的预测模型。 根据输入量的不同,现有预测模型可以分为物理模型、统计模型和物理统计混合模型:物理模型[8-9]使用如气象学(数值天气预报等)、地质学(山岳形态等)和风电机组的技术特征(轮毂高度、功率曲线和推力系数)等信息作为模型输入量,目的是得到本地风速的最佳估计值,进而利用模式输出统计(MOS)方法减小预测残差;统计模型[10-12]使用解释变量和在线测量方法,通常使用如递归最小二乘法和人工神经网络ANN(Artificial Neural Network)法等递归技术;混合模型[13-15]作为最优模型,先得到风电机组区域内的气流等物理量,再使用先进的统计模型补充物理模型得到的信息,因而能够得到更为精确的预测值。 由于物理方法中风电场周围的物理信息等对预测结果的精确度有很大影响,而统计方法能够根据风电场自身的特点和位置,随时修改预测模型参数,可以得到比较高的准确度。

国内对风电功率测报与预测提出了最新的技术要求,2009年国家电网公司发布企业标准Q/GDW 392—2009《风电场接入电网技术规定实施细则(试行)》中明确规定,风电功率预测系统应能通过专网向调度机构上报相关数据,应至少具备日前预报功能和超短期预报功能,每日12:00之前向调度机构申报次日日前风电功率预测曲线,根据超短期预测结果,滚动调整2 h以后的风电功率预测曲线[16]。 国外主要集中在预测模型的研究,比国内起步早,技术也相对成熟。 为得到提前0.5~36 h的预测值,丹麦科技大学(DTU)提出一种计及遗忘因子的自适应递归最小二乘估计法的风电功率预测模型;马德里卡洛斯第三大学提出Sipreolico模型,该模型由9个自适应非参数统计模型组成,使用递归最小二乘算法或者卡尔曼滤波算法循环计算;TrueWind公司提出一种EWIND模型,该模型使用一次性参数设计方法研究顺风向NWP模型输出量的局部效应[17]。

现有的风电功率预测建模方法中,很少有考虑原始风电功率序列的非平稳特性的方法,ANN是一种应用广泛的风电功率预测建模方法,但由于其自适应训练的收敛性受步长、隐含层神经元个数、隐含层输出函数和输出层输出函数等因素影响,训练时间较长,往往不能完全映射风电功率的非平稳特性。 因此本文采用一种具有很强的非平稳信号跟踪、预测能力的信号处理新方法———原子稀疏分解ASD (Atomic Sparse Decomposition) 法, 作为ANN的前置分解手段。 现实中的风电功率具有很强的非平稳性,可看作具有多个不同参数的原子分量和残差分量的叠加,其非平稳性导致原子参数不断变化。 对原子分量进行自预测,残差分量进行ANN预测,叠加后得到最终预测结果。 较之常规的ANN预测方法,所提预测方法具有更好的处理非平稳特性的能力。 将该方法应用于国内外某几个实际风电场的风电功率预测中,取得了令人满意的结果,预测精度也符合要求。

1 ASD理论

1.1基本概念

近年来ASD技术在信号处理领域成为热点。 该方法源于文献[18]提出的信号在过完备原子库上分解的思想,在信号建模、压缩、特征提取等方面具有重要作用。 ASD采取的是一种贪婪的自适应分解策略,其原子库是高度冗余(过完备)的,以保证任意信号都可以从中自适应地选择一组最佳的原子来表示,使得分解结果非常稀疏,此过程则为稀疏分解。 核心问题是如何自适应地寻找最佳匹配原子及其系数。

当前,ASD法在电力系统中的应用才刚刚起步,且主要集中在电能质量扰动信号特征提取[19-20]和低频振荡模式识别[21]的研究。

1.2字典集的构造

在ASD中,原子通常由一般性的核函数表示, 在信号处理领域,多种核函数可用来表示原子,例如正弦函数[22]、Chirp函数[23]。 本文采用的核函数为高斯函数,如下式所示:

其中,g(o)为高斯核函数;c和 σ 分别为中心和尺度参数。 选择不同的中心和尺度参数,可以构造一系列不同的原子,这些原子的集合称为字典集。 图1中列举了3个不同的原子,其中,g1表示中心为0、 尺度为2的原子;g2表示中心为2、尺度为2的原子; g3表示中心为0、尺度为3的原子。

1.3双字典集的ASD

在迭代分解过程中,每次迭代的待选择原子可以分为2类:先前已经被选择过的旧原子和尚未被选择过的新原子。 因此,过完备字典集可以划分为2个分离的字典集:一个由旧原子构成,一个由新原子构成。

在开始阶段,所有的原子均属于新字典集,在分解过程的前面几次迭代中,大部分被选择的最优原子属于新字典集,随着迭代的继续,旧字典集慢慢增大。 当旧字典集足够大时,被选择的最优原子大部分属于旧字典集。 从稀疏性的角度看,新原子的选取不利于分解的稀疏性,为达到稀疏分解的目的应尽可能从旧字典集中选取原子。 因此,文献[24-25]提出了一种有利于分解稀疏性的最优原子选取流程, 具体第k步迭代描述如下。

根据第k步以前的分解结果,过完备字典集已被分为旧和新2个字典集。 由于2个信号的内积(各时刻2个信号乘积之和)描述了它们的线性相关性:内积的绝对值越大,2个信号的相关性越强;内积为0,2个信号线性无关。 所以,分别计算残差信号和2个字典集中各个原子的内积,并选出每个字典集中最大的内积cold和cnew,其对应的旧和新字典集中的原子分别用 Φold和 Φnew表示。

如果cold≥cnew,选Φold作为该次迭代中的最优原子,即Φopt=Φold,cold作为该最优原子的迭代系数,即copt=cold。显然,第k步之前的分解使得旧字典集中的每个原子都已有一个分解系数。因此,需将copt添加到被选原子Φold的分解系数上进行累加。最后,更新第k步的残差信号,即R(k)=R(k-1)-coptΦopt。

如果|cold|<|cnew|,么|k步|的最优原子按照如下步骤进行选择。

a. 分别计算旧和新字典集上的残差:

b. 计算相对误差re:

其中,‖·‖表示信号的欧氏范数。

c.通过给定的阈值T决定最优原子:如果re≤T,选Φold作为该次迭代中的最优原子,后续的计算过程与cold≥cnew情况下的相同;如果re>T,选Φnew作为该次迭代中的最优原子,更新变量,即令copt=cnew,Φopt=Φnew,R(k)=Rnew,把该原子添加到旧字典集中,并从新字典集中删除,系数copt作为该原子的分解系数。

d. 更新阈值。 通过给定阈值T,在旧和新字典集中选择最优原子。 为保证收敛和稳定性,T是一个随着迭代步数而递减的函数,本文采用的是模拟退火算法中的退火函数:

其中,0.7≤α<1;T0为初始温度,并且设定小于1;k为当前迭代步数;N为退火速度因子。 随着迭代步数的增加,T趋于0。 相关算法流程如图2所示。

2 ASD-ANN滑动预测模型

2.1数据预处理

在进行模型预测时,当输入或输出向量的各个分量量纲不同或大小相差很大时,应对不同的分量在其取值范围内分别进行归一化处理。 考虑到本文所采用各分量的物理意义相同且为同一量纲,故采用先在整个数据范围内确定最大值和最小值再进行统一的归一化变换处理,进而将模型输入输出变换为[0,1]区间的值,具体归一化公式如下:

其中,xput为模型的输入或输出分量;x軇put为经过归一化处理后的输入或输出分量;xmax和xmin分别为模型输入或输出量的最大值和最小值。

2.2 ASD预测模型建立

风机的发电量由于受到平均风速、主导风向、最大风速、极大风速、气温、气压、空气密度、雷暴、电线覆冰等因素的影响,使得风电功率本质上具有很强的非平稳性。 ANN能较好地拟合输入输出数据间的非线性关系,但离线训练的ANN无法自适应输入信号的非平稳性,不能完全映射其特性。 因此,本文采用ASD算法对风电功率进行滑动分解,并用残差信号取代原始信号作为ANN的输入量预测下个时刻的残差信号。 由于残差信号的能量(相对于原始信号)很小,这样能极大地避免具有主导能量的信号成分(原子的线性组合)的非平稳性对ANN预测产生影响。 模型结构见图3,具体的预测过程描述如下。

a. 对风电功率数据进行n次ASD,见式(6)。其中,r(t)为残差分量;aj(t)为第j个原子分量,等于原子与其分解系数的乘积。

b. 将残差信号作为ANN残差预测模型的输入量进行下一时刻的残差预测。

c. 根据ASD的表达式自预测下一时刻的原子分量值。

d. 将ANN残差预测结果与ASD预测结果叠加即得到下一时刻的原子稀疏分解-人工神经网络(ASD-ANN)模型预测值,如图3所示。 图中,x(t)表示风电功率的时间序列,x(t+1)表示模型预测值。

2.3 ANN残差预测模型

文献[26-28]中已说明ANN可以用来拟合能量较小的风电功率残差分量。 因此,本文采用常用的ANN作为残差预测模型,输入节点数设定为15,根据Kolmogorov定理,中间节点数设定为31,输出节点数设定为1,对应于前15个时刻的残差分量,中间层采用正切Sigmoid传递函数,输出层采用对数Sigmoid传递函数,输出节点对应于下一时刻的风电功率残差分量预测值,模型网络结构如图4所示。

2.4滑动预测模型

一般而言,一个确定的风电功率序列经过ASD- ANN预测模型前置分解后,必然得到一组参数确定、 平稳且具有主导能量的原子分量和非平稳、随机性强但能量小的残差分量。 由于原子分量占有主导作用,这种方法从本质上讲还是片面地把风电功率数据当作平稳序列来处理。 因此,本文提出一种滑动预测方法,通过2.2节建立50个最佳的预测模型。 采用最新的输入变量和对应的不同模型来对下一个15 min风电功率进行滑动预测。 每个模型对应的输入变量不同:模型Mi以第i点前的400点风电功率作为输入变量,对第i点预测;Mi+1利用第i点的实测值和第i点前的399点风电功率作为输入变量,依此类推。

采用时间尺度确定的时间窗进行下一时刻的预测,得到预测值后时间窗滑动向前推进一个时刻,继续类似的预测。 该方法的优点在于,虽然时间窗内的分解得到的是平稳的原子分量,但随着时间窗不断地向前滑动,时间窗与时间窗之间的原子分量参数是不断变化的,当训练样本足够大时,时间窗个数很多,可视为对风电功率数据进行了非平稳化处理; 同时,随着时间窗的滑动,原子分量自适应地调节自身的参数,以适应非平稳的风电功率数据,极大地增强了预测模型的泛化能力。

2.5预测结果校正

目前国内外的风电功率预测模型中,风电功率预测精度较差,使用滑动模型进行风电功率预测时有必要对预测结果进行校正。 采用线性回归方法进行校正[29],校正模型如下:

其中,PASD，t为采用ASD-ANN模型的t时刻的风电功率预测值;PcASD，t为校正后的t时刻预测值;eASD，t= a + bPASD，t为t时刻的预测误差,a和b为参数,可采用最小二乘法计算得到,由历史风电功率及其误差样本数据进行估计,方法如式(8)、(9)所示。

其中,Nc为历史样本数量;eASD,i = PASD,i - Pmeas,i为历史风电功率预测误差,Pmeas，i为风电场实测风电功率数据,PASD，i为实测数据的模型训练值。

综上所述,可以得到本文提出的ASD-ANN风电功率滑动预测模型结构如图5所示。

3算例分析

3.1算例1

为验证本文采用的基于双字典集的ASD算法的有效性,使用MATLAB 7.10进行算法编程,构造一个典型非平稳的测试信号,如下所示:

其中,i∈ [-15,15],从中均匀抽取600个采样点,e为随机白噪声信号。 根据2.3节所述,取前400个点进行ASD,残差信号序列作为ANN的训练样本, 后200个点作为测试样本。

根据文献[24-25],计算过程中的迭代分解参数 α = 0.935,初始温度T0= 0.09,退火速度因子N = 2.5, 迭代终止阈值为1.0×10-6。 图6为构造的测试信号, 图7为预测结果。

当i∈[-10,10]时,分析比较单字典集和双字典集的ASD算法仿真效果如表1所示。 定义测试信号的绝对平均误差eNMAE和均方根误差eNRMSE的公式分别如式(11)、(12)所示。

其中,y(i)为测试信号;y赞(i)为预测结果;M为预测样本个数。

由图6、图7和表1可见:

a. 相比基于单字典集的ASD算法,基于双字典集的ASD算法预测结果的绝对平均误差和均方根误差均有一定程度的降低;

b. 在无噪声和加噪声的情况下,基于双字典集的ASD算法的预测效果均更加精确;

c. 完成仿真计算的时间仅为3.485 s,约减少1 s的时间。

3.2算例2

选取2006年5月10日至2006年5月24日内国内某风电场58台风电机组中的某一台风电机组输出功率数据,时间分辨率为15 min,总共1440个数据。 选取其中前450点作为样本数据进行实验, 如图8所示。

选取图8中前400个数据点分解的残差信号作为ANN的训练样本,后50个数据点作为测试样本。 预测结果如图9所示,参数如表2所示。

作为对比,采用常规ANN预测风电功率[27]。 由于预测误差大小和风机容量有直接关系,为了定量评价预测效果,采用国际上普遍的归一化绝对平均误差eNMAE和归一化均方根误差eNRMSE为依据,定义如下:

其中,x(i)为实际值;赞x(i) 为预测值;Pcap为风机的额定容量。 表3分别列出了常规ANN法、校正前后ASD-ANN风电功率的绝对平均误差和均方根误差。

由图9和表3可见:

a. 校正后得到的风电功率预测结果相对于校正前得到较大的改善,而常规ANN法不能较好地自适应风电功率的强非平稳性,导致其预测误差最大;

b. ASD-ANN滑动预测的误差总体上小于ANN预测的误差;

c. 对预测结果的校正只能有限地提高预测精度,而采用合适的预测方法、充分考虑风电功率的非平稳性才是提高预测精度的主要方法。

3.3算例3

为验证所提预测方法的鲁棒性,采用统计学方法分析单台风机在不同时段和不同风机在同一时段的预测效果。

3.3.1单台风机在不同时段预测误差统计分析

选取2001年1月1日至2008年6月23日内国外某风机输出功率数据作为样本。 统计过程中,每600个采样点作为一个时间段记录一次预测结果, 总共统计50个时间段,即30000个采样点。 计算每个时间段的绝对平均误差并统计出现的频率,结果如图10所示。 图中,频率直方图指样本观测值在各区间单位长度内的频数与总数的比值。

采用统计学常用的正态分布拟合方法分别拟合图10中2种预测方法的绝对平均误差(置信区间取值均为95%),得到均值的估计值分别为16.13%和20.72 %,标准差的估计值分别为1.26 % 和1.46 %。 可见,就单台风机在不同时段的预测误差统计分析而言,统计结果有效验证了本文预测方法的鲁棒性, 绝对平均误差由常规ANN预测方法得到的20.72% 降至本文所提方法的16.13%,精度提高约30%。

3.3.2不同风机在同一时段的预测结果对比分析

选取某风电场4台风电机组输出功率数据,时间为2006年5月10日至2006年5月16日,分辨率为15 min,每台机组选取600个采样点,前400个点作为训练样本,后200个点作为测试样本[30]。 不同风机在同一时段的预测结果如表4所示。

由图10和表4可得出以下结论。

a. 根据文献[31]提出的相关系数判断随机变量相关性的准则:ANN预测方法的预测结果与实测结果相关系数约为0.36~0.49,二者具有低度相关;本文所提预测方法的预测结果与实测结果相关系数约为0.50~0.57,二者具有中度相关。

b. 预测结果与实测结果的偏差分析可知,ANN预测方法的绝对平均误差约为17%~19%,均方根误差约为22%~24%,预测数据合格率(eNRMSE< 10 %) 约为57%~65%;本文所提预测方法的绝对平均误差约为15%~18%,均方根误差约为20%~23%,预测数据合格率(eNRMSE<10%)约为68 % ~ 74 %。

可以看出,就不同风机在同一时段的预测结果对比分析而言,所提预测方法提高了预测结果与实测结果的相关度和预测数据合格率的统计区间,降低了绝对平均误差、均方根误差计算值的统计区间。

4结论

ASD是近年来信号稀疏分解领域的新热点,可实现对信号更加灵活、简洁和自适应的表示。 针对风电功率的非平稳性,本文将双字典集ASD和ANN相结合的方法引入到风力发电功率时间序列的预测中。 通过滑动的ASD,信号中能量占主导地位的非平稳部分能较准确地被预测。 同时,ANN可以较好地映射出较平稳的残差信号输入输出之间的非线性关系。 将该方法应用于某几个实际风电场的功率预测中,结果表明,该方法可以有效地处理风电功率的非平稳性,能够产生更为稀疏的分解效果,降低绝对平均误差、均方根误差计算值的统计区间。

摘要：采用一种具有很强的非平稳信号跟踪、预测能力的原子稀疏分解(ASD)法,作为人工神经网络(ANN)的前置分解方法,将风电功率序列分解为原子分量和残差分量,对原子分量进行自预测,残差分量进行ANN预测,再通过追加最新的风电功率实时数据来更新ASD的结果,进而滑动预测下一个时刻的风电功率。以实际风电场数据进行验证,结果证明了该模型可以有效地处理风电功率非平稳性,产生更为稀疏的分解效果,显著地降低了绝对平均误差、均方根误差计算值的统计区间。

原子稀疏分解 第2篇

一、信号的表示

一组线性独立的矢量φ={φ1}1∈T组成基, 它能张开整个空间。任意一个在空间中的矢量S, 可以通过线性组合来展开。

s在基矢量φ1上的展开系数是αi=〈s, φi〉。因为基具有独立性, 所以它是唯一的一种展开方式。如果φi⊥φj, i≠j, 则基ψ为空间S的一组正交基。式 (1.1) 可转化成如下矩阵:S=准α (1.2)

基矢量组成φ∈RM×N矩阵。如果N<M, 矢量空间S的有些矢量将不能完整的表示, 则基集合{αi}为非完备基;如果N=M, 矢量空间S的每个矢量都能完整表示, 而且表示方式是唯一的, 称基集合{αi}为一组完备基;如果M<N, 泽矢量空间S中的每个矢量都可以用αi线性组合展开, 且展开形式有无穷方式, 称基集合{αi}是超完备基。

信号也是采用一组矢量{ti}Ni=1来进行表示的。给定信号f∈RN, 通过基矢量{ti}Ni=1的线性组合展开。f在基矢量ti上的展开系数是f=∑Ni=1αiti=Tα, αi=〈f, ti〉, T为合成矩阵。目前常用的基有余弦基、傅里叶基和小波基, 这些基对应的合成矢量形成完备正交基。采用超完备的冗余基函数代替正交基函数来使表示形式多样化。此时, 基矢量ti组成一个超完备集{ti}Ki=1, 且N<K。因为ti不是线性相关的, 因此集合{ti}Ki=1是个字典, 不是个基。信号采用冗余字典对表示时, 合成矩阵D不再是方阵, 而是一个N×K维的矩阵, 且N≤K。设{ti}Ki=1对应合成矩阵D, 则信号f∈RN的超完备表示为:f=Dα (1.3)

给定f和D, 式 (1.3) 是一个欠定问题, 有无穷多个解α, 因此信号的稀疏表示是可行的。

二、MP算法

MP算法通过迭代方式从超完备元子库中选出与信号或者是残余信号最匹配的原子, 从而将信号表示为这些原子的线性组合, 实现信号的稀疏分解。设D={gm}0≤m<K是由K (K>N) 个原子构成的字典, 该字典包含K个线性无关的原子, 构成长度为N的信号空间的一组基。匹配跟踪算法对信号s进行分解的过程如下。设定初值R0s=s, 信号首先被分解为:

由于g0与R1s正交, 因此:

为了极小化||R1s||, 即让残余能量最小化, 必须选定原子g0, 使|<g0, R0s>|最大化。

用同样的方法, 得到:

经过M次迭代后, 信号s可以表示为:

其中RMs为M项近似残差, 并且满足:

三、一种基于互相关的快速匹配算法

(一) 算法思想

在信号稀疏表示中, 一般使用Gabor原子库。一个Gabor原子由尺度因子s、原子频率v、原子相位w和位移因子u共同决定, 且s、v、w三个参数决定原子形状, u只是决定原子的中心位置。如果对超完备原子库中的某一个原子gi (参数为si、vi、wi) 的选取方法保持不变, 信号长度为N, 取u=N/2, 通过平移就可以得到参数为 (参数为si、vi、wi、ui) 的原子 (ui≠N/2) 。为了不影响信号稀疏分解的效果, 尽量使ui的值在区间[0, N-1]上。这相当扩大了元子库中原子的数目。由于ui从0到N-1连续取值, 所有的N次内积运算<Rks, gk>都可以转化成两个列向量Rks和gk的一次互相关运算。内积运算的过程实际上是找互相关运算最大值的过程。如果提高了互相关运算的速度, 就提高了信号稀疏分解的速度。

(二) 算法描述

对两列长为N的数字信号, 其互相关函数为:

两列数字信号在互相关运算的结果附近存在着一个明显的峰值特征, 为了快速求出互相关运算的最大值, 可以间隔一定距离进行互相关运算, 确定最大值出现的区间, 然后再在这个区间内进行相关运算, 求出最大值和最大值点所在的位置。

假设每次间隔m个点进行计算, 式3.1可表达为:

将得到的个点互相关运算的值进行比较, 找出最大值 (此最大值可能不是实际的最大值) 。根据峰值特性, 最大值点应在区间[m-1, m+1]上。如果N比较大, 改进后互相关运算的计算量是原互相关运算的计算量的1/m。改进后互相关流程如图1。

四、实验结果与仿真

利用Matlab进行算法仿真, 取不同长度的信号和步进点数进行互相关运算, 求得的最大值与MP算法求得的内积最大值的误差, 算法速度与MP算法速度的比值。从实验结果可以看出, 当步进点数为2、4时, 误差很小, 当步进点数大于或等于6时, 误差开始变大。因此, 在实际运用中, 步进点数不能太大, 即使步进点数为1, 互相关算法比MP算法的速度要快。

本文针对信号稀疏分解中计算量大、分解速度缓慢等问题, 提出了一种信号稀疏分解的快速算法。它利用一种互相关的快速运算来代替稀疏分解中计算量巨大的内积运算, 提高了运算速度。通过理论仿真和实验表明, 在保证稀疏精度的情况下, 可以将信号稀疏分解的速度提高至普通MP算法速度的10倍以上。

参考文献

[1]张春梅, 尹忠科, 肖明霞.基于冗余字典的信号超完备表示与稀疏分解[J].科学通报, 2006, (6) :628-633.

[2]何虹丽, 嵇银辉, 等.基于差分进化算法的交通图像稀疏分解[J].电子元器件应用, 2011, (3) :35-37.

原子稀疏分解 第3篇

关键词：稀疏分解,快速计算,互相关运算

0 引言

近几年来,突破传统奈奎斯特采样定理的压缩感知(CS)理论框架得到了研究者们越来越多的关注。CS理论主要由3个部分组成:信号的稀疏分解;投影,观测;信号重构。由此可见,CS理论实现的前提就是要求被采样信号具有稀疏性,因此,对信号稀疏分解算法进行研究,有着极其重要的理论意义和广泛的应用价值。而其中以贪婪算法为核心的匹配追踪(MP)算法在计算复杂度和逼近效果方面要优于其他的算法,是目前研究者们最常用的算法。但是,MP算法的不足之处是基于MP的稀疏分解过程中,每一步都要计算残余信号在过完备原子库中每一个原子上的投影,这就造成了分解过程中相当巨大的计算量。高复杂度的计算量使得信号的稀疏分解在实际的研究应用中难以推广。针对这种情况,目前研究者们都在针对如何提高信号稀疏分解的速度、降低运算量这两个方面的问题进行研究。

本文提出一种基于互相关运算的信号稀疏分解的快速算法来改进信号的稀疏分解算法,该算法的思想是通过将基于MP的信号稀疏分解中花费大部分时间的内积运算转换成一次互相关运算,并且通过一种互相关运算的快速算法来对内积运算进行改进,这样一来,大大地提高了计算速度,在不增大计算误差的情况下缩短了计算时间,具有良好的应用前景。

1 匹配追踪算法的基本思想

我们假设待分解信号f的长度为N。其中,f∈H,H为有限维Hibert空间,D为过完备原子库(DH)。过完备原子库的形成方法引用参考资料[2],这种由Gabor原子组成的过完备原子库是颇具代表性的,它在多个文献中被研究者引用。原子gγ与待分解信号的长度相同并做归一化处理,即‖gγ‖=1。一个Gabor原子由一个经过调制的高斯窗函数来构建:

其中,g(t)=e-πt2是高斯窗函数,时频参数

MP算法是一种贪婪算法,该算法通过迭代从过完备原子库中选出与信号或者是残余信号最匹配的原子(也称为最佳原子),从而将信号表示为这些最佳原子的线性组合,即实现了信号的稀疏分解。其算法的具体流程图如图1所示。

信号的稀疏分解目前面临的主要问题就是分解速度十分缓慢,当信号长度比较长时,将需要花费大量的存储计算时间,而其中内积的运算占用了分解时间的一个很大部分。如何简化内积运算,提高运算速度和效率是亟待解决的关键问题。

由文献[1]可知,在稀疏分解的过程中,如果对原子库中的某一个原子(参数:si,vi,wi)的选取方法保持不变,我们可以取u=N/2,通过平移就可以得到参数为(si,vi,wi,ui)的原子(ui≠N/2)。但是为了不影响信号稀疏分解的效果,不妨让ui取所有可能的值[0,N-1]。这样一来,对于具有参数(si,vi,wi)的一个原子gγ,则这个原子要和残余信号作N次内积运算。由于ui从0到N-1连续取值,从理论上讲,所有的N次内积〈Rkf,gγk〉运算可以转换成两个列向量Rkf和gγk的一次互相关运算,内积运算的过程实际是找互相关运算的最大值的过程,如果我们能够提高互相关运算的速度,我们就提高了信号稀疏分解的速度。

2 一种互相关运算的快速方法

我们知道,对两列长为N的数字信号,其互相关函数可以用下式来表示:

由文献[5]可知,两列数字信号在互相关运算的结果附近存在着一个明显的峰值特性,那么如果我们需要得到互相关运算结果中的最大值,我们可以间隔一定点数进行互相关运算,也就是进行多个点的跳跃计算,这样我们就可以先找出最大值存在的一个比较小的区间,然后再在这个小区间内进行互相关运算,这样来找出最大值点位置和最大值,这样一来,就明显减少了计算的点数,加快了计算的速度。

假设我们取每次跳跃m个点进行计算,表达式将为:

将得到的个点进行比较,得到的最大值点可能不是实际的最大值点,根据它的峰值特性,实际的最大值点肯定就在这个最大值点左边第m-1个点到它右边第m+1个点共2 m-1个点中某个。如果N比较大的话,这样一来计算量就可以减小到原来的1/m。具体计算流程如图2所示。

3 实验结果与仿真分析

利用MATLAB进行算法仿真,将不同长度的信号和不同的步进点数的条件下得到互相关运算最大值与MP算法求得的内积最大值进行比较,得到的正确率如表1所示:

从表1可以看出,当步进点数为2,4时,基本没有误差或是误差很小,能够保证比较好的分解效果,而当步进点数大等于6时,误差开始变大,从而不能满足实际要求。因此实际运算中步进点数的选取可以参考计算精度和计算速度的要求进行合理选取。

试验中采用长度不同的实际数字信号,选取不同的步进点数的环境下进行仿真。由于信号稀疏分解的速度在不同的条件下(不同的计算机以及软硬件配置)上是不同的,所以给出信号稀疏分解实际的时间是没有多大参考价值的(几乎没有文献给出信号稀疏分解的时间)。所以我们以MP基本算法的速度为参照物,将其稀疏分解的速度设为1,表2中给出的是一个长度为400的信号用本文提出的算法进行稀疏分解的速度基于MP基本算法的倍数。

从表2可以看出,即使是步进点数为1时,本文的算法也比MP基本算法速度快。虽然从理论上来讲,所有N次内积运算转化成一次互相关运算的计算量几乎没有变化,但是计算效率却大大地提高了。而采用步进点数为4时,本文改进算法的速度在保证相对精度的情况下与MP基本算法相比得到了很大的提高,是MP基本算法的20多倍。

参考文献

[1]YIN Z K,WANG J Y,SHAO J.Sparse decomposition based onstructural properties of atom dictionary[J].Journal of SouthwestJiao tong University,2004(2).

[2]ARTHUR P L,PHILIPOS C L.Voiced/unvoiced speech discrimi-nation in noise using gabor atomic decomposition[A].Proc.OfIEEE ICASSP[C].Hong Kong,April,2003.

[3]KUSH R.VARSHNEY,MJDATETIN,JOHN W.FISHER,and ALAN S.WILLSKY.Sparse representation in structured dic-tionaries with application to synthetic aperture radar[J].IEEETransactions on Signal Processing,2008(8).

[4]ZHANG W Y.Study on matching pursuit based low-bit rate speechcoding[D].Beijing:Institute of Software,Chinese Academy of Sci-ence,2002

原子稀疏分解 第4篇

随着汽车数量的增 大,高速公路 收费站及 停车场增多,迫切需要智能收费系统来代替人工收费,减少人力开支。视频识别受环境因素影响较大,并且对摄像头位置有一定要求。

针对此问题,本文提出了采用基于发动机声音进行汽车识别的方法,从特征选取及识别机方面进行分析。针对识别实时性要求较高的特点,通过实验,采用多帧进行训练,采用单帧进行识别的方法可以减少运算量,提高运算速度并较好地保证识别准确率。

1发动机声音信号的 MP稀疏分解

1.1MP算法原理及分解方法

首先将原子参数离散化构成完备的原子库,然后通过信号与原子做内积迭代运算以寻找最优原子,通过反复迭代,达到设定阈值后退出循环,原始信号能量被一系列原子所取代。

(1)从原子库中寻找满足式(1)的原子,f(t)为原始信号,gr为原子:

(2)原始信号分解为内积和原子之积加上残余信号,如式(2),其中残余信号为:

(3)由式(2),根据式(3),对残差R1f(t)进行同样的迭代运算,如式(4)、式(5)所示。迭代过程如图1所示,直至残余信号小于设定阈值,分解结束。

(4)最终分解信号如式(6),其中grk为分解所得的时频原子,l为分解的原子个数,< Rkf(t),grk>为分解系数,Rnf(t)为残余信号。

1.2发动机声音信号分解与重构

针对发动机声音相对平稳的特点,采用256点为一帧可有效反映信号特征。Gabor原子是自适应信号分解中较常用的一种原子,在时频平面上的时频积很 小[1,2,3],能准确表示时频信号的细节特征。因此,选用Gabor原子形成原子库,采用20个原子进行分解重构,图2为原信号,图3为重构信号,可以看出重构信号可较好地反映原信号的特征[4]。

2发动机声音识别

2.1MP稀疏分解在发动机声音识别中的应用

本文所研究的发动机信号相对平稳,采用Gabor原子能有效反应发动机信号特征。Gabor原子公式见式(7)。

其中g(t)为高斯窗,发动机声音分解过程中,每次迭代都会得到相应的原子grk,最终信号分解为式(6)中20个原子组合的形式,原子由时 频参数γ = (s,u,v,w)表示,其中s反映原子尺度、u反映位移、v和w反映原子频率。时频参数γ反映发动机信号的重要特征,因此选用发动机声音信号稀疏分解所得参数γ作为发动机声音信号特征进行识别。

2.2BP神经网络在发动机声音识别中的应用

本文研究的发动机声音信号相对平稳,采用BP(BackPropagation)神经网络作为训练识别网络。BP神经网络实现原理如图4所示。

实现方法为:

(1)初始化各层网络权值为较小的随机数,依次输入P个学习样本[5]。设当前输入为第P个样本,依次计算各层输出。

(2)求各层反传误差δ,期望输出d ,实际的输出y,输出层的节点数m ,中间层节点数n1[6],输入层节点数n2,如下式[5]:

(3)记录已学习过的样本个数P。如果p<P,转到步骤(1)继续计算[6],如果p=P,按权值修正公式修正各层的权值ω。η为学习速率[5]。为加快收敛,引入动量项α。

(4)按新的权值计算,若对每个P和l都满足|dl(p)yl(p)|<ε,ε为一个给定正小数 ,或达到最大训练次数,则终止学习[5],否则进行新一轮学习。

3识别结果分析

本实验采用汽车收费站附近受环境干扰较少的机车发动机声音作为训练和识别样本。由于发动机声音信号在频域内相对平稳,本实验采用1024点为一帧进行实验。信号MP稀疏分解时采用Gabor原子形成原子库进行分解[7],分解为20原子,采用时频参数γ作为发动机声音信号特征进行识别[6]。

进行BP训练时采用matlab进行模拟仿真,将车型分为轿车、面包车、客车、重型汽车4类,分别采用100辆作为训练集,100辆作为识 别集,共800辆车作为 样本[8]。采用期望最小均方误差为0.001,以时频参数γ= (s,u,v,w)作为特征参数,输入采用4节点,4维时频参数s、u、v、w作为特征参数进行输入,输出为4个节点,分别设定(0、0、0、1)作为轿车输出,(0、0、1、0)作为面包车输出,(0、1、0、0)作为客车输出,(1、0、0、0)作为重型汽车输出[6],隐层为50个,动量因子为0.85,学习速率为0.1。

由于发动机声音识别实时性要求较高,因此多帧进行识别与采用单帧识别在时间上有一定影响。表1为单帧识别和3帧识别对比,可以看出,单帧识别在时间上有较大优势,识别效果相差不大,因此实际应用中选择单帧识别。

4结语

本文针对发动机声音相对平稳的特点,采用MP稀疏分解方法,用Gabor原子作为完备原子库分解发动机声音信号9,通过分解重构,反映原始信号特征[9]。

原子稀疏分解 第5篇

大型变带宽正定对称矩阵的分解具有重要的工程应用价值,局部分解在现代科技生产研究中有很大的用处,合理利用已有矩阵分解结果,对变动矩阵采用局部分解算法,可以缩短整体矩阵的分解时间,提高矩阵分解效率。

1 矩阵局部分解的原理

Cholesky分解法是求解对称正定线性方程组最常用的方法之一。

定理:若A∈Rn×n对称正定,则存在一个对角元为正数的下三角矩阵L∈Rn×n,使得A=LLT成立,其中L∈(lij)n×n是下三角矩阵。即:

从式(1)可以看出,对角线值aii的分解结果lii只依赖于第i行的已分解值,而第i列的其它值aki的分解值lki依赖于k行的已有分解值以及第i行分解值。模型的矩阵信息重用与局部化分解原理如图1所示,在矩阵局部分解中,下三角矩阵某一列的矩阵分解值仅与该列的对角线上元素所在行以及该行之后所有行的分解结果有关。根据式(1),可以推导信息重用的局部分解原理。

本文基于cholesky分解法,设计并实现了大型稀疏矩阵的局部分解算法,通过列优先的方式进行分解,矩阵某元素变化可以归结到其列对应得对角线元素。通过重用在此点分解之前的分解结果,减少整个矩阵分解产生的时间代价。

2 矩阵读取与分解结果的存储

在对矩阵进行cholesky分解时虽然只需要读取一半的数据,但是对于大型的稀疏矩阵,仍然需要对内存进行优化。待分解的原矩阵存放在txt文件中,然后利用文件指针来进行矩阵数据的读取操作。

在计算分解时,首先要读取文件中的矩阵数据,若用二维数组把矩阵的数据一次遍历完完全放在其中,则会消耗大量的内存空间,此时可以一边打开文件、读取文件中的数据一边计算,再将分解的结果写入另一文件中。分解的结果存储在文件时,可以进行压缩存储,利用三元组的方法将非零结果和列号、行号存储在文件中。利用此种方法可以解决大型矩阵数据占用大量内存的问题。

在分解时同时打开两个txt文件,一个读取原矩阵数据另一个存储分解的结果,这样既节省内存又可以在局部分解时方便地利用先前已经分解的结果进行局部分解计算,优化了大型稀疏矩阵计算时内存的问题。

3 分解结果的重用分析

对于图1变化值,根据Cholesky分解列优先的分解原则,可知道此点变化只会影响下面阴影部分三角形的值变化,而不会影响前面值的变化,所以只重新分解矩阵元素(图1的阴影部分)。即从矩阵元素开始有变化的k值起用Cholesky分解公式计算一遍。但下面小阴影部分三角形重新计算要重用前面已经分分解过梯形部分的数据,所以当局部变化时,不再需要进行整体矩阵的重新分解计算。

只考虑对角线上元素值的变化,对于特殊点的变化可以归结为对角线上元素值的变化,因为列行所对应的对角线上的值变化会影响此下面一列的变化,其实质还是一样,都可以从其列行对应的对角线上的点为起始点开始进行局部分解。

信息重用的关键代码如下:

其中,num是对角线改变点行/列号。由代码可以看出,变化点出现时,是从变化点处开始局部计算的,而不是整体矩阵的分解。

4 局部分解与整体分解的效率对比

以100阶矩阵的局部分解为例,设矩阵改变元素所在的行对角线位置为30~100之间,则整体矩阵分解和局部分解运算速度如表1所示。

由表1可见,整体分解,由于每次分解次数相同,除去改变点数值变化引起的运算复杂度增加或降低而导致分解时间少量变化外,总体分解时间在(720~860)×10-6s之间,变化不大;重用分解结果数据的局部分解,分解时间在变化点行/列数较小时,运算时间接近整体分解,随着变化点行/列数变大,由于重用的运算结果数据越来越多,本身需要运算的数据越来越少,效率大幅提升。最显著就是100阶的元素值改变时,整体分解需要重算99阶,而局部分解只需要分解1阶数据。

5 结语

提出了一种能在分解大型变带宽正定稀疏矩阵中,局部点发生时变化利用分解结果数据快速分解的方法。与原始整体分解算法相比,局部分解大大提高了分解效率,为分解百万千万阶矩阵节约了大量时间,这种信息重用的方法,为以后实际应用到网格化模型刚度变化快速响应提供了可能。不足之处在于内存的优化,对于大量原数据,边读边算边存储,让内存持续保持较低占用率。

后续的研究中,一方面要继续简化局部分解,深入研究信息重用和内存优化的方法;另一方面,要利用局部分解的高效性,进行网格化模型局部变化的仿真模拟。

参考文献

[1]唐建国,曹魏娟,高曙明,等.面向CAE的简化模型误差评估与信息重用[J].中国机械工程,2012,23(8):961-967

[2]阮百尧,熊彬.大型对称变带宽方程组的Cholesky分解法[J].物探化探计算技术,2000,22(4):361-363,368

[3]谢晓峰,李代平,陈璟华,等.大型稀疏线性方程组的一种压缩求解算法[J].计算机工程与应用,2001,37(5):110-111

[4]陈建平,Jerzy Wasniewski.Cholesky分解递归算法与改进[J].计算机研究与发展,2001,38(8):923-926

[5]袁晖坪.行(列)对称矩阵的LDU分解与Cholesky分解[J].华侨大学学报(自然科学版),2007,28(1):88-91

[6]刘巧华,魏木生.LU和Cholesky分解的向前舍入误差分析[J].高等学校计算数学学报,2006,28(4):358-366

[7]宋滔,王绪本.稀疏矩阵快速回代的Cholesky分解法[J].物探化探计算技术,2013,35(3):293-296

原子稀疏分解 第6篇

局部放电(Partial Discharge ,PD)是电气设备绝缘故障的早期表现形式,是评定电气设备绝缘状况的重要指标之一,倍受重视[1,2]。由于局部放电信号微弱,而电气设备运行环境中普遍存在载波通信和高频保护等窄带干扰信号,同时存在各种随机干扰信号,这些干扰信号易将局部放电信号淹没,以致于不易提取出真实、有效的局部放电信号。

国内外学者在局部放电信号的抗噪、滤波和干扰抑制方面做了大量而深入的研究,将先进的信号处理方法应用于抑制局部放电干扰信号中,主要的方法有基于傅里叶变换的数字滤波法[3,4]、基于小波变换的滤波法[5,6]和基于经验模 态分解 (Empirical Mode Decomposition,EMD)的时空滤波法[7,8,9]以及相关的组合算法。在这些方法中,基于傅里叶变换的滤波法受干扰信号的中心频率影响较大,自适应能力不足;基于小波变换的滤波法能抑制窄带干扰信号并去除白噪声,但存在频带混叠、泄漏效应、原始信号能量损失严重、稳定性差、参数较难统一等缺陷;基于EMD的时空滤波法在模态分解过程中容易出现虚假频率,且存在端点效应。参考文献[10] 提出了最优离散谐波小波包变换法,在抑制混频随机周期性窄带干扰方面取得了较好的效果。

Selesnick提出的信号共振稀疏分解方法[11,12,13,14,15,16], 在国内已成功应用于滚珠轴承的故障诊断[17,18,19],取得了良好的检测效果,但该方法在电气设备故障诊断方面的应用尚未见报道。

本文提出了基于共振稀疏分解的电气设备局部放电信号窄带干扰抑制新方法,该方法根据局部放电信号和窄带干扰信号的特点,分析了品质因子Q、 冗余因子r、分解层数L和权重系数λ 对共振稀疏分解的影响;针对局部放电信号的特征,通过合理选择分解参数,增强了小波基选择的灵活性,实现了对含窄带干扰的局部放电信号的稀疏分解,将其分解为高共振分量、低共振分量和残余分量,它们分别与窄带信号、局部放电信号和随机干扰相对应,从而提取出局部放电信号。分别对仿真数据和试验数据进行共振稀疏分解,并与小波软阈值滤波法进行对比, 处理结果验证了本文方法的有效性。

1共振稀疏分解理论

信号共振稀疏分解方法的基本思想是用2种品质因子不同的小波基函数稀疏表示复杂的信号,主要包括品质因子可调的小波变换和信号的稀疏分解2个部分。

1.1品质因子Q可调的小波变换

品质因子是表征信号的频率聚集程度的量度, 品质因子越高,信号的频率越集中,时域上振荡次数越多。可调品质因子小波变换[10,11]通过直接指定品质因子Q(定义为中心频率与频率带宽的比值)和冗余因子r来设计小波,增强了品质因子选择的灵活性。

可调品质因子小波变换由2个不同的滤波库 (高共振分量滤波库、低共振分量滤波库)组成,具体变换流程如图1所示。

图1中,x(n)为待分析信号;H0(ω)为低通滤波频响函数;H1(ω)为高通滤波频响函数;α,β分别为低通滤波库和高通滤波库的尺度参数;x1(n)为提取出的高共振分量;x2(n)为提取出的低共振分量。

带通滤波的品质因子Q定义为中心频率fc和带宽BW的比值,即Q=fc/BW。当品质因子Q和冗余因子r的值选定后,对应的滤波库尺度参数α 和β 由式(1)确定 。

品质因子Q的小波变换的3层分解如图2所示,第1层分解的输入信号为x(n),第2层分解的输入信号为第1层分解输出的低通滤波信号CA1, 第3层分解的输入信号为第2层分解输出的低通滤波信号CA2;CD1,CD2,CD3分别为各层输出的高通滤波信号;CA3为第3层分解输出的低通滤波信号。

根据信号自身特征自适应调节Q,使得小波基的持续振荡程度与待分析信号振荡程度匹配良好,从而达到降噪的目的,增强信号的稀疏表示效果。

在可调品质因子小波变换中,由确定的Q和r生成的小波基函数库中,所有的小波基函数都具有相同的品质因子Q。Q=3,r=3,L=12时的各级小波时域波形和频率响应曲线如图3所示。

在同一个小波基函数库中,不同级数(层数)的小波的中心频率fc和带宽BW不同。一个信号中由第j级小波表示的部分称为该信号的第j个子带。每个子带的中心频率fc和带宽BW由式(2)、 式(3)计算。

式中:fs为待分解信号的采样频率;j为小波所在的级数。

1.2高共振分量与低共振分量的分离

信号共振稀疏分解方法利用形态分量分析[10,11]将信号中各分量按振荡特性进行非线性分离,建立起高共振分量和低共振分量各自的最佳稀疏表示形式。

由品质因子可调的小波变换得到2种品质因子为Q1和Q2的小波基函数库s1和s2后,将待分解信号x表示为

式中:w1,w2为小波变换系数;w1s1,w2s2分别为待分解信号x中的高共振分量和低共振分量;n为待分解信号y中不能用w1s1和w2s2表示的部分,称为残余分量。

如果系数矩阵w1和w2中元素数值极小的个数多,且残余分量n很小,则式(4)是待分解信号y的稀疏表示。

形态分量分析的目的是从待分解信号x中分离出源信号,为了达到将高共振分量与低共振分量有效分离的目的,形态分量分析的目标函数可表示为

式(5)中的权重系数λ1和λ2决定了最优系数w1*和w2*,进而决定了最终分解结果中高共振分量和低共振分量的能量。

式(5)综合考虑了式(4)中残余分量n的大小以及w1,w2小波变换系数的稀疏程度。目标函数的值越小,表示分解越稀疏。

利用分裂 增广拉格 朗日收缩 算法[12,13](SALSA)获得最优的w1和w2,使得目标函数的值最小,稀疏分解的过程就是寻找能使该函数达到极小值时的w1*和w2*的过程,求取出的高共振分量和低共振分量的估计值分别表示为

1.3共振稀疏分解的实现步骤

(1)根据实际信号,选择合适的品质因子Q,通过品质因子可调小波获得相应的小波基函数库。

(2)建立目标函数,用SALSA收缩算法获取最优的系数w1*和w2*,使目标函数值最小。

(3)由w1*和w2*对原信号进行重构,分别获取原信号中的高共振分量

2局部放电信号分析

2.1仿真分析

现场检测到的局部放电信号为衰减振荡型信号,可用单指数和双指数衰减振荡来模拟。

单指数衰减振荡波形表达式如式(7)所示:

双指数衰减振荡波形表达式如式(8)所示:

式中 : A为信号幅值 ; τ 为衰减系数 。

局部放电仿真信号参数见表1 。

现场监测结果表明,局部放电信号受窄带干扰及白噪声干扰影响较为严重。为模拟现场真实局部放电信号,给理想的局部放电信号叠加窄带信号和随机白噪声,其中窄带干扰数学表达式如式(9)所示,所加窄带干扰的频率为150 ,500kHz和1.75, 2.5 MHz,且幅值均为1mV;随机白噪声满足高斯分布N(0,0.022)。

理想的局部放电信号s(t)和被窄带、随机白噪声污染的局部放 电信号x(t)波形如图4所示,显然,图4(b)中的理想局部放电信号s(t)已完全被窄带干扰信号所淹没。对s(t)和x(t)进行频谱分析, 幅值谱如图5所示。将图5(a)、(b)两图对比可知, 理想局部放电信号的幅值小,而含干扰的局部放电信号的最大幅值约为理想局部放电信号最大幅值的40倍,且呈离散谱线。

对信号x(t)进行共振稀疏分解,根据信号特点和实际经验[15,16,17]设置品质因子Q、冗余因子r、分解层数L和权重系数λ,参数见表2。

将x(t)分解为高共振分量x1(t)和低共振分量x2(t)以及残余分量n(t),各分量的时域波形如图6所示。由图6可知,高共振分量以窄带干扰信号为主,低共振分量以衰减振荡信号为主,残余分量中包含随机白噪声。

对高共振分量x1(t)和低共振分量x2(t)进行幅值谱分析,幅值谱如图7所示。由图7可知,高共振分量的频率与所加窄带干扰信号的频率分量基本一致,说明分解得到的高共振信号主要是窄带干扰信号;低共振分量的幅值谱与理想的局部放电幅值谱基本吻合,说明分解得到的低共振信号中主要是局部放电信号。

将含干扰的局部放电信号x(t)用Db44小波分解,并进行软阈值滤波,滤波后重构得到的局部放电信号如图8所示。

为比较本文方法与小波软阈值滤波法抑制窄带干扰的能力,以信噪比SNR、均方根误差e和互相关系数c作为评价指标,3个指标的定义为

式中:s(i)和s^(i)分别为理想和消噪后的局部放电信号;N为信号长度。

本文方法与小波软阈值滤波法抑制窄带干扰的评价指标见表3。由表3可见,小波软阈值滤波法提高了局部放电信号的信噪比,而本文方法提高信噪比的能 力更强、均方根误 差更小、波形相似 度更高。

在运用共振稀疏分解法时,调节高共振和低共振滤波库的品质因子,使高共振滤波库小波基的时域振荡特征与窄带信号的特征相匹配,且Q越大, 各级(层)小波的频率响应重叠度越小,分辨率越高;使共振滤波库小波基的时域振荡特征与局部放电信号的特征相匹配,提取出局部放电信号。减小冗余因子r,则降低各级小波频率响应间的重叠度,可以区分不同频率的信号;分解层数L越大,信号分解的计算 量越大;在最优稀 疏分解的 目标函数 式 (式 (5))中,权重系数λ 越大,分解越稀疏,同时使待分解信号分解到高共振和低共振分量中的成分减少,而分解到残余分量中的可能性增大,所以权重系数λ的取值需要综合考虑。

2.2实测局部放电信号的信号分析

用本文方法对实测的某发电机中性点局部放电信号进行处理,并与小波软阈值滤波法对比。图9为实测信号的波形和幅值谱。图10、图11分别为小波软阈值滤波法和本文方法的处理结果。

由图9可知,实测局部放电信号已被窄带和随机干扰信号淹没。对比图10和图11可知,经小波软阈值滤波和小波重构后得到的局部放电信号波形与理 想局部放 电信号所 具有的衰 减振荡特 征 (图4(a))相比,有较为明显的失真;而用本文方法提取的局部放 电信号失 真较小,且幅值谱 峰值在2.5 MHz处,与理想局部放电信号幅值谱(图5(a)) 特征一致。可见,采用共振稀疏分解法可有效消除窄带干扰,并保留局部放电信号的主要特征。

3结语

原子稀疏分解 第7篇

心电信号 (ECG) 是心脏电生理活动的反映, 它蕴涵着关于心脏节律及其电传导活动的丰富生理病理信息, 可为心脏类疾病诊断、监护和评价提供重要的数据参考, 特别是因其无侵入、成本低、携带方便、可靠性较高等优点, 受到近年兴起的远程医疗、移动健康等领域的重视[1,2]。信号稀疏分解理论为心电信号的波形检测与识别提供了新的思路[3,4,5]。文[6]将该思想应用于ECG波形的特征点检测及识别, QRS波群检测率宣称达到99%, 但P波、T波检测率仅为85%。

本文在稀疏分解理论的基础上, 提高了算法对任意输入的ECG片段的自适应能力和P波、T波特征点的正确识别率。

1 稀疏分解基本原理

稀疏分解理论认为很多现实的信号可被分解为少数基本信号。匹配追踪算法 (MP) 可用来寻找目标信号基于特定字典的稀疏表示, 其基本思路是在每一次的迭代过程中, 从过完备字典D中选择与信号最为匹配的原子来构建稀疏逼近, 并且求出信号表示残差, 即可得到信号的稀疏表示之一。

本文应用该思想和MP分解, 通过将目标ECG分解为若干典型的特征波分量, 可实现更准确的特征点检测与识别。

2 冗余原子集的构建

一个完整周期的ECG信号可以表示为P波、一个QRS波群和一个T波的组合。因此, 在原子集中, 至少应包含表示P波、QRS波群和T波自身特征的成分。下面根据心电信号各特征波的特点, 构造出P波、QRS波群、T波的子超完备字典。

2.1 P原子

P波的幅值比较小, 波形圆钝, 时间宽度不大于0.11秒, 幅值电压 (高度) 小于等于0.25毫伏, 如图1所示三个典型的离散序列。对其按公式gγ=gγ/norm (gγ) 进行归一化处理。

对这三个离散序列通过在坐标系上进行拉伸、平移变换, 构造出P波的子超完备字典。令s为尺度因子, u为平移因子, 将原始离散数列按照s做拉伸变换可以获得不同波宽的P波, 按照u做平移变换可以获得在横轴不同位置的P波, 从而得到P波的子超完备字典。

2.2 QRS原子

对于正常成年人, 其心电图中的QRS复合波时间宽度为0.06~0.10s。Q波的幅值变化量一般不超过0.3mv。正常S波在标准导联, 其幅值变化量应在0.6mv以内。

根据QRS波的特征给出图2所示的离散序列。将归一化后的离散序列根据尺度因子s和平移因子u在坐标系上进行拉伸、平移变换, 构造出QRS波的子超完备字典。

2.3 T原子

在心电图QRS波之后的一个圆钝波形就是反应心室快速除极过程的T波。其时间宽度0.05~0.25s。T波的幅值与R波幅值的比例应大于1/10。

根据T波特征, 提取到图3所示的离散序列, 并对其进行归一化处理。将归一化后的离散序列根据尺度因子s和平移因子u在坐标系上进行拉伸、平移变换, 构造出T波的子超完备字典。

3 基于改进MP的特征点提取与辨识算法

设原始信号为f, 长度为N, 本算法所得特征点在原始信号中的位置为t。和为原始信号首尾可能存在的不完整特征波的匹配原子。P、QRS、T波对应的原子集分别为DP、DQRS、DT, 原子集中特征点的位置为。该算法的具体步骤为:

1) 基于P、QRS、T波形特征生成原子集DP、DQRS、DT;

2) 对原始信号做稀疏分解和重构, 得到匹配原子gbest, , 其中gγ∈DP、DQRS、DT;

4) 记录匹配原子gbest对应的稀疏系数a, 尺度因子向量s, 平移因子向量u;

5) 判断找到的是否对应了原始信号中的某一个特征波, 如果稀疏系数a>0.8*sum (a_before) , 则表示找到一个原始信号中的特征波g;

7) 根据4) 和6) 中获得的匹配原子对应的参数a、s、u, 计算原始信号f特征点的位置, t=t′*s+u, 特征点的类型由对应原子的类型直接获得。

4 实验结果及分析

4.1 数据来源说明

NIF ECG数据库是Physionet数据库中的一个, 全称是Non-Invasive Fetal Electrocardiogram Database。从NIF ECG数据库提供的55组心电数据中随机选取20组胸导联ECG信号, 并从每组信号中随机选取10段长为2s的片段进行测试。

4.2 评价指标

该算法检测得到的特征点是否正确可以通过与原始ECG的人工标定信息相比较而得出。设某个特征点标定的位置信息为, 算法测得的该特征点的位置信息为, 令该算法在某个特征点处的检测时间误差为τ, 则

令τ0=10毫秒, 若某个特征点处的检测时间误差为τ≤τ0, 则认为该特征点被正确识别。

我们通过该算法识别特征点的灵敏度 (Detection Sensitivity) 和预测力 (Positive Predicitivity) 来评价该算法, 其中灵敏度和预测力的计算公式分别为:

其中TP表示算法正确识别的特征点的个数, FN表示算法没能识别到的特征点的个数, FP表示算法错误识别的特征点的个数。

4.3 测试与结果讨论

经统计, QRS波群、T波和S波峰值点的识别灵敏度分别达到了98.5%、93.2%和83.8%, 具体统计结果如表1所示。

图4展示了NIF ECG数据库中“ecgca595”的特征点识别。

图5给出了NIF ECG数据库中“ecgca595”重构后的ECG图。

从以上的图形可以看出, 即使输入的ECG信号也包含了不完整心拍, 本文的算法也能进行正确重构, 并准确地得到特征波的位置信息。

5 结论

本文算法经验证能够以98%以上的准确率识别ECG信号的QRS波群位置、幅度等信息。在稀疏分解原子集构建中, 考虑到不一定能找到合适的生成函数来表示ECG信号的特征波, 而若采用训练的方法则会增加算法复杂度, 无法保证算法的实时性, 故本文中采用直接从ECG信号中获取离散序列的方法组成原子集, 亦获得了较好的效果。在基于原子集的信号分解和重构过程中, 对于P波、QRS波群和T波子原子集的重构顺序没有特别的要求。即P波和T波的识别完全根据ECG信号波形自身的特点来完成, 不依赖于QRS波群的识别, 和传统算法相比大大提高了P波和T波的识别率。

参考文献

[1]Cuiwei L, Chongxun Z, Changfeng T.Detection of ECG characteristic points using wavelet transforms[J].Biomedical Engineering, IEEE Transactions on, 1995, 42 (1) :21-28.

[2]Sayadi O, Shamsollahi M B.Model-Based Fiducial Points Extraction for Baseline Wandered Electrocardiograms[J].Biomedical Engineering, IEEE Transactions on, 2008, 55 (1) :347-351.

[3]Babaie-Zadeh M, Mehrdad B, Giannakis G B.Weighted sparse signal decomposition[C].Kyoto:2012.3425-3428.

[4]Shen Y, Ye C, Li H, et al.An MP sparse decomposition algorithm with2-step search method based on partitioning the over-complete dictionary[C].Ningbo:2011.2539-2542.

[5]Ahani S, Ghaemmaghami S.Image steganography based on sparse decomposition in wavelet space[C].Beijing:2010.632-637.