CN104299265B - 一种流体环境影响下的群体行为控制方法 - Google Patents

Abstract

本发明提供了一种在流体环境影响下的群体行为控制方法。在基于物理的流体仿真方面，通过求解计算流体动力学方程（Navier‑Stokes方程）中的平流项、外力项和压力项，更新流体的水平集、速度等物理量，基于GPU（Graphic Processing Unit）完成流体的仿真计算。在基于势能场的群体行为规划方面，对于局部势能场和全局势能场进行有效的融合。从而在保证路径最优的同时，提高运算效率。在流体环境对于群体行为的影响方面，主要是通过对于个体（Agent）产生作用力来影响最终的导航力，以达到对于群体行为控制的目的。本发明能够有效地实现流体环境影响下的群体行为控制的真实效果。

Description

一种流体环境影响下的群体行为控制方法

技术领域

本发明属于虚拟现实技术领域，尤其涉及一种流体环境影响下的群体行为控制方法。

背景技术

群体动画在计算机游戏、影视动漫、城市建筑规划等很多方面都有着广泛的应用。然而，对大规模群体行为进行仿真尤其是实时仿真是一件复杂而艰巨的工作。因为对群体行为的仿真不仅涉及到高层的决策过程，还同时要考虑底层事物的表示以及事物之间的交互计算。更为麻烦的是，个体与个体之间以及个体与环境之间的复杂约束关系决定了在群体规模增大的时候，仿真运算复杂度将呈非线性增长。

在计算机仿真研究领域，流体仿真是其中一个重要的研究方向，它就像是一种有用的建筑积木，是仿真很多自然现象的基础。而基于物理学的流体(如烟、水、火)动画提供了计算机图形中最让人叹为观止的视觉效果。然而由于其巨大的计算成本，它在历史上一直属于高质量离线渲染的领域。在实时应用中，3D流体的效果往往都是依赖于粒子系统实现。尽管粒子系统能够产生很具有吸引力的结果，但是它们并不符合流体的真实外观和行为。基于物理的实时流体仿真仍然是一种挑战，尤其是流体与各种各样的障碍物进行交互的仿真模拟。

流体环境影响下的群体行为控制很少有人研究，这是因为其本身的复杂度和计算机软硬件条件的限制。但是，群体行为控制中添加流体的影响可以增加模拟的真实度。同时，在游戏、动画等实际应用方面也有较大的需求。

发明内容

为了使群体行为更接近真实效果，本发明提供了一种流体环境影响下的群体行为控制方法。本发明分别实现了基于物理的流体仿真、基于势能场(势能使用一个数值进行表示，群体运动的趋势是从高势能的地方到低势能的地方，对于场景中的个体来说，会判断周围势能的大小，从而向势能低的地方运动)的群体行为规划和流体环境对于群体行为的影响。在基于物理的流体仿真方面，通过求解计算流体动力学方程中的平流项、外力项和压力项，更新流体的水平集、速度等物理量，完成流体的仿真计算。对于流体中求出的物理量，将其存储在纹理(一种数据结构，类似于CPU上的数组。在GPU上对于数据进行读写操作，大部分都是基于纹理这种数据结构。对于纹理而言，是由纹理元素简称纹素构成，它可以设定为单个通道或者多个通道)中。然后，流体的仿真计算就可以在GPU(图形处理器，它是相对于CPU的一个概念，是一个专门的图形的核心处理器)下进行，大大提高了仿真的效率。在基于势能场的群体行为规划方面，将目的地设置为势能最低的地方，对于个体周围和障碍物周围则设置较高的势能。从而能够使得群体向着目的地运动，对于周围的障碍物和其他个体能够进行很好的规避。对于局部势能场和全局势能场进行有效的融合，从而在保证路径最优的同时，提高运算效率。将整个群体划分为多个小群体，这些小群体具有共同的特征。通过这种方式，可以减少势能的计算，提高效率。在流体环境对于群体行为的影响方面，主要是通过流体环境对个体产生作用力，从而影响个体的速度，然后达到对于群体行为控制的目的。本发明能够有效的实现流体环境影响下的群体行为控制的真实效果。

本发明一种流体环境影响下的群体行为控制方法，该方法包括步骤：

步骤1，初始化流体的水平集、速度等物理量，对流体的水平集和速度进行平流；

步骤2，处理静态场景物体，生成静态局部势能网格；

步骤3，根据“静态局部势能”网格，针对每个网格点，生成对应的全局势能网格；

步骤4，处理动态物体，更新周围网格；

步骤5，对于路径规划所需要考虑的因素进行计算，这里的因素主要是指目的地的位置和个体所受的作用力；路径规划是指个体为了到达目的地所选择的路径；

对于路径规划所需要考虑的因素进行计算，这里的因素主要是指目的地的位置和个体所受到的作用力；目的地的位置使用全局的最小势能表示；如果全局中有若干个目的地，则令这些目的地的势能具有同样的大小，即同时保持全局最小；当群体中个体的目的地信息确定之后，接下来最重要的因素就是群体中个体所受到的作用力；通过个体所受到的作用力，可以得到加速度，从而对于个体的速度产生影响；假设个体是以环境所能允许的最大速度行进，这里的环境是指上下坡的坡度以及人流量的大小因素；计算个体速度如公式(1)所示：

v＝u(x,θ)n_θ (1)

其中，u表示最大速度场，x为个体所在位置，θ表示速度的方向，n_θ＝[cosθ,sinθ]^T表示在θ方向上的单位向量。

上述的个体速度v没有考虑周围环境的影响，这里使用g表示周围环境的影响，即是个体受到周围环境的作用力，对于g来说，它是可以动态变化的。通过动态变化的环境作用力可以帮助个体规避动态的障碍物，例如汽车等。

将上述因素结合起来考虑，下面描述个体进行路径规划需要考虑的因素。总的来说，个体会选择到达目的地的最短路径。假定D为从地点x到目的地的路径集合，一个位于x的个体选择的路径P∈D，按照如下公式(2)所示：

L+T+g (2)

其中，L表示路径的长度，T表示所消耗的时间，g表示环境作用力。这里的公式(2)描述了路径规划所需要考虑的三个因素，即路径长度、通过路径的时间和所受的环境作用力。通过这三个因素的权衡，可以得到最终所选择的路径。

步骤6，处理流体环境，计算流体环境对于周围个体产生的作用力，并计

算最终的导航力；所谓导航力，是通过对于个体的受力进行分析之后，得到的最终的合力，它会影响个体如何运动。

处理流体环境，计算流体环境对于周围个体产生的作用力，并计算个体最终的导航力时，又进一步包括：

步骤6.1，对流体区域所在的网格进行切片，在每个切片上区分是个体还是流体。对于每一个切片来说，个体区域和流体区域分别使用不同数值进行标记，这样就能区分每一个切片上的个体区域和流体区域。然后，做一次最终遍历，根据每一个切片，把数值复制到3D纹理(这个3D纹理称之为流固3D纹理)中去。

步骤6.2，与步骤6.1的做法类似，对于流体网格进行切片。然后，将网格切片与个体网状模型进行相交操作。在将切片与网状模型相交时，通过将个体网状模型划分成三角形图元来实现。在得到三角形图元与切片的交点之后，通过插值即可得到流体与个体相交的边界点。然后，创建一个新的3D纹理(这个3D纹理称之为边界3D纹理，同流固3D纹理的大小和维度一致)用于存放边界点的信息。在得到边界3D纹理之后，通过对于边界点周围的单元进行判断，可以得到边界点处的法向。将边界点的法向存储于一个新的3D纹理(这个3D纹理称之为法向3D纹理，同边界3D纹理的大小和维度一致)。

步骤6.3，考虑流体环境对于个体施加的作用力，通过将个体所受的全部作用力进行矢量合成之后，即可得到最终的导航力。流体环境对于个体施加的作用力如公式(3)所示：

其中，Div表示流体网格上的离散化的散度求解符号， 对于A和b来说，c是一个3×m的雅克比矩阵(m代表个体模型自由度(degreesoffreedom)的数量)，令[n₁...n_k]为边界点(流体与个体相交点)的法线向量，则M为个体模型的质量矩阵。S为一个k×n的选择矩阵，用来选择耦合处的压力p，k表示流体与个体的相交点的数量，n表示流体单元的数量。向量S_p则是包含了个体周围的所有压力值。S_p中的每一个分量p_i都对个体施加了一个为(Δx)²p_in_i的压力。其中的是矩阵c对于时间的导数，s^*为个体的中间速度。u^*为流体的中间速度场。

至此，根据上述公式(3)，就可以求出流体环境对于个体产生的压力p，也即流体环境对于个体施加的作用力。通过对于个体所受作用力进行合成之后，即可得到最终的个体导航力。

步骤7，根据最终的导航力产生的加速度，可以得到个体的速度，通过速度来实现个体的路径规划以及规避障碍物；

根据最终的导航力所产生的加速度，可以得到个体的速度，通过速度来实现个体的路径规划和规避障碍物。在公式(2)描述的路径规划所需的开销基础之上，下面给出最优路径的选择方法。假定有一个函数φ，这个函数描述了一条路径。该路径上任意点的开销与最优路径都是一样的。直觉上来说，对于任何个体，最优的策略就是沿着这个函数梯度的相反方向移动，这是降低路径开销最为迅速的方式。事实上，函数φ可以按照下述方式精确定义：在目标点φ＝0，其他的地方φ满足如下方程(4)：

||▽φ(x)||＝C (4)

其中，单位开销C是由梯度▽φ的方向来得到。

令a(x)为位于位置x的个体受到的导航力所产生的加速度，结合公式(1)可以得到速度的计算如公式(5)所示：

其中，v表示速度，u(x,θ)由个体的位置和运动方向而得的最大速度场。由于u与群体密度ρ密切相关。因此，在求解u的时候，根据群体密度的不同，将密度划分为高密度、低密度和中等密度这三种情况来分别进行求解。

步骤8，求解流体环境的压力项，得到流体的最终速度，并调整边界点处的水平集信息；

求解流体环境的压力项，得到流体的最终速度，并调整边界点处的水平集信息。在步骤1中对于流体速度的平流所得到的是中间的散度速度，这个中间速度需要减去压力项的梯度才能得到最终的无散度速度。因此，需要求解压力项。对于压力项的求解采用了Jacobi迭代的方法。这里需要注意的一点是，在步骤1中并没有考虑流体环境与个体的相互作用。由于流体环境对于个体会施加一个作用力，同时个体对于流体环境也会施加一个反作用力(大小相同，方向相反)。这个反作用力会影响边界点处流体的压力项，因此边界点处流体的速度也会发生改变。同时，需要根据边界点处的流体速度对于水平集进行修正，才能得到最终的水平集信息。

步骤9，对于流体环境与个体的边界进行处理，使两者在边界处满足自由滑动边界条件。

本发明的有益效果是：在基于物理的流体仿真方面，将流体的物理信息保存在纹理当中，使得仿真在GPU上得以实现，提高了仿真效率。在更新水平集时，采取了平流与调整相结合的方法，获得了更高的效率。在基于势能场的群体行为规划方面，将局部势能场和全局势能场进行有效的融合。由于局部势能场方式存在局部极小值问题，不能保证到达最终目标点。而针对特定目标点的全局势能场计算又太费时。因此，将二者进行融合：全局势能建立在局部势能之上。同时，全局势能网格仅取低精度，在保证大致全局导航最优的同时，将局部的避让细节交由局部势能处理，从而产生平滑的导航路径。在流体环境对于群体行为的影响方面，需要进行三步操作：区分流体网格中的个体和流体部分、得到个体的边界点以及个体边界点处的法向、计算流体环境对于个体产生的作用力。然后，根据作用力得到个体的导航力。导航力产生的加速度对于个体的速度产生影响，从而影响到个体的路径规划以及规避障碍物。

附图说明

图1示出了本发明一种流体环境影响下的群体行为控制方法流程图；

图2示出了逆向追踪的平流方法；

图3示出了规则几何体的势能作用力示意图；

图4示出了双线性插值；

图5示出了梯度计算示意图；

图6示出了局部势能分布图、路径拓扑图和求解到的全局势能(从左到右)；

图7示出了流体网格切片的示意图；

具体实施方案

下面结合附图和实施实例对本发明优先实施方式进一步说明：

图1所示的流程图给出了本发明整个实施的具体过程：

步骤1，初始化流体的水平集、速度等物理量，对流体的水平集和速度进行平流。在初始化流体的物理量之前，需要初始化流体的仿真网格。然后根据网格创建3D纹理，保存流体的初始状态(水平集、速度等信息)。在流体仿真的过程中，速度是最重要的量，因为速度决定了流体如何移动其本身，以及在流体中的东西。由速度传送流体本身和在流体中的其他量，这个过程就是平流。本发明中计算平流是通过逆向处理和使用隐式积分法。不是通过计算粒子在当前时间段移动到哪里来得到移动量，而是从每个网格单元逆时追踪粒子轨道，得到它以前的位置，并到开始的网格单元，复制那个位置的量。如图2所示，为了计算以双圈标记的单元的平流，逆时追踪速度场到X。对最靠近X的4个网格值(在图中连接成一个正方形)进行双线性插值，而对出发的网格单元写下结果。将这种平流方法应用到速度和水平集上，得到中间的速度和水平集信息。

步骤2，处理静态场景物体，生成“静态局部势能”网格。具体包括的步骤如下：

步骤2.1，为了避免额外的梯度计算，本发明直接采用二维向量数组进行存储，数组值代表地理网格顶点上的势能值。向量方向代表势能梯度下降方向，大小代表势能对周围物体的作用力大小。如图3所示。

为了求得二维空间任意点上的平滑势能作用力，这里采用双线性插值进行计算。如图4所示。分别代表网格顶点p上的势能；L为网格大小；X和Y分别代表距格子边的距离。则势能如公式(6)所示。

步骤2.2，对于规则物体，采用其边沿及其平行扩展作为势能等高线，因此势能方向垂直几何边缘。势能大小采用类似库伦力的方式计算。为了方便描述，引入两个参数---影响范围δ和影响力度λ。影响范围δ表示几何体对其边界以外周围δ单位以内的点存在影响：影响力度λ表示几何体对影响范围以内的点的影响程度比例参数。则几何体外部且距离边缘的最短距离向量为上的点所受的势能如公式(7)计算求得。

当点处于几何体边沿时，势能大小为无穷大；当处于作用范围边缘时，势能大小为0。为了简化考虑，所有处于几何体内部的点所受力大小为无穷大；所有处于物体作用范围δ外的点所受力大小为0。

当一个点受到多个几何体的势能作用时，采用矢量叠加方式处理，从而可以方便地进行几何体的添加和删除操作。图3分别显示了方形、圆形、椭圆形的势能作用方式。

步骤2.3，本发明使用采样方式近似计算地形梯度。如图5所示，v₁…v₈为网格顶点P周围的八个顶点。将该九个顶点垂直映射到地形上，求得对应8个三维空间向量中与水平面夹角最大的向量V_max。取V_max方向作为梯度V的方向：V_max与水平面夹角的正切值作为梯度V的大小。

地形产生的势能将叠加到已求得的势能当中。假设在地形某点处的地形曲面梯度为V，行走难度为λ。则该点的势能使用公式(8)求得。

F＝-λ*V (8)

步骤3，根据“静态局部势能”网格，针对每个“网格点”，生成对应的“全局势能”网格。如图6所示，左图是已经构建好的局部势能。如前所述，局部势能是采用作用力的形式存储起来的，这里只关心其大小而不必考虑方向。现要求解对目标点A的全局导航势能。

这里将全局网格精度设为局部网格的1/3，构建一张路径拓扑图，如图6中间所示，图中每条线段的权值采用公式(9)计算。

其中F(i)表示i点的局部势能，意思为从顶点A开始，以局部格子大小为采样步长，直到顶点B，所有采样点上的局部势能作用力大小之和。Len(AB)表示顶点A到顶点B的距离。分别代表局部势能和距离对全局势能的影响权值。

当完成对每条线段的权值计算后，就可以求出整个全局势能场。整个运算过程并不求出每条到达A点的单独路径，而是在每个网格点上存储一个指向其最近邻接点的向量作为全局路径导航作用力。如图6右所示。这里需要注意的是，此时的全局路径导航力是只考虑了静态局部势能的结果。

步骤4，处理动态物体，更新周围网格。动态物体的势能作用力计算与静态物体非常类似。动态物体在一定时间内运行的路径会形成一个几何图形(例如，运动的汽车在一定时间内的轨迹会形成一个几何图形)，利用该几何图形的静态计算模型即为该物体的动态计算模型。

处理动态物体中的行为个体(Agent)，为了准确反映行为个体间的相互避让作用关系，个体信息被放入周围网格。网格大小一般取行为个体的感知半径(在个体运动的过程中，需要设定一个数值。当物体和个体的距离小于这个数值时，个体能知道这个物体的存在，从而在运动中将这个物体的影响考虑进去。这个数值就称之为感知半径)大小。行为个体根据其当前所在位置，将其势能信息(一般取较大值，即个体表示为障碍物)放入相应格子中。每个行为个体直接从其附近的格子中获悉周围是否有其他个体，从而进行避让。

行为个体从其所在格子及周围的九个格子中获得行为个体信息。当行为个体从一个格子走入另外一个格子后，行为个体从前一个格子中注销，然后注册到另外一个格子中。

步骤5，对于路径规划需要考虑的因素进行计算。这里的因素主要是指目的地的位置和个体所受到的作用力；路径规划的首要因素就是目的地的位置。目的地的位置使用全局的最小势能表示。如果全局中有若干个目的地，则令这些目的地的势能具有同样的大小，即同时保持全局最小。当群体中个体的目的地信息确定之后(大量个体可能具有同样的目的地)，接下来最重要的因素就是群体中个体所受的作用力。通过个体所受的作用力可以得到加速度，从而影响到个体的速度。总的来说，假设个体是以环境所能允许的最大速度进行。环境可以通过以下方式来影响个体的速度：1、地形因素。例如，当遇到下坡时，个体会加快移动的速度。2、物理边界。例如，当遇到障碍物等物理边界的时候，个体是不能穿过障碍物的。3、个体之间的相互影响。本实例中假定当群体的密度较大的时候，个体的移动速度会降低。一个极端的情况是，两个个体之间不能互相碰撞。计算速度如公式(1)所示。这里的环境是指上下坡的坡度以及人流量的大小因素；计算个体速度如公式(1)所示：

v＝u(x,θ)n_θ (1)

其中，u表示最大速度场，x为个体所在位置，θ表示速度的方向，n_θ＝[cosθ,sinθ]^T表示在θ方向上的单位向量。

上述的个体速度v没有考虑周围环境的影响，这里使用g表示周围环境的影响，即是个体受到周围环境的作用力，对于g来说，它是可以动态变化的。通过动态变化的环境作用力可以帮助个体规避动态的障碍物，例如汽车等。

将上述因素结合起来考虑，下面描述个体进行路径规划需要考虑的因素。

即使个体可以自由移动，也会选择特定的路径。因为周围环境会对个体产生一定的影响，使用g表示周围环境的影响，即个体受到周围环境的作用力。在后续步骤中所讨论的流体环境，也是属于周围环境的范畴，即流体环境属于g的一部分。关于流体环境对于个体的影响，在步骤6中进行详细介绍。

综合考虑以上因素，下面描述个体如何进行路径规划。总的来说，个体会选择到达目的地的最短路径。假定D为从地点x到目的地的路径集合，一个位于x的个体选择的路径P∈D，按照如下公式(2)所示：

L+T+g (2)

其中，L表示路径的长度，T表示所消耗的时间，g表示环境作用力。这里的公式(2)描述了路径规划所需要考虑的三个因素，即路径长度、通过路径的时间和所受的环境作用力。通过这三个因素的权衡，可以得到最终所选择的路径。路径规划的开销如公式(2)所示。

步骤6，处理流体环境，计算流体环境对于周围个体产生的作用力，并计算最终的导航力，具体的步骤如下：

步骤6.1，对流体区域所在的网格进行切片，在每个切片上区分是个体还是流体。其思想很简单：使用矩阵投影把输入的三角形网状结构在目的3D纹理结构的每一个切片里渲染一次。这里实现的方式就是将流体区域所在的网格进行切片，对于每一个切片来说，个体区域和流体区域分别使用不同数值进行标记，这样就能区分每一个切片上的个体区域和流体区域。如图7所示，图中每一个切片中黑色部分即是个体区域，其余白色部分为流体区域。然后，做一次最终遍历，根据每一个切片，把数值复制到3D纹理(这个3D纹理称之为流固3D纹理)中去。

步骤6.2，与步骤6.1的做法类似，首先需要对于流体网格进行切片，然后将切片与个体网状模型进行相交操作。由于个体网状模型可以视作由三角形图元构成，因此切片与网状模型的相交即可转化为切片与三角形图元的相交操作。在具体的实施过程中，通过将三角形图元的每条边与切片进行相交来得到交点。然后，根据交点进行插值来得到流体与个体相交的边界点。创建一个新的3D纹理(这个3D纹理称之为边界3D纹理，与流固3D纹理的大小和维度保持一致)用于存放边界点的信息。在得到边界3D纹理之后，根据边界点周围的单元是否存在流体单元来得到边界点的法向信息。具体来说，如果边界点相邻单元是流体单元，则该流体单元对于边界点产生一个法向分量。当边界点周围的单元都判断完毕之后，将这些法向分量进行合成，即可得到最终的法向。创建一个新的3D纹理(这个3D纹理称之为法向3D纹理，同边界3D纹理的大小和维度一致)存储边界点的法向信息。

步骤6.3，考虑流体环境对于个体施加的作用力，通过将个体所受的全部作用力进行矢量合成之后，即可得到最终的导航力。流体环境对于个体施加的作用力如公式(3)所示：

其中，Div表示流体网格上的离散化的散度求解符号， 对于A和b来说，c是一个3×m的雅克比矩阵(m代表个体模型自由度(degreesoffreedom)的数量)，令[n₁...n_k]为边界点(流体与个体相交点)的法线向量，则M为个体模型的质量矩阵。S为一个k×n的选择矩阵，用来选择耦合处的压力p，k表示流体与个体的相交点的数量，n表示流体单元的数量。向量S_p则是包含了个体周围的所有压力值。S_p中的每一个分量p_i都对个体施加了一个为(Δx)²p_in_i的压力。其中的是矩阵c对于时间的导数，s^*为个体的中间速度。u^*为流体的中间速度场。

在得到流体环境和个体相交的边界点之后，通过公式(3)可以得到流体环境对个体施加的作用力。下面介绍得到公式(3)的过程。

通过求解下面的方程来计算作用力产生的加速度，如公式(10)所示：

其中，c是一个3×m的雅克比矩阵(m代表个体模型自由度的数量)，令[n₁...n_k]为边界点(流体与个体相交点)的法线向量，则M为个体模型的质量矩阵。S为一个k×n的选择矩阵，k表示流体与个体的边界点的数量，n表示流体单元的数量。向量Sp则是包含了个体周围的压力值。Sp中的每一个分量p_i都对个体施加了一个为(Δx)²p_in_i的压力。其中的是矩阵c对于时间的导数，s^*为个体的中间速度。

注：在统计模型中，自由度是指样本中可以自由变动的变量的个数，即能独立变化的数据数目。这里的自由度是指个体模型可以划分成相互独立部分的个数。

由于流体要满足不可压缩条件，因此流体速度在下一步骤开始前应当保持散度为0，如公式(11)所示：

▽·uⁿ⁺¹＝▽·(u^*+Δta)＝0 (11)

通过对于公式(11)的变形，可以得出公式(12)：

将上述的推导过程汇总，即可得到流体环境对个体施加的作用力如公式(3)所示。将个体所受的作用力进行合成之后，可以得到最终的个体导航力。

步骤7，根据最终的导航力所产生的加速度，可以得到个体的速度，通过速度来实现个体的路径规划以及规避障碍物。在公式(2)描述的路径规划所需的开销基础之上，下面给出最优路径的选择方法。假定有一个函数φ:R²-＞R，(R²-＞R表示将一个二维数据转化为一维数据)，这个函数描述了一条路径。该路径上任意点的开销与最优路径都是一样的。直觉上来说，对于任何个体，最优的策略就是沿着这个函数梯度的相反方向移动，这是降低路径成本最为迅速的方式。事实上，函数φ可以按照下述方式精确定义：在目标点φ＝0，其他的地方φ满足如下方程(4)：

||▽φ(x)||＝C (4)

其中，单位开销C是由梯度▽φ的方向来得到。

令a(x)为位于位置x的个体受到的导航力所产生的加速度，结合公式(1)可以得到速度的计算如公式(5)所示：

其中，v表示速度，u(x,θ)由个体的位置和运动方向而得的最大速度场。由于u与群体密度ρ密切相关。因此，在求解u的时候，根据群体密度的不同，将密度划分为高密度、低密度和中等密度这三种情况来分别进行求解。

计算势能来寻找最优路径似乎是一件比较麻烦的事情。不过，可以对于势能计算进行简化。假设一群个体有着相同的速度场、目的地以及周围环境。针对这种情况，对于群体而言只需要计算一次势能即可，群体中的成员共享这个结果。事实上，不同的个体有着不同的速度、处在不同的环境，并且目的地也不尽相同。因此，本实例将群体按照目的地的不同划分成不同的小组，每个小组的个体拥有相同的目的地。

下面介绍公式(5)中速度场u的求解。速度是群体密度相关的变量。在群体密度较低的地方，速度由地形来控制。在平坦的地形，速度保持为一个常量。但是，速度会随着斜率的改变而改变。在群体密度较高的地方，速度是由周围的群体运动来控制的。当个体运动方向是逆着群体运动的方向时，个体运动受到阻碍；当个体运动的方向是顺着群体运动的方向时，个体运动不受影响。

由于速度是群体密度相关的，引入密度场ρ来表示群体密度。然后把个体转化到群体的密度场中，令ρ_i表示第i个个体所在的密度场。这个密度场的密度在个体所在的位置达到最高峰，然后向四周进行衰减。群体的密度场就是将个体的密度场进行求和而得的。当计算群体的密度场时，可以同时计算平均速度v，平均速度就是通过速度对于密度的平均而得到的。如公式(13)所示：

和

其中，v_i表示第i个个体的速度。

下面具体描述群体密度如何影响速度。在低密度的区域(ρ≤ρ_min)，速度u等于地形速度u_T。令地形的最大和最小斜率为s_min和s_max，速度与斜率成反比。如公式(14)所示：

其中，▽h(x)·n_θ是沿着方向θ高度场h的斜率。

在高密度的区域(ρ≥ρ_max)，速度u等于群体的速度u_v，如公式(15)所示：

群体速度本质上就是从当前位置x移动距离r后，所得到的平均速度

在中等密度的区域(ρ_min≤ρ≤ρ_max)，在地形速度和群体速度直接进行线性插值。如公式(16)所示：

需要注意的是，就像群体的速度一样，密度ρ通过rn_θ的偏移来进行计算的。原因是不应该让个体对于密度场的作用反过来成为阻碍自己运动的力量。

步骤8，求解流体环境的压力项，得到流体的最终速度，并调整边界点处的水平集信息。在步骤1中对于流体的速度进行了平流，不过这个速度是中间的散度速度。想要得到最终的无散度速度，就需要减去压力项的梯度。

本实例中是通过Jacobi迭代求解泊松压力方程计算压力，泊松压力方程如公式(17)所示：

▽·▽p＝▽·w (17)

其中，w表示流体的中间速度，▽·▽p中间表示先对压力求梯度，然后求散度。

针对泊松压力方程的求解，需要将方程离散化，重写为公式(18)的形式：

其中，α和β是常数。在泊松压力方程中，x代表p，b代表▽·w，α＝-(δx)²，而β＝6。上标k和k+1表示迭代次数。

引入散度和拉普拉斯算子的有限微分形式将公式(17)离散化。如公式(19)和公式(20)所示：

假设网格单元的三维是相同的，即δx＝δy＝δz。对于公式(19)和公式(20)整理可得公式(21)：

其中，p_i＝p_i+1,j,k+p_i-1,j,k，p_j＝p_i,j+1,k+p_i,j-1,k，p_k＝p_i,j,k+1+p_i,j,k-1。

方程离散化成目标形式之后，简单地运行若干次迭代就可以求解方程。即把前一个迭代结果作为下一个迭代的输入(x^(k+1)变成x^(k))，在每个网格单元应用方程(21)。

这里需要注意的是，对于流体环境和个体相交的边界处流体，压力项通过步骤6中所求得的作用力来进行修正。即相交边界处的流体压力项由流体环境对于个体作用力的反作用力来求得。

在得到流体的压力项之后，计算其梯度值，采用Helmholtz-Hodge分解定理(任何的矢量场都可以分解为两个其他矢量场之和：一个是无散度矢量场，另一个是标量场的梯度)将流体的中间速度投影到其无散度场上去，得到流体的最终速度，保证仿真前后都满足流体的不可压缩条件。如公式(22)所示：

u＝u^*-▽p (22)

其中，u表示散度为零的速度(流体最终速度)；u^*表示散度非零的中间速度。

由于流体使用了水平集方法(一种数值方法，主要用于解决曲线演化问题。它具有计算稳定，适应任意维空间的特点，这种方法在图像处理领域中大量使用，尤其是图像分割领域)，现在给出水平集基本方程的推导：

用x(t)表示水平集，用φ表示水平集函数，根据两者的对应关系有公式(23)：

φ(x(t),t)＝0 (23)

其中，t表示时间。将公式(23)两端对于时间求偏导，则有公式(24)：

其中，▽符号为梯度算子。

假设F为边界外法向方向的速度，则有公式(25)：

其中，n＝▽φ/|▽φ|。

将公式(25)中所求得的▽φ代入公式(24)中可得公式(26)：

由公式(26)可以看出，水平集的变化只与运动速度(即公式(26)中的边界外法向方向的速度F有关)。

根据水平集的基本方程可知，可以在得到流体的最终速度之后对于水平集进行求解。但是，本实例中并没有采用这种方法。考虑到流体环境需要与个体进行交互，本实例中先对水平集进行平流(水平集的平流方法与速度的平流方法一致)。然后，在求得流体的最终速度之后，对于流体环境和个体相交点处的水平集进行调整。通过这种方法，可以获得比直接根据最终速度求解水平集更高的效率。这是因为在流体环境与个体进行交互时需要对于流体仿真网格(由于在计算机上对Navier-Stokes方程进行求解需要将方程离散化，因此需要对于流体限定一个仿真网格。这里使用的网格是MAC网格，它是一种交错网格，将流体的速度存储在网格的边界上，而将流体其他的物理量存储在网格的中心)进行切片。如果有了平流之后的水平集信息，可以排除水平集之上的切片(这部分只有个体，没有流体)，从而减少了计算量，提高效率。

步骤9，对于流体环境与个体的边界进行处理，使两者在边界处满足自由滑动边界条件。所谓自由滑动边界条件，简单来说，就是流体不能流进个体和从个体中流出，但是允许沿着个体的表面自由流动。这个条件表示流体的速度与个体的速度在边界法向方向上的分量是一样的，如公式(27)所示：

u_fluid·n＝u_solid·n (27)

在本实例中，为了满足自由滑动边界条件，使用的方法是：首先，计算出个体的速度在边界法向方向上的分量；然后，使用这些分量替换相邻流体单元中的法向方向上的分量。

通过使用这种方法，在流体环境与个体的边界处满足了自由滑动边界条件。

Claims (1)

1.一种流体环境影响下的群体行为控制方法，其特征在于，该方法包括以下步骤：

步骤1，初始化流体的水平集、速度，对流体的水平集和速度进行平流：

在初始化流体的物理量之前，需要初始化流体的仿真网格，然后根据网格创建3D纹理，保存流体的初始状态；通过逆向处理和使用隐式积分法计算平流，从每个网格单元逆时追踪粒子轨道，得到它以前的位置，并到开始的网格单元，复制那个位置的量，将这种平流方法应用到速度和水平集上，得到中间的速度和水平集信息；

步骤2，处理静态场景物体，生成静态局部势能网格：

步骤2.1，采用二维向量数组进行存储，数组值代表地理网格顶点上的势能值，向量方向代表势能梯度下降方向，大小代表势能对周围物体的作用力大小；

采用双线性插值进行计算二维空间任意点上的平滑势能作用力；设 和分别代表网格顶点p上的势能；L为网格大小；X和Y分别代表距格子边的距离；则势能如公式(6)所示：

$\stackrel{\→}{F}=(L-X)(L-Y){\stackrel{\→}{p}}_1+X(L-Y){\stackrel{\→}{p}}_2+(L-X)Y{\stackrel{\→}{p}}_3+XY{\stackrel{\→}{p}}_4---\left(6\right)$

步骤2.2，对于规则物体，采用其边沿及其平行扩展作为势能等高线，因此势能方向垂直几何边缘，势能大小采用类似库伦力的方式计算，为了方便描述，引入两个参数---影响范围δ和影响力度λ，影响范围δ表示几何体对其边界以外周围δ单位以内的点存在影响：影响力度λ表示几何体对影响范围以内的点的影响程度比例参数，则几何体外部且距离边缘的最短距离向量为上的点所受的势能如公式(7)计算求得：

$\stackrel{\→}{S}=\mathrm{\λ}(\mathrm{\δ}-|\stackrel{\→}{V}\left|\right)/\left(\right|\stackrel{\→}{V}|\·|\stackrel{\→}{V}\left|\right)\·\stackrel{\→}{V}---\left(7\right)$

当点处于几何体边沿时，势能大小为无穷大；当处于作用范围边缘时，势能大小为0；为了简化考虑，所有处于几何体内部的点所受力大小为无穷大；所有处于物体作用范围δ外的点所受力大小为0；

当一个点受到多个几何体的势能作用时，采用矢量叠加方式处理，从而可以方便地进行几何体的添加和删除操作；

步骤2.3，使用采样方式近似计算地形梯度，设v₁…v₈为网格顶点p周围的八个顶点，将该九个顶点垂直映射到地形上，求得对应8个三维空间向量中与水平面夹角最大的向量V_max；取V_max方向作为梯度V的方向：V_max与水平面夹角的正切值作为梯度V的大小；

地形产生的势能将叠加到已求得的势能当中，假设在地形某点处的地形曲面梯度为V，行走难度为λ，则该点的势能使用公式(8)求得：

F＝-λ*V (8)；

步骤3，根据静态局部势能网格，针对每个网格点，生成对应的全局势能网格：

求解对目标点A的全局导航势能，将全局网格精度设为局部网格的1/3，构建一张路径拓扑图，图中每条线段的权值采用公式(9)计算：

其中F(i)表示i点的局部势能，意思为从顶点A开始，以局部格子大小为采样步长，直到顶点B，所有采样点上的局部势能作用力大小之和，Len(AB)表示顶点A到顶点B的距离，分别代表局部势能和距离对全局势能的影响权值；

当完成对每条线段的权值计算后，就可以求出整个全局势能场，整个运算过程并不求出每条到达A点的单独路径，而是在每个网格点上存储一个指向其最近邻接点的向量作为全局路径导航作用力；

步骤4，处理动态物体，更新周围网格：

动态物体在一定时间内运行的路径会形成一个几何图形，利用该几何图形的静态计算模型即为该物体的动态计算模型；处理动态物体中的行为个体，为了准确反映行为个体间的相互避让作用关系，个体信息被放入周围网格，网格大小一般取行为个体的感知半径大小，行为个体根据其当前所在位置，将其势能信息放入相应格子中，每个行为个体直接从其附近的格子中获悉周围是否有其他个体，从而进行避让；行为个体从其所在格子及周围的九个格子中获得行为个体信息，当行为个体从一个格子走入另外一个格子后，行为个体从前一个格子中注销，然后注册到另外一个格子中；

步骤5，对于路径规划所需要考虑的因素进行计算：

因素主要是指目的地的位置和个体所受的作用力；路径规划是指个体为了到达目的地所选择的路径；目的地的位置使用全局的最小势能表示；如果全局中有若干个目的地，则令这些目的地的势能具有同样的大小，即同时保持全局最小；当群体中个体的目的地信息确定之后，接下来最重要的因素就是群体中个体所受到的作用力；通过个体所受到的作用力，可以得到加速度，从而对于个体的速度产生影响；假设个体是以环境所能允许的最大速度行进，这里的环境是指上下坡的坡度以及人流量的大小因素；计算个体速度如公式(1)所示：

v＝u(x,θ)n_θ (1)

其中，u表示最大速度场，x为个体所在位置，θ表示速度的方向，n_θ＝[cosθ,sinθ]^T表示在θ方向上的单位向量，

上述的个体速度v没有考虑周围环境的影响，这里使用g表示周围环境的影响，即是个体受到周围环境的作用力，对于g来说，它是可以动态变化的，通过动态变化的环境作用力可以帮助个体规避动态的障碍物；

将上述因素结合起来考虑，下面描述个体进行路径规划需要考虑的因素；总的来说，个体会选择到达目的地的最短路径，假定D为从地点x到目的地的路径集合，一个位于x的个体选择的路径P∈D，按照如下公式(2)所示：

L+T+g (2)

其中，L表示路径的长度，T表示所消耗的时间，g表示环境作用力；这里的公式(2)描述了路径规划所需要考虑的三个因素，即路径长度、通过路径的时间和所受的环境作用力，通过这三个因素的权衡，可以得到最终所选择的路径；

步骤6，处理流体环境，计算流体环境对于周围个体产生的作用力，并计算最终的导航力：

导航力是通过对于个体的受力进行分析之后，得到的最终的合力，它会影响个体如何运动，

处理流体环境，计算流体环境对于周围个体产生的作用力，并计算个体最终的导航力时，又进一步包括：

步骤6.1，对流体区域所在的网格进行切片，在每个切片上区分是个体还是流体；对于每一个切片来说，个体区域和流体区域分别使用不同数值进行标记，这样就能区分每一个切片上的个体区域和流体区域；然后，做一次最终遍历，根据每一个切片，把数值复制到流固3D纹理中去；

步骤6.2，对于流体网格进行切片；然后，将网格切片与个体网状模型进行相交操作；在将切片与网状模型相交时，通过将个体网状模型划分成三角形图元来实现；在得到三角形图元与切片的交点之后，通过插值即可得到流体与个体相交的边界点；然后，创建一个新的边界3D纹理，用于存放边界点的信息；在得到边界3D纹理之后，通过对于边界点周围的单元进行判断，得到边界点处的法向，将边界点的法向存储于一个新的3D纹理；

步骤6.3，考虑流体环境对于个体施加的作用力，通过将个体所受的全部作用力进行矢量合成之后，即可得到最终的导航力，流体环境对于个体施加的作用力如公式(3)所示：

$DivAp=Div(-\frac{u^{*}}{\mathrm{\Δ}t}+b)---\left(3\right)$

其中，Div表示流体网格上的离散化的散度求解符号， 对于A和b来说，c是一个3×m的雅克比矩阵，m代表个体模型自由度的数量，令[n₁...n_k]为边界点即流体与个体相交点的法线向量，则M为个体模型的质量矩阵；S为一个k×n的选择矩阵，用来选择耦合处的压力p，k表示流体与个体的相交点的数量，n表示流体单元的数量；向量S_p则是包含了个体周围的所有压力值；S_p中的每一个分量p_i都对个体施加了一个为(Δx)²p_in_i的压力；其中的是矩阵c对于时间的导数，s^*为个体的中间速度；u^*为流体的中间速度场；

根据上述公式(3)，求出流体环境对于个体产生的压力p，也即流体环境对于个体施加的作用力；通过对于个体所受作用力进行合成之后，即可得到最终的个体导航力；

步骤7，根据最终的导航力产生的加速度，得到个体的速度，通过速度来实现个体的路径规划以及规避障碍物：

根据最终的导航力所产生的加速度，可以得到个体的速度，通过速度来实现个体的路径规划和规避障碍物；在公式(2)描述的路径规划所需的开销基础之上，下面给出最优路径的选择方法；假定有一个函数φ，这个函数描述了一条路径；该路径上任意点的开销与最优路径都是一样的；直觉上来说，对于任何个体，最优的策略就是沿着这个函数梯度的相反方向移动，这是降低路径开销最为迅速的方式；事实上，函数φ可以按照下述方式精确定义：在目标点φ＝0，其他的地方φ满足如下方程(4)：

||▽φ(x)||＝C (4)

其中，单位开销C是由梯度▽φ的方向来得到；

令a(x)为位于位置x的个体受到的导航力所产生的加速度，结合公式(1)可以得到速度的计算如公式(5)所示：

$v=-u(x,\mathrm{\θ})\frac{\▿\mathrm{\φ}\left(x\right)}{\left|\right|\▿\mathrm{\φ}\left(x\right)\left|\right|}+a\left(x\right)---\left(5\right)$

其中，v表示速度，u(x,θ)由个体的位置和运动方向而得的最大速度场；由于u与群体密度ρ密切相关；因此，在求解u的时候，根据群体密度的不同，将密度划分为高密度、低密度和中等密度这三种情况来分别进行求解；

步骤8，求解流体环境的压力项，得到流体的最终速度，并调整边界点处的水平集信息：

求解流体环境的压力项，得到流体的最终速度，并调整边界点处的水平集信息；在步骤1中对于流体速度的平流所得到的是中间的散度速度，这个中间速度需要减去压力项的梯度才能得到最终的无散度速度；因此，需要求解压力项；对于压力项的求解采用了Jacobi迭代的方法；这里需要注意的一点是，在步骤1中并没有考虑流体环境与个体的相互作用；由于流体环境对于个体会施加一个作用力，同时个体对于流体环境也会施加一个反作用力；这个反作用力会影响边界点处流体的压力项，因此边界点处流体的速度也会发生改变；同时，需要根据边界点处的流体速度对于水平集进行修正，才能得到最终的水平集信息；

步骤9，对于流体环境与个体的边界进行处理，使两者在边界处满足自由滑动边界条件：

自由滑动边界条件就是流体不能流进个体和从个体中流出，但是允许沿着个体的表面自由流动，这个条件表示流体的速度与个体的速度在边界法向方向上的分量是一样的，如公式(27)所示：

u_fluid·n＝u_solid·n (27)

为了满足自由滑动边界条件，使用的方法是：首先，计算出个体的速度在边界法向方向上的分量；然后，使用这些分量替换相邻流体单元中的法向方向上的分量；通过使用这种方法，在流体环境与个体的边界处满足了自由滑动边界条件。