CN111460741A - 一种基于数据驱动的流体模拟方法 - Google Patents

Abstract

一种基于数据驱动的流体模拟方法，其较准确地匹配了两个流体的表面,同时避免了误差的不断积累，能够更好地恢复流体表面之间的动作，能获得更自然的流体表面，在自动匹配和获得适定的对应关系方面均比用户参与的配准方法以及第三方数据点的方法更好。该方法包括：(1)通过两组不同的设置参数生成两个流体序列，然后转化成两组时空流体表面φ₁和φ₂作为输入数据；(2)进行流体动态模型的配准，添加基于L0范数的光滑约束条件，通过限制变形场中所包含的非零梯度的数量实现变形场的光滑化处理；(3)进行变形的对齐，通过残差迭代以及分辨率的分层逐渐优化；(4)进行表面投影；(5)进行插值运算并输出结果。

Description

一种基于数据驱动的流体模拟方法

技术领域

本发明涉及空间动画变形的技术领域，尤其涉及一种基于数据驱动的流体模拟方法。

背景技术

利用计算机模拟具有真实感的现实世界场景受到众多研究者重视，其中关于流体运动的模拟在电影制作、游戏的实时渲染、虚拟现实等方面有着越来越广泛的应用。流体模拟算法主要分为两大类：基于物理的和数据驱动的模拟方法。

基于流体的物理属性去模拟流体运动受流体动力学中经典的纳维-斯托克斯(Navier-Stokes)方程组的启发，依据流体的实际运动规律，能获得任意初始设置下的流体速度场，生成富有真实感的流体运动过程。目前的流体模拟方法都包含了许多重要的设置参数，如果想要模拟具体的流体效果则需要不断地调整一系列的参数：当生成的流体动画不符合目标时，设置参数进一步被修改以使新的流体动画更接近目标。基于物理的流体模拟通常需要处理大量的流体数据，而且经过大量的计算之后还要对流体数据进行渲染，任何参数设置出现问题都需要再重新执行整个模拟算法。

由于那些不符合要求的流体数据被轻易地丢弃了，Thuerey等人提出了一种基于数据驱动的流体模拟方法将这些数据重新利用起来。该模拟方法只要匹配两个流体动画，获得两者之间的对应关系，就能高效地生成以两个流体动画为边界的一系列中间运动状态，并且不需要重新运行算法。光流匹配模型可实现源流体表面至目标流体表面的形态转换：首先要计算合适的变形场，寻找源流体表面与目标流体表面之间的对应关系从而对齐两组流体表面，然后才能利用变形场来融合流体表面并模拟新的流体运动。

下面对光流技术进行简单的介绍。包含运动物体的三维场景投影到二维平面上形成了一系列连续变化的图像，光流描述的正是图像序列中像素点的瞬时变化速度，其中就包含了目标物体的运动信息。具体来讲，图像上每一点随着时间的变化而移动，因此二维平面上的像素点具有速度场，可以近似地表示空间中运动物体的运动场。用像素点的亮度变化表示像素点的运动，光流场可以用来反映图像序列中每个点亮度的变化。光流技术利用亮度随时间的变化信息以及相邻帧之间的对应关系计算空间运动物体的运动情况。亮度恒定是光流技术的基本假设，该假设认为图像序列中某一像素点在不同帧之间可能存在位置变化，但该点亮度值的大小是不会变化的。以(x(t),y(t))表示空间中的一点投影到图像平面的运动轨迹，这一点在时间t的亮度用I(x(t),y(t),t)来表示，亮度一致假设可以用以下公式表示：dI/dt＝0。用速度场u(x,y)来表示图像中每个点的运动状态，应用链式法则可以推出

公式(1)中，

这就是光流的约束方程。

仅利用光流约束方程无法得到唯一的解，因此添加一些正则项约束变形场则能够获得更加合适的对应关系。在匹配源流体表面与目标流体表面的光流约束的基础上，L2范数的光滑正则化和吉洪诺夫正则化可以作为新的约束项来更好地匹配流体表面。不过，虽然基于L2范数的光滑正则化作用于梯度值上能够压制变形场相邻点之间的差异，但是对于有较大差异的源流体表面和目标流体表面却不一定能获得理想的变形。液滴落入水池的例子，初始状态为不同形状的液滴从不同的位置落下。其中，方形液滴的对应点应该位于球形液滴而不是水池液面，然而基于L2的光滑正则项的约束并没有很好地优化未对齐的表面，得到了不自然的结果。这是由于匹配模型属于最小化问题，主要采用的是容易受限于局部最小化的最近距离思想。

由于表面某些区域上的点通常与周围的点有着相似的运动变化，所以对应的变形场呈现出了稀疏特性。利用这种稀疏特性最直接的方法就是将基于L0范数的正则项加到光流匹配模型中。L0范数的最小化问题是典型的凸优化问题，不能直接进行求解。图像光滑处理中运用了基于L0的梯度最小化方法，其提出的交替方向迭代法将最小化问题分解为两个子问题，能很好地对这一凸问题近似求解。实际上，基于L0范数的正则化已经在图像处理、计算机视觉以及图形学等领域中有着广泛的应用。光流估计以及动作重建等课题中也很重视稀疏先验的作用。在Shen等人中，利用稀疏表示作为先验信息估计流场，光流估计模型结合了这种针对流场梯度值的稀约束项来处理测量噪声。Chen等人为了抑制过稀疏的问题而将运动场分为稀疏与非稀疏两部分。Guo等人在动作跟踪中将基于L0范数的动作正则项用于关节处的动作估计，从而提出了基于L2-L0的非刚性动作跟踪算法，并且有效地减少了误差的传播。

发明内容

数据驱动流体模拟方法在匹配流体表面时的变形场存在误差,为克服现有技术的缺陷，本发明要解决的技术问题是提供了一种基于数据驱动的流体模拟方法,该方法利用L0正则化约束，能较准确地匹配了两个流体表面,同时避免了误差的不断积累，能够更好地恢复流体表面之间的动作，能获得更自然的流体表面，在自动匹配和获得适定的对应关系方面均比用户参与的配准方法以及第三方数据点的方法更好。

本发明的技术方案是：这种基于数据驱动的流体模拟方法，该方法包括以下步骤：

(1)通过两组不同的设置参数生成两个流体序列，然后转化成两组时空流体表面φ₁和φ₂作为输入数据；

(2)进行流体动态模型的配准，添加基于L0范数的光滑约束条件，通过限制变形场中所包含的非零梯度的数量实现变形场的光滑化处理；

(3)进行变形的对齐，通过残差迭代以及分辨率的分层逐渐优化；

(4)进行表面投影，为了恢复两组时空流体表面的细节，当光流匹配算法计算出表面SDF的变形场之后，通过线性查找的方式来获取更多细节的变形；

(5)进行插值运算并输出结果，利用变形场在流体表面φ₁与流体表面φ₂之间建立对应关系，将该变形场应用于插值运算生成不同流体表面间光滑且流畅的转换过程，插值运算结合两组流体表面SDF的并集，从而在结果中保留小液滴以及薄膜结构。

本发明通过两组不同的设置参数生成两个流体序列，然后转化成两组时空流体表面φ₁和φ₂作为输入数据，通过结合L0范数正则化的光流匹配模型，限制变形场中所包含的非零梯度的数量实现变形场的光滑化处理，然后对齐已生成流体动态模型的表面，进行表面投影，通过插值运算将两者融合，进一步生成多种新的流体运动，优化后的匹配模型较准确地匹配了两个流体表面,同时避免了误差的不断积累，能够更好地恢复流体表面之间的动作，能获得更自然的流体表面，在自动匹配和获得适定的对应关系方面均比用户参与的配准方法以及第三方数据点的方法更好。

附图说明

图1是根据本发明的基于数据驱动的流体模拟方法的整体流程图。

图2是根据本发明的基于数据驱动的流体模拟方法的一个详细实施例的流程图。

图3是根据本发明的分层方法的示意图。

图4是根据本发明的变形场对齐的示意图。

具体实施方式

如图1所示，这种基于数据驱动的流体模拟方法，该方法包括以下步骤：

(1)通过两组不同的设置参数生成两个流体序列，然后转化成两组时空流体表面φ₁和φ₂作为输入数据；

(2)进行流体动态模型的配准，添加基于L0范数的光滑约束条件，通过限制变形场中所包含的非零梯度的数量实现变形场的光滑化处理；

(3)进行变形的对齐，通过残差迭代以及分辨率的分层逐渐优化；

(4)进行表面投影，为了恢复两组时空流体表面的细节，当光流匹配算法计算出表面SDF的变形场之后，通过线性查找的方式来获取更多细节的变形；

(5)进行插值运算并输出结果，利用变形场在流体表面φ₁与流体表面φ₂之间建立对应关系，将该变形场应用于插值运算生成不同流体表面间光滑且流畅的转换过程，插值运算结合两组流体表面SDF的并集，从而在结果中保留小液滴以及薄膜结构。

本发明通过两组不同的设置参数生成两个流体序列，然后转化成两组时空流体表面φ₁和φ₂作为输入数据，通过结合L0范数正则化的光流匹配模型，限制变形场中所包含的非零梯度的数量实现变形场的光滑化处理，然后对齐已生成流体动态模型的表面，进行表面投影，通过插值运算将两者融合，进一步生成多种新的流体运动，优化后的匹配模型较准确地匹配了两个流体表面,同时避免了误差的不断积累，能够更好地恢复流体表面之间的动作，能获得更自然的流体表面，在自动匹配和获得适定的对应关系方面均比用户参与的配准方法以及第三方数据点的方法更好。

优选地，所述步骤(1)中，通过流体模拟算法获得的流体动画由一系列连续的流体表面序列组成，流体的每个微元在某一序列帧的时刻具有相应的瞬时速度信息，一个流体动画经过变形之后得到的结果仍然保持时空相关性,根据流体运动的速度场将流体动态模型转化成包含时间和空间维度的表面数据，使得在对动态数据进行变形优化时有更多的参考信息，流体模拟算法将利用不同的设置参数生成两组流体序列，然后将流体数据转化成两组时空表面φ₁和φ₂作为输入数据。

优选地，所述步骤(2)中，通过计算4D稠密变形场u将流体表面φ₁的点映射到流体表面φ₂，以使φ₁的表面形态变形为φ₂，通过光流技术计算两组流体动态模型之间的动作，流体模拟法的设置参数变量r作为导致两组流体表面之间出现区别的自变量，表示一种外部的控制参数，由光流的亮度恒定假设得出dφ/dr＝0，应用链式法则推出流体动态模型的配准方程表示为公式(2)：

其中，流体表面上点的位置变化dx/dr表示表面φ₁到φ₂的变形u，用表面φ₂对时间和空间的导数来近似表示

优选地，所述步骤(2)中，基于L0范数的光滑约束条件的光流匹配模型的能量方程为公式(9)：

其中，参数β_T用于控制吉洪诺夫正则项的作用在整个方程中所占的影响，

自适应参数β_A用于控制点p的梯度

与其对应的辅助变量v^p之间的相似性，参数β_S用于控制变形场的光滑程度。

优选地，所述步骤(2)中，能量方程包括两种不同性质的子问题，分别对其建立了闭合解的模型；

子问题1：v的最小化

将未涉及到变量v的项去掉，根据公式(10)得到关于v的最小化模型

如果把点p的梯度

作为已知量，对应变量v^p可以通过条件方程来进行确定，该方程的表示为公式(11)：

计算出流体表面上每个点p的最小能量，然后根据公式(12)求和就得到了子问题1的最小能量

子问题2：u的最小化，如果变量v已经确定，变形u的计算模型通过公式(13)表示：

该公式属于二次方程的最小化问题，利用梯度下降法直接求解，求E(u)关于变形量u的导数并令其为0，得到的线性方程表示为公式(14)：

其中，用两组流体表面SDF的差φ₁-φ₂来表示梯度φ_r，作为变形目标的流体表面的空间梯度表示

L为拉普拉斯算子。

优选地，所述步骤(3)中的残差迭代，通过光流技术求解变形场，然后对流体表面φ₁进行变形，并判断其是否与流体表面φ₂的形态更相似；如果相似那么再次执行算法来改进变形场，重复多次直到算法收敛。

优选地，所述步骤(3)中的分辨率的分层逐渐优化，采用金字塔结构，其中五角星代表流经某一流体微元的表面，他的位置从流体表面φ₁的左下方变形到流体表面φ₂的右上角，在红色方框标注的小窗口范围内寻找当前表面移动位置，对流体数据进行下采样得到粗分辨率的数据；在粗分辨率这一层得到的变形场记作u¹，然后传递至下一层并在此基础上继续修正变形得到u²，直至最后一层从而实现变形场由粗到精的分层迭代运算。

优选地，所述步骤(4)中，变形表面φ₁已经非常接近目标φ₂的情况下，沿着目标表面上某一点的法线方向找到在表面φ₁中最近的交点，然后通过缩减两者差值的方式改进变形场；当光流匹配计算的变形场使流体表面之间的差距越小，这种投影方法越稳定地恢复表面的小变形。

优选地，所述步骤(5)中，流体表面φ₁由变形场u_1→2变形得到φ′₁，由变形场u_2→1得到表面φ′₂，已知插值权重系数为α，插值运算的表达式为公式(15)、(16)：

φ′_1∪2＝min(φ′₁,φ′₂) (15)

φ_output＝ω₁φ′₁+ω_1∪2φ′_1∪2+ω₂φ′₂ (16)

φ_output表示插值运算生成的流体表面，插值的权重系数表示为公式(17)-(19)，其中clamp用于约束参数范围至0到1之间：

ω₁＝clamp(1-2α) (17)

ω₂＝clamp(2α-1) (18)

ω_1∪2＝1-ω₁-ω₂ (19)。

以下更详细地说明本发明的内容。

给定的3D流体模型能够基于流体速度场得到具有连续运动效果的流体动画，动画序列在每一帧的水体形态构成了随时间变化的流体动态模型。数据驱动的流体模拟算法采用表面匹配的方法，建立两种流体动态模型的表面在同一时间的对应关系，并计算合适的变形场去对齐流体表面，然后利用插值运算进一步扩展出新的流体运动。为了对齐不同的流体表面需要添加正则化惩罚项，从而改善变形场的过拟合问题。本发明提出了一种结合L0范数正则化的光流匹配模型，该模型采用基于L0的梯度最小化算法解决光流估计中的光滑化处理问题，使得由于错误的对应关系造成的变形量误差累积以及流体运动不连续等情况得到有效缓解。为了实现L0范数在匹配模型中的稳健求解，目标函数的最小化问题被转化成两个子问题的最小化，并运用交替最小化得到变形量的最优解。

基于物理的流体模拟算法能够根据给定的流体模型获得相应的流体速度场，从而模拟流体运动。但是这种方法对于复杂的流体模拟无法达到实时效果，所以为了快速模拟流体运动可以采用数据驱动的流体模拟。其中基于模型数据的方法是对齐已生成流体动态模型的表面，然后通过插值运算将两者融合，进一步生成多种新的流体运动。

图2是根据本发明的基于数据驱动的流体模拟方法的一个详细实施例的流程图。

1.作为输入的流体动画由基于物理的流体模拟算法所生成。我们首先通过设置流体的初始状态，调整流体模拟的设置参数如液体的粘性系数、液滴的位置、对应场景下的重力大小等生成两组相似的流体运动序列，使得生成的两组流体之间存在一些差异但又具有相似性。匹配具有相似运动的两组流体动画，与光流技术中提到的运动目标检测的思路相类似。

2.采用有向距离场描述流体表面。光流技术中依靠的仅仅是图像中的亮度信息，如果当前区域的像素点亮度信息一致，那么该区域当中的某一像素点就存在多种可能的速度矢量，这就是亮度恒定假设引起的孔径问题。这种运动模糊的问题很大程度上是由于场景中的各种复杂因素给亮度信息带来的影响所造成的，因此降低输入信息的不确定性才能从根本上解决孔径问题。有向距离场是通过流体表面上每个点与已知的零等值面之间的最短距离来表示流体表面的，因此整个流体表面上不存在信息不充分或者模糊的情况。运用光流技术进行动作估计是利用输入数据的梯度实现的，与直接利用流体数据的密度等信息相比有向距离场能够获得光滑且准确的空间梯度。为了利用光流技术获得高质量的变形场，将流体表面转化为采用有向距离场来描述是非常关键的一步。

3.整个光流匹配算法需要完成多次迭代，逐渐优化变形场。第一次迭代的变形场初始化为0。

4.利用当前的变形场u，对流体表面序列φ₁进行变形，得到新的流体表面φ′₁。

5.由于能量方程中包含有基于L0范数的正则项，所以能量方程的求解被分成两个子问题。辅助变量v被引入能量方程中，其中v与变形场梯度

相对应，利用闭合公式(11)计算v。

6.子问题二是通过梯度下降法最小化包含变形场u的能量方程(13)，然后求解能量最小时的变形量u⁽ⁱ⁾。

7.将上一步得到的变形量u⁽ⁱ⁾与上一次迭代中得到的变形量对齐，生成一个光滑的变形量u^tmp(i)。

8.利用变形场u^tmp(i)对流体表面φ₁进行变形得到新的流体表面，然后与流体表面φ₂进行比较；同时利用上一次迭代得到的变形流体表面与流体表面φ₂进行比较。分别采用公式(20)提到的误差测量方程得到二者的误差值。

9.如果变形场u^tmp(i)能够使得最终的误差值更小，那么用κ＝2乘以参数β_A以加速算法收敛，然后继续执行步骤4。

10.如果变形场u^tmp(i)无法使最终的误差值更小，那么就得到了最优的变形场，然后保存变形场。

11.基于步骤10获得的变形场生成了新的流体表面φ′₁，然后通过线性查找的方式来对齐φ′₁与φ₂之间的细节变形。

12.通过光流匹配以及投影之后得到优化的变形场u，并生成合适的对应关系，然后利用插值公式(16)融合两组流体表面序列，生成新的流体动画序列。

通过流体模拟算法获得的流体动画由一系列连续的流体表面序列组成，流体的每个微元在某一序列帧的时刻具有相应的瞬时速度信息。如果根据流体运动的速度场将输入的流体表面SDF按照时间的顺序连接起来，那么流体动态模型将被转化为包含时间与空间维度的表面数据。这种时空流体表面数据包含了每个点在时间维度上的运动信息，因此在对流体动画这种动态数据进行变形优化的时候提供了更多的参考信息来获得最佳匹配。流体模拟算法将利用两组不同的设置参数生成两个流体序列，然后转化成时空表面作为融合模型的输入数据。融合模型主要包括四个步骤：流体动态模型的配准、多组连续欧拉变形的对齐方法、表面投影以及插值运算。

流体动态模型的配准。在这里将作为输入的两组时空流体表面分别用φ₁和φ₂来表示，配准的目标主要是通过计算4D稠密变形场u将流体表面φ₁的点映射到流体表面φ₂，以使φ₁的表面形态变形为φ₂。接下来的内容都是以流体表面φ₂为目标，计算流体表面φ₁的变形场，当然也可以改变两组输入的角色计算反方向的变形。光流技术通常用于恢复在不同时间获得的图像中动态物体的动作，同样地可以用于计算两组由不同模拟生成的时空表面之间的变形。两组流体表面之间的变化可以统一地用设置参数r来表示，作用类似于光流技术中的时间变量t，属于一种外部的控制变量。对于一组关于设置参数r发生运动变化的时空表面φ(x(r),r)，上述的亮度恒定假设约束表面的运动变化可以表示为dφ/dr＝0，应用链式法则推出配准时空表面的约束方程表示如下：

式中，时空表面的位置变化dx/dr表示表面φ₁到φ₂的变形u，用表面φ₂关于时间和空间的导数来表示

由于一组时空表面φ之间的变化源于设置参数r的变化，因此利用

表示时空表面关于参数r的偏导数。如果改变输入的流体数据那么相应的设置参数r也随之发生变化，从而说明这个方程高效地参照了五个维度的变化来计算变形场。

将光流技术应用于数据驱动流体模拟能有效地计算时空表面的变形场并对齐表面，不过仅依靠由光流亮度恒定假设推导的约束方程不能求得唯一解，因为配准时空表面的约束方程(1)具有非线性特性，此时必须引入新的约束项。

众所周知，场景中某一物体的运动变化应该是连续且平滑的，这也正是光流亮度恒定假设的依据。以时空表面上某一点为例，该点与邻域内的这些点应该具有相似的速度矢量，从而在给定的空间范围内，速度场保持一种光滑的变化状态。因此，基于前面描述的光流约束方程，可以添加一个约束项用于惩罚“非光滑解”。常见的基于L2范式的光滑约束是限制时空变形场中某一点的变形量与周围变形量之间的差值，即通过最小化该点的变形量沿着各维度的偏导数平方和来确保变形场的光滑性。然而，对于两组形态并不相似的流体表面序列，由基于L2范数的光滑约束优化后的光流匹配模型所得到的变形场依然不能准确地匹配两组表面，因为L2范数采用的是最近邻思想，流体表面φ₁上的某个点在查找关于流体表面φ₂的对应点时只在该点邻近的位置处而不会考虑整个流体表面φ₂。

实际上在对齐两个时空表面时，由于每个点的变换不是独立进行的，除了某些差别较大的区域以外其他区域的变形量梯度为零，这就表明时空变形场应该具有稀疏性。因此，本文提出了一种基于L0范数的光滑约束条件，通过限制变形场中所包含的“非零梯度”的数量实现变形场的光滑化处理。这种梯度计数方法可以从整体上平滑化变形场并且很好地保留细节性变形，从而防止流体表面上的点被移动到错误的地方。接下来，通过最小化能量方程的方法求解时空变形场u，该能量方程包括光流约束项以及关于变形场梯度的光滑约束项：

式中参数β_S用于控制时空变形场的光滑程度。下面先利用基于L2范数的光滑约束求解E_smooth(u)，并利用最小二乘思想重新表述该能量方程，其具体表示如下：

E_total(u)中的第一项是光流约束项与公式(1)相对应，u是由空间上的三个维度以及时间维度这四个部分组成的稠密时空变形场，公式中的m表示这四个维度。那些一定范围内的表面通常具有相似的运动形变，因此对应的变形场梯度呈现出了稀疏的特性，本文利用变形场的稀疏性这一优势，采用了L0范数求解能量方程中的光滑约束，用E′_smooth(u)代替E_smooth(u)实现基于L0范数的光滑约束，能量最小化表达式如下所示：

光滑化通常约束的是变量的梯度，时空变形场中每个点p与关于时间和空间方向相邻的点之间的差值用梯度

表示，用E′_smooth(u)测量梯度

的数量，其具体的测量表达式可写为：

E′_smooth(u)＝S(u) (6)

用

表示变形场中点p的梯度值，S(u)统计了其中非零梯度的个数。另外，为了限制变形向量的大小，方程中还添加了一项有利于小变形存在的吉洪诺夫正则化。这样，在光流约束的基础上加入光滑正则化和吉洪诺夫正则化，通过最小化能量方程求解变形场u。本文提出的基于L0范数正则化的光流匹配模型具体表示如下：

式中，参数β_T用于控制吉洪诺夫正则项的作用在整个方程中所产生的影响。

通常来说，复杂的变形场中除了变化相似的区域外还存在类似于运动变化的边界或者非连续性运动这样的区域，这种地方的变形梯度比较大，也就是说变形场不是完全光滑的，因此合适的变形场正则化对于运动估计的准确性至关重要。很多运动估计模型采用的是常见的基于L2范数的光滑约束来正则化变形场，但是对这种不连续效果的优化并不好。由此可以看出，计算变形场时在以亮度恒定假设作为约束的基础上还应该加入稀疏的光滑正则项作为先验信息。这样，与单独对那些不连续变形效果建立对应关系相较而言，利用稀疏的光滑约束能有效地光滑化运动一致的区域同时还能够优化那些具有不连续性动作的区域，从而获得更准确的变形场。

通过观察这个能量方程可以发现，整个模型中同时包括了两种不同类型的解决办法，前两项的计算是以时空表面上的点为单位进行求解的，而最后一项的L0范数是一种统计型的计数模型。这种凸优化问题是无法直接求解的，因此本文通过引入对应于梯度

的辅助变量v，运用基于半二次型分解的交替方向迭代优化方法近似求解最小化问题。重新描述模型并表示如下：

上式中

自适应参数β_A用于控制点p的梯度

与其对应的辅助变量v^p之间的相似性。该模型中包括了两种不同性质的子问题，本文分别对其建立了闭合解的模型。

子问题1：v的最小化将未涉及到变量v的项去掉，得到关于v的最小化模型

如果把点p的梯度

作为已知量，对应变量v^p可以通过条件方程来进行确定，该方程的表示如下：

这样，计算出时空表面上每个点p的最小能量，然后求和就得到了子问题1的最小能量：

子问题2：u的最小化同样的，如果变量v已经确定，变形u的计算模型可以用一个二次方程来表示：

该公式属于二次方程的最小化问题，利用梯度下降法可以直接求解。求E(u)关于变形量u的导数并令其为0，得到关于变形量u的线性方程表示如下：

其中，用两组时空表面SDF的差φ₁-φ₂来表示梯度φ_r，作为变形目标的时空表面的空间梯度用来表示

L为拉普拉斯算子。

变形的对齐。流体运动的变形场不可能仅通过一次求解就得到，还需要通过残差迭代以及分辨率的分层逐渐优化。首先，通过光流技术求解变形场，然后对流体表面φ₁进行变形，并判断其是否与流体表面φ₂的形态更相似。如果相似那么应该再次执行算法来求解变形场，这一过程需要重复多次直到算法收敛。传统的光流法对于估计变化较小的运动有较高的准确性，但是对于短时间内运动剧烈的物体就很难检测到了。类似地，如果流体运动的范围非常小，仅有1到2个流体微元的距离，变形场能很好地被光流技术恢复，但实际上流体运动的复杂性说明运动距离不仅限于几个微元，因此实现分辨率的分层是很有必要的。分辨率分层就是通过重采样当前的流体数据获得更低分辨率的流体数据，这样在原先的分辨率时无法检测到的大变形可能变成了所谓的小变形。

图3描述了分层方法采用的金字塔结构，其中五角星代表流经某一流体微元的表面，他的位置从流体φ₁的左下方变形到流体φ₂的右上角，但是由于亮度恒定假设导致的孔径问题，目前只能在加粗的方框标注的小窗口范围内寻找流体φ₁的表面是如何移动的，流体φ₁的实际运动路径就无法正确地估计。在这种情况下，对流体数据进行下采样获得低分辨率的流体数据，能更明确时空表面变形量的大小以及方向。由低分辨率的流体数据得到的变形场记作u¹，然后传递至更高分辨率的数据并在此基础上继续修正变形场得到u²，直至最后一层从而实现变形场由粗到精的分层迭代优化。

为了得到最优的变形场而进行了多次计算，得到与迭代顺序相关的变形场序列u^1…n。如果每迭代一次都要将变形场存储起来并结合前几次计算的变形场依次变形流体φ₁，运算量比较大。因此获得一个光滑的线性变形场不仅能够更好地配合迭代运算，而且有利于算法的模块化运行。但是由于本文在描述流体运动时采用的是欧拉描述，变形向量位于每个流体微元体上而不在流体的表面上，只能通过变形向量反推出哪个位置的表面需要移动到变形向量所处微元体的位置，所以合并变形场序列不能通过直接相加来获得，否则就会产生错误的结果，如图4第二行所示。本文假设图4所描述的例子中，变形矢量的大小均为1，仅用箭头表示其方向信息，其它的微元体中没有标注箭头是因为没有对当前示例中的流体表面产生影响。第二行描述的错误结果中忽略了变形矢量u²的作用，所以应该在具有变形矢量u²的微元体位置结合u¹，从而获得一个正确的位于表面目标位置的变形矢量。

表面投影。经过光流匹配模型的计算获得了两个时空表面序列之间的变形场，不过此时的变形场只是简单地恢复了表面φ₂的大致形态，而流体表面的一些细节性变形还没有被检测到。为了从细节上恢复时空表面，在表面SDF由光流匹配得到的变形场处理之后，可以通过线性查找的方式来获取更多细节性的变形。

形变表面φ₁已经非常接近目标φ₂的情况下，沿着目标表面上某一点的法线方向，即SDF的梯度方向找到与表面φ₁最近的交点，然后通过缩减两者差值的方式改进变形场。光流匹配计算的变形场越能使流体表面之间的差距缩小，这种投影方法就越是能稳定地恢复表面的小变形。

插值运算。计算流体序列之间的变形场并逐渐对其进行优化之后，就可以建立流体表面φ₁与流体表面φ₂之间的对应关系，将该变形场应用于插值运算则可以生成不同流体时空表面间光滑且流畅的变换过程。由于流体运动中其表面形态常伴随着小液滴以及薄膜结构的特点，这些细小的形态特征在插值过程中很容易丢失，所以插值运算还要结合两组时空表面SDF的并集从而在结果中保留这些细节。流体表面φ₁由变形场u_1→2变形得到φ′₁，同样地由变形场u_2→1得到表面φ′₂，已知插值权重系数为α，插值运算的表达式如下：

φ′_1∪2＝min(φ′₁,φ′₂) (15)

φ_output＝ω₁φ′₁+ω_1∪2φ′_1∪2+ω₂φ′₂ (16)

φ_output表示插值运算生成的流体时空表面，插值的权重系数表示如下：

ω₁＝clamp(1-2α) (17)

ω₂＝clamp(2α-1) (18)

ω_1∪2＝1-ω₁-ω₂ (19)

式中的clamp用于约束权重系数的范围至0到1之间。

下面具体描述基于L0正则化的光流匹配算法，如表1所示。算法的输入由两组时空表面SDFφ₁和φ₂以及初始化为0的变形场u组成。函数deform(φ,u,α)的作用是采用变形场u来变换流体φ的形态，表面的变形相当于流体模拟中的水平对流过程，参数α用来控制变形的强度。最终的变形场u需要经过多次迭代运算的优化，流体数据融合模型中介绍的对齐方法alignVelocity(u)在这里发挥了很大的作用，极大地降低了算法的复杂性。

表1

本发明的有益效果如下：

利用流体模拟框架mantaflow进行了流体动画的生成、流体表面数据的匹配以及流体的插值运算。通过表面匹配获得的变形量与输入的两组流体表面共同作为插值运算的输入，然后生成新的流体动画。其中，源流体表面与目标流体表面序列是由基于物理的流体模拟算法生成的。插值运算生成了一组介于源与目标流体表面之间的过渡状态的流体表面动画。分别用本发明提出的基于L0的流体模拟方法以及传统的基于L2的流体模拟方法计算变形量，然后插值获得新的流体动画。接下来，将对比两种优化模型生成的流体表面，验证本发明的优化的准确性。两种典型的流体运动场景“液滴下落”和“溃坝”将作为实验的输入，首先通过插值结果的视觉效果来评估变形的质量，为了能更直观地观察，采用了单向最近邻插值，比如计算流体表面φ₁变换为φ₂的变形量来运行插值结果。

液滴场景的初始状态设置为两个形状与位置均不一致的液滴落入水池的过程。他们会在不同时间落入水池，因此水面出现了不一致的波纹。我们的实验结果均在插值权重为1的条件下生成，这时的源流体表面完全变形为目标流体表面，从而能够直接比较目标流体表面与插值结果的相似性。由源流体表面变形为目标流体表面的插值结果中生成流体的第二行表示由基于L2范数优化的匹配方法生成流体形态，第一行表示由基于L0范数优化的匹配方法所生成的结果。为了更直观地观察其中的差异，还可以验证两组结果的截面部分。通过结果可知，与基于L2范数优化的光流法匹配模型相比，本发明提出的优化方法生成的结果与目标流体表面更接近，改善了由未对齐表面造成的错误对应关系，因此流体表面的状态更自然。从流体表面φ₂变换为φ₁的插值结果，基于L2的变形方法得到了不准确的对应关系，使得水面上“凭空”出现了一些波纹。因此提出的优化模型能更好地对齐两个输入流体表面，得到的结果与目标流体表面更接近。另外，由基于L2方法生成的那些不自然的运动也能通过基于L0的优化模型得到缓解。

为了验证优化模型的适用性，通过“溃坝”这一场景来验证基于L0的优化模型是否能更准确地匹配流体表面。首先，将一部分水置于水池的一侧，然后“打开”大坝的开关，让水从一边流下并流过整个水池。为两组输入的流体动画设置了不同的宽度，这样生成的大坝将实现一系列的介于当前这两种水坝宽度的流体动画，即随着插值权重的增加对其应的宽度逐渐发生变化。这个实验中，设置的插值权重分别是是1和0.3。同样的，比较了基于L2范数优化的光流法变形模型与基于L0范数的优化模型所生成的结果。与目标流体表面相比，基于L0的优化模型的结果中出现的错误变形相对更少。由于基于L2正则化的模型主要针对流体表面的局部区域，而基于L0的正则化则能够考虑流体表面的整体。因此提出的优化模型能够减少在基于L2正则化的模型中产生的不自然的变形，使得结果与目标流体表面更接近。

以上的场景的流体表面是运动初期的状态，基于L0的匹配方法生成的变形量在捕获两个水花整体动作的基础上还能够恢复一些小的细节。

对于变形量质量的判断，除了以视觉效果来判断变形场的准确性以外，还能以什么定量的标准作为检验变形场质量的依据呢？有一些简单直接的方法，例如计算变形后的φ₁与目标表面φ₂之间的差值，或者比较最小化能量方程后得到的能量值的大小。但是实际上，即使变形量有让两组表面更相似的趋势，也依然被认为是不够好的。当时空表面SDF中的距离值较大时，可能对形变表面的整体形态造成影响，因此直接进行比较后很容易出现偏高的误差值。本文采用的变形场质量判断方法更多地关注了流体的表面区域，建立时空表面的误差统计模型表示如下：

error(φ₁,φ₂)＝∑_x h(φ₁(x),φ₂(x))V (20)

式中的V表示空间体积，其中指示函数h(φ₁(x),φ₂(x))用于检测表面是否对齐，具体的定义如下：

式中的s₁和s₂分别代表时空表面SDFφ₁和φ₂。每次迭代都需要用误差统计模型来判断当前的变形场与上一次的变形场相比是否使得流体表面φ₁和φ₂进一步地接近了。随着每次迭代运算的执行，利用参数κ乘以自适应参数β_A加速算法的收敛，本文在下一节中详细讨论了有关参数κ的选取问题。

表2展示了两种场景在两个变形方向计算得到的变形量分别得到的变形源流体与目标流体的误差值。基于L2的匹配方法和基于L0的匹配方法均在多分辨率下迭代计算直到所生成的流体表面与目标流体表面的误差不再降低。

表2

以上所述，仅是本发明的较佳实施例，并非对本发明作任何形式上的限制，凡是依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属本发明技术方案的保护范围。

Claims (9)

1.一种基于数据驱动的流体模拟方法，其特征在于：该方法包括以下步骤：

(1)通过两组不同的设置参数生成两个流体序列，然后转化成两组时空流体表面φ₁和φ₂作为输入数据；

(2)进行流体动态模型的配准，添加基于L0范数的光滑约束条件，通过限制变形场中所包含的非零梯度的数量实现变形场的光滑化处理；

(3)进行变形的对齐，通过残差迭代以及分辨率的分层逐渐优化；

(4)进行表面投影，为了恢复两组时空流体表面的细节，当光流匹配算法计算出表面SDF的变形场之后，通过线性查找的方式来获取更多细节的变形；

(5)进行插值运算并输出结果，利用变形场在流体表面φ₁与流体表面φ₂之间建立对应关系，将该变形场应用于插值运算生成不同流体表面间光滑且流畅的转换过程，插值运算结合两组流体表面SDF的并集，从而在结果中保留小液滴以及薄膜结构。

2.根据权利要求1所述的基于数据驱动的流体模拟方法，其特征在于：所述步骤(1)中，通过流体模拟算法获得的流体动画由一系列连续的流体表面序列组成，流体的每个微元在某一序列帧的时刻具有相应的瞬时速度信息，一个流体动画经过变形之后得到的结果仍然保持时空相关性，根据流体运动的速度场将流体动态模型转化成包含时间和空间维度的表面数据，使得在对动态数据进行变形优化时有更多的参考信息，流体模拟算法将利用不同的设置参数生成两组流体序列，然后将流体数据转化成两组时空表面φ₁和φ₂作为输入数据。

3.根据权利要求2所述的基于数据驱动的流体模拟方法，其特征在于：所述步骤(2)中，通过计算4D稠密变形场u将流体表面φ₁的点映射到流体表面φ₂，以使φ₁的表面形态变形为φ₂，通过光流技术计算两组流体动态模型之间的动作，流体模拟法的设置参数变量r作为导致两组流体表面之间出现区别的自变量，表示一种外部的控制参数，由光流的亮度恒定假设得出dφ/dr＝0，应用链式法则推出流体动态模型的配准方程表示为公式(2)：

其中，流体表面上点的位置变化dx/dr表示表面φ₁到φ₂的变形u，用表面φ₂对时间和空间的导数来近似表示

4.根据权利要求3所述的基于数据驱动的流体模拟方法，其特征在于：所述步骤(2)中，基于L0范数的光滑约束条件的光流匹配模型的能量方程为公式(9)：

其中，参数β_T用于控制吉洪诺夫正则项的作用在整个方程中所占的影响，

自适应参数β_A用于控制点p的梯度

与其对应的辅助变量v^p之间的相似性，参数β_S用于控制变形场的光滑程度。

5.根据权利要求4所述的基于数据驱动的流体模拟方法，其特征在于：所述步骤(2)中，能量方程包括两种不同性质的子问题，分别对其建立了闭合解的模型；

子问题1：v的最小化

将未涉及到变量v的项去掉，根据公式(10)得到关于v的最小化模型

如果把点p的梯度

作为已知量，对应变量v^p可以通过条件方程来进行确定，该方程的表示为公式(11)：

计算出流体表面上每个点p的最小能量，然后根据公式(12)求和就得到了子问题1的最小能量

子问题2：u的最小化，如果变量v已经确定，变形u的计算模型通过公式(13)表示：

该公式属于二次方程的最小化问题，利用梯度下降法直接求解，求E(u)关于变形量u的导数并令其为0，得到的线性方程表示为公式(14)：

其中，用两组流体表面SDF的差φ₁-φ₂来表示梯度φ_r，作为变形目标的流体表面的空间梯度表示

L为拉普拉斯算子。

6.根据权利要求5所述的基于数据驱动的流体模拟方法，其特征在于：所述步骤(3)中的残差迭代，通过光流技术求解变形场，然后对流体表面φ₁进行变形，并判断其是否与流体表面φ₂的形态更相似；如果相似那么再次执行算法来改进变形场，重复多次直到算法收敛。

7.根据权利要求6所述的基于数据驱动的流体模拟方法，其特征在于：所述步骤(3)中的分辨率的分层逐渐优化，采用金字塔结构，其中五角星代表流经某一流体微元的表面，他的位置从流体表面φ₁的左下方变形到流体表面φ₂的右上角，在红色方框标注的小窗口范围内寻找当前表面移动位置，对流体数据进行下采样得到粗分辨率的数据；在粗分辨率这一层得到的变形场记作u¹，然后传递至下一层并在此基础上继续修正变形得到u²，直至最后一层从而实现变形场由粗到精的分层迭代运算。

8.根据权利要求7所述的基于数据驱动的流体模拟方法，其特征在于：所述步骤(4)中，变形表面φ₁已经非常接近目标φ₂的情况下，沿着目标表面上某一点的法线方向找到在表面φ₁中最近的交点，然后通过缩减两者差值的方式改进变形场；当光流匹配计算的变形场使流体表面之间的差距越小，这种投影方法越稳定地恢复表面的小变形。

9.根据权利要求8所述的基于数据驱动的流体模拟方法，其特征在于：所述步骤(5)中，流体表面φ₁由变形场u_1→2变形得到φ′₁，由变形场u_2→1得到表面φ′₂，已知插值权重系数为α，插值运算的表达式为公式(15)、(16)：

φ′_1∪2＝min(φ′₁，φ′₂) (15)

φ_output＝ω₁φ′₁+ω_1∪2φ′_1∪2+ω₂φ′₂ (16)

φ_output表示插值运算生成的流体表面，插值的权重系数表示为公式(17)-(19)，其中clamp用于约束参数范围至0到1之间：

ω₁＝clamp(1-2α) (17)

ω₂＝clamp(2α-1) (18)

ω_1∪2＝1-ω₁-ω₂ (19)。

Images