变尺度算法（共7篇）
变尺度算法 篇1
图像分割质量的好坏对后续的图像处理具有重要的影响[1], 而进化计算作为一种具有良好全局寻优性能的智能随机搜索技术, 已广泛应用于图像的分割、检索、复原和压缩等方面[2,3]。但进化计算通常采用单一的目标函数来评价图像处理的质量, 存在评价标准过于片面的不足[4]。而多目标进化优化技术作为近些年进化计算领域新兴的一个研究方向, 致力于同时处理多个相互制约的目标函数的优化问题, 已经取得了大量的研究成果[5,6]。这使得将多目标进化优化技术应用于图像处理领域具有了可行性。目前已有一些学者提出了各自的多目标进化图像处理算法[7]。如Shirakawa等人[8]提出了用于图像分割的多目标进化算法, 而Bandyopadhyay等人[9]则提出了用于遥感图像分类的多目标进化聚类算法。上述方法均直接以单个像素点作为样本来进行数据聚类, 但图像中包含的像素点往往很多, 这影响了算法的处理效率, 限制了多目标进化算法的应用范围。图像分割的目的是将图像中不同种类的对象区分开, 而图像中像素之间的差异并非都属于不同类对象的差异, 在同类对象内部实际也存在些许细节差异[10]。但这些内部差异对图像分割并无太大贡献, 反而会让算法耗费时间去计算这些细节差异, 影响了优化的速度。为了提高处理效率, 可以在早期将若干相差不大的邻域像素合并, 便于算法快速找出不同类对象的大致轮廓, 再将合并的像素逐步拆分, 最终准确确定不同对象的边界。为此, 提出了一种变尺度进化的多目标图像分割算法 (Multi-objective Evolutionary Image Segmentation with Varying-scale, MEISV) , 采用了变尺度进化的方法, 对各尺度过度的时机进行了控制, 还采用了多目标进化的框架, 并改进了连接度目标函数。对Sailboat图像的分割则表明了算法改进的有效性。
1 总体架构变尺度进化框架的设计
本文所提变尺度进化包含合并、进化、判别和分解4个步骤。框架如图1所示。
为了便于描述, 设图像I包含m×n个像素, 其中Iij (R, G, B) (i=1, 2, …, m, j=1, 2, …, n) 表示第i行第j列所在像素的R, G, B取值。合并步骤是将图像I按照所处坐标位置划分为p×q个区域, 如图2所示。其中, floor (m/p) 表示小于但最接近m/p值的整数, mod (m/p) 则表示m除以p所得的余数。然后将每个区域看作一个像素, 称为“区域像素”Igij, 其中g表示当前所在尺度的序号, 初始序号为0, 每分解1次序号加1。区域像素的R, G, B值为区域中所有像素的均值, 记为Igij (R, G, B) , i=1, 2, …, p, j=1, 2, …, q。
区域合并步骤只在整个算法的初始阶段执行一次, 通过将图像等分而得到粗粒度的像素, 此后进化算法对这些粗粒度像素进行分割, 达到一定判别准则后即进行粗粒度像素的分解, 再对分解后的像素执行进化算法, 直到所有像素均分解为单一像素且达到进化算法的迭代停止准则。
判别步骤是为了保证图像分割的准确度, 在进化到一定时期以后通过关键解的变化情况来决定是否需要对粗粒度像素进行逐步分解。所谓关键解是指相对于周围其他解, 它的某些目标函数取值仅有小幅变化, 而另外一些目标函数取值则出现大幅变化, 如图3所示。
这些解是解集中的拐点位置, 表明图像分割结果有重要的变化, 是感兴趣的解。因此当这些解变化不大时, 说明当前尺度下的进化已经足够充分, 需要分解区域以执行更小尺度下的进化。关键解的检测方法为:对多目标优化所得解集中的每一个解, 计算其与左侧第k个邻域解的斜率bk1, 以及右侧第k个邻域解与该解的斜率bk2, k=1, 2, …, k0, 然后将所有k0个bk1-bk2值的均值作为当前解的评估结果;评估结果越大说明所在位置的突变程度越大, 本文将所有解对应的评估结果降序排列, 选取评估结果最大的前z个解作为关键解, 其中z为设定的参数, 最小取值为1, 最大取值为解集的所有解。而进行区域分解的判别标准则为:在某一尺度的进化过程中, 累积t0次连续两代关键解的目标函数值的欧氏距离小于设定的阈值ε。因为小于一定的阈值则说明当前粒度下聚类所得的结果已经趋于稳定, 需要再更细一层像素粒度下进行聚类。
分解步骤是对于当前尺度下的每个图像区域, 分别以横向和纵向的中间位置为边界, 将其划分为4个子区域。然后计算每个子区域像素的R, G, B均值, 再对各个子区域相邻边界位置的像素进行局部调整。调整的方式为:对于边界位置的每一个像素, 计算其与各个子区域R, G, B均值的欧氏距离, 按照距离最近原则划归新的子区域, 最后更新各子区域的R, G, B均值。经过分解以后图像中的区域数量将变为原先的4倍, 当某一区域的像素个数等于或少于4个时, 则直接分解为单个像素。
2 多目标图像分割算法的设计
2.1 多目标函数的设计
本文所设计的多目标函数为2个, 称为紧致度指标和连接度指标。其中紧致度指标与文献[8]所提类似, 如式 (1) 所示
式中:x为进化算法的编码, 对应某一图像分割结果;wij, k表示区域像素Igij与第k个类的隶属关系, 若属于第k类则为1, 否则为0;C (x) 表示编码x所对应的总分割类数;vk则表示属于第k类的区域元素的均值[11]。而连接度指标则如式 (2) 、式 (3) 所示
式中:L=l1, l2, l3, l4分别表示区域像素Igij在其上、下、左、右4个方向上的邻域区域像素个数;表示不存在。
所提连接度指标与文献[8]的区别在于dij, nn (k) 取值的条件不同。已有方法仅仅判别区域像素与其邻域是否属于同类, 若不是则按照距离远近做出惩罚。但图像通常包含有噪声或者一些由极少像素构成的对象。这些像素无关主要对象的识别, 若将其作为一类, 则会影响分割结果的连接度指标判别, 使其取值偏大, 导致算法误以为此分割结果不合适。为此, 对惩罚条件作了修改:若区域像素Igij与其邻域Ignn (k) 不属于同类, 并且Ignn (k) 的左侧和右侧邻域区域像素中不同时存在和Igij属于同类的情况, 则做出相应惩罚, 如图4所示。
图4中“·”、“□”和“#”分别表示3类对象, “△”表示噪声像素, 原先连接度倾向于将区域 (1) 和 (3) 作为两类对象, 而改进后的连接度可以去除噪声的干扰, 正确地将其划分为一类对象。
2.2 算法描述
多目标图像分割算法是基于经典的多目标聚类算法MOCK[10]改进而得的, 用于图1所示框架的进化步骤中。其采用整数型编码x= (x11, x12, …, x1q, …, xp1, xp2, …, xpq) , 其中xij与图像中的区域像素Igij相对应。xij的取值范围为1到p×q, 按照从左至右、从上至下的顺序分别对应图像中的p×q个区域。若xij=r×s, 则表示Igij与Igrs存在连接边。通过上述编码, 图像中的不同区域像素构成了一个众多边连接的网络图。在初始化阶段, 编码的生成方式为:利用最小生成树算法得到一个网络图, 其中边的长度值为边两端R, G, B值的欧氏距离, 然后将前C个最长的边断开, 从而将所有区域像素分为C类, 其中C随机从[1, Cmax]中选取。在进化过程中, 其采用的交叉和变异算子与MOCK算法相同。这些算子通过不断的断开和建立新的连接边来对网络图进行调整, 从而形成新的子图, 使得图像分割的类数以及区域发生变化。此外, 在每次图像的尺度发生变化时, 需要对上一尺度下所得解集中各个解进行重新编码, 并重新计算相应的目标函数值。
3 实验结果与分析
为了验证算法的有效性, 选取标准图像Sailboat和Terra做测试, 如图5所示。
设置了两种对比算法, 分别为:1) 文献[8]提出的多目标聚类图像分割算法, 记为MOCKIS;2) 连接度指标未做改进, 其余与本文算法相同, 记为CMEISV。
3.1 参数设置
相关参数的设置为:MOCKIS的所有参数设置与文献[8]相同;MEISV和CMEISV的种群大小为50, 存储多目标优化所得解集的集合大小也为50, 合并步骤中p和q的取值均为8, 判别步骤中参数k0, z, t0以及ε的取值分别为3, 3, 10和10-3, 区域像素在4个邻域方向上的近邻数量l1, l2, l3, l4均取值为5, 最大类数Cmax则设为50, 迭代停止条件为所有区域分解为单个像素并且连续t0代关键解的变化量小于阈值ε。所有算法中像素的相似性度量均是以R, G, B值之间的欧氏距离作为评价标准的。
3.2 分割效果对比与分析
3种算法在Sailboat和Terra图像上的测试结果, 如图6和图7所示。图6和图7中, 子图a~c为本文所提算法MEI-SV在优化结束时所得3个关键解的分割结果, 子图d~f为对比算法CMEISV所得关键解的分割结果, 子图g则为文献[8]所提供的MOCKIS的最佳分割结果。
从图6结果显示, 本文提出的MEISV所得的3个关键解均可以将树木 (近处) 、草地、树林 (远处) 、天空和湖泊这些景物区分开, 其中子图b和c还识别出湖泊中的船, 子图c还识别出天空中的部分云彩。相比而言, 对比算法CMEISV虽然也将大部分景物识别, 但是由于连接度指标对噪声和少量像素构成的对象过于敏感, 使得图像分割以后显示出大量小的分割区域, 反而不利于主要信息的判别。对比算法MOCKIS的结果则更为粗糙一些, 其未能将图中的树林与草地、天空和云彩、湖泊和船只区分开来。此外, 由子图a~c还可以发现MEISV所得到的几个关键解随着分割类数的增加, 其细节信息也逐渐增多, 颜色相近的两个物体 (如天空和云彩) 也逐渐区分。这可以为技术人员提供不同层次下有价值的图像分割信息。
从图7中结果也可以得出类似结论:本文算法可以将远近不同的树林、田野、山峦等区分开, 而对比算法CMEISV则得到的结果过于杂乱, 缺乏层次性, MOCKIS则划分得过于简单, 一些信息被忽略 (如近处草地中包含有黄色和绿色两种草) 。
3.3 运算时间对比
除了分析算法所得图像分割结果的质量, 本文还在相同运行条件下测试了所提算法MEISV、MOCKIS和MEIS在运行时间上的差别。其中MEIS与本文算法MEISV采用相同的编码、交叉和变异以及目标函数, 只是不进行变尺度操作, 直接对各个像素点进行多目标图像分割。3种算法的迭代停止条件均设为函数调用次数达到预先设定的15 000次。测试结果如表1所示。
s
从表1结果可以看出本文算法在所用时间最短, 在图像Sailboat上为MEIS的48.08%, 为MOCKIS的48.86%, 在图像Terra上则分别为55.53%和54.80%。可见变尺度进化框架在保持求解质量的前提下, 可以有效减少优化过程所消耗的时间。
4 小结
提出的变尺度进化的多目标图像分割算法设计了先将局部像素合并再在进化过程中逐步分解的框架, 有效降低了算法的运行时间, 同时保证了图像分割的质量。其中设计的改进连接度指标对噪声和少量像素所构成的对象不敏感, 从而有效降低了算法所受到的干扰。标准图像Sailboat的实验结果表明了改进算法的有效性, 其所得分割结果能识别出图像中的主要内容, 具有层次清晰的特点, 而且算法消耗的时间最多减少达一半以上。
参考文献
[1]江成顺, 汪先超.两相图像变分分割凸松弛模型快速算法[J].计算机学报, 2013, 36 (5) :1086-1096.
[2]靳璐, 付梦印.基于遗传模糊核聚类的图像分割[J].模式识别与人工智能, 2013, 26 (1) :205-210.
[3]任通, 程江华, 金阳.一种基于白像素增量比的字幕图像分割算法[J].电视技术, 2014, 38 (5) :190-193.
[4]李永乐, 娄静涛, 张茂军, 等.基于全向总变分最小化的折反射散焦模糊图像复原方法[J].光学学报, 2013 (8) :99-105.
[5]王璐, 胡晰远, 彭思龙.基于分块的空间变化抖动模糊图像的全局模糊去除[J].计算机辅助设计与图形学学报, 2012, 24 (6) :766-774.
[6]王冬冬, 张炜, 金国锋, 等.尖点突变理论在红外热波检测图像分割中的应用[J].红外与激光工程, 2014, 43 (3) :1009-1015.
[7]李磊, 金良海, 宋恩民.基于区域合并和图割的非监督彩色图像分割算法[J].华中科技大学学报:自然科学版, 2013, 41 (12) :44-47.
[8]SHIRAKAWA S, NAGAO T.Evolutionary image segmentation based on multiobjective clustering[C]//Proc.2009 IEEE Congress on Evolutionary Computation.NY, USA:IEEE Press, 2009:2466-2473.
[9]BANDYOPADHYAY S, MAULIK U, MUKHOPADHYAY A.Multiobjective genetic clustering for pixel classification in remote sensing imagery[J].IEEE Trans.Geoscience and Remote Sensing, 2007, 45 (5) :1506-1511.
[10]HANDL J, KNOWLES J.An evolutionary approach to multiobjective clustering[J].IEEE Trans.Evolutionary Computation, 2007, 11 (1) :56-76.
[11]陈恺, 陈芳, 戴敏, 等.基于萤火虫算法的二维熵多阈值快速图像分割[J].光学精密工程, 2014, 22 (2) :517-523.
变尺度算法 篇2
微网是一种由负荷、微电源(分布式电源)和储能装置共同组成的有机系统。微网可以有效地整合各种新能源发电技术,充分发挥新能源发电所带来的经济效益和环境效益;可以更好地满足用户对电能质量和供电可靠性更高的要求;可以实现电网的削峰填谷,提高能源的利用率[1,2,3]。
微网的优化运行是微网研究的重点和难点问题,已经引起了人们的广泛关注。微网的优化运行主要有如下几个特点:由于微网主要接入配电网中,电压等级低,系统中输电线路的电阻起主导作用,线路损耗较大,不能忽略[4];由于风力发电、光伏发电等新能源发电技术受自然环境的影响较大,输出功率具有随机性、波动性、间歇性,微网的优化运行要求实时性比较强,一般在秒级;由于微网内的微电源与传统的火电机组有很大的区别,不能简单应用等耗量微增率准则,一般采用智能优化算法来进行优化。
混沌优化算法具有遍历性、随机性、规律性的特点,能在一定的范围内按照自身的规律不重复地遍历所有的状态[5]。混沌优化算法能避免陷入局部极小,比随机搜索更具有优越性,易于跳出局部最优解[6]。文献[7]将混沌优化算法用于电力系统经济负荷分配,并与遗传算法相比较,证明了混沌优化算法的优越性。文献[8]首次提出了变尺度混沌优化算法,并通过数值算例验证了所提算法的有效性。文献[9]对变尺度混沌优化算法进行了改进,并将改进后的算法应用于大规模电力系统经济负荷分配的算例中,与其他算法相比取得了更优的结果。
本文研究了变尺度混沌优化算法的初值、变量空间缩小系数、“二次搜索”调节系数对优化结果的影响,并提出了相应的改进措施,通过数值算例验证了改进措施的有效性。结合微网中微电源的特点建立了微网优化运行的数学模型,应用改进后的算法进行优化,取得了预期的效果。
1 对变尺度混沌优化算法的改进
1.1 参数的影响
变尺度混沌优化方法的具体算法步骤详见文献[8]。从理论上来讲,由于混沌变量的遍历性特点,不受初值的影响,但是迭代步数是有限的,混沌变量也就不可能遍历所有的状态,不同的初值会出现不同的结果。针对这一问题,文献[9]提出了将初值选为随机数和并行搜索的方法,但是该方法可能会使计算时间明显增加。
变量空间缩小系数和“二次搜索”调节系数选择得不恰当,可能会造成无法找到全局最优点,而只能找到局部最优的情况。同时,变量空间缩小系数和“二次搜索”调节系数的选择与目标函数有很大的关系,目标函数不变,改变缩小系数或调节系数,可以得到不同的优化结果。
1.2 算法改进
本文在文献[8-9]的基础上对变尺度混沌优化算法中变量空间缩小系数、“二次搜索”调节系数的选择进行了改进,减少了这些参数对不同目标函数优化结果的影响。
1.2.1 变量空间缩小系数的选取
在变尺度混沌优化算法中,需要根据搜索进程不断地缩小优化变量的搜索空间[8],变量空间缩小系数就是表征优化变量的搜索空间在每次“二次搜索”过程中缩小程度的参数,本文用t来表示。
文献[8]给出了变量空间变化的公式,如式(1)所示:
其中,i表示第i个变量,r表示“二次搜索”的次数,xir表示第r次“二次搜索”得到的第i个变量的最优值,bir、air分别表示第i个变量在第r次“二次搜索”时的上限和下限,bir+1、air+1分别表示第i个变量在第r+1次“二次搜索”时的上限和下限。
文献[8]指出t的范围是(0,0.5),同时随着“二次搜索”次数r的不断增加,变量的搜索范围也在不断缩小。由式(1)可以看出变量的搜索空间在当前最优解附近进行调整,搜索空间的调整程度与变量搜索空间的大小有关,当搜索空间较大时,t应取较大的值,以保证搜索的速度;当搜索空间较小时,t应取较小的值,以保证搜索的精度。
为此本文提出了参数t的确定公式,如式(2)所示:
由式(2)可以看出,t的范围在(0,0.5)之间,满足文献[8]的要求。同时,随着“二次搜索”次数r的增加,变量搜索空间不断减小,t的数值也将不断减小。
1.2.2“二次搜索”调节系数的选取
文献[8]指出在变尺度混沌优化算法中,还需要根据搜索的进程不断改变“二次搜索”调节系数。“二次搜索”调节系数是指在粗搜索得到的次优点的基础上进行微调,得到新的混沌变量,用新的混沌变量进行“二次搜索”,本文中,“二次搜索”调节系数用α表示。文献[8]给出了变尺度的公式,如式(3)所示:
其中,k表示第k次混沌搜索;xir,k+1表示在第r+1次“二次搜索”中,第i个变量在第k次混沌搜索中用到的混沌变量;xir表示第r次搜索得到的第i个变量的次优点;xir,k表示在第r次“二次搜索”中,第i个变量在第k次混沌搜索中用到的混沌变量。
由式(3)可以看出α应该是一个与“二次搜索”次数r相关的数,并且α的取值应该是一个较小的数,以保证在次优点附近进行微调。同时,随着“二次搜索”次数的增加,寻优结果不断向真值靠近,α应该不断地减小以保证寻优结果的精度。
为此本文提出了参数α的确定公式,如式(4)所示:
1.3 数值仿真
本文采用2个常用的测试函数F1、F2对算法进行测试[8,9],所有的数值仿真和实例计算都在MAT-LAB中编程实现。
测试函数F1如下:
其中,变量的取值范围为:-2.048≤x1≤2.048,-2.048≤x2≤2.048。此测试函数的理论最优解为min(1,1)=0。
测试函数F2如下:
其中,变量的取值范围为:-2≤x1≤2,-2≤x2≤2。此测试函数的理论最优解为min(0,-1)=3。
本文在进行测试的过程中,选定了相同的迭代步数,即算法耗时是一定的。表1是本文算法和文献[8]算法的对比结果。
由表1可以看出本文所选取的变量空间缩小系数t、“二次搜索”调节系数α,在迭代步数一定的情况下,能够得到更优的解,证明了本文算法的寻优效果更好。
1.4 实例验证
以文献[7]中的3机6母线系统为例,总负荷为500 MW,不考虑阀点效应和网损的情况下,惩罚因子选为10,在机组之间进行负荷的优化分配。各单元机组的参数如表2所示。本文算法运行3次得到的可行解,与文献[7]、文献[8]中算法的结果进行比较,如表3所示。由表3可见,本文算法得到的负荷分配结果比文献[7]和文献[8]的算法得到的结果费用更低,证明了本文算法的有效性。
2 微网优化运行的数学模型
2.1 微网参数
本文采用简化的微网模型,该微网共有10个节点,具体结构如图1所示。微电源有微型燃气轮机(MTG)、柴油发电机(DEG)、蓄电池(BAT)、光伏电池(PV)、超级电容器(SC)。其中PV和SC的输出功率作为已知量处理,不作为优化变量。PV的实际出力设定为3 kW,SC的实际出力设定为2 kW。微电源容量如下:MTG为20 kW,DEG为15 k W,BAT为15 kW,PV为5 kW,SC为5 kW。负荷容量如下:负荷1为10+j0.1 kW,负荷2为15+j0.17 kW,负荷3为7.5+j0.05 kW,负荷4为7.5+j0.08 kW。本文在做微网优化运行时,只考虑了微网孤岛运行的情况,所以设定并网节点9的功率为0。
2.2 目标函数
微网的优化运行是一个多目标、多约束条件的复杂优化问题[10]。本文同时考虑了微网的经济成本最小和网损最小作为目标函数,其中经济成本主要考虑了燃料成本、运行维护成本、环境折算成本。给不同的子目标函数赋予不同的权重,进行线性加权,将多目标问题转化为单目标问题,同时采用罚函数的方法对约束条件进行处理。
2.2.1 微网经济成本
2.2.1. 1 燃料成本
a.MTG的燃料成本与自身的工作效率有关,表达式如式(7)所示:
其中,FMTG为MTG的燃料成本;C为MTG采用的燃料气体的单价,本文取2元/m3;VLH为天然气的低热热值,本文取9.7 kW·h/m3;PMTG为MTG的输出功率;ηMTG为MTG的效率,其大小与MTG输出功率的大小有关[11]。
b.DEG的燃料成本就是它的耗量特性函数,如式(9)所示[12]:
其中,参数a、b、c的大小一般由生产厂家给定,本文选取a=6,b=0.012,c=8.5×10-4。
c.根据BAT的特性可知,BAT不消耗燃料,不存在燃料费用。
2.2.1. 2 运行维护成本
微电源的运行维护成本可以用微电源输出功率乘以相关的系数来表示,如式(10)所示:
其中,E为微电源总的运行维护成本,ki为第i个微电源的运行维护成本系数,Pi为第i个微电源的输出功率,N为微电源的数目。ki的大小按照文献[12]选取,具体如式(11)所示:
2.2.1. 3 环保折算成本
MTG和DEG在运行的过程中会产生氮氧化物(NOx)、二氧化硫(SO2)、二氧化碳(CO2)等空气污染物。考虑到微网的环境效益,将这些污染物按照一定的成本进行折算,作为微网优化运行的目标。具体计算公式如式(12)所示:
其中,C1为微网的环保折算成本,Pi为第i个微电源的输出功率,N为微电源的数目,aij为第i个微电源排放的第j种污染物的量,q为污染物的种类,cj为第j种污染物的折算成本。
不同种类的污染物折算成本以及MTG、DEG的排放因子如表4所示[13]。蓄电池由于其自身的运行特性,不产生污染物。
综合考虑以上因素,微网的经济成本如式(13)所示:
其中,Vcost表示微网的经济成本。
2.2.2 网损
由于微网一般都接在电压等级较低的配电网中,而配电网中R/X的值一般较大,在5到几十之间[14],因此微网中传输线的网损比传统大电网的网损明显加大,不可忽略不计。网损通过潮流计算的方法得到,如式(14)所示:
其中,Pk、Qk为第k条支路传输的有功、无功功率,M为支路总数,Rk为支路k电阻,|Uk|为支路电压幅值。
2.3 约束条件
功率平衡约束:
其中,Pi为第i个微电源输出的功率,N为微电源的数目,Pload为总负荷。
微电源输出功率约束:
其中,Pimin、Pimax分别为第i个微电源输出功率的下限和上限。
节点电压约束:
其中,Ui为第i个节点的电压,U imin、Uimax分别为第i个节点的电压下限和上限。
3 算例仿真
本文选取的微网电压等级为380 V,线路选择L J-16型导线,线路阻抗为R=1.98Ω/km,X=0.358Ω/km[15,16]。采用本文提出的改进变尺度混沌优化算法进行计算。目标函数为F=λ1Vcost+λ2Ploss,λ1、λ2分别为多目标的权重系数,并且满足λ1+λ2=1。
算法步骤如下。
a.初始化。输入优化变量的维数N,各个微电源的输出功率上下限bi、ai,“一次搜索”最大迭代步数,“二次搜索”最大迭代步数,权重λ1、λ2的取值,惩罚因子等参数,并随机生成N个混沌变量。
b.一次搜索。将混沌变量映射到待优化变量的取值范围内,用混沌变量进行搜索。
c.判断“一次搜索”迭代步数是否满足“一次搜索”的最大迭代步数。若不满足则继续迭代搜索,若满足则进行步骤d。
d.变尺度。对“一次搜索”得到的当前最优解进行变尺度得到新的优化变量,并调整各变量的搜索空间。
e.二次搜索。用新的混沌变量在调整后的搜索空间内执行步骤b、c。
f.重复执行步骤d、e。
g.判断“二次搜索”迭代步数是否满足“二次搜索”的最大迭代步数。若不满足则继续迭代搜索,若满足则输出微网优化运行结果。
对λ1、λ2赋予不同的值,计算结果如表5所示。
由表5的结果可以看出,考虑不同的权重,将得到不同的优化结果,各个微电源的输出功率也不同。在进行微网的优化运行时,应结合不同用户的需要,综合考虑不同的影响因素。
4 结论
本文对变尺度混沌优化算法进行了改进,并利用典型的数值算例进行了验证,证明了改进的有效性。将该方法应用到电力系统负荷优化分配中,取得了比传统的变尺度混沌优化算法更好的效果。
变尺度算法 篇3
随着经济的快速发展, 人们对能源的需求与日俱增。面对不可再生能源日趋减少和全球变暖等环境污染问题, 迫切需要一种可稳定持续开发的可再生能源。由地球表面空气流动所产生的风能无疑是最好的选择之一, 这就使风力发电成为一个重要的研究课题[1]。其中, 短期功率预测[2]是风力发电系统研究的一个重要组成部分, 它对发电系统的功率控制和稳定运行起着至关重要的作用。
目前, 在风力发电系统功率预测方面, 国内外学者都做了一定的研究工作。由于国外对这方面的研究起步较早, 很多国家都已拥有比较完善的风电功率预测模型和系统。而国内的风电产业起步较晚, 在此的研究还局限于风电场的规划、发电机组的控制以及风电场以最小影响并入电网等方面, 对风力发电的短期预测研究成果甚少, 预测效果不佳[3]。
本文在现有研究的基础上, 结合国内数值天气预报不完备等特点, 提出一种将历史数据与物理模型相结合的插值方法。并针对风力发电系统短期功率预测因预测信息的不同对预测周期的不同需求, 提出一种在已有插值方法基础上的变时间尺度短期功率预测方法。
2 短期功率预测方法
物理模型与历史数据结合基础上的变时间尺度功率预测方法由发电机组功率预测输入预处理模型、基于历史数据的天气预测模块以及功率预测模块组成[4]。
2.1 风力发电功率预测模型输入处理
依据物理模型对风力发电短期功率预测的输入处理步骤可总结如下: (1) 通过数据校正模型对已有数值天气预报数据进行校正, 减小预测误差; (2) 通过地形变化模型修整因风电场地表凹凸不平给预测精度带来的影响; (3) 建立基于不同地貌的粗糙度分类模型, 计算受粗糙度影响后的风速; (4) 通过风轮背面风速模型求解风功曲线对应的风速。
2.2 变时间尺度的天气预测模型
数值天气预报风速经过数据校正、地形和粗糙度模型修正后可作为距离最近风力发电机组的风速值。但鉴于短期功率预测的预测周期存在小于当前所研究风电场天气预报信息的分辨周期, 达不到功率预测周期的要求的情况, 提出一种基于历史数据的插值方法, 提高天气预报信息的分辨率。又根据控制系统对短期功率预测预测周期的不同要求, 提出一种基于历史数据的变时间尺度的天气信息数据预测模型。
2.3 功率预测模块
在系统执行能量分配时, 根据功率预测模块对预测周期的要求设置插值周期, 再按此周期从天气信息历史数据中提取数据。若根据控制系统要求需要不同的预测周期, 则对应设置不同的插值周期, 提取不同插值的历史数据, 至此实现变时间尺度的插值运算, 完成变时间尺度的风力发电机组功率预测物理模型。
执行预测时, 引入相似度法则对历史风速变化趋势与预测点的预测风速和实际风速组成的综合数据做相似度分析计算。其具体方法如下:
以日特征向量反应风速变化趋势, 假设为预测日的风速综合数据, 为历史第j天的风速综合数据, 为预测日的特征向量, 其中, k表示预测日k时刻相邻时间点的风速变化趋势量, , k表示历史数据第j日k时刻相邻时间点的风速变化趋势量, 则预测日数据与历史第j日数据变化趋势相似度可由如下表达式求解。
通过对预测日风速变化趋势与历史风速变化趋势做相似度结果的比较, 从中选取相似度较高的m天的数值天气信息历史数据, 并以求加权平均的形式对预测日当天的风速进行插值。求解获得预测日下一预测周期的预测风速, 达到为接下来基于历史数据的变时间尺度风力发电功率预测提供可靠天气信息支持的目的。
2.4 风力发电机组预测功率
数值天气预报信息经过各输入预处理模型, 结合变时间尺度功率预测模块, 成为距离最近风力发电机组的输入风速值。在无干扰情况下, 根据风电场空气密度、发电机组相关参数以及风速值可推到出理想输出功率如下。
式中, P代表发电机组输出功率, R代表发电机组桨叶半径, 这δ代表风电场空气密度, Dp代表风能利用系数, vr、vin和vout分别代表发电机组的额定风速、起始风速和终止风速。
在现实发电工程系统中, 发电系统受环境等诸多因素影响, 很难达到理想状态, 故需采用实际风电场功率风速曲线进行功率输出预测。
3 结论
在本文中, 考虑到目前国内数值天气预报分辨率不足以满足短期功率预测周期要求的情况, 提出一种将物理模型与历史数据结合的相似度插值方法, 通过与线性插值法的对比试验, 验证了该方法的可行性;并根据风力发电系统对功率预测周期的不同要求, 提出一种基于历史数据的变时间尺度方法。
参考文献
[1]王常贵.新能源发电技术[M].北京:中国电力出版社, 2003.
[2]杨秀媛, 肖洋, 陈树勇.风电场风速和发电功率预测研究[J].中国电机工程学报, 2005 (11) :1-5.
[3]李炎, 高山.风电功率短期预测技术综述[A].中国农业大学.中国高等学校电力系统及其自动化专业第二十四届学术年会论文集 (下册) [C].中国农业大学:2008:5.
变尺度算法 篇4
Marching Cubes算法(简称MC算法)是由William E.Lorensen和Harvey E.Cline在1987年的Computer Graphics中第一次提出来的,它是面显示算法中的经典算法,也被称为“等值面提取”(Isosurface Extraction)算法,它根据阈值在一个3维数据场中抽取等值面实现3维重建。自从MC算法提出以来,通过不断改进,已广泛地应用在可视化的各个领域,特别是在医学图像三维重建、立体地图测绘等领域。MC算法提出以后,一直就被认为是迄今为止应用最广泛的面重建算法之一。
2. MC算法简介
在医学图像三维重建中,一张CT图像为一个断层切面,两张CT图像组成上下两层切面,每张断层切面取相邻的4个点,共8个点组成了一个立方体(Cube),如图1所示。Marching Cubes算法先单独计算一个立方体与等值面相交的情况,得到若干个三角面,然后再遍历所有的立方体,得到所有的相交三角面,最后将这些三角面拟合为等值面输出。
在原始数据(比如CT、MRI三维图像)中,它们的三维空间数据通常是离散的,沿着立方体边的方向数据场呈连续线性变化,如果立方体的一条边的两个顶点分别大于、小于等值面的值,则在该边上有且仅有一个交点。如果边上的两个顶点都大于或都小于等值面的值,则这条边和等值面没有交点。它的基本算法是逐个处理数据场中的立方体,找出等值面与这些立方体的交点,进而找出包含在这些体元中的等值面。
Lorensen在他们的论文中第一次提出MC算法时,具体分析了Cube与等值面相交的各种情况:因为一个立方体有8个顶点,每个顶点有2种状态,即表面内和表面外,因而Cube与表面的关系有2的8次方,即256种情况,通过枚举这256种情况,可得到一个状态表以得到等值面与Cube各边相交的信息以及Cube内表面的拓扑结构。因为Cube是立方体,利用2种不同的对称性,即表面与Cube上每个顶点与表面的关系同时逆转,不会影响Cube内三角化表面的拓扑结构,也就是说256种情况是两两互补的,所以只需要考虑0-4种顶点在表面内的情况。另外,根据立方体3个方向旋转上的对称性,即Cube向3个方向任意旋转90度的倍数,Cube内表面的三角形拓扑结构不变。这样经过分析,最终将256种情况减少为15种。这15种情况如图2所示:
Marching Cubes具体的算法步骤:
(1)根据对称关系构建一个256种相交关系的索引表。该表指明等值面与体素的哪条边相交。
(2)取相邻两层图片中相邻的8个象素,构成一个体素并把这8个象素编号。
(3)比较每个象素与阈值确定该象素是l还是0。
(4)把这8个象素构成的01串组成一个8位的索引值。
(5)用索引值在上边的索引表里查找对应关系,并求出与立方体每条边的交点。
(6)交点构成三角形面片。
(7)遍历三维图像的所有体素,重复执行(2)到(6)。
(8)将得到的所有的三角面片输出。
3. MC算法的缺点
Marching Cube算法的优点是明显的,但它也存几个缺点:
(1)MC算法局限于通过阈值分割来提取等值面,而医学图像由于成像机理不同,所适用的分割方法也不同,这就造成了MC算法适用的局限性;
(2)MC算法中关于等值面与立方体交点的计算是采用插值法的,计算量比较大;
(3)MC算法存在二义性。如果在体元的一个边界面上,标号为0和标号为1的角点分别位于对角线的两端,则该表面上存在四个等值点,那么就会有两种可能的等值线的连接方式,这就是二义性;
(4)MC算法的输出是三角形网格,由于得到的数据分布十分密集,由此建立的三维网格数量庞大。
4. MC算法在不同尺度上的改进
根据上面的分析,MC算法由于计算量非常大,所以显示三维面时会适成延时。本文根据MC算法可以选择不同的像素点组成立方体进行求交点提出了一种多尺度的MC算法。其原理如下:
设每相隔N个点取一个顶点,上下两层共8个顶点组成一个大间距的立方体,从这个大间距立方体开始计算。下图3为同一个方向上的两个顶点a1,a2。求这个两个顶点组成的边与等值面的交点K。
根据MC算法求交点的公式,K的坐标为:
K=(a2-a1)(阈值-v1)/(v2-v1)(1)
遍历完所有大间距的立方体后,就得到了一个大尺度的交点集,由这个交点集所组成的三角面片会明显减少,计算量下降,使得三维图像绘制速度加快。
计算完上述大尺度的立方体后,缩小顶点之间的间距,将立方体的边长变为原大尺度的一半。下图4为缩小一半间距所同一方向上的三个顶点。
则以a1,a3组成边的两个顶点,a3,a2组成另一组边的两个顶点。则这两条边和等值面的交点求法如下:
K1=(a3-a1)(阈值-v1)/(v 3-v1)(2)
K2=(a2-a3)(阈值-v3)/(v 2-v3)(3)
观察(2)式和(3)式有
K=(a2-a1)(阈值-v1)/(v2-v1)
K1=(a3-a1)(阈值-v1)/(v3-v1)
其中a3-a1=(a2-a1)/2
所以有:
K(v2-v1)=(a2-a1)(阈值-v1)
K1(v3-v1)=(a2-a1)(阈值-v1)/2
即
K1(v3-v1)=K(v2-v1)/2
K1=K(v2-v1)/2(v3-v1)
同样K2=K(v2-v1)/2(v2-v3)
根所上述理论,计算一次小尺度的交点K1,K2,都可以利用到上一次大尺度的计算结果K,而v2,v1,v3都是已知值,所以上述计算过程将会非常之简单。这样就从理论上找到了一个可以提高MC算法计算速度的方法。以下表1是本文根据这个理论所做的实验,实验的数据是来自田捷的中科院自动化研究所医学影像研究室提供的一个老鼠的CT图,该数据的大小是512层,每层是512×512。每个像素点是32bit的值。
由表1可以看出,采用不同尺度的方向后,生成的三解形面片少了,速度也快了。对比图5和图6可以看出,用两个尺度相差数倍组成立方体的MC算法计算出的三维图像的效果差别是不明显的。
5. 结论
本文根据MC算法计算量大,难以实时绘图的缺点,提出了一个多尺度的MC算法,经过实验对比,新的算法比原算法速度更快,生成三维面的效率更高,而且生成三维面的质量并没有明显的下降。
参考文献
[1]徐毅,李晓梅.对体可视化Marching Cube算法的改进[J].计算机工程,1999,25(11):52—54.
[2]谢小棉,李树祥.基于MC的医学三维等值面的平滑与归并[J].中国图像图形学报,2001,6(8):806—809.
[3]William E.Lorensen,Harvey E.Cline,Marching Cubes:A High Resolution3D Surface Construction Algorithm,1987
[4]Doi A,Koide A.An Efficient Method of Triangulating Equi-Valued Surfaces hy Using Tetrahedral Cells,1988.
变尺度算法 篇5
关键词：小波变换,MALLAT算法,多尺度分解
小波变换是一种进行信号分析的时频分析方法, 除被广泛应用于图片、音频和视频分析处理外, 还被广泛应用于计算机科学、小波降噪和薄层划分等领域, 是一种实用性强且方便高效的信号处理工具。但当前多数情况下小波分析方法大都是直接应用MATLAB软件上面的小波工具箱, 不利于包装成具有特定功能的分析工具, 而采用C语言实现的小波分析算法缺乏交互界面, 可用性不强。
本文MALLAT算法原理, 采用当前较为流行的、跨平台性较好的JAVA语言编程实现了小波多尺度快速分解算法, 为开发面向对象的小波多尺度分析算法提供参考。
1 多尺度及MALLAT分解算法
离散小波变换是将输入信号通过与高低频滤波矩阵进行信号分解, 其数学表达式如下:
式中, an+1 (2m) 为an (2m) 经低频滤波器h (2m) 数字滤波之后的偶抽取样本数;dn+1 (m) 为an (m) 经高频滤波器g (2m) 数字滤波之后的偶抽取样本数;h (2m) , g (2m) 分别为低频滤波器H和高频滤波器G的冲击响应序列, K为位移。
2 相关算法的JAVA编程实现及代码
2.1 数据初始化
MALLAT算法所具有的数组之间的运算关联性, 有利于在进行数据初始化时将数据读入到指定的动态数组中, 便于使用JAVA中的动态数组实现数组内成员的自主动态管理。
2.2 滤波器组
因小波滤波器之间存在镜像与倒序的关联, 所以在进行小波分解时, 只要低频分解滤波器矩阵, 就能根据其倒序和镜像算法计算出其他的滤波器。
2.3 边界延拓
MALLAT算法假定输入的数据信号是无限长的, 而实际中的数据信号往往是有限长的, 所以在处理实际数据信号时, 为降低边界信号突变的影响, 一般都涉及边界信号的延拓问题。信号延拓的常用方法有周期延拓、对称延拓和零延拓。其数学表达式为:
式中, f′ (x) 为f (x) 经对称边界延拓之后得到的延拓序列, filter_LEN为对应滤波器长度, s_LEN为信号长度。
2.4 小波分解
MALLAT算法的实质就是用实际信号与低频和高频信号进行卷积运算, 然后进行偶抽取, 从而得到相应的细节信号和逼近信号, 故而在编程实现过程中只要给定相应的滤波器和数据信号, 就能够运用统一算法进行分解。其算法分解实现过程如下:
3 结语
本文讨论了一维离散小波多尺度变换MALLAT算法及JAVA语言编程实现的原理和方法。并结合MALLAT快速分解原理及小波重构原理制作了独立的小波多尺度分析软件, 给出了进行多尺度分解主要步骤和算法, 能够为通过面向对象语言开发小波能谱分析相关的算法提供指导。
参考文献
[1]侯正信, 王成优, 杨爱萍.有限长度信号Mallat算法的边界延拓方法[J].数据采集与处理, 2009 (6) :714-720.
[2]岳文正, 陶果.小波变换在识别储层流体性质中的应用[J].地球物理学报, 2003 (6) :863-869.
变尺度算法 篇6
1 Contourlet 变换
Contourlet变换先进行多分辨率分解然后在进行多方向分解,首先用拉普拉斯金字塔(LP,Laplacian Pyramid)变换进行多尺度分解,并捕捉图像信号的奇异点,接着使用方向性滤波器组(DFB,Directional Filter Bank)进行多方向分解,将同个方向上的奇异点合成,捕捉二维信号的高频成分。图1给出了Contourlet变换的流程图。LP变换能够有效的避免由方向性滤波器所造成的低频分量的“泄漏”。
Contourlet变换的核心LP和DFB,因此Contourlet变换也称“塔式方向滤波器组”。而Contourlet变换的最终是用类似于轮廓段 (contour segment)的基结构来逼近表示二维图像信息。图2表示一个Contourlet变换可能分解的频率。由于拉普拉斯金字塔和方向性滤波器组具备的完全重构特性,因此其组合PDFB也具备能完全重构。但是由于拉普拉斯金字塔本身具有冗余性,因此Con-tourlet变换具有冗余性。
图像经Contourlet分解之后绝大多数系数幅值接近零,系数表现相当稀疏,而在原图像的边缘细节以及纹理特征部分稀疏幅值比较大,且在他们在各尺度上具有有一定的相关性。图3给出了lena图的Contourlet分解的例子,图像经3级LP分解,以及8个方向的分解,从图我们可以看到最大多系数都接近零,而在图像的轮廓和边缘部分系数的幅值比较大。
2 基于层结构的多阈值Contourlet降噪算法
阈值降噪是图像降噪处理中最常用的方法,其关键是选取一个适当的阈值把图像噪声系数和信号系数区分开来,然后对噪声系数进行置零或者修改其值,随后通过反变换重构图像,达到图像降噪目的。
文献[8]中提出了基于层结构的Contourlet多阈值图像降噪算法。文献中提出“算法使用Contourlet变换来代替小波变换,Con-tourlet的基函数能有效地对分段光滑的线段进行表示,而对奇异点则影响不大”,在降噪处理的过程中,即使将噪声系数误认为图像细节系数给予保留,或者相反,将图像细节系数误认噪声给予删除,也不会使降噪重构的图像出现显著数值较大的孤立点,从而造成图像边缘细节或纹理特征稀疏的丢失。其方法如下:
1)首先,对带噪声的图像进行CT变换,确定对图像进行分解的级数J和各个级数J上所要分解的方向数s。
2)根据下面式子计算出各个尺度和各个方向上噪声方差估计值,
wgi,j(s) 为带噪声图像高频系数。
3)对图像CT变换后的系数方差进行估计
4) 阈值调整:对低分解层的高频系数进行调整。
5) Contourlet反变换:将处理后的Contourlet变换系数进行反变换,可以重构经过降噪之后的图像。
3 基于Contourlet的相关尺度降噪算法
文献[10]提出了一种基于Contourlet变换的相关降噪新算法,其思想主要是将图像的Contourlet系数跟含噪声图像的Contourlet变换系数进行比较,随着分辨率越高,图像的细节部分对应的Contourlet变换系数越稳定,而噪声则衰减越大,即代表图像纹理特征或边缘细节的Contourlet系数相关比较强,而代表噪声的Contourlet变换系数相关比较弱或不相关。采用了利用小波系数相关量来确定图像边缘信息,并且与阈值函数相关结合。其过程如下:
4 基于Contourlet相关尺度多阈值降噪算法
采用硬阈值函数进行降噪处理后的信号虽然可以比较好地保留图像边缘纹理等细节信息,都由于这种理想滤波器会造成图像出现振铃或者伪吉布斯效应等视觉失真,而基于相关尺度的降噪方法,有时候易将过多噪声系数误判为有用信号,为了取得比较好的降噪效果,该文对以上两种方法进行综合,采用了相关尺度多阈值降噪方法,将阈值函数与相关尺度结合起来。在基于层机构的阈值基础上,对于判断为噪声的Contourlet系数进行相关量的计算,具体算法如下:
1)首先对带有噪声的图像进行CT变换,一般情况下确定对图像进行分解的级数J为3级和各个级数J上所要分解的方向数s为8。
2)运用文献[8]的方法,计算3个尺度以及8个方向上噪声方差估计值,并对对图像的Contourlet系数的方差进行估计最后得到阈值函数:
5 结果及分析
本文采用标准的大小为512*512像素测试图象Lena, Peppers为实验图象,进行仿真实验,比较了多阈值降噪、相关降噪和本文提出的多阈值相关降噪效果的比较,分解级数为4,4,8;表1为各个降噪方法后的信噪比(SNR),图4是lena图各种方法降噪后的结果。
从表1可以看出虽然尺度相关多阈值降噪的SNR不是最高,但在局部图上可以看出此方法在边缘细节的保持上效果更好。
6 结论
变尺度算法 篇7
图像去噪是图像处理领域基本方法之一,其目的是提高信噪比、改善主观视觉效果。在图像去噪处理研究中,寻找一种即能去除图像中的噪声又能较好地保持图像的边缘细节信息是图像处理学、计算机视觉和模式识别研究中的热点问题[1]。小波分析兴起于20世纪80年代,小波变换不但可以提供比较准确的时域定位,而且能够提供较精确的频域定位,而傅里叶分析虽然具有很强的频域局部化能力,但并不具有时间局部化能力。一般认为,小波的兴起得益于Yves Meyer和Mallat等人引入的多分辨分析的概念、紧支撑小波的构造及快速小波算法的实现,另外一个重要的原因是小波对具有点奇异的目标函数具有最稀疏的逼近性能。对于二维图像而言,图像的主要信息由边缘和纹理特征所描述,边缘和纹理是图像函数中的“线奇异”。常用的二维小波由一维小波张成,这种小波只有“水平”、“垂直”、“对角”3个有限方向的空间辨析能力。研究表明,理想的高维函数表示方法其基本元素至少应该是多方向的,或者说“各向异性”,正是空间方向辨析能力的不足使得二维张量小波对含“线奇异”的函数不能达到最优的非线性逼近阶。在高维情况下,小波分析不可以充分应用数据本身所具有的几何特征,小波变换在高维情况下,还不是最好的或者“最稀疏”的表示方法[2,3]。多尺度几何分析发展的目的和动力正是要努力寻找一种新的高维函数的表示方法。
曲线波变换是在脊波变换基础上发展起来的一种新的多尺度变换,笔者在研究Curvelet变换的基本理论及其在数字图像处理应用基础上,提出新型方向及尺度乘积的曲线波去噪算法。实验表明,笔者提出的新算法在图像去噪方面具有更好的视觉效果,峰值信噪比(PSNR)不但达到较高值,同时能明显减轻划痕现象的产生。与传统的小波去噪和Curvelet阈值去噪方法进行相比较,新算法在主观和客观评价指标上,体现了明显的优越性。
2 Curvelet变换及其离散实现
2.1 Curvelet变换的产生背景
在高维情况下,因为小波分析有限的方向性,由一维小波张成的可分离小波并不可以“最优逼近”具有线或者面奇异的函数。为了能够更好地逼近具有直线奇异的多变量函数,由E.J.Candès在1998年提出了脊波理论和基本理论框架,然而对于含曲线奇异的多变量函数来说,它的逼近性能仅和小波函数相当,不具备最好的非线性逼近误差衰减阶。对于多变量函数f∈L1∩L2(Rn),脊波变换的定义为:
若函数Ψ满足特定容许条件,则称Ψ是容许神经元激励函数,并称
为脊波,脊波变换为
R(f)(捃)=
脊波非线性逼近误差为
2.2 Curvelet变换
脊波是一个非常优秀的函数分析工具,然而相对于普通的目标函数(例如图像)来说,边缘是图像的不连续所在,通常情况下,边缘不是直线型的。在二维情况下,标准正交的脊波分析等价于Randon域中的非正交小波分析。针对具有曲线奇异性的目标函数来说,脊波分析并不是最稀疏的表示方法。
鉴于此,E.J.Candès和D.L.Donoho等人提出了曲线波(Curvelet)变换,本质上是多尺度局部的脊波变换形式。它的主要思想是,对曲线进行无限分割时,每一小段大体上能够看作为直线段,这时,即能够采用脊波分析理论。当对所有的尺度s≥0进行脊波变换时,对于具有光滑奇异性的目标函数,曲线波提供近似最优的表示。
曲线波分解的基本元素是脊波金字塔。曲线波变换的过程如下:
1)子带分解。首先定义子带滤波器组P0,二维目标函数f(x1,x2)经P0滤波后被分解为子带(△s,s≥0):f|→(P0,△1f,△2f,…)。
上的光滑窗ωQ(x1,x2)的集合中,用对应的窗函数乘上函数将产生局部化的Q。对每个尺度的Q执行同样的操作,也就是对所有s固定而k1,k2变化的Q=Q(s,k1,k2)。将窗分割算法应用于步骤(1)的每个子带:△sf|→(ωQ△sf)Q∈Qs。
3)归一化。设算子(TQf)(x1,x2)=2sf(2sx1-k1,2sx2-k2),把定义在Q上的输入部分变成定义在[0,1]2的输出部分。由步骤(2)得到的每个正方形再归一化为单位尺度
4)脊波分解。每个剖分块都在基元素为ρλ的离散正交脊波系统中进行分解,这是由ρλ产生L2(R2)上的正交基:
2.3 快速离散Curvelet变换实现
为了使曲线波变换更加容易使用和理解,并加快变换的速度以及减少冗余度,D.L.Donoho等人提出了快速离散曲线波变换(Fast Discrete Curvelet Transform,FDCT),即第2代曲线波。由于楔形窗不是标准的矩形,因此不能直接运用快速傅里叶算法实现,本文采用基于非均匀采样Fourier变换(Unequally-spaced Fast Fourier Transform,USFFT)的快速离散曲线波变换。
频率域的Curvelet变换如下
式中:j为尺度;l为角度;L1,j,L2,j为四边形Pj,l的边长;且
基于上述思想,可得快速离散曲线波算法如下:
1)对图像f(n1,n2)进行二维傅里叶变换,得到F(k,l),-n/2≤k,l≤n/2。
2)在F(k,l),中对每一个尺度角度(j,l)组合再取样(或插值),得到新取样函数:F(k,l,-ktanθl),(j,l)∈Pj。
3)用窗口Uj(k,l)乘以新的F(k,l,-ktan∈θl)得到λFj,l(k,l),(j,l)∈Pj。
4)对每个λFj,l(k,l)作二维傅里叶逆变换,得到Curvele系数c(j,l,k)。
以Lena图为源图像(图1a),快速离散曲波变换后系数显示图如图1b所示。
3 新型方向及尺度乘积的曲线波去噪算法
3.1 Curvelet硬阈值去噪算法假设含噪声图像模型为[4]
式中:I为原始图像;N为标准高斯白噪声。
去噪算法曲线波系数的硬阈值处理为
式中:是含噪声图像Curvelet系数;是经过硬阈值处理后的Curvelet系数;为子带系数的标准差;k为系数。
3.2 新型方向及尺度乘积的曲线波去噪算法
小波分解后各子带系数满足如下关系:边缘信号的特征和噪声所对应的小波系数都具有较大的值,但两者在尺度上的表现有很大不同。边缘信号的系数随着尺度传递下去,而噪声系数随着尺度会有较大的衰减,每层衰减约一半,因此小波系数在不同尺度上同样位置和方向具有相关性。将这个思想引入到Curvelet域中,同小波变换类似,Curvelet分解中不同尺度的系数也具有相同的特点[5,6]。在高频尺度,边缘信息的特征和噪声特征所对应的系数都有较大的值,但前者沿着尺度有较大的传承,而后者在传递过程中有较大的衰减[7]。基于上述思想,笔者的算法步骤如下:
1)对带噪声的原图像进行5层的Curvelet分解。
2)由于分解后的第3层系数是第4层系数的父系数,将分解后系数的第3层和第4层运用方向及尺度乘积的关系进行处理,处理方式如下:先判断系数是否满足
式中:Ri,j表示层相关系数,l=4表示尺度,s表示方向,T(l,s)表示l层s方向的阈值,如满足上式,则为边缘信号。将不满足式(9)的系数再用式(10)判断
式中:Ci,j(l,s)表示处理前的曲线波系数,α为调整系数。如果满足式(10)也判定为边缘信号。将剩余系数判定为噪声。
3)第1层为最粗尺度,不进行处理。
4)对第2和第5层采用改进阈值的方法处理。
曲线波传统硬阈值法去噪容易“过扼杀”变换系数,重构图像损失较多图像细节信息;用软阈值法去噪后的图像相对柔和,但往往残留不少轻微噪声,残留的噪声表现为长条状,严重影响图像的视觉效果。鉴于传统的硬阈值和软阈值方法去噪效果不够理想,笔者采用改进的阈值函数
式中:表示阈值处理后的曲线波系数,T(l,s)表示阈值,β是调节系数(0≤β≤1),可由实验确定,通过调节β,可以获得更好的去噪效果。
5)将处理后得到的5层系数进行Curvelet重构,得到去噪图像。
4 仿真及分析
为了验证本文算法的有效性,仿真采用512×512像素的Lena图像作为测试图像,叠加均值为零,标准差σ分别为10,20,25,30的高斯白噪声。图2a为含噪图像,图2b至图2d分别为用小波硬阈值、Curvelet变换阈值和笔者所提出的算法进行去噪后的图像,测试实验的去噪效果通过去噪图像的主观视觉效果和客观评价标准峰值信噪比(PSNR)两方面进行评价。
表1给出了不同方法去噪结果的对比。
从以上的去噪仿真可以看出,对于轮廓清晰、对比度明显的图像,所提算法具有更好的去噪效果,尤其在大噪声的情况下,也能得到比较清晰的重构图像。将本算法与小波变换进行比较,本算法的PSNR值提高很多,具有明显的效果,而且视觉上效果感觉明显地改善了。所提算法与Curvelet相比,由于Curvelet阈值方法在噪声比较大的时候,仅通过阈值无法区分边缘与噪声系数,而所提算法的PSNR值不但提高,而且视觉效果明显改善,其中的划痕现象明显减轻。
5 结论
笔者提出的新型方向及尺度乘积曲线波去噪算法,能将噪声信号比较正确地从边缘信号中区分出来,对于边缘信息较明显的图像,该算法比小波变换及Curvelet阈值去噪法具有更好的去噪效果,不但PSNR值得到提高,而且能明显减轻划痕现象的产生。要更进一步地消除划痕现象的产生,还需要结合其他图像处理技术,使主观视觉和客观评价指标达到最好的效果。Curvelet变换作为一种新的多尺度变换在图像去噪应用上显示出巨大的潜力,虽存在一些亟待解决的问题,但是由于其崭新的理论面貌和独到应用特点,Curvelet将在图像处理中发挥越来越大的作用。
摘要：小波对处理具有点奇异的目标函数是一种最优的表示方法。对于二维图像而言,图像的主要特征由具有线奇异的边缘所刻画。曲线波变换是一种多尺度变换,对于具有光滑曲线奇异性的目标函数,曲线波提供了稳定的、高效的和近于最优的表示。本文在第二代曲线波的基础上,提出了新型方向及尺度乘积的曲线波去噪算法。实验结果表明,本文算法明显优于小波图像去噪方法,也优于曲线波的阈值方法。
关键词：小波,多尺度,曲线波变换,图像去噪
参考文献
[1]冈萨雷斯.数字图像处理[M].北京:电子工业出版社,2007:128-132.
[2]焦李成,侯彪,王爽,等.图像多尺度几何分析理论与应用[M].西安:西安电子科技大学出版社,2008:101-132.
[3]李传珍,王晖,王静冬,等.基于BayesShrink阈值估计的Curvelet图像去噪[J].电视技术,2007,31(6):14-16.
[4]朱为,李国辉,涂丹.一种基于第2代曲波变换的尺度相关图像去噪方法[J].中国图象图形学报,2008,13(12):2282-2284.
[5]STARCK J L,CANDS E J,DONOHO D L.The curvelet transform forimage denoising[J].IEEE Trans Im.Proc.,2002(6):670-684.
[6]CANDS E J,DONOHO D L.Continuous curvelet transform:i.resolution of the wavefront Set[J].Applied and Computational Harmonic Analysis,2005,19(2):162-197.