结合空间插值的水污染浓度边界跟踪*

张谦，裴海龙，邬依林

1. 广东第二师范学院计算机科学系，广东 广州 510303

2. 华南理工大学/自主系统与网络控制教育部重点实验室/广东省无人机系统工程技术研究中心，广东 广州 510640

水体污染是水资源保护的主要问题之一，在常规水环境保护上，各国在废水排放、水质监测以及质量标准等方面均已建立其行之有效的水环境监控保护体系，可以有效保护水资源［1］。然而随着经济社会的发展，工业排放、水上作业、农业生产等活动，均有可能造成一些不可预见的突发性水污染事故，进而对污染监控提出了更高的要求，需要快速、准确以及廉价的水体污染监控方法［2］。一般而言，对于水体污染监测，可以采用卫星遥感［2-4］、传感器网络、数值模拟等方式进行监测。然而，对于突发性水污染事件，一般具有局部性的特点，而针对突发性水污染扩散的数值模拟，是水污染事故应急处理中的重要手段［5］。考虑到在不同水体环境下水污染浓度场的数值模拟条件，文献［6］模拟了急流条件下突发水污染的污染物二维水流-输运模型；文献［7］结合深度平均流速的横向指数分布和横向扩散系数的二维变化关系式，提出岸边排放污染混合区的形状特征与分类；文献［8］提出了倾斜岸水面污染源下角形域污染物二维扩散浓度分布；这些模型为水污染控制与监测以及应急处理中的快速决策提供了很好的依据。但在实际的水污染突发事件的应急处置中，更关注的是污染物扩散的边界、浓度分布及其运动过程［8］。不同于常规的水环境监测，水利部门可以通过固定观测点构建监测网络，然而对于突发性水污染区域的监测，因为污染源位置的不确定性，更方便的是使用移动设备对污染区域进行监测，建立污染区域的浓度场分布，进而确定浓度边界。然而，考虑到地形等因素的影响，很难使用移动监测设备对突发性水污染扩散浓度场进行完全的实时测量，因此依靠移动采样设备获取测量样本，采用空间插值重建浓度场是一种有效的手段。地统计空间插值方可以估计任何地点的数据，其中克里金（Kriging）插值法是最常用的空间插值法［9-11］。如果能够通过移动检测设备采样并使用克里金插值法，实现水污染扩散区域的浓度场重构，那么将为移动机器人进行水污染浓度场边界的确立提供便利。浓度场边界追踪属于机器人轨迹跟踪问题，在这方面很多学者使用模型预测控制方法开展研究［12-16］，例如文献［12］使用模型预测控制方法设计了一个水下机器人轨迹跟踪系统；文献［13］考虑实际工作约束和限制，提出了一种基于模型预测控制轮式机器人轨迹规划方法；文献［14］针对水下动态目标跟踪，提出了基于模型预测控制的水下机器人跟踪控制方法等。然而，这些跟踪的方法，一般都是由状态方程直接收敛到目标轨线。在水污染浓度边界跟踪上，因为无法事先获取目标点的位置，必须依靠移动机器人在移动过程中的采样，来预测目标点的位置，同时决定移动机器人的移动方向，达到跟踪的目的。

因此，针对上述问题，提出了使用模型预测控制与空间插值算法相结合，进行水污染扩散浓度边界跟踪的方法，通过空间插值重构的水污染扩散浓度场分布结合系统的状态模型，再加入符合设定环境的控制量约束，寻找局部最优点作为移动机器人跟踪的浓度场边界最优轨迹，从而确立污染区域浓度场边界。

1 浓度边界跟踪方案

如图1 所示框架，为结合空间插值与模型预测控制的浓度场边界跟踪方案。为了验证方案的可行性，首先对水污染扩散区域的浓度场进行了数值模拟，并将数值模拟的结果作为假设真实发生水污染的浓度场，然后通过对真实浓度场的随机采样，对比不同条件下的空间插值结果，选择合适的参数拟合结果进行污染区域浓度场的重构。最后根据重构的浓度场以及拟合参数，结合水面移动机器人状态模型，对移动机器人的移动位置浓度值进行预测，并以模型预测控制方法对移动机器人所在区域进行局部寻优，寻找并确定下一步的移动目标，以此迭代实现对浓度场边界的跟踪并最终确立污染扩散区域的浓度边界。

图1 浓度边界跟踪方案Fig.1 Scheme of boundary tracking for water pollution concentration field

2 污染区域浓度场分布与重构

2.1 污染源扩散模型

污染物在水体环境下的浓度扩散，是一个非常复杂的过程，需要综合考虑水流速度以及扩散系数变化因素。近年来，国内外很多学者对污染物扩散浓度场进行了数值模拟［6-8］，其中文献［7］提出考虑河流流速和横向扩散系数变化的河流岸边排放污染混合区的估算方法，描述了河流岸边排放污染区的形状特征，较好地表征了排污等值线形状，对水污染排放区域的数值描述有较强的代表性。本文主要研究的是对污染源在水面扩散时的浓度场边界追踪问题，浓度场数值模型［7］为

式（1）中x为自排污口沿河流流向的纵向坐标；y为垂直于x轴从排污口指向河心的横向坐标，u1表示离岸横坐标y1的流速，Ṁ表示岸边污染源单位时间的排放质量，H表示河道平均水深，m和n为正常数指数，γ和α为正常数，p= 2 +m-n＞1，φ=(1+m)/p，Γ(φ)为完全伽马函数。

如图2 所示为简化的模拟岸边污染源扩散浓度场情况，图2（a）为浓度场扩散三维视图，图2（b）为浓度场等高线，仿真参数如表1 所示，这里假设单位时间排放量Ṁ = 1000 kg/s，平均水深H= 5 m，河流流速为常数0.1 m/s，横向分布指数m和n都为0，p= 2+m-n，γ= 5.96，α= 0.01，φ=(1+m)/p，Γ(φ)= 1.

图2 浓度场模拟Fig.2 Concentration field simulation

表1 水污染浓度场数值模拟参数Table1 Parameters for concentration field simulation

从图2（b）可以看出，在近排放点位置的地方，浓度值比较大，等浓度线也比较密集，随着位置逐渐变化，浓度场变化符合文献［7］中所归纳的岸边排放污染源混合区标准型类别，其形状特征近似于半椭圆形。

2.2 浓度场重建

通过式（1）可以模拟出污染源扩散的浓度场情况，但在实际应用过程中，采用的是单一移动采样设备对污染源的浓度场边界进行追踪，也很难对浓度场分布进行建模后得到精确的分布情况后再进行追踪，因此，更直接的方式是通过采样设备对浓度场分布范围内进行随机浓度采样，然后通过地统计学方式，对采样得到的浓度值进行空间插值，重构污染区域的浓度场。在地统计学中，最常用的方法是Kriging 插值法。Kriging 方法从变量相关性和变异性出发，在有限区域内对区域化变量的取值进行无偏、最优估计的一种方法。很多文献对Kriging插值法进行了详细阐述［9-11］，其基本算法思路如下所述。

假设x是浓度场的任一点，Z(x)是该点的测量值，如果在整个浓度场内总共有n个观测点x1，x2，…，xn，那么对于任意待估计点的实测浓度值Zv(x)，其估计的浓度值Z*v(x)，可以通过该待估点影响范围内的n个有效样本值Zv(xi)(i= 1，2，…，n)的线性组合来表示，即

为了拟合观测数据值与空间距离之间的关系，一般变异函数有多种模型，常见的有高斯模型、指数模型、球状模型等。例如，高斯模型可以表示为

综上所述，在进行污染源扩散的浓度场重建的时候，就是通过采样得到的观测值，进行变异函数拟合，再求解式（8）得到权重系数λi，即可对不同位置点进行估计，从而最终重建出整个浓度场。

为了验证使用Kriging 插值重构岸边污染源扩散浓度场的效果，图3 以不同的随机样本数拟合球形变差函数模型，对图2模拟的浓度场进行了插值测仿真，以不同的采样样本数，拟合球形变差函数的结果如表2所示，仿真结果如图3所示。

表2 拟合球形变差函数结果Table 2 Fitness statistics of spherical variation function models

图3 不同样本数拟合球形变差函数的浓度场重构Fig.3 Reconstruction of concentration field by fitting spherical variation function with different sample numbers

如图3 所示，图3（a）使用100 个样本、图3（b）使用200 个样本、图3（c）使用500 个样本，分别对图2所示模拟的浓度场进行Kriging 插值所得的浓度场预测等高线、变差函数拟合效果以及预测结果与实际分布之间的误差，从表2 以及图3 可知，样本数越多，变差函数的拟合效果越好，对浓度场的重建效果越好，此外，预测误差分布主要集中在离污染源排放位置较近的区域，越远离排放源，误差越接近于0。因为本文主要是对污染源扩散边界进行跟踪，跟踪的轨线离排放源较远，因此经过插值重构的污染源浓度场边界，可以反映真实的浓度扩散，是可以用来作为追踪轨线的。

图4为采用500个采样点时，变差函数分别拟合高斯模型和指数模型时重建浓度场分布的效果。从图4（a）拟合高斯模型的效果来看，重构的浓度场边界不能完成呈现半椭圆形状，且重建的浓度场与原浓度场之间的误差有较大的波动，尤其是离排放源较远的地方，误差变动仍然较大，不利于后续开展边界追踪操作。图4（b）为拟合指数模型效果，从重构的浓度场分布情况看，等浓度线重构效果优于图3（b）采用200 采样点的球状模型，但在x方向较小时，y方向曲线效果不如图3（c），误差波动比较大，且变差函数拟合的效果不理想。因此，综上分析，在进行边界追踪重建浓度场分布的时候，将采用500个采样点拟合球状变差函数模型进行空间插值。

图4 变差函数拟合高斯模型和指数模型的浓度场重构Fig.4 Reconstruction of concentration field by fitting Gaussian and exponential variation function

3 浓度边界追踪

3.1 轨迹跟踪

对污染物扩散浓度边界进行追踪，事实上是控制移动机器人的位置和方向，使得移动机器人自起始位置开始，其移动轨迹快速靠近浓度边界等浓度线，并沿着等浓度线进行移动。假设移动机器人是一个在水面移动的采样设备，那么它的3自由度动态模型如下式所示［13-14］

式（10）中，状态向量[x，y，θ]T分别表示位置和移动的方向角，v表示移动的速度，可以设为常数。移动机器人在进行浓度边界追踪时，移动机器人首先要根据随机采样得到样本进行当前位置的Kriging 插值浓度预测，根据预测值与重构的边界浓度值进行比较，寻找局部最优点，从而控制移动机器人的移动位置和方向，如此循环操作达到边界最终的目的，这个过程与模型预测控制进行移动机器人轨迹跟踪相匹配。因此，将以模型预测控制的方式对浓度边界进行追踪。

在模型预测控制中，对于连续的系统状态方程，通过选取合适的采样频率，将系统离散化为

式（17）中，Np和Nu分别表示预测步长与控制步长。K(x(k+i)，y(k+i))表示根据移动机器人预测所在位置，结合拟合参数进行Kriging 插值所得到的预测位置的污染浓度值，cr表示追踪的参考浓度值，Δω为控制量的变化，Q= diag[q1，q2，…，qNp]和R表示惩罚权重。通过最小化目标函数J，就可以得到控制量ω。

3.2 边界追踪仿真

根据式（10）和式（17），结合第二节所述浓度场重构，水面移动机器人跟踪浓度边界进行跟踪仿真。由式（17），假设参考轨迹为污染源浓度场分布中浓度值为50 的等浓度线，轨迹预测值为当前位置下，通过500个采样点拟合球形变差函数后的Kriging插值预测值，预测步长和控制步长均设置为3，系统的控制量u=ω，移动速度设置为常数，仅对角速度变化进行约束，约束条件为|ω|＜π/2，也就是角度θ=ωt不能在单位时间内进行大方向的转动。

如图5 所示，图5（a）所模拟的真实的浓度场分布数值设置参数与图2（b）一致，图中黑色的点为在浓度场范围内随机采样点的位置。图5（b）是变差函数拟合球状模型的结果，从图上可知，基本能够完全拟合。图5（c）是通过采样模拟的真实浓度场，进行空间插值而重构的浓度场，这个重构的浓度场就是将要进行轨迹跟踪的边界。图5（d）即为浓度边界的轨迹跟踪结果，红色点表示移动机器人从起始位置向指定浓度值所在区域靠拢时的每一次计算所得到的移动轨迹，红色的星形图标表示移动机器人的位置移动到了预定浓度边界，并将从该位置开始进行等浓度线边界追踪，当到达岸边时跟踪停止。图中假设移动机器人起始位置为任意设置，并将跟踪浓度值为50 的等浓度线，移动机器人的速度v= 1m/s 设置为常数。对于移动机器人的每一次迭代位置，均需要通过Kriging 插值计算预测步长决定的移动机器人所在位置的浓度值，并根据预测的浓度值与等浓度线的差，以及控制量增量，对最小化式（17）的目标函数求最优值，从而得到所需的最优控制量ω，以此控制量修正移动机器人的位置，当移动到下一个位置时，再重复上述预测求最优的步骤，从而最终实现水面移动机器人跟踪污染物扩散浓度边界。从图5（d）可以看到，浓度边界跟踪的效果非常好，能沿着等浓度线进行移动，这也验证了通过空间插值预测的模型预测控制浓度边界轨迹跟踪方法可行。

图5 浓度边界追踪Fig.5 Boundary tracking for water pollution concentration field

4 结 论

本文通过数值模拟河流岸边排放污染区浓度扩散情况，并比较不同的Kriging 空间插值变差函数拟合方法，重构污染区域的浓度场分布，通过模型预测控制方式，实现使用水面移动机器人设备跟踪浓度扩散边界轨迹的目的。仿真实验结果表明，该方法能够稳定、快速、准确地进行浓度边界跟踪，可以为实际应用过程中的区域跟踪提供有价值的借鉴。