基于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 direc 1(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：捕食者的位置；
3.12、i pos ：第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
direc 1(t 1)*direc 1(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 β-+=
t=20s t=30s
这里
max v 0.1=，β为一可调参数，当0β=时，速度为原速度，当0β
，个体
的运动速度比原模型快，系统更易趋于同步。

这样，速率不仅具有改变下一时刻的作用，而且是携带信息的载体。

这种信息就是个体的一致性序数。

为了是所有的个体尽快达到同步，下一时刻速度方向的计算中，我们就利用这一信息，以加快收敛速度。

当个体超出活动范围时，它会向相反的方向运动，即：
i i
direc 1(t 1)direc 1(t)+=- 当在活动范围内时，对各个原则的影响设置权重，在此更多地考虑聚集性和排列性。

下一时刻的方向：
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
direc5 direc6
小鱼direc6 direc5 4.2、问题二
当个体和捕食者的距离较短时，该个体迅速逃逸，暂时不考虑对群体的影响。

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

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 direc 1(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(t1)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(t1)pos(t)*v
direc1(t)
+=+
运用MATLAB模拟，可得：
t=0s t=10s
t=20s t=30s
4.3、问题三
采用Lead-follower模型，假设鱼群中有一部分个体是信息丰富者,设置含有食物的环境,这一部分个体最先感知到食物的存在,并以一定速率接近食物.建立在第一问的基础上,由于个体之间有运动的联系,靠近信息丰富者的个体会感知到邻居的运动，这部分个体会首先向信息丰富者靠拢，并使自身的运动方向向信息丰富者的运动方向转移，使得有更多的个体向食物运动。

再有了多数个体向食物运动后，剩下的少部分离信息丰富较远的个体，由于必须向整体的中心和平均方向靠拢的原则，所以这部分个体会先向整体靠近，然后感知到向食物运动的邻居，接着使自己的位置和速度方向向食物转移。

对于信息丰富者，几乎不受群体的影响，只考虑排斥性，但是速度不能太快，否则，
就不能达到带动群体的效果，以最大速度的1/8做为信息丰富者的速度。

信息丰富者下一时刻的方向：
i i
i
direc1(t1)0.5
*direc4(t)0.5*(food pos(t))
+=+-
信息丰富者下一时刻的位置：
i
i i
i
direc1(t)
pos(t1)pos(t)0.125*v*
direc1(t)
+=+
当其他鱼接近食物时，它们的速度也会发生变化：
i i i i
i i i i direc1(t1)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*direc 7{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)<norm(direcpre)
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
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)<norm(direcpre)
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;
direc1{i}=0.1*direc1{i}+0.1*direc2{i}+0.1*direc3{i}+0.1*direc4{i}+0.3*direc
7{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*(foo d-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。