[发明专利]一种接近湖泊岸边的三维污染源定位方法

摘要

本发明涉及一种接近湖泊岸边的三维污染源定位方法。其技术方案为将tk(k＝1,...,S)时刻测量的水体污染物浓度值tk‑1时刻待求的估计向量X的后验估计和tk‑1时刻的误差协方差矩阵用无迹卡尔曼滤波数据处理方法进行数据处理，获取tk时刻的待求向量X的后验估计和tk时刻的误差协方差矩阵重复上述方法直至得到k＝S时tS时刻的待求向量X的后验估计即得污染源的横坐标、纵坐标和深度坐标依次为和本发明具有定位周期短、精度高且稳定性强的特点。

主权项

一种接近湖泊岸边的三维污染源定位方法，其特征在于所述三维污染源定位方法的步骤为：步骤1、在接近湖泊岸边的三维水体中，以岸边方向为x轴、垂直于岸边方向为y轴、水体深度为z轴和岸边水体底部任意一点为原点，建立三维直角坐标系；假设在接近湖泊岸边(ε,η,ξ)位置处有一连续污染源从τ时刻开始以速率为Q均匀地向三维水体中排放浓度为C0的污染物，所述污染物的质量流率M＝QC0；三维水体中设有N个传感器节点，N个传感器节点共计监测了S个采样时刻的水体污染物浓度，第i(i＝1,...,N)个传感器节点(xi,yi,zi)在tk(k＝1,...,S)时刻监测的水体污染物浓度为：C~(xi,yi,zi,tk)=C(xi,yi,zi,tk,X)+Wi(k)---(1)]]>式(1)中：tk(k＝1,...,S)表示第k个采样时刻，h，X表示以污染源位置(ε,η,ξ)、初始扩散时间τ和所述污染物的质量流率M所组成的待求向量,X＝(ε,η,ξ,τ,M)T，Wi(k)表示第i(i＝1,...,N)个传感器节点(xi,yi,zi)在tk(k＝1,...,S)时刻的监测噪声，kg/m3；Wi(k)的均值为0，Wi(k)的方差为Vi(k)，(kg/m3)2，C(xi,yi,zi,tk,X)表示第i(i＝1,...,N)个传感器节点(xi,yi,zi)在tk时刻的理论水体污染物浓度，kg/m3，C(xi,yi,zi,tk,X)=M2πD[1Rierfc(Ri2D(tk-τ))+1R‾ierfc(R‾i2D(tk-τ))]---(2)]]>式(2)中：Ri表示第i(i＝1,...,N)个传感器节点(xi,yi,zi)到污染源(ε,η,ξ)的距离，m，Ri=(xi-ϵ)2+(yi-η)2+(zi-ξ)2---(3)]]>Ri表示第i(i＝1,...,N)个传感器节点(xi,yi,zi)到镜像污染源(ε,‑η,ξ)的距离，m，R‾i=(xi-ϵ)2+(yi+η)2+(zi-ξ)2---(4)]]>D表示同向扩散系数，m2/h，M表示所述污染物的质量流率，kg/h，tk(k＝1,...,S)表示第k个采样时刻，h，τ表示初始扩散时间，h，erfc(x)表示余补误差函数，erfc(x)＝1‑erf(x)                                        (5)式(5)中：erf(x)表示误差函数，其值可以从专业文献中查表获取；步骤2、初始化t0时刻待求向量X的后验估计为X^0=(ϵ^(0),η^(0),ξ^(0),τ^(0),M^(0))T---(6)]]>式(6)中：和依次表示t0时刻污染源的估计横坐标、污染源的估计纵坐标和污染源的估计深度坐标，m，表示t0时刻污染源的估计初始扩散时间，h，表示t0时刻污染源的估计质量流率，kg/h，初始化t0时刻误差协方差矩阵为P^0=diag(σ^ϵ(0)2,σ^η(0)2,σ^ξ(0)2,σ^τ(0)2,σ^M(0)2)---(7)]]>式(7)中：和依次表示t0时刻污染源的横坐标估计方差、污染源的纵坐标估计方差和污染源的深度坐标估计方差，m2，表示t0时刻污染源的初始扩散时间τ估计方差，h2，表示t0时刻污染源的质量流率M估计方差，(kg/h)2；将采样时刻tk置为第1个采样时刻t1，即令k＝1；步骤3、将tk时刻的N个传感器节点监测的水体污染物浓度tk‑1时刻的待求向量X的后验估计和tk‑1时刻的误差协方差矩阵用无迹卡尔曼滤波数据处理方法进行处理；用所述无迹卡尔曼滤波数据处理方法的处理过程如下：步骤3.1、根据tk‑1时刻的待求向量X的后验估计得到tk时刻待求向量X的第1个Sigma粒子χk|k-10=X^k-1---(8)]]>再根据tk‑1时刻的误差协方差矩阵得到tk时刻第p+1个Sigma粒子和第q+1个Sigma粒子χk|k-1p=X^k-1+[((n+λ)P^k-1)p]T---(9)]]>χk|k-1q=X^k-1-[((n+λ)P^k-1)q]T---(10)]]>式(9)、(10)中：λ表示缩放比例参数，λ＝1，n表示待求向量X的维数，n＝5，p表示Sigma粒子的个数，p＝1,...,n，q表示Sigma粒子的个数，q＝n+1,...,2n，表示开方矩阵的第p行，表示开方矩阵的第q行；计算2n+1个Sigma粒子的权值W0=λn+λ---(11)]]>Wp=12(n+λ)---(12)]]>Wq=12(n+λ)---(13)]]>式(11)、(12)、(13)中：λ表示缩放比例参数，λ＝1，n表示待求向量X的维数，n＝5，p表示Sigma粒子的个数，p＝1,...,n，q表示Sigma粒子的个数，q＝n+1,...,2n；步骤3.2、预测tk时刻各传感器节点(xi,yi,zi)(i＝1,...,N)处的水体污染物浓度(Ck-)i=W0C(xi,yi,zi,tk,χk|k-10)+Σp=1nWpC(xi,yi,zi,tk,χk|k-1p)+Σq=n+12nWqC(xi,yi,zi,tk,χk|k-1q)---(14)]]>式(14)中：W0表示第1个sigma粒子的权值，同式(11)，Wp表示第p+1个sigma粒子的权值，同式(12)，Wq表示第q+1个sigma粒子的权值，同式(13)，n为待求向量X的维数，n＝5，表示待求向量X在tk时刻的第1个sigma粒子时传感器节点(xi,yi,zi)的理论水体污染物浓度，kg/m3，表示待求向量X在tk时刻的第p+1(p＝1,...,n)个sigma粒子时传感器节点(xi,yi,zi)的理论水体污染物浓度，kg/m3，表示待求向量X在tk时刻的第q+1(q＝n+1,...,2n)个sigma粒子时传感器节点(xi,yi,zi)的理论水体污染物浓度，kg/m3；步骤3.3、根据tk时刻N个传感器节点监测的水体污染物浓度更新待求向量的后验估计X^k=Xk-+KkC~(x1,y1,z1,tk)-(Ck-)1C~(x2,y2,z2,tk)-(Ck-)2...C~(xN,yN,zN,tk)-(Ck-)N---(15)]]>式(15)中：表示tk时刻的待求向量X的先验估计，表示tk‑1时刻的待求向量X的后验估计，表示tk时刻第1个传感器节点处预测的水体污染物浓度，同式(14)，kg/m3，表示tk时刻第2个传感器节点处预测的水体污染物浓度，同式(14)，kg/m3，表示tk时刻第N个传感器节点处预测的水体污染物浓度，同式(14)，kg/m3，表示tk时刻第1个传感器节点处监测的水体污染物浓度，同式(1)，kg/m3，表示tk时刻第2个传感器节点处监测的水体污染物浓度，同式(1)，kg/m3，分别表示tk时刻第N个传感器节点处监测的水体污染物浓度，同式(1)，kg/m3，Kk表示tk时刻的卡尔曼增益矩阵，Kk=PXCSk-1---(16)]]>式(16)中：PXC表示tk时刻的估计最小均方误差协方差矩阵，PXC=W0(χk|k-10-Xk-)Q0T+Σp=1nWp(χk|k-1p-X^k-)QpT+Σq=n+12nWq(χk|k-1q-X^k-)QqT---(17)]]>Sk表示tk时刻测量噪声协方差与估计测量噪声最小均方误差协方差之和组成的测量总误差矩阵，Sk＝PCC+Rk                                          (18)式(17)、(18)中：p表示Sigma粒子的个数，p＝1,...,n，q表示Sigma粒子的个数，q＝n+1,...,2n，n表示待求向量X的维数，n＝5，W0表示第1个sigma粒子的权值，同式(11)，Wp表示第p+1个sigma粒子的权值，同式(12)，Wq表示第q+1个sigma粒子的权值，同式(13)，表示第1个sigma粒子，同式(8)，表示第p+1个sigma粒子，同式(9)，表示第q+1个sigma粒子，同式(10)，表示tk时刻的待求向量X的先验估计，表示tk‑1时刻的待求向量X的后验估计，Q0表示待求向量X在tk时刻的第1个sigma粒子时N个传感器节点理论水体污染物浓度与监测的水体污染物浓度的差值所成的偏差矩阵，Q0=C(x1,y1,z1,tk,χk|k-1(0))-(Ck-)1C(x2,y2,z2,tk,χk|k-1(0))-(Ck-)2...C(xN,yN,zN,tk,χk|k-1(0))-(Ck-)N---(19)]]>Qp表示待求向量X在tk时刻的第p+1个sigma粒子时N个传感器节点理论水体污染物浓度与监测的水体污染物浓度的差值所成的偏差矩阵，Qp=C(x1,y1,z1,tk,χk|k-1(p))-(Ck-)1C(x2,y2,z2,tk,χk|k-1(p))-(Ck-)2...C(xN,yN,zN,tk,χk|k-1(p))-(Ck-)N---(20)]]>Qq表示待求向量X在tk时刻的第q+1个sigma粒子时N个传感器节点理论水体污染物浓度与监测的水体污染物浓度的差值所成的偏差矩阵，Qq=C(x1,y1,z1,tk,χk|k-1(q))-(Ck-)1C(x2,y2,z2,tk,χk|k-1(q))-(Ck-)2...C(xN,yN,zN,tk,χk|k-1(q))-(Ck-)N---(21)]]>PCC表示tk时刻的估计测量噪声最小均方误差协方差矩阵，PCC=W0Q0Q0T+Σp=1nWpQpQpT+Σq=n+12nWqQqQqT---(22)]]>Rk表示tk时刻的测量噪声协方差矩阵，为N×N的对角矩阵，Rk=V1(k)0...00V2(k)...0.........00...VN(k)N×N---(23)]]>式(19)、(20)、(21)中：表示第1个sigma粒子，同式(8)，表示第p+1个sigma粒子，同式(9)，表示第q+1个sigma粒子，同式(10)，表示待求向量X在tk时刻的第1个粒子时第1个传感器节点(x1,y1,z1)的理论水体污染物浓度，kg/m3，表示待求向量X在tk时刻的第1个粒子时第2个传感器节点(x2,y2,z2)的理论水体污染物浓度，kg/m3，表示待求向量X在tk时刻的第1个粒子时第N个传感器节点(xN,yN,zN)的理论水体污染物浓度，kg/m3，表示待求向量X在tk时刻的第p+1个粒子时第1个传感器节点(x1,y1,z1)的理论水体污染物浓度，kg/m3，表示待求向量X在tk时刻的第p+1个粒子时第2个传感器节点(x2,y2,z2)的理论水体污染物浓度，kg/m3，表示待求向量X在tk时刻的第p+1个粒子时第N个传感器节点(xN,yN,zN)的理论水体污染物浓度，kg/m3，表示待求向量X在tk时刻的第q+1个粒子时第1个传感器节点(x1,y1,z1)的理论水体污染物浓度，kg/m3，表示待求向量X在tk时刻的第q+1个粒子时第2个传感器节点(x2,y2,z2)的理论水体污染物浓度，kg/m3，表示待求向量X在tk时刻的第q+1个粒子时第N个传感器节点(xN,yN,zN)的理论水体污染物浓度，kg/m3，表示tk时刻第1个传感器节点处预测的水体污染物浓度，同式(14)，kg/m3，表示tk时刻第2个传感器节点处预测的水体污染物浓度，同式(14)，kg/m3，表示tk时刻第N个传感器节点处预测的水体污染物浓度，同式(14)，kg/m3；式(22)中：W0表示第1个sigma粒子的权值，同式(11)，Wp表示第p+1个sigma粒子的权值，同式(12)，Wq表示第q+1个sigma粒子的权值，同式(13)；式(23)中：表示第1个传感器节点在tk(k＝1,...,S)时刻的监测噪声方差，(kg/m3)2，表示第2个传感器节点在tk(k＝1,...,S)时刻的监测噪声方差，(kg/m3)2，表示第N个传感器节点在tk(k＝1,...,S)时刻的监测噪声方差，(kg/m3)2；步骤3.4、计算tk时刻的后验估计误差协方差矩阵P^=P^k--KkSkKkT---(24)]]>式(24)中：表示tk时刻先验估计误差协方差矩阵，Kk表示tk时刻的卡尔曼增益矩阵，同式(16)，Sk表示tk时刻测量噪声协方差与估计测量噪声最小均方误差协方差之和组成的测量总误差矩阵，同式(18)；得到tk时刻的待求向量X的后验估计和tk时刻的误差协方差矩阵步骤4、将采样时刻tk置为第k+1个采样时刻tk+1，即令k＝k+1，重复步骤3，直至得到k＝S时的tS时刻的待求向量X的后验估计和误差协方差矩阵所述tS时刻的待求向量X的后验估计中的和依次为所求的污染源的横坐标、纵坐标和深度坐标。