双频毫米波雷达反演云微物理参数的模拟研究

杨敏 黄兴友 赵泉钦 张正 杨恒祥 周红根

(1 江苏省气象探测中心，南京 210019；2 南京信息工程大学 气象灾害预报预警与评估协同创新中心，南京 210044)

引 言

云是地球—大气系统的重要组成部分，它不仅是产生降水的前期物，也是影响地表辐射能量收支和全球气候变化的关键因素。毫米波雷达通过接收云粒子群对雷达电磁波的后向散射信号实现对云的遥感探测，并采用多频雷达组合探测，成为获取云滴谱的新技术。如果选择适当组合雷达的频率或波长，使得云滴粒子对其中一个频段的雷达电磁波产生瑞利散射，对另一个产生非瑞利散射(米散射)，那么双频雷达探测就可较准确地反演出云滴谱参数并得到其他微物理特性[1-2]。利用云滴群对双频雷达的衰减差异、垂直速度差异等，也能得到云滴群的粒子大小和数量的总体特性，有助于反演云滴谱和微物理参数。

在地基雷达的双频反演方案早期研究中，任何距离库上的滴谱参数都是计算两个波段的差分衰减得到的[3-4]。Sekelsky，et al[5]提出一种利用33 GHz和95 GHz的反射率因子差和平均多普勒速度差分来估算云粒子中体积直径的模型。Rose，et al[6]介绍了星载双频雷达回波反演的系统方法，通过滴谱分布参数修正控制后向散射和前向散射的电磁波的必要积分方程，并使用这些方程重新评估向后迭代算法。LIAO，et al[7]利用了X、W波段雷达反射率因子和多普勒速度的测量结果来估算雪与雨的微物理特性。Tridon，et al[8]基于双波长多普勒频谱理论，研究了Ka-W双波段反演大气降水滴谱分布特征和空气状态参数。吴琼等[9]使用星载雷达模拟器输出的模拟数据，假设冰晶为球形，利用双频比对星载双频云雷达最佳频点组合以及联合反演冰云微物理参数的算法进行相关研究。

研究表明，双频测量能够更好地估计滴谱分布参数，提高云参数反演的精度[10]。国内目前监测云要素的地基云雷达多为35 GHz。2014年，安徽38所和中国航天科技集团九院704所成功研制了94 GHz毫米波云雷达。此前也开展了35 GHz和94 GHz云雷达联合观测试验，且我国也计划发展星载94-220 GHz的双频测云雷达。因此针对这3个频段，对35-94 GHz、35-220 GHz和94-220 GHz三种组合的双频毫米波雷达反演云微物理参数进行模拟研究。此外，探讨双频算法对滴谱形状、粒子形状的敏感性，比较两种方案反演水云和冰云的优劣势，并结合模拟数据试验了三种频段组合的反演结果。

1 双频反演算法的理论基础

两部云雷达的频率分别为Flow和Fhigh(low和high分别表示双频组合中的低频和高频)，等效雷达反射率因子Ze与云滴谱n(D)、后向散射截面α(D,T,λ)以及雷达波长λ有关，具体表达式如下：

(1)

其中：Ki为介电因子，等于(m2-1)/(m2+2)；m为粒子的复折射指数，与频率和温度有关；D为云滴直径。

(2)

(3)

其中：Vf(D)为直径为D的云滴粒子在静止空气中的下落末速度。

在瑞利散射条件下，粒子的雷达截面σ为：

(4)

云粒子的下落末速度与云滴直径的关系可近似表示为[11]：

Vf(D)=avDbv，

(5)

其中：对于球形液态云滴，D<90 μm时，av=3.26×107m-1·s-1，bv=2[12]。随着粒径的增长(90

对于粒径更大的雨滴(D≥600 μm)，Atlas， et al[13]给出经验公式：

Vf(D)=9.65-10.30 exp(-6D) ,

(6)

对于非球形冰晶粒子，其粒子落速则参考Mitchell， et al[14]和Heymsfield, et al[15]的结果：av=211，bv=0.416。

在云滴谱符合伽马分布的假设下[16-17]，云滴谱为：

(7)

其中：N0为滴谱截距参数(单位：个·cm-3·μm-μ+1)；μ为滴谱形状参数，为无量纲参数；D0为中值直径(单位：μm)。

描述两个频段雷达反射率因子差异的双频比定义为：

(8)

其中：Zlow(R)和Zhigh(R)分别为距离雷达R处的粒子群在频率为Flow和Fhigh时的雷达反射率因子。这样，在μ一定时，DFR的大小与N0无关，只与雷达频率F、复折射指数m和粒子直径D有关。不考虑衰减以及信号的测量误差，那么DFR大小取决于非瑞利后向散射效应差异。

图1 球形水滴在35、94以及220 GHz的后向散射截面(a)和衰减截面(b)Fig.1 Back-scattering cross section (a) and attenuation cross section(b) of spherical water droplets at 35， 94 and 220 GHz

类似地，多普勒速度差为两个频段雷达测的多普勒速度差，定义为：

(9)

其中：V(Flow)和V(Fhigh)分别为距离雷达R处的粒子在频率为Flow和Fhigh时所测的多普勒速度。结合公式(2)、(3)，多普勒速度差DDV的大小与N0无关且不受空气垂直运动速度ω的影响。

根据散射理论和云滴粒子的复折射指数数据，可以计算出不同大小的云滴粒子群的双频比DFR和多普勒速度差DDV，这样DFR和DDV分别和粒子中值直径D0成一一对应的关系，即可建立DFR-D0和DDV-D0数据库。在对雷达系统进行精确标定后，通过双频雷达系统获取的雷达反射率因子和多普勒速度，可以得到DFR和DDV，再根据DFR和DDV与D0的关系，就能反演得到云滴粒子的D0。

2 云参数的双频反演

一般对于液相云滴做(椭)球形处理，而冰晶在自然界中呈现出各种各样的形状：单晶、聚合物以及霰。本文采用的冰云模型为：D<60 μm时，100%滴晶；602 000 μm，97%子弹花，3%聚合物[18]。

根据Ray[19]以及Mätzler[20]的理论，三种频段下水和冰的复折射指数的计算结果如表1所示，其中n是普通的折射指数，k是吸收系数。

在计算雷达后向散射截面时，利用DDA算法[21-23]计算液相云滴以及冰晶云滴的雷达后向散射截面，该模型允许在计算粒子散射特性时改变形状、尺寸和方向。图1为利用DDA计算得到的当温度为0 ℃时球形液相云滴在3个频段下的后向散射截面和衰减截面。图2则是非球形冰晶(以六棱柱为例)当温度为-20 ℃时在3个频段下的后向散射截面和衰减截面。

表1 不同温度与频率时水和冰的n、k值Table 1 The n and k values of water and ice at different temperatures and frequencies

图2 六棱柱冰晶在35、94以及220 GHz的后向散射截面(a)和衰减截面(b)Fig.2 Back-scattering cross section (a) and attenuation cross section(b) of hexagonal prism ice crystal at 35， 94 and 220 GHz

图3 球形液相云(雨)滴双频比DFR(a)和多普勒速度差DDV(b)与中值直径D0的关系Fig.3 The relationship between the dual frequency ratio DFR (a) and the Doppler velocity difference DDV(b) ofspherical liquid cloud (rain) drop and the median diameter D0

2.1 水云的DFR-D0和DDV-D0关系

利用公式(8)、(9)，得到35-94 GHz、35-220 GHz和94-220 GHz三种频段组合下液相云滴DFR和D0以及DDV和D0的关系曲线(图3)。在计算中，采取控制变量法，D0为唯一自变量，温度、滴谱形状因子和液滴形状为固定值，分别为0 ℃、0和球形。

图3表明，对于液相云(雨)滴，当D0从0～2 mm变化时，雷达频点差异越大，同一D0所对应的双频比DFR也较大，35-220 GHz的DFR动态范围是最大的，35-94 GHz的DFR 动态范围最小。D0<0.3 mm时，94-220 GHz的DFR最大，35-94 GHz的DFR最小且为负值；D0>0.3 mm时，35-220 GHz组合的DFR反超94-220 GHz，35-94 GHz的DFR最小且由负值增加为正值。35-220 GHz 的DDV动态范围最大，94-220 GHz 的DDV的动态范围最小。三种频段组合的DFR随着D0的增加而增加，DDV则先增后减。在液相小云滴范围内(液相云的最大液滴直径通常小于0.5 mm)，35-94 GHz的DFR对D0敏感性很低，不适合反演云滴谱。35-220 GHz相较于94-220 GHz，DFR与DDV的动态范围都较大，且35 GHz一般对液相云滴以及弱降水的探测较为有效，而94 GHz则对小粒子和冰晶更为敏感。虽然在液相小云粒范围内，35-94 GHz的DDV对D0敏感性很低，但是考虑到高频的强衰减，该组合最适合反演水云以及弱降水粒子。

2.1.1 水云的DFR-D0和DDV-D0关系对μ的敏感性

当中值直径D0一定时，对于不同的粒子谱形状因子μ，云滴粒子大小占比不同，所以粒子谱形状因子μ影响了DFR-D0和DDV-D0关系。为了说明DFR和DDV的结果对μ的敏感性，图4提供了DFR-D0和DDV-D0随μ变化的一些计算结果。在计算中，μ从-1变到2，粒子形状设置为球形，温度为0 ℃。结果表明，当D0从0～2 mm变化时，在D0较小时，形状因子μ对DFR-D0以及DDV-D0关系基本上没有影响，但是随着D0的增大，μ对DFR和DDV的影响也逐渐明显，μ对DDV-D0的影响大于对DFR-D0的影响。由图4a、c、e，当D0较小时，μ对DFR几乎无影响，随着D0增大，DFR与μ先成反比后成正比；类似的，由图4b、d、f，当D0较小时，μ对DDV几乎无影响，随着D0增大，DDV与μ成反比。当D0从0～2 mm变化时，μ导致的35-94 GHz、35-220 GHz和94-220 GHz三组频段下的DFR最大差异分别达到了3.2 dB、3.5 dB 和2.4 dB，μ导致的35-94 GHz、35-220 GHz和94-220 GHz三组频段下的DDV最大差异则分别为1.7 、2.3 和1.3 m·s-1。形状因子μ越小，DFR和DDV对D0的微小变化就越敏感。针对液相云滴，μ对DFR-D0以及DDV-D0的影响较小。

2.1.2 水云的DFR-D0和DDV-D0关系对γ的敏感性

轴比γ造成的DFR-D0和DDV-D0关系差异主要是由于当粒子体积保持恒定时，对于不同轴比的云滴粒子，其后向散射截面的差异。为了说明液相云滴粒子的形状的改变对DFR-D0和DDV-D0关系的影响，模拟扁球体的粒子，考虑雷达垂直发射电磁波，线型偏振的特性，假定扁球形粒子的对称轴沿垂直方向，它们的长轴在水平面上，长轴与云雷达的偏振方向一致。轴比γ从1(球形)到0.5(扁椭球形)变化，固定形状因子μ=0，图5为计算结果，结果表明，对于三组频段组合，当μ固定时，轴比γ对DFR-D0的影响小于对DDV-D0的影响，球形液滴产生最大动态范围的DFR和DDV。由图5a、c、e，当D0从0～2 mm变化，D0较小时，DFR与轴比成反比，D0较大时，DFR与轴比成正比；轴比导致的35-94 GHz、35-220 GHz和94-220 GHz三组频段下的DFR最大差异分别达到了1.8 dB、3.6 dB 和2.0 dB。所以当使用双频雷达的反射率和多普勒速度探测雨的参数时，γ对DFR-D0以及DDV-D0关系的影响不可忽略，而对于液相云滴，γ对DFR-D0以及DDV-D0几乎无影响。由图5b、d、f，当D0从0～2 mm变化时，DDV整体上呈现先增后减的趋势，球形液滴产生最大动态范围的DDV，而最偏心的液滴(γ=0.5)产生最小动态范围的DDV；轴比导致的35-94 GHz、35-220 GHz和94-220 GHz三组频段下的DDV最大差异分别达到了0.9、0.7和-0.7 m·s-1。因为低频段的小粒子的散射主要是瑞利散射，所以轴比对小粒子的后向散射截面几乎没有影响。对于随机取向的非球形颗粒，DFR和DDV与D0关系的结果类似于球体[7]。

图4 不同粒子谱形状因子下的球形液相云滴DFR(a、c、e)和DDV(b、e、f)与D0的关系Fig.4 The relationship between DFR (a， c， e) and DDV (b， e， f) of spherical liquid cloud dropletsunder different particle spectrum shape factors and D0

图5 不同粒子形状下的球形液相云滴DFR(a、c、e)和DDV(b、d、f)和D0的关系Fig.5 The relationship between DFR (a， c， e) and DDV (b， d， f) of spherical liquid cloud droplets under different particle shapes and D0

图6 冰晶双频比DFR(a)和多普勒速度差DDV(b)与中值直径D0的关系Fig.6 The relationship between dual frequency ratio DFR (a) and Doppler velocity difference DDV(b) of ice crystal and median diameter D0

2.2 冰云的DFR-D0和DDV-D0关系

类似的，计算得到35-94 GHz、35-220 GHz和94-220 GHz三种频段组合下冰晶DFR和DDV与D0的关系曲线(图6)。在计算中，温度和滴谱形状因子设置为固定值，分别为-20 ℃、0，冰晶形状基于Baum， et al[18]的结果。由图6a，当D0在0～2 mm范围内时，35-220 GHz的DFR动态范围最大，约为16.5 dB，94-220 GHz的DFR 动态范围最小，约7 dB，35-94 GHz的DFR对小粒子最不敏感，且35-94 GHz和94-220 GHz都出现了双解。由图6b，多普勒速度差DDV随着D0的增大波动较大，且动态范围都很小。原则上，DFR和DDV的测量可以用于独立地估计D0，相关参数也可以调整使其一致性最大化。然而实际DDV测量时不可避免会出现种种噪声，这种噪声相对小动态范围的DDV将会导致D0估计值与实际值相比会出现非常大的偏差，即从DDV推导出来的D0是没有指导意义的。

2.2.1 冰云的DFR-D0和DDV-D0关系对μ的敏感性

图7为冰云DFR-D0和DDV-D0随粒子谱形状因子μ变化的一些计算结果。在计算中，μ从-1变到2，温度为-20 ℃。结果表明，当D0从0～2 mm变化时，μ对DDV-D0的影响大于对DFR-D0的影响。冰晶对毫米波的衰减非常有限，在没有云水的情况下回波可以轻松订正。因此，对于冰云，更适合利用DFR反演云参数。

图9 冰云反演廓线和模拟廓线的对比：(a)衰减后的数据反演结果；(b)衰减订正后的反演结果Fig.9 Comparison of ice cloud retrieval profile and simulated profile：(a) retrieval profile after attenuation；(b) retrieval profile after attenuation correction

3 双频反演模拟评估检验

为了检验算法的精度，数据仿真两组云滴谱数据，水云模拟的观测方式为地基双频雷达，冰云模拟的观测方式为星载双频雷达。毫米波雷达在衰减订正时不仅要考虑云雨滴粒子对电磁波的衰减，还要考虑大气气体的吸收衰减[24]。计算中，假定μ=0，即云滴谱函数为指数分布。模拟的云层距离分辨率为100 m，云厚1.2 km，每个距离库的D0固定为0.4 mm，其中模拟的液相云云层高度为1.0～2.2 km，冰云云层高度则为4.0～5.2 km(N0固定为8×103mm-1·m-3)[25]，模拟冰云的高度较低是为了最大化模拟大气对电磁波的衰减(标准大气中，0 ℃层约2.3 km)。不同通道组合对不同粒径的云滴反演效果是不同的，接下来就三组频段组合对模拟数据的反演效果分别进行评估检验。

3.1 DDV反演水云效果的检验

实际两个频段测量的雷达反射率的差异，表示为DFRm，是两个频段的衰减差异、非瑞利后向散射效应差异以及与测量值相关的测量误差的综合结果。干的冰晶对雷达信号的衰减非常小，而云水对雷达信号的衰减则是不可忽视的，而且对流层水汽吸收也很重要。实现两个频段的雷达反射率因子的完全订正是准确提取DFR估算D0的关键，就目前来说这项工作很难做到。与基于DFR的方法相比，DDV技术与路径衰减无关，如果不考虑测量误差，那么用DDV估算液相云滴D0将是直接有效的。

实际探测中雷达系统噪声、定标误差以及随机误差不可避免，所以，为了更好地评估DDV反演水云的效果，取两部雷达的DDV相对测量误差分别在±10%、±15%、±20%、±25%范围内。

图8为反演后得到的D0廓线，对其相应的均值以及标准偏差的统计结果如表2所示。从图8和表2可以看出，当D0=0.4 mm，相对测量误差一定时，35-94 GHz反演精度最高，效果更稳定，35-220 GHz次之，94-220 GHz反演结果波动最大。不同频段组合在相对测量误差一定时，反演效果不同。结合图3b，随着粒径逐渐增大到0.5 mm，94-220 GHz的DDV对D0的敏感性开始降低，这时微小的测量误差也会带来较大的反演误差，这在模拟数据中表现得也很明显。所以同一频段组合在相对测量误差一定时，反演效果也是不同的，要想提高反演精度，不仅要尽量降低各种噪声，更要针对反演目标特性，选择对目标粒径敏感性较强的频段组合。

表2 DDV相对测量误差对反演结果的影响Table 2 The influence of DDV relative measurement error on retrieval results

3.2 DFR反演冰云效果的检验

应用DFR反演冰云D0时，反演的依据是等效雷达反射率因子Ee，而双频反演的实际输入值则是实测雷达反射率因子Zm，两者关系为[26]：

(10)

(11)

其中：R为雷达与云层的距离；kc为云滴衰减系数；ka为大气吸收系数。

模拟仿真时，不考虑雷达系统噪声、定标误差以及随机误差等，通过DDA[21-23]计算云滴衰减截面Qt可由公式(11)得到云滴衰减系数kc，利用Liebe模式[27]计算大气对电磁波的吸收系数ka，其中大气廓线数据取自美国(1976)标准大气。由公式(10)就模拟得到每个距离库经过云滴和大气衰减后的雷达反射率因子Zm。

类似地，为了更好地评估DFR反演冰云的效果，不考虑衰减取两部雷达的DFR相对测量误差分别在±10%、±15%、±20%、±25%范围内。

图10为反演后得到的D0廓线，对其相应的均值以及标准偏差的统计结果如表3所示。

图10 不同相对测量误差下的DFR反演结果：(a)相对测量误差(-10%，10%)；(b)相对测量误差(-15%，15%)；(c)相对测量误差(-20%，20%)；(d)相对测量误差(-25%，25%)Fig.10 DFR retrieval results under different relative measurement errors：(a) relative measurement error (-10%， 10%)；(b) relative measurement error (-15%， 15%)；(c) relative measurement error (-20%， 20%)；(d) relative measurement error (-25%， 25%)

表3 DFR相对测量误差对反演效果的影响Table 3 The influence of DFR relative measurement error on retrieval results

从图10和表3可以看出，当D0=0.4 mm，DFR相对测量误差一定时，35-94 GHz反演效果更稳定，35-220 GHz次之，94-220 GHz反演结果波动最大。在DFR相对测量误差一定时，不同频段组合对同一粒径的反演效果是不同的。

4 结论

选取3个频段组合研究了双频反演水云和冰云滴谱参数的算法，并利用模拟数据分别进行了反演效果检验。主要结论如下：

(1)双频雷达提供了一种推断云滴谱信息的方法，虽然DFR和DDV方法都受到假定粒子滴谱分布模型可能不准确性的影响，但是水云中基于DFR反演的主要误差来源是雷达信号的衰减，而径向速度可以近似认为不受衰减影响。且水云一般高度较低，当星载雷达对其测量时，不匹配的波束宽度等测量误差将会对测量结果产生很大的系统误差，而双频精确反演的基础之一是两个频段测量目标的一致性。如果使用地基雷达的测量数据，在一定程度上可以减小这些测量不确定性。考虑到高频段雷达回波信号对水云以及对流层丰富的水汽衰减订正相关的不确定性，理论上地基35-94 GHz组合的DDV技术更适合反演水云或者弱降水粒子。

(2)在冰云中，DDV的动态范围很小，系统和测量误差相对于信号的动态范围是大的，因此应用DDV反演冰云有限。对于干的冰云粒子，采取高空观测模式，那么其衰减几乎可以忽略(高空大气稀薄，衰减也非常有限)，基于DFR的方法比DDV更准确。由于雷达频率越高，雷达对小粒子的探测更为敏感，星载94-220 GHz组合的DFR技术更适合高空卷云的探测。

(3)在双频反演中，随着粒径的增大，DFR以及DDV出现了双解的问题。可以考虑引用雷达反射率因子和多普勒速度作为第二约束条件，可以在很大程度上解决这种不确定性。