余 欣，李其江，劉希勝，閔 敏，楊 陽(yáng)，時(shí) 璐，馬 蕊

(1.黃河水利委員會(huì) 黃河水利科學(xué)研究院，河南 鄭州 450003;2.青海省水利廳，青海 西寧 810001;3.青海省水文水資源測(cè)報(bào)中心，青海 西寧 810001）

湟水流域位于我國(guó)西北地區(qū)的生態(tài)過(guò)渡帶,干流西寧站以下的黃土丘陵區(qū)水土流失嚴(yán)重[1］。近年來(lái),受人類活動(dòng)和氣候變化的影響,湟水流域水沙關(guān)系發(fā)生了很大變化。

本文以湟水流域民和站以上區(qū)域?yàn)檠芯繉?duì)象,基于水沙長(zhǎng)時(shí)間序列資料,利用趨勢(shì)分析、突變檢驗(yàn)等方法,從侵蝕性降水量(日降水量大于12 mm)、NDVI指數(shù)和水土保持措施等三個(gè)方面,分析了水沙情勢(shì)演變特征及其驅(qū)動(dòng)因素,以期為水資源開發(fā)利用、水旱災(zāi)害防治、水土流失治理等提供參考。

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

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

收集了湟水干流湟源、西寧和樂都3個(gè)青海省屬水文站實(shí)測(cè)徑流、泥沙數(shù)據(jù)以及黃河水利委員會(huì)水文局所屬湟水干流民和水文站實(shí)測(cè)徑流、泥沙數(shù)據(jù)。數(shù)據(jù)序列為1956—2020年,資料質(zhì)量較好、序列完整。

1.2 研究方法

1.2.1 趨勢(shì)檢驗(yàn)

趨勢(shì)檢驗(yàn)采用Sen斜率估計(jì)法和Mann-Kendall[1］檢驗(yàn)法(M-K法)。Sen斜率估計(jì)法是一種穩(wěn)健的非參數(shù)統(tǒng)計(jì)趨勢(shì)計(jì)算方法,該方法計(jì)算效率高,對(duì)測(cè)量誤差和離群數(shù)據(jù)不敏感,常被用于長(zhǎng)時(shí)間序列數(shù)據(jù)的趨勢(shì)分析中。Sen斜率β公式為

式中:median為中值函數(shù)；Xi、Xj為序列中數(shù)據(jù)值。

β值為正表示上升趨勢(shì),β值為負(fù)表示下降趨勢(shì)。

M-K[3-5］法是一種非參數(shù)統(tǒng)計(jì)檢驗(yàn)方法,變量可以不具有正態(tài)分布特征,適用于水文變量的趨勢(shì)檢驗(yàn)。假設(shè)時(shí)間序列(X1,X2,…,Xn),統(tǒng)計(jì)量S為

方差Var(S):

檢驗(yàn)指標(biāo)Z:

在雙邊趨勢(shì)檢驗(yàn)中,對(duì)于給定的置信水平α,若則在置信水平α上時(shí)間序列數(shù)據(jù)存在明顯的上升或下降趨勢(shì)。Z為正值表示上升趨勢(shì),Z為負(fù)值表示下降趨勢(shì)。表示其通過(guò)了置信度為95%的顯著性檢驗(yàn)。

以上兩種方法結(jié)合可以增強(qiáng)研究方法的抗噪性,并在一定程度上提高檢驗(yàn)成果的準(zhǔn)確性。

1.2.2 突變檢驗(yàn)

突變檢驗(yàn)采用Pettitt檢驗(yàn)法[2,6］。Pettitt檢驗(yàn)法是非參數(shù)檢驗(yàn)法,該方法在檢驗(yàn)連續(xù)序列時(shí)是高效的,能夠較好地識(shí)別水文時(shí)間序列的突變點(diǎn),在突變點(diǎn)檢驗(yàn)中應(yīng)用較多,且物理意義清晰。該檢驗(yàn)基于Mann-Whitney的統(tǒng)計(jì)函數(shù)認(rèn)為兩個(gè)序列(X1,X2,…,Xt)和(Xt＋1,Xt＋2,…,XT)均來(lái)自于同一序列,對(duì)于連續(xù)的序列,Ut,T和Vt,T由下式計(jì)算:

式中:Ut,T、Vt,T為不同時(shí)段檢驗(yàn)值統(tǒng)計(jì)量。

當(dāng)｜Ut,T｜最大時(shí),對(duì)應(yīng)的Xt是可能的突變點(diǎn)。當(dāng)突變點(diǎn)Ut,T＞0時(shí),該序列具有向下突變趨勢(shì),反之則具有向上突變的趨勢(shì)?？赡芡蛔凕c(diǎn)的顯著性水平由下式計(jì)算:

當(dāng)POA(t)≤0.5時(shí),該點(diǎn)為有效突變點(diǎn)。

1.2.3 周期檢驗(yàn)

Morlet小波[7-9］作為一種分解工具和檢測(cè)長(zhǎng)時(shí)間序列中多尺度突變的方法,具有時(shí)域和頻域的良好局部化特征,而且隨著信號(hào)不同頻率成分在時(shí)域取樣的疏密,調(diào)節(jié)其放大倍數(shù),可以清楚地分析信號(hào)在各個(gè)層次上的變化,進(jìn)而對(duì)源信號(hào)在不同時(shí)間尺度上的變化趨勢(shì)進(jìn)行判斷和甄別。

2 結(jié)果與分析

2.1 水沙變化分析

2.1.1 趨勢(shì)性

采用湟水干流湟源、西寧、樂都、民和4個(gè)水文站的徑流、泥沙數(shù)據(jù),分析1956—2020年湟水流域從上游到中上游、中游、中下游主要水文站的天然年水沙時(shí)序,見圖1、圖2。

圖1 湟水干流各水文站天然年徑流量過(guò)程線

圖2 湟水干流各水文站天然年輸沙量過(guò)程線

通過(guò)Sen斜率估計(jì)法和M-K趨勢(shì)檢驗(yàn)法,分析各站徑流量、輸沙量變化趨勢(shì),結(jié)果見表1。1956—2020年各水文站徑流量呈增加趨勢(shì),增加速率在0.2億～0.5億m3/10 a之間,M-K檢驗(yàn)指標(biāo)Z在2.04～3.40之間,呈顯著上升趨勢(shì)。同期各站輸沙量呈減小趨勢(shì),減少速率在3.3萬(wàn)～386.2萬(wàn)t/10 a之間,M-K檢驗(yàn)指標(biāo)Z在-6.8～-2.2之間,呈顯著下降趨勢(shì)。

表1 湟水干流水沙趨勢(shì)分析結(jié)果

通過(guò)Pettitt突變檢驗(yàn)可得,各站徑流量突變發(fā)生在1982年和2002年左右,見表2。1982—2001年與1956—1981年相比多年平均徑流量增加3.9%～11.6%,2002—2020年與1982—2001年相比多年平均徑流量增加7.5%～19.7%。各站輸沙量在2000年左右發(fā)生突變,2001—2020年與1956—2000年相比多年平均輸沙量減少45.6%～79.3%。

表2 湟水干流水沙突變分析結(jié)果

2.1.2 周期性

用小波分析法識(shí)別湟水流域1956—2020年徑流量、輸沙量序列周期成分。湟水徑流量上下游呈現(xiàn)相對(duì)一致的周期變化,在1956—2020年主要存在中心時(shí)間尺度為9、22、34、47 a的周期分布,在34 a的周期尺度上對(duì)應(yīng)小波方差最大,說(shuō)明周期震蕩最強(qiáng),為主要周期,見圖3。

圖3 湟水干流各水文站徑流量小波周期分析

從湟水干流各站輸沙量小波周期分析可以看出,上游湟源站與下游3站周期變化不同(見圖4)。上游的湟源站以30 a尺度為主周期；中下游西寧、樂都、民和站的小波變換周期類似,存在11～25 a尺度的震蕩,尤其在11、14、21 a尺度的小波震蕩能量較強(qiáng)。

圖4 湟水干流各水文站輸沙量小波周期分析

2.2 水沙關(guān)系分析

2.2.1 水沙異源

按照水文站區(qū)間控制單元,將流域劃分為湟源以上、湟源—西寧、西寧—樂都、樂都—民和4個(gè)區(qū)間。經(jīng)過(guò)分析,流域產(chǎn)水、產(chǎn)沙存在區(qū)間差異性,水沙存在明顯異源現(xiàn)象。

流域水量主要來(lái)源于湟源至西寧區(qū)間,該區(qū)間貢獻(xiàn)水量和沙量分別占民和站實(shí)測(cè)的41.4%和19.1%；而沙量主要來(lái)源于樂都至民和區(qū)間,該區(qū)間貢獻(xiàn)水量和沙量分別占民和站實(shí)測(cè)的16.1%和42.7%,兩個(gè)區(qū)間水量和沙量貢獻(xiàn)率大小相反,見表3。

表3 湟水流域各區(qū)間水沙量對(duì)比

考慮各水文站控制面積不同,分別計(jì)算各區(qū)間的產(chǎn)水模數(shù)與產(chǎn)沙模數(shù),分析可知,4個(gè)區(qū)間的產(chǎn)水模數(shù)接近,在9.48萬(wàn)～11.31萬(wàn)m3/km2之間。因?yàn)楫a(chǎn)水模數(shù)差異不大,所以區(qū)間面積的大小是產(chǎn)水量多少的主導(dǎo)因素,由于湟源—西寧區(qū)間面積最大,因此其成為主要產(chǎn)水區(qū)間。

從上游至下游,4個(gè)區(qū)間的產(chǎn)沙模數(shù)成倍增加,由湟源以上的134.16 t/km2增加到樂都—民和區(qū)間的2289.86 t/km2,在產(chǎn)水模數(shù)接近的前提下,4個(gè)區(qū)間的水土流失程度呈現(xiàn)自上游至下游成倍遞增趨勢(shì)。

2.2.2 上下游水沙關(guān)聯(lián)性

湟水干流上下游相鄰兩站年徑流量相關(guān)性較好,由上游至下游兩兩相關(guān)關(guān)系增強(qiáng),說(shuō)明年徑流自上游至下游具有良好的協(xié)同關(guān)系,見圖5。

圖5 湟水干流各站實(shí)測(cè)年徑流量上下游相鄰兩站相關(guān)關(guān)系(單位:億m3)

在相鄰兩站實(shí)測(cè)年輸沙量相關(guān)關(guān)系中,除湟源與西寧站的相關(guān)關(guān)系較差外,中下游各站相關(guān)性較好,說(shuō)明湟水流域中下游主要產(chǎn)沙地區(qū)各站的協(xié)同關(guān)系好于上游,見圖6。

圖6 湟水干流各站實(shí)測(cè)年輸沙量上下游相鄰兩站相關(guān)關(guān)系(單位:萬(wàn)t)

2.2.3 水沙時(shí)序變化協(xié)同作用

湟水干流各站徑流量與輸沙量的關(guān)系散亂,相關(guān)性較差,呈現(xiàn)不同變化趨勢(shì)。時(shí)序上,各站年徑流量整體呈現(xiàn)上升態(tài)勢(shì),輸沙量則呈現(xiàn)下降態(tài)勢(shì)。

綜合徑流輸沙突變時(shí)間、歷年過(guò)程以及下墊面變化等情況,將長(zhǎng)程時(shí)序分為1956—1980年、1981—2002年以及2003—2020年3個(gè)時(shí)段進(jìn)行縱深分析。各站實(shí)測(cè)年徑流量總體呈現(xiàn)上升態(tài)勢(shì),在2002年以后,各站均值、中位數(shù)、四分位值(上、下四分位值)、最大最小值均有所上升。3個(gè)階段,各站實(shí)測(cè)年輸沙量呈現(xiàn)下降態(tài)勢(shì),尤其是西寧、樂都、民和站輸沙量持續(xù)下降,其均值、中位數(shù)、四分位值(上、下四分位值)、最大最小值均有所下降,見圖7。

圖7 湟水干流各站實(shí)測(cè)年徑流量、輸沙量不同時(shí)段對(duì)比

2.3 驅(qū)動(dòng)力分析

2.3.1 降水指標(biāo)

降水是湟水流域水資源的主要補(bǔ)給來(lái)源,但侵蝕性降水同時(shí)會(huì)導(dǎo)致土壤水蝕。本文選用日降水量大于12 mm的年侵蝕性降水量(PRCPTOT12)、侵蝕性降水天數(shù)(R12)、侵蝕性降水強(qiáng)度(SDII12)三個(gè)降水指標(biāo),分析湟水1981—2020年P(guān)RCPTOT12、R12、SDII12的長(zhǎng)期變化趨勢(shì)(見圖8)。湟水流域PRCPTOT12、R12趨勢(shì)表現(xiàn)基本一致,年代際呈現(xiàn)“先增后減”態(tài)勢(shì)。2001—2010年年平均PRCPTOT12為研究時(shí)段最大,達(dá)到211.0 mm,PRCPTOT12以9.00 mm/10 a速率呈不顯著上升趨勢(shì),最高值為2009年的264.2 mm,最低值為1982年的109.8 mm。湟水流域研究時(shí)段內(nèi)多年平均R12以0.43 d/10 a速率呈不顯著上升趨勢(shì),總體波動(dòng)較大,變化范圍在6.0～14.2 d之間,1991—2000年時(shí)段均值上升,至2000年后轉(zhuǎn)為偏多趨勢(shì),在2001—2010年為研究時(shí)段最大,R12達(dá)到11.5 d,2001—2010年、2011—2020年兩時(shí)期分別較均值偏多0.8、0.3 d。湟水流域研究時(shí)段內(nèi)多年平均SDII12為18.3 mm/d,變化范圍在16.2～20.1 mm/d之間,總體波動(dòng)不大,年代際呈現(xiàn)“增—減—增”特征。黃土高原地區(qū)降水量顯著減少,減少幅度自東南向西北逐漸變小,侵蝕性降水量與年降水量變化相似,但湟水流域降水量和侵蝕性降水量呈增加趨勢(shì),不同于黃土高原變化趨勢(shì)[10］,這與本文研究結(jié)果類似。

圖8 侵蝕性降水指標(biāo)變化趨勢(shì)

采用相關(guān)分析法進(jìn)一步分析湟水干流民和站徑流量、輸沙量與降水指標(biāo)的響應(yīng)關(guān)系(見表4)。從相關(guān)性分析結(jié)果看,湟水流域PRCPTOT12、R12與民和站徑流量呈極顯著正相關(guān)關(guān)系,SDII12與輸沙量呈顯著正相關(guān)關(guān)系,其余相關(guān)關(guān)系不顯著。這表明PRCPTOT12、R12的增加是湟水流域徑流量增加的主要因素,SDII12是引起輸沙量變化的重要因素。這也意味著盡管增加趨勢(shì)不顯著,但隨著極端天氣事件的增多,由水蝕導(dǎo)致的土壤侵蝕總體上還將繼續(xù)增多,這必然對(duì)湟水流域高質(zhì)量發(fā)展產(chǎn)生深遠(yuǎn)影響。

表4 民和站徑流量、輸沙量與降水指標(biāo)的相關(guān)系數(shù)

2.3.2NDVI指數(shù)

植被變化受人類活動(dòng)、氣候條件等因素共同作用,不同下墊面條件導(dǎo)致水熱條件空間分布不均,進(jìn)而造成NDVI指數(shù)空間異質(zhì)性[11］。植被覆蓋也是衡量水沙變化的關(guān)鍵因素之一[12］。湟水流域NDVI年均值與5—9月均值空間分布見圖9,NDVI年均值介于0.13～0.68之間,區(qū)域性差異明顯,表現(xiàn)為由西北向東南方向遞減,從上游向中下游遞減,表明西北部植被覆蓋率高于東南部,這與王兮之等[13-14］研究結(jié)論相似。流域西北部多為海拔較高、氣候寒冷區(qū)域,植被以林地和草地為主,區(qū)域內(nèi)氣溫低,熱量少,蒸發(fā)量小,城鄉(xiāng)居民建設(shè)用地較少,受人類活動(dòng)干擾小,植被覆蓋度相對(duì)較高,降水的侵蝕產(chǎn)沙能力較弱。被淺山區(qū)包圍的中下游河谷地帶,地勢(shì)相對(duì)較低,水熱條件較好,人類活動(dòng)頻繁,沿河兩岸遍布城鄉(xiāng)居民建設(shè)用地、工礦企業(yè)和耕地,NDVI較上游偏低,下墊面截留固沙能力較弱。

圖9 湟水流域NDVI指數(shù)空間分布

湟水流域降水主要集中在5—9月,在此期間NDVI均值較年均值有所提高,尤其中上游區(qū)域變化最為明顯。夏季降水較多,上游海拔較高,氣溫相對(duì)下游偏低,降水蒸發(fā)較慢,雨水經(jīng)下滲進(jìn)入土壤促進(jìn)植物生長(zhǎng),因此在一定程度上可以減少雨水侵蝕、增強(qiáng)地表滲透性與抗蝕性[15］,但局部強(qiáng)降水易沖刷植被覆蓋度較低的湟水谷地坡面和溝道,造成輸沙量增大,進(jìn)一步印證中下游是輸沙的主要來(lái)源。

2.3.3 水土保持措施

湟水流域人口占青海省總?cè)丝诘?0%左右,人類活動(dòng)頻繁,下墊面條件復(fù)雜。自2000年以來(lái),青海省水利部門通過(guò)新建基本農(nóng)田、營(yíng)造水保林和經(jīng)濟(jì)林、封育治理等措施有效改善湟水流域水土流失狀況,重點(diǎn)實(shí)施了小流域綜合治理、淤地壩壩系建設(shè)、坡耕地水土流失綜合防治、鞏固退耕還林區(qū)基本口糧田工程等一系列生態(tài)修復(fù)和預(yù)防保護(hù)工程,對(duì)于減少水土流失、增加土壤入滲、攔蓄坡溝泥沙和地表徑流在一定程度上發(fā)揮了重要作用。根據(jù)全國(guó)第一次水利普查數(shù)據(jù),湟水水利工程數(shù)量增長(zhǎng)不明顯,但隨著黑泉水庫(kù)建成累計(jì)總庫(kù)容陡增,這也是西寧—民和區(qū)間輸沙量減少的原因之一。累計(jì)建設(shè)淤地壩661座,小型蓄水工程62864個(gè)。湟水中下游黃土丘陵區(qū)水土流失得到抑制,2000年之后輸沙量顯著下降,治理效果、生態(tài)效益顯著提高。

3 結(jié)論

(1)湟水干流4站年徑流量和輸沙量分別呈現(xiàn)不同變化趨勢(shì)。各站徑流量均呈增加趨勢(shì),速率為0.2億～0.5億m3/10 a,突變發(fā)生在1982年和2002年左右。各站輸沙量則呈下降趨勢(shì),下游減小幅度較大,在2000年左右發(fā)生突變。

(2)湟水流域徑流量上下游呈現(xiàn)相對(duì)一致的周期變化,主要存在中心時(shí)間尺度為9、22、34、47 a的周期分布；輸沙量湟源站存在30 a尺度的主周期,而下游西寧、樂都、民和站較為相似,在11、14、21 a尺度的小波方差較為突出。

(3)湟水流域徑流縱向空間具有連貫協(xié)同性,中下游泥沙協(xié)同關(guān)系好于上游。水沙趨勢(shì)各異,產(chǎn)水、產(chǎn)沙存在明顯異源現(xiàn)象,湟源—西寧區(qū)間面積最大致使其成為主要來(lái)水區(qū)間,樂都—民和區(qū)間為主要產(chǎn)沙區(qū)間。

(4)PRCPTOT12、R12對(duì)徑流量變化具有促進(jìn)作用,SDII12變化引起輸沙量的變化。得益于湟水流域生態(tài)修復(fù)和預(yù)防保護(hù)工程的實(shí)施,流域內(nèi)植被覆蓋度提高,水土流失程度降低、土壤入滲增加、攔蓄坡溝泥沙和地表徑流成效顯著。