喬文佑，高志山，袁群，朱丹，許寧晏，倫旭磊，車嘯宇

（南京理工大學(xué) 電子工程與光電技術(shù)學(xué)院，南京 210094）

0 引言

在空間相機、衛(wèi)星激光通信、引力波探測等高精度空間平臺光機系統(tǒng)中，雜散光對系統(tǒng)信號的影響是不容忽視的［1-3］。光機系統(tǒng)的元器件是雜散光的重要來源，其表面光散射特性直接影響著系統(tǒng)內(nèi)雜散光分布情況，針對系統(tǒng)元器件表面散射特性進行光線追跡和雜散光分析，是光學(xué)設(shè)計及仿真中重要的研究內(nèi)容。以蒙特卡洛法（Monte Carlo Method，MCM）模擬抽樣為基礎(chǔ)，基于描述表面散射特性的雙向反射分布函數(shù)（Bidirectional Reflectance Distribution Function，BRDF）來進行散射光線追跡，是目前科研實驗及雜散光分析軟件的常用算法。

應(yīng)用MCM 實現(xiàn)散射光線準(zhǔn)確追跡的關(guān)鍵在于由BRDF 模型建立正確的散射光線追跡概率模型，常用的概率模型求解方案主要為反變換法［4］。以反變換法設(shè)計概率模型進行光線追跡的方法，主要核心內(nèi)容：1）選用現(xiàn)有的BRDF 模型解析式對實測BRDF 數(shù)據(jù)進行擬合解析，選出合理的BRDF 模型；2）將BRDF 轉(zhuǎn)化為概率密度函數(shù)（Probability Density Function，PDF），通過求解PDF 的累積分布函數(shù)（Cumulative Distribution Function，CDF）的反函數(shù)，完成合理采樣并以此確定散射光線方向，最終實現(xiàn)光線追跡。然而，系統(tǒng)元器件表面由于其微觀形貌特點和紋理分布的特殊性，表面散射特性顯示出復(fù)雜多樣化，且伴隨加工技術(shù)創(chuàng)新及新材料表面的出現(xiàn)，商用軟件包內(nèi)建的或現(xiàn)有的BRDF 模型或相關(guān)修正模型，可用于精確表征表面散射屬性的模型種類有限。BRDF 模型參數(shù)解析過程繁雜，存在擬合誤差，模型適用條件存在限制等問題［5-6］。同時，反變換法雖然可高效產(chǎn)生服從指定分布的隨機樣本，但是多數(shù)散射模型的BRDF 受光線坐標(biāo)變量θ與φ調(diào)制，在復(fù)雜概率模型設(shè)計過程中，對形式復(fù)雜的PDF，其CDF 存在無解析解的情況，反變換法失效［7］。

針對以上問題，為簡化表面散射特性的建模分析過程，增強散射光線追跡方法的適用性，本文提出了一種直接使用散射表面BRDF 離散數(shù)值的散射光線追跡方法。對隨散射角離散變化的各向同性散射表面BRDF 離散測量數(shù)據(jù)組，由散射角半球空間轉(zhuǎn)換到方向余弦空間，通過等間距賦值插等值方式，得到方向余弦空間內(nèi)完整分布的BRDF 數(shù)據(jù)。利用舍選抽樣法不受限于CDF 求解過程的優(yōu)點，設(shè)計新的散射概率模型，即計算方向余弦空間內(nèi)BRDF 數(shù)值比表示離散光線的概率分布，設(shè)定檢驗條件篩選出散射光線空間坐標(biāo)，實現(xiàn)散射光線追跡。

1 理論基礎(chǔ)

1.1 雙向反射分布函數(shù)表征模型

雙向散射分布函數(shù)（Bidirectional Scattering Distribution Function，BSDF）是描述材料表面對入射光在半球空間內(nèi)散射情況的度量函數(shù)［8］，表征結(jié)構(gòu)表面散射特性，被廣泛應(yīng)用于光學(xué)遙感等測量領(lǐng)域［2，9］。BSDF 實際包括反射和透射兩個部分，如對透射部分的散射光感興趣，需要定義雙向透射分布函數(shù)（Bidirectional Transmittance Distribution Function，BTDF）。因篇幅有限，不失一般性，以表面散射光線追跡算法為例，討論與反射部分BRDF 有關(guān)的散射光線追跡方法，方法同樣適用于透射情形。BRDF 被定義為空間某一散射方向光線的輻射亮度與入射光照在材料表面上的輻射照度比值，即

式中，光線入射角由θi與φi兩個分量表示，分別為入射光線的天頂角與方位角；散射光線的散射角也由θs與φs兩個分量表示，分別為散射光線天頂角與方位角；dLs(θi，φi，θs，φs)為散射方向的輻射亮度；dEi(θi，φi)為入射方向的輻射照度，物理模型如圖1 所示。

圖1 BRDF 物理模型Fig.1 BRDF physical model

自然界中，具有各向同性散射特征的表面是普遍存在的，工程應(yīng)用中，加工工藝決定了光機系統(tǒng)的很多介質(zhì)表面散射可近似看作是各向同性的，對此類介質(zhì)表面散射特性表征的經(jīng)典BRDF 模型有ABg 模型［5］、Harvey 模型［10］等。該類模型解析式中，可將變量(θi，φi，θs，φs)轉(zhuǎn)換成變量為R的函數(shù)，

式中，ri與rs分別為入射光線與散射光線在方向余弦x-y坐標(biāo)系中投影單位向量，如圖1 所示。xi和yi與xs和ys分別為入射光與散射光在方向余弦坐標(biāo)系中的坐標(biāo)，其中xi=sinθicosφi，yi=sinθisinφi，xs=sinθscosφs，ys=sinθssinφs，即為半球空間內(nèi)光線的天頂角與方位角在方向余弦空間的變量關(guān)系。在方向余弦空間（x-y坐標(biāo)系）中，光線映射成點，與鏡像反射光線映射點等間距位置的BRDF 數(shù)值始終相等?？筛鶕?jù)入射光線和鏡像反射光線映射點所在平面（如圖2 方框切面xs=0）內(nèi)的BRDF 數(shù)據(jù)，通過等間距賦值插等值方式，得到散射表面在方向余弦空間內(nèi)完整分布的BRDF 數(shù)據(jù)，即散射角半球空間內(nèi)任意散射角度光線的BRDF 值。

圖2 方向余弦空間介質(zhì)表面BRDF 數(shù)值分布Fig.2 Cosine-space BRDF distribution on medium surface

1.2 基于蒙特卡洛-舍選抽樣法的光線追跡模型

蒙特卡洛法（MCM）是一種隨機統(tǒng)計的概率模擬方法，其思想是通過對隨機變量進行統(tǒng)計試驗，進而求得近似解的數(shù)值方法［11］。在雜光分析仿真中，MCM 是非序列光線追跡的基礎(chǔ)算法，對光線在系統(tǒng)傳播過程出現(xiàn)的散射等問題，通常采取構(gòu)建散射光線分布的概率模型，將光線的空間坐標(biāo)作為抽樣對象，經(jīng)過模擬抽樣確定散射光線方向，實現(xiàn)散射能量在空間分布，轉(zhuǎn)變?yōu)閿y帶等能量的離散光線分布，以散射空間光線密度表征能量分布，最終完成散射光線追跡［12］。

如前所述，基于MCM 實現(xiàn)散射光線追跡的關(guān)鍵，在于對其概率模型的正確構(gòu)建，而概率模型求解方法常使用反變換法，但反變換法時常失效［7，13］。對于已知光機表面BRDF 離散實測數(shù)據(jù)的情形，可以舍棄由實測數(shù)據(jù)經(jīng)模型擬合、PDF 和CDF 獲取、求解CDF 的反函數(shù)等復(fù)雜繁瑣的技術(shù)途徑，引入舍選抽樣法，直接應(yīng)用BRDF 離散實測數(shù)據(jù)組，得到用于散射光線追跡的光線抽樣概率函數(shù)。

舍選抽樣法也稱為拒絕采樣法［14］，旨在由簡單樣本完成對復(fù)雜樣本的抽樣。以一維函數(shù)為例，設(shè)由BRDF 實測數(shù)據(jù)獲取的PDF 函數(shù)記為P(x)，通過建立檢驗條件，對產(chǎn)生的隨機數(shù)樣本完成服從P(x)分布的舍選抽樣，具有采樣靈活、操作方便等優(yōu)點。舍選抽樣法既可以用于連續(xù)分布對象，也可以用于離散分布對象，以舍選抽樣法設(shè)計概率模型時對多維等復(fù)雜數(shù)學(xué)分布對象同樣適用。如圖3，舍選抽樣時，會選取密度函數(shù)c?M(x)，c?M(x)也被稱為擠壓函數(shù)，使得c?M(x)能夠覆蓋待抽樣對象的密度函數(shù) ，其中c為拒絕常數(shù)且通常取c=max {M(x)/P(x)}。為保證高效的抽樣效率，選取的c?M(x)應(yīng)盡量貼近P(x)曲線且表達式簡單。

圖3 舍選抽樣模型Fig.3 Rejection sampling model

舍選抽樣過程中，首先生成密度函數(shù)為c?M(x)的一個偽隨機數(shù)x0，再從均勻分布（0，1）中產(chǎn)生一個偽隨機數(shù)u，當(dāng)u滿足關(guān)系式（3）檢驗條件時，x0為接受樣本，否則拒絕，重新生成隨機數(shù)。重復(fù)以上過程，即可得到許多概率密度為P(x)的偽隨機數(shù)，最終完成對樣本x的抽樣。

舍選抽樣算法的抽樣效率η，實則為c?M(x)的樣本被接受成為P(x)的樣本的概率，即有

可見，抽樣效率由選取的函數(shù)c?M(x)與P(x)共同決定，兩者越貼近，抽樣效率越大。若c?M(x)函數(shù)在極小值處，存在函數(shù)值與P(x)最貼近，抽樣效率最大，若c?M(x)函數(shù)在極大值處，存在函數(shù)值與P(x)最貼近，抽樣效率最大。當(dāng)函數(shù)c?M(x)越接近P(x)，舍選抽樣的效率越高，但相應(yīng)的抽樣難度也越來越大。由于P(x)形式復(fù)雜多樣，若一一選取貼合度高的擠壓函數(shù)，不僅選取困難且對其他形式分布的抽樣對象存在適用性問題。為了降低抽樣難度及提升通用性，對函數(shù)c?M(x)取常數(shù)Z處理，一般Z值可取P(x)的最大值，即

對于常數(shù)Z的處理方式，雖然在一定程度上降低了抽樣效率，但可以做到對任意分布形式的P(x)函數(shù)，都可進行模擬抽樣，同時規(guī)避了反變換法過程中積分求解難題，通用性好。

接下來，需要解決的問題是，如何從半球空間內(nèi)任意散射角度光線的BRDF 值（由BRDF 離散測量值得到）產(chǎn)生密度函數(shù)P(x)。從概率論角度看，密度函數(shù)P(x)的本質(zhì)，是描述隨機變量的輸出值在某個確定的取值點附近的可能性。因此，最簡捷的產(chǎn)生方法，是將方向余弦空間內(nèi)的個體BRDF 離散數(shù)值與總體BRDF 離散數(shù)值總和的比值，作為表征各抽樣點相應(yīng)密度函數(shù)值P(x)。對于二維情況，由BRDF 實測數(shù)據(jù)計算的P(xi，yj)表達式如式（6）。

式中，設(shè)定方向余弦空間BRDF 的變量(xi，yj)取隨機數(shù)，即為待舍選的抽樣樣本，且需要滿足x2i+yj2≤1。 此時，Z=max(P(xi，yj))，式（3）改寫為u?max(P(xi，yj))≤P(xi，yj)，則在舍選抽樣過程中，檢驗條件“u?max(P(xi，yj))≤P(x?i，yj?)且(x?i)2+(yj?)2≤1”對樣本(x?i，yj?) 進行舍選，完成服從P(xi，yj)分布的抽樣，以此來確定散射光線方向，為直接利用BRDF 離散數(shù)值的散射光線追跡方法奠定基礎(chǔ)，算法流程如圖4 所示。

圖4 舍選抽樣算法流程Fig.4 Flow chart of rejection sampling method

2 仿真驗證

2.1 仿真設(shè)置

為驗證基于舍選抽樣BRDF 離散數(shù)值的散射光線追跡方法的準(zhǔn)確性，設(shè)置相同的入射角、追跡光線數(shù)量等仿真參數(shù)，利用Matlab編寫了本文光線追跡仿真程序，并選擇與商用軟件LightTools相同的概率分光原理設(shè)計概率模型［7］，獲取散射光線追跡結(jié)果的數(shù)據(jù)群，與相同條件的LightTools軟件散射光線計算結(jié)果進行比對。

不失一般性，為了方便與LightTools 軟件結(jié)果的比對，光機表面BRDF 的離散實測值，應(yīng)選擇LightTools 軟件中具有的內(nèi)建光機表面模型，在相同入射角、追跡光線數(shù)量情況下，通過對散射角離散獲得BRDF 數(shù)值，即獲取其半球空間中入射光線與鏡像反射光線所在平面內(nèi)的BRDF 離散測量數(shù)據(jù)，作為本文方法Matlab 追跡程序中散射表面BRDF 實測值。

Harvey 模型是光機表面BRDF 建模的常用散射模型，常用于表征粗糙度在納米量級的光學(xué)表面散射特性，式（7）為Harvey 模型BRDF 表達式。在LightTools 中設(shè)置Harvey 散射表面屬性，可以表征不同光機表面的模型參數(shù)，構(gòu)建光線追跡仿真系統(tǒng)。

仿真中，對于反射鏡表面，Harvey 模型參數(shù)為b0=2.43，L=1.8×10?2，s=?2.32［15］；以及Harvey 模型參數(shù)為b0=1.185×10?2，L=3.10×10?2，s=?1.9［16］的光學(xué)面，仿真設(shè)置入射光線能量為1 W，追跡光線數(shù)為107，散射能量分析面為邊長0.5 mm 的正方形，界面內(nèi)網(wǎng)格數(shù)為100×100，圖5 為LightTools 中散射幾何空間仿真示意圖，表1 為仿真元件的空間坐標(biāo)。

表1 仿真系統(tǒng)坐標(biāo)Table 1 Simulation system coordinates

圖5 散射仿真系統(tǒng)Fig.5 Scattering simulation system

2.2 驗證分析

基于Harvey 模型，按照2.1 所述參數(shù)設(shè)置方法，構(gòu)建散射光線追跡仿真系統(tǒng)，仿真結(jié)果能量分布如圖6所示，其中（a）為本文方法Matlab 追跡程序的計算結(jié)果，（b）為LightTools 軟件的計算結(jié)果。

圖6 Harvey 仿真模型散射能量分布圖，b0=1.185×10?2，L=3.10×10?2，s=?1.9Fig.6 Scattering energy distribution of Harvey simulation model， b0=1.185×10?2，L=3.10×10?2，s=?1.9

由圖6 可以看出，仿真結(jié)果是與追跡光線數(shù)量匹配的數(shù)據(jù)群，數(shù)據(jù)量較大。為對比仿真結(jié)果的偏差，本文引入了通用質(zhì)量指數(shù)（Universal Quality Index， UQI）。UQI 是通用的圖像質(zhì)量多元素同步評價指標(biāo)，其定義包含了相關(guān)性損失、亮度失真和對比度失真三個元素［17］，可以量化圖像之間的相對失真程度。因為成像本質(zhì)與能量分布相關(guān)，因此，本文采用UQI 作為仿真結(jié)果相似性的評價值，通過計算UQI 來量化散射能量分布的相似度。UQI 表達式為

由于蒙特卡洛模擬過程中采樣具有隨機性，造成連續(xù)散射光線追跡仿真結(jié)果并不完全相同，為驗證仿真的精確度，對上述Harvey 模型及參數(shù)進行了多次同步仿真。通過觀察Matlab 追跡程序第n次與LightTools 第n次仿真散射能量之間的UQI 值動態(tài)范圍，驗證本文光線追跡方法的合理性。

圖7 為不同Harvey 散射模型參數(shù)表征的材料表面散射光線追跡仿真結(jié)果，對每種參數(shù)模型分別進行了10 次關(guān)于本文光線追跡方法與LightTools 同步仿真對比?？梢?，本文光線追跡方法與LightTools 同步仿真結(jié)果的UQI 值動態(tài)范圍分別在0.999 4～0.999 5 和0.999 2～0.999 3 區(qū)間，且UQI 值的波動范圍小，說明本文光線追跡方法與LightTools 仿真得到的散射能量分布高度一致，因此，本文基于舍選抽樣BRDF 離散數(shù)值的散射光線追跡方法具有較高地準(zhǔn)確度與穩(wěn)定度。對于不同材料表面散射仿真結(jié)果評價的UQI 值出現(xiàn)的差異，是由BRDF 模型種類及其模型參數(shù)決定。

圖7 Harvey 模型仿真UQI 動態(tài)范圍Fig.7 The UQI range of Harvey model simulation

為了進一步驗證本文基于舍選抽樣BRDF 離散數(shù)值的散射光線追跡方法的適用性，同上述仿真對比過程，選擇能夠?qū)Υ罅抗鈾C結(jié)構(gòu)表面建模的ABg 模型。其中，選擇了用于表征材質(zhì)為BK7 平凸型準(zhǔn)直物鏡［18］表面散射特性和發(fā)黑鋁合金［5］表面散射特性的ABg 模型參數(shù)，進行多次同步光線追跡仿真實驗，模型參數(shù)見表2。

表2 仿真系統(tǒng)坐標(biāo)Table 2 Simulation system coordinates

不同參數(shù)的ABg 模型所表征的材料表面散射光線追跡仿真結(jié)果如圖8 所示?？梢?，本文光線追跡方法與LightTools 同步仿真結(jié)果的UQI 值動態(tài)范圍分別在0.999 0～0.999 1 和0.998 7～0.998 8 區(qū)間，UQI 值的波動范圍小，說明兩者仿真所得到的散射能量分布高度一致，本文基于舍選抽樣BRDF 離散數(shù)值的散射光線追跡方法不僅準(zhǔn)確度高，對不同表面散射特性的表征同樣具有良好地適用性。

圖8 ABg 模型仿真UQI 動態(tài)范圍Fig.8 The UQI range of ABg model simulation

2.3 復(fù)合散射表面應(yīng)用

由于系統(tǒng)元器件或涂層表面微觀形貌特點和紋理分布的特殊性，其表面散射呈現(xiàn)復(fù)雜多樣化，單一BRDF 解析模型很難對測定的散射數(shù)據(jù)進行高精度擬合，通常需要多個BRDF 解析模型進行復(fù)合擬合，并確定各BRDF 模型的權(quán)重系數(shù)，復(fù)合BRDF 模型如式（9）所示。

式中，BRDF1、BRDF2分別為各種類BRDF 模型，k1、k2分別為各類BRDF 模型的權(quán)重系數(shù)。

為驗證本文基于舍選抽樣BRDF 離散數(shù)值的散射光線追跡方法對復(fù)合散射表面光線追跡的適用性，在LightTools 軟件中構(gòu)建了包含朗伯散射與高斯散射的復(fù)合散射表面光線追跡模型。其中，朗伯散射模型表達式為

式中，ρ為總反射比。

高斯散射模型表達式為

式中，K為歸一化常數(shù)，，σx與σy為高斯分布的標(biāo)準(zhǔn)差。

LightTools 仿真模型中設(shè)置朗伯散射與高斯散射的權(quán)重系數(shù)k1與k2分別為0.3 與0.7，朗伯散射的總反射比ρ為1，高斯散射分布的標(biāo)準(zhǔn)差σx與σy均為15°。通過以上復(fù)合BRDF 模型及參數(shù)，獲取其半球空間中入射光線與鏡像反射光線所在平面內(nèi)的BRDF 離散測量數(shù)據(jù)，作為本文方法Matlab 追跡程序中散射表面BRDF 實測值，利用Matlab 追跡程序與LightTools 軟件進行多次同步仿真對比。

圖9 為復(fù)合散射表面基于本文光線追跡方法與LightTools 軟件建模進行10 次同步仿真對比結(jié)果?？梢?，本文光線追跡方法與LightTools 同步仿真結(jié)果的UQI 值動態(tài)范圍在0.998 8～0.998 9 區(qū)間，UQI 值的波動范圍小，說明兩者仿真得到的散射能量分布高度一致。因此，本文基于舍選抽樣BRDF 離散數(shù)值的散射光線追跡方法對于復(fù)合散射表面的光線追跡同樣具有良好的適用性。

圖9 復(fù)合散射表面模型仿真UQI 動態(tài)范圍Fig.9 The UQI range of ABg model simulation

相較于LightTools 軟件無法直接構(gòu)建介質(zhì)表面為Harvey、ABg、高斯等復(fù)雜的復(fù)合散射表面仿真模型，本文提出的基于舍選抽樣BRDF 離散數(shù)值的散射光線追跡方法具有良好的光學(xué)仿真優(yōu)勢，可以準(zhǔn)確、便捷地完成此類復(fù)雜的復(fù)合散射表面光線追跡。

3 結(jié)論

為簡化表面散射特性建模分析過程，增強散射光線追跡方法的適用性，本文提出了一種直接使用散射表面BRDF 離散數(shù)值的散射光線追跡方法?；诿商乜宸M散射光線分布思想，以方向余弦空間內(nèi)BRDF 數(shù)值比表示離散光線的概率分布，結(jié)合舍選抽樣法篩選出散射光線的空間坐標(biāo)，實現(xiàn)散射光線追跡。該方法不僅避免了對BRDF 模型選用的繁雜過程，同時也規(guī)避了反變換法過程中積分求解難題。為了驗證本文光線追跡方法的準(zhǔn)確性與適用性，在相同入射角、追跡光線數(shù)量情況下，選用光學(xué)軟件LightTools 設(shè)置表征各類光機元器件表面散射特性的Harvey、ABg 以及復(fù)合散射表面的BRDF 模型及參數(shù)，進行建模仿真。利用Matlab 編制本文光線追跡方法的仿真程序，通過LightTools 設(shè)置的BRDF 模型及參數(shù)，獲取其半球空間中入射光線與鏡像反射光線所在平面內(nèi)的BRDF 離散測量數(shù)據(jù)，作為本文方法Matlab 追跡程序中散射表面BRDF 離散實測值，進行同步光線追跡仿真。采用UQI 作為LightTools 軟件仿真與本文光線追跡方法仿真結(jié)果比對的評價值。結(jié)果顯示，UQI 值均在0.998 以上且波動范圍小，說明本文基于舍選抽樣BRDF 離散數(shù)值的散射光線追跡方法準(zhǔn)確度高且具有良好地適用性。隨著加工制造技術(shù)和光學(xué)涂層工藝的不斷更新，為滿足精密光學(xué)系統(tǒng)研制中對新型機械加工材料和涂層表面散射特性的分析工作，本文散射表面光線追跡方法可為此項工作提供一種簡捷準(zhǔn)確地分析途徑與思路，具有重要的實際應(yīng)用和推廣價值。