矩阵奇异值（共7篇）
矩阵奇异值 篇1
图像处理的内容十分广泛, 具体而言, 可以分为:图像获取、图像增强、图像复原、图像压缩、图像分割等。这些内容都是基于矩阵的处理得到的。图像处理领域中, 物体的形状复原技术应用也十分广泛。在模式识别、三维技术、物体识别等方面都有十分广泛的应用。采用经验的照度差的方法固然可以获得物体形状复原的信息, 但是往往需要一个和物体本身反射率以及物质特征一致的参照球才能解决问题。但也有一些不利用参照物, 直接从物体本身就能进行形状复原的办法。这些方法都需要在形状复原的时候使用多个光源来完成。本文探讨奇异值分解法来解决形状复原的问题。
1 奇异值分解
数字图像处理当中:一幅图像可以表示成为一个二维函数, 也就是空间笛卡尔坐标, 坐标的值称为该点图像的灰度。模拟图像是不方便用计算机处理, 利用取样和量化把模拟图像转换为数字图像是比较可取的方法。这样数字图像的结果是一个实际的矩阵。
更一般的矩阵表达方式为:
奇异值分解 (singular value decomposition) 是一种矩阵的分解方法, 奇异值分解是现代数值的最基本和最重要的工具之一。简单地说, 我们可以把m行n列的数字图像A分解成如公式 (1) 这样的表达式:
其中, U是n行n列的矩阵。Σ是n行m列的矩阵, 它除了对角线的成分之外, 其余都是0, 对角线的数据是非负数。V*是m行m列的V矩阵的伴随矩阵。
我们具体在做形状复原的时候, 往往只需要获得物体表面的法线向量矩阵, 也就是数字图像中每个像素点对应的向量x, y, z就可以得到物体的形状信息。
2 使用偏光镜分离镜面反射成分
数字图像实际是光的一种数字表现, 众所周知, 光是一种电磁波, 分为横和纵两个方向。偏光镜片的性质之一就是只允许一个方向的光通过透镜传送。从观测角度看的话, 如果相机镜头前安装了偏光镜, 随着偏光镜片的旋转, 观测值也就是摄影得到的数字图像的值也会有相应的变化。我们在做数字图像处理的时候, 往往镜面反射强烈的地方不利于物体的观测和形状复原。
实际上得到的数字图像的灰度值如公式 (2) 所示, 它分为扩散反射成分Idif和镜面反射成分Isc和Isv。其中, f (θ) 是一个系数, θ是偏光镜片旋转的角度。可见, 旋转角度θ变化, 函数f (θ) 也发生变化, 从而引起f (θ) 和扩散反射成分Isv的乘积也发生变化。
如果, 把Idif和Isc表示成I′的话, 公式 (2) 可以简单地表示成:
把公式 (3) 变成
这样的形式的话, 也就是说, 如果把Is看做I′和Isv是一个相对应n列的矩阵。A是1和f (θ) , 是一个相对应n行的矩阵。利用公式 (4) 可以把其中一部分镜面反射成分Isv从矩阵中分离出来。
在具体做实验的时候, 固定一个光源下拍摄到的数字图像往往无法避免镜面反射的问题, 特别对于一些镜面反射特性很强的物体来说尤为如此。这就需要利用奇异值分解法分解出镜面反射成分, 具体一副图像无法完成矩阵的构成, 我们需要每隔一定角度旋转偏光镜片拍摄照片, 从而得到数枚图像构成的矩阵来提供分解。
3 数字图像的奇异值分解
一副数字图像一共有n个点, 其中我们一共有m张图像, 也就是可以组成一个 (n×m) 灰度值矩阵I。I可以表示成 (n×n) 的反射率矩阵R, (n×3) 的法线向量矩阵N, (3×m) 的光源方向矩阵S, 和 (m×m) 的光源照度矩阵L。如公式 (5) 所示。
如果把反射率矩阵R和法线矩阵N一起表示成面特性矩阵N, 光源照度矩阵L和光源方向矩阵S一起组成光源特性矩阵S。这样公式就可以简单化的表示成为:I=NS这样的形式。
利用公式 (1) 和公式 (5) 可以把数字图像灰度值I像I=U∑V这样利用奇异值分解法分解。这里, 如果U的前3列表示为U′, ∑是一个3行3列的∑′, V的前3列表示成V′的话。就会变成I′=U′∑′V′这样的公式。利用面特性矩阵N和光源特性矩阵S, 可以做公式 (6) 那样的奇异值分解。
其中, ∑可以通过最小二乘法解出来。实际在图像分解的时候, 因为光源的不确定性, 利用具体光源方向可以判断出具体的向量矩阵。
这里, 从不同光源角度观测的图像都会由前一章的奇异值分解法分解出镜面反射成分, 去除镜面反射成分的灰度图像可以被利用在本部分的奇异值分解法当中。这时候也需要多个光源角度下拍摄数组图片, 从而得到丰富的矩阵数据来提供分解。
4 实验结果
为了获得大量的实验数据, 具体分别使用16个方向的光源。也就是在16个不同角度的光源下, 分别拍摄16组数字图像。每组图像都采用偏光镜拍摄, 其中一组图像一共有16张, 每张对应偏光镜片的一个角度, 也就是从0到360度, 偏光镜片每旋转15度拍摄一组图像。每组采用偏光镜拍摄的图像16枚, 这16枚图像可以利用第二章的奇异值分解法将镜面反射成分从数字灰度图像中分解出来。也就是说1组光源一共得到1枚没有镜面反射的图像。将利用奇异值分解法分解的16组光源一共16枚没有镜面反射的数字图像, 利用第三章的奇异值分解公式, 把物体表面向量分别出来, 从而获得形状复原的结果。如图2和3所示。图1是目标物体原图之一。图2是物体x和y方向的向量, 图3是z方向的向量。
在这里, 其实有些地方有光泽的镜面反射成分并没有完全去除, 这是因为Isc这一部分的镜面反射成分还可以透过偏光镜片过滤。从而导致奇异值分解法也无法去除这部分数据, 也给之后的计算带来了一些误差。
参考文献
[1]周波, 陈健.基于奇异值分解的、抗几何失真的数字水印算法[J].中国图象图形学报A, 2004 (04) :506-512.[2]李海峰, 王树勋, 温泉.基于SVD和ICA的鲁棒水印算法[J].中山大学学报 (自然科学版) , 2004 (zk) :53-57.[3]陈永红, 黄席樾.基于混沌映射和矩阵奇异分解的公开数字水印技术[J].计算机仿真, 2005 (01) :138-141.
矩阵奇异值 篇2
潜在语义分析LSA(Latent Semantic Analysis)是一种用于知识获取和展示的计算理论和方法,它使用统计的方法对大量的文本集进行分析,从而提取和表示出词语的语义,通过这种语义能克服传统的信息检索中使用词汇匹配带来的两个问题:词语的歧义性和同义性。
LSA需要通过数学的方法建立潜在语义空间模型,矩阵奇异值分解SVD(singular value decomposition)是目前普遍使用的LSA空间的构造方法。SVD实现中存在的一个关键问题就是词汇—文本矩阵的SVD计算量。本文主要针对潜在语义分析中SVD计算的特点设计并实现了一种并行算法。
目前矩阵SVD计算的方法主要分为两类:基于矩阵QR分解的方法和基于Jacobi迭代的方法。基于矩阵QR分解的算法比基于Jacobi迭代的算法计算速度更快,但是基于Jacobi迭代的算法计算得到的结果精度更高[3]。文献[4]采用双边Jacobi算法实现了奇异值分解的并行化。由于单边Jacobi算法具有更高的效率,而且精确性也不比双边Jacobi算法差[5],所以本文采用单边Jacobi算法实现了奇异值分解的并行化。
1 潜在语义空间及Hestenes算法简介
1.1 潜在语义空间模型
定理1(奇异值分解定理) 任何一个矩阵X(m*n),设X的秩记为r,均可分解为两个正交矩阵和一个对角矩阵的乘积:X=TSDT 其中T(m*r)=(t1,t2,…,tr)为正交矩阵,其中t1,t2,…,tr被称为X的左奇异向量;S(r*r)=diag(σ1,σ2,…,σr)为对角矩阵,σ1,σ2,…,σr为X的所有奇异值,并且满足如下的关系σ1≥σ2≥…≥σr>0;D(n*r)=(d1,d2,…,dr)为正交矩阵,其中d1,d2,…,dr为X的右奇异向量。
1.2 Hestense算法介绍
Hestense算法是一种利用一系列Jacobi平面旋转矩阵在计算机上实现矩阵奇异值分解的方法。该算法的基本思想是构造一系列的Jacobi平面旋转矩阵V(θ1),V(θ2),…,V(θm),使得矩阵A变为正交矩阵H。即
A V(θ1)V(θ2),…,V(θm)=H
式中假设需经过S次变换才能使m*n维实数矩阵A的任意两列正交。其中每一个Jacobi平面旋转矩阵V(θi)使乘积A V(θi)中变换的两个列向量正交,即,
其中Jacobi平面旋转矩阵V(θi)是把一个n*n的单位矩阵中的第i行,第j行,第i列,第j列交叉的四个元素用c和s表示。c=cosθ,s=sinθ,ai和aj表示矩阵A的第i列和第j列。我们选择θ使得a
一般称n(n-1)/2次Jacobi变换为一次扫描,每一次扫描使得矩阵A中的任意两列都恰好正交一次。经过有限次的扫描,可以使矩阵A的任意两列都正交,即变为一个正交矩阵。
1.3 Hestense算法的改进
本文对Hestense算法进行了以下两方面的改进,使得可以在矩阵奇异值分解的同时对矩阵列向量的模有序排列。这样可以在SVD计算结束时,直接获得有序排列的奇异值。
① 当每个Jacobi平面旋转矩阵使矩阵A的两列正交的同时对矩阵两列的模进行排序;
② 当矩阵两列的内积相对于所有列向量对的平均内积较小的时候,则他们已经接近正交化时,不再对其正交化处理。
2 词汇文本矩阵奇异值分解的并行化
2.1 划分策略
定义1 A-sweep扫描是指采用Hestense算法对矩阵的所有列两两正交一次。
定义2 AB-supersweep扫描是指采用改进的Hestense算法把一个矩阵(Xi)中的所有列向量与另一个矩阵(Xj)中的所有列向量都分别正交化一次。
在应用中,可用的CPU数量是有限的,而问题的规模是无限的。为了使并行算法具有很好的可扩展性,本文采用了一种任务划分策略,把矩阵的计算分配到多个节点上进行计算。其基本思想是设矩阵为X(m*n),进程的个数为p,则把矩阵按列分为2p个小矩阵,每个小矩阵q=n/2p列(不能整除时补0向量),即 X(m*n)= (X1,X2,…,X2p),其中每个小矩阵Xi (i=1,…,2p)看作一个超级列,然后每个进程上分配两个超级列。
采用这种划分策略后,当采用A-sweep和AB-supersweep两种扫描方法对每个进程上的超级列扫描时,每次扫描同样产生n*(n-1)/2个不同的列向量对。
2.2 并行扫描策略
并行扫描策略定义了矩阵列向量之间的两两正交的顺序,
及各个进程之间如何进行通信。本文采用了文献[1]提出的扫描策略。该策略的优点是每次扫描能够产生n(n-1)/2个不同的Jacobi对,同时把矩阵的n个列向量按摸排列成递增或递减。该策略(如图1所示)包括两个过程:前向扫描和后向扫描。其中前向扫描使得矩阵列向量的模按递增排序,后向扫描使得矩阵列向量的模按递减排序。在计算过程中这两个过程交替地执行。
前向扫描和后向扫描(如图1所示)都是将矩阵的n列分成两行。在同一列的两个列向量组成一个列向量对,每个进程处理一个列向量对,这样n/2个进程可以对n个列向量并行地处理。图1中上下方向的箭头表示同一个进程的两个列交换位置,左右方向的箭头表示进程间的通信方式。
2.3 并行算法的框架描述
Step1 初始化 根据上述任务划分策略把矩阵按列划分成多个超级列,每个进程分别读取两个超级列。
Step2 中止条件 判断是否满足中止条件,若是则转到step5。
Step3 前向扫描 首先对每个进程上的两个超级列分别执行一次A-sweep。每个进程再分别执行q-1次AB-supersweep,每执行一次AB-supersweep各个进程要执行一次同步,并根据前向扫描的方法在各个进程之间通信。
Step4 后向扫描 首先对每个进程上的两个超级列分别执行一次A-sweep扫描,再分别执行q-1次AB-supersweep,每执行一次AB-supersweep各个进程要执行一次同步,并根据后向扫描的方法在各个进程之间通信。转到Step2。
Step5 结束 0号进程收集所有进程上的计算结果。
2.4 并行算法的框架图
3 实验结果及分析
3.1 实验平台及实验数据
本实验在自强3000高性能计算机上测试,集群共有192个节点,394个CPU。实验中用的矩阵是随机生成的n*n的稠密矩阵。并行算法使用c+MPI实现。为了程序的可扩展性和负载平衡,并行程序中使用了对等模式,各个节点平均分配数据。
3.2 实验分析
该算法考虑了潜在语义分析中词汇-文本矩阵SVD计算时要求奇异值排序的特点,在计算SVD的同时得到了从大到小排列的奇异值。从表1和图3得出以下结论:
① 从表1可以看出,并行算法大大缩短了程序运行的时间。
② 当矩阵规模比较小(n<400)的时候,由于通信量占用了较多的时间,所以加速比比较小。随着矩阵规模的逐渐变大,加速比也逐渐变大,当达到一定值(n>400)时趋于稳定。
③ 算法具有很好的扩展性,随着问题规模和CPU数量的增大,加速的变化更加明显。
4 结束语
本文根据潜在语义分析中SVD计算的特点,设计并实现了一种矩阵SVD计算的并行算法,该算法在矩阵正交化的同时能够对矩阵列向量按模排序。同时为了提高算法的可扩展性,本文采用了一种任务划分策略。这样不仅得到了有序排列的奇异值,而且大大缩短了程序运行的时间。实验结果表明了上述算法的有效性。潜在语义分析中词汇-文本矩阵中的词汇和文本经常需要更新,因此下一步工作的目标是针对潜在语义空间更新的具体特点设计一种并行算法。
参考文献
[1]Zhou B B,Brent R P.Aparallel ring ordering algorithmfor efficient one-sided Jacobi SVD computations[C].In Proc.Sixth IASTED/ISMMInternat.Conf.Washington,USA,October,1994:369- 372.
[2]Brent R P,Luk F T.The solution of singular-value and symmetric ei-genvalue problems on multiprocessor arrays[J].SIAMJ.Sci.and Stat.Comput.,1985,6:69 -84.
[3]Demmel J.Jacobi s method is more accurate than QR[J].SIAMJ.Sci.Stat.Comput.,1992,11:1204 -1246.
[4]Becka M,Oksa G,Vajtersic M.Dynamic ordering for a parallel block-Jacobi SVD algorithm[J].Parallel Computing 2002,28:243 -262.
矩阵奇异值 篇3
学生成绩是学籍管理中的重要内容之一, 它是学校评价学生的最主要的指标之一[1~3]。学生成绩恢复是指学生成绩由于漏输或者学生没有选修该课程等原因导致学生成绩数据不完整, 通过一定的技术将该成绩恢复出来。目前, 学校鼓励开设任选课以提高学生的综合素质, 学生具有较大的选择课程的自主权。由于学生选择的课程不同, 任课教师也就不可能一样, 那么教师在给学生评分时尺度也就不一样。同时, 所选修的课程也难度不一样, 所以导致学生最终的分数也存在差别。但学生的得分却是学生评比中最主要的因素。为了克服这些原因导致学生得分不一致, 许多学校在对学生成绩评比中, 不考虑任选课, 而仅考虑全班所有同学都学了的课程。这样导致学生不重视选修课的学生。同时, 对于那些选修课学得好的同学, 也是不公平的。为了克服这些缺点, 本文利用学生成绩表可以表示为一个矩阵, 并利用该矩阵具有低秩的特点, 将没有选修该课程的成绩恢复出来。
基于低秩矩阵缺失元素的恢复在计算机视觉、医学图像分析等领域中具有广泛的用途[4~6]。Tomasi等人利用子矩阵法对矩阵缺失元素进行恢复[7], 但是该方法恢复结果依赖于所选取的子矩阵。为了克服该缺点, Ma等人通过矩阵的行列变换, 将所有的已知元素变换到矩阵的左上角, 再利用这些已知的元素, 一个一个地将缺失元素恢复[8]。但是该方法在恢复缺失元素时, 利用了已恢复元素的信息, 这样必然会导致误差的积累。有些学者采用进化的方法对缺失元素进行恢复[9], 但该方法运算量大, 而且易出现早熟现象。有些学者采用非线性规划的方法进行求解[10], 但是这些方法收敛速度比较慢, 尤其是到了算法的后期。
为了克服上述缺点, 本文采用迭代奇异值分解的方法对缺失元素进行恢复, 该方法对矩阵进行SVD分解, 得到一个低秩矩阵, 将低秩矩阵中的元素代替未知元素, 再循环进行SVD分解。
该方法的优点是在矩阵恢复过程中, 利用了所有已知元素信息, 而且将所有已知元素平等地对待, 这样必然恢复精度高。
1 迭代矩阵奇异值分解的缺失数据恢复方法
近年来, 基于低秩矩阵的恢复受到了广泛地关注[11~12], 它是一个求矩阵最小秩的优化问题, 其求解非常困难, 但可通过求其最紧凸闭包的矩阵核范数问题, 它的解和原始问题的解非常接近[13], 即:
式中:Ω 表示已知元素下标的矩阵集合, A矩阵为已知数据矩阵, ||∙||*表示核范数, 即矩阵的奇异值之和。
在学生成绩的获取过程中, 由于各种原因, 可能会出现一些数据的缺失, 即矩阵A含有缺失元素。矩阵缺失元素恢复问题可表示为如下形式[14]:
式中:M是为非负的缺失矩阵, 它与矩阵A大小相同, 当矩阵A中元素的值已知时, 对应的矩阵M处的值为1, 否则为0; (·) M表示在矩阵M上进行投影。
可通过拉格朗日求解式 (2) , 即:
式中Z为引入的矩阵。
由于教师在批改试卷时, 会受到许多因素的影响, 必然会出现判定给分时尺度不一, 因此, 学生成绩数据可以认为含有噪声。当矩阵中含有噪声时, 式 (2) 可重写为:
式中 τ 和 θ 均为参数。
利用奇异值阈值方法, 可采用式 (5) 迭代形式求解式 (3) :
式中:Hτ (⋅) 表示对矩阵进行SVD分解;表示矩阵M的补集;τ 为一个非负阈值, 用于控制矩阵秩的取值。
本文方法是一个线性迭代过程, 首先利用初值对矩阵进行SVD分解, 提取矩阵中部分主成分, 得到一个与原矩阵近似的低秩矩阵, 降低了噪声的影响;提取恢复出的缺失元素, 并结合已知元素再进行SVD分解, 直到收敛。
1.1 计算过程
本文方法是一个线性迭代过程, 首先利用初值对矩阵进行SVD分解, 提取矩阵中部分主成分, 得到一个与原矩阵近似的低秩矩阵, 降低了噪声的影响;提取恢复出的缺失元素, 并结合已知元素再进行SVD分解, 直到收敛。
1.1.1 计算矩阵Z
本文方法是一个线性迭代过程, 首先利用初值对矩阵进行SVD分解, 提取矩阵中部分主成分, 得到一个与原矩阵近似的低秩矩阵, 降低了噪声的影响;提取恢复出的缺失元素, 并结合已知元素再进行SVD分解, 直到收敛。
首先, 对矩阵Y进行SVD分解, 即:
利用阈值 τ , 求到的主成分的个数rN, 得到矩阵Z为:
1.1.2 计算矩阵Y
根据矩阵Z值, 对矩阵Y进行更新, 即:
同时, 为了有效性地度量恢复矩阵与原始矩阵之间的差异, 本文定义恢复误差为:
可采用下式作为迭代终止条件, 即:
式中 ε 为一个给定的非负阈值。
1.2 算法总结
可采用下式作为迭代终止条件, 即:
(1) 初始化Z0、ε , 给矩阵Z0中的缺失数据赋初始值;
(2) 利用式 (6) 和式 (7) 计算矩阵Z值;
(3) 利用式 (9) 计算恢复误差e , 若满足式 (10) 的终止条件则转至第 (6) 步;否则转第 (4) 步;
(4) 按照式 (8) 计算矩阵Y的值;
(5) 转第 (2) 步;
(6) 输出Yk=Yopt和恢复误差值。
2模拟实验与真实实验
图1 为在不同噪声下恢复误差随迭代次数的变化, 图2 为在不同的缺失率下恢复误差随迭代次数变化。为了检验本文方法的收敛性能, 在Matlab环境下, 通过计算机模拟产生一个大小为80×100 秩为2 的矩阵, 其中每个元素在 (1, 100) 范围内, 随机缺失30%的元素, 同时在矩阵中加入均值为0, 方差分别为0, 1, 2 的高斯噪声。利用这些模拟产生的数据进行缺失元素恢复, 并利用式 (10) 求取恢复误差, 模拟结果如图1 所示。
为了研究本节方法的收敛性与遮挡率之间的关系, 在上述实验中加入均值为0, 方差为2 的高斯噪, 并分别随机缺失10%, 20%及30%的矩阵元素, 实验结果如图2所示。
从图1 和图2 可以看出, 本节方法最多迭代10 次便能够达到收敛, 因此它具有良好的收敛性。从图2 中还可以看出, 矩阵元素缺失率越低, 算法收敛速度越快。原因是由于缺失率越低, 提供的约束信息就越多, 因此算法收敛速度越快。
为了比较本文方法和Tomasi的子矩阵法[7]的收敛精度, 在矩阵元素缺失率为20%的情况下, 并在矩阵中加入高斯噪声, 分别用本文方法和Tomasi的子矩阵法进行缺失元素恢复, 在每种噪声水平下各进行100 次实验, 并分别求取平均恢复误差, 实验结果如图3 所示。
从图3 可以看出, 本节方法的误差比较小。原因是由于本文方法在恢复矩阵元素时, 利用了所有的矩阵元素信息, 且将所有的元素都平等地对待, 并没有倚重某些元素, 因此算法的误差比较小。而子矩阵法在恢复缺失元素时, 仅利用了子矩阵中的元素, 并没有利用全部已知元素的信息, 其恢复精度依赖于所选取的子矩阵, 因此, 子矩阵法的恢复精度比较低。
为了验证本方法的正确性, 本文利用陕西师范大学计算机科学与技术专业某毕业班的所有课程成绩组成一个矩阵, 该班级共有53 名学生, 共有115 门课程, 由于有些课程只有部分学生选修, 因此, 矩阵中含有大量的缺失数据, 缺失率为44%, 其中部分成绩如表1 中所示, 表中“*”表示该学生没有选修该课程。
采用本文方法对学生成绩进行恢复, 恢复结果如表2 所示。
为了检验本文方法的正确性, 本文首先采用我院现的方法计算学生平均成绩及排名, 即仅考虑学生的必选课的平均成绩进行排名, 因为必选课所有的学生都必须选修, 这样必选课所有的学生都有成绩。而所有的任选课都不考虑, 因为任选课有些学生没有选修该课程, 即没有该门课程的成绩, 学生平均成绩及排名如表3 所示。
最后利用本文方法对学生的成绩进行恢复, 利用恢复后的成绩求平均成绩, 并根据平均成绩进行排名, 恢复后平均成绩及排名如表3 所示。
从表3 可看出, 经过恢复后的平均成绩和原始平均成绩非常接近。从排名上看, 两者排名会稍微有点变化, 但变化不会太大, 这是由于本文方法在对缺失元素进行恢复时, 还是利用了基于该生的各科课程成绩及该生在该课程在所有选修学生中的排名。但是本文的排名更加合理, 原因是由于本文的方法考虑了学生的任选课成绩, 而现有的方法并没有考虑任选课的成绩。
3 结语
本文为了有效地对学生成绩数据进行恢复, 提出了一种迭代奇异值分解的学生成绩恢复方法。该方法将学生成绩表示为一个矩阵, 利用该矩阵具有低秩的特性, 在给定缺失元素的初始值后, 采用迭代奇异值分解的方法, 最终求取成绩表中缺失元素的真实值。模拟实验和真实实验结果表明该方法能够快速、正确地恢复出学生的真实成绩。
摘要：为了有效地对学生成绩数据进行恢复, 提出一种迭代奇异值分解的学生成绩恢复方法。该方法采用矩阵表示学生成绩, 利用该矩阵具有低秩的特性, 在给定缺失元素的初始值后, 利用奇异值分解得到缺失元素的近似值, 而该近似值比初始值更加接近真实值。再将求到的近似值代替初始值, 经过多次迭代, 最终可求到成绩表中缺失元素的真实值。该方法的优点是在缺失元素恢复过程中, 利用了所有已知元素信息, 并将所有已知元素平等地对待。模拟实验和真实实验结果 表明该方法能够快速、精确地恢复出学生的真实成绩。
奇异值数字水印算法的数学基础 篇4
数字图像的奇异值分解。数值分析中的奇异值分解 (SVD) 是一种将矩阵对角化的数值算法。在图像处理中应用的主要理论背景是: (1) 图像奇异值的稳定性非常好, 当图像被施加小的扰动时, 图像的奇异值不会有大的变化; (2) 奇异值所表现的是图像的代数特性。
从线性代数的角度看, 一幅数字图像可以看成是一个非负矩阵。用A∈Rm×n来表示这样一个图像矩阵, 其中R表示实数域。矩阵A的奇异值分解如下:A=UDVT
其中都是正交矩阵, 是A或AT的奇异值, 其中R是A的秩, 是AAT或ATA特征值的平方根。U、V分别是A的左奇异向量和右奇异向量。因为
所以U的列向量是AAT的特征向量, V的列向量是ATA的列向量, 并且它们所对应的特征值都是A的奇异值的平方。
基于HVS的图像块分类。根据HVS的特性, 人眼对不同类型图像区域的敏感程度依次为:对图像平滑区域的噪声最敏感;对图像的边缘信息较为敏感;对纹理区域的噪声不敏感。根据这个特性, 我们希望在图像的纹理区域进行水印的嵌入。文献基于图像块的灰度均值、熵和方差等统计特性将图像块进行划分。
图像处理与信息论的相关理论表明, 方差反应各像素灰度值与图像平均灰度值的总的离散程度, 而熵则表征图像所含信息的多少, 图像纹理区域熵值较大而方差较小。在这里, 我们引入, 其中H表示熵, 为标准差, 显然熵值越大, 方差越小, T值越大。通过引入T值, 我们在对图像的各分块进行比较时, 就不必对方差和熵值分别进行比较, 而统一对T值大小比较即可, 这样对于各分块的比较变得简单明了, 同时也利于程序的编写。
水印的嵌入与提取。算法步骤如下:水印嵌入:首先将图像进行分块, 计算每块图像的T值, 选取T值最大的块, 进行奇异值分解, 将水印嵌入到进行奇异值分解的图像块中。方法如下:水印是由服从高斯分布的实数序列组成, 其中Wi相互独立, n为水印的长度。挑选出除了d1以外的n个数值最大的奇异值组成的单调递减序列, 之所以排除d1是因为d1远远大于其它的奇异值。对d1的改变会导致图像质量的下降。根据将水印w嵌入到S中, 得到一个嵌入后的奇异值序列。其中
a是一个比例参数, 用来调节w改变S的程度, 从而控制嵌入水印后的图像质量。将S'回代入D中, 生成D', 然后再UD'VT就得到嵌入水印后的图像块A'。由于调整后的奇异值序列可能不符合单调递减的顺序, 而在奇异值分解后奇异值总是递减排列的, 所以根据d'i的递减顺序对S中对应的di排序, 生成序列。用它来标明A的奇异值及其原始顺序。用它来标明A的奇异值及其原始顺序。
提取水印时需要:嵌入水印的图像 (可能已经遭受到攻击) 、参考水印w、奇异值序列和比例参数a。
将进行分块, 确定原先嵌入水印的图像块, 设为A*, 并从中提取出可能已遭到破坏的水印w*并进行检测:
(1) A*进行奇异值分解, 得到按递减顺序排列的奇异值向量D*。
(2) 将按单调递减顺序排列, 还原成。
(3) 根据S和, 使D*中的奇异值与S的相应元素对齐, 并从中D*提取出奇异值序列。
(4) 提取水印。
(5) 计算标准相关系数:
利用corr2 (w, w*) 来测量w与w*相似程度。
实验结果:我们使用256级灰度图像Lena.bmp为例, 来测试本水印算法对各种几何攻击的稳健程度。Lena.bmp的大小为256×256, 先将图像进行64×64分块, 再嵌入水印。本文只给出算法对放大、旋转等几何攻击以及对剪切攻击的结果。
分块的大小为64×64, 取水印长度N=50。实验结果如下:原图像 (图1) 。
未攻击时, 在原图像的下半部分加入水印, a=0.3, 图像为 (图2) 。
当a=0.3时, 峰值信噪比psnr=46.9693;
从上面各幅图像中提取的水印并给出相关系数r=0.9991
下面我们对嵌入水印的图像施加各种攻击
旋转90度, a=0.3, 相关系数r=0.9991
横向拉伸2倍, a=0.3, 相关系数r=0.9991
纵向拉伸2倍, a=0.3, 相关系数r=0.9991
可见本文的算法对常见的几何攻击是非常稳健的。经过几何攻击后的图像可以非常稳健的提取水印。
现在我们施加其他攻击:
裁减左半幅图像, a=0.3, 从中提取的水印与原始水印的相关系数r=0.9991
裁减下半幅图像, a=0.3, 从中提取的水印与原始水印的相关系数r=0
实验结果表明, 本文算法对剪切有一定的抵抗性, 关键是剪切的位置, 若把嵌入水印的图像块剪切掉, 那不能提取出水印;若没有剪切掉嵌入水印的块, 则可以非常稳定的提取出水印。
本文提出了一种能够抵抗一般几何攻击的数字水印算法, 由服从高斯分布的随机数组成的水印被嵌入到数字图像的奇异值之上, 并考虑到人类视觉系统特性的影响, 引入了基于统计量熵值和方差的参数T, 作为我们划分不同种类块的标准。实验结果表明本水印算法对这些几何攻击具有很好的稳定性。
摘要：本文提出了一种利用奇异值分解的图像数字水印算法, 针对的是在水印提取之前遭受了未知几何攻击的图像。实验结果表明, 本文的算法对放大、旋转等几何攻击以及对剪切, 噪声攻击是稳健的, 而且此水印在受到了上述的几何攻击后仍能够被有效的提取。
关键词：奇异值分解,数字水印,几何攻击
参考文献
[1]刘瑞祯, 谭铁牛.基于奇异值分解的数字图像水印方法[J].电子学报, 2001, 29 (2) :pp.168~171.
[2][美]G..H.戈卢布C.F范洛恩[J]《.矩阵计算》, 8.1&8.6.
矩阵奇异值 篇5
滚动轴承是机械设备中最常用且最容易损坏的零件之一,当其发生故障时损失巨大,因此轴承故障诊断技术得到了人们的重视。轴承的内圈、外圈以及滚动体等出现损伤,会在轴承特征频率上产生一系列冲击振动,使得原来的平稳振动信号变成非平稳振动信号。小波变换由于具有“数学显微镜”和多分辨率的特性,在机械故障诊断领域得到了广泛应用[1,2,3]。但是,小波基函数的选择问题是小波分析的一个难题。
Huang等[4]提出的希尔伯特-黄变换(HHT)是一种新的具有自适应的时频分析方法,已经成为机械故障诊断领域研究的热点[5,6,7,8]。目前,基于HHT变换的故障诊断方法大部分是先利用经验模态分解(EMD)方法分解得到内蕴模式函数,再对内蕴模式函数进行分析提取机械故障特征。张超等[5]利用齿轮振动信号EMD分解后得到的内蕴模式函数的能量熵作为特征,对齿轮进行裂纹和断齿故障诊断。高强等[6]利用EMD对轴承振动信号进行分解,将局部损伤轴承产生的高频调幅信号成分作为内蕴模式函数分离出来,然后用Hilbert变换得到其包络信号,通过包络谱提取轴承故障特征频率。程军圣等[7]提出了一种基于内蕴模式函数奇异值分解和支持向量机的故障诊断方法,该方法利用内蕴模式函数形成向量矩阵,然后对该矩阵进行奇异值分解,提取其奇异值作为故障特征向量。相比较而言,利用Hilbert谱进行故障诊断的文献很少,于德介等[8]提出了一种基于Hilbert谱时频熵的齿轮故障诊断方法,该方法将Hilbert谱的时-频平面等分为N个面积相等的时频块,对每块进行能量归一化,然后仿照信息熵的方式定义时频熵,实验表明,不同齿轮状况的时频熵不同。由于Hilbert谱中包含了机械故障中蕴含的非常丰富的信息,因此针对Hilbert谱进行特征提取以便对机械故障进行诊断具有重要的研究价值。
奇异值分解(singular value decomposition,SVD)是一种重要的矩阵分解,奇异值个数反映了原矩阵中独立行(列)矢量的个数。目前,奇异值分解技术已获得了广泛应用,如应用于数据压缩[9]、信号降噪[10]、机器状态监测[11]等。
本文提出一种基于Hilbert谱奇异值分解的特征提取方法,首先利用EMD将轴承振动信号分解为若干个内蕴模式函数(IMF)分量之和,然后对每个IMF分量进行Hilbert变换得到瞬时频率和瞬时幅值,从而得到轴承振动信号的Hilbert谱,Hilbert谱表示了信号完整的时间-频率分布,接着对Hilbert谱进行奇异值分解,得到的奇异值作为轴承故障诊断的特征向量,最后利用支持向量机(support vector machine,SVM)进行故障分类,验证上述特征提取方法的有效性。
1 Hilbert谱
Huang等[4]提出的EMD方法目的是将待分析信号分解为一系列表征时间尺度的IMF分量,使得各IMF分量是窄带信号,满足Hilbert变换的要求,即IMF分量必须满足两个条件:采样数据的极大点和极小点数之和与过零点的个数之差不超过1;在任意点,由局部极大值和局部极小值定义的包络均值必须为零。
EMD的分解过程其实是一个“筛分”过程,在“筛分”的过程中,不仅消除了模态波形的叠加,而且使波形轮廓更加对称。EMD方法通过下面的步骤对信号x(t)进行分解[4,12]:
(1)确定信号所有的局部极值点,然后用三次样条将所有的局部极大值点连接起来形成上包络线。
(2)再用三次样条线将所有的局部极小值点连接起来形成下包络线,上下包络线应该包络所有的数据点。
(3)上下包络线的平均值记为m1,求出:
理想地,如果h1是一个IMF,那么h1就是x(t)的第1个IMF分量。
(4)如果h1不满足IMF的条件,把h1作为原始数据,重复步骤(1)~步骤(3),得到上下包络线的平均值m11,再判断h11=h1-m11是否满足IMF的条件,如不满足,则继续循环,直到得到的h1k满足IMF的条件为止。记c1=h1k,则c1为信号x(t)的第1个满足IMF条件的分量。
(5)将c1从x(t)中分离出来,得到
将r1作为原始数据重复步骤(1)~步骤(4),得到x(t)的第2个满足IMF条件的分量c2,重复循环n次,得到信号x(t)的n个满足IMF条件的分量ci(i=1,2,…,n)。记
当rn成为一个单调函数不能再从中提取满足IMF条件的分量时,循环结束。
式(3)中,rn称为残余函数,代表信号的平均趋势。
将EMD分解得到的每个内蕴模式函数ci(t)作Hilbert变换,得
构造解析信号:
于是得到幅值函数:
相位函数:
可以求出瞬时频率:
这样,可以得到:
其中,RP表示取实部。展开式(9)即称为Hilbert谱[8],记作H(ω,t)。
H(ω,t)精确地描述了信号的幅值在整个频率段上随时间和频率的变化规律。
2 奇异值分解
由Hilbert谱得到的振动信号的时频分布矩阵记为A,假设A为m×n阶矩阵,则由奇异值分解理论可知,存在正交矩阵
使得UTAV=diag[σ1σ2…σp]=S,p=min(m,n)。即
则该式称为A的奇异值分解。其中σ1≥σ2≥…≥σp≥0,σi(i=1,2,…,p)为A的奇异值,U和V中分别是A的奇异向量。
矩阵A的奇异值分解提供了一些关于矩阵A的重要信息,利用奇异值和奇异向量可以刻画矩阵的本身结构,其中大的奇异值包含了矩阵中所蕴含的模式的更多信息。
对机械故障振动信号进行Hilbert-Huang变换得到的Hilbert谱时频矩阵,维数比较高,数据量大,且包含大量的相互关联的特征,不利于下一步的分析计算。可以利用奇异值分解的方法提取故障信号Hilbert谱矩阵的特征,通过对Hilbert谱矩阵进行奇异值分解而得到的奇异值代表了振动信号Hilbert谱的时频结构特征。机械在出现故障的情况下,时频结构会出现变化,因此,Hilbert谱的奇异值可以用来作为机械故障诊断的特征。
3 轴承信号的Hilbert谱分析
试验采用6205-2RS JEM SKF型轴承,轴承的内直径为25.0012mm,外直径为51.9989mm,厚度为0.5906mm,节径为39.0398mm,滚动体直径为15.0012mm,其中内圈故障、滚动体故障、外圈故障直径为0.1778mm,电动机负载0HP,电机转速为1797r/min。试验装置由一个1.5kW电动机、一个扭矩传感器/译码器、一个功率测试计与相应电器控制装置组成。将振动加速度传感器安装在带有磁力基座的机架上,振动信号由16通道数据记录仪采集得到,采样频率为12kHz。选择其中时间长度为0.25s的故障数据,时域信号如图1所示。
对轴承振动信号进行EMD分解,得到内蕴模式函数IMF,然后对得到的IMF进行Hilbert变换,得到Hilbert谱,其中轴承正常、内圈故障、滚动体故障、外圈故障4种工况信号频率范围为0~1200Hz的Hilbert谱如图2所示。从图中可以看出,轴承振动信号的Hilbert谱体现出明显的非平稳特征,其中包含了丰富的信息,能量主要集中在低频0~300Hz范围。从图2还可以看到4种工况振动信号的Hilbert谱有不同之处,但并不容易对4种工况进行区分。
Hilbert谱是轴承振动信号的一种时频域表示方法,为了提取轴承振动信号在时频域的特征,对得到的Hilbert谱使用奇异值分解技术。其中较大的奇异值含有振动信号结构的更多信息。轴承正常状态、内圈故障、滚动体故障、外圈故障振动信号的Hilbert谱的前50个奇异值均值如图3所示。从图3中可以看到,轴承不同工况振动信号的Hilbert谱的奇异值明显不同,轴承在正常情况下的奇异值较小,当故障出现后,其Hilbert谱的奇异值增大,并且,不同故障的奇异值不同,外圈故障时的奇异值最大。从图中可以看出,轴承振动信号Hilbert谱的奇异值可以作为轴承故障诊断的特征。
4 应用
故障诊断试验所用的实测轴承振动加速度数据来自于Case Western Reserve University(CWRU)[13]。利用其中的“Drive End Bearing Fault Data”,故障类型有外圈故障、滚动体故障和内圈故障,采样频率为12kHz,对载荷为0的不同损伤程度的10组数据进行9种不同的故障诊断测试,测试数据集的选择参考文献[14]的方法,其中包含了不同故障状态的组合,数据集的详细描述如表1所示。数据集中每个样本的采样点数为4096。其中轴承正常情况的样本数为59个,随机选择其中29个样本作为训练样本集,其余30个样本作为测试样本集;其余每种故障情况的样本数为29个,随机选择其中14个样本作为训练样本集,其余15个样本集作为测试样本集。故障数据集名称中的3组数字分别代表轴承的3种故障状态所对应的故障直径的大小,如对于D070707数据集,代表内圈、滚动体和外圈的故障直径都为177.8μm(7mils),对于DBALL、DINN和DOUT数据集分别包含了滚动体故障、内圈故障和外圈故障3种故障程度的数据。
支持向量机是在统计学习理论基础上发展起来的一种新型机器学习方法,具有泛化能力强、维数不敏感等优点,适于求解高维、小样本、非线性情况下的模式分类和回归分析等问题。因此,这里选择SVM作为轴承故障诊断的分类器。
试验中选用“一对一”的SVM多类识别方法,核函数选择最常使用的高斯核函数形式,核函数形式为
其中,σ′为控制核函数高宽的参数。对于线性不可分情况,引入惩罚因子C来控制错误分类。参数取值为σ′=1,惩罚因子C=100。
应用Hilbert谱前20个奇异值作为特征对表1的数据集进行训练与识别,故障分类结果如表2所示。从表2中可以看出,在9个故障诊断测试集上都取得了比较高的识别率,说明了基于Hilbert谱奇异值方法用于轴承故障诊断的有效性。
5 结论
本文提出一种基于Hilbert谱奇异值的机械故障特征提取方法。轴承信号的Hilbert谱包含了时间、频率两方面的信息,是分析非平稳信号的有力工具,通过结合奇异值分解方法,较好地实现了Hilbert谱时频特征的有效提取。通过轴承正常、内圈故障、滚动体故障、外圈故障4种工况实测信号的识别试验,验证了本文提出的方法可以准确提取轴承故障特征。本文提供了一种机械故障诊断的新思路。
摘要：针对机械故障振动信号时频特征提取问题,提出一种基于Hilbert谱奇异值的特征提取方法,并将其应用于轴承故障诊断。该方法首先利用经验模式分解方法将振动信号分解为若干个内蕴模式函数之和,接着对每个内蕴模式函数进行Hilbert变换得到振动信号的Hilbert谱,然后对Hilbert谱进行奇异值分解,得到反映机械状态特征的奇异值序列,最后利用奇异值作为特征向量,使用支持向量机进行轴承故障诊断。轴承正常、内圈故障、滚动体故障、外圈故障实测信号实验结果表明,该方法能有效地提取轴承故障振动信号特征。
矩阵奇异值 篇6
盲源分离(blind source separation,BSS)［１－６］作为一种预处理手段,可以从机械设备测量信号中分离(或恢复)出各个机械部件的振动源信号。目前,盲源分离方法,如独立成分分析(independ-ent component analysis,ICA)法,已在机械设备故障诊断领域取得了初步的应用。
在机械设备盲信号处理中,目前的方法大多是基于分离源信号来提取故障特征信息并进行故障诊断［３－６］,而很少从分离矩阵出发来处理。但是从信息的角度来看,分离源信号和分离矩阵是从两个不同的角度和层面来阐释机械设备故障特征信息的,它们各有侧重。在一定意义上,分离矩阵包含的信息与源信号的信息是互补的,它们共同表征了机械设备的状态。若将分离源信号和分离矩阵包含的信息互相融合,则必然能够更全面、更深刻地表征其状态。因此,本文将盲信号处理、包络分析、奇异值分解和信息熵等信号处理方法相结合,提出基于奇异值融合的故障特征信息提取方法,从分离源信号和分离矩阵融合的角度来提取机械设备的故障特征信息,并将其应用于液压齿轮泵的故障特征信息提取中。
1 基本概念和理论
1.1 包络分析
包络分析(envelope analysis,EA)［７－１１］,又称解调分析,目的是提取载附在高频调制信号上的与故障有关的低频信号,从时域上看,就是提取时域信号波形的包络轨迹。包络分析是目前最常用、最有效的机械设备故障诊断方法之一。调制信号的包络解调方法主要有三种:Hilbert幅值解调法、检波- 滤波解调法和高通绝对值解调法。Hilbert幅值解调法由于简单、有效,故在齿轮和滚动轴承等机械信号处理中已应用较多［１２］,本文也利用Hilbert幅值解调法来提取源信号的上下包络信号。
1.2 奇异值分解
设有M行N列的实矩阵A,它可以作奇异值分解(singular value decomposition,SVD)［１３－１４］,若有
则U和V分别称为实矩阵A的左右奇异阵。Λ 为矩阵[diag(σ１,σ２,…,σｐ):0] 或其转置,Λ ∈RＭ×Ｎ,p= min(M,N),σ１≥σ２≥…≥σｐ≥0,σ１,σ２,…,σｐ称为实矩阵A的奇异值。
可进行奇异值分解是矩阵的固有特征,矩阵奇异值具有比较高的稳定性,即当矩阵元素发生小的变动时,矩阵的奇异值变化很小;同时矩阵的奇异值还具有比例不变性、旋转不变性和可降维压缩等特性。因此,矩阵的奇异值符合机械故障诊断学对于故障特征信息的基本要求。
1.3 奇异值熵
为了定量地描述矩阵奇异值的变化程度,引入信息熵理论,根据信息熵的定义来提出并构造奇异值熵。
假设奇异值统一编号为σ１、σ２、…、σｎ,对其进行归一化处理,即可得到:
其中,E =σ１+σ２+ … +σｎ,从而有
这符合计算信息熵的初始归一化条件,根据信息熵理论,构造奇异值熵为
2 基于奇异值融合的盲信息提取方法
2.1 盲信息提取策略
2.1.1 源信号包络矩阵奇异值
在机械故障诊断中,齿轮和滚动轴承等部件的振动信号是调制类型信号,而且传感器测量的信号是这些振动信号的混合信号,这就大大增加了故障特征信息提取的难度,使得故障特征信息的准确度降低。盲源分离技术可以从已知的混合信号(测量信号)中分离(或恢复)出各个振动信号(源信号),这些源信号提供了较“纯”的信息载体,可以从中提取更为准确的故障特征信息。进一步,振动源信号是调制信号,其上下包络更集中地携带了机械设备的运行状态信息,因此,分离源信号进行包络解调可以得到上下包络信号,进而,上下包络信号可以分别组成上下包络信号矩阵,称为初始特征向量矩阵,它们可以准确地刻画振动源信号的特征信息。若将初始特征向量矩阵进行奇异值分解,则可提取其奇异值作为机械设备的初始故障特征信息。
一般,盲源分离得到的源信号有多个,因此其上下包络信号矩阵的维数也比较大,这使得包络矩阵的奇异值比较多,导致后续模式识别的复杂性和计算量也大幅增加。为了减少奇异值数量和计算量,同时又不损失或少损失故障特征信息,考虑对奇异值进行降维压缩,将上下包络矩阵奇异值的均值和奇异值熵,作为降维的故障特征信息。
2.1.2 盲源分离矩阵奇异值
在盲源分离中,分离(或恢复)源信号相当于估计分离矩阵;当分离矩阵确定了,源信号也就很容易得到。若从信息的角度来看,机械设备的部分故障特征信息不仅包含在分离源信号中,也包含在分离矩阵中,源信号和分离矩阵是从两个不同的角度和层面来阐释机械设备故障特征信息的。因此,可将分离矩阵直接作为初始特征向量矩阵进行奇异值分解,并从中提取机械设备的故障特征信息。
2.1.3 基于奇异值融合的特征信息提取
在一定意义上,源信号是从宏观现象来表征机械设备的故障特征信息的,而分离矩阵则是从微观结构来刻画机械设备的故障特征信息的,它们是从两个不同的角度和层面来阐释同一机械设备的故障特征信息的。若从信息融合的角度来考虑,可将源信号包络矩阵奇异值(均值和奇异值熵)和分离矩阵奇异值融合共同作为机械设备的故障特征信息,这必然能够更全面和深刻地表征机械设备的运行状态信息。
2.2 盲信息提取流程和步骤
综上所述,本文提出的基于奇异值融合的机械盲信息提取方法,即将源信号包络矩阵奇异值均值、奇异值熵和分离矩阵奇异值进行特征层信息融合,共同作为机械设备的故障特征信息。它的流程如图1所示。具体步骤为:
(1)机械信号测量。从待检测机械设备上测得多通道传感器观测信号(混合信号)。
(2)盲源分离。混合信号经中心化和白化等预处理,利用盲源分离技术(例如FastICA算法等)计算分离矩阵,并分离调制源信号。
图1基于奇异值融合的盲信息提取方法
(3)源信号包络矩阵奇异值。利用Hilbert幅值解调法提取调制源信号的上下包络信号,进而组成初始特征向量矩阵,进行奇异值分解,计算奇异值均值和奇异值熵。
(4)分离矩阵的奇异值分解。分离矩阵进行奇异值分解,得到它的奇异值。
(5)奇异值融合。源信号包络矩阵奇异值(均值和奇异值熵)和分离矩阵奇异值进行融合,共同作为机械设备的故障特征信息。
3 液压齿轮泵试验
为了验证该方法的有效性,将其应用于国产CB-Kp63型液压齿轮泵中。齿轮泵试验台架和加速度计(传感器)的设置如图2所示,其中,泵轴转速为定速1480r/min。
图2液压齿轮泵试验台架及加速度计设置
在密闭实验室环境下,设定齿轮泵故障模式类型包括:正常状态、齿面磨损和轴承故障;在每个故障模式下各测量64 组传感器观测信号x(t)=(x１(t),x２(t),x３(t),x４(t))Ｔ,即每个故障模式均为64组4通道混合信号,其中,泵壳振动信号采样频率设为10kHz,采样时间设为1s。图3为其中一组混合信号的时域波形。
图３液压齿轮泵的混合信号时域波形
混合信号x(t)经中心化和白化等预处理,再利用FastICA算法进行盲源分离,可以得到分离矩阵W和源信号s (t),s (t) =(s１(t),s２(t),s３(t),s４(t))Ｔ。图4a、图5a、图6a即为图3混合信号x(t)对应源信号s(t)的时域波形。进一步,源信号进行包络解调可以提取其上下包络信号(图4~图6)。
在得到源信号上下包络矩阵,即初始特征向量矩阵后,再进行奇异值分解,可以得到对应的奇异值向量σ（１）=(σ１（１）,σ２（１）,σ３（１）,σ４（１）)和σ（２）=(σ１（２）,σ２（２）,σ３（２）,σ４（２）),其中,上标“1”和“2”分别代表上下包络矩阵。进一步,由σ（１）和σ（２）计算其奇异值均值和奇异值熵。
同时,分离矩阵W进行奇异值分解也可以提取故障特征信息。例如,由图3混合信号x(t)可以得到它的分离矩阵:Wnormal(正常状态)、Wgear(齿面磨损)和Wbearing(轴承故障),它们再进行奇异值分解可以得到对应的奇异值向量:
分析和比较上述奇异值数值,可以得出如下结论:
(1)与源信号上下包络矩阵的奇异值向量(初始特征向量)相比,奇异值均值和奇异值熵的聚类划分特性更为典型,其数值稳定性也更为鲁棒。
(2)源信号包络矩阵的奇异值熵可以定量地表征齿轮泵的运行状态信息,在同一泵轴转速下,齿轮泵奇异值熵值从小到大的排列顺序依次为:正常状态、轴承故障、齿面磨损。
(3)盲源分离矩阵的奇异值同样也具有类间差异显著而类内聚类集中的良好特性。
(4)与单独的源信号包络矩阵奇异值或分离矩阵奇异值相比,融合的奇异值特征信息向量可以更全面、更准确地表征齿轮泵的故障特征信息,因而具有更为优良的故障特征信息刻画性能。
奇异值特征信息向量的分布如图7所示。由于奇异值特征信息向量是一个9维向量,不易直观观察,因此,这里仅给出了它的2D平面分布(第5维和第6维)和3D空间分布(第1维、第5维和第6维)。
由图7可知,基于奇异值融合的机械盲信息提取方法是可行的,也是有效的。
图4液压齿轮泵的源信号和上下包络信号时域波形(正常状态)
图5液压齿轮泵的源信号和上下包络信号时域波形(齿面磨损)
图6液压齿轮泵的源信号和上下包络信号时域波形(轴承故障)
4 结语
在机械盲信号处理中,分离源信号和分离矩阵是从两个不同的角度和层面来阐释机械设备故障特征信息的,它们各有侧重,互为补充。本文将分离源信号和分离矩阵包含的故障特征信息相融合,提出基于奇异值融合的机械盲信息提取方法,并利用液压齿轮泵试验验证了该方法的可行性和有效性。
矩阵奇异值 篇7
表面形貌是指零件在加工过程中残留下来的微观几何形态,可用粗糙度、波纹度、形状误差等来表征。表面形貌客观地反映了表面生成机制,影响工程表面的功能特性[1]。为了有效地分离粗糙度、波纹度、形状误差等表面形貌特征,国内外许多学者对表面滤波技术进行了深入研究[2,3]。目前,表面滤波的方法主要有多项式拟合法、2RC滤波法、高斯滤波法、稳健高斯回归滤波法、样条滤波法和小波滤波法等。多项式拟合法能保证数学处理速度,但受函数形式和多项式次数的制约,其拟合精度有限;2RC滤波器能有效分离出高频和低频部分,但存在非线性相移的缺陷;高斯滤波法能有效分离表面成分且无相移,但存在边界效应和稳健性差的问题;稳健高斯回归滤波法解决了高斯滤波法存在的边界效应和稳健性差的问题[4];样条滤波法能有效提取大曲率曲线且没有边界效应,但其频率传输特性较差[5,6,7]。上述滤波法都缺乏多尺度性能,小波滤波法则具有多尺度性能但存在最优小波函数选择的问题[8]。针对目前表面滤波法所存在的问题,本文提出了一种基于多分辨率奇异值分解(MSVD)的表面滤波方法。该方法不存在边界效应和最优小波函数选择的问题,并且具有多尺度性。
1 一维多分辨率奇异值分解原理
对于任何一维信号X=(x1,x2,…,xN),构造矩阵
对矩阵H进行奇异值分解[9](singular value decomposition,SVD)得到:
其中,δ1≥δ2,u1、u2、v1、v2分别为SVD分解后得到的列向量,H=δ1u1v1T+δ2u2v2T。记H1=δ1u1v1T和H2=δ2u2v2T,H1对应的是大奇异值,反映的是H的主体成分;H2对应的是小奇异值,反映的是H的细节成分。令Hj为j尺度下的构造矩阵,Hj,1为Hj的主体成分,Hj,2为Hj的细节成分;Aj为j尺度下信号的主体成分,Dj为j尺度下信号的细节成分。由式(2)可知:
一维MSVD递推分解方式如图1所示[10]。
其中,Aj→Aj+1、Dj+1算法如图2所示。
由式(2)~式(4)可知:Aj=Aj+1+Dj+1,Aj和Aj+1都为N维向量,所以一维MSVD不会存在边界效应。对一维离散信号进行MSVD可知,MSVD与小波分析具有一定的差别,小波分析的分解层数是有限的,而MSVD不受分解层数的限制。文献[11]证明了{Dj}构成近似公比为0.5的等比数列,所以{Aj}是收敛的,则信号通过SVD迭代分解最终得到的结果是稳定的。
2 二维多分辨率奇异值分解的构造
根据一维MSVD分析理论,本文结合二维小波分析中Mallat算法来构造二维MSVD,其算法如图3所示。其中,Ajx为对矩阵Aj进行行SVD分解得到的主体成分,Djx为对矩阵Aj进行行SVD分解得到的细节成分,AjxAjy为对矩阵Ajx进行列SVD分解得到的主体成分,AjxDjy为对矩阵Ajx进行列SVD分解得到的细节成分,DjxAjy为对矩阵Djx进行列SVD分解得到的主体成分,DjxDjy为对矩阵Djx进行列SVD分解得到的细节成分。
Aj→Aj+1、D1j+1、D2j+1、D3j+1的算法如图4所示。其中,Aj为j尺度下信号的主体成分,通过一维H矩阵的构造方法来构造Aj进行行递推分解的构造矩阵Hjx和列递推分解的构造矩阵Hjy;Hjx,1、Hjx,2分别为Hjx的主体成分和细节成分。Hjy1、Hjy2分别为通过行递推主体部分和细节部分的列构造矩阵。
令
式中,M为二维离散信号的行数;N为二维离散信号的列数。
构造Aj进行行递推分解的构造矩阵Hjx,并进行SVD分解如下:
其中,int()为向下取整函数。且i=1,2,…,2 M;k=1,2,…,N-1。
对Hjx进行SVD分解如下:
令
令
其中,i=1,2,…,2 M;k=1,2,…,N-1。
通过式(3)和式(4)得出Ajx和Djx:
根据Axj和Dxj构造进行列递推分解的,并分别进行SVD分解如下:
分别令
其中,i=1,2,…,M-1;k=1,2,…,2 N。
分别令
其中,i=1,2,…,M-1;k=1,2,…,2 N。
通过式(3)和式(4)得出Aj+1、D1j+1、D2j+1和D3j+1如下:
由式(2)~式(4)和图3、图4可知:
由对一维离散信号的MSVD分解结果可知,一维离散信号通过SVD迭代分解最终得到的结果是稳定的,且不存在边界效应。二维离散信号的SVD迭代分解实质上是通过行分解和列分解迭代的形式进行的,因此二维离散信号MSVD最终得到的结果是稳定的,且不存在边界效应。
3 基于MSVD的表面滤波法
3.1 表面形貌滤波的MSVD模型
令z(x)为二维表面轮廓信号,z1(x)为二维表面形状误差,z2(x)为二维表面波纹度,z3(x)为二维表面粗糙度;z(x,y)为三维表面轮廓信号,z1(x,y)为三维表面形状误差,z2(x,y)为三维表面波纹度,z3(x,y)为三维表面粗糙度。二维表面形貌滤波的MSVD模型可表示为
其中,j为分解层数,z1(x)+z2(x)为表面评定基准线,Aj、Dk的求解见图2。三维表面形貌滤波的MSVD模型可表示为
其中,z1(x,y)+z2(x,y)为表面评定基准面,Aj、Dk1、Dk2、Dk3的求解见图4。
3.2 MSVD模型分解层数的确定
MSVD具有多尺度性,因此可在不同尺度下对表面特征进行分离和提取。对于具有多尺度性的小波分析应用于提取表面形貌评定基准时,存在小波分解次数的确定问题,同理,MSVD用于提取表面形貌评定基准也存在分解层数确定的问题。为此,本文提出了一种基于最大临近层差值条件的方法来确定分解层数。定义二维表面和三维表面最大一维相邻层差值:
其中,Δ(j)为二分递推的第j层Aj与第j+1层Aj+1的最大差值。设定最大邻层差值为Δ,当j≥k时,满足Δ(j)≤Δ,则确定递推分解终止层数为k。因为{Aj}是收敛的,所以必然存在k,使得j≥k时,Δ(j)≤Δ。
3.3 基于MSVD滤波法的计算复杂度
设一维测量信号为X1×N,T11、T12、T13分别为Aj→Hj、Hj→Hj,1和Hj,1→Aj+1的计算复杂度,k为分解层数,则一维MSVD滤波法的计算复杂度T1为
设二维测量信号为XM×N,T21、T22、T23、T24、T25、T26分别为Aj→Hjx、Hjx→Hjx,1、Hjx,1→Ajx、Ajx→Hjy1、Hjy1→Hjy1,1、Hjy1,1→Aj+1的计算复杂度,k为分解层数,则二维MSVD滤波法的计算复杂度T2为
4 仿真计算和实验验证
4.1 一维MSVD的实测实验
通过光切显微镜测量长方形铜块表面轮廓,该实验的数据评定长度为ln=7.5mm,采样间隔为Δx=0.05mm,采样点数为150,如图5所示。对该实验数据进行一维MSVD,设定迭代终止条件Δ=0.05,由图6可知分解层数可选17。
为了验证该方法的可行性与正确性,本文分别采用了高斯滤波法和高斯回归滤波法来提取实验数据的表面基线(图5、图7)。图8为基于一维MSVD滤波法得到的表面基线。从图5、图7、图8可知三种滤波方式得到的表面基线有较好的一致性,并且本文算法不会存在高斯滤波法所存在的边界数据丢失的问题。
4.2 二维MSVD仿真实验
为了说明二维MSVD对三维表面滤波法的有效性,本文以多峰函数peaks(200)作为一个采样点数为200×200的三维基准面,同时在该基准面上加白噪声作为模拟表面,如图9所示。运用基于二维MSVD滤波法对模拟表面进行滤波处理,设置迭代终止条件Δ=0.05,得到的基准面如图10所示。同时本文还分别采取了高斯回归滤波法和双树复小波滤波法来提取基准面,如图11、图12所示。
通过对仿真数据进行对比发现:从基准面的形貌来说,三种滤波结果具有较好的一致性。不过滤波性能还是有一定的差异:采用高斯回归滤波时,得到的基准面不够自然平滑,与理想基准面有一定的差异。采用双树复小波滤波法和基于二维MSVD滤波法得到的基准面自然光滑且不存在边界效应。
为了进一步评价不同滤波法对仿真信号的滤波效果,本文采用信噪比和最大误差两个评价指标来评定。其中二维信噪比RSNR为理想基准面功率和滤波后基准面的噪声功率之比:
最大误差Δmax为理想基准面与滤波后基准面的最大误差:
其中,z为标准基准面的高度矩阵,z1为滤波后基准面的高度矩阵。通过不同滤波法提取基准面的评价指标如表1所示。
从表1可以看出,高斯回归滤波法信噪比最小和最大误差最大,从而评价结果最差;基于二维MSVD滤波模型得到基准面的信噪比最大并且最大误差最小,从而评价结果最优。
4.3 二维MSVD实测数据分析
采用英国Taylor-Hobson公司的CCI型白光干涉式表面测量仪进行数据分析,其最小采样间距为0.078μm,垂直分辨率为0.01nm,实验样品为某磨削工件,采样点数为256×256,采样间距Δx=1μm,截止波长选择为λc=80μm,实测表面如图13所示。为了验证方法的可行性与正确性,本文分别采取高斯回归滤波法、双树复小波滤波法和二维MSVD滤波法对磨削样品进行滤波处理,结果如图14~图16所示。
由图14~图16可知,三种滤波法能有效地提取表面基准面,并且基准面的形貌具有较好的一致性。但高斯回归滤波结果受异常信号的影响有明显的凸峰和凹谷,而双树复小波滤波法和基于二维MSVD滤波法具有多尺度性,可通过选取较大尺度来抑制凸峰和凹谷,但双树复小波滤波法相邻尺度滤波结果差异性较大,并且分解尺度有限。为保证滤波结果的准确性,双树复小波分解尺度可根据截止波长或能量守恒法来确定[11,12,13],由采样间距和截止波长计算双树复小波分解层数为5,其滤波结果见图15。基于二维MSVD滤波法因细节成分是公比小于1的等比数列,所以相邻尺度滤波结果差异随分解尺度的增大而减小。可通过改变最大差值终止条件来微调表面基准面,这是双树复小波滤波法所不具有的。本文设置最大差值终止条件为0.05,其滤波结果如图16所示,从图16可知基于二维MSVD滤波得到基准面存在不明显凸峰和凹谷,在一定程度上能抑制异常信号的影响。
5 结束语