珠江磨刀門河口日均水位變化及影響因子辨識(shí)

楊 昊，歐素英，傅林曦，劉 鋒，蔡華陽，楊清書

（1.中山大學(xué) 海洋工程與技術(shù)學(xué)院 河口海岸研究所 廣東 廣州 510275；2.河口水利技術(shù)國家地方聯(lián)合工程實(shí)驗(yàn)室 廣東 廣州 510275；3.廣東省海岸與島礁工程技術(shù)研究中心 廣東 廣州 510275；4.南方海洋科學(xué)與工程廣東省實(shí)驗(yàn)室（珠海），廣東 珠海 519000）

1 研究背景

河口日均水位作為全球氣候變化的重要指標(biāo)，其時(shí)空演變對河口三角洲系統(tǒng)的可持續(xù)發(fā)展具有重要指示意義。河口日均水位變動(dòng)的態(tài)勢由河道地形和動(dòng)力邊界（流量、海平面）等因子共同驅(qū)動(dòng)，全球氣候變暖導(dǎo)致海平面上升[1]將影響河口日均水位的時(shí)空分布；同時(shí)，隨著經(jīng)濟(jì)社會(huì)的高速發(fā)展，強(qiáng)人類活動(dòng)已成為河口三角洲系統(tǒng)演變的第三驅(qū)動(dòng)力[2]。當(dāng)氣候變化與強(qiáng)人類活動(dòng)干預(yù)的累積效應(yīng)達(dá)到甚至超過水位閾值（即自然恢復(fù)能力所對應(yīng)的水位變化極值）時(shí)，河口日均水位將會(huì)發(fā)生異變（本文所提“異變”是指強(qiáng)人類活動(dòng)驅(qū)動(dòng)，對河口地形及動(dòng)力邊界的改變，無論是變化幅度和頻度均超過自然因子所引起的變化，變化特征有別于自然過程，下同）。由于水位變化直接影響河口區(qū)的防洪和供水安全，因此，河口日均水位的異變態(tài)勢及其影響因子辨識(shí)，不僅是河口動(dòng)力學(xué)研究的重要科學(xué)問題，而且是河口治理和港口航道資源開發(fā)的關(guān)鍵[3-5]。

河口日均水位是徑潮動(dòng)力耦合的典型結(jié)果，其時(shí)空演變過程是探究河口復(fù)雜徑潮相互作用過程及機(jī)制的重要基礎(chǔ)。在徑流作用占主導(dǎo)河段，河口日均水位的變化主要取決于上游來水量[6-7]，而在感潮河段，隨著潮汐作用的增強(qiáng)，日均水位變化還受徑潮動(dòng)力非線性作用的影響。徑潮動(dòng)力耦合下單一感潮河段的水位梯度主要與非線性摩擦項(xiàng)相平衡[8-9]，因此，感潮河段沿河流方向形成正的水位梯度，河口日均水位由口門至上游沿程抬升，并具有顯著的大小潮和季節(jié)性變化[10-11]。采用功率譜分析、小波變換等分解方法能夠有效統(tǒng)計(jì)分析河口日均水位的時(shí)空分布[3，12]，基于數(shù)值模型與解析模型相結(jié)合的方法亦可對河口日均水位分布的形成機(jī)制進(jìn)行探討[13]。近年來，基于數(shù)據(jù)驅(qū)動(dòng)模型率定參數(shù)反演自然過程下水溫及水位的變化過程，從而辨識(shí)強(qiáng)人類活動(dòng)以及自然變化對河口水文要素的影響程度得到廣泛應(yīng)用[14-15]。

珠江三角洲（珠三角）是我國國家重大戰(zhàn)略——粵港澳大灣區(qū)建設(shè)的重要依托，其面積僅占全國國土面積的0.5%，創(chuàng)造出全國GDP的14.02%[16]。西江下游磨刀門作為珠三角入海徑流量和輸沙量最大的河口，是珠江河口最主要的泄洪和取水通道，探究磨刀門河口沿程日均水位的時(shí)空演變對保障珠三角的防洪和供水安全起著極其重要的作用。從1950年代開始，磨刀門河口先后經(jīng)歷了大規(guī)模的灘涂圍墾、聯(lián)圍筑閘和無序河道采沙等強(qiáng)人類活動(dòng)。在各種人類活動(dòng)的累積影響下，對磨刀門河口日均水位的影響進(jìn)行辨識(shí)研究具有典型性及必要性。本文聚焦于珠江磨刀門河口，基于數(shù)據(jù)驅(qū)動(dòng)模型對磨刀門河口長時(shí)間序列的日均水位資料進(jìn)行分析，辨識(shí)地形與動(dòng)力邊界變化對日均水位時(shí)空演變的影響，探究強(qiáng)人類活動(dòng)驅(qū)動(dòng)下日均水位的季節(jié)性異變特征。研究成果可為珠江河網(wǎng)區(qū)的防洪防澇、河口治理、航道整治、堤圍防護(hù)及水資源開發(fā)利用等提供科學(xué)依據(jù)。

2 研究區(qū)域概況

珠江河口具有三江（西江、北江和東江）匯流，八口（虎門、蕉門、洪奇門、橫門、磨刀門、雞啼門、虎跳門、崖門）入海的特點(diǎn)。根據(jù)徑流動(dòng)力與潮汐動(dòng)力的相對強(qiáng)弱，八大口門可劃分為河優(yōu)型和潮優(yōu)型河口。磨刀門河口的輸水輸沙量大、潮差小，屬于典型的徑流強(qiáng)、潮流弱的河優(yōu)型河口。作為西北江河網(wǎng)最主要的泄洪通道，磨刀門多年平均徑流量為923億m3，約占馬口站徑流量的三分之一，約占西、北江徑流量的四分之一，居珠江八大口門之首[17]。磨刀門潮汐屬于不正規(guī)半日潮，日不等現(xiàn)象顯著，口門內(nèi)潮流為往復(fù)流，流向與河道走向基本一致，口外主要潮流流向?yàn)槲鞅?東南向，漲潮流為西北向，落潮流為東南向。本文研究范圍為西江主干入?？凇サ堕T河口，上游流量控制站位為馬口水文站，沿程涵蓋甘竹、竹銀、燈籠山和三灶4個(gè)潮位站，具體站位如圖1（a）所示。據(jù)馬口水文站1959—2016年月均流量資料統(tǒng)計(jì)，其多年平均流量為7078 m3/s，最大月均流量達(dá)29 000 m3/s，最小月均流量為1210 m3/s，洪枯季差異顯著。

珠江磨刀門河口的典型強(qiáng)人類活動(dòng)主要包括：上游水庫建設(shè)、河道無序采沙、聯(lián)圍筑閘工程及口門圍墾整治，具體人類活動(dòng)作用區(qū)域見圖1（b）。珠江流域自1950年代開始已興建超過9000座水庫大壩，總庫容達(dá)到650億m3，約占珠江流域年徑流量的23 %[18]，其中西江流域水庫建設(shè)強(qiáng)度在1990—2000年間有較大的提高，期間庫容超過1億m3的水庫數(shù)達(dá)到10座，總庫容達(dá)到162億m3[19]，由于水庫工程擔(dān)任著“削洪調(diào)枯”的任務(wù)，會(huì)對下泄流量的季節(jié)性分布產(chǎn)生直接影響。磨刀門河口自1980年代起出現(xiàn)大規(guī)模河道（馬口-燈籠山河段）采沙活動(dòng)，至1990年代初達(dá)到高峰。高強(qiáng)度的采沙導(dǎo)致磨刀門馬口水文站的“水位-流量”關(guān)系發(fā)生顯著改變，多年平均流量條件下水位下降約達(dá)2 m[20-22]。聯(lián)圍筑閘旨在控支強(qiáng)干，聯(lián)圍并流?？陂T圍墾整治（燈籠山-三灶河段）導(dǎo)致河口快速向海延伸，至1990年代初磨刀門河口延伸長度達(dá)到13 km，河口水位逐漸抬升[23-24]。

總體而言，上游水庫調(diào)蓄導(dǎo)致馬口水文站流量（上游動(dòng)力邊界）呈現(xiàn)季節(jié)性異變，河道挖沙使得地形下切以及口門圍墾使得河道變窄（地形邊界）均導(dǎo)致磨刀門河口地形發(fā)生異變，而海平面的自然變化驅(qū)使下游動(dòng)力邊界發(fā)生改變?；诘匦闻c動(dòng)力邊界的時(shí)空變化不均勻性的特征，磨刀門河口日均水位具有明顯的階段性異變特征。

圖1 研究區(qū)域與強(qiáng)人類活動(dòng)概化示意圖

3 數(shù)據(jù)及方法

3.1 數(shù)據(jù)處理本文收集了1959—2016年馬口、甘竹、竹銀、燈籠山和三灶5個(gè)站的逐日高、低潮位，以及對應(yīng)時(shí)間段馬口站的日均流量數(shù)據(jù)，具體信息見表1（三灶站起始時(shí)間為1966年，其余站點(diǎn)起始時(shí)間均為1959年）。所用數(shù)據(jù)來源于《廣東省水文年鑒》和廣東省水文局。潮位原始數(shù)據(jù)高程基面為凍結(jié)基面，已統(tǒng)一轉(zhuǎn)換至珠江基面。將珠江磨刀門河口沿程站點(diǎn)的逐日高、低潮位數(shù)據(jù)插值得到逐時(shí)潮位，取24 h 平均統(tǒng)計(jì)得到各站點(diǎn)的日均水位。日均水位坡度S 定義為日均水位隨距離的變化率：

式中：Δz為兩站點(diǎn)之間的日均水位之差，m；Δx為兩站點(diǎn)之間的距離，m。

3.2 RTHA模型理論本文采用徑潮耦合的數(shù)據(jù)驅(qū)動(dòng)模型[25]（River-Tidal Harmonic Analysis，簡稱RTHA）用于反演徑流驅(qū)動(dòng)下的日均水位時(shí)空演變。該模型基于Matte 等[26]提出的非線性調(diào)和分析思路，在傳統(tǒng)調(diào)和分析理論的基礎(chǔ)上，假設(shè)河口任意位置的分潮振幅及相位變化主要受上游流量和地形變化的非線性調(diào)制影響，適用于河優(yōu)型河口，與本文研究區(qū)域的河口類型一致。

當(dāng)潮波沿河道上溯，潮水位的變化可由一個(gè)余水位項(xiàng)和多個(gè)分潮波（天文分潮、淺水分潮及氣象分潮等）的余弦曲線進(jìn)行線性疊加：

表1 磨刀門河口沿程主要潮位（水文）站點(diǎn)數(shù)據(jù)

假定河道平均水深為常數(shù)且不隨時(shí)間發(fā)生變化，引入潮波衰減系數(shù)公式：

式中：p0、p1、γ均為待定常數(shù)；Q為流量，m3/s，即河口潮汐的非線性衰減主要由上游徑流控制。

假設(shè)河口區(qū)潮波非線性衰減主要受徑流流量（Q）大小控制，結(jié)合式（2）和式（3），河口區(qū)水位z（tj）可由下式表示：

式中bi、ci為待定常數(shù)。

式（4）表明水位Z可由兩個(gè)隨時(shí)間變化的部分表示：

其中：

式中：d0、d1為狀態(tài)模型的待定參數(shù)；v0、v1、r0、r1為各分潮的待定參數(shù)；s(tj)為河床底高程、水深、徑流流量、海平面變化等引起的平均海平面（MSL）的變化，其隨徑流流量、外海潮汐的改變而改變，稱為平均水面狀態(tài)函數(shù)；F(tj)用來描述不同徑流條件對分潮的調(diào)制影響。

式（5）即為徑流驅(qū)動(dòng)下的RTHA模型。

本文根據(jù)實(shí)測水位數(shù)據(jù)，基于式（5）—式（7），建立多元回歸方程組，求解得出待定參數(shù)d0、d1及隱含潮汐調(diào)和常數(shù)信息的、各分潮待定參數(shù)v0、v1、r0、r1。但在最終解方程組之前，流量的冪次方（γ）是必須先給定的參數(shù)。本文采用粒子群優(yōu)化算法[27]（Particle Swarm Optimization），以水位的均方根誤差RMSE值最小或相關(guān)指數(shù)R2最大為目標(biāo)，迭代求解出流量的冪次方（γ）。均方根誤差RMSE和相關(guān)系數(shù)R2定義如下：

式中：m為各階段的時(shí)間序列長度；y為各階段水位實(shí)測值，m；yˉ為各階段水位實(shí)測值的均值，m。

本文雖采用分離徑潮信號(hào)的數(shù)據(jù)驅(qū)動(dòng)模型分離出徑流信號(hào)和徑流調(diào)制下的潮汐信號(hào)，但結(jié)果只選取由s(tj)與F(tj)之和所重構(gòu)出的日均水位z(tj)與實(shí)測日均水位進(jìn)行比較（均統(tǒng)一至日均水位尺度），并不分別討論s(tj)與F(tj)各自的變化，及各分潮“調(diào)和常數(shù)”在時(shí)間上的變化。

3.3 地形與動(dòng)力邊界對日均水位影響的辨識(shí)方法本文中馬口與三灶站日均水位作為邊界條件不參與模型計(jì)算，僅由馬口站日均流量驅(qū)動(dòng)對甘竹、竹銀、燈籠山站日均水位進(jìn)行模擬。累積距平法是顯示單一因素在其長時(shí)間序列中突變年份的常用手段[28]。本文采用歸一化的月均水位坡度（歸一化是指對原始時(shí)間序列進(jìn)行線性變換），即歸一化的月均水位坡度=（原始月均水位坡度-最小值）（/最大值-最小值）作為河口系統(tǒng)的特征動(dòng)力指標(biāo)，其構(gòu)成的時(shí)間序列為S，則任意時(shí)刻S的累積距平表示為：

式中：St為任意時(shí)刻t的月均水位坡度累積距平值；Si為第i月的月均水位坡度；Sˉ為n個(gè)月的月均水位坡度平均值。

通過繪制累積距平曲線將磨刀門河口日均水位序列分為緩變期、過渡期和調(diào)整期3個(gè)階段。將人類活動(dòng)作用前（地形未發(fā)生異變）自然狀態(tài)下的日均水位記為緩變期_實(shí)測值、強(qiáng)人類活動(dòng)改造地形邊界階段內(nèi)的日均水位記為過渡期_實(shí)測值、人類活動(dòng)作用減弱河口進(jìn)入恢復(fù)調(diào)整階段的日均水位記為調(diào)整期_實(shí)測值?；诰徸兤赺實(shí)測值率定出自然演變條件下日均水位與流量的函數(shù)關(guān)系與參數(shù)（緩變期_計(jì)算值），并以過渡期及調(diào)整期的日均流量驅(qū)動(dòng)預(yù)報(bào)出對應(yīng)階段的日均水位（即過渡期_預(yù)報(bào)值和調(diào)整期_預(yù)報(bào)值）。根據(jù)RTHA模型重構(gòu)出各站點(diǎn)日均水位的計(jì)算值（預(yù)報(bào)值），并結(jié)合各站點(diǎn)實(shí)測值，可通過下式計(jì)算分別得到流量邊界驅(qū)動(dòng)值ΔDIS、地形與海平面邊界驅(qū)動(dòng)值ΔGEO_MSL和總變化值ΔTOT：

ΔDIS等效于僅由上游流量邊界驅(qū)動(dòng)下的水位變化值，主要反演上游水庫建設(shè)對磨刀門沿程日均水位的影響（其中也包括流域的自然變化），而ΔGEO_MSL等效于剔除上游流量由地形改變（挖沙與口門圍墾）和海平面邊界影響下的日均水位變化值。采用緩變期所率定的參數(shù)來預(yù)測未受強(qiáng)人類活動(dòng)改造的自然水位演變，通過預(yù)報(bào)值與實(shí)測值的對比來達(dá)到辨識(shí)目的。結(jié)合式（11）—式（13）可知：ΔTOT=ΔGEO_MSL+ ΔDIS+ε，式中ε=z（緩變期_計(jì)算值）- z（緩變期_實(shí)測值），即采用緩變期日均水位的平均偏差來表示RTHA模型的誤差。

4 日均水位的階段性演變及影響辨識(shí)

4.1 日均水位的階段性演變月均水位坡度是研究河口區(qū)日均水位時(shí)空變化的切入點(diǎn)，本文采用三灶-馬口河段月均水位坡度的累積距平曲線對磨刀門河口日均水位的異變過程進(jìn)行階段性劃分。統(tǒng)計(jì)得到三灶-馬口段月均水位坡度并繪制累積距平曲線，結(jié)果如圖2所示。累積距平曲線呈現(xiàn)先增大、再波動(dòng)上升、后減小的階段性變化特征，結(jié)合強(qiáng)人類活動(dòng)作用于磨刀門的時(shí)間綜合判斷，可將1959—2016年分為3個(gè)階段，即1959—1990年為緩變期（該階段以自然變化為主，人類活動(dòng)影響的累積效應(yīng)甚微），1991—2000年為過渡期（該階段為強(qiáng)人類活動(dòng)累積效應(yīng)開始作用時(shí)期），2001—2016年為調(diào)整期（強(qiáng)人類活動(dòng)在該階段的強(qiáng)度減弱，磨刀門河口進(jìn)入自適應(yīng)調(diào)整階段），其中圖2中陰影部分為過渡期。從圖2中可知，在1990年之前曲線持續(xù)上升為正距平即水位坡度大于整體數(shù)據(jù)的平均值，經(jīng)歷過渡期后曲線下降，表明2000年后為負(fù)距平即水位坡度小于整體數(shù)據(jù)平均值。累積距平曲線呈現(xiàn)的結(jié)果表明自然調(diào)節(jié)與強(qiáng)人類活動(dòng)驅(qū)動(dòng)下磨刀門河口水位坡度變緩。

圖2 三灶-馬口段歸一化月均水位坡度累積距平線

根據(jù)累積距平曲線對時(shí)間節(jié)點(diǎn)的劃分，做出甘竹、竹銀及燈籠山站緩變期（1959—1990年）和調(diào)整期（2001—2016年）的多年日平均水位過程線如圖3所示。日均水位的季節(jié)分布表現(xiàn)為洪季大于枯季即上游流量越大日均水位抬升越明顯，緩變期洪峰位于7月而進(jìn)入調(diào)整期后洪峰轉(zhuǎn)移至6月。當(dāng)磨刀門河口進(jìn)入調(diào)整期，甘竹、竹銀站日均水位明顯下降，由緩變期至調(diào)整期，甘竹站日均水位平均下降0.49 m，竹銀站日均水位平均下降0.14 m。甘竹站調(diào)整期日均水位的季節(jié)分布均低于緩變期，竹銀站2月份水位略高于緩變期，而其余時(shí)間段分布規(guī)律與甘竹站一致。燈籠山站調(diào)整期與緩變期季節(jié)分布差異較為復(fù)雜，調(diào)整期內(nèi)夏、秋兩季日均水位低于緩變期，而春、冬兩季日均水位高于調(diào)整期，由緩變期至調(diào)整期，燈籠山站日均水位平均上升0.01 m。

圖3 磨刀門河口站點(diǎn)不同階段日均水位的季節(jié)分布

4.2 RTHA模型的率定與驗(yàn)證基于磨刀門河口實(shí)測日均水位的階段性演變特征，采用RTHA模型反演緩變期階段的日均水位時(shí)空變化，為模型效果評價(jià)提供依據(jù)。將緩變期分成2個(gè)階段，緩變期前2/3（1959—1980年）時(shí)段用于參數(shù)率定并預(yù)報(bào)緩變期后1/3（1981—1990年）時(shí)段作進(jìn)一步驗(yàn)證。模型率定及驗(yàn)證通過式（8）均方根誤差（RMSE）和式（9）相關(guān)系數(shù)（R2）來顯示其模型效果。

圖4為模型率定驗(yàn)證結(jié)果。表2顯示的是緩變期以及率定與驗(yàn)證階段的RMSE、R2及利用PSO算法得到的參數(shù)γ，而各分潮參數(shù)由于數(shù)目眾多并未在表2中顯示。模型相關(guān)系數(shù)大于0.78，表明模擬效果較好，但由圖4（a）可見該模型在大流量條件下存在較大誤差；甘竹站緩變期整體日均水位的RMSE為0.19 m，竹銀站緩變期整體日均水位的RMSE為0.19 m，燈籠山站整體緩變期日均水位的RMSE為0.20 m。RTHA數(shù)據(jù)驅(qū)動(dòng)模型反演日均水位效果良好，可用于后續(xù)預(yù)報(bào)分析，后續(xù)預(yù)報(bào)過渡與調(diào)整期日均水位統(tǒng)一采用緩變期（1959—1990年）率定的參數(shù)γ（具體參數(shù)見表2）。

圖4 磨刀門河口各站點(diǎn)RTHA模型計(jì)算值與實(shí)測值的對比

表2 率定及驗(yàn)證階段RTHA模型參數(shù)及效果

4.3 日均水位的影響因子辨識(shí)圖5是磨刀門河口甘竹、竹銀與燈籠山站過渡及調(diào)整期日均ΔGEO_MSL與ΔDIS的季節(jié)分布。根據(jù)甘竹、竹銀與燈籠山站在磨刀門河口的空間位置其動(dòng)力條件可分為3種類型，甘竹站為徑流優(yōu)勢型、竹銀站為徑潮相互作用型、燈籠山站為潮流優(yōu)勢型。從ΔGEO_MSL季節(jié)分布來看，地形與海平面共同邊界驅(qū)動(dòng)下磨刀門河口3個(gè)站點(diǎn)由過渡期至調(diào)整期存在明顯的累積現(xiàn)象，甘竹、竹銀站累積效應(yīng)主要體現(xiàn)在夏季，而燈籠山站的累積效應(yīng)主要體現(xiàn)在秋、冬兩季，說明人類活動(dòng)改變地形以及海平面變化的影響都具有累積效應(yīng)且呈現(xiàn)明顯的季節(jié)特征。從過渡期至調(diào)整期，春夏兩季的ΔGEO_MSL主要由挖沙導(dǎo)致的河床下切引起，累積效應(yīng)量值持續(xù)增大方向?yàn)樨?fù)（即日均水位下降），秋冬兩季由海平面及口門圍墾主導(dǎo)，ΔGEO_MSL持續(xù)增大方向?yàn)檎慈站惶蠲黠@的月份為10—11月。

三灶站日均水位的季節(jié)分布如圖6（b）所示，其近似等效于海平面的季節(jié)分布，對比發(fā)現(xiàn)海平面在秋季（9—11月）有明顯上升，與竹銀、燈籠山站ΔGEO_MSL大于0的時(shí)間段吻合。海平面以及口門圍墾對日均水位的影響越往上游其影響程度越弱，其影響主要體現(xiàn)在潮流優(yōu)勢及徑潮相互作用河段（甘竹站調(diào)整期的ΔGEO_MSL僅在11月份略微大于0），而挖沙對日均水位的影響主要體現(xiàn)在徑流優(yōu)勢與徑潮相互作用河段。

僅由上游流量邊界驅(qū)動(dòng)引起的日均水位變化ΔDIS在過渡期呈現(xiàn)洪季抬升、枯季下降的變化趨勢，這與過渡期內(nèi)上游洪季流量大于緩變期，而枯季流量小于緩變期相對應(yīng)，并未體現(xiàn)出水庫的調(diào)蓄作用。西江磨刀門上游水庫建設(shè)始于1964年，但庫容小，大庫容水庫投入使用的時(shí)間為2000年（天生橋一級(jí)，庫容量為103.60億m3）其庫容量超過之前上游水庫建設(shè)的總庫容[18]，因此水庫調(diào)蓄作用主要體現(xiàn)在調(diào)整期。當(dāng)磨刀門河口進(jìn)入調(diào)整期后，ΔDIS在7—11月出現(xiàn)明顯下降，而1、2、12月出現(xiàn)明顯抬升（圖5（d）—圖5（f）），這與圖6（a）中馬口站流量逐月分布因水庫調(diào)蓄而發(fā)生洪季減小、枯季增大的季節(jié)性異變相對應(yīng)。

圖5 過渡及調(diào)整階段ΔGEO_MSL和ΔDIS的季節(jié)分布

圖6 磨刀門河口馬口站流量與三灶站水位的季節(jié)分布

地形與動(dòng)力邊界驅(qū)動(dòng)下磨刀門河口日均水位調(diào)整期與緩變期的年均差值與季節(jié)差值如表3所示。由表3的年均差值可知，當(dāng)磨刀門河口進(jìn)入調(diào)整期后，甘竹、竹銀的日均水位均有不同程度下降，ΔTOT反映出越往上游（徑流優(yōu)勢河段）下降幅度越大（甘竹站ΔTOT= -0.47 m），而下游燈籠山站ΔTOT為正值即日均水位抬升（燈籠山站ΔTOT= 0.01 m）。竹銀與甘竹站ΔGEO_MSL為負(fù)值，其中甘竹站ΔGEO_MSL達(dá)到-0.37 m，竹銀站為-0.10 m，而燈籠山站是挖沙主要區(qū)域，但ΔGEO_MSL為正值（ΔGEO_MSL= 0.02 m），說明燈籠山站受到與挖沙影響趨勢相反的人類活動(dòng)（如口門圍墾）及海平面的影響，且對日均水位的影響要大于挖沙活動(dòng)。各站點(diǎn)ΔDIS均為負(fù)值，說明受上游水庫調(diào)蓄影響下馬口站日均流量整體減小導(dǎo)致磨刀門河口日均水位普遍下降。從量值上看ΔGEO_MSL對日均水位的影響明顯大于ΔDIS。

對比各站的季節(jié)差值，甘竹與竹銀站四季ΔTOT變化趨勢一致，均表明日均水位下降，呈現(xiàn)洪季大于枯季，夏季最大，冬季最小，春秋兩季交替的現(xiàn)象。燈籠山站四季ΔTOT呈現(xiàn)秋冬季抬升（秋季ΔTOT= 0.01 m、冬季ΔTOT= 0.06 m）而春夏季減小的變化趨勢，與海平面的季節(jié)性變化相對應(yīng)（圖6（b））。各站ΔDIS的四季變化規(guī)律一致，均呈現(xiàn)冬季抬升、春夏秋三季下降，這與水庫洪季蓄水、枯季放水有良好的對應(yīng)關(guān)系（圖6（a））。但ΔGEO_MSL的季節(jié)性差異出現(xiàn)不同的變化趨勢，甘竹站四季均值為負(fù)說明該區(qū)域水位變化主要受挖沙導(dǎo)致的地形下切所控制；竹銀站春夏冬三季變化規(guī)律與甘竹站一致，但在秋季ΔGEO_MSL為正值（ΔGEO_MSL= 0.03 m），表明海平面以及口門圍墾在竹銀站秋季日均水位中占主導(dǎo)作用；燈籠山站的ΔGEO_MSL與其ΔTOT規(guī)律一致，由于燈籠山站受到兩種人類活動(dòng)的影響且影響的趨勢剛好相反，即在口門圍墾和海平面綜合作用下使得日均水位抬升，結(jié)果表現(xiàn)為燈籠山站秋冬兩季的ΔGEO_MSL為正值（秋季ΔGEO_MSL= 0.16 m、冬季ΔGEO_MSL= 0.02 m），而挖沙使得日均水位下降導(dǎo)致燈籠山站春夏兩季ΔGEO_MSL為負(fù)值（春季ΔGEO_MSL= -0.05 m、夏季ΔGEO_MSL= -0.02 m）。

表3 磨刀門沿程站點(diǎn)緩變期至調(diào)整期的日均水位差（單位：m）

5 討論

磨刀門河口主要潮位站日均水位發(fā)生階段性變化，勢必影響沿程日均水位坡度的時(shí)空分布。圖7為燈籠山-甘竹河段日均水位坡度的季節(jié)分布。

將式（11）和式（12）中日均水位替換成日均水位坡度，則藍(lán)、黑線圍成的面積代表不考慮地形因素，僅由上游流量驅(qū)動(dòng)的日均水位坡度變化值即ΔS-DIS，紅、黑線圍成的面積代表地形及海平面邊界共同引起的日均水位坡度變化值即ΔS-GEO_MSL（圖7陰影面積）。整體而言，燈籠山-甘竹河段日均水位坡度中ΔS-GEO_MSL所圍的面積明顯大于ΔS-DIS所圍面積，表明磨刀門河口地形與海平面邊界對日均水位坡度的影響大于上游流量邊界的影響。磨刀門河口燈籠山-甘竹河段日均水位坡度各參數(shù)如表4所示。由統(tǒng)計(jì)的日均水位坡度ΔS-TOT的變化可知，整體而言強(qiáng)人類活動(dòng)以及海平面上升使得磨刀門燈籠山-甘竹河段日均水位坡度變緩（即為負(fù)值），夏季降幅最大（量值為3.43×10-6），其次是春季和秋季，冬季變緩幅度最小。日均水位坡度的ΔS-GEO_MSL四季變化趨勢與ΔS-TOT一致，只是量值上有所差別；而日均水位坡度的ΔS-DIS四季變化規(guī)律有所不同，受水庫調(diào)蓄影響，秋季降幅大于春季（與ΔS-TOT、ΔS-GEO_MSL規(guī)律不一致），而在冬季日均水位坡度的ΔS-DIS為正值（量值為1.13×10-7），即水庫在冬季放水使得燈籠山-甘竹河段日均水位坡度變陡。

表4 燈籠山-甘竹河段緩變期至調(diào)整期的日均水位坡度差

圖7 甘竹-燈籠山段日均水位坡度的季節(jié)分布

綜上所述，磨刀門主要的三類強(qiáng)人類活動(dòng)（上游水庫建設(shè)疊加上無序挖沙和口門圍墾）均導(dǎo)致沿程日均水位坡度變緩。根據(jù)簡化的一維動(dòng)量守恒過程可知，平均水位坡度項(xiàng)主要與有效摩擦項(xiàng)相平衡[11]，當(dāng)平均水面坡降減緩，外海潮波傳播所受的有效摩擦減弱，因此磨刀門河口潮差增大，咸潮上溯加劇。這表明珠江磨刀門河口進(jìn)入調(diào)整期后的防洪納潮、潮波傳播及沖淤演變等正面臨著水位坡降變緩和潮汐動(dòng)力增強(qiáng)的影響。

6 結(jié)論

本文基于珠江磨刀門河口沿程5個(gè)站位1959—2016年的日均水位數(shù)據(jù)以及馬口站日均流量數(shù)據(jù)，探究地形與動(dòng)力邊界異變格局下，甘竹、竹銀、燈籠山站日均水位以及燈籠山-甘竹河段日均水位坡度的季節(jié)性異變，量化僅由上游流量驅(qū)動(dòng)、地形變化與海平面共同驅(qū)動(dòng)的日均水位（坡度）變化值，得到以下主要結(jié)論：

（1）由三灶-馬口河段月均水位坡度的累積距平曲線可知，磨刀門河口的日均水位呈現(xiàn)出階段性變化，1959—1990年為自然變化主導(dǎo)的緩變期、1991—2000年為強(qiáng)人類活動(dòng)主導(dǎo)的過渡期、2001—2016年為自適應(yīng)調(diào)整期。甘竹與竹銀站日均水位由緩變期至調(diào)整期逐漸下降（分別下降0.49 m和0.13 m），燈籠山站日均水位總體上升0.01 m。因此造成磨刀門河口沿程水位坡降變緩，由緩變期至調(diào)整期，燈籠山-甘竹河段日均水位坡度平均下降6.75×10-6。

（2）地形和海平面邊界共同驅(qū)動(dòng)下（ΔGEO_MSL），磨刀門河口徑流優(yōu)勢河段以及徑潮相互作用河段的日均水位變化，均以河道下切作用主導(dǎo)，海平面和口門圍墾作用僅在徑潮相互作用河段的秋季成為主導(dǎo)因子。而在潮流優(yōu)勢河段中，洪季（春、夏）以河道下切因子為主，枯季（秋、冬）則替換成海平面以及口門圍墾因子主導(dǎo)；其次ΔGEO_MSL具有累積效應(yīng)，徑流優(yōu)勢與徑潮相互作用河段的累積效應(yīng)均體現(xiàn)在夏季，而潮流優(yōu)勢河段則體現(xiàn)在秋、冬兩季。上游流量邊界驅(qū)動(dòng)下（ΔDIS），磨刀門河口春夏秋三季的日均水位變化以水庫蓄水主導(dǎo)，冬季則以水庫調(diào)枯作用為主。整體而言，ΔGEO_MSL對磨刀門河口日均水位的影響明顯大于ΔDIS。

（3）磨刀門河口上游水庫建設(shè)導(dǎo)致僅由流量邊界驅(qū)動(dòng)下，燈籠山-甘竹河段冬季日均水位坡度變陡（冬季ΔS-DIS為1.13×10-7）；強(qiáng)人類活動(dòng)改造地形邊界疊加海平面變化導(dǎo)致燈籠山-甘竹河段日均水位坡度變緩，夏季ΔS-GEO_MSL變緩幅度最大（夏季ΔS-GEO_MSL為-3.43×10-6），其次為春季、秋季，冬季變緩幅度最小（冬季ΔS-GEO_MSL為-6.95×10-7）。磨刀門河口日均水位坡度變緩導(dǎo)致潮波沿河口傳播所受的有效摩擦減弱，磨刀門河口將面臨潮汐動(dòng)力增強(qiáng)導(dǎo)致咸潮與風(fēng)暴潮強(qiáng)度加劇等自然災(zāi)害的影響。