灞河流域水環(huán)境數(shù)值模擬研究

陳 軍，吳 程，李 爽，王 穎，龔 毅

(1 中國電建集團市政規(guī)劃設計研究院有限公司，廣東 珠海 519000； 2 中國電建集團昆明勘測設計研究院有限公司，云南 昆明 650032；3 西安理工大學 水利水電學院，陜西 西安 710048)

我國水資源時空分布不均，隨著全球氣候變化的加劇，水資源短缺與水環(huán)境惡化問題愈顯突出。城市水系周邊人口密集、工企業(yè)聚集、污染源眾多，人口快速增長與城市化進程加快造成了一系列城市水環(huán)境問題[1-2]。城市水系承擔著調(diào)蓄防洪、維持區(qū)域生物多樣性、氣候微調(diào)和景觀文娛等功能[3]，是城市經(jīng)濟發(fā)展和生態(tài)平衡維系的重要依托。如何對城市水系水環(huán)境進行高效管理是目前研究熱點之一[4-5]。

作為水環(huán)境問題快速識別與高效管理的重要技術手段之一，水環(huán)境模型得到了廣泛應用[6-10]。水環(huán)境模型按照計算域尺度可分為水域模型與流域模型，其中水域模型最具代表性的有DHI MIKE21/MIKE3、EFDC(environmental fluid dynamics code)和Delft3D等，流域模型最具代表性的有SWAT(soil & water assessment tool)和HSPF(hydrological simulation program-fortran)[11-13]等。目前相關研究以水質(zhì)指標空間分布與典型控制斷面(點位)的指標變化特征以及水環(huán)境對各類工程響應特性的研究為主[14-16]，但利用水環(huán)境模型對水質(zhì)數(shù)據(jù)進行空間擴展后再進行綜合水質(zhì)評價的研究較少。

灞河是西安市第二大河流，下游城區(qū)段受鄉(xiāng)村、城鎮(zhèn)等典型人類活動的影響，在浐灞橡膠群壩工程實施后，河道內(nèi)形成了較為穩(wěn)定的水域面積，水體流速減緩，使得河道流態(tài)更加接近湖泊，形成了典型的兼具河道與湖泊特征的湖泊型城市河道。前期已有學者利用水環(huán)境模型針對灞河做了如下研究，Li等[17]利用SWAT模型研究了灞河流域非點源污染空間分布特征且識別了重點源區(qū)；胡勝等[18]基于數(shù)字濾波法和SWAT模型，借助空間分析與全局趨勢分析方法對灞河流域年均基流空間分布特征進行了解析，為流域地下水資源開發(fā)與保護提供了決策依據(jù)。綜上，目前關于水環(huán)境模型在灞河層面的應用研究多以灞河流域水文和非點源污染特征為主，較少涉及到灞河和浐河水域尺度。

本研究以灞河流域下游城區(qū)段灞河和浐河水域為對象，利用DHI MIKE21FM模型中HD模塊耦合Ecolab模塊，選取溶解氧(dissolved oxygen,DO)、化學需氧量(chemical oxygen demand,COD)、葉綠素a(chlorophyll a,Chl-a)、氨氮(ammonia nitrogen,NH3-N)和總磷(total phosphorus,TP)5個水質(zhì)指標作為狀態(tài)變量，參考WQ模板中各狀態(tài)變量的生化反應方程對其修改后進行水環(huán)境模型構建，利用監(jiān)測數(shù)據(jù)對模型校準后，開展灞河和浐河水環(huán)境數(shù)值模擬研究，分析灞河口水質(zhì)指標變化特性?；跀?shù)值模擬結果利用BP(back propagation)神經(jīng)網(wǎng)絡算法對研究區(qū)時均綜合水質(zhì)類型進行評價，解析基準年、豐水年、平水年與枯水年研究區(qū)時均綜合水質(zhì)類型空間分布特征，利用主成分分析(principal component analysis,PCA)法識別了研究區(qū)水質(zhì)類型的主導性指標，以期為湖泊型城市河道水環(huán)境數(shù)值模擬研究提供參考。

1 研究區(qū)域概況

灞河流域位于陜西省中部，西安市東南部，發(fā)源于秦嶺北麓藍田縣灞源鄉(xiāng)九道溝，自東南流向西北，最終于灞橋區(qū)三郎村附近匯入渭河。灞河總長104.1 km，河床平均比降6.2‰，流域面積2 581 km2，流域內(nèi)現(xiàn)有羅李村和馬渡王2座水文站[19]。灞河主要支流有浐河、藍橋河與輞峪河。浐河系灞河下游最大支流，全長64 km，河床平均比降8.9‰，流域面積752.8 km2。灞河流域?qū)贉貛О霛駶櫞箨懶约撅L氣候，多年平均氣溫13.3 ℃，多年平均降水量720 mm，多年平均蒸發(fā)量946.6 mm，降水時空分布不均，降水量由南向北遞減，降水主要集中在6-10月。灞河流域內(nèi)有5區(qū)1縣，分別為商州區(qū)、長安區(qū)、雁塔區(qū)、灞橋區(qū)、浐灞生態(tài)區(qū)和藍田縣[20]。流域上游地類以林地為主，中游以耕地為主，下游以建設用地為主。灞河是西安市重要的飲用水源地之一[19]，也是西安市經(jīng)濟與生態(tài)發(fā)展的重要載體。研究區(qū)以灞河流域下游城區(qū)段內(nèi)灞河和浐河為主，南起隴海線北至灞河入渭河口，沿水流方向自東南向西北貫穿浐灞生態(tài)區(qū)核心區(qū)。灞河流域及研究區(qū)域位置如圖1所示。

圖1 灞河流域位置及各監(jiān)測點位的分布Fig.1 Distribution of Ba River basin and monitoring point

2 研究數(shù)據(jù)來源

(1)水文數(shù)據(jù)。1958-2012年水文數(shù)據(jù)由西安理工大學圖書館館藏《黃河流域水文資料》(紙質(zhì)版)整理得到。(2)氣象數(shù)據(jù)。氣象數(shù)據(jù)由中國氣象數(shù)據(jù)網(wǎng)(http://data.cma.cn)獲取。(3)水陸邊界數(shù)據(jù)。利用2016年Sentinel-2兩景影像(影像數(shù)據(jù)來源于美國地質(zhì)勘探局USGS網(wǎng)站https://earthexplorer.usgs.gov，影像ID：S2A_OPER_MSI_L1C_TL_SGS__20160429T034347_20160429T084621_A004448_T49SCT_N02_01_01、S2A_OPER_MSI_L1C_TL_SGS__20160429T034347_20160429T084621_A004448_T49SCU_N02_01_01)，采用開源GIS工具中最大似然法解譯后再根據(jù)Google衛(wèi)片調(diào)整后獲取。(4)水下地形數(shù)據(jù)。將現(xiàn)場實測水深、水位數(shù)據(jù)導入GIS軟件后，將水深與水位數(shù)據(jù)轉(zhuǎn)換成河底高程數(shù)據(jù)，再與河道岸線高程數(shù)據(jù)合并后，得到整個研究區(qū)水下地形數(shù)據(jù)。(5)水工建筑物數(shù)據(jù)。來源于文獻[21]。(6)水質(zhì)數(shù)據(jù)。于2012年7-10月開展了5期水質(zhì)監(jiān)測(2012-07-10，2012-07-25，2012-08-24，2012-09-25和2012-10-15)，其中DO與水溫2項指標由便攜式多參數(shù)水質(zhì)監(jiān)測儀進行檢測，COD、NH3-N、TP和Chl-a 4項指標均參照《水和廢水監(jiān)測分析方法(第4版)》[22]進行檢測。(7)點源污染數(shù)據(jù)。研究區(qū)共有6個污染點源，其流量(于研究點源處多次測量10 L容器注滿點源廢水所需時間，對時間取算術平均值后換算得到流量數(shù)據(jù))和水質(zhì)(DO、水溫、COD、NH3-N、TP和Chl-a)與水質(zhì)數(shù)據(jù)(6)同期檢測。(8)水深數(shù)據(jù)。由便攜式超深波水深儀測定，與水質(zhì)數(shù)據(jù)(6)同期測定。

3 水環(huán)境模型構建與數(shù)據(jù)處理

3.1 水環(huán)境模型構建

利用DHI MIKE21FM模型中的HD模塊耦合Ecolab模塊，選取DO、COD、Chl-a、NH3-N和TP共5個水質(zhì)指標作為狀態(tài)變量，參考WQ模板中各狀態(tài)變量的生化反應方程并對其修改后進行水環(huán)境模型構建，修改后各狀態(tài)變量的生化反應控制方程為[23]：

(1)

(2)

(3)

(4)

(5)

式中：C1表示DO的質(zhì)量濃度；A表示大氣富氧率；S表示光合作用富氧率；B表示呼吸作用耗氧率；CSC2表示COD降解過程中耗氧率；CSC4表示NH3-N降解過程中耗氧率；C2表示COD的質(zhì)量濃度；DC2表示COD降解率；C3表示Chl-a的質(zhì)量濃度；PC3表示Chl-a凈產(chǎn)率；DEC3表示Chl-a死亡率；C4表示NH3-N的質(zhì)量濃度；PC2-C4表示COD降解過程中NH3-N產(chǎn)率；DC4表示NH3-N降解率；C5表示TP的質(zhì)量濃度；PC2-C5表示COD降解過程中TP產(chǎn)率；DC5表示TP降解率。

3.2 水環(huán)境模型設置

研究區(qū)計算域2個起點分別為灞河6號壩下游100 m處和浐河4號壩下游50 m處，終點為灞河1號壩下游10 m處(圖1)。采用非結構化網(wǎng)格進行計算域網(wǎng)格劃分，經(jīng)網(wǎng)格無關性分析后，本次數(shù)值模擬最優(yōu)網(wǎng)格數(shù)量為9 254個。模擬時間設置為2012-07-10-10-15，時間步長取10 s。上游水動力邊界統(tǒng)一采用流量邊界，下游采用水位邊界。上游水質(zhì)邊界均采用各指標實測濃度輸入，下游水質(zhì)邊界采用第二類邊界條件(諾依曼邊界條件)。模型中點源均以簡單源進行概化，水動力輸入條件統(tǒng)一采用流量輸入，水質(zhì)輸入條件統(tǒng)一采用濃度輸入。由于監(jiān)測數(shù)據(jù)期數(shù)限制，故本研究未將實測數(shù)據(jù)分為率定與驗證兩個階段，而是將率定與驗證階段進行合并。

3.3 計算工況設置

本研究共設置基準年、豐水年、平水年與枯水年4個計算工況，基準年水動力邊界條件采用2012年《黃河流域水文資料》(紙質(zhì)版)中灞河和浐河的流量與水位數(shù)據(jù)，水質(zhì)邊界條件采用2012年DO、COD、Chl-a、NH3-N和TP質(zhì)量濃度的實測數(shù)據(jù)，豐水年、平水年與枯水年的水質(zhì)邊界條件與基準年保持一致，水動力邊界條件分別為豐水年、平水年與枯水年相應流量。豐水年、平水年和枯水年相應流量由以下步驟計算：對馬渡王水文站1958-2012年的年平均流量進行P-Ⅲ曲線擬合后，選取Cs=2Cv(其中Cs表示偏態(tài)系數(shù)，Cv表示變差系數(shù))曲線確定豐水年(保證率P=10%)、平水年(保證率P=50%)和枯水年(保證率P=90%)年平均流量，選取典型年后對年平均流量進行年內(nèi)分配，分別得到豐水年、平水年與枯水年的水動力邊界條件。

3.4 BP神經(jīng)網(wǎng)絡訓練

為對研究區(qū)水質(zhì)類型進行準確評估，引入BP神經(jīng)網(wǎng)絡法算法，依據(jù)地表水環(huán)境質(zhì)量標準(GB 3838-2002)選取DO、COD、NH3-N和TP的標準限值，利用python3語言numpy包中的random模塊下的uniform函數(shù)生成4×300的數(shù)組作為BP神經(jīng)網(wǎng)絡輸入層，并將1～6的十進制數(shù)值進行二進制轉(zhuǎn)化后生成1×300的數(shù)組，作為BP神經(jīng)網(wǎng)絡輸出層，進行BP神經(jīng)網(wǎng)絡訓練，其中隱藏層層數(shù)采用黃金分割法[24]最終確定為13，輸入層與隱藏層之間的激勵函數(shù)采用tansig函數(shù)(雙曲正切S型傳輸函數(shù))，隱藏層與輸出層之間的激勵函數(shù)采用purelin函數(shù)(線性傳輸函數(shù))，樣本訓練過程中DO中Ⅰ類水標準限值定為0 ℃下的水體飽和溶解氧質(zhì)量濃度(14.64 mg/L)。

3.5 數(shù)據(jù)處理

數(shù)據(jù)統(tǒng)計分析、圖形可視化表達以及BP神經(jīng)網(wǎng)絡訓練均采用python3語言。計算結果云圖采用DHI MIKE自帶的Plot Composer進行處理。

4 結果與分析

4.1 模型精度驗證

由于監(jiān)測數(shù)據(jù)完整性的限制，因此選取圖1中的S1、S4點位的水深及S1、S2、S3和S4點位的DO、COD、NH3-N、TP與Chl-a數(shù)據(jù)對水環(huán)境模型模擬值與實測值進行對比分析，為方便后續(xù)統(tǒng)計，本研究將研究時段2012-07-10-10-15轉(zhuǎn)化為以天計數(shù)，即2012-07-10為第1天，共98 d，水環(huán)境模型模擬值與實測值的對比結果見圖2。從圖2可以觀察到，各指標的實測值基本分布于模擬值兩側，且模擬值的變化趨勢與實測值基本保持一致。

圖2 灞河流域S1、S2、S3和S4監(jiān)測點位水環(huán)境模型模擬值與實測值的對比Fig.2 Comparison of simulated and measured values of water environment model at S1,S2,S3 and S4 monitoring points of Ba River basin

為了對灞河流域水環(huán)境模型精度進行量化分析，引入百分比偏差系數(shù)(PBIAS)[25]。當PBIAS<0時表示模擬值較實測值偏大，PBIAS>0時表示模擬值較實測值偏小。PBIAS的計算結果如表1所示。

表1 灞河流域水環(huán)境模型百分比偏差系數(shù)(PBIAS)的計算結果Table 1 Percent bias of water environment model at Ba River basin %

表1顯示，除了S1點位的NH3-N(-25.41%)和S3點位的TP(25.06%)外，各項指標的PBIAS絕對值均小于25%，這表明模型模擬值與實測值擬合趨勢良好[26]，已建立的模型較好地反映了研究區(qū)水質(zhì)指標的變化情況，可用于該區(qū)域不同工況下水質(zhì)指標的模擬計算。

4.2 灞河口水質(zhì)分析

根據(jù)《陜西省水功能區(qū)劃》[27]規(guī)定，研究區(qū)水質(zhì)目標為地表水Ⅳ類(GB 3838-2002)，相對渭河而言，灞河口水質(zhì)需滿足Ⅳ類目標。本研究計算域未包含灞河1號壩下游至灞河口段，考慮到未納入計算域的河段長度僅約為1.7 km，且該河段基本無外源污染擾動，因此選取該區(qū)域出口斷面水質(zhì)代表灞河口水質(zhì)進行分析。不同工況下灞河口水質(zhì)指標變化趨勢如圖3所示。

圖3 不同工況下灞河口斷面DO、COD、NH3-N和TP質(zhì)量濃度的變化Fig.3 Variation of DO,COD,NH3-N and TP concentrations at the Bahekou section for different scenarios

圖3表明，TP與DO在研究時段內(nèi)均滿足Ⅳ類水質(zhì)標準，COD與NH3-N部分時段未達Ⅳ類水質(zhì)標準，這與文獻[28]的論述一致。COD和NH3-N的超標時段主要集中在1～31 d(7月至8月初)。

對灞河口斷面COD與NH3-N 2項指標超標情況進行統(tǒng)計，結果見表2。

表2 不同工況下灞河口斷面COD和NH3-N質(zhì)量濃度超標情況統(tǒng)計結果Table 2 Statistics of COD and NH3-N concentrations exceeding standard at the Bahekou section for different scenarios

從表2可以觀察到，2012-07-10-10-15共計98 d內(nèi)，豐水年與平水年COD超標時間少于基準年與枯水年。豐水年與平水年NH3-N超標時間分別為22和27 d，也小于基準年與枯水年。豐水年COD與NH3-N最大超標率大于平水年與基準年，這是因為在水質(zhì)邊界條件相同時，由于豐水年水動力邊界的流量在部分時段小于基準年與平水年，且該時段點源流量相對較入口流量大。基準年NH3-N的超標時間明顯大于COD。分別對基準年計算時段內(nèi)研究區(qū)域的COD與NH3-N取算術平均值，由這2項指標的均值明確的水質(zhì)類型與西安市生態(tài)環(huán)境局公布的《2012年第三季度環(huán)境質(zhì)量監(jiān)測分析》中的水質(zhì)類型一致。

4.3 水質(zhì)類型空間特性分析

研究水質(zhì)空間分布特性是識別水域重點保護區(qū)的重要途徑之一。鑒于上述原因，提取水環(huán)境模型計算結果中DO、COD、NH3-N和TP共4項指標后，對各指標在時間尺度(2012-07-10-10-15)上取算術平均值，然后利用前期訓練好的BP神經(jīng)網(wǎng)絡對研究區(qū)綜合水質(zhì)類型進行計算，該計算結果稱為時均綜合水質(zhì)類型。不同工況下研究區(qū)時均綜合水質(zhì)類型的空間分布如圖4所示。從圖4可以看出，浐河和灞河在交匯前，灞河水質(zhì)整體優(yōu)于浐河水質(zhì)，灞河水質(zhì)類型主要為Ⅱ類，部分區(qū)域為Ⅲ類；浐河水質(zhì)類型主要為Ⅳ類，部分區(qū)域為Ⅲ類。灞河3號壩下游水質(zhì)類型除豐水年為Ⅲ類外，其余工況下均為Ⅳ類。不同工況下，盡管2012年7-10月灞河4號壩下游部分區(qū)域時均綜合水質(zhì)類型超過Ⅳ類，鑒于其影響范圍較小，且該范圍內(nèi)水質(zhì)下降的主要原因在于河道左岸點源污染的影響，西安浐灞生態(tài)區(qū)水生態(tài)系統(tǒng)保護與修試點工作中的污水治理與中水回用工程實施后，已消除了該范圍內(nèi)的點源影響，因此認為研究區(qū)時均綜合水質(zhì)類型在空間范圍內(nèi)于基準年、豐水年、平水年和枯水年基本都滿足Ⅳ類水質(zhì)目標要求。

圖4 不同工況下灞河時均綜合水質(zhì)類型的空間分布Fig.4 Spatial distribution of time-average comprehensive water quality at Ba River

由圖4可知，受浐河水體擾動，灞河水體在與浐河水體交匯后，交匯區(qū)水質(zhì)類型出現(xiàn)平面梯度分區(qū)，水質(zhì)類型梯度分區(qū)從左岸到右岸依次為Ⅳ類-Ⅲ類-Ⅱ類。交匯區(qū)水質(zhì)類型分布特征主要受灞河與浐河交匯前的動量比率以及兩水體各項水質(zhì)指標影響。交匯區(qū)水質(zhì)類型分布特征在不同工況下表現(xiàn)出明顯差異，灞河入交匯區(qū)前，豐水年動量大于基準年，因此豐水年條件下交匯區(qū)Ⅲ類水的帶狀面積較基準年條件下有所增加，且Ⅲ類水的條帶明顯向左岸推進?？菟陾l件下灞河入交匯區(qū)前動量比率小于基準年條件下，交匯區(qū)Ⅲ類水條帶明顯向右岸推進。

為分析不同工況下研究區(qū)內(nèi)水質(zhì)類型的主導性指標，提取研究區(qū)內(nèi)DO、COD、NH3-N和TP共4項指標，對各指標在時間尺度(2012-07-10-10-15)上取算術平均值，分別對不同工況進行PCA分析。由于KMO(Kaiser-Meyer-Olkin)值均大于0.5且Bartlett顯著性概率均為0.000，滿足PCA分析條件。根據(jù)特征值λi>1的原則提取的1個主成分反映了原始4個指標57%以上的信息。該結果表明，不同工況下COD與NH3-N對第一主成分的貢獻率為63%～64%，可以認為不同工況下研究區(qū)時均綜合水質(zhì)類型的主導性指標主要為COD與NH3-N。

5 結 論

將DHI MIKE21FM中HD模塊與Ecolab模塊耦合后構建了適用于灞河流域的水環(huán)境模型，研究了灞河口水質(zhì)指標的變化情況，基于BP神經(jīng)網(wǎng)絡算法結合模擬結果對研究區(qū)時均綜合水質(zhì)類型進行了評價，借助PCA識別了研究區(qū)時均綜合水質(zhì)類型的主導性指標，得到的研究結論如下：

1)灞河口水質(zhì)超標因子主要為COD與NH3-N，2012年7-10月NH3-N的超標時間與最大超標率明顯大于COD。豐水年與平水年COD與NH3-N的超標時間較基準年明顯下降，但最大超標率增加。

2)灞河和浐河交匯前，灞河水質(zhì)整體優(yōu)于浐河，灞河水質(zhì)類型主要為Ⅱ類，部分區(qū)域為Ⅲ類；浐河水質(zhì)類型主要為Ⅳ類，部分區(qū)域為Ⅲ類。除豐水年灞河3號壩下游水質(zhì)類型為Ⅲ類外，基準年、平水年與枯水年該區(qū)域水質(zhì)類型均為Ⅳ類。研究區(qū)時均綜合水質(zhì)類型基本滿足Ⅳ類水質(zhì)目標，水質(zhì)類型主要受COD與NH3-N 2項指標影響。

3)所構建的水環(huán)境模型，經(jīng)驗證后具有較高的精度，可適用于類似區(qū)域的水質(zhì)預測及評價。

需要指出的是，由于水質(zhì)監(jiān)測數(shù)據(jù)時長限制，基準年僅計算了7-10月水質(zhì)指標的變化情況，計算時段并未覆蓋1個水文年。因此，本研究得出的結論主要適用于灞河流域雨季的水環(huán)境特征，時均綜合水質(zhì)類型并不能代表全年平均的綜合水質(zhì)類型。水環(huán)境數(shù)值模擬是一個多情境與多維度的過程，本研究在水質(zhì)模型構建過程中，以實測數(shù)據(jù)表征了上游的點源與非點源的污染特征，但同時也忽略了下游非點源輸入，后續(xù)還需進一步研究。