2018, Vol. 30 
(2: 江西师范大学地理与环境学院, 南昌 330022)
(3: 江西省赣州市气象局, 赣州 341000)
(2: School of Geography and Environment, Jiangxi Normal University, Nanchang 330022, P. R. China)
(3: Ganzhou Weather Bureau of Jiangxi Province, Ganzhou 341000, P. R. China)
农林系统是非点源污染形成和发展的重要生态系统，非点源污染涉及到的气象、土壤、水文、植物生长、农事活动和管理等多种过程，都与农林生态系统密不可分.运用分布式机理性非点源模型进行时间和空间尺度上的模拟，是定量研究和评估农林系统对非点源污染影响的有效方法.在众多的流域分布式非点源模型中，SWAT(Soil and Water Assessment Tool)模型是目前具代表性、使用广泛、应用前景广阔的分布式机理性模型[1]，是流域尺度、时间连续、基于过程的模型.它不仅可以在水文响应单元的空间尺度上模拟地表径流、入渗、侧流、地下水流、回流、融雪径流、土壤温湿度、蒸散发、产沙输沙、作物生长、营养盐输移、流域水质等多种自然过程，还可以模拟包括耕作、灌溉、施肥、收割、用水调度等农业管理措施对这些过程的影响[2].该模型自1990s由USDA-ARS开发以来，在国内外已广泛应用于水质和水量的模拟评估、非点源污染负荷估算及形成机制探讨、情景分析与预测、环境变化及农业管理措施对水文水质的影响、气候变化对区域水循环和作物生长的影响等多方面[3]；应用的流域尺度从几百平方千米到几十万平方千米不等，甚至大到国家级或大陆级空间尺度，如Jayakrishnan等[4]将SWAT模型用于宏观分析和评估美洲大陆或美国全国的水质、水量等水资源管理措施的效应；应用的时间尺度从某次降水过程的数小时、数天到几十年不等[5].
SWAT模型虽然在国内外得到了广泛的应用[2]，但模型中单一的植物生长模式难于表达变化密度、多林种混交的森林植被状态以及不同的种植制度与农耕方式，因而这些模型的应用难于考虑植被景观、种植制度和农耕方式对非点源污染的影响.针对这些问题，国内外学者对该模型进行了不同程度的修正.如Watson等[6-7]引入了森林生长的机理模型3-PG来代替SWAT的植物生长模型；Kiniry等[8]通过引进并修正ALMANAC植物模型以改善SWAT模型在模拟森林生态系统有关水和生物地球化学循环方面不足的方法.国内学者在这方面，近年来也作了大量探索，代俊峰等[9]针对中国南方丘陵水稻灌区的水文特点，改进了SWAT的灌溉水运动模块、稻田水分循环模块、稻田水量平衡各要素和产量模拟的计算方法，增加了渠系渗漏模拟模块及其对地下水的补给作用、塘堰的灌溉模块等.谢先红等[10]在稻田蒸发蒸腾、控制灌溉排水、塘堰实时灌溉等方面对SWAT模型进行了改进，可以模拟“间歇”、“淹灌”和“薄浅湿晒”等灌溉模式的水稻灌区水分运动过程.仕玉治等[11]针对灌溉资料相对缺乏流域，以SWAT模型为基础，修改了其水稻田的自动灌溉制度，并构建了多水源综合灌溉模块.郑捷等[12]针对平原型灌区人工-自然复合的水文循环特点，考虑平原灌区灌溉渠道、排水沟和河道等人工干扰，在沟渠河网的提取方法、子流域与水文响应单元的划分以及作物耗水量计算模块等方面对SWAT模型进行了改进.
在非点源模型的研究与应用过程中，针对亚热带季风湿润区的森林植被景观和我国南方水稻种植区种植制度与农耕方式的特征，开展非点源污染建模及模拟方法的研究，并应用遥感技术获取植被混杂比例、复种和套种等相关参数，对于应用现有的非点源模型开展非点源污染的研究具有实际意义.因此，本文以农林系统的非点源污染模拟为目标，通过研究建立变化密度及多种类混杂的森林生长模型，修正SWAT采用平均森林植被密度和单一植物生长模式估算生物累积量的问题，并建立与之相适应的森林优势组份丰度遥感反演模型、叶面积指数和消光系数遥感反演模型以获取森林生长模型的相关参数.同时，根据间作套种下的辐射能利用Keating方程，引入间作套种指数变量，修正模型原有的单一生物量日积累模型，并以亚热带季风湿润区红壤背景下的鄱阳湖流域为试验区，以野外实测数据为基础，对SWAT修正模型的有效性进行探讨.
1 SWAT模型植物生长模式的应用局限与修正 1.1 SWAT模型植物生长模式及应用中存在的问题分析 1.1.1 SWAT模型原有植物生长模式的缺陷SWAT植物生长模型源自简化的EPIC(Environmental Policy-Integrated Climate)模型，它利用单一的植物生长模型来模拟所有植物的生长，SWAT模型主要适用于高度农业区的非点源模拟[3].农林系统的差异导致了多个研究者对SWAT模型进行了修正[6-7].然而，Watson[6]和MacDonald等[7]引进的3-PG模型需要大量森林生长数据的支持，增加了模型的复杂性，使模型的应用受到一定的局限；而ALMANAC模型原本来源于EPIC模型，为多作物通用模型，其作物生长模型不再单一，但其生物量累积过程并没有考虑到间作套种这种复合群体的增益效能[8].在最新的SWAT 2012版本中，虽然作了适当的修正，但其生物量累积过程中的两个重要参数：消光系数和叶面积指数仍然还是简化形式，所有植物的消光系数都被缺省地定为-0.65，而叶面积指数则基于平均的植物密度进行估算，且由统一的理想叶面积发育模型来控制.这种单一的植物生长模型难于描述森林系统和农业系统中多林种混交、间作套种的植被景观和种植方式.而农林系统是流域生态系统的重要子系统，也是氮(N)和磷(P)循环的重要界面，模型的这种简化局限了该模型在农林系统中的一些应用.
1.1.2 模拟间作套种农耕方式的局限SWAT模型主要用于估算大流域的非点源负荷，预测、评估不同土壤、土地利用/覆盖和农业管理对水、泥沙沉积、农业化学物质的长期影响.但由于其植物生长模型的单一，难于表达复杂的种植方式，且在农事管理上，虽然考虑了多种农事活动和农事管理，但这种考虑主要基于美洲和欧洲等国家粗放型的农耕方式[6]，对我国常见的间作套种等种植制度和农耕方式基本上无法处理，而这些耕作方式极大地改变了农作物的光能利用效率[13-14]，改变了作物对N和P的吸收，影响了营养物质的循环.
1.2 SWAT模型原有的植物生长模式SWAT的植物生长模型为单一的植物生长模型，它主要由两部分组成，一是生物量累积模型，另一个是理想叶面积发育模型.具体内容如下：
(1) 日生物量最大增量由植物叶面积截取的光合有效辐射来估算：
| $ \Delta bio = D \cdot RUE \cdot Reg \cdot {H_{{\rm{phosyn}}}} $ | (1) |
式中，Δbio为逐日的总生物量潜在增量; RUE为植物的辐射利用率; D为空间单元中的平均密度; Reg为温度、水分和氮磷养分的综合胁迫因子(其值在0~1之间)，如温度的胁迫以植物(作物)生长的下限温度、最适温度和上限温度等来表达，以植物(作物)生长的S型曲线来描述，植物的三基点温度及其他胁迫因子数据都存储在SWAT的植物生长数据库中.
而Hphosyn为给定时间内植物叶面积截取的光合有效辐射，该变量由下式确定：
| $ {H_{{\rm{phosyn}}}} = 0.5{H_{{\rm{day}}}}\left( {1-\exp \left( {{k_1} \cdot LAI} \right)} \right) $ | (2) |
式中，Hday为在给定时间内入射的总辐射，0.5Hday为光合有效辐射，kl是消光系数(缺省值为-0.65)，LAI是叶面积指数，由理想叶面积发育模型估算.
(2) 理想叶面积发育模型.植物生长初期，在叶面积达到最大叶面积之前，冠层高度和叶面积发育由理想叶面积发育模型(曲线)所控制：
| $ f{r_{LAI{\rm{mx}}}} = \frac{{f{r_{PHU}}}}{{f{r_{PHU}} + \exp \left( {{l_1}-{l_2} \cdot f{r_{PHU}}} \right)}} $ | (3) |
式中，frLAImx为对应于所给植物潜热单位比(frPHU)的植物最大叶面积指数比, frPHU为该植物在生长季节内给定时间的潜热单位比, l1和l2为形状系数. frPHU的计算公式为：
| $ f{r_{PHU}} = \frac{{\sum\limits_{i = 1}^d {HU} }}{{PHU}} $ | (4) |
式中，HU为第i天的累积热量单位，PHU为该植物的总热量单位.
形状系数通过2个已知点(frLAI, 1，frPHU, 1)和(frLAI, 2，frPHU, 2)解方程(3)来计算，其公式如下：
| $ {l_1} = \ln \left[{\frac{{f{r_{PHU, 1}}}}{{f{r_{LAI, 1}}}}-f{r_{PHU, 1}}} \right] + {l_2} \cdot f{r_{PHU, 1}} $ | (5) |
| $ {l_2} = \frac{{\left\{ {\ln \left[{\frac{{f{r_{PHU, 1}}}}{{f{r_{LAI, 1}}}}-f{r_{PHU, 1}}} \right]- \ln \left[{\frac{{f{r_{PHU, 2}}}}{{f{r_{LAI, 2}}}}-f{r_{PHU, 2}}} \right]} \right\}}}{{f{r_{PHU, 2}}-f{r_{PHU, 1}}}} $ | (6) |
某天的冠层高度由下式计算：
| $ {h_{\rm{c}}} = {h_{{\rm{c, mx}}}} \cdot \sqrt {f{r_{LAI{\rm{mx}}}}} $ | (7) |
式中，hc为给定某天的冠层高度(m)，hc, mx为该植物的最大冠层高度.
1.3 基于多植物生长模式的SWAT模型的修正 1.3.1 森林生长模型的调整与修正原始的SWAT植物生长模型用于森林生长，存在如下问题：一是由于叶面积季节变化不大，常绿树木类叶面积如采用理想叶面积发育模型来估算，会存在一定的误差；二是所有的森林植被冠层的消光系数都设置为-0.65；三是森林密度采用平均密度来表示.
基于上述考虑，对SWAT模型有关森林植被的生长模型，提出如下修正方法：
(1) 对于落叶类树木，仍然采用原来的植物生长模型；(2)对于常绿树木类，植物生长模型的理想叶面积发育模型由叶面积指数与遥感植被指数统计模型来代替；(3)对于消光系数，通过试验建立研究区优势植物类型遥感消光系数模型来代替统一的消光系数-0.65.研究区优势植物类型遥感消光系数模型由实测的消光系数与遥感植被指数的统计关系来确定；(4)落叶常绿针阔混交林的植物生长模型由下式确定：
| $ {H_{{\rm{phosyn}}}} = \sum\limits_{i = 1}^n {{P_i} \cdot {H_{{\rm{phosyn, }}i}}} $ | (8) |
式中，Pi为i类优势树种在混交林中的比例，采用流域森林优势植被组份丰度来表示，采用遥感方法来获取. Hphosyn, i为i类优势树种给定时间内植物叶面积截取的光合有效辐射，n为混交林中优势树种的数量. (5)森林密度问题的修正，SWAT植物生长的密度体现在日生物量最大增量公式中(公式(1))，它对于农作物种植来说，问题不大，但对于森林系统来说，问题比较大.在SWAT模型原有的输入变量中，森林的密度主要体现在土地利用/覆盖中的林地类的2级分类，即有林地、灌木林地等.因此，在式(1)中引入流域森林植被覆盖度来代替其中的平均密度D，流域森林植被覆盖度采用遥感方法提取.
1.3.2 间作套种方式下农作物生长模型的建立间作套种是一种在时间和空间上实现种植集约化的种植方式.根据Keating和Carberry[14-15]的方程可以得出间作套种的高秆作物在给定时间内植物叶面积截取的光合有效辐射为：
| $ \begin{array}{l} {H_{{\rm{phosyn, H}}}} = 0.5{H_{{\rm{day}}}}\\ \left\{ {\left[{1-\exp \left( {{k_{1, {\rm{H}}}} \cdot LA{I_{{\rm{H, 1}}}}} \right)} \right] + \frac{{{k_{1, {\rm{H}}}} \cdot LA{I_{{\rm{H, 2}}}}}}{{{k_{1, {\rm{H}}}} \cdot LA{I_{{\rm{H, 2}}}} + {k_{1, {\rm{L}}}} \cdot LA{I_{\rm{L}}}}}\left[{1-\exp \left( {{k_{1, {\rm{H}}}} \cdot LA{I_{{\rm{H, 2}}}} + {k_{1, {\rm{L}}}} \cdot LA{I_{\rm{L}}}} \right)} \right]} \right\} \end{array} $ | (9) |
而矮秆作物在给定时间内植物叶面积截取的光合有效辐射为：
| $ \begin{array}{l} {H_{{\rm{phosyn, L}}}} = 0.5{H_{{\rm{day}}}}\\ \left\{ {\frac{{{k_{1, {\rm{L}}}} \cdot LA{I_{\rm{L}}}}}{{{k_{1, {\rm{H}}}} \cdot LA{I_{{\rm{H, 2}}}} + {k_{1, {\rm{L}}}} \cdot LA{I_{\rm{L}}}}}\left[{1-\exp \left( {{k_{1, {\rm{H}}}} \cdot LA{I_{{\rm{H, 2}}}} + {k_{1, {\rm{L}}}} \cdot LA{I_{\rm{L}}}} \right)} \right]} \right\} \end{array} $ | (10) |
式中, LAIH, 1、LAIH, 2和LAIL分别为高秆作物上层和下层的叶面积指数和矮秆作物的叶面积指数，kl, H和kl, L分别为高矮搭配作物类型的消光系数，上层与下层以矮秆作物的冠层高度为分界点. LAIH, 1、LAIH, 2的计算由下式表示：
| $ \left\{ \begin{array}{l} LA{I_{{\rm{H, 1}}}} = \frac{{{h_{\rm{H}}}-{h_{\rm{L}}}}}{{{h_{\rm{H}}}}}LA{I_{\rm{H}}}\\ LA{I_{{\rm{H, 2}}}} = \frac{{{h_{\rm{L}}}}}{{{h_{\rm{H}}}}}LA{I_{\rm{H}}} \end{array} \right. $ | (11) |
式中, hH和hL分别为高秆和矮秆作物的冠层高度，由作物性状和生长习性数据库提供.因此，间作套种产生的日最大生物量累积为：
| $ \Delta bio = RUE\left( {{I_{\rm{H}}} \cdot {H_{{\rm{phosyn, H}}}} + {I_{\rm{L}}} \cdot {H_{{\rm{phosyn, L}}}}} \right) $ | (12) |
式中, IH、IL分别为高秆和矮秆作物的间作套种指数.所谓间作套种指数，是指在单位面积中某种作物在间作套种中所占的比重.不同作物的间作套种指数将采用遥感技术通过建立间作套种指数遥感模型来获取.
1.4 SWAT模型源代码的修正及编译SWAT模型由Fortran语言编制而成，其执行框架由主程序及大量的子程序构成，共有307个程序文件.
针对公式(8)和公式(12)，构造了3个新变量P(森林优势植被组份丰度)、IH(高杆作物的间作套种指数)和IL(矮秆作物的间作套种指数)，由于叶面积指数和消光系数两个变量在原始的SWAT模型中就有，因此无需增设.由于SWAT模型的叶面积指数(LAI)是在水文响应单元(HRU)中的一个变量，即在HRU中植被的LAI是均一的，因此，公式(8)和公式(12)的P、IH和IL均在水文响应单元级别上提取其数值，并增加了相应的数据读取子程序ReadNewVar.f.该子程序完成修正后模型所需要的变量包括LAI、植被冠层消光系数(K，即(2)式的kl，与后面的EC意义相同)等原有变量和P、IH和IL等新增变量.
由于SWAT模型对于模拟间作套种存在一定的局限(见2.1.3)，对SWAT模型农事管理模块的修正，主要增加了一个判断某HRU内是否存在间作套种的逻辑变量InterPlant，并增设了新的间作套种规模模块InterPlantMgt.f，该模块参照原有的管理模块，增加了允许多种作物并存管理的功能.
SWAT模型的修正与编译是采用Intel Visual Fortran Compiler 11.1.067编译器在Visual Studio 2005平台上修正及编译完成.
2 试验区概况及模型的数据准备 2.1 试验区概况选取梅江流域(26°17′34″~27°08′54″N，115°44′28″~116°16′06″E)为研究区，该流域集水面积约3304 km2.境内以中低山为主，高程分布在170~1440 m，平均高程为379.1 m(图 1).贯穿全流域的梅江发源于江西省宁都县肖田乡北缘的王陂障，梅江是赣江的水源区，也是鄱阳湖的主要源头之一，属赣江的上支.流域内气候温和湿润，季风明显，雨量充沛，年平均气温17.3℃，极端最高气温37.9℃，极端最低气温-6.2℃，年降雨量1706 mm，年平均相对湿度80 %以上.土壤多为变质岩、花岗岩发育而成的山地红壤和黄红壤，质地为壤土和轻砂壤，土层深厚肥沃.流域内地带性植被为中亚热带常绿阔叶林，森林覆盖率达70 %，天然林面积占森林面积的87 %.流域内的植被群系主要为针叶林、常绿落叶阔叶混交林、常绿阔叶林、落叶阔叶林、竹林、矮(曲)林、灌丛、草丛.由于研究区处于亚热带湿润季风区，其森林分布特征以混交林为主，即使是人工林，也难见单一的林木.天然植被群系中，常绿阔叶林和针叶林优势明显，其中针叶林的优势树种是马尾松(P. massoniana Lamb.)和杉木(Cunninghamia lanceolata(Lamb.)Hook.)，主要分布在中低丘陵岗地；常绿阔叶林有17种树种，但优势树种为丝栗栲(C. fargesii Fr.)和苦槠栲(C. sclerophylla (Lindl.) Schott.)，主要分布在人类活动较少且海拔高于600~700 m的高山上[16-17]，低丘陵以樟(Cinnamomum camphora)和木荷(Schima superba)为主，主要分布在居民地附近；落叶阔叶林的垂直地带性分布表现出低矮岗地以杨树(Populus)、桉树(Eucalyptus spp.)、苦楝(Melia azedarach)、乌桕(Sapium sebiferum)等为主，中高山地以枫(Liquidambar formosana Hance)、栎(Q. aliena var. acu-teserrata)为主；竹林以毛竹(P. pubescens Mazel)为主，但总量不多，主要分布在低矮岗地.此外，流域内还有零星分布的果茶园林.
|
图 1 研究区的地形(a)、主要河流及雨量站分布(b) Fig.1 Distribution of terrain(a), domain rivers and precipitation stations in the study area(b) |
流域内主要水系为梅江(干流)、琳池河(支流)、黄陂河(支流)、会同河(支流)、竹坑河(支流)、固厚河(支流)6条主干河流(图 1).其中5条支流的水量占整个梅江干流水量的92.3 %.
按2010年人口统计资料，三种产业结构比为24.7 :39.7 :35.6.农业以水稻种植为主，兼有脐橙、莲子等作物种植.行政区涉及广昌县、乐安县、宁都县、石城县、兴国县、宜黄县共6个县，其中宁都县占研究区面积的96.4 %，其余县仅占3.6 %.
2.2 模型的数据输入及其处理SWAT模拟过程中涉及的数据包括2部分：1)SWAT模型需要的基本数据，主要有：地形DEM、河网、土地利用与覆盖、土壤类型及其理化属性、植物与植物生长数据、肥料数据、水文与水质数据、气象数据和农业养殖数据与人口排放数据[18]，其中水文与水质数据主要用于模型的建立、模型参数的率定和验证；2)修正后SWAT模型所需的数据，包括植被覆盖度、森林优势植被组份丰度、森林优势植被冠层叶面积指数及消光系数、农耕方式间作套种指数及复种指数等森林植被与间作套种参数.
2.2.1 森林植被与间作套种参数的获取方法(1) 森林优势组份丰度的遥感估算方法流域森林优势组份丰度是反映流域中某森林组份的单位面积比重，是修正后SWAT模型的一个重要变量.本文基于线性混合模型方法通过遥感技术来获取研究区森林优势组份丰度的空间分布数据.线性混合模型假定混合像元的光谱是该瞬时视场内各类地物光谱的线性组合[19-20].根据研究区的实际地表覆盖，将研究区地表覆盖分成阔叶林、针叶林、草地与稀疏灌丛、裸地、水体、耕地、居民地共7种地表覆盖组份，其中森林植被包括阔叶林、针叶林和草地与稀疏灌丛3种基本组份.作为研究成果的一部分，地表覆盖组份丰度的遥感提取方法见文献[21-23]，提取的地表覆盖丰度是以遥感影像的像元为空间单元，一个像元中的7种地物覆盖丰度累加值为1.研究区阔叶林、针叶林和草地与稀疏灌丛3种森林植被组份丰度见图 2.
|
图 2 研究区3种主要森林植被组份丰度 Fig.2 Three domain forest vegetation component abundance of the study area |
(2) 森林优势植被冠层叶面积指数及消光系数的遥感反演叶面积指数和消光系数是表征植被群体冠层结构及光能利用的地球表层下垫面参量，也是修正后SWAT模型的重要参数.我国南方地处中亚热带湿润地区，雨量丰沛，植被多呈混交状态，且疏密程度不均，少见单一大面积的均匀分布.因此，以目前的遥感技术及遥感信息源难于获取单一树种或单一植被的叶面积指数和消光系数.以遥感技术反演的叶面积指数或消光系数只能是一定空间范围内森林优势植被类型的叶面积指数或消光系数.为更好地匹配叶面积指数和光合有效辐射(用于估算消光系数)的实测数据，反映植被混交和疏密不均的状态，本文以美国Landsat ETM作为遥感信息源，通过Gram-Schmidt(GS)图像融合方法，将ETM融合成空间分辨率为15 m的多光谱数据.在实测森林优势植被冠层叶面积指数和光合有效辐射的基础上，利用植被指数法经验公式法反演流域的叶面积指数，并根据Beer-Lambert定律，建立了流域优势植被冠层消光系数的反演模型，具体方法参见文献[24]. 图 3给出了利用该反演模型估算的研究区森林优势植被冠层叶面积指数和消光系数.
|
图 3 研究区森林优势植被冠层的叶面积指数和消光系数 Fig.3 The LAI and EC of forest domain species vegetation canopy of the study area |
(3) 农耕方式间作套种指数及复种指数的遥感信息提取农业系统中的农耕方式在原有的SWAT模型中已有多种方式体现，集中在农事管理模块中，如播种、施肥/杀虫、深耕、收割等，但间作套种的农耕方式在原有的SWAT模型中体现的不够.为此本文修正了SWAT模型作物生长模型的日最大生物量累积方程(公式(12)).式中的IH、IL将通过建立间作套种指数遥感模型获取.
复种问题在SWAT模型中已有所体现，可以通过农事管理模块来实现有关复种的相关参数的设置及数据库建立，复种指数的遥感反演也有比较成熟的方法[25-27].本文所涉及的试验区复种指数的遥感反演参见文献[28].
2.2.2 模型的其他输入数据(1) DEM数据和河网数据. DEM和河网的GIS数据来源于1 :5万的地形图，主要用于SWAT模型的流域离散化过程.
(2) 土地利用数据.土地利用/覆盖数据是SWAT模型的主要输入变量，同时也是人口及牲畜养殖排放估算、森林优势组份丰度模型建立、农作物复种与间作套种指数模型建立和叶面积指数与消光系数模型建立的重要参考数据.为了获得研究区的土地利用/覆盖数据，采用了高级陆地成像仪(Advanced Land Imager，ALI)数据作为遥感数据源(来源于中国科学院计算机网络信息中心国际科学数据服务平台(http://datamirror.csdb.cn)的Level L1Gst数据).通过多光谱数据与全色波段的Gram-Schmidt(GS)数据融合，生成了10 m分辨率的多光谱假彩色合成图像，通过目视解译方法对研究区土地利用/覆盖进行了解译，并对解译结果进行了多次调研并验证.形成的最终类别为：农村居住用地、城镇及交通建设用地、旱地、未利用地、水田、低矮植被(包括草地、灌丛)、疏林地、有林地、河流、湖泊水库坑塘10个类别(图 4).
|
图 4 土地利用与覆盖 Fig.4 Land useand cover |
(3) 土壤数据.土壤数据主要包括土壤类型及其理化属性，来源于江西省宁都县2009年完成的《宁都县耕地地力评价》成果和1986年完成的土壤普查成果[29]，其中土壤类型空间分布地图比例尺为1 :5万(图 5). SWAT的土壤用户数据库(usersoil)包括17个物理属性，其中SOL_CLAY(layer #)、SOL_SILT(layer #)、SOL_SAND(layer #)、SOL_ROCK(layer #)和SOL_Z(layer #)按参考文献[29]进行转换，SOL_BD(layer #)、SOL_AWC(layer #)、SOL_K(layer #)、SOL_CBN(layer #)、USLE_K等变量则按照参考文献[30]进行转换并计算；而化学属性存放在SWAT的土壤输入文件(.chm)中，共有7个属性.
|
图 5 土壤类型空间分布 Fig.5 Spatial distribution of soil type |
(4) 植物/作物生长数据.植物/作物的生长数据主要包括植物生长发育所需的有效积温、发育的最低温度、最高温度和最适温度等农业气象参数，以及最大冠层高度、最大潜力叶面积指数、消光系数、根深、收获指数等品种性状参数.部分植物参数取自SWAT模型自带的参数，如松树、杉树、水稻、大豆等，研究区特有的作物参数来自《宁都县耕地地力评价》所收集的宁都县作物品种特征数据.
(5) 肥料数据与各地施肥数据.肥料数据来自于《宁都县耕地地力评价》所收集的肥料信息，包括所施用肥料的名称、含氮量、含磷量、含钾量等，利用这些数据，建立了SWAT模型运行所需肥料数据库；研究区不同区域农业施肥数据来源相同，包括了24个乡镇共21062户农户施肥品种和每亩施肥量、亩产等内容.
(6) 水文和水质数据.水文和水质数据主要用于SWAT模型的参数率定和模型验证，包括河流流量数据和河流主要营养盐浓度数据.为了获得研究试验区水文水质的实测数据，对研究区空间离散化了87个子流域，设置了9个断面的测点(图 6)，这些测点大约控制了70个子流域，占研究区子流域总数的80.5 %；对流量与水质进行了6次测定，其中水质采样要素包括TN、TP、氨氮、硝酸盐氮、亚硝酸盐氮.
|
图 6 水质采样点的分布 Fig.6 Distribution of the water quality sampling sites |
(7) 气象数据.采用了研究区周边9个气象站的数据，包括宁都县、宜黄县、瑞金市、石城县、于都县、南丰县、广昌县、泰宁县(福建)、长汀县(福建)1957-2011年共54 a的气温(最高和最低)、湿度、风速和降水数据，利用这些数据，按照SWAT气象用户数据库的要求，整理成SWAT的气象用户数据库[31].此外，在流域内还采用了33个雨量站2005-2011年共7 a的实测降水数据，作为模型的输入，用于模型的建立、参数率定和验证.
(8) 农业养殖数据与人口排放数据.研究区不仅包含林业、农业和少量的工业，同时还有一定规模的养殖业，为了能较好地利用SWAT模型模拟流域内的非点源污染，农业养殖和人口排放不能不考虑.本文对农业养殖和人口排放采用了GIS的空间分析方法，建立了在中小尺度流域农业养殖和人口排放空间化方法，具体可参见文献[32-33].
3 模型的调试及参数率定采用了汇流累积区面积为3000 hm2的阈值进行了空间离散化，生成了87个子流域，并采用土地利用/覆盖、土壤、坡度比例阈值为20 :10 :20的水文响应单元(HRU)划分阈值，得到700个HRU单元.
选取了2005-2011年作为模型调试运行、参数率定、验证的时间，其中2005-2007年为调试运行时间，2008-2010年为模型参数率定时间，2011年为模型验证时间.
模型参数率定先对水文参数进行率定，率定完水文参数以后，再对营养盐负荷参数进行率定，这期间不再对水文参数进行率定.使用的模型均为原始模型.
修正模型的模拟参数采用原始模型率定的参数结果，除了模型进行了修正及输入参数进行了相应调整外，其余变量及参数均控制在与原始模型的变量与参数一致，并使用修正模型的模拟结果进行不同研究目标的分析.
3.1 流量的参数率定与敏感性分析参数率定采用SWAT-CUP软件进行，该软件采用t检验的方法，对变量进行敏感性评价，一般认为P < 0.05为敏感性变量.流量的参数敏感性试验包括CN2、ALPHA_BF、GW_DELAY、GWQMN、GW_REVAP、ESCO、SOL_AWC等13个参数，各参数名称、物理意义和率定结果见表 1，其中参数名称中的扩展名表示所在数据库名称.
|
表 1 流量的率定参数及率定结果 Tab.1 Calibration parameters and results of the flow |
图 7给出了这13个参数的率定敏感性程度，从P值可以看出，CH_K2、ALPHA_BF和GW_REVAP 3个参数比较敏感，其中以CH_K2最敏感[34].
|
图 7 流量参数率定的敏感程度t值和P值 Fig.7 The flow parameters calibration of sensitive degree of the t and P |
氮营养盐负荷的参数率定选择了NPERCO、CMN、ERORGN、SPCON、BIOMIX、SHALLST_N等8个参数，其中NPERCO、CMN 2个参数对氮营养盐的排放具有一定的敏感性(图 8a).磷营养盐负荷的参数率定选择了PSP、RSDCO、BIOMIX、PHOSKD和PPERCO等10个参数，其中PSP和SPEXP 2个参数反应敏感，以PSP参数最为敏感(图 8b).各参数的名称、物理意义及率定结果见表 2.
|
图 8 氮(a)和磷(b)相关参数率定的敏感程度t值和P值 Fig.8 The sensitive degree of t and P of relevant parameters of nitrogen and phosphorus |
|
|
表 2 氮和磷相关参数的率定结果 Tab.2 The calibration results of relevant parameters of nitrogen and phosphorus |
在氮和磷参数率定基础上，对上述反应敏感的参数，再用SWAT-CUP进行微调率定，得到最优的参数，利用这些率定好的参数进行模拟.在实测氮和磷负荷数据基础上对模拟结果进行验证，并分析模拟的有效性程度.
4 模型的验证与有效性分析 4.1 原始模型的验证与有效性分析在参数率定结果的基础上，用修正模型对2011年进行模拟.利用6次9个采样点的水质采样数据和同期的宁都县水文站的实测流量，与模型同期同断面输出的流量及营养盐负荷进行有效性分析.采用相关系数(R2)和有效性系数Nash-Sutcliff (NS)系数来衡量模型的有效性.其中NS系数的计算公式如下：
| $ NS = 1-\frac{{\sum\limits_{i = 1}^n {{{\left( {{\xi _{{\rm{o}}i}}-{\xi _{{\rm{m}}i}}} \right)}^2}} }}{{\sum\limits_{i = 1}^n {{{\left( {{\xi _{{\rm{o}}i}}-\overline {{\xi _{\rm{o}}}} } \right)}^2}} }} $ | (13) |
式中，ξo、ξo和ξm分别为观测值、观测值平均值及模拟值，n为统计的样本数. NS值在-1~1之间变化，其值越接近1，说明模型的有效性越好.
原始模型模拟的流量与实测流量及流域平均降雨量的对照图可以看出，模拟与实测的流量基本趋势一致，但模拟的峰值与实测的峰值存在较大的出入，主要表现在模拟的极大峰值偏小、中小峰值偏大(图 9).与流域平均降雨量比较可看出，流量模拟值与流域降雨量的变化趋势有比较好的一致性，模型的地表产流调蓄模拟方面不及实际情况.在研究区不同部位有多个不同规模的水库，由于没有获取到相关的水库数据，因而模型没有考虑水库对河道的调蓄作用.模型结果中基流和降雨量小的峰值模拟相对较差可能与此有关.
|
图 9 原始模型模拟的流量与实测流量及流域平均降雨量对照 Fig.9 Comparison of the original model simulation and measured flow and the average rainfall of basin |
原始模型模拟流量的有效性分析可以看出，模拟值与实测值的相关系数为0.753，达到极显著水平.此外模型的NS系数为0.686，说明模型存在较好的有效性(表 3).
|
|
表 3 原始的模型模拟流量的有效性分析 Tab.3 The effectiveness analysis of the original model of the flow simulation |
模型在模拟TP和TN营养盐负荷方面，时间流与空间流趋势都基本一致(图 10). 表 4给出了原始模型模拟营养盐方面的有效性分析，从营养盐负荷的平均值、最大值和最小值来看，模拟值总体比实测值偏高.另外，从R2值和NS值来看，模拟效果略差于流量的模拟，且TN的模拟效果比TP的模拟效果差.
|
图 10 原始模型模拟的TP、TN与实测的TP、TN浓度对照 Fig.10 Comparison of the original model simulation and measured TP, TN concentrations |
|
|
表 4 原始模型模拟营养盐方面的有效性分析 Tab.4 The effectiveness analysis of the original model of the nutrient simulation |
修正模型模拟的流量与实测流量及流域平均降雨量的对照示意图可以看出，修正模型在流量方面的模拟与原始模型基本一致(图 11)，但对比图 7，在流量的峰值模拟方面，修正模型有较大的改善，说明修正模型能够较好地反映地表径流方面的实际情况.
|
图 11 修正模型模拟的流量与实测流量及流域平均降雨量对照 Fig.11 Comparison of the modified model simulation and measured flow and the average rainfall of basins |
从模拟的流量平均值、最大值、最小值、流量总量来看，修正模型的模拟结果更接近实测值(表 5).相关系数由原始模型的0.753增加到修正模型的0.863，而NS系数则由0.686增加到0.740，分别增加了14.6 %和7.8 %.说明修正模型的流量模拟效果比原始模型的流量模拟效果要好.
|
|
表 5 修正模型的流量模拟有效性分析 Tab.5 The effectiveness analysis of the modified model of the nutrient loading simulation |
在营养盐负荷模拟方面，将9个水质采样点的6次(分别为2011年2月11日、4月2日、5月8日、7月11日、9月10日和11月21日)实测TP和TN浓度与模拟的TP和TN浓度见图 12和图 13.修正模型在模拟营养盐负荷方面的有效性分析可以看出，修正模型在模拟TP时效果好于模拟TN，R2与NS系数分别为0.82、0.78和0.71、0.69(表 6).对比原始模型，修正模型在模拟TP和TN方面的R2和NS系数，分别增加了3.8 %、4.4 %和6.4 %、6.1 %.
|
图 12 修正模型在模拟TP浓度的时空流与实测的时空流对照 Fig.12 Comparison of the temporal-spatial flow between the modified model simulation TP concentration and the measurement |
|
图 13 修正模型在模拟TN浓度的时空流与实测的时空流对照 Fig.13 Comparison of the temporal-spatial flow between the modified model simulation TN concentration and the measurement |
|
|
表 6 修正模型在模拟营养盐负荷方面的有效性分析 Tab.6 The effectiveness analysis of the modified model of the nutrient loading |
从模型对流量和营养盐的模拟结果来看，修正后的SWAT模型在模拟峰值流量和TN、TP时，效果优于原始的模型.说明相较于原始SWAT模型采用平均森林植被密度和单一的植物生长模式估算生物累积量，修正的SWAT模型采用变化密度、多种类和多种类混杂的森林生长模型，更能反映森林生长的真实情况，因而在进行森林植被景观对非点源污染影响的模拟中，可以有效利用植被覆盖度、森林组份丰度和叶面积指数等变量从不同角度来描述森林植被景观的不同状态.
5 结语以农林系统的非点源污染模拟为目标，通过研究并建立变化密度及多种类混杂的森林生长模型，修正了SWAT模型采用平均森林植被密度和单一植物生长模式估算生物累积量的问题.同时，根据间作套种下的辐射能利用Keating方程，引入间作套种指数变量，修正SWAT模型原有的单一生物量日积累模型.以亚热带季风湿润区红壤背景下的鄱阳湖流域为试验区，以野外实测数据为基础，探讨修正SWAT模型的有效性.结果表明：
1) 在流量的峰值模拟方面，修正模型有较大的改善，说明修正模型能够较好地反映地表径流方面的实际情况；从模拟的流量平均值、最大值、最小值、流量总量来看，修正模型的模拟结果更接近实测值.相关系数由原始模型的0.753增加到修正模型的0.863，而NS系数则由0.686增加到0.740，分别增加了14.6 %和7.8 %.说明修正模型的流量模拟效果比原始模型的流量模拟效果要好.
2) 在模拟营养盐负荷方面，修正模型在模拟TP浓度时效果好于模拟TN，相关系数与NS系数分别为0.82、0.78和0.71、0.69.对比原始模型，修正模型在模拟TP和TN浓度的相关系数和NS系数，分别增加了3.8 %、4.4 %和6.4 %和6.1 %.
3) 修正的SWAT模型采用植被覆盖度、针叶林、阔叶林及稀疏灌丛等不同森林组份丰度和植被叶面积指数等变量，比原始SWAT模型更真实和更有效地反映了森林植被景观的宏观与微观特征，因而修正的SWAT模型在模拟流量和营养盐方面比原始模型有更好的有效性.
| [1] |
Krysanova V, Arnold JG. Advances in ecohydrological modelling with SWAT-A review. Hydrological Sciences Journal-des Sciences Hydrologiques, 2008, 53(5): 939-947. DOI:10.1623/hysj.53.5.939 |
| [2] |
Lai GY, Wu DY, Zhong YX et al. Progress in development and applications of SWAT model. Journal of Hohai University:Natural Sciences, 2012, 40(3): 243-251. [赖格英, 吴敦银, 钟业喜等. SWAT模型的开发与应用进展. 河海大学学报:自然科学版, 2012, 40(3): 243-251.] |
| [3] |
Gassman PW, Reyes MR, Green CH et al. The soil and water assessment tool:Historical development, applications and future research directions. Transactions of the ASABE, 2007, 50(4): 1211-1250. DOI:10.13031/2013.23637 |
| [4] |
Jayakrishnan R, Srinivasan R, Santhi C et al. Advances in the application of the SWAT model for water resources management. Hydrol Process, 2005(19): 749-762. |
| [5] |
Muttiah RS, Wurbs RA. Modeling the impacts of climate change on water supply reliabilities. Water Resources Assoc, 2002, 27(3): 407-419. |
| [6] |
Watson B, Coops N, Selvalingam S. Integration of 3-PG into SWAT to simulate growth of evergreen forests. Hydrol Earth Syst Sci, 2005, 19(3): 827-838. |
| [7] |
MacDonald DJ, Kiniry J, Arnold R et al. Developing parameters to simulate trees with SWAT. 3rd International SWAT Conference, 2005.
|
| [8] |
Kiniry JR, MacDonald JD, Kemanian AR. Plant growth, and simulation for landscape-scale hydroplogical modeling. Hydrological Sciences Journal, 2008, 53(5): 1030-1042. DOI:10.1623/hysj.53.5.1030 |
| [9] |
Dai JF, Cui YL. Distributed hydrological model for irrigation area based on SWATⅠ. Principle and method. Journal of Hydraulic Engineering, 2009, 40(2): 145-152. [代俊峰, 崔远来. 基于SWAT的灌区分布式水文模型——Ⅰ.模型构建的原理与方法. 水利学报, 2009, 40(2): 145-152.] |
| [10] |
Xie XH, Cui YL. Analysis of scaling variation of irrigation water use efficiency under typical irrigation regimes. Engineering Journal of Wuhan University, 2009, 42(5): 653-660. [谢先红, 崔远来. 典型灌溉模式下灌溉水利用效率尺度变化模拟. 武汉大学学报:工学版, 2009, 42(5): 653-660.] |
| [11] |
Shi YZ, Zhang C, Zhou HC et al. Development and application of SWAT model to paddy district in watershed scale. Water Resources and Power, 2010, 28(7): 18-22. [仕玉治, 张弛, 周惠成等. SWAT模型在水稻灌区的改进及应用研究. 水电能源科学, 2010, 28(7): 18-22.] |
| [12] |
Zheng J, Li GY, Han ZZ et al. Application of modified SWAT model in plain irrigation district. Journal of Hydraulic Engineering, 2011, 42(1): 88-97. [郑捷, 李光永, 韩振中等. 改进的SWAT模型在平原灌区的应用. 水利学报, 2011, 42(1): 88-97.] |
| [13] |
Tsubo M, Walker SR, Ogindo HO. A simulation model of cereal-legume intercropping systems for semi-arid regions I. Model development. Field Crops Research, 2005, 93(1): 10-22. DOI:10.1016/j.fcr.2004.09.002 |
| [14] |
Gao Y, Duan AW, Liu ZG. Effect of monoculture and intercropping on radiation use efficiency and yield of maize and soybean. Chinese Journal of Eco-Agricultuer, 2009, 17(1): 7-12. [高阳, 段爱旺, 刘祖贵等. 单作和间作对玉米和大豆群体辐射利用率及产量的影响. 中国生态农业学报, 2009, 17(1): 7-12.] |
| [15] |
Keating BA, Carberry PS. Resource capture and use in intercropping:solar radiation. Field Crops Res, 1993, 34: 273-301. DOI:10.1016/0378-4290(93)90118-7 |
| [16] |
Zhang HX, Fang HY, Tu XB et al. Investigation and assessment of biodiversity. Nanchang: Jiangxi Science and Technology Press, 2010. [张海星, 方红亚, 涂晓斌等. 江西生物多样性调查与评估. 南昌: 江西科学技术出版社, 2010.]
|
| [17] |
Administrative Office of Ganzhou in Jiangxi Province Forestry Bureau of Reclamation ed. Gannan trees. Ganzhou: Administrative Office of Ganzhou in Jiangxi Province Forestry Bureau of Reclamation, 1981. [江西省赣州地区行署林垦局. 赣南树木. 赣州: 江西省赣州地区行署林垦局, 1981. ]
|
| [18] |
Neitsch SL, Arnold JG, Kiniry JR et al. Soil and Water Assessment Tool Input/Output File Documentation Version 2005, 2004.
|
| [19] |
Hill JH. Monitoring 20 years of increased grazing impact on the Greek Island of Crete with earth observation satellites. Journal of Arid Environments, 1998, 39: 165-178. DOI:10.1006/jare.1998.0392 |
| [20] |
Adams JB. Spectral mixture modeling:A new analysis of rock and soil types at the Viking Lander I Site1. Geophys Res, 1986, 91: 8098-8112. DOI:10.1029/JB091iB08p08098 |
| [21] |
Chen XZ, Lai GY, Pan RX et al. Extract forest information using the mixed end-members from images and field survey data. Lanzhou, Gansu, China: 2012 International Symposium on Geomatics for Integrated Water Resources Management, 2012: 2134-2139. http://ieeexplore.ieee.org/document/6349580/
|
| [22] |
Chen XZ, Lai GY. Abundance extraction of end-members of forest based on linear mixed model-A case study of Meijiang Basin. Nanjing, China: The 2nd International Conference on Remote Sensing, Environment and Transportation Engineering (RSETE 2012), 2012: 203-206. http://ieeexplore.ieee.org/xpls/icp.jsp?arnumber=6260377
|
| [23] |
Chen XZ, Lai GY, Pan RX et al. Vegetation abundance remote sensing information extraction based on linear spectral mixed model in the Hilly Area. Jiangxi Science, 2012, 30(4): 473-479. [陈绪志, 赖格英, 潘瑞鑫等. 基于线性混合模型的丘陵山区植被丰度遥感信息提取研究——以江西梅江流域为例. 江西科学, 2012, 30(4): 473-479.] |
| [24] |
Lai GY, Zeng XG, Liu Y et al. Retrieving leaf area index and extinction on coefficient of dominated vegetation canopy cover in Meijiang Watershed of China using image-fusion and Landsat ETM Data. Remote Sensing Technology And Application, 2013, 28(4): 697-706. [赖格英, 曾祥贵, 刘影等. 基于ETM和图像融合的优势植被冠层叶面积指数和消光系数的遥感反演. 遥感技术与应用, 2013, 28(4): 697-706. DOI:10.11873/j.issn.1004-0323.2013.4.697] |
| [25] |
Peng DL, Huang JF, Jin HM. The Monitoring for sequential cropping index of arable land in Zhejiang Province using MODIS-NDVI. Scientia Agricultura Sinica, 2006, 39(7): 1352-1357. [彭代亮, 黄敬峰, 金辉民. 基于MODIS-NDVI的浙江省耕地复种指数监测. 中国农业科学, 2006, 39(7): 1352-1357.] |
| [26] |
Yan HM, Liu JY, Cao MK. Remotely sensed multiple cropping index variations in China during 1981-2000. Acta Geographica Sinica, 2005, 50(6): 559-566. [闫慧敏, 刘纪远, 曹明奎等. 近20年中国耕地复种指数的时空变化. 地理学报, 2005, 50(6): 559-566.] |
| [27] |
Yan HM, Cao MK, Liu JY et al. Characterizing spatial patterns of multiple cropping system in China from multi-temporal remote sensing images. Transactions of the CSAE, 2005, 21(4): 85-90. [闫慧敏, 曹明奎, 刘纪远等. 基于多时相遥感信息的中国农业种植制度空间格局研究. 农业工程学报, 2005, 21(4): 85-90.] |
| [28] |
Pan RX, Lai GY, Li XJ et al. A methodology for multiple cropping index extraction of Mei Jiang Basin based on MODIS/NDVI. Jiangxi Science, 2012, 30(6): 764-768. [潘瑞鑫, 赖格英, 李秀娟等. 基于MODIS-NDVI的梅江流域耕地复种指数提取研究. 江西科学, 2012, 30(6): 764-768.] |
| [29] |
Soil Survey Office Nanjing County in Jiangxi Province, Nanjing County soil. Nanchang:Jiangxi People's Publishing House, 1988. Nanchang: Jiangxi People's Publishing House,, 1988. [江西省宁都县土壤普查办公室. 宁都县土壤. 南昌: 江西人民出版社, 1988.]
|
| [30] |
Wei HB, Zhang ZP, Yang JP. Establishing method for soil database of SWAT model. Water Resources and Hydropaver Engineering, 2007, 38(6): 15-18. [魏怀斌, 张占庞, 杨金鹏. SWAT模型土壤数据库建立方法. 水利水电技术, 2007, 38(6): 15-18.] |
| [31] |
Pang JP, Xu ZX, Liu CM. Weather generator and database in the SWAT Model. Journal of China Hydrology, 2007, 27(5): 25-30. [庞靖鹏, 徐宗学, 刘昌明. SWAT模型中天气发生器与数据库构建及其验证. 水文, 2007, 27(5): 25-30.] |
| [32] |
Yi FZ, Lai GY, Zhang L et al. Estimation and spatial distribution of the potential load of Three-Yellow Chicken Feces in Meijiang Watershed in China. Lanzhou, Gansu, China: 2012 International Symposium on Geomatics for Integrated Water Resources Management, 2012.
|
| [33] |
Zeng XG, Lai GY, Yi FZ. A study population data spatialization of small river basin based on GIS:Take Meijiang River Basin as the example. Geography and Geo-Information Science, 2013, 29(6): 40-44. [曾祥贵, 赖格英, 易发钊等. 基于GIS的小流域人口数据空间化研究——以梅江流域为例. 地理与地理信息科学, 2013, 29(6): 40-44.] |
| [34] |
SWAT-CUP4: SWAT Calibration and Uncertainty Programs-A User Manual. Eawag: Swiss Federal Institute of Aquatic Science and Technology, 2009.
|