張 帥，彭福利，季雨來(lái)，張 京，張奇謀，李 琪，錢(qián) 瑞，齊凌艷，黃佳聰**

(1：安徽師范大學(xué)地理與旅游學(xué)院，蕪湖 241003) (2：中國(guó)科學(xué)院南京地理與湖泊研究所, 中國(guó)科學(xué)院流域地理學(xué)重點(diǎn)實(shí)驗(yàn)室，南京 210008) (3：中國(guó)環(huán)境監(jiān)測(cè)總站，北京 100012) (4：資源環(huán)境與地理信息工程安徽省工程技術(shù)研究中心，蕪湖 241003)

湖庫(kù)富營(yíng)養(yǎng)化與藻類(lèi)水華是我國(guó)面臨的突出水環(huán)境問(wèn)題，2020年《中國(guó)環(huán)境狀況公報(bào)》數(shù)據(jù)顯示，在我國(guó)110個(gè)重點(diǎn)監(jiān)測(cè)的湖庫(kù)中有29%處于富營(yíng)養(yǎng)化狀態(tài)，其中三湖(太湖、滇池、巢湖)全部處于富營(yíng)養(yǎng)化狀態(tài)[1]，造成藍(lán)藻水華頻發(fā)[2]；2018年夏季，巢湖湖心區(qū)和西北岸暴發(fā)了面積約121.38 km2的藍(lán)藻水華，嚴(yán)重影響周邊居民日常生活和用水安全[3]. 提前預(yù)測(cè)藻類(lèi)未來(lái)變化，并采取應(yīng)急措施能夠有效緩解藍(lán)藻水華造成的危害[4]；近年來(lái)，我國(guó)在巢湖等210個(gè)重點(diǎn)湖庫(kù)設(shè)置了349個(gè)監(jiān)測(cè)站，積累了長(zhǎng)時(shí)間歷史系列的人工與自動(dòng)監(jiān)測(cè)數(shù)據(jù)，并在此基礎(chǔ)上，構(gòu)建了包含白洋淀[5]、太湖[6]、巢湖[7]等湖泊的藻類(lèi)動(dòng)態(tài)模型. 但湖泊藻類(lèi)動(dòng)態(tài)變化過(guò)程十分復(fù)雜[8]，數(shù)值模擬難度大，亟需依托與日俱增的海量監(jiān)測(cè)數(shù)據(jù)，發(fā)展模型-數(shù)據(jù)融合算法，提升模型模擬精度[9].

數(shù)據(jù)同化是在模型和觀測(cè)的誤差估計(jì)基礎(chǔ)上，在模型運(yùn)行過(guò)程中，實(shí)時(shí)融合觀測(cè)數(shù)據(jù)的方法[10]；是融合模型與數(shù)據(jù)、提升模擬精度的重要技術(shù)之一[11]. 最初應(yīng)用于數(shù)值天氣預(yù)報(bào)，為天氣預(yù)報(bào)提供初始場(chǎng)的數(shù)據(jù)處理技術(shù)[12]；目前已廣泛應(yīng)用于大氣[13]、海洋[14]等領(lǐng)域，在生態(tài)環(huán)境領(lǐng)域也有相關(guān)應(yīng)用[15]，包括流域水文模型、湖泊藻類(lèi)模型等，并取得了良好應(yīng)用效果[16]. 集合卡爾曼濾波(Ensemble Kalman Filter, EnKF)是數(shù)據(jù)同化領(lǐng)域主流的算法之一，其通過(guò)蒙特卡洛方法計(jì)算狀態(tài)的預(yù)報(bào)誤差協(xié)方差，用集合的思想解決了實(shí)際應(yīng)用中背景誤差協(xié)方差矩陣的估計(jì)和預(yù)報(bào)困難的問(wèn)題，能夠有效降低數(shù)據(jù)同化計(jì)算量，具有程序設(shè)計(jì)相對(duì)簡(jiǎn)單、可應(yīng)用于非線(xiàn)性系統(tǒng)、易于實(shí)現(xiàn)并行運(yùn)算等優(yōu)點(diǎn)[11,17]. 但集合卡爾曼濾波應(yīng)用于湖泊藻類(lèi)動(dòng)態(tài)模型中，缺乏模型敏感參數(shù)動(dòng)態(tài)變化的現(xiàn)實(shí)[18-20]，導(dǎo)致模型模擬精度提升效果有限，還存在可能增加模型不確定性的不足[21]. 參數(shù)敏感性分析能夠有效識(shí)別復(fù)雜模型的關(guān)鍵參數(shù)，是實(shí)時(shí)跟蹤敏感參數(shù)動(dòng)態(tài)變化的有效方法[22]，有望與數(shù)據(jù)同化耦合，實(shí)現(xiàn)對(duì)模型敏感參數(shù)的靶向更新，精準(zhǔn)提升模型模擬精度.

本研究通過(guò)集合卡爾曼濾波與參數(shù)敏感性分析的耦合，開(kāi)發(fā)了一種能夠?qū)崟r(shí)識(shí)別敏感參數(shù)的新型數(shù)據(jù)同化算法，為湖泊藻類(lèi)動(dòng)態(tài)模型模擬精度提升提供技術(shù)支持. 新型數(shù)據(jù)同化算法將參數(shù)敏感性分析耦合到集合卡爾曼濾波，在模型運(yùn)行每一階段，根據(jù)敏感性分析定位敏感參數(shù)，采用集合卡爾曼濾波同步同化狀態(tài)變量和敏感參數(shù). 由此可精確定位并調(diào)整每階段敏感參數(shù)，實(shí)時(shí)校正模型運(yùn)行軌跡，減少模型模擬不確定性，提升模型模擬精度. 為驗(yàn)證新算法，以巢湖為研究對(duì)象，依據(jù)水質(zhì)自動(dòng)監(jiān)測(cè)數(shù)據(jù)，將新型數(shù)據(jù)同化算法應(yīng)用于湖泊藻類(lèi)動(dòng)態(tài)模型[23]，檢驗(yàn)新算法對(duì)模型精度的提升效果，研究成果也為其它湖庫(kù)水質(zhì)模型數(shù)據(jù)同化提供技術(shù)參考.

1 材料與方法

1.1 研究區(qū)概況

巢湖(31°26′～31°42′N(xiāo)，117°17′～117°51′E)位于安徽省中部、合肥市內(nèi)，由巢湖、肥東、肥西和廬江一市三縣環(huán)抱，湖長(zhǎng)54.5 km，最大寬度21 km，平均水深3 m，面積約787 km2，容積約20億m3(圖1)，是沿岸居民生活生產(chǎn)的主要水源[24]. 巢湖入湖河流約33條，其中6條主要入湖河流為杭埠河、白石天河、兆河、柘皋河、南淝河和派河；出湖河流主要為裕溪河，連接長(zhǎng)江. 近年來(lái)隨著流域內(nèi)人口增加，工農(nóng)業(yè)迅速發(fā)展，城鎮(zhèn)工業(yè)廢水、生活污水排入，導(dǎo)致水體營(yíng)養(yǎng)鹽增加，藍(lán)藻水華事件頻發(fā)，水質(zhì)惡化，影響周?chē)用耧嬘盟踩玔25]. 此外，藍(lán)藻水華還分泌毒素，危害人類(lèi)健康[26]. 因此，巢湖富營(yíng)養(yǎng)化已成為突出的生態(tài)、環(huán)境和社會(huì)問(wèn)題[27].

圖1 巢湖與水質(zhì)自動(dòng)監(jiān)測(cè)站分布Fig.1 Locations of Lake Chaohu and its automatic water quality monitoring stations

1.2 數(shù)據(jù)

巢湖水質(zhì)監(jiān)測(cè)數(shù)據(jù)來(lái)自8個(gè)自動(dòng)監(jiān)測(cè)站的逐日監(jiān)測(cè)(圖1)，監(jiān)測(cè)指標(biāo)包括水溫(WT, ℃)、總磷(TP, mg/L)、總氮(TN, mg/L)及葉綠素a濃度(Chl.a, μg/L)等，監(jiān)測(cè)時(shí)間為2020年1-5月. 氣象數(shù)據(jù)來(lái)自國(guó)家氣象信息中心(http://data.cma.cn/)的巢湖站(圖1)，位于巢湖東部沿岸，代表巢湖氣象條件；氣象指標(biāo)包括最低氣溫(℃)、最高氣溫(℃)、平均氣溫(℃)、日照時(shí)長(zhǎng)(h)、平均風(fēng)速(m/s)及降水量(mm)等. 以上數(shù)據(jù)為湖泊藻類(lèi)動(dòng)態(tài)模型運(yùn)行與校驗(yàn)數(shù)據(jù).

1.3 巢湖藻類(lèi)動(dòng)態(tài)模型

目前已開(kāi)發(fā)了系列湖泊藻類(lèi)模擬模型，如澳大利亞西澳大學(xué)水研究中心的ELCOM-CAEDYM模型[28]、威廉瑪麗大學(xué)海洋學(xué)院維吉尼亞海洋科學(xué)研究所的Environmental Fluid Dynamics Code(EFDC)模型等[29]. 但此類(lèi)模型多需要海量、高時(shí)空分辨率的數(shù)據(jù)輸入，部分?jǐn)?shù)據(jù)難以實(shí)時(shí)獲取，不適用于湖泊藻類(lèi)的短期(1～3 d)模擬預(yù)測(cè)[30]. 為解決模型復(fù)雜數(shù)據(jù)輸入問(wèn)題，國(guó)內(nèi)外結(jié)合研究區(qū)特征，開(kāi)發(fā)了少量數(shù)據(jù)輸入的湖泊藻類(lèi)動(dòng)態(tài)模型[31-32]. 本研究采用黃佳聰?shù)葮?gòu)建的湖泊藻類(lèi)動(dòng)態(tài)模型[23](圖2).

圖2 湖泊藻類(lèi)動(dòng)態(tài)模型概念圖[23]Fig.2 Conceptual diagram of phytoplankton dynamic model for lakes[23]

該模型采用葉綠素a濃度表征藻類(lèi)生物量，藻類(lèi)動(dòng)態(tài)變化模擬分為兩個(gè)模塊(表1、表2)：①藻類(lèi)生消模塊：用于模擬藻類(lèi)在太陽(yáng)輻射、水溫、營(yíng)養(yǎng)鹽等條件下的生長(zhǎng)，以及在排泄、沉降、呼吸、死亡、被浮游動(dòng)物捕食等機(jī)理過(guò)程中的消亡；②基于二維水動(dòng)力的藻類(lèi)運(yùn)移模塊：用于模擬藻類(lèi)在水動(dòng)力作用下的運(yùn)移過(guò)程. 其它機(jī)理過(guò)程的數(shù)學(xué)表達(dá)可以參考文獻(xiàn)[23]. 為了兼顧模型運(yùn)行效率與算法耦合效率，本研究?jī)H考慮了點(diǎn)位的藻類(lèi)模擬，即模型為零維模型，僅考慮藻類(lèi)生長(zhǎng)與消亡等過(guò)程導(dǎo)致的藻類(lèi)動(dòng)態(tài)變化，不考慮水動(dòng)力輸移對(duì)藻類(lèi)的影響. 相比于其它模型，該模型輸入數(shù)據(jù)易實(shí)時(shí)獲取，能滿(mǎn)足國(guó)內(nèi)湖泊藻類(lèi)模擬預(yù)測(cè)的需求.

表1 湖泊藻類(lèi)生消過(guò)程表達(dá)[23]

表2 藻類(lèi)模型參數(shù)與變量[23]

Tab.2 Parameters and variables of the phytoplankton model[23]

類(lèi)型簡(jiǎn)稱(chēng)描述單位初值參數(shù)Umax藻類(lèi)最大生長(zhǎng)速率d-11.145Topt藻類(lèi)生長(zhǎng)的最佳水溫℃27α水及非藻類(lèi)吸收短波輻射的平均消光系數(shù)m-10.45β藻類(lèi)吸收短波輻射的平均消光系數(shù)μg Chl/(L·m)0.016h水深m3.0Iopt飽和光強(qiáng)MJ/(m2·d)12Gsc太陽(yáng)常數(shù)MJ/(m2·d)1367φ緯度(°)31.5KP影響藻類(lèi)磷吸收的米氏常數(shù)μg/L10KN影響藻類(lèi)氮吸收的米氏常數(shù)μg/L22K藻類(lèi)沉降速率m/d0.0864GRmax浮游動(dòng)物最大捕食速率d-10.09Fmin可供捕食藻類(lèi)最小濃度μg/L100Fs可供捕食藻類(lèi)米氏常數(shù)μg/L500km藻類(lèi)死亡引起葉綠素a減少速率系數(shù)d-10.027kr藻類(lèi)呼吸引起葉綠素a減少速率系數(shù)d-10.17ke藻類(lèi)排泄引起葉綠素a減少速率系數(shù)d-10.01?溫度因子1.08待解變量T水體溫度℃r光衰減系數(shù)mChl葉綠素a濃度μg/LI水體表面光強(qiáng)MJ/(m2·d)S日照時(shí)數(shù)hS0晝長(zhǎng)hk離心校正因子U藻類(lèi)生長(zhǎng)速率d-1RA藻類(lèi)呼吸速率d-1SA藻類(lèi)沉降速率d-1MA藻類(lèi)死亡速率d-1GA藻類(lèi)被捕食速率d-1EA藻類(lèi)排泄速率d-1δ太陽(yáng)赤緯°ωs日落時(shí)角°DP水體可溶性磷μg/LDN水體可溶性氮μg/LRAmax藻類(lèi)最大呼吸速率d-1F可供捕食藻類(lèi)濃度μg/L

1.4 新型數(shù)據(jù)同化算法

基于集合卡爾曼濾波的數(shù)據(jù)同化在模型運(yùn)行的每一階段會(huì)對(duì)狀態(tài)變量和模型參數(shù)進(jìn)行調(diào)整[18]. 算法運(yùn)行分兩步驟：預(yù)測(cè)和更新. 預(yù)測(cè)：利用模型模擬狀態(tài)變量先驗(yàn)值集合；更新：根據(jù)實(shí)測(cè)值統(tǒng)計(jì)卡爾曼增益并估計(jì)狀態(tài)變量和參數(shù)，作為下一時(shí)刻的模型輸入[33]. 算法遞推式運(yùn)行，每階段的更新只依賴(lài)于上一時(shí)刻集合空間，運(yùn)行效率高[11]，但難以應(yīng)對(duì)模型參數(shù)敏感性動(dòng)態(tài)變化的問(wèn)題.

敏感性分析能夠在多參數(shù)模型中有效識(shí)別關(guān)鍵參數(shù)，從而減少敏感參數(shù)對(duì)模型的影響，使模型模擬更加精確[22]. 本研究決定采用的擾動(dòng)分析法簡(jiǎn)單實(shí)用且應(yīng)用普遍[34]，采用此方法對(duì)湖泊藻類(lèi)動(dòng)態(tài)模型的參數(shù)進(jìn)行敏感性分析. 擾動(dòng)分析法基本思路如下：選定一個(gè)參數(shù)作為敏感性分析目標(biāo)，其它參數(shù)都保持不變，將目標(biāo)參數(shù)在一定范圍內(nèi)取值擾動(dòng)n次并計(jì)算模型輸出，依次求出相鄰兩次擾動(dòng)的敏感程度，取平均值作為目標(biāo)參數(shù)的相對(duì)敏感度S(式(6))(表3). 該方法簡(jiǎn)單高效，既方便耦合模型又能兼顧運(yùn)行效率. 根據(jù)本研究湖泊藻類(lèi)動(dòng)態(tài)模型的機(jī)理和結(jié)構(gòu)，將模型的參數(shù)敏感程度劃分為5個(gè)等級(jí)：不敏感(Ⅰ，|S|<0.25)；弱敏感(Ⅱ，0.25≤|S|<0.5)；一般敏感(Ⅲ，0.5≤|S|<0.75)；比較敏感(Ⅳ，0.75≤|S|<1.0)；極度敏感(Ⅴ，|S|≥1.0). 將敏感程度|S|>0.5作為敏感參數(shù)用于調(diào)整，敏感度閾值0.5是根據(jù)模型的機(jī)理結(jié)構(gòu)及模擬效果綜合率定的，以此為閾值調(diào)整參數(shù)能夠獲得的最優(yōu)模擬結(jié)果.

本研究開(kāi)發(fā)的新型數(shù)據(jù)同化算法將參數(shù)敏感性分析耦合到集合卡爾曼濾波過(guò)程中，使其能夠精確定位并調(diào)整每階段模型敏感參數(shù)，實(shí)時(shí)有效校正模型運(yùn)行軌跡，提升模型模擬精度. 綜合考慮集合卡爾曼濾波對(duì)參數(shù)的調(diào)整過(guò)程以及每一階段參數(shù)敏感性分析的時(shí)效性，本研究提出的新型數(shù)據(jù)同化算法運(yùn)行分3個(gè)步驟：預(yù)測(cè)、分析與更新(表3).

①預(yù)測(cè)：與集合卡爾曼濾波數(shù)據(jù)同化的預(yù)測(cè)步相同，根據(jù)上一時(shí)刻狀態(tài)變量(即葉綠素a濃度)后驗(yàn)集合，使用湖泊藻類(lèi)動(dòng)態(tài)模型模擬該時(shí)刻狀態(tài)變量先驗(yàn)集合(式(1))，對(duì)先驗(yàn)集合求期望(式(5))即可得狀態(tài)變量的模擬值，該模擬值在該時(shí)刻實(shí)測(cè)數(shù)據(jù)可用之前都有效.

②分析：使用擾動(dòng)分析法定位該時(shí)刻敏感參數(shù). 根據(jù)擾動(dòng)分析法流程，本研究對(duì)湖泊藻類(lèi)動(dòng)態(tài)模型每個(gè)待分析參數(shù)的擾動(dòng)范圍為：±1%、±5%、±10%、±15%、±20%，共擾動(dòng)10次. 使用(式(6))計(jì)算出該時(shí)刻每個(gè)參數(shù)的敏感度S，將敏感度|S|>0.5作為敏感參數(shù)輸出到更新步.

③更新：當(dāng)該時(shí)刻實(shí)測(cè)數(shù)據(jù)可用時(shí)即進(jìn)入更新步，使用集合卡爾曼濾波更新?tīng)顟B(tài)變量和敏感參數(shù)集合. 先更新?tīng)顟B(tài)變量，根據(jù)(式(7))更新?tīng)顟B(tài)變量集合，狀態(tài)變量先驗(yàn)集合被更新為后驗(yàn)集合，對(duì)其求期望(式(10))即可得該時(shí)刻狀態(tài)變量的同化值；再更新敏感參數(shù)，由于耦合了參數(shù)敏感性分析，每一步需要更新的參數(shù)維度不同，因此將參數(shù)維度固定為一維，然后依次調(diào)整，如有多個(gè)參數(shù)處于敏感狀態(tài)則依次執(zhí)行調(diào)整參數(shù)步(式(11))，如沒(méi)有參數(shù)處于敏感狀態(tài)則不執(zhí)行調(diào)整參數(shù)步. 考慮到集合卡爾曼濾波的假設(shè)是各參數(shù)誤差相互獨(dú)立且服從高斯分布[11]，因此涉及到多個(gè)敏感參數(shù)需要調(diào)整時(shí)，依次調(diào)整每個(gè)參數(shù)效果等同于同步調(diào)整；然后使用(式(13))對(duì)所有參數(shù)集合求期望即可得該時(shí)刻模型參數(shù)的同化值. 該時(shí)刻狀態(tài)變量和參數(shù)的同化值將作為模型下一時(shí)刻預(yù)測(cè)步的初值，依此遞推，運(yùn)行流程如圖3所示.

圖3 新型數(shù)據(jù)同化算法流程圖[18]Fig.3 New data assimilation algorithm flow chart[18]

表3 新型數(shù)據(jù)同化算法公式*

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

本研究提出的新型數(shù)據(jù)同化算法以集合卡爾曼濾波為主體改進(jìn)而成. 為測(cè)試其在傳統(tǒng)集合卡爾曼濾波基礎(chǔ)上的提升效果，同步設(shè)計(jì)基于傳統(tǒng)集合卡爾曼濾波數(shù)據(jù)同化算法方案用于對(duì)比分析. 目前集合卡爾曼濾波在水文水質(zhì)模型中主要有2種同化方案，一是僅同化狀態(tài)變量[35]，二是同步同化狀態(tài)變量和模型參數(shù)[18]. 為詳細(xì)對(duì)比差異，分別將這兩種方案應(yīng)用到模型查看其對(duì)模型的提升效果. 同時(shí)，還需一種完全憑借湖泊藻類(lèi)動(dòng)態(tài)模型的模擬方案作為評(píng)價(jià)上述3種數(shù)據(jù)同化方案效果的基準(zhǔn). 因此，本研究共有了4種方案，包括1個(gè)純模型方案和3個(gè)數(shù)據(jù)同化方案(僅同化狀態(tài)變量、同步同化狀態(tài)變量和模型參數(shù)及新型數(shù)據(jù)同化算法方案)，下面詳細(xì)介紹4種方案的策略及運(yùn)行步驟.

1)Reference：不使用數(shù)據(jù)同化算法，僅依據(jù)湖泊藻類(lèi)動(dòng)態(tài)模型模擬葉綠素a濃度. 將2020年1月1日實(shí)測(cè)的水質(zhì)數(shù)據(jù)和2020年1月2日的氣象數(shù)據(jù)作為初值輸入，將該階段模擬的葉綠素a濃度作為下一階段模型的初值，依此遞推，持續(xù)模擬至2020年6月1日. 該方案旨在了解湖泊藻類(lèi)動(dòng)態(tài)模型本身的模擬精度，為3種數(shù)據(jù)同化方案的模擬值提供對(duì)比基準(zhǔn).

2)EnKF_A：使用集合卡爾曼濾波數(shù)據(jù)同化算法同化狀態(tài)變量(即葉綠素a濃度)，不同化模型參數(shù). 輸入初值時(shí)使用蒙特卡洛方法分別對(duì)狀態(tài)變量和邊界條件添加誤差并生成集合，模擬出下一階段葉綠素a濃度集合，對(duì)集合求期望即可獲得模擬值. 當(dāng)獲得葉綠素a濃度實(shí)測(cè)值時(shí)，使用集合卡爾曼濾波同化狀態(tài)變量集合，對(duì)集合求期望獲得同化值，作為模型下一階段狀態(tài)變量輸入，依此遞推.

3)EnKF_B：使用集合卡爾曼濾波數(shù)據(jù)同化算法同步同化狀態(tài)變量和模型參數(shù). 考慮到使用的湖泊藻類(lèi)動(dòng)態(tài)模型參數(shù)眾多，僅同步同化一個(gè)參數(shù)，即Umax(藻類(lèi)最大生長(zhǎng)速率)，其為模型主控方程參數(shù)，是模型初值階段最敏感的參數(shù). 初值輸入時(shí)還需要同步生成待同化參數(shù)集合，模擬過(guò)程同上方案相同，當(dāng)獲得葉綠素a濃度實(shí)測(cè)值時(shí)，使用集合卡爾曼濾波同步同化狀態(tài)變量和參數(shù)集合，分別對(duì)集合求期望獲得狀態(tài)變量和參數(shù)的同化值，作為模型下一階段狀態(tài)變量和該參數(shù)輸入，依此遞推.

4)EnKF_C：使用本研究開(kāi)發(fā)的新型數(shù)據(jù)同化算法同步同化狀態(tài)變量和敏感參數(shù). 考慮到湖泊藻類(lèi)動(dòng)態(tài)模型參數(shù)較多(18個(gè)參數(shù))，全部用于敏感性分析影響同化效率，因此選取其中8個(gè)總體較為敏感的參數(shù)，用于新型數(shù)據(jù)同化算法方案中模型每一階段的敏感性分析，表4列出了用于敏感性分析的8個(gè)參數(shù)及其初值. 其中p4(藻類(lèi)最大生長(zhǎng)速率)即為同步同化模型參數(shù)方案(EnKF_B)的同化參數(shù)，該參數(shù)是模型主控方程的參數(shù)，在初始狀態(tài)是模型最為敏感的參數(shù)，其余7個(gè)參數(shù)也是模型各模塊的主要控制部分. 因此，通過(guò)分析上述8個(gè)參數(shù)的敏感性能表征模型的總體敏感程度. 新算法初值輸入時(shí)同步生成待分析敏感度參數(shù)集合，模擬過(guò)程依然相同，當(dāng)獲得葉綠素a濃度實(shí)測(cè)值時(shí)，先使用擾動(dòng)分析法定位該階段模型的敏感參數(shù)，然后使用集合卡爾曼濾波依次同化狀態(tài)變量和敏感參數(shù)集合，分別對(duì)集合求期望獲得狀態(tài)變量和敏感參數(shù)的同化值，作為模型下一階段狀態(tài)變量和敏感參數(shù)輸入，依此遞推.

表4 敏感性分析的8個(gè)參數(shù)[29]

Tab.4 Eight parameters for sensitivity analysis[29]

編號(hào)參數(shù)描述單位初值p1?溫度因子1.08p2α水及非藻類(lèi)吸收短波輻射的平均消光系數(shù)m-10.45p3K藻類(lèi)沉降速率m/d0.0864p4Umax藻類(lèi)最大生長(zhǎng)速率d-11.145p5GRmax浮游動(dòng)物最大捕食速率d-10.09p6kr藻類(lèi)呼吸引起葉綠素a減少速率系數(shù)d-10.17p7Topt最適藻類(lèi)生長(zhǎng)的水溫℃27p8ke藻類(lèi)排泄引起葉綠素a減少速率系數(shù)d-10.01

以上4種方案中，Reference為純模型模擬方案，其模擬值作為評(píng)價(jià)3種數(shù)據(jù)同化方案對(duì)模型精度提升效果的基準(zhǔn)；EnKF_A和EnKF_B為傳統(tǒng)集合卡爾曼濾波數(shù)據(jù)同化算法2種主要應(yīng)用方案；EnKF_C即為本研究提出的新型數(shù)據(jù)同化算法.

1.6 數(shù)據(jù)同化效果評(píng)價(jià)

為評(píng)價(jià)模擬精度，采用納什系數(shù)(Nash-sutcliffe model efficiency,NSE)衡量模擬值和實(shí)測(cè)值的符合程度；采用均方根誤差(root mean square error,RMSE)、平均絕對(duì)誤差(mean absolute error,MAE)和平均絕對(duì)百分比誤差(mean absolute percentage error,MAPE)衡量模擬值和實(shí)測(cè)值的誤差[36].

2 結(jié)果與分析

由于集合卡爾曼濾波中測(cè)量誤差和模型誤差兩參數(shù)無(wú)法顯式給出[11]. 為保證數(shù)據(jù)同化效果，先使用數(shù)值模擬實(shí)驗(yàn)測(cè)試；然后依次實(shí)施4種方案并根據(jù)性能指標(biāo)評(píng)價(jià)模擬精度和誤差，對(duì)比分析4種方案模擬效果.

2.1 數(shù)據(jù)同化算法參數(shù)率定

集合卡爾曼濾波中模型誤差和測(cè)量誤差兩參數(shù)相互影響，需同步率定[11]. 模型誤差主要來(lái)源于模型結(jié)構(gòu)、參數(shù)及邊界條件誤差等[37]，為合理評(píng)估模型誤差，使用模型模擬值與實(shí)測(cè)值計(jì)算8個(gè)自動(dòng)監(jiān)測(cè)站MAPE，取其中最小值(10%)和最大值(50%)作為模型誤差率定區(qū)間邊界；測(cè)量誤差主要來(lái)源于測(cè)量?jī)x器及環(huán)境(氣溫、風(fēng)力等)條件等，綜合考慮儀器精度及測(cè)量過(guò)程，測(cè)量誤差率定區(qū)間邊界取定為1%～10%. 同步率定時(shí)使用RMSE評(píng)價(jià)數(shù)據(jù)同化算法模擬誤差，RMSE取值越小說(shuō)明同化效果越佳. 根據(jù)數(shù)值模擬結(jié)果所示(圖4)，當(dāng)模型誤差和測(cè)量誤差分別取18%和3%時(shí)，數(shù)據(jù)同化效果最佳，由此即確定了集合卡爾曼濾波模型誤差和測(cè)量誤差取值. 研究表明集合卡爾曼濾波過(guò)程中集合個(gè)數(shù)越多提升效果越好[17]，但考慮到效率問(wèn)題，本研究中集合卡爾曼濾波的集合個(gè)數(shù)取定為100，此集合數(shù)已被證明在大多數(shù)模型中都足夠[38]. 集合卡爾曼濾波的初值對(duì)同化結(jié)果影響不大，因此將模型初始狀態(tài)的誤差設(shè)置為20%.

圖4 不同模型誤差與測(cè)量誤差背景下的同化結(jié)果Fig.4 Assimilation results in case of different models and measurement errors

根據(jù)上述討論，3種數(shù)據(jù)同化方案中集合卡爾曼濾波的參數(shù)如下：模型誤差為18%，測(cè)量誤差為3%，集合個(gè)數(shù)為100，初始狀態(tài)誤差為20%.

2.2 4種方案結(jié)果對(duì)比分析

圖5展示了8個(gè)自動(dòng)監(jiān)測(cè)站2020年1月1日-6月1日逐日葉綠素a濃度實(shí)測(cè)值和4種方案的模擬曲線(xiàn)，并給出了4種方案模擬曲線(xiàn)的一項(xiàng)擬合指標(biāo)(NSE)和3項(xiàng)誤差指標(biāo)(RMSE、MAE、MAPE). 從圖5葉綠素a濃度實(shí)測(cè)值看出，時(shí)間上，多數(shù)自動(dòng)監(jiān)測(cè)站葉綠素a濃度峰值出現(xiàn)在4月份；空間上，巢湖西部和南部水域的葉綠素a濃度總體高于東部和北部水域；從葉綠素a濃度極大值模擬看，數(shù)據(jù)同化后的模擬值大部分略低于實(shí)測(cè)值，因多數(shù)藻類(lèi)實(shí)測(cè)的極大值梯度大、回落迅速，模型往往無(wú)法及時(shí)反應(yīng)，其模擬值具有滯后性且難以達(dá)到峰值高度，經(jīng)過(guò)濾波后會(huì)取到一個(gè)中間值，因此濾波模擬結(jié)果多數(shù)略低于實(shí)測(cè)極大值而又高于模型模擬結(jié)果；另外，還需注意藍(lán)藻水華暴發(fā)期藻類(lèi)上浮堆積在水面可能會(huì)遮住自動(dòng)監(jiān)測(cè)探頭，監(jiān)測(cè)數(shù)據(jù)不夠準(zhǔn)確，存在出現(xiàn)異常峰值的可能. 以葉綠素a濃度實(shí)測(cè)值為依據(jù)，下面詳細(xì)分析4種方案模擬效果.

1)Reference：模型模擬的葉綠素a濃度曲線(xiàn)僅能夠反映實(shí)測(cè)值大致變化趨勢(shì)，且在模擬預(yù)測(cè)葉綠素a濃度峰值上難以發(fā)揮作用，模擬結(jié)果誤差較大. 從性能評(píng)價(jià)指標(biāo)得出模型模擬精度在8個(gè)自動(dòng)監(jiān)測(cè)站差距較大，新河入湖區(qū)站模擬值的NSE僅為0.1，而東半湖湖心站模擬值NSE則達(dá)到了0.69. 8個(gè)自動(dòng)監(jiān)測(cè)站平均NSE為0.49，平均RMSE為7.31 μg/L，平均MAPE為44%，可見(jiàn)模型總體模擬效果并不理想.

2)EnKF_A：從集合卡爾曼濾波同化狀態(tài)變量方案模擬曲線(xiàn)看出，經(jīng)過(guò)同化狀態(tài)變量的模型模擬曲線(xiàn)更加貼近實(shí)測(cè)值. 從性能評(píng)價(jià)指標(biāo)得出，8個(gè)自動(dòng)監(jiān)測(cè)站模擬值的平均NSE提升到了0.63，平均RMSE降低了2.47 μg/L，平均MAPE也降低了9%，說(shuō)明同化狀態(tài)變量能夠提升模型模擬精度.

3)EnKF_B：從集合卡爾曼濾波同步同化狀態(tài)變量和模型參數(shù)方案模擬曲線(xiàn)看出，相較于僅同化狀態(tài)變量，同步同化模型參數(shù)后模型模擬精度提升有限. 盡管8個(gè)自動(dòng)監(jiān)測(cè)站平均NSE提升了0.08，但注意到其中3個(gè)站(湖濱站、東半湖湖心、兆河入湖區(qū))的NSE不增反減，而平均NSE的提升主要來(lái)自新河入湖區(qū)站，從0.26提升到0.53. 這反映不考慮模型參數(shù)敏感性動(dòng)態(tài)變化而對(duì)固定參數(shù)調(diào)整難以有效校正模型運(yùn)行軌跡，反而會(huì)因不當(dāng)調(diào)整參數(shù)使同化效果更差，存在給模型模擬帶來(lái)額外不確定性的隱患.

4)EnKF_C：從本研究提出的耦合參數(shù)敏感性分析的新型數(shù)據(jù)同化方案模擬曲線(xiàn)看出，新型數(shù)據(jù)同化方案能夠更加準(zhǔn)確地模擬預(yù)測(cè)葉綠素a濃度，在上述2種數(shù)據(jù)同化方案上有了進(jìn)一步提升. 從性能評(píng)價(jià)指標(biāo)得出該方案的模型模擬值平均NSE達(dá)到了0.76，平均RMSE為4.27 μg/L，平均MAPE下降到23%. 無(wú)論總體上還是在每個(gè)自動(dòng)監(jiān)測(cè)站，考慮了參數(shù)敏感性的新型數(shù)據(jù)同化算法對(duì)傳統(tǒng)集合卡爾曼濾波數(shù)據(jù)同化效果都有穩(wěn)定的提升，對(duì)葉綠素a濃度峰值有更加準(zhǔn)確的模擬. 考慮到本研究設(shè)計(jì)的3種數(shù)據(jù)同化方案參數(shù)條件都相同，因此能夠確定模型模擬精度的提升來(lái)源于新型數(shù)據(jù)同化算法對(duì)敏感參數(shù)的實(shí)時(shí)定位和調(diào)整. 這說(shuō)明耦合了參數(shù)敏感性分析的新型數(shù)據(jù)同化算法能夠有效定位模型運(yùn)行每階段的敏感參數(shù)，調(diào)整模型運(yùn)行軌跡，增加模型穩(wěn)定性，提高模型葉綠素a濃度模擬精度.

3 討論

案例測(cè)試結(jié)果表明：巢湖藻類(lèi)模型的參數(shù)敏感性存在顯著的動(dòng)態(tài)變化過(guò)程，提出新型數(shù)據(jù)同化算法能夠?qū)崟r(shí)跟蹤模型敏感參數(shù)，實(shí)現(xiàn)模型模擬精度的最大化提升. 在新型數(shù)據(jù)同化算法方案中，根據(jù)每階段模型模擬值，將待分析敏感性的8個(gè)參數(shù)依次輸入擾動(dòng)分析模塊即可獲得其敏感度S. 以湖濱站1-5月參數(shù)敏感性的動(dòng)態(tài)變化(圖6)為例，p1(溫度因子)和p4(藻類(lèi)最大生長(zhǎng)速率)在多數(shù)時(shí)段敏感度等級(jí)都在Ⅲ級(jí)(|S|≥0.5，一般敏感)及以上，敏感性總體較高. 尤其是p4，在模型運(yùn)行時(shí)間段有62 d的敏感度都在Ⅴ級(jí)(|S|≥1.0，極度敏感)以上，是模型最敏感的參數(shù)，但也有39 d處于不敏感狀態(tài). 表明：受邊界條件等要素影響，參數(shù)敏感度在模型運(yùn)行各階段不斷變化，敏感參數(shù)在模型動(dòng)態(tài)運(yùn)行過(guò)程中存在轉(zhuǎn)為不敏感的可能. 其它參數(shù)雖多不敏感，但隨著模型的運(yùn)行也有敏感時(shí)刻. 例如，p3(藻類(lèi)沉降速率)在模型運(yùn)行前20 d都不敏感，但在1月21日變?yōu)闃O度敏感參數(shù)(|S|≥1.0)，其在當(dāng)天敏感程度超過(guò)了p4(1.0≥|S|≥0.75)，表明隨著模型的動(dòng)態(tài)運(yùn)行，原本不敏感參數(shù)也可能受邊界條件影響轉(zhuǎn)為敏感參數(shù).

圖6 基于擾動(dòng)分析法的參數(shù)敏感性動(dòng)態(tài)變化過(guò)程(湖濱站，參數(shù)p1～p8含義見(jiàn)表4)Fig.6 The parameter sensitivity dynamic change process based on the disturbance analysis method (Hubin Station, the meaning of parameters p1-p8 are shown in Tab.4)

通過(guò)對(duì)比分析可知，新型數(shù)據(jù)同化算法能夠顯著提升模型模擬精度的原因在于耦合了參數(shù)敏感性分析，使其能夠?qū)崟r(shí)定位模型敏感參數(shù)并有效調(diào)整，校正模型運(yùn)行軌跡，使模型模擬更加精確. 在本研究的湖泊藻類(lèi)動(dòng)態(tài)模型應(yīng)用中，參數(shù)敏感性分析方法選擇了成熟的擾動(dòng)分析法，從模型中選取了8個(gè)具有代表性參數(shù)(表4)用于敏感性分析，定位到敏感參數(shù)后使用集合卡爾曼濾波更新敏感參數(shù)，實(shí)時(shí)校正模型運(yùn)行軌跡. 而傳統(tǒng)集合卡爾曼濾波數(shù)據(jù)同化算法則僅能夠調(diào)整固定的參數(shù)，難以應(yīng)對(duì)模型運(yùn)行中的敏感參數(shù)動(dòng)態(tài)變化，無(wú)法有效校正模型運(yùn)行軌跡，甚至不當(dāng)?shù)膮?shù)調(diào)整還會(huì)增加模型不確定性. 通過(guò)兩種模型同化方案(EnKF_B、EnKF_C)中參數(shù)值動(dòng)態(tài)變化，能夠看出新算法相較于傳統(tǒng)算法對(duì)模型參數(shù)的優(yōu)化調(diào)整過(guò)程，考慮到站點(diǎn)多、時(shí)序長(zhǎng)，圖7僅展示湖濱站1月份兩種方案中參數(shù)值的動(dòng)態(tài)變化率，選擇湖濱站是為了方便從圖6中同步觀察各參數(shù)1月份的敏感度. 從圖7a看出，由于固定同化參數(shù)p4(藻類(lèi)最大生長(zhǎng)速率)，EnKF_B方案僅能且不斷的調(diào)整該參數(shù)，即使是在該參數(shù)處于不敏感狀態(tài)(如1月5日)依然調(diào)整，而其它參數(shù)即使處于敏感狀態(tài)也不會(huì)進(jìn)行調(diào)整. 而從圖7b看出，EnKF_C方案則能夠根據(jù)模型敏感參數(shù)動(dòng)態(tài)變化靈活改變對(duì)參數(shù)的調(diào)整，合理應(yīng)對(duì)敏感參數(shù)動(dòng)態(tài)變化. 以p4為例，當(dāng)其處于敏感狀態(tài)時(shí)(如1月2-4日)會(huì)及時(shí)進(jìn)行調(diào)整，當(dāng)其處于非敏感狀態(tài)時(shí)(如1月6日)則不會(huì)進(jìn)行調(diào)整. 由此可見(jiàn)，本研究提出的新型數(shù)據(jù)同化算法正是通過(guò)耦合參數(shù)敏感性分析，實(shí)時(shí)動(dòng)態(tài)調(diào)整敏感參數(shù)，實(shí)現(xiàn)模型模擬精度的穩(wěn)定提升，具有廣泛應(yīng)用前景.

圖7 不同數(shù)據(jù)同化方案(EnKF_B(a)和EnKF_C(b))的參數(shù)動(dòng)態(tài)變化(湖濱站，2020年1月，參數(shù)p1～p8含義見(jiàn)表4)Fig.7 The change of parameter values for the assimilation strategies of EnKF_B (a) and EnKF_C (b)(Hubin Station, January 2020, parameters (p1-p8) descriptions can be found in Tab.4)

本研究模擬時(shí)段為2020年1-5月，主要是由于在5-10月的藍(lán)藻水華暴發(fā)期間，藻類(lèi)水華堆積可能影響探頭準(zhǔn)確測(cè)定，自動(dòng)監(jiān)測(cè)數(shù)據(jù)質(zhì)量受影響；同時(shí)，根據(jù)孔繁翔等提出的湖泊藍(lán)藻水華形成四階段理論[39]，本研究模擬時(shí)間段(1-5月)主要處于藻類(lèi)復(fù)蘇和生物量增加的初期，也是決定全年藍(lán)藻水華情勢(shì)的重要時(shí)期，開(kāi)展此時(shí)段的藍(lán)藻水華預(yù)測(cè)也有一定意義.

4 結(jié)論

本研究通過(guò)耦合數(shù)據(jù)同化算法與敏感參數(shù)識(shí)別算法，實(shí)現(xiàn)了模型敏感參數(shù)的實(shí)時(shí)識(shí)別與同化，解決了傳統(tǒng)數(shù)據(jù)同化算法難以適應(yīng)模型敏感參數(shù)動(dòng)態(tài)變化的問(wèn)題；改進(jìn)算法應(yīng)用于巢湖藻類(lèi)動(dòng)態(tài)模擬，測(cè)試結(jié)果表明：改進(jìn)算法實(shí)現(xiàn)了湖泊藻類(lèi)動(dòng)態(tài)模型與高頻自動(dòng)監(jiān)測(cè)數(shù)據(jù)的實(shí)時(shí)融合，顯著提升了模型模擬精度(模擬誤差減少42%)，對(duì)藻類(lèi)峰值的模擬精度提升尤為明顯；算法與模型采用松散耦合模式，可移植性強(qiáng)，可應(yīng)用于其它水生態(tài)環(huán)境模型的數(shù)據(jù)同化，隨著我國(guó)河湖多源監(jiān)測(cè)數(shù)據(jù)的劇增，算法具有廣闊應(yīng)用前景.