降水是水循环的重要环节， 雷达是准确迅速地探测降水的一种重要而有效方式。目前广泛使用的多普勒天气雷达估测降水， 其主要原理是利用Z-I关系: Z=AIb， 其中A和b均由统计关系给出。由于雨滴谱分布随季节、 地域、 降水类型等变化， 使这种方式的估测结果存在很大的不确定性(张培昌等， 2000; 吴亚昊等， 2016; 牛生杰等， 2002)。双频雷达通过联合雷达的两个频率通道的返回信号， 能够更准确地提供云、 降水的微物理参数及降水率的估测（Liao et al， 2005; Mardiana et al， 2004）。

国外利用双频雷达反演云、 降水的微物理参数分布及降水率的研究已经有40多年的历史。1991年6月， 在CaPE（Convection and Precipitation Electrification Experiment）试验中， 美国国家航空航天局的喷气式飞机T-39上搭载了日本通信实验室发展的工作波长为X波段和Ka波段的双频雷达， 对层云降水进行了有组织的飞行观测试验， 采用后向反演算法（BRM）反演出层云降水的雨滴谱及垂直廓线(Meneghini et al， 1992)。2014年全球降水观测卫星（GPM）发射成功， 作为TRMM的后续卫星降水计划， 其核心观测平台（GPMCO）搭载了全球首个星载双频降水雷达（DPR）， 工作频率分别为Ku波段和Ka波段（分别为13 GHz和35 GHz）（唐国强等， 2015）。针对GPM上的DPR， Adhikari et al（2007）发展了基于GPM的双频降水雷达的差分衰减（DA）修订方法来反演降水粒子的滴谱分布及降水率。目前， 国内已经开始发展双频雷达， 航天部北京无线电测量研究所和青海大学已经联合研制了Ku和Ka波段的双频雷达， 并于2016年底进行了测试验收。中国科学院大气物理研究所正在研制三波长雷达， 中国气象科学研究院灾害天气国家重点实验室对已有的Ka波段云雷达升级改造为双频雷达（Ka/Ku）， 并已测试验收， 于2019年4月初开始在广东开展外场观测试验。

双频测雨雷达通过估计雨滴谱参数来计算降水率， 雨滴谱分布通常采用指数分布、 Gamma分布、 或者对数分布等函数来表示。国内外许多学者很早对实际雨滴谱分布进行过研究， 雨滴在下落过程中不断发生碰并和破碎最终达到一个稳定的状态， 得到一个稳定的雨滴谱分布(Srivastava， 2010)。Marshall et al（1948)最早提出指数形式的M~P分布。Ulbrich（1983）则认为Gamma分布使用三种参数， 能更好地描述实际滴谱分布。M~P分布只适用于稳定的层云降水， 而于对流云等波动较大的降水， 滴谱两端偏差较大（陈宝君等， 1998; Carbone et al， 1977; Joss et al， 1978）。相比较而言Gamma分布不仅能反映层状云雨滴谱分布， 对积雨云和混合云雨滴谱也有较为理想的描述， 故用Gamma分布可以提高整个滴谱范围的拟合精度。但是由于真实的雨滴谱分布本身随时空变化， 双频技术反演的雨滴谱参数并不能完全真实描述实际的雨滴谱分布， 而且反演方法本身存在一定误差， 用双频技术反演得到的降水也存在一定的不确定性（Liao et al， 2014）。

Willis提出的Gamma分布有三个未知参量: 标准化截距参数N_w(与粒子数浓度有关)、 形状因子μ、 质量加权的粒子直径D_m。用双频雷达技术反演雨滴谱分布参数及降水率时需要先假定Gamma分布的形状因子μ， 关于μ的取值， 国外有很多学者曾做过相关研究。Le et al(2014)提出， 利用双频理论反演雨滴谱参数时， 在液态降水中Gamma模型的形状因子μ取为0最为合适， 在混合态和固态降水中μ一般取值为3。Liao et al（2014）利用美国爱荷华州洪水试验中的雨滴谱数据分析了采用不同μ值， 对利用DPR估测降水产生的不确定性， 认为当μ取3~6时反演效果最好。或者当μ和斜度因子有某种关系时， 反演的相对误差也很小， 都维持在10%以内。然而张鸿发等(1989）曾指出由于我们国家所处的地理位置、 降水机制和粒子下落过程与国外不同， 雨滴谱与国外也有很大差别。因而利用国内的雨滴谱仪资料开展雨滴谱参数对双频技术反演效果的影响研究， 对促进双频雷达技术在我国的应用和发展具有重要意义。

本研究在利用双频雷达技术反演降水的微物理参数时， 采用Gamma分布来描述雨滴谱分布， 并对其在双频雷达估测降水中的不确定性进行了评估。研究中利用雨滴谱仪测量的DSD数据来计算Ka和Ku波段的雷达反射率因子。在假定的DSD模型下， 这两个波段的雷达反射率因子被用来估测D_m、 衰减系数k及降水率R。双频雷达通常采用两个波长雷达反射率因子的差值（差分频率比DFR）及较长波长的雷达反射率因子来反演D_m。将雨滴谱仪计算得到的D_m、 衰减系数和降水率作为真值， 和双频雷达技术反演得到的估计值进行比较， 考察采用不同的μ值时， 双频技术反演的D_m、 衰减系数、 降水率的估计值偏离真值的程度， 从而衡量双频反演技术中DSD参数造成的影响。研究数据来自广州龙门和西藏那曲， 这两个地区一个位于华南区域， 一个属于青藏高原， 气候地形都有很大差异， 对于探究雨滴谱参数对反演效果影响的地区差异性具有很好的代表性。

本研究以架设在广东龙门和西藏那曲气象局观测场的两个雨滴谱仪的连续2~4年的夏季雨滴谱资料为基础。观测站的雨滴谱仪采用德国OTT公司设计生产的PARSIVEL(Particle Size and Velocity)激光雨滴谱仪。

PARSIVEL激光雨滴谱仪的基本原理是通过发射平行激光束同时测定粒子直径和下落速度。激光测量的优越性表现在可以测量微小物体。PARSIVEL雨滴谱仪可以测量最小直径0.2 mm的降水粒子。其中降水粒子直径大小的测量是由降水粒子通过水平光束时， 相应直径遮住部分光束产生的输出电压的大小确定。而用电子信号的持续时间即可推算出降水粒子的下降速度， 其中持续时间指的是降水粒子从开始进入光束到完全离开所经历的时间。

雨滴谱仪的架设高度距地面1.4 m， 取样面积54 cm²(18×3.0 cm²)， 采样间隔60 s。仪器测量的数据共有1024种粒子等级， 分为32个直径通道和32个速度通道， 分别对应的数据范围是0.2～25.0 mm以及0.2～20.0 m·s^-1。雨滴直径的有效观测范围是0.2～6.0 mm（濮江平等， 2010）， 因为第1、 2个直径通道信噪比低， 不能使用， 因此研究中剔除了观测中出现的小于0.2 mm的数据。此外， 王鹏飞等（1989）认为自然界中甚少存在>6 mm的雨滴， 濮江平等（2012）认为仪器测得的>6 mm的雨滴是由于雨滴重叠引起的， 因此， 也剔除了观测记录中出现的大于6 mm的数据（22~32个直径通道）。为保证滴谱结构的稳定性， 选取10为阈值， 剔除每分钟雨滴总个数小于10的数据文件。

使用雨滴谱数据可直接计算D_m和R:

式中: i、 j分别表示速度通道和直径通道; n(D_j)为直径在D_j到D_j+ΔD间隔内的雨滴数浓度(单位: m^-3); v_i为粒子下落速度(单位: m·s^-1) 。

图1给出了广东龙门2016年6月一次降水过程的18 h的雨滴谱仪测量结果。其中图1(a)~(b)为6月13日19:00（北京时， 下同）至14日13:00龙门降水的DSD及D_m变化情况。图1(a)表示不同直径(D)范围内的粒子数浓度随时间的变化， 色标中N(D)代表粒子数浓度（单位: m^-3·m^-1）。从图1(a)中可以看到， 13日19:00至14日01:00， 降水较少， 雨滴数浓度和D_m变化也比较平稳; 强降水发生在14日04:00 -10:00， 粒子最大直径达到6 mm， 粒子数浓度明显增加， 计算的D_m最大值超过4 mm。

为了检验雨滴谱数据的可靠性， 研究中将雨滴谱资料计算的降水率与龙门观测场雨量计直接测量的降水率进行对比[图1（c）]。从图1（c）中可以看到, 雨滴谱资料的计算结果和雨量计测得的降水率变化趋势是一致的， 在个别降水率大的时间点出现雨滴谱资料的计算结果略高于雨量计的测量结果， 但总体来看， 一致性良好， 雨滴谱仪的资料可信。

雷达波长λ_i （i=1， 2， 分别代表Ku和Ka波段）， 距离r处的观测到的雷达反射率因子Z_m可以通过雷达等效反射率因子Z_e和衰减系数k来表达:

其中Z_{e， i}和k_i可以分别表示为:

式中: σ_b、 σ_e分别是直径为D的降水粒子的后向散射截面和消光截面; K_w是介电因子， 与水的负折射指数m有关。关于N(D)， Willis（1984）提出一种“标准化”Gamma分布形式， 其中N_w不受μ的影响:

其中D_m可以表示为N(D)的4阶矩量和3阶矩量之比:

将式（6）代入式（4）和式（5）， 计算反射率因子和衰减系数。为方便计算定义I_b和I_e:

定义差分频率比DFR(单位: dB):

DFR不依赖于N_w， 因而也可以用I_b表示为:

从式（12）可以看出， 若已知两个波长的大小， 假设环境温度为20 ℃， Nw=1 m^-3·cm^-1， 给定一个μ值， 即可根据粒子米散射理论计算得到D_m与DFR的关系， 并进行存储。根据双频雷达测得的DFR及D_m与DFR的关系， 可以反演得到D_m， 将得到的D_m代入式（9）和式（10）得到I_b和I_e， 再将I_b、 I_e代入式（13）和式（14）， 即可计算出N_w和k。

并可根据下式计算降水率R:

这里N_w和D_m均由前面的反演得到， 龙门地区v(D)采用Atlas et al（1973）提出的标准大气压和温度情况下静止大气雨滴的下落末速度， 那曲海拔高（4500 m左右）， 空气密度低， 根据空气密度调整下落末速度(即在龙门地区下落末速度基础上乘以ρ0ρ0.4作为校正因子， 其中ρ0是标准大气压和温度下的空气密度， ρ是那曲地区的空气密度)。

为了从理论上计算DFR与D_m之间的对应关系， 需要计算不同降水粒子在不同电磁波长下的散射截面和衰减截面。对于非球形粒子的散射， 目前还没有一个统一的解决办法。许丽生等（2014）认为扩展边界条件法(EBCM)是严格求解具有任意大小粒子电磁散射最有效的方法之一。因此， 采用EBCM计算不同降水粒子对Ku和Ka波段电磁波的散射截面和衰减截面， EBCM的详细过程请参见Barber et al（1975）的文献。

假设温度为20 ℃， 形状因子μ分别取-2， 0， 2， 3， 4和6时， 利用公式（12）得到DFR和D_m之间的对应关系[图2（a）]。可以看到， 无论μ取何值， 在DFR>0 dB时， DFR的增大或减小完全依赖于D_m的增大或减小， 因而可以利用这种一一对应关系， 来反演D_m。但在DFR<0 dB部分， 可以发现同一个DFR的值会有两个D_m与之对应， 这样的不确定性给D_m的反演带来困难， 需要增加反射率因子的信息来帮助找到合适的D_m。由于粒子的散射特性及介电系数受温度影响， 同样利用公式（12）式， 得到μ=3时不同温度下DFR与D_m的函数对应关系[图2（b）]。可以看到， 不同温度会造成DFR-D_m关系的改变， 但这种影响主要集中在小粒子部分， DFR>0 dB部分几乎不受温度影响。因此， 为避免DFR<0时造成的D_m的反演误差， 在研究中只选取了DFR>0 dB的雨滴谱仪资料， 龙门和那曲筛选出符合要求的雨滴谱数据分别为3754和1503个。

根据雨滴谱仪的测量原理， 在某一分钟的雨滴谱样本中， 测量不同直径的粒子数浓度N(D_j)可表示为:

式中: A_ij是雨滴谱数据中第i个速度通道、 第j个直径通道的粒子数目; S为雨滴谱仪的取样面积（54 cm²）; T为取样时间（60 s）。得到粒子数浓度之后， 利用式（1）和式（2）即可用雨滴谱仪观测的DSD数据计算D_m(单位: mm)和R(单位: mm·h^-1)。衰减系数k (单位: dB·km^-1)可用衰减截面表示为:

将D_m、 R和k的计算结果作为选取数据样本的理论真值。再计算双频雷达Ka和Ku波段（λ₁和λ₂）的等效反射率因子Z_e（单位: mm⁶·m^-3）:

其中:

式中: σ_{b， i}表示波长为λ_i时粒子的后向散射截面（单位: cm²）; m表示水的负折射指数。设定雨滴谱Gamma分布中μ为不同值， 用前面介绍的双频技术的原理， 反演得到D_m、 R、 k的雷达估计值， 将雨滴谱仪计算的真值和双频雷达的反演结果进行比较， 以评估DSD参数对双频雷达估测降水的不确定性。研究中分别以广东龙门和西藏那曲的雨滴谱仪观测数据为基础进行了计算及分析， 比较广东龙门和西藏那曲地区DSD参数对双频雷达反演效果影响的差异性。具体的流程如图3所示。

为了定量比较Gamma分布模型对双频反演算法的影响， 还计算了相对误差随各参量的变化。其中， 为了尽可能的降低样本数量不足带来的误差的不稳定性， 定义绝对误差δ(x)为在(x-Δx/2,x+Δx/2)范围内的真值和估计值的平均差异， 即:

式中: N_x表示在(x-Δx/2， x+Δx/2)范围内数据点数; x_est和x_true分别代表雷达估计值和参考真值; 且(x-Δx/2， x+Δx/2）范围内的点都满足:

当x分别为D_m、 k_Ka、 k_Ku和R时， Δx分别取0.1 mm、 0.2 dB·km^-1、 0.2 dB·km^-1和2 mm·h^-1。利用δ(x)即可计算x处的相对误差:

龙门位于广东省中部， 属南亚热带季风气候， 夏季降雨量大。根据上文介绍的研究方法， 利用广东龙门地区筛选的雨滴谱仪观测资料计算了4个降水参数: D_m、 k_Ku、 k_Ka和R。并将得到的计算结果作为广东龙门地区数据样本的参考真值。图4给出的是设定μ=3时D_m、 k_Ku、 k_Ka、 R的参考真值和双频雷达反演结果的散点图。每张图有3754个点， 这些点越接近中间对角线， 表示双频雷达估计值和参考真值越接近， 说明双频雷达反演效果越好。图4表明采用μ=3的Gamma分布， 双频雷达反演结果与雨滴谱仪数据直接计算的参考真值一致性很好。随着粒子增大， D_m的数量减少变得集中， 并在大于3 mm时稍偏离对角线， 出现高估， 说明μ=3时双频雷达可能对中值区D_m反演效果更好; Ka波段衰减率在小粒子区很集中而在大粒子区发散; Ku波段的衰减率大部分点都集中在对角线上; 从降水率来看， 当μ=3时， 双频算法对小到中雨的反演效果很好， 在大雨时， 反演的降水率略低于真值。

图5给出了广东龙门地区， Gamma分布取不同μ时， 双频雷达反演的 D_m、 k_Ka、 k_Ku、 R的相对误差的计算结果。如图5所示, 当μ=-2和μ=0， 在D_m<3 mm时， 很明显被低估; 当μ=2、 3、 6和μ=4， D_m>1.5 mm时， 分别有10%左右的高估和低估， 其中μ=2、 3在D_m为1.5~3 mm时相对误差在0~10%; 对k_Ku的反演， 相对误差都在0附近， 只有μ=-2时对衰减有高估; 在对k_Ka的反演中， 除了μ=-2对雷达反演有很明显的高估， 其他μ值反演结果的相对误差都在±20%以内， μ=2、 3、 4的相对误差在±5%以内， 其中μ=3时误差相对最小， 在0值附近波动; 在降水率反演中， μ=-2降水有严重高估， μ=0和μ=6对降水分别有10%左右的高估和低估， μ=2、 3、 4时， 对降水率估计的相对误差基本在±10%以内， 其中μ=3时对5~50 mm·h^-1的降水率估计误差最小且在0值附近波动。

利用龙门地区的3754个降水雨滴谱样本， 将雨滴谱仪直接计算的降水率作为真值， 分别与利用双频技术和Z-I方法得到的降水率进行比较。其中， 当Z_e>40 dBZ时， 认为是对流性降水， 采用Z=300I^1.4; 当Z_e≤40 dBZ时， 认为是层状云降水（Steiner et al， 1995）， 采用Z=200I^1.6。计算两种方法得到的降水率的相对误差平均值δ¯。发现采用Z-I方法， δ¯ =0.47; 采用双频技术， δ¯ =0.16。表明利用双频雷达测量降水率可以提高降水率估测的准确度。

总的来说， 从广东龙门的数据结果来看， 当μ=-2时， 双频雷达对四个参数的反演效果都很差， 当μ取2、 3、 4时， 双频雷达对四个参数的反演都能得到较为稳定准确的结果， 其中μ=3时， 双频雷达对k_Ka和降水率的反演效果最优。可以看到双频雷达的反演效果与D_m、 降水率的大小也有关， 当μ=2、 3、 4， D_m在1.5~3 mm之间， 降水强度为5~50 mm·h^-1时， 反演的D_m相对误差最小、 反演的降水率最准确。Liao et al（2014）认为μ值的选择可能会受地区差异的影响。因此在研究中利用架设在西藏那曲雨滴谱仪资料继续进行了上述研究。

青藏高原作为世界上最高的高原， 被称作“世界屋脊”。独特的地理条件， 产生的动力和热力强迫对水汽和能量循环都产生巨大的影响（陈月等， 2019; 常姝婷等， 2019; 王玉琦等， 2019）， 青藏高原南部来自印度洋的暖湿气流北上过程中受高山阻留， 动力抬升， 水汽凝结， 使得青藏高原从北到南降水率差异很大。其中那曲位于青藏高原北部， 平均海拔4500 m以上， 寒冷干燥， 只有五到九月较为温暖降水量较多。Wu et al（2017）分析了高原那曲地区和我国南方的雨滴谱特征， 认为二者有很多不同之处。因此研究中利用那曲夏季雨滴谱仪观测数据， 分析了青藏高原那曲地区DSD参数对双频雷达反演降水的影响。

图6给出了那曲地区当μ=6时D_m、 k_Ku、 k_Ka、 R的反演结果和参考真值的散点分布。图6表明采用μ=6的Gamma分布来描述雨滴谱参数时， 在那曲地区， 双频雷达反演结果与雨滴谱仪观测数据直接计算的参考真值一致性较好。

为了定量比较Gamma分布模型对双频雷达反演结果的影响， 同样计算了那曲地区各降水参量反演值相对误差随各参量的变化。由于在那曲地区的观测时间较短， 计算的大滴部分相对误差波动太大， 故在计算相对误差时剔除Δx中样本数小于5个的数据。

图7给出了那曲地区Gamma分布取不同μ时， 双频雷达反演的D_m、 k_Ka、 k_Ku、 R的相对误差。在那曲地区， μ=-2对D_m严重低估， μ=0对D_m有10%~30%的低估， μ=2、 3、 4、 6反演的D_m在1.5~3 mm范围内相对误差都在10%之内， 其中μ=3、 4时相对误差最小， 在0值附近波动; k_Ku的反演效果与龙门基本没有差异; 对于那曲地区k_Ka和5～40 mm·h^-1降水率的反演， 都呈现出随着μ值的增加相对误差减小的趋势， 其中当μ=4、 6时， k_Ka相对误差最小， 在5%以内; 当μ=6时， 5～40 mm·h^-1降水率相对误差最小， 误差值在5%以内。

比较龙门和那曲双频雷达的反演效果可知， 从降水的微物理参数来看， 当μ=3时， 龙门和那曲对D_m在1.5~3 mm的反演效果都较为稳定准确; 而从宏观上看， 当μ=3时龙门地区k_Ka和5~50 mm·h^-1降水率反演的误差最小， 而那曲地区则是μ=4、 6时k_Ka相对误差最小， μ=6在5~40 mm·h^-1的降水率反演中表现更优。综上所述， 地区差异对用双频雷达反演1.5~3 mm的D_m造成的影响较小， 而利用双频算法反演Ka波段衰减和降水率时， 要达到最小的反演误差， 那曲则需要比龙门更大的μ值。

本研究的目的是探究双频雷达反演算法中DSD模型对反演结果造成的影响， 具体评价方案是基于地面雨滴谱仪观测的DSD数据， 直接计算得到不同波长的雷达等效反射率因子， 根据双频雷达理论得到反演的D_m、 衰减系数、 降水率等。将这些反演结果与用雨滴谱仪观测的DSD数据直接计算的结果进行比较， 评估DSD模型对双频雷达反演结果的影响。研究中利用广东龙门和青藏高原那曲两个气象观测场架设的雨滴谱仪获得的2~4年的数据资料， 讨论Gamma分布函数中形状因子μ取不同值时双频雷达反演D_m、 k_Ku、 k_Ka、 R 4个参数的反演误差， 并比较了其在那曲和龙门地区的异同， 得到以下结论:

（1） 当μ=-2时， 双频雷达在两个地区对所有参数的反演效果都很差。当μ取2、 3、 4时， 双频雷达都能对龙门降水四个参数的反演都得到稳定准确的结果。其中μ=2、 3时， 在D_m为1.5~3 mm时相对误差在0~10%; μ=2、 3、 4时， k_Ka反演的相对误差在±5%以内; 除μ=-2时， k_Ku反演的相对误差值都在0值附近; μ=3时对5~50 mm·h^-1的降水率估计误差基本在0值附近波动。这说明在华南地区μ取2、 3、 4时双频反演技术都能达到良好的反演效果。

（2） 对龙门降水反演结果的误差分析中发现: 双频雷达的反演效果与D_m、 降水率的大小范围也有关系。

（3） 对比双频雷达技术在那曲和龙门的反演效果: μ取3时两个地区反演1.5~3 mm D_m的误差都很小， 但μ=3在龙门反演D_m的相对误差在0~10%浮动， 相应的μ=3在那曲反演D_m的相对误差则在0值附近波动; μ取3时龙门地区k_Ka和5~50 mm·h^-1降水率反演的误差都最小， 而那曲地区则μ=4、 6时和μ=6在k_Ka和5~40 mm·h^-1的降水率反演中相对误差最小。在k_Ka和降水率的反演中， 那曲需要比龙门更大的μ值， 这说明双频雷达技术在青藏高原的应用有其特殊性， 需要与平原地区区别对待。

事实上， 要提高降水率的反演准确度， 更需要假设的Gamma分布函数贴近当地实际雨滴谱分布， 而之前有学者对青藏高原雨滴谱分布的研究中发现青藏高原雨滴谱分布和平原地区有很多不同之处(Wu et al， 2017; Porcù et al， 2014; 常祎等， 2016)。因而在今后的工作中我们将对当地实际雨滴谱分布进行拟合, 将拟合结果与Gamma分布模型做对比, 寻找双频雷达反演降水不准确的根本原因。此外我们也将对龙门的双频雷达探测一致性进行对比, 并利用雷达实测数据开展液水含量和雨滴谱反演的个例研究。