段 景 耀

(上海交通大學電子信息與電氣工程學院 上海 200240)

0 引 言

流體模擬是計算機圖形學的重要分支之一，被廣泛應用在電影、游戲、動畫、虛擬現實場景以及科學工程等諸多領域?；谖锢淼牧黧w模擬過程中需要求解納維爾-斯托克斯(Navier-Stokes，NS)方程，在計算機領域NS方程的求解方式一般分為兩種：歐拉法和拉格朗日法。相比于拉格朗日法，歐拉法在高精度的模擬中具有求解更為穩(wěn)定、表面重構和渲染更為簡便等優(yōu)點，因而被廣泛使用。

在流體模擬的規(guī)模增大時，歐拉法模擬所需的計算量也迅速增加。因此，大場景中的流體模擬需要良好的加速算法來保證效率。在模擬中，泊松方程的求解是歐拉法的關鍵步驟，同時也是計算開銷最大的步驟。為了提高流體模擬的運行速度，許多研究都致力于優(yōu)化泊松方程的求解乃至整個歐拉法的求解中。

最為直接的優(yōu)化效率方法是優(yōu)化泊松方程的求解過程。例如，Bridson等[1]提出的預處理共扼梯度法，通過減少求解大型稀疏線性方程組的迭代次數，來達到加速計算的目的。McAdams等[2]在共扼梯度法的預處理過程中使用了多重網格來進行加速，在模擬較大規(guī)模的場景時具有非常良好的性能提升。

另外一種優(yōu)化效率的方法是對模擬的空間進行區(qū)域劃分，優(yōu)化計算資源的分配，在更加關鍵的區(qū)域中增加計算資源和精度。例如文獻[3-4]提出通過八叉樹數據結構對空間進行劃分。文獻[5]的方法在液體表面采用精細的立方體網格，而在液體內部使用長方體網格。文獻[6]的方法在流體的來源位置應用精細的網格，在其他區(qū)域則使用較為粗粒度的網格。文獻[7]將整個空間進行粗粒度劃分，并將每個區(qū)域進行抽象，在宏觀和局部分別求解，最后再合成流體的運動狀態(tài)，通過減小計算規(guī)模的方法減小計算量。這些方法都是在精度和性能之間進行平衡，提升性能的同時不可避免地會引入誤差。

機器學習(Machine Learning,ML)作為一個新興的學科，被廣泛應用于多個計算機領域，能夠實現自動化地在數據中提取特征進行分類或擬合的操作，具有速度快、非線性等諸多優(yōu)點。卷積神經網絡(Convolutional Neural Network,CNN)作為機器學習的重要分支之一，同樣被應用在流體模擬中。例如，文獻[8-9]利用CNN通過數據訓練來求解NS方程，來達到加速的目的。文獻[10]等人通過CNN來對流體的局部特征進行學習，再通過合成得到整個場景的流體運動狀態(tài)，也達到了良好的模擬效果。

本文提出一種基于二級網格和CNN的流體模擬方法，首先通過二級網格進行局部和宏觀的流體狀態(tài)的模擬，再通過訓練好的CNN將局部細節(jié)狀態(tài)和宏觀運動趨勢進行整合，得到最終的流體模擬結果。

相比于其他方法，本文的方法具有以下的優(yōu)點。首先，基于流體的連續(xù)性和受局部狀態(tài)影響較大的特性，可以將流體進行二級網格的劃分。這樣的劃分能夠有效減少計算規(guī)模，并允許同級網格之間進行并行計算。其次，CNN的引入能夠有效地對高低精度網格之間的非線性關系進行擬合，從而通過數據驅動的方式，減少由于二級網格分割計算帶來的誤差，在加速流體模擬的同時，仍然保持較好的數據精確性。

1 算法描述

本文的算法主要分為兩個部分，分別是二級網格的流體模擬，CNN的學習和應用。二級網格的模塊主要實現了分塊的流體模擬算法，在此階段同時生成二級網格流體的狀態(tài)數據，作為CNN學習的訓練數據。在CNN部分中，對上述數據進行學習，并將得到的訓練好的模型應用在二級網格的模擬中，得到最終的流體模擬結果。算法流程如圖1所示。

圖1 包含二級網格與CNN的流體模擬算法流程圖

1.1 Navier-Stokes方程

在歐拉方法中，流體被視為不可壓縮的。流體所在的空間被均勻分割成大小相同的網格。記流體的速度場為a，流體的密度為ρ，流體的壓強為p，f為所有的外力的合力(重力、浮力等)，則Navier-Stokes方程以如下的方式給出：

(1)

▽·u=0

(2)

式中：▽為梯度算符。若已知第n幀的速度場un,為求解最終的速度un+1，首先求解中間速度u*，滿足：

(3)

進一步，求解第n+1幀速度：

(4)

壓強p通過下列方程給出：

(5)

通過上述5個方程，便可以根據當前幀的流體狀態(tài)，求得下一幀中流體的狀態(tài)。

1.2 二級網格下的流體模擬

首先，我們將網格均勻分割為粗網格，再將每個粗網格均勻分割成為若干細網格。如圖2所示，若流體空間為正方體形狀，記每一邊上的粗網格數目為rc，每個粗網格中每一邊上細網格的數目為rf，空間中的精細網格的總數為(rc×rf)3。若此時每個精細網格的速度記為u，則我們最終要求解的是下一幀中的每個精細網格中的速度u′。

圖2 二級網格的空間劃分

(6)

平均值反映了此區(qū)域中整體的流動情況。由于流體中的速度是連續(xù)的，在較小的局部區(qū)域內流體的速度差異并不明顯，因此我們可以選取平均值作為粗網格中的速度代表值。

(7)

(8)

我們使用CNN通過大量數據的訓練來得到較為理想的函數f。

(3) 計算復雜度。在計算壓強時需要求解方程組：

(9)

式(9)可以等價為AP=b的形式，其中：A為系數矩陣，P為待求解的壓強構成的向量，b為常數矩陣。由于壓強p的長度為n3，可知矩陣A的大小為n3×n3，在實際計算中A較為稀疏，但求解A的逆矩陣A-1的復雜度必定大于O(n3)。由經驗得出，此步驟中的實際復雜度約為O(n4)。

在我們的問題中，若直接應用NS計算求解，所需時間為：

T=k·(rc·rf)4

(10)

式中：k為復雜度系數。

若進行二級網格的劃分，則在局部精細網格和全局粗網格上，NS方程求解所花費的時間總共為：

(11)

式中：t為整合粗細網格所需的時間，與rc和rf的三次方成正比?？梢钥闯龇指罘椒ǖ臅r間復雜度要小于直接求解NS方程的時間復雜度，即T′

綜上所述，二級網格劃分本身能夠降低計算量，又由于各個粗網格區(qū)域內的計算相互獨立，因此可以進行多線程或者分布式計算進行進一步加速。如將空間分割為16個粗網格，則在求解速度，可以同時計算16個網格中的速度，使此步驟的性能至少提升16倍。

1.3 卷積神經網絡

由于二級網格的區(qū)域分割的計算，將完整的流體分割成為了相互獨立的精細網格區(qū)域，這樣在提高效率的同時也不可避免地引入了誤差。卷積神經網絡(CNN)能夠較好地處理平滑曲線的擬合問題，而流體恰好具有良好的連續(xù)和平滑的特性。鑒于粗網格、細網格分別代表了局部運動和宏觀運動，但二者的對應關系未知，我們采用CNN來通過數據來擬合二者的關系。

(1) 數據生成和預處理。在二級網格的模擬過程中，我們將每一幀的粗網格和細網格中的速度數據進行記錄，作為訓練數據。同時在同樣的場景不使用二級網格，而是直接求解NS方程得到下一幀的網格中的速度，作為對照參考的真實數據。而我們的預測數據就是下一幀中細網格中的速度數據，目的是讓預測數據盡可能接近真實數據。

對于粗網格的數據，其網格大小與全局網格大小不一致，其分辨率為rc。我們采用擴充的方法，將數據集大小擴充為全局數據網格大小rc×rf，記錄每個精細網格所對應的粗網格的數據，以便于在CNN中對數據進行卷積處理。

由于流體的連續(xù)性，當時間步長較小時，相鄰兩幀之間的速度等各項流體狀態(tài)數據差異較小，因此我們在時間上進行抽樣，每N幀之中選擇1幀數據進行記錄，從而在數據集大小不變的情況下，增加了數據集的多樣性。同時，我們針對多個不同的場景進行模擬，進一步提升數據集的多樣性。

(2) CNN結構。卷積神經網絡具有這樣的特性，每個神經元都會考慮其附近區(qū)域的數據狀態(tài)，符合流體的局部區(qū)域影響的特點，與我們二級網格劃分的思想高度相似。鑒于我們的輸入數據(當前幀的粗網格和細網格的速度場)和輸出數據(下一幀的細網格的速度場)的形狀相同，我們采用兩層帶有邊界補齊padding的卷積層對數據進行學習。

圖3 卷積神經網絡結構

在此神經網絡結構中，要訓練的對象是每一層中的卷積核conv。例如在我們的實際應用中，每個卷積核寬度為7，深度為6，即每一個卷積層中卷積核中的數據個數為7×7×6=294?？偣残栌柧毜木矸e核變量的數目為588個。通過對大量訓練數據進行學習，使預測的誤差不斷減小，從而找到預期的卷積核的值。

(3) Loss函數。Loss函數作為衡量預測結果的標準，直接影響到模型的準確性。Loss函數的選取應該具有以下幾點特點：(1) 模型預測越準確，Loss函數的返回值越小，反之，Loss函數的返回值越大;(2) Loss函數的選取應盡量簡潔，使其自身的求解和梯度下降法反向傳播梯度的時候，不會帶來過大的運算量，保證機器學習的效率。

基于這兩點要求，我們采用每個維度上，預測數據與真實數據的L1距離差異比例作為模型的Loss，即：

(11)

在實際操作中，由于我們的模型是遞歸進行的，即上一幀的輸出作為下一幀的輸入，倘若在神經網絡的計算過程中，流體速度的比例產生縮放，即計算得到的u′與真實的utruth的關系為u=αutruth，則隨著模擬的進行，u′與真實的utruth之間的差異將會呈指數級增長。如果α小于1，u′將會迅速趨近于0;如果α小于1，u′將會迅速增大，以至于模擬場景“爆炸”。因此我們增加一項Loss，來衡量u′與真實的utruth的比例關系，即:

Lossscale=w∑(convweight-1)

(13)

式中：convweight為卷積核的權重之和，w為Lossscale所占總的Loss的權重。我們將此項與式(12)中的Loss相加，得到最終模型的Loss函數。

2 實驗結果

2.1 實驗效果分析

首先在數據上，我們通過Loss函數對我們方法進行評估。為了避免過擬合，我們對數據集訓練的Epoch數目為2 000。通過訓練，我們的方法訓練的最終Loss值約為0.34(見圖4)。此數值在具體流體場景中的意義為：在對速度場求解結果的逐網格的比較中，我們的方法與PCG相比的逐網格的差異平均為34%左右。作為對比，在未使用CNN而采用線性插值時，根據場景的不同，此數值在0.6～2.1之間變動，即在最壞的情況下相對于PCG誤差達到了209%左右。

圖4 CNN訓練EPOCH與Loss的關系

如圖5、圖6所示，在視覺效果上，我們的模型在保持流體宏觀走勢的同時仍然可以展示流體的微觀細節(jié)。

(a) 本文方法(b) 直接求解NS方程圖5 模擬結果

(a) 本文方法 (b) 直接求解NS方程圖6 模擬結果

2.2 性能分析

我們使用的機器的處理器為Intel i7-3770@CPU 3.40 GHz×8，內存為15.6 GB,顯卡為NVIDIA GTX 970。

在實際操作中，取N=5即每5幀取一幀數據，粗網格和細網格的分辨率均為16，即場景中的分辨率為256，同時在450個不同的場景中進行模擬，每個場景模擬100幀。在模擬過程中，總計耗時約為1.3 h，我們總共得到了450×100/5=9 000組數據，數據總大小約為1.8 GB。

訓練過程中，我們對每個訓練集訓練的Epoch數目為3 000，批處理的大小為10組數據，訓練時間長為3.1 h。在保持流體形態(tài)的同時，我們的方法能在分辨率較高時夠大幅縮減模擬的時間，提高幀率。每幀所需的時間與分辨率的關系如表1所示。在網格規(guī)模為1 024×1 024時，我們的方法的運算效率僅為PCG方法的一半左右。

表1 兩種方法每幀耗時和分辨率的關系

3 結 語

本文提出了一種基于二級網格和機器學習的流體模擬方法，通過二級網格的區(qū)域劃分來加速流體狀態(tài)的求解，并通過卷積神經網絡對劃分的區(qū)域進行整合，同時通過數據驅動的方法來消除誤差，從而在加速流體計算的同時基本上保持原本的流體形態(tài)特征，達到了良好的模擬效果。

二級網格的使用建立在流體連續(xù)性較好，流體的運動狀態(tài)受局部影響較大的前提下。但當每個粗網格包含的細網格數目(rf)較大時，二級網格的區(qū)域劃分將會對空間造成較大的離散程度，由此產生較大的誤差。

二級網格劃分減少了計算量，因而減少了整個流體場景中所包含的信息量，我們的做法本質上是利用CNN，通過數據驅動的方法外推因計算量減小而丟失的信息，但當場景本身網格數量太少時，其包含的信息量便不足以支持外推，此時本方法將會失效。此外在網格分辨率較小時，PCG方法本身的耗時較小，但本文方法由于需要處理粗細網格間的調度，同時還需調用CNN進行預測，因此造成在低分辨率下運算效率低于PCG方法。因此本文方法更適用于大規(guī)模的流體場景。

此外，由于使用了數據驅動的方法，本文算法結果依賴于大量的真實實驗數據。數據的數量不足或者多樣性不足，會使得此方法在特殊邊界條件或一些未訓練過的場景中失效。由于我們的訓練模型依賴于流體模型，當流體模型發(fā)生改變時(如流體的空間結構、大小和劃分方式等)，訓練結果仍然會失效，因此我們模型運行的前提條件是流體模型的一致性。

本文的工作提供了很多可擴展性。首先，每個劃分的區(qū)域都可以進行并行或分布式計算。其次，可以進行多于二級的網格劃分，并使用同樣的方法進行訓練，進一步優(yōu)化運算性能并保持流體的形態(tài)。除此之外，可以在CNN模型中增加對其他流體變量的分析學習，如對流體壓強和密度的分析，從而進一步增加整合后流體的精確程度。