精确算法

精确算法（通用7篇）

精确算法 篇1

随着21世纪信息技术的迅猛发展,微电子机械系统(MEMS)越来越广泛地用于光纤通信系统中,使得光纤通信器件不断集成化、小型化,大幅提升了器件的功能和可靠性。封装过程的自动化是降低光纤器件制造成本的关键,其中最重要的工序之一是光纤对准。

平面光波导器件的封装,是利用亚微米精度的活动平台,将波导芯片与阵列光纤进行对准耦合并固定,是目前光电子器件封装研究的热点问题之一[1]。目前国外技术可实现阵列波导器件的单通道功率损耗<0.3 dB,波导与阵列光纤的对准误差<1 μm。下一代阵列波导器件已提出功率损耗<0.15 dB的性能要求,平面波导与光纤的对准误差也要求<1 μm,对准技术面临更加严峻的挑战[2]。

精确对准过程中的位置搜索法是对准耦合的关键。本文就爬山法和模式搜索法进行仿真,通过仿真对这两种方法进行了比较,然后再确定可靠的搜索方式。

1 对准耦合原理

波导与光纤之间的连接损耗主要有对准偏差损耗(包括横向位错损耗、纵向间距损耗和轴向角度损耗)、模场匹配损耗和传输损耗等。前一项损耗取决于波导与光纤对准是否精确,后两项损耗取决于波导与光纤自身的结构。

1.1 对准耦合模型

当平面光波导与单根光纤对接时,常用的方法是将平面光波导固定,通过移动光纤实现对准。如图1所示,对准时共5个自由度,分别是横向位错x、y,纵向间距z,轴向角度α(俯仰角度)和β(横摆角度)。

1.2 对准耦合理论

波导与光纤的模场为高斯分布,如图2所示。波导与光纤的耦合效率由波导模场与光纤模场之间的重叠积分确定[3]。

考虑光波导的模场半径对称分布(ωx0=ωy0=ωg0),则对波导模场与光纤模场重叠部分积分可得平面光波导与光纤的耦合效率为

$\begin{array}{l}\eta =\eta {}_x\cdot \eta {}_y,(1)\\ \eta {}_x=k{}_x\exp \{-k{}_x[\frac{x{}^2}{2}\left(\frac{1}{\omega {}_{\text{g}0}^2}+\frac{1}{\omega {}_{\text{f}0}^2}\right)+\\ \text{π}{}^2\alpha {}^2(\omega {}_{\text{g}}^2(z)+\omega {}_{\text{f}0}^2)/2\lambda {}^2-\frac{x\cdot \alpha \cdot z}{\omega {}_{\text{g}0}^2}]\},(2)\\ k{}_x=4\omega {}_x^2\omega {}_{\text{f}0}^2/[(\omega {}_x^2+\omega {}_{\text{f}0}^2){}^2+\lambda {}^{2}z{}^2/\text{π}{}^2],(3)\\ \omega {}_{\text{g}}(z)=\omega {}_{\text{g}0}\sqrt{1+(\lambda z/\text{π}\omega {}_{\text{g}0}^2){}^2},(4)\end{array}$

式中,x为横向位错;z为纵向位错;ωg0、ωf0分别为z=0时光波导和光纤的模场半径;ωg(z)为光纤到光波导距离为z时的模场半径;λ为波长。将式(2)、式(3)中的x换作y,α换作β,即得ηy的计算公式。

$\text{功}\text{率}\text{损}\text{耗}(\text{d}\text{B})=-10\log {}_{10}({\rm P}{}_{\text{o}\text{u}\text{t}}/{\rm P}{}_{\text{i}\text{n}}),(5)$

式中,Pin为输入光功率;Pout为输出光功率。

2 算法原理

2.1 爬山法原理

有源自动对准的传统控制算法称为爬山法,该算法是当今发展最成熟、应用厂商最多的一种算法。

爬山法最早由美国Newport公司提出。该算法是一种局部择优的方法,是一维搜索算法,它采用启发式方法,每走一步都会有启发信息进行引导,确定下一步前进的方向以及前进的距离,也就是说在当前节点的时候,会比较各个方向与当前节点的高度差,然后选择与当前节点高度差最大的方向以一定的步长前进。当走到下一步的时候再进行选择,如此下去,一直走到山顶,找到最优解。如图3所示,从当前的节点开始,与邻节点的值进行比较。 如果当前节点是最大的,那么返回当前节点,作为最大值(即山峰最高点);反之就用最高的邻节点来替换当前节点,从而实现向山峰的高处攀爬的目的。如此循环直到达到最高点[4]。

如图4(a)所示,设定前进的步长,以一定的方向搜寻到峰值最高处,为了避开只找到局部的峰值(Local Peak),在峰值开始下降时,继续向前再搜寻N(checkpoints)个步长。如果在搜寻的N个步长中有比刚刚找到的峰值更高的光强度值,则判定刚刚找到峰值是局部峰值,需要再往下搜索,直到找到最大峰值;找到峰值后,步长减半,向反方向搜索,同样是在搜寻到峰值且开始下降时,再往前搜索N个步长,确定没有更高峰值后反向,如图 4(b)所示;然后如图4(c)所示,不停地重复以上两个步骤,便可找到一个最佳的峰值位置;最后如图4(d)所示,完成整个搜索。

2.2 模式搜索法原理

模式搜索法是由Hooke和Jeeves于1961年提出的,因此又称为Hooke-Jeeves方法[5]。与爬山法一样,模式搜索法不需要目标函数的导数信息,仅通过比较目标函数值的大小来移动迭代点,因此可用于光纤对准问题。与爬山法不同的是,它能以多个自由度同时进行寻优迭代过程,并通过模式移动保证迭代方向始终朝着“山峰”的走向逼近极大点,适合于多自由度搜索。

模式搜索法通过探测移动和模式移动完成对目标函数极大值的搜索,算法实现的关键在于确定基点、沿基点方向的移动以及后退策略[6]。如图5所示,初始点x(1)为一个基点,从该点开始沿X和Y方向进行探测性移动,得到第2个基点x(2),条件是x(2)处的功率值比x(1)位置更大,然后,沿着x(2)-x(1)的矢量方向向前跨越一步,即进行模式移动,步长是x(1)到x(2)的距离,随后从模式移动到达的位置x1再次开始探测移动,如此迭代,直到基点的功率值开始减小(图5中的x3),此时,这一次的模式移动和探测移动都是无效的,应该退回到前一基点x(4),减小步长,重新开始搜索,直到步长达到探索精度要求。

3 仿真实验

3.1 仿真参数

为了比较两种算法的差异,拟从以下两方面入手进行仿真实验比较:(1) 设定同一初始条件,比较两种算法的搜索速度和搜索步数。其中光纤的传输波长λ=1.33 μm,光纤的模场直径ωf0=4.6 μm,波导的模场直径ωg0=4.5 μm。(2) 两种算法各仿真10次,求平均值和标准差,比较两种算法的可靠性和搜索精度。

(1) 爬山法仿真参数定义

步长(stepsize=1 μm):在各个轴上,驱动器行走时需检测的最小长度;

检测点(checkpoint=4):为了防止陷入局部峰值,当爬山法搜索到峰值后要继续向前搜索的步数;

迭代次数(iteration=4):当爬山法搜索到峰值后,步长将要减半,反复搜索,更加细致地靠近峰值位置,重复次数即反复减小步长搜索的次数;

终止条件:达到设定的阈值或者已完成设定的重复次数,即可停止此次搜索。

(2) 模式搜索法仿真参数定义

缩减因子(α=0.5):当任何方向的探测都没有更大的功率值时,需缩减步长,小范围搜索峰值;

加速因子(β=1):用于移动方向时,方向探索正确的情况下加快搜索速度;

初始步长(stepsize=1 μm):初始时,每步移动的最小距离;

误差条件(eps1=1E-4):经过多次缩减步长后,步长减小到的最小值。

3.2 仿真实验

当初始位置为x=3 μm、y=2 μm、θx=0.1°、θy=0.5°、z=20 μm,目标效率为0.99时,爬山法和模式搜索法的仿真结果分别如图6和图7所示。

3.3 仿真结果分析

从自动耦合精密对准阶段的仿真结果对比(见表1)中不难发现:两种算法都能搜索到耦合效率>0.99的最大光强度值,其中模式搜索法的平均耦合效率和标准差都优于爬山法,此外,模式搜索法在搜索速度上也有一定优势。

4 结束语

爬山法是一维算法,多自由度搜索时需要在几个自由度之间反复循环,采样点数较多,搜索时间比较长,而且容易受到交叉耦合等各种噪声的影响从而导致搜索失败。

模式搜索法中探测移动的作用是确定基点和保证功率值增大的方向,而模式移动则沿着“山峰”使功率值增加得更快。由于该算法在X和Y方向同时探测,然后沿着两个基点矢量方向前进,所以能够克服两个自由度交叉耦合的影响。

仿真实验证明,模式搜索法可以更快地搜索到最大光功率值,算法收敛性明显高于爬山法。在光波导器件的封装制造中,阵列光纤精密对准是影响器件性能的决定性因素。从提高对准速度和算法稳定性出发,在光波导与阵列光纤对准中,运用模式搜索法可以更快地搜索到最大光功率值。

摘要：文章基于单模阵列光纤与平面光波导对准的理论模型,对在精确对准过程中常见的两种位置搜索法(爬山法和模式搜索法)进行了仿真模拟实验。仿真实验证明模式搜索法可以更快地搜索到最大光功率值,算法收敛性明显好于爬山法。

关键词：精确对准,爬山法,模式搜索法

参考文献

[1]Murakawa M,Nosato H,Higuchi T.Automatic opti-cal fiber alignment system using genetic algorithms[C].Marseilles:the 6th International Conference onArtificial Evolution,2003.129-140.

[2]沙慧军,陈抱雪.光波导-光纤耦合对接自动化系统的研究[J].光子学报,2005,34(12):1 773-1 777.

[3]Elsey J,Law S,Poladian L.Model based alignment ofoptical waveguides[C].New Orleans:53rd ElectronicComponents and Technology Conference,2003.1 565-1 567.

[4]周甘霖.TO-CAN封装镭射焊接机的寻光对位法则之研究[D].台湾:台湾省成功大学,1994.

[5]MOBARHAN K S,HAGENBUECHLE M,HEY-LER R.Fiber to waveguide alignment algorithm[J].Newport Application note,2000,(3):1-6.

[6]陈宝林.最优化理论与算法[M].北京:清华大学出版社,2000.

任意多面体重心及积分的精确算法 篇2

在工程应用和力学模拟中经常涉及到三维领域的多项式积分问题,比如有限元模拟中单元刚度矩阵(龙驭球,1988)[1]、非连续变形分析(DDA)和数值流形方法(NMM)中的惯性矩阵(Chen,1997;Shi,1984,1992,1996,2001;姜冬茹,2002;骆少明,2005;石根华,1997;郑榕明,2004)[2—10]。考虑到实际问题中涉及的积分区域多为复杂的多面体,通常的解决办法是把复杂形状的多面体划分成规则的几何单元后,采用数值积分方法求解;比如高斯积分方法(封建湖,2001)[1]。因此,能否用一种简单的方法直接从整体上求解复杂区域的多项式积分是值得研究的。石根华(1997)[9]给出了N维单纯形积分的公式,在此基础上林绍忠(2005)[12]给出了单纯形积分的递推公式。但两位学者给出的公式表述形式较为复杂,需要有具体的算法实现才有利于单纯形积分方法的推广。本文通过对经典单纯形积分公式(石根华,1997)的分析[9],给出了三维单纯形积分的伪代码,在此基础上给出了任意多面体重心及积分算法。利用上述算法,通过对中空规则形状和非规则形状积分区域下体积和重心理论值与计算值的对比分析,讨论了单纯形积分结果的精确性和普适性。同时,通过对比不同边长比(10-4—104)情况下计算结果的分析,讨论了积分区域图形条件对积分结果的影响。

1三维单纯形积分的基本原理及算法实现

石根华(1997)[9]给出了N维单纯形积分的通用公式,参考林绍忠(2005)[12]对单纯形积分公式的描述,公式(1)直接给出了三维单纯形积分表达式。

公式(1)中iw=ixw+iyw+iyw(w=0,1,2,3);J为雅克比行列式(石根华,1997)[9],在三维情况下有如下形式分别对应三维单纯形(四面体)的第i顶点的坐标。

从表现形式上看公式(1)比较复杂,致使一些数值流形方法的研究者选择其它积分方法,如姜冬茹等(2002)[7]、骆少明等(2003)[8]采用了Hammer积分方法,郑榕明等(2004)[10]采用了高斯积分方法。为了解决骆少明等(2003)[8]提到的“单纯形积分的形式十分复杂,不便于编制通用程序”的问题,表1给出了三维单纯形积分的伪代码实现,该代码经过简单变化并增加4阶行列式和阶乘函数即可方便转化成为C、C++或Fortran代码,进而实现多项式在三维单纯形上的积分运算。表1给出的算法仅对三维单纯形有效。要想实现任意形状积分区域的单纯形积分需要明确具体的积分流程。通过向石根华先生请教得到了任意形状多面体三维单纯形积分过程,结合图1凹多面体说明其流程。图1中多面体积分的过程如下:(1)把多面体每个面的顶点按外法线方向逆时针排列,即P6P5P4P3P2P1、P1P7P6、P6P7P5、P5P7P4、P4P7P3、P3P7P2和P2P7P1;(2)任意选一点P0,此点可以位于多面体内部或外部,如选在坐标原点,则可简化雅克比行列式J的计算;(3)以多面体的面为单位进行单纯形积分,比如针对面P6P5P4P3P2P1,分别以P0P1P2P3、P0P1P3P4、P0P1P4P5和P0P1P5P6四个单纯形进行积分,累加后完成该面的积分操作;(4)按(3)的步骤完成其他面的积分操作,全部累加后完成整个积分过程。

2中空规则形状单纯形积分的精度分析

根据表1的算法和图1的计算流程,我们可以计算出任意形状多面体的积分结果。为了分析积分结果的可靠性,下面通过对图2中空规则多面体的体积和重心的计算结果与理论结果对比分析进行具体说明。公式(2)给出了体积和重心的计算公式,式中V为多面体体积,X、Y和Z多面体重心的三维坐标。

针对中空多面体的三维单纯形积分问题(图2),对于外表面形成的多面体采用图1的流程正常积分,内表面形成的多面体每个面的顶点按该面的内法线方向逆时针排列并进行积分,最后累加即可得到最终结果。图2中外表面形成的正六面体的长、宽、高分别为100、50和100,内表面形成的六面体为边长为50的正方体。根据上述信息可以计算出两个六面体的的理论体积和重心,进而计算出填充部分的体积、重心X和Y。根据表1的算法和图1的计算流程,表2给出了图2中外表面形成的六面体、内表面形成的六面体和填充部分的体积和重心结果,其中V、X、Y和Z的含义与公式(2)一致。

表2结果表明理论结果与计算结果的差异很小,二者的差异出现在第16位有效数字。由于图2的规则图形的边长差异不大,下面通过调整边长比来研究图形条件对计算结果的影响(表3)。

表3结果表明图形条件对中空规则形状下三维单纯形积分结果的影响很小。由于只对Z坐标进行了尺度缩放,因此重心X、Y坐标应该保持不变,体积及重心Z随缩放尺度变化。表3中重心X、Y如果取14位有效数字,结果为V=375 000×scale、X=50, Y=175和Z=53.333 333 333 333×scale,相对误差约为10-15。另外,重心X值、Y值不随Z坐标尺度的缩放呈系统性变化,即当X和Y尺度远大于Z尺度或Z尺度远大于X和Y尺度时,重心X和Y均保持稳定。面积和重心Z的结果也具有上述特性。由此可见,图形条件对三维单纯形积分结果的影响不具有系统性,其差异主要是由计算机的计算误差引起的。因此,可以初步确定三维单纯形积分具有高精度的特性。

3 非规则形状单纯形积分的普适性分析

前面通过中空规则形状下单纯形积分结果理论值与计算值的对比分析,初步讨论了三维单纯形积分算法的有效性和精确性。为了分析该算法的普适性,下面对图3给出的非规则形状凹多面体进行体积和重心分析。由于图3的凹多面体的重心和体积没有理论结果可供对比,我们采用类似表3的分析思路来验证积分结果的可靠性。利用公式(2)及表1算法,可以方便的计算出图3凹多面体的体积和重心坐标,分别为V=76 727.166 666 666 7、重心X=49.503 344 208 809 1、重心Y=-1.961 582 381 729 2、 重心Z=28.775 000 000 000 0。

通过对图3的Y进行尺度缩放,表4给出了体积及重心的计算结果,结果表明在非规则积分区域上述参数均表现出高精度特性。由于仅对Y坐标进行了尺度缩放,因此重心X、Z应该保持稳定,体积V和重心Y应该随Y坐标的缩放尺度发生变化。如果取15位有效数字,结果为

V=76 727.166 666 666 7×scale,

X=49.503 344 208 809 1,

Y=-1.961 582 381 729 20×scale和

Z=28.775 000 000 000 0,相对误差约为10-15。另外,重心X值、Z值不随Y坐标尺度的缩放呈系统性变化,变化表现出随机特性。因此,上述结果表明即使积分区域为非规则的凹多面体,积分结果依然保持高稳定性和精确性,并且图形条件对积分结果基本没有影响。

4 讨论与结论

本文主要采用多面体体积和重心为参数对三维单纯形积分的精确性进行了分析,其中体积涉及0阶积分问题、重心涉及1阶积分问题,体积和重心的高精确性反映了积分结果的精确性。对于高阶积分问题,笔者也进行了相关研究(限于篇幅未在本文给出相应图表),结果表明高阶积分结果同样保持高精确性。

(1) 通过本文给出的三维单纯形积分算法和流程可以实现高阶多项式在任意复杂多面体上的积分。该积分算法可以直接应用于三维DDA(非连续变形分析)和三维NMM(数值流形方法)的程序设计中(Chen,1997;Shi,1984,1992,1996,2001)[2]。

(2) 本文给出的多面体的体积和重心算法可以精确求解任意复杂多面体的体积和重心问题。该算法适合凹多面体、中空多面体,具有高度的普适形和高精度性(相对误差为10-15量级)。

(3) 通过图形条件对三维单纯形积分结果影响的分析,结果表明即使尺度比从10-4变化到104,体积和重心结果依然维持稳定状态,误差分布具有随机特性,不依赖于尺度变化。由此表明结果误差主要来源于计算机的计算误差,算法本身是高度稳定的。

致谢:本文的研究得到了石根华先生的悉心指导,在此表示感谢。

参考文献

[1]龙驭球.有限元法概论.北京:人民教育出版社,1988

[2] Chen G Q,Ohnishi Y,Ito T.Development of high order manifoldmethod.ICADD-2,Kyoto,Japan,1997

[3] Shi G H,Goodman.Discontinuous deformation analysis.In:Proceed-ings of the 25th U.S.Symposium on Rock Mechanics,Evanston,25-27 June 1984:269—277

[4] Shi G H.Discontinuous deformation analysis:a new numerical modelfor the statistics and dynamics of deformable block structures.Engi-neering Computations,1992;9(2):157—168

[5] Shi G.H.Manifold method.Proc of the 1st 1nt.Forum on DDASimulation of Discontinuous Media.Bekerley:[S.n.],1996:52—204

[6] Shi G H.Three dimensional discontinuous deformation analysis.In:Rock Mechanics in the National Interest.Elsworth,Tinucci and Hea-sley,2001:1421—1428

[7]姜冬茹,骆少明.三维数值流形方法及其积分区域的确定算法.汕头大学学报(自然科学版),2002;17(3):29—36

[8]骆少明,张湘伟,吕文阁,等.三维数值流形方法的理论研究.应用数学和力学,2005;26(9):1027—1032

[9]石根华.数值流形方法与非连续变形分析.装觉民,译.北京:清华大学出版社,1997

[10]郑榕明,张勇慧.基于六面体覆盖的三维数值流形方法的理论探讨与应用.岩石力学与工程学报,2004;23(10):1745—1754

[11]封建湖,车刚明,聂玉峰.数值分析原理.北京:科学出版社,2001

精确算法 篇3

随着我国电源与负荷的快速增长,直流输电在我国有了飞跃发展。由于高压直流输电换流器的非线性特性,导致在换流站的直流侧和交流侧产生谐波电压和谐波电流。这些谐波严重影响输电效率,影响设备寿命和安全,同时线路上的谐波也会对周边通信产生干扰,对直流系统的谐波治理主要是加装滤波器。

在直流系统滤波器设计时,准确分析直流系统的谐波成份,是合理科学配置滤波器参数的基础。同时,对于有源滤波器的运行,需要实时获得直流系统相关谐波参数。然而由于实际运行中,直流系统的谐波是不断变化的,这给谐波的动态测量增加了很大难度。因此,谐波的监测仍没有得到很好的改进[1]。傅里叶分解算法需要一定时间的谐波参数,且需进行两次变换,计算量大,计算时间长,从而使得检测时间较长,检测结果实时性较差,在频域分析和降噪方面,其精度不是太高。近来许多神经网络模型被推荐用于电力系统谐波估算分析,文献[2]中推荐ANN算法用于电力谐波分析,该算法采用ADALINE模型和LMS学习机制,该方法可以满足一些用相对较少设备进行的测量工作。文献[3]中推荐了用加权多点插值方法改善多频率信号的幅值估计。

Prony算法在现代频谱分析中发挥了重要的作用,然而,将其用于电力系统信号处理和分析的时间还不长。本文详细介绍了基于扩展Prony算法的谐波分析方法,该方法解决了快速傅里叶分解算法在计算谐波幅值和相角方面的弱点。

1 基于扩展Prony算法谐波分析模型

1.1 Prony算法描述

Prony算法是由法国数学家提出的,用指数函数的线性组合来描述间距采样数据的数据模型[4,5]。形如

或者

假设f(n)由等时间间距Δt的N点组成,从到。由式(1)有以下一系列等式成立:

在式(1)中:ci和μi都是复数,ci=Aiexp(jθi),μi=exp[(αi+j2πfiΔt)]。其中Ai是幅值;θi是相角;αi是衰减因子;fi是频率;Δt是采样间隔;N是采样点数;f(^n)表示f(n)的估计值。将式(2)转换成矩阵形式如下:

其中:

1.2 Prony模型参数求解

构造一个价值函数如下:

为了求得参数iA,θi,αi,fi,需要求解非线性最小二乘方程式。方程组(2)是常系数线性差分方程均匀解,这也是Prony方法的关键。因此,采用所测数据通过求解差分方程可以得到相关系数。为了推导出线性差分方程,定义一个特征多项式如下:

其中:α0=1。

从方程组(2),可以得到

将上述等式两边同时乘以αk并求和可得

考虑μin-k=μin-pμip-k,可得

式(5)第二个和项恰好是特征多项式(4)根μ的位置,因此式(5)为0,从式(5)可知,f^(n-k)满足如下等式:

定义σ(n)为f^(n)的实测与估计值之间的偏差,则

将式(6)代入式(7)可得

显然,该模型是特殊的序列协方差矩阵求解过程,其含有白噪声并带有相同的AR和MA参数。为了避免求解非线性方程式,将偏差定义如下:

式(8)可以变换为

显然,式(9)是个AR模型,采用整体最小二乘方法可以求得使偏差最小的相应系数a1,…,ap。从而可以求得特征多项方程式代入式(3),最小二乘的解是

最后,可以得到谐波的参数如下

2 仿真实例与结果分析

为了验证Prony算法在直流谐波测量方面的有效性和精确性,本文采用PSCAD自带的直流输电模型作为测试对象,该模型采用单12脉动换流器作为交直流转换装置,用T型线路模型来近似表示直流线路。该模型中不含直流滤波器,这样更易于测试直流线路谐波特性。具体电路图如图1。

图2给出了仿真结果直流线路谐波录波图,表1给出了分别采用FFT算法和扩展Prony算法对直流线路谐波测量结果对比表。由于直流线路谐波主要是3倍次谐波,表中只给出了3倍次谐波的对比结果。由于本文所采用的直流输电模型是简化的直流输电模型,未考虑背景谐波等不确定因素的影响,表1中所给出的直流线路谐波电流理论幅值和相角是通过理论公式计算所得[6,7]。

由表1的对比结果可知,扩展Prony算法在直流谐波测量上精度高于傅里叶分解算法,尤其在相角的估计测量上,其精度高的优势更加明显。

3 结论

扩展Prony算法以最小二乘法为基础,采用指数矩阵线性组合可以直接求解信号幅值、相角、阻尼比和频率,具有较高的实时性和精确度。本文通过理论分析和仿真测试对比,证明了扩展Prony算法应用在直流谐波测量方面是有效的,并且具有更高的精度。扩展Prony算法的这一特点为直流输电谐波分析提供了更好的测量工具。

摘要：介绍了一种更精确的,适用于直流输电系统谐波分析的扩展Prony算法,该算法应用指数方程线性组合估算信号,给出了其在直流谐波分析上的应用。同时,基于PSCAD平台,利用其自定义模型功能构造了采用扩展Prony算法的谐波分析元件,通过对直流系统谐波仿真测试并与PSCAD自带傅里叶算法对比,证明了该方法的有效性和精确性。

关键词：直流系统,谐波分析,Prony,PSCAD,FFT

参考文献

[1]Arrillaga J,Bradley D A,Bodger P S.Power System Harmonics[M].New York:John Wiley&s Sons,1985.

[2]金明,刘远龙.用于电力系统谐波分析的ANN算法[J].电网技术,1997,21(5):52-54.JIN Ming,LIU Yuan-long.A New ANN Harmonic Analysis Algorithm For Power System[J].Power System Technology,1997,21(5):52-54.

[3]Dusan Agrez.Weighted Multipoint Interpolated DFT to Improve Amplitude Estimation of Multi-frequency Signal[J].IEEE Trans on Instrumentation and Measurement,2002,51(2):287-291.

[4]Hildedrand F B.Introduction to Numerical Analysis[M].New York:McGraw-Hill,1956.

[5]张贤达.现代信号处理[M].北京:清华大学出版社,1995.ZHANG Xian-da.Modern Signal Processing[M].Beijing:Tsinghua University Press,1995.

[6]王峰,游广增,黄莹,等.高压直流输电3脉动谐波电压源计算[J].高电压技术,2009,35(10):264-268.WANG Feng,YOU Guang-zeng,HUANG Ying,et al.Calculation of Three-pulse Harmonic Voltage Source for HVDC Transmission[J].High Voltage Engineering,2009,35(10):264-268.

精确算法 篇4

关键词：超宽带,信道模型,室内定位,多层神经网络,TDOA,AOA

无线通信技术[1]的发展使得对于无线定位的需求日趋增加,传统的室外定位技术已被广泛应用,人们对复杂室内环境精确定位的需求越来越强烈,比如消防、医疗等。超宽带(Ultra - WideBand,UWB)[2]因为是一种无载波通信技术,以及GHz量级的带宽,所以可有效地降低发射信号功率、系统复杂度的同时提供厘米级的定位精度。

根据目前的研究状况,超宽带室内精确定位算法可以分为5 类:基于信号到达时间(Direction of Arrival,TOA) 的定位算法[3]、基于信号到达时间差( Time Difference of Arrival,TDOA) 的定位算法[4]、基于信号到达角度(Angle of Attack,AOA)的定位算法[5]、基于接收信号强度( Received Signal Strength Indication,RSSI)的定位算法以及混合定位算法[6]。考虑到实现的可能性以及复杂程度,基于信号到达时间(Time of Arrival,TOA)的定位算法被排除,同时考虑到定位的精度,基于接收信号强度(RSSI) 的定位算法也被排除。

随着对定位精度要求的不断提高,单一定位方法无法满足需求,融合两种及以上的混合定位方法就成了一个发展趋势。文献[7]讨论了TDOA/AOA混合定位方法,获得了比TDOA、AOA单一定位算法更高的定位精度。因此,本文提出了一种基于多层神经网络的超宽带混合定位算法,利用多层神经网络对到达时间差(TDOA)和到达角度(AOA)测量值中的非视距传播误差进行修正,来提高超宽带的定位精度,与传统TDOA / AOA定位算法相比,定位精度有所提高。

1 仿真信道模型

工作组提出了3 不同的模型,本文研究的方向为室内环境下的超宽带精确定位,基于这个目的选择2 ~ 10 GHz的超宽带信道模型作为仿真的信道模型,着重研究住宅环境。2 ~10 GHz的超宽带信道模型是建立在路径增益模型(Pathgain)、功率延时分布模型(Power Delay Profile)和小尺度衰落模型(Small - scale Fading)的基础上。

1. 1 路径增益模型

路径增益[8,9]与频率的关系函数为,k是频率对于路径增益的影响因子。

路径增益与距离的关系函数为,d0是基准距离,设定为1 m,G0是距离d0=1 m处的路径增益,n的值取决于实际环境是视距(Line-of-Sight,LOS)还是非视距(Nonline-of-Sight,NLOS)。可得到路径增益、频率及距离之间的函数为

式中,ηRX - ant(f)是接收天线的增益;ηTX - ant(f)是发射天线的增益。

如果考虑阴影(Shadow) 的影响,S为服从均值为0;方差为 σ 的高斯分布的随机变量,不同环境表示为不同的具体数值。

1. 2 功率延时分布

本文主要研究室内环境下的超宽带精确定位,因此只考虑信号到达率(Ray Arrival Rates)。信号到达率由综合两个泊松过程得到

τ0,l= 0,β 为混合概率,λ1,λ2为到达率。

IEEE802. 15 工作组借鉴了S - V模型来拟合信号的分簇到达时间(Cluster Arrival Times),可用冲击响应hdiscr(t)来描述

Øk,l为相对相位,在一个带通系统中其服从[0,2π]中均匀分布,ak,l为第l个簇,第k个分量的加权系数,Tl是第一簇的到达时间,τk,l是第l个簇的k个分量相对于Tl的延时。

由于簇的数量L是满足概率密度函数的泊松分布,因此均值可以用来描述概率密度函数

泊松分布函数 ρ(TlTl - 1) = ∧lexp[- ∧l(Tl-Tl - 1)],l >0,∧l第l簇的到达率。

1. 3 小尺度衰落(Small - scale Fading)

IEEE802. 15. 4a的小尺度衰落采用Nakajami信道模型

式中,m≥1/2是Nakagami信道的m因子;Γ(m)为伽马函数;Ω 是平均能量。

2 超宽带TDOA / AOA混合算法的改进

2. 1 TDOA / AOA定位原理

TDOA的几何原理是双曲线测量定位,AOA的几何原理是三角测量定位,二维空间原理如图1 和图2所示。

如图1 所示(x1,y1)、(x2,y2)、(x3,y3)为3 个已知节点位置,则待测节点(xm,ym)与已知节点的位置关系如下

如图2 所示(x1,y1)和(x2,y2)为两个已知节点位置,角度值 α1、α2为通过定向天线测得的相对角度值。则待测节点(xm,ym)与已知节点的位置关系如下

2. 2 改进的TDOA / AOA混合算法

研究需要获取三维位置信息,根据上述TDOA/AOA定位原理,定位基站至少需要4 个,基站数量增加可以提高定位精度,但相应的成本也随之增加,目前市场上每个基站的售价约5 000 元,以下是在4 个基站的情况下建立的多层神经网络模型[10,11,12]。

(1)输入层由参与定位的4 个基站提供的3 个TDOA测量值和4 个AOA测量值所组成,因此,输入向量可表示为

(2)隐含层的层数、隐含层各层神经元数目以及各层激励函数的选择,目前还没有理论上的指导。隐含层的层数、神经元数目增加可以提高定位精度,但计算的复杂性也随之增加。根据仿真数据如图4 所示,隐含层的层数选为3 层。根据基站布设模型,隐含层各层神经元数目选为9。

各层激励函数

(3)输出层由3 个神经元构成,输出向量O[X,Y,Z]。基于神经网络的定位算法具体步骤: 1 ) 利用Matlab神经网络工具包,建立本文提出的神经网络模型。2)训练神经网络模型,如图5 所示。

第一层,定位范围1. 8 m × 1. 8 m,每200 组数据训练一个神经元,一共训练9 个,各神经元间距<0. 5 m。

第二层,定位范围0. 9 m × 0. 9 m,每150 组数据修正一个神经元,将第一层的9 个神经元全部修正一次,各神经元间距<0. 4 m。

第三层,定位范围0. 4 m ×0. 4 m,每100 组数据修正一个神经元,将第二层的9 个神经元全部修正一次,各神经元间距<0. 3 m。

3 仿真与分析

为检验该算法的实际性能,用Matlab软件对本文定位算法进行仿真,对本文算法和单纯采用AOA、TDOA、TDOA / AOA算法进行了比较,仿真结果如图6所示。

如图6 所示,在相同信噪比环境下,利用多层神经网络的TDOA/AOA算法能够明显降低均方根误差。如果当信噪比降低时,采用多层神经网络的TDOA/AOA算法比原TDOA / AOA算法效果更明显。

4 结束语

精确算法 篇5

在物流配送中, 配送路径的规划对行业的利益有着重大影响。针对于最优路径规划问题有着很多算法可以选择, 如模拟退火算法、蚁群算法、禁忌搜索算法等。而在这些算法中, Dijkstra算法随着整个配送路径网络的增大, 算法内部的二重循环将使其执行效率严重下降。模拟退火算法能够获取优秀的解决方案, 但是冷却参数的设置难以把握。蚁群算法虽然具有很强的实用性, 但易于陷入局部最优解并且收敛速度慢。禁忌搜索算法通俗易懂、易于实现、具有较强的局部搜索能力, 且能够避免陷于局部最优。

在现有的一些对禁忌搜索算法改进的优化方案中, 如配送中K辆配送车辆服务M个客户, 每车均有最大行驶距离, 当有车辆达到最大配送路径, 则由新车辆加入配送服务[1]、配送中心拥有同种车辆多辆, 有车辆载重限制[2]、限定车辆装载质量为Q, 每辆车辆一次配送的最大行驶距离一定[3], 都是在配送路径规划问题上, 将整体路径规划划分为子问题 (多辆车与多用户的配送路径规划问题) , 如根据配送车辆的载货量, 分配各配送车辆的责任范围, 协调配送车辆的调度。但是在子问题上 (每辆配送车辆的路径规划) 仍是使用禁忌搜索算法。即是这些方案都是对配送任务做出了细化, 为每辆配送车辆指定任务量, 通过算法协调车辆调度, 以达到优化配送的目的, 并未对禁忌搜索算法做出优化。而像利用遗传算法变异算子能产生具有新特性的个体的优点, 以遗传算法与禁忌搜索相结合, 利用两种算法的特性互补, 提高了算法的优化性[4], 但无疑增大了算法的难度与计算强度。

本文针对物流配送路径规划问题为禁忌搜索算法设计了伪精确优化算法, 在给定优化强度参数 (优化强度可调节) 的情况下, 对禁忌搜索算法的结果进行再次优化以达到获取最优路径的目的。

禁忌搜索算法的原理与实现步骤

禁忌搜索 (Taboo Search) 是人工智能的一种体现, 是局部领域搜索的一种扩展。禁忌搜索最重要的思想是标记对应已搜索的局部最优解的一些对象, 并在进一步的迭代搜索中尽量的避开这些对象 (并不是绝对的禁止循环) , 从而保证对不同的有效搜索途径的探索。禁忌搜索算法通用、易于实现。

算法流程[5]:

(1) 给定算法参数, 随机产生初始解x, 置禁忌表为空。

(2) 判断算法终止条件是否满足。若是, 则结束算法并输出优化结果;否则, 继续一下步骤。

(3) 利用当前解的领域函数产生其所有 (或若干) 邻域解, 并从中确定若干候选解。

(4) 对候选解判断藐视准则是否满足。若成立, 则用满足藐视准则的最佳状态y替代x成为新的当前解, 并以与y对应的禁忌对象替换最早进入禁忌表的对象, 以y替换最优状态, 然后转步骤 (2) , 否则, 继续以下步骤。

(5) 判断候选解对应的各对象的禁忌属性, 选择候选解集中非禁忌对象对应的最佳状态为新的当前解, 同时以与之对应的禁忌对象替换最早进入禁忌表的禁忌对象元素。

具体描述如下:

输入P为有序数组, 禁忌表D, 优化强度S。

1.初始化D, 循环次数S=0, 令P′=P, 计算P的权值W (P) 。

2.当S≤s时, 输出P, 退出程序。

4.当d<f (s) 时, 按m对随机数交换P′中的数, 将其优化后得到新数列P″, 计算W (P″) , 否则返回第二步。f (s) 为禁忌限制强度函数。

5.当W (P) >W (P″) , P=P″, W (P) =W (P″) , s=0, 初始化D;否则s=s+1, 更新D:

一、伪精确优化算法设计原理与实现步骤

伪精确优化是针对在禁忌搜索算法给定优化强度的情况下, 禁忌搜索算法的响应速度相应的提高, 但是其结果相对的需要进一步地优化来保证与提高。伪精确优化算法对优化结果进行再次优化而设计的算法。伪精确优化算法由移动和翻转两个部分组成:

二、移动

简述:在当前解的路径点序列中, 随机移动其中任意一点的位置, 计算移动后的序列路径权值。若新序列的权值低于旧路径的权值, 则接受新路径, 并更新当前路径权值。

具体描述如下:

输入完全图G (V, E) , P为V的有序数组。

1.初始化, 将P赋值于P1, 计算P及P1的权值:

2.选取P1中任意一位置i上的点, 再任意选取中任意一位置j, 且|i-j|>i;将i上的点移动到j, 其中间的数列 (含j) 水平向i平移一格, 得到P1′, 并计算W (P1′)

4.返回第二步, 遍历P1中每一对相距大于1的位置。

5.当W (P) >W (P1) , 则P=P1, W (P) =W (P1) , 并返回第四步, 否则输出P。

三、翻转

简述:在当前解的路径序列中, 随机翻转其中一段子序列, 计算翻转后的新路径序列的权值。若新路径权值低于旧路径权值, 则接受新路径, 并更新路径权值。

具体描述如下:

输入完全图G (V, E) , P为V的有序数组。

1.初始化, 将P赋值于P1, 计算P及P1的权值:

2.选取P1中任意两位置i, j且i<j。将i与j之间的数列颠倒 (含i, j) , 得到P1′。

4.返回第二步, 遍历P1中每一对的位置。

5.当W (P) >W (P1) , 则P=P1, W (P) =W (P1) , 并返回第四步, 否则输出P。

算法流程简述:

(1) 以最短路径算法生成初始解, 初始化禁忌表, 并用伪精确优化算法对初始解进行首次优化。

(2) 判断是否达到优化强度, 若是, 则输出最优结果;否则, 继续以下步骤。

(3) 产生随机数对组M。

(4) P经过M变化为P’

(5) 判断是否达到循环条件, 若是, 则输出优化结果并继续后续步骤;否则转到 (2) 。

(6) 针对输出的当前优化结果进行伪精确优化, 转到 (2) 。

五、测试结果

1.生成配送路线

测试内容:生成路径所需的响应时间。

测试条件:仓库 (一个) 和客户 (多个) 的地址信息, 计算强度。

输出:优化配送路径。

测试结果:

从表1可以看出, 对于选址数小于等于80的情况, 计算强度3比较好。对于选址数位100的情况, 虽然计算强度大于3可以将总里程优化将近1%, 但却将平均响应时间增大了1~2倍。对于大型物流企业, 需要在总体运输规划时间成本和最终路径运输成本之间进行仔细权衡, 应避免规划时间过长而造成系统不可用。

2.动态规划路线

测试内容:动态规划路线响应时间 (路径生成, 随意改变路况信息生成的配送路线) 。

测试条件:路况更改信息, 计算强度1。

输出:更改路况后的后续运输点之间的路径。

测试结果:

六、结论

本文在禁忌搜索算法的基础上设计了物流配送中的最优路径规划优化算法, 在给定优化强度的条件下提高了算法响应速度, 并通过伪精确优化算法实际的对禁忌搜索的优化能力进行了提升。算法整体继承了禁忌搜索算法通用、易于实现, 且容易理解的特性, 比使用更加复杂、更加高级的智能算法而言, 设计结果更具有普适性, 并且算法复杂度为O (nk) 。本文设计出的算法亦可用于以问题着手进行算法优化设计的方案, 解决每一条的配送路径规划。

摘要：本文针对物流配送中最优路径规划问题设计出了一种结合禁忌搜索和伪精确优化的路径规划算法, 本算法规划出的配送路径具有比较好的整体优化能力与快速响应能力。算法经物流配送模拟软件实际运行验证效果较好。

关键词：物流配送,路径规划,禁忌搜索,伪精确优化

参考文献

[1]戴东, 王江晴.一种基于禁忌搜索算法的车辆路径问题的改进算法[J].中南民族大学学报, 2007, 26 (1) :64-66.

[2]陈晓眯, 孟志青, 徐杰.基于混合禁忌搜索算法的动态车辆路径研究[J].浙江工业大学学报, 2009, 37 (5) :580-585.

[3]李松, 李瑞彩, 刘兴.基于改进禁忌搜索算法的车辆路径优化[J].铁道运输与经济, 2008, (5) :91-94.

[4]余丽, 陆锋, 杨林.交通网络旅行商路径优化的遗传禁忌搜索算法[J].测绘学报, 2014, 44 (11) :1197-1203.

精确算法 篇6

在解决机器人的路径规划问题中,通常用多面体模型来模拟机器人与障碍物,多面体又分为凸多面体和凹多面体两种。那么,机器人与障碍物之间的碰撞检测情况就转换为检测多个几何多面体之间的碰撞情况[1]。目前,在三维空间中,科学技术人员已经提出了计算两个多面体的Minkowski和的多种不同方法,主要都是以计算Minkowski和的边界值为目标,并利用一些方法来表示它[2]。其中,凸多面体的Minkowski和的计算方法已经比较成熟,而凹多面体的Minkowski和只能获得近似值,用抽取等值面的方法来高度近似Minkowski和的边界值[3]。

移动立方体MC(Marching Cubes)算法[4]作为一种有效的构造等值面的算法,使用三角面片作为中间几何图元的基本表达元素。MC算法本身原理简单 、容易实现 ,因此得到了广泛的应用 ,被认为是目前最为流行的等值面抽取算法之一[5]。对科学计算可视化领域的研究产生了深远的影响。然而,算法也存在一些不足之处,主要集中在提取曲面正确的拓扑结构、提高绘制精度、减少MC存储空间和提高处理速度等[6],本文则主要针对减少储存空间和提高处理速度进行改进。

1 GPU加速等值面生成

MC算法是一种基于立方体的等值面重建算法,基本思想是逐个处理数据场中的体素,分类出与等值面相交的体素,采用插值计算出等值面与体素棱边的交点。根据体素中每一顶点与等值面的相对位置,将等值面与立方体边的交点按一定方式连接生成等值面,作为等值面在该立方体内的一个逼近表示。尽管算法本身能够很好地实现近似,但其中还是存在着一些问题,比如在构造等值面时,是用三角形面片拟合的方式来实现的,由于形成的三维数据场密度是任意的,构造的立方体大,势必会造成很大的误差,影响最终的近似效果;构造的立方体小,近似得比较精确,但三维数据场密度就会很大,计算量极大,响应速度相应变慢。为了保证计算的实时性和高精度,针对这些问题,在传统MC算法的基础上提出了一种改进的MC算法——基于GPU的快速MC算法。

GPU是目前计算机上普遍采用的图形图像处理专用器件,具有强大的并行处理能力和高效的数据传输能力,是一个功能强大的并行计算单元。科研人员用它来加速科学计算和可视化应用程序,取得了令人鼓舞的研究成果[7]。

由MC算法可知,近似的等值面是由小的三角面片组成的,小的三角面片是从体素中提取出来的,且只需要知道一个体素八个定点就可以确定,用MC算法来近似等值面的计算具有很强的可并行性。

1.1 白立方体的排除

定义1 白立方体就是将三维多面体模型体素化后形成三维数据场,三维数据场由数以万计的立方体组成,其中与三维多面没有交点的立方体叫白立方体即白体素,其他的叫作边界体素,等值面就在这些边界体素中存在。2002年徐晓玲等提出了区域增长边界体素搜索策略。区域增长边界体素搜索算法可将三维数据场中的白立方体排除,保留下来的体素必是图1中的一种,其中0号是白立方体,这样可避免对整个体素空间的搜索,大大减少后续工作的计算量。

1.2 三次线性插值确定交点

把白立方体排除之后,剩余的立方体体素都是与等值面有交点的,并且这些交点出现在立方体体素的棱边上,现有的确定交点的方法中,线性插值比较具有代表性,也有为了实时性用中点来代替的。那么为了能又快又不影响精度,这里提出了三次线性插值的方法来确定交点,将误差降低到棱边长度的0.125倍以下。如图2所示。

M是Vi和Vj的中点,P是M和Vi的中点,Q是M和Vj的中点。P1,P2,P3,P4分别是中点。采样点与其中哪个点的距离比较小就用那个点来代替。计算快,而且提高重建精度。

1.3 将边界立方体数据映射到纹理单元

一个立方体数据要存储八个顶点坐标和立方体连接关系两部分,按顶点顺序将所有的顶点坐标存储下来,并将顺序号作为顶点索引,然后以顶点索引的方式存储立方体中顶点的连接关系,在实际操作中,用三个二维纹理来实现,第一个纹理采用32位浮点RGBA格式存储顶点坐标和函数值,一个顶点对应一个纹理但愿,顶点的三维坐标和函数值分别和R、G、B、A四个分量一一对应。第二个纹理采用8个8位的RGB格式存储顶点对应的索引信息,方便查找顶点坐标。第三个纹理单元同样采用8个8位的RGB格式存储边表的对应关系。

1.4 基于GPU加速的MC

当输入数据映射到纹理并设计好输出数据结构以后,就可以在片段着色程序中实现等值面提取算法了,顶点的函数值vi 可以根据索引在顶点纹理中查询得到,再按照阈值F可以将顶点进行分类,阈值的计算方法为:其中t1,t2,…,tn,是我们选择的边缘点,f(t1),f(t2) ,…,f(tn)是各点的灰度值,则:

$F=\frac{{\displaystyle \sum _{i=1}^{n}f}(t{}_i)}{n}$

αi表示掩码的第i位,根据比特掩码和边表就可决定立方体与等值面相交的边,利用线性插值方法即可算出顶点,连接即是我们所求的等值面。

$\alpha {}_i=\{\begin{array}{cc}0& (v{}_i<F)\\ 1& (v{}_i\ge F)\end{array}i=(0,1,\cdots ,7)$

利用GPU的并行浮点运算的性能,实现实时提取等值面[8,9,10]。如下所示是基于GPU加速的MC实时等值面的算法如图3所示。

2 能量函数优化等值面

2.1 能量函数公式

等值面生成技术所产生三角网格中的三角面片数目通常数以万计,且三角面片均是平面图形。连接而成的等值面不够圆滑,不够接近实物,为了能更好地模拟实际模型,提出一种能量函数[11,12],在保持原有网格结构不变得前提下,对等值面的三角网格进行优化。

定义的能量函数如式(1)所示:

$\begin{array}{l}E({\rm M})={\displaystyle \sum _{i=1}^{n}d}{}^2(p{}_i,\text{ϕ}{}_i(\left|{\rm K}\right|))+{\displaystyle \sum k}d{}^2(v{}_i-v{}_j)+\\ E{}_{scalar}({\rm M})+E{}_{disc}({\rm M})(1)\end{array}$

其中M表示三角网格,P={p1,p2,…,pi}表示点集,$d{}^2(p{}_i,\text{ϕ}{}_i(\left|k\right|))$表示点集P到网格M距离的平方和,d2(vi-vj)表示网格所有边长度的平方和,其中k表示边的弹性系数,vi表示网格的顶点。

2.2 寻找分裂节点

上述算法得到的三维网格模型是由三角面片组成的,三角面片的大小不一,对于其中面积较大的三角面片分裂成更小的三角面片,能使相似度大大得到提高,那么寻找分裂节点将成为一个关键,用三角面片的三个顶点坐标来表示,计算算法比较复杂,我们用三角面片的重心T来做分裂节点,计算方法比较简单,一个网格的重心坐标向量可以通过下述方法得到,假设,S表示三角形面积,Vi,Vj,Vk表示三角面片的三个顶点的坐标。

bi=S△TVjVk/S△ViVjVk

bj=S△TViVk/S△ViVjVk

bk=S△TViVj/S△ViVjVk

则T的坐标可表示为:

T=biVi+bjVj+bkVk (2)

则T的重心坐标向量为:

b=biei+bjej+bkek (3)

2.3 网格优化

通过线性最小二乘数值算法使能量函数式(1)的值达到最小时,网格能得到很好的优化。对于${\displaystyle \sum _{i=1}^{n}d}{}^2(p{}_i,\text{ϕ}(\left|k\right|))$中,每个Pi所对应的$\text{ϕ}(\left|k\right|)$都可以用${\displaystyle \sum _{i=1}^{n}d}{}^2(p{}_i,\text{ϕ}(\left|k\right|))$来表示。式(1)中的所有变量均可用顶点坐标V(x,y,z)来表示,那么就会得到:

E(M)=E1(x,y,z) (4)

欲使式(4)的值达到最小,根据极值理论可知,只有在E(M)=E1(x,y,z)的导数为零的时候才能取得,如下:

这样,根据求每个变量的偏导数,即可求得每个网格顶点坐标的最小值,这样就得到使用能量函数优化网格后所有顶点的新值,从而得到网格被优化后各顶点的新位置。

3 实验分析

本文的算法是在酷睿双核2.4GHz、内存4GB、VC++6.0和OpenGL环境下实现,以重建一个兔子的模型为例,对重建后的兔子模型用能量函数进行优化,表1给出了加速前后重建时间的区别,图4为重建后兔子的模型,图5为用能量函数优化后的兔子网格结构模型,优化前后有关数据的分析如表2所示。

通过上述实验数据可以看出,利用GPU加速后,重建速度明显提高,尤其是在数据量比较大的情况下,效果更加明显。优化前的三角网格中,三角面片的相对面积大小很不均匀,经过能量函数的优化,大部分三角面片的相对面积大小比较接近,兔子网格的空间形状进一步得到了改善,使重建的模型更逼近实际形状。使本文的算法得到了很好的验证。

4 结 语

利用GPU的浮点运算能力和并行运算能力,提高了MC算法的重建速度,与CPU相对,GPU的并行运行运算性能突出,尤其是当需要处理的数据量很大时,提升的效果更加明显。为了能更好地接近实物,运用能量函数优化生成的三角网格,在不改变原有结构的同时,使重建的模型与实物的相似度到达最大。但是实际需求对实时性和精确度要求越来越高,还有待更好的方法来解决问题。

摘要：Minkowski和算法能够精确地检测出若干子凸多面体之间是否发生碰撞,而凹多面体Minkowski和的边界只能用等值面来近似表示,等值面的生成通常需要计算数以百万计的体素,计算量比较大,耗时比较长,利用GPU的高速并行运算能力和浮点运算能力,加速等值面的形成,以达到实时生成等值面的效果。提出一种能量函数,在不改变原有结构的同时优化等值面,使其更接近实际形状,从而实现凹多面体实时的和高精度的碰撞检测。实验表明,提出的方法提高了等值面的绘制速度和绘制精度,从而提高了碰撞检测的精确性,尤其是在数据量比较大的情况下,效果更明显。

关键词：Minkowski和,移动立方体,等值面,能量函数

参考文献

[1]Guo Xijuan,Gao Yanli,Liu Yong.Optimization Algorithm for Compu-ting Exact Minkowski Sum of 3D Convex Polyhedra[J].InternationalJournal of Innovative Computing,Information and Control,2008,4(6):1401-1410.

[2]Jia Tianqi,Feng Xue,Guo Xijuan.Improved Marching Cubes algorithmbased on sign determination[C].The 2nd International Workshop Edu-cation Technology and Computer Science,2010(5):194-197.

[3]Guo Xijuan,Xie Lei.Optimal Accurate Minkowski Sum Approximation[C]//Advanced Intelligent Computing Technology and Applications-ICIC2008,Lecture Notes in Computer Sciences(LNCS),ICIC 2008,LNCS 5226:197-206.

[4]Varadhan G,Manocha D.Accurate Minkowski Sum Approximation ofPolyhedral Models[C]//Proc.Comput.Graph and Appl,12th PG’04,IEEE Comput.Sci,2005:392-401.

[5]孙伟,张彩明,杨兴强.Marching Cubes算法研究现状[J].计算机辅助设计与图形学学报,2007,19(7):947-952.

[6]徐晓玲,李现民,李桂清,等.体素重建中的快速移动立方体方法[J].系统仿真学报,2002,14(4):509-513.

[7]袁友伟.采用GPU加速的三维实体模型绘制[J].电子学报,2008,36(12A):144-146.

[8]陈鹏,杨超,吴玲达.硬件加速的等值面提取与绘制[J].小型微型计算机系统,2008,8(8):1538-1541.

[9]王世东,张佑生,偶春生,等.基于能量最小化的网格优化算法[J].计算机工程与设计,2007,28(2):267-269.

[10]吴恩华,柳有权.基于图形处理器(GPU)的通用计算[J].计算机辅助设计与图形学学报,2004,16(5):601-602.

[11]Wang Shidong,Zhang Yousheng,Ou Chunsheng,et al.Algorithm formesh optimization based on energy minimization[J].Engineering andDesign,2007,28(2):267-269.

精确算法 篇7

心理学研究[1]表明, 眼睛、嘴等器官在人脸感知过程中起着重要作用, 其中眼睛的作用尤为重要。绝大多数人脸识别算法[2,3]在提取人脸特征之前需要根据器官位置进行人脸的几何仿射校正等预处理;一些识别方法还要特别提取器官部位的特征以及器官间的距离等信息作为识别特征。

眼睛定位是一个困难的问题。即使是在已知人脸粗略位置的情况下, 眼睛的精确定位仍很困难。这主要是由外界干扰和器官本身复杂的变化造成的。除了由光照、尺寸、姿态、平面旋转、图像质量等因素造成的变化外, 眼睛的开闭、眼镜的反光、头发和镜框的遮挡等也给眼睛的外观造成较大的变化;特别是在眼睛闭合的情况下, 眉毛以及粗镜框会给眼睛定位带来较大干扰。目前的眼睛定位算法可以分为基于启发式规则的方法和基于样本统计学习的方法。基于启发式规则的方法一般结合图像处理方法, 根据眼睛分布、形状、颜色、边缘等先验知识来定位眼睛。而基于样本统计学习的方法一般需要收集大量样本, 对模型的参数进行训练获取, 并采用训练获得的模型对未知样本进行定位。

基于启发式规则的方法, 由于其假设模型一般无法代表较多情况下的外观, 而且参数获取通常会很繁复, 需要大量的人工干预, 还容易受主观影响, 可扩展性很差, 对于一些特殊条件无法处理, 对于新的没有考虑的情况也缺乏可扩展性, 因而应用受到很多限制。而以往的基于样本统计学习的方法, 常常采用与人脸相似的处理方式, 忽略了人脸外观模式强于眼睛外观模式这一特点。由于, 眼睛的外观比人脸变化程度更大, 模式更加不确定, 因而, 采用普通样本统计方法获得的定位算法或者由于限制太死, 无法定位到特殊情况下的器官位置, 或者由于限制太少, 使得检测结果存在大量虚警, 造成定错, 或者受虚警影响而偏离真实位置。这些不可克服的缺点表现为或者无法定位到侧面、戴眼镜等特殊模式下的特征点, 或者会错把眉毛当作眼睛发生定位错误。

现采用经典的Real Adaboost算法进行人眼检测, 这种算法不仅稳定, 而且能处理多种模式。而且其已在人脸检测中得到成功应用, 成为人脸检测的一种主流算法, 也将证明其在人眼检测中体现很好的性能。

1 Haar-like特征及积分图像

在进行目标检测时, 往往用一些特征值来表达图像, 而非基于图像的像素值。Viola[4]首先将Haar-like特征引到人脸检测中去, 并体现出了非常好的性能。现将证明Haar-like特征在眼睛定位中也有很好的效果。

1.1 Haar-Like特征

Haar-like特征是一种矩形特征 (见图1) 。其中A和B是两个矩形框的特征, C是三个矩形框的特征, D是四个矩形框的特征, 而且每个特征内所包含的矩形框是全等的。每个Haar-like特征的特征值是白色矩形框的所有像素的和减去灰色矩形框内的所有像素的和。

在实际计算中, 每个矩形特征在子窗口中进行滑动, 并在每个位置该特征又有多个尺度, 因此在一子图像中将包含大量的矩形特征。比如在我们实验程序中, 对于24×12的子图像, 都要计算12 900个haar-like 特征。再考虑到Adaboost训练时往往需要几万个, 甚至几十万个24×12的子图像, 显然这是一个庞大的计算量。为解决这个问题, 引入了积分图像的概念。

1.2 积分图像

为了解决快速计算Haar-like 特征的问题, 采用计算积分图像[5]的概念。它可以在多种尺度下, 用相同的时间来计算不同的特征。在积分图像中, 在位置 (x, y) 的值表示在它左上方所有像素的和, 即

$i{}_i(x,y)={\displaystyle \sum _{x^{\prime }\le x,y^{\prime }\le y}i}(x^{\prime },y^{\prime })$ (1)

式 (1) 中ii (x, y) 表示积分图像, i (x, y) 是原始图像。而积分图像可通过以下两公式计算出来

s (x, y) =s (x, y-1) +i (x, y) (2)

ii (x, y) =ii (x-1, y) +s (x, y) (3)

式 (3) 中s (x, y) 表示该列像素的行积分和。而且规定, s (x, -1) =0, ii (-1, y) =0。

这样可以计算出图像的积分图像, 在积分图像里计算任意尺度的矩形框内的像素和的时间都是固定的。比如图2, 在D里面的像素的和是4+1- (2+3) 。

2 Real Adaboost算法介绍

Adaboost算法通过不断地调用弱学习机来得到一系列弱分类器, 然后把这些弱分类器合并成一个强分类器。Adaboost有很多变种, 比如Discrete Adaboost, Gentle Adaboost及本文要研究的Real Adaboost等。这些算法都来自于相同的思想:就是由很多弱分类器来组合成强分类器。他们的主要不同在于强分类器的构造方式和弱分类器的构造方式。

Real Adaboost的弱分类器是把样本空间χ映射到实数空间R。

2.1 Real Adaboost强分类器的构造算法

Real Adaboost算法的强分类器训练算法流程如下。

(1) 给定训练集L={ (xi, yi) }, i=1, …, n, ;yi∈{+1, -1}是样本标号, xi∈X是样本特征。

(2) 样本的初始权重$D{}_1(i)=\frac{1}{n},i=1,\cdots ,n$;

(3) 迭代次数t=1, …, T:

1) 在该轮分布Dt上, 采用弱分类器构造选择算法, 获取此轮最佳的弱分类器, 得到

ht:X→{-∞, +∞};

2) 更新样本的权重

$D{}_{t+1}(i)=\frac{D{}_t(i)\exp (-y{}_{i}h{}_t(x{}_i))}{{\rm Z}{}_t}$ (4)

${\rm Z}{}_t={\displaystyle \sum _{i}D}{}_t(i)\exp (-y{}_{i}h{}_t(x{}_i))$ (5)

其中Zt是归一化因子;

(4) 输出的强分类器:

H (x) =sign$({\displaystyle \sum _{i=1}^{{\rm T}}h}{}_t(x)-b)$ (6)

2.2 弱分类器的构造算法

在每轮迭代中, 对于弱分类器构造算法, 我们采用区域分割弱分类器构造方法来构造弱分类器, 并从中选取使得分类误差上界最小的弱分类器作为本轮输出。

弱分类器构造算法如下。

(1) 在分布Dt上, 对于候选特征空间H中的每个特征构造其对应的弱分类器如下:

1) 将样本空间分为n个不同的区间X1, …, Xn, X1∪X2∪…Xn=X且$X{}_i{\displaystyle \underset{i\ne j}{\cap }}X{}_j=\Phi$;

2) 在分布Dt上, 计算:

$W{}_l^j=p(x{}_i\in X{}_j,y{}_i=l)={\displaystyle \sum _{x{}_i\in X{}_j,y{}_i=l}D}{}_t(i)$

其中l=±1;

3) 对于Xj中的每个x, 设定其对应弱分类器输出为:

$\forall x\in X{}_{j,}h(x)=\frac{1}{2}\ln \left(\frac{W{}_{+1}^j+\epsilon }{W{}_{-1}^j+\epsilon }\right)$ (7)

其中ε<<1/2N, 引入ε是为了避免分母为零时引起的数字溢出;

4) 计算

${\rm Z}=2{\displaystyle \sum _j\sqrt{W{}_{+1}^j}}W{}_{-1}^j$ (8)

(2) 从构造的所有弱分类器中, 选择使得Z最小的h作为此轮最终选择输出的弱分类器, 即$h{}_t=\underset{h\in {\rm H}}{{\displaystyle \text{a}\text{r}\text{g}\text{m}\text{i}\text{n}}}{\rm Z}$。

2.3 Cascade结构

Viola提出Cascade结构 (见图3) 可显著加快检测速度, 因此层次型检测器分为多层, 每一层都是上述Real Adaboost 算法训练得到的一个强分类器。经过阈值调整, 使得每一层都能让几乎全部人脸样本通过, 而拒绝很大一部分非人脸样本。这样, 靠近前面的各层, 采用少量的特征却拒绝了大部分的非人脸样本, 靠近后面的各层, 虽然采用大量的特征来排除近似人脸的非人脸候选图像的干扰, 但是由于需要处理的窗口数目很少, 对于整体运算时间的耗费很小, 就在保证检测效果的同时, 提高了人脸检测速度。

3 Nest结构的分类器

在Cascade结构中, 每个节点都是Real Adaboost分类器, 他们之间的相关性很小。对于后一层来说, 所有通过前一层的样本都是目标。实际上, 由Real Adaboost算法可知, 每层分类器的输出为

${\rm H}(X)={\displaystyle \sum _{t=1}^{{\rm T}}h}{}_t(X)$ (9)

H (X) 是各个弱分类器输出的和, 是一个实数值。

为了观察H (X) 对正样本和负样本的的不同响应, 设计了这样一个实验:首先训练一个单层的Real Adaboost分类器, 然后记录了通过该分类器的500个正样本和500个负样本的输出值, 见图4。

由图4可以明显看出, 虽然这500个正样本和500个负样本都通过了单层Real Adaboost分类器, 也就是他们的输出值都大于该层分类器的门限, 但该分类器对正样本的输出值明显大于对负样本的输出值, 而且很分明。

通过上面的分析, 知道前一层分类器的输出值很重要, 在训练下一层分类器的时候, 应该利用这一宝贵信息。把所有前一个节点的输出值作为每一个样本的特征值来训练一个新的弱分类器, 这个新分类器我们叫做Nest弱分类器[7]。

修正后的Real Adaboost算法如下。

(1) 给定训练集L={ (xi, yi) }, i=1, …, n, ;yi∈{+1, -1}是样本标号, xi∈X是样本特征。

(2) 样本的初始权重$D{}_1(i)=\frac{1}{n},i=1,\cdots ,n$;

(3) 迭代次数t=1, …, T:

1) 如果这是层次型结构的第二个或以后的节点, 而且t等于1, 则利用层次型结构前一个节点的输出作为特征值训练Nest弱分类器并更新权重。

2) 在该轮分布Dt上, 采用弱分类器构造选择算法, 获取此轮最佳的弱分类器, 得到ht:X→{-∞, +∞};

3) 更新样本的权重

$D{}_{t+1}(i)=\frac{D{}_t(i)\exp (-y{}_{i}h{}_t(x{}_i))}{{\rm Z}{}_t}$,

${\rm Z}{}_t={\displaystyle \sum _{i}D}{}_t(i)\exp (-y{}_{i}h{}_t(x{}_i))$,

其中Zt是规一化因子。

4) 输出的强分类器:

H (x) =sign$({\displaystyle \sum _{i=1}^{{\rm T}}h}{}_t(x)-b)$。

与普通层次型结构分类器不同, Nest结构分类器充分利用了前一分类器的信息。通过用训练的Nest弱分类器, 后一个分类器就继承了前一个分类器的信息。Nest层次型分类器结构见图5。

4 实验结果

为了比较Nest结构和普通层次型分类器的性能, 设计了一个实验:首先收集了10 000个左眼的正样本和200 000负样本, 然后采用相同的参数, 分别训练了一个Nest层次型结构的分类器和普通层次型结构的分类器。下面是部分样本。

正样本集示例:

负样本集示例:

Nest层次型结构的分类器和普通层次型结构的分类器有着相同的参数:每层的最小检测率为0.999 9, 最大误检率为0.4, 总体目标误检率为0.001。

4.1 ROC曲线

ROC曲线即检测率～虚警率曲线。检测率和虚警率分别定义如下

检测率$=\frac{\text{正}\text{确}\text{检}\text{测}\text{到}\text{的}\text{眼}\text{睛}\text{数}\text{目}}{\text{测}\text{试}\text{集}\text{中}\text{包}\text{含}\text{的}\text{眼}\text{睛}\text{数}\text{目}}\times 100\%$ (9)

误检率$=\frac{\text{虚}\text{警}\text{数}\text{目}}{\text{总}\text{搜}\text{索}\text{窗}\text{口}\text{中}\text{非}\text{眼}\text{睛}\text{的}\text{数}\text{目}}\times 100\%$ (10)

一般而言, 可以采用虚警率作为横坐标, 检测率作为纵坐标构建ROC曲线。不断改变检测器阈值的方式来获取不同阈值下的检测率和对应虚警率, 从而获取ROC曲线。对于在相同条件下训练的Nest结构和普通层次型结构分类器的ROC曲线对比如图6。

由ROC曲线对比可知, 在相同的误检率下, Nest层次型分类器有更高的检测率;而在相同的检测率下, Nest层次型分类器有更低的误检率。也就是说Nest层次型分类器有的检测性能要好于普通层次型的分类器。

4.2 分类器弱特征器数目比较

除了比较分类器的检测性能外, 又比较了分类器选择的弱特征数, 也就是弱分类器的个数。分类器选择的弱特征个数直接关系到分类的训练速度和检测速度, 其中检测速度是衡量分类器的一个很重要的指标。

所训练的Nest层次型分类器和普通分类器都含有12层Real Adaboost分类器, 但每层分类器所用的弱特征数差别很大, 如图7所示。

从图7可以看到, Nest层次型分类器每层所用的弱特征数明显少于普通层次型结构每层所用的弱特征数, 而Nest层次型分类器所用的总的弱特征数更远少于普通层次型结构所用的弱特征数。也就是说Nest层次型结构的分类器有更快的检测速度。实验表明, 对于24×12的样本来说, 基本层次分类器每个需要大概138.2 μs, 而Nest层次型结构的分类器平均仅需要77.8 μs。

4.3 检测结果示例

按照上述的方法, 同样又训练了一个Nest层次型的右眼分类器。为了衡量Nest层次型分类器的性能, 做了以下定义。

对于一张人脸, 假定手工标定的左眼、右眼位置分别为 (P${}_{\text{l}\text{e}}^{\text{l}\text{a}\text{b}\text{e}\text{l}}$, P${}_{\text{r}\text{e}}^{\text{l}\text{a}\text{b}\text{e}\text{l}}$) , 自动定位的左眼、右眼位置分别为 ( (P${}_{\text{l}\text{e}}^{\text{D}\text{e}\text{t}}$, P${}_{\text{r}\text{e}}^{\text{D}\text{e}\text{t}}$) , 自动定位结果与手工标定结果之间的欧氏距离对于左眼、右眼分别为 (dle, dre) 。手工标定左右眼之间的欧式距离为deye。

则眼睛定位误差定义为:

$\epsilon =\frac{\max (d{}_{le},d{}_{re})}{d{}_{eye}}$ (11)

当眼睛定位误差ε<0.15时, 认为双眼的定位是准确的。

按照上述定义, 我们在Feret Person数据库上测试了我们的算法, 结果是我们的左眼检测率为98.57%, 右眼检测率为99.31%。下面是部分检测结果 (图8) 。

5 结论

在他人工作的基础上, 实现了一种准确鲁棒的眼睛特征点定位算法。该算法不仅能够检测到正面正常光照条件下的眼睛中心点, 而且对于戴眼镜、闭眼等情况也有非常好的性能。

本文贡献主要在以下几个点:

(1) 将Haar-like特征、积分图像、Real Adaboost分类器等方法应用于眼睛定位。之前Viola等人已经证明这些方法能很好地解决人脸检测的问题, 本文证明了这些算法也能很好地应用于眼睛定位;

(2) 仔细分析了普通层次型分类器的不足之处, 并通过实验给出了直观的说明;

(3) 针对普通层次型分类器的不足, 介绍了Nest结构的层次型分类器。在相同的条件下训练了一个普通层次型分类器和Nest结构分类器, 通过比较了它们的ROC曲线和使用的弱特征的个数, 证明了Nest层次型分类器能更快和更准确地检测眼睛。

参考文献

[1]Keating C F, Keating E G.Monkeys and mug shots:cues used by Rhesus monkeys (Macaca mulata) to recognize a human face.Jour-nal Comparative Psychology, 1993;107 (2) :131—139

[2] Belhumeur P N, Hespanha J P. Eigenface vs Fisherfaces: recognition using class specific linear projection. IEEE Trans PAMI, 1997; 19:711—720

[3] Brunelli R, Poggio T. Face recognition: features versus templates. IEEE Trans on PAMI, 1993; 15 (10) : 1042—1052

[4]Viola P, Jones M.Robust real time object detection.IEEE ICCV Workshop on Statistical and Computational Theories of Vision, Van-couver, Canada, July13, 2001

[5] Jones M, Viola P. Fast multi-view face detection. Technical Report 2003-96, Mitsubishi Electric Research Lab, July 2003

[6]Schapire E, Singer Y.Improved boosting algorithms using confi-dence-rated predictions.Proceedings of the Eleventh Annual Confer-ence on Computational Learning Theory, 1998:80—91