戈壁是指地表由砾石覆盖的一类荒漠景观，也是我国西北干旱区的典型景观类型之一^[1]。戈壁地区大规模的开发建设破坏戈壁风蚀面后，下伏可蚀性颗粒上翻，风沙活动出现并加剧^[2]，使戈壁同沙漠一样为西北和华北地区的沙尘源区^[3]。研究戈壁地表颗粒粒度特征及空间分布特征可用来评价戈壁风蚀面的空气动力学稳定性、解释沙尘源、分析戈壁造成自然灾害的原因，进而指导防沙治沙工程建设^[4-5]。地表颗粒粒度组成在很大程度上决定了土壤抗风蚀性的强弱，是影响风蚀发生与发展的关键因子^[6-12]。因此，地表颗粒粒度特征与风蚀强度的关系一直是国际风沙学界的研究热点^{[7, 9, 13-14]}。准确获取地表不同粒级颗粒的百分比是进行准确的风蚀监测与评价的首要条件^[6-9]。以往的研究主要是通过现地取样，然后对样品筛分称重计算各粒级颗粒的百分含量，费时费力，破坏原始地表，且不适于大量、长期定点监测。随着图像处理技术的发展，一些学者开始应用数字图像处理技术来计算风蚀地表颗粒的百分含量^[15-20]，这样可以不破坏原地表，实现现地非接触、长期连续监测。

吉兰泰盐湖位于乌兰布和沙漠西南缘的贺兰山与巴音乌拉山之间的冲积洪积扇上，是内蒙古重要的盐矿生产基地。由于区域气候干旱，加之盐湖外围过度的人为干扰，使盐湖周边生态环境出现整体退化，风沙灾害威胁着盐湖盐业生产和盐湖寿命。为使盐湖正常生产，国家及地方政府也颁布了一些保护性政策，开展了一些研究，同时也建立了一系列防护措施，使盐湖周边风蚀得到了一定遏制^[21]。然而，为更进一步保护盐湖生产与运营寿命，亟待了解盐湖周边风蚀地表(戈壁、沙地等)颗粒粒度特征及其空间分布，进而建立更加科学、合理的保护与治理措施。本文以吉兰泰盐湖西北方向(主害风方向)的堆积型细砾质戈壁为研究对象，基于数字图像处理技术提取地表不同粒级(> 0.84 mm、0.84 ~ 0.42 mm和 < 0.42 mm)颗粒的质量分数，并结合地统计学分析探讨颗粒的空间分布格局，以期为盐湖周边的防沙治沙工程建设提供理论基础与科学数据。

1 材料与方法 1.1 研究区概况

研究区行政区划隶属内蒙古阿拉善左旗吉兰泰镇。据吉兰泰气象站1955—2000年的气象数据显示：区域年均气温8.96 ℃，降水量107.8 mm，蒸发量2 956.8 mm；风速3.4 m/s，最大24 m/s；大风日数33.9 d，扬沙天数85.4 d，风沙流出现频率为112.9次/年。常年以东北风出现次数最多，西南风次之；10 m/s以上大风出现次数最多的是西风和西北风。因此，区域主风向为东北风和西南风，而主害风为西风与西北风^[21]。

天然植被以旱生和超旱生的荒漠植被和盐生植被为主，如：白刺(Nitraria tangutorum Bobror.)、梭梭(Haloxylon ammodendron (C.A. Mey.) Bunge.)、沙冬青(Ammopiptanthus mongolicus (Maxim.)Cheng f.)、阿拉善沙拐枣(Calligonum alaschanicum A. Los.)、刺叶柄棘豆(Oxytropis aciphylla Ledeb.)、骆驼蓬(Peganum harmala L.)、油蒿(Artemisia ordosica krasch.)和盐爪爪(Kalidium foliatum (Pall.) Moq.)等。

1.2 样地选择与采样点布设

在吉兰泰盐湖西北方向(主害风方向)距离湖边约5 km的地方设置长期监测样地(样地中心地理坐标39°47′41″ N，105°40′04″ E，海拔1 040 m)，样地面积为250 000 m²(500 m × 500 m)。地貌类型为巴音乌拉山山前堆积型戈壁，地势为西北高于东南。样地内不同区域植被生长状况差异明显，东部生长有白刺、沙冬青等，且形成了灌丛堆；而西部仅散生草本和小灌木，如骆驼蓬、刺叶柄棘豆等(图 1)。

样地选定后，以样地西南角作为样地的原点位(0, 0)，然后向正东方向布设500 m长的样线，在样线上每隔50 m设一个采样点；从南往北依次每50 m布设1条样线。此外，为满足地统计学分析中尺度分析需要，在从西向东、南向北方向居中的样线上加密采样点，以组成不同尺度的样点距离，采样点共计141个(图 2)。在布设采样点的过程中，为了保证每条样线的走向端直，采用经纬仪和100 m卷尺相结合的方式来精确布设采样点。

1.3 图像采集

数字图像的采集使用机身为Canon Eos 5D-Mark Ⅱ，镜头为Canon EF 24 ~ 105 mm f/4L ISUSM的数码相机，分辨率为5 616 × 3 744，有效像素为21 026 304 pixel。在选定样点上放置已贴标签的标尺框；将相机固定在三脚架上，调节三脚架高度；调焦使拍摄区域恰好对准标尺框，用遮光布遮蔽阳光；全自动拍摄，每个样点拍摄3张图像^[15-16]。

1.4 图像处理及颗粒大小换算

应用Erdas Imagine 9.0遥感图像处理软件对图像处理，提取图像中颗粒的基本信息(面积、周长等)。具体提取流程参考课题组及笔者前期研究结果^[15-19]。地表颗粒大小的分级参照Chepil划分风蚀颗粒的标准：不可蚀(> 0.84 mm)、半可蚀(0.84 ~ 0.42 mm)和高度可蚀颗粒(< 0.42 mm)^[7]。

用等效直径来表示颗粒大小^[22-23]，通过公式(1)和表 1中的具体参数来计算。

$ d = 2 \times {\left( {S/{\rm{ \mathsf{ π} }}} \right)^{1/2}} $

(1)

式中：d：图像中单个颗粒的等效直径(mm)，S：图像中单个颗粒的面积(mm²)。

不同粒级颗粒质量分数采用公式(2)计算：

$ {d_{\left( X \right)}} = ({d_{(X)}}所占像元数/图像总像元数) \times 100 $

(2)

式中：d_(X)：某一粒级颗粒在图像中的像元总数。

1.5 数据处理

采用SAS 9.0软件进行回归分析和经典统计学分析，GS+ 7.0软件进行变异函数分析和Moran’s I系数分析，Surfer 8.0软件进行克立格空间插值绘图。

2 结果与分析 2.1 地表颗粒经典统计特征分析

对粒径 > 0.84 mm、0.84 ~ 0.42 mm和 < 0.42 mm颗粒的质量分数进行分析，其基本统计特征值如下(表 2)。

从理论上讲，一幅地表图像中不同粒级颗粒的质量分数之和应接近于100%，但图像处理过程中，由于在区分颗粒主体与其边缘颜色灰度值时会造成一定程度颗粒面积损失，而使颗粒总和未达到100%^[20]。本研究用图像处理技术提取的3个粒级颗粒的质量分数之和为93.12%(表 2)。不同粒级颗粒质量分数的均值可反映出样区内地表不同粒级颗粒的集中趋势，> 0.84 mm、0.84 ~ 0.42 mm和 < 0.42 mm颗粒质量分数的均值分别为69.81%、14.99%和8.32%，表明采样区地表物质组成以粒径 > 0.84 mm的颗粒为主。3个粒级颗粒均为中等变异，> 0.84 mm、0.84 ~ 0.42 mm和 < 0.42 mm颗粒的变异系数分别为14.12%、40.02%和55.06%，说明地表颗粒越小，其空间分布的变异性越大，即采样区内细颗粒的空间异质性远大于粗颗粒。

2.2 地表颗粒分布特征的变异函数分析

经K-S检验发现，3个粒级颗粒的质量分数数据均呈一定的偏态分布，但考虑到其原始数据偏态值较小，可近似视其为正态分布，可进行变异函数分析^[24]。

表 3表明，3个粒级颗粒的空间变异拟合模型相关性均显著(R²：0.60 ~ 0.95)，说明地表颗粒的空间分布具有明显的空间结构特征。> 0.84 mm和0.84 ~ 0.42 mm颗粒的变异函数理论模型为高斯模型，有效变程(A₀)=${\sqrt 3 }$a，分别为632 m和691 m；< 0.42 mm颗粒的变异函数理论模型为指数模型，A₀=3a，为1 875 m。三类颗粒的A₀均大于取样尺度，说明取样尺度设计合理。> 0.84 mm颗粒的C₀/(C₀+C) > 25%，空间相关性中等，而其他二者的C₀/(C₀+C) < 25%，空间相关性强。三者的C₀/(C₀+C)和C/(C₀+C)值表明地表颗粒的空间变异是由随机因素和结构因素共同作用的结果。其中，随机因素引起的变异所占组分相对较少，且主要表现在25 m以下尺度上；而结构因素(空间自相关部分)所占组分比例相对较多，且作用尺度范围大。

2.3 地表颗粒空间自相关分析

Moran's I系数可定量描述变量在空间上的依赖关系。I值为-1 ~ 1。I > 0表示变量在空间上正相关，I < 0表示负相关，I = 0表示空间依赖性小或空间随机变异较大。其对应的横坐标表示空间自相关范围^[25-26]。

图 3表明，> 0.84 mm、0.84 ~ 0.42 mm和 < 0.42 mm颗粒在采样尺度上I值最大值与最小值分别为0.54和-0.22、0.50和-0.22、0.42和-0.18，三者差异不显著(P > 0.05)。3个粒级颗粒的I值随间隔距离增大整体上均呈线性下降趋势(R²分别为0.93、0.88、0.92)，表明研究区地表颗粒呈现一定的空间结构性。在0 ~ 230 m范围内均表现出正相关；随着距离增大，其值逐渐向负方向增长，从正空间自相关到达负空间自相关，但负相关增加幅度较小。自相关图的变化趋势表明地表不同粒级颗粒的空间分布格局均呈斑块状。这与在研究区所调查的情况相一致，采样区西部地表植被分布较少，粗颗粒较多，为一种性质的斑块；采样区东部植被覆盖度增加，粗颗粒减少而细颗粒相对增多，为另一种性质斑块。此外，自相关图的变化趋势说明了研究区地表颗粒的变异主要由结构性因素引起，即地表颗粒的空间异质性主要与成土过程及区域气候因素等有密切的关系。

2.4 戈壁地表颗粒空间分布格局分析

根据变异函数模型，采用普通克里格法空间插值，用Surfer 8.0软件绘制颗粒空间分布等值线图(图 4)。

从图 4中可知，3个粒级颗粒均呈斑块状分布。0.84 ~ 0.42 mm和 < 0.42 mm颗粒的空间分布状况较一致，> 0.84 mm颗粒分布与其他二者相反。> 0.84 mm颗粒的含量在采样区北部较南部高，西部较东部高，西北部显著高于东南部，而其他二者的高值区与低值区的分布状况正好与之相反。这主要是因为在地表同一个采样斑块中(同一图像中)，3个粒级颗粒的质量分数总和理论上约为100%，因此三者之间存在此消彼长的现象，即同一图像中粗颗粒物含量高的地方细颗粒含量则相对低。此外，等值线图中的斑块状很好地解释了上述自相关图中的I系数值的波动现象。地表颗粒的斑块状分布可能主要与地表微地形及灌木空间分布格局有关。即大风发生时，表层较轻的细颗粒物质会被气流吹蚀而发生位移，而当遇到植被时则会沉积在植被背风区域，使区域内地表颗粒的分布格局成斑块状^[27]。

此外，为了定量揭示采样区内3个粒级颗粒的空间分布及变化情况，将采样区内南北走向的11条样线，每条样线上所有样点的同一粒级颗粒质量分数平均，分析地表不同粒级颗粒从东向西的分布及变化情况；同样的方法，分析不同粒级颗粒从南向北的分布及变化情况(图 5)。从采样区西部到东部，> 0.84 mm颗粒极显著减少(P < 0.01)，而0.84 ~ 0.42 mm和 < 0.42 mm颗粒极显著增加(P < 0.01)；从南到北，> 0.84 mm颗粒增加，0.84 ~ 0.42 mm和 < 0.42 mm颗粒减少，但均未达到显著水平(P > 0.05)。

3 讨论

研究区地表物质组成中以 > 0.84 mm颗粒为主，< 0.42 mm颗粒的变异系数大于其他二者，说明颗粒越小，其空间分布变异性越大。这主要是由于地表细颗粒物质更易受到外界条件(风蚀、践踏、碾轧等)的影响，其分布格局更易发生变化而使其空间变异性增大，所占比例的稳定性也降低^{[10, 28]}。

不同粒级土壤颗粒空间异质性是由结构性因素(土壤母质类型、地形、降水、热量等)和随机性因素(人为因素与土壤微变异等)共同作用的结果^[29-30]。本研究也得到同样的结果，> 0.84 mm、0.84 ~ 0.42 mm和 < 0.42 mm颗粒的C₀/(C₀+C)分别为27.21%、20.63%和24.43%。由此说明吉兰泰堆积型戈壁地表颗粒的空间变异特征也是由结构因素(戈壁形成、地形地势、主风向及主害风向等)和随机因素(人工挖沙、樵采、车轮碾轧及牲畜践踏等)共同作用的结果。其中，自然成土过程中的相关因子(地形地势、土壤母质、土壤类型)和气候因子(主风向及主害风向等)等非人为的区域因素是导致地表颗粒空间变异的主要内在驱动力，它有利于地表颗粒属性空间变异结构性的加强和相关性的提高，尤其是在较大尺度水平上表现更明显；而人为过程(放牧、樵采、挖沙、车轮碾轧等)则是地表颗粒特性空间变异的外在影响因素，表现为较大的随机性，它可以削弱颗粒空间变异的结构性和相关性，使颗粒特性的空间分布较均一化，这种情况在小尺度水平上表现得尤为强烈。

干旱区斑块植被区表层土壤颗粒的分布格局与植被的分布格局均为斑块状分布，且在研究空间上具有明显的差异性。植被较多且覆盖面积大的区域内，植被通过改变近地表气流，使细砂和粉砂沉积在植被周围，进而影响了样区内土壤颗粒组成的分布格局^[27]。我们的研究结果也类似，样地内地表颗粒空间分布与灌丛分布均为斑块状分布，简单对比可发现二者具有较好的吻合性，因此认为植被分布状况对地表颗粒的空间分布产生了很重要的影响。研究区大风日数较多，主害风为西风与西北风，土壤风蚀强烈。但是，由于灌丛能够有效地降低风速，使细粒物质在灌丛周围出现沉积，进而使地表不同位置不同粒级颗粒的含量发生变化，进而促使地表颗粒空间分布发生变化。地表颗粒空间分布的变化进一步导致养分、水分发生变化，最终形成了以灌木为优势种的荒漠灌丛堆(“沃岛”)^[31]，随着灌丛堆的自身发育，灌丛堆内部相对于其外部对于土壤养分及细颗粒物质的富集程度均不断增强^[32-33]，即在多风地区，灌丛“沃岛”效应可能与大风的作用结合，共同促进沙化进一步发展^[34]。此外，课题组在前期的研究中对该区域白刺和沙冬青灌丛堆周边的地表颗粒进行了探讨^{[18, 20]}，认为白刺和沙冬青灌丛堆周边地表 < 0.42 mm颗粒的水平空间异质性尺度分别为8.41 m和14.91 m，结果阐明了灌丛堆对地表颗粒空间分布所产生的影响，也可以很好地解释本研究中地表不同粒级颗粒出现斑块状分布的原因。

4 结论

基于图像处理技术提取戈壁地表颗粒信息的效果较好，可实现在现地非接触、无破坏、长期性、定点监测戈壁地表风蚀情况。研究区3个不同粒级(> 0.84 mm、0.84 ~ 0.42 mm和 < 0.42 mm)颗粒的空间分布格局是由结构因素和随机因素共同作用的结果，戈壁形成过程、地势、风蚀及植被分布状况是造成颗粒分布格局出现差异的主要内在驱动力，人为过程在一定程度上影响了颗粒的再分布。地表细颗粒由于更易受到外界条件(风蚀、践踏、碾轧等)的影响，分布格局更容易发生变化而使其空间变异性增大。