輸氣管道泄漏的波達(dá)時差交叉定位方法*

王 強(qiáng) 鄭曉亮 薛 生 袁宏永 付 明

(1 安徽理工大學(xué)電氣與信息工程學(xué)院 淮南 232001)

(2 安徽理工大學(xué)能源與安全學(xué)院 淮南 232001)

(3 清華大學(xué)合肥公共安全研究院 合肥 230601)

0 引言

隨著燃?xì)庑枨罅康闹鹉晏岣?，由輸氣管道泄漏帶來的潛在風(fēng)險不可避免地隨之增加。因此，對管道泄漏的安全監(jiān)測已經(jīng)成為公共安全領(lǐng)域的重要內(nèi)容。目前常用的大范圍監(jiān)測方法[1-4]可實現(xiàn)管道泄漏位置的粗略定位，但還需要較大范圍的搜索定位才能準(zhǔn)確找到泄漏位置，整個過程工作量大、耗時長、不利于快速排查修復(fù)漏點。泄漏發(fā)生后，管內(nèi)高壓氣體快速噴出并與泄漏孔摩擦振動產(chǎn)生泄漏聲源，聲源位置與泄漏位置高度重合，對該聲源進(jìn)行定位即可實現(xiàn)對管道泄漏的精確溯源。實際工況下必須對泄漏位置進(jìn)行三維定位即同時獲取其方位和深度信息，才具有應(yīng)用價值。

基于傳感器陣列的聲源定位方法主要包括：(1)高分辨率譜估計法[5-6]，可實現(xiàn)對波達(dá)方向的超分辨率估計，計算量較大；(2)波束形成法[7-9]，理論上可對近場聲源進(jìn)行三維定位，但計算量大幅提高，丁浩等[10]提出一種可實現(xiàn)近場聲源三維定位的波束形成方法，對2 m 以內(nèi)的3000 Hz 以上中高頻聲源定位誤差低于10%；(3)波達(dá)時差(Time difference of arrival,TDOA)法[11-15]，可實現(xiàn)對近場聲源的三維定位且計算量小，楊祥清等[16]提出一種基于球形差值的隨機(jī)梯度下降算法，實現(xiàn)了對三維空間點聲源的精確定位。對管道泄漏聲波信號特性的研究表明泄漏聲波的能量集中在低頻波段[17-18]，低頻聲波較小的傳播衰減以及較強(qiáng)的障礙穿透能力有利于信號的檢測與定位。根據(jù)近場聲源判據(jù)r＜2d2/λ(r為聲源與陣列距離，d為陣元間距，λ為波長)，由于低頻聲波波長λ較長，且實際定位中陣元間距d不能過大，因此較遠(yuǎn)距離的低頻聲源很難滿足近場判據(jù)，其更趨近于遠(yuǎn)場聲源。而平面陣列只能觀測遠(yuǎn)場聲源的空間方位信息，因此基于平面陣列的定位方法無法實現(xiàn)對遠(yuǎn)場聲源的三維定位。綜上所述，3 種方法中TDOA 法雖然能以較低的計算量進(jìn)行近場聲源三維定位，但不適用于對較遠(yuǎn)距離低頻泄漏聲源的三維定位。

在雷達(dá)定位中，常使用多站測向交叉[19-20]的方法，對不同位置基站所觀測信源方向進(jìn)行交叉以計算信源距離。為了將TDOA法應(yīng)用到管道泄漏三維定位中，本文提出一種基于TDOA的交叉定位方法，將陣列設(shè)置在泄漏區(qū)域不同位置并獲取兩組空間方位信息，對其進(jìn)行交叉求取空間交點從而完成定位。

1 管道泄漏定位原理

TDOA 法的基本原理是通過延時估計算法估計聲波信號到達(dá)不同傳感器陣元的時差，再結(jié)合陣列與聲源的幾何關(guān)系聯(lián)立方程解出聲源坐標(biāo)。由泄漏源所產(chǎn)生的聲波可分為近場球面波和遠(yuǎn)場平面波，理論上平面陣列可實現(xiàn)對近場聲源的三維定位和遠(yuǎn)場聲源的空間方位定位(相當(dāng)于二維定位)。圖1為平面陣列與兩種聲源的位置關(guān)系，當(dāng)傳感器陣列與聲源足夠近時，聲波沿不同方向傳播至4 個傳感器，設(shè)聲波到達(dá)傳感器1～3相對于參考傳感器的延時分別為τ1、τ2、τ3。結(jié)合延時τ1、τ2、τ3以及傳感器與聲源的位置關(guān)系，建立3 個不相關(guān)的方程即可解出近場聲源坐標(biāo)。

圖1 平面陣列與聲源位置關(guān)系示意圖Fig.1 Relation of position between plane array and acoustic sources

為提高定位范圍，無法始終保持傳感器陣列與泄漏位置足夠近。一方面，由遠(yuǎn)場判據(jù)r＞2d2/λ可知，將泄漏聲源視為遠(yuǎn)場聲源更有利于縮小陣列孔徑；另一方面，管道泄漏所產(chǎn)生低頻聲波的波長更長，這進(jìn)一步加劇了泄漏聲源趨向遠(yuǎn)場聲源的趨勢。因此在實際工況下，將泄漏源作為遠(yuǎn)場聲源更具有可操作性。當(dāng)泄漏源被視為遠(yuǎn)場聲源時，決定延時值的變量變?yōu)槁曉聪鄬τ陉嚵械目臻g方位而不再是三維坐標(biāo)。設(shè)傳感器1～3 與參考傳感器的距離分別為r1、r2、r3，兩傳感器連線與x軸夾角為θ1、θ2、θ3。以傳感器1為例，圖2為遠(yuǎn)場平面聲波到達(dá)傳感器1和參考傳感器的波程差示意圖，θ為聲源方位角，φ為仰角。

圖2 波程差示意圖Fig.2 Illustration of wave path difference

由圖2 中幾何關(guān)系可知，遠(yuǎn)場聲源到達(dá)傳感器1～3與參考傳感器的波程差可表示為

延時τ1、τ2、τ3由算法估計得到，結(jié)合聲速c可得3個包含變量θ、φ的獨立方程：

求解式(2)即可得出空間方位解。為進(jìn)一步實現(xiàn)對遠(yuǎn)場泄漏聲源的三維定位，本文提出一種基于TDOA法的交叉定位方法，其步驟為：(1)將傳感器陣列先后設(shè)置在泄漏區(qū)域的兩個不同位置；(2)通過TDOA 法進(jìn)行泄漏聲源空間定向；(3)結(jié)合兩組空間方位，使用交叉定位法獲取泄漏聲源的三維坐標(biāo)，完成泄漏定位。

2 延時估計

空間定向步驟的關(guān)鍵是延時估計。廣義互相關(guān)法從信號的相關(guān)分析基礎(chǔ)上發(fā)展而來，首先計算兩個傳感器所接收相關(guān)信號的互功率譜，對互功率譜進(jìn)行加權(quán)處理后經(jīng)傅里葉反變換即可得到兩個信號的互相關(guān)函數(shù)。設(shè)不同位置的兩個傳感器采集到的相關(guān)信號分別為x1(t)和x2(t)，求兩個信號的互功率譜得

其中，W(ω)表示加權(quán)函數(shù)，取1 時不進(jìn)行任何加權(quán)，此時R12(τ)表示基本互相關(guān)法的傅里葉變換形式。目前常用的加權(quán)函數(shù)有PHAT、Roth、SCOT、Echart、ML 等[21]，其中Echart 加權(quán)函數(shù)需要噪聲的先驗知識，因此本文只對基本互相關(guān)法和PHAT、Roth、SCOT、ML 加權(quán)廣義互相關(guān)法進(jìn)行時延估計效果對比。

在實驗室環(huán)境下(無強(qiáng)噪聲干擾)模擬管道泄漏，所使用管道為DN50 鋼管，內(nèi)壓0.6 MPa，泄漏孔形狀為圓形，孔徑1.5 mm。分別在距離泄漏孔0.4 m 和2 m 處設(shè)置聲波傳感器，圖3為時域連續(xù)信號的波形和頻譜圖。相較于距離泄漏孔0.4 m 處信號，2 m 處時域信號衰減幅度達(dá)到一個數(shù)量級，但500 Hz以下頻段衰減較慢?？紤]到所使用聲波傳感器頻響范圍下限為8 Hz，而高頻聲波衰減較快且穿透障礙物的能力較差，不利于提高傳感器陣列有效定位范圍，因此本文僅截取信號的10～500 Hz分量進(jìn)行分析和定位。

使用基本互相關(guān)法和多種加權(quán)廣義互相關(guān)法對兩個傳感器信號的10～500 Hz 分量進(jìn)行延時估計，保持室內(nèi)溫度為23°C，管道內(nèi)壓和泄漏孔規(guī)格不變，改變信號到達(dá)兩傳感器的距離差，每組距離進(jìn)行5 次實驗。由表1 估計結(jié)果可得，PHAT、Roth、SCOT 以及ML加權(quán)廣義互相關(guān)法的估計結(jié)果基本保持為0，基本互相關(guān)法的估計結(jié)果則隨著距離的增大而增大。使用15 組數(shù)據(jù)的基本互相關(guān)分析結(jié)果計算聲速，得到聲速估計均值為349.5 m/s，與一個標(biāo)準(zhǔn)大氣壓和15°C 條件下的340 m/s 空氣介質(zhì)聲速典型值較為接近。分析本文條件下加權(quán)廣義互相關(guān)法估計效果不佳的原因為：本文采用低通濾波截取信號10～500 Hz 分量，導(dǎo)致信號互功率譜G12(ω)高頻部分能量大幅降低。傅里葉反變換實際是對頻域信號再次進(jìn)行傅里葉變換，根據(jù)式(4)對G12(ω)進(jìn)行傅里葉反變換后，其高頻部分(可視為幅度接近0 的時域直流信號)本不應(yīng)該被反映在反變換后的互相關(guān)系數(shù)上(可視為G12(ω)的頻譜)。而多種加權(quán)函數(shù)[21]的分母均包含信號的互功率譜或自功率譜，對G12(ω)的高頻部分起到了放大作用，則反變換后G12(ω)的頻譜出現(xiàn)了直流分量，即對應(yīng)互相關(guān)系數(shù)在0 時刻出現(xiàn)了峰值。為保證延時估計精度，采樣率必須足夠高，同時為避免高頻信號干擾，有必要對信號進(jìn)行低通濾波。因此加權(quán)廣義互相關(guān)法不適合本文高采樣率低頻信號的延時估計，故采用基本互相關(guān)法進(jìn)行延時估計并完成定位。

圖3 信號波形和頻譜Fig.3 Wave form and spectrum of signal

表1 多種互相關(guān)法延時估計結(jié)果Table 1 Results of delay estimation using multiple methods

3 交叉定位

交叉定位是本文所提出定位方法的關(guān)鍵。由于TDOA空間定向必然存在誤差，空間中的兩條直線幾乎不可能相交。為解決不相交直線的交點求取問題，將兩條空間直線投影到多個平面上并在平面上求取兩條直線的交點，再將多個平面交點聚焦到空間中某一點上，從而完成偽交點的求取。圖4為基于投影法的空間直線偽交點求取過程原理圖。泄漏點g位于xOy面下方，坐標(biāo)為(x0,y0,z0)。陣列p1、p2位于xOy平面，坐標(biāo)分別為(x1,y1,0)、(x2,y2,0)。利用TDOA法可獲得泄漏點g相對于陣列p1、p2的兩組空間方位角(θ1,φ1)、(θ2,φ2)，并根據(jù)該方位角形成兩條空間直線l1、l2。為方便計算，將φ1、φ2定義為直線l1、l2與z軸負(fù)方向的夾角，θ1、θ2仍為直線l1、l2在xOy面投影與x軸正方向夾角。由于誤差的存在，圖4 中直線l1、l2均不與泄漏點g相交。將直線l1、l2投影到xOz面得交點s1，投影到xOy面得交點s2，過交點s1做平行于xOy面的平面C，最后將s2投影到平面C即為所求偽交點s。

圖4 投影法原理圖Fig.4 Schematic diagram of projection method

由p1(x1,y1,0)、p2(x2,y2,0)以及(θ1,φ1)、(θ2,φ2)得直線l1、l2：

分別將直線l1、l2投影到xOz面和xOy面得直線：

由式(7)即可得經(jīng)投影后再聚焦的偽交點坐標(biāo)s(x,y,z)。

4 管道泄漏定位實驗

4.1 實驗平臺

根據(jù)本文所提出的基于TDOA 法的管道泄漏交叉定位方法，如圖5所示搭建定位實驗平臺。該實驗平臺包含以下部分：

(1)氣體壓縮裝置：包含空氣壓縮機(jī)、儲氣罐、減壓閥以及球閥若干，當(dāng)儲氣罐充滿并關(guān)閉空壓機(jī)時，氣體壓縮裝置可向模擬測試管道提供最高1 MPa的持續(xù)穩(wěn)定壓力輸出；

(2)泄漏管道：包含一根長度為6 m 的DN50鍍鋅鋼管，預(yù)留圓形泄漏孔直徑為1.5 mm；

(3)定位陣列：包含4 個聲波傳感器和信號放大器、采集儀和配套分析軟件以及傳感器陣列支架。

圖5 定位實驗平臺Fig.5 Experimental platform for determining leakage localization

使用上述實驗平臺進(jìn)行定位實驗，調(diào)節(jié)減壓閥保持管道內(nèi)壓為0.6 MPa，關(guān)閉空壓機(jī)以避免其產(chǎn)生干擾噪聲；采集儀采樣率為3 MHz，采樣時間為10 s；4 個聲波傳感器按十字交叉排列，圖6為傳感器陣列的實物圖和示意圖。由圖3所示信號時域波形可知泄漏信號為持續(xù)穩(wěn)定信號，保持管內(nèi)壓力、管道規(guī)格和泄漏孔規(guī)格等條件不變，則兩次定位的信源條件可視為不變。傳感器陣列采用支架固定，實驗室地面平坦且支架高度固定，改變位置后可保證陣列高度、陣元相對位置等條件基本不變。按圖4所示以陣列位置連線為y軸建立空間直角坐標(biāo)系，保持陣列平面與xOz面重合(即陣列平面保持豎直)。預(yù)先設(shè)定陣列兩次布放位置的坐標(biāo)，計算陣列與泄漏點的距離和相對角度，并按照此相對位置布放陣列。因此，兩次采樣的實驗條件、陣列安裝、陣列相對位置等參數(shù)得以保證。綜上所述，實驗步驟為：儲氣罐充滿氣后關(guān)閉空壓機(jī)，調(diào)節(jié)減壓閥為管道提供0.6 MPa 壓力輸出；使用傳感器陣列先后在兩個不同位置對泄漏信號進(jìn)行采集；導(dǎo)出不同位置所采集信號，使用本文定位算法完成定位。

圖6 傳感器陣列Fig.6 Sensor array

4.2 聲源模型對比

陣列第一次布放的位置p1坐標(biāo)為(1.90 m,0,0)，第二次布放位置p2坐標(biāo)為(-1.90 m,0,0)，泄漏點g的坐標(biāo)為(0,3.10 m,-0.84 m)，則聲源與兩陣列的距離r均為3.73 m。根據(jù)本文第2 節(jié)多次延時估計結(jié)果得到的聲速平均值349.5 m/s，以最高頻率500 Hz 計算所截取信號最小波長λ為787.8 mm。如圖6(b)陣列示意圖所示可得傳感器1和傳感器4 的間距(最大陣元間距)d為54.4 cm，由判據(jù)r＞2d2/λ計算得到遠(yuǎn)場聲源與陣列之間距離的臨界值為0.75 m，因此本文實驗條件下泄漏聲源滿足遠(yuǎn)場聲源判據(jù)。根據(jù)本文第3 節(jié)對空間方位角(θ,φ)的定義，計算得泄漏點g相對于布放位置p1、p2的空間方位角坐標(biāo)為(121.5°,77.07°)和(58.5°,77.07°)。將傳感器1 設(shè)為參考傳感器，傳感器2、3、4 相對于參考傳感器延時τ2-1、τ3-1、τ4-1的估計結(jié)果如表2所示。根據(jù)陣列中各傳感器和泄漏點的位置，計算出遠(yuǎn)近場聲源的理論延時值并在表2 中一同給出。對比可知，實驗條件下的泄漏聲源更趨近于遠(yuǎn)場聲源。如圖2所示以陣列平面為xOy面建立空間直角坐標(biāo)系，根據(jù)圖6(b)所示傳感器位置關(guān)系，將延時估計結(jié)果帶入式(2)得到兩組空間方位坐標(biāo)。再將空間方位計算結(jié)果轉(zhuǎn)換至圖4所示坐標(biāo)系，得到轉(zhuǎn)換后的坐標(biāo)為(122.0°,78.1°)和(58.1°,77.0°)，與理論值較為接近。將陣列位置坐標(biāo)以及方位角帶入式(7)得偽交點s的坐標(biāo)為(-0.004 m,3.047 m,-0.791 m)，相較于與泄漏點g(0,3.10 m,-0.84 m)的誤差為0.07 m。將延時估計結(jié)果帶入近場聲源模型進(jìn)行空間定位，此時無法得出有效解，進(jìn)一步說明實驗條件下泄漏源為遠(yuǎn)場聲源。

表2 延時估計值與理論值對比Table 2 Comparison between estimated and theoretical values of time delay

4.3 定位實驗

上述定位實驗將傳感器陣列對稱放置在泄漏點兩側(cè)，但實際泄漏點位置未知，無法保證傳感器對稱放置。為模擬實際場景，仍以陣列兩次布放的位置連線作為y軸、連線中點作為原點建立空間直角坐標(biāo)系，但泄漏點分布放置。本文僅考慮陣列前向泄漏點的定位，因此實際泄漏點的分布范圍應(yīng)為xOz面前半空間?？紤]到管道多以直線鋪設(shè)，泄漏點高度變化有限，并且xOz面前半空間關(guān)于yOz面對稱、陣列位置關(guān)于原點對稱，討論yOz面右半空間的定位情況即可推廣至左半空間。因此，定位實驗劃定泄漏點區(qū)域為x、y坐標(biāo)取[0,10 m]范圍、z坐標(biāo)取-0.84 m。圖7為陣列位置與泄漏位置的俯視圖，陣列兩次布放的位置分別為(M1,N1)、(M2,N2)和(M3,N3)，A1～C3為9個泄漏位置(每次實驗僅一個位置發(fā)生泄漏)，各點x和y坐標(biāo)如圖7所示，陣列位置z坐標(biāo)為0，泄漏位置z坐標(biāo)為-0.84 m。為了探究陣列孔徑對定位精度的影響，還需要改變陣列模型孔徑參數(shù)進(jìn)行實驗。記圖6(b)所示陣列模型孔徑為孔徑1，進(jìn)一步縮小孔徑使傳感器1、3 與十字交點的距離為35 cm，傳感器2、4與交點距離為15 cm，記為孔徑2；傳感器1、3 與交點的距離為20 cm，傳感器2、4 與交點距離為15 cm，記為孔徑3。

圖7 陣列與泄漏位置俯視圖Fig.7 Position of arrays and leakage points

使用3 種孔徑的陣列，在3 組陣列布放位置下分別對9 個泄漏位置進(jìn)行定位實驗，每個條件下進(jìn)行5 次實驗并計算平均定位誤差。表3為定位結(jié)果，對于最大的孔徑1，陣列布放位置為M1、N1(即陣列間距2 m)時，A 組泄漏點定位精度相對較高，B、C兩組明顯下降；陣列布放位置為M2、N2(即陣列間距6 m)時，整體定位誤差有所降低；陣列布放位置為M3、N3(即陣列間距10 m)時，A、B 兩組泄漏點定位誤差明顯降低，C 組有所下降但仍較高。同樣條件下，A 組定位精度高于B、C 兩組，A 組內(nèi)距離原點最遠(yuǎn)的A1點定位誤差最大。隨著陣元間距的減小，定位誤差呈增大趨勢，且孔徑3 下誤差明顯增大。分析原因為

(1)A 組泄漏點位于y軸上而B、C 兩組偏離y軸，則A 組各點相對于兩個陣列位置的方位角θ分別屬于(0,90°)和(90°,180°)范圍且關(guān)于θ=90°對稱，B、C 兩組各點兩個方位角θ均小于90°且距離原點越遠(yuǎn)兩方位角大小越接近，由式(7)中x、y坐標(biāo)表達(dá)式可知計算A 組各點坐標(biāo)時分母tanθ1-tanθ2的絕對值更大，此時分子值的上下浮動所引起的誤差更小，反之計算B、C兩組各點坐標(biāo)時分母絕對值較小，則分子波動對誤差影響增大；

(2)對于A 組各點，距離原點越遠(yuǎn)或陣列間距越小則方位角越接近90°，由正切函數(shù)變化率可知此時方位角θ的誤差對定位誤差的影響更大，因此表3 中同一條件下A1點的定位誤差相對于A2、A3更大，陣列布放位置為M1、N1(陣列間距最小)時A組定位誤差同樣更大；

(3)對于B、C兩組，陣列布放位置為M1、N1(陣列間距最小)時，各點兩個方位角θ的差值較小，由原因(1)可知計算坐標(biāo)時分子波動對誤差影響較大，擴(kuò)大陣列間距則誤差得以降低；

(4)陣列孔徑?jīng)Q定空間方位計算的分別率，孔徑較小時信源方向的改變所引起的延時變化較小，不利于互相關(guān)法對延時進(jìn)行準(zhǔn)確估計，從而導(dǎo)致定位誤差增大。

對比本文方法與現(xiàn)有基于TDOA 的聲源空間定位方法。楊祥清等[16]使用SI-LMS算法在辦公室環(huán)境進(jìn)行語音信號定位，定位結(jié)果表明：SI-LMS 算法對距離陣列1 m 以內(nèi)的近場聲源定位誤差小于0.05 m，對2～3 m 的中場聲源定位誤差小于0.1 m，對4 m 以上的遠(yuǎn)場聲源定位誤差則顯著提高，最大誤差可達(dá)到0.57 m。本文方法定位結(jié)果顯示：陣列布放于M3、N3時，距離陣列或原點較近的泄漏點(A3、B3點)定位誤差與SI-LMS 法基本一致，但距離超過4 m 且靠近y軸的A1、A2、B1和B2點定位誤差仍保持較低，而偏離y軸的C1、C2和C3點定位誤差則明顯增加。相較于現(xiàn)有適用于近場聲源三維定位的TDOA 法，本文方法對y軸附近較遠(yuǎn)距離的泄漏點定位精度有較明顯提升。

表3 定位結(jié)果Table 3 Results of localization

綜上所述，使用本文方法進(jìn)行泄漏定位時，在有效信號檢測范圍內(nèi)，應(yīng)增大陣列布放間距；在不影響遠(yuǎn)場聲源判據(jù)成立的前提下，應(yīng)使用較大孔徑陣列；當(dāng)定位結(jié)果偏離y軸時應(yīng)適當(dāng)調(diào)整陣列位置，以保證泄漏點位于y軸附近，提高定位精度。

5 結(jié)論

為實現(xiàn)對管道泄漏遠(yuǎn)場聲源的三維定位，本文使用交叉定位法對現(xiàn)有TDOA 定位方法進(jìn)行了改進(jìn)。將泄漏聲源視為遠(yuǎn)場聲源并使用TDOA 法進(jìn)行不同位置的兩次定向，提出空間不相交直線的偽交點求取方法，對兩組空間方位進(jìn)行交叉求取交點，從而完成定位。針對本文信號采樣率高、頻率低的特點，選取基本互相關(guān)法進(jìn)行延時估計。建立管道泄漏定位實驗平臺，對陣列孔徑、布放間距以及泄漏位置等因素對定位精度的影響進(jìn)行了實驗分析。定位結(jié)果表明：適當(dāng)提高陣列間距、增大陣列孔徑以及調(diào)整陣列指向，能夠提高定位精度；與現(xiàn)有TDOA法相比，基于交叉定位改進(jìn)的TDOA法對距離原點4 m 以上遠(yuǎn)場泄漏聲源的定位精度有明顯提升。

由于實驗條件所限，本文采用陣列先后兩次布放的方式進(jìn)行定位，且陣元數(shù)和陣列孔徑無法進(jìn)一步擴(kuò)大。實際定位中，應(yīng)采用雙陣列同時布放、采集的方式提高定陣列參數(shù)的一致性和可控性，同時應(yīng)進(jìn)一步探索更多陣元數(shù)、更復(fù)雜陣列模型以及更大陣列孔徑對定位精度的有益影響。