资助项目:西北旱区生态水利国家重点实验室“黄河中游典型流域植被变化对蓝水—绿水转化的影响机制”(2019KFKT-3)
第一作者:杜婷婷(1987—),女,河北沧州人,硕士,讲师,研究方向:水文水资源。E-mail:498921822@qq.com 通信作者:郭梦京(1986—),男,山西运城人,博士,副教授,研究方向:流域水循环与水环境保护。E-mail:147433250@qq.com
(1.绵阳职业技术学院,四川 绵阳621000; 2.西安理工大学 水利水电学院,西安 710077; 3.江西省水利科学研究院,南昌330029; 4.江西省鄱阳湖流域农业资源与生态重点实验室,南昌 330045)
(1.Mianyang Vocational and Technical College,Mianyang,Sichuan 621000,China; 2.College of Water Conservancy and Hydropower,Xi'an University Of Technology,Xi'an 710077,China; 3.Jiangxi Institute of Water Sciences,Nanchang 330029,China; 4.Key Laboratory of Poyang Lake Watershed Agricultural Resources and Ecology of Jiangxi Province,Nanchang 330045,China)
Xijiang River Basin; VIC model; model application; hydrological elements; change trend analysis
水文模型[1]是用数学公式和物理性质来描述流域自然水循环和能量交换过程,随着计算机科学的不断发展,水文模型被广泛的应用在水文模拟预测、旱涝灾害预警、流域水资源管理等各个方面。得益于不断进步的地理信息系统(GIS)[2]技术与遥感技术(RS)[3],结合了GIS与RS技术的分布式水文模型不但考虑了传统水文模型所需要的降雨数据还重点考虑到了下垫面以及气候的变化对于流域水文生态的影响,使得分布式水文模型在流域自然水循环模拟过程中更具优势。
珠江流域主要支流有三条,分别是东江、北江以及西江,其中西江是珠海、澳门等地区城市用水最重要的上游来水源[4]。纵观前人的研究,对于珠江流域的水文研究主要集中在流域面积较小的北江和东江,李建庆,罗显刚等人基于不同土地利用情景利用SWAT模型对北江流域进行了水文模拟,结果表明:城市化以及退耕还林还草等土地利用变化都会对流域自然水文过程造成较大影响[5],綦昕瑶,刘贵花等利用IHACRES水文模型对东江流域径流变化的原因进行分析时发现:气候变化是引起东江流域径流变化的主要因素,人类活动是次要因素[6]。而在西江流域水文方面的研究由于研究资料的匮乏以及其他原因相对较少,廖卫红,雷晓辉等利用分布式水文模型EasyDHM对西江流域径流进行模拟,结果表明该模型在西江流域具有良好的适用性[7]; 曾凌,熊立华等利用DDRM(DEM-based distributed rainfall-runoff model)模型模拟了西江流域的土壤湿度的时空分布,结合ASCAT卫星遥感反演土壤湿度产品作为对比,结果表明DDRM模型和卫星遥感土壤湿度指数具有良好的时间一致性和空间一致性[8]; 赵胤懋,廖卫红等以西江为研究区利用地面雨量站的数据对比分析了CMORPH降水产品精度,分析结果表明:CMORPH降水产品具有较高的产品精度,能够替代地面雨量站的数据[9]; 吴志勇,林青霞等选取1951—2010年水文站径流数据,利用SRDI干旱指数研究西江流域水文干旱的空间变化特征,结果表明:西江流域轻旱发生范围有扩大趋势,极旱发生范围有减小趋势[10]。综合来看针对西江的水文研究大都集中在某一点或者某一产品的数据精度上,对于西江流域整体的系统性水文研究则相对较少。本研究利用遥感技术获取水文分析的基础数据、利用GIS空间分析技术对水文数据进行转换与处理,以西江作为目标研究流域,探究分布式水文模型VIC模型在西江流域的适用性情况,以期为西江流域水资源科学管理及洪水灾害防治提供一定的科学依据。
西江[11]地处北纬21°31'—26°49'、东经102°14'—114°48',流域横跨云南、广西以及贵州,集水面积约35万km2,占整个珠江流域总面积超过七成,年均径流量超过2 300亿m3,主要支流包括左江、右江、黔江、红水河、郁江、桂江以及柳江等。
西江流域在气候上的垂直差异和水平差异相当明显; 不同海拔,不同经纬度地区之间差异明显。流域气温差异明显,年平均气温常年保持在14~23℃。流域年均降水无论在空间上还是时间上分布都极不均匀。时间上:汛期(4—9月)水量占到全年总径流量的70%; 空间上则呈现出东南多西北少的趋势,年均降雨量在1 000~2 200 mm。
最初版本的VIC[12-14]模型又叫VIC-2L模型,其中的2L代表的是对土壤的分层,即将土壤分成了上下两层。随着越来越多的水文工作者对VIC模型研究水平不断提高,VIC版本随之不断更迭,将第二层土壤进行了分割,划出顶薄层,新增了对于土壤表层的动态描述,提高了土壤湿度的模拟精度,最终形成了如今被广大水文研究工作者广泛采用的VIC-3L模型。VIC-3L模型考虑了裸土以及不同植被覆盖类型对于降雨的动态响应。
VIC模型在产流方式上参考新安江模型同时考虑了两种产流方式(蓄满产流和超渗产流[15-17]),在模拟历时任意一个Δt内,该段时间内的降雨量会被VIC模型产流过程拆分为蓄满产流R1、超渗产流R2以及下渗到土壤的降雨ΔW3个部分。Wt是t时刻的土壤含水量,而不同性质土壤的空间分布的不均匀性导致降雨的入渗产流的方式也随着空间变化而变化,从蓄满产流变成超渗产流,网格内降雨在土壤中的入渗能力也随土壤性质的变化而变化,具体如下式:
f=fm[1-(1-C)1/B] (1)
式中:f表征不同土壤的降雨入渗能力在空间上的变化; fm代表流域内网格降雨量的最大入渗能力大小; C表征的是流域内降雨的入渗能力值小于或者等于f的网格面积占流域总网格面积的比例; B的意义与上式类似,是表征流域内土壤入渗能力参数。
P=R1(y)+R2(y)+ΔW(y) (2)
式中:P为在历时内的流域降雨量总和; R1为蓄满产流量; R2为超渗产流量; ΔW为入渗到土壤的雨量,且y=R1(y)+ΔW(y)。其中蓄满产流和超渗产流加上入渗到土壤的总水量ΔW可以用y式来进行表示,具体如下:
式中:fmm为流域土壤平均入渗能力,其他参数的意义与前面公式说明相同,经过上式的计算后就可以得到地下基流、地表径流以及三层土壤的含水量。
VIC模型的一大特点是通过全局参数文件的参数设置可以自由决定是否开启能量平衡模拟从而在研究流域内同时进行能量平衡模拟以及水量平衡模拟。此外VIC模型将流域内气候、土壤属性、地形地貌以及植被的综合作用以数学公式的方式进行集成处理,使得VIC模型具备扎实的数学意义和较高的物理基础。最后VIC模型源码开源,方便后续研究人员做模型的改进与发展。
VIC模型在应用过程中,模型的预处理是一个必不可少的环节,借助Arcmap平台,基于STRM 90 m高程数据对西江流域进行划分,以2 km×2 km模型分辨率在Arcmap软件中利用西江流边界将西江流域划分为723个网格单元。
本研究所使用的部分VIC模型输入数据来源如下:(1)地表覆盖数据来源于马大1 km分辨率全球地表覆盖数据集。该数据集的处理方法为:在流域所划分好的所有网格中,统计网格内所有的地表覆盖类型,同时将所有类型的地表覆盖属性信息赋值到对应网格上。(2)土壤质地数据采集自世界土壤数据库; 土壤质地数据集的处理原则与地表覆盖数据大体类似,区别之处在于网格内占比最大的土壤质地的理化性质数据代表整个网格的理化性质属性数据。(3)气象驱动数据来源中国地面累年日值数据集。该数据集包含了全国699个气象站点的70 a日值气象数据。本研究针对西江流域提取了流域周围共88个气象站点的气象数据进行研究,其数据处理原则为:根据流域周围共88个气象站点的气象要素数据,从中提取降水、最高气温、最低气温、风速4个子要素组成新的气象要素数据集,作为运行VIC模型的气象驱动数据。(4)用于VIC模型模拟结果的率定和验证所需要的实际径流量数据来自于武宣水文站1981—1989年逐日径流量数据。
根据已经收集的数据可知,西江流域土壤主要以沙壤土、黏土以及沙质黏土为主,其中又以沙壤土占据最主要部分。土地覆盖类型主要草地、耕地和林地3种基本覆盖类型
M-K[18-19](Mann-Kendell)趋势检验方法是分析水文气象要素在时间上变化特点常用的一种分析方法,该方法的优势在于用于检验的时间要素序列样本无须遵循某一特定的分布规律,在进行变化特点分析时极少会受到异常值的影响,而且计算过程简单且易于分析。M-K趋势检验方法的原理是假设有一个时间序列Xn,将每一个时间序列的子元素进行对比,根据其数值的大小来判断该时间序列是否具有某种变化趋势,具体见表1。
表1中Z值的正负表示该时间要素序列变化趋势为增加趋势还是减小趋势,当值大于零时为增加趋势,小于零表示该要素序列表现出减小趋势,Z的值的大小用来定量表征该时间序列变化趋势是否显著,参考他人研究[20]当Z>|±1.64|,则该要素序列具有某一显著的变化趋势,正值具有显著增加趋势,反之则有减小趋势。
VIC模型土壤参数众多,绝大部分的参数有其固定物理意义[21],通过不同土壤属性数据的获取,得到流域不同网格的土壤理化性质参数,而参与率定的土壤参数[22]共6个,其中控制研究流域地下基流量的土壤参数变量包括:Dmax,Ds,Ws,表征流域不同区域土壤蓄水能力大小的形状变量b以及由VIC模型所划分出来的土壤第二层厚度D2和第三层土壤厚度D3。本研究借助Matlab科学计算语言将待率定的VIC模型土壤参数设置成土壤参数变量,每次用不同的土壤参数变量的运行VIC模型,得到对应的模型模拟日值数据,将该次模拟数据与水文站实测径流数据进行对比,对比完成后,更改土壤参数变量的值,再次运行VIC模型,重复上述过程,直到所有的土壤变量组合方式都运行完毕。
根据现有的武宣水文站实测日值径流数据以及模型模拟日值径流数据通过评价指标纳什效率系数(Nash-Sutcliffe,NS)[23]、相关性系数(r)[24]以及相对偏差(BIAS)[25]对VIC模型的在对应土壤参数变量条件下的率定结果进行评估,评价指标具体数学公式见表2。
在上述3个公式中Qi,c为第i天的日尺度模拟径流量; Qi,o是第i天的日尺度实测径流量; n为总模拟历时,单位为天; Q^-o为实际日值观测径流量的平均值; Q^-m为模拟日值径流量的平均值。
由于西江流域面积广阔其数据量相应的也会更大,导致每次模拟过程所需时间较长(30 min以上),若将每种可能的组合方式均运行一次,在时间成本上是不可接受的,因此本研究参考Rosenbrock法[26]通过确定参数变化步长的方法来减少模型运行所需时间成本,即对VIC模型土壤参数变量的变化范围进行划分,每次只使用一个变化步长范围内的土壤参数变量,输入到VIC模型中,运行VIC模型,得到一组模拟径流数据,重复上述过程直到所有的土壤参数变量组合都运行完完毕,最后根据评价指标的最优值找出最优土壤参数变量组合。
本研究将整个VIC模型模拟过程划分为率定期(1981—1986年)和验证期(1987—1989年)两个阶段。土壤参数变量的物理意义及率定工作完成后得到的土壤参数变量最优值结果见表3。
以武宣站站率定期(1981—1986年)和验证期(1987—1989年)实测日值径流作为参考,VIC模型模拟结果见图1—2。
通过图1率定期与图2验证期实测与模拟日流量过程对比以及表4西江流域VIC模型的评价指标在率定期和验证期的结果可以看出在率定期:NS系数在日尺度上的值为0.71,在月尺度上的值为0.84。在相关性系数r的值上,日尺度的值为0.84,与NS系数的表现类似,月尺度的值较日尺度有所上升为0.95,表现优异。在相对偏差BIAS上,月尺度较日尺度也有所下降,分别为9%以及8%。在验证期内:3个评价指标总体上较率定期均有不同程度的下降,其中NS系数在日尺度和月尺度上的值分别为0.65,0.85。相关性系数r的变化情况与NS系数的变化情况类似,在月尺度上的表现要优于日尺度,分别为0.81,0.94,在相对偏差上,日尺度的值为13%而月尺度的值为9%较率定期内的表现均有不同程度的增大。
为了进一步分析VIC模型在西江流域年内丰水期(4—9月)和枯水期(10—12月、1—3月)的模拟效果是否存在差别,本研究将整个模拟历时进行拆分,将每年的丰水期以及枯水期分别提取出来单独组成一个径流时间序列,分析VIC模型在年内不同时期的模拟效果是否存在差别。
图3和图4为率定期丰水期以及枯水期的水文过程线,由图3可以看出在1981—1986年的丰水期中虽然水文过程曲线的变化趋势保持一致,但在多个年份出现了洪水峰值模拟不足。但VIC模型在枯水期(图4)的表现上,其模拟效果无论是从变化趋势上看,还是对于峰值的模拟,其效果都要好于丰水期的模拟效果,部分情况下实测径流量曲线甚至与模拟径流量曲线处于重合状态。
由表5中可知:日尺度下,率定期VIC模型在丰水期的NS的值为0.60,BIAS的值为13%,月尺度NS的值为0.63,BIAS的值为11%,而枯水期的NS的值为0.89,BIAS的值为4%。在月尺度下,枯水期的NS的值为0.91,BIAS的值略差于日尺度为6%。验证期的模拟结果与率定期类似,具体为在日尺度验证期的模拟结果中; 丰水期NS的值为0.55,略优于率定期日尺度NS的值,BIAS的值为17%。在月尺度上NS的值为0.58,BIAS的值为13%。在枯水期的模拟中,日尺度上NS的值为0.78,BIAS的值为6%。月尺度上NS的值为0.87,BIAS的值为8%。因此在不同时间尺度条件下,VIC模型在枯水期的模拟效果无论是在率定期还是验证期都要优于丰水期的模拟效果。
由率定期和验证期的水文过程线可以看出,VIC模型在年内不同时期的模拟精度上枯水期的模拟精度都要优于丰水期的模拟精度,但从总体上年际模拟效果上看,对比廖卫红等人利用EasyDHM模型对西江流域的径流模拟结果,虽然VIC模型的输入数据种类更多、数据要求也更高,但该模型也有着更好的模拟效果,同时也能够对流域热量的传输过程进行模拟,这是EasyDHM模型所不具备的。综合来说VIC模型在西江流域表现良好,可以为西江流域水资源合理分配、旱涝灾害预警与防治提供一定的科学依据。
利用经过率定与验证得出的最优土壤参数组合,对VIC模型进行更深层次的应用。基于已经建立好的西江流域VIC模型框架,单独制作了1971—2010年共四十年的气象数据作为VIC模型的气象驱动数据以驱动VIC模型运行,获取1971—2010年的VIC模型产流过程输出的西江流域多个水文要素数据。
VIC模型产流过程的输出文件参数众多,内容涵盖了经过率定过程得出的最优水量性质参数以及最优能量性质参数。本研究提取了水量性质参数中的P(多年平均降水)、R(地表径流)、E(蒸散发)、B(地下基流)4个水文要素,以3个月为一个时间步长,结合M-K趋势检验方法基于网格探究西江流域P,R,E,B这4个水文要素的近40 a的变化规律。
本文以基于网格的形式,对VIC模型所划分出的西江流域723个网格进行逐网格变化趋势分析,其原理是根据VIC模型产流过程所生成西江流域723个网格1971—2010年共40 a的水文要素,每一个网格单独进行M-K变化趋势检验分析,计算各个网格的Z值,在所有网格的完成M-K变化趋势分析后,将所有的Z值进行汇总,制作成ASCII格式,导入ArcMap中进行可视化处理。
研究表明:在这40 a中降雨在西江流域约一半的地区有增加趋势,而Z>1.64的网格共213个占西江流域总网格的29%,主要分布在西江流域的西部云南和贵州境内以及东部广西境内,降雨表现出下降趋势的网格主要集中在西江流域中部,但其Z均没有大于1.64,即增加趋势不明显。蒸散发具有增加趋势的网格集中在西江流域的中部贵州境内,而Z>1.64的网格较降水有一定的增加,共266个网格,占流域总网格数的37%; 与降雨类似,蒸散发也具有部分网格具有下降趋势,但其Z值均没有大于1.64。西江流域地表经流要素具有增加趋势的网格最多,其Z>1.64的网格共476个,在流域总网格数中占比超过65%,主要集中在广西境内流域部分。地下基流相较于其他3个水文要素变化最不明显,具有显著增加趋势的网格共10个,主要集中在广西境内,占比不到2%。
流域自然水循环过程的模拟需要综合考虑多方面的因素,尤其是大尺度的流域水文物理过程尤其复杂。本研究依托RS技术和GIS空间分析技术完成了分布式水文模型VIC模型的资料收集和数据预处理工作,在西江流域构建了VIC模型数据库,利用3个评价指标(NS,r,BIAS)评估VIC模型在西江流域径的适用性,结果表明:VIC模型在不同的时间尺度的径流模拟过程中均表现出了良好的适用性。文章最后对VIC模型的输出文件做了进一步的应用,从中提取了P,E,R,B共4个水文要素,利用M-K变化趋势分析方法,对这西江流域这4个水文要素在1971—2010年共四十年内的变化趋势进行了逐网格分析,分析结果表明:西江流域有超过65%的流域面积的地表径流R具有显著增加趋势,年均蒸散发E和年均降雨P具有显著增加趋势的面积占流域总面积约35%。地下基流B则相对稳定,具有显著变化趋势的面积占流域总面积不到2%。
尽管大尺度VIC模型在流域自然水文物理过程的模拟研究中具有众多优点,但是由于研究区数据资料的限制以及模型参数的不确定性等方面的影响,会对实际的模拟过程产生不良后果。如本研究在西江流域所采用的VIC模型分辨率为2 km,较大的模型分辨率会实测流域内很多的细节被忽略、概化,从而影响模型的模拟精度,因此在未来的研究中可以考虑采用更加精细的数据精度以提高模型的模拟精度从而更好的了解流域自然水循环过程的变化特点。