藏北高原砾石粒径空间异质性研究_徐涛.pdf

第40卷第2期2023年2月Vol.40No.2Feb.2023干 旱 区 研 究ARIDZONERESEARCHhttp:/DOI：10.13866/j.azr.2023.02.13藏北高原砾石粒径空间异质性研究徐涛1，于欢1，孔博2，邱霞1,3，胡孟珂1，凌鹏飞1（1.成都理工大学地球科学学院，四川 成都610059；2.中国科学院水利部成都山地与灾害研究所，四川 成都610041；3.四川省不动产登记中心，四川 成都610014)摘要：砾石是各种水文和侵蚀等过程综合作用的产物，是草地和土壤退化、生态系统恶化的一个标志，反过来这些砾石也影响到侵蚀的各个过程。研究藏北高原地表砾石的空间分异对区域生态环境恢复具有重要意义。本文以地表砾石粒径大小和空间位置为研究对象，通过Moran s I指数、空间变异函数、地理探测器、回归分析等方法对其空间异质性进行系统性分析。结果表明：（1）全局Moran s I指数值为0.481，呈显著的正相关，局部Moran s I指数显示，研究区东部砾石聚集模式为高-高聚集，中部为低-低聚集，其余区域多呈随机分布。（2）砾石空间异质性由结构因素主导，但是变异函数最佳拟合模型与特征参数值均存在一定差异性，即存在一定的各向异性特征。（3）地理探测器结果显示，NDVI、土地利用类型为影响研究区砾石粒径空间异质性的主要因素，人口密度、植被类型、年均降水为次要因素。（4）回归分析结果显示，最优尺度回归为最佳回归模型，NDVI对砾石粒径影响最大，其后依次为土地利用类型、年均降水、植被类型。关键词：藏北高原；砾石粒径；空间异质性；地理探测器；回归分析平均海拔4500 m以上的藏北高原，是我国的重要生态安全屏障1。但是，近年来由于气候变暖以及人类活动的加剧，藏北高原地区高寒草地退化严重，植被覆盖度下降，生态环境被严重破坏2，这直接影响着藏北地区农牧民的生产生活,也对青藏高原乃至全国生态造成严重威胁3。因此,藏北高原生态环境恢复成为我国生态工作的重要内容4。而在藏北高原生态领域研究中，多数学者以草原植被和土壤组分为主要研究对象5-10，以高原砾石作为对象的研究相对较少。一方面，砾石是各种水文和侵蚀等过程综合作用的产物，是草地和土壤退化，生态系统恶化的一个标志；另一方面砾石反过来又影响水文和侵蚀过程，如入渗、蒸发、径流和水蚀等11。因此，本文以地表砾石粒径大小和空间位置为研究对象，深入分析砾石分异空间格局的效应关系，以期为藏北高原生态修复提供参考。在对砾石粒径空间分异的研究中，一方面在研究对象上多是针对沙丘、濒海（河）的沙粒以及戈壁地区砾石的空间分异，如Okin等12在加利福利亚东南部的马尼克斯（Manix）研究了羽状沙丘沙粒空间分异现象；董玉祥等13对河北黄金海岸沙粒的垂直分布模式进行研究；陈西庆等14通过对长江入河口底沙粒径的变化进行研究，结果表明研究区河口底沙粒径存在继续粗化的趋势。另一方面在研究方法上多数利用经典统计学分析海拔、粒度参数与粒径的相关关系，或利用地统计学中的空间变异函数模型分析砾石粒径的空间异质性，如Quade15通过研究位于莫哈维戈壁表面的砾石分布随海拔的变化，构建了砾石覆盖度随海拔变化的函数；曹晓阳等16利用经典统计学方法分析了噶顺戈壁洪积扇区地表砾石的粒度特征；王利兵等17采用地统计学的理论和方法研究了浑善达克沙地沙粒粒径组成的空间异质性。上述研究从统计学及空间分析等方面对砾石粒径空间分异展开相应研究，但是，对高原砾石粒径空间分异以及不同因子叠加贡献率的研究报道收稿日期：2022-07-09；修订日期：2022-11-27基金项目：国家自然科学基金项目（41971226，41871357）；中国科学院战略性先导科技专项（XDA19030303，XDA28110503）；国家重点基础研究发展规划项目：矿山环境地质灾害协同监测预警技术与装备（2017YFC1503103）作者简介：徐涛（1996-），男，硕士研究生，生态地理信息系统.E-mail:通讯作者：于欢.E-mail:292302页2期徐涛等：藏北高原砾石粒径空间异质性研究较少。地理探测器模型是一种基于空间异质性的研究方法，可以定量探测表达某一时空现象的主要驱动因子以及不同驱动因子之间交互作用，与其他空间异质性探测工具相比地理探测器具有更高的解释效率18。因此，本文在借鉴前人研究方法的基础上，加入地理探测器模型开展藏北高原砾石粒径空间分异的多因子定量归因，并结合回归分析建立砾石粒径的空间预测模型，以期为藏北高原生态修复及环境保护提供科学依据。1研究区与研究方法1.1 研究区概况由于藏北高原面积较大，且部分地区砾石分布不具代表性，为系统分析藏北高原砾石的空间异质性，故研究区选择砾石分布特征显著的那曲地区南部地区（295521335322N，87222197312E），包括安多县南部（除色务乡外其余安多县境内区域）、班戈县、申扎县、双湖县南部（含多玛乡、措折罗玛镇、巴岭乡、协德乡）以及色尼区。研究区总面积约18.31 km2，平均海拔4100 m以上，地势四周高、中间低，气候寒冷干燥，大部分地区年均气温低于0。年均降水量从东部的548.8 mm左右，逐渐减至西部的216.3 mm，且降水主要集中在69月。研究区土壤以寒钙土和草毡土为主，植被类型较为单一，90%左右均为草甸。受到自然地理条件的影响，研究区经济发展水平较差，经济发展以放牧业与旅游业为主，人口密度较低，仅在东部地区有少量人口分布，地广人稀。1.2 数据来源砾石粒径空间分布数据来自野外采样与室内反演。野外采样点位于藏北高原的那曲南部地区（图1），海拔每升高100 m设一条样带，每个样带设45个样区，每个样区的4个边角及中间位置设置1个0.5 m0.5 m的样方（图1d），由于区域情况复杂，部分地区受现实条件的限制，无法测得5个样方，总计测得30个样区（图1b），145个样方。每个样方内选取占主体粒级10个以内的砾石，并用游标卡尺对单个砾石的Feret直径进行测量（图1e），并记录海拔、坡度、坡向、砾石数量等基本信息。室内反演主要将145个样方分为130个反演样方和15个验证样方，通过采用ImageJ软件计算130个反演样方中经形貌轮廓提取后的野外采样砾石的粒径值，并对研究区Landsat 8原始影像、光谱指数因子、DEM数据等进行了预处理并去除因子间多重共线性，分析粒注：底图采用自然资源部标准地图制作，审图号为GS（2019）3333号，对底图边界无修改。下同。图1 野外采样与反演成果Fig.1 Field sampling and inversion results29340卷干旱区研究径值与上述因子间的相关性，筛选显著变量因子，结合多元线性逐步回归与基于PCA的多元回归分析构建砾石特征参数反演模型，进而得到砾石粒径反演数据19。该数据经15个验证样方验证，均方根误差为1.4，相对误差在30%以内的样区数为11个，约占总样区数的73%，平均误差为22.75%，反演精度良好。根据Wentworth20提出的砾石分级理论，将本次采集的砾石粒径分为 3 种等级。其中粒径在 24mm之间的砾石分为极细砾，粒径在48 mm之间的分为细砾，粒径在816 mm的砾石分为中砾。本研究反演得到的砾石粒径分布如图1c所示，由图可知，研究区粒径主要为48 mm的细砾，分布于研究区的中西部，占总面积的76.99%，其次为816 mm的中砾，主要分布于研究区东部，占总面积的19.32%，而极细粒仅占总面积的3.69%。环境变量数据有自然因素、人文因素两大类，共8个代理变量，数据信息如表1所示。1.3 研究方法1.3.1 数据预处理方法由于地理探测器的正确运行需要分类数据，故对土地利用类型、土壤类型、植被类型采用原始分类，其余影响因子采用自然间断点分类法进行离散化，经分异及因子探测器不断调试，最终确定NDVI、年均降水分9类，夜间灯光强度分8类，高程分7类，人口密度分4类。利用ArcGIS渔网工具将研究区划分为10 km10 km的研究单元，然后采用双线性插值法提取每个研究单元内砾石粒径值Y与分类后环境变量的类型量X，从而将两种数据的属性进行匹配，最后将提取的数据带入地理探测器模型进行计算21。具体流程如图2所示。1.3.2 空间异质性方法（1）空间局域异质性 Moran s I指数可以分析空间局域异质性，即分析某区域内砾石粒径与周围砾石粒径的空间自相关程度。Moran s I指数介于（-1，1），I值越接近于1，表示研究单元正相关性越强；I值越接近于-1，表示研究单元负相关性越强；I值趋近于0时，表示研究单元呈随机分布22。I=i=1nj=1nWij(Yi-Y)(Yj-Y)S2i=1nj=1nWij（1）式中：I为Moran s I指数；n为研究砾石单元个数；Yi和Yj分别为第i个砾石单元和第j个砾石单元的粒径值；Wij为空间权重矩阵；S2为观测值的方差;Y为粒径值的平均值23。（2）空间结构异质性 为进一步分析空间自相关的范围及空间结构异质性，引入空间变异函数探究其空间结构异质性，并计算各向同性及各向异性下砾石粒径的最佳拟合模型及特征参数值，其公式为：R(h)=12N(h)i=1N(h)z(xi)-z(xi+h)2（2）式中：z(xi)和z(xi+h)分别为区域化随机变量z在空间中xi和xi+h处的值；N(h)为抽样间隔等于h时的点对数；R(h)为变异函数值24。（3）空间分层异质性 砾石的空间异质性极易受环境条件的影响，为此本文引入地理探测器探测影响藏北高原砾石粒径空间分异的主导因子。地理探测器可以将变量进行分层，用来探究空间分层异质性，并揭示变量间存在某种关联的空间分析模型25。该模型被广泛应用于土地利用、区域经济、精准扶贫等领域，是一种新型空间分异性分析工表1 数据来源及介绍Tab.1 Data source and introduction数据类型自然因素人文因素数据名称DEM（X1）年NDVI（X2）土地利用类型（X3）土壤类型（X4）植被类型（X5）年均降水（X6）人口密度（X7）夜间灯光强度（X8）数据来源国家地球系统科学数据中心（http:）中国科学院资源环境科学与数据中心（http:）国家地球系统科学数据中心（http:）国家地球系统科学数据中心（http:）中国科学院资源环境科学与数据中心（http:）中国科学院资源环境科学与数据中心（http:）中国科学院资源环境科学与数据中心（http:）美国国家地球物理数据中心(https:www.ngdc.noaa.gov)分辨率30 m1000 m30 m1000 m1000 m500 m100 m500 m2942期徐涛等：藏北高原砾石粒径空间异质性研究具26。在本研究中用以探测砾石粒径的空间分异以及环境变量因子对砾石粒径的解释程度，用q值度量，表达式如下：q=1-Nh2hN2（3）式中：N为全区的单元数；2为指标的方差；h为变量的分区，h=1，2，3，L；L表示分区数目。q的大小反映了空间分异的程度，q值越大，表示空间分层异质性越强，反之则空间分布的随机性越强。2结果与分析2.1 基于Moran s I指数的空间局域异质性Moran s I指数可以用来探究研究区砾石粒径的局域异质性，研究区全局Moran s I指数为0.481，Z值高于2.5且通过1%的置信检验，表明研究区砾石粒径呈显著的正相关，空间聚集性较强，存在异质性。局部Moran s I指数结果显示（图3），研究区砾石粒径聚集模式主要为高-高聚集模式、低-低聚集模式与随机模式。在空间分布上，研究区东部砾石多呈现高-高聚集模式，研究区南部砾石则以随机分布为主，在研究区西部和北部地区，随机分布和低-低聚集呈现相间分布。2.2 基于空间变异函数的空间结构异质性通过空间变异函数计算出研究区砾石粒径的一系列变异参数，并选择最佳空间变异函数模型进行拟合，从而得到研究区砾石粒径的空间变异曲线（图4a）。考虑到砾石粒径在不同方向上的变化存在一定差异性，本文从各向异性的角度出发，