1 引言

近年来，X波段双线偏振多普勒天气雷达(简称X波段偏振雷达)作为业务天气雷达的有益补充，可以提供比业务天气雷达更加细微的云、雨滴、冰雹等水成物的微物理过程信息，使得X波段偏振雷达在龙卷、下击暴流、暴雨和冰雹等灾害性天气过程的监测和预警以及人工影响天气作业指挥和效果评估等方面有着巨大的应用潜力(刘黎平等, 1996, 1997；王致君等，2002；陈羿辰等，2017)。北京市人工影响天气办公室、青岛市气象局、上海市气象局等业务单位通过引进X波段偏振雷达，分别在北京奥运会、青岛奥帆赛、上海世博会等重大气象保障工作中发挥了重要作用(马建立等，2008；刘宝君等，2008；吴志根等，2010)。

差分传播相移(ϕ_DP)是X波段偏振雷达非常重要的一个探测参量。若X波段偏振雷达采用双发双收模式，则雷达探测的偏振参量包括反射率、差分反射率、差分传播相移、差分传播相移率以及相关系数(Bringi, et al，2001)。Ryzhkov等(1996)指出，差分传播相移率(K_DP)受雷达波束充塞系数的影响较小，而且不受降水粒子衰减的影响，因此K_DP常被用于定量降水估测和衰减订正。偏振参量差分传播相移率是由差分传播相移计算得到，差分传播相移数据质量的好坏将决定差分传播相移率的准确性，如果反射率、差分反射率采用相位方法进行订正(Park et al, 2005a, 2005b；何宇翔等，2009)，则ϕ_DP数据质量将影响反射率、差分反射率进而影响差分传播相移率这3个偏振参量的数据质量，因此，差分传播相移数据质量控制显得尤为重要。由于X波段偏振雷达波长较短，雷达发射的电磁波经过冰雹等大粒子区域时，米散射引起的后向差分传播相移(δ)会使差分传播相移出现“陡增”现象，后向差分传播相移和差分传播相移混杂在一起是雷达探测的ψ_DP，因此，在利用雷达探测得到的ψ_DP计算差分传播相移率时必然会产生误差，导致利用差分传播相移率进行降水估测或粒子相态识别时准确率降低。因此，在利用差分传播相移计算差分传播相移率时需要消除米散射引起的后向差分传播相移的影响。

对差分传播相移数据质量控制方法已开展了大量的研究，丁青兰等(2003)比较和分析了不同差分传播相移率计算方法的精度；何宇翔等(2009)引入卡尔曼滤波开展差分传播相移的质量控制，并用质量控制后的差分传播相移进行反射率的衰减订正；Hubbert等(1993)通过设计低通滤波器对回波的差分传播相移进行滤波，发现该方法不能有效抑制后向差分传播相移的变动，随后Hubbert等(1995)采用低通FIR滤波器，利用循环迭代的方法，消除后向差分传播相移的影响，但该方法并不能保证电磁波经过雨区后差分传播相移的累积递增属性；毕永恒等(2012)也同样采用该方法，对中国科学院大气物理研究所的X波段偏振雷达，先进行差分传播相移数据预处理，然后采用自适应约束算法对雷达反射率进行衰减订正；杜牧云等(2012)采用小波分析方法，对差分传播相移进行滤波处理，指出小波滤波后，差分传播相移具有更好的连续性和平滑度，随后杜牧云等(2013a, 2013b)对C波段偏振雷达开展了数据质量分析, 指出ϕ_DP经分类处理后拟合的K_DP对降水粒子的特性反映比较符合实际情况；曹俊武等(2011)分析了影响差分传播相移探测精度的因素，并指出对差分传播相移进行3 km的平滑时，可以满足差分传播相移的探测精度要求；肖艳娇等(2012)使用地物杂波抑制、退折叠、初始值调整、滤波等方法对中国气象局武汉暴雨研究所X波段双线偏振雷达的资料做了质量控制，并给出各个径向获取初始相位的方法；郭凤霞等(2014)利用714XDP偏振雷达观测的偏振参量对中国西北地区夏季雷暴云内的水凝物进行了分类；王洪等(2016, 2018)先后开展了S波段偏振雷达偏振参量的模拟研究和超级单体的偏振参量表现特征分析。

上述差分传播相移质量控制方法在消除后向差分传播相移和差分传播相移数据波动方面有一定的效果，但是得出的差分传播相移率均存在负值，或是质量控制后差分传播相移并不能保证是累积递增属性(累积递增属性，即要求在降水区ϕ_DP是累积单调增加的，K_DP>0，而在非降水区，ϕ_DP不增加的区域，ϕ_DP增量为0，K_DP=0，ϕ_DP径向上保持经过雨区最后一个库的值)，理论上与差分传播相移应该具备的累积递增属性不符。

文中将采用数学统筹方法——线性规划(LP)对北京地区组网的4部X波段偏振雷达探测的差分传播相移进行数据质量控制, 有效剔除X波段偏振雷达在米散射条件下产生的后向差分传播相移，以保证差分传播相移率的非负属性和差分传播相移的累积递增属性。

2 线性规划 2.1 线性规划标准型

线性规划是研究线性约束条件下线性目标函数极值问题的数学理论和方法，是合理利用有限资源制订最佳决策的有力工具，广泛用于军事作战、经济分析、经营管理和工程技术等方面。线性规划问题常用形式是标准型，包括以下3个部分：

(1) 一个需要极大化的线性函数，例如:c₁x₁+c₂x₂；

(2) 以下形式的线性约束条件，例如:

(3) 非负变量，例如x₁≥0，x₂≥0。

从标准型可以看出，“线性”要求变量(x₁，x₂)的幂值都是一次的，线性规划问题常用矩阵形式表达，即目标表达式C ^T X，在条件AX ≤b和X ≥0约束条件下取得最小值。其他类型的问题，例如极大化问题、不同形式的约束问题和有负变量的问题，均可以改写成其等价问题的标准型。

2.2 线性规划在差分传播相移质量控制中的应用

假定X波段偏振雷达探测的某条径向数据中有n个距离库，每个距离库上的差分传播相移构成的矩阵向量为ψ =(ψ₁，ψ₂，…，ψ_n)，经过滤波或处理后的差分传播相移构成的矩阵向量为ϕ =(ϕ₁，ϕ₂，…，ϕ_n)。因此，差分传播相移线性规划问题转化为使得矩阵ϕ和ψ的差值最小，即最小，选取中间矩阵矢量为z =(z₁, z₂, …, z_n)，令|ϕ_i-ψ_i|≤z_i，要使得最小，即是在z_i-ϕ_i≥-ψ_i和z_i+ϕ_i≥ψ_i两个限制条件下，使得价值函数最小，令矩阵ψ_c=(Z, ϕ)^T，要使价值函数最小，即使C·ψ_c最小，其中C =(1₁, …, 1_n, 0_n+1, …, 0_2n)，差分传播相移的线性规划问题转化为在条件Aψ_c≥ ψ限制下，矩阵向量C·ψ_c取得最小值，其中，I_n为n×n单位矩阵，ψ_c=(Z, ϕ)^T，T表示矩阵转置。如果没有其他限制条件，当ψ＝ϕ时，价值函数=0；显然这个起不到对差分传播相移处理的效果，为此，需要增加限定函数。由于差分传播相移的物理性质是单调增加的，因此，单位距离内差分传播相移的变化率是正值，即K_DP≥0, 为此引入Savitzky-Golay (SG)二阶多项式差分滤波器(Madden，1978)，其m点的表达式为

(1)

m点差分滤波构成矩阵M_{n-m+1, n}，矩阵K_DP≥0可以表示为M_{n-m+1, n} ϕ ^T≥ Z_n-m+1，Z_n-m+1为0矩阵。上述的线性规划问题可转换为：在限制条件A_AUG ψ_c≥ ψ_AUG下，使C·ψ_c最小，其中，

(2)

假定选取5点差分滤波，C_SG={-0.2, -0.1, 0, 0.1, 0.2}，则有

(3)

Giangrande等(2013)指出，在线性规划时，用了m点的差分滤波，还需用k=(m-1)/2点数的平滑滤波。平滑滤波器形式为

(4)

且要求

(5)

为了要求差分传播相移率的非负条件，则滤波要求

(6)

下面利用北京顺义X波段偏振雷达2017年7月7日21时09分(北京时，下同)2.4°平面位置显示(PPI)，325°方位径向差分传播相移数据，开展5点和25点差分滤波线性规划效果分析(图 1)，图 1中Raw-ϕ_DP表示雷达输出的原始差分传播相移数据，Lp-ϕ_DP-5(25)表示对原始差分传播相移数据开展5点(25点)平滑的线性规划结果，Lp-K_DP-5(25)表示按5点(25点)平滑的差分传播相移线性规划结果计算的差分传播相移率。差分传播相移率是根据相邻两个距离库线性规划后差分传播相移的值按式(7)计算得到的。

(7)

式中，ϕ_DP(r_n)和ϕ_DP(r_n+1)分别为在降水区内相邻距离r_n和r_n+1处测得的双程传播相移。

从图 1可以看出，5点和25点线性规划都能满足K_DP≥0，但是25点线性规划得到差分传播相移更加平滑，差分传播相移率波动性小于5点线性规划。由于雷达距离库长75 m，5点差分滤波相当于0.375 km滑动窗口内保证K_DP≥0，25点差分滤波相当于1.875 km滑动窗口内K_DP≥0。在后续内容中，线性规划选择25点差分滤波器。

3 线性规划数据处理 3.1 数据源和观测设备

线性规划可以有效剔除X波段偏振雷达在米散射条件下产生的后向差分传播相移，为了验证剔除后向差分传播相移的效果，选取北京4部组网X波段偏振雷达2017年7月7日冰雹天气过程的观测数据进行分析。当日夜间，受东移强对流云团的影响，北京地区出现冰雹、大风等强对流天气，延庆、昌平、怀柔、顺义、门头沟、海淀出现冰雹；全市大部分地区出现8级以上短时大风，其中延庆、昌平、怀柔、门头沟、平谷出现10—11级以上短时大风。虽然雷达定量探测可以到150 km，但此次天气过程较强，X波段偏振雷达电磁波衰减严重，故给出的雷达回波(图 2)只显示了75 km范围。4部X波段偏振雷达主要性能指标见表 1。

3.2 非气象回波对线性规划的影响

从图 2可见，房山、顺义和通州雷达的站点周围存在很多低于15 dBz的弱化回波，显然不是降水回波，应归为非气象回波类。为了分析该弱回波对差分传播相移的影响，以顺义雷达为例进行分析，可见强天气过程回波还处于雷达站点的西北方向(方位240°—60°)，未到达雷达站点，站点东南方向(方位60°—240°)基本上是小于15 dBz的弱回波，可以判断为非气象回波，但从顺义雷达差分传播相移数据(图 3)可以看出，东南方向弱回波却导致差分传播相移增大(约从40°增至120°)迅速，这显然不是雷达电磁波经过雨区导致差分传播相移增大。因此，可以判断雷达站点东南方向差分传播相移增大是由非气象回波导致，可能是地物或是生物目标导致，而地物和生物目标通常具有相关系数小的特点。由于非气象回波会导致差分传播相移的增大，且存在波动大的特点(图 4)，0—3 km处，差分反射率(Z_DR)波动范围在-8—0.5 dB，主要表现为负值，有可能是受地物影响；3—17 km处差分反射率主要表现为正值，均值在0附近，波动范围约在-8—6 dB，有可能是生物目标或是晴空回波。差分传播相移在0—17 km处波动范围在0°—360°，而其线性规划是按径向所有距离库进行规划的，因此，非气象回波必然影响线性规划结果，为此需要剔除非气象回波。

由图 4可知，非气象回波差分传播相移和差分反射率波动大，因此，对气象回波和非气象回波差分反射率和差分传播相移的纹理特征开展分析，纹理计算表达式为，x为观测参量，x_i为某距离库观测参量值，x 为相邻N个距离库参量x_i平均值，N表示选取的径向距离库数，此处取N=7，表示选取7个径向距离库计算纹理。图 5是图 3中方位120°和325°纹理变化，图 5a为差分反射率纹理，图 5b为差分传播相移纹理，图中A_z=120表示120°方位非气象目标(地物或生物目标)纹理特征，A_z=325是325°方位气象目标纹理特征，可见气象目标和非气象目标纹理特征分界明显，差分反射率分界线大致为SD(Z_DP)=1，差分传播相移分界线大致为SD(ϕ_DP)=5，因此，可以利用纹理特征开展非气象目标的剔除。

通过上述分析，气象回波和非气象回波在相关系数、差分反射率纹理、差分传播相移纹理特征有明显的差异，因此可以考虑用相关系数、差分反射率纹理、差分传播相移纹理剔除非气象回波。文中后续剔除非气象回波综合利用了相关系数、差分反射率纹理、差分传播相移纹理特征，取ρ_HV < 0.9，SD(Z_DP)>1，SD(ϕ_DP)>5作为剔除地物约束条件。

3.3 数据平滑对线性规划结果的影响分析

利用图 2中房山X波段偏振雷达324°径向差分传播相移数据进行分析，采用3.2节中相关系数、差分反射率纹理、差分传播相移纹理信息剔除非气象回波，然后对该径向原始差分传播相移数据进行线性规划，以及对原始差分传播相移采取3、5、7点平滑，看平滑后线性规划的效果，图 6中ϕ_DP-3表示差分传播相移进行3点平滑结果，ϕ_DP-3-Lp表示3点平滑后差分传播相移线性规划结果，其他含义类推。图 7为原始差分传播相移和3、5、7点平滑后的差分传播相移线性规划结果比较，从不同平滑点数看，除了在线性规划起始的41—42 km距离内，原始差分传播相移线性规划与平滑的线性规划结果有一点差异外，42 km以外的原始差分传播相移线性规划结果与3、5、7点差分传播相移平滑后线性规划结果具有非常好的一致性，线性规划结果均方根误差接近0。但是平滑后，后向差分传播相移要比没有平滑前特征明显，图 8给出了利用差分传播相移线性规划结果对反射率和差分反射率订正的结果，订正方法采用何宇翔等(2009)的相位订正方法，从订正的效果看，Z_H和Z_DR得到了明显改善；随着差分传播相移的增大，相关系数在反射率(未订正前)的最大值63.7 dBz(对应距离46.875 km，距离库=625)附近从约0.97开始逐渐减小，说明电磁波进入了冰雹区域，在约48 km处相关系数达到最小值(0.82)，然后逐渐增大至0.90左右，并在此波动；在相关系数逐渐减小的过程中，后向差分传播相移增大明显，雷达发射的电磁波经过冰雹区域后，发生了米散射，在48—55 km径向距离范围内，差分传播相移与线性规划后的差分传播相移并不重合，结合相关系数基本小于0.9，反射率在60—40 dBz变化，可以判断在48—55 km径向距离内，主要应该是冰雹区域或是雨夹雹集中区域。为了分析线性规划前、后后向差分传播相移对差分传播相移率计算结果的影响，将房山X波段偏振雷达2.4° PPI 324°方位原始差分传播相移7点平滑结果和7点平滑后线性规划结果，采用相邻距离库分别计算差分传播相移率(图 9)，图 9a中ϕ_DP-7表示原始差分传播相移7点平滑，K_DP-7-LSF表示用Wang等(2009)最小二乘法拟合(LSF)的差分传播相移率，图 9b中ϕ_DP-7-Lp表示原始差分传播相移7点平滑后线性规划结果，K_DP-7-Lp表示原始差分传播相移7点平滑后，根据线性规划结果计算的差分传播相移率。线性规划方法得到的差分传播相移率未出现负值，且波动小，最小值0°/km，最大值5.2°/km，而常用的最小二乘法拟合方法得到的差分传播相移率出现负值，且波动范围大，最小值-6.3°/km，最大值12.1°/km；而通常K_DP < 8°/km(Giangrande, et al，2013)，显然线性规划明显优于最小二乘法拟合。

3.4 线性规划效果

利用上述剔除非气象回波方法对北京房山、顺义、昌平、通州X波段偏振雷达2017年7月7日观测的2.4° PPI差分传播相移开展线性规划数据质控(图 10)，可以看出，质量控制前各个雷达各个径向方位上差分传播相移均有起伏，并不满足累积递增属性，质量控制后各个雷达差分传播相移均为单调增大，且利用质量控制后的差分传播相移数据计算得到的差分传播相移率均为正值(图 11)。

图 11中线性规划前的差分传播相移率是雷达基数据中给出的值，K_DP < 0的值对应色标为黑色。线性规划后K_DP≥0，与理论上X波段偏振雷达发射的电磁波经过雨区差分传播相移是单调增大的物理属性一致。

对X波段偏振雷达而言，冰雹区域容易产生后向差分传播相移效应，且冰雹区域对应的相关系数偏小，图 12a—h分别给出了房山、顺义、通州、密云雷达2.4° PPI剔除非气象回波后对应的相关系数以及利用模糊逻辑算法(曹俊武等，2005)识别出的冰雹结果，可见在识别出冰雹的区域(红色区域)，对应的相关系数较小，基本上ρ_HV < 0.9，衰减订正后的反射率最强值的地区不一定和冰雹区对应(房山雷达表现出最强回波与冰雹区域不完全一一对应)，这点与刘黎平等(1993)结论一致。为了验证线性规划对冰雹区域的后向差分传播相移抑制作用，图 12e₁—h₁为原始差分传播相移经过7点平滑后得到的ϕ_DP′在识别的冰雹区域与线性规划后ϕ_LP的差值(Δϕ_DP(≈δ)=ϕ_DP′-ϕ_LP)。可以看出，冰雹区域δ值最大可超过15°，如果按相邻两个距离库计算差分传播相移率，δ=15°，K_DP=93.3°/km，而通常K_DP < 8°/km，因此，δ会对差分传播相移率的计算产生较大误差，影响利用差分传播相移率估测降水的精度。

3.5 线性规划并行计算

由3.4节内容可以看出，差分传播相移线性规划可以起到非常好的差分传播相移滤波作用，以及保证差分传播相移率的非负性质，但是由于线性规划是对雷达各个体扫数据各条径向所有数据的统一规划，因此，各条径向数据的多少将决定该条径向数据的计算量。北京参与组网的5部X波段偏振雷达(1部2017年7月仍在建设中)，最远探测距离230 km，对应距离库数3066，由于X波段雷达在有天气过程时衰减较强，假定取1000库计算(对应75 km)，一共9层PPI，每层360条径向，按2.2节介绍内容，要完成一个体扫数据差分传播相移线性规划，若采用5点差分滤波器，需要计算9×360个大小为2995×1000矩阵数据，计算量非常庞大，当北京5部X波段偏振雷达数据组网时，每部雷达都要进行线性规划，这个计算量就更加庞大。因此，需要考虑提高差分传播相移线性规划算法的计算速度。文中线性规划算法是在vs2010平台中实现的，在程序中引入OpenMP并行命令。OpenMP是一种用于共享内存并行系统的多线程程序设计方案，支持的编程语言包括C、C++和Fortran。采用并行计算后，线性规划方法运算速度大幅度提高。

4 结论和讨论

利用线性规划方法对差分传播相移进行质量控制，可以有效保证差分传播相移是一条累积递增光滑的曲线，与物理期望的意义一致。不论是在降水区还是非降水区域，线性规划对差分传播相移质量控制都是适用的，原因是线性约束条件(式(3))保证了在降水区K_DP>0，即差分传播相移的递增属性，而在非降水区，差分传播相移不增大的区域K_DP=0，即差分传播相移增量为0，因此式(3)线性约束条件非负性保证了线性规划方法在降水和非降水区域都适用，只是在利用线性规划后的差分传播相移订正反射率和差分反射率时应注意是在雨区进行，因为反射率和差分反射率衰减订正系数是在雨区得到的，并不是在非雨区得到的，非雨区的订正系数获取非常困难，目前尚无有效方法。

采用线性规划方法，对北京参与组网的其中4部X波段偏振雷达探测的差分传播相移进行质量控制，结果表明，线性规划可以有效降低因冰雹等大粒子导致的后向差分传播相移对差分传播相移的影响，利用线性规划后的差分传播相移数据对降水区域反射率、差分反射率进行衰减订正，订正后得到的偏振参量利用模糊逻辑方法识别的冰雹区域中后向差分传播相移效应明显，线性规划后可以有效降低后向差分传播相移对差分传播相移率计算的影响。

线性规划方法用于差分传播相移质量控制，计算量大，在实际应用过程中需要考虑用并行计算以提高计算速度，尤其是北京5部组网的X波段偏振雷达都采用线性规划方法进行差分传播相移质量控制时，并行计算的应用就显得非常必要。

下一步工作，增加线性规划约束条件，利用线性规划得到的差分传播相移率，开展降水估测，利用降水估测结果，与其他计算差分传播相移率的方法进行比较，检验线性规划得到的差分传播相移率在降水估测中的效果。