CN110377869A - 一种海洋浮标污染源寻找方法 - Google Patents

Abstract

本发明公开了一种海洋浮标污染源寻找方法，其步骤是：(1)对污染水域的水污染源数据进行去噪滤波处理；(2)建立基于高斯分布的污染源运动模型：(3)利用参数估计方法建立误差分布函数；(4)利用遗传算法对污染源的中心位置和污染源的扩散强度进行参数估计。本发明提出的海洋浮标污染源寻找方法，使得浮标系统通过海洋水质水文数据能够精确寻污溯源，海洋污染源数据通过基于自相关衰减率和小波变化的滤波方法能够消除由随机噪声产生的毛刺、不稳定波动等现象，使海洋水质水文数据的精度提高。

Description

一种海洋浮标污染源寻找方法

技术领域

本发明涉及一种污染源寻找方法，具体是一种海洋浮标污染源寻找方法。

背景技术

海洋污染监测一直以来备受关注，以海洋水产养殖为例，水域环境的优劣直接影响到鱼、虾、蟹、海参等海洋水产品质量的高低。海洋污染除了赤潮等大规模污染，基本是以源头为中心向外辐射污染(如化学药剂泄露)，以往的浮标水质监测，仅仅是对单点区域的水质要素进行记录，没有对数据进行分析以寻污溯源，并且由于海洋波浪等原因，数据精度往往较低。

在目前研发中的国外水质监测系统，如美国Heliosware公司的EMNET系统和澳大利亚CSIRO公司的Fleck系统，通信速率低，体积大，成本昂贵。国内的水质监测系统研究，如刘兴国等提出的基于水质监测技术的水产养殖安全保障系统和王骥等设计的基于无线传感器网络海洋环境监测系统，实现了水质多参数实时传输，但其测得数据的精确度不高，且仅仅局限在采集数据阶段，并没有进一步寻找污染源位置，没有达到从源头治理海洋污染的效果。同时由于海洋水质水文数据往往会伴随大量的误差干扰，数据的滤波没有有效处理。

发明内容

发明目的：本发明的目的是提供一种海洋浮标污染源寻找方法，通过建立基于高斯分布的污染源模型精确确定污染源点源的位置。

技术方案：一种海洋浮标污染源寻找方法，包括以下步骤：

(1)取水流上的封闭曲面所围成污染水域，对其污染源数据进行去噪滤波处理；

(2)设第i个观测点的污染物浓度指数为：

C_i＝C(x_i,y_i),i＝1,2,…,n

式中x_i，y_i是观测点的水平位置，由于在无限大介质中扩散浓度集中在点源，设定污染扩散属于多元的正态分布，若扩散模式解具有高斯解形式，即当有N个污染源时，(x_i,y_i)处的浓度分布为：

其中，k_j表示污染源的权重系数，Σ_j表示污染源的扩散强度，污染源的中心位置

(3)设定目标函数为：表示估计的污染分布与实际分布偏差最小；

(4)利用遗传算法对污染源的中心位置N个污染源的扩散强度σ_i进行参数估计，设置待估计参数的估计区间，设置遗传算法初始种群数和遗传代数，进行遗传交叉和变异等算法，以为适应度函数，选择算法执行过程中适应度函数最低的解，作为最终的污染源参数估计。

有益效果：(1)本发明提出的基于高斯分布的污染源寻找方法，使得浮标系统通过海洋水质水文数据精确寻污溯源，试验结果显示污染源位置的误差与研究区域的空间尺寸之比小于1％。

(2)本发明在海洋污染源数据在EEMD的基础上建立基于自相关衰减率和小波变化的滤波消除由随机噪声产生的毛刺、不稳定波动等现象，使海洋水质水文数据的精度提高。其滤除噪声效果明显比传统的低通滤波和中值滤波要好。

附图说明

图1是本发明的方法流程图；

图2是本发明水文水质数据的原始信号图；

图3是本发明的经验模态分解图；

图4是本发明基于自相关衰减率和小波变换的EEMD滤波方法滤波后信号图；

图5是本发明低通滤波信号图；

图6是本发明中值滤波信号图；

图7是本发明的静态水面污染源定位图；

图8是本发明的遗传算法定位污染源图；

图9是本发明的动态水面污染源定位图。

具体实施方式

如图1所示，一种海洋浮标污染源寻找方法，包括以下步骤：

(1)取水流上的封闭曲面所围成污染水域，对其污染源数据进行去噪滤波处理；

(2)设第i个观测点的污染物浓度指数为：

C_i＝C(x_i,y_i),i＝1,2,…,n

式中x_i，y_i是观测点的水平位置，由于在无限大介质中扩散浓度集中在点源，设定污染扩散属于多元的正态分布，若扩散模式解具有高斯解形式，即当有N个污染源时，(x_i,y_i)处的浓度分布为：

其中，k_j表示污染源的权重系数，Σ_j表示污染源的扩散强度，污染源的中心位置

(3)设定目标函数为：表示估计的污染分布与实际分布偏差最小；

(4)利用遗传算法对污染源的中心位置N个污染源的扩散强度σ_i进行参数估计，设置待估计参数的估计区间，设置遗传算法初始种群数和遗传代数，进行遗传交叉和变异等算法，以为适应度函数，选择算法执行过程中适应度函数最低的解，作为最终的污染源参数估计。

在步骤(1)中，对污染数据进行滤波处理的方法有很多种，如emd、eemd滤波方法，本发明提出一种基于自相关衰减率和小波变换的EEMD滤波方法。首先EMD方法将一般信号s(t)分解成若干个本征信号函数imf和一个残余信号r(t)组成，即：

而为了弥补EMD方法的模态混叠问题，现有技术中还提出了EEMD，一种噪声辅助信号处理方法。在EMD方法中，得到合理IMF的能力取决于信号极值点的分布情况，如果信号极值点分布不均匀，会出现模态混叠的情况。为此，EEMD方法将白噪声加入待分解信号，利用白噪声频谱的均匀分布，当信号加在遍布整个时频空间分布一致的白噪声背景上时，不同时间尺度的信号会自动分布到合适的参考尺度上，并且由于零均值噪声的特性，经过多次平均后，噪声将相互抵消,集成均值的结果就可作为最终结果。EEMD滤波方法的具体步骤如下:(1)向信号加入正态分布白噪声；(2)将加入白噪声的信号分解成各IMF分量；(3)重复步骤(1)和(2)，每次加入新的白噪声序列；(4)将每次得到的IMF集成均值作为最终结果。

在EEMD滤波方法的基础上，本发明提出了一种基于自相关衰减率和小波变化的滤波方法：基于的去噪方法主要目的是寻找到噪声主导模态和信号主导模态的分界点k。首先,此处引入随机信号的自相关函数统计特性。随机信号的自相关函数是信号时域特性的一种平均度量,它反映了信号在不同时刻,和取值的相关程度。自相关的直观含义就是：把一个信号平移一段距离，跟原来有多相似，则有表达式如下：

为了准确表现随机信号在不同时刻取值的相关程度,本文使用归一化自相关函数,即：

其中，R_x(τ)表示信号x(t)的自相关函数。

对于随机噪声由于其各个时刻的弱关联性及随机性,决定了其在零点处自相关函数值最大,在其它点处自相关函数值迅速衰减到很小的特点。对于理想高斯白噪声而言,其归一化自相关函数值在零点处为一,在其它点为零。对于一般信号而言,其自相关函数在零点处取得最大,由于信号间存在着关联性,因此在其他点处,自相关函数并没有迅速衰减到很小的值,而是随着时间差的变化而变化,变化规律明显有别于噪声的自相关函数的变化情况。

基于以上的特点,本发明首先根据各的自相关函数的特点,提出了一种基于自相关函数衰减率的分选有用信号与噪声分量分界点的准则，再对于我们所判定的噪声起主导作用的模态中,其组成成分除了噪声外,也会有少量有用信号的高频部分包含其中,此时,对噪声模态分量选小波变换去噪方法将要去除的分量进行滤波最后将所有的经过处理及未处理的分量进行重构得到去噪后的信号。

本发明所提基于自相关衰减率和小波变化的EEMD滤波方法步骤如下：

(1.1)对含噪声的信号x(t)进行EEMD分解，得到N个模态分量imf_i(i＝1,2,…,N)

(1.2)分别计算各个模态分量imf_i的归一化后的自相关函数

(1.3)定义衰减率β，其中ω根据实际信号取定，并根据实际信号设定阈值W，基于β和W判断噪声与信号起主导作用的模态分界点k；

(1.4)对判定为噪声起主导的模态imf₁～imf_k，采用小波软阈值去噪，得到去噪后的各分量imf₁’^w～imf_k’^w；

(1.5)重构信号， 即为滤波后的信号

其中，小波去噪软阈值函数为：

其中，ε的取值为，当在第j层时，N为信号长度，σ_j时信号在第j层的标准差。

将发明中模态分量带入，可得：

其中，ε_j为第j个imf_j的阈值。

图2-4是以标准正弦信号叠加随机信号为输入信号的仿真结果，模态分量总数为8，取k＝3，将自相关系数最大的三个模态分量进行叠加，其余进行小波软阈值去噪再叠加，最后得到滤波后的输出信号。可以看到原始信号中由于随机噪声产生的毛刺、不稳定波动等现象，在滤波后得到很好的滤除。

图5和图6表示与传统的低通滤波和中值滤波的效果，可以看出低通滤波虽然有效的抑制了高频噪声，但是对原始信号也有一定影响，而中值滤波对于噪声的滤除效果明显不如本文方法。

基于以上方法，可以自适应的从经验模态分解中提取合适信号，组合为滤波后的信号。

在步骤(2)中的所设立的污染源寻找模型是通过以下推论得出：

首先对水域模型进行模型简化处理，任取水流上以封闭曲面S所围成的水域Ω，C(x,y,z,t)表示t时刻(x,y,z)位置点的污染物浓度；通常污染物因为河流流动和分子自由运动而发生扩散，D_x，D_y，D_z分别表示x，y，z方向的扩散系数；因为水流都有吐故纳新和自我净化的能力，假设自我降解系数为K；θ(x,y,z,t)是(x,y,z)点在t时刻单位体积单位时间污染物的排放量。

由多重积分物理意义可知，通过封闭曲面S从时刻t到时刻t+dt流入Ω的污染物质量为

其中，cosα，cosβ，cosγ为S的外法向余弦。

由高斯定理可知：

由于河水自净作用，从t到t+Δt时刻Ω水域的污染物减少为：

由于河水不断流动，污染物从t到t+Δt时刻流出S的质量为：

同理，由高斯定理知：

封闭曲面S内从t到t+Δt时刻污染物的排放量为：

从另一个角度看，由于浓度的变化引起Ω内污染物的质量的增加量为：

由质量守恒定律得到：

M₅＝M₁-M₂-M₃+M₄

所以河流水质污染的4D数学模型为：

初始条件为：

C(x,y,z,0)＝λ(x,y,z)

对四维水质模型左右两边同时进行傅里叶变换得到：

合并同类项得：

解此线性微分方程得：

再对上式左右同时取傅里叶逆变换得：

从上式中可以看到C(x,y,z,t)的表达式十分复杂，实际上，现实水流中绝大部分情况下用不到这样面面俱到的计算模型，只需要抓住问题的本质即可。考虑瞬时污染点源四维扩散模型：

其中，-∞<x<+∞，-∞<y<+∞，-∞<z<+∞，t>0。

初始条件为：

C(x,y,z)＝Mδ(x)δ(y)δ(z)

边界条件为：

该模型在三维水体中，于原点(0,0,0)处瞬时投入质量为M的污染源，污染源在x，y，z方向上进行扩散，其解析解为：

z_i即为湖面深度，恒为0。假设河水的垂直对流速度远小于水平平流速度，则上式可继续简化为二维问题：

由于在无限大介质中扩散浓度集中在点源，这类扩散问题表达式都具有高斯分布的形式，即“高斯解”，由此，假设污染扩散是多元的正态分布。

设一个复合函数g(μ)，根据扩散模式呈高斯分布形态，g(μ)符合正态分布。多元正态分布是一元正态分布的直接推广。设随机向量X(x₁,x₂,…,x_m)服从p维正态分布，则有：

在研究污染问题时，g(μ)可简化为：

式中k为污染源的权重系数。

设第i个观测点的污染物浓度指数为：

C_i＝C(x_i,y_i),i＝1,2,…,n

式中x_i，y_i是观测点的水平位置。

若扩散模式解具有高斯解形式，扩散函数应该符合多元正态分布函数形式，即当只有一个污染源时，(x_i,y_i)处的浓度分布为：

当有N个污染源时，(x_i,y_i)处的浓度分布为N个污染源在此处污染效果的叠加：

其中，k_j，∑_j，都是待估计的参数。

在步骤(4)目标函数为：表示估计的污染分布与实际分布偏差最小。利用遗传算法，可以对以上参数进行估计。首先设定污染源个数N,作为混合二元正态分布的变量个数，其中待估计的参数有，N个污染源的中心位置N个污染源的扩散强度σ_i，σ_i是污染源本身强度与时间的函数，但在某一固定时间下，或者当时间t→+∞时，σ_i仅仅时污染源本身强度的函数。因此，有3N个待估计参数。其次设置待估计参数的估计区间，由问题本身研究的区域决定，如在上述模拟过程中，0≤σ_i≤1。最后设置遗传算法初始种群数和遗传代数，进行遗传交叉和变异等算法。以为适应度函数，选择算法执行过程中适应度函数最低的解，作为最终的污染源参数估计。

如图7和8所示，图7为污染源分布图，图8中正方形标识出的是利用遗传算法定位出的污染源位置，而三角形标识出的是设定的污染源位置，污染源位置的误差与研究区域的空间尺寸之比小于1％。

对于河流存在一定流速的情况下，也可以实现上述算法。其理论为，对于连续污染源，河流中任意位置的污染值可以由以下积分表示

如果当时间足够长的时候，污染物浓度可以认为不再随时间而变化，仅仅是空间位置的函数，其稳态模型的解析解为

因此，可模拟连续污染源在流动河水中的污染情况。如图9所示模拟在水平方向上河流由一定流速时，连续污染源在经过一定时间后达到的稳定污染情况，虽然这时候污染情况显然符合高斯解，但仍然可以利用参数估计的方法，以最小化为目标函数，利用遗传算法寻找污染源位置。最终，寻找到的污染源用正方形标识，设定的污染源实际位置用三角形标识在图中，污染源位置的误差与研究区域的空间尺寸之比小于2％。

Claims (5)

1.一种海洋浮标污染源寻找方法，其特征在于包括以下步骤：

(1)取水流上的封闭曲面所围成污染水域，对其污染源数据进行去噪滤波处理；

(2)设第i个观测点的污染物浓度指数为：

C_i＝C(x_i,y_i),i＝1,2,…,n

式中x_i，y_i是观测点的水平位置，由于在无限大介质中扩散浓度集中在点源，设定污染扩散属于多元的正态分布，若扩散模式解具有高斯解形式，即当有N个污染源时，(x_i,y_i)处的浓度分布为：

其中，k_j表示污染源的权重系数，Σ_j表示污染源的扩散强度，污染源的中心位置

(3)设定目标函数为：表示估计的污染分布与实际分布偏差最小；

(4)利用遗传算法对污染源的中心位置N个污染源的扩散强度σ_i进行参数估计，设置待估计参数的估计区间，设置遗传算法初始种群数和遗传代数，进行遗传交叉和变异算法，以为适应度函数，选择算法执行过程中适应度函数最低的解，作为最终的污染源参数估计。

2.根据权利要求1所述的海洋浮标污染源寻找方法，其特征在于：步骤(1)所述对污染源数据进行去噪滤波处理的方法为基于自相关衰减率和小波变化的EEMD滤波方法，具体包括以下步骤：

(1.1)对含噪声的信号x(t)进行EEMD分解，得到N个模态分量imf_i(i＝1,2,…,N)；

(1.2)分别计算各个模态分量imf_i的归一化后的自相关函数

(1.3)定义衰减率β:其中ω根据实际信号取定，并根据实际信号设定阈值W，基于衰减率β和阈值W判断噪声与信号起主导作用的模态分界点k；

(1.4)对判定为噪声起主导的模态imf₁～imf_k，采用小波软阈值去噪，得到去噪后的各分量imf₁’^w～imf_k’^w；

(1.5)重构信号， 即为滤波后的信号。

3.根据权利要求2所述的海洋浮标污染源寻找方法，其特征在于：步骤(1.4)所述小波去噪软阈值函数为：

其中，ε的取值为，当在第j层时，N为信号长度，σ_j时信号在第j层的标准差。

4.根据权利要求1所述的海洋浮标污染源寻找方法，其特征在于：步骤(1)中污染水域的水流雷诺数达到大于2000Rep的阈值，水流的流动方式为紊流而不是层流。

5.根据权利要求1所述的海洋浮标污染源寻找方法，其特征在于：步骤(1)中污染水域的水流沿(x,y,)两方向的流动速度u_x，u_y在一定时间内是恒定的。