逆正向蒙特卡羅方法及應(yīng)用

施國棟

(上海財經(jīng)大學(xué)金融學(xué)院,上海 200433)

0 引言

蒙特卡羅方法在確定性問題、粒子運(yùn)輸、大氣環(huán)境模擬、數(shù)理統(tǒng)計(jì)學(xué)、金融經(jīng)濟(jì)學(xué)和科學(xué)實(shí)驗(yàn)?zāi)M等領(lǐng)域有著廣泛應(yīng)用。蒙特卡羅方法計(jì)算結(jié)果精確,當(dāng)所研究的問題沒有解析解時,蒙特卡羅方法的結(jié)果往往可以作為基準(zhǔn)解[1]。蒙特卡羅方法計(jì)算比較靈活,在解決介質(zhì)的輻射特性問題時,根據(jù)光路可逆性原理,可以分別使用正向或逆向蒙特卡羅方法模擬光線的傳播與能量分配[2,3]。然而對于一些不可逆的問題,即需要根據(jù)事物當(dāng)前的狀態(tài)值計(jì)算下一個狀態(tài),則不可以用逆向蒙特卡羅方法計(jì)算,因此需要對傳統(tǒng)的單向蒙特卡羅方法進(jìn)行改進(jìn)。

目前考慮偏振效應(yīng)的輻射傳輸研究中,Lau與Watson[4]導(dǎo)出了梯度折射率介質(zhì)中偏振光輻射傳輸方程。Zhao等[5,6]完整地推導(dǎo)了梯度折射率介質(zhì)中矢量輻射傳輸方程。Ben等[7]采用正向蒙特卡羅法求解了矢量輻射傳輸問題,提出了一種模擬梯度介質(zhì)中偏振輻射傳輸?shù)姆椒ā?/p>

在非均勻折射率介質(zhì)中光線軌跡通常是彎曲的,光線在傳輸過程中會出現(xiàn)會聚或發(fā)散,即使在透明介質(zhì)中,光線傳播也會引起偏振橢圓旋轉(zhuǎn),從而引起能量的重新分配。橢圓偏振的傳輸是由Rytov場矢量旋轉(zhuǎn)法描述[8],所以光的偏振狀態(tài)在梯度折射率介質(zhì)中會不斷變化,并影響下一時刻的結(jié)果,這使得對光線反向跟蹤變得不可能。反向蒙特卡羅法可以從某個特定角度出發(fā)對光線進(jìn)行跟蹤,從而極大減少了跟蹤光束的數(shù)量[9-11]?？紤]偏振狀態(tài)的輻射傳輸問題,如果可以通過反向蒙特卡羅法求解,可以有效減少計(jì)算時間和避免角度離散,從而獲得更高精度的結(jié)果。

近幾十年來,一些學(xué)者試圖采用反向計(jì)算方法來解決矢量輻射傳輸?shù)膯栴}[10,12-14]。但是光子的散射方向受偏振狀態(tài)影響較大,這些近似方法可能無法提供準(zhǔn)確的結(jié)果[15]。為了解決這類復(fù)雜問題,本文發(fā)展了更為靈活的逆正向蒙特卡羅法(RFMC),打破了傳統(tǒng)單向蒙特卡羅方法計(jì)算的連貫性,模擬計(jì)算了一維梯度折射率介質(zhì)自身的Stokes矢量。

1 逆正向蒙特卡羅法

1.1 物理模型

圖1為一維梯度折射率半透明介質(zhì)的多層模型,假設(shè)介質(zhì)為無限大的平板,厚度為L,折射率的分布沿著z軸方向,下表面處折射率為n0,上表面處折射率為nL,這種模型與大氣層的梯度折射率分布類似。介質(zhì)被假定為吸收、發(fā)射、散射或非散射性介質(zhì),吸收系數(shù)為κa,介質(zhì)的溫度場為等溫分布,灰色半透明介質(zhì)和純鏡面反射邊界符合Snell定律和Fresnel定律。研究的重點(diǎn)為梯度折射率介質(zhì)模型以及上下表面的表觀發(fā)射率。

圖1 一維梯度折射率半透明介質(zhì)多層模型Fig.1 Multilayer model of one-dimensional gradient refractive index translucent media

光線在一維梯度折射率介質(zhì)中軌跡的跟蹤可以使用簡單的多層近似方法[16],每個子層的折射率近似等于該子層中心點(diǎn)折射率。基于該假設(shè),每一層都有固定的折射率,光線在每個子層中沿直線傳播,兩層之間虛擬的界面被假定為折射或全反射界面[17]。利用蒙特卡羅方法對光線模擬時,將光線跟蹤的過程分成發(fā)射、傳播、散射、反射、吸收等一系列子過程。

1.2 Stokes矢量計(jì)算子過程

光線在非散射一維梯度折射率介質(zhì)中的軌跡為一個平面,從文獻(xiàn)[4,5]可以得到

式中S為Stokes參數(shù),S=(I,Q,U,V)[18],I是輻射強(qiáng)度,Q是垂直和水平方向線偏振強(qiáng)度,U是±45°方向上偏振強(qiáng)度,V是圓偏振強(qiáng)度;s為光線軌跡弧長;n為折射率。

光線在經(jīng)過介質(zhì)反射、折射和散射后的Stokes矢量為Sref、Srefr和Ss,其計(jì)算公式分別為

式中Si是入射Stokes矢量,R(θi)為反射矩陣,T(θi)為折射矩陣,θi為入射角,M(Θ)為單次散射的Mueller矩陣,Θ為散射角,L表示沿著光線傳播方向的旋轉(zhuǎn)矩陣,i1、i2是圖2中旋轉(zhuǎn)的角度。下面分別列出R(θi)、T(θi)和M(Θ)的計(jì)算方法。

式(2)中的反射矩陣R(θi)可以表示為[19-21]

式中θi為入射角,rp和rs分別表示平行和垂直偏振的振幅反射率,Re()和Im()分別表示復(fù)數(shù)的實(shí)部和虛部,*表示共軛復(fù)數(shù)。rp和rs分別用菲涅耳方程來計(jì)算,即

式中ni和nt分別為入射邊的折射率和折射邊的折射率。

式(3)中的折射矩陣T(θi)可以表示為

式中θt為折射角,tp和ts分別為p和s偏振的透射系數(shù),其計(jì)算公式分別為

矢量輻射傳輸中散射比較復(fù)雜,由Van de Hulst定義[22],M(Θ)的形式為

圖2為入射光與散射光線的幾何關(guān)系。Si(θ1,φ1)和Ss(θ2,φ2)分別是光線散射前后的方向,Θ為散射角,即入射和散射方向之間的夾角。矩陣L表示沿著光線傳播方向順時針旋轉(zhuǎn)的旋轉(zhuǎn)矩陣,其表達(dá)式為[23]

式中Φ(等于i1或者i2)是圖中旋轉(zhuǎn)的角度。

相關(guān)的一些角度可以由圖中幾何關(guān)系得到。圖2中入射Stokes矢量Si和z軸確定平面aoz,散射后的Stokes矢量Ss和z軸確定平面boz,Si和Ss確定了散射平面。i1是平面aoz與散射平面的夾角,i1∈[0,2π),服從均勻概率分布,i2是平面boz與散射平面旋轉(zhuǎn)的夾角。在不考慮矢量輻射傳輸?shù)臉?biāo)量模擬計(jì)算中,光線的散射角通常采用反向插值方法計(jì)算[24],而在考慮矢量輻射傳輸時,可采用拒絕法確定光線的散射方向。

圖2 入射光與散射光的幾何關(guān)系Fig.2 Geometric relation of incident light and scattered light

1.3 拒絕法確定散射方向

拒絕法是一種產(chǎn)生給定分布函數(shù)的方法。在圖2中已知入射光線Si(θ1,φ1)的天頂角和圓周角,散射方向i1和Θ通過拒絕法來確定[25,26],其他角度i2、θ2和φ2可以用空間余弦定理來確定[27]。

采用拒絕法確定散射方向i1和Θ的步驟為:

1)散射角Θ和i1在局部坐標(biāo)系中定義,其計(jì)算公式分別為

式中RΘ和Rφ是均勻分布的隨機(jī)數(shù),隨機(jī)范圍在[0,1]區(qū)間。

2)計(jì)算Is′和Is,max′。

在散射平面定義Stokes矢量,散射后Stokes矢量可以表示為

將式(12)展開,得到

將步驟1)中生成的散射角Θ和i1代入式(12),得到Ss′中I部分的分量,即為Is′。

根據(jù)三角函數(shù)公式將式(13)變形為

得到式(14)中散射角Θ對應(yīng)的最大值Is,max′[7]。

3)采用拒絕法來判斷生成的散射方向Θ和i1是否有效。

比較步驟2)中的Is′和ζIs′,max,其中ζ是[0,1]區(qū)間內(nèi)均勻分布的隨機(jī)數(shù)。如果Is′≥ζIs′,max,那么步驟1)中生成的散射方向是有效的;反之則步驟1)中生成的散射方向是無效的,重復(fù)步驟1)直到生成新的散射方向被接受。

1.4 逆正向蒙特卡羅法計(jì)算流程

Stokes矢量的第一部分分量I表示光子所攜帶的能量,而其他部分(Q,U,V)表示光線的偏振狀態(tài)。采用蒙特卡羅法模擬時,為了確保能量守恒,需要對Stokes矢量進(jìn)行歸一化處理,即在光線軌跡上,Stokes矢量中第一部分I保持為1。Stokes矢量的歸一化公式為

在一維非散射梯度折射率介質(zhì)中光線軌跡是一個平面,根據(jù)式(1)和式(15),光線在介質(zhì)中傳播的子過程中Stokes矢量雖然不斷變化,但是在歸一化后其S′值不變,即在一維非散射梯度折射率介質(zhì)的子層傳播過程中Stokes矢量不發(fā)生變化。

但不能使用傳統(tǒng)標(biāo)量反向蒙特卡羅法,主要原因是:

1)在矢量輻射傳輸問題中,計(jì)算反射系數(shù)的公式與標(biāo)量輻射傳輸問題中不同。當(dāng)光線從介質(zhì)1入射到介質(zhì)2時,考慮偏振的反射率ρf為

式中θ為入射角;ρ是標(biāo)量方法中的反射率;Q是入射Stokes矢量Si中的第二部分分量;rp和rs分別是平行和垂直偏振的振幅反射率,采用菲涅耳公式(6)來計(jì)算。雖然Q分量最終引起的反射率變化不是很大,但是考慮偏振的反射率公式和標(biāo)量方法的反射率公式不一樣。

2)在考慮偏振時確定散射方向與標(biāo)量方法不同?？紤]偏振時散射方向Θ和i1是否有效,可以使用拒絕法來判斷。在式(13)中,考慮偏振時的散射方向與Ii、Qi、Ui、M1(Θ)、M2(Θ)這些量都有關(guān);而只考慮標(biāo)量問題時散射方向的確定與散射相函數(shù)有關(guān)。在各向異性散射的問題中,Qi、Ui不等于0,因此各向異性散射介質(zhì)考慮散射偏振問題的結(jié)果與只計(jì)算標(biāo)量問題的結(jié)果不同。

根據(jù)以上兩點(diǎn),不能直接使用標(biāo)量反向蒙特卡羅法。因此提出了逆正向蒙特卡羅法,即在光線跟蹤時采用逆向跟蹤方法,通過反向光線跟蹤找到光線在介質(zhì)內(nèi)的終點(diǎn),然后在Stokes矢量計(jì)算時則采用正向計(jì)算。在正向Stokes矢量計(jì)算時式(2)不變,因?yàn)槿肷浣呛头瓷浣窍嗟?式(3)原折射角變?yōu)槿肷浣?i1和i2相互交換,坐標(biāo)旋轉(zhuǎn)方向由順時針方向變化為逆時針方向。式(4)則變?yōu)?/p>

逆正向蒙特卡羅法的理論框架應(yīng)遵循“正向”蒙特卡羅法,逆向光線跟蹤的過程只是在介質(zhì)中找到光的終點(diǎn),然后該終點(diǎn)變成下一步正向計(jì)算時光線的起點(diǎn)。

圖3為逆正向蒙特卡羅法的光線跟蹤示意圖。圖中黑色實(shí)線為光線反向跟蹤的軌跡,灰色實(shí)線是光線正向Stokes矢量計(jì)算過程。利用逆正向蒙特卡羅法求解考慮偏振狀態(tài)的表觀發(fā)射率的算法,可以概括為如下步驟:

1)從反向跟蹤的起點(diǎn)開始,模擬光線的總數(shù)量為N,假設(shè)光線全部進(jìn)入介質(zhì)內(nèi)而不發(fā)生反射。跟蹤光線并記錄各子過程的相關(guān)信息,如入射方向、折射方向、折射率、散射方向等,以用于正向Stokes矢量計(jì)算。

2)每一束模擬光線在被散射和吸收之前所傳播的距離為

式中κ為介質(zhì)的吸收系數(shù);Rs∈[0,1],是一個均勻分布的隨機(jī)數(shù)。

3)光線在傳播過程中如果發(fā)生散射,散射角Θ和i1通過式(11)確定,并假設(shè)所有生成散射方向都是有效的,繼續(xù)跟蹤散射后的光線。

4)圖3中光線(黑色實(shí)線)在傳播過程中遇到邊界時,不考慮光線折射透過介質(zhì),即假設(shè)上下界面為鏡面反射,繼續(xù)跟蹤直到光線被介質(zhì)所吸收。

圖3 逆正向蒙特卡羅法光線跟蹤軌跡示意圖Fig.3 Ray tracing track diagram of reverse forward Monte Carlo method

以上四個步驟是光線的“特殊”反向跟蹤過程,下面是正向Stokes矢量計(jì)算過程。

5)在逆向光線跟蹤完成后,根據(jù)光路的可逆性,光線可以沿原來路徑返回。光線的Stokes矢量計(jì)算時使用之前已存儲的數(shù)據(jù),而無需重復(fù)跟蹤光線。

6)在正向計(jì)算的光線到達(dá)界面處,使用矢量輻射傳輸問題的反射率公式(16)判斷光線是透射還是反射,如果在表面透射則停止跟蹤。同樣,當(dāng)光線到達(dá)散射的位置時,根據(jù)拒絕法來判斷之前存儲的散射方向是否有效,如果無效則停止這束光線的跟蹤,并記錄被拒絕的總次數(shù)。最后在模擬光線的總數(shù)量中減去拒絕法所拒絕的總次數(shù),這相當(dāng)于在標(biāo)準(zhǔn)的拒絕法中重新生成隨機(jī)數(shù)的過程。

7)最終可以沿著原來入射方向射出介質(zhì)光線的比例就是表觀發(fā)射率,表觀Stokes矢量S(θ,φ)公式為

式中i為圖1中的i子層;K為介質(zhì)子層數(shù);Ib(Ti)是i子層溫度的黑體輻射強(qiáng)度;Ni→s1為從i子層發(fā)射的光線可以沿著介質(zhì)表面的(θ,φ)方向發(fā)射出去的總數(shù);Sm(θ,φ)是第m根光線歸一化的Stokes矢量;N為光線追蹤的總數(shù);Nrject為各向異性散射介質(zhì)中被拒絕法拒絕的總次數(shù),對于非散射問題和各向同性散射類問題,Nrject=0。

2 一維梯度折射率介質(zhì)矢量輻射傳輸

2.1 非散射性介質(zhì)

設(shè)模型的介質(zhì)折射率分布為n(z)=n0+(nL-n0)z/L,z為折射率變化方向,上下兩個邊界處的折射率分別為n0和nL,n0=1.2,nL=1.8,厚度L=0.1。該模型是一個無限大平行吸收發(fā)射的平板,為非散射灰色半透明介質(zhì),其邊界被認(rèn)為是純菲涅耳反射。介質(zhì)的吸收系數(shù)κa=10m-1,光線的初始Stokes矢量為S0=(1,0,0,0)T。研究的重點(diǎn)是介質(zhì)上下表面某點(diǎn)的Stokes矢量,并與文獻(xiàn)[28]中的標(biāo)量方法結(jié)果進(jìn)行了對比。

圖4是采用逆正向蒙特卡羅法計(jì)算上表面處I、Q兩個分量的結(jié)果,而U和V分量等于零沒有在圖中繪制,其計(jì)算條件為n(z)=1.2+0.6z/L,L=0.1,κa=10m-1,S0=(1,0,0,0)T。反向跟蹤模擬光線的總數(shù)量為106,當(dāng)模型被劃分為100層和200層時,計(jì)算結(jié)果十分接近,而10、20、50層之間的偏差較大。為了保證其結(jié)果的精度,本研究采用將介質(zhì)劃分為200個子層的近似方案進(jìn)行計(jì)算。

圖4 不同子層數(shù)量時上表面某點(diǎn)的Stokes矢量。(a)I分量;(b)Q分量Fig.4 Stokes vector of a point on the upper surface with different sublayers.(a)I component;(b)Q component

εp和εs分別是表觀定向平行發(fā)射率和表觀定向垂直發(fā)射率,其與Stokes矢量的換算關(guān)系為

式中Ip是P極化波的輻射強(qiáng)度,Ib是黑體輻射強(qiáng)度,Is是S極化波的輻射強(qiáng)度。

圖5是采用逆正向蒙特卡羅法所計(jì)算的等溫條件下介質(zhì)上表面某點(diǎn)的方向發(fā)射率與文獻(xiàn)[28]結(jié)果的對比,計(jì)算條件為n(z)=1.2+0.6z/L,L=0.1,κa=10m-1,S0=(1,0,0,0)T。采用逆正向蒙特卡羅法所得到的結(jié)果和文獻(xiàn)[28]標(biāo)量方法結(jié)果一致,其表觀定向平行發(fā)射率εp相對于表觀定向垂直發(fā)射率εs要大一些。對于均勻分布的等溫介質(zhì),其上下表面的發(fā)射特性一致。而對于梯度折射率分布介質(zhì),介質(zhì)上下表面的發(fā)射率則存在明顯的區(qū)別,梯度折射率介質(zhì)與均勻折射率介質(zhì)表現(xiàn)出不一樣的輻射特性。

圖5 介質(zhì)上下表面某點(diǎn)的發(fā)射率對比結(jié)果。(a)I分量;(b)Q分量Fig.5 Emissivity comparison of a point on the upper and lower surface of a medium.(a)I component;(b)Q component

圖6 為不同光學(xué)厚度介質(zhì)τ在上表面的Stokes矢量,計(jì)算條件為n(z)=1.2+0.6z/L,L=0.1,κa=10m-1,S0=(1,0,0,0)T,200子層,模擬光線的總數(shù)量為106。隨著光學(xué)厚度的減小,在大多數(shù)角度中更多的光線會直接穿過介質(zhì),使得I和Q部分的分量隨著光學(xué)厚度變小而顯著減少。但當(dāng)光學(xué)厚度在0.7、角度在75°以上時,I和Q的分量沒有減小。此時光線厚度雖然有所減少,但是在入射角度較大時,界面反射概率增加和全發(fā)射的發(fā)生使得光線繼續(xù)在介質(zhì)中傳播,并產(chǎn)生線性偏振。

圖6 不同光學(xué)厚度介質(zhì)上表面的Stokes矢量。(a)I分量;(b)Q分量Fig.6 Stokes vectors of upper surfaces on media with different optical thickness.(a)I component;(b)Q component

2.2 各向異性散射性介質(zhì)

當(dāng)介質(zhì)散射類型為瑞利散射時,散射光的強(qiáng)度和入射光波長的四次方成反比,其散射相矩陣M(Θ)為

當(dāng)散射反照率ω=0.5時,采用逆正向蒙特卡羅法得到的介質(zhì)上表面的Stokes矢量如圖7所示,其計(jì)算條件為n(z)=1.2+0.6z/L,L=0.1,κa=10m-1,S0=(1,0,0,0)T,200子層,模擬光線的總數(shù)量為106。在圖7(a)、(b)中I和Q分量數(shù)值相對較大,而在圖7(c)、(d)中U和V分量非常小,接近0。與圖5相比,散射對Stokes矢量的影響是明顯的,光線發(fā)生了散射,這使得I和Q分量都有所減少。

圖7 反照率為ω=0.5時介質(zhì)上表面Stokes矢量。(a)I分量;(b)Q分量;(c)U分量;(d)V分量Fig.7 Stokes vector of upper surface on medium when ω=0.5.(a)I component;(b)Q component;(c)U component;(d)V component

圖8是在散射反照率ω=0.5時,采用逆正向蒙特卡羅法和標(biāo)量反向蒙特卡羅法求得的上表面和下表面某處發(fā)射率的對比,圖8的計(jì)算條件為n(z)=1.2+0.6z/L,L=0.1,κa=10m-1,S0=(1,0,0,0)T,200子層,模擬光線的總數(shù)量為106,此時兩種方法的計(jì)算結(jié)果并不完全一致。在散射的條件下,介質(zhì)中光線被散射的概率增加了,考慮偏振狀態(tài)時光線散射方向的分布概率與標(biāo)量不同,因此逆正向蒙特卡羅法發(fā)射率的結(jié)果與標(biāo)量方法也存在一定的差異。

圖8 介質(zhì)上下表面某點(diǎn)的發(fā)射率結(jié)果。(a)上表面;(b)下表面Fig.8 Emissivity results for a point on the upper and lower surfaces of a medium.(a)Upper surface;(b)lower surface

3 結(jié)論

采用逆正向蒙特卡羅法求解了考慮光線偏振情況時一維梯度折射率介質(zhì)的表觀發(fā)射特性。分別研究了一維線性折射率分布非散射性介質(zhì)和各向異性散射介質(zhì)的表觀矢量,得到了Stokes矢量圖像,并與已有文獻(xiàn)中計(jì)算的標(biāo)量表觀發(fā)射率結(jié)果進(jìn)行了對比驗(yàn)證。主要內(nèi)容和結(jié)論有:

1)傳統(tǒng)的蒙特卡羅法各個子過程之間相互獨(dú)立且連貫,一個子過程計(jì)算結(jié)束后才進(jìn)入下一個子過程。例如在界面反射的子過程中,生成了隨機(jī)數(shù)就立即去判斷光線是反射還是透射。而逆正向蒙特卡羅法中,并沒有把光線軌跡上事件保持在一個連續(xù)的計(jì)算過程。各種事件是否發(fā)生的判斷過程,可以發(fā)生在其他子過程中。這使得蒙特卡羅法的應(yīng)用更為靈活,為解決其他更為復(fù)雜問題提供了新的思路。

2)提出了逆正向蒙特卡羅法,實(shí)現(xiàn)Stokes矢量反向跟蹤正向計(jì)算。通過與標(biāo)量方法的對比和分析,逆正向蒙特卡羅法可以很好解決此類復(fù)雜問題,其結(jié)果表明考慮光線偏振的表觀發(fā)射特性與標(biāo)量方法的結(jié)果存在一定的區(qū)別。

3)由于介質(zhì)中梯度折射率的分布,使得介質(zhì)上下表面的發(fā)射率存在了差別。隨著光學(xué)厚度減小,介質(zhì)表觀Stokes矢量也隨之減小,但在較大觀測角度呈現(xiàn)復(fù)雜結(jié)果。在一維梯度折射率介質(zhì)表面的Stokes矢量圖像中,I和Q分量數(shù)值較大,而U和V分量非常小。