孙 跃，肖 辉，冯 亮，杨慧玲，陈 程，周筠王君

(1.成都信息工程大学大气科学学院，四川成都610225;2.中国科学院大气物理研究所云降水物理与强风暴重点实验室，北京100029)

差分传播相移ΦDP是双偏振雷达的一个重要的偏振参量，其不仅可以被用于X-band双偏振雷达的衰减订正[1－3]当中，由 ΦDP得到订正的差分传播相移率KDP，还可以被应用到双偏振雷达降水估算[4]、相态识别[5]当中。但是，ΦDP原始信号中相对于上述用途而言包含了由雷达本身造成的系统噪声和地物可能引起的观测噪声，因此在使用之前需要先对其进行质量控制。对于ΦDP的质量控制，国内外学者分别采用过滑动平均滤波[6－7]、中值滤波[6－8]、Kalman 滤波[2，6，8]、小波滤波[6，8－9]、迭代滤波[10]、分类处理法[7]、FIR 低通滤波[6，8]等多种算法。综合上述相关文献所列举的滤波个例看，当质量控制目标对应的是连续的云体，且差分传播相移ΦDP中只存在零星地物噪声和幅度变化不大的系统噪声时，用上述算法便能取得较为理想的质量控制结果。但是，当扫描径向上存在不只一个云体，尤其是对雷暴群扫描时，ΦDP在云体的边缘位置会产生异常的震荡且在云体间隔的位置不连续，或者当观测噪声较大以及系统噪声较大时，使用上述方法均不能得到理想的ΦDP质量控制效果，其表现是大幅度杂波只被控制成了小幅度杂波，而没有完全被滤除。这主要是因为，上述方法大都无法在算法层面将有利用价值的ΦDP变化信息与脉冲式的杂波完全区分开。杜牧云等[7]提出分类处理法虽然对此问题有了一定的改进，但是其剔除大幅度杂波的算法首先是基于连续性检查和大量历史资料统计得到的，当径向数据不连续或者大幅度杂波连续出现时，基于统计意义上的连续性检查可能无法正确地从中分离出有价值的ΦDP信号。因此，进一步研究更实用的ΦDP质量控制算法是有必要的。

值得注意的是，以往的ΦDP质量控制方法均采用一般性的信号处理算法，但ΦDP有别于一般模拟信号/数字信号的特点。简而言之，在ΦDP信号中对气象研究与应用有重要价值的信息不是波动信号或脉冲信号，多数情况下也不能用单一的解析方式描述。从上述研究中成功的滤波个例[2，6－8，10]来看，有使用价值的ΦDP信号在其扫描经向上更像是一个导数不为负的连续分段函数，而这种描述也非常符合ΦDP信号的定义。作为一个距离积分量，ΦDP在粒子直径较大的探测距离段上增长较快，而在粒子直径较小的距离段上增长缓慢，甚至有小的波动。从双偏振雷达的衰减订正、雨量估算、相态识别等用途看，也需要使ΦDP曲线的斜率尽量不为负。同时，在以往成功的ΦDP质量控制个例中，ΦDP的连续上升段呈现出明显的线性增长特征。因此，ΦDP的质量控制可以由一个单纯的降噪滤波问题转化为一组导数不为负的连续线段的确定问题。基于以上论述，设计了一种针对ΦDP特点的“线性拟合-递推法”的质量控制方法。

1 线性拟合-递推法介绍

为了说明线性拟合-递推法的具体应用，下面以中国科学院大气物理研究所车载X波段双偏振雷达2008年6月23日在首都机场观测的多单体强对流天气为例[11]，选取14:59 分方位角 296°RHI扫描 0.34°仰角的雷达水平偏振反射率ZH与差分传播相移ΦDP作为研究对象(图1a)(径向库长150 m)。从图1(a)的ZH曲线可以看出，在扫描线上前400个距离库上有两个主要云体对应的两段回波，而两段回波并非是连续的。因此需要先将对应两段的ΦDP进行拼接(合并)，即填补中间不连续的数据段，再按照介绍的方法对ΦDP进行质量控制。具体步骤如下:

(1)滑动线性平滑

采用的滑动线性平滑的原理类似于经典的多项式平滑滤波，如五点二次滤波、五点三次滤波[12]，即在滑动窗口D范围内通过最小二乘法求解该段的线性拟合系数aj与bj，进而求得点j处的线性平滑结果 ^yj。

其中xj代表第j个距离库距雷达原点的距离，系数aj与bj的计算公式即在滑动窗口内利用典型的最小二乘法[13]求解一次多项式系数的公式

xl与yl分别代表第l个点距雷达原点的距离和观测值。

与一般多项式滤波的不同之处是采用的多项式为一次多项式且滑动窗口D更长，以求得原ΦDP信号各点的线性特征估计值(图1b)，并以每次线性拟合的残差平方和RSS的算式定义该j点ΦDP的线性拟合偏差σj(图1c)，该值越小表明该点符合线性变化特征的程度越好，同时记录下该点线性拟合的斜率k(图1d)。单从线性滑动平滑的结果看，其平滑出的图形与一般性的滑动平均类似，即波动只在幅度上存在不同程度的减小而并没有消失，但滑动线性平滑过程中得到的线性拟合偏差和斜率可以在下面的步骤中为递推求得目标点提供依据。

其中，σj为线性拟合偏差，D为滑动窗口，D=21。RSS为每次线性拟合的残差平方和。

图1 D=21时的线性滑动平滑

(2)对线性拟合偏差进行排序

从经验的角度来说，针对图1(b)中的原始信号，一般先确定图中波动幅度较小且线性特征明显的数据段(如虚线框确定的两段位置)，再依据波动幅度和斜率不为负的原则依次确定其他数据段。对比图1(b)和图1(c)中虚线框位置的对应关系，可以将“先确定波动幅度较小且线性特征明显的数据段”这一目标转化为“先确定线性拟合偏差小的点”。这就需要将上一步骤中得到的线性拟合偏差序列进行升序排序，并求得排序索引，因为希望先确定噪声小且线性趋势明显的数据段，再依据“ΦDP质量控制后斜率不明显为负”的实用目标去确定噪声大且线性趋势不明显的数据段，而不是像一般性的滑动平均、中值滤波、低通滤波等方法那样按照距离库的探测顺序确定质量控制后的最终计算值，因为上述这些按照探测顺序计算的方法显然没有一种机制是专门保留正斜率而抑制负斜率的，即无法考虑后未计算的远端距离库是否还有更可信的但是比近端已经确定的数据还要显著低的数据。因此对于反应噪声和线性趋势的线性拟合偏差σj进行排序并按照重新排序的顺序进行逐库递推计算是非常必要的。

(3)依据排序结果递推目标点

为了确定导数不为负的连续线段，给出逐点递推判断的规则为

依据上述步骤求得排序索引，按照式(1)的规则逐点递推判断该点是否应该为最终求得的点，其中Yi代表将^yj升序排序所得的n个排序索引中第i个索引对应的距离库上线性滑动平滑值，Yleft代表Yi左侧第一个已被赋值的数据值，Yright代表Yi右侧第一个已被赋值的数据值，－999和999代表初始的左右边界。ξ代表松弛变量或称损失函数，表征递推过程中容忍波动的程度。Y0和Yn+1分别代表递推起始的左右边界，k代表线性拟合斜率。在每次递推中，Yleft和Yright中至少有一个会改变，以此判断新的Yi是否应该属于最终求得的“导数不为负的连续线段”上的点。表1给出D=21、ξ=1时，前15次递推的过程列表结果。

表1 部分递推过程举例

为了结合递推规则说明参数ξ与D在该方法中的选择依据，下面首先给出D=21时ξ分别为0、0.5、1、10时的递推结果(图2)，可以看出，随着ξ的增大，递推判断后得到的点的数量增多，但是，波动的幅度也随之增大，尤其是当ξ=10时，150距离库附近开始出现明显的导数值正负波动，并出现导数为负值的点，这时递推过程逐渐泛化为滑动线性平滑;而ξ=0递推出的符合条件的点过少，最终以这些点做线性插值可能不具有代表性。因此，对“尽可能抑制波动及负斜率”与“尽可能多的反应数据细节特征”两个目标折中，经过一系列实践，选取ξ=1。但ξ究竟取什么值才具有绝对的普适性，是一个很难论证的问题。不过此处可以设计成按照经验给定一个ξ初值，当极少数计算达不到某些特定要求时，迭代增加1到2此ξ即可。

图2 取不同松弛函数ξ值的递推结果

而对于D的选取依据，总的来说其不宜过小或过大。因为当D过小时，某些云体边缘带有趋势性的异常噪声可能会被错误的识别为线性趋势，从而导致线性递推错误。图3以另一条扫描径向为例，给出了ξ=1，D分别为11和21的情况，可以看到，700至800距离库间存在一处有连续上升趋势的异常噪声，过短的滑动窗口会使该段计算出的线性拟合偏差很小，从而错误的将该段识别为正常的线性趋势，这将影响它前后位置的递归过程，譬如在该例中，700距离库之前的数据点由于线性平滑后的值以及线性拟合斜率均比该处大，按照递推关系将不被保留，从而导致数据断裂(图3a)。当D=21时，滑动窗口足够大到贯穿这个云体边缘的异常噪声区，便不存在上述问题(图3b)。而为了尽可能多的保留数据趋势变化的细节，D也不宜过大。保守起见，选取D=21(3.15 km)以保证这个滑动窗口能贯穿大部分云体边缘的噪声异常区域。

图3 不同滑动窗口D对递推结果的影响

(4)对递推得到的目标点以外的点进行线性插值

将给定ξ时得到的递推结果中的空缺值进行线性插值，即可得到完整连续的ΦDP质量控制结果。图4(a)为D=21，ξ=1时经过最终质量控制得到的ΦDP结果与原信号的比较。由图4(a)可见，质量控制后ΦDP图形基本滤去了所有噪声，并尽可能地反映了原始信号中的线性变化趋势，且没有明显的负斜率波动，实现了预设的“一组导数不为负的连续线段的确定”的目标。图4(b)为将质量控制后的数据按照最初拼合数据段填补空白的方式再展开，以和ZH进行对比。从图中可以看到，在整个距离径向上，当雷达波束经过连续的ZH大值的较强回波区时，ΦDP的增长是非常显著的，这符合ΦDP变化的基本原理。

图4 D=21，ξ=1时的质量控制结果、原始ΦDP和ZH

为了体现上述质量控制结果的优越性，下面分别列举滑动平均、中值滤波、FIR低通滤波3种具有简单参数的质量控制方法进行对比(图5)。从图中可以看到，这些传统的滤波方法对于云区内幅度不大且较为规则的噪声或者零星噪声都起到了一定的滤除效果，但在云体边缘仍然有异常噪声没有被滤除。而参数更加复杂的滤波方法如卡尔曼滤波和小波滤波虽然在此无法通过个例来说明其质量控制的必然结果，但这些传统的信号处理方法对ΦDP质量控制的局限性是可以预见的。

图5 其他质量控制方法得到的结果

2 在衰减订正中的应用与分析

衰减订正是差分传播相移(ΦDP)的主要用途之一。X-band双偏振雷达由于波长较短，雷达波经过强降水区后反射率会出现明显的衰减，因此，需要对雷达探测到的反射率进行衰减订正。从双偏振雷达的原理上看，ΦDP是不因探测粒子尺度的大小而减弱的，因此，许多双偏振雷达反射率衰减订正方法[1，2，9]都是基于ΦDP进行的。这里列举一种双偏振雷达衰减订正的自适应算法[1，14－15]，该方法通过计算每个距离库上的衰减量AH并积分，再与观测值Z'h相加得到衰减订正后的反射率 Zh，10lg[Zh(r)] =10lg[Z'h(r)] +，衰减量 AH为

其中计算参数I

式中，r0、r1是扫描径向上雨区的起止距离库，ΔΦDP是r1、r0位置上的ΦDP之差。所谓“自适应”算法指的是通过枚举参数α和b可能的取值，得到许多组AH(r，α，b)并利用其重构出 ΦcalDP(r，α，b):

再通过雨区径向范围内 ΦcalDP(r，α，b)与 ΦDP(r)的差值最小作为约束条件确定最优的参数α和b。

以第一节中ΦDP质量控制个例的扫描径向数据为例，参考毕永恒等[15]同样是针对X波段双偏振雷达的衰减订正研究，选取α取值范围0.13～0.35，b值固定为0.8，对于第242至357距离库上的雨区，采用上述衰减订正方法对反射率ZH进行衰减订正，得到图6结果。从图中可以看到，在300距离库以后ΦDP增长明显的部分，ZH在被订正以后也明显增大了。

图6 采用经过质量控制后ΦDP进行雷达回波强度的衰减订正前后对比

为了进一步说明新的ΦDP质量控制方法在衰减订正中的优势，下面列举一个典型的在ΦDP质量控制效果并不理想的情况下衰减订正失败的例子。图7为采用传统的滑动平均滤波方法(这里取13点滑动平均)滤波后的ΦDP代入与图6相同的自适应算法得到的衰减订正结果。由图7可见，订正后的ZH反而减小了，这是因为以滑动平均做ΦDP质量控制并不能有效地滤除ΦDP在云体边缘异常振幅的噪声，因而在计算242到357距离库雨区内的 ΦDP变化时，ΔΦDP的符号为负，这将使式4不能正确地反映衰减量，从而完全影响了该算法的衰减订正结果。针对这一情况，虽然可以通过改变计算雨区范围的大小减小这种影响，但只要质量控制后的ΦDP仍然存在较大量级的噪声或者较多负斜率的情况，这种问题就无法完全避免。

图7 采用滑动平均后ΦDP根据采用自适应算法进行衰减订正

同时还应注意到，其他与ΦDP有关的衰减订正方法，比如以 KDP计算 AH的方法[1，9]、利用 ΦDP与总衰减量成线性关系的方法[2]等，对ΦDP的质量控制结果是同样敏感的。因此，能否进行有效的ΦDP质量控制是双偏振气象雷达衰减订正结果好坏的重要因素。

另外，文中所介绍的线性拟合-递推法在其他观测个例中的ΦDP质量控制与ZH衰减订正工作中也取得了较好的结果，在这里就不再全面介绍了。

3 小结与讨论

提出了一种双参数的线性拟合-递推法来解决双偏振雷达差分传播相移ΦDP的质量控制问题，并通过多个扫描径向上的ΦDP质量控制应用说明该新算法的有效性，结论如下:

(1)提出了ΦDP质量控制的线性拟合-递推法。该方法通过滑动线性平滑得到每个距离库上ΦDP的线性拟合偏差与线性拟合斜率，并根据ΦDP斜率非负性的质量控制目标构造递推关系式，依据线性拟合偏差的排序索引逐点递推，再经过线性插值，最终得到斜率基本非负、可滤除所有噪声且保留有价值变化趋势的ΦDP质量控制结果。

(2)多个实例应用表明，该线性拟合-递推法均取得了较好的ΦDP质量控制结果，不同量级幅度的噪声基本得到了有效的滤除，基本抑制了ΦDP的负斜率，且体现了ΦDP增长的部分与ZH连续大值区的对应关系，说明该方法物理意义明确，而且操作简单。但对于该方法以及两个参数D和ξ的普适性，应在今后的工作中开展更多的实践与讨论。

[1] Park S G，V N Bringi，V Chandrasekar，et al.Correction of Radar Reflectivity for Rain Attenuation at X Band.Part I:Theoretical and Empirical Basis[J].J.Atmos.Oceanic Techno.，2005，22:1621－1632.

[2] 何宇翔，吕达仁，肖辉，等.X波段双线极化雷达反射率的衰减订正[J].大气科学，2009，33(5):1027－1037.

[3] 何宇翔，吕达仁，肖辉.X波段双线极化雷达差分反射率的衰减订正[J].高原气象，2009，28(3):607－616.

[4] 刘黎平，葛润生，张沛源.双线偏振多普勒天气雷达遥测降水强度和液态含水量的方法和精度研究[J].大气科学，2002，26(5):709－719.

[5] 郭凤霞，马学谦，王涛，等.基于X波段双线偏振天气雷达的雷暴云粒子识别[J].气象学报，2014，72(6):1－15.

[6] 魏庆，胡志群，刘黎平.双偏振雷达差分传播相移的五种滤波方法对比分析[J]，成都信息工程学院学报，2014，29(6):596－602.

[7] 杜牧云，刘黎平，胡志群，等.双线偏振雷达差分传播相移的质量控制[J]，应用气象学报，2012，23(6):710－720.

[8] Hu Z Q，Liu L P.Applications of wavelet analysis in differential propagation phase shift data denoising[J].Adv.Atmos.Sci.，2014，31(4):824－834.

[9] 杜牧云，刘黎平，胡志群，等.双线偏振雷达差分传播相移的小波滤波初探[J].暴雨灾害，2012，31(3):248－254.

[10] Hubbert J，Bringi V N.An iterative filtering technique for the analysis of copolar differential phase and dual-frequency radar measurements[J].J.Atmos.Oceanic Technol.，1995，12(3):643－648.

[11] 陈程，肖辉，冯亮，等.北京地区强冰雹风暴的双偏振特征观测分析[J].安徽农业科学，2014，42(22):7511－7515，7516.

[12] 徐士良，马尔妮.常用算法程序集(C/C++描述)(第五版)[M].北京:清华大学出版社，2013:591.

[13] 贾俊平.统计学(第四版)[M].北京:中国人民大学出版社，2011:170－173.

[14] Bringi V N，Keenan T D，Chandrasekar V.Correcting C-band radar reflectivity and differential reflectivity data for rain attenuation:a self-consistent method with constraints[J].IEEE Transactions on Geosience and Remote Sensing，2001，39(9):1906－1915.

[15] 毕永恒，刘锦丽，段树，等.X波段双线偏振气象雷达反射率的衰减订正[J].大气科学，2012，36(3):495－506.