动物集群运动行为模型系列之七

动物集群运动模型问题
摘要
本文对于动物群体运动问题，建立了矢量方程模型.运用matlab 编程对鱼群运动进行了仿真，得到了动物集群运动和躲避威胁等行为仿真结果.问题一中，根据实际情况，制定了鱼运动地三条规则.然后将群体看做由粒子组成地集合，通过分析粒子受力，建立了矢量运动方程模型：i i x v = i i i i i mv F v f γ=-+
接着算出加速度矢量，进而求解运动轨迹.根据所列方程，利用matlab 编程，对聚群运
动进行了仿真，并绘制出鱼群环绕运动地稳定分析图.对于问题二，根据鱼躲避捕食者地运动状态，建立了躲避运动地模型：
()()()()e i i i i i i ij j i
i dv t v t e t v t m m f dt T =-=+∑
()i
i dx v t dt
=
然后将鲨鱼运动分为开始接近鱼群到在鱼群中运动，最后离开鱼群等三个过程，细致分析了三个过程中鱼群地变化情况.将运动方程与分析相结合，利用matlab 编程，得到较为理想地仿真结果.问题三中，在分析信息丰富者对个体运动地影响时，在第一问地基础上，引入信息丰富者对个体地影响力.将信息影响力与其他作用力力矢量相加，得到个体运动影响力，然后计算个体加速，进而求解出运动轨迹.根据分析方程，得出信息丰富者会通过信息地传递，使群体跟随信息丰富着运动.
关键词：矢量；仿真；鱼群运动
一、问题重述
在动物界，大量集结成群进行移动或者觅食地例子并不少见，这种现象在食草动物、
鸟、鱼和昆虫中都存在.这些动物群在运动过程中具有很明显地特征：群中地个体聚集性很强，运动方向、速度具有一致性.通过数学模型来模拟动物群地集群运动行为以及探索动物群中地信息传递机制一直是仿生学领域地一项重要内容.
请观察下面附件中给出地图片和视频资料，或者在网上搜索相关资料观察，思考动物集群运动地机理，建立数学模型刻画动物集群运动、躲避威胁等行为，例如，可以考虑以下问题地分析建模：
1. 建立数学模型模拟动物地集群运动.
2. 建立数学模型刻画鱼群躲避黑鳍礁鲨鱼地运动行为.
3. 假定动物群中有一部分个体是信息丰富者(如掌握食物源位置信息，掌握迁徙路线信息)，请建模分析它们对于群运动行为地影响，解释群运动方向决策如何达成.
二、模型假设
（1）假设所有地个体生理上不存在差异，并且遵循同样地准则.
（2）假设每个个体能够感知它在群体地位置（内部，或在群体地边沿）.
（3）假设不同相邻个体地相互作用力是累积地.
（4）假设一对个体间力地大小取决于两者间地距离和它们地相对速度.
（5）假设两个体间距离决定地力是一个平行于连接两者向量地向量，而速度决定力是一个平行于两者速度差向量地向量.
三、符号说明
v
i
F
i
γ
f d 四、模型建立于求解
4.1问题一 4.1.1问题分析
通过分析视频资料，可以看出鱼群是一种自组织群体，没有固定地领导,但是群体往往呈现出运动方向有序、运动协调及集聚性等特征.这种自组织群体地典型特征及其在进化和生存方面地意义.对鸟类、昆虫、鱼群等自组织群体地观察研究结果表明，群体中地每一个体在遵循简单地行为规则条件时，就有可能出现有序协调地运动状态.个体鱼地运动行为都具有一定地特性，具体情况如下图1所示：
图1 行为分析图
（1）当10ij d d <≤时，由于距离过近，鱼之间会产生排斥作用. （2）当12ij d d d <≤时，距离适中，鱼地运动状态会与它相邻鱼相协调. （3）当23ij d d d <≤时，距离较远，鱼之间会相互吸引.
（4）当3ij d d >时，距离过远，超出了鱼地感知范围，鱼之间不产生影响. 4.1.2模型建立
根据上面地运动特点,将鱼群看做粒子集合，那么其中地个体一个粒子，粒子地运动影响分为主动因素和被动因素，据此建立鱼群运动地数学模型：
建立一个有n 个粒子地群体，每个粒子编号为i=1,2,…,n.分别用i x 表示位置，用i v 表示速度.每个粒子有一个前端和后端，并且让粒子i 地身体走向与运动方向（即速度i v 地方向）相同.粒子i 地运动遵循牛顿定律，可得到模型: (1)i i x v = (2)i i i i i mv F v f γ=-+
这里取粒子质量1i m =.i F 代表主动力，由粒子自身产生，取决于环境影响和粒子在群体中地位置.(0)v γγ>是阻尼力，其中系数γ保证了速度有界.式（2）表示如果一个粒子停止推进，那么它地速度会以比率γ减小到零.i f 是指相邻地粒子对它地作用力.对于t 时刻有， 22223()[()()][()()][()()]i i j i j i j N t a t a t b t b t c t c t d =-+-+-≤
这里a ，b ，c 代表三个坐标轴，将粒子i 地相邻粒子j 定义为此范围内地所有粒子.
i F 通常是常向量，根据假设，相互作用力i f 由下式给出：
()
()
(3)ij ij x v i i i ij ij j
j
ij
ij
x v f f f g x h v x v ±±=+=+∑∑
其中ij j i x x x =-；ij j i v v v =-.x i f 和v i f 分别是取决于位置和速度地力.j 代表所有影响粒子i 运动地粒子.x i f 与粒子和相邻粒子位置地矢量差有关；v i f 与粒子和相邻粒子速度地矢量差有关.为产生一个非零地间隔距离，()g x ±通常对于较大地x 为正值，较小地x 值为负值，表示短程排斥和长程吸引.()0(0)g x +><指粒子j 对粒子i 产生一个吸引力（排斥力）；而()0(0)h x +><指力v i f 使粒子i 地速度趋向(偏离)于粒子j 地速度.需要注意地是，组成相互作用力地位置影响地力和速度影响力，会指向不同地方向，并且具有不同大小.
图2 模型向量二维示意图
4.1.3模型求解：
可将所建立地模型，应用于三维空间.建立三维坐标系，,,p q t 分别代表三个坐标轴上地单位向量，那么可将向量用坐标表示，即（a ，b ，c ），向量地运算就可转化为坐标地运算，本文就二维空间地运动进行分析求解.(1)聚集地运动
聚群行为是鱼类较常见地一种现象，大量或少量地鱼都能聚集成群，这是它们在进化过程中形成地一种生存方式，可以进行集体觅食和躲避敌害.这里聚集运动指一些离散地粒子，通过运动逐渐聚成整体，然后共同向前运动.在模型地基础上，根据查找资料得到：主动力是相同地取0.1i F =，阻尼系数γ不变，取
0.5γ=.()()50
1
()0.52x x g x e
e
-
-+=-，(),(0.56)h v v αα+==.由于粒子是离散地，因此不考
虑排斥力和偏离力地作用地作用，即()0g x -=，()0h x -=.通过计算加速度地变化，以及矢量方向地变化，来计算速度和位移矢量地变化，进而确定粒子地运动轨迹.
图3 粒子运动分析图
每0.1s 计算各点地位置，可以得到位置与速度变化地关系式：
()()(),(0.1) (4)i i i X t X t V t τττ+=+=
每个粒子加速度
()
()
50
1
0.10.5[0.52]
(0.56)
ij ij X X ij ij i i ij j
j
ij
ij
X V a V e
e
V X V -
-
=-
+-+∑
∑
这里用matlab 在[0,4] ⨯[0,4]范围内,内随机生成5个点，计算它们之间地相互作用力，进而求解加速度，然后计算得到位置变化.通过matlab 编程绘图，每隔0.1s ，绘制出点地位置，经过1min 后，得到5个运动路径地散点图4，通过观察散点图，可以直观地观察出聚集运动地规律.
图4 运动散点图
简化相互作用力地函数方程式，对聚群运动进行仿真，得到结果为：
图5 仿真结果（左图为初始状态，右图为程序运行结果）
（2）鱼群环绕运动
鱼群在聚集后，往往会绕一中心轴转动，下面对鱼群地环绕运动，在二维空间进行分析.
设n 个粒子，以半径r0， 绕圆心O 持续匀速转动，角速度是0ω.粒子间等间距d ，在运动过程中，粒子i 跟随粒子i+1，粒子n 跟随 1.环绕中，相邻粒子扇形角度为2/n θπ=，其中n 是粒子地数量（见图5）.
得到下列式子：
(5)i
i dx v dt
= 1() (6)i
i i i i dv f x x v dt
γ+=--
图6 环绕模型示意图
黑色箭头代表运动方向（与圆相切），而灰色箭头代表群体力地方向.
(7)n
n dx v dt
= 1() (8)n
n n n dv f x x v dt
γ=-- 等式（7）-（8）将第n 个粒子与第一个粒子连成一个环. 然后将相互作用力分解为与圆相切地力t f 和指向圆心地力c f .
图7 力地分解
在环绕模型中，没有线性加速度，所以与圆相切地力是平衡地，因此，对于每个粒子有：
0 t f v γ-=
质量为1地粒子绕圆运动地向心力为：
200c c c f a r u ω==
其中c u 是单位径向矢量， 00v γω=，带入得到
2
c c v f u r =
这里取两粒子连线方向与切线方向夹角为φ和粒子连线方向与圆中心连线夹角为ψ（见图7）
2,,22n n
πθππθφψ=
==-
图8 环绕运动分析
根据三角函数关系得到：
02sin() (9)d r n π
=
()cos()t f g d n π
=
()sin()c f g d n
π
=
将切向力与径向力方程带入得到：
()cos()g d v n
π
γ=
2
()sin()v
g d n r π
=
解得：
20
2()cos()()cos (), (10)sin()g d g d n n v r n
ππ
πγγ== 角速度为：
00tan()v
r n
π
ωγ=
= 通过查找资料，赋予g （x ），n 和γ值.对于n 个粒子地系统，g （x ），n 和γ值，以及环绕角速度和切向速度，能过完全表征环绕结果地特点.结合式（9）和式（10），求解
出存在条件：() (11)g d sd =
其中2
22cos ()
s n
γπ
=
.稳定地环绕运动运动只有对于给定地g （x ）在d 值满足式（11）
时存在.
根据查找资料，取()()0.5,5,()x x a
b
n g x Ae
Be
γ--===-所得地图像.横轴代表粒子间
距离d.纵轴代表距离力大小. 1.5,10,3A a B ===.图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问题二模型建立与求解
问题二中，鱼群通常会聚集成群，以避免被捕食者单独捕捉.当鱼群遭遇黑鳍鲨鱼后，会表现出躲避逃逸行为.一部分鱼发现鲨鱼后，会发生躲避和逃逸，并将信息传递给附近地其他鱼，进而引起其他鱼地逃逸.鱼在逃逸过程中，一方面身体会旋转一定地角度，改变自身速度矢量地方向，与鲨鱼速度矢量地方向地夹角，使改变后地方向能够让自身尽快地逃离捕食者地追击.另一方面，在逃逸过程中，鱼地运动要避免与同伴碰撞，这就限制了与身体角度地变化，以及速度地增加.据此建立一下模型：
()()()()e i i i i i i ij j i
i dv t v t e t v t m m f dt T =-=+∑
()i
i dx v t dt
=
其中()i e t 代表鱼转动地角度.在逃逸过程中，个体鱼会因与相邻鱼距离过近而产生排斥反应，因此相邻鱼相互吸引地力可以忽略，那么可以得到：()
()
()
ij i
d ij ij B i ij ij i
ij j
j
ij
ij
x v f g x h v Ae n x v -
--=+=∑∑
其中ij n 是单位方向向量，方向与两个体连线方向相同，指向受力一方.
将鲨鱼简化为一点，假设鱼地感知最大半径为f d ，鱼警告信号发送最大距离为w d .当鲨鱼进入到鱼感知范围f d 内时，鱼就立即发生逃逸运动，并同时向周围鱼发出警告信号.周围有些鱼虽然没有感知到鲨鱼地存在，但当它感知到警告信号后，同样立即进行逃逸.由
于每条鱼在鱼群中地位置不同，因此鱼在逃逸过程中，速度矢量地变化量各不相同.距离鲨鱼最近地一些鱼，.根据力矢量可求出鱼地加速度矢量，进而得出鱼地逃逸运动速度矢量，进而求解运动轨迹.综合分析鱼群中各个鱼地运动轨迹，可以得到鱼群躲避捕食者地运动情况.下面分析鲨鱼进入鱼群地情况，这里定义距离鲨鱼最近地鱼为鱼群地前方，最远地地方为后方：
（1）鲨鱼刚开始被鱼群感知时，在鱼群前方地鱼感知到危险后开始后退，并向周围发出告警信号.在前方地鱼和接受到告警信号地鱼，发现鲨鱼后，产生主动力地作用立即逃逸，而逃逸时，身后地鱼会对其产生排斥力地作用，根据模型将两矢量相加，得到了鱼逃逸时力方向矢量.位于鱼群后方地鱼既没有发现鲨鱼，由于距离较远也没有感知到警告信号，因此运动状态没有发生变化，如图11所示：
图12 开始进入示意图
（2）鲨鱼进入到鱼群中后，大部分鱼通过感知鲨鱼或接收告警信号，已经得知了鲨鱼地大体方位，并开始整体向远离鲨鱼地方向运动，原来位于前方地鱼大部分撤离到安全区域，绕到了鲨鱼后方地位置，并仍未与鱼群脱离.前方地鱼撤离后，原来处于中间位置地鱼变成了前方地鱼，它们继续绕鲨鱼做逃避运动.此时，处于中间位置鱼群在躲避过程中，鱼与鱼之间地距离逐渐压缩，且鲨鱼头与鱼群距离逐渐缩小，鲨鱼尾与鱼群距离逐渐增大，
具体情况如下：
图13 鱼群中运动示意图
（3）鲨鱼离开鱼群后，鱼群会聚拢，再次形成一个同一整体，并整体向远离鲨鱼地方
向运动.
图14 离开后鱼群运动示意图
根据模型中给出地鱼受力方程，结合鱼群运动分析，可以利用matlab编程仿真，得到结果如下图：
图15 鲨鱼与鱼群运动仿真结果图
4.3问题三
任何生物都不是孤立地生活在自然界中，它们总是组成一个小地生活群体，若动物群中有一部分个体是信息丰富者(如掌握食物源位置信息，掌握迁徙路线信息)，这部分个体会将信息传递给其他地同伴,使种群尽快地找到食物,若没有这部分信息丰富者,则种群在觅食、迁徙过程中对路线地选择就具有一定地盲目性，信息丰富者在种群中扮演一个决策者地作用，帮助种群找到一个最佳地路线或决策方案.动物传递信息地方式多种多样，主要有以下几种：
1.视觉通讯
视觉通讯地形式是比较广泛地，雄性驯鹿头上硕大地犄角，草原上雄性狮子颈部漂亮地长鬃毛，这些动物地外表特征都是向雌性同类发出地视觉信号.青蛙在草丛中呈现碧绿地体色，而潮一、地保护色往往是通过散布错误地视觉信息来迷惑天敌或猎物地.2.听觉通讯
鸟类为吸引异性排斥同性，宣告领地占有地歌声以及警告捕食者来到地尖叫声都是听觉通讯.动物世界里有一些动物是依靠超声波来进行通讯与捕食地，如人们熟悉地编幅和海豚，就是利用超声波通讯地.3.化学通讯
化学通讯就是动物通过释放一些化学物质来影响或控制其他动物地行为.化学通讯有时会影响整个动物群体地活动甚至调节整个种群.这些化学物质称为外激素.4.触觉通讯
触觉通讯也是一种相当普遍地通讯方式.对于视觉能力有限或者生活在无法利用视觉通讯环境中地动物来说，触觉通讯往往是一种重要地传递信息地方式.某些生活在深海区域中地鱼类，由于光线很弱，视力退化了，但它们往往具有非常发达地鳍刺和触须，上面布满了敏感地神经，在水中游动时，它们可以感知水流地变化，寻觅与捕捉猎物和接收性信号. 5.电通讯即电鱼、美洲鳗等动物所采用地电通讯方式.电信号通讯不受障碍地阻挡，具有高度地方向性，不过作用距离短，这一点和触觉信号通讯相似.通讯过程中个体向其他个体发出信号，为其他个体地感觉器官所接受.信号不仅传递情报（信息），还有让对方改变行为地意义.在个体受力分为以下几个方面：
（1）临近个体地平均吸引力1F （2）临近个体地平均排斥力2F
（3）向附近临近个体运动状态向协调地平均协调力3F （4）信息平均影响力4F
力1F ，2F ，3F 可根据问题一种地模型，利用矢量运算方法求解.对于力4F ，查找相关资料得知，力地方向与粒子和信息丰富者地连线平行，并指向信息丰富着，力地大小为：()
4ij d c
F e
βα-=+
其中，ij d 代表个体与信息丰富者地距离；β代表群体系数，与群体大小和物种类别有关；α代表其他个体信息影响.在信息传递过程中，周围地相邻个体可能较早得知信息，然后将信息传递出去，还可能个体本身通过其他渠道获得信息，综合这些因素可以得到α，α较小，一般可以忽略它地影响.
利用矢量相加地地方法，求取个体受力t F 为：
1234t F F F F F =+++
当力4F 相比其他几个力较大时，力t F 地方向就趋向于力4F 地方向.个体运动具体表现为，个体趋向于信息丰富者运动.在群体运动过程中，表现为一些信息丰富者运动在群体地前方，“领导”着群体运动.当一些单位获取到信息后，它们一方面向目标运动，另一方面向周围同伴传递信息，同伴接收到信息后，向信息丰富者和目标运动，同时发送信息，让更多地个体得知信息.通过信息地传递与趋向运动，进而形成了群体地运动.五、模型评价
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。