CN102508961A - 一种高分辨率的全分布式水文模型topx的设计方法 - Google Patents

Abstract

本发明公开了一种高分辨率的全分布式水文模型TOPX的设计方法，属于全分布式水文模型的设计方法领域。本发明结合新的地形指数计算方法，建立具有地形指数尺度转换机制高分辨率的全分布式水文模型TOPX方法，其步骤为：A）高分辨率的全分布式水文模型TOPX模型数据准备；B)高分辨率的全分布式水文模型TOPX地形指数的计算；C）高分辨率的全分布式水文模型TOPX模型产流计算；D）高分辨率的全分布式水文模型TOPX模型土壤湿度计算；E）高分辨率的全分布式水文模型TOPX模型的汇流计算；F）高分辨率的全分布式水文模型TOPX模型测试结果及其分析。本发明能够简单高效的处理DEM数据，提高地形指数的效率，有效的提高了TOPX模型在在高分辨率的DEM区域的应用。

Description

一种高分辨率的全分布式水文模型TOPX的设计方法

技术领域

本发明涉及一种结合新的地形指数的全分布式水文模型设计方法，具体的说是一种借助于即考虑地形指数概念和新安江模型的水量平衡思想的TOPX模型，结合新的地形指数计算方法，建立耦合新的地形指数计算方法的具有地形指数尺度转换机制的高分辨率的全分布式水文模型TOPX的设计方法。

背景技术

水文过程是一种非常复杂的大气、陆面及地下水循环现象，流域水文模型是水文科学研究的有效工具。所谓流域水文模型是用数学的方法描述和模拟水文循环的过程，即将流域概化成一个系统，降水量是系统的输入，流量是系统的输出，同样，蒸发和壤中流也可以被认为是输出。水文模型中典型代表有TOPMODEL和新安江模型等。

TOPMODEL模型和新安江模型(徐精文.基于蓄水容量曲线和TOPMODEL概念的流域水文模型构建及其与Noah LSM的耦合研究[D].北京：中国科学院大气物理所，2009.)都是基于变动产流源面积理论的分布式概念性模型，这两个模型结构简单，优选参数少，尤其是TOPMODEL模型的优选参数更少，而且与观测的物理水文过程有密切联系，模拟效果好，已被应用到水科学研究的各个方面。然而，新安江模型只考虑了土壤中含水量的变化，无法有效表达地下水位的动态变化，而地下水位对大气模式和气候模式中水分和能量平衡有重要的影响。TOPMODEL模型对地下水位的刻画上存在一个严重问题，即TOPMODEL模型刻画地下水动态变化时存在与实际情况偏差的问题(平均地下水水位线高于地面)，而且仅仅通过改进该模式很难解决这一问题。具有地形指数尺度转换机制的大尺度陆面水文过程模型TOPX(雍斌.陆面水文过程模型TOPX构建及其与区域气候模式RIEMS的耦合应用[D].南京：南京大学，2008.)很好的解决这些问题，该模型能够简单高效的表达陆地表面地形地貌特征的，同时又能较好描述二维水文过程，提高了水文模型的模拟效果，对大尺度陆面水文过程模拟有良好的稳定性。

在TOPX模型中，由于涉及大尺度，DEM网格大，地形指数的计算量小，应用方便。在对于高分辨率DEM的研究区域，面临以下问题：

涉及高分辨率的研究流域，DEM网格分辨率小，DEM网格数量多，DEM数据量大时，地形指数的计算工作量则会非常巨大，处理大的DEM数据时非常耗时，导致死机，严重影响TOPX模型在高分辨率的DEM区域的应用。

发明内容

1.发明要解决的技术问题

本发明的目的是借助于考虑地形指数概念和新安江模型的水量平衡思想的TOPX模型，结合新的地形指数计算方法，建立构建具有耦合新的地形指数方法的地形指数尺度转换机制高分辨率的全分布式水文模型TOPX(见图1)，通过利用当地的数字高程模型数据、土壤类型矢量数据、土地利用类型栅格数据和水文气象数据，对高分辨率的全分布式水文模型TOPX模型参数进行率定，以便高分辨率的全分布式水文模型TOPX模型能准确的模拟当地的水文过程。该模型能够简单高效的处理DEM数据，提高地形指数的效率，解决TOPX模型在高分辨率研究流域情况下，DEM网格数量多，数据量大，地形指数的计算工作量非常巨大，处理大的DEM数据非常耗时，甚至导致死机的问题，有效的提高了TOPX模型在在高分辨率的DEM区域的应用。

2.技术方案

本发明目的通过以下方案实现：

本发明的一种高分辨率的全分布式水文模型TOPX的设计方法，其步骤为：

A)高分辨率的全分布式水文模型TOPX模型数据准备；

B)高分辨率的全分布式水文模型TOPX地形指数的计算；

C)高分辨率的全分布式水文模型TOPX模型产流计算；

D)高分辨率的全分布式水文模型TOPX土壤湿度计算；

E)高分辨率的全分布式水文模型TOPX模型的汇流计算；

F)模型测试结果及其分析。

步骤A：高分辨率的全分布式水文模型TOPX模型数据准备包括：

1)高分辨率的全分布式水文模型TOPX模型需要的气象数据包括：逐日的最高气温、最低气温、降水量、相对湿度、风速、日照时数和20cm小型蒸发皿蒸发量；

2)高分辨率的全分布式水文模型TOPX模型需要的水文数据：逐日径流量数据和次洪数据；

3)高分辨率的全分布式水文模型TOPX模型需要的基础地理数据：数字高程模型数据、土壤类型矢量数据和土地利用类型栅格数据。

步骤B高分辨率的全分布式水文模型TOPX模型地形指数的计算：

1)填洼算法：

填洼算法先建立两个栈(栈A和栈B)，把临时DEM中分位数最小的那些栅格的位置信息(i行j列)存储在栈A或栈B中，分位数最小的那一类中每个栅格(用c表示)进行循环处理，如果栅格单元属于洼地则令其等于其相邻栅格中最低的一个；若不是洼地则令其等于它对应的初始矩阵的值，并且将此栅格从栈中剔出，将没有被剔出的栅格存放到另一个栈中(栈B或栈A)。然后再对栈B或栈A进行类似的处理。重复此过程，直到栈中只剩下洼地为止，这些洼地的位置信息最终会保存在两个栈中的一个里。第一类洼地填充和平坦区处理的C++伪码如表3-1所示。

表3-1中符号δ(n)要么为0，要么代表微小高程。如果δ(n)取值为0，洼地被填成平地；如果δ(n)取一个微小高程，比如对于对角线方向取0.000141，非对角线方向取0.0001，则初始DEM被新方法处理后，没有平坦区的存在。

表3-1对分位数最小的一类进行填洼的算法伪码

第一类处理完后，在第二类处理之前还有一关键步骤：检查栅格中是否还存在数值等于该类分位数的栅格，如果存在，则将其赋值为下一类(即第二类)分位数，因为这些栅格虽然是洼地，但其在空间上与填洼后的第一类不相邻，需要进一步处理。之后进行第二类，第三类等填洼处理，直到所有类全部处理完，这时的CDEM就是经过填洼和平坦区处理后的DEM，此时所有的洼地对应的位置信息存放在栈A或栈B两个中的一个之中，这非常有利用随后的水文模拟与分析。

2)地形指数提取方法：

地形指数提取方法首先用分位数法将DEM分成8类，然后将DEM处理成无洼地及平坦区的DEM。其流程如图2所示。然后，由于DEM中无洼地及平坦区，可直接根据D8算法指定水流方向(JENSON S K，DOMINGUE J O.Extractingtopographic structure from digital elevation data for geographic information systemanalysis[J].Photogrammetric Engineering and Remote Sensing，1988，54(11)：1593-1600.)。最后进行地形指数的计算，其具体计算方法参照Wolock and McCabe(2000)的方法(Wolock D M，McCabe G J.Differences in topographic characteristicscomputed from 100-and 1000-m resolution digital elevation model data[J].Hydrological Processes，2000，14(6)：987-1002.)。

步骤C：高分辨率的全分布式水文模型TOPX模型产流计算：

借鉴TOPX模型的产流计算方法(雍斌.陆面水文过程模型TOPX构建及其与区域气候模式RIEMS的耦合应用[D].南京：南京大学，2008.)，建立模型的产流方案。

1)在BATS模型基础上(Bonan G B，Oleson K W，Vertenstein M，et al.The LandSurface Climatology of the Community Land Model Coupled to the NCARCommunity Climate Model.2002，15(22)：3123-3149.Dai Yongjiu，Zeng Xubin，Dickinson R E，et al.The common land model[J].Bulletin of the AmericanMeteorological Society，2003，84(8)：1013-1023.)并借鉴TOPMODEL的产流思想推出了一个混合式产流方案，即基线径流方案，其产流公式如下：

R_s＝F_satQ_wat+(1-F_sat)w_s ⁴Q_wat    (3-1)

$R_{\mathrm{sb}}=F_{\mathrm{sat}}{l}_b{e}^{-D}+(1-F_{\mathrm{sat}})K_{\mathrm{sat}}\left(z_b\right){w_{\mathrm{bm}}}^{2B+3}+\underset{i=1}{\overset{10}{\mathrm{\Σ}}}\max [0,(({\mathrm{\θ}}_i-{\mathrm{\θ}}_e)\mathrm{\Δ}{z}_i/\mathrm{\Δt})]---(3-2)$

式中，F_sat为饱和面积百分比，Q_wat为净降水，w_s为土壤层深权重指数，l_b为最大地下径流系数，K_sat(z_b)是土壤底层的饱和水力传导率，θ_i为第i层的体积土壤水含量，θ_e是土壤有效空隙度，Δz_i和Δt是第i层的土壤深度和时间步长，D为缺水深，B为土壤质地构造等参数，W_bm为蓄水容量参数；

2)对CLM2.0的地表径流公式进行改进(Chen Ji，Kumar P.Topographic influence onthe seasonal and inter-annual variation of water and energy balance of basins in North America[J].Journal of Climate，2001，14：1989-2014.)，而对地下径流公式则进行了简化，得到如下公式：

R_s＝F_satQ_wat+(1-F_sat)max(0，(Q_wat-I_max))    (3-3)

$R_{\mathrm{sb}}=\frac{\mathrm{\α}{K}_{\mathrm{sat}}\left(0\right)}{f}{e}^{-{\mathrm{\λ}}_m}{e}^{-{\mathrm{fz}}_{\▿}}---(3-4)$

式中，I_max是土壤下渗能力，α各向异性系数，K_sat(0)为饱和水力传导度，

为计算格网的土壤平均缺水深，F_sat为饱和面积百分比；

3)在上述两种产流方案的基础上(Niu G Y，Yang Z L.The versatile integrator ofsurface and atmosphere processes(VISA)Part II：Evaluation of three topography-based runoffschemes[J].Global Planet Change，2003，38：191-208.Stieglitz M，Rind D，Famiglietti J，et al.An efficient approach to modeling the topographic control of surface hydrology for regional andglobal modeling[J].Journal of Climate，1997，10：118-137.)，对产流的方案简化如下：

R_s＝F_satQ_wat                      (3-5)

$R_{\mathrm{sb}}=R_{\mathrm{sb},\max }{e}^{-{\mathrm{fz}}_{\▿}}---(3-6)$

式中，R_sb，max是计算格网平均土壤缺水深为0时的最大地下径流量，f是土壤衰减系数，F_sat是由计算格网的地形特性和土壤湿度状态所决定的，通过双参数幂指数函数法来建立F_sat和地形信息以及土壤湿度之间的联系，其表达式如下：

$F_{\mathrm{sat}}=F_{\max }{e}^{-C_s{\mathrm{fz}}_{\▿}}---(3-7)$

式中，F_max为计算格网的最大饱和面积百分比，即流域计算格网上地形指数大于或等于均值的象素占整个流域总象素的百分比，可通过地形指数的空间累积分布曲线CDF和地形指数均值来推求，而CDF曲线由IMFD地形指数计算新方法得到的，地形指数均值是通过基于IMFD算法的地形指数尺度转换公式获取的(雍斌.陆面水文过程模型TOPX构建及其与区域气候模式RIEMS的耦合应用[D].南京：南京大学，2008.)，C_s为一系数，通过地形指数均值λ_m来获得，即C_s＝0.04λ_m+0.21。

步骤D：高分辨率的全分布式水文模型TOPX土壤湿度计算：

1)借鉴新安江三层土壤蒸散发计算模型(赵人俊.流域水文模型—新安江模型与陕北模型[M].北京：水利电力出版社，1984.)以及TOPX模型的土壤湿度计算方法(雍斌.陆面水文过程模型TOPX构建及其与区域气候模式RIEMS的耦合应用[D].南京：南京大学，2008.)按照土壤垂向分布的不均匀性将土壤分为三层，主要参数有上层、下层及深层土壤张力水容量WUM、WLM、WDM(WM＝WUM+WLM+WDM)，流域蒸散发折算系数K以及深层蒸散发系数C。计算需要输入的数据是蒸发皿实测水面蒸发量，计算输出的变量为上层、下层及深层实时土壤蒸发量EU、EL、ED(E＝EU+EL+ED)和实时土壤含水量WU、WL、WD(W＝WU+WL+WD)。三层土壤湿度和蒸散发动态变化计算原则是：当降水事件发生时，先补充上层最大蓄水容量，满足WUM值后再补充下层最大蓄水容量，依次类推当中层WLM满足后最后都补充到深层。同理，对蒸散发来说，上层按蒸散发能力蒸发，上层含水量不满足蒸发能力时，剩余蒸散发能力从下层蒸发，下层的蒸发量与蒸散发能力及下层蓄水量成正比，不足部分由下层蓄水量补给，当下层蓄水量不够补给时，用深层蓄水量补给。另外需要考虑土壤湿度变化模块，久旱之后，WL已经很小，而仍然没有下雨的情况。由于植物根系作用，深层水分将供给蒸散发。因此在下层土壤湿度蒸散发判断时，引入一个下层蒸发土壤植被根系影响系数C。当WL＞C×WLM表明下层土壤含水丰富，则按照正常情况蒸发；而当WL≤C×WLM时，则又分为两种情况，一种是下层仍然有水可以蒸发，但此时比较难以蒸发，用C值来控制蒸发的量，另一种则是下层已无水可蒸发，则到深层进行蒸发。这样模型就可以把植被根系对土壤湿度的影响有效考虑进去了。具体的计算见图4的蒸散发计算流程图。

2)产流量R计算，可以通过新安江模型中的蓄水容量曲线来获得。一般来说，流域内各点的蓄水容量并不相同，新安江三水源模型把流域内各点的蓄水容量概化为一条抛物线，即张力水蓄水容量曲线：

$\frac{f}{F}=1-{(1-\frac{\mathrm{WM}}{\mathrm{WMM}})}^B---(3-8)$

流域平均蓄水容量：

$\mathrm{WM}={\∫}_0^{\mathrm{WMM}}(1-\frac{f}{F})\mathrm{dWM}=\frac{\mathrm{WMM}}{(1+B)}\×(1-\mathrm{IMP})---(3-9)$

与流域初始平均蓄水量的相应的纵坐标A为：

$A=\mathrm{WMM}\×(1-{(1-\frac{W}{\mathrm{WM}})}^{\frac{1}{1+B}})---(3-10)$

当PE＝P-E＞0时，产流，否则不产流。产流时：

当PE+A＜WMM：

$R=\mathrm{PE}-\mathrm{WM}+W+\mathrm{WM}\×{(1-\frac{A+\mathrm{PE}}{\mathrm{WMM}})}^{1+B}---(3-11)$

当PE+A＞＝WMM：

R＝PE-(WM-W)          (3-12)

式中，WMM为流域内蓄水容量最大值；WM为流域内某点的蓄水容量；f为流域蓄水容量小于WM的流域面积；F为流域面积；B为蓄水容量曲线方次，反应蓄水容量分布的不均匀性。

通过三层土壤蒸散发模型的计算，获得土壤湿度的时空动态变量，根据平均蓄水容量和土壤湿度的差值求的流域网格土壤缺水深进而根据公式3-5和3-6计算出地表径流和地下径流。

步骤E)高分辨率的全分布式水文模型TOPX模型的汇流计算：

采用栅格水系滞时演算法汇流来进行汇流处理，通过考虑不同径流类型的沿程传输损失，获得流域出口的流量过程。计算步骤和方法如下：

1)从流域出口栅格向上游按广度优先算法遍历所有对出口有流量贡献的栅格，压入队列；

2)按队列的先后顺序计算栅格的水流路径长度。首先根据栅格DEM高程差，采用常用的D8算法确定出栅格的水流流向，接着按照栅格水流流向，计算本栅格的水流流经长度；然后按照队列的先后顺序，后面的栅格(上游栅格)的水流路径长度为栅格本身的水流流经长度与下游栅格水流路径长度之和；

3)遍历完队列中的所有栅格，完成流域的栅格水流路径长度计算；

4)通过栅格水流路径首尾栅格高程差除以栅格水流路径长度，得到本栅格的水流路径平均坡度；

5)当已知了栅格的水流路径长度和水流路径平均坡度后，栅格汇流时间就可以按下式计算得出：

${\mathrm{RT}}_{i,j}=\frac{L_{i,j}}{V_{i,j}}=\frac{L_{i,j}}{K_v\×\sqrt{S_{i,j}}}---(3-13)$

式中，i和j为栅格的行列号，RT为栅格平均汇流时间(hr)；L为栅格水流路径长度(m)；V为栅格水流平均流速；S为栅格水流路径平均坡度。由于V随降雨量、降雨强度、土壤下垫面性质以及时间而变，难以确定，因此利用S^0.5*K_v来代替计算(杨涛，张鹰，陈界仁等.基于数字平台的黄河多沙粗沙区分布式水文模型研究—以黄河岔巴沟流域为例[J].水利学报，2005，36(4)：456-459.)。K_v为速度参数，包含了糙率、水力半径等因素对水流的影响，运算时需要通过率定得到。根据降雨径流数据率定出流域地表径流、壤中径流和地下径流的速度参数，根据率定得到的速度参数，分别计算出临沂流域模拟中的地表径流、壤中径流和地下径流的栅格汇流时间；

6)在水源沿汇流路径到达(子)流域出口的途中会产生沿程的传输损失，包括沿程的下渗和蒸散发，传输损失按以下公式计算：

$Q_{\mathrm{out}}=\frac{Q_{i,j}}{K\·{\mathrm{RT}}_{i,j}^n}---(3-14)$

式中，Q_out为(子)流域出口处的汇流流量(m³·s^-1)；Q_i，j为分水源后的栅格径流流量(m³)；K为汇流传输损失系数，通过率定得到；n为汇流时间的指数，指数越大，沿程损失越厉害，通常取n＝1，表示汇流流量与汇流时间成倒数关系。

步骤F：模型测试结果及其分析：

为了检验模型的模拟结果的精度，通常采用Nash模型效率系数(确定性系数)(Nash J E，Sutcliffe J V，1970.River flow forecasting through conceptual model：adiscussion of principles[J].Journal of Hydrology，10：282-290.)来对模拟的径流量进行分析和评价。Nash系数的定义如下：

$\mathrm{NE}=1-\frac{{\mathrm{\Σ}}_{i=1}^n{(Q_{\mathrm{obs},i}-Q_{\mathrm{sim},i})}^2}{{\mathrm{\Σ}}_{i=1}^n{(Q_{\mathrm{obs},i}-{\stackrel{\‾}{Q}}_{\mathrm{obs}})}^2}---(3-15)$

式中，Q_obs，i为实测流量(m³·s^-1)，Q_sim，i为模拟流量(m³·s^-1)，

为实测流量均值(m³·s^-1)，n为模拟的时段数。确定性系数越接近1，表明模型效率越高。

3.有益效果

本发明与现有技术中的模型设计方法相比，其显著进步在于：

(1)本发明能够简单高效的处理DEM数据，提高地形指数的效率，解决TOPX模型在高分辨率研究流域情况下，DEM网格数量多，数据量大，地形指数的计算工作量非常巨大，处理大的DEM数据非常耗时，导致死机的问题，有效的提高了TOPX模型在在高分辨率的DEM区域的应用。

(2)本发明的模型能够完成土壤湿度过程模拟，定量分析土壤湿度过程及其动态变化，为人类减少旱涝灾害提供科学的参考依据，科学管理水资源的调度，为合理开发和利用水资源提供科学依据和合理建议；

(3)本发明的模型能够完成洪水过程模拟，定量分析和预测洪水过程及其影响，为人类减少洪水灾害提供科学的参考依据；

(4)本发明的模型可以模拟长期的降雨径流过程，定量分析和预测不同时间尺度下流域的水文水资源动态变化，为合理开发和利用水资源提供科学的参考依据。

附图说明

图1为本发明的高分辨率的全分布式水文模型开发原理框架图；

图2 DEM洼地及平坦区处理新算法；

图3新旧方法提取的流向及地形指数对比；

图4为本发明的高分辨率全分布式水文TOPX模型蒸散发计算原理框图；

图5为临沂流域栅格水流路径长度提取结果；

图6为临沂流域栅格的水流路径平均坡度提取结果；

图7为本发明对临沂流域2001-2003年率定期模拟的日降雨-径流过程线图示；

图8为本发明对临沂流域2001-2003年逐月径流量和月平均流量的结果图示；其中：(a)逐月径流量图；(b)月平均流量图；

图9为本发明对临沂流域2006-2007年验证期模拟的日降雨-径流过程线图示；

图10为本发明对临沂流域2006-2007年逐月径流量和月平均流量的结果图示；其中：(a)逐月径流量图；(b)月平均流量图。

具体实施方式

下面结合实施例和附图，对本发明做进一步的解释说明。

实施例

临沂流域，经度117.35°～118.74°，纬度34.98°～36.30°。地形起伏较大，海拔跨度从57m到1125m。临沂流域气候属于典型的暖温带大陆性季风气候。冬季寒冷干燥且多风，降水稀少；夏季受副热带高压影响，炎热多雨；春秋两季气候温和。据多年资料统计，流域内多年平均气温11.8℃～13.3℃，多年平均蒸发量839mm，多年平均降水量830mm，年降水量主要集中在夏季，暴雨多而强度大，是山东省的暴雨中心之一，年际变化也较大。多年平均天然河川径流量约为30.16亿立方米。流域的山地、丘陵、平原面积各占总面积的三分之一。土壤类型主要有淋溶褐土、棕壤性土、棕壤、潮土和水稻土。选取临沂流域为研究实施区域，利用本发明的高分辨率的全分布式TOPX模型，利用DEM数据驱动地形指数算法，利用临沂流域的水文气象数据，土壤类型矢量数据，土地利用类型栅格数据和数字高程模型数据，驱动高分辨率的全分布式水文模型TOPX模型，模拟临沂流域水文过程。

步骤一：高分辨率的全分布式水文模型TOPX模型数据准备

1)高分辨率的全分布式水文模型TOPX模型需要数字高程数据，包括临沂流域数字高程模型数据等；

2)高分辨率的全分布式水文模型TOPX模型需要临沂流域气象数据，临沂、费县等10个气象台站，2001-2003年，2006-2007年逐日气象数据，气象数据包括日最高气温、日最低气温、日降水量、日相对湿度、日均风速、日照时数和20cm小型蒸发皿蒸发量；

3)高分辨率的全分布式水文模型TOPX模型需要临沂流域土壤类型矢量数据，土地利用栅格数据；

4)高分辨率的全分布式水文模型TOPX模型需要临沂流域水文数据，水文数据包括临沂水文站2001-2003年，2006-2007年逐日径流量数据和次洪数据；

5)高分辨率的全分布式水文模型TOPX模型输入DEM数据进行地形指数提取；

6)高分辨率的全分布式水文模型TOPX模型输入临沂流域1)～4)准备数据，对模型地形指数提取，然后对模型进行率定及验证。

步骤二：高分辨率的全分布式水文模型TOPX模型地形指数提取

1)填充洼地的新算法

新填洼算法先建立两个栈(栈A和栈B)，把临时DEM中分位数最小的那些栅格的位置信息(i行j列)存储在栈A或栈B中，分位数最小的那一类中每个栅格(用c表示)进行循环处理，如果栅格单元属于洼地则令其等于其相邻栅格中最低的一个；若不是洼地则令其等于它对应的初始矩阵的值，并且将此栅格从栈中剔出，将没有被剔出的栅格存放到另一个栈中(栈B或栈A)。然后再对栈B或栈A进行类似的处理。重复此过程，直到栈中只剩下洼地为止，这些洼地的位置信息最终会保存在两个栈中的一个里。

第一类处理完后，在第二类处理之前还有一关键步骤：检查栅格中是否还存在数值等于该类分位数的栅格，如果存在，则将其赋值为下一类(即第二类)分位数，因为这些栅格虽然是洼地，但其在空间上与填洼后的第一类不相邻，需要进一步处理。之后进行第二类，第三类等填洼处理，直到所有类全部处理完，这时的CDEM就是经过填洼和平坦区处理后的DEM，此时所有的洼地对应的位置信息存放在栈A或栈B两个中的一个之中，这非常有利用随后的水文模拟与分析。

2)方法评价

提出的新方法用1°×1°USGS DEM数据(分辨率为90m)进行评价。选用的DEM大小为1×1042441个栅格到8×1042441个栅格(在空间上覆盖范围大约为10000km²到80000km²)，用这样大小的DEM已足够用来测试不同方法的效率。使用普通的900MHz Intel PentiumⅣ处理器、512M RAM内存的台式计算机对Jenson和Domingue(JENSON S K，DOMINGUE J O.Extracting topographicstructure from digital elevation data for geographic information system analysis[J].Photogrammetric Engineering and Remote Sensing，1988，54(11)：1593-1600.)的方法，Planchon和Darboux(Planchon O，Darboux F.A fast，simple and versatilealgorithm to fill the depressions of digital elevation models[J].Catena，2002，46(2-3)：159-176.)的方法及新方法进行效率对比测试图3所示。本文采用以下三个指标对算法的效率进行评价：

扫描栅格总数：完成DEM洼地和平坦区处理所需要扫描的总栅格个数。DEM中每个栅格要被扫描多次，扫描次数的多数与算法有密切关系。一般说来，扫描栅格总数愈大，算法效率愈低。

比较总次数：即栅格与它周围相邻的8个栅格进行比较的总次数。使比较总次数降低的算法往往是较优的算法。

执行时间：尽管这一指标会受到具体的编译器及计算机性能的影响，但它可作为算法效率高低的一个很有用的辅助评价指标。

新方法比传统的Jenson和Domingue(JENSON S K，DOMINGUE J O.Extracting topographic structure from digital elevation data for geographicinformation system analysis[J].Photogrammetric Engineering and Remote Sensing，1988，54(11)：1593-1600.)的方法效率提高了很多倍。比如，对于大小为8×1042441个栅格的DEM，用传统的Jenson和Domingue(JENSON S K，DOMINGUE J O.Extracting topographic structure from digital elevation data forgeographic information system analysis[J].Photogrammetric Engineering andRemote Sensing，1988，54(11)：1593-1600.)的算法需要15.2个小时，而用新方法只需要19秒！至于Planchon和Darboux(2001)的两个算法(直接算法和改进算法)，新方法仍然显示出了优越性：平均而言，新算法要比Planchon和Darboux(Planchon O，Darboux F.A fast，simple and versatile algorithm to fill thedepressions of digital elevation models[J].Catena，2002，46(2-3)：159-176.)的直接算法快几百倍以上，比其提高算法快7倍以上。

3)新地形指数提取方法

本文提出的新地形指数提取方法首先用分位数法将DEM分成8类，然后将DEM处理成无洼地及平坦区的DEM。然后，由于DEM中无洼地及平坦区，可直接根据D8算法指定水流方向(JENSON S K，DOMINGUE J O.Extracting topographicstructure from digital elevation data for geographic information system analysis[J].Photogrammetric Engineering and Remote Sensing，1988，54(11)：1593-1600.)。最后进行地形指数的计算，其具体计算方法参照Wolock and McCabe(Wolock D M，McCabe G J.Differences in topographic characteristics computed from 100-and1000-m resolution digital elevation model data[J].Hydrological Processes，2000，14(6)：987-1002.)的方法。

为了方便起见，挑选一大小为15×15栅格的DEM进行新旧地形指数计算方法的对比分析。如图3(b)&(c)所示，新旧方法提取的流向存在微小差别(以粗体箭头表示)，且差别主要产生在洼地及平坦区。相应地，新旧方法计算的地形指数也存在微小差别，且差别主要产生在洼地及平坦区，见图3(d)&(e)。这表明，虽然新旧方法在提取地形指数怕耗时间上存在很大差异，但提取的地形指数差别非常小，且差别主要产生在洼地及平坦区。

步骤三：高分辨率的全分布式水文模型TOPX模型产流计算

1)借鉴TOPX模型的产流计算方法(雍斌.陆面水文过程模型TOPX构建及其与区域气候模式RIEMS的耦合应用[D].南京：南京大学，2008.)，建立模型的产流方案。计算临沂流域产流，其产流公式(Niu G Y，Yang Z L.The versatileintegrator of surface and atmosphere processes(VISA)Part II：Evaluation of threetopography-based runoff schemes[J].Global Planet Change，2003，38：191-208.Stieglitz M，Rind D，Famiglietti J，et al.An efficient approach to modeling thetopographic control of surface hydrology for regional and global modeling[J].Journalof Climate，1997，10：118-137.)如下：

R_s＝F_satQ_wat                         (5-5)

$R_{\mathrm{sb}}=R_{\mathrm{sb},\max }{e}^{-{\mathrm{fz}}_{\▿}}---(5-6)$

式中，R_sb，max是计算格网平均土壤缺水深为0时的最大地下径流量，f是土壤衰减系数，F_sat是由计算格网的地形特性和土壤湿度状态所决定的。通过双参数幂指数函数法来建立F_sat和地形信息以及土壤湿度之间的联系，其表达式如下：

$F_{\mathrm{sat}}=F_{\max }{e}^{-C_s{\mathrm{fz}}_{\▿}}---(5-7)$

式中，F_max为计算格网的最大饱和面积百分比，即流域计算格网上地形指数大于或等于均值的象素占整个流域总象素的百分比，可通过地形指数的空间累积分布曲线CDF和地形指数均值来推求，而CDF曲线由IMFD地形指数计算新方法得到的，地形指数均值是通过基于IMFD算法的地形指数尺度转换公式获取的(雍斌.陆面水文过程模型TOPX构建及其与区域气候模式RIEMS的耦合应用[D].南京：南京大学，2008.)。C_s为一系数，通过地形指数均值λ_m来获得，即C_s＝0.04λ_m+0.21。

步骤四：高分辨率的全分布式水文模型TOPX土壤湿度计算

1)借鉴新安江三层土壤蒸散发计算模型(赵人俊.流域水文模型—新安江模型与陕北模型[M].北京：水利电力出版社，1984.)，以及TOPX模型的土壤湿度计算方法(雍斌.陆面水文过程模型TOPX构建及其与区域气候模式RIEMS的耦合应用[D].南京：南京大学，2008.)，建立模型的产流方案。

按照土壤垂向分布的不均匀性将土壤分为三层，主要参数有上层、下层及深层土壤张力水容量WUM、WLM、WDM(WM＝WUM+WLM+WDM)，流域蒸散发折算系数K以及深层蒸散发系数C。计算需要输入的数据是蒸发皿实测水面蒸发量，计算输出的变量为上层、下层及深层实时土壤蒸发量EU、EL、ED(E＝EU+EL+ED)和实时土壤含水量WU、WL、WD(W＝WU+WL+WD)。在下层土壤湿度蒸散发判断时，引入一个下层蒸发土壤植被根系影响系数C。当WL＞C×WLM表明下层土壤含水丰富，则按照正常情况蒸发；而当WL≤C×WLM时，则又分为两种情况，一种是下层仍然有水可以蒸发，但此时比较难以蒸发，用C值来控制蒸发的量，另一种则是下层已无水可蒸发，则到深层进行蒸发。这样模型就可以把植被根系对土壤湿度的影响有效考虑进去了。具体的计算见蒸散发计算流程图4；

2)产流量R计算。一般来说，流域内各点的蓄水容量并不相同，新安江三水源模型把流域内各点的蓄水容量概化为一条抛物线，即张力水蓄水容量曲线：

$\frac{f}{F}=1-{(1-\frac{\mathrm{WM}}{\mathrm{WMM}})}^B---(5-8)$

流域平均蓄水容量：

$\mathrm{WM}={\∫}_0^{\mathrm{WMM}}(1-\frac{f}{F})\mathrm{dWM}=\frac{\mathrm{WMM}}{(1+B)}\×(1-\mathrm{IMP})---5-9)$

与流域初始平均蓄水量的相应的纵坐标A为：

$A=\mathrm{WMM}\×(1-{(1-\frac{W}{\mathrm{WM}})}^{\frac{1}{1+B}})---(5-10)$

当PE＝P-E＞0时，产流，否则不产流。

产流时：

当PE+A＜WMM：

$R=\mathrm{PE}-\mathrm{WM}+W+\mathrm{WM}\×{(1-\frac{A+\mathrm{PE}}{\mathrm{WMM}})}^{1+B}---(5-11)$

当PE+A＞＝WMM：

R＝PE-(WM-W)          (5-12)

式中，WMM为流域内蓄水容量最大值；WM为流域内某点的蓄水容量；f为流域蓄水容量小于WM的流域面积；F为流域面积；B为蓄水容量曲线方次，反应蓄水容量分布的不均匀性。

通过三层土壤蒸散发模型的计算，我们可以获得土壤湿度的时空动态变量，根据平均蓄水容量和土壤湿度的差值求的流域网格土壤缺水深

进而根据公式5-5和5-6计算出地表径流和地下径流。然后根据KSS壤中流出流系数结合计算出的径流得出壤中流。

步骤五：高分辨率的全分布式水文模型TOPX模型的汇流计算

采用栅格水系滞时演算法汇流来进行汇流处理，通过考虑不同径流类型的沿程传输损失，获得流域出口的流量过程。计算步骤和方法如下：

1)从流域出口栅格向上游按广度优先算法遍历所有对出口有流量贡献的栅格，压入队列；

2)按队列的先后顺序计算栅格的水流路径长度。首先根据栅格DEM高程差，采用常用的D8算法确定出栅格的水流流向，接着按照栅格水流流向，计算本栅格的水流流经长度；然后按照队列的先后顺序，后面的栅格(上游栅格)的水流路径长度为栅格本身的水流流经长度与下游栅格水流路径长度之和；

3)遍历完队列中的所有栅格，完成流域的栅格水流路径长度计算(见图5)；

4)通过栅格水流路径首尾栅格高程差除以栅格水流路径长度，得到本栅格的水流路径平均坡度(见图6)；

5)当已知了栅格的水流路径长度和水流路径平均坡度后，栅格汇流时间就可以按下式计算得出：

${\mathrm{RT}}_{i,j}=\frac{L_{i,j}}{V_{i,j}}=\frac{L_{i,j}}{K_v\×\sqrt{S_{i,j}}}---(5-13)$

式中，i和j为栅格的行列号，RT为栅格平均汇流时间(hr)；L为栅格水流路径长度(m)；V为栅格水流平均流速；S为栅格水流路径平均坡度。由于V随降雨量、降雨强度、土壤下垫面性质以及时间而变，难以确定，因此利用S^0.5*K_v来代替计算(杨涛，张鹰，陈界仁等.基于数字平台的黄河多沙粗沙区分布式水文模型研究—以黄河岔巴沟流域为例[J].水利学报，2005，36(4)：456-459.)。K_v为速度参数，包含了糙率、水力半径等因素对水流的影响，运算时需要通过率定得到。根据降雨径流数据率定出流域地表径流、壤中径流和地下径流的速度参数，根据率定得到的速度参数，分别计算出临沂流域模拟中的地表径流、壤中径流和地下径流的栅格汇流时间；

6)在水源沿汇流路径到达(子)流域出口的途中会产生沿程的传输损失，包括沿程的下渗和蒸散发。传输损失按以下公式计算：

$Q_{\mathrm{out}}=\frac{Q_{i,j}}{K\·{\mathrm{RT}}_{i,j}^n}---(5-14)$

式中：Q_out为(子)流域出口处的汇流流量(m³·s^-1)；Q_i，j为分水源后的栅格径流流量(m³)；K为汇流传输损失系数，通过SQLoss地表径流沿程传输损失系数，LQLoss壤中径流沿程传输损失系数，BQLoss地下径流沿程传输损失系数获得，通过率定得到；n为汇流时间的指数，指数越大，沿程损失越厉害。通常取n＝1，表示汇流流量与汇流时间成倒数关系。

步骤六：模型测试结果及其分析

采用Nash模型效率系数(确定性系数)(Nash J E，Sutcliffe J V，1970.Riverflow forecasting through conceptual model：a discussion of principles[J].Journal ofHydrology，10：282-290.)来对模拟的径流量进行分析和评价，利用临沂水文站2001-2003年，2006-2007年逐日径流量数据对高分辨率的全分布式水文模型TOPX模型临沂水文过程的模拟结果进行详细分析，以评判该模型的模拟能力。

根据临沂流域2001-2003年水文、气象数据，率定出模型运行参数如表5-1所示。

根据表5-1率定的模型参数，模拟得到临沂流域率定期日降雨-径流过程如图7所示。

表5-1：临沂流域2001-2003年径流模拟参数率定结果

从图7模拟的日降雨-径流过程线可以看出：模型能够对日径流量进行很好的模拟，能够较好的捕捉洪峰出现的时间，但是对于洪峰量上的模拟有一定的误差，且总体偏小。模型对降水的丰沛程度比较敏感，在丰水年的模拟效果无论是洪峰流量还是峰现时间都要明显好于枯水年。在丰水年中，单峰值的丰水年(2001年)效果又要好于多峰值的丰水年(2003年)；在2002年多峰值的枯水年，模型模拟的流量过程线较差，除了对7月底的最大洪峰和4月底的次大洪峰捕捉的较好外，其它的次大洪峰(3月底、5月初和5月中下旬)都没有捕捉到；在流量过程线中，还有一些值得注意的现象：模型模拟的流量对降水有很好的响应，特别是在雨季响应较好，但对久旱后的降雨事件反应不是很好，比如枯水年的3月底、5月初和5月中下旬以及2003年丰水年的6月底的四次降雨过程都没有在流量过程线中体现出来，分析其主要原因为：久旱之后，土壤中的含水量较低，降雨主要用来满足三层土壤水分的蒸散发，短时间内没有形成径流。同时存在有降雨事件发生，但实测流量没有相应的显示出来，例如2001年的夏季初期6月份，主要原因为：一是短时间强度不大的降雨过程使得雨量不能马上形成汇流；二是人类活动对水文过程的影响，大型水库在夏季蓄水将对降雨时间后的流量产生较大的影响。

图8为2001-2003年逐月径流量和月平均流量得到的降雨-径流过程线结果图示。由图8(a)的逐月流量图和图8(b)的月平均流量图可以看出：模拟的流量和实测流量之间存在着较高的一致性，对洪峰流量和峰现时间都有较高的捕捉能力，模拟的流量对降水有非常好的响应，当有大的降雨过后，就有较高的洪峰流量出现。

进一步对长时段降雨-径流过程的模拟精度进行定量分析和评价。

表5-2：临沂流域率定期日、月降雨径流过程误差分析表(2001-2003)

模拟时间段	确定性系数	相关系数
			日模拟	0.790	0.890
月模拟	0.853	0.929

表5-3：临沂流域率定期月流量相对误差分析表(2001-2003)

表5-2的日、月降雨径流过程误差分析表可知：高分辨率的全分布式水文模型TOPX模型在日模拟上的确定性系数和相关系数分别为0.790和0.890，在月模拟上的确定性系数和相关系数分别为0.853和0.929；说明高分辨率的全分布式水文模型TOPX模型完全能够在日和月尺度上达到径流过程模拟和预测的要求。

为了验证模型在半湿润半干旱区域的适应性，采用率定期得到的模型率定参数，对临沂流域2006-2007年的数据进行了验证。验证期得到的日降雨-径流过程如图9所示。对验证期日模拟结果按月进行统计平均，得到如图10(a)所示的逐月径流量图，如图10(b)所示的月平均流量图。

从图9验证期模拟的日降雨-径流过程线可以看出：模型模拟的结果与验证期结果相似，对洪峰流量和峰现时间都有较好的捕捉能力。在验证期，单峰值的丰水年(2006年)效果次于多峰值的丰水年(2007年)。其中2006年，最大洪峰出现的时间比观测的洪峰有滞后现象，分析其原因可能为：由于流域上不同地方的降雨汇至流域出口需要一定的时间，因此洪峰出现的时间相对于降雨事件都有一滞后现象，而8月初的这次降雨事件，模型模拟峰现时间比较合理，出现在降雨事件之后，而实测峰现时间则和降雨出现在同一天，不符合常规，造成模型模拟效率系数较低，分析其过程线，其前期降雨并不是很多，强度也不大，因此最大可能原因是人为因素的影响，水库开闸放水造成。

图10的验证期逐月径流量和月平均流量得到的降雨-径流过程线也与验证期相似，模拟的流量和实测流量之间存在着较高的一致性，对洪峰流量和峰现时间都有较高的捕捉能力，模拟的流量对降水有非常好的响应。

表5-4：临沂流域验证期日、月降雨径流过程误差分析表(2006-2007)

模拟时间段	确定性系数	相关系数
			日模拟	0.724	0.869
月模拟	0.905	0.973

表5-5：临沂流域验证期月流量相对误差分析表(2006-2007)

验证期长时段降雨—径流过程的模拟精度进行定量分析和评价如表5-4和表5-5。由表可知：高分辨率的全分布式水文模型TOPX在验证期日模拟上的确定性系数和相关系数分别为0.724和0.869，在月模拟上的确定性系数和相关系数分别为0.905和0.973，总体来看和率定期接近，月尺度上的模拟效果稍好于日尺度上。月流量数据的相对误差统计分析表明：在冬春季节，实际流量值比较小，平均相对误差较大；在雨季(7-10月)，实测流量接近模拟流量，模拟精度较高，平均相对误差在整年中最小，约为16％左右；验证期的日、月误差的量值和误差分布规律也和率定期基本一致，体现出所率定的参数在临沂流域有一定的代表性，同时模型对半湿润半干旱区的水文过程模拟有良好的稳定性。

Claims (7)

1.一种高分辨率的全分布式水文模型TOPX的设计方法，其步骤为：

A）高分辨率的全分布式水文模型TOPX模型数据准备；

B)  高分辨率的全分布式水文模型TOPX模型地形指数计算；

C）高分辨率的全分布式水文模型TOPX模型产流计算；

D）高分辨率的全分布式水文模型TOPX土壤湿度计算；

E）高分辨率的全分布式水文模型TOPX模型的汇流计算。

2.根据权利要求1所述的一种高分辨率的全分布式水文模型TOPX的设计方法，其特征在于步骤A）中，所述高分辨率的全分布式水文模型TOPX模型的数据准备包括气象数据、水文数据和基础地理数据，其中，气象数据包括日最高气温、日最低气温、日降水量、日相对湿度、日均风速、日照时数和20cm小型蒸发皿蒸发量；水文数据包括逐日径流量数据和次洪数据；基础地理数据包括土壤类型矢量数据、土地利用类型栅格数据和数字高程模型数据。

3.根据权利要求2所述的一种高分辨率的全分布式水文模型TOPX的设计方法，其特征在于步骤B）中模型地形指数计算，计算步骤如下：

（1）输入DEM栅格数据( 除边界以外的所有点存储在二维数组A中 ) ；

计算出高程数据的8个分位数；

创建一个大小和A一样的数组B；

对于数组A中的每一个元素，如果其值小于第1分位数，则将数组B中对应位置的值全部赋值为第一分位数；如果其值小于第二分位数且大于第一分位数，则将数组B中对应位置的值全部赋值为第二分位数；… ；如果其值小于第八分位数且大于第七分位数，则将数组B中对应位置的值全部赋值为第八分位数；

创建两个完全一样的栈：栈A 和 B；

并将数组B中值为第一分位数的位置信息存储在栈A中；

（a）令 flag = False；

对于栈A 或栈B中的每一个元素；

对于任一栅格C，如果其在数组A中的值≥其在数组B中相邻的8 个栅格的值+δ(n)？

将栅格C在数组B中的值赋值为数组A对应位置的值, 并令 flag=True, 同时从栈A或栈B中删除栅格C；

（b）令 flag =ture

对于栈A 或栈B中的每一个元素

flag = False？

结束 

数组 B 中第一类填洼完成，其他7类方法相同 

（c）令 flag = False

对于栈A 或栈B中的每一个元素

对于任一栅格C，如果其在数组A中的值≦其

在数组B中相邻的8 个栅格的值+δ(n)？

对于任一栅格C，如果其在数组B中的值>其

在数组B中相邻的8 个栅格的值+δ(n)？ 

将栅格C在数组B中的值赋值为这一相邻栅格的值, 并令 flag=True, 同时将栈A/B中元素移到栈B/A中

（d）令 flag = False

对于栈A 或栈B中的每一个元素

对于任一栅格C，如果其在数组A中的值≦其

在数组B中相邻的8 个栅格的值+δ(n)？

对于任一栅格C，如果其在数组B中的值﹤其

在数组B中相邻的8 个栅格的值+δ(n)？

flag = False？

结束 

数组 B 中第一类填洼完成，其他7类方法相同 

（2）直接根据D8算法指定水流方向，计算地形指数:

T=ln(α/tanβ)，

式中 ln 是指自然对数，α 代表单位等高线长度上游集水面积，tanβ 表示坡度。

4.根据权利要求3所述的一种高分辨率的全分布式水文模型TOPX的设计方法，其特征在于步骤B）中高分辨率的全分布式水文模型TOPX模型产流计算，推出一个混合式产流方案，即基线径流方案，其产流公式如下：

                                                 

式中，F _sat为饱和面积百分比，为净降水，w _s 为土壤层深权重指数，l _b 为最大地下径流系数，K _sat (z _b )是土壤底层的饱和水力传导率，θ _i 为第i层的体积土壤水含量，θ _e 是土壤有效空隙度，

和

是第i层的土壤深度和时间步长，D为缺水深，B为土壤质地构造参数，W_bm蓄水容量参数。

5.根据权利要求4所述的一种高分辨率的全分布式水文模型TOPX的设计方法，其特征在于步骤C）建立土壤蒸散发计算模型，模型公式如下：

a.当PE>0

（1）WU+PE-R>WUM

WU+WL+PE-R-WUM>WLM

WUM=WU

WLM=WL

W+PE-R-WU-WL=WD

K*EM=EU

0=ED=EL

EU+EL+ED=E

WU+WL+WD=W

（2）WU+PE-R>WUM

WU+WL+PE-R-WUM<WLM

WU+WL+PE-R-WUM=WL

WUM=WU

K*EM=EU

0=ED=EL

EU+EL+ED=E

WU+WL+WD=W

（3）WU+PE-R<WUM

WU+PE-R=WU

K*EM=EU

0=ED=EL

EU+EL+ED=E

WU+WL+WD=W

b.当PE<0

（1）WU>PE

K*EM=EU

0=EL=ED

WU+PE=WU

EU+EL+ED=E

WU+WL+WD=W

（2）WU<PE

WU+P=EU

0=WU

WL>C*WLM

(K*EM-EU)*WL/WLM=EL

WL-EL=WL

0=ED

EU+EL+ED=E

WU+WL+WD=W

（3）WU<PE

WU+P=EU

0=WU

WL<C*WLM

WL>C*(K*EM-EU)=EL

WL-EL=WL

0=ED

EU+EL+ED=E

WU+WL+WD=W

（4）WU<PE

WU+P=EU

0=WU

WL<C*WLM

WL<C*(K*EM-EU)

WL=EL

0=WL

C*(K*EM-EU)=EL=ED

WD-ED=WD

EU+EL+ED=E

WU+WL+WD=W

式中：WUM、WLM、WDM(WM=WUM+WLM+WDM)为上层、下层及深层土壤张力水容量；WU、WL、WD为各层土壤含水量；W为总土壤含水量；EU、EL、ED(E=EU+EL+ED)为上层、下层及深层实时土壤蒸发量；EM为实际蒸散发能力；PE为超渗降雨；R为流域产流量；K为流域蒸散发折算系数；C为下层蒸发土壤植被根系影响系数。

6.根据权利要求1~4中任一项所述的一种高分辨率的全分布式水文模型高分辨率的全分布式水文模型TOPX的设计方法，其特征在于步骤D）采用栅格水系滞时演算法汇流来进行汇流处理，栅格汇流时间公式如下：

 

式中，i和j为栅格的行列号，RT为栅格平均汇流时间（hr）；L为栅格水流路径长度(m)；V为栅格水流平均流速；S为栅格水流路径平均坡度；利用V=S^0.5*K _v 来计算V；K _v 为速度参数，运算时通过率定得到。

7.根据权利要求1~5中任一项所述的一种高分辨率的全分布式水文模型TOPX的设计方法，其特征在于步骤E）通过考虑不同径流类型的沿程传输损失，获得流域出口的流量过程，按以下公式计算：

 

式中，

为（子）流域出口处的汇流流量（m³·s^-1）；

为分水源后的栅格径流流量（m³）；K为汇流传输损失系数，通过率定得到；n为汇流时间的指数，指数越大，沿程损失越严重，通常取n＝1，表示汇流流量与汇流时间成倒数关系。

Images