基于智能体和元胞自动机的鱼道内鱼类行为模拟

2021-09-04汪靖阳郄志红孙双科

农业工程学报 2021年11期

汪靖阳，郄志红，孙双科，张超

（1.河北农业大学城乡建设学院，保定 071001；2.中国水利水电科学研究院流域水循环模拟与调控国家重点试验室，北京 100038；3.国家北方山区农业工程技术研究中心，保定 071001）

0 引言1

鱼道是改善河道上下游连通性、保护物种多样性的有效措施之一。鱼道的过鱼效果不但受水流因素的影响，同时又与鱼的行为特性直接相关，二者又相互耦合，十分复杂。目前的研究中：第一类侧重鱼的行为特性和规律，以指导鱼道设计。郑铁刚等[1]根据过鱼对象的游泳能力和生活习性，对鱼道进口位置选择进行了阐述；Chen等[2]对目标鱼的游泳能力指标进行测试，并根据流速划分鱼类行为区；Silva等[3]对鱼类行为学的主要研究成果进行了分析和综合，提出了鱼道设计的建议。第二类通过过鱼试验优化鱼道的体型结构，工作量大、成本较高，且由于不同鱼类的习性不尽相同，而鱼道试验一般只针对典型鱼类，故难以克服“众口难调”的矛盾[4]。第三类是通过计算机模拟鱼道内鱼的游泳行为，如：Goodwin等[5-6]建立鱼类模型与流体总应变的响应关系，认为过鱼对象的上溯轨迹遵循Weber定律；Zhu等[7]通过数值模拟技术，采用基于欧拉-拉格朗日的智能体方法对目标鱼上溯轨迹进行模拟，其结果与Rodriguez等[8]的试验鱼轨迹相近；Tan等[9]运用ELSR模型预测异侧竖缝式鱼道内鱼类的洄游轨迹，且与试验轨迹对比发现两者吻合度较高，初步验证了鱼类轨迹预测模型的适用性；Padgett等[10]分别基于个体模型和元胞自动机模型模拟了幼年鳗鱼在鱼道的上溯过程，发现运用基于个体模型来模拟幼年鳗鱼上溯效果更优；陈求稳等[11]运用基于个体模型模拟了鱼类生长和分布的变化。国内外一般采用智能体模型（即基于个体模型）来模拟鱼类行为，而智能体模型存在计算量大，实现困难等缺点，尤其是在鱼类运动行为规则复杂时，仿真效率较低[12]。因此，本文考虑元胞自动机模型在研究复杂系统方面的优越性，拟将其融入鱼类行为模型中演化鱼类运动系统，提高鱼类行为模型模拟复杂规则时的仿真效率。

鱼道试验过鱼对象的选择需考虑易获性与代表性，相关资料对中国 2000年以后的鱼道过鱼对象进行了统计，发现裂腹鱼等冷水性鱼类和四大家鱼等鱼类占主要部分[13]，其中四大家鱼的草鱼在市面上容易获取，且国内以草鱼为过鱼对象的鱼道研究资料较多[14-16]，有助于建立鱼类行为模型。因此，本文选择草鱼幼鱼为主要研究鱼类，研究草鱼幼鱼在竖缝式鱼道中上溯行为特征，并结合智能体和元胞自动机模型的特点建立鱼类行为模型，预测草鱼幼鱼在上溯时的运动情况，为鱼道中的鱼类行为模拟研究奠定基础。

1 基本理论

1.1 元胞自动机

如图1所示，元胞自动机（Cellular Automata，CA）系统是一个由元胞空间、元胞状态、邻居和变化函数规则等元素组成的离散动力学系统[17-20]。元胞是元胞自动机模型将空间按一定规则划分后得到的最小单位，所有元胞相互离散，构成一个元胞空间；在任一时刻每个元胞都只能有一种元胞状态；邻居是元胞周围按一定形状划定的有限集合，它们影响中心元胞下一个时刻的状态；变化函数规则定义了元胞状态转换的规则。用集合的语言来描述元胞自动机模型，即：

CA=(Ω,N,S,F) （1）

式中CA是元胞自动机；Ω是离散的元胞空间；S代表元胞的状态空间，即任一时刻元胞处在的状态集合；N表示元胞的邻域；F表示元胞状态演变的局部规则，即按照前一时刻元胞的相邻状态组合来确定当前时刻状态值的函数。

1.2 基于Agent的元胞自动机建模理论

智能体（Agent）是指能够感知环境并反作用于环境的客观的或虚拟的个体。其功能主要包括：1）能感应周围环境和其他Agent的信息；2）各个Agent都拥有相互独立的规则；3）Agent能够对感知的信息作出反馈。

2 研究方法

2.1 数值模拟计算

2.1.1 数学模型

按照鱼道设计相关规范[21]，设计竖缝式鱼道模型。模型尺寸依据重力相似准则确定，具体尺寸如下：鱼道长6.675 m，宽0.6 m，坡比1:100，由7个尺寸相同的池室组成，单个池室宽为0.6 m，长度为0.675 m，竖缝宽度为0.075 m，该研究拟通过改变导板相对长度G/B以实现不同水流结构之间的转变（G为导板长度，B为鱼道池室宽度），进而对比各种水流结构条件下Agent-CA模型的模拟效果，使模型更具有可靠性。该研究选取导板相对长度为G/B＝0.1、0.25及 0.5的 3种计算工况，图2为3种工况竖缝式鱼道内部结构图。

采用Gambit软件划分网格，水池内采用网格尺寸为7×10-3m，竖缝内网格尺寸为 5×10-3m，垂向网格尺寸2.5×10-2m。采用计算流体力学软件Fluent模拟计算流场，进流段边界条件设置为压力进口，出流段边界条件设为压力出口，顶部边界设置为压力进口。

2.1.2 控制方程

控制方程采用RNGk-ε紊流方程，包括连续方程、动量方程、湍动能方程与湍动能耗散率方程[22-23]。

连续方程：

动量方程：

湍动能k方程：

湍动能耗散率ε方程：

式中下标i和j取值1、2、3；u、v、w表示各方向的流速分量，m/s；P表示作用在水体的压力，Pa；ρ=1 000 kg/m3为水流密度；k为湍动能，m2/s2，ε表示湍动耗散率，(kg·m2)/s2；μ和μt为紊动黏性系数和涡粘系数，N·s/m2；Cμ=0.0845；σk、σε为湍动能和耗散率的 Prandtl数，分别取 1.3和 1.0；η是湍流时间尺度和平均流时间尺度的比值；Cε1取1.44、Cε2取1.92。Gk是k的产生项，其控制方程如下：

2.2 物理模型设计与试验

2.2.1 物理模型的设计

物理模型依据鱼道设计导则和水力设计基本要点制作，根据重力相似准则确定模型具体尺寸。试验模型主要由鱼道试验段、摄像头、回水池等组成。鱼道试验段长度为10.2 m，宽度为1.3 m，坡比为1∶100。试验模型采用自循环系统，以水泵为动力，摄像头用于记录鱼道中间池室的过鱼情况，试验布置如图3所示。试验水槽内部布置竖缝式鱼道，鱼道体型与数学模型保持一致。3种鱼道体型其余结构尺寸参数均相同。

2.2.2 试验方案

试验装置制作完成后进行对比性过鱼试验：首先，考虑到鱼道尺寸等因素，选用 85条体长范围为（10±2）cm的草鱼作为试验用鱼，并在试验之前对草鱼均进行空腹处理，以确保试验时过鱼对象的活力；其次，试验选择第五级池室作为观察区域，将摄像头架设至中间 5#水池顶部中心处，以实时记录过鱼对象在鱼道水池内的上溯行为表现；最后选择试验鱼放置鱼道下游放鱼池中，等待15 min后打开拦鱼格栅，并同步启动摄像系统以记录鱼类的上溯信息。

2.3 基于智能体和元胞自动机的鱼道内鱼类行为模型

模型将所研究的鱼道区域划分成m×n个相同的正方形网格，每个网格点代表一个元胞，将所有连接在一起的网格点组成一个二维空间，即鱼道的二维元胞空间。每一个元胞都有4种属性特征，即障碍（如鱼道壁面）、鱼个体、鱼道区域和上游区域。对于任意一个处在(Xi,Yj)位置的元胞，以Ui,j、ui,j及vi,j分别表示该位置的流速大小、x轴方向的流速和y轴方向的流速，m/s。其状态集合S的表现形式为：

基于 Agent理论将元胞空间中每个状态为鱼个体的网格看作是Agent，这些网格在几何意义上也代表模拟鱼个体的质心。

2.3.1 鱼类感应范围

鱼类洄游是依靠侧线系统感知周围的水动力特征，从而选择上溯方向[24-25]。本试验为模拟鱼类的侧线系统感应水流的能力，对每个Agent赋予一个感应范围（如图4所示）：以元胞空间中Agent所在网格为中心，采用边长为2R的正方形平面来表示Agent的感知区域（即图4中阴影部分），其中感知平面的邻居元胞半径为R，由于关于鱼类侧线系统的感应范围尚无定论，该研究结合实际观察与相关文献中智能系统对鱼类感应范围的研究计算[5-7]，假设鱼的感应半径的取值在一个体长Lf附近波动，公式如下：

式中Lf表示鱼的体长，cm；RN是介于0～1之间的随机数。

2.3.2 游动方向和步长

传统Agent-CA模型常采用如图5a所示8个方向来代表Agent运动方向，应用上述运动方式在人员疏散模拟中取得了较好的效果[18]，但应用到鱼类行为模拟时，由于鱼个体的游动方向更加灵活，仅仅用 8个运动方向并不能准确地描述鱼类行为，因此本试验采用鱼个体运动方向和步长相结合的方式来描述鱼道内的鱼类行为：如图5b是以鱼个体的质心所占据的网格为中心，采用边长为 2Δs正方形平面的边界来表示Agent的运动区域（即图6中阴影部分），采用该运动区域来描述鱼类运动能够在增加运动方向的同时，将鱼个体能够运动的方向数与步长联系起来，提高鱼类行为模拟的准确性。在本试验的模型中鱼个体能够运动的方向数w与步长Δs可用如下公式表示：

式中a代表网格的边长。相关文献表明，鱼在鱼道内上溯时除了进入局部高流速区域（如竖缝式鱼道的竖缝处），其游泳速度变幅较小，且鱼类的克流速度主要与体长有关[26-27]。因此，为了满足模拟的准确性以及方便程序运行，做出以下假设：1）以鱼的相对临界游泳速度U′crit作为该模型中鱼类正常游泳速度vf；2）以 Δt=0.3s代表鱼类每步所用时间。参考鱼类克流能力的相关文献[14-16]设置鱼类正常游泳速度vf=4Lf（cm/s）。鱼类的上溯步长（Δs）用公式（10）计算。

2.3.3 模型假设

过鱼试验中发现，上溯鱼会根据流速、流速梯度、紊动能等因素采用对应的行为模式。在数值模拟计算中发现试验鱼道的最大紊动能TKEmax=0.025 m2/s2，属于低紊流区[28]，因此在本试验暂不考虑紊动能对鱼类运动的影响。参考鱼类行为学的相关研究成果[16,28-34]，按流速与鱼类游泳行为的特性指标间关系将鱼类上溯分为 3种场景，如表1所示。

表1 不同场景的鱼类行为模式Table 1 Fish behavior patterns in different scenes

在模型中当Agent处于以上假设的3种场景中时，会模拟鱼进行对应的运动行为，这也是该模型控制 Agent运动的主要方式，本研究结合鱼类行为资料将模型假设的3种场景对应的鱼类运动行为归结如下：

1）避开障碍行为：在洄游过程中鱼类可以通过测线感知水的流动从而规避障碍物。在模型中表示为当状态为障碍的网格出现在Agent的感知区域内时，Agent将不会向靠近障碍物的方向运动。

2）随机行为：因为该模型未考虑鱼类行为的所有影响因素，以及鱼类运动具有一定的随机性，所以在鱼类活动全过程中，都有可能向任意方向游动。在模型中表示为Agent拥有向各个方向运动的概率。

3）寻找主流行为：鱼类能通过水流来感知主流方向，并向主流游动。在模型中表示为Agent具有一定概率向感知区域内流速最大的网格移动。

4）逆流向前行为：鱼类选择水流的来流方向上溯。在模型中表示为 Agent具有一定概率沿感知区域的平均流速的反方向运动。

5）逆流后退行为：鱼类的游泳速度小于流速，鱼类会沿水流方向后退。在模型中表示为Agent具有一定概率沿感知区域的平均流速方向运动。

6）冲刺行为：鱼从高流速区域迅速进入低流速区域时，会随着惯性向前冲刺一段距离。在模型中表示为当Agent上一时间步的感知区域平均流速与当前感知区域平均流速之差大于冲刺指数时，该Agent在下一个时间步的上溯步长提高。

2.3.4 鱼个体移动规则

文献[5,7]的移动规则为 Agent下一个时间步移动到对其吸引力最大的邻居网格，但在模拟的过程中发现在固定位置放鱼时，过鱼对象均会沿相同的轨迹上溯，这是不符合实际情况的。因此本文参考陈求稳等[11]考虑鱼类行为随机性的思想，采用随机过程的方法建立移动规则：根据建立的鱼对水环境因子的响应关系，计算Agent（鱼）下一个时间步移动到运动区域（如图6）的概率。根据计算的概率，采用随机过程的方法，确定Agent（鱼）在下一个时间步最终所处的网格。

Agent向第K个方向移动的概率由公式（11）计算。

式中prob_dir(K)表示向第K个方向移动的概率；n表示网格属性系数，当网格属性为障碍时，n=0，反之，n=1；(K)为t场景下第K个方向的网格吸引力概率，其中，t=1时，认为Agent处于感应场景；t=2时，Agent处于上溯场景；t=3时，Agent处于冲回场景。

该研究的网格吸引力概率的计算公式如下：

式中At表示t场景下各方向网格吸引力权重因子之和；(K)表示t场景下I行为在第K个方向的网格吸引力权重因子，I的取值范围为1～5，分别代表避开障碍行为、随机行为、寻找主流行为、逆流向前行为和逆流后退行为。当t=1时，(K)、(K)均为0；当t=2时，(K)=0；当t=3时，(K)、(K)和(K)均为0。

Agent各个方向的网格吸引力概率主要由不同场景下的运动行为规则决定，包括避开障碍行为、随机行为、寻找主流行为、逆流向前行为和逆流后退行为，相关公式如下：

1）避开障碍行为对各方向网格吸引力权重因子的影响

避开障碍行为在模型中表示为当状态为障碍的网格出现在Agent的感知区域内时，Agent将不会向靠近障碍物的方向运动，其表达式如下：

式中Kbar表示靠近障碍物的方向；(Kbar)表示t场景下避开障碍行为使 Agent向第Kbar方向运动的网格吸引力权重因子。

2）随机行为对各方向网格吸引力权重因子的影响

随机行为在模型中表示为 Agent拥有向各个方向运动的概率，其表达式如下：

式中w表示鱼类运动方向的个数；αran表示随机行为的网格吸引力权重值；表示t场景下随机行为使Agent向第1～w方向运动的网格吸引力权重因子。

3）寻找主流行为对各方向网格吸引力权重因子的影响

寻找主流行为在模型中表示为 Agent具有一定概率向其感知区域内流速最大的网格移动，其表达式如下：

式中Kmain表示靠近主流运动的方向；αmain表示寻找主流行为的网格吸引力权重值；(Kmain)表示t场景下寻找主流行为使Agent向第Kmain方向运动的网格吸引力权重因子。

4）逆流向前行为对各方向网格吸引力权重因子的影响

逆流向前行为在模型中表示为 Agent具有一定概率沿感知区域的平均流速的反方向运动，其表达式如下：

式中Kre表示逆流向前运动的方向；αre表示逆流向前行为的网格吸引力权重值；(Kre)表示t场景下逆流向前行为使 Agent向第Kre方向运动的网格吸引力权重因子。

5）逆流后退行为对各方向网格吸引力权重因子的影响

逆流后退行为在模型中表示为 Agent具有一定概率沿感知区域的平均流速方向运动，其表达式如下：

式中Kback表示逆流后退运动的方向；αback表示逆流后退行为的网格吸引力权重值；(Kback)表示t场景下逆流后退行为使Agent向第Kback方向运动的网格吸引力权重因子。

6）冲刺行为对鱼个体移动规则的影响

冲刺行为在模型中表示为当 Agent上一时间步的感知区域平均流速与当前感知区域平均流速之差大于冲刺指数Tgra时，该Agent在下一个时间步的上溯步长提高，其表达式如下：

式中ΔSsprint表示冲刺行为移动的步长；Δv表示Agent两个时间步的流速差值；K′表示冲刺系数，经检验发现，K′取0.833时，更符合鱼类实际运动情况。

3 实例模拟与对比分析

3.1 流速分布规律

相关研究[35]表明竖缝式鱼道水流结构具有典型的二元特性，水力因子沿垂向变化较小，并且草鱼在上溯时较倾向于中底层活动，因此该研究选取水深为 0.5H（H为鱼道水深）平面的流场作为研究区域。

图6给出了导板相对长度为G/B＝0.1、0.25及0.5下3种工况下竖缝式鱼道中间池室0.5H处流速分布情况。由图6可以看出，G/B＝0.1时鱼道主流蜿蜒性较弱，主流左侧回流区域大于右侧；G/B＝0.25时鱼道主流蜿蜒性相对增加，使主流集中于池室中心，两侧回流区域尺度相近；G/B＝0.5时鱼道主流蜿蜒性进一步提升，使左侧回流区进一步被压缩，尺度明显小于右侧。

3.2 过鱼试验对比分析

为了便于观察分析鱼类行为与流场水力因子间的响应关系，将3种工况下草鱼的上溯轨迹线与图7的流速场进行叠加，并将轨迹分布宽度W标注在图中，如图7a、图7b和图7c所示。由于本试验选用的试验用鱼数目较多，如果将所有轨迹绘制在一个池室内会造成轨迹十分散乱，因此该研究基于鱼类所有轨迹点坐标作为基础数据，对原始上溯轨迹进行特征化处理[33]，如图9d、图9e和图9f:1）先将鱼道区域划分为a×b的网格，网格内各点位置用坐标(x,y)表示；2）针对相同的xi值，记录所有轨迹M(xi,yj)横向位置yj，并对所有的yi值从小到大进行排序；3）针对相同的xi值，查找排序为25%、50%、75%的位置点处所对应的yj值；4）将各xi值对应的3个频率的yj点绘出，简记为轨迹T25、T50、T75。

由图7可知，G/B=0.1的鱼道中试验鱼主要沿主流上溯，上溯路程较短；而在G/B=0.25的鱼道中大部分以冲刺的方式通过竖缝处进入池室右侧，而后在上游竖缝附近寻找主流，通过竖缝游至上一级水池；在G/B=0.5的鱼道中，草鱼的上溯轨迹与G/B=0.25的鱼道中相似，先冲游通过竖缝处，再在池室中央寻找主流进入上一级池室。

3.3 模拟与分析

进行鱼类行为模拟前，先确定模型的重要参数值：随机行为主要反映了鱼类活动的随机性，并不是鱼类上溯过程的主要行为，因此该行为的网格吸引力权重值较小，该研究设置aran=5；由图7的特征轨迹线可以看出，寻找主流行为和逆流向前行为是指导鱼类上溯的主要行为，发生在整个上溯过程的40%～60%，因此定义其网格吸引力权重值分别为amain=45和are=50；逆流后退行为主要发生在冲回的场景，此时的流速对鱼的胁迫作用非常明显，其网格吸引力权重值aback=50；由图9b和图9c观察可知，过鱼对象的冲刺行为主要发生在竖缝附近，当单位时间步内流速差在0.07～0.1 m/s时，极易发生冲刺行为，因此该研究设定冲刺指数Tgra=0.08 m/s。

根据上述数值模拟获得的流场数据和竖缝式鱼道几何特征建立基于Agent-CA的鱼道内鱼类行为模型：先将竖缝式鱼道划分为一个834×75二维网格，大小为8 mm×8 mm；再将Fluent软件计算所得的水深为0.5H的流场数据（流速、x轴方向流速分量和y轴方向流速分量）导入元胞空间中；最后在第7级池室中生成1条体长在（10±2）cm范围Agent（鱼）进行鱼类行为模拟，重复模拟85次，即可得到85个Agent（鱼）的原始轨迹分布图。

模型模拟得到G/B=0.1、0.25及0.5鱼道Agent的原始轨迹分布图，如图8a、图8b和图8c所示。

在上述模拟的轨迹图中，黑色部分代表竖缝式鱼道内的导板或者隔板，池室内各种颜色的散点代表模型随机生成的各个Agent的运动轨迹，鉴于模拟Agent数目较多，为了便于观察与分析，按3.2节方法将G/B=0.1、0.25及0.5的Agent轨迹线进行特征化处理，得到如图8d、图8e和图8f的特征轨迹分布图。

为了反映模型的模拟效果，本文选用轨迹范围（W/B）的相对误差作为评价模型的指标（W为轨迹宽度，B为鱼道池室宽度）。表2给出了模拟轨迹和模拟特征轨迹的相对误差，结合图7和图8中试验轨迹与模拟轨迹分析：通过比较图7和图8的轨迹形态，可以看出模拟特征轨迹与试验轨迹较为接近，且模拟特征轨迹范围相对误差较低，最大值为23.5%，说明本文建立的Agent-CA模型能够较为准确地描述鱼类主要运动范围。但对比模拟轨迹和试验轨迹发现模拟轨迹分布更为集中，与试验轨迹分布有所差异，且模拟轨迹范围相对误差较高，最大值为42.5%，造成这一现象的原因可能是由于模型中鱼类运动随机性较弱，使Agent（鱼）不易偏离主流上溯，因此提高模型中鱼类运动的随机性仍是今后要考虑的问题。

表2 3种导板相对长度下试验轨迹范围、模拟轨迹范围和相对误差Table 2 Experimental track scope, simulated track scope and relative error of three types of relative length of guide plates