松花江是我国的七大河流之一，其流域面积为55.68万km²，约占东北地区总面积的44.8%，约占全国面积的5.8%.当前松花江流域水资源开发利用率为21.7%(根据水利部编制的《全国中型水库建设规划》报告)，远低于国际上通常认为水资源开发利用的生态警戒线30%~40%(李云生，2007)，水资源总体开发利用程度偏低.同时，近年来随着流域内经济的快速发展和城镇化进程的不断加快，松花江流域内水污染物排放量逐步增多，水环境状况逐渐恶化，松花江流域面临着越来越严重的水污染问题，其主要表现为入河污染量大，不仅涉及工业污水和生活污水排放，还有日趋突出的农业非点源污染，严重地制约着松花江流域内经济社会的可持续发展(李玮等，2010).针对松花江流域日益突出的水资源和水环境问题，开展对其水资源和水环境系统特征及时空分布规律的研究显得尤为重要，许多学者在这方面也开展了大量工作(Li et al., 2009; He et al., 2011;Wang et al., 2012a)，如针对松花江流域水环境质量的特征分析(许琳娟等，2012)、松花江流域工业污染物结构的特征研究(杨丽花等，2012)、松花江流域面源污染特征分析(烟贯发等，2013)、松花江流域水功能区限制纳污控制研究(周娜等，2014)等.已有的研究成果从水环境承载力、水环境污染控制等方面开展了理论、模型和方法的研究(於方等，2009；Wang et al., 2012b; Wang et al., 2012c;曾维华等，2012; Trainor et al., 2013)，为解决松花江流域水污染总量分配问题提供了参考依据.

如何在保证流域内社会经济快速发展的同时，既能高效合理地开发利用水资源，又能有效地控制水体污染并改善流域水质，实现经济社会与水资源、水环境的可持续发展，已成为松花江流域决策管理中一个亟待解决的重大问题.水资源的开发利用、水环境的污染与控制是一个复杂的、有机联系的系统整体，涉及到社会、资源、环境、法规、技术等诸多方面的因素，各个因素实时动态变化并相互作用反馈，进而形成了一个庞大复杂的流域动态循环系统.系统动力学(System Dynamic，SD)将定性分析与定量分析相结合，采用计算机仿真技术，分析系统内部结构与其动态行为特征，预测未来事物发展趋势，适合研究非线性、不确定性、动态性的复杂大系统问题(解海静等，2012;刘婧尧等，2014).因此，本文采用系统动力学方法对松花江流域水资源消耗与水污染物排放的动态变化特征及联合控制过程的动态模拟展开研究，目的是为了能够对松花江流域水资源与水环境的动态变化机理过程进行系统刻画，研究其未来可能的发展趋势和量化其动态过程的变化行为，提高对松花江流域的决策能力和管理水平.

松花江流域水资源-水污染联合控制动态模拟系统是一个多层次、多目标、多变量、动态性的综合系统.首先，要分析水资源消耗过程和水污染物排放过程中的动态变化特征并建立二者之间的有机联系；其次，从水污染物排放源头来划分水资源和水污染联合控制的子系统，筛选关键、典型的并能够对不同子系统中的变化状态进行综合刻画的指标变量；第三，选择合理的数学模型(或公式)来对不同变化状态特征进行量化；第四，采用系统论观点从水资源消耗和水污染物排放到水资源和水污染联合控制整体过程进行综合考虑，构建水资源消耗和水污染物排放之间的耦合关系，注重系统之间和系统内部各因素之间的动态变化特征和相互关系刻画；最后，应用系统动力学工具软件Vensim PLE，建立松花江流域水资源-水污染联合控制动态模拟模型，并在一定的假设条件下对不同的情景决策方案进行仿真，模拟各种情景决策方案条件下水资源消耗的动态变化对流域内经济社会的影响，以及通过水污染控制对水环境的改善效果，进而为未来松花江流域水资源优化调控和水污染减排制定出科学、合理的决策方案.具体流程见图 1.

图 1(Fig.1)

图 1 水资源-水污染控制联合控制系统流程 Fig.1 The system flow of combined control of water resources and water pollution

系统研究边界范围选取为松花江流域及其汇水区域.为了能够对松花江流域水资源消耗和水污染物排放的动态变化过程进行全方位、多层次的综合反映，揭示其对流域内经济社会发展的反馈作用，数据必须要能够支撑对模型的历史性验证.本文以2011年为基准年，时间跨度为2012—2020年，模拟时间间隔为1年.模拟历史数据资料主要来源于《中国统计年鉴—2012》、《黑龙江省统计年鉴—2012》、《吉林省统计年鉴—2012》、《内蒙古统计年鉴—2012》及《重点流域水污染防治规划(2011—2015)》报告.

从水污染物产生的源头考虑，将松花江流域水资源与水污染划分为4个子系统，分别为农业子系统、畜禽养殖子系统、工业子系统、生活(城镇和农村)子系统.其中，农业、畜禽养殖业、工业、生活(城镇生活和农村生活)子系统共同决定了松花江流域水资源消耗量与水污染物排放量，影响着松花江流域水环境系统，同时，水环境系统也对其他子系统有相互反馈作用.

松花江流域是国家粮食基地的重点建设地区，农业用水和总用水需求增长较快.长期以来为保障流域内农业生产，化肥、农药施用总量大，但有效利用率低，大型灌区农田退水污染问题突出，农业面源污染有加重趋势.农业子系统是为了对松花江流域的面源污染特征进行系统反映.根据灌溉面积的增长变化，合理控制农业施肥，提高灌溉用水量，降低农业废水排放量用总氮和总磷的排放量，有效控制农业面源污染.农业子系统选择了种植业灌溉面积为关键状态变量，构建了详细的农业子系统SD流图(图 2).

图 2(Fig.2)

图 2 农业子系统SD结构流图 Fig.2 Flow chart of agricultural subsystem SD model

松花江流域畜禽养殖近几年增长较快，根据污染源普查结果，畜禽养殖业也是松花江流域面源污染的主要来源，不同畜禽类型和畜禽数量的污染物排放量差异较大.畜禽养殖子系统的构建充分考虑了松花江流域畜禽养殖现状，选择了猪、肉牛、奶牛、肉鸡、蛋鸡和羊6类畜禽养殖数量作为主要状态变量，具体的畜禽养殖子系统SD流图见图 3.

图 3(Fig.3)

图 3 畜禽养殖子系统SD结构流图 Fig.3 Flow chart of livestock and poultry breeding subsystem SD model

松花江流域是我国重要的工业基地，随着工业经济的复苏和发展，造成了流域内水环境逐渐恶化，成为制约振兴东北老工业基地经济发展的瓶颈.工业子系统以工业总产值为状态变量，以工业增加值作为流速变量，通过工业增加值产值的取水量现状和历史变化趋势，设定不同阶段工业增加值取水量指标，推测目前或将来为实现某一工业增加值目标所需的工业用水量、新鲜水取水量.同时，结合工业用水重复利用率和工业废水排放处理情况，确定工业子系统最终的废水及COD、NH₄⁺-N污染物排放量.具体工业子系统SD流图详见图 4.

图 4(Fig.4)

图 4 工业子系统SD结构流图 Fig.4 Flow chart of industrial subsystem SD model

随着松花江流域内人口的不断增长，导致流域内生活用水量和生活污水排放量增加，从而加剧了流域内水资源紧张并增大水污染负荷.生活子系统主要考虑了松花江流域人口数量变化、人口结构、城镇化水平等因素对区域人口总数的影响，而人口数量的发展直接决定了松花江流域内生活源污染物的排放量.生活子系统选择了以总人口数量作为状态变量，人口增长率作为流率变量构建生活子系统，生活子系统SD结构流图详见图 5.

图 5(Fig.5)

图 5 生活子系统SD结构流图 Fig.5 Flow chart of life(urban and rural)subsystem SD model

以国家“十二五”流域污染物减排目标和《重点流域水污染防治规划》为依据，根据本文所建立的农业子系统SD模型、畜禽养殖子系统SD模型、工业子系统SD模型和生活子系统SD模型，分别通过模拟仿真计算农业、畜禽养殖业、工业和生活对水资源的消耗量、废水的排放量、水污染物(COD、NH₄⁺-N)的排放量;然后根据建立的松花江流域水资源-水污染联合控制数学模型(公式(1))，计算松花江流域整个区域内水资源消耗总量、废水排放总量、COD排放总量及NH₄⁺-N排放总量.

由于所建立的水资源-水污染联合控制动态模拟模型结构复杂且涉及到的参数众多，而模型参数设定的合理性直接决定了最终仿真结果的准确性和合理性.对于定量化参数变量的取值，可直接通过测量方法、公式计算和模型模拟来确定，而对于半定量化参数变量的取值，可借助统计和预测技术手段来参数估计，如趋势推算法、取平均值法、灰色系统方法及回归分析法等.此外，对于缺乏历史数据的参数变量的取值，主要根据专家实际经验做出合理评估确定.本研究主要参数取值见表 1.

表 1(Table 1)

表 1 主要表函数参数取值 Table 1 Parameters for the main table function

表 1 主要表函数参数取值 Table 1 Parameters for the main table function 主要参数 单位 取值 2011年(基准年) 2015年 2020年 注：规模化畜禽养殖比例为猪 ∶ 肉牛 ∶ 奶牛 ∶ 肉鸡 ∶ 蛋鸡 ∶ 羊，湿法比例为猪 ∶ 肉牛 ∶ 奶牛 ∶ 肉鸡 ∶ 蛋鸡 ∶ 羊. 单位面积化肥施用量 kg·hm-2 34.9 32.6 29.8 农业消耗系数 0.65% 0.60% 0.58% 单位面积灌溉用水量 m3·亩-1 420.0 392.6 364.3 人口自然增长率 0.10% 0.15% 0.20% 农村生活耗水系数 35% 30% 25% 沼气化率 35% 40% 45% 城镇居民耗水系数 30% 33% 35% 生活废水占管网废水总产生量比例 88% 90% 95% 工业增加值 亿元 6635.8 8795.4 11494.9 工业用水重复利用率 65% 75% 85% 规模化畜禽养殖比例 35.0%∶25.0%∶42.0%∶ 53.0%∶40.0%∶25.0% 42.5%∶35.0%∶49.0%∶ 60.0%∶50.0%∶30.0% 50.0%∶45.0%∶55.0%∶ 65.0%∶60.0%∶35.0% 湿法比例 61.5%∶59.5%∶ 59.5%∶ 39.0%∶ 39.0%∶20.0% 30.0%∶ 30.0%∶ 30.0%∶ 15.0%∶ 15.0%∶10.0% 20.0%∶20.0%∶20.0%∶ 10.0%∶ 10.0%∶5.0%

选取2009—2011年的总人口、耗水总量、工业GDP、废水排放量、COD排放量和NH₄⁺-N排放量的历史数据进行模拟验证，将得到的模拟结果与历史数据进行比较，以验证所建立的模型系统结构的合理性和模拟结果的准确性.从表 2检验结果可知，所选变量的模拟值与历史值的相对误差(|模拟值-历史值|/模拟值)均小于10%，说明所建立模型的仿真结果与历史检验值拟合程度较好，可用于对松花江流域水资源消耗与水污染排放联合控制的预测模拟.

表 2(Table 2)

表 2 模拟值与实际数据的对比 Table 2 Comparison of simulation and actual data

表 2 模拟值与实际数据的对比 Table 2 Comparison of simulation and actual data 年份 检验结果 总人口/ 万人 耗水总量/ 亿m3 工业GDP/ 万亿元 废水排放 总量/亿t COD排放 总量/万t NH4+-N排放 总量/万t 2009 历史值 6017 317.1 1.58 29.4 110.2 11.1 模拟值 6023 319.7 1.69 31.5 109.5 12.0 误差 0.09% 0.8% 6.5% 6.67% - 0.63% 7.5% 2010 历史值 6252 323.1 1.76 32.6 107.1 11.4 模拟值 6263 335.7 1.89 34.1 109.7 11.8 误差 0.06% 3.8% 6.9% 4.39% 2.37% 3.39% 2011 历史值 6583 362.9 1.94 36.7 332.2 20.8 模拟值 6593 377.2 2.11 38.5 341.6 22.6 误差 0.15% 3.8% 8.1% 4.67% 2.75% 7.9%

以当前松花江流域经济社会发展现状为基准，应用本文所提出的模型在不同情景和约束条件下进行仿真模拟.为便于对模型仿真模拟得到的结果进行比较分析，在此仅选出了具有代表性的4种情景决策方案.

按照当前松花江流域经济社会发展的速度，保持水资源消耗、人口增长速率及污染物排放量等指标值的现状，模拟松花江流域未来水资源消耗与水污染物排放的趋势.同时，该方案作为其他情景方案的比较基础.

开源、节流提高工业和农业节水能力，工业再生水回用率在原有基础上以每年2%的速度递增，万元工业产值耗水量在原有基础上以每年10%的速度递减；大力发展农业节水灌溉技术，农业产值耗水量在原有基础上以每年1%的速度递减.

提高废水的处理率，在现有40%的基础上逐年提高至80%；提高畜禽养殖粪便处理工程，降低畜禽养殖排放污染物入河量；提高化肥的利用率，减少化肥、农药施用量，降低农业面源污染.按照国家污染物总量削减要求，到2020年工业产排污系数在现状基础上削减60%，沼气入户率达到45%，农村污水分散处理率达到25%以上，其他参数取值不变.

综合方案二和方案三，在保持当前经济稳定增长的基础上，以开源、节流为主，兼顾治污、减排，实现水资源和水环境的可持续发展.该方案综合考虑了松花江流域水资源开发利用和水环境改善，是松花江流域社会经济和水环境系统协调发展的最优策略.

应用本文所建立的松花江流域水资源-水污染联合控制系统模型，分别对上述涉及的4种情景方进行了仿真模拟.松花江流域2012—2020年水资源消耗总量、废水排放总量、COD排放总量及NH₄⁺-N排放总量仿真模拟结果如图 6所示.

图 6(Fig.6)

图 6 耗水总量、废水排放总量、COD排放总量和NH4+-N排放总量模拟结果对比 Fig.6 Comparison of simulation results for the total volume of water consumption，total volume of waste water discharge，total emission COD and total emission NH4+-N

方案一是维持当前松花江流域经济社会发展的现状，来对松花江流域未来水资源消耗与水污染物排放进行模拟预测.由图 6可知，在保持松花江流域当前发展速度条件下，方案一消耗的水资源量最大，2012年的耗水总量是389亿m³，到2020年增加到475.3亿m³，增长了18.15%.尽管“十一五”和“十二五”期间对流域污染减排总量目标进行了控制，废水排放量、COD排放量和NH₄⁺-N排放量也逐渐下降，但按照当前松花江流域环境状况和水环境质量来看，方案一对于松花江流域未来污染控制和减排的压力很大，对流域内经济的发展有很大的制约，仍然不能够改变当前松花江流域水资源利用率低下和水环境恶化的现状.

方案二是通过开源、节流措施，有效提高松花江流域的水环境承载力.通过提高农业灌溉节水技术，提高农业灌溉水资源的利用率，降低农业灌溉废水的流失.同时改进工业中水利用的工艺技术，提高中水的回用率.根据仿真模拟结果(图 6)可知，通过实施方案二，2012年耗水量、废水排放量、COD排放量和NH₄⁺-N排放量分别为388.2亿m³、35.3亿t、280.5万t和18.6万t，与方案1相比分别下降了0.21%、3.5%、5.0%和2.5%；到2020年耗水量、废水排放量、COD排放量和NH₄⁺-N排放量分别为468.9亿m³、23.6亿t、233.6万t和16.7万t，与方案一相比分别下降了1.34%.2.59%、1.0%和8.5%.由此可见，开源、节流是提高松花江流域水环境承载力的关键所在，可以提高松花江流域水资源的利用率.

方案三是通过治污、减排措施，有效提高松花江流域水环境承载力和改善流域的水环境质量.根据仿真模拟结果(图 6)可知，通过实施方案三，2012年耗水量、废水排放量、COD排放量和NH₄⁺-N排放量分别为365.8亿m³、34.1亿t、274.9万t和18.2万t，与方案一相比分别下降了5.7%、6.7%、6.9%和4.9%；到2020年耗水量、废水排放量、COD排放量和NH₄⁺-N排放量分别为443.1亿m³、21.1亿t、224.2万t和16.3万t，与方案一相比分别下降了6.7%、12.6%、4.9%和8.8%.由此说明通过采取一系列污染物治理和减排措施，能够缓解松花江流域水环境承载力的压力和水功能区纳污的负荷.

方案四是方案二和方案三的综合，通过实施开源、节流、治污、减排的综合措施，松花江流域水环境承载力逐年上升，而且通过开展水环境污染治理和减排措施，流域内的水环境也得到有效的改善.根据仿真模拟结果(图 6)可知，通过实施方案四，2012年耗水量、废水排放量、COD排放量和NH₄⁺-N排放量分别为335.6亿m³、30.4亿t、258.4万t和17.0万t，与方案一相比分别下降了13.7%、17.0%、12.5%和11.1%；到2020年耗水量、废水排放量、COD排放量和NH₄⁺-N排放量分别为405.4亿m³、20.6亿t、204.1万t和13.9万t，与方案1相比分别下降了14.7%、15.0%、13.5%和24.2%.通过开源节流和治污减排的联合控制，松花江流域废水排放量和水污染物排放量远低于方案一、方案二和方案三中的排放量.因此，方案四是松花江流域经济社会与水环境系统协调发展的最佳策略.

根据当前松花江流域水环境承载力设计了上述的4种情景方案，对松花江流域水资源消耗量和水污染物排放量的情景进行了模拟，从大尺度上对松花江流域的水资源消耗和水污染物排放联合调控进行初步性探索研究，因此，本文研究侧重在水资源消耗总量和水污染物排放总量控制的系统刻画.由于模型当中涉及到众多半定量化和经验型的参数(如污染物去除率、废水处理率、污染物流失系数等)，且不能够对国家法规、经济发展速度及污染治理措施等因素进行完全模拟，使得在情景模拟过程中充满了诸多不确定性，最终模拟结果也存在一定的误差.因此，针对不同的对象在模拟中需要结合具体的环境来对参数进行率定和合理调整.

1)结合定性描述和定量分析的方法构建了松花江流域水资源消耗和水污染物排放二者联合控制的耦合关系，应用系统动力学方法建立了松花江流域水资源-水污染联合控制的动态模拟模型.通过制定松花江流域4种不同的发展策略方案，应用所建立的动态模拟模型进行了情景模拟并对仿真模拟结果综合分析.结果表明，本文所提出的联合控制动态模拟模型具有合理性和可行性，能够为松花江流域水资源管理和水污染控制提供决策支撑.

2)仿真模拟的结果表明，针对松花江流域单纯追求经济的高速发展而不注重对流域水资源合理开发和水环境保护，或是仅考虑水环境保护而减缓经济发展速度的策略都是不可取的，唯有兼顾流域经济社会发展和水环境保护的协调发展，才是实现松花江流域社会、经济和环境可持续发展的最佳策略.