1、收稿日期：2021-05-17；修订日期：2022-07-08基金项目：国家自然科学基金面上基金(41975027)；国家重点研发计划重点专项(2017YFC1501401)共同资助通讯作者：寇蕾蕾，女，湖北省人，副教授，主要从事气象雷达数据处理和应用研究。E-mail：第38卷 第6期2022年12月热 带 气 象 学 报JOURNAL OF TROPICAL METEOROLOGYVol.38,No.6Dec.,2022陈垚，寇蕾蕾，蒋银丰，等.基于新型观测算子的双偏振雷达雨滴谱变分反演J.热带气象学报，2022，38(6)：854-869.文章编号：1004-4965(2022)06-0
2、854-16基于新型观测算子的双偏振雷达雨滴谱变分反演陈垚1,2,4,5，寇蕾蕾1,3，蒋银丰1，杨春生1，林正健1，楚志刚1,3(1.南京信息工程大学大气物理学院，江苏 南京 210044；2.福建省平潭综合实验区气象局，福建 平潭 350400；3.南京信息工程大学气象灾害预警与评估协同创新中心，江苏 南京 210044；4.平潭海洋气象野外科学观测研究站，福建 平潭 350400；5.福建省灾害天气重点实验室，福建 福州 350001）摘要：提出一种基于变分理论的双偏振多普勒天气雷达反演雨滴谱方法，最优反演雨滴谱的同时可实现衰减订正。反演过程中使用一种新型观测算子，利用滴谱仪实测数据计算
3、状态变量和双偏振参量，使观测算子更能代表本地降水特性。将新型观测算子、误差协方差矩阵和状态变量的先验估计用于代价函数中，基于高斯-牛顿迭代方法求解代价函数得到最优雨滴谱。利用理想模拟试验和南京信息工程大学C波段双偏振多普勒天气雷达实测个例对算法进行验证和评估。结果表明：算法反演得到的状态变量(液态水含量LWC和质量权重平均直径Dm)与滴谱仪数据计算结果的相关系数达到了0.96和0.80，相对偏差为25.19%和10.63%，均方根误差比常规反演结果改善了50%左右，比基于模拟观测算子的变分反演结果改善了30%左右，最优反演雨滴谱得到的降雨率R和雨量计数据相关系数达到0.89，相对偏差为14.7
4、8%。关键词：变分反演；观测算子；双偏振多普勒天气雷达；雨滴谱反演；衰减订正中图分类号：P406文献标志码：ADoi：10.16032/j.issn.1004-4965.2022.0771 1 引引言言雨是最常见的降水形式。雨滴谱(RaindropSize Distribution，DSD)是描述降雨的一个重要参量，它与雨滴形状、取向以及下落末速度一起，能够提供较完整的降雨信息，并可以计算雨滴散射振幅、双偏振雷达参量和降雨参数。降雨参数主要包括降雨率(R)、液态水含量(LWC)、质量权重平均直径(Dm)、体积中值直径(D0)、归一化截距参数(NW)和总数密度(Nt)等，可用于雷达定量降水估测、
5、全球水循环分析和数值预报雷达资料同化等研究。通过分析不同类型降水的降雨参数，可更加了解降水微物理过程，改善降水微物理参数化，因此雨滴谱和降雨参数的精确反演对准确监测和预报降水十分重要。双偏振雷达通过发射垂直和水平方向的电磁波来实现探测，可获得两个通道的回波信息。它不仅能获得常规多普勒雷达的参数如水平反射率因子(ZH)、径向速度(v)和谱宽(SW)，还可得到双偏振观测量如差分反射率因子(ZDR)、比差分传播相移(KDP)、总差分相移(DP)、协相关系数(hv)等，这些观测量能够反映降水系统内部水凝物大小、形状、相态、浓度等信息，所以双偏振天气雷达的观测量可用于雨滴谱反演，进而计算降雨参数，有助于
6、更好地研究我国的降水微物理特性并定量估测降水1-4。目前，雨滴谱反演主要有两类方法。第一类是常规的雨滴谱反演，主要通过先建立雨滴谱模型，再推导双偏振参量和雨滴谱参数的关系式，实现雨滴谱反演。Seliga 等5在研究中建立了指数雨滴谱模型，并由ZH、ZDR和KDP确定模型参数后得到雨滴谱。Ulbrich等6提出了 gamma分布模型，并将双偏振雷达参量和降雨特性联系在一起。第6期陈垚等：基于新型观测算子的双偏振雷达雨滴谱变分反演Zhang等7提出利用-关系约束gamma分布，得到两参数模型(constrained-gamma model，简称C-G模型)8-12，由ZH和ZDR计算C-G模型参数
7、，成功反演雨滴谱。国内也有许多学者利用常规方法反演雨滴谱，如张鸿发等13通过统计分析得到 ZH、ZDR与雨滴谱参数和雨强的关系，用于反演雨滴谱。李宗飞等14通过对雨滴模型的散射模拟以及gamma分布的拟合，建立雨滴谱参数与双偏振参量之间的关系反演雨滴谱。Liu等15用广东惠州龙门站滴谱仪数据计算-关系，并基于适用于华南地区的C-G模型，用ZH和ZDR反演雨滴谱。第二类方法是基于最优化理论的雨滴谱反演，主要将最优化理论和常规的雨滴谱反演方法相结合，或通过建立观测算子，基于最优化算法求解包含观测 误 差 在 内 的 非 线 性 方 程 组 反 演 雨 滴 谱。Vulpiani等16基于神经网络算法
8、由S波段双偏振雷达数据计算雨滴谱。Cao等17-18先利用贝叶斯理论反演雨滴谱，用雨滴谱的先验估计降低观测误差，之后又基于C-G模型建立观测算子，利用二维变分理论反演雨滴谱。Yoshikawa等19基于极大似然法由X波段雷达双偏振观测量反演雨滴谱。Wen等20提出了一个由双偏振观测量反演雨滴谱的模型，首先基于C-G模型建立观测算子，其次构造一个伪训练集，利用最邻近算法反演雨滴谱。Vivek等21和Huang等22分别基于指数模型和C-G模型建立观测算子，利用S波段和C波段雷达数据变分反演得到最优雨滴谱。与观测算子同理，Sun等23建立了雨滴谱模型参数和雷达观测量之间的逆映射表，由X波段偏振雷达
9、数据反演雨滴谱。无论是常规反演还是变分反演，都需要在反演前建立雨滴谱模型，利用模型参数表示雨滴谱。但在实际情况中，雨滴谱是复杂且多变的，某一种雨滴谱模型并不能代表所有的雨滴谱分布。前人在变分反演中所用的观测算子多为建立雨滴谱模型后模拟得到，本研究不建立雨滴谱模型，使用滴谱仪实测数据计算更能反映当地降水特点的观测算子，将其应用于双偏振雷达雨滴谱最优化反演中，进一步提高了雨滴谱反演的准确性。为反演得到更精确的雨滴谱，本文提出了一种基于新型观测算子的双偏振雷达变分反演雨滴谱方法。首先，介绍本文所用到的C、S波段双偏振多普勒天气雷达数据和Parsivel滴谱仪数据；其次，描述基于新型观测算子的变分反演
10、方法，并阐述观测算子的建立；随后，阐述基于新型观测算子变分反演雨滴谱的详细过程；最后，利用一个理想模拟试验和NUIST-CDP观测到的实测个例验证该方法的可行性，并在文末分析了结果优化的原因。2 2 数据介绍数据介绍2.1 C波段双偏振多普勒天气雷达本研究基于南京信息工程大学 C 波段双偏振多普勒天气雷达(简称 NUIST-CDP)数据最优反演雨滴谱。雷达站点坐标为 118.717 2 E，32.206 9 N，天线直径11 m，扫描半径150 km，脉冲宽度0.5 s，距离库长75 m，波束宽度0.54，完成一次体扫所需时间为78 min，可探测得到水平反射率因子(ZH)、径向速度(v)、速
11、度谱宽(SW)、差分反射率因子(ZDR)、线性退偏振比(LDR)、差分传播相移(DP)、差分传播相移率(KDP)、协相关系数(hv)等参量。本研究在反演之前对雷达数据进行如下预处理：基于地物杂波识别算法24去除雷达数据中的地物杂波。利用 21 点中值滤波器平滑 ZH、ZDR、DP，消除原始数据中出现的异常跳变值。对DP进行退模糊，当相邻两点的 DP值之差大于320，若前一点大于0，则将后一点加上360，反之减360，模糊处理后的DP再用21点中值滤波器平滑，最后用最小二乘法计算KDP。因为信噪比和hv较小时雷达数据质量较差25，所以在实际应用中剔除了信噪比小于20 dB、hv小于0.9、ZDR
12、和KDP为负值的雷达数据。本研究在使用常规方法反演雨滴谱时，对雷达数据进行了如下衰减订正26：当 0.1/kmKDP3/km时，使用 AH和 KDP的关系进行订正(该关系由滴谱仪数据拟合得到，下同)，订正公式为AH=0.092KDP；当KDP0.1/km或KDP3/km时，使用AH和ZH的关系订正，订正公式为AH=1.43610-5ZH0.855(ZH的单位为 mm6/m3)。2.2 S波段多普勒天气雷达由于本研究的最优雨滴谱反演可同时实现衰减订正，南京龙王山S波段多普勒天气雷达(简称SA雷达)反射率因子 PPI(Plane Position Indicator)图像主要用于检验衰减订正的结果
13、。SA雷达站点坐标为 118.696 9 E，32.190 8 N，距离库长 1km，波束宽度约为1，完成一次体扫所需时间为855热 带 气 象 学 报第38卷6 min，其距离NUIST-CDP约2.5 km，二者高度相差约50 m，对比时选取二者体扫起始时刻最接近的数据。2.3 OTT Parsivel滴谱仪本研究使用滴谱仪数据计算结果来建立观测算子、模拟雷达径向廓线和检验反演结果，数据来自于三部德国 OTT 公司生产的 Parsivel 滴谱仪。第一部滴谱仪位于南京信息工程大学大气综合观测基地，站点坐标为118 42 E，32 12 N，距离NUIST-CDP 约 1 800 m，相对雷
14、达的方位角为250。第二部滴谱仪位于高淳站点，坐标为118.903 9 E，31.333 3 N，距离NUIST-CDP约98km，相对雷达的方位角为170。第三部滴谱仪位于江宁站点，坐标为118.899 7 E，31.931 4 N，距离 NUIST-CDP 约 35 km，相对雷达的方位角为151。在反演之前对滴谱仪数据进行如下预处理：用6点中值滤波器平滑数据，当观测通道的雨滴个数小于10个或计算出的降雨率小于0.5 mm/h，则剔除该点数据。本研究从滴谱仪数据计算双偏振雷达观测量和降雨参数，先从滴谱仪数据中读取出滴谱信息即N(D)，再利用T矩阵方法27求得C波段粒子散射系数，进而计算双偏
15、振雷达观测量和降雨参数。2.4 雨量计本文所用的雨量计数据来源于南京市地面雨量自动站，共计85个雨量站点。雨量计数据时间分辨率为 1 小时，雨量分辨率为 0.1 mm。本文以雨量计数据为实际雨强来验证最优雨滴谱计算降雨率的准确性。NUIST-CDP、SA雷达、三部滴谱仪以及雨量计站点的空间分布如图1所示。3 3 雨滴谱反演方法雨滴谱反演方法3.1 代价函数基于变分理论的双偏振雷达雨滴谱反演本质上是求状态变量的先验估计、背景误差、观测误差、双偏振观测量和观测算子预测结果组成的非线性方程组最优解的过程。因为雨滴谱参数和Dm、LWC之间可相互转化，且Dm和LWC符合高斯分布21，本研究使用Dm和LW
16、C作为状态变量，建立双偏振观测量和雨滴谱之间的关系。将变分反演雨滴谱的代价函数定义为：J(x)=(x-xb)TB-1(x-xb)+y-H(x)TR-1 y-H(x)(1)式(1)中，x 为每条雷达径向上 n 个距离库 Dm和LWC的向量形式：xDm=Dm(1),Dm(2),Dm(n)T(2)xLWC=LWC(1),LWC(2),LWC(n)T(3)x=xDm;xLWC(4)式(2)(4)中，Dm(n)和 LWC(n)代表第 n个距离库上的状态变量，当状态变量确定时，可用观测算子预测双偏振观测量。xb为状态变量的先验估计(背景场)，由状态变量和双偏振参量之间的经验关系计算得到。B为背景场的误差协方差矩阵，通过对状态变量进行南京地区气候统计分析后，将背景误差标准差 Dm和 LWC分别设置为 1 mm 和0.707 g/cm3(方差分别为1 mm2和0.5 g2/cm6)，本研究假设了一个高斯相关模型来计算协方差矩阵28：Bij=b2exp-12(rijrL)2(5)式(5)中，b2是背景误差协方差(Dm2和 LWC2)，rij是第i个和第j个雷达距离库之间的距离，rL是空间去相关长度，会影
