數(shù)字孿生丹江口水質(zhì)模型關(guān)鍵技術(shù)與應(yīng)用

摘要：數(shù)字孿生技術(shù)作為近年來(lái)的創(chuàng)新技術(shù)，通過(guò)構(gòu)建精確的虛擬模型實(shí)現(xiàn)了對(duì)實(shí)際物理世界的高度映射，為解決傳統(tǒng)水質(zhì)安全管理中的難題提供了全新的視角和方法。水質(zhì)模型是實(shí)現(xiàn)數(shù)字孿生湖庫(kù)水質(zhì)安全精準(zhǔn)映射和虛實(shí)互動(dòng)的核心。圍繞南水北調(diào)中線工程的水源地——丹江口水庫(kù)的水質(zhì)安全管理需求，針對(duì)傳統(tǒng)模型存在的模擬失真、誤差累積、計(jì)算速度慢及無(wú)法滿足突發(fā)污染事件實(shí)際決策要求等各種技術(shù)難題，提出了一種新型三維水動(dòng)力水質(zhì)模型，并利用實(shí)地流場(chǎng)觀測(cè)和水質(zhì)監(jiān)測(cè)數(shù)據(jù)進(jìn)行了模型校驗(yàn)。結(jié)果表明：該模型主要改進(jìn)了組合σ/z網(wǎng)格、總磷本構(gòu)方程、突發(fā)污染投藥處置模擬、OpenMP并行計(jì)算、水動(dòng)力預(yù)計(jì)算和雙時(shí)間步長(zhǎng)模式、污染團(tuán)擴(kuò)散軌跡示蹤等關(guān)鍵技術(shù)，顯著提升了模型的計(jì)算效率和模擬精度；通過(guò)在2023年漢江秋汛期間對(duì)丹江口水庫(kù)水質(zhì)安全管理的實(shí)際應(yīng)用，證明了該模型能有效支持水質(zhì)安全“四預(yù)”功能（預(yù)報(bào)、預(yù)警、預(yù)演、預(yù)案），模擬結(jié)果為漢江秋汛防御及汛后蓄水水質(zhì)安全保障提供了重要技術(shù)支撐。

關(guān) 鍵 詞：數(shù)字孿生；水質(zhì)安全；水質(zhì)模型；水動(dòng)力水質(zhì)模型；丹江口水庫(kù)

中圖法分類號(hào)：X832；TP391.985

文獻(xiàn)標(biāo)志碼：BDOI：10.16232/j.cnki.1001-4179.2024.08.030

0 引 言

丹江口水庫(kù)是南水北調(diào)中線工程水源地、國(guó)家一級(jí)水源保護(hù)區(qū)，保障其水質(zhì)安全是支撐京津冀豫社會(huì)經(jīng)濟(jì)高質(zhì)量發(fā)展、助力華北地區(qū)河湖生態(tài)環(huán)境復(fù)蘇、提升受水區(qū)人民獲得感、幸福感和安全感的重要抓手。為此，丹江口水利樞紐被列為水利部《“十四五”智慧水利建設(shè)規(guī)劃》確定的12個(gè)數(shù)字孿生先行先試重點(diǎn)建設(shè)水利工程之一，是國(guó)內(nèi)首個(gè)以水利工程水質(zhì)安全為業(yè)務(wù)需求的數(shù)字孿生水利工程^［1-2］。

作為數(shù)字孿生丹江口工程的核心模型之一，水質(zhì)數(shù)學(xué)模型是明晰庫(kù)區(qū)水環(huán)境時(shí)空演變特征、預(yù)測(cè)上游突發(fā)污染物遷移轉(zhuǎn)化過(guò)程的重要技術(shù)手段^［3-4］。近幾十年來(lái)，隨著計(jì)算機(jī)技術(shù)和應(yīng)用數(shù)學(xué)的快速發(fā)展，水質(zhì)模型已經(jīng)取得了顯著進(jìn)展。國(guó)外在水質(zhì)模型的開(kāi)發(fā)和應(yīng)用方面起步較早，諸如丹麥MIKE^［4］、美國(guó)WASP^［5］、EFDC^［6］和CE-QUAL-W2 ^［7］、荷蘭Delft3D^［8］，已被廣泛應(yīng)用于河流、湖泊、水庫(kù)等不同水體的水質(zhì)評(píng)估和管理。這些模型綜合考慮了水動(dòng)力學(xué)、化學(xué)反應(yīng)動(dòng)力學(xué)以及生物過(guò)程，能夠模擬多種污染物的遷移轉(zhuǎn)化過(guò)程。然而，這些模型在應(yīng)用時(shí)往往需要經(jīng)過(guò)大量的參數(shù)調(diào)優(yōu)和驗(yàn)證，且在處理復(fù)雜水體時(shí)的模擬時(shí)間較長(zhǎng)，難以實(shí)現(xiàn)實(shí)時(shí)預(yù)測(cè)。

國(guó)內(nèi)關(guān)于水質(zhì)模擬的研究起步較晚，但發(fā)展迅速，已經(jīng)形成了一批自主研發(fā)的水質(zhì)模型，如CJK3D^［9］等。這些模型在借鑒國(guó)外模型的基礎(chǔ)上，更加注重適應(yīng)中國(guó)特有的水環(huán)境特點(diǎn)和污染問(wèn)題，如對(duì)農(nóng)業(yè)面源污染、工業(yè)污染的特別考慮，在規(guī)劃環(huán)評(píng)、環(huán)境容量計(jì)算等方面得到了較多應(yīng)用^［10-13］。然而，現(xiàn)有模型在實(shí)現(xiàn)精確模擬的同時(shí)，仍面臨著模型復(fù)雜度高、計(jì)算效率低等問(wèn)題，特別是在處理突發(fā)水質(zhì)污染事件時(shí)，缺乏足夠的靈活性和實(shí)時(shí)性。

盡管國(guó)內(nèi)外現(xiàn)有的水質(zhì)模型已經(jīng)能夠支持一定程度上的水質(zhì)管理和決策，但現(xiàn)有技術(shù)和方法在應(yīng)對(duì)復(fù)雜和多變的水環(huán)境挑戰(zhàn)時(shí)，仍然面臨一系列的問(wèn)題：① 現(xiàn)有水質(zhì)模型在處理復(fù)雜的水體系統(tǒng)時(shí)往往難以達(dá)到高模擬精度。特別是對(duì)于具有復(fù)雜地形、變化水文條件和多種污染源的大型水體，模型往往難以準(zhǔn)確反映水質(zhì)變化趨勢(shì)，導(dǎo)致模擬結(jié)果與實(shí)際情況存在偏差。② 高精度的水質(zhì)模擬往往伴隨著高計(jì)算成本。特別是對(duì)于需要進(jìn)行長(zhǎng)期模擬或?qū)崟r(shí)預(yù)測(cè)的應(yīng)用場(chǎng)景，現(xiàn)有模型往往難以滿足快速響應(yīng)的需求。③ 實(shí)用化水平低和自主可控性差。盡管模型理論研究及開(kāi)發(fā)取得了一定進(jìn)展，但將模型轉(zhuǎn)化為實(shí)際應(yīng)用的工具仍然存在諸多障礙，且大多數(shù)先進(jìn)的水質(zhì)模型都是由國(guó)外研究機(jī)構(gòu)或公司開(kāi)發(fā)，在缺乏源碼自主可控的情況下，針對(duì)特定水體的模型優(yōu)化和定制開(kāi)發(fā)變得更加困難。

針對(duì)上述問(wèn)題，本文提出了基于數(shù)字孿生技術(shù)的丹江口水質(zhì)模型DJKWQM（danjiangkou water quality model），該模型通過(guò)改進(jìn)三維水質(zhì)模型垂向網(wǎng)格的分層方法、分區(qū)設(shè)置空間異質(zhì)的生化參數(shù)、并行化大規(guī)模線性方程組的求解過(guò)程、動(dòng)態(tài)解耦水動(dòng)力模型與水質(zhì)模型的數(shù)據(jù)流，顯著提升了模型計(jì)算效率和模擬精度，為全面構(gòu)建數(shù)字孿生丹江口水質(zhì)安全“四預(yù)”功能體系奠定了基礎(chǔ)，發(fā)揮了試點(diǎn)項(xiàng)目的引領(lǐng)示范作用。

1 現(xiàn)有水質(zhì)模型不足以支撐數(shù)字孿生“四預(yù)”體系

目前，盡管國(guó)內(nèi)外眾多水質(zhì)模型已被廣泛采用，但它們?cè)诶碚摶A(chǔ)或算法設(shè)計(jì)上存在不足，限制了其在數(shù)字孿生系統(tǒng)中的應(yīng)用。這些局限主要表現(xiàn)在以下5個(gè)關(guān)鍵方面：

（1）垂向網(wǎng)格坐標(biāo)系會(huì)產(chǎn)生梯度誤差。丹江口水庫(kù)最大水深167 m，是典型的深水型水庫(kù)，其流場(chǎng)和污染物濃度場(chǎng)具有極明顯的三維特征。三維水動(dòng)力水質(zhì)數(shù)學(xué)模型必須采用適當(dāng)?shù)拇瓜蜃鴺?biāo)系才能精確描述垂向變化特征。早期的z垂向坐標(biāo)系雖便于離散方程，但在淺水區(qū)域僅有一層垂向網(wǎng)格，垂向分辨率過(guò)低，近岸區(qū)域模擬效果差。后續(xù)發(fā)展出的垂向σ坐標(biāo)系可使深水區(qū)和淺水區(qū)具有相同的垂向分層數(shù)，保證了全域網(wǎng)格的垂向分辨率一致，但在陡峭深水區(qū)會(huì)產(chǎn)生明顯壓力梯度誤差^［14］，在應(yīng)用于庫(kù)底起伏大的丹江口庫(kù)區(qū)時(shí)，需要改進(jìn)上述垂向σ坐標(biāo)系以避免梯度誤差。

（2）模型過(guò)簡(jiǎn)或過(guò)繁。一部分模型對(duì)空間維度和關(guān)鍵過(guò)程進(jìn)行過(guò)多簡(jiǎn)化，無(wú)法精確刻畫污染物在湖庫(kù)內(nèi)的遷移轉(zhuǎn)化過(guò)程，實(shí)用性較弱，且該類模型多是將率定參數(shù)用于整個(gè)水域，未考慮環(huán)境變化和參數(shù)的空間異質(zhì)性，精確度較低。另一部分模型直接采用國(guó)外商業(yè)軟件的復(fù)雜水質(zhì)模塊進(jìn)行建模，如EFDC的water quality模塊和MIKE的ecolab模塊^［11］，此類模型涉及數(shù)百個(gè)參數(shù)，盡管這些模型在國(guó)外有諸多成功案例，但越復(fù)雜的水質(zhì)模型對(duì)輸入水質(zhì)數(shù)據(jù)的要求越精細(xì)，收集數(shù)據(jù)費(fèi)用也越高，模擬效果卻并未因考慮了更精細(xì)的生物地球化學(xué)循環(huán)過(guò)程而變得更好，性價(jià)比較低。因此，在模擬丹江口水庫(kù)水質(zhì)的時(shí)空變化過(guò)程時(shí)，需要針對(duì)本流域水污染特點(diǎn)和具體污染物，開(kāi)發(fā)出概化得當(dāng)、簡(jiǎn)繁適宜的水質(zhì)模塊。

（3）模型計(jì)算速度不能滿足應(yīng)對(duì)突發(fā)污染事件的決策需求。作為中線水源地，丹江口水庫(kù)及其上游如遇突發(fā)污染等緊急事件，需要在短時(shí)間內(nèi)預(yù)測(cè)分析污染擴(kuò)散范圍，以支撐水質(zhì)實(shí)時(shí)滾動(dòng)推演、快速預(yù)演決策，這對(duì)模型計(jì)算效率和精度都提出了較高的要求，即“不損失精度的條件下提高計(jì)算效率”。然而，現(xiàn)有三維水動(dòng)力水質(zhì)模型需要對(duì)大規(guī)模線性方程組進(jìn)行迭代計(jì)算以求解水位和流速分量，且多數(shù)模型是基于串行計(jì)算模式，在模擬水域面積達(dá)1 050 km²的丹江口水庫(kù)時(shí)，計(jì)算耗時(shí)明顯變長(zhǎng)，遠(yuǎn)不能滿足應(yīng)急決策的時(shí)限要求。雖可通過(guò)LU矩陣分解得到預(yù)條件子（preconditioner）以加快穩(wěn)定雙共軛梯度法（Bi-CGStab，Bi-conjugate gradient stabilized method）^［15］、廣義最小殘差法（GMRES，generalized minimum residual method）^［16］、共軛梯度法（CG，conjugate gradient method）^［17］等求解器的迭代速度，但提升效果較為有限，模型計(jì)算速度已成為阻礙數(shù)字孿生系統(tǒng)走向?qū)嶋H應(yīng)用的瓶頸。

（4）模型實(shí)用性不足。一方面是通用化程度低，模型前后處理沒(méi)有實(shí)現(xiàn)標(biāo)準(zhǔn)化，不同用戶根據(jù)自身需求定制不同模型，模型之間無(wú)法通用，移植到其他系統(tǒng)時(shí)仍需重新編制模型^［18］；另一方面，國(guó)外商業(yè)軟件僅支持離線計(jì)算，交互性差，無(wú)法滿足水質(zhì)滾動(dòng)推演、突發(fā)污染在線快速模擬等實(shí)時(shí)決策的業(yè)務(wù)場(chǎng)景需求。

（5）模型自主可控性差。目前稍復(fù)雜的水質(zhì)模型軟件大部分不開(kāi)源，少數(shù)開(kāi)源軟件僅提供包含基礎(chǔ)功能的早期代碼，對(duì)于大規(guī)模加速計(jì)算、動(dòng)態(tài)網(wǎng)格等關(guān)鍵底層技術(shù)仍實(shí)行封鎖。例如，EFDC模型僅公開(kāi)了基礎(chǔ)版源碼，必須支付昂貴的軟件費(fèi)才能使用包含MPI/OMP 混合多線程等核心技術(shù)的增強(qiáng)版，加上商業(yè)軟件的保護(hù)機(jī)制，其模型不能嵌入數(shù)字孿生系統(tǒng)“為我所用”，因此亟需針對(duì)核心“卡脖子”技術(shù)加強(qiáng)自主研發(fā)。

2 數(shù)字孿生丹江口水質(zhì)模型關(guān)鍵技術(shù)

針對(duì)現(xiàn)有水質(zhì)模型在應(yīng)用于丹江口水庫(kù)水質(zhì)模擬時(shí)的諸多問(wèn)題，利用開(kāi)源計(jì)算流體力學(xué)（CFD）類庫(kù)進(jìn)行二次開(kāi)發(fā)，構(gòu)建了丹江口水質(zhì)模型DJKWQM。DJKWQM交叉融合了水利科學(xué)與工程、環(huán)境科學(xué)、計(jì)算機(jī)可視化與并行計(jì)算技術(shù)等多個(gè)學(xué)科，在計(jì)算精度、效率、實(shí)用性等方面創(chuàng)新性地突破了多項(xiàng)關(guān)鍵技術(shù)，有力支撐了數(shù)字孿生丹江口水質(zhì)“四預(yù)”功能實(shí)現(xiàn)。

在開(kāi)發(fā)過(guò)程中，采用了基于有限差分的模型數(shù)值解法，同時(shí)對(duì)組合σ/z垂向網(wǎng)格技術(shù)和總磷本構(gòu)方程進(jìn)行了創(chuàng)新改進(jìn)，從而提升了水動(dòng)力及磷等關(guān)鍵水質(zhì)指標(biāo)的模擬精度。通過(guò)開(kāi)發(fā)OpenMP并行計(jì)算、水動(dòng)力預(yù)計(jì)算和雙時(shí)間步長(zhǎng)模式，提升了模型的計(jì)算效率。通過(guò)開(kāi)發(fā)突發(fā)污染投藥處置模擬技術(shù)、污染團(tuán)擴(kuò)散軌跡示蹤技術(shù)等，實(shí)現(xiàn)任意位置突發(fā)污染擴(kuò)散與處置快速模擬，以及污染物運(yùn)動(dòng)軌跡的實(shí)時(shí)追蹤監(jiān)控，滿足突發(fā)污染事件應(yīng)對(duì)的實(shí)際需求。最后通過(guò)開(kāi)發(fā)模型標(biāo)準(zhǔn)化接口，滿足數(shù)字孿生的模型實(shí)時(shí)演算需求。

模型總體框架和各功能模塊見(jiàn)圖1。主要包括8個(gè)功能模塊。其中，邊界條件讀取模塊負(fù)責(zé)從外部數(shù)據(jù)源（如氣象數(shù)據(jù)、水文數(shù)據(jù)等）中讀取邊界條件，包括流域的地形、流速、水位等信息。前處理模塊用于初始化模擬所需的各項(xiàng)變量，包括水位、流速、濃度等，并根據(jù)外部氣象數(shù)據(jù)計(jì)算水面上的風(fēng)應(yīng)力。紊流計(jì)算模塊采用適當(dāng)?shù)奈闪髂Ｐ?，?jì)算流場(chǎng)中的湍流運(yùn)動(dòng)。求解器模塊包括用于求解線性方程組的迭代法（如PCG法）和求解三對(duì)角矩陣的TMDA法。水動(dòng)力求解模塊（動(dòng)量方程）通過(guò)差分處理后求解動(dòng)量方程，得到水深平均的流場(chǎng)分布。水動(dòng)力求解模塊（三維流場(chǎng)、水位更新）在水深平均流場(chǎng)的基礎(chǔ)上，求解三維流場(chǎng)，并根據(jù)求解結(jié)果更新水位信息。物質(zhì)輸運(yùn)求解模塊考慮流體中物質(zhì)的對(duì)流和擴(kuò)散運(yùn)動(dòng)，根據(jù)流場(chǎng)信息計(jì)算物質(zhì)在流體中的傳輸過(guò)程（對(duì)流擴(kuò)散項(xiàng)），同時(shí)考慮水體中生物、化學(xué)反應(yīng)對(duì)物質(zhì)濃度的影響，以模擬水體中物質(zhì)的生物、化學(xué)變化（生化反應(yīng)項(xiàng)）。迭代控制管理模塊用于監(jiān)控模擬計(jì)算的進(jìn)行情況，判斷模擬是否達(dá)到設(shè)定的條件，如收斂性判定、時(shí)間步長(zhǎng)控制等，控制模型迭代計(jì)算流程。

2.1 模型數(shù)值計(jì)算流程

DJKWQM以完整形式的三維水動(dòng)力學(xué)方程（RANS方程）和污染物輸移方程為控制方程，利用有限差分法將上述偏微分方程組轉(zhuǎn)化為代數(shù)方程組。在每個(gè)時(shí)間步，首先將動(dòng)量方程與水位方程聯(lián)立，采用預(yù)處理共軛梯度法（PCG）求解關(guān)于所有單元格水位的大規(guī)模五對(duì)角線性方程組，得到各單元格的水位和水深平均流速，并將其傳遞至關(guān)于各層單元格垂向流速的三對(duì)角線性方程組，采用追趕法（TMDA）求解此方程組后得到各層單元格的垂向和水平流速分量，完成三維水動(dòng)力學(xué)模型的求解。以所得流場(chǎng)為基礎(chǔ)，考慮污染物的理化生過(guò)程，基于算子分裂技術(shù)^［19］分兩步求解污染物輸移方程，模擬污染物的遷移-擴(kuò)散-轉(zhuǎn)化全過(guò)程，第一步采用Smolarkiewicz和Clark高階迎風(fēng)差分格式^［20］求解對(duì)流項(xiàng)和反應(yīng)項(xiàng)（源匯項(xiàng)），第二步采用隱式迎風(fēng)差分格式求解擴(kuò)散項(xiàng)，完成水質(zhì)模型的求解。具體過(guò)程見(jiàn)圖2。

2.2 改進(jìn)的組合σ/z網(wǎng)格技術(shù)

部分國(guó)外先進(jìn)模型（如EFDC、MIKE3）采用了混合σ/z垂向網(wǎng)格方法^［21］，即從自由表面到指定深度使用傳統(tǒng)σ坐標(biāo)，指定深度以下則使用z坐標(biāo)，并在每個(gè)時(shí)間步內(nèi)根據(jù)新水深重新確定每個(gè)水平單元的垂向網(wǎng)格分層數(shù)，更新各變量所在的垂向坐標(biāo)。此法根據(jù)水深變化實(shí)時(shí)更新垂向分層總數(shù)，減小了梯度誤差，但也明顯增加了計(jì)算成本。

丹江口水庫(kù)蓄水運(yùn)行后水位相對(duì)穩(wěn)定，絕大部分區(qū)域的水深年內(nèi)變幅較小，水平單元的垂向網(wǎng)格分層數(shù)比較穩(wěn)定，使用混合σ/z垂向網(wǎng)格法在每個(gè)時(shí)間步都計(jì)算垂向網(wǎng)格分層數(shù)的必要性較弱。因此，對(duì)混合σ/z網(wǎng)格做了進(jìn)一步改進(jìn)（圖3），即在設(shè)置初始水深時(shí)確定計(jì)算域內(nèi)各水平單元的垂向分層數(shù)，且在模擬過(guò)程中保持不變，而考慮垂向各層厚度隨水深動(dòng)態(tài)變化而變化。這與現(xiàn)有混合σ/z垂向網(wǎng)格法需每步確定垂向分層數(shù)相比，大大減少了計(jì)算量，提高了計(jì)算效率。

2.3 重構(gòu)總磷模型本構(gòu)方程

磷是丹江口水質(zhì)管理中的重點(diǎn)關(guān)注指標(biāo)。考慮磷在水體中的降解沉降、底泥釋放、吸附解吸等關(guān)鍵動(dòng)力學(xué)過(guò)程（圖4），重構(gòu)總磷模型本構(gòu)方程，將磷的吸附沉降表達(dá)式、底泥磷釋放強(qiáng)度公式添加到方程反應(yīng)項(xiàng)（源匯項(xiàng)）。重構(gòu)后的磷模型包括對(duì)流項(xiàng)、擴(kuò)散項(xiàng)、反應(yīng)項(xiàng)（降解沉降項(xiàng)、底泥釋放項(xiàng)、泥沙與磷吸附解吸項(xiàng)），其中，降解沉降項(xiàng)采用一級(jí)反應(yīng)動(dòng)力學(xué)描述（即降解沉降速率給定一個(gè)經(jīng)過(guò)率定的固定常數(shù)），底泥釋放項(xiàng)采用修正的Elovich方程描述^［22］，泥沙與磷吸附解析項(xiàng)采用Langmuir動(dòng)力學(xué)方程描述^［23］，所有涉及的動(dòng)力學(xué)參數(shù)可通過(guò)原位觀測(cè)或室內(nèi)物理試驗(yàn)結(jié)果校正后獲得，改寫后的模型能更真實(shí)地描述磷在水體中的各類理化生過(guò)程，實(shí)現(xiàn)磷模擬精度的提高。

2.4 突發(fā)污染投藥處置模擬技術(shù)

在曾經(jīng)發(fā)生的丹江口水庫(kù)上游老灌河銻升高事件中，采取了筑橡膠壩與投放絮凝劑（硫酸亞鐵）的組合方式進(jìn)行應(yīng)急處置。銻污染物在水體中絮凝沉淀的快慢可采用絮凝速率來(lái)表征，傳統(tǒng)水質(zhì)模型通常將整個(gè)計(jì)算域的動(dòng)力學(xué)參數(shù)設(shè)置為一個(gè)固定常數(shù)，但在實(shí)際處置過(guò)程中，銻污染物的絮凝去除過(guò)程非常復(fù)雜，會(huì)隨著投藥量、投放時(shí)間、投放地點(diǎn)的不同而發(fā)生變化。

為滿足銻污染投藥處置的實(shí)際應(yīng)用需求，提出動(dòng)力學(xué)參數(shù)的動(dòng)態(tài)分區(qū)設(shè)置方法，首先對(duì)模型代碼中銻絮凝速率的對(duì)應(yīng)變量進(jìn)行擴(kuò)維；在此基礎(chǔ)上，結(jié)合有關(guān)文獻(xiàn)或機(jī)理試驗(yàn)研究得到投藥量、投藥時(shí)間與銻絮凝速率回歸關(guān)系的數(shù)學(xué)表達(dá)式，通過(guò)數(shù)字孿生系統(tǒng)獲取用戶輸入的投藥時(shí)間、地點(diǎn)、投藥量，根據(jù)投藥地點(diǎn)匹配對(duì)應(yīng)的計(jì)算網(wǎng)格單元后，計(jì)算出隨時(shí)間變化的絮凝速率，實(shí)現(xiàn)銻污染投藥處置過(guò)程的精準(zhǔn)模擬。假設(shè)丹江發(fā)生銻濃度升高，在武當(dāng)山水電站投藥點(diǎn)投放一定的水質(zhì)凈化絮凝劑進(jìn)行處置，圖5展示了投放絮凝劑后的丹江干流銻濃度分布。

由于丹江口水庫(kù)下墊面類型復(fù)雜、區(qū)域差異顯著，此次還將提出的動(dòng)力學(xué)參數(shù)動(dòng)態(tài)分區(qū)設(shè)置方法用于污染降解系數(shù)、底泥釋放速率等動(dòng)力學(xué)參數(shù)的分區(qū)設(shè)置，在模型中可以調(diào)整動(dòng)力學(xué)選項(xiàng)以及是否使用動(dòng)力學(xué)分區(qū)，更靈活地模擬復(fù)雜水質(zhì)遷移轉(zhuǎn)化過(guò)程。

2.5 OpenMP并行計(jì)算技術(shù)

為了支撐水質(zhì)實(shí)時(shí)滾動(dòng)推演、突發(fā)污染事件快速預(yù)演決策，開(kāi)發(fā)基于熱點(diǎn)分析的OpenMP并行算法，首先使用性能分析工具（如gprof或valgrind）分析程序中的瓶頸，識(shí)別占用大部分計(jì)算時(shí)間的相應(yīng)函數(shù)（熱點(diǎn)），確定待優(yōu)化函數(shù)，使用OpenMP進(jìn)行共享內(nèi)存的并行化處理，使多個(gè)CPU核心并行執(zhí)行計(jì)算任務(wù)。為避免并行化過(guò)程中的競(jìng)態(tài)條件和性能瓶頸，對(duì)部分模型數(shù)據(jù)結(jié)構(gòu)進(jìn)行重構(gòu)，使其更好地支持共享內(nèi)存操作。對(duì)比測(cè)試并行程序與原始串行程序，結(jié)果表明兩種程序計(jì)算結(jié)果幾乎一致（相對(duì)誤差小于3%），但并行程序的計(jì)算效率提升了30%。

2.6 水動(dòng)力預(yù)計(jì)算和雙時(shí)間步長(zhǎng)模式技術(shù)

傳統(tǒng)水動(dòng)力水質(zhì)模型的計(jì)算模式在每個(gè)時(shí)間步內(nèi)先后完成水動(dòng)力和水質(zhì)計(jì)算，再進(jìn)行下一時(shí)間步的計(jì)算。由于三維水動(dòng)力模型計(jì)算量大，并不能滿足突發(fā)污染過(guò)程預(yù)演及處置對(duì)模擬速度的要求。為此，開(kāi)發(fā)了水動(dòng)力預(yù)計(jì)算和雙時(shí)間步長(zhǎng)技術(shù)。該技術(shù)重塑了三維水動(dòng)力水質(zhì)模型的計(jì)算流程，重定向了模型數(shù)據(jù)流，即提前計(jì)算不同情景的庫(kù)區(qū)流場(chǎng)（預(yù)計(jì)算技術(shù)），在模擬突發(fā)污染時(shí)直接讀取流場(chǎng)信息驅(qū)動(dòng)水質(zhì)模型計(jì)算，不再執(zhí)行極為耗時(shí)的水動(dòng)力模型計(jì)算，從而實(shí)現(xiàn)了分鐘級(jí)地快速模擬污染物動(dòng)態(tài)變化。

該技術(shù)實(shí)現(xiàn)了水動(dòng)力預(yù)計(jì)算和雙時(shí)間步長(zhǎng)模式兩大技術(shù)突破。具體包括：首先將水動(dòng)力計(jì)算和水質(zhì)計(jì)算進(jìn)行動(dòng)態(tài)解耦，同時(shí)對(duì)水動(dòng)力模型的輸出數(shù)據(jù)結(jié)構(gòu)進(jìn)行重構(gòu)，避免生成的文件過(guò)大降低文件讀取速度；然后，由數(shù)字孿生系統(tǒng)為水質(zhì)計(jì)算提供各種不同場(chǎng)景的流場(chǎng)信息，例如在線推演功能提供當(dāng)日最新流場(chǎng)演算結(jié)果，歷史場(chǎng)景知識(shí)庫(kù)可提供歷史典型水文條件下庫(kù)區(qū)流場(chǎng)演算結(jié)果；最后，待模型解耦后，允許水動(dòng)力和水質(zhì)計(jì)算分別采用不同時(shí)間步長(zhǎng)（即雙時(shí)間步長(zhǎng)），在保證結(jié)果收斂的前提下（滿足CFL穩(wěn)定判據(jù)）^［24］，水質(zhì)模型可設(shè)置合適的時(shí)間步長(zhǎng)，大大加快模型整體計(jì)算效率。經(jīng)測(cè)試，某模擬時(shí)長(zhǎng)為5 d的突發(fā)污染情景（8萬(wàn)個(gè)計(jì)算單元）的計(jì)算時(shí)間為1～2 min，可有效支撐實(shí)際突發(fā)污染的快速模擬與決策。

2.7 污染團(tuán)演進(jìn)擴(kuò)散軌跡示蹤技術(shù)

由于突發(fā)性水污染事件具有隨機(jī)性及不可預(yù)見(jiàn)性的特點(diǎn)，根據(jù)實(shí)際管理需求，模型要能實(shí)現(xiàn)任意位置突發(fā)污染快速模擬并能隨時(shí)追蹤污染團(tuán)的位置，動(dòng)態(tài)識(shí)別污染團(tuán)與陶岔等關(guān)鍵點(diǎn)位的距離、到達(dá)時(shí)間等。

為此，開(kāi)發(fā)污染團(tuán)演進(jìn)擴(kuò)散軌跡示蹤技術(shù)，在每個(gè)計(jì)算時(shí)間步內(nèi)，通過(guò)水質(zhì)模型計(jì)算的全域濃度值進(jìn)行聚類分析判斷，識(shí)別污染團(tuán)前鋒位置所在網(wǎng)格，計(jì)算出其與目標(biāo)點(diǎn)位（如陶岔）所在網(wǎng)格的距離，結(jié)合流速情況估算其到達(dá)目標(biāo)點(diǎn)位的時(shí)間。該技術(shù)避免采用復(fù)雜拉格朗日示蹤機(jī)理模型進(jìn)行求解，為污染團(tuán)演進(jìn)動(dòng)態(tài)追蹤提供一種簡(jiǎn)單實(shí)用的方法，同時(shí)結(jié)合GIS地圖定位技術(shù)與模型進(jìn)行交互，實(shí)現(xiàn)任意突發(fā)污染源位置情景下的遷移擴(kuò)散計(jì)算。圖6展示了2023年4月，泗河排污口偷排高濃度氨氮事件被曝光后，通過(guò)模型對(duì)偷排事件污染擴(kuò)散過(guò)程進(jìn)行模擬分析的成果。

2.8 模型標(biāo)準(zhǔn)化和通用化

為使機(jī)理模型能夠滿足實(shí)時(shí)計(jì)算、動(dòng)態(tài)更新的需求，開(kāi)發(fā)了基于Service的獨(dú)立低耦合、高可配模型通用化封裝技術(shù)及標(biāo)準(zhǔn)化接口，實(shí)現(xiàn)系統(tǒng)發(fā)送計(jì)算指令傳遞參數(shù)-調(diào)用模型計(jì)算-計(jì)算結(jié)果三維可視化仿真的數(shù)據(jù)流接合和全流程貫通（圖7）。在分析水質(zhì)模型與數(shù)字孿生系統(tǒng)數(shù)據(jù)交互關(guān)系的基礎(chǔ)上，梳理明確數(shù)據(jù)流向，并指定模型輸入輸出的標(biāo)準(zhǔn)數(shù)據(jù)格式，數(shù)字孿生系統(tǒng)接到計(jì)算指令后，采集用戶設(shè)置的相關(guān)計(jì)算參數(shù)（如模型參數(shù)、水文水質(zhì)邊界條件），生成模型運(yùn)行所需的標(biāo)準(zhǔn)輸入文件并通過(guò)接口程序傳遞給模型，模型完成計(jì)算后，再通過(guò)接口程序?qū)⒔Y(jié)果反饋給數(shù)字孿生系統(tǒng)進(jìn)行展示和仿真渲染。同時(shí)按照統(tǒng)一標(biāo)準(zhǔn)對(duì)構(gòu)建好的各類專業(yè)模型（如洪水預(yù)報(bào)模型、一維水質(zhì)模型、三維水質(zhì)模型）進(jìn)行模塊化封裝和整合集成，并留出相關(guān)數(shù)據(jù)接口，便于規(guī)范化調(diào)用不同異構(gòu)計(jì)算模型。

3 模型校驗(yàn)

3.1 水動(dòng)力模型驗(yàn)證

為驗(yàn)證水動(dòng)力模型的有效性，課題組于2023年6月4日和9月27日，針對(duì)丹江口水庫(kù)肖川—龍口斷面和柳陂鎮(zhèn)斷面開(kāi)展流場(chǎng)現(xiàn)場(chǎng)觀測(cè)。丹江口肖川龍口斷面和柳陂鎮(zhèn)斷面的流速模擬值與監(jiān)測(cè)值見(jiàn)圖8。驗(yàn)證結(jié)果顯示，肖川—龍口斷面和柳陂鎮(zhèn)斷面的平均誤差分別為14.96%和10.18%，流速模擬誤差基本在15%以內(nèi)，表明丹江口水庫(kù)水動(dòng)力模型能比較準(zhǔn)確地反映丹江口水庫(kù)的流速變化情況。

3.2 水質(zhì)模型驗(yàn)證

課題組收集了2021～2022年丹江口水庫(kù)的逐月水質(zhì)監(jiān)測(cè)數(shù)據(jù)，對(duì)水質(zhì)模型進(jìn)行驗(yàn)證測(cè)試，總磷等關(guān)鍵指標(biāo)模擬與實(shí)測(cè)的過(guò)程曲線比較吻合，相對(duì)誤差在20%以內(nèi)，能較準(zhǔn)確反映庫(kù)區(qū)的水質(zhì)變化過(guò)程和分布，部分驗(yàn)證結(jié)果如圖9所示。

4 實(shí)例應(yīng)用

2023年9月下旬，受強(qiáng)降雨和漢江上游水庫(kù)調(diào)度影響，丹江口水庫(kù)水位持續(xù)上漲，漢江2023年第1號(hào)洪水在漢江上游形成，漢江汛情緊迫，洪水匯入將大量上游污染物帶入水庫(kù)，數(shù)字孿生系統(tǒng)發(fā)出漢江總磷升高告警。系統(tǒng)從9月29日起至事件結(jié)束，每日開(kāi)展?jié)L動(dòng)推演，采用當(dāng)日最新的水文、水質(zhì)監(jiān)測(cè)數(shù)據(jù)以及丹江口水庫(kù)預(yù)報(bào)調(diào)洪演算成果，調(diào)用數(shù)字孿生丹江口的三維水動(dòng)力水質(zhì)模型，實(shí)時(shí)研判庫(kù)區(qū)水質(zhì)安全形勢(shì)。

以10月7日的推演為例，當(dāng)日水文預(yù)報(bào)成果顯示，庫(kù)區(qū)預(yù)計(jì)在6 d后（10月13日）蓄水到170.00 m，基于9月29日至10月7日實(shí)測(cè)的庫(kù)區(qū)水文、水質(zhì)數(shù)據(jù)和10月7日的水文預(yù)報(bào)數(shù)據(jù)，推演9月29日至10月13日期間的庫(kù)區(qū)總磷分布和變化趨勢(shì)（即9月29日至10月7日為歷史復(fù)演，10月7～13日為預(yù)測(cè)預(yù)報(bào)），分析水庫(kù)蓄水到170.00 m過(guò)程中總磷在庫(kù)區(qū)內(nèi)的擴(kuò)散演進(jìn)過(guò)程及其對(duì)丹庫(kù)和陶岔水質(zhì)的影響。模型工況條件設(shè)置如下：漢江白河站的入庫(kù)總磷為0.17 mg/L；丹江磨峪灣的入庫(kù)總磷為0.106 mg/L；采用預(yù)報(bào)水文數(shù)據(jù)，6 d內(nèi)漢江入庫(kù)平均流量為1 901.56 m³/s，丹江的入庫(kù)平均流量為415 m³/s，大壩下泄平均流量為1 705 m³/s，陶岔平均引水流量為245 m³/s，水庫(kù)水位由169.37 m增加至170.00 m。

推演結(jié)果表明：自9月27日以來(lái)，高濃度總磷從漢庫(kù)和丹庫(kù)向庫(kù)中推進(jìn)，漢庫(kù)推進(jìn)很快，約在10月1日達(dá)到肖川—龍口；10月4日大約擴(kuò)散到壩前，整個(gè)漢庫(kù)總磷濃度比本底值有一定提高，大部分污染物隨大壩下泄，未進(jìn)入丹庫(kù)；4號(hào)以后水庫(kù)下泄流量逐漸減少，庫(kù)水位逐漸抬升，部分污染物逐漸向丹庫(kù)移動(dòng)；10月9日大約擴(kuò)散到?jīng)鏊印_(tái)子山，該處的總磷濃度模擬值約為0.038 mg/L（10月9月開(kāi)展了應(yīng)急監(jiān)測(cè)，總磷濃度監(jiān)測(cè)值為0.041 mg/L，模擬誤差在8%以內(nèi)，推演情況與實(shí)測(cè)數(shù)據(jù)趨勢(shì)基本一致）。根據(jù)對(duì)10月13日總磷的預(yù)測(cè)（圖10和圖11），水庫(kù)已蓄水到170.00 m，此時(shí)漢庫(kù)總磷濃度總體平穩(wěn)，丹江入庫(kù)總磷從河口大約擴(kuò)散至老城鎮(zhèn)，整個(gè)蓄水過(guò)程中，丹庫(kù)中心和陶岔水質(zhì)穩(wěn)定在地表水Ⅱ類標(biāo)準(zhǔn)。

每日根據(jù)最新實(shí)測(cè)數(shù)據(jù)，對(duì)前期推演結(jié)果進(jìn)行復(fù)盤和模型驗(yàn)證，并動(dòng)態(tài)優(yōu)化和調(diào)整模型參數(shù)。以10月7日為例，陶岔、馬蹬、青山、1號(hào)船陶岔界站、2號(hào)船倉(cāng)房香花鎮(zhèn)、3號(hào)船壩前站的模擬值與實(shí)測(cè)值見(jiàn)表1，從9月27日以來(lái)的推演結(jié)果來(lái)看，模型模擬值與實(shí)測(cè)值的相對(duì)誤差總體在12%以內(nèi)，基本可滿足實(shí)際模擬推演應(yīng)用的需求。

數(shù)字孿生丹江口工程為中國(guó)首個(gè)在大型水庫(kù)滿蓄中深度運(yùn)用的數(shù)字孿生工程，基于數(shù)字孿生丹江口水質(zhì)模型，實(shí)現(xiàn)了洪水期水質(zhì)滾動(dòng)推演與跟蹤分析，為保障漢江秋汛防御與丹江口水庫(kù)170 m蓄水水質(zhì)安全提供科學(xué)支撐。

5 結(jié)論與展望

本文詳細(xì)分析了現(xiàn)階段數(shù)字孿生丹江口水質(zhì)模型關(guān)鍵技術(shù)的技術(shù)難點(diǎn)，重點(diǎn)闡述了以三維水動(dòng)力水質(zhì)模型為核心的模型技術(shù)研發(fā)成果，最后介紹了數(shù)字孿生水質(zhì)模型在2023年漢江秋汛防御與丹江口水庫(kù)170 m蓄水過(guò)程中的實(shí)際演算應(yīng)用。

目前，數(shù)字孿生丹江口水質(zhì)模型技術(shù)在建設(shè)中已經(jīng)取得了一定實(shí)效，形成了一定的實(shí)踐積累，然而，離真正的“虛實(shí)孿生融合”的目標(biāo)仍有一定的距離，今后將以下5個(gè)方面作為技術(shù)突破方向：

（1）加強(qiáng)核心專業(yè)模型研發(fā)，將三維水動(dòng)力水質(zhì)模型1.0版本升級(jí)到2.0版本，所有模塊實(shí)現(xiàn)全面自主研發(fā)，同時(shí)引入AI智能模型，通過(guò)耦合機(jī)理模型與智能模型，實(shí)現(xiàn)參數(shù)自動(dòng)率定，提升水質(zhì)預(yù)報(bào)精度、延長(zhǎng)預(yù)見(jiàn)期。

（2）構(gòu)建自主可控、具有高保真度的流域水文面源-水動(dòng)力-水質(zhì)-生態(tài)動(dòng)力學(xué)多維耦合模型，進(jìn)一步完善水文-水質(zhì)-水生態(tài)全鏈條模擬推演功能。

（3）集成漢江、堵河、丹江、老灌河4條主要入庫(kù)河流水文水質(zhì)高頻同步監(jiān)測(cè)數(shù)據(jù)，同時(shí)推進(jìn)其余12條入庫(kù)支流的水文站和水質(zhì)自動(dòng)站建設(shè)，確保覆蓋所有關(guān)鍵入庫(kù)支流的實(shí)時(shí)數(shù)據(jù)獲取，從而使模型能夠更準(zhǔn)確地反映實(shí)際水質(zhì)狀況。

（4）研發(fā)非靜壓熱-流耦合模型、污染溯源模型等特殊水流和污染運(yùn)動(dòng)模擬模型，提高模擬精度、增加模型功能和應(yīng)用場(chǎng)景。同時(shí)加強(qiáng)并行加速算法自主開(kāi)發(fā)力度，推進(jìn)國(guó)家超算中心的水利共享應(yīng)用，增強(qiáng)高性能計(jì)算算力，提高模型演算速率。

（5）進(jìn)一步梳理與細(xì)化模型標(biāo)準(zhǔn)化需求，建立開(kāi)放兼容、穩(wěn)定成熟的模型平臺(tái)，并形成相應(yīng)的核心技術(shù)標(biāo)準(zhǔn)。

參考文獻(xiàn)：

［1］蔡陽(yáng).以數(shù)字孿生流域建設(shè)為核心構(gòu)建具有\(zhòng)"四預(yù)\"功能智慧水利體系［J］.中國(guó)水利，2022（20）：2-6.

［2］林莉，李全宏，曹慧群，等.數(shù)字孿生丹江口水質(zhì)安全建設(shè)挑戰(zhàn)與舉措［J］.中國(guó)水利，2023（11）：32-36.

［3］胡春宏，郭慶超，張磊，等.數(shù)字孿生流域模型研發(fā)若干問(wèn)題思考［J］.中國(guó)水利，2022（20）：7-10.

［4］舒長(zhǎng)莉，李林，馮韜.基于MIKE21的河道飲用水源地突發(fā)污染事故模擬：以贛江南昌段為例［J］.人民長(zhǎng)江，2019，50（3）：73-77.

［5］李一平，唐春燕，龔然.水環(huán)境數(shù)學(xué)模型原理及應(yīng)用［M］.北京：科學(xué)出版社，2021.

［6］HAMRICK J M.A three-dimensional environmental fluid dynamic computer code：theoretical and computational aspects［R］.Virginia：Virginia Institute of Marine Science，1992.

［7］李光浩，楊霞，陳和春，等.三峽水庫(kù)香溪河庫(kù)灣水華水動(dòng)力調(diào)控初步研究［J］.人民長(zhǎng)江，2017，48（10）：18-23，29.

［8］張繼民，張新周，湯紅亮，等.Delft3D在海灣電廠溫排水?dāng)?shù)值模擬中的應(yīng)用［J］.人民長(zhǎng)江，2009，40（1）：59-62.

［9］白一冰，張成剛，羅小峰.呂泗漁場(chǎng)人工魚礁群流場(chǎng)效應(yīng)及穩(wěn)定性研究［J］.人民長(zhǎng)江，2018，49（8）：25-30.

［10］賈鵬，王慶改，周俊，等.地表水環(huán)評(píng)數(shù)值模擬精細(xì)化研究［J］.環(huán)境影響評(píng)價(jià)，2015，37（1）：51-54.

［11］李一平.太湖生態(tài)動(dòng)力學(xué)模型研究［M］.北京：中國(guó)水利水電出版社，2014.

［12］JING Z，KANG L，YAO H M，et al.Simulation of water temperature by using urban lake temperature model［J］.Journal of Hydraulic Engineering，2018，144（1）：06017025.

［13］KANG L，GUO X M.Hydrodynamic effects of reconnecting lake group with Yangtze River in China［J］.Water Science and Engineering，2011，4（4）：405-420.

［14］黑鵬飛，周剛，假冬冬，等.σ-數(shù)值切割單元法（σ-SIBM）：解決σ-坐標(biāo)在海洋陡坡區(qū)域數(shù)值模擬中壓強(qiáng)梯度誤差的新方法［J］.中國(guó)科學(xué)（地球科學(xué)），2014，44（4）：753-765.

［15］郭曉明，田治宗，李書霞，等.基于Keller-Box格式的三維非靜壓水動(dòng)力學(xué)模型［J］.華中科技大學(xué)學(xué)報(bào)（自然科學(xué)版），2015（7）：29-33.

［16］SAAD Y F.Iterative methods for sparse linear systems［M］.Philadelphia：Society for Industrial and Applied Mathematics，2003.

［17］GUO X M，LING K.A new efficient depth-integrated non-hydrostatic model for free surface flows［J］.Science China Technological Science，2013，56（4）：824-830.

［18］黃艷.以數(shù)字孿生長(zhǎng)江支撐流域治理管理［J］.中國(guó)水利，2022（8）：30-35.

［19］賈宏恩，李開(kāi)泰，鐘賀.算子分裂法求解對(duì)流-擴(kuò)散-反應(yīng)方程［J］.工程數(shù)學(xué)學(xué)報(bào)，2012，29（1）：89-95.

［20］申霞，洪大林，王鵬.基于POM的物質(zhì)輸運(yùn)方程數(shù)值格式研究［J］.海洋科學(xué)進(jìn)展，2009（4）：452-459.

［21］孫先，楊文俊，王國(guó)棟.基于σ坐標(biāo)變換和水位函數(shù)法的立面二維水動(dòng)力學(xué)模型［J］.長(zhǎng)江科學(xué)院院報(bào)，2012，29（7）：6-9，26.

［22］盧鑫.硫化物對(duì)電鍍廠鉻污染土壤穩(wěn)定化處理及其長(zhǎng)期穩(wěn)定性研究［D］.南昌：華東理工大學(xué)，2017.

［23］王穎，郭世花.湖泊底泥中磷釋放的隨機(jī)動(dòng)力學(xué)模型［J］.水利學(xué)報(bào)，2003，34（11）：71-77.

［24］柳尚斌，孫懷鳳，李文翰，等.瞬變電磁BEDS-FDTD三維正演新算法及穩(wěn)定性驗(yàn)證［J］.地球物理學(xué)報(bào)，2023，66（2）：841-853.

（編輯：劉 媛）

Development and application of key technologies of water quality model in

digital twin danjiangkou Project

JING Zheng¹，CAO Huiqun¹，LIN Li¹，QIN He²，GUO Xiaoming³

（1.Water Environment Department，Changjiang Scientific Research Institute，Wuhan 430010，China; 2.South-to-North Water Diversion Middle Line Water Source Co.，Ltd.，Shiyan 442700，China; 3.Yellow River Civilization and Sustainable Development Research Center，Henan University，Zhengzhou 450046，China）

Abstract：Digital twin technology，as an innovative advancement in recent years，achieves high-precision mapping of the real physical world by constructing accurate virtual models.This provides a new perspective and approach to addressing the issues in traditional water quality safety management.Water quality models are central to realizing precise mapping and virtual-real interaction of water quality in digital twin lakes and reservoirs.To address the water quality safety management needs of Danjiangkou Reservoir，the source of the Middle Route of the South-to-North Water Diversion Project，a novel three-dimensional hydrodynamic water quality model was proposed.This model overcame various technical challenges of traditional models，including simulation distortion，error accumulation，slow computing speed，and the inability to meet decision-making requirements for sudden pollution events.The model was calibrated using in-situ flow field observations and water quality monitoring data.Results demonstrated that the model significantly enhances computational efficiency and simulation accuracy through key technologies such as combined σ/z grid，total phosphorus constitutive equation，simulation of sudden pollution treatment，OpenMP parallel computing，hydrodynamic pre-computation，dual time-step mode，and pollution plume diffusion trajectory tracing.Practical application during the Hanjiang River autumn flood in 2023 proves the model′s effectiveness in supporting the \"four-pre\" functions of water quality safety （forecast，pre-warning，rehearsal，and pre-emergency planning）.The simulation results provide crucial technical support for flood prevention and post-flood water quality safety assurance in the Hanjiang River.

Key words：digital twin; water quality safety; water quality model; hydrodynamic water quality model; Danjiangkou Reservoir