南京大学C波段双偏振雷达_用Matlab将所有IQ数据全部重新计算_2014年9月11日

南京大学C波段双偏振雷达_用Matlab将所有IQ数据全部
重新计算_2014年9月11日
南京大学C波段双偏振雷达用Matlab将所有IQ数据
全部重新计算
南京大学
2014年9月11日
目录
1 概
述 ..................................................................... (1)
2 整个的计算步
骤 ..................................................................... (2)
3 Matlab和C语言联合编
程 ............................................................... 3 3.1 概述...................................................................... ............................................ 3 3.2 编程方法...................................................................... (3)
3.2.1 编译器的设置...................................................................... . (3)
3.2.2 IPP的编
译 ..................................................................... ....................... 4 3.3 在64bit 操作系统下的编
译 ..................................................................... ..... 4 3.4 速度测试...................................................................... (5)
3.4.1 计算机配置...................................................................... .. (5)
3.4.2 处理一个完整的体扫的IQ数据的速
度 (5)
3.4.3 mexFun_Radar_Data_DSP 函数中，各个子函数的速
度 (6)
3.4.4 多个Matlab进程同时运
行 (6)
3.5 经验和教训...................................................................... .. (7)
4 IQ计算时的参数设
置 ......................................................................
7 4.1 雷达的性能参数...................................................................... ........................ 7 4.2 信号处理算法设定的参
数...................................................................... .. (8)
5 信号处理算法介
绍 ..................................................................... ..... 9 5.1 第1个子过程的算法介绍...................................................................... .. (9)
5.1.1
Read_MSD_Radar_IQData_And_Process_Main.m (9)
5.1.2
Fun_Read_MSD_Radar_IQData_And_DSP.m (10)
5.1.3
Fun_Radar_Data_DSP.m ................................................... .. (10)
5.1.3.1
Fun_PreProcess.m........................................................ (11)
5.1.3.2
Fun_ClutterFilter.m..................................................... . (11)
5.1.3.3
Fun_CorrCalc.m.......................................................... . (12)
5.1.3.4
Fun_ParameterCalc.m..................................................... (12)
5.1.3.5
Fun_ChannelCorrect.m ................................................... (13)
5.1.3.6
Fun_ZDR_KDP_Calc.m ..................................................... .. (13)
5.1.3.7
Fun_FormatChange.m ..................................................... .. (14)
5.1.3.8
Fun_Thresholding.m ..................................................... .. (14)
5.1.3.9
Fun_SpeckleFilter.m .................................................... ............ 14 5.2 第2个子过程的算法介绍...................................................................... ...... 15 5.3 第3个子过程的算法介绍(未完成).. (16)
5.3.1
Read_MSD_Radar_BaseData_PHIDP_Correct_ALL (16)
5.3.2
Fun_PHIDP_Initial_Modify ............................................... . (16)
5.3.3
Fun_PHIDP_Reject_Singular_Data ......................................... . (16)
5.3.4
Fun_PHIDP_Remove_Fold .................................................. .. (16)
5.3.5
Fun_PHIDP_Filter ....................................................... .. (16)
6 计算结果的初步对比(未完成).................................................. 16 6.1 对各个IOP 强度起伏的对
比 ..................................................................... .. 17 6.2 PPI的对
比 ..................................................................... (20)
6.2.1 层状云降水...................................................................... (20)
6.2.2 对流型降水......................................................................
.................. 22 6.3 RHI的对比...................................................................... . (24)
6.3.1 层状云降水...................................................................... (24)
6.3.2 对流型降水...................................................................... .................. 27 6.4 地物抑制效果的对比...................................................................... .. (29)
7 下一步工
作 ..................................................................... .............. 30 7.1 将用matlab重新计算出的基数据分发给相关的研究人员，终于可以开展下一步的科学研究了...................................................................... ..................... 30 7.2 按照在第2章中谈到的第2个阶段，把目前基于Matlab的GMAP和解距离模糊等算法，改用C语言来实
现 .............................................................. 30 7.3 召集人员，进行和RVP9的深入对比分析，包括两者的地物抑制效果、参数的随机差的对比，看看两者算法各自的优势和缺陷. (30)
1 概述
2014年6、7月份，由北京敏视达公司制造的南京大学C波段双偏振雷达在安徽长丰站进行了持续的降水观测，录取了大量的基数据和IQ数据。

但是，在对雷达自身的RVP9信号处理器给出的基数据的分析过程中，发现该基数据存在一些较大的问题，主要是:强度值和金属球对比有3dB的偏差、强度经常会突然发生最高达6dB的变化、地物抑制只有20多dB 等问题(见《用金属球对敏视达C波段双极化雷达标定结果.doc》、《南京大学C波段双偏振雷达_回波强度发生突然变化_分析报告.doc》和《南京大学C波段双偏振雷达_地物抑制效果比较差的问题_分析总结报告.doc》等文档)。

因此，我们决定采用自己的、基于Matlab编写的信号处理算法，对IQ数据进行全面的计算，重新生成基数据，以便于后续的科学研究。

于是，从2014年9月7日开始，一直到9月底才结束，对IQ数据进行了第一次的全面重新计算。

这次重新计算的基数据，经过和金属球标定、仪表测试、接收到的太阳能量等多方面的对比，我们认为其给出的dBZ、ZDR是比较准确的，dBZ的误差不超过
1dB，ZDR的误差不超过0.2dB，可以满足后续科学研究的要求了。

并且，还和位于长丰站北面30km处的一维雨滴谱仪进行了对比，也表明dBZ
和ZDR和雨滴谱计算的结果是高度一致的(相差约,dB)。

但需注意的是，有以下任务是不归这次计算负责的:
, 利用PHIDP订正ZDR、dBZ降水衰减的系数优化问题;
, 利用基数据中的各个参数，如:谱宽，识别和去除二次回波;
, 利用风场连续性识别和去除速度模糊;
, 利用基数据中的各个参数，如:CC、ZDR等，对非降水回波，包括:地
物、昆虫、鸟类等回波的识别和去除;
, 以及其它涉及到科学问题(也就是说，不是雷达的问题)上的各种算法，
如:所谓共振对PHIDP的影响、所谓共振对CC的影响等。

上述任务是应该由其他研究人员，在后续的科学研究中，负责完成的。

本文将计算过程中的混合编程方法、信号处理算法、经验教训、和RVP9结果的对比等进行详细介绍，为今后下一步的算法改进打下基础。

注意:由于时间关系，第5.3章和第6章尚未完成。

1
9月26日，补充了{'IOP1';'IOP2';'IOP4';'IOP7';'IOP10';'IOP11'}的强度起伏的结果。

10月13日，把其它几个IOP的结果页补充完整了。

2 整个计算步骤
考虑到我们目前的信号处理水平还在很低的层次，也不可能等到信号处理算法全面完善之后，再计算给出基数据(那样的话，后续的研究人员要等的花儿也谢了)。

因此，我们准备将整个的IQ重新计算的任务，分为3个阶段，分步实施: , 第1个阶段(就是本次9月份重新计算的过程):利用常规的FFT地物
抑制算法、MultiLag相关系数计算方法，但没有考虑解距离模糊的问题，将IQ数据先计算一遍，给出基数据给后续的研究人员使用。

, 第2个阶段:(准备在今年底、明年初开始计算):采用基于GMAP或
BGMAP的地物抑制算法，并考虑解距离模糊的问题，将IQ数据重新
计算一遍。

, 第3个阶段:(准备到明年6、7月份的观测试验前完成):采用新一代
的算法，主要是地物的频谱分析和识别、地物反演折射率等算法(还尚
未开始研究呢,)，并对天线解码延迟造成的角度误差进行补偿修正，
将IQ数据重新再度计算一遍，其结果将作为本次课题观测试验的最终
数据。

每个阶段又再分为3个子过程:
, 第1个子过程:对IQ数据进行地物抑制、相关计算、参数估计，得出双
偏振的各个参数，形成初步的基数据文件。

这个子过程耗时最长，计算
机需要连续计算20天左右(总共约20TB的IQ数据)。

, 第2个子过程:对上面的基数据进行垂直90度分析，给出其ZDR、PHIDP 的系统偏差随时间的变化情况，并绘制成曲线，写成文档记录下来(主
要是由于雷达硬件性能的变化造成的)。

这个子过程的计算需耗时2天
左右。

, 第3个子过程:利用分析出每一个体扫的ZDR、PHIDP系统偏差情况，
对基数据进行修正，形成最终的基数据结果。

同时，对PHIDP和KDP
进行剔除异常值、去折叠、自适应平滑滤波等处理，给出一个更高质量
的KDP。

这个子过程的计算需耗时3天左右。

2
在计算之前，已经将所有的IQ数据都拷贝至DS1813磁盘阵列中了(具备28TB 的容量)，该磁盘阵列的读写速度可以达到100MByte/s以上，满足了IQ计算的要求。

IQ数据保存在磁盘阵列的“NJU_CPol\IQData”目录下，用Matlab重新计算的基数据保存在“NJU_CPol\BaseData_DSP”目录下。

3 Matlab和C语言混合编程方法
3.1 概述
为了提高Matlab对IQ数据计算的速度，决定采用Matlab调用C语言编写的核心算法的方法，来提高速度。

在编程中，主要参考了《精通Matlab与C/C++混合程序设计(第3版)》，刘维编著，第80页~第180页的内容。

C程序在“H:\QXLD\NJU_CPol\Matlab\敏视达雷达IQ数据的信号处理\C” 目录下，其主程序是mexFun_Radar_Data_DSP.c，其它几个C文件是子函数。

程序是参照了Fun_Radar_Data_DSP.m 和以前写过的PDRadarDSP.cpp。

为了提高速度，程序大量采用了Intel 的IPP 算法库和32bit的浮点运算。

程序中的各种信号处理算法，以及其中的各种系数，是完全按照
Fun_Radar_Data_DSP.m 程序写的(但没有数据显示的功能，并且双重频处理、高阶的MultiLag算法等，暂时没有实现)。

为了确保C程序计算结果的正确性，在C程序编写完成之后，将各个体扫下、各个仰角层面、各个径向、各个距离单元的处理结果和基于Matlab程序的处理结果进行了仔细的对比，发现两者是完全一致的(由于C语言是采用32bit浮点计算的，而matlab是采用64bit浮点计算的，因此会有百万分之一以下的微小差别，这是完全可以忽略的)。

考虑到充分发挥Matlab和C语言各自的优势，因此并没有把全部的Matlab程序都改用C语言来实现，而仅仅把耗时最长的Fun_Radar_Data_DSP.m用C语言实现了。

今后，可以将GMAP等先进的算法也用C语言来实现，因为这些算法用Matlab 编写时，难以实现“向量化计算”，造成计算速度很慢。

而用C语言编写就不存在“向量化”的难题，计算速度应该很快。

3.2 编程方法
3.2.1 编译器的设置
3
先在matlab下执行:
mex -setup
进行C语言编译器的设定。

注意要选择VC2010作为编译器(不要采用Matlab自带的编译器)。

并且要注意:如果要把编译好的程序在没有安装VC2010的电脑上运行，则要安装VC2010的可发行包。

然后，运行:
mex C\mexFun_Radar_Data_DSP.c
就可以将 C\ 目录下的“mexFun_Radar_Data_DSP.c”文件，以及包含的各个子函数文件，编译为“mexFun_Radar_Data_DSP.mexw32”，就可以被matlab直接调用了。

如果采用了-g选项(可以进行断点调试)，即运行:
mex C\mexFun_Radar_Data_DSP.c -g
就不但可以生成 mexFun_Radar_Data_DSP.mexw32文件，还能生成
“mexFun_Radar_Data_DSP.mexw32.pdb” 文件，用于断点调试。

注意:如果在Win7 64bit下运行，则需要重新编译一下，生成
“mexFun_Radar_Data_DSP.mexw64”，才能在64bit 的matlab下调用。

3.2.2 Intel公司IPP算法库的使用
需要先安装IPP的开发包，在“Y:\SETUP\DevTools\各种软件开发工具包
\Intel算法库\IPP6.1” 目录下，安装到C盘的默认目录下。

需要用到intel_TBE (2).lic文件进行安装。

在C程序的开头，需要加入:
#include "C:\\Program
Files\\Intel\\IPP\\6.1.2.041\\ia32\\include\\ipp.h"
#pragma comment( lib,"C:\\Program
Files\\Intel\\IPP\\6.1.2.041\\ia32\\lib\\ippcorel.lib")
#pragma comment( lib,"C:\\Program
Files\\Intel\\IPP\\6.1.2.041\\ia32\\lib\\ippsemerged.lib") #pragma comment( lib,"C:\\Program
Files\\Intel\\IPP\\6.1.2.041\\ia32\\lib\\ippsmerged.lib") 注意:经测试发现，使用支持多线程并行计算的静态库，没有加速效果(实际运行时，仍然是单线程)。

因此，就决定只采用单线程的库了。

具体使用方法详见“C:\Program Files\Intel\IPP\6.1.2.041\ia32\doc”目录下的《userguide_win_ia32.pdf》。

3.2.3 IPP在64bit 操作系统下的编译
4
注意:在64bit操作系统下，需要安装IPP 6.0的64bit版本(暂时没有6.1的64bit 版本)。

然后，在C程序的开头，需要加入:
#include "C:\\Program
Files\\Intel\\IPP\\6.0.0.062\\em64t\\include\\ipp.h"
#pragma comment( lib,"C:\\Program
Files\\Intel\\IPP\\6.0.0.062\\em64t\\lib\\ippcoreem64tl.lib") #pragma comment( lib,"C:\\Program
Files\\Intel\\IPP\\6.0.0.062\\em64t\\lib\\ippsmergedem64t.lib") #pragma comment( lib,"C:\\Program
Files\\Intel\\IPP\\6.0.0.062\\em64t\\lib\\ippsemergedem64t.lib") 具体使用方法详见“doc”目录下的《userguide_win_ em64t.pdf》。

3.3 速度测试
3.3.1 计算机配置
计算机的配置如下:
, CPU:I5
, 主频:3.2GHz
, 内存:4GB
, 操作系统:Win7 32bit
3.3.2 处理一个体扫的IQ数据的耗时
基数据在“H:\QXLD\NJU_CPol\BaseData_Select\用于C程序测试的数据_20140813”目录下，基数据为:NJU.20140704.203525.AR2。

对应的IQ数据在“H:\QXLD\NJU_CPol\IQData\用于C程序测试的数据_20140813”目录下，共15个IQ文件。

结果对比如下:
全部耗时扣除了IQ数据
读取的耗时
560.25s 408.81s Matlab的
Fun_Radar_Data_DSP 函
数的处理时间
277.75 126.31 用C语言编写的
mex_Fun_Radar_Data_DSP
函数的处理时间(采用-g
编译)
256.33s 104.89 用C语言编写的
mex_Fun_Radar_Data_DSP
函数的处理时间(不采用-g
编译)
5
从上表可以看出，采用C语言编写的函数，其执行时间仅为Matlab耗时的
1/3~1/4(扣除了IQ数据读取的耗时)。

在实际观测中，天线一个体扫模式(共15个仰角层面)，共耗费400s左右。

因此，从表中可以看出，采用C语言对Matlab程序进行优化之后，就能实时实现各种处理算法了。

3.3.3 mexFun_Radar_Data_DSP 函数中，各个子函数的速度
测试数据为64个相关脉冲(即方位上的积累数)，距离库为1933个，库长为75m。

测试方法是通过1000次循环，然后平均得到的。

函数功能耗时(采用-g 耗时(不采用
编译)(ms) -g 编译)(ms)
1.79 1.00 只有数据的输入和输出，即只执行了
Get_DSP_Setup、Get_CorrelationData、
Set_RayRecord 3个函数，用于和Matlab空
间的数据交互
4.05 Fun_PreProcess.c 函数 4.48
10.61 8.27 Fun_ClutterFilter.c函数
Fun_CorrCalc 4.23 4.14 Fun_ParameterCalc 1.18 0.64
Fun_ChannelCorrect 0.1 0.1 Fun_ZDR_KDP_Calc 1.16 0.73 Fun_FormatChange 0.32 0.31 Fun_Thresholding 0.1 0.1 Fun_SpeckleFilter 0.18 0.34
从上表看出，Fun_ClutterFilter.c(对地物进行抑制)的计算最耗费时间。

这是正常的，因为该函数涉及到频谱的变换，运算量最大。

另外，经测试，对于同样的IQ数据，Matlab的Fun_Radar_Data_DSP 函数的执行时间为:99.86ms，是C语言耗时的3~4倍。

3.3.4 多个Matlab进程同时运行
经测试，多个Matlab同时运行，全局变量不会在不同进程间共享，从而不会对计算结果造成影响。

也就是说，可以同时运行几个matlab进程进行计算，分别计算不同目录下的IQ数据。

6
在正式的计算中，采用了配备有12GB内存的I7 四核计算机，主频为
3.4GHz，安装了Windows7 64bit，同时运行了4个matlab进行计算，将计算速度提高了3倍左右。

需要注意:如果硬盘的读写速度比较慢，就会造成IQ数据读取需要耗费大量的时间，就体现不出多个Matlab同时运行的好处了。

3.4 经验和教训
在编程中，有以下一些经验和教训要吸取:
, Matlab提供的(int nlhs, mxArray *plhs[],int nrhs, const mxArray
*prhs[])几
个参数，用于和C语言进行数据交互，非常方便，非常好用，需要彻底
掌握该编程方法。

参见mxArrayToStruct.c程序。

, C语言中，所有的变量要在函数的一开头进行定义，而不能像C++ 语言
一样(只要在用之前定义即可)。

, 占用内存比较大的变量(如保存IQ数据的变量)，必须要定义成全局变
量。

否则，编译会报错(堆栈的空间不够)。

, 注意:要用fabs，而不能用abs! 因为C语言中的abs是整数取绝对值～, 经过Matlab和IPP中的窗函数对比发现，两者有点差别。

因此，这里要将matlab的窗函数用公式写出来，而不能直接采用matlab提供的，如
hamming、hanning的函数等，以便和IPP中的完全一致
, 在反异步的程序中，原来的C语言代码有缺陷。

当判断出当前距离单元
是异步干扰信号时，不能将Abs的结果清零，而只能先将IQ清零。

等
到反异步结束之后，再从IQ重新计算出 Abs的结果。

4 IQ计算时的参数设置
4.1 雷达的性能参数
Matlab信号处理中的雷达性能参数，是根据2014年5月27日，对雷达进行了仪表测试的结果而设定。

测试结果详见《南大C波段雷达标校情况_2014年5月30日.doc》。

参数名称值备注
-28.31dBz 水平通道的dBZ0 这两个结果是根据5月27日用仪表
-28.58dBz 测试的结果，计算出来的垂直通道的dBZ0
-75.96dBm 水平通道的噪声功率接收机输入端接噪声源，但不通电时
-76.40 dBm 测得垂直通道的噪声功率
7
0 deg 两个接收通道的相位差 2014年5月24日，由于对波导进行
了拆装，因此初相发生了变化。

于是
在数字接收机中，对初相进行了调整
(大约是5月24日UTC时间4点左
右)。

调整之后，ReceiverPhaseDiff=0
下面一些参数，是从RVP9给出的IQ数据中，读取出来的(和雷达的工作模式有关):
参数名称值备注
PolarizationMode 3 雷达的极化工作方式 h=0,v=1 ,hv=3 PulseWidth 5e-07 发射脉冲宽度，单位:s CarrierFrequency 5.6235e9 雷达中心载频，单位:Hz
注意:实际的载频为5.625e9 TrigDataLength 1933 每个触发下的距离库的数目 CorrelationNumber 364 天线一圈的径向数
PRT1 0.001 每个触发的重复间隔，单位:s SamplePeriod 5e-07 每个距离库的库长，单位:s
4.2 信号处理算法设定的参数
以下是Fun_Read_MSD_Radar_IQData_And_DSP.m函数中，设定信号处理算法所需的一些门限参数:
参数名称值备注
InterferenceFilterGate 2 在反异步(抗同频干扰)的算法中，
相邻周期的幅度相差多少倍，就认
为是异步干扰信号 ClutterFilterSel 根据仰角不对杂波滤除的滤波器的凹口的大
同而设定小设置
仰角<1度，则ClutterFilterSel=4;
1度?仰角<2度，则
ClutterFilterSel=3;
2度?仰角<3度，则
ClutterFilterSel=2;
其它仰角则ClutterFilterSel=1
8
KDP_Average_Distance 2000m 定义KDP计算中，参与平均计算
的距离范围
ThresholdLOG -2dB 信噪比门限，单位:dB。

经过和
RVP9给的结果对比，认为-2比较
合适(另外，由于高仰角(>5度)
的噪声会低1.5dB，因此对于高仰
角回波而言，相当于门限为
-0.5dB)
ThresholdCCOR -125.0dB 杂波对消的门限注意:是负数。

注意:设定这么小，就是为了让该
功能无效
USE_MULTILAG_DBZ false dBZ的计算不用Multi-Lag USE_MULTILAG_W false W的计算也不用Multi-Lag USE_MULTILAG_CORR true 相关系数采用Multi-Lag的方法计
算
USE_MULTILAG_NUMBER 1 采用Multi-Lag中的 Lag1的方法
计算
ClutterFilter_USE_GMAP 0 暂时不采用GMAP算法来进行地
物的滤除
5 信号处理算法介绍
5.1 第1个子过程的算法介绍
第1个子过程是指:对IQ数据进行地物抑制、相关计算、参数估计，得出双偏振的各个参数，形成初步的基数据文件。

主要分为以下各个函数:
5.1.1 Read_MSD_Radar_IQData_And_Process_Main.m
这是整个计算的主程序，其计算流程如下:
1. 首先，设定要进行计算的文件的目录，包括:RVP9给出的基数据的目
录、IQ数据的目录、保存计算结果的目录;
2. 然后，遍历RVP9给出的基数据目录下的每一个文件。

检查该基数据文
件是否已经重新用信号处理算法，计算过IQ数据了。

如果已经生成了，
则根据用户设定，决定是否重新生成;
9
3. 从RVP9给出的基数据文件中，读取出需要的数据结构，如Header、Site、
Task、Cut等数据，赋值给BaseData_DSP(该结构中保存的就是信号处
理的结果)
4. 然后，遍历所有的层，查出这个仰角层面的时刻所对应的IQ文件名。

注
意:由于基数据中的时间信息和IQ数据的文件名，可能会存在一点点
的差异，所以需要检查前后5秒钟内，有无对应的IQ数据文件; 5. 找到对应的IQ文件之后，读取记录在该文件中的雷达特性，对dBZ0、
dBZ0_2、RecieverNoise_dBm、RecieverNoise_dBm_2、ReceiverPhaseDiff 等5个参数进行修改(这5个参数根据5月27日用仪表标定的结果而
设定);
6. 调用信号处理函数Fun_Read_MSD_Radar_IQData_And_DSP，处理该IQ
文件;
7. 将信号处理函数生成的结果(已经保存为mat文件)，读取到
BaseData_DSP 结构中;
8. 不断循环，遍历每个仰角层面;
9. 将BaseData_DSP结构，保存至计算结果的基数据的目录中。

注意:这
里比RVP9 输出的多一个“相位”的数据，可用于研究大气折射率的变
化。

10. 不断循环，遍历目录下的每个文件。

5.1.2 Fun_Read_MSD_Radar_IQData_And_DSP.m
该函数的计算流程如下:
1. 设定信号处理算法所需的一些门限参数;
2. 预先分配空间，以加快程序运行速度;
3. 遍历每个相关脉冲组的IQ数据;
4. 调用信号处理函数mexFun_Radar_Data_DSP(C语言版)，或者是
Fun_Radar_Data_DSP(Matlab语言版)，进行计算;
5. 根据需要，对比C语言版和Matlab版两种的计算结果;
6. 根据需要，可以对某一个距离库进行深入分析，绘制其FFT结果;
7. 最后，将结果保存成MAT 文件。

5.1.3 Fun_Radar_Data_DSP.m 该函数的计算流程如下(依次调用9个计算子函数):
1. 调用Fun_PreProcess子函数，对输入的数据进行反异步处理;
10
2. 调用Fun_ClutterFilter子函数，对地物杂波进行滤波，采用FFT变换至
频域，滤除零速附近的通道，再IFFT至时域的方法;
3. 调用Fun_CorrCalc子函数，对相关脉冲的数据进行相关计算，得出T0、
R0、T1、R1等;
4. 调用Fun_ParameterCalc子函数，对强度、速度、谱宽等各种参数进行计
算;
5. 调用Fun_ChannelCorrect子函数，对两个通道的幅度和相位进行补偿;
6. 调用Fun_ZDR_KDP_Calc，进行ZDR和KDP的计算;
7. 调用Fun_FormatChange子函数，将结果转换成特定的格式输出;
8. 调用Fun_Thresholding子函数，对数据进行门限控制;
9. 最后，调用Fun_SpeckleFilter 子函数，进行异常点滤除。

5.1.3.1 Fun_PreProcess.m
该子函数的功能是:反异步(即抗同频异步雷达干扰)。

计算过程如下: 1. 首先对第1个脉冲进行反异步处理，计算该脉冲和第2个脉冲的幅度差; 2. 对幅度差进行判断，若超过门限(InterferenceFilterGate与当前回波幅度相乘+噪声电压*3)，则认为是异步干扰
3. 将是异步干扰信号的距离单元的IQ值置为0
4. 然后，对其余相关脉冲，与前1个比较，进行反异步处理
5. 双偏振方式，则对另一个通道的数据也要进行处理
5.1.3.2 Fun_ClutterFilter.m 该子函数的功能是:对地物杂波进行滤波。

采用FFT变换至频域，滤除零
速附近的通道，再IFFT至时域的方法。

计算过程如下:
1. 对位于同一个距离库、相邻触发下的IQ数据进行加权处理(目的是为了
抑制FFT变换的泄露)。

经过试验，这里选用hanning窗，效果比较好。

注意:经过对比Matlab和IPP算法库中的窗函数发现，两者有点细微的
差别。

因此，需要将matlab的窗函数用公式写出来，以便和IPP中的完
全一致。

2. 对窗函数加权之后的IQ数据进行FFT变换;
3. 根据不同的仰角，选择不同的处理算法。

如果仰角<1度，则
ClutterFilterSel=4;1度?仰角<2度，则ClutterFilterSel=3;2度?仰角
<3度，则ClutterFilterSel=2;其它仰角则ClutterFilterSel=1;
11
4. 根据ClutterFilterSel的值，去除零速及附近几个通道的值，这些通道的
值采用相邻通道的值来代替;
5. 进行IFFT，变换回时域;
6. 由于进行了窗函数加权，因此要乘以修正系数，以便补偿幅度的加权损
失;
7. 最后，如果是双偏振方式，则对另一个通道的数据也要进行处理。

5.1.3.3 Fun_CorrCalc.m
功能:对相关脉冲的数据进行相关计算，得出T0、R0、T1、R1等计算过程如下:
1. 首先用abs、mean、angle等函数，计算出T0、Phase、R0;
2. R1(经过地物对消后的IQ的一阶相关)、T1(未经过地物对消后的IQ
的一阶相关)采用for循环，即可计算出来;
3. 如果是双偏振方式，则对另一个通道的数据也要进行处理;
4. 然后，采用for循环，对经过地物对消后的IQ进行水平、垂直之间的相
关计算，得到R0*conj(R0_2);对未经过地物对消后的IQ进行水平、垂
直之间的互相关计算，得到T0*conj(T0_2);
5. 根据Multi-Lag算法的要求，增加计算fcT0_MUL_CONJ_T0_2_LAG1、
fcR0_MUL_CONJ_R0_2_LAG1等参数;
6. 根据需要，可以绘制特定距离上、不同LAG的结果，以便深入研究
Multi-Lag的原理和效果。

5.1.3.4 Fun_ParameterCalc.m 功能:对强度、速度、谱宽等各种参数进行计算
计算过程如下:
一、计算强度
, 计算20*log10(r/r0)，r0=1km，进行距离订正
, 计算大气损耗a*r修正根据雷达手册第2版第2章图2.22 ~ 2.23 , 计算信噪比，10*log10((T0-N)/N)，采用RVP8 用户手册第5.2.7章的算法
, 计算杂波修正CCOR =10*log10(R0/T0) , 计算dBZ=10*log10((T0-N)/N)
+dBZ0+ 20*log10(r/r0) + a*r +
10*log10(R0/T0)
dBT=10*log10((T0-N)/N) +dBZ0+ 20*log10(r/r0) + a*r;
12
二、计算速度(这里只考虑单重频的方式，因为在观测试验中，没有采用DPRF 模式)
, 计算最大速度fVMax =fRadarWavelength/DSP_Setup.PRT1/4.0
, 计算R1的相位/PI* fVMax，得出速度
三、计算谱宽，。

, 计算SQI，采用RVP8 用户手册第5.2.9章的算法，用R0、R1计算。

, 采用Lag0和Lag1，即标准方法来计算谱宽。

注意:要将T0 减去噪声
才行。

注意:这里的 Lag1是采用的未经过地物抑制的IQ数据来进行计算的(还是应该用经过地物抑制的IQ数据来计算?)。

四、计算双极化的一些参数
, 计算LDR，即垂直和水平通道的dBZ之差;
, 计算PHIDP，即fcR0_MUL_CONJ_R0_2的角度
, 根据“A Multi-lag Correlation Estimator for Polarimetric Radar Variables in
the Presence of Noise .pdf”论文中的公式，采用Multi-Lag中的Lag1公
式来进行相关系数的计算(因为在前一阶段的研究中发现，采用Lag1
的结果比较可信，而采用更多阶数计算的结果存在问题，特别是当回波
频谱比较宽的时候)。

计算公式如下:
注意:这里的 Lag1是采用的未经过地物抑制的IQ数据来进行计算的(还是应该用经过地物抑制的IQ数据来计算?)。

5.1.3.5 Fun_ChannelCorrect.m
功能:对两个通道的幅度和相位进行补偿。

代码很简单，如下:
m_pstResultDataStruct.fDiffPhase = m_pstResultDataStruct.fDiffPhase +
DSP_Setup.ReceiverPhaseDiff /180.0*pi;
m_pstResultDataStruct.fDiffPhase=mod(m_pstResultDataStruct.fDiffPhas e,2*pi
);
5.1.3.6 Fun_ZDR_KDP_Calc.m
功能:进行ZDR和KDP的计算
13
计算过程如下:
1. 通过水平和垂直通道的dBZ之差，得到ZDR;
2. 然后，遍历每个距离库，计算KDP。

首先，计算出要用多少个距离库的
PHIDP进行计算(根据KDP_Average_Distance而设定，这里采用前后
2km的距离库);
3. 对该距离单元前后2km内的PHIDP值进行挑选，根据门限(如果
SQI<0.5，则不选择该组数据)，挑选比较好的PHIDP的值; 4. 如果有效数据太少，则该距离库的KDP无效;
5. 以第一个数据的相位为基准，后续的数据依次与其相比较，对相位进行
修正，防止相位出现折叠;
6. 最后，用最小二乘法，计算相位的斜率，并将单位换算成度/km。

注意:
要将斜率/2，因为KDP是双程定义的。

5.1.3.7 Fun_FormatChange.m
功能:将结果转换成特定的格式输出。

将m_pstResultDataStruct中的参数拷贝给m_stRayRecord，并将PHIDP转换成以度为单位(原来是以弧度为单位)。

5.1.3.8 Fun_Thresholding.m
功能:对数据进行门限控制。

计算过程如下:
1. 遍历每个距离库，对每个库单独进行门限的控制;
2. 利用LOG(即SNR)和CCOR(即地物杂波对消的门限)，对dBZ进
行判断，如果小于门限，则置为无效值;
3. 如果dBZ为无效值，则ZDR或LDR(单发双收模式下)也置为无效值;
4. 利用LOG，对dBT、V、W进行判断，如果小于门限，则置为无效值;
5. 利用LOG，对双偏振参数CC、PHIDP、KDP进行判断，如果小于门限，
则置为无效值;
6. 利用LOG，对回波的Phase(该参数可用于研究大气折射率的变化)进
行判断，如果小于门限，则置为无效值;
7. 注意:SNR和SQI两个结果，不受门限的控制。

5.1.3.9 Fun_SpeckleFilter.m 功能:进行异常点滤除，采用RVP8 用户手册第5.3.3章的算法。

计算过程如下:
14
1. 若当前距离单元的值无效，则若前后1个距离单元的值都有效，则将2
个值平均后作为当前距离单元的值;
2. 若当前距离单元的值是有效的，则若前后1个距离单元的值都无效，则
将当前距离单元的值也设为无效;
3. 注意:第一个距离点和最后一个点不进行异常点滤除
5.2 第2个子过程的算法介绍
第2个子过程是指:对第1个子过程给出的基数据进行垂直90度分析，给
出其ZDR、PHIDP的系统偏差随时间的变化情况，并绘制成曲线，写成文档记录下来(系统偏差的变化主要是由于雷达硬件性能的变化造成的)。

详见《南京大学C波段双偏振雷达_垂直90度观测_分析总结报告.doc》。

Matlab程序为:
, Read_NJU_CPol_Radar_BaseData_90deg_Analysis.m
, Read_NJU_CPol_Radar_90deg_Analysis_And_Display.m
主要分析过程如下:
1. 从雷达记录的基数据文件中，读取各个数据，如果是RHI模式，则不进
行后续的计算;如果该文件是单发双收模式，该文件也不进行分析;
接下来要判断90度是否有合适的降水，要有合适的降水，才适合进行下2.
面的各种分析。

判断的准则是:选取零度层以下，且不含地物的相关系
数(CC)的数据来进行判断。

选取距离(即高度)从1000m~3000m的
CC值，总共有约2500个数据。

如果在这些数据中，有超过5%的CC
小于0.95(包括没有回波的数据)，则说明该文件不是降水回波，该文
件就不再进行后续的分析了;
3. 先对ZDR进行分析。

分析从0~6000m不同高度上的分布，各个参数的
随机误差和系统偏差，再绘制ZDR随方位变化的波形图(选取
1000~2000m的高度)、直方图、PPI图;
4. 再对PHIDP随方位的变化进行分析，绘制垂直90度扫描时PHIDP随方
位变化的波形图(选取1000~2000m的高度)、PPI图。

5. 将结果保存为TXT文件和mat文件。

6. 最后，从90度分析结果的文件中，读取雷达的各个体扫的分析结果，选
取距离(即高度)为2500m的数据，作为该体扫下，雷达自身的ZDR、
PHIDP的系统偏差(选取2500m的距离，是因为经观察回波发现，在
这个距离处，地物的回波比较弱)。

15
另外，编写了mean_PHIDP.m和std_PHIDP.m两个函数，用于对PHIDP进行平均和方差的计算。

注意:对于PHIDP，不能简单的用Matlab提供的mean函数进行平均、用std函数进行方差计算，否则对于接近180度和-180度附近的数据，就会出现严重的错误～
5.3 第3个子过程的算法介绍(未完成)
第3个子过程是指:利用分析出每一个体扫的ZDR、PHIDP系统偏差情况，对基数据进行修正，形成最终的基数据结果。

同时，对PHIDP和KDP进行剔除异常值、去折叠、自适应平滑滤波等处理，给出一个更高质量的KDP。

参见《对敏视达雷达回波进行基于PHIDP的dBZ和ZDR订正_2014年4月5日,18日.doc》。

主要分为以下各个函数:。