珠江河口磨刀門枯水期鹽度統計模型

王 彪，朱建榮，呂愛琴

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

珠江是我國第二大河，磨刀門為珠江河口八大入?？陂T之一，是主要的徑流和泥沙出口，也是珠海和澳門主要水源地。20世紀90年代以來，珠江三角洲地區咸潮活動越來越頻繁，持續時間增加，上溯影響范圍越來越大，強度趨于嚴重，發展態勢不容樂觀[1]，嚴重影響了珠海和澳門的供水安全。

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

1 研究數據和方法

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

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

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

(a) 時段1流量

(b) 時段1潮位潮差

(c) 時段1鹽度

(d) 時段2流量

(e) 時段1潮位、潮差

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

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

2 模型構建

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

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

(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)中的曲線即為擬合曲線，潮差-鹽度的擬合優度R2為0.90?；诖顺辈?鹽度擬合關系式，可計算出潮差隨時間變化對應的潮差擬合鹽度變化，其與實測鹽度對比見圖4。

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

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

SSD=SD/SFT0

(2)

其中

SD=S-SFT0

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

通過分析歸一化實測鹽度差與徑流之間的關系，發現二者具有反比例函數關系。因此，采用反比例函數進行歸一化實測鹽度差與徑流擬合，得到徑流-鹽度差擬合函數：

(3)

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

SFQ0=SQSD0SFT0

(4)

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

(5)

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

3 修正鹽度統計模型

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

SFT1=0.001 127e4.372HTR

(6)

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

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

(a) 鹽度與潮差

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

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

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

(7)

SFQ1=SQSD1SFT1

(8)

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

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

(9)

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

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

4 誤差分析和討論

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

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

5 結 語

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