基于線性規(guī)劃模型的河流水環(huán)境容量分配研究

荊海曉， 李小寶， 房懷陽(yáng)， 劉長(zhǎng)根

(1.西安理工大學(xué) 省部共建西北旱區(qū)生態(tài)水利國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710048；2.環(huán)境保護(hù)部華南環(huán)境科學(xué)研究所，廣東 廣州 510535； 3.天津大學(xué) 機(jī)械工程學(xué)院，天津 300350)

1 研究背景

水環(huán)境容量是指水體在規(guī)定的環(huán)境目標(biāo)下所能容納的最大污染物量[1]，其反映了污染物在環(huán)境中的遷移及轉(zhuǎn)化規(guī)律，是衡量水體的水環(huán)境承載能力的一個(gè)重要指標(biāo)。一般情況下，某一水體的水環(huán)境容量的大小主要取決于水體特性、水質(zhì)目標(biāo)和污染物的特性等因素[2]。通過研究水體的水環(huán)境容量來實(shí)現(xiàn)對(duì)水體的水環(huán)境質(zhì)量的控制已經(jīng)成為現(xiàn)代水環(huán)境保護(hù)和規(guī)劃的重要依據(jù)。目前，對(duì)于水環(huán)境容量計(jì)算方法以及在不同水域的應(yīng)用，國(guó)內(nèi)外學(xué)者已進(jìn)行了很多研究，并取得了一系列的成果[1,3-6]。李珊等[1]在一維穩(wěn)態(tài)水質(zhì)模型的基礎(chǔ)上推導(dǎo)出了太子河水環(huán)境容量計(jì)算公式，對(duì)不同頻率年，該河段COD、氨氮的環(huán)境容量進(jìn)行了計(jì)算評(píng)價(jià)，并提出相應(yīng)的排污措施。王子軒等[4]通過建立淀山湖流域一維平原河網(wǎng)水量水質(zhì)模型，對(duì)該流域的COD、氨氮環(huán)境容量進(jìn)行了計(jì)算評(píng)估，設(shè)計(jì)了相應(yīng)的斷面水質(zhì)達(dá)標(biāo)方案。熊鴻斌等[5]基于MIKE11模型研究了引江濟(jì)淮工程渦河段水環(huán)境容量。韓梓流等[6]采用一維水量水質(zhì)數(shù)學(xué)模型研究了京杭運(yùn)河五牧河斷面達(dá)標(biāo)時(shí)各概化排污口的允許排污量，得到了該研究區(qū)域的水環(huán)境容量，為改善京杭運(yùn)河水環(huán)境質(zhì)量提供了科學(xué)依據(jù)。閻非等[7]研究了綜合考慮多個(gè)點(diǎn)源排污口、支流、取水口以及面源污染，并結(jié)合排污口權(quán)重，提出了一種實(shí)用化的水環(huán)境容量計(jì)算方法和排污控制方案。李兵等[8]通過對(duì)蘇南運(yùn)河蘇錫常段沿岸的工業(yè)污染源、農(nóng)業(yè)污染源和城鎮(zhèn)生活污水排污分布情況進(jìn)行調(diào)查，對(duì)污染源進(jìn)行了評(píng)價(jià)，并根據(jù)水質(zhì)目標(biāo)，計(jì)算出了各段水體的水環(huán)境容量、排污控制量和削減率。董一博等[9]采用一維水動(dòng)力水質(zhì)模型并結(jié)合線性規(guī)劃模型研究了太子河遼河段水環(huán)境容量。

北運(yùn)河是海河北系四大河流之一，發(fā)源于北京市昌平縣燕山南麓，先后流經(jīng)北京市海淀、昌平、順義、通州各區(qū)、河北省廊坊地區(qū)、天津市武清區(qū)、北辰區(qū)和紅橋區(qū)。其流域面積為6 166 km2，其中山區(qū)面積為952 km2，平原面積為5 214 km2。多年平均徑流量為4.81×108m3，其中山區(qū)年均凈流量為1.29×108m3、平原為3.52×108m3。通州區(qū)北關(guān)閘以上稱溫榆河，北關(guān)閘以下稱北運(yùn)河，沿途的支流包括通惠河、涼水河、鳳港減河等平原河道。溫榆河主要支流包括藺溝河、清河、壩河、小中河等。北運(yùn)河主要支流包括通惠河、涼水河、鳳港減河等，其中通惠河主要負(fù)擔(dān)北京市城區(qū)、西郊地區(qū)及東郊部分地區(qū)的排水任務(wù)。北運(yùn)河(含溫榆河)作為北京市最重要的排水河道，承擔(dān)著北京城區(qū)90%的排洪任務(wù)，此外，還負(fù)責(zé)天津部分地區(qū)的排洪排污任務(wù)，其水質(zhì)的好壞直接關(guān)系到這些地區(qū)社會(huì)經(jīng)濟(jì)的發(fā)展，對(duì)沿岸地區(qū)實(shí)施可持續(xù)發(fā)展的重要性不言而喻。由于城市工業(yè)的發(fā)展及人口的增加，與日俱增的污水排入溫榆河、北運(yùn)河，其干支流水體已受到嚴(yán)重污染[10-13]，這已經(jīng)成為制約京津地區(qū)社會(huì)發(fā)展的重要因素之一。因此，對(duì)于北運(yùn)河水環(huán)境容量的研究可為該河流水質(zhì)污染治理提供科學(xué)的依據(jù)，對(duì)于京津地區(qū)的經(jīng)濟(jì)、社會(huì)和生態(tài)環(huán)境可持續(xù)發(fā)展也有著重要的意義。

本文以北運(yùn)河為例，基于控制斷面達(dá)標(biāo)的原則，采用水動(dòng)力及水質(zhì)模型建立該河段各支流及點(diǎn)源與控制斷面之間的響應(yīng)關(guān)系，在此基礎(chǔ)上，考慮各點(diǎn)源及支流排污或水質(zhì)現(xiàn)狀，建立各支流及點(diǎn)源排放削減的線性規(guī)劃優(yōu)化模型，對(duì)對(duì)該河段水環(huán)境容量進(jìn)行優(yōu)化分配，為北運(yùn)河水環(huán)境規(guī)劃治理決策的制定提供科學(xué)的依據(jù)。

2 研究方法

2.1 水動(dòng)力及水質(zhì)模型

河流水動(dòng)力及水質(zhì)模型的控制方程為Saint-Venant方程[14-15]和一維對(duì)流擴(kuò)散方程[16-17]，如下所示：

(1)

(2)

(3)

方程(1)～(2)描述了河道中水體的質(zhì)量守恒和動(dòng)量守恒，方程(3)表示水體中污染物量的守恒，其右邊第一項(xiàng)為擴(kuò)散項(xiàng)，反映了污染物在水中隨著水分子的擴(kuò)散而發(fā)生的變化；右邊第二項(xiàng)為對(duì)流項(xiàng)，反映了水中污染物隨著水體運(yùn)動(dòng)發(fā)生的濃度變化；右邊第三項(xiàng)為污染物的降解項(xiàng)，主要反映污染物本身隨時(shí)間變化的降解等作用。通過綜合考慮水中污染物在以上各種作用下的變化，得出水體中污染物濃度的變化規(guī)律。

對(duì)于上述方程的求解有很多方法，主要包括有限差分法、有限體積法等。本文采用有限差分法對(duì)上述方程進(jìn)行求解，其中對(duì)水動(dòng)力方程的求解采用Preissmenn格式[15]，而對(duì)對(duì)流擴(kuò)散方程的求解采用二階隱式差分格式[16]。

2.2 排放源與控制斷面的響應(yīng)矩陣

響應(yīng)矩陣反映了在給定水動(dòng)力條件下各控制斷面對(duì)各支流或點(diǎn)源排污的響應(yīng)關(guān)系，如對(duì)于一個(gè)M×N的矩陣A，其中的元素Aij表示控制斷面i對(duì)點(diǎn)源j的響應(yīng)系數(shù)。控制斷面濃度、響應(yīng)矩陣和點(diǎn)源排放濃度的關(guān)系如下：

Ck=ACp

(4)

式中：Ck為控制斷面的濃度，其為M×1 的矩陣;Cp為支流或點(diǎn)源的排放污染物濃度，為N×1的矩陣;M、N分別表示河流控制斷面和支流或點(diǎn)源的個(gè)數(shù)。

2.3 線性規(guī)劃(LP)優(yōu)化模型

為了對(duì)河段水環(huán)境容量進(jìn)行優(yōu)化，得到各支流及點(diǎn)源的排放量分配方案，需建立河段水環(huán)境容量分配優(yōu)化模型，本文采用線性規(guī)劃模型對(duì)河段水環(huán)境容量的分配進(jìn)行優(yōu)化，如下。

目標(biāo)函數(shù)(總納污量最大原則)：

(5)

約束條件：

(6)

Cpj≥0j=1,2，…,N

(7)

式中：TC為排入干流的污染物的總量，mg；Qi為點(diǎn)源i的流量值，m3/s;Cpj為點(diǎn)源j的排放污染物濃度值，mg/L。Aij為響應(yīng)矩陣;Csi為控制斷面i的水質(zhì)目標(biāo)濃度值，mg/L。

上述線性規(guī)劃模型中目標(biāo)函數(shù)，即方程(5)為所有點(diǎn)源和支流排放污染物量總量最大原則下的優(yōu)化目標(biāo)函數(shù)，在實(shí)際的優(yōu)化過程中也可根據(jù)實(shí)際情況建立此函數(shù)，如總經(jīng)濟(jì)效益最大原則等。此外，約束條件(6)、(7)分別表示每個(gè)控制斷面都必須達(dá)到水質(zhì)標(biāo)準(zhǔn)和各點(diǎn)源的排放濃度不能為負(fù)數(shù)，和目標(biāo)函數(shù)一樣，若實(shí)際問題中有其他對(duì)點(diǎn)源或控制斷面的約束條件，均可添加。

3 計(jì)算條件

采用本文所建模型研究北運(yùn)河環(huán)境容量及分配，需給定相關(guān)水動(dòng)力條件、點(diǎn)源及支流水質(zhì)現(xiàn)狀、河段地表水功能區(qū)劃、控制斷面選取等資料，本節(jié)將分別介紹。

3.1 計(jì)算區(qū)域

研究區(qū)域?yàn)楸边\(yùn)河沙河閘至老米店閘之間河段，河道總長(zhǎng)151 km，如圖1所示。本文研究中，共考慮概化點(diǎn)源或支流共12個(gè)。其中包括沙河閘上游河道(P1)和6條一級(jí)支流，從上游到下游依次為清河(P4)、壩河(P7)、小中河(P9)、通惠河(P10)、涼水河(P11)和龍鳳河(P12)。此外，經(jīng)調(diào)查沙河閘上游河段共有5個(gè)點(diǎn)源，其中沙河閘至清河之間有2個(gè)點(diǎn)源(P2和P3)，清河至壩河之間有2個(gè)點(diǎn)源(P5和P6)，壩河到小中河之間有1個(gè)點(diǎn)源(P8)。

圖1 北運(yùn)河干、支流及控制斷面分布圖

3.2 水環(huán)境功能區(qū)劃及控制斷面的選取

水功能區(qū)劃是指根據(jù)流域區(qū)域的水資源狀況，并考慮水資源開發(fā)利用現(xiàn)狀和經(jīng)濟(jì)社會(huì)發(fā)展對(duì)水量和水質(zhì)的需求，在響應(yīng)的水域劃定的具有特定功能，有利于水資源的合理開發(fā)利用和保護(hù)，能夠發(fā)揮最佳效益的區(qū)域。根據(jù)北京市和天津市的水功能區(qū)劃，北運(yùn)河干流段的水質(zhì)目標(biāo)如表1所示。

本文以化學(xué)需氧量(COD)和氨氮(NH3—N)兩個(gè)代表性污染物指標(biāo)為例，研究其水環(huán)境容量分配，表2為國(guó)家地表水環(huán)境質(zhì)量標(biāo)準(zhǔn)規(guī)定兩個(gè)指標(biāo)的濃度值。

對(duì)于控制斷面的劃分，本文中綜合考慮了水環(huán)境功能區(qū)劃和行政區(qū)劃，設(shè)置了3個(gè)控制斷面，見圖1，控制斷面位置及水質(zhì)標(biāo)準(zhǔn)見表3。其中斷面1為不同水功能區(qū)交界斷面，斷面2為北京市和天津市的分界斷面，斷面3為研究區(qū)域邊界斷面。

表1 北運(yùn)河干流水功能區(qū)段劃分

表2 地表水環(huán)境質(zhì)量標(biāo)準(zhǔn)[18] mg/L

表3 控制斷面位置及水質(zhì)標(biāo)準(zhǔn)

3.3 水動(dòng)力條件和污染物降解系數(shù)的確定

本文主要計(jì)算年平均情況下(設(shè)計(jì)條件下)北運(yùn)河各支流及點(diǎn)源排入干流的污染物水環(huán)境容量?jī)?yōu)化分配，其中水動(dòng)力條件包括計(jì)算邊界條件及支流或點(diǎn)源的入流量等，本文利用流域產(chǎn)匯污模型，根據(jù)北運(yùn)河流域典型年(2008年)的降水等水文資料，計(jì)算出了北運(yùn)河流域各支流及點(diǎn)源的入河流量及污染物的量。

污染物降解系數(shù)是河流水質(zhì)模型中的一個(gè)關(guān)鍵參數(shù)，該參數(shù)取值直接影響到水質(zhì)模型的準(zhǔn)確性。該系數(shù)的確定采用實(shí)測(cè)值，實(shí)驗(yàn)測(cè)得北運(yùn)河不同河段四個(gè)季度的COD和NH3—N降解系數(shù)，其中COD降解系數(shù)取值范圍為0.1145～0.2076 d-1，NH3—N降解系數(shù)取值范圍為0.0369～0.2232 d-1。

4 計(jì)算結(jié)果及分析

采用本文所建立的水動(dòng)力及水質(zhì)模型，計(jì)算得COD和NH3—N的響應(yīng)系數(shù)，如圖2～3所示。其中，第1個(gè)點(diǎn)源為上游邊界，即沙河閘。從圖中可以發(fā)現(xiàn)，響應(yīng)矩陣中控制斷面1對(duì)點(diǎn)源4～12的響應(yīng)系數(shù)為零，控制斷面2對(duì)點(diǎn)源11、12的響應(yīng)系數(shù)為零，表明這些點(diǎn)源的排放污染物量對(duì)對(duì)應(yīng)的控制斷面濃度沒有影響，其原因是點(diǎn)源位于控制斷面的下游，而北運(yùn)河屬于非感潮河流，其水流不會(huì)從下游流向上游。其次，通過對(duì)比發(fā)現(xiàn)，相同水動(dòng)力條件下，兩種污染物濃度的響應(yīng)矩陣差別不大，由于污染物降解時(shí)間較長(zhǎng)，因此，控制斷面對(duì)點(diǎn)源或支流的水質(zhì)響應(yīng)主要由水動(dòng)力過程決定。

圖2 COD響應(yīng)系數(shù)圖3 NH3—N響應(yīng)系數(shù)

根據(jù)前期監(jiān)測(cè)及水文模型結(jié)果，現(xiàn)狀情況下北運(yùn)河3個(gè)控制斷面的濃度均不能滿足水功能區(qū)劃標(biāo)準(zhǔn)要求，因此，本文采用所建立線性規(guī)劃模型對(duì)點(diǎn)源、支流排污削減量進(jìn)行了分配，結(jié)合現(xiàn)狀排污情況，設(shè)置兩種排污優(yōu)化方案，如表4所示。其中，優(yōu)化目標(biāo)函數(shù)為所有支流、點(diǎn)源排放量總和最大。3個(gè)約束條件，分別為各控制斷面滿足水功能區(qū)劃要求的水質(zhì)濃度；各點(diǎn)源、支流最大排放量不超過現(xiàn)狀；各支流、點(diǎn)源的最大削減量。兩種方案的不同之處在于其對(duì)點(diǎn)源或支流的最大削減量要求不同，方案一最大削減量為現(xiàn)狀的50%，方案二最大削減量為現(xiàn)狀的70%。

圖4～5分別給出了兩種方案下12個(gè)支流及點(diǎn)源COD和NH3—N的年分配容量。從圖4～5中可以看出，方案1情況下，各支流及點(diǎn)源的COD的分配容量均為負(fù)值，表明均需要削減。而方案2情況下，支流3和支流8的分配容量為正值，表明還有一定的容量，而其他支流及點(diǎn)源仍需削減。兩種方案下，點(diǎn)源12的COD年削減量較大，均超過15 000 t。結(jié)合調(diào)查統(tǒng)計(jì)COD年排放量，可進(jìn)一步計(jì)算出各斷面的COD削減率，兩種方案情況下，北運(yùn)河各支流的COD平均削減率在30%～50%之間。兩種方案情況下，NH3—N的容量并沒有很明顯的變化，均為負(fù)值，需要進(jìn)一步削減。這表明點(diǎn)源最大削減量的變化對(duì)其優(yōu)化分配沒有產(chǎn)生明顯的影響。其中，點(diǎn)源11、12的年削減量最大，均超過3 000 t。結(jié)合現(xiàn)調(diào)查統(tǒng)計(jì)排放量，兩種方案情況下，各點(diǎn)源或支流NH3—N平均削減率在40%～70%之間。若要達(dá)到水環(huán)境功能區(qū)劃要求，仍需進(jìn)一步提高各支流及點(diǎn)源的最大削減量比例。

表4 兩種排污優(yōu)化方案

圖4兩種方案下各點(diǎn)源或支流的COD年分配量圖5兩種方案下各點(diǎn)源或支流的NH3—N年分配量

5 結(jié) 論

水環(huán)境容量分配的合理性對(duì)實(shí)現(xiàn)水系水環(huán)境規(guī)劃治理有著非常重要的作用，本文通過水動(dòng)力水質(zhì)模型計(jì)算得到了點(diǎn)源與控制斷面水質(zhì)響應(yīng)矩陣，在此基礎(chǔ)上，建立了水環(huán)境容量分配的線性規(guī)劃優(yōu)化模型。采用所建立的模型對(duì)北運(yùn)河干流段水環(huán)境容量分配進(jìn)行了研究，得出以下結(jié)論：

(1)除個(gè)別支流或點(diǎn)源外(點(diǎn)源3和8)，兩種方案下北運(yùn)河各支流及點(diǎn)源的COD和NH3—N分配容量均為負(fù)值。與監(jiān)測(cè)統(tǒng)計(jì)排放量比較，各點(diǎn)源或支流COD平均削減率為30%～50%，NH3—N平均削減率在40%～70%之間。

(2)線性規(guī)劃模型作為一種有效的優(yōu)化方法可以用于水環(huán)境容量的分配研究中，通過改變不同的優(yōu)化約束條件，實(shí)現(xiàn)對(duì)不同影響因素影響的考慮，最終達(dá)到對(duì)水環(huán)境容量的優(yōu)化分配。

(3)在水環(huán)境容量?jī)?yōu)化分配中，可能會(huì)無法找到最優(yōu)結(jié)果，這時(shí)就需要對(duì)約束條件進(jìn)行修改，比如，考慮到經(jīng)濟(jì)效益等影響，一些支流或點(diǎn)源無法實(shí)現(xiàn)削減，會(huì)導(dǎo)致優(yōu)化模型找不到最優(yōu)解。