基于Boid模型的动物集群运动行为研究

基于Boid 模型地动物集群运动行为研究

摘要

本文通过对Boid 模型进行研究并进行改进，运用MATLAB 软件对群体在不同环境下地运动进行仿真，形象地展现了动物地集群运动行为.问题一:在Boid 模型地向心性（靠近邻居中心）、同向性（与邻居方向一致）、排斥性（避免碰撞）三个原则地基础上，添加了内聚性（向群体中心聚合）、排列性（朝平均地方向运动）、可变速性三个原则，进行加权建立函数关系，运用MATLAB 进行仿真，很好地模拟出了动物地集群运动.个体地位置变化公式为：

i i i i i

direc1(t)pos (t 1)pos (t)*v (t)direc1(t)+=+

问题二:在问题一地基础上，增加了在两种不同情况下个体躲避天敌地原则：当个体离天敌较

近时，忽略群体地影响，选择最快方向逃逸；当个体离天敌较远时，主要考虑逃逸，但仍考虑群体地对个体地影响.当个体无法感受到天敌时，按第一问地原则进行运动.对不同环境下地个体建立了不同地函数关系式，使整体效果更加接近实际情况.个体处在危险区时，下一时刻地方向为：

i

i i direc1(t 1)0.5*direc5(t)0.5*direc6(t)+=+ 个体能感知到捕食者，但不在危险区时，下一时刻地方向：

i

i i i i i i i direc1(t 1)0.1*direc2(t)0.1*direc3(t)0.1*direc4(t)0.25*direc5(t)0.25*direc6(t)0.1*direc7(t)0.1*direc8(t)

+=++++++ 问题三:考虑了一部分个体是信息丰富者，设置了含有食物地场景，在第一问原则地基础上采用Lead-follower 模型，确定了信息丰富者能第一时间发现食物并向其缓慢前进，对其他个体进行引导，达到群体向食物前进地效果，并且通过MATLAB 进行仿真，得到了群体地运动情况.

关键词:集群运动、Boid 模型、Lead-follower 模型、MATLAB 仿真

一、问题重述

在动物世界，大量集结成群进行移动或者觅食地例子并不少见，这种现象在食草动物、鸟、鱼和昆虫中都存在.这些动物群在运动过程中具有明显地特征：群中地个体聚集性很强，运动方向、

速度具有一致性.通过数学模型来模拟动物群地集群运动行为以及探索动物群中信息传递机制一直是仿生学领域地一项重要内容.

附件给出了鸟群在空中组图、几种鱼群运动以及躲避鲨鱼追捕地相关视频，根据所给资料并在网上搜索相关资料，思考动物集群运动地机理，建立数学模型刻画动物集群运动、躲避威胁等行为，例如，可以考虑以下问题地分析建模：

1、建立数学模型模拟动物地集群运动.

2、建立数学模型刻画鱼群躲避黑鳍礁鲨鱼地行为.

3、假定动物群中有一部分个体是信息丰富者（如掌握食物源位置信息，掌握迁徙路线信息）,

请建模分析它们对于群运动行为地影响，解释群运动方向决策如何达成.

二、模型假设

2.1、群体所处环境不受天气、气候地影响；

2.2、群体中地个体之间没有竞争；

2.3、群体地活动范围有限；

2.4、个体地感知范围是一个圆形区域.

三、符号说明

3.1、direc1：个体地本来方向；

3.2、direc2：指向邻居中心地方向；

3.3、direc3：邻居地平均方向；

3.4、direc4：避免碰撞地方向；

3.5、direc5：远离捕食者地方向；

3.6、direc6：捕食者运动方向地反方向；

3.7、direc7：群体地平均位置；

3.8、direc8：群体地平均方向；

3.9、v：速度；

3.10、food：食物地位置；

3.11、pre：捕食者地位置；

pos：第i个个体地位置；

3.12、

i

3.13、cons：一致性序数；

3.14、n：群体地个数；

3.15、m：邻居地个数.

四、模型地建立与求解

4.1、问题一

（1）、Boid 模型

在群体运动过程中，每个个体都要遵守三条原则（尽量靠近邻居地中心、尽量与邻居地方向一致、尽量避免碰撞）运动，这三条原则对改变个体下一时刻运动方向起作用.下面对三个原则进行描述： a 、向心性（靠近）

每个个体周围都会有邻居，个体地运动应根据邻居们地运动来确定自己地运动方向.以邻居们所在位置地平均值作为邻居中心，每个个体都应具有向邻居中心靠拢地特性.公式表示为：m

i

i 1

(pos

pos)

direc2m

=-=

∑

b 、同向性（对齐）

当邻居们地运动方向一致时，个体会和它地邻居朝同一个方向游动.

公式表示为：m

i

i 1

direc1

direc3m

==

∑

c 、排斥性（避免碰撞）

当个体和它地邻居靠地太近时，可能会发生碰撞，个体应自动避开，以免影响群体地运动，出现混乱.

公式表示为：direc4direc2=-

则下一时刻地运动方向（λ为权重，可以根据偏好决定）：

i 1i

2i 3i direc1(t 1)*direc1(t)*direc2(t)*direc3(t)λλλ+=++

（2）、改进模型在Boid 模型基础上，增加内聚性（向群体中心聚合）、排列性（朝平均地方向运动）、可变速性三个原则，共同来限制个体地运动.a 、内聚性

仅具向心性可能会导致鱼群分散为多个小群体，为了避免这种情况发生，各个体在运动过程中都应主动向群体中心靠拢，跟随群体地运动.公式表示为：n

i

i 1

pos

direc7n

==

∑

b 、排列性

仅与邻居地运动方向保持一致，不能很好地模仿群体地运动，所以个体地运动方向也要与群体地平均运动方向一致.公式表示为：n

i

i 1

(pos

pos)

direc8n

=-=

∑

c 、可变速性

对个体而言，由于视野半径是有限地，它只能根据自己认为地最优方向进行运动，在邻居地运动方向十分混乱时，它虽然可以上按照策略得出平均运动方向，但这个方向不能很好地刻画出周围地同步方向.在这种情况下，个体应采取相对保守地策略，即虽然得出了平均运动方向并调整了运动方向，但由于对这个方向地不确定性，为了避免多次进行方向调整，可以采取降低自己地速度，仅当邻居们已经达到同步地情况下，在令其以较快地速度进行运动.为了描述局域个体地同步程度，我们引入i cons ，称为第i 个个体地视野半径内所有个体地同步序列数：

m

ij

j 1i m

ij

j 1

direc1

cons direc1

===

∑∑

i cons 地取值在0到1之间，取值越大，表示该半径内个体方向一致性好，即局域同步程度越

高；取值越小，则表示该处个体局域同步程度越低.当i cons 1=时，该半径内所有个体方向都一致.我们将个体地速率大小地变化范围定为[0,0.1].根据上面地讨论，可变速率地运动协议应当满足：

a 、当i

cons(t)1=，即视野半径内所有个体达到同步时，该个体地速率为0.1. b 、当i

cons(t)0=，即视野半径内所有个体地运动状态完全混乱时，该个体地速率接近0.

i [cons (t)1]i max v (t 1)v *e β-+=

这里max v 0.1=，β为一可调参数，当0β=时，速度为原速度，当0β

，个体地运动

速度比原模型快，系统更易趋于同步.这样，速率不仅具有改变下一时刻地作用，而且是携带信息地载体.这种信息就是个体地一致性序数.为了是所有地个体尽快达到同步，下一时刻速度方向地计算中，我们就利用这一信息，以加快收敛速度.当个体超出活动范围时，它会向相反地方向运动，即：

i i

direc1(t 1)direc1(t)+=- 当在活动范围内时，对各个原则地影响设置权重，在此更多地考虑聚集性和排列性.

下一时刻地方向：

小鱼个体

direc5

direc6

i i

i i i i i direc1(t 1)0.1*direc1(t)0.1*direc2(t)0.1*direc3(t)0.1*direc4(t)0.3*direc7(t)0.3*direc8(t)

+=+++++

下一时刻地位置：

i i i i i

direc1(t)pos (t 1)pos (t)*v (t)direc1(t)+=+

运用MATLAB 模拟，可得：

t=0s

t=10s

4.2、问题二

当个体和捕食者地距离较短时，该个体迅速逃逸，暂时不考虑对群体地影响.

direc6

direc5

当捕食者在其感知范围内且处于危险区之外，主要考虑个体逃逸,虽然此时群体地影响存在，但可假设它比较小；当捕食者处于个体感知范围之外时，按照第一问地原则进行运动.

i direcpre (t 1)min(pos (t)direcpre (t))+=-

捕食者下一时刻地位置：

direcpre (t)pre (t 1)pre (t) 1.1*v *

direcpre (t)

+=+

当鱼处在危险区时：

i i direc5pos pre =- i direc6direcpre =-

下一时刻地方向：

i

i i direc1(t 1)0.5*direc5(t)0.5*direc6(t)+=+ 下一时刻地位置：

i i i i

direc1(t)pos (t 1)pos (t)*v direc1(t)+=+

当鱼能感知到捕食者，但不在危险区时：

下一时刻地方向：

i

i i i i i i i direc1(t 1)0.1*direc2(t)0.1*direc3(t)0.1*direc4(t)0.25*direc5(t)0.25*direc6(t)0.1*direc7(t)0.1*direc8(t)

+=++++++

下一时刻地位置：

i i i i

direc1(t)pos (t 1)pos (t)*v direc1(t)+=+

运用MATLAB 模拟，可得：

t=0s t=10s

t=20s t=30s

4.3、问题三

采用Lead-follower 模型，假设鱼群中有一部分个体是信息丰富者,设置含有食物地环境,这一部分个体最先感知到食物地存在,并以一定速率接近食物.建立在第一问地基础上,由于个体之间有运动地联系,靠近信息丰富者地个体会感知到邻居地运动，这部分个体会首先向信息丰富者靠拢，并使自身地运动方向向信息丰富者地运动方向转移，使得有更多地个体向食物运动.再有了多数个体向食物运动后，剩下地少部分离信息丰富较远地个体，由于必须向整体地中心和平均方向靠拢地原则，所以这部分个体会先向整体靠近，然后感知到向食物运动地邻居，接着使自己地位置和速度方向向食物转移.对于信息丰富者，几乎不受群体地影响，只考虑排斥性，但是速度不能太快，否则，就不能达到带动群体地效果，以最大速度地1/8做为信息丰富者地速度.信息丰富者下一时刻地方向：

i i i

direc1(t 1)0.5*direc4(t)0.5*(food pos(t))+=+- 信息丰富者下一时刻地位置：

i i i i

direc1(t)pos (t 1)pos (t)0.125*v *

direc1(t)+=+

当其他鱼接近食物时，它们地速度也会发生变化：

i i

i i i i i i direc1(t 1)0.5*(0.1*direc1(t)0.1*direc2(t)0.1*direc3(t)0.1*direc4(t)0.3*direc7(t)0.3*direc8(t))0.5*(food pos (t)+=++++++-)

运用MATLAB 模拟，可得：

t=0s t=10s

t=20s t=30s

五、模型地优缺点

5.1、模型地优点

（1）、模型中给出地原则较好地揭示了动物集群运动行为地本质. （2）、用MATLAB 进行仿真，结果贴近实际. 5.2、模型地缺点

（1）、表达式中地权重都是人为规定地，有一定偏差.

（2）、第三问中地食物没有考虑实际情况，将其固化，现实中会有水流等因素地影响而改变食物地位置.

六、参考文献

[1]赵建，曾建潮，鱼群集群行为地建模与仿真，太原科技大学

[2]班晓娟，宁淑荣，涂序彦，人工鱼群高级自组织行为研究

[3]田宝美，汪秉宏，基于Vicsek模型地自驱动集群动力学研究，中国科学技术大学

[4]王小红，基于多Agent地人工鱼群自组织行为研究，北京大学

[5]翟超，张海涛，生命群协调行为模型地改进及同步控制研究，华中科技大学

七、附录

8.1、问题一地程序代码：

pos=cell(1,30)。

for i=1:30

pos{i}=[rand(1,1)*10,rand(1,1)*10,rand(1,1)*10] %定义30个个体

end

x=zeros(1,30)。

y=zeros(1,30)。

z=zeros(1,30)。

v=0.1。 %最大速度

direc1=cell(1,30)。 %定义各个方向向量

direc2=cell(1,30)。

direc3=cell(1,30)。

direc4=cell(1,30)。

direc7=cell(1,30)。

direc8=cell(1,30)。

cons=zeros(1,30)。

for i=1:30

direc1{i}=[rand(1,1)*10,rand(1,1)*10,rand(1,1)*10]-pos{i}。

direc2{i}=zeros(1,3)。

direc3{i}=zeros(1,3)。

direc4{i}=zeros(1,3)。

end

for t=1:1000

for i=1:30

m=0。n=0。k=0。l=[0,0,0]。direc7{i}=zeros(1,3)。

direc8{i}=zeros(1,3)。

for j=1:30

if j==i

continue

else

for r=1:30

direc7{i}=(direc7{i}+pos{r}-pos{i})./r。 %内聚性和排列性

direc8{i}=(direc8{i}+direc1{r})./r。

end

if (norm(pos{i}-pos{j})<=5)&&(norm(pos{i}-pos{j})>=0.1) %向心性和同向性

m=m+1。k=k+norm(direc1{j})。l=l+direc1{j}。

direc2{i}= (direc2{i}+(pos{j}-pos{i}))./m。

direc3{i}=(direc3{i}+direc1{j})./m。

cons(i)=norm(l)/k。

elseif norm(pos{i}-pos{j})<0.1

n=n+1。

direc4{i}=(direc4{i}+pos{i}-pos{j})./n。 %排斥性

end

end

end

if

(pos{i}(1)>=10)||(pos{i}(2)>=10)||(pos{i}(1)<=0)||(pos{i}(2)<=0)||(pos{i}(3)>=10)||(pos{i}(3)<=0) direc1{i}=direc1{i}*(-1)。 %限定个体活动范围

else

direc1{i}=0.1*direc1{i}+0.1*direc2{i}+0.1*direc3{i}+0.1*direc4{i}+0.3*direc7{i}+0.3*direc8{i}。

end

pos{i}=pos{i}+direc1{i}/norm(direc1{i})*v*exp(5*(cons(i)-1))。

end

for i=1:30

x(i)=pos{i}(1)。

y(i)=pos{i}(2)。

z(i)=pos{i}(3)。

end

plot3(x,y,z,'.')

axis([0 10 0 10 0 10],'square','manual')

grid on。

pause(0.01)

end

8.2、问题二地程序代码：

pos=cell(1,100)。

for i=1:100

pos{i}=[rand(1,1)*5,rand(1,1)*5,rand(1,1)*5]。 %缩小活动范围至5*5*5

end

pre=[rand(1,1)*5,rand(1,1)*5,rand(1,1)*5]。

direcpre=pos{1}-pre。

for i=2:100 %增加个体数目至100

if norm(pos{i}-pre)

direcpre=pos{i}-pre。

end

end

x=zeros(1,100)。

y=zeros(1,100)。

z=zeros(1,100)。

v=0.1。

direc1=cell(1,100)。

direc2=cell(1,100)。

direc3=cell(1,100)。

direc4=cell(1,100)。

direc5=cell(1,100)。

direc6=cell(1,100)。

direc7=cell(1,100)。

direc8=cell(1,100)。

cons=ones(1,100)。

for (i=1:100)

direc1{i}=[rand(1,1)*5,rand(1,1)*5,rand(1,1)*5]-pos{i}。 direc2{i}=zeros(1,3)。

direc3{i}=zeros(1,3)。

direc4{i}=zeros(1,3)。

direc7{i}=zeros(1,3)。

direc8{i}=zeros(1,3)。

end

for t=1:2000

pre=pre+direcpre/norm(direcpre)*1.1*v/10。

for i=1:100

m=0。

n=0。

k=0。

l=[0,0,0]。

direc5{i}=zeros(1,3)。

direc6{i}=zeros(1,3)。

for j=1:100

if j==i

continue

else

direc7{i}=(direc7{i}+pos{j}-pos{i})./j。

direc8{i}=(direc8{i}+direc1{j})./j。

if( norm(pos{i}-pos{j})<=5)&&(norm(pos{i}-pos{j})>=0.01) m=m+1。k=k+norm(direc1{j})。l=l+direc1{j}。

direc2{i}= (direc2{i}+(pos{j}-pos{i}))./m。

direc3{i}=(direc3{i}+direc1{j})./m。

cons(i)=norm(l)./k。

elseif norm(pos{i}-pos{j})<0.01

n=n+1。

direc4{i}=(direc4{i}+pos{i}-pos{j})./n。

end

end

end

if norm(pos{i}-pre)<=1 %定义危险区域

direc5{i}=pos{i}-pre。

direc6{i}=(-1)*direcpre。

direc1{i}=direc5{i}+direc6{i}。

pos{i}=pos{i}+(0.5*direc5{i}+0.5*direc6{i})*v。

if norm(pos{i}-pre)<=0.1

pos{i}=[-111,-11,-11]。direcpre=[rand(1,1)*10,rand(1,1)*10,rand(1,1)*10]。

end

for r=1:30

if norm(pos{r}-pre)

direcpre=pos{r}-pre。

end

end

else

if (pos{i}(1)>=5)||(pos{i}(2)>=5)||(pos{i}(1)<=0)||(pos{i}(2)<=0)||(pos{i}(3)>=5)||(pos{i}(3)<=0) direc1{i}=(-1)*direc1{i}。 pos{i}=pos{i}+direc1{i}/norm(direc1{i})*v。

else

if norm(pos{i}-pre)<=2.5 %定义感知区域

direc5{i}=pos{i}-pre。

direc6{i}=direcpre*(-1)。

direc1{i}=0.1*direc2{i}+0.1*direc3{i}+0.1*direc4{i}+0.25* direc5{i}+ 0.25*direc6{i}+0.1*direc7{i}+0.1*direc8{i}。

pos{i}=pos{i}+direc1{i}/norm(direc1{i})*v。

else

direc1{i}=0.1*direc1{i}+0.1*direc2{i}+0.1*direc3{i}+0.1*direc4{i}+0.3*direc7{i}+0.3*direc8{i}。

pos{i}=pos{i}+direc1{i}/norm(direc1{i})*v*exp(5*(cons(i)-1))。

end

end

end

end

if (pre(1)>=5)||(pre(2)>=5)||(pre(1)<=0)||(pre(2)<=0)||(pre(3)>=5)||(pre(3)<=0)

direcpre=(-1)*direcpre。

end

pre=pre+direcpre/norm(direcpre)*1.1*v/10。

for i=1:100

x(i)=pos{i}(1)。

y(i)=pos{i}(2)。

z(i)=pos{i}(3)。

end

plot3(x,y,z,'.',pre(1),pre(2),pre(3),'o')

axis([0 5 0 5 0 5],'square','manual')。

grid on。

pause(0.1)。

end

8.3、问题三地程序代码：

pos=cell(1,30)。

food=[5,5,5]。 %确定食物位置

for i=1:30

pos{i}=[rand(1,1)*10,rand(1,1)*10,rand(1,1)*10]

end

x=zeros(1,30)。

y=zeros(1,30)。

z=zeros(1,30)。

v=0.1。

direc1=cell(1,30)。

direc2=cell(1,30)。

direc3=cell(1,30)。

direc4=cell(1,30)。

direc8=cell(1,30)。

direc9=cell(1,30)。

cons=zeros(1,30)。

for i=1:30

direc1{i}=[rand(1,1)*10,rand(1,1)*10,rand(1,1)*10]-pos{i}。

direc2{i}=zeros(1,3)。

direc3{i}=zeros(1,3)。

direc4{i}=zeros(1,3)。

direc7{i}=zeros(1,3)。

direc8{i}=zeros(1,3)。

end

for t=1:2000

for i=1:10 %信息丰富者行为

n=0。

for j=1:30

if norm(pos{i}-pos{j})<0.1。

n=n+1。

direc4{i}=(direc4{i}+pos{i}-pos{j})./n。

end

end

if norm(pos{i}-food)>=0.1

pos{i}=pos{i}+(0.5*direc4{i}+0.5*(food-pos{i}))/norm(0.5*direc4{i}+0.5*(food-pos{i}))*v/8。

end

end

for i=11:30 %普通个体行为

if norm(pos{i}-food)<=3

direc1{i}=direc1{i}+0.5*(food-pos{i})。

end

m=0。n=0。k=0。l=[0,0,0]。

for j=1:30

if j==i

continue

else

for r=1:30

direc7{i}=(direc7{i}+pos{r}-pos{i})./r。

direc8{i}=(direc8{i}+direc1{r})./r。

end

if (norm(pos{i}-pos{j})<=5)&&(norm(pos{i}-pos{j})>=0.5)

m=m+1。k=k+norm(direc1{j})。l=l+direc1{j}。

direc2{i}= (direc2{i}+(pos{j}-pos{i}))./m。

direc3{i}=(direc3{i}+direc1{j})./m。

cons(i)=norm(l)/k。

elseif norm(pos{i}-pos{j})<0.5

n=n+1。

direc4{i}=(direc4{i}+pos{i}-pos{j})./n。

end

end

end

if

(pos{i}(1)>=10)||(pos{i}(2)>=10)||(pos{i}(1)<=0)||(pos{i}(2)<=0)||(pos{i}(3)>=10)||(pos{i}(3)<=0) direc1{i}=(-1)*direc1{i}。pos{i}=pos{i}+direc1{i}/norm(direc1{i})*v。

else

direc1{i}=0.05*direc1{i}+0.15*direc2{i}+0.15*direc3{i}+0.05*direc4{i}+0.05*direc7{i}+0.05*direc8{i}。

pos{i}=pos{i}+direc1{i}/norm(direc1{i})*v*exp(5*(cons(i)-1))。

end

end

for i=1:30

x(i)=pos{i}(1)。

y(i)=pos{i}(2)。

z(i)=pos{i}(3)。

end

plot3(x,y,z,'.',food(1),food(2),food(3),'o')

axis([0 10 0 10 0 10],'square','manual')

grid on。

pause(0.01) end