成渝地区双城经济圈人类活动氮磷输入时空演变及其驱动因素

邓陈宁，张泽乾，徐 睿，李 虹，邹天森，杨鹊平*，刘伟玲，聂 冲，赵艳民

1. 中国环境科学研究院，北京 100012

2. 国家长江生态环境保护修复联合研究中心，北京 100012

3. 北京师范大学水科学研究院，北京 100875

随着社会经济的飞速发展，伴生的水污染已成为备受关注的热点问题[1]. 高强度的人类活动显著改变了氮磷等养分的地球化学循环，导致水体富营养化等一系列水生态环境问题[2-3]. 河流养分通量与人类活动养分输入之间存在较强的相关性，氮磷养分输出的比例分别为输入的20%～30%和3%～6%[4-6]. 人为氮磷输入主要来源于农业活动和食物生产消费[7]. 在农业集约程度较高的区域，作物施肥和畜禽养殖可能引入过量的养分输入[8]；在人口密度较高的城市群，食物生产消费方式和土地利用的改变也直接影响着养分循环[9]. 此外，气候条件、土地覆盖和社会经济发展水平的差异，也会导致各地区养分输入与环境响应过程不尽相同[10].

随着我国城市化进程的加快，诸多经济区和城市群等高复杂度的社会-经济-自然复合系统正在加速形成[11]. 目前，我国城市群大多处于发展的初级阶段，以水体富营养化为代表的环境问题尚未完全解决.城市群内跨市社会经济活动相互作用，使区域环境问题放大升级为跨区域问题[12-14]，亟待找寻科学的解决方案. 成渝城市群是我国城市群的典型代表，其作为长江经济带战略支撑之一，是西部大开发的重要平台和高质量发展的重要增长极之一. 《成渝地区双城经济圈建设规划纲要》明确提出，要共筑长江上游生态屏障，推动生态共建共保，加强污染跨界协同治理，探索绿色转型发展新路径. 因此，开展成渝地区双城经济圈人为氮磷输入研究对绿色可持续发展至关重要.

模型是研究区域氮磷输入和水环境响应的重要技术手段. 用于定量分析区域氮磷输入和输出的数学模型主要分为机理模型和统计模型. 机理模型以HSPF (hydrological simulation program-fortran)[15]和SWAT (soil and water assessment tool)[16]等为代表，统计模型以SPARROW (spatially referenced regressions on watershed)[17]等为代表. 然而，由于对输入数据的精度要求高，上述模型在数据缺乏或质量不佳地区的应用十分受限. 基于活动水平数据的准物料平衡模型净人为氮输入(NANI)和净人为磷输入(NAPI)可定量评估区域人为氮磷的输入[18-19]，为有限数据区域养分输入及环境响应估算提供了途径. 目前，NANI和NAPI已广泛应用于多尺度不同类型的区域[7,20]，成为定量评估河流氮磷输入的有效指标[21]. 为此，该研究基于NANI和NAPI模型，以成渝地区双城经济圈人类活动氮磷输入为研究对象，定量辨识其时空格局及环境响应过程，以期为城市群环境养分管理提供参考.

1 数据与方法

1.1 研究区域

成渝地区双城经济圈地处长江经济带上游地区的四川盆地，位于长江经济带与“一带一路”的交汇处，是我国西部地区人口最密集、产业基础最雄厚、发展水平最高的城镇化区域[22]. 成渝地区双城经济圈生态资源、水能资源、矿产资源丰富，总面积18.5×104km2，占全国的1.9%. 2019年，常住人口约9 600×104人，地区生产总值近6.3×1012元，分别占全国的6.9%和6.3%. 成渝地区双城经济圈以成都市和重庆市为中心，包括重庆市的中心城区及万州、涪陵、綦江、大足、黔江、长寿、江津、合川、永川、南川、璧山、铜梁、潼南、荣昌、梁平、丰都、垫江、忠县等27个区(县)以及开州、云阳的部分地区，四川省的成都、自贡、泸州、德阳、绵阳(除平武县、北川县)、遂宁、内江、乐山、南充、眉山、宜宾、广安、达州(除万源市)、雅安(除天全县、宝兴县)、资阳等15个市. 长江在区域内流经高海拔生态脆弱区、平原人口稠密区及工业聚集区，生态环境风险叠加交织.

1.2 研究方法

该研究基于张汪寿等[5,23]改进的NANI和NAPI模型对2011-2019年成渝地区双城经济圈人类活动氮磷输入进行估算. 该方法将养分输入分为非点源输入(NANIn、NAPIn)和点源输入(NANIp、NAPIp)，计算公式：

1.2.1非点源输入

人为非点源氮输入(NANIn)主要由农村地区食品/饲料净氮输入量(Nr-im)、氮肥施用量(Nfert)、作物固氮量(Nfix)和大气NOx沉降量(Ndep)组成. 人为非点源磷输入(NAPIn)主要由农村地区食品/饲料净磷输入量(Pr-im)、磷肥施用量(Pfert)和非食品磷输入量(Pnon)组成.

式中，Nhc和Phc分别为农村地区人类N、P消费量，Nlc和Plc分别为畜禽N、P消费量，Ncp和Pcp分别为作物产品中供人类食用的N、P量，Nlp和Plp分别为畜禽产品中供人类食用的N、P量(含水产品). 上述变量单位均为kg/(km2·a)，分别按N和P计.

由于食物产品在加工、运输、储存和零售等过程中存在损耗，因此作物以及畜禽、水产品不能被完全利用，取折损系数为10%[4]；化肥施用量根据单一肥料和复合肥施用量(折纯)估算[24]；作物固氮量根据共生固氮作物种植面积(豆类、花生)和非共生固氮作物种植面积(稻谷、其他作物)估算[5]；大气NOx沉降量根据氮氧化物排放量估算[25]；非食品磷输入量根据常住人口数估算.

1.2.2点源输入

点源氮磷输入(NANIp和NAPIp)主要为经处理后排放的生活污水和工业废水中的N、P量，计算公式：

式中：Nurb和Purb分别为通过生活污水输入的N、P量，通过城市地区常住人口数和人均氮磷消耗量估算[26-27]，单位均为kg/(km2·a)，分别按N和P计；Nind和Pind分别为通过工业废水输入的N、P，单位均为kg/(km2·a)，分别按N和P计；Isew为城市污水处理率，%；Irem-tn和Irem-tp分别为污水处理工艺对于N、P的平均去除率，分别取0.6和0.7[5,28].

1.3 数据来源

所用数据主要包括活动水平数据和相应参数. 活动水平数据包括2011-2019年研究区域涉及地市的城镇和农村常住人口、人均食品消费量(包括粮食、蔬菜、肉类、禽类、水产品、蛋类、奶类、干鲜瓜果)、作物产量及面积(稻谷、小麦、玉米、高粱、豆类、薯类、花生、油菜籽、芝麻、蔬菜、水果)、农用化肥施用量(氮肥、磷肥、复合肥)、畜禽养殖数量(牛、马、驴、骡、羊、猪、家禽、肉兔)、水产品产量、工业废水排放量、城市污水处理率、氮氧化物排放量等，主要来源于国家和地方统计局及其他政府部门发布的统计年鉴等. 作物产品和畜禽产品的氮磷含量、作物固氮率等相关参数来源于长江流域相关研究[5-6,24,29-34].由于研究区域中重庆市的开州、云阳以及四川省的绵阳、达州、雅安均为部分区域，故采用面积加权[35]的方法将市域尺度的数据进行分配.

2 结果与讨论

2.1 人类活动氮磷输入总量及其组成结构

2011-2019年成渝地区双城经济圈平均NANI和NAPI分别为13 063和2 291 kg/(km2·a)，高于长江经济带尺度[7]和全国尺度的平均水平[25,34-38]，但低于长三角和长江中游城市群[10]. 如图1所示，研究区域中，非点源和点源NANI分别为12 514和549 kg/(km2·a)，分别占NANI总量的95.8%和4.2%. 非点源NANI中氮肥施用量占比最高，为非点源氮输入量的60.3%，为NANI总量的57.8%. 非点源和点源NAPI分别为2 210和82 kg/(km2·a)，分别占NAPI总量的96.4%和3.6%. 与非点源NANI组分占比结构相似，非点源NAPI中磷肥施用量占比最高，为非点源磷输入58.8%，为NAPI总量的56.7%. 从NANI和NANI的组成结构来看，化肥施用量是主要组分，食品/饲料养分净输入量次之，与长江经济带尺度和全国尺度的氮磷输入结构相似[9,27,37]. 结果表明，成渝地区农业生产活动对于人为氮磷输入影响显著. 农业土地利用方式不合理、种养模式不科学、农用化肥施用强度过高等均会引起人为非点源氮磷的过量输入，造成氮磷在土壤中滞留，并伴随降雨引起的冲刷、淋滤等过程进入水环境，增大了环境污染风险[39]. 因此，可将“源头-过程-末端”相结合的养分管理模式作为深入研究方向和有效措施途径，实现经济效益与环境效益的双赢.

图1 人类活动氮磷输入总量及其组成特征Fig.1 Anthropogenic nitrogen and phosphorus inputs and their composition

2.2 人类活动氮磷输入时间轨迹及空间格局

2011-2019年成渝地区双城经济圈NANI和NAPI总量均呈下降趋势(见图2). 由图2可见，与2011年相比，2019年NANI和NAPI总量分别降低了43.2%和25.9%. 这主要是由非点源NANI和NAPI下降引起的，二者分别降低了44.9%和26.8%，但点源NANI和NAPI则分别上升了7.0%和3.3%. 对于非点源NANI，农村地区食品/饲料净氮输入量、氮肥施用量以及大气氮沉降量均呈明显下降趋势，分别下降了105.8%、13.4%和43.5%. 农村地区食品/饲料净氮输入量降为负值，表明该区域以出口食品/饲料为主，氮素随食品/饲料的输出流出研究区域. 氮肥施用量以及大气氮沉降量的下降与《到2020年化肥使用量零增长行动方案》和《大气污染防治行动计划》等的实施密不可分[12,40]. 对于点源NANI，经处理后生活污水输入量增加了12.9%，这是城市化进程加快，城市常住人口数量急剧增长的结果. 2019年成都市和重庆市城市常住人口数量比2011年分别增加了30.9%和30.0%. 经处理后工业废水输入量降低了45.5%，与工业产业升级和结构调整等有关. 对于非点源NAPI，磷肥施用量下降了13.4%，农村地区食品/饲料净磷输入量、非食品磷输入量也因农村常住人口减少等分别下降了45.4%和18.6%. 对于点源NAPI，其趋势与点源NANI相似，经处理后生活污水输入量增加了5.7%，经处理后工业废水输入量降低了49.4%.

图2 NANI和NAPI及其组分的年际变化Fig.2 Interannual variation of NANI, NAPI and their components

空间尺度上，2011-2019年NANI(见图3)和NAPI(见图4)平均强度较高的地区主要集中于区域中部，包括遂宁市、内江市、自贡市、资阳市，以及与成都临近的德阳市，与重庆市临近的广安市等. 德阳市NANI和NAPI均为最高，分别为30 451和455 kg/(km2·a)，遂宁市、内江市、广安市、自贡市、资阳市依次降低且均高于成渝地区双城经济圈NANI和NAPI总量. 雅安市NANI和NAPI均为最低，分别为4 801和546 kg/(km2·a). 非点源NANI和NAPI空间异质特性与NANI和NAPI总量的高度相似，N、P肥施用量空间格局与非点源NANI和NAPI高度相似.不难发现，德阳市、遂宁市、内江市、广安市、自贡市、资阳市等地区化肥施用强度均较高. 已有研究[1,7]表明，NANI和NAPI与耕地面积均呈正相关. 2019年上述各市耕地面积占比分别为42.0%、50.7%、50.8%、48.5%、49.4%和55.8%，单位面积粮食总产量分别达到330.5、267.5、317.0、283.9、316.7和289.1 t/km2，均处于较高水平. 点源NANI和NAPI最高的为成都市，分别为1 646和243 kg/(km2·a)，主要由城市生活污水主导. 成都市作为成渝城市群双核心之一，城市常住人口数远高于其他地级市. 与非点源NANI和NAPI最低值相似，雅安市点源NANI和NAPI也最低. 这与其环境底数、产业结构、人口分布等有关. 截至2019年，雅安市森林覆盖率高达67.38%，而NANI和NAPI与森林面积均成反比[2,9]. 诸如此类因素造成雅安市NANI和NAPI均呈现较低水平.

图3 NANI及其组分的空间异质性Fig.3 Spatial heterogeneity of NANI and its components

图4 NAPI及其组分的空间异质性Fig.4 Spatial heterogeneity of NAPI and its components

2.3 人类活动氮磷输入驱动因素

基于NANI和NAPI模型结构，人类活动氮磷输入各组分对其有直接影响. 由于研究区域内各地区社会经济发展水平不同，为量化不同组分对NANI和NAPI变化的相对贡献度，引入Lindeman-Merenda-Gold(LMG)模型. LMG模型基于多元回归模型，在考虑因变量之间相关关系和顺序效应的基础上，通过平均对变量所有可能的边际贡献，以此分解因变量的方差，进而对因素的相对贡献度进行计算[41-43]. 该模型基于Rstudio的“relaimpo”包[44]实现. LMG模型结果显示，2011-2019年研究区域内各地区农村地区食品/饲料氮输入量、氮肥施用量、作物固氮量、大气氮沉降量、经处理生活污水氮输入量、经处理工业废水氮输入量对于NANI的相对贡献度分别为24.5%、25.2%、17.3%、17.1%、6.1%和9.7%. 2011-2019年研究区域内各地区农村地区食品/饲料磷输入量、磷施用量、非食品磷输入量、经处理生活污水磷输入量、经处理工业废水磷输入量对于NAPI的相对贡献度分别为21.0%、30.4%、24.2%、14.0%和10.5%. 其中农用化肥施用相对重要性最高，食品/饲料养分输入量等次之.

此外，NANI和NAPI也受社会经济或自然因素的间接影响. 基于文献调研[1,27,45]，选取社会因素(总人口密度、城市人口密度、农村人口密度)、经济因素(人均GDP、单位面积GDP)及土地利用因素(耕地面积、森林面积)进行影响因素分析. Pearson相关性分析结果(见图5)表明，NANI和NAPI均与农村人口密度及耕地面积呈显著正相关，与森林面积呈显著负相关(由于土地利用变化速率较人口和GDP缓慢，因此土地利用选取2019年进行分析)，但其与人均或单位面积GDP均无显著关系.

图5 NANI和NAPI与影响因素的Pearson相关性分析Fig.5 Pearson correlation analysis of NANI, NAPI and driving factors

随着农村人口密度的增加，NANI和NAPI也呈上升趋势〔见图6(a)〕. 近年来，随着城市化进程的加快，部分农村人口转变为城市人口，使得NANI和NAPI中农村地区食品/饲料养分输入量降低，直接影响氮磷输入总量. 相关研究[1]表明，人口密度存在阈值，当人口密度高于100人/km2时，其对人为氮磷输入变异解释的能力趋于稳定. 此外，随着耕地面积占比的上升，人为氮磷输入强度也会上升〔见图6(b)〕，这也是德阳等市NANI和NAPI较高的原因. 然而，森林面积占比则与人为氮磷输入强度呈显著负相关〔见图6(c)〕，雅安市NANI和NAPI值较低也可反映这一特点.

图6 NANI和NAPI与主要影响因素的相关关系Fig.6 The correlation between NANI-NAPI and the main driving factors

3 结论与展望

a) 成渝地区农业生产活动对于人为氮磷输入影响显著. 2011-2019年，成渝地区双城经济圈平均NANI和NAPI分别为13 063和2 291 kg/(km2·a)，高于长江经济带尺度和全国尺度的平均水平，但低于长三角和长江中游城市群的强度水平. 其中，非点源NANI和NAPI中农用化肥输入量为主要来源，食品/饲料养分输入量次之.

b) 人为氮磷输入年际变化特征由城乡人口结构演变和农业生产措施优化驱动. 在时间尺度上，2011-2019年成渝地区双城经济圈总NANI和NAPI均整体呈下降趋势，其中NANI和NAPI总量分别降低了43.2%和25.9%，非点源NANI和NAPI分别降低了44.9%和26.8%，但点源NANI和NAPI分别上升了7.0%和3.3%.

c) 人为氮磷输入空间异质性受城市间农业发展水平和人口密度差异影响明显. 在空间尺度上，2011-2019年NANI和NAPI平均强度较高的地区大部分集中于区域中部以及成都市、重庆市的临近地市等.其中德阳市NANI和NAPI均为最高，分别为30 451和4 553 kg/(km2·a)，遂宁市、内江市、广安市、自贡市、资阳市依次递减，雅安市NANI和NAPI均为最低，分别为4 801和546 kg/(km2·a).

d) 人为氮磷输入模式受社会经济因素和自然因素的间接影响. NANI和NAPI与农村人口密度、耕地面积占比均呈显著正相关，与森林覆盖率均呈显著负相关.

e) 在未来的研究中，算法方面，可结合养分流动和养分输出，改进模型、优化参数、完善体系，从全链条对氮磷的社会经济循环进行剖析，基于经济与环境效益共赢视角，削减其对大气、土壤和水环境的直接与间接影响；管理方面，加强区域环境的协同治理研究和政策制定，可为城市群、跨界流域等的可持续发展提供保障和动能.