GPU加速的水體輻射傳輸Monte Carlo模擬模型研究*

杜克平，薛　坤

(1:北京師范大學(xué)地理學(xué)與遙感科學(xué)學(xué)院，遙感科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100875)(2:中國(guó)科學(xué)院南京地理與湖泊研究所，湖泊與環(huán)境國(guó)家重點(diǎn)實(shí)驗(yàn)室，南京 210008)(3:中國(guó)科學(xué)院大學(xué)，北京100049)

?

GPU加速的水體輻射傳輸Monte Carlo模擬模型研究*

杜克平1，薛坤2,3

(1:北京師范大學(xué)地理學(xué)與遙感科學(xué)學(xué)院，遙感科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100875)(2:中國(guó)科學(xué)院南京地理與湖泊研究所，湖泊與環(huán)境國(guó)家重點(diǎn)實(shí)驗(yàn)室，南京 210008)(3:中國(guó)科學(xué)院大學(xué)，北京100049)

水體輻射傳輸方程是復(fù)雜的微積分方程，只能利用數(shù)值方法求解，如Monte Carlo光線追蹤法、不變嵌入法、離散坐標(biāo)法等，其中，Monte Carlo方法是目前解決水體水下光場(chǎng)三維問題的唯一有效方法. 根據(jù)輻射傳輸理論，開發(fā)了水下光場(chǎng)的Monte Carlo模擬模型，主要包含大氣、水-氣界面、層化水體和水底邊界4個(gè)模塊. 實(shí)現(xiàn)了模擬任意太陽角度、不同水體固有光學(xué)屬性和任意深度條件下，考慮大氣、粗糙水面和水底邊界的水下光場(chǎng)，能夠獲取輻亮度、輻照度等輻射量的空間分布. 該模型暫不考慮Raman散射、偏振、內(nèi)部光源的影響. 實(shí)現(xiàn)了GPU加速水下光場(chǎng)Monte Carlo模擬，并用Mobley等提出的海洋光學(xué)標(biāo)準(zhǔn)問題中的問題1～6進(jìn)行驗(yàn)證. 在兩種計(jì)算環(huán)境下，通過對(duì)不同邊界條件下的CPU、GPU運(yùn)行時(shí)間及加速比的對(duì)比，發(fā)現(xiàn)GPU計(jì)算可以達(dá)到幾百至上千倍的加速比.

GPU；水體輻射傳輸方程；Monte Carlo模擬；水下光場(chǎng)

水下光場(chǎng)計(jì)算模型主要通過計(jì)算不同水深和太陽天頂角、方位角等條件下的輻照度、輻亮度、漫衰減系數(shù)等輻射量與表觀光學(xué)量來反映水下光場(chǎng)的分布. 精確模擬水下光場(chǎng)對(duì)研究海洋初級(jí)生產(chǎn)力、生物地球化學(xué)模型以及全球循環(huán)模型中的熱量收支平衡有重要意義.

水體輻射傳輸方程描述了電磁波在水中傳播時(shí)，受介質(zhì)的吸收、散射等作用的影響發(fā)生衰減的過程. 水體光學(xué)模型的核心就是求解水體輻射傳輸方程，由于方程本身的復(fù)雜性及邊界條件的限制，只能利用數(shù)值分析的方法得到數(shù)值解. 求輻射傳輸方程數(shù)值解的方法有Monte Carlo光線追蹤法[1-6]、不變嵌入法[7]、離散坐標(biāo)法[8]和矩陣算法[9-12]等.

Monte Carlo光線追蹤法直觀地表達(dá)輻射傳輸過程中光束被水中粒子(水分子、懸浮顆粒物等)吸收和散射過程，物理結(jié)構(gòu)清晰. 它可以對(duì)任意入射光線、散射相函數(shù)和固有光學(xué)屬性分布情況下的水下光場(chǎng)進(jìn)行三維模擬，可以解決一維模型所不能解決的問題，更加符合客觀世界的實(shí)際情況. 它具有較強(qiáng)的靈活性，可以模擬野外測(cè)量難以達(dá)到的極限條件[13]. 在水體光學(xué)研究領(lǐng)域，Kattawar[1-2]、Gordon[3]、Kirk[4]、Morel[5]、Mobley[6]等很多學(xué)者已經(jīng)用Monte Carlo方法解決了不同的水體光學(xué)問題. 國(guó)內(nèi)，曹文熙等運(yùn)用Monte Carlo方法模擬了浮標(biāo)浮體陰影及安裝儀器自陰影對(duì)水下光輻射測(cè)量的影響[14]. 唐軍武等[15]和凌在盈[13]利用Monte Carlo方法模擬了水下光場(chǎng)，但沒有考慮大氣和水底邊界狀況.

Monte Carlo最大的缺點(diǎn)是計(jì)算速度慢，即便有提高計(jì)算速度的方法(如方差縮減技術(shù)[2,16]、BMC[3]等)，但這些方法也只在一定程度上提高了光子的利用效率. Monte Carlo方法具有計(jì)算數(shù)據(jù)相對(duì)獨(dú)立、串行數(shù)據(jù)處理少等優(yōu)良的并行化特性，這種特性很適合用GPU加速運(yùn)算，基于GPU加速M(fèi)onte Carlo方法可以幾倍、幾十倍甚至上百倍地提高計(jì)算速度[17-20]，在圖形渲染、真實(shí)感場(chǎng)景模擬等應(yīng)用領(lǐng)域中應(yīng)用廣泛.

本文主要研究GPU加速的三維水體輻射傳輸Monte Carlo模擬模型，描述了水體輻射傳輸Monte Carlo模擬模型的建立過程，并實(shí)現(xiàn)了GPU加速. 比較分析了兩種計(jì)算環(huán)境下不同海洋光學(xué)標(biāo)準(zhǔn)問題及復(fù)雜水下光場(chǎng)的加速性能.

圖1 水體輻射傳輸Monte Carlo模擬模型框圖Fig.1　The flow chart of Monte Carlo simulation model of water

1　水體輻射傳輸 Monte Carlo模擬模型

水體輻射傳輸Monte Carlo模擬模型，主要包含大氣部分、水氣界面、層化水體和水底邊界4大模塊，其中大氣部分包括大氣模型和入射光場(chǎng)分布模型(圖1). 輸入?yún)?shù)包括光子數(shù)、入射角度、固有光學(xué)量、水深、風(fēng)速、大氣參數(shù)、底邊界等. 輸出結(jié)果包括輻亮度(L)、輻照度(E)、遙感反射率(Rrs)、遙感反射比(R)、平均余弦(μ)、漫衰減系數(shù)(K)等輻射量與表觀光學(xué)量. 該模型實(shí)現(xiàn)了模擬任意太陽角度(天頂角和方位角)、不同水體固有光學(xué)屬性(吸收系數(shù)、散射系數(shù)和散射相函數(shù))、任意深度條件下，考慮大氣、粗糙水面、層化水體和水底邊界的水下光場(chǎng)，可以得到輻亮度、輻照度等輻射量的空間分布. 該模型暫不考慮Raman散射、偏振和內(nèi)部光源的影響.

1.1大氣模型及入射光場(chǎng)分布

本文采用的大氣模型為RADTRAN-X模型[21-22]. 該模型需要的大氣參數(shù)包括：太陽角度、大氣壓強(qiáng)、大氣類型、相對(duì)濕度、降雨量、風(fēng)速、能見度、臭氧含量、氣溶膠參數(shù)等，得到特定條件下的漫射天空光輻照度(Esky)、太陽直射光輻照度(Edir)和水面之上總?cè)肷漭椪斩?Etot=Esky+Edir).

光子到達(dá)水-氣界面之前，與大氣中粒子發(fā)生相互作用，傳播方向發(fā)生變化. 實(shí)際情況下，很難得到天空漫射光的空間分布，需要假設(shè)入射光場(chǎng)分布. 入射光場(chǎng)分布模型有理想模型和半分析模型，其中常用的理想模型有各向同性模型和心型分布模型[22]. 本文使用由Harrison等提出的輻亮度分布半分析模型[23]，得到上半球空間的入射輻亮度分布. 為了方便統(tǒng)計(jì)、記錄光子信息，把2π立體角上半球空間分為天頂角間隔10°、方位角間隔15°的格網(wǎng)[7]，入射光的初始方向?yàn)楦窬W(wǎng)的中心對(duì)應(yīng)的天頂角和方位角.

1.2粗糙水面

光子到達(dá)水-氣界面會(huì)發(fā)生折射、反射. 平靜水面情況下，根據(jù)Snell定律和菲涅爾反射定律得到折射角和反射率. 然而，實(shí)際的海面受風(fēng)浪、波動(dòng)等影響，幾乎不可能是水平的，需要考慮粗糙水面的情況.

粗糙水面的模擬使用Cox等[24]提出的小波面斜率正態(tài)分布模型，不考慮重力波. 在模擬過程中，考慮小波面的斜率的概率分布計(jì)算、相互遮擋問題以及多次反射等問題. 此處，存在兩個(gè)坐標(biāo)系統(tǒng)：系統(tǒng)坐標(biāo)系和小波面坐標(biāo)系，需要進(jìn)行坐標(biāo)系轉(zhuǎn)換. 系統(tǒng)坐標(biāo)系以水平面垂直向下為Z軸的正方向，X-Y平面是水平面. 小波面坐標(biāo)系是以小波面的法線作為Z軸，X-Y平面在小波面所在的平面內(nèi).

θn是小波面的法線與水平面的法線的夾角，代表小波面的傾斜角度，θn=0時(shí)，小波面水平.θn的方差為σ2=0.003+0.00512U，U為風(fēng)速.

1.3水體內(nèi)光子與粒子的相互作用

光子與水中粒子發(fā)生相互作用時(shí)，一部分被吸收，為了減小Monte Carlo統(tǒng)計(jì)噪聲，本文采用變權(quán)重的光子跟蹤技術(shù)，即每次光子的權(quán)重乘以單次散射反照率(ω0)作為被吸收后光子的權(quán)重. 其中ω0的值越接近于1表示水體散射能力越強(qiáng)，越接近于0表示水體吸收能力越強(qiáng).

相函數(shù)是光子與粒子相互作用后散射方向的概率密度函數(shù). 實(shí)際模擬中，需要根據(jù)水體的類型和水體中組分的含量選擇合適的散射相函數(shù). 水中粒子的散射主要分為純水的散射和懸浮顆粒物的散射. 純水的散射相函數(shù)采用近似瑞利散射相函數(shù)，懸浮顆粒物的散射相函數(shù)可使用Petzold在1972年測(cè)量的不同水體的體散射函數(shù)的平均值，或者使用其它相函數(shù)公式(如：Henyey-Greenstein相函數(shù)、Fournier-Forand相函數(shù))[25]等.

針對(duì)垂向不均一的水體，把野外實(shí)際測(cè)量的IOP隨深度的廓線或者水體層化模型[26]輸入水體輻射傳輸Monte Carlo模擬模型，可以得到垂向異質(zhì)水體的水下光場(chǎng)分布.

1.4水底邊界

光子到達(dá)水底邊界部分被水底物質(zhì)吸收，部分被反射，反射的方向和反射率與水底狀況有關(guān)，一般用二向反射函數(shù)(BRDF)定量描述入射光線與出射光線的關(guān)系. 由于水底BRDF測(cè)量存在一定的困難，在研究中通常忽略水底邊界的二向性，假設(shè)水底表面是反射率固定的朗伯體，反射方向?yàn)楦飨蛲?

2　GPU加速的水體輻射傳輸Monte Carlo模擬

最初，GPU(Graphic Processing Unit，圖形處理器)是一個(gè)高度并行的用于圖形渲染任務(wù)的計(jì)算設(shè)備，現(xiàn)在已經(jīng)發(fā)展成為常用的計(jì)算處理器. 在數(shù)據(jù)并行運(yùn)算方面，GPU已經(jīng)得到了廣泛的應(yīng)用，也變得越來越高效.

水體輻射傳輸Monte Carlo(MC)模擬模型基于概率統(tǒng)計(jì)的方法，模擬光線在水體中與水中粒子相互作用的過程，獲取水體的表觀光學(xué)量. 要模擬得到精確的水下光場(chǎng)，需要初始設(shè)置大量光子(107)，傳統(tǒng)利用CPU計(jì)算時(shí)需要依次計(jì)算每個(gè)光子在水體中的傳播過程，耗時(shí)較長(zhǎng)，運(yùn)行效率低. 然而，根據(jù)MC模擬過程的可并行計(jì)算的特點(diǎn)，引入GPU計(jì)算可以同時(shí)并行計(jì)算大量光子，大幅提高計(jì)算效率. 在基于CPU的水體輻射傳輸Monte Carlo模型的基礎(chǔ)上，按照并行計(jì)算的要求，把可以并行計(jì)算的部分放到GPU上運(yùn)行，條件判斷、邊界確定、結(jié)果保存等不具備并行特性的部分在CPU上計(jì)算.

對(duì)CPU-MC程序進(jìn)行改寫可實(shí)現(xiàn)MC模擬模型的GPU并行運(yùn)算. 表1列出了兩種計(jì)算環(huán)境的CPU、GPU配置及對(duì)應(yīng)的軟件環(huán)境.

3　結(jié)果驗(yàn)證與加速性能分析

3.1海洋光學(xué)標(biāo)準(zhǔn)問題驗(yàn)證

Mobley等提出了海洋光學(xué)的7個(gè)標(biāo)準(zhǔn)問題[26]，包括理想的簡(jiǎn)單問題、基本問題、層化水體問題、大氣效應(yīng)問題、粗糙水面問題、水底邊界問題和拉曼散射問題. 他們用7種輻射傳輸模型(Hydrolight、離散坐標(biāo)法和Monte Carlo方法)，得到特定條件下的輻射量(輻亮度和輻照度)的平均值，經(jīng)常被國(guó)內(nèi)外學(xué)者作為檢驗(yàn)海洋輻射傳輸模型精度的標(biāo)準(zhǔn)[9,15].

表1　計(jì)算環(huán)境描述

海洋光學(xué)7個(gè)標(biāo)準(zhǔn)問題中的基本問題假設(shè)為：(1)水體表面水平；(2)水體水平與垂向均一且無限水深；(3)不考慮天空光；(4)太陽光為點(diǎn)光源，天頂角為60°；(5)太陽光在其垂直表面上的輻照度為1 W/(m2·nm)；(6)水體內(nèi)部無Raman散射及內(nèi)部光源存在；(7)假設(shè)水體為高散射(ω0=0.9)或者高吸收(ω0=0.2)水體兩種情況；(8)水體散射包括水分子和顆粒物散射兩部分. 表2是對(duì)問題1～6的簡(jiǎn)單描述，參照文獻(xiàn)[26]中的表3，本文沒有考慮問題7(Raman散射).

表2　海洋光學(xué)標(biāo)準(zhǔn)問題中問題1～6的描述[26]

利用GPU加速的水體輻射傳輸Monte Carlo模擬模型，對(duì)表2的6個(gè)問題分別進(jìn)行水下光場(chǎng)模擬. 表3是入射光子數(shù)為107，不同標(biāo)準(zhǔn)問題、光學(xué)深度、ω0下的模擬結(jié)果(下行輻照度(Ed)、上行輻照度(Eou)和向上輻亮度(Lu))與Mobley等[26]的相應(yīng)模擬結(jié)果的平均值(Ed_aver、Eou_aver和Lu_aver)，以及它們之間的平均相對(duì)誤差(δ)(式1)，x分別為Ed、Eou和Lu.

(1)

從對(duì)6個(gè)標(biāo)準(zhǔn)問題驗(yàn)證的結(jié)果(表3)可以得出，水體輻射傳輸模擬模型的精度基本達(dá)到了Mobley等的研究[26]中多個(gè)數(shù)值計(jì)算模型計(jì)算得到的平均值的精度范圍.Ed、Eou比Lu的相對(duì)誤差低，是因?yàn)長(zhǎng)u是由落入某個(gè)角度范圍內(nèi)少數(shù)光子計(jì)算得到的，受光子數(shù)目影響，穩(wěn)定性差，而Eou與整個(gè)半球空間中的光子數(shù)有關(guān)，結(jié)果相對(duì)穩(wěn)定. 在光學(xué)深度較大或強(qiáng)吸收水體存在誤差較大的情況下(問題3，在60 m處模擬結(jié)果誤差達(dá)到20%以上)，是由于Monte Carlo方法本身存在噪聲，復(fù)雜水體條件下假設(shè)的隨機(jī)數(shù)較多，同樣光子數(shù)的情況下，在強(qiáng)吸收水體噪聲較大. 因此，欲得到較好的結(jié)果，在強(qiáng)吸收水體需要增加模擬的光子數(shù)以保證模擬精度.

3.2不同模擬條件的GPU加速性能分析

對(duì)海洋光學(xué)基本問題及其組合邊界條件下的CPU計(jì)算時(shí)間、GPU加速的運(yùn)行時(shí)間和加速比進(jìn)行統(tǒng)計(jì)，總光子數(shù)為106.

通過GPU加速，在問題1和問題3中加速比為300～500，而問題2和問題6達(dá)到了1000多(表4)，這表明GPU加速可以明顯提高M(jìn)onte Carlo算法的運(yùn)行效率. 問題2(ω0=0.9)CPU計(jì)算需要2 d(173604 s/3600 s=48.22 h)，而GPU加速之后只需要加入大氣模型及天空光之后，GPU加速?zèng)]有明顯提高運(yùn)算速度，2～3 min(142 s). 從CPU的每個(gè)光子依次計(jì)算，到GPU的許多光子同時(shí)在GPU中計(jì)算，節(jié)省了時(shí)間，大幅提高了運(yùn)行效率. 與文獻(xiàn)[26]中的表4中5個(gè)MC模型對(duì)問題3的計(jì)算時(shí)間進(jìn)行比較，運(yùn)行速度平均提高了236倍.

表3　本文模擬結(jié)果與Mobley等[26]研究結(jié)果的平均值及平均相對(duì)誤差(入射光子數(shù)為107)

表4 兩種計(jì)算環(huán)境下不同問題的CPU和GPU計(jì)算加速比

z表示ω0隨深度變化.

加入大氣模型及天空光之后，GPU加速?zèng)]有明顯提高運(yùn)算速度，反而在一些情況下，計(jì)算效率降低. 問題4僅能達(dá)到幾倍到十幾倍的加速比，在考慮各種邊界條件的復(fù)雜水下光場(chǎng)(問題3、4、5、6)模擬時(shí)，GPU運(yùn)算反而降低了運(yùn)行效率. 原因是在大氣入射的上半球空間里，總?cè)肷涔庾訑?shù)根據(jù)天空光的分布分配到每個(gè)quad(9*24)里，而在非直射方向quad里的光子數(shù)量較少，應(yīng)用GPU加速反而降低了運(yùn)行時(shí)間.

3.3不同水體散射特性下的加速性能

從以上分析可以看出，在同一個(gè)標(biāo)準(zhǔn)問題中，ω0=0.2和ω0=0.9的運(yùn)行時(shí)間相差很大. 因此，針對(duì)不同散射特性條件(ω0=0.1～0.9)下的Monte Carlo模擬的運(yùn)行時(shí)間和效率進(jìn)行了比較. 結(jié)果(圖2)表明，兩種計(jì)算環(huán)境下，CPU運(yùn)行時(shí)間相差不大，而計(jì)算環(huán)境1的GPU運(yùn)行時(shí)間大約為計(jì)算環(huán)境2的2倍，也就是計(jì)算環(huán)境2的加速比是計(jì)算環(huán)境1的2倍. 在同一計(jì)算環(huán)境下，水體從高吸收到高散射的過程中，CPU和GPU的運(yùn)行時(shí)間呈指數(shù)增加. 在不同散射程度的水體，光子的壽命不同，在高吸收水體，光子存活時(shí)間較短. 因此，在總?cè)肷涔庾訑?shù)一定的情況下，高吸收水體的運(yùn)行時(shí)間會(huì)短，光子的利用率下降，導(dǎo)致模擬結(jié)果的噪聲大于高散射水體.

圖2 問題2 (a)、5 (b)和6 (c)在兩種計(jì)算機(jī)環(huán)境下的CPU和GPU運(yùn)行時(shí)間及加速比(總光子數(shù)為106)隨ω0的變化Fig.2　Computing time and speedup ratio vary with ω0 under two computer environments of problem 2 (a), problem 5 (b) and problem 6 (c)

4　結(jié)論

根據(jù)輻射傳輸理論，利用Monte Carlo模擬方法，開發(fā)了包括大氣、水-氣界面、層化水體、水底邊界4大模塊的水下光場(chǎng)模擬模型. 并利用Mobley等發(fā)表的海洋光學(xué)標(biāo)準(zhǔn)問題中的問題1～6對(duì)模型的精度進(jìn)行驗(yàn)證. 引入GPU加速計(jì)算技術(shù)，使得水下光場(chǎng)Monte Carlo模擬的加速比得到大幅提高，一些簡(jiǎn)單問題的加速比可以達(dá)到幾百到上千. 對(duì)于不同散射特性的水體，因?yàn)楣庾永眯什煌?，需要加大模擬光子數(shù)來保證運(yùn)行結(jié)果的準(zhǔn)確性. 但是，對(duì)于包含大氣模型和入射光場(chǎng)的水下光場(chǎng)Monte Carlo模擬方面，因?yàn)槟Ｐ托枰A(yù)先判斷每個(gè)quad里的光子數(shù)，導(dǎo)致大量光子一起運(yùn)算的并行特性被削弱，GPU并行加速的優(yōu)勢(shì)沒有得到很好的體現(xiàn). 以后研究中，將進(jìn)一步解決大氣入射光場(chǎng)耦合水下光場(chǎng)的Monte Carlo模擬模型的GPU加速問題.

致謝：感謝中國(guó)科學(xué)院南京地理與湖泊研究所馬榮華研究員提供的建議和幫助，以及可以應(yīng)用GPU加速的工作站作為本文的計(jì)算環(huán)境2. 感謝美國(guó)University of Massachusetts Boston的Zhongping Lee教授對(duì)本研究的建議與幫助.

[1]Plass GN, Kattawar GW. Monte Carlo calculations of light scattering from clouds.AppliedOptics, 1968, 7(3): 415-419.

[2]Plass GN, Kattawar GW. Radiative transfer in an atmosphere-ocean system.AppliedOptics, 1969, 8(2): 455-466.

[3]Gordon HR. Ship perturbation of irradiance measurements at sea. 1: Monte Carlo simulations.AppliedOptics, 1985, 24(23): 4172-4182.

[4]Kirk JTO. Monte Carlo study of the nature of the underwater light field in, and the relationships between optical properties of, turbid yellow waters.MarineandFreshwaterResearch, 1981, 32(4): 517-532.

[5]Morel A, Gentili B. Diffuse reflectance of oceanic waters: its dependence on Sun angle as influenced by the molecular scattering contribution.AppliedOptics, 1991, 30(30): 4427-4438.

[6]Mobley CD, Sundman LK. Effects of optically shallow bottoms on upwelling radiances: Inhomogeneous and sloping bottoms.LimnologyandOceanography, 2003, 48(1): 329-336.

[7]Mobley CD. Light and water: Radiative transfer in natural waters. New York: Academic Press, 1994.

[8]Jin Z, Stamnes K. Radiative transfer in nonuniformly refracting layered media: atmosphere-ocean system.AppliedOptics, 1994, 33(3): 431-442.

[9]何賢強(qiáng), 潘德爐, 白雁等. 基于矩陣算法的海洋-大氣耦合矢量輻射傳輸數(shù)值計(jì)算模型. 中國(guó)科學(xué): D輯: 地球科學(xué), 2006, (9): 860-870.

[10]張鑒, 何曉雄, 趙鳳生. 利用大氣-海洋系統(tǒng)輻射傳輸模擬水色遙感信息量的變化特性. 量子電子學(xué)報(bào), 2003, (5): 623-628.

[11]何賢強(qiáng), 潘德爐, 白雁等. 基于輻射傳輸數(shù)值模型PCOART的大氣漫射透過率精確計(jì)算. 紅外與毫米波學(xué)報(bào), 2008, (4): 303-307.

[12]何賢強(qiáng), 潘德爐, 白雁等. 海洋-大氣耦合矢量輻射傳輸粗糙海面模型. 光學(xué)學(xué)報(bào), 2010, (3): 618-624.

[13]凌在盈. 基于Monte Carlo方法的水體二向反射分布函數(shù)(BRDF)模擬[學(xué)位論文]. 杭州： 浙江大學(xué), 2010.

[14]曹文熙, 吳廷芳, 楊躍忠等. 光學(xué)浮標(biāo)陰影效應(yīng)的蒙特卡洛模擬. 高技術(shù)通訊, 2003,(3): 80-84.

[15]唐軍武, 田國(guó)良, 陳清蓮. 離水輻射非朗伯特性的Monte Carlo模擬及分析. 海洋學(xué)報(bào): 中文版, 2000, (2): 48-57.

[16]Plass GN, Kattawar GW, Guinn JJA. Radiative transfer in the earth’s atmosphere and ocean: influence of ocean waves.AppliedOptics, 1975, 14(8): 1924-1936.

[17]Fang Q, Boas DA. Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units.OpticsExpress, 2009, 17(22): 20178-20190.

[18]Preis T, Virnau P, Paul Wetal. GPU accelerated Monte Carlo simulation of the 2D and 3D Ising model.JournalofComputationalPhysics, 2009, 228(12): 4468-4477.

[19]Ren N, Liang J, Qu Xetal. GPU-based Monte Carlo simulation for light propagation in complex heterogeneous tissues.OpticsExpress, 2010, 18(7): 6811-6823.

[20]Wang Y, Li P, Jiang Cetal. GPU accelerated electric field Monte Carlo simulation of light propagation in turbid media using a finite-size beam model.OpticsExpress, 2012, 20(15): 16618-16630.

[21]Gregg WW, Carder KL. A simple spectral solar irradiance model for cloudless maritime atmospheres.LimnologyandOceanography, 1990, 35(8): 1657-1675.

[22]Mobley CD, Sundman LK. HydroLight 5-EcoLight 5 Technical Documentation. Sequoia Scientific, Inc, 2008.

[23]Harrison AW, Coombes CA. An opaque cloud cover model of sky short wavelength radiance.SolarEnergy, 1988, 41(4): 387-392.

[24]Cox C, Munk W. Measurement of the roughness of the sea surface from photographs of the Suns Glitter.JournaloftheOpticalSocietyofAmerica, 1954, 44(11): 838-850.

[25]Mobley CD, Sundman LK, Boss E. Phase function effects on oceanic light fields.AppliedOptics, 2002, 41(6): 1035-1050.

[26]Mobley CD, Gentili B, Gordon HRetal. Comparison of numerical models for computing underwater light fields.AppliedOptics, 1993, 32(36): 7484-7504.

Accelerating Monte Carlo radiative transfer simulation of water using GPU technique

DU Keping1& XUE Kun2,3

(1:StateKeyLaboratoryofRemoteSensingScience,SchoolofGeography,BeijingNormalUniversity,Beijing100875,P.R.China)(2:StateKeyLaboratoryofLakeScienceandEnvironment,NanjingInstituteofGeographyandLimnology,ChineseAcademyofSciences,Nanjing210008,P.R.China)(3:UniversityofChineseAcademyofSciences,Beijing100049,P.R.China)

Radiative Transfer Equations (RTE) of water body are complex integro-differential equations, which can be solved by different numerical methods, e.g., Monte Carlo ray tracing, invariant imbedding, and discrete ordinates. Monte Carlo method is a powerful technique, which can be used in any water body, even those whose boundary conditions and inherent optical properties (IOPs) vary in three dimensions. However, the Monte Carlo method is computationally costly, which limits the use for many problems in optical oceanography. In this paper, a new kind of acceleration technology to accelerate the ocean radiative transfer simulation, using the CUDA-enabled graphics processing unit (GPU) is presented. With the approach’s help, it is easy to code on NVIDIA GPUs and there is no need to worry about the hardware details of a specific GPU. Firstly, some basic ideas of ocean radiative Monte Carlo simulation are introduced, then GPU programs for ocean radiative transfer simulation are implemented. Finally, the performances of the two GPUs (NVIDIA GTX 670 GPU and NVIDIA Quadro 6000 GPU) with their CPU counterparts are compared. From our numerical results, the speedup over hundreds of times for solving the issues is achieved compared with that obtained using CPU.

Graphics processing unit (GPU); radiative transfer equation; Monte Carlo simulation; underwater light field

J.LakeSci.(湖泊科學(xué))， 2016， 28(3)： 654-660

10.18307/2016.0322

?2016 byJournalofLakeSciences

*國(guó)家自然科學(xué)基金項(xiàng)目(41471284，41431176)、國(guó)家高技術(shù)研究發(fā)展計(jì)劃“863”項(xiàng)目(2012AA12A303)和國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃“973”項(xiàng)目(2013CB733403)聯(lián)合資助. 2015-06-03收稿；2015-09-17收修改稿. 杜克平(1974～)，男，博士，副教授；E-mail：kpdu@bnu.edu.cn.