(北京理工大學(xué)自動化學(xué)院， 北京 100081)

引言

現(xiàn)代工業(yè)和生活中廣泛使用了用于存儲和輸送壓縮氣體的壓力容器以及管道等， 氣密性是這些設(shè)備質(zhì)量和安全的重要指標(biāo)之一，因此氣體泄漏檢測非常重要[1-2]。一些傳統(tǒng)的泄漏檢測方法，比如壓差法、絕對壓力法等，操作非常復(fù)雜并且對技術(shù)人員要求較高，而且一般不具有實時性[3]。現(xiàn)代工業(yè)中廣泛利用泄漏產(chǎn)生超聲波的原理來進(jìn)行泄漏檢測， 通過分析高壓氣體從密封容器或者管道的泄漏孔噴射產(chǎn)生的超聲波，達(dá)到泄漏檢測與泄漏點(diǎn)定位的目的[4]。國內(nèi)外學(xué)者根據(jù)不同研究重點(diǎn)開發(fā)了多種超聲波氣體泄漏檢測儀， 采用的檢測定位算法主要包括以時延技術(shù)為基礎(chǔ)的到達(dá)時間差定位算法(TDOA)和以信號傳播特性為基礎(chǔ)的能量衰減算法(ED)。然而到達(dá)時間差定位算法誤差不穩(wěn)定：當(dāng)傳感器與泄漏點(diǎn)距離較大時，能量衰減定位算法所得的聲壓比值接近1，這使得該方法對誤差非常敏感，故采用基于擴(kuò)展Kalman濾波的數(shù)據(jù)融合方法，利用傳感器陣列，將TDOA定位法和ED定位法的計算結(jié)果進(jìn)行聯(lián)合計算出泄漏點(diǎn)與陣列中各換能器之間的距離，最終通過幾何方法計算出泄漏點(diǎn)與傳感器平面的垂直距離和相應(yīng)的投影坐標(biāo)。通過實驗證明了該算法的有效性及其與傳統(tǒng)算法相比的優(yōu)越性。

1 氣體泄漏模型

對于各種液壓、氣壓密閉容器，管道或是焊接加工接口等由于加工質(zhì)量不過關(guān)、安裝不合理，或者是長期連續(xù)使用等原因都會導(dǎo)致細(xì)微的孔隙和裂縫產(chǎn)生，因此在壓力系統(tǒng)的作用下就會產(chǎn)生流體泄漏[5]。當(dāng)孔隙的尺寸足夠小，而容器的內(nèi)外壓力差足夠大的時候，從孔隙中泄漏產(chǎn)生的氣體的流速會很大，而且泄漏的氣體的雷諾數(shù)通常比較高，因此就產(chǎn)生了湍流射流。對于湍流射流來說，從靠近狹縫處的初始段至主段以后很長的區(qū)域內(nèi)將會不斷地產(chǎn)生旋渦，而這些漩渦又在不斷地發(fā)展破裂產(chǎn)生新的漩渦，氣體泄漏模型如圖1所示。關(guān)于這些漩渦，1952年Lighthill就給出了論述[6]，即漩渦產(chǎn)生的渦流實際上就是流體的聲音。

圖1 氣體泄漏模型

實際工程中，容器或管道的泄漏孔很小，其孔徑大多在0.1 mm量級，甚至更小。超聲波的產(chǎn)生區(qū)域在噴口以外5倍孔徑的區(qū)域，即該區(qū)域長度在1 mm量級。通常的泄漏檢測中，傳感器與泄漏孔的距離在100 mm量級上，遠(yuǎn)大于超聲波的產(chǎn)生區(qū)域，故可以將噴注噪聲產(chǎn)生區(qū)域看做是點(diǎn)聲源，聲波以球面向外傳播。泄漏產(chǎn)生的聲波的頻率譜范圍一般是從十幾千赫茲到接近百千赫茲，但是能量大都會集中在10～50 kHz內(nèi)，屬于超聲波頻譜范圍之內(nèi)，因此利用超聲波傳感器陣列對其進(jìn)行采集、處理，可以估計出點(diǎn)聲源的位置，即泄漏點(diǎn)的位置。

2 到達(dá)時間差定位算法

到達(dá)時間定位方法的原理是測量信號從信號源到達(dá)觀測點(diǎn)的時間，進(jìn)而得到信號源與觀測點(diǎn)之間的距離，多觀測點(diǎn)的信息聯(lián)合解算出信號源位置。TDOA定位法可以消除對時間基準(zhǔn)的依賴性，因而可以降低成本并仍然保證一定的定位精度，但是需要有較好的時延估計方法，才能保證較高的時延估計精度[7-8]。

圖2為泄漏定位算法原理圖，假設(shè)4個超聲波換能器A，B，C，D以正方形平面陣列的形式放置，其坐標(biāo)為A(x1,y1,z1)，B(x2,y2,z2)，C(x3,y3,z3)，D(x4,y4,z4)，信號到達(dá)換能器的時間為t1，t2，t3，t4，信號到達(dá)換能器的時間差為Δt12，Δt23，Δt34(下標(biāo)1對應(yīng)換能器A，2對應(yīng)換能器B，以此類推)。

圖2 泄漏定位算法原理圖

則旋轉(zhuǎn)雙曲面方程組如下：

(1)

式中，ξij為到達(dá)時間差計算值的誤差。方程組為非線性恰定方程組，但是到達(dá)時間差誤差為隨機(jī)數(shù)，因此上述方程組很難存在解。上述雙曲面方程組聯(lián)合其他數(shù)學(xué)估計方法，例如最小二乘法或牛頓迭代法等，可以解出方程組的近似解，該解即為泄漏點(diǎn)位置[9]。圖3表示泄漏點(diǎn)位于(-10,10,450)、初始迭代值為(10,-10,350)、到達(dá)時間差誤差在5 μs內(nèi)時，使用牛頓迭代法計算泄漏點(diǎn)位置的迭代過程?？梢钥闯?，迭代過程中系統(tǒng)狀態(tài)收斂到局部極值，與泄漏點(diǎn)真實位置相差巨大。

圖3 到達(dá)時間差算法仿真結(jié)果

3 能量衰減定位算法

在使用能量衰減定位算法時，主要利用的是信號的能量在傳播過程中隨傳播距離的增大而逐漸衰減的原理，導(dǎo)致不同換能器接收到的來自同一泄漏聲源的超聲波信號在強(qiáng)度上的不同，從而構(gòu)成約束條件，與泄漏點(diǎn)和超聲波換能器陣列構(gòu)成的幾何位置關(guān)系，聯(lián)立方程組得出泄漏孔所在的位置坐標(biāo)。能量衰減定位方法因為硬件簡單、低成本等特點(diǎn)，在無線傳感器網(wǎng)絡(luò)領(lǐng)域的節(jié)點(diǎn)定位中使用較多[10-11]。

超聲波可被視為球面波，將聲波的衰減模型描述為：

(2)

式中，pi(t)為t時刻換能器i處聲壓值；p(t)為t時刻聲源的聲壓值；rs與ri分別為聲源與換能器的位置向量；npi(t)為聲壓噪聲，可以用高斯分布近似。

本研究選擇的壓電式超聲波換能器對聲壓敏感，換能器的輸出值正比于聲波的聲壓，即：

ui(t)=gipi(t)

(3)

式中，gi為換能器i的接收系數(shù)。

考慮不存在噪聲的情況下，換能器i,j處聲壓比值kij為：

(4)

則實際檢測時，由平面陣中傳感器測得信號得到的定位方程組為:

(5)

簡單地說，到達(dá)時間差定位算法是根據(jù)泄漏點(diǎn)與任兩傳感器距離的差值做出以傳感器位置為焦點(diǎn)、距離差為長軸的空間旋轉(zhuǎn)雙曲面，多個雙曲面的交點(diǎn)即是泄漏點(diǎn)位置。能量衰減算法是根據(jù)泄漏點(diǎn)與任兩傳感器距離的比值得出一球面，多個球面的交點(diǎn)為泄漏點(diǎn)位置，如圖4所示。

圖4 面型傳感器陣列定位示意圖

4 基于Kalman濾波的融合定位算法

數(shù)據(jù)融合就是將多傳感器測得信息或者多種算法得到的最終或中間結(jié)果聯(lián)合來提高最終結(jié)果或估計的精確度的方法[12-13]。

數(shù)據(jù)融合可分為3個層次，數(shù)據(jù)層融合、特征層融合與決策層融合。數(shù)據(jù)層融合是將多傳感器測得的最底層數(shù)據(jù)送至數(shù)據(jù)融合中心進(jìn)行處理，得到需要的結(jié)果，計算量大時效性差；特征層融合是先抽象出多傳感器或者多算法的特征向量，將各特征向量送至融合中心進(jìn)行關(guān)聯(lián)，根據(jù)融合結(jié)果做出判斷決策；決策層融合是先利用各傳感器或算法分別計算出結(jié)果，再將結(jié)果融合得到最終值，精度較低。

本研究是兩種算法的融合，無法進(jìn)行數(shù)據(jù)層融合。決策層融合相當(dāng)于根據(jù)兩種算法分別求出泄漏點(diǎn)位置，再對各自的估計結(jié)果進(jìn)行融合。由于單獨(dú)的到達(dá)時間差算法與能量衰減算法估計結(jié)果誤差大，不適合決策層融合。本研究中，到達(dá)時間差算法(TDOA)得到泄漏點(diǎn)與傳感器的距離差,能量衰減算法(ED)得到泄漏點(diǎn)與傳感器的距離比值。特征層融合是根據(jù)距離差與距離比值送到融合處理中心進(jìn)行融合計算，本研究采用特征層融合。

數(shù)據(jù)融合常用方法包括概率論方法、D-S證據(jù)理論、神經(jīng)網(wǎng)絡(luò)與Kalman濾波等，本研究采用Kalman濾波方法，該方法適用于實時處理動態(tài)環(huán)境中的信息。Kalman濾波屬于遞歸濾波器，通常被用于動態(tài)系統(tǒng)中估計含有噪聲的系統(tǒng)狀態(tài)值。然而Kalman濾波理論只能應(yīng)用于線性系統(tǒng)，本研究的測量方程涉及泄漏點(diǎn)到傳感器的距離比值，含有非線性成分，所以本研究采用擴(kuò)展Kalman濾波方法(EKF)進(jìn)行數(shù)據(jù)融合。

本研究采用的數(shù)據(jù)融合的結(jié)構(gòu)框圖如圖5所示。將到達(dá)時間差算法與能量衰減算法視為兩子系統(tǒng)，分別得到到達(dá)時間差與聲壓比值。圖5中特征抽取表示將到達(dá)時間差與聲壓比值轉(zhuǎn)換為到達(dá)距離差與到達(dá)距離比，關(guān)聯(lián)即將到達(dá)距離差與到達(dá)距離比聯(lián)合作為數(shù)據(jù)融合的觀測值。數(shù)據(jù)融合部分采用擴(kuò)展Kalman濾波理論進(jìn)行泄漏點(diǎn)位置的迭代計算，最終得到位置估計結(jié)果。

圖5 數(shù)據(jù)融合結(jié)構(gòu)框圖

關(guān)于數(shù)據(jù)融合部分的模型，本研究將超聲波泄漏檢測類比為一個控制系統(tǒng)，假設(shè)泄漏點(diǎn)位置坐標(biāo)x，y，z為系統(tǒng)狀態(tài)，r=[x,y,z]T為狀態(tài)向量，根據(jù)不同傳感器測得的時間差和聲壓比計算距離差Δsij與距離比值δsij，作為系統(tǒng)的特征值，將數(shù)據(jù)融合看作對系統(tǒng)狀態(tài)的估計。

首先，建立系統(tǒng)的狀態(tài)轉(zhuǎn)移方程與觀測方程。本研究考慮泄漏檢測時泄漏點(diǎn)相對于傳感器陣列靜止的情況，則該離散控制過程的狀態(tài)方程可以描述為:

r(k+1)=Φr(k)+ω(k)

(6)

式中，本研究中泄漏點(diǎn)坐標(biāo)的單位為mm。Φ為狀態(tài)轉(zhuǎn)移矩陣，Φ=I3為單位陣。ω(k)∈R3為過程噪聲，可以將其元素當(dāng)做均值為0的高斯白噪聲?？梢詫⒌竭_(dá)時間差算法與能量衰減算法看做系統(tǒng)的兩子系統(tǒng)，前者的特征向量為泄漏點(diǎn)到傳感器的距離差值Δs，后者的特征向量為距離比值δs。本研究采用特征層的數(shù)據(jù)融合，即將兩子系統(tǒng)的特征向量關(guān)聯(lián)，則將系統(tǒng)的觀測方程表示為：

y(k)=h(r(k))+υ(k)

(7)

式中，y(k)=[Δs12(k), Δs23(k), Δs34(k),δs12(k),δs23(k),δs34(k)]T為系統(tǒng)的觀測向量； Δsij=c·Δtij，δsij=kij，v(k)∈R6表示觀測噪聲，其元素可用均值為0的高斯白噪聲近似。其中，h(r(k))為系統(tǒng)的非線性觀測函數(shù):

(8)

其中，si(k) (i=1,2,3,4)表示泄漏點(diǎn)與傳感器i之間的距離：

(9)

過程噪聲主要由檢測過程中傳感器陣列與泄漏點(diǎn)的相對移動造成，觀測噪聲由信號調(diào)理采集環(huán)節(jié)以及計算過程引入，故認(rèn)為過程噪聲與觀測噪聲沒有相關(guān)性。

根據(jù)擴(kuò)展Kalman濾波理論，系統(tǒng)的先驗狀態(tài)估計為:

(10)

P(k+1|k)=ΦP(k|k)ΦT+Q

(11)

由k+1時刻的觀測值，得到k+1時刻實際觀測值與先驗估計的觀測值間的差，也稱為新息：

(12)

系統(tǒng)的后驗估計為:

(13)

其中，K(k+1)代表kalman系數(shù)，也稱為kalman增益，為:

K(k+1)=P(k+1|k)HT(k+1)

[H(k+1)P(k+1|k)HT(k+1)+R]-1(14)

H(k+1)表示根據(jù)k+1時刻先驗狀態(tài)估計局部線性化得到的觀測矩陣，為h(r(k))的雅克比矩陣。后驗狀態(tài)估計的誤差協(xié)方差矩陣為:

P(k+1|k+1)=[I3-K(k+1)H(k+1)]

P(k+1|k)

(15)

給定迭代的初始狀態(tài)估計及其誤差協(xié)方差為:

P(0|0)=P0

(16)

根據(jù)以上的濾波模型進(jìn)行擴(kuò)展Kalman濾波的迭代，直至相鄰兩次系統(tǒng)狀態(tài)向量差的模的平方小于某一固定值時，輸出融合結(jié)果。

5 實驗證明

為驗證本研究提出的泄漏定位算法的性能，建立了一套由4個傳感器組成的面陣、檢測機(jī)箱、顯示器和存在0.1 mm直徑小孔的尼龍罐構(gòu)成的實驗裝置，實驗裝置如圖6所示，檢測機(jī)箱內(nèi)包含前置調(diào)理電路、數(shù)據(jù)采集卡、單板電腦、電源和風(fēng)扇。傳感器采集的信號經(jīng)過前置調(diào)理電路的放大和濾波、數(shù)據(jù)采集卡的數(shù)模轉(zhuǎn)換后，輸入單板電腦進(jìn)行定位計算，最終的定位結(jié)果顯示在顯示屏中，使用的軟件平臺為LabVIEW。

圖6 實驗裝置

4個傳感器按照四邊形邊長為60 mm的正方形排列。實驗時，尼龍罐內(nèi)空氣相對壓力為100 kPa。設(shè)4個傳感器的坐標(biāo)分別為A(-30,30,0)，B(30,30,0)，C(30,-30,0)，D(-30,-30,0)(單位為mm)。將泄漏點(diǎn)固定在距離傳感器面不同距離的不同位置，然后分別使用TDOA、ED與本研究提出的基于數(shù)據(jù)融合的定位算法估計泄漏點(diǎn)位置。實驗將泄漏點(diǎn)固定于S1(-10,70,300)，S2(-80,50,500)，S3(-20, -20,800), 3個位置。本研究將泄漏點(diǎn)位置計算結(jié)果的誤差定義為結(jié)果與實際泄漏點(diǎn)位置的三坐標(biāo)差值的平方和根，也即是：

(17)

式中，e為位置結(jié)果誤差;xe,ye,ze為定位算法計算得到的泄漏點(diǎn)位置，不同算法的計算結(jié)果誤差如表1所示。

表1 不同算法的實驗誤差

如表1所示，基于Kalman濾波的數(shù)據(jù)融合定位算法的精度明顯高于ED定位算法、TDOA定位算法，這主要是因為ED算法在容器壓力較低的情況下，聲壓比接近于1，因此該算法對誤差極為敏感，定位精度一直不高，為計算帶來了困難；而TDOA算法的實驗誤差在一般情況下是較小的，但有時會突然變得非常大，實驗誤差很不穩(wěn)定。這是由于該算法主要是通過計算2個信號的相位差從而得到時間差，當(dāng)測試距離過大時，不能準(zhǔn)確得到相位差為θ或是(θ+360°·n)(n為任意自然數(shù))，當(dāng)估計結(jié)果與真實結(jié)果相差整整一個周期時，實驗誤差達(dá)到了340.28 mm?；贙alman濾波的數(shù)據(jù)融合定位算法克服了兩種算法的局限性，結(jié)果誤差小于20 mm。

6 結(jié)論

本研究針對基于超聲波的氣體泄漏檢測與定位，提出了一種基于Kalman濾波的的泄漏點(diǎn)定位算法，設(shè)計了一套基于該算法的泄漏定位裝置并利用該裝置驗證了所提出的算法。本研究的創(chuàng)新點(diǎn)在于將平面四元陣應(yīng)用于超聲波泄漏檢測定位中；設(shè)計到達(dá)時間差算法和能量衰減算法在平面陣上的定位計算模型；關(guān)聯(lián)兩種算法的特征向量進(jìn)行泄漏點(diǎn)的位置估計，即將數(shù)據(jù)融合應(yīng)用于泄漏點(diǎn)定位。實驗結(jié)果表明，本研究提出的定位算法能計算出泄漏點(diǎn)位置，結(jié)果誤差小于20 mm。

同時，該方法也存在一定的局限性。首先，在泄漏信號較小的情況下，定位系統(tǒng)受到環(huán)境噪聲的影響較大，難以檢測到氣體泄漏，因此，需要改進(jìn)放大濾波電路。此外，我們正在計劃增加傳感器陣列上的傳感器數(shù)量，并研究傳感器陣列的布局模式對定位結(jié)果的影響，以找到進(jìn)一步提高定位系統(tǒng)精度的方法。