孫博聞，宗慶志，楊晰淯，張?jiān)瑢?，高學(xué)平

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

1 研究背景

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

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

數(shù)學(xué)模型是研究深水水庫水動(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)行大量的代碼改寫，這也是造成當(dāng)前大多數(shù)數(shù)學(xué)模型并未廣泛采用GPU并行計(jì)算的重要原因。在提升模型計(jì)算精度方面，有學(xué)者提出利用數(shù)據(jù)同化的方法將觀測數(shù)據(jù)與模型計(jì)算結(jié)果相融合，結(jié)合現(xiàn)場觀測和數(shù)值模擬的優(yōu)點(diǎn)，實(shí)現(xiàn)模型狀態(tài)和參數(shù)不斷更新，以有效提高模型的模擬精度[19]。但受限于初始條件、邊界條件等不確定性的影響，目前數(shù)據(jù)同化方法在深水水庫中短期水質(zhì)預(yù)報(bào)的應(yīng)用依然較少。

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

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

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

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

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

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

2.4 基于OpenDA的數(shù)據(jù)同化與預(yù)報(bào)OpenDA(Open data assimilation)是一個(gè)可快速實(shí)現(xiàn)參數(shù)校準(zhǔn)和濾波測試等功能的數(shù)據(jù)同化工具箱。它提供了一個(gè)允許算法和模型進(jìn)行簡單交換的平臺(tái)，通過將數(shù)據(jù)同化算法作為單獨(dú)組件，降低數(shù)據(jù)同化的應(yīng)用成本，近年來廣泛應(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ù)測模型系統(tǒng)與觀測系統(tǒng)，同化算法選用集合卡爾曼濾波算法(EnKF，Ensemble Kalman filter algorithm)。在同化階段預(yù)測模型都要運(yùn)行集合數(shù)次，當(dāng)有觀測數(shù)據(jù)時(shí)，通過EnKF算法進(jìn)行數(shù)據(jù)同化。數(shù)據(jù)同化后的結(jié)果被稱為分析值集合，該集合的均值被認(rèn)為是最優(yōu)估計(jì)值，分析值集合將作為下個(gè)時(shí)間步預(yù)報(bào)的初始值，直到下一個(gè)有觀測數(shù)據(jù)的時(shí)刻進(jìn)行數(shù)據(jù)同化。

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

3 結(jié)果與分析

3.1 GPU加速效果

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

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

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)語組合二作為本實(shí)驗(yàn)的OpenACC導(dǎo)語形式，實(shí)驗(yàn)結(jié)果如圖3(b)所示。在計(jì)算網(wǎng)格數(shù)小于600時(shí)，CPU單精度的計(jì)算模式耗時(shí)最少，但隨著網(wǎng)格數(shù)的增加，GPU并行計(jì)算耗時(shí)低短的優(yōu)勢開始逐漸顯現(xiàn)出來。而當(dāng)計(jì)算網(wǎng)格數(shù)增加到950時(shí)，雖然仍存在數(shù)據(jù)傳輸?shù)暮臅r(shí)問題，但由于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)語組合或網(wǎng)格數(shù)量條件下的計(jì)算用時(shí)對比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) 采用大黑汀水庫壩前斷面水質(zhì)監(jiān)測站點(diǎn)2019年實(shí)測5個(gè)月的溫度和溶解氧數(shù)據(jù)(采樣間隔30 min)，對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用來表征模擬值和觀測值的一致性，取值在0～1之間，越接近1表示二者之間的一致程度越高，PBIAS結(jié)果以百分比定量化給出模擬值比觀測值整體被低估或高估的平均趨勢，PBIAS為0表示模擬效果最佳，PBIAS>0表示模擬值傾向于低估，反之PBIAS<0表示模型傾向于高估。

(1)

(2)

(3)

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

圖4 W2模型模擬結(jié)果與數(shù)據(jù)同化系統(tǒng)模擬結(jié)果相比實(shí)測水溫和溶解氧的偏差值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)行水庫水溫預(yù)報(bào)。參考Zhang等[28]的研究結(jié)果，以單日內(nèi)是否所有剖面數(shù)據(jù)都監(jiān)測出溫躍層為依據(jù)判斷分層狀態(tài)，即5月27日—10月4日為分層穩(wěn)定期，此前為分層形成期，此后為分層削弱期。為分析數(shù)據(jù)同化系統(tǒng)的水溫預(yù)報(bào)性能，繪制分層期、混合期與全時(shí)段下各水深處未來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)來看，隨著預(yù)報(bào)時(shí)段由1天延長至10天，預(yù)報(bào)誤差由0.22～0.35 ℃增大至0.77～1.09 ℃，即預(yù)報(bào)精度會(huì)隨著預(yù)報(bào)時(shí)段延長逐漸降低。其他學(xué)者的研究結(jié)果也表明，預(yù)報(bào)時(shí)長與預(yù)報(bào)精度顯著相關(guān)[29-31]，預(yù)報(bào)時(shí)長的增加會(huì)使預(yù)報(bào)結(jié)果的不確定性增大。從IOA指標(biāo)來看，分層期7天預(yù)報(bào)期內(nèi)各深度的IOA值均超過0.65，具有較好的一致性；預(yù)報(bào)期延長至10天后，7 m處的IOA由0.78降至0.57，降幅遠(yuǎn)超1 m與13 m深度。從PBIAS指標(biāo)來看，1 m和13 m處預(yù)報(bào)均傾向于高估，7 m處則傾向于低估，但高/低估均在±3%以內(nèi)。近年來引灤工程常采用潘家口水庫大流量調(diào)度抑制大黑汀水庫缺氧發(fā)展[19-20，32]，水庫調(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)水庫1～10天內(nèi)水溫變化。

圖5 分層期、混合期與全時(shí)段各水深處未來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)對溶解氧的預(yù)報(bào)性能，在水庫分層狀態(tài)和混合狀態(tài)下對壩前斷面各水深處未來1～10天的溶解氧預(yù)報(bào)結(jié)果進(jìn)行評(píng)價(jià)，結(jié)果如圖6。隨著預(yù)報(bào)時(shí)段由1天延長至10天，RMSE由0.06～0.25 mg/L增大至0.45～1.09 mg/L，可見隨著預(yù)報(bào)時(shí)段延長，不確定因素的累計(jì)效應(yīng)對預(yù)報(bào)會(huì)產(chǎn)生一定影響。分層期7天預(yù)報(bào)期內(nèi)各深度的IOA值均超過0.65，具有較好的一致性；預(yù)報(bào)期延長至10天后，13 m處的IOA由0.8降至0.64，降幅遠(yuǎn)超1 m與7 m深度。從PBIAS指標(biāo)來看，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)水庫1～10天內(nèi)溶解氧變化。

圖6 各水深處未來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模型建立水庫水動(dòng)力水質(zhì)數(shù)學(xué)模型，基于OpenACC的GPU并行方法提升模型計(jì)算效率，構(gòu)建大黑汀水庫水溫與溶解氧的數(shù)據(jù)同化系統(tǒng)，在中短期時(shí)間尺度上實(shí)現(xiàn)對水庫水溫與溶解氧動(dòng)態(tài)變化的高精度、高效率預(yù)報(bào)，主要結(jié)論包括：(1)GPU并行計(jì)算中，不同導(dǎo)語組合的計(jì)算耗時(shí)差異明顯，網(wǎng)格數(shù)量為1000時(shí)，計(jì)算耗時(shí)相差可達(dá)83 min。當(dāng)OpenACC采用kernels導(dǎo)語加copy和copyin子語的組合時(shí)，GPU并行加速的表現(xiàn)最好。當(dāng)網(wǎng)格數(shù)超過600后，網(wǎng)格數(shù)量越多，GPU并行計(jì)算的效率優(yōu)勢越大。(2)采用EnKF方法進(jìn)行數(shù)據(jù)同化，同時(shí)同化狀態(tài)變量與模型參數(shù)，不同深度處水溫模擬結(jié)果較無同化條件提升了41.2%～68.8%，溶解氧模擬結(jié)果提升了54.4%～71.2%。有助于提升對水溫和溶解氧的中短期預(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é)果符合水庫水溫與溶解氧的中短期變化規(guī)律，能夠?yàn)榇蠛谕∷畮斓墓┧c生態(tài)安全提供支撐。