==========
沈阳雷达质控方案
==========

.. toctree::
   :maxdepth: 2
   :caption: Contents:

因为雷达在运行过程基本不保存IQ信号，该方案主要在基数据（即Z，V，和W）的基础上开展雷达资料的质量控制工作。针对国内雷达资料和辽宁雷达的数据特征，质量控制主要涉及包括以下几个方面的工作：地物杂波剔除、异常传播杂波的识别、晴空回波的识别与剔除、径向干扰回识别与剔除、径向速度退模糊，从基数据角度对雷达信号质量进行控制，利用杂波的统计特性将其从有效气象信号中剔除。

地物杂波处理
==========================

利用移动目标检测（MTD，Sawyer等于1985提出）的方案对地物杂波进行抑制并清除，从而得到仅包含降水信号的有效数据，在此方案中，将同时考虑地物杂波的统计特性以及对径向速度（V）的不敏感性。MTD为一种抑制低速杂波的数字运动目标指示器雷达系统，其利用新的技术识别径向速度小于预定值的目标物雷达回波。该方法通过计算作为相干相位检测器输出信号幅度和径向目标速度的函数的滤波器输出，并将这些值存储在一个可由相干相位输出信号幅度来处理的只读存储器中，然后可以进行比较以确定来自滤波器的哪个信号是来自径向速度等于或小于径向目标速度的目标，这样的滤波器输出信号就可以被拒绝。MTD流程如下：

.. figure:: images/MTDflow.png
    :height: 220px
    :align: center
    
    MTD流程图。


径向干扰杂波的质量控制
==========================

综合雷达径向干扰造成的回波污染特征，充分利用径向干扰信号的径向延伸特征，径向速度特征和反射率因子的垂直连续性特征，构建相应的启发式处理方案，一方面可以从回波单元象素的孤立性上考虑，对其进行处理滤波，这对于雷达受外界同频雷达出现的螺旋带状回波和麻点干扰相当有效；另一方面，需要通过其径向的连续性和方位上的不连续性的检查来识别和剔除：

第一、对雷达资料进行预处理，即对取样的雷达资料进行排序。由于CINRAD/SA雷达是在一系列固定的仰角上扫描360°进行取样的，因此，获取的资料中每个PPI扫描面的径向个数不固定，每次取样的方位角度也不是固定对应的。为此，我们需要最低两层的雷达资料进行径向处理，按照严格的1°的间隔顺次排列这些资料，以满足每个PPI的资料都是360个径向，同一方位回波强度、径向速度和速度谱宽资料的对应，以及不同层次雷达资料方位上的对应。

第二、需要对相对孤立的干扰回波进行滤波。这些回波大多以孤立点或细线的方式出现，使用的方法是在PPI 扫描面上, 移动一个 NxN (5×5个点) 的窗口, 如果该窗口内中心点周围的有效值的点少于某一阈值（缺省值为75%）时, 则将其视作噪声回波，予以剔除。 

第三、检查雷达回波在径向上有效回波的连续性，对每个方位上的有效回波的总库数进行统计，找出最大的有效回波总库数的方位 :math:`N`，并记下最大值 :math:`R_{MAX}` 和最小值 :math:`R_{MIN}`。

    .. math::     
        R_i=\sum_{j=1}^{N_R}Z_j   &\\
        Z_j=
        \begin{cases} 
        1  &Z_{i,j}>Spval& \\
        0  &Z_{i,j}\leq Spval
        \end{cases}        
        
 其中， :math:`R` 为计算的每个方位的有效回波库数的总值；:math:`i` 为方位，最大方位为 :math:`N_A`，值为360（即 :math:`1\leq i\leq N_A`）； :math:`N_R` 为距离库的个数，对于SA雷达的回波强度来说，为460个库； :math:`Spval` 为无效回波值，为基数据中读入的无效回波的标示值； :math:`Z` 表示为回波强度值在径向上有效回波的连续性。找出最多有效库数的方位。        
        
第四、对找到的最多有效库数的方位进行判定，确定其是否为干扰回波。先通过几个简单的判断值来判定。   
      
    .. math::     
        R_D=\frac{R_{MAX}-R_{MIN}}{R_{MIN}} \times 100\% 

 其中， :math:`R_D` 是根据最大有效回波的库数与最小有效库数计算的一个判断值。当 :math:`R_D` 的值小于10%时，则认为不是干扰回波， :math:`R_D` 是根据经验取值，一般取10%～20%。并且同时当满足 :math:`R_{MIN}` 大于10（一般可取10～20）， :math:`R_{MAX}` 大于300（一般根据该雷达的库数决定）时，判定找的该径向为出现干扰回波的径向。

第五、检查最大有效径向上的连续性。即对找到的最大有效径向 :math:`N` 进行方位上的连续性检查，找出其周围方位上是否也有干扰回波。采取方法是依靠该方位与周围方位( :math:`N\pm 10` )有效回波总库数的差值来判定，给定阈值(一般选取20～40)，小于该阈值的方位也会标记为干扰回波方位。若大于阈值，则只有该方位有干扰回波。这样，即可找到干扰回波的起始方位 :math:`N1` 和终止方位 :math:`N2` 。


第六、将找出来的干扰回波方位块做为一个整体，进行剔除处理。分别对 :math:`N1-1` 和 :math:`N2+1` 的方位， :math:`N1-2` 和 :math:`N2+2` 的方位进行检查，当对于 :math:`N1-1` 和 :math:`N2+1` 两个方位中的同一库中都没有有效回波时，则认为从 :math:`N1\sim N2`的方位中相同的的库中是没有有效回波的；同理，检查 :math:`N1-2` 和 :math:`N2+2` 的方位同一库中是否同时没有有效回波，若无，则认为 :math:`N1\sim N2` 方位的该库中是没有有效回波的。

经过以上6个步骤的处理，就可以把找到的连续成块状分布的干扰回波标记并剔除。

能流泄露回波的质量控制
==========================

能流杂波目前可以确定与雷达硬件环境紧密联系在一起，主要是雷达接收支路不正常，或测试信号灌入，也有可能来自外界的稳定性信号漏入，导致雷达回波大面积污染的一种非气象杂波现象。能流泄漏严重时，被污染的回波会充斥雷达探测范围的全部360个方位和460 km范围的全部距离库。这类回波具有以下特点：（1）出现具有非常强的随机性，仰角和方位角均不确定；（2）很强的虚假回波；（3）径向速度也可能被污染，所以径向速度不可参考。（4）沿径向方向，观测数据呈递增趋势但变化范围很小。（5）沿方位方向，观测数据变化不大。这类回波的处理相对较为复杂，目前可以采用基于模糊逻辑的方法处理这类回波：

针对能流泄露时，雷达回波布满整个PPI平面，相比于其他各种回波（降水、晴空、其它），回波强度有效探测像素点数最多，采用回波强度有效探测比率，计算回波强度有效探测总像素点数占整个PPI平面像素点的比率：

    .. math:: 
        R_z=\frac{\sum_{i=1}^{Na_z}\sum_{j=1}^{Nr_z}M_z}{Na_z\times Nr_z}  &\\
        M_z=
        \begin{cases} 
        1  &Z_{i,j}\neq Spval& \\
        0  &Z_{i,j}\leq Spval
        \end{cases} 

上式中， :math:`Z_{i,j}` 为回波点 :math:`(i,j)` 强度值； :math:`Spval` 为无效探测值； :math:`Na_z` ， :math:`Nr_z` 为回波强度的有效探测范围的方位数和距离库数。

针对径向速度的有效探测会得到正常的观测值和距离折叠的污染值RF，当速度资料被能流泄露污染时，污染值会标示为RF，计算RF值像素点数占PPI上径向速度有效探测像素点数的比例：

    .. math:: 
        R_{rf}=\frac{\sum_{i=1}^{Na_v}\sum_{j=1}^{Nr_v}M_{rf}}{\sum_{i=1}^{Na_v}\sum_{j=1}^{Nr_v}M_{v}}  &\\
        M_{rf}=
        \begin{cases} 
        1  &V_{i,j}=R_f& \\
        0  &V_{i,j}\neq R_f &\\
        \end{cases} &\\
        M_v=
        \begin{cases} 
        1  &V_{i,j}= Spval& \\
        0  &V_{i,j}\neq Spval
        \end{cases} 

上式中， :math:`V_{i,j}`为 :math:`(i,j)` 回波点的径向速度值； :math:`R_f` 为距离折叠污染值； :math:`Na_v` ， :math:`Nr_v` 为径向速度的有效探测范围的方位数和距离库数。

针对多普勒天气雷达进行观测时，获取的回波强度资料和径向速度资料在同一位置是一一对应的。即在有回波强度值对应的位置也对应有观测到的速度值，特别是对于观测到的正常回波资料，有效回波的速度和强度的不匹配比率将会很小。

    .. math:: 
        R_{vz}=\frac{\sum_{i=1}^{Na_v}\sum_{j=1}^{Nr_v}M_{v}}{\sum_{i=1}^{Na_v}\sum_{j=1}^{Nr_v}M_{z}} 

该式中的 :math:`M_z` 和 :math:`M_v` 和之前所指相同。

计算每个距离圈上统计回波强度沿径向变化小的像素点，得到一圈中变化小的像素点占每圈总数的比率：

    .. math:: 
        Ra=\frac{\sum_{i=1}^{Na_z}M_{Zv}}{\sum_{i=1}^{Na_z}M_{Na}}  &\\
        M_{Zv}=
        \begin{cases} 
        1  &Zv_{i,j}\leq 1& \\
        0  &Zv_{i,j}> 1 &\\
        \end{cases} &\\
        M_{Na}=
        \begin{cases} 
        1  &Zv_{i,j} \neq Spval& \\
        0  &Zv_{i,j} =  Spval
        \end{cases} 


    .. math:: 
        R_{Nrz}=\frac{\sum_{i=1}^{Nr_z}M_{Ra}}{\sum_{i=1}^{Nr_z}M_{a}}  &\\
        M_{Ra}=
        \begin{cases} 
        1  &Ra\geq 60\%& \\
        0  &Ra < 60\% &\\
        \end{cases} &\\
        M_{a}=
        \begin{cases} 
        1  &Ra \neq Spval& \\
        0  &Ra  =  Spval
        \end{cases} 

计算每个方位统计每根径向上沿方位的回波强度库间变化比率 :math:`R_{Naz}` ：

    .. math:: 
        Rr=\frac{\sum_{i=1}^{Nr_z}M_{Zv}}{\sum_{i=1}^{Nr_z}M_{Nr}}  &\\
        M_{Nr}=
        \begin{cases} 
        1  &Zv_{i,j} \neq Spval& \\
        0  &Zv_{i,j} =  Spval
        \end{cases} 

    .. math:: 
        R_{Naz}=\frac{\sum_{i=1}^{Na_z}M_{Rr}}{\sum_{i=1}^{Na_z}M_r}  &\\
        M_{Rr}=
        \begin{cases} 
        1  &Rr\geq 50\%& \\
        0  &Rr < 50\% &\\
        \end{cases} &\\
        M_{r}=
        \begin{cases} 
        1  &Rr \neq Spval& \\
        0  &Rr  =  Spval
        \end{cases} 

在这些参量的基础上，采用基于模糊逻辑方法，根据如下图所示的有效回波强度比率，RF值比率，有效回波匹配比率，以及沿径向回波强度库间变化和沿方位回波强度库间变化的隶属函数，识别和剔除能流泄露杂波的影响。

.. figure:: images/MBF.png
    :height: 550px
    :align: center
    
    识别能流泄露杂波的5个参量的隶属函数图。


雷达反射率因子与地面雨滴谱的比对与标校
==========================

将雷达的所有部件的输出作为一个整体进行校准，首先需要基于T-matrix散射矩阵方法和本地雨滴谱数据计算雷达相关参量，将计算结果与雷达实测数据进行大规模比较，从而得到雷达测量的数据在不同降水类型下的系统测量误差，并将该系统误差应用到雷达的实测数据中。

T-Matrix方法是由Waterman 于1965年正式介绍的，在随后相关领域研究和文章中被广泛应用和推广。T-Matrix法是计算非球形粒子散射、解释非球形粒子散射特性的有效方法之一，已成为处理电磁、声波和弹性波散射等广阔科学领域的核心技术。在T-Matrix法中，入射和散射的电场在一系列合适的(矢量)球谐函数中扩展，并且通过过渡矩阵(T-Matrix)来建立各个展开系数的列之间的关系。入射和散射场矢量球谐函数展开形式如下(Bo Peterson et al.，1974; 郭丽君等，2014)：

    .. math:: 
      \overline{E^1(\bar r)} = \sum_{v=1}^\infty D_v[a_vM_v^1(k\bar r) + b_vN_v^1(k\bar r)]  &\\
      \overline{E^s(\bar r)} = \sum_{v=1}^\infty D_v[f_vM_v^3(k\bar r) + g_vN_v^3(k\bar r)] 

其中 :math:`Dv` 是标准化系数，:math:`M_v^1` 、:math:`N_v^1` 、 :math:`M_v^3` 、 :math:`N_v^3` 是矢量球谐函数，:math:`a_v` ， :math:`b_v` ， :math:`f_v` ， :math:`g_v` 为矢量球谐函数的展开系数，则散射场展开系数可由入射场展开系数通过下式确定：

    .. math:: 
      \left[ \begin{array}{c}
      f_v \\
      g_v  \end{array} \right]=[T] \cdot \left[ \begin{array}{c}
      a_v \\
      b_v  \end{array} \right]

其中 :math:`T` 为转换矩阵，通过在散射体表面的二维数值积分求得（Barber et al.，1975），从而得到散射场。 :math:`T` 矩阵算法的核心是计算矩阵 :math:`T` ，求得 :math:`T` 就能求得散射场，进而求得任意方向上的散射截面 :math:`\sigma (\theta_s)` ，

    .. math:: 
      \sigma (\theta_s)=4\pi |\overline F|^2 
      
其中 :math:`\overline F` 为散射场远场矢量振幅， :math:`\theta_s` 为散射角，当 :math:`\theta_s = \pi` 时，得到后向散射截面。

理想情况下，雷达反射率因子与雨滴谱模拟反演的反射率因子应具有非常高的一致性。在雷达标定不完善的情况下，二者之间会存在一个明显的数据偏差，只需要在雷达资料处理之后，将这个偏差叠加在雷达反射率因子上即可得到质量控之后的数据资料。

具体而言，该研究方案首先利用激光雨滴谱观测资料反演得得一个反射率因子时间序列结果，如下图中黑色虚线所示；同时与雨滴谱站点匹配的雷达观测得到的反射率因子用红色虚线标识。由此对比可以看出雷达测得的原始反射率因子比雨滴谱模拟的反射率因子时间序列明显偏弱，但订正后的结果与雨滴谱模拟结果显著增强。

.. figure:: images/DSD_compare1.png
    :height: 250px
    :align: center
    
    雷达测得的反射率因子和雨滴谱模拟的反射率因子时间序列（原始观测偏弱）。

.. figure:: images/DSD_compare2.png
    :height: 250px
    :align: center
    
    雷达测得的反射率因子和雨滴谱模拟的反射率因子（质量控制和订正后）。

在业务化过程中，可以用雷达测值和雨滴谱反演结果的平均差异 :math:`\triangle Z` 进行定量评估：

    .. math:: 
      \triangle Z =\frac{1}{n}\sum_{i=1}^{n}(Z_H^R-Z_H^{DSD}) 

其中，:math:`Z_H^R` 是雷达观测的反射率因子，:math:`Z_H^{DSD}` 是地面雨滴谱反演得反射率因子。下图中原始的雷达观测资料的平均偏差 :math:`\triangle Z` 在-4.875dBZ，经过订正后平均偏差 :math:`\triangle Z` 减小为-0.134dBZ。

当然，雷达测得的反射率因子也可能因为雷达系统误差出现偏大的情况。 如下图所示，雷达观测资料的平均偏差 :math:`\triangle Z` 在+3.533dBZ；经过订正和质量控制，雷达反射率因子的偏差得到有效抑制，并使雷达观测资料的平均偏差 :math:`\triangle Z` 减小为-0.103dBZ。
    
.. figure:: images/DSD_compare3.png
    :height: 250px
    :align: center
    
    雷达测得的反射率因子和雨滴谱模拟的反射率因子时间序列（原始观测偏强）。    

.. figure:: images/DSD_compare4.png
    :height: 250px
    :align: center
    
    雷达测得的反射率因子和雨滴谱模拟的反射率因子时间序列（质量控制和订正后）。
     
值得一提的是，通过雷达观测（包括质控前后数据）与地面雨滴谱的比对与标校分析，我们不仅能够实现对单偏振雷达反射率因子的评估，标定和误差修正，也能多双偏振雷达参量进行定量评估，质控分析。这讲为雷达双偏振升级后的数据质控提供保障。因此，该项目推进过程中，我们也对双偏振参量进行了研究分析。下图演示了雨滴谱模拟的雷达反射率因子（ZH）需要满足与差分传播相移率（KDP）之间良好的非线性关系（下图中的拟合曲线）。如图a和b所示的两个散点对比结果，椭圆内的订正前的ZH明显偏离了理论上KDP-ZH关系，这些偏离的反射率因子数据就是需要订正的数据。图c和d描述了订正后的ZH观测数据与KDP之间的散点关系，结果表明订正后的反射率因子数据明显更加趋近于理论的模拟结果。
    
.. figure:: images/DSD_compare5.png
    :height: 600px
    :align: center
    
    雷达测得的反射率因子与雨滴谱反演的反射率因子与差分相移率（KDP）的对比结果。   


径向速度退模糊
==========================

径向速度的质量控制方案包括预处理和退模糊两部分，主要计算流程如下图所示。预处理部分主要包括插值和剔除噪音两个步骤。插值的目的是把雷达资料转化为规则格点资料，方便与分析场格点资料进行比较。首先根据Mohr和Vaughan (1979)提出的插值方法对雷达径向风进行插值，从原始极坐标系转化为三维笛卡尔坐标系，插值后资料格点水平分辨率为1 km，垂直分辨率为 500 m。后续的处理步骤，都是在已插值后的笛卡尔坐标系的资料上进行。

径向速度的质量控制主要也是对径向速度进行退模糊，尤其在高强度风的情况下进行退模糊，利用科罗拉多州立大学（CSU）提出的自适应算法寻找速度模糊的距离库，并将该距离库观测到的风场和附近距离库的风场数据进行比较分析，如果该算法确认该径向速度值折叠（模糊），则利用雷达系统最大不模糊速度信息对该观测进行退模糊。 

.. figure:: images/Radial.png
    :align: center
    
    雷达观测资料径向速度的预处理和质量控制流程。    

简要来说，可用 :math:`x(k)` 表示输入信号， :math:`y(k)` 表示输出信号， :math:`d(k)` 表示期望信号， :math:`F` 是输入、输出和参考信号的一个函数，即 :math:`F=F[e(k)]=F[x(k), y(k), d(k)]` ，自适应算法试图使函数 :math:`F` 达到最小化，从而 :math:`y(k)` 与 :math:`d(k)` 近似相等，使得参数 :math:`θ(k)` 收敛到 :math:`\theta _0` ， :math:`\theta _0` 为导致目标函数 :math:`F` 最小化的最优系数构成的集合(Diniz，2014)。

观测误差的确定：对比使用Hollisworth 1986方法、Desroziers 2001方法、Desroziers 2005方法，诊断雷达反射率和径向速度资料的观测误差，在对同化预报效果进行对比检验的基础上，选择合适的观测误差诊断方法进行应用。