藏北高原砾石粒径空间异质性研究

徐 涛， 于 欢， 孔 博， 邱 霞,3， 胡孟珂， 凌鹏飞

（1.成都理工大学地球科学学院，四川 成都 610059；2.中国科学院水利部成都山地与灾害研究所，四川 成都610041；3.四川省不动产登记中心，四川 成都 610014)

平均海拔4500 m以上的藏北高原，是我国的重要生态安全屏障［1］。但是，近年来由于气候变暖以及人类活动的加剧，藏北高原地区高寒草地退化严重，植被覆盖度下降，生态环境被严重破坏［2］，这直接影响着藏北地区农牧民的生产生活,也对青藏高原乃至全国生态造成严重威胁［3］。因此,藏北高原生态环境恢复成为我国生态工作的重要内容［4］。而在藏北高原生态领域研究中，多数学者以草原植被和土壤组分为主要研究对象［5-10］，以高原砾石作为对象的研究相对较少。一方面，砾石是各种水文和侵蚀等过程综合作用的产物，是草地和土壤退化，生态系统恶化的一个标志；另一方面砾石反过来又影响水文和侵蚀过程，如入渗、蒸发、径流和水蚀等［11］。因此，本文以地表砾石粒径大小和空间位置为研究对象，深入分析砾石分异空间格局的效应关系，以期为藏北高原生态修复提供参考。

在对砾石粒径空间分异的研究中，一方面在研究对象上多是针对沙丘、濒海（河）的沙粒以及戈壁地区砾石的空间分异，如Okin等［12］在加利福利亚东南部的马尼克斯（Manix）研究了羽状沙丘沙粒空间分异现象；董玉祥等［13］对河北黄金海岸沙粒的垂直分布模式进行研究；陈西庆等［14］通过对长江入河口底沙粒径的变化进行研究，结果表明研究区河口底沙粒径存在继续粗化的趋势。另一方面在研究方法上多数利用经典统计学分析海拔、粒度参数与粒径的相关关系，或利用地统计学中的空间变异函数模型分析砾石粒径的空间异质性，如Quade［15］通过研究位于莫哈维戈壁表面的砾石分布随海拔的变化，构建了砾石覆盖度随海拔变化的函数；曹晓阳等［16］利用经典统计学方法分析了噶顺戈壁洪积扇区地表砾石的粒度特征；王利兵等［17］采用地统计学的理论和方法研究了浑善达克沙地沙粒粒径组成的空间异质性。

上述研究从统计学及空间分析等方面对砾石粒径空间分异展开相应研究，但是，对高原砾石粒径空间分异以及不同因子叠加贡献率的研究报道较少。地理探测器模型是一种基于空间异质性的研究方法，可以定量探测表达某一时空现象的主要驱动因子以及不同驱动因子之间交互作用，与其他空间异质性探测工具相比地理探测器具有更高的解释效率［18］。因此，本文在借鉴前人研究方法的基础上，加入地理探测器模型开展藏北高原砾石粒径空间分异的多因子定量归因，并结合回归分析建立砾石粒径的空间预测模型，以期为藏北高原生态修复及环境保护提供科学依据。

1 研究区与研究方法

1.1 研究区概况

由于藏北高原面积较大，且部分地区砾石分布不具代表性，为系统分析藏北高原砾石的空间异质性，故研究区选择砾石分布特征显著的那曲地区南部地区（29°55′21″～33°53′22″N，87°22′21″～97°3′12″E），包括安多县南部（除色务乡外其余安多县境内区域）、班戈县、申扎县、双湖县南部（含多玛乡、措折罗玛镇、巴岭乡、协德乡）以及色尼区。研究区总面积约18.31 km2，平均海拔4100 m以上，地势四周高、中间低，气候寒冷干燥，大部分地区年均气温低于0 ℃。年均降水量从东部的548.8 mm 左右，逐渐减至西部的216.3 mm，且降水主要集中在6—9月。研究区土壤以寒钙土和草毡土为主，植被类型较为单一，90%左右均为草甸。受到自然地理条件的影响，研究区经济发展水平较差，经济发展以放牧业与旅游业为主，人口密度较低，仅在东部地区有少量人口分布，地广人稀。

1.2 数据来源

砾石粒径空间分布数据来自野外采样与室内反演。野外采样点位于藏北高原的那曲南部地区（图1），海拔每升高100 m 设一条样带，每个样带设4～5个样区，每个样区的4个边角及中间位置设置1个0.5 m×0.5 m 的样方（图1d），由于区域情况复杂，部分地区受现实条件的限制，无法测得5个样方，总计测得30个样区（图1b），145个样方。每个样方内选取占主体粒级10个以内的砾石，并用游标卡尺对单个砾石的Feret 直径进行测量（图1e），并记录海拔、坡度、坡向、砾石数量等基本信息。室内反演主要将145个样方分为130个反演样方和15个验证样方，通过采用ImageJ 软件计算130 个反演样方中经形貌轮廓提取后的野外采样砾石的粒径值，并对研究区Landsat 8 原始影像、光谱指数因子、DEM 数据等进行了预处理并去除因子间多重共线性，分析粒径值与上述因子间的相关性，筛选显著变量因子，结合多元线性逐步回归与基于PCA 的多元回归分析构建砾石特征参数反演模型，进而得到砾石粒径反演数据［19］。该数据经15个验证样方验证，均方根误差为1.4，相对误差在30%以内的样区数为11个，约占总样区数的73%，平均误差为22.75%，反演精度良好。

图1 野外采样与反演成果Fig.1 Field sampling and inversion results

根据Wentworth［20］提出的砾石分级理论，将本次采集的砾石粒径分为3 种等级。其中粒径在2～4 mm之间的砾石分为极细砾，粒径在4～8 mm之间的分为细砾，粒径在8～16 mm 的砾石分为中砾。本研究反演得到的砾石粒径分布如图1c 所示，由图可知，研究区粒径主要为4～8 mm的细砾，分布于研究区的中西部，占总面积的76.99%，其次为8～16 mm的中砾，主要分布于研究区东部，占总面积的19.32%，而极细粒仅占总面积的3.69%。

环境变量数据有自然因素、人文因素两大类，共8个代理变量，数据信息如表1所示。

表1 数据来源及介绍Tab.1 Data source and introduction

1.3 研究方法

1.3.1 数据预处理方法 由于地理探测器的正确运行需要分类数据，故对土地利用类型、土壤类型、植被类型采用原始分类，其余影响因子采用自然间断点分类法进行离散化，经分异及因子探测器不断调试，最终确定NDVI、年均降水分9类，夜间灯光强度分8类，高程分7类，人口密度分4类。

利用ArcGIS 渔网工具将研究区划分为10 km×10 km 的研究单元，然后采用双线性插值法提取每个研究单元内砾石粒径值Y与分类后环境变量的类型量X，从而将两种数据的属性进行匹配，最后将提取的数据带入地理探测器模型进行计算［21］。具体流程如图2所示。

图2 数据处理流程Fig.2 Data processing flow

1.3.2 空间异质性方法

（1）空间局域异质性Moran’sI指数可以分析空间局域异质性，即分析某区域内砾石粒径与周围砾石粒径的空间自相关程度。Moran’sI指数介于（-1，1），I值越接近于1，表示研究单元正相关性越强；I值越接近于-1，表示研究单元负相关性越强；I值趋近于0时，表示研究单元呈随机分布［22］。

式中：I为Moran’sI指数；n为研究砾石单元个数；Yi和Yj分别为第i个砾石单元和第j个砾石单元的粒径值；Wij为空间权重矩阵；S2为观测值的方差;Yˉ为粒径值的平均值［23］。

（2）空间结构异质性 为进一步分析空间自相关的范围及空间结构异质性，引入空间变异函数探究其空间结构异质性，并计算各向同性及各向异性下砾石粒径的最佳拟合模型及特征参数值，其公式为：

式中：z(xi)和z(xi+h)分别为区域化随机变量z在空间中xi和xi+h处的值；N(h)为抽样间隔等于h时的点对数；R(h)为变异函数值［24］。

（3）空间分层异质性 砾石的空间异质性极易受环境条件的影响，为此本文引入地理探测器探测影响藏北高原砾石粒径空间分异的主导因子。地理探测器可以将变量进行分层，用来探究空间分层异质性，并揭示变量间存在某种关联的空间分析模型［25］。该模型被广泛应用于土地利用、区域经济、精准扶贫等领域，是一种新型空间分异性分析工具［26］。在本研究中用以探测砾石粒径的空间分异以及环境变量因子对砾石粒径的解释程度，用q值度量，表达式如下：

式中：N为全区的单元数；σ2为指标的方差；h为变量的分区，h=1，2，3，…，L；L表示分区数目。q的大小反映了空间分异的程度，q值越大，表示空间分层异质性越强，反之则空间分布的随机性越强。

2 结果与分析

2.1 基于Moran’s I指数的空间局域异质性

Moran’sI指数可以用来探究研究区砾石粒径的局域异质性，研究区全局Moran’sI指数为0.481，Z值高于2.5 且通过1%的置信检验，表明研究区砾石粒径呈显著的正相关，空间聚集性较强，存在异质性。局部Moran’sI指数结果显示（图3），研究区砾石粒径聚集模式主要为高-高聚集模式、低-低聚集模式与随机模式。在空间分布上，研究区东部砾石多呈现高-高聚集模式，研究区南部砾石则以随机分布为主，在研究区西部和北部地区，随机分布和低-低聚集呈现相间分布。

图3 研究区砾石粒径的空间聚集特征Fig.3 Spatial clustering characteristics of gravel size in the study area

2.2 基于空间变异函数的空间结构异质性

通过空间变异函数计算出研究区砾石粒径的一系列变异参数，并选择最佳空间变异函数模型进行拟合，从而得到研究区砾石粒径的空间变异曲线（图4a）。考虑到砾石粒径在不同方向上的变化存在一定差异性，本文从各向异性的角度出发，将研究区砾石粒径划分为0°、22.5°、45°、67.5°、90°、112.5°、135°及157.5°共8 个方向，各个方向上的点线图如图4b所示。

图4 砾石粒径空间变异曲线Fig.4 Spatial variation curve of gravel size

由图4可知，从各向同性的角度出发，总的砾石粒径空间变异函数类型为线性，而从各向异性的角度来看，0°、22.5°、45°、157.5°方向上变异函数类型为指数模型，90°、112.5°、135°方向上变异函数类型为高斯模型，67.5°方向上为线性模型。

根据变异函数曲线，得到研究区砾石粒径空间结构异质性相关特征参数如表2所示。

表2 砾石粒径空间变异函数特征参数Tab.2 Characteristic parameters of spatial variogram of gravel size

从各向同性角度分析，块金值C0为0.391，具有一定的块金效应，亦说明随机因素引起的空间异质性不可忽略。偏基台值代表结构因素引起的空间变异，为0.708，大于块金值，表明研究区砾石粒径的空间异质性主要由结构因素引起。基台值为1.099，代表了由随机因素及结构因素共同引起的空间变异，亦表明了研究区砾石粒径的空间异质性不可忽略。块基比为0.356，表明研究区砾石粒径空间异质性的因素中，随机因素占比0.356，结构因素占比0.644，结构因素占主导地位。变程为26.178 km，表明研究区砾石粒径的空间自相关范围为26.178 km。

从各向异性的角度分析，各方向特征参数值存在差异性，但在任意方向上，偏基台值均大于块金值，说明研究区砾石粒径空间异质性均由结构因素主导。空间变异函数变程的大小反映砾石粒径沿某一方向变化速度的快慢，变程值越大，表明该参数在该方向上变化慢，即非匀质性越弱，反之越强。由表可知，157.5°方向变程值最大，即该方向变异性最小，而90°方向变程值最小，即该方向变异性最大。故研究区砾石粒径各向异性特征非常明显，各向异性主轴为157.5°方向。

2.3 基于地理探测器的空间分层异质性

因子探测器探测各因子对砾石粒径的解释度（q值）大小。探测结果如表3所示，通过显著性检验的因子中，按对砾石粒径解释度强弱排序为：X2（NDVI）＞X3（土地利用类型）＞X7（人口密度）＞X5（植被类型）＞X6（年均降水），对应的q值分别为0.41、0.27、0.16、0.15、0.13。其中，NDVI 的因子解释力最大（0.41），表明藏北高原砾石粒径空间分布受NDVI的控制作用最为强烈，即NDVI 与砾石粒径在空间分布上具有最强的一致性。次要因素为土地利用类型，因子解释力为0.27，说明土地利用类型也是影响砾石粒径空间分布的重要因素。

表3 砾石粒径因子探测结果Tab.3 Gravel size factor detection results

交互探测主要分析各因子对砾石粒径空间分布是否存在交互作用。结果如表4所示，由表可知，除X3∩X6表现为非线性增强外，其余交互作用类型均表现为双因子增强。其中X2（NDVI）∩X3（土地利用类型）交互作用值最强，解释力达到0.56，其次为X2（NDVI）∩X7（人口密度）与X3（土地利用类型）∩X6（年均降水），解释力均为0.47，最小的为X6（年均降水）∩X7（人口密度），解释力为0.27。各个因子间的交互作用共同促成了研究区砾石粒径的空间分异格局。这说明任意两个环境因子交互后对砾石粒径空间分布的因子解释力均会显著提升,也就是说砾石粒径空间分布是由各维度环境因子的共同作用形成。

表4 砾石粒径影响因子间交互作用Tab.4 Interaction between impact factors of gravel size

风险探测器探测各因子不同分层对砾石粒径空间分布是否存在显著性差异及各分层中砾石粒径均值大小。探测结果如图5所示，NDVI 中，0.79～0.92 分区中粒径均值最大，0.1～0.18 分区粒径均值最小，总体粒径均值与NDVI呈正相关，除0.24～0.32与0.61～0.79 与0.79～0.92 分层外，其余区间均存在显著性差异。在土地利用类型中，疏林地中的粒径均值最大，盐碱地中粒径均值最小，绝大多数分层间均存在显著性差异。植被类型中高山垫状植被中的粒径均值最大，除高山垫状植被分层外，其余均不存在显著性差异。年降水量中，年降水量在470.4～548.8 mm 区间内粒径均值最大，261.3～281.4 mm区间内粒径均值最小，总体两者呈正相关。

图5 影响因子分区粒径大小Fig.5 Influence factor partition particle size

人口密度因子中，人口密度与砾石粒径两者呈现一定的正相关。但是，NDVI 与人口密度均呈现“东高西低”的分布态势，两者具有空间分布上的相似性，且该区人口稀少，人类活动对砾石粒径的影响相对有限，结合因子探测器探测结果，NDVI 对砾石粒径的解释力远大于人口密度，NDVI 仍为影响砾石粒径分布的主导因素。NDVI 一定程度上代表了植被覆盖度，NDVI 大，则植被覆盖度相对较高，水土保持能力强，流水、风力等外力侵蚀作用较弱，则砾石粒径值较大；反之，砾石粒径越大，则土壤水分的入渗能力越强，则植被覆盖度越高。因此，人口分布与砾石粒径可能存在某种假相关，即两者空间分布的相似性可能是由于植被等外部环境引起的。

2.4 多元线性回归、随机森林回归与最优尺度回归验证

由地理探测器的因子探测器及交互作用探测器可知，上述8 种影响因子均对砾石粒径有不同程度的影响，但是DEM、土壤类型、夜间灯光强度P值过大，未通过0.05 的显著性检验，故将NDVI、土地利用、植被类型、年均降水、人口密度利用多元线性回归进行拟合，拟合结果显示各因子的VIF值均远＜5，说明各因子间不存在多重共线性，模型拟合调整后如图6a所示，R2为0.426。将上述8个影响因子加入随机森林回归模型进行拟合，结果如图6b 所示，由图可知，随机森林回归模型拟合R2为0.508，随机森林的R2高于多元线性回归。

图6 多元线性回归与随机森林回归拟合Fig.6 Multiple linear regression and random forest regression fitting

考虑到影响因素中有土地利用类型、植被类型等名义变量，利用多元线性回归与随机森林回归模型进行拟合，无法准确地反映分类变量不同取值的距离，故本文引入最优尺度回归进行拟合。将砾石粒径、NDVI、年均降水、人口密度作为数字型变量，土地利用、植被类型作为名义变量输入模型，结果显示最优尺度回归调整后R2为0.491，拟合结果好于多元线性回归，但次于随机森林回归，且人口密度未通过0.05的显著性检验，原因可能为人口密度与砾石粒径存在某种假相关，考虑剔除人口密度再次进行回归，回归结果如表5所示。

表5 最优尺度回归模型Tab.5 Optimal scale regression model

改进后模型拟合精度由0.491 上升到0.518，通过0.01的显著性检验，拟合结果好于多元线性回归与随机森林回归，且各因子均通过0.05的显著性检验，但该模型仅用于分析研究区的砾石粒径，模型为：

式中：Y为砾石粒径的标准分；X2、X3、X5、X6分别为NDVI、土地利用类型、植被类型、年均降水的标准分。由地理探测器筛选出的因子除人口密度外均有较显著的水平，通过公式（4）系数绝对值大小可知，在研究区影响砾石粒径大小的因子中，NDVI、土地利用类型起主要作用，植被类型、年均降水起次要作用，结果与地理探测器吻合。

3 讨论

由于自然环境与人类活动的影响，砾石粒径存在空间异质性。本研究发现，藏北高原砾石粒径存在一定程度的空间分异，且这种空间分异是由结构因素主导。NDVI、土地利用类型是影响研究区砾石粒径空间分异的主要因素，年均降水、植被类型、人口密度是影响研究区砾石粒径空间分异的次要因素。

地理探测器探测结果显示，NDVI 与土地利用类型是影响藏北高原砾石粒径空间分异的主要因素，这主要是由当地的生态环境决定的（图7）。在研究区中，降水、NDVI 呈现“东高西低”的态势，植被类型、土地利用类型以高寒禾草、苔草草原，高寒嵩草、杂草类草甸为主，四者具有空间分布上的一致性。降水影响植被覆盖度和生长状况，植被类型反映了植被的背景信息，二者均对NDVI 与土地利用类型产生影响。NDVI 反映了植被覆盖度与植被生长状态，其与土地利用类型又与水流、风力等外力侵蚀息息相关（如高覆盖度草地与低覆盖度草地相比，水流、风力等外力侵蚀作用较弱）。而砾石一方面是各种水文和侵蚀等过程综合作用的产物，另一方面砾石又通过入渗、改变土壤储水量等反过来影响植被和土壤肥力，进而影响NDVI 与土地利用类型。

图7 环境变量对砾石粒径的影响机制Fig.7 The influence mechanism of environmental variables on gravel size

综上所述，可以从砾石粒径角度分析草地退化趋势，为草原牧草退化提供动态监测。人类活动导致藏北高原天然草地植被急剧下降，中高覆盖草地向低覆盖草地转变，居民和工业用地面积增加，土地利用类型发生变化，导致流域水土保持能力下降，土壤破坏更为严重。植被退化与土地荒漠化的相互作用，恶化了区域生态环境，严重浪费了区域水资源。由于供水不足和缺水，灌溉草地的牲畜负重过重，草地退化、荒漠化和水土流失更加严重，形成了恶性循环。上述变化将会使水土保持能力下降，外力侵蚀日益严重，砾石数量增多，粒径变小。但受研究资料限制，缺乏长时间序列的砾石粒径分布数据，因此无法确定砾石粒径变化状态，后期将在时间尺度上进一步研究砾石粒径的空间分异，以期为藏北高原生态环境修复提供科学的参考依据。

4 结论

（1）研究区砾石粒径呈显著的正相关，空间聚集性较强，存在异质性；在研究区东部砾石多呈现高-高聚集，研究区南部砾石则为随机分布，在研究区西部和北部地区，随机分布和低-低聚集呈现相间分布。

（2）研究区砾石粒径具有明显的各向异性特征，且由结构因素主导。从变异函数模型来看，总的变异函数模型为线性，其余方向的模型为指数、高斯、线性；从变程来看，157.5°方向变异性最小，90°方向变异性最大。

（3）因子探测器显示NDVI、土地利用类型为主要影响因素，植被类型、年均降水与人口密度为次要影响因素；交互探测器显示各因子存在明显的交互增强效应，NDVI∩土地利用类型的交互作用最高，年均降水∩人口密度的交互作用最小；风险探测器表明不同分级区间各因子对砾石粒径空间分布具有一定的差异性，其中土地利用类型具有最强的区间影响差异。

（4）回归分析结果显示，最优尺度回归精度最高，人口密度未通过显著性检验，相关系数由大到小为：NDVI、土地利用类型、年均降水、植被类型，且均与砾石粒径呈正相关，除年均降水与植被类型与地理探测器略有不同外，其余结果均与地理探测器结果相吻合。