基于HEC－RAS和反問題算法的突發(fā)污染溯源方法研究

文章編號：1674-6139（2025）07-0015-06

中圖分類號：X522文獻標(biāo)志碼：A

Identification of Sudden Pollution Sources in Rivers Using HEC-RAS and Inverse Problem Algorithm

Lin Yiyuan

（FujianProvicialKyboratoryfnldingogujadefingschoduzh

Abstract：Amethodologicalframeworkwasestablishedtorealizetheall-weatheranddynamicsupervisionofsuddenwaterpolution sources.Theproposedmodel integrates theconceptsofaBayesian-multi MarkovchainMonteCarlomethodandHEC-RAShydrody namicmodel.Themethodologicalframeworkwastestedusingtwohypotheticalcasesofapotentialsudenastewaterspillincidentand multi-pointsourceshichthelocationinformationisknown.Thestudyresultsshowthat：Forthetrackingofsuddenwastewaterdischarge，heinverseproblemmodelcanefectivelyestimatesourceparametersincludingsourcelocation，sourceflowrate，startingand endingtieofdischarge.Forthdyamicsupervisionofmultiplepolutionsouces，theinverseproblemmodelcanfctivelyidentifythe dischargeamoutsandtarting/endingtimeofmultiplepolutionsources.Thedevelopedinversemodelcouldprovidetechicalsolutiofor dynamic monitoring of sudden water pollution with the support of on-line data in the future.

Keywords：polltion source tracking；inverse problem；hydrodynamic model；HEC-RAS；Bayesianalgorithm

前言

突發(fā)水污染事件具有不確定性、隱蔽性和非連續(xù)性等特點，監(jiān)管難度大，確定污染排放信息是開展應(yīng)急處置和風(fēng)險防控的前提。目前在突發(fā)污染溯源領(lǐng)域的研究，主要可以分為示蹤劑實驗法和數(shù)值模擬法兩大類。清華大學(xué)[]基于三維熒光光譜與水樣一一對應(yīng)的原理，研制出污染預(yù)警溯源儀，用于異常情況污染溯源，但需同步建立工業(yè)、生活污染源的水質(zhì)指紋圖譜，目前多數(shù)地區(qū)暫未達到該條件。與示蹤劑實驗法不同，數(shù)學(xué)模擬法采用反問題的思想，依據(jù)環(huán)境水力學(xué)監(jiān)測數(shù)據(jù)來反推污染源信息，求解污染物源項識別問題[2]

但現(xiàn)有研究多是針對單點源污染排放的溯源，基于污染物輸移、稀釋和轉(zhuǎn)化機制，根據(jù)監(jiān)測斷面的污染物濃度數(shù)據(jù)來反演源項參數(shù)的解[3]?？紤]到河道側(cè)向入流是污染排放的間接表征，會引起河道流量/水位的變化，本論文從水動力角度出發(fā)，基于水動力模型和反問題算法展開單點非固定源/多點固定源的突發(fā)污染溯源的方法學(xué)研究。與水質(zhì)模型相比，水動力模型涉及的不確定參數(shù)少，且水位可方便全天候在線監(jiān)測。研究為水污染溯源工作提供新的參考思路。

1模型的基本原理

1. 1 一維非恒定流模型

由于河道縱向長度通過遠大于其寬度和深度，可將河道水流流動簡單概化為一維問題，水動力過程可用一維圣維南方程組表示[4]，即

式（1）中，A為橫斷面面積， m²;t 為時間，s; Q 為斷面過流流量， m³/s;x 為河道縱向距離， m;h 為斷面水深， m;g 為重力加速度， m/s²;S₀ 和 S_f 分別為河道底坡度和摩阻坡度，無量綱； q_l 為單位長度河道上的旁側(cè)入流量， m²/s 。

1.2水動力模型的二次開發(fā)

HEC－RAS軟件在水文水力計算中應(yīng)用廣泛[5]，最大的優(yōu)勢在于可以自由下載以及與其他模型/軟件耦合使用。研究在前人的基礎(chǔ)上[對HEC-RAS（版本：HEC-RAS4.1.0）進行二次開發(fā)，基于Matlab語言調(diào)用HEC－RAS的非穩(wěn)態(tài)流文件（*.u##），修改計算節(jié)點入流模塊的相關(guān)參數(shù)來開展水污染溯源模擬研究。

在Matlab語言環(huán)境下調(diào)用HEC-RAS進行水力計算，主要分為打開和運行程序、更改輸入文件、讀取、比較結(jié)果和輸出4個功能模塊。（1）打開和運行程序：涉及函數(shù)主要包括創(chuàng)建HEC－RAS軟件對應(yīng)的COM服務(wù)器：actxserver函數(shù)、打開項目文件（h.Project_Open 函數(shù)）、運行計劃文件（h.Compute_CurrentPlan函數(shù)）；（2）更改輸入文件：借助strrep函數(shù)，將非穩(wěn)態(tài)文件（*.u##）中目標(biāo)數(shù)值進行替換;借助fget1函數(shù)逐行讀取文件內(nèi)容；（3）讀取結(jié)果：借助h.OutputDSS_GetStageFlow函數(shù)讀取指定河道斷面的流量/水位數(shù)據(jù)；（4）比較和輸出：通過之前已經(jīng)整理的實測水位數(shù)據(jù)（obs.mat），以似然函數(shù)為判斷依據(jù)，輸出相對理想的參數(shù)后驗概率密度分布。

2 反問題算法

2.1基于貝葉斯－蒙特卡洛算法的不確定性分析

研究采用在不確定性量化方面具有優(yōu)勢的改進貝葉斯-蒙特卡洛算法（Bayesian－MCMC）進行污染反演分析，貝葉斯推理基于以下公式[7]：

式（2）中， θ 為模型參數(shù)； 為觀測數(shù)據(jù)； 為參數(shù)的后驗概率密度函數(shù)； p（θ） 為參數(shù)的先驗概率密度函數(shù)； 為似然函數(shù) 是 的規(guī)范化常數(shù)。

似然函數(shù) 通過水動力模型模擬值 Y= {y₁，y₂，…，y_t，…，y_n} 與觀測數(shù)值 ， 的之差的模數(shù)來表征模擬值與實測值之間的擬合程度。定義觀測數(shù)據(jù)噪聲為白噪聲，其誤差服從均值為0，標(biāo)準(zhǔn)偏差為 σ 的正態(tài)分布， θ ）可表示為：

式（3）中， .y_t（θ） 為第 Ψ_tΨ_t 個節(jié)點的模擬值； 為第 χ_t 個節(jié)點的觀測值； n 為觀測數(shù)據(jù)個數(shù)。

2.2突發(fā)污染源模型框架

研究采用動態(tài)多鏈的采樣方法，引入空間采樣及異常鏈檢驗相關(guān)算法，對經(jīng)典的隨機游走Metropolis 算法進行改進[8]。污染溯源模型系統(tǒng)運行框架見圖1，以下為運行流程：

（1）確定模型參數(shù)維度 d ，根據(jù)參數(shù)先驗分布信息隨機產(chǎn)生 N 個樣本；

（2）調(diào)用水動力模型計算出樣本初始狀態(tài)對應(yīng)的監(jiān)測斷面水位值，根據(jù)式（3）計算各樣本的后驗概率密度 ：

（3） N 條Markov鏈平行進化。在第 Ψ_t 次（ ，…，T} ）迭代中，根據(jù)交叉概率 CR 選擇參數(shù)子空間A進行交叉操作，判斷建議分布中各分量 θ_p^i，j（i= {1，…，d} ）是否進行交叉替換；

（4）計算第j條Markov鏈 Θ={θ_t-1¹，…，θ_t-1^N} 在第 χ_t 次迭代的跳轉(zhuǎn)步長 dθ_ιⁱ 及按建議分布產(chǎn)生的候

選狀態(tài) θ_p^i，j

（5）計算建議分布在候選樣本狀態(tài)下的后驗概率密度函數(shù) ，通過接受該判斷是否接受候選樣本 θ_pⁱ ，否則仍保留 θ_t-1ⁱ 5

（6）在非平穩(wěn)期采用 IQR （Inter-Quartile-Range，IQR）[9]進行異常鏈的統(tǒng)計與剔除；

（7）依據(jù)Gelmanamp;Rubin方法計算各參數(shù)當(dāng)前迭代次數(shù)的 統(tǒng)計值，判斷Markov鏈?zhǔn)欠衿椒€(wěn)。

（8）重復(fù) 直至滿足設(shè)定的迭代次數(shù) T 同時結(jié)合平穩(wěn)期樣本參數(shù)迭代值進行后驗概率密度統(tǒng)計分析。

3仿真數(shù)值實驗

3.1案例I：單個工業(yè)點源突發(fā)污染排放

3.1.1 問題描述

設(shè)計案例I如下：一條河流的某河段均勻，在距起始端 4km 斷面安裝了水位在線監(jiān)測設(shè)備，監(jiān)測頻次為每分鐘一次。在距離起始端 1.2km 處某工業(yè)企業(yè)產(chǎn)生了突發(fā)性排放，排放起始時間為 59min （從0：00起開始記錄，下同），終止時間為 180min 。排放流量為 0.2m³/s 。對應(yīng)的反問題為根據(jù)斷面水位在線監(jiān)測數(shù)據(jù) 反演某突發(fā)污染排放的位置 x 流量 排放起始時刻 T_?1 、結(jié)束時間刻 T₂ 。河道水文參數(shù)及待反演參數(shù)

信息見表1。

3.1.2 反演分析

在改進的Bayesian－MCMC算法中，取 N=5 ，并設(shè)定迭代次數(shù) T=10 000 ，利用拉丁超立方抽樣法從先驗分布中抽取生成排放位置、水量和排放起始時刻的初始樣本。觀測數(shù)據(jù) 定義為河道下游監(jiān)測斷面的水位連續(xù)監(jiān)測數(shù)據(jù)。參數(shù) 的搜索空間根據(jù)監(jiān)測位置而定，該工況下設(shè)置為［0，4000]；參數(shù)q 的搜索空間設(shè)置為[0，0.1];根據(jù)水位監(jiān)測數(shù)據(jù)可以判斷出工業(yè)企業(yè)排放時間一定在0：00一5：00a.m 區(qū)間內(nèi)（從 開始以分鐘記錄），故參數(shù) ^T1，T2 的搜索空間設(shè)置為[0，500]。

圖2展示了4個參數(shù)邊緣后驗概率分布，并給出了最大后驗密度概率估計（MaximumAPosterioriEstimation，MAP），其中MAP在圖中用藍色的 × 表示?？梢钥闯?，Markov鏈進入穩(wěn)定收斂區(qū)域后，4個參數(shù)的MAP均與真值完全吻合。

結(jié)果表明，針對 4km 長河道中單個點源突發(fā)污染排放的設(shè)計算例，反演參數(shù)包括污染源排放位置、排放水量、排放開始時刻、排放結(jié)束時刻。溯源模型反演值與真實值完全吻合，污染溯源精度可以達200m ，模型的可靠性得到驗證。

3.2案例Ⅱ：多個工業(yè)點源動態(tài)監(jiān)管

3.2.1 問題描述

設(shè)計案例Ⅱ如下：一條河流的某河段均勻，在距起始端 4km 斷面安裝了水位在線監(jiān)測設(shè)備，監(jiān)測頻次為每分鐘一次。在監(jiān)測設(shè)備上游存在3家工業(yè)企業(yè)，在一天內(nèi)3家企業(yè)均存在污水排放現(xiàn)象。企業(yè)A在距離起始端 1.2km 處，污水排放起始時間為59min （從 起開始記錄，下同），終止時間為

180min ，排放水量為 0.2m³/s ;企業(yè) B 在距離起始端1.8km 處，污水排放起始時間為 209min ，終止時間為270min ，排放水量為 0.1m³/s ;企業(yè)C在距離起始段2.6km 處，污水排放起始時間為 1019min ，終止時間為1 140min ，排放水量為 0.05m³/s 。則對應(yīng)的反問題為根據(jù)斷面在線監(jiān)測數(shù)據(jù) （h_ι）^～，…，（h_n）^～} 反演上游3家工業(yè)企業(yè)排污口排污的起始時間 T_g01，T_g02，T_g03 ，結(jié)束時間 T_g11?T_g12?T_g13 及對應(yīng)的排放水量qg1qg2qg3o

3.2.2 反演分析

在改進的Bayesian－MCMC算法中，取 N=7 ，并設(shè)定迭代次數(shù) T=20 000 ，利用拉丁超立方抽樣法從先驗分布中抽取生成各工業(yè)污染點源的排放水量和對應(yīng)的排放起始時刻的初始樣本。觀測數(shù)據(jù) 定義為河道下游監(jiān)測斷面的水位連續(xù)監(jiān)測數(shù)據(jù)。參數(shù) q_gl?q_g2?q_g3 的搜索空間均設(shè)置為 [0，0.1];T_g01 、T_g02?T_g03?T_g11?T_g12?T_g13 的搜索空間均設(shè)置為[0，1440]（從 0：00a.m 開始以分鐘記錄）。

圖3展示了9個參數(shù)邊緣后驗概率分布，并給出了最大后驗密度概率估計（MAP）。與案例I工況不同，由于存在多解的情況，多條Markov鏈經(jīng)過2萬次迭代后仍不能進入穩(wěn)定收斂階段，如參數(shù) q_g1 的搜索軌跡會在 0.2，0.1，0.05 三個值間跳躍， T_g01 的搜索軌跡在60、1020間跳躍，即會出現(xiàn)三個排放點源信息錯位的情況。但是通過足夠多次的迭代，得到的各參數(shù)概率密度曲線仍能提供準(zhǔn)確的反演信息。由反演結(jié)果（見圖3和表2）可以看出，9個參數(shù)的MAP均與真值幾乎完全吻合，最大誤差也不到 1% 。反演結(jié)果說明即使在多點源排放，反演組合爆炸的情況下，改進的MCMC模型仍能反演得到準(zhǔn)確的污染源排放信息。

4結(jié)論

針對突發(fā)污染具有不確定性、隱蔽性，難以通過常規(guī)手段進行日常監(jiān)管的特點，研究構(gòu)建了基于反問題算法的污染溯源模型，應(yīng)用于突發(fā)污染排放動態(tài)反演中。文中給出了兩類溯源仿真設(shè)計算例，一類是針對單個點源大流量、突發(fā)性排放的溯源追蹤;另一類則是針對多個涉污工業(yè)企業(yè)排放的動態(tài)監(jiān)管問題。案例I、Ⅱ均表明，針對高維、不適定性的環(huán)境水力學(xué)反問題，構(gòu)建的溯源模型能夠給出精度較高的反演結(jié)果。特別是案例Ⅱ突破了以往多數(shù)研究僅針對單個點源污染進行溯源的技術(shù)瓶頸，后續(xù)可首先在有條件的工業(yè)園區(qū)開展該方法體系的實例驗證與推廣應(yīng)用，必要時可與水質(zhì)監(jiān)測配合，助力河道水污染智慧監(jiān)管。

參考文獻：

[1]熊秋燃，沈鑒，胡遠，等.基于水質(zhì)熒光指紋技術(shù)的復(fù)合污染河道溯源研究[J].光譜學(xué)與光譜分析，2024，44（6）：1773-1780.