基于傳感器網(wǎng)絡(luò)的分布式生化氣體源參數(shù)測定算法

張 勇 ，孟慶浩，吳玉秀，曾 明

(1. 天津大學(xué)電氣與自動(dòng)化工程學(xué)院，天津 300072；2. 天津商業(yè)大學(xué)信息工程學(xué)院，天津 300134)

隨著傳感器網(wǎng)絡(luò)技術(shù)的不斷發(fā)展，其在環(huán)境監(jiān)測領(lǐng)域的應(yīng)用也越來越廣泛，其中對(duì)生化氣體源參數(shù)測定問題[1-2]的研究在近些年受到越來越多的關(guān)注．生化氣體源參數(shù)測定是指運(yùn)用概率估計(jì)的方法確定監(jiān)測環(huán)境中的生化氣體源的位置、釋放率、擴(kuò)散系數(shù)等參數(shù)．該問題也被認(rèn)為是對(duì)生化氣體在環(huán)境中擴(kuò)散建模問題的逆問題[1]．

文獻(xiàn)[3]提出了一個(gè)完整的檢測和估計(jì)理論框架并運(yùn)用靜態(tài)傳感器陣列采集信號(hào)實(shí)現(xiàn)蒸汽源定位. 文獻(xiàn)[4]對(duì)問題進(jìn)行了擴(kuò)展，采用一個(gè)移動(dòng)傳感器節(jié)點(diǎn)取代文獻(xiàn)[3]中靜態(tài)傳感器陣列，實(shí)現(xiàn)了蒸汽源定位，同時(shí)對(duì)移動(dòng)節(jié)點(diǎn)在信號(hào)采集過程中的最優(yōu)路徑規(guī)劃問題展開研究．文獻(xiàn)[5]把應(yīng)用延伸到對(duì)地雷的監(jiān)測中，對(duì)雷場中傳感器節(jié)點(diǎn)的部署和成功檢測到地雷的概率之間的關(guān)系進(jìn)行研究．以上研究均基于相同的生化氣體物理擴(kuò)散模型且模型簡單，參數(shù)測定只針對(duì)源位置坐標(biāo)一個(gè)參數(shù)展開研究，采用的是極大似然估計(jì)算法．文獻(xiàn)等[6]提出一種基于最小二乘算法的 2步氣體源定位算法，該算法第1步先使用固定的傳感器節(jié)點(diǎn)采集濃度信號(hào)，并在給定的氣體擴(kuò)散模型的基礎(chǔ)上對(duì)氣體源的位置進(jìn)行預(yù)估運(yùn)算，給出一些類似氣體源坐標(biāo)點(diǎn)的集合；第2步對(duì)這些不同傳感器節(jié)點(diǎn)通過預(yù)估運(yùn)算得到的氣體源坐標(biāo)點(diǎn)集合求交集，其交集可以用來最后判定是否為氣味源．在該算法中運(yùn)用了氣體釋放率來對(duì)預(yù)估氣體源的位置坐標(biāo)進(jìn)行確認(rèn)．文獻(xiàn)[7]對(duì)一個(gè)放射源進(jìn)行參數(shù)測定工作，主要包括源的位置預(yù)估和釋放率的估計(jì)．

以上的生化氣體源參數(shù)測定算法均是通過傳感器網(wǎng)絡(luò)完成信號(hào)采集并將采集的信號(hào)發(fā)送到數(shù)據(jù)融合中心，生化氣體源的參數(shù)估計(jì)在數(shù)據(jù)融合中心實(shí)現(xiàn)．當(dāng)傳感器節(jié)點(diǎn)數(shù)目不太多時(shí)，這種方法實(shí)現(xiàn)簡單，隨著傳感器節(jié)點(diǎn)數(shù)目的增加，其運(yùn)算量和通信負(fù)擔(dān)將急劇增大，采用此方法將嚴(yán)重影響系統(tǒng)的性能．

筆者基于分布式卡爾曼濾波(distributed Kalman filter，DKF)算法[8]，提出了分布式擴(kuò)展卡爾曼濾波(extended Kalman filter，EKF)和無跡卡爾曼濾波(unscented Kalman filter，UKF)迭代估計(jì)算法用于氣體源的位置和釋放率的估計(jì)，并通過計(jì)算機(jī)仿真對(duì)這2種方法進(jìn)行了分析比較．

1 問題描述

生化氣體源的氣體擴(kuò)散過程可以描述為一個(gè)概率的隨機(jī)過程，因此，采用概率估計(jì)的方法對(duì)氣體源參數(shù)進(jìn)行測定應(yīng)是一種可行的解決方法．如圖 1所示，假設(shè)在二維區(qū)域中存在 1個(gè)氣體源，坐標(biāo)用(xs, ys)表示，同時(shí)在該環(huán)境區(qū)域內(nèi)隨機(jī)部署 m個(gè)傳感器節(jié)點(diǎn)，各個(gè)節(jié)點(diǎn)的坐標(biāo)已知或者可以通過節(jié)點(diǎn)自身定位算法獲得，用(xi, yi)表示第 i個(gè)節(jié)點(diǎn)的坐標(biāo)，i = 1 ,2,… ,m ．在實(shí)際應(yīng)用中，每個(gè)節(jié)點(diǎn)配置一個(gè)可采集氣體或氣味濃度信息的傳感器．假定 zk= { zk:k=1,2,…,N}為位于(xi, yi)節(jié)點(diǎn)在不同時(shí)間周期內(nèi)所采集到的濃度測量值．則氣體源測定問題可以描述為運(yùn)用不同的濃度測量值對(duì)氣體源參數(shù)向量= [,,?]T進(jìn)行估計(jì)的過程．其中，?為氣體釋放率估計(jì)值，[,]為氣體源坐標(biāo)估計(jì)值．

圖1 生化氣體源參數(shù)測定問題描述Fig.1 Geometry of chemical odor source determination problem

1.1 氣體分布模型

運(yùn)用概率估計(jì)的方法解決生化氣體源的參數(shù)測定問題需要給定氣體分布的物理模型．準(zhǔn)確地建立實(shí)際監(jiān)測環(huán)境下的氣體分布模型是非常困難的，而且對(duì)于突發(fā)氣體泄漏事件來說也是不現(xiàn)實(shí)的，因?yàn)闅怏w的擴(kuò)散通常受到湍動(dòng)氣流影響，其隨機(jī)性很大．但是，通過湍流擴(kuò)散理論可以推導(dǎo)相對(duì)簡單的用于描述時(shí)均氣體分布的模型．筆者采用文獻(xiàn)[9]中提出的時(shí)均氣體分布模型．此種模型可以描述在時(shí)均風(fēng)速恒定且均勻(homogeneous)、各向同性(isotropic)的湍動(dòng)氣流作用下的氣體分布狀況．假設(shè)生化氣體點(diǎn)源位于地平面上 ( xs,ys)處，則氣體分布模型表達(dá)式為

式中：c(xi, yi) 為監(jiān)測區(qū)域中傳感器節(jié)點(diǎn)(xi, yi)處的濃度值；q為氣體釋放率；K是湍流擴(kuò)散系數(shù)；u為風(fēng)速；d是區(qū)域中傳感器節(jié)點(diǎn)(xi, yi) 到氣體源 ( xs,ys)的距離

如果假設(shè)風(fēng)向沿著x軸，則式(1)可以簡化為

根據(jù)上述的模型可以得出在一定區(qū)域內(nèi)氣體濃度分布．

1.2 測量模型

系統(tǒng)的測量模型，即傳感器節(jié)點(diǎn)(xi, yi)的濃度測量值 zk可以描述為

式中：x =[xs, ys,q]T；非線性函數(shù) hk(x)為位置(xi, yi)處觀測到的氣體理想濃度值；vk為環(huán)境噪聲．

2 分布式氣體源參數(shù)測定算法

根據(jù)貝葉斯濾波理論[10]，所有關(guān)于狀態(tài)變量的信息都可以從其概率密度分布獲得．從統(tǒng)計(jì)方法的角度，生化氣體源參數(shù)測定可以看作在給定濃度觀測值的情況下，求解包含氣體源狀態(tài)參數(shù)信息的后驗(yàn)概率密度p ( x z1,z2,… ,zk)的狀態(tài)估計(jì)問題．

當(dāng)狀態(tài)空間模型是線性高斯模型，同時(shí)狀態(tài)噪聲和觀測噪聲都是加性高斯噪聲，通常可以采用卡爾曼濾波(Kalman filter，KF)預(yù)估狀態(tài)．而環(huán)境中生化氣體的擴(kuò)散過程是一種非線性隨機(jī)過程，因此采用針對(duì)非線性系統(tǒng)的EKF[11]和UKF[12]來完成氣體源參數(shù)向量x的估計(jì)．

EKF和UKF都把后驗(yàn)概率近似描述為一種符合N(xk; xkk, Pkk)的高斯分布，其中 xk為k時(shí)刻的狀態(tài)參量，xkk和Pkk分別為均值和方差，其值可以通過迭代的方法得到．

2.1 EKF算法

針對(duì)非線性模型，利用泰勒展開來線性化狀態(tài)方程和觀測方程，用高斯分布來近似狀態(tài)的后驗(yàn)分布，然后使用KF進(jìn)行估計(jì)．EKF算法表示為

式(8)和式(9)分別是 k-1時(shí)刻的狀態(tài)均值和狀態(tài)方差預(yù)測方程．下標(biāo)kk?1含義為在 k ?1時(shí)刻的測量值基礎(chǔ)上對(duì)k時(shí)刻狀態(tài)的估計(jì)，kk含義為得到 k時(shí)刻測量值后對(duì)狀態(tài)進(jìn)行更新的過程．

式中：Sk為估計(jì)值與真值殘差方差，Rk；Kk為Kalman濾波增益，Kk=Pkk?1；Rk為測量噪聲的方差；和分別為非線性函數(shù)fk?1(x)和 hk(x)局部線性化值．

式(10)和式(11)表述的是得到k時(shí)刻測量值 zk后系統(tǒng)的狀態(tài)均值和狀態(tài)方差的更新過程．

由于是對(duì)釋放率恒定的靜態(tài)氣體源進(jìn)行位置和釋放率估計(jì)，因此，氣體源的位置和釋放率是不變的，即向量x是恒定的．轉(zhuǎn)換方程可以簡化為

2.2 UKF算法

與EKF方法用泰勒級(jí)數(shù)近似非線性轉(zhuǎn)換函數(shù)不同，UKF用一組采樣點(diǎn)表示系統(tǒng)狀態(tài)的概率分布函數(shù)．確切地說，UKF利用這組采樣點(diǎn)通過無跡變換近似得到后驗(yàn)概率分布．無跡變換是 UKF算法的核心和基礎(chǔ)[10]，UKF利用無跡變換代替泰勒級(jí)數(shù)展開，遞推估計(jì)高斯隨機(jī)變量的均值和方差．假設(shè)系統(tǒng)采樣點(diǎn)的權(quán)重系數(shù)為，i = 1 ,2,… ,M ．權(quán)重系數(shù)可以通過無跡變換進(jìn)行更新．狀態(tài)預(yù)估均值和狀態(tài)預(yù)估方差的計(jì)算式為

測量值的預(yù)估方程為

得到 k時(shí)刻測量值kz后系統(tǒng)的狀態(tài)均值和狀態(tài)方差的更新過程為

矩陣 Pxz和 Pzz的計(jì)算式分別為

2.3 分布式迭代估計(jì)算法

分布式迭代估計(jì)算法由算法初始化、迭代運(yùn)算和算法結(jié)束判定3個(gè)部分構(gòu)成．算法流程如圖2所示．

圖2 分布式估計(jì)算法流程Fig.2 Flow chart of distributed estimation algorithm

在算法初始化階段，需要選定一個(gè)節(jié)點(diǎn)作為初始節(jié)點(diǎn)．初始節(jié)點(diǎn)可以是部署在環(huán)境中的任意節(jié)點(diǎn)，但若初始節(jié)點(diǎn)離氣體源較近，則估計(jì)過程相對(duì)容易，因此為了驗(yàn)證本文分布式迭代估計(jì)算法的可行性，初始節(jié)點(diǎn)均在距離氣味源較遠(yuǎn)的節(jié)點(diǎn)集合中選擇．初始節(jié)點(diǎn)確定以后，該節(jié)點(diǎn)在一個(gè)采樣周期內(nèi)對(duì)所處的環(huán)境進(jìn)行濃度信息采集．同時(shí)在采樣過程中該節(jié)點(diǎn)需要和鄰居節(jié)點(diǎn)完成信息共享以獲得盡可能多的測量值．初始節(jié)點(diǎn)發(fā)布廣播信息給單跳路由的鄰居節(jié)點(diǎn)，接收到廣播信息的鄰居節(jié)點(diǎn)返回其對(duì)應(yīng)的坐標(biāo)信息和相應(yīng)濃度測量值，根據(jù)這些信息由初始節(jié)點(diǎn)完成對(duì)目標(biāo)參數(shù)的預(yù)估運(yùn)算，并給出估計(jì)結(jié)果．初始節(jié)點(diǎn)得到的估計(jì)結(jié)果需要傳遞到鄰居節(jié)點(diǎn)中的某個(gè)傳感節(jié)點(diǎn)，然后進(jìn)一步完成迭代運(yùn)算，其目標(biāo)是不斷地向接近氣味源的節(jié)點(diǎn)傳遞氣體源參數(shù)預(yù)估信息．

下一個(gè)節(jié)點(diǎn)的選擇通常是由節(jié)點(diǎn)坐標(biāo)與氣體源坐標(biāo)估計(jì)值距離函數(shù)來判定，?d值最小的某個(gè)鄰居節(jié)點(diǎn)通常被選為下一個(gè)路由節(jié)點(diǎn)．該節(jié)點(diǎn)再通過自身及其鄰居節(jié)點(diǎn)提供的測量值對(duì)接收到的參數(shù)預(yù)估信息進(jìn)行更新，然后以此類推完成迭代運(yùn)算．

這種迭代估計(jì)算法所進(jìn)行的“預(yù)估-更新”過程實(shí)際上是一種卡爾曼濾波的實(shí)現(xiàn)過程，算法的最終目的是通過盡可能少的路由節(jié)點(diǎn)把初始節(jié)點(diǎn)的信息傳遞到氣味源附近的某個(gè)節(jié)點(diǎn)，并完成氣體源的定位和釋放率估計(jì)．由于不需要把測量信息傳遞到融合中心，且融合計(jì)算不集中在融合中心實(shí)現(xiàn)，而是分布到整個(gè)監(jiān)測環(huán)境中的某些節(jié)點(diǎn)上實(shí)現(xiàn)，降低了系統(tǒng)的運(yùn)算和通信消耗．

3 仿真結(jié)果與分析

為了驗(yàn)證算法的可行性，在 MATLAB平臺(tái)上進(jìn)行了仿真實(shí)驗(yàn)．仿真所用的計(jì)算機(jī) CPU主頻為2.4,GHz，內(nèi)存為 2,GB．實(shí)驗(yàn)區(qū)域?yàn)?500,m×400,m的一個(gè)二維空間，如圖3所示．

圖3 氣體源和100個(gè)傳感器節(jié)點(diǎn)分布Fig.3 Deployment of chemical odor source and 100 random sensor nodes

假設(shè)氣味源真實(shí)坐標(biāo)為(50，0)，通過文獻(xiàn)[9]中的氣體擴(kuò)散模型生成相應(yīng)的濃度數(shù)據(jù)，其分布由等高線給出，其中釋放率q選擇為 30,mL/min，湍流擴(kuò)散系數(shù)K為 30,cm2/s．100個(gè)傳感器節(jié)點(diǎn)(用黑色圓點(diǎn)表示)隨機(jī)的部署在監(jiān)測環(huán)境區(qū)域內(nèi)，假設(shè)每個(gè)節(jié)點(diǎn)的位置信息已知，在仿真實(shí)驗(yàn)中各個(gè)節(jié)點(diǎn)均可以實(shí)時(shí)采集并具有時(shí)空一致性．考慮到通常所用的金屬氧化物半導(dǎo)體氣體傳感器的相對(duì)較長的響應(yīng)和恢復(fù)，仿真中節(jié)點(diǎn)的濃度采樣周期設(shè)為 5,s．對(duì)每個(gè)傳感器節(jié)點(diǎn)而言系統(tǒng)模型噪聲采用高斯白噪聲(均值μ=0，方差 σi=3×10-6)．

對(duì)基于 EKF和 UKF的分布式迭代估計(jì)算法分別進(jìn)行了研究和實(shí)現(xiàn)，并通過多次仿真實(shí)驗(yàn)對(duì)其性能進(jìn)行了分析和比較．首先對(duì)2種算法在10種不同節(jié)點(diǎn)部署條件下的成功率以及成功實(shí)現(xiàn)的算法中參數(shù)估計(jì)誤差均值進(jìn)行了比較，然后在相同實(shí)驗(yàn)環(huán)境和條件下對(duì) 2種算法的估計(jì)結(jié)果與所需傳感器節(jié)點(diǎn)數(shù)量之間的關(guān)系進(jìn)行了分析．

在每一種節(jié)點(diǎn)部署條件下，分別選擇 10個(gè)與氣體源距離大于 400,m的不同傳感器節(jié)點(diǎn)作為初始節(jié)點(diǎn)．每個(gè)初始節(jié)點(diǎn)運(yùn)用EKF和UKF 2種不同的算法分別完成 10次迭代估計(jì)運(yùn)算，記錄 2種算法能夠成功得到估計(jì)結(jié)果的次數(shù)，然后更換節(jié)點(diǎn)部署方案重復(fù)上面的運(yùn)算．

算法完成運(yùn)算的執(zhí)行時(shí)間小于5,min則視為1次成功的迭代估計(jì)，否則視為 1次失敗的運(yùn)算(即沒有收斂)．成功率定義為成功的次數(shù)除以總的實(shí)驗(yàn)次數(shù)(每種算法的總實(shí)驗(yàn)次數(shù)為1,000次)．最后在成功的迭代運(yùn)算中，求估計(jì)誤差Δx的算術(shù)平均值．

圖4給出了不同節(jié)點(diǎn)部署條件下UKF算法和EKF算法的氣體源定位誤差均值、釋放率估計(jì)誤差均值和成功率3種統(tǒng)計(jì)結(jié)果．從圖中可以看出：UKF算法的成功率高于EKF算法，表明UKF估計(jì)算法的收斂性能好于EKF算法；UKF算法估計(jì)誤差均值比EKF小，說明UKF算法比EKF算法具有更高的估計(jì)精度．

圖4 不同節(jié)點(diǎn)部署條件下成功的迭代算法統(tǒng)計(jì)結(jié)果Fig.4 Results of successful iterative algorithm with different sensor deployment

圖 5給出了 2種估計(jì)算法時(shí)其中一次成功估計(jì)的過程．圖中，黑色圓點(diǎn)代表參與了數(shù)據(jù)融合運(yùn)算的節(jié)點(diǎn)；以黑色節(jié)點(diǎn)為圓心的黑色圓圈表示節(jié)點(diǎn)測量的濃度值大小，圓圈半徑越大，表示濃度越高．為便于性能分析和比較，EKF算法和UKF算法都選擇灰色節(jié)點(diǎn)(坐標(biāo)為(475，112))作為其初始節(jié)點(diǎn)．EKF和UKF算法最終估計(jì)的氣體源位置分別為(73，－47)和(55，11)．

圖 5(a)中的星號(hào)和圖 5(b)中的三角形符號(hào)分別表示EKF和UKF算法最終估計(jì)的氣體源位置．由圖中結(jié)果可以判定 UKF算法在定位過程中比 EKF算法使用的節(jié)點(diǎn)數(shù)目少，且路由路徑短．這也表明UKF在運(yùn)算過程中的計(jì)算量較小，從而定位速度比較快，同時(shí)通信消耗也會(huì)相應(yīng)地低．

為了更清晰地看出 EKF和 UKF算法的估計(jì)精度與所使用傳感器節(jié)點(diǎn)數(shù)量之間的關(guān)系，基于圖5所示傳感器節(jié)點(diǎn)分布及氣味源位置，圖6和圖7分別給出了氣體源位置和釋放率估計(jì)誤差隨傳感器節(jié)點(diǎn)數(shù)量的變化情況．從圖6和圖7中可以看出UKF算法定位精度和釋放率估計(jì)精度均好于 EKF算法．對(duì)于UKF算法，收斂比較快，經(jīng)過大約 20個(gè)節(jié)點(diǎn)后估計(jì)結(jié)果基本不再發(fā)生大的變化，而對(duì)于EKF算法，收斂過程相對(duì)緩慢，使用的節(jié)點(diǎn)也比較多，而且隨著節(jié)點(diǎn)的增加其最終定位誤差也比較大．

圖5 2種不同算法分布式估計(jì)算法估計(jì)實(shí)例Fig.5 Example of two different distributed estimation algorithms

圖6 氣體源位置估計(jì)誤差與傳感器節(jié)點(diǎn)數(shù)量關(guān)系Fig.6 Position estimation error of chemical odor source with different sensor node numbers

圖7 氣體源釋放率估計(jì)誤差與傳感器節(jié)點(diǎn)數(shù)量關(guān)系Fig.7 Emission rate estimation errors of chemical odor source with different sensor node numbers

4 結(jié) 語

研究了基于傳感器網(wǎng)絡(luò)技術(shù)實(shí)現(xiàn)生化氣體源參數(shù)測定問題．針對(duì)傳統(tǒng)的集中式算法中存在的問題提出了一種基于貝葉斯框架的分布式迭代估計(jì)算法，并分別運(yùn)用分布式擴(kuò)展卡爾曼濾波和無跡卡爾曼濾波方法進(jìn)行了實(shí)現(xiàn)．初步的仿真實(shí)驗(yàn)表明：分布式卡爾曼濾波算法可用于生化氣體源參數(shù)測定；無跡卡爾曼濾波算法在生化氣體源參數(shù)估計(jì)的成功率和估計(jì)誤差均值2個(gè)方面均好于擴(kuò)展卡爾曼濾波算法；在相同的仿真條件下，無跡卡爾曼濾波算法比擴(kuò)展卡爾曼濾波算法能夠更快地收斂，即使用的傳感器節(jié)點(diǎn)數(shù)目更少．

［1］ Keats W A，Yee E，Lien F S. Bayesian inference for source determination with applications to a complex urban environment[J]. Atmospheric Environment，2007，41(1)：465-479.

［2］ Oispuu M，Nickel U. Direct detection and position determination of multiple sources with intermittent emission[J]. Signal Processing，2010，90(5)：3056-3064.

［3］ Nehorai A，Porat B，Paldi E. Detection and localization of vapor-emitting sources[J]. IEEE Transactions on Signal Processing，1995，43(1)：234-253.

［4］ Porat B，Nehorai A. Localizing vapor-emitting sources by moving sensors[J]. IEEE Transactions on Signal Processing，1996，44(4)：1018-1021.

［5］ Jeremic J，Nehorai A. Landmine detection and localization using chemical sensor array processing[J]. IEEE Transactions on Signal Processing，2000，48(5)：386-395.

［6］ Matthes J，Groll L，Keller H B. Source localization by spatially distributed electronic noses for advection and diffusion[J]. IEEE Transactions on Signal Processing，2005，53(5)：1711-1719.

［7］ Pudykiewicz J A. Application of adjoint tracer transport equations for evaluating source parameters[J]. Atmospheric Environment，1998，32(9)：3039-3050.

［8］ Spanos D P，Olfati-Saber R，Murray R M. Approximate distributed Kalman filtering in sensor networks with quantifiable performance[C]//Proceedings of the 4,th international Symposium on Information Processing in Sensor Networks. Piscataway，USA，2005：133-139.

［9］ Ishida H，Nakamoto T，Moriizumi T. Remote sensing of gas/odor source location and concentration distribution using mobile system[J]. Sensors and Actuators B：Chemical，1998，49(6)：52-57

［10］ 匡興紅. 無線傳感器網(wǎng)絡(luò)中定位跟蹤技術(shù)的研究[D].上海：上海交通大學(xué)電子信息與電氣工程學(xué)院，2008.Kuang Xinghong. Research on Localization and Target Tracking in Wireless Sensor Networks[D]. Shanghai：School of Electronic Information and Electrical Engineering，Shanghai Jiao Tong University，2008(in Chinese).

［11］ Bar-Shalom Y，Li X R，Kirubarajan T. Estimation with Applications to Tracking and Navigation[M]. USA：John Wiley and Sons，2001.

［12］ Julier S，Uhlmann J，Durrant-White H F. A new method for the nonlinear transformation of means and covariances in filters and estimators[J]. IEEE Transactions on Automatic Control，2000，45(3)：477-482.