河流系统是社会生产、生活与生态环境可持续发展的重要基础^[1]。近年来, 以城市化为代表的人类活动导致了区域不透水面积扩张、水系数量锐减、连通性下降、河道干渠化、鱼类洄游受阻等现象^[2-5], 使得区域汇流时间缩短, 洪涝调蓄能力减弱, 易形成“小水大灾”、水质下降、生态环境恶化等问题^[6-10]。城市化导致的水系变化及其所引发的水文生态效应成为了国内外学者关注的热点问题。

不同的城市化发展阶段由于对河流水系功能存在不同需求, 人类活动对其进行改造所产生的影响也存在差异。缓慢城市化时期, 为保障农业高产, 人们对原有河道进行疏浚、裁弯取直、围湖造田等, 最终导致河网密度、水面率和河流曲度剧烈衰减^[11-12]；快速城市化时期, 为了保障城市防洪安全, 主干河道被拓宽, 加之城镇用地扩张导致大量末级河道被填埋, 河流水系复杂性急剧下降^[13]；而新时期人类改造河流还会受到水生态、水环境、水景观、水文化等理念因素影响。

目前, 针对城市化下水系变化的定量分析, 多是采用与河流地貌学相关的指标^[14]。Sear等发现受到城市化强烈干扰, 世界上60%的河流其发育过程都发生了变化^[15]；Chin等将河流对城市化的适应过程划分为反应阶段、张弛阶段和平衡阶段^[16]；Gregory等总结出在城市化的影响下英国Monk河流域河槽特征发生显著改变^[17]；袁雯等提出了土地利用类型变化带来的平原河网演变可能具有“自然型-井型-干流型”的发展趋势^[18]；吴雷等采用格网化的方式分析了城市化发展与河流水系演变的空间异质性^[19]；张凤等刻画了城镇体系和水系结构演化特征及其时空关系, 证实两者具有不同的时空演化方向^[20]。

分形理论成功揭示了水系结构特征中的由整体到部分、由宏观到微观之间的关系^[21-25]。然而流域系统十分复杂, 可以划分出许多地形、构造、气候等自然发育条件各异的子流域, 各个子流域有其单独的分形, 因此大型水系属于多重分维范畴^[26]。近年来, 学者针对多重分形方法展开许多有益尝试^[27-28], 在流域地形地貌刻画^[29]、水文过程模拟^[30]及流域景观特征分析^[31]等方面取得良好进展。

根据上述背景, 人们对城市化下水系演变特征有了更为深刻的理解。但现有多重分形研究大都用来刻画流域地貌形态特征, 对水系结构的描述也多集中在山区或者重点在揭示水系多分维对水文过程的影响等方面, 在平原河网地区应用较少。杭嘉湖地区地处太湖平原南部, 河道密布, 社会经济发展快速, 水系变化深刻, 加剧了区域洪涝灾害, 部分学者从水利片区^[12-13]、城区-郊区^[11]等宏观、微观尺度分析了城市化对河流水系变化的影响, 而以县级行政区为分析单元的研究仍较为少见。为此, 本文基于1 ∶ 50000地形图获取的1980s与2010s水系数据, 在分析水系数量和形态特征的基础上, 融合多重分形理论与地理加权回归模型, 探究平原河网地区城市化和水系演变之间的关系, 以期为保护我国城市化地区的河流水系、实现人类与河流系统协调发展提供科学参考。

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

杭嘉湖地区(30°09′—31°02′N、119°52′—121°16′E)位于太湖流域南部, 钱塘江及杭州湾以北, 黄浦江以西, 浙西丘陵区以东, 总面积约为7607万km², 是太湖流域面积最大的平原区, 且整个地势呈现西南高东北低的格局。自更新世以来该区经受过多次海浸、海退, 对地貌产生很大影响, 湖相、河相、海相交叉, 加之人类有史以来造田耕作, 形成水田和旱地镶嵌、河浜纵横交错的现状。研究区地处我国东部季风区, 属于亚热带季风气候, 雨热同期, 多年平均温度16—18℃, 多年平均降水在1100—1500 mm, 受季风强弱影响, 降水年际变化较大。区内河道纵横交错, 发育程度高, 是典型的平原河网区(图 1)。

1.2 数据与方法 1.2.1 数据获取

基于1 ∶ 50000地形图获取的1980s和2010s的河流水系资料, 分析杭嘉湖地区近30年来的河流水系演变特征。其中1980s水系数据由1 ∶ 50000纸质地形图经数字化处理获得, 2010s水系数据源于1 ∶ 50000数字线划图。由于研究区地势平坦, 河道密集且纵横交错, 难以按照传统的Strahler、Shreve等规则进行河流分级。综合考虑到河流的自然特性和社会功能, 本文一般将宽度大于20 m的河流划为一级河流, 即为主干河流；将宽度为10—20 m的河流划为二级河流, 同时也是一级支流；将宽度小于10 m的河流划为三级河流(二级支流)。其中, 干流在河网中主要起泄洪排涝的作用, 而支流主要起调蓄功能。遥感影像选取1991年、2015年的Landsat卫星数据, 分辨率均为30 m, 对其进行影像剪裁、拼接等处理, 根据研究区土地利用类型特点与研究需要将土地利用类型划分为城镇、水田、旱地、林草地和水域五类。经调查验证, 两期遥感影像解译精度均超过80%, Kappa系数达到0.80, 解译结果可信。采用1990年与2015年的社会经济数据, 数据源自中国科学院资源环境科学与数据中心(http://www.resdc.cn/)。其中, 人口空间分布网格数据集精度为1 km, 单位为人/km², 该数据集反映了人口数据在全国范围内的空间分布状况；国内生产总值(GDP)空间分布网格数据集精度为1 km, 单位为万元/km²。

杭嘉湖地区城市化发展较快, 人口密集, 河网水系受城市化影响较为深刻。本文以人口密度、GDP密度及空间城市化水平3个指标表征杭嘉湖地区的城市化发展(图 2)。1990—2015年研究区人口密度出现了较大幅度(101%)的增长, 东北和西南地区为近几十年来人口集聚速率最大, 集聚程度最高的区域。GDP密度总体增长幅度为33%, 各县区经济虽出现不同程度的增长, 但贫富差距趋于扩大化。空间分布上, 经济快速增长主要集中在东北部沿海开放区及西南地区。城市发展主要通过城镇用地扩张实现, 1990—2015年杭嘉湖地区空间城市化水平上升幅度达31%, 有明显提高, 但城市化发展的速率差异越来越大。城市化水平越高, 城镇用地扩张越慢；相反, 在城市化水平较低的区域, 城镇用地增幅较大。依据2015年空间城市化水平, 本文将杭嘉湖地区各县级行政单元划分为低度城市化区(≤30%)、中度城市化区(30%—50%)和高度城市化区(>50%)三类^[12]。

1.2.2 水系演变特征参数

河流水系演变特征主要体现在水系的数量、结构和形态等方面。多重分形特征不仅能够反映出水系结构整体的复杂程度, 更能够有效刻画水系分布的均匀性与集中度。本文综合国内外关于水系指标的基础上, 选取了河网密度R_d^[13]、水面率W_p^[11]描述水系数量特征, 结构特征通过计算水系多分维谱得到的奇异指数分布范围Δα及多重分形谱高差Δf表征, 形态特征选取河流曲度S_r^[11]、干流面积长度比R_m^[13]进行定量刻画。

分形维数的计算方法有两种：一种基于Horton定律, 其实质是分形递归, 但仅基于现象表述的角度开展研究^{[21, 32-33]}；另一种则是计盒法, 也是最为简单且应用最广的水系分形维数计算方法^[34]。本文以水系长度为目标量, 通过计盒法统计研究区水系长度分布概率。用边长为ε×ε的盒子覆盖整个河网水系, 统计第i个盒子中的水系长度总和M_i(ε), M_i(ε)与全流域所有盒子内水系长度值的总和M之比即为第i个盒子内的水系分布概率Q_i(ε)。公式如下：

(3)

式中, Q_i(ε)为第i 个盒子内的水系分布概率, M_i(ε)为第i 个盒子中的水系长度和, M为流域内的水系总长度。

在水系分布概率Q_i(ε)的基础上, 分别通过式4、式5与式6计算得到配分函数X_q(ε)、质量指数τ(q)。

(4)

式中, N(ε)为覆盖水系的盒子中所有非空盒子的数目；q∈(-∞, +∞), 为权重因子。q值不同反映水系分布概率Q_i(ε)在配分函数X_q(ε)中的具体作用不同。本文结合研究区实际情况, 设置q∈[-4, +4], 计算步长为0.5。

(5)

如果研究区内水系分布概率Q_i(ε)具有多重分形特征, 则配分函数X_q(ε)与盒子边长ε具有式5所表示的关系。对配分函数X_q(ε)与盒子边长ε分别取自然对数(式6), 其自然对数之比即为质量指数τ(q)。

(6)

通过对质量指数函数τ(q)-q进行Legendre变换, 得到多重分形谱函数f(α)-α(式7)。

(7)

式中, α(q)为奇异指数函数, 表征水系分布的不规则程度；f(α)为多重分形谱函数。

(1) 多重分形谱奇异指数分布范围Δα

在多重分形谱函数f(α)-α中, 奇异指数分布范围Δα定量地表征了分形体内最大概率子集和最小概率子集的对比关系, 反映了区域内部水系结构的差异化程度。也就是说, Δα数值越大, 各子集概率的两极化趋势越明显, 区域内部水系结构疏密差异越大, 分布越不均匀；Δα数值越小, 水系分布的均匀性越好。其计算公式如下：

(8)

式中, α_max与α_min分别表示奇异指数的最大值和最小值。

(2) 多重分形谱高差Δf

多重分形谱高差Δf可根据以下公式进行计算。

(9)

多重分形谱高差Δf主要表征流域内具有相同奇异指数α的水系结构变化趋势, 即水系分布的集中度。Δf > 0, 即最小概率子集的数目小于最大概率子集的数目, 区域内水系分布较为集中；若Δf < 0, 则研究区水系呈分散分布；若Δf=0, 则表示大小概率子集的数目相等, 二者对水系分布的影响相当。

1.2.3 地理加权回归模型

地理加权回归模型(Geographically Weighted Regression, GWR)是一种对区域内自变量和因变量之间的变化关系进行建模的非参数局部空间回归模型^[35-36]。该模型利用随着空间位置不同而发生变化的参数来估计结果, 量化反映数据之间的空间响应关系。其公式如下：

(10)

式中, y_i为i点的因变量值, x_ik(k=1, 2, …m)为影响因变量y_i变化的自变量, (u_i, v_i)为i点的坐标, β₀(u_i, v_i)为截距, β_k(u_i, v_i)(k=1, 2, …m)为斜率, 即回归系数, ε_i则为误差项^[37]。

2 结果分析 2.1 杭嘉湖地区水系特征及其演变规律 2.1.1 数量特征演变规律

1980s—2010s研究区平均河网密度减少态势明显, 衰减率为10%(图 3)。杭州市的衰减率达20.5%, 近30年该市河网密度降低了1/5, 可见城市的快速发展对自然水系的巨大冲击。就空间分布而言, 河网密度在1980s基本呈现自东向西逐渐减小的格局, 2010s分布则相对随机。且高值区(4.01—5.00 km/km²)从地区东部的集中分布演化为以平湖为核心, 在上海都市圈快速城市化影响下, 松江和金山的河网密度受到不同程度的影响而逐渐下降。太湖平原与钱塘江作为推动长三角经济发展的载体与纽带, 其沿岸及周边地区水系受城市化影响较大, 河网密度低值区逐渐趋于沿太湖和钱塘江分布。

水面率分布的总体格局较为稳定(图 3)。两个时期水面率中、高值区主要沿太湖分布, 而水面率较小的区域则集中在地区东部和南部, 大致沿钱塘江分布。从数值上看, 1980s研究区平均水面率为9.72%, 最大值和最小值分别为20.4%和3.4%；而2010s平均水面率为8.43%, 最大与最小分别为18.9%和1.2%。可见, 近30年杭嘉湖地区平均水面率降幅达13%, 衰减幅度较大, 最低水面率均在南湖区, 其降幅更加剧烈(64.7%)。总体上, 研究区水面率的空间分布格局改变不大；就变化率与数值而言, 区内湖泊等水体缩减明显。

2.1.2 多重分形特征演变规律

1980s—2010s杭嘉湖地区大部分区域的奇异指数分布范围Δα大于1, 相对较大, 说明水系分布不规则, 疏密程度差异性较高(图 4), 这是由于平原河网区的水系纵横交错, 分布情况比较复杂。从Δα变率看, 研究区Δα增大趋势明显, 表明杭嘉湖地区水系结构的均匀性持续变差, 分布更不规则, 以钱塘江沿岸最为突出。两个时期中以2010s吴兴区Δα最大, 为1.547, 与1980s相比增加了103%。中等大小的Δα(1.11—1.30)分布格局变化比较明显, 由1980s连片分布于地区南部和中西部演化为2010s基本沿钱塘江和以青浦、金山区为次级核心分布, 表明钱塘江沿岸以及青浦、金山区近30年来内部水系趋于破碎化。而Δα低值区(小于1.10)呈现出沿太湖分布的趋势。

1980s和2010s水系多重分形谱高差Δf平均值分别为1.04和1.05, 均大于1, 说明研究区的河网水系空间分布较为集中。就数值来看, 1980s和2010sΔf最小值为0.61和0.71, 最大值分别为1.27和1.34, 均呈增大趋势, 表明20世纪80年代以来水系集中分布趋势明显。空间分布上, 2010sΔf的高值区呈现往东南方向转移的趋势(图 4), 钱塘江沿岸Δf变率普遍较大, 且为正值, 说明沿岸区域水系分布空间集聚度较1980s更大。

2.1.3 形态特征演变规律

1980s—2010s杭嘉湖各区域的河流曲度均小于1.30, 同时仅1980s海宁市的河流曲度超过了1.13(图 5)。相反, 两个时期中大部分区域的河流曲度均介于1.05—1.10之间, 尤其是1980s, 16个区域中有14个其河流曲度低于1.10, 占比87.5%, 说明平原河网地区绝大部分河流属于顺直型河流^[38]。从河流曲度变化率来看, 研究区平均河流曲度减小了0.5%, 除嘉善、德清、吴兴、南浔四个区域的河流曲度在2010s出现微小的增长, 其余区域的河流曲度均呈略微下降的趋势。空间上, 河流曲度的分布格局变化主要表现为高值区(1.14—1.16)的消失与中值区(1.08—1.10)从较为集中的分布演变为相对随机的分布。

与河流曲度相比, 干流面积长度比的分布格局变化更大(图 5)。空间上, 干流面积长度比高值区(大于60.00)从以松江为核心演变为双核心分布(松江区和德清县), 较高值区(50.01—60.00)从西北部集中分布向以吴江区为核心演化。此外, 低值区(30.01—40.00)从随机分布变为集中分布在东南部沿钱塘江和杭州湾一带。总体看来, 中、高度城市化区域的干流面积长度比增加幅度较大, 如德清县与杭州市的干流面积长度比增幅分别为18.4%和15.2%。这是因为干流是城市发展规划的主要保护对象, 为了缓解区域洪涝, 拓宽主干河道与修建新的行洪河道成为城市可持续发展的必要选择。

2.2 城市化对水系演变的影响

杭嘉湖地区水系变化的空间分异现象较为显著, 传统统计模型(如线性回归、逐步回归、主成分分析等)在研究过程中无法真实反映城市化对水系影响的空间分异特征。由于人口与GDP数据存在较强共线性, 且城市化进程对河流水系造成直接影响主要通过城市土地扩张实现。因此, 本文分别构建各水系指标的变化率和空间城市化水平变化率之间的地理加权回归模型, 定量探讨城市化和水系变化之间的空间响应关系。

2.2.1 城市化对水系数量特征的影响

中度城市化区域, 河网密度显著降低。在高度与低度城市化区域, 水面率与城市化的空间响应更为敏感。R_d-Urban表示以河网密度变化率为因变量、空间城市化水平变化率为自变量构建的地理加权回归模型, 以此类推。从GWR模型回归系数的空间分布可以看出, 河网密度、水面率对空间城市化水平的响应关系并非完全表现为负相关。R_d-Urban关系在高度城市化与低度城市化区域为正相关, 在中度城市化区域则表现为负相关(图 6), 说明在中度城市化区域, 河网密度衰减较为明显；W_p-Urban关系在高度城市化与低度城市化区域为负相关, 在中度城市化区域表现为正相关(图 6)。类似于R_d-Urban响应关系, 呈负相关则表明在高度与低度城市化区域, 水面率受城市扩张影响更易发生缩减。

2.2.2 城市化对水系多重分形特征的影响

杭嘉湖地区城市化的快速发展对水系多重分形特征影响的空间异质性表现较为明显。在Δα-Urban和Δf-Urban模型(图 7)中, 回归系数的空间分布均表现为：在高度和中度城市化区域为负相关, 而在低度城市化区域则为正相关关系。这说明, 随着城市发展, 低度城市化区域水系多重分形特征表现出复杂性增强、疏密程度差异扩大的规律；在中度和高度城市化区域, 水系不规则性减弱, 趋于均匀化分布, 主要是由于城市建设用地快速扩张, 干扰并破坏了河网的自然发育过程, 导致大量末级支流被填埋, 复杂多元的水系结构趋于简单化。从回归系数大小来看, 城镇用地扩张能使杭嘉湖地区多重分形谱高差变化率最大增加0.1%—14.7%, 即水系分布的空间集聚度变化幅度较大。

2.2.3 城市化对水系形态特征的影响

杭嘉湖地区河流曲度S_r及干流面积长度比R_m的变化均与城市化发展呈正相关关系。从回归系数大小来看, 研究区城镇用地扩张能使局部区域河流曲度变化率最大增加0.1%—0.9%, 干流面积长度比变化率最大增加0.1%—4.5%(图 8)。就空间分布而言, S_r-Urban及R_m-Urban模型回归系数均存在从地区西南向东北方向逐渐增加的趋势。两个模型的回归系数从大到小均排列为：低度城市化区域>中度城市化区域>高度城市化区域, 表明城市化水平越低, 形态特征演变越剧烈。

3 讨论

近30年杭嘉湖地区的水系变化较为剧烈。以杭州市为例, 河网密度和水面率分别衰减20.5%、18.7%；河流曲度下降1.8%, 干流面积长度比增幅达15.2%；奇异指数分布范围减小11.1%, 表明水系分布趋于均匀化。然而, 从河网密度和水面率的角度看, 其河网水系呈衰减状态。这可能是由于在城市化等人类活动的影响下, 杭州市水系受扰动剧烈, 河道数量缩减, 河网密度与水面率有所减小。但由于被填埋的多为各级支流, 加之区域防洪需要, 主干河道保留较为完整且有所拓宽, 水系复杂性下降, 趋于简单化。

为进一步探究不同城市化水平在转变过程中对水系的影响, 基于地理加权回归模型统计得到“低-低、低-中、低-高、中-高”4种不同城市化水平转变类型下各水系指标的变化率和空间城市化水平变化率之间的局部回归系数(图 9), 从而反映不同区域城市化水平的变化对于水系的影响。在“低-低”类型下, 城市化发展水平较低, 水面率衰减以及形态特征演变更为剧烈, 且水系空间分布局部集中程度(Δf)增加趋势更为明显；而在“低-中”及“低-高”情况下, 水系分布存在向均匀化发展的变化趋势, 且城市化水平转变幅度越大, 河网密度变化越剧烈, 水系均匀化分布趋势越明显。在“中-高”类型下, 城市发展水平较高, 受水面率控制及城市景观生态发展需要等政策影响, 加之末级河道被填埋, 因此除河网密度受扰动变化明显外, 其余水系特征值变化较小。

4 结论

对太湖平原杭嘉湖地区水系变化及其与城市化的空间响应进行分析, 得到如下结论：

(1) 快速城市化背景之下, 杭嘉湖地区水系变化较为剧烈。数量特征上, 河网密度与水面率均呈减小的趋势, 近30年衰减幅度分别为10%和13%。多重分形特征上, 2010s相比于1980s水系分布均匀性降低, 空间集中程度增加。形态特征上, 研究区多为顺直型河流, 且河流曲度总体呈略微下降的趋势, 干流面积长度比在中、高度城市化区域增幅较大。

(2) 运用地理加权回归模型分析发现, 杭嘉湖地区水系变化对城市化的响应在空间上存在较为明显的异质性。城市化水平越低, 水系形态演变越剧烈。高度城市化区域, 水面率受城市化影响显著, 水系分布均匀性增强；中度城市化区域, 河网密度对城市化的响应更为明显, 水系复杂性减小；而在低度城市化区域, 以水面率与城市扩张响应更为敏感, 且水系分布呈现局部集中程度增加、不规则性增强的趋势。

(3) 杭嘉湖地区存在“低-低”、“低-中”、“低-高”、“中-高”四种不同的城市化水平转变类型。转变类型不同, 对水系的影响也存在差异。城市化水平未发生转变时, 水面率衰减和形态特征演变更为剧烈, 且水系空间分布集中程度呈增加趋势；而在“低-中”及“低-高”情况下, 水系有向均匀化分布的变化趋势, 且城市化水平转变幅度越大, 均匀化分布趋势越明显。“中-高”类型发展过程中, 受政策影响仅河网密度变化较为明显。

河流水系在人类活动的干扰下不断演变, 其中城市建设大量侵占各级河流, 河道缩减使城市洪涝问题加剧。为此, 人们又不断疏浚、拓宽主干河道, 使得水系结构趋于简单化。可见, 河流系统与人类活动是相互反馈与制约的, 二者在这种过程中相互协调, 达到河网水系演变适应城市化发展的平衡态。然而, 城镇体系(城市群)的发展使得河流水系演变更加复杂, 城镇体系的空间结构变化与水系结构的时空演变关系以及适应城市发展的各水系特征值的合理阈值区间仍待进一步探究。