降雨粒子群散射特性的高效分析

丁大志 沈 鵬 陳如山 尚 社 范曉彥

（1.南京理工大學(xué)通信工程系，江蘇 南京 210094；2.空間微波技術(shù)國(guó)防科技重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710000）

引 言

全球降雨分布及其測(cè)量是目前全球環(huán)境氣候變化研究中的薄弱環(huán)節(jié)，當(dāng)前只有少數(shù)國(guó)家可以利用熱帶降雨測(cè)量（TRMM）衛(wèi)星獲得可靠的全球降雨分布數(shù)據(jù)［1-2］。TRMM降雨雷達(dá)工作在Ku波段，中心頻率為13.8GHz，可檢測(cè)的最小降雨量約為0.6mm／h.為了進(jìn)一步提高星載降雨雷達(dá)對(duì)小雨、雪等氣候現(xiàn)象的測(cè)量靈敏度，新一代星載降雨雷達(dá)將增加K和Ka工作頻段［3］，使得降雨雷達(dá)能夠?qū)π∮暄┑葰夂蜻M(jìn)行更加精確的測(cè)量。因此，對(duì)K和Ka波段降雨粒子散射特性及不同頻率對(duì)降雨粒子散射特性的影響進(jìn)行研究對(duì)改進(jìn)天氣和洪水預(yù)報(bào)、水資源管理以及監(jiān)測(cè)全球變暖與氣候變化具有重要意義。

傳統(tǒng)分析降雨粒子群散射的方法可以概括為兩大類:解析解法和數(shù)值近似解法。其中，解析解法主要是分離變量（SVD）法，最早由 Asano等人提出［4］，但是該方法對(duì)坐標(biāo)系的依賴性限制了其對(duì)復(fù)雜形狀粒子散射特性的研究。隨著計(jì)算機(jī)硬件環(huán)境的發(fā)展，數(shù)值類方法得到快速發(fā)展。用于解決粒子散射問(wèn)題的計(jì)算方法主要包括離散偶極子近似（DDA）法、T矩陣（T-Matrix）法、基于體積分方程（VIE）的矩量法、稀疏矩陣規(guī)則網(wǎng)格（SMCG）法等。DDA方法最早由Purcell和Pennypacker在1973年提出，用于分析星際塵埃物質(zhì)的散射和吸收問(wèn)題［5］。但是其對(duì)存儲(chǔ)量、計(jì)算時(shí)間的高需求使其發(fā)展受到很大限制。T矩陣方法最早由P.C.Waterman［6］于1965年首次引入計(jì)算電磁學(xué)領(lǐng)域，最顯著的特點(diǎn)是它只與粒子自身的屬性有關(guān)，對(duì)于某一固定的散射粒子，計(jì)算一次它的T矩陣，便可以得到任意方位下該散射體的散射特征，但是在分析非球?qū)ΨQ類模型時(shí)，其計(jì)算復(fù)雜度會(huì)大幅度增加。利用體積分方程（VIE）結(jié)合矩量法［7］求解宏觀粒子的散射問(wèn)題，可以精確地分析目標(biāo)的電磁特性，但是其未知量大、計(jì)算時(shí)間長(zhǎng)等缺點(diǎn)使其只能對(duì)較小點(diǎn)尺寸的目標(biāo)進(jìn)行分析。SMCG方法［8-9］將傳統(tǒng)矩量法中的阻抗矩陣分為強(qiáng)相互作用矩陣和弱相互作用矩陣，對(duì)弱相互作用矩陣進(jìn)行Taylor級(jí)數(shù)展開(kāi)，利用FFT變換計(jì)算弱相互作用矩陣與待求向量的乘積，降低了內(nèi)存需求和計(jì)算復(fù)雜度。

簡(jiǎn)要闡述了電場(chǎng)積分方程和SMCG方法的基本原理。對(duì)K和Ka波段下均勻分布的介質(zhì)降雨粒子的散射特性進(jìn)行了對(duì)比分析，為了更好地模擬實(shí)際降雨過(guò)程，分析了隨機(jī)分布的降雨粒子群的散射特性，研究了一定區(qū)域粒徑隨機(jī)分布的降雨雨團(tuán)的散射特性。

1.SMCG方法基本理論

1.1 電場(chǎng)積分方程

對(duì)于介質(zhì)體散射問(wèn)題，通?？梢酝ㄟ^(guò)求解表面積分方程和體積分方程來(lái)解決。但是，由于降雨粒子為非均勻的介質(zhì)散射體，一般選擇體積分方程求解散射場(chǎng)。主要針對(duì)體積分方程進(jìn)行介紹。

假設(shè)散射體的介電常數(shù)和電導(dǎo)率分別為ε（r）、σ（r），r為場(chǎng)點(diǎn)，r′為源點(diǎn)入射平面電磁波在散射體上產(chǎn)生的等效感應(yīng)電流Jeq可表示為［10］

式中:E表示介質(zhì)散射體上的總電場(chǎng)強(qiáng)度。此等效電流在背景均勻的介質(zhì)中產(chǎn)生的散射場(chǎng)可以表示為

式中:k20=ω2μ0ε0，且μ0為自由空間中的磁導(dǎo)率，pv表示主值積分。

由于介質(zhì)體的總場(chǎng)強(qiáng)度可以表示為入射場(chǎng)和散射場(chǎng)總和的形式，即

將式（2）代入式（5）可得

上式便是介質(zhì)散射體的體積分方程表達(dá)式。對(duì)于脈沖基點(diǎn)匹配方法，E（r′）·=G（r，r′）可表示為

簡(jiǎn)化起見(jiàn)，令x1=x x2=y x3=z，則

式（7）中的（r，r′）是對(duì)稱矩陣，可以將式（6）表示為

假設(shè)介質(zhì)散射體被剖分為n個(gè)小體積單元，每個(gè)介質(zhì)體單元內(nèi)的E（r）和τ（r）都是常量。第m個(gè)體積元的體積表示為Vm，rm表示在Vm內(nèi)選取的任意點(diǎn)。在每個(gè)點(diǎn)rm處，式（6）表示成矩陣方程為

則式（10）可以表示為

令［Gxpxq］為N×N階的矩陣，其中的元素可以定義為

［Exp］和［Eixp］均為N維向量，表示成如下形式

式（6）可以表示成如下矩陣形式

將上式簡(jiǎn)化表示為

式中:［G］是3 N×3 N 階的矩陣，［E］和［Ei］為3 N維的向量。通過(guò)求解上式可以確定N個(gè)選定點(diǎn)處的電場(chǎng)。

1.2 SMCG方法

SMCG方法的主要思想是根據(jù)剖分單元之間相互作用的強(qiáng)弱將阻抗矩陣劃分成強(qiáng)相互作用矩陣［GS］和弱相互矩陣［GW］

式中:［GS］為稀疏矩陣；［GW］是稠密矩陣。對(duì)于［GS］矩陣可以采用壓縮存儲(chǔ)的形式節(jié)省內(nèi)存，而弱相互矩陣矢量乘的計(jì)算及存儲(chǔ)耗費(fèi)了大部分的CPU時(shí)間和內(nèi)存。為加速弱相互矩陣矢量乘的計(jì)算，通常采用SMCG快速算法來(lái)加速計(jì)算。SMCG方法是將Green函數(shù)在規(guī)則網(wǎng)格點(diǎn)上進(jìn)行Taylor級(jí)數(shù)展開(kāi)，則弱相互矩陣［GW］可以進(jìn)一步表示為

式中:K為Taylor級(jí)數(shù)總的展開(kāi)項(xiàng)數(shù)。則式（16）相應(yīng)的零階和高階解的迭代計(jì)算式為

圖1 弱相互作用矩陣元素的直接和間接積分路徑

圖1從二維角度示意了兩種積分路徑，三維情況可以類似推導(dǎo)出。假定介質(zhì)單元rm=（xm，ym，zm）與rn=（xn，yn，zn）為弱作用的兩個(gè)單元，則其阻抗元素G（rm，rn）的積分路徑有兩種選擇:其一是從點(diǎn)rm到點(diǎn)rn的直接路徑；其二是從點(diǎn)rm到其對(duì)應(yīng)規(guī)則網(wǎng)格點(diǎn)rmo=（xmo，ymo，zmo），再?gòu)狞c(diǎn)rno到與點(diǎn)rn對(duì)應(yīng)規(guī)則網(wǎng)格點(diǎn)rno= （xno，yno，zno），最后，從點(diǎn)rno到點(diǎn)rno.間接計(jì)算可以通過(guò)格林函數(shù)關(guān)于規(guī)則網(wǎng)格點(diǎn)的Taylor級(jí)數(shù)展開(kāi)實(shí)現(xiàn)，即

式中:xdm= （xm-xmo）；xdn= （xno-xn）；xmn=（xmo-xno），同理，對(duì)應(yīng)的ydm、ydk、ymn、zdm、zdn、zmn有類似的表示。弱相互作用矩陣與待求向量的乘積可以寫成

式中:［Tsi］相當(dāng)于將基元素中心點(diǎn)移到與其最近的規(guī)則網(wǎng)格點(diǎn)上；［Gi］為Toeplitz矩陣，該矩陣元素反映了規(guī)則網(wǎng)格點(diǎn)間的相互作用；矩陣［Tti］相當(dāng)于將規(guī)則網(wǎng)格點(diǎn)移回到檢測(cè)中心點(diǎn)。Toeplitz矩陣的行之間滿足循環(huán)平移關(guān)系，進(jìn)行矩陣矢量乘時(shí)可通過(guò)快速傅里葉變換（FFT）實(shí)現(xiàn)，SMCG算法能夠?qū)?nèi)存需求從O（N2）降至O（N），計(jì)算復(fù)雜度從O（N2）降至O（NlogN）.

在對(duì)弱相互作用矩陣元素進(jìn)行Taylor級(jí)數(shù)展開(kāi)時(shí)，需要對(duì)近場(chǎng)作用與遠(yuǎn)場(chǎng)作用的臨界距離rmax和展開(kāi)級(jí)數(shù)所取階數(shù)（對(duì)應(yīng)項(xiàng)數(shù)）進(jìn)行折中選取。rmax選取過(guò)小、級(jí)數(shù)展開(kāi)項(xiàng)過(guò)少都會(huì)直接造成計(jì)算精度的下降；rmax選取過(guò)大、級(jí)數(shù)展開(kāi)項(xiàng)過(guò)多則會(huì)占用大量的計(jì)算時(shí)間和空間。一般選取弱作用距離在Δrg＜rmax＜2Δrg，Δrg表示規(guī)則網(wǎng)格間距，Taylor展開(kāi)階數(shù)一般取為2，規(guī)則網(wǎng)格點(diǎn)取每波長(zhǎng)10～20個(gè)點(diǎn)。

對(duì)于分析由大量降雨粒子構(gòu)成的雨團(tuán)的電磁散射特性，SMCG方法相比于傳統(tǒng)矩量法將有效地降低計(jì)算復(fù)雜度和內(nèi)存需求。通過(guò)一些典型的算例來(lái)分析降雨粒子群的散射情況。由于星載降雨雷達(dá)工作在K和Ka頻段，雨滴模型不能直接近似成質(zhì)點(diǎn)，還需進(jìn)一步對(duì)雨滴模型進(jìn)行細(xì)分成更小的單元，利用正六面體單元對(duì)雨滴模型進(jìn)行了剖分細(xì)化。

2.實(shí)驗(yàn)結(jié)果分析

Pruppacher等人［11］通過(guò)對(duì)大量實(shí)驗(yàn)數(shù)據(jù)分析指出，降雨粒子半徑在小于140μm時(shí)可以等效為球形，當(dāng)半徑大于140μm雨滴的形狀通常呈橢球形。

降雨粒子的數(shù)學(xué)模型通常等效為橢球，因此，在研究降雨粒子散射特性時(shí)，主要針對(duì)橢球形粒子模型的電磁特性進(jìn)行分析。假設(shè)等效后橢球的旋轉(zhuǎn)軸為z軸，對(duì)應(yīng)軸長(zhǎng)為b，在x，y方向上的橢球軸長(zhǎng)為a，則直角坐標(biāo)系下的橢球方程可寫為

當(dāng)a＞b時(shí)，稱之為扁橢球；a＜b時(shí)，稱之為長(zhǎng)橢球；特殊情況下，a=b時(shí)為球模型。

在實(shí)際降雨過(guò)程中，降雨雷達(dá)觀測(cè)到的都是大量降雨粒子構(gòu)成的雨團(tuán)信息。大量的實(shí)際觀測(cè)數(shù)據(jù)表明，構(gòu)成雨團(tuán)的降雨粒子大都滿足某種粒徑分布（如:L-P，M-P，Gamma分布）的形式。粒徑分布是指在一定觀測(cè)空間中的大量粒子，其粒子半徑的數(shù)學(xué)統(tǒng)計(jì)分布形式。針對(duì)粒徑分布呈隨機(jī)分布的模型進(jìn)行了散射特性的分析，隨機(jī)是指粒子半徑大小在一定范圍內(nèi)是隨機(jī)的。

針對(duì)仿真分析，入射電磁波的工作頻率是基本參數(shù)。例如:現(xiàn)有的TRMM型降雨雷達(dá)入射電磁波的工作頻率選取在Ku波段，這主要是考慮到Ku波段波長(zhǎng)較長(zhǎng)，具有很好的探測(cè)能力。而K和Ka波段的入射波波長(zhǎng)較短，與降雨粒子的等效半徑基本保持在一個(gè)數(shù)量級(jí)。粒子尺寸、相態(tài)等參數(shù)的變化對(duì)其散射特性的影響很明顯。新一代的降雨雷達(dá)（PR2，Precipitation Radar 2）就將 K和Ka工作頻段下降雨粒子的電磁特性作為重要的觀測(cè)對(duì)象。因此，從接近實(shí)際情況的角度考慮，本文選取K和Ka波段作為工作頻率。

降雨粒子模型對(duì)電磁波的折射參數(shù)選取可以參考水和冰的復(fù)折射指數(shù)數(shù)據(jù)［12］。Ka波段下水的復(fù)折射參數(shù)在1.784+0.0013i到4.954+2.79i的范圍內(nèi)。

2.1 固定方位下球形降雨模型的散射特性

圖2所示為直徑1mm的27個(gè)球形降雨粒子，均勻分布在空間中，球心距為1.6mm，復(fù)折射參數(shù)選取為（2.0，0.01），入射波工作頻率為37.5GHz，波長(zhǎng)λ=8mm，沿z方向入射，且只有x方向的電場(chǎng)。SMCG方法將粒子剖分后未知量為5859，強(qiáng)弱相互作用的臨界距離取為rmax=1.4Δrg；泰勒級(jí)數(shù)展開(kāi)的階數(shù)為γ=2.

從圖2可以看出:SMCG方法所得計(jì)算結(jié)果與VIE結(jié)合矩量法計(jì)算所得角分布曲線吻合，從而可以驗(yàn)證本文所采用SMCG方法的正確性。表1給出了SMCG方法和VIE方法在內(nèi)存消耗和計(jì)算時(shí)間上的比較。從表中可以看出:在保證計(jì)算精度的前提下，SMCG方法大大節(jié)省了計(jì)算時(shí)間和內(nèi)存消耗。

圖2 降雨粒子群在固定方位下SMCG法和VIE結(jié)合矩量法計(jì)算所得水平極化散射特性

表1 SMCG和VIE方法耗時(shí)和內(nèi)存消耗對(duì)比

2.2 降雨粒子群在均勻／隨機(jī)分布下的散射特性

對(duì)于微波遙感而言，降雨粒子在空間的分布通常是隨機(jī)的，僅僅研究均勻分布下粒子的散射特性便不足以滿足實(shí)際情況的需求。我們對(duì)比分析了球形降雨粒子在均勻／隨機(jī)分布下的散射特性和不同軸比的橢球近似雨團(tuán)模型均勻分布對(duì)散射特性的影響。降雨粒子團(tuán)在空間的隨機(jī)分布采用蒙特卡羅法來(lái)模擬，且粒子間相互分開(kāi)不接觸。以125個(gè)球形降雨粒子模型隨機(jī)分布在1λ×1λ×1λ空間中為例，圖3描述了該模型的空間分布情況。對(duì)比分析了均勻／隨機(jī)分布和不同軸比的降雨粒子對(duì)散射特性的影響。

圖3 降雨粒子團(tuán)在1λ×1λ×1λ空間中隨機(jī)分布模型

分析球形近似雨團(tuán)模型空間分布對(duì)散射特性的影響。參數(shù)為:直徑為1mm的125個(gè)球形降雨粒子，粒子均勻／隨機(jī)分布在1λ×1λ×1λ空間中，復(fù)折射參數(shù)選取為（2.0，0.01），入射波工作頻率為37.5 GHz，沿x方向極化，波長(zhǎng)λ=8mm，沿z方向入射。SMCG方法將粒子剖分后未知量為27125，每一維方向上規(guī)則網(wǎng)格數(shù)Ng=16（Ng為規(guī)則網(wǎng)格數(shù)），強(qiáng)弱相互作用的臨界距離取為rmax=1.4Δrg；泰勒級(jí)數(shù)展開(kāi)的階數(shù)為γ=2.圖4描述了125個(gè)同樣尺寸的球形近似降雨粒子模型在空間中均勻分布和隨機(jī)分布兩種情況下的電磁特性。從圖中可以看出:當(dāng)散射角較?。ń咏跋蛏⑸洌r(shí)，兩種分布對(duì)散射結(jié)果的影響并不明顯；隨著散射角的增大，兩種分布下的散射強(qiáng)度存在著明顯的差異。換言之，入射波入射方向的不同對(duì)散射體的前向散射影響不大。

圖4 直徑為1mm的125個(gè)球形降雨粒子均勻／隨機(jī)分布于1λ×1λ×1λ空間中的粒子團(tuán)模型在Ka波段下的散射強(qiáng)度

分析不同軸比的橢球形近似雨團(tuán)模型均勻分布對(duì)散射特性的影響。參數(shù)為:等效球直徑r≈1mm，軸比分別為1.0、1.5和2.0，z軸為旋轉(zhuǎn)軸的125個(gè)扁橢球形降雨粒子，粒子均勻分布在1λ×1λ×1λ空間中，復(fù)折射參數(shù)選取為（2.0，0.01），沿x方向極化的入射波波長(zhǎng)λ=8mm，沿z方向入射。粒子剖分密度不變，以上不同軸比粒子剖分后對(duì)應(yīng)的未知量分別為27125，27375，27875，Ng=16，rmax=1.4Δrg和γ=2.

圖5描述了不同軸比的橢球降雨粒子模型的散射特性，模型具有相同的等效球半徑。從圖中可以看出，軸比的不同對(duì)模型的前向散射和后向散射影響均不是很大，散射強(qiáng)度的第一個(gè)谷點(diǎn)隨著降雨粒子軸比的增大而逐漸減小，且谷點(diǎn)的位置也向左偏移，這是由于軸比越大，模型的非球形效應(yīng)越明顯。

圖5 等效球直徑為1mm的125個(gè)軸比分別為1.0、1.5和2.0的橢球形降雨粒子均勻分布于1λ×1λ×1λ空間中的粒子團(tuán)模型在Ka波段下的散射強(qiáng)度

2.3 降雨粒子群在不同頻段下的散射特性

實(shí)際現(xiàn)象表明，降雨粒子在不同頻段下的散射特性存在著明顯的不同。這是因?yàn)楫?dāng)電磁波投射到氣體分子或液態(tài)、固態(tài)的云和降雨粒子上時(shí)，一部分能量被散射，另一部分能量被吸收而轉(zhuǎn)變?yōu)闊崮芑蚱渌问降哪芰?。?duì)隨機(jī)分布在一定空間中的大量降雨粒子在不同頻率點(diǎn)處的散射特性進(jìn)行了研究。為了建立盡可能逼近現(xiàn)實(shí)的降雨模型，獲得所選雨滴模型的理想極化參數(shù)，雨滴的軸比選擇特別重要。Pruppacher等人［11］提出的雨滴軸比與等效直徑的關(guān)系為

Beard等人［13］通過(guò)對(duì)試驗(yàn)數(shù)據(jù)的擬合得到的多項(xiàng)式關(guān)系為

Andsager等人［14］對(duì)諧振雨滴群的形狀進(jìn)行平均，得到的軸比公式為

Keenan等人［15］對(duì)文獻(xiàn)中一些雨滴模型的參數(shù)擬合結(jié)果為

式中:eps=a／b代表軸比；D為橢球等體積球的直徑。a為長(zhǎng)半軸長(zhǎng)，b為短半軸長(zhǎng)，eps＜1表示粒子為長(zhǎng)橢球粒子，eps＞1表示粒子為扁橢球粒子，eps=1表示粒子為球形粒子。

圖6所示為橢球雨滴軸比與其等效半徑R的關(guān)系。所示的四種模型在對(duì)真實(shí)雨滴進(jìn)行近似時(shí)，每一模型只有在某一特定雨滴直徑范圍上好于其它三者。參考圖7所示的不同大小雨滴的橫截面可以看出，當(dāng)雨滴半徑小于0.675mm時(shí)，采用Keenan等人的模型，大于0.675mm小于2.2mm的雨滴采用Andsage等人的模型，其它值時(shí)采用Beard等人的模型。當(dāng)雨滴變大以后，其形狀更加扁平。由軸比圖示可以知道，當(dāng)雨滴等效直徑大于2.2mm時(shí)，Beard等人的模型的軸比值更小，這表明雨滴形狀更扁，更貼近實(shí)際情況。

在500mm×500mm×500mm空間內(nèi)分布了500個(gè)隨機(jī)分布的橢球形降雨粒子，粒子的短半軸長(zhǎng)b小于2mm.考慮軸比eps與等效球半徑R關(guān)系，采用Keenan等人的模型［15］建立降雨粒子模型。

圖8 500mm×500mm×500mm空間內(nèi)隨機(jī)分布500個(gè)半短軸0～2mm橢球形降雨群在不同頻段下的散射特性

對(duì)降雨模型進(jìn)行剖分細(xì)化，未知量為27498.入射波入射方向與z軸夾θ=30°，φ=0°.在10～23 GHz的頻率范圍內(nèi)每0.01GHz取一個(gè)采樣點(diǎn)，雨滴粒子的介電常數(shù)隨著頻率變化，介于4+0.04i和3.3495+0.105i之間，SMCG方法參數(shù)為:Ng=16，rmax=1.4Δrg和γ=2.圖8給出了隨機(jī)分布500個(gè)橢球形降雨群在不同頻段下的散射特性，散射角為0°，是前向散射。從圖中可以看出:在頻率較低時(shí)散射強(qiáng)度曲線變化比較緩慢，隨著工作頻率的提高，散射強(qiáng)度呈現(xiàn)明顯的抖動(dòng)，且抖動(dòng)幅度劇烈。這些都是由于頻率提高，降雨粒子電尺寸與波長(zhǎng)可比擬，對(duì)電磁波的散射更加劇烈導(dǎo)致的。

3.結(jié) 論

從理論研究的角度出發(fā)，針對(duì)新一代星載降雨雷達(dá)新增K和Ka工作頻段測(cè)量的特點(diǎn)，對(duì)均勻／隨機(jī)分布的降雨雨團(tuán)在K和Ka波段下的電磁散射特性進(jìn)行了研究。從中可以看出，散射角較小時(shí)，降雨分布對(duì)散射特性的影響是可以忽略不計(jì)的，而當(dāng)散射角較大時(shí)，這種分布因素的影響就必須加以考慮；隨著雷達(dá)工作波段的提高，粒子的電磁散射特性變化亦更加劇烈，對(duì)較小降雨進(jìn)行測(cè)量具有較高的測(cè)量靈敏度，從而進(jìn)一步提高星載降雨雷達(dá)探測(cè)中、高緯度地區(qū)小雨或小雪的能力。

綜上所述，新一代降雨雷達(dá)新增加K和Ka波段觀測(cè)，可以擴(kuò)展降雨的檢測(cè)范圍，有助于區(qū)分不同類型的降雨，對(duì)提高星載降雨雷達(dá)測(cè)量精度具有重要的意義。

［1］KUMMEROW C，BARNES W，KOZU T.The tropical rainfall measuring mission（TRMM）sensor package［J］.Journal of atmospheric and oceanic technology，1998:809-817.

［2］KOZU T，KAWANISHI T，KUROIWA H，et al.Development of precipitation radar onboard the tropical rainfall measuring mission （TRMM）satellite ［J］.IEEE transactions on geoscience and remote sensing，2001，39（1）:102-116.

［3］唐治華，等.國(guó)外星載氣象雷達(dá)技術(shù)研究［J］.航天器工程.2008，17（5）:88-94.TANG Zhihua，DU Hui.Study on Overseas Spaceborne Meteorological Radar［J］.Spacecraft Engineering，2008，17（5）:88-94.（in Chinese）

［4］BARTON J P.Interal and near-surface electromagnetic fields for an absorbing spheroidal particle with arbitrary illumination［J］.Applied Optics，1995，34（36）:6471-8473.

［5］PURCELL E M，PENNYPACKER C R.Scattering and absorption of light by nonspherical dielectric grains［J］.Astrophys.J.，1973，186:705-714.

［6］WATERMAN P C.Symmetry，unitary，and geometry in electromagnetic scattering［J］.Phys Rev D，1971，3:825-839.

［7］丁大志，劉金權(quán)，胡云琴，等.基于多層快速多極子算法分析平面多層結(jié)構(gòu)［J］.電波科學(xué)學(xué)報(bào)，2009，24（3）:422-426.DING Dazhi，LIU Jinquan，HU Yunqin，et al.Multilayered plane structure using multilevel fast multipole algorithm［J］.Chinese Journal of Radio Science，2009，24（3）:422-426.（in Chinese）

［8］孫玉發(fā)，陳志豪.分析三維隨機(jī)介質(zhì)目標(biāo)散射問(wèn)題的SMCG方法［J］.中國(guó)科學(xué)技術(shù)大學(xué)學(xué)報(bào)，2003，33（3）:345-350.SUN Yufa，CHAN C H.Sparse matrix canonical grid method for scattering from 3Drandom dielectric scatterers［J］.Journal of University of Science and Technology of China，2003，33（3）:345-350.（in Chinese）

［9］CHAN C H.TSANG L..A sparse-matrix／canonical grid method for scattering by many scatterers［J］.Microwave Opt.Technol.Lett.，1995，8（2）:114-118.

［10］丁大志.復(fù)雜電磁問(wèn)題的快速分析和軟件實(shí)現(xiàn)［D］.南京:南京理工大學(xué)，2007.

［11］PRUPPACHER H R，BEARD K V.A wind tunnel investigation of the internal circulation and shape of water drops falling at terminal velocity in air ［J］.The Quarterly Journal of the Royal Meteorological Society，1970，96（408）:247-256.

［12］石廣玉.大氣輻射學(xué)［M］.北京:科學(xué)出版社，2007.

［13］BEARD K V，CHUANG C.A new model for the equilibrium shape of raindrops.Journal of Atmospheric Sciences，1987，44（11）:1509-1524.

［14］ANDSAGER K，BEARD K V，LAIRD N F.Laboratory measurements of axis ratios for large raindrops.Journal of Atmospheric Sciences，1999，56（15）:2673-2683.

［15］KEENAN T D，CAREY L D ，ZRNIC D S，et al.Sensitivity of 5-cm wavelength polarimetric radar variables to raindrop axial ratio and drop size distribution.Journal of Applied Meteorology，2001，40（3）:526-545.