黎 坤，林凱榮，江 濤，陳曉宏

中山大學(xué)水資源與環(huán)境研究中心，廣東 廣州 510275

數(shù)字濾波法在點源和非點源污染負荷分割中的應(yīng)用

黎 坤，林凱榮，江 濤，陳曉宏

中山大學(xué)水資源與環(huán)境研究中心，廣東 廣州 510275

以分辨流域點源和非點源污染負荷為研究目標，在傅立葉分析的基礎(chǔ)上，提出了從污染負荷時間序列中分辨點源和非點源污染負荷的新方法——數(shù)字濾波法，探討了數(shù)字濾波方程的參數(shù)與濾波次數(shù)的關(guān)系，并將該方法應(yīng)用到東江流域CODMn負荷的分割上.應(yīng)用實例表明，分割的點源和非點源負荷系列曲線符合點源和非點源負荷的產(chǎn)生特點，使用者能夠在污染負荷分割過程中通過濾波參數(shù)的選取比較方便地加入自己的經(jīng)驗.一般取較大濾波參數(shù)時，只需要3次濾波就能得出較滿意的結(jié)果.

非點源;污染負荷;數(shù)字濾波;解析

由于非點源具有隨機性、廣泛性和難監(jiān)測性，在我國對非點源的研究還處于起步階段.目前的水質(zhì)監(jiān)測主要是在可控制的河流、湖泊、水庫的進出口斷面上，無法區(qū)分點源和非點源污染負荷.如何在現(xiàn)有資料的基礎(chǔ)上來分辨點源和非點源負荷，為污染源的解析提供依據(jù)，成為目前應(yīng)該解決的一個關(guān)鍵問題.筆者在傅立葉分析的基礎(chǔ)上，提出了從污染負荷時間序列中分辨點源和非點源污染負荷的新方法——數(shù)字濾波法.

數(shù)字濾波法是近年來國際上研究最多的分割基流的方法［1-7］.LYNE［8］于 1979 年首次提出用該方法來分割水文過程中徑流對降雨的快速響應(yīng)和慢速響應(yīng)成分.NATHAN等［9］研究認為，該方法是一種快速而且客觀的基流分割方法，同時發(fā)現(xiàn)濾波參數(shù)在0.900～0.950內(nèi)，能夠得到較好的分割結(jié)果.ARNOLD等［10］把數(shù)字濾波法與其他基流分割方法進行了分析比較，結(jié)果表明，采用3次濾波且濾波參數(shù)取為0.925能夠得到較好的分割結(jié)果;同時認為該方法很容易使用，而且能夠再現(xiàn)分割的結(jié)果;使用者能夠在分割基流的過程中比較方便地加入自己的經(jīng)驗.MAU 等［11］則發(fā)現(xiàn)，濾波參數(shù)取為 0.850，進行多次濾波可以得到最好的結(jié)果.SPONGBERG［12］在傅立葉分析的基礎(chǔ)上，闡述數(shù)字濾波分割基流的基本原理，并在時間域和頻率域上對數(shù)字濾波法分割基流的結(jié)果進行了分析驗證.但鮮見應(yīng)用數(shù)字濾波法來研究非點源的解析，筆者選擇應(yīng)用數(shù)字濾波法分辨點源和非點源污染負荷通量.

1 數(shù)字濾波法基本原理

數(shù)字濾波法源自傅立葉分析，所以在陳述數(shù)字濾波法基本原理時有必要了解傅立葉分析原理.

1.1 傅立葉分析

傅立葉分析的研究與應(yīng)用至今已經(jīng)歷了100余年.傅立葉分析主要是通過傅立葉變換由時域分析轉(zhuǎn)入變換域分析.傅立葉變換是在以時間為自變量的“信號”與以頻率為自變量的“頻譜”函數(shù)之間的某種變換關(guān)系，即建立了時間函數(shù)和頻譜函數(shù)之間的轉(zhuǎn)換關(guān)系.

1.1.1 傅立葉變換

如果x(t)是定義在整個實軸上的實值或復(fù)值函數(shù)，則其傅立葉變換可由下式給出：

若對任意參數(shù)f，上述積分都存在，則式(1)確定了函數(shù)X(f)，稱為X(f)的傅立葉變換.如果已知X(f)，則利用如下的傅立葉逆變換，還可復(fù)原x(t)：

若x(t)和X(f)同時滿足式(1)和(2)，則稱它們是一個傅立葉變換對，記為 x(t)?X(f).通常X(f)是一個復(fù)函數(shù)，因此可以寫成如下兩部分：

式中，R(f)，I(f)分別是 X(f)的實部和虛部.將式(3)表示為指數(shù)形式：

其中，

工程技術(shù)中，常將x(t)看成時間信號，相應(yīng)的空間，稱為時間域和空域;將其傅立葉變換 X(f)看成頻率函數(shù)，相應(yīng)的空間稱為頻域. X(f )稱為x(t)的傅立葉譜，而 φ(f)稱為其相角，這在物理上是有良好背景的.該頻率的含義可以這樣來理解：應(yīng)用歐拉公式可將指數(shù)項表示成正弦－余弦的形式，如果把式(1)解釋成離散項和的極限，則顯然X(f)是包含了無限項正弦－余弦的和，而且 f的每一個值確定了所對應(yīng)的正弦－余弦的頻率.

1.1.2 離散傅立葉變換

實際上，水文時間序列常以離散的形式給出，因此常采用的是離散傅立葉變換(DFT).離散傅立葉變換的一般形式如下：

給定 N 個 實 或 復(fù) 的 數(shù) 列{x(0)，x(1)，…，x(N － 1) } ，定義

X(n)為{x(k)}的離散傅立葉變換.

離散的頻率可表示為：

如果不考慮式(7)中的復(fù)指數(shù)部分的運算，則求解式(7)共需要N×N次乘法和N×(N－1)次加法.顯然當N很大時，其工作量是相當可觀的.因此，通常采用凱萊和卡柯［13］提出的專門用于處理DFT的快速算法(FFT)，可減少DFT的計算時間.

1.2 濾波原理

數(shù)字濾波系統(tǒng)的基本輸入輸出關(guān)系為：

式中，yk為濾波輸出項;ck為濾波系數(shù);xk為輸入項;*表示卷積.

卷積原理表明，在時域上的卷積等價于在頻域上的乘積，故有：式中，T(n)為ck的離散傅立葉變換.

式(10)表明，濾波輸出項(時間序列)的離散傅立葉變換，變換過程如1.1節(jié)所示，可以表示為2個復(fù)數(shù)的頻率序列的乘積，用極數(shù)的形式可表示為：

式中，z為復(fù)數(shù);r為復(fù)數(shù)的模(幅值);θ為相位角.

從式(10)和(11)可以看出，初始序列的濾波輸出項的大小取決于T(n)的幅值，為此，把T(n)稱為濾波系數(shù)ck的頻率轉(zhuǎn)移函數(shù)，簡稱FTF［14］.如果FTF為復(fù)數(shù)，其相位角不為零，那么它所表示的濾波將會改變初始時間序列的相位角.

線性遞歸濾波的一般形式為：

式中，dm為濾波系數(shù);J和M分別為cj和dm個數(shù)的總和.

由式(12)表示濾波的頻率轉(zhuǎn)移函數(shù)(FTF)可表示為［15］：

1.3 頻譜分析

污染負荷是指區(qū)域或某環(huán)境要素對污染物的負載量［16］，目前對污染負荷的監(jiān)測是通過水量水質(zhì)同步監(jiān)測來實現(xiàn)的，即采集水樣的同時進行流量觀測，其監(jiān)測的是污染負荷隨時間的非連續(xù)變化過程，是某斷面的污染負荷的通量，流域某斷面的污染負荷通量是點源負荷和非點源負荷的總和［17］.

點源污染是指有固定排放點的污染源，如工業(yè)廢水及城市生活污水等，由排放口集中匯入水體引起的污染;非點源污染是指溶解性或固體污染物在大面積降水和徑流沖刷作用下匯入受納水體而引起的污染.點源受降水等自然因素的影響較小，其排放過程相對穩(wěn)定.而非點源污染受降水等因素影響，排放呈現(xiàn)隨機性，波動較大.

一般某流域出口斷面點源污染負荷的出流時間序列是相對穩(wěn)定的序列.對點源污染負荷的出流時間序列進行傅立葉變換，其變換后的幅值波形圖見圖1.由圖1可見，點源污染負荷的出流時間序列頻域下的信號主要是低頻信號.

圖1 點源污染負荷的時間序列傅立葉變換后的幅值波形圖Fig.1 Power spectrum of point pollutant series after Fourier transform analysis

非點源污染負荷的出流時間序列由于受降水等因素影響波動很大，對非點源污染負荷的出流時間序列進行傅立葉變換，其變換后的幅值波形圖見圖2.由圖2可見，非點源污染負荷的出流時間序列頻域下的信號明顯比點源污染負荷出流的信號要大得多，且主要是高頻信號，這也就是可以利用分離高低頻信號進行點源和非點源污染負荷系列分割的理論依據(jù).

但是，從圖2還可以看出，非點源污染負荷的時間序列在頻域下的信號并不能完全看作是高頻信號，因為它還包括了部分低頻信號，這就意味著要想完全地分割點源和非點源污染負荷是不可能的，因為有一部分頻率是交疊的.污染負荷系列中低頻率信號的分離會部分地削弱點源污染負荷.但是，濾波的頻率轉(zhuǎn)移函數(shù)(FTF)提供了能夠得到最好分割點源和非點源污染負荷時間序列的途徑.通過頻率轉(zhuǎn)移函數(shù)選擇最佳的濾波系統(tǒng)使得在最小削弱點源負荷的同時最大限度地削減非點源污染負荷，而且盡可能地消除相位失真.

圖2 非點源污染負荷的時間序列傅立葉變換后的幅值波形圖Fig.2 Power spectrum of non-point pollutant series after Fourier transform analysis

1.4 數(shù)字濾波算法

筆者將數(shù)字濾波法引入到流域點源和非點源污染負荷的分割上.參考LYNE［8］于1979年提出的濾波方程，點源和非點源污染負荷分割濾波方程為：

式中，Pt為t時段內(nèi)過濾后的污染負荷出流量(非點源負荷出流量)，kg，計算中 P1設(shè)為0;Bt為 t時段的污染負荷出流總量，kg;β為無量綱參數(shù)，稱為濾波參數(shù).已知污染負荷出流總量和非點源污染負荷出流量后，點源污染負荷出流量(Nt)可由下式求出：

由式(14)和(15)可試算出點源和非點源污染負荷的出流過程，從而對污染負荷的出流過程線進行點源和非點源出流負荷分割.

濾波參數(shù)(β)需要優(yōu)化率定.對于濾波參數(shù)和濾波次數(shù)，不同的學(xué)者有著不同的認識.下面通過頻率轉(zhuǎn)移函數(shù)來分析LYNE［8］提出的數(shù)字濾波法的濾波參數(shù)和濾波次數(shù)的關(guān)系.

把式(14)轉(zhuǎn)化為式(12)的形式，可得：

且有 J=1，M=1.

根據(jù)式(12)，可得 LYNE［8］數(shù)字濾波的頻率轉(zhuǎn)移函數(shù)為：

由式(17)通過Matlab軟件轉(zhuǎn)換，可以得到頻率轉(zhuǎn)移函數(shù)的幅度譜和相位譜，如圖3，4所示.

圖3 不同濾波參數(shù)下LYNE濾波的FTF的幅值波形圖Fig.3 Power spectrum curves of FTF based LYNE filter with differentβ

圖4 不同濾波參數(shù)下LYNE濾波的FTF的相位譜Fig.4 Phase curves of FTF based LYNE filter with differentβ

從理論上來說，非點源污染負荷信號將全部通過濾波，而點源污染負荷信號則不能通過.從圖3可以看出，對所有比較常用的濾波參數(shù)來說，在最高頻率范圍中，F(xiàn)TF的幅值為1.0，包含全部非點源污染負荷信號.而多重濾波并不能進一步衰減高頻率信號，因為在第一次濾波中它們已經(jīng)完全衰減了.增大濾波參數(shù)將會增加點源污染負荷和非點源污染負荷信號中低頻率部分的衰減量.使用較小的參數(shù)將會減少該頻率范圍內(nèi)的衰減量.雖然這在FTF的振幅范圍內(nèi)不是很明顯，但在實際應(yīng)用中可以發(fā)現(xiàn)，如果選用比較小的參數(shù)，以上一次濾波結(jié)果作為下一次濾波輸入，通常要濾波5～6次才能得到合理的點源和非點源負荷的分割.當使用較大濾波參數(shù)時，只需濾波2～3次就可以得到相同的分割量.

圖4顯示的是 β為 0.850，0.925和 0.975的LYNE濾波的FTF的相位譜.由于相位角不為零，那么它所表示的濾波將會改變初始時間序列的相位，所以實際上在濾波處理過程中，相位失真總是不可避免的.由于正向和反向濾波的相位相反可以起到相互抵消的作用，從而達到修正相位失真的目的.這也是多重濾波時常使用正反交替的原因.但是由于在每次濾波過程中通過濾波器的信號需要從污染總負荷中減去以得到點源污染負荷，因此采用反向濾波也不能夠完全修正由于正向濾波所產(chǎn)生的相位失真［17］.

2 實例計算

以東江流域為研究對象，采用數(shù)字濾波法進行點源與非點源污染負荷分割，并與平均濃度法［19］的計算結(jié)果進行比較，探討濾波參數(shù)和濾波次數(shù)的選擇對于點源和非點源污染負荷分割結(jié)果的影響.

東江發(fā)源于江西省尋烏縣椏髻缽，上游稱尋烏水，南流入廣東省境內(nèi)，至龍川合河壩匯安遠水后稱東江.東江是珠江流域的主要支流之一，東江與西江、北江和珠江三角洲組成珠江.東江流域是香港、深圳、東莞、廣州、惠州、河源等城市的重要水源，其水環(huán)境質(zhì)量影響到上述城市的供水質(zhì)量.由于東江流域的重要性，從2000年開始廣東省水文局在博羅水文站對東江流域水質(zhì)水量進行同步監(jiān)測.根據(jù)博羅水文站2000—2005年同步水質(zhì)水量月系列資料，對東江流域污染物負荷進行點源和非點源污染負荷的分割，分析東江流域下游控制站博羅水文站的點源和非點源污染負荷，為東江流域的水環(huán)境整治提供理論支持.

筆者收集到博羅站2000—2005年水量月系列資料及氨氮，CODMn，總磷，BOD5等4種污染物同步監(jiān)測資料.僅以CODMn負荷月系列資料為例，應(yīng)用數(shù)字濾波法進行點源和非點源污染負荷的分割計算.

從1.4節(jié)分析知道，不同的濾波參數(shù)和濾波次數(shù)對點源和非點源分割結(jié)果有很大的影響.以東江流域為研究對象，根據(jù)經(jīng)驗，選擇 β分別為0.800，0.850，0.900，0.925和0.950的不同濾波次數(shù)(1～3次)，采用式(14)和(15)對 CODMn負荷月系列資料進行點源和非點源分割.表1為數(shù)字濾波法分割CODMn負荷點源的月均值與平均濃度法［19］的比較.從表1可以看出，數(shù)字濾波法分割CODMn負荷時，不同的濾波參數(shù)和不同濾波次數(shù)可以達到基本相同的效果.

對β為0.800進行4次濾波的CODMn負荷點源分割曲線，與 β為0.950進行2次濾波的 CODMn負荷點源分割曲線進行比較(見圖5).由圖5可見，兩曲線基本吻合，說明不同濾波參數(shù)經(jīng)過一定的濾波次數(shù)后，其曲線線形基本能夠吻合.

表1 數(shù)字濾波法和平均濃度法COD Mn負荷點源分割的月均值比較Table 1 Comparison of differentiation results between digital filter method and average concentration with respect to point CODMn pollutants t/月

為了提高濾波的效率，采用較大濾波參數(shù)(β=0.925)，用式(14)和(15)對博羅水文站 CODMn負荷月系列資料進行點源和非點源污染負荷的5次濾波分割計算，并與平均濃度法［9］的計算結(jié)果進行比較，結(jié)果見表2.由表2可見，采用比較大的濾波參數(shù)(β=0.925)，只需要3次濾波就可以近似達到平均濃度法的平均結(jié)果.

β為0.925的3次濾波CODMn負荷的分割結(jié)果見圖6和7.

圖5 β為0.950的2次濾波和β為0.800的4次濾波CODMn點源分割結(jié)果對比Fig.5 Comparison of differentiation of CODMn pollutants after two filtering practices withβof 0.950 and four filtering practices withβof 0.800

表2 COD Mn負荷點源分割的數(shù)字濾波法(β=0.925)和平均濃度法比較Table 2 Comparison between digital filter differentiation method(β=0.925)and average concentration method in terms of differentiation of non-point CODMn pollutants t

圖6 CODMn點源和非點源污染負荷數(shù)字濾波分割結(jié)果(β=0.925)Fig.6 Final results of filter-based differentiation of point and non-point CODMn pollutants(β=0.925)

圖7 CODMn點源負荷數(shù)字濾波最終分割曲線(β=0.925)Fig.7 Final results of filter-based differentiation of point CODMn pollutants(β=0.925)

由圖6可見，數(shù)字濾波法分割出來的CODMn點源負荷曲線趨勢平滑，波動幅度不大;而非點源負荷曲線由于受降水徑流影響，波動劇烈.從圖6看出，數(shù)字濾波法能將點源和非點源污染進行很好地分割，且分割出來的點源負荷和非點源負荷曲線比較符合其各自的產(chǎn)生特點，因此認為數(shù)字濾波法分割出來的點源負荷和非點源負荷是合理的.

由圖7可見，CODMn點源負荷從2000年1月到2005年3月(第63個月)是逐漸減少的，波動幅度基本維持在920～2 200 t/月，但從2005年初(第64個月)開始，CODMn點源月負荷開始增加，呈上升趨勢.這主要是廣東省在2004—2005年實施的產(chǎn)業(yè)轉(zhuǎn)移政策使部分產(chǎn)業(yè)從珠江三角洲轉(zhuǎn)移到了東江的中上游地區(qū)，造成了東江流域 CODMn點源污染在2005年出現(xiàn)上升趨勢.

根據(jù)數(shù)字濾波的計算結(jié)果，2000—2005年東江流域CODMn負荷主要來自于非點源污染(見表3).

表3 COD Mn負荷數(shù)字濾波的最終計算結(jié)果及點源和非點源所占比例Table 3 Differentiation results of point and non-point CODMn pollutants

3 結(jié)論

a.數(shù)字濾波法是根據(jù)污染負荷監(jiān)測系列資料進行污染負荷解析的有效方法，可用于點源負荷和非點源負荷的分割.

b.數(shù)字濾波法分割的點源和非點源負荷系列曲線符合點源負荷和非點源負荷產(chǎn)生特點.

c.數(shù)字濾波法具有簡單方便以及可操作性強等優(yōu)點，而且能夠再現(xiàn)分割的結(jié)果，使用者能夠在污染負荷分割的過程中通過濾波參數(shù)的選取比較方便地加入自己的經(jīng)驗.一般取較大濾波參數(shù)只需要3次濾波就能得出較滿意的結(jié)果.

d.實例計算表明，東江流域 2000—2005年CODMn負荷主要來自于非點源污染.

［1］ FUREY P R，GUPTA V K.A physically based filter for separating baseflow from streamflow time series ［J］.Water Resour Res，2001，37(11)：2709-2722.

［2］ FUREY P R，GUPTA V K.Tests of two physically based filters for baseflow separation［J］.Water Resour Res，2003，39(10)：1297-1307.

［3］ LIN K，GUO S，ZHANG W.A new baseflow separation method based on analytical solutions of the Horton infiltration capacity curve［J］.Hydrological Processes，2007，21(13)：1719-1736.

［4］ 林凱榮.數(shù)字水文模擬與基流分割方法研究［D］.武漢：武漢大學(xué)，2007.

［5］ 林凱榮，郭生練，張文華.基于霍頓下滲能力曲線的流量過程線連續(xù)分割方法研究［J］.水文，2008，28(1)：10-14.

［6］ 楊桂蓮，郝芳華，劉昌明，等.基于 SWAT模型的基流估算及評價：以洛河流域為例［J］.地理科學(xué)進展，2003，22(5)：463-471.

［7］ 陳利群，劉昌明，楊聰，等.黃河源區(qū)基流估算［J］.地理研究，2006，25(4)：659-665.

［8］ LYNE V.Stochastic time-variable rainfall-runoff modeling［M］//HOLLICK M.Hydrology and Water Resources Symposium.Perth，Australia： National Committee on Hydrology and Water Resources of the Institute of Engineering，1979：89-93.

［9］ NATHAN R J，MCMAHON T A.Evaluation of automated techniques for baseflow and recession analyses［J］.Water Resour Res，1990，26(7)：1465-1473.

［10］ ARNOLD J G，ALLEN P M，MUTTIAH R，et al.Automated baseflow separation and recession analysis techniques［J］.Ground Water，1995，33：1010-1018.

［11］ MAU D P，WINTER T C.Estimating ground-water recharge from streamflow hydrographs for a small mountain watershed in a temperate humid climate，New Hampshire，USA［J］.Ground Water，1997，35：291-304.

［12］ SPONGBERG ME.Spectral analysis of baseflow separation with digital filters［J］.Water Resour Res，2000，36(3)：745-752.

［13］ BURRUS C S，T W PARKS T W.DFTxA1/FFT andconvolut algorithms［M］.New York： Wiley，1985：48-60.

［14］ NEWTON H J.TIMESLAB： a time series analysis laboratory［M］.California：Wadsworth，Belmont，1988：75-80.

［15］ PRESS W H，TEUKOLSKY S A，VETTERING W T，et al.Numerical recipes in FORTRAN［M］.New York：Cambridge University Press，1990：465-1473.

［16］ 王麗萍，朱玉麗，王春，等.江蘇省環(huán)境負荷預(yù)測分析［J］.環(huán)境科學(xué)研究，2008，21(1)：208-212.

［17］ 袁宇，朱京海，侯永順，等.污染物入海通量非點源貢獻率的分析方法［J］.環(huán)境科學(xué)研究，2008，21(5)：169-172.

［18］ 林凱榮，陳曉宏，江濤，等.數(shù)字濾波進行基流分割的應(yīng)用研究［J］.水力發(fā)電，2008(6)：28-30.

［19］ 李懷恩.估算非點源污染負荷的平均濃度法及其應(yīng)用［J］.環(huán)境科學(xué)學(xué)報，2000，20(4)：397-400.

Application of a Digital Filter in the Differentiation of Pollutant Loads from Point and Non-Point Sources

LI Kun，LIN Kai-rong，JIANG Tao，CHEN Xiao-hong
Center for Water Resources and Environment Research，Sun Yat-sen University，Guangzhou 510275，China

This study aims to differentiate pollutant loads originating from point and non-point pollutant sources.Based on Fourier analysis techniques，a new method，digital filter，was developed with the aim to distinguish whether the pollutants are of point source or of non-point source through the time sequence of pollutant load.The relationship between parameters of filter equations and number of filter computations was discussed.This technique was successfully applied to the differentiation of pollutant loads in terms of CODMnin the East River basin.The results of this case study indicate that this newly-developed method has high maneuverability and is easy to use with a simple algorithm.The differentiated pollutants can well match the actual conditions.Additionally，the parameters of the newly-developed technique can be altered to well capture the changing properties of pollutants originating from point or non-point sources，which indicates larger flexibility of this technique in practice.Generally，satisfactory results can be obtained after filtering only three times by selecting larger filter parameters.

non-point source;pollutant loads;digital filter;apportionment

X32

A

1001－6929(2010)03－0298－06

2009－05－30

2009－11－26

國家自然科學(xué)基金青年科學(xué)項目(50809078);水資源與水電工程科學(xué)國家重點實驗室開放基金項目(2008B043)

黎坤(1969－)，男，廣西玉林人，工程師，博士，主要研究水資源與環(huán)境，非點源污染，eeslk@mail.sysu.edu.cn.

(責(zé)任編輯：孔 欣)