基于聲波衰減模型對液體管道泄漏位置的極大似然估計

池招招，蔣軍成,，刁 旭，倪 磊，王志榮，沈國棟，郝永梅

(1.南京工業大學 安全科學與工程學院，南京 211816；2.常州大學 環境與安全工程學院，江蘇 常州 213164)

隨著我國工業化進程的加快和能源結構優化的推進，我國油氣管道建設正迎來一個大的發展機遇期。雖然油氣管道運輸為經濟發展帶來諸多便利，但是油氣管道安全影響因素多，部分因素偶發性強，運行風險高，管理難度大[1-3]。油氣管道泄漏可能導致爆炸等嚴重事故，造成巨大的生命財產損失和環境問題。因此泄漏檢測方法的研究將有助于提升我國油氣集輸系統安全保障的技術水平，具有重要的經濟價值和社會意義。

針對管道泄漏定位研究，Li等[4]通過對信號使用時頻分析譜(crosstime-frequency spectrum，CTFS)去噪，然后通過互相關分析確定上下游傳感器接收到泄漏信號的時差，使用時差法對泄漏點進行定位，所得結果相較單純時差法誤差降低六倍。Liu等[5]研究了聲波在管道上的傳播衰減公式，并對該公式進行修正，通過實驗獲得管道上下游處聲波幅值，使用聲波衰減公式計算得到泄漏點位置。該方法與時差法相比，無需獲得聲波傳播波速和上下游傳感器接收到泄漏信號的時差，且Liu等將該方法成功應用到實際油氣管道上。Sun等[6]仍舊使用時差法進行泄漏定位，但不同點是使用了威格納-維爾分布(wigner-ville distribution，WVD)計算時頻參數，并通過時頻參數獲得時差。Sun等是將時差法用于聲發射(acoustic emission, AE)信號，Liu等[7]則是在壓力信號中使用時差法。Liu等[8]使用了改進之后的小波包對動態壓力傳感器獲取的信號進行去噪，并對比時差法與幅值衰減法的定位情況，結果表明幅值衰減法定位結果優于時差法，且不受時間測量與波速測量裝置的限制。

由于管道中液體流動的影響，上下游衰減參數是有差異的；使用聲波衰減公式定位時，若工況發生變化，則需要重新計算衰減參數。本文針對上述問題改進了衰減參數的設置，分別計算上下游參數值；且為了避免每次計算衰減參數前需要確定管道運行參數，改進了傳感器的布設方式，直接通過實驗數據獲得衰減參數值，然后使用極大似然估計(maximum likelihood estimation, MLE)法對改進后的衰減公式進行變換并求解，從而獲得預測的泄漏位置。

1 原 理

1.1 聲波衰減模型

泄漏狀態下，流體從泄漏口噴射至管外時，液體與管壁相互作用產生聲發射源，聲發射源向外輻射能量形成壓力波并以波速沿著管壁進行傳播。通過在管壁上安裝AE傳感器則可記錄得到聲波信號。

聲波信號隨傳播距離的增加而發生的衰減公式[9]為

p=p0e-σαx

(1)

式中：x為聲波傳播距離，m；p0為泄漏點處的聲波幅值，Pa；p為相應傳播距離x處的聲波幅值，Pa；σ為修正系數，用于修正管道運輸不同介質時所產生的差異，對于水管，取值3～4，而對于氣體管道則是0.5～1.5；α為黏熱吸收系數，用于表征由于黏性吸收和熱傳導所造成的聲波強度衰減，根據公式(1)可知，α值越大，黏性吸收和熱傳導對聲波的衰減作用越強，α計算公式為

α=αη′+αη″+αχ=

(2)

式中[10]：αη′為切變黏滯吸收系數;αη″為容變黏滯吸收系數；αχ為熱傳導吸收系數；ω為角頻率，rad/s，ω=2πf；f為聲波頻率，Hz；r為管道內徑，m；ρ0是介質密度，kg/m3；c0為聲波在管道中的傳播速度，m/s；η為動力黏滯系數，Pa·s；η′為切變黏滯系數，Pa·s；η″為容變黏滯系數，Pa·s；χ為熱傳導系數,W/(m·K)；Cv為定容比熱；Cp為定壓比熱,kJ/(kg·K)。

若是按照傳統的聲波衰減定位方法，則需要根據管道運行狀態確定上述所有參數，從而獲得黏熱吸收系數α。當所需檢測的管道發生變化時，其對應的r，c0，η′，η″，χ相應地也需要做出調整，而當管道內介質發生變化時，需要相應調整ρ0，Cv和Cp，這就造成了計算上的繁瑣以及變換應用場合之后的重復計算。

因此本文在該基礎上，增加了上下游衰減系數的差異性；且改變了傳統的傳感器布置方式，以便于σ和α的計算，傳感器布置示意圖如圖1所示。

圖1 水管泄漏聲波檢測的傳感器布置示意圖

按照式(1)可知，對于水管修正系數σ取值為3～4，只給出取值范圍而沒有確切的數值，所以本文為了避免修正系數取值的差異對結果的影響，將修正系數σ和黏熱吸收系數α合并，在本文中兩系數乘積σα即為聲波衰減系數。

根據式(1)，聲波向上游點1和點2的傳播公式為式(3)和(4)，上游的衰減系數為σ1α1

p1=p0e-σ1α1x1

(3)

p2=p0e-σ1α1x2

(4)

同理，聲波向下游點3和點4的傳播公式如式(5)和式(6)，下游的衰減系數為σ2α2

p3=p0e-σ2α2x3

(5)

p4=p0e-σ2α2x4

(6)

式中，p0為泄漏點處的聲波信號幅值，mV。(式(1)的聲波為壓力值，而實驗設備所得幅值為電信號的電壓值，本文直接使用電信號的電壓值代替聲波壓力值)；x1，x2，x3和x4分別為傳感器1，2，3，4與泄漏點之間的距離，m；p1，p2，p3和p4分別為傳感器1，2，3，4所獲取的聲波信號幅值，mV；σ1α1可由上游兩個傳感器的幅值比值獲得，即式(3)除以式(4)

(7)

變換公式獲得

(8)

同理可得

(9)

(10)

從上式可以看出，只要知道聲發射傳感器在測量點1、2、3、4處的聲波信號幅值p以及同一端傳感器之間(1和2之間，3和4之間)的距離，即可由式(8)和(10)獲得衰減參數的數值。

1.2 聲波衰減模型的極大似然估計

公式(8)和(10)僅能得到σα的值，如果使用衰減公式(1)計算x，則需要知道泄漏點處的聲波幅值，但是實際工程應用中，泄漏點處的聲波幅值是無法得到的，僅僅用該公式并不能得到泄漏位置。所以對于泄漏位置的預測，本文在聲波衰減公式的基礎上采用MLE。MLE的估計流程圖如圖2。

圖2 MLE流程圖

按照上述流程，首先定義一個上下游的幅值比值公式為

(11)

由于泄漏檢測過程中，信號常常夾雜著噪聲，所以本文在公式(11)的基礎上加入白噪聲ω[11]

(12)

假設白噪聲ωi滿足N(μi,εi)的正態分布，對式(12)標準化處理，得

(13)

(14)

對式(14)聯合概率函數取對數可得似然函數，如式(15)所示

(15)

式中：μi為噪聲均值；εi為噪聲方差；L為測量段長度，m，即傳感器1和4之間的距離。

上式中，Ki，σ1α1和σ2α2可通過去噪后的純泄漏信號進行幅值計算，μi和εi由噪聲信號獲得，將實驗信號所得到的數值代入式(15)，使用MLE即可對該公式進行求解，求得x1。

2 實驗設備與操作

由于測量油品泄漏會有較大危險性，所以本文采用水代替油品進行實驗。

實驗裝置包括管道系統和數據采集系統，如圖3。該水管由法蘭連接，法蘭連接處有橡膠墊圈，管長67 m，閥門安裝在距離上游水流入口33.5 m的地方。管道外徑108 mm，壁厚3 mm。運行時壓力0.1 MPa，流量為32～34 m3/h。管線為碳鋼材質，焊縫寬度為5.3 mm，焊縫系數為0.85。

圖3 水管泄漏聲波檢測與定位裝置

閥門開啟時假定為泄漏狀態，為了探討不同泄漏口徑下，本文所用定位方法的準確性，設置了不同大小的泄漏孔；泄漏孔的直徑分別為3,6,8,10,12,15,20和27 mm。不同的測量距離會影響傳感器接收到的信號強度，測量距離越遠，信號越弱，也就影響了泄漏定位的精度。因此，本研究通過控制變量，在相同的泄漏尺寸下設置了不同的測量距離。

由于聲發射泄漏檢測是一種無損檢測方法，所以傳感器可以沿著外壁管道的任意位置固定，以捕獲聲壓信號[12]。本實驗進行了兩組不同距離的測量，分別對應圖1傳感器1,2,3,4的不同位置，實驗中傳感器位置設置如表1，表中距離為傳感器到泄漏口的距離。

表1 傳感器位置布置表

實驗設備如圖3，在PCI-2系統中，R15a諧振傳感器非侵入性地獲得包含噪聲的泄漏信號，該傳感器采樣率范圍為1～3 000 kHz；利用40 dB的前置放大器放大微弱的輸入信號；18位a/D轉換的聲發射處理卡，即PCI-2采集卡進行數據處理；數據在AEwin軟件中顯示。在閥門開啟，管道內液體流動穩定之后以200 kHz的采樣頻率進行數據采集，每次采集持續時間為25 s，每組工況采集3次數據。

3 實驗結果

3.1 數據處理

由數據采集系統直接獲得的原始AE信號如圖4所示。上述實驗獲得的信號，首先用VMD(variational mode decomposition)方法將噪聲信號從實驗獲得的信號中分離出來，以便獲得噪聲參數；原始信號去除噪聲信號后被認定為純泄漏信號，純泄漏信號被用于計算式(8)、式(10)和式(12)中的幅值。

圖4 數據采集系統所得的原始AE信號圖

VMD在運行前，需要預設二次懲罰項ξ，分解模態分量個數k，噪聲容忍度τ以及收斂準則容差值φ，根據相關計算結果，最終預設值如表2所示。

表2 VMD相關運行參數設置

分解出的7個模態分別標記為u1～u7。隨機選取實驗獲得的聲發射信號作為VMD輸入信號，得到原始信號及VMD分解之后的各模態頻譜圖如圖5所示，各模態有其各自的中心頻率。計算各模態與原始信號之間的相關系數，將相關系數大于0.3的模態進行重構[13]。經過計算，u2～u7與原始信號之間相關系數大于0.3，即對u2～u7重構，記為去噪之后的信號；u1被設定為噪聲信號。

圖5 原始信號頻譜與VMD分解后各模態頻譜圖

分離出的噪聲信號與純泄漏信號需要白化[14]。白化后的噪聲信號用于計算噪聲方差與均值。白化后的純泄漏信號經過特征值提取，首先得到能量值，幅值則用能量的開方獲得。能量值計算公式為

(16)

式中:xi為重構信號;M為每個信號中采樣點數。實驗設定采樣時間為25 s，但是實際操作時，采樣的時間會由于操作人員的反應速度以及軟件反應速度等因素產生微小差異。考慮到不同泄漏口徑進行實驗時，采樣時間的微小差異會影響能量計算結果，進而影響幅值結果，因此在原有能量值的基礎上，將按照采樣時間與能量值的比例，計算5 s采樣時間的能量值。第一組傳感器布置距離在不同泄漏尺寸下聲波信號幅值計算結果在表3中顯示(三次實驗的平均值)；第二組距離在表4中顯示(三次實驗的平均值)。將表3和4的數據畫圖6和圖7，以求更直觀地分析數據。

表3 第一組距離布設下不同泄漏尺寸聲波信號幅值表

表4 第二組距離布設下不同泄漏尺寸聲波信號幅值表

圖6 第一組傳感器布設下不同泄漏口徑幅值對比圖

圖7 第二組傳感器布設下不同泄漏口徑幅值對比圖

圖6中，3,6,8,10,12和15 mm之間幅值差異較小，原因在于泄漏尺寸之間差距較小；而20 mm和27 mm之間幅值差異就比較明顯。隨著泄漏口徑的增加，幅值逐漸增大，12 mm的幅值除外，原因推斷是，在該條件下當泄漏孔小到一個臨界值時，由于出現管路膨脹振動形式[15]，增強了傳感器接收到的聲發射信號。

圖7中，3 mm的信號幅值大于6 mm；8 mm的信號幅值大于10 mm；與第一組距離的規律不同。由于聲波沿管道傳播過程中不同的頻率衰減程度是不同的，也就意味著不同探測距離，傳感器所探測到的信號中占主導的頻率是不同的。在第二組距離中，3 mm信號幅值大于6 mm的原因推斷是3 mm泄漏口處液體汽化，從而發生氣固耦合[15]，增強了AE信號中的某一頻率，而在該探測距離下，這一頻率在聲信號中占主導，所以造成了第二組距離的實驗現象。8 mm信號幅值大于10 mm原因是8 mm管徑在該探測距離下，管路發生膨脹振動形式，且對應的頻率在聲信號中占主導。

本文所用方法只需要確定每一組實驗中不同傳感器之間的信號幅值比，定位結果的準確性主要取決于幅值比值。因此，不同泄漏口徑之間幅值的大小對比不影響結果。

3.2 定位結果

(17)

在兩組不同測量距離下，用基于衰減模型的極大似然估計得到的定位結果如表5。由表5中三組實驗數據得到的誤差平均值做不同泄漏口徑及不同傳感器布設下定位誤差對比，如圖8所示。

表5 不同測量距離處泄漏定位結果

圖8 不同泄漏口徑及不同傳感器布設下定位誤差對比圖

表5中可見，第一組距離，3 mm泄漏孔的實驗1和實驗3所得結果較差，推斷是由于泄漏口徑較小，導致傳感器接收的泄漏信號較弱，被湮沒于噪聲信號中，從而影響幅值計算結果。

由圖8可知，第一組傳感器布設的條件下，3 mm的泄漏口徑誤差最大。而在第二組傳感器布設的條件下，15 mm和20 mm的誤差較大。除此之外，其他預測值誤差都在6%以內，所有實驗數據中最大的誤差也不超過15%，說明該方法能較好地對泄漏點進行定位。也說明在本文的實驗中，傳感器的探測距離對定位結果的影響較小。

該實驗的誤差來源于：

(1)在測量傳感器距泄漏口的距離時，由于測量工具精度以及測量方法產生的誤差。

(2)當泄漏點距離焊縫、法蘭等管道特征較近時，由于反射波的影響，會使定位出現較大誤差[16]。

(3)傳感器的采樣率不夠，導致測量得到的信號不能較完全地反映泄漏情況。

式中，j和l表示φ(jΔx,lΔy)所形成二維矩陣的第j行第L列，Δκx=2π/Dx和Δκx=2π/Dx分別為x方向和y方向頻率間隔，Dx=NxΔx和Dy=NyΔy分別為相位屏兩個方向上的長度，Nx和Ny分別為兩個方向上的網格數，Δx和Δy分別為兩個方向上的網格間距.

(4)在對聲發射信號的處理過程中，使用的VMD去噪方法分離的噪聲與純泄漏信號，與實際的噪聲和泄漏信號有所差別。

(5)極大似然估計法只能估計出最優位置，但是結果可能與實際泄漏位置之間會有誤差。

為了佐證該方法的有效性，對現有設備選取一段較長直管進行試驗。在管段上選取一點進行敲擊，共實驗三次，且每次敲擊強度不同。實驗中敲擊點位于傳感器2與3之間，傳感器與敲擊點之間距離如表6所示。

表6 傳感器位置布置表

實驗所得聲信號經過VMD處理，計算幅值，如表7所示。

表7 敲擊實驗聲波信號幅值表

使用基于衰減模型的極大似然估計所得預測結果如表8。結果表明，該方法在直管段定位能力良好。

表8 敲擊實驗衰減模型定位結果

時差法由于其運算的簡易而被廣泛運用，下文將基于衰減模型的方法與時差法進行比較。

時差法定位的公式為

(18)

式中:L為上下游傳感器之間的距離，此處為3.080 m；v為聲波在管道內傳播速度，此處為4 692 m/s；t為上下游傳感器之間接收信號的時間差，s；x為泄漏點與上游傳感器之間距離，m。

對傳感器1和4得到的敲擊聲信號進行VMD處理，再進行互相關，得到上下游傳感器接收信號的時間差t如表9所示。

表9 敲擊實驗時差法定位結果

由表9看出，時差法所得的定位結果的精確度差距較大，數據處理結果的穩定性較差，且定位誤差大于基于衰減模型的極大似然估計。

4 結 論

時差法作為管道泄漏定位的主要方法，受限于管道中聲速以及時差的確定，所以本文主要研究聲波衰減模型，在原模型的基礎上進行改進，并通過實驗驗證該方法的可靠性。和原衰減模型定位方法相比，該方法的優勢在于：

(1)上下游的衰減參數假定為不同數值，保證了后續公式計算的準確性；

(2)無需確定α在公式(2)中所涉及的參數，α值可直接由實驗獲得；當應用場合變換之后，原衰減模型定位方法需要重新確定α的計算所涉及的參數，本文所提方法避免了定位前對管道運行參數的確定，減少了檢測設備運行前的準備工作。

(3)在原衰減模型的基礎上，增加了噪聲的影響，使得衰減模型更符合實際的管道運行狀況，提高了定位模型的準確性。

本文所得結論如下：

(1)定位結果的準確性表明聲波傳播規律滿足公式(1)指數函數形式。泄漏口徑從3 mm增大到27 mm，聲波幅值逐漸增加，對于幾組異常值，推斷為特定泄漏口徑產生管道膨脹振動形式或者氣固耦合現象，在一定探測距離下，膨脹振動形式或氣固耦合所產生的頻率占主導，導致聲發射信號增強。

(2)采用改進的聲波衰減模型并且結合極大似然估計，所得到的定位方法在不同測量距離與不同泄漏口徑的實驗中，誤差較小；在敲擊實驗中，本文所提模型與時差法比較，基于衰減模型的極大似然估計定位結果優于時差法，說明該方法切實可行。

(3)對于3 mm及以下的泄漏口徑，本文的方法檢測到的泄漏信號較弱，從而導致定位結果誤差較大，所以接下來的研究中應對該方法進行改進，以求對微小泄漏也能準確定位。