基于數(shù)據(jù)同化的深水湖庫(kù)水溫中短期預(yù)報(bào)

孫博聞，楊晰淯，暴 柱，劉曉波，劉 暢，高學(xué)平

(1.天津大學(xué) 水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300350；2.水利部海河水利委員會(huì)引灤工程管理局 海河流域?yàn)春铀|(zhì)監(jiān)測(cè)中心，河北 唐山 064309；3.中國(guó)水利水電科學(xué)研究院 水生態(tài)環(huán)境研究所，北京 100038)

1 研究背景

深水湖庫(kù)通常存在季節(jié)性溫度分層現(xiàn)象，水溫不僅影響水體的運(yùn)動(dòng)規(guī)律，也影響著水體中生化反應(yīng)以及水生生物的新陳代謝過(guò)程，包括溶解氧飽和度[1]、微生物和藻類生長(zhǎng)[2]、沉積物向水體釋放化學(xué)成分[3-5]以及魚(yú)類對(duì)棲息地適宜性[2，6]；由溫度分層引起溶解氧等環(huán)境指標(biāo)的分層還會(huì)誘發(fā)水環(huán)境水生態(tài)問(wèn)題[7-8]，進(jìn)而威脅供水安全及下游生態(tài)系統(tǒng)健康。在氣候變化與水庫(kù)調(diào)控的共同影響下，上述變化更加趨于復(fù)雜[9-10]，因此研究水溫變化對(duì)湖庫(kù)水質(zhì)管理與生態(tài)安全十分必要。

當(dāng)前對(duì)水庫(kù)水溫的研究大多關(guān)注其長(zhǎng)期演變規(guī)律[11-13]，中短期預(yù)報(bào)方面的研究相對(duì)較少。事實(shí)上，在未來(lái)10天這一時(shí)間尺度上預(yù)報(bào)湖庫(kù)水溫變化對(duì)管理者具有重要意義[14]，如Weber等[15]利用水溫-深度變化數(shù)據(jù)動(dòng)態(tài)確定Grosse Dhuenn水庫(kù)取水深度，在優(yōu)化下泄水溫同時(shí)還抑制了水庫(kù)缺氧的發(fā)展；Baracchini等[16]在日內(nèi)瓦湖的觀測(cè)發(fā)現(xiàn)，強(qiáng)風(fēng)驅(qū)動(dòng)產(chǎn)生的劇烈涌升流使部分湖區(qū)的上層水溫在2天內(nèi)驟降10 ℃，直接影響下游的供水與河流生態(tài)。此外水溫垂向分布還決定了湖庫(kù)熱分層強(qiáng)度，進(jìn)而決定湖庫(kù)的分層或混合狀態(tài)[17]。秋季湖庫(kù)由分層狀態(tài)過(guò)渡到混合狀態(tài)時(shí)，庫(kù)底積累的營(yíng)養(yǎng)物質(zhì)和重金屬在整個(gè)水體混合，此時(shí)冷雨等氣象條件誘發(fā)的翻庫(kù)過(guò)程則會(huì)使垂向混合在短時(shí)間內(nèi)發(fā)生[18]，這類由剖面水溫驟變引發(fā)的水動(dòng)力過(guò)程也會(huì)對(duì)湖庫(kù)生態(tài)安全產(chǎn)生較大沖擊。通過(guò)上述分析可見(jiàn)，中短期湖庫(kù)水溫預(yù)報(bào)能夠在調(diào)控湖庫(kù)下泄水溫這類常規(guī)管理中發(fā)揮作用，還能在應(yīng)對(duì)極端天氣(如風(fēng)暴)等誘發(fā)的翻庫(kù)、涌升流這類威脅供水與生態(tài)安全的應(yīng)急管理工作中提供決策支撐。

數(shù)學(xué)模型是湖庫(kù)水環(huán)境管理與決策分析的重要工具，初始條件與模型參數(shù)的確定是提升模型預(yù)報(bào)準(zhǔn)確性面臨的挑戰(zhàn)[19-20]，隨著實(shí)時(shí)觀測(cè)數(shù)據(jù)的獲取和傳輸?shù)玫斤w速發(fā)展，使得將原型觀測(cè)數(shù)據(jù)引入水環(huán)境預(yù)報(bào)，以提高預(yù)報(bào)精度成為可能[21]。數(shù)據(jù)同化融合觀測(cè)數(shù)據(jù)與模型計(jì)算結(jié)果，綜合考慮模型結(jié)構(gòu)誤差、邊界條件誤差和觀測(cè)誤差，對(duì)模型狀態(tài)變量進(jìn)行最優(yōu)估計(jì)，實(shí)現(xiàn)了模型狀態(tài)和參數(shù)不斷更新[19，22]，目前已應(yīng)用于提高河流磷遷移估計(jì)和湖泊富營(yíng)養(yǎng)化等方面的模擬精度[23-24]。集合卡爾曼濾波是應(yīng)用最為廣泛的數(shù)據(jù)同化算法，該方法能夠提供狀態(tài)量的均值及其相應(yīng)的誤差協(xié)方差，同時(shí)由于引入了集合的思想且不需要伴隨或線性算子，因此具有適用于非線性系統(tǒng)、程序設(shè)計(jì)相對(duì)簡(jiǎn)單等特點(diǎn)[25]。

本文以引灤入津工程源頭水庫(kù)—大黑汀水庫(kù)為研究對(duì)象，建立大黑汀水庫(kù)的二維水動(dòng)力模型，基于集合卡爾曼濾波算法構(gòu)建可綜合考慮模型參數(shù)、邊界條件以及觀測(cè)數(shù)據(jù)不確定性的湖庫(kù)水溫?cái)?shù)據(jù)同化系統(tǒng)。系統(tǒng)利用未來(lái)水庫(kù)調(diào)度數(shù)據(jù)與氣象數(shù)據(jù)等作為預(yù)報(bào)條件，能夠進(jìn)行未來(lái)1～10 d內(nèi)的湖庫(kù)水溫預(yù)報(bào)，高精度的水溫中短期預(yù)報(bào)方法可以為深水湖庫(kù)供水與生態(tài)安全提供理論與技術(shù)支撐。

2 研究方法

2.1 CE-QUAL-W2模型CE-QUAL-W2是一個(gè)考慮縱向和垂向的二維水動(dòng)力水質(zhì)模型(簡(jiǎn)稱W2模型)[26]，模型假設(shè)水體橫向平均，適宜模擬河道型水庫(kù)的水動(dòng)力與水質(zhì)變化，模型主要控制方程如下。

連續(xù)性方程：

(1)

X方向動(dòng)量方程：

(2)

Z方向動(dòng)量方程：

(3)

狀態(tài)方程：

ρ=f(Tw，ΦTDS，ΦSS)

(4)

水面線方程：

(5)

質(zhì)量/熱量守恒方程：

(6)

式中：U與W為河道縱向與垂向流速，m/s；q為單寬流量，m3/s；B為網(wǎng)格寬度，m；g為重力加速度，m2/s；α為河底線與水平面夾角， °；η為水面線高程，m；Bη為水面寬度，m；ρ為密度，kg/m3；τxx與τxz與分別縱向流速的垂直梯度在x方向和z方向所產(chǎn)生的單位面積的剪切應(yīng)力，Pa；P為大氣壓，kPa；f(Tw，ΦTDS，Φss)為考慮水溫、總?cè)芙夤腆w、無(wú)機(jī)懸浮物的密度函數(shù)；h為沿河道底坡法線方向的水深，m；Φ為各組分濃度，g/m3；Dx為縱向擴(kuò)散系數(shù)，m/s2；Dz為垂向擴(kuò)散系數(shù)，m/s2；qΦ為入流組分濃度，g/(m3·s)；SΦ為源匯項(xiàng)，g/(m3·s)。

圖1 基于數(shù)據(jù)同化的湖庫(kù)水溫中短期預(yù)報(bào)工作流程

(7)

(8)

(9)

(10)

(11)

(12)

(13)

2.3 湖庫(kù)水溫中短期預(yù)報(bào)流程基于數(shù)據(jù)同化的湖庫(kù)水溫中短期預(yù)報(bào)包括四個(gè)階段(圖1)：(1)建模階段，根據(jù)獲取的地形、氣象等數(shù)據(jù)建立W2數(shù)學(xué)模型；(2)參數(shù)優(yōu)化階段，對(duì)影響W2水溫模擬的敏感參數(shù)增加噪聲擾動(dòng)，運(yùn)行N組W2模型，將水溫模擬值與觀測(cè)值輸入EnKF算法，得到關(guān)鍵參數(shù)的優(yōu)化結(jié)果，此階段模擬時(shí)長(zhǎng)不少于20 d；(3)變量同化階段，在得到的關(guān)鍵參數(shù)后驗(yàn)分布基礎(chǔ)上，為邊界條件增加噪聲擾動(dòng)以體現(xiàn)其不確定性，再次運(yùn)行N組W2模型，將水溫模擬值與觀測(cè)值二次輸入EnKF算法，得到水溫同化結(jié)果；(4)預(yù)報(bào)階段，將同化結(jié)果作為本階段的湖庫(kù)水溫初始條件與參數(shù)組合，在獲取未來(lái)氣象數(shù)據(jù)、水庫(kù)入流條件以及調(diào)度計(jì)劃作為預(yù)報(bào)邊界條件后，運(yùn)行CE-QUAL-W2模型實(shí)現(xiàn)對(duì)湖庫(kù)水溫在未來(lái)1～10 d的中短期預(yù)報(bào)。

3 研究實(shí)例

3.1 研究區(qū)域與數(shù)據(jù)來(lái)源大黑汀水庫(kù)位于河北省唐山市遷西縣灤河干流上，是引灤樞紐工程體系中的骨干水庫(kù)工程，通過(guò)承接其上游來(lái)水與潘家口水庫(kù)泄水為天津、唐山等地供水。水庫(kù)總庫(kù)容約3.37億m3，為大Ⅱ型水庫(kù)，庫(kù)區(qū)回水長(zhǎng)度約23 km，具有河道型水庫(kù)狹長(zhǎng)的地形特征(圖2)。壩前水深可達(dá)25 m，全庫(kù)表現(xiàn)出典型的季節(jié)性水溫分層特征(圖2(d))。水庫(kù)溫度分層直接影響著水體理化過(guò)程，是目前大黑汀水庫(kù)供水所面臨的難題之一[29]。對(duì)大壩下游而言，大黑汀水庫(kù)下泄水溫是控制河道水溫的主要因素，進(jìn)而影響著下游生態(tài)系統(tǒng)及耕種活動(dòng)；對(duì)庫(kù)區(qū)而言，水溫分層誘發(fā)的庫(kù)底缺氧環(huán)境使得沉積物中的污染物發(fā)生厭氧分解，加速水質(zhì)惡化并威脅供水安全[30]。本文水溫觀測(cè)數(shù)據(jù)包括長(zhǎng)期定點(diǎn)觀測(cè)與定期走航觀測(cè)(圖2)，定點(diǎn)觀測(cè)位于大黑汀水庫(kù)壩前500 m，在該斷面1、7和13 m等深度布置自動(dòng)水溫記錄儀，觀測(cè)頻率30 min/組；氣象數(shù)據(jù)源自國(guó)家氣象中心以及布設(shè)于水庫(kù)的在線氣象站；入庫(kù)水溫、流量及水庫(kù)調(diào)度數(shù)據(jù)由水利部海河水利委員會(huì)引灤工程管理局提供。

圖2 大黑汀水庫(kù)概況與2019年水溫觀測(cè)數(shù)據(jù)

3.2 大黑汀水庫(kù)水動(dòng)力模型本文將大黑汀水庫(kù)縱向網(wǎng)格間距(DLX)設(shè)為300 m，垂向?qū)雍穸?H)設(shè)為1 m，模型底部高程設(shè)為103 m，頂部高程為133 m，包含虛擬網(wǎng)格共劃分69個(gè)單元(segment)、最多33層(layer)，考慮到大黑汀水庫(kù)水面波動(dòng)較小，因此設(shè)水面坡度(slope)為0。W2模型中影響水溫計(jì)算的關(guān)鍵參數(shù)及取值列于表1。

表1 模型關(guān)鍵參數(shù)取值

對(duì)于垂向擴(kuò)散系數(shù)W2模型要求用戶選定計(jì)算模式后，由模型內(nèi)部計(jì)算得出，本文選定的計(jì)算模式為W2N。率定過(guò)程中發(fā)現(xiàn)，AX和DX對(duì)水溫結(jié)構(gòu)影響較小，因此選取默認(rèn)值1，這與李娟等[31]的研究結(jié)果一致。采用W2模型中的TERM算法計(jì)算水面熱量交換，該算法考慮了日照短波與長(zhǎng)波輻射的直接入射和反射、水面-大氣反向熱輻射、空氣-水面熱傳導(dǎo)以及蒸發(fā)潛熱幾方面的影響。率定結(jié)果表明模型模擬的水動(dòng)力過(guò)程與實(shí)測(cè)值符合較好，可用于模型計(jì)算。模型參數(shù)率定與驗(yàn)證過(guò)程參照姚嘉偉[32]研究成果，本文不再贅述。

3.3 數(shù)據(jù)同化方案設(shè)計(jì)合理的同化方案與參數(shù)設(shè)置是保證同化效果的前提，EnKF算法利用集合思想解決了實(shí)際應(yīng)用中背景協(xié)方差矩陣估計(jì)和預(yù)報(bào)困難的問(wèn)題，集合數(shù)越多，集合均值協(xié)方差與真值協(xié)方差越接近，但計(jì)算負(fù)荷會(huì)顯著增加，集合數(shù)過(guò)少，則可能導(dǎo)致模型誤差協(xié)方差估計(jì)錯(cuò)誤[33]。為此本文設(shè)置7組集合數(shù)(3、10、20、50、100、200和400)進(jìn)行數(shù)值試驗(yàn)，尋找可兼顧模擬精度與效率的集合數(shù)。數(shù)據(jù)同化通常包括只同化狀態(tài)變量以及同時(shí)同化狀態(tài)變量與模型參數(shù)兩種模式，根據(jù)徐興亞等[24]等的研究，后者表現(xiàn)通常優(yōu)于前者，因此本文也采用該模式。參考相關(guān)研究，選擇底部熱交換系數(shù)(CBHE)、太陽(yáng)輻射吸收系數(shù)(BETA)和風(fēng)遮蔽系數(shù)(WSC)這三個(gè)對(duì)W2模型水溫模擬最敏感參數(shù)[34-35]，作為參數(shù)優(yōu)化階段考慮的模型參數(shù)。雖然數(shù)據(jù)同化方法能夠在理論上實(shí)現(xiàn)模型模擬結(jié)果和觀測(cè)數(shù)據(jù)的最佳融合，但這種“最佳”與模型的模擬誤差協(xié)方差和觀測(cè)誤差協(xié)方差的精度密切相關(guān)，因此選擇合適的模擬誤差和觀測(cè)誤差對(duì)提高同化精度同樣至關(guān)重要。參考李港等[36]研究，將觀測(cè)誤差與模擬誤差分別設(shè)置為1%、10%、20%、30%，篩選觀測(cè)誤差和模擬誤差的最優(yōu)組合。

在此基礎(chǔ)上，選擇W2模型的入流流量以及氣象數(shù)據(jù)中的太陽(yáng)輻射熱量和風(fēng)速作為增加噪聲擾動(dòng)的邊界條件，以壩前測(cè)點(diǎn)1、7和13 m水深處的水溫觀測(cè)為同化數(shù)據(jù)(觀測(cè)頻率12 h/次)，對(duì)W2模型的參數(shù)和水溫進(jìn)行數(shù)據(jù)同化。首先將4月28日至5月9日的數(shù)據(jù)進(jìn)行同化模擬，再自5月10日起開(kāi)展基于數(shù)據(jù)同化結(jié)果和未來(lái)氣象數(shù)據(jù)、水庫(kù)入流條件以及調(diào)度計(jì)劃作為預(yù)報(bào)邊界條件的水溫預(yù)報(bào)，每次預(yù)報(bào)未來(lái)1～10 d的水溫，直至10月27日(預(yù)報(bào)流程詳見(jiàn)圖1)。采用均方根誤差(RMSE)、一致性系數(shù)(IOA)和百分比偏差系數(shù)(PBIAS)評(píng)價(jià)數(shù)據(jù)同化系統(tǒng)的表現(xiàn)。IOA用來(lái)表征模擬值和觀測(cè)值的一致性，IOA取值在0～1之間，越接近1表示二者之間的一致程度越高。PBIAS結(jié)果以百分比定量化給出模擬值比觀測(cè)值整體被低估或高估的平均趨勢(shì)[37]，PBIAS為0表示模擬效果最佳，PBIAS>0表示模擬值傾向于低估，反之PBIAS<0表示模型傾向于高估。IOA和PBIAS的計(jì)算公式分別為：

(14)

(15)

4 結(jié)果與分析

4.1 數(shù)據(jù)同化方案確定對(duì)壩前斷面處不同集合數(shù)下的數(shù)據(jù)同化水溫模擬結(jié)果與實(shí)測(cè)水溫進(jìn)行對(duì)比(圖3)，可以看出隨著集合數(shù)增大RMSE逐漸減小，當(dāng)集合數(shù)達(dá)到100后RMSE減小幅度大幅降低，同化結(jié)果趨于穩(wěn)定。將觀測(cè)誤差與模擬誤差分別設(shè)置為1%、10%、20%、30%共4種情況，圖3展示了不同的誤差組合下的同化結(jié)果?？梢钥闯?，當(dāng)觀測(cè)誤差取1%，模型誤差取10%、20%時(shí)，數(shù)據(jù)同化結(jié)果的RMSE均較小。綜上，本文的數(shù)據(jù)同化方案為：集合數(shù)取100，模擬誤差與觀測(cè)誤差分別取10%和1%。

圖3 不同同化方案模擬結(jié)果對(duì)比

4.2 數(shù)據(jù)同化效果圖4為有無(wú)數(shù)據(jù)同化情況下，壩前500 m斷面1、7和13 m水深處W2模型模擬結(jié)果與數(shù)據(jù)同化系統(tǒng)模擬結(jié)果相比實(shí)測(cè)水溫的偏差值?？梢钥闯觯诮?jīng)過(guò)EnKF算法對(duì)模型參數(shù)和狀態(tài)變量同時(shí)進(jìn)行更新后，數(shù)據(jù)同化系統(tǒng)在1、7和13 m水深處的平均模擬偏差分別為1.01、0.31和0.21 ℃，較無(wú)數(shù)據(jù)同化的模擬結(jié)果分別提升了68.8%、51.6%和41.2%，模擬精度顯著提升。本文參考Zhang等[38]提出的方法，以垂向溫度梯度1.5 ℃/m為閾值識(shí)別溫躍層，當(dāng)某日所有觀測(cè)水溫剖面全部出現(xiàn)溫躍層時(shí)，認(rèn)定進(jìn)入分層期，否則為混合期。針對(duì)分層期(6月5日至10月3日)、混合期和全時(shí)段分別評(píng)價(jià)W2模型與數(shù)據(jù)同化系統(tǒng)的模擬結(jié)果，結(jié)果列于表2?？梢钥闯觯瑹o(wú)論是分層期還是混合期，W2模型的均方根誤差(RMSE)明顯大于同化系統(tǒng)，最大可達(dá)到1.37 ℃；W2模型的一致性系數(shù)(IOA)相比1的偏差值也均大于同化系統(tǒng)；W2模型的百分比偏差系數(shù)(PBIAS)同樣大于同化系統(tǒng)?？梢?jiàn)數(shù)據(jù)同化方法可以有效的提升模型的模擬效果，為提高湖庫(kù)水溫的中短期預(yù)報(bào)精度奠定了基礎(chǔ)。

圖4 W2模型模擬結(jié)果與數(shù)據(jù)同化系統(tǒng)模擬結(jié)果相比實(shí)測(cè)水溫的偏差值(2019年)

表2 W2模型與數(shù)據(jù)同化系統(tǒng)的模擬結(jié)果評(píng)價(jià)

4.3 數(shù)據(jù)同化系統(tǒng)預(yù)報(bào)性能數(shù)據(jù)同化系統(tǒng)可為湖庫(kù)水溫預(yù)報(bào)提供可靠的模型參數(shù)與初始狀態(tài)，在此基礎(chǔ)上利用未來(lái)氣象數(shù)據(jù)、水庫(kù)入流條件以及調(diào)度計(jì)劃作為預(yù)報(bào)邊界條件，在考慮邊界條件的不確定性前提下，運(yùn)行W2模型進(jìn)行湖庫(kù)水溫預(yù)報(bào)。為分析數(shù)據(jù)同化系統(tǒng)的水溫預(yù)報(bào)性能，對(duì)分層期、混合期與全時(shí)段下各水深處未來(lái)1～10 d的水溫預(yù)報(bào)結(jié)果進(jìn)行評(píng)價(jià)，結(jié)果如圖5所示?？梢钥闯觯?與13 m處分層期的預(yù)報(bào)性能優(yōu)于混合期，7 m處兩階段預(yù)報(bào)性能各有優(yōu)劣。從RMSE指標(biāo)來(lái)看，隨著預(yù)報(bào)時(shí)段由1 d延長(zhǎng)至10 d，預(yù)報(bào)誤差由0.22～0.35 ℃增大至0.77～1.09 ℃，可見(jiàn)隨著預(yù)報(bào)時(shí)段延長(zhǎng)，不確定因素的累計(jì)效應(yīng)對(duì)預(yù)報(bào)會(huì)產(chǎn)生一定影響。從IOA指標(biāo)來(lái)看，分層期7 d預(yù)報(bào)期內(nèi)各深度的IOA值均超過(guò)0.65，具有較好的一致性；預(yù)報(bào)期延長(zhǎng)至10 d后，7 m處的IOA由0.78降至0.57，降幅遠(yuǎn)超1與13 m深度。近年來(lái)引灤工程常采用潘家口水庫(kù)大流量調(diào)度抑制大黑汀水庫(kù)缺氧發(fā)展[39]，該現(xiàn)象可能與此種調(diào)度有關(guān)。從PBIAS指標(biāo)來(lái)看，1和13 m處預(yù)報(bào)均傾向于高估，7 m處則傾向于低估，但高/低估均在±3%以內(nèi)。從各統(tǒng)計(jì)指標(biāo)結(jié)果可以看出，構(gòu)建的數(shù)據(jù)同化系統(tǒng)能夠較好的預(yù)報(bào)水庫(kù)1～10 d內(nèi)水溫變化。

圖5 不同預(yù)測(cè)時(shí)段水溫預(yù)報(bào)效果對(duì)比

以10 d預(yù)報(bào)期為例，在整個(gè)模擬周期內(nèi)選擇湖庫(kù)水溫分層發(fā)展、穩(wěn)定和減弱的三個(gè)典型時(shí)段進(jìn)一步分析，對(duì)應(yīng)時(shí)段分別5月10—19日、8月2—11日和10月18—27日，預(yù)報(bào)結(jié)果如圖6所示，圖6中虛線為95%預(yù)報(bào)置信區(qū)間。由圖6(a)可以看出，熱分層形成階段受氣象條件的影響，庫(kù)區(qū)水溫整體緩慢上升、表底溫差逐漸增大，1與13 m的溫差由4.3上升至8.6 ℃。沿水深方向來(lái)看，1 m處的預(yù)報(bào)誤差始終保持在1.3 ℃以內(nèi)，7和13 m處預(yù)報(bào)誤差也始終保持在0.8和0.3 ℃以內(nèi)。7和13 m處水溫受氣象條件影響較小，受入流條件和水庫(kù)取水等的影響較大，5月18日預(yù)報(bào)得到的水溫上升現(xiàn)象即由水庫(kù)取水所致，這也表明同化預(yù)報(bào)系統(tǒng)能夠較好的預(yù)報(bào)湖庫(kù)調(diào)度對(duì)水庫(kù)水溫的影響。由圖6(b)可以看出，此時(shí)水庫(kù)熱分層已完全形成，入庫(kù)水溫和氣溫均維持在全年的較高水平，整體而言7和13 m的預(yù)報(bào)精度依然高于1m。數(shù)據(jù)同化系統(tǒng)預(yù)報(bào)出了8月6日受短期強(qiáng)降雨和降溫引起的水庫(kù)在1和7 m處1.6～2.3 ℃的降溫過(guò)程，13 m處水溫波動(dòng)不大。由圖6(c)可以看出，水庫(kù)水溫整體逐漸降低且1和7 m處水溫已基本相等，此時(shí)較弱的熱分層使10月24日時(shí)在大風(fēng)降溫驅(qū)動(dòng)下，水庫(kù)可能產(chǎn)生較強(qiáng)的垂向摻混，發(fā)生“翻庫(kù)”等危害供水安全事件的風(fēng)險(xiǎn)較大，在運(yùn)行管理中值得關(guān)注。整體而言，數(shù)據(jù)同化系統(tǒng)在受氣象條件及水庫(kù)調(diào)度等內(nèi)外部因素影響下，能夠在10 d的預(yù)報(bào)期內(nèi)維持較高準(zhǔn)確性，為湖庫(kù)調(diào)度管理提供依據(jù)。

圖6 2019年3個(gè)典型時(shí)段未來(lái)10 d水庫(kù)水溫預(yù)報(bào)結(jié)果

5 討論

5.1 觀測(cè)數(shù)據(jù)密度的影響觀測(cè)數(shù)據(jù)在數(shù)據(jù)同化過(guò)程中扮演著關(guān)鍵作用，觀測(cè)數(shù)據(jù)的密度與同化模型的預(yù)報(bào)精度息息相關(guān)。但較高的觀測(cè)頻次往往需要投入更多的設(shè)備成本，此外由于設(shè)備、通訊等軟硬件限制，并非所有觀測(cè)點(diǎn)都能保證有時(shí)間序列上的完整數(shù)據(jù)，因此有必要分析觀測(cè)數(shù)據(jù)密度對(duì)數(shù)據(jù)同化效果的影響。在對(duì)本文預(yù)報(bào)期內(nèi)9個(gè)時(shí)段進(jìn)行10 d的同化預(yù)測(cè)中，分別設(shè)置3個(gè)深度的觀測(cè)數(shù)據(jù)密度為12、24和36 h/個(gè)，圖7為不同深度水溫預(yù)報(bào)結(jié)果(RMSE)。

圖7 不同觀測(cè)數(shù)據(jù)密度下各深度預(yù)測(cè)RMSE結(jié)果對(duì)比

從圖7可以看出，在本文的數(shù)據(jù)同化模型中，當(dāng)預(yù)報(bào)點(diǎn)位與W2模型模擬時(shí)間一致的條件下，觀測(cè)數(shù)據(jù)密度從12 h/個(gè)降低至36 h/個(gè)時(shí)，水溫預(yù)報(bào)的精度逐漸降低。但觀測(cè)數(shù)據(jù)密度并非越高越好，崔凱鵬等[40]的研究表明，過(guò)高的觀測(cè)時(shí)間密度可能會(huì)產(chǎn)生信息冗余，Thomas等[14]的研究也發(fā)現(xiàn)由于設(shè)備故障產(chǎn)生的數(shù)據(jù)損失并未對(duì)同化表現(xiàn)產(chǎn)生持續(xù)影響。因此在兼顧效率和模型預(yù)測(cè)準(zhǔn)確性的同時(shí)，確定合適的觀測(cè)數(shù)據(jù)密度，并在模擬期內(nèi)明確哪些時(shí)段的觀測(cè)數(shù)據(jù)能帶來(lái)最大價(jià)值，對(duì)選取觀測(cè)資料有著重要意義。

5.2 同化方案設(shè)定的影響基于EnKF算法的同化方案設(shè)定中，集合數(shù)(N)、觀測(cè)誤差、模擬誤差等的選取對(duì)模型精度有較大影響，對(duì)于上述指標(biāo)選取均需要進(jìn)行大量的模擬測(cè)試。集合數(shù)方面，Houtekamer等[41]的研究表明N在100這一量級(jí)對(duì)于數(shù)據(jù)同化研究是足夠的；Evensen[42]的研究表明模型誤差與N-1/2相關(guān)，當(dāng)該值更接近于0后會(huì)有利于系統(tǒng)的誤差控制。相關(guān)研究中N的取值通常在20～500之間[14，23，33]，本文根據(jù)調(diào)整集合數(shù)后模型的RMSE表現(xiàn)確定N取值為100，結(jié)果表明該值能夠較好地兼顧同化系統(tǒng)的模擬精度與計(jì)算效率。觀測(cè)誤差和模擬誤差也會(huì)對(duì)模型預(yù)報(bào)精度有較大的影響，通常認(rèn)為觀測(cè)數(shù)據(jù)精度較高且更接近真實(shí)值，因此觀測(cè)誤差取值往往小于模擬誤差[23，43]。從本文16組誤差組合下的模擬結(jié)果看，當(dāng)模擬誤差一定時(shí)，RMSE隨觀測(cè)誤差的增大而增大；當(dāng)觀測(cè)誤差一定時(shí)，RMSE隨模擬誤差的增大而增大的趨勢(shì)并不明顯，劉卓等[23]認(rèn)為觀測(cè)誤差對(duì)同化精度的影響較模擬誤差大是產(chǎn)生這一現(xiàn)象的主要原因，因此設(shè)置較小的觀測(cè)誤差是合理的。

6 結(jié)論與展望

本文基于集合卡爾曼濾波算法與CE-QUAL-W2模型，構(gòu)建可綜合考慮模型參數(shù)、邊界條件以及觀測(cè)數(shù)據(jù)不確定性的湖庫(kù)水溫?cái)?shù)據(jù)同化系統(tǒng)，將其應(yīng)用于大黑汀水庫(kù)進(jìn)行1～10 d的中短期水溫預(yù)報(bào)，主要結(jié)論包括：(1)同化方案設(shè)計(jì)對(duì)數(shù)據(jù)同化系統(tǒng)具有重要影響，當(dāng)集合數(shù)為100、模擬誤差和觀測(cè)誤差分別為10%和1%時(shí)，本文構(gòu)建的同化系統(tǒng)能夠兼顧較高的計(jì)算效率與模擬精度。(2)采用雙重假設(shè)方法，經(jīng)過(guò)兩次集合卡爾曼濾波算法對(duì)W2模型的水溫敏感參數(shù)和水溫狀態(tài)變量進(jìn)行同步更新后，數(shù)據(jù)同化系統(tǒng)在1、7和13 m水深處的水溫模擬準(zhǔn)確性較無(wú)數(shù)據(jù)同化的模擬結(jié)果分別提升了68.8%、51.6%和41.2%，模擬精度顯著提升。(3)利用數(shù)據(jù)同化系統(tǒng)進(jìn)行水溫預(yù)報(bào)時(shí)，隨著預(yù)報(bào)期由1 d延長(zhǎng)至10 d，不同水深的預(yù)報(bào)誤差由0.22～0.35 ℃增大至0.77～1.09 ℃。無(wú)論水庫(kù)處于分層期或混合期，數(shù)據(jù)同化系統(tǒng)均能夠在預(yù)報(bào)期內(nèi)的氣象條件及水庫(kù)調(diào)度等內(nèi)外部因素驅(qū)動(dòng)下維持較高準(zhǔn)確性，為水庫(kù)中短期調(diào)度管理提供理論支撐。此外本文是針對(duì)W2模型開(kāi)展數(shù)據(jù)同化研究的首個(gè)案例，所構(gòu)建的同化系統(tǒng)對(duì)于其他應(yīng)用W2模型開(kāi)展水動(dòng)力水質(zhì)模擬與預(yù)報(bào)的研究具有參考價(jià)值。

本項(xiàng)研究還能在以下四方面進(jìn)行提升：(1)不確定性分析方面。本文的預(yù)報(bào)效果除了與模型參數(shù)和邊界條件有關(guān)，還會(huì)受觀測(cè)數(shù)據(jù)密度與同化方案的影響，因此系統(tǒng)分析上述要素對(duì)預(yù)報(bào)結(jié)果的不確定性很有必要。此外本文只選取了3個(gè)對(duì)W2模型水溫模擬最敏感的參數(shù)進(jìn)行數(shù)據(jù)同化，這也可能造成一定的不確定性低估。(2)觀測(cè)數(shù)據(jù)獲取方面。本文采用同一點(diǎn)位、不同深度的原位觀測(cè)數(shù)據(jù)進(jìn)行同化系統(tǒng)構(gòu)建，但在對(duì)空間異質(zhì)性較強(qiáng)的湖庫(kù)進(jìn)行預(yù)報(bào)時(shí)，增大觀測(cè)數(shù)據(jù)的空間密度十分必要。布設(shè)多個(gè)觀測(cè)站點(diǎn)，提升觀測(cè)數(shù)據(jù)的空間密度較容易實(shí)現(xiàn)，但也意味著需要投入成倍的觀測(cè)成本。遙感數(shù)據(jù)具有較強(qiáng)的空間代表性，利用該方法獲取觀測(cè)數(shù)據(jù)已有研究在開(kāi)展[33，36]，多源觀測(cè)數(shù)據(jù)的融合將會(huì)進(jìn)一步提升數(shù)據(jù)同化系統(tǒng)的預(yù)報(bào)性能。(3)從信息融合角度而言，數(shù)據(jù)同化與機(jī)器學(xué)習(xí)在數(shù)學(xué)原理上采用的均為貝葉斯估計(jì)策略[44]，融合機(jī)器學(xué)習(xí)和與機(jī)理模型的水溫預(yù)測(cè)研究已經(jīng)逐漸開(kāi)展[45]，數(shù)據(jù)同化方法可能再次扮演重要角色。(4)數(shù)字孿生是當(dāng)前水利工程領(lǐng)域的重要研究方向，建設(shè)數(shù)字孿生水利工程就是為了服務(wù)預(yù)報(bào)、預(yù)警、預(yù)演、預(yù)案四項(xiàng)基本功能，其中預(yù)報(bào)功能要實(shí)現(xiàn)對(duì)水安全要素的高精度短期預(yù)報(bào)和中期預(yù)測(cè)，本文也可為開(kāi)展數(shù)字孿生水利工程的中短期預(yù)報(bào)研究提供思路和方法。