气象雷达探测的有效立体空域细分逼近计算方法与流程

本发明涉及气象雷达应用技术领域，具体为气象雷达探测的有效立体空域细分逼近计算方法。

背景技术：

以往，气象雷达建设前期选址工作中基本上是根据阻挡角图和等射束高度图粗略估测气象雷达周边的探测环境状况，很难客观准确掌握立体空域的阻挡情况，因为没有人研究如何计算阻挡空域问题。然而气象工作人员更关注气象雷达探测的有效空域和扫描盲区，关于雷达盲区的精确计算方法,教科书上介绍的内容极少,可找到的相关资料也非常有限，雷达探测的有效立体空域分析和计算方法更是未见介绍，所以本发明重点解决这个问题。

通常气象雷达安装在四周开阔、阻挡较小的地方。天气过程绝大多数发生在地面到40km高空范围内，选择站址的时候，通常会选在四周阻挡相对小的高山顶上，但是总会存在不同程度的阻挡。考虑到低空也有许多重要天气过程，所以在阻挡容许范围内，站点海拔高度也不是越高越好。气象雷达很少使用负仰角扫描。因此我们重点关注相对于雷达0°仰角以上的空域就可以了。

气象雷达探测的有效空域受到各个方位角方向的高山高楼等阻挡物阻挡影响有关，也与雷达能够扫描的最大仰角有关，还与地球是圆形有关。雷达图像在实际业务使用中，绝大多数应用人员根本不知道也不清楚探测范围内哪些位置被阻挡掉了，扫描垂直高度上有没有阻挡，诸如此类影响回波数据质量的问题，直接影响气象预报人员决策的科学性，因此必须搞清楚探测空域的有效性。目前我国业务使用的新一代天气雷达功能规格书给出的定量探测范围是半径230km，定性探测范围是半径460km，扫描仰角通常设定在0°到20°之间，扫描层数一般是9至17层；最新发展起来的x波段双极化相控阵天气雷达，定量探测范围是半径45km左右，它的扫描仰角范围通常设定在0°到37°之间。由于地物的阻挡也是非常复杂的，很难简单地应用几何计算有效空域，所以设计气象雷达探测的有效立体空域细分逼近计算方法是很有必要的。

技术实现要素：

本发明的目的在于提供气象雷达探测的有效立体空域细分逼近计算方法，以解决上述背景技术中提出的问题。

为了解决上述技术问题，本发明提供如下技术方案：气象雷达探测的有效立体空域细分逼近计算方法，包括以下步骤：步骤一，高程数据的获取；步骤二，设定高度范围和层次；步骤三，打开文件；步骤四，一次计算；步骤五，绘画；步骤六，读出数据；步骤七，二次计算；步骤八，三次计算；步骤九，四次计算；步骤十，五次计算；步骤十一，计算角度；步骤十二，计算体积；步骤十三，计算半径；步骤十四，计算占比；

其中在上述步骤一中，通过地理信息系统获取dem数据，以雷达站为中心，按照一周360°每0.25°方位角从dem中读取l1范围内的阻挡物高程值，且阻挡物高程值单位为千米，经过地球曲率修正，筛选出对雷达发射波束阻挡最大高度值和对应的距离数据，生成阻挡角数据文件“data.txt”；

其中在上述步骤二中，根据雷达业务实际运行情况，假定最小扫描仰角为0°，测站馈源上空0-l2高度范围每隔100m一层，且l2的取值单位为千米；考虑100m高度变化对等射束高度图的面积变化较小，如果分为50m一层则等射束高度图面积变化更小，实际上细分高度越小对有效探测空域计算结果越精准，也就是立体空域体积细分逼近计算方法；

其中在上述步骤三中，打开data.txt文件，读出第一组高程数据和距离数据，计算相对应方位的阻挡角；

其中在上述步骤四中，计算第n高度层雷达可以探测到的最远距离distant，在第n高度层上大于distant的该方位后面的空域是扫描不到的；

其中在上述步骤五中，画出该测站的等射束高度图，以雷达站为原点，方位角0.25°夹角，以这个最大探测距离为等腰三角形的腰长distant，构成的等腰三角形，计算出这个等腰三角形的面积，面积计算公式为s＝0.5*distant*distant*sin0.25°；且等腰三角形的面积单位为平方千米；一周总共1440个不同的等腰三角形，求面积之和，就是逼近该高度层的有效探测面积；循环往复，总共计算0-l2之间的各层有效面积，保存到各层面积数据文件s.txt；

其中在上述步骤六中，打开s.txt文件，逐条读出400组面积数据，按照每层100m高度近似计算该层可探测空间体积，且可探测空间体积计算公式为v＝面积si*高度＝si*0.1，且可探测空间体积单位为立方千米，一共高度层体积累加之和得到有效空域体积v1，也就是整个扫描仰角0°到20°之间、探测距离最大l1、测站馈源上空0-l2高空范围雷达能够探测到的空域，这个体积就是雷达扫描没有受到阻挡的那部分空域；

其中在上述步骤七中，计算静椎区体积v2：由于目前我国业务使用的s波段气象雷达观测模式设定扫描仰角范围是0°到20°，所仰角20°以上、测站馈源上空0-l2高空范围是业务雷达探测不到的静锥区空域；根据圆锥体积公式是v＝1/3sh1，s是圆柱的底面积，h1是圆柱的高，那么静锥区体积v2＝1/3*π*(l2/tan20°)*(l2/tan(20°)*l2，且静锥区体积v2的单位为立方千米；

其中在上述步骤八中，假如没有任何阻挡，那么整个扫描仰角0°到20°之间、探测距离最大l1、测站馈源上空0-l2高空范围雷达应探测空域体积，可以根据圆柱体公式计算，应探测空域体积＝π*l1*l1*l2，且应探测空域体积的单位为立方千米；默认扫描空域体积v10＝应探测空域体积-静椎区体积v2，且默认扫描空域体积v10的单位为立方千米；

其中在上述步骤九中，计算扫描盲区体积v3：v3＝天线海拔高度*雷达扫描覆盖区面积＝h2*π*l1*l1；且h2为雷达天线馈源安装在海拔的高度，且单位为千米；

其中在上述步骤十中，计算雷达站下面海拔高度0m到地球中心与l1探测半径范围形成的圆锥柱体体积v5：根据圆锥体积公式是v＝1/3sh，s是圆柱的底面积，h1是圆柱的高，那么v5＝1/3*π*探测半径*探测半径*地球半径＝π*l1*l1*6371.393/3，且v5的单位为立方千米；

其中在上述步骤十一中，首先计算θ＝atan(l1/6371.393)；球缺的高h＝地球半径r(1-cosθ)；根据球缺体积公式v＝πh²(r-h/3)，从而计算得到球缺体积v6，其中r是地球的半径,h是球缺的高；球缺下面圆锥体体积v7＝1/3*(6371.393*cosθ)*π*(6371.393*sinθ)*(6371.393*sinθ)；

其中在上述步骤十二中，计算地曲盲区体积v4＝圆锥柱体体积v5-球缺体积v6-球缺下面圆锥体体积v7；

其中在上述步骤十三中，计算雷达探测半径l1、测站馈源上空l2以下的总空域体积v0＝默认扫描空域体积v10+静椎区体积v2+扫描盲区体积v3+地曲盲区体积v4；

其中在上述步骤十四中，通过以上各体积的计算结果，可以计算出各部分体积的占比，其中有效探测立体空域占总空域体积的百分比：percentage1＝有效空域体积v1/总空域体积v0；且有效探测立体空域占总空域体积的百分比雷达站l1探测半径、馈源上空l2以下范围内计算；有效探测立体空域占默认扫描空域的百分比：percentage2＝有效空域体积v1/默认扫描空域体积v10，且有效探测立体空域占默认扫描空域的百分比在雷达站0°到20°仰角、雷达馈源上空0-l2高度范围内计算；静椎区空域体积占总空域体积的百分比：percentage3＝静椎区体积v2/总空域体积v0；低空扫描盲区占总空域体积的百分比：percentage4＝(扫描盲区体积v3+地曲盲区体积v4)/总空域体积v0，且低空扫描盲区占总空域体积的百分比在天线馈源以下、仰角0°以下范围内计算。

根据上述技术方案，所述步骤一中，l1的取值范围为0-230。

根据上述技术方案，所述步骤一中，文件“data.txt”中包含1440组挡角数据。

根据上述技术方案，所述步骤五中，细分面积运用了等腰三角形，且腰长distant的单位为千米。

根据上述技术方案，所述步骤十中，π的取值为3.14159。

与现有技术相比，本发明所达到的有益效果是：该气象雷达探测的有效立体空域细分逼近计算方法，通过该计算方法能够准确计算出雷达扫描空域的实际情况，通过各个占比情况，结合雷达阻挡图和等射束高度图，能够较好地评价拟选雷达站条件的好坏，具有非常重要的实际意义。

附图说明

附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。在附图中：

图1是本发明的方法流程图；

图2是本发明中高程数据获取流程图；

图3是本发明中等射束高度图面积计算流程图；

图4是本发明中有效探测空域计算流程图；

图5是本发明中不同高度层的最大探测距离图；

图6是本发明中等射束高度图垂直投影示意图；

图7是本发明中雷达扫描区域示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

请参阅图1-7，本发明提供一种技术方案：气象雷达探测的有效立体空域细分逼近计算方法，包括以下步骤：步骤一，高程数据的获取；步骤二，设定高度范围和层次；步骤三，打开文件；步骤四，一次计算；步骤五，绘画；步骤六，读出数据；步骤七，二次计算；步骤八，三次计算；步骤九，四次计算；步骤十，五次计算；步骤十一，计算角度；步骤十二，计算体积；步骤十三，计算半径；步骤十四，计算占比；

其中在上述步骤一中，通过地理信息系统获取dem数据，以雷达站为中心，按照一周360°每0.25°方位角从dem中读取230km范围内的阻挡物高程值，经过地球曲率修正，筛选出对雷达发射波束阻挡最大高度值和对应的距离数据，生成阻挡角数据文件“data.txt”，文件中包含1440组挡角数据；

其中在上述步骤二中，根据雷达业务实际运行情况，假定最小扫描仰角为0°，测站馈源上空0-40km高度范围每隔100m一层，分为共400层；考虑100m高度变化对等射束高度图的面积变化较小，如果分为50m一层则等射束高度图面积变化更小，实际上细分高度越小对有效探测空域计算结果越精准，也就是立体空域体积细分逼近计算方法；

其中在上述步骤三中，打开data.txt文件，读出第一组高程数据和距离数据，计算相对应方位的阻挡角；

其中在上述步骤四中，计算第n高度层雷达可以探测到的最远距离distant，在第n高度层上大于distant的该方位后面的空域是扫描不到的；

其中在上述步骤五中，画出该测站的等射束高度图，以雷达站为原点，方位角0.25°夹角，以这个最大探测距离为等腰三角形的腰长distant，构成的等腰三角形，计算出这个等腰三角形的面积，s＝0.5*distant*distant*sin(0.25°)；且细分面积运用了等腰三角形，且腰长distant的单位为千米，腰三角形的面积单位为平方千米，一周总共1440个不同的等腰三角形，求面积之和，就是逼近该高度层的有效探测面积；循环往复，总共计算0km至40km之间的400层有效面积，保存到各层面积数据文件s.txt；

其中在上述步骤六中，打开s.txt文件，逐条读出400组面积数据，按照每层100m高度近似计算该层可探测空间体积v＝面积si*高度＝si*0.1，且可探测空间体积v的单位为立方千米，400个高度层体积累加之和得到有效空域体积v1，也就是整个扫描仰角0°到20°之间、探测距离最大230km、测站馈源上空0m至40km高空范围雷达能够探测到的空域，这个体积就是雷达扫描没有受到阻挡的那部分空域；

其中在上述步骤七中，计算静椎区体积v2：由于目前我国业务使用的s波段气象雷达观测模式设定扫描仰角范围是0°到20°，所仰角20°以上、测站馈源上空0m至40km高空范围是业务雷达探测不到的静锥区空域；根据圆锥体积公式是v＝1/3sh，s是圆柱的底面积，h是圆柱的高，那么静锥区体积v2＝1/3*π*(40/tan20°)*(40/tan(20°)*40＝505899.3(km³)；

其中在上述步骤八中，假如没有任何阻挡，那么整个扫描仰角0°到20°之间、探测距离最大230km、测站馈源上空0m至40km高空范围雷达应探测空域体积，可以根据圆柱体公式计算，应探测空域体积＝π*230*230*40＝6647414(km³)；默认扫描空域体积v10＝应探测空域体积-静椎区体积v2＝6647414-505899.3＝6141514.7(km³)；

其中在上述步骤九中，计算扫描盲区体积v3：v3＝天线海拔高度x雷达扫描覆盖区面积＝h*π*230*230；举例，假设雷达天线馈源安装在海拔高度h＝200m，则扫描盲区体积v3＝3.14159*230*230*0.2＝33237.1(km³)；

其中在上述步骤十中，计算雷达站下面海拔高度0m到地球中心与230km探测半径范围形成的圆锥柱体体积v5：根据圆锥体积公式是v＝1/3sh，s是圆柱的底面积，h是圆柱的高，那么v5＝1/3*π*探测半径*探测半径*地球半径＝π*230*230*6371.393/3＝352944058.6(km³)；

其中在上述步骤十一中，首先计算θ＝atan(230/6371.393)＝2.0674°；球缺的高h＝地球半径r(1-cosθ)＝6371.393*(1-cos2.0674°)＝4.1473km；根据球缺体积公式v＝πh²(r-h/3)，其中r是地球的半径,h是球缺的高，那么球缺体积v6：v6＝π*4.1473*4.1473*(6371.393-4.1473/3)＝344197.8(km3)；球缺下面圆锥体体积v7＝1/3*(6371.393*cosθ)*π*(6371.393*sinθ)*(6371.393*sinθ)＝1/3*(6371.393*cos2.0674°)*π*(6371.393*sin2.0674°)*(6371.393*sin2.0674°)＝6364.5*3.14159*229.8*229.8/3＝351949334.9(km³)；

其中在上述步骤十二中，计算地曲盲区体积v4＝圆锥柱体体积v5-球缺体积v6-球缺下面圆锥体体积v7＝352944058.6-344197.8-352250387.3＝349473.5(km³)；

其中在上述步骤十三中，雷达探测半径230km、测站馈源上空40km以下的总空域体积v0＝默认扫描空域体积v10+静椎区体积v2+扫描盲区体积v3+地曲盲区体积v4；举例，假设雷达天线馈源安装在海拔高度200m，那么总空域体积v0＝6141514.7+505899.3+33237.1+349473.5＝7030124.6(km³)；

其中在上述步骤十四中，通过以上各体积的计算结果，可以计算出各部分体积的占比，其中有效探测立体空域占总空域体积的百分比：percentage1＝有效空域体积v1/总空域体积v0；且有效探测立体空域占总空域体积的百分比雷达站l1探测半径、馈源上空l2以下范围内计算；有效探测立体空域占默认扫描空域的百分比：percentage2＝有效空域体积v1/默认扫描空域体积v10，且有效探测立体空域占默认扫描空域的百分比在雷达站0°到20°仰角、雷达馈源上空0-l2高度范围内计算；静椎区空域体积占总空域体积的百分比：percentage3＝静椎区体积v2/总空域体积v0；低空扫描盲区占总空域体积的百分比：percentage4＝(扫描盲区体积v3+地曲盲区体积v4)/总空域体积v0，且低空扫描盲区占总空域体积的百分比在天线馈源以下、仰角0°以下范围内计算。

基于上述，本发明的优点在于，本发明，通过计算得到各个占比情况，从而根据得到的雷达扫描空域的实际情况，结合雷达阻挡图和等射束高度图，能够较好地评价拟选雷达站条件的好坏，具有非常重要的实际意义。

需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。

最后应说明的是：以上所述仅为本发明的优选实施例而已，并不用于限制本发明，尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。