胡帥 高太長李浩 楊波 江志東 陳鳴 李書磊

1)(解放軍理工大學氣象海洋學院,南京 211101)

2)(解放軍理工大學,電磁環(huán)境效應與電光工程國家級重點實驗室,南京 210007)

3)(海軍航空工程學院,青島 266041)

基于時域多分辨算法的非球形氣溶膠散射特性仿真模擬?

胡帥1)高太長1)2)?李浩1)楊波2)江志東3)陳鳴1)李書磊1)

1)(解放軍理工大學氣象海洋學院,南京 211101)

2)(解放軍理工大學,電磁環(huán)境效應與電光工程國家級重點實驗室,南京 210007)

3)(海軍航空工程學院,青島 266041)

(2016年9月19日收到;2016年11月30日收到修改稿)

非球形氣溶膠的散射特性是影響輻射傳輸模擬準確性的重要因素.為實現(xiàn)非球形、非均質(zhì)氣溶膠散射特性的模擬,基于MRTD(multi-resolution time-domain)方法建立了一個新的氣溶膠散射模型.采用MRTD技術(shù)實現(xiàn)了近場電磁場的計算;考慮氣溶膠的特殊性,推導了基于體積積分方法的近遠場外推方法,實現(xiàn)了粒子散射振幅矩陣和穆勒矩陣的仿真;構(gòu)建了粒子吸收和消光截面的計算模型,實現(xiàn)了粒子積分散射特性的高精度模擬.將MRTD散射模型的結(jié)果與Mie理論、T矩陣法進行了對比,驗證了模型的準確性;討論了空間網(wǎng)格粗細對模擬精度的影響,并定量分析了模型的運行效率.結(jié)果表明,MRTD散射模型的相函數(shù)模擬誤差小于8%,其中前向散射方向小于4%;當粒徑與入射光波長相當時,消光和散射效率因子的相對誤差小于0.1%;空間網(wǎng)格粗細對模擬精度影響顯著,當粒子尺度參數(shù)小于20時,在相同模擬精度要求下,所需網(wǎng)格尺寸隨尺度參數(shù)呈先增大后減小的特征.

非球形氣溶膠,散射特性,時域多分辨算法,穆勒矩陣

1 引 言

非球形氣溶膠散射特性的不確定性一直是制約輻射傳輸和氣候模擬精度的重要因素[1,2].對比多年的IPCC報告可發(fā)現(xiàn)[3],氣溶膠直接輻射強迫的評估精度雖然有所提高,但仍存在較大的不確定性(-0.9—-0.1 W/m2),其中氣溶膠和卷云中冰晶等非球形粒子的散射特性就是制約其評估精確性的一個重要原因[4,5].氣溶膠光學遙感及大氣校正等均需要以輻射傳輸模式為基礎(chǔ)[6,7],而無論是標量還是矢量輻射傳輸模式,均需要氣溶膠散射特性作為輸入;目前考慮多次散射的輻射傳輸求解模型如DISORT(discrete-ordinate-method),RT3/PolRadtran及SHDOM(spherical harmonics discrete ordinate method)等已得到廣泛檢驗[8,9],但作為模式輸入?yún)?shù)的氣溶膠散射特性,如穆勒矩陣和散射相函數(shù)等[10],仍主要通過Mie散射計算得到,并未充分考慮氣溶膠及冰晶的非球形效應,而實際非球形粒子與等效球形粒子的散射特性存在較大差異[11,12],這將影響遙感正演模型的準確性,進而影響反演精度,因此準確獲取非球形氣溶膠的散射特性,提高正演模型模擬精度已成為發(fā)展大氣遙感技術(shù)的關(guān)鍵環(huán)節(jié)之一[13,14].

為模擬非球形氣溶膠的散射特性,目前已經(jīng)建立了許多氣溶膠散射模型,大致可分為兩大類：一類是近似計算模型,主要針對粒子尺度參數(shù)較為極端的情形(當入射光波長遠大于粒徑或遠小于粒徑的情形),主要的模型包括Rayleigh近似、Rayleigh-Gans-Stevenson近似(RGA)、異常衍射近似和幾何光學近似GOA(geometric optics approximation)等[15],其中GOA方法廣泛用于大尺度冰晶粒子的散射特性計算[14,16];第二類是精確求解模型,這類模型通過在一定邊界條件下數(shù)值求解Maxwell方程組或亥姆霍茲方程,獲取模擬粒子電磁散射特性[15].按照計算原理,精確求解模型可大致分為三類：基于基函數(shù)展開的散射模型、基于體積積分方程的散射模型及基于微元法的散射模型.基于基函數(shù)展開的散射模型的特點是基于矢量波函數(shù)對電磁入射場、內(nèi)場及散射場進行展開,然后求解微分方程或積分方程建立入射場與散射場展開系數(shù)間的關(guān)系,進而實現(xiàn)對散射過程的求解.這一類模型主要包括T矩陣法(T matrix method)[17],EBCM(extended boundary condition method)[18],SVM(separation of variables method)[19]和PMM(point matching method)[20]等.這類方法具有計算精度高、速度快、內(nèi)存消耗小的優(yōu)點,但當粒子尺度參數(shù)較大、形狀較為極端或復折射率較大時,易出現(xiàn)數(shù)值求解的不穩(wěn)定;散射場及內(nèi)部場采用基函數(shù)展開時,展開系數(shù)存在截斷誤差;此外,受求解邊界條件的限制,這類模型只能模擬一些形狀簡單的粒子,如橢球粒子、圓柱粒子等.體積積分方程法是通過數(shù)值求解亥姆霍茲方程的格林函數(shù)形式解來模擬散射場,典型代表是矩量法MoM(method of moments)[21],DDA(discrete dipole approximation)[22,23]和FIEM(fredholm integral equation method)積分方程法.DDA是MoM的改進,其本質(zhì)差異在于自身感應電場的計算,可適用于任意形狀、非均質(zhì)的氣溶膠粒子,但是模擬過程涉及大型矩陣方程的求解,存在不穩(wěn)定性;同時,計算精度取決于偶極子的數(shù)目,計算精度與內(nèi)存消耗、收斂速度存在突出矛盾.因此,該方法主要適合尺度參數(shù)較小的粒子的散射特性模擬,且計算精度不高.FIEM積分方程法的優(yōu)點是積分過程僅和粒子尺度參數(shù)和形狀有關(guān),而與入射光的入射角和偏振狀態(tài)無關(guān),但只適用于尺度和復折射率較小的規(guī)則粒子,且模型計算量很大[15].基于微元法思想的散射模型主要包括時域有限差分法(finite difference time domain,FDTD)[15,16]和有限元法(finite element method,FEM)[24].其中FDTD已廣泛應用于氣溶膠散射特性模擬[16],由于該方法避免了積分方程法的奇異核問題,因此更適用于復雜形狀、非均質(zhì)氣溶膠的光散射特性模擬.但與DDA相似,該方法存在計算精度與計算量的矛盾.為解決上述問題,Liu等[25]將時域偽譜法(pseudo spectral time domain,PSTD)引入氣溶膠散射過程模擬中,由于該方法取相對較少的采樣點便可獲得較高模擬精度,因此可顯著降低內(nèi)存消耗與運行時間.目前,Liu等[26]已采用該方法實現(xiàn)了尺度參數(shù)達200的粒子散射特性模擬,但受原理限制,該方法只能采用純散射場技術(shù)引入入射波,因此在計算不同形狀粒子的散射特性時,相應程序也要進行修改,過程繁瑣.MRTD是矩量法和FDTD的結(jié)合體[27,28],是Maxwell方程的高階差分展開,具有良好數(shù)值色散特性,也廣泛用于大目標的散射特性的模擬[29].基于此,本文將其引入非球形氣溶膠的散射特性的模擬中,結(jié)合研究對象和學科需求的特殊性,建立了一套新的非球形氣溶膠散射模型.

本文思路如下：首先介紹了一些物理概念和MRTD散射模型的基本框架,然后簡單介紹了基于MRTD近場計算方案,包括基本離散迭代格式、吸收邊界設(shè)計及總場散射場技術(shù);針對氣溶膠的特殊性,推導了基于體積積分的近遠場外推方案和粒子散射特性計算方法;最后對模型的準確性和性能進行了分析和驗證.

2 基本概念

2.1 Maxwell旋度方程組及電磁參數(shù)轉(zhuǎn)換

電磁場傳播規(guī)律可用Maxwell旋度方程組描述.無源條件下,該方程的矢量表達形式如下：

上式中,E和H分別表示電場和磁場的強度矢量;ε表示介質(zhì)的介電常數(shù);μ表示磁導系數(shù);σ和σm分別表示介質(zhì)的電導率和磁導率.

在大氣科學領(lǐng)域,粒子的光學屬性通常采用復折射率m表示,其實部和虛部分別表征了粒子散射效應和吸收效應的強弱.粒子復折射率m與復介電常數(shù)εc存在以下轉(zhuǎn)換關(guān)系：若假設(shè)復折射率m=mr-jmi,則有關(guān)系式成立. 采用以上關(guān)系便可將粒子的復折射率轉(zhuǎn)化為電磁散射模擬所需的電磁參數(shù).

2.2 穆勒矩陣和Stokes矢量

穆勒矩陣F(?rsca,?rinc)是輻射傳輸模式的重要參數(shù),它不僅表征了散射能量的空間分布,同時也定量描述了入射光Stokes矢量Sinc和散射光Stokes矢量Ssc轉(zhuǎn)換關(guān)系.穆勒矩陣F的定義如(4)式所示[30],其中F11也稱為散射相函數(shù).

3 MTRD非球形氣溶膠散射模型

3.1 模型的基本框架

MRTD非球形氣溶膠散射模型包括三個模塊：電磁場近場計算模塊、近遠場外推模塊及氣溶膠散射特性計算模塊.電磁場近場計算模塊主要實現(xiàn)包含散射體空間的電磁場模擬,其中關(guān)鍵技術(shù)包括Maxwell方程離散迭代格式、吸收邊界的設(shè)計以及基于總場散射場技術(shù)的入射波引入;近遠場外推模塊的功能是將近場電磁場變換至頻域,并實現(xiàn)近遠場場變換;氣溶膠散射特性計算模塊主要實現(xiàn)粒子吸收、消光和散射截面及穆勒矩陣等參數(shù)的模擬.由于氣溶膠一般為電介質(zhì),本文基于介質(zhì)中的Maxwell方程組推導了基于體積積分遠場變換方法;在粒子散射特性計算方面,本文從大氣輻射學需求出發(fā),推導了粒子的吸收截面和消光截面計算模型.

圖1 MTRD非球形氣溶膠散射模型基本框架Fig.1.The basic frame of MRTD scattering model for nonspherical aerosol particles.

MRTD散射模型的計算域是在直角坐標系下進行離散的,電磁場空間離散格式采用Yee元胞.由于散射體形狀不規(guī)則,在不同介質(zhì)的交界處,單個Yee元胞內(nèi)可能存在多種介質(zhì)混合的現(xiàn)象,若直接對其進行截斷,則可造成較大的“階梯誤差”[4],因此需要對元胞內(nèi)電磁參數(shù)進行等效平均.本模型采用Maxwell-Garnett公式進行計算[4],具體計算方法如(5)式所示.

式中,εr(x,y,z)表示位置(x,y,z)處的相對介電常數(shù),表示復介電常數(shù)εc與ε0之比,r是元胞的平均相對介電常數(shù)值.

3.2 基本迭代格式的建立和穩(wěn)定性條件

在直角坐標系下,分別采用Daubechies尺度函數(shù)φi(x)與矩形脈沖函數(shù)hn(t)在時間和空間上對電磁場分量進行展開[31]：式中r=(x,y,z)為場分量的位置向量;為電場分量的展開系數(shù);和表示磁場分量的展開系數(shù);i,j,k分別表示場分量的空間坐標;n表示時間迭代步數(shù).hn(t)滿足正交性,其定義如下：

Daubechies尺度函數(shù)φi(x)定義如下：

式中M1為Daubechies函數(shù)φ(x)的一階矩.尺度函數(shù)φi(x)不僅滿足正交性,同時還近似滿足位移內(nèi)插特性,即φi(x)≈δi,0,δi,i′為狄拉克函數(shù).

簡單而不失一般性,以電場分量Ex為例,對Maxwell離散迭代格式進行簡單說明.將(6)式入Maxwell旋度方程,可得

采用Galerkin法則,在方程兩側(cè)同時乘以φi+1/2(x)φj(y)φk(z)hn+1/2(t),并求取對應項的內(nèi)積,結(jié)合函數(shù)hn(t)和φi(x)的正交性,(15)式可化簡為

式中,a(l)稱為連接系數(shù),l的變化范圍為-Ls—Ls-1,Ls=2N-1稱為基函數(shù)的有效支撐尺寸,N為小波函數(shù)消失矩,文中取N=2.a(l)可由尺度函數(shù)在傅里葉頻域的內(nèi)積計算得到：

式中,m=(i+1/2,j,k)為場分量的網(wǎng)格離散坐標,CA和CB的計算方法如下：

由于函數(shù)φi(x)的位移內(nèi)插特性,空間r0的電磁場值就等于該處的場展開系數(shù),也即有成立,因此可避免對場分量的重構(gòu),這是采用Daubechies尺度函數(shù)作為基函數(shù)展開的優(yōu)勢所在.

同理,電場分量Ey,Ez和磁場分量的MRTD迭代步進方程也有類似形式.由于MRTD是Maxwell方程的一種顯式差分格式,因此為保證迭代過程的收斂性,時間差分間隔Δt與網(wǎng)格差分步長必須滿足Courant穩(wěn)定性條件,其中對于基于Daubechies尺度函數(shù)的MRTD算法,時間迭代步長Δt必須滿足以下關(guān)系式[31]：

式中,c為光速;Δx,Δy和Δz分別為x,y和z方向的離散步長.

3.3 ADE-PML吸收邊界

采用基于輔助微分方程的PML邊界(perfectly matched layer with auxiliary differential equation,ADE-PML)作為吸收邊界.本文在劉亞文等[32]提出的二維ADE-PML邊界基礎(chǔ)上,將其推廣至三維情形.下面仍以Ex分量為例,介紹吸收邊界的處理過程.在頻域情形下,坐標伸縮的Maxwell旋度方程組如下：

求取si的倒數(shù)可得其 中和為非負實數(shù),且κi≥1[32].將其代入(21)式Ex分量的計算公式中,并引進輔助變量ψexy和ψexz可得：

進一步將上式改寫為時域形式：

參照3.2節(jié)的處理方法,采用Galerkin法對三個偏微分方程進行離散,可得場分量的迭代步進方程：

以上各式中,ai和bi(i=x,y,z)滿足以下公式：

ADE-PML吸收邊界通常采用PEC邊界進行截斷.在邊界處,由于MRTD迭代過程涉及邊界外的場分量,因此為保證邊界附近場分量的正常迭代,需要采用鏡像原理對邊界外的場分量進行拓展.

3.4 MRTD的總場-散射場技術(shù)

與FDTD模型類似,采用總場散射場技術(shù)引入平面入射波.在該方案中,計算域被分為總場區(qū)和散射場區(qū),總場區(qū)電磁分量為散射場Esca和Hsca與入射場Einc和Hinc之和(如(32)式),散射場區(qū)僅包含散射場Esca和Hsca[33].

由MRTD的步進迭代方程可知,在計算連接邊界附近節(jié)點的場分量時,需要同時使用散射場和總場區(qū)的切向分量,因此需要通過引入入射波對相應計算節(jié)點進行訂正,考慮各場分量迭代公式的相似性,下面僅以電場分量Ex為例進行簡單說明.

圖2所示為Ex分量的總場和散射場節(jié)點訂正區(qū)域,其中圖2(a)為散射場需要進行訂正的節(jié)點區(qū)域,共4個;圖2(b)為總場區(qū)需要訂正的節(jié)點區(qū)域,共8個.考慮各區(qū)域訂正方法的基本一致,僅以個散射場區(qū)域1和總場區(qū)域5作為例子進行解釋.

區(qū)域1的范圍為：i∈[Imin,Imax-1],j∈[Jmin,Jmax],k∈[Kmax+1,Kmax+Ls-1].由于該區(qū)域為散射場區(qū),在迭代計算過程中,所涉及的總場節(jié)點需要扣除入射場的值,因此,該區(qū)域內(nèi)各元胞Ex的計算方法如下：

圖2 Ex分量的總場和散射場節(jié)點訂正區(qū)域Fig.2.Correction area for the Yee cell distributed in scattering field and toal field.

以上各式中,下標s,tol和inc分別表示該場分量為散射場、總場和入射場分量.連接邊界附近各節(jié)點的入射波切向分量采用一維MRTD投影至三維總場邊界的方案計算.該方法基本流程與FDTD相似,其中一維MRTD和三維MRTD的時間步長必須相同.

3.5 基于體積積分的遠場計算模型

氣溶膠和冰晶一般為非導磁電介質(zhì),因此在介質(zhì)內(nèi),存在電極化矢量P(r),而磁極化矢量M(r)=0.取電磁波隨時間的變化因子為exp(jωt),則在頻域中,介質(zhì)中的Maxwell方程可改寫為[34]

上式中P(r)=χE(r)=(ε(r)/ε0-1)E(r),χ為電極化率.

對(35)式第一個方程兩端取旋度,并將第二個方程代入可得

采用矢量恒等式?×?×E=?(?·E)-ΔE進一步化簡可得到

式中Δ為拉普拉斯算子,由數(shù)學知識,該方程的格林函數(shù)形式解應滿足(38)式和(39)式.

其中格林函數(shù)的具體形式如下：

在遠場,取近似|r-r′|=r-r′·er,并將上式寫成頻域形式,可得散射場的表達形式如下：

由于MRTD計算區(qū)域是離散的,Esca需要采用數(shù)值計算方案求解,本文采用以下方法計算：

式中(θob,φob)為觀測方位角;為Yee元胞的平均介電常數(shù);和為電場分量在坐標(i′,j′,k′)的平均值(采用相鄰兩個Yee元胞平均);k=ker=(kx,ky,kz)為波矢量.Ts為中間變量,如下式所示：

進一步將直角坐標系下的場分量變換至球坐標系下,可得

3.6 粒子散射特性的計算

3.6.1 散射振幅矩陣與穆勒矩陣的計算

通過(50)式便可獲得散射振幅矩陣,在此基礎(chǔ)上,通過穆勒矩陣F與散射振幅矩陣S的轉(zhuǎn)換關(guān)系,便可計算穆勒矩陣,詳細請參見文獻[15,30].

3.6.2 吸收截面的推導

吸收截面表征的是粒子吸收電磁波的能力,實質(zhì)是損耗電磁能而將其轉(zhuǎn)化為熱能.記電磁波損耗功率功率為Pab,則Pab可采用下式表示：

上式中,ns為積分面的外法向矢量,S(r)為電磁場的坡印廷矢量.采用矢量恒等式?·(E(r)×H?(r))=H?(r)·?×E(r)-E(r)·?×H?(r)對體積積分部分進行化簡可得

代入(51)式,可得電磁波損耗功率Pab的形式為

上式表明粒子吸收的能量轉(zhuǎn)化為焦耳熱,與吸收截面的定義相符.進一步,入射場的坡印廷矢量Sinc可表示為

式中Z0為波阻抗,為歸一化的入射波波矢量.根據(jù)吸收截面的定義,粒子吸收截面Cab可表示為損耗功率Pab與入射功率Pinc之比,可采用下式計算：

3.6.3 消光截面的推導

頻域條件下,總場的坡印廷矢量可分解為以下三部分：入射坡印廷矢量Sinc,散射坡印廷矢量Ssca和消光坡印廷矢量Sext[30].其中消光坡印廷矢量Sext可表示為[30]

對坡印廷矢量Sext曲面積分并取實部計算粒子消光功率,如下式所示：采用矢量恒等式?·(A×B)=B·?×A-

A·?×B,結(jié)合麥克斯韋方程組對上式進行化簡可以得到：

對上式進一步化簡,可以得到

在獲取消光功率的基礎(chǔ)上,采用入射波功率對其歸一化,即可得到消光截面Cext,具體計算方法如下：

在獲得消光截面Cext和吸收截面Cab的基礎(chǔ)上,求取兩者之差便可得到粒子的散射截面Csc.進一步,采用粒子幾何截面對粒子消光截面、吸收截面及散射截面進行歸一化,即可得到消光、吸收和散射效率因子,具體公式請參見文獻[4].

4 模型驗證與性能分析

基于FORTRAN實現(xiàn)了MRTD散射模型的編碼,下面分別將該模型模擬的穆勒矩陣與積分散射參數(shù)(消光、吸收及散射效率因子)和Mie理論、T矩陣法進行對比以驗證模型的準確性,并對模型的性能進行定量評估.

4.1 與Mie散射理論的對比驗證

4.1.1 穆勒矩陣的對比驗證

取入射光波長λ為0.55μm,粒子半徑為0.5μm,復折射率為1.53-0.008i,設(shè)置離散網(wǎng)格步長為λ/50,采用MRTD模擬粒子的穆勒矩陣,并將其與Mie理論的模擬結(jié)果進行對比,結(jié)果如圖3所示.左圖為穆勒矩陣隨散射角的變化,其中F12,F34和F44采用F11進行了歸一化處理;右側(cè)為兩模型的相對偏差分布,由于F12和F34在散射角0?和180?附近值趨于,因而相對誤差極大,無實際意義,因此僅給出散射角20?—170?區(qū)間的相對誤差.

圖3 (網(wǎng)刊彩色)MRTD與Mie散射模型模擬的穆勒矩陣及相對誤差分布Fig.3.(color online)Comparison of Müeller matrix between MRTD and Mie scattering model and the distribution of relative simulation error.

如圖3所示,MRTD模型模擬結(jié)果與Mie散射模型一致性較好,兩者曲線基本重合.對于散射相函數(shù)F11,兩模型的相對誤差在7%以內(nèi),其中在散射角0?—120?區(qū)間,兩者相對誤差小于3%;對于F12/F11和F34/F11,其相對誤差隨散射角波動顯著,但整體而言,其相對誤差在10%以內(nèi),F44/F11的模擬相對誤差小于F12/F11和F34/F11,其中前向散射方向的相對誤差在3.5%以內(nèi).從相對誤差分布上看,前向散射方向的相對誤差整體小于后向散射方向,其原因是由于采用網(wǎng)格對光滑球形粒子近似存在“階梯誤差”,即形狀存在差異,而這種形狀的差異對后向散射的影響強于前向散射[36].

4.1.2 散射、消光及吸收效率因子的對比驗證

入射光波長、粒子復折射率及網(wǎng)格設(shè)置與穆勒矩陣驗證實驗相同,在不同粒徑條件下采用MRTD模擬粒子的消光、吸收和散射效率因子,并將其與Mie理論結(jié)果對比,結(jié)果列于表1.

表1 MRTD與Mie理論模擬的消光、吸收和散射效率因子結(jié)果對比Table 1.Extinction,absorption and scattering efficiency simulated by MRTD and Mie scattering model.

由表1可知,不同粒徑條件下,MRTD與Mie理論模擬的消光、吸收和散射效率因子相對偏差均在6%以內(nèi),驗證了兩模型模擬結(jié)果的一致性.MRTD與Mie理論模擬結(jié)果的相對誤差多數(shù)為負值,說明該MRTD模擬結(jié)果整體小于Mie散射理論的模擬值;隨著粒子半徑的增大,消光、吸收和散射效率因子三個參數(shù)的模擬誤差整體呈先減小后增大的特征,其中當粒徑和入射光波長相當時,MRTD模型的計算精度最高,三個參數(shù)的模擬誤差均小于0.5%.

4.2 與T矩陣法的對比驗證

T矩陣法是計算非球形粒子散射特性的有力工具,其準確性已得到廣泛驗證.本文將MRTD散射模型與T矩陣法對比,進一步在非球形條件下驗證MRTD模型的準確性,其中模擬的粒子形狀包括橢球和圓柱粒子.如圖4所示,橢球形狀采用橫縱軸比a/b表示,圓柱粒子采用直徑高度比D/L表示,粒子的空間取向采用旋轉(zhuǎn)軸(rotational axis)的天頂角α和方位角β表示.

圖4 (網(wǎng)刊彩色)橢球及圓柱粒子的參數(shù)說明Fig.4.(color online)The illustrations of the parameters for spheriods and cylinders.

4.2.1 橢球粒子的情形

4.2.1.1 穆勒矩陣的對比驗證

取入射光波長λ為0.55μm,計算域離散網(wǎng)格間距為λ/5,設(shè)置橢球粒子的等效球形半徑(等體積法)為0.5μm,橫縱軸比a/b=0.5,粒子取向(α,β)=(0,0),復折射率為1.53-0.008i,采用MRTD模型與T矩陣法分別計算粒子的穆勒矩陣并進行對比,結(jié)果如圖5所示.

由圖5可知,MRTD散射模型與T矩陣法模擬的穆勒矩陣一致性較好,其中對于F11,在前向散射方向,兩模型的相對偏差小于3.3%;在后向散射方向,平均相對偏差僅為4.6%.與球形粒子結(jié)果類似,對于F12/F11,F34/F11和F12/F11,散射角0?—90?區(qū)間的模擬誤差明顯小于90?—180?方向,其中在側(cè)向散射方向(120?—140?),由于曲線變化劇烈,兩模型計算值相差略大,但整體變化趨勢是一致的,最大相對誤差為9.7%,這可能是由于對粒子形狀建模存在誤差所致.

4.2.1.2 散射、消光及吸收效率因子的對比驗證

入射光、網(wǎng)格粗細、粒子復折射率及取向的設(shè)置與上面穆勒矩陣驗證實驗相同;取粒子橫縱軸比a/b取0.6,等效粒徑為0.1,0.5,1.0和2.0μm,分別采用MRTD散射模型與T矩陣法模擬橢球粒子消光效率因子、吸收效率因子和散射效率因子,并求取兩者相對偏差,結(jié)果如表2所列,其中,粒子的幾何截面采用等效球形半徑計算.

由表2可知,橢球粒子情形下,MRTD與T矩陣法模擬的粒子積分散射參數(shù)誤差較小,其中散射效率因子和消光效率因子的誤差小于5%,吸收效率因子誤差小于7.5%.與球形粒子相似,當粒子尺度與入射光波長相當時,散射參數(shù)的模擬精度最高,其中當粒徑為0.5μm時,消光效率因子的模擬誤差小于0.01%,吸收效率因子的模擬誤差優(yōu)于0.25%.從誤差分布上看,MRTD計算的橢球粒子的散射參數(shù)略小于T矩陣法,說明MRTD模型可能略低估實際粒子的散射吸收能力.

圖5 (網(wǎng)刊彩色)MRTD與T矩陣模型模擬的橢球粒子穆勒矩陣Fig.5.(color online)Comparison of Müeller matrix of spheriod simulated by MRTD and T matrix method.

表2 MRTD與T矩陣模型模擬的消光、吸收和散射效率因子結(jié)果對比Table 2.Extinction,absorption and scattering efficiency simulated by MRTD and T matrix method.

4.2.2 圓柱粒子情形

4.2.2.1 穆勒矩陣的對比驗證

在圓柱粒子情形下,比較MRTD和T矩陣法模擬的穆勒矩陣,結(jié)果如圖6所示.模擬過程中,入射光波長、網(wǎng)格粗細及粒子取向設(shè)置與橢球粒子情形相同,粒子復折射率為1.53-0.008i,直徑高度比為D/L=1.0.

由圖6可知,圓柱粒子情形下,MRTD模擬的穆勒矩陣變化曲線與T矩陣結(jié)果一致性程度較高,其中對于F11,在散射角0?—90?區(qū)間,兩者的相對偏差在3%以內(nèi),在后向散射方向,兩者最大相對誤差僅為7.6%.對于F12/F11,F34/F11和F44/F11,在0?—140?散射角區(qū)間,MRTD與T矩陣法模擬結(jié)果基本重合,最大相對偏差小于10%,在F12/F11和F34/F11隨散射角劇烈變化的區(qū)間,兩者的一致性程度稍差,特別是在散射角160?—170?附近.整體而言,在前向散射方向,兩模型的一致性程度高于后向散射方向.

圖6 (網(wǎng)刊彩色)MRTD與T矩陣模型模擬的圓柱形粒子穆勒矩陣Fig.6.(color online)Comparison of Müeller matrix of cylinder simulated by MRTD and T matrix method.

4.2.2.2 散射、消光及吸收效率因子的對比驗證

參數(shù)設(shè)置與穆勒矩陣驗證實驗相同,在不同等效球形粒徑(等體積法)條件下,分別采用MRTD模型與T矩陣法模擬圓柱粒子消光、吸收和散射效率因子,并計算兩者相對偏差,結(jié)果如表3所列,與橢球粒子情形類似,粒子幾何截面采用等效球形半徑計算.

由表3可知,MRTD模擬的粒子積分散射特性與T矩陣法模擬值相符,其中消光效率因子的模擬相對偏差小于5.5%,吸收效率因子最大偏差僅為8.745%.相比而言,吸收效率因子的相對偏差略大于消光和散射效率因子,其原因是由于吸收截面值相對較小所致.與球形和橢球粒子情形類似,當粒子等效球形半徑與入射光波長相當時,粒子散射特性模擬的精度最高,隨著粒子半徑偏離該尺寸,模擬誤差隨之增大.從誤差符號分布上看,MRTD模型整體上也是略微低估了圓柱形粒子的散射吸收特性.

表3 MRTD與T矩陣理論模擬的消光、吸收和散射效率因子結(jié)果對比Table 3.Extinction,absorption and scattering efficiency simulated by MRTD and T matrix method.

4.3 網(wǎng)格粗細對模擬精度的影響

網(wǎng)格粗細決定了粒子形狀描述的精度,也決定了粒子散射特性模擬的準確性.網(wǎng)格精細,粒子形狀特征描述越好,但同時計算量與內(nèi)存消耗隨之迅速增長,因此針對不同尺度粒子選擇適當粗細的網(wǎng)格,對緩解平衡計算量與精度的矛盾是十分有必要的.基于此,本節(jié)主要討論網(wǎng)格粗細對不同尺度粒子散射特性模擬精度的影響.

取入射光波長λ為0.55μm,粒子形狀為球形,復折射率為1.33—0.010i,分別取粒子半徑為0.1,0.5,1.0和2.0μm,設(shè)置網(wǎng)格間距δ為λ/10,λ/20,λ/30,λ/40,λ/50和λ/60,分別采用MRTD和Mie散射模型模擬粒子的消光、吸收效率因子及穆勒矩陣.為定量評估MRTD模型的模擬精度,消光和吸收效率因子的精度采用MRTD與Mie散射模型模擬結(jié)果的相對偏差表征;散射光空間分布特性的模擬精度采用F11的均方根誤差RSE表征,其定義如下：

計算得到不同粒徑條件下MRTD模型與Mie散射模型的散射參數(shù)模擬相對偏差如表4所列.

由表4可知,隨著計算網(wǎng)格的細化,粒子散射特性的模擬精度迅速提升,其中對于粒徑為1.0μm的粒子,當δ從λ/10變化至λ/60時,Qex和Qab的模擬誤差分別由13.5122%和10.2059%減小至1.0634%和1.1583%.當粒子尺度參數(shù)小于22時,相同網(wǎng)格粗細條件下,粒子半徑與入射光波長越相近,MRTD方法模擬的粒子散射特性精度越高;對于尺度參數(shù)x約為1,5,10和20的粒子,空間離散網(wǎng)格步長分別取λ/60,λ/30,λ/40和λ/40即可滿足精度需求.

4.4 模型運行效率分析

進一步分析MRTD模型的運行效率.所用計算機的CPU是Intel(R)Core i5,主頻3.1 GHz,內(nèi)存8 GB,軟件系統(tǒng)為Win7操作系統(tǒng)和Fortran 90 Compiler.根據(jù)4.3節(jié)中確定的保精度網(wǎng)格間隔,分別模擬大折射率實部、大折射率虛部及一般復折射率三種情形的粒子散射特性,并統(tǒng)計所需運行時間,如表5所列.

由表5可知,隨著粒徑的增大,相應的散射特性模擬時間迅速增長,其中當粒徑為0.1μm時,運行時間僅僅為285 s,當粒徑為3.0μm時,模擬時間高達2.8855×105s;折射率實部和虛部的變化并不會顯著地改變程序運行的時間.

表4 計算域網(wǎng)格間距對粒子散射特性模擬精度的影響Table 4.The influence of grid size on the simulation accuracy of aerosol scattering properties.

表5 不同粒徑粒子散射特性模擬所需的時間統(tǒng)計Table 5.Computational time for particles with different radius and space resolution.

5 結(jié) 論

非球形氣溶膠的散射特性是影響輻射傳輸過程的重要因素.為提高復雜非球形氣溶膠的散射特性模擬能力,將計算電磁學中的MRTD技術(shù)引入至粒子近場散射場的計算,推導了僅基于電場分量體積積分的近遠場外推方案,建立了基于體積積分的粒子吸收截面、消光截面等散射特性計算方案,將MRTD散射模型與Mie散射模型、T矩陣法的模擬結(jié)果進行了對比驗證,并分析網(wǎng)格大小對模擬精度的影響.主要結(jié)論如下：

1)MRTD散射模型與Mie理論、T矩陣法的模擬結(jié)果一致性較好,其中相函數(shù)的最大相對誤差小于8%,前向散射方向小于4%;當粒徑與入射光波長相當時,消光和散射效率因子的相對誤差小于0.1%;

2)空間網(wǎng)格粗細對模擬精度影響顯著,當粒子尺度參數(shù)小于20時,在相同精度要求下,所需網(wǎng)格間距隨粒子尺度參數(shù)呈先減小后增加的特征,其中當入射光波長與粒子半徑相當時,需要的Yee元胞數(shù)量最少;

3)與FDTD及PSTD相似,本模型的運行時間也隨著粒子半徑的增大而迅速增加,下一步需要對模型進行并行化處理以提高模型的運行效率.

感謝解放軍理工大學野戰(zhàn)工程學院電磁環(huán)境效應與電光工程國家級重點實驗室陳彬教授及劉亞文博士在電磁場理論方面提供的指導.

參考文獻

[1]Dou T,Xiao C,Shindell D T,Liu J,Ming J,Qin D 2012Atmos.Chem.Phys.12 7995

[2]Liou K N,Takano Y 1994Atmos.Res.31 271

[3]Intergovernmental Panel of Global Climate Change 2007IPCC:Climate Change

[4]Liou K N 2003An Introduction to Atmospheric Radiation(San Diego：Academic Press)

[5]Liou K N,Takano Y,Yang P 2013J.Quantit.Spectrosc.Radiat.Transfer127 149

[6]Rao R Z 2012Modern Optics(Beijing：Scientific Express)p31[饒瑞中2012現(xiàn)代大氣光學(北京：科學出版社)第31頁]

[7]Hu S,Gao T C,Li H,Liu L,Liu X C,Zhang T,Cheng T J,Li W T,Dai Z H,Su X J 2016J.Geophys.Res.121 doi：101002/2015JD024105

[8]Evans K F,Stephens G L 1991J.Quant.Spectrosc.Radiat.Transfer46 413

[9]Evans K F 1998J.Atmos.Sci.55 429

[10]Han Y,Wang T J,Rao R Z,Wang Y J 2008Acta Phys.Sin.57 7396(in Chinese)[韓永,王體健,饒瑞中,王英儉2008物理學報57 7396]

[11]Hu S,Gao T C,Liu L 2014J.Meteorolog.Sci.34 612(in Chinese)[胡帥,高太長,劉磊2014氣象科學34 612]

[12]Zhang X L,Huang Y B,Rao R Z 2013High Power Laser and Particle Beams25 1675(in Chinese)[張小林,黃印博,饒瑞中2013強激光與粒子束25 1675]

[13]Herman M,Deuzé J L,Marchand A,Roger B,Lallart P 2005J.Geophys.Res.110 D10S02

[14]Yang P,Liou K N,Bi L,Liu C,Yi B,Baum B A 2015Adv.Atmos.Sci.32 32

[15]Mishchenko M I,Hovenier J W,Travis L D 2000Light Scattering by Nonspherical Particles,Thoery,Measurements,and Application(New York：Academic Press)

[16]Yang P,Liou K N 1995J.Opt.Soc.Am.A12 12

[17]Mishchenko M I,Travis L D 1998J.Quant.Spectrosc.Radiat.Transfer60 309

[18]Bi L,Yang P 2016J.Quant.Spectrosc.Radiat.Transfer178 93

[19]Voshchinnikov N V,Farafonov V G 1993Astrophys.Space Sci.204 19

[20]Al-Rizzo H M,Tranquilla J M 1995J.Computat.Phys.119 356

[21]Harrington R F 1968Field Computation by Moment Methods(New York：Macmillan)

[22]Draine B T 1988Astrophys.J.333 848

[23]Draine B T,Flatau P J 1994J.Opt.Soc.Am.A11 1491

[24]Morgan M A,Mei K K 1979IEEE Trans.Antenn.Propagat.27 202

[25]Liu C,Panetta R L,Yang P 2013J.Quant.Spectrosc.Radiat.Transfer129 169

[26]Liu C,Panetta R L,Yang P 2012J.Quant.Spectrosc.Radiat.Transfer113 1728

[27]Liu Y,Chen Y,Zhang P 2013Prog.Electromagn.Res.143 223

[28]Cheong Y W,Lee Y M,Ra K H,Kang J G,Shin C C 1999IEEE Microw.Guided Wave Lett.9 297

[29]Dai S Y,Wu Z S 2008Acta Phys.Sin.57 7636(in Chinese)[代少玉,吳振森 2008物理學報 57 7636]

[30]Bohren C F,Huffman D R 1983Absorption and Scattering of Light by Small Particles(New York：John Wiley&Sons)

[31]Tentzeris E M,Cangellaris A,Katehi L P B,Harvey J 2003IEEE Trans.Microw.Tech.50 501

[32]Liu Y W,Chen Y W,Xu X,Liu Z X 2013Acta Phys.Sin.62 034101(in Chinese)[劉亞文,陳亦望,徐鑫,劉宗信2013物理學報62 034101]

[33]Gao Q,Cao Q,Zhou J 2009International Forum on Information Technology and Applications2009 359-362

[34]Griffiths D J 2014Introduction to Electrodynamics(third edtion)(New York：Pearson Education)

[35]van der Hulst H C 1981Light Scattering by Small Particles(New York：Dover Publications)

[36]Curtis D B,Meland B,Aycibin M 2008J.Geophys.Res.113 D08210

PACS:42.68.Mj,42.68.Ay,95.30.Jx,42.25.Bs DOI:10.7498/aps.66.044207

Simulating scattering properties of nonspherical aerosol particles using multiresolution timedomain method?

Hu Shuai1)Gao Tai-Chang1)2)?Li Hao1)Yang Bo2)Jiang Zhi-Dong3)Chen Ming1)Li Shu-Lei1)

1)(College of Meteorology and Oceanography,PLA University of Science and Technology,Nanjing 211101,China)
2)(National Key Laboratory on Electromagnetic Environment and Electro-optical Engineering,PLA University of Science and Technology,Nanjing 210007,China)
3)(Navy Aeronautical and Astroautical University,Qingdao 266041,China)

19 September 2016;revised manuscript

30 November 2016)

Scattering process of aerosol particles plays an important role in atmospheric radiative transfer since it can modify the transmission,reflection and absorption ability of atmospheric system.Owning to the uncertainty of aerosol particles’scattering properties,which results from their complicated geometries and inhomogeneous compositions,there still exists a considerable uncertainty in the radiative transfer numerical simulation,and simulating the scattering properties of aerosol with irregular shapes has become a hotspot in meteorological study.To this end,a new aerosol scattering model is developed based on multi-resolution time-domain(MRTD),by which the scattering processes of nonspherical and inhomogeneous particles can be simulated.In this model,the near electromagnetic field is calculated by MRTD technique.Considering the particularity of aerosol medium,a transformation technique from near field to far field is derived based on volume integration method,and then the scattering amplitude matrix and Müeller matrix can be calculated by the obtained far electric field as well.The models for particle extinction and absorption cross section are derived from Maxwell’s curl equations in the frequency domain,by which the integration scattering properties can be simulated accurately.The MRTD scattering model is validated by comparing with Mie theory and T matrix method for spherical particle,ellipsoidal particle and cylindrical particle,and the influence of grid size on the simulation accuracy is analyzed subsequently.In the last part,the efficiency of the MRTD scattering model is quantitatively discussed.The simulation results show that the relative errors of scattering phase function simulated by our model are less than 8%,and the errors in forward scattering direction are much smaller,which are less than 4%.The precisions for extinction and absorption efficiency are much higher than the results from the scattering phase function,and the relative errors can reduce to 0.1%for particles with their radii comparable to the wavelength of incident light.The gird size has a significant influence on model precision;to achieve the same accuracy,the grid size first increases with increasing particle radius,and then decreases as a function of particle size for particles with size parameter less than 20.In the next step,we will try to establish the scattering property database of nonspherical particles based on the MRTD scattering model developed here.

nonspherical aerosol,scattering properties,multi-resolution time-domain,Müeller matrix

:42.68.Mj,42.68.Ay,95.30.Jx,42.25.Bs

10.7498/aps.66.044207

?國家自然科學基金(批準號：41575025,41575024)資助的課題.

?通信作者.E-mail：gaotc@gmail.com

*Project supported by the National Natural Science Foundation of China(Grant Nos.41575025,41575024).

?Corresponding author.E-mail：gaotc@gmail.com