CN103870699A - 基于双层异步迭代策略的水动力学洪水演进模拟方法 - Google Patents

Abstract

本发明提供了基于双层异步迭代策略的水动力学洪水演进模拟方法，用和数字高程模型同分辨率的栅格数据表示洪水数据，采用双层异步迭代算法模拟实际的水流过程；其中，外层迭代过程通过给定的迭代步长确定洪水演进时刻，内层迭代过程通过对比栅格流向相邻栅格所用的时间和栅格迭代步长的大小关系，确定栅格各自的内层迭代步长，从而实现迭代次数的自适应选择，迭代计算各栅格的断面水量，并通过设定输出洪水形态的时间间隔，动态绘制洪水演进图。采用FloodArea软件、均一化迭代算法以及本发明提供的方法分别对福建省万安流域暴雨洪水历史数据进行模拟，结果表明本发明模拟的结果与实际灾情的吻合度最高，能有效提升洪水演进模拟精度。

Description

基于双层异步迭代策略的水动力学洪水演进模拟方法

技术领域

本发明涉及基于双层异步迭代策略的水动力学洪水演进模拟方法，属于地理信息系统的虚拟地理环境研究领域。

背景技术

由于洪涝灾害频繁发生、破坏力强且难以预测，模拟和预报洪水演进过程成为了一项重要的研究任务。目前，模拟和预报洪水演进过程的方法主要有水文模型法、数值模拟法以及水动力学建模法等，前两类方法都需要以复杂的数学模型为基础，而水动力学建模算法则较为简单，但需要进行大量计算，而这一问题随着计算机性能的逐步提高得到了解决，因此，基于GIS栅格数据的水动力学洪水演进模型越来越成为学者们的研究热点。目前该方面的研究主要集中在两个方向：一是基于特定的地形地貌特征，对水动力学洪水演进模型做精简，使其能够更加便捷地运用到该类地形区域当中。如Batesand De Roo针对河道内洪水和河漫滩洪水演进的不同特征分别利用一维运动波和二维扩散波对模型进行简化；Hagen et al基于阿富汗的地形特征提出了一种简化的水动力学洪水演进模型，使其能够在灾害突发时更好地运用。二是完善模型，建立更为精准的洪水演进模型。例如李大鸣等为建立适应河道、滞洪区复杂情况的洪水演进一、二维衔接数学模型，采用有限体积法进行实现；Manoj and Vijay针对地表径流对圣维南方程组做扩散波近似法解算。

上述研究内容的实现，都需要与具体的迭代策略相匹配。目前，较为通用的方法为均一化迭代法。均一化迭代是指在一个迭代过程中赋予所有栅格相同的迭代步长，其具体迭代过程如图7所示：首先利用数字高程模型（DigitalElevation Model）获得水文特征参数信息，如水面比降、水流方向等；然后计算所有栅格的断面流量q、栅格水深WL等信息；当一次迭代计算完毕，在总模拟时长T中减去Δt并进行下一次迭代计算；当T小于等于0时，迭代计算结束，此时的栅格数据即为淹没模拟时长T的洪水演进图，栅格值即淹没水深。

由于数字高程模型（Digital Elevation Model）各个栅格单元的最大水面比降不同，导致栅格中水的流速和流向邻近栅格所需时间也不同。而均一化迭代过程赋予了所有栅格相同的迭代步长，这就造成了迭代步长和实际所需时间的不一致，影响了栅格单元水位、流量的计算精度。通过计算某时刻万安流域的数字高程模型（Digital Elevation Model）中每个栅格单元分辨率和流速的比值，统计流向临近栅格单元的所需时间，统计结果如图8所示。根据图8可知，流向邻近栅格单元所需的时间主要集中在5s—129s之间，在129s—253s之间的栅格也仍有分布。因此，在大时间跨度的区间内，均一化迭代法在一次迭代过程中用相同的迭代步长来计算栅格单元间的断面流量，势必会影响洪水演进模拟的精度。在复杂的地形条件下，迭代步长均一化带来的洪水演进误差尤为明显。

发明内容

为了解决现有技术的不足，本发明提供了基于双层异步迭代策略的水动力学洪水演进模拟方法，通过分析迭代过程中不同栅格单元洪水的水动力学特征，采用地理信息系统（Geographic Information System）栅格分析技术，在数字高程模型（Digital Elevation Model）地形因子分析的基础上，提出双层异步迭代算法，其中通过内层迭代过程实现迭代次数的自适应选择，迭代计算各格栅的断面水量，从而更为精确地动态模拟洪水的演进过程。

本发明为解决其技术问题所采用的技术方案是：

本发明提供了基于双层异步迭代策略的水动力学洪水演进模拟方法，用和数字高程模型同分辨率的栅格数据表示洪水数据，采用双层异步迭代算法模拟实际的水流过程；其中，外层迭代过程通过给定的迭代步长确定洪水演进时刻，内层迭代过程通过对比栅格流向相邻栅格所用的时间和栅格迭代步长的大小关系，确定栅格各自的内层迭代步长，从而实现迭代次数的自适应选择，迭代计算各栅格的断面水量，并通过设定输出洪水形态的时间间隔，动态绘制洪水演进图。

所述的基于双层异步迭代策略的水动力学洪水演进模拟方法具体包括以下步骤：

（1）洪水数据用和数字高程模型同分辨率的栅格数据表示，包括栅格总数N、各栅格的栅格水深WL、各栅格的栅格高程值Elevation、各栅格与相邻栅格的栅格间距d；设置一次完整模拟过程的模拟时长T；内层迭代次数指示BoolFirst初始化为true；变量ΔpauseT初始化为0；

（2）进入外层迭代过程，通过以下公式计算进行一次外层迭代的整幅迭代步长Δt：

Δt=∑Tflow_k·P_k……………………………………….….（1）

其中，Tflow_k表示当前栅格流向相邻栅格的时间，P_k表示N个栅格中出现时间为Tflow_k值的频率，k表示出现Tflow_k的栅格个数，P_k通过以下公式计算：

P_k=k/N……………………………………...…………….….（2）

公式（1）中，Tflow_k∈(μ-2σ,μ+2σ)，μ表示Tflow_k的期望，σ表示方差；每个栅格各自的Tflow通过以下公式计算：

Tflow=d/V…………….…………...…………..（3）

其中，V表示当前栅格的栅格流速；栅格流速V通过曼宁公式计算：

$V=\frac{1}{n}\·{R_{\mathrm{hy}}}^{2/3}\·\sqrt{I}\·\·\·\left(4\right)$

其中，R_hy为水力半径，水力半径为栅格水深WL，n为曼宁系数，I为水面比降；水面比降I通过以下公式计算：

$I=\frac{h-h^{\′}}{d}\·\·\·\left(5\right)$

其中，h表示当前栅格的自由液面高度，h'表示相邻栅格的自由液面高度，d表示栅格间距；自由液面高度H通过以下公式计算：

H=Elevation+WL….…………………………..（6）

（3）设当前栅格为主栅格，进入该主栅格的内层迭代过程；初始化变量，AdjectNum赋值为0，AdjectNum2赋值为0，bstop赋值为false；

（4）若bstop为true，则迭代结束，进入步骤（11），否则进入步骤（5）；

（5）更新AdjectNum2，使AdjectNum2赋值为AdjectNum，初始化变量j为0；采用D8算法，通过公式（5）分别计算主栅格与周围相邻8个栅格单元的水面比降，取水面比降最大的相邻栅格为主栅格的流向栅格；判断BoolFirst的取值：

（a）若BoolFirst为true，则表示本次计算是主栅格第一次进入内层迭代，并通过以下公式计算主栅格进行本次内层迭代的栅格迭代步长ΔcellT：

ΔcellT=Δt+ΔpauseT.…………………………(7）

其中，Δt为步骤（2）所述的整幅迭代步长，ΔpauseT为主栅格进行上一次内层迭代计算后栅格余下的步长；更新BoolFirst为false，通过公式（3）计算主栅格到其流向栅格所用时间Tflow；

（b）若BoolFirst为false，则更新ΔcellT，将ΔcellT减去上次迭代计算时主栅格到其流向栅格的时间Tflow，Tflow通过公式（3）计算得到；

（6）比较Tflow与栅格迭代步长ΔcellT：若ΔcellT＜0.5Tflow，则主栅格水量滞留，更新ΔpauseT，使ΔcellT累加到ΔpauseT，bstop更新为true，进入步骤（8）；若0.5Tflow≤ΔcellT≤Tflow，则计算ΔcellT时间内主栅格到其流向栅格的水量，更新主栅格及其流向的栅格水深WL和栅格高程值Elevation，ΔpauseT更新为0，bstop更新为true，进入步骤（8）；若ΔcellT＞Tflow，则计算Tflow时间内主栅格到其流向栅格的水量，更新主栅格及其流向的栅格水深WL和栅格高程值Elevation，ΔpauseT更新为0，此时主栅格的该流向栅格成为主栅格的附属栅格j，进入步骤（7）；其中，主栅格到其流向栅格的水量通过以下公式计算：

Q=V·t·A…………………………….………..（8）

其中，Q为流向相邻栅格的水量，V为栅格流速，t为时间间隔，t根据Tflow与ΔcellT比较的结果赋值为ΔcellT或Tflow，A为栅格截面面积；其中，栅格流速V通过公式（4）计算，栅格截面面积A通过以下公式计算：

A=d·h_flow…….………………………………….（9）

其中，d为栅格间距，h_flow表示两个栅格间的断面水深；两个栅格间的断面水深h_flow用主栅格的栅格水深减去主栅格的栅格高程值和附属栅格的栅格高程值之间的较大值计算得到；

（7）若附属栅格j已被标记过是主栅格的附属栅格，则进入步骤（8）；否则，更新AdjectNum，使AdjectNum累加1，进入步骤（8）；

（8）若j≥AdjectNum2，则表示当前主栅格的所有附属栅格水量已计算过，进入步骤（4）以判断主栅格的本次内层迭代过程是否结束；否则进入步骤（9）；

（9）采用D8算法，通过公式（5）分别计算附属栅格j与其周围相邻8个栅格单元的水面比降，取水面比降最大的相邻栅格为附属栅格j的流向栅格；通过以下公式计算附属栅格j的栅格迭代步长ΔcellTj：

ΔcellT_j=Tflow+ΔpauseT….………………….（11）

其中，ΔpauseT根据步骤（5）所述的Tflow与ΔcellT比较的结果赋值，Tflow为主栅格流向附属栅格j所用的时间；

（10）通过公式（3）计算得到附属栅格j到其流向栅格所用的时间Tflow_j，比较Tflow_j与附属栅格j的栅格迭代步长ΔcellT_j：若ΔcellT_j＜0.5Tflow_j，则附属栅格j水量滞留，更新ΔpauseT使ΔcellT_j累加到该附属栅格j的滞留迭代步长ΔpauseT，更新j，使j累加1，返回步骤（8）；若0.5Tflow_j≤ΔcellT_j≤Tflow_j，则通过公式（8）计算ΔcellT_j时间内附属栅格j到其流向栅格的水量，更新附属栅格j及其流向的栅格水深WL和栅格高程值Elevation，更新ΔpauseT为0，更新j，使j累加1，返回步骤（8）；若ΔcellT_j＞Tflow_j，则附属栅格j的流向栅格成为主栅格的附属栅格，通过公式（8）计算Tflow_j时间内附属栅格j到其流向栅格的水量，更新附属栅格j及其流向的栅格水深WL和栅格高程值Elevation，更新ΔpauseT为0，更新j，使j累加1，返回步骤（7）；

（11）返回步骤（2）以进入下一主栅格的内层迭代过程，直到所有N个栅格都完成内层迭代过程，进入步骤（12）；

（12）更新T，将T减去Δt，则完成了一次外层迭代过程；若T≤0，进入（13）；否则返回（2），以进行下一次外层迭代过程；

（13）判断是否查看当前洪水形态，若是，则输出当前时刻各栅格的水量Q，返回步骤（2）以进行下一次外层迭代过程；否则直接进入步骤（2）以进行下一次外层迭代过程；

（14）输出T时间段内各栅格的水量Q，得到洪水演进图。

本发明相比于现有技术具有的有益效果是：

（1）本发明提供的方法解决了基于均一化迭代的洪水演进模拟中迭代步长和单个栅格水流演进时间不一致的问题，通过外层迭代控制洪水演进的时间尺度，同时，使用内层迭代过程通过对比栅格流向相邻栅格所用的时间和栅格迭代步长的大小关系，确定栅格各自的内层迭代步长，从而实现迭代次数的自适应选择，迭代计算各栅格的断面水量；

（2）同传统的均一化迭代法和FloodArea软件模拟的结果相比，本发明提供的方法模拟结果与灾情点实际水深平均误差更小，证明本发明提供的方法模拟精度更高。

附图说明

图1是本发明的双层异步迭代过程；

图2是万安流域调查受灾点与模拟淹没水深比较图；

图3是降雨2小时洪水淹没形态模拟分布对比；

图4是降雨6小时洪水淹没形态模拟分布对比；

图5是降雨10小时洪水淹没形态模拟分布对比；

图6是降雨13小时洪水淹没形态模拟分布对比；

图7是均一化迭代算过程；

图8是流向邻近栅格单元时间统计图。

具体实施方式

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

本发明用和数字高程模型同分辨率的栅格数据表示洪水数据，采用双层异步迭代算法模拟实际的水流过程；其中，外层迭代过程通过给定的迭代步长确定洪水演进时刻，内层迭代过程通过对比栅格流向相邻栅格所用的时间和栅格迭代步长的大小关系，确定栅格各自的内层迭代步长，从而实现迭代次数的自适应选择，迭代计算各栅格的断面水量，并通过设定输出洪水形态的时间间隔，动态绘制洪水演进图。

参照图1，本发明提供的基于双层异步迭代策略的水动力学洪水演进模拟方法的步骤如下：

（1）洪水数据用和数字高程模型同分辨率的栅格数据表示，包括栅格总数N、各栅格的栅格水深WL、各栅格的栅格高程值Elevation、各栅格与相邻栅格的栅格间距d；设置一次完整模拟过程的模拟时长T；内层迭代次数指示BoolFirst初始化为true；变量ΔpauseT初始化为0；

（2）进入外层迭代过程，通过以下公式计算进行一次外层迭代的整幅迭代步长Δt：

Δt=∑Tflow_k·P_k……………………………………….….（1）

其中，Tflow_k表示当前栅格流向相邻栅格的时间，P_k表示N个栅格中出现时间为Tflow_k值的频率，k表示出现Tflow_k的栅格个数，P_k通过以下公式计算：

P_k=k/N……………………………………...…………….….（2）

公式（1）中，Tflow_k∈(μ-2σ,μ+2σ)，μ表示Tflowk的期望，σ表示方差；每个栅格各自的Tflow通过以下公式计算：

Tflow=d/V…………….…………...…………..（3）

其中，V表示当前栅格的栅格流速；栅格流速V通过曼宁公式计算：

$V=\frac{1}{n}\·{R_{\mathrm{hy}}}^{2/3}\·\sqrt{I}\·\·\·\left(4\right)$

其中，R_hy为水力半径，水力半径为栅格水深WL，n为曼宁系数，I为水面比降；水面比降I通过以下公式计算：

$I=\frac{h-h^{\′}}{d}\·\·\·\left(5\right)$

其中，h表示当前栅格的自由液面高度，h'表示相邻栅格的自由液面高度，d表示栅格间距；自由液面高度H通过以下公式计算：

H=Elevation+WL….…………………………..（6）

（3）设当前栅格为主栅格，进入该主栅格的内层迭代过程；初始化变量，AdjectNum赋值为0，AdjectNum2赋值为0，bstop赋值为false；

（4）若bstop为true，则迭代结束，进入步骤（11），否则进入步骤（5）；

（5）更新AdjectNum2，使AdjectNum2赋值为AdjectNum，初始化变量j为0；采用D8算法，通过公式（5）分别计算主栅格与周围相邻8个栅格单元的水面比降，取水面比降最大的相邻栅格为主栅格的流向栅格；判断BoolFirst的取值：

（a）若BoolFirst为true，则表示本次计算是主栅格第一次进入内层迭代，并通过以下公式计算主栅格进行本次内层迭代的栅格迭代步长ΔcellT：

ΔcellT=Δt+ΔpauseT.…………………………(7）

其中，Δt为步骤（2）所述的整幅迭代步长，ΔpauseT为主栅格进行上一次内层迭代计算后栅格余下的步长；更新BoolFirst为false，通过公式（3）计算主栅格到其流向栅格所用时间Tflow；

（b）若BoolFirst为false，则更新ΔcellT，将ΔcellT减去上次迭代计算时主栅格到其流向栅格的时间Tflow，Tflow通过公式（3）计算得到；

（6）比较Tflow与栅格迭代步长ΔcellT：若ΔcellT＜0.5Tflow，则主栅格水量滞留，更新ΔpauseT，使ΔcellT累加到ΔpauseT，bstop更新为true，进入步骤（8）；若0.5Tflow≤ΔcellT≤Tflow，则计算ΔcellT时间内主栅格到其流向栅格的水量，更新主栅格及其流向的栅格水深WL和栅格高程值Elevation，ΔpauseT更新为0，bstop更新为true，进入步骤（8）；若ΔcellT＞Tflow，则计算Tflow时间内主栅格到其流向栅格的水量，更新主栅格及其流向的栅格水深WL和栅格高程值Elevation，ΔpauseT更新为0，此时主栅格的该流向栅格成为主栅格的附属栅格j，进入步骤（7）；其中，主栅格到其流向栅格的水量通过以下公式计算：

Q=V·t·A…………………………….………..（8）

其中，Q为流向相邻栅格的水量，V为栅格流速，t为时间间隔，t根据Tflow与ΔcellT比较的结果赋值为ΔcellT或Tflow，A为栅格截面面积；其中，栅格流速V通过公式（4）计算，栅格截面面积A通过以下公式计算：

A=d·h_flow…….………………………………….（9）

其中，d为栅格间距，h_flow表示两个栅格间的断面水深；两个栅格间的断面水深h_flow用主栅格的栅格水深减去主栅格的栅格高程值和附属栅格的栅格高程值之间的较大值计算得到；

（7）若附属栅格j已被标记过是主栅格的附属栅格，则进入步骤（8）；否则，更新AdjectNum，使AdjectNum累加1，进入步骤（8）；

（8）若j≥AdjectNum2，则表示当前主栅格的所有附属栅格水量已计算过，进入步骤（4）以判断主栅格的本次内层迭代过程是否结束；否则进入步骤（9）；

（9）采用D8算法，通过公式（5）分别计算附属栅格j与其周围相邻8个栅格单元的水面比降，取水面比降最大的相邻栅格为附属栅格j的流向栅格；通过以下公式计算附属栅格j的栅格迭代步长ΔcellT_j：

ΔcellT_j=Tflow+ΔpauseT….………………….（11）

其中，ΔpauseT根据步骤（5）所述的Tflow与ΔcellT比较的结果赋值，Tflow为主栅格流向附属栅格j所用的时间；

（10）通过公式（3）计算得到附属栅格j到其流向栅格所用的时间Tflow_j，比较Tflow_j与附属栅格j的栅格迭代步长ΔcellT_j：若ΔcellT_j＜0.5Tflow_j，则附属栅格j水量滞留，更新ΔpauseT使ΔcellT_j累加到该附属栅格j的滞留迭代步长ΔpauseT，更新j，使j累加1，返回步骤（8）；若0.5Tflow_j≤ΔcellT_j≤Tflow_j，则通过公式（8）计算ΔcellT_j时间内附属栅格j到其流向栅格的水量，更新附属栅格j及其流向的栅格水深WL和栅格高程值Elevation，更新ΔpauseT为0，更新j，使j累加1，返回步骤（8）；若ΔcellT_j＞Tflow_j，则附属栅格j的流向栅格成为主栅格的附属栅格，通过公式（8）计算Tflow_j时间内附属栅格j到其流向栅格的水量，更新附属栅格j及其流向的栅格水深WL和栅格高程值Elevation，更新ΔpauseT为0，更新j，使j累加1，返回步骤（7）；

（11）返回步骤（2）以进入下一主栅格的内层迭代过程，直到所有N个栅格都完成内层迭代过程，进入步骤（12）；

（12）更新T，将T减去Δt，则完成了一次外层迭代过程；若T≤0，进入（14）；否则进入步骤（13）；

（13）判断是否查看当前洪水形态，若是，则输出当前时刻各栅格的水量Q，返回步骤（2）以进行下一次外层迭代过程；否则直接进入步骤（2）以进行下一次外层迭代过程；

（14）输出T时间段内各栅格的水量Q，得到洪水演进图。

分别利用均一化迭代法、FloodArea软件以及本发明提供的基于双层异步迭代策略的水动力学洪水演进模拟方法模拟万安流域的暴雨洪涝历史过程，并针对7个典型受灾点进行实验对比分析，并将结果和实际检测结果进行对比，对比结果如参照图2所示。整理数据后通过以下公式计算模拟结果和实际检测结果的平均误差：

$\mathrm{Error}=\underset{i=0}{\overset{n}{\mathrm{\Σ}}}|h_i-{\mathrm{hreal}}_i|/n$

其中，h_i表示第i个灾情考察点的模拟水深，hreal_i表示第i个灾情考察点的实际水深，n表示灾情考察点的个数。

三种模拟方法的平均误差如下表所示，单位为m：

	均一化迭代法	FloodArea	异步迭代法
				平均误差	1.014	0.721	0.360

可知本发明提供的基于双层异步迭代策略的水动力学洪水演进模拟方法与灾情点实际水深平均误差最小，该方法模拟精度最高。但是同时也发现，三种方法的模拟误差均较大，这是由于1:50000的DEM对中小流域地形表现力有限，未考虑蒸发渗透等水文模型，以及未考虑人工建筑对洪水走向改变的影响等多方面原因造成的。

在洪水淹没形态分布对比方面，分别选取降雨2小时、6小时、10小时、13小时这四个时次的模拟结果对比，如图3至图6所示。其中，a表示均一化迭代法的模拟结果，b表示FloodArea的模拟结果，c表示本发明提供的基于双层异步迭代策略的水动力学洪水演进模拟方法的模拟结果。

观察各时段模拟对比情况，其中降雨2小时三个模拟结果的淹没状况较为一致；降雨6小时后，本发明提供的基于双层异步迭代策略的水动力学洪水演进模拟方法0.5-1米的淹没范围较其他两个略大；在降雨10小时和13小时后，本发明提供的基于双层异步迭代策略的水动力学洪水演进模拟方法和FloodArea的淹没状况较为一致，而均一化迭代模拟出现了较小的“河网分叉”。从总体分布上看，三种模拟方法的结果在淹没形态分布上基本一致。

通过淹没水深和洪水形态分布的对比分析可知，本发明提供的基于双层异步迭代策略的水动力学洪水演进模拟方法的模拟结果在淹没形态上和其他两种方法基本一致，并且在与实际灾情点淹没水深匹配上是最好的。

Claims (2)

1.基于双层异步迭代策略的水动力学洪水演进模拟方法，其特征在于：用和数字高程模型同分辨率的栅格数据表示洪水数据，采用双层异步迭代算法模拟实际的水流过程；其中，外层迭代过程通过给定的迭代步长确定洪水演进时刻，内层迭代过程通过对比栅格流向相邻栅格所用的时间和栅格迭代步长的大小关系，确定栅格各自的内层迭代步长，从而实现迭代次数的自适应选择，迭代计算各栅格的断面水量，并通过设定输出洪水形态的时间间隔，动态绘制洪水演进图。

2.根据权利要求1所述的基于双层异步迭代策略的水动力学洪水演进模拟方法具体包括以下步骤：

（1）洪水数据用和数字高程模型同分辨率的栅格数据表示，包括栅格总数N、各栅格的栅格水深WL、各栅格的栅格高程值Elevation、各栅格与相邻栅格的栅格间距d；设置一次完整模拟过程的模拟时长T；内层迭代次数指示BoolFirst初始化为true；变量ΔpauseT初始化为0；

（2）进入外层迭代过程，通过以下公式计算进行一次外层迭代的整幅迭代步长Δt：

Δt=∑Tflow_k·P_k………………………………………….….（1）

其中，Tflow_k表示当前栅格流向相邻栅格的时间，Pk表示N个栅格中出现时间为Tflow_k值的频率，k表示出现Tflow_k的栅格个数，P_k通过以下公式计算：

P_k=k/N………………………………………...…………….….（2）

公式（1）中，Tflow_k∈(μ-2σ,μ+2σ)，μ表示Tflow_k的期望，σ表示方差；每个栅格各自的Tflow通过以下公式计算：

Tflow=d/V………………………………………...……………（3）

其中，V表示当前栅格的栅格流速；栅格流速V通过曼宁公式计算：

$V=\frac{1}{n}\·{R_{\mathrm{hy}}}^{2/3}\·\sqrt{I}\·\·\·\left(4\right)$

其中，R_hy为水力半径，水力半径为栅格水深WL，n为曼宁系数，I为水面比降；水面比降I通过以下公式计算：

$I=\frac{h-h^{\′}}{d}\·\·\·\left(5\right)$

其中，h表示当前栅格的自由液面高度，h'表示相邻栅格的自由液面高度，d表示栅格间距；自由液面高度H通过以下公式计算：

H=Elevation+WL….…………..………………..（6）

（3）设当前栅格为主栅格，进入该主栅格的内层迭代过程；初始化变量，AdjectNum赋值为0，AdjectNum2赋值为0，bstop赋值为false；

（4）若bstop为true，则迭代结束，进入步骤（11），否则进入步骤（5）；

（5）更新AdjectNum2，使AdjectNum2赋值为AdjectNum，初始化变量j为0；采用D8算法，通过公式（5）分别计算主栅格与周围相邻8个栅格单元的水面比降，取水面比降最大的相邻栅格为主栅格的流向栅格；判断BoolFirst的取值：

（a）若BoolFirst为true，则表示本次计算是主栅格第一次进入内层迭代，并通过以下公式计算主栅格进行本次内层迭代的栅格迭代步长ΔcellT：

ΔcellT=Δt+ΔpauseT.……………………………(7）

其中，Δt为步骤（2）所述的整幅迭代步长，ΔpauseT为主栅格进行上一次内层迭代计算后栅格余下的步长；更新BoolFirst为false，通过公式（3）计算主栅格到其流向栅格所用时间Tflow；

（b）若BoolFirst为false，则更新ΔcellT，将ΔcellT减去上次迭代计算时主栅格到其流向栅格的时间Tflow，Tflow通过公式（3）计算得到；

（6）比较Tflow与栅格迭代步长ΔcellT：若ΔcellT＜0.5Tflow，则主栅格水量滞留，更新ΔpauseT，使ΔcellT累加到ΔpauseT，bstop更新为true，进入步骤（8）；若0.5Tflow≤ΔcellT≤Tflow，则计算ΔcellT时间内主栅格到其流向栅格的水量，更新主栅格及其流向的栅格水深WL和栅格高程值Elevation，ΔpauseT更新为0，bstop更新为true，进入步骤（8）；若ΔcellT＞Tflow，则计算Tflow时间内主栅格到其流向栅格的水量，更新主栅格及其流向的栅格水深WL和栅格高程值Elevation，ΔpauseT更新为0，此时主栅格的该流向栅格成为主栅格的附属栅格j，进入步骤（7）；其中，主栅格到其流向栅格的水量通过以下公式计算：

Q=V·t·A…….………………………….………..（8）

其中，Q为流向相邻栅格的水量，V为栅格流速，t为时间间隔，t根据Tflow与ΔcellT比较的结果赋值为ΔcellT或Tflow，A为栅格截面面积；其中，栅格流速V通过公式（4）计算，栅格截面面积A通过以下公式计算：

A=d·h_flow…….………………………………….（9）

其中，d为栅格间距，h_flow表示两个栅格间的断面水深；两个栅格间的断面水深h_flow用主栅格的栅格水深减去主栅格的栅格高程值和附属栅格的栅格高程值之间的较大值计算得到；

（7）若附属栅格j已被标记过是主栅格的附属栅格，则进入步骤（8）；否则，更新AdjectNum，使AdjectNum累加1，进入步骤（8）；

（8）若j≥AdjectNum2，则表示当前主栅格的所有附属栅格水量已计算过，进入步骤（4）以判断主栅格的本次内层迭代过程是否结束；否则进入步骤（9）；

（9）采用D8算法，通过公式（5）分别计算附属栅格j与其周围相邻8个栅格单元的水面比降，取水面比降最大的相邻栅格为附属栅格j的流向栅格；通过以下公式计算附属栅格j的栅格迭代步长ΔcellT_j：

ΔcellT_j=Tflow+ΔpauseT…….……….………..（11）

其中，ΔpauseT根据步骤（5）所述的Tflow与ΔcellT比较的结果赋值，Tflow为主栅格流向附属栅格j所用的时间，Tflow已在步骤（6）求得；

（10）通过公式（3）计算得到附属栅格j到其流向栅格所用的时间Tflow_j，比较Tflow_j与附属栅格j的栅格迭代步长ΔcellT_j：若ΔcellT_j＜0.5Tflow_j，则附属栅格j水量滞留，更新ΔpauseT使ΔcellT_j累加到该附属栅格j的滞留迭代步长ΔpauseT，更新j，使j累加1，返回步骤（8）；若0.5Tflow_j≤ΔcellT_j≤Tflow_j，则通过公式（8）计算ΔcellT_j时间内附属栅格j到其流向栅格的水量，更新附属栅格j及其流向的栅格水深WL和栅格高程值Elevation，更新ΔpauseT为0，更新j，使j累加1，返回步骤（8）；若ΔcellT_j＞Tflow_j，则附属栅格j的流向栅格成为主栅格的附属栅格，通过公式（8）计算Tflow_j时间内附属栅格j到其流向栅格的水量，更新附属栅格j及其流向的栅格水深WL和栅格高程值Elevation，更新ΔpauseT为0，更新j，使j累加1，返回步骤（7）；

（11）返回步骤（2）以进入下一主栅格的内层迭代过程，直到所有N个栅格都完成内层迭代过程，进入步骤（12）；

（12）更新T，将T减去Δt，则完成了一次外层迭代过程；若T≤0，进入（14）；否则进入步骤（13）；

（13）判断是否查看当前洪水形态，若是，则输出当前时刻各栅格的水量Q，返回步骤（2）以进行下一次外层迭代过程；否则直接进入步骤（2）以进行下一次外层迭代过程；

（14）输出T时间段内各栅格的水量Q，得到洪水演进图。

Images