基于波環粒子的實時水波仿真方法

顧浩杰，張 軍，2*

（1.江南大學 人工智能與計算機學院，江蘇 無錫 214122；2.江蘇省媒體設計與軟件技術重點實驗室，江蘇 無錫 214122）

0 引言

流體模擬一直是計算機圖形學中極具挑戰的研究問題，其計算復雜度遠超現有通用計算機的實時計算能力。與此同時，由于自然場景的虛擬仿真中不可避免涉及流體對象，因此流體模擬吸引著包括互動游戲，媒體藝術以及增強現實（Augmented Reality，AR）與虛擬現實（Virtual Reality，VR）等領域的大量開發人員。

水波模擬是流體模擬研究中一個較為特殊的研究點，直接采用流體模擬通用方法求解納維-斯托克斯方程（Navier-Stokes equations，NS），很難實現高逼真度的實時計算性能。伴隨著游戲領域出現大量戶外開放環境仿真的需求，水波的實時仿真和渲染技術越來越受到行業人員的重視。

自1980 年以來，計算機圖形學領域的學者已經提出過多種模擬水波方法。這些方法的主要策略是采用微幅波假設［1］對水波問題進行線性化處理，以近似求解滿足不可壓縮條件的NS。根據具體模擬方法的不同，可以將這些方法細分為歐拉網格法［2-3］和拉格朗日粒子法［4-5］。

歐拉網格法計算開銷相對較少，但會受CFL（Courant-Friedrichs-Lewy）條件限制。該方法將空間分割為離散的元胞，然后追蹤流過元胞的被模擬材料，并將振幅等模擬信息記錄在每個網格單元里。2015 年，Jeschke 等［6］提出一種高效的波前追蹤算法模擬水波與復雜環境的交互效果，并利用多值函數插值的方法避開了奈奎斯特限制。這種模擬水波的方法需要進行大量的預計算，因此不適用于實時交互的可變場景。為了模擬出真實的頻散效果，Canabal 等［7］在求解偏微分方程數值解的基礎上使用頻散核，表現了豐富的水波動態細節，并通過金字塔過濾器和陰影卷積運算有效提高了計算效率以及與障礙物碰撞現象的幾何計算效率。Jeschke等［8］基于小波的水波離散化，使水面波的模擬既有豐富的水波細節，又可以高效處理碰撞關系。

相較于歐拉網格法，拉格朗日粒子法避免了網格分辨率的限制，可以在任意不規則區域實現水波仿真，并避免了無水波區域的“空計算”問題。拉格朗日方法將水波分割為離散、獨立的元素，然后分別追蹤其中每個元素的運動過程。每個離散的元素在內存空間中有一個獨立的數據結構體，負責攜帶相關水波的模擬數據。2007 年，Yuksel 等［9-10］通過拉格朗日粒子法模擬水波的反射效果，并根據體積守恒的規則推算出物體入水與出水過程中水波在物體周圍的位置與振幅。該方法渲染出的水波具有一定視覺逼真性，計算性能也達到了實時交互水平，是基于粒子系統實時模擬水波交互場景的經典算法。

Yuksel 等［9-10］的方法在水波物理特性的細節（比如頻散現象）方面仿真精度非常有限，經不起仔細觀察。對此，Jeschke 等［11］引入“波包”的概念，不僅能有效模擬出水波的高頻波細節，而且使水波具有物體上更加精確的波速。最近，Skrivan 等［12］基于波包的概念，提出水波曲線方法，通過插值函數與多個控制點代表曲線的方式，成功地在三維動態流體表面上模擬出豐富的水波細節。雖然這種方法模擬出的水波更為細致，但其計算復雜度較大，無法達到實時繪制的要求。

拉格朗日粒子方法面臨的主要問題是粒子數量在水波發展初期和后期的差異非常大，導致計算量不穩定。初始水波較為簡單，少量粒子即可完成仿真需求；而后期的水波環半徑急劇增大，且發生復雜碰撞后會產生更多新水波，必須使用大量粒子才能疊加產生視覺上無裂痕的仿真效果。隨著大量新粒子被添加到仿真系統，計算量會顯著提升，渲染幀率出現大幅下降，嚴重影響開發人員對幀率的控制能力。

本文將有效減少水波模擬中粒子數量作為解決上述問題的途徑。而減少粒子數量的一個直接方法就是增大粒子的覆蓋面積，以避免海量微小粒子堆積在一起（如圖1 所示）。為此，本文提出以波環作為基礎粒子的方法（如圖2（b）所示），動態提高粒子覆蓋面積，從而顯著減少水波模擬過程中所需的粒子數量。為解決由此帶來的水波反射時碰撞分析的復雜性問題，本文采用鏡像波源的方法。該方法有效解決了波環粒子作為圓環無法模擬碰撞反射的問題，并同時解決了傳統小粒子處理反射可能出現的水波斷裂問題。實踐表明，本文方法在保持與波包方法近似的逼真度的同時，大幅減少了粒子使用數量，并有效控制了粒子增長規模，能在低端PC 機、移動平臺等低性能硬件上實現交互級水波仿真應用。

圖1 微小粒子堆積現象Fig.1 Accumulation phenomenon of small particles

1 算法理論

由于不同頻率、不同方向傳播的水波彼此不相干，所以水面可認為是大小不等、方向與頻率各異的水波的疊加組合——波包（如式（1））。相較于傳統的采用鄰域傳播思想對水波進行動態造型的方法［13］，使用拉格朗日粒子法并引入波包的概念，往往可以使模擬的水波更接近真實的水波。波包是局限在空間的某有限范圍區域內的波動，波包整體按群速度cg在空間移動，內部按相速度cp移動。

在小粒子上引用波包，需要在二維條件下進行水面高度的計算（如圖2（a）的矩形塊）。本文利用波環粒子表示波包的圓環結構（如圖2（b）），有效地將水面高度的計算問題維持在一維條件下，能大幅度減少模擬時所需的粒子數量，降低計算開銷。

如圖2（a）所示，代表相同大小的水波時，所需矩形粒子的數量遠多于所需波環粒子的數量。隨著水波繼續發展，其波環半徑持續增大，周長也相應增加。為保持波包與波包相互重疊（按照香農采樣定理，必須重疊一半才能保持視覺連續性），新的波包必須持續加入到場景中，大幅延長了計算時間；而本文的波環只需修正其內部參數，波環數量始終保持不變，計算量保持穩定。

圖2 矩形粒子與波環粒子模擬水波對比Fig.2 Comparison of rectangular particle and wave annulus particle in simulating water waves

依據波環粒子作為基礎粒子的方式，本文算法結構如圖3 所示。

圖3 本文算法結構Fig.3 Framework of the proposed algorithm

首先，水面網格點的高度計算可以簡化成網格點上波環粒子高度的疊加，因此高度函數可以近似表達為：

其中：N是網格點上對應的波環粒子數量，aj是波環粒子的振幅，kj是波環粒子的代表波數，cg是波環粒子的半徑擴張速度，cp是波環粒子的內部相速度。?(x)代表核函數，因為波包需要一個核函數來描述其在空間中起局部作用，并且其傅里葉變換Φ(k)在波數坐標系中起局部作用。Φ(k)越寬，波包則越狹窄，一般使用高斯函數作為波包的核函數。此外，波環粒子的寬度l選擇為3 倍波環粒子的代表波長λj=2π/kj。

其次，由于波環繼承了波包的概念，水波的頻散變得易于模擬。通常使用ω(k)來描述頻散關系，相速度由此表示為ω/k，群速度則表示為?ω/?k。角頻率ω與波數k的關系一般表示為：

其中：g是重力，σ是表面張力，ρ是水的密度，k是波數，h是水深。由于波環粒子無法直接反映各處的水深帶來的影響，因此需利用波環上多點的水深取平均的方式來近似水深對波環的影響。

波的相速度cp是波的相位在空間中傳遞的速度，描述波包中某一單頻波的相位移動速度，其計算公式為：

式（3）正確描述了當水深h很大（深水波）或者波數k很大（表面張力波）時，相速度與水深無關；其他情況，相速度受水深影響。

群速度cg是能量或信息順著波動傳播的速度，其計算公式為：

根據式（4），當不同波長的水波表現出不同的傳播速度時，造成混合各種波長的波環漸漸地在空間中分離開來，從而實現頻散的效果。

式（5）即為最終水面高度的計算公式。此外，為了便于計算，本文將局部坐標范圍控制在[-1，1]。通過對不同頻段的水波進行疊加，可以得到較為真實的水波，如圖4、5 分別展示了多頻段一維波函數疊加的效果與三維下多頻段波環疊加的效果。

圖4 多頻段一維波函數疊加Fig.4 Multiband one-dimensional wave function stacking

圖5 多頻段波環疊加Fig.5 Multiband band wave annulus particles stacking

最后，在傳統的小粒子模擬水波的方法中，面對水波的碰撞反射，往往讓入射水波和反射水波遵守光學中的鏡面反射定律［14］，并由此更新碰撞后的粒子行進方向與位置。相較于小粒子，本文的波環粒子是代表一個圓環上的水波，而非水波中的局部，所以無法直接通過更新傳播方向進行處理。綜合考慮性能和精確性后，本文采取鏡像虛擬波環的形式處理水波反射現象。當發生波環與容器壁發生碰撞后，以碰撞邊界為鏡面，在鏡像波源處產生相同大小的新波環粒子，來模擬碰撞后的反射效果。此外，由于頻散模擬會產生大量的同波源的波環，因此波環的數量遠遠大于波源的數量。使用波環檢測是否與障礙物邊界發生碰撞，需要大量成本來計算距離dki(k=1，2，…，p，i=1，2，…，n)，其中k為波環數量，n為邊界數量。為了減少計算成本，本文采用障礙物邊界檢測碰撞的方法。如圖6 所示，當且僅當出現新波源時，需更新各邊界的動態表；邊界依據波環粒子的波源ID 查找應該發生碰撞的距離，并與波環粒子當前時刻的半徑進行對比以檢測碰撞，從而達到避免重復計算的目的。

圖6 邊界動態表的更新Fig.6 Update of boundary dynamic table

2 水波環粒子的具體實現

本文算法的計算流程如圖7 所示，在CPU 上處理水波的行為狀態，在GPU 上處理水面網格的高度計算以及最終的水面繪制。

圖7 本文算法具體實現流程Fig.7 Implementations flow of the proposed algorithm

水波行為狀態處理需要每一“時步”調整波環粒子的屬性，以保證水波物理特性以及位置的正確。水面高度的計算則需要水面網格點根據不同波環粒子的波源、外環的半徑以及波環寬度來判斷自身在哪些波環內，并疊加高度得到最終改點高度，如圖8 所示，其中，Oi是波源，ri是波環粒子外環半徑，li是波環寬度。

圖8 水面網格的高度計算Fig.8 Height calculation of water surface grid

2.1 水波頻散效果模擬

因為水波物理屬性的復雜性以及可視性，水波頻散效果對于水波的模擬十分重要。水波波長的不同，會使其傳播速度也不相同，因此同一波源所產生的各種不同頻率波在向外傳播時會發生分離的現象。為了使水波具有較為真實的頻散效果，本文采用頻散細分方法，為每個波環粒子分配一個具體的波數范圍[kmin，kmax]，使波環粒子的代表波數kj=(kmax-kmin)/2。根據式（4）可知，群速度是關于波數k的函數，群中的一些波比其他的波傳播得更快，會將數據包的一部分推到平均群速度之前；而較慢的波則會將數據包的另一部分拖到后面。因此當分別以cg(kmin)與cg(kmax)行進的虛擬距離相差大于0.3 倍的波環寬度時，生成一個同波源的新波環。兩個波環的波數范圍分別調整為[kmin，kj]與[kj，kmax]，調整后兩個波環的代表波數不同，傳播速度也隨之發生改變，從而模擬出頻散效果。

此外，當發生頻散后波環粒子細分為兩個新波環粒子的情況，遵循能量守恒定律，需要對波環粒子的振幅與能量進行重新計算。振幅與能量滿足：

其中：E為波環粒子能量，a為波環粒子振幅，S為波環面積。

假設發生頻散后，兩個新波環陡度是一樣，則可以得到兩個新波環粒子的能量之比：

其中：Ei是對應新波環粒子的能量；li是對應新波環粒子的波環寬度；ki是對應新波環粒子的代表波數；r是波環粒子的外環半徑，因為細分是對波數的劃分，因此兩個波環粒子的外環半徑相同。

圖9、10 分別展示了有無頻散對水波模擬的影響。通過對比可以明顯看出，圖10 中不同波長的水波是以不同的傳播速度行進的，這使模擬的水波更為真實。

圖9 沒有頻散現象的單波源模擬Fig.9 Simulation of single wave source without dispersion

圖10 具有頻散現象的單波源模擬Fig.10 Simulation of single wave source with dispersion

2.2 水波環反射計算

在復雜環境中，水波與障礙物發生碰撞會發生反射。當發生碰撞時，為了體現障礙物的表面粗糙性，可以將碰撞部分的水波的頻譜拉伸到高頻。針對不同計算平臺的性能差異，本文通過設定全局碰撞變量γ控制波環的反射次數。γ值越大，反射次數越多，但計算量也會隨之增加。由于本文為減少重復計算，記錄各障礙物邊界到各波源距離，所以邊界檢測波環粒子時，利用波環上一“時步”的外環半徑與當前“時步”的外環半徑進行處理。碰撞處理的算法如下所示。

為驗證頻散導致波環粒子增多而不會增加波源數目，表1 對不同時間段兩處水波產生至不斷擴散與頻散的過程進行部分記錄，并由表2 給出所使用的邊界動態表的并集。

表1 波環粒子的半徑變化Tab.1 Radius variation of wave annulus particles

表2 邊界動態表的并集Tab.2 Union of boundary dynamic tables

圖11 展示了本文方法所模擬的碰撞效果，此外，也可以清楚看到不同的γ造成的幀數與效果上的差別，相較于γ=1，當γ=2 時，有明顯的二次反射，但粒子數量也明顯增多，導致幀數下降。

圖11 本文方法的碰撞效果模擬Fig.11 Collision simulation of the proposed method

2.3 水波擴散效果模擬

水波向外擴散的過程中，波環粒子的面積會發生改變，根據能量守恒定律，波環的振幅需要隨之改變；所以，每個“時步”需要使用式（6）重新計算每個波環粒子的振幅。

此外，為了增加水波模擬的真實性，本文在對水波擴散的模擬中，對波環能量引入阻尼項。雖然現實中存在各種因素會使水波損失能量，但很多因素難以在較小的計算成本下模擬出；因此，本文主要考慮粘性勢流、表面污染以及水波破碎所引起的能量耗散。粘性勢流與表面污染可以通過對能量的衰減進行表示［15-16］，公式為：

其中：ν是水的粘度，一般約為10-6m2/s；a是波環的振幅；ξ≥0 為阻尼系數，用于調整數值耗散的程度。

水波的穩定性的一個重要參考指標是波陡，當波陡達一定限度，波峰就會破碎。由于無法直接模擬出水波的破碎的效果，因此通過對波陡大于0.142［17］的水波移除其一部分能量的方式進行近似的模擬。

為驗證水波擴散過程中阻尼帶來的影響，通過追蹤單個波環粒子的振幅變化，從表3 可以看出振幅衰減速度的不同。為體現在視覺效果上的差別，圖12 和圖13 分別展示無阻尼和使用本文的阻尼項所模擬的水波擴散效果。通過對比可以看到，增加阻尼后，有效讓振幅得到了漸進的指數衰減，并且讓短波相較于長波得到較快衰減，使水波在擴散的過程中更為真實。

表3 不同阻尼下單粒子的振幅對比Tab.3 Comparison of amplitudes of single particles with different damping

圖12 無阻尼的水波擴散效果Fig.12 Wave diffusion without damping

圖13 本文方法的水波擴散效果Fig.13 Wave diffusion of proposed method

3 實驗與結果分析

本文實驗平臺為普通筆記本計算機，主要硬件配置為Intel Core i7-6700HQ，顯卡為GeForce GTX 1060M。

為體現水波模擬的效果與實時性，本文將實驗結果與主流的文獻［11］方法進行橫向對比，其中本文的全局碰撞變量γ設定為3，阻尼系數ξ設定為1。

通過對粒子邊界的繪制（如圖14）可以看出，在模擬相同大小的水波時，本文所需的波環粒子數量遠遠小于文獻［11］所需的矩形粒子數量。

圖14 水波模擬所需粒子數對比Fig.14 Comparison of number of particles required for water wave simulation

圖15 展示了本文方法和文獻［11］方法水波碰撞反射后的擴散效果。可以從圖中長方形框中看到，兩種方法都能產生較為清晰的二次反射。此外，由于文獻［11］中采用矩形粒子，所以在發生碰撞反射后，改變方向后已細分的矩形粒子無法快速再次左右細分，會產生如圖15（b）橢圓形框中水波斷裂的現象；而本文是以鏡像波源產生新波環粒子的方式處理碰撞反射，因此避免了水波的斷裂。

圖15 水波碰撞反射效果對比Fig.15 Comparison of water wave collision reflection effect

圖16 展示了文獻［18］、本文以及文獻［11］對雨滴漣漪的模擬效果，可以觀察到，本文方法與文獻［11］方法相較于文獻［18］方法具有更為真實的頻散特性。

圖16 雨滴漣漪模擬效果對比Fig.16 Comparison of raindrop ripple simulation effect

為了對比本文方法和文獻［11］方法在更新粒子數據上的成本，本文追蹤在復雜環境下單波源水波產生到消逝的過程，記錄其在CPU 端處理水波行為以及各種狀態更新所需的計算成本。由于模擬計算時波包和波環粒子數量是動態變化的，造成CPU 端消耗時間波動較大，對消耗時間采用滑動平均，再取對數的方式進行對比分析。如圖17 所示，可以明顯看出，相較于文獻［11］方法，本文方法在這方面的計算成本一直維持在2 ms 以下，對粒子數據更新的開銷低且穩定。

圖17 單波源水波CPU 端計算時間對比Fig.17 Comparison of calculation time of single wave source on CPU

在上述追蹤在復雜環境下單波源水波產生到消逝的過程中，還同時記錄了模擬時的幀率差別，圖18 展示了本文方法和文獻［11］方法模擬時的幀率。可以發現，本文方法幀數始終高于文獻［11］方法，且維持在60 幀以上。

圖18 復雜環境下單波模擬幀率對比Fig.18 Comparison of single wave simulation FPS in complex environment

為了進一步定量分析本文方法的實時性，表4 給出在復雜環境下同時產生多波源水波的情況進行追蹤的時間統計結果。通過對比可以發現，本文方法在同時有多個波源水波產生至消逝的過程中，能保持流暢穩定，而文獻［11］方法則會在中間時段出現明顯的卡頓現象。

表4 復雜環境下多波源模擬幀率對比Tab.4 Comparison of multi-wave source simulation FPS in complex environment

最后，圖19 給出本文方法實際水池場景中的具體應用效果。

圖19 水池場景中多波源水波模擬Fig.19 Simulation of multi-wave source water wave in pool scene

4 結語

本文利用水波動態過程可以近似為中心平靜且呈環形擴展的特點，將波環作為基礎粒子用于水波實時仿真模擬。該方法繼承了波包算法的物體特性逼真優點，并顯著降低了水波模擬過程中對粒子數據更新的計算成本。通過使用鏡像波源的方法，新波環算法可以正確模擬出水波碰撞后的反射波，有效避免復雜碰撞檢測步驟。通過大量對比實驗表明，本文方法的渲染結果在保持足夠水波動態細節的情況下，使用更少的粒子數量，且不會出現波紋斷裂現象。最終，本文提出的方法在普通PC 上可以達到更高的渲染幀數，從而更加適合移動、廉價硬件平臺部署虛擬仿真應用。

由于波環在碰撞邊界后會繼續向外傳播，故本文方法不能對水面漂浮的小型障礙物實現水波的阻擋效果，這是算法今后的改進方向。