[发明专利]基于双层异步迭代策略的水动力学洪水演进模拟方法

权利要求书

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=1n·Rhy2/3·I···(4)]]>

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

I=h-h′d···(5)]]>

其中，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，得到洪水演进图。