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

孫博聞，宗慶志，楊晰淯，張?jiān)瑢帲邔W(xué)平

(天津大學(xué) 水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300350)

1 研究背景

深水湖庫(kù)通常存在季節(jié)性溫度分層現(xiàn)象，水溫影響著水體中生化反應(yīng)[1-5]，由溫度分層引起溶解氧分層與季節(jié)性缺氧還會(huì)誘發(fā)水環(huán)境水生態(tài)問(wèn)題[6-8]，進(jìn)而威脅供水安全及下游生態(tài)系統(tǒng)健康。在氣候變化與水庫(kù)調(diào)控的共同影響下，上述變化更加趨于復(fù)雜[9-10]，因此研究水溫和溶解氧的變化對(duì)湖庫(kù)生態(tài)環(huán)境安全十分必要。

當(dāng)前對(duì)水庫(kù)水溫和溶解氧的研究大多關(guān)注其長(zhǎng)期演變規(guī)律[11-13]，但在未來(lái)1～10天的中短期時(shí)間尺度上預(yù)報(bào)湖庫(kù)水溫與溶解氧變化對(duì)管理者同樣具有重要意義[14]。如Weber等[15]利用水溫-深度變化數(shù)據(jù)動(dòng)態(tài)確定Grosse Dhuenn水庫(kù)取水深度，在優(yōu)化下泄水溫同時(shí)還抑制了水庫(kù)缺氧的發(fā)展；強(qiáng)風(fēng)驅(qū)動(dòng)產(chǎn)生的劇烈涌升流使日內(nèi)瓦湖部分湖區(qū)的上層水溫驟降，直接影響湖區(qū)溶解氧分布與下游供水[16]；余曉等[17]對(duì)潘家口水庫(kù)全年的水溫和溶解氧濃度進(jìn)行時(shí)空分布特征和演化成因的研究，得出水溫分層控制溶解氧分層的結(jié)論；劉暢等[18]發(fā)現(xiàn)大流量調(diào)度過(guò)程引起的強(qiáng)對(duì)流條件也可以產(chǎn)生短期削弱甚至完全消除由水庫(kù)分層導(dǎo)致的缺氧現(xiàn)象。可見(jiàn)對(duì)湖庫(kù)水溫與溶解氧在中短期時(shí)間尺度進(jìn)行預(yù)報(bào)，能夠在調(diào)控湖庫(kù)下泄水溫這類常規(guī)管理中發(fā)揮作用，還能在應(yīng)對(duì)極端天氣等威脅供水與生態(tài)安全的應(yīng)急管理工作中提供決策支撐。

數(shù)學(xué)模型是研究深水水庫(kù)水動(dòng)力水質(zhì)演化規(guī)律的重要方法，但受制于模型計(jì)算速度以及模型初始條件與參數(shù)取值等原因，模型在中短期時(shí)間尺度上的預(yù)報(bào)效率與精度還有較大的提升空間。在提升模型計(jì)算效率方面，并行計(jì)算是當(dāng)前主流解決方案。如能夠采用合適的優(yōu)化方法利用CPU(Central Processing Unit)處理數(shù)學(xué)模型邏輯運(yùn)算部分，利用GPU(Graphic Processing Unit)執(zhí)行并行運(yùn)算部分，有可能大幅縮減模型運(yùn)行時(shí)間。但利用GPU運(yùn)算通常需要進(jìn)行大量的代碼改寫(xiě)，這也是造成當(dāng)前大多數(shù)數(shù)學(xué)模型并未廣泛采用GPU并行計(jì)算的重要原因。在提升模型計(jì)算精度方面，有學(xué)者提出利用數(shù)據(jù)同化的方法將觀測(cè)數(shù)據(jù)與模型計(jì)算結(jié)果相融合，結(jié)合現(xiàn)場(chǎng)觀測(cè)和數(shù)值模擬的優(yōu)點(diǎn)，實(shí)現(xiàn)模型狀態(tài)和參數(shù)不斷更新，以有效提高模型的模擬精度[19]。但受限于初始條件、邊界條件等不確定性的影響，目前數(shù)據(jù)同化方法在深水水庫(kù)中短期水質(zhì)預(yù)報(bào)的應(yīng)用依然較少。

本文以引灤入津工程—大黑汀水庫(kù)為研究對(duì)象，利用水庫(kù)長(zhǎng)時(shí)段、高頻次自動(dòng)監(jiān)測(cè)數(shù)據(jù)為同化數(shù)據(jù)源，采用集合卡爾曼濾波算法作為同化方法，基于CE-QUAL-W2模型(the two-dimensional，laterally averaged，hydrodynamic and water quality model)建立水庫(kù)水動(dòng)力水質(zhì)數(shù)學(xué)模型，選用GPU并行計(jì)算方法提高模型計(jì)算效率，構(gòu)建大黑汀水庫(kù)水溫與溶解氧的數(shù)據(jù)同化系統(tǒng)，以期實(shí)現(xiàn)對(duì)水庫(kù)垂向水溫與溶解氧的動(dòng)態(tài)變化的高精度、高效率的中短期預(yù)報(bào)，為水庫(kù)水環(huán)境管理提供可靠的技術(shù)支撐。

2 深水湖庫(kù)溫-氧中短期預(yù)報(bào)計(jì)算

2.1 研究區(qū)域大黑汀水庫(kù)位于唐山市遷西縣城北5 km的灤河干流，地處北溫帶，四季分明，控制流域面積35 100 km2，總庫(kù)容3.37億m3。大黑汀水庫(kù)是引灤入津樞紐工程的源頭水庫(kù)，承擔(dān)著向唐山及天津供水的重要任務(wù)。自1990年代以來(lái)，大黑汀水庫(kù)開(kāi)始出現(xiàn)大規(guī)模的網(wǎng)箱養(yǎng)殖活動(dòng)(圖1)，人工投入的養(yǎng)殖飼料、魚(yú)類排泄物等外源有機(jī)污染持續(xù)被投入，對(duì)水質(zhì)造成了嚴(yán)重污染。水庫(kù)管理局自2016年開(kāi)始大面積清理養(yǎng)殖網(wǎng)箱，污染負(fù)荷被清除。但幾十年養(yǎng)殖活動(dòng)產(chǎn)生的污染物依舊沉積在水庫(kù)底泥當(dāng)中。對(duì)水庫(kù)水溫和溶解氧的監(jiān)測(cè)數(shù)據(jù)表明，庫(kù)區(qū)底部存在明顯的缺氧區(qū)，水環(huán)境質(zhì)量明顯下降。

圖1 大黑汀水庫(kù)概況圖Fig.1 General situation of Daheiting Reservoir

2.2 CE-QUAL-W2模型本文采用CE-QUAL-W2構(gòu)建大黑汀水庫(kù)水動(dòng)力水質(zhì)數(shù)學(xué)模型。CE-QUAL-W2是一個(gè)考慮縱向和垂向的二維水動(dòng)力水質(zhì)模型(簡(jiǎn)稱W2模型)，模型假設(shè)水體橫向平均，因此適宜模擬河道型水庫(kù)的水動(dòng)力與水質(zhì)變化。模型邊界條件包括上游流量、水溫及水質(zhì)、下游水位等，氣象條件數(shù)據(jù)包括氣溫、露點(diǎn)溫度、風(fēng)速、風(fēng)向和云量等。模型建立和參數(shù)率定過(guò)程參考姚嘉偉研究成果[20]，本文不再贅述。

2.3 基于OpenACC的GPU并行計(jì)算OpenACC(Open Accelerators)是一種新型的GPU并行方式，它利用導(dǎo)語(yǔ)向編譯器提供編譯和優(yōu)化代碼的相關(guān)信息，由于其執(zhí)行效率與實(shí)現(xiàn)方法綜合較優(yōu)，近年來(lái)開(kāi)始應(yīng)用于水動(dòng)力數(shù)學(xué)模型的GPU并行計(jì)算[21-23]。在對(duì)程序進(jìn)行并行計(jì)算時(shí)，由于外部邏輯循環(huán)包括前后處理及串行代碼很難實(shí)現(xiàn)在GPU上的并行計(jì)算，所以選擇適當(dāng)?shù)难h(huán)代碼進(jìn)行并行計(jì)算有利于提升程序并行加速效果。其次，選出循環(huán)復(fù)雜的代碼部分，在代碼的開(kāi)始和結(jié)尾處添加OpenACC導(dǎo)語(yǔ)代碼，對(duì)程序進(jìn)行編譯，查看編譯器的反饋信息并記錄初步的計(jì)算用時(shí)。此時(shí)有可能會(huì)發(fā)生計(jì)算用時(shí)變長(zhǎng)的情況，這是因?yàn)殚_(kāi)始GPU并行計(jì)算后，增加了GPU和CPU之間的數(shù)據(jù)傳輸?shù)暮臅r(shí)，如果傳輸?shù)臄?shù)據(jù)量過(guò)大，就會(huì)增加計(jì)算耗時(shí)，進(jìn)而影響程序的計(jì)算性能。本研究基于OpenACC采用CPU-GPU協(xié)同計(jì)算結(jié)構(gòu)實(shí)現(xiàn)W2模型的并行計(jì)算(圖2(a))，通過(guò)OpenACC計(jì)算導(dǎo)語(yǔ)將W2模型水動(dòng)力計(jì)算中較復(fù)雜的部分分配給GPU處理器，CPU處理器完成剩余計(jì)算任務(wù)及邏輯控制，選用Tesla K40c處理器作為GPU加速設(shè)備。由于數(shù)據(jù)同化過(guò)程需要執(zhí)行大量模型計(jì)算，因此通過(guò)比較不同計(jì)算導(dǎo)語(yǔ)和模型網(wǎng)格數(shù)量對(duì)GPU加速效果的影響，確定數(shù)據(jù)同化階段采用的GPU并行計(jì)算方案。

2.4 基于OpenDA的數(shù)據(jù)同化與預(yù)報(bào)OpenDA(Open data assimilation)是一個(gè)可快速實(shí)現(xiàn)參數(shù)校準(zhǔn)和濾波測(cè)試等功能的數(shù)據(jù)同化工具箱。它提供了一個(gè)允許算法和模型進(jìn)行簡(jiǎn)單交換的平臺(tái)，通過(guò)將數(shù)據(jù)同化算法作為單獨(dú)組件，降低數(shù)據(jù)同化的應(yīng)用成本，近年來(lái)廣泛應(yīng)用于水文模型與水動(dòng)力模型的數(shù)據(jù)同化領(lǐng)域[24-27]。本研究基于OpenDA進(jìn)行水溫與溶解氧的數(shù)據(jù)同化和預(yù)報(bào)，流程如圖2(b)所示。其中數(shù)據(jù)同化系統(tǒng)包括預(yù)測(cè)模型系統(tǒng)與觀測(cè)系統(tǒng)，同化算法選用集合卡爾曼濾波算法(EnKF，Ensemble Kalman filter algorithm)。在同化階段預(yù)測(cè)模型都要運(yùn)行集合數(shù)次，當(dāng)有觀測(cè)數(shù)據(jù)時(shí)，通過(guò)EnKF算法進(jìn)行數(shù)據(jù)同化。數(shù)據(jù)同化后的結(jié)果被稱為分析值集合，該集合的均值被認(rèn)為是最優(yōu)估計(jì)值，分析值集合將作為下個(gè)時(shí)間步預(yù)報(bào)的初始值，直到下一個(gè)有觀測(cè)數(shù)據(jù)的時(shí)刻進(jìn)行數(shù)據(jù)同化。

圖2 計(jì)算流程Fig.2 Computing process

3 結(jié)果與分析

3.1 GPU加速效果

3.1.1 不同導(dǎo)語(yǔ)選擇對(duì)GPU優(yōu)化效果影響 表1列出了兩種導(dǎo)語(yǔ)組合格式，相應(yīng)計(jì)算用時(shí)對(duì)比如圖3(a)所示?？梢钥闯鲈贕PU并行計(jì)算后，由于增加了主機(jī)和GPU之間的數(shù)據(jù)傳輸?shù)暮臅r(shí)，因此不同導(dǎo)語(yǔ)組合方式的計(jì)算時(shí)間有較明顯的差異。在不同計(jì)算網(wǎng)格數(shù)下，導(dǎo)語(yǔ)組合二的計(jì)算耗時(shí)始終小于導(dǎo)語(yǔ)組合一。隨著網(wǎng)格數(shù)的增加，兩組導(dǎo)語(yǔ)組合的計(jì)算耗時(shí)差逐漸增大，網(wǎng)格數(shù)為1000時(shí)耗時(shí)相差最大，達(dá)到了83 min。因?yàn)閷?dǎo)語(yǔ)組合二在導(dǎo)語(yǔ)組合一的基礎(chǔ)上，增加了copy子語(yǔ)和copyin子語(yǔ)，其中copy子語(yǔ)確保了數(shù)據(jù)由主機(jī)到GPU之間的雙向傳輸?shù)目尚行?，而copyin子語(yǔ)則可以將目標(biāo)計(jì)算數(shù)據(jù)復(fù)制到顯卡內(nèi)存中，并且在GPU完成第一部分并行計(jì)算后不把數(shù)據(jù)傳輸回主機(jī)，所以可以有效的減少數(shù)據(jù)傳輸消耗時(shí)間，提升GPU并行效率。因此選擇適當(dāng)?shù)腛penACC導(dǎo)語(yǔ)，可以有效減少數(shù)據(jù)傳輸所消耗的時(shí)間，提升GPU并行加速效果。

表1 OpenACC的不同導(dǎo)語(yǔ)組合

3.1.2 不同網(wǎng)格數(shù)量下GPU加速效果分析 本研究網(wǎng)格數(shù)設(shè)定為300、350、400、450、500、550、600、650、700、750、800、850、900、950和1000。選取3.1.1節(jié)中GPU并行效率較高的導(dǎo)語(yǔ)組合二作為本實(shí)驗(yàn)的OpenACC導(dǎo)語(yǔ)形式，實(shí)驗(yàn)結(jié)果如圖3(b)所示。在計(jì)算網(wǎng)格數(shù)小于600時(shí)，CPU單精度的計(jì)算模式耗時(shí)最少，但隨著網(wǎng)格數(shù)的增加，GPU并行計(jì)算耗時(shí)低短的優(yōu)勢(shì)開(kāi)始逐漸顯現(xiàn)出來(lái)。而當(dāng)計(jì)算網(wǎng)格數(shù)增加到950時(shí)，雖然仍存在數(shù)據(jù)傳輸?shù)暮臅r(shí)問(wèn)題，但由于GPU強(qiáng)大的并行計(jì)算能力，CPU-GPU協(xié)同計(jì)算比CPU串行計(jì)算節(jié)省用時(shí)可達(dá)46 min。為保證GPU加速效果，同時(shí)考慮模擬精度要求，最終確定模型縱向劃分為69個(gè)單元(Segment)，垂向最多33層(Layer)。

圖3 不同導(dǎo)語(yǔ)組合或網(wǎng)格數(shù)量條件下的計(jì)算用時(shí)對(duì)比Fig.3 Comparison of computational time under different combinations of directives or grid quantity

3.2 水溫及溶解氧數(shù)據(jù)同化

3.2.1 同化和預(yù)報(bào)性能判別標(biāo)準(zhǔn) 采用大黑汀水庫(kù)壩前斷面水質(zhì)監(jiān)測(cè)站點(diǎn)2019年實(shí)測(cè)5個(gè)月的溫度和溶解氧數(shù)據(jù)(采樣間隔30 min)，對(duì)CE-QUAL-W2模型的狀態(tài)變量和參數(shù)同時(shí)進(jìn)行同化。以均方根誤差(RMSE)、一致性系數(shù)(IOA)和百分比偏差系數(shù)(PBIAS)作為評(píng)價(jià)標(biāo)準(zhǔn)(式(1)—式(3))，RMSE越接近0表示模擬效果越佳，IOA用來(lái)表征模擬值和觀測(cè)值的一致性，取值在0～1之間，越接近1表示二者之間的一致程度越高，PBIAS結(jié)果以百分比定量化給出模擬值比觀測(cè)值整體被低估或高估的平均趨勢(shì)，PBIAS為0表示模擬效果最佳，PBIAS>0表示模擬值傾向于低估，反之PBIAS<0表示模型傾向于高估。

(1)

(2)

(3)

3.2.2 同化效果分析 圖4為有無(wú)數(shù)據(jù)同化的情況下，壩前斷面水下1 m、7 m和13 m處的模擬水溫與實(shí)測(cè)水溫、模擬溶解氧和實(shí)測(cè)溶解氧間的偏差值。可以看出，在經(jīng)過(guò)EnKF方法對(duì)模型參數(shù)和狀態(tài)變量同時(shí)進(jìn)行更新后，模型的模擬精度有了很大的提升，數(shù)據(jù)同化模型在壩前斷面水下1 m、7 m和13 m處的水溫平均模擬偏差分別為1.01 ℃、0.31 ℃和0.21 ℃，較無(wú)數(shù)據(jù)同化的模擬結(jié)果分別提升了68.8%、51.6%和41.2%；溶解氧平均模擬偏差分別為0.38 mg/L、0.22 mg/L和0.18 mg/L，較無(wú)數(shù)據(jù)同化的模擬結(jié)果分別提升了71.2%、63.6%和54.4%。

圖4 W2模型模擬結(jié)果與數(shù)據(jù)同化系統(tǒng)模擬結(jié)果相比實(shí)測(cè)水溫和溶解氧的偏差值Fig.4 The water temperature errors and the dissolved oxygen (DO) errors

3.3 水溫及溶解氧中短期預(yù)報(bào)

3.3.1 水溫中短期預(yù)報(bào) 以數(shù)據(jù)同化系統(tǒng)為基礎(chǔ)，考慮邊界條件等不確定性前提下，運(yùn)行W2模型進(jìn)行水庫(kù)水溫預(yù)報(bào)。參考Zhang等[28]的研究結(jié)果，以單日內(nèi)是否所有剖面數(shù)據(jù)都監(jiān)測(cè)出溫躍層為依據(jù)判斷分層狀態(tài)，即5月27日—10月4日為分層穩(wěn)定期，此前為分層形成期，此后為分層削弱期。為分析數(shù)據(jù)同化系統(tǒng)的水溫預(yù)報(bào)性能，繪制分層期、混合期與全時(shí)段下各水深處未來(lái)1～10天的水溫預(yù)報(bào)結(jié)果，如圖5(第1—3列)所示?？梢钥闯?，在1 m與13 m處分層期的預(yù)報(bào)性能優(yōu)于混合期，7 m處兩階段預(yù)報(bào)性能各有優(yōu)劣。從RMSE指標(biāo)來(lái)看，隨著預(yù)報(bào)時(shí)段由1天延長(zhǎng)至10天，預(yù)報(bào)誤差由0.22～0.35 ℃增大至0.77～1.09 ℃，即預(yù)報(bào)精度會(huì)隨著預(yù)報(bào)時(shí)段延長(zhǎng)逐漸降低。其他學(xué)者的研究結(jié)果也表明，預(yù)報(bào)時(shí)長(zhǎng)與預(yù)報(bào)精度顯著相關(guān)[29-31]，預(yù)報(bào)時(shí)長(zhǎng)的增加會(huì)使預(yù)報(bào)結(jié)果的不確定性增大。從IOA指標(biāo)來(lái)看，分層期7天預(yù)報(bào)期內(nèi)各深度的IOA值均超過(guò)0.65，具有較好的一致性；預(yù)報(bào)期延長(zhǎng)至10天后，7 m處的IOA由0.78降至0.57，降幅遠(yuǎn)超1 m與13 m深度。從PBIAS指標(biāo)來(lái)看，1 m和13 m處預(yù)報(bào)均傾向于高估，7 m處則傾向于低估，但高/低估均在±3%以內(nèi)。近年來(lái)引灤工程常采用潘家口水庫(kù)大流量調(diào)度抑制大黑汀水庫(kù)缺氧發(fā)展[19-20，32]，水庫(kù)調(diào)度提高了水溫預(yù)報(bào)的不確定性，在溫躍層附近的水溫預(yù)報(bào)精度依然存在提升空間。從各項(xiàng)統(tǒng)計(jì)指標(biāo)結(jié)果可以看出，構(gòu)建的數(shù)據(jù)同化系統(tǒng)能夠較好地預(yù)報(bào)水庫(kù)1～10天內(nèi)水溫變化。

圖5 分層期、混合期與全時(shí)段各水深處未來(lái)1～10天的水溫預(yù)報(bào)結(jié)果評(píng)價(jià)Fig.5 Evaluation of water temperature forecast results for the next 1-10 days at different depths

3.3.2 溶解氧中短期預(yù)報(bào) 為分析數(shù)據(jù)同化系統(tǒng)對(duì)溶解氧的預(yù)報(bào)性能，在水庫(kù)分層狀態(tài)和混合狀態(tài)下對(duì)壩前斷面各水深處未來(lái)1～10天的溶解氧預(yù)報(bào)結(jié)果進(jìn)行評(píng)價(jià)，結(jié)果如圖6。隨著預(yù)報(bào)時(shí)段由1天延長(zhǎng)至10天，RMSE由0.06～0.25 mg/L增大至0.45～1.09 mg/L，可見(jiàn)隨著預(yù)報(bào)時(shí)段延長(zhǎng)，不確定因素的累計(jì)效應(yīng)對(duì)預(yù)報(bào)會(huì)產(chǎn)生一定影響。分層期7天預(yù)報(bào)期內(nèi)各深度的IOA值均超過(guò)0.65，具有較好的一致性；預(yù)報(bào)期延長(zhǎng)至10天后，13 m處的IOA由0.8降至0.64，降幅遠(yuǎn)超1 m與7 m深度。從PBIAS指標(biāo)來(lái)看，1 m和13 m處預(yù)報(bào)均傾向于高估，7 m處則傾向于低估，但高/低估均在±3%以內(nèi)。從各項(xiàng)統(tǒng)計(jì)指標(biāo)結(jié)果可以看出，構(gòu)建的數(shù)據(jù)同化系統(tǒng)能夠較好的預(yù)報(bào)水庫(kù)1～10天內(nèi)溶解氧變化。

圖6 各水深處未來(lái)1～10天的溶解氧預(yù)報(bào)結(jié)果評(píng)價(jià)Fig.6 Evaluation of DO forecast results for the next 1-10 days at different depths

4 結(jié)論

本文采用集合卡爾曼濾波算法作為同化方法，基于CE-QUAL-W2模型建立水庫(kù)水動(dòng)力水質(zhì)數(shù)學(xué)模型，基于OpenACC的GPU并行方法提升模型計(jì)算效率，構(gòu)建大黑汀水庫(kù)水溫與溶解氧的數(shù)據(jù)同化系統(tǒng)，在中短期時(shí)間尺度上實(shí)現(xiàn)對(duì)水庫(kù)水溫與溶解氧動(dòng)態(tài)變化的高精度、高效率預(yù)報(bào)，主要結(jié)論包括：(1)GPU并行計(jì)算中，不同導(dǎo)語(yǔ)組合的計(jì)算耗時(shí)差異明顯，網(wǎng)格數(shù)量為1000時(shí)，計(jì)算耗時(shí)相差可達(dá)83 min。當(dāng)OpenACC采用kernels導(dǎo)語(yǔ)加copy和copyin子語(yǔ)的組合時(shí)，GPU并行加速的表現(xiàn)最好。當(dāng)網(wǎng)格數(shù)超過(guò)600后，網(wǎng)格數(shù)量越多，GPU并行計(jì)算的效率優(yōu)勢(shì)越大。(2)采用EnKF方法進(jìn)行數(shù)據(jù)同化，同時(shí)同化狀態(tài)變量與模型參數(shù)，不同深度處水溫模擬結(jié)果較無(wú)同化條件提升了41.2%～68.8%，溶解氧模擬結(jié)果提升了54.4%～71.2%。有助于提升對(duì)水溫和溶解氧的中短期預(yù)報(bào)性能。(3)1 m、7 m和13 m處水溫預(yù)報(bào)的平均偏差分別在0.60 ℃、0.32 ℃和0.28 ℃以內(nèi)；1 m、7 m和13 m處溶解氧預(yù)報(bào)平均偏差分別在0.54 mg/L、0.32 mg/L和0.24 mg/L以內(nèi)。預(yù)報(bào)結(jié)果符合水庫(kù)水溫與溶解氧的中短期變化規(guī)律，能夠?yàn)榇蠛谕∷畮?kù)的供水與生態(tài)安全提供支撐。