水资源作为生产生活的重要资源，与人类生存环境和社会发展密切相关。然而随着我国社会经济的快速发展，水资源短缺和水污染严重已成为中国水资源面临的主要问题。中国水资源总量丰富，但人均占有量仅为世界平均水平的四分之一。据《2013年中国环境公报》数据显示，全国59.6%的地下水监测点水质属较差或极差级别，地表水中河流Ⅳ类及以上水质达到28.3%，水环境压力较大。传统研究大多把水资源数量与水污染单独研究，虽有Azeved^[1]、王西琴^[2]、夏星辉^{[3, 4]}等学者将二者相结合进行水资源综合评价，但依旧无法定量分析二者之间的关系，在此背景下灰水足迹的概念应运而生。

Hoekstra和Chapagain提出“灰水足迹”的概念，经水足迹网络的灰水足迹小组进一步完善，将灰水足迹定义为以水环境质量标准为基准，将污染负荷稀释至高于特定环境水质标准所需淡水的体积^{[5, 6]}。灰水足迹将水污染从水量角度进行评价，为水污染与水资源数量的综合研究提供了新的方法。国内外对水足迹的研究已经取得较大的研究成果，但单独针对灰水足迹的研究，至今尚未全面展开。目前国外研究主要集中在特定产品灰水足迹核算评价和区域整体灰水足迹核算评价两方面。其中，特定产品灰水足迹核算评价主要集中在：(1)农作物生产灰水足迹核算评价^{[7, 8, 9, 10, 11]}，如：Chapagain等^[7]从生产和消费领域研究了水稻的绿水足迹、蓝水足迹、灰水足迹。(2)畜牧产品灰水足迹核算评价^{[12, 13]}，如Mekonnen等^[12]从全球范围评价了农场动物的水足迹(其中包括灰水足迹)。(3)工业产品灰水足迹核算评价，如Ercin^[14]等探索了一种含糖碳酸饮料的水足迹(其中包括灰水足迹)。区域整体灰水足迹核算评价方面^{[15, 16]}如Mekonnen等^[15]从生产与消费两个方面计算了1995—2006年全球大部分国家的水足迹(其中包括灰水足迹)。同时，国内的研究也取得了探索性进展。国内灰水足迹研究主要聚集在农作物灰水足迹的核算与评价方面^{[17, 18, 19, 20, 21]}，如盖力强等^[17]、秦丽杰等^[21]基于化肥污染对区域小麦、玉米等粮食生产的水足迹(其中包括灰水足迹)进行了计算和评价。另外，王来力等^[22]从初始灰水足迹、工序灰水足迹等方面核算了纺织产品的灰水足迹。孙才志等^[23]与曾昭等^[24]分别对中国31个省市(自治区)水足迹(其中包括污染水足迹)和北京市灰水足迹进行了核算评价。我国灰水足迹研究尚属初级阶段，灰水足迹测算细节方面还有待进一步完善、区域综合角度鲜有涉及，同时未能从全国角度对灰水足迹的空间关联规律进行有效分析研究。

空间自相关分析是一系列空间数据分析方法和技术的集合，用来定量分析事物在空间上的依赖关系并将其可视化表达，目前主要集中应用在区域经济、土地利用、生态保护等方面，近年来在水资源领域也得到一定运用。鉴于此，本文运用灰水足迹相关理念，结合前人研究的基础上，选取1998—2012年中国大陆31个省市(自治区)为研究对象，在测度其灰水足迹的基础上，构建灰水足迹荷载系数并基于时间、空间两个维度运用空间自相关模型对中国省际灰水足迹荷载系数的动态演变趋势及其空间关联格局进行分析，旨在为解决区域水资源问题提供合理依据。

空间自相关分析是ESDA(Exploratory spatial data analysis)模型的中心内容之一，指研究范围中同一观测值在不同空间位置上的相关性，是研究特定空间范围聚集程度的一种度量^[25]。其主要包括空间权重矩阵、空间自相关系数等内容。空间自相关系数包含全局空间自相关系数和局部空间自相关系数。

在研究全局空间自相关的系数中，文章选取全局Moran′s指数做进一步研究。全局Moran′s指数可以辨别研究范围内观测值的聚集或分散特征。定义如下：

点击浏览公式

全局Moran′s指数只能判断观测值在整个研究范围的关联与差异程度，但不能判断研究范围内部的具体空间集聚特征及其显著性。局部Moran′s指数的提出解决了这一问题。其定义如下：

点击浏览公式

根据灰水足迹的概念，借鉴Arjen Y.Hoekstra等^{[5, 6]}编著的《水足迹评价手册》中的灰水足迹计算公式并结合前人研究的成果，将灰水足迹按人类生产生活的不同方面计算：

农药、化肥的大量施用及畜禽粪便的随意排放都会造成不同程度的水污染。据《第一次全国污染源普查》资料，农业已成为水污染物排放的主要源头。根据农业生产的主要方面，将农业灰水足迹分为种植业灰水足迹和畜禽养殖业灰水足迹两大类。

(1)种植业灰水足迹

农作物在生长过程中施加化肥、喷洒农药等除被植物利用外，部分在降水、灌溉等作用下，经过地表径流和地下淋溶进入水体，对水体造成污染。区别于工业与生活部门的水污染，种植业主要是面源污染，它以扩散方式进入水体。综合考虑种植业中的主要水污染物，选取氮肥作为水污染物进行灰水足迹计算。参照《水足迹评价手册》的常用模型，假设施用氮肥中固定比例(氮肥淋失率)的氮进入水体做简单的估算。其公式如下：

点击浏览公式

(2)养殖业灰水足迹

作为农业生产的重要组成部分，随着人口的增长、饮食结构的改变等原因，中国养殖业迅速发展，但畜舍冲洗、粪便排放等造成了严重水污染，据《第一次全国污染源普查》，畜禽养殖业已成为农业水污染的主要来源。选取具有代表性的畜禽(包括猪、牛、羊、家禽)养殖排污作为考量对象，通过养殖数量、饲养周期、排泄系数、单位粪便水污染物含量和进入水体流失率计算出畜禽养殖业污染负荷。为避免重复计算，饲养周期小于1年的猪与家禽数量取年末出栏量；饲养周期大于等于1年的牛、羊的数量取年末存栏量^[27]。一般而言污水中包含多种形式的污染物，灰水足迹由其中最关键的污染物决定即造成灰水足迹最大的决定^{[5, 6]}。结合单位粪便水污染物含量，选取COD、总氮作为指标评价畜禽养殖业的灰水足迹。由于水体可以同时对COD和氮元素进行稀释，因此，选择以两者为不同衡量指标的灰水足迹中的最大值作为畜禽养殖业灰水足迹。计算公式表达如下：

点击浏览公式

点击浏览公式

点击浏览公式

(3)农业总灰水足迹计算

通过以上对种植业灰水足迹和养殖业灰水足迹的计算，结合计算的特殊性，农业总灰水足迹的计算公式如下：

点击浏览公式

工业作为国民经济的重要组成部分，它的发展程度已成为衡量一个国家综合实力的主要标志，但是大力发展工业的同时也引起了大量的污染。较农业部门灰水足迹，工业源废水是直接排放进入水体，属点源污染，故可直接测算废水中主要污染物排放量，而COD和氨氮是排放污水中的主要污染物，所以，工业部门采用COD和氨氮作为衡量指标。其计算公式如下：

点击浏览公式

点击浏览公式

由于生活污水与工业污水同属于点源污染且排放污水中均是COD和氨氮为主要污染物，所以生活灰水足迹(WF_dom-grey)的计算同工业灰水足迹。

区域总灰水足迹由农业灰水足迹、工业灰水足迹、生活灰水足迹三部分组成。总结可得区域总灰水足迹计算公式如下：

点击浏览公式

灰水足迹并不能完全反应出一个地区的水环境压力，如部分地区虽灰水足迹相对较小，但地区水资源缺乏从而导致其水环境压力较大，基于此本文构建了灰水足迹荷载系数。灰水足迹荷载系数是稀释水污染到一定环境水质标准所用水量与水资源量的比值，建立起灰水足迹与水资源总量之间的关系，用于表示水污染压力的系数，系数越大，水污染压力越大。

点击浏览公式

本文农业部门的氮肥施用量和畜禽养殖数量数据来源于《中国农村统计年鉴》、《中国畜牧业年鉴》和《新中国五十年农业统计资料》；氮肥淋失率选取全国平均氮肥淋失率7%^[28]；畜禽饲养周期、排放系数、单位粪便COD和总氮排放量、进入水体流失率数据取自于《全国规模化畜禽养殖业污染情况调查技术报告》；工业废水和生活污水的废水排放量及其COD和氨氮量数据来源于《中国环境统计年鉴》；水资源总量数据选取《中国水资源公报》及各省市(自治区)统计年鉴。

文中污染物浓度达标排放标准采用《污水综合排放标准》(GB8978—1996)中一级排放标准，COD和氨氮的排放达标浓度分别为60、15 mg/L(综合文中总氮也采用氨氮标准考虑)。受纳水体的自然本底浓度指自然条件、无人为影响下水体中某污染物的浓度设为0^{[5, 6]}。

通过以上各个部门的灰水足迹计算，得到1998—2012年中国31个省市(自治区)的灰水足迹见表 1(限于页面宽度仅列偶数年份)。如图 1所示，1998—2006年全国灰水足迹呈现波动趋势，由1998年的5078.58×10⁸ m³上升至2006年的5238.53×10⁸ m³，并达到研究期间内全国灰水足迹最大值；2007年开始，全国灰水足迹呈现下降趋势，下降至2012年的4400.85×10⁸ m³同时达到研究期间内全国灰水足迹最小值，说明随着经济发展、技术水平的进步、环保意识的增强等原因，我国灰水足迹开始有所减小。

图1 中国1998—2012年灰水足迹结构与变化图 Fig.1 The composition and change of grey water footprint from 1998 to 2012

图选项

由图 1可知，全国灰水足迹中农业灰水足迹所占比例最大、生活次之、工业最小。同时随着时间的变化，农业灰水足迹虽在2007年降低，但2007年之后灰水足迹又出现上升趋势；生活灰水足迹相对较稳定，虽然随着人口的增长、人民生活水平的提高在研究期间开始有所上升，但随着环保意识的提高、技术水平的进步，生活灰水足迹开始下降；工业灰水足迹在研究期间呈现明显的下降趋势，这得益于国家可持续发展观念的提升、工业结构的调整等原因。由以上分析可知，我国工业、生活部门的点源污染得到有效控制、农业部门的面源污染呈现扩大趋势。

全国灰水足迹在地区上呈现出人口大省、农业大省灰水足迹较大的特征，河南、山东、四川、河北等地区灰水足迹明显高出全国平均水平，其中河南省灰水足迹为全国灰水足迹量最大的地区，平均达到387.47×10⁸ m³；北京、天津、上海、海南等地区灰水足迹较小，其中北京为全国灰水足迹最小的地区，平均仅为22.27×10⁸ m³。

由公式计算得到1998—2012年全国31个省市(自治区)灰水足迹荷载系数见表 2(限于页面宽度仅列偶数年份)。31个省市(自治区)15a灰水足迹荷载系数变化不大，整体呈现小幅波动趋势。在全国内部也存在着明显的地区差异，将31个省市(自治区)15年平均灰水足迹荷载系数分为五类(图 2)，第一类为高荷载区包括天津、宁夏两地，灰水足迹荷载系数远高于全国平均水平，水污染压力巨大，两地区虽然灰水足迹较低，但宁夏研究期间平均水资源量位列全国末位、天津其次，导致其水污染压力巨大；第二类为较高荷载区包括北京、河北、河南、山东、山西、上海6省市，除上海外均位于华北平原、黄土高原缺水地区，灰水足迹荷载系数均大于1，水污染情况不容忽视；第三类为中度荷载区包括辽宁、江苏、甘肃、吉林、内蒙古、陕西、安徽、黑龙江、湖北9省市(自治区)；较低荷载区包括重庆、湖南、广西、四川、贵州、新疆、广东、浙江、海南9个省市(自治区)，主要集中于长江以南水资源丰富地区；第五类为低荷载区包括福建、江西、云南、青海、西藏5地区，主要分布于青藏高原与南方水资源丰富地区，水污染压力较小。

图2 全国灰水足迹荷载系数分类图 Fig.2 The classification map of grey water footprint load coefficient in China

图选项

灰水足迹荷载系数受灰水足迹与水资源总量共同决定，通过分析可知灰水足迹荷载系数高值主要集中在华北与西北部分地区，低值主要集中于青藏高原与江南地区。灰水足迹荷载系数表现出一定的空间集聚性，即相邻省份的灰水足迹荷载系数有存在空间相关的可能性，故本文运用空间自相关分析方法进一步探讨中国省际灰水足迹荷载系数的空间关联格局。

文章将中国31个省市(自治区)研究期间的灰水足迹荷载系数运用全局空间自相关模型分析得到其全局Moran′s指数I值(表 3)。经检验，2007年之后全局Moran′s指数I值未通过显著性检验即不具有统计学意义上的全局自相关关系，2007之前(除2000年)全局Moran′s指数I值的正态统计量Z值均大于在0.05置信水平的下的临界值(1.96)，即通过显著性检验。说明中国省际灰水足迹荷载系数高值和低值分别呈现一定程度的空间集聚。动态分析可知，15a全局Moran′s指数虽然出现波动，但整体呈下降趋势，直至2008年呈现不显著相关。表明中国省际灰水足迹荷载系数全局空间相关性减弱，省际灰水足迹荷载系数的逐渐从空间集聚向分散转变。

通过以上分析可知：①区域水资源总量直接影响其灰水足迹荷载系数，相邻省份水资源量受相同地形、气候等因素共同作用下，呈现一定的相似特征从而使灰水足迹荷载系数出现空间集聚性：②同时灰水足迹也是区域灰水足迹荷载系数的决定因素，受经济发展水平、人口、产业结构等社会经济因素影响，临近省份社会经济特征也呈现一定的相似性，故导致灰水足迹荷载系数在空间上集聚；③随着区域发展水平的差异，区域间差异逐渐变大，影响下灰水足迹荷载系数的空间集聚现象开始减弱。

在全局自相关分析的基础上，利用Matlab软件计算出中国31个省市(自治区)灰水足迹荷载系数局部Moran′s指数I_i值，并做出1998年与2012年的中国灰水足迹荷载系数的LISA集聚地图(图 3)。如图 3所示：高高集聚(H-H)与低低集聚(L-L)地区较多、高低集聚(H-L)与低高集聚(L-H)地区相对较少，中国灰水足迹荷载系数在空间上集聚现象明显，省际水污染压力存在明显的空间关联性；但2012年H-H集聚与H-L集聚地区较1998年减少，说明中国灰水足迹荷载系数空间聚集现象开始减弱；省际水污染压力联系减弱。

图3 全国31个省市(自治区)灰水足迹荷载系数的LISA集聚分布 Fig.3 LISA cluster map of the grey water footprint load coefficient in China

图选项

(1)H-H集聚地区

稳定存在H-H集聚的地区分别为：天津、河北、山西、山东、河南，主要集中在华北地区。天津虽然在产业结构、土地面积、国家政策等影响下灰水足迹小，但地处华北地区水资源短缺，从而导致其灰水足迹荷载系数大，水污染压力巨大。河北、山东、河南三省作为我国的人口与农业大省，粮食生产过程中农药化肥的过量施用、能源原材料产业比重大、污水处理不足等因素影响下致使其灰水足迹远远高于全国平均水平，研究期间平均灰水足迹分列全国第三、四、一位，导致现有水资源无法有效稀释水污染。山西位于我国缺水地区，研究期间平均水资源总量属全国倒数第五，同等污染排放下水污染压力大。

(2)H-L集聚地区

宁夏稳定位于H-L集聚地区，属高度荷载地区。宁夏位于西北缺水地区，是全国水资源总量最少的省份，研究期间平均水资源量仅为9.9×10⁸ m³，水污染稀释可用水量少，灰水足迹荷载系数大。

(3)L-L集聚地区

吉林、黑龙江、浙江、福建、江西、湖北、湖南、广东、广西、海南、重庆、四川、贵州、云南、青海、西藏、新疆17省市(自治区)稳定位于在L-L集聚，大部分位于南方、青藏高原等地区。南方与青藏高原区是中国水资源丰富地区，新疆虽位于我国西北地区，但区域内水资源分布极端不均衡，北疆地区受地形等自然条件影响，水资源丰富。此类型由于可用于水污染稀释的水量大，故其相对缺水省市(自治区)灰水足迹荷载压力较小。同时新疆、青海、西藏地区由于自然条件等限制下经济发展相对落后、人口密度小，导致灰水足迹也较小。

(4)L-H集聚地区

L-H集聚区较为不稳定，内蒙古、陕西、辽宁、甘肃稳定位于该区。内蒙古畜牧业发达、养殖业灰水足迹量大，影响下其1998—2012年平均灰水足迹位列全国15位，灰水系数荷载系数位列全国13位。但内蒙古临近河北、山西、宁夏等灰水足迹系数更大的地区，相比较而言为低值区。陕西、甘肃位于西部地区，工农业发展相对落后，水污染相对较小，灰水足迹荷载系数较小，同时临近宁夏等H-H集聚区，故存在L-H集聚。辽宁省属中度荷载地区，临近河北省等较高荷载区，故存在L-H集聚。

北京市由H-H集聚区落入L-H集聚区，北京市位于华北缺水地区，且人口密度大、工业发达等原因影响下水环境压力较大，但在发达的经济水平、先进的技术水平、政策驱动等因素下研究期间灰水足迹呈现明显下降趋势，相较临近的天津、河北等地灰水足迹荷载系数变低；安徽省由L-L集聚地区变为L-H集聚区，由于其位于中国灰水足迹荷载系数H-H集聚区与L-L集聚区过渡地区，易受二者影响；江苏省由L-H集聚区变为H-L集聚区，与安徽相同其也位于中国灰水足迹荷载系数H-H集聚区与L-L集聚区过渡地区，易受二者影响，同时江苏省灰水足迹荷载系数由1998年的0.45增加到2012年的0.50，导致其由L-H集聚区变为H-L集聚区；上海由H-L集聚区落入L-L集聚区，与北京类似，上海作为我国的经济中心，在经济高速增长的同时，工业废水和生活污水的排放量大，水污染压力大。同时邻近南方水资源丰富地区，故位于H-L集聚区，但在发达的经济水平、先进的技术水平、政策驱动等因素下研究期间灰水足迹呈现明显下降趋势，使其落入L-L集聚区。

通过本文研究，可以得出如下主要结论：

(1)1998—2012年全国灰水足迹呈现波动趋势，1998年至2006年的灰水足迹呈现上升趋势；2007年开始，全国灰水足迹呈现下降趋势；农业在总灰水足迹的贡献率最高、工业最低；河南是灰水足迹最大的省份，北京是灰水足迹最小的地区。

(2)31个省市(自治区)15年灰水足迹荷载系数整体呈现小幅波动趋势。在全国内部也存在着明显的地区差异，大体分为五类，第一类为高荷载区包括天津、宁夏两地；第二类为较高荷载区包括北京、河北、河南、山东、山西、上海6省市；第三类为中度荷载区包括辽宁、江苏、甘肃、吉林、内蒙古、陕西、安徽、黑龙江、湖北9省市(自治区)；较低荷载区包括重庆、湖南、广西、四川、贵州、新疆、广东、浙江、海南9个省市(自治区)；第五类为低荷载区包括福建、江西、云南、青海、西藏5地区。

(3)借助全局与局部空间自相关对全国31个省市(自治区)15年的灰水足迹荷载系数进行时空演变与空间分布特征的分析可知，中国省级灰水足迹存在空间集聚现象且集聚现象逐渐减弱，依据计算结果得到H-H集聚区主要包括天津、河北等省市(自治区)，H-L集聚区包括宁夏，L-L集聚区主要集中在我国南方与青藏高原区，L-H集聚区较为不稳定，主要包括内蒙古、陕西等省市(自治区)绘制LISA集聚地图并进行分析。

本文通过灰水足迹及灰水足迹荷载系数的计算与空间关联格局分析，为水资源利用和水污染防治提供了一定的理论基础，可以通过对不同地区、不同污染物灰水足迹的计算，分析水污染原因并结合水资源量制定合适有效的污染物控制指标，以进行有效水污染防治，实现水资源的合理利用才能实现区域的可持续发展。

但本研究还存在一定不足，对灰水足迹的研究只取简单和代表性的污染物，许多方面估算可能偏小。种植业灰水足迹不只受氮肥及农药施用等的影响，且氮肥淋失率只是简要整理而得未充分考虑土地利用类型、气候地貌基础等对氮肥淋失的影响；养殖业灰水足迹计算中为了便于计算只简单考虑了主要畜禽粪便的水污染；整个面源污染考虑不足，如：农业地膜、水土流失、城市生活垃圾、工业废弃物都可能会随降水进入水体造成污染^[29]。