数字滤波法在点源和非点源污染负荷分割中的应用

黎 坤，林凯荣，江 涛，陈晓宏

中山大学水资源与环境研究中心，广东 广州 510275

数字滤波法在点源和非点源污染负荷分割中的应用

黎 坤，林凯荣，江 涛，陈晓宏

中山大学水资源与环境研究中心，广东 广州 510275

以分辨流域点源和非点源污染负荷为研究目标，在傅立叶分析的基础上，提出了从污染负荷时间序列中分辨点源和非点源污染负荷的新方法——数字滤波法，探讨了数字滤波方程的参数与滤波次数的关系，并将该方法应用到东江流域CODMn负荷的分割上.应用实例表明，分割的点源和非点源负荷系列曲线符合点源和非点源负荷的产生特点，使用者能够在污染负荷分割过程中通过滤波参数的选取比较方便地加入自己的经验.一般取较大滤波参数时，只需要3次滤波就能得出较满意的结果.

非点源;污染负荷;数字滤波;解析

由于非点源具有随机性、广泛性和难监测性，在我国对非点源的研究还处于起步阶段.目前的水质监测主要是在可控制的河流、湖泊、水库的进出口断面上，无法区分点源和非点源污染负荷.如何在现有资料的基础上来分辨点源和非点源负荷，为污染源的解析提供依据，成为目前应该解决的一个关键问题.笔者在傅立叶分析的基础上，提出了从污染负荷时间序列中分辨点源和非点源污染负荷的新方法——数字滤波法.

数字滤波法是近年来国际上研究最多的分割基流的方法［1-7］.LYNE［8］于 1979 年首次提出用该方法来分割水文过程中径流对降雨的快速响应和慢速响应成分.NATHAN等［9］研究认为，该方法是一种快速而且客观的基流分割方法，同时发现滤波参数在0.900～0.950内，能够得到较好的分割结果.ARNOLD等［10］把数字滤波法与其他基流分割方法进行了分析比较，结果表明，采用3次滤波且滤波参数取为0.925能够得到较好的分割结果;同时认为该方法很容易使用，而且能够再现分割的结果;使用者能够在分割基流的过程中比较方便地加入自己的经验.MAU 等［11］则发现，滤波参数取为 0.850，进行多次滤波可以得到最好的结果.SPONGBERG［12］在傅立叶分析的基础上，阐述数字滤波分割基流的基本原理，并在时间域和频率域上对数字滤波法分割基流的结果进行了分析验证.但鲜见应用数字滤波法来研究非点源的解析，笔者选择应用数字滤波法分辨点源和非点源污染负荷通量.

1 数字滤波法基本原理

数字滤波法源自傅立叶分析，所以在陈述数字滤波法基本原理时有必要了解傅立叶分析原理.

1.1 傅立叶分析

傅立叶分析的研究与应用至今已经历了100余年.傅立叶分析主要是通过傅立叶变换由时域分析转入变换域分析.傅立叶变换是在以时间为自变量的“信号”与以频率为自变量的“频谱”函数之间的某种变换关系，即建立了时间函数和频谱函数之间的转换关系.

1.1.1 傅立叶变换

如果x(t)是定义在整个实轴上的实值或复值函数，则其傅立叶变换可由下式给出：

若对任意参数f，上述积分都存在，则式(1)确定了函数X(f)，称为X(f)的傅立叶变换.如果已知X(f)，则利用如下的傅立叶逆变换，还可复原x(t)：

若x(t)和X(f)同时满足式(1)和(2)，则称它们是一个傅立叶变换对，记为 x(t)⇐X(f).通常X(f)是一个复函数，因此可以写成如下两部分：

式中，R(f)，I(f)分别是 X(f)的实部和虚部.将式(3)表示为指数形式：

其中，

工程技术中，常将x(t)看成时间信号，相应的空间，称为时间域和空域;将其傅立叶变换 X(f)看成频率函数，相应的空间称为频域. X(f )称为x(t)的傅立叶谱，而 φ(f)称为其相角，这在物理上是有良好背景的.该频率的含义可以这样来理解：应用欧拉公式可将指数项表示成正弦－余弦的形式，如果把式(1)解释成离散项和的极限，则显然X(f)是包含了无限项正弦－余弦的和，而且 f的每一个值确定了所对应的正弦－余弦的频率.

1.1.2 离散傅立叶变换

实际上，水文时间序列常以离散的形式给出，因此常采用的是离散傅立叶变换(DFT).离散傅立叶变换的一般形式如下：

给定 N 个 实 或 复 的 数 列{x(0)，x(1)，…，x(N － 1) } ，定义

X(n)为{x(k)}的离散傅立叶变换.

离散的频率可表示为：

如果不考虑式(7)中的复指数部分的运算，则求解式(7)共需要N×N次乘法和N×(N－1)次加法.显然当N很大时，其工作量是相当可观的.因此，通常采用凯莱和卡柯［13］提出的专门用于处理DFT的快速算法(FFT)，可减少DFT的计算时间.

1.2 滤波原理

数字滤波系统的基本输入输出关系为：

式中，yk为滤波输出项;ck为滤波系数;xk为输入项;*表示卷积.

卷积原理表明，在时域上的卷积等价于在频域上的乘积，故有：式中，T(n)为ck的离散傅立叶变换.

式(10)表明，滤波输出项(时间序列)的离散傅立叶变换，变换过程如1.1节所示，可以表示为2个复数的频率序列的乘积，用极数的形式可表示为：

式中，z为复数;r为复数的模(幅值);θ为相位角.

从式(10)和(11)可以看出，初始序列的滤波输出项的大小取决于T(n)的幅值，为此，把T(n)称为滤波系数ck的频率转移函数，简称FTF［14］.如果FTF为复数，其相位角不为零，那么它所表示的滤波将会改变初始时间序列的相位角.

线性递归滤波的一般形式为：

式中，dm为滤波系数;J和M分别为cj和dm个数的总和.

由式(12)表示滤波的频率转移函数(FTF)可表示为［15］：

1.3 频谱分析

污染负荷是指区域或某环境要素对污染物的负载量［16］，目前对污染负荷的监测是通过水量水质同步监测来实现的，即采集水样的同时进行流量观测，其监测的是污染负荷随时间的非连续变化过程，是某断面的污染负荷的通量，流域某断面的污染负荷通量是点源负荷和非点源负荷的总和［17］.

点源污染是指有固定排放点的污染源，如工业废水及城市生活污水等，由排放口集中汇入水体引起的污染;非点源污染是指溶解性或固体污染物在大面积降水和径流冲刷作用下汇入受纳水体而引起的污染.点源受降水等自然因素的影响较小，其排放过程相对稳定.而非点源污染受降水等因素影响，排放呈现随机性，波动较大.

一般某流域出口断面点源污染负荷的出流时间序列是相对稳定的序列.对点源污染负荷的出流时间序列进行傅立叶变换，其变换后的幅值波形图见图1.由图1可见，点源污染负荷的出流时间序列频域下的信号主要是低频信号.

图1 点源污染负荷的时间序列傅立叶变换后的幅值波形图Fig.1 Power spectrum of point pollutant series after Fourier transform analysis

非点源污染负荷的出流时间序列由于受降水等因素影响波动很大，对非点源污染负荷的出流时间序列进行傅立叶变换，其变换后的幅值波形图见图2.由图2可见，非点源污染负荷的出流时间序列频域下的信号明显比点源污染负荷出流的信号要大得多，且主要是高频信号，这也就是可以利用分离高低频信号进行点源和非点源污染负荷系列分割的理论依据.

但是，从图2还可以看出，非点源污染负荷的时间序列在频域下的信号并不能完全看作是高频信号，因为它还包括了部分低频信号，这就意味着要想完全地分割点源和非点源污染负荷是不可能的，因为有一部分频率是交叠的.污染负荷系列中低频率信号的分离会部分地削弱点源污染负荷.但是，滤波的频率转移函数(FTF)提供了能够得到最好分割点源和非点源污染负荷时间序列的途径.通过频率转移函数选择最佳的滤波系统使得在最小削弱点源负荷的同时最大限度地削减非点源污染负荷，而且尽可能地消除相位失真.

图2 非点源污染负荷的时间序列傅立叶变换后的幅值波形图Fig.2 Power spectrum of non-point pollutant series after Fourier transform analysis

1.4 数字滤波算法

笔者将数字滤波法引入到流域点源和非点源污染负荷的分割上.参考LYNE［8］于1979年提出的滤波方程，点源和非点源污染负荷分割滤波方程为：

式中，Pt为t时段内过滤后的污染负荷出流量(非点源负荷出流量)，kg，计算中 P1设为0;Bt为 t时段的污染负荷出流总量，kg;β为无量纲参数，称为滤波参数.已知污染负荷出流总量和非点源污染负荷出流量后，点源污染负荷出流量(Nt)可由下式求出：

由式(14)和(15)可试算出点源和非点源污染负荷的出流过程，从而对污染负荷的出流过程线进行点源和非点源出流负荷分割.

滤波参数(β)需要优化率定.对于滤波参数和滤波次数，不同的学者有着不同的认识.下面通过频率转移函数来分析LYNE［8］提出的数字滤波法的滤波参数和滤波次数的关系.

把式(14)转化为式(12)的形式，可得：

且有 J=1，M=1.

根据式(12)，可得 LYNE［8］数字滤波的频率转移函数为：

由式(17)通过Matlab软件转换，可以得到频率转移函数的幅度谱和相位谱，如图3，4所示.

图3 不同滤波参数下LYNE滤波的FTF的幅值波形图Fig.3 Power spectrum curves of FTF based LYNE filter with differentβ

图4 不同滤波参数下LYNE滤波的FTF的相位谱Fig.4 Phase curves of FTF based LYNE filter with differentβ

从理论上来说，非点源污染负荷信号将全部通过滤波，而点源污染负荷信号则不能通过.从图3可以看出，对所有比较常用的滤波参数来说，在最高频率范围中，FTF的幅值为1.0，包含全部非点源污染负荷信号.而多重滤波并不能进一步衰减高频率信号，因为在第一次滤波中它们已经完全衰减了.增大滤波参数将会增加点源污染负荷和非点源污染负荷信号中低频率部分的衰减量.使用较小的参数将会减少该频率范围内的衰减量.虽然这在FTF的振幅范围内不是很明显，但在实际应用中可以发现，如果选用比较小的参数，以上一次滤波结果作为下一次滤波输入，通常要滤波5～6次才能得到合理的点源和非点源负荷的分割.当使用较大滤波参数时，只需滤波2～3次就可以得到相同的分割量.

图4显示的是 β为 0.850，0.925和 0.975的LYNE滤波的FTF的相位谱.由于相位角不为零，那么它所表示的滤波将会改变初始时间序列的相位，所以实际上在滤波处理过程中，相位失真总是不可避免的.由于正向和反向滤波的相位相反可以起到相互抵消的作用，从而达到修正相位失真的目的.这也是多重滤波时常使用正反交替的原因.但是由于在每次滤波过程中通过滤波器的信号需要从污染总负荷中减去以得到点源污染负荷，因此采用反向滤波也不能够完全修正由于正向滤波所产生的相位失真［17］.

2 实例计算

以东江流域为研究对象，采用数字滤波法进行点源与非点源污染负荷分割，并与平均浓度法［19］的计算结果进行比较，探讨滤波参数和滤波次数的选择对于点源和非点源污染负荷分割结果的影响.

东江发源于江西省寻乌县桠髻钵，上游称寻乌水，南流入广东省境内，至龙川合河坝汇安远水后称东江.东江是珠江流域的主要支流之一，东江与西江、北江和珠江三角洲组成珠江.东江流域是香港、深圳、东莞、广州、惠州、河源等城市的重要水源，其水环境质量影响到上述城市的供水质量.由于东江流域的重要性，从2000年开始广东省水文局在博罗水文站对东江流域水质水量进行同步监测.根据博罗水文站2000—2005年同步水质水量月系列资料，对东江流域污染物负荷进行点源和非点源污染负荷的分割，分析东江流域下游控制站博罗水文站的点源和非点源污染负荷，为东江流域的水环境整治提供理论支持.

笔者收集到博罗站2000—2005年水量月系列资料及氨氮，CODMn，总磷，BOD5等4种污染物同步监测资料.仅以CODMn负荷月系列资料为例，应用数字滤波法进行点源和非点源污染负荷的分割计算.

从1.4节分析知道，不同的滤波参数和滤波次数对点源和非点源分割结果有很大的影响.以东江流域为研究对象，根据经验，选择 β分别为0.800，0.850，0.900，0.925和0.950的不同滤波次数(1～3次)，采用式(14)和(15)对 CODMn负荷月系列资料进行点源和非点源分割.表1为数字滤波法分割CODMn负荷点源的月均值与平均浓度法［19］的比较.从表1可以看出，数字滤波法分割CODMn负荷时，不同的滤波参数和不同滤波次数可以达到基本相同的效果.

对β为0.800进行4次滤波的CODMn负荷点源分割曲线，与 β为0.950进行2次滤波的 CODMn负荷点源分割曲线进行比较(见图5).由图5可见，两曲线基本吻合，说明不同滤波参数经过一定的滤波次数后，其曲线线形基本能够吻合.

表1 数字滤波法和平均浓度法COD Mn负荷点源分割的月均值比较Table 1 Comparison of differentiation results between digital filter method and average concentration with respect to point CODMn pollutants t/月

为了提高滤波的效率，采用较大滤波参数(β=0.925)，用式(14)和(15)对博罗水文站 CODMn负荷月系列资料进行点源和非点源污染负荷的5次滤波分割计算，并与平均浓度法［9］的计算结果进行比较，结果见表2.由表2可见，采用比较大的滤波参数(β=0.925)，只需要3次滤波就可以近似达到平均浓度法的平均结果.

β为0.925的3次滤波CODMn负荷的分割结果见图6和7.

图5 β为0.950的2次滤波和β为0.800的4次滤波CODMn点源分割结果对比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负荷点源分割的数字滤波法(β=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点源和非点源污染负荷数字滤波分割结果(β=0.925)Fig.6 Final results of filter-based differentiation of point and non-point CODMn pollutants(β=0.925)

图7 CODMn点源负荷数字滤波最终分割曲线(β=0.925)Fig.7 Final results of filter-based differentiation of point CODMn pollutants(β=0.925)

由图6可见，数字滤波法分割出来的CODMn点源负荷曲线趋势平滑，波动幅度不大;而非点源负荷曲线由于受降水径流影响，波动剧烈.从图6看出，数字滤波法能将点源和非点源污染进行很好地分割，且分割出来的点源负荷和非点源负荷曲线比较符合其各自的产生特点，因此认为数字滤波法分割出来的点源负荷和非点源负荷是合理的.

由图7可见，CODMn点源负荷从2000年1月到2005年3月(第63个月)是逐渐减少的，波动幅度基本维持在920～2 200 t/月，但从2005年初(第64个月)开始，CODMn点源月负荷开始增加，呈上升趋势.这主要是广东省在2004—2005年实施的产业转移政策使部分产业从珠江三角洲转移到了东江的中上游地区，造成了东江流域 CODMn点源污染在2005年出现上升趋势.

根据数字滤波的计算结果，2000—2005年东江流域CODMn负荷主要来自于非点源污染(见表3).

表3 COD Mn负荷数字滤波的最终计算结果及点源和非点源所占比例Table 3 Differentiation results of point and non-point CODMn pollutants

3 结论

a.数字滤波法是根据污染负荷监测系列资料进行污染负荷解析的有效方法，可用于点源负荷和非点源负荷的分割.

b.数字滤波法分割的点源和非点源负荷系列曲线符合点源负荷和非点源负荷产生特点.

c.数字滤波法具有简单方便以及可操作性强等优点，而且能够再现分割的结果，使用者能够在污染负荷分割的过程中通过滤波参数的选取比较方便地加入自己的经验.一般取较大滤波参数只需要3次滤波就能得出较满意的结果.

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］ 林凯荣.数字水文模拟与基流分割方法研究［D］.武汉：武汉大学，2007.

［5］ 林凯荣，郭生练，张文华.基于霍顿下渗能力曲线的流量过程线连续分割方法研究［J］.水文，2008，28(1)：10-14.

［6］ 杨桂莲，郝芳华，刘昌明，等.基于 SWAT模型的基流估算及评价：以洛河流域为例［J］.地理科学进展，2003，22(5)：463-471.

［7］ 陈利群，刘昌明，杨聪，等.黄河源区基流估算［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］ 王丽萍，朱玉丽，王春，等.江苏省环境负荷预测分析［J］.环境科学研究，2008，21(1)：208-212.

［17］ 袁宇，朱京海，侯永顺，等.污染物入海通量非点源贡献率的分析方法［J］.环境科学研究，2008，21(5)：169-172.

［18］ 林凯荣，陈晓宏，江涛，等.数字滤波进行基流分割的应用研究［J］.水力发电，2008(6)：28-30.

［19］ 李怀恩.估算非点源污染负荷的平均浓度法及其应用［J］.环境科学学报，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

国家自然科学基金青年科学项目(50809078);水资源与水电工程科学国家重点实验室开放基金项目(2008B043)

黎坤(1969－)，男，广西玉林人，工程师，博士，主要研究水资源与环境，非点源污染，eeslk@mail.sysu.edu.cn.

(责任编辑：孔 欣)