[發明專利]一種接近湖泊岸邊的三維污染源定位方法

在審

• 申請號: 201710254491.6
• 申請日: 2017-04-18
• 公開（公告）號: CN107423467A
• 公開（公告）日: 2017-12-01
• 發明（設計）人: 柴利;馮政;楊君;羅旭;吳晨暉
• 申請（專利權）人: 武漢科技大學
• 主分類號: G06F17/50
• 分類號: G06F17/50;G01N33/18;G06Q50/06
• 代理公司: 武漢科皓知識產權代理事務所(特殊普通合伙)42222
• 代理人: 張火春

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的后驗估計中的和依次為所求的污染源的橫坐標、縱坐標和深度坐標。