对敏视达雷达回波进行基于PHIDP的dBZ和ZDR订正_2014年4月5日～18日.doc

对敏视达C波段雷达进行基于PHIDP的dBZ和ZDR订正研究2014年4月5日18日目 录1 雷达自身PHIDP的变化11.1 一个月以来的变化11.2 相连体扫之间的变化31.3 小结52 PHIDP的初始值的修正52.1 从垂直90度的回波数据中提取PHIDP的初始值（Fun_PHIDP_Initial_Calc.m）52.1.1 计算方法52.1.2 计算结果62.1.3 小结92.2 初始值修正（Fun_PHIDP_Initial_Modify.m）92.2.1 计算方法92.2.2 计算结果93 对PHIDP进行各种处理153.1 剔除异常的PHIDP值(Fun_PHIDP_Reject_Singular_Data.m)153.1.1 剔除准则的研究153.1.2 计算结果173.2 去折叠（Fun_PHIDP_Remove_Fold.m）203.2.1 计算方法203.2.2 计算结果203.3 平滑滤波(Fun_PHIDP_Filter.m)223.3.1 确定参与平滑滤波的点数223.3.2 平滑滤波计算243.3.3 计算结果253.3.4 小结283.4 对一些问题的说明283.4.1 为何不采用FIR滤波器对PHIDP进行平滑滤波283.4.2 一元线性回归和鲁棒线性回归的比较283.4.3 对前后多少距离范围，再相邻10点的平滑的说明314 利用PHIDP进行dBZ和ZDR修正354.1 计算方法(Fun_PHIDP_Modify_dBZ_dBT_ZDR.m)354.2 NJU.20140329.000632.AR2.bz2，仰角=1.5度的计算结果354.2.1 PHIDP处理前后的对比354.2.2 dBZ修正前后对比374.2.3 ZDR修正前后对比384.2.4 RVP9给出的KDP和重新计算的KDP 对比394.3 NJU.20140329.030007.AR2.bz2，仰角=2.4度的计算结果414.3.1 PHIDP处理前后的对比424.3.2 dBZ修正前后对比434.3.3 ZDR修正前后对比444.3.4 RVP9给出的KDP和重新计算的KDP 对比455 总结465.1 基于PHIDP修正ZDR和PHIDP很重要，否则dBZ和ZDR的结果会有很大的偏差465.2 由于雷达采用接收机恒温技术，自身的PHIDP比较稳定，因此可以在进行体扫的时候，利用垂直90度观测的数据（如果有合适的降水回波数据的话），来进行初始值的计算，不需要每个径向上都进行初始值的统计465.3 通过垂直90度的回波计算的PHIDP初始值非常精确，不受四周地物等的影响465.4 通过对PHIDP进行精心的处理（剔除异常值、去折叠、自适应的平滑滤波等），可以很好的去除PHIDP受地物等非降水回波的影响，并能大大降低PHIDP的随机误差475.5 根据PHIDP，重新计算的KDP，看起来效果比RVP9信号处理器给出的KDP值，要好得多。接下来终于可以用KDP的值进行定量降水的研究了。476 下一步工作476.1 将上面的一整套处理方法，全面应用到该雷达至今所观测到的所有回波数据之中，进行全面的验证，主要是要看对于各个不同情况的回波，算法不能出现错误或奇异的结果。476.2 根据PHIDP订正之后的dBZ、ZDR，和S波段的雷达进行对比476.3 重新计算得到的KDP，和雨量计进行对比476.4 从IQ层面，开展对地物等非气象回波的研究，看能否找到更加合理、科学的剔除异常PHIDP的准则4731 雷达自身PHIDP的变化1.1 一个月以来的变化由于这部雷达的双通道接收机采用了半导体恒温设计，因此PHIDP随时间、环境温度的变化应该是比较小的。下面是最近一个月以来，PHIDP的变化（垂直90度扫描观测）图 2月18日降雪回波的PHIDP图 2月28日降水回波的PHIDP图 3月18日降水回波的PHIDP图 3月29日降水回波的PHIDP注意：这几张图的PHIDP的显示范围比较精细，从-10+30度，因此虽然从画面上看，有很多随机的斑点，但实际上数据的方差是很小的。从上面的四张图可以看出，在一个月以来，雷达自身PHIDP初始值的变化为： 约8度-4度-12度-14度。这应该说是比较小的。1.2 相连体扫之间的变化下面三张是相邻体扫之间PHIDP的值：可以看出，PHIDP的初始值在间隔几分钟的时间内，几乎没有变化。1.3 小结从上面的回波图可以看出，这部雷达自身PHIDP的变化是比较小的。这也就说明，我们在进行PHIDP初始值的提取的时候，不需要每个径向都进行一次，也不需要每个天线周期都进行一次，而是可以在进行体扫的时候，利用垂直90度观测的数据（如果有合适的降水回波数据的话），来进行初始值的计算。2 PHIDP的初始值的修正2.1 从垂直90度的回波数据中提取PHIDP的初始值（Fun_PHIDP_Initial_Calc.m）2.1.1 计算方法首先，从零距离开始，去除最近的3个距离单元（一般是发射机的漏信号），连续10个距离库（注意：如果距离的分辨率比较大，则不需要10个，5个即可。总之，高度不要超过1500m）的CC0.9，且SNR20dB，则将这10个库的PHIDP进行平均（起到平滑的作用，降低随机起伏），作为该径向上的初始值。如果该径向上没有符合的结果，则为NaN（无效值）。然后，将全方位上所有的径向上算出的PHIDP值，去除无效值，再去除最大和最小的10%，求出平均值。如果样本数足够多(100)，则将平均值作为最终的结果输出。注意：对相位进行平均时，不能简单的将所有的值加起来求平均，而是要考虑到相位折叠的问题。计算方法见matlab程序。由于前面提到了：雷达自身的PHIDP初始值变化很小，因此在这里进行计算的时候，我们就可以宁缺毋滥，去掉所有可疑的回波。2.1.2 计算结果对几个不同时间的基数据进行计算的结果如下（从文件名中可以看出数据的日期）：注意：上面的最后两张图片是相邻2个体扫周期的图，可见：相邻2个体扫周期的PHIDP 的初始值，几乎没有变化。2.1.3 小结l 通过垂直90度的回波计算的PHIDP初始值非常精确，不受四周地物等的影响；l PHIDP初始值随着天线的转动有起伏（这是这部雷达一直存在的现象，估计与波导旋转铰链机械结构有关），起伏量为正负3度。l 起伏量为3度，会造成dBZ的订正误差为：0.11*3=0.33dB；ZDR的订正误差为：0.033*3=0.1dB。可见，误差比较小，一般可以忽略。如果为了精确修正，则需要考虑这个起伏量。但目前我们就不考虑这个起伏了，只取平均值来作为PHIDP的初始值。2.2 初始值修正（Fun_PHIDP_Initial_Modify.m）2.2.1 计算方法用matlab计算非常简单，代码如下：PHIDP.Data=mod( PHIDP.Data -PHIDP_Initial_Result ,360); %先转换到0360度之间PHIDP.Data( PHIDP.Data 180 ) =PHIDP.Data( PHIDP.Data 180 ) -360; %将凡是180度的，减去360度，换算到-180180之间PHIDP.Data( PHIDP.Data -180 ) =PHIDP.Data( PHIDP.Data 180度的，减去360度，换算到-180180之间2.2.2 计算结果下面对各个仰角层面的PHIDP值，利用前面计算的PHIDP的初始值进行修正。图 21 垂直90度的回波（修正前）图 22 垂直90度的回波（初相修正后）图 23 仰角=10度的回波（修正前）图 24 仰角=10度的回波（初相修正后）图 25 仰角=1.5度的回波（修正前）图 26仰角=1.5度的回波（初相修正后）从上面对比的6张图可以看出，利用前面垂直90度扫描得到的PHIDP的初始值，可以很精确的将雷达在各个层面上扫描得到的PHIDP的初始相位修正正确。下面我们用A显的方式，仔细观察一下各个不同的仰角层面，用垂直90度扫描得到的PHIDP的初始值进行修正之后的初始相位是否为零度。图 27仰角=0.5度的回波（初相修正后）图 28仰角=1.5度的回波（初相修正后）图 29仰角=4.3度的回波（初相修正后）图 210仰角=10.0度的回波（初相修正后）从上面几张A显图可以看出，不同仰角层面的PHIDP的初始相位是基本相同的。有一点差别，这是由于仰角的波导铰链，也会像方位波导铰链一样，会对PHIDP有点影响。当然，对于那些将固态发射机安装在天线背面的固态体制雷达而言，由于没有了方位和仰角的波导铰链，因此PHIDP肯定就不再受到方位和仰角的影响了。3 对PHIDP进行各种处理3.1 剔除异常的PHIDP值(Fun_PHIDP_Reject_Singular_Data.m)3.1.1 剔除准则的研究下面是几张不同文件中，在不同的径向上，PHIDP和CC、SNR的关系：图 31 在某条径向上，PHIDP、CC、SNR的关系（1）图 32 在某条径向上，PHIDP、CC、SNR的关系（2）图 33 在某条径向上，PHIDP、CC、SNR的关系（3）通过对各种不同情况下的PHIDP和CC、SNR的波形对比，我们决定采用以下的准则：l 如果当前距离单元和前一个距离单元的PHIDP相差如果超过30度，则剔除； l 如果当前距离单元的SNR3dB，或者CC40dB且CC0.9的数据且dBZ50dB，属于数据比较好的点，采取前后1km内的点参与计算； l 对于其它的数据（低SNR、低CC），属于数据比较差的点，采取前后更多点的点（4km）进行滤波；l 对于处于中间的点，则采用前后2km内的点参与计算；但是，我们发现采用分为几种类型的方法，一方面难以控制分类的门限准则；而且在回波数据处于分类的交接区的时候，会造成计算出来的PHIDP有一些不合理的振荡现象。因此，经过多次试验，我们最终采用了一种新的办法来得到需要前后多少点，该方法是根据CC、SNR和dBZ的值，自动计算出需要前后多少距离范围内的点参与计算。关键代码如下： if (CC0.8) if (SNR20 ) if dBZ50 DataLength_km=1; elseif dBZ20 DataLength_km=5; elseif SNR0 DataLength_km=10; elseif isnan(SNR) DataLength_km=10-SNR/4; else DataLength_km=10; end end else DataLength_km=10; end;经过上述计算得出的“前后多少距离范围”，还要进行相邻10点的平滑，以便得出一个更稳健的距离范围（详见第3.4.3章）。下图是根据某条径向上的SNR和dBZ的值，计算出来的需要前后多少距离范围内的点参与后续的PHIDP平滑滤波计算：图 312 自动计算出来的、需要前后多少距离范围内的点参与计算图 313 与前一张图对应的 PHIDP与SNR的图从上图可以看出，如果dBZ或SNR比较高，则参与平滑滤波计算的点数比较少；反之，则参与计算的点数比较多。3.3.2 平滑滤波计算下面要遍历每个距离点，用该距离点前后一定区间内的PHIDP值，进行平滑滤波计算。注意：这段代码比较耗费时间，以后需要优化。但如果用C+实现，则速度应该很快的。在进行计算的时候，需要注意如果参与计算的点数20，则说明有效的点数太少了，为了计算的精确，需要将参与运算的前后距离范围增加1倍。另外，如果发现在当前距离点的左边或者右边没有足够的数据（至少要有4个数据才行），则还要将参与运算的前后距离范围再增加1倍。平滑滤波的算法有两种：l 用Matlab的polyfit函数实现的常用的一元线性回归；l 用Matlab提供的robustfit实现的鲁棒线性回归；一元线性回归的主要代码如下： p=polyfit(KDP_Calc_Index_No_NaN,KDP_Calc_PHIDP_No_NaN,1); PHIDP_Select_Temp(ii)=polyval(p,ii); KDP_Select_Temp(ii)=p(1)/PHIDP.Resolution*1e3;鲁棒线性回归的主要代码如下： brob = robustfit(KDP_Calc_Index_No_NaN,KDP_Calc_PHIDP_No_NaN); PHIDP_Select_Temp(ii)=brob(1)+brob(2)*ii; KDP_Select_Temp(ii)=brob(2)/PHIDP.Resolution*1e3;经试验，这两种算法的效果都很好。但是，一元线性回归的计算速度快得多。3.3.3 计算结果图 314 在某条径向上，PHIDP的平滑滤波前后对比图 315 在某条径向上，PHIDP的平滑滤波前后对比（近距离处地物干扰）图 316 在某条径向上，PHIDP的平滑滤波前后对比（强回波，高斜率）图 317 在某条径向上，PHIDP的平滑滤波前后对比（远处弱信号）3.3.4 小结通过PHIDP在不同距离段上、平滑滤波的表现，我们可以发现：通过合理、科学的确定需要前后多少距离范围内的点参与平滑滤波计算，然后采用线性回归的方法，可以很好的实现PHIDP的滤波任务。同时，线性回归还能将前面剔除的、异常PHIDP值所在的距离单元上的值，用线性内插的方式填充起来。3.4 对一些问题的说明3.4.1 为何不采用FIR滤波器对PHIDP进行平滑滤波我们认为，FIR滤波器用于PHIDP的平滑滤波存在以下的缺陷：l 由于FIR滤波器不具备对奇异的点的抵抗能力，因此万一有一个奇异点（由于地物抑制不干净、非气象回波等原因），则会严重影响计算的结果。而线性回归的算法就具备这个抵抗奇异点的能力。l 由于前面进行了“剔除异常的PHIDP值”的处理，那些异常值的距离单元的PHIDP，也要想办法给恢复出来（如果这些距离单元的PHIDP值直接扔掉了，最终的dBZ和ZDR的画面上就会出现比较多的缺失）。而FIR滤波器难以完成这个插值的功能，只有本文中提到的线性回归才适合完成这件事情；l 由于需要根据不同的数据质量，选择不同的参与平滑滤波的点数，而FIR滤波器难以做到根据点数而变化（因为FIR的阶数一般都是预先设计好的）；l 还有一条：线性回归的计算量也要比FIR滤波器要少。3.4.2 一元线性回归和鲁棒线性回归的比较我们挑选了一组PHIDP起伏很大的数据，来比较两种方法的效果。用Matlab的polyfit函数实现的常用的一元线性回归的效果如下：图 318 polyfit函数实现的常用的一元线性回归的效果（1）用Matlab提供的robustfit实现的鲁棒线性回归的效果如下；图 319 robustfit实现的鲁棒线性回归的效果（1）再换一组PHIDP起伏很大的数据，再比较两种方法的效果：图 320 polyfit函数实现的常用的一元线性回归的效果（2）图 321 robustfit实现的鲁棒线性回归的效果（2）对比两种方法发现，这两种算法的效果都很好，能很好的起到平滑滤波的作用，基本不受异常点的影响。但是，一元线性回归的计算速度快得多，而且容易用C+ 实现。因此，今后的基于C+的信号处理程序，将采用polyfit的一元线性回归来实现。3.4.3 对前后多少距离范围，再相邻10点的平滑的说明在确定参与平滑滤波的点数的时候，我们提到：自动计算出需要前后多少距离范围之后，还需要对前后多少距离范围，进行相邻10点的平滑，以便得出一个更稳健的距离范围。那么，为何还要仅有相邻10点的平滑呢？先绘制一张不经过相邻10点平滑，直接就是自动计算出的需要前后多少距离范围的图（横轴是距离，竖轴是需要前后多少距离范围）：图 322 计算出的需要前后多少距离范围的图下图是此时，经过PHIDP滤波后的图：这里也有一点振荡PHIDP出现不合理的振荡现象图 323 根据前后多少距离范围，对PHIDP进行滤波后的图从上面两张图可以看出，由于回波的SNR、dBZ天生就存在起伏，从而就造成了“前后多少距离的点参与后续计算”也会发生很大的起伏。而这个起伏就会造成对PHIDP的平滑滤波效果大大下降，会出现不合理的振荡现象。注意：在国内外的文献中，一般采用根据数据的质量，将参与计算的前后距离的长度分为几种不同的情况。由于回波天生存在的起伏，当回波数据处于分类的交接区的时候，也会造成平滑滤波后的PHIDP出现不合理的振荡现象。下面是对自动计算出的前后距离，进行10点平滑之后的图（横轴是距离，竖轴是需要前后多少距离范围）：图 324 计算出的需要前后多少距离范围，再进行10点平滑的图下图是此时，经过PHIDP滤波后的图：图 325 根据10点平滑后的距离范围，对PHIDP进行滤波后的图从上面两张图可以看出，经过对“前后多少距离范围”经过10点平滑之后，可以得出一个更稳健的距离范围。然后再以这个距离范围内的PHIDP数据点，进行平滑滤波，效果就非常好了。当然，对“前后多少距离范围再进行10点平滑，得到更稳健的距离范围”这件事情并不是非常必需的。也就是说，即使不采取这个措施，经过平滑滤波之后的PHIDP，以及推算出的KDP，已经比信号处理器直接给出的PHIDP和KDP要好的多了。4 利用PHIDP进行dBZ和ZDR修正4.1 计算方法(Fun_PHIDP_Modify_dBZ_dBT_ZDR.m)这里的计算就很简单了。主要代码如下：dBZ.Data=dBZ.Data+PHIDP *0.11; ZDR.Data=ZDR.Data+PHIDP *0.033; 在该程序中，0.11和0.033分别是dBZ和ZDR相对于PHIDP的修正系数，见有关论文。该系数和雷达波段有关。我们先直接选用国外的研究成果，以后再研究更加准确的系数值。下面以两组回波为例，展示一下处理的效果。更多的处理效果，直接用以下Matlab程序，进行计算吧：Read_MSD_Radar_BaseData_PHIDP_Correct.m4.2 NJU.20140329.000632.AR2.bz2，仰角=1.5度的计算结果4.2.1 PHIDP处理前后的对比图 41 雷达信号处理器直接给出的PHIDP结果图 42 经过初相修正、数据剔除、去折叠、平滑滤波之后的PHIDP结果从上面两张图对比可以看出，经过初始相位修正、数据剔除、去折叠、平滑滤波之后的PHIDP的数据质量，要远远好于雷达信号处理器直接给出的PHIDP结果。4.2.2 dBZ修正前后对比图 43 雷达信号处理器直接给出的dBZ结果dBZ增加了图 44 经过PHIDP订正之后的dBZ结果4.2.3 ZDR修正前后对比图 45 雷达信号处理器直接给出的ZDR结果ZDR增加了图 46 经过PHIDP订正之后的ZDR结果4.2.4 RVP9给出的KDP和重新计算的KDP 对比图 47 RVP9给出的KDP结果强回波区强回波区图 48 从PHIDP重新计算的KDP结果下面我们对方位为298度的PHIDP和KDP结果绘制A显如下，可以更明显的看出两者的差异：通过人工核算，此处的斜率约为：5deg/km图 49 方位为298度的PHIDP的值RVP9信号处理给出的KDP波动大RVP9信号处理给出的KDP有偏差图 410 方位为298度的两种KDP结果的对比注意：这个偏差的原因很简单：RVP9处理器给出的KDP是单程的，而Matlab计算的KDP是双程的。单程的KDP是双程的1/2。4.3 NJU.20140329.030007.AR2.bz2，仰角=2.4度的计算结果4.3.1 PHIDP处理前后的对比图