侯 琨，王秀茹?，王計(jì)平，袁普金，程 復(fù)

(1.北京林業(yè)大學(xué)水土保持學(xué)院，100083，北京；2.中國(guó)林業(yè)科學(xué)研究院國(guó)家林業(yè)和草原局鹽堿地研究中心，100091，北京；3.水利部水土保持監(jiān)測(cè)中心，100053，北京)

由景觀要素的空間分布和組合特征構(gòu)成的景觀格局與生態(tài)過(guò)程的相互聯(lián)系是開(kāi)展景觀生態(tài)學(xué)研究的重要切入點(diǎn)[1]。格局決定過(guò)程，過(guò)程也會(huì)推進(jìn)格局的演進(jìn)。作為一種典型的面狀物質(zhì)能量流動(dòng)過(guò)程，水土流失的發(fā)生和發(fā)展與景觀格局的演變密不可分[2-3]。在人類干擾較強(qiáng)的景觀格局中，劇烈波動(dòng)變化的土地利用/覆被情況是景觀格局演進(jìn)的直接表現(xiàn)[4]，其引起的景觀結(jié)構(gòu)和功能變化，可對(duì)水土流失等生態(tài)過(guò)程產(chǎn)生重大影響[5]。

從景觀類型角度，土地利用作用于水土流失過(guò)程主要體現(xiàn)在用地類型改變、用地方式差異、產(chǎn)沙源相對(duì)入河口的距離等方面[6]。研究表明，要達(dá)到切實(shí)有效控制土壤侵蝕的目的，除了改變景觀類型，還要實(shí)現(xiàn)景觀類型在空間上的合理配置[7-8]。異質(zhì)景觀的空間分布格局作用于水文結(jié)構(gòu)和產(chǎn)沙系統(tǒng)，會(huì)造成景觀破碎化程度、水文連通性等指標(biāo)的改變，增強(qiáng)或削弱景觀斑塊產(chǎn)生或攔截水土流失的能力，最終直觀表現(xiàn)在河流斷面泥沙輸出數(shù)量的變化[9]。

當(dāng)前，關(guān)于景觀格局與水土流失過(guò)程關(guān)系的研究大多數(shù)是基于小區(qū)或坡面試驗(yàn)[10]，集中于單一或少數(shù)幾種土地利用類型、植被分布等因素對(duì)侵蝕產(chǎn)沙過(guò)程的獨(dú)立或綜合影響等[11]，并在此基礎(chǔ)上推演更大尺度上的水沙過(guò)程，而從流域、區(qū)域尺度考慮景觀空間配置因素開(kāi)展的研究稍顯薄弱。隨著尺度上升至流域?qū)用?，流域?nèi)部地表徑流和泥沙的產(chǎn)生機(jī)制及運(yùn)移過(guò)程將會(huì)發(fā)生明顯改變，在大尺度上應(yīng)用有限小尺度模擬預(yù)測(cè)得出的模型可能會(huì)造成信息遺漏和丟失[12]；因此，在流域尺度上引入“源”“匯”理論，探究景觀格局的時(shí)空分布特征變化與水沙過(guò)程的定量關(guān)系，可為揭示有利于流域可持續(xù)發(fā)展的合理景觀格局配置提供新的角度和理論依據(jù)。

筆者以位于黃河多沙粗沙國(guó)家級(jí)水土流失重點(diǎn)治理區(qū)的延河流域作為研究對(duì)象，基于遙感與地理信息系統(tǒng)等技術(shù)，識(shí)別流域“源”“匯”景觀類型，采取水文統(tǒng)計(jì)分析、景觀指數(shù)提取、建立多元回歸模型等方法，構(gòu)建景觀格局與平均輸沙率的量化模型，以期為延河流域未來(lái)的生態(tài)恢復(fù)和流域科學(xué)治理工作探索新的方向。

1 研究區(qū)概況

延河是黃河中游的一級(jí)支流，地理坐標(biāo)范圍大致為E 108°37′48″～110°28′48″、N 36°21′～37°19′12″，全長(zhǎng)約286.9 km，流域面積約7 725 km2[13]，95%為黃土丘陵溝壑地貌，土地利用方向主要為草地、林地和耕地。多年平均降水量514 mm，降水多集中在6—9月。土壤類型以黃綿土為主；植被類型主要為森林草原植被，分帶明顯。

為反映該流域的真實(shí)情況，嚴(yán)格控制變量，筆者針對(duì)甘谷驛水文站控制斷面以上區(qū)域(以下簡(jiǎn)稱“研究區(qū)”)開(kāi)展工作；該水文站位于延安市寶塔區(qū)，為延河流域的把口水文站，建站時(shí)間早，長(zhǎng)時(shí)間序列資料較完整，控制流域面積約5 891 km2。

2 數(shù)據(jù)與方法

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

土地利用等地理信息數(shù)據(jù)來(lái)源于中國(guó)黃土高原資源與環(huán)境遙感系列圖編委會(huì)以及中國(guó)科學(xué)院生態(tài)環(huán)境研究中心；甘谷驛站水文泥沙數(shù)據(jù)來(lái)源于水利部水文局(1969—1996年)和黃河水利委員會(huì)中游水文水資源局(1997—2019年)存檔的逐日實(shí)測(cè)資料；降水量數(shù)據(jù)通過(guò)向國(guó)家氣象科學(xué)數(shù)據(jù)中心申請(qǐng)獲得，具體站點(diǎn)為位于研究區(qū)中部區(qū)域的安塞氣象站(編號(hào)53841，地理位置E 109°19′12″、N 36°52′48″)。

2.2 研究方法

2.2.1 泥沙過(guò)程時(shí)變趨勢(shì) 以年為單位對(duì)1969—2019年的逐日輸沙實(shí)測(cè)資料進(jìn)行整理統(tǒng)計(jì)，采用Mann-Kendall趨勢(shì)檢驗(yàn)法判斷其變化趨勢(shì)[14]，采用單位降水距平累積曲線法判斷降水因素對(duì)研究區(qū)泥沙過(guò)程的影響程度[15]。

2.2.2 “源”“匯”景觀類型分類 以GB/T 21010—2017《土地利用現(xiàn)狀分類》為參考，結(jié)合研究區(qū)實(shí)際情況及本研究的需要，將土地利用類型劃分為耕地、林地、草地、建設(shè)用地、水域和其他用地。根據(jù)“源”和“匯”理論，將耕地、建設(shè)用地和其他用地歸類為“源”景觀類型，將林地、草地和水域歸類為“匯”景觀類型(受遙感影像分辨率所限，將園地納入林地一并考慮)。

2.2.3 景觀指數(shù)篩選及提取 景觀指數(shù)是體現(xiàn)景觀幾何特性和配置結(jié)構(gòu)的重要指標(biāo)，本次分別從類型水平和景觀水平選取部分代表性較強(qiáng)的常見(jiàn)指數(shù)(表1)。景觀指數(shù)基于景觀結(jié)構(gòu)數(shù)量化軟件Fragstats 4.2進(jìn)行提取，指數(shù)計(jì)算公式及其生態(tài)學(xué)意義詳見(jiàn)文獻(xiàn)[16]。

表1 選取的景觀指數(shù)及尺度Tab.1 Selected landscape indexes and scales

2.2.4 線性相關(guān)分析 對(duì)于固定面積的流域，“源”“匯”景觀在理論上應(yīng)呈現(xiàn)出相反的變化趨勢(shì)，即在類型景觀水平上，“源”“匯”景觀指數(shù)與流域產(chǎn)沙量若分別表現(xiàn)為反向相關(guān)，可認(rèn)為是有效指標(biāo)。在類型水平上，對(duì)同屬于“源”或“匯”景觀類型的各類景觀要素，若景觀指數(shù)與平均輸沙率均表現(xiàn)為同向相關(guān)，則可認(rèn)為是有效指標(biāo)[17]，納入多元回歸模型的解釋變量集合。

景觀格局指數(shù)與產(chǎn)沙過(guò)程的相關(guān)性借助皮爾遜簡(jiǎn)單相關(guān)系數(shù)R進(jìn)行判斷

(1)

式中：X、Y為變量；Cov(X,Y)為X與Y的協(xié)方差；Var[X]為X的方差；Var[Y]為Y的方差。

2.2.5 多元回歸模型 通過(guò)構(gòu)建多元回歸模型，確定景觀格局指數(shù)與平均輸沙率的定量相關(guān)關(guān)系。模型分析基于SPSS Statistics 21軟件包進(jìn)行，一般如式(2)所示：

y=β0+β1x1+β2x2+…+βpxp+ε。

(2)

式中：y為響應(yīng)變量；x1，x2，…，xp為解釋變量；β0為回歸常數(shù)；β1，β2，…，βp為偏相關(guān)系數(shù)；ε為隨機(jī)誤差。

3 結(jié)果與分析

3.1 泥沙過(guò)程時(shí)變特征

3.1.1 泥沙變化趨勢(shì)性 1969—2019年間，研究區(qū)年及汛期(6—9月)平均輸沙率分別為996.58和2 930.92 kg/s(圖1)，整體波動(dòng)幅度較大，最大值(1977年)與最小值(2015年)相差約276倍，汛期輸沙量占全年比例達(dá)到97.47%，年內(nèi)及年際差異均極其顯著。

圖1 1969—2019年及汛期平均輸沙率變化趨勢(shì)Fig.1 Variation trend of annual and flood season average sediment transport rate from 1969 to 2019

年及汛期平均輸沙率的Mann-Kendall趨勢(shì)檢驗(yàn)結(jié)果表明，統(tǒng)計(jì)量Z值分別為-4.76和-9.03，說(shuō)明研究區(qū)年及汛期平均輸沙率均呈現(xiàn)下降趨勢(shì)，而|Z|>Z1-α/2=2.58，表明其通過(guò)了給定顯著性水平α為0.01的檢驗(yàn)，下降趨勢(shì)極其顯著。

3.1.2 泥沙變化階段性 為判斷降水因素對(duì)泥沙過(guò)程的影響程度，分別繪制輸沙率(圖2)和單位降水輸沙率的距平累積曲線(圖3)。通過(guò)對(duì)比可以看出，2條曲線均存在一定程度的波動(dòng)，但波動(dòng)規(guī)律和趨勢(shì)幾近一致，說(shuō)明降水因素并不是造成研究區(qū)產(chǎn)沙輸沙量波動(dòng)變化的主要因素。

圖2 年平均輸沙率距平累積曲線Fig.2 Cumulative departure curve of annual sediment transport rate

圖3 單位降水輸沙率距平累積曲線Fig.3 Cumulative departure curve of annual sediment transport rate by unit precipitation

3.2 “源”“匯”景觀的時(shí)變趨勢(shì)

3.2.1 “源”“匯”景觀類型的時(shí)變趨勢(shì) 1985—2015年，研究區(qū)“源”“匯”景觀各類型以2000年為界呈現(xiàn)較為明顯的分時(shí)段變化趨勢(shì)(圖4)：1985—2000年，景觀結(jié)構(gòu)基本維持穩(wěn)定，耕地、草地是優(yōu)勢(shì)景觀要素；2000年之后，“源”“匯”景觀類型結(jié)構(gòu)呈現(xiàn)顯著變化，草地、林地面積大幅增加，逐漸成為研究區(qū)內(nèi)的優(yōu)勢(shì)景觀要素，耕地面積急劇減少。

圖4 1985—2015年各類“源”“匯”景觀比例Fig.4 Percentage of the various “source” and “sink” landscape from 1985 to 2015

3.2.2 “源”“匯”景觀整體時(shí)變趨勢(shì) 整體上看，研究時(shí)段內(nèi)研究區(qū)“匯”景觀面積比例始終高于“源”景觀(圖5和圖6)，其中，1985—2000年，“源”景觀與“匯”景觀分布呈現(xiàn)相對(duì)均勢(shì)且較為穩(wěn)定，“匯”景觀的面積比例在55%左右浮動(dòng)，略高于“源”景觀；2000年之后，“源”景觀與“匯”景觀的均勢(shì)分布狀態(tài)開(kāi)始發(fā)生傾斜，“匯”景觀面積大幅增加，至2015年，“匯”景觀面積比例達(dá)到91.0%。

圖5 1985—2015年 “源”“匯”景觀比例Fig.5 Percentage of the “source” and “sink” landscape from 1985 to 2015

圖6 1985—2015年研究區(qū)“源”“匯”景觀結(jié)構(gòu)Fig.6 “Source” and “sink” landscape structure of the study region from 1985 to 2015

3.3 “源”“匯”景觀格局變化對(duì)泥沙過(guò)程的影響

3.3.1 最優(yōu)響應(yīng)時(shí)段的選取 鑒于遙感影像解譯獲取的景觀結(jié)構(gòu)信息僅能代表研究區(qū)一段時(shí)間內(nèi)相對(duì)穩(wěn)定且分布平均的景觀格局特征，而景觀格局結(jié)構(gòu)發(fā)生變化之后需要經(jīng)歷一段時(shí)間才能在水土資源循環(huán)系統(tǒng)中發(fā)揮作用，因此，研究以生態(tài)建設(shè)工程為主要驅(qū)動(dòng)力的“源”“匯”景觀格局變化與泥沙過(guò)程的關(guān)系，需要尋求最優(yōu)響應(yīng)時(shí)段。

如表2所示，不同年份耕地面積和林草地面積與若干不同時(shí)間段內(nèi)的年均輸沙率均有較強(qiáng)的相關(guān)性，其中，對(duì)應(yīng)年份以來(lái)3年的平均輸沙率與耕地面積、林草地面積相關(guān)性最為顯著，由此確定本研究的最優(yōu)響應(yīng)時(shí)段為獲取土地利用信息對(duì)應(yīng)年份以來(lái)3年。

表2 不同年份耕地面積(X)、林草地面積(Y)與對(duì)應(yīng)年份以來(lái)不同時(shí)段平均輸沙率(Zi)決定系數(shù)Tab.2 Determination coefficient between cultivated land area (X), forest and grassland area (Y) in different years and average sediment transport rate (Zi) in different periods since corresponding year

3.3.2 泥沙過(guò)程對(duì)類型水平景觀指數(shù)變化的響應(yīng) 以1985、1995、2000、2008和2015年各“源”“匯”景觀斑塊類型以及合并后的“源”“匯”景觀的類型水平景觀指數(shù)為解釋變量，以最優(yōu)響應(yīng)時(shí)段3年的平均輸沙率為響應(yīng)變量，分別計(jì)算皮爾遜簡(jiǎn)單相關(guān)系數(shù)，判斷景觀指數(shù)指示流域泥沙過(guò)程的有效性。

對(duì)于將各類景觀類型合并后的“源”、“匯”景觀，除了景觀形狀指數(shù)，其余參與統(tǒng)計(jì)的8個(gè)類型水平景觀指數(shù)均呈現(xiàn)出反向相關(guān)的態(tài)勢(shì)(圖7)，表明“源”“匯”景觀整體類型水平景觀指數(shù)的波動(dòng)變化作用于泥沙過(guò)程的機(jī)制較為復(fù)雜，需進(jìn)一步分析各種“源”“匯”景觀類型水平景觀指數(shù)的指示能力。

圖7 “源”“匯”景觀類型水平景觀指數(shù)與平均輸沙率相關(guān)性Fig.7 Correlation between the class level index of the “source” and “sink” landscape and average sediment transport rate

對(duì)于各類“源”景觀斑塊類型，斑塊密度、平均斑塊面積、相似臨近比例、聚合度與年均輸沙率表現(xiàn)出同向相關(guān)性(圖8)，這4個(gè)指標(biāo)可有效指示“源”景觀斑塊變化對(duì)流域產(chǎn)沙量的影響。對(duì)于各類“匯”景觀斑塊類型，斑塊密度、平均臨近指數(shù)與年均輸沙率表現(xiàn)出同向相關(guān)性(圖9)，均為負(fù)相關(guān)，這4個(gè)指標(biāo)可有效指示“匯”景觀斑塊變化對(duì)流域產(chǎn)沙量的影響。

圖8 各類“源”景觀類型水平景觀指數(shù)與平均輸沙率相關(guān)性Fig.8 Correlation between the class level index of the various “source” landscape and average sediment transport rate

圖9 各類“匯”景觀類型水平景觀指數(shù)與平均輸沙率相關(guān)性Fig.9 Correlation between the class level index of the various “sink” landscape and average sediment transport rate

3.3.3 泥沙過(guò)程對(duì)景觀水平景觀指數(shù)變化的響應(yīng) 將最優(yōu)響應(yīng)時(shí)段的平均輸沙率與研究區(qū)景觀水平景觀指數(shù)進(jìn)行相關(guān)分析(圖10)，結(jié)果顯示，個(gè)別指數(shù)表現(xiàn)出較強(qiáng)的相關(guān)性，但相關(guān)系數(shù)整體較類型水平景觀指數(shù)偏低，指示能力稍遜。

圖10 景觀水平景觀指數(shù)與平均輸沙率相關(guān)性Fig.10 Correlation between the class level landscape index and average sediment transport rate

此外，景觀水平景觀指數(shù)的衡量對(duì)象是流域整體景觀，“源”景觀和“匯”景觀的類型組成和結(jié)構(gòu)特征變化均可引起景觀水平景觀指數(shù)的波動(dòng)，因而難以識(shí)別對(duì)流域產(chǎn)沙量變化起決定作用的景觀指數(shù)。

3.4 “源”“匯”景觀格局與泥沙過(guò)程的量化關(guān)系

分別以有效的“源”“匯”景觀類型水平景觀指數(shù)為解釋變量，以最優(yōu)響應(yīng)時(shí)段的平均輸沙率為響應(yīng)變量，建立多元線性回歸模型，定量分析景觀指數(shù)對(duì)泥沙過(guò)程的指示程度。

對(duì)于“源”景觀，提出假設(shè)模型為

y=β0+β1x1+β2x2+β3x3+β4x4+ε。

(3)

式中：y為最優(yōu)響應(yīng)時(shí)段的平均輸沙率，kg/s；x為“源”景觀類型水平景觀指數(shù)，量綱為1，其中，x1為PD，x2為AREA_MN，x3為PLADJ，x4為AI。

在SPSS Statistics 21中輸入數(shù)據(jù)進(jìn)行多元回歸分析，初步回歸結(jié)果顯示：x3變量被排除；x1與x4的方差膨脹因子(variance inflation factor，VIF)遠(yuǎn)大于10，存在顯著的多重共線性，剔除x1，x4參與進(jìn)一步回歸。將保留的x2和x4與y再次進(jìn)行回歸，結(jié)果為

y=2 584.923+7.164x2-32.207x4+ε。

(4)

(2 878.665) (1.696) (35.008)

n=5R2=0.939F=15.304

回歸方程復(fù)相關(guān)系數(shù)為0.969，擬合優(yōu)度達(dá)到0.939，F(xiàn)檢驗(yàn)在0.1的水平上顯著，由此判斷，回歸方程整體較顯著，在“源”景觀類型中，平均斑塊面積均與年均輸沙率正相關(guān)，聚合度與年均輸沙率負(fù)相關(guān)。研究時(shí)段內(nèi)，“源”景觀類型的平均斑塊面積大幅減小，聚合度略有波動(dòng)減小，這與多元回歸的結(jié)果并不完全一致，表明2個(gè)因素并非獨(dú)立作用于泥沙過(guò)程，而是共同影響年均輸沙率的波動(dòng)變化，平均斑塊面積的影響程度相較更高。

對(duì)于“匯”景觀，采取相同的分析步驟，提出假設(shè)模型為：

y=β0+β1x1+β2x2+ε。

(5)

式中：y為最優(yōu)響應(yīng)時(shí)段的平均輸沙率，kg/s；x為“匯”景觀類型水平景觀指數(shù)，量綱為1，其中，x1為PD；x2為PROX_MN?；貧w結(jié)果如下：

y=371.439+5 765.753x1-0.001x2+ε。

(6)

(291.113) (2 130.384) (0.000)

n=5R2=0.945F=17.029

回歸方程復(fù)相關(guān)系數(shù)為0.972，擬合優(yōu)度達(dá)到0.945，F(xiàn)檢驗(yàn)在0.1的水平上顯著，由此判斷，回歸方程整體較顯著，在“匯”景觀類型中，斑塊密度與年均輸沙率正相關(guān)，平均臨近指數(shù)與年均輸沙率負(fù)相關(guān)。研究時(shí)段內(nèi)，“匯”景觀類型的斑塊密度略有波動(dòng)減小，平均臨近指數(shù)成倍大幅增加，表明“匯”景觀斑塊破碎化程度降低，斑塊間鄰近度逐漸升高，景觀連接性趨好，攔蓄泥沙的能力增強(qiáng)，利于削弱泥沙的產(chǎn)生和匯集過(guò)程，減少入河泥沙量。

4 結(jié)論

1)1969—2019年間，研究區(qū)汛期輸沙量占全年比例超過(guò)97%；年及汛期平均輸沙率通過(guò)了給定顯著性水平α為0.01的Mann-Kendall趨勢(shì)檢驗(yàn)，下降趨勢(shì)極其顯著。對(duì)比輸沙率和單位降水輸沙率的距平累積曲線，排除了降水因素主導(dǎo)研究區(qū)產(chǎn)沙輸沙量波動(dòng)變化的可能性。

2)1985—2019年間，研究區(qū)“匯”景觀逐漸占據(jù)主導(dǎo)格局，引起這種波動(dòng)變化的主因是退耕還林還草等生態(tài)建設(shè)工程導(dǎo)致耕地面積大幅減少，林、草地面積大幅增加。在提出最優(yōu)響應(yīng)時(shí)段的基礎(chǔ)上，通過(guò)計(jì)算皮爾遜簡(jiǎn)單相關(guān)系數(shù)，選取可有效指示研究區(qū)泥沙過(guò)程的景觀指數(shù)，建立多元回歸模型。結(jié)果表明：“源”景觀類型的平均斑塊面積和聚合度共同影響年均輸沙率的波動(dòng)，“匯”景觀類型的斑塊密度和平均臨近指數(shù)共同影響年均輸沙率的波動(dòng)，模型擬合優(yōu)度分別達(dá)到0.939和0.945，影響作用較為顯著。

對(duì)于景觀水平景觀指數(shù)，相關(guān)系數(shù)整體較類型水平景觀指數(shù)偏低，且無(wú)法區(qū)分引起指數(shù)波動(dòng)的主導(dǎo)要素是“源”景觀類型還是“匯”景觀類型，因而難以識(shí)別對(duì)流域產(chǎn)沙量變化起決定作用的景觀指數(shù)。

3)多元回歸的結(jié)果表明，影響泥沙過(guò)程的景觀格局要素主要為景觀的破碎化程度和景觀的聚合連接性，當(dāng)“源”景觀的破碎化程度降低，“匯”景觀的鄰近度升高、連接性增強(qiáng)時(shí)，形成的景觀格局結(jié)構(gòu)利于流域尺度水文循環(huán)系統(tǒng)的泥沙攔蓄和調(diào)控。