1、1,粒子群优化算法求优 及LBG算法的运行,深圳大学信息工程学院 黄彩玲,2,一. 粒子群优化算法求最优解,初始化一群随机粒子（随机解）每次迭代中，粒子通过跟踪两个极值更新自己：粒子本身找到的历史最后解（个体极值点pbest）整个种群目前找到的最好解（全局极值点gbest）需要计算粒子的适应值，以判断粒子位置距最优点的距离。 每次迭代中，根据适应度值更新pbest和gbest。 迭代中止条件：设置最大迭代次数或全局最优位置满足预定最小适应阈值。,3,粒子群优化算法求最优解,D维空间中，有N个粒子；粒子i位置：xi=(xi1,xi2,xid)，将xi代入适应函数f(xi)求适应值；粒子i速度：v
2、i=(vi1,vi2,vid)粒子i个体极值点位置：pbesti=(pi1,pi2,pid)种群的全局极值点位置：gbest=(g1,g2,gd)粒子i的第n维速度和位置更新公式：vin=w*vin+c1*r1*(pbestin-xin)+c2*r2*(gbestn-xin)xin=xin+vinc1,c2学习因子，经验值取c1=c2=2，调节学习最大步长r1,r2两个随机数，取值范围（0，1），以增加搜索随机性w 惯性因子，非负数，调节对解空间的搜索范围,4,基本粒子群优化算法流程图,5,程序分析,主要数据结构：种群大小(PopSize) 空间维数(NDim)矢量的边界(Bound) 最大迭
3、代次数(MaxIter) C1、C2、W、R1、R2 各粒子当前适应度值(fvalue)各粒子位置(population) 各粒子速度(velocity) 各粒子的最佳位置(pbest) 全局最佳粒子位置(gbest) 全局最佳粒子序号(index) 更新前各粒子适应度值(fpbest)得到相近适应度值的迭代次数(samecounter)临时适应度值 (oldbestval),6,初始化各主要数据（设三维的Sphere函数求最优）,flag=0; %停止程序标志 oldbestval=0; %记录旧的适应度值 samecounter=0; %记录得到相同适应度值的迭代次数 iteration
4、= 0; %迭代次数 MaxIter =100; %最大迭代次数 PopSize=20; %种群大小 c1 = .5; %学习因子 c2 = .5; %学习因子 w=0.8; %惯性因子 Bound=-100 100;-100 100;-100 100; %粒子的坐标范围 NDim = length(Bound); %空间维数 for i=1:PopSize %定义粒子上下边界lowerbound(:,i)=Bound(:,1);upperbound(:,i)=Bound(:,2); end,7,初始化各主要数据, for i=1:Ndim %初始化各粒子初始位置，在有效范围内随机选数 pop
5、ulation(i,:)=rand(1, PopSize)*(Bound(i,2)-Bound(i,1) + Bound(i,1); end for i=1:Ndim %初始化各粒子最大速度，使粒子不能越出边界vmax(i,:)=(Bound(i,2)-Bound(i,1)/2; end velocity = vmax.*rand(NDim, PopSize); for i = 1:PopSize %计算各粒子的适应度值fvalue(i) = population(1,i)2+population(2,i)2+population(3,i)2; endpbest = population; %
6、记录各粒子的个体极值点位置 fpbest = fvalue; %记录最佳适应度值 fbestval,index = min(fvalue); % 找出全局极值和相应的序号,8,主程序,while(flag = 0) % 阻止逸出,9,主程序,for i = 1:PopSize % 更新各粒子适应度值fvalue(i) = population(1,i)2+population(2,i)2+population(3,i)2;endchangeColumns = fvalue fpbest; % 更新后的适应度值优于更新前的，记录序号pbest(:, find(changeColumns) = p
7、opulation(:, find(changeColumns); % 更新个体极值点位置fpbest = fpbest.*( changeColumns) + fvalue.*changeColumns; %更新个体极值fbestval, index = min(fvalue); %更新全局极值和相应的序号if floor(fbestval*1e30)=oldbestval %比较更新前和更新后的放大的适应度值;samecounter=samecounter+1; %相等时记录加一;elseoldbestval=floor(fbestval*1e30); %不相等时更新放大的适应度值，并记录
8、清零;samecounter=0;end,10,主程序,if samecounter = 20 %多次迭代的适应度值相近时程序停止flag=1;endBest(iteration) =fbestval; % 输出及描出找到的全局极值plot(Best,ro);xlabel(generation); ylabel(f(x);text(0.5,0.95,Best = , num2str(Best(iteration),Units,normalized); drawnow; end,11,运行结果(添加三维图等，使之直观点， 达优范围要查),12,二. LBG算法矢量量化码书设计,LBG算法以初始码
9、字开始，不断迭代直至收敛，每个迭代过程包括对训练矢量分类和更新码字。缺点：计算量大，生成的码书无序，码书质量受初始码书影响。基于两条优化准则：最近邻域准则。对于给定码书，训练矢量集的最优分类可通过把每个矢量映射为离它最近的码字而得到。质心条件。对于给定的训练矢量分类，其对应的最优码书中各码字是通过求各胞腔的中心矢量而获得。,13,LBG矢量量化码书程序流程图,14,程序分析,主要函数：main_lbg(): 主程序cal_ed_ad()：由矢量和码书得到矢量和码字距离，矢量索引号，空胞腔数目 LBG()：由训练矢量、旧码书，旧ed、矢量索引号等进 行LBG计算，得到新码书，新的矢量索引号、新的
10、矢量失真、空胞腔数目,15,main_lbg()主要程序分析,V=50; % 迭代次数 epsilon=0.001; % 出错率极限 Line,Col=size(img_lena); % 读出图像的大小 n=4; % 设定向量维数为n*n M=Line*Col/n/n; % 训练矢量大小 tv=zeros(Line*Col/n/n,n*n); for i=0:Line/n-1 % 分割图像，得到矢量，tv(4096*16)for j=0:Col/n-1 tv(i*Col/n+1+j,:)=img_lena(i*n+1,j*n+1:j*n+n),img_lena(i*n+2,j*n+1:j*n+
11、n),img_lena(i*n+3,j*n+1:j*n+n),img_lena(i*n+4,j*n+1:j*n+n); (令f=i*n使程序速度提高)end end,16,main_lbg()主要程序分析,cbook_new=tv(round(rand(1,256)*M),:); cbook_old=cbook_new; % cbook_old为更新前码书，cbook_new为更新后码书，初始化 N,L=size(cbook_old); % N=256,L=16，码书大小 img_new=zeros(M,L); % 训练矢量更新, 4096 x 16 ety_cv=zeros(1,V+1);
12、% No_of_empty_cv 每次迭代后空胞腔数目，共V+151次迭代ed_old,ad,ind_tv_old,no_ecv=cal_ed_ad(tv,cbook_old,0); % 由矢量和码书得到矢量和码字距离，矢%量索引号，空胞腔数目 D(1)=ad; %第一次迭代后的平均失真 ety_cv(1)=no_ecv; %第一次迭代后的空胞腔数目 img_ini=zeros(M,L); for i = 1:M,img_ini(i,:)=cbook_old(ind_tv_old(i),:); % 重新划分后的矢量 End psnr=10*log10(GL*GL/mean(mean(img_i
13、ni-tv).2) %第一次迭代后的信噪比 PSNR(1)=psnr;,17,main_lbg()主要程序分析,for v = 1:V %迭代50次if stop_flag=0 cbook_new,ind_tv_new,ed_new,img_new,ad,psnr,no_ecv.= LBG(tv,cbook_old,ed_old,ind_tv_old,L,M,N,GL); %由训练矢量、旧码书，旧ed、矢量索引号等进 行LBG计算，得到新码书，新的矢量索引号、新的矢量失真、空胞腔数目no_ecv=no_ecv %空胞腔数目ety_cv(v+1)=no_ecv; %第v1次迭代后的空胞腔数目ad
14、r=abs(D(v+1)-D(v)/D(v) % 平均失真率if adr=epsilon % 如果平均失真率达到要求就停止运算stop_flag=1; elsecbook_old=cbook_new; % 否则更新码书ind_tv_old=ind_tv_new; % 更新矢量索引号ed_old=ed_new; % 更新所有训练矢量与所有码字的距离endend end,18,cal_ed_ad()主要程序分析,for i = 1:M,for j = 1:N,ed(i,j)=norm(tv(i,:)-cb(j,:); % 计算每个训练矢量与每个码字的距离end end if cho=0 y_ed,
15、ind_tv=min(ed,2); % 找每个矢量对应最近码字及相应的索引号for i = 1:M,ad=ad+(norm(tv(i,:)-cb(ind_tv(i),:)2/M; %计算平均失真率endno_ecv=0; %记录空胞腔数目identified_cb=zeros(1,N); % 值为1 是非空, 0 是空胞腔 for i = 1:M,identified_cb(ind_tv(i)=1; % 非空矢量作标记1endno_ecv=N-sum(identified_cb); % 计算空胞腔数目 end,19,LBG()主要程序分析,for i = 1:M %产生想得到新的码书用，码书为
16、胞腔内矢量之和mf_lbg(i,ind_tv_old(i)=1; end cb_new=mf_lbg*tv; %码书为胞腔内矢量之和for j = 1:N %更新码书no_tv=sum(mf_lbg(:,j); if no_tv=0cb_new(j,:)=cb_old(j,:); % 为空胞腔时这个码字不变elsecb_new(j,:)=cb_new(j,:)/no_tv; %非空时取中间值产生新码字end Endcb_new_r=round(cb_new); % 四舍五入取整，新码书ed_new,ad,ind_tv_new,no_ecv=cal_ed_ad(tv,cb_new_r,0); %更新ed和ad,20,LBG()主要程序分析,for i = 1:M %由新的码书产生新的训练矢量img_new(i,:)=cb_new_r(ind_tv_new(i),:); endpsnr=10*log10(GL*GL/mean(mean(img_new-tv).2);,21,运行结果,Lena原图 LBG算法译码图像,22,运行结果,PSNR 28.7641 耗时：59.3350s 平均失真：1.3829e+003,23,谢谢！,
