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

敖亮挺

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

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

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

1 研究資料與方法

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

圖1 研究區域地理位置

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

1.1 滑動平均法

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

1.2 Mann-Kendall檢驗法

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

1.3 有序聚類分析法

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

1.4 小波分析法

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

2 研究結果

2.1 榕江流域水沙分布特性

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

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

根據東橋園水文站泥沙資料(1955—2017年)統計，榕江南河多年平均含沙量為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 徑流泥沙的趨勢性分析

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

a)南河徑流量

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

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

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

2.3 徑流泥沙的突變性分析

2.3.1泥沙序列突變分析

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

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

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

a)輸沙量

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

a)輸沙量

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

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

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

2.3.2沙量突變成因分析

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

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

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

c)航道疏浚的作用。榕江流域屬于平原地帶，防洪要求較高，航道整治通常采用挖槽進行[14]。20世紀90年代初，為滿足運移量不斷增長的需求，榕江流域進行了長距離的航道疏浚，榕江河道水流的挾沙能力隨著河道水深的增大和流速的減小而減小，泥沙受水流挾帶而運動，水流的挾沙能力與流速成正比。

2.4 徑流泥沙周期性分析

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

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

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

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

a)南河徑流量

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

a)南河徑流量

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

a)南河徑流量

d)含沙量

2.5 降雨量和徑流量分析

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

表2 榕江流域代表站多年平均降水量年內分配統計

2.6 水文潮汐分析

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

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

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

3 結論分析

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

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

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

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

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