动物集群运动行为模型系列之七-毕业论文.doc

动物集群运动行为模型
动物集群运动模型问题
摘要
本文对于动物群体运动问题，建立了矢量方程模型。运用matlab编程对鱼群运动进行了仿真，得到了动物集群运动和躲避威胁等行为仿真结果。
问题一中，根据实际情况，制定了鱼运动的三条规则。然后将群体看做由粒子组成的集合，通过分析粒子受力，建立了矢量运动方程模型：

接着算出加速度矢量，进而求解运动轨迹。根据所列方程，利用matlab编程，对聚群运动进行了仿真，并绘制出鱼群环绕运动的稳定分析图。
对于问题二，根据鱼躲避捕食者的运动状态，建立了躲避运动的模型：

然后将鲨鱼运动分为开始接近鱼群到在鱼群中运动，最后离开鱼群等三个过程，细致分析了三个过程中鱼群的变化情况。将运动方程与分析相结合，利用matlab编程，得到较为理想的仿真结果。
问题三中，在分析信息丰富者对个体运动的影响时，在第一问的基础上，引入信息丰富者对个体的影响力。将信息影响力与其他作用力力矢量相加，得到个体运动影响力，然后计算个体加速，进而求解出运动轨迹。根据分析方程，得出信息丰富者会通过信息的传递，使群体跟随信息丰富着运动。

关键词：矢量；仿真；鱼群运动

一、问题重述
在动物界，大量集结成群进行移动或者觅食的例子并不少见，这种现象在食草动物、鸟、鱼和昆虫中都存在。这些动物群在运动过程中具有很明显的特征：群中的个体聚集性很强，运动方向、速度具有一致性。通过数学模型来模拟动物群的集群运动行为以及探索动物群中的信息传递机制一直是仿生学领域的一项重要内容。
请观察下面附件中给出的图片和视频资料，或者在网上搜索相关资料观察，思考动物集群运动的机理，建立数学模型刻画动物集群运动、躲避威胁等行为，例如，可以考虑以下问题的分析建模：
1. 建立数学模型模拟动物的集群运动。
2. 建立数学模型刻画鱼群躲避黑鳍礁鲨鱼的运动行为。
3. 假定动物群中有一部分个体是信息丰富者(如掌握食物源位置信息，掌握迁徙路线信息)，请建模分析它们对于群运动行为的影响，解释群运动方向决策如何达成。
二、模型假设
（1）假设所有的个体生理上不存在差异，并且遵循同样的准则。
（2）假设每个个体能够感知它在群体的位置（内部，或在群体的边沿）。
（3）假设不同相邻个体的相互作用力是累积的。
（4）假设一对个体间力的大小取决于两者间的距离和它们的相对速度。
（5）假设两个体间距离决定的力是一个平行于连接两者向量的向量，而速度决定力是一个平行于两者速度差向量的向量。
三、符号说明
符号
含义

相邻粒子距离

环绕半径

粒子位置矢量

粒子速度矢量

主动力

阻尼系数

环绕角速度

斜率

相互作用力

鱼最大感知距离

鱼最大告警距离
四、模型建立于求解
4.1问题一
4.1.1问题分析
通过分析视频资料，可以看出鱼群是一种自组织群体，没有固定的领导,但是群体往往呈现出运动方向有序、运动协调及集聚性等特征。这种自组织群体的典型特征及其在进化和生存方面的意义。对鸟类、昆虫、鱼群等自组织群体的观察研究结果表明，群体中的每一个体在遵循简单的行为规则条件时，就有可能出现有序协调的运动状态。
个体鱼的运动行为都具有一定的特性，具体情况如下图1所示：

图1 行为分析图
（1）当时，由于距离过近，鱼之间会产生排斥作用。
（2）当时，距离适中，鱼的运动状态会与它相邻鱼相协调。
（3）当时，距离较远，鱼之间会相互吸引。
（4）当时，距离过远，超出了鱼的感知范围，鱼之间不产生影响。
4.1.2模型建立
根据上面的运动特点,将鱼群看做粒子集合，那么其中的个体一个粒子，粒子的运动影响分为主动因素和被动因素，据此建立鱼群运动的数学模型：
建立一个有n个粒子的群体，每个粒子编号为i=1,2,…,n。分别用表示位置，用表示速度。每个粒子有一个前端和后端，并且让粒子i的身体走向与运动方向（即速度的方向）相同。粒子i的运动遵循牛顿定律，可得到模型:

这里取粒子质量。代表主动力，由粒子自身产生，取决于环境影响和粒子在群体中的位置。是阻尼力，其中系数保证了速度有界。式（2）表示如果一个粒子停止推进，那么它的速度会以比率减小到零。是指相邻的粒子对它的作用力。对于t时刻有，

这里a，b，c代表三个坐标轴，将粒子i的相邻粒子j定义为此范围内的所有粒子。
通常是常向量，根据假设，相互作用力由下式给出：

其中；。和分别是取决于位置和速度的力。j代表所有影响粒子i运动的粒子。与粒子和相邻粒子位置的矢量差有关；与粒子和相邻粒子速度的矢量差有关。
为产生一个非零的间隔距离，通常对于较大的x为正值，较小的x值为负值，表示短程排斥和长程吸引。指粒子j对粒子i产生一个吸引力（排斥力）；而指力使粒子i的速度趋向(偏离)于粒子j的速度。需要注意的是，组成相互作用力的位置影响的力和速度影响力，会指向不同的方向，并且具有不同大小。

图2 模型向量二维示意图
4.1.3模型求解：
可将所建立的模型，应用于三维空间。建立三维坐标系，分别代表三个坐标轴上的单位向量，那么可将向量用坐标表示，即（a，b，c），向量的运算就可转化为坐标的运算，本文就二维空间的运动进行分析求解。
(1)聚集的运动
聚群行为是鱼类较常见的一种现象，大量或少量的鱼都能聚集成群，这是它们在进化过程中形成的一种生存方式，可以进行集体觅食和躲避敌害。这里聚集运动指一些离散的粒子，通过运动逐渐聚成整体，然后共同向前运动。在模型的基础上，根据查找资料得到：
主动力是相同的取，阻尼系数不变，取。，。由于粒子是离散的，因此不考虑排斥力和偏离力的作用的作用，即，。通过计算加速度的变化，以及矢量方向的变化，来计算速度和位移矢量的变化，进而确定粒子的运动轨迹。

图3 粒子运动分析图
每0.1s计算各点的位置，可以得到位置与速度变化的关系式：

每个粒子加速度

这里用matlab在[0,4] [0,4]范围内,内随机生成5个点，计算它们之间的相互作用力，进而求解加速度，然后计算得到位置变化。通过matlab编程绘图，每隔0.1s，绘制出点的位置，经过1min后，得到5个运动路径的散点图4，通过观察散点图，可以直观的观察出聚集运动的规律。

图4 运动散点图
简化相互作用力的函数方程式，对聚群运动进行仿真，得到结果为：

图5 仿真结果（左图为初始状态，右图为程序运行结果）
（2）鱼群环绕运动
鱼群在聚集后，往往会绕一中心轴转动，下面对鱼群的环绕运动，在二维空间进行分析。
设n个粒子，以半径r0， 绕圆心O持续匀速转动，角速度是。粒子间等间距d，在运动过程中，粒子i跟随粒子i+1，粒子n跟随1。环绕中，相邻粒子扇形角度为，其中n是粒子的数量（见图5）。
得到下列式子：

图6 环绕模型示意图
黑色箭头代表运动方向（与圆相切），而灰色箭头代表群体力的方向。

等式（7）-（8）将第n个粒子与第一个粒子连成一个环。
然后将相互作用力分解为与圆相切的力和指向圆心的力。

图7 力的分解
在环绕模型中，没有线性加速度，所以与圆相切的力是平衡的，因此，对于每个粒子有：

质量为1的粒子绕圆运动的向心力为：

其中是单位径向矢量， ，带入得到

这里取两粒子连线方向与切线方向夹角为和粒子连线方向与圆中心连线夹角为（见图7）

图8 环绕运动分析
根据三角函数关系得到：

将切向力与径向力方程带入得到：

解得：

角速度为：

通过查找资料，赋予g（x），n和值。对于n个粒子的系统，g（x），n和值，以及环绕角速度和切向速度，能过完全表征环绕结果的特点。结合式（9）和式（10），求解出存在条件：

其中。稳定的环绕运动运动只有对于给定的g（x）在d值满足式（11）时存在。
根据查找资料，取所得的图像。横轴代表粒子间距离d。纵轴代表距离力大小。。图8中取b=1.5，两曲线存在交点，；而在图9中取b=2，两曲线不存在交点。

图9 b=1.5时函数图像

图10 b=2时函数图像
g（d）与直线sd的交点是环绕运动点。在图8中，交点是环绕运动存在的d值；图9中没有交点，表明环绕运动不存在。坡度s影响交点是否存在，小的坡度增加了相交的可能性。所以减小阻尼系数或增大n值，增大了环绕运动稳定存在可能性。当一个或多个粒子离开群体后，会打破磨盘运动的平衡。图10表明对于n的不同值的情况。在这个结果中，少于五个粒子，环绕运动结果不存在。

图11 不同n值对应图像结果
4.2问题二模型建立与求解
问题二中，鱼群通常会聚集成群，以避免被捕食者单独捕捉。当鱼群遭遇黑鳍鲨鱼后，会表现出躲避逃逸行为。一部分鱼发现鲨鱼后，会发生躲避和逃逸，并将信息传递给附近的其他鱼，进而引起其他鱼的逃逸。鱼在逃逸过程中，一方面身体会旋转一定的角度，改变自身速度矢量的方向，与鲨鱼速度矢量的方向的夹角，使改变后的方向能够让自身尽快的逃离捕食者的追击。另一方面，在逃逸过程中，鱼的运动要避免与同伴碰撞，这就限制了与身体角度的变化，以及速度的增加。据此建立一下模型：

其中代表鱼转动的角度。在逃逸过程中，个体鱼会因与相邻鱼距离过近而产生排斥反应，因此相邻鱼相互吸引的力可以忽略，那么可以得到：

其中是单位方向向量，方向与两个体连线方向相同，指向受力一方。
将鲨鱼简化为一点，假设鱼的感知最大半径为，鱼警告信号发送最大距离为。当鲨鱼进入到鱼感知范围内时，鱼就立即发生逃逸运动，并同时向周围鱼发出警告信号。周围有些鱼虽然没有感知到鲨鱼的存在，但当它感知到警告信号后，同样立即进行逃逸。由于每条鱼在鱼群中的位置不同，因此鱼在逃逸过程中，速度矢量的变化量各不相同。
距离鲨鱼最近的一些鱼，。根据力矢量可求出鱼的加速度矢量，进而得出鱼的逃逸运动速度矢量，进而求解运动轨迹。综合分析鱼群中各个鱼的运动轨迹，可以得到鱼群躲避捕食者的运动情况。
下面分析鲨鱼进入鱼群的情况，这里定义距离鲨鱼最近的鱼为鱼群的前方，最远的地方为后方：
（1）鲨鱼刚开始被鱼群感知时，在鱼群前方的鱼感知到危险后开始后退，并向周围发出告警信号。在前方的鱼和接受到告警信号的鱼，发现鲨鱼后，产生主动力的作用立即逃逸，而逃逸时，身后的鱼会对其产生排斥力的作用，根据模型将两矢量相加，得到了鱼逃逸时力方向矢量。位于鱼群后方的鱼既没有发现鲨鱼，由于距离较远也没有感知到警告信号，因此运动状态没有发生变化，如图11所示：

图12 开始进入示意图
（2）鲨鱼进入到鱼群中后，大部分鱼通过感知鲨鱼或接收告警信号，已经得知了鲨鱼的大体方位，并开始整体向远离鲨鱼的方向运动，原来位于前方的鱼大部分撤离到安全区域，绕到了鲨鱼后方的位置，并仍未与鱼群脱离。前方的鱼撤离后，原来处于中间位置的鱼变成了前方的鱼，它们继续绕鲨鱼做逃避运动。此时，处于中间位置鱼群在躲避过程中，鱼与鱼之间的距离逐渐压缩，且鲨鱼头与鱼群距离逐渐缩小，鲨鱼尾与鱼群距离逐渐增大，具体情况如下：

图13 鱼群中运动示意图
（3）鲨鱼离开鱼群后，鱼群会聚拢，再次形成一个同一整体，并整体向远离鲨鱼的方向运动。

图14 离开后鱼群运动示意图
根据模型中给出的鱼受力方程，结合鱼群运动分析，可以利用matlab编程仿真，得到结果如下图：

图15 鲨鱼与鱼群运动仿真结果图
4.3问题三
任何生物都不是孤立地生活在自然界中，它们总是组成一个小的生活群体，若动物群中有一部分个体是信息丰富者(如掌握食物源位置信息，掌握迁徙路线信息)，这部分个体会将信息传递给其他的同伴,使种群尽快地找到食物,若没有这部分信息丰富者,则种群在觅食、迁徙过程中对路线的选择就具有一定的盲目性，信息丰富者在种群中扮演一个决策者的作用，帮助种群找到一个最佳的路线或决策方案。动物传递信息的方式多种多样，主要有以下几种：
1.视觉通讯
视觉通讯的形式是比较广泛的，雄性驯鹿头上硕大的犄角，草原上雄性狮子颈部漂亮的长鬃毛，这些动物的外表特征都是向雌性同类发出的视觉信号。青蛙在草丛中呈现碧绿的体色，而潮一、的保护色往往是通过散布错误的视觉信息来迷惑天敌或猎物的。
2.听觉通讯
鸟类为吸引异性排斥同性，宣告领地占有的歌声以及警告捕食者来到的尖叫声都是听觉通讯。动物世界里有一些动物是依靠超声波来进行通讯与捕食的，如人们熟悉的编幅和海豚，就是利用超声波通讯的。
3.化学通讯
化学通讯就是动物通过释放一些化学物质来影响或控制其他动物的行为。化学通讯有时会影响整个动物群体的活动甚至调节整个种群。这些化学物质称为外激素。
4.触觉通讯
触觉通讯也是一种相当普遍的通讯方式。对于视觉能力有限或者生活在无法利用视觉通讯环境中的动物来说，触觉通讯往往是一种重要的传递信息的方式。某些生活在深海区域中的鱼类，由于光线很弱，视力退化了，但它们往往具有非常发达的鳍刺和触须，上面布满了敏感的神经，在水中游动时，它们可以感知水流的变化，寻觅与捕捉猎物和接收性信号。
5.电通讯
即电鱼、美洲鳗等动物所采用的电通讯方式。电信号通讯不受障碍的阻挡，具有高度的方向性，不过作用距离短，这一点和触觉信号通讯相似。
通讯过程中个体向其他个体发出信号，为其他个体的感觉器官所接受。信号不仅传递情报（信息），还有让对方改变行为的意义。
在个体受力分为以下几个方面：
（1）临近个体的平均吸引力
（2）临近个体的平均排斥力
（3）向附近临近个体运动状态向协调的平均协调力
（4）信息平均影响力
力，，可根据问题一种的模型，利用矢量运算方法求解。对于力，查找相关资料得知，力的方向与粒子和信息丰富者的连线平行，并指向信息丰富着，力的大小为：

其中，代表个体与信息丰富者的距离；代表群体系数，与群体大小和物种类别有关；代表其他个体信息影响。在信息传递过程中，周围的相邻个体可能较早得知信息，然后将信息传递出去，还可能个体本身通过其他渠道获得信息，
综合这些因素可以得到，较小，一般可以忽略它的影响。
利用矢量相加的的方法，求取个体受力为：

当力相比其他几个力较大时，力的方向就趋向于力的方向。个体运动具体表现为，个体趋向于信息丰富者运动。在群体运动过程中，表现为一些信息丰富者运动在群体的前方，“领导”着群体运动。
当一些单位获取到信息后，它们一方面向目标运动，另一方面向周围同伴传递信息，同伴接收到信息后，向信息丰富者和目标运动，同时发送信息，让更多的个体得知信息。通过信息的传递与趋向运动，进而形成了群体的运动。
五、模型评价
5.1模型优点
（1）从分析受力的角度，建立了矢量模型。定量的计算了鱼的运动状态，结果准确可靠。
（2）将鱼群中鱼当做有前后端的粒子，简化了问题，减少了计算量。
（3）引入了空间坐标，利用坐标进行运算，使得运算更加便捷，结果更加准确。
（3）分析鱼群运动较为细致全面，仿真结果较为准确。
5.2模型缺点
（1）方程数量较多，使仿真程序较复杂。
（2）查找的数据有限，结果可能会存在一定的误差。
参考文献
[1] 柳玲飞,周应祺.红鼻鱼群体结构的数学建模与仿真可视化，上海海洋大学海洋科学学院,2012.12.
[2] 程代展,胨翰馥.从群集到社会行为控制[J],科技导报, 2004.8.
[3] 赵建，曾建潮.鱼群集群行为的建模与仿真[J],太原科技大学学报.
[4] 肖人彬,陶振武.群集智能研究进展[J].管理科学学报, 2007.10
[5] 郑毅,吴斌.由鸟群和蚂蚁想到的—基于主体的仿真与群集智能的研究[J].微电脑世界,2001.1.

附录
程序
问题一：
%仿真
clear;clc
t=0.5;n=100;
x(1,:)=40*rand(1,n);y(1,:)=40*rand(1,n);
draction(1,:)=(rand(1,n)-0.5)*2*pi;
d=[];
for k=1:100
%%距离
for i=1:n
for j=1:n
d(i,j)=sqrt((x(k,i)-x(k,j))^2+(y(k,i)-y(k,j))^2);
end
end
for i=1:n
d(i,i)=inf;
end
%%速度方向
for i=1:n
A=0;B=0.1;
for j=1:n
if d(i,j)<0.3
A=A-draction(k,j)/d(i,j);
B=B-1/d(i,j);
if A~=0
draction(k+1,i)=(1-t)*draction(k,i)-t*A/B;
else draction(k+1,i)=draction(k,i);
end
elseif d(i,j)<5
A=A-draction(k,j)/d(i,j);
B=B-1/d(i,j);

if A~=0
draction(k+1,i)=(1-t)*draction(k,i)+t*A/B;
else draction(k+1,i)=draction(k,i);
end
end
end

v(k+1,i)=0.5;%速度大小的更新
end
%%坐标的更新
for i=1:n
x(k+1,i)=x(k,i)+v(k+1,i)*cos(draction(k+1,i))*1;
y(k+1,i)=y(k,i)+v(k+1,i)*sin(draction(k+1,i))*1;
if x(k+1,i)>40
x(k+1,i)=x(k+1,i)-40;
elseif x(k+1,i)<0
x(k+1,i)=x(k+1,i)+40;
end
if y(k+1,i)>40
y(k+1,i)=y(k+1,i)-40;
elseif y(k+1,i)<0
y(k+1,i)=y(k+1,i)+40;
end
end

end

for i=1:k
pause(0.2)
plot(x(i,:),y(i,:),'*')
axis([0 40 0 40])
getframe;
end
问题二：
clear
clc
t=0.5;n=300;
x=[];y=[];xs=[];
x(1,:)=40*rand(1,n);y(1,:)=40*rand(1,n);
draction(1,:)=(rand(1,n)-0.5)*2*pi;
d=[];
%初始化鲨鱼坐标
xs(1,:)=40*rand(1,1);ys(1,:)=40*rand(1,1);
dractions(1,1)=(rand(1,1)-0.5)*2*pi;
for k=1:200
dractions(k+1,1)=dractions(k,1)+(rand(1,1)-0.5);
vs(k+1,1)=0.7;
xs(k+1,1)=xs(k,1)+vs(k+1,1)*cos(dractions(k+1,1))*1;
ys(k+1,1)=ys(k,1)+vs(k+1,1)*sin(dractions(k+1,1))*1;
if xs(k+1,1)>40
xs(k+1,1)=xs(k+1,1)-40;
elseif xs(k+1,1)<0
xs(k+1,1)=xs(k+1,1)+40;
end
if ys(k+1,1)>40
ys(k+1,1)=ys(k+1,1)-40;
elseif ys(k+1,1)<0
ys(k+1,1)=ys(k+1,1)+40;
end

for i=1:n
ds(i)=sqrt((x(k,i)-xs(k,1))^2+(y(k,i)-ys(k,1))^2);
for j=1:n
d(i,j)=sqrt((x(k,i)-x(k,j))^2+(y(k,i)-y(k,j))^2);
end
end
for i=1:n
d(i,i)=inf;
end
%%方向
for i=1:n
A=0;B=0.1;
for j=1:n
if d(i,j)<0.2
A=A-draction(k,j)/3;
B=B-1/d(i,j);
if A~=0
draction(k+1,i)=(1-t)*draction(k,i)-t*A/B;
else draction(k+1,i)=draction(k,i);
end
elseif d(i,j)<5
A=A-draction(k,j)/d(i,j);
B=B-1/d(i,j);
if A~=0
draction(k+1,i)=(1-t)*draction(k,i)+t*A/B;
else draction(k+1,i)=draction(k,i);
end
end
end
if ds(i)<5
if x(k,i)>xs(k,1);
draction(k+1,i)=atan((y(k,i)-ys(k,1))./(x(k,i)-xs(k,1)));
end
if x(k,i)<xs(k,1);
draction(k+1,i)=atan((y(k,i)-ys(k,1))./(x(k,i)-xs(k,1)))+pi;
end
end
v(k+1,i)=0.5;

end
%%坐标的更新
for i=1:n
x(k+1,i)=x(k,i)+v(k+1,i)*cos(draction(k+1,i))*1;
y(k+1,i)=y(k,i)+v(k+1,i)*sin(draction(k+1,i))*1;
if x(k+1,i)>40
x(k+1,i)=x(k+1,i)-40;
elseif x(k+1,i)<0
x(k+1,i)=x(k+1,i)+40;
end
if y(k+1,i)>40
y(k+1,i)=y(k+1,i)-40;
elseif y(k+1,i)<0
y(k+1,i)=y(k+1,i)+40;
end
end

end

for i=1:k
pause(0.2)
plot(x(i,:),y(i,:),'.',xs(i,:),ys(i,:),'ro')
axis([0 40 0 40])
getframe;
end
20