珠江河口磨刀門枯水期鹽度統(tǒng)計模型

王 彪，朱建榮，呂愛琴

(1.上海市環(huán)境科學(xué)研究院，上海 200233； 2.上海水環(huán)境模擬與水生態(tài)修復(fù)工程技術(shù)研究中心，上海 200233；3.長三角區(qū)域生態(tài)環(huán)境聯(lián)合研究中心，上海 200233； 4.華東師范大學(xué)河口海岸學(xué)國家重點實驗室，上海 200062； 5.廣東省水文局佛山分局，廣東 佛山 528000)

珠江是我國第二大河，磨刀門為珠江河口八大入海口門之一，是主要的徑流和泥沙出口，也是珠海和澳門主要水源地。20世紀(jì)90年代以來，珠江三角洲地區(qū)咸潮活動越來越頻繁，持續(xù)時間增加，上溯影響范圍越來越大，強(qiáng)度趨于嚴(yán)重，發(fā)展態(tài)勢不容樂觀[1]，嚴(yán)重影響了珠海和澳門的供水安全。

鹽水入侵是河口地區(qū)普遍存在的一種徑流、潮汐和風(fēng)力等動力因子綜合作用形成的一種自然現(xiàn)象。以往對珠江河口鹽水入侵的大量研究表明，潮汐和徑流是影響珠江河口鹽水入侵的主要因子[2-6]。楊干然[2]根據(jù)海、河兩大動力體系的強(qiáng)弱程度將磨刀門劃分為徑流型河口，認(rèn)為鹽度分布在洪季呈弱混合型，而枯季為緩混合型。李素瓊等[3]則指出磨刀門無論枯、汛期咸淡水高度分層，鹽水楔活動明顯。樂燦等[4]基于珠江河口2016 年枯季同步觀測鹽度資料，總結(jié)了枯季八大口門同步鹽度垂向分布和鹽淡水混合特征，指出一般潮汐動力越強(qiáng)，鹽度垂向分布的拐點位置越高。徐俊峰等[7]對珠江三角洲網(wǎng)河區(qū)一維水沙、含氯度輸移數(shù)學(xué)模型和河口區(qū)二維水沙、含氯度輸移數(shù)學(xué)模型進(jìn)行聯(lián)解，建立了珠江三角洲及河口區(qū)整體數(shù)學(xué)模型。為了更深入研究珠江河口的水流運動和鹽度輸移的三維特征，不少學(xué)者基于動力學(xué)過程的三維數(shù)學(xué)模型開展了大量研究[8-13]。Wong等[8]基于POM(Princeton ocean model)分析了珠江河口鹽度梯度形成的羽狀鋒面的動力機(jī)制以及與冬季枯水期羽狀鋒面相關(guān)的小尺度環(huán)流。林若蘭等[11]采用MIKE3 構(gòu)建珠江東四口門三維水動力數(shù)值模型，對枯水期各風(fēng)向下珠江口的水動力進(jìn)行模擬，分析北風(fēng)、東北風(fēng)和東風(fēng)對河口漲落潮流速、鹽度分布、潮通量的影響。呂紫君等[12]基于SCHISM (semi-implicit cross-scale hydroscience integrated system model)建立了磨刀門三維鹽度數(shù)值模型，利用該模型并結(jié)合鹽通量分解方法，研究了磨刀門水道鹽淡水混合特征及各驅(qū)動力對鹽度通量變化的貢獻(xiàn)作用。此外，陳水森等[14]基于經(jīng)驗?zāi)Ｐ头治隽四サ堕T水道含氯度及咸潮入侵最大距離與上游徑流、河口鹽度的關(guān)系；路劍飛等[15]利用引入滯后因子的徑向基函數(shù)神經(jīng)網(wǎng)絡(luò)模型，以珠江口磨刀門水道為例,研究建立了鹽度多步預(yù)測模型；周凡涵等[16]基于小波-神經(jīng)網(wǎng)絡(luò)構(gòu)建了珠江口門區(qū)磨刀門水道的鹽水入侵預(yù)報模型，對磨刀門鹽水入侵進(jìn)行了預(yù)見期1～3 d 的預(yù)測應(yīng)用?？傮w上采用統(tǒng)計學(xué)經(jīng)驗?zāi)Ｐ蛯χ榻涌邴}水入侵的研究相對較少，原因是統(tǒng)計學(xué)模型主要針對特定站點、特定指標(biāo)(如日均鹽度、入侵距離等)進(jìn)行統(tǒng)計學(xué)意義上的反演和預(yù)測，而基于動力方程的數(shù)學(xué)模型則能深入反演和預(yù)測河口鹽度的時空變化。但統(tǒng)計學(xué)經(jīng)驗?zāi)Ｐ蛯A(chǔ)數(shù)據(jù)要求較少，易于快速實現(xiàn)，在實際管理應(yīng)用中仍有較強(qiáng)實用性。本文針對潮汐和徑流兩個主要動力因子，采用最小二乘擬合方法研究構(gòu)建枯水期磨刀門鹽度統(tǒng)計模型。

1 研究數(shù)據(jù)和方法

珠江河口磨刀門地理位置如圖1所示，影響磨刀門水道鹽水入侵最主要的兩個動力因子是潮汐和徑流。本文基于該動力機(jī)制，應(yīng)用平崗泵站2005年1月1日至3月31日(時段1)和2005年10月1日至2006年2月23日(時段2)兩個時段的鹽度和潮位資料，以及同期的磨刀門上游西江的高要水文站和北江的石角水文站的合計徑流量，建立鹽度與潮汐和徑流的統(tǒng)計關(guān)系模型。

圖1 珠江河口磨刀門地理位置Fig.1 Location of Modaomen in Pearl River Estuary

基于統(tǒng)計學(xué)方法構(gòu)建鹽水入侵統(tǒng)計模型，主要針對鹽度的特征值，如鹽度日平均值、日最大值和日最小值，挖掘其與相關(guān)動力因子之間的量化響應(yīng)關(guān)系。本文主要構(gòu)建日平均鹽度與潮汐、徑流動力因子之間的量化響應(yīng)關(guān)系，日平均鹽度通過對平崗泵站日內(nèi)逐時鹽度進(jìn)行算術(shù)平均得到。潮汐動力因子以潮差進(jìn)行表征，珠江河口潮汐為不規(guī)則潮，大潮期間為不規(guī)則半日潮，小潮期間為不規(guī)則全日潮，本文取日內(nèi)最大的潮差作為當(dāng)日潮汐動力強(qiáng)度。由于徑流對鹽度的影響不僅與近期的徑流量大小有關(guān)，其前期的徑流量大小對河口鹽度也存在累積影響，因而分析徑流動力因子對鹽度的影響時，應(yīng)考慮其前期的徑流量情況，考慮潮汐的半月周期，本文取前15 d的平均徑流量作為徑流動力因子強(qiáng)度。處理前后的鹽度、潮汐、徑流如圖2所示，原始鹽度、潮位數(shù)據(jù)時間分辨率為1 h，徑流量數(shù)據(jù)時間分辨率為1 d。原始數(shù)據(jù)表明，平崗泵站鹽度和水位除了有日內(nèi)變化特征外，具有顯著的半月變化特征；徑流量與鹽度之間總體上呈現(xiàn)徑流量大、鹽水入侵弱和徑流量小、鹽水入侵強(qiáng)的特征。

(a) 時段1流量

(b) 時段1潮位潮差

(c) 時段1鹽度

(d) 時段2流量

(e) 時段1潮位、潮差

(f) 時段2鹽度圖2 徑流、潮汐和鹽度隨時間變化Fig.2 Temporal variations of runoff, tide, and salinity

由于外海高鹽水入侵河口是一個復(fù)雜的動力過程，潮汐強(qiáng)度峰值與鹽度峰值往往不會同步出現(xiàn)，而是存在一定的相位差。圖2表明，日均鹽度和潮差雖然都有明顯的半月周期變化，但二者峰值確實存在明顯的相位差。尹小玲等[17]參考珠江河口磨刀門宏觀幾何條件，建立了概化河口數(shù)值模型，分析指出由于混合過程時間效應(yīng)，咸水運動調(diào)整過程呈現(xiàn)非線性特征，鹽分向陸輸運量在小潮后2～3 d達(dá)到最大。Wang等[18]通過珠江河口三維鹽水入侵?jǐn)?shù)值模型得出，該相位差動力成因主要是由于下游河口處洪灣水道汊道(圖1中通向澳門的水道)的特殊地形與風(fēng)和潮汐動力相互作用形成的。為使構(gòu)建的鹽度統(tǒng)計模型具有更好的適用性，需要對鹽度和潮差進(jìn)行相位校正處理。通過潮差和鹽度的相關(guān)分析，發(fā)現(xiàn)將潮差進(jìn)行3.5 d的相位校正處理后，二者相關(guān)系數(shù)達(dá)到最大。因此，下文統(tǒng)計建模中當(dāng)日鹽度對應(yīng)的潮差為當(dāng)日后第3和第4天的潮差均值?；谌掌骄}度、第3和第4天的平均潮差和前15 d平均徑流量數(shù)據(jù)(下文分別簡稱為鹽度、潮差和徑流)，本文采用統(tǒng)計學(xué)方法構(gòu)建珠江河口磨刀門鹽度統(tǒng)計模型。

2 模型構(gòu)建

珠江河口鹽水入侵主要受潮汐和徑流的作用，且潮汐引起的半月時間尺度變化更為明顯，故采用最小二乘法先擬合出鹽度和潮差之間的關(guān)系式，據(jù)此預(yù)估出潮差對鹽度影響；在此基礎(chǔ)上，觀測鹽度去除潮差影響貢獻(xiàn)，可得到一個鹽度差，這個差值可認(rèn)為是徑流變化引起，通過擬合方法可找出鹽度差與徑流的關(guān)系式。將潮差-鹽度擬合關(guān)系式疊加徑流-鹽度差擬合關(guān)系式，最終可得到一個綜合關(guān)系式，此即本文構(gòu)建鹽度統(tǒng)計模型的基本思路。

圖3(a)給出了平崗泵站的鹽度-潮差散點分布情況，可見二者為非線性關(guān)系，總體可認(rèn)為是指數(shù)關(guān)系。因此采用指數(shù)關(guān)系對鹽度和潮差進(jìn)行擬合，通過最小二乘法得出擬合關(guān)系為

(a) 鹽度與潮差

(b) 歸一化鹽度差與徑流圖3 鹽度與潮差及歸一化鹽度差與徑流的最小二乘擬合Fig.3 Least square fitting curves of salinity and tidal range as well as normalized salinity difference and runoff

SFT0=0.000 429 4e4.996HTR

(1)

式中：SFT0為鹽度受到潮差影響的潮差擬合鹽度；HTR為潮差。圖3(a)中的曲線即為擬合曲線，潮差-鹽度的擬合優(yōu)度R2為0.90?；诖顺辈?鹽度擬合關(guān)系式，可計算出潮差隨時間變化對應(yīng)的潮差擬合鹽度變化，其與實測鹽度對比見圖4。

圖4 實測鹽度、鹽度差與擬合結(jié)果對比Fig.4 Comparison of measured and fitting results of salinity and salinity difference

實測鹽度中包含了潮差和徑流的共同作用，在得到式(1)后，實測鹽度減去潮差擬合鹽度得到的鹽度差可認(rèn)為是徑流的作用。若對此鹽度差與徑流直接進(jìn)行擬合建模，會發(fā)現(xiàn)擬合效果不理想。分析其原因發(fā)現(xiàn)，此鹽度差還受到潮形的影響，相同徑流下大潮期間由于河口總體鹽度大，對應(yīng)的鹽度差也相對較大；小潮期間對應(yīng)的河口鹽度總體較低，鹽度差也相對較小。為了去除潮形不同引起的徑流對應(yīng)的鹽度差差異，本文通過對鹽度差進(jìn)行潮形歸一化處理，得到歸一化鹽度差：

SSD=SD/SFT0

(2)

其中

SD=S-SFT0

式中：SSD為歸一化實測鹽度差；SD為徑流作用對應(yīng)鹽度差；S為鹽度實測值。

通過分析歸一化實測鹽度差與徑流之間的關(guān)系，發(fā)現(xiàn)二者具有反比例函數(shù)關(guān)系。因此，采用反比例函數(shù)進(jìn)行歸一化實測鹽度差與徑流擬合，得到徑流-鹽度差擬合函數(shù)：

(3)

式中：SQSD0為徑流作用對應(yīng)的歸一化擬合鹽度差；Q為徑流。擬合曲線如圖3(b)所示，徑流-歸一化實測鹽度差的擬合優(yōu)度R2為0.71。在此基礎(chǔ)上，去歸一化便可得到徑流擬合鹽度差：

SFQ0=SQSD0SFT0

(4)

實測徑流作用對應(yīng)鹽度差SD與徑流擬合鹽度差SFQ0對比如圖4所示。綜合潮差-鹽度擬合關(guān)系(式(1))和徑流-鹽度差之間關(guān)系(式(3)和(4))，可得到河口鹽度統(tǒng)計模型表達(dá)式為

(5)

式中SiniF為鹽度統(tǒng)計模型計算的鹽度，其與實測鹽度對比情況見圖4。

3 修正鹽度統(tǒng)計模型

鹽度統(tǒng)計模型中進(jìn)行潮差-鹽度關(guān)系擬合時，實測鹽度資料中包含了徑流的貢獻(xiàn)。通過式(3)和(4)可以計算出徑流對鹽度的貢獻(xiàn)，將實測鹽度S減去徑流對鹽度的貢獻(xiàn)SFQ0，得到濾去徑流影響的實測鹽度SRFQ0，這個鹽度與潮差關(guān)系更為緊密。基于SRFQ0重新與潮差進(jìn)行擬合，能進(jìn)一步提高模型精度。因此基于指數(shù)關(guān)系式用最小二乘法進(jìn)行重新擬合，可得到潮差-鹽度修正擬合關(guān)系式：

SFT1=0.001 127e4.372HTR

(6)

式中SFT1為修正潮差擬合鹽度。

擬合曲線如圖5(a)所示，兩者擬合優(yōu)度R2為0.92。濾去徑流影響的實測鹽度與修正潮差擬合鹽度的時間序列對比見圖6。由于時段2第130天附近徑流擬合鹽度偏高，因此濾去該徑流影響的實測鹽度出現(xiàn)負(fù)值，導(dǎo)致該時段擬合效果不好，但整體上修正潮差擬合鹽度與實測鹽度更為一致。

(a) 鹽度與潮差

(b) 歸一化鹽度差與徑流圖5 鹽度與潮差及歸一化鹽度差與徑流的修正擬合Fig.5 Revised fitting curves of salinity and tidal range as well as normalized salinity difference and runoff

圖6 實測鹽度、鹽度差與修正擬合結(jié)果對比Fig.6 Comparison of measured and revised fitting results of salinity and salinity difference

由于對潮差-鹽度擬合關(guān)系進(jìn)行了修正，相應(yīng)的徑流-鹽度差擬合關(guān)系式也需要重新修正。采用與上文同樣的方法，可得到修正歸一化鹽度差SQSD1與徑流Q的擬合關(guān)系，進(jìn)而去歸一化便可得到修正徑流擬合鹽度差SFQ1：

(7)

SFQ1=SQSD1SFT1

(8)

修正歸一化鹽度差SQSD1與徑流Q擬合曲線如圖5(b)所示，兩者的擬合優(yōu)度R2為0.71。實測徑流作用對應(yīng)鹽度差與修正的徑流擬合鹽度差的時間序列對比見圖6。

通過潮差-鹽度修正擬合關(guān)系(式(6))以及徑流-鹽度差的修正關(guān)系(式(7)和(8))，便可得到修正鹽度統(tǒng)計模型：

(9)

式中SrevF為修正鹽度統(tǒng)計模型計算的鹽度，其與實測鹽度對比情況見圖6。

比較圖6與圖4，可以看出修正鹽度統(tǒng)計模型比鹽度統(tǒng)計模型更好地再現(xiàn)了鹽度的實際變化。理論上，修正方法可以反復(fù)進(jìn)行以提高模型精度，但實際應(yīng)用中，第一次修正效果較為明顯，此后修正效果的提高相對有限，建議采用一次修正為宜。

4 誤差分析和討論

本文實測鹽度數(shù)據(jù)序列的均值為0.94‰。鹽度統(tǒng)計模型計算的鹽度序列均值為1.16‰，均值偏差為0.22‰，相對誤差23.4%；修正模型計算的鹽度序列均值為0.96‰，均值偏差約0.02‰，均值相對誤差為2.1%。需要指出的是，雖然均值偏差、相對誤差結(jié)果較好，但對于具體某天的鹽度，尤其是實際鹽度很低的情況，模型誤差比較大。這主要是因為鹽度很低，幾乎為0時，河口不受鹽水入侵影響，導(dǎo)致期間鹽度與潮汐、徑流量的相關(guān)性差，因而統(tǒng)計模型對此期間數(shù)據(jù)識別能力相應(yīng)較差。即便統(tǒng)計模型計算出的鹽度較小，由于實測鹽度幾乎為0，其相對誤差也可能極大，沒有代表性。實際管理應(yīng)用中，往往關(guān)心鹽水入侵較強(qiáng)導(dǎo)致的高鹽度情況。

實際河口鹽水入侵動力機(jī)制非常復(fù)雜，其影響因素除徑流和潮汐外，還受河口地形變化、海平面上升以及風(fēng)力等因素影響[19]。有研究指出，珠江磨刀門口門攔門沙的蝕退、口門向“河-波型”類型轉(zhuǎn)化以及磨刀門水道受采砂形成河床下切均會加劇磨刀門水道鹽水入侵[19]；不同風(fēng)向?qū)χ榻涌邴}水入侵影響不一，北風(fēng)、東北風(fēng)促進(jìn)虎門鹽水的入侵，東風(fēng)則相反[11]。此外，短期特殊氣象過程也會對河口鹽水入侵過程造成影響，如臺風(fēng)“納沙”期間，由于臺風(fēng)增水效應(yīng)使得磨刀門鹽度分布呈現(xiàn)出“雙峰”特征。本文主要關(guān)注對河口影響較大的徑流和潮汐兩個動力因子，基于對徑流和潮汐影響機(jī)理的認(rèn)識和統(tǒng)計學(xué)方法提出了枯水期磨刀門鹽水入侵統(tǒng)計模型構(gòu)建方法。本文提出的鹽水入侵統(tǒng)計模型與建模采用的數(shù)據(jù)樣本緊密相關(guān)，若樣本數(shù)據(jù)代表性不足則可能影響模型精度；考慮到河口地區(qū)水文地貌情況往往容易發(fā)生變化，進(jìn)行模型預(yù)測應(yīng)用時，還應(yīng)注意建模數(shù)據(jù)的時效性。建議在數(shù)據(jù)積累較為豐富的區(qū)域應(yīng)用本文提出的方法時，盡量采用代表性較好的近期數(shù)據(jù)進(jìn)行建模，并采用不同批次的數(shù)據(jù)進(jìn)行應(yīng)用驗證，以確保模型的可靠性。

5 結(jié) 語

基于日均鹽度、日最大潮差、半月平均徑流量等統(tǒng)計數(shù)據(jù)，分別采用指數(shù)函數(shù)和反比例函數(shù)進(jìn)行最小二乘擬合得到潮差-鹽度擬合關(guān)系式和徑流-鹽度差擬合關(guān)系式，并綜合得到基于對潮差和徑流擬合的珠江河口磨刀門鹽度統(tǒng)計模型。針對鹽度統(tǒng)計模型中潮差擬合的徑流影響問題，通過濾去徑流影響后對潮差擬合的方式進(jìn)行了修正，得到修正鹽度統(tǒng)計模型。構(gòu)建的鹽度統(tǒng)計模型和修正鹽度統(tǒng)計模型能反映珠江河口磨刀門枯水期日特征鹽度的變化過程。誤差分析表明，鹽度統(tǒng)計模型和修正鹽度統(tǒng)計模型總體均值偏差分別約為0.22‰和0.02‰，相對誤差分別為23.4%和2.1%，修正鹽度統(tǒng)計模型比鹽度統(tǒng)計模型具有更高的精度。河口地區(qū)鹽度、徑流和潮位資料較容易獲取，基于本文提出的分離短時間尺度潮差影響和長時間尺度徑流影響的鹽度統(tǒng)計擬合方法，可快速建立河口鹽度與潮差、徑流的統(tǒng)計學(xué)模型，并預(yù)測河口鹽水入侵程度，具有較好的實用性。