中国水资源总量可观，但人均占有量低，部分地区存在水资源短缺、水质恶化的现象.水资源在中国天然分布不均，主要体现在时间上的夏季多、冬季少，空间上的南方多、北方少.随着国内外对水资源重视程度的提高，关于水资源的相关研究越来越广泛，研究内容主要包括水资源的利用效率、时空分布格局、可持续发展、水资源分配等^[1].

随着研究的深入，中国已经在法律层面重视大型流域的分段化管理，2017年修正的《中华人民共和国水污染防治法》提出：“省、市、县、乡建立河长制，分级分段组织领导本行政区域内江河、湖泊的水资源保护，水域岸线管理，水污染防治，水环境治理等工作.”实施河长制，可以因地制宜地对每段河流采用有针对性的管理措施，有利于河流的管控和可持续发展^[2].此外，以流域水资源分布情况进行生态功能划分也是目前研究的热点，早期的生态区划始于美国，Bailey等从生态系统的角度提出了美国的生态区划图^[3].国内生态区划分研究起步相对较晚，傅伯杰等进行了首个全国性的生态区划研究，将中国大陆地区分为3个生态大区、13个生态地区和57个生态区^[4].目前的研究划分结果多为三级功能区划，对于典型区域四级生态功能区划的研究还较少.

黑龙江省是全国粮食主产区，水资源储量丰富，地表径流发达，但对水资源的开发利用效率较低.哈尔滨作为黑龙江省的省会城市，是东北的老工业基地，其中阿城区更是省内重要的粮食生产基地，城镇、工业及农业用水的合理配置是当前急需解决的问题.哈尔滨市区段流经的主要河流为松花江，因此研究哈尔滨优控单元的水资源分布对城市发展、工业配置、农业种植有重要意义.《重点流域水污染防治规划》(2011—2015年)中将黑龙江省划分为16个控制单元，其中的松花江哈尔滨市市辖区控制单元(以下简称哈尔滨优控单元)是黑龙江省内的16个控制单元之一，也是黑龙江省控制区内4个优先控制单元之一.

在此背景下，选择哈尔滨优控单元作为研究区域，考虑采用模型对控制单元内的水资源时空分布情况进行研究.SWAT(soil and water assessment tool)模型是美国农业部开发的流域尺度分布式水文模型，目前在国内已经有广泛应用，包括对流域水量、悬移质、非点源等进行模拟预测分析等^[5-7].SWAT模型多用于大尺度流域的模拟预测，本文创新性地建立SWAT模型对特定研究区域哈尔滨优控单元范围内的水资源分布情况进行研究，并在目前三级生态区划的基础上以亚流域水资源敏感性作为标准继续划分了四级生态功能区，为实施河长制分段管理和规划哈尔滨优控单元的水资源提供一定的依据.

1 数据及研究方法 1.1 研究区概况

研究所选范围为哈尔滨优控单元，区域范围内的流域总面积为6 534.27 km²，优控单元所在位置如图 1，优控单元具体流域范围见图 2.哈尔滨优控单元位于松花江流域中段，朱顺屯及大顶子山两个国控监测断面之间，是中国重要的老工业基地和农牧场生产基地.流域内海拔高度73~98 m，属于中温带的大陆性季风气候.其春季多风干燥，夏季多雨温热；秋季较短，冬季寒冷漫长，寒潮频繁.区域多年平均气温较低，在3~5 ℃，年内温度变化较大，8月份日平均气温可达21~26 ℃，2月份可达-20 ℃以下.

研究区域主要的土地利用类型为农田和林地(见图 3)，其中水田主要为沿河分布，旱地主要集中在哈尔滨市区松花江干流以及支流阿什河下游周边，林地主要分布在支流阿什河的上游和松花江干流沿岸.

旱地是哈尔滨优控单元最主要的土地利用方式，占比为48.15%，其次是林地，占比为31.25%，其他各种土地利用总计占比20.6%.

1.2 数据来源

DEM来自中国科学院计算机网络信息中心国际科学数据镜像网站，土地利用图来自中国科学院资源环境科学数据中心，土壤类型图来自联合国粮农组织(FAO)，气象数据来自寒区旱区科学数据中心(SWAT模型中国大气同化驱动集CMADS V1.0)^[8]，水文数据来自哈尔滨市水文水资源勘测总站，点源数据来自环保部哈尔滨环境监测中心站，各数据精度及年限见表 1.

1.3 研究方法

SWAT模型包含水文过程模型、土壤侵蚀模型和污染负荷模型，本研究主要采用其中的水文过程模型，根据获得的基础数据建立SWAT模型，经过径流率定和验证，最终得到适用于研究区域的模型，对研究区域水资源分布情况进行分析^[9-10].

MK检验.Mann-Kendall(MK)法主要用于降水、径流等水文要素的变化趋势检验.计算公式如下^[11-12]：

$ \mathrm{MK}_{\mathrm{s}}=\left\{\begin{array}{cc}{\frac{S-1}{\sqrt{\operatorname{var}(S)}}, } & {S>0}, \\ {0, } & {S=0}, \\ {\frac{S+1}{\sqrt{\operatorname{var}(S)}}, } & {S<0}.\end{array}\right. $

(1)

式中：$S=\sum\limits_{i=1}^{n-1} \sum\limits_{k=i+1}^{n} \operatorname{sgn}\left(X_{k}-X_{i}\right), \quad \operatorname{var}(S)= $$ \frac{n(n-1)(2 n+5)}{18}, \operatorname{sgn}\left(X_{k}-X_{i}\right)=\left\{\begin{array}{cc}{1} & {X_{k}-X_{i}>0}, \\ {0} & {X_{k}-X_{i}=0}, \\ {-1} & {X_{k}-X_{i}<0}.\end{array}\right.$n为样本数，X_k，X_i为样本值.MK_s＞0，表示呈上升趋势，MK_s＜0，表示呈下降趋势.显著性水平α=0.05时，MK_s统计量临界值为1.96；显著性水平α=0.01时，MK_s统计量临界值为2.58.

2 结果与讨论 2.1 SWAT模型的建立

SWAT模型建立过程主要包括基础数据库的收集及本地化、亚流域及HRU的划分、敏感性分析及参数的调整、模型的率定及验证.针对哈尔滨优控单元区域进行研究，由于SWAT模型不能直接按照行政区域进行流域划分，选择大于哈尔滨优控单元的MASK进行圈定，然后根据SWAT模型自然划分的流域结果进行计算，最后选择哈尔滨优控单元边界内的区域进行研究.SWAT模型选择三家子断面作为上游入水口、大顶子山断面做为流域出水口，亚流域的划分积水面积阈值为13 070 hm^2[13-14]，最终将哈尔滨优控单元流域划分为25个亚流域^[15-16]，划分结果见图 4.

研究获得的数据时间尺度为2008—2014年，其中2008—2009年作为模型预热期，用来减少模型运行过程中初始化产生的问题，2010—2012年作为率定期，2013—2014年作为验证期.率定采用SWAT-CUP中的SUFI-2算法^[17-18]，得到的径流参数敏感性分析结果见表 2，径流模拟结果见图 5.

模型适用效果评价采用相对误差(E_R)、决定系数(R²)、Nash-Sutcliffe效率系数(E_NS)3个指标.采用月径流数据对模型进行率定和验证后，得到径流率定期E_R =1.42%，R²=0.77，E_NS =0.76；验证期E_R =-6.01%，R²=0.82，E_NS=0.74.根据得到的率定期及验证期的结果，对照Moriasi等给出的模型评价等级^[19]，0.65＜E_NS≤0.75评价等级属于“好”范围，E_R≤±10评价等级属于“非常好”范围，综合模拟结果评价等级属于“好”范围，因此，认为模型质量良好，本次研究建立的模型可以应用于松花江流域哈尔滨区段.

2.2 哈尔滨优控单元水资源时空变化特征

1) 时间变化特征.对哈尔滨优控单元2008—2014年的月径流数据进行分析，研究流域年内变化规律如图 6所示.可以看出，汛期(6—9月)径流总量达8 873.8 m³/s，占全年径流量的56%，其中8月份径流最多，平均径流量为2 805.4 m³/s，占全年径流量的17.7%.枯水期(11—2月)径流总量为2 542.6 m³/s，仅占全年径流量的16.1%.4月份开始，流域温度达到0 ℃以上，松花江开始融化，径流量增加，形成“桃花水”，因此，4月份流量比3月份增加幅度达49.3%.

对哈尔滨优控单元2008—2014年的年径流数据进行分析，研究流域年际变化规律见图 7.可以看出，年径流量范围为725~2 872 m³/s，多年平均径流量为1 320 m³/s.年径流量极大值2 872 m³/s出现在2013年，极小值725 m³/s出现在2008年，年径流量的极值比K_a为3.96，极值相差2 147 m³/s，可以看出，2008—2014年间哈尔滨优控单元降水量年际变化较大.

经计算，多年径流量数据的标准差为667 m³/s，年径流量的变差系数(C_v)为0.48，表明哈尔滨优控单元径流在年际间变化程度较大，其水源主要来自降水，地下水对河流的补给作用相对较小.

采用MK法检验哈尔滨优控单元径流的变化趋势及其显著性，选用的数据是2008—2014年间的月均径流量，样本总数为7，根据式(1)计算结果见表 3.可以看出，研究区在该时段内年均径流呈现不显著上升趋势；逐月分析的结果表明，除6月和8月外，研究时段内月均径流均呈现不显著上升趋势，6月和8月上升显著，上升速率分别为2.103和2.403 (m³/s)/10 a.

2) 空间变化特征.根据降水频率，2013年、2014年、2012年分别作为丰水年、平水年和枯水年的典型年，各年径流在亚流域的贡献分布见图 8.

在丰水年2013年，亚流域径流贡献率较大的为20，23，22，4，贡献率分别为10.64%，8.15%，7.94%和7.44%；在平水年2014年，亚流域径流贡献较大的为20，8，14，23，贡献率分别为8.2%，8.03%，8.02%和6.42%；在枯水年2012年，亚流域径流贡献较大的为20，23，22，14，贡献率分别为11.31%，8.64%，8.19%和5.64%.亚流域中的20，在丰水年、平水年和在枯水年其水文径流贡献量较其他区域都要大，属于丰水区域，参照土地利用图和坡度图，亚流域20土地类型主要是林地，区域坡度相对较大，汇水能力相对更强，因此，在各个年份均有较大贡献.亚流域8和14在平水年2014年内贡献比例较其他年份高，主要是由于2014年该区域降水量有所增加.亚流域3，5，9，12在丰水年、平水年和在枯水年径流占比贡献均较低，水资源相对缺乏.其中亚流域3为城镇，亚流域5主要为城镇及林地，亚流域9和12主要为旱地及水田.

可以看出，哈尔滨优控单元水资源分布具有时空分异性，合理地进行水资源调配利用可以在最大程度上整合区域水资源，充分发挥其作用.水资源与生态环境密切相关，基于一定的生态功能区划进行水资源调控，可以更科学地划分水资源功能区.

2.3 哈尔滨优控单元生态功能区划分

2015年发布的《全国生态功能区划》(修编版)^[20]，按照生态区、生态亚区、生态功能区依次进行了三级划分，哈尔滨优控单元按照现行三级生态功能划分可以分为阿什河上游水源涵养与土壤保持生态功能区(I-04-02-05)、哈尔滨市区城市与城郊农业生态功能区(I-05-02-01)、拉林—阿什河流域农业与土壤保持生态功能区(I-05-02-02)，具体见图 9.

哈尔滨优控单元中主要包括哈尔滨市区及周边农、林地区，属于水资源依赖性强的区域，目前该区域按照三级生态功能分区仅包含3个类别，三级划分的结果对其社会发展的规划、农业林业管理而言过于粗放，难以有效指导区域内的水资源综合利用，且哈尔滨优控单元水资源具有明显的异质性，根据三级生态功能划分结果，不能体现其水资源分布特点，基于此，本研究将以上3个功能区按照水资源敏感性特征，以亚流域作为水文生态功能的均质性单元，继续划分四级生态功能区^[21].水资源敏感性特征按照亚流域水资源占整个流域的百分比划分为3个级别，占比在6.42%~11.3%的为高水资源区、占比在2.3%~6.04%的为中等水资源区、占比在0.59%~1.94%的为低水资源区，将各区域按照水资源禀赋的高中低程度继续划分四级生态功能区^[22]，结果见图 10.

I-04-02-05-01为阿什河上游水源涵养保护区，该区域是哈尔滨优控单元内水资源最丰富的地区，目前主要土地利用类型为林地和旱地，可以考虑将该区域的旱田等量置换为其他区域的水田，充分利用区域水资源进行自然灌溉，同时兼顾水源涵养和农业种植需求；I-04-02-05-02为阿什河上游土壤保持生态功能区，该区域目前主要土地利用类型为林地和旱地，但森林覆盖率低，对土壤的侵蚀较敏感，因此，考虑在此区域进行适当的退耕还林，并在重点地区加密种植森林植被带，以达到防治水土流失和水土保持的目的.

I-05-02-02-01为拉林河—阿什河流域水田种植区，该区域水资源相对丰富，目前主要土地利用方式为旱地和水田，可以在此区域进行种植业调整，扩大水田面积，逐步减少旱地面积，充分利用自然水资源灌溉；I-05-02-02-02为拉林河—阿什河流域旱地种植区，该区域水资源相对较低，目前主要土地利用方式为水田和旱地，该区域更适合种植旱地，可以与I-05-02-02-01区域结合进行合理的调配，扩大区域内的旱地面积，逐步减少水田面积；I-05-02-02-03为拉林河—阿什河流域土壤保持生态功能区，该区域水资源低，目前主要土地利用方式为旱地和水田，建议将此区域进行适当的退耕还林，并在重点地区加密种植森林植被带，以达到防治水土流失和水土保持的目的.

I-05-02-01-01为哈尔滨市区城市与城郊核心区，该区域包括哈尔滨市城镇及周边部分区域，该区域应逐步降低农业耕作占用面积，逐步完善城市周边的防护林；I-05-02-01-02为哈尔滨市区城市与城郊农业种植区，该区域主要土地利用方式为旱地，可以在维持现有农业种植的基础上，适当增加林地植被，作为林地与耕地的交汇区；I-05-02-01-03为哈尔滨市区城市与城郊水源涵养林保护区，该区域紧邻哈尔滨市区，水资源较丰富，现有土地利用中已经有部分林地，建议增大林地密度，在松花江干流两侧建立水源涵养林带.

2.4 规划管理建议

研究范围内农田面积较大，根据实际情况，北方地区农业大面积实行旱改水，水田所占比例逐步增加，农业用水的需求量已经出现短缺，本研究根据水资源分布特点划定了四级生态功能分区，在水资源较丰富的区域多种植水田，等量替换其他区域的旱地，水资源相对缺乏的地区种植玉米、大豆等需水量少的经济作物，并在适宜的区域实行退耕还林.管理部门可以结合本次研究的结果，在区域内实行水资源统一调度与管理，充分利用哈尔滨优控单元内的水资源.

3 结论

1) 建立了适用于哈尔滨优控单元的SWAT模型，径流率定期E_R=1.42%，R²=0.77，E_NS=0.76；验证期E_R=-6.01%，R²=0.82，E_NS =0.74，模拟结果评价等级属于“好”范围.因此，认为模型质量良好，建立的模型可以应用于松花江流域哈尔滨区段.

2) 哈尔滨优控单元2008—2014年间年内径流在汛期(6—9月)最大，占全年径流量的56%，8月份径流最多，平均径流量2 805.4 m³/s，占全年径流量的17.7%，枯水期(11—2月)仅占全年径流量的16.1%.4月份流量与3月份对比增加幅度较大，达49.3%，径流量年际变化范围为725~2 872 m³/s，多年平均径流量为1 320 m³/s；2008—2014年间哈尔滨优控单元降水量年际变化较大，MK检验结果表明，年均径流呈现不显著上升趋势，逐月分析的结果表明，除6月和8月外，研究时段内月均径流均呈现不显著上升趋势，6月份和8月份上升显著，上升速率分别为2.103和2.403 (m³/s)/10 a.

3) 哈尔滨优控单元2008—2014年径流贡献量在空间上分布不均匀，丰/平/枯三期各亚流域水资源分布情况不同，水资源贡献较大的区域主要集中在亚流域20和23，其中亚流域20在丰水年、平水年和在枯水年径流贡献量上均为最大，亚流域9和12水资源相对缺乏.根据水资源空间分布情况，按照水资源敏感性特征，在三级生态功能分区的基础上，继续进行了四级分区，将哈尔滨优控单元一共划分了8个四级生态功能分区，水资源丰富的地区主要用于水源涵养、种植水田，水资源较少的地区用作旱地农业区，水资源少的地区用作土壤保持区.本次研究结果可为实施河长制分段管理和规划哈尔滨优控单元的水资源提供一定的依据.