基于線性規(guī)劃模型的水環(huán)境容量計(jì)算研究

洪夕媛，孫明東，程全國(guó)

(1.沈陽大學(xué)環(huán)境學(xué)院，沈陽 110044；2.中國(guó)環(huán)境科學(xué)研究院，北京 100000)

0 引 言

隨著中國(guó)經(jīng)濟(jì)的持續(xù)發(fā)展以及城市化進(jìn)程的不斷加快，水體受污染的程度也越來越嚴(yán)重，水環(huán)境容量保護(hù)領(lǐng)域的研究熱點(diǎn)逐漸轉(zhuǎn)向污染物排放總量控制和水環(huán)境容量的計(jì)算，在控制總污染物的情況下，對(duì)水環(huán)境容量的研究成為當(dāng)前研究的熱點(diǎn)和重要課題。

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

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

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

1 研究區(qū)概況

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

文章將東洋河-洋河作為干流，支流設(shè)為清水河和桑干河進(jìn)行模型的搭建工作，水質(zhì)目標(biāo)依據(jù)《河北省水環(huán)境功能區(qū)劃》。研究區(qū)示意圖1，見圖1。

圖1 研究區(qū)示意圖

2 研究方法

2.1 水動(dòng)力水質(zhì)模型的構(gòu)建

2.1.1 模型原理

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

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是指點(diǎn)源和非點(diǎn)源直接匯入單元i的總流量，m3·d-1；Qout，i是指單元i的點(diǎn)源和非點(diǎn)源直接流出的總流量，m3·d-1。

模型的水質(zhì)方程是一維平流-擴(kuò)散物質(zhì)遷移方程，考慮了水質(zhì)組分如水的對(duì)流擴(kuò)散、稀釋和外源輸入或變化的響應(yīng)關(guān)系，對(duì)于任意一種水質(zhì)組分，有：

(2)

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

2.1.2 水質(zhì)評(píng)價(jià)

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

(3)

(4)

(5)

(6)

(7)

式中：wi為各污染因子的權(quán)重；smax為m個(gè)污染因子中計(jì)算標(biāo)準(zhǔn)的最大值，如所選的6個(gè)污染因子中，COD的三類標(biāo)準(zhǔn)最大，則smax=20；ri為第i種污染因子的相關(guān)性比值；si為第i個(gè)污染因子的計(jì)算標(biāo)準(zhǔn)；Fi為第i個(gè)污染因子的污染指數(shù)；ci為第i個(gè)污染因子的實(shí)測(cè)濃度；si為第i個(gè)污染因子的計(jì)算標(biāo)準(zhǔn)；Fw為權(quán)重值最大的污染因子的Fi值；Fi,max為第i個(gè)污染因子中污染指數(shù)最大值；P′即為改進(jìn)的內(nèi)梅羅污染指數(shù)。

通過對(duì)研究區(qū)內(nèi)9控制斷面的水質(zhì)數(shù)據(jù)對(duì)各河流水質(zhì)進(jìn)行評(píng)價(jià)，改進(jìn)的內(nèi)梅河綜合污染指數(shù)評(píng)價(jià)結(jié)果，見表1。通過水質(zhì)評(píng)價(jià)的結(jié)果顯示，該流域的水質(zhì)基本可以滿足水環(huán)境功能區(qū)要求，桑干河水質(zhì)要優(yōu)于洋河和清水河，洋河水質(zhì)最差；9個(gè)監(jiān)測(cè)斷面中，位于下游邊界的“八號(hào)橋”控制斷面水質(zhì)最差，該斷面為洋河與桑干河匯合后的監(jiān)測(cè)斷面同樣也作為進(jìn)入官?gòu)d水庫(kù)的水質(zhì)控制斷面，其改進(jìn)的內(nèi)梅羅指數(shù)為0.995已接近四類水質(zhì)等級(jí)，因此對(duì)于洋河水質(zhì)的改善還有待進(jìn)一步提高。

表1 改進(jìn)的內(nèi)梅羅綜合污染指數(shù)評(píng)價(jià)結(jié)果

結(jié)合“十三五”期間水污染物總量控制指標(biāo)為CODcr和NH3-N以及永定河污染現(xiàn)狀，本研究選擇CODcr和NH3-N作為水環(huán)境容量核算的評(píng)價(jià)因子。

2.1.3 水質(zhì)參數(shù)的確定

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

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

明確控制斷面是計(jì)算響應(yīng)系數(shù)的前提，由于不是所有河流都有數(shù)量完整且距離合適的省控或國(guó)控?cái)嗝?，因此需要在河道中合理選取控制斷面，這對(duì)進(jìn)行污染源與控制斷面響應(yīng)系數(shù)的計(jì)算起著至關(guān)重要的作用。

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

(8)

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

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

表2 排污口概化細(xì)則

2.3 水環(huán)境容量計(jì)算

2.3.1 線性規(guī)劃模型

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

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

(9)

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

圖2 河流概化模型示意圖

2.3.2 響應(yīng)系數(shù)矩陣

污染負(fù)荷的輸入與受納水體的水質(zhì)響應(yīng)關(guān)系反映了一定排污負(fù)荷條件下對(duì)河流控制斷面水質(zhì)貢獻(xiàn)率，可以通過響應(yīng)系數(shù)加以表達(dá)，如對(duì)于一個(gè)M×N的矩陣A，其中Aij表示控制斷面i對(duì)點(diǎn)源j的響應(yīng)系數(shù)，響應(yīng)系數(shù)矩陣為各斷面與概化后排污口的響應(yīng)系數(shù)構(gòu)成的矩陣。文章選用經(jīng)過實(shí)測(cè)值驗(yàn)證過的QUAL2K一維水動(dòng)力水質(zhì)模型計(jì)算污染源與斷面水質(zhì)的響應(yīng)系數(shù)，最終得到排污口排放量與斷面水質(zhì)的響應(yīng)系數(shù)矩陣結(jié)果，CODcr排放量與斷面水質(zhì)響應(yīng)關(guān)系矩陣見表2和NH3-N排放量與斷面水質(zhì)響應(yīng)關(guān)系矩陣見表3。

2.3.3 設(shè)計(jì)流量

對(duì)于設(shè)計(jì)流量的選取，關(guān)乎到水環(huán)境容量的大小，但至今為止沒有絕對(duì)的規(guī)定?！端蚣{污能力計(jì)算規(guī)程》中規(guī)定，河流的設(shè)計(jì)流量應(yīng)采用近10a最枯月平均流量或90%保證率最枯月平均流量。根據(jù)中國(guó)環(huán)境規(guī)劃研究院給出的水環(huán)境容量計(jì)算設(shè)計(jì)流量的分析，對(duì)于北方河流，90%保證率下的流量往往為零，可選擇近10a最枯季75%保證率或選擇50%保證率同等條件進(jìn)行對(duì)比計(jì)算[9]。

表2 CODcr排放量與斷面水質(zhì)響應(yīng)關(guān)系矩陣

表3 NH3-N排放量與斷面水質(zhì)響應(yīng)關(guān)系矩陣

永定河流域?yàn)榈湫偷谋狈降貐^(qū)河流，枯水期流量小甚至斷流，河流徑流量的季節(jié)性變化較大[10]。文章選取2002-2017年≠0的最小月逐日平均流量作為樣本，采用水文保證率法，借助Execl工具，求得水文變量設(shè)計(jì)值與頻率之間的定量關(guān)系，確定干流及其主要支流75%保證率條件下的設(shè)計(jì)流量，其中友誼水庫(kù)并非全年放水，因此設(shè)計(jì)流量采用近10a最枯月平均流量，干流及支流設(shè)計(jì)流量，見表4。

表4 干流及支流設(shè)計(jì)流量 m3/s

4)水動(dòng)力學(xué)參數(shù)：水動(dòng)力學(xué)參數(shù)的計(jì)算目的是明確河道的水力特征，根據(jù)流域內(nèi)水文站實(shí)測(cè)數(shù)據(jù)確定河流流速、水深、河寬與流量的關(guān)系系數(shù)，其經(jīng)驗(yàn)公式如下：

水深—流量關(guān)系：H=a1Qb1

河寬—流量關(guān)系：B=a2Qb2

流速—流量關(guān)系：U=a3Qb2

式中：a、b為經(jīng)驗(yàn)系數(shù)，a1、a2和a3三者的乘積為1，b1、b2和b3三者之和為1。通過擬合得到干流及支流的水動(dòng)力參數(shù)，各河流水動(dòng)力學(xué)參數(shù)，見表5。

3 結(jié)果與分析

通過模型建立的污染源與水質(zhì)的響應(yīng)關(guān)系，基于河北省水環(huán)境功能區(qū)劃設(shè)定的水質(zhì)目標(biāo)，計(jì)算各排污口允許納污量，再結(jié)合各斷面污染物的現(xiàn)狀負(fù)荷計(jì)算各區(qū)域的削減量。水環(huán)境容量計(jì)算結(jié)果，見表6。

表5 各河流水動(dòng)力學(xué)參數(shù)

表6 水環(huán)境容量計(jì)算結(jié)果 t/a，%

由表6可知，永定河流域張家口段CODcr水環(huán)境容量為4724.7 t/a，NH3-N水環(huán)境容量為377.34 t/a，研究區(qū)域內(nèi)剩余水環(huán)境容量為負(fù)值，表明河流的納污能力已經(jīng)超過最大允許納污量。桑干河雖水環(huán)境容量最大，但其污染負(fù)荷量也是最多的，通過CODcr、NH3-N剩余水環(huán)境容量的計(jì)算結(jié)果可以看出，僅有上游尚義縣和萬全縣存在剩余水環(huán)境容量，因?yàn)槎咦鳛樯嫌蝸硭?，水質(zhì)較好，因此水環(huán)境容量較大，不需要進(jìn)行削減；自清水河入洋河之后，剩余水環(huán)境容量均為負(fù)值，表明現(xiàn)狀負(fù)荷均超過水環(huán)境容量。CODcr及NH3-N平均削減量為25.8%、6.3%；其中，宣化縣削減任務(wù)最為繁重，CODcr和NH3-N分別需要削減906.05 t/a、105.12 t/a，分別占總削減任務(wù)的76.09%、83.27%。

通過對(duì)永定河流域張家口段的污染負(fù)荷核定，CODcr和NH3-N污染負(fù)荷入河量分別為6575.33 t/a、418 t/a，其中點(diǎn)源CODcr負(fù)荷為2902.81 t/a，NH3-N負(fù)荷為174.33 t/a，分別占總?cè)牒恿康?4.1%和41.7%，由于直排企業(yè)僅剩兩家，因此點(diǎn)源污染的主要來源為城鎮(zhèn)生活源，且污水處理廠多集中于宣化縣河流沿岸，該縣CODcr和NH3-N點(diǎn)源入河量占比分別為66.76%和53.06%。研究區(qū)內(nèi)非點(diǎn)源污染中CODcr的入河量為3675.52 t/a，其中分散畜禽養(yǎng)殖為CODcr非點(diǎn)源污染的主要來源，占CODcr非點(diǎn)源污染總?cè)牒恿康?2.8%；NH3-N非點(diǎn)源入河量為244.14 t/a，其中農(nóng)田徑流污染為NH3-N非點(diǎn)源污染的主要來源，占NH3-N非點(diǎn)源污染總?cè)牒恿康?3.6%。對(duì)于清水河以及桑干河污染主要以農(nóng)田徑流及分散畜禽養(yǎng)殖污染為主。

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

4 結(jié) 論

通過QUAL2K水動(dòng)力水質(zhì)模型計(jì)算污染源與控制斷面的水質(zhì)響應(yīng)關(guān)系，并利用線性規(guī)模型計(jì)算得到永定河流域張家口段洋河干流各排污口的最大允許納污量，得到以下結(jié)論：

1)基于改進(jìn)的內(nèi)梅羅污染指數(shù)法對(duì)水質(zhì)進(jìn)行的評(píng)價(jià)分析，流域內(nèi)主要河流水質(zhì)現(xiàn)狀依次為桑干河>(優(yōu)于)清水河>(優(yōu)于)洋河，下游邊界“八號(hào)橋”斷面水質(zhì)最差。因此為達(dá)到水質(zhì)目標(biāo)，永定河上游流域應(yīng)首要改善洋河水質(zhì)。

2)通過對(duì)洋河干流進(jìn)行QUAL2K水動(dòng)力水質(zhì)模型的構(gòu)建，最終驗(yàn)證結(jié)果均在20%以內(nèi)，表明該模型在永定河流域張家口段的適用性良好。

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