小波包算法

2024-10-25

小波包算法（精选7篇）

小波包算法篇1

0 引言

数字水印作为一种新兴的防止盗版的技术日益受到人们的关注。它是将数字签名、商标和安全码等信息作为水印嵌入到图像中,同时要求不引起原始图像质量的明显下降,而且对于常见的图像处理操作应具有稳健的特性[1]。

目前,常见的数字水印技术是在信号的空域[2]和变换域[3,4],而现在研究的主要方向是变换域嵌入水印的方法。基于离散小波变换的数字水印技术主要是利用小波分解的局部特性和多尺度分析的特性来对原始宿主图像的低频子带进行进一步分解,再结合人眼视觉特性(Human Visual System,HVS)把水印嵌入到不同层次的低频系数上。由于在低频系数嵌入水印不会提高噪声水平,而且对低通滤波、有损压缩等具有较强的鲁棒性,但是仅在小波域的低频系数中嵌入水印具有嵌入水印容量有限和嵌入水印后可视性较差等缺点。

为了解决小波域低频嵌入水印有限和嵌入后图像不可见性差的问题,提出了一种基于小波包分解的自适应水印嵌入方案。算法嵌入的主要思想是对原始图像进行二级小波包分解,将待嵌入的二值水印进行一级小波分解,将原始图像的低频区域进行分块,并且判断各子块是否为纹理块,采用不同的嵌入强度将水印的低频信息自适应的嵌入到原始图像的低频频带中。对于水印的中高频分量,采用位置自适应的加性嵌入方式嵌入到原始图像的中高频分量的低频部分中。在整个嵌入过程中充分的利用了Hilbert扫描曲线对图像处理的优势。实验结果表明,本算法很好地解决了鲁棒性和不可见性之间的矛盾,对常规图像处理具有较强的鲁棒性。

1 小波包分解和Hilbert扫描曲线

1.1 小波包分解

小波包分解能够为图像提供一种比小波分解更加精细的方法,具有更广泛的应用价值。通过一级小波变换,原始图像被分解为4个子带图像,若分别对4个子带图像进行进一步的小波分解,又可得到更高分辨率的4个子带图像。如此反复,便可对图像进行多级小波包分解。图1和图2分别为二级小波包分解的树形图和Barbara图像二级小波包分解的示意图。

1.2 Hilbert扫描曲线

在文献[5]中,Hilbert曲线的定义是通过正方形逐次分割的标号次序给出的。例如,一单位正方形(如图3所示),将其逐级分割为4个相等的小正方形,图3(a)、图3(b)分别为一、二级分割,以此类推。对于每个自然数k,对k-1次分割再做细分割,即得到4k个k级正方形,每个正方形的边长为1/2k,记这些正方形为Q $_{i}^{(k)}$ (i=1,2,…,4k),它们的标号次序满足以下条件:

① 每个Q $_{i}^{(k)}$ 包含初始点O(0,0);

② k级正方形Q $_{4 i - 3}^{(k)}$ ,Q $_{4 i - 2}^{(k)}$ ,Q $_{4 i - 1}^{(k)}$ ,Q $_{4 i}^{(k)}$ 包含在k-1级的正方形Q(k-1)i中(k>1,i=1,2,…,4(k-1));

③ 任何下标为连续整数的正方形Q $_{i}^{(k)}$ 与Q $_{i + 1}^{(k)}$ 具有公共边。

将以上标号按从小到大的顺序用折线连接起来,就得到了Hilbert曲线。C.Gostman和M.Lindenbaum证明它是所有扫描曲线中最好的保持空间局部邻接性的曲线,特别适用于图像处理。

2 水印的嵌入及提取算法

在介绍嵌入算法前,先介绍一下图像纹理块的确定依据[7]。

① 将宿主图像的低频频带划分为8×8的子块,并计算子块的平均灰度值;

② 用子块中每个像素点的灰度值减去该子块的平均灰度值,并将差值的绝对值做和,记为sum;

③ 根据①和②计算出各子块sum的最大值max和最小值min;

④ 将t=(max-min)×3/5+min作为判定子块是否为纹理块的阈值,若子块的sum大于t,则该子块是纹理块,用T表示。否则为非纹理块,用 $\tilde{Τ}$ 表示。

2.1 水印的嵌入算法

① 水印信号的生成。将水印图像映射成一维向量为:

Wt $= {W t (i) | W t (i) \in {0, 1}, 0 \leq i \leq n \times n}$ ,

对Wt进行一级小波分解。

② 对原始图像I进行二级小波包分解,将低频频带nd(5)分割成互不重叠、大小为8×8的子块,记为S1k(k=1,2,…,N),N为分块的总数。对S1k进行Hilbert扫描,得到Hilbert序列B1k。

③ 将Wt的低频信息M自适应地嵌入到Hilbert序列B1k中。

令 $F (i, j) = \max | B 1_{k} |$

${\begin{matrix} F^{'} (i, j) = F (i, j) + α_{1} \times Μ ‚ B 1_{k} \in Τ \\ F^{'} (i, j) = F (i, j) + α_{2} \times Μ ‚ B 1_{k} \in \tilde{Τ} \end{matrix}$

。

式中,F(i,j)为第k(k=1,2,…,N)个分块中绝对值最大的小波包系数;F′(i,j)为嵌入水印信息后的小波包系数;α1,α2为嵌入强度。

④ Wt的中频信息采用位置自适应嵌入方式,分别嵌入到二级小波包分解的中频分量的低频信息中,具体方法如下:

分别将nd(9),nd(13)分割成互不重叠、大小为8×8的子块,分别记为S2k,S3k。并对以上得到的子块序列进行Hilbert扫描,得到Hilbert序列B2k,B3k。令

$F (m) {i, j} = \max | B (m)_{k} |$ ,m=2,3。

F′(m){i,j}=F(m){i,j}+α3×N(m)。

式中,F(m){i,j}为第k个分块中绝对值最大的小波包系数;F′(m){i,j}为嵌入水印信息后的小波包系数;α3为嵌入强度;N(m)为Wt的中频信息。

⑤ Wt的高频信息采用步骤④的算法,直接嵌入到原始图像经一级小波包分解的高频信息nd(4)中,这是因为原始图像经一级小波包变换后高频信息量本来就很少,不宜再进行细分。

⑥ 对各嵌入水印的分块图像分别进行Hilbert扫描的逆变换,得到正常序列,并进行二级小波包逆操作,生成最终含有水印信息的合成图像IWt。

2.2 水印的提取算法

水印的提取是水印嵌入的逆过程。本文算法中,由于水印的低频部分和中高频部分采用不同的嵌入算法,所以提取过程中要分别对水印的低频和中高频进行提取,步骤如下:

① 对含水印图像IWt做二级小波包分解,对低频频带进行分块,采用2.1节水印的嵌入算法的逆算法提取水印的低频信息;

② 对含水印图像IWt做二级小波包分解,对中高频频带进行分块,采用2.1节水印的嵌入算法的逆算法提取水印的中高频信息;

③ 对提取出的水印的低频信息和中高频信息进行一级小波逆变换,生成水印信息Wt;

④ 将得到的序列Wt按一定的顺序排列,即得到提取的水印图像。

3 仿真实验及结果

仿真实验所用原始图像以512×512的标准灰度图像“Barbara”为载体图像,以32×32的二值图像“任职教育”作为水印。对水印的嵌入、提取和攻击进行仿真实验。图4(a)、图4(b)和图4(c)分别是Barbara原始图像、二值水印图像和Barbara图像嵌入水印后的图像。从主观上看,嵌入水印前的图像和嵌入水印后的图像的视觉效果并没有什么差别,说明了该算法具有较好的不可见性。实验中取α1=0.1,α2=0.08,α3=0.001。

在参考文献[8]中,以峰值信噪比(Peak Signal Noise Ratio,PSNR)评价原始图像与含水印图像之间的差别,主观上可以容忍的图像PSNR值都在20 dB以上。提取出的水印图像 $\tilde{W} t$ 与原始水印图像Wt之间的相似性,用NC值来表示,它的值越大表示2个图像Wt与 $\tilde{W} t$ 越相像。PSNR和NC的定义如下:

$Ρ S Ν R (d B) = 10 \lg \frac{D^{2} Ν Ν}{\sum_{x = 1}^{Ν} \sum_{y = 1}^{Ν} (Ι (x, y) - Ι_{W t} (x, y))^{2}}$ ;

$Ν C = \frac{\sum_{i = 1}^{n} \sum_{j = 1}^{n} W \tilde{t} (i, j) \oplus \bar{W} t (i, j)}{n \times n}$ 。

式中, $\bar{W} t$ 为Wt的逻辑非运算。

对其他一些纹理相对复杂的标准测试灰度图像进行了水印嵌入的不可见性实验,取得了很好的效果,并与文献[3]的算法得到的峰值信噪比(PSNR)的实验结果进行了比较,结果如表1所示。

最后是以Barbara图像为载体对含水印图像进行攻击后提取的水印信息,如图5(a)、(b)、(c)、(d)、(e)、(f)所示,分别为加高斯噪声(0.01)、加椒盐噪声(0.01)、进行中值滤波(3×3)和JPEG压缩(压缩因子分别为70,50,30)。其NC值分别为:0.840 8,0.842 8,0.614 3,0.996 1,0.953 1和0.859 4。

从实验结果来看,本文算法对常规的图像处理具有良好的鲁棒性,解决了水印的不可见性和鲁棒性之间的矛盾,且能嵌入较大容量的水印。

4 结束语

利用图像小波包变换系数的特点,将经一级小波分解后水印的低频信息和中高频信息采用不同嵌入算法分别嵌入到原始图像的不同位置,在嵌入过程中充分考虑了图像的自身内容和人类视觉特性。而Hilbert扫描可以很好地保持空间局部邻接性,本文充分结合了Hilbert曲线扫描和小波包分解的特性。实验结果表明,本算法确实是有效的,很好地解决了水印的不可见性和鲁棒性之间的矛盾,且能嵌入较大容量的水印,具有较强的实用意义。算法需要进一步解决的任务是加强算法抗几何攻击的研究。 

摘要：为了解决日益严重的盗版问题,提出了一种基于小波包变换的自适应水印算法。对宿主图像进行小波包分解,对水印图像进行小波分解。根据视觉感知特性将水印的低频分量自适应地嵌入到宿主图像的低频分量中,嵌入强度根据图像的内容自适应地决定。水印的中高频分量用位置自适应的方法嵌入到图像的中高频分量中。嵌入过程中充分利用Hilbert扫描曲线的良好特性。实验结果表明,该算法对常见的图像处理操作具有较强的鲁棒性。

关键词：数字水印,小波包,Hilbert扫描,人类视觉特性

参考文献

[1]蒋天发.网络信息安全及数字水印技术的研究[J].武汉理工大学学报(交通科学与工程版),2003,27(6):826-828.

[2]COX I J,MILLER M L,BLOOM J A.Digital Watermarking[M].Morgan Kaufmann Publishers,San Francisco,CA,2001:89-152.

[3]张建伟,鲍政,王顺风.图像小波域分块奇异值分解的自适应水印方案[J].中国图象图形学报,2007,12(5):811-818.

[4]张旭,张贵仓.一种基于小波包分解的自适应数字水印算法[J].计算机辅助设计与图形学学报,2007,19(7):931-934.

[5]齐东旭.分形及其计算机生成[M].北京:科学出版社,1994:15-18.

[6]王笋,徐小双.Hilbert曲线扫描矩阵的生成算法及其MATLAB程序代码[J].中国图象图形学报,2006,11(1):119-122.

[7]刘卓,王相海.一种基于提升方案小波和纹理块的图像鲁棒水印算法[J].微电子学与计算机,2006,23(6):20-23.

[8]孙圣和,陆哲明,牛夏牧.数字水印技术及应用[M].北京:科学出版社,2004:510-512.

小波包算法篇2

近年来,数字水印技术作为一种新兴的数字产品版权保护技术,成为网络信息安全研究领域的一个热点。基于小波域的数字水印技术因具有良好的多分辨表示、时频局部分析等特性,且易于与JPEG2000、MPEG4压缩标准兼容等特点,得到了广泛的关注。但是,小波变换域数字水印算法也存在缺陷,主要是在处理一些纹理图像时存在透明性与鲁棒性不甚理想等问题。其根本原因在于传统的塔式小波仅仅递归分解了低频带,未能对纹理细节丰富的高频子带实施分解处理。

本文提出一种基于小波包变换和视觉模型的数字水印算法。该算法采用小波包变换对原始图像进行分解,因为小波包变换比小波变换更符合信号的频域特征,而且只有知道图像的小波包分解结构才能提取出水印,因而能够提高水印的安全性。水印图像经置乱后被分块嵌入原始图像中不同的中频子带,并结合人类视觉特性确定不同区域的嵌入强度,以实现水印的自适应嵌入,采用双极性量化的方法修改选定的小波包系数,使算法具有良好的透明性和鲁棒性。

2 图像的小波包分解

图像的小波包分解可以用一个四叉小波包树来表示,一个三级的小波包分解树如图1所示。在对原始图像进行分解时,是按小波包分解算法逐层进行分解,即不仅对低频段进行分解,而且对高频段也进行分解,可以提高频率分辨率,因此它不仅是一种比小波多分辨分析更加精细的分解方法,而且具有更好的时频特性。图像的三级小波包分解如图2所示。

为了能够自适应地嵌入水印强度,本文采用Lewis[1]和Barni[2]提出基于小波域的视觉模型,利用从视觉模型导出的JND[3]来确定在图像的各个部分所能容忍的数字水印信号的最大强度,从而避免破坏视觉质量。根据视觉掩盖特性,水印的不可感知性受信号频率、背景亮度、纹理复杂度的影响。在该模型中,DWT系数的掩蔽特性qθl(i,j)的计算公式为:

其中,θ和l分别表示小波子带的方向和小波分解的级数,Θ(l,θ)为各级子带的频率掩蔽特性,Λ(l,i,j)为各级子带的亮度掩蔽特性,Ξ(l,i,j)为各级子带的纹理掩蔽特性。

3 水印的嵌入

水印嵌入算法的基本思想是:用推广的Arnold变换对水印图像进行置乱,以提高算法的安全性。对原始图像进行三级小波包分解,选择四个中频子带分别嵌入水印图像的1/4。将选定的子带进行分块,每块嵌入一位水印信息,嵌入时采用改进的双极性量化的方法修改小波系数。本文在量化方法上进行了改进,先按照视觉模型计算图像局部区域的JND阈值,再结合区域的系数特征计算嵌入强度α,根据α确定量化步长△的大小,实现了水印的自适应嵌入,有利于提高水印的不可见性。

设原始图像是大小为M×M的灰度图,水印图像是大小为N×N的二值图像,水印的嵌入过程按如下几个步骤进行:

3.1 图像置乱

对水印图像采用Arnold变换进行置乱,水印密钥为P和T。

3.2 小波包分解

对原始图像做三级小波包分解,以得到一个最低频子带和若干高频子带。传统的金字塔式小波分解,图像仅在低频部分不断地被分割,而小波包算法可对任意的子带进行进一步分解,因此提取水印时必须知道子带分解结构,提高了水印的安全性。

3.3 划分小波子块并确定嵌入强度

选择小波树中的四个中频子带作为嵌入区域。设所选的中频子带的大小为M×M,水印图像是大小为N×N的二值图像。将水印图像分为四部分,分别嵌入一个中频子带中。先将选中的中频子带分成N2/4块,每块大小为(2M/N)×(2M/N),并为其设定一个唯一的编号C,(C=0,1,2,…,N2/2-1),每个子块嵌入一位水印信息。对子块中的系数进行Z字扫描,形成一个一维序列,选择其中第K个系数嵌入水印,K=Cmod(4M2/N2),设选中的系数在子带中的坐标为(k1,k2)。对于处于Iθ1,θ2,θ3子带的编号为Cs的子块的嵌入强度α,其计算公式如下:

其中α0为预先设定的全局嵌入强度,qθ1,θ2,θ3(k1,k2)为该子块的JND阈值,mean(|Cs|)表示该子块内所有小波系数绝对值的平均值。

3.4 修改小波系数

对每一块子图作如下处理:设选中的小波系数为X,嵌入水印信息后X被修改为X',修改小波系数的过程如下:

1)划分区间集:选取△将坐标轴划分成如图3所示的A区间集和B区间集。区间宽度△的大小由下式决定:

△=α·△0

其中α为嵌入强度,△0为预先选定的全局量化步长。

2)修改小波系数。根据待嵌入的水印比特W和选定的小波系数X所在的区间集对X进行修改,修改的结果X'可表示为:

根据选定的区间宽度△对X进行取整数商,设求得的整数商m,则有:

m=[X/△]

具体的修改方法如下:

这样,对小波系数X修改之后,水印比特W包含的信息由修改后的小波系数X'所在的区间集唯一确定:如果X'处在A区间集内,则X'代表水印比特信息“1”;反之X'处在B区间集内,则X'代表水印比特信息“0”。区间集的宽度△由嵌入强度α调整,α越大则△也越大,水印的鲁棒性也就越好。

3.5 用修改后的系数作小波包逆变换,得到嵌入水印的图像

由于将水印图像分块嵌入不同的小波包子带中,因此大大提高了算法的水印容量。另外,提取水印时必须知道图像的小波包子带分解结构,从而增加了水印的安全性。

4 水印的检测和提取

提取水印时先将图像做3级小波包变换,选择嵌入水印的中频子带,对选择的中频子带分块并编号,然后根据各子块的小波系数平均值及该子块的JND阈值确定该区域的嵌入强度α,再由α计算出相应的量化步长△,利用△针对每个分块划分区间集A和B。再由分块的编号和大小找出嵌入水印的小波系数X',则提取出的水印比特信息W’由下式决定:

也就是说,如果X'位于A区间集,则水印比特为“1”,若X'位于B区间集,则水印比特为“0”。这样,在每个子带中分别恢复出一幅水印子图,最后将四幅子图拼接,根据置乱密钥对其进行反置乱变换得到水印图像。

5 实验结果分析

5.1 性能测试

为了验证本文算法的有效性,本文在Matlab下进行数字水印的嵌入和提取实验。原始图像为512×512的标准灰度图像,如图4(a)所示,水印图像为64×64的二值图像,如图4(b)所示,采用本文算法把水印嵌入原始图像后,得到嵌入水印后的图像,如图5(a)所示,再用提取算法从含水印图像中提取出二值水印图像,如图5(b)所示。

原始图像和嵌入水印后图像的大小、亮度、灰度范围、熵等性能参数如表1所示,从表1可以看出,嵌入水印前和嵌入水印后的图像在大小上没有任何变化,图像的亮度值以及灰度分布范围也没有明显的变化,图像亮度均值相对误差为0.1%,而灰度相对误差为5.2%,熵值变化相对误差为3.6%。

除了进行主观视觉效果评价外,还要进行水印嵌入与提取的定量评价,水印的隐蔽性评价常采用信号处理中的峰值信噪比PSNR,它是在原始图像X和含水印图像X'间进行的,其定义为:

单位是dB,PSNR值越大,水印隐蔽性越好。水印的鲁棒性评价在原始水印W与提取水印W'间进行。通常采用归一化相关系数ρ对水印和提取水印的相似性进行定量评价,它定义为:

相关系数值在0~1之间,相关系数越大,水印鲁棒性越好[4]。

在本次实验中,峰值信噪比PSNR=48.37,相关系数ρ=0.99。

5.2 抗攻击能力测试

为了检测水印的鲁棒性,对含水印图像进行一些攻击。常见的攻击方式包括叠加噪声、几何剪切、JPEG压缩、平滑滤波、旋转复位、图像增强等。表2给出含水印图像受攻击后提取出的水印的相关系数值。

由表2可见,除几何剪切攻击外,其他攻击发生时,相关系数均大于0.9。而且,在45%的情况下,相关系数均不小于0.95。这说明,本文的数字水印算法具有理想的鲁棒性。

6 结束语

本文提出了一种基于小波包变换和人类视觉特性的自适应盲水印算法。算法利用更符合信号频域特性的小波包变换对原始图像进行分解,将水印图像置乱后分块嵌入不同的中频子带中,并结合人类视觉特性自适应地确定每个区域的嵌入强度。本文创新之处在于采用改进的双极性量化的方法修改选定的小波包系数,提高了水印的不可见性。实验结果表明,本文提出的数字水印算法不仅具有较好的透明性,而且可以有效地抵抗常见的图像处理手段,具有理想的水印鲁棒性。

参考文献

[1]Lewis A S and Knowles G.Image compression using the 2-D wavelet transform[J].IEEE Transactions on Image Processing,1992,1(4):244-250.

[2]Barni M,Bartolini F and Piva A.Improved Wavelet-based Watermarking Through Pixel-wise Masking[J].IEEE Transaction on ImageProcessing,2001,10(5):783-791.

[3]X.M.Niu,S.H.Sun.Wavelet-based Gray-level Digital Watermarking[J].Journal of Harbin Institute of Technology,2000,8(2):111-115.

[4]黄达人,刘九芬,黄继武.小波变换域图像水印嵌入对策和算法[J].软件学报,2002,13(7):1290-1295.

[5]王向阳,杨红颖,陈利科.基于人眼视觉系统的自适应量化数字水印算法研究[J].小型微型计算机系统,2005,26(9):1525-1529.

[6]陈中孝,王曦.基于最优小波包基的自适应数字水印算法[J].微电子学与计算机,2005,22(2):114-116.

小波包算法篇3

关键词：遥感图像,小波包,父冲突,扩展方向树,SPIHT

1 引言

遥感,是20世纪60年代兴起并迅速发展起来的一种信息高技术,已在城市规划、环境保护、军事侦察等领域得到了广泛的应用[1]。遥感图像的特征是数据量庞大,图像尺寸高达30000×30000像素,谱段包含数百个,因此,要保证遥感数据能够在有限带宽的信道中传输,就必须进行压缩。

与常规图像相比,遥感图像的纹理更多、分辨率更高,空间相关性很小,采用普通的压缩方法可能造成纹理、边缘的损失,因此,如何在提高压缩比的同时又能较好地恢复细节也就成为了遥感压缩的重点。基于此,我们利用小波包优良的高频分析能力[2],提出了一种结合SPIHT的小波包编码算法—扩展方向树小波包压缩算法(Extended SOT Wavelet Packet Coding)。该方法借鉴SPIHT算法中的零树结构[3,4],通过重新定义其方向树,即扩展方向树(Extended SOT),解决零树结构中的“父冲突”(parental conflict)[7,8,9]问题,最终通过SPIHT算法实现整个编解码。

2 小波包变换

小波包变换[2]是小波变换的推广,它通过对小波变换的高频部分作进一步分解,使信号频带划分更加精细,并且能够自适应地选择小波包基,以达到对原始信号的最优逼近。小波包定义如下所示:

称由式(1)所定义的函数族{un(n=0,1,2···,n∈Z+)}为由φ生成的小波包。其中u0(t)=φ(t),u1(t)=ψ(t),φ(t)和ψ(t)分别为尺度函数和小波函数。

同时,将un进行二进伸缩和平移可以得到:

定义由式(2)组成的函数库{uj,n,k(t),j,k∈Z,n∈Z+}为由φ(t)导出的小波包库。由式(2)不难看出,小波包库由许多小波包组成,而不同小波包又具有不同性质,反映不同的信号特性,小波包必然具有分解形式多样化的特点。因此,本文采用Shannon熵[5,9]作为代价函数来得到小波包变换的最优基。

3 扩展方向树

传统小波编码如SPIHT算法,成功地利用了小波变换的零树结构和子带间相关性,通过方向树有效地表示了有效值映射,其方向树定义如图1所示。

图中每个子带用黑色正方形表示,其边长与尺度成反比。除最低频子带LLM和最高频的三个子带外,每个系数在它的同方向相邻更高频子带的相同位置上有四个子女,最高频子带没有子女,同方向相邻子带间尺度满足二进关系。

然而,由于小波包有可能对高频子带作进一步分解,相邻子带间差异可能不再满足二进关系,其父子关系会更复杂,子系数可能出现在比父系数尺度更粗糙的子带中。这样就出现了如下两类“父冲突”:

1)一个子系数对应多个父系数。2)不同节点的父系数对应同一个子带的子系数。

以三层小波包分解为例(图2(a)所示),R表示最低频子带LL3,T1,T2,T3表示水平,垂直,对角三个不同方向树的根节点,每个方向在下一精细尺度的子带都是其子节点。这样,可以定义如下四种父子关系:

a)子节点位于父节点同方向相同位置的更粗糙尺度子带。如图2(b),树T2:14与15、16、17、18所示。

b)子节点位于父节点同方向相同位置的相同尺度子带。如图2(b),树T2:14与19、20、21所示。

c)子节点位于父节点同方向相同位置的下一个精细尺度子带。如图2(b),树T3:7与11所示。

d)子节点位于父节点同方向相同位置的更精细尺度子带。如图2(b),树T1:1、2、3、4与5所示。

可以看出,A、D两种父子关系与小波变换父子关系不同,产生了“父冲突”问题,我们称A为第一类“父冲突”,D为第二类“父冲突”。对于前者,文献[8]通过将产生冲突的粗糙尺度的子节点上移,使之与更粗糙或者相同尺度下的节点产生父子关系而解决冲突。而对于第二类“父冲突”,由于四个父节点都对应于同一个子节点,我们将其拆分,产生冲突的子节点只与其中一个父节点产生父子关系。综合考虑C、D两种父子关系,一个父节点就不是仅仅对应一个子节点,一个父系数也不是仅仅对应4个子系数,而是应该对应4N,N表示父子之间的尺度差异。因此本文采用的方向树不同于SPIHT算法,我们称为扩展方向树,在实际编码时需对子节点个数和子孙数目做相应调整。图3所示解决“父冲突”后的方向树图。

4 零树假设验证

对于一个给定的小波包分解,其方向树的定义如上所述。而零树结构是建立在一个假设的基础上的:在小波分解系数中,对于给定的阈值,如果父系数是不重要的,那么很有可能子孙系数也是不重要的[3,4]。正是利用这种假设,才能够有效表示“无效值映射”,也就等价于有效地表示了有效值映射,这就是零树结构的基础。因此,本文编码算法是否成功,首先就取决于零树假设在小波包分解中是否成立。我们引入小波分解零树结构检验的经验统计方法,该方法同样适用于小波包分解。

我们采取两种方式统计变换系数的幅值:1)按照子带所属频带,以频率递增的方向画出所有系数幅值;2)对不同的方向树,从根节点开始,按照先父节点后子节点的顺序分别画出水平,垂直以及对角方向的系数幅值。以一幅2048×2048的遥感图像为例,小波及小波包分解后子带扫描顺序分别如图4(a)、图4(b)所示,扫描按照子带标号递增方向进行,子带内部则以行为单位,按照从上到下的顺序扫描。图5则表示在当前小波包分解下所生成的扩展方向树。

图6和图7分别表示小波分解和小波包分解后按频率增加顺序和方向树顺序画出的幅度值。在图6(b)和图7(b)中,小波分解和小波包分解都能够按照方向树将系数分为三个部分,而且系数基本呈现由大到小逐渐变化的趋势。对于小波包分解,在每一个树内,还存在若干子树(图8),每一个子树内部同样呈现父节点系数幅度较大,子节点系数幅度较小的趋势。换言之,对同一个树,父系数总是比子系数大,如果父系数不重要,子系数也很有可能不重要,因此,小波包扩展方向树满足零树假设。

5 实验

为验证ESOT-WPC算法的准确性,我们选择3幅2048×2048×8bit的遥感图像进行实验,并与SPIHT算法比较,三幅图像分别命名为RS1、RS2和RS3。小波包分解选用双正交Db9/7小波,5层小波分解,并通过代价函数的计算确定最优基,解决方向树中的“父冲突”问题,最优基结构如图9所示。SPIHT算法同样选用双正交Db9/7小波,5层小波分解。

图10、图11和图12分别是RS1、RS2和RS3在0.5bpp压缩率下,使用SPIHT算法和ESOT-WPC算法压缩后得到的图像。不难看出,图10(a)与图10(b)和图10(c)在视觉观察上基本无异,这说明如果图像冗余度较高,包含大量平坦区域,本文算法和S PIHT算法都能得到较好的效果;对于RS2和RS3,图11(a)中边缘在图11(b)中虽有所恢复,但已经发生扭曲,而在图11(c)中则基本保持了边缘的线性。图12(a)白色柱状物之间的斜向纹理在图12(c)中基本得到恢复,但在图12(b)中则是模糊一片。表1是三幅图像使用SPIHT算法和ESOT-WPC算法压缩后的PSNR(峰值信噪比)的对比,对于细节变化很少的RS1图像,由于图像中包含大量的低频信息,过多的高频分解反而会降低图像的压缩效率,因此本文算法的P SNR低于SPIHT。而对于富含纹理和细节的图像RS2和RS3,本文算法在无论在PSNR还是在视觉效果上都明显优于SPIHT算法,更好的恢复了遥感图像的细节信息。

6 结论

本文针对普通压缩算法容易造成遥感图像高频信息丢失的问题,提出一种结合SPIHT的小波包编码算法—ESOT-WPC。借鉴SPIHT算法中的零树结构,通过重新定义子带之间的父子关系,即扩展方向树,解决其“父冲突”问题,并成功实现编解码。实验表明,该算法在PSNR和视觉效果上均优于SPIHT算法,是一种有效的遥感图像压缩方法。但目前算法计算复杂度较高,需进一步改善算法性能,降低复杂度。

参考文献

[1]彭望琭.遥感概论[M].北京:高等教育出版社,2002.PENG Wang-lu.The Overview of remote sensing[M].Beijing:Higher Education Press,2002.

[2]Ramchandran Kannan,Vetterli Martin.Best wavelet packet bases in a rate distortion sense[J].IEEE Trans.Image Processing,1993,2(4):160-175.

[3]Said Amir,Pearlman William A..A new fast and efficient image codec based on set portioning in hierarchical trees[J].IEEE Trans.on CASVT,1996,6(3):243-250.

[4]Shapiro Jerome M.Embedded image coding using zero-trees of wavelet coefficients[J].IEEE Transactions on Signal Process-ing,1993,41(12):3445-3462.

[5]Coifman Ronald R,Wickerhauser Mladen Victor.Entropy based algorithm for best basis selection[J].IEEE Trans.Information Theory,1992,38(2):713-718.

[6]Xiong Zixiang,Ramchandran Kannan,Orchard Michael I.Wavelet packet image coding using space-frequency quantization[J].IEEE Trans.Image Processing,1998,7(6):892-898.

[7]Xiong Zixiang,Ramchandran Kannan,Orchard Michael I d.Space-frequency quantization for wavelet image coding[J].IEEE Trans.Image Processing,1997,6(5):677-693.

[8]Rajpoot Nasir M,Wilson Roland G,Meyer Francois G,et al.Adaptive wavelet packet basis selection for zerotree image coding[J].IEEE Trans.Image Processing,2003,12(12):1460-1472.

小波包算法篇4

电缆故障双端行波测距的准确度很大程度上取决于行波到达的时间和行波传播的速度。当电缆发生单相接地故障时将经历一个复杂的暂态过渡过程,暂态信号中含有丰富的高频分量。不同频率的行波信号具有不同的传播速度及衰减率,即产生色散,导致电缆故障行波波头信号在传播过程中发生畸变,影响了对行波准确到达时间的判别和对行波波速的确定,因此降低了电缆故障行波测距的准确度[1]。

近年来,小波分析方法在故障行波测距领域等方面得到广泛的应用[2,3,4,5]。但是行波传播色散特性导致的行波波头达到时间差以及波速的选择仍没得到很好的解决。文献[6]提出了采用相关分析法实现输电线路双端行波测距,通过相关分析方法确定两端波形到达时间差。但是双端行波信号所含频率成分不同,各频率成分的衰减也不一致,波形差异较大,难以得到较高的相关系数。

本文将小波包提取算法和相关分析方法相结合,利用小波包算法提取电缆双端较窄频段的信号来重构波形,对重构波形进行相关分析,进一步增大了相关系数,从而得到更加精确的行波到达时间差,行波波速由该频段频率特性决定,可近似为常数,解决了难以选择波速的问题,进一步提高了故障测距的准确度。

1 测距算法原理

1.1 小波包提取算法

小波包分析能将频带进行多层次划分,为信号提供一种更为精细的分析方法。设hk和gk是相互正交的镜像滤波器,则小波包定义为由正交尺度函数Φ(t)确定的函数族。

其中:W0(t)=Φ(t),W1(t)=φ(t),即W0(t)为尺度函数Φ(t),W1(t)为小波函数φ(t);n为采样点总数,k为采样点序号,k=1,2,…,n。

小波包提取算法就是利用小波包原理将离散信号按不同频段正交分解,在不同的尺度下分解出不同频段的信号,保留分解序列中感兴趣的频段,舍去其他频段,然后进行重构。重构后的信号频段比原始信号要窄,信号长度与原始信号相等。小波包提取算法的步骤如下[7]。

1)选取正交滤波器hk,令gk=(-1)k-1h1-k。

2)确定分解层数M,M>0。设置采样频率fs,第M层每个序列的带宽为fs/2M+1。

3)进行逐层小波包分解,提取有效频段,剔除无效频段。

4)利用小波包的重构公式对有效频段进行重构,重构后的信号与原始信号长度相等,频段宽度为fs/2M+1。

1.2 小波包提取信号的相关分析

故障点处产生的信号可以等效为各个频段信号的叠加,在向电缆线路两侧线路传播过程中,不同频段的信号会产生不同程度的衰减和相移,但各个频段内的波形不会发生大的变化,因此,双端信号小波包提取后具有很大的相似性。可采用相关分析法来求取到达检测点的行波波头之间的时间差。

相关分析是研究随机变量之间相关性的统计分析方法。确定性信号可以看作是平稳的,且具有遍历性的随机信号的特例,因此可以用相关分析方法来判断两个波形的相似程度[8]。设x(t),y(t)为两个能量有限的连续函数,则它们相关运算定义为

式中,Rxy为x(t),y(t-τ)在时延τ上的互相关函数,也称相关系数。

相关运算的过程相当于是将一个函数平移后,对另一个函数进行扫描,反映的是两个函数在不同偏移量下的相似程度。求τ使相关系数Rxy最大,此时的波形x(t)平移τ后恰好与y(t)对齐,也即求得了两个波形的时间差。

1.3 有效频段的选取

信号进行小波包分解后,有些频段信号微弱,并不能代表原信号,对其进行相关分析,可能出现伪最大相关系数。为此,须选取能量相对集中的频段进行小波包重构。方法如下:首先,信号经小波包分解后,得到一系列频段,分别记为:C1,C2,C3,…,Cn,代表从低频到高频的信号。然后求各个频段的总能量,设Cj(j=1,2,…,n)对应的能量Ej(j=1,2,…,n),则有

式中,xjk(j=1,2,…,n;k=1,2,…,m)为Cj离散点的值。最后选取能量相对集中的频段进行小波包重构。

同时,为了使电缆双端信号都能代表原信号,以远离故障端为准选取频谱中能量集中的频段,进行小波包提取,近故障端根据相应频段进行提取。若以近故障端为准,则远故障端信号在该频段可能衰减幅度较大而不能代表原信号。

2 测距方法的实现

基于小波包提取算法和相关分析的电缆双端行波测距方法步骤如下。

1)求取线模分量。对故障后电缆双端的三相电流信号进行采样,相模变换,选用公认的色散较小的行波线模分量进行测距。暂态行波信号是叠加在工频信号上的,因此在对暂态行波信号进行分析之前应采用高通或带通滤波器滤掉工频信号和各次谐波,保留暂态行波信号。

2)小波包提取。确定远故障端和近故障端。对双端电流行波线模分量傅里叶变换,得到能量频谱,确定远故障端和近故障端。根据行波传播理论,行波传播过程中,由于高频分量衰减快,故障行波中的有效频率范围反比于故障距离[9],容易确定远故障端和近故障端。以远离故障端为准选取频谱中能量集中的频段,进行小波包提取,近故障端根据相应频段进行提取。

3)相关分析。利用式(3)对提取的信号进行相关分析,得到最大相关系数,确定行波到达送端和受端时间差Δt。

4)测距。根据双端测距原理,可得到故障距离为

式中:x为故障点距送电端的距离;v为选定频段的波速;L为电缆全长。

3 实验及数据分析

3.1 实验模型

为了验证方法的正确性,本文对实际电缆进行了实验,实验电路如图1所示。

由于三相线路故障时,可以将故障量分解为α,β,0这3个独立的量,同样可以按单相线路的方法计算,所以实验采用单相线路。

电源为单相10 V工频电压,电缆长度为323.9m,故障点F接地电阻为5Ω,距离M点61 m。P1,P2为同规格的高频电流传感器,分别连接至示波器的通道1和通道2,以达到时间同步。

3.2 实验数据分析

F点发生接地故障时示波器采集到的传感器P1和P2的波形如图2所示,其中纵轴U为传感器输出电压,横轴n为示波器采样点数,示波器的采样频率设置为100 MHz。取P1和P2采集到的第一个电流行波进行分析。

根据实测数据,传感器获得的电流行波频率主要分布在0~10 MHz,确定上限频率fs=10 MHz,将电流行波进行3层小波包分解,则整个频带被分成8段,每个频带宽度为fs/2M=1.25 MHz,得到能量频谱图如图3所示。从图2中容易确定传感器P2端为远故障端,传感器P1端为近故障端。

然后提取P2能量最集中的C1频段(0~1.25MHz)进行重构,P1端也以0~1.25MHz频段作为有效频段进行重构,双端重构波形如图4所示。

对图4波形进行相关分析,得到最大相关系数时延,即行波到达M和N端的时间差Δt=1.18µs。取故障行波的波速作为C1频段的行波速度,故障行波速度由实测方法确定,采集电缆尾端发生故障时的双端传感器数据,测得故障行波速度约为v=1.73×108m/s。利用式(5)计算得x=59.88 m,与实际故障距离相差1.12 m。

对电缆在不同故障距离、不同接地电阻的电缆故障进行实验,利用小波模极大值法和小波包相关分析法分别对实验结果进行分析,结果如表1和表2所示。由数据可表明所得测距结果有较高的精确性,满足故障测距的要求。

4 结论

根据小波包分解和重构原理,提取故障行波中能量相对集中的频段进行相关分析,为电缆双端行波故障测距提供了一种新的方法。该方法能准确地求取波头到达时间差和波速,有效地消除了行波传播色散对行波测距准确度的影响。通过建立实验模型,对实验数据进行分析,验证了本文提出方法的有效性,提高了电缆故障测距精度。

参考文献

[1]曾祥君,张小丽,马洪江,等.基于小波包能量谱的电网故障行波定位方法[J].高电压技术,2008,34(11): 2311-2316. ZENG Xiang-jun,ZHANG Xiao-li,MA Hong-jiang, et al.Traveling wave fault location method for power grids based on wavelet packet energy spectra[J].High Voltage Engineering,2008,34(11):2311-2316.

[2]马丹丹,王晓茹.基于小波模极大值的单端行波故障测距[J].电力系统保护与控制,2009,37(3):55-59. MA Dan-dan,WANG Xiao-ru.Single terminal methods of traveling wave fault location based on wavelet modulus maxima[J].Power System Protection and Control,2009,37(3):55-59.

[3]成乐祥,李扬,唐瑜.基于改进形态Haar小波在输电线路故障测距中的应用研究[J].电力系统保护与控制,2010,38(6):30-34. CHENG Le-xiang,LI Yang,TANG Yu.Research and application of improved morphological Haar wavelet in transmission line fault location[J].Power System Protection and Control,2010,38(6):30-34.

[4]何正友,钱清泉.电力系统暂态信号分析中小波基的选择原则[J].电力系统自动化,2003,27(10):45-49. HE Zheng-you,QIAN Qing-quan.Mother wavelet option method in the transient signal analysis of electric power systems[J].Automation of Electric Power Systems, 2003,27(10):45-49.

[5]周鑫,吴飞,张默霓,等.基于小波变换的T型线路故障测距[J].电力系统保护与控制,2010,38(2):8-11. ZHOU Xin,WU Fei,ZHANG Mo-ni,et al.A new fault location method for T-connection transmission lines based on wavelet transform[J].Power System Protection and Control,2010,38(2):8-11.

[6]DU Lin,PANG Jun,SIMAWen-xia.Correlation analysis algorithm for transmission line fault location based on travelling wave[J].High Voltage Engineering,2008,34 (21):2637-2641.

[7]马强,荣命哲,贾申利.基于振动信号小波包提取和短时能量分析的高压断路器合闸同期性的研究[J].中国电机工程学报,2005,25(13):149-154. MA Qiang,RONG Ming-ze,JIA Shen-li.Study of switching synchronization of high voltage breakers based on the wavelet packets extraction algorithm and short time analysis method[J].Proceedings of the CSEE, 2005,25(13):149-154.

[8]王清亮,刘军良.基于高频暂态分量相关性的选择性漏电保护[J].电力自动化设备,2007,27(9):59-61. WANG Qing-liang,LIU Jun-liang.Selective leakage protection based on correlation between high-frequency transient components[J].Electric Power Automation Equipment,2007,27(9):59-61.

小波包算法篇5

基于小波变换的小波包变换,不仅在时频域都具有良好的局部化性质,而且对信号的高频部分提供了更精细的分解,从而得到大量高频信息[1,2,3]。小波包变换在电力系统多方面的应用都有明显的优势,例如在数据压缩领域,小波包变换可以针对电力系统故障录波信号中蕴含丰富暂态特征信号的特点,将故障信号分解成各个频带的信息,选择有重要信息的频带采用各种压缩算法进行压缩。近些年,在小波包算法基础之上,又提出了与嵌入式零数编码相结合[4]、与快速傅里叶变换相结合[5]、混合小波包[6]等各种提高压缩精度的压缩算法。这些算法虽然在一定程度上使压缩性能得到提升,但却增加了算法复杂度。随着电力系统数据记录装置的飞速发展、电力系统高采样的逐步演变、滤波装置的不断更新,电力系统电能质量数据逐渐呈现海量化趋势,传输、存储等环节对实时性的要求又进一步提高,因此,电力系统数据压缩效率问题得到越来越多的关注。以上各种压缩算法所基于的小波包变换,其核心部分的卷积运算过程复杂、运算量大、实时性差,严重阻碍了其在实际工程中的应用。为了提高小波包变换的计算速度,解决目前电力系统海量数据的压缩效率问题,本文提出将小波包变换并行化处理。

目前,针对这一问题,国内外学者提出利用各种并行技术来加快小波包变换的运算速度[7,8,9,10,11,12,13],大部分文献都是利用多机网络环境进行并行化处理,该方法不可避免地受到系统维护及网络传输环节所带来的不利影响。随着计算机的发展,单机多核处理器应用越来越广泛,因此基于多核技术的并行小波包研究具有十分重要的意义。

本文利用Pthreads[14]和Open MP[15,16,17]2种并行编程环境,针对提高小波运算速度问题,提出小波包Mallat分解重构算法不同的并行策略。在Pthreads环境下,提出小波包分解过程数据级并行、重构过程任务级并行;在Open MP环境下,分别提出小波包分解重构特征嵌套和非嵌套并行策略。并利用C语言在单机双核处理器平台上进行了海量数据压缩实验研究,取得了明显的效果,为解决电力系统数据压缩效率问题提供一种新的方法。

1 小波包Mallat算法并行性分析

小波包是小波概念的推广。它是在小波变换(V0=V1⊕W1=V2⊕W2⊕W1=…)只将V(尺度)空间进行分解的基础上,进一步对Wj(小波空间)进行分解,使正交小波变换中随j的增大而变宽的频谱窗口进一步变细,达到最适合于待分析信号的时频窗口或最优基[18]。其分解过程如图1所示,图中h0与h1分别表示小波包分解后的高频系数与低频系数,a表示尺度空间数据,d表示小波空间数据。

小波包分解重构过程的核心是滤波器与待处理数据的卷积过程。从卷积计算的公式可以了解到,卷积计算结果与待处理数据序列直接相关,不同位置的待处理数据与滤波器的卷积过程得到不同的卷积结果。因此不同位置的卷积结果之间并没有相关性,可以将待处理数据进行分配并行卷积。

小波包Mallat重构原理如图2所示,图中g表示由小波空间数据与尺度空间数据重构后得到的上一层的重构数据。可以看出,小波包变换中对尺度空间和小波空间的处理也不存在关联性,该部分同样可以作为并行任务执行。

2 基于Pthreads与Open MP的并行小波包实现

2.1 Pthreads与Open MP

Pthreads是Linux操作系统的多线程接口标准,在Windows操作系统上,可以利用Pthreads-win32开放源代码版本实现并行编程。与Pthreads库不同,Open MP是通过采用指导语句、库函数调用和环境变量来给用户提供在共享存储计算机上创建和管理可移植并行程序的方法。它是一个在多处理机上用以编写可移植的多线程应用程序的API库,采用ForkJoin并行模式。

2.2 基于Pthreads的并行小波包实现

小波包对信号数据处理的流程,需要经过数据加载、延拓、分解、二抽取、二插值、重构这一过程。利用Pthreads实现分解过程数据级并行分解、重构过程任务级重构。其中,延拓、二抽取与二插值的过程与分解、重构过程融合为同一个流程。

a.数据延拓。为了执行数据延拓,将每组与后一组重复NFilter-1个数据,其中NFilter是滤波器系数个数,本文采取周期延拓,最后一组与第1组数据重复前NFilter-1个数据,假设有p个线程参与执行并行计算,则每组有N/p+NFilter-1个数,N为原始数据个数,卷积时滤波器第1位只对应每组前N/p个数进行计算。

b.数据级并行分解。由于小波包分解时各数据序列之间没有相关性,可以将各个数据分配给各个线程同时执行。为保证各个执行线程负载平衡,本文根据小波包分解层数将近似系数与细节系数的数据分别分为N/p组,每层N的大小为上一层的1/2,数据分组如图3所示。

将分解的前q组数据分给q个线程,每个线程同时执行分解过程,待线程挂起后,主线程继续向这些子线程派发新的q组数据。分解一层后,先对近似系数进行分组并行分解,再对细节系数分组执行。此时需要注意的是,由于延拓方式的影响,对应于每一层中近似与细节系数的延拓值将会不同,因此在分组完成后先要计算此次延拓后的余值项,再采用新的余值进行延拓与卷积计算。

c.二抽取。为了减少小波包变换计算时间,可以将分解中二抽取过程与分解过程相结合[20],即不再计算所有数据,而是选择奇数位数据计算,偶数位数据直接跳过。图4采用滤波器系数为4说明此过程。图中当数据延拓分组后,不同的线程针对本组数据进行卷积计算时,每个线程执行第1、3、5、…位为首的运算。

图5为小波包Pthreads并行分解流程,图中Ti为不同的线程,di为不同组的数据,主线程完成数据载入与延拓后,将数据分组,分出一组数据后创建一个线程执行该组的分解运算,完成近似系数的分解,接着继续分解第i层细节系数,直到所有系数分解完成。

由于数据回存需用到上一步线程的结果,并且是在同一堆栈中进行读写,为了保证数据回存的正确性,各线程完成卷积计算以后,需在分解数据回存之前与主线程合并,主线程完成分组后在此等待,创建的线程得到释放。

d.二插值与任务级重构。重构的过程中可以发现,每一层对应的近似系数和细节系数个数是相同的,并且不同的系数与不同的滤波器进行运算,不同类型系数之间没有相关性,读取数据实现卷积的过程就不会出现同步的问题,因此可以将近似系数和细节系数重构的过程看成是不同的任务,实现任务分解方式及任务大小相同,这也满足负载平衡要求。

小波包重构并行过程,先根据每一层循环与分解层数的关系,计算数据延拓后的系数值,二插值时不进行补零操作,而是在非补零位置计算待处理数据与滤波器奇数位的系数卷积,补零位置计算待处理数据与滤波器偶数位的系数卷积。而整个过程中尺度空间与小波空间同时进行重构运算,其任务级并行重构流程如图6所示。

2.3 基于Open MP的并行小波包实现

Open MP标准可以直接在串行程序的基础上调整程序的并行性,因此利用Open MP实现小波包并行必须首先编写小波包分解的串行程序,本文在Visual studio 2008平台上利用C语言完成串行编程。Open MP最显著的用处在于,可以将小波包for循环分配到不同处理器上同时执行,即根据处理器的个数p将循环次数分为p组,每组循环由一个线程执行,这样由p个线程同时完成不同次循环体中卷积、二抽取计算,如图7所示。

小波包分解串行程序核心部分是一个嵌套循环过程。嵌套循环的含义是把小波包分解层次与计算低频、高频系数的次数相联系。为处理边界问题,对边界进行周期延拓,共享数据作为私有副本分配给各个线程。根据图2,假设分解n层(n≥1),外层循环需要进行2n-1次。假设外层分解至第j(j=0,1,2,…,2n-2)次,内层循环则从j N/(2n-1)开始至(j+1)N/(2n-1)次结束。

a.嵌套方案。在Open MP内,嵌套是一个默认不使能的布尔属性。当并行区域嵌套出现在线程已经运行在并行区域又遇到另一个并行区域时,若嵌套使能,那么程序就会按照之前动态线程的规则生产一个新的线程组。可以通过omp_set_nested(i)来设置嵌套使能状态,i=0表示嵌套未使能,i=1表示嵌套使能。

小波包分解过程利用嵌套使能并行执行嵌套循环。当主线程遇到外层循环时,派生出一组线程,这些线程同时继续向内层循环执行下一层并行,这时派生出的线程和原来的主线程成为针对内存循环的主线程,继续按照之前的方式派生新的线程从而并行执行内层循环,当内层循环执行结束时,这些线程逐级挂起,最终回到最初的主线程上。嵌套循环并行执行线程运行过程如图8所示,m为实际计算机核数。

b.非嵌套方案。另外一种方法即类似于小波分解并行执行过程。对于嵌套循环的程序,可以选择只对一层循环作并行处理,由于小波包分解中外层循环次数不多,而内层循环复杂度较高,本文选择只对内层循环过程作并行处理,即利用多核处理器同时对不同尺度的系数做卷积和二抽取过程。非嵌套循环并行执行线程运行状态如图9所示。

3 数据压缩分析实验

3.1 实验基本过程

为了检测并行小波包在数据压缩中的性能表现,本文对数据进行并行小波包压缩解压处理,由于小波包变换可以使信号经过小波包变换后能量集中在少数的系数上,并且这些系数的能量能够远大于其他多数小波包的系数,因此可采用阈值方法,现有的文献阈值选取包括全局阈值与局部阈值,考虑到全局阈值能降低计算的复杂度,而本文分析重点在于研究算法的并行性,因此选取全局阈值压缩。

数据压缩方法采用阈值压缩算法,阈值处理后需要记录非零小波系数的位置以及系数的取值,此过程采用文献[22]所提出的位图压缩法。

小波分解的最大级数J定义[23]如下:

其中,ff是信号的基础频率,fs是信号的采样频率。

实验主要步骤如下。

a.载入信号后,对信号作并行小波包分解;分解层数根据式(1)得出,本实验分解层数为3层,分解过程采用第2节提出的小波包并行分解方案。

b.将分解后的小波系数,经过全局阈值处理,利用位图压缩算法将剩余后的小波系数依次排列,完成压缩过程。解压过程根据位图压缩法还原各系数原始位置,对其余位置进行补零操作。

c.采用第2节提出的并行小波包重构方案,对解压后系数作小波包重构,还原原始数据。

实验主要用于检验并行小波包的计算效率,在保证分解重构的结果与串行计算结果相一致的前提下进行。

实验时间记录从主线程执行并行小波包分解开始,直到小波包重构结束。由于实验针对小波包分解压缩重构的并行过程,因此不记录载入信号与写出信号的时间。

仿真信号源采用的是:

3.2 实验环境

本文采用的实验平台为:处理器为Intel(R)Core(TM)2 Duo CPU T5750,内存为3 GB。所有数据为5次计算时间平均值。

3.3 实验方案

为比较并行算法性能,采用本文所提出基于Pthreads及Open MP 2种并行编程环境下的小波包并行算法与串行小波包算法分别进行数据压缩。实验在210、215、216、220、221这5种数据量下,采用db4小波包,将信号分解为3层。实验结果见表1、2。

为了显示并行算法相对串行算法的性能改变,采用加速比描述:加速比=串行算法计算时间/并行算法计算时间。

表3中S1、S2、S3分别表示Pthreads环境下小波包并行压缩解压加速比、Open MP环境下小波包嵌套并行压缩解压加速比、Open MP环境下小波包非嵌套并行压缩解压加速比。

3.4 实验结论

a.从实验可以看出基于Pthreads与Open MP的并行小波包算法,在计算速度上有着较大的优势,这种性能的提升随着数据量的增大而增大。相对于串行算法有着接近2倍的速度提升,甚至在用Open MP进行并行处理时出现超线性加速比的情况。

b.在数据量较少时,由于并行计算不可避免的开销要比并行效果带来的节约时间要大,并行计算的性能不如串行计算,在数据量为215时,各种方法的加速比表现较为一致,且都大于1,说明在此数据量附近,并行带来的节约只略大于并行开销,一旦计算机有其他程序的任务,就会严重影响此时加速比的表现。这一问题将随着海量数据处理的发展而逐渐被忽略。

c.随着数据量的增大,基于Open MP的并行小波包加速比比基于Pthreads的并行小波包要好,这是由于Pthreads是直接控制线程调度,可以由程序员自主控制程序的底层运行,一旦代码优化效果好,就可以达到较高的加速比,反之,由于优化存在难度,在代码复杂度较高的情况下,可能达不到Open MP的并行效果。

d.嵌套循环的并行性能比非嵌套循环较弱,这是由于嵌套使能开启后,内部循环仍旧按照之前的方式派生线程,由于外层循环根据处理器个数已经产生最大可并行执行的线程,因此,在内层循环中所派生出的新线程与其主线程并没有同时执行程序而是交替进行。此外创建和退出线程的过程又带来一定的时间消耗,使得性能下降。但是在处理器个数大于2时,通过动态分配线程,改变外层和内层线程分配方式,可以使嵌套并行循环的性能更优。

4 结语

小波包理论与图像小波包分解篇6

小波包变换是一种提取虹膜纹理特征的有效方法。对虹膜图像的纹理特征的处理, 主要集中在中高频能量部分, 可以对图像进行二级小波包分解, 并提取出第一级和第二级的中高频部分的能量分量作为特征向量。

1 小波包理论

小波包分析是一种精细的分析方法, 它将频带进行多层次划分, 对高频部分进一步分解, 并能够根据被分析信号的特征, 自动选择相应频带, 进行匹配, 提高了时频分辨率。

1.1 小波包的定义

通常, 在多分辨率的分析中, 利用不同的尺度因子j, 把Hilbert空间L2 (R) 分解为子空间Wj (j∈Z) 的正交和, 即其中, Wj为小波函数ψ (t) 的小波子空间。为达到提高频率分辨率的目的, 应对小波子空间Wj进行频率的细分。

先用一个子空间Ujn来统一表征尺度空间Vj和小波子空间Wj, 设:

Hilbert空间的正交分解Vj+1=Vj⊕Wj, 于是可以把Ujn的分解统一为以下形式:

设Ujn是函数Un (t) 的闭包空间, Un (t) 是函数U2n (t) 的闭包空间, 令Un (t) 满足双尺度方程:

其中, g (k) = (-1) kh (1-k) 。当n=0时, 上式变为:

在多分辨率分析中, φ (t) 和ψ (t) 满足双尺度方程:

将式 (4) 和式 (5) 相比较发现, u0 (t) 和u1 (t) 分别退化为尺度函数φ (t) 和小波基函数ψ (t) 。式 (4) 是式 (1) 的等价表示, 将它推广到n∈Z+ (非负整数) 的情况, 则式 (3) 可等价表示为:

小波包的定义:式 (3) 构造的序列{un (t) } (其中n∈Z+) 成为由基函数u0 (t) =φ (t) 确定的正交小波包。n=0, 即为式 (3) 的情况。

因为φ (t) 是由hk唯一确定, 所以又称{un (t) }n∈Z为关于序列{hk}的正交小波包。

1.2 小波包的空间分解

设{un (t) }n∈Z是hk的小波包族, 对式 (1) 作迭代分解, 生成子空间族如下 (n=1, 2, …;j=1, 2, …) :

于是, 有:

…

空间Wj分解的子空间序列记为Uj-12j+m, m=0, 1, …, 2j-1;l=1, 2, …。而U2j-1j+m的标准正交基是{2- (j-1) /2u2j+m (2j-ltk) :k∈Z}。当l=0和m=0时, U2j-1j+m为Ujl=Wj, 它的正交基也简化为2-jl2ul (2-jt-k) =2-jl2ψ (2-jt-k) , 即是标准正交小波族{ψj, k (t) }。

如果n是一个倍频呈细化的参数, 不妨令n=2l+m, 则小波包可以简记为ψj, k, n (t) =2-j/2ψn (2-jt-k) , 其中, ψj, k, n (t) 是有尺度指标j、位置指标k和频率指标n的小波包。而小波ψj, k (t) 只有参数离散尺度j和离散平移k, 小波包多一个频率参数n=2l+m。正因如此, 小波包克服了小波变换时频分辨率低的缺陷。参数n表示函数的零交叉数, 即其波形的震荡次数。

由尺度函数ψn (t) 生成的函数族ψj, k, n (t) (n∈Z+;j, k∈Z) 称为由ψ (t) 构造的小波库。

如果对于任意j=0, 1, 2, …, 有:

则对应的族{uj, k, un (t-k) |j=…, -1, 0;n=2, 3, …;k∈Z}是L2 (R) 空间的一个正交基。

尺度j增大, 相应正交小波基函数的空间分辨率越来越高, 而其频率分辨率越来越低, 这是正交小波基难以克服的缺陷。相比较之下, 小波包却有将随j增大而变宽的频谱窗口, 可进一步分割变细的优良性质, 克服了小波变换的不足。

2 图像的小波包分解

先用小波包将图像分解为四个子图像, 即为水平与垂直方向上的低频分量LL、水平方向上高频垂直方向上低频的分量HL、水平方向上低频垂直方向上高频的分量LH和水平与垂直方向上的高频分量HH。接下来可以仅对低频分量继续递归分解, 也可以包括对中频分量继续递归分解, 或进一步包括对高频分量继续递归分解。这三种分解情况所得到的子图像的数目是不相同的, 设分解级数为j, 这三种分解得到的子图像数分别为显然, 当分解的级数j比较大时, 三种子图像数目的差距也会很大。

而第三种分解方法, 也称为完全树结构分解, 就是小波包分解。从频谱空间的角度来看, 小波包分解指对下一级的U空间和各个V空间都继续分解。1-D小波包变换的三级频谱空间划分示意图见图1, 其中A表示近似, D表示细节。

3 结语

一般来说, 小波分析是将信号分解成低频的粗略部分和高频的细节部分, 然后仅对低频信息继续分解, 得到低频信息和高频信息, 依此下去, 而不对高频信息做分解处理。小波包分解则对各频带进行分解, 既对低频信息分解, 也对高频信息分解。用小波包分解分析纹理图像时, 对低频信息分析的同时, 对图像高频信息有选择地正交分解, 这样可以得到比较全面的图像特征信息。

参考文献

[1]Yang w L, Lu G M, Wang K Q.Iris recognition based on location of key points[J].Biometric Authentication, 2004

[2]邵宇, 马义德.基于小波变换与差分矩阵的虹膜识别技术[J].甘肃科技, 2007, (6) :27-29

小波包算法篇7

图像内容的主要特征包括形状、颜色、纹理等,纹理是其中重要的特征之一。目前纹理还没有统一的定义,纹理可以认为是灰度或颜色在空间以一定的形式变化而产生的图案[1]。在纹理分析中,最重要的是纹理特征的提取。纹理特征提取现有的方法基本上可归纳为统计法、模型法、空间-频率域分析法和结构法几类。由于纹理的多样性和复杂性,要实现对他们全面的描述是很困难的,这些分析方法都只能从某些方面反映出纹理的特征,而缺乏对纹理整体特征的描述。

基于上面的事实,本文提出一种小波包和分形相结合的纹理特征提取的算法。通过实验证明,小波包和分形相结合提取纹理特征的方法,在纹理分类中取得良好的效果。

2 小波包

小波分析是在Fourier(傅里叶)分析的基础上发展起来的,作为时-频分析方法,小波分析比Fourier分析有着许多本质上的进步。小波分析提供了一种自适应的时域和频域同时局部化的分析方法,在局部时-频分析中具有很强的灵活性,能聚焦到信号时段和频段的任何细节,被喻为时-频分析的显微镜[2]。

但是在小波分析中,每次都是对低频部分进行分解,所以在高频部分分辨率很低。纹理图像利用空间-频率域分析法时,纹理特征主要集中在中、高频段。所以利用小波包来提取纹理特征就给我们提供了一种更为灵活的方式[3]。

小波包和小波在纹理特征提取的时候,最重要的区别在于小波包不仅对低频部分进行了分解,而且对高频部分也进行了详细的分解。

小波包空间分解结构如图1所示。

由于正交小波分解中算子H(低通滤波算子)和G(高通滤波算子)的作用,在第一层分解中,有U3=U $_{2}^{1}$ ⊕U $_{2}^{2}$ ,U $_{2}^{1}$ ⊥U $_{2}^{2}$ 。

由于C $_{2}^{1}$ =HC3,所以U $_{2}^{1}$ 是U3中的低频分量组成的子空间;同理,C $_{2}^{2}$ =GC3,所以U $_{2}^{2}$ 是U3中高频分量组成的子空间;这样,子空间U3被分解为不同频带的两个子空间U $_{2}^{1}$ 和U $_{2}^{2}$ ;由于正交分解的缘故,所以子空间U $_{2}^{1}$ 和U $_{2}^{2}$ 是相互正交的。

类比可知,第2层分解中有:

$\begin{array}{l} U_{2}^{1} = U_{1}^{1} \oplus U_{1}^{2} ‚ U_{1}^{1} ⊥ U_{1}^{2} \\ U_{2}^{2} = U_{1}^{3} \oplus U_{1}^{4} ‚ U_{1}^{3} ⊥ U_{1}^{4} \end{array}$

同样类比可知,在第3层分解中有:

$U_{1}^{n} = U_{0}^{2 n - 1} \oplus U_{0}^{2 n}, U_{0}^{2 n - 1} ⊥ U_{0}^{2 n}, n = 1, 2, 3, 4$

从空间分解关系来看,他把正交小波分解的子空间作进一步细分;从频域划分来看,他将有限频带细分为若干更细频带的组合[2]。

3 分形

分形的概念是数学家B.B.Mandelbrot于1975年提出的,他把分形定义为“一种由许多个与整体有某种相似性的局部所构成的形体”。他提出了分形这个区别于传统的、超越尺寸的新概念,即不把微小的变化与宏观的、大的变化分离开来,而是把他们紧密联系起来。也就是说分形无所谓大尺度与小尺度,而是超越一切尺度,并存在着大小尺寸的相似性。分形理论的一个主要概念是分形维数,即维数可以是分数。

分形维数基本上刻画了分形体的复杂程度和占据空间的规模。常见的关于维数的概念有:Hausdorff维数,盒维数,关联维数,容量维数,信息维数,李亚普诺夫维数,广义维数等。本实验中,我们选用盒维数。

盒维数:用圆、球、线段和正方形以及立方体等具有特征长度的基本图形去近似分形图形。例如用长度为ε的线段集合近似海岸线那样的复杂曲线,把所得到的线段总数记为L(ε),如果该曲线有L(ε)∝ε1-D的关系,即可称曲线具有D维数[4]。

4 小波包和分形的纹理特征提取

某些只含有小面积纹理特征的子图在小波包分解中被舍弃,但在纹理的边缘,其频率成分比较丰富,这对基于小波包框架变换或者别的基于多通道滤波的纹理分析会造成一定的影响,特别是当纹理面积小时,这种影响更加不可忽略。纹理图像总体上并不一定呈现分形特征,因此分形只能从一个有限的区域和尺度来描述纹理。对于以上情况,本文章采取小波包框架和分形相结合,首先利用小波包分解,进行特征提取,某些只含有小面积纹理特征的部分因能量小而容易被忽略,对这部分利用分形提取特征值。把这两次提取的特征值共同作为图像的纹理特征值。

(1) 首先对纹理图像进行小波包分解,我们采用的是用‘db3’和‘shannon’熵标准,第一级小波包分解后,可以得到4副子图像。第二级可以得到16副子图像[5]。

(2) 求每个子图像的能量。纹理特征主要集中在中、高频段部分,所以最低频段的两幅子图的能量我们可以不考虑。求其余18个图像各自的平均能量[6]。

$E = \frac{1}{Μ Ν} \sum_{i = 1}^{Μ} \sum_{j = 1}^{Ν} x^{2} (i, j)$

M×N是小波包分解后子图像的尺寸,i,j表示子图像的行值和列值,x是小波包分解的系数[5]。

(3) 对于18幅图像中,依次排序,得到3个平均能量最低的3幅图像,求这3幅图像的分形维数。

(4) 特征向量:把得到的15副图像的平均能量,以及能量最低的3副图像的分形维数一起构成一个18维的特征向量。

(5) 采用SVM分类器[7]对纹理图像进行分类。

5 实验结果

我们利用Brodatz纹理库[8]的纹理图像作为原始图像,一共112幅。选取每个原始图像进行旋转等操作,生成尺度为64×64的子图像,一共得到2 240幅子图像。实验结果与尺度矩阵法和只用小波包分解图像提取特征值进行纹理分类的结果比较。测试结果如表1所示。从表1可以看出,本文章采取的方法有更好的平均正确识别率。

注:训练样本:1 008个;测试样本:1 232个。

6 结语

在图像纹理分析中,无论是纹理分类还是纹理分割等,最重要的就是纹理特征的提取[7]。现今的纹理特征提取的方法,都有各自的优势和缺点,因此,纹理特征提取方法的相互结合,是纹理特征提取方法的发展趋势之一。本文方法结合了小波包和分形的各自的优势,很好地进行了纹理特征的提取,但算法实现的时候,由于计算量大,计算速度还有待提高。下一步的工作还需要大量的实验,对存在的不足进行改善。

摘要：纹理是图像内容的重要特征之一,纹理图像分类中最重要的步骤是纹理特征的提取。纹理图像利用空间-频率域分析法时,纹理特征主要集中在中、高频段,纹理图像边缘包含着丰富的纹理特征信息。利用小波包和分形相结合的方法提取纹理特征值进行纹理图像分类,实验结果表明此算法是合理有效的,纹理分类的正确识别率较高。

关键词：纹理特征提取,小波包,分形,纹理分类

参考文献

[1]王丽亚,李小平,方凯.纹理图像的特征提取和分类[J].微电子学与计算机,2005,22(9):96-102.

[2]徐长发,李国宽.实用小波方法[M].武汉:华中科技大学出版社,2004.

[3]王首勇,朱光喜,唐远炎.应用最优小波包变换的特征提取方法[J].电子学报,2003,31(7):1 035-1 038.

[4]王克奇,谢永华,陈立君.基于分形理论的木材纹理特征研究[J].林业机械与木工设备,2005,33(7):19-20.

[5]王琪,费耀平.基于小波包分析的虹膜特征提取方法[J].计算机工程与应用,2006(13):60-62.

[6]刘曙光,屈萍鸽.基于小波包的织物纹理分类[J].纺织学报,2004,25(4):47-48,8.

[7]陈武凡.小波分析及其在图像处理中的应用[M].北京:科学出版社,2002.

【小波包算法】推荐阅读：

小波阈值算法09-17

小波域算法10-24