一种基于CPU和GPU混合的高速并行地形遮蔽计算方法与流程

本发明属于数据分析领域，涉及一种基于cpu和gpu混合的高速并行地形遮蔽计算方法。

背景技术：

三维场景模拟仿真与地形数据息息相关，而对地形遮蔽的计算能够很好的分析一定范围数据与地形之间的关系，在军事、国土、测绘等方向具有很好的实际意义。

传统的地形遮蔽算法的实现过程常使用cpu编程的方法。cpu并行计算最常用的方法是多线程并行运算，通常利用openmp对循环进行并行化处理，openmp提供了一系列应用程序接口，通过在源代码中加入“pragma”来标明并行化，并通过系统级的抽象降低了并行编码的复杂性。openmp以fork-join模型执行。但是cpu编程的方法在大量数据情况下计算效率不高。

目前，在大规模的并行计算中，越来越多地应用图形处理器（gpu）进行并行计算。运算平台cuda(computeunifieddevicearchitecture)是一种操作gpu计算的硬件和软件架构，它将gpu视作一个数据并行计算设备，而且无需把这些计算映射到图形api。但在使用gpu进行并行计算的过程中，经常会遇到服务器或者工作站没有显卡或者显卡不支持cuda的情况。

技术实现要素：

本发明的目的在于采用普通cpu并行算法和cuda并行计算框架混合的方法，解决地形遮蔽计算过程中对环境适应能力不足、计算效率低下的问题。

为实现上述技术目的，达到上述技术效果，本发明通过以下技术方案实现：

一种基于cpu和gpu混合的高速并行地形遮蔽计算方法，该方法包括以下步骤：

步骤1）建立网格化数字高程模型，以高程矩阵方式来表示地面高程，根据需要计算的区域和高程的读取范围，将该区域划分成若干网格，对网格进行编码，记录每个网格的高程值，计算网格化的高程矩阵，使得区域内任一位置有对应的网格编码，从而得到该位置的高程值；

步骤2）采用遮蔽分析算法对区域的地形数据进行分析，将两点在高程矩阵中的编码拟合成一条直线，通过计算直线与高程网格矩阵的交点，得到视点到观察点的地形剖面，根据点与点之间的夹角判断是否通视，计算出最远通视点；

步骤3）高程数据校正，将影响高程数据的地球曲率和大气衰减因素作为校正参数，构建校正公式对高程数据进行校正。

进一步的，所述步骤1）中，将需要计算的区域划分成网格m×n，网格编码以起始点的地理坐标为网格的原点(0，0)，水平方向为0到n-1，垂直方向为0到m-1。

进一步的，所述步骤1）中，计算网格化的高程矩阵，包括以下步骤：

步骤1.1）根据需要计算的区域的网格大小（m×n），最小点的位置pmin(lon,lat),其中lon和lat分别表示经度和纬度，最大点的位置pmax(lon,lat),计算点p(lon,lat)对应的网格编码（x,y）,其计算公式如下：

单位网格的经纬度分辨率为：

r(∆lon,∆lat)=pmax-pmin/(n,m)；

步骤1.2）根据单位网格的经纬度分辨率和起始点的经纬度计算出对应点的网格编码，计算公式如下：

（x,y）=(p-pmin)/r。

进一步的，所述步骤2）中，采用地形遮蔽算法计算具体位置，包括以下步骤：

步骤2.1）计算视点p0(x0,y0,z0)在高程矩阵中的编码(m0,n0)，得到该位置的高程值z0，计算观察点p1(x1,y1,z1)在高程矩阵中的编码(m1,n1)，得到该位置的高程值z1；

步骤2.2）将p0点和p1点在高程矩阵中的编码拟合成一条直线h[x]，该直线以高程矩阵的行数为纵坐标y轴，以高程矩阵的列数为横坐标x轴，其直线公式为：

h(x)-m0=(m1-m0)/(n1-n0)×(x-n0)；

步骤2.3）直线h[x]为视点到观察点在高程矩阵上的投影，通过计算直线与高程网格矩阵的交点，得到对应的网格的高程数组z[k]，得到视点到观察点的地形剖面{x,y,z[k]}；

步骤2.4）计算视点与高程数组z[k]中每个点的与x轴的夹角beta，若tan(beta)大于直线h[x]的斜率则通视，否则不通视，该点为需要计算的最远通视点，根据最远通视点在直线h[x]上的位置，计算出其具体的地理位置。

进一步的，所述步骤3）中，采用以下校正公式对高程数据进行校正：

z=z0+d²(r-1)/d，

其中：

z为计入大气折射影响后的校正高程，

z0为观察位置的表面高程，

d为观测要素和观测位置之间的平面距离，

d为地球的直径，

r为光的折射系数。

进一步的，光的折射系数r默认值取0.13，适用于白天天空晴朗条件的标准气压下、高程范围为40到100米的位置，采用不同的折射值以模拟各种大气条件和高程对可见性的影响。

本发明的有益效果是:

本发明采用cpu和gpu混合的模式，通过对现有地形遮蔽算法的改造与优化，实现了对多种应用场景下的支持，如：低空、超低空和山地等等，大大提高了算法的计算效率。

附图说明

图1为地形遮蔽示意图；

图2为本发明的高程网格图。

具体实施方式

下面将参考附图并结合实施例，来详细说明本发明。

一种基于cpu和gpu混合的高速并行地形遮蔽计算方法，该方法包括以下步骤：

步骤1）建立网格化数字高程模型，以高程矩阵方式来表示地面高程，数字高程模型简称dem，是通过一组有序数值阵列形式表示地面高程的一种地面模型，本发明读取高程模型采用的是高程矩阵的方法，根据需要计算的区域和高程的读取范围，将该区域划分成若干网格，对网格进行编码，记录每个网格的高程值，计算网格化的高程矩阵，使得区域内任一位置有对应的网格编码，从而得到该位置的高程值；

步骤2）采用遮蔽分析算法对区域的地形数据进行分析，通过高程矩阵将点对点的地形遮蔽问题转化为离散空间中的直线与某一地形剖面的相交问题，本实施例中，将两点在高程矩阵中的编码拟合成一条直线，通过计算直线与高程网格矩阵的交点，得到视点到观察点的地形剖面，根据点与点之间的夹角判断是否通视，计算出最远通视点；

步骤3）高程数据校正，由于地球是一不规则椭圆，简单起见可将地球简化为一个半径为6371km的圆；其中，水准面是指处处与铅垂线垂直的连续封闭曲面；若使用水平面代替水准面，受地球曲率影响，必然会导致误差；所以，在进行地形遮蔽计算时，若计算的点为低空或超低空目标点时，则需要考虑地球曲率对计算结果的影响；而受天气的气压，密度、湿度、温度等变化，也会对计算结果产生误差；因此，本发明利用公式对高程数据进行校正，从而减小地球曲率和大气衰减对通视分析的影响，本实施例中，将影响高程数据的地球曲率和大气衰减因素作为校正参数，构建校正公式对高程数据进行校正。

所述步骤1）中，将需要计算的区域划分成网格m×n，网格编码以起始点的地理坐标为网格的原点(0，0)，水平方向为0到n-1，垂直方向为0到m-1。

所述步骤1）中，计算网格化的高程矩阵，包括以下步骤：

步骤1.1）根据需要计算的区域的网格大小（m×n），最小点的位置pmin(lon,lat),其中lon和lat分别表示经度和纬度，最大点的位置pmax(lon,lat),计算点p(lon,lat)对应的网格编码（x,y）,其计算公式如下：

单位网格的经纬度分辨率为：

r(∆lon,∆lat)=pmax-pmin/(n,m)；

步骤1.2）根据单位网格的经纬度分辨率和起始点的经纬度计算出对应点的网格编码，计算公式如下：

（x,y）=(p-pmin)/r。

具体计算结果如图2所示。

所述步骤2）中，采用地形遮蔽算法计算具体位置，如图1所示，包括以下步骤：

步骤2.1）计算视点p0(x0,y0,z0)在高程矩阵中的编码(m0,n0)，得到该位置的高程值z0，计算观察点p1(x1,y1,z1)在高程矩阵中的编码(m1,n1)，得到该位置的高程值z1；

步骤2.2）将p0点和p1点在高程矩阵中的编码拟合成一条直线h[x]，该直线以高程矩阵的行数为纵坐标y轴，以高程矩阵的列数为横坐标x轴，其直线公式为：

h(x)-m0=(m1-m0)/(n1-n0)×(x-n0)；

步骤2.3）直线h[x]为视点到观察点在高程矩阵上的投影，通过计算直线与高程网格矩阵的交点，得到对应的网格的高程数组z[k]，得到视点到观察点的地形剖面{x,y,z[k]}；

步骤2.4）计算视点与高程数组z[k]中每个点的与x轴的夹角beta，若tan(beta)大于直线h[x]的斜率则通视，否则不通视，该点为需要计算的最远通视点，根据最远通视点在直线h[x]上的位置，计算出其具体的地理位置。

所述步骤3）中，采用以下校正公式对高程数据进行校正：

z=z0+d²(r-1)/d，

其中：

z为计入大气折射影响后的校正高程，

z0为观察位置的表面高程，

d为观测要素和观测位置之间的平面距离，

d为地球的直径，本实施例中d取值为12740千米，

r为光的折射系数。

光的折射系数r默认值取0.13，适用于白天天空晴朗条件的标准气压下、高程范围为40到100米的位置，采用不同的折射值以模拟各种大气条件和高程对可见性的影响。

本实施例中，采用三种方法混合进行地形遮蔽计算，来实现相同地区，同一分辨率下，相同观察点个数进行遮蔽分析。方法一为使用cpu进行遮蔽计算；方法二使用cpu线程进行遮蔽计算；方法三使用cuda进行遮蔽计算。

实验过程中，计算需要计算地形遮蔽的点为1024个，计算区域的高程网格大小为2048*2048，三种方法的计算效率比较如下:

由上述结果进行对比分析可以看出，在地形遮蔽计算中，优先使用gpu进行并行运算，当gpu不支持cuda时或无gpu时，自动退化至cpu并行计算，当底层系统不支持openmp时，自动退化至cpu单线程计算。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。