MIT 6.837：Particle System 粒子系统

在接下来的若干篇博客中，我们将介绍用于实现计算机动画的若干技术。粒子系统是一种基于有限元思想的技术，将复杂的系统转化为有限个粒子模型进行模拟，再利用模拟结果重建出近似的系统模拟结果。

计算机动画的实现方法

关键帧动画

动画师通过指定动画中的部分帧中物体的状态，再通过插值方法生成状态间的过渡帧。

程序生成动画

将物体的运动描述为算法或函数，例如圆周运动用圆的方程表示、球的弹跳用带阻尼衰减的三角函数表示等。

基于物理的动画

记录物体的物理属性，利用物理规律模拟物体接下来的运动状态。这种方法与现实的物理规律相符，但难以控制其动画效果，模拟过程计算量开销也很大。

动作捕捉

在真实的物体上粘贴光学标记，再使用高速摄像机拍摄真实世界中的动画，利用这些数据在计算机中重建运动过程。常用于制作需要自然感或细致的动画，但需要在现实中利用设备捕捉动作，成本和难度较高。

粒子系统

模拟物理的过程中，常用的物理元有质点、刚体和变形物体。粒子系统就将要模拟的系统看做一系列质点组成，并模拟这些质点之间的物理规律，从而得到系统结果的近似解。

粒子系统通常生成许多粒子，用空间中的力场描述外界对粒子施加的力，并应用物理学中的规律求解例子状态。这些状态可以包括空间位置、速度、加速度、质量、生命周期等。

粒子的运动

根据经典的牛顿定律，运动涉及的参数包括质量 $m$、位置 $\boldsymbol{x}$、速度 $\boldsymbol{v}$ 和受力 $\boldsymbol{F}$。通过牛顿第二定律，它们建立起了这样的一组关系

$$ \begin{aligned}
& \frac{\mathrm{d}}{\mathrm{d}t} \boldsymbol{x}(t) = \boldsymbol{v}(t) \\
& \frac{\mathrm{d}}{\mathrm{d}t} \boldsymbol{v}(t) = \frac{1}{m} \boldsymbol{F}(\boldsymbol{x}, \boldsymbol{v}, t)
\end{aligned} $$

这是一组常微分方程（Ordinary Differential Equation，ODE），记 $\boldsymbol{X}(t) = \begin{pmatrix} \boldsymbol{x} \\ \boldsymbol{v} \end{pmatrix}, \boldsymbol{f}(\boldsymbol{X}, t) = \begin{pmatrix} \boldsymbol{v} \\ \dfrac{1}{m} \boldsymbol{F} \end{pmatrix}$，则可以表示为规范化的写法

$$ \frac{\mathrm{d}\boldsymbol{X}}{\mathrm{d}t} = \boldsymbol{f}(\boldsymbol{X}, t) $$

ODE 的解是一个困难的问题，有些 ODE 可以利用通解公式简单求值，但也有许多 ODE 没有解析解。而对于动画中的粒子运动，我们其实只是想知道某些 $t$ 值上的解，因此不需要关注它是否存在解析解，而只需要一个数值解，甚至一个近似解，就可以实现动画效果了。

微分方程的数值解

对于计算机动画来说，我们只需要知道在一些离散时刻状态的值，再将其渲染为帧即可。因此，ODE 的解析解不是必要的，我们只需要有办法得到这些离散时刻的数值解即可。而由于帧间时间间隔很小，利用导数估计下一帧的状态是一个可能的思路。

Euler 方法

欧拉方法直接将时间间隔与导数相乘，即

$$ \begin{aligned}
& t_1 = t_0 + h \\
& \boldsymbol{X}_1 = \boldsymbol{X}_0 + h\boldsymbol{f}(\boldsymbol{X}_0, t_0)
\end{aligned} $$

此处时间间隔 $h$ 取得越小，估计的错误则越小。由于在整个迭代过程中错误会积累，如果 $h$ 的值取的不合适，得到的结果可能与实际结果相去甚远。另外，Euler 方法存在一些无论 $h$ 取得多小都会出错的情况，例如圆周运动。

Euler 方法的收敛性也无法保证，例如对于下面这个例子：

这些错误是由于 Euler 方法是一阶估计方法，其错误为 $O(h^2)$，且在整个迭代过程中错误会积累起来。

二阶方法

一阶方法的准确度不佳，因此我们可以使用二阶估计方法，降低 $h$ 取值对估计准确度的影响，使误差缩小到 $O(h^3)$。

改进 Euler 方法：利用一阶估计终点处导数的值补偿估计误差，令 $\boldsymbol{f}_0 = \boldsymbol{f}(\boldsymbol{X}_0, t_0), \boldsymbol{f}_1 = \boldsymbol{f}(\boldsymbol{X}_0+h\boldsymbol{f}_0, t_0+h)$，则 $\boldsymbol{X}(t_0+h) = \boldsymbol{X}_0 + \dfrac{h}{2}(\boldsymbol{f}_0+\boldsymbol{f}_1)$。

中点方法：利用步长中点处导数的值补偿误差，令 $\boldsymbol{f}_0 = \boldsymbol{f}(\boldsymbol{X}_0, t_0), \boldsymbol{f}_m = \boldsymbol{f} \left( \boldsymbol{X}_0+\dfrac{h}{2}\boldsymbol{f}_0, t_0+\dfrac{h}{2} \right)$，则 $\boldsymbol{X}(t_0+h) = \boldsymbol{X}_0 + h\boldsymbol{f}_m$。

改进 Euler 方法使用两个全步长估计的均值，而中点方法则是先取半步长导数再估计全步长，两者的结果不一定相同，但具有同阶数的精确度。

粒子动画的流程

伪代码如下：

AnimateParticles(n, y0, t0, tf) {
    y = y0
    t = t0
    DrawParticles(n, y)
    while(t != tf) {
        f = ComputeForces(y, t)
        dydt = AssembleDerivative(y, f)
            //there could be multiple force fields
        {y, t} = ODESolverStep(6n, y, dy/dt)
        DrawParticles(n, y)
    }
}

粒子的外部受力

重力：$\boldsymbol{F} = \begin{pmatrix} 0 \\ 0 \\ -mg \end{pmatrix}$，很轻的粒子（烟尘等）和火焰可以忽略重力作用

阻力：$\boldsymbol{F} = -d \boldsymbol{v}$，可以利用阻力让数值解稳定下来，也可以模拟在粘稠液体中的运动效果

空域力场：力与空间坐标相关，例如风推力/阻力、电磁力、漩涡等，有时需要改变物体的能量则速度也作为参数

相互作用力：与两粒子的空间坐标相关，例如用于模拟液体内部的作用力 $\boldsymbol{F} = \dfrac{k_1}{\left| \boldsymbol{x}^{(i)} – \boldsymbol{x}^{(j)} \right|^m} – \dfrac{k_2}{\left| \boldsymbol{x}^{(i)} – \boldsymbol{x}^{(j)} \right|^n}$（Lennard-Jones 力）

粒子的碰撞

粒子间的碰撞通常忽略不做处理，一般只做与环境的碰撞检测。碰撞处理大致可以分为两种，一种为模拟真实的反弹，一种只是将其移出物体表面。与平面的求交可以类似之前射线的方法，计算平面方程的值即可。

在碰撞后，其速度变化规律有多种处理方式。理想反弹时，其切向速度不变，而法向速度反向，也即 $\boldsymbol{v}’ = \boldsymbol{v}_t – \boldsymbol{v}_n$。也可以引入系数 $\varepsilon$ 反应非理想反弹时法向速度的损失，即 $\boldsymbol{v}’ = \boldsymbol{v}_t – \varepsilon \boldsymbol{v}_n$。

粒子的生命周期

利用生命周期，我们可以控制粒子实例的数量，移除超出屏幕或世界范围的粒子等。此外，还可以用生命值来实现颜色、透明度等的改变，实现变色、淡出等效果。

粒子的应用

除粒子模拟火焰、粉尘等效果之外，粒子还可用于实现其他一些效果。

粒子建模

模拟鸟群、人群

游戏中的系统和效果

天气、爆炸、草丛等。