CN103971005A - 一种靠近湖泊岸边的污染源定位方法 - Google Patents

Abstract

本发明涉及一种靠近湖泊岸边的污染源定位方法。其技术方案是：将t_k时刻的N个传感器节点测量的水体离子浓度时刻的待求向量X的后验估计和t_k-1时刻的误差协方差矩阵用动态数据处理模块进行处理，得到t_k时刻的待求向量X的后验估计和t_k时刻的误差协方差矩阵k：＝k+1，再调用动态数据处理模块，直至计算出k＝L时的t_L时刻的待求向量X的后验估计则为所求的污染源的横坐标，

Description

一种靠近湖泊岸边的污染源定位方法

技术领域

本发明属于污染源定位技术领域。具体涉及一种靠近湖泊岸边的污染源定位方法。

背景技术

湖泊是工农业用水、居民饮用水的重要水源。近年来，水污染事件的频发已引起人类的广泛关注，湖泊水体的污染给社会、经济和人们的生活带来了不可估量的损失。由于湖泊水体的污染主要来自工业废水、废渣和生活垃圾的沿岸排放，因此，靠近湖泊岸边污染源的定位对保护水资源具有重要意义。

现有靠近湖泊岸边的污染源定位方法有：近似函数法和非线性最小二乘法。近似函数法是在已知污染源扩散模型的基础上，用分段函数近似污染源扩散模型中的非线性函数，得到分段扩散模型。分段扩散模型将扩散过程简化为不稳定阶段和稳定阶段。采用该方法时需提前给定浓度阈值以判定扩散处于哪个阶段，然后选择处于同一阶段的两个测量值计算出污染源的位置。然而在实际应用中该浓度阈值难以预先确定，该方法不易实现。并且该方法忽略了过渡阶段的浓度信息，影响了定位精度。该方法的定位效果严重依赖于被选取的少量传感器节点的测量值，所以此方法易受测量噪声的干扰，定位结果的稳定性差。

非线性最小二乘法是在实际的约束条件下，满足测量的浓度值与理论的浓度值之差最小时的参数估计。该方法在每个测量时刻都有历史时刻的测量值和当前时刻的测量值参与最小二乘数值运算，然而随着测量次数的增加，数据量会急剧增长，在不同的测量时刻有重复而繁杂的最小二乘数值计算，因此该方法数据存储量大、工作效率低和耗时长。

发明内容

本发明旨在克服现有技术中的不足，目的是提供一种易实现、定位精度高、工作效率高和耗时短的靠近湖泊岸边的污染源定位方法。

为实现上述目的，本发明采用的技术方案的具体步骤是：

步骤1、以岸边为纵轴，以垂直于岸边为横轴，在水平面上建立直角坐标系。设污染源在(α,β)处从τ时刻连续释放污染物，污染源的质量流率为Q，以污染源的横坐标α、污染源的纵坐标β、污染源的初始扩散时间τ和污染源的质量流率Q为未知量组成待求向量X，待求向量X＝(α,β,τ,Q)^T；传感器节点的总个数为N，传感器节点共测量了L个时刻的水体离子浓度，第i(i＝1,...,N)个传感器节点(x_i,y_i)在t_k(k＝1,...,L)时刻测量的水体离子浓度为：

$\tilde{Z}(x_i,y_i{t}_k)=Z(x_i,y_i,t_{k}X)+V_i^{\left(k\right)}---\left(1\right)$

式(1)中：是第i个传感器节点(x_i,y_i)在t_k时刻的测量噪声，kg/m³；

的均值为0，的方差为D_i，(kg/m³)²；

Z(x_i,y_i,t_k,X)表示t_k时刻第i个传感器节点(x_i,y_i)处的理论水体离子浓度，kg/m³；

$Z(x_i,y_i,t_k,X)=\frac{\mathrm{\ξQ}}{2f\sqrt{\mathrm{\πd}}}[\frac{1}{R_i}\mathrm{erfc}\left(\frac{R_i}{2\sqrt{d(t_k-\mathrm{\τ})}}\right)+\frac{1}{{\stackrel{\‾}{R}}_i}\mathrm{erfc}\left(\frac{{\stackrel{\‾}{R}}_i}{2\sqrt{d(t_k-\mathrm{\τ})}}\right)]---\left(2\right)$

式(2)中：ξ为补量纲常量，tunit为时间单位，h；

R_i为第i个传感器节点(x_i,y_i)到污染源(α,β)的距离，m；

$R_i=\sqrt{{(x_i-\mathrm{\α})}^2+{(y_i-\mathrm{\β})}^2}---\left(3\right)$

为第i个传感器节点(x_i,y_i)到镜像污染源(-α,β)的距离，m；

${\stackrel{\‾}{R}}_i=\sqrt{{(x_i+\mathrm{\α})}^2+{(y_i-\mathrm{\β})}^2}---\left(4\right)$

erfc(x)＝1-erf(x) (5)

f为水体平均深度，m；

d为同向扩散系数，m²/h。

步骤2、设t₀时刻的待求向量X的后验估计为设t₀时刻的误差协方差矩阵为 ${\hat{P}}_0=\mathrm{diag}({\widehat{\mathrm{\σ}}}_{\mathrm{\α}\left(0\right)}^2,{\widehat{\mathrm{\σ}}}_{\mathrm{\β}\left(0\right)}^2,{\widehat{\mathrm{\σ}}}_{\mathrm{\τ}\left(0\right)}^2,{\widehat{\mathrm{\σ}}}_{Q\left(0\right)}^2),$ 其中： 依次为t₀时刻的污染源的横坐标的估计方差、t₀时刻的污染源的纵坐标的估计方差、t₀时刻的污染源的初始扩散时间的估计方差和t₀时刻的污染源的质量流率的估计方差；令k＝1。

步骤3、将t_k时刻的N个传感器节点测量的水体离子浓度时刻的待求向量X的后验估计和t_k-1时刻的误差协方差矩阵用动态数据处理模块进行处理，得到t_k时刻的待求向量X的后验估计和t_k时刻的误差协方差矩阵

步骤4、k：＝k+1，跳转到步骤3，直至计算出k＝L时的t_L时刻的待求向量X的后验估计 ${\hat{X}}_L={({\widehat{\mathrm{\α}}}_{\left(L\right)},{\widehat{\mathrm{\β}}}_{\left(L\right)},{\widehat{\mathrm{\τ}}}_{\left(L\right)},{\hat{Q}}_{\left(L\right)})}^T$ 和误差协方差矩阵

步骤5、在t_L时刻的待求向量X的后验估计中，为所求的污染源的横坐标，为所求的污染源的纵坐标。

2、根据权利要求1所述的靠近湖泊岸边的污染源定位方法，其特征在于所述的用动态数据处理模块进行处理的步骤是：

步骤(1)、由t_k-1时刻的待求向量X的后验估计和t_k-1时刻的误差协方差矩阵得到t_k时刻的第1个Sigma点第p+1个Sigma点和第q+1个Sigma点 ${\mathrm{\λ}}_{k|k-1}^q(q=n+1,...,2n):$

${\mathrm{\λ}}_{k|k-1}^0={\hat{X}}_{k-1}---\left(6\right)$

${\mathrm{\λ}}_{k|k-1}^p={\hat{X}}_{k-1}+{\left[{\left(\sqrt{(n+1){\hat{P}}_{k-1}}\right)}_p\right]}^T---\left(7\right)$

${\mathrm{\λ}}_{k|k-1}^q={\hat{X}}_{k-1}-{\left[{\left(\sqrt{(n+1){\hat{P}}_{k-1}}\right)}_{(q-n)}\right]}^T---\left(8\right)$

式(7)中：表示开方矩阵的第p行；

式(8)中：表示开方矩阵的第(q-n)行。

计算t_k时刻第1个Sigma点第p+1个Sigma点和第q+1个Sigma点依次对应的权值W₀、W_p和W_q：

$W_0=\frac{1}{n+1}---\left(9\right)$

$W_p=\frac{1}{2(n+1)}---\left(10\right)$

$W_q=\frac{1}{2(n+1)}---\left(11\right)$

式(7)(8)(9)(10)(11)中：n为待求向量X的维数，n＝4。

步骤(2)、预测t_k时刻各传感器节点(x_i,y_i)(i＝1,...,N)处的水体离子浓度

${\left(Z_k^{-}\right)}_i=W_{0}Z(x_i,y_i,t_k,{\mathrm{\λ}}_{k|k-1}^0)+\underset{p=1}{\overset{n}{\mathrm{\Σ}}}{W}_{p}Z(x_i,y_i,t_k,{\mathrm{\λ}}_{k|k-1}^p)+\underset{q=n+1}{\overset{2n}{\mathrm{\Σ}}}{W}_{q}Z(x_i,y_i,t_k,{\mathrm{\λ}}_{k|k-1}^q)---\left(12\right)$

式(12)：W₀、W_p和W_q依次同式(9)、式(10)和式(11)；

和依次同式(6)、式(7)和式(8)；

n为待求向量X的维数，n＝4；

表示当待求向量X为t_k时刻的第1个Sigma点时，t_k时刻第i个传感器节点(x_i,y_i)处的理论水体离子浓度；

表示当待求向量X为t_k时刻的第p+1个Sigma点时，t_k时刻第i个传感器节点(x_i,y_i)处的理论水体离子浓度；

表示当待求向量X为t_k时刻的第q+1个Sigma点时，t_k时刻第i个传感器节点(x_i,y_i)处的理论水体离子浓度；

步骤(3)、计算t_k时刻的待求向量X的后验估计和t_k时刻的误差协方差矩阵首先计算t_k时刻的待求向量X的后验估计

${\hat{X}}_k={\hat{X}}_k^{-}+G_k\left(\begin{array}{c}\tilde{Z}(x_1,y_1,t_k)-{\left(Z_k^{-}\right)}_1\\ \tilde{Z}(x_2,y_2,t_k)-{\left(Z_k^{-}\right)}_2\\ \·\\ \·\\ \·\\ \tilde{Z}(x_N,y_N,t_k)-{\left(Z_k^{-}\right)}_N\end{array}\right)---\left(13\right)$

式(13)中：为t_k时刻的待求向量X的先验估计，

${\hat{X}}_k^{-}={\hat{X}}_{k-1}---\left(14\right)$

G_k是t_k时刻的增益矩阵，G_k用于修正t_k时刻的待求向量X的先验估计得到t_k时刻的待求向量X的后验估计

$G_k=P_{\mathrm{XZ}}{S}_k^{-1}---\left(15\right)$

和依次表示第1个传感器节点、第2个传感器节点和第N个传感器节点在t_k时刻测量的水体离子浓度；

和依次表示预测的t_k时刻第1个传感器节点(x₁,y₁)处、第2个传感器节点(x₂,y₂)处和第N个传感器节点(x_N,y_N)处的水体离子浓度；

式(14)中：为t_k-1时刻的待求向量X的后验估计。

式(15)中：

$P_{\mathrm{XZ}}=W_0({\mathrm{\λ}}_{k|k-1}^0-{\hat{X}}_k^{-}){\left(\mathrm{\Δ}{Z}_0\right)}^T+\underset{p=1}{\overset{n}{\mathrm{\Σ}}}{W}_p({\mathrm{\λ}}_{k|k-1}^p-{\hat{X}}_k^{-}){\left(\mathrm{\Δ}{Z}_p\right)}^T+\underset{q=n+1}{\overset{2n}{\mathrm{\Σ}}}{W}_q({\mathrm{\λ}}_{k|k-1}^q-{\hat{X}}_k^{-}){\left(\mathrm{\Δ}{Z}_q\right)}^T---\left(16\right)$

S_k＝R_k+P_ZZ (17)

式(16)：W₀、W_p和W_q依次同式(9)、式(10)和式(11)；

和依次同式(6)、式(7)和式(8)；

为t_k时刻的待求向量X的先验估计，同式(14)；

n为待求向量X的维数，n＝4；

$\mathrm{\Δ}{Z}_0=\left(\begin{array}{c}Z(x_1,y_1,t_k,{\mathrm{\λ}}_{k|k-1}^0)-{\left(Z_k^{-}\right)}_1\\ Z(x_2,y_2,t_k,{\mathrm{\λ}}_{k|k-1}^0)-{\left(Z_k^{-}\right)}_2\\ \·\\ \·\\ \·\\ Z(x_N,y_N,t_k,{\mathrm{\λ}}_{k|k-1}^0)-{\left(Z_k^{-}\right)}_N\end{array}\right)---\left(18\right)$

$\mathrm{\Δ}{Z}_p=\left(\begin{array}{c}Z(x_1,y_1,t_k,{\mathrm{\λ}}_{k|k-1}^p)-{\left(Z_k^{-}\right)}_1\\ Z(x_2,y_2,t_k,{\mathrm{\λ}}_{k|k-1}^p)-{\left(Z_k^{-}\right)}_2\\ \·\\ \·\\ \·\\ Z(x_N,y_N,t_k,{\mathrm{\λ}}_{k|k-1}^p)-{\left(Z_k^{-}\right)}_N\end{array}\right)---\left(19\right)$

$\mathrm{\Δ}{Z}_q=\left(\begin{array}{c}Z(x_1,y_1,t_k,{\mathrm{\λ}}_{k|k-1}^q)-{\left(Z_k^{-}\right)}_1\\ Z(x_2,y_2,t_k,{\mathrm{\λ}}_{k|k-1}^q)-{\left(Z_k^{-}\right)}_2\\ \·\\ \·\\ \·\\ Z(x_N,y_N,t_k,{\mathrm{\λ}}_{k|k-1}^q)-{\left(Z_k^{-}\right)}_N\end{array}\right)---\left(20\right)$

式(18)中：同式(6)；

表示当待求向量X为t_k时刻的第1个Sigma点时，t_k时刻第1个传感器节点(x₁,y₁)处的理论水体离子浓度；

表示当待求向量X为t_k时刻的第1个Sigma点时，t_k时刻第2个传感器节点(x₂,y₂)处的理论水体离子浓度；

表示当待求向量X为t_k时刻的第1个Sigma点时，t_k时刻第N个传感器节点(x_N,y_N)处的理论水体离子浓度；

和依次表示预测的t_k时刻第1个传感器节点(x₁,y₁)处、第2个传感器节点(x₂,y₂)处和第N个传感器节点(x_N,y_N)处的水体离子浓度。

式(19)中：同式(7)；

表示当待求向量X为t_k时刻的第p+1个Sigma点时，t_k时刻第1个传感器节点(x₁,y₁)处的理论水体离子浓度；

表示当待求向量X为t_k时刻的第p+1个Sigma点时，t_k时刻第2个传感器节点(x₂,y₂)处的理论水体离子浓度；

表示当待求向量X为t_k时刻的第p+1个Sigma点时，t_k时刻第N个传感器节点(x_N,y_N)处的理论水体离子浓度；

和依次表示预测的t_k时刻第1个传感器节点(x₁,y₁)处、第2个传感器节点(x₂,y₂)处和第N个传感器节点(x_N,y_N)处的水体离子浓度。

式(20)中：同式(8)；

表示当待求向量X为t_k时刻的第q+1个Sigma点时，t_k时刻第1个传感器节点(x₁,y₁)处的理论水体离子浓度；

表示当待求向量X为t_k时刻的第q+1个Sigma点时，t_k时刻第2个传感器节点(x₂,y₂)处测量的理论水体离子浓度；

表示当待求向量X为t_k时刻的第q+1个Sigma点时，t_k时刻第N个传感器节点(x_N,y_N)处的理论水体离子浓度；

和依次表示预测的t_k时刻第1个传感器节点(x₁,y₁)处、第2个传感器节点(x₂,y₂)处和第N个传感器节点(x_N,y_N)处的水体离子浓度。

式(17)中：R_k是测量噪声的协方差矩阵，所述测量噪声的协方差矩阵为对角矩阵，

$R_k={\left(\begin{array}{cccc}{D}_1& 0& \·\·\·& 0\\ 0& D_2& \·\·\·& 0\\ \·& \·& \mathrm{}& \·\\ \·& \·& \mathrm{}& \·\\ \·& \·& \mathrm{}& \·\\ 0& 0& \·\·\·& D_N\end{array}\right)}_{N\×N}---\left(21\right)$

$P_{\mathrm{ZZ}}=W_0\left(\mathrm{\Δ}{Z}_0\right){\left(\mathrm{\Δ}{Z}_0\right)}^T+\underset{p=1}{\overset{n}{\mathrm{\Σ}}}{W}_p\left(\mathrm{\Δ}{Z}_p\right){\left(\mathrm{\Δ}{Z}_p\right)}^T+\underset{q=n+1}{\overset{2n}{\mathrm{\Σ}}}{W}_q\left(\mathrm{\Δ}{Z}_q\right){\left(\mathrm{\Δ}{Z}_p\right)}^T---\left(22\right)$

式(22)中：W₀、W_p和W_q依次同式(9)、式(10)和式(11)；

ΔZ₀、ΔZ_p和ΔZ_q依次同式(18)、式(19)和式(20)；

n为待求向量X的维数，n＝4。

然后计算误差协方差矩阵

${\hat{P}}_k=P_k^{-}-G_k{S}_k{\left(G_k\right)}^T---\left(23\right)$

式(23)中： $P_k^{-}={\hat{P}}_{k-1}---\left(24\right)$

G_k和S_k依次同式(15)和式(17)；

式(24)中：为t_k-1时刻的误差协方差矩阵。

由于采用上述技术方案，本发明与现有技术相比具有如下经济效果：

本发明无需提前给定阈值，用每个时刻的测量值去更新待求向量X的先验估计得到待求向量X的后验估计，在实际定位操作中简单易行；其次，本发明中测量的浓度数据可以包含整个扩散过程的信息，避免了因忽略过渡阶段影响定位精度的问题；然后，本发明在每个测量时刻均处理了所有传感器节点(x_i,y_i)(i＝1,...,N)的测量数据，由于在每个测量时刻都考虑到了所有传感器节点(x_i,y_i)(i＝1,...,N)的测量噪声V_i ^(k)(i＝1,...,N)，该方法的定位效果相对稳定。

本发明在定位计算过程中只需保存当前时刻的测量值、上一时刻的待求向量X的后验估计和上一时刻的误差协方差矩阵，用当前时刻的测量值去更新待求向量X的先验估计，计算过程中不需要保存和计算历史测量数据。故本发明具有工作效率高和耗时短的优点。

因此，本发明具有易实现、定位精度高、工作效率高和耗时短的特点。

具体实施方式

下面结合具体实施方式对本发明作进一步的描述，并非对保护范围的限制：

实施例1

一种靠近湖泊岸边的污染源定位方法。该方法的具体步骤是：

步骤1、以岸边为纵轴，以垂直于岸边为横轴，在水平面上建立直角坐标系。设污染源在(α,β)处从τ时刻连续释放污染物，污染源的质量流率为Q，以污染源的横坐标α、污染源的纵坐标β、污染源的初始扩散时间τ和污染源的质量流率Q为未知量组成待求向量X，待求向量X＝(α,β,τ,Q)^T；传感器节点的总个数为8，传感器节点共测量了5个时刻的水体离子浓度，8个传感器节点(x_i,y_i)(i＝1,...,8)在t_k(k＝1,...,5)时刻测量的水体离子浓度见表1.1。

表1.1

第i(i＝1,...,8)个传感器节点(x_i,y_i)在t_k(k＝1,...,5)时刻测量的水体离子浓度为：

$\tilde{Z}(x_i,y_i{t}_k)=Z(x_i,y_i,t_{k}X)+V_i^{\left(k\right)}---\left(1\right)$

式(1)中：是第i个传感器节点(x_i,y_i)在t_k时刻的测量噪声，kg/m³；

的均值为0，的方差为D_i，(kg/m³)²；8个传感器节点(x_i,y_i)(i＝1,...,8)在t_k(k＝1,...,5)时刻的测量噪声的方差见表1.2；

Z(x_i,y_i,t_k,X)表示t_k时刻第i个传感器节点(x_i,y_i)处的理论水体离子浓度，kg/m³；

$Z(x_i,y_i,t_k,X)=\frac{\mathrm{\ξQ}}{2f\sqrt{\mathrm{\πd}}}[\frac{1}{R_i}\mathrm{erfc}\left(\frac{R_i}{2\sqrt{d(t_k-\mathrm{\τ})}}\right)+\frac{1}{{\stackrel{\‾}{R}}_i}\mathrm{erfc}\left(\frac{{\stackrel{\‾}{R}}_i}{2\sqrt{d(t_k-\mathrm{\τ})}}\right)]---\left(2\right)$

式(2)中：ξ为补量纲常量，tunit为时间单位，h；

R_i为第i个传感器节点(x_i,y_i)到污染源(α,β)的距离，m；

$R_i=\sqrt{{(x_i-\mathrm{\α})}^2+{(y_i-\mathrm{\β})}^2}---\left(3\right)$

为第i个传感器节点(x_i,y_i)到镜像污染源(-α,β)的距离，m；

${\stackrel{\‾}{R}}_i=\sqrt{{(x_i+\mathrm{\α})}^2+{(y_i-\mathrm{\β})}^2}---\left(4\right)$

erfc(x)＝1-erf(x) (5)

f为水体平均深度，f＝10m；

d为同向扩散系数，d＝0.5m²/h。

表1.2

步骤2、设t₀＝0h时刻的待求向量X的后验估计为设t₀＝0h时刻的误差协方差矩阵为其中：4、4、4和100依次为t₀＝0h时刻的污染源的横坐标的估计方差、t₀＝0h时刻的污染源的纵坐标的估计方差、t₀＝0h时刻的污染源的初始扩散时间的估计方差和t₀＝0h时刻的污染源的质量流率的估计方差；令k＝1。

步骤3、将t₁＝6h时刻的8个传感器节点测量的水体离子浓度时刻的待求向量X的后验估计和t₀＝0h时刻的误差协方差矩阵用动态数据处理模块进行处理。

步骤3.1、由t₀＝0h时刻的待求向量X的后验估计和t₀＝0h时刻的误差协方差矩阵得到t₁＝6h时刻的第1个Sigma点第p+1个Sigma点和第q+1个Sigma点 ${\mathrm{\λ}}_{1|0}^q(q=5,...,8):$

${\mathrm{\λ}}_{1|0}^0={\hat{X}}_0---\left(6\right)$

${\mathrm{\λ}}_{1|0}^p={\hat{X}}_0+{\left[{\left(\sqrt{(4+1){\hat{P}}_0}\right)}_p\right]}^T---\left(7\right)$

${\mathrm{\λ}}_{1|0}^q={\hat{X}}_0-{\left[{\left(\sqrt{(4+1){\hat{P}}_0}\right)}_{(q-4)}\right]}^T---\left(8\right)$

式(7)中：表示开方矩阵的第p行；

式(8)中：表示开方矩阵的第(q-4)行。

t₁＝6h时刻的9个Sigma点如下：

$\begin{array}{c}{\mathrm{\λ}}_{1|0}^0={\left(\mathrm{2,8,3,80}\right)}^T;{\mathrm{\λ}}_{1|0}^1={\left(\mathrm{6.483302,8,3,80}\right)}^T;{\mathrm{\λ}}_{1|0}^2={\left(\mathrm{2,12.4833,3,80}\right)}^T;{\mathrm{\λ}}_{1|0}^3={\left(\mathrm{2,8,7.483302,80}\right)}^T;\\ {\mathrm{\λ}}_{1|0}^4={\left(\mathrm{2,8,3,102.3629}\right)}^T;{\mathrm{\λ}}_{1|0}^5={(-\mathrm{2.4833,8,3,80})}^T;{\mathrm{\λ}}_{1|0}^6={\left(\mathrm{2,3.516698,3,80}\right)}^T;{\mathrm{\λ}}_{1|0}^7={(\mathrm{2,8},-\mathrm{1.4833,80})}^T;\\ {\mathrm{\λ}}_{1|0}^8={\left(\mathrm{2,8,3,57.63708}\right)}^T.\end{array}$

计算t₁＝6h时刻第1个Sigma点第p+1个Sigma点和第q+1个Sigma点依次对应的权值W₀、W_p和W_q：

$W_0=\frac{1}{4+1}---\left(9\right)$

$W_p=\frac{1}{2(4+1)}---\left(10\right)$

$W_q=\frac{1}{2(4+1)}---\left(11\right)$

t₁＝6h时刻的9个Sigma点依次对应的权值如下：

W₀＝0.2；W₁＝0.1；W₂＝0.1；W₃＝0.1；W₄＝0.1；W₅＝0.1；W₆＝0.1；W₇＝0.1；W₈＝0.1。

步骤3.2、预测t₁＝6h时刻各传感器节点(x_i,y_i)(i＝1,...,8)处的水体离子浓度

${\left(Z_1^{-}\right)}_i=W_{0}Z(x_i,y_i,t_1,{\mathrm{\λ}}_{1|0}^0)+\underset{p=1}{\overset{4}{\mathrm{\Σ}}}{W}_{p}Z(x_i,y_i,t_1,{\mathrm{\λ}}_{1|0}^p)+\underset{q=5}{\overset{8}{\mathrm{\Σ}}}{W}_{q}Z(x_i,y_i,t_1,{\mathrm{\λ}}_{1|0}^q)---\left(12\right)$

式(12)：W₀、W_p和W_q依次同式(9)、式(10)和式(11)；

和依次同式(6)、式(7)和式(8)；

表示当待求向量X为t₁＝6h时刻的第1个Sigma点时，t₁＝6h时刻第i个传感器节点(x_i,y_i)处的理论水体离子浓度：

当待求向量X为t₁＝6h时刻的第1个Sigma点时，t₁＝6h时刻各传感器节点处的理论水体离子浓度如下：

$\begin{array}{c}Z(x_1,y_1,t_1,{\mathrm{\λ}}_{1|0}^1)=0.304573;Z(x_2,y_2,t_1,{\mathrm{\λ}}_{1|0}^0)=0.169127;\\ Z(x_3,y_3,t_1,{\mathrm{\λ}}_{1|0}^0)=0.141195;Z(x_4,y_4,t_1,{\mathrm{\λ}}_{1|0}^0)=0.497168;\\ Z(x_5,y_5,t_1,{\mathrm{\λ}}_{1|0}^0)=0.3217;Z(x_6,y_6,t_1,{\mathrm{\λ}}_{1|0}^0)=0.722341;\\ Z(x_7,y_7,t_1,{\mathrm{\λ}}_{1|0}^0)=0.34804;Z(x_8,y_8,t_1,{\mathrm{\λ}}_{1|0}^0)=0.162404.\end{array}$

表示当待求向量X为t₁＝6h时刻的第p+1个Sigma点时，t₁＝6h时刻第i个传感器节点(x_i,y_i)处的理论水体离子浓度；

当待求向量X为t₁＝6h时刻的第2个Sigma点时，t₁＝6h时刻各传感器节点处的理论水体离子浓度如下：

$\begin{array}{c}Z(x_1,y_1,t_1,{\mathrm{\λ}}_{1|0}^1)=0.003248;Z(x_2,y_2,t_1,{\mathrm{\λ}}_{1|0}^0)=0.002208;\\ Z(x_3,y_3,t_1,{\mathrm{\λ}}_{1|0}^1)=0.001944;Z(x_4,y_4,t_1,{\mathrm{\λ}}_{1|0}^1)=0.006499;\\ Z(x_5,y_5,t_1,{\mathrm{\λ}}_{1|0}^1)=0.014043;Z(x_6,y_6,t_1,{\mathrm{\λ}}_{1|0}^1)=0.025606;\\ Z(x_7,y_7,t_1,{\mathrm{\λ}}_{1|0}^1)=0.008697;Z(x_8,y_8,t_1,{\mathrm{\λ}}_{1|0}^1)=0.005376;\end{array}$

当待求向量X为t₁＝6h时刻的第3个Sigma点时，t₁＝6h时刻各传感器节点处的理论水体离子浓度如下：

$\begin{array}{c}Z(x_1,y_1,t_1,{\mathrm{\λ}}_{1|0}^2)=0.002315;Z(x_2,y_2,t_1,{\mathrm{\λ}}_{1|0}^2)=0.000847;\\ Z(x_3,y_3,t_1,{\mathrm{\λ}}_{1|0}^2)=0.000627;Z(x_4,y_4,t_1,{\mathrm{\λ}}_{1|0}^2)=0.0045;\\ Z(x_5,y_5,t_1,{\mathrm{\λ}}_{1|0}^2)=0.001852;Z(x_6,y_6,t_1,{\mathrm{\λ}}_{1|0}^2)=0.005346;\\ Z(x_7,y_7,t_1,{\mathrm{\λ}}_{1|0}^2)=0.002229;Z(x_8,y_8,t_1,{\mathrm{\λ}}_{1|0}^2)=0.005376.\end{array}$

当待求向量X为t₁＝6h时刻的第4个Sigma点时，t₁＝6h时刻各传感器节点处的理论水体离子浓度如下：

$\begin{array}{c}Z(x_1,y_1,t_1,{\mathrm{\λ}}_{1|0}^3)=0.107885;Z(x_2,y_2,t_1,{\mathrm{\λ}}_{1|0}^3)=0.042568;\\ Z(x_3,y_3,t_1,{\mathrm{\λ}}_{1|0}^3)=0.031761;Z(x_4,y_4,t_1,{\mathrm{\λ}}_{1|0}^3)=0.245627;\\ Z(x_5,y_5,t_1,{\mathrm{\λ}}_{1|0}^3)=0.14569;Z(x_6,y_6,t_1,{\mathrm{\λ}}_{1|0}^3)=0.467618;\\ Z(x_7,y_7,t_1,{\mathrm{\λ}}_{1|0}^3)=0.154867;Z(x_8,y_8,t_1,{\mathrm{\λ}}_{1|0}^3)=0.04712.\end{array}$

当待求向量X为t₁＝6h时刻的第5个Sigma点时，t₁＝6h时刻各传感器节点处的理论水体离子浓度如下：

$\begin{array}{c}Z(x_1,y_1,t_1,{\mathrm{\λ}}_{1|0}^4)=0.38984;Z(x_2,y_2,t_1,{\mathrm{\λ}}_{1|0}^4)=0.216404;\\ Z(x_3,y_3,t_1,{\mathrm{\λ}}_{1|0}^4)=0.180664;Z(x_4,y_4,t_1,{\mathrm{\λ}}_{1|0}^4)=0.636144;\\ Z(x_5,y_5,t_1,{\mathrm{\λ}}_{1|0}^4)=0.411627;Z(x_6,y_6,t_1,{\mathrm{\λ}}_{1|0}^4)=0.924261;\\ Z(x_7,y_7,t_1,{\mathrm{\λ}}_{1|0}^4)=0.44533;Z(x_8,y_8,t_1,{\mathrm{\λ}}_{1|0}^4)=0.207802.\end{array}$

表示当待求向量X为t₁＝6h时刻的第q+1个Sigma点时，t₁＝6h时刻第i个传感器节点(x_i,y_i)处的理论水体离子浓度；

当待求向量X为t₁＝6h时刻的第6个Sigma点时，t₁＝6h时刻各传感器节点处的理论水体离子浓度如下：

$\begin{array}{c}Z(x_1,y_1,t_1,{\mathrm{\λ}}_{1|0}^5)=0.205871;Z(x_2,y_2,t_1,{\mathrm{\λ}}_{1|0}^5)=0.124709;\\ Z(x_3,y_3,t_1,{\mathrm{\λ}}_{1|0}^5)=0.105074;Z(x_4,y_4,t_1,{\mathrm{\λ}}_{1|0}^5)=0.357685;\\ Z(x_5,y_5,t_1,{\mathrm{\λ}}_{1|0}^5)=0.288039;Z(x_6,y_6,t_1,{\mathrm{\λ}}_{1|0}^5)=0.665988;\\ Z(x_7,y_7,t_1,{\mathrm{\λ}}_{1|0}^5)=0.281433;Z(x_8,y_8,t_1,{\mathrm{\λ}}_{1|0}^5)=0.135239.\end{array}$

当待求向量X为t₁＝6h时刻的第7个Sigma点时，t₁＝6h时刻各传感器节点处的理论水体离子浓度如下：

$\begin{array}{c}Z(x_1,y_1,t_1,{\mathrm{\λ}}_{1|0}^6)=0.360151;Z(x_2,y_2,t_1,{\mathrm{\λ}}_{1|0}^6)=0.581824;\\ Z(x_3,y_3,t_1,{\mathrm{\λ}}_{1|0}^6)=0.653197;Z(x_4,y_4,t_1,{\mathrm{\λ}}_{1|0}^6)=0.254905;\end{array}$

$\begin{array}{c}Z(x_5,y_5,t_1,{\mathrm{\λ}}_{1|0}^6)=0.487653;Z(x_6,y_6,t_1,{\mathrm{\λ}}_{1|0}^6)=0.216877;\\ Z(x_7,y_7,t_1,{\mathrm{\λ}}_{1|0}^5)=0.426542;Z(x_8,y_8,t_1,{\mathrm{\λ}}_{1|0}^6)=0.917164.\end{array}$

当待求向量X为t₁＝6h时刻的第8个Sigma点时，t₁＝6h时刻各传感器节点处的理论水体离子浓度如下：

$\begin{array}{c}Z(x_1,y_1,t_1,{\mathrm{\λ}}_{1|0}^7)=0.62234;Z(x_2,y_2,t_1,{\mathrm{\λ}}_{1|0}^7)=0.427979;\\ Z(x_3,y_3,t_1,{\mathrm{\λ}}_{1|0}^7)=0.383123;Z(x_4,y_4,t_1,{\mathrm{\λ}}_{1|0}^7)=0.846294;\\ Z(x_5,y_5,t_1,{\mathrm{\λ}}_{1|0}^7)=0.601173;Z(x_6,y_6,t_1,{\mathrm{\λ}}_{1|0}^7)=1.048705;\\ Z(x_7,y_7,t_1,{\mathrm{\λ}}_{1|0}^7)=0.64922;Z(x_8,y_8,t_1,{\mathrm{\λ}}_{1|0}^7)=0.399828.\end{array}$

当待求向量X为t₁＝6h时刻的第9个Sigma点时，t₁＝6h时刻各传感器节点处的理论水体离子浓度如下：

$\begin{array}{c}Z(x_1,y_1,t_1,{\mathrm{\λ}}_{1|0}^8)=0.219506;Z(x_2,y_2,t_1,{\mathrm{\λ}}_{1|0}^8)=0.12185;\\ Z(x_3,y_3,t_1,{\mathrm{\λ}}_{1|0}^8)=0.101726;Z(x_4,y_4,t_1,{\mathrm{\λ}}_{1|0}^8)=0.358191;\\ Z(x_5,y_5,t_1,{\mathrm{\λ}}_{1|0}^8)=0.231773;Z(x_6,y_6,t_1,{\mathrm{\λ}}_{1|0}^8)=0.52042;\\ Z(x_7,y_7,t_1,{\mathrm{\λ}}_{1|0}^8)=0.25075;Z(x_8,y_8,t_1,{\mathrm{\λ}}_{1|0}^8)=0.117006.\end{array}$

预测的t₁＝6h时刻各传感器节点处的水体离子浓度如下：

$\begin{array}{c}{\left(Z_1^{-}\right)}_0=W_{0}Z(x_1,y_1,t_1,{\mathrm{\λ}}_{1|0}^0)+\underset{p=1}{\overset{4}{\mathrm{\Σ}}}{W}_{p}Z(x_1,y_1,t_1,{\mathrm{\λ}}_{1|0}^p)+\underset{q=5}{\overset{8}{\mathrm{\Σ}}}{W}_{q}Z(x_1,y_1,t_1,{\mathrm{\λ}}_{1|0}^q)=0.25305;\\ {\left(Z_1^{-}\right)}_2=W_{0}Z(x_2,y_2,t_1,{\mathrm{\λ}}_{1|0}^0)+\underset{p=1}{\overset{4}{\mathrm{\Σ}}}{W}_{p}Z(x_2,y_2,t_1,{\mathrm{\λ}}_{1|0}^p)+\underset{q=5}{\overset{8}{\mathrm{\Σ}}}{W}_{q}Z(x_2,y_2,t_1,{\mathrm{\λ}}_{1|0}^q)=0.185664;\\ {\left(Z_1^{-}\right)}_3=W_{0}Z(x_3,y_3,t_1,{\mathrm{\λ}}_{1|0}^0)+\underset{p=1}{\overset{4}{\mathrm{\Σ}}}{W}_{p}Z(x_3,y_3,t_1,{\mathrm{\λ}}_{1|0}^p)+\underset{q=5}{\overset{8}{\mathrm{\Σ}}}{W}_{q}Z(x_3,y_3,t_1,{\mathrm{\λ}}_{1|0}^q)=0.174051;\\ {\left(Z_1^{-}\right)}_4=W_{0}Z(x_4,y_4,{t_1,\mathrm{\λ}}_{1|0}^0)+\underset{p=1}{\overset{4}{\mathrm{\Σ}}}{W}_{p}Z(x_4,y_4,t_1,{\mathrm{\λ}}_{1|0}^p)+\underset{q=5}{\overset{8}{\mathrm{\Σ}}}{W}_{q}Z(x_4,y_4,t_1,{\mathrm{\λ}}_{1|0}^q)=0.370418;\\ {\left(Z_1^{-}\right)}_5=W_{0}Z(x_5,y_5,t_1,{\mathrm{\λ}}_{1|0}^0)+\underset{p=1}{\overset{4}{\mathrm{\Σ}}}{W}_{p}Z(x_5,y_5,t_1,{\mathrm{\λ}}_{1|0}^p)+\underset{q=5}{\overset{8}{\mathrm{\Σ}}}{W}_{q}Z(x_5,y_5,t_1,{\mathrm{\λ}}_{1|0}^q)=0.282525;\end{array}$

$\begin{array}{c}{\left(Z_1^{-}\right)}_6=W_{0}Z(x_6,y_6,t_1,{\mathrm{\λ}}_{1|0}^0)+\underset{p=1}{\overset{4}{\mathrm{\Σ}}}{W}_{p}Z(x_6,y_6,t_1,{\mathrm{\λ}}_{1|0}^p)+\underset{q=5}{\overset{8}{\mathrm{\Σ}}}{W}_{q}Z(x_6,y_6,t_1,{\mathrm{\λ}}_{1|0}^q)=0.53195;\\ {\left(Z_1^{-}\right)}_7=W_{0}Z(x_7,y_7,t_1,{\mathrm{\λ}}_{1|0}^0)+\underset{p=1}{\overset{4}{\mathrm{\Σ}}}{W}_{p}Z(x_7,y_7,t_1,{\mathrm{\λ}}_{1|0}^p)+\underset{q=5}{\overset{8}{\mathrm{\Σ}}}{W}_{q}Z(x_7,y_7,t_1,{\mathrm{\λ}}_{1|0}^q)=0.291515;\\ {\left(Z_1^{-}\right)}_8=W_{0}Z(x_8,y_8,t_1,{\mathrm{\λ}}_{1|0}^0)+\underset{p=1}{\overset{4}{\mathrm{\Σ}}}{W}_{p}Z(x_8,y_8,t_1,{\mathrm{\λ}}_{1|0}^p)+\underset{q=5}{\overset{8}{\mathrm{\Σ}}}{W}_{q}Z(x_8,y_8,t_1,{\mathrm{\λ}}_{1|0}^q)=0.215504.\end{array}$

步骤3.3、计算t₁＝6h时刻的待求向量X的后验估计和t₁＝6h时刻的误差协方差矩阵首先计算t₁＝6h时刻的待求向量X的后验估计

${\hat{X}}_1={\hat{X}}_1^{-}+G_1\left(\begin{array}{c}\tilde{Z}(x_1,y_1,t_1)-{\left(Z_1^{-}\right)}_1\\ \tilde{Z}(x_2,y_2,t_1)-{\left(Z_1^{-}\right)}_2\\ \·\\ \·\\ \·\\ \tilde{Z}(x_8,y_8,t_1)-{\left(Z_1^{-}\right)}_8\end{array}\right)---\left(13\right)$

式(13)中：为t₁＝6h时刻的待求向量X的先验估计，

${\hat{X}}_1^{-}={\hat{X}}_0---\left(14\right)$

${\hat{X}}_1^{-}={\left(\mathrm{2,8,3,80}\right)}^T;$

G₁是t₁＝6h时刻的增益矩阵，G₁用于修正t₁＝6h时刻的待求向量X的先验估计得到t₁＝6h时刻的待求向量X的后验估计

$G_1=P_{\mathrm{XZ}}{S}_1^{-1}---\left(15\right)$

和依次表示第1个传感器节点、第2个传感器节点和第8个传感器节点在t₁＝6h时刻测量的水体离子浓度；

和依次表示预测的t₁＝6h时刻第1个传感器节点(x₁,y₁)处、第2个传感器节点(x₂,y₂)处和第8个传感器节点(x₈,y₈)处的水体离子浓度。

式(14)中：为t₀＝0h时刻的待求向量X的后验估计。

式(15)中：

$P_{\mathrm{XZ}}=W_0({\mathrm{\λ}}_{1|0}^0-{\hat{X}}_1^{-}){\left(\mathrm{\Δ}{Z}_0\right)}^T+\underset{p=1}{\overset{n}{\mathrm{\Σ}}}{W}_p({\mathrm{\λ}}_{1|0}^p-{\hat{X}}_1^{-}){\left(\mathrm{\Δ}{Z}_p\right)}^T+\underset{q=5}{\overset{8}{\mathrm{\Σ}}}{W}_q({\mathrm{\λ}}_{1|0}^q-{\hat{X}}_1^{-}){\left(\mathrm{\Δ}{Z}_q\right)}^T---\left(16\right)$

S₁＝R₁+P_ZZ (17)式(16)：W₀、W_p和W_q依次同式(9)、式(10)和式(11)；

和依次同式(6)、式(7)和式(8)；

为t₁＝6h时刻的待求向量X的先验估计，同式(14)；

$\mathrm{\Δ}{Z}_0=\left(\begin{array}{c}Z(x_1,y_1,t_1,{\mathrm{\λ}}_{1|0}^0)-{\left(Z_1^{-}\right)}_1\\ Z(x_2,y_2,t_1,{\mathrm{\λ}}_{1|0}^0)-{\left(Z_1^{-}\right)}_2\\ \·\\ \·\\ \·\\ Z(x_8,y_8,t_1,{\mathrm{\λ}}_{1|0}^0)-{\left(Z_1^{-}\right)}_8\end{array}\right)---\left(18\right)$

ΔZ₀＝(0.051623,-0.01654,-0.03286,0.12675,0.039175,0.19039,0.056525,-0.0531)^T。

$\mathrm{\Δ}{Z}_p=\left(\begin{array}{c}Z(x_1,y_1,t_1,{\mathrm{\λ}}_{1|0}^p)-{\left(Z_1^{-}\right)}_1\\ Z(x_2,y_2,t_1,{\mathrm{\λ}}_{1|0}^p)-{\left(Z_1^{-}\right)}_2\\ \·\\ \·\\ \·\\ Z(x_8,y_8,t_1,{\mathrm{\λ}}_{1|0}^p)-{\left(Z_1^{-}\right)}_8\end{array}\right)---\left(19\right)$

ΔZ₁＝(-0.2498,-0.18346,-0.17211,-0.36392,-0.26848,-0.50634,-0.28282,-0.21013)^T；

ΔZ₂＝(-0.25074,-0.18482,-0.17342,-0.36592,-0.28067,-0.5266,-0.28929,-0.2148)^T；

ΔZ₃＝(-0.14517,-0.14309,-0.14229,-0.12479,-0.13684,-0.06433,-0.13665,-0.16838)^T；

ΔZ₄＝(0.13679,0.03074,0.006614,0.265726,0.129102,0.392311,0.153815,-0.0077)^T。

$\mathrm{\Δ}{Z}_q=\left(\begin{array}{c}Z(x_1,y_1,t_1,{\mathrm{\λ}}_{1|0}^q)-{\left(Z_1^{-}\right)}_1\\ Z(x_2,y_2,t_1,{\mathrm{\λ}}_{1|0}^q)-{\left(Z_1^{-}\right)}_2\\ \·\\ \·\\ \·\\ Z(x_8,y_8,t_1,{\mathrm{\λ}}_{1|0}^q)-{\left(Z_1^{-}\right)}_8\end{array}\right)---\left(20\right)$

ΔZ₅＝(-0.03718,-0.06095,-0.06898,-0.01273,0.005514,0.134037,-0.01008,-0.08027)^T；

ΔZ₆＝(0.1071,0.396159,0.479146,-0.11551,0.205128,-0.31507,0.135028,0.70166)^T；

ΔZ₇＝(0.36929,0.242315,0.209072,0.475876,0.318648,0.516755,0.357705,0.184324)^T；

ΔZ₈＝(-0.03354,-0.06381,-0.07232,-0.01223,-0.05075,-0.01153,-0.04076,-0.0985)^T。

式(18)中：同式(6)；

表示当待求向量X为t₁＝6h时刻的第1个Sigma点时，t₁＝6h时刻第1个传感器节点(x₁,y₁)处的理论水体离子浓度；

表示当待求向量X为t₁＝6h时刻的第1个Sigma点时，t₁＝6h时刻第2个传感器节点(x₂,y₂)处的理论水体离子浓度；

表示当待求向量X为t₁＝6h时刻的第1个Sigma点时，t₁＝6h时刻第8个传感器节点(x₈,y₈)处的理论水体离子浓度；

和依次表示预测的t₁＝6h时刻第1个传感器节点(x₁,y₁)处、第2个传感器节点(x₂,y₂)处和第8个传感器节点(x₈,y₈)处的水体离子浓度。

式(19)中：同式(7)；

表示当待求向量X为t₁＝6h时刻的第p+1个Sigma点时，t₁＝6h时刻第1个传感器节点(x₁,y₁)处的理论水体离子浓度；

表示当待求向量X为t₁＝6h时刻的第p+1个Sigma点时，t₁＝6h时刻第2个传感器节点(x₂,y₂)处的理论水体离子浓度；

表示当待求向量X为t₁＝6h时刻的第p+1个Sigma点时，t₁＝6h时刻第8个传感器节点(x₈,y₈)处的理论水体离子浓度；

和依次表示预测的t₁＝6h时刻第1个传感器节点(x₁,y₁)处、第2个传感器节点(x₂,y₂)处和第8个传感器节点(x₈,y₈)处的水体离子浓度。

式(20)中：同式(8)；

表示当待求向量X为t₁＝6h时刻的第q+1个Sigma点时，t₁＝6h时刻第1个传感器节点(x₁,y₁)处的理论水体离子浓度；

表示当待求向量X为t₁＝6h时刻的第q+1个Sigma点时，t₁＝6h时刻第2个传感器节点(x₂,y₂)处测量的理论水体离子浓度；

表示当待求向量X为t₁＝6h时刻的第q+1个Sigma点时，t₁＝6h时刻第8个传感器节点(x₈,y₈)处的理论水体离子浓度；

和依次表示预测的t₁＝6h时刻第1个传感器节点(x₁,y₁)处、第2个传感器节点(x₂,y₂)处和第8个传感器节点(x₈,y₈)处的水体离子浓度。

式(17)中：R₁是测量噪声的协方差矩阵，所述测量噪声的协方差矩阵为对角矩阵，

$R_1={\left(\begin{array}{cccc}1& 0& \·\·\·& 0\\ 0& 1& \·\·\·& 0\\ \·& \·& \mathrm{}& \·\\ \·& \·& \mathrm{}& \·\\ \·& \·& \mathrm{}& \·\\ 0& 0& \·\·\·& 1\end{array}\right)}_{8\×8}---\left(21\right)$

$P_{\mathrm{ZZ}}=W_0\left(\mathrm{\Δ}{Z}_0\right){\left(\mathrm{\Δ}{Z}_0\right)}^T+\underset{p=1}{\overset{n}{\mathrm{\Σ}}}{W}_p\left(\mathrm{\Δ}{Z}_p\right){\left(\mathrm{\Δ}{Z}_p\right)}^T+\underset{q=5}{\overset{8}{\mathrm{\Σ}}}{W}_q\left(\mathrm{\Δ}{Z}_q\right){\left(\mathrm{\Δ}{Z}_p\right)}^T---\left(22\right)$

式(22)中：W₀、W_p和W_q依次同式(9)、式(10)和式(11)；

ΔZ₀、ΔZ_p和ΔZ_q依次同式(18)、式(19)和式(20)。

$P_{\mathrm{XZ}}=\left(\begin{array}{cccccccc}-0.09533& -0.05492& -0.04624& -0.15745& -0.12284& -0.2871& -0.12228& -0.05822\\ -0.16043& -0.26047& -0.29257& -0.11226& -0.2178& -0.09484& -0.19023& -0.41088\\ -0.23065& -0.17279& -0.15753& -0.2693& -0.20421& -0.26052& -0.22163& -0.15813\\ 0.380918& 0.211451& 0.176529& 0.621584& 0.402205& 0.903106& 0.435137& 0.203046\end{array}\right);$

$P_{\mathrm{ZZ}}=\left(\begin{array}{cccccccc}0.032074& 0.025176& 0.023816& 0.041445& 0.032015& 0.049368& 0.03382& 0.027376\\ 0.025176& 0.031323& 0.033458& 0.022733& 0.028476& 0.019815& 0.027115& 0.043768\\ 0.023816& 0.033458& 0.036543& 0.018318& 0.028083& 0.012637& 0.025872& 0.048821\\ 0.041445& 0.022733& 0.018318& 0.062476& 0.039021& 0.081824& 0.043628& 0.016947\\ 0.032015& 0.028476& 0.028083& 0.039021& 0.033554& 0.045947& 0.03438& 0.034181\\ 0.049368& 0.019815& 0.012637& 0.081824& 0.045947& 0.114865& 0.052762& 0.007166\\ 0.03382& 0.027115& 0.025872& 0.043628& 0.03438& 0.052762& 0.036034& 0.030289\\ 0.027376& 0.043768& 0.048821& 0.016947& 0.034181& 0.007166& 0.030289& 0.066679\end{array}\right);$

$S_1=\left(\begin{array}{cccccccc}0.132074& 0.025176& 0.023816& 0.041445& 0.032015& 0.049368& 0.03382& 0.027376\\ 0.025176& 0.131323& 0.033458& 0.022733& 0.028476& 0.019815& 0.027115& 0.043768\\ 0.023816& 0.033458& 0.136543& 0.018318& 0.028083& 0.012637& 0.025872& 0.048821\\ 0.041445& 0.022733& 0.018318& 0.162476& 0.039021& 0.081824& 0.043628& 0.016947\\ 0.032015& 0.028476& 0.028083& 0.039021& 0.133554& 0.045947& 0.03438& 0.034181\\ 0.049368& 0.019815& 0.012637& 0.081824& 0.045947& 0.214865& 0.052762& 0.007166\\ 0.03382& 0.027115& 0.025872& 0.043628& 0.03438& 0.052762& 0.136034& 0.030289\\ 0.027376& 0.043768& 0.048821& 0.01647& 0.034181& 0.007166& 0.030289& 0.166679\end{array}\right);$

$G_1=\left(\begin{array}{cccccccc}-0.05306& -0.02173& -0.02051& -0.23301& -0.35318& -1.08792& -0.26062& -0.13862\\ -0.26911& -0.88585& -1.07656& 0.04595& -0.63915& 0.051271& -0.45221& -1.66658\\ -0.83763& -0.52728& -0.42212& -0.81684& -0.54311& -0.35225& -0.6505& -0.22125\\ 0.657762& 0.198536& 0.13095& 1.592045& 0.954352& 2.960267& 0.991823& 0.354588\end{array}\right);$

${\hat{X}}_1={\left(\mathrm{1.783544,8.199684,3.343394,80.13257}\right)}^T.$

然后计算误差协方差矩阵

${\hat{P}}_1=P_1^{-}-G_1{S}_1{\left(G_1\right)}^T---\left(23\right)$

式(23)中： $P_1^{-}{\hat{P}}_0---\left(24\right)$

$P_1^{-}=\mathrm{diag}\left(\mathrm{4,4,4,100}\right);$

G₁和S₁依次同式(15)和式(17)；

式(24)中：为t₀＝0h时刻的误差协方差矩阵；

${\hat{P}}_1=\left(\begin{array}{cccc}3.560445& -0.33296& -0.5172& 1.439375\\ -0.33296& 2.511154& -0.85327& 1.197235\\ -0.5172& -0.85327& 3.047394& 1.877354\\ 1.439375& 1.197235& 1.877354& 95.13391\end{array}\right).$

故：t₁＝6h时刻的待求向量X的后验估计 ${\hat{X}}_1={\left(\mathrm{1.783544,8.199684,3.343394,80.13257}\right)}^T;$

t₁＝6h时刻的误差协方差矩阵 ${\hat{P}}_1=\left(\begin{array}{cccc}3.560445& -0.33296& -0.5172& 1.439375\\ -0.33296& 2.511154& -0.85327& 1.197235\\ -0.5172& -0.85327& 3.047394& 1.877354\\ 1.439375& 1.197235& 1.877354& 95.13391\end{array}\right).$

步骤4、k＝2，同理步骤3，得到t₂＝10h时刻：

待求向量X的后验估计 ${\hat{X}}_2={\left(\mathrm{1.683838,8.038214,3.561119,82.43101}\right)}^T;$

误差协方差矩阵 ${\hat{P}}_2=\left(\begin{array}{cccc}3.270747& -0.60335& -0.25665& 1.583841\\ -0.60335& 0.85556& -0.40224& 1.233375\\ -0.25665& -0.40224& 2.746954& 1.749259\\ 1.583841& 1.233375& 1.749259& 94.46743\end{array}\right).$

……，

k＝5，同理步骤3，得到t₅＝18h时刻：

待求向量X的后验估计 ${\hat{X}}_5={\left(\mathrm{2.559427,7.21534,4.065863,82.51426}\right)}^T;$

误差协方差矩阵 ${\hat{P}}_5=\left(\begin{array}{cccc}0.744604& -0.23172& 0.054246& 1.655629\\ -0.23172& 0.213794& -0.1765& 1.163522\\ 0.054246& -0.1765& 2.397213& 2.187048\\ 1.655629& 1.163522& 2.187048& 92.45414\end{array}\right).$

步骤5、在t₅＝18h时刻的待求向量X的后验估计 ${\hat{X}}_5={\left(\mathrm{2.559427,7.21534,4.065863,82.51426}\right)}^T$ 中，2.559427为所求的污染源的横坐标，7.21534为所求的污染源的纵坐标。

实施例2

一种靠近湖泊岸边的污染源定位方法。该方法的具体步骤是：

步骤1、以岸边为纵轴，以垂直于岸边为横轴，在水平面上建立直角坐标系。设污染源在(α,β)处从τ时刻连续释放污染物，污染源的质量流率为Q，以污染源的横坐标α、污染源的纵坐标β、污染源的初始扩散时间τ和污染源的质量流率Q为未知量组成待求向量X，待求向量X＝(α,β,τ,Q)^T；传感器节点的总个数为8，传感器节点共测量了9个时刻的水体离子浓度，8个传感器节点(x_i,y_i)(i＝1,...,8)在t_k(k＝1,...,9)时刻测量的水体离子浓度见表2.1。

表2.1

第i(i＝1,...,8)个传感器节点(x_i,y_i)在t_k(k＝1,...,9)时刻测量的水体离子浓度为：

$\tilde{Z}(x_i,y_i{t}_k)=Z(x_i,y_i,t_{k}X)+V_i^{\left(k\right)}---\left(1\right)$

式(1)中：是第i个传感器节点(x_i,y_i)在t_k时刻的测量噪声，kg/m³；

的均值为0，的方差为D_i，(kg/m³)²；8个传感器节(x_i,y_i)(i＝1,...,8)在t_k(k＝1,...,9)时刻的测量噪声的方差见表2.2；

Z(x_i,y_i,t_k,X)表示t_k时刻第i个传感器节点(x_i,y_i)处的理论水体离子浓度，kg/m³；

$Z(x_i,y_i,t_k,X)=\frac{\mathrm{\ξQ}}{2f\sqrt{\mathrm{\πd}}}[\frac{1}{R_i}\mathrm{erfc}\left(\frac{R_i}{2\sqrt{d(t_k-\mathrm{\τ})}}\right)+\frac{1}{{\stackrel{\‾}{R}}_i}\mathrm{erfc}\left(\frac{{\stackrel{\‾}{R}}_i}{2\sqrt{d(t_k-\mathrm{\τ})}}\right)]---\left(2\right)$

式(2)中：ξ为补量纲常量，tunit为时间单位，h；

R_i为第i个传感器节点(x_i,y_i)到污染源(α,β)的距离，m；

$R_i=\sqrt{{(x_i-\mathrm{\α})}^2+{(y_i-\mathrm{\β})}^2}---\left(3\right)$

为第i个传感器节点(x_i,y_i)到镜像污染源(-α,β)的距离，m；

${\stackrel{\‾}{R}}_i=\sqrt{{(x_i+\mathrm{\α})}^2+{(y_i-\mathrm{\β})}^2}---\left(4\right)$

erfc(x)=1-erf(x) (5)

f为水体平均深度，f=10m；

d为同向扩散系数，d=0.5m²/h

表2.2

步骤2、同实施例1的步骤2。

步骤3、同实施例1的步骤3。

步骤4、k＝2，同理步骤3，得到t₂＝10h时刻：

待求向量X的后验估计 ${\hat{X}}_2={\left(\mathrm{1.683838,8.038214,3.561119,82.43101}\right)}^T;$

误差协方差矩阵 ${\hat{P}}_2=\left(\begin{array}{cccc}3.270747& -0.60335& -0.25665& 1.583841\\ -0.60335& 0.85556& -0.40224& 1.233375\\ -0.25665& -0.40224& 2.746954& 1.749259\\ 1.583841& 1.233375& 1.749259& 94.46743\end{array}\right).$

……。

k＝9，同理步骤3，得到t₉＝30h时刻：

待求向量X的后验估计 ${\hat{X}}_9={\left(\mathrm{2.93044,6.062932,5.648708,85.69018}\right)}^T;$

误差协方差矩阵 ${\hat{P}}_9=\left(\begin{array}{cccc}-0.00139& -0.00785& 0.03308& 0.70294\\ -0.00785& 0.044006& -0.03705& 0.702772\\ 0.03308& -0.03705& 2.159088& 3.178541\\ 0.70294& 0.702772& 3.178541& 60.2314\end{array}\right).$

步骤5、在t₉＝30h时刻的待求向量X的后验估计 ${\hat{X}}_9={\left(\mathrm{2.93044,6.062932,5.648708,85.69018}\right)}^T;$ 中，2.93044为所求的污染源的横坐标，6.062932为所求的污染源的纵坐标。

本具体实施方式与现有技术相比具有如下经济效果：

本具体实施方式无需提前给定阈值，用每个时刻的测量值去更新待求向量X的先验估计得到待求向量X的后验估计，在实际定位操作中简单易行；其次，本具体实施方式中测量的浓度数据可以包含整个扩散过程的信息，避免了因忽略过渡阶段影响定位精度的问题；然后，本具体实施方式在每个测量时刻均处理了所有传感器节点(x_i,y_i)(i＝1,...,N)的测量数据，由于在每个测量时刻都考虑到了所有传感器节点(x_i,y_i)(i＝1,...,N)的测量噪声该方法的定位效果相对稳定。

本具体实施方式在定位计算过程中只需保存当前时刻的测量值、上一时刻的待求向量X的后验估计和上一时刻的误差协方差矩阵，用当前时刻的测量值去更新待求向量X的先验估计，计算过程中不需要保存和计算历史测量数据。故本具体实施方式具有工作效率高和耗时短的优点。

因此，本具体实施方式具有易实现、定位精度高、工作效率高和耗时短的特点。

Claims (2)

1.一种靠近湖泊岸边的污染源定位方法，其特征在于所述污染源定位方法的步骤是：

步骤1、以岸边为纵轴，以垂直于岸边为横轴，在水平面上建立直角坐标系；设污染源在(α,β)处从τ时刻连续释放污染物，污染源的质量流率为Q，以污染源的横坐标α、污染源的纵坐标β、污染源的初始扩散时间τ和污染源的质量流率Q为未知量组成待求向量X，待求向量X＝(α,β,τ,Q)^T；传感器节点的总个数为N，传感器节点共测量了L个时刻的水体离子浓度，第i(i＝1,...,N)个传感器节点(x_i,y_i)在t_k(k＝1,...,L)时刻测量的水体离子浓度为：

$\tilde{Z}(x_i,y_i{t}_k)=Z(x_i,y_i,t_{k}X)+V_i^{\left(k\right)}---\left(1\right)$

式(1)中：是第i个传感器节点(x_i,y_i)在t_k时刻的测量噪声，kg/m³，

的均值为0，的方差为D_i，(kg/m³)²，

Z(x_i,y_i,t_k,X)表示t_k时刻第i个传感器节点(x_i,y_i)处的理论水体离子浓度，

kg/m³，

$Z(x_i,y_i,t_k,X)=\frac{\mathrm{\ξQ}}{2f\sqrt{\mathrm{\πd}}}[\frac{1}{R_i}\mathrm{erfc}\left(\frac{R_i}{2\sqrt{d(t_k-\mathrm{\τ})}}\right)+\frac{1}{{\stackrel{\‾}{R}}_i}\mathrm{erfc}\left(\frac{{\stackrel{\‾}{R}}_i}{2\sqrt{d(t_k-\mathrm{\τ})}}\right)]---\left(2\right)$

式(2)中：ξ为补量纲常量，tunit为时间单位，h，

R_i为第i个传感器节点(x_i,y_i)到污染源(α,β)的距离，m，

$R_i=\sqrt{{(x_i-\mathrm{\α})}^2+{(y_i-\mathrm{\β})}^2}---\left(3\right)$

为第i个传感器节点(x_i,y_i)到镜像污染源(-α,β)的距离，m，

${\stackrel{\‾}{R}}_i=\sqrt{{(x_i+\mathrm{\α})}^2+{(y_i-\mathrm{\β})}^2}---\left(4\right)$

erfc(x)＝1-erf(x) (5)

f为水体平均深度，m，

d为同向扩散系数，m²/h；

步骤2、设t₀时刻的待求向量X的后验估计为设t₀时刻的误差协方差矩阵为 ${\hat{P}}_0=\mathrm{diag}({\widehat{\mathrm{\σ}}}_{\mathrm{\α}\left(0\right)}^2,{\widehat{\mathrm{\σ}}}_{\mathrm{\β}\left(0\right)}^2,{\widehat{\mathrm{\σ}}}_{\mathrm{\τ}\left(0\right)}^2,{\widehat{\mathrm{\σ}}}_{Q\left(0\right)}^2),$ 其中：和依次为t₀时刻的污染源的横坐标的估计方差、t₀时刻的污染源的纵坐标的估计方差、t₀时刻的污染源的初始扩散时间的估计方差和t₀时刻的污染源的质量流率的估计方差；令k＝1；

步骤3、将t_k时刻的N个传感器节点测量的水体离子浓度t_k-1时刻的待求向量X的后验估计和t_k-1时刻的误差协方差矩阵用动态数据处理模块进行处理，得到t_k时刻的待求向量X的后验估计和t_k时刻的误差协方差矩阵

步骤4、k：＝k+1，跳转到步骤3，直至计算出k＝L时的t_L时刻的待求向量X的后验估计 ${\hat{X}}_L={({\widehat{\mathrm{\α}}}_{\left(L\right)},{\widehat{\mathrm{\β}}}_{\left(L\right)},{\widehat{\mathrm{\τ}}}_{\left(L\right)},{\hat{Q}}_{\left(L\right)})}^T$ 和误差协方差矩阵

步骤5、在t_L时刻的待求向量X的后验估计中，为所求的污染源的横坐标，为所求的污染源的纵坐标。

2.根据权利要求1所述的靠近湖泊岸边的污染源定位方法，其特征在于所述的用动态数据处理模块进行处理的步骤是：

步骤(1)、由t_k-1时刻的待求向量X的后验估计和t_k-1时刻的误差协方差矩阵得到t_k时刻的第1个Sigma点第p+1个Sigma点和第q+1个Sigma点 ${\mathrm{\λ}}_{k|k-1}^q(q=n+1,...,2n):$

${\mathrm{\λ}}_{k|k-1}^0={\hat{X}}_{k-1}---\left(6\right)$

${\mathrm{\λ}}_{k|k-1}^p={\hat{X}}_{k-1}+{\left[{\left(\sqrt{(n+1){\hat{P}}_{k-1}}\right)}_p\right]}^T---\left(7\right)$

${\mathrm{\λ}}_{k|k-1}^q={\hat{X}}_{k-1}-{\left[{\left(\sqrt{(n+1){\hat{P}}_{k-1}}\right)}_{(q-n)}\right]}^T---\left(8\right)$

式(7)中：表示开方矩阵的第p行；

式(8)中：表示开方矩阵的第(q-n)行；

计算t_k时刻第1个Sigma点第p+1个Sigma点和第q+1个Sigma点依次对应的权值W₀、W_p和W_q：

$W_0=\frac{1}{n+1}---\left(9\right)$

$W_p=\frac{1}{2(n+1)}---\left(10\right)$

$W_q=\frac{1}{2(n+1)}---\left(11\right)$

式(7)(8)(9)(10)(11)中：n为待求向量X的维数，n＝4；

步骤(2)、预测t_k时刻各传感器节点(x_i,y_i)(i＝1,...,N)处的水体离子浓度

${\left(Z_k^{-}\right)}_i=W_{0}Z(x_i,y_i,t_k,{\mathrm{\λ}}_{k|k-1}^0)+\underset{p=1}{\overset{n}{\mathrm{\Σ}}}{W}_{p}Z(x_i,y_i,t_k,{\mathrm{\λ}}_{k|k-1}^p)+\underset{q=n+1}{\overset{2n}{\mathrm{\Σ}}}{W}_{q}Z(x_i,y_i,t_k,{\mathrm{\λ}}_{k|k-1}^q)---\left(12\right)$

式(12)：W₀、W_p和W_q依次同式(9)、式(10)和式(11)，

和依次同式(6)、式(7)和式(8)，

n为待求向量X的维数，n＝4,表示当待求向量X为t_k时刻的第1个Sigma点时，t_k时刻第i个传感器节点(x_i,y_i)处的理论水体离子浓度，

表示当待求向量X为t_k时刻的第p+1个Sigma点时，t_k时刻第i个传感器节点(x_i,y_i)处的理论水体离子浓度，

表示当待求向量X为t_k时刻的第q+1个Sigma点时，t_k时刻第i个传感器节点(x_i,y_i)处的理论水体离子浓度，

步骤(3)、计算t_k时刻的待求向量X的后验估计和t_k时刻的误差协方差矩阵首先计算t_k时刻的待求向量X的后验估计

${\hat{X}}_k={\hat{X}}_k^{-}+G_k\left(\begin{array}{c}\tilde{Z}(x_1,y_1,t_k)-{\left(Z_k^{-}\right)}_1\\ \tilde{Z}(x_2,y_2,t_k)-{\left(Z_k^{-}\right)}_2\\ \·\\ \·\\ \·\\ \tilde{Z}(x_N,y_N,t_k)-{\left(Z_k^{-}\right)}_N\end{array}\right)---\left(13\right)$

式(13)中：为t_k时刻的待求向量X的先验估计，

${\hat{X}}_k^{-}={\hat{X}}_{k-1}---\left(14\right)$

G_k是t_k时刻的增益矩阵，G_k用于修正t_k时刻的待求向量X的先验估计得到t_k时刻的待求向量X的后验估计

$G_k=P_{\mathrm{XZ}}{S}_k^{-1}---\left(15\right)$

和依次表示第1个传感器节点、第2个传感器节点和第N个传感器节点在t_k时刻测量的水体离子浓度，

和依次表示预测的t_k时刻第1个传感器节点(x₁,y₁)处、第2个传感器节点(x₂,y₂)处和第N个传感器节点(x_N,y_N)处的水体离子浓度，

式(14)中：为t_k-1时刻的待求向量X的后验估计，

式(15)中：

$P_{\mathrm{XZ}}=W_0({\mathrm{\λ}}_{k|k-1}^0-{\hat{X}}_k^{-}){\left(\mathrm{\Δ}{Z}_0\right)}^T+\underset{p=1}{\overset{n}{\mathrm{\Σ}}}{W}_p({\mathrm{\λ}}_{k|k-1}^p-{\hat{X}}_k^{-}){\left(\mathrm{\Δ}{Z}_p\right)}^T+\underset{q=n+1}{\overset{2n}{\mathrm{\Σ}}}{W}_q({\mathrm{\λ}}_{k|k-1}^q-{\hat{X}}_k^{-}){\left(\mathrm{\Δ}{Z}_q\right)}^T---\left(16\right)$

S_k＝R_k+P_ZZ (17)式(16)：W₀、W_p和W_q依次同式(9)、式(10)和式(11)，

和依次同式(6)、式(7)和式(8)，

为t_k时刻的待求向量X的先验估计，同式(14),

n为待求向量X的维数，n＝4；

$\mathrm{\Δ}{Z}_0=\left(\begin{array}{c}Z(x_1,y_1,t_k,{\mathrm{\λ}}_{k|k-1}^0)-{\left(Z_k^{-}\right)}_1\\ Z(x_2,y_2,t_k,{\mathrm{\λ}}_{k|k-1}^0)-{\left(Z_k^{-}\right)}_2\\ \·\\ \·\\ \·\\ Z(x_N,y_N,t_k,{\mathrm{\λ}}_{k|k-1}^0)-{\left(Z_k^{-}\right)}_N\end{array}\right)---\left(18\right)$

$\mathrm{\Δ}{Z}_p=\left(\begin{array}{c}Z(x_1,y_1,t_k,{\mathrm{\λ}}_{k|k-1}^p)-{\left(Z_k^{-}\right)}_1\\ Z(x_2,y_2,t_k,{\mathrm{\λ}}_{k|k-1}^p)-{\left(Z_k^{-}\right)}_2\\ \·\\ \·\\ \·\\ Z(x_N,y_N,t_k,{\mathrm{\λ}}_{k|k-1}^p)-{\left(Z_k^{-}\right)}_N\end{array}\right)---\left(19\right)$

$\mathrm{\Δ}{Z}_q=\left(\begin{array}{c}Z(x_1,y_1,t_k,{\mathrm{\λ}}_{k|k-1}^q)-{\left(Z_k^{-}\right)}_1\\ Z(x_2,y_2,t_k,{\mathrm{\λ}}_{k|k-1}^q)-{\left(Z_k^{-}\right)}_2\\ \·\\ \·\\ \·\\ Z(x_N,y_N,t_k,{\mathrm{\λ}}_{k|k-1}^q)-{\left(Z_k^{-}\right)}_N\end{array}\right)---\left(20\right)$

式(18)中：同式(6)，

表示当待求向量X为t_k时刻的第1个Sigma点时，t_k时刻第1个传感器节点(x₁,y₁)处的理论水体离子浓度，

表示当待求向量X为t_k时刻的第1个Sigma点时，t_k时刻第2个传感器节点(x₂,y₂)处的理论水体离子浓度，

表示当待求向量X为t_k时刻的第1个Sigma点时，t_k时刻第N个传感器节点(x_N,y_N)处的理论水体离子浓度，

和依次表示预测的t_k时刻第1个传感器节点(x₁,y₁)处、第2个传感器节点(x₂,y₂)处和第N个传感器节点(x_N,y_N)处的水体离子浓度；

式(19)中：同式(7)，

表示当待求向量X为t_k时刻的第p+1个Sigma点时，t_k时刻第1个传感器节点(x₁,y₁)处的理论水体离子浓度，

表示当待求向量X为t_k时刻的第p+1个Sigma点时，t_k时刻第2个传感器节点(x₂,y₂)处的理论水体离子浓度，

表示当待求向量X为t_k时刻的第p+1个Sigma点时，t_k时刻第N个传感器节点(x_N,y_N)处的理论水体离子浓度，

和依次表示预测的t_k时刻第1个传感器节点(x₁,y₁)处、第2个传感器节点(x₂,y₂)处和第N个传感器节点(x_N,y_N)处的水体离子浓度；

式(20)中：同式(8)，

表示当待求向量X为t_k时刻的第q+1个Sigma点时，t_k时刻第1个传感器节点(x₁,y₁)处的理论水体离子浓度，

表示当待求向量X为t_k时刻的第q+1个Sigma点时，t_k时刻第2个传感器节点(x₂,y₂)处的理论水体离子浓度，

表示当待求向量X为t_k时刻的第q+1个Sigma点时，t_k时刻第N个传感器节点(x_N,y_N)处的理论水体离子浓度，

和依次表示预测的t_k时刻第1个传感器节点(x₁,y₁)处、第2个传感器节点(x₂,y₂)处和第N个传感器节点(x_N,y_N)处的水体离子浓度；

式(17)中：R_k是测量噪声的协方差矩阵，所述测量噪声的协方差矩阵为对角矩阵，

$R_k={\left(\begin{array}{cccc}{D}_1& 0& \·\·\·& 0\\ 0& D_2& \·\·\·& 0\\ \·& \·& \mathrm{}& \·\\ \·& \·& \mathrm{}& \·\\ \·& \·& \mathrm{}& \·\\ 0& 0& \·\·\·& D_N\end{array}\right)}_{N\×N}---\left(21\right)$

$P_{\mathrm{ZZ}}=W_0\left(\mathrm{\Δ}{Z}_0\right){\left(\mathrm{\Δ}{Z}_0\right)}^T+\underset{p=1}{\overset{n}{\mathrm{\Σ}}}{W}_p\left(\mathrm{\Δ}{Z}_p\right){\left(\mathrm{\Δ}{Z}_p\right)}^T+\underset{q=n+1}{\overset{2n}{\mathrm{\Σ}}}{W}_q\left(\mathrm{\Δ}{Z}_q\right){\left(\mathrm{\Δ}{Z}_p\right)}^T---\left(22\right)$

式(22)中：W₀、W_p和W_q依次同式(9)、式(10)和式(11)，

ΔZ₀、ΔZ_p和ΔZ_q依次同式(18)、式(19)和式(20),

n为待求向量X的维数，n＝4；

然后计算误差协方差矩阵

${\hat{P}}_k=P_k^{-}-G_k{S}_k{\left(G_k\right)}^T---\left(23\right)$

式(23)中： $P_k^{-}={\hat{P}}_{k-1}---\left(24\right)$

G_k和S_k依次同式(15)和式(17)；

式(24)中：为t_k-1时刻的误差协方差矩阵。