藏北高原礫石粒徑空間異質(zhì)性研究

徐 濤， 于 歡， 孔 博， 邱 霞,3， 胡孟珂， 凌鵬飛

（1.成都理工大學(xué)地球科學(xué)學(xué)院，四川 成都 610059；2.中國(guó)科學(xué)院水利部成都山地與災(zāi)害研究所，四川 成都610041；3.四川省不動(dòng)產(chǎn)登記中心，四川 成都 610014)

平均海拔4500 m以上的藏北高原，是我國(guó)的重要生態(tài)安全屏障［1］。但是，近年來(lái)由于氣候變暖以及人類活動(dòng)的加劇，藏北高原地區(qū)高寒草地退化嚴(yán)重，植被覆蓋度下降，生態(tài)環(huán)境被嚴(yán)重破壞［2］，這直接影響著藏北地區(qū)農(nóng)牧民的生產(chǎn)生活,也對(duì)青藏高原乃至全國(guó)生態(tài)造成嚴(yán)重威脅［3］。因此,藏北高原生態(tài)環(huán)境恢復(fù)成為我國(guó)生態(tài)工作的重要內(nèi)容［4］。而在藏北高原生態(tài)領(lǐng)域研究中，多數(shù)學(xué)者以草原植被和土壤組分為主要研究對(duì)象［5-10］，以高原礫石作為對(duì)象的研究相對(duì)較少。一方面，礫石是各種水文和侵蝕等過(guò)程綜合作用的產(chǎn)物，是草地和土壤退化，生態(tài)系統(tǒng)惡化的一個(gè)標(biāo)志；另一方面礫石反過(guò)來(lái)又影響水文和侵蝕過(guò)程，如入滲、蒸發(fā)、徑流和水蝕等［11］。因此，本文以地表礫石粒徑大小和空間位置為研究對(duì)象，深入分析礫石分異空間格局的效應(yīng)關(guān)系，以期為藏北高原生態(tài)修復(fù)提供參考。

在對(duì)礫石粒徑空間分異的研究中，一方面在研究對(duì)象上多是針對(duì)沙丘、瀕海（河）的沙粒以及戈壁地區(qū)礫石的空間分異，如Okin等［12］在加利福利亞?wèn)|南部的馬尼克斯（Manix）研究了羽狀沙丘沙?？臻g分異現(xiàn)象；董玉祥等［13］對(duì)河北黃金海岸沙粒的垂直分布模式進(jìn)行研究；陳西慶等［14］通過(guò)對(duì)長(zhǎng)江入河口底沙粒徑的變化進(jìn)行研究，結(jié)果表明研究區(qū)河口底沙粒徑存在繼續(xù)粗化的趨勢(shì)。另一方面在研究方法上多數(shù)利用經(jīng)典統(tǒng)計(jì)學(xué)分析海拔、粒度參數(shù)與粒徑的相關(guān)關(guān)系，或利用地統(tǒng)計(jì)學(xué)中的空間變異函數(shù)模型分析礫石粒徑的空間異質(zhì)性，如Quade［15］通過(guò)研究位于莫哈維戈壁表面的礫石分布隨海拔的變化，構(gòu)建了礫石覆蓋度隨海拔變化的函數(shù)；曹曉陽(yáng)等［16］利用經(jīng)典統(tǒng)計(jì)學(xué)方法分析了噶順戈壁洪積扇區(qū)地表礫石的粒度特征；王利兵等［17］采用地統(tǒng)計(jì)學(xué)的理論和方法研究了渾善達(dá)克沙地沙粒粒徑組成的空間異質(zhì)性。

上述研究從統(tǒng)計(jì)學(xué)及空間分析等方面對(duì)礫石粒徑空間分異展開相應(yīng)研究，但是，對(duì)高原礫石粒徑空間分異以及不同因子疊加貢獻(xiàn)率的研究報(bào)道較少。地理探測(cè)器模型是一種基于空間異質(zhì)性的研究方法，可以定量探測(cè)表達(dá)某一時(shí)空現(xiàn)象的主要驅(qū)動(dòng)因子以及不同驅(qū)動(dòng)因子之間交互作用，與其他空間異質(zhì)性探測(cè)工具相比地理探測(cè)器具有更高的解釋效率［18］。因此，本文在借鑒前人研究方法的基礎(chǔ)上，加入地理探測(cè)器模型開展藏北高原礫石粒徑空間分異的多因子定量歸因，并結(jié)合回歸分析建立礫石粒徑的空間預(yù)測(cè)模型，以期為藏北高原生態(tài)修復(fù)及環(huán)境保護(hù)提供科學(xué)依據(jù)。

1 研究區(qū)與研究方法

1.1 研究區(qū)概況

由于藏北高原面積較大，且部分地區(qū)礫石分布不具代表性，為系統(tǒng)分析藏北高原礫石的空間異質(zhì)性，故研究區(qū)選擇礫石分布特征顯著的那曲地區(qū)南部地區(qū)（29°55′21″～33°53′22″N，87°22′21″～97°3′12″E），包括安多縣南部（除色務(wù)鄉(xiāng)外其余安多縣境內(nèi)區(qū)域）、班戈縣、申扎縣、雙湖縣南部（含多瑪鄉(xiāng)、措折羅瑪鎮(zhèn)、巴嶺鄉(xiāng)、協(xié)德鄉(xiāng)）以及色尼區(qū)。研究區(qū)總面積約18.31 km2，平均海拔4100 m以上，地勢(shì)四周高、中間低，氣候寒冷干燥，大部分地區(qū)年均氣溫低于0 ℃。年均降水量從東部的548.8 mm 左右，逐漸減至西部的216.3 mm，且降水主要集中在6—9月。研究區(qū)土壤以寒鈣土和草氈土為主，植被類型較為單一，90%左右均為草甸。受到自然地理?xiàng)l件的影響，研究區(qū)經(jīng)濟(jì)發(fā)展水平較差，經(jīng)濟(jì)發(fā)展以放牧業(yè)與旅游業(yè)為主，人口密度較低，僅在東部地區(qū)有少量人口分布，地廣人稀。

1.2 數(shù)據(jù)來(lái)源

礫石粒徑空間分布數(shù)據(jù)來(lái)自野外采樣與室內(nèi)反演。野外采樣點(diǎn)位于藏北高原的那曲南部地區(qū)（圖1），海拔每升高100 m 設(shè)一條樣帶，每個(gè)樣帶設(shè)4～5個(gè)樣區(qū)，每個(gè)樣區(qū)的4個(gè)邊角及中間位置設(shè)置1個(gè)0.5 m×0.5 m 的樣方（圖1d），由于區(qū)域情況復(fù)雜，部分地區(qū)受現(xiàn)實(shí)條件的限制，無(wú)法測(cè)得5個(gè)樣方，總計(jì)測(cè)得30個(gè)樣區(qū)（圖1b），145個(gè)樣方。每個(gè)樣方內(nèi)選取占主體粒級(jí)10個(gè)以內(nèi)的礫石，并用游標(biāo)卡尺對(duì)單個(gè)礫石的Feret 直徑進(jìn)行測(cè)量（圖1e），并記錄海拔、坡度、坡向、礫石數(shù)量等基本信息。室內(nèi)反演主要將145個(gè)樣方分為130個(gè)反演樣方和15個(gè)驗(yàn)證樣方，通過(guò)采用ImageJ 軟件計(jì)算130 個(gè)反演樣方中經(jīng)形貌輪廓提取后的野外采樣礫石的粒徑值，并對(duì)研究區(qū)Landsat 8 原始影像、光譜指數(shù)因子、DEM 數(shù)據(jù)等進(jìn)行了預(yù)處理并去除因子間多重共線性，分析粒徑值與上述因子間的相關(guān)性，篩選顯著變量因子，結(jié)合多元線性逐步回歸與基于PCA 的多元回歸分析構(gòu)建礫石特征參數(shù)反演模型，進(jìn)而得到礫石粒徑反演數(shù)據(jù)［19］。該數(shù)據(jù)經(jīng)15個(gè)驗(yàn)證樣方驗(yàn)證，均方根誤差為1.4，相對(duì)誤差在30%以內(nèi)的樣區(qū)數(shù)為11個(gè)，約占總樣區(qū)數(shù)的73%，平均誤差為22.75%，反演精度良好。

圖1 野外采樣與反演成果Fig.1 Field sampling and inversion results

根據(jù)Wentworth［20］提出的礫石分級(jí)理論，將本次采集的礫石粒徑分為3 種等級(jí)。其中粒徑在2～4 mm之間的礫石分為極細(xì)礫，粒徑在4～8 mm之間的分為細(xì)礫，粒徑在8～16 mm 的礫石分為中礫。本研究反演得到的礫石粒徑分布如圖1c 所示，由圖可知，研究區(qū)粒徑主要為4～8 mm的細(xì)礫，分布于研究區(qū)的中西部，占總面積的76.99%，其次為8～16 mm的中礫，主要分布于研究區(qū)東部，占總面積的19.32%，而極細(xì)粒僅占總面積的3.69%。

環(huán)境變量數(shù)據(jù)有自然因素、人文因素兩大類，共8個(gè)代理變量，數(shù)據(jù)信息如表1所示。

表1 數(shù)據(jù)來(lái)源及介紹Tab.1 Data source and introduction

1.3 研究方法

1.3.1 數(shù)據(jù)預(yù)處理方法 由于地理探測(cè)器的正確運(yùn)行需要分類數(shù)據(jù)，故對(duì)土地利用類型、土壤類型、植被類型采用原始分類，其余影響因子采用自然間斷點(diǎn)分類法進(jìn)行離散化，經(jīng)分異及因子探測(cè)器不斷調(diào)試，最終確定NDVI、年均降水分9類，夜間燈光強(qiáng)度分8類，高程分7類，人口密度分4類。

利用ArcGIS 漁網(wǎng)工具將研究區(qū)劃分為10 km×10 km 的研究單元，然后采用雙線性插值法提取每個(gè)研究單元內(nèi)礫石粒徑值Y與分類后環(huán)境變量的類型量X，從而將兩種數(shù)據(jù)的屬性進(jìn)行匹配，最后將提取的數(shù)據(jù)帶入地理探測(cè)器模型進(jìn)行計(jì)算［21］。具體流程如圖2所示。

圖2 數(shù)據(jù)處理流程Fig.2 Data processing flow

1.3.2 空間異質(zhì)性方法

（1）空間局域異質(zhì)性Moran’sI指數(shù)可以分析空間局域異質(zhì)性，即分析某區(qū)域內(nèi)礫石粒徑與周圍礫石粒徑的空間自相關(guān)程度。Moran’sI指數(shù)介于（-1，1），I值越接近于1，表示研究單元正相關(guān)性越強(qiáng)；I值越接近于-1，表示研究單元負(fù)相關(guān)性越強(qiáng)；I值趨近于0時(shí)，表示研究單元呈隨機(jī)分布［22］。

式中：I為Moran’sI指數(shù)；n為研究礫石單元個(gè)數(shù)；Yi和Yj分別為第i個(gè)礫石單元和第j個(gè)礫石單元的粒徑值；Wij為空間權(quán)重矩陣；S2為觀測(cè)值的方差;Yˉ為粒徑值的平均值［23］。

（2）空間結(jié)構(gòu)異質(zhì)性 為進(jìn)一步分析空間自相關(guān)的范圍及空間結(jié)構(gòu)異質(zhì)性，引入空間變異函數(shù)探究其空間結(jié)構(gòu)異質(zhì)性，并計(jì)算各向同性及各向異性下礫石粒徑的最佳擬合模型及特征參數(shù)值，其公式為：

式中：z(xi)和z(xi+h)分別為區(qū)域化隨機(jī)變量z在空間中xi和xi+h處的值；N(h)為抽樣間隔等于h時(shí)的點(diǎn)對(duì)數(shù)；R(h)為變異函數(shù)值［24］。

（3）空間分層異質(zhì)性 礫石的空間異質(zhì)性極易受環(huán)境條件的影響，為此本文引入地理探測(cè)器探測(cè)影響藏北高原礫石粒徑空間分異的主導(dǎo)因子。地理探測(cè)器可以將變量進(jìn)行分層，用來(lái)探究空間分層異質(zhì)性，并揭示變量間存在某種關(guān)聯(lián)的空間分析模型［25］。該模型被廣泛應(yīng)用于土地利用、區(qū)域經(jīng)濟(jì)、精準(zhǔn)扶貧等領(lǐng)域，是一種新型空間分異性分析工具［26］。在本研究中用以探測(cè)礫石粒徑的空間分異以及環(huán)境變量因子對(duì)礫石粒徑的解釋程度，用q值度量，表達(dá)式如下：

式中：N為全區(qū)的單元數(shù)；σ2為指標(biāo)的方差；h為變量的分區(qū)，h=1，2，3，…，L；L表示分區(qū)數(shù)目。q的大小反映了空間分異的程度，q值越大，表示空間分層異質(zhì)性越強(qiáng)，反之則空間分布的隨機(jī)性越強(qiáng)。

2 結(jié)果與分析

2.1 基于Moran’s I指數(shù)的空間局域異質(zhì)性

Moran’sI指數(shù)可以用來(lái)探究研究區(qū)礫石粒徑的局域異質(zhì)性，研究區(qū)全局Moran’sI指數(shù)為0.481，Z值高于2.5 且通過(guò)1%的置信檢驗(yàn)，表明研究區(qū)礫石粒徑呈顯著的正相關(guān)，空間聚集性較強(qiáng)，存在異質(zhì)性。局部Moran’sI指數(shù)結(jié)果顯示（圖3），研究區(qū)礫石粒徑聚集模式主要為高-高聚集模式、低-低聚集模式與隨機(jī)模式。在空間分布上，研究區(qū)東部礫石多呈現(xiàn)高-高聚集模式，研究區(qū)南部礫石則以隨機(jī)分布為主，在研究區(qū)西部和北部地區(qū)，隨機(jī)分布和低-低聚集呈現(xiàn)相間分布。

圖3 研究區(qū)礫石粒徑的空間聚集特征Fig.3 Spatial clustering characteristics of gravel size in the study area

2.2 基于空間變異函數(shù)的空間結(jié)構(gòu)異質(zhì)性

通過(guò)空間變異函數(shù)計(jì)算出研究區(qū)礫石粒徑的一系列變異參數(shù)，并選擇最佳空間變異函數(shù)模型進(jìn)行擬合，從而得到研究區(qū)礫石粒徑的空間變異曲線（圖4a）?？紤]到礫石粒徑在不同方向上的變化存在一定差異性，本文從各向異性的角度出發(fā)，將研究區(qū)礫石粒徑劃分為0°、22.5°、45°、67.5°、90°、112.5°、135°及157.5°共8 個(gè)方向，各個(gè)方向上的點(diǎn)線圖如圖4b所示。

圖4 礫石粒徑空間變異曲線Fig.4 Spatial variation curve of gravel size

由圖4可知，從各向同性的角度出發(fā)，總的礫石粒徑空間變異函數(shù)類型為線性，而從各向異性的角度來(lái)看，0°、22.5°、45°、157.5°方向上變異函數(shù)類型為指數(shù)模型，90°、112.5°、135°方向上變異函數(shù)類型為高斯模型，67.5°方向上為線性模型。

根據(jù)變異函數(shù)曲線，得到研究區(qū)礫石粒徑空間結(jié)構(gòu)異質(zhì)性相關(guān)特征參數(shù)如表2所示。

表2 礫石粒徑空間變異函數(shù)特征參數(shù)Tab.2 Characteristic parameters of spatial variogram of gravel size

從各向同性角度分析，塊金值C0為0.391，具有一定的塊金效應(yīng)，亦說(shuō)明隨機(jī)因素引起的空間異質(zhì)性不可忽略。偏基臺(tái)值代表結(jié)構(gòu)因素引起的空間變異，為0.708，大于塊金值，表明研究區(qū)礫石粒徑的空間異質(zhì)性主要由結(jié)構(gòu)因素引起。基臺(tái)值為1.099，代表了由隨機(jī)因素及結(jié)構(gòu)因素共同引起的空間變異，亦表明了研究區(qū)礫石粒徑的空間異質(zhì)性不可忽略。塊基比為0.356，表明研究區(qū)礫石粒徑空間異質(zhì)性的因素中，隨機(jī)因素占比0.356，結(jié)構(gòu)因素占比0.644，結(jié)構(gòu)因素占主導(dǎo)地位。變程為26.178 km，表明研究區(qū)礫石粒徑的空間自相關(guān)范圍為26.178 km。

從各向異性的角度分析，各方向特征參數(shù)值存在差異性，但在任意方向上，偏基臺(tái)值均大于塊金值，說(shuō)明研究區(qū)礫石粒徑空間異質(zhì)性均由結(jié)構(gòu)因素主導(dǎo)?？臻g變異函數(shù)變程的大小反映礫石粒徑沿某一方向變化速度的快慢，變程值越大，表明該參數(shù)在該方向上變化慢，即非勻質(zhì)性越弱，反之越強(qiáng)。由表可知，157.5°方向變程值最大，即該方向變異性最小，而90°方向變程值最小，即該方向變異性最大。故研究區(qū)礫石粒徑各向異性特征非常明顯，各向異性主軸為157.5°方向。

2.3 基于地理探測(cè)器的空間分層異質(zhì)性

因子探測(cè)器探測(cè)各因子對(duì)礫石粒徑的解釋度（q值）大小。探測(cè)結(jié)果如表3所示，通過(guò)顯著性檢驗(yàn)的因子中，按對(duì)礫石粒徑解釋度強(qiáng)弱排序?yàn)椋篨2（NDVI）＞X3（土地利用類型）＞X7（人口密度）＞X5（植被類型）＞X6（年均降水），對(duì)應(yīng)的q值分別為0.41、0.27、0.16、0.15、0.13。其中，NDVI 的因子解釋力最大（0.41），表明藏北高原礫石粒徑空間分布受NDVI的控制作用最為強(qiáng)烈，即NDVI 與礫石粒徑在空間分布上具有最強(qiáng)的一致性。次要因素為土地利用類型，因子解釋力為0.27，說(shuō)明土地利用類型也是影響礫石粒徑空間分布的重要因素。

表3 礫石粒徑因子探測(cè)結(jié)果Tab.3 Gravel size factor detection results

交互探測(cè)主要分析各因子對(duì)礫石粒徑空間分布是否存在交互作用。結(jié)果如表4所示，由表可知，除X3∩X6表現(xiàn)為非線性增強(qiáng)外，其余交互作用類型均表現(xiàn)為雙因子增強(qiáng)。其中X2（NDVI）∩X3（土地利用類型）交互作用值最強(qiáng)，解釋力達(dá)到0.56，其次為X2（NDVI）∩X7（人口密度）與X3（土地利用類型）∩X6（年均降水），解釋力均為0.47，最小的為X6（年均降水）∩X7（人口密度），解釋力為0.27。各個(gè)因子間的交互作用共同促成了研究區(qū)礫石粒徑的空間分異格局。這說(shuō)明任意兩個(gè)環(huán)境因子交互后對(duì)礫石粒徑空間分布的因子解釋力均會(huì)顯著提升,也就是說(shuō)礫石粒徑空間分布是由各維度環(huán)境因子的共同作用形成。

表4 礫石粒徑影響因子間交互作用Tab.4 Interaction between impact factors of gravel size

風(fēng)險(xiǎn)探測(cè)器探測(cè)各因子不同分層對(duì)礫石粒徑空間分布是否存在顯著性差異及各分層中礫石粒徑均值大小。探測(cè)結(jié)果如圖5所示，NDVI 中，0.79～0.92 分區(qū)中粒徑均值最大，0.1～0.18 分區(qū)粒徑均值最小，總體粒徑均值與NDVI呈正相關(guān)，除0.24～0.32與0.61～0.79 與0.79～0.92 分層外，其余區(qū)間均存在顯著性差異。在土地利用類型中，疏林地中的粒徑均值最大，鹽堿地中粒徑均值最小，絕大多數(shù)分層間均存在顯著性差異。植被類型中高山墊狀植被中的粒徑均值最大，除高山墊狀植被分層外，其余均不存在顯著性差異。年降水量中，年降水量在470.4～548.8 mm 區(qū)間內(nèi)粒徑均值最大，261.3～281.4 mm區(qū)間內(nèi)粒徑均值最小，總體兩者呈正相關(guān)。

圖5 影響因子分區(qū)粒徑大小Fig.5 Influence factor partition particle size

人口密度因子中，人口密度與礫石粒徑兩者呈現(xiàn)一定的正相關(guān)。但是，NDVI 與人口密度均呈現(xiàn)“東高西低”的分布態(tài)勢(shì)，兩者具有空間分布上的相似性，且該區(qū)人口稀少，人類活動(dòng)對(duì)礫石粒徑的影響相對(duì)有限，結(jié)合因子探測(cè)器探測(cè)結(jié)果，NDVI 對(duì)礫石粒徑的解釋力遠(yuǎn)大于人口密度，NDVI 仍為影響礫石粒徑分布的主導(dǎo)因素。NDVI 一定程度上代表了植被覆蓋度，NDVI 大，則植被覆蓋度相對(duì)較高，水土保持能力強(qiáng)，流水、風(fēng)力等外力侵蝕作用較弱，則礫石粒徑值較大；反之，礫石粒徑越大，則土壤水分的入滲能力越強(qiáng)，則植被覆蓋度越高。因此，人口分布與礫石粒徑可能存在某種假相關(guān)，即兩者空間分布的相似性可能是由于植被等外部環(huán)境引起的。

2.4 多元線性回歸、隨機(jī)森林回歸與最優(yōu)尺度回歸驗(yàn)證

由地理探測(cè)器的因子探測(cè)器及交互作用探測(cè)器可知，上述8 種影響因子均對(duì)礫石粒徑有不同程度的影響，但是DEM、土壤類型、夜間燈光強(qiáng)度P值過(guò)大，未通過(guò)0.05 的顯著性檢驗(yàn)，故將NDVI、土地利用、植被類型、年均降水、人口密度利用多元線性回歸進(jìn)行擬合，擬合結(jié)果顯示各因子的VIF值均遠(yuǎn)＜5，說(shuō)明各因子間不存在多重共線性，模型擬合調(diào)整后如圖6a所示，R2為0.426。將上述8個(gè)影響因子加入隨機(jī)森林回歸模型進(jìn)行擬合，結(jié)果如圖6b 所示，由圖可知，隨機(jī)森林回歸模型擬合R2為0.508，隨機(jī)森林的R2高于多元線性回歸。

圖6 多元線性回歸與隨機(jī)森林回歸擬合Fig.6 Multiple linear regression and random forest regression fitting

考慮到影響因素中有土地利用類型、植被類型等名義變量，利用多元線性回歸與隨機(jī)森林回歸模型進(jìn)行擬合，無(wú)法準(zhǔn)確地反映分類變量不同取值的距離，故本文引入最優(yōu)尺度回歸進(jìn)行擬合。將礫石粒徑、NDVI、年均降水、人口密度作為數(shù)字型變量，土地利用、植被類型作為名義變量輸入模型，結(jié)果顯示最優(yōu)尺度回歸調(diào)整后R2為0.491，擬合結(jié)果好于多元線性回歸，但次于隨機(jī)森林回歸，且人口密度未通過(guò)0.05的顯著性檢驗(yàn)，原因可能為人口密度與礫石粒徑存在某種假相關(guān)，考慮剔除人口密度再次進(jìn)行回歸，回歸結(jié)果如表5所示。

表5 最優(yōu)尺度回歸模型Tab.5 Optimal scale regression model

改進(jìn)后模型擬合精度由0.491 上升到0.518，通過(guò)0.01的顯著性檢驗(yàn)，擬合結(jié)果好于多元線性回歸與隨機(jī)森林回歸，且各因子均通過(guò)0.05的顯著性檢驗(yàn)，但該模型僅用于分析研究區(qū)的礫石粒徑，模型為：

式中：Y為礫石粒徑的標(biāo)準(zhǔn)分；X2、X3、X5、X6分別為NDVI、土地利用類型、植被類型、年均降水的標(biāo)準(zhǔn)分。由地理探測(cè)器篩選出的因子除人口密度外均有較顯著的水平，通過(guò)公式（4）系數(shù)絕對(duì)值大小可知，在研究區(qū)影響礫石粒徑大小的因子中，NDVI、土地利用類型起主要作用，植被類型、年均降水起次要作用，結(jié)果與地理探測(cè)器吻合。

3 討論

由于自然環(huán)境與人類活動(dòng)的影響，礫石粒徑存在空間異質(zhì)性。本研究發(fā)現(xiàn)，藏北高原礫石粒徑存在一定程度的空間分異，且這種空間分異是由結(jié)構(gòu)因素主導(dǎo)。NDVI、土地利用類型是影響研究區(qū)礫石粒徑空間分異的主要因素，年均降水、植被類型、人口密度是影響研究區(qū)礫石粒徑空間分異的次要因素。

地理探測(cè)器探測(cè)結(jié)果顯示，NDVI 與土地利用類型是影響藏北高原礫石粒徑空間分異的主要因素，這主要是由當(dāng)?shù)氐纳鷳B(tài)環(huán)境決定的（圖7）。在研究區(qū)中，降水、NDVI 呈現(xiàn)“東高西低”的態(tài)勢(shì)，植被類型、土地利用類型以高寒禾草、苔草草原，高寒嵩草、雜草類草甸為主，四者具有空間分布上的一致性。降水影響植被覆蓋度和生長(zhǎng)狀況，植被類型反映了植被的背景信息，二者均對(duì)NDVI 與土地利用類型產(chǎn)生影響。NDVI 反映了植被覆蓋度與植被生長(zhǎng)狀態(tài)，其與土地利用類型又與水流、風(fēng)力等外力侵蝕息息相關(guān)（如高覆蓋度草地與低覆蓋度草地相比，水流、風(fēng)力等外力侵蝕作用較弱）。而礫石一方面是各種水文和侵蝕等過(guò)程綜合作用的產(chǎn)物，另一方面礫石又通過(guò)入滲、改變土壤儲(chǔ)水量等反過(guò)來(lái)影響植被和土壤肥力，進(jìn)而影響NDVI 與土地利用類型。

圖7 環(huán)境變量對(duì)礫石粒徑的影響機(jī)制Fig.7 The influence mechanism of environmental variables on gravel size

綜上所述，可以從礫石粒徑角度分析草地退化趨勢(shì)，為草原牧草退化提供動(dòng)態(tài)監(jiān)測(cè)。人類活動(dòng)導(dǎo)致藏北高原天然草地植被急劇下降，中高覆蓋草地向低覆蓋草地轉(zhuǎn)變，居民和工業(yè)用地面積增加，土地利用類型發(fā)生變化，導(dǎo)致流域水土保持能力下降，土壤破壞更為嚴(yán)重。植被退化與土地荒漠化的相互作用，惡化了區(qū)域生態(tài)環(huán)境，嚴(yán)重浪費(fèi)了區(qū)域水資源。由于供水不足和缺水，灌溉草地的牲畜負(fù)重過(guò)重，草地退化、荒漠化和水土流失更加嚴(yán)重，形成了惡性循環(huán)。上述變化將會(huì)使水土保持能力下降，外力侵蝕日益嚴(yán)重，礫石數(shù)量增多，粒徑變小。但受研究資料限制，缺乏長(zhǎng)時(shí)間序列的礫石粒徑分布數(shù)據(jù)，因此無(wú)法確定礫石粒徑變化狀態(tài)，后期將在時(shí)間尺度上進(jìn)一步研究礫石粒徑的空間分異，以期為藏北高原生態(tài)環(huán)境修復(fù)提供科學(xué)的參考依據(jù)。

4 結(jié)論

（1）研究區(qū)礫石粒徑呈顯著的正相關(guān)，空間聚集性較強(qiáng)，存在異質(zhì)性；在研究區(qū)東部礫石多呈現(xiàn)高-高聚集，研究區(qū)南部礫石則為隨機(jī)分布，在研究區(qū)西部和北部地區(qū)，隨機(jī)分布和低-低聚集呈現(xiàn)相間分布。

（2）研究區(qū)礫石粒徑具有明顯的各向異性特征，且由結(jié)構(gòu)因素主導(dǎo)。從變異函數(shù)模型來(lái)看，總的變異函數(shù)模型為線性，其余方向的模型為指數(shù)、高斯、線性；從變程來(lái)看，157.5°方向變異性最小，90°方向變異性最大。

（3）因子探測(cè)器顯示NDVI、土地利用類型為主要影響因素，植被類型、年均降水與人口密度為次要影響因素；交互探測(cè)器顯示各因子存在明顯的交互增強(qiáng)效應(yīng)，NDVI∩土地利用類型的交互作用最高，年均降水∩人口密度的交互作用最??；風(fēng)險(xiǎn)探測(cè)器表明不同分級(jí)區(qū)間各因子對(duì)礫石粒徑空間分布具有一定的差異性，其中土地利用類型具有最強(qiáng)的區(qū)間影響差異。

（4）回歸分析結(jié)果顯示，最優(yōu)尺度回歸精度最高，人口密度未通過(guò)顯著性檢驗(yàn)，相關(guān)系數(shù)由大到小為：NDVI、土地利用類型、年均降水、植被類型，且均與礫石粒徑呈正相關(guān)，除年均降水與植被類型與地理探測(cè)器略有不同外，其余結(jié)果均與地理探測(cè)器結(jié)果相吻合。