湖泊藻類動態模型數據同化模式的改進

李 港，陳 誠，何欣霞,2，何夢男，鄧 躍，陳求穩,2

(1.南京水利科學研究院生態環境研究所，江蘇 南京 210029； 2.重慶交通大學河海學院，重慶 400074；3.中建四局第五建筑工程有限公司，廣東 深圳 518000)

湖泊是重要的淡水資源庫，近些年來，由于營養鹽的大量排放以及氣候變化的影響，湖泊富營養化及水華問題日益嚴重[1]。據2018年中國生態環境狀況公報顯示，我國約有1/3的湖庫達到了富營養化狀態[2]，其中太湖、滇池、巢湖等重要水源地富營養化形勢嚴峻，嚴重威脅飲用水安全和人體健康[3]。水華的精確模擬預測是富營養化湖泊風險管控的基礎，對于富營養化湖泊的治理具有重要意義。

富營養化模型是進行湖泊水華模擬預測的重要工具。國內外學者開發了EFDC(environmental fluid dynamics code)、WASP(water quality analysis simulation program)、Delft3D-Eco和三維水生態動力學(three-dimensional hydro-ecological dynamics, 3DHED)模型等一系列湖泊水質模型[4-7]，已成功應用于滇池[8]、Langelmavesi-Roine湖[9]、維多利亞湖[10]等富營養化湖泊的藻類模擬預測。但由于大范圍藍藻生物量的空間分布獲取困難，模型使用站點觀測數據作為輸入數據時，會產生較大的初值誤差，同時，模型運行過程中模擬誤差也會不斷積累，使得模擬精度無法得到有效保證[11]。數據同化(data assimilation, DA)技術通過將觀測數據融入模型計算過程，校正模型模擬結果，同步更新模型參數，降低了輸入數據和模型結構不確定性的影響[12]。該技術包含連續數據同化和順序數據同化兩大類，其中連續順序同化算法利用同化時間窗口內的所有觀測數據和模型狀態值進行最優估計，不斷調整模型初始場，最終將模型軌跡擬合至同化窗口周期內的觀測值上，主要以變分算法為主[13]。相較于變分算法，以集合卡爾曼濾波為代表的順序數據同化算法能夠提供狀態量的均值及其相應的誤差協方差，其預報誤差也可以隨著模型動態變化，同時由于引入了集合的思想，不需要伴隨或線性算子，具有適用于非線性系統、程序設計相對簡單以及易實現并行運算等特點[14]，近年來得到了廣泛的應用。Loos等[15]結合HSFP、EFDC以及WQFS(water quality forecast system)，建立了基于EnKF (ensemble Kalman filter)的Yeongsan河藻類數值模型和數據同化系統。劉卓等[16]采用2009—2011年太湖站點監測數據對Delft3D-BLOOM模型進行EnKF數據同化，提高了Delft3D-BLOOM模型對太湖藍藻的模擬精度。對于大型湖泊而言，數據同化技術雖然能將初始場的使用時間減少到一個觀測周期內，降低了模型誤差的累積，但加入的同化數據是稀疏站點數據插值得到的藍藻生物量數據，忽略了藍藻生物量的空間異質性，存在較大的初值誤差。遙感數據由于具有較好的空間代表性[17]，尤其對大型湖區的藍藻生物量表征能力優于站點觀測數據，所以將其與數據同化技術相結合，以降低初值誤差，有效提高模型模擬與預測精度。但某些時刻遙感數據質量受云霧等氣象條件的影響，觀測誤差較大，此時加入低質量的遙感數據進行數據同化，反而會造成一個觀測周期內的模擬精度降低。

基于此，本文以太湖為研究區域，采用2014—2016年水環境監測數據率定的3DHED模型進行藻類時空模擬，加入藍藻生物量遙感反演數據進行模型的EnKF數據同化，并提出一種數據同化改進策略，通過比較兩種方式模擬結果與實測數據的均方根誤差(RMSE)，進行數據同化時遙感數據的有效篩選，避免因部分時刻遙感數據誤差較大給模型帶來的影響，提升藍藻生物量模擬精度，以期為富營養化湖泊水華的精確模擬預測和風險管控提供技術支撐。

1 研究區概況與數據資料

1.1 研究區概況

太湖是中國第三大淡水湖泊，總面積為 2 338 km2，平均水深約2 m[18]，分為梅梁灣、竺山湖、貢湖、西北湖、湖心區、東部湖區、西南湖和東太湖8個湖區(圖1)。1980年以來，隨著太湖周邊城市化進程的不斷加快，湖泊的不合理利用開發加劇[19]，如工業生產、圍湖造田以及沿岸漁業發展等，導致大量污染物匯入太湖，使得太湖生態系統逐步退化，水體自凈能力下降，水質日趨惡化，藍藻水華頻發[20]。2007年藍藻水華暴發造成了無錫市近百萬人的飲用水危機，引起全社會廣泛關注[21]。2017年，太湖再次暴發大面積藍藻水華，嚴重影響周邊群眾生產生活[22]，太湖富營養化形勢依然嚴峻。

圖1 研究區域及采樣點分布

1.2 數據資料

2 模型與研究方法

2.1 3DHED模型

3DHED模型通過耦合SELFE模型與SALMO模型，考慮了動力條件對藻類分布的影響，可以有效模擬藻類時空變化[7]。其耦合原理為：①通過SELFE模型將湖泊劃分為若干個三角網格，并將SALMO模型應用到其中的每個網格；②通過SELFE模型計算得到全湖所有網格的水動力數據，然后將其提供給SALMO模型；③根據水動力數據與模型水質初始場信息，利用SALMO模型在每個網格上進行水質模擬。3DHED模型參數率定參考郭靜等[24]的前期研究成果，采用2014—2016年的水環境生態監測數據對模型影響較大的20個參數進行率定。選定的參數及其率定結果如表1所示。

2.2 EnKF數據同化

EnKF計算過程包括預測和分析兩個步驟：設有N個集合，在k=0時刻對每個集合進行初始化。式(1)為預測階段，式(2)～(6)為分析階段。

(1)

式中：Xai,k為k時刻第i集合狀態變量分析值；Xfi,k+1為k+1時刻狀態變量預測值；Mk,k+1為模型算子，本文選取3DHED模型；wi,k為模型結構不確定性引起的模擬誤差，服從均值為0、協方差矩陣為Qk的正態分布。

Xai,k+1=Xfi,k+1+Kk+1[Y0k+1-Hk+1(Xfk+1)]

(2)

(3)

表1 模型參數率定結果

(4)

(5)

(6)

同時進行狀態變量和模型參數同化，對模型模擬精度的提高有明顯效果[25]。本文以太湖優勢藻種藍藻生物量作為狀態變量，以敏感性最高的藍藻最佳生長溫度作為模型參數[16]，同時進行狀態變量和模型參數數據同化。采用EnKF算法進行數據同化時，合理的同化方案和合適的參數設置是保證同化效果的前提。EnKF算法利用集合思想解決了實際應用中背景協方差矩陣估計和預報困難的問題，集合數越多，集合均值協方差與真值協方差越接近，但計算負荷會顯著增加；集合數過少，則可能導致模型誤差協方差估計錯誤。本研究將集合數設置為30、50、70、100、200和500進行數值試驗，結果發現集合數為100時可較好地兼顧模擬精度和計算效率。采用數據同化方法將觀測數據與模擬數據相融合，發現融合的優劣性與模擬誤差協方差和觀測誤差協方差緊密相關[26]，所以選擇合適的模擬誤差和觀測誤差對提高同化精度至關重要。本文參考Chen等[27-28]的研究結果，將觀測誤差與模擬誤差分別設置為1%、10%、20%、30%，篩選觀測誤差和模擬誤差的最優組合，進行同化參數設置。

2.3 改進數據同化

同化模型中背景場和觀測數據的誤差均會導致模型的預報誤差，而背景場與觀測誤差均與實測數據相關，因此本研究提出一種改進的數據同化(modified data assimilation, mDA)策略。首先在有遙感數據的時刻，加入遙感數據進行模型的數據同化計算，同時進行模型模擬計算，直到有站點實測數據的下一時刻，以站點實測數據為參考，計算并比較該時刻同化結果和模型模擬結果的RMSE，從而進行同化時刻數據源的篩選。提出的mDA能夠有效剔除低質量的遙感數據，避免遙感數據質量較差所導致的同化結果精度降低的問題，從而提高同化模型的魯棒性。改進數據同化流程如圖2所示。

圖2 改進的數據同化流程

2.4 模型精度評價指標

本文采用RMSE和一致性系數IOA(index of agreement)定量評價藍藻生物量的觀測數據與模擬結果的擬合度。RMSE表示模擬結果與站點實測數據的偏差，其值越小說明模擬結果與站點實測數據越接近。IOA表示模擬結果與站點實測站點數據的吻合程度，描述了模擬結果與站點實測數據在變化趨勢上的一致性，IOA值介于0～1之間，越接近1表示模擬結果與實測數據吻合度越高。

(7)

(8)

3 結果與分析

3.1 數據同化參數確定

表2為16種不同觀測誤差和模擬誤差組合下模擬結果與觀測數據的RMSE值。從表2中可以看出，當觀測誤差設置為1%且模擬誤差設置為10%時，數據同化結果的RMSE最小，因此本文以觀測誤差1%、模擬誤差10%作為EnKF數據同化方法的參數設置。

表2 不同觀測誤差與模擬誤差下的RMSE值

3.2 數據同化結果分析

針對太湖8個湖區各湖區選擇一個代表性站點進行3DHED模擬結果和DA模擬結果的比較，圖3為8個代表性站點藍藻生物量觀測值、3DHED模擬結果以及DA模擬結果的對比。從圖3可以看出，3DHED模擬結果雖能在一定程度上反映藍藻生物量的變化趨勢，但對藍藻生物量峰值的捕捉能力欠佳，而加入遙感數據進行同化后，同化時刻初值誤差有所降低，模型誤差累積周期縮短，使得DA模擬結果的藍藻生物量變化趨勢與峰值的模擬效果均較3DHED模擬結果有顯著提高。8個代表性站點的3DHED模擬結果、DA模擬結果與站點實測數據的RMSE和IOA及其平均值見表3。從表3可以看出，8個代表性站點3DHED模擬結果的RMSE均值為 1.44 mg/L，IOA均值為0.43；DA模擬結果的RMSE均值為1.34 mg/L，IOA均值為0.60。相較于3DHED模擬結果，8個代表性站點IOA平均值提升了39.5%，RMSE平均值降低了6.9%，表明數據同化方法能夠有效提高模型的模擬精度。

圖3 藍藻生物量觀測值、3DHED模擬結果和DA模擬結果對比

3.3 改進數據同化結果分析

圖4為8個湖區代表性站點藍藻生物量觀測值、DA模擬結果以及mDA模擬結果的對比。從圖4可以看出，DA模擬結果和mDA模擬結果均能較為準確地模擬藍藻生物量的變化趨勢，但mDA捕捉到藍藻生物量峰值的次數有一定增加，表明采用DA方法進行數據同化時，部分同化時刻遙感數據存在較大誤差，mDA方法可進行高質量遙感數據的篩選，有效提升模擬精度。mDA模擬結果的RMSE值介于0.79～1.68 mg/L，8個站點RMSE均值為 1.21 mg/L，在DA模擬結果基礎上降低了9.7%，mDA模擬結果的IOA值介于0.53～0.89，均值為0.68，在DA模擬結果基礎上提高了13.3%。

圖4 藍藻生物量觀測值、DA模擬結果和mDA模擬結果對比

圖5為太湖30個采樣點3DHED模擬結果、DA模擬結果和mDA模擬結果的IOA值及RMSE值對比。3種模擬方式下30個站點的IOA均值分別為0.43、0.64和0.71，RMSE均值為1.42 mg/L、1.27 mg/L和1.16 mg/L。各站點IOA值由大到小順序為mDA、DA、3DHED，RMSE值由小到大順序為mDA、DA、3DHED。總體而言，mDA的藍藻生物量模擬精度最高，但少部分站點存在mDA模擬結果的RMSE高于DA模擬結果的現象。

(a) IOA值

表3為各個湖區3DHED模擬結果、DA模擬結果以及mDA模擬結果的IOA值和RMSE值的對比情況。從表3可以看出，DA模擬結果的IOA均高于3DHED模擬結果的IOA，而mDA模擬結果的IOA值在DA模擬結果的基礎上又得到了進一步提升。對于3種模擬結果的RMSE，mDA模擬結果的RMSE值最低，模擬效果最好，整體上太湖各湖區RMSE值表現為西北部湖區大于東南部湖區。

4 討 論

在真實的環境下，與藍藻相關的一些物理參數(如藍藻沉降率和被捕食率等)是難以精確模擬和確定的，并且參數會隨著時間而變化，其過程遠比模型模擬復雜得多[29]。本文采用3DHED模型模擬2014—2016年太湖全湖藍藻生物量的時空變化時，雖然能夠在一定程度上反映太湖藍藻的變化趨勢，但總體模擬效果欠佳，峰值模擬效果有待進一步提高(圖3)。藍藻從生長到形成水華會經歷下沉—休眠—復蘇—上浮4個階段[30]，冬季沉于底泥中的藻類細胞在春季適宜條件下迅速生長，于夏季聚集上浮形成水華[31]。水動力條件也是決定水華暴發的關鍵因子[32]，當藍藻生物量聚集到一定程度后，風浪的擾動使藍藻細胞通過碰撞形成大群體，風浪作用一旦減弱，便上浮形成表面可見水華[33]。結合本文模擬結果，3DHED模擬過程中未能考慮冬季沉底部分的藍藻，同時對于風力作用下藍藻的遷移過程描述不夠清晰，這可能是導致模型對于藍藻生物量峰值捕捉不夠理想的重要原因[7]。采用遙感數據對3DHED進行EnKF數據同化后，由于各個觀測時刻的藍藻生物量得到及時更新，使得DA模式下的藍藻生物量模擬精度得到明顯提升(圖3)，其IOA均值較3DHED模擬結果增加了48.8%，RMSE均值降低了10.4%。

表3 8個湖區3DHED、DA和mDA模擬結果的IOA值和RMSE值的對比

模擬誤差和觀測誤差是數據同化的關鍵因素[16]，對同化性能有重要影響。設置較小的觀測誤差時，觀測數據精度更高，同化結果會更接近觀測值，反之，同化結果會更接近模擬值。一般情況下，觀測數據的誤差更小，更接近于真實值。李淵等[34-35]在進行不同同化試驗時采用了較小的觀測誤差，取得了較為理想的試驗結果。但遙感數據受到云覆蓋等因素的影響，在某些時刻、某些區域的數據誤差較大，在進行數據同化時，由于模型采用的仍是較小的觀測誤差，會導致遙感數據誤差被低估，同化結果趨近于遙感數據，反而會更加偏離真實值。本文采用改進的數據同化模式后，結果顯示mDA的藍藻生物量模擬效果較DA整體上有了一定的提升，對藍藻生物量峰值的捕捉頻次也有所增加(圖4)，全湖IOA均值在DA基礎上增加10.9%，RMSE均值降低了8.6%，說明未改進前遙感數據在某些時刻存在誤差被低估的現象，模型模擬結果趨向于誤差被低估的遙感數據，從而偏離真實值。

太湖藍藻水華總體呈現東南低西北高的趨勢，這與入湖河流污染負荷輸入密切相關。西南部的天目山水系，城鎮較少，人口密度低，水質較好。西北部的南溪水系水流多來自工農業發達的平原和城鎮地區，入湖污染負荷較大[36]。同時太湖夏季盛行東南風，受到太湖風浪的影響，藍藻水華易在西北部聚集和堆積[37]，因此，重度水華主要集中在太湖西北部沿岸(如竺山湖、梅梁灣等)[38]。本文采用mDA方法進行太湖藍藻生物量的模擬時，除了在一定程度上減小遙感數據誤差對同化結果的影響外，還能較好地模擬藍藻生物量的空間分布。結合圖1和表3可以看出，RMSE值總體表現為西北部湖區大于東南部湖區，這是由于藍藻生物量越大的區域，其RMSE值也越大。

湖泊藻類動態模型采用常微分方程描述湖泊生態系統的食物鏈動態和營養物質循環導致的自身結構不確定性[25]，使得模擬誤差會隨著模擬時間的增加而不斷積累。數據同化技術通過實時加入觀測數據，可及時更新初始場數據[39]，將一個初始場的使用時間縮短到一個觀測周期內，有效降低了模擬誤差的積累，此時初始場數據的誤差對模型模擬精度起到決定性作用。遙感數據具有時間上較為連續、監測范圍廣等特點[40]，但若在觀測時刻受到云霧等因素影響，提供的初始場數據誤差可能大于模型模擬得到的初始場數據。本文以兩種模式下提供的初始場數據模擬到下一觀測時刻，比較得到的模擬結果與實測數據的RMSE，再篩選出與真實情況誤差較小的初始場數據進行模型模擬，從而進一步提高模擬精度。近年來，國內學者針對太湖藍藻生物量模擬開展了一系列研究。張成成[41]利用三維水生態模型模擬了藍藻生物量的變化，驗證站點的IOA值介于0.23～0.50。Huang等[28]采用EnKF方法通過更新模型狀態變量和模型參數以提高太湖富營養化模型對浮游植物的模擬精度，其最佳模擬效果下IOA值達到了0.67。本研究mDA模式下的藍藻生物量模擬結果改善效果明顯(圖5)，全湖平均IOA值為0.71，表明本文提出的mDA同化方法具有較高的模擬精度，尤其對于監測站點稀疏的大型湖庫水華的精確模擬預測具有廣闊的應用前景。值得注意的是，存在少部分站點mDA結果精度低于DA的現象(圖5中拖山、洑東等站點)，這主要是由于本研究是通過計算整個湖區36個站點處同化結果與實測數據的RMSE值，來判斷是否加入遙感數據進行同化，少部分站點遙感數據質量較差，但在整個湖區范圍內精度較高，使得部分站點模擬效果不佳。因此后續的研究可以通過劃分不同湖區來進行遙感數據分區的篩選，從而進一步提高數據同化精度。同時，太湖面積廣闊，不同湖區的藍藻水華嚴重程度亦有差別，進行不同湖區模型參數的分區率定對模型精度提高具有一定的效果。

5 結 論

a.采用遙感數據進行藍藻生物量EnKF數據同化后，同化結果的IOA均值較3DHED模擬結果的IOA均值提升了48.8%，RMSE均值降低了10.4%，DA模擬精度較3DHED模型模擬精度有顯著提高。

b.本文提出的mDA方法對藍藻生物量的模擬精度最高，總體精度由大到小為mDA、DA、3DHED，其中mDA的全湖IOA均值達到0.71，RMSE均值降低至1.16 mg/L，不僅能較好反映藍藻生物量時空變化趨勢，還顯著提高了對藍藻生物量的峰值捕捉能力。