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

权利要求书

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

步骤1、在接近湖泊岸边的三维水体中，以岸边方向为x轴、垂直于岸边方向为y轴、水体深度为z轴和岸边水体底部任意一点为原点，建立三维直角坐标系；假设在接近湖泊岸边(ε,η,ξ)位置处有一连续污染源从τ时刻开始以速率为Q均匀地向三维水体中排放浓度为C₀的污染物，所述污染物的质量流率M＝QC₀；三维水体中设有N个传感器节点，N个传感器节点共计监测了S个采样时刻的水体污染物浓度，第i(i＝1,...,N)个传感器节点(x_i,y_i,z_i)在t_k(k＝1,...,S)时刻监测的水体污染物浓度为：

C~(xi,yi,zi,tk)=C(xi,yi,zi,tk,X)+Wi(k)---(1)]]>

式(1)中：

t_k(k＝1,...,S)表示第k个采样时刻，h，

X表示以污染源位置(ε,η,ξ)、初始扩散时间τ和所述污染物的质量流率M所组成的待求向量,X＝(ε,η,ξ,τ,M)^T，

W_i^(k)表示第i(i＝1,...,N)个传感器节点(x_i,y_i,z_i)在t_k(k＝1,...,S)时刻的监测噪声，kg/m³；

W_i^(k)的均值为0，W_i^(k)的方差为V_i^(k)，(kg/m³)²，

C(x_i,y_i,z_i,t_k,X)表示第i(i＝1,...,N)个传感器节点(x_i,y_i,z_i)在t_k时刻的理论水体污染物浓度，kg/m³，

C(xi,yi,zi,tk,X)=M2πD[1Rierfc(Ri2D(tk-τ))+1R‾ierfc(R‾i2D(tk-τ))]---(2)]]>

式(2)中：

R_i表示第i(i＝1,...,N)个传感器节点(x_i,y_i,z_i)到污染源(ε,η,ξ)的距离，m，

Ri=(xi-ϵ)2+(yi-η)2+(zi-ξ)2---(3)]]>

R_i表示第i(i＝1,...,N)个传感器节点(x_i,y_i,z_i)到镜像污染源(ε,-η,ξ)的距离，m，

R‾i=(xi-ϵ)2+(yi+η)2+(zi-ξ)2---(4)]]>

D表示同向扩散系数，m²/h，

M表示所述污染物的质量流率，kg/h，

t_k(k＝1,...,S)表示第k个采样时刻，h，

τ表示初始扩散时间，h，

erfc(x)表示余补误差函数，

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

式(5)中：

erf(x)表示误差函数，其值可以从专业文献中查表获取；

步骤2、初始化t₀时刻待求向量X的后验估计为

X^0=(ϵ^(0),η^(0),ξ^(0),τ^(0),M^(0))T---(6)]]>

式(6)中：

和依次表示t₀时刻污染源的估计横坐标、污染源的估计纵坐标和污染源的估计深度坐标，m，

表示t₀时刻污染源的估计初始扩散时间，h，

表示t₀时刻污染源的估计质量流率，kg/h，

初始化t₀时刻误差协方差矩阵为

P^0=diag(σ^ϵ(0)2,σ^η(0)2,σ^ξ(0)2,σ^τ(0)2,σ^M(0)2)---(7)]]>

式(7)中：

和依次表示t₀时刻污染源的横坐标估计方差、污染源的纵坐标估计方差和污染源的深度坐标估计方差，m²，

表示t₀时刻污染源的初始扩散时间τ估计方差，h²，

表示t₀时刻污染源的质量流率M估计方差，(kg/h)²；

将采样时刻t_k置为第1个采样时刻t₁，即令k＝1；

步骤3、将t_k时刻的N个传感器节点监测的水体污染物浓度t_k-1时刻的待求向量X的后验估计和t_k-1时刻的误差协方差矩阵用无迹卡尔曼滤波数据处理方法进行处理；

用所述无迹卡尔曼滤波数据处理方法的处理过程如下：

步骤3.1、根据t_k-1时刻的待求向量X的后验估计得到t_k时刻待求向量X的第1个Sigma粒子

χk|k-10=X^k-1---(8)]]>

再根据t_k-1时刻的误差协方差矩阵得到t_k时刻第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)]]> 2

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、预测t_k时刻各传感器节点(x_i,y_i,z_i)(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)中：

W₀表示第1个sigma粒子的权值，同式(11)，

W_p表示第p+1个sigma粒子的权值，同式(12)，

W_q表示第q+1个sigma粒子的权值，同式(13)，

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

表示待求向量X在t_k时刻的第1个sigma粒子时传感器节点(x_i,y_i,z_i)的理论水体污染物浓度，kg/m³，

表示待求向量X在t_k时刻的第p+1(p＝1,...,n)个sigma粒子时传感器节点(x_i,y_i,z_i)的理论水体污染物浓度，kg/m³，

表示待求向量X在t_k时刻的第q+1(q＝n+1,...,2n)个sigma粒子时传感器节点(x_i,y_i,z_i)的理论水体污染物浓度，kg/m³；

步骤3.3、根据t_k时刻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)中：

表示t_k时刻的待求向量X的先验估计，

表示t_k-1时刻的待求向量X的后验估计，

表示t_k时刻第1个传感器节点处预测的水体污染物浓度，同式(14)，kg/m³，

表示t_k时刻第2个传感器节点处预测的水体污染物浓度，同式(14)，kg/m³，

表示t_k时刻第N个传感器节点处预测的水体污染物浓度，同式(14)，kg/m³，

表示t_k时刻第1个传感器节点处监测的水体污染物浓度，同式(1)，kg/m³，表示t_k时刻第2个传感器节点处监测的水体污染物浓度，同式(1)，kg/m³，

分别表示t_k时刻第N个传感器节点处监测的水体污染物浓度，同式(1)，kg/m³，

K_k表示t_k时刻的卡尔曼增益矩阵，

Kk=PXCSk-1---(16)]]>

式(16)中：

P_XC表示t_k时刻的估计最小均方误差协方差矩阵，

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)]]>

S_k表示t_k时刻测量噪声协方差与估计测量噪声最小均方误差协方差之和组成的测量总误差矩阵，

S_k＝P_CC+R_k(18)

式(17)、(18)中：

p表示Sigma粒子的个数，p＝1,...,n，

q表示Sigma粒子的个数，q＝n+1,...,2n，

n表示待求向量X的维数，n＝5，

W₀表示第1个sigma粒子的权值，同式(11)，

W_p表示第p+1个sigma粒子的权值，同式(12)，

W_q表示第q+1个sigma粒子的权值，同式(13)，

表示第1个sigma粒子，同式(8)，

表示第p+1个sigma粒子，同式(9)，

表示第q+1个sigma粒子，同式(10)，

表示t_k时刻的待求向量X的先验估计，

表示t_k-1时刻的待求向量X的后验估计，

Q₀表示待求向量X在t_k时刻的第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)]]>

Q_p表示待求向量X在t_k时刻的第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)]]>

Q_q表示待求向量X在t_k时刻的第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)]]>

P_CC表示t_k时刻的估计测量噪声最小均方误差协方差矩阵，

PCC=W0Q0Q0T+Σp=1nWpQpQpT+Σq=n+12nWqQqQqT---(22)]]>

R_k表示t_k时刻的测量噪声协方差矩阵，为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在t_k时刻的第1个粒子时第1个传感器节点(x₁,y₁,z₁)的理论水体污染物浓度，kg/m³，

表示待求向量X在t_k时刻的第1个粒子时第2个传感器节点(x₂,y₂,z₂)的理论水体污染物浓度，kg/m³，

表示待求向量X在t_k时刻的第1个粒子时第N个传感器节点(x_N,y_N,z_N)的理论水体污染物浓度，kg/m³，

表示待求向量X在t_k时刻的第p+1个粒子时第1个传感器节点(x₁,y₁,z₁)的理论水体污染物浓度，kg/m³，

表示待求向量X在t_k时刻的第p+1个粒子时第2个传感器节点(x₂,y₂,z₂)的理论水体污染物浓度，kg/m³，

表示待求向量X在t_k时刻的第p+1个粒子时第N个传感器节点(x_N,y_N,z_N)的理论水体污染物浓度，kg/m³，

表示待求向量X在t_k时刻的第q+1个粒子时第1个传感器节点(x₁,y₁,z₁)的理论水体污染物浓度，kg/m³，

表示待求向量X在t_k时刻的第q+1个粒子时第2个传感器节点(x₂,y₂,z₂)的理论水体污染物浓度，kg/m³，

表示待求向量X在t_k时刻的第q+1个粒子时第N个传感器节点(x_N,y_N,z_N)的理论水体污染物浓度，kg/m³，

表示t_k时刻第1个传感器节点处预测的水体污染物浓度，同式(14)，kg/m³，

表示t_k时刻第2个传感器节点处预测的水体污染物浓度，同式(14)，kg/m³，

表示t_k时刻第N个传感器节点处预测的水体污染物浓度，同式(14)，kg/m³；

式(22)中：

W₀表示第1个sigma粒子的权值，同式(11)，

W_p表示第p+1个sigma粒子的权值，同式(12)，

W_q表示第q+1个sigma粒子的权值，同式(13)；

式(23)中：

表示第1个传感器节点在t_k(k＝1,...,S)时刻的监测噪声方差，(kg/m³)²，

表示第2个传感器节点在t_k(k＝1,...,S)时刻的监测噪声方差，(kg/m³)²，

表示第N个传感器节点在t_k(k＝1,...,S)时刻的监测噪声方差，(kg/m³)²；

步骤3.4、计算t_k时刻的后验估计误差协方差矩阵

P^=P^k--KkSkKkT---(24)]]>

式(24)中：

表示t_k时刻先验估计误差协方差矩阵，

K_k表示t_k时刻的卡尔曼增益矩阵，同式(16)，

S_k表示t_k时刻测量噪声协方差与估计测量噪声最小均方误差协方差之和组成的测量总误差矩阵，同式(18)；

得到t_k时刻的待求向量X的后验估计和t_k时刻的误差协方差矩阵

步骤4、将采样时刻t_k置为第k+1个采样时刻t_k+1，即令k＝k+1，重复步骤3，直至得到k＝S时的t_S时刻的待求向量X的后验估计和误差协方差矩阵

所述t_S时刻的待求向量X的后验估计中的和依次为所求的污染源的横坐标、纵坐标和深度坐标。