艾祖亮， 張立民， 于文龍

（海軍航空工程學(xué)院，山東 煙臺 264001）

多重散射的天空光照效果建模與實時繪制

艾祖亮， 張立民， 于文龍

（海軍航空工程學(xué)院，山東 煙臺 264001）

依據(jù)大氣散射的物理原理，提出了一種考慮多重散射的天空光照效果建模與實時繪制方法。該方法首先以太陽和天空光作為光源建立了多重散射的天空光照效果模型，然后綜合多種大氣粒子密度, 采用合理的分段采樣策略，對天空顏色模型的積分進行簡化，以減少積分近似計算所帶來的誤差；通過對簡化后的模型進行分析提出了采用二維紋理與三維紋理對光學(xué)深度預(yù)計算的方法，避免了運行時計算光學(xué)深度積分的問題；最后該算法在GPU的片段處理器上執(zhí)行，實現(xiàn)了天空光照效果模型的實時繪制，可以滿足各種實時應(yīng)用需求。

多重散射；天空光照；光學(xué)深度；預(yù)計算；實時繪制

太陽光在大氣中的多重散射對于天空光照效果的逼真模擬具有重要的意義。大多數(shù)大氣散射模型[1-5]在模擬天空光照效果時只考慮單散射，而不考慮多重散射光。由于單散射模型忽略了天空體本身作為天空光光源的影響，并不適合模擬黎明或黃昏時的天空光照效果，這主要是因為此時太陽位于水平線下方，僅僅離太陽方位近的區(qū)域才能夠收到直接的太陽光，而大氣層其他區(qū)域只能收到非直接的多重散射的天空光。

因此，許多學(xué)者提出了多重散射模擬方法，Preetham等[6]通過一個分析模型去擬合多散射蒙特卡羅（Monte-Carlo）仿真的結(jié)果，但是他們的模型要求視點必須在地面上才有效；Nishita等[7]和Haber等[8]用體積輻射度方法計算多重散射，但是他們的方法很難達(dá)到實時性；Bruneton和Neyret[9]提出了一個實時大氣散射繪制方法，并考慮了多重散射，該方法顯著提高了多重散射模型繪制速度，但在實時性上距離工程實際應(yīng)用仍然存在一定的距離。

國內(nèi)，劉世光和彭群生[10]及錢徽和朱淼良[11]分別對大氣效果的建模和繪制技術(shù)進行了綜述；王長波等[12]提出一種考慮大氣折射的天空光照模型, 綜合考慮大氣壓強、溫度及水汽壓等因素,實現(xiàn)了日、月升落中形狀及其周圍光暈色彩的變化等效果，但是他們的方法非常復(fù)雜，不能達(dá)到實時性；黃雷等[13]提出一種適用于雨霧天氣的天空光照的多粒子散射模型，實現(xiàn)了雨霧場景中大氣散射效果及雨后彩虹的真實感繪制，但是該模型僅僅只考慮了單散射；杜芳等[14]和馮玉康等[15]分別基于GPU對Nishita等[16]的大氣單散射模型進行了實現(xiàn)，但對于多重散射并沒涉及。

本文提出了一種考慮多重散射的天空光照效果建模與實時繪制方法。該方法首先以太陽和天空光作為光源建立了多重散射的天空光照效果模型，然后根據(jù)多種大氣粒子密度采用合理的分段采樣策略對天空顏色模型的積分進行簡化；通過對簡化后的模型進行分析提出了采用二維紋理與三維紋理對光學(xué)深度預(yù)計算的方法，避免了運行時計算光學(xué)深度積分的問題；最后該算法在 GPU的片斷處理器上執(zhí)行，實現(xiàn)了多重散射的天空光照效果的實時繪制，可以滿足各種實時應(yīng)用需求。

1 天空光照效果模型

1.1 太陽光與天空光模型

多重散射模型以太陽光和天空光為光源，通過借鑒Nishita模型的思想[7]對天空光進行建模，以大氣中任意一點p為中心，按照方位角 φj和高度角ψj采樣把整個天空體分割成若干個發(fā)光帶，其中，單位為度。發(fā)光帶的強度由采樣線的強度確定，為了簡化計算，模型不考慮地面方向的多重散射光，如圖1所示。

圖1 太陽光與天空光

入射到大氣中任意一點p的光強是太陽光和天空體各發(fā)光帶入射到該點處光強的總和Ip(λ)，公式如下：

其中，λ為光線波長，天空體發(fā)光帶數(shù)目為 Nsky，θj為天空光發(fā)光帶采樣線入射方向與太陽光入射方向的夾角，也稱散射角。

首先計算太陽光 Isun,p(λ)，s為太陽光經(jīng)過點p與天空體圓頂?shù)慕稽c，Is(λ)為太陽光在進入天空體圓頂前的初始入射光強，根據(jù)大氣光衰減模型[17]可知：

式(2)中， t(s p,λ) 為光學(xué)深度，定義為大氣粒子的消光系數(shù) βex（λ） 與密度比ρ的乘積沿著光線傳輸路徑的積分[16]，公式如下：

其中， H0為大氣粒子均值高度，不同類型大氣粒子的消光系數(shù)與密度比是不同的，對于空氣分子主要發(fā)生Rayleigh散射，消光系數(shù) β(λ)等于散射系數(shù)λ)；對于其他大氣懸浮粒子主要發(fā)生Mie散射，消光系數(shù)(λ)等于粒子的散射系數(shù)β(λ) 與吸收系數(shù) β之和。公式(3)中的βex(λ)ρ實質(zhì)上是等于大氣中所有存在的粒子消光系數(shù)與密度比乘積的和。

然后計算天空光發(fā)光帶 Isky,p(λ,φj,ψj,θj)， pj為發(fā)光帶采樣線 ajp上任意一點， sj為經(jīng)過 pj的太陽光與天空體圓頂?shù)慕稽c，則點 pj處的光強為：

其中， F(θj,λ)為散射相位函數(shù)[18]，描述了入射光線在某一點散射后進入特定散射方向的光能占整體入射光能的比例，它取決于發(fā)生散射時的光線波長、粒子尺寸和散射角。

光強 Ipj(λ,θj)經(jīng)過衰減到達(dá)p處的光強為：

天空光發(fā)光帶Isky,p(λ,φj, ψj,θj)為 光 強Ipjp(λ,θj)沿著采樣線 ajp的積分：

1.2 天空顏色模型

天空光照計算是天空顏色（sky color）的計算，主要針對天空體圓頂（sky dome）的光照效果模型，如圖2所示。

圖2 天空顏色模型

圖2中， pv為視點， pa為視線與天空體圓頂?shù)慕稽c，θ為視線方向與太陽光入射方向的夾角，θsky,j為視線方向與天空光發(fā)光帶采樣線入射方向的夾角，θj為天空光發(fā)光帶采樣線入射方向與太陽光入射方向的夾角。點p處的太陽光與天空光總的入射光強 Ip（λ）可以通過上一節(jié)公式得出，則光強 Ip（λ）在點p處發(fā)生散射進入視線方向的光強為：

視線方向總的入射光強為 Ippv(λ,θ,θsky,j)沿著視線 pvpa的積分：

上式即為天空顏色模型，通過該公式可以計算出大氣中任意一點所有視線方向的天空顏色。

2 天空光照效果的實時繪制

2.1 模型的積分簡化

天空顏色模型 Iv(λ,θ,θsky,j)的計算涉及太陽光和天空光沿著視線方向的積分，天空光發(fā)光帶Isky,p(λ,φj,ψj,θj)也需要計算太陽光沿著發(fā)光帶采樣線方向的積分，而上述積分并無解析表達(dá)且計算非常耗時，因此需對該式進行簡化。

積分的近似計算通常采用對積分路徑分段采樣，然后通過矩形、梯形或者拋物線法來求得積分結(jié)果，積分路徑的分段間隔往往是相同的。但是對于大氣散射模型來說，如果采用常量間隔的分段采樣方式，積分近似的誤差會非常大，這主要是因為大氣粒子密度是隨著海拔高度以指數(shù)方式遞減的，理想的情況應(yīng)該是分段間隔與大氣粒子密度成反比，小的間隔對應(yīng)于低的高度，大的間隔對應(yīng)于高的高度。

為此，按照海拔高度把大氣層離散為具有一定厚度的一層層球殼和 Hi,max分別為每層球殼的內(nèi)表面與外表面高度，則有）。通過這些球殼對視線方向上的積分路徑進行分段采樣，每一段內(nèi)的粒子密度采樣值取該段內(nèi)中間點的密度值，每一層球殼的厚度就等于分段間隔，如圖3所示。

圖3 大氣層球殼

大氣層的厚度大約在1000km以上，并沒有明顯的界限，但在飛行模擬視景仿真中，通常只考慮離地球海拔高度為35km以下的大氣層，即假設(shè)大氣層最大的厚度 Hmax= HN,max= 35km 。之所以這樣假設(shè)是因為在此高度之上的大氣粒子已經(jīng)非常稀薄，對大氣散射光照效果模型計算結(jié)果的影響完全可以忽略。

積分路徑的分段大小要綜合地考慮大氣粒子密度的分布情況，使其在段中近似恒定，分段數(shù)N可根據(jù)繪制的精度來調(diào)整。為了得到合理的分段間隔，假設(shè)每層球殼中的粒子密度衰減都相同[8]：

由于不同粒子的均值高度 H0不同，空氣粒子均值高度通常為7.994km，灰塵粒子的均值高度為1.2km，其他霧粒子、雪粒子、雨粒子的均值高度與具體氣象條件有關(guān)，為了使得分段間隔更加合理，公式(10)實際采用的均值高度為上述所有存在粒子均值高度的平均值 Ha。

求解公式(10)得到如下：

其中， H1,min= 0, Hmax= 35km , i= 1,… ,N -1。

通過該遞推公式可以計算出所有球殼內(nèi)外表面高度以及球殼厚度，然后用這些球殼對天空顏色模型和天空光模型的積分路徑分段采樣，每一段內(nèi)的采樣點取該段的中心點，顯然每段中心點的高度即是相應(yīng)球殼的中心高度 (Hi,min+ Hi,max)/2，每一段內(nèi)的粒子密度就是該中心高度的密度。通過上述分段采樣方法，公式(7)與公式(9)可以簡化為：

其中，Nsky,j表示天空光發(fā)光帶采樣線方向積分路徑分段采樣數(shù)目， Nv表示天空顏色模型視線方向積分路徑分段采樣數(shù)目， pj,k與 pi分別為相應(yīng)的采樣點。

2.2 模型的預(yù)計算與繪制

經(jīng)過上節(jié)的簡化解決了天空顏色模型與天空光發(fā)光帶的視線方向積分計算問題，但是光學(xué)深度積分的計算仍然不可避免，O'Neil[19]針對光學(xué)深度的計算提出了預(yù)計算的方法，用一個二維紋理存儲光學(xué)深度的計算結(jié)果，這張表將光學(xué)深度的積分式替換成一張二維紋理查找表，表的兩個維度分別是視點高度與視線方向的天頂角。二維紋理查找表方法利用太陽光是平行光這一近似假設(shè)，通過查找表就可以獲得大氣中任何一點到大氣層頂部的光學(xué)深度，避免了運行時計算光學(xué)深度的問題。采用二維紋理查找表對于單散射天空光照效果模型的計算取得了實時繪制效果，但是對于多重散射天空光照效果模型，由于光學(xué)深度的計算量巨大，只采用二維紋理查找表仍然難以達(dá)到實時繪制要求。因此，針對這個問題，本文提出了一種采用二維紋理和三維紋理對光學(xué)深度預(yù)計算的方法，并基于 GPU實現(xiàn)了多重散射天空光照效果模型的實時繪制。

下面對天空光照顏色模型公式(13)的計算進行分析，公式第一項任意采樣點 pi來自太陽光的光強貢獻為：

公式第二項任意采樣點 pi來自天空光某一方向發(fā)光帶的光強貢獻為：

t(papv,λ)光學(xué)深度的值與視點的高度、視線方向的天頂角有關(guān)，可以預(yù)先計算為以視點的高度和視線方向的天頂角為參數(shù)的一個二維紋理查找表。

F(θ,λ)、 F(θj,λ)散射相位函數(shù)的值與散射角有關(guān)，可以預(yù)先計算為以散射角為參數(shù)的一個一維紋理查找表。

ρpi、 ρpj,k大氣粒子密度比的值與高度相關(guān)，可以預(yù)先計算以高度為參數(shù)的一個一維紋理查找表。

由于大氣中不同種粒子的光學(xué)深度 t(S,λ)、散射相位函數(shù) F(θ,λ)、密度比ρ都是不同的，因此每個預(yù)計算的紋理查找表都需要存儲多個數(shù)值。本文的模型主要考慮了5種類型的大氣粒子，它們是空氣粒子、灰塵粒子、霧粒子、雨粒子、雪粒子，其中空氣分子發(fā)生Rayleigh散射，后面4種粒子發(fā)生Mie散射。每個紋理查找表的值都包括4個通道，可以存儲4個float值，前3個通道分別存儲與空氣分子、灰塵粒子、霧粒子相關(guān)的預(yù)計算結(jié)果，第4個通道存儲與雨粒子或者雪粒子相關(guān)的預(yù)計算結(jié)果，即雨粒子、雪粒子同時只能有一個有效。

通過上述4個預(yù)計算的紋理查找表，天空顏色模型表達(dá)式的每一項值都可以在運行中通過查表獲得，然后只需要通過簡單加減乘積運算就可以計算出天空顏色模型的值，因此該方法非常適合在GPU上執(zhí)行。算法的繪制在OpenGL片段處理器執(zhí)行，通過簡單地調(diào)用 texure1D、texure2D和texure3D等幾個著色器指令函數(shù)就可以直接從已經(jīng)預(yù)計算好的4個紋理查找表中獲取光學(xué)深度、散射相位函數(shù)、密度比的值，從而實現(xiàn)天空顏色模型的實時計算與繪制。

3 實驗及結(jié)果分析

本文算法的測試環(huán)境如下：硬件平臺為Intel酷睿i7 3770 CPU、4 GB內(nèi)存、NVIDIA GeForce GTX 660顯卡（2GB顯存），軟件平臺為Windows 7、OpenGL 4.1、GLEW 1.9。

以下所有實驗只考慮空氣分子和灰塵兩種大氣粒子，不考慮霧粒子、雨粒子、雪粒子。天空光發(fā)光帶數(shù)目 Nsky為8，大氣球殼分層數(shù)目N為20。空氣分子發(fā)生Rayleigh散射，均值高度 HR為7.994km，對應(yīng)于光線3個波長(680nm, 550nm, 440nm)的散射系數(shù)為(0.0058, 0.0135, 0.0331)；灰塵粒子發(fā)生Mie散射。

下面對3種不同方法在實際應(yīng)用中的繪制效率進行實驗分析。表1為相同場景環(huán)境下采用3種不同方法所消耗的繪制時間對比，生成圖像分辨率為1024×768。為了使仿真環(huán)境更加真實，場景地形數(shù)據(jù)庫為臺灣地區(qū)，地形圖像分辨率為5m，數(shù)據(jù)格式為TerraPage(.txp)，對txp格式的多線程頁面調(diào)度支持采用 OSG開源圖形引擎代碼。從繪制各階段消耗時間以及幀數(shù)結(jié)果可以看出，本文的多重散射方法相比Bruneton和Neyret[9]的方法在繪制速度上具有明顯的優(yōu)勢，基于顏色插值的方法[4]雖然繪制速度最快，但由于只是一個簡單經(jīng)驗方法，沒有任何基于大氣物理模型的仿真計算，不能模擬出真實的天空光照效果。

下面針對太陽以及灰塵粒子不同的高度繪制相應(yīng)的天空光照效果。

圖4為太陽不同高度下的天空光照效果圖，灰塵粒子均值高度 H 為1.2km，散射系數(shù)為

M0.02，消光系數(shù)為 0.0222，從繪制的效果可以看出本文的算法能夠?qū)崟r模擬不同太陽高度的真實天空光照效果。

圖5為灰塵粒子不同均值高度下的天空光照效果圖，灰塵粒子均值高度從左至右依次為0.6、1.2、3.0、4.0、5.0、6.0，單位為km，散射系數(shù)為0.02，消光系數(shù)為0.0222。從繪制的效果可以看出隨著粒子均值高度的增加，粒子的密度越來越大，相應(yīng)對太陽光的吸收和散射效果也越明顯。

下面針對灰塵粒子不同的散射系數(shù)和消光系數(shù)繪制相應(yīng)的天空光照效果。

圖6所示為灰塵粒子不同散射系數(shù)下的天空光照效果，灰塵粒子均值高度為1.2km，散射系數(shù)從左至右依次為0.005、0.010、0.020、0.030、0.040、0.050。從繪制的效果可以看出隨著散射系數(shù)的增加，太陽光在空氣中散射的效果越來越明顯。

圖7所示為灰塵粒子不同消光系數(shù)下的天空光照效果，灰塵粒子均值高度為1.2km，散射系數(shù)為 0.02，消光系數(shù)從左至右依次為0.4000、0.2500、0.2000、0.0667、0.0333、0.0222。從繪制的效果可以看出隨著消光系數(shù)的減少，光線在大氣中的衰減越來越小，相應(yīng)天空越來越明亮，顏色也越來越藍(lán)。

表1 不同方法繪制時間的對比

圖4 太陽不同高度下的天空光照效果

圖5 灰塵粒子不同均值高度下的天空光照效果

圖6 灰塵粒子不同散射系數(shù)下的天空光照效果

圖7 灰塵粒子不同消光系數(shù)下的天空光照效果

4 結(jié) 論

本文提出了一種多重散射的天空光照效果建模與實時繪制方法，實驗結(jié)果證明該方法能夠?qū)崟r模擬不同時間、不同大氣條件下的真實感天空光照效果，可以滿足各種實時應(yīng)用需求。

[1] Dobashi Y, Nishita T, Kaneda K, Yamashita H. A fast display method of sky color using basis functions [J]. The Journal of Visualization and Computer Animation, 1997, 8(2): 115-127.

[2] Dobashi Y, Yamamoto T, Nishita T. Interactive rendering of atmospheric scattering effects using graphics hardware [C]//HWWS 2002 Proceedings of the ACM SIGGRAPH/EUROGRAPHICS Conference on Graphics Hardware, 2002: 99-107.

[3] O'Neil S. Accurate atmospheric scattering [M]. GPU Gems 2: Programming Techniques for High-Performance Graphics and General-Purpose Computation. Addison-Wesley, 2005: 253-268.

[4] Abad J A. A fast, simple method to render sky color using gradients maps[EB/OL]. (2006-10-09) [2013-04-03]. http://www.geocities.com/ngdash/whitepapers/skydom ecolor.html.

[5] Schafhitzel T, Falk M, Ertl T. Real-time rendering of planets with atmospheres [C]//WSCG International Conference in Central Europe on Computer Graphics, Visualization and Computer Vision, 2007: 91-98.

[6] Preetham A J, Shirley P, Smits B. A practical analytic model for daylight [C]//SIGGRAPH 99, 1999: 91-100.

[7] Nishita T, Dobashi Y, Kaneda K. Display method of the sky color taking into account multiple scattering [C]//Pacific Graphics 1996, 1996: 117-132.

[8] Haber J, Magnor M, Seidel H P. Physically based simulation of twilight phenomena [J]. ACM Transactions on Graphics, 2005, 24(4): 1353-1373.

[9] Bruneton E, Neyret F. Precomputed atmospheric scattering[J]. Computer Graphics Forum, 2008, 27(4): 1079-1086.

[10] 劉世光, 彭群生. 氣象景觀的真實感模擬技術(shù)綜述[J]. 計算機輔助設(shè)計與圖形學(xué)學(xué)報, 2008, 20(4): 409-416.

[11] 錢 徽, 朱淼良. 大氣效果的繪制技術(shù)綜述[J]. 計算機輔助設(shè)計與圖形學(xué)學(xué)報, 2002, 14(5): 475-482.

[12] 王長波, 王章野, 曾 運, 彭群生. 考慮大氣折射的天空場景真實感繪制[J]. 計算機學(xué)報, 2005, 28(6): 939-949.

[13] 黃 雷, 王章野, 王長波. 雨霧天氣下光線散射效果的實時繪制[J]. 軟件學(xué)報, 2006, 17(11): 126-137.

[14] 杜 芳, 張 炎, 晁建剛, 王金坤. 基于GPU的地球大氣散射現(xiàn)象可視化仿真[J]. 系統(tǒng)仿真學(xué)報, 2009, 21(2): 147-150.

[15] 馮玉康, 周圣川, 馬純永, 韓 勇, 陳 戈. 基于GPU的地球大氣層和三維體積云仿真[J]. 計算機工程, 2012, 38(19): 218-225.

[16] Nishita T, Sirai T, Tadamura K. Display of the earth taking into account atmospheric scattering [C]// SIGGRAPH 93, 1993: 175-182.

[17] Hoffman N, Preetham A J. Rendering outdoor light scattering in real time [C]//Game Developer Conference, 2002: 337-352.

[18] Cornette W M, Shanks J G. Physical reasonable analytic expression for the single scattering phase function [J]. Applied Optics, 1992, 31(16): 3152-3160.

[19] O'Neil S. Real-time atmospheric scattering[EB/OL]. (2004-05-03)[2013-04-03]. http://www.gamedev.net/ page/resources/_/technical/graphics-programming-an d-theory/real-time-atmospatmos-scattering-r2093.

Modeling and Real-time Rendering of Sky Illumination Effects Taking Account of Multiple Scattering

Ai Zuliang, Zhang Limin, Yu Wenlong
(Naval Aeronautical Engineering Institute, Yantai Shandong 264001, China)

This paper presents a method of modeling and real-time rendering of the sky illumination effects taking account of multiple scattering with atmospheric scattering based on physical models. Firstly, with the sun and the sky as light sources, the model of the sky illumination effects taking account of multiple scattering is built. Then a reasonable segmented sampling strategy based on a variety of atmospheric particle densities is adopted to simplify the integral of the sky color model to reduce the integral approximation calculation error. By analyzing the simplified model, the pre-computed 2D texture and 3D texture for optical depth are proposed to avoid a run-time calculation of optical depth. Finally, the algorithm is executed on the fragment processor of GPU so that the sky illumination effects taking account of multiple scattering are rendered with high frame rate for real-time applications.

multiple scattering; sky illumination; optical depth; pre-compute; real-time rendering

TP 391.9

A

2095-302X (2014)02-0181-07

2013-04-03；定稿日期：2013-06-28

艾祖亮（1980-），男，湖北隨州人，博士研究生。主要研究方向為自然現(xiàn)象模擬、實時繪制技術(shù)。E-mail：aizuliang@live.cn