说明书一种靠近湖泊岸边的污染源定位方法
技术领域
本发明属于污染源定位技术领域。具体涉及一种靠近湖泊岸边的污染源定位方法。
背景技术
湖泊是工农业用水、居民饮用水的重要水源。近年来，水污染事件的频发已引起人类的广泛关注，湖泊水体的污染给社会、经济和人们的生活带来了不可估量的损失。由于湖泊水体的污染主要来自工业废水、废渣和生活垃圾的沿岸排放，因此，靠近湖泊岸边污染源的定位对?；に试淳哂兄匾庖?。
现有靠近湖泊岸边的污染源定位方法有：近似函数法和非线性最小二乘法。近似函数法是在已知污染源扩散模型的基础上，用分段函数近似污染源扩散模型中的非线性函数，得到分段扩散模型。分段扩散模型将扩散过程简化为不稳定阶段和稳定阶段。采用该方法时需提前给定浓度阈值以判定扩散处于哪个阶段，然后选择处于同一阶段的两个测量值计算出污染源的位置。然而在实际应用中该浓度阈值难以预先确定，该方法不易实现。并且该方法忽略了过渡阶段的浓度信息，影响了定位精度。该方法的定位效果严重依赖于被选取的少量传感器节点的测量值，所以此方法易受测量噪声的干扰，定位结果的稳定性差。
非线性最小二乘法是在实际的约束条件下，满足测量的浓度值与理论的浓度值之差最小时的参数估计。该方法在每个测量时刻都有历史时刻的测量值和当前时刻的测量值参与最小二乘数值运算，然而随着测量次数的增加，数据量会急剧增长，在不同的测量时刻有重复而繁杂的最小二乘数值计算，因此该方法数据存储量大、工作效率低和耗时长。
发明内容
本发明旨在克服现有技术中的不足，目的是提供一种易实现、定位精度高、工作效率高和耗时短的靠近湖泊岸边的污染源定位方法。
为实现上述目的，本发明采用的技术方案的具体步骤是：
步骤1、以岸边为纵轴，以垂直于岸边为横轴，在水平面上建立直角坐标系。设污染源在(α,β)处从τ时刻连续释放污染物，污染源的质量流率为Q，以污染源的横坐标α、污染源的纵坐标β、污染源的初始扩散时间τ和污染源的质量流率Q为未知量组成待求向量X，待求向量X＝(α,β,τ,Q)T；传感器节点的总个数为N，传感器节点共测量了L个时刻的水体离子浓 度，第i(i＝1,...,N)个传感器节点(xi,yi)在tk(k＝1,...,L)时刻测量的水体离子浓度为：
Z~(xi,yitk)=Z(xi,yi,tkX)+Vi(k)---(1)]]>
式(1)中：是第i个传感器节点(xi,yi)在tk时刻的测量噪声，kg/m3；
的均值为0，的方差为Di，(kg/m3)2；
Z(xi,yi,tk,X)表示tk时刻第i个传感器节点(xi,yi)处的理论水体离子浓度，kg/m3；
Z(xi,yi,tk,X)=ξQ2fπd[1Rierfc(Ri2d(tk-τ))+1R‾ierfc(R‾i2d(tk-τ))]---(2)]]>
式(2)中：ξ为补量纲常量，tunit为时间单位，h；
Ri为第i个传感器节点(xi,yi)到污染源(α,β)的距离，m；
Ri=(xi-α)2+(yi-β)2---(3)]]>
为第i个传感器节点(xi,yi)到镜像污染源(-α,β)的距离，m；
R‾i=(xi+α)2+(yi-β)2---(4)]]>
erfc(x)＝1-erf(x) (5)
f为水体平均深度，m；
d为同向扩散系数，m2/h。
步骤2、设t0时刻的待求向量X的后验估计为设t0时刻的误差协方差矩阵为P^0=diag(σ^α(0)2,σ^β(0)2,σ^τ(0)2,σ^Q(0)2),]]>其中：依次为t0时刻的污染源的横坐标的估计方差、t0时刻的污染源的纵坐标的估计方差、t0时刻的污染源的初始扩散时间的估计方差和t0时刻的污染源的质量流率的估计方差；令k＝1。
步骤3、将tk时刻的N个传感器节点测量的水体离子浓度时刻的待求向量X的后验估计和tk-1时刻的误差协方差矩阵用动态数据处理?？榻写?，得到tk时刻的待求向量X的后验估计和tk时刻的 误差协方差矩阵
步骤4、k：＝k+1，跳转到步骤3，直至计算出k＝L时的tL时刻的待求向量X的后验估计X^L=(α^(L),β^(L),τ^(L),Q^(L))T]]>和误差协方差矩阵
步骤5、在tL时刻的待求向量X的后验估计中，为所求的污染源的横坐标，为所求的污染源的纵坐标。
2、根据权利要求1所述的靠近湖泊岸边的污染源定位方法，其特征在于所述的用动态数据处理?？榻写淼牟街枋牵?
步骤(1)、由tk-1时刻的待求向量X的后验估计和tk-1时刻的误差协方差矩阵得到tk时刻的第1个Sigma点第p+1个Sigma点和第q+1个Sigma点λk|k-1q(q=n+1,...,2n):]]>
λk|k-10=X^k-1---(6)]]>
λk|k-1p=X^k-1+[((n+1)P^k-1)p]T---(7)]]>
λk|k-1q=X^k-1-[((n+1)P^k-1)(q-n)]T---(8)]]>
式(7)中：表示开方矩阵的第p行；
式(8)中：表示开方矩阵的第(q-n)行。
计算tk时刻第1个Sigma点第p+1个Sigma点和第q+1个Sigma点依次对应的权值W0、Wp和Wq：
W0=1n+1---(9)]]>
Wp=12(n+1)---(10)]]>
Wq=12(n+1)---(11)]]>
式(7)(8)(9)(10)(11)中：n为待求向量X的维数，n＝4。
步骤(2)、预测tk时刻各传感器节点(xi,yi)(i＝1,...,N)处的水体离子浓度
(Zk-)i=W0Z(xi,yi,tk,λk|k-10)+Σp=1nWpZ(xi,yi,tk,λk|k-1p)+Σq=n+12nWqZ(xi,yi,tk,λk|k-1q)---(12)]]>
式(12)：W0、Wp和Wq依次同式(9)、式(10)和式(11)；
和依次同式(6)、式(7)和式(8)；
n为待求向量X的维数，n＝4；
表示当待求向量X为tk时刻的第1个Sigma点时，tk时刻第i个传感器节点(xi,yi)处的理论水体离子浓度；
表示当待求向量X为tk时刻的第p+1个Sigma点时，tk时刻第i个传感器节点(xi,yi)处的理论水体离子浓度；
表示当待求向量X为tk时刻的第q+1个Sigma点时，tk时刻第i个传感器节点(xi,yi)处的理论水体离子浓度；
步骤(3)、计算tk时刻的待求向量X的后验估计和tk时刻的误差协方差矩阵首先计算tk时刻的待求向量X的后验估计
X^k=X^k-+GkZ~(x1,y1,tk)-(Zk-)1Z~(x2,y2,tk)-(Zk-)2···Z~(xN,yN,tk)-(Zk-)N---(13)]]>
式(13)中：为tk时刻的待求向量X的先验估计，
X^k-=X^k-1---(14)]]>
Gk是tk时刻的增益矩阵，Gk用于修正tk时刻的待求向量X的先验估计得到tk时刻的待求向量X的后验估计
Gk=PXZSk-1---(15)]]>
和依次表示第1个传感器节点、第2个传感器节点和第N个传感器节点在tk时刻测量的水体离子浓度；
和依次表示预测的tk时刻第1个传感器节点(x1,y1)处、第2个传感器节点(x2,y2)处和第N个传感器节点(xN,yN)处的水体 离子浓度；
式(14)中：为tk-1时刻的待求向量X的后验估计。
式(15)中：
PXZ=W0(λk|k-10-X^k-)(ΔZ0)T+Σp=1nWp(λk|k-1p-X^k-)(ΔZp)T+Σq=n+12nWq(λk|k-1q-X^k-)(ΔZq)T---(16)]]>
Sk＝Rk+PZZ (17)
式(16)：W0、Wp和Wq依次同式(9)、式(10)和式(11)；
和依次同式(6)、式(7)和式(8)；
为tk时刻的待求向量X的先验估计，同式(14)；
n为待求向量X的维数，n＝4；
ΔZ0=Z(x1,y1,tk,λk|k-10)-(Zk-)1Z(x2,y2,tk,λk|k-10)-(Zk-)2···Z(xN,yN,tk,λk|k-10)-(Zk-)N---(18)]]>
ΔZp=Z(x1,y1,tk,λk|k-1p)-(Zk-)1Z(x2,y2,tk,λk|k-1p)-(Zk-)2···Z(xN,yN,tk,λk|k-1p)-(Zk-)N---(19)]]>
ΔZq=Z(x1,y1,tk,λk|k-1q)-(Zk-)1Z(x2,y2,tk,λk|k-1q)-(Zk-)2···Z(xN,yN,tk,λk|k-1q)-(Zk-)N---(20)]]>
式(18)中：同式(6)；
表示当待求向量X为tk时刻的第1个Sigma点时，tk时刻第1个传感器节点(x1,y1)处的理论水体离子浓度；
表示当待求向量X为tk时刻的第1个Sigma点时，tk时刻第2个传感器节点(x2,y2)处的理论水体离子浓度；
表示当待求向量X为tk时刻的第1个Sigma点时，tk时刻第N个传感器节点(xN,yN)处的理论水体离子浓度；
和依次表示预测的tk时刻第1个传感器节点(x1,y1)处、第2个传感器节点(x2,y2)处和第N个传感器节点(xN,yN)处的水体离子浓度。
式(19)中：同式(7)；
表示当待求向量X为tk时刻的第p+1个Sigma点时，tk时刻第1个传感器节点(x1,y1)处的理论水体离子浓度；
表示当待求向量X为tk时刻的第p+1个Sigma点时，tk时刻第2个传感器节点(x2,y2)处的理论水体离子浓度；
表示当待求向量X为tk时刻的第p+1个Sigma点时，tk时刻第N个传感器节点(xN,yN)处的理论水体离子浓度；
和依次表示预测的tk时刻第1个传感器节点(x1,y1)处、第2个传感器节点(x2,y2)处和第N个传感器节点(xN,yN)处的水体离子浓度。
式(20)中：同式(8)；
表示当待求向量X为tk时刻的第q+1个Sigma点时，tk时刻第1个传感器节点(x1,y1)处的理论水体离子浓度；
表示当待求向量X为tk时刻的第q+1个Sigma点时，tk时刻第2个传感器节点(x2,y2)处测量的理论水体离子浓度；
表示当待求向量X为tk时刻的第q+1个Sigma点时，tk时刻第N个传感器节点(xN,yN)处的理论水体离子浓度；
和依次表示预测的tk时刻第1个传感器节点(x1,y1)处、 第2个传感器节点(x2,y2)处和第N个传感器节点(xN,yN)处的水体离子浓度。
式(17)中：Rk是测量噪声的协方差矩阵，所述测量噪声的协方差矩阵为对角矩阵，
Rk=D10···00D2···0·········00···DNN×N---(21)]]>
PZZ=W0(ΔZ0)(ΔZ0)T+Σp=1nWp(ΔZp)(ΔZp)T+Σq=n+12nWq(ΔZq)(ΔZp)T---(22)]]>
式(22)中：W0、Wp和Wq依次同式(9)、式(10)和式(11)；
ΔZ0、ΔZp和ΔZq依次同式(18)、式(19)和式(20)；
n为待求向量X的维数，n＝4。
然后计算误差协方差矩阵
P^k=Pk--GkSk(Gk)T---(23)]]>
式(23)中：Pk-=P^k-1---(24)]]>
Gk和Sk依次同式(15)和式(17)；
式(24)中：为tk-1时刻的误差协方差矩阵。
由于采用上述技术方案，本发明与现有技术相比具有如下经济效果：
本发明无需提前给定阈值，用每个时刻的测量值去更新待求向量X的先验估计得到待求向量X的后验估计，在实际定位操作中简单易行；其次，本发明中测量的浓度数据可以包含整个扩散过程的信息，避免了因忽略过渡阶段影响定位精度的问题；然后，本发明在每个测量时刻均处理了所有传感器节点(xi,yi)(i＝1,...,N)的测量数据，由于在每个测量时刻都考虑到了所有传感器节点(xi,yi)(i＝1,...,N)的测量噪声Vi(k)(i＝1,...,N)，该方法的定位效果相对稳定。
本发明在定位计算过程中只需保存当前时刻的测量值、上一时刻的待求向量X的后验估计和上一时刻的误差协方差矩阵，用当前时刻的测量值去更新待求向量X的先验估计，计算过程中不需要保存和计算历史测量数据。故本发明具有工作效率高和耗时短的优点。
因此，本发明具有易实现、定位精度高、工作效率高和耗时短的特点。
具体实施方式
下面结合具体实施方式对本发明作进一步的描述，并非对?；し段У南拗疲?
实施例1
一种靠近湖泊岸边的污染源定位方法。该方法的具体步骤是：
步骤1、以岸边为纵轴，以垂直于岸边为横轴，在水平面上建立直角坐标系。设污染源在(α,β)处从τ时刻连续释放污染物，污染源的质量流率为Q，以污染源的横坐标α、污染源的纵坐标β、污染源的初始扩散时间τ和污染源的质量流率Q为未知量组成待求向量X，待求向量X＝(α,β,τ,Q)T；传感器节点的总个数为8，传感器节点共测量了5个时刻的水体离子浓度，8个传感器节点(xi,yi)(i＝1,...,8)在tk(k＝1,...,5)时刻测量的水体离子浓度见表1.1。
表1.1
第i(i＝1,...,8)个传感器节点(xi,yi)在tk(k＝1,...,5)时刻测量的水体离子浓度为：
Z~(xi,yitk)=Z(xi,yi,tkX)+Vi(k)---(1)]]>
式(1)中：是第i个传感器节点(xi,yi)在tk时刻的测量噪声，kg/m3；
的均值为0，的方差为Di，(kg/m3)2；8个传感器节点(xi,yi)(i＝1,...,8)在tk(k＝1,...,5)时刻的测量噪声的方差见表1.2；
Z(xi,yi,tk,X)表示tk时刻第i个传感器节点(xi,yi)处的理论水体离子浓度，kg/m3；
Z(xi,yi,tk,X)=ξQ2fπd[1Rierfc(Ri2d(tk-τ))+1R‾ierfc(R‾i2d(tk-τ))]---(2)]]>
式(2)中：ξ为补量纲常量，tunit为时间单位，h；
Ri为第i个传感器节点(xi,yi)到污染源(α,β)的距离，m；
Ri=(xi-α)2+(yi-β)2---(3)]]>
为第i个传感器节点(xi,yi)到镜像污染源(-α,β)的距离，m；
R‾i=(xi+α)2+(yi-β)2---(4)]]>
erfc(x)＝1-erf(x) (5)
f为水体平均深度，f＝10m；
d为同向扩散系数，d＝0.5m2/h。
表1.2
步骤2、设t0＝0h时刻的待求向量X的后验估计为设t0＝0h时刻的误差协方差矩阵为其中：4、4、4和100依次为t0＝0h时刻的污染源的横坐标的估计方差、t0＝0h时刻的污染源的纵坐标的估计方差、t0＝0h时刻的污染源的初始扩散时间的估计方差和t0＝0h时刻的污染源的质量流率的估计方差；令k＝1。
步骤3、将t1＝6h时刻的8个传感器节点测量的水体离子浓度时刻的待求向量X的后验估计和t0＝0h时刻的误差协方差矩阵用动态数据处理?？榻写?。
步骤3.1、由t0＝0h时刻的待求向量X的后验估计和t0＝0h时刻的误差协方差矩阵 得到t1＝6h时刻的第1个Sigma点第p+1个Sigma点和第q+1个Sigma点λ1|0q(q=5,...,8):]]>
λ1|00=X^0---(6)]]>
λ1|0p=X^0+[((4+1)P^0)p]T---(7)]]>
λ1|0q=X^0-[((4+1)P^0)(q-4)]T---(8)]]>
式(7)中：表示开方矩阵的第p行；
式(8)中：表示开方矩阵的第(q-4)行。
t1＝6h时刻的9个Sigma点如下：
λ1|00=(2,8,3,80)T;λ1|01=(6.483302,8,3,80)T;λ1|02=(2,12.4833,3,80)T;λ1|03=(2,8,7.483302,80)T;λ1|04=(2,8,3,102.3629)T;λ1|05=(-2.4833,8,3,80)T;λ1|06=(2,3.516698,3,80)T;λ1|07=(2,8,-1.4833,80)T;λ1|08=(2,8,3,57.63708)T.]]>
计算t1＝6h时刻第1个Sigma点第p+1个Sigma点和第q+1个Sigma点依次对应的权值W0、Wp和Wq：
W0=14+1---(9)]]>
Wp=12(4+1)---(10)]]>
Wq=12(4+1)---(11)]]>
t1＝6h时刻的9个Sigma点依次对应的权值如下：
W0＝0.2；W1＝0.1；W2＝0.1；W3＝0.1；W4＝0.1；W5＝0.1；W6＝0.1；W7＝0.1；W8＝0.1。
步骤3.2、预测t1＝6h时刻各传感器节点(xi,yi)(i＝1,...,8)处的水体离子浓度
(Z1-)i=W0Z(xi,yi,t1,λ1|00)+Σp=14WpZ(xi,yi,t1,λ1|0p)+Σq=58WqZ(xi,yi,t1,λ1|0q)---(12)]]>
式(12)：W0、Wp和Wq依次同式(9)、式(10)和式(11)；
和依次同式(6)、式(7)和式(8)；
表示当待求向量X为t1＝6h时刻的第1个Sigma点时，t1＝6h时刻第i个传感器节点(xi,yi)处的理论水体离子浓度：
当待求向量X为t1＝6h时刻的第1个Sigma点时，t1＝6h时刻各传感器节点处的理论水体离子浓度如下：
Z(x1,y1,t1,λ1|01)=0.304573;Z(x2,y2,t1,λ1|00)=0.169127;Z(x3,y3,t1,λ1|00)=0.141195;Z(x4,y4,t1,λ1|00)=0.497168;Z(x5,y5,t1,λ1|00)=0.3217;Z(x6,y6,t1,λ1|00)=0.722341;Z(x7,y7,t1,λ1|00)=0.34804;Z(x8,y8,t1,λ1|00)=0.162404.]]>
表示当待求向量X为t1＝6h时刻的第p+1个Sigma点时，t1＝6h时刻第i个传感器节点(xi,yi)处的理论水体离子浓度；
当待求向量X为t1＝6h时刻的第2个Sigma点时，t1＝6h时刻各传感器节点处的理论水体离子浓度如下：
Z(x1,y1,t1,λ1|01)=0.003248;Z(x2,y2,t1,λ1|00)=0.002208;Z(x3,y3,t1,λ1|01)=0.001944;Z(x4,y4,t1,λ1|01)=0.006499;Z(x5,y5,t1,λ1|01)=0.014043;Z(x6,y6,t1,λ1|01)=0.025606;Z(x7,y7,t1,λ1|01)=0.008697;Z(x8,y8,t1,λ1|01)=0.005376;]]>
当待求向量X为t1＝6h时刻的第3个Sigma点时，t1＝6h时刻各传感器节点处的理论水体离子浓度如下：
Z(x1,y1,t1,λ1|02)=0.002315;Z(x2,y2,t1,λ1|02)=0.000847;Z(x3,y3,t1,λ1|02)=0.000627;Z(x4,y4,t1,λ1|02)=0.0045;Z(x5,y5,t1,λ1|02)=0.001852;Z(x6,y6,t1,λ1|02)=0.005346;Z(x7,y7,t1,λ1|02)=0.002229;Z(x8,y8,t1,λ1|02)=0.005376.]]>
当待求向量X为t1＝6h时刻的第4个Sigma点时，t1＝6h时刻各传 感器节点处的理论水体离子浓度如下：
Z(x1,y1,t1,λ1|03)=0.107885;Z(x2,y2,t1,λ1|03)=0.042568;Z(x3,y3,t1,λ1|03)=0.031761;Z(x4,y4,t1,λ1|03)=0.245627;Z(x5,y5,t1,λ1|03)=0.14569;Z(x6,y6,t1,λ1|03)=0.467618;Z(x7,y7,t1,λ1|03)=0.154867;Z(x8,y8,t1,λ1|03)=0.04712.]]>
当待求向量X为t1＝6h时刻的第5个Sigma点时，t1＝6h时刻各传感器节点处的理论水体离子浓度如下：
Z(x1,y1,t1,λ1|04)=0.38984;Z(x2,y2,t1,λ1|04)=0.216404;Z(x3,y3,t1,λ1|04)=0.180664;Z(x4,y4,t1,λ1|04)=0.636144;Z(x5,y5,t1,λ1|04)=0.411627;Z(x6,y6,t1,λ1|04)=0.924261;Z(x7,y7,t1,λ1|04)=0.44533;Z(x8,y8,t1,λ1|04)=0.207802.]]>
表示当待求向量X为t1＝6h时刻的第q+1个Sigma点时，t1＝6h时刻第i个传感器节点(xi,yi)处的理论水体离子浓度；
当待求向量X为t1＝6h时刻的第6个Sigma点时，t1＝6h时刻各传感器节点处的理论水体离子浓度如下：
Z(x1,y1,t1,λ1|05)=0.205871;Z(x2,y2,t1,λ1|05)=0.124709;Z(x3,y3,t1,λ1|05)=0.105074;Z(x4,y4,t1,λ1|05)=0.357685;Z(x5,y5,t1,λ1|05)=0.288039;Z(x6,y6,t1,λ1|05)=0.665988;Z(x7,y7,t1,λ1|05)=0.281433;Z(x8,y8,t1,λ1|05)=0.135239.]]>
当待求向量X为t1＝6h时刻的第7个Sigma点时，t1＝6h时刻各传感器节点处的理论水体离子浓度如下：
Z(x1,y1,t1,λ1|06)=0.360151;Z(x2,y2,t1,λ1|06)=0.581824;Z(x3,y3,t1,λ1|06)=0.653197;Z(x4,y4,t1,λ1|06)=0.254905;]]>
Z(x5,y5,t1,λ1|06)=0.487653;Z(x6,y6,t1,λ1|06)=0.216877;Z(x7,y7,t1,λ1|05)=0.426542;Z(x8,y8,t1,λ1|06)=0.917164.]]>
当待求向量X为t1＝6h时刻的第8个Sigma点时，t1＝6h时刻各传感器节点处的理论水体离子浓度如下：
Z(x1,y1,t1,λ1|07)=0.62234;Z(x2,y2,t1,λ1|07)=0.427979;Z(x3,y3,t1,λ1|07)=0.383123;Z(x4,y4,t1,λ1|07)=0.846294;Z(x5,y5,t1,λ1|07)=0.601173;Z(x6,y6,t1,λ1|07)=1.048705;Z(x7,y7,t1,λ1|07)=0.64922;Z(x8,y8,t1,λ1|07)=0.399828.]]>
当待求向量X为t1＝6h时刻的第9个Sigma点时，t1＝6h时刻各传感器节点处的理论水体离子浓度如下：
Z(x1,y1,t1,λ1|08)=0.219506;Z(x2,y2,t1,λ1|08)=0.12185;Z(x3,y3,t1,λ1|08)=0.101726;Z(x4,y4,t1,λ1|08)=0.358191;Z(x5,y5,t1,λ1|08)=0.231773;Z(x6,y6,t1,λ1|08)=0.52042;Z(x7,y7,t1,λ1|08)=0.25075;Z(x8,y8,t1,λ1|08)=0.117006.]]>
预测的t1＝6h时刻各传感器节点处的水体离子浓度如下：
(Z1-)0=W0Z(x1,y1,t1,λ1|00)+Σp=14WpZ(x1,y1,t1,λ1|0p)+Σq=58WqZ(x1,y1,t1,λ1|0q)=0.25305;(Z1-)2=W0Z(x2,y2,t1,λ1|00)+Σp=14WpZ(x2,y2,t1,λ1|0p)+Σq=58WqZ(x2,y2,t1,λ1|0q)=0.185664;(Z1-)3=W0Z(x3,y3,t1,λ1|00)+Σp=14WpZ(x3,y3,t1,λ1|0p)+Σq=58WqZ(x3,y3,t1,λ1|0q)=0.174051;(Z1-)4=W0Z(x4,y4,t1,λ1|00)+Σp=14WpZ(x4,y4,t1,λ1|0p)+Σq=58WqZ(x4,y4,t1,λ1|0q)=0.370418;(Z1-)5=W0Z(x5,y5,t1,λ1|00)+Σp=14WpZ(x5,y5,t1,λ1|0p)+Σq=58WqZ(x5,y5,t1,λ1|0q)=0.282525;]]>
(Z1-)6=W0Z(x6,y6,t1,λ1|00)+Σp=14WpZ(x6,y6,t1,λ1|0p)+Σq=58WqZ(x6,y6,t1,λ1|0q)=0.53195;(Z1-)7=W0Z(x7,y7,t1,λ1|00)+Σp=14WpZ(x7,y7,t1,λ1|0p)+Σq=58WqZ(x7,y7,t1,λ1|0q)=0.291515;(Z1-)8=W0Z(x8,y8,t1,λ1|00)+Σp=14WpZ(x8,y8,t1,λ1|0p)+Σq=58WqZ(x8,y8,t1,λ1|0q)=0.215504.]]>
步骤3.3、计算t1＝6h时刻的待求向量X的后验估计和t1＝6h时刻的误差协方差矩阵首先计算t1＝6h时刻的待求向量X的后验估计
X^1=X^1-+G1Z~(x1,y1,t1)-(Z1-)1Z~(x2,y2,t1)-(Z1-)2···Z~(x8,y8,t1)-(Z1-)8---(13)]]>
式(13)中：为t1＝6h时刻的待求向量X的先验估计，
X^1-=X^0---(14)]]>
X^1-=(2,8,3,80)T;]]>
G1是t1＝6h时刻的增益矩阵，G1用于修正t1＝6h时刻的待求向量X的先验估计得到t1＝6h时刻的待求向量X的后验估计
G1=PXZS1-1---(15)]]>
和依次表示第1个传感器节点、第2个传感器节点和第8个传感器节点在t1＝6h时刻测量的水体离子浓度；
和依次表示预测的t1＝6h时刻第1个传感器节点(x1,y1)处、第2个传感器节点(x2,y2)处和第8个传感器节点(x8,y8)处的水体离子浓度。
式(14)中：为t0＝0h时刻的待求向量X的后验估计。
式(15)中：
PXZ=W0(λ1|00-X^1-)(ΔZ0)T+Σp=1nWp(λ1|0p-X^1-)(ΔZp)T+Σq=58Wq(λ1|0q-X^1-)(ΔZq)T---(16)]]>
S1＝R1+PZZ (17)式(16)：W0、Wp和Wq依次同式(9)、式(10)和式(11)；
和依次同式(6)、式(7)和式(8)；
为t1＝6h时刻的待求向量X的先验估计，同式(14)；
ΔZ0=Z(x1,y1,t1,λ1|00)-(Z1-)1Z(x2,y2,t1,λ1|00)-(Z1-)2···Z(x8,y8,t1,λ1|00)-(Z1-)8---(18)]]>
ΔZ0＝(0.051623,-0.01654,-0.03286,0.12675,0.039175,0.19039,0.056525,-0.0531)T。
ΔZp=Z(x1,y1,t1,λ1|0p)-(Z1-)1Z(x2,y2,t1,λ1|0p)-(Z1-)2···Z(x8,y8,t1,λ1|0p)-(Z1-)8---(19)]]>
ΔZ1＝(-0.2498,-0.18346,-0.17211,-0.36392,-0.26848,-0.50634,-0.28282,-0.21013)T；
ΔZ2＝(-0.25074,-0.18482,-0.17342,-0.36592,-0.28067,-0.5266,-0.28929,-0.2148)T；
ΔZ3＝(-0.14517,-0.14309,-0.14229,-0.12479,-0.13684,-0.06433,-0.13665,-0.16838)T；
ΔZ4＝(0.13679,0.03074,0.006614,0.265726,0.129102,0.392311,0.153815,-0.0077)T。
ΔZq=Z(x1,y1,t1,λ1|0q)-(Z1-)1Z(x2,y2,t1,λ1|0q)-(Z1-)2···Z(x8,y8,t1,λ1|0q)-(Z1-)8---(20)]]>
ΔZ5＝(-0.03718,-0.06095,-0.06898,-0.01273,0.005514,0.134037,-0.01008,-0.08027)T；
ΔZ6＝(0.1071,0.396159,0.479146,-0.11551,0.205128,-0.31507,0.135028,0.70166)T；
ΔZ7＝(0.36929,0.242315,0.209072,0.475876,0.318648,0.516755,0.357705,0.184324)T；
ΔZ8＝(-0.03354,-0.06381,-0.07232,-0.01223,-0.05075,-0.01153,-0.04076,-0.0985)T。
式(18)中：同式(6)；
表示当待求向量X为t1＝6h时刻的第1个Sigma点时，t1＝6h时刻第1个传感器节点(x1,y1)处的理论水体离子浓度；
表示当待求向量X为t1＝6h时刻的第1个Sigma点时，t1＝6h时刻第2个传感器节点(x2,y2)处的理论水体离子浓度；
表示当待求向量X为t1＝6h时刻的第1个Sigma点时，t1＝6h时刻第8个传感器节点(x8,y8)处的理论水体离子浓度；
和依次表示预测的t1＝6h时刻第1个传感器节点(x1,y1)处、第2个传感器节点(x2,y2)处和第8个传感器节点(x8,y8)处的水体离子浓度。
式(19)中：同式(7)；
表示当待求向量X为t1＝6h时刻的第p+1个Sigma点时，t1＝6h时刻第1个传感器节点(x1,y1)处的理论水体离子浓度；
表示当待求向量X为t1＝6h时刻的第p+1个Sigma点时，t1＝6h时刻第2个传感器节点(x2,y2)处的理论水体离子浓度；
表示当待求向量X为t1＝6h时刻的第p+1个Sigma点时，t1＝6h时刻第8个传感器节点(x8,y8)处的理论水体离子浓度；
和依次表示预测的t1＝6h时刻第1个传感器节点(x1,y1)处、第2个传感器节点(x2,y2)处和第8个传感器节点(x8,y8)处的水体离子浓度。
式(20)中：同式(8)；
表示当待求向量X为t1＝6h时刻的第q+1个Sigma点时， t1＝6h时刻第1个传感器节点(x1,y1)处的理论水体离子浓度；
表示当待求向量X为t1＝6h时刻的第q+1个Sigma点时，t1＝6h时刻第2个传感器节点(x2,y2)处测量的理论水体离子浓度；
表示当待求向量X为t1＝6h时刻的第q+1个Sigma点时，t1＝6h时刻第8个传感器节点(x8,y8)处的理论水体离子浓度；
和依次表示预测的t1＝6h时刻第1个传感器节点(x1,y1)处、第2个传感器节点(x2,y2)处和第8个传感器节点(x8,y8)处的水体离子浓度。
式(17)中：R1是测量噪声的协方差矩阵，所述测量噪声的协方差矩阵为对角矩阵，
R1=10···001···0·········00···18×8---(21)]]>
PZZ=W0(ΔZ0)(ΔZ0)T+Σp=1nWp(ΔZp)(ΔZp)T+Σq=58Wq(ΔZq)(ΔZp)T---(22)]]>
式(22)中：W0、Wp和Wq依次同式(9)、式(10)和式(11)；
ΔZ0、ΔZp和ΔZq依次同式(18)、式(19)和式(20)。
PXZ=-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.158130.3809180.2114510.1765290.6215840.4022050.9031060.4351370.203046;]]>
PZZ=0.0320740.0251760.0238160.0414450.0320150.0493680.033820.0273760.0251760.0313230.0334580.0227330.0284760.0198150.0271150.0437680.0238160.0334580.0365430.0183180.0280830.0126370.0258720.0488210.0414450.0227330.0183180.0624760.0390210.0818240.0436280.0169470.0320150.0284760.0280830.0390210.0335540.0459470.034380.0341810.0493680.0198150.0126370.0818240.0459470.1148650.0527620.0071660.033820.0271150.0258720.0436280.034380.0527620.0360340.0302890.0273760.0437680.0488210.0169470.0341810.0071660.0302890.066679;]]>
S1=0.1320740.0251760.0238160.0414450.0320150.0493680.033820.0273760.0251760.1313230.0334580.0227330.0284760.0198150.0271150.0437680.0238160.0334580.1365430.0183180.0280830.0126370.0258720.0488210.0414450.0227330.0183180.1624760.0390210.0818240.0436280.0169470.0320150.0284760.0280830.0390210.1335540.0459470.034380.0341810.0493680.0198150.0126370.0818240.0459470.2148650.0527620.0071660.033820.0271150.0258720.0436280.034380.0527620.1360340.0302890.0273760.0437680.0488210.016470.0341810.0071660.0302890.166679;]]>
G1=-0.05306-0.02173-0.02051-0.23301-0.35318-1.08792-0.26062-0.13862-0.26911-0.88585-1.076560.04595-0.639150.051271-0.45221-1.66658-0.83763-0.52728-0.42212-0.81684-0.54311-0.35225-0.6505-0.221250.6577620.1985360.130951.5920450.9543522.9602670.9918230.354588;]]>
X^1=(1.783544,8.199684,3.343394,80.13257)T.]]>
然后计算误差协方差矩阵
P^1=P1--G1S1(G1)T---(23)]]>
式(23)中：P1-P^0---(24)]]>
P1-=diag(4,4,4,100);]]>
G1和S1依次同式(15)和式(17)；
式(24)中：为t0＝0h时刻的误差协方差矩阵；
P^1=3.560445-0.33296-0.51721.439375-0.332962.511154-0.853271.197235-0.5172-0.853273.0473941.8773541.4393751.1972351.87735495.13391.]]>
故：t1＝6h时刻的待求向量X的后验估计X^1=(1.783544,8.199684,3.343394,80.13257)T;]]>
t1＝6h时刻的误差协方差矩阵P^1=3.560445-0.33296-0.51721.439375-0.332962.511154-0.853271.197235-0.5172-0.853273.0473941.8773541.4393751.1972351.87735495.13391.]]>
步骤4、k＝2，同理步骤3，得到t2＝10h时刻：
待求向量X的后验估计X^2=(1.683838,8.038214,3.561119,82.43101)T;]]>
误差协方差矩阵P^2=3.270747-0.60335-0.256651.583841-0.603350.85556-0.402241.233375-0.25665-0.402242.7469541.7492591.5838411.2333751.74925994.46743.]]>
……，
k＝5，同理步骤3，得到t5＝18h时刻：
待求向量X的后验估计X^5=(2.559427,7.21534,4.065863,82.51426)T;]]>
误差协方差矩阵P^5=0.744604-0.231720.0542461.655629-0.231720.213794-0.17651.1635220.054246-0.17652.3972132.1870481.6556291.1635222.18704892.45414.]]>
步骤5、在t5＝18h时刻的待求向量X的后验估计X^5=(2.559427,7.21534,4.065863,82.51426)T]]>中，2.559427为所求的污染源的横坐标，7.21534为所求的污染源的纵坐标。
实施例2
一种靠近湖泊岸边的污染源定位方法。该方法的具体步骤是：
步骤1、以岸边为纵轴，以垂直于岸边为横轴，在水平面上建立直角坐标系。设污染源在(α,β)处从τ时刻连续释放污染物，污染源的质量流率为Q，以污染源的横坐标α、污染源的纵坐标β、污染源的初始扩散时间τ和污染源的质量流率Q为未知量组成待求向量X，待求向量X＝(α,β,τ,Q)T；传感器节点的总个数为8，传感器节点共测量了9个时刻的水体离子浓度，8个传感器节点(xi,yi)(i＝1,...,8)在tk(k＝1,...,9)时刻测量的水体离子浓度见表2.1。
表2.1
第i(i＝1,...,8)个传感器节点(xi,yi)在tk(k＝1,...,9)时刻测量的水体离子浓度为：
Z~(xi,yitk)=Z(xi,yi,tkX)+Vi(k)---(1)]]>
式(1)中：是第i个传感器节点(xi,yi)在tk时刻的测量噪声，kg/m3；
的均值为0，的方差为Di，(kg/m3)2；8个传感器节(xi,yi)(i＝1,...,8)在tk(k＝1,...,9)时刻的测量噪声的方差见表2.2；
Z(xi,yi,tk,X)表示tk时刻第i个传感器节点(xi,yi)处的理论水体离子浓度，kg/m3；
Z(xi,yi,tk,X)=ξQ2fπd[1Rierfc(Ri2d(tk-τ))+1R‾ierfc(R‾i2d(tk-τ))]---(2)]]>
式(2)中：ξ为补量纲常量，tunit为时间单位，h；
Ri为第i个传感器节点(xi,yi)到污染源(α,β)的距离，m；
Ri=(xi-α)2+(yi-β)2---(3)]]>
为第i个传感器节点(xi,yi)到镜像污染源(-α,β)的距离，m；
R‾i=(xi+α)2+(yi-β)2---(4)]]>
erfc(x)=1-erf(x) (5)
f为水体平均深度，f=10m；
d为同向扩散系数，d=0.5m2/h
表2.2
步骤2、同实施例1的步骤2。
步骤3、同实施例1的步骤3。
步骤4、k＝2，同理步骤3，得到t2＝10h时刻：
待求向量X的后验估计X^2=(1.683838,8.038214,3.561119,82.43101)T;]]>
误差协方差矩阵P^2=3.270747-0.60335-0.256651.583841-0.603350.85556-0.402241.233375-0.25665-0.402242.7469541.7492591.5838411.2333751.74925994.46743.]]>
……。
k＝9，同理步骤3，得到t9＝30h时刻：
待求向量X的后验估计X^9=(2.93044,6.062932,5.648708,85.69018)T;]]>
误差协方差矩阵P^9=-0.00139-0.007850.033080.70294-0.007850.044006-0.037050.7027720.03308-0.037052.1590883.1785410.702940.7027723.17854160.2314.]]>
步骤5、在t9＝30h时刻的待求向量X的后验估计X^9=(2.93044,6.062932,5.648708,85.69018)T;]]>中，2.93044为所求的污染源的横坐标，6.062932为所求的污染源的纵坐标。
本具体实施方式与现有技术相比具有如下经济效果：
本具体实施方式无需提前给定阈值，用每个时刻的测量值去更新待求向量X的先验估计得到待求向量X的后验估计，在实际定位操作中简单易行；其次，本具体实施方式中测量的浓度数据可以包含整个扩散过程的信息，避免了因忽略过渡阶段影响定位精度的问题；然后，本具体实施方式在每个测量时刻均处理了所有传感器节点(xi,yi)(i＝1,...,N)的测量数据，由于在每个测量时刻都考虑到了所有传感器节点(xi,yi)(i＝1,...,N)的测量噪声该方法的定位效果相对稳定。
本具体实施方式在定位计算过程中只需保存当前时刻的测量值、上一时刻的待求向量X的后验估计和上一时刻的误差协方差矩阵，用当前时刻的测量值去更新待求向量X的先验估计，计算过程中不需要保存和计算历史测量数据。故本具体实施方式具有工作效率高和耗时 短的优点。
因此，本具体实施方式具有易实现、定位精度高、工作效率高和耗时短的特点。