粒子群优化算法求优 及LBG算法的运行 深圳大学信息工程学院 黄彩玲

一. 粒子群优化算法求最优解 初始化一群随机粒子（随机解） 每次迭代中，粒子通过跟踪两个极值更新自己： －粒子本身找到的历史最后解（个体极值点pbest） －整个种群目前找到的最好解（全局极值点gbest） 需要计算粒子的适应值，以判断粒子位置距最优点的距离。 每次迭代中，根据适应度值更新pbest和gbest。 迭代中止条件：设置最大迭代次数或全局最优位置满足预定最小适应阈值。

粒子群优化算法求最优解 D维空间中，有N个粒子； 粒子i位置：xi=(xi1,xi2,…xid)，将xi代入适应函数f(xi)求适应值； 粒子i速度：vi=(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+vin c1,c2—学习因子，经验值取c1=c2=2，调节学习最大步长 r1,r2—两个随机数，取值范围（0，1），以增加搜索随机性 w —惯性因子，非负数，调节对解空间的搜索范围

根据适应度更新pbest、gbest，更新粒子位置速度 基本粒子群优化算法流程图 开始 初始化粒子群 计算每个粒子的适应度 根据适应度更新pbest、gbest，更新粒子位置速度 结束 no yes 达到最大迭代次数或 全局最优位置满足最小界限？

程序分析 主要数据结构： 种群大小(PopSize) 空间维数(NDim) 矢量的边界(Bound) 最大迭代次数(MaxIter) C1、C2、W、R1、R2 各粒子当前适应度值(fvalue) 各粒子位置(population) 各粒子速度(velocity) 各粒子的最佳位置(pbest) 全局最佳粒子位置(gbest) 全局最佳粒子序号(index) 更新前各粒子适应度值(fpbest) 得到相近适应度值的迭代次数(samecounter) 临时适应度值 (oldbestval)

初始化各主要数据（设三维的Sphere函数求最优） flag=0; %停止程序标志 oldbestval=0; %记录旧的适应度值 samecounter=0; %记录得到相同适应度值的迭代次数 iteration = 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

初始化各主要数据 … for i=1:Ndim %初始化各粒子初始位置，在有效范围内随机选数 population(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; velocity = vmax.*rand(NDim, PopSize); for i = 1:PopSize %计算各粒子的适应度值 fvalue(i) = population(1,i)^2+population(2,i)^2+population(3,i)^2; pbest = population; %记录各粒子的个体极值点位置 fpbest = fvalue; %记录最佳适应度值 [fbestval,index] = min(fvalue); % 找出全局极值和相应的序号

主程序 while(flag == 0) & (iteration < MaxIter) %寻找最优主程序 iteration = iteration +1; %迭代次数递增 for i=1:PopSize %更新全局极值点位置 gbest(:,i) = population(:,index); end R1 = rand(NDim, PopSize); %重新设置随机性 R2 = rand(NDim, PopSize); R3 = rand(NDim, PopSize); velocity = w*velocity + c1*R1.*(pbest-population) + c2*R2.*(gbest-population); % 更新各粒子速度 population = population + velocity; % 更新各粒子位置 OutFlag = population<=lowerbound | population>=upperbound; % 逸出标志 population = population - OutFlag.*velocity; % 阻止逸出

主程序 for i = 1:PopSize % 更新各粒子适应度值 fvalue(i) = population(1,i)^2+population(2,i)^2+population(3,i)^2; end changeColumns = fvalue <fpbest; % 更新后的适应度值优于更新前的，记录序号 pbest(:, find(changeColumns)) = population(:, find(changeColumns)); % 更新个体极值点位置 fpbest = fpbest.*( ~changeColumns) + fvalue.*changeColumns; %更新个体极值 [fbestval, index] = min(fvalue); %更新全局极值和相应的序号 if floor(fbestval*1e30)==oldbestval %比较更新前和更新后的放大的适应度值; samecounter=samecounter+1; %相等时记录加一; else oldbestval=floor(fbestval*1e30); %不相等时更新放大的适应度值，并记录清零; samecounter=0;

主程序 if samecounter >= 20 %多次迭代的适应度值相近时程序停止 flag=1; end Best(iteration) =fbestval; % 输出及描出找到的全局极值 plot(Best,'ro');xlabel('generation'); ylabel('f(x)'); text(0.5,0.95,['Best = ', num2str(Best(iteration))],'Units','normalized'); drawnow;

运行结果(添加三维图等，使之直观点， 达优范围要查) 函数名 表达式 维数 初值范围 理想达优值 运行次数 最大迭代次数 达优次数 平均值 Sphere ，-100<x<100 3 [-100,100] 50 46 0.053 Rosenbrock 5 [-15,30] 200 26 21.9756 六峰驼背函数 ，0<x<1 1 (0,1) 0.1989 41 18.8867

二. LBG算法矢量量化码书设计 LBG算法以初始码字开始，不断迭代直至收敛，每个迭代过程包括对训练矢量分类和更新码字。 缺点：计算量大，生成的码书无序，码书质量受初始码书影响。 基于两条优化准则： －最近邻域准则。对于给定码书，训练矢量集的最优分类可通过把每个矢量映射为离它最近的码字而得到。 －质心条件。对于给定的训练矢量分类，其对应的最优码书中各码字是通过求各胞腔的中心矢量而获得。

LBG矢量量化码书程序流程图 开始 初始化矢量和码书 找最接近码字，矢量重新分类 算失真和信噪比 更新码书 no 平均失真率达到要求？ 结束 输出图像 no yes

程序分析 主要函数： main_lbg(): 主程序 cal_ed_ad()：由矢量和码书得到矢量和码字距离，矢 量索引号，空胞腔数目 LBG()：由训练矢量、旧码书，旧ed、矢量索引号等进 行LBG计算，得到新码书，新的矢量索引号、 新的矢量失真、空胞腔数目

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+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

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); % No_of_empty_cv 每次迭代后空胞腔数目，共V+1＝51次迭代 [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_ini-tv).^2))) %第一次迭代后的信噪比 PSNR(1)=psnr;

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; %第v＋1次迭代后的空胞腔数目 adr=abs(D(v+1)-D(v))/D(v) % 平均失真率 if adr<=epsilon % 如果平均失真率达到要求就停止运算 stop_flag=1; else cbook_old=cbook_new; % 否则更新码书 ind_tv_old=ind_tv_new; % 更新矢量索引号 ed_old=ed_new; % 更新所有训练矢量与所有码字的距离 end

cal_ed_ad()主要程序分析 for i = 1:M, for j = 1:N, ed(i,j)=norm(tv(i,:)-cb(j,:)); % 计算每个训练矢量与每个码字的距离 end … if cho==0 [y_ed,ind_tv]=min(ed,[],2); % 找每个矢量对应最近码字及相应的索引号 ad=ad+(norm(tv(i,:)-cb(ind_tv(i),:)))^2/M; %计算平均失真率 no_ecv=0; %记录空胞腔数目 identified_cb=zeros(1,N); % 值为‘1’ 是非空, ‘0’ 是空胞腔 identified_cb(ind_tv(i))=1; % 非空矢量作标记1 no_ecv=N-sum(identified_cb); % 计算空胞腔数目

LBG()主要程序分析 for i = 1:M %产生想得到新的码书用，码书为胞腔内矢量之和 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==0 cb_new(j,:)=cb_old(j,:); % 为空胞腔时这个码字不变 else cb_new(j,:)=cb_new(j,:)/no_tv; %非空时取中间值产生新码字 End cb_new_r=round(cb_new); % 四舍五入取整，新码书 [ed_new,ad,ind_tv_new,no_ecv]=cal_ed_ad(tv,cb_new_r,0); %更新ed和ad

LBG()主要程序分析 for i = 1:M %由新的码书产生新的训练矢量 img_new(i,:)=cb_new_r(ind_tv_new(i),:); end psnr=10*log10(GL*GL/mean(mean((img_new-tv).^2)));

运行结果 Lena原图 LBG算法译码图像

运行结果 PSNR＝ 28.7641 耗时：59.3350s 平均失真：1.3829e+003

谢谢！