珠江磨刀門河口水位與海平面、上游流量 的聯(lián)合分布關(guān)系異變研究*

唐啟邦 歐素英 蔡華陽 楊清書

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

河口是陸海相互作用的過渡區(qū)域， 其水位變化一方面是徑潮動力相互作用的直接結(jié)果， 既受自下而上的河口潮波、海平面上升等海洋動力的影響， 也受自上而下的上游徑流作用(Caiet al， 2016， 2019a)； 另一方面強人類活動引起的地形邊界改變， 直接或間接導(dǎo)致河口三角洲水位、潮差、鹽度等動力結(jié)構(gòu)的變化， 引起河口沉積作用過程、化學及生態(tài)過程的調(diào)整改變， 進而對河口水資源規(guī)劃、防洪、排澇等河口治理策略提出了新的挑戰(zhàn)。作為我國人口集聚程度極大的珠江河口地區(qū)， 隨著其經(jīng)濟的快速發(fā)展及城市規(guī)模的迅速擴張， 人類活動對河口環(huán)境的影響已遠超環(huán)境自我修復(fù)能力， 導(dǎo)致河口環(huán)境急劇惡化。研究表明， 近幾十年來珠江三角洲河道大規(guī)模無序挖砂、口門圍墾、航道整治及上游水庫建設(shè)等強人類活動導(dǎo)致河口三角洲的地形邊界、水文環(huán)境已經(jīng)發(fā)生明顯異變， 如口門水道縮窄延長、河道下切、過水斷面面積增大， 珠江來沙量顯著減少， 三角洲內(nèi)馬口分流分沙比減少、三水分流分沙比增加等( 燊陳子 等， 2015； 倪培桐等， 2016； Wuet al， 2016； 陸永軍等， 2019； Liuet al， 2019)， 相同珠江流量條件下三角洲中上段的水位明顯下降， 水位出現(xiàn)由緩慢上升轉(zhuǎn)為急劇下降的異變等(蔣陳娟等， 2012)； 而水位異變導(dǎo)致洪季河網(wǎng)區(qū)水位坡降變緩， 泄洪排澇效率降低， 枯季潮汐動力明顯增強， 河口鹽水入侵加劇(黃洪城等， 2014； Caiet al， 2019b)。珠江三角洲水文的自然演變規(guī)律在強人類活動干預(yù)下甚至出現(xiàn)自然過程難以達到的異變， 而河口水文過程及其規(guī)律是水資源管理利用的重要依據(jù)， 因此研究陸海相互作用下強人類活動對河口三角洲區(qū)水位變化的影響， 既是河口動力過程研究的熱點問題， 也對探討人類活動驅(qū)動下河口區(qū)域洪潮水位預(yù)報、水資源的高效利用以及河口治理等具有重要的現(xiàn)實意義， 為今后磨刀門河口的綜合整治規(guī)劃、有序開發(fā)航運資源、防洪排澇和調(diào)水壓咸等迫切社會經(jīng)濟需求提供重要技術(shù)支撐。

近20 年來， 眾多研究者采用小波變換、理論解等不同的方法辨識河口地區(qū)不同時間尺度的水位變化及機制， 比如潮周期變化、水位的長期變化趨勢、極端低水位及設(shè)計水位的變化等(歐素英等， 2004a； 歐素英等， 2004b； 蔣陳娟等， 2012； 燊陳子 等， 2013， 蔡華陽等， 2018； 楊昊等， 2019； Yanget al， 2019a)。但如何從多因素共同影響下的水位變化中辨識出強人類活動的作用仍有待于進一步研究。近年來廣泛應(yīng)用于水文要素研究中的Copula 聯(lián)合概率分布方法， 在研究多邊界聯(lián)合約束下的水文要素分布上具有明顯的 優(yōu) 勢(陳 子 燊 等， 2015， 2018； 陳 玲 舫 等， 2017； Bacopoulos， 2017； Liuet al， 2018)。本文嘗試從概率統(tǒng)計分布角度， 采用Copula 聯(lián)合概率分布方法， 選取受強人類活動影響顯著的珠江磨刀門河口為研究區(qū)域， 以珠江流量為上游邊界， 下邊界不考慮潮汐的短周期變化， 以河口月均海平面為月均態(tài)下海洋動力的綜合驅(qū)動因子， 擬研究磨刀門月均水位與上下邊界條件即三灶海平面(即平均海平面水位)和馬口流量之間的聯(lián)合依賴關(guān)系， 比較研究上游流量和外海海平面共同約束下強人類活動前后的水位概率分布特征及極端水位風險的變化過程， 對比研究其聯(lián)合依賴關(guān)系異變， 進而定量計算強人類活動對水位變化的影響程度。

1 研究區(qū)域概況

本文研究區(qū)域為珠江磨刀門河口， 包括西江干流水道-磨刀門水道(如圖1 所示)。西北江水沙經(jīng)思賢滘口分水分沙后， 約78%—84%的徑流量通過西江干流水道進入西北江三角洲。磨刀門作為珠江八大口門中入海徑流量和輸沙量最大的河口， 約36%的珠江徑流量通過磨刀門進入南海， 是珠江最重要的泄洪通道。磨刀門河口潮差較小， 口門三灶站多年平均潮差小于1m。因此， 磨刀門河口徑流動力強而潮汐動力弱， 是典型的河優(yōu)型河口。基于上游馬口水文控制站1965—2016 年的月均流量統(tǒng)計資料， 可得洪季(4—9 月)歷史最大流量為 29000m3/s， 平均流量為10810m3/s； 枯季(10 月—次年3 月)歷史最大流量為14500m3/s， 平均流量為3428m3/s； 年均歷史最大流量為9978m3/s， 平均流量為7119m3/s。口門三灶站1965—2016 年的月均海平面數(shù)據(jù)顯示， 其洪季歷史最大海平面為0.90m， 平均海平面為0.63m； 枯季歷史最大海平面為0.87m， 平均海平面為0.56m； 年平均歷史最大海平面為0.70m， 平均海平面為0.59m。

圖1 研究區(qū)域圖 Fig.1 Map of the study area

從20 世紀50 年代開始， 磨刀門河口先后經(jīng)歷了聯(lián)圍筑閘、灘涂圍墾、航道整治規(guī)劃以及大規(guī)模的河道無序采砂等強人類活動。特別是20 世紀90 年代磨刀門完成口門整治后， 口門人為向海快速延伸約16km， 使得原本具有泄洪納潮功能的寬闊內(nèi)海灣消失， 河口環(huán)境轉(zhuǎn)變?yōu)槭苋斯?dǎo)堤控制的河道， 口門迅速向海延伸， 延伸速度超越自然過程變化， 口門邊界條件發(fā)生明顯異變。西江上游流域方面， 建有大藤峽等骨干水源工程， 以西江龍灘及大藤峽等水庫為骨干構(gòu)建水資源調(diào)配體系， 直接改變磨刀門河口的來水來沙條件(張旭等， 2013)。而對磨刀門河口地形影響最為明顯的是大規(guī)模、無序的河道采砂。根據(jù)調(diào)查分析(韓龍喜等， 2005)， 1991—1994 年(采砂高峰期)3年內(nèi)從珠海大橋至竹銀近口段僅24.4km 的范圍內(nèi)， 共挖河沙2.4×107m3， 平均年挖砂量8×106m3； 1994—2000 年6 年內(nèi)共挖砂 1.92×107m3， 平均年挖砂量3.2×106m3； 1997 年因在洪灣水道中段建“洪灣北保稅區(qū)”， 需挖砂吹填， 從橫州水道北口左側(cè)河床采砂， 采砂量約為1×107m3。大規(guī)模挖砂以及航道整治引起磨刀門河口上游河段河床大幅度下切， 河床深泓平均下切0.59—2.25m， 河床沖淤演變趨勢由原來淤積為主轉(zhuǎn)變?yōu)轱@著沖刷(賈良文等， 2002)。自禁沙令頒布后， 2000 年后河道采砂才逐漸減少。

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

2.1 數(shù)據(jù)及處理

本文以1965—2016 年的馬口月均流量作為上游河流邊界， 不考慮外海潮汐的短周期變化影響， 以三灶站的月均海平面為外海動力邊界， 分別采用燈籠山、竹銀、江門、甘竹站的月均水位序列進行聯(lián)合依賴關(guān)系研究， 站位位置如圖1 所示。根據(jù)人類活動歷程(羅憲林等， 2002； 劉鋒等， 2011； 袁菲等， 2018)， 自1960 年以來， 磨刀門河道整體上橫斷面平均河寬縮窄， 水深增加， 寬深比減小， 主槽斷面面積增加； 20 世紀80 年代后期磨刀門開始大規(guī)模圍墾整治， 直至20 世紀90 年代中期才基本結(jié)束。西北江水沙經(jīng)思賢滘 馬口三水多年分流、分沙比大幅度變化發(fā)生在1993 年， 其中洪季同多年變化一樣， 也發(fā)生在1993 年， 枯季比多年變化提早4 年， 發(fā)生在1989 年(胡德禮等， 2009)。基于此， 本文將1965—1985 年劃分為強人類活動前自然演變階段， 2000—2016 年作為強人類活動劇烈影響后階段， 1985—2000 年為相應(yīng)的過渡階段， 比較研究在馬口流量和三灶海平面影響下磨刀門河口沿程各站水位對強人類活動的響應(yīng)過程及演變規(guī)律。

2.2 Copula 方法

自然系統(tǒng)中水文事件大部分都是多因素共同影響作用的結(jié)果。Sklar(1959)提出Copula 理論， 是研究多變量相關(guān)性及關(guān)聯(lián)性的重要數(shù)學工具， 能夠有效建立多變量之間的非線性關(guān)系， 且其邊緣分布靈活， 形式不受邊緣分布限制。目前Copula 方法已被廣泛應(yīng)用于洪水、暴雨、干旱等多變量水文事件的聯(lián)合概率分布研究( 燊陳子 等， 2013； Yanget al， 2019b)。

根據(jù) Sklar 定理， 若F是一個p維隨機變量X=(X1，???，XP)的累積分布函數(shù)， 其邊緣分布是連續(xù)函數(shù)Fi(i=1，???，p)， 則存在唯一一個Copula函數(shù)C， 使得：

一般來說， 如果x1出現(xiàn)超過某個極大值或低于某個極小值的極端情況時， 則事件x1可定義為危險或處于危險之中(Salvadoriet al， 2004)。例如， 當河口沿程各站水位低于該站10%水位累積概率所對應(yīng)的值時定義為低水位事件， 記為F(x1)low_WL≤10%， 類似的， 高水位事件定義為水位高于第90 分位數(shù)所對應(yīng)的水位值， 記為F(x1)high_WL≥90%。高水位及低水位事件的條件概率可由式(2)或(3)計算。

由于河口三角洲的月均水位既受到上游流量的影響， 同時也受到外海海平面變化的作用， 本文通過Copula 函數(shù)來構(gòu)建水位-海平面-流量的聯(lián)合依賴關(guān)系模型， 不僅可以描述水位與海平面、流量的非線性依賴關(guān)系， 而且能夠通過比較強人類活動前后聯(lián)合依賴關(guān)系的變化， 從概率統(tǒng)計角度定量研究強人類活動對磨刀門河口沿程各站水位的影響程度。具體步驟為： 1) 首先通過相關(guān)分析初步研究強人類活動影響下磨刀門河道沿程各站水位與海平面、流量的線性關(guān)系是否發(fā)生明顯變化。2) 確定強人類活動前后磨刀門河口沿程各站月均水位、馬口流量及三灶海平面的最優(yōu)概率分布， 即選擇四種理論概率分布(即高斯分布、伽馬分布、對數(shù)正態(tài)分布和威布爾分布)對各變量數(shù)據(jù)進行擬合， 采用卡方檢驗進行擬合優(yōu)度比選， 確定最優(yōu)邊緣分布， 并初步研究人類活動前后單變量概率分布的變化。3) 確定強人類活動前后磨刀門河口沿程各站月平均水位與上游流量、月平均海平面的最優(yōu)Copula 函數(shù)， 分析變量間相互依賴關(guān)系的變化?；诔喑匦畔⒘繙蕜t(Akaike information criterion， AIC)， 從四個Copula 函數(shù)即Gaussian 函數(shù)和阿基米德Copula 家族中Clayton、Gumbel 和Frank 函數(shù)中選擇最優(yōu)Copula 函數(shù)來連接各變量， 構(gòu)建雙變量和三變量的聯(lián)合分布模型， 比較分析人類活動前后兩個階段水位與流量、海平面之間聯(lián)合依賴結(jié)構(gòu)的變化。4) 根據(jù)構(gòu)建的水位與海平面、流量的聯(lián)合依賴模型， 通過式(2)或(3)計算任意給定流量或(和)海平面條件下出現(xiàn)某個水位事件(特別是高水位及低水位事件)的條件概率。通過比較強人類活動前后相同條件概率所對應(yīng)的水位事件的改變， 或同樣水位事件所對應(yīng)條件概率的變化， 定量分析強人類活動對水位的影響程度。

3 結(jié)果分析

3.1 人類活動前后月均水位與海平面、流量線性關(guān)系的變化

強人類活動前(1965—1985 年)和強人類活動后(2000—2016 年)各站水位與三灶海平面及馬口流量的相關(guān)系數(shù)如表1 所示。由表1 可見， 自甘竹至燈籠山段， 其月均水位與海平面的線性相關(guān)系數(shù)隨距三灶站距離的縮小逐漸增大， 接近外海三灶的燈籠山站， 其相關(guān)系數(shù)最大。強人類活動后， 外海海平面對沿程各站水位的線性影響有所增強， 水位-海平面的線性相關(guān)系數(shù)較人類活動前均略有增加， 其中， 竹銀站的增加幅度最大(由0.76 增大到0.82)。反之， 上游馬口流量對磨刀門沿程各站水位的影響在強人類活動影響后則有所減弱。其中， 甘竹和江門站水位-流量的線性相關(guān)系數(shù)分別由0.98、0.98 分別減小至0.97、0.95， 上游流量影響強度略有減小， 但強人類活動后磨刀門水道的竹銀和燈籠山站的水位-流量相關(guān)系數(shù)分別由0.94、0.74 大幅減少為0.81、0.56， 即強人類活動影響下水位與馬口流量的線性關(guān)系顯著減弱。

表1 強人類活動前后水位-海平面、水位-流量的線性相關(guān)系數(shù) Tab.1 Coefficients of water level - sea level and water level - river discharge before and after intensive human interventions

3.2 月均水位與海平面-流量的聯(lián)合依賴關(guān)系模型的構(gòu)建

圖2 各站不同變量邊緣分布的卡方統(tǒng)計量 Fig.2 Chi-square statistics of marginal distributions for different variables at different stations

單變量邊緣分布的卡方檢驗結(jié)果如圖2 所示。各站月均水位的最優(yōu)概率分布有所差異， 其中， 竹銀站更傾向于伽馬分布， 燈籠山站服從廣義極值分布， 而 馬口流量更符合威布爾分布。采用圖2 中各變量的最優(yōu)概率分布， 可得強人類活動前后的累積頻率分布曲線， 如圖3 所示。受強人類活動影響后， 沿程各站低累積概率(＜0.18)所對應(yīng)的水位有所增大(曲線右移)， 反之， 高累積概率(＞0.18)所對應(yīng)的水位有所減小(曲線左移)。對于海平面變化來說， 強人類活動影響后的概率密度分布曲線整體右移， 表明高水位事件概率明顯增大， 海平面變化有較明顯的上升趨勢。對馬口流量而言， 強人類活動前后的概率分布曲線形態(tài)比較一致， 但因思賢 滘分流比的減小和上游流量峰值的降低， 平均流量所對應(yīng)的概率有所增大， 而相同累積概率條件下所對應(yīng)的流量有所降低(累積概率曲線整體左移)。由圖3 可見擬合的邊緣分布與觀測數(shù)據(jù)吻合較好， 能夠用于后續(xù)依賴關(guān)系的異變研究。

圖3 強人類活動前后磨刀門河口沿程各站變量邊緣分布累積頻率曲線 Fig.3 Cumulative frequency curves of marginal distribution of different variables at stations along the Modaomen estuary before and after intensive human interventions

由2.2 節(jié)可知， 根據(jù)AIC 準則選擇最佳的Copula函數(shù)， 結(jié)果如表2 所示。選取AIC 值最小的Gaussian函數(shù)作為水位-海平面和水位-流量的最佳擬合Copula聯(lián)合分布函數(shù)用于后續(xù)依賴關(guān)系的異變研究。

3.3 強人類活動前后水位-海平面和水位-流量的聯(lián)合依賴關(guān)系變化

圖4 和圖5 分別給出各站強人類活動前后水位-海平面和水位-流量的聯(lián)合累積頻率關(guān)系曲線， 從中可以更直觀的看出強人類活動影響下水位-海平面和水位-流量依賴關(guān)系異變情況。相同海平面影響下， 強人類活動影響后甘竹、竹銀和燈籠山站的低水位-低海平面聯(lián)合累積概率曲線明顯抬升， 而高水位-高海平面事件則下降(圖4)， 表明低海平面條件下發(fā)生低水位事件的聯(lián)合累積概率明顯提升， 發(fā)生低水位事件的風險增大， 反之， 高水位事件風險降低； 同流量條件下， 低水位-低流量變化不明顯， 而強人類活動后甘竹至燈籠山站高水位-高流量的聯(lián)合累積概率曲線顯著下降(圖5)。這表明水位與海平面、流量的聯(lián)合依賴關(guān)系發(fā)生明顯改變， 低水位與海平面關(guān)系的敏感度增加， 而高水位與海平面關(guān)系的敏感度下降。低水位與流量關(guān)系的敏感度基本不變， 而高水位與海流量關(guān)系的敏感度明顯下降。

表2 各站強人類活動前、后， 水位-海平面和水位-流量對應(yīng)不同copula 家族的赤池信息量準則(Akaike information criterion, AIC)值 Tab.2 The AIC values for different copulas fitted to water-sea level and water level-river discharge before and after intensive human interventions at different stations

圖4 強人類活動前后海平面-水位聯(lián)合累積頻率曲線 Fig.4 Combined sea-water level recurrence frequency before and after intensive human interventions

圖5 強人類活動前后流量-水位聯(lián)合累積頻率曲線 Fig.5 Combined river discharge-water level recurrence frequency before and after intensive human interventions

由3.1 節(jié)可知， 海平面或流量和水位的線性相關(guān)性較強， 因此， 海平面或流量和水位同時出現(xiàn)的條件概率會存在一個閾值范圍(即圖4、圖5 中等值線彎曲處)， 表明低水位條件僅匹配出現(xiàn)一定范圍的低海平面或低流量， 相應(yīng)地， 磨刀門河口沿程各站高水位出現(xiàn)時， 三灶站的海平面或馬口站的流量也必定在某一個范圍內(nèi)變動， 大于或小于此海平面或流量范圍的值基本不會出現(xiàn)。

3.4 強人類活動前后月均水位異變的風險評估

在給定海平面或流量概率事件的條件下， 基于構(gòu)建的聯(lián)合分布依賴模型可探討水位事件發(fā)生的條件概率， 特別是極端水位事件發(fā)生的概率。采用式(2)計算極端海平面或流量條件下各站水位的條件概率， 即F(x2)low_(SL/D)≤10%和F(x2)high_(SL/D)≥90%， 同時計算海平面事件或流量事件影響下水位閾值的累積概率，即高(F(x1)high_WT≥90%)和低(F(x1)low_WT≤10%)水位事件， 計算結(jié)果如圖6 和圖7 所示， 其中較小的數(shù)值意味著發(fā)生極端水位事件的風險較低的。圖6 和7中兩條實線(紅、藍)分別代表強人類活動前極端低海平面事件(對應(yīng)累積概率小于10%的海平面閾值)和極端高海平面事件(對應(yīng)累積概率大于90%的海平面閾值)條件下水位的累積概率， 而兩條虛線(紅， 藍)對應(yīng)強人類活動后極端海平面事件條件下的水位累積概率。兩種條件分布曲線之間的概率差(即圖中陰影區(qū)域)反映水位對海平面或流量變化敏感程度， 陰影區(qū)域面積的大小表明水位對海平面和流量的敏感性大小因位置(不同站位)和時段(強人類活動前后)的不同而有所差異。

對于各站水位-海平面的依賴性(圖6)， 強人類活動影響后概率差區(qū)域面積明顯縮小， 表明各站月均水位對海平面的依賴性有所減弱。對于甘竹站(圖6a)水位與海平面的遭遇概率， 在高海平面(90%)條件下， 當累計頻率大于0.18 時， 水位曲線左移， 表明水位下降； 當累計頻率小于0.18 時， 水位曲線右移， 表明水位上升。而在低海平面(10%)條件下， 同樣出現(xiàn)這種情況， 當累計頻率大于0.5 時， 水位曲線左移， 表明水位下降； 當累計頻率小于0.5 時， 水位曲線右移， 表明水位上升。三角洲中段的江門站(圖6b)， 在同一海平面條件下， 強人類活動后水位曲線整體左移， 表明水位整體下降。三角洲中下段的竹銀站(圖6c)， 強人類活動后在高海平面(90%)條件下， 當累計頻率大于0.18 時， 水位曲線左移， 同概率下水位下降； 當累計頻率小于0.18時， 水位曲線右移， 表明水位上升。而在低海平面(10%)條件下， 同樣出現(xiàn)這種情況， 當累計頻率大于 0.6 時， 水位曲線左移， 表明水位下降； 當累計頻率小于0.6 時， 水位曲線右移， 表明水位上升。三角洲下段的燈籠山站(圖6d)， 在同一海平面條件下， 強人類活動后水位曲線整體左移， 表明水位整體下降。由此可見， 強人類活動影響后， 各站水位變化幅度變窄。其中甘竹站和竹銀站在同一海平面概率下， 較高的累計頻率對應(yīng)的水位下降， 較低的累計頻率對應(yīng)的水位上升； 江門站和燈籠山站則整體水位下降。

圖6 強人類活動前后海平面極端情況下的水位條件分布 Fig.6 Conditional distribution of water level under extreme cases of sea level before and after intensive human interventions

類似的， 圖7 顯示強人類活動對磨刀門河口沿程各站水位-流量聯(lián)合依賴關(guān)系影響。與水位-海平面聯(lián)合依賴關(guān)系相同， 強人類活動影響后各站的概率差區(qū)域面積明顯縮小， 表明各站月均水位對流量的依賴性均有所減弱。對于水位與流量的遭遇概率， 甘竹站(圖7a)在高流量(90%)條件下， 當累計頻率大于0.18 時， 水位曲線左移， 同條件概率下水位下降； 當累計頻率小于0.18 時， 水位曲線右移， 水位上升。而在低流量(10%)條件下， 水位曲線整體右移， 表明水位上升。三角洲中段的江門站(圖7b)， 在同一海平面條件下， 強人類活動后水位曲線整體左移， 表明水位整體下降。三角洲中下段的竹銀站(圖7c)， 在高流量(90%)條件下， 當累計頻率大于0.2 時， 水位曲線左移， 表明水位下降； 當累計頻率小于0.2 時， 水位曲線右移， 表明水位上升。而在低流量(10%)條件下， 水位曲線整體右移， 表明水位整體抬升。三角洲下段的燈籠山站(圖7d)， 在高流量(90%)條件下， 強人類活動后水位曲線整體左移， 表明水位整體下降； 在低流量(10%)條件下， 水位曲線整體右移， 表明水位整體抬升?？梢钥闯?， 各站在高流量(90%)條件下， 水位基本上呈下降趨勢， 但甘竹站和竹銀站在低累計頻率下對應(yīng)的水位上升； 而低流量(10%)條件下， 甘竹， 竹銀和燈籠山站整體水位上升， 但江門站整體水位呈下降趨勢。

前文已經(jīng)分別分析了海平面或流量對水位的影響， 然而不同的極端水位事件可能是單因素影響， 也可能是海平面與流量共同作用引起的(如復(fù)合極低海平面和低流量事件)， 即三角洲水位同時受海平面和流量的共同影響， 圖8 即為磨刀門河口四個站點的月均水位、海平面和流量耦合條件下， 采用式(3)所計算的各站發(fā)生極端水位事件的條件概率。

圖7 強人類活動前后流量極端情況下的水位條件分布 Fig.7 Conditional distribution of water level under extreme cases of river discharge before and after intensive human interventions

圖8 海平面和流量極端情況下的水位條件分布 Fig.8 Conditional distribution of water level under extreme cases of sea level and river discharge

由圖8 可見， 在上、下邊界的共同影響下， 三角洲大規(guī)模強人類活動之后， 各站概率差區(qū)域面積均有所減小， 表明各站水位對流量及海平面的敏感性均有不同程度的減弱。各站在高海平面和流量(90%)條件下， 水位基本上呈下降趨勢， 但甘竹站和竹銀站在低累積頻率下(＜18%)對應(yīng)的水位上升； 而低海平面和流量(10%)條件下， 甘竹站和竹銀站整體水位上升， 但江門站和燈籠山站整體水位呈下降趨勢。結(jié)合前文， 可以看出圖6 中竹銀、燈籠山站的圖形與圖8中兩站相對應(yīng)的圖形相似， 說明下游站位水位與海平面的依賴關(guān)系更加密切。而圖7 中甘竹、江門站的圖形與圖8 中兩站相對應(yīng)的圖形相似， 說明上游站位流量與水位的依賴關(guān)系更加密切。

給定遭遇概率下， 當水位不斷增大到一定值時海平面或流量將達到一個極小值， 這個水位值定義為有效最大水位(劉曾美等， 2013)。為定量計算強人類活動的影響， 規(guī)定累計頻率趨近于1 時的水位為有效最大水位。計算給定不同海平面和流量遭遇概率下的有效最大水位， 結(jié)果如表3 所示。由表3 可見， 對于相同的水位與海平面、流量的遭遇概率， 強人類活動后其所對應(yīng)的水位值較人類活動前有明顯下降， 特別是在高海平面和高流量概率事件驅(qū)動下， 磨刀門河口沿程水位下降越為顯著。如給定水位與海平面、流量的遭遇概率為0.9 時， 甘竹站降幅最大， 下降了1.24m； 江門站次之， 下降幅度為0.63m； 竹銀站和燈籠山站分別下降了0.25m、0.09m。當給定水位與海平面、流量的遭遇概率為0.5 時， 甘竹站和江門站的水位降幅分別為0.88m 和0.65m， 竹銀站和燈籠山站降幅偏小， 分別下降了0.07m 和0.06m； 而當水位與海平面、流量的遭遇概率為0.1 時， 甘竹站和江門站月平均水位分別下降了0.21m 和0.26m， 竹銀站略為上升0.05m， 燈籠山站下降0.01m。從遭遇概率的角度看， 低海平面和低流量條件下各站月平均水位的下降幅度明顯小于高海平面和高流量條件下月平均水位的變幅。從空間上看， 上游站位月平均水位的下降幅度要明顯大于下游站位。

從河口動力的角度來看， 河口水位受徑潮相互作用影響， 同時也受強人類活動引起的地形邊界改變影響。河口不同區(qū)段水位的主控因素有明顯的洪枯季變化， 上游段枯季水位受徑流和海平面共同控制， 洪季受徑流控制；下游段枯季水位受海平面控制， 而洪季水位主要受徑流和海平面共同控制(蔡華陽等， 2018)。多年大規(guī)模的挖掘河床泥沙， 在很大程度上改變了河道河床演變的自然沖淤特性， 即磨刀門河道河床由總體緩慢淤積轉(zhuǎn)變?yōu)榇蠓葲_刷。大規(guī)模采砂引起磨刀門河床河道大幅度下切和過水面積增大， 必將引起水位發(fā)生明顯變化。由表3 可知， 近幾十年來磨刀門河口的水位普遍下降， 表明挖砂已經(jīng)使得水位變化由過去自然演變的緩慢上升轉(zhuǎn)變?yōu)樗幌陆担?這是性質(zhì)不同、方向不同的變化， 且在很長一段時間內(nèi)不可恢復(fù)。因此， 未來河口三角洲水資源規(guī)劃調(diào)用和河口治理都必須在新水位格局的基礎(chǔ)上進行， 比如基于人類活動后階段的Copula 聯(lián)合概率分布， 為引水口、排污口的位置， 防洪工程及堤防的極端水位設(shè)計、航道開發(fā)、防洪排澇等等提供新的設(shè)計水位參數(shù)。同時， 在不同強度的上下游邊界驅(qū)動下， 強人類活動影響下河口三角洲系統(tǒng)內(nèi)部不同概率水位變化的定量研究結(jié)果， 可為未來河口三角洲高效開發(fā)利用和治理比如綜合整治規(guī)劃、防洪排澇和調(diào)水壓咸策略等用提供重要的基礎(chǔ)參數(shù)。

表3 給定海平面和流量遭遇概率下水位變化范圍 Tab.3 Variation in water level under given sea level and discharge encounter probability

4 結(jié)論

本文選取受強人類活動(如大規(guī)模采砂及圍墾等)影響顯著的珠江磨刀門河口為主要研究區(qū)域， 根據(jù)人類活動歷程將數(shù)據(jù)資料分成兩個時段， 即強人類活動前(1965—1985 年)和強人類活動后(2000—2016年)， 在以馬口站流量、三灶站海平面分別作為上游河流邊界、外海邊界的雙邊界驅(qū)動下， 分別對燈籠山、竹銀、江門、甘竹站的水位序列在月均尺度上采用Copula 方法研究水位與口外海平面、上游流量之間的聯(lián)合分布依賴關(guān)系， 深入分析強人類活動對其聯(lián)合分布依賴關(guān)系的影響， 指出強人類活動后磨刀門河口沿程水位與口外海平面、上游流量之間的依賴關(guān)系減弱， 主要結(jié)論如下：

1) 強人類活動后， 河道地形下切顯著， 外海海平面對甘竹、竹銀、燈籠山站月均水位的線性影響有所增強， 而上游馬口流量對甘竹等站月平均水位的線性影響減弱， 竹銀站和燈籠山站的水位-流量線性相關(guān)系數(shù)由強人類活動前的0.94、0.74 分別減少為0.81、0.56。

2) 強人類活動前后月均水位、海平面和流量的概率分布特征變化明顯。強人類活動后甘竹至燈籠山站水位概率密度曲線明顯左移， 同概率的水位事件較強人類活動前明顯降低， 同時高累級頻率(＞0.18)對應(yīng)的水位下降， 低累級頻率(＜0.18)對應(yīng)的水位上升； 三灶海平面概率密度曲線于強人類活動后右移， 同概率海平面呈上升趨勢； 強人類活動前后流量的威布爾分布變化較小， 流量分布基本不受強人類活動影響， 但累計頻率曲線整體左移， 表明流量整體下降。

3) 強人類活動后各站水位變化幅度變窄， 水位與海平面、流量的聯(lián)合依賴關(guān)系發(fā)生明顯改變， 低水位與海平面關(guān)系的敏感度增加， 而高水位與海平面關(guān)系的敏感度下降。低水位與流量關(guān)系的敏感度基本不變， 而高水位與海流量關(guān)系的敏感度明顯下降。在相同海平面影響下， 強人類活動后發(fā)生低水位事件的風險增大， 反之， 高水位事件風險降低； 而在相同流量影響下， 強人類活動影響后低水位事件的風險基本不變， 而高水位事件的風險顯著下降S。

4) 基于三變量聯(lián)合概率分布分析， 與二元情況比較下發(fā)現(xiàn)下游站位水位與海平面的依賴關(guān)系更加密切， 而上游站位與流量的依賴關(guān)系更加密切。對于不同的水位與海平面、流量遭遇概率， 強人類活動后各站有效最大水位均有明顯下降， 且遭遇概率越大， 相應(yīng)的水位下降幅度越高， 遭遇概率為0.9 時， 甘竹站月均有效最大水位下降幅度達到1.24m， 且上游站位的降幅大于下游站位。