基于幾何光學近似模型的大氣泡粒子散射光場分布計算

呂且妮 ，徐 暢 ，靳文華

(1. 天津大學精密儀器與光電子工程學院，天津300072；2. 天津大學光電信息技術教育部重點實驗室，天津300072)

在基于粒子光散射理論的各種檢測方法中，粒子散射光特性分析是其主要研究內容之一．散射光強計算常用的理論模型有 Mie理論[1-2]、Debye級數展開模型[3-4]、幾何光學近似模型(geometrical optics approximation model，GOM)[2,4-13]．GOM 是基于光線理論計算每一束光線對散射場的貢獻，并可以給出每一束光線的物理意義．在 GOM 模型中，由于散射角通常是入射角的多值函數，因此，利用 GOM 計算粒子散射光強分布，就必須求解出每個散射角所對應的入射角，也就是說，基于 GOM 模型的數值計算首先是要求解散射角所對應的所有入射角，這是基于GOM 理論計算散射光場分布必須首先要解決的問題．其方法之一是利用三角函數關系構建一個方程，求解該方程的解得到滿足條件的入射角，如 Zhou等[11]利用cos(·)函數構造了散射角與入射角的多項式，采用牛頓迭代法和分段法相結合的方法求解多項式方程得到入射角，并開發(fā)了相對折射率大于1的大顆粒的 GOM 計算程序．文獻[4]和文獻[12]利用sin(·/2)函數構造了散射角與入射角的多項式，基于MATLAB語言編寫了相對折射率小于1的大氣泡的GOM 計算程序．呂且妮等[13]已經提出利用sin(·/2)和 MATLAB中的 fzero(·)函數的遍歷取值的求解方法．本文分析了基于不同三角函數方程，利用MATLAB 中的 fsolve(·)函數/fzero(·)函數，采用遍歷取值求解入射角的準確性和難度，給出了一種簡便和準確的求解方法，計算了球形粒子的散射光強分布和散射相位分布，與 Mie理論計算進行了比對，并與文獻[8]和文獻[10]給出的計算結果進行了比較，結果吻合很好．

1 GOM理論模型

1.1 散射光強函數

圖 1中，設波長為λ的平面平行光Ei入射于半徑為a的球形粒子，散射光場分布包括反射光E0、折射光E1、E2、…、Ep( p=0,1,2,… ,為弦數)和衍射光?，其復振幅分布為

衍射光場Ediff的平行和垂直分量的復振幅分布分別為

式中：θ為散射角；J1為第一類Bessel函數；x為粒子尺寸的無因次參數，

圖1 GOM理論模型Fig.1 GOM theoretical model

反射光場E0和折射光場Ep的平行和垂直分量的復振幅[2]分別為

式中：p=0時為反射光，即E0；p=1時表示入射光折射進入粒子后再直接折射(透射)出粒子，即E1；p＞ 1表示入射光折射入粒子后，在粒子中發(fā)生(p?1)次反射，再折射(透射)出粒子．經粒子后的出射光Ep由入射角θi和弦數p唯一決定，其偏向角θp為

式中：s inθi=msinθr，θr為折射角，m=n2n1為粒子相對于介質的折射率；偏向角θp與散射角θ關系可表示為

式中q=±1，分別表示入射光從光軸的上半方和下半方入射．如圖2所示，

圖2 q=±1時入射光的入射位置Fig.2 Incident ray position for q=±1

波前擴展因子為

式中

式中：Ei為入射光振幅；r⊥和r‖為費涅耳反射系數[14]，分別為

對m＜1，當入射角θi＞θC=arcsinm (臨界角)時產生全反射，這時沒有光入射到粒子內部．全反射引入的相位變化[14]為

E‖和 E⊥分別為

散射光場總的相位分布[8]為

散射光場的強度分布

則非偏振光的散射光強分布為

1.2 散射相位函數

散射光能量分布的另一種描述是散射光的相位函數，它決定了光散射能量的角分布．非偏振光的散射光相位函數[9-11]為

則其相位函數P()θ為

相位函數P()θ為

2 入射角-散射角的數值求解

由式(15)～式(21)可知，散射光強I()θ和散射相位P()θ是散射角θ的函數，由式(4)和式(5)可知，散射角θ是入射角θi的函數，且是入射角θi的多值函數．對一散射角θ，其散射光強分布可能是由不同入射角的光線所共同作用的結果，因此，基于 GOM 模型的數值計算關鍵是已知弦數p和散射角θ，求滿足式(4)和式(5)的所有入射角θi．分析式(4)可知，根據三角函數的特性，采用不同的函數將有不同形式的求解方法．

2.1 函數的構造

2.1.1 tan(θp/2)函數

根據正切函數特性，對式(4)有

設

2.1.2 sin(θp/2)函數

根據正弦函數的特性，由式(4)得

設

求解方程f(θi)=0的根，即可得到入射角θi．

2.1.3 sinpθ函數

對式(4)，不再對偏向角pθ取半角，直接采用周期為2π的sinpθ函數求解，即有

設

同樣，求解方程 f(θi)=0的根，即可得到相應的入射角θi．

2.2 多項式方程的求解

構造多項式之后，將式(22)、(24)和(26)中的角度求解轉化為求方程 f(θi)根的問題．本文采用MATLAB 中 fsolve(·)函數/fzero(·)函數求解方程其步驟為：對一給定的散射角θ和 p，設作為 fsolve(·)函數的求解估計值，其中Δθi為所設定的入射角取值間隔，n為正整數．對q=±1時，遍歷所有估計值，返回求出的入射角θi，再將θ、p、q=±1以及求出的入射角θi代入式(5)，判斷l(xiāng)是否為整數，若l為整數，保存入射角θi，否則說明給定的散射角θ無解．當遍歷的較小時，通過 fsolve(·)函數可準確地求出每一個入射角θi．這種求解方法與文獻[11]相比，不需要預判拐點．

2.3 結果及分析

圖3 m＝0.75時散射角與入射角關系Fig.3 Relationship between scattering angle and incidence angle for m = 0.75

分析式(24)可知，當l為偶數時，sin( 2)pθ=當l為奇數時，，因此，式(24)可表示為．由于則式(24)為．數學上，l的奇偶性并不影響式(24)，利用該函數同樣可得到正確的入射角θi．對m＝0.75，其散射角所對應的入射角如圖3所示，但所求出的入射角θi對應的q值將可能不正確．如對散射角θ=53°，求解得到的時，l=?1．此時，應取q=1，才能滿足l為整數，進而得到正確的散射光強分布和散射相位分布．因此，在利用求解時，需要根據l的值，對q值進行修正．

對式(24)，l的奇偶性影響了 q的求解，但不影響入射角θi．為了避免l對q的求解影響，利用sinθp函數進行求解．圖 4(a)為對m=0.75，利用該求解算法得到的 p＝4階的入射角與散射角關系圖．與圖3(e)相比，在圖4(a)中，一個入射角θi對應兩個散射角θ′和這顯然不對．對散射角θ′和π?θ′，由式(5)，則有

圖4 m＝0.75、p＝4時散射角與入射角關系Fig.4 Relationship between scattering angle and incidence angle for m＝0.75,and p＝4

若設正確的散射角為θ′，則而將再不滿足整數的限制條件，故散射角 π?θ′將不正確，應予以剔除．修正后的入射角與散射角關系如圖 4(b)所示，與圖 3(e)相同．

圖5 m＝1.333時不同階數p的散射角與入射角關系Fig.5 Relationship between scattering angle and incidence angle of various components p for m＝1.333

對上述 3種函數，采用 fzero(·)函數也可得到正確的θi和q值．文獻[13]采用sin(θp2)函數和fzero(·)函數求解得到了相應的入射角．但對p≥4，散射角θ≈180°時，出現了入射角θi的缺值現象，主要原因是：fzero(·)函數是采用二分法或分割法等算法進行零點的求解，文獻[13]在利用 fzero(·)函數計算時，將中間值作為估計值代入 fzero(·)函數進行求解，將區(qū)間值作為輸入值即可．MATLAB中的fzero(·)函數和fsolve(·)函數都可求解給定初值附近的零值點，但fzero(·)函數的運算速度稍快，而 fsolve(·)函數可減小由估計值選擇不當引起的缺值現象．

圖6 m＝0.75時GOM和Mie理論計算的總散射光強比較Fig.6 Comparison of total scattering intensity calculated by GOM and Mie theories for m＝0.75

3 散射光場分布計算

利用上述的入射角-散射角求解方法以及式(15)～式(21)，分別計算了不同尺寸和不同折射率的球形粒子散射光強和散射相位函數分布．圖6給出了m＝0.75、x取不同值時，基于 GOM 和 Mie理論計算的不同粒徑的總散射光強分布，其中，截止階數，基于GOM計算光強時，已計算了由全反射引入的相移量，圖 6與文獻[8]中的圖 3完全吻合．圖7為m＝1.333時，基于GOM和Mie理論計算的 x＝100、500的總散射光強分布，同樣pmax=10.圖6和圖7中，隨著粒子直徑d的增大即x的增大，GOM 與 Mie計算的散射光強分布吻合越好．對圖6，在θ=80°區(qū)域，GOM與Mie計算的散射光強分布有差異，是因為當θi=θC時，p=1階的散射光劇烈減少，I→0．根據幾何光學理論，當θ=π?2θC=82.82°時，沒有光線進入到粒子內部產生高階數的散射光，也就說散射角θ≤82.82°，即散射臨界角θC=π?2θC=82.82°(可參看文獻[13]中的圖 6(b)).對圖 7，在θ=130°區(qū)域，GOM 理論計算的光強分布出現了尖峰，該角度位置為光學彩虹角位置，此時擴散因子 D(p)(p, m,θi)→∞ ．圖 8(a)所示為m＝0.75、m＝1.333時，基于 GOM 理論計算的球形粒子近場散射相位分布．由式(19)可知，近場散射相位分布與粒子尺寸和波長無關，圖 8(a)與文獻[10]中的圖 2.7完全吻合．圖 8(b)給出了m＝0.75和m＝1.333，基于 GOM 理論計算的不同粒徑的球形粒子的遠場相位圖，其中截止階數pmax=10．由圖 8(b)可以看出，隨著粒徑的增大，散射相位函數分布趨于一致．這是因為當粒子尺寸不斷增大時，夫瑯和費衍射光場Ediff的貢獻將逐漸減小，J1( x sinθ)→0．圖 9(a)和圖9(b)分別為波長λ＝0.55,μm，粒子半徑a＝15,μm時，m＝0.75和m＝1.333的球形粒子遠場相位分布，并與 Mie理論計算的非均勻粒子體系的散射相位分布進行了對比，其中對于 Mie理論，非均勻粒子分布服從半寬系數μ＝6的伽馬分布、粒子有效半徑a＝15,μm，圖9與文獻[10]中的圖2.6相同．計算結果驗證了該算法的正確性．

圖7 m＝1.333時GOM和Mie理論計算的總散射光強比較Fig.7 Comparison of total scattering intensity computed by GOM and Mie theories for m＝1.333

圖8 GOM理論計算的相位函數分布Fig.8 Scattering phase function computed by GOM theory

圖9 GOM和Mie理論計算的散射相位函數分布Fig.9 Calculations of scattering phase function computed by GOM and Mie theories

4 結 語

本文對 GOM 基本理論進行了研究，并對 GOM中的入射角與散射角關系進行了分析，對比了不同三角函數的求解方法．sin(θp2)函數和sinθp函數兩種方法均由于l的影響，使計算得到的 q值或入射角與散射角關系有所偏差．tan(θp/2)函數，由于其π的周期，無論l取何值時，式(22)將始終成立．相比前兩種方法，無需加入任何約束條件以驗證入射角θi和q值的求解是否正確，改善了程序的復雜度和運算量．利用本文所提出的多值散射角求解方法，計算了球形大氣泡粒子和水粒子的總散射光強和散射相位分布，并與 Mie理論計算進行了比較，計算結果與 Mie理論計算結果吻合很好，驗證了本算法的可行性和正確性.

［1］ Gogoi A，Choudhury A，Ahmed G A. Mie scattering computation of spherical particles with very large size parameters using an improved program with variable speed and accuracy[J].JMod Opt，2010，57(10)：2192-2202.

［2］ Van de Hulst H C.Light Scattering by Small Parti-cles[M]. New York：Wiley，1957.

［3］ Shen Jianqi，Wang Huarui. Calculation of Debye series expansion of light scattering [J].Appl Opt，2010，49(13)：2422-2428.

［4］ Wu L，Yang H R，Li X D，et al. Scattering by large bubbles：Comparisons between geometrical-optics theory and Debye series [J].J Quant Spectrosc Radiat Transf，2007，108(1)：54-64.

［5］ Xu F，Cai X S，Ren K F. Geometrical-optics approximation of forward scattering by coated particles [J].Appl Opt，2004，43(9)：1870-1879.

［6］ Li X Z，Han X E，Li R X，et al. Geometrical-optics approximation of forward scattering by gradient-index spheres [J].Appl Opt，2007，46(22)：5241-5247.

［7］ Li W，Yang K C，Xia M，et al. Computation of the scattering intensity distribution for natural light scattered by an air bubble in water [J].J Opt A：Pure Appl Opt，2006，8(1)：93-99.

［8］ Yu H T，Shen J Q，Wei Y H. Geometrical optics approximation of light scattering by large air bubble [J].Particuology，2008，6(5)：340-346.

［9］ Kokhanovsky A A. Optical properties of bubbles [J].J Opt A：Pure Appl Opt，2003，5(1)：47-52.

［10］ Kokhanovsky A A.Light Scattering Media Optics：Problems and Solutions[M]. Praxis，Chichester，UK：Springer，2004.

［11］ Zhou X B，Li S S，Stamnes K. Geometrical-optics code for computing the optical properties of large dielectric spheres [J].Appl Opt，2003，42(21)：4295-4306.

［12］ 李旭東，楊鴻儒，張 郁，等. 大氣泡散射的幾何物理模型數值計算[J]. 應用光學，2006，27(6)：539-542.Li Xudong，Yang Hongru，Zhang Yu，et al. Numerical calculation of light scattering caused by large spherical bubbles with geometrical-physical model[J].J Appl Opt，2006，27(6)：539-542(in Chinese).

［13］ 呂且妮，靳文華，葛寶臻，等. 粒子散射的幾何光學模型中的多值散射角函數計算[J]. 光電子·激光，2010，21(11)：1677-1682.Lü Qieni，Jin Wenhua，Ge Baozhen，et al. Calculation of the multi-value scattering angle function in geometrical optics model for particle scattering[J].J Opt·Laser，2010，21(11)：1677-1682(in Chinese).

［14］ Born M，Wolf E. 光學原理[M]. 楊霞蓀，譯. 北京：電子工業(yè)出版社，2009.Born M，Wolf E.Principles of Optics[M]. Yang Xiasun，Trans. Beijing：Publishing House of Electronics Industry，2009(in Chinese).