唐啟邦 歐素英 蔡華陽(yáng) 楊清書

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

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

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

1 研究區(qū)域概況

本文研究區(qū)域?yàn)橹榻サ堕T河口， 包括西江干流水道-磨刀門水道(如圖1 所示)。西北江水沙經(jīng)思賢滘口分水分沙后， 約78%—84%的徑流量通過西江干流水道進(jìn)入西北江三角洲。磨刀門作為珠江八大口門中入海徑流量和輸沙量最大的河口， 約36%的珠江徑流量通過磨刀門進(jìn)入南海， 是珠江最重要的泄洪通道。磨刀門河口潮差較小， 口門三灶站多年平均潮差小于1m。因此， 磨刀門河口徑流動(dòng)力強(qiáng)而潮汐動(dòng)力弱， 是典型的河優(yōu)型河口?；谏嫌务R口水文控制站1965—2016 年的月均流量統(tǒng)計(jì)資料， 可得洪季(4—9 月)歷史最大流量為 29000m3/s， 平均流量為10810m3/s； 枯季(10 月—次年3 月)歷史最大流量為14500m3/s， 平均流量為3428m3/s； 年均歷史最大流量為9978m3/s， 平均流量為7119m3/s?？陂T三灶站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 世紀(jì)50 年代開始， 磨刀門河口先后經(jīng)歷了聯(lián)圍筑閘、灘涂圍墾、航道整治規(guī)劃以及大規(guī)模的河道無序采砂等強(qiáng)人類活動(dòng)。特別是20 世紀(jì)90 年代磨刀門完成口門整治后， 口門人為向?？焖傺由旒s16km， 使得原本具有泄洪納潮功能的寬闊內(nèi)海灣消失， 河口環(huán)境轉(zhuǎn)變?yōu)槭苋斯?dǎo)堤控制的河道， 口門迅速向海延伸， 延伸速度超越自然過程變化， 口門邊界條件發(fā)生明顯異變。西江上游流域方面， 建有大藤峽等骨干水源工程， 以西江龍灘及大藤峽等水庫(kù)為骨干構(gòu)建水資源調(diào)配體系， 直接改變磨刀門河口的來水來沙條件(張旭等， 2013)。而對(duì)磨刀門河口地形影響最為明顯的是大規(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， 河床沖淤演變趨勢(shì)由原來淤積為主轉(zhuǎn)變?yōu)轱@著沖刷(賈良文等， 2002)。自禁沙令頒布后， 2000 年后河道采砂才逐漸減少。

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

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

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

2.2 Copula 方法

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

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

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

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

3 結(jié)果分析

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

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

表1 強(qiáng)人類活動(dòng)前后水位-海平面、水位-流量的線性相關(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)計(jì)量 Fig.2 Chi-square statistics of marginal distributions for different variables at different stations

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

圖3 強(qiáng)人類活動(dòng)前后磨刀門河口沿程各站變量邊緣分布累積頻率曲線 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 準(zhǔn)則選擇最佳的Copula函數(shù)， 結(jié)果如表2 所示。選取AIC 值最小的Gaussian函數(shù)作為水位-海平面和水位-流量的最佳擬合Copula聯(lián)合分布函數(shù)用于后續(xù)依賴關(guān)系的異變研究。

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

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

表2 各站強(qiáng)人類活動(dòng)前、后， 水位-海平面和水位-流量對(duì)應(yīng)不同copula 家族的赤池信息量準(zhǔn)則(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 強(qiáng)人類活動(dòng)前后海平面-水位聯(lián)合累積頻率曲線 Fig.4 Combined sea-water level recurrence frequency before and after intensive human interventions

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

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

3.4 強(qiáng)人類活動(dòng)前后月均水位異變的風(fēng)險(xiǎn)評(píng)估

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

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

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

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

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

圖7 強(qiáng)人類活動(dòng)前后流量極端情況下的水位條件分布 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ī)模強(qiáng)人類活動(dòng)之后， 各站概率差區(qū)域面積均有所減小， 表明各站水位對(duì)流量及海平面的敏感性均有不同程度的減弱。各站在高海平面和流量(90%)條件下， 水位基本上呈下降趨勢(shì)， 但甘竹站和竹銀站在低累積頻率下(＜18%)對(duì)應(yīng)的水位上升； 而低海平面和流量(10%)條件下， 甘竹站和竹銀站整體水位上升， 但江門站和燈籠山站整體水位呈下降趨勢(shì)。結(jié)合前文， 可以看出圖6 中竹銀、燈籠山站的圖形與圖8中兩站相對(duì)應(yīng)的圖形相似， 說明下游站位水位與海平面的依賴關(guān)系更加密切。而圖7 中甘竹、江門站的圖形與圖8 中兩站相對(duì)應(yīng)的圖形相似， 說明上游站位流量與水位的依賴關(guān)系更加密切。

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

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

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

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

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

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

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