劉 卓，李志杰，胡柳明，林育青，陳求穩(wěn),

(1：三峽大學(xué)水利與環(huán)境學(xué)院，宜昌 443002)(2：南京水利科學(xué)研究院生態(tài)環(huán)境研究中心，南京 210029)(3：中國科學(xué)院生態(tài)環(huán)境研究中心，北京100085)

基于集合卡爾曼濾波的湖泊富營養(yǎng)化模型Delft3D-BLOOM數(shù)據(jù)同化

劉 卓1，李志杰3，胡柳明2，林育青2，陳求穩(wěn)2,3

(1：三峽大學(xué)水利與環(huán)境學(xué)院，宜昌 443002)(2：南京水利科學(xué)研究院生態(tài)環(huán)境研究中心，南京 210029)(3：中國科學(xué)院生態(tài)環(huán)境研究中心，北京100085)

富營養(yǎng)化模型是進(jìn)行湖泊水環(huán)境質(zhì)量預(yù)測和管理的重要工具，然而模型客觀存在的誤差一直是應(yīng)用者關(guān)心的重要問題. 數(shù)據(jù)同化作為連接觀測數(shù)據(jù)與數(shù)值模型的重要方法，可以有效提高模型的準(zhǔn)確性. 集合卡爾曼濾波(EnKF)是眾多數(shù)據(jù)同化算法中應(yīng)用最為廣泛的一種，可進(jìn)行非線性系統(tǒng)的數(shù)據(jù)同化，并能有效降低數(shù)據(jù)同化的計(jì)算量. 本研究以太湖作為具體實(shí)例，選擇Delft3D-BLOOM作為富營養(yǎng)化模型，在數(shù)值實(shí)驗(yàn)確定EnKF集合數(shù)為100、觀測誤差方差為1%、模擬誤差方差為10%的基礎(chǔ)上分別進(jìn)行模型狀態(tài)變量同化以及狀態(tài)變量與關(guān)鍵參數(shù)同步同化. 結(jié)果顯示，僅同化狀態(tài)變量時(shí)，模型預(yù)測精度有所增加；同時(shí)同化狀態(tài)變量和關(guān)鍵參數(shù)時(shí)，可顯著提升模型在湖泊水環(huán)境質(zhì)量預(yù)測中的精度. 該研究為應(yīng)用集合卡爾曼濾波以提高復(fù)雜的湖庫富營養(yǎng)化模型模擬精度提供了有效的方法.

集合卡爾曼濾波；富營養(yǎng)化模型；數(shù)據(jù)同化；湖泊；太湖

水是人類社會重要的基礎(chǔ)性自然資源，而湖泊是世界上重要的淡水資源庫之一. 然而，隨著社會經(jīng)濟(jì)的快速發(fā)展，湖泊水環(huán)境已遭受嚴(yán)重的污染，眾多湖泊呈現(xiàn)出嚴(yán)重的富營養(yǎng)化問題，降低了湖泊功能，加劇了水資源短缺. 無論是經(jīng)濟(jì)發(fā)達(dá)還是發(fā)展中國家與地區(qū)，都面臨著不同程度的湖泊富營養(yǎng)化問題[1]. 加拿大的溫尼伯湖，受集水區(qū)域內(nèi)大量農(nóng)業(yè)生產(chǎn)的影響，每年有高達(dá)5900 t的總磷和56300 t的總氮凈輸入該湖泊[2]，導(dǎo)致沉積物中氮、磷濃度迅速增加，并在夏季出現(xiàn)越來越嚴(yán)重的藍(lán)藻水華. 非洲最大的內(nèi)陸漁業(yè)基地維多利亞湖，受沿岸生活污水、工業(yè)廢水排放等影響，在1990s魚類能夠生長繁殖的區(qū)域面積減少了70%，藍(lán)藻水華暴發(fā)時(shí)，水中藍(lán)藻濃度可達(dá)34000 cells/ml，水體透明度下降到1 m左右[3]. 歐洲湖泊同樣面臨嚴(yán)重的富營養(yǎng)化問題：如匈牙利的巴拉頓湖[4]、位于瑞士和德國交界的Constance湖[5]等都出現(xiàn)過嚴(yán)重的富營養(yǎng)化現(xiàn)象. 我國由于近幾十年經(jīng)濟(jì)高速發(fā)展，湖泊資源的不合理開發(fā)利用，使湖泊富營養(yǎng)化進(jìn)程大大加快，如滇池、太湖、巢湖等都已處于富營養(yǎng)化狀態(tài).

湖泊富營養(yǎng)化和藻類水華嚴(yán)重威脅著生態(tài)安全和飲用水安全，在湖泊生態(tài)系統(tǒng)未根本性恢復(fù)到健康狀態(tài)之前，及時(shí)準(zhǔn)確進(jìn)行湖泊水華預(yù)測對湖泊水環(huán)境管理具有重要作用. 富營養(yǎng)化數(shù)學(xué)模型是進(jìn)行湖泊水環(huán)境預(yù)測、治理和管理的重要工具. 從1970s至今，富營養(yǎng)化模型有了很大的發(fā)展，經(jīng)歷了由單一元素向多種元素、局部過程向系統(tǒng)過程、時(shí)空均質(zhì)性向時(shí)空異質(zhì)性的發(fā)展[6]，大致可分為4類：單一營養(yǎng)物質(zhì)負(fù)荷模型[7]、浮游植物與營養(yǎng)鹽相關(guān)模型[8]、生態(tài)水動力模型[9-13]以及人工智能模型[14-18]. 由于對富營養(yǎng)化及水華暴發(fā)的機(jī)理尚未完全清楚，導(dǎo)致模型過程概化、參數(shù)率定、邊界條件設(shè)定等存在誤差，從而降低模型的精度. 隨著觀測技術(shù)的快速發(fā)展，高質(zhì)量的監(jiān)測數(shù)據(jù)量迅速增加，如何有效利用不斷增長的監(jiān)測數(shù)據(jù)和持續(xù)改進(jìn)的數(shù)值模型，是當(dāng)前研究的熱點(diǎn).

數(shù)據(jù)同化是將觀測數(shù)據(jù)與數(shù)值模型計(jì)算結(jié)果相融合，對模型狀態(tài)進(jìn)行最優(yōu)估計(jì)，實(shí)現(xiàn)模型狀態(tài)和參數(shù)的不斷更新，從而提高模型的模擬精度. 國內(nèi)外學(xué)者在數(shù)據(jù)同化方面開展了大量研究，已成功應(yīng)用于數(shù)值天氣預(yù)報(bào)[19]、海洋科學(xué)[20]、水文預(yù)報(bào)[21]以及陸面過程[22]等多個(gè)領(lǐng)域. 1990s后，數(shù)據(jù)同化發(fā)展為連續(xù)數(shù)據(jù)同化和順序數(shù)據(jù)同化兩大類. 連續(xù)數(shù)據(jù)同化算法定義一個(gè)同化的時(shí)間窗口T，利用該同化窗口內(nèi)的所有觀測數(shù)據(jù)和模型狀態(tài)值進(jìn)行最優(yōu)估計(jì)，通過迭代以不斷調(diào)整模型初始場，最終將模型軌跡擬合到同化窗口周期內(nèi)獲取的所有觀測上，如三維變分和四維變分算法等[23]. 三維變分進(jìn)行的是三維空間的全局分析，避免了分析不是全局最優(yōu)的問題；但三維變分無法利用后面時(shí)刻的資料來訂正前面的結(jié)果，同化的解在時(shí)間上不連續(xù). 四維變分是三維變分的推廣，考慮了時(shí)間維，成為目前最常用的變分方法，在很多數(shù)值天氣預(yù)報(bào)業(yè)務(wù)系統(tǒng)(如歐洲中期天氣預(yù)報(bào)中心、美國國家大氣研究中心等)中使用. 四維變分可同化多時(shí)刻的資料，也可在代價(jià)函數(shù)中加上其他的弱約束項(xiàng). 但對于復(fù)雜模式而言，由正模型的離散形式直接得到伴隨方程的離散形式，仍然是一個(gè)比較艱巨的任務(wù)，因此，在計(jì)算方面，四維變分比三維變分代價(jià)大得多. 近年來，以卡爾曼濾波為代表的順序數(shù)據(jù)同化方法得到了廣泛應(yīng)用[24]，其獨(dú)特優(yōu)勢在于：①與傳統(tǒng)最優(yōu)插值算法相比，它使得預(yù)報(bào)誤差隨模型動態(tài)發(fā)展，而最優(yōu)插值中的預(yù)報(bào)誤差是給定的，脫離了模型狀態(tài)的動態(tài)變化；②與變分相比，卡爾曼濾波提供了狀態(tài)量的均值及其相應(yīng)的誤差協(xié)方差，而變分無法提供狀態(tài)量的協(xié)方差矩陣；③卡爾曼濾波無需編寫模型的伴隨模式，更容易實(shí)現(xiàn). 由于卡爾曼濾波需要分析和存儲隨時(shí)間變化的協(xié)方差，因此計(jì)算量大，在狀態(tài)線性變化和高斯分布假設(shè)條件下能獲得最優(yōu)解；但是在數(shù)據(jù)同化中狀態(tài)變量通常都是非線性變化的，且高斯分布條件也很難滿足，從而導(dǎo)致其使用受到極大的制約. 基于集合論和統(tǒng)計(jì)理論的集合卡爾曼濾波(EnKF)[25]具有程序設(shè)計(jì)相對簡單、不需要伴隨或切線性算子、可應(yīng)用于非線性系統(tǒng)、易于實(shí)現(xiàn)并行運(yùn)算等優(yōu)點(diǎn)，彌補(bǔ)了傳統(tǒng)卡爾曼濾波的不足.

本研究選擇成熟的Delft3D-BLOOM作為富營養(yǎng)化模型，以太湖為研究對象，建立了基于EnKF的藍(lán)藻濃度數(shù)值模型與數(shù)據(jù)同化系統(tǒng)，通過分析確定了初始場擾動、觀測場擾動、集合數(shù)等同化系統(tǒng)的關(guān)鍵參數(shù)，提高了太湖藍(lán)藻水華預(yù)測精度，為集合卡爾曼濾波在湖庫富營養(yǎng)化模型數(shù)據(jù)同化中的應(yīng)用提供了參考.

1 材料與方法

1.1 研究區(qū)域及采樣點(diǎn)分布

太湖(30°55′40″～31°32′58″N， 119°52′32″～120°36′10″E)是我國第二大淡水湖，是典型的淺水湖泊，湖面面積2338.11 km2，南北長68.5 km，東西平均寬34 km，最大水深2.6 m，平均水深1.89 m，湖盆形態(tài)呈淺碟形[26]. 太湖自1990s中期大部分區(qū)域已達(dá)到中營養(yǎng)到富營養(yǎng)狀態(tài)，尤其是梅梁灣處于重度富營養(yǎng)水平[27]. 針對太湖水體富營養(yǎng)化現(xiàn)象，國家采取了一些相應(yīng)措施改善水質(zhì)，雖然取得一定成效，但未能從根本上解決太湖富營養(yǎng)化問題[28]. 1980s太湖水體中藻類以綠藻、硅藻、藍(lán)藻為總體優(yōu)勢群，分別占40%、28%和20%. 而1990s以來，太湖藍(lán)藻不僅呈全湖性分布，而且出現(xiàn)時(shí)間長，幾乎全年出現(xiàn)，藍(lán)藻門最高時(shí)約占總量的91.6%[29]. 2003年以來，藍(lán)藻水華開始向湖心擴(kuò)散，嚴(yán)重時(shí)幾乎覆蓋整個(gè)太湖的非水生植被區(qū)[30]. 2007年暴發(fā)了大面積藍(lán)藻水華，無錫貢湖水源地受到嚴(yán)重污染，引發(fā)供水危機(jī). 2008年太湖藻類平均密度為1.488×107cells/L，其中藍(lán)藻密度為1.199×107cells/L[29]. 從2009-2011年太湖水質(zhì)特性(表1)可見，藍(lán)藻門生物量遠(yuǎn)大于其他門類. 因此，本研究選擇太湖為對象，開展藍(lán)藻濃度數(shù)值模型的數(shù)據(jù)同化分析.

1.2 數(shù)據(jù)收集

圖1 太湖監(jiān)測點(diǎn)位分布(1：三號標(biāo)、2：拖山、3：直湖港、4：小灣里、5：竺山湖、6：龍頭、7：大貢山、8：沙墩港、9：漁業(yè)村、10：貢湖、11：東太湖、12：廟港、13：戧港外、14：烏龜山、15：平臺山、16：焦山、17：14號燈標(biāo)、18：湖心南、19：橫山、20：大浦口、21：伏東、22：夾浦、23：新塘、24：小梅口、25：大錢閘、26：西山、27：胥口、28：漫山、29：胥湖)Fig.1 Distribution of monitoring sites in Lake Taihu

1.3 富營養(yǎng)化模型簡介

Delft3D-BLOOM模型是廣泛應(yīng)用于計(jì)算初級生產(chǎn)力、葉綠素a濃度、浮游植物功能群結(jié)構(gòu)的二維水生態(tài)模型，由荷蘭WL | Delft Hydraulics開發(fā)，現(xiàn)已得到廣泛應(yīng)用. 它可以模擬包括藍(lán)藻、綠藻和硅藻在內(nèi)多達(dá)15類浮游植物的生消. Los等[31]對荷蘭North Sea的葉綠素a 濃度進(jìn)行了模擬；Chen等[32]利用Delft3D模擬了淀山湖不同季節(jié)水動力水質(zhì)條件下藻類的生長情況；Kuang等[33]在香港大鵬灣應(yīng)用Delft3D模型研究其水體干濕兩季的水文環(huán)境和水華的運(yùn)移規(guī)律；Lee等[34]在珠江入?？诨贒elft3D建立了三維模型來研究赤潮對流. 為了模擬藻類在不同環(huán)境條件下的適應(yīng)性，通常將藻類分為光能限制型(E-型)、氮限制型(N-型)和磷限制型(P-型)3種類型. 模型主要包括營養(yǎng)物質(zhì)輸移和初級生產(chǎn)力等主要過程. 在計(jì)算過程中與水動力模塊Delft3D-FLOW耦合，由于太湖水動力以風(fēng)生流為主，模擬時(shí)風(fēng)場條件設(shè)置參考許旭峰等[26]的研究成果.

表1 2009-2011年太湖水質(zhì)特性

對于淺水湖泊，水體中物質(zhì)的輸移采用包括源匯項(xiàng)及反應(yīng)項(xiàng)的二維對流-擴(kuò)散方程：

(1)

式中，C為濃度(kg/m3)；Dx和Dy分別為x和y方向上的擴(kuò)散系數(shù)(m2/s)；S為源項(xiàng)；t為時(shí)間(s)；vx和vy分別為x和y方向上的渦流粘度(m2/s)；fR(C,t)為反應(yīng)項(xiàng)，考慮了氮的硝化、反硝化、藻類吸收和礦化過程，磷的吸附、解吸、藻類吸收和礦化過程，溶解氧的耗氧、大氣復(fù)氧、藻類光合作用泌氧和呼吸作用耗氧等過程，具體參考Delft3D-WAQ技術(shù)手冊[35].

BLOOM中初級生產(chǎn)力主要取決于溫度修正的特定生長速率、死亡速率和維持性呼吸速率，其函數(shù)為：

(2)

(3)

(4)

1.4 集合卡爾曼濾波(EnKF)

EnKF的計(jì)算過程包括預(yù)報(bào)和分析兩部分，其主要計(jì)算步驟為：

1)預(yù)報(bào)階段：在預(yù)報(bào)階段，根據(jù)公式(5)來預(yù)測集合中每個(gè)成員在T+ΔT時(shí)刻的狀態(tài)變量.

(5)

(6)

(7)

(8)

2)分析階段:在分析階段，根據(jù)公式(9)對參數(shù)進(jìn)行更新，然后通過公式(12)求得更新后的參數(shù)集合估計(jì)最優(yōu)值.

(9)

(10)

(11)

(12)

狀態(tài)變量根據(jù)公式(13)進(jìn)行更新，通過公式(15)求得更新后狀態(tài)變量集合的平均值作為最優(yōu)的狀態(tài)估計(jì).

(13)

(14)

(15)

1.5 數(shù)據(jù)同化方案設(shè)置

采用EnKF 方法對實(shí)際觀測數(shù)據(jù)進(jìn)行同化之前，需要進(jìn)行數(shù)值實(shí)驗(yàn)，以確定合理的同化方案. 首先恰當(dāng)?shù)膮?shù)設(shè)置有利于EnKF方法獲得更令人滿意的結(jié)果[37]. EnKF是利用各成員在集合平均值上的離散度來統(tǒng)計(jì)出預(yù)報(bào)誤差協(xié)方差和分析誤差協(xié)方差. 要準(zhǔn)確估計(jì)一個(gè)概率密度，需要一定的集合數(shù)，如果集合數(shù)太少，模型的誤差協(xié)方差可能會估計(jì)錯(cuò)誤；如果集合數(shù)過多，則會不可避免地帶來巨大的計(jì)算量，并且當(dāng)集合數(shù)超過一定數(shù)量后，同化結(jié)果將趨于穩(wěn)定，而不再有大幅度的改進(jìn)[38]. 在參考前人研究的基礎(chǔ)上，為有效地平衡同化精度和計(jì)算量，本研究對不同的集合數(shù)取值對同化結(jié)果的影響進(jìn)行了研究.

其次，雖然理論上將觀測數(shù)據(jù)與預(yù)報(bào)模型結(jié)合，通過數(shù)據(jù)同化方法能夠?qū)崿F(xiàn)模型模擬結(jié)果和觀測數(shù)據(jù)的最佳融合，但這種“最佳”與模型的模擬誤差協(xié)方差和觀測誤差協(xié)方差的精度密切相關(guān). 預(yù)測誤差協(xié)方差和觀測誤差協(xié)方差的估計(jì)準(zhǔn)確與否，是決定EnKF 方法同化精度的關(guān)鍵. 然而，實(shí)際應(yīng)用中很難精確得到這些統(tǒng)計(jì)量信息[39]，因此對模擬誤差協(xié)方差和觀測誤差協(xié)方差進(jìn)行優(yōu)化估計(jì)是EnKF 方法應(yīng)用過程中的必要步驟. 不同的應(yīng)用實(shí)際其同化方案不同，數(shù)據(jù)同化方案的設(shè)計(jì)應(yīng)兼顧同化精度和計(jì)算負(fù)荷. 因此，在具體應(yīng)用之前，應(yīng)進(jìn)行數(shù)值實(shí)驗(yàn)以確定合理的同化方案.

1.6 模型性能分析

分別采用均方根誤差(root mean square error,RMSE)和一致性指數(shù)(Index of Agreement,IOA)來定量衡量模型的模擬值與觀測值的一致程度.

(16)

(17)

式中，yo為藍(lán)藻生物量觀測值，ys為藍(lán)藻生物量模型模擬值.

2 結(jié)果

選擇藍(lán)藻生物量作為狀態(tài)變量，選擇敏感性最高的E-型藍(lán)藻最大生長速率溫度系數(shù)作為待同化參數(shù)[40]，以太湖29個(gè)站點(diǎn)每月1次的藍(lán)藻生物量實(shí)測值作為觀測值，對模型的狀態(tài)變量和參數(shù)進(jìn)行同化. 首先將模型計(jì)算1年(2009-01-01-2009-12-31)，使模型達(dá)到穩(wěn)定后，在2010-01-01-2011-12-31時(shí)段開展數(shù)據(jù)同化.

圖2 不同集合數(shù)下的均方根誤差Fig.2 Root mean square error (RMSE) of different ensemble sizes

表3 不同觀測誤差和模擬誤差下的均方根誤差

Tab.3 RMSE under different error variance settings for observation and simulation

模型誤差1%10%20%30%觀測誤差1%1.53830.27010.27071.364010%1.36400.61381.05741.675420%0.68920.53211.10991.649830%1.64981.64981.64982.2180

2.1 數(shù)據(jù)同化方案確定

研究設(shè)置了6種集合數(shù)(30、50、70、100、200和500)，并用RMSE來表征不同集合數(shù)下同化后的狀態(tài)變量與觀測值的一致程度. 從圖2可以看出隨著集合數(shù)的增大，RMSE逐漸減小，當(dāng)集合數(shù)達(dá)到100后，RMSE減小速度大幅降低，同化結(jié)果趨于穩(wěn)定.

本研究將觀測誤差和模擬誤差均分別設(shè)置了1%、10%、20%、30% 4種情況. 表3列出了16種不同的觀測誤差和模擬誤差組合下同化結(jié)果的RMSE. 從表中可以看出，當(dāng)觀測誤差為1%、模擬誤差為10%和20%時(shí)，同化結(jié)果的RMSE均較小，且差別不大.

根據(jù)以上分析，確定適合本研究實(shí)例和富營養(yǎng)化模型的集合卡爾曼濾波數(shù)據(jù)同化設(shè)置：集合數(shù)為100，模擬誤差為10%，觀測誤差為1%.

2.2 數(shù)據(jù)同化分析

2.2.1 只同化狀態(tài)變量 采用直接插值法對藍(lán)藻生物量進(jìn)行更新，即參數(shù)值保持不變，不對觀測值和狀態(tài)變量進(jìn)行擾動，直接以觀測值作為狀態(tài)變量的最優(yōu)估計(jì). 將太湖29個(gè)觀測站藍(lán)藻生物量的觀測值插值到整個(gè)研究區(qū)域，在模型按設(shè)置的時(shí)間步長(dt)進(jìn)行計(jì)算的過程中，在有觀測值的時(shí)刻，將插值后得到的觀測場作為下一時(shí)刻藍(lán)藻生物量的初始場對模型繼續(xù)進(jìn)行計(jì)算. 考慮到監(jiān)測站點(diǎn)眾多，計(jì)算結(jié)果無法一一列出，故從8個(gè)湖區(qū)各取1個(gè)代表性站點(diǎn)進(jìn)行結(jié)果展示. 從圖3可以看出，在每次遇到觀測值時(shí)刻之后，數(shù)據(jù)同化模型結(jié)果比Delft3D-BLOOM直接模擬結(jié)果更為接近觀測值，然后又趨近于Delft3D-BLOOM模擬值. 當(dāng)采用狀態(tài)變量同化模式時(shí)，大部分站點(diǎn)模擬結(jié)果的IOA值較Delft3D-BLOOM 直接模擬結(jié)果的IOA值有所上升(圖4)，說明狀態(tài)變量同化對提高模型的預(yù)測精度有一定效果.

2.2.2 同時(shí)同化狀態(tài)變量和模型參數(shù) 采用太湖29個(gè)觀測站點(diǎn)實(shí)測的藍(lán)藻生物量數(shù)據(jù)，對Delft3D-BLOOM模型的狀態(tài)變量和參數(shù)同時(shí)進(jìn)行同化. 結(jié)果顯示經(jīng)過EnKF方法對模型參數(shù)和狀態(tài)同時(shí)進(jìn)行更新后，模型的預(yù)測精度有了明顯的提高(圖5). 模擬結(jié)果不僅可以很好地反映藍(lán)藻生物量的變化趨勢，而且對藍(lán)藻生物量峰值的模擬結(jié)果有顯著改善. 采用EnKF 方法進(jìn)行數(shù)據(jù)同化后，模擬結(jié)果的IOA值有顯著提高，絕大部分達(dá)到0.9以上(圖6)，說明同時(shí)對狀態(tài)變量和參數(shù)進(jìn)行同化對提高模型的預(yù)測精度有明顯效果.

2.3 兩種同化模式的結(jié)果比較

通過對比Delft3D-BLOOM 模型本身模擬結(jié)果的IOA值、只同化狀態(tài)變量、同時(shí)同化狀態(tài)變量和參數(shù)時(shí)的IOA值大小發(fā)現(xiàn)，采用只同化狀態(tài)變量的模式時(shí)，模擬結(jié)果的IOA值較BLOOM 模型本身模擬結(jié)果的IOA值略有上升；采用同時(shí)同化狀態(tài)變量和參數(shù)的模式時(shí)，模擬結(jié)果的IOA值較BLOOM 模擬結(jié)果有了明顯的提高(圖7). 這說明同時(shí)同化狀態(tài)變量和參數(shù)的同化模式較只同化狀態(tài)變量的同化模式更優(yōu).

圖3 藍(lán)藻生物量觀測值(Obs)、Delft3D-BLOOM 直接模擬值(D3D)以及只同化狀態(tài)變量時(shí)的模擬值(DA1)Fig.3 Observations (Obs), simulation results of Delft3D-BLOOM (D3D) and state variables assimilated model (DA1) of blue-green alga biomass

圖4 Delft3D-BLOOM 直接模擬結(jié)果和只同化狀態(tài)變量時(shí)模擬結(jié)果的IOA值Fig.4 IOA values of Delft3D-BLOOM model and state variables assimilated model

3 討論

本研究通過數(shù)據(jù)同化方法，利用監(jiān)測數(shù)據(jù)，持續(xù)動態(tài)地提升了太湖藍(lán)藻濃度模擬模型計(jì)算精度. Evensen等[41]和Houtekamer等[42]的研究結(jié)果表明當(dāng)集合數(shù)為100左右時(shí)即可獲得令人滿意的結(jié)果. 本文數(shù)值實(shí)驗(yàn)表明當(dāng)集合數(shù)<100時(shí)，模型的RMSE隨集合數(shù)的增大而快速減小，集合數(shù)的大小對RMSE影響較大；當(dāng)集合數(shù)>100時(shí)，集合數(shù)對RMSE的影響不明顯,主要原因是EnKF統(tǒng)計(jì)誤差相關(guān)場存在虛假的相關(guān)，而集合數(shù)的增加會減少虛假的相關(guān). 因此，集合數(shù)設(shè)置為100既能達(dá)到較理想的模擬精度，又能兼顧同化系統(tǒng)的計(jì)算效率.

系統(tǒng)噪聲是數(shù)據(jù)同化的一個(gè)關(guān)鍵要素，對同化性能有重要影響[43]. 當(dāng)設(shè)置的觀測誤差較小時(shí)，認(rèn)為觀測值的精度較高，同化結(jié)果就會更接近觀測值；同理，當(dāng)設(shè)置的模擬誤差較小時(shí)，認(rèn)為模型模擬值的精度較高，則同化結(jié)果會較為接近模擬值. 一般情況下，我們認(rèn)為觀測值的誤差較小，相較模擬值，觀測值更接近真實(shí)值. Rens[43]研究了狀態(tài)變量系統(tǒng)噪聲對模型輸出葉綠素a濃度的影響，發(fā)現(xiàn)當(dāng)模擬誤差為1%、5%、10%、15%、20%和25%時(shí)，葉綠素a濃度的誤差均穩(wěn)定在3%左右. 師春香等[44]進(jìn)行單點(diǎn)陸面土壤濕度同化試驗(yàn)時(shí)，根據(jù)經(jīng)驗(yàn)選取模擬誤差3%，觀測誤差分別取1%和3%，同化后進(jìn)行對比得到當(dāng)觀測誤差取1%時(shí)結(jié)果更優(yōu). 從本研究進(jìn)行的16種不同觀測誤差和模擬誤差組合下同化結(jié)果的RMSE來看，當(dāng)模擬誤差一定時(shí)，RMSE隨觀測誤差的增大而增大；當(dāng)觀測誤差一定時(shí)，RMSE隨模擬誤差的增大而增大的趨勢并不明顯. 所以，觀測誤差對同化精度的影響較模擬誤差大，李淵等[24]在太湖葉綠素a模型分析中也得到了相似的結(jié)論. 對本研究而言，最合適的觀測誤差為1%，模擬誤差為10%.

從圖4和圖6的IOA值可以看出兩種同化模式都不同程度提高了模型的模擬精度，但是在個(gè)別時(shí)刻數(shù)據(jù)同化后的結(jié)果比未同化的結(jié)果差，可能此時(shí)模型誤差被低估，以至于模型狀態(tài)更新的時(shí)候忽略了數(shù)據(jù)同化的預(yù)報(bào)集合. 當(dāng)只同化狀態(tài)變量時(shí)，在有觀測值的時(shí)刻數(shù)據(jù)同化模擬結(jié)果接近觀測值，這是因?yàn)椴捎弥苯硬逯捣▽顟B(tài)變量進(jìn)行同化時(shí)，遇到觀測值的時(shí)刻對初始場重新賦值，使得模擬結(jié)果瞬間重合觀測值；但由于初始場對模型計(jì)算結(jié)果的影響隨著時(shí)間步長的增加逐漸消失，模型的計(jì)算結(jié)果逐漸趨近Delft3D-BLOOM直接模擬的結(jié)果. 當(dāng)同時(shí)同化狀態(tài)變量和參數(shù)時(shí)，模型的預(yù)測精度有了很大的提高，模擬結(jié)果不僅可以很好地模擬藍(lán)藻生物量的變化趨勢，還使藍(lán)藻生物量峰值的模擬效果有了很大的提升. 同時(shí)對狀態(tài)變量和參數(shù)進(jìn)行同化時(shí)，一方面對狀態(tài)變量進(jìn)行同化為下一時(shí)段的模型計(jì)算提供了更精確的初始場；另一方面，根據(jù)站點(diǎn)實(shí)測值同化后得到的參數(shù)更適合當(dāng)前時(shí)段的系統(tǒng)，這兩方面共同改善了模型的模擬精度[45]. 因此，采用

圖5 藍(lán)藻生物量觀測值、Delft3D-BLOOM直接模擬值以及同時(shí)同化狀態(tài)變量和參數(shù)時(shí)的模擬值(DA2)Fig.5 Observations, simulation results of Delft3D-BLOOM, simulation results state variables and parameter assimilated model (DA2) of blue-green alga biomass

圖6 Delft3D-BLOOM直接模擬結(jié)果以及同時(shí)同化狀態(tài)變量和參數(shù)時(shí)模擬結(jié)果的IOA值Fig.6 IOA values of Delft3D-BLOOM model, state variables and parameter assimilated model

圖7 Delft3D-BLOOM直接模擬結(jié)果、只同化狀態(tài)變量時(shí)模擬結(jié)果以及同時(shí)同化狀態(tài)變量和參數(shù)時(shí)模擬結(jié)果的IOA 值Fig.7 IOA values of Delft3D-BLOOM, state variables assimilated model and state variables and parameter assimilated model

監(jiān)測數(shù)據(jù)動態(tài)更新狀態(tài)變量和模型參數(shù)可有效提升復(fù)雜的湖泊富營養(yǎng)化數(shù)值模型的模擬精度.

4 結(jié)論

1)同化方案的設(shè)置對模型效率和計(jì)算成本具有顯著影響. 在本研究中，當(dāng)集合數(shù)為100時(shí)，模型數(shù)據(jù)同化的計(jì)算效率和模擬精度均可得到令人滿意的結(jié)果；當(dāng)模擬誤差設(shè)置為10%、觀測誤差設(shè)置為1%時(shí)，同化后的模擬精度最高；觀測誤差較模擬誤差對系統(tǒng)的影響更大.

2)采用EnKF方法進(jìn)行數(shù)據(jù)同化后，模型的模擬結(jié)果IOA值均得到提高. 只同化狀態(tài)變量時(shí)，大部分站點(diǎn)的IOA值略有上升；同時(shí)同化狀態(tài)變量和參數(shù)時(shí)，絕大部分站點(diǎn)的IOA值有顯著提高，且達(dá)到0.9以上.

3)同化狀態(tài)變量對提高模型預(yù)測精度有一定效果，但是效果維持的時(shí)間很短；同時(shí)同化狀態(tài)變量和參數(shù)比只同化狀態(tài)變量時(shí)的結(jié)果更優(yōu)，其模擬結(jié)果不僅能很好地反映藍(lán)藻生物量的變化趨勢，還顯著提高了藍(lán)藻生物量峰值的模擬精度，有助于提升對藍(lán)藻水華的預(yù)測能力.

[1] Qin Boqiang, Xu Hai, Dong Baili eds. Eutrophication of lakes governance theory and practice. Beijing: Higher Education Press, 2011: 161-162. [秦伯強(qiáng), 許海, 董百麗. 富營養(yǎng)化湖泊治理的理論與實(shí)踐. 北京: 高等教育出版社, 2011: 161-162.]

[2] Board LWS. Reducing nutrient loading to Lake Winnipeg and its watershed: Our collective responsibility and commitment to action. Lake Winnipeg Stewardship Board, 2006: 30.

[3] Ochumba PBO, Kibaara DI. Observations on blue-green algal blooms in the open waters of Lake Victoria, Kenya.AfricanJournalofEcology, 1989, 27(1): 23-34.

[4] Wang Debin, Nie Ming, Geng Kaiyou. Comparative study on pollution control between dianchi lake and balaton lake.JournalofKunmingUniversity, 2013, 35(6): 51-53. [王德斌, 聶銘, 耿開友. 滇池與巴拉頓湖污染治理的比較研究. 昆明學(xué)院學(xué)報(bào), 2013, 35(6): 51-53.]

[5] Straile D, Geller W. Crustacean zooplankton in Lake Constance from 1920 to 1995: Response to eutrophication and re-oligotrophication.AdvancesinLimnology, 1998, 53: 255-274.

[6] Liu Yuanbo, Chen Weimin. Review on simulation study on algal dynamics.JLakeSci, 2000, 12(2): 171-177. DOI: 10.18307/2000.0212. [劉元波, 陳偉民. 湖泊藻類動態(tài)模擬. 湖泊科學(xué), 2000, 12(2): 171-177.]

[7] Vollenweider RA. Input-Output models with special reference to the phosphorus loading concept in limnology.SchweizerischeZeitschriftfürHydrologie, 1975, 37(1): 53-84.

[8] Han Fei, Chen Yongcan, Liu Zhaowei. Advance in the eutrophication models for lakes and reservoirs.AdvancesinWaterScience, 2003, 14(6): 785-791. [韓菲, 陳永燦, 劉昭偉. 湖泊及水庫富營養(yǎng)化模型研究綜述. 水科學(xué)進(jìn)展, 2003, 14(6): 785-791.]

[9] J?rgensen SE, Mejer H, Friis M. Examination of a lake model.EcologicalModelling, 1978, 4(2): 253-278.

[10] Liang Li, Deng Yun, Zheng Meifangetal. Predicting of eutrophication in the LongChuan River based on CE-QUAL-W2 model.ResourcesandEnvironmentintheYangtzeBasin, 2014, (s1): 103-111. [梁俐, 鄧云, 鄭美芳等. 基于CE-QUAL-W2模型的龍川江支庫富營養(yǎng)化預(yù)測. 長江流域資源與環(huán)境, 2014, (s1): 103-111.]

[11] Chen Faxing. Numercial simulation of water quality in Minghu reservoir based on EFDC and water quality prediction.YellowRiver, 2016, 38(1): 85-87. [陳法興. 基于EFDC模型的明湖水庫水質(zhì)數(shù)值模擬及預(yù)測. 人民黃河, 2016, 38(1): 85-87.]

[12] Zhong Xiaohang, Wang Feier, Yu Jieetal. Application of WASP model and gini coefficient in the total amount control of water pollutant: A case study in Tiaoxi estuary of south Taihu Lake.JournalofZhejiangUniversity:ScienceEdition, 2015, 42(2): 181-188. [鐘曉航, 王飛兒, 俞潔等. 基于WASP水質(zhì)模型與基尼系數(shù)的水污染物總量分配——以南太湖苕溪入湖口區(qū)域?yàn)槔? 浙江大學(xué)學(xué)報(bào): 理學(xué)版, 2015, 42(2): 181-188.]

[13] Wei Xingyao, Wang Chao, Wang Peifang. AQUATOX model for study of eutrophication of tributaries.WaterResourcesandPower, 2016, (3): 44-48. [魏星瑤, 王超, 王沛芳. 基于AQUATOX模型的入湖河道富營養(yǎng)化模擬研究. 水電能源科學(xué), 2016, (3): 44-48.]

[14] Whigham PA, Recknagel F. Predicting chlorophyll-a in freshwater lakes by hybridising process-based models and genetic algorithms.EcologicalModelling, 2001, 146(1): 243-251.

[15] Zhang Chengcheng, Shen Aichun, Zhang Xiaoqingetal. Application of support vector machine to evaluate the eutrophication status of Taihu Lake.ActaEcologicaSinica, 2013, 33(23): 7563-7569. [張成成, 沈愛春, 張曉晴等. 應(yīng)用支持向量機(jī)評價(jià)太湖富營養(yǎng)化狀態(tài). 生態(tài)學(xué)報(bào), 2013, 33(23): 7563-7569.]

[16] Wang Xuelian, Song Yuzhi, Kong Fanfanetal. Estimation of chlorophyll-a concentration in Taihu Lake by using Back Propagation (BP) neural network forecast model.ChineseJournalofAgrometeorology, 2016, 37(4): 408-414. [王雪蓮, 宋玉芝, 孔繁璠等. 利用BP神經(jīng)網(wǎng)絡(luò)模型對太湖水體葉綠素a含量的估算. 中國農(nóng)業(yè)氣象, 2016, 37(4): 408-414.]

[17] Xiang Xianquan, Tao Jianhua. Eutrophication model of Bohai bay based on GA-SVM.JournalofTianjinUniversity, 2011, 44(3): 215-220. [向先全, 陶建華. 基于GA-SVM的渤海灣富營養(yǎng)化模型. 天津大學(xué)學(xué)報(bào), 2011, 44(3): 215-220.]

[18] Xu Xiaoyi, Han Lianghua, Luo Guyuanetal. Prediction of chlorophyll-a by Genetic BP neural network.ChinaWaterandWastewater, 2012, 28(1): 68-71. [許曉毅, 韓亮華, 羅固源等. 遺傳-BP神經(jīng)網(wǎng)絡(luò)法預(yù)測葉綠素a濃度變化. 中國給水排水, 2012, 28(1): 68-71.]

[19] Panofsky R. Objective weather-map analysis.JournalofMeteorology, 1949, 6(6): 386-392.

[20] Pham DT, Verron J, Roubaud MC. A singular evolutive extended Kalman filter for data assimilation in oceanography.JournalofMarineSystems, 1998, 16(3): 323-340.

[21] Weerts AH, El Serafy GY. Particle filtering and ensemble Kalman filtering for state updating with hydrological conceptual rainfall-runoff models.WaterResourcesResearch, 2006, 42(9): 123-154.

[22] Galantowicz JF, Entekhabi D, Njoku EG. Tests of sequential data assimilation for retrieving profile soil moisture and temperature from observed L-band radiobrightness.IEEETransactionsonGeoscienceandRemoteSensing, 1999, 37(4): 1860-1870.

[23] Ma Jianwen, Qing Sixian. Recent advances and development of data assimilation algorithms.AdvancesinEarthScience, 2012, 27(7): 747-757. [馬建文, 秦思嫻. 數(shù)據(jù)同化算法研究現(xiàn)狀綜述. 地球科學(xué)進(jìn)展, 2012, 27(7): 747-757.]

[24] Li Yuan, Li Yunmei, Lü Hengetal. Sensitivity analysis on Lake Taihu data assimilation scheme of chlorophyll-a concentration.JLakeSci, 2015, 27(1): 163-174. DOI: 10.18307/2015.0119. [李淵, 李云梅, 呂恒等. 太湖葉綠素 a 同化系統(tǒng)敏感性分析. 湖泊科學(xué), 2015, 27(1): 163-174.]

[25] Evensen G. Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics.JournalofGeophysicalResearchAtmospheres, 1994, 99(C5): 10143-10162.

[26] Xu Xufeng, Liu Qingquan. Numerical study on the characteristics of wind-induced current in Taihu lake.ChineseJournalofHydrodynamics, 2009, 24(4): 512-518. [許旭峰, 劉青泉. 太湖風(fēng)生流特征的數(shù)值模擬研究. 水動力學(xué)研究與進(jìn)展, 2009, 24(4): 512-518.]

[27] Guo Jing, Chen Qiuwen, Li Weifeng. Application of lake model SALMO to the Meiliang Bay of Taihu lake.ActaScientiaeCircumstantiae, 2012, 32(12): 3119-3127. [郭靜, 陳求穩(wěn), 李偉峰. 湖泊水質(zhì)模型SALMO在太湖梅梁灣的應(yīng)用. 環(huán)境科學(xué)學(xué)報(bào), 2012, 32(12): 3119-3127.]

[28] Qin Boqiang, Wang Xiaodong, Tang Xiangmingetal. Drinking water crisis caused by eutrophication and cyanobacterial bloom in lake Taihu: Cause and measurement.AdvancesinEarthScience, 2007, 22(9): 896-906. [秦伯強(qiáng), 王小冬, 湯祥明等. 太湖富營養(yǎng)化與藍(lán)藻水華引起的飲用水危機(jī)——原因與對策. 地球科學(xué)進(jìn)展, 2007, 22(9): 896-906.]

[29] Gu Suli, Chen Fang, Sun Jiangling. Analysis of cyanobacteria monitoring and algal blooms in Taihu Lake.WaterResourcesProtection, 2011, 27(3): 28-32. [顧蘇莉, 陳方, 孫將陵. 太湖藍(lán)藻監(jiān)測及暴發(fā)情況分析. 水資源保護(hù), 2011, 27(3): 28-32.]

[30] Wang Changyou, Yu Yang, Sun Yunkunetal. The discussion of the early forecasting of cyanobacteria bloom in the Lake Taihu based on ELCOM-CAEDYM model.ChinaEnvironmentalScience, 2013, 33(3): 491-502. [王長友, 于洋, 孫運(yùn)坤等. 基于ELCOM-CAEDYM模型的太湖藍(lán)藻水華早期預(yù)測探討. 中國環(huán)境科學(xué), 2013, 33(3): 491-502.]

[31] Los FJ, Villars MT, Tol MWMVD. A 3-dimensional primary production model (BLOOM/GEM) and its applications to the (southern) North Sea (coupled physical-chemical-ecological model).JournalofMarineSystems, 2008, 74(1/2): 259-294.

[32] Chen Y, Lin W, Zhu Jetal. Numerical simulation of an algal bloom in Dianshan Lake.ChineseJournalofOceanologyandLimnology, 2016, 34(1): 231-244.

[33] Kuang C, Lee JHW. Physical hydrography and algal bloom transport in Hong Kong waters.ChinaOceanEngineering, 2005, 19(4): 539-556.

[34] Lee JHW, Qu B. Hydrodynamic tracking of the massive spring 1998 red tide in Hong Kong.JournalofEnvironmentalEngineering, 2004, 130(5): 535-550.

[35] WL | Delft hydraulics. Delft3D-WAQ technical reference manual. Delft, The Netherlands, 2009. http: ∥www.delftsoftware.com.

[36] Li Z, Chen Q, Xu Q. Modeling algae dynamics in Meiliang Bay of Taihu Lake and parameter sensitivity analysis.JournalofHydro-environmentResearch, 2014, 9(2): 216-225.

[37] Eknes M, Evensen G. An Ensemble Kalman filter with a 1-D marine ecosystem model.JournalofMarineSystems, 2002, 36(1): 75-100.

[38] Huang Chunlin, Li Xin. Sensitivity analysis on land data assimilation scheme of soil moisture.AdvancesinWaterScience, 2006, 17(4): 457-465. [黃春林, 李新. 土壤水分同化系統(tǒng)的敏感性試驗(yàn)研究. 水科學(xué)進(jìn)展, 2006, 17(4): 457-465.]

[39] Simon E, Bertino L. Application of the Gaussian anamorphosis to assimilation in a 3-D coupled physical-ecosystem model of the North Atlantic with the EnKF: a twin experiment.OceanScience, 2009, 5(4): 495-510.

[40] Huang J, Gao J, H?rmann G. Hydrodynamic-phytoplankton model for short-term forecasts of phytoplankton in Lake Taihu, China.Limnologica-EcologyandManagementofInlandWaters, 2012, 42(1): 7-18.

[41] Evensen G, Van LPJ. Assimilation of Geosat altimeter data for the Agulhas current using the ensemble Kalman filter with a quasigeostrophic model.MonthlyWeatherReview, 1996, 124(1): 85-96.

[42] Houtekamer PL, Mitchell HL. A sequential ensemble Kalman filter for atmospheric data assimilation.MonthlyWeatherReview, 2001, 129(1): 123-137.

[43] Rens EG. Improving the algal bloom prediction in the North Sea by Ensemble Kalman Filtering in the GEM/BLOOM Model[Dissertation]. Delft: Delft University of Technology, 2013.

[44] Shi Chunxiang, Xie Zhenghui, Qian Huietal. Based on satellite remote sensing data of China’s regional soil moisture EnKF data assimilation.ScienceChinaEarthScience, 2011, 41(3): 375-385. [師春香, 謝正輝, 錢輝等. 基于衛(wèi)星遙感資料的中國區(qū)域土壤濕度 EnKF 數(shù)據(jù)同化. 中國科學(xué): 地球科學(xué), 2011, 41(3): 375-385.]

[45] Mao J, Lee JH, Choi K. The extended Kalman filter for forecast of algal bloom dynamics.WaterResearch, 2009, 43(17): 4214-4224.

Ensemble Kalman filter based data assimilation in the Delft3D-BLOOM lake eutrophication model

LIU Zhuo1, LI Zhijie3, HU Liuming2, LIN Yuqing2& CHEN Qiuwen2,3**

(1:CollegeofHydraulicandEnvironmentalEngineering,ThreeGorgesUniversity,Yichang443002,P.R.China)(2:CenterforEco-EnvironmentalResearch,NanjingHydraulicResearchInstitute,Nanjing210029,P.R.China)(3:ResearchCenterforEco-EnvironmentSciences,ChineseAcademyofSciences,Beijing100085,P.R.China)

Numerical eutrophication model is an important tool to predict and manage the ecosystem of lakes and reservoirs. However, the objective errors of the model are always vital problems the users concerned. Data assimilation, which connects observations and model simulations, can effectively improve the accuracy of models. Ensemble Kalman filter (EnKF), which is one of the most widely used methods for data assimilation, is suitable for nonlinear system and has high computation efficiency. In this research, the Delft3D-BLOOM was taken as the eutrophication model, and Lake Taihu was taken as the study case. After numerical testing, the ensemble size was set to 100, the observation error variance was set to 1%, and the simulation error variance was set to 10%. Two data assimilation modes, assimilation of model state variables and synchronous assimilation of both state variables and key parameters, were examined. The results showed that the fitness between model simulation and observation was slightly improved when the state variable was updated. When both the state variables and parameters were assimilated, the fitness was significantly improved. The study provides a promising approach in using EnKF to improve the simulation accuracy of complex eutrophication models.

Ensemble Kalman filter; eutrophication model; data assimilation; lake; Lake Taihu

國家自然科學(xué)基金項(xiàng)目(51579149，51609142)和江蘇省水利科技項(xiàng)目(2016021)聯(lián)合資助. 2016-09-21收稿；2016-12-20收修改稿. 劉卓(1991～)，男，碩士研究生；E-mail: rrzliu@foxmail.com.

；E-mail: qchen@rcees.ac.cn.

DOI 10.18307/2017.0505