污染物在水体环境下的浓度扩散，是一个非常复杂的过程，需要综合考虑水流速度以及扩散系数变化因素。近年来，国内外很多学者对污染物扩散浓度场进行了数值模拟^［

式（1）中$x$为自排污口沿河流流向的纵向坐标；$y$为垂直于$x$轴从排污口指向河心的横向坐标，$u_{\mathrm{1}}$表示离岸横坐标$y_{\mathrm{1}}$的流速，$\dot{M}$表示岸边污染源单位时间的排放质量，$H$表示河道平均水深，$m$和$n$为正常数指数，$\gamma$和$\alpha$为正常数，$p=\mathrm{2}+m-n>\mathrm{1}$，$\phi =(\mathrm{1}+m)/p$，$\mathrm{\Gamma }(\phi )$为完全伽马函数。

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

图2  浓度场模拟

Fig.2  Concentration field simulation

下载:  原图 | 高精图 | 低精图

表1  水污染浓度场数值模拟参数

Table1  Parameters for concentration field simulation

参数	意义	值	单位
$\dot{M}$	岸边稳定点源单位时间排放质量	1 000	kg/s
$H$	平均水深	5	m
$u_1$	离岸横坐标为$y_1$的流速	0.1	m/s
$y_1$	离岸横坐标	0~80	m
$x_{\mathrm{m}\mathrm{a}\mathrm{x}}$	最大长度	150	m
$y_{\mathrm{m}\mathrm{a}\mathrm{x}}$	最大宽度	80	m
$m$	横向分布指数	0
$n$	横向分布指数	0
$p$	分类条件	2
$\gamma$	变横向扩散系数中的正常数系数	5.96
$\alpha$	变横向扩散系数中的正常数指数	0.01
$\phi$	正常数指数	0.5

下载:  CSV 下载:  表格图片

从图2（b）可以看出，在近排放点位置的地方，浓度值比较大，等浓度线也比较密集，随着位置逐渐变化，浓度场变化符合文献［

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

假设$x$是浓度场的任一点，$Z(x)$是该点的测量值，如果在整个浓度场内总共有$n$个观测点$x_{\mathrm{1}}$，$x_{\mathrm{2}}$，$\cdots$，$x_n$，那么对于任意待估计点的实测浓度值$Z_v\left(x\right)$，其估计的浓度值$Z_v^{*}\left(x\right)$，可以通过该待估点影响范围内的$n$个有效样本值$Z_v\left(x_i\right)(i=\mathrm{1,2},\cdots ,n)$的线性组合来表示，即

其中，权重系数${\lambda }_i$是各已知样本$Z(x_i)$在估计$Z_v^{*}\left(x\right)$时影响大小的系数，而估计$Z_v^{*}\left(x\right)$的好坏主要取决于怎样计算或选择权重系数。在求取权重系数时必须满足两个条件：

（i） 估计是无偏的，即偏差的数学期望为零

（ii） 估计是最优的，即使估计值$Z_v^{*}\left(x\right)$和实际值$Z_v\left(x\right)$之差的平方和最小

为了使得估计方差最小，可以通过拉格朗日乘数原理进行计算，令

对式（5）求对拉格朗日乘子$\mu$以及${\lambda }_i$偏导，令其偏导数为0，可以得到方程组

式中$c\left(x_i,x_j\right)$表示协方差，求解方程组（6）即可得到${\lambda }_i$，代入式（2）就可得到待估计点$x$的估计值。当变异函数存在时，协方差函数与变异函数存在如下关系

因此，式（6）也可以写成

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

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

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

图3  不同样本数拟合球形变差函数的浓度场重构

Fig.3  Reconstruction of concentration field by fitting spherical variation function with different sample numbers

下载:  原图 | 高精图 | 低精图

表2  拟合球形变差函数结果

Table 2  Fitness statistics of spherical variation function models

随机样本数	块金常数	基台值	变程
100	0	1 677.2	74.087 2
200	0	1 406.1	55.7
500	0	973.8	48.9

下载:  CSV 下载:  表格图片

如图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自由度动态模型如下式所示^［

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

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

其中$\mathit{x}$表示状态向量，$\mathit{u}$表示控制输入向量，$\widehat{X}$ 和$\widehat{U}$和表示状态和控制输入约束集，约束状态和输入的最大值与最小值，那么模型预测控制问题可以描述为如式（12）所示的带约束条件下的最优化问题。

式中$\mathit{Q}$和$\mathit{P}$是一个半正定矩阵，$\mathit{R}$为正定矩阵，${\mathit{r}}_k$表示参考轨迹。为了求最优化问题，可以通过参数向量$\mathit{U}=[{\mathit{u}}_{\mathrm{0}},{\mathit{u}}_{\mathrm{1}},\dots ,{\mathit{u}}_{N-\mathrm{1}}$］将最优化问题转化为多参数规划问题，即通过当前的状态向量值${\mathit{x}}_{\mathrm{0}}$，获得预测的状态向量${\mathit{x}}_{\mathrm{1}},{\mathit{x}}_{\mathrm{2}},\dots ,{\mathit{x}}_{N-\mathrm{1}}$以及控制输入向量${\mathit{u}}_{\mathrm{0}},{\mathit{u}}_{\mathrm{1}},\dots ,{\mathit{u}}_{N-\mathrm{1}}$的表达式，如式（13）所示

将式（13）写为

代价函数（12）也可以写为

式中，$\overline{\mathit{Q}}=\mathrm{b}\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{k}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}([\mathit{P},\mathit{P},\dots ,\mathit{Q}])$，$\overline{\mathit{R}}=\mathrm{b}\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{k}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}([\mathit{R},\dots ,\mathit{R}])$，$\mathit{r}$为参考轨迹。将式（14）代入式（15）得

其中$\mathit{H}={\mathit{H}}_{\mathit{u}}^{\mathrm{T}}\overline{\mathit{Q}}{\mathit{H}}_{\mathit{u}}+\overline{\mathit{R}}$，$\mathit{F}=\mathrm{2}\left({\mathit{x}}_{\mathrm{0}}^{\mathrm{T}}{\mathit{H}}_{\mathit{x}}^{\mathrm{T}}-\mathit{r}\right)\overline{\mathit{Q}}{\mathit{H}}_{\mathit{u}}$，$\mathit{M}\left({\mathit{H}}_{\mathit{u}},{\mathit{H}}_{\mathit{x}},{\mathit{x}}_{\mathrm{0}},\mathit{r}\right)$是一个常数，因此，最优化问题（15）是式（16）表示的关于$U$的二次函数，最优化问题可以转化为二次规划问题。

根据式（10）、式（15）和式（16），使用移动机器人结合空间插值的浓度边界追踪的目标函数如式（17）所示。

式（17）中，$N_p$和$N_u$分别表示预测步长与控制步长。$\mathrm{𝒦}\left(\mathit{x}(k+i),\mathit{y}(k+i)\right)$表示根据移动机器人预测所在位置，结合拟合参数进行Kriging插值所得到的预测位置的污染浓度值，$c_r$表示追踪的参考浓度值，$∆\omega$为控制量的变化，$\mathit{Q}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}[q_{\mathrm{1}},q_{\mathrm{2}},\dots ,q_{N_p}]$和$\mathit{R}$表示惩罚权重。通过最小化目标函数$J$，就可以得到控制量$\omega$。

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

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

图5  浓度边界追踪

Fig.5  Boundary tracking for water pollution concentration field

下载:  原图 | 高精图 | 低精图