岷江上游流域分形地貌对农村居民点分布的影响

项 清，黄 弘，于 欢，阚瑷珂2，①，何 杰

(1.成都理工大学地球科学学院，四川 成都 610059；2.青藏高原及其东缘人文地理研究中心，四川 成都 610059；3.成都理工大学旅游与城乡规划学院，四川 成都 610059)

流域地貌是地貌类型中的重要一类，流域地貌的形成主要是地表流水在流动过程中侵蚀地面，形成各种形态的侵蚀河谷，同时又将被侵蚀的物质沿途堆积，形成各类堆积地貌[1]。农村居民点作为乡村聚落的主要形态和重要人文景观，其空间分布状况是人地之间关系演变规律的一种重要表现形式[2]。地貌形态要素对于农村居民点的布局具有基础性影响，而流域地貌区聚落分布与其特殊的流水冲刷地貌也密切相关[3]。掌握流域地貌形态影响下的农村居民点布局规律，对流域尺度下的农村居民点整理规划和乡村振兴具有重要的指导意义。

目前国内外学者对农村居民点分布与地貌环境的关系做了大量研究，研究区域涉及丘陵区、平原地区、喀斯特槽谷区、山区等。如焦贝贝等[4]对黄土高原低山丘陵区农村居民点分布与地形因子关系进行了研究。李珊珊等[5]以秦巴山区为例，分析了高程、坡度、地形起伏度、河流4个因子对农村居民点分布格局的影响。另外，现有研究多选择网格法、景观格局指数、GIS空间分析等方法。如徐羽等[6]利用网格法刻画了江西省农村居民点的时空特征，并用地理探测器测度了地形地貌对农村居民点分布的影响。姜转芳等[7]基于GIS分析得出甘肃河西地区的居民点集中分布在海拔1 200～2 400 m、坡度<5°的走廊平原和盆地以及坡度较小的缓坡地。综合当前已有研究发现，当前在流域尺度下对流域地貌与农村居民点分布关联性的研究成果较少，并且已有研究大多分析农村居民点随着各地形因子等级的增加而变化的规律，较少定量分析比较各类地形因子对农村居民点分布的相关性强度大小。根据分形地貌学理论，分形维数能够解释水系结构的拓扑关系，并揭示流域地貌侵蚀发育阶段和规律[8]。因此，通过以分形维数作为特征值来表征内外营力作用下流域地貌形态，并通过地理探测器进一步分析地形因子对农村居民点的影响特征，则可在宏观层面和中微观层面上定量计算流域地貌形态与农村居民点分布的关联性。

岷江上游流域由于岷江水系长期的深切割作用，在低海拔区域形成以“Ⅴ”型河谷为主要特点的流域地貌特征[9]。同时青藏高原东缘的隆升加强了岷江水系主河道及其支流的侵蚀作用，进一步形成沟谷纵横的河流侵蚀地貌[10-11]。另外，沿岷江及其支流河谷地带为研究区人口密集区，且自然灾害发生率较高[12-13]。研究该地区农村居民点空间分布与流域自然环境的空间关联性，有利于探索青藏高原东缘农村人居生态单元的特征及模式，对地区农村居民点布局安全优化具有一定意义。因此根据分形理论，通过水系分维数来描述岷江上游流域地貌形态，计算岷江上游流域地貌环境与农村居民点分布的数量关系，并基于地理探测器进一步探究流域地貌各地形因子与农村居民点分布的相关性，揭示流域地貌环境对农村居民点分布的影响机制，以期为实现多民族聚居山区乡村振兴的可持续发展提供一定的参考。

1 研究区概况

岷江上游地处四川盆地丘陵山地向川西北高原的过渡地带，属青藏高原东缘高山峡谷区。流域平均海拔3 440 m，流域干流全长337 km，流域面积约为2.12万km2(图1)。

图1 研究区位及范围

流域范围包括四川省的汶川、理县、黑水、茂县、松潘5个县。岷江上游地区自古以来民族迁移频繁，被称为“民族迁移走廊”。在多民族不断迁徙与融合中，最终形成羌-藏-汉-回多民族聚居区。

2 数据来源与研究方法

2.1 数据来源与预处理

研究采用30 m空间分辨率的数字高程模型(DEM)，数据来源于中国科学院计算机网络信息中心地理空间数据云镜像网站(http:∥www.gscloud.cn/)。青藏高原边界数据来源于国家青藏高原科学数据中心(http:∥data.tpdc.ac.cn/)。1∶250 000居民点数据、行政区界线和道路等数据来源于国家基础地理信息中心(http:∥www.ngcc.cn/)，时间为2017年。主要对数据进行如下预处理：

(1)由于国家基础地理信息中心下载的数据是分幅下载，空间存储单元为15′(经差)×10′(纬差)，因此通过arcpy工具包按照研究区范围对数据进行合并裁剪。

(2)对DEM数据进行预处理，通过模拟地表径流在地表的流动来产生水系，并采用Strahler河网分级方法确定沟道等级，同时以分水岭为界将岷江上游流域沟道水系划分为8个二级支流流域和97个小流域。

2.2 研究方法

2.2.1分形维数计算

提取水系后采取盒维数法测算分形维数。对于测度为r的方格网，假定测量对象占据的边长r的网格数(非空网格数)为Nr。那么，测量对象与测度之间服从负幂律关系。

Nr∝r-d。

(1)

对式(1)取对数可得：

lgNr=-dlgr+c。

(2)

式(1)～(2)中,d为标度指数；c为常数。如果测度r与非空盒子数Nr在双对数坐标轴上呈线性关系，则该物体满足分形规律。对r与Nr采取最小二乘法进行线性回归，回归直线斜率D的绝对值则为该物体的分形维数[14]。

(3)

2.2.2农村居民点分布影响因素探测

运用地理探测器分析农村居民点分布与与驱动因子之间的空间关联关系。地理探测器主要由风险探测、因子探测、生态探测和交互探测组成。模型的公式表达如下：

(4)

式(4)中，q为某指标的空间异质性，q∈[0,1]；N为研究区全部样本数；σ2为指标的方差；i为分区，i=1，2，…，L；L为分区数目。q的大小反映了空间分异的程度，q值越大，表示空间分层异质性越强，反之则空间分布的随机性越强。当q=0时，指示研究对象不存在空间异质性；当q=1时，指示完美的空间异质性[15]。

3 结果分析

3.1 流域分形地貌与农村居民点分布关联特征

已有研究表明，按照计盒法进行水系分维数计算时，计算结果会受到网格边长r的变化倍数影响。当r变化倍数大于2时出现计算波动，当变化倍数小于1.4时，计算值则非常稳定。因此，为了保证计算的可靠性，网格的变化倍数取1.2。将200 m作为网格边长r的基数，以1.2作为网格边长的变化倍数，获取20组r与Nr的数列。

通过lnNr与lnr双对数线性关系拟合回归结果来看，两者相关系数达0.98(P<0.005)，相关性显著。进一步利用相关系数来确定相应的无标度区间(分形维数有意义且相关程度最好的区间)，发现网格边长r在200～3 698 m时，岷江上游总体水系形态特征具有分形性质。由图2可知，拟合方程为lnNr=-1.229lnr+17.712，回归线斜率为-1.222 9，即水系分形维数D为1.222 9。已有研究证明，水系分维值高低表征地貌发育状况，当D≤1.6时，流域地貌处于侵蚀发育阶段的幼年期，因此岷江上游流域地貌总体上处于侵蚀发育阶段的幼年期。该阶段主要表现为水系尚未充分发育，河网密度小，河流深切侵蚀剧烈，河谷多呈“V”形。核密度分析结果显示，研究区农村居民点分布呈现狭长带状分布。这种分布形态与“V”形河谷谷坡陡峭，河谷狭窄具有直接关系，由于河谷狭窄，河道几乎占据整个河谷，居民点只能集聚分布于高等级河道两侧的谷坡上。由图3可知，岷江上游流域农村居民点整体的分布形态与分形地貌呈现正相关的叶状分支特征。

岷江上游流域主要有鱼子溪、杂谷脑河、黑水河、大姓沟、寿溪5大支流。同时根据岷江主干道与断层的相对位置，将岷江上游主干道分为3段。对8个二级支流流域分形维数进行计算，相关系数均在0．98以上，满足相关性检验，说明无标度区间的选择在各二级支流流域尺度也具备较好的标度不变性，不同区域在尺度上皆存在良好的分形特征。8个二级支流流域的横向比较中，分维值由大到小依次为干流流域上游(lnNr=-1.266 6lnr+15.503)、干流流域下游(lnNr=-1.265 6lnr+15.638)、寿溪流域(lnNr=-1.257 1lnr+14.571)、大兴沟流域(lnNr=-1.251 9lnr+15.871)、干流流域中游(lnNr=-1.227 8lnr+15.520)、渔子溪流域(lnNr=-1.220 9lnr+15.086)、杂谷脑流域(lnNr=-1.219 6lnr+15.946)和黑水河流域(lnNr=-1.197 1lnr+16.199)。

r—网格边长，Nr—网格数。

结合流域地貌来看，岷江上游主干河道流向自北向南，河谷深切，支流在干流的西侧发育，这些支流向高原腹地溯源侵蚀，形成了沟谷纵横的山地侵蚀地貌。各支流发育有明显的不均衡差异。岷江上游自北向南分布的大姓沟、黑水河、杂谷脑河、鱼子溪、寿溪这5条支流地表切割逐渐加剧，而各支流流域尺度的分形维数也反映了这一特点，基本呈现自北向南逐渐增大的趋势。如表1所示，总体上分维值较高的干流流域、寿溪流域农村居民点分布也较集中。同时岷江上游流域河谷与山脊之间相对高差达3 000 m以上，流域地貌的垂直差异也影响了农村居民点的民族类型。从各支流流域的农村居民点分布来看，614个藏族农村居民点主要分布在干流流域上段、大兴沟流域、黑水河流域、杂谷脑河流域、渔子溪流域的高半山区带(>2 000～3 000 m)；307个羌族农村居民点主要分布在干流流域中段、下段以及黑水河流域、杂谷脑河流域的低半山区和干旱河谷区带(>1 200～2 000 m)；16个回族农村居民点分布在干流流域上段且嵌于藏族农村居民点内部；96个汉族农村居民点主要分布在干流流域下段的河谷地带。

从97个小流域尺度流域的分形维数来看，地貌分形维数大约处于0.86～1.11的区间，偏态系数为-0.42，峰度系数为3.31。按照自然断点法对97个小流域尺度的分形维数进行分类，其中分形维数为0.86～0.95的小流域尺度23个，分布有156个农村居民点，平均每个单元分布7个农村居民点；分形维数为>0.95～1.01的小流域尺度45个，分布有447个农村居民点，平均每个流域单元分布10个农村居民点；分形维数为>1.01～1.11的小流域尺度31个，分布了430个农村居民点，平均每个流域单元分布14个农村居民点。因此从不同梯度流域分形维数的分级来看，呈现出流域分形维数区间越高，分布的农村居民点数量越多的特点。以居民点空间坐标为对象，提取每个农村居民点所在的小流域尺度单元水系分形维数，利用空间聚类分析进一步判断小流域尺度上城乡居民点分布对流域地貌环境发育的响应规律。从空间聚类的结果来看，聚类模式为高分形维数-高密度(H-H)的城乡居民点最多，数量为211个，占4种聚类模式的51%；其次是低分形维数-低密度(L-L)聚类模式，数量为165个，占比为40%；水系分形维数与农村居民点密度分布同效应的H-H和L-L模式总占比91%；L-H和H-L聚类模式的城乡居民点较少，占比分别为2%和7%。这也表明在小流域尺度上流域分形维数与城乡居民点分布呈正相关。

图3 全局流域尺度内农村居民点核密度分析

表1 支流流域尺度农村居民点分布特征

3.2 流域地貌因子对农村居民点分布的影响

岷江上游流域水系分形维数反映了流域地貌侵蚀发育阶段，并通过水系分形维数与农村居民点的数量关系定量揭示了流域地貌发育与农村居民点分布的空间关联性。为进一步探索流域地貌形态对农村居民点分布的影响特征，以小流域单元中的居民点密度值变化为因变量，从流域侵蚀地貌有关的微观坡面因子中选取坡度、坡长、坡度变率、剖面曲率，从宏观坡面因子中选取地表粗糙度、地形起伏度、地表切割深度以及沟壑密度、海拔等作为自变量。以上指标均以小流域尺度单元进行统计，分别计算小流域内海拔、坡度、坡度变率、剖面曲率、地表粗糙度、地形起伏度、地表切割深度的均值，以及坡长和沟壑密度。如图4～8所示，采用自然断点法将指标因素的值划分为5类，并叠加每个小流域单元内农村居民点的空间分布情况。

获取数据后，基于R语言的GD包从因子探测、风险探测、交互探测3个方面对影响因素进行分析[16]。其中因子探测器计算各地形地貌环境因子的PD值排序；风险探测器用来计算某个地形地貌因子不同分类范围农村居民点分布的差异性。交互作用探测器通过对比单个指示因子PD值与双因子PD值，判断指示因子之间是独立还是交互作用，以及其交互表现为协同还是抑制作用。

图4 支流流域尺度农村居民点民族特征分布

图5 小流域水系分形维数与农村居民点密度

L—低，H—高。

通过地理探测器的分异及因子探测器模块，评估岷江上游流域地形地貌因子与农村居民点空间分布的相关性(q值)。如图8所示，沟壑密度(0.397)>海拔(0.365)>地形起伏度(0.289)>地表切割深度(0.286)>剖面曲率(0.224)>坡度(0.22)>坡长(0.219)>地表粗糙度(0.216)>坡度变率(0.214)。由此可知，岷江上游流域沟壑密度与农村居民点分布的相关性最强，其次是海拔高程。其余因子中，地形起伏度与地表切割深度较为接近且高于其他因子。剖面曲率、坡度、坡长、地表粗糙度、坡度变率等地形因子相关性较为接近。

通过风险探测模块揭示在不同分级区间各因子对农村居民点分布的显著性差异。以q值最高的沟壑密度(X9)来看，在0.633～0.764 km·km-2的沟壑密度区间内对农村居民点分布解释力最强，当沟壑密度处于最低值区间时，其对农村居民点分布的解释力明显减少。海拔在750～1 390 m的区间内对农村居民点分布的解释力最高，表明分布在干旱河谷区的汉族农村居民点对海拔的敏感性更高。另外坡长因子对农村居民点分布的解释主要集中在坡长高值区间内。坡度、剖面曲率、地表粗糙度、地形起伏度、地表切割深度、坡度变率等因子各区间的解释能力也具有一定差异，但总体上差异较小。

图7 岷江上游流域农村居民点密度与影响因素分布

X1—坡度，X2—坡长，X3—剖面曲率，X4—地表粗糙度，X5—地形起伏度，X6—地表切割深度，X7—坡度变率，X8—海拔，X9—沟壑密度。

通过地理探测器的交互作用探测模块得出，各影响因子两两之间均存在双因子增加效应，且多数因子之间交互作用后具有非线性增强效应(表2)。其中相关性最强的交互作用影响因子是沟壑密度∩海拔(0.744)，其次是沟壑密度∩地形起伏度(0.716)、沟壑密度∩地表切割深度(0.685)、海拔∩地表切割深度(0.663)、沟壑密度∩坡度变率(0.663)。其中海拔与任意因子结合的交互作用强度均高于单项作用强度之和。而沟壑密度与任意因子的交互作用强度在所有交互结果中处于较高水平。结合因子探测器和交互作用探测器的结果，表明岷江上游农村居民点的分布受地表切割深度、地形起伏度、坡度变率、沟壑密度等地表侵蚀度指标与海拔的影响较大。

表2 岷江上游农村分布影响因素交互探测结果

4 讨论

从对岷江上游流域总体尺度、支流和小流域尺度上水系分形维数与农村居民点数量关系的研究结果来看，无论是总体尺度上农村居民点的分布形态与岷江上游流域地貌处于侵蚀发育阶段幼年期的地貌特征相吻合，还是支流和小流域尺度上农村居民点密度分布、民族特征与侵蚀地貌差异的紧密联系，均显示农村居民点与流域地貌发育具有重要的关联性(图9)。这主要是由于流域地貌发育特征所提供的生产资源、产生的自然灾害以及形成的景观特征也具有明显差异，这与农村居民点生产、生活以及居民的心理需求相耦合，从而影响农村居民点选址、迁移、消亡和扩张，形成了流域地貌环境影响农村居民点分布的内在机制。

流域地貌发育阶段性特征决定了区域可开发的土地资源、交通条件以及水利等农村居民点必需的生产生活要素，如发育相对成熟的沟道才能产生面积较大的平坦河谷地，从而为居民点扩张提供实体空间。而发育处于相对初级阶段的流域由于河谷狭窄，沟谷深切，可利用的土地资源相对稀少，导致农村居民点土地生产力低，从而使得居民点规模难以扩大。同时较大河流冲刷形成的河谷地上延伸出二级沟道，在二级沟道上又延伸出支毛沟，而交通道路往往沿河流布局，在河流等级网络的基础上形成了岷江上游流域的交通道路系统[17]，从而使得沿河流沟道布局的农村居民点具有不同的交通和水利条件。这也一定程度上解释了地理探测器探测结果中农村居民点对沟壑密度较为敏感的现象。

图9 河谷地貌发育对农村居民点分布的影响机制

另外在地理探测器的探测结果中，坡度、地表切割深度等表征地表侵蚀程度的因子对农村居民点分布的影响较强，这与沟谷发育阶段性特征所产生的自然灾害有很大关系。以泥石流为例，在河谷地貌发育的初级阶段，流水侵蚀作用强烈，沟道长度较短，山坡坡度和沟床纵比降较大，这些特征都成为泥石流灾害发生的条件。但随着沟谷发育逐渐成熟，其相对切割程度和沟床纵比降也随之降低，而沟道长度的增加使得固体物质的移动距离也将增大，这些均不利于泥石流的发生，泥石流活动将逐渐趋于衰退甚至停止[14]。总的来说，发育成熟的沟道相对稳定，发生滑坡、水土流失、坍塌、泥石流等自然灾害的几率较低，为居民点发展提供了相对安全的外部条件，而发育不成熟的沟道由于自然灾害原因易导致农村居民点迁移和消亡。

岷江上游流域的农村居民点是由汉族、藏族、回族、羌族等多个民族选址聚居而形成。在聚落选址时，由于不同民族具有明显的环境心理需求倾向，居民点会在不同的地貌环境进行选址。如藏族农村居民点选址时讲究“东有路、西有树、南有水、北有山”的择址观念，汉族农村居民点在选址时强调坐北朝南、背山面水。因此岷江上游各支流流域地貌发育的差异性所形成的不同地貌景观也对农村居民点的分布产生影响。

此外，岷江上游地区人口活动也会对流域自然环境造成一定影响，如地区人工渠系的大量修建也对该地区水系结构及地貌发育造成了不同程度的影响。受限于资料收集不足，未能在该研究中展开讨论，将在后续工作中进一步开展聚落人口活动对流域自然环境的影响研究，以期为该流域人地关系研究和村庄布局优化提供更科学的参考依据。

5 结论

为探究岷江上游流域地貌环境与农村居民点分布的关联性，通过分形理论计算了岷江上游流域总体尺度、8个支流尺度、97个小流域尺度的水系分形维数，并通过分形维数与农村居民点分布的数量关系揭示了流域地貌发育对农村居民点分布影响的空间关系。同时基于地理探测器论证了农村居民点分布与流域地貌不同因子的相关性，并进一步讨论了岷江上游流域地貌环境对农村居民点分布的影响机制。得出以下结论：

(1)在总体、支流及小流域尺度上，岷江上游流域均存在显著的分形特征，并且在总体尺度上，岷江上游流域分形地貌的幼年发育阶段促使农村居民点整体的分布形态与分形地貌呈现正相关的叶状分支特征。在支流和小流域尺度上，农村居民点密度分布和民族特征均与流域地貌发育的差异具有一致性。

(2)因子探测器表明，沟壑密度是影响农村居民点分布的主要因素，其次是海拔高程、地形起伏度、地表切割深度等因子。风险探测器表明，不同分级区间各因子对农村居民点分布具有一定的差异性，其中沟壑密度、海拔高程和坡长均具有显著的区间影响差异。交互探测器表明各地形因子之间存在明显的交互增强效应。沟壑密度∩海拔高程的交互效应最高，其次是沟壑密度∩地形起伏度。

(3)流域地貌发育所形成的生产资源、自然灾害以及景观特征通过与农村居民点生产、生活以及居民的心理需求相耦合，从而影响农村居民点选址、迁移、消亡和扩张，形成了流域自然环境影响农村居民点分布的内在机制。