榕江流域多年水沙特性分析

敖亮挺

(廣東省水利電力勘測設(shè)計研究院有限公司，廣東 廣州 510635)

榕江俗稱南河，又稱揭陽江，位于廣東省東部，瀕臨南海，南海水系外河流，其東北和西南方向分布有韓江、練江。榕江是粵東地區(qū)的第二大河流，僅次于韓江，也是廣東省著名深水河，僅次于珠江，可進(jìn)出3 000～5 000 t級貨輪，直航香港和廣州、上海、湛江等地，歷史上有“黃金水道”和“狀元港”的美譽。榕江上游由南、北兩河組成，于砲臺鎮(zhèn)雙溪咀匯合形成榕江干流，經(jīng)汕頭港內(nèi)的牛田洋注入南海。榕江全長175 km，多年平均徑流量61.4億m3，河床平均坡降0.49‰。其流域面積4 408 km2，山區(qū)面積占47.8%，丘陵占16.2%，平原占36%。

南河為榕江主流，發(fā)源于汕尾市陸河縣鳳凰山南麓，發(fā)展到揭陽市三洲攔河閘以下進(jìn)入感潮區(qū)，每天出現(xiàn)兩次高潮和兩次低潮，相鄰兩次高潮或低潮的潮位不等，漲落潮歷時也不等，屬不規(guī)則半日潮。北河是榕江最大的支流，發(fā)源于梅州市豐順縣桐子洋，發(fā)展到北河橋閘以下為感潮河段。榕江南河上游河道坡降較陡,洪水傳播迅速，而中游河槽平衍淤淺，沙洲礙流，河道彎曲，易受潮汐沙頂托，洪水宣泄不暢，下游則河面寬，吹程大，風(fēng)暴潮威脅嚴(yán)重。20世紀(jì)90年代以前河槽淤淺嚴(yán)重,沙洲不斷增多。榕江北河河道彎曲狹窄，坡降平緩，水土流失嚴(yán)重。1957—1959年，由于新嶺礦場采砂沖洗泥沙大量入河，河槽淤積較為嚴(yán)重，1960年以后浚河固堤，河道有所疏通。1970年大洪水，新西河水庫泄洪道二、三級跌水被沖垮，大量泥沙流入北河，致使赤坎河段淤高1 m左右。榕江河段洪水期受暴雨影響來水較豐，枯水期基本沒有徑流匯入，由于受到潮流影響，泥沙運動特點較為復(fù)雜[1]。榕江流域歷年來水土流失比較嚴(yán)重，河道沖淤變化明顯，并且水流受到潮汐的影響較大，漲落急時刻流速可達(dá)1.5 m/s[2]。在現(xiàn)代經(jīng)濟快速發(fā)展的過程中，人類不合理的開發(fā)、利用自然資源和工業(yè)無限發(fā)展對生態(tài)環(huán)境造成的破壞，導(dǎo)致了河道污染、水質(zhì)下降、水資源短缺、水土流失等系列系統(tǒng)性問題。地表植被的水源涵養(yǎng)力下降，氣候和季節(jié)的突變，生態(tài)重構(gòu)對水域變異等問題都迫切地等待人們?nèi)ゲ粩嘌芯亢吞接憽?/p>

1 研究資料與方法

本文研究區(qū)域見圖1，南河至東橋園水文站，北河至赤坎水文站，兩河于雙溪咀匯合，其中南河為榕江主流，為主要研究區(qū)域。徑流量分析資料取自榕江南河?xùn)|橋園站1953—2017年的流量數(shù)據(jù)，榕江北河赤坎站1967—2017年的流量數(shù)據(jù)；泥沙資料取自榕江南河?xùn)|橋園站1955—2017年的輸沙量與含沙量數(shù)據(jù)，北河流量資料缺少而與南河流量資料不同步，因此僅可作為輔助分析。因北河與南河水土流失情況近似，故將歷年泥沙演變與南河作同等程度對待。

圖1 研究區(qū)域地理位置

趨勢性、突變性和周期性是水文時間序列的重要特征[3-4]。本文采用5 a滑動平均法、Mann-Kendall檢驗等方法對水沙序列的趨勢性與突變性進(jìn)行分析，利用Morlet小波分析進(jìn)行周期性分析。下文對這些方法的主要特點和適用性作簡要的介紹。

1.1 滑動平均法

滑動平均法是以遞推的形式，不斷滑動取一定量的相鄰數(shù)據(jù)進(jìn)行算術(shù)平均[5]，以削弱序列中的短周期波動，能有效抑制頻繁波動起伏的隨機誤差，最后得到的滑動平均值隨時間的變化可反映序列的變化趨勢[6]。滑動平均法的主要特點是簡單便捷，但存在一定的主觀性和隨意性。

1.2 Mann-Kendall檢驗法

Mann-Kendall檢驗法是一種非參數(shù)檢驗方法，該方法適用范圍廣、避免人為干擾、高度定量化，是水文要素評估常會用到的方法。該方法可以檢測分析時間序列的變化趨勢，并且可以分析出序列發(fā)生突變的位置[7]。M-K方法采用原時間序列得到秩序列，由秩序列計算得到隨時間變化的統(tǒng)計量過程線UF與UB，其交點通常為突變點。當(dāng)UF與UB過程線超過臨界線時，表明變化趨勢顯著，并根據(jù)UF與UB的正負(fù)情況判斷序列的升降[6]；當(dāng)其交點位于臨界線外時，表明時間序列顯著變異。M-K檢驗法的顯著水平一般取為α=0.05，此時臨界值為±1.96。

1.3 有序聚類分析法

有序聚類分析法是通過遞推的方式設(shè)分割點，以此點將序列分為兩部分，求得突變點前后序列的離差平方和，再求出總的離差平方和S，S最小時的分割點為最優(yōu)分割點，即為序列的突變點[8]。有序聚類分析法簡易直觀，計算結(jié)果較精確。有序聚類分析法識別突變點簡單有效，但在水文序列的突變點識別過程中可能產(chǎn)失漏和偏差而導(dǎo)致計算結(jié)果不準(zhǔn)確。

1.4 小波分析法

小波分析方法源于對信號分析的伸縮與平移[9]，在水文中主要用于具有多時間尺度和振蕩特征的時間序列，識別其時間-頻率尺度，能夠清楚的顯示出隱含在時間序列中隨時間變化的尺度周期，也能反映時間序列的變化趨勢，并對未來的演變趨勢作出定性估計[10]。另外，由小波系數(shù)得到的小波方差圖也可反映出該序列演變過程的主周期。小波分析法通過小波變換系數(shù)，能快速獲得任意長度的水沙序列的變化周期，且計算結(jié)果精確全面，適用性好。

2 研究結(jié)果

2.1 榕江流域水沙分布特性

根據(jù)東橋園水文站流量資料(1953—2017年)統(tǒng)計，榕江南河多年平均徑流量27.62億m3，多年平均流量87.44 m3/s；1961年對應(yīng)年徑流量和年均流量最大，分別為48.60億m3和154.12 m3/s；2004年對應(yīng)年徑流量和年均流量最小，分別為12.04億m3和38.17 m3/s。

根據(jù)赤坎水文站流量資料(1968—2017年)統(tǒng)計，榕江北河多年平均徑流量為8.17億m3，多年平均流量26.10 m3/s；1983年對應(yīng)年徑流量和年均流量最大，分別為12.99億m3和41 m3/s；2004年對應(yīng)年徑流量和年均流量最小，分別為4.83億m3和15.40 m3/s。

根據(jù)東橋園水文站泥沙資料(1955—2017年)統(tǒng)計，榕江南河多年平均含沙量為0.157 kg/m3，最大年平均含沙量為0.321 kg/m3(1959年)，最小年平均含沙量為0.034 kg/m3(2014、2017年)。多年平均輸沙率為14.543 kg/s，最大年均輸沙率為37.8 kg/s(1973年)，最小年均輸沙率為1.288 kg/s(2004年)。因北河與南河水土流失情況近似，故可作同等程度對待。

2.2 徑流泥沙的趨勢性分析

趨勢性分析主要是用于分析時間序列總的走勢，結(jié)論分別僅有增加、減少或不變3種[11]。趨勢分析對預(yù)測未來流域演變的趨勢有著重要意義。榕江流域多年徑流量過程線見圖2a、2b，采用5 a滑動平均輔助分析，由圖可見，近60年來的南、北河徑流量整體上均呈現(xiàn)為年際波動的交替曲線，5 a滑動平均值在線性趨勢附近上下波動，表明徑流量無明顯增加或減少的趨勢。

a)南河徑流量

多年來沙統(tǒng)計過程線見圖2c、2d，由圖可見，近60年來的輸沙量及含沙量均呈年際波動，且隨時間推移均呈整體下降趨勢；2005—2013年期間，兩者的年際波動現(xiàn)象較為劇烈，2013年后較為穩(wěn)定。對來沙資料進(jìn)行線性擬合，得到輸沙量以0.914萬t/a的趨勢減小，含沙量以0.003 16 kg/(m3·a)的趨勢減小。

采用Mann-Kendall趨勢檢驗法進(jìn)一步對水沙資料進(jìn)行趨勢性分析，計算結(jié)果見表1，由表可見，東橋園站1953—2017年徑流量的統(tǒng)計量為0.47，赤坎站1967—2017年徑流量的統(tǒng)計量為0.655，兩站統(tǒng)計量|Z|均小于顯著性水平α=0.05時的臨界值1.96，即表明南河與北河的多年徑流量均有增加的趨勢，但趨勢不顯著；輸沙量與含沙量的統(tǒng)計量|Z|均大于顯著性水平，表明輸沙量與含沙量均有顯著減小的趨勢。

表1 榕江流域徑流量、輸沙量與含沙量逐年變化趨勢M-K檢驗結(jié)果

2.3 徑流泥沙的突變性分析

2.3.1泥沙序列突變分析

水文要素的突變性是指：突變位置的前后序列的變化趨勢不同[12]。即時間序列從一個均值狀態(tài)或者趨勢突然跳躍到另一個均值狀態(tài)或趨勢[13]。據(jù)此對造成突變的自然因素與人為因素進(jìn)行分析，有利于對將來流域的管理治理。

由圖2及表1可知，南河與北河多年徑流量無明顯趨勢性變化，時間序列不存在顯著突變，而泥沙序列存在明顯的減小趨勢，因此，下文僅對輸沙量與含沙量進(jìn)行顯著性突變分析。在分析榕江流域近60年來水沙趨勢的基礎(chǔ)上，通過Mann-Kendall突變檢驗法確定榕江流域歷史輸沙量與含沙量的突變時間。

圖3a為輸沙量的M-K突變檢驗統(tǒng)計量的曲線。圖中可看出在0.05顯著水平內(nèi)，輸沙量的統(tǒng)計曲線UF與UB曲線在1961、1975—1980年出現(xiàn)交點，而在顯著區(qū)間外無交點，表明輸沙量的突變點可能在1961年或1975—1980年，但均不顯著。含沙量的統(tǒng)計曲線見圖3b，可看出在1960、1968—1971、1980—1987、2013—2015年UF與UB曲線均出現(xiàn)交點，而1980—1987年、2013—2015年的交點在顯著區(qū)間外，為顯著性突變可能發(fā)生的時間段。

a)輸沙量

以上僅僅確定了泥沙序列突變可能發(fā)生的時間段，為進(jìn)一步確定輸沙量與含沙量的突變點，采用有序聚類分析法對輸沙量與含沙量序列進(jìn)行的突變檢驗。檢驗結(jié)果見圖4，輸沙量與含沙量離差平方和的最小值均在1987年，表明榕江流域的輸沙量與含沙量可能在1987年發(fā)生了突變。

a)輸沙量

為進(jìn)一步直觀展示輸沙量與含沙量的突變情況，采用累計距平法對泥沙序列的變化規(guī)律進(jìn)行分析。圖5為輸沙量與含沙量的累計距平曲線，結(jié)合圖中信息可知：輸沙量與含沙量的累計距平曲線趨勢大致相同，均呈先升后降，其由升轉(zhuǎn)為降的年份出現(xiàn)在1987年，表明在1987年前，榕江流域整體處于多沙期，1987—2017年間整個階段處于少沙期。

圖5 輸沙量與含沙量累計距平曲線

因此，在結(jié)合多種突變分析方法的基礎(chǔ)上，得到榕江流域的輸沙量在1987年變化不大，而含沙量在1987年變化顯著。

2.3.2沙量突變成因分析

a)持續(xù)枯水年的作用。根據(jù)東橋園站徑流資料作差積曲線分析，見圖6，東橋園站與赤坎站呈現(xiàn)豐枯水年交替變化，而20世紀(jì)80年代末開始，接近20年幾乎為連續(xù)枯水時間段，河段的輸沙量與含沙量因此受到一定的影響。

圖6 東橋園站年徑流差積曲線

b)建筑工程的大量采砂作用。20世紀(jì)90年代前，榕江南河河槽淤淺嚴(yán)重,沙洲不斷增多；榕江北河河道彎曲狹窄，坡降平緩，水土流失嚴(yán)重。20世紀(jì)80年代后期，榕江流域的整治疏浚開始受到重視，20世紀(jì)90年代后，由于建筑工程的大量采砂，河槽開始逐年下切，沙洲面積減小，上游來沙落淤至深槽內(nèi)，使得流域的輸沙量與含沙量均減小。

c)航道疏浚的作用。榕江流域?qū)儆谄皆貛В篮橐筝^高，航道整治通常采用挖槽進(jìn)行[14]。20世紀(jì)90年代初，為滿足運移量不斷增長的需求，榕江流域進(jìn)行了長距離的航道疏浚，榕江河道水流的挾沙能力隨著河道水深的增大和流速的減小而減小，泥沙受水流挾帶而運動，水流的挾沙能力與流速成正比。

2.4 徑流泥沙周期性分析

水文序列的多年變化周期十分復(fù)雜，提取出來的一般為近似周期。周期性跟隨研究尺度的不同而發(fā)生相應(yīng)的變化，這種變化一般表現(xiàn)為大尺度的變化周期嵌套著小時間尺度的變化周期[15-16]。得到周期性規(guī)律能夠更好地預(yù)測序列未來的變化規(guī)律與未來年份的水文特征，而且可以為流域的防洪與灌溉等工程的設(shè)計提供理論幫助[17]。

本文采用Morlet復(fù)小波變換，以榕江流域近60 a的水沙序列為基本資料，來探討榕江流域的水沙周期性變化特征。圖7為各序列距平值小波變換的實部等值線圖，圖中為實線且顏色較淺的部分為正，即為豐水期(多沙期)，虛線且顏色較深的部分為負(fù)，即為枯水期(少沙期)。等值線為0(實線與虛線相互過渡的區(qū)域)的地方表示豐枯交替(多沙-少沙交替)的轉(zhuǎn)折點。

圖7a可看出：南河徑流量在1953—2017年的演變過程中可分為3個時間尺度周期，分別為5 a以下、5～8 a與10～15 a。其中10～15 a在全局振蕩最強烈，并從1953—2017年經(jīng)歷了枯-豐交替的11個階段，且2017年剛好出現(xiàn)正值等值線，表明2017年后將進(jìn)入豐水期；5～8 a則是隨機出現(xiàn)。圖7b為北河徑流量在1967—2017年的小波變換實部等值線，可看出北河徑流量在1967—2017年的演變過程中尺度周期為5a以下、5～8 a與10～15 a，其中10～15 a從1967—2017年經(jīng)歷了枯-豐交替的9個階段，同樣2017年后迎來豐水期；80年代末，10～15 a的尺度周期振蕩減弱，5～8 a開始活躍。由圖7a、7b可得出：榕江流域徑流的演變過程中，10～15 a的尺度周期在整個時間段內(nèi)的發(fā)展最穩(wěn)定，為徑流的主周期。

圖7c、7d為泥沙序列小波變換的實部等值線，可見圖中左側(cè)等值線密度均大于右側(cè)，表明隨時間推移榕江流域的輸沙量與含沙量均減小，與前文整體趨勢性分析的結(jié)論一致，而進(jìn)入21世紀(jì)后輸沙量與含沙量的等值線較密，表明輸沙量與含沙量有所回升，這與前文泥沙序列逐年變化過曲線的狀態(tài)相符。圖中可看出輸沙量與含沙量的振蕩周期主要為5 a以下、5～8 a與10～15 a，主振蕩周期為10～15 a，從1956—2017年經(jīng)歷了“少沙-多沙”交替的11個階段，與徑流演變相對應(yīng)，而2017年等值線為負(fù)值，表明2017年后的一段時間仍為少沙期。由圖7c、7d可得出：輸沙量與含沙量的演變過程具有與徑流演變相似的特點，即10～15 a的尺度周期在整個時間段內(nèi)的發(fā)展最穩(wěn)定，為演變過程的主周期。但進(jìn)入21世紀(jì)后10～15 a的振蕩程度減弱，5～8 a的尺度周期逐漸活躍。

a)南河徑流量

為更加直觀表現(xiàn)水沙時間序列演變過程的主周期，以小波方差對小波變換進(jìn)行分析，小波方差的大小反應(yīng)各時段時間尺度周期振蕩的強弱，即方差越大，該時間對應(yīng)的尺度周期振蕩越強[18]。小波方差圖的峰值的對應(yīng)周期尺度即為序列的主周期。從圖8可看出，所有序列的方差峰值均出現(xiàn)在時間尺度周期為10～15 a內(nèi)，約為12 a，表明該尺度周期的振蕩最強，為序列的第一主周期，與上述小波實部等值線分析結(jié)論一致，圖中其他極值點也表明各序列還存在小周期。

a)南河徑流量

結(jié)合小波系數(shù)實部等值線與小波方差分析，取水沙序列的第一主周期12 a附近的小波系數(shù)進(jìn)行趨勢分析，見圖9，縱坐標(biāo)為第一主周期的小波系數(shù)，其為正值時表示徑流量、沙量在此時段偏多，對應(yīng)為豐水期(多沙期)，為負(fù)值時表示徑流量、沙量在此時段偏少，對應(yīng)為枯水期(少沙期)。圖中可直觀看出南河、北河徑流量分別經(jīng)歷了11次與9次豐枯交替，且2017年后小波系數(shù)出現(xiàn)正值，表明2017年后將進(jìn)入豐水期，這與上述徑流量小波系數(shù)等值線分析結(jié)論相互印證；輸沙量與含沙量均經(jīng)歷了11次交替，但在2017年的小波系數(shù)為負(fù)值，表明2017年后一段時間仍為少沙期，與上述泥沙序列的小波系數(shù)等值線分析結(jié)論相互印證。

a)南河徑流量

d)含沙量

2.5 降雨量和徑流量分析

降雨量是榕江流域地表徑流的唯一來源，地表徑流量的變化與降雨量基本一致(表2)[19]。榕江上游地勢陡峻，山高谷深，溪潤密布，河流落差大，多急流、瀑布，降水強度大，洪水匯流快；揭西縣河婆鎮(zhèn)以下，河谷逐漸開闊，比降較為平緩，地勢平坦；揭陽市三洲攔河閘以下為潮感區(qū)，河道更為平緩，兩岸農(nóng)田高程多在3 m以下，不同程度地受臺風(fēng)暴雨、風(fēng)暴潮影響，歷史上洪(潮)澇災(zāi)較為嚴(yán)重。

表2 榕江流域代表站多年平均降水量年內(nèi)分配統(tǒng)計

2.6 水文潮汐分析

流入海洋的河流稱為外流河，流入內(nèi)陸湖泊或消失于沙漠之中的這類瞎尾河稱為內(nèi)流河。榕江由汕頭市牛田洋入海，是南海熱帶氣旋常影響和登陸的地區(qū)之一，每年受其影響，少則1～2次，多則5～6次，熱帶氣旋(包括熱帶低壓、熱帶風(fēng)暴、強熱帶風(fēng)暴、臺風(fēng))帶來狂風(fēng)暴雨暴潮，這是榕江流域的一個重要氣候特征。

榕江多年平均水位3.32 m(珠江基面，下同)，歷史最高洪水位9.92 m(1970年9月14日)，最低水位2.29 m(1955年3月22日)。據(jù)東橋園水文站1954—1979年統(tǒng)計，榕江年平均流量87.4 m3/s，歷史上最大洪峰流量4 830 m3/s(1970年9月14日)，最枯流量為0(1954年3月31日斷流)。年平均含沙量0.224 kg/m3，年均輸沙量63.9萬t。北河與南河水土流失情況相似，但北河雖苦旱數(shù)月，也從沒斷流。

榕江受潮汐影響較大，漲潮時北河回水線(感潮區(qū))至瑯山上邊的北良。南河漲潮時回水線至三洲，楓江漲潮的回水線至潮州市浮崗。

3 結(jié)論分析

榕江上游由南、北兩河組成，于砲臺鎮(zhèn)雙溪咀匯合形成榕江干流，經(jīng)汕頭港內(nèi)的牛田洋注入南海。對榕江流域水沙過程的進(jìn)行趨勢性、突變性以及周期性的分析研究，以期對該流域的水資源管理與河道治理方面提供理論參考。研究主要得到以下幾點認(rèn)識和結(jié)論。

a)采用Mann-Kendall趨勢檢驗法對榕江流域的水沙資料進(jìn)行趨勢性分析，得出東橋園站與赤坎站兩站統(tǒng)計量均為正數(shù)，但絕對值小于顯著性水平α=0.05，表明南河與北河的多年徑流量均有增加的趨勢，但并不顯著，徑流量整體上圍繞一定值呈年際波動；輸沙量與含沙量的統(tǒng)計量均為負(fù)數(shù)，但絕對值大于顯著性水平，表明輸沙量與含沙量均有顯著減小的趨勢。

b)通過組合Mann-Kendall突變檢驗法、有序聚類分析法與累計距平法對有顯著變化的輸沙量與含沙量序列進(jìn)行突變性分析，經(jīng)檢驗分析，輸沙量與含沙量的突變點均為1987年，進(jìn)一步確定榕江流域的輸沙量在1987年發(fā)生不顯著突變，而含沙量在1987年發(fā)生顯著突變；以1987年為轉(zhuǎn)折點，榕江來沙過程整體上經(jīng)歷了多沙、少沙2個過程。

c)采用Morlet復(fù)小波變換分析榕江流域的水沙周期性變化特征，以小波系數(shù)方差圖與其相互印證，得到徑流量、輸沙量與含沙量的第一主周期約為12 a。以12 a為尺度周期，徑流量從1953—2017年經(jīng)歷了枯-豐交替的11個階段，并在2017年后將進(jìn)入豐水期；來沙過程從1955—2017年經(jīng)歷了少沙-多沙交替的11個階段，在2017年后一段時間仍為少沙期。

綜上幾種分析方法和數(shù)據(jù)可以看出，全球變暖，氣候變化，頻繁的人類活動等各類生產(chǎn)經(jīng)濟的建設(shè)，會直接對江河流域的趨勢性、突變性以及周期性帶來不可恢復(fù)的變化。在一定的空間尺度上，這種變化具有一定的重復(fù)性，也就是周期性。Mann-Kendall趨勢檢驗法，Morlet 小波分析法，在某周期內(nèi)能較為準(zhǔn)確分析計算出流域的趨勢和流域的周期突變，并開展流域水沙趨勢的預(yù)測。