山西省煤矿区地表水文过程的空间模拟研究

2018-08-02张吴平贾若男

中国农村水利水电 2018年7期

侯帅，张吴平，贾若男

(山西农业大学资源环境学院，山西太谷 030801)

0 引言

长期以来，煤炭资源都是中国工业发展的重要支撑力量，特别是20世纪80年代以来，煤炭资源为我国的经济建设和发展做出了巨大的贡献，但长期的煤炭开采造成水循环途径变化、土地和植被损毁、生态环境退化[1]。煤炭开采改变了地表形态，破坏了地表原有的水文过程，研究煤炭开采对地表水文过程的影响，可为煤矿区的环境保护、生态重建、科学开采等提供重要的理论依据。

自20世纪50年代以来，水文模型的研究取得了巨大的进展，其中分布式水文模型可更好地分析自然和人为因素对流域水环境的影响[2]，通常水文模型对水文过程的模拟主要包括植被截留、土壤入渗、径流汇流等过程，需要大量的水文、气象观测资料，而无观测站且受采煤扰动的长河流域难以满足这些要求。

植被截留降雨是水文循环中极其重要的一部分，植被截留受到降雨强度、降雨历时、植被类型、环境因素等诸多因素的影响，是一个复杂的混合过程[3]。国外的Horton、Rutter和Gash模型的应用较多，国内研究和应用较多的多为半经验性理论模型，如李崇巍、李玉霞等以仪垂祥模型为基础，建立了岷江上游林冠截持降雨遥感模型，结果较为理想[4-8]。但对农作物冠层截留的研究相对于林冠截留的研究较少且其结果差异较大[9]。SWAP模型中的作物截留计算采用了由Von Hoyningen-Hune和Braden提出的农田作物冠层截留量的一般方程，并且也被应用到其他研究水分运动的模型中去，如张圣微采用该方法建立了农田水分运动模型的降雨截留子模块[10]。大多数水文模型中的植被截留过程模拟以经验模型为主，且没有区分天然植被和种植植被，对此过程的研究较为简单。

土壤入渗指降雨或灌溉的水分从地表渗入土壤的运动过程，土壤入渗在水文循环过程中起着重要的纽带作用。土壤入渗的研究模型主要有经验模型，通过土壤入渗曲线得到模型参数如Kostiakov模型；半经验模型，采用简化的连续方程如Hortan、Holtan等模型；物理模型，基于Darcy定律和Richard方程，具有明确的物理意义，如Green-Ampt、Philip等模型,其中Green-Ampt模型广泛应用于水土侵蚀、降雨入渗等领域的研究中，如Mein 等将其改进应用于恒定降雨条件下的土壤入渗研究，也被应用到SWAT模型中计算下渗量[11,12]。

降雨经植被截留、土壤入渗后，在地形作用下形成地表径流。而数字高程模型(DEM)可反映地形特征，借助一定的算法可获取水系特征，为径流模拟提供了基础数据，因此已被应用到SWAT、SHE、TOPMODEL等诸多水文模型中[13]。

笔者以30 m×30 m的栅格为计算单元，基于3S技术构建了长河流域地表水文过程模型，分析了煤炭开采下模型参数的空间分布特征，并模拟了不同降雨强度下的地表水文过程。

1 材料与方法

1.1 研究区概况

山西省长河流域位于晋城市泽州县(图1)，地理坐标北纬35°30′～35°38′，东经112°37′～112°46′，包括47个村庄，面积约为113 km2，海拔725～1 167 m，是典型的黄土丘陵区。气候属暖温带大陆性季风气候，年平均气温11.2 ℃，年平均降雨量576 mm，长河于中部自北向南流过，长河以西分布着年产均在45 万t以上的12座煤矿，长期的煤炭开采造成土地破坏、水资源污染、空气污染等一系列问题，已在很大程度上影响了当地人民的生产生活和社会经济的健康发展。

图1 研究区地理位置及高程Fig.1 Geographical location and elevation of the study area

1.2 地表水文过程模型

研究区地处黄土高原，进行短时间径流模拟可以不考虑蒸散发和壤中流的影响，因此本文重点模拟地表水文过程中的植被截留降雨、土壤入渗、径流汇流过程。地表水文过程模型建立的条件为：以分水岭划定研究区边界，研究区内居民以井水为主要水源，在情景模拟中无外界河流汇入研究区。

本文构建的地表水文过程模型包括植被截留模型、土壤入渗模型、径流汇流模型，模型的总体结构及各模型间的联系如图2所示。由给定的降雨类型，首先根据各类型植被的截留降雨模型可得到截留的降雨量，落到地表的降雨再由土壤入渗模型得到入渗量，最后地表净雨量经地表汇流过程模拟得到汇流累积量的空间分布。本文首次在长河流域构建地表水文过程模型，由于研究区内缺乏水文观测资料，模型具体应用于下述降雨情景模拟中。

图2 地表水文过程模型结构示意图Fig.2 Schematic of the surface hydrological process model structure

1.2.1 植被截留模型

(1)天然植被截留模型。本部分采用仪垂祥提出的植被截留模型，对于给定的植被类型，植被冠层存在一个最大截留量，当降雨量P(mm)大于最大截留量时，最大截留降雨量PJ(mm)为[7]：

PJ=αVFCLAI

(1)

式中：α为叶表面最大持水深度。mm；VFC为植被覆盖度；LAI为叶面积指数，m2/m2。

当降雨量小于最大截留量时，植被截留量由植被覆盖度决定，其截留量PJ为：

PJ=VFCP

(2)

(2)种植植被截留模型。根据SWAP模型中计算作物截留的一般方程来模拟作物的截留量，随着降雨量的增加，其截留量也逐渐达到最大值aLAI[14]，截留量可由下式计算：

(3)

式中：PJ为降雨截留量，mm；a为经验系数；b为植被覆盖度。

1.2.2 土壤入渗模型

采用以下由Mein和Larson将Green-Ampt模型推广至降雨入渗情况的模型来对入渗过程进行模拟[15]。

假定有稳定雨强p，在降雨的开始阶段，降雨全部渗入土壤，当p大于土壤的入渗能力时，地表开始形成积水，由模型可得出积水时的累积入渗量Fp(cm)：

Fp=Sf(θs-θi)/(p/Ks-1)

(4)

由此可得开始积水时间tp=Fp/p，只有当p>Ks时才可能积水，所以整个过程的入渗率f(cm/d)为：

(5)

式中：Ks为饱和导水率，cm/d；Sf为湿润峰处的土壤水吸力值，cm；θs、θi分别为土壤的饱和含水率和初始含水率，cm3/cm3；F为开始积水以后的累积入渗量，cm。F的计算采用下面修正后的公式：

t>tp

(6)

式中：ts表示假设从t=0时开始积水，到入渗量为Fp时所需的时间，计算方法如下：

(7)

在实际应用中，首先根据降雨强度和降雨时间判断地表是否可以形成积水，积水前的累积入渗量由式(4)得到，积水后的累积入渗量由式(6)、(7)得出。

1.2.3 径流汇流模型

建立基于DEM的径流汇流模型，DEM填洼处理后，采用单流向算法中应用较多的D8法确定水流方向，根据合适的阈值提取出河网，应用Strahler法对河网分级，并将整个流域划分成若干子流域，从而建立研究区的数字水系[16]。

以植被截留、土壤入渗过程得到的地表净雨量为赋值栅格，对生成的数字水系赋值可得到地表汇流量的空间分布及各子流域的汇流累积量。

1.3 数据来源与处理

1.3.1 植被截留模型参数

(1)遥感数据。来源于美国地质勘探局(http:∥www.usgs.gov/)的2000-2015年的Landsat影像，首先需要对原始影像进行预处理，主要是辐射校正和几何校正。

(2)植被类型。采用研究区遥感影像结合当地土地利用规划图及《泽州县林地保护利用规划》对遥感影像进行解译，并利用研究区1∶2 000正射影像对地物判别和修正得到。

(3)叶表面最大持水深度。根据各乡镇的统计资料及文献[17]，并结合实地调查确定。

(4)植被覆盖度。利用获取的遥感影像来提取归一化植被指数(NDVI)来估算[18]。

(5)叶面积指数。通过遥感影像结合植被冠层多次散射模型来获取，具体步骤见文献[19]。

1.3.2 土壤入渗模型参数

模型的关键参数土壤饱和导水率由土壤理化性质通过van Genuchten模型得到，湿润峰处的土壤水吸力值根据张光辉等的研究得出[20,21]。所需的土壤的各项理化性质通过GPS确定1 km×1 km的采样网格对土壤采样后测定，共得到117组采样点数据。

1.3.3 径流汇流模型参数

DEM数据来源于泽州县国土资源局，分辨率为10 m，在ArcGIS中经填洼、确定流向、提取河网、划分子流域后可生成研究区的数字水系，其中含184个子流域，各径流于研究区中部汇流，总体从东北流向西南，经流域西南出口流出。

2 结果与讨论

2.1 参数的空间分布特征

2.1.1 植被覆盖度的空间分布特征

植被覆盖度是影响植被截留的重要因素，也反映了煤炭开采对植被的影响程度。2015年研究区的植被类型的空间分布如图3，天然植被类型主要为暖温带落叶阔叶林和针阔混交林，农作物主要是玉米、小麦。草地、农作物的分布基本相等，灌木林地、阔叶林、针阔混交林在矿区的面积远大于非矿区。

图3 研究区植被类型分布图Fig.3 Distribution map of vegetation types in the study area

将2000-2015年的植被覆盖度分为5级：①无覆盖，VFC≤0.1；②低覆盖，0.1

图4 2000-2015年研究区植被覆盖度等级的空间分布图Fig.4 Spatial distribution of vegetation coverage grades in the study area during 2000-2015

图5 2000-2015年矿区各植被覆盖度等级的面积Fig. 5 Area of the vegetation coverage grades in mining area during 2000-2015

表1 2000-2015年矿区植被覆盖度等级转移概率矩阵Tab.1 Transfer probability matrix of vegetation coveragegrades in mining area during 2000-2015

2.1.2 土壤饱和导水率的空间分布特征

土壤饱和导水率反映了土壤的渗透能力，是土壤水分运动的重要特征之一。表2是采样点土壤饱和导水率Ks的基本统计特征，矿区Ks的最小值、最大值、均值均小于非矿区，且Ks属中等变异，变异性较强，但矿区变异系数大于非矿区，说明受采煤影响，矿区Ks的变化情况大于非矿区。经Kolmogorov-Smirnov法检验，表明土壤饱和导水率符合正态分布(P=0.121>0.05)，通过地统计学方法对数据进行半变异函数的拟合，并通过交叉验证检验其合理性和拟合精度，选取决定系数较大(0.568)、残差较小(1.108×10-3)的最优拟合函数模型为指数模型，根据Kriging插值方法可得到土壤饱和导水率的空间分布(图6)[23]。研究区土壤饱和导水率的范围为4.337～139.321 cm/d，西部矿区的平均土壤饱和导水率34.375 cm/d小于东部非矿区的44.426 cm/d，土壤饱和导水率总体上呈现出西部小于东部的趋势。说明煤炭开采改变了原有的土壤性质，土壤饱和导水率变小。

表2 矿区、非矿区土壤饱和导水率的基本统计特征Tab.2 Basic statistics of soil saturated hydraulic conductivityin mining area and non-mining area

图6 研究区土壤饱和导水率的空间分布图Fig.6 Spatial distribution of soil saturated hydraulic conductivity in the study area

2.2 基于多年降雨资料的降雨情景模拟

2.2.1 长河流域降雨趋势分析

分析来源于泽州县气象中心的流域1986-2015年间的降雨资料可知，降雨主要集中在7、8月，两月降雨量占年总降雨量比例的均值为41.524%。按照中国气象局降雨量等级的划分标准，降雨类型主要为小雨(0.1～9.9 mm)和中雨(10～24.9 mm)，分别占67.982%、20.175%，其余大雨(25～49.9 mm)占9.357%，暴雨(50～99.9 mm)、大暴雨(100～249.9 mm)共占2.486%。

2.2.2 情景模拟分析

在上述降雨资料分析的基础上，以降雨强度分别为2.5、8、16 mm/h，降雨时间为1 h的降雨情景对研究区2015年7月的地表水文过程进行空间模拟。对于有植被覆盖的栅格，首先根据植被截留模型判断3种降雨情景下的降雨量是否达到各类型植被的最大截留降雨量：截留量未达最大截留量的栅格，降雨将不会落到地表；超过最大截留量的栅格根据土壤入渗模型中地表积水形成的时间，判断落到地表的降雨在1 h内是否会形成地表积水，形成积水的栅格再由径流汇流模型进行汇流累积量的计算。对于无植被覆盖的栅格，降雨直接落到地表，同理根据土壤入渗模型判断降雨在1 h内是否形成地表积水，是否进行径流汇流的计算。计算结果表明：

(1)2.5 mm/h降雨强度下，降雨大部分被植被截留，较少降雨渗入土壤，如图7(a)、(b)植被冠层平均截留量为0.160 mm，土壤平均入渗量分别为1.955 mm；

(2)8 mm/h降雨强度下，1 h内植被冠层截留量达到最大，其中农作物平均截留量为0.194 mm；阔叶林平均截留量为0.656 mm；针阔混交林平均截留量为0.442 mm；灌木林平均截留量为0.331 mm；草地平均截留量为0.197 mm。由于降雨强度小于土壤饱和入渗率，其余降雨经土壤入渗后难以形成地表径流，如图7(c)，(d)植被冠层平均截留量为0.198 mm，土壤平均入渗量分别为6.672 mm；

(3)16 mm/h降雨强度下，1 h内植被冠层截留量达到最大后，其余降雨一部分渗入土壤，一部分形成积水，植被截留量和土壤入渗量的空间分布如图7(e)、(f)，土壤平均入渗量为13.457 mm。积水继而形成径流，经径流汇流模拟可得到研究区地表径流汇流量的空间分布(图8)，汇流量由两侧向中部、由北向南逐渐增加，在流域西南即第184号子流域出口达到最大值。

图7 不同降雨强度下研究区植被截留量、土壤入渗量的空间分布图Fig.7 Spatial distribution of vegetation interception and soil infiltration in the study area under different rainfall intensities