基于EMD方法的X波段双偏振雷达衰减订正

李 海， 冯兴寰, 孟凡旺

(1. 中国民航大学天津市智能信号与图像处理重点实验室， 天津 300300；2. 中国航空工业集团公司雷华电子技术研究所， 江苏无锡 210403)

0 引言

天气雷达是监测和预警突发灾害性天气最有效的手段，常规的天气雷达通过反射率与降雨率的经验关系实现对降雨的探测和预报，然而此类雷达无法满足对大范围降水测量准确度的要求，双偏振雷达能提高天气雷达的探测能力，包括降水估算能力和粒子相态识别能力。与C波段和S波段相比，X波段双偏振雷达对目标的定位更准确，并且具有天线尺寸小、易于移动等优点。当双偏振雷达探测气象环境时，回波信号会被路径中降雨区域吸收，使接收端的反射率产生衰减，特别是X波段(中心波长3 cm)，由于其波长较短，反射率衰减问题更为严重，使实测的反射率与真实的反射率之间存在差异，导致双偏振雷达对降雨估计不准确，因此必须进行衰减订正。

Zhang等根据雷达气象方程和-关系，提出了雷达径向逐库算法进行衰减订正，该算法减小了过度订正问题，但没有解决衰减订正中不稳定的-关系。随着双偏振雷达的发展，Rahimi提出差分传播相移与衰减率之间可以近似为线性关系，在使用差分传播相移进行衰减订正时需消除后向散射以及环境噪声的影响，得到准确的差分传播相移。常见的平滑滤波，中值滤波等方法虽然能够改善差分传播相移数据的平滑性，但失去了对原有数据变化趋势的反映。Hubbert等提出采用FIR和IIR低通滤波器估计差分传播相移，该方法可得到差分传播相移的平均走势，但随着距离库的增加，不能有效抑制差分传播相移波动。何宇翔等提出的卡尔曼滤波方法能够对差分传播相移进行较精确的估计，但这类方法对于距离门之间数据的关联性较强，导致运算速度较慢。杜牧云等提出了小波分析法，通过小波分析估计得到的差分传播相移具有良好的平滑度，并减少了差分传播相移率的负值，但是此方法需要选择合适的小波基，影响运行速度。Bringi等提出了迭代滤波方法，该方法能够达到剔除干扰的目的，但是迭代次数难以确定，数据处理时间较长。文献[11]利用MCMC方法对差分传播相移进行处理，该方法能对差分传播相移进行滤波处理，但MCMC方法存在一个显著问题就是需要计算接受率，导致计算量大，并且由于接受率的原因导致算法收敛时间变长。目前，上述方法均是将降水区域视为连续的区域进行衰减订正处理，且将天气雷达不同偏振参量的关系理想化为线性关系，当雷达测量值出现不连续或缺测点时将影响整个降雨区域的衰减订正效果，因而这些方法具有一定的局限性。

本文提出了一种基于EMD方法的X波段双偏振雷达衰减订正，该方法首先利用EMD方法对总差分传播相移进行自适应分解，获得由高频到低频分布的多个IMF，其次通过皮尔逊相关系数准则对IMF进行筛选，将有用的IMF进行重构获得差分传播相移，且将差分传播相移采用最小二乘法拟合差分传播相移率，然后将求得的差分传播相移与差分传播相移率采用自适应约束方法进行反射率衰减订正，最后对所提方法进行了仿真实验。实验结果表明所提方法可以降低观测数据带来的误差，并极大程度上保证差分传播相移的递增性，从而实现更加准确的衰减订正。

1 基于EMD的衰减订正方法

本文使用EMD方法进行衰减订正，首先采用EMD方法估计差分传播相移，之后再对差分传播相移采用最小二乘法估计差分传播相移率，最后对求得的差分传播相移与差分传播相移率采用自适应约束方法进行反射率衰减订正。下面将分别进行详细描述。

1.1 基于EMD的差分传播相移估计

假设降水区域包含个距离库，将第个距离库的水平和垂直偏振波相位分别记为()、()，则该距离库的差分传播相移()为

()=()-()

(1)

它表示水平、垂直偏振波传播到第个距离库之后散射回来的信号相位差。实际雷达观测到的总差分传播相移()由差分传播相移与后向散射差分相移构成，表示为

()=()+()

(2)

式中，()表示该距离库后向散射差分相移，是散射过程中降雨粒子本身引起的相位差。从频率上看，()属于高频噪声。对于不同强度的降雨区，()会随着雨滴直径的增大而增大，导致()在距离廓线上表现出短距离内的波动。由式(2)可以看出，当()近似为常数或能被忽略时，可以将()作为()的估计值，否则就需要滤除()，估计准确的()。

1.1.1 EMD分解

EMD算法是一种自适应的信号去噪方法，该算法可将总差分传播相移自适应分解为多个IMF和残余项之和，且每个IMF都应满足以下2个条件：1)在整个曲线中，极值点和过零点的数目相等或至多相差1个；2)在任意局部区间，曲线的局部极大值包络线和局部极小值包络线的平均值为0。EMD分解包括提取分量、筛选IMF、计算余项三个部分：

(a) 提取分量

确定原始序列={(1),(2),…,(),…,()}的局部极值序列，记为={()|⊂}，其中，表示所有极值点的集合。局部极值序列分为局部极大值序列与局部极小值序列(例如，局部极大值定义为序列中的某个距离门所对的总差分传播相移值，其前一距离门的值比它小，后一距离门的值也比它小)。

从序列中提取相邻极值点()、()之间的局部序列设为={(),…,(),…,()}，其中，≤≤。根据文献[14]，求的局部均值()为

(3)

=-

(4)

(b) 筛选IMF

若符合IMF条件，就是的第一个IMF；若不满足IMF的条件，继续将作为新的信号，重复步骤(a)进行分解。在实际中，通常无法满足IMF条件，假设经过(一般小于10)次分解后，有以下式(5)成立，则认为其满足IMF条件：

(5)

式中，为门限，一般取值为02～03。

那么原始信号的第一个IMF可记为

(6)

(c) 计算余项

将得到的从中分离出来，得到余项=-，将作为新的信号，重复步骤(a)～(b)继续进行分解，依次可得到第2，第3，…，第个，…，第个IMF，记为

(7)

当余项变成一个单调信号时，EMD分解过程停止。

通过上述分解过程，可将分解为个IMF与余项的和，记为

(8)

基于EMD的总差分传播相移分解流程如图1所示。

图1 总差分传播相移分解流程图

1.1.2 EMD重构

传统的EMD凭借经验选取前几个IMF重构差分传播相移，这会导致重构的差分传播相移中混入多余的后向散射或漏掉某些有用信息，为了克服这一问题，本文提出利用皮尔逊相关系数准则(||)分析原始信号与EMD分解得到的个IMF相关性，确定最佳IMF数目。

首先计算=((1),(2),…,(),…,())与=((1),(2),…,(),…,())的值，定义式为

=1,2,…,

(9)

表1 皮尔逊相关系数评估标准

通过计算EMD分解得到的个IMF与原始信号的||值，找到合适的分界点∈(0,1)，确保除去噪声较大的前(∈(1,2,…,))个，其他IMF的||值稳定在评估标准范围内。

经过EMD分解后，总差分传播相移可以表示为

(10)

式中，～表示噪声分量，+1～表示信号分量，表示余项。差分传播相移由信号IMF与余项共同构成，即

(11)

基于EMD的差分传播相移重构流程如图2所示。

图2 差分传播相移重构示意图

1.2 基于最小二乘法的差分传播相移率估计

在降水区距离内含有个距离库，第个距离库对应值为()，用最小二乘法拟合这个库对应的，根据拟合曲线的斜率即可得到雨区距离的值,表示为

(12)

1.3 基于自适应约束方法的反射率衰减订正

衰减率对接收的反射率影响较大，且随着探测距离的增加而增加，造成反射率估计误差，因此本文在自适应约束算法的基础上，根据EMD方法处理得到的与最小二乘法拟合的，提出了适用于双偏振雷达的自适应约束方法进行反射率衰减订正。该方法首先对雨区径向数据划分区间，根据-的关系，计算所有区间的衰减系数，完成整个径向数据的衰减订正。

由文献[18]可知，雷达反射率衰减订正原理如式(13)所示：

(13)

式中，()与()分别为订正前后的第个距离库反射率，()为第距离库的衰减率。

(14)

(15)

(16)

最后将最优()代入式(13)中得到订正后的()。

对反射率处理的具体流程如图3所示。

图3 衰减订正处理流程图

2 实验结果及分析2.1 差分传播相移估计效果分析

本文使用的是ARM(Atmospheric Radiation Measurement Climate Research Facility)探测到的X波段天气雷达数据。该雷达采用双偏振体制，同时在水平和垂直偏振中传输，其主要包括偏振参数，，等。下面通过仿真试验分析了不同滤波方法的估计效果。

图4为X-SAPR雷达于0.5°俯仰角、30°方位角的雷达观测数据，以及经过不同滤波方法的径向距离廓线图。从图中可以看出，原始信号(图4(a))由于外界噪声存在较大的波动，但基本保持递增趋势。由均值滤波(图4(b))与FIR滤波(图4(c))后的处理效果可以看出，这两种方法对的波动起伏有较好的平滑效果，但随着平滑库数增大，不能有效抑制波动。卡尔曼滤波(图4(d))由于缺乏先验数据，造成初始数据失真。采用EMD方法对原始信号进行分解得到7个IMF，从低到高依次计算各IMF与原信号的||值分别为0.07，0.13，0.15，0.20，0.23，0.36，0.58，可知前3个IMF的||值均在0.00～0.19极弱相关这一范围内，为噪声分量，后4个IMF的||值均大于0.19，为信号分量，对后4个IMF进行有效重构，其结果如图4(e)所示，可以看出，数据的毛刺和波动均得到很好抑制，数据的连续性和平滑度有了显著提升，即去除了后向散射的影响，这也有利于对数据进行后续处理及应用。综上可见，EMD方法较之常规方法具有更明显的处理效果。

(a) 原始数据

图5为X-SAPR雷达于0.5°俯仰角观测数据的PPI图及经过EMD方法后的PPI图。从图中可见，在原始数据PPI图(图5(a))中，大部分径向都具有较为明显的层次性，颜色上表现为分层渐进递增形式，但也存在许多“麻点”(黑框所选区域为麻点)，即为波动数据点；经过EMD方法后(图5(b))，数据整体层次性得到明显加强，颜色分布均匀，且原始数据的有效的气象回波也得到了很好的保留。

(a) Ψdp原始数据

为了进一步分析不同滤波方法对差分传播相移估计的影响，分别通过平均波动指数、相对误差、互相关系数进行降雨估计性能分析，计算公式为

(17)

(18)

(19)

式中，()表示实测值，()表示经过处理后的值，通过来比较距离廓线的波动情况，的值越大，说明径向距离廓线的波动性越大；通过来比较数据可靠性，值越小数据的可信度越高；通过来比较数据变化趋势的一致程度，值越大表明线性一致程度越高。

表2统计了()经本文提出的方法与常规方法处理并应用于降水估计()的、和。根据统计结果可以看出，()经常规方法处理后的和较大，而EMD方法处理的()的与较小，且值相对较大，说明EMD方法对()数据的处理波动较小，可靠性高，应用于降水估测的精度最高。

表2 不同滤波方法对QPE影响的对比图

2.2 差分传播相移率估计效果分析

图6为X-SAPR雷达于0.5°俯仰角、30°方位角的雷达观测数据，以及经过不同滤波方法的径向距离廓线图。结果表明，经过均值滤波、FIR滤波、卡尔曼滤波和EMD方法处理后估计的的负值数量分别为182，185，167和148。说明EMD方法能够有效地减少的负值，保证数据的基本趋势。

(a) 原始数据

2.3 反射率衰减订正效果分析

X波段双偏振雷达探测的径向距离为40 km，一个距离库表示100 m。图7为同一径向的反射率衰减订正结果对比图，从图中可以看到，变化与的变化趋势是相对应的。在初始距离库，雨区的衰减较少，反射率基本在20dBz左右，而在第200个距离库处有一强回波区域，雨区的衰减较多，使用EMD方法进行订正后的反射率因子能够在保持原始数据变化趋势的同时，对其进行有效的衰减订正。

图7 反射率衰减订正对比图

图8为反射率因子订正前后PPI对比图。由图8(a)可知，强回波的区域分布离散，主要集中在图中心位置，雷达强回波反射率因子一般在30～40 dBz，由图8(b)可知，经过衰减订正后，雷达回波反射率因子的订正值一般在1～5 dBz。强回波区域扩大，回波中心更为明显，订正值达到了5～10 dBz。

由于S波段雷达基本上不存在雨区衰减，除非在一些湿雹区。因此将同时间同一区域内的反射率因子订正后的效果与S波段数据进行对比，其效果与S波段探测结果越相似，则认为订正效果越好。S波段双偏振雷达的探测距离范围为230 km，而X波段的探测距离范围为40 km，所以X波段气象雷达包含在S波段双偏振雷达的探测范围内，图8(c)为S波段双偏振雷达回波反射率因子，黑色方框中的区域作为X波段双偏振雷达衰减订正后的效果参照图。对比图8(a)、图8(c)可以看到，X波段双偏振雷达衰减更为严重。对比图8(b)、图8(c)可以看到，经过EMD方法处理后，X波段双偏振雷达的反射率因子中心处强回波有所加强，与S波段雷达回波反射率因子差异显著减小，表明EMD方法对X波段双偏振雷达衰减订正具有一定的效果。

(a) X波段ZH原始数据

3 结束语

本文提出了基于EMD方法的X波段双偏振雷达衰减订正，首先使用经验模式分解方法将总差分传播相移按照频率分布自适应地分解为多个IMF，使得原始序列平稳化，从而提高差分传播相移的预测精度，接着利用最小二乘法估计差分传播相移率，最后根据得到的差分传播相移与差分传播相移率采用自适应约束算法对反射率因子进行衰减订正。实验结果表明，利用EMD方法能够实现对差分传播相移滤波处理，并且保留原始信号的基本趋势；利用自适应约束的方法能够对反射率进行有效的衰减订正，这对雷达降水量估计以及降水分类等有重要的现实意义，所以本方法具有一定的实际应用价值。