張 萍,蔡強(qiáng)國,鄭明國,何天樂

(1.中國科學(xué)院 地理科學(xué)與資源研究所 陸地水循環(huán)及地表過程重點(diǎn)實驗室，北京 100101；2.中國科學(xué)院大學(xué) 資源環(huán)境學(xué)院，北京 10049；3.廣東省科學(xué)院生態(tài)環(huán)境與土壤研究所 華南土壤 污染控制與修復(fù)國家地方聯(lián)合工程研究中心 廣東省農(nóng)業(yè)環(huán)境綜合治理重點(diǎn)實驗室 廣東省面源污染防治工程技術(shù)研究中心,廣東 廣州 510650；4.中國科學(xué)院 地理科學(xué)與資源研究所 資源與環(huán)境信息系統(tǒng)國家重點(diǎn)實驗室，北京 100101)

流域徑流產(chǎn)沙的空間變異性問題的核心是尺度[1]，不同時空尺度上的降雨及徑流產(chǎn)沙過程變化較大，這種現(xiàn)象在黃河流域尤為突出。然而，目前多年時間尺度上對于影響流域產(chǎn)沙的降雨空間變異性研究較少[2]。黃土高原地區(qū)次降雨的空間變異性較大，使得次降雨下侵蝕產(chǎn)沙過程的空間變異大[3]。但有研究[4]表明，多年時間尺度上無定河流域產(chǎn)沙量的空間變異性確很小。降雨是影響流域產(chǎn)沙的主要動力因素之一[5]，特別是汛期降雨影響著黃河流域的侵蝕產(chǎn)沙[6-7]，降雨的空間變化對流域產(chǎn)沙空間分布有著重要影響[8-9]。降雨量空間變異性越大，不同時段流域產(chǎn)沙就越集中，相同降雨量時產(chǎn)沙量就越多[10]。降雨特征空間分布的研究需要在多重時空尺度上廣泛開展[11]。然而，相對于降雨和人類活動導(dǎo)致侵蝕產(chǎn)沙時間變化的大量研究而言，影響侵蝕產(chǎn)沙空間變化的降雨特征的空間分布相關(guān)研究目前不僅數(shù)量上很少[12-14]，這在一定程度上影響了侵蝕產(chǎn)沙空間分布式模型預(yù)報精度的提高[15]。河口—龍門區(qū)間(以下簡稱“河龍區(qū)間”)是黃河泥沙的主要來源區(qū)[16]，但自20世紀(jì)70年代以來該流域水沙特征發(fā)生了明顯的變化，產(chǎn)沙量呈現(xiàn)劇烈減少的趨勢[17]，1972—1998年期間有21 a黃河出現(xiàn)斷流[18]。造成黃河徑流和產(chǎn)沙減少的主要因素是降雨等氣象要素發(fā)生顯著變化以及黃河流域內(nèi)劇烈的人類活動[19]。因此，本研究以黃土高原典型多沙粗沙區(qū)河龍區(qū)間為研究對象，分析了該區(qū)域降雨特征的空間分布形式，以及在大規(guī)模水土流失治理措施實施前后降雨特征對流域產(chǎn)沙影響的變化。

1 材料和方法

1.1 研究區(qū)域概況

河龍區(qū)間位于黃河中游上段，包括陜西、山西省和內(nèi)蒙古自治區(qū)的50個縣，面積約為1.13×105km2。該區(qū)地貌上表現(xiàn)為黃土丘陵溝壑區(qū)、黃土峁?fàn)钋鹆隃羡謪^(qū)、風(fēng)沙黃土丘陵溝壑區(qū)、森林黃土丘陵溝壑區(qū)、濕潤類黃土丘陵溝壑區(qū)、黃土殘垣溝壑區(qū)、風(fēng)沙草原區(qū)和高原土石山區(qū)。該區(qū)屬于溫帶大陸性氣候，年均溫度6～14 ℃；植被以草原為主，自西北向東南依次分布著荒漠草原、干旱草原和森林草原。區(qū)間有83.45%的土地為水土流失區(qū)，多年平均輸沙量高達(dá)9.08×105t，占黃河全年粗沙輸沙量的72%[20]。

1.2 數(shù)據(jù)來源

為研究降雨特征的空間分布及降雨特征對流域產(chǎn)沙的影響，本文引入汛期降雨，汛期降雨侵蝕力和侵蝕產(chǎn)沙模數(shù)指標(biāo)。降雨和侵蝕產(chǎn)沙數(shù)據(jù)來源于《中華人民共和國水文年鑒》、黃土高原科學(xué)數(shù)據(jù)中心(http:∥loess.geodata.cn/)和中國氣象數(shù)據(jù)網(wǎng)(http:∥data.cma.cn/)。利用陜西、山西和內(nèi)蒙古地區(qū)3個省的降雨資料，得到河龍區(qū)間1959—2015年(其中1990—2007年的數(shù)據(jù)資料缺失)有效的雨量站點(diǎn)為79個，除去降雨缺測的年份共計有40 a的有效降雨資料。為了探究河龍區(qū)間不同時間尺度降雨的空間變異特征，本文從半月、月以及1 a，2 a，4 a，5 a，8 a，10 a，20 a，40 a等時間尺度上對降雨特征的空間變異性進(jìn)行分析。

1.3 研究方法

1.3.1 降雨特征的計算 本文所涉及的降雨特征包括汛期降雨量(Pflood,mm)和汛期降雨侵蝕力〔Rflood，MJ·mm/(hm2·h)〕。通過計算河龍區(qū)間1959—2015年多年平均月降雨量占多年平均降雨量的比例，得出該流域年內(nèi)70%以上的降雨量都集中于6—9月。研究表明河龍區(qū)間年內(nèi)的產(chǎn)沙性降雨[21]和96%的輸沙量主要集中在汛期[22]。因此，本文將針對6—9月的降雨和產(chǎn)沙資料進(jìn)行降雨特征和產(chǎn)沙模數(shù)的統(tǒng)計計算。降雨侵蝕力是指降雨引起土壤侵蝕的潛在能力[23]。有研究者對比了國內(nèi)外常用的降雨侵蝕力簡易計算模型在中國東部地區(qū)的適用性[24]，發(fā)現(xiàn)全國第一次水利普查采用的簡易降雨侵蝕力計算模型[25]可直接用于降雨侵蝕力的計算〔公式(1)—(2)〕。

(1)

(2)

式中：Rk為某半月時段的降雨侵蝕力值〔MJ·mm/(hm2·h)〕；n為半月時段內(nèi)的天數(shù)；(Pd≥12)m為半月時段內(nèi)第md≥12 mm的侵蝕性日雨量；α為模型參數(shù)。

1.3.2 面雨量及空間變異系數(shù) 面雨量是指某一時段內(nèi)一定面積上的平均雨量，準(zhǔn)確地計算面雨量對流域產(chǎn)沙的研究具有重要意義[26]。面雨量的計算方法有很多，主要有逐步訂正格點(diǎn)法、三角形法、等雨量法和泰森多邊形法等。泰森多邊形插值是根據(jù)雨量站點(diǎn)降雨資料估算流域面雨量的方法[8](公式3)。但由于其對流域降雨空間變異性的描述能力較差，本文選取泰森多邊形加權(quán)變異系數(shù)(Cv-Thiessen)來計算降雨特征的空間變異性[8](公式5)。降雨特征值及流域產(chǎn)沙的空間變異性用變異系數(shù)指標(biāo)來表示。泰森多邊形加權(quán)的變異系數(shù)Cv-Thiessen越大，離散程度越大，Cv-Thiessen越接近于0，離散程度越小。

產(chǎn)沙模數(shù)表示一個流域范圍內(nèi)單位面積的產(chǎn)沙量，本身就是一個面指標(biāo)，所以產(chǎn)沙模數(shù)的變異系數(shù)不需要再做面積加權(quán)(公式6)。

(3)

(4)

(5)

(6)

1.3.3K均值聚類法(K-means) 聚類分析是將數(shù)據(jù)按照一定規(guī)則進(jìn)行分類的數(shù)據(jù)處理分析過程，同一類內(nèi)的對象相似性較高，而類與類之間的差異較大[27]。K-means聚類法原本是基于距離特征的算法[28]，有研究表明聚類分析法可用于氣候分類區(qū)劃[29]。本文中采用汛期降雨量和汛期降雨侵蝕力值，對降雨特征的空間分布進(jìn)行聚類分析，可將研究區(qū)根據(jù)降雨量和降雨侵蝕力的大小分為多個類型區(qū)。K-means 算法[30]為：

在數(shù)據(jù)集中隨機(jī)選取k個雨量站作為k個簇(類群)的中心，分別計算其它雨量站的汛期降雨量和汛期降雨侵蝕力值到各個簇中心的差值，將這些雨量站點(diǎn)分別劃分到最近的簇中，差值的大小以歐氏距離的計算方法為基礎(chǔ)：

(7)

式中：d(i,j)為兩個雨量站間降雨特征的差值；i=(xi1,xi2,…,xin)，j=(xj1,xj2,…,xjn)為兩個n維的數(shù)據(jù)對象。

根據(jù)聚類結(jié)果，通過計算簇中數(shù)據(jù)對象的算術(shù)平均值重新確定簇中心，然后按照新中心進(jìn)行重新聚類，直到下列函數(shù)收斂(簇中心不再變化)：

(8)

式中：SSE為綜合平方誤差，SSE值越小，表明越接近簇中心，聚類效果越好；p為空間點(diǎn)；mi為簇Ci的平均值。

2 結(jié)果與分析

2.1 不同時間和空間尺度上降雨的空間變異性

本文在不同時間尺度上對汛期降雨和汛期降雨侵蝕力的空間變異性進(jìn)行了分析，分析結(jié)果可以得出，在半月、月以及1 a等年內(nèi)時間尺度上，汛期降雨和汛期降雨侵蝕力的空間變異系數(shù)隨著時間尺度的增大而銳減。汛期降雨和汛期降雨侵蝕力在半月時間尺度上空間變異系數(shù)的平均數(shù)分別為0.66和1.68，相差2.6倍；汛期降雨和汛期降雨侵蝕力在1 a時間尺度上空間變異系數(shù)的平均數(shù)分別為0.25和0.54，相差2.2倍。在年內(nèi)和月際時間尺度上汛期降雨的空間變異性要明顯小于汛期降雨侵蝕力。在2，4，5，8，10，20 a以及40 a等年際時間尺度上，汛期降雨和汛期降雨侵蝕力的空間變異系數(shù)隨著時間尺度的增大而緩慢減小。汛期降雨和汛期降雨侵蝕力在2 a時間尺度上空間變異系數(shù)的平均數(shù)分別為0.44和0.20，相差2.2倍；汛期降雨和汛期降雨侵蝕力在40 a時間尺度上空間變異系數(shù)的平均數(shù)分別為0.20和0.13，相差1.54倍。在年際時間尺度上汛期降雨的空間變異性要明顯大于汛期降雨侵蝕力。在河龍區(qū)間選取了偏關(guān)河、皇甫川、孤山川、朱家川、嵐漪河、窟野河、蔚汾河、禿尾河、佳蘆河、湫水河、三川河、無定河、清澗河、延水、昕水河、汾川河、州川河以及仕望川等18個主要支流，探究了河龍區(qū)間降雨特征在支流上的空間變異特征。

河龍區(qū)間各個支流汛期降雨的空間變異系數(shù)的最大值為0.22，最小為0.11，相差2倍；汛期降雨侵蝕力的空間變異系數(shù)的最大值為0.51，最小為0.25，相差2倍。支流間汛期降雨和汛期降雨侵蝕力的空間變異性相差不大，并且汛期降雨的空間變異性大的支流相應(yīng)的汛期降雨侵蝕力的空間變異性也大。

2.2 降雨對流域輸沙模數(shù)的影響

為分析降雨和流域產(chǎn)沙在1959—2015年大規(guī)模水土保持措施實施前后的相關(guān)性，本文利用線性回歸法，對河龍區(qū)間9個主要支流偏關(guān)河、皇甫川、朱家川、窟野河、蔚汾河、禿尾河、三川河、無定河以及昕水河1959—2015年汛期降雨、汛期降雨侵蝕力等和產(chǎn)沙模數(shù)進(jìn)行回歸分析。結(jié)果發(fā)現(xiàn)(見圖1—2)，在1959—1970年期間，產(chǎn)沙模數(shù)與汛期降雨和汛期降雨侵蝕力呈顯著冪函數(shù)關(guān)系(p<0.01)，決定系數(shù)都在0.7左右，并且在9個支流都表現(xiàn)出產(chǎn)沙模數(shù)隨降雨特征的增加而急劇增加的趨勢。相比之下，1971—2015年期間，產(chǎn)沙模數(shù)與汛期降雨和汛期降雨侵蝕力之間無明顯相關(guān)關(guān)系，這是由于該區(qū)域在20世紀(jì)70年代起實施了大規(guī)模水土保持措施，例如將坡地改造為梯田，使局部微地形改變，田面坡度大大減小，有效增加了降雨的入滲，使產(chǎn)流減少，因而侵蝕產(chǎn)沙減少；黃河中游地區(qū)的溝谷侵蝕和重力侵蝕十分劇烈，大量修筑淤地壩對侵蝕的泥沙進(jìn)行攔截，增加流域中泥沙的存貯，從而減少了產(chǎn)沙量[17]。由于以上人類活動的干預(yù)導(dǎo)致降雨對產(chǎn)沙模數(shù)的貢獻(xiàn)率極大減小。

圖1 河龍區(qū)間主要支流1959-2015年P(guān)flood對流域產(chǎn)沙的影響

圖2 河龍區(qū)間主要支流1959-2015年Rflood對流域產(chǎn)沙的影響

2.3 降雨的空間分布特征

3 結(jié) 論

(1) 1959—2015年河龍區(qū)間汛期降雨和汛期降雨侵蝕力的空間分布特征一致，都表現(xiàn)出從東南先西北區(qū)域減小的趨勢。汛期降雨的空間變異性為12%，汛期降雨侵蝕力的空間變異性為20%，汛期降雨侵蝕力的空間變異性要遠(yuǎn)遠(yuǎn)大于汛期降雨。

汛期降雨和汛期降雨侵蝕力在所有時間尺度上的空間分布特征一致，具有顯著的緯度地帶性(p<0.01)，且汛期降雨侵蝕力的緯度地帶性要強(qiáng)于汛期降雨。汛期降雨和汛期降雨侵蝕力在年內(nèi)時間尺度上空間變異系數(shù)隨時間尺度的增大而銳減，并且汛期降雨的空間變異性小于汛期降雨侵蝕力；在年際時間尺度上汛期降雨和汛期降雨侵蝕力的空間變異性隨時間尺度的增大緩慢減小，汛期降雨的空間變異性大于汛期降雨侵蝕力。

(2) 河龍區(qū)間18個主要支流汛期降雨和汛期降雨侵蝕力的空間變異程度差異不大，汛期降雨空間變異性大的支流相應(yīng)的汛期降雨侵蝕力的空間變異性也大，支流間汛期降雨和汛期降雨侵蝕力空間變異性的差異與其面積大小無關(guān)。

(3) 1959—1970年期間，河龍區(qū)間SSY與Pflood和Rflood呈顯著正相關(guān)關(guān)系(p<0.01)；1971—2015年，SSY與Pflood和Rflood無明顯相關(guān)關(guān)系，降雨對流域產(chǎn)沙的貢獻(xiàn)率大大減小。