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

丁大志 沈 鹏 陈如山 尚 社 范晓彦

（1.南京理工大学通信工程系，江苏 南京 210094；2.空间微波技术国防科技重点实验室，陕西 西安 710000）

引 言

全球降雨分布及其测量是目前全球环境气候变化研究中的薄弱环节，当前只有少数国家可以利用热带降雨测量（TRMM）卫星获得可靠的全球降雨分布数据［1-2］。TRMM降雨雷达工作在Ku波段，中心频率为13.8GHz，可检测的最小降雨量约为0.6mm／h.为了进一步提高星载降雨雷达对小雨、雪等气候现象的测量灵敏度，新一代星载降雨雷达将增加K和Ka工作频段［3］，使得降雨雷达能够对小雨雪等气候进行更加精确的测量。因此，对K和Ka波段降雨粒子散射特性及不同频率对降雨粒子散射特性的影响进行研究对改进天气和洪水预报、水资源管理以及监测全球变暖与气候变化具有重要意义。

传统分析降雨粒子群散射的方法可以概括为两大类:解析解法和数值近似解法。其中，解析解法主要是分离变量（SVD）法，最早由 Asano等人提出［4］，但是该方法对坐标系的依赖性限制了其对复杂形状粒子散射特性的研究。随着计算机硬件环境的发展，数值类方法得到快速发展。用于解决粒子散射问题的计算方法主要包括离散偶极子近似（DDA）法、T矩阵（T-Matrix）法、基于体积分方程（VIE）的矩量法、稀疏矩阵规则网格（SMCG）法等。DDA方法最早由Purcell和Pennypacker在1973年提出，用于分析星际尘埃物质的散射和吸收问题［5］。但是其对存储量、计算时间的高需求使其发展受到很大限制。T矩阵方法最早由P.C.Waterman［6］于1965年首次引入计算电磁学领域，最显著的特点是它只与粒子自身的属性有关，对于某一固定的散射粒子，计算一次它的T矩阵，便可以得到任意方位下该散射体的散射特征，但是在分析非球对称类模型时，其计算复杂度会大幅度增加。利用体积分方程（VIE）结合矩量法［7］求解宏观粒子的散射问题，可以精确地分析目标的电磁特性，但是其未知量大、计算时间长等缺点使其只能对较小点尺寸的目标进行分析。SMCG方法［8-9］将传统矩量法中的阻抗矩阵分为强相互作用矩阵和弱相互作用矩阵，对弱相互作用矩阵进行Taylor级数展开，利用FFT变换计算弱相互作用矩阵与待求向量的乘积，降低了内存需求和计算复杂度。

简要阐述了电场积分方程和SMCG方法的基本原理。对K和Ka波段下均匀分布的介质降雨粒子的散射特性进行了对比分析，为了更好地模拟实际降雨过程，分析了随机分布的降雨粒子群的散射特性，研究了一定区域粒径随机分布的降雨雨团的散射特性。

1.SMCG方法基本理论

1.1 电场积分方程

对于介质体散射问题，通常可以通过求解表面积分方程和体积分方程来解决。但是，由于降雨粒子为非均匀的介质散射体，一般选择体积分方程求解散射场。主要针对体积分方程进行介绍。

假设散射体的介电常数和电导率分别为ε（r）、σ（r），r为场点，r′为源点入射平面电磁波在散射体上产生的等效感应电流Jeq可表示为［10］

式中:E表示介质散射体上的总电场强度。此等效电流在背景均匀的介质中产生的散射场可以表示为

式中:k20=ω2μ0ε0，且μ0为自由空间中的磁导率，pv表示主值积分。

由于介质体的总场强度可以表示为入射场和散射场总和的形式，即

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

上式便是介质散射体的体积分方程表达式。对于脉冲基点匹配方法，E（r′）·=G（r，r′）可表示为

简化起见，令x1=x x2=y x3=z，则

式（7）中的（r，r′）是对称矩阵，可以将式（6）表示为

假设介质散射体被剖分为n个小体积单元，每个介质体单元内的E（r）和τ（r）都是常量。第m个体积元的体积表示为Vm，rm表示在Vm内选取的任意点。在每个点rm处，式（6）表示成矩阵方程为

则式（10）可以表示为

令［Gxpxq］为N×N阶的矩阵，其中的元素可以定义为

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

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

将上式简化表示为

式中:［G］是3 N×3 N 阶的矩阵，［E］和［Ei］为3 N维的向量。通过求解上式可以确定N个选定点处的电场。

1.2 SMCG方法

SMCG方法的主要思想是根据剖分单元之间相互作用的强弱将阻抗矩阵划分成强相互作用矩阵［GS］和弱相互矩阵［GW］

式中:［GS］为稀疏矩阵；［GW］是稠密矩阵。对于［GS］矩阵可以采用压缩存储的形式节省内存，而弱相互矩阵矢量乘的计算及存储耗费了大部分的CPU时间和内存。为加速弱相互矩阵矢量乘的计算，通常采用SMCG快速算法来加速计算。SMCG方法是将Green函数在规则网格点上进行Taylor级数展开，则弱相互矩阵［GW］可以进一步表示为

式中:K为Taylor级数总的展开项数。则式（16）相应的零阶和高阶解的迭代计算式为

图1 弱相互作用矩阵元素的直接和间接积分路径

图1从二维角度示意了两种积分路径，三维情况可以类似推导出。假定介质单元rm=（xm，ym，zm）与rn=（xn，yn，zn）为弱作用的两个单元，则其阻抗元素G（rm，rn）的积分路径有两种选择:其一是从点rm到点rn的直接路径；其二是从点rm到其对应规则网格点rmo=（xmo，ymo，zmo），再从点rno到与点rn对应规则网格点rno= （xno，yno，zno），最后，从点rno到点rno.间接计算可以通过格林函数关于规则网格点的Taylor级数展开实现，即

式中:xdm= （xm-xmo）；xdn= （xno-xn）；xmn=（xmo-xno），同理，对应的ydm、ydk、ymn、zdm、zdn、zmn有类似的表示。弱相互作用矩阵与待求向量的乘积可以写成

式中:［Tsi］相当于将基元素中心点移到与其最近的规则网格点上；［Gi］为Toeplitz矩阵，该矩阵元素反映了规则网格点间的相互作用；矩阵［Tti］相当于将规则网格点移回到检测中心点。Toeplitz矩阵的行之间满足循环平移关系，进行矩阵矢量乘时可通过快速傅里叶变换（FFT）实现，SMCG算法能够将内存需求从O（N2）降至O（N），计算复杂度从O（N2）降至O（NlogN）.

在对弱相互作用矩阵元素进行Taylor级数展开时，需要对近场作用与远场作用的临界距离rmax和展开级数所取阶数（对应项数）进行折中选取。rmax选取过小、级数展开项过少都会直接造成计算精度的下降；rmax选取过大、级数展开项过多则会占用大量的计算时间和空间。一般选取弱作用距离在Δrg＜rmax＜2Δrg，Δrg表示规则网格间距，Taylor展开阶数一般取为2，规则网格点取每波长10～20个点。

对于分析由大量降雨粒子构成的雨团的电磁散射特性，SMCG方法相比于传统矩量法将有效地降低计算复杂度和内存需求。通过一些典型的算例来分析降雨粒子群的散射情况。由于星载降雨雷达工作在K和Ka频段，雨滴模型不能直接近似成质点，还需进一步对雨滴模型进行细分成更小的单元，利用正六面体单元对雨滴模型进行了剖分细化。

2.实验结果分析

Pruppacher等人［11］通过对大量实验数据分析指出，降雨粒子半径在小于140μm时可以等效为球形，当半径大于140μm雨滴的形状通常呈椭球形。

降雨粒子的数学模型通常等效为椭球，因此，在研究降雨粒子散射特性时，主要针对椭球形粒子模型的电磁特性进行分析。假设等效后椭球的旋转轴为z轴，对应轴长为b，在x，y方向上的椭球轴长为a，则直角坐标系下的椭球方程可写为

当a＞b时，称之为扁椭球；a＜b时，称之为长椭球；特殊情况下，a=b时为球模型。

在实际降雨过程中，降雨雷达观测到的都是大量降雨粒子构成的雨团信息。大量的实际观测数据表明，构成雨团的降雨粒子大都满足某种粒径分布（如:L-P，M-P，Gamma分布）的形式。粒径分布是指在一定观测空间中的大量粒子，其粒子半径的数学统计分布形式。针对粒径分布呈随机分布的模型进行了散射特性的分析，随机是指粒子半径大小在一定范围内是随机的。

针对仿真分析，入射电磁波的工作频率是基本参数。例如:现有的TRMM型降雨雷达入射电磁波的工作频率选取在Ku波段，这主要是考虑到Ku波段波长较长，具有很好的探测能力。而K和Ka波段的入射波波长较短，与降雨粒子的等效半径基本保持在一个数量级。粒子尺寸、相态等参数的变化对其散射特性的影响很明显。新一代的降雨雷达（PR2，Precipitation Radar 2）就将 K和Ka工作频段下降雨粒子的电磁特性作为重要的观测对象。因此，从接近实际情况的角度考虑，本文选取K和Ka波段作为工作频率。

降雨粒子模型对电磁波的折射参数选取可以参考水和冰的复折射指数数据［12］。Ka波段下水的复折射参数在1.784+0.0013i到4.954+2.79i的范围内。

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

图2所示为直径1mm的27个球形降雨粒子，均匀分布在空间中，球心距为1.6mm，复折射参数选取为（2.0，0.01），入射波工作频率为37.5GHz，波长λ=8mm，沿z方向入射，且只有x方向的电场。SMCG方法将粒子剖分后未知量为5859，强弱相互作用的临界距离取为rmax=1.4Δrg；泰勒级数展开的阶数为γ=2.

从图2可以看出:SMCG方法所得计算结果与VIE结合矩量法计算所得角分布曲线吻合，从而可以验证本文所采用SMCG方法的正确性。表1给出了SMCG方法和VIE方法在内存消耗和计算时间上的比较。从表中可以看出:在保证计算精度的前提下，SMCG方法大大节省了计算时间和内存消耗。

图2 降雨粒子群在固定方位下SMCG法和VIE结合矩量法计算所得水平极化散射特性

表1 SMCG和VIE方法耗时和内存消耗对比

2.2 降雨粒子群在均匀／随机分布下的散射特性

对于微波遥感而言，降雨粒子在空间的分布通常是随机的，仅仅研究均匀分布下粒子的散射特性便不足以满足实际情况的需求。我们对比分析了球形降雨粒子在均匀／随机分布下的散射特性和不同轴比的椭球近似雨团模型均匀分布对散射特性的影响。降雨粒子团在空间的随机分布采用蒙特卡罗法来模拟，且粒子间相互分开不接触。以125个球形降雨粒子模型随机分布在1λ×1λ×1λ空间中为例，图3描述了该模型的空间分布情况。对比分析了均匀／随机分布和不同轴比的降雨粒子对散射特性的影响。

图3 降雨粒子团在1λ×1λ×1λ空间中随机分布模型

分析球形近似雨团模型空间分布对散射特性的影响。参数为:直径为1mm的125个球形降雨粒子，粒子均匀／随机分布在1λ×1λ×1λ空间中，复折射参数选取为（2.0，0.01），入射波工作频率为37.5 GHz，沿x方向极化，波长λ=8mm，沿z方向入射。SMCG方法将粒子剖分后未知量为27125，每一维方向上规则网格数Ng=16（Ng为规则网格数），强弱相互作用的临界距离取为rmax=1.4Δrg；泰勒级数展开的阶数为γ=2.图4描述了125个同样尺寸的球形近似降雨粒子模型在空间中均匀分布和随机分布两种情况下的电磁特性。从图中可以看出:当散射角较小（接近前向散射）时，两种分布对散射结果的影响并不明显；随着散射角的增大，两种分布下的散射强度存在着明显的差异。换言之，入射波入射方向的不同对散射体的前向散射影响不大。

图4 直径为1mm的125个球形降雨粒子均匀／随机分布于1λ×1λ×1λ空间中的粒子团模型在Ka波段下的散射强度

分析不同轴比的椭球形近似雨团模型均匀分布对散射特性的影响。参数为:等效球直径r≈1mm，轴比分别为1.0、1.5和2.0，z轴为旋转轴的125个扁椭球形降雨粒子，粒子均匀分布在1λ×1λ×1λ空间中，复折射参数选取为（2.0，0.01），沿x方向极化的入射波波长λ=8mm，沿z方向入射。粒子剖分密度不变，以上不同轴比粒子剖分后对应的未知量分别为27125，27375，27875，Ng=16，rmax=1.4Δrg和γ=2.

图5描述了不同轴比的椭球降雨粒子模型的散射特性，模型具有相同的等效球半径。从图中可以看出，轴比的不同对模型的前向散射和后向散射影响均不是很大，散射强度的第一个谷点随着降雨粒子轴比的增大而逐渐减小，且谷点的位置也向左偏移，这是由于轴比越大，模型的非球形效应越明显。

图5 等效球直径为1mm的125个轴比分别为1.0、1.5和2.0的椭球形降雨粒子均匀分布于1λ×1λ×1λ空间中的粒子团模型在Ka波段下的散射强度

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

实际现象表明，降雨粒子在不同频段下的散射特性存在着明显的不同。这是因为当电磁波投射到气体分子或液态、固态的云和降雨粒子上时，一部分能量被散射，另一部分能量被吸收而转变为热能或其它形式的能量。对随机分布在一定空间中的大量降雨粒子在不同频率点处的散射特性进行了研究。为了建立尽可能逼近现实的降雨模型，获得所选雨滴模型的理想极化参数，雨滴的轴比选择特别重要。Pruppacher等人［11］提出的雨滴轴比与等效直径的关系为

Beard等人［13］通过对试验数据的拟合得到的多项式关系为

Andsager等人［14］对谐振雨滴群的形状进行平均，得到的轴比公式为

Keenan等人［15］对文献中一些雨滴模型的参数拟合结果为

式中:eps=a／b代表轴比；D为椭球等体积球的直径。a为长半轴长，b为短半轴长，eps＜1表示粒子为长椭球粒子，eps＞1表示粒子为扁椭球粒子，eps=1表示粒子为球形粒子。

图6所示为椭球雨滴轴比与其等效半径R的关系。所示的四种模型在对真实雨滴进行近似时，每一模型只有在某一特定雨滴直径范围上好于其它三者。参考图7所示的不同大小雨滴的横截面可以看出，当雨滴半径小于0.675mm时，采用Keenan等人的模型，大于0.675mm小于2.2mm的雨滴采用Andsage等人的模型，其它值时采用Beard等人的模型。当雨滴变大以后，其形状更加扁平。由轴比图示可以知道，当雨滴等效直径大于2.2mm时，Beard等人的模型的轴比值更小，这表明雨滴形状更扁，更贴近实际情况。

在500mm×500mm×500mm空间内分布了500个随机分布的椭球形降雨粒子，粒子的短半轴长b小于2mm.考虑轴比eps与等效球半径R关系，采用Keenan等人的模型［15］建立降雨粒子模型。

图8 500mm×500mm×500mm空间内随机分布500个半短轴0～2mm椭球形降雨群在不同频段下的散射特性

对降雨模型进行剖分细化，未知量为27498.入射波入射方向与z轴夹θ=30°，φ=0°.在10～23 GHz的频率范围内每0.01GHz取一个采样点，雨滴粒子的介电常数随着频率变化，介于4+0.04i和3.3495+0.105i之间，SMCG方法参数为:Ng=16，rmax=1.4Δrg和γ=2.图8给出了随机分布500个椭球形降雨群在不同频段下的散射特性，散射角为0°，是前向散射。从图中可以看出:在频率较低时散射强度曲线变化比较缓慢，随着工作频率的提高，散射强度呈现明显的抖动，且抖动幅度剧烈。这些都是由于频率提高，降雨粒子电尺寸与波长可比拟，对电磁波的散射更加剧烈导致的。

3.结 论

从理论研究的角度出发，针对新一代星载降雨雷达新增K和Ka工作频段测量的特点，对均匀／随机分布的降雨雨团在K和Ka波段下的电磁散射特性进行了研究。从中可以看出，散射角较小时，降雨分布对散射特性的影响是可以忽略不计的，而当散射角较大时，这种分布因素的影响就必须加以考虑；随着雷达工作波段的提高，粒子的电磁散射特性变化亦更加剧烈，对较小降雨进行测量具有较高的测量灵敏度，从而进一步提高星载降雨雷达探测中、高纬度地区小雨或小雪的能力。

综上所述，新一代降雨雷达新增加K和Ka波段观测，可以扩展降雨的检测范围，有助于区分不同类型的降雨，对提高星载降雨雷达测量精度具有重要的意义。

［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］唐治华，等.国外星载气象雷达技术研究［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］丁大志，刘金权，胡云琴，等.基于多层快速多极子算法分析平面多层结构［J］.电波科学学报，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］孙玉发，陈志豪.分析三维随机介质目标散射问题的SMCG方法［J］.中国科学技术大学学报，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］丁大志.复杂电磁问题的快速分析和软件实现［D］.南京:南京理工大学，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］石广玉.大气辐射学［M］.北京:科学出版社，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.