阮曉晗，白一茹，王幼奇，高小龍

（1.寧夏大學(xué)地理科學(xué)與規(guī)劃學(xué)院，銀川 750021；2.寧夏大學(xué)生態(tài)環(huán)境學(xué)院，銀川 750021）

礫石是粒徑≥2 mm 的礦物質(zhì)顆粒，在土壤中廣泛分布[1]。研究表明，礫石類型、粒徑大小、礫石配比、形狀、含量及其在土壤基質(zhì)中的分布等能夠顯著影響土壤物理性質(zhì)[2]、水分運(yùn)動(dòng)[3]、碳儲(chǔ)存[4]及養(yǎng)分循環(huán)[5]，并且能制約地表結(jié)皮[6]及侵蝕[7]等一系列土壤過程。由于區(qū)位、氣候和管理措施等不同，土體中礫石的空間分布在自然及人為作用的綜合影響下存在明顯的異質(zhì)性，因此探究土體中礫石的空間分布特征能夠?yàn)楹侠砝猛恋刭Y源、明確灌溉定額及恢復(fù)區(qū)域生態(tài)環(huán)境提供理論依據(jù)。

近些年已有不少學(xué)者對(duì)土體中礫石的空間分布進(jìn)行了研究，如Simanton 等[8]在美國亞利桑那州東南部的坡面研究發(fā)現(xiàn)，坡度梯度與地表巖石碎片含量顯著相關(guān)（P＜0.01）。Poesen 等[9]對(duì)地中海半干旱丘陵和橫斷面的地表礫石進(jìn)行了研究，發(fā)現(xiàn)礫石粒徑大小和覆蓋的空間變化主要受丘陵坡度梯度控制，同時(shí)坡向與坡位等也會(huì)顯著影響礫石的空間分布（P＜0.05）。Govers 等[10]在西班牙東南部半干旱丘陵地區(qū)研究表明，土壤表面覆蓋的礫石量不僅與坡度梯度有關(guān)，同時(shí)與曲率線性相關(guān)。朱元駿等[11]在黃土高原水蝕風(fēng)蝕交錯(cuò)帶研究了坡面表土礫石覆蓋度和粒徑分布特征，發(fā)現(xiàn)坡面礫石分布是侵蝕與坡度因素疊加作用的結(jié)果。馮雪瑾等[12]測定了太行山南麓地區(qū)不同坡位、不同土層深度的礫石含量及礫石粒徑的分布特征，結(jié)果表明不同粒徑礫石體積分?jǐn)?shù)均呈現(xiàn)從坡頂?shù)狡履_逐漸減小的趨勢(shì)。Nyssen 等[13]的研究表明，耕地中礫石的空間分布受地形地貌條件、水文過程及人為作用共同影響。綜上可知，在不同區(qū)域和尺度條件下礫石分布均存在明顯空間異質(zhì)性而主要影響因素有所不同，當(dāng)前研究大多基于不同區(qū)域自然因素主導(dǎo)的礫石空間異質(zhì)性研究，針對(duì)人為影響較大的土體含礫石區(qū)域研究較少，而壓砂地等人為形成的特殊礫石-土壤聚集體中礫石的空間異質(zhì)性及其影響因素研究相對(duì)缺乏。

壓砂地是我國西北地區(qū)充分利用當(dāng)?shù)氐[石資源，應(yīng)對(duì)干旱氣候的一種典型耕作模式[14]，其表面覆蓋的礫石能夠有效促進(jìn)土壤水分入滲、增溫保墑[15]。然而隨著種植年限增加，壓砂地的礫石與土壤層混合、連通，導(dǎo)致其土壤理化性質(zhì)[16]、保水保墑功能變差[17]，干燥化程度加劇[16]。當(dāng)前壓砂地出現(xiàn)大量經(jīng)長期種植后撂荒的現(xiàn)象，約有1/3退化嚴(yán)重，近年來寧夏各地區(qū)逐步開展壓砂地有序退出及區(qū)域生態(tài)修復(fù)等工作，因此明確壓砂地不同粒徑礫石空間變異特征，并進(jìn)一步分析其影響因素是研究壓砂地退化機(jī)制的基礎(chǔ)問題。為此本研究對(duì)壓砂地不同粒徑礫石進(jìn)行統(tǒng)計(jì)分析，采用空間自相關(guān)和半方差函數(shù)的分析方法，從不同角度揭示研究區(qū)壓砂地不同粒徑礫石空間結(jié)構(gòu)特征和分異規(guī)律，并通過地理探測器深入探討種植年限、坡向和地表粗糙度等因子對(duì)不同粒徑礫石空間分異的影響程度，選取影響因子較高的環(huán)境變量協(xié)同地理加權(quán)回歸克里格法進(jìn)行插值，獲取了不同粒徑礫石空間分布特征，以期為壓砂地逐步退出、修復(fù)重建和分類施策等提供理論支撐。

1 材料與方法

1.1 研究區(qū)概況

研究區(qū)位于中衛(wèi)市沙坡頭區(qū)興仁鎮(zhèn)紅溝（36°56′32″～36°57′45″N，105°22′22″～105°21′29″E），地處寧夏中部干旱帶，屬干旱半干旱氣候，年平均溫度6.9 ℃，≥10 ℃積溫2 320 ℃，無霜期163 d 左右。多年平均降水量不足200 mm，年蒸發(fā)量2 390 mm左右，全年日照3 800 h 左右，光熱資源十分豐富，晝夜溫差大，土壤為灰鈣土，是硒砂瓜主產(chǎn)區(qū)。覆蓋礫石來自香山風(fēng)化碎石，覆蓋厚度為20 cm 左右，利用ArcGIS 10.6制作樣點(diǎn)分布及高程圖（圖1）。

圖1 研究區(qū)樣點(diǎn)分布及高程Figure 1 Sampling point distribution and elevation of the study area

1.2 采樣方法

2020 年9 月20 日至10 月5 日進(jìn)行不同年限壓砂地礫石分選。根據(jù)前期實(shí)地調(diào)研，明確了研究區(qū)范圍內(nèi)所有壓砂地的種植年限，并將所有壓砂地的種植年限劃分為：2、5、10、20、30 a 和40 a。研究區(qū)總面積為2.54 km2，種植2、5、10、20、30 a 和40 a 的壓砂地面積占比分別為10.75%、30.42%、32.37%、14.69%、5.92%和5.85%。在研究區(qū)壓砂地集中分布區(qū)域根據(jù)各年限壓砂地所占面積比，同時(shí)結(jié)合網(wǎng)格法均勻布點(diǎn)，樣點(diǎn)數(shù)分別為11、31、33、15、6 和6，共計(jì)102 個(gè)。在每一個(gè)樣點(diǎn)內(nèi)選取50 cm×50 cm 范圍內(nèi)礫石覆蓋層全部砂石，風(fēng)干后用50、31.5、25、16、10、5、2 mm 篩子逐層過篩，選取不同粒徑的礫石測定質(zhì)量并記錄，壓砂地礫石粒徑配比情況見表1。

表1 壓砂地不同粒徑礫石配比（%）Table 1 Proportion of gravel with different particle sizes in gravel-sand mulched fields（%）

1.3 影響因子采集與分類

礫石空間分布受自然因素與人文因素共同影響，有研究表明土體中礫石的空間分布與地形因子之間存在緊密聯(lián)系[8-10]，地形對(duì)土體內(nèi)部礫石空間變異的影響不容忽視，這些地形因子能夠通過DEM 獲取[18]。本研究以研究區(qū)DEM 影像為基礎(chǔ)，提取坡度、坡向、地表粗糙度、地表起伏度、剖面曲率及平面曲率6 種地形因子，并通過對(duì)研究區(qū)的前期調(diào)查獲取壓砂地種植年限數(shù)據(jù)，其中DEM 數(shù)據(jù)為30 m 的ASTER GDEM數(shù)據(jù)，來源于地理空間數(shù)據(jù)云網(wǎng)站。為充分考慮多因子對(duì)礫石空間分布的影響，本文選取自然因素和人文因素共7 個(gè)影響因子，并利用自然斷點(diǎn)法，對(duì)數(shù)值型驅(qū)動(dòng)因子進(jìn)行等級(jí)劃分。

1.4 統(tǒng)計(jì)分析方法

1.4.1 空間自相關(guān)分析

空間自相關(guān)是指地理對(duì)象的某一屬性值在相鄰空間位置處取值之間的相互關(guān)系，是對(duì)該屬性空間聚集或分散程度的度量[19]?？臻g自相關(guān)現(xiàn)象廣泛存在于地理空間數(shù)據(jù)中，因此，進(jìn)行空間自相關(guān)分析對(duì)深入理解地理對(duì)象的空間分布特征具有重要意義。Moran′sI指數(shù)在地理對(duì)象的空間自相關(guān)性研究中具有廣泛應(yīng)用，其不僅能在總體水平上對(duì)空間自相關(guān)現(xiàn)象做出評(píng)價(jià)，而且能夠?qū)Φ乩韺傩灾档木奂碗x散程度作進(jìn)一步的區(qū)域可視化。本研究中全局空間自相關(guān)采用全局Moran′sI指數(shù)描述：

式中：N為采樣點(diǎn)總數(shù)；Wij為空間權(quán)重；Xi和Xj分別為i處和j處的礫石含量；為礫石含量均值；Z為標(biāo)準(zhǔn)化統(tǒng)計(jì)量值；E（I）為觀測變量自相關(guān)性的期望值；Var（I）為方差。Moran′sI指數(shù)的取值范圍為[-1，1]，I＞0 代表屬性值呈現(xiàn)空間正相關(guān)，趨于空間聚集；I＜0代表屬性值呈現(xiàn)空間負(fù)相關(guān)，趨于空間分散；I=0代表屬性值趨于空間隨機(jī)分布。采用Z值進(jìn)行顯著性檢驗(yàn)，當(dāng)Z≥1.96或Z≤-1.96時(shí)，表示空間相關(guān)性顯著（P＜0.05）。

1.4.2 空間變異性分析

采用地統(tǒng)計(jì)學(xué)分析壓砂地不同粒徑礫石占比的空間變異特征，其公式為：

式中：γ（h）代表半方差函數(shù)值；N（h）代表相距為h的點(diǎn)對(duì)數(shù)；Z（xi）代表x=xi處變量Z的實(shí)測值；Z（xi+h）代表x=xi+h處變量Z的實(shí)測值。以滯后距h為橫軸、半方差函數(shù)值為縱軸可以繪制半方差圖。根據(jù)半方差圖可以得到塊金值（C0）、基臺(tái)值（C0+C）、變程（A）三個(gè)重要特征值。C0表示隨機(jī)部分引起的空間變異性；C0+C表示變量的最大變異程度；A表示變程，即變量自相關(guān)范圍。通過公式C0/（C0+C）可以計(jì)算得到一個(gè)描述空間變量依賴性的重要指標(biāo)，即空間異質(zhì)比，根據(jù)Cambardella 等[20]的劃分標(biāo)準(zhǔn)，C0/（C0+C）≤25%表示強(qiáng)的空間依賴性，25%＜C0/（C0+C）＜75%表示中等空間依賴性，C0/（C0+C）≥75%表示弱的空間依賴性。

1.4.3 地理探測器

地理探測器是探測空間分異并揭示某種地理屬性與其解釋因子之間關(guān)系的一種空間分析模型[21]，能夠探測單因子對(duì)因變量空間分異的影響，并對(duì)其顯著性進(jìn)行統(tǒng)計(jì)學(xué)檢驗(yàn)。本研究采用因子探測功能來分析各影響因子對(duì)壓砂地土體中不同粒徑礫石空間分異的影響。其計(jì)算公式[22]如下：

式中：Nh為子類型區(qū)h的單元數(shù)；L為變量的分層；N為全區(qū)的單元數(shù)；為不同粒徑礫石占比在子類型h區(qū)的方差；σ2為不同粒徑礫石占比在全區(qū)的方差。q值為某因子對(duì)不同粒徑礫石空間分異的解釋程度，其取值范圍為[0，1]，值越大表示該因子對(duì)不同粒徑礫石空間分異的解釋力越強(qiáng)，反之則越弱。

1.4.4 地理加權(quán)回歸克里格法

地理加權(quán)回歸克里格法（GWRK）是通過對(duì)目標(biāo)變量和輔助變量進(jìn)行地理加權(quán)回歸擬合（GWR），得到局部回歸的殘差項(xiàng)，然后使用普通克里格插值法（OK）對(duì)所得殘差項(xiàng)進(jìn)行插值[23]。GWR模型結(jié)構(gòu)如下：

式中：YGWR（s0）為因變量Y在s0處的模擬值；β0（s0）為截距項(xiàng)；βi（s0）為第k個(gè)協(xié)變量的回歸系數(shù)；Xk（s0）為第i個(gè)樣點(diǎn)上第k個(gè)變量的測量值；s0為局部估計(jì)的系數(shù)，用加權(quán)最小二乘法求得，是空間位置的函數(shù)，其權(quán)重矩陣的估計(jì)采用高斯函數(shù)；k為自變量數(shù)。

GWRK 法是用OK 法對(duì)GWR 模型擬合后的殘差進(jìn)行插值，該方法是對(duì)GWR法的一種深入推廣，并增加了GWR法擬合的趨勢(shì)[24]，其表達(dá)式為：

式中：ε（s0）為GWR 模型在點(diǎn)s0處擬合后所得的殘差，采用OK法進(jìn)行插值。

1.4.5 精度評(píng)價(jià)

隨機(jī)生成的80 個(gè)用于模型擬合的測試點(diǎn)參與所有分析過程，其余22 個(gè)測試點(diǎn)用于驗(yàn)證空間插值結(jié)果。用驗(yàn)證點(diǎn)估計(jì)值和實(shí)測值的平均相對(duì)誤差（MRE）和均方根誤差（RMSE）評(píng)價(jià)不同方法的預(yù)測準(zhǔn)確性。

1.5 數(shù)據(jù)處理

利用SPSS 20.0 進(jìn)行經(jīng)典統(tǒng)計(jì)分析，在GS+9.0 軟件中完成半方差函數(shù)分析，運(yùn)用GeoDa 1.16軟件進(jìn)行全局空間自相關(guān)分析，在GeoDetector 中完成地理探測器的影響因子分析，在GWR 4.0中完成GWR分析，用ArcGIS 10.6 軟件分別對(duì)研究區(qū)不同粒徑礫石進(jìn)行GWRK 插值，常規(guī)統(tǒng)計(jì)分析和繪圖處理分別在Excel 2016和Origin 2017中進(jìn)行。

2 結(jié)果與討論

2.1 礫石粒徑占比的描述性統(tǒng)計(jì)

表2是壓砂地102個(gè)采樣點(diǎn)礫石覆蓋層礫石粒徑占比的描述性統(tǒng)計(jì)特征值，各樣點(diǎn)礫石粒徑＞50、31.5～50、25～31.5、16～25、10～16、5～10、2～5 mm 和≤2 mm 平均占比分別為3.40%、3.98%、4.76%、7.33%、11.56%、15.09%、22.36%和31.52%。經(jīng)典統(tǒng)計(jì)學(xué)通常用變異系數(shù)（CV）的大小描述變量的變異程度，當(dāng)CV≤10%時(shí)為弱變異，當(dāng)10%＜CV＜100%時(shí)為中等變異，當(dāng)CV≥100%時(shí)為強(qiáng)變異[25]。表2結(jié)果表明礫石粒徑占比的變異程度隨粒徑的減小逐漸減小，＞50 mm礫石CV最大，為59.62%，≤2 mm 礫石CV最小，為19.82%，所有粒徑礫石占比均屬中等變異，說明不同粒徑礫石含量受隨機(jī)性因素和結(jié)構(gòu)性因素共同影響。通過K-S檢驗(yàn)可發(fā)現(xiàn)僅＞50、5～10 mm和≤2 mm粒徑礫石屬正態(tài)分布，其余5 個(gè)粒徑范圍的礫石均不符合正態(tài)分布，說明不同粒徑礫石含量受人為干擾較大。為滿足地統(tǒng)計(jì)分析要求，對(duì)不符合正態(tài)分布的粒徑進(jìn)行對(duì)數(shù)正態(tài)轉(zhuǎn)換，使其符合正態(tài)分布。

表2 不同粒徑礫石占比的描述性統(tǒng)計(jì)Table 2 Descriptive statistics of gravel proportion with different particle size

2.2 基于地統(tǒng)計(jì)學(xué)的不同粒徑礫石空間變異特征

在GS+9.0 中對(duì)壓砂地不同粒徑礫石占比進(jìn)行半方差函數(shù)模型擬合，結(jié)果如表3 所示。參照半方差函數(shù)理論模型的選取標(biāo)準(zhǔn)（決定系數(shù)R2接近于1、殘差RSS趨近于0），發(fā)現(xiàn)粒徑＞50、31.5～50、25～31.5、16～25、10～16、5～10、2～5 mm 和≤2 mm 的礫石含量適宜的模型分別為指數(shù)模型、線性模型、純塊金模型、指數(shù)模型、線性模型、高斯模型、線性模型、球型模型。礫石的空間變異性受結(jié)構(gòu)性因素和隨機(jī)性因素的共同影響，結(jié)構(gòu)性因素（礫石母質(zhì)、地形等）會(huì)導(dǎo)致空間相關(guān)性增強(qiáng)，隨機(jī)性因素（耕作制度、種植年限等）使礫石空間相關(guān)性減弱[26]。在8 個(gè)粒徑礫石中2～5、10～16 mm 的礫石空間異質(zhì)比均大于75%，說明其受隨機(jī)性因素影響劇烈；25～31.5 mm 的最佳模型為純塊金模型，說明其空間上呈隨機(jī)性分布；其余粒徑的空間異質(zhì)比均小于25%，說明結(jié)構(gòu)性因素是空間分布的主導(dǎo)因素。變程結(jié)果顯示10～16、2～5 mm 粒徑的變程均遠(yuǎn)大于平均采樣間距（250 m），說明相較于其他粒徑礫石，10～16、2～5 mm 礫石的空間連續(xù)性的尺度范圍更大。

表3 壓砂地不同粒徑礫石占比的半方差模型參數(shù)Table 3 Semivariance model parameters of gravel proportion with different particle size in gravel-sand mulched field

2.3 Moran′s I指數(shù)

圖2 為通過GeoDa 軟件獲取的不同粒徑礫石全局Moran 散點(diǎn)圖，全局Moran′sI指數(shù)為Moran 散點(diǎn)圖的斜率。對(duì)壓砂地不同粒徑礫石進(jìn)行全局空間自相關(guān)分析發(fā)現(xiàn)，礫石粒徑范圍為＞50、31.5～50、25～31.5、16～25、10～16、5～10、2～5 mm 和≤2 mm 的全局Moran′sI指數(shù)分別為0.892、0.717、0.928、0.924、0.923、0.674、0.944和0.779，各粒徑礫石的全局Moran′sI值均高于0.674，P均小于0.01，表明壓砂地不同粒徑礫石均表現(xiàn)出很強(qiáng)的空間正相關(guān)，且Z值均大于1.96，表示空間相關(guān)性顯著，說明樣點(diǎn)與周圍樣點(diǎn)的值相似。在Moran散點(diǎn)圖中，數(shù)據(jù)處在第一、第三象限表示正的空間自相關(guān)性，數(shù)據(jù)處于第二、第四象限表示負(fù)的空間自相關(guān)性。由圖2可知，各粒徑礫石的散點(diǎn)分布狀況不同，但整體均主要分布在第一、第三象限，說明其在空間分布上并非完全隨機(jī)，而是具有極顯著的空間依賴特征，即各粒徑礫石在空間上呈顯著聚集性。

圖2 不同粒徑礫石占比的空間自相關(guān)分析Figure 2 Spatial autocorrelation analysis of gravel proportion with different particle size

2.4 壓砂地不同粒徑礫石空間分異的因素識(shí)別

本研究利用地理探測器的因子探測器得到種植年限、坡度、坡向、地表粗糙度、地表起伏度、剖面曲率及平面曲率7 個(gè)可能影響壓砂地不同粒徑礫石空間分異因素的q值并對(duì)其進(jìn)行了排序（表4），各個(gè)因子的解釋力均通過了1%的顯著性水平檢驗(yàn)。對(duì)各因子的平均q值綜合分析可知，影響壓砂地不同粒徑礫石空間分異的因子按影響大小依次為種植年限＞坡向＞地表粗糙度＞坡度＞剖面曲率＞平面曲率＞地表起伏度。

表4 壓砂地不同粒徑礫石占比空間分異的因子探測結(jié)果Table 4 Factor detection results of spatial differentiation of gravel proportion with different particle size in gravel-sand mulched fields

2.5 壓砂地不同粒徑礫石的空間分布特征及精度評(píng)價(jià)

GWR 模型通過計(jì)算回歸模型的局部參數(shù)來消除地理數(shù)據(jù)中存在的空間自相關(guān)性及空間非平穩(wěn)性問題，克服了回歸模型中不同空間位置的環(huán)境因子權(quán)重一致的不足，實(shí)現(xiàn)了局部空間信息的捕獲，提高了插值精度，同時(shí)其在空間制圖方面表現(xiàn)出更豐富的局部細(xì)節(jié)和更光滑的空間漸變特征[27]，因此可以更加直觀揭示出壓砂地礫石配比的空間分布特征。

本研究在地理探測器識(shí)別壓砂地不同粒徑礫石空間分異因子的基礎(chǔ)上，選取q值最高的三個(gè)因子（種植年限、坡向和地表粗糙度）作為GWR 模型的環(huán)境變量，并通過地理加權(quán)回歸克里格法（GWRK）插值獲取研究區(qū)壓砂地礫石分布圖（圖3）?！? mm的礫石占比高值主要集中于研究區(qū)中部，西南部含量較低。2～5 mm 礫石占比高值區(qū)和低值區(qū)均呈破碎化、分散化分布。5～10 mm 和10～16 mm 礫石占比分布呈現(xiàn)相似的分布格局，高值均分布于東北部，低值集中分布于中部；16～25 mm 礫石占比高值區(qū)主要分布于西部和西南部，低值區(qū)主要集中于東北部；25～31.5 mm 礫石占比高值主要分布于西南部，東北部有零散分布，低值主要分布于中部；31.5～50 mm礫石占比的低值主要分布于東北部，其余部分均有明顯高值分布；＞50 mm礫石占比主要集中分布于中部、西部及南部。

圖3 不同粒徑礫石占比空間分布狀況Figure 3 Spatial distribution of gravel proportion with different particle size

本研究通過計(jì)算22 個(gè)驗(yàn)證點(diǎn)的MRE和RMSE來檢驗(yàn)的GWRK 的插值精度，評(píng)價(jià)結(jié)果如表5 所示。＞50、31.5～50、25～31.5、16～25、10～16、5～10、2～5 mm和≤2 mm 粒徑礫石占比的MRE和RMSE均趨近于0，說明引入種植年限、坡向和地表粗糙度三個(gè)影響因子的GWRK 模型能夠較為精準(zhǔn)地模擬壓砂地不同粒徑礫石的空間分布狀況，同時(shí)也可以驗(yàn)證種植年限、坡向和地表粗糙度為壓砂地不同粒徑礫石的主要影響因子。

表5 GWRK精度評(píng)價(jià)Table 5 Accuracy evaluation of GWRK

2.6 不同粒徑礫石占比的空間分布規(guī)律的成因分析

地統(tǒng)計(jì)分析表明不同粒徑礫石的空間異質(zhì)性受隨機(jī)性因素和結(jié)構(gòu)性因素共同影響，通過地理探測器進(jìn)一步分析發(fā)現(xiàn)，種植年限、坡向、地表粗糙度、坡度、剖面曲率、平面曲率、地表起伏度均對(duì)礫石占比的空間分布產(chǎn)生不同程度的影響。

隨種植年限增加，壓砂地礫石層中大粒徑礫石占比逐漸降低而小粒徑礫石占比逐漸提高。相關(guān)研究表明影響礫石分布的原因主要有地表徑流、風(fēng)蝕風(fēng)化及人為活動(dòng)干擾[13，28]，研究區(qū)處于干旱半干旱帶，降雨稀缺，自然條件下難以形成地表徑流，因此僅有小粒徑礫石和細(xì)土可能在農(nóng)業(yè)灌溉作用下向下遷移。風(fēng)蝕風(fēng)化作用對(duì)不同粒徑礫石的分布有顯著影響，坡向會(huì)影響太陽輻射的強(qiáng)度和時(shí)間，地表粗糙度表示地表對(duì)風(fēng)速和風(fēng)沙活動(dòng)的影響，坡度、剖面曲率、平面曲率和地表起伏度均反映了區(qū)域微地形狀況，均可影響礫石的風(fēng)蝕風(fēng)化進(jìn)程[29]，導(dǎo)致不同位置的礫石在流水、風(fēng)力及重力等外營力作用下礫石的空間分布有所差異。一些學(xué)者研究表明，坡度是礫石分布的重要影響因子[8-10]，壓砂地不同粒徑礫石空間分異的主導(dǎo)因素識(shí)別表明，本研究中坡度對(duì)不同粒徑礫石分布的作用相對(duì)較弱，這主要是由于壓砂地一般選取地表平整、地表起伏度較小的區(qū)域進(jìn)行鋪砂及后續(xù)農(nóng)業(yè)耕作，因此研究區(qū)內(nèi)各樣點(diǎn)的坡度、剖面曲率、平面曲率和地表起伏度變化較小，對(duì)于礫石分布的影響相對(duì)較低。坡向的q值高于其他幾個(gè)地形因子，這是因?yàn)檠芯繀^(qū)地貌類型為黃土丘陵，區(qū)內(nèi)溝壑縱橫，雖然壓砂地所處位置一般為坡底，但其坡向存在差異，因此研究區(qū)內(nèi)各樣點(diǎn)所處位置的坡向成為影響侵蝕和風(fēng)化的重要因子，對(duì)礫石的空間分布產(chǎn)生較大作用。種植年限的q值明顯高于其他因子，說明不同粒徑礫石的空間異質(zhì)性受種植年限影響更為顯著。這可能是由于種植年限增加，耕作過程中機(jī)械碾壓和耖砂翻耕會(huì)加速礫石由大粒徑向小粒徑轉(zhuǎn)化的過程，導(dǎo)致出現(xiàn)種植年限越大礫石破碎化越明顯，以及土體中小粒徑礫石占比提高的現(xiàn)象，譚軍利等[30]也發(fā)現(xiàn)隨著砂田種植年限延長，壓砂地大粒徑的礫石量減少，小粒徑的礫石增多。

本研究從多個(gè)角度分析了壓砂地土體中不同粒徑礫石的空間分異及其影響因素，但不同因子的具體作用仍有待進(jìn)一步探究。同時(shí)有必要針對(duì)不同種植年限壓砂地土壤理化性質(zhì)、土壤水分運(yùn)動(dòng)及溶質(zhì)運(yùn)移狀況開展進(jìn)一步研究，以期為實(shí)現(xiàn)壓砂地農(nóng)業(yè)可持續(xù)發(fā)展、恢復(fù)區(qū)域生態(tài)功能提供理論支撐。

3 結(jié)論

（1）不同粒徑礫石空間變異系數(shù)為19.82%～59.62%，屬中等變異，其空間變異性受隨機(jī)性因素和結(jié)構(gòu)性因素共同影響。

（2）不同粒徑礫石占比的空間分布狀況不同，各粒徑礫石的全局Moran′sI值均高于0.674，并全部通過極顯著性檢驗(yàn)（P＜0.01），表明在全局尺度上壓砂地不同粒徑礫石均表現(xiàn)出很強(qiáng)空間正相關(guān)，即均存在極顯著的空間依賴特征，且呈聚集性分布。

（3）壓砂地不同粒徑礫石空間分異的解釋力由大到小依次為種植年限＞坡向＞地表粗糙度＞坡度＞剖面曲率＞平面曲率＞地表起伏度，種植年限、坡向和地表粗糙度為礫石空間分異的主要影響因子。