呂且妮 ，徐 暢 ，靳文華

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

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

1 GOM理論模型

1.1 散射光強函數(shù)

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

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

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

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

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

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

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

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

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

波前擴展因子為

式中

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

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

E‖和 E⊥分別為

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

散射光場的強度分布

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

1.2 散射相位函數(shù)

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

則其相位函數(shù)P()θ為

相位函數(shù)P()θ為

2 入射角-散射角的數(shù)值求解

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

2.1 函數(shù)的構(gòu)造

2.1.1 tan(θp/2)函數(shù)

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

設(shè)

2.1.2 sin(θp/2)函數(shù)

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

設(shè)

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

2.1.3 sinpθ函數(shù)

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

設(shè)

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

2.2 多項式方程的求解

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

2.3 結(jié)果及分析

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

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

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

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

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

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

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

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

圖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理論計算的相位函數(shù)分布Fig.8 Scattering phase function computed by GOM theory

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

4 結(jié) 語

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

［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］ 李旭東，楊鴻儒，張 郁，等. 大氣泡散射的幾何物理模型數(shù)值計算[J]. 應(yīng)用光學(xué)，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］ 呂且妮，靳文華，葛寶臻，等. 粒子散射的幾何光學(xué)模型中的多值散射角函數(shù)計算[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. 光學(xué)原理[M]. 楊霞蓀，譯. 北京：電子工業(yè)出版社，2009.Born M，Wolf E.Principles of Optics[M]. Yang Xiasun，Trans. Beijing：Publishing House of Electronics Industry，2009(in Chinese).