雷达回波外推自学笔记，希望大家多多指导，谢谢（待续）

小莹子 · 发表于 2019-11-28 19:58:25

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

本帖最后由小莹子于 2019-12-27 11:24 编辑

   国内天气雷达经历了模拟天气雷达、数字化天气雷达、多普勒天气雷达三个阶段，并在美国WSR-88D的基上上，研发了新一代天气预警雷达CINRAD。多普勒天气雷达基于多普勒效应，即物体的频率会随运动发生变化。云、雨等散射体与雷达进行相对运动是，会形成多普勒效应。多普勒天气雷达不仅能探测散射体的形状、位置和分布情况，还能探测散射体内部质点的运动状态。多普勒天气雷达已经成为气象业务应用中最重要、最常用的雷达设备，新一代多普勒天气雷达发射脉冲形式的电磁波，然后雷达接受散射回来的能量，并根据接受回波的返回时间、来向、强度、相位、偏振等信息，确定气象米表的空间位置和性质。多普勒天气雷达除了测量回波强度之外，还可以用于监测和警戒强对流天气、灾害性天气等等。传统天气雷达只能测量回波强度（即反射率因子）。    我国新一代天气雷达网包含了S波段和C波段共8个型号的雷达。S波段波长为10cm（波长7.5cm~15cm），具有较好的穿透力，多布设在沿海和多强降水地区（湿润多雨、云层较厚），主要用于监测和预警灾害性大风、短时暴洪等强对流天气。C波段波长为5cm（波长3.75cm~7.5cm）,具有较高的分辨率，主要分布在内陆地区（干燥少雨、云层较薄），用于监测局部暴雨、冰雹等灾害性天气等。S波段包括SA、SB、SC型雷达；C波段包括CB、CC、CD、CCJ型雷达。雷达扫描方式包括PPI、RHI、CAPPI、VCS。
一、基数据结构
   天气雷达的基数据以二进制文本进行存储，主要由径向的运动速度（径向速度）、回波强度、速度谱宽（速度脉动程度的度量）等数据组成。基数据分为公共数据块和径向数据块两部分，公共数据块用于提供数据站点信息、配置任务等公共信息；径向数据块用于存储天气雷达的探测资料，包括径向头、径向数据头以及径向数据。
   基数据文件格式主要有SAB、CB、88D、SCD、CC1.0、CC2.0。其中SAB格式基数据文件，无文件头，使用固定长度径向数据（2432个字节），VCP21扫描层数9层。CB格式基数据文件，无文件头，使用固定长度径向数据（4132个字节），按主机字节库存储，VCP21扫描层数11层。88D格式文件头（16个字节），基数据文件使用固定长度径向数据（2432个字节），VCP21扫描层数9层。SCD格式基数据文件头（1024个字节）+固定长度径向数据（4000个字节），VCP扫描层数9层。CC1.0格式基数据文件使用文件头（1024个字节），固定长度径向数据（3000个字节），按主机字节序存储。CC2.0格式基数据文件使用文件头（2060个字节），径向数据是变长的，VCP21扫描层数9层。
   基数据：Z_RADR_I_Z9230_20140417155600_0_DOR_SA_CAP.bin
文件命名规则：Z：固定编码，表示国内交换文件；RADR：固定编码，表示雷达资料；I_Z9230表示雷达区站号；yyymmddhhMMss为体扫结束后文件生成时间（年、月、日、时、分、秒，用UTC世界时间表示）；O，固定编码，表示观测资料；DOR，表示多普勒雷达；Bin，表示二进制文件。
   反射率因子随时间的变化常被用于预测降水回波的变化趋势，回波强度变化，强降雨（雪）带，风暴的结构和强弱等信息。反射率因子Z的定义是单位体积中降水粒子直径的6次方总和，常用单位为mm6/m3，通常用dBZ表示;常用来表示气象目标的强度，反映了降水粒子的尺度和疏密度。
   基本反射率因子（R）用途：估计系统的强度，冰雹的潜在性，风暴结构；平均径向速度产品（V）用途，估计风向、风速，确定水平和垂直切变（不连续边界），识别强天气系统；基本径向速度谱宽产品（SW）用途：提供由于风切变、湍流或样本至两引起的径向速度变化幅度的度量，可用来帮助确定边界位置、TBSS估计湍流大小及检查径向速度是否可靠；组合反照率因子（CR）揭示了一个体积扫中的最高反射率因子的分布。（摘自传说中的谁，CINRAD-SA雷达系统主要产品与算法）
二、CINRAD雷达产品显示系统读取雷达基数据
下载安装操作参考：
“CINRAD雷达产品显示系统读取雷达基数据”http://bbs.06climate.com/forum.p ... hlight=%C0%D7%B4%EF”
“与大家分享CINRAD雷达产品显示系统http://bbs.06climate.com/forum.p ... %C0%D7%B4%EF&page=1”
1.安装成功并打开后，点【设置】-【站点信息】-【数据】，选择自己的文件格式如SA。
2.点【变量】，里面的基本产品包括“反射率、速度和谱宽”，这里以反射率为例；
3.然后在界面右侧有所要打开的文件信息，其中可以选择雷达的扫描层数即仰角，图中上侧可以看到雷达扫描时间，对应的仰角等信息；
4.如果需要将数据导出进行其他应用，点【文件】-【导出】-【当前产品】，会保存成txt文件，这时点开数据文件会看到方位角、库数和显示产品。
5.【分析】中可以选择RHI距离高度显示或VCS任意垂直剖面，任选两点即可。
注：如若更换了站点，需改变地图。地图文件名为province.map，位于CINRAD雷达产品显示系统安装目录（如：C:\Program Files\Yuhan\Radar）。程序附带的是北京市地图文件，如果要显示其他地区的地图，您需要使用该地区的map文件替换掉程序自带的文件（map文件可以在PUP软件中找到）。将您的地图文件命名为province.map，覆盖系统安装目录的province.map文件即可。（如果替换地图文件后，显示的地图没有即时更新，请点击工具栏上的"复位"按钮）
三、Matlab读取雷达基数据
代码参考又是那只猫“matlab读取SA雷达基数据文件并进行显示”http://bbs.06climate.com/forum.php?mod=viewthread&tid=14628
我自己添加了一些注释，便于理解，如有错误和理解不足之处，还请大家指教！
雷达基数据格式（SA、SB、CD）在论坛或是百度上都可以搜索到
四、雷达回波外推
   雷达回波外推技术是进行临近天气预报的主要技术手段。雷达回波外推是指根据天气雷达探测到的回波数据，确定回波的强度分布以及回波体（如风暴单体、降水区等）的移动速度和方向，通过对回波体进行线性或者非线性的外推，预报一定时间段后的雷达回波状态。
   雷达回波外推目前使用最广泛的传统外推算法为单体质心法和交叉相关法。光流法是从计算机视觉领域引入的算法，在预报变化比较快的强对流天气过程中更显优势。BP神经网络模型利用了神经网络的自主学习能力，可以获得回波变化的一些隐含特征，展现了用神经网络模型进行回波外推的潜力。
   单体质心法将目标简化为一点，适合对孤立的、大而强的回波单体或单体群进行跟踪和预报，当回波比较零散或者回波发生合并、分裂时，预报准确率大大降低，因此不能用于层状云降水系统。
   交叉相关法使用相邻两个时刻的空间优化相关系数建立最佳拟合关系，可以跟踪层状云降雨系统，但对于回波变化很快的强对流天气过程，难以保证追踪的准确性，效果会明显降低。

参考文献：

[1]滕志伟. 基于深度学习的多普勒气象雷达回波外推算法研究[D].湖南师范大学,2017.

[2]梁启君,梁军,黄骞.雷达回波外推算法的并行化及实现[J].地理空间信息,2013,11(06):47-50+9.

[3]曾小团,梁巧倩,农孟松,冯业荣,许向春,陈业国.交叉相关算法在强对流天气临近预报中的应用[J].气象,2010,36(01):31-40.

[4]朱平,李生辰,肖建设,徐亮,靳世强.天气雷达回波外推技术应用研究[J].气象,2008(07):3-9+129-132.

[5]乔春贵,郑世林,杨立志,霍锐,李周,鲁坦.质心法雷达回波外推的原理及应用[J].河南气象,2006(03):29-30.

爱到花开 · 发表于 2019-12-3 09:21:39

下一步呢?准备采用什么方法进行外推。另外，我看现在的做法都是用的BR或CREF外推，这种不知道合不合理，因为好多外推的限定条件都是在二维平面上。

小莹子 · 发表于 2019-12-3 09:30:23

本帖最后由小莹子于 2019-12-3 09:45 编辑

爱到花开发表于 2019-12-3 09:21
下一步呢?准备采用什么方法进行外推。另外，我看现在的做法都是用的BR或CREF外推，这种不知道合不合理，因 ...

我是自学，也不是气象背景，还没有学习到这一步，其实也是不知道往下怎么走。正在学习它的存储格式，并尝试用用matlab将其成功读取，只是还没有成功。还请多多指教，谢谢

爱到花开 · 发表于 2019-12-3 15:46:19

小莹子发表于 2019-12-3 09:30
我是自学，也不是气象背景，还没有学习到这一步，其实也是不知道往下怎么走。正在学习它的存储格式，并尝 ...

好吧，实际上我也仅是做到第二步，算出了风场，不知道怎么进行外推

小莹子 · 发表于 2019-12-3 15:50:33

爱到花开发表于 2019-12-3 15:46
好吧，实际上我也仅是做到第二步，算出了风场，不知道怎么进行外推

你好，那你第一步是怎么做的呀？用CINRAD还是Matlab？还是什么的呀？我现在正在装PUP，谢谢

爱到花开 · 发表于 2019-12-3 16:05:27

小莹子发表于 2019-12-3 15:50
你好，那你第一步是怎么做的呀？用CINRAD还是Matlab？还是什么的呀？我现在正在装PUP，谢谢

我用的是雷达CREF或者CAPPI数据，CREF就用Mosaic那个c程序就行，读出来是6201*4201的格点，筛选区域即可；CAPPI文件中我记得好像就是一般二进制文件吧，前14个字节是文件头

小莹子 · 发表于 2019-12-3 16:17:54

爱到花开发表于 2019-12-3 16:05
我用的是雷达CREF或者CAPPI数据，CREF就用Mosaic那个c程序就行，读出来是6201*4201的格点，筛选区域即可 ...

我的是SA文件，然后我用matlab读取，在网上找的代码，不过人家没有太全乎的注释，所以我图出来了，但不明白中间过程啥意思。我看论文上是主要是质心法和交叉相关法。

爱到花开 · 发表于 2019-12-4 08:05:17

/////////////////////////////////////////////////////////////
// //
// 读取全部雷达基数据最终形式，并实用 //
// //
/////////////////////////////////////////////////////////////
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//#define RGates 460 //S BAND 反射率距离库数
//#define VGates 920 //S BAND 速度距离库数
//#define WGates 920 //S BAND 谱宽距离库数
void pause()
{ char c;
printf("\nPlease enter any key; then go on...\n");
c=getchar();
}
int main(void)
{
//雷达信息头
//第一段定义传输信息类型说明，共14类
char Unused1[14]; //1~14未用文件数据变量的定义
// char DigitalRadarData； //1数字化雷达数据（综合了反射率，径向风速或速度谱宽）
// char RadarDataCollectionStatus //2RDA雷达数据采集状态数据
// char DisplayOrMaintenanceData //3显示/维护数据
// char ConsoleSignal //4控制台信号 RDA—RPG雷达产品生成
// char MaintenanceData //5维护记录数据
// char ControlInstruction //6RDA控制指令
// char VolumeScanMode //7体扫模式
// char ErrorCorrection //8纠错
// char DataRequest //9数据请求
// char ConsoleInformation //10控制台信息 RPG—RDA
// char LoopTestRDA //11环路测试 RDA—RPG
// char LoopTestRPG //12环路测试 RPG-RDA
// char ErrorMessageFilter //13错误信息过滤分流 RPG-RDA
// char EditedErrorMessageFilter //14编辑后的错误信息分流 RPG-RDA
//第二段定义报头数据项名称及内容描述
// unsigned short int MessageSize; //12~13数据报长度，表示记录长度，双字节单位表示
// char ChannelID; //14 RDA冗余通道，RDA信号传输通道号(=0,1,2)，对冗余通道配置
unsigned short int MessageType; //15~16报文类型记录数据类型,具体说明见6.1
// char Unused2[12]; //17~28未用
short int IDSequence; //17~18报文顺序号记录顺序号，0~0x7FFF之间循环表示
unsigned short int MessageGenerationDate; //19~20 记录产生日期,以1/1/1970以来
// 的Julian日期表示
unsigned long int MessageGenerationTime; //21~24 记录产生的GMT时间(毫秒)
char Unused2[4];//25~28未用
// short int MessageSegmentsNumber; //25~26报文段数记录传输分段段数
// short int MessageSegmentNumber; //27~28报文段号记录传输段号
//第三段定义第一类报文
unsigned long int RadialCollectionTime; //29~32 Zulu参考时间，
// 径向资料采集的GMT时间(毫秒)
unsigned short int RadialCollectionDate; //33~34 修正的Julian日期，
// 径向资料采集日期,以1/1/1970以来的Julian日期表示
unsigned short int UnambiguousRange; //35~36 无模糊区，不模糊距离,单位：0.1Km
unsigned short int AzimuthAngle; //37~38 方位角,实际值=(value>>3)*0.043945(度)
unsigned short int AzimuthNumber; //39~40 方位值，当前仰角内径向数据序号
unsigned short int RadialStatus; //41~42 径向数据状态，径向位置参数(=0,1,2,3,4)
short int ElevationAngle;//43~44 仰角, 实际值
//第四段定义
short int ElevationNumber; //45~46 仰角编号，体积扫描仰角数目
short int FirstGateRangeOfRef;//47~48 反射率数据第一个距离库的实际距离，
// 第一个强度库的距离(米)
short int FirstGateRangeOfDoppler;//49~50 多普勒数据第一个距离库的实际距离，
// 第一个速度\谱宽库的距离(米)
short int ReflectivityGateSize; //51~52 警戒距离，反射率数据距离库长\强度库长(米)
short int DopplerGateSize; //53~54 多普勒数据距离库长，速度\谱宽库长(米)
short int ReflectivityGates;//55~56 反射率距离\强度库数
short int DopplerGates;//57~58 多普勒距离库数，速度\谱宽库数
short int SectorNumber;//59~60 cut内的扇区号
float CalibrationConstant;//61~64 系统增益标校常数（dB）
unsigned short int RefPointer;//65~66 反射率数据指针（偏离雷达数据头的字节数），
// 表示第一个反射率数据的位置，从雷达数据头到强度数据开始的字节数
unsigned short int VelPointer;//67~68 速度数据指针（偏离雷达数据头的字节数），
// 表示第一个速度数据的位置，从雷达数据头到速度数据开始的字节数
unsigned short int SWPointer;//69~70 谱宽数据指针（偏离雷达数据头的字节数），
// 表示第一个谱宽数据的位置，从雷达数据头到谱宽数据开始的字节数
short int VelResolution;//71~72 多普勒速度分辨率:2=0.5m/s,4=1.0m/s
short int VCP;//73~74体扫VCP方式（11、21、31、32）
char Unused3[8];//75~82 未用。
unsigned short int RefPlaybackPointer;//83~84 用于回放的反射率数据指针，同33，
// 强度数据指针（用于RDA回放）
unsigned short int VelPlaybackPointer;//85~86 用于回放的速度数据指针，同34，
// 速度数据指针（用于RDA回放）
unsigned short int SWPlaybackPointer;//87~88 用于回放的谱宽数据指针，同35，
// 谱宽数据指针（用于RDA回放）
short int NyquistVelocity;//89~90 Nyquist速率(0.01m/s)不模糊速度
char Unused4[38];//91~128 未用。
// short int AAF;//91~92 Atmospheric Attenuation Factor(0.001 dB/Km)大气衰减因子
// short int TOVER;//93~94没有距离折叠两个回波功率之间最小差值(0.1瓦）
unsigned char BaseData[2300]; //129~2428基数据部分2300字节
// char *pData;
char Unused5[4]; //2429~2432 说明记录中三个基本数据是否存在，或者某个数据的库数是否有损失的情况。
// 出现这些情况时，仍然按2432个字节传输，有最后四个字节说明。
double Ldata[2300]; //******yuanlai char******
float UnRange;
float AzAngle;
float ElAngle;
float NyVelocity;
int i,k;
/* FILE *fnamelist
char filename1[80],filename2[80],filename3[80];
int fi,f2,f3,f4;
for(f1=0,f1<=9,f1++)
{
for(f2=0,f2<=9,f2++)
{
for(f3=0,f3<=9,f3++)
{
for(f4=0,f4<=9,f4++)
{
filename1=
// int a=0,b=0,c=0,d=0,e=0,f=0,g=0,h=0,sum;
*/
FILE *fp,*fpout,*fpnumber;
if((fp=fopen("2009060312.32A","rb"))==NULL)
{
printf("%s\n","2009060312.32A不存在，请确认其位置。");
exit(1);
}
if((fpout=fopen("2009060312.32A.out","w"))==NULL)
{
printf("%s\n","2009060312.32A.out不存在，请确认其位置。");
exit(2);
}
/* if((fpgates=fopen("E:\\data\\2009060315.10ADopplerGates.txt","w"))==NULL)
{
printf("%s\n","E:\\readdata\\DopplerGates.txt不存在，请确认其位置。");
exit(3);
}
*/
if((fpnumber=fopen("2009060312.32ADopplerGatesNumber.txt","w"))==NULL)
{
printf("%s\n","E:\\readdata\\DopplerGates.txt不存在，请确认其位置。");
exit(4);
}
//程序段的开始部分
k=1;
// while(cha=getc(fp)!=EOF)
// for(k=1;k<5885&&RadialStatus!=4;) //k++) 5885
do
{
// fread("%s%d%s",ch,&MessageSize,ChannelID);
// if(cha=getc(fp)==EOF&&fp>12l){printf("k=%d\n",k);break;}
// printf("This is the %d number\n",k);
fprintf(fpout,"\nThis is the %d number\n",k);
for(i=1;i<=14;i++)
fread( &Unused1[i-1],sizeof(char),1,fp);
// printf("%u ",ch[i-1]);
// printf("\n");
//原来有的 fread(&MessageSize,sizeof(unsigned short int),1,fp);//12~13 表示记录长度，双字节单位表示
// printf("This is massageSize:%u\n",MessageSize);
// fprintf(fpout,"This is massageSize:\n%u\n",MessageSize);
// pause();
//原来有的 fread(&ChannelID,sizeof(char),1,fp);//14 信号传输通道号(=0,1,2)
// printf("This is ChannelID:%d\n",ChannelID);
fread( &MessageType,sizeof(unsigned short int),1,fp);//15~16 记录数据类型,具体说明见6.1
// printf("MessageType:%u\n",MessageType);
// for(i=1;i<=12;i++)
// fread( &Unused2[i-1],sizeof(char),1,fp);//17~28 未用。
fread(&IDSequence,sizeof(short int),1,fp);//17~18 记录顺序号，0~0x7FFF之间循环表示
// printf("IDSequence:%X\n",IDSequence);
fread(&MessageGenerationDate,sizeof(unsigned short int),1,fp);
//19~20 记录产生日期,以1/1/1970以来的 Julian日期表示
// fprintf(fpout,"This is MessageGenerationDate:\n%u\n",MessageGenerationDate);
// printf("MessageGenerationDate:%u\n",MessageGenerationDate);
fread(&MessageGenerationTime,sizeof(unsigned long int),1,fp);//21~24 记录产生的GMT时间(毫秒)
// fprintf(fpout,"This is MessageGenerationTime:\n%u\n",MessageGenerationTime);
// printf("MessageGenerationTime:%u\n",MessageGenerationTime);
for(i=1;i<=4;i++)
fread( &Unused2[i-1],sizeof(char),1,fp);//25~28 未用。
// fread(&MessageSegmentsNumber, sizeof(short int),1,fp);//25~26 记录传输分段段数
// printf("MessageSegmentsNumber:%u\n",MessageSegmentsNumber);
// fread(&MessageSegmentNumber, sizeof(short int),1,fp); //27~28 记录传输段号
// printf("MessageGenerationTime:%u\n",MessageSegmentNumber);
// pause();
fread( &RadialCollectionTime, sizeof(unsigned long int),1,fp); //29~32 径向资料采集的GMT时间(毫秒)
fprintf(fpout,"This is RadialCollectionTime:\n%u\n",RadialCollectionTime);
// printf("RadialCollectionTime:%u\n",RadialCollectionTime);
fread( &RadialCollectionDate, sizeof(unsigned short int),1,fp);
//33~34 径向资料采集日期,以1/1/1970以来的Julian日期表示
fprintf(fpout,"This is RadialCollectionDate:\n%u\n",RadialCollectionDate);
// printf("RadialCollectionDate:%u\n",RadialCollectionDate);
fread( &UnambiguousRange, sizeof(unsigned short int),1,fp);//35~36 不模糊距离,单位：0.1Km
UnRange=(float)(UnambiguousRange*0.1);
fprintf(fpout,"This is UnambiguousRange:\n%f\n",UnRange);
// printf("UnambiguousRange:%u\n",UnambiguousRange);
fread( &AzimuthAngle, sizeof(unsigned short int),1,fp);//37~38 方位角,实际值=(value>>3)*0.043945(度)
// AAngle=(short int)((AzimuthAngle>>3)*0.043945/360*366+0.5);//**********重大不同*************
// AzAngle=(float)((AzimuthAngle/8.)*(180./4096.));
AzAngle=(float)((AzimuthAngle>>3)*(180./4096.));
if(AzAngle<0) AzAngle=0.;
if(AzAngle>=360) AzAngle=360.;
// printf("%d",AzimuthAngle);
// pause();
fprintf(fpout,"This is AzimuthAngle:\n%f\n",AzAngle);
// printf("AzimuthAngle:%u\n",AzimuthAngle);
fread( &AzimuthNumber, sizeof(unsigned short int),1,fp);//39~40 本层仰角扫描的径向数
fprintf(fpout,"This is AzimuthNumber:\n%u\n",AzimuthNumber);
// printf("AzimuthNumber:%u\n",AzimuthNumber);
fread( &RadialStatus, sizeof(unsigned short int),1,fp);//41~42 径向位置参数(=0,1,2,3,4)
fprintf(fpout,"This is RadialStatus:\n%u\n",RadialStatus);
// printf("RadialStatus:%d\n",RadialStatus);
fread( &ElevationAngle, sizeof(short int),1,fp);
ElAngle=(float)((ElevationAngle>>3)*(180./4096.));
// pause();
fprintf(fpout,"This is ElevationAngle:\n%f\n",ElAngle); //43~44 仰角, 实际值
// printf("ElevationAngle:%d\n",ElevationAngle);
fread( &ElevationNumber,sizeof(short int),1,fp); //45~46 体积扫描仰角数编号
fprintf(fpout,"This is ElevationNumber:\n%d\n",ElevationNumber);
// printf("ElevationNumber:%d\n",ElevationNumber);
fread( &FirstGateRangeOfRef,sizeof(short int),1,fp);//47~48第一个强度库的距离(米)
fprintf(fpout,"This is FirstGateRangeOfRef:\n%d\n",FirstGateRangeOfRef);
// printf("FirstGateRangeOfRef:%d\n",FirstGateRangeOfRef);
fread( &FirstGateRangeOfDoppler,sizeof(short int),1,fp);//49~50 第一个速度\谱宽库的距离(米)
fprintf(fpout,"This is FirstGateRangeOfDoppler:\n%d\n",FirstGateRangeOfDoppler);
// printf("FirstGateRangeOfDoppler:%d\n",FirstGateRangeOfDoppler);
fread( &ReflectivityGateSize,sizeof(short int),1,fp); //51~52 强度库长(米)
fprintf(fpout,"This is ReflectivityGateSize:\n%d\n",ReflectivityGateSize);
// printf("ReflectivityGateSize:%d\n",ReflectivityGateSize);
fread( &DopplerGateSize,sizeof(short int),1,fp); //53~54速度\谱宽库长(米)
fprintf(fpout,"This is DopplerGateSize:\n%d\n",DopplerGateSize);
// printf("DopplerGateSize:%d\n",DopplerGateSize);
fread( &ReflectivityGates,sizeof(short int),1,fp);//55~56 强度库数
fprintf(fpout,"This is ReflectivityGates:\n%d\n",ReflectivityGates);
// printf("ReflectivityGates:%d\n",ReflectivityGates);
fread( &DopplerGates,sizeof(short int),1,fp);//57~58 速度\谱宽库数
fprintf(fpout,"This is DopplerGates:\n%d\n",DopplerGates);
// printf("DopplerGates:%d\n",DopplerGates);
// if(0<DopplerGates<=320) a=a+1;
/* if(320<DopplerGates<=420) b=b+1;
else if(420<DopplerGates<=520) c=c+1;
else if(520<DopplerGates<=620) d=d+1;
else if(620<DopplerGates<=720) e=e+1;
else if(720<DopplerGates<=820) f=f+1;
else if(820<DopplerGates<=920) g=g+1;
else h=h+1;
sum=a+b+c+d+e+f+g;
fprintf(fpgates,"a=%d\nb=%d\nc=%d\nd=%d\ne=%d\nf=%d\ng=%d\nh=%d\nsum=%d\n",a,b,c,d,e,f,g,h,sum);
*/
fread( &SectorNumber,sizeof(short int),1,fp);//59~60 cut内的扇区数
fprintf(fpout,"This is SectorNumber:\n%d\n",SectorNumber);
fread( &CalibrationConstant,sizeof(float),1,fp);//61~64 系统增益标校常数（dB）
fprintf(fpout,"This is CalibrationConstant:\n%f\n",CalibrationConstant);
fread( &RefPointer,sizeof(unsigned short int),1,fp);//65~66从雷达数据头到强度数据开始的字节数
fprintf(fpout,"This is RefPointer:\n%u\n",RefPointer);
fread( &VelPointer,sizeof(unsigned short int),1,fp);//67~68从雷达数据头到速度数据开始的字节数
fprintf(fpout,"This is VelPointer:\n%u\n",VelPointer);
fread( &SWPointer,sizeof(unsigned short int),1,fp);//69~70从雷达数据头到谱宽数据开始的字节数
fprintf(fpout,"This is SWPointer:\n%u\n",SWPointer);
fread( &VelResolution,sizeof(short int),1,fp);//71~72 速度分辨率:2=0.5m/s,4=1.0m/s
fprintf(fpout,"This is VelResolution:\n%d\n",VelResolution);
// printf("VelResolution:%d\n",VelResolution);
fread( &VCP,sizeof(short int),1,fp);//73~74体扫VCP方式（11、21、31、32）
fprintf(fpout,"This is VCP:\n%d\n",VCP);
// printf("VCP:%d\n",VCP);
for(i=1;i<=8;i++)
fread( &Unused3[i-1],sizeof(char),1,fp);//75~82 未用。
// for(i=1;i<=4;i++)
// printf("Unused[%1d]:%d\n",i-1,Unused[i-1]);
fread( &RefPlaybackPointer,sizeof(unsigned short int),1,fp);//83~84强度数据指针（用于RDA回放）
// fprintf(fpout,"This is RefPlaybackPointer:\n%u\n",RefPlaybackPointer);
// printf("RefPlaybackPointer:%u\n",RefPlaybackPointer);
fread( &VelPlaybackPointer,sizeof(unsigned short int),1,fp);//85~86 速度数据指针（用于RDA回放）
// fprintf(fpout,"This is VelPlaybackPointer:\n%u\n",VelPlaybackPointer);
// printf("VelPlaybackPointer:%u\n",VelPlaybackPointer);
fread( &SWPlaybackPointer,sizeof(unsigned short int),1,fp);//87~88 谱宽数据指针（用于RDA回放）
// fprintf(fpout,"This is SWPlaybackPointer:\n%u\n",SWPlaybackPointer);
// printf("SWPlaybackPointer:%u\n",SWPlaybackPointer);
fread( &NyquistVelocity,sizeof(short int),1,fp);//89~90 Nyquist不模糊速率(0.01m/s)
NyVelocity=(float)(NyquistVelocity*0.01);
fprintf(fpout,"This is NyquistVelocity:\n%f\n",NyVelocity);
// printf("NyquistVelocity:%d\n",NyquistVelocity);
for(i=1;i<=38;i++)
fread( &Unused4[i-1],sizeof(char),1,fp);//91~128 未用。
// fread( &AAF,sizeof(short int),1,fp);//91~92 Atmospheric Attenuation Factor(0.001 dB/Km)
// printf("AAF:%d\n",AAF);
// fread( &TOVER,sizeof(short int),1,fp);//93~94没有距离折叠两个回波功率之间最小差值(0.1瓦）
// printf("TOVER:%d\n",TOVER);
fprintf(fpnumber,"%d %d %d\n",AzimuthNumber,ElevationNumber,DopplerGates);
for(i=1;i<=2300;i++)
fread( &BaseData[i-1],sizeof(unsigned char),1,fp);
// fprintf(fpout,"%u\t",BaseData[i-1]);
// for(i=1;i<=2300;i++)
// {
// printf("%d ",Mydata[i-1]);
// if(i%20==0)
// fprintf(fpout,"\n");
if(ReflectivityGates!=0)
{
fprintf(fpout,"This is the %d ReflectivityData:\n",k);
for(i=RefPointer-99;i<=RefPointer-100+ReflectivityGates;i++)
{
if(BaseData[i-1]==(char)0)
{
// BaseData[i-1]='0'; //无回波
Ldata[i-1]=-999.0;
}
else if(BaseData[i-1]==(char)1) //*******(char *)1 **********
{
// BaseData[i-1]='1'; //回波信号低于信噪比,数据是模糊的
Ldata[i-1]=999.0;
}
else
{
Ldata[i-1]=((BaseData[i-1]-2.0)/2.-32.0);
}
fprintf(fpout,"%10.4f\t",Ldata[i-1]);
}
}
else
fprintf(fpout,"This is the %d ReflectivityData:\n",k);
if(DopplerGates!=0)
{
fprintf(fpout,"\nThis is the %d DopplerVelocityData:\n",k);
for(i=VelPointer-99;i<=VelPointer-100+DopplerGates;i++)
{
if(BaseData[i-1]==(char)0) //*******(char *)0 **********
{
// BaseData[i-1]='0'; //无回波
Ldata[i-1]=-999.0;
}
else if(BaseData[i-1]==(char)1) //*******(char *)1 **********
{
// BaseData[i-1]='1'; //回波信号低于信噪比,数据是模糊的
Ldata[i-1]=999.0;
}
else
{
if(VelResolution==2) //0.5m/s
Ldata[i-1]=((BaseData[i-1]-2.0)/2.0-63.5);
else if(VelResolution==4) //1.0m/s
Ldata[i-1]=((BaseData[i-1]-2.0)-127.0);
// else
// Ldata[i-1]=((Mydata[i-1]-2)/2.0-63.5);
}
fprintf(fpout,"%10.4f\t",Ldata[i-1]);
}
}
else
fprintf(fpout,"\nThis is the %d DopplerVelocityData:\n",k);
if(DopplerGates!=0)
{
fprintf(fpout,"\nThis is the %d SWData:\n",k);
for(i=SWPointer-99;i<=SWPointer-100+DopplerGates;i++)
{
if(BaseData[i-1]==(char)0) //*******(char *)0 **********
{
// BaseData[i-1]='0'; //无回波
Ldata[i-1]=-999.0;
}
else if(BaseData[i-1]==(char)1) //*******(char *)1 **********
{
// BaseData[i-1]='1'; //回波信号低于信噪比,数据是模糊的
Ldata[i-1]=999.0;
}
else
{
Ldata[i-1]=((BaseData[i-1]-2.0)/2.0-63.5);
}
fprintf(fpout,"%10.4f\t",Ldata[i-1]);
}
}
else
fprintf(fpout,"\nThis is the %d SWData:\n",k);
// printf("\n输出保留位数：:");
// pause();
// }
for(i=1;i<=4;i++)
{
fread( &Unused5[i-1],sizeof(char),1,fp);
// fprintf(fpout,"%c",Unused5[i-1]);
}
k++;
}
//while (k<5885&&RadialStatus!=4);
while (RadialStatus!=4);
// while (k<2000);
}

复制代码

小莹子 · 发表于 2019-12-4 10:26:55

爱到花开发表于 2019-12-4 08:05

非常感谢，虽然不想矫情，但真的看着好难。我先摸索，再次谢谢您

		自动登录	找回密码
密码			立即注册