黄河流域是中华民族最主要的发源地之一, 也是我国重要的生态屏障和经济带^[1]。由于生态本底脆弱, 经济发展方式与生态保护要求相矛盾, 黄河流域也是我国生态环境敏感地区。自1978年以来, 我国陆续在黄河流域实施了“三北防护林工程”“退耕还林还草工程”等重大生态工程, 流域内生态环境得到较大改善。然而, 伴随社会生产、居民生活的快速发展, 流域超载的面积和人口比重分别为14.23%、17.22%^[2], 城市群污染排放占比超过70%^[3]。城市作为经济发展和生态治理的载体, 城市生态问题逐渐成为黄河流域生态保护和高质量发展的现实阻碍。城市生态治理的核心任务被认为应由减轻脆弱性向增强韧性转变^[4], 在已有研究中, 学者们主要关注黄河流域城市生态安全^[5]、生态脆弱性^[6]、生态系统服务^[7]、生态效率^[8]、绿色发展^[9]等, 较少从韧性视角探讨城市生态系统对不确定威胁的吸收和转化能力。科学测度黄河流域城市生态韧性、探究时空演变特征、剖析其发展的异质性对于实现黄河流域全面推进生态保护与建设“韧性城市”具有重要意义。

目前学者们主要从城市生态韧性概念界定、指标体系构建与方法选择、影响因素等方面展开探索。(1)概念内涵方面。韧性内涵经历了从工程韧性到生态韧性再到演化韧性两次重要的嬗变。“工程韧性”和“生态韧性”均属于均衡论视角^[10]; 演化韧性基于演化论视角, 将韧性看作是系统的固有属性, 无论是否受到干扰, 其始终具备干扰发生前进行防御、干扰发生时及时响应和不断学习创新的能力^[11]。基于此, 学者们对于城市生态韧性的概念界定大致可以分为两类：①均衡论视角：城市生态韧性指城市生态系统对风险干扰的预防、响应能力, 及在遭受干扰后返回干扰前状态的的反弹恢复能力^[12—14]。②演化论视角：学者不再强调城市生态系统恢复到干扰前状态的反弹能力, 更强调系统通过调整结构、改变路径实现转型发展的能力^[15—16]。(2)指标体系构建与方法选择。鉴于目前学界对城市生态韧性的内涵界定仍未达成共识, 其在指标体系构建方面也各有差异。国外研究多从生态系统服务^[17]、自然环境可持续发展^[18]、社区抗灾^[19]等视角进行分析; 国内研究侧重构建多维指标体系进行量化测度, 具体包括规模-密度-形态^[14]、抵抗-适应-恢复^[20]等多个目标层。部分学者结合压力-状态-响应模型, 通过核心变量法、适应循环模型、综合评价等探析区域生态韧性的动态演化^[13]。(3)影响因素。目前对城市生态韧性影响因素的研究主要分为自然因素和社会因素两类。自然因素的研究侧重自然生态系统本身, 例如, 气候、土壤、水文、植被和地形地貌等^[21]; 社会因素主要包括城市化、土地开发、科技创新等^[22]。部分研究认为城市生态韧性是社会活动与自然生境共同作用的结果^[23]。

总结既有研究发现：(1)现有文献对城市生态韧性的研究多是侧重均衡论, 较为缺乏从演化韧性视角对城市生态韧性进行探讨, 难以客观认识城市生态系统的现实水平与发展转型能力。(2)虽然已有研究成果基于地理学角度对城市生态韧性的时空格局进行了分析, 但关于城市生态韧性空间关联特征及动态转移规律的研究鲜有涉及。(3)在对城市生态韧性影响因素的研究中, 已有研究多基于障碍度模型^[15]、地理加权回归模型^[20]等方法, 以上方法虽然可以探究城市生态韧性的关键驱动因素, 但均忽略了影响因素对不同分位点上城市生态韧性的差异性作用。鉴于以上分析, 本文依托2011—2020年黄河流域七大城市群59个地级市, 基于演化韧性视角, 从“抵抗-响应-转型”三个能力构建城市生态韧性指标体系, 探究时空演变特征和动态转移规律, 利用面板分位数回归探究主要因素对城市生态韧性的影响差异性。以期全面把握黄河流域城市生态韧性的时空特征和影响机理, 为处于不同韧性等级的地级市推进生态保护和建设“韧性城市”提供一定的科学依据和理论参考。

1 研究区域概况与数据来源 1.1 研究区域概况

黄河流域横跨青海、四川、甘肃、宁夏、内蒙古、陕西、山西、河南及山东9省, 包括3个区域级城市群(山东半岛城市群、中原城市群和关中平原城市群)和4个地区性城市群(兰西城市群、晋中城市群、呼包鄂榆城市群和宁夏沿黄城市群)。流域内上中游城市群, 面临经济发展缓慢与生态本底脆弱的双重瓶颈, 下游城市群产业倚能倚重、环境污染问题突出, 城市污染排放约占沿黄九省区的80%^[24]。因此, 本文依据《黄河流域生态保护与高质量发展规划纲要》, 选取黄河流域七大城市群内地级及以上城市为研究对象。考虑到数据的可得性, 剔除了青海、甘肃的部分自治州(盟), 最终确定了59个地级及以上城市(图 1)。

1.2 数据来源

本文统计数据主要来源于《中国城市统计年鉴(2012—2021年)》《中国统计年鉴(2012—2021年)》以及各地级市的统计年鉴和统计公报。部分处理数据(人均水资源量、人均公园绿地面积、环境规制强度指数)主要由指标间的运算得出; 归一化植被指数(NDVI)为30m分辨率NDVI数据, 来源于中国科学院资源环境科学数据中心(http://www.resdc.cn/); 底图来源于自然资源部标准地图服务网站(https://www.mnr.gov.cn)审图号为GS(2019)1822号的标准地图, 底图无修改; 部分缺失数据通过线性插值法补足。

2 研究方法与模型构建 2.1 研究方法 2.1.1 熵权-逼近理想解排序法(熵权-TOPSIS法)

本文采用熵权-TOPSIS法对城市生态韧性作出评价。熵权-TOPSIS法在数据标准化的基础上依据数据的离散程度对指标进行赋权, 再采用TOPSIS法计算出各评价对象所有指标值与最优方案和最劣方案的距离从而得出历年的综合评价指数^[25]。具体步骤如下：

(1) 熵权法确定标准化指标权重, 公式为：

(1)

(2)

式中, w_j为第j项指标的权重; e_j表示第j项指标的信息熵; m为评价年数; y_ij为第i年第j项指标的原始数值; k为波尔茨曼常量。

(2) 确定正理想解与负理想解, 计算欧式距离：

(3)

(4)

式中，f_j⁺与f_j^-分别为第j个目标到最优目标及最劣目标的距离，f_ij是第i个目标第j个评价指标的权重规格化值；S_i⁺为目标到理想解的距离，S_i^-为目标到反理想解的距离。

(3) 计算各指标与理想解的贴近程度C值, 公式为：

(5)

式中, C_j取值范围为[0, 1], 其值越大, 则研究单元的城市生态韧性值越高, 反之越低。

2.1.2 传统马尔科夫链

传统马尔科夫链是一种时间和状态均为离散的马尔科夫过程, 能够反映出区域内每个个体的状态及其向上或向下转移的可能性。本文首先将测度出的城市生态韧性数据离散化, 变成k种类型, 接着计算每种类型状态的概率分布及转移变化, 最后得出城市生态韧性的整体演进规律。具体计算步骤如下^[26]：

(6)

式中, M表示k×k的马尔科夫转移概率矩阵; m_ij元素表示t年份i类型的地级市在t+1年份转移到类型j的概率, 并采用m_ij=n_ij/n_i进行估计。

2.1.3 空间马尔科夫链

空间马尔科夫链是在传统马尔科夫链基础上加入“空间滞后”的一种创新模型, 能够深入分析面板数据在时间演变和空间背景下的演变规律^[14]。本文将上述得到的传统k×k阶状态转移概率矩阵M以初始年份的空间滞后类型为条件, 分解成k个k×k转移条件概率矩阵, 分析不同邻域背景下, 城市生态韧性发生转移的概率^[27]。具体计算步骤如下：

(7)

式中, Lag_a是地级市a的空间滞后值; Y_b为地级市b的观测值; W_ab为空间权重矩阵(相邻为1, 反之为0)。

2.1.4 面板分位数回归

面板分位数回归通过分析影响因素对不同分位点城市生态韧性的影响, 提供观测变量之间更多的信息, 可以避免均值回归的理想化模式, 使研究结论更为可靠^[26]。具体步骤如下：

(8)

式中, Q_Yit为城市生态韧性的条件分位函数; Y_it为城市生态韧性值; a_i为常数项; X_it为解释变量矩阵; β(τ) 为T分位下的影响系数; T为本文设定的分位点。

(9)

式中, β(θ) 为影响系数; q为分位数数组个数; k为分位数数组的第k组; ρ_τk为分位数损失函数; w_k为第k分位数的权重系数; β(τ_k) 为第k分位数的影响系数。

2.2 城市生态韧性指标体系构建

“演化韧性”的本质源于适应循环理论, 该理论指出系统在不确定环境中会经历“利用-保存-释放-重组“四个阶段, 在每个阶段系统内部呈现出的不同潜能、应对能力及重组转型能力使得系统韧性高低各异^[28]。依据这个框架, Walker等提出韧性是系统为回应压力和限制条件而激发的抵抗、应对和改变的能力^[29]; 李阳力等认为水生态系统韧性指水生态系统具有的防御、响应和转换能力^[16]; 丁建军从抵抗与恢复能力、适应与调整能力、创新与转型能力三个维度构建了城市经济韧性指标体系^[30]。因此, 本文基于“演化韧性”视角, 将城市生态韧性视为城市生态系统的基础属性, 具体指城市生态系统通过内部抵抗、综合响应及发展转型, 最大限度预防和吸收干扰, 实现可持续发展的能力; 并参照已有相关研究^{[15, 30]}, 建立了包含抵抗能力、响应能力、转型能力三个方面的城市生态韧性指标体系, 见表 1。其中, 抵抗能力表示系统依托自然条件禀赋, 所具备的承受干扰、控制功能和结构不变的能力; 响应能力指系统在受干扰的过程中, 多样化应对冲击的能力; 转型能力体现系统在应对干扰的基础上, 通过学习、创新实现全新发展的能力。

2.3 城市生态韧性影响因素

本文参考相关研究^{[20, 23]}, 综合考虑黄河流域城市发展特点和数据可得性, 分析经济发展水平、人口密度、产业结构、科技创新、环境规制对城市生态韧性的影响程度。具体见表 2：(1)经济发展水平。一方面经济发展是城市生态文明建设的物质支撑; 另一方面, 经济发展水平的提高意味着生产规模和消费需求的扩大, 间接带来能源消耗和污染排放。(2)人口密度。依据新增长理论, 人力资本为城市生态建设提供了充足的劳动力和高素质智库; 另一方面, 人口扩张挤占城市生态资源, 加剧生态承载压力。(3)产业结构。产业结构是经济活动和生态环境的重要连接体。(4)科技创新。技术创新水平的提高能够为生态保护和破解资源消耗提供新途径。(5)环境规制。“波特假说”认为科学的环境规制可以产生显著“创新效应”; “遵从成本”理论提出环境规制会增加生产成本, 不利于企业绿色生产。本文使用环境规制强度指数作为代理变量, 具体计算公式为：

式中i为城市, t为年份, j∈ 1, 2, 3为污染物类型, 分别取工业废水、工业二氧化硫、工业烟(粉)尘, Polu为污染物排放量, Y为第二产业产值, ER_{i, t}为i城市t年的环境规制强度, ER越大则污染排放量越小, 环境规制强度越大。

3 黄河流域城市生态韧性的时空特征 3.1 黄河流域城市生态韧性的时序演变特征

图 2报告了2011—2020年黄河流域城市生态韧性的均值与变异系数。研究期内, 城市生态韧性由0.092增至0.125, 年均增长0.33%, 表现为缓慢增长态势; 变异系数先由0.307(2011年)增至0.545(2017年)再降至0.482(2020年), 表现为倒“U”型变化特征。十年间, 济南、青岛、烟台、西安和郑州始终高于黄河流域平均水平, 其均值分别为0.374、0.375、0.330、0.193和0.182;与之相反, 白银、中卫、阳泉和临汾始终低于黄河流域平均水平, 其均值分别为0.043、0.050、0.051和0.052。从城市群来看, 2011—2020年, 山东半岛城市群一直遥遥领先, 城市生态韧性处于(0.138, 0.217)范围内, 表现为显著提升态势; 中原城市群呈现持续稳定增长趋势, 处于(0.063, 0.094)范围内; 呼包鄂榆和关中平原城市群较为稳定, 变化范围分别为(0.071.0.082)、(0.067, 0.089);晋中和兰西城市群时序演变分别呈现“V”和倒“U”型变化特征; 研究期内, 宁夏沿黄城市群由0.089降至0.071, 呈现持续降低态势。

3.2 黄河流域城市生态韧性的空间演变特征 3.2.1 空间分布特征

借助Arc GIS技术平台将2011年、2014年、2017年和2020年黄河流域城市群城市生态韧性进行空间可视化分析(图 3)。

从地级市尺度看, 研究期内, 城市生态韧性主要以兰州、郑州、西安、太原、呼和浩特等单核心省会城市向外围地级市递减, 呈现出明显的“中心—外围”的空间格局。2011—2020年, 兰州、郑州、西安、太原、呼和浩特的城市生态韧性均值分别为0.119、0.182、0.193、0.124、0.105, 显著高于黄河流域平均水平(0.092)。其中, 郑州变化最为显著, 十年间增幅为208.5%;银川由0.082降至0.066, 是唯一一个出现负增长的省会城市。研究期内, 山东半岛城市群各地级市显著领先于其他地级市, 形成了以济南和青岛为中心的“双核领跑”格局, 其均值分别为0.374和0.375。

从城市群尺度看, 研究期内, 黄河流域城市生态韧性总体表现为“下游强、上中游弱”的分布格局, 山东半岛城市群始终最优, 均值为0.185, 呈现“一枝独放”的态势。2011年, 宁夏与晋中城市群仅次于山东半岛城市群, 分别为0.089、0.086;同年, 兰西和中原城市群相对脆弱, 分别为0.064、0.063。2014年, 呼包鄂榆与兰西城市群表现出增强趋势, 分别增长至0.074、0.075;晋中城市群则由0.086降低至0.066。2017年, 七大城市群城市生态韧性均有所提升, 其中兰西城市群增幅最大, 为33.3%;中原、晋中和关中平原城市群较低, 分别为0.072、0.070和0.071;至2020年, 以上三大城市群表现为显著增长局面, 增幅分别为30.5%、22.8%和29.6%;兰西与宁夏沿黄城市群分别下降至0.087和0.071。

3.2.2 动态演变特征

基于上文城市生态韧性的时空分布特征, 将城市生态韧性分为低韧性(0—0.100)、中等韧性(0.101—0.150)、较高韧性(0.151—0.300)、高韧性(0.301—0.555)4种类型。

表 3报告了未纳入地理空间因素的黄河流域城市生态韧性的动态演化特征：(1)马尔科夫转移概率矩阵中对角线上的数值均大于非对角线上的数值, 反映出城市生态韧性具有“路径依赖”与“自身锁定”的发展特征。(2)低韧性、中等韧性和较高韧性表现为向相邻高等级转移, 高韧性表现为向相邻低等级转移, 四个韧性等级在其他非对角线上的转移概率均为0, 说明城市生态韧性转型主要发生在相邻等级之间, 仍未实现跨等级转移。(3)低韧性等级的地级市维持现状的概率为0.9486, 向中等韧性转移的概率仅为0.0514;中等韧性向较高韧性转移概率为0.4545;由较高韧性向高韧性转移的概率仅为0.0815;由高韧性向较高韧性下移的概率为0.0145。

表 4报告了纳入地理空间因素的黄河流域城市生态韧性的动态演化特征：(1)邻域背景对城市生态韧性转移变化具有重要影响。不考虑空间邻域背景时, 低韧性地级市维持现状的概率为0.9486, 而与低韧性、中等韧性、较高韧性和高韧性的地级市相邻时, 其维持现状的概率分别为0.9487、0.8397、0.8797和0.7642。(2)不同的邻域背景对城市生态韧性转移的影响各不相同。中等韧性的地级市与低韧性和中等韧性地级市相邻时, 其向上转移的概率分别为0.0370和0.1667, 与较高韧性和高韧性的地级市相邻时, 其转移为较高韧性的概率分别为0.4545和0.7344;较高韧性的地级市与低韧性、中等韧性和较高韧性的地级市相邻时, 其向上转移的概率分别为0.0211、0.0200和0.0769, 而与高韧性地级市为邻时, 其转移概率升至0.5475, 表现为非对角线元素数值大于对角线上的概率。(3)高韧性地级市与低韧性地级市相邻, 维持现状的概率为0.2193, 向下转移的概率高达0.7807, 与中等韧性、较高韧性和高韧性地级市相邻时, 维持现状的概率分别为0.9630、0.9746和0.9848。

3.3 黄河流域城市生态韧的影响因素分析

参考穆学青的研究成果^[26], 设置10%、25%、50%、75%和90%为面板分位数回归模型的分位点。同时, 为避免数据出现内生性问题, 本文亦进行普通最小二乘法(OLS)回归进行比较。

表 5和图 4分别报告了面板分位数回归结果及可视化结果。经济发展水平(ED)对城市生态韧性具有显著正向提升作用, 在75%分位数水平下达到最高。人口密度(PD)对城市生态韧性的影响具有明显差异性。在10%、25%分位数水平下, 人口密度表现为显著正向促进作用; 在50%、75%分位数水平下, 人口密度的正向影响未通过显著性检验; 在90%分位数水平下, 人口密度表现为显著负向抑制影响。产业结构(IS)对城市生态韧性表现为显著负向影响, 说明黄河流域以第二产业为主导的产业结构对城市生态韧性提升具有抑制作用。科技创新(TI)对城市生态韧性具有显著正向提升作用。环境规制(ER)对城市生态韧性的影响作用存在水平差异。在10%、25%、50%分位数水平上, 环境规制对城市生态韧性表现为不显著的负向影响, 在75%、90%分位数水平上, 环境规制表现为显著正向影响。综上, 经济发展水平和科技创新是提升城市生态韧性的重要驱动因素; 当前, 不合理的产业结构抑制了城市生态韧性的增强; 人口密度的增大会阻碍高韧性地级市城市生态韧性的提升; 高韧性地级市应继续重视环境规制的正向推动作用。

4 讨论

通过上文的分析, 黄河流域城市生态韧性呈现“下游强、上中游弱”的空间布局, 山东半岛城市群始终最优, 主要原因是山东半岛城市群逐渐成长为北方地区重要的增长极, 通过“腾笼换鸟”等战略推进了产业绿色转型, 加之南山北水优越的自然条件, 共同筑起了生态城市建设的领先地位^[11]。至2017年, 七大城市群城市生态韧性均得到提升, 兰西城市群增幅最大, 这得益于国家在该区域长期积极实施退耕还林还草、水土流失治理等生态工程, 城市生态系统抵抗能力和响应能力增强^[31]。因此, 该阶段城市生态韧性的变异系数在2017年开始逐渐减小, 呈现倒“U”型变化特征; 至2020年, 中原、晋中和关中平原城市群增幅分别为30.5%、22.8%和29.6%, 这主要得益于《促进中部地区崛起规划(2016—2025年)》, 为恰处于中心地带的城市群建设生态文明城市提供了政策倾斜^[32]。尽管如此, 城市生态韧性的提升非一日之功, 在加大城市生态系统保护、修复投入力度的同时, 注重对边缘城市的政策支持。近十年来, 城市群内部形成了以省会城市为核心向外围地级市递减的“中心—外围”空间格局, 即城市群交界地带容易形成城市生态韧性“洼地”。究其原因, 城市群间生态保护和治理存在“以邻为壑”的现象, 边缘城市是最容易接受到污染转嫁的地方。

由于空间格局和历史发展的特殊性, 黄河流域跨越青藏高寒、西北干旱半干旱、华北湿润半湿润等若干个自然分区, 呈现显著的自然地理分异格局和东西发展梯度差异^[33]。因此, 受限于自然和社会双重复杂环境, 各因素的影响作用存在显著差异。科技创新表现为正向促进作用, 与“科技创新是改善黄河流域生态环境的重要动力”的研究论断一致^[34]。产业结构对各分位点的城市生态韧性均存在显著的负向抑制作用, 与已有研究“工业活动是造成地区生态压力的显著影响因素”的观点相符^[35]。究其原因, 黄河流域约50%以上的地级市属于资源依赖型和重工业主导型城市, “高污染、高耗能、高排放”的工业发展模式影响城市生态建设和可持续发展^[36]。在75%、90%分位数水平下, 环境规制表现为显著正向影响, 证实了“环境规制对生态韧性存在正向贡献”的研究论断^[20]。伴随城市生态韧性从分位点10%上升为90%, 人口密度的影响逐渐由正向作用转变为负向抑制, 可能的原因是, 城市生态韧性提升初期对劳动力需求量大, 注重人口“量”的供给, 当城市生态韧性提升至一定水平, 对劳动力的需求从“众”向“优”转变, 大量人口流入可能会造成空间拥挤、资源抢占、污染加剧等问题。

5 结论

(1) 城市生态韧性由0.092增至0.125, 年均增长0.33%, 表现为缓慢增长态势; 变异系数由0.307增至0.545再降至0.482, 表现为倒“U”型变化特征。空间上, 城市生态韧性表现为“下游强、上中游弱”; 城市群主要以兰州、郑州、西安、太原、呼和浩特等省会城市为核心向外围地级市递减, 呈现出明显的“中心—外围”的空间格局。

(2) 城市生态韧性类型转移具有“路径依赖”与“自身锁定”效应, 其转型主要在相邻等级之间进行递次转移, 仍未实现跨等级转移; 地理空间背景对城市生态韧性转移变化具有重要影响, 高韧性地级市的正向溢出作用更突出, 低韧性地级市更倾向抑制相邻地级市城市生态韧性的向上转移。

(3) 从城市生态韧性的影响因素来看, 经济发展水平的显著正向影响在75%分位数水平上达到最高; 科技创新能够有效促进各分位数水平城市生态韧性的提升; 不合理的产业结构严重抑制了城市生态韧性的增强; 人口密度在10%、25%分位数水平上表现为显著正向促进作用, 对高韧性等级地级市城市生态韧性提升存在抑制影响; 不合理的产业结构严重抑制了城市生态韧性的增强; 环境规制仅在75%、90%分位数水平上表现为显著正向影响。