基于蒙特卡羅方法的水體光場傳輸特性研究?

（海軍駐武漢地區(qū)第七軍事代表室 武漢 430223）

1 引言

進(jìn)入21世紀(jì)，隨著經(jīng)濟(jì)和技術(shù)的快速發(fā)展，對(duì)海洋的研究已經(jīng)成為一個(gè)新的領(lǐng)域。海洋資源開發(fā)、海洋自然環(huán)境保護(hù)以及海洋軍事應(yīng)用等方面的研究越來越引起人們的關(guān)注［1］。

光電探測成像技術(shù)是海洋探測及環(huán)境保護(hù)的重要手段之一。然而，光在水下傳播過程中，受到水分子和其他懸浮粒子的作用，將發(fā)生吸收和散射現(xiàn)象，使得光線能量逐漸衰減，且成像質(zhì)量急劇下降。因此，對(duì)水下光場傳輸特性的研究具有重要的意義［7］。

研究水下光場傳播規(guī)律的主要方法有幾何光學(xué)方法、輻射傳輸方程方法［3］和蒙特卡羅模擬方法等［4～5］。幾何光學(xué)方法主要是采用光在水空界面的折射定律和光在純凈水中的直線傳播規(guī)律，無法深入研究散射背景條件下的光傳播特性。輻射傳輸方程雖然可以描述光在強(qiáng)散射水體中的傳播特性，可詳細(xì)地考慮散射、吸收、輻射等光與物質(zhì)相互作用的現(xiàn)象，但是對(duì)于復(fù)雜邊界條件和初始條件下方程求解比較困難。而蒙特卡羅模擬方法［6～9］作為研究散射介質(zhì)中光傳輸特性的一種適用廣泛而嚴(yán)謹(jǐn)?shù)姆椒?。本文基于該方法仿真分析了在?qiáng)散射背景條件下水體的光場傳輸規(guī)律。

2 方法

光子在水體中的傳播規(guī)則采用概率分布來描述。如圖1是蒙特卡羅方法水體散射仿真流程圖，單個(gè)光子首先被初始化統(tǒng)一權(quán)重，設(shè)置光子步長，并在坐標(biāo)系中移動(dòng)它。

圖1 蒙特卡羅方法水體散射仿真流程圖

若光子到達(dá)水體邊緣，檢查其內(nèi)部反射的幾率，并決定其是繼續(xù)留在水體內(nèi)還是離開水體。若其在內(nèi)部反射，則需修改其運(yùn)動(dòng)方向，繼續(xù)運(yùn)動(dòng)；否則，光子離開水體，記錄為反射（透射）。在每運(yùn)動(dòng)一個(gè)步長，計(jì)算光子由于水體的吸收造成的權(quán)重衰減，并記錄在該處水體吸收矩陣中。剩下的光子權(quán)重代入到散射后移動(dòng)方向傳輸步長計(jì)算中。若光子權(quán)重低于一定小的閾值，則執(zhí)行輪盤賭算法來決定忽略此光子或繼續(xù)計(jì)算它的移動(dòng)。

2.1 光子的初始化

每一個(gè)光子在最初統(tǒng)一設(shè)定一個(gè)權(quán)重W。當(dāng)光子在水體邊界遇到折射率改變界面時(shí)，會(huì)發(fā)生鏡面反射。設(shè)入射空氣折射率為n1，水體折射率為n2，則鏡面反射率 Rsp為

則光子的權(quán)重被減去Rsp。

2.2 光子步長

根據(jù)吸收系數(shù)和散射系數(shù)的定義，光子與水體介質(zhì)在單位步長［s1，s1+ds1］內(nèi)相互作用（包括吸收與散射）的幾率為μ1ds1。由水體的比耳朗伯定律［4］可得:

其中，μa為吸收系數(shù)，μs為散射系數(shù)，μt是總衰減系數(shù)，單位是m-1。對(duì)上式進(jìn)行變換，在［0，s1］內(nèi)積分，得到指數(shù)分布:

表示光子步長大于 s1的幾率，即P｛s≥s1｝。對(duì)于自由路徑s的累積分布:

變換上式求解得到:

設(shè)定ζ=P｛s＜s1｝，且在（0，1］區(qū)間均勻分布，因此ζ與1-ζ的分布相同，則-ln（1-ζ）與-ln（ζ）等價(jià)，故:

依此式，在（0，1］區(qū)間內(nèi)生成隨機(jī)數(shù)ζ，便得到步長s=s1。

2.3 光子移動(dòng)

當(dāng)s被確定以后，光子就可以在坐標(biāo)系中移動(dòng)。設(shè)光子的瞬時(shí)位置為（x，y，z），瞬時(shí)方向單位向量r，可表示為（μx，μy，μz）:

在移動(dòng)后光子的位置為（x’，y’，z’），則:

2.4 邊界處理

在光子移動(dòng)一個(gè)步長過程中，光子可能會(huì)遇到水體邊界。當(dāng)光子遇到邊界時(shí)，光子可能反射（或透射）逃逸出水體或者向內(nèi)部反射。光子向內(nèi)反射的幾率取決于入射到邊界的角度（θi）?？捎?jì)算得到:

由Snell公式表示入射角θi與折射角θt的關(guān)系。向內(nèi)反射的幾率R（θi），可由Fresnel公式得出:

因此，1-R（θi）表示光子離開水體成為反射部分，計(jì)入反射矩陣。

在光子離開介質(zhì)前，只行走了步長s的一部分。因此，實(shí)際離開位置應(yīng)該是基于減短的步長s':

其中τ表示水體的厚度。若光子內(nèi)反射，則光子權(quán)重計(jì)算為

剩下的光子有新的位置和方向。x，y坐標(biāo)可直接用前述公式計(jì)算，z坐標(biāo)更新為

2.5 光子吸收

當(dāng)光子移動(dòng)一定步長時(shí)，光子由于水體介質(zhì)的吸收，會(huì)使得權(quán)重減小。權(quán)重減小的幅度ΔQ:

光子權(quán)重W更新:

由于（μa/μt+μs/μt）=1，因此能量是平衡的。

2.6 光子散射

當(dāng)光子發(fā)生散射時(shí)，光子的軌道偏轉(zhuǎn)角度θ分布在［0，π］。Henyey與Greenstein［10］最早提出利用散射相函數(shù)來計(jì)算散射偏轉(zhuǎn)角的概率密度方程:

其中，g為各向異性因子，表示散射的角度分布。令μ=cosθ，分布在［-1，1］?？梢酝茖?dǎo)出:

因此，g表示介質(zhì)中散射角余弦的平均值。各向異性因子g取值在-1～1之間。

計(jì)算出偏轉(zhuǎn)角和方位角，新的光子軌道μ′x，）則可從舊的軌道（μx，μy，μz）和偏轉(zhuǎn)角θ與方位角ψ計(jì)算得到:

2.7 光子終止

若光子的權(quán)重衰減到足夠小，并低于一定閾值時(shí)，那么它的傳輸對(duì)計(jì)算結(jié)果的影響將可忽略。為了減小忽略權(quán)重而帶來的誤差。采用輪盤賭方法對(duì)權(quán)重小于閾值W的光子進(jìn)行處理。在一定數(shù)量m的光子中讓其中一個(gè)光子能夠存活下來，并使權(quán)重增為mW，其他的光子則為0。

3 結(jié)果

3.1 水體參數(shù)設(shè)置

根據(jù)實(shí)際水體所滿足的參數(shù)和成像探測范圍，選取深度為30m，半徑為5m的圓柱形水體區(qū)域，網(wǎng)絡(luò)劃分的精度為0.01m，進(jìn)行蒙特卡羅程序仿真實(shí)驗(yàn)，如圖2所示。設(shè)從水體的吸收系數(shù)為0.2 m-1，散射系數(shù)為 0.1 m-1，各項(xiàng)異性因子為 0.8［11～12］。為了保證一定的仿真精度，選擇入射光子包數(shù)為100萬。

圖2 水體介質(zhì)模型示意圖

3.2 單光子散射路徑

根據(jù)仿真設(shè)計(jì)的參數(shù)，首先分析單光子的散射路徑，并對(duì)單光子的散射路徑和散射方向進(jìn)行描述，如圖3所示。部分單光子經(jīng)過了水體的折射、散射和吸收后，最終又在水面出射（如圖3中上面三圖所示）；部分單光子經(jīng)過水體折射、散射和吸收后，最終被吸收掉（如圖3中左下圖所示）；部分單光子進(jìn)入水體后，直接在通過水體出射（如圖3中右下圖所示）。

3.3 水下光子密度空間分布

將100萬個(gè)光子在規(guī)定的水域中心位置，準(zhǔn)直入射到水體中，光子經(jīng)過水體的反射、折射、散射和吸收，并對(duì)100萬個(gè)光子能量進(jìn)行統(tǒng)計(jì)，可以得到不同位置的光子密度空間分布和不同位置的光子密度空間分布等高線圖，如圖4所示。Colorbar采用jet顏色映射［13］，表示以10為底的光子數(shù)對(duì)數(shù)大小。

圖3 單個(gè)光子在水體中的散射路徑

圖4 光子密度空間分布仿真圖和高線圖

3.4 不同深度下徑向光子密度分布

根據(jù)仿真數(shù)據(jù)，進(jìn)一步對(duì)比統(tǒng)計(jì)分析了在水下不同深度時(shí)，光子密度（或光子數(shù)）的徑向分布，如下圖所示。分別統(tǒng)計(jì)了水下深度為5m、15m和25m處的光子密度，其對(duì)數(shù)分布曲線如圖所示，可以看到光子密度隨著深度增加，光子密度整體上成遞減趨勢，但是在距離入射位置處光子密度比遠(yuǎn)離入射位置處的值要大。在深度為25m時(shí)，光子密度已經(jīng)非常低，達(dá)到10-4。

3.5 不同深度下散射光子分析

透過水體的光子可以分為彈道光子和散射光子，散射光子可以分為蛇形光子和擴(kuò)散光子。一般地，認(rèn)為彈道光子在傳統(tǒng)透鏡式成像中起到了主要貢獻(xiàn)，而散射光子是被水中散射子散射的光子，對(duì)成像起到干擾作用，為噪聲成分。

圖5 不同深度下光子數(shù)徑向分布

圖6 不同深度下，彈道光子去除和未去除透射光照度對(duì)比

為了對(duì)比分析不同水體深度下，彈道光子去除和未去除時(shí)對(duì)透射光子密度分布影響，分別對(duì)水深度為2m、6m、10m和15m的水體仿真研究，其對(duì)比圖如圖6所示。在每個(gè)水體深度下，左上圖是未去除彈道光子的透射光子密度分布圖像，左下圖是未去除彈道光子的透射光子密度分布圖像對(duì)應(yīng)y為0位置的光子密度分布（Colorbar以光子密度的對(duì)數(shù)形式表示出），右上圖是去除彈道光子的透射光子密度分布圖像，右下圖是去除彈道光子的透射光子密度分布圖像對(duì)應(yīng)y為0位置的光子密度分布。其他水體深度，同理可得到如圖4所示。

由圖中可以看出，未去除彈道光子時(shí)x為0位置處的光子密度數(shù)遠(yuǎn)大于去除彈道光子時(shí)，但是隨著深度增加，其比值會(huì)逐漸減小，說明散射光子在深度比較小時(shí)對(duì)成像影響比較小，但是在水體深度比較大時(shí)，彈道光子較少，散射光子也會(huì)減小，值得注意的是，彈道光子減小的比例比散射光子減小的比例要大。所以在大深度區(qū)域成像時(shí)，主要貢獻(xiàn)的直行光子會(huì)被噪聲淹沒，無法進(jìn)行成像。

4 結(jié)語

本文基于蒙特卡羅仿真思想設(shè)計(jì)了水下光場傳輸?shù)拿商乜_程序算法，選取實(shí)際水體參數(shù)，對(duì)水下單光子散射路徑、水體光子密度分布、不同深度徑向光子密度分布和不同深度下散射光子影響進(jìn)行了仿真分析。單光子路徑主要有三種情況:直接透射、最終出射水面、最終吸收。光子包垂直進(jìn)入水中時(shí)，彈道光子密度數(shù)遠(yuǎn)大于散射光子時(shí)，但是隨著深度增加，其比值會(huì)逐漸減小，說明散射光子在深度比較小時(shí)對(duì)成像影響比較小，但是在水體深度比較大時(shí)，彈道光子減小的比例比散射光子減小的比例要大。所以在大深度區(qū)域成像時(shí)，起主要貢獻(xiàn)的直行光子會(huì)被噪聲淹沒，無法進(jìn)行成像。