相控阵天气雷达站环境参数智能分析系统研发

吕雪芹，何艳丽，敖振浪

(广东省气象探测数据中心，广州 510080)

0 引言

天气雷达是探测台风、暴雨、龙卷风等强对流天气的重要技术手段之一。雷达探测效果不仅要求雷达设备本身具有精准的探测性能，更是要求雷达站周边环境条件好，在综合考虑工程建设基础条件可行的前提下，要求阻挡角尽量越小越好。一般情况下，雷达站站址的拔海高度越高，四周阻挡角越小。但是拔海高度不是越高越好，因为天气雷达一般不使用负仰角扫描，所以0°仰角以下的区域没有扫描到，就会出现雷达扫描盲区。而在低空范围里，往往又是强对流天气经常发生的关键区域。目前广东省气象部门快速推进X波段双极化相控阵天气雷达建设，正是弥补原有S波段新一代天气雷达探测网的低空探测盲区的不足。在原有雷达探测网的基础上，如何更好的优化布局X波段双极化相控阵天气雷达，站点选择显得尤为重要。按照以往的定性分析雷达站探测环境的常规做法已经不能完全满足精细化探测业务需要，因此如何更好地定量化计算、分析、判断雷达站环境条件参数是非常重要的。近些年来有不少关于雷达站选址方面的文章，但主要是分析四周阻挡角和等射束高度图的介绍，如伍洋等[1]借助GIS数据分析新一代天气雷达站选址,并直观展现站点的遮蔽角图和等射束高度图。景号然[2]等利用SRTM高程数据作为选址基础数据，计算得到天气雷达在0.5°，1.0°，2.4°仰角上地物遮挡情况；利用高程格点数据获得3个仰角的地物剖面数据，获得站点上空1 km,海拔3 km和6 km等射束高度图，分析了四川省天气雷达网探测环境，并给出了评估结果。李明凤[3]等利用数字高程模型(DEM)对广东省已建成的12部新一代天气雷达分别进行地形阻挡分析，主要是通过制作遮蔽角图和等射束高度图分析地形阻挡，根据阻挡分析探测效果优劣。叶飞[4]通过获取DEM高程数据、SRTM数据和雷达周边高大建筑物的遮挡数据，分别绘制天气雷达单站遮蔽角图和内蒙全区新一代天气雷达网等射束高度拼图，对内蒙全区新一代天气雷达探测环境进行了评估分析，通过评估分析可以对雷达覆盖率不够的区域提出补偿机制。高永芹、刘锋、吴焕萍、杨磊等[5-8]介绍了基于组件式GIS技术的业务应用系统建设的总体框架,并阐述了系统的建设内容,主要包括电子地图建设,制图产品输出,实时与历史资料检索，评估以及参数化制图等功能，对GIS应用于决策服务系统的若干关键技术问题进行了讨论与分析。这些研究和系统，基本上是基于GIS地理信息获取高程数据制作雷达站四周遮蔽角图和等射束高度图，再人工目视评估探测环境，费时费力，缺少定量化分析，评价完整性也存在局限性，不能很好地满足精密探测业务需求。本系统重点解决定量化、自动化评估探测环境内容。

1 架构及界面设计

1.1 总体架构

整个系统的架构设计分为四个部分：第一个部分是站点基本信息的录入，包括站点名称、经纬度、海拔高度、雷达天线计划离地面高度等参数；第二部分是通过GIS获取雷达站地理位置360°范围内距离站点一定半径内的高程数据[9-11]；第三部分是制作遮蔽角图和等射束高度图，进而制作多个站点的等射束高度图的拼图；第四部分是自动综合分析探测环境状况，输出综合分析报告。流程框图如图1所示。

图1 系统架构图

1.2 界面设计

采用C#编程语言构建WINDOWS交互主界面，菜单栏分为站点基本参数、获取高程数据、探测环境分析及显示等4个部分，功能设计着重考虑了各个部分在实际使用中的关联性和灵活性。使用广东省连州雷达站作为背景图，功能标签和按键采取图形化凸凹立体轮廓，立体感强，对约定可选内容采取可编辑文本框，方便输入。界面力求美观，交互友好便捷，如图2所示。

图2 主界面

在实际应用中，对于某个拟选站点，可使用仪器测量它的经纬度和海拔高度，根据雷达塔楼和雷达基座计算天线离地面的高度，输入到对应的“站点基本参数”菜单栏。

设计可编辑文本框的预设半径范围分别为50 km、75 km、100 km和150 km。对于X波段双极化相控阵雷达来说，有效探测范围一般是50 km以内，选取数字高程数据DEM的范围一般选择50 km即可，如果考虑为S波段雷达所用，则选择100 km就足够了，因为地球曲率原因，100 km以外只有超过8km以上的高山才可能产生阻挡[12-13]。读取的数字高程数据存入对应站点的TXT文件中，TXT文件还包含基本参数等信息，作为后续分析和评估所需基本数据。

考虑使用的灵活性，通过“获取高程数据”菜单栏中一键读取数据，系统自动搜索指定半径范围的高程数据，生成TXT文件保存于磁盘中用于任何时候使用。用户可以在TXT文件中根据实际状况增加或者修改高程数据，例如附近有建筑物阻挡而GIS里又没有标示的点。

设计“探测环境分析”菜单栏一次性完成阻挡角图、等射束高度图和多站拼图的制作，生成对应的JPG图片，同时，自动分析不同仰角的连排阻挡物的方位角以及全方位占比，进一步自动计算静锥盲区、低空盲区、扫描面积和扫描立体空域完备度等，再自动分析区域内的环境条件情况，输出综合分析结论。

2 高程数据读取方法

一般来说，从网络上商用的GIS信息系统在线读取高程数据是有限制的，不方便且速度慢，不太适合雷站址选择的应用场景需求[14-15]。本文重点介绍离线读取的方法，首先从地理信息系统下载相应区域的TIF格式高程DEM文件，TIF格式高程文件很大，涵盖一个省的矩形区域的文件大小至少2 GB，占用太多计算机资源，影响运行性能，所以选择区域应该适当，够用就行了。数字高程的分辨率当然越高越好，本设计下载30米距离分辨率1米高度分辨率的数据，基本满足业务需求。通过文件读取的方法离线读取数字高程数据，这样速度非常快[16]，而且使用起来方便，完全不依赖于互联网即可随时随地使用，这对于野外现场选址的时候至关重要。下面是实现的主要代码。在菜单点击“获取高程”按键，直接进入事件代码。

private void buttonScanning_Click(object sender, EventArgs e)

Gdal.AllRegister(); Gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); //支持中文

Dataset ds = Gdal.Open("D:雷达站选址智能测评程序GIS基础信息文件(不可删除)GIS_DEM.tif", Access.GA_ReadOnly);

int rasterX = ds.RasterXSize; //影像宽度

int rasterY = ds.RasterYSize; //影像高度

int bandCount = ds.RasterCount; //波段数

double[] tmpD = new double[6];

ds.GetGeoTransform(tmpD); //影像坐标变换参数

string proj = ds.GetProjection(); //影像坐标系信息(WKT格式字符串)

backgroundWorker1=new System.ComponentModel.BackgroundWorker();

backgroundWorker1.DoWork+= ackgroundWorker1_DoWork;//注册后台操作

backgroundWorker1.RunWorkerCompleted+= backgroundWorker1_RunWorkerCompleted;

Control.CheckForIllegalCrossThreadCalls= false;//.NET禁止了跨线程调用控件 backgroundWorker1.RunWorkerAsync(); //启动后台操作

buttonScanning.Enabled = false;// 在获取资料期间禁止再次按该按钮

textBox7.ForeColor = Color.Black;

textBox7.Text = "正在读取高程数据......";

button1.Enabled = true;

众所周知，从地图读取数据是需要花费一定时间的。如果前台长时间读数据的时候不能进行其他操作，则非常乏味且不可接受。采取的方法就是将前台读取数据的任务通过委托的方法放在后台运行，前台就能够响应用户操作。使用委托语句backgroundWorker1.DoWork += backgroundWorker1_DoWork； 进行注册后台操作，把读取任务放在后台，前台无感知。后台读取事件由如下主要代码实现。

private void backgroundWorker1_DoWork(object sender, DoWorkEventArgs e)

string filePathOne = "D:雷达站选址智能测评程序数据资料demo.TXT" StreamWriter Rad_fp = new StreamWriter(filePathOne);//实例化打开一个文本文件

//扫描四周的阻挡物的海拔高度值

for (loop = 0; loop < LOOP; loop++)

float AVG = 0, MaxAVG = 0;

Distant = 0;

NewDistant = 0;

MaxHight = 0;

for (i = 0; Distant < DataAskRange; i++)//寻找每条扫描线上的最高点

Distant = Distant + 50; //按照每50米距离间隔查询一个点的海拔高度

double[] LatLon = new double[2];

LatLon = ComputerNewPointLatLon.computerThatLonLat(Lon0, Lat0, Distant, Angle);//计算新的坐标的经纬

AVG = ProcessGisTiffFile.GetMultifyElevation(ds, LatLon[0], LatLon[1]);

hight = (float)(0.5 * Distant * Distant / radiusEarthMetres);

hight = AVG - hight;// 地球曲率订正

HigthTemp = (float)(MaxAVG - StationHeight - AntennaHeight);

BlockingAngle = (float)(Math.Asin(((HigthTemp-(0.5*NewDistant* NewDistant/radiusEarthMetres))/ NewDistant)) * 180 / Math.PI); //计算某方位角点的阻挡角

//保存数据到TXT文件

Rad_fp.Write((DataAskRange / 1000).ToString().PadLeft(5, ' ') + " ");//测评范围

Rad_fp.Write(Angle.ToString().PadLeft(5, ' ') + " ");//方位角

Rad_fp.Write((NewDistant / 1000.0).ToString("F1").PadLeft(6, ' ') + " ");//最高点的距离

Rad_fp.Write((MaxAVG).ToString("F1").PadLeft(6, ' ') + " ");//最高点的高度

Rad_fp.Write((MaxHight).ToString("F1").PadLeft(6, ' ') + " ");//订正后的高度

Rad_fp.Write(BlockingAngle.ToString("F2").PadLeft(6, ' ') + " ");//阻挡角

Angle = Angle + AngleInterval;//增加一个间隔的方位角

if (backgroundWorker1.CancellationPending) //检查是否有中止后台继续执行的按钮操作，有则退出

e.Cancel = true;

break;

Rad_fp.Flush();//清除内存释放资源

Rad_fp.Close();//关闭资料文件

为了便于后续处理，扫描到的高程数据存储到TXT文件的格式如表1所示。

表1 存储格式

需要注意的是，以站点为中心扫描四周的阻挡物，究竟以多少方位角间隔和距离间隔进行步进扫描，应该由实际需求决定[17]。本设计以0.25°间隔和50 m距离间隔在扫描范围内扫描地形的海拔高度值，从实际效果来看还是比较理想的。

3 遮蔽角和等射束高度图及拼图

3.1 遮蔽角和等射束高度图的制作

评价雷达探测能力的三种图[18]的制作流程如图3所示。

图3 阻挡图等制作流程图

首先打开某一个单站的txt文件，逐行读取全部数据，分别对每一个方位对应的原始高度值和距离值进行提取，作进一步的距离订正。利用三角函数计算方位角对应点的阻挡角，HigthTemp=hight-StationHeight-AntennaHeight;BlockingAngle=(float)(Math.Atan(HigthTemp/(NewDistant*1000))*180/Math.PI);其中HigthTemp是雷达天线订正后的高度。NewDistant是阻挡物与雷达站之间的距离。

在极坐标系绘制阻挡图和在直角坐标系上绘制阻挡图。从图4可以看出，制作的阻挡角图既具有同圆极坐标形式，也具有对应的直角坐标形式，比传统只有极坐标形式更加形象直观，容易理解，如图4所示。

图4 阻挡角图

天气雷达探测能力不仅受雷达性能、各种衰减、电磁波折射和降水云性质等因素的影响，也往往会受到雷达站四周的高大建筑物、山体的阻挡，有的影响相当严重[19-20]。为了具体、定量地分析掌握雷达站在各个方向的探测能力，需要制作表征探测能力的等射束高度图。

等射束高度图是指在一定折射条件下，测站四周由于地物阻挡，绘制的各个方向上、各种斜距下，雷达波束中心轴线能够到达的最低高度等值线图。等射束高度图考虑了雷达天线高度的情况，地球球面及大气折射的影响．特别是它考虑了周围地物、地形对雷达探测能力的影响[21]，这对分析雷达探测环境条件是非常有用的工具。应用测高公式：

H=h+1000Rsinδ+0.058 9R2

(1)

式(1)中,h为天线高度(m)，δ为遮蔽角(°)，H为观测高度(m) ，R是最大探测距离以km为单位。根据遮蔽角δ，业务上常常分别设置H为1 km、3 km和6 km，可以分别计算雷达最大观测距离R。

以0.25°方位间隔的各个方向上，应用测高公式可以分别求出各个方向上1 km高度的雷达可以达到的最大探测距离，采用极坐标的形式将探测距离依次标出，顺序连线形成等射束高度图，如图5 所示。

图5 等射束高度图

3.2 多站拼图的制作

多站点拼图实际上是将每个雷达站点的位置分别在所处的经纬度上绘制等射束高度图，单站使用不同颜色绘制，形成一幅多站拼图。必须非常注意：低纬度地区应该使用麦卡托投影坐标转换变成真实的地理坐标[22]！以减少位置误差。等高度通常取1 km、3 km和6 km。由于参与拼图的站点相距的距离可能不一样，有些相距很近，有些相距很远，相距太近绘出来拼图容易密集重叠不好看，如果相距太远又会超出屏幕范围。因此，在屏幕固定区域范围里绘制多个站点的拼图，为了尽可能的利用屏幕可显示范围，增强拼图的清晰度，采用了自动适应屏幕大小的技术。这里通过计算最左、最右、最上、最下的站点经纬度，计算出拼图的中心位置，确定向四边延伸或者缩小一定比例，进行归一化处理自适应屏幕大小。海拔高度1 km、单站半径45 km的拼图如图6所示。

图6 拼图

4 自动分析功能实现

打开单站txt文件，逐行读取各个参数。统计和分析三个设定门限的阻挡面，即阻挡角0.5°以上、1.0°以上和1.5°以上。假定如果阻挡角0.5°以上连续有10°方位阻挡，那么就视作连排阻挡。同理，阻挡角1.0°以上连续7°以上方位阻挡，就视作连排阻挡，阻挡角1.5°以上连续5°以上方位阻挡，就视作连排阻挡。根据阻挡严重程度分为优、良、一般或者差等4个等次，作为阻挡优劣的判据在输出报告中给出建议。下面是主要代码实现。

TempSaveFileName ="客观分析报告.txt";

StreamWriter fp1 = new StreamWriter(TempSaveFileName);//实例化打开一个文本文件

fp1.Write("{0}雷达站探测环境净空条件客观分析报告： ", StationName);

float sum = 0, percent;

int flagCnt = 0;

bool Startflag = false;

int More10BlockingAngle = 0;

for (loop = 0; loop < LOOP; loop++)

if (Angle_BlockingAngle[loop, 1] >= 0.5)//阻挡角>=0.5的占比统计

{ ++sum;

if (Startflag == true)

{ flagCnt++;

} else

{ Startflag = true; flagCnt = 1;

else

{ Startflag = false;

if (flagCnt >= 20)//步进0.5°，方位连排0.5X20=10°以上

{ fp1.Write("方位从{0}度到{1}度的区段，共有{2}度范围，并且挡角高于0.5°的连排(连续多于10°方位角)山体阻挡。 ", Angle_BlockingAngle[loop - flagCnt, 0], Angle_BlockingAngle[loop - 1, 0], flagCnt/2);

flagCnt = 0;

More10BlockingAngle++;

percent = sum * 100 / LOOP;

fp1.Write("本站四周高于0.5°阻挡角一共有{0}个点(方位0.5°间隔一个点),", sum);

fp1.Write("即是高于0.5°的阻挡面占比{0}%。 ", percent.ToString("#0.0"));

探测净空面积的计算方法是采用作者提出的细分逼近法，限于篇幅，细分逼近算法不在这里介绍。使用这种算法计算出X波动双极化相控阵雷达45 km扫描半径范围内，1 km高度探测净空投影面积、净空占总投影面积比例、投影阻挡面积比例、有效探测空域及其占比等。对于其他类型天气雷达，不同扫描范围和3 km、6 km高度探测净空面积的计算方法类同。通过自动分析功能，给出对应的综合客观评价。下面是独树岗雷达站探测环境净空条件客观分析报告实例。

方位从82.25°～98°的区段，共有32°范围，并且挡角高于0.5°的连排(连续多于10°方位角)山体阻挡。

本站四周高于0.5°阻挡角一共有515个点(方位0.25°间隔一个点),即是高于0.5°的阻挡面占比35.8%。

方位从86.75～98°的区段，共有23°范围，并且挡角高于1.0°的连排(连续7度方位角以上)山体阻挡。

本站四周高于1.0°阻挡角一共有118个点(方位0.25°间隔一个点),即是高于1.0°的阻挡面占比8.2%。

本站四周有1个方向存在高于1.0°阻挡的连排(连续多于7°方位角)区域，反映探测环境条件不太理想。

方位从89.5°～97°的区段，共有15°范围，并且挡角高于1.5°的连排(连续5°方位角以上)山体阻挡。

本站四周高于1.5°阻挡角一共有42个点(方位0.25°间隔一个点),即是高于1.5°的阻挡面占比2.9%。

本站四周有1个方向存在高于1.5°阻挡的连排(连续多于5°方位角)区域，说明探测环境条件不理想。

本站45 km扫描半径、馈源上空25 km以下范围内，有效探测立体空域占总空域体积的百分比(雷达站45 km探测半径、馈源上空25 km以下范围内):80.64%。

本站45 km扫描半径、馈源上空25 km以下范围内，有效探测立体空域占默认扫描空域的百分比(雷达站0～37°仰角、雷达馈源上空0～25 km高度范围内):98.85%。

本站45 km扫描半径、馈源上空25 km以下范围内，静椎区空域体积占总空域体积的百分比:18.05%。

本站45 km扫描半径、馈源上空25 km以下范围内，低空扫描盲区占总空域体积的百分比(天线馈源以下、仰角0°以下):0.38%。

5 结束语

以往雷达选址的通常做法，是使用经纬仪在拟选点实地测量四周的遮蔽角，这种人工方法相当费时费力和粗糙，非常容易受地理位置和天气状况影响，较远地方阻挡物看不见而忽视了其阻挡影响。也有采用人工在地图上作业，寻找各个方位角上最大海拔高度的点，再使用尺子测量与雷达站之间的距离，然后计算阻挡角。所有这些方法很难保证客观地反映实际探测环境状况，而且效率非常缓慢，往往需要几天时间才能完成，难以满足现代业务需要。

本文基于离线的GIS地理信息系统数字高程模型DEM获取雷达站探测范围内的山脉、森林、高楼等阻挡情况，非常适合现场初步寻找天气雷达最佳候选站址，灵活性便利性不言而喻。自动画出雷达探测范围内的阻挡图、测站上空1 km和海拔3 km、6 km高度的等高度射束图、多站拼图，自动计算探测净空投影面积、有效探测空域等覆盖率指标，并进行智能分析，给出最优合理化评估报告。系统能够在数十秒内完成整个计算、画图及分析过程，其客观性和效率凸显。引入有效探测空域的分析和计算功能是传统雷达选址方法不具有的，定量化的分析极大地提高了天气雷达站选址决策的科学性，提高全网探测能力，为雷达监测精密提供良好基础，该系统在广东省26部X波段双极化相控阵雷达选址中得到很好的应用。

为了便于读者参考借鉴，给出有关的实现代码以期对读者有用。