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

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

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

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

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

1 一維、二維水動力與水質(zhì)模型

1.1 一維河網(wǎng)水動力模型

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

(1)

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

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

1.2 一維河網(wǎng)水質(zhì)模型

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

(2)

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

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

1.3 平面二維水動力模型

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

(3)

(4)

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

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

1.4 平面二維水質(zhì)模型

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

(5)

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

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

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

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

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

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

Ψ(y)=az1y2+az2y+az3

(6)

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

(7)

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

(8)

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

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

2.2 改進的耦合模型的實現(xiàn)步驟

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

圖1 重疊-投影示意圖

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

以訓練穩(wěn)定后的分布函數(shù)作為模型一維、二維的連接條件進行耦合計算，進而構(gòu)建改進的重疊投影區(qū)域。改進的重疊投影關系見圖2。

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

2.2.1參數(shù)訓練步驟

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

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

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

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

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

(9)

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

2.2.2模型耦合步驟

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

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

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

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

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

g. 重復上述步驟，根據(jù)實時訓練出的面Ω1的分布函數(shù)，實現(xiàn)一維、二維交接面物理量的轉(zhuǎn)化，直至模擬結(jié)束。

3 模型檢驗

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

3.1 研究區(qū)域概況

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

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

圖3 河網(wǎng)各河段典型斷面位置概化示意圖

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

3.2 定解條件

3.2.1初始條件

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

3.2.2邊界條件

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

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

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

3.3 參數(shù)計算方法

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

3.4 模擬場景

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

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

圖4 斷面流量的時空變化

圖5 斷面水位的時空變化

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

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

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

(a) t=19.33 h

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

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

表3 首次訓練穩(wěn)定后各物理量分布函數(shù)的系數(shù)

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

(a) 平均縱向流速差值

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

4 結(jié) 語

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