对敏视达雷达回波进行基于phidp的dbz和zdr订正.docx
《对敏视达雷达回波进行基于phidp的dbz和zdr订正.docx》由会员分享，可在线阅读，更多相关《对敏视达雷达回波进行基于phidp的dbz和zdr订正.docx（52页珍藏版）》请在冰豆网上搜索。
对敏视达雷达回波进行基于phidp的dbz和zdr订正
对敏视达C波段雷达进行基于PHIDP的dBZ和ZDR订正研究
2014年4月5日～18日
1 雷达自身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的初始值（）
2.1.1 计算方法
首先，从零距离开始，去除最近的3个距离单元（一般是发射机的漏信号），连续10个距离库（注意：
如果距离的分辨率比较大，则不需要10个，5个即可。
总之，高度不要超过1500m）的CC>，且SNR>20dB，则将这10个库的PHIDP进行平均（起到平滑的作用，降低随机起伏），作为该径向上的初始值。
如果该径向上没有符合的结果，则为NaN（无效值）。
然后，将全方位上所有的径向上算出的PHIDP值，去除无效值，再去除最大和最小的10%，求出平均值。
如果样本数足够多（>100），则将平均值作为最终的结果输出。
注意：
对相位进行平均时，不能简单的将所有的值加起来求平均，而是要考虑到相位折叠的问题。
计算方法见matlab程序。
由于前面提到了：
雷达自身的PHIDP初始值变化很小，因此在这里进行计算的时候，我们就可以宁缺毋滥，去掉所有可疑的回波。
2.1.2 计算结果
对几个不同时间的基数据进行计算的结果如下（从文件名中可以看出数据的日期）：
注意：
上面的最后两张图片是相邻2个体扫周期的图，可见：
相邻2个体扫周期的PHIDP的初始值，几乎没有变化。
2.1.3 小结
●通过垂直90度的回波计算的PHIDP初始值非常精确，不受四周地物等的影响；
●PHIDP初始值随着天线的转动有起伏（这是这部雷达一直存在的现象，估计与波导旋转铰链机械结构有关），起伏量为正负3度。
●起伏量为3度，会造成dBZ的订正误差为：
*3=；ZDR的订正误差为：
*3=。
可见，误差比较小，一般可以忽略。
如果为了精确修正，则需要考虑这个起伏量。
但目前我们就不考虑这个起伏了，只取平均值来作为PHIDP的初始值。
2.2 初始值修正（）
2.2.1 计算方法
用matlab计算非常简单，代码如下：
=mod（-PHIDP_Initial_Result,360）;%先转换到0~360度之间
（>180）=（>180）-360;%将凡是>180度的，减去360度，换算到-180~180之间
（<-180）=（<-180）+360;%将凡是>180度的，减去360度，换算到-180~180之间
2.2.2 计算结果
下面对各个仰角层面的PHIDP值，利用前面计算的PHIDP的初始值进行修正。
图21垂直90度的回波（修正前）
图22垂直90度的回波（初相修正后）
图23仰角=10度的回波（修正前）
图24仰角=10度的回波（初相修正后）
图25仰角=度的回波（修正前）
图26仰角=度的回波（初相修正后）
从上面对比的6张图可以看出，利用前面垂直90度扫描得到的PHIDP的初始值，可以很精确的将雷达在各个层面上扫描得到的PHIDP的初始相位修正正确。
下面我们用A显的方式，仔细观察一下各个不同的仰角层面，用垂直90度扫描得到的PHIDP的初始值进行修正之后的初始相位是否为零度。
图27仰角=度的回波（初相修正后）
图28仰角=度的回波（初相修正后）
图29仰角=度的回波（初相修正后）
图210仰角=度的回波（初相修正后）
从上面几张A显图可以看出，不同仰角层面的PHIDP的初始相位是基本相同的。
有一点差别，这是由于仰角的波导铰链，也会像方位波导铰链一样，会对PHIDP有点影响。
当然，对于那些将固态发射机安装在天线背面的固态体制雷达而言，由于没有了方位和仰角的波导铰链，因此PHIDP肯定就不再受到方位和仰角的影响了。
3 对PHIDP进行各种处理
3.1 剔除异常的PHIDP值
3.1.1 剔除准则的研究
下面是几张不同文件中，在不同的径向上，PHIDP和CC、SNR的关系：
图31在某条径向上，PHIDP、CC、SNR的关系
（1）
图32在某条径向上，PHIDP、CC、SNR的关系
（2）
图33在某条径向上，PHIDP、CC、SNR的关系（3）
通过对各种不同情况下的PHIDP和CC、SNR的波形对比，我们决定采用以下的准则：
●如果当前距离单元和前一个距离单元的PHIDP相差如果超过30度，则剔除；
●如果当前距离单元的SNR<3dB，或者CC<，则剔除；
完成上面两条准则之后，接下来还需要对PHIDP进行判断，要求前后的PHIDP都不是无效值，当前距离点才有效。
也就是说，若当前距离单元的值是有效的，则若前后1个距离单元的有一个无效，则将当前距离单元的值也设为无效。
另外，对对于第一个点和最后一点，要单独判断。
3.1.2 计算结果
处理结果如下：
图34在某条径向上，原始PHIDP和剔除虚假数据之后的对比
（1）
图35在某条径向上，原始PHIDP和剔除虚假数据之后的对比
（2）
图36在某条径向上，原始PHIDP和剔除虚假数据之后的对比（3）
图37在某条径向上，原始PHIDP和剔除虚假数据之后的对比（4）
可以看出，通过前面的判断准则，可以去除绝大多数地物、低SNR造成的错误的PHIDP值。
3.2 去折叠（）
3.2.1 计算方法
去折叠的方法如下（详见Matlab程序）：
●先选出所有不是NaN的数据，即有效的数据；
●然后计算相邻两个数据之差，依次判断差值，让差值在+-180度之间；
●然后采用了Matlab的一个巧妙的函数:
cumsum，就可以将相邻两个数据之差不断进行累加，就能完成去折叠的任务了
●最后，要将结果重新填到原来的数组中（原来数组中的NaN数据不变）
3.2.2 计算结果
下面我们故意将PHIDP加上一个初始值（这里是170度），以便能对去折叠的效果进行测试。
如果不加初始值的话，则看不出去折叠的效果。
（因为这几组数据的回波强度不够强，回波的PHIDP还没有达到180度的量级，因此PHIDP根本就没有折叠）
原始的PHIDP随距离的变化如下：
图38未经过去折叠处理的PHIDP
经过去折叠处理之后的图如下：
图39经过了去折叠处理的PHIDP
更换了另一组数据进行测试，原始的PHIDP（故意加了150度的初始值）随距离的变化如下：
图310未经过去折叠处理的PHIDP
经过去折叠处理之后的图如下：
图311经过了去折叠处理的PHIDP
可见，这种方法是能完成去折叠的处理的，并且即使数据中存在无效数据（无效数据是由于回波的SNR太低等原因，被剔除掉的），也不会影响去折叠的计算。
3.3 平滑滤波
3.3.1 确定参与平滑滤波的点数
总所周知，参与平滑滤波的点数越多，则平滑效果越好，但是会影响对小范围强回波的计算。
因此，首先要合理、科学的确定需要前后多少距离范围内的点参与平滑滤波计算。
在国内外的文献中，普遍采用根据数据的质量，将参与计算的前后距离的长度分为几种类型，例如：
●对于SNR>40dB且CC>的数据且dBZ>50dB，属于数据比较好的点，采取前后1km内的点参与计算；
●对于其它的数据（低SNR、低CC），属于数据比较差的点，采取前后更多点的点（4km）进行滤波；
●对于处于中间的点，则采用前后2km内的点参与计算；
但是，我们发现采用分为几种类型的方法，一方面难以控制分类的门限准则；而且在回波数据处于分类的交接区的时候，会造成计算出来的PHIDP有一些不合理的振荡现象。
因此，经过多次试验，我们最终采用了一种新的办法来得到需要前后多少点，该方法是根据CC、SNR和dBZ的值，自动计算出需要前后多少距离范围内的点参与计算。
关键代码如下：
if（CC>
if（SNR>20）
ifdBZ>50
DataLength_km=1;
elseifdBZ<20
DataLength_km=5;
elseif~isnan（dBZ）
DataLength_km=5-（dBZ-20）/;
else
DataLength_km=5;
end
else
ifSNR>20
DataLength_km=5;
elseifSNR<0
DataLength_km=10;
elseif~isnan（SNR）
DataLength_km=10-SNR/4;
else
DataLength_km=10;
end
end
else
DataLength_km=10;
end;
经过上述计算得出的“前后多少距离范围”，还要进行相邻10点的平滑，以便得出一个更稳健的距离范围（详见第Error!
Referencesourcenotfound.章）。
下图是根据某条径向上的SNR和dBZ的值，计算出来的需要前后多少距离范围内的点参与后续的PHIDP平滑滤波计算：
图312自动计算出来的、需要前后多少距离范围内的点参与计算
图313与前一张图对应的PHIDP与SNR的图
从上图可以看出，如果dBZ或SNR比较高，则参与平滑滤波计算的点数比较少；反之，则参与计算的点数比较多。
3.3.2 平滑滤波计算
下面要遍历每个距离点，用该距离点前后一定区间内的PHIDP值，进行平滑滤波计算。
注意：
这段代码比较耗费时间，以后需要优化。
但如果用C++实现，则速度应该很快的。
在进行计算的时候，需要注意如果参与计算的点数<20，则说明有效的点数太少了，为了计算的精确，需要将参与运算的前后距离范围增加1倍。
另外，如果发现在当前距离点的左边或者右边没有足够的数据（至少要有4个数据才行），则还要将参与运算的前后距离范围再增加1倍。
平滑滤波的算法有两种：
●用Matlab的polyfit函数实现的常用的一元线性回归；
●用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）/*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）/*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的平滑滤波存在以下的缺陷：
●由于FIR滤波器不具备对奇异的点的抵抗能力，因此万一有一个奇异点（由于地物抑制不干净、非气象回波等原因），则会严重影响计算的结果。
而线性回归的算法就具备这个抵抗奇异点的能力。
●由于前面进行了“剔除异常的PHIDP值”的处理，那些异常值的距离单元的PHIDP，也要想办法给恢复出来（如果这些距离单元的PHIDP值直接扔掉了，最终的dBZ和ZDR的画面上就会出现比较多的缺失）。
而FIR滤波器难以完成这个插值的功能，只有本文中提到的线性回归才适合完成这件事情；
●由于需要根据不同的数据质量，选择不同的参与平滑滤波的点数，而FIR滤波器难以做到根据点数而变化（因为FIR的阶数一般都是预先设计好的）；
●还有一条：
线性回归的计算量也要比FIR滤波器要少。
3.4.2 一元线性回归和鲁棒线性回归的比较
我们挑选了一组PHIDP起伏很大的数据，来比较两种方法的效果。
用Matlab的polyfit函数实现的常用的一元线性回归的效果如下：
图318polyfit函数实现的常用的一元线性回归的效果
（1）
用Matlab提供的robustfit实现的鲁棒线性回归的效果如下；
图319robustfit实现的鲁棒线性回归的效果
（1）
再换一组PHIDP起伏很大的数据，再比较两种方法的效果：
图320polyfit函数实现的常用的一元线性回归的效果
（2）
图321robustfit实现的鲁棒线性回归的效果
（2）
对比两种方法发现，这两种算法的效果都很好，能很好的起到平滑滤波的作用，基本不受异常点的影响。
但是，一元线性回归的计算速度快得多，而且容易用C++实现。
因此，今后的基于C++的信号处理程序，将采用polyfit的一元线性回归来实现。
3.4.3 对前后多少距离范围，再相邻10点的平滑的说明
在确定参与平滑滤波的点数的时候，我们提到：
自动计算出需要前后多少距离范围之后，还需要对前后多少距离范围，进行相邻10点的平滑，以便得出一个更稳健的距离范围。
那么，为何还要仅有相邻10点的平滑呢？
先绘制一张不经过相邻10点平滑，直接就是自动计算出的需要前后多少距离范围的图（横轴是距离，竖轴是需要前后多少距离范围）：
图322计算出的需要前后多少距离范围的图
下图是此时，经过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 计算方法
这里的计算就很简单了。
主要代码如下：
=+PHIDP*;
=+PHIDP*;
在该程序中，和分别是dBZ和ZDR相对于PHIDP的修正系数，见有关论文。
该系数和雷达波段有关。
我们先直接选用国外的研究成果，以后再研究更加准确的系数值。
下面以两组回波为例，展示一下处理的效果。
更多的处理效果，直接用以下Matlab程序，进行计算吧：
4.2 NJU..，仰角=度的计算结果
4.2.1 PHIDP处理前后的对比
图41雷达信号处理器直接给出的PHIDP结果
图42经过初相修正、数据剔除、去折叠、平滑滤波之后的PHIDP结果
从上面两张图对比可以看出，经过初始相位修正、数据剔除、去折叠、平滑滤波之后的PHIDP的数据质量，要远远好于雷达信号处理器直接给出的PHIDP结果。
4.2.2 dBZ修正前后对比
图43雷达信号处理器直接给出的dBZ结果
图44经过PHIDP订正之后的dBZ结果
4.2.3 ZDR修正前后对比
图45雷达信号处理器直接给出的ZDR结果
图46经过PHIDP订正之后的ZDR结果
4.2.4 RVP9给出的KDP和重新计算的KDP对比
图47RVP9给出的KDP结果
图48从PHIDP重新计算的KDP结果
下面我们对方位为298度的PHIDP和KDP结果绘制A显如下，可以更明显的看出两者的差异：
图49方位为298度的PHIDP的值
图410方位为298度的两种KDP结果的对比
注意：
这个偏差的原因很简单：
RVP9处理器给出的KDP是单程的，而Matlab计算的KDP是双程的。
单程的KDP是双程的1/2。
4.3 NJU..，仰角=度的计算结果
4.3.1 PHIDP处理前后的对比
图411雷达信号处理器直接给出的PHIDP结果
图412经过初相修正、数据剔除、去折叠、平滑滤波之后的PHIDP结果
4.3.2 dBZ修正前后对比
图413雷达信号处理器直接给出的dBZ结果
图414经过PHIDP订正之后的dBZ结果
4.3.3 ZDR修正前后对比
图415雷达信号处理器直接给出的ZDR结果
图416经过PHIDP订正之后的ZDR结果
4.3.4 RVP9给出的KDP和重新计算的KDP对比
图417RVP9给出的KDP结果
图418从PHIDP重新计算的KDP结果
5 总结
5.1 基于PHIDP修正ZDR和PHIDP很重要，否则dBZ和ZDR的结果会有很大的偏差
5.2 由于雷达采用接收机恒温技术，自身的PHIDP比较稳定，因此可以在进行体扫的时候，利用垂直90度观测的数据（如果有合适的降水回波数据的话），来进行初始值的计算，不需要每个径向上都进行初始值的统计
5.3 通过垂直90度的回波计算的PHIDP初始值非常精确，不受四周地物等的影响
5.4 通过对PHIDP进行精心的处理（剔除异常值、去折叠、自适应的平滑滤波等），可以很好的去除PHIDP受地物等非降水回波的影响，并能大大降低PHIDP的随机误差
5.5 根据PHIDP，重新计算的KDP，看起来效果比RVP9信号处理器给出的KDP值，要好得多。
接下来终于可以用KDP的值进行定量降水的研究了。
6 下一步工作
6.1 将上面的一整套处理方法，全面应用到该雷达至今所观测到的所有回波数据之中，进行全面的验证，主要是要看对于各个不同情况的回波，算法不能出现错误或奇异的结果。
Matlab程序为：
计算后生成的PHIDP、dBZ、ZDR、KDP的图片，保存在“xxx_New”目录下。
6.2 根据PHIDP订正之后的dBZ、ZDR，和S波段的雷达进行对比
6.3 重新计算得到的KDP，和雨量计进行对比
6.4 从IQ层面，开展对地物等非气象回波的研究，看能否找到更加合理、科学的剔除异常PHIDP的准则
升级会员 