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

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

中圖分類號：X522文獻標志碼：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

前言

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

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

1模型的基本原理

1. 1 一維非恒定流模型

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

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

1.2水動力模型的二次開發

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

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

2 反問題算法

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

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

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

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

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

2.2突發污染源模型框架

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

（1）確定模型參數維度 d ，根據參數先驗分布信息隨機產生 N 個樣本；

（2）調用水動力模型計算出樣本初始狀態對應的監測斷面水位值，根據式（3）計算各樣本的后驗概率密度 ：

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

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

選狀態 θ_p^i，j

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

（6）在非平穩期采用 IQR （Inter-Quartile-Range，IQR）[9]進行異常鏈的統計與剔除；

（7）依據Gelmanamp;Rubin方法計算各參數當前迭代次數的 統計值，判斷Markov鏈是否平穩。

（8）重復 直至滿足設定的迭代次數 T 同時結合平穩期樣本參數迭代值進行后驗概率密度統計分析。

3仿真數值實驗

3.1案例I：單個工業點源突發污染排放

3.1.1 問題描述

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

信息見表1。

3.1.2 反演分析

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

圖2展示了4個參數邊緣后驗概率分布，并給出了最大后驗密度概率估計（MaximumAPosterioriEstimation，MAP），其中MAP在圖中用藍色的 × 表示。可以看出，Markov鏈進入穩定收斂區域后，4個參數的MAP均與真值完全吻合。

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

3.2案例Ⅱ：多個工業點源動態監管

3.2.1 問題描述

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

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

3.2.2 反演分析

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

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

4結論

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

參考文獻：

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