CN110598159A - 基于有效计算量的局部型栅格空间分析并行计算方法 - Google Patents

Abstract

本发明涉及一种基于有效计算量的局部型栅格空间分析并行计算方法，包括以下步骤：并行环境初始化，读取源栅格图层数据，将并行进程分为主进程和若干个从进程，主进程根据有效计算量将栅格图层划分个多个栅格分块并执行并行计算过程，最后输出计算结果。本发明能够充分利用并行计算资源、提高算法并行效率。

Description

基于有效计算量的局部型栅格空间分析并行计算方法

技术领域

本发明涉及一种基于有效计算量的局部型栅格空间分析并行计算方法，属于并行计算技术领域。

背景技术

局部型栅格数据空间分析通常是在栅格图层中对单个或多个栅格单元开辟矩形分析窗口，并在该窗口内进行诸如极值、均值等一系列统计计算，具有典型的数据密集性的算法特征。局部型栅格空间分析通常包括求取均值、极值、中间值、众数、整体插值、空间统计等，从而实现图像增强、特征信息提取等应用。局部型栅格空间分析具体实现原理不同，但具有相同的算法特征，即其基本计算单元为单个栅格单元或规则分析窗口包含的栅格单元，并循环对不同基本单元调用相同规则进行计算。因此，该类型算法具有良好的并行性，主要表现在：算法原理简单、重复性高，适合使用并行计算技术进行处理；此外，计算过程限于规则、有限的分析窗口，因而可将栅格图层划分成计算独立、不需通讯的栅格分块，有利于并行计算。

对局部型栅格数据空间分析算法，传统的并行方法主要包括数据划分方法和并行调度方法。通常数据划分方法包括规则格网划分及四叉树划分等；这些方法主要基于空间范围、以保证面积相等为准则，实现对栅格数据的规则划分。尽管已有较多研究应用上述并行方法实现了对局部型栅格空间分析算法的并行加速，但仍存在一定的效率提升空间，主要表现在：现有的数据划分方法以面积相同为标准划分不同的栅格分块，对栅格图层的划分均较为粗略，未考虑划分后不同栅格分块间的实际计算量差异，从而容易引起并行过程中的数据倾斜。

发明内容

本发明要解决技术问题是：提供一种能够充分利用并行计算资源、提高算法并行效率的局部型栅格空间分析计算方法。

为了解决上述技术问题，本发明提出的技术方案是：一种基于有效计算量的局部型栅格空间分析并行计算方法，包括以下步骤：

步骤一、并行环境初始化，接收并行参数，所述并行参数包括并行进程个数p、两个栅格分块之间的计算量差异比例阈值s；

步骤二、读取源栅格图层数据，并获取该栅格图层的栅格行数w、栅格列数u和栅格尺寸；

步骤三、将并行进程分为主进程和p-1个从进程，主进程根据预设的具有nr×nc个栅格单元的分析窗口按栅格行的方向或栅格列的方向进行数据划分，将栅格图层划分为p-1个具有相同非空值栅格单元数量的栅格分块；其中按栅格行的方向进行数据划分的方法具体如下：

S301、统计栅格图层中的非空值栅格单元数量，即

式中，N_t为栅格图层中的非空值栅格单元数量，N_i为第i行的非空值栅格单元个数；

这样每个栅格分块包含的理论非空值栅格单元数量C_t＝N_t/(p-1)；

S302、按照栅格行的方向根据面积相同的原则确定各栅格分块的接缝线初始位置，任一栅格分块包含的栅格行数为v＝w/(p-1)、栅格列数为u；

S303、依次统计每个栅格分块内部各栅格行的非空值栅格单元数量，对其进行累加得到该栅格分块的实际计算量C_p；

S304、若C_p<C_t，则将接缝线向着使栅格行数增加的方向移动，直到满足以下条件时移动终止：

若C_p>C_t，则将接缝线向着使栅格行数减少的方向移动，直到满足以下条件时移动终止：

其中，step为预设的移动步长，即接缝线每次移动的栅格行数；C_step为待移动的step个栅格行包含的非空值栅格单元数量；

而按栅格列的方向进行数据划分的方法，除接缝线移动方向不同以外，均与按栅格列的方向进行数据划分的方法一致；

步骤四、主进程将p-1个栅格分块依次分配给p-1个从进程，各从进程并行处理各自分配的栅格分块，并在处理完后将计算结果反馈回主进程；

步骤五、主进程接收并行从进程的反馈信息，将计算结果写入目标文件，从而得到完整的并行计算结果。

首先需要说明的是，“计算量”和“有效计算量”的含义相同，本发明中以非空值栅格单元数量作为度量栅格有效计算量的标准。为了便于理解，本发明将非空值栅格单元数量视为计算量或有效计算量。

在并行计算时，传统的规则数据划分方法以面积相同为标准划分不同的栅格分块，从而忽略了各栅格分块的有效计算量，这使得划分后的结果并不能代表实际的计算量，从而不能实现并行计算过程中的负载均衡。本发明通过考虑计算过程中的有效计算量，并通过动态调整划分接缝线位置，以保证划分后的各栅格分块包含的有效计算量大致相当，从而实现负载均衡。

上述技术方案的进一步改进是：若步骤S304中某个栅格分块通过移动接缝线无法满足|C_p-C_t|<s×N_t，则对单个接缝线进行精确调整，方法如下：

1)当C_p>C_t时，需要减少栅格分块中的栅格单元个数N_a，当C_p<C_t时，需要增加栅格分块中的栅格单元个数N_a，其中N_a的取值范围为：

|C_t-C_p|-s×N_t＜N_a＜|C_t-C_p|+s×N_t，

在精确调整过程中需要调整的最少栅格单元个数为其取值应满足：

需要调整以step×step为大小的窗口个数N_dg为：

2)若C_p>C_t，则在栅格行内从首个非空值栅格单元开始依次向栅格分块内部减少step²个栅格单元，直至减少的窗口个数达到N_dg为止；

若C_p<C_t，则在栅格行内从首个非空值栅格单元开始依次向栅格分块内部增加包含step²个的栅格单元的规则窗口，直至增加的窗口个数达到N_dg为止。

这样本发明就突破了数据划分时栅格分块接缝线为直线段的限制，通过将接缝线设置为不规则曲线段，以更有效地实现不同栅格分块间有效计算量的平衡，以实现对并行计算资源的充分利用。

上述技术方案的再进一步改进是：步骤三中主进程在数据划分完成后，若step<nr或step<nc，对各划分后的栅格分块建立数据缓冲区，即在每个栅格分块基础上预读若干栅格行和栅格列，以独立地完成本栅格分块内的计算，从而避免在并行计算过程中与其他栅格分块进行通讯；其中，位于尾端的栅格分块不需要建立数据缓冲区，对于位于其他位置的栅格分块，在分块底部建立由宽度为nr–step个栅格行、长度为u个栅格列组成的不规则数据缓冲区。

上述技术方案的更进一步改进是：在步骤S302中划分栅格分块后，对任一栅格分块，若v％step＝0，则保持接缝线位置不变；否则，调整接缝线的位置，使接缝线向着栅格行数减少的方向移动v％step行，从而使得栅格行数为step的整数倍

附图说明

下面结合附图对本发明作进一步说明。

图1是本发明实施例中数据源的示意图。

图2是本发明实施例中不规则数据划分方法示意图。

图3是栅格分块数据缓冲区的原理示意图。

图4是并行计算结果融合的原理示意图。

图5是本发明实施例的并行k-means算法分类前后的对比图。

具体实施方式

实施例

本实施例的基于有效计算量的局部型栅格空间分析并行计算方法，包括以下步骤：

步骤一、并行环境初始化，接收并行参数，所述并行参数包括并行进程个数p、两个栅格分块之间的计算量差异比例阈值s。

步骤二、读取源栅格图层数据，并获取该栅格图层的栅格行数w、栅格列数u和栅格尺寸。

步骤三、将并行进程分为主进程和p-1个从进程，主进程根据预设的具有nr×nc个栅格单元的分析窗口按栅格行的方向或栅格列的方向进行数据划分，将栅格图层划分为p-1个具有相同非空值栅格单元的栅格分块。

本实施例在数据划分时，首先需要确定不同算法类型包含的数据粒度。根据不同算法原理，可将局部型栅格数据空间分析中的数据粒度分为单个栅格单元及分析窗口中包含的多个栅格单元。通常，计算过程中的分析窗口大小可认为是包含nr×nc个栅格单元的矩形窗口；因此，当数据粒度为单个栅格单元时，可认为其分析窗口大小为1×1。在本例中，包含nr×nc个栅格单元的分析窗口即为并行计算中的数据粒度。进而，局部型栅格数据空间分析过程可概括为：从栅格数据图层左上角开始，按照先栅格行、后栅格列的顺序逐数据粒度移动；在各数据粒度内按照计算规则对nr×nc个栅格单元进行处理，从而完成计算。针对不同的算法类型，数据粒度的移动步长step也不相同。当移动步长满足step＜nr或step＜nc时，两个相邻数据粒度包含重叠的栅格单元；当移动步长满足step≥nr且step≥nc时，两个相邻数据粒度包含的栅格单元均互不相同。

本实施例在确定了数据粒度、移动步长后，即可实现对源栅格数据图层进行数据划分。针对栅格数据的划分，其本质是确定各栅格分块的接缝线位置；不同的接缝线位置决定了栅格分块包含计算量的大小。因此，本例主要围绕确定栅格分块接缝线展开，包含三个步骤：确定接缝线初始位置、粗略调整接缝线位置及精确调整接缝线位置。考虑到数据划分后栅格分块接缝线边界往往引起不同分析窗口的割裂，从而所有接缝线的长度应尽可能短，以减少被割裂栅格单元的数目。因此，对于包含w个栅格行、u个栅格列的栅格图层，在划分过程中，若栅格图层中栅格行数大于栅格列数(w>u)，则按照栅格行的方向进行划分；若栅格图层中栅格列数大于栅格行数(u>w)，则按照栅格列的方向进行划分。本文以按照栅格行方向进行划分的方式(w>u)为例，说明不规则划分方法的主要过程。数据划分过程中，栅格分块数设置为p-1、接缝线可用sl表示，则sl_i(i>0且i<n)即表示第i个栅格分块所属的接缝线；两个栅格分块之间的计算量差异阈值比例设置为s，且以非空值的栅格单元个数作为度量栅格有效计算量的标准。

其中，按栅格行的方向进行数据划分的方法具体如下：

S301、统计栅格图层中的非空值栅格单元数量，即

式中，N_t为栅格图层中的非空值栅格单元数量，N_i为第i行的非空值栅格单元个数；

每个栅格分块包含的理论非空值栅格单元数量C_t＝N_t/n。

S302、按照栅格行的方向根据面积相同的原则确定各栅格分块的接缝线初始位置，任一栅格分块包含的栅格行数为v＝w/(p-1)、栅格列数为u。这样，使得任一栅格分块包含的面积为整个栅格图层的1/(p-1)，即确定了接缝线的初始位置如图2a所示。

S303、依次统计每个栅格分块内部各栅格行的非空值栅格单元数量，对其进行累加得到该栅格分块的实际计算量C_p；通常，栅格分块的实际计算量与理论计算量相差较大，从而需要循环调整接缝线的位置，以使得实际计算量与理论计算量相当。

S304、如图2b所示，若C_p<C_t，则将接缝线向着使栅格行数增加的方向移动，直到满足以下条件时移动终止：

若C_p>C_t，则将接缝线向着使栅格行数减少的方向移动，直到满足以下条件时移动终止：

其中，step为预设的移动步长，即接缝线每次移动的栅格行数；C_step为待移动的step个栅格行包含的非空值栅格单元数量。

若步骤S304中某个栅格分块通过移动接缝线无法满足|C_p-C_t|<s×N_t，则对单个接缝线进行精确调整，方法如下：

1)当C_p>C_t时，需要减少栅格分块中的栅格单元个数N_a，当C_p<C_t时，需要增加栅格分块中的栅格单元个数N_a，其中N_a的取值范围为：

|C_t-C_p|-s×N_t＜N_a＜|C_t-C_p|+s×N_t，

在精确调整过程中需要调整的最少栅格单元个数为其取值应满足：

需要调整以step×step为大小的窗口个数N_dg为：

2)若C_p>C_t，则在栅格行内从首个非空值栅格单元开始依次向栅格分块内部减少step²个栅格单元，直至减少的窗口个数达到N_dg为止；

若C_p<C_t，则在栅格行内从首个非空值栅格单元开始依次向栅格分块内部增加包含step²个的栅格单元的规则窗口，直至增加的窗口个数达到N_dg为止。

当完成上述过程后，即完成了接缝线的精确调整过程，如图2c所示。在栅格分块接缝线精确调整前，接缝线的形状为直线段，在精确调整后，接缝线的形状将变成不规则的弯曲线段，如图2d所示。

若step<nr或step<nc，对各划分后的栅格分块建立数据缓冲区，即在每个栅格分块基础上预读若干栅格行和栅格列，以独立地完成本栅格分块内的计算，从而避免在并行计算过程中与其他栅格分块进行通讯；其中，位于尾端的栅格分块不需要建立数据缓冲区，对于位于其他位置的栅格分块，在分块底部建立由宽度为nr–step个栅格行、长度为u个栅格列组成的不规则数据缓冲区。

这是因为移动步长与数据粒度大小满足step≥nr且step≥nc时，则在完成不规则数据划分后，各栅格分块间的计算相互独立；若数据粒度移动步长与数据粒度大小满足step<nr或step<nc，则位于栅格分块接缝线边界处的数据粒度将不可避免地被割裂，需要获取其他栅格分块的栅格单元以完成计算。对于满足step<nr或step<nc的算法类型，需要在数据划分完成后对各划分后的栅格分块建立数据缓冲区，即在源栅格分块基础上预读若干栅格行和栅格列，以独立地完成本栅格分块内的计算，从而避免在并行计算过程中与其他栅格分块进行通讯，如图3所示。其中，位于尾端的栅格分块不需要建立数据缓冲区；对于位于其他空间位置的栅格分块，需在分块底部建立由宽度为nr–step个栅格行、长度为u个栅格列组成的不规则数据缓冲区。这样各栅格分块均保留了一定的数据缓冲区，从而避免了并行计算过程中的数据传输。

步骤四、主进程将栅格图层中p-1个栅格分块，依次分配给p-1个从进程，各从进程并行处理各自分配的栅格分块，并在处理完后将计算结果反馈回主进程。

步骤五、主进程接收并行从进程的反馈信息，将计算结果写入目标文件，从而形成完整的并行计算结果。

在并行处理过程完成后，各栅格分块的结果相互独立，从而需要进行各栅格分块的结果融合，从而形成最终的完整计算结果，如图4所示。考虑到各栅格分块分布相对规整，且相邻栅格分块的接缝线完全重合，因而可以利用各栅格分块左上角点坐标作为融合过程中该栅格分块在整个栅格图层中的标识点。在完成各栅格分块的计算后，通过计算该栅格分块的左上角点坐标并定位其存储位置，将计算结果顺序写入目标文件中。

为了验证本实施例的方法，发明人采用3幅中国江苏省南京市江宁区无人机遥感影像(如图1所示，分别为数据1、数据2和数据3)作为数据源进行并行计算实验。其中，数据1为江宁区完整遥感影像，数据2和数据3分别为江宁区部分区域的遥感影像。数据1数据量为53.1GB，包含115,209栅格行×122,425栅格列；数据2数据量为26.7GB，包含94,001栅格行×94,001栅格列；数据3数据量为6.9GB，包含45,001栅格行×50,401栅格列。各数据中均包含3个波段，空间分辨率为0.5米。数据源包含的土地利用类型主要有水体、裸地、建设用地、林地及耕地。

k-means算法是典型的局部型栅格数据空间分析算法类型，其功能是实现遥感影像的非监督分类，主要通过反复迭代并更新遥感影像的分类中心，从而达到最低的分类误差平方和。该算法的基本原理可表述为：对一个待分类的数据集X＝(x₁,x₂,…,x_n)，数据集中的每个待分类数据均是m维向量，k-means算法即是将这个数据集分为t(t≤n)个子集的过程。分类后的子集可表达为S＝(s₁,s₂,…,s_k)。算法的执行过程包含以下步骤：(1)随机从数据集中选取t个数据作为初始的遥感影像分类中心；(2)在第h次迭代中，对任一X中的i≠j(i＝1,2,…,t)，若则其中，是以为中心的类；(3)由步骤(2)得到的类即为新的分类中心可表达为

公式中，N_j为C_j类中的样本数目；通常是按照使得J取值最小的原则确定，J可表达为

(4)对所有的i＝1,2,…,h，均执行上述判断操作，直至则分类结束。

在k-means分类算法中，数据粒度为单个栅格单元，需逐个对非空值栅格单元判断其类别归属，具有典型的数据密集型算法特征。

本例首先实现了k-means算法的并行化，在该算法中，计算复杂度主要取决于遥感影像中的迭代次数L、分类数目t、非空值栅格单元个数n和波段数b。若分类中的基本时间单元为T，则单次迭代耗费的时间主要包括分类中心更新时间和迭代处理时间。其中，更新时间可表示为tb；迭代处理时间由像元个数、分类数目及波段数决定，可表示为f(n,k,b)。基于以上分析，k-means的完整计算时间T_s可计算为

T_s＝(f(n,k,b)+kb)LT

在并行算法中，执行时间主要由数据划分(T_gd)、类别归属计算(T_dc)、分类中心更新(T_cu)、类别汇总(T_sc)和数据读写(T_rw)组成，可进一步表示为

T_p＝T_gd+T_dc+T_cu+T_sc+T_rw。

为了验证本例实现的并行算法的分类精度，利用本例并行算法调用不同的进程，对实验数据3进行并行k-means分类，生成数据格式为ERDAS Imagine(*.img)格式的栅格结果数据。原始试验区的遥感影像部分放大区域如图5a所示，调用8个并行进程执行该并行算法，根据土地利用类型种类将待分类遥感影像分成5个类别，分类后的栅格结果图如图5b所示。对比分类前遥感影像和本例并行算法分类结果可以看出，两者在不同土地利用类型的对比上均具有较高的相似性，因而可看出本文提出的并行分类算法结果具有一定的正确性。

为了进一步验证本文实现的基于k-means遥感影像分类并行算法的正确性，将本例算法的分类结果与遥感影像分类专业软件ENVI中k-means算法的分类结果进行对比。具体来说，利用ENVI中的“k-means Classification”功能模块同样将本实验中的源数据分为五类，并将其计算结果与本例并行算法的分类结果进行对比，通过比较各类别栅格单元个数差异以验证本例方法的有效性，对比结果如表1所示。通过对比可以看出，本例算法与ENVI软件算法运行结果差异较小，取得了精度较高的分类结果。

在实验中，主要采用运行时间、加速比和负载均衡指数三个指标来评价算法的并行效率。本例实现的并行k-means算法主要采用主从式的并行模式实现对计算数据的并行动态调度，因此并行算法的最少进程数为2。实验中改变并行进程数从2增加至120，调用并行k-means算法执行数据3；同时，测试并行算法在不同进程数时的运行时间、加速比及对应的负载均衡指数。

从计算结果可以看出，并行算法的串行执行时间为2400.28秒。随着并行进程数的增加，并行运行时间逐渐减少；在该过程中，并行算法加速比逐渐增加，并与进程数呈近似线性变化关系。当进程数达到108时，并行时间达到最少，为118.42秒；此时，并行加速比达到峰值，为20.27。此后，尽管并行进程数进一步增加，但此时并行环境内的计算资源达到饱和，更多的进程数加剧了并行环境内的资源竞争；因此，并行时间达到稳定状态，并行加速比开始逐渐下降。上述结果表明，应用本文并行方法实现的并行算法可大大减少算法运行时间，取得了良好的并行加速比。

本发明不局限于上述实施例所述的具体技术方案，除上述实施例外，本发明还可以有其他实施方式。对于本领域的技术人员来说，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等形成的技术方案，均应包含在本发明的保护范围之内。

Claims (4)

1.一种基于有效计算量的局部型栅格空间分析并行计算方法，包括以下步骤：

步骤一、并行环境初始化，接收并行参数，所述并行参数包括并行进程个数p、两个栅格分块之间的计算量差异比例阈值s；

步骤二、读取源栅格图层数据，并获取该栅格图层的栅格行数w、栅格列数u和栅格尺寸；

步骤三、将并行进程分为主进程和p-1个从进程，主进程根据预设的具有nr×nc个栅格单元的分析窗口按栅格行的方向或栅格列的方向进行数据划分，将栅格图层划分为p-1个具有相同非空值栅格单元数量的栅格分块；其中按栅格行的方向进行数据划分的方法具体如下：

S301、统计栅格图层中的非空值栅格单元数量，即

式中，N_t为栅格图层中的非空值栅格单元数量，N_i为第i行的非空值栅格单元个数；

这样每个栅格分块包含的理论非空值栅格单元数量C_t＝N_t/(p-1)；

S302、按照栅格行的方向根据面积相同的原则确定各栅格分块的接缝线初始位置，任一栅格分块包含的栅格行数为v＝w/(p-1)、栅格列数为u；

S303、依次统计每个栅格分块内部各栅格行的非空值栅格单元数量，对其进行累加得到该栅格分块的实际计算量C_p；

S304、若C_p<C_t，则将接缝线向着使栅格行数增加的方向移动，直到满足以下条件时移动终止：

若C_p>C_t，则将接缝线向着使栅格行数减少的方向移动，直到满足以下条件时移动终止：

其中，step为预设的移动步长，即接缝线每次移动的栅格行数；C_step为待移动的step个栅格行包含的非空值栅格单元数量；

而按栅格列的方向进行数据划分的方法，除接缝线移动方向不同以外，均与按栅格列的方向进行数据划分的方法一致；

步骤四、主进程将p-1个栅格分块依次分配给p-1个从进程，各从进程并行处理各自分配的栅格分块，并在处理完后将计算结果反馈回主进程；

步骤五、主进程接收并行从进程的反馈信息，将计算结果写入目标文件，从而得到完整的并行计算结果。

2.根据权利要求1所述的基于有效计算量的局部型栅格空间分析并行计算方法，其特征在于：若步骤S304中某个栅格分块通过移动接缝线无法满足|C_p-C_t|<s×N_t，则对单个接缝线进行精确调整，方法如下：

1)当C_p>C_t时，需要减少栅格分块中的栅格单元个数N_a，当C_p<C_t时，需要增加栅格分块中的栅格单元个数N_a，其中N_a的取值范围为：

|C_t-C_p|-s×N_t＜N_a＜|C_t-C_p|+s×N_t，

在精确调整过程中需要调整的最少栅格单元个数为其取值应满足：

需要调整以step×step为大小的窗口个数N_dg为：

2)若C_p>C_t，则在栅格行内从首个非空值栅格单元开始依次向栅格分块内部减少step²个栅格单元，直至减少的窗口个数达到N_dg为止；

若C_p<C_t，则在栅格行内从首个非空值栅格单元开始依次向栅格分块内部增加包含step²个的栅格单元的规则窗口，直至增加的窗口个数达到N_dg为止。

3.根据权利要求1所述的基于有效计算量的局部型栅格空间分析并行计算方法，其特征在于：步骤三中主进程在数据划分完成后，若step<nr或step<nc，对各划分后的栅格分块建立数据缓冲区，即在每个栅格分块基础上预读若干栅格行和栅格列，以独立地完成本栅格分块内的计算，从而避免在并行计算过程中与其他栅格分块进行通讯；其中，位于尾端的栅格分块不需要建立数据缓冲区，对于位于其他位置的栅格分块，在分块底部建立由宽度为nr–step个栅格行、长度为u个栅格列组成的不规则数据缓冲区。

4.根据权利要求1所述的基于有效计算量的局部型栅格空间分析并行计算方法，其特征在于：在步骤S302中划分栅格分块后，对任一栅格分块，若v％step＝0，则保持接缝线位置不变；否则，调整接缝线的位置，使接缝线向着栅格行数减少的方向移动v％step行，从而使得栅格行数为step的整数倍。