CN115690339A - 一种柱坐标空间下的高精度空域通视分析方法 - Google Patents

Abstract

本发明公开了一种柱坐标空间下的高精度空域通视分析方法，包括以下步骤：S1，选择需要分析的地域，采集待分析的DEM高程数据；S2，将所述DEM高程数据转换到柱坐标空间；S3，计算所述DEM高程数据中的网格方位信息；S4，利用所述网格方位信息，将所述DEM高程数据按照方位划分为M个子区域；S5，计算某一网格上空高度Δh处的通视情况；S6，利用步骤5计算的所有网格的通视情况，获取该地域上高度Δh的空域通视情况。本方法提出了用最大方位角、最小方位角来描述一个DEM网格的方位信息，更加准确地描述了网格信息，从而可以获得更高的分析精度，提高了分析处理效率和精度，避免了大量数值计算，易于实现，适合大规模的数据分析。

Description

一种柱坐标空间下的高精度空域通视分析方法

技术领域

本发明涉及通视分析领域，特别涉及一种柱坐标空间下的高精度空域通视分析方法。

背景技术

利用通视分析方法分析某一区域的空域通视情况，在雷达站选址、通信基站选址和军事作战仿真等领域有着重要应用。

通视分析研究的一个重要方向就是利用DEM(Digital Elevation Model)地形高程数据(以下简称DEM数据)，分析观测点周围的空域通视情况。DEM数据实际上是在地里坐标系下描述地形的高程数据，DEM数据把地形按照经纬度剖分成网格，某一DEM数据值就标识了某个经纬度网络范围内的地形高程值，网格大小反应了DEM数据的精细程度。

现有的通视分析方法中，大部分算法假设DEM划分的网格是规则的，然后利用这些网格的高程数据判断目标是否通视。(陈龙，王永军，李大军.R3视域算法与参考面算法的对比研究[J].测绘与空间地理信息，2014,37(5)：34-38；Chuanjun Wu.et al.PDERL:anaccurate and fast algorithm with a novel perspective on solving the oldviewshed analysis problem[J].Earth Science Information,2021,14:619-623)。事实上由于地球纬度线长度不同，不同地点的经纬度网格大小并不相同，地域面积越大，网格大小差别越明显。特别是把经纬度网格投影到观测点所在水平面时，会进一步发生形变。因此，如果假设地形网格是规则的，必然存在分析误差。另外，有的通视算法利用地理坐标系和地心坐标系之间的转换关系寻找目标视线经过DEM数据网格，这种分析算法虽然不存在地形误差，但是需要大量的数值计算，因此不适合大规模计算(李易依，贾维敏，姚敏立.一种改进的DYNTACS通视性算法[J].计算机仿真,2020,37(5):218-221)。

通过上述分析可见，寻找一种计算简单、且高精度的通视分析方法是通视分析领域研究和发展的一个重要方向。

发明内容

为了解决现有问题，本发明提供了一种柱坐标空间下的高精度空域通视分析方法，具体方案如下：

一种柱坐标空间下的高精度空域通视分析方法，包括以下步骤：

S1，选择需要分析的地域，采集待分析的DEM高程数据；

S2，将所述DEM高程数据转换到柱坐标空间；

S3，计算所述DEM高程数据中的网格方位信息；

S4，利用所述网格方位信息，将所述DEM高程数据按照方位划分为M个子区域；

S5，计算某一网格上空高度Δh处的通视情况；

S6，利用步骤5计算的所有网格的通视情况，获取该地域上高度Δh的空域通视情况。

优选的，所述DEM高程数据存储的高度值H_ij表示网格G_ij的高程，其经纬度坐标为(B_i，L_j)，约定左下角顶点V_ij的(B_i，L_j，H_ij)坐标代表网格G_ij的坐标，所述步骤2将待分析地域的所述DEM数据经过三次坐标转换，转换到以地理坐标为(B_P，L_P，H_P)的观测点V_P为中心的柱坐标空间下，具体要经过三次坐标转换，步骤包括：

S21，将所述DEM高程数据从地理坐标系转换到空间直角坐标系，转换方法见公式(1):

其中，N＝a/W，N为地球卯酉圈曲率半径，a为地球长半轴，a＝6378.137km，

e²＝(a²-b²)/a²，e为地球的第一偏心率，b为地球的短半径，b＝6356.7523141km，

转换到空间直角坐标系后，观测点V_P的坐标则变换为(X_P，Y_P，Z_P)；

S2.2，将所述DEM高程数据从空间直角坐标系转换到以V_P为中心的站心坐标系，转换方法见公式(2)：

其中，(B_P，L_P，H_P)为观测点在地理坐标系下的坐标，(x，y，z)是所述DEM高程数据在站心坐标系下的坐标，此时观测点V_P的坐标为(0，0，0)；

S2.3，将所述DEM高程数据转换到柱坐标空间(R，θ，h)上，R表示某一网格距离观测点的水平距离，θ表示网格相对于观测点方位，以正北为方位0点，h表示网格相对于观测点水平面的高度，经过坐标转换后，h已包含地球曲率的影响，转换方法见公式(3)：

h＝z (3)。

优选的，所述步骤3中的所述网格方位信息包括所述DEM高程数据中的每一个网格相对于观测点的最大方位角θ(G)_max、最小方位角θ(G)_min和俯仰角β(G)。

优选的，第i行第j列网格中的所述网格方位信息分别为θ(G_ij)_max、θ(G_ij)_min、和β(G_ij)，其定义见公式(4)：

θ(G_ij)_max＝max([θ_i,jθ_i,j+1θ_i+1,jθ_i+1,j+1])

θ(G_ij)_min＝min([θ_i,jθ_i,j+1θ_i+1,jθ_i+1,j+1])

β(G_ij)＝arctan(h_ij/R_ij) (4)

其中，θ_i,j，θ_i,j+1，θ_i+1,j，θ_i+1,j+1表示网格G_ij的四个顶点的方位角。

优选的，步骤4中分区的方法包括以下步骤：

S4.1，计算每个子区域的起始方位角和结束方位角，见公式(5)

zone(i)_start＝(i-1)*(360/N)

zone(i)_end＝i*(360/N) (5)

其中，zone(i)代表第i个区域，zone(i)_start表示该区域的起始方位，zone(i)_end表示该区域的结束方位；

S4.2，根据区域的起始、结束方位角zone(i)_start、zone(i)_end，和网格的最大方位角θ(G)_max、最小方位角θ(G)_min，判断网格G是否属于该区域，当网格的最大、最小方位角满足公式(6)(7)(8)中的任一组条件，则判定该网格属于第i个子区域zone(i)：

优选的，所述步骤5中网格G_ij上空Δh处的通视情况的方法包括以下步骤：

S5.1，在柱坐标空间下的所述DEM高程数据中，选择一个网格G_ij，其坐标为(R_Gij，θ_Gij，h_Gij)，选择该网格上空某一点作为目标点V_T(R_Gij，θ_Gij，h_T)，其中h_T＝h_Gij+Δh，Δh表示目标点相对G_ij高度；

S5.2，根据G_ij的方位角，选择一个所述DEM高程数据子区域，利用该子区域内的所述DEM高程数据进行通视分析，如果zone(i)_start≤θ_Gij≤zone(i)_end，则选择第i个子区域；

S5.3，寻找影响目标通视性的网格集合{G}，V_T在水平面投影为V_T'，线段V_PV_T'在柱坐标下表示为满足{R≤R_T,θ＝θ_T}点的集合，因此{G}就是涵盖V_T'所在方位，且水平距离小于R_T的网格的集合，即G满足条件(9)：

公式(9)利用了网格最大方位角、最小方位角和距离信息，准确找到影响通视性的网格集合{G}，避免了其它算法大量的数值计算，也不会受到网格大小不一致的影响；

S5.4，通过步骤S5.3找到的影响通视的所有网格{G}，比较目标点俯仰角β(V_T)与{G}中所有网格的俯仰角大小，根据公式(10)判断是否通视：

see(V_T)表示目标点是否通视，通视记为1，不通视记为0，n表示集合{G}中包含了n个网格，如果分析空域高度是相对地形的高度，则β(V_T)＝arctan((h_Gij+Δh)/R_T)，Δh表示目标点相对G_ij高度；如果分析的空域高度是相对观测点高度，则β(V_T)＝arctan(h_T/R_T)，h_T表示目标点相对V_P的高度。

本发明还揭示了一种计算机可读存储介质，介质上存有计算机程序，计算机程序运行后，执行上述的柱坐标空间下的高精度空域通视分析方法。

本发明还揭示了一种计算机系统，包括处理器、存储介质，存储介质上存有计算机程序，处理器从存储介质上读取并运行计算机程序以执行上述的柱坐标空间下的高精度空域通视分析方法。

本发明的有益效果在于：

(1)该方法将DEM数据转换到柱坐标空间，首次提出了用最大方位角、最小方位角来描述一个DEM网格的方位信息，因此本方法中用水平距离、方位、高度、最大方位、最小方位五个参数，更加准确地描述了网格信息，从而可以获得更高的分析精度。

(2)本方法提出了以观测点为中心将DEM数据分成M个子区域，在子区域内进行通视分析的思路，该思路可以极大降低一次分析处理的网格单元数量，提高分析处理效率；

(3)该方法通过比较网格最大方位角、最小方位角、网格水平距离、目标点方位、目标点水平距离等信息，能准确获取观测点与目标点投影经过的所有地形网格，避免了网格大小不一致引起的误差问题，通视分析的精度较高；

(4)本方法在完成坐标转换等数据预处理后，仅通过比较运算即可实现通视分析，避免了大量数值计算，易于实现，适合大规模的数据分析。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图。

图1：一种规则网格DEM地形数据；

图2：柱坐标下DEM数据的定义(R，θ，h)；

图3：网格G_ij的俯仰角及四个顶点方位角；

图4：DEM数据子区域的划分；

图5：寻找影响目标通视性的网格集合{G}；

图6：网格G的β(V_T)的计算；

图7：经纬度坐标下的DEM地形数据；

图8：柱坐标下的DEM地形数据；

图9：影响影响V_T通视性的网格集合{G}；

图10：观测点周围相对地表200米高度空域通视情况。

具体实施方式

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

一种柱坐标空间下的高精度空域通视分析方法，包括以下步骤：

S1，选择需要分析的地域，下载采集待分析的DEM高程数据。一种常用的DEM数据格，如图1所示，阴影部分表示实际的地形。这种基于网格的DEM数据根据经纬度把地形划分成小网格。DEM高程数据存储的高度值H_ij表示网格G_ij的高程，其经纬度坐标为(B_i，L_j)，约定左下角顶点V_ij的(B_i，L_j，H_ij)坐标代表网格G_ij的坐标。

本实施例从中科院地理空间数据云下载了ASTER GDEM2格式的DEM数据，地域涵盖北纬43-45度，东经87-89度，将该数据读入计算机中，经过数据拼接和间隔抽样后，最终DEM数据中包含了size＝3600×3600个网格高程数据。经纬度坐标下的三维地形，如图7所示。每一个地形网格的坐标为(B，L，H)。每个网格大小为2″×2″，精度大约60m。网格的坐标包含三个变量，因此，可以用一个3×size矩阵存储DEM网格数据，本实施例中定义一个矩阵变量BLH(3×size)，其中BLH(1，:)表示网格的经度，BLH(2，:)表示网格的纬度，BLH(3，:)表示网格高度。

S2，将DEM高程数据转换到柱坐标空间。

将待分析地域的DEM数据经过三次坐标转换，转换到以地理坐标为(B_P，L_P，H_P)的观测点V_P为中心的柱坐标空间下，具体要经过三次坐标转换，步骤包括：

S21，将DEM高程数据从地理坐标系转换到空间直角坐标系，转换方法见公式(1)：

其中，N＝a/W，N为地球卯酉圈曲率半径，a为地球长半轴，a＝6378.137km，

e²＝(a²-b²)/a²，e为地球的第一偏心率，b为地球的短半径，b＝6356.7523141km，

转换到空间直角坐标系后，观测点V_P的坐标则变换为(X_P，Y_P，Z_P)。

S2.2，将DEM高程数据从空间直角坐标系转换到以V_P为中心的站心坐标系，转换方法见公式(2)：

其中，(B_P，L_P，H_P)为观测点在地理坐标系下的坐标，(x，y，z)是DEM高程数据在站心坐标系下的坐标，此时观测点V_P的坐标为(0，0，0)。

S2.3，将DEM高程数据转换到柱坐标空间(R，θ，h)上，R表示某一网格距离观测点的水平距离，θ表示网格相对于观测点方位，以正北为方位0点，h表示网格相对于观测点水平面的高度，经过坐标转换后，h已包含地球曲率的影响，如图2所示。转换方法见公式(3)：

h＝z (3)。

具体地，选择观测点VP(N44°26'40"，E88°，473m)，利用公式(1)、(2)和(3)，把DEM数据转换到以观测点VP为中心的柱坐标系下，转换后的DEM数据如图8所示，每一个网格数据的坐标为(R，θ，h)，其中h表示了网格相对于VP平面的高度。定义三维矩阵变量RDE(3×size)，存储柱坐标下的DEM数据，其中RDE(1，:)表示网格与VP的距离，RDE(2，:)表示网格的方位，RDE(3，:)表示网格相对VP的高度。

通过对比图7与图8右侧的高度条数据，不难发现坐标转换后的高度h已经消除了地球曲率的影响。

在柱坐标空间下进行空域通视分析，用距离、方位和高度等信息描述DEM数据，更便于实现目标点的通视分析。

S3，计算DEM高程数据中的网格方位信息。

步骤3中的网格方位信息包括DEM高程数据中的每一个网格相对于观测点的最大方位角θ(G)max、最小方位角θ(G)min和俯仰角β(G)。

DEM网格数据中，每个网格代表着一个小范围区域，仅通过网格顶点(R，θ，h)坐标无法全面反应网格的信息。因此，这里定义网格的最大方位角θ(G)max、最小方位角θ(G)min和俯仰角β(G)三个参数。其中，俯仰角β的大小反应了网格对观测点遮蔽的大小。第i行第j列网格中的网格方位信息分别为θ(G_ij)_max、θ(G_ij)_min、和β(G_ij)，其定义见公式(4)：

θ(G_ij)_max＝max([θ_i,jθ_i,j+1θ_i+1,jθ_i+1,j+1])

θ(G_ij)_min＝min([θ_i,jθ_i,j+1θ_i+1,jθ_i+1,j+1])

β(G_ij)＝arctan(h_ij/R_ij) (4)

根据图3可知，网格G_ij的四个顶点的方位(θ_i,j，θ_i,j+1，θ_i+1,j，θ_i+1,j+1)，实际上代表了四个网格的方位(G_i,j，G_i,j+1，G_i+1,j，G_i+1,j+1)。

因此，通过索引值可以在RDE(2，:)中直接获取四个顶点的方位值。

其中，G_i,j的索引值为IndexG_i,j＝j*3600+i；

G_i,j+1的索引值为IndexG_i,j+1＝(j+1)*3600+i；

G_i+1,j的索引值为IndexG_i,j+1＝j*3600+(i+1)；

G_i+1,j+1的索引值为IndexG_i+1,j+1＝(j+1)*3600+(i+1)；

由此，可知：

θ(G)_max＝max([RDE(2,IndexG_i,j),RDE(2,IndexG_i,j+1,RDE(2,IndexG_i+1,j)],RDE(2,IndexG_i+1,j+1))，

θ(G)_min＝min([RDE(2,IndexG_i,j),RDE(2,IndexG_i,j+1),RDE(2,IndexGi+1,j)],RDE(2,IndexGi+1,j+1))，

俯仰角计算方法：

β(G_i,j)＝hG_ij/RG_ij＝RDE(3,IndexG_i,j)/RDE(1,IndexG_i,j)，

计算完成所有网格的最大方位角θ(G)max、最小方位角θ(G)min和俯仰角β(G)后，将数据分别存储在AngleMax(1，size)，AngleMin(1，size)和Beta(1，size)三个矩阵变量中。

本方法首次提出了网格最大方位角和最小方位角的概念，基于这两个方位角，可更加准确描述网格的方位信息，消除了传统算法中假设网格大小一致带来的误差，也为后续更加准确的进行子区域划分、确定影响通视性的网格集合等操作提供了基础。

S4，利用网格方位信息，将DEM高程数据按照方位划分为M个子区域，每个子区域包含网格的数量远小于总区域的数据量，在子区域内进行通视分析可以大大提高分析效率。

其中，分区的方法包括以下步骤：

S4.1，计算每个子区域的起始方位角和结束方位角，见公式(5)

zone(i)_start＝(i-1)*(360/N)

zone(i)_end＝i*(360/N) (5)

其中，zone(i)代表第i个区域，zone(i)_start表示该区域的起始方位，zone(i)_end表示该区域的结束方位，如图5所示。

S4.2，根据区域的起始、结束方位角zone(i)_start、zone(i)_end，和网格的最大方位角θ(G)_max、最小方位角θ(G)_min，判断网格G是否属于该区域，当网格的最大、最小方位角满足公式(6)(7)(8)中的任一组条件，则判定该网格属于第i个子区域zone(i)：

本实施例中取N＝5760,根据公式(6)(7)(8)将网格进行分区，定义矩阵变量ZoneIndex(i)，ZoneIndex(i)中存储了第i个子区域中所有网格的索引值。通过索引值就可以访问到该子区域内所有网格的信息，例如，RDE(1,ZoneIndex(6))表示了第6个子区域网格的距离信息，其它子区域网格信息可以用同样方法获取。

本方法提出了一种根据网格方位角信息进行自由划分区域的方法，在划分的子区域内进行通视分析，将提高通视分析的速度。

S5，计算某一网格上空高度Δh处的通视情况。

步骤5中网格G_ij上空Δh处的通视情况的方法包括以下步骤：

S5.1，在柱坐标空间下的DEM高程数据中，选择一个网格G_ij，其坐标为(R_Gij，θ_Gij，h_Gij)，选择该网格上空某一点作为目标点V_T(R_Gij，θ_Gij，h_T)，其中h_T＝h_Gij+Δh，Δh表示目标点相对G_ij高度；

S5.2，根据G_ij的方位角，选择一个DEM高程数据子区域，利用该子区域内的DEM高程数据进行通视分析，如果zone(i)_start≤θ_Gij≤zone(i)_end，则选择第i个子区域；

S5.3，寻找影响目标通视性的网格集合{G}，如图5所示。V_T在水平面投影为V_T'，线段V_PV_T'在柱坐标下表示为满足{R≤R_T,θ＝θ_T}点的集合，因此{G}就是涵盖V_T'所在方位，且水平距离小于R_T的网格的集合，即G满足条件(9)：

公式(9)利用了网格最大方位角、最小方位角和距离信息，准确找到影响通视性的网格集合{G}，避免了其它算法大量的数值计算，也不会受到网格大小不一致的影响；

S5.4，通过步骤S5.3找到的影响通视的所有网格{G}，比较目标点俯仰角β(V_T)与{G}中所有网格的俯仰角大小，根据公式(10)判断是否通视：

see(V_T)表示目标点是否通视，通视记为1，不通视记为0，n表示集合{G}中包含了n个网格，如果分析空域高度是相对地形的高度，则β(V_T)＝arctan((h_Gij+Δh)/R_T)，Δh表示目标点相对G_ij高度；如果分析的空域高度是相对观测点高度，则β(V_T)＝arctan(h_T/R_T)，h_T表示目标点相对V_P的高度，如图6所示。

具体地，令G_ij上空Δh处V_T的坐标为(RGij，θGij，hT)，其中h_T＝hG_ij+Δh，Δh表示目标点相对网格G_ij高度。以网格G(3600)(3600)为例，其坐标为(100400m，51.7072°，-1263m)。假设分析地形上空200m的空域的通视性，令Δh＝200m。则目标点V_T的坐标为(100400m，51.7072°，-1063m)。

根据V_T的方位θ_T＝51.7072°，利用不等式zone(i)_start≤θ_Gij≤zone(i)_end确定其所在子区域为第828个子区域。该子区域网格点索引值存在ZoneIndex(828)的矩阵变量中，该子区域包含了5040个网格。

根据公式(9)在第828个子区域中进行筛选，寻找影响该点通视性的所有网格集合{G}。影响VT通视性的网格集合{G}中包含了3003个网格，这3003网格地形绘制成柱坐标下三维图形，如图9所示。找到{G}后，计算目标V_T的俯仰角β(V_T)，β(V_T)＝arctan(h_T/R_Gij)＝-0.6131°，计算max{β(G)}＝0.1031°。比较β(V_T)和max{β(G)}大小，通过公式(10)可以判断，see(3600)(3600)＝0，即网格G(3600)(3600)上空200m处不可通视。

通过比较网格点方位和距离信息，可以快速便捷地找到影响目标点通视性的所有网格，算法实现简单且精度较高。

S6，利用步骤5计算的所有网格的通视情况，获取该地域上高度Δh的空域通视情况。

具体地，分析完所有网格上空Δh＝200m处的通视情况后，就完成了该区域地面200m高度空域的通视情况分析。本实施例中200m高度通视效果图，如图10所示。

该方法将DEM数据转换到柱坐标空间，首次提出了用最大方位角、最小方位角来描述一个DEM网格的方位信息，因此本方法中用水平距离、方位、高度、最大方位、最小方位五个参数，更加准确地描述了网格信息，从而可以获得更高的分析精度。本方法提出了以观测点为中心将DEM数据分成M个子区域，在子区域内进行通视分析的思路，该思路可以极大降低一次分析处理的网格单元数量，提高分析处理效率。该方法通过比较网格最大方位角、最小方位角、网格水平距离、目标点方位、目标点水平距离等信息，能准确获取观测点与目标点投影经过的所有地形网格，避免了网格大小不一致引起的误差问题，通视分析的精度较高。本方法在完成坐标转换等数据预处理后，仅通过比较运算即可实现通视分析，避免了大量数值计算，易于实现，适合大规模的数据分析。

本发明还揭示了一种计算机可读存储介质，介质上存有计算机程序，计算机程序运行后，执行上述的柱坐标空间下的高精度空域通视分析方法。

本发明还揭示了一种计算机系统，包括处理器、存储介质，存储介质上存有计算机程序，处理器从存储介质上读取并运行计算机程序以执行上述的柱坐标空间下的高精度空域通视分析方法。

本领域技术人员将进一步领会，结合本文中所公开的实施例来描述的各种解说性逻辑板块、模块、电路、和算法步骤可实现为电子硬件、计算机软件、或这两者的组合。为清楚地解说硬件与软件的这一可互换性，各种解说性组件、框、模块、电路、和步骤在上面是以其功能性的形式作一般化描述的。此类功能性是被实现为硬件还是软件取决于具体应用和施加于整体系统的设计约束。技术人员对于每种特定应用可用不同的方式来实现所描述的功能性，但这样的实现决策不应被解读成导致脱离了本发明的范围。

结合本文所公开的实施例描述的各种解说性逻辑板块、模块、和电路可用通用处理器、数字信号处理器(DSP)、专用集成电路(ASIC)、现场可编程门阵列(FPGA)或其它可编程逻辑器件、分立的门或晶体管逻辑、分立的硬件组件、或其设计成执行本文所描述功能的任何组合来实现或执行。通用处理器可以是微处理器，但在替换方案中，该处理器可以是任何常规的处理器、控制器、微控制器、或状态机。处理器还可以被实现为计算设备的组合，例如DSP与微处理器的组合、多个微处理器、与DSP核心协作的一个或多个微处理器、或任何其他此类配置。

结合本文中公开的实施例描述的方法或算法的步骤可直接在硬件中、在由处理器执行的软件模块中或在这两者的组合中体现。软件模块可驻留在RAM存储器、闪存、ROM存储器、EPROM存储器、EEPROM存储器、寄存器、硬盘、可移动盘、CD-ROM、或本领域中所知的任何其他形式的存储介质中。示例性存储介质耦合到处理器以使得该处理器能从/向该存储介质读取和写入信息。在替换方案中，存储介质可以被整合到处理器。处理器和存储介质可驻留在ASIC中。ASIC可驻留在用户终端中。在替换方案中，处理器和存储介质可作为分立组件驻留在用户终端中。

在一个或多个示例性实施例中，所描述的功能可在硬件、软件、固件或其任何组合中实现。如果在软件中实现为计算机程序产品，则各功能可以作为一条或更多条指令或代码存储在计算机可读介质上或藉其进行传送。计算机可读介质包括计算机存储介质和通信介质两者，其包括促成计算机程序从一地向另一地转移的任何介质。存储介质可以是能被计算机访问的任何可用介质。作为示例而非限定，这样的计算机可读介质可包括RAM、ROM、EEPROM、CD-ROM或其它光盘存储、磁盘存储或其它磁存储设备、或能被用来携带或存储指令或数据结构形式的合意程序代码且能被计算机访问的任何其它介质。任何连接也被正当地称为计算机可读介质。例如，如果软件是使用同轴电缆、光纤电缆、双绞线、数字订户线(DSL)、或诸如红外、无线电、以及微波之类的无线技术从web网站、服务器、或其它远程源传送而来，则该同轴电缆、光纤电缆、双绞线、DSL、或诸如红外、无线电、以及微波之类的无线技术就被包括在介质的定义之中。如本文中所使用的盘(disk)和碟(disc)包括压缩碟(CD)、激光碟、光碟、数字多用碟(DVD)、软盘和蓝光碟，其中盘(disk)往往以磁的方式再现数据，而碟(disc)用激光以光学方式再现数据。上述的组合也应被包括在计算机可读介质的范围内。

提供对本公开的先前描述是为使得本领域任何技术人员皆能够制作或使用本公开。对本公开的各种修改对本领域技术人员来说都将是显而易见的，且本文中所定义的普适原理可被应用到其他变体而不会脱离本公开的精神或范围。由此，本公开并非旨在被限定于本文中所描述的示例和设计，而是应被授予与本文中所公开的原理和新颖性特征相一致的最广范围。

尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (8)

1.一种柱坐标空间下的高精度空域通视分析方法，其特征在于，包括以下步骤：

S1，选择需要分析的地域，采集待分析的DEM高程数据；

S2，将所述DEM高程数据转换到柱坐标空间；

S3，计算所述DEM高程数据中的网格方位信息；

S4，利用所述网格方位信息，将所述DEM高程数据按照方位划分为M个子区域；

S5，计算某一网格上空高度Δh处的通视情况；

S6，利用步骤5计算的所有网格的通视情况，获取该地域上高度Δh的空域通视情况。

2.根据权利要求1所述的方法，其特征在于，所述DEM高程数据存储的高度值H_ij表示网格G_ij的高程，其经纬度坐标为(B_i，L_j)，约定左下角顶点V_ij的(B_i，L_j，H_ij)坐标代表网格G_ij的坐标，所述步骤2将待分析地域的所述DEM数据经过三次坐标转换，转换到以地理坐标为(B_P，L_P，H_P)的观测点V_P为中心的柱坐标空间下，具体要经过三次坐标转换，步骤包括：

S21，将所述DEM高程数据从地理坐标系转换到空间直角坐标系，转换方法见公式(1):

其中，N＝a/W，N为地球卯酉圈曲率半径，a为地球长半轴，a＝6378.137km，

e²＝(a²-b²)/a²，e为地球的第一偏心率，b为地球的短半径，b＝6356.7523141km，

转换到空间直角坐标系后，观测点V_P的坐标则变换为(X_P，Y_P，Z_P)；

S2.2，将所述DEM高程数据从空间直角坐标系转换到以V_P为中心的站心坐标系，转换方法见公式(2)：

其中，(B_P，L_P，H_P)为观测点在地理坐标系下的坐标，(x，y，z)是所述DEM高程数据在站心坐标系下的坐标，此时观测点V_P的坐标为(0，0，0)；

S2.3，将所述DEM高程数据转换到柱坐标空间(R，θ，h)上，R表示某一网格距离观测点的水平距离，θ表示网格相对于观测点方位，以正北为方位0点，h表示网格相对于观测点水平面的高度，经过坐标转换后，h已包含地球曲率的影响，转换方法见公式(3)：

h＝z (3)。

3.根据权利要求1所述的方法，其特征在于：所述步骤3中的所述网格方位信息包括所述DEM高程数据中的每一个网格相对于观测点的最大方位角θ(G)_max、最小方位角θ(G)_min和俯仰角β(G)。

4.根据权利要求3所述的方法，其特征在于：第i行第j列网格中的所述网格方位信息分别为θ(G_ij)_max、θ(G_ij)_min、和β(G_ij)，其定义见公式(4)：

θ(G_ij)_max＝max([θ_i,jθ_i,j+1θ_i+1,jθ_i+1,j+1])

θ(G_ij)_min＝min([θ_i,jθ_i,j+1θ_i+1,jθ_i+1,j+1])

β(G_ij)＝arctan(h_ij/R_ij) (4)

其中，θ_i,j，θ_i,j+1，θ_i+1,j，θ_i+1,j+1表示网格G_ij的四个顶点的方位角。

5.根据权利要求3所述的方法，其特征在于，步骤4中分区的方法包括以下步骤：

S4.1，计算每个子区域的起始方位角和结束方位角，见公式(5)

zone(i)_start＝(i-1)*(360/N)

zone(i)_end＝i*(360/N) (5)

其中，zone(i)代表第i个区域，zone(i)_start表示该区域的起始方位，zone(i)_end表示该区域的结束方位；

S4.2，根据区域的起始、结束方位角zone(i)_start、zone(i)_end，和网格的最大方位角θ(G)_max、最小方位角θ(G)_min，判断网格G是否属于该区域，当网格的最大、最小方位角满足公式(6)(7)(8)中的任一组条件，则判定该网格属于第i个子区域zone(i)：

6.根据权利要求3所述的方法，其特征在于，所述步骤5中网格G_ij上空Δh处的通视情况的方法包括以下步骤：

S5.1，在柱坐标空间下的所述DEM高程数据中，选择一个网格G_ij，其坐标为(R_Gij，θ_Gij，h_Gij)，选择该网格上空某一点作为目标点V_T(R_Gij，θ_Gij，h_T)，其中h_T＝h_Gij+Δh，Δh表示目标点相对G_ij高度；

S5.2，根据G_ij的方位角，选择一个所述DEM高程数据子区域，利用该子区域内的所述DEM高程数据进行通视分析，如果zone(i)_start≤θ_Gij≤zone(i)_end，则选择第i个子区域；

S5.3，寻找影响目标通视性的网格集合{G}，V_T在水平面投影为V_T'，线段V_PV_T'在柱坐标下表示为满足{R≤R_T,θ＝θ_T}点的集合，因此{G}就是涵盖V_T'所在方位，且水平距离小于R_T的网格的集合，即G满足条件(9)：

公式(9)利用了网格最大方位角、最小方位角和距离信息，准确找到影响通视性的网格集合{G}，避免了其它算法大量的数值计算，也不会受到网格大小不一致的影响；

S5.4，通过步骤S5.3找到的影响通视的所有网格{G}，比较目标点俯仰角β(V_T)与{G}中所有网格的俯仰角大小，根据公式(10)判断是否通视：

see(V_T)表示目标点是否通视，通视记为1，不通视记为0，n表示集合{G}中包含了n个网格，如果分析空域高度是相对地形的高度，则β(V_T)＝arctan((h_Gij+Δh)/R_T)，Δh表示目标点相对G_ij高度；如果分析的空域高度是相对观测点高度，则β(V_T)＝arctan(h_T/R_T)，h_T表示目标点相对V_P的高度。

7.一种计算机可读存储介质，其特征在于：介质上存有计算机程序，计算机程序运行后，执行如权利要求1至6中任一项所述的柱坐标空间下的高精度空域通视分析方法。

8.一种计算机系统，其特征在于：包括处理器、存储介质，存储介质上存有计算机程序，处理器从存储介质上读取并运行计算机程序以执行如权利要求1至6中任一项所述的柱坐标空间下的高精度空域通视分析方法。

Images