珠江流域下游1960—2019水沙異變?cè)\斷及其控制因素

王永紅, 彭 錦, 黃 暢, 吳宏旭, 蔡斯龍, 劉 鋒

(1.中國海洋大學(xué) 海洋地球科學(xué)學(xué)院, 海底科學(xué)與探測技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室, 山東 青島 266100;2.廣東省水文局, 廣州 510150; 3.中山大學(xué) 海洋工程與技術(shù)學(xué)院, 廣東 珠海 519082)

隨著社會(huì)經(jīng)濟(jì)的發(fā)展，人類活動(dòng)對(duì)于河口三角洲水沙變化的影響日益增強(qiáng)。強(qiáng)人類活動(dòng)導(dǎo)致了水沙情勢(shì)的異變，引起水道的洪潮水位、徑流量和輸沙率以及分流比等水文要素的變化，容易造成洪水災(zāi)害以及河段沖淤逆轉(zhuǎn)等，從而直接影響河口三角洲社會(huì)經(jīng)濟(jì)的可持續(xù)發(fā)展。在珠江三角洲，強(qiáng)烈人類活動(dòng)尤為顯著，例如中上游水庫建設(shè)，森林砍伐，水土保持，河網(wǎng)采砂，口門圍墾以及河道疏浚等，都會(huì)造成河口三角洲區(qū)域水沙的異變，影響大灣區(qū)的水安全。因此各個(gè)相關(guān)部門都極為重視珠江流域以及河口三角洲的水沙變化及其原因，特別是自然和人類活動(dòng)貢獻(xiàn)的定量化，以及水沙發(fā)生異變的時(shí)間節(jié)點(diǎn)以及控制因素。

目前的研究手段基本可以定量區(qū)分自然因素和人類活動(dòng)因素對(duì)于水沙變化的影響，認(rèn)為目前人類活動(dòng)對(duì)于水沙變化的貢獻(xiàn)率在80%以上。例如在對(duì)珠江流域1954—2018年的水文氣象數(shù)據(jù)以及高要、石角和博羅站水沙數(shù)據(jù)分析的基礎(chǔ)上，發(fā)現(xiàn)人類活動(dòng)特別是大壩建設(shè)造成了96%的輸沙量變化，而4%是由于氣候變化[1]。在對(duì)1950—2000年的氣候變化和人類活動(dòng)對(duì)珠江流域(西北東江)徑流和泥沙進(jìn)行定量計(jì)算中，發(fā)現(xiàn)人類活動(dòng)對(duì)三江的貢獻(xiàn)值不同，貢獻(xiàn)值從高到低依次為西東北江，其中西江的人類活動(dòng)貢獻(xiàn)值最大為96%[2]。對(duì)1956—2013年珠江流域(西北東江)輸沙率進(jìn)行研究得出2000年后人類活動(dòng)對(duì)輸沙率減少作用更加顯著，貢獻(xiàn)率達(dá)到80%[3]。通過對(duì)高要(西江)、石角(北江)和博羅(東江)水文站1954—2011年的徑流和輸沙率的相關(guān)關(guān)系計(jì)算得出氣候變化和人類活動(dòng)對(duì)輸沙率的貢獻(xiàn)值，認(rèn)為人類活動(dòng)的貢獻(xiàn)值從30%增加到80%[4]。

對(duì)于珠江流域以及河口三角洲水沙的階段性變化以及影響因素也有了一定的認(rèn)識(shí)，例如一些學(xué)者在對(duì)西江流域(梧州站和高要站)的輸沙進(jìn)行研究時(shí)發(fā)現(xiàn)，輸沙變化的主控因素是水土保持和水利建設(shè)等人類活動(dòng)[5]，如西江流域(梧州站)輸沙量可以分為兩個(gè)階段，其發(fā)生階段性變化主要受流域大型水利樞紐工程的建設(shè)攔沙等控制因素影響[6]。另外，有研究認(rèn)為珠江流域(西北東江)輸沙率分為3個(gè)階段，人類活動(dòng)(大壩建設(shè))是造成輸沙率變化的主要原因[7]。在對(duì)東江流域(博羅站)的輸沙量的研究中發(fā)現(xiàn)，其輸沙變化分為3個(gè)階段，其控制因素主要為水庫和采砂等人類活動(dòng)[8]。三水、馬口站泥沙變化分為兩個(gè)階段，其變化的控制因素主要為上游大型水庫建設(shè)和森林覆蓋率的提高[9]。

對(duì)于珠江流域的水沙變化的研究還有很多，但大都集中在珠江流域和河口三角洲區(qū)域的水沙階段性變化趨勢(shì)研究。雖然也有一些研究指明水沙的異變，但是只是在某一個(gè)站位，例如只談了西江[10]，或者馬口三水等[11]。由于珠江流域和河網(wǎng)的復(fù)雜性，西江、北江和東江流域的氣候以及人類活動(dòng)程度不同，目前還沒有給出每一個(gè)江的具體水沙異變時(shí)間節(jié)點(diǎn)以及明確的異變機(jī)制。主要的原因是研究者使用的數(shù)據(jù)時(shí)間長度序列不一樣，或者研究的水文站位比較單一，空間分布不足，分析方法也不相同，從而造成對(duì)于水沙數(shù)據(jù)分段以及異變研究結(jié)果的不一致性。本文收集珠江流域下游5個(gè)典型水文站的1960—2019年最新的水沙數(shù)據(jù)，利用水沙突變檢測和雙累積曲線的方法，對(duì)珠江流域和網(wǎng)河區(qū)的水沙情勢(shì)異變分別進(jìn)行判斷，并且使用異變分級(jí)的方法，給出三江不同流域的水沙異變點(diǎn)和顯著變化點(diǎn)，并闡明異變的主控因子。

1 試驗(yàn)數(shù)據(jù)與方法

1.1 水沙數(shù)據(jù)來源

本文收集了高要、石角、博羅、馬口和三水5個(gè)典型水文站(圖1)1960—2019年年度流量與輸沙率數(shù)據(jù)，數(shù)據(jù)來源于廣東省水文局(歷史測量數(shù)據(jù))，中山大學(xué)(歷史保存數(shù)據(jù))、《中國河流泥沙公報(bào)》(2003—2019版) (http:∥www.mwr.gov.cn/)。統(tǒng)計(jì)年份見表1。珠江流域主要大壩建設(shè)數(shù)據(jù)來源于中國大壩工程學(xué)會(huì)網(wǎng)站(http:∥www.chincold.org.cn/)。

表1 五大水文站歷史水沙統(tǒng)計(jì)年份

1.2 分析方法

1.2.1 Mann-Kendall檢驗(yàn)法 M-K趨勢(shì)分析法(Mann-Kendall檢驗(yàn)法)是一種非參數(shù)統(tǒng)計(jì)檢驗(yàn)方法，其優(yōu)點(diǎn)是樣本不需要遵循一定的分布，也不受少數(shù)異常值的干擾，更適用于類型變量與順序變量，計(jì)算也比較簡便，最初由曼和肯德爾提出原理并發(fā)展了這一方法。現(xiàn)已成為世界氣象組織推薦并已廣泛使用的非參數(shù)檢驗(yàn)方法，適用于水文、氣象等非正態(tài)分布的數(shù)據(jù)。其計(jì)算公式如公式(1—4) ，設(shè)長時(shí)間序列(流量、輸沙率)為x1，x2，…，xn;Sk表示第i個(gè)樣本xi>xj(1≤j≤i)的累計(jì)數(shù)：

圖1 珠江流域下游5個(gè)典型水文站和上游水庫位置

(1)

其中：

(2)

在時(shí)間序列隨機(jī)獨(dú)立的假設(shè)下，Sk的均值與方差分別為：

(3)

將Sk標(biāo)準(zhǔn)化為：

(4)

式中：UF1=0,給定顯著性水平α,若|UF|>Uα,表明序列存在明顯的趨勢(shì)變化。將此方法引用到反序列,序列順序變?yōu)閤n,xn-1,…,x1,反序列標(biāo)準(zhǔn)化由UBk表示,UBk,UFk可組成兩條曲線,在這里臨界值U0.05為±1.96,設(shè)定α值為0.05。當(dāng)UBk,UFk兩條曲線出現(xiàn)交點(diǎn),則交點(diǎn)處可能是變化點(diǎn)[11]。

1.2.2 Pettitt檢驗(yàn)法 Pettitt檢驗(yàn)法依據(jù)長時(shí)間序列趨勢(shì)性的變化為基礎(chǔ)，確定序列突變點(diǎn)。假設(shè)序長度為T，變化點(diǎn)為t，把假設(shè)序列劃分為前后兩段，即兩樣本為x1,x2,…,xt和xt+1,xt+2,…,xT, Pettitt突變檢驗(yàn)的本質(zhì)是檢驗(yàn)序列中的兩個(gè)樣本是否來自于同一個(gè)樣本的Mann-Whitney統(tǒng)計(jì)量。

統(tǒng)計(jì)量Ut的計(jì)算公式為：

(5)

if(Xi-Xj)>0 sgn(Xi-Xj)=1

if(Xi-Xj)=0 sgn(Xi-Xj)=0

if(Xi-Xj)<0 sgn(Xi-Xj)=-1

(6)

在檢驗(yàn)序列時(shí)，根據(jù)順序統(tǒng)計(jì)量理論引用了K-S兩樣本檢驗(yàn)，檢驗(yàn)兩個(gè)樣本累積分布函數(shù)的最大差值而求得的累積概率為P，即突變點(diǎn)位置。

Kt=max1≤t≤T|Ut|

(7)

(8)

如果累積概率越接近于1，說明存在突變點(diǎn)的趨勢(shì)越顯著；令α為置信度，若p>Pα，則突變點(diǎn)趨勢(shì)顯著。當(dāng)α<0.5時(shí)，表示存在有效突變點(diǎn)；當(dāng)α<0.05時(shí)，表示存在顯著突變點(diǎn)；當(dāng)α<0.01時(shí)，表示存在極顯著突變點(diǎn)[12-13]。

1.2.3 判別水沙異變的程度 本研究根據(jù)Mann-Kendall檢驗(yàn)和Pettitt檢驗(yàn)的結(jié)果，對(duì)水沙變化程度進(jìn)行層級(jí)劃分。(1) 當(dāng)Mann-Kendall檢驗(yàn)的UF曲線超過置信度區(qū)間且在0.05顯著水平內(nèi)與UB曲線交點(diǎn)，此時(shí)Pettitt檢驗(yàn)曲線超過0.01顯著水平線時(shí)定義為異變點(diǎn); (2) 當(dāng)Mann-Kendall檢驗(yàn)的UF曲線在置信度區(qū)間且在0.05顯著水平內(nèi)與UB曲線交點(diǎn)，同時(shí)Pettitt檢驗(yàn)曲線僅超過0.5顯著水平線時(shí)定義為顯著變化點(diǎn); (3) 當(dāng)Mann-Kendall檢驗(yàn)的UF曲線在置信度區(qū)間且在0.05顯著水平內(nèi)與UB曲線交點(diǎn)，同時(shí)Pettitt檢驗(yàn)曲線未超過0.5顯著水平線時(shí)只表示為一般趨勢(shì)性變化點(diǎn)。本文結(jié)合水沙比(流量與輸沙率的比值，m3/kg)和水沙雙累積曲線的變化趨勢(shì)，重點(diǎn)分析水沙的異變點(diǎn)和顯著變化點(diǎn)。

2 結(jié)果與分析

根據(jù)Mann-Kendall檢驗(yàn)、Pettitt檢驗(yàn)以及雙累積曲線的方法，對(duì)西江的高要站和馬口站、北江的石角和三水站，以及東江的博羅站進(jìn)行了異變點(diǎn)和顯著變化點(diǎn)的確定(圖2，表2)。異變點(diǎn)表明此點(diǎn)前后水沙發(fā)生了異變，顯著變化點(diǎn)前后表明水沙發(fā)生了顯著的變化，但不是異變(表2)。

2.1 西江高要站和干流馬口站水沙異變的主控因子

根據(jù)西江高要站徑流量的Mann-Kendall檢驗(yàn)和Pettitt檢驗(yàn)，Mann-Kendall檢驗(yàn)中UF曲線在置信度區(qū)間且在0.05顯著水平內(nèi)與UB曲線有交點(diǎn)，同時(shí)Pettitt檢驗(yàn)曲線未超過0.5顯著水平線，因此可以認(rèn)為高要站的流量沒有異變點(diǎn)(圖2)。多年的徑流量變化范圍為3 367～10 198 m3/s，平均值為6 912 m3/s(圖3，表3)。西江高要站輸沙率Mann-Kendall檢驗(yàn)顯示UF曲線超過置信度區(qū)間，且在0.05顯著水平內(nèi)與UB曲線交點(diǎn)，表明其異變點(diǎn)為2001年。Pettitt高要站輸沙率檢驗(yàn)曲線超過0.01顯著水平線，確定的異變點(diǎn)為2000年。結(jié)合輸沙率過程圖發(fā)現(xiàn)在2000年后輸沙率開始進(jìn)入一個(gè)新階段，因此判斷2000年為西江水沙異變點(diǎn)(圖3，表3)。異變發(fā)生之前，1960—2000年的輸沙率變化范圍為521～4 167 kg/s，平均值為2 249 kg/s，水沙比為3.4 m3/kg(表3)。異變之后2000—2019年的輸沙率變化范圍為255～1 823 kg/s，平均值為764 kg/s，水沙比為11 m3/kg。異變之后輸沙率約為異變前的三分之一。說明輸沙率在量級(jí)上已經(jīng)發(fā)生了根本的變化。從雙累積曲線也可看出在2000年前后斜率發(fā)生明顯變化，斜率從0.33變?yōu)?.10(表3)。

高要站在2000年前后發(fā)生異變的主要原因是西江上游水庫的建設(shè)。1960—2000年，輸沙率變化不大，主要受到氣候的控制。在1960—1985年，氣候和人類活動(dòng)貢獻(xiàn)值分別為83%和17%[2]，此時(shí)輸沙率由于森林大量砍伐導(dǎo)致流域水土流失加劇而增大。此時(shí)廣西的水土流失面積在20世紀(jì)50—60年代為12 000 km2，而80年代為30 600 km2，是50—60年代的近3倍。而同期廣東省水土流失面積則從7 444 km2增加到了 17 070 km2，是原來的約2.3倍[14-15]。由于西江流域在廣西境內(nèi)集水面積共計(jì)20.24萬km2，占全流域集水面積的85.7%，因此廣西壯族自治區(qū)和廣東省的水土流失面積增加，使得高要站泥沙含量呈現(xiàn)上升趨勢(shì)。

而在1985—2000年，人類活動(dòng)開始增強(qiáng)，但輸沙率仍然以氣候控制為主。氣候和人類活動(dòng)貢獻(xiàn)值分別為63%和37%[2]。此階段水土流失程度減弱，廣西壯族自治區(qū)20世紀(jì)90年代的水土流失面積為28 100 km2，比80年代下降了8.2%，變化不大[14-15]。相反此時(shí)水庫大壩的影響增強(qiáng)，1985—2000年，西江上游主要修建的水庫為魯布革水庫、巖灘水庫和天生橋水庫(表4，圖1)。1992年巖灘水庫建成后高要站輸沙率則由1981—1991年的7 658萬t/a降到了1992—2002年的6 310萬t/a，降低了18%[16-17]。巖灘水庫的泥沙沉積速率為35 Mt/a，高要站的泥沙量減少了10 Mt/a[18]。因此可以看出巖灘水庫對(duì)于泥沙減少起到重要的調(diào)節(jié)作用[17]。不過巖灘水庫庫容較少，因此對(duì)于泥沙的調(diào)節(jié)作用仍然有限，沒有在1992年前后引起輸沙率的異變。

2000—2019年，雖然前期巖灘水庫調(diào)節(jié)作用有限，但在1997年建成的天生橋水庫的庫容量是巖灘水庫的3倍多(表4)，建成后經(jīng)過一定時(shí)間的蓄水和攔沙過程，導(dǎo)致2000年后的輸沙率大幅下降，致使輸沙率發(fā)生異變，之后在水庫的調(diào)節(jié)作用下，輸沙率保持穩(wěn)定波動(dòng)。其中天峨站(巖灘水庫上游150 km處的一個(gè)水文站)在天生橋水庫建成后的輸沙量減少了30 Mt/a，發(fā)生了較大的變化，也可以幫助說明天生橋水庫的巨大調(diào)節(jié)作用[17]。值得注意的是2006年在西江上游建設(shè)了更大型的龍灘水庫(表4，圖1)，其庫容量約為天生橋水庫的2.5倍。水庫建成后使得2006—2019年輸沙減小之后保持基本穩(wěn)定(圖3)，但是并沒有引起輸沙率的異變。這是因?yàn)樘焐鷺蛩畮煳挥谖鹘嫌?圖1)，而龍灘水庫位于天生橋水庫下游，天生橋水庫已經(jīng)攔截大部分泥沙，從而導(dǎo)致在2000年輸沙率發(fā)生異變，而龍灘水庫庫容量大，但建成后對(duì)下游的輸沙影響不如天生橋水庫，沒有引起輸沙率發(fā)生異變。

圖2 1960-2019年珠江流域五站位年流量與年輸沙率MK檢驗(yàn)和Pettitt檢驗(yàn)

表2 珠江流域五站點(diǎn)異變的判定

輸沙率異變后森林覆蓋面積持續(xù)上升，例如廣西壯族自治區(qū)森林覆蓋率由2000年的39%增長到2005年的53%[11]，因此水土保持使得輸沙率持續(xù)小幅度降低。另外，1995—2014年高要站附近也發(fā)生了挖沙的行為[19]，也對(duì)輸沙率減少起到一定的作用。此階段人類活動(dòng)的貢獻(xiàn)值大幅上升，氣候和人類活動(dòng)貢獻(xiàn)值分別為4%和96%[2]。

據(jù)馬口站流量的Mann-Kendall檢驗(yàn)和Pettitt檢驗(yàn)，發(fā)現(xiàn)在1986年Mann-Kendall檢驗(yàn)的UF曲線在置信度區(qū)間且在0.05顯著水平內(nèi)與UB曲線交點(diǎn)，同時(shí)Pettitt檢驗(yàn)曲線僅超過0.5顯著水平線，因此可以認(rèn)為馬口站的流量在1986年存在顯著變化點(diǎn)(圖2)。1960—1986年的流量變化范圍為3 818～9 978 m3/s，平均值為7 498 m3/s，水沙比為3.3 m3/kg，水沙變化的相關(guān)系數(shù)為0.52。雙累積曲線斜率為0.33。顯著變化之后1986—2019年的流量變化范圍為4 403～9 602 m3/s，平均值為6 777 m3/s，較異變前減少9.6%。水沙比為8.0 m3/kg，水沙變化的相關(guān)系數(shù)為0.08，雙累積曲線斜率為0.17(表3)。馬口在1986年發(fā)生的徑流顯著變化，1986年初期發(fā)生的徑流量變化和當(dāng)時(shí)劇烈的氣候動(dòng)蕩變化有關(guān)[4]。

西江馬口站輸沙率的Mann-Kendall檢驗(yàn)線在0.05顯著水平之外，顯示的異變點(diǎn)為1999年；Pettitt檢驗(yàn)線0.01顯著水平線之外，異變點(diǎn)為1999年。因此判斷1999年為西江干流馬口站水沙異變點(diǎn)(圖3)。異變發(fā)生之前，1960—1999年的輸沙率變化范圍為506～4 186 kg/s，平均值為2 306 kg/s，水沙比為3.5 m3/kg，水沙變化的相關(guān)系數(shù)為0.36(表3)。異變之后1999—2019年的輸沙率變化范圍為202～1 680 kg/s，平均值為745 kg/s，水沙比為10.8 m3/kg，水沙變化的相關(guān)系數(shù)為0.46，輸沙率約為異變前的三分之一。從雙累積曲線可看出在1999年前后斜率發(fā)生明顯變化，斜率為從0.34變?yōu)?.10(表3)。因此馬口站的水沙變化基本和高要站相似。其不同階段輸沙異變的原因也和高要站一致。因此馬口的輸沙率異變主要受到上游水庫的建設(shè)誘發(fā)。不過馬口站的輸沙率異變的時(shí)間略微早于高要站，可能和馬口站和三水站的水流分水比在1988—1997年期間的顯著變化有關(guān)[11]。馬口站的采砂活動(dòng)主要集中在1995—2006年[19]，因此采砂會(huì)對(duì)輸沙率有一定的影響，但并不是主要因素。

2.2 北江的石角站和干流三水站水沙異變和主控因子

根據(jù)徑流和輸沙率異變和顯著變化的診斷標(biāo)準(zhǔn)，和高要以及馬口站的判別方式一樣，確定北江石角站徑流量沒有異變點(diǎn)。多年的徑流量變化范圍為515～2 280 m3/s，平均值為1 335 m3/s(圖3)。石角站輸沙率沒有異變點(diǎn)，但是存在顯著變化點(diǎn)，結(jié)合輸沙率過程圖發(fā)現(xiàn)在1999年后輸沙率開始進(jìn)入一個(gè)新階段，因此判斷1999年為北江石角輸沙率顯著變化點(diǎn)(圖3，表3)。1960—1999年，輸沙率整體平穩(wěn)，呈現(xiàn)上升趨勢(shì)，輸沙率平均值為190 kg/s，水沙比為7.9 m3/kg，水沙變化的相關(guān)系數(shù)為0.68，雙累積曲線斜率為0.15(表3)。1965年之前，幾乎不受人類活動(dòng)影響，氣候貢獻(xiàn)90%，人類活動(dòng)貢獻(xiàn)10%[2]。該階段降雨量減少，造成徑流量和輸沙率的降低[2]。1965年后人類活動(dòng)逐漸增強(qiáng)。例如1971年修建的南水水庫，庫容量約為12億m3(表4)，引起了隨后的輸沙率的急劇減少(圖2)。

注：虛線為顯著變化點(diǎn)，實(shí)線為異變點(diǎn)。

表3 珠江流域五站位流量與輸沙率之間的相關(guān)關(guān)系

表4 珠江流域較大規(guī)模水庫建設(shè)情況

雖然繼續(xù)有水庫的建設(shè)，但是水庫庫容量都小于2億m3，對(duì)泥沙的攔截作用有限。同時(shí)森林砍伐引起的水土流失作用超過了水庫攔沙的影響[11]，覆蓋北江流域的廣東省水土流失面積從20世紀(jì) 50—60年代的7 444 km2增加到了80年代的17 070 km2，增加了1.3倍[14-15]，比1950s初期增加了2倍多[9]，使得河口三角洲區(qū)域輸沙率增加。而水庫的建設(shè)所攔截的泥沙量小于水土流失的量，因此泥沙含量呈現(xiàn)上升趨勢(shì)。1985—1993年珠江的降雨量、徑流量和輸沙率處于劇烈的減小階段[4]，本階段輸沙率受到氣候影響的貢獻(xiàn)大致為88%，人類活動(dòng)影響的貢獻(xiàn)為12%[2]。因而輸沙率在該階段前期呈下降趨勢(shì)，但是輸沙率在此階段后期(1989年)呈現(xiàn)一定上升，主要還是受到氣候和徑流的控制。

1999—2019年，在1999年建設(shè)的飛來峽大壩(庫容約為20億m3)，造成了石角站輸沙率的顯著變化，這一階段的平均輸沙率為134 kg/s，水沙比為12.3 m3/kg，水沙變化的相關(guān)系數(shù)為0.67，雙累積曲線斜率為0.11(表3)。因此1999年的水庫建設(shè)造成了石角輸沙率的顯著變化。本階段氣候貢獻(xiàn)為36%，人類活動(dòng)貢獻(xiàn)為64%[2]。這一階段國家水土保持政策使得森林覆蓋面積有一定增長，20世紀(jì)80年代中期開始，在國家水土保持政策的引導(dǎo)下，北江森林覆蓋率有一定的提高，北江區(qū)域面積為4.8萬 km2，森林覆蓋率由20世紀(jì)80年代3.5萬 km2增加到20世紀(jì)90年代3.7萬 km2，森林覆蓋率由73%上升為77%[7]，上升速度并不快，因此水土流失得到一定控制，對(duì)輸沙率的減少有一定的作用，但不是主控因子。另外，石角站從2003—2014年有持續(xù)的挖沙活動(dòng)，也對(duì)輸沙率的減小有一定貢獻(xiàn)。雖然水沙存在顯著性變化，但是不是異變，水沙相關(guān)性一直相對(duì)較好。

根據(jù)異變和顯著變化的診斷標(biāo)準(zhǔn)，北江干流三水站徑流量的異變點(diǎn)為1991年。異變發(fā)生之前，1960—1991年的流量變化范圍為298～1 910 m3/s，平均值為1 213 m3/s。異變之后1991—2019年的流量變化范圍為1 139～2 938 m3/s，平均值為1 952 m3/s，較異變前增加61%。水沙比為原來的2倍，水沙變化的相關(guān)系數(shù)變差，雙累積曲線斜率減小(表3)。

北江三水站異變點(diǎn)為2003年(圖2)。異變發(fā)生之前，1960—2003年輸沙率變化范圍為18～578 kg/s，平均值為292 kg/s(表3)。異變之后2003—2019年的輸沙率變化范圍為44～296 kg/s，平均值為155 kg/s，水沙比為原來的2.3倍，水沙變化的相關(guān)性變差。泥沙量約為異變前的二分之一，說明輸沙率在量級(jí)上已經(jīng)發(fā)生了根本的變化。從雙累積曲線可看出在1991年和2003年前后斜率發(fā)生顯著變化，斜率從0.21變?yōu)?.09(表3)。

雖然三水站主要的徑流來自北江，但是其異變現(xiàn)象并不和石角站一致，石角站徑流量并不存在異變點(diǎn)和顯著變化點(diǎn)。三水站1991年前徑流量基本保持穩(wěn)定，徑流量突然增加的原因是由于三水站在1989—2006年發(fā)生了大量的采砂活動(dòng)[19]，而北江河網(wǎng)區(qū)的采砂量遠(yuǎn)遠(yuǎn)大于西江網(wǎng)河區(qū)，致使北江網(wǎng)河區(qū)的河床出現(xiàn)大量下切，導(dǎo)致三水站分流比升高，更多的徑流流入到三水站，引起了徑流量的突然增大，是以前的1.6倍，因而發(fā)生異變。而此時(shí)我們發(fā)現(xiàn)馬口站并沒有同時(shí)發(fā)生徑流異變，主要是因?yàn)轳R口站的徑流量較大，例如20世紀(jì)80年代之后，馬口站年徑流量減少12%，而三水站年徑流量增加54%[11]，因此由于和三水站分水分沙引起的變化量并不足以使得馬口站產(chǎn)生徑流量的異變。

同樣，三水站的輸沙率并不和徑流量異變時(shí)間一致，而是在2003年發(fā)生異變，1960—2003年水沙變化的相關(guān)系數(shù)為0.56，異變后2003—2019年水沙相關(guān)性變差，相關(guān)系數(shù)為0.30(表3)。這是因?yàn)殡m然在三水站附近的挖沙活動(dòng)集中于1989—2006年[19]，加上1999年的飛來峽水庫，實(shí)際上已經(jīng)造成了三水站輸沙率的持續(xù)降低，但三水站的輸沙率在1999—2003年仍然保持一定的水平(圖3)。而2003年，三水站的輸沙率有了斷崖式下降，并且在2003年之后輸沙率都在較低水平震蕩，說明2003年后的異變主要是這段時(shí)間的上游來沙量減少所致。2000年左右，廣東省將水土保持面上資金納入省財(cái)政專項(xiàng)，每年固定投入1 500萬元。近幾年來，廣東省每年綜合治理水土流失面積約400 km2[20]。根據(jù)時(shí)間點(diǎn)變化，判斷上游的水土保持是三水站輸沙率發(fā)生異變的主要原因。另外2003—2004年遭遇極其干旱的天氣，降雨量只有1 423～1 315 mm，鄰近的2001年、2002年、2005年、2006年的降雨量都在1 766～2 110 mm[21]，加劇了上游來沙量的減少。

2.3 東江博羅站的水沙異變和主控因子

根據(jù)異變和顯著變化的診斷標(biāo)準(zhǔn)，東江博羅站徑流量不存在異變點(diǎn)。多年的流量變化范圍為284～1 313 m3/s，平均值為741 m3/s，水沙比為16.4 m3/kg。東江博羅利用輸沙率異變點(diǎn)為1988年(圖2)。異變發(fā)生之前，1960—1988年的輸沙率變化范圍為10～161 kg/s，平均值為87 kg/s，水沙比為9.7 m3/kg。異變之后1988—2019年的輸沙率變化范圍為7～128 kg/s，平均值為43 kg/s，水沙比為22.6 m3/kg。泥沙量約為異變前的二分之一。說明輸沙率在量級(jí)上已經(jīng)發(fā)生了根本的變化。從雙累積曲線可看出在1988年前后斜率發(fā)生顯著變化，斜率從0.11變?yōu)?.06(表3)。

輸沙率異變之前的1960—1988年，此階段輸沙率變化較為平穩(wěn)，輸沙率有略微下降。水沙變化的相關(guān)系數(shù)為0.60(表3)。此階段主要由人類活動(dòng)主導(dǎo)，人類活動(dòng)的貢獻(xiàn)為89%[2]。河流輸沙隨著流域水庫累計(jì)庫容的增加而減少。主要受到上游新豐江水庫(1962年)、楓樹壩水庫(1973年)和白盆珠水庫(1985年)的影響(表4)。3個(gè)水庫的共同作用下雖然使輸沙率減小，但是總的庫容量并不大，對(duì)于白盆珠水庫來說，其多年平均流量為37.5 m3/s，而博羅站本身徑流量為741 m3/s，遠(yuǎn)遠(yuǎn)大于其調(diào)節(jié)流量[22-23]。從輸沙率的變化過程圖來看(圖3)，輸沙率在2004—2006年較高，也可以證明水庫的調(diào)節(jié)作用有限，因此水庫建設(shè)對(duì)于博羅站輸沙率異變影響不大，不是異變的主要原因。而在1988年后，1988—2003年人類活動(dòng)的貢獻(xiàn)為從89%降為75%[2]。在國家水土保持政策引導(dǎo)，大范圍的治山治水、開發(fā)“四荒”以及退耕還林，水土流失得到較好控制。自20世紀(jì)90年代以來，共治理東江中上游水土流失區(qū)域961.3 km2，植被覆蓋率從治理前的30%提高到65%，入河泥沙則從治理前800～1 000萬t/a降到350～400萬t/a[23]。2000年左右，廣東省每年綜合治理水土流失面積約400 km2[20]。此階段并沒有大型水庫建設(shè)，根據(jù)時(shí)間點(diǎn)變化，水土保持是造成博羅站的徑流量的異變的主要控制因素。由于在博羅站附近1998—2008年階段才有大規(guī)模持續(xù)的挖沙行為[19]，特別是2002—2005年間，新角水閘—譚公廟河段(博羅站上游)河床采砂量高達(dá)1 800萬 m3，河床平均下切5.7 m；受此影響博羅站同流量級(jí)水位下降嚴(yán)重，潮汐動(dòng)力顯著增強(qiáng)[24]，但博羅站的輸沙率在1988年就已經(jīng)發(fā)生了異變，因此挖沙對(duì)博羅站輸沙率發(fā)生異變的影響很小。

3 討 論

從研究結(jié)果看，長序列的最新的水沙數(shù)據(jù)，以及對(duì)前人大量的工作成果總結(jié)，更加全面準(zhǔn)確地認(rèn)識(shí)珠江流域以及河口三角洲發(fā)生的水沙異變及其原因，基本解決了由于數(shù)據(jù)時(shí)間長度序列不一樣、站位不同以及分析方法不同造成的計(jì)算和判別結(jié)果不一致的問題。特別是西江和北江及其干流、以及東江的同步數(shù)據(jù)分析，更加清晰地理解這些地區(qū)的水沙異變的時(shí)空變化及其原因。

水沙異變主要是在徑流量和輸沙率持續(xù)變化后由于大型的事件導(dǎo)致兩者發(fā)生了數(shù)量上本質(zhì)性的改變，除了三水的徑流量發(fā)生異變后徑流量大幅上升，異變后徑流量平均值是異變前的1.6倍，其他各站徑流量保持平穩(wěn)。各站的輸沙率發(fā)生異變或者顯著變化后普遍減少，例如高要站和馬口站輸沙率發(fā)生異變后，輸沙率只有異變之前的三分之一。三水和博羅站發(fā)生異變后，輸沙率只有異變之前的二分之一。

近70 a來，珠江流域和河口三角洲的徑流量除了個(gè)別站位，基本仍然受到氣候的控制，而泥沙則受到人類活動(dòng)的控制，上游水庫建設(shè)、水土保持、以及河網(wǎng)區(qū)挖沙仍然是影響不同流域河流泥沙異變的主要原因。不過在珠江流域及河口三角洲區(qū)域，三江及其干流的水沙異變情勢(shì)并不相同。在西江高要、北江石角、東江博羅，以及西江干流馬口和北江干流三水，只有三水站的徑流在1991年發(fā)生異變，主要是由于人類挖沙導(dǎo)致三水河槽剖面迅速下降，經(jīng)過思賢滘的分水比增大引起了徑流量異變。這和前人根據(jù)1954—2011年徑流量和輸沙率的數(shù)據(jù)確定三水站徑流1990年發(fā)生突變相差不大[11]。但是前人認(rèn)為馬口站在2002年前后徑流量減少12%，發(fā)生異變，但本文根據(jù)判別標(biāo)準(zhǔn)，認(rèn)為這種量級(jí)的變化并不構(gòu)成異變，因此馬口站徑流只有顯著變化，并沒有發(fā)生異變，主要是因?yàn)轳R口站的徑流量大，減少的水量不至于引起其異變。

5個(gè)站位有4個(gè)站位的輸沙率量發(fā)生了異變，而每個(gè)站位的異變時(shí)間和原因并不相同。高要站和馬口站在2000年和1999年發(fā)生了輸沙率的異變，主控制因子為上游水庫和大壩的建設(shè)。這與前人的研究結(jié)果相差不大。例如在西江流域(梧州站)以1960—2015年年輸沙量為研究對(duì)象，采用R/S法和Mann-Kendall檢驗(yàn)法分析確定2000年為異變點(diǎn)[6]。對(duì)西江高要站1957—2016年的輸沙數(shù)據(jù)通過滑動(dòng)T檢驗(yàn)、有序聚類法和累積距平法綜合判斷西江高要站輸沙量的異變點(diǎn)為1999年[25]。

北江干流的三水站和東江博羅站的輸沙率分別在2003年和1988年發(fā)生異變，主要是由于三水站附近挖沙以及東江上游水土保持造成的。只有石角站沒有發(fā)生異變，但是在1999年發(fā)生了顯著變化，1999年的顯著減小是由于上游大壩的建設(shè)引起的。前人根據(jù)1957—2000年輸沙量計(jì)算確立了石角站1982年為輸沙量的明顯變化點(diǎn)[25]，而1983年和1993年為兩個(gè)水沙變化節(jié)點(diǎn)[5]，認(rèn)為其主要控制因素是采砂活動(dòng)和水庫建設(shè)，和本文的研究結(jié)果有所差異，這些節(jié)點(diǎn)只能是一般趨勢(shì)性變化節(jié)點(diǎn)。前人研究表明博羅站輸沙量(1954—2011年)在20世紀(jì)80年代以來下降是一突變現(xiàn)象，并確定1984年前后輸沙量發(fā)生異變，異變后輸沙量減少46%[11]，這與本文結(jié)果相差不大。

另外，馬口三水的徑流量和輸沙率的異變和其分流分沙比的異變并不相同。前人研究表明由于北江下游河道下切幅度大于西江河道下切幅度，使馬口三水站分流比于1988年發(fā)生異變、分沙比1992年發(fā)生異變[11]。不過三水的徑流量異變發(fā)生在1991年，而輸沙率異變發(fā)生在2003年。因此徑流量發(fā)生異變是在馬口三水的分流比異變后發(fā)生，主要和挖沙導(dǎo)致的河道下切有關(guān)。但是輸沙率在2003年發(fā)生異變，卻和分沙比1992年的異變沒有關(guān)聯(lián)，因?yàn)槿据斏陈实淖兓饕艿缴嫌嗡帘３值膭×易兓嘘P(guān)，其分沙比引起的輸沙率變化還不足以引起異變。

4 結(jié) 論

(1) 根據(jù)Mann-Kendall檢驗(yàn)法和Pettitt檢驗(yàn)法兩種診斷方法，將水沙變化分為異變和顯著變化兩個(gè)層級(jí)。當(dāng)Mann-Kendall檢驗(yàn)的UF曲線超過置信度區(qū)間且在0.05顯著水平內(nèi)與UB曲線交點(diǎn)，此時(shí)Pettitt檢驗(yàn)曲線超過0.01顯著水平線時(shí)定義為異變點(diǎn)。而當(dāng)Mann-Kendall檢驗(yàn)的UF曲線在置信度區(qū)間且在0.05顯著水平內(nèi)與UB曲線交點(diǎn)，同時(shí)Pettitt檢驗(yàn)曲線僅超過0.5顯著水平線時(shí)定義為顯著變化點(diǎn)。

(2) 對(duì)于流量來說，馬口站流量在1986年發(fā)生顯著變化，主要受氣候控制；顯著變化前水沙比為3.3 m3/kg，顯著變化后水沙比為8.0 m3/kg，水沙相關(guān)系數(shù)從0.52變?yōu)?.08；雙累積曲線斜率也從0.33變?yōu)?.17。三水站的流量在1991年發(fā)生異變，異變后徑流量較異變前增加61%，主要是由于人類挖沙導(dǎo)致三水剖面迅速下降，經(jīng)過思賢滘的分水比增大引起了徑流量異變。異變前水沙比為5.2 m3/kg，異變后水沙比為10.5 m3/kg，水沙相關(guān)系數(shù)從0.61變?yōu)?.44；雙累積曲線斜率也從0.22變?yōu)?.12。其他站位沒有發(fā)生徑流異變以及顯著變化；

(3) 5個(gè)站位的輸沙量只有石角站沒有發(fā)生異變，但是在1999年發(fā)生了顯著變化；石角站發(fā)生輸沙率顯著變化后，較第一階段減少29.5%。水沙比從7.9 m3/kg變?yōu)?2.3 m3/kg；雙累積曲線斜率從0.15變?yōu)?.11。受到上游水庫的建設(shè)控制。

(4) 除石角外，其余4個(gè)站位的輸沙率量都發(fā)生了異變，而每個(gè)站位的異變時(shí)間和原因并不相同。高要站和馬口站在1999—2000年發(fā)生了輸沙率的異變，輸沙率只有異變之前的三分之一，異變前后兩站水沙比變化幅度相似，從3.4～3.5 m3/kg異變?yōu)?1.0～10.8 m3/kg；雙累積曲線斜率兩站變化幅度相同從0.34變?yōu)?.10，主要受到上游水庫和大壩的建設(shè)的控制。北江干流的三水站和東江博羅站的輸沙率分別在2003年和1988年發(fā)生異變，主要是由于上游水土保持造成的。異變后輸沙率只有異變之前的二分之一，異變前后兩站水沙比變化幅度相似，從5.6～9.7 m3/kg異變?yōu)?3.4～22.6 m3/kg；雙累積曲線斜率也從0.11～0.21變?yōu)?.06～0.09。