Mie 散射系數的新算法王少清 任中京　張希明　何芳 江海鷹(山東建筑材料工業(yè)學(xué)院，廈門(mén)　250022)提要：介紹了一種計算Mie散射系數的新方法，給出了計算實(shí)例。關(guān)鍵詞：米氏散射，Mie系數，Mie計算New algorithm of Mie scattering parametersWang Shaoqing Ren Zhongjing　Zhang Ximing 　HeFang Jiang Haiying(Shandong Institute o f Building Materials，Jinan　250022)Abstract：A new algorithm of calculating Mie Parameters was introduced. Some calculation results done with this algorithm were given.Key words: Mie scattering ，Mie parameters，Mie calculation1　引言Mie 理論是球形顆粒對單色光的散射場(chǎng)分布的嚴格解析解^[1]，目前在環(huán)保、動(dòng)力、氣象、天文、兩相流及粉體顆粒尺寸分布測量等工程技術(shù)領(lǐng)域中有廣泛的應用。利用單一顆?；蝾w粒群光散射場(chǎng)的測量數據?？梢苑赐频蒙⑸漕w?；蝾w粒群的很多物理性質(zhì)，如顆粒的尺寸、顆粒的折射率等^[2]。但反推必須事先計算出各種尺寸的顆粒在各種復折射率下的散射場(chǎng)分布數據。1968 年Dave^[3]zui先發(fā)表了完整的Mie散射計算方法，以后Lentz^[4]和Wiscomb^[5]又針對部分計算提出了新算法。國內也有人發(fā)表了他們自己的算法^{[6、7、8]}。但總的看來(lái)，這些算法均有各自的局限性。尤其是當顆粒尺寸或折射率的虛部值較大時(shí)，往往計算速度過(guò)慢或產(chǎn)生溢出和不收斂等現象。本文介紹作者發(fā)展的Mie 散射新算法。該算法的特點(diǎn)是不受顆粒尺寸及折射率的限制，不會(huì )產(chǎn)生溢出和不收斂的現象，且具有較快的計算速度。2　Mie散射系數的計算公式Mie 散射計算的中心問(wèn)題是計算Mie 散射系數a_n 和b_n，其表達式為^[9]
其中為顆粒的尺寸參數，定義為α=πd/λ，d為顆粒直徑，λ為入射光在顆粒周?chē)橘|(zhì)中的波長(cháng)，而m為顆粒在周?chē)橘|(zhì)中的相對復折射率，即m= m¹ + im² ( m²< 0) , 式中i 為虛數單位。而Ψ_n( Z) 和ξ_n( Z) ( Z 既表示α，又表示mα)的表達式為Ψ_n( Z)=（πZ/2）J_n+1(Z)ξ_n( Z)=Ψ_n( Z)+iΧ_n(Z)Χ_n(Z)= ? (πZ/2)( ? 1)_n-1J_-(n-1)(Z) = ? (πZ/2)N_n+1(Z)3 計算中產(chǎn)生溢出的原因計算Mie 散射系數須先計算Ψ_n 和Χ_n。一般采用遞推的方法。遞推又分為向前遞推( 即從n= 0 開(kāi)始) 與向后遞推( 即從n= N 開(kāi)始至n=0, N 為預先設定值) 。實(shí)驗表明，向前遞推總是快于向后遞推。Ψ_n 和Χ_n的初值為
分析以上兩式可知，當m₂≠0 時(shí), 若顆粒尺寸d很大, 或復折射率的虛部值m₂很大，將使得乘積m₂d很大，可使兩式中的項exp( - m₂ α) =exp ( - πm₂d/ λ) 的值超過(guò)計算機的數據限，從而產(chǎn)生溢出，這是產(chǎn)生溢出的重要原因。另外，在遞推過(guò)程中，不恰當的算法也可能造成溢出。4　a_n和b_n的新算法為解決上述問(wèn)題，作者提出了新的算法。將a_n和b_n的公式變形如下: 令
其中Lnr、Lnj分別表示Ln(m)的實(shí)部與虛部。將（2）式代入（1）式，并用a_nr、a_nj和b_nr、b_nj分別表示的實(shí)部與虛部。如此可推得
在上述四個(gè)公式中采用比值的形式是非常重要的，這樣可避免遞推過(guò)程中當a_i與b_i 較大時(shí)乘法運算可能產(chǎn)生的溢出。這是本算法的一個(gè)重要特點(diǎn)。在以上四式中
由于均為實(shí)變量函數，計算不產(chǎn)生溢出。關(guān)鍵是和的算法如何處理，才能保證計算中不產(chǎn)生溢出。在Lentz的算法中是采用連分式計算Ln的值，其精度的保證是由在大量計算基礎上得出一個(gè)截斷項數N與參數a與m的經(jīng)驗公式而實(shí)現的。這樣的經(jīng)驗公式，一則有實(shí)用上的局限性，再則也會(huì )帶來(lái)截斷誤差。文獻（6）對此經(jīng)驗公式做了改進(jìn)，但仍陷于a=1~100，m1=1~2，m2=0~1的范圍。下面介紹本文作者發(fā)展的關(guān)于Ln的新算法。該算法的特點(diǎn)是不受a及m值的限制，不會(huì )產(chǎn)生溢出或不收斂等病態(tài)現象，且具有較快的計算速度。令


上面導出的(3) - (20)諸式，構成了Mie系數an和bn的完整算法。由于an和bn是從n= 1開(kāi)始計算，利用初值公式(16) - (20)即可算得任意級數的an和bn的值，故沒(méi)有舍入誤差的問(wèn)題。從(16)式可見(jiàn)，因為y= m₂ɑ≤0, 故無(wú)論m₂和ɑ取何值均不會(huì )產(chǎn)生溢出，再加上(3)的各式中采用了比的形式，又避免了計算過(guò)程中的溢出，這就從根本上解決了溢出的問(wèn)題。5　計算實(shí)例利用以上算法編制了計算顆粒散射場(chǎng)強度和消光系數的計算機程序。當波長(cháng)為λ的單位振幅的平面自然光入射顆粒時(shí)，顆粒的散射光強為^[9]

Ξ為計算機在雙精度下的zui小數據限。圖1示出散射光強的一組計算實(shí)例。其中取m₁=1.33，m₂=-0.4，λ=0.6328三圖分別對應于顆粒的直徑為d=0.001,1.0和30μm。而d圖為顆粒直徑d=100μm時(shí)散射花樣的局部放大?？梢?jiàn)，隨著(zhù)顆粒尺寸的增大，前向散射迅速加強，并且出現了復雜的旁瓣。
實(shí)部(a)與虛部(b)的變化情況?？梢?jiàn)隨著(zhù)m₁和m₂的增大，雖然顆粒的尺寸保持不變，但散射也加強，且后向散射隨著(zhù)m₁和m₂的增大而加強。
圖3給出有關(guān)消光系數的計算結果，其中a)和b)分別表示消光系數隨折射率實(shí)部與虛部的變化情況?？梢?jiàn)隨著(zhù)顆粒直徑的增大，消光系數趨近于2; 折射率的增大，尤其是折射率虛部的增大，使這一趨近變得更快和更明顯。另外當折射率的虛部m₂= 0時(shí)，消光系數隨顆粒直徑的增大而振蕩；但當m2≠0 時(shí)，振蕩迅速消失。

光強zui大值處所對應的F_M(Z)的值即由上式確定。由上式也可見(jiàn)：此極限情況下的焦移大小主要由S₀/f及N_a所決定。
參考文獻：[1]M. Bor n and E. Wolf, Principles of optics ( sixt hedition ) , Pergaman Press ( OXFORD, NEW YORK, PARIS) , pp611(1980)[2]Peter Chylek, V. Ramasw amy , A. Ashkin and J.M . Dziedzic,“Simultaneous determination o f refractiv endex and size o f spherical dielectric particles from light scattering data",Appl. Opt . ,Vo l. 22, No . 15, p2302- 2307(1983)[3]J. V. Dave, Repo rt 320 - 3237 ( IBM Scientific center , 1968)[4]W. J. Lentz, Appl. Opt. , 15, 668( 1976)W. J. Wiscome, Appl. Opt, Vol. 19, No. 9, pp1505( 1980)[6]顧冠亮等, 有關(guān)光散射物理量的數值計算, 上海機械學(xué)院學(xué)報, 1984 年4 期pp21[7]余其錚等, Mie 散射算法的改進(jìn), 哈爾濱工業(yè)大學(xué)學(xué)報, 1987 年3 期, pp21[8]鄭剛等, M ie 散射的數值計算, 應用激光, Vol.12, No . 5, pp220( 1992)[9]H. C. v an de Hulst , Lig ht scatter ing by small particles,Do ver Publica tion lnc. , New Yor k, chp. 13( 1981)