目前在皖南山区中小流域, 或是由于现有水文站点的分布密度和它们所记录的数据信息并不足以满足建立高精度水文模型的要求, 或是由于各种自然条件与人为因素的限制, 一些基础性的水文数据都无法取得^[1].面对这样的事实, 水文工作者不得不开展无资料地区的水文模型研究以获得能够满足水资源管理需求的相关数据, 实现对广大区域水资源的科学调度管理^[2-5].

为了解决上述问题, 系统分析的许多概念、方法被逐步引入, 主要有单位线地区综合法、地貌瞬时单位线法和推理公式法等^[6].本文通过回归分析的手段建立基于DEM的Nash单位线参数推求公式, 并通过通径分析的手段揭示对流域汇流过程有影响的地貌因子及其作用机制, 为Nash单位线在中小河流无资料地区的应用提供参考.

1 数据与方法 1.1 数据来源

为了率定出Nash单位线参数, 本文选取皖南山区29个中小流域的实测水文资料, 资料来源于安徽省水文局; 为了提取流域的地形地貌因子和地表覆盖类型, 本文选用皖南山区中小流域的DEM数据, 数据来源于中国科学院计算机网络信息中心国际科学数据镜像网站, 精度为90 m×90 m, 数据采集时间为2000年; LandCover数据来源于国家基础地理信息中心, 精度为30 m×30 m, 数据采集时间为2010年.

1.2 模型原理及分析方法 1.2.1 Nash单位线

Nash^[7]提出流域对地面净雨的调蓄作用可用n个串联的线性水库的调节作用来模拟假设, 并由此推导出Nash瞬时单位线的数学表达式:

$ u\left( t \right) = \frac{1}{{k\left( {n - 1} \right)!}}{\left( {\frac{t}{k}} \right)^{n - 1}} \cdot {\rm{e}^{ - \frac{t}{k}}} $

(1)

式中, u(t)为瞬时单位线; n为线性水库个数; k为线性水库的蓄量常数, 表征流域平均汇流时间.

瞬时单位线转化为时段单位线采用S曲线, 瞬时单位线转换为时段为Δt、净雨量为10 mm的时段单位线计算公式为:

$ S\left( t \right) = \int_0^t {u\left( {t, 0} \right) \cdot {\rm{d}}t = \int_0^t {\frac{1}{{k\left( {n - 1} \right)!}}{{\left( {\frac{t}{k}} \right)}^{n - 1}} \cdot {{\rm{e}}^{ - \frac{t}{k}}} \cdot {\rm{d}}t} } $

(2)

$ u\left( {\Delta t, t} \right) = S\left( t \right) - S\left( {t - \Delta t} \right) $

(3)

$ q\left( {\Delta t, t} \right) = \frac{{10F}}{{3.6\Delta t}} \cdot u\left( {\Delta t, t} \right) $

(4)

式中, q(Δt, t)为时段单位线, 单位为m³/s; F为流域面积, 单位为km².

由公式(1)~(4)可以看出, 推求Nash单位线的关键是参数n和k的确定.

1.2.2 多元线性回归分析

多元线性回归分析是确定两种或两种以上变量间相互依赖的定量关系的一种统计分析方法.它基于观测数据建立变量间适当的依赖关系, 分析数据内在规律, 并可用于预报、控制等问题.在大多数实际问题中, 影响因变量的因素不是一个而是多个, 可以建立因变量y与各自变量x_i之间的多元线性回归模型^[8]:

$ y = {b_0} + {b_1} \cdot {x_1} + {b_2} \cdot {x_2} + \cdots +{b_k} \cdot {x_k} + e $

(5)

式中, b为回归系数, e为随机误差.

除了线性关系, 进行多元线性回归分析还需要其他3个假设:

a) 独立性:自变量不能有多重共线性; 在多元回归分析中, 由于涉及多个自变量, 所以需要考虑多重共线性.多重共线性是指线性回归模型中自变量之间由于存在较高的相关而使模型估计不准确或难以估计.因此, 在进行分析解释之前, 需要诊断共线性^[9].常用的参数及其标准有:①容差(tolerance) < 0.1, 则可能存在共线性问题; ②方差膨胀因子(variance inflate factor, VIF)实际上就是容差的倒数, < 10, 说明多重共线性弱; ③条件指数(Condition Index, CI) < 10, 表明存在的多重共线性弱.

b) 正态性:自变量任意线性组合, 残差服从正态分布.

c) 方差齐性:自变量任意线性组合, 因变量方差不变, 残差齐性.

建立回归方程之后, 还需要检验方程的有效性^[10].通过方差分析(F检验)来检验回归方程的显著性, 代表的是对进行回归的所有自变量的回归系数的总体检验; 通过t检验来说明偏回归系数的显著性; 通过确定性系数R²来说明回归方程对因变量的解释程度, 调整R²是针对自变量的增多会不断增强预测力的一个矫正, >0.6表示模拟效果很不错.

1.2.3 通径分析

通径分析是数量遗传学家Sewall Wright于1921年提出来的一种多元统计技术.通径分析可以通过对自变量与因变量之间表面直接相关性的分解, 来研究自变量对因变量的直接重要性和间接重要性, 从而为统计决策提供可靠的依据, 在众多领域得到广泛应用^[11].

通径分析在多元回归的基础上将相关系数r_iy分解为直接通径系数P_iy(某一自变量对因变量的直接作用)和间接通径系数(该自变量通过其他自变量对因变量的间接作用)^[12-13].通径分析的理论已证明, 任一自变量X_i与因变量y之间的简单相关系数(r_iy)=P_iy+所有X_i与y的间接通径系数, 任一自变量X_i与因变量y的间接通径系数=r_iy·P_iy.在通径分析过程中, 一般认为最难计算的就是通径系数.事实上, 通过软件进行线性回归计算, 计算结果给出的线性回归方程的标准系数(standardized coefficients)即为我们需要的通径系数.

以Nash单位线参数为因变量, 地形地貌因子和地表覆盖类型为自变量, 建立多元线性回归模型, 并就该回归模型进行通径分析.

2 结果与讨论 2.1 Nash单位线参数推求

选取皖南山区29个中小流域, 根据实测水文资料率定出Nash单位线参数.率定方法采用矩法初估, 人工微调.皖南山区29个中小流域的参数n均为3(表 1).早在1977年, Singh^[14]根据大量的实测资料分析, 曾得出不论流域面积多大, n可近似取作3的结论.因此, 本文对参数n不做过多的分析, 主要对参数k进行分析.

2.2 地形地貌因子和地表覆盖类型提取

利用ArcGIS 10.1提取皖南山区29个中小流域的11个地形地貌因子和10种地表覆盖类型. 11个地形地貌因子分别是:流域面积、面积坡度、河道平均坡度、河网总长度、河网密度、主河道长、流域平均高程、形状系数、流域长度、流域平均坡度和河道几何参数; 10种地表覆盖类型分别是:耕地、森林、草地、灌丛地、湿地、水体、苔原、人造森林、裸地、冰川与永久积雪.部分提取结果如表 2.

2.3 参数k的多元回归分析

从29个中小流域中选取26个进行多元回归分析, 得出回归方程.用剩下的3个流域(三元、黄尾河、杨山岭)对所得方程进行检验.

2.3.1 多元回归方程的建立与分析

根据芮孝芳利用地形地貌资料确定Nash模型参数的研究^[15]可知, k可以表示为地形地貌因子的幂函数形式.为便于做多元线性回归分析, 将k值、地形地貌因子和地表覆盖类型数值均取ln值, 这样可以满足自变量与因变量之间具有线性关系的条件, 这是进行线性回归分析的前提.

为了满足独立性要求, 本文采用的变量筛选方法为逐步回归.逐步回归的基本思想是将变量逐个引入模型, 每引入一个解释变量后都要进行F检验, 并对已经选入的解释变量逐个进行t检验, 当原来引入的解释变量由于后面解释变量的引入变得不再显著时, 则将其删除.这是一个反复的过程, 直到既没有显著的解释变量选入回归方程, 也没用不显著的解释变量从回归方程中剔除为止.依据上述思想, 可利用逐步回归筛选并剔除引起多重共线性的变量.

关于正态性和方差齐性的检验, 在回归分析后进行.对回归模型显著效果解释如下:

1) 该模型调整R²=0.796>0.6, 表明该回归模型对因变量的解释程度为79.6%, 模拟效果很好.

2) 回归方程的显著性水平sig < 0.05, 在F检验中拒绝原假设(所有自变量对因变量不能产生显著性影响), 说明至少存在1个自变量对因变量产生显著性影响, 即该模型能够有效预测因变量.

3) 各自变量VIF < 10(表 3), 故认为该回归模型不存在多重共线性的问题.

4) 除个别观测值外, 残差基本服从正态分布(图 1); 除金家流域外, 标准化残差均在[-2,2]之间波动(图 2), 与自变量无关, 满足方差齐性假设.

综上分析, 该回归模型服从上文提到的4个假设, 且回归效果显著.根据系数B(表 3), 可得k的推求公式:

$ \ln \left( k \right) = 1.909 - 0.206\ln \left( {{S_{{\rm{area}}}}} \right){\rm{ + }}0.157\ln \left( {{F_{{\rm{shape}}}}} \right){\rm{ + }}0.096\ln \left( {{G_{{\rm{river}}}}} \right){\rm{ + }}0.271\left( {{C_{{\rm{forest}}}}} \right) $

(6)

式(6)两边同时取e底, 得:

$ k = \frac{{6.746{F_{{\rm{shape}}}}^{0.157} \cdot {G_{{\rm{river}}}}^{0.096} \cdot {C_{{\rm{forest}}}}^{0.271}}}{{{S_{{\rm{area}}}}^{0.206}}} $

(7)

式中, F_shape为形状系数, G_river为河道几何参数, C_forest为森林覆盖率, S_area为面积坡度.

流域汇流研究的是流域对注入其中的净雨的再分配作用.近代对于这种再分配机理可归结为由于两种扩散作用^[16]:地貌扩散作用和水动力扩散作用.地貌扩散作用与流域的面积、形状和水系分布特点等有关.水动力扩散对地面径流汇流流域的空间分布不均匀主要表现在流速沿水流方向分布的不均匀性, 主要取决于流域的地形坡度沿水流方向的不均一性.式(7)中的面积坡度、形状系数、河道几何参数以及森林覆盖率这4个因子反映了流域的面积、坡度、形状和水系分布特点, 通过影响地貌扩散和水动力扩散两种作用来影响该流域净雨的再分配作用, 从而影响Nash单位线参数k的大小.因此, 从物理机制上分析, 所得k的推求公式合理.

第26个流域(金家流域)的率定k值和推求k值相差较大(图 3).与图 2所得结果一致, 所得回归方程对金家流域的拟合效果不好.一方面, 金家流域洪水场次资料太少, 大洪水只有一场, 利用API模型对该流域k值的率定存在误差; 另一方面, 金家流域下游有东方红水库, 对径流有拦蓄作用.由于在利用实测资料对k值进行率定时, 实测洪水资料已包含了水库的调蓄作用, 而在用根据地形地貌因子和地表覆盖类型得出的回归方程对该流域k值进行推求时未考虑该水库的影响, 从而导致率定值比推求值大很多.

2.3.2 多元回归方程的检验

将三元、黄尾河、杨山岭3个流域的面积坡度、形状系数、河道集合参数及森林覆盖率的值带入公式(7), 得到3个流域预测的k值.

推求的k值减去率定的k值为预测误差, 其绝对值为绝对误差.许可误差是依据预测精度的使用要求和实际预测技术水平等综合确定的误差允许范围. k值的预测以率定k值的20%作为许可误差.本研究中3个流域的预测结果均在允许误差之内, 检验结果均为合格(表 4).该结果表明模型的模拟效果较好, 该回归方程适用于皖南山区中小河流域, 具有较强的实用性, 可为其他无资料地区的水文预报提供借鉴.

2.4 参数k的通径分析

通径分析是在回归分析的基础上, 将简单相关系数分解为直接通径系数(某一自变量对因变量的直接作用)和间接通径系数(该自变量通过其他自变量对因变量的间接作用).

面积坡度对Nash单位线参数k的直接影响最大, 森林覆盖率对k的直接影响最小, 形状系数和河道几何参数对k的直接影响较小(表 5).但形状系数和森林覆盖率通过面积坡度对k产生较明显的间接影响.值得一提的是, 由于植被覆盖对流域汇流过程有阻滞作用, 会导致k增大, 所以森林覆盖率本身对k的直接影响为正方向; 但是在山区, 森林覆盖率高的流域一般在上游, 面积坡度大, 流域汇流过程更快, 导致k减小, 所以森林覆盖率对k的最终影响为负方向.这也就解释了森林覆盖率对k的直接通径系数为正, 而简单相关系数为负的现象(表 5).

3 小结

本文为解决Nash模型在无资料地区的应用问题提出了新思路, 即通过数学分析方法定量揭示地形地貌因子和地表覆盖类型对Nash模型参数k的影响机制, 同时建立模型参数与地貌因子之间的数学公式, 并进行检验.

利用多元线性回归方法研究基于地貌因子的参数k的数学推求公式.为消除多重共线性的影响, 采用了逐步回归的变量筛选方式.回归分析结果表明预测残差服从正态分布, 且不受自变量影响, 方差齐性.回归模型各项指标良好, 回归效果显著, 具有统计学意义.

利用通径分析方法定量研究各地貌因子对参数k的影响机制.基于回归分析得到的简单相关系数矩阵和回归系数矩阵, 计算各自变量对k的直接通径系数和间接通径系数.通径分析结果揭示了面积坡度对参数k的直接影响最大, 同时它通过其他因子(森林覆盖率等)对参数k有着明显间接影响.