CN104392147A - 面向区域尺度土壤侵蚀建模的地形因子并行计算方法 - Google Patents

Abstract

本发明公开了一种面向区域尺度土壤侵蚀建模的地形因子并行计算方法，属于计算机和信息技术领域。该方法包括以下过程：(1)确定数据划分策略，根据并行节点数目，采用行带缓冲区划分方法对原始地形数据进行等行划分，同时每个子数据块包含了数据区、缓冲区以及控制点三个要素；(2)按照数据划分策略，基于预处理后的DEM数据，进行坡度并行计算、水流方向计算、汇流累积计算、河网提取计算、坡长计算以及坡度坡长因子计算六步。本发明在并行算法设计过程中充分考虑了串行算法的计算特征，针对局部型算法和全局型算法分别设计了不同的并行化策略，有效提高了算法的并行效率，实现了基于海量地形数据的地形因子快速准确计算。

Description

面向区域尺度土壤侵蚀建模的地形因子并行计算方法

技术领域

本发明属于计算机和信息技术领域，特别涉及一种面向区域土壤侵蚀建模的地形因子并行计算方法。

背景技术

土壤侵蚀是全人类共同面对的重大环境问题，而地形要素对土壤侵蚀发生发展具有重要的控制作用，也是现有土壤侵蚀模型中不可或缺的因子之一。当前学者普通采用坡度坡长因子以表达地形要素对土壤侵蚀的影响。坡度坡长因子包含了坡度因子和坡长因子两部分，其中坡长的提取是一个难点。Hickey,Van Remortel，张宏鸣等学者基于数字高程模型产品(DigitalElevation Models，缩写为DEM)提出并发展了面向侵蚀过程的坡长提取方法，并已经在坡面尺度和流域尺度的土壤侵蚀预测中得到广泛应用。

近年来为服务于区域水土保持宏观决策、生态环境综合治理、全球变化等问题，在已有的坡面尺度、流域尺度土壤侵蚀研究的基础上，区域尺度土壤侵蚀相关研究得到越来越多的重视，相关学者开发了一系列的区域土壤侵蚀模型。在区域尺度土壤侵蚀预测与模拟中，由于大区域高精度数据较难获取，同时由于计算机存储、计算等能力的限制，现有研究通常采用以下两种方式求取适用于区域尺度土壤侵蚀定量评价的地形指标。一是根据侵蚀地貌学理论拟订替代指标,间接表示坡度坡长因子,如地形起伏度等。二是对基于中低分辨率DEM提取的地形因子进行尺度变换,以满足土壤侵蚀建模的需要。然而这两种方法都只是权宜之计，目前在区域土壤侵蚀建模中，基于高精度地形数据的坡度坡长因子高效计算已成为制约其应用与发展的关键问题。

随着地理数据采集方式的进步，特别是LiDAR、INSAR等技术的应用，大区域高精度DEM数据越来越丰富。数据来源的不断扩展更加凸显了现有计算方式的落后，运用并行计算方法实现对区域尺度的坡度坡长因子提取是解决这一问题的有效方法。目前并行计算在数字地形分析领域已经得到了广泛的应用，特别是在流域分析方面，基于并行计算提取河网、进行流域划分已有大量研究。然而，目前现有针对坡度坡长因子的计算还采用串行方法，并行计算在流域分析中的成功运用可以借鉴到坡度坡长因子计算中。

发明内容

本发明所要解决的技术问题在于，针对现有串行算法无法高效提取面向区域尺度土壤侵蚀建模的地形因子，提出了一种多核集群环境下坡度坡长因子并行计算方法，为区域尺度的土壤侵蚀预测与模拟提供技术支撑。

针对上述发明目的，本发明所采用的技术方案如下：

面向区域尺度土壤侵蚀建模的地形因子并行计算方法，包括以下过程：

(1)确定数据划分策略，其划分方式为：根据并行节点数目，采用行带缓冲区划分方法对原始地形数据进行等行划分，同时每个子数据块包含了三个要素：a)数据区：本计算节点所分配的数据区；b)缓冲区：每个子数据块附加了上下两个缓冲区，对于第一子数据块只有下缓冲区，对于最后一个数据块只有上缓冲区；c)控制点：用于定位每个数据块的起始读写位置；

(2)按照步骤(1)的数据划分策略，基于预处理后的DEM数据，进行坡度并行计算；

(3)按照步骤(1)的数据划分策略，基于预处理后的DEM数据，进行坡长并行计算；具体包含了以下几个子算法：a)基于预处理后的DEM数据，进行水流方向并行计算；b)基于水流方向计算结果进行汇流累积并行计算；c)基于汇流累积计算结果和设定的河网阈值，进行河网提取并行计算；d)基于水流方向计算结果进行坡长并行计算，同时此步骤需要坡度数据用于坡度截断栅格判断，河网水系数据用于河网截断栅格判断；

(4)按照步骤(1)的数据划分策略，基于坡度和坡长计算结果，进行坡度坡长因子并行计算。

上述步骤(2)具体的计算流程如下：

步骤21.根据数据划分策略，对输入的DEM数据进行数据划分，并创建结果数据文件；

步骤22.各子节点遍历输入的DEM数据，并执行坡度算法，计算完成后将计算结果写入结果文件。

上述步骤(3)具体的计算流程如下：

步骤31.根据数据划分策略，对输入的DEM数据进行数据划分，并创建结果数据文件；各子节点遍历输入的DEM数据，并执行水流方向算法，计算完成后将计算结果写入结果文件；

步骤32.根据数据划分策略，对水流方向矩阵进行数据划分，并创建汇流累积结果文件；各子节点遍历输入的水流方向矩阵，将源点栅格赋值为0同时入栈，非源点栅格赋值为-m，m为邻接汇入栅格数目，遍历完成后依次取出栈内栅格并沿流向计算其下游栅格汇流累积值，当栈内元素为空时计算结束；

步骤33.数据块间更新缓冲区信息，各子数据块已更新后的缓冲区栅格为源点栅格依次更新本数据区汇流累积值；一次更新计算完成后，遍历各子数据块边界行数据，若全部栅格完成计算则计算结束，否则需要重复步骤33；计算完成后各子节点将计算结果写入已创建的汇流累积文件；

步骤34.根据数据划分策略，对输入的汇流累积数据进行数据划分，并创建结果数据文件；各子节点遍历输入数据，并执行河网提取算法，计算完成后将计算结果写入结果文件；

步骤35.根据数据划分策略，分别对河网矩阵、坡度矩阵、水流方向矩阵进行数据划分，并创建坡长结果文件；各子节点遍历输入数据，根据侵蚀学坡长算法定义分别执行栅格点标识、非累积坡长计算、累积坡长计算步骤；

步骤36.数据块间更新缓冲区信息，各子数据块以更新后的缓冲区栅格为源点栅格依次更新本数据区的坡长值；一次更新计算完成后，遍历各子数据块边界行数据，若全部栅格完成计算则计算结束，否则需要重复步骤36；计算完成后各子节点将计算结果写入已创建的坡长文件。

上述步骤(4)具体的计算流程如下：

步骤41.根据数据划分策略，对输入的坡度、坡长数据进行数据划分，并创建结果数据文件；

步骤42.各子节点遍历输入数据，并执行坡度坡长因子算法，计算完成后将计算结果写入结果文件。

本发明具有以下两点优势：

(1)将现有的小流域尺度的土壤侵蚀模型地形因子计算方法扩展到了区域尺度，运用并行技术克服了现有串行算法面对大区域海量高精度数据处理能力的瓶颈，为区域尺度土壤侵蚀模拟提供有力支持。

(2)并行算法设计过程中充分考虑了串行算法的计算特征，针对局部型算法和全局型算法分别设计了不同的并行化策略，有效提高了算法的并行效率。实现了基于海量地形数据的地形因子快速准确计算。

附图说明

图1是本发明中数据划分示意图；

图2是全局型算法数据依赖示意图，图2(a)是水流方向矩阵，图2(b)是计算源点栅格示意图，图2(c)是数据分块导致的计算未完成示意图；

图3是本发明中缓冲区更新计算策略示意图，图3(a)是进程间更新缓冲区示意图，图3(b)是缓冲区更新后各进程分别计算示意图；

图4是本发明中算法并行化策略示意图；

图5是本发明中并行化流程示意图；

图6是实施例中输入DEM数据；

图7是实施例中输出的坡度坡长因子结果。

具体实施方式

下面结合附图和具体实施例，对本发明作进一步详细说明。

一、数据划分策略

本发明采用的是一种条带划分方式，其基本思想是根据并行算法的进程数，将原始数据划分为等数据行的子数据块，每个子数据块主要包含三个元素：(1)数据区：本计算节点所能分配的数据区域；(2)控制点：用于定位每个数据块的起始读写位置；(3)缓冲区：每个子数据块附加了上下两个缓冲区(对于第一子数据块只有下缓冲区，最后一个数据块只有上缓冲区)。缓冲区一方面可解决跨数据访问数据的问题，同时也可用于进程间计算结果的传递。(如图1所示)

二、算法并行策略

本发明涉及的六个子算法，根据其计算特征可分为局部型算法和全局型算法两类。对于局部型算法，每一个栅格的计算都是相对独立的，即计算过程中栅格之间的计算结果不具有依赖性。对于局部型算法并行化，数据划分后各子数据块之间不存在数据依赖，因此不需要数据块间的通信。而对于全局型算法，数据划分后各子数据块之间需要信息传递以消除依赖关系。

局部型算法可进一步划分为：(1)点运算，如河网提取和因子计算两类算法，每个栅格计算是完全独立的；(2)窗口运算，如坡度和水流方向两类算法，每个栅格计算时需要开辟一个分析窗口，在本发明中采用的是3*3窗口。具体实施时，局部型算法的并行化策略如下所示：

(1)程序初始化，获取进程数及算法相关参数。

(2)根据进程数拆分DEM数据。

(3)读取数据，各进程根据定位信息读取本节点数据及缓冲区数据。

(4)算法计算，根据具体算法对各进程的数据进行计算。

(5)写出数据，根据定位点信息，各进程分别将计算结果写入结果文件。

(6)释放内存，并行计算任务完成。

局部型算法的并行化策略的伪代码如表1所示：

表1 局部型算法的并行化策略的伪代码

不同于局部算法，全局型算法各栅格之间存在依赖关系，这种依赖关系导致全局型算法中每一个栅格无法保证独立计算。针对全局型算法的特征，在计算过程中，本发明定义了计算阀与计算源两个概念。其中计算阀是指待计算栅格可以开始计算的条件，满足了计算阀则表明待计算栅格已经满足了依赖关系。计算源是指在初始状态下即满足计算阀的栅格。如图2所示，当根据图2(a)水流方向矩阵计算汇流累积时，该计算任务的计算阀可定义为待计算栅格的上游栅格均已计算完成，则图2(b)中的灰色栅格在初始状态下即满足计算阀，可视为计算源。从计算源开始顺流计算，当其他栅格满足计算阀时，即可赋值相应的汇流累积值。然而，当此类全局型算法并行计算时，数据划分会导致具有依赖关系的栅格被分配在不同的子数据块，因此单纯的各进程分块执行算法会造成部分栅格无法计算。如图2(c)所示，分块计算后部分栅格无法满足计算阀，计算无法完成。

为此，本发明设计了缓冲区更新计算策略，其基本步骤如下：

(1)每个子数据块初始计算。每个栅格是否能一次完成计算取决于该栅格以及它的依赖栅格是否在相同的数据域内。正常情况下，各进程初始计算结束后仍有一定栅格未完成计算。

(2)进程间通信更新计算结果。考虑到数据间的计算依赖性具有传递性，因此本发明通过各子数据块的边界行状态判断该数据块是否全部完成计算。当边界行有栅格没有完成计算，则每个进程间更新缓冲区以传递边界行信息。(如图3(a))

(3)每个子数据块更新计算。各进程搜索更新后的缓冲区，如果缓冲区中有栅格已取得最终值且未被被访问，则该栅格被入栈bufferSource且其访问标识设置为1。搜索完成后，栈bufferSource中的每个栅格依次出栈，此时出栈栅格可视为计算源，并由此栅格开始更新未计算栅格。(如图3(b))

(4)更新计算完成后，如果各子数据块的边界行有还未计算完成的栅格，则重复步骤2和步骤3直到全部栅格均已计算完成。

全局型算法并行化的伪代码如表2所示：

表2 全局型算法的并行化策略的伪代码

根据以上讨论，本发明中的并行化策略如4所示。局部型算法和全局型算法的并行化策略的主要差异在于，全局型算法需要进程间的通信以传递各自计算结果，以消除数据依赖的影响，而其他部分均一致。

三、主要类设计

本发明共设计了8个类：RasterClass用于定义栅格DEM的数据基本格式，数据的统计信息等。ParatitionClass用于定义并行计算相关的的基本函数，包含了数据划分、数据融合、数据读写、进程间数据通信等；SlopeClass用于实现坡度因子的并行计算功能；FlowDirectionClass用于实现水流方向算法的并行计算功能；FlowAccumulationClass用于实现汇流累积算法的并行计算功能；StreamClass用于实现河网提取算法的并行计算功能；SlopelengthClass用于实现坡长并行计算功能；FactorClass用于实现坡度坡长因子并行计算功能。

本发明中具体技术流程如下：

第一步、坡度计算：在本发明中，坡度计算结果用于坡长计算中坡度截断栅格判断以及最终坡度坡长因子计算。本发明中采用的是三阶不带权算法，其计算函数如公式(1)-(3)所示。根据坡度算法算法的特征，采用局部型算法并行化策略进行并行计算。

$f_x=\frac{Z_{i-1,j+1}+Z_{i,j+1}+Z_{i+1,j+1}-Z_{i-1,j-1}-Z_{i,j-1}-Z_{i+1,j-1}}{6g}---\left(1\right)$

$f_y=\frac{Z_{i+1,j+1}+Z_{i+1,j}+Z_{i+1,j-1}-Z_{i-1,j-1}-Z_{i-1,j}-Z_{i-1,j+1}}{6g}---\left(2\right)$

$\mathrm{Slope}(i,j)=\arctan \sqrt{{f_x}^2+{f_y}^2}---\left(3\right)$

式中Z_i,j表示的是(i，j)栅格处的高程值，g为单位栅格大小。

第二步、水流方向计算：水流方向是指流水离开格网时的流向。流向算法目前主要有单流向和多流向两类。在本发明中采用的是单流向算法中的最大坡降算法，其基本原理是首先给定一个水流方向矩阵，规定了各方向的标识。然后计算中心栅格与领域的八个栅格的坡降，最大坡降方向即确定为水流方向，并赋予相应的方向标识。水流方向算法也属于局部型算法，因其并行化方法采用局部型算法并行化策略。

第三步、汇流累积计算：汇流累积矩阵是指流向该格网的所有上游格网单元的水流累积量，它是基于水流方向矩阵确定的。具体计算时，每一个栅格需要等待汇入该栅格的所有相邻栅格均已完成计算时才可以进行计算,其计算公式为：

$\mathrm{Fac}\left(c\right)=\mathrm{Num}+\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}\mathrm{Fac}\left(i\right)---\left(4\right)$

式中Num是汇入该栅格的个数，i是每个汇入该栅格的编号

汇流累积算法属于全局型算法，采用全局型算法并行策略，其核心计算环节如下：

(1)根据流向栅格判断每一个栅格的相邻栅格汇入其的个数n，如果n等于0则表明该栅格为计算源，则该点的汇流累积赋值为0，同时该点将入栈sourceStack。如果n大于0，则该点赋值为-n。

(2)将栈sourceStack中的栅格依次出栈，当出栈栅格mCell的下一个栅格nCell的值小于-1，则将其值加1同时等待下一次计算。如果nCell的值等于-1，则表明该栅格满足计算阀，将根据公式(4)计算的结果赋值给该栅格。以nCell为计算源，依次计算其下游栅格汇流累积值，直到有不满足计算条件的栅格。当栈sourceStack为空，则表明一次计算结束。

第四步、河网提取计算：河网是在汇流累积矩阵基础上形成的，它是通过所设定的阈值，将汇流累积值高于此阈值的格网连接起来，从而形成流域网络。具体算法为：遍历汇流累积数据，如果栅格的汇流累积值大于等于给定阈值，则将该栅格赋值为1，否则赋值为0。河网提取算法属于局部型算法，本发明对河网提取算法并行化采用局部型算法并行化策略。

第五步、坡长计算：坡长计算时坡度坡长因子计算的核心部分。根据USLE中的定义，坡长被定为从地表汇流产生的那点到出现明显渠道或者发生坡度减小导致发生沉积的那一点的水平距离。坡长计算的输入参数分别为：水流方向矩阵flowDirection,坡度矩阵slope,河网矩阵stream，坡度截断因子cutoff1和cutoff2。在计算过程中，全部栅格被定义为四类：(1)源点栅格，一般为局部高值点或位于数据块的边缘区域；(2)坡度截断栅格，即坡度骤降发生沉积区域，当栅格坡度小于2.87以及大于2.87度时，分别采用公式(5)和公式(6)进行坡段截断栅格判断。

$\mathrm{cutoffCell}(i,j)=\left\{\begin{array}{cc}1& \frac{\mathrm{slope}(\mathrm{prei},\mathrm{prej})-\mathrm{slope}(i,j)}{\mathrm{slope}(i,j)}>\mathrm{cutoffFactor}1\\ 0& \frac{\mathrm{slope}(\mathrm{prei},\mathrm{prej})-\mathrm{slope}(i,j)}{\mathrm{slope}(i,j)}\≤\mathrm{cutoffFactor}1\end{array}\right.---\left(5\right)$

$\mathrm{cutoffCell}(i,j)=\left\{\begin{array}{cc}1& \frac{\mathrm{slope}(\mathrm{prei},\mathrm{prej})-\mathrm{slope}(i,j)}{\mathrm{slope}(i,j)}>\mathrm{cutoffFactor}2\\ 0& \frac{\mathrm{slope}(\mathrm{prei},\mathrm{prej})-\mathrm{slope}(i,j)}{\mathrm{slope}(i,j)}\≤\mathrm{cutoffFactor}2\end{array}\right.---\left(6\right)$

式中，slope(i,j)和slope(prei,prej)分别是是栅格(i,j)和其上游栅格(prei,prej)处的坡度值，cutoffFactor1和cutoffFactor2分别是截断因子，由用户定义。

(3)河网截断栅格，即河网上的栅格；(4)普通栅格，即不属于以上三类的栅格。本发明中坡长计算的核心算法如下：

(1)栅格点标识：遍历流向矩阵、河网矩阵、坡度矩阵判断每一个栅格的类型。如果栅格mCell是河网截断栅格，则标识矩阵T(m)赋值为1并且坡长值赋值为0；如果mCell是源点栅格，则将其入栈sourceStack同时T(m)赋值为2；如果是坡度截断栅格，则T(m)赋值为3；如果是其他栅格，T(m)赋值为-n，n为mCell八领域栅格中流入它的栅格数。

(2)非累积坡长计算。非累积坡长即只考虑计算栅格自身的坡长值，根据栅格类型不同，具体计算公式如下：

式中NSL为非累积坡长值，cellsize为栅格分辨率大小

(3)累积坡长计算。对于栈sourceStack中的每个栅格依次出栈，如果出栈栅格iCell的流入栅格jCell的栅格标识T(jCell)不等于-1，则T(jCell)加1并等待下一次计算；如果T(jCell)等于-1,则表明jCell已满足计算条件，根据公式(8)计算其坡长。同时以jCell栅格为计算源顺水流方向依次判断其下游栅格，直至遇到不满足计算条件栅格。

${\mathrm{\λ}}_{i,j}=\underset{x=0,y=0}{\overset{x=i,y=j}{\mathrm{\Σ}}}\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_{x,y}---\left(8\right)$

式中k是坐标(x,y)处栅格的相邻栅格的编号，λ_x,y是每一个栅格的坡长值。

第六步、坡度坡长因子计算。坡度坡长因子计算是本流程的最后一步，其计算函数如公式(9)-(13)所示。具体计算时，对每一个栅格分别计算坡度坡长值，并将值写入结果文件。坡度坡长因子计算属于局部型算法，并行计算时采用本发明中的局部型算法并行化策略。

LS＝L*S   (9)

$L={\left(\frac{l}{22.13}\right)}^m---\left(10\right)$

$m=\frac{\mathrm{\β}}{1+\mathrm{\β}}---\left(11\right)$

β＝(sinθ/0.089)/[3*sinθ^0.8+0.56]   (12)

式中L是坡长因子，S是坡度因子，l是坡长值，m坡长指数变量，β是沟蚀与沟间侵蚀的比例，θ是坡度值。

以上步骤是本发明的具体的计算流程。此外，本发明采用的是从主式并行架构，主节点主要负责管理从节点、数据划分以及管理输入输出数据；从节点则主要负责输入数据读取、算法执行和输出数据写出。(如图4)

以下采用黄土高原数据为例，说明本发明。

黄土高原具有独特的地理景观，其黄土分布连续而深厚，水土流失严重为世界罕见，对于黄土高原土壤侵蚀的研究具有重要意义。现有研究中，由于单机计算资源受限，相关学者在对该区域土壤侵蚀建模时往往采用的较低分辨率DEM提取相关地形因子，这也造成了当前高精度DEM数据在相关研究中无法有效利用。本实例采用了覆盖黄土高原的Aster数据，数据分辨率为30米，数据量约6GB，数据已进行了洼地处理。具体实验基于IBM高性能计算集群，集群采用11台单机组成cluster结构，其中有一个主控节点(CPU：Xeon E5620主频：2.4GHz；内存：32G)，10台计算节点(CPU：Xeon E5650主频：2.4GHz；内存：16G)，各节点通过1000Mbps快速交换式以太网连接，系统环境为：Debain 6.0(64位)，应用软件为:OpenMPI 1.4.2、GDAL1.8.1。

坡度坡长因子并行计算具体流程如下：

第一步：坡度计算：调用SlopeClass完成坡度因子并行计算，其中输入数据为无洼地DEM数据。

第二步：流向计算：调用FlowDirection完成水流方向并行计算，其中输入数据为无洼地DEM数据。

第三步：汇流累积计算：调用FLowAccumulationClass完成汇流累积并行计算，其中输入数据为水流方向矩阵。

第四步：河网提取计算：调用StreamClass完成河网提取并行计算，其中输入数据为汇流累积矩阵。

第五步：坡长计算：调用SlopelengthClass完成坡长并行计算，其中输入数据为水流方向矩阵、坡度矩阵、河网矩阵，同时坡度截断因子分别设置为0.5和0.7。

第六步：坡度坡长因子计算：调用FactorClass完成坡度坡长因子并行计算，其中输入数据为坡度矩阵和坡长矩阵。

图6是黄土高原DEM数据，图7为坡度坡长因子计算结果。

Claims (4)

1.面向区域尺度土壤侵蚀建模的地形因子并行计算方法，其特征在于，包括以下过程：

(1)确定数据划分策略，其划分方式为：根据并行节点数目，采用行带缓冲区划分方法对原始地形数据进行等行划分，同时每个子数据块包含了三个要素：a)数据区：本计算节点所分配的数据区；b)缓冲区：每个子数据块附加了上下两个缓冲区，对于第一子数据块只有下缓冲区，对于最后一个数据块只有上缓冲区；c)控制点：用于定位每个数据块的起始读写位置；

(2)按照步骤(1)的数据划分策略，基于预处理后的DEM数据，进行坡度并行计算；

(3)按照步骤(1)的数据划分策略，基于预处理后的DEM数据，进行坡长并行计算；其中包含了水流方向计算、汇流累积计算、河网提取计算以及坡长计算四个子计算流程；

(4)按照步骤(1)的数据划分策略，基于坡度和坡长计算结果，进行坡度坡长因子并行计算。

2.根据权利要求1所述面向区域尺度土壤侵蚀建模的地形因子并行计算方法，其特征在于，步骤(2)具体的计算流程如下：

步骤21.根据数据划分策略，对输入的DEM数据进行数据划分，并创建结果数据文件；

步骤22.各子节点遍历输入的DEM数据，并执行坡度算法，计算完成后将计算结果写入结果文件。

3.根据权利要求1所述面向区域尺度土壤侵蚀建模的地形因子并行计算方法，其特征在于，步骤(3)具体的计算流程如下：

步骤31.根据数据划分策略，对输入的DEM数据进行数据划分，并创建结果数据文件；各子节点遍历输入的DEM数据，并执行水流方向算法，计算完成后将计算结果写入结果文件；

步骤32.根据数据划分策略，对水流方向矩阵进行数据划分，并创建汇流累积结果文件；各子节点遍历输入的水流方向矩阵，将源点栅格赋值为0同时入栈，非源点栅格赋值为-m，m为邻接汇入栅格数目，遍历完成后依次取出栈内栅格并沿流向计算其下游栅格汇流累积值，当栈内元素为空时计算结束；

步骤33.数据块间更新缓冲区信息，各子数据块已更新后的缓冲区栅格为源点栅格依次更新本数据区汇流累积值；一次更新计算完成后，遍历各子数据块边界行数据，若全部栅格完成计算则计算结束，否则需要重复步骤33；计算完成后各子节点将计算结果写入已创建的汇流累积文件；

步骤34.根据数据划分策略，对输入的汇流累积数据进行数据划分，并创建结果数据文件；各子节点遍历输入数据，并执行河网提取算法，计算完成后将计算结果写入结果文件；

步骤35.根据数据划分策略，分别对河网矩阵、坡度矩阵、水流方向矩阵进行数据划分，并创建坡长结果文件；各子节点遍历输入数据，根据侵蚀学坡长算法定义分别执行栅格点标识、非累积坡长计算、累积坡长计算步骤；

步骤36.数据块间更新缓冲区信息，各子数据块以更新后的缓冲区栅格为源点栅格依次更新本数据区的坡长值；一次更新计算完成后，遍历各子数据块边界行数据，若全部栅格完成计算则计算结束，否则需要重复步骤36；计算完成后各子节点将计算结果写入已创建的坡长文件。

4.根据权利要求1所述面向区域尺度土壤侵蚀建模的地形因子并行计算方法，其特征在于，步骤(4)具体的计算流程如下：

步骤41.根据数据划分策略，对输入的坡度、坡长数据进行数据划分，并创建结果数据文件；

步骤42.各子节点遍历输入数据，并执行坡度坡长因子算法，计算完成后将计算结果写入结果文件。