奇异值分解

奇异值分解（通用8篇）

奇异值分解 篇1

0 引言

学生成绩是学籍管理中的重要内容之一, 它是学校评价学生的最主要的指标之一[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 结语

本文为了有效地对学生成绩数据进行恢复, 提出了一种迭代奇异值分解的学生成绩恢复方法。该方法将学生成绩表示为一个矩阵, 利用该矩阵具有低秩的特性, 在给定缺失元素的初始值后, 采用迭代奇异值分解的方法, 最终求取成绩表中缺失元素的真实值。模拟实验和真实实验结果表明该方法能够快速、正确地恢复出学生的真实成绩。

摘要：为了有效地对学生成绩数据进行恢复, 提出一种迭代奇异值分解的学生成绩恢复方法。该方法采用矩阵表示学生成绩, 利用该矩阵具有低秩的特性, 在给定缺失元素的初始值后, 利用奇异值分解得到缺失元素的近似值, 而该近似值比初始值更加接近真实值。再将求到的近似值代替初始值, 经过多次迭代, 最终可求到成绩表中缺失元素的真实值。该方法的优点是在缺失元素恢复过程中, 利用了所有已知元素信息, 并将所有已知元素平等地对待。模拟实验和真实实验结果 表明该方法能够快速、精确地恢复出学生的真实成绩。

关键词：学生成绩,缺失元素,低秩矩阵,奇异值分解

奇异值分解 篇2

奇异二阶四点边值问题的正解

研究了如下奇异二阶四点边值问题{un(t)+h(t)f(t.u(t))=0,0＜t＜1,u(0)=au(ξ).u(1)=βu(η),其中0＜ξ＜η＜1.0≤α.β＜1,f(t,u)在u=0点具有奇性,h(t)在t=0.1点具有奇性.应用Green函数的.性质和存在性原则,得到了正解的存在性.

作 者：苗春梅 葛渭高 MIAO Chun-mei GE Wei-gao  作者单位：苗春梅,MIAO Chun-mei(北京理工大学理学院,数学系,北京,100081;长春大学,理学院,吉林,长春,130022)

葛渭高,GE Wei-gao(北京理工大学理学院,数学系,北京,100081)

刊 名：数学的实践与认识  ISTIC PKU英文刊名：MATHEMATICS IN PRACTICE AND THEORY 年，卷(期)：2008 38(13) 分类号：O1 关键词：奇异   四点边值问题   正解   Green函数  

奇异值分解 篇3

随着多媒体和互联网技术的发展,数字图像信息的应用越来越广泛,要求大量存储和传输图像。未压缩的数字图像信息量过于庞大,通过网络传输图像信息会占用大量的网络资源,往往要求在保证质量的前提下以较小的空间存储图像,以较少的比特率传输图像,因此采用合适的方法对图像进行压缩,以便图像的存储和传输,具有重要的意义。

奇异值分解(SVD,Singular Value Decomposition)是线性代数最有用的工具之一,广泛应用于应用数学、信号处理和模式识别等诸多领域。它是通过对矩阵的分解和近似,有效地把任意矩阵降低到一个较小的可逆方阵,然后再利用可逆矩阵进行数据分析,从而达到简化分析的目的。

1 图像压缩原理

图像数据压缩的目的是在满足一定图像质量条件下,用尽可能少的比特数来表示原始图像,以提高图像传输的效率和减少图像存储的容量。压缩编码节省了数据的存储空间,这样不论是在传输数据还是在处理数据的时候都会带来很大的便利,图像数据压缩在信息论中称为信源编码[1]。

去除多余数据以数学的观点来看,这一过程实际上就是将二维像素阵列变换为一个在统计上无关联的数据集合,图像压缩是指以较少的比特有损或无损地表示原来的像素矩阵的技术,也称图像编码。

图像数据之所以能被压缩,就是因为数据中存在着冗余,图像数据的冗余主要表现为:图像中相邻像素间的相关性引起的空间冗余;图像序列中不同帧之间存在相关性引起的时间冗余;不同彩色平面或频谱带的相关性引起的频谱冗余。数据压缩的目的就是通过去除这些数据冗余来减少表示数据所需的比特数。由于图像数据量的庞大,在存储、传输、处理时非常困难,因此图像数据的压缩就显得非常重要[2,3]。

2 奇异值分解

奇异值分解是线性代数中一种重要的矩阵分解,是矩阵分析中正规矩阵酉对角化的推广。在信号处理、统计学等领域有重要应用。

奇异值分解在某些方面与对称矩阵或Hermite矩阵基于特征向量的对角化类似。然而这两种矩阵分解尽管有其相关性,但还是有明显的不同。对称阵特征向量分解的基础是谱分析,而奇异值分解则是谱分析理论在任意矩阵上的推广。

在目前的图像编码研究领域中,虽然存在着种类繁多的编码方法,但是他们都或多或少地遭遇到其自身发展的瓶颈,于是不断开拓新的图像压缩编码方法自然成为这个信息时代的必然选择。SVD方法具有编码效率不确定性、数据矩阵复杂性等等特点。SVD方法在图像压缩领域有巨大的应用潜力[4,5]。

3 奇异值分解的图像压缩原理

矩阵奇异值分解在图像处理中有着重要应用。假设一幅图像有n×n个像素,如果将这个n2数据一起传送,数据量往往会很大。因此,希望能够改为传送一些较少的数据,并且在接收端还能够利用这些传送的数据重构原图像。

假设用n×n维矩阵A表示要传送的原始图像。假定对矩阵A进行奇异值分解,便得到

A=U∑VT. (1)

其中,奇异值按照从大到小的顺序排列。如果从中选择k个大奇异值以及与这些奇异值相对应的左右奇异向量逼近原图像,便可以使用k(2n+1)个数值代替原来的n×n个图像数据。这k(2n+1)个被选择的新数据是矩阵A的前k个奇异值、n×n左奇异向量矩阵U的前k列和n×n右奇异向量矩阵V的前k列的元素。

undefined

式(2)称为图像的压缩比。显然,被选择的大奇异值的个数k应该满足条件k(2n+1)

图1是有效压缩曲线图,图中,当选定n后,对应的k值只有选择在曲线下方才有效,所以图1也是k值的有效性边界曲线。

显而易见,若k值偏小,则压缩比偏大,重构的图像质量较差。反之,过大的k值又会导致压缩比过小,降低图像压缩的效率。

由上述分析可知,在传送图像的过程中,就无需传送n2个原始数据,而只需要传送k(2n+1)个有关奇异值和奇异向量的数据即可。另外,在接收端接收到奇异值σ1,σ2,…,σk以及左奇异向量u1,u2,…,uk和右奇异向量v1,v2,…,vk后,即可通过公式(2)重构出原图像。

undefined

图2给出了SVD图像压缩方法的流程图。

4 仿真分析

仿真中,将一幅卡通图像变成矩阵的形式,然后对矩阵进行奇异值分解,在不影响原图像品质的情况下,缩小图像的数据量,选取适当的k值,然后再将分解后的矩阵转换成图像,来达到压缩图像的目的。图3是一个大小为5.9 kB的jpg图像,利用提出的SVD方法对其进行压缩。其压缩结果如图3～6所示。

k=5压缩为3.38 kB,k=20为4.36 kB,k=30为4.52 kB,可以看出k越小,压缩后的图像越小,但图像不清晰,随着k的变大,图像变清晰,但压缩率也变小,满足奇异值的大L曲线。

5 结论

图像压缩的最终目的是达到减少存储资源、提高传输效率。为此,本文利用矩阵的奇异值分解这一数学工具,实现了图像的压缩与恢复,在选择忽略部分较小奇异值的基础上,可以达到很高的压缩比,同时不产生图像失真。仿真分析表明,矩阵的奇异值分解压缩方法具有较好的压缩效率。

参考文献

[1]DaVidSalomon.数据压缩原理与应用[M].第2版.吴乐南,等译.北京:电子工业出版社,2003.

[2]刘榴娣,刘明奇,党长民.实用数字图像处理[M].北京:北京理工大学出版社,2001.

[3]冯玉氓,邵玉明,张星.数据图像压缩编码[M].北京:中国铁道出版社,1993.

[4]胡乡峰,卫金茂.基于奇异值分解(SVD)的图像压缩[J].东北师大学报,2006,38(3):36-39.

奇异值分解 篇4

一、奇异值分解

奇异值分解是线性代数中非常重要的一种矩阵分解。该理论的诞生已经有百余年的历史, 随着信息工程的需求和计算机技术发展, 它被广泛地应用到统计分析、信号与图像处理、系统理论和控制等领域中[2]。

二、基于奇异值分解的非带限信号采样算法

下面通过对一种非带限信号——周期狄拉克流的分析来推导基于奇异值分解的非带限信号采样算法。

设x (t) 是周期为τ的狄拉克流信号, 其表达式为:

式中Ck为狄拉克流的权值, tk为狄拉克流的位置, 则x (t) 的傅里叶变换为

对该狄拉克流信号x (t) , 取hB (t) =Bsinc (Bt) 作为采样核, 经过采样核后得到采样值yn, n=0, 1, ..., N-1, 则采样值yn是原信号x (t) 的一个充分的描述。

利用yn重建原信号x (t) 的过程如下。

(1) 通过yn确定x (t) 的傅里叶变换X[m], |m|≤M;

式中HB是hB (t) 的傅里叶变换。

(2) 利用X[m]构建Hankel矩阵X, 然后对其进行奇异值分解, 求得狄拉克流的K个位置信息

利用X[m]构建一个P×Q (P, Q≥K) 维的Hankel矩阵, 矩阵X的秩为rank (X) =K, 因此通过对Hankel矩阵X的秩的判别, 可以确定出每周期内狄拉克流信号的个数K。

根据奇异值分解的理论, X可以被分解成如下形式

其中US、SS和VS是矩阵X的K个最大奇异值分解三对组阵, 它们包含了有用信号的主要信息;Un、Sn和Vn是矩阵X的小奇异值分解三对组阵, 对于含有噪声的信号, 其加性噪声主要集中在这些小奇异值项上。

假设φ=diag (Zk) k×k, φh=φH, 显然矩阵US和VS都满足平移不变空间特性, 即

(3) 求取狄拉克流的K个权值

三、结论

本文以周期狄拉克流信号为例, 研究了一种基于奇异值分解的非带限信号采样与重建算法。它突破了Shannon采样定理中Nyquist率的限制, 是传统采样定理的一个非常有益的补充, 在宽带通信领域尤其是超宽带通信中有着广阔的应用前景。

摘要：本文借助奇异值分解技术, 研究了对狄拉克流这类特殊非带限信号的采样与重建算法。先对信号的采样点进行离散傅里叶变换, 并用生成的系数形成Hankel矩阵, 然后对Hankel矩阵进行奇异值分解, 求取狄拉克流的位置信息, 再解范德蒙方程组求得狄拉克流的权信息, 重建出原信号。该算法具有计算量小、信号恢复精确率高和抗噪声能力强的特点。

关键词：奇异值分解,非带限信号,狄拉克流,采样,重建

参考文献

[1]Chang C C, Tsai P Y, Lin C C.SVD-based digital image watermarking scheme[J].Pattern Recognition Letters, 2005, 10:1577-1586.

奇异值分解 篇5

表面形貌是指零件在加工过程中残留下来的微观几何形态,可用粗糙度、波纹度、形状误差等来表征。表面形貌客观地反映了表面生成机制,影响工程表面的功能特性[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 结束语

奇异值分解 篇6

在图像处理中, 常常需要对图像进行滤波分解处理, 使感兴趣的目标信息与其他信息分离, 以突出目标信息, 获取特征。传统的彩色图像依据灰度图像理论, 将图像分成3个分量, 分别进行滤波。然而, 彩色图像的3分量原本是有机的统一整体, 相互之间具有较强的相关性, 因此, 如果人为地将其分开处理, 将会对图像整体结构信息造成影响。

1996年, Pei第一次提出了彩色图像的四元数模型。即:

q (x, y) =r (x, y) i+g (x, y) j+b (x, y) k

其中r (x, y) 、g (x, y) 和b (x, y) 分别表示图像中在 (x, y) 点处的红、绿、蓝色的灰度值 (或者是H、I、S的灰度值) 。基于此模型, 近年来有大量的研究将四元数及四元数矩阵作为一种重要的彩色图像处理的工具, 在图像的压缩、去噪和增强等领域得到广泛应用。

1 四元数矩阵分解

定理1 (SVDQ) 对于任一个秩为r的四元数矩阵A∈QN×M, 必有两个四元数酉矩阵U, V使得

其中, H表示共轭转秩:

Σ=diag (σ1, σ2, ..., σr, 0, ..., 0) , 为四元数矩阵A的奇异值且σ1≥σ2≥...≥σr>σr+1=...=0

U= (u1, u2, ..., uN) ∈QN×N, 为四元数矩阵A的左奇异特征向量

V= (v1, v2, ..., vM) ∈QM×M, 为四元数矩阵A的右奇异特征向量

Ψi=uiviH, 所有的Ψi构成A的正交完备基, σrΨi称为A的第i级特征子图像

根据矩阵的相关知识, 一个矩阵A的能量E可以由它的Frobenius范数表示:即

同时, Σ是实对角矩阵故有, 即EA=Σσi

可见, 图像大部分能量集中在前面较大特征值所对应的特征图像当中。

2 基于四元数奇异值分解的彩色图像重构

通过四元数奇异值分解, 可以得到一幅彩色图像的一系列奇异值。在重构过程中, 拥有较大数量级的奇异值表征的是图像的轮廓, 而较小的奇异值所表征的则是图像的细节信息。式 (1) 意味着, 如果选择全部奇异值和特征图像, 则图像可以完整重建;如果选择部分奇异值和特征图像, 则图像部分重建。基于四元数奇异值分解的彩色图像重构可由下式表明:

其中, Σk是保留Σ的前k个奇异值而将后面n-k奇异值置为0所得到的对角矩阵。但是, 怎样选取奇异值及其特征图像, 往往凭经验而定。还有利用相邻特征图像之间的差异性, 通过某一度量指标进行度量确定k值, 但对于拥有大量特征图像的高维图像, 这样做处理量大、速度慢。由于前k个子图重建的图像能量为:

同时σi与Ek之间具有分形规律, 即是有:

式 (4) 在双对数图上表现为lnα, 如果图像包含多个分形规律, 则双对数图上有相应条数β的直线, 如图 (c) 。为求得不同直线的拐点, 将 (logσk, log Ek) 有序分组, 用最小均方误差 (LS) 拟合多段直线, 使总误差平方和最小 (即总残差平方和最小) , 从而求得重建图像的k值。

3 研究实验

我们在MATLAB R2008a的平台下, 以Lena (140×140) 为例, 采用矩阵算法分解四元数矩阵进行处理。

说明:图 (a) 是加入均值为0, 方差为15的高斯噪声后的图像;图 (b) 是奇异值累计百分比图;图 (c) 是奇异值与能量的双对数图, I段是细节信息即高频信息, II段是概略信息即低频信息;图 (d) 是能量累计70%的重建图, 此时K=20;图 (e) 是能量累计80%的重建图, 此时K=37;图 (f) 是根据双对数图似合直线求得拐点的重建图像, 此时拐占K=24。

传统上对于k值的选择, 都是凭借个人经验及感观, 反复尝试选择某一个k值, 很难说清楚k值选择的合理性;运用本文的方法, 可以较好解释k值的选择:图 (d) 和图 (f) 比较, 尽管两者k值相近, 但图 (f) 更清晰些;图 (e) 和图 (f) 比较, 两者相差不大, 但图 (e) 的k值远比图 (f) 大了很多。文献[3]也从能量的角度确定k值, 但需要对每个子图再次进行分块求奇异值, 得到相邻子图的奇异值之间的全局误差序列图, 从而选出拐点确定k值, 这样多次进行分解效率很低, 在图像很大时特明显, 本文一次分解后, 迅速用双对数图圈出拐点确定k值, 效率更高些。

4 结束语

将四元数奇异值分解理论应用数字图像处理领域, 能较好保持彩色图像的信息结构完整性;尽管传统方法及其它确定k值的方法有其特有优点, 但实验表明, 本文能较快确定拐点, 拐点位置也有较明确的物理含义, 对于确定k值有很大的借鉴意义。

摘要：基于四元数矩阵彩色图像奇异值分解, 得到表征彩色图像的不同分量的奇异值, 运用分形快速确定图像的拐点。该方法具有快速和简单可行的优点。以受噪声污染图像为例, 该方法针对彩色图像去噪具有较好的效果。

关键词：彩色图像,四元数奇异值分解,分形

参考文献

[1]PEI S-C, CHENG C-M.A novel block truncation coding of color images by using quaternion moment preserving principle.IEEE In-ternational Symposium on Circuit s and Systems, Atlanta, 1996.

[2]ZHANG F.Quaternion and matrices of quaternion.Linear algebra and its applications, 1997, 251 (1) :21-570.

[3]雷印杰, 余艳梅, 周激流, 等.四元数奇异值分解与彩色图像去噪[J].四川大学学报 (自然科学版) , 2007 (6) .

[4]黄锐.基于矩阵分解的遥感图像处理方法[D].武汉:中国地质大学, 硕士学位论文, 2008.

[5]李庆谋, 成秋明.分形奇异 (特征) 值分解方法与地球物理和地球化学异常重建[J].地球科学, 2004 (1) .

奇异值分解 篇7

关键词：混沌序列,小波变换,奇异值分解,鲁棒性

1、引言

近年来, 数字化技术和Internet的飞速发展, 在最大限度的方便人类信息交流的同时, 也带来了版权的危机[1]。数字水印具有双重安全性, 即水印的添加与否具有不可知性以及水印提取受密钥的保护, 因而非常适用于信息安全问题。该技术为知识产权保护提供了一种新的手段[2,3]。根据水印加载方法的不同, 可以分为:空间域水印和变换域水印。本文主要讨论的是变换域上的水印技术, 与空间域水印方法比较, 变换域水印方法具有如下优点:1) 在变换域中嵌入的水印信号能量可以散布到载体的所有位置上, 有利于保证水印的不可察觉性;2) 在变换域, 人类视觉系统和听觉系统的某些特性 (如频率掩蔽效应) 可以更方便地结合到水印加载过程中;3) 一些变换域的方法可与数据压缩标准相兼容, 从而实现在压缩域内的水印算法, 同时也能抵抗相应的有损压缩。

2、混沌序列数字水印

混沌系统内在的随机性和对初始条件十分敏感的特性, 使得混沌理论在近十年被广泛应用于数字水印领域。混沌序列数字水印是指依靠混沌序列生成水印或者加密水印的方式。混沌序列是一种非线性序列, 结构复杂, 难以预测, 它具有比较好的特性:生成形式简单, 对初始条件极其敏感和类白噪声[4]。利用初值可以精确的重构混沌序列。我们借助混沌的这一特性, 利用Logistic混沌映射来进行预处理, 将进一步提高秘密信息的安全性。

2.1 Logistic映射

Logistic映射, 也称为虫口模型。Logistic混沌序列的遍历统计特性近似于零均值白噪声, 具有良好的随机性、相关性和复杂性 (其具体描述及相关特性详见参考文献[5]) 。

2.2 广义猫映射

混沌系统具有对初始条件的极端敏感性和类随机性, 使之具有很好的密码学特征;因此, 近年来混沌系统也被应用到密码学中。猫映射矩阵形式如下所示:

将上述猫映射作如下推广, 首先将相空间推广为{0, 1, 2, …, N-1}×{0, 1, 2, …, N-1}, 即只取0到N-1的正整数;其次, 将方程推广为最一般的二维可逆保面积方程:

式中a, b, c, d为正整数, 其保面积要求|C|=ad-bc=1。在此要求下a, b, c, d中四个参数只有三个是独立的, 比如我们可以让a, b, c独立, 而d由保面积条件确定。推广的猫映射仍然具有混沌映射的特性。因此, 我们利用广义猫映射来产生水印各像素点嵌入到载体图像小波系数的对角矩阵中的位置, 这样, 水印的全部像素将被随机而均匀地置乱到载体子图的对角矩阵的整个系数空间。

3、小波变换

小波变换的基本思想就是对信号进行细致的频率分离即多分辨率分解。通过一级小波变换, 原始图像被分解为4个一级子图。若对一级近似子图再进行小波分解又可得更低分辨率的4个二级子图;如此反复可对数字图像进行多层小波分解。其中最深层的低频子图集中了被分解图像的绝大部分信息, 刻画了图像的主体特征, 所以称为被分解图像的近似子图;低频子图抗外来影响的能力好。而各层高频子图则分别保持了被分解图像各方向的边缘细节, 刻画了被分解图像的边缘细节特征, 故统称为被分解图像的细节子图;但高频子图这些边缘细节易受外来噪声、常规图像处理等因素影响, 其稳定性较差。

4、奇异值分解

奇异值分解是线性代数中一种重要的矩阵分解, 它在最优化问题、特征值问题、最小二乘方问题、广义逆矩阵问题、信号处理和统计学等领域有重要应用。奇异值分解是通过svd函数实现, 具体定义如下所示:[U, S, V]=svd (A) 返回一个与A同大小的对角矩阵S, 两个矩阵U和V, 且满足=U*S*V'。

5、水印的嵌入提取及检测

5.1 水印的嵌入

Step1:将图像读取到相应矩阵中;

Step2:使用混沌Logistic映射对水印信息进行加密处理;

Step3:指定广义猫映射的独立参数a, b, c, 迭代次数n;

Step4:小波变换, 首先对原始数字图像进行L级小波分解, 得到不同分辨率级下的多个细节子图和第L级的一个逼近子图;

Step5:对第L级的一个逼近子图进行奇异值分解, 得到对角矩阵;

Step6:扫描水印图像一遍找出水印优值b0 (优值b0为像素最多的那个) , 确定小数下、上阈值dc1、dc2及广义猫映射中的参数d, 本算法中, 取dc1=4, dc2=6;

Step7:对水印图像各像素位置 (x=1~Mm, y=1~Nm, Mm=Nm) 做如下操作:

(1) 首先由猫映射公式及给定的独立参数a、b、c和迭代次数n, 确定该水印像素在对角矩阵中的嵌入位置 (x’, y’) ;然后求出对角矩阵 (x’, y’) 位置处的个位数部分值 (image假设为对角矩阵) c (x’, y’) ;

(2) 如果水印像素值message (x, y) =b0, 则修改相应对角矩阵的个位数部分到[dc1, dc2]范围内;

如果水印像素值为劣值message (x, y) =1-b0, 则修改相应对角矩阵使其个位数变为零;

Step8:对嵌入水印的对角矩阵再次进行奇异值分解, 通过奇异值的逆运算得到载体图像的低频系数;

Step9:小波逆变换, 将修改过的逼近子图结合各细节子图进行小波逆变换, 得到嵌有水印的隐秘载体图像;

嵌入算法的运行效果如图1所示:

5.2 水印提取

水印的提取是水印嵌入的逆过程:

Step1:使用小波变换得到水印的L级低频子图;

Step2:利用奇异值分解得到对角矩阵, 通过对角矩阵求得一个中间矩阵;

Step3:应用混沌猫映射确定中间矩阵中水印的对应位置, 利用统计学方法统计出此像素点上的值, 从而提取出加密的水印图像;

Step4:通过Logistic映射解密, 得到我们提取的最后水印图像;

提取算法的运行效果

5.3 水印攻击与检测

6、结论

本文研究了一种基于奇异值分解及小波变换的二值图像水印算法, 同时水印的恢复基于改正的二值算法, 可以免疫于待测图像在有限范围内失真造成的影响, 这类似于畸变的离散数字信号通过重新抽样可以完全恢复的特点, 因而具有满意的鲁棒性。从实验结果可以看出:该算法对JPEG压缩, 椒盐噪声, 剪切等具有比较好的鲁棒性。而且, 可以通过调整个位数阈值dc1和dc2可以在水印的鲁棒性和不可感知性之间进行平衡或向某一指标倾斜。水印嵌入位置被混沌置乱, 增加了水印的安全性。水印的嵌入和提取基于相同的密钥 (即混沌参数、迭代次数和个位数阈值) , 使水印的恢复不需要二值逻辑矩阵的辅助;算法简单, 更具有实用性。

参考文献

[1]李明, 廖晓峰.结合混沌的小波变换数字水印技术[J].计算机科学, 2007, 34 (8) :245-247.

[2]许红山.置乱技术在信息隐藏中的应用。广州大学学报, 2004, 3 (2) :P134-136.

[3]黄这人, 刘九芬。黄继武.小波变换城图像水印嵌入对策和算法软件学报.2002.13 (7) :1291-1297.

[4]朱从旭, 陈志刚.基于DWT域的混沌置乱二值图像数字水印新算法[M].湖北:小型微型计算机系统, 2005, 26 (7) :2-3

奇异值分解 篇8

单相接地故障作为矿井电网最主要的故障类型,约占电气故障总数的70% 左右[3]。单相接地故障后,非故障相的对地电压升高为原来的倍。尽管煤矿井下同样采用小电流接地方法,可以在故障发生后持续运行1 ~ 2 h,但是要迅速找到故障点并排除故障,以免造成更严重后果。由于故障存在的过度电阻及电弧接地情况等不确定性及配电网的多分支等方面的影响,传统的故障定位方法很难准确识别故障点[4]。文献[5]提出的利用电压、电流行波的综合行波极性判别,并在行波法失效时辅以阻抗测距。

1 行波的基本原理

按照基本原理,故障定位可以基本分为两大类:阻抗法和行波法。阻抗法受线路结构不对称、过渡电阻、线路参数沿走廊分布不均匀以及互感器变换误差等因素的影响,测距误差较大[6]。行波法拥有诸多优点克服了阻抗法的不足,提升了定位精度。行波法又分为单端和双端行波之分。

1. 1 双端行波法

双端行波法原理简单,在故障发生后,只需分别测出由故障点产生的行波首波到达故障线路的两端的时刻即可求出故障点信息。如图1 所示为双端测距的原理图,假设F点发生接地故障,距离R,S两端的分别为d1和d2,各方向首波到达R,S的时间为t1和t2。

式( 1) 中 Δt = t1- t2,这样就得到故障点的位置。

GPS作为同步授时单元,随着算法及处理器等不断进步,授时精度达到了纳秒级,但是由于两端都采用授时单元会可能造成误差积累。行波传输速度之快接近与光速,些许时间误差会对测距结果产生很大的影响。文献[7]提出首先利用双端行波测距判断故障大致距离,然后基于该故障距离及线路自身参数关系来准确识别故障点反射波和对端母线反射波达到测量端的时刻,从而精确测得故障距离。此方法虽然克服了双端测距GPS授时精度影响,但是没有考虑速度变化带来影响。

作为文献[8]中,采用一种改进的双端定位法,在线路中间增加第三个采集点,从而消去速度,只含有各点的时间函数进行定位。此方法前提条件是速度不变,可是由于各频率行波衰减程度不同,速度肯定随着改变。即使不计成本,增加第三个点也增加了信号采集出现的误差,对最终结果的精确性造成一定的影响。

1. 2 单端行波法

单端行波法根据故障点产生行波及经过故障点反射产生行波到达一端的时间差定位。原理图如图2 所示。

当F点靠近R端发生接地故障的时候,会分别产生向R,S方向运动的行波。行波到达R端时刻为t1,因为R端的波阻抗发生变换,行波会发生折射与反射。由R端反射的行波到达故障点F,又发生折射与反射,当F点反射的行波到达R端的时刻为t2,这样就可以算得故障点F的位置信息:

若故障点F靠近S,则经过S点折射回来的行波会早于故障点再次折射的行波。其他条件不变,波形经过F-S-F-R记录到达R端时间为t'2,从而求得故障距离为:

文献[9,10]等文章已经解决如何判别第二次透射过来的波形是来自与对端母线还是故障点。文献[11]根据第二个到达测量点的反向行波零模与线模的极性相同则为故障点反射; 若极性相反则为对端母线透射过来得。此方法更为简单,计算量更小。

目前还没有准确适用于煤矿井下的单相接地故障定位方法。采用的双单端定位法减小了只依靠单端定位时的误差。也不需要双端法需要精确的授时单元,此外行波色散对波速变化的考虑也提高了传统测距的测量精度。

2 Hankel矩阵下的信号奇异性检测

行波波头到达母线的时刻是故障测距的关键之一,行波波头到达母线的信号会引起突变,或者其某阶导数的突变。发生突变的点称为奇异点,用Lipschitz指数来描述。一维信号构造Hankel矩阵,经过奇异值分解SVD( sigular value decomposition) 后获得相似于小波变换得到的奇异值信息。而小波变换一旦选定了小波基,奇异性指数则不会再发生变换[11]。每一端进行单端测距时,通过构造Hankel矩阵进行SVD分解,这样能获得比小波模极大值法更准确的首波与反射波到达的时间。

SVD是指一个实矩阵A ∈ Rm ×n,存在两个正交矩阵U = [u1,u2,…,um]∈ Rm ×m和V = [v1,v2,…,vn]∈ Rn ×n,满足式( 4)

当m > n时,∑ = [diag( σ1,σ2,…,σq) ,0]m ×n;m < n时候,∑ = [diag( σ1,σ2,…,σq) ,0]Tm ×n。其中σi被称为奇异值,而且存在 σ1≥ σ2≥ … ≥ σq> 0。则可以将A写成特征值与特征向量乘积形式:

式( 5) 中,q = min( m,n) 。令Ai= σiuiviT,则

假设一组离散信号X = [x( 1) ,x( 2) ,…,x( N) ],构造Hankel矩阵( 6)

式( 6) 中,1 < n < N,m + n - 1 = N。再由式( 5) 得到图3。

图3中Ri,1与Ci,n分别是矩阵Ai的第一行及最后一列除去首元素所剩下的部分。让Pi=[Ri,1,CTi,n];Ri,1∈R1×n,Ci,n∈R(m-1)×1。这样,原离散信号X就能用来表示。而Pi从第二层开始,就具备指示奇异值位置信息的能力。相对于小波极值法的优点是:保持了各分量信号在原信号中的相位不变,具有零相位偏移的特性[8]。在经过仿真对比发现,小波变换需要对数据进行截取,否则在信号初始阶段与结束阶段会出现大的越变导致所需信号很难辨认,给故障测距增加了难度。采用图6所示的系统进行仿真,若对故障时刻0.035~0.040 s母线端采集的电流信号进行小波分析,不进行信号截取,就会出现如图4所示现象。而Hankel矩阵下的SVD分解不存在这样的现象。

3 行波波速的确定

3. 1 波速与频率关系

煤矿井下使用的都是电力电缆,而电缆波速受到多种因素影响。不同的频带更是影响波速的大小,因此确定行波的波速是故障测距的另一个关键点。

在电力传输的过程中,线路的三相之间存在耦合现象,需要进行相应的解耦变换。经过变换后的0、1、2 模成为独立的单相系统,可以分别进行研究。常用的变换为Karenbaner和Clarke变换。根据文献[12],本文选取 α 模电流信号进行研究。

根据行波的相关理论可知,传播系数。其中,R,L,G,C分别为单位长度的电阻,电感,电导和电容,均为频率的函数。式中 α 为衰减系数,决定行波的衰减程度; β为相速度,与 ω 一起决定行波波速。行波波速可以用公式( 7) 表示。

可见,行波波速v随着频率f变化而变化。由于行波在线路传输过程中存在衰减,高频部分波速最快,衰减也最快。在给定条件下,得到了频率值就确定了行波波速值。电缆中,速度与频率的关系如图5 所示。

3. 2 中心频率

文中仿真用的一个简单的煤矿高压供电系统如图6 所示,供电端为10 k V电压,三条线路馈出,故障初相角- 90°,接地电阻为50 Ω,断路器闭合,系统经消弧线圈接地。

假设图5 所示系统,线路3 在t = 0. 035 s得时候发生故障,线路总长为15 km,故障点距离左端母线长度为9 km。故障时刻起,线路3 两端测量到的行波分别经过Clarke变换,得到母线端与负荷端的α 模电流行波。

由前文分析看出,频率对波速有很大的影响,因此频率的正确选取也对精确定位起到了至关重要的作用。小波分解能将原包含混叠在一起的信号分解为不同频带的信号,便于对原始信号在不同频带进行分析。经故障点折射或透射的电流行波会有很大的相似性,若在某尺度分解下通过计算首波与第二个波形的相关程度最大,则选取该尺度的计算中心频率[13]。

在对比两种在信号奇异性检测中,常用的db小波与Gauss小波,不同尺度下分别计算首波与第二个波形的相关系数。如表1 所示。

文献[14]也给出了在电力系统暂态信号分析中小波基的选择原则。本文根据两个波形相关度最大为标准,选择Gauss小波相应尺度的中心频率作为行波波速计算的频率。

4 双-单端定位方法

4. 1 误差修正

经过SVD与中心频率的确定,母线端与负荷端能够分别求得测量值d1与d2,假设故障点距两端的实际值为l1与l2。l为线路的总长度,则有l = l1+l2。那么测量值d1与d2的测量误差为式( 8) 所示。

将式( 8) 中上下两式相加得到式( 9) 。

令,这样得经过修正的故障点距离母线端d1',距离负荷端d2',有d1'+d2'=(d1-Δl/2)+(d2-Δl/2)=l。从整体上降低了误差,提高了故障测距的精度。

4. 2 故障测距步骤

应用的故障测距方法,总结起来分为以下三个步骤:

( 1) 母线端与负荷端分别采集电流信号,构造Hankel矩阵,采用SVD分析信号,用P3分量进行奇异值点描述。通过最大值搜索得到首波与第二个波形到达的时间。采用文献[15]中的方法判断第二个波形来自故障点直接折射还是对端透射过来的。从而确定是近端故障还是远端故障。

( 2) 对采集的信号在以尺度间隔0. 5 的gaus1小波上,对首波及第二个波形进行相关性分析,采用相关度最大的尺度为准,求取中心频率来计算行波波速。

( 3) 利用4. 1 节中的方法进行结果修正,得到更为准确的故障距离。

5 仿真分析

以图6 所示的模型进行仿真,仿真条件也如前文所述。在gaus1 小波2. 0 尺度上的前两个波形相似度最大,这时该尺度下小波的中心频率为0. 100 ×106Hz,根据公式( 7 ) 求得行波波速为v = 1. 989 ×108m / s。根据图4 可以求得母线端前两个波形到达的时间差为( 106 - 46) × 10- 6s。这样求得的故障点距离母线端长度为d1= 9. 032 km。

负荷端采用相同的方法,求得故障点距离负荷端长度为d2= 6. 069 km,Δl = 0. 101 km。这样求得母线端到故障点的修正距离为d'1= d1- Δl /2 =8. 981 5 km。修正前相对误差为0. 356% ,经过修正后的相对误差为0. 206% ,精度提高了0. 15% 。

传统的小波极值法在没有考虑色散对波速影响的情况下,仿真后测得母线到故障点的距离为:8. 929 km,相对误差为0. 789% 。不仅高于本文未经修正结果的相对误差0. 356% ,更高于经过修正的相对误差值。

还有一种概率非常低的情况,即当两端产生的误差相互抵消时,就会出现零修正,从前面的例子看出,即使没有修正的结果准确性也高于传统的小波极值法。

两个终端测距,还能对最终结果可靠性进行判断。实际应用时,现场设备以及煤矿井下的环境复杂性带来了很多未知因素。以母线端结果为参考,修正前相对误差为eo,修正后的相对误差ec,可以设置相对误差阈值er。若|ec- eo|≥ er则至少某一端出现了判断失误,需要用其他的测距方法来补充。

经过验证,当断路器打开时,系统以中性点不接地方式运行时,本文方法判断的故障距离跟中性点接地结果一样。

其他条件不变,只改变故障点位置,可以看到表2 的数据统计与传统小波极值法的比较。

从表2 可以看出本文提出的方法在精度上要高于传统的小波模极值法。

接地时的故障电阻也会对行波阻抗造成一定影响,采用3. 2 节中的仿真条件,只更改接地电阻大小,故障点到母线端实际距离还是9 km,仿真测量值如表3。

本文方法以及传统的小波模极值法都不受接地故障初相角的影响。因为初相角改变并没有使影响行波波速的线路参数或频率产生变化。

6 结论

本文从影响行波准确性因素的时间与波速入手。首先构造Hankel矩阵进行SVD分解,得到优于小波变换的极值点时间。然后分析不同频率对波速的影响,从而采用小波变换后,首波与第二个行波波头相似度最大的尺度中心频率作为计算波速的频率。母线端与负荷端都经过相同的办法得到初值,最后再用误差纠正法得到更优的修正结果。仿真表明,本文方法较传统的小波极值法有更准确的结果,不受接地电阻、故障初相角及接地点位置等影响。

摘要：行波测距最重要的是时间和波速的确定。双端测距需要GPS授时系统,两端的授时误差增加了首波到达时刻的准确性。对电缆线路来说,行波的色散更加严重。而且不管单端还是双端测距,波速的确定对结果的影响都非常大。故障线路两端同时利用Hankel矩阵方式下对到达两端的电流行波进行奇异性分解(SVD)。每端分别得到首波及第二个波形到达的时间。SVD在不同的分解层上的奇异点结果不会发生偏移,克服小波变换随着分解层数增加带来的奇异点偏移缺点。波速的大小与频带关系紧密,行波信号经过小波分解再重构为多个频带的时域信号。在计算每个频带下前两个波头之间的相似度后,选取相似度最大的频带得到电缆的波速。在考虑线路时频特征对行波波速v的影响同时,提出双-单端故障测距。最后两端分别得到的故障点位置经过误差纠正算法,使得结果更加准确。