基于線性規劃模型的水環境容量計算研究

洪夕媛，孫明東，程全國

(1.沈陽大學環境學院，沈陽 110044；2.中國環境科學研究院，北京 100000)

0 引 言

隨著中國經濟的持續發展以及城市化進程的不斷加快，水體受污染的程度也越來越嚴重，水環境容量保護領域的研究熱點逐漸轉向污染物排放總量控制和水環境容量的計算，在控制總污染物的情況下，對水環境容量的研究成為當前研究的熱點和重要課題。

永定河是海河流域面積最大、最重要的水系，同時也是京津冀區域重要的水源涵養區。但目前永定河流域存在水資源過渡開發、環境承載力差、河段斷流等突出問題，嚴重制約了地區的經濟社會健康發展。為了解決永定河流域的水環境污染問題，明確其水環境容量可以為河流水質管理提供依據，對實現水環境規劃治理有著重要作用[1]。

水環境容量是指特定的水體在給定的設計水文條件、水功能區目標水質、排污口位置及排放方式條件下，滿足水域某一水質目標的污染物在排污口的最大排放量，也可以定量解釋為水體對該污染物的承載力[2]。目前應用最廣泛的水環境容量計算方法是解析公式法和線性規劃法[3]，線性規劃法較解析法來說對排污口的總量控制更具有針對性，并且自動化程度高、精度高、對邊界條件及設計條件具備更強的適應性。董飛等基于綜合水動力水質模型，提出了水質響應系數矩陣的高效計算方法；荊曉海等[5]通過線性規劃模型對北運河水環境容量進行優化分配，并提供兩種排污優化方案。

基于2017年河北省第二次污染普查數據顯示，張家口市的直排工業企業還未清零，沿岸污水處理廠的數量也比較多，作為官廳水庫上游來水河段，改善其水質、明確水體的負載現狀及納污能力尤為重要。因此文章采用線性規劃模型對永定河上游張家口段的水環境容量進行優化求解，明確各排污口的最大允許納污量，為后續的污染物總量控制提供技術支持。

1 研究區概況

永定河內蒙古高原和山西高原，由西至東流經內蒙古、山西、河北、北京、天津5個省(直轄市)，其中河北省流經張家口市和廊坊市。文章研究范圍為張家口市境內永定河上游流域，流域面積為17713 km2，境內主要包含洋河、桑干河、清水河等。地勢西北高東南低，陰山山脈橫貫中部，將張家口劃分為壩上高原區和壩下低中山盆地，流域內上游海拔高1500m以上，屬溫帶大陸性季風氣候，平均降水量約為400mm/a，多年平均氣溫8.5℃。

文章將東洋河-洋河作為干流，支流設為清水河和桑干河進行模型的搭建工作，水質目標依據《河北省水環境功能區劃》。研究區示意圖1，見圖1。

圖1 研究區示意圖

2 研究方法

2.1 水動力水質模型的構建

2.1.1 模型原理

QUAL2K模型為綜合一維水動力水質模型，主要模擬河流干流。該模型適用于混合良好的枝狀河流中的BOD、N、P等污染物的遷移轉化模擬[6]。永定河流域張家口段屬于中小型河流，因此采用該模型進行流域的水質模擬研究。該模型以計算單元為最小模擬單位，將穩態條件下的流量平衡方程應用于每個計算單元，流量平衡方程如下：

Qi=Qi-1+Qin,i-Qout,i

(1)

式中：Qi是單元i到單元i+1的出流流量，m3·d-1；Qi-1是單元i-1到單元i的入流流量，m3·d-1；Qin，i是指點源和非點源直接匯入單元i的總流量，m3·d-1；Qout，i是指單元i的點源和非點源直接流出的總流量，m3·d-1。

模型的水質方程是一維平流-擴散物質遷移方程，考慮了水質組分如水的對流擴散、稀釋和外源輸入或變化的響應關系，對于任意一種水質組分，有：

(2)

式中：C為污染物濃度，mg·L-1；x為所計算的距離，m；t為計算的時間，s；Ax為距離為x處的河流斷面面積，m2；DL為縱向彌散系數，m2·s-1；u為平均流速，m·s-1；s為組分外部的點源和非點源污染物濃度，mg·L-1；V為計算單元體積，m3。方程右側4項分別代表擴散、平流、組分反應和組分的外部源和匯，并采用有限差分法進行求解。

2.1.2 水質評價

文章采用改進的內梅羅污染指數法進行流域的水質評價。傳統的內梅羅指數法突出了污染指數最大的污染因子對環境的影響而忽略了一些實測濃度較低(如重金屬離子、總磷等)卻對水環境有很大影響的污染因子，而改進的內梅羅指數法為克服這種缺陷，通過賦予將每個污染因子賦予權重量化每種因子對水體的影響程度，使得評價結果更為準確[7]，計算公式為：

(3)

(4)

(5)

(6)

(7)

式中：wi為各污染因子的權重；smax為m個污染因子中計算標準的最大值，如所選的6個污染因子中，COD的三類標準最大，則smax=20；ri為第i種污染因子的相關性比值；si為第i個污染因子的計算標準；Fi為第i個污染因子的污染指數；ci為第i個污染因子的實測濃度；si為第i個污染因子的計算標準；Fw為權重值最大的污染因子的Fi值；Fi,max為第i個污染因子中污染指數最大值；P′即為改進的內梅羅污染指數。

通過對研究區內9控制斷面的水質數據對各河流水質進行評價，改進的內梅河綜合污染指數評價結果，見表1。通過水質評價的結果顯示，該流域的水質基本可以滿足水環境功能區要求，桑干河水質要優于洋河和清水河，洋河水質最差；9個監測斷面中，位于下游邊界的“八號橋”控制斷面水質最差，該斷面為洋河與桑干河匯合后的監測斷面同樣也作為進入官廳水庫的水質控制斷面，其改進的內梅羅指數為0.995已接近四類水質等級，因此對于洋河水質的改善還有待進一步提高。

表1 改進的內梅羅綜合污染指數評價結果

結合“十三五”期間水污染物總量控制指標為CODcr和NH3-N以及永定河污染現狀，本研究選擇CODcr和NH3-N作為水環境容量核算的評價因子。

2.1.3 水質參數的確定

QUAL2K模型水質參數的選擇包括COD氧化速率kdc和氨氮硝化速率kna，模型參數的率定選用2017年1-12月的水質數據，模型驗證選用2018年1-12月實測水質數據。最終驗證結果均在20%以內，表明該模型在永定河流域張家口的適用性良好。最終確定COD氧化速率kdc的取值范圍為0.1-0.35 d-1，氨氮硝化速率kna的取值范圍為0.1-0.6 d-1。

2.2 控制斷面選取及排污口概化

明確控制斷面是計算響應系數的前提，由于不是所有河流都有數量完整且距離合適的省控或國控斷面，因此需要在河道中合理選取控制斷面，這對進行污染源與控制斷面響應系數的計算起著至關重要的作用。

控制斷面的選擇遵循以下原則[8]：①水環境功能區的上界和下界、支流匯入口、排污口等造成河道流量水質發生突變的點；②當斷面距離較大時，為減少斷面的各參數的差距，可適當插入一個或幾個斷面；③直排企業概化后的排污口也可以作為控制斷面。根據2017年河北省第2次污染普查數據顯示，洋河上游主要為面源污染，下游排污口主要是城鎮污水處理廠、各支流入河口以及工業企業直排口，支流清水河和桑干河面源污染較集中。文章將流域內一級支流視為排污口，面源污染比較集中的沿江村鎮概化為一個排污口，對于污水處理廠分布較密集的區域可將其概化成一個集中的排污口，合并后的排污口與所在河段上游斷面的距離，計算公式為：

(8)

式中：L為概化后集中排放口到控制斷面的距離，km；Qi為第i個排污口(支流口)水量，m3/s；Ci為第i個排污口排放污染物的濃度，mg/L；Li為第i個排污口到下一個控制斷面的距離，km。

根據張家口市永定河流域的污染源排放情況，將主要區縣及其重要支流作為其排放排污口，排污口概化細則，見表2；河流概化模型示意圖，見圖2。

表2 排污口概化細則

2.3 水環境容量計算

2.3.1 線性規劃模型

線性規劃模型計算水環境容量的優點在于，能夠考慮到各個排污口的分布，方便對各個排污口進行總量分配研究。該方法計算水環境容量的基本思路是：①在干流建立以污染物排放量為變量的控制斷面；②目標函數設為所有排污口排放量總量最大；③響應濃度≤實際約束濃度設置為約束方程；④各排污口的最大允許負荷量為決策變量，通過Execl中規劃求解工具進行線性規劃求解，最后通過對所求區段內的各污染源允許排放量加和，即可得到全流域的最大允許負荷量即為水環境容量。

Xj≥0,j=1,2,…,N

(9)

式中：決策變量xj為第j污染源的排放量，g/s；aij為第j污染源對第i控制斷面的響應系數，可由水質模型計算得到；ci為控制斷面的實際約束濃度，mg/L；aijxij稱為響應濃度，mg/L。

圖2 河流概化模型示意圖

2.3.2 響應系數矩陣

污染負荷的輸入與受納水體的水質響應關系反映了一定排污負荷條件下對河流控制斷面水質貢獻率，可以通過響應系數加以表達，如對于一個M×N的矩陣A，其中Aij表示控制斷面i對點源j的響應系數，響應系數矩陣為各斷面與概化后排污口的響應系數構成的矩陣。文章選用經過實測值驗證過的QUAL2K一維水動力水質模型計算污染源與斷面水質的響應系數，最終得到排污口排放量與斷面水質的響應系數矩陣結果，CODcr排放量與斷面水質響應關系矩陣見表2和NH3-N排放量與斷面水質響應關系矩陣見表3。

2.3.3 設計流量

對于設計流量的選取，關乎到水環境容量的大小，但至今為止沒有絕對的規定。《水域納污能力計算規程》中規定，河流的設計流量應采用近10a最枯月平均流量或90%保證率最枯月平均流量。根據中國環境規劃研究院給出的水環境容量計算設計流量的分析，對于北方河流，90%保證率下的流量往往為零，可選擇近10a最枯季75%保證率或選擇50%保證率同等條件進行對比計算[9]。

表2 CODcr排放量與斷面水質響應關系矩陣

表3 NH3-N排放量與斷面水質響應關系矩陣

永定河流域為典型的北方地區河流，枯水期流量小甚至斷流，河流徑流量的季節性變化較大[10]。文章選取2002-2017年≠0的最小月逐日平均流量作為樣本，采用水文保證率法，借助Execl工具，求得水文變量設計值與頻率之間的定量關系，確定干流及其主要支流75%保證率條件下的設計流量，其中友誼水庫并非全年放水，因此設計流量采用近10a最枯月平均流量，干流及支流設計流量，見表4。

表4 干流及支流設計流量 m3/s

4)水動力學參數：水動力學參數的計算目的是明確河道的水力特征，根據流域內水文站實測數據確定河流流速、水深、河寬與流量的關系系數，其經驗公式如下：

水深—流量關系：H=a1Qb1

河寬—流量關系：B=a2Qb2

流速—流量關系：U=a3Qb2

式中：a、b為經驗系數，a1、a2和a3三者的乘積為1，b1、b2和b3三者之和為1。通過擬合得到干流及支流的水動力參數，各河流水動力學參數，見表5。

3 結果與分析

通過模型建立的污染源與水質的響應關系，基于河北省水環境功能區劃設定的水質目標，計算各排污口允許納污量，再結合各斷面污染物的現狀負荷計算各區域的削減量。水環境容量計算結果，見表6。

表5 各河流水動力學參數

表6 水環境容量計算結果 t/a，%

由表6可知，永定河流域張家口段CODcr水環境容量為4724.7 t/a，NH3-N水環境容量為377.34 t/a，研究區域內剩余水環境容量為負值，表明河流的納污能力已經超過最大允許納污量。桑干河雖水環境容量最大，但其污染負荷量也是最多的，通過CODcr、NH3-N剩余水環境容量的計算結果可以看出，僅有上游尚義縣和萬全縣存在剩余水環境容量，因為二者作為上游來水，水質較好，因此水環境容量較大，不需要進行削減；自清水河入洋河之后，剩余水環境容量均為負值，表明現狀負荷均超過水環境容量。CODcr及NH3-N平均削減量為25.8%、6.3%；其中，宣化縣削減任務最為繁重，CODcr和NH3-N分別需要削減906.05 t/a、105.12 t/a，分別占總削減任務的76.09%、83.27%。

通過對永定河流域張家口段的污染負荷核定，CODcr和NH3-N污染負荷入河量分別為6575.33 t/a、418 t/a，其中點源CODcr負荷為2902.81 t/a，NH3-N負荷為174.33 t/a，分別占總入河量的44.1%和41.7%，由于直排企業僅剩兩家，因此點源污染的主要來源為城鎮生活源，且污水處理廠多集中于宣化縣河流沿岸，該縣CODcr和NH3-N點源入河量占比分別為66.76%和53.06%。研究區內非點源污染中CODcr的入河量為3675.52 t/a，其中分散畜禽養殖為CODcr非點源污染的主要來源，占CODcr非點源污染總入河量的62.8%；NH3-N非點源入河量為244.14 t/a，其中農田徑流污染為NH3-N非點源污染的主要來源，占NH3-N非點源污染總入河量的73.6%。對于清水河以及桑干河污染主要以農田徑流及分散畜禽養殖污染為主。

綜上，通過對污染負荷的分析，能夠更加明確各排污口的污染物來源方向，接下來便可采用更有針對性的減排措施進行污染物的削減措施。對于工業排污口集中的宣化縣，應著重提升污水處理能力、強化污水的收集以及提高中水使用率；清水河和桑干河，應對畜禽養殖污染實施相應的管理政策，將畜禽糞便資源化，與種植業密切結合發展生態農業、有機農業，對于農田徑流的防治可以通過控制施肥量、調整施肥時間和方法、施用有機肥等措施進行改善。

4 結 論

通過QUAL2K水動力水質模型計算污染源與控制斷面的水質響應關系，并利用線性規模型計算得到永定河流域張家口段洋河干流各排污口的最大允許納污量，得到以下結論：

1)基于改進的內梅羅污染指數法對水質進行的評價分析，流域內主要河流水質現狀依次為桑干河>(優于)清水河>(優于)洋河，下游邊界“八號橋”斷面水質最差。因此為達到水質目標，永定河上游流域應首要改善洋河水質。

2)通過對洋河干流進行QUAL2K水動力水質模型的構建，最終驗證結果均在20%以內，表明該模型在永定河流域張家口段的適用性良好。

3)CODcr及NH3-N的水環境容量分別為4724 t/a、377.34 t/a。自清水河入洋河之后，河段均無剩余水環境容量，其中宣化縣超標最為嚴重，清水河次之，桑干河的納污能力與現狀污染負荷量大致相當；流域內CODcr及NH3-N污染負荷平均削減量為25.8%、6.3%。