一种基于地质统计学的河流污染溯源方法

摘要

本发明公开了一种基于地质统计学的河流污染溯源方法，应用于环境污染技术领域，具体步骤如下：获得区面上的多时段水质与流速数据；利用地质统计学方法构建污染源溯源的反演数学模型；根据多时段水质与流速数据初步判断污染源所处的区域，并假定可能的污染源位置；利用流速数据建立分析河段的水动力及水质模拟模型，并通过水动力及水质模拟模型计算传递函数；采用GPU并行加速方法改进原最小化问题的计算方法求解污染源溯源的反演数学模型，通过比较模型计算结果与实测的水质与流速数据拟合程度，确定污染源的位置与排放过程。该方法为了提高计算效率采用GPU并行加速方法改进原最小化问题的计算方法求解污染源溯源的反演数学模型。

说明书

技术领域

本发明涉及环境保护技术领域，更具体的说是涉及一种基于地质统计学的河流污染溯源方法。

背景技术

由水污染导致的水质型缺水已经成为制约国民经济可持续发展的重要因素。水污染的来源主要有城市生活污水、农业面源污染、工业污水等。目前，我国不断加大对环境保护的力度，在水污染治理方面已经开展了大量工作，水环境质量在整体上有较大的提升。然而，仍存在部分偷排或污染物意外泄漏的情况。这些污染排放存在隐蔽、随机的特点。一旦发生，会对水环境产生较大的影响，但是由于其排放不规律等原因，难以快速确定污染源，为水环境的管理与治理带来了一定困难。

传统的污染源调查工作方法主要依据水质监测数据来判定污染源的位置与强度。水质监测数据主要来自于固定的监测站点与人工取样，这些方法虽然精确度高，但仅能提供有限点位的水质数据，耗时耗力，调查效率低，无法反应区面上的水质变化，也难以查找隐蔽、随机的污染源。基于无人机平台的高光谱或多光谱水质分析技术，侦测范围广、使用条件灵活、数据分辨率高，近年来已出现了大量的相关研究，为高效地获取较大区面上水质数据提供了技术基础。根据多时段、大区面上的水质数据，利用地质统计学方法反演污染源位置与排放过程，能够极大的节省甄别隐蔽污染源所需的人力物力，有利于提高污染源排查和相关环境保护工作的效率和准确性。

进一步，在溯源问题较为复杂时，由于计算负荷较大，需要的计算时间较长，对于计算溯源问题时，计算量比较大的问题是本领域技术人员亟需解决的问题。

发明内容

有鉴于此，本发明提供了一种针对偷排或污染物意外泄漏时污染源具有隐蔽、随机、难以确定的特点，提出一种新的污染源溯源方法。该方法先利用无人机平台的高光谱或多光谱水质分析技术为该方法获得区面上的多时段水质数据，再利用地质统计学方法反演出隐蔽污染源的位置与排放过程，同时为了提高计算效率采用GPU并行加速方法改进原最小化问题的计算方法求解污染源溯源的反演数学模型。

为了实现上述目的，本发明提供如下技术方案：

一种基于地质统计学的河流污染溯源方法，具体步骤如下：

获得区面上的多时段水质与流速数据；

利用地质统计学方法构建污染源溯源的反演数学模型；

根据多时段水质与流速数据初步判断污染源所处的区域，并假定可能的污染源位置；

利用流速数据建立分析河段的水动力及水质模拟模型，并通过水动力及水质模拟模型计算传递函数；

采用GPU并行加速方法改进原最小化问题的计算方法求解污染源溯源的反演数学模型，通过比较模型计算结果与实测的水质与流速数据拟合程度，确定污染源的位置与排放过程。

可选的，在上述的一种基于地质统计学的河流污染溯源方法中，构建污染源溯源的反演数学模型具体步骤如下：

假定河流中污染物未发生化学反应或发生一级化学反应，将污染物排放过程与河流中污染物浓度观测值间的关系概化为式1：

z＝Hs+v 式1；

式中：z为一个由观测数据组成的m×1维向量；s为待求未知函数经离散后形成的一个n×1维向量，代表污染物排放过程；v为测量误差组成的m×1维向量；H为由传递函数构成的灵敏度矩阵，代表了污染物在河水中的迁移过程；

假定v的均值为0，其协方差矩阵为R，s为一个随机向量，其数学期望和协方差分别由式2和式3表示：

E[s]＝Xβ 式2；

Q(θ)＝E[(s-Xβ)(s-Xβ)

式中：X为一个已知的n×p的矩阵；β为p个漂移系数组成的向量；Q(θ)为关于未知参数θ的已知函数，通常取高斯分布函数。

可选的，在上述的一种基于地质统计学的河流污染溯源方法中，根据污染物浓度的非负性，对未知函数s进行如下变换：

将式4代入式1得到

变换后，

可选的，在上述的一种基于地质统计学的河流污染溯源方法中，通过水动力及水质模拟模型计算传递函数具体步骤如下：

由于假定河流中污染物未发生化学反应或发生一级化学反应，污染物迁移过程用线性偏微分方程表示：

式中c为污染物浓度，u为流速，D为弥散系数，t为时间，x为迁移距离，r(c)为反应项；

此时，污染物浓度由传递函数方程表示：

其中s(t)为实际排放过程，f(·)为传递函数。

可选的，在上述的一种基于地质统计学的河流污染溯源方法中，通过差分方法计算传递函数，从而得到含有一级化学反应过程的传递函数；

通过变量替换将式8变为式9：

当式9中t>0时，s(T-t)＝1，t≤0时，s(T-t)＝0，得到穿透曲线方程：

使式10对时间t求导，即可得到传递函数在坐标x时刻t的值；

根据假定的污染源位置，利用正向的河流污染物迁移模型计算s(t)＝1条件下各观测点的浓度曲线，然后采用向后差分方法，式12计算各观测点在各个时刻的传递函数值，从而得到由传递函数所组成的灵敏度矩阵H，最后利用地质统计学方法反演出相应污染源污染物的排放过程；

可选的，在上述的一种基于地质统计学的河流污染溯源方法中，待求未知函数离散形式s的求解过程可分成两个步骤，首先是确定未知参数θ，而后估计未知函数s；

当式13表示的概率取最大值时，θ的取值即为待求的参数值；

Σ＝HQH

Ξ＝Σ

式13描述的极值问题，等同于求式16的最小值；

L(θ)对θ的导数取0值时，L(θ)取最小值；当迭代计算收敛后，再利用θ求Q，并代入式17，解得到n×m的系数矩阵Λ和p×n的拉格朗日乘子阵M；

对未知函数s的最佳估计

V＝-XM+Q-QH

得到假定污染源对应的排放过程。

可选的，在上述的一种基于地质统计学的河流污染溯源方法中，采用GPU并行加速方法改进原最小化问题的计算方法具体步骤如下：

要求L(θ

可选的，在上述的一种基于地质统计学的河流污染溯源方法中，将GPU中Grid中的若干个Block作为一个Block组，每个Block组执行一个λ

可选的，在上述的一种基于地质统计学的河流污染溯源方法中，将假定污染源与排放过程输入正向模型中，计算各离散网格的水质与流速，根据计算结果与测量结果的拟合效果确定污染源的位置与排放过程；采用下式描述拟合效果：

式中

经由上述的技术方案可知，与现有技术相比，本发明公开提供了一种基于地质统计学的河流污染溯源方法，充分利用了无人机平台的灵活性，以及光谱分析技术便于取得区面数据的优势，克服了区域水质取样密度低及取样频次低的不足，在水质与流速数据的基础上，利用地质统计学方法反演污染物在河水中的迁移过程，通过对比分析计算结果与测量数据的拟合程度，确定隐蔽污染源的位置与排放过程。与传统水质分析方法或人工探查方法相比，能够反映区面上的水质信息，便于查找隐蔽污染源，同时提高污染源调查工作的效率，有利于提高河流污染治理及管理的技术水平。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图获得其他的附图。

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

图2为本发明的并行加速示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明实施例公开了一种基于地质统计学的河流污染溯源方法，该方法先利用无人机平台的高光谱或多光谱水质分析技术为该方法获得区面上的多时段水质数据，再利用地质统计学方法反演出隐蔽污染源的位置与排放过程。具体步骤如下：

(1)利用无人机平台的高光谱或多光谱水质分析技术为该方法获得区面上的多时段水质与流速数据。

(2)利用地质统计学方法构建污染源溯源的反演数学模型。

(3)根据已知数据与相关调查初步判断污染源所处的区域，并假定可能的污染源位置。

(4)利用流速数据建立分析河段的水动力及水质模拟模型，并得通过模型计算传递函数。

(5)求解污染源溯源的反演数学模型，通过比较模型计算结果与实测的水质与流速数据拟合程度，分析确定污染源的位置与排放过程。

本发明实施例的的主要内容如下：

1.获取区面上的水质与流速数据

无人机平台和高光谱或多光谱水质分析技术的作用是为本发明提供区面上的水质参数与流速数据。目前，已有大量利用无人机平台和高光谱或多光谱手段分析水质与流速的研究，本发明此处仅简述水质与流速数据获取的过程。

(1)水质参数

首先，使用搭载有高光谱或多光谱成像设备的无人机，在拟分析河段上按一定航迹飞行，并获取飞行路线上的高光谱或多光谱低空遥感影像。其次，利用实测水质数据和遥感影像，确定水质参数与光谱曲线间的定量关系。再次，利用机器学习技术构建水质参数反演模型。最终，利用水质参数反演模型分析不同时段的影像数据，从而获得分析河段上相应时段的水质参数空间分布数据。

(2)流速参数

利用获得的影像数据，分析河水表面水流特征和漂浮物交错轨迹的流动时间计算流动速度，然后分析其分布并间接地计算水流速度，最终得到不同时段的流速空间分布数据。

(3)水质与流速参数的空间分配

根据拟使用的河流水质模拟模型，采用结构性或非结构性网格将分析河段离散化。根据离散后的空间信息，将水质与流速参数分配到各网格中。

2.建立溯源数学模型

(1)地质统计学溯源模型

假定河流中污染物未发生化学反应或发生一级化学反应，将污染物排放过程与河流中污染物浓度观测值间的关系概化为式1：

z＝Hs+v 式1

式中：z为一个由观测数据组成的m×1维向量；s为待求未知函数经离散后形成的一个n×1维向量，代表污染物排放过程；v为测量误差组成的m×1维向量；H为由传递函数构成的灵敏度矩阵，代表了污染物在河水中的迁移过程。

假定v的均值为0，其协方差矩阵为R，s为一个随机向量，其数学期望和协方差分别由式2和式3表示。

E[s]＝Xβ 式2

Q(θ)＝E[(s-Xβ)(s-Xβ)

式中：X为一个已知的n×p的矩阵；β为p个漂移系数组成的向量；Q(θ)为关于未知参数θ的已知函数，通常取高斯分布函数。

(2)污染物浓度非负性约束

考虑污染物浓度的非负性，对未知函数s进行如下变换：

将式4代入式1得到

实施上述变换后，

3.计算传递函数

为便于说明问题，此处用一维问题说明传递函数的计算方法。由于假定河流中污染物未发生化学反应或发生一级化学反应，污染物迁移过程可用线性偏微分方程表示。

式中c为污染物浓度，u为流速，D为弥散系数，t为时间，x为迁移距离，r(c)为反应项。

此时，污染物浓度可由传递函数方程表示：

其中s(t)为实际排放过程，f(·)为传递函数。在迁移问题较复杂的情况下，难以得到传递函数的解析形式。因此，本发明通过差分方法计算传递函数，从而得到含有一级化学反应过程的传递函数。

首先，通过变量替换将式8变为式9：

当式9中t>0时，s(T-t)＝1，t≤0时，s(T-t)＝0，得到穿透曲线方程：

使式10对时间t求导，即可得到传递函数在坐标x时刻t的值。

因此，可先根据假定的污染源位置，利用正向的河流污染物迁移模型计算s(t)＝1条件下各观测点的浓度曲线，然后采用向后差分方法(式12)计算各观测点在各个时刻的传递函数值，从而得到由传递函数所组成的灵敏度矩阵H，最后利用地质统计学方法反演出相应污染源污染物的排放过程。

4.求解溯源数学模型

待求未知函数离散形式s的求解过程可分成两个步骤，首先是确定未知参数θ，而后估计未知函数s。当式13表示的概率取最大值时，θ的取值即为待求的参数值。

Σ＝HQH

Ξ＝Σ

式13描述的极值问题，等同于求式16的最小值。

L(θ)对θ的导数取0值时，L(θ)取最小值。可采用Gauss-Newton迭代法、最速下降法等方法求解此极值问题，也可以结合遗传算法等智能优化方法求解此问题。当迭代计算收敛后，再利用θ求Q，并代入式17，解得到n×m的系数矩阵Λ和p×n的拉格朗日乘子阵M。

对未知函数s的最佳估计

V＝-XM+Q-QH

此时，得到假定污染源对应的排放过程。

5.确定污染源的位置与排放过程

将假定污染源与排放过程输入正向模型中，计算各离散网格的水质与流速，根据计算结果与测量结果的拟合效果确定污染源的位置与排放过程。本发明采用下式描述拟合效果：

式中

6.并行加速

采用地质统计学溯源时，求解式16的最小化问题需要进行大量迭代计算。在溯源问题较为复杂时，由于计算负荷较大，需要的计算时间较长。为提高计算效率，采用GPU并行加速方法改进原最小化问题的计算方法。

以最速下降法为例，该方法法要求每次迭代后的目标函数值应是下降的，即要求L(θ

针对这一特点，本文将GPU中Grid(网格)中的几个Block(线程块)作为一个Block组，每个Block组执行一个λ

由于遗传算法等智能优化算法也可以采用GPU并行的方式来提高计算效率。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的装置而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。

对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。