黃銳貞，陳曉宏，何艷虎, 于海霞

(1.中山大學水資源與環(huán)境研究中心，廣州 510275；2.華南地區(qū)水循環(huán)和水安全廣東普通高校重點實驗室，廣州 510275； 3. 中山市火炬高技術產(chǎn)業(yè)開發(fā)區(qū)水利所，廣東 中山 528437)

河口位于海洋和河流的交接地帶，受海洋動力和河流動力的雙重影響，其區(qū)域通常交通便利，是人類生活、生產(chǎn)和貿(mào)易活動頻繁的場所[1]。而咸潮上溯對于河口生態(tài)系統(tǒng)和河口地區(qū)的正常取水都會帶來較大的影響。影響咸潮上溯的因素有很多，其中上游徑流和海洋潮汐這兩個因素尤為重要[2]。而政府間氣候變化[3]專門委員會(IPCC)公布的第五次評估報告(AR5)中的第一工作組報告《氣候變化2013：自然科學基礎》指出，氣候變暖十分明顯，2003-2012年的平均溫度比1850-1900年的平均溫度上升了0.78 ℃。氣溫的上升使得冰川融化和海水溫度升高，從而導致了海平面的上升，進而使得河口地區(qū)的潮位隨之增高，從而加劇了咸潮上溯的嚴重性[4]。河口地區(qū)咸潮上溯一方面令地下水水質(zhì)變壞，影響地下水的抽取利用，另一方面可能引起土地鹽漬化，影響周邊農(nóng)業(yè)生產(chǎn)[5,6]；同時使得河水咸度增加，水質(zhì)變壞，嚴重影響陸地淡水資源。因此，分析咸潮上溯的規(guī)律以及上游來水對河口含氯度的影響及其模擬預測對該地區(qū)的供水安全具有重要意義。

就河口區(qū)咸潮上溯問題，國內(nèi)外已有較為豐富的研究成果。當前國內(nèi)外主要采取一維或多維咸潮入侵數(shù)值模型模擬河口含氯度變化，并耦合河口水動力模型研究鹽度輸送的驅動機制[7-11]。咸潮上溯的影響因素復雜多變，主要包含上游來水、潮汐、地形、風以及海平面變化等[12-14]。在缺乏系統(tǒng)的實測含氯度數(shù)據(jù)及水文資料和大范圍的河口地形下，采用一維或多維數(shù)值模型來研究咸潮上溯對河口飲用水源地和生態(tài)環(huán)境的影響較為困難。而人工神經(jīng)網(wǎng)絡因其具有能夠捕獲要素中有關的非線性關系的能力，可解決因數(shù)據(jù)缺失所帶來的計算問題，因此可用于部分資料缺失情況下的河口區(qū)含氯度的模擬預測。

珠江三角洲地區(qū)是我國重要的經(jīng)濟中心區(qū)域，珠三角河口區(qū)是周邊城市工農(nóng)業(yè)生產(chǎn)和人民生活的重要水源地。磨刀門水道作為珠江的重要入?？谥?，其鹽度變化周期特征與潮汐周期特征基本一致，自有監(jiān)測記錄以來，最嚴重咸潮發(fā)生在2011年12月，受上游來水少和潮周期的影響，磨刀門水道取水口自12月4日起全線24 h含氯度超標，其中珠海的平崗泵站連續(xù)10 d含氯度超標，連續(xù)半月無法取水[15]，咸潮的發(fā)生嚴重影響到人類的基本生活和健康，也是學術界關注的熱點[16]。為此，鑒于珠三角磨刀門河口地區(qū)部分實測含氯度及河口地形數(shù)據(jù)的缺失，本文將運用人工神經(jīng)網(wǎng)絡中的BP神經(jīng)網(wǎng)絡建立珠江磨刀門河口含氯度的預測模型，并分析模型的可靠性；利用訓練好的含氯度預測模型模擬計算90%、97%上游來水頻率與下游5年一遇、10年一遇潮位和20年一遇潮位不同組合下的河口含氯度變化，分析上游來水和下游潮位變化對河口含氯度的影響，進而為河口地區(qū)的供、取水安全以及為“壓咸”而進行的上游水庫水量調(diào)度提供參考價值。

1 研究區(qū)概況及數(shù)據(jù)來源

1.1 研究區(qū)概況

珠江口地區(qū)東起深圳、西迄臺山、北至廣州、南達萬山群島，包括珠江三角洲的大部分地區(qū)和內(nèi)陸架水域。其中磨刀門水道(見圖1)是珠江出海門口之一，是珠江口咸潮活躍的典型河段，其位于廣東省珠海市洪灣企人石，是西江徑流的主要出?？陂T。磨刀門水道全場33.35 km，縣境內(nèi)流域面積177.8 km2，過境流量523 億m3，枯水期漲潮最大流量為9 370 m3/s潮最大，洪水期漲潮最大流量為1.0 萬m3/s，年輸沙量為2 700 萬t。

圖1 磨刀門水道示意圖Fig.1 Schematic diagram of the Modaomen waterway

1999年之前，咸潮主要對中山市西江出?？诘奶怪捩?zhèn)造成影響；1999年以后，咸潮大幅上溯，開始對城區(qū)及周邊地區(qū)造成影響[17]；到了2005年咸潮是近30年來最嚴重的，在中山市24個鎮(zhèn)區(qū)中，供水受咸潮不同程度影響的鎮(zhèn)區(qū)已達到18個[18]；而到了2011年，下游咸潮兇猛，平崗泵站已破2005年的咸情記錄[15]。

1.2 數(shù)據(jù)及來源

本文采用磨刀門水道中三水站+馬口站的前一天日平均徑流量、燈籠山站日平均潮位和平崗泵站日平均含氯度作為研究數(shù)據(jù)。數(shù)據(jù)均在2005年1、2、3、9、10、11和12月的實測資料中選取，總樣本個數(shù)為173個，用于神經(jīng)網(wǎng)絡模型的訓練和驗證。其中1、2、3、9、10月的數(shù)據(jù)(某些天數(shù)的數(shù)據(jù)有缺失)作為第1組數(shù)據(jù)，樣本個數(shù)為134個，其用于模型的訓練；而11和12月的數(shù)據(jù)(某些天數(shù)的數(shù)據(jù)有缺失)則作為第2組數(shù)據(jù)，樣本個數(shù)為39個，其用于模型的驗證。

以上述三水站+馬口站的前一天日平均徑流量和燈籠山日平均潮位作為本文神經(jīng)網(wǎng)絡模型中的輸入，平崗泵站日平均含氯度作為模型的輸出。而由于各輸入輸出數(shù)據(jù)之間所屬范圍差別較大，直接使用觀測數(shù)據(jù)是網(wǎng)絡收斂緩慢甚至無法收斂，造成模型無法模擬，因此，為了使程序運行時加快收斂，對數(shù)據(jù)中的平均徑流量數(shù)據(jù)和含氯度數(shù)據(jù)均除以1 000，而后使用premnmx函數(shù)對所有數(shù)據(jù)進行歸一化處理；模型訓練和驗證后，對輸出數(shù)據(jù)運用postmnmx函數(shù)進行反歸一化處理[19]。

2 BP神經(jīng)網(wǎng)絡方法

近年來人工神經(jīng)網(wǎng)絡(ANN)在水文分析和水文預測中的應用越發(fā)廣泛[20]，而且取得良好效果。而在這些應用中，BP神經(jīng)網(wǎng)絡是較為有效的方法，且大部分采用三層網(wǎng)絡結構以解決水文研究中的問題[21]。本文采用的3層神經(jīng)網(wǎng)絡結構如圖2所示，該網(wǎng)絡結構帶有1個輸入層、1個帶非線性轉換函數(shù)的隱含層和1個帶線性轉換函數(shù)的輸出層。本文對單隱含層的BP神經(jīng)網(wǎng)絡進行說明[22]，假設BP神經(jīng)網(wǎng)絡的輸入矢量x∈Rn，其中x=(x0,x1,…,xn-1)T；設隱含層有n1個神經(jīng)元，其輸出為x′∈Rn1，其中x′=(x′0,x′1,…,x′n1-1)T；輸出矢量y∈Rm，其中y=(y0,y1,…,ym-1)T。從輸入層到隱含層，其權值為wij，閾值為θj；從隱含層到輸出層，其權值為w′jk，閾值為θ′k。各層神經(jīng)元的輸出為：

(1)

圖2 三層BP神經(jīng)網(wǎng)絡Fig.2 Three layers of BP neural network

訓練中神經(jīng)網(wǎng)絡會隨機指定每個神經(jīng)元的初始權值和閥值，然后把隨機權重下得到的模型模擬值與目標值作比較，通過不斷調(diào)整各層的權值和閥值，使得網(wǎng)絡輸出的模擬值與目標值之間的均方誤差達到所給定的期望值。

3 結果分析

3.1 BP神經(jīng)網(wǎng)絡的模擬與驗證

運用MATLAB軟件對模型進行編程，模型中的一層隱含層采用雙曲正切S形傳輸函數(shù)，輸出層則使用線性傳輸函數(shù)，其中訓練函數(shù)為traingdm算法。在模型訓練階段，根據(jù)第1組數(shù)據(jù)的輸入，運用BP神經(jīng)網(wǎng)絡模型計算的平崗泵站日平均含氯度與實測值進行比較，在未達到理想的結果前不斷調(diào)整模型中的各個參數(shù)，直到計算值與實測值之間的均方誤差達到所設定的期望值。之后采用第2組數(shù)據(jù)對模型進行驗證。

設定模型的訓練步數(shù)為50 000步，期望訓練精度為0.05。經(jīng)過多次的參數(shù)調(diào)整，最后確定隱含層的節(jié)點數(shù)為16，模型訓練后在第48 562步時就達到所設計的最小均方差(0.05)，可見隱含層節(jié)點數(shù)在這個設定下有利于網(wǎng)絡的收斂。

圖3為設定好的網(wǎng)絡模型計算的日平均含氯度與實測數(shù)據(jù)之間的比較，可以看出計算值與實測值符合良好，其相關系數(shù)為0.867 8。

圖3 BP神經(jīng)網(wǎng)絡模型訓練后的結果Fig.3 Curves of observed chlorinity and its calibration results

而后運用第2組數(shù)據(jù)對訓練后的模型進行仿真驗證，結果如圖4所示?？梢钥闯瞿Ｐ偷尿炞C效果良好，其相關系數(shù)達到0.84。為進一步驗證模型的適用性，利用訓練好的模型對2009年10-12月的含氯度進行模擬，實測值與模擬值之間的相關系數(shù)為0.82，說明了模型較好的適用性。

圖4 BP神經(jīng)網(wǎng)絡模型計算與實測數(shù)據(jù)的比較Fig.4 Curves of observed chlorinity and its validation results

利用BP神經(jīng)網(wǎng)絡訓練和驗證的河口含氯度計算結果與實際含氯度之間的相關系數(shù)均在0.82以上， BP神經(jīng)網(wǎng)絡能夠較好地反映出上游徑流量、河口咸潮期潮位與含氯度三者之間的非線性關系。因此，根據(jù)徑流量和潮位數(shù)據(jù)，即能夠較準確地模擬計算出相應的含氯度。

3.2 徑流量的變化對含氯度的影響

運用上述基于BP神經(jīng)網(wǎng)絡的珠江河口含氯度預測模型研究受調(diào)節(jié)的上游水電站的徑流量變化對河口含氯度的影響。對上游徑流量分別按升序和降序進行排列，由于潮汐變化具有顯著的周期性，不同潮周期間潮位過程相近，因此計算時仍將潮位作為下游邊界條件[21]。利用上述訓練好的BP神經(jīng)網(wǎng)絡模型進行模擬計算，結果如圖5所示。

圖5 上游徑流量升序與降序時含氯度的變化Fig.5 Chlorinity changes under upstream runoff sequences of ascending and descending

由圖5可以看出，上游徑流量與河口含氯度呈負相關關系：當徑流量按升序排列時，含氯度呈下降趨勢；當徑流量按降序排列時，含氯度呈上升趨勢。這反映了上游徑流量越大，河口含氯度越小，體現(xiàn)出上游徑流量的“壓咸”作用，與Liu[12]的研究成果一致，印證了咸潮入侵時期上游來水對減輕河口含氯度確保供水安全具有至關重要作用。表1為上游徑流量的原始值、升序、降序排列對河口含氯度影響的對比。

表1 不同排序的徑流量對含氯度的影響Tab.1 Effects of different sort of runoff to chlorinity

由表1可知，當徑流量以降序排列時，最大含氯度為1 225 mg/L，有34 d的濃度是低于250 mg/L，即水質(zhì)達標[12]，這比原始徑流量下的達標天數(shù)多了8 d，也低于其最大含氯度2 029 mg/L。

對第2組數(shù)據(jù)(驗證數(shù)據(jù)，共39 d)的上游徑流量分別按升序和降序進行排列，下游潮位不變，運用上述訓練好的BP神經(jīng)網(wǎng)絡進行計算?？梢钥闯?，對上游徑流量進行排序泄水時，其水質(zhì)達標情況均比實際徑流量好；當上游徑流量以降序形式進行泄水時，其結果最好，徑流量降序時水質(zhì)達標的天數(shù)比實際情況的要多8 d，其最大含氯度為1 225 mg/L，也比實際情況的最大含氯度2 029 mg/L要小，這可能是因為對于受咸潮入侵影響的河口區(qū)，前期徑流量較大時，會使河口段本底含氯度相對較小，這在一定程度上會減輕前期含量度對后期的影響。

3.3 不同來水頻率和潮位設計方案下的含氯度模擬預測

基于1989-2005年上游三水站+馬口站前一天徑流序列和1971-2004年燈籠山逐時潮位數(shù)據(jù)，設計以下6種邊界條件情形：90%來水頻率、5年一遇潮位；90%來水頻率、10年一遇潮位；90%來水頻率、20年一遇潮位；97%來水頻率、5年一遇潮位；97%來水頻率、10年一遇潮位；97%來水頻率、20年一遇潮位。運用訓練好的BP神經(jīng)網(wǎng)絡模型對平崗泵站日平均含氯度進行模擬計算，以分析不同來水頻率和潮位邊界情形下的含氯度變化，結果如表2所示。

表2 模型計算結果 mg/L

由表2可知，以上游90%來水頻率為例，在5年一遇潮位情形下，方案設計的結果(366.1 mg/L)較2005年實際河口含氯度(140.25 mg/L)[15]大，這是因為實際數(shù)據(jù)的潮位(-0.215 m)比方案預測中的低，海水入侵程度較低，所以其含氯度比方案設計計算的結果要??；在十年一遇潮位(0.11 m)情形下，在2005年的實測數(shù)據(jù)中找出與其接近的數(shù)據(jù)，得到某天的平均潮位為0.112 m，平均徑流量為2 412 m3/s，河口含氯度362.17 mg/L，這組實測數(shù)據(jù)與方案90%來水頻率+10年一遇潮位的數(shù)據(jù)與計算結果(366.1 mg/L)相似，說明本文神經(jīng)網(wǎng)絡模型的模擬效果較好?？梢钥闯?，河口設計潮位越高，上游來水越小，河口含氯度越大。

當上游來水頻率為90%時，5年一遇潮位下的含氯度與10年一遇的含氯度相差10 mg/L左右；而當上游來水頻率為97%時，5年一遇潮位下的含氯度與10年一遇的含氯度相差15 mg/L。由此可以看出當上游徑流量較小時，下游潮位對河口含氯度的影響會比徑流量較大的影響要大。10年一遇潮位與20年一遇潮位亦是如此，表明當上游徑流量較小時，下游潮位對于河口含氯度的影響程度較大。章文[14]等的研究也指出，磨刀門水道的含氯度在極端枯水期很大程度上受下游潮位變化幅度的影響。

同一設計潮位方案下，特枯來水年份(97%來水頻率)平崗泵站日平均含氯度較正常來水年份(90%來水頻率)平均增幅達到3.3倍。尤其是特枯來水年份，平崗泵站10年一遇潮位下日平均含氯度達到了1 500 mg/L。說明枯水期上游來水的減少較大增加了咸潮上溯的風險，進而造成了河口地區(qū)正常供取水的困難。

4 結 語

(1)本文建立的基于BP神經(jīng)網(wǎng)絡的含氯度模擬模型對平崗泵站的日平均含氯度模擬訓練和驗證的結果與實際含氯度之間的相關系數(shù)均在0.82以上，模型體現(xiàn)了較好的適用性，可用于磨刀門河口區(qū)含氯度的模擬預測。

(2)不同流量過程條件下日平均含氯度的模擬表明，采用水庫泄水來抑制潮汐河口水源地咸潮上溯，在泄水總量不變的情況下，按照徑流量降序來泄水更能有效地緩解咸潮上溯對水源地的影響。這對于河口地區(qū)的供水和取水安全及水庫泄水調(diào)度均具有重要的意義。

(3)不同來水頻率和潮位設計方案下的含氯度預測模擬結果表明：上游徑流量較小時，下游潮位對河口含氯度的影響會比徑流量較大的影響要大；同一設計潮位方案下，特枯來水年份(97%來水頻率)平崗泵站日平均含氯度較正常來水年份(90%來水頻率)平均增幅達到3.3倍。

值得指出的是，本文僅選取了影響河口含氯度的兩個主要因素：上游來水和下游潮位，以分析上游來水變化及其與不同潮位的組合對網(wǎng)河區(qū)河口含氯度的影響。實際上，除了上游來水和下游潮位，河口含氯度還受到風向，河道地形變化，下游支流的交互作用以及其他因素的影響。例如，磨刀門水道下游的洪灣和賀洲水道關閉可導致20%咸潮入侵的減少[23]；在小潮或是落潮階段，更多咸潮由洪灣水道涌向磨刀門水道，這是此階段磨刀門水道上游咸潮入侵達到最大化的內(nèi)在動力機制[24]。因此，河口含氯度的變化呈現(xiàn)出高度的非線性[12]，咸潮入侵的動力機制及其對河口含氯度的影響機理較為復雜，仍需進一步深入研究。

□

[1] 楊莉玲. 河口鹽水入侵的數(shù)值模擬研究[D]. 上海：上海交通大學,2007.

[2] Lian JJ, He W, Ma C, et al. Guarantee rate of freshwater in a river mouth intruded by saltwater with respect to the joint impact of runoff and tide [J]. Journal of Hydroinformatics, 2015,17(6):917-929.

[3] IPCC: Climate Change 2013. The Physical Science Basis[Z].

[4] Ketabchi H, Mahmoodzadeh D, Ataie-Ashtiani B et al. Sea-level rise impacts on seawater intrusion in coastal aquifers: review and integration[J]. Journal of Hydrology, 2016,535:235-255.

[5] Colon-Rivera RJ, Feagin RA, West JB et al. Hydrological modification, saltwater intrusion, and tree water use of a Pterocarpus officinalis swamp in Puerto Rico [J]. Estuarine Coastal and Shelf Science, 2014,147:156-167.

[6] Vijay R, Mohapatra PK. Hydrodynamic assessment of coastal aquifer against saltwater intrusion for city water supply of Puri, India [J]. Environmental Earth Sciences, 2016,75:588.

[7] 孔 蘭,陳曉宏,杜 建,等. 基于數(shù)學模型的海平面上升對咸潮上溯的影響[J]. 自然資源學報,2010,(7):1 097-1 104.

[8] 史英標,李若華,姚凱華. 錢塘江河口一維鹽度動床預報模型及應用[J]. 水科學進展,2015,(2):212-220.

[9] Kernell G R. Selected streamflow statistics for stream gaging stations in Northeastern Maryland[R]. U.S. Department of the Interior. U.S. Geological Survey, Open-file Report, 2006.

[10] Kernell G R. The national streamflow statisics program: a computer program for estimating streamflow stistics for ungaged sites[M]. U.S. Geological Survey, 2007.

[11] BMJ Langland, JP Raffensperger, DL Moyer, et al. Changes in streamflow and water quality in Selected Nontidal basins in the Chesapeake Bay Watershed, 1985-2004[R]. U.S. Geological Survey Scientific Investigations Report, 2006.

[12] Liu BJ, Yan SL, Chen XH, et al. Wavelet analysis of the dynamic characteristics of saltwater intrusion - a case study in the Pearl River Estuary of China[J]. Ocean & Coastal Management, 2014,95:81-92.

[13] Yuan R, Zhu JR. The effects of dredging on tidal range and saltwater intrusion in the Pearl River Estuary [J]. Journal of Coastal Research, 2015,31(6):1 357-1 362.

[14] 章 文,劉丙軍,辛彥博，等. 珠江河口區(qū)鹽度變化周期特征分析----以磨刀門水道為例[J]. 熱帶地理,2013,33(1):28-33.

[15] 黃建華,黃庭輝,陳治家. 中山珠海咸潮兇猛[N]. 廣州日報,2011-12-20(A14).

[16] 林偉波. 甌江口鹽水入侵預測模式研究[J]. 中國農(nóng)村水利水電,2013,(3): 24-33.

[17] 胡 溪. 珠江口磨刀門水道咸潮入侵數(shù)值模擬研究[D]. 北京：清華大學, 2010.

[18] 周 琳,汪 進,萬 飛,等. 面對咸潮----市供水有限公司努力抗咸保供水[N]. 中山日報畫刊,2006-12-01(第6期27版).

[19] 朱 凱. 精通MATLAB神經(jīng)網(wǎng)絡[M]. 北京:電子工業(yè)出版社,2010.

[20] 胡鐵松,袁 鵬,丁 晶. 人工神經(jīng)網(wǎng)絡在水文水資源中的應用[J]. 水科學進展,1995,(1): 76-82.

[21] 劉國東,丁 晶. BP網(wǎng)絡用于水文預測的幾個問題探討[J]. 水利學報,1999,(1): 66 -71.

[22] 黃 麗. BP神經(jīng)網(wǎng)絡算法改進及應用研究[D]. 重慶：重慶師范大學,2008.

[23] 胡 溪,毛獻忠. 珠江口磨刀門水道咸潮入侵規(guī)律研究[J]. 水利學報,2012,(5): 529 -536.

[24] Gong WP, Wang YP, Jia JJ. The effect of interacting downstream branches on saltwater intrusion in the Modaomen Estuary, China[J]. Journal of Asian Earth Sciences， 2012,45:223-238.