基于改進重疊投影法的空間耦合水污染評估模型

張 影，楊曉東，賈子璇，章 胤，王 晶

(1.燕山大學里仁學院，河北 秦皇島 066004； 2.燕山大學理學院，河北 秦皇島 066004；3.燕山大學經濟管理學院，河北 秦皇島 066004； 4.燕山大學區域經濟發展研究中心，河北 秦皇島 066004)

水污染評估是研究水污染問題的重要途徑，利用水動力水質數值模型模擬計算水污染的輸移擴散過程，可為水污染風險評估和相關部門采取應對措施提供重要參考信息。水庫作為水資源重要存儲地，庫區上游突發性水污染的研究一直備受關注。針對庫區上游河床特點，往往建立一維、二維耦合水動力水質模型模擬污染物的擴散，于庫區中上游河道窄長區域構建一維河網模型[1]，下游河道開闊水域構建二維模型，并在交接面處進行合理的耦合，既可較準確地反映污染物濃度的時空變化規律，同時又可控制模型的復雜度以滿足實際需求[2]。

近年來，眾多學者針對一維、二維耦合模型開展了大量研究工作。諸裕良等[3]在河口一維、二維連接處通過接口斷面法傳遞水力因子，建立了一種一維、二維全隱河網海灣水動力聯網數學模型；Lai等[4]提出了適用于大尺度水動力模擬的一維、二維耦合方法，并對長江中游流域的河流-湖泊洪水演進進行了模擬；陳文龍等[5]構造并求解Riemann問題實現了一維、二維模型耦合，有效克服了傳統堰流公式缺點，并提出時間步長自適應匹配方法解決了一維、二維模型時間步長不一致的問題；王秀杰等[6]建立了復雜條件下天然河道漫潰堤洪水在防洪保護區的一維、二維水動力模型；顧杰等[2]基于MIKE FLOOD建立了入海河流及近岸海域一維、二維耦合河流-海岸水動力和水質模型；田福昌等[7]建立了山洪溝道潰堤洪水演進一維、二維水動力耦合數值模擬模型以分析評估潰堤山洪淹沒風險。目前，一維、二維耦合模型研究多為基于堰流公式和數值通量的側向型聯解耦合，構建水動力的一維、二維耦合模型以研究漫潰堤洪水問題，而針對一維、二維水動力水質耦合模型的研究較少，多為MIKE軟件的實際應用。此外，針對耦合交接面的連接問題，對參數及連接條件的過度概化，往往造成模擬結果的不準確。

針對上述問題，本文在傳統重疊投影法[8]的基礎上，于重疊區域的上邊界引入流速等物理量的二次橫向分布函數，基于神經網絡訓練思想，采用最速下降法構建訓練算法，對分布函數進行優化訓練，將收斂后的分布函數作為模型的連接條件，進而建立了基于改進重疊投影法的空間耦合模型，有效地降低了耦合處物理量的突變，實現了一維、二維水動力和水質模型精準耦合。

1 一維、二維水動力與水質模型

1.1 一維河網水動力模型

采用描述明渠非恒定流的水動力運動方程：

(1)

式中：x為沿程距離；t為時間；z為水位；Q為流量；qt為旁側入流流量；B為寬度；A為斷面面積；K為流量模數；g為重力加速度。

結合河網汊點連接條件，利用Preissmann四點加權隱式格式[9]離散控制方程，然后采用河網三級聯解法求解(詳見文獻[10])。

1.2 一維河網水質模型

河道污染物的一維對流擴散方程為

(2)

式中：c為水流輸送的污染物濃度；Ex為污染物縱向擴散系數；K1為污染物綜合降解系數。

對方程(2)采用前差分離散時間項、隱式迎風格式離散對流項、中心差分離散擴散項，整理得各斷面濃度遞推關系式，并結合汊點質量平衡方程建立方程組求解，回代至各河段得各斷面污染物濃度(詳見文獻[11])。

1.3 平面二維水動力模型

二維控制方程有水流連續性方程和水流運動方程組：

(3)

(4)

式中：H為水深；u、v為x、y方向的流速；f為阻力系數；Ω為科氏力系數；ν為紊流渦黏性系數；λ為風應力系數；wx、wy為風速在x、y方向上的分量。

采用貼體坐標法[12-13]，通過Poison方程對控制方程進行坐標變換，并采用交替方向隱格式法(ADI)求解(詳見文獻[14])。

1.4 平面二維水質模型

二維對流擴散方程用以描述污染物在水體中的輸移擴散規律：

(5)

式中Ex、Ey為x、y方向的擴散系數。

對流擴散方程經坐標變換后，基于水動力模型計算的流速、水位等結果，采用交替方向隱格式法求解(詳見文獻[14])。

2 一維、二維模型的改進耦合模型

采用一維、二維耦合的水動力和水質模型模擬污染物的擴散時，在耦合交接面的連接區域上對參數及連接條件的過度概化，往往造成模擬結果的不準確。一般情況下，沿河寬方向各物理量——如水位、流速及污染物濃度等——不是均勻分布的，為體現連接后二維區域各物理量的真實分布情況，在傳統重疊投影法的基礎上，于一維、二維區域引入水位、流速、污染物濃度等物理量的二次橫向分布函數，進而構建訓練算法，對各分布函數中的未知參數進行訓練，將收斂后的分布函數作為模型的連接條件，建立改進的一維、二維水動力和水質空間耦合模型。

2.1 耦合模型連接條件的改進

考慮在一般情況下沿河寬方向各物理量不是均勻分布，在一維、二維交接面設置物理量沿河寬方向的分布函數Ψ(y)，并采用二次函數形式來擬合其在邊界的分布：

Ψ(y)=az1y2+az2y+az3

(6)

式中：y為河寬方向坐標；az1、az2、az3為待定系數。將分布函數Ψ(y)于二維邊界離散化，設其在控制體i(i=1，2，…，N,N為二維邊界處控制體個數)的取值為Ψi，有：

(7)

由于二維邊界物理量隨時間變化的模擬計算過于復雜，因而本文代之以適用于一定時段的二維邊界各控制體的波動比值ωi，有:

(8)

此時，通過訓練得到波動比值向量ω，在已知交接面處的平均值后就可得該物理量的分布情況，能準確體現其在交接面的連接關系，從而提高耦合準確度。

據此，在二維邊界控制體上設置水位、流速、污染物濃度分布函數，需要注意的是重疊區域二維邊界的流速需要縱向流速分布函數u(y)以及流速矢量方向與縱向的夾角分布函數θ(y)兩個分布函數以表示其分布情況，分布函數皆采用二次函數形式擬合，通過訓練可得收斂的水位波動比值向量ωz、流速波動比值向量ωu、夾角分布向量θ，進而得到污染物濃度波動比值向量ωc，作為污染物濃度在交接面上的連接條件。

2.2 改進的耦合模型的實現步驟

考慮到二維區域數值波的傳播受柯朗-弗里德里西斯-列維(Courant-Fredrich-Lewy， CFL)條件限制，進行水域延伸構造虛擬重疊區域[8]，并考慮在一維準確計算水域中進行二維區域物理量的細化分布函數的訓練，將二維計算水域的邊界向一維計算水域延伸Courant個網格點，見圖1。此時，定義左側虛擬邊界為Ω1，右側一維、二維交接面為Ω2，其中，Ω1為一維、二維轉化平面，Ω2為比較訓練平面；圖1中Φ為一維區域內的物理量，Ψ為二維區域的物理量，L表示重疊區域中虛擬的一維斷面數。

圖1 重疊-投影示意圖

訓練時，面Ω2采用滯后條件，即下一時步虛擬重疊區域的物理量通過該時間步長內的迭代計算得到含參精確值。在一維準確計算水域中進行二維區域物理量的細化分布函數的訓練，利用該精確值與原滯后值之差構建目標函數以訓練分布函數相關系數。

以訓練穩定后的分布函數作為模型一維、二維的連接條件進行耦合計算，進而構建改進的重疊投影區域。改進的重疊投影關系見圖2。

圖2 重疊區域一、二維連接關系示意圖

2.2.1參數訓練步驟

a. 先計算水動力模型，面Ω2的時間滯后條件取為ΔΦn+1=0(n為模型已完成迭代的次數)，記此時的邊界值為ΦΩ2。

b. 求解一維隱格式，得面Ω1的物理量值ΦΩ1，并乘以波動比值向量ωz、ωu將面Ω1上一維物理量值轉化為二維邊界條件ΨΩ1(P1)。

c. 計算二維區域。取面Ω2上流速等物理量投影得到一維精確邊界條件Φ′Ω2。

d. 確定目標函數p(P1):

p(P1)=(Φ′Ω2-ΦΩ2)2

(9)

目標函數越小，一維、二維連接效果越好。

2.2.2模型耦合步驟

a. 面Ω2的時間滯后條件取為ΔΦn+1=0，記此時的邊界值為ΦΩ2。

b. 將滯后邊界值ΦΩ2作為一維水動力模型的下邊界，采用三級聯解法得到各水力要素在一維計算區域各斷面的取值。

d. 求解二維計算區域，代入二維邊界條件ΨΩ1(P1)，采用交替方向隱格式法結合追趕法得水力要素于二維計算區域的分布。

e. 取面Ω2的水力要素分布，投影得到下一時刻一維邊界條件Φ′Ω2，即ΦΩ2,n+2。

f. 將水動力模型相關水力參數輸入水質模型。同水動力模型訓練耦合步驟，并由面Ω1的污染物濃度波動比值向量ωc實現污染物濃度一維、二維區域的連接，最終得到整個計算區域內的污染物濃度分布。

g. 重復上述步驟，根據實時訓練出的面Ω1的分布函數，實現一維、二維交接面物理量的轉化，直至模擬結束。

3 模型檢驗

以秦皇島桃林口庫區上游青龍河流域為研究區域，模擬青龍縣某污水泄漏事件，對上述改進的一維、二維耦合河網水動力水質模型進行檢驗和分析。

3.1 研究區域概況

桃林口水庫位于河北省東北部的灤河主要支流青龍河上，總庫容8.59億m3，區域內地勢北高南低，降雨主要集中于每年的7、8月。庫區上游河流以青龍河水系為主，河長265 km，控制面積3 431.5 km2，河床主要為砂卵石，中游河床寬50～70 m，下游段河床較寬，為400～1 000 m，平均縱坡0.155 6%。

通過分析區域的河段組成，確定一維河網研究范圍為：上邊界取自小柳條溝北莊S358縣道，下邊界取至小菜峪。徑流量較大支流包括沙河、起河、都源河以及星干河。一維模型共布設了505個斷面，典型斷面如圖3所示，斷面平均距離Δx=500 m。為滿足CFL條件和控制方程的微分性質，設定一維模型的時間步長Δt=120 s。

圖3 河網各河段典型斷面位置概化示意圖

二維模型研究范圍包括小菜峪斷面至入庫口。為適應復雜河道邊界，本文建立貼體坐標系將物理計算區域進行轉化[15]，最終生成了3 451個計算網格。基于CFL條件選擇時間步長為4 s。一維、二維計算結果的輸出間隔均為120 s。

3.2 定解條件

3.2.1初始條件

以青龍河水系監測站最近(2019年7月15日)的一次檢測數據為基礎，通過數據處理，得到各典型斷面的流量、水位、自然狀態下污染物的濃度等數據。采用一維、二維插值方法，分別得到一維、二維區域各斷面、控制體的流量、水位等數據作為水動力模型的初始條件；各斷面、控制體的自然狀態下污染物的濃度作為水質模型的初始條件。

3.2.2邊界條件

一維水動力模型的上邊界條件是各干支流河段上游實際的流量、水位變化，下邊界條件為一維、二維耦合界面處二維模型計算得到的平均流量和平均水位。

二維水動力模型的耦合界面為流量-水位邊界，其余邊界為固壁邊界，本文根據流域特點和模擬時長，固壁邊界采用閉邊界條件，沿邊界法線(nΓ0)方向流速為零。

河網水質模型的邊界條件是單源點污染源，位于青龍縣河南村斷面的滿源污水處理廠的生活污水均勻地排入所在的河流系統中。

3.3 參數計算方法

以秦皇島市水文局等部門提供的部分數據為基礎，通過Google Earth遙感與圖像處理技術[16-17]設計提取算法，提取河網各斷面和控制體的基礎數據，并結合河流縱向擴散系數[18]等水力學經驗計算公式，計算得到模型必要的基本參數。對于某些需要率定的參數，按斷面采用試錯法調試有較大任意性，河網規模較大時調試工作量十分巨大，本文采用分級法，即同一等級的河道按一定的參數取值，保證模型參數的相對準確性及可行性，并通過實地勘察和咨詢水力學專業人員等途徑，綜合考量，對一維、二維每個經典斷面和經典控制體的水力參數進行打分，采用插值方法計算所有斷面和控制體的水力參數。

3.4 模擬場景

假定滿源污水處理廠突發污水泄漏事故，有 9 400 m3的生活污水均勻地排入研究水系中，本文主要研究濃度為1.606×10-4mol/L的總磷(TP)在未來26.7 h內的變化情況[19]，并從秦皇島市水文局監測站監測時間(2019-07-15)開始計時。

圖4和圖5為河網水動力模型在一維區域的計算結果。在整個模擬期間，最大流量是876.23 m3/s，在一個時間步長(120 s)內變化最大為14.17 m3/s，最高水位是216.34 m，一個時間步長內變化最大為0.23 m。

圖4 斷面流量的時空變化

圖5 斷面水位的時空變化

圖6為河網水質模型在一維區域的計算結果，可見污染事件發生后561個時間步長(即18.7 h)污染物到達了二維區域(小菜峪斷面)。值得注意的是，污染物匯入汊點時未污染水流匯入稀釋，導致輸出斷面的污染物濃度出現突變。

圖6 污染物TP濃度的時空變化

圖7為河網模型在二維區域的計算結果，即在污染事件發生后581、601個時間步長(即19.33 h、20.23 h)時，二維區域中各控制體上污染物濃度的分布，可以看出，當以二次函數作為交接面污染物濃度分布的擬合函數時，后續的二維區域同一橫截面的控制體上的TP濃度，仍然呈中間高、兩端低近似二次函數趨勢的分布，從而驗證了用二次函數作為擬合函數的合理性。

(a) t=19.33 h

由以上計算結果可知，通過最速下降法訓練改進后的重疊投影法實現一維、二維模型的耦合，在河段推演、汊點連接、一維和二維耦合等過程均未出現解的“奇點”，方法的改進并未引起模型解的突變，表明改進后的模型穩定性較高，達到了預期目標。

首次訓練穩定后的一維和二維交接面處縱向流速、流速矢量方向與縱向的夾角、水位和污染物濃度的分布函數中的參數如表3所示，可以看到，4個分布函數的二次項系數均較小，表明河流同一截面上的該4個變量數值相差不會太大；夾角θ的二次系數為負，河流中央夾角小，兩端夾角大，其他3個變量正好相反；各分布函數的一次項系數很小，表明河流中該4個水力要素近似呈對稱分布。以上結果近似符合青龍河在二維區域的河流特征，進而驗證了參數訓練算法的正確性。

表3 首次訓練穩定后各物理量分布函數的系數

圖8為通過傳統的重疊投影法、改進的重疊投影法計算的一維、二維交接面與其下游相鄰20 m斷面處的平均縱向流速差值和TP平均濃度差值的變化，傳統重疊投影法計算的平均縱向流速差值的平均值為0.091 m/s，改進的重疊投影法計算的平均值為0.040 m/s，實測的平均值為0.046 m/s，可見改進的重疊投影法使一維向二維過渡時流速等水動力參量的突變性大大減小，保證了水動力耦合模型的精確性。由于流速和TP濃度的非均勻分布(河中央流速大處TP濃度高)，改進后的水質耦合模型污染物TP擴散速度稍微加快，更加符合實際情況。

(a) 平均縱向流速差值

由于傳統的重疊投影法在處理一維、二維數值交換時，是將交接面上各物理量的一維數值賦予該截面所有二維控制體，即交接面上所有控制體上各物理量的數值均相等，而河流的實際情況是各物理量的數值在橫截面上呈某種分布，因而導致了河網模型的解在一維、二維耦合處出現較大的突變。針對庫區上游河流，本文通過引入分布函數的方法可實現減小此種突變的目的。

4 結 語

本文于一維、二維交接面處引入流速、水位等物理量的二次分布函數，通過最速下降法對分布函數中的系數進行訓練，將訓練得到的分布函數作為物理量從一維向二維過渡的連接，提高了耦合的精確性。桃林口庫區上游青龍河流域實例驗證結果表明，模型耦合的精確性得到了很大的提高，且計算的穩定性較高。相比較傳統的重疊投影法，由于增加了最速下降法訓練參數的過程，計算的時間復雜度有了較大提升，但隨著并行計算算法和超級計算機的快速發展，改進的重疊投影法適用性會得到更大的提高。改進的模型適用于上中游河道狹窄、下游河面開闊，且在一維、二維交接面河道寬度變化較緩的河流發生的突發性污染場景。本文選用二次函數作為各物理量在交接面分布的擬合函數，尋找更精確的分布函數是下一步的研究方向。