引言

上海地处欧亚大陆东岸，长江三角洲前缘，濒临太平洋，季风气候特征明显。受中纬度系统与副热带系统交替影响，上海市暴雨、雷暴大风、冰雹、台风等气象灾害频发(曹晓岗等，2014；陈雷等，2015)。天气雷达作为短临监测和预警的重要手段之一，在业务中的作用日益凸显(严红梅等，2019)。上海市气象局于2014年完成南汇WSR-88D雷达的双线偏振升级改造工作，相较单偏振雷达，升级后的双偏振雷达能够获取差分反射率因子Z_dr，差分相位ϕ_dp、相关系数ρ_hv等更多独立变量。丰富的偏振量在回波的相态识别、暴雨和冰雹等灾害性天气预警及定量估测降水等方面具备重大的应用价值(程周杰等，2009；尹春光等，2011；孙丝雨等，2013；杨祖祥等，2019)。

差分相移率K_DP作为差分相位ϕ_dp的二次产品是双偏振雷达的重要变量之一，由于其不受地物阻挡和回波衰减影响的特性，在粒子相态识别和定量估测降水等方面具备显著优势。差分相移率K_DP的数据质量主要由差分相位ϕ_dp决定，关于差分相位ϕ_dp的质量控制，国内外已有众多学者进行了相关研究工作。Hubbert等(1993, 1995)针对C波段双线偏振雷达，采用低通滤波器对差分相位ϕ_dp进行滤波，并在此基础上设计了有限冲击响应(FIR)滤波器通过多次迭代滤波进行差分相位ϕ_dp的质控；Ryzhkov(2007)针对双偏振雷达非均匀充塞(Non-Uniform Beam Filling，简称NBF，下同)的误差推导了WSR-88D雷达适用的订正公式；杜牧云等(2011, 2012)利用一次降水过程对C波段双线偏振进行了偏振参量的评估，并利用小波分析法进行了差分相位ϕ_dp质量控制的初探；肖艳姣等(2012)在退相位折叠和初相位调整后，利用FIR滤波器进行X波段双线偏振雷达差分相位质量控制工作。曹杨等(2016)利用质控后的C波段双偏振雷达差分相位ϕ_dp计算差分相移率K_DP，并在一次暴雨过程中进行了验证。

但上述方法的适用领域多集中于C波段和X波段双线偏振雷达，对目前业务主要应用的S波段双线偏振雷达涉及较少。本文在WSR-88D雷达差分相位ϕ_dp质控算法的基础上(Bringi and Chandrasekar, 2010；Branch，2012；Kumjian，2013)，针对影响差分相位ϕ_dp数据质量非均匀充塞(NBF)问题，提出了针对NBF径向进行识别并订正径向相关系数，进而根据设定阈值进行回波差分相移率分类计算的自适应方法，以期对业务雷达的差分相位以及差分相移率数据质量进行控制，提高其在降水估测中的可用性。

1 雷达参数

上海市气象局业务用S波段双偏振雷达为洛克希顿-马丁公司生产的WSR-88D型雷达, 于1997年从美国引进，布设于浦东新区滨海镇(原南汇滨海镇)，2014年进行双偏振改造升级(以下简称南汇双偏振雷达)。升级后的南汇双偏振雷达LevelⅡ数据(基数据)提供反射率因子、径向速度、速度谱宽、差分反射率因子、相关系数、差分相位共6个基本产品。WSR-88D雷达主要技术指标见表 1。

2 算法介绍

本文设计的算法主要包括差分相位质量控制和差分相移率计算两部分。具体实现如下：首先进行退相位折叠操作，其次针对是否NBF径向选择不同判据区分径向上距离库内的气象回波和非气象回波，并对非气象回波进行插值，最后利用9点和25点最小二乘法自适应计算差分相移率，整体算法流程如图 1所示。本节将具体说明算法的详细计算步骤。

2.1 差分相位质量控制

差分相位的质量控制工作主要包括退相位折叠、非均匀充塞径向识别、非气象回波识别和相位滤波四个方面。

2.1.1 退相位折叠

本文所用的WSR-88D双线偏振雷达最大不模糊相位为360°，退相位折叠主要包括如下步骤：首先将径向上连续30个距离库视作一个数据块，以同时满足相关系数大于等于0.9和存在有效差分相位数据以及差分相位标准差小于120°三个条件作为判别标准，统计数据块内有效点的个数。若有效点个数大于等于15认为该部分为有效数据块，该数据块内所有有效点的中值为ϕ_{dp_median}。如公式(1)—(3)所示，计算中值滤波差值和单次、二次折叠差值(分别用Δϕ_dp(i)、Δϕ_dp1(i)、Δϕ_dp2(i)表示)

$ \Delta \phi_{\mathrm{dp}}(i)=\left|\phi_{\mathrm{dp}}(i)-\phi_{\mathrm{dp}_{-} \mathrm{median}}\right| $

(1)

$ \Delta \phi_{\mathrm{dp} 1}=\left|\phi_{\mathrm{dp}_{\_\mathrm{median}}}-\left(\phi_{\mathrm{dp}}(i)+360\right)\right| $

(2)

$ \Delta \phi_{\mathrm{dp} 2}=\left|\phi_{\mathrm{dp}_{\_\text {median }}}-\left(\phi_{\mathrm{dp}}(i)+360 \times 2\right)\right| $

(3)

若Δϕ_dp(i)> Δϕ_dp1(i)，则认为发生单次差分相位折叠，利用单次相位折叠退折叠（ϕ_{dp_unfolded}）公式(4)进行计算

$ \phi_{\mathrm{dp}{\_\text {unfold }}}=\phi_{\mathrm{dp}}(i)+360 $

(4)

若Δϕ_dp(i)> Δϕ_dp2(i)，则认为发生了差分相位二次折叠，退相位折叠公式如公式(5)所示

$ \phi_{\mathrm{dp}\_{\text {unfold }}}=\phi_{\mathrm{dp}}(i)+360 \times 2 $

(5)

2.1.2 非均匀充塞(NBF)径向确定

在距离雷达较远处存在垂直发展较高回波单体时，雷达波束展宽会造成取样体积内非均匀充塞(NBF)现象，即取样体积内存在不同相态的降水粒子(如较高高度处为冰相的冰晶粒子，较低高度处为液态的降水粒子)，从而不满足取样体积内粒子均匀分布且相态单一的假设条件。进而造成雷达探测物理量的偏差，而其中又以反射率和相关系数对其最为敏感。故在进行径向上距离库内的气象回波和非气象回波识别前，需对是否NBF径向分情况讨论，本文确定NBF径向的标准为：自径向180个距离库开始(距离雷达距离较近处发生NBF情况较少)，计算单根径向满足反射率在30~50 dBz之间，速度绝对值大于1 m·s^-1并且相关系数小于0.7的点的个数。若满足条件的点总数大于10则认为对应径向为NBF径向。

2.1.3 气象回波识别及差分相位滤波

由于NBF径向的相关系数可靠性降低，故NBF径向利用信噪比和差分相位参量识别气象回波和非气象回波；非NBF径向使用相关系数和差分相位信息。由于本文使用的WSR-88D雷达LevelII数据中不包含信噪比参量，对于NBF径向首先利用均值滤波后(本文选用径向3点均值滤波)的反射率因子数据计算信噪比参量(SNR_hc(i, jc))，推导过程中进行了无量纲处理。如公式(6)所示

$ S N R_{\mathrm{hc}}(i, j c)=Z_{{\rm c}}(i, j c)-20 \log (R(i))-C $

(6)

Z_c(i, jc)为均值滤波后的反射率数据(单位: dBz)，R为径向距离(单位:km)，C为雷达常量(基数据头文件提供，单位:dB)

气象回波和非气象回波的识别判据为：针对NBF径向，若信噪比大于5 dB且对应差分相位为有效值，则认为该距离库为气象回波；针对非NBF径向，若对应距离库内相关系数大于等于0.9且存在差分相位数据，则认为该距离库为气象回波。

由于差分相移率反映的是差分相位在一定距离上的变化情况，为方便下一步的差分相移率计算，对径向上的非气象回波进行插值填补。在对径向差分相位数据进行5点中值滤波的基础上，分别利用9点和25点平滑及插值方法进行非气象回波的插值工作。以9点平滑和插值为例，首先对径向差分相位进行9点均值滤波(对应位置前4个和后4个及其本身有效相位数据进行滤波)，再选取9点内的有效数据拟合线性曲线并进行非气象回波差分相位的插值。

2.2 差分相移率计算

差分相移率(K_DP)反映的是差分相位在单位距离上的相位变化，其形式如公式(7)所示。差分相位质量控制完成后，根据径向单点的相关系数和反射率因子值，分别选择9点和25点最小二乘法计算差分相移率。反射率大于40 dBz并且相关系数大于等于0.9的点选用9点最小二乘法计算；反射率小于40 dBz且相关系数大于等于0.9的点选用25点最小二乘法计算，其计算如公式(8)

$ K_{\mathrm{DP}}=\frac{\phi_{\mathrm{dp}}\left(r_{2}\right)-\phi_{\mathrm{dp}}\left(r_{1}\right)}{2\left(r_{2}-r_{1}\right)} $

(7)

$ a=\frac{\sum x y-\frac{1}{N} \sum x \sum y}{\sum x^{2}-\frac{1}{N}\left(\sum x\right)^{2}} $

(8)

由于NBF径向的相关系数受粒子相态非均匀分布的影响，相关系数会受到低估。本文选用Ryzhkov针对WSR-88D型雷达推导的相关系数订正公式对NBF径向的相关系数进行订正，而后再进行相关系数阈值判断，如公式(9)所示

$ \frac{\left|\rho_{\mathrm{hv}}^{(m)}\right|}{\rho_{\mathrm{hv}}}=\xi_{1}=\exp \left\{-1.37 \times 10^{-5} \Omega^{2}\left[\left(\frac{\mathrm{d} \phi_{\mathrm{DP}}}{\mathrm{d} \theta}\right)^{2}+\left(\frac{\mathrm{d} \phi_{\mathrm{DP}}}{\mathrm{d} \phi}\right)^{2}\right]\right\} $

(9)

其中ρ_hv^(m)为探测值，ρ_hv为实际值(订正值)，Ω为3 dB波束宽度的实测值(WSR-88D型双线偏振雷达为0.97°)，θ和ϕ分别对应雷达的仰角和方位角。

3 效果检验

2018年9月20日，受蒙古气旋东移以及南支槽北抬的影响，华东地区出现一次大范围降水天气。本节以该次过程双偏振雷达的观测资料为例说明本算法在实际过程中的应用。

3.1 退相位折叠

退相位折叠是差分相位质量控制中的首要任务，直接影响到后续的相位插值平滑以及差分相移率的计算工作。图 2为原始差分相位PPI与退相位折叠差分相位PPI的效果对比。由图 2a中可知，椭圆圈内，沿径向的差分相位数值由接近360°突变到20°左右。结合相邻径向的差分相位分布，可以确定椭圆圈内发生沿径向方向的差分相位折叠。需要根据前文所述的内容进行退相位折叠处理，其结果如图 2b，经以上处理后超过360°的差分相位获得了有效的订正，满足差分相位在沿径向方向上相位增加的经典概念模型。

下文以206°径向详细阐述本文的退相位折叠效果，如图 3所示，在径向160 km到190 km距离内差分相位发生了单次相位折叠(差分相位差为360°)。退折叠后的差分相位沿径向整体上呈递增趋势。由于差分相移率主要表现径向距离上相位的变化，若不进行退折叠处理，对于后续的差分相移率计算影响巨大。

3.2 非均匀充塞径向识别

在距离雷达中心径向距离较远处，雷达取样体积内会观测到不同相态的降水粒子，不同相态降水粒子位于同一取样体积内会间接造成相关系数的低估。图 4为NBF识别所需的径向特征量及识别出的NBF径向。

如图 4a所示的0.5°仰角反射率PPI，黑色椭圆内的回波位于强回波(回波强度50 dBz以上)单体后部，反射率强度值位于30~50 dBz之间。对应位置处的相关系数数值大多低于0.8，较远距离处的相关系数值甚至低于0.7。结合4c径向速度图PPI黑色椭圆内径向速度普遍大于2 m·s^-1的实际，可以初步认定黑色椭圆部分为气象回波，但由于空间粒子非均匀充塞(NBF)引起的反射率值和相关系数值的低估。

通常针对非NBF径向，以相关系数0.9作为阈值能够区分绝大多数的气象和非气象回波，而针对NBF径向，相关系数的低估导致其不能作为区分降水和非降水回波的有效判据。在本文的研究中，对于NBF径向的回波数据，参照美国的业务试验研究成果(Ryzhkov, 2014)，选用信噪比(SNR，阈值5 dB)代替相关系数进行气象和非气象回波的初步区分。相较NBF径向，均匀填充径向较为普遍，因此本文选择NBF径向数据阐述初步识别径向气象和非气象回波的基本原理。以206°径向(NBF径向)为例，识别后的气象回波如图 5所示。

径向距离约140~200 km范围内，受非均匀充塞影响，相关系数抖动较大，并且数值多低于0.9，选用一般的判别标准(相关系数阈值0.9)会导致气象回波的过度剔除，进而影响后续的差分相位插值并最终导致差分相移率数据的缺失。针对NBF径向选用信噪比作为判别条件能够在气象回波和非气象回波筛选过程中尽可能保留更多的回波信息，为后续的差分相位插值和差分相移率的计算提供更多的数据支持。

3.3 差分相位插值

对退相位折叠后的差分相位数据进行5点中值滤波，再分别使用9点和25点滤波窗对上述中值滤波结果进行均值滤波。随后根据气象回波和非气象回波的初步识别结果，对非气象回波数据块进行径向9点和径向25点的线性插值。即通过拟合径向距离库前后4个点(12个点)有效差分相位与径向距离的线性关系式插值补充非气象回波数据块。径向差分相位的插值结果如图 6所示。

由图 6可见，在径向140-180 km范围内，9点和25点插值相位曲线趋势保持一致，证明插值方法能够正确反映电磁波通过雨区后相位变化的规律。9点插值曲线在强回波位置处的变化较为明显，在随后的差分相移率计算中不至于平滑掉过多局部变化的趋势，保留强回波处相位变化的真实信息。25点插值曲线则较为平滑，可滤除弱降水回波处差分相移的个别变化较大的不合理点，有利于弱回波区域的差分相移率计算。对于滤波窗范围均为无效相位点的情况，利用本根径向最近距离处的有效相位点代替。对于整根径向均为无效相位点的情况，采用系统初相位(本文选用的WSR-88D雷达系统初相位值为60°)代替整根径向的差分相位值。

3.4 差分相移率计算

差分相移率的最终计算结果是结合相关系数和反射率数值选择9点和25点最小二乘法得到。针对NBF径向，需要对相关系数进行订正，从而在噪声较多数据中保留更多真实气象回波信息。本文选用Ryzhkov提出的针对WSR-88D型号雷达相关系数订正方法进行NBF径向相关系数的订正工作，订正结果如图 7所示。

在径向发生非均匀充塞的区域(140~200 km)，相关系数得到了适度的订正。经订正后，黑色框内原本位于相关系数位于0.9以下的点经订正后相关系数位于0.9以上(黑色方框内)。虽出现了过度订正情况(黑色圆圈内)，但也仅限于个别情况。综合而言，该订正算法对于非均匀充塞波束内的气象回波差分相位真实信息的保留起到了有利的作用，在径向140~200 km范围内保留了更多的气象信息。

关于差分相移率的计算，针对NBF径向，选择订正后的相关系数作为差分相移率是否进行计算的判别条件，即相关系数大于等于0.9的点视为有效点进行差分相位计算。针对非NBF径向，选取5点均值滤波后的相关系数作为判据。再结合反射率因子是否大于40 dBz，选择9点或25点最小二乘法计算差分相移率。

图 8给出了单根径向(206°)差分相移率的计算结果以及对应的反射率因子信息。可见，在径向150 km以内的范围内，差分相移率曲线的波动趋势与径向反射率因子的波动趋势较为一致。75 km以内由于回波强度不强，差分相位的值变化较小，选择25点最小二乘法计算得到较为平滑的差分相移率曲线。而75~ 150 km范围内选用9点最小二乘法计算保留差分相移率的变化趋势。150 km以外受粒子非均匀充塞影响，通过对相关系数进行订正尽量保留有效的差分相位信息计算差分相移率。

图 9a给出最终计算得出的差分相移率PPI，图 9b为WSR-88D雷达的K_DP产品(LevelⅢ数据)。对比图 9a中NBF径向相关系数订正后的0.5°差分相移率PPI与同仰角反射率因子PPI(图 4a)，可以发现，反射率强中心的区域在差分相移率PPI上均表现为较大的值与之对应。对于NBF径向部分(图 4d)，与WSR-88D的差分相移率产品局部放大图对比(图 9b，与图 9a黑色框内部分对应)可知，经相关系数订正后的差分相移率计算结果(图 9c)在径向上尽可能还原了WSR-88D算法直接设为无效值的部分(原始相关系数小于0.9的部分)，扩充了差分相移率K_DP的可用性。

4 结论与讨论

(1) 在业务用S波段双偏振雷达中，差分相位折叠是影响差分相位数据质量及后续差分相移率计算的重要因素之一，特别是在大范围强降水区域后侧，其相位折叠明显，需要进行退相位折叠，才能进一步应用该数据。

(2) 非均匀充塞(NBF)径向会导致相关系数的低估，影响气象回波和非气象回波的判别，可能导致该类径向上远距离的气象回波数据被误判而剔除，因此在进行差分相位质量控制时需要首先识别出NBF径向，并利用信噪比代替原有不合理的相关系数参数，进行NBF径向气象和非气象回波的识别。

(3) 针对NBF径向，采用合理的订正公式可以有效解决NBF径向中相关系数低估的问题，从而能够保留更多原来被错误剔除的气象回波信息，使之能参与其他物理量的计算。

(4) 根据一定的反射率及相关系数阈值，合理区分强回波区和弱回波区，再根据区分情况选用9点和25点进行插值和最小二乘计算差分相移率，能够在弱回波区域内减少差分相位变化过大的不合理点影响，而在强回波区域体现相位的快速变化，从而使得计算的差分相移率更加接近真实情况。

双偏振雷达由于采样体积内的粒子相态非均匀充塞会造成反射率及其他偏振量的异常，但同时非均匀充塞现象能否定量说明冰相转换层的大致高度？未来需要收集更多非均匀充塞的个例并与对应时次的探空数据进行对比。