基于貝葉斯公式的不確定性水環境容量研究
——以北運河為例

曹若馨,曾維華,李 晴,馬鵬宇 (北京師范大學環境學院,北京 100875)

水環境容量的概念最早是指在滿足水體水功能目標的前提下,水體可容納污染物的最大負荷量.歐美學者對環境容量這一稱謂較少使用,而一般將其稱為同化容量[1]或納污能力、水體允許排污水平[2-3]等.

通常,水環境容量核算是基于確定設計條件的,即最枯周(7Q10)或最枯月(30Q10)設計流量條件下的最大容許排放量[4-5].基于不同的實際情況建立水質模型,并對水環境容量進行計算[6-7].周孝德等[8]提出了一維穩態設計條件下計算河流水環境容量的三種方法,即段首控制法,段尾控制法和功能區段尾控制法,并對各種計算方法適用條件進行了分析討論.美國等發達國家在確定條件下的允許排污負荷方面的研究起步較早,經過近五十年的發展,美國已從最初的最枯月、周最大允許排污負荷,發展到目前的最大日容許排污負荷(TMDL),兼顧面源污染與安全裕量[9].國外很多學者基于系統優化原理將水體允許排污負荷計算與污染負荷分配統一考慮,即水質規劃[10-11].

但是,由于決定水環境容量研究涉及到的流量與背景濃度等因素的隨機不確定性,導致水環境容量同樣具有隨機不確定性特征.確定條件下的水環境容量核算的不足之處逐漸顯露出來,由此引起一些學者對基于設計流量的確定性水環境容量核算方法產生了質疑.因此,從20 世紀80 年代末開始,就有學者著手研究不確定水環境容量.曾維華等[12]從水環境的隨機不確定性出發,提出了風險水環境容量的隨機估算模式.除了運用隨機理論外,之后還有學者從水系統的灰色性[13]、模糊性[14]角度出發,計算水體納污能力.相應的研究方法主要有蒙特卡洛法[15]、區間數理論[16]、盲數理論[17-18]等.

常用的傳統確定性水環境容量核算方法與TMDL 技術體系中對污染負荷容量的核算方法,都屬于確定性容量核算方法體系,得到的都是確定的數值,不符合水系統的隨機特性.而隨機不確定性水環境容量核算方法,如隨機估算模式等,往往需要依靠大量監測數據,這就導致研究過程中通常因為難以滿足數據量和數據質量而影響結果的準確性甚至難以進行實際應用.考慮到現有方法的不足之處,本文引進概率論中的貝葉斯公式,將歷史經驗,包括傳統核算知識,作為先驗知識,通過利用實際監測數據對先驗知識進行校正,最終得到水環境容量的后驗分布.這一方面能夠通過使用驗前知識,對實際監測樣本進行合理補充,在一定程度上保證結果的可靠性;另一方面能夠得到容量概率分布而非確定性數值,從而對水環境容量有一個更為全面準確的描述.

本文以北運河流域為研究對象,以月為研究尺度,基于流域自然匯水邊界將研究區劃分為一系列控制單元,同時為了進行比較分析,并且保證不確定性水環境容量核算的合理性與準確性,首先對研究流域進行傳統確定性COD 水環境容量核算,然后通過分析各個參數作為隨機變量的概率分布,利用貝葉斯方法對其不確定性COD 水環境容量進行核算,最終得到各計算單元的水環境容量分布,并與傳統方法核算的結果進行比較.

1 研究方法

1.1 確定性水環境容量核算方法

1.1.1 單元劃分與斷面設定 計算單元的劃分通常以保證重要水域水體功能和保持單元內穩定計算條件為原則,以城市及工業區、生活區等重要敏感的區域或斷面為劃分點進行節點劃分.一維穩態計算條件下,水環境容量研究通常采用段末濃度控制[8]的方法:污染物濃度在研究河段的段末達到水質目標濃度,此時全河段污染物濃度高于該段水質目標濃度.

1.1.2 排污口調查與概化 為了更便捷的進行研究,通常將排污口在河段內的分布進行概化.重心概化法是一種常用的排污口概化方法,根據排污口的排放量與排放位置,將所有排污口概化到其重心位置,即認為河段內污染源排污口在河段排污口重心斷面排放[19].重心位置根據式(1)計算,排污量根據式(2)計算.

式中:Xi為第i 河段概化排污口到河段段末的距離,m;Mni為第i 河段第n 個排污口的排污量,g/s;Xni為第i 河段第n 個排污口到河段段末的距離,m.

1.1.3 水環境容量水質模型選取 水質模型的選取參考《水域納污能力計算規程》[20]中的相關規定.一維水質模型適用于寬深比不大的中小河流,污染物質基本在河段內均勻混合,斷面污染物濃度橫向變化不大,主要適用于Q＜150m3的中小型河流.

1.1.4 模型參數設定 確定性水環境容量核算涉及的參數包括:流量、流速、水質目標以及污染物綜合衰減系數等.

(1)設計流量與流速.設計流量的設定一般采用近10 年最枯月平均流量或90%保證率最枯月平均流量,集中式飲用水源地采用95%保證率最枯月平均流量.對于北方地區部分河流,可根據實際情況適當調整設計保證率,也可選取偏枯典型年的枯水期流量作為設計流量.設計流速是與設計流量是相對應的,當實測資料不夠充足時,通常可通過經驗公式,基于實測數據推導得出.式(3)是常采用的流速-流量經驗關系式[21]:

式中:U 為流速,m/s;Q 為流量,m3/s;a、b 是水力學參數,即經驗參數,由實測資料確定.

(2)水質目標.各單元水質目標濃度可根據《國家地表水環境質量標準(GB3838-2002)》[22]確定.水功能區水質目標與水質目標濃度值對應關系如表1 所示.進行確定性水環境容量核算時,水質目標濃度也作為下游斷面的背景濃度進行計算.

(3)污染物綜合衰減系數.污染物綜合衰減系數的設定方法采用實測資料反推法[20]:選取順直、水流穩定、無支流匯入、無入河排污口的河段,在上下游分別對污染物濃度和流速進行監測,選擇合適的水質模型計算污染物降解系數.一維模型中可通過以下公式計算K 值:

式中:u 為河流平均流速,m/s;x 為上下游斷面間距離,m;CA為上游斷面污染物濃度,mg/L;CB為下游斷面污染物濃度,mg/L.

表1 水質目標與COD 目標值對應關系Table 1 Water quality goals and COD concentration goal values

1.2 基于貝葉斯公式的不確定性水環境容量核算方法

基于貝葉斯公式的不確定性水環境容量核算與傳統確定性水環境容量核算方法的基本步驟一致,然而不同之處在于,基于貝葉斯公式的水環境容量核算的參數都是服從某一概率分布的隨機變量,而非某一確定的數值,因此需要確定各個參數相對應的分布.

1.2.1 流量分布 根據歷史經驗,大部分水文參數服從皮爾遜III 型(P-III)分布.本文根據各水文站歷年逐月的流量數據,采用常用的矩法對P-III 型曲線進行參數估計,并基于優化適線法比較理論頻率與經驗頻率,選取最優的參數,從而最終得到流量的概率分布.

確定P-Ⅲ分布的參數包括形狀參數α、尺度參數β和位置參數α0.這三個參數可由樣本均值Ex、變差系數Cv和偏態系數Cs推導得出[23-24].其計算公式為:

式中:xi為序列中的值;n 為序列長度;σ為序列標準差;Sc為偏斜度.

兩組參數間的關系為:

為了檢驗得到的理論頻率曲線是否較為合理,需要將其與經驗頻率進行比較.在對水文變量序列由大到小排序后,得到的新序列中每一項對應的經驗頻率p 為:

式中:m 為序數; p 為經驗頻率.

基于優化適線法,采用離差平方和最小準則(OLS),使經驗點和同頻率的頻率曲線縱坐標之差的平方和達到最小.在選擇候選參數時,通常設置Cs為Cv的某一倍數.對于P-III 型曲線,就是使下列目標函數式取最小[25]:

即

式中:Q 為參數(Ex,Cv;Cs), Q’為參數Q 的估計值;Pi為頻率;n 為系列長度;f(Pi,Q)為頻率曲線縱坐標.

1.2.2 流速分布 與傳統確定性水環境容量核算一致,流速與流量是對應的關系,通過流速-流量經驗關系式推導得出.得到經驗參數后,流速分布也就是隨機變量函數的分布,函數關系式為式(3).

1.2.3 背景濃度分布及水質目標 背景濃度分布根據各河段段首實際監測水質數據來進行分析,可利用編程軟件進行數據觀察及數據擬合分析,從而確定其分布.每個河段的水環境質量目標取各河段背景水質濃度的數學期望值.

1.2.4 污染物綜合衰減系數分布 基于貝葉斯公式[26]對污染污染物綜合衰減系數k 進行率定:

式中:k 為污染污染物綜合衰減系數,是一隨機變量;p(k|D)為參數k 的后驗分布,也就是參數k 的率定結果分布,表示在獲得觀測數據之后模型參數的分布規律;D 為樣本信息;p(k)為參數k 的先驗分布函數,表示在未收集數據前關于參數的認識;p(D|k)為似然函數,表示模型參數擬合實測數據的程度,值越大擬合程度越好.由于p(D)是常數項,因此,參數的后驗概率正比于似然函數與先驗概率的乘積:

同樣選取順直、水流穩定、無支流匯入、無入河排污口的河段,利用實際污染物濃度數據(樣本)D和所選擇的水質模型,基于平方誤差損失函數計算對比段末水質的計算值L′2與實際值L2,計算得到觀測樣本對于參數k 的似然函數p(D|k):

這里暫時無須要求似然函數滿足歸一性,關注的是參數變化時函數的變化而非其具體取值.最后得到似然函數與先驗概率的乘積(式9),進行歸一化后得到參數k 的后驗概率p(k|D).

1.2.5 水環境容量后驗分布 基于貝葉斯公式對水環境容量W 進行核算:

為保證結果的可靠性以及對最后的合理性驗證,先驗分布p(W)根據不同傳統設計條件,也就是不同保證率下的一系列確定性水環境容量核算結果值來進行設置.

基于已得到的各參數的分布,對各參數進行相應的隨機抽樣.最后基于平方誤差損失函數計算對比段末水質的計算值C′與水質目標濃度Cs,從而得到觀測樣本對于W 的似然函數p(D|W):

同樣這里無須要求似然函數滿足歸一性.最后得到似然函數與先驗概率的乘積(式11),進行歸一化后得到W 的后驗分布p(W|D).

2 北運河水環境容量核算案例研究

2.1 研究區概況

北運河發源于北京市昌平區燕山南麓,西界永定河,東臨潮白河,干流從西北方向依次經過北京市順義區、朝陽區、通州區、廊坊市香河縣、天津市武清區,最后在天津市大紅橋匯入海河.在北京境內,北運河干流總長90km,境內流域面積為4348km2,是北京市內平原流域最廣的水系,北京有70%以上的人口在該流域內生活工作.北運河沿線排污口眾多,污水排放量大,其作為北京市最重要的排水河道,承擔著北京城區90%的排洪任務,此外,北運河還承擔天津部分地區的排污排洪任務,可見北運河的水環境狀況與這些地區的社會經濟發展息息相關.

圖1 北運河流域河網分布Fig.1 Beiyunhe river network

圖2 北運河流域河段概化Fig.2 Generalized network of Beiyunhe river basin

2.2 數據來源

本研究數據資料主要來自2008～2017 年的《北京區域統計年鑒》、《天津統計年鑒》、《廊坊經濟統計年鑒》、《河北農村統計年鑒》等統計數據,北京、天津、廊坊市工業企業排污等統計核算數據,以及水文站監測站數據與水環境監測數據等.

2.3 研究流域單元劃分

表2 計算單元劃分Table 2 Division of calculation units

北運河流域計算單元劃分采用控制單元區劃點為主,水力條件急劇變化點為輔的劃分方式.單元劃分結果如表2 所示.

2.4 水環境容量計算模型

北運河流域降水偏少,流量遠小于150m3/s.根據《水域納污能力計算規程》(SL348-2006)[20],北運河適用于一維穩態水質模型.采用國內通常使用的現狀-段末濃度控制的方法計算研究河段的水環境容量.在考慮排污口流量時,河段納污能力公式為:

式中:W 為河段水環境容量,g/s;Q 為河段流量,m3/s;Qi為河段中i排污口的排口流量,m3/s;C0為河段初始濃度,mg/L;Cs為河段水質目標,mg/L;k 為綜合降解系數,1/d;u 為河段平均流速,m/s;l1為河段概化的排污口到計算河段段首處的距離,m;l2為河段概化的排污口到計算河段段末處的距離,m.

2.5 排污口概化

根據統計的污染排放數據采用重心概化法對北運河流域各計算單元沿岸的排污口進行概化,對于干流河段,其支流作為穩定排放源一同進行概化.概化結果如表2.

由于本研究選取段末濃度控制現狀水環境容量核算模型,因此只需考慮存在排污源的10 個計算單元.

2.6 確定性水環境容量核算

2.6.1 參數設定 根據北運河流域各水文站流量,以及排污源數據,得到各計算單元的設計流量如表3所示.

表3 北運河流域各計算單元設計流量Table 3 Design water flow of each calculation unit in Beiyunhe river basin

表4 各計算單元水力系數及水質目標Table 4 Hydraulic coefficients and water quality goals of each calculation unit

基于式(3),水力學參數a、b 根據北運河各水文站流量-流速曲線推導,鄰近河段根據就近原則,結合自身水力條件選取就近水文站水力參數值,最終確定各河段水力學參數.各計算單元水力系數及水質目標如表4 所示.

在研究北運河流域確定性水環境容量時,綜合衰減系數采用實測資料反推法得出.根據北運河流域監測斷面的分布,選取順直、水流穩定、無支流匯入、無入河排污口的河段榆林莊-王家擺,基于式(4)進行推導,經推算,污染物綜合衰減系數K= 0.255d-1.2.6.2 確定性水環境容量結果 根據以上設定參數,各河段在不同保證率下的確定性水環境容量如表5 所示.

表5 北運河流域COD 確定性水環境容量Table 5 COD water environmental capacity in Beiyunhe river basin

2.7 基于貝葉斯公式的不確定性水環境容量核算

為了高效的處理大量數據,并進行可視化分析,本文利用R 語言編程,得到各參數的分布.

2.7.1 流量與流速 由于流速與流量是對應的關系,其關系式與經驗參數與傳統確定性水環境容量核算一致,因此只需要得到流量的分布.

為了得到P-III 曲線對應的形狀參數α、尺度參數β和位置參數α0,首先需要得到樣本均值Ex、變差系數Cv和偏態系數Cs.在R語言中均值Ex可由mean函數計算得到,標準差σ可由sd 函數計算得到,偏斜度Sc可由skewness 函數計算得到.在計算出α、β、α0后,可通過調用pgamma 函數與qgamma 函數求理論頻率值P 與水文變量值xp.

圖3 各計算單元流量概率密度分布Fig.3 Probability density distribution of water flow in each calculation unit

圖4 各計算單元背景濃度概率密度分布Fig.4 Probability density distribution of background pollutant concentration in each calculation unit

通過R 語言編程得到經驗頻率后,利用基于離差平方和的優化適線法,適當調整偏態系數Cs,比較不同參數的擬合效果,最后選取最優的一組參數.得到各河段的流量分布如圖3 所示.2.7.2 背景濃度 為得到未知的背景濃度分布,針對各首斷面連續實測逐月監測數據,首先要對數據進行觀察分析.利用plotdist 函數觀察數據的直方圖或密度圖,另外利用descdist 函數得到對數據的一系列描述性統計,包括最大最小值、均值、方差、偏度、峰度以及偏度-峰度圖等.通過對數據的觀察分析,選取幾個較為合理的候選分布,通過fitdist 函數分別進行擬合分析.使用fitdist 函數中提供的最大似然估計法(mle)進行擬合分析,分別得到各候選函數對應擬合結果的參數估計值、標準差、似然度等.最后通過denscomp 函數和cdfcomp 函數將各候選函數的擬合結果曲線圖與原始數據的密度分布及累積概率分布分別進行比較,擬合結果最好的則為最終選定的理論分布[27].大部分水質分布都服從gamma 分布.

2.7.3 污染物綜合衰減系數 根據北運河實際情況及歷史文獻資料[28-29],本文設置COD 降解系數K的先驗分布為0.01～0.8d-1上的均勻分布p(k).選取順直、水流穩定、無支流匯入、無入河排污口的河段榆林莊-王家擺,利用實際污染物濃度數據(樣本)D 和所選擇的水質模型,基于平方誤差損失函數計算對比段末水質的計算值L2′與實際值L2,計算得到觀測樣本對于參數k 的似然函數p(D|k).最后得到似然函數與先驗概率的乘積,從而得到參數k的后驗概率,如圖5 所示.k 值在0.255d-1處概率密度最大.

2.7.4 水環境容量后驗概率分布 根據各河段在不同保證率下的確定性水環境容量值,相對應的設置每個計算單元的水環境容量的正態先驗分布p(W).利用已得到的各個參數的分布,對每個參數進行隨機抽樣.

對于常見分布如gamma 分布,可直接利用R 語言中的抽樣函數rgamma 對其進行抽樣.對于沒有抽樣函數的分布,如參數k 的分布,可利用舍選抽樣法(Acceptance-rejection Sampling)進行抽樣.舍選抽樣法的抽樣規則為:隨機變量k 的取值區間為[a,b],其概率密度為f(k),且max{f(k)|a≤k≤b}=c,首先產生[a,b]區間內均勻分布的隨機數x,然后產生[0,c]區間內均勻分布的隨機數y,當y≤f(x)時,接受x 為所需要的隨機數,否則返回到第一步重新抽取一對(x,y).

圖5 污染物綜合衰減系數概率密度分布Fig.5 Probability density distribution of attenuation coefficient of pollutants

基于平方誤差損失函數計算對比段末水質的計算值C′與水質目標濃度Cs,從而得到觀測樣本對于參數W 的似然函數p(D|W).最后得到似然函數與先驗概率的乘積,從而得到W 的后驗分布p(W|D).

3 結果與討論

3.1 水環境容量核算結果

得到有排口的10 個計算單元的水環境容量分布如圖6 所示,具體參數如表6.

圖6 北運河流域COD 水環境容量概率密度分布Fig.6 Probability density distribution of COD water environmental capacity in Beiyunhe river basin

圖7 北運河流域COD 水環境容量累積概率分布Fig.7 Cumulative distribution of COD water environmental capacity in Beiyunhe river basin

表6 不確定性COD 水環境容量核算結果參數Table 6 Result values of uncertain COD water environmental capacity

為了更直觀的看到觀測數據對先驗信息進行的“校正”,將水環境容量概率密度分布轉化為累積概率分布,同時將確定性水環境容量值(表5)在累積概率分布圖上用圓點標記出來(圖7).由于流量保證率體現的是不低于該設計流量的月均流量百分比,因此在累積概率分布上所對應的概率應該是1-保證率,90%保證率對應的累積概率應該是1-90%= 10%.

3.2 結果分析

從各河段的水環境容量分布來看,除了筐兒港-新老米店閘河段為左偏態,其余河段均為右偏態,相比較于峰度為3 的標準正態分布,各個分布的峰度都小于3,都屬于瘦尾分布.

北運河干流河段有污染物濃度更高、流量更大、水環境容量更大的特征,干流的5 個河段分別為土溝橋-溫榆河順義區、溫榆河順義區-榆林莊、王家擺-土門樓、秦營揚水站-筐兒港、筐兒港-新老米店閘,水環境容量期望值大多在1萬t/a以上,包括全流域水環境容量期望值最大的溫榆河順義區-榆林莊河段(4.56 萬t/a),其中期望值最小的是筐兒港-新老米店閘河段(6881.5t/a).支流河段相對而言水環境容量較小,除了清河閘-沙子營河段流量較大,大紅門閘上-涼水河口河段流量大且河段長度長,導致這兩段的水環境容量較大,分別為1.08 萬t/a 和2.15 萬t/a,其余河段容量都較小,包括全流域水環境容量期望值最小的南沙河入昌平-南沙河口河段(782.1t/a).

對于河段長度較長、河段流量較大的河段,水環境容量往往較大,如水環境容量期望值最大的溫榆河順義區-榆林莊河段,該河段最長,長度為45.10km,匯入支流較多;相反對于流量小、長度短的河段,其水環境容量較小,如南沙河入昌平-南沙河口河段的水環境容量期望值最小,長度10.23km,流量小,集中在3m3/s 以下,相較而言,王家擺-土門樓河段長度雖然較短,為 8.92km,但流量較大,集中在 20m3/s～30m3/s,其水環境容量也較大.背景濃度對水環境容量的分布也存在影響,但各河段背景濃度的分布差別相對較小.

比較傳統設計條件下的確定性水環境容量值,從圖7 的累積概率分布圖上可以看出,不同保證率下的確定性水環境容量基本符合不確定性水環境容量的概率分布趨勢,設計條件下的確定性水環境容量體現的是先驗信息,貝葉斯方法基于實際監測數據對先驗信息進行了“修正”,一方面,得到的后驗累積分布的基本形狀和趨勢仍然比較符合先驗信息,可見結果的合理性,另一方面,這也進一步顯示流量對水環境容量核算的顯著影響.

陳丁江等[30]分析指出設計流量是對水環境容量最靈敏的輸入參數,其次是綜合衰減系數.從水環境容量的影響因素上看,根據研究結果,研究區的水環境容量主要受河段流量與長度的影響,其次是背景濃度與水質目標濃度.

4 結論

4.1 各河段水環境容量呈瘦尾偏態分布,其中筐兒港-新老米店閘河段為左偏態,其余河段均為右偏態.

4.2 水環境容量的分布受河段流量和河段長度影響最大,整個流域內干流河段水環境容量較大,支流相對較小,全流域水環境容量期望值最大的是溫榆河順義區-榆林莊河段,為4.56 萬t/a,該河段長度最長,而水環境容量期望值最小的是南沙河入昌平-南沙河口河段,為782.1t/a.

4.3 傳統設計條件下的確定性水環境容量核算結果與基于貝葉斯公式的不確定性水環境容量核算結果基本趨勢一致,后者不僅體現了前者所包含的信息,并且針對各參數的隨機特性,利用實際監測數據對前者進行了“修正”.