文章編號(hào)：1674-6139（2025）07-0015-06

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

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

1模型的基本原理

1. 1 一維非恒定流模型

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

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

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

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

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

2 反問題算法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3仿真數(shù)值實(shí)驗(yàn)

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

3.1.1 問題描述

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

信息見表1。

3.1.2 反演分析

在改進(jìn)的Bayesian－MCMC算法中，取 N=5 ，并設(shè)定迭代次數(shù) T=10 000 ，利用拉丁超立方抽樣法從先驗(yàn)分布中抽取生成排放位置、水量和排放起始時(shí)刻的初始樣本。觀測數(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è)排放時(shí)間一定在0：00一5：00a.m 區(qū)間內(nèi)（從 開始以分鐘記錄），故參數(shù) ^T1，T2 的搜索空間設(shè)置為[0，500]。

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

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

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

3.2.1 問題描述

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

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

3.2.2 反演分析

在改進(jìn)的Bayesian－MCMC算法中，取 N=7 ，并設(shè)定迭代次數(shù) T=20 000 ，利用拉丁超立方抽樣法從先驗(yàn)分布中抽取生成各工業(yè)污染點(diǎn)源的排放水量和對(duì)應(yīng)的排放起始時(shí)刻的初始樣本。觀測數(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個(gè)參數(shù)邊緣后驗(yàn)概率分布，并給出了最大后驗(yàn)密度概率估計(jì)（MAP）。與案例I工況不同，由于存在多解的情況，多條Markov鏈經(jīng)過2萬次迭代后仍不能進(jìn)入穩(wěn)定收斂階段，如參數(shù) q_g1 的搜索軌跡會(huì)在 0.2，0.1，0.05 三個(gè)值間跳躍， T_g01 的搜索軌跡在60、1020間跳躍，即會(huì)出現(xiàn)三個(gè)排放點(diǎn)源信息錯(cuò)位的情況。但是通過足夠多次的迭代，得到的各參數(shù)概率密度曲線仍能提供準(zhǔn)確的反演信息。由反演結(jié)果（見圖3和表2）可以看出，9個(gè)參數(shù)的MAP均與真值幾乎完全吻合，最大誤差也不到 1% 。反演結(jié)果說明即使在多點(diǎn)源排放，反演組合爆炸的情況下，改進(jìn)的MCMC模型仍能反演得到準(zhǔn)確的污染源排放信息。

4結(jié)論

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

參考文獻(xiàn)：

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