陈青松,鲁亚雄,罗 勇,张洪吉,谭小琴,朱 琛,费建波,杨钰赐
(1.四川省自然资源科学研究院,四川 成都 610041;2.四川农业大学资源学院,四川 成都 611130;3.四川省自然资源厅,四川 成都 610072;4.四川省林业和草原调查规划院,四川 成都 610084;5.四川省国土空间规划研究院,四川 成都 610081)
资源环境承载力是衡量资源环境状况受到人类生产生活干扰能力的一个重要指标,相关研究也由来已久[1],总体上经历了从生物生态学到资源环境科学的综合认识过程[2,3]。目前,通过构建指标体系开展资源环境承载力评价是常用方法[4,5],而运用熵值法对指标进行权重赋值能有效排除主观因素的干扰[4,5];基于上述方法,学者开展了不同研究尺度的资源环境承载力的综合评价研究工作,如国家尺度[6-8]、省域尺度[9-11]、区域尺度[12-16]、市县尺度[17-19]等,在分析资源环境承载力基础上借助变量深入分析其驱动因素研究还相对较少[16,20];地理加权回归模型(Geographically Weighted Regression,GWR)是常用的空间异质影响因素分析模型[21],可用于驱动因子的筛选[22];相较于大尺度的资源环境承载力评价,县域内的评价更需要精细的基础资料才能有针对性地指导资源利用和环境调控[23],因此在县域开展资源环境承载力评价研究更具有理论与实践探索必要性,有利于完善和改进现有的评价框架和承载力评价方法,将空间管制和空间治理相结合,为城镇空间扩展、生态空间保护、生产要素和资源的空间配置以及空间规划编制等研究提供参考。
仁寿县地理坐标为103°55'—104°30'E、29°38'—30°20'N,隶属于眉山市,位于四川盆地南部,人口密度高、人均耕地少、人地矛盾突出,周边与6市区接壤,北接成都市双流区,是西南经济百强县,全域进入天府新区辐射区、影响区,经济发展较快。全县幅员面积2 716.86km2,东西长73km,南北长55km,地势由西北向东南逐渐倾斜,整体南北高、中间低、呈谷峰型,海拔高度470.6—1 383m。辖32个乡镇(街道),县城建成区面积超过40km2,据统计2020 年仁寿县常住人口数量111.1 万人,县内分布有彝族、藏族、羌族等16 个少数民族,地区生产总值为457.4 亿元。仁寿县自然资源丰富,拥有煤炭、石灰石、石英砂等矿产资源;同时也是产粮大县,果蔬之乡;拥有野生动物100 余种,野生有利用价值的植物560 余种;地表水沿龙泉山山脉东西分流,县内流域面积可达2 200 余km2,年均径流量约10 亿m3;县内全域呈现南北高中间低的地势,境内地质构造单元处于川西台陷龙泉褶皱与川中台拱、威远穹隆的接合部位,属地质灾害易发县。仁寿县属亚热带湿润季风气候区,年均气温17.3℃,年均降水量1 018.7mm,年均日照时数1 179.4h,年均相对湿度77%。
选取仁寿县2016 年开始融入成都天府新区到2020 年发展成型的两个年份,搜集了相关自然和社会经济的研究数据,包括指标层中25 个指标(表1),数据收集整理以乡镇为单元,数据来源于仁寿县规划和自然资源局、生态环境局、水利局及农业农村局等部门,经过数据申请和走访得到,同时查阅仁寿县统计年鉴、眉山市统计年鉴等对数据进行补充。
图1 仁寿县区位Figure 1 Location of Renshou County
1.2.1 基于PSR模型的指标体系构建
“压力—状态—响应”(Pressure - State - Response,PSR)被广泛用于资源环境承载力的评价,该模型最早有加拿大统计学家David J.Rapport 和Tony Friend两人提出,后经经济合作与发展组织(OECD)和联合国环境规划署(UNEP)逐渐发展为研究环境问题的框架[24],该模型是结合自然生态环境和社会经济相互内在互动关系,以一定的逻辑关系,对区域资源环境承载力进行评价,其动因则是人类的各类生产和生活活动给生态环境所带来的各种扰动,人类的扰动又使得生态环境出现相应的不平衡状态,从而影响人类生存和发展,一段时间持续的扰动下生态环境无法通过自我调节维持这种平衡状态,从而人类社会必须根据生态环境的状况做出相应的应对措施,对其进行评估后采取相应的调整措施,以帮助生态环境维持平衡状态。其中压力是人类带来的扰动,如资源消耗、环境污染等,这种扰动甚至可称之为对资源环境的胁迫,本文分别选取4个方面的内容构建10 种压力型指标,为负指标。状态,指的是在一定时间下生态环境状态的变化,本文构建9 种状态指标,为正指标;响应则是人类社会所采取的减少扰动应对措施,能够体现为生态环境提供保障的指标,共计6 种,为正指标。
以PSR模型为基准构建与仁寿县人口、资源、环境和社会经济等要素和资源环境承载力密切相关的因素,同时在参阅相关文献的基础上[24,25],遵循科学性、综合性、代表性和层次性的原则,构建仁寿县资源环境承载力评价体系。评价指标体系包括目标层、准则层、分准则层、指标层等4 个层次。基于PSR模型及层次划分原则,结合典型县域资源环境承载力综合指数分为压力评价指标、状态评价指标和相应评价指标,共25 个指标因子构建仁寿县资源环境承载力评价指标体系(表1)。
1.2.2 熵权法构建资源环境承载力指数
随着权重确定方法研究的日趋成熟,国内外学者发展了层次分析法、特尔菲法、回归系数法、等差法和熵权法等权重确定方法。相对于其他权重确定方法,熵权法能够避免主观色彩对评价结果的影响[20],是一种客观的赋权方法。
首先,根据各指标评价数据建立矩阵Rij=n ×m,一般认为所统计的指标为一般分布;利用极值标准化(Min—Max)法对评价的指标进行标准化处理,从而消除指标间的内部差异[26],其标准化方法如下:
式中:Rij为第j个指标第i 个项目的初始值;为对应第j个指标标准化后的值;Rj-max和Rj-min分别为第i个指标的最大值和最小值。
其次,计算各指标的比重Pij:
式中:n为项目个数。
最后,在获得项目的熵值后,计算第j个指标的熵权:
式中:wj为项目的熵权值;m为指标数量。
本文采用多指标加权函数模型来描述县域资源环境承载力指数。即将单一指标标准化后的数值乘以与之对应的权重再相加求和的方法得到一个综合性指标数值,这个数值即为综合承载力指数[27],取值范围在0—1 之间,分值越高,表明资源环境承载力水平越高,具体计算公式如下:
式中:RECC为资源承载力综合评价指数;Wj为第j个指标的权重;为第i年的无量纲标准化值;n为指标个数。
1.2.3 全局空间自相关
为探究县域资源环境承载力的空间变异情况,运用GIS的空间统计分析,全局空间自相关分析,在给定的空间要素即相关属性的值,用以评估所表达的模式为聚类模式、离散模式或者是随机模式,Z和P 值的得分表示统计学的显著意义,而莫兰指数(Moran's I)的值则指示着空间上的聚集和离散,即莫兰指数为正值,则得到空间聚集的统计学意义,若莫兰指数为负值时则指示空间上为离散的趋势。本研究以空间自相关分析资源环境承载力的空间关联和空间差异,在空间统计学中,常用的度量空间自相关程度的统计指标为莫兰指数,其计算公式为:
式中:n为样本量;xi、xj为空间位置i 和j 的观察量;Wij表示空间位置i 和j 的邻近关系,当i 和j邻近时,Wij=1;反之为0。全局Moran's I指数的取值范围为[-1,1],大于0 为空间正相关,小于0 为负相关,等于0 为不相关。
1.2.4 地理加权回归
地理加权回归(GWR)是一种广泛应用的空间分析方法,该方法通过建立空间局部点的局部回归方程,可以探究所研究的对象在一定空间尺度下的变化和相关驱动因素[22,28]。因此本研究选择GWR模型评估相关驱动因子在空间上对县域资源环境承载力指数的影响。同时可以进行预测,而一般的回归模型不考虑空间性,但在实际问题中,回归参数往往是以从不同地理位置上采集的,有位置上的变化,为了考虑回归参数在不同地理位置上的变化,引入GWR模型进行不同指标的回归分析,以得到资源环境承载力的空间驱动力。GWR 模型是普通的回归模型的扩展,加入了空间特性,并以距离加权的方式引入到模型中,具体模型如下[29]。
式中:yi为第i 个空间位置上的被解释变量;(ui,vi)为第i 个样点的空间坐标;βio为第i 个样本点的常数估计值;βij(ui,vi)为第i个样本点第j个自变量系数,与其空间位置关系有关;xij为第j 个自变量在样本i的值;εi为残差。
通常假定其服从的是独立正太分布,在空间上个点回归系数的球解如下:
式中:X、y为各样本点的自变量与因变量矩阵;W(ui,vi)=diag(Wi1,Wi2,…,Win)为样本点i 的空间权重矩阵;=(xi1,xi2,…,xin)。其中,X =为由解释变量观测值构成的设计矩阵;为由被解释变量观测值构成的列向量;元素dij,表示空间位置点j 到回归点i 的欧氏距离;D 是从回归点i到其第M个最近邻居的距离;M是最近邻居的最优值。空间核函数类型选择高斯函数法,高斯函数用一个连续单调衰减函数表示权重Wij与距离dij之间的连续单调递减函数,可以克服上述空间权函数不连续的缺点,函数形式如下所示:
带宽则是由Akaike信息准则,可以通过使该值达到最小来确定[30],计算公式如下:
通过熵值法确定的各资源环境承载力指数情况,运用“压力—状态—响应”模型构建了32 个乡镇2016 年和2020 年的PSR 各子系统的承载指数值,结合县内的资源环境概况,对不同的子系统进行综合分析评价。
2016 年仁寿县乡镇的整体情况如图2 所示。从图2 可见,各乡镇压力指数与状态指数较为接近,都表现出较高的水平,但总体的压力指数水平小于状态指数,突兀的峰值出现在县城所在的街道,即怀仁、普宁、文林及视高街道,上述乡镇均表现为经济发展迅速,土地供应量大,人口不断涌入,人地矛盾较突出,资源消耗和环境压力较大,因而表现出较高的压力指数;响应指数方面,与压力和状态指数不同,存在压力和状态指数高而响应指数低的状况,如宝马镇;其中黑龙滩镇的响应指数最高,说明在该时期有相关减轻资源索取消耗和环境整治的措施,通过响应年份的仁寿县统计年鉴发现,2016 年对黑龙滩镇中的黑龙滩水库进行了整治,且对位于仁寿城区内的中央水体公园、仁寿城市湿地公园、响水六坊公园开园投用,减轻了资源与环境压力。
图2 2016 年仁寿县各乡镇PSR指数Figure 2 PSR index of each township in Renshou County in 2016
2020 年各乡镇的PSR指数如图3 所示。从图3可见,压力和状态指数在县城所在乡镇任然表现较高,但与2016 年相比压力指数的变化较大,整体呈现下降的趋势;状态指数较2016 年相比在个别乡镇上表现较明显,整体有一个增加的趋势,具体在方家镇、怀仁、普宁、视高及文林街道等乡镇,其中变化突出的方家镇和视高街道,方家镇作为撤乡并镇的排头兵,行政区划的调整使得镇内资源利用和产业发展提升较大;响应指数方面整体呈现减小的趋势,其中县城所在的怀仁街道、普宁街道以、文林街道及视高镇表现出较高的响应指数,上述街道的响应指数和压力指数呈现同减的趋势。
图3 2020 年仁寿县各乡镇PSR指数Figure 3 PSR index of each township in Renshou County in 2020
通过熵值法确定2016—2020 年仁寿县资源环境承载力综合评价指数变化情况,结果如图4 所示。从图4 可见,2016—2020 年仁寿县各乡镇的资源环境承载力整体上增加的乡镇多,达到23 个,整体上减小的乡镇数仅9 个。但存在研究期内个别乡镇资源环境承载力剧烈变化的情况,如黑龙滩镇;人口集中、资源消耗较大的县城所在街道的资源环境承载力均出现增加,这表明人口用地高度集中的城市在面对资源环境承载压力不断增加的同时,注意提升城市自身对资源环境承载力的响应能力,合理调配资源,注重保护城市环境,其资源环境承载力是可以优化提升的。
图4 2016—2020 年仁寿县资源环境承载力指数的变化情况Figure 4 Changes in the resource and environment carrying capacity index of Renshou County,2016—2020
通过空间自相关统计,两个年份的莫兰指数统计如表2 所示。由表2 所示,两个年份的的莫兰指数均大于0,说明研究中各乡镇的资源环境承载力在空间上表出聚集的趋势。在Z 值的得分上,当Z值得分大于2.58 时,随机产生此高聚类模式的可能性要小于1%,而Z 值得分大于1.96 时,则随机产生此高聚类模式的可能性要小于5%,本研究中2016 年的Z值为2.164,该值大于1.96 小于2.58,则本研究资源环境承载力的空间表现指示随机产生高聚类模式的可能性小于5%,这说明在研究中各乡镇的资源环境承载力值的空间分布表现聚集性,且高聚集模式概率超过95%。2020 年的Z 得分为1.55,小于1.96,没达到显著性要求,聚集的概率小于2016 年。
表2 莫兰指数统计表Table 2 Moran's I index statistics
各乡镇的空间聚类分析如图5 所示,可以看出,2016 年资源环境承载空间上无显著聚集性的乡镇数为27 个,有显著聚集性意义的乡镇数为5 个,约占总乡镇的15.6%,其中高高聚类集中在县城辖区的街道,如文林街道、怀仁街道;而2020 年的聚类P值为0.121,大于0.05,未达到显著水平,则表明该年度资源环境承载力的空间分布是随机空间过程的结果,不存在明显的聚类现象。与2020 年相比2016 年的县域内资源环境承载力变异度降低。
图5 仁寿县资源环境承载力空间变异特征分布Figure 5 Distribution of spatial variation characteristics of resource and environment carrying capacity
从高低聚集程度来说,2016 年出现了一定的空间聚集现象,其中高—高聚集在文林街道和怀仁街道,为县城所在区域,呈现高承载力,周边也是高承载力现象;两个年度中富家镇2016 年还呈现高低聚集的现象,而2020 年则表现为低低聚集,资源环境承载力的下降;大化镇两个年度均表现出低高聚集,即周围表现高资源环境承载,本镇表现为低资源环境承载力;低低聚集现象的乡镇在2020 年有所增加,对整体资源环境承载力的下降起到指示作用。
2.3.1 驱动力变量选择
为了弄清资源环境承载力变化的驱动因素,选择适当的变量进行研究。在选择研究变量时,主要是根据自变量与因变量之间的相关关系以及相关描述性统计分析结果来反映,将驱动力变量指标作描述性统计分析,结果能够反映原始指标体系数据的整体情况,本研究借助描述性统计和回归分析对2016 年和2020 年的指标变量进行了分析,统计得到了25 个指标的整体情况,2016 年的指标中各指标的Pearson相关系数值在-0.565—0.783 之间,25 个指标变量与资源环境承载力显著性统计显示,除了A13和A19,其余指标显著性均在95%置信区间,其中指标A16的统计量被排除在指标之外;2020年的指标中Pearson 相关系数值在-0.732—0.899之间,显著性统计显示,全部满足显著性在95%置信区间,其中指标A5的统计量被排除。
在模型构建之前,引入方差膨胀因子(VIF)对各指标的共线性进行判断,VIF 值越大则表明变量之间共线性越强,一般而言,VIF 值越小越能满足地理加权回归的要求,对于VIF 大于5 的变量则认为不满足共线性要求[21]。从表5 中选取VIF 小于5且显著性水平在95%置信区间的变量指标,从而筛选出构建地理加权回归模型的变量(表3)。
表3 驱动变量筛选结果Table 3 Screening results of driving variables
2.3.2 驱动力因素分析
研究对比了最小二乘法(OLS)与地理加权回归(GWR)的结果(表4),由表4 可知,OLS 的回归结果拟合度调整R2小于GWR 的拟合度,因此研究采用GWR的回归结果进行驱动力分析。基于筛选后的驱动因子,确定相关参数后,输入地理加权回归模型得到相应的驱动因素的拟合结果,整体结果包括AICc和调整R2,其反映了地理加权回归的可信度。其中2016 年筛选了3 个因子,2020 年筛选了5 个因子,两个年份的AICc值均小于100,表明模型拟合较好。将拟合结果输入空间可视化软件后,获得仁寿县资源环境承载力GWR 拟合结果空间分布,其中渲染得到2016 年局部拟合度(图6)。
表4 地理加权回归分析结果统计表Table 4 Statistical table of geographically weighted regression analysis results
图6 2016、2020 年地理加权回归的局部R2 渲染Figure 6 Local R2 rendering and standard error plots of geographically weighted regressions in 2016 and 2020
从图6 可见,2016 年的全局拟合度即响应的拟合效果,整体拟合效果较好,R2均超过了0.5,仁寿县北部乡镇的拟合度高于南部乡镇,北部乡镇出现了聚集的趋势。2020 年的拟合度整体也超过0.5,呈现自北向南逐渐降低的趋势,标准误差显示仅有视高街道、文林街道的误差相对较大,拟合效果不太理想,其余地区的回归拟合效果都比较理想,整体误差还是偏小的。
将地理加权回归所得系数分别输入到空间可视化软件进行可视化展示(图7),分别对GDP增长率(A2)、万元工业产值标准能耗(A23)和农药使用量(A25)进行渲染。由图7a可知,在汪洋镇、禄加镇和禾加镇的GDP增长率对资源环境承载力表现出负向作用,3 个镇处在仁寿县的最南端,产业结构和布局处在边缘区域,可能是表现负向作用的原因;而其余镇对整个县的资源环境承载力均表现正向作用,正向作用力从北向南逐渐增强;而万元工业产值标准能耗对资源环境承载力均表现为负向作用,自西北向东南的负向作用逐渐增强;农药使用量的回归结果表明,农药使用量对资源环境承载力整体表现为正向的驱动作用,但驱动力较小,呈现出东西方向的变化趋势,自东向西逐渐减小,农药化肥是农业生产的必须投入品,合理施用农药是资源环境有一定的正向作用。
图7 2016 年地理加权回归系数渲染Figure 7 Rendering of geographically weighted regression coefficients in 2016
同样对2020 年回归系数结果进行渲染,结果如图8 所示。由图8a 可知,2020 年的第二产业、第三产业占GDP的比例(A3)对资源环境承载力的均表现为负向的驱动力,二三产业对GDP 的贡献率越大,对资源环境承载力的压力也就越大,压力趋向降低资源环境承载力的大小,分布的趋势呈现从西南向东北逐渐增强的的趋势,即高家镇向东南方向至北斗镇表现较高的负向驱动力,该区域内是环天经济区的主要区域,与成都天府新区、东部新区接壤,二三产业的发展尤其是三产业的发展是能源消耗和环境压力的主要因素,因而表现较强的负向能力。由图8b可知,水资源开发利用率(A18)同样整体表现出对资源环境承载力负向的影响能力,强弱从西北向东南逐渐增强,最强的依然出现在汪洋镇、禄加镇及禾加镇等,该区域是县域内水环境工业污染重点控制区域,水资源的开发利用需要更注重防止工业污染。湿地是动植物多样性的集中表现区域,能够调节气候,缓冲环境污染的压力,是资源环境承载力正向的驱动力,图8c显示结果表明县区内湿地保有量(A22)对资源环境承载力的驱动均表现为正向的驱动影响力,正向趋势的影响力差异较小,但依然呈现自北向南逐渐增强的趋势,这种格局可能是因为北部靠近成都的部分乡镇的湿地对资源环境承载力的正向驱动力要小于经济发展所带来的压力。对比图8d与图7b可以发现,2016 年的万元工业产值标准能耗(A23)相比2020 年,表现有所不同,整体呈现一种正向的驱动作用力,分布趋势呈现自西向东逐渐增强的趋势,说明产业升级,产业结构的调整,包括发展理念的转变,使得资源节约、环境友好型工业发展对资源环境承载力的驱动由负转正。同2016 年的农药使用量相比,化肥的使用表现除了不同的趋势(图8e),农用化肥使用量(A24)对资源环境承载力整体都表现出了负向的影响力,表现出自西南向东北逐渐增强的趋势,东北部分乡镇的农用化肥使用量对资源环境承载力的驱动作用更强,化肥更容易引起资源环境承载力的负向变化。
图8 2020 年驱动因素地理加权回归系数渲染Figure 8 GWR coefficient rendering of driving factors in 2020
本文以县域为研究尺度,通过构建指标体系,运用PSR模型和熵值法确定了资源环境评价指数,结合地理加权回归分析,开展了县域尺度资源环境承载力的时空变化及其驱动因素的实证研究。主要结论如下:①仁寿县2016 年资源环境压力和状态指数均处于较高水平,部分乡镇的环境响应指数高,人地矛盾突出的县城和经济增长迅速的乡镇街道资源环境承载力的压力指数明显大于其他乡镇;2020 年人口聚集和城市开发集中的乡镇压力指数有所缓解,同时响应指数和状态的增加明显,因此在城市开发过程中更加注重对环境的正向响应,有助于提高整体资源环境状态指数;2016—2020 年虽然各状态指数波动明显,但主要人口聚集和城市开发的乡镇资源环境承载力呈现增加的趋势,表明提升城市自身对资源环境的响应,合理调配资源,注重保护城市环境,其资源环境承载力是可以优化提升的。②仁寿县资源环境承载力指数在2016 年呈现出一定的聚集分布,主要为县城所在的文林街道、怀仁街道,呈现高高聚集;2020 年未呈现聚集分布,与2016 年相比低低聚集乡镇有所增加,资源环境承载力空间相关性下降。③仁寿县资源环境承载力驱动因素主要集中在经济和生态子系统内,拟合度2016 年优于2020 年,整体误差偏小;2016 年万元工业产值标准能耗是主要负向驱动力,GDP增长率和农药使用量是主要正向驱动力;2020 年二三产业占比、农用化肥施用和水资源开发是主要负向驱动力,湿地保有和万元工业产值标准能耗是主要正向驱动力;从整体的回归系数看负向驱动力强于正向驱动力,今后因在产业结构调整上多下功夫,合理水资源开发和农药使用,大力发展湿地生态,有助于提升资源环境承载力。
仁寿县作为成都市的正南方,独特的地理分布和自然资源禀赋为该区域的发展提供者源源不断的支撑作用,因此县域的未来高质量发展必须建立在合理保护的基础上,保护和开发都需要提升和优化,以提高资源环境承载力水平。依据资源环境承载力的分析结果,首先,应建立资源环境承载力底线约束机制。县域内生产生活与城镇建设是影响资源环境承载力的关键因素,一方面要明确三线一单制度,即资源利用上限、环境保护底线、城市开发红线以及资源环境负面清单,一方面要构建相应的保障制度确保制度的落实和实施;同时约束机制的建立还应考虑不同发展的历史时期下人们对资源环境总需求不降低的情况下,提高现有资源环境的开发效率。其次,持续优化资源环境承载正负向驱动。资源环境承载力的变化受到各类社会、经济和环境因素的影响,加强正向驱动的促进作用、减小负向驱动的抑制作用是维持资源环境承载力的途径之一,对于仁寿县,加强水资源开发利用的管控以及大力建设生态湿地无疑是维持资源环境承载力的重要方面。