mie散射物理參量的一種改進數(shù)值算法

mie散射物理參量的一種改進數(shù)值算法

1考慮復折射率對計算結(jié)果的影響主要表現(xiàn)首先，可以嚴格地分解mie和debie運動，但mie運動的計算并不復雜，因為mie分散計算非常復雜。另一方面，直接編程mie理論可能會導致嚴重的計算錯誤，尤其是在計算復變量bessel函數(shù)時。到1968年Dave最先發(fā)表了完整的算法,以后Lentz和Wiscombe又針對部分計算提出了改進算法,但Dave的算法不能應用到大尺度參數(shù)的情況下,Wiscombe的方案詳細討論了計算時間,循環(huán)次數(shù),得到了準確的結(jié)果,但是他沒有考慮復折射率較大的粒子,安慶師范學院的張杰等人對大氣中具有復折射率氣溶膠粒子的散射光學特性,余其錚等人利用連分式進行計算,通過對計算次數(shù)進行擬合,提高了計算的速度,但是對于尺度參數(shù)在100以外的擬合公式?jīng)]有討論。王式民等利用后向遞推辦法大大提高了計算速度,但是文中提到的遞推表達式是和粒子的復折射率相關的表達式。朱震等人在以往Mie散射計算的基礎上提出一種新的算法,此算雖然區(qū)別于傳統(tǒng)的遞推算法,精度高,但仍然使用了遞推算法,不可避免的會出現(xiàn)積累誤差。前面所說的改進都是在同一種算法上改進,而姬豐等人利用Matlab對散射參量進行數(shù)值模擬,得到了正確的結(jié)果,但是它沒有考慮折射率虛部較大和粒子尺度參數(shù)較大對計算結(jié)果的影響。而我們的數(shù)值計算完全不依賴連分式或者后向遞推算法,直接利用Matlab中內(nèi)置的半整數(shù)階Besselj函數(shù),半整數(shù)階Besselh(由半整數(shù)階Besselj和Bessely函數(shù)組合得到)函數(shù),連帶Legendre函數(shù)通過簡單有效的編程得到了準確的結(jié)果,由于Matlab中基于矩陣的運算方式,比起用后向遞推或者改進后向遞推的編程方式有很大的優(yōu)勢。2半個數(shù)階的貝塞爾函數(shù)法若散射體為均勻球體,如圖1所示,照射光為線偏振平面波,振幅為E,光強為I0,沿z軸傳播,其電場矢量沿x軸振動。散射體位于坐標原點O,P為觀測點。散射光方向(OP方向)與照射光方向(z軸)所組成的平面稱為散射面,照射光方向至散射光方向之間的夾角θ稱為散射角,而x軸至OP在xy平面上投影線(OP′)之間的夾角φ稱為極化角。觀測點與散射體相距r。根據(jù)經(jīng)典的Mie散射理論,散射粒子的尺度參數(shù)為α=2πa/λ,其中a為球形粒子的半徑,散射粒子相對周圍介質(zhì)的折射率為m=m1+i*m2。則散射光垂直于散射面和平行于散射面的兩個分量的振幅函數(shù)為:S1=∞∑n=12n+1n(n+1)(anπn+bnτn)(1)S1=∑n=1∞2n+1n(n+1)(anπn+bnτn)(1)S2=∞∑n=12n+1n(n+1)(anτn+bnπn)(2)S2=∑n=1∞2n+1n(n+1)(anτn+bnπn)(2)以上兩式中:an=Ψn(α)Ψ′n(mα)-mΨ′n(α)Ψn(mα)ξn(α)Ψ′n(mα)-mξ′n(α)Ψn(mα)(3)bn=Ψn(α)Ψ′n(mα)-Ψ′n(α)Ψn(mα)mξn(α)Ψ′n(mα)-ξ′n(α)Ψn(mα)(4)πn=Ρ(1)n(cosθ)sinθ(5)τn=ddθΡ(1)n(cosθ)(6)an=Ψn(α)Ψ′n(mα)?mΨ′n(α)Ψn(mα)ξn(α)Ψ′n(mα)?mξ′n(α)Ψn(mα)(3)bn=Ψn(α)Ψ′n(mα)?Ψ′n(α)Ψn(mα)mξn(α)Ψ′n(mα)?ξ′n(α)Ψn(mα)(4)πn=P(1)n(cosθ)sinθ(5)τn=ddθP(1)n(cosθ)(6)以上式中:ψn(z)=(zπ2)1/2Jn+1/2(z)(7)ξn(z)=(zπ2)1/2?[Jn+1/2(z)-i*Yn+1/2(z)](8)ψ′n(z)=ψn-1(z)-nzψn(z)(9)ξ′n(z)=ξn-1(z)-nzξn(z)(10)ψn(z)=(zπ2)1/2Jn+1/2(z)(7)ξn(z)=(zπ2)1/2?[Jn+1/2(z)?i*Yn+1/2(z)](8)ψ′n(z)=ψn?1(z)?nzψn(z)(9)ξ′n(z)=ξn?1(z)?nzξn(z)(10)Jn+1/2(z)和Yn+1/2(z)分別為半整數(shù)階的第一類,第二類貝塞爾函數(shù)。P(1)n(1)n(cosθ)為一階n次第一類締合勒讓德函數(shù);Pn(cosθ)為第一類勒讓德函數(shù)。在數(shù)值模擬過程中選取初始條件如下:{ξ1(z)=1z(sinz+icosz)-(cos-isinz)ψ1(z)=1zsinz-coszτ1=cosθπ1=1(11)通過設定符號變量x和兩個符號函數(shù)ψn(x)和ξn(x),對于ψn(x)的計算直接根據(jù)式(7)運用Matlab自帶函數(shù)Besselj來實現(xiàn);對于ξn(x)的計算直接根據(jù)式(8)運用Matlab自帶函數(shù)Besselj和Bessely的組合來實現(xiàn);由于初值已給出,所以兩個導數(shù)項ψ′n(x)和ξ′n(x)分別用前面已經(jīng)求出的ψn(x)和ξn(x)根據(jù)遞推公式(9)和(10)得出;現(xiàn)在只需要運用Matlab自帶函數(shù)Inline和Vectorize將符號變量和符號函數(shù)轉(zhuǎn)化為數(shù)值變量α和mα數(shù)值函數(shù)。代入數(shù)值變量和就可以分別計算出ψn(α),ψn(mα),ψ′n(α),ψ′n(mα),ξn(α),ξn(mα),ξ′n(α),ξ′n(mα),代入(3)式和(4)式就分別可以求出an和bn,系數(shù)πn和τn可以由Matlab自帶Legendre函數(shù)得出,本文算法的具體程序流程如圖2所示:在這4個系數(shù)確定的情況下就可以計算粒子的散射系數(shù)和消光系數(shù)以及散射相位函數(shù)等物理參量。值得注意的是:Matlab自帶Legendre函數(shù)偶數(shù)行和實際值互為相反數(shù),這一點在本文作者編程的過程中深入分析才得以發(fā)現(xiàn)。3改進的反演算法在第1節(jié)中已經(jīng)將計算Mie散射參量所需的四個最基本的系數(shù)通過編程計算得到,還需對其準確性進行檢驗,根據(jù)散射系數(shù),消光系數(shù)計算公式:ksca=2α2∞∑n=1(2n+1)(|an|2+|bn|2)(12)kext=2α2∞∑n=1(2n+1)Re(an+bn)(13)我們計算了部分粒子的散射和消光系數(shù)與Dave和Wiscombe的結(jié)果作為對比如表1(程序Appendix1.),其中例1～8是Dave(D)的計算結(jié)果,例9～18是Wiscombe(W)的計算結(jié)果。(表中畫橫線的表示沒有提供相應數(shù)據(jù))。由表可以看出,我們所采用的簡單的Matlab程序得出了任意折射率且尺度參數(shù)在10-4-104的球形粒子散射參量的準確計算結(jié)果,與傳統(tǒng)的計算結(jié)果吻合得很好。充分證明了系數(shù)an和bn算法的正確性。而且我們的計算除了(9),(10)兩式用到了簡單遞推以外,其余根本不出現(xiàn)遞推,而傳統(tǒng)的后向遞推在處理粒子尺寸較大或者折射率虛部較大時,往往出現(xiàn)計算速度慢或產(chǎn)生溢出和不收斂的現(xiàn)象,本文的算法在根本上避免了以上問題。為了更充分的論證本算法的正確性和適用性,本文又從尺度參數(shù)連續(xù)改變的角度畫出了消光系數(shù)與尺度參數(shù)及折射率實部乘積的關系如圖3、4所示(程序Appendix2.):以上兩圖與參考文獻完全吻合,而上述兩圖中靠近坐標原點所對應的粒子尺度參數(shù)趨于無窮小,即便在此情況下,本文的算法得到的結(jié)果仍然與前人算法吻合,說明了本算法在粒子尺度參數(shù)很小的情況下仍然適用,對于尺度參數(shù)較大的情況,沒有可對照的參考文獻,但是從表1中可以看到在尺度參數(shù)為10000時仍然給出了準確的計算結(jié)果,至少說明在這個較大尺度范圍內(nèi)是準確的。4散射的量化見表2散射相位函數(shù)是描述光經(jīng)過粒子散射以后能量空間分布的物理量,對于微粒散射特性的研究具有重要意義,由參考文獻中關于非偏振狀態(tài)下散射相位函數(shù)的表達式:Ρ(θ)=|S1(θ)|2+|S2(θ)|2∞∑n=1(2n+1)(|an|2+|bn|2)(14)其中各個物理量的含義及其表達式在上文中已經(jīng)提到,由于粒子的散射是粒子本身的折射率以及尺度參數(shù)決定的,與其它因素無關,只須給出這兩個物理量的具體值即可。本文選取折射率為m=1.5+0.005×i,尺度參數(shù)分別為5、15如圖5、6所示(程序Appendix3.):以上兩圖與參考文獻符合得很好,由公式(1),(2)可以看出振幅函數(shù)種含有系數(shù)an、bn、πn、τn,而系數(shù)an和bn的準確性已經(jīng)通過散射系數(shù)和消光系數(shù)得以證明,由散射相位函數(shù)計算的準確性可知系數(shù)πn和τn計算的準確性。5算法的收斂次數(shù)在討論收斂速時,原則上講應該是與前人算法在相同狀態(tài)下進行比較,但由于前人算法比較繁瑣,況且程序執(zhí)行效率低,因此本文闡述的收斂速度是從收斂次數(shù)來考慮,因為在傳統(tǒng)的后向遞推和改進后向遞推中都用到了遞推次數(shù)的經(jīng)驗公式,各個文獻中所給的表達式不盡相同,但是得出的收斂次數(shù)相差不大,本文采用的后向遞推的上限次數(shù):N=1.5m1α+10(15)其中m1為粒子折射率的實部,a為粒子的尺度參數(shù)。由此可以大致確定程序的循環(huán)次數(shù),然后在此基礎上減少求和次數(shù),直到小數(shù)點后14位為穩(wěn)定輸出為止,以此作為收斂的標準,而Dave或者Wiscombe等人是將小數(shù)點后面4位或5位穩(wěn)定輸出作為最終計算結(jié)果。現(xiàn)將本文計算的級數(shù)求和次數(shù)N1與后向遞推的上限次數(shù)N進行對比如表2.由表2可知:本文的算法即使是在保證小數(shù)點后面14位輸出穩(wěn)定作為循環(huán)終止的前提下,循環(huán)次數(shù)也比后向遞推所用的上限次數(shù)少,假如是在同樣5位輸出穩(wěn)定的情況下,我們的算法在求和次數(shù)上將會更少,計算速度將會更快,本算法在計算粒子尺度參數(shù)或者折射率實部較大的粒子方面更具有優(yōu)勢。6davd/wiso模擬算法本文拋棄了傳統(tǒng)的關于經(jīng)典Mie散射參量的后向遞推和連分式法,充分利用Matlab內(nèi)部集成的函數(shù)模塊,對Mie散射消光系數(shù),散射系數(shù)進行了準確的計算,計算結(jié)果和Dave或者Wiscombe的結(jié)果吻合得很好,本算法的一大優(yōu)點在對于尺