基于希爾伯特-黃變換的雙馬赫-曾德分布式光纖傳感振動(dòng)定位方法

楊文晨，秦增光，劉兆軍，徐演平，李 釗，渠 帥，叢振華，王澤群

（1.山東大學(xué) 信息科學(xué)與工程學(xué)院，山東 青島 266237；2.山東省激光技術(shù)與應(yīng)用重點(diǎn)實(shí)驗(yàn)室，山東 青島 266237；3.山東大學(xué) 激光與紅外系統(tǒng)集成技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，山東 青島 266237）

1 引言

近年來(lái)，光纖通信技術(shù)發(fā)展迅速并且得到廣泛應(yīng)用，光纖傳感技術(shù)也受到了廣泛的關(guān)注，在越來(lái)越多的領(lǐng)域得以應(yīng)用，大型建筑、橋梁、管道的安全監(jiān)測(cè)，飛機(jī)、高鐵的運(yùn)行安全監(jiān)測(cè)，電纜的溫度監(jiān)測(cè)，周界安全的防護(hù)監(jiān)測(cè)，海底地質(zhì)災(zāi)害預(yù)測(cè)，甚至石油產(chǎn)業(yè)、電力行業(yè)、機(jī)械加工及醫(yī)療行業(yè)等領(lǐng)域都與其息息相關(guān)[1-6]。作為光纖傳感領(lǐng)域最重要的一類，分布式光纖傳感主要分為基于后向散射光的分布式光纖傳感系統(tǒng)和基于干涉儀結(jié)構(gòu)的分布式光纖傳感系統(tǒng)?；诤笙蛏⑸涞姆植际焦饫w傳感系統(tǒng)又根據(jù)后向散射光的不同，可以分為拉曼光時(shí)域反射計(jì)（Raman Optical Time Domain Reflectometer，R-OTDR）、布里淵光時(shí)域反射計(jì)（Brillouin Optical Time Domain Reflectometer，B-OTDR）和相位敏感型光時(shí)域反射計(jì)（Phasesensitive Optical Time Domain Reflectometer， φ-OTDR）?；诟缮鎯x結(jié)構(gòu)的分布式光纖傳感系統(tǒng)根據(jù)其裝置的不同可分為薩格納克（Sagnac）干涉儀、馬赫-曾德（Mach-Zehnder，M-Z）干涉儀和邁克爾遜（Michelson）干涉儀。

基于后向散射型的分布式光纖傳感技術(shù)原理是當(dāng)光纖中的光受到外界擾動(dòng)時(shí)會(huì)產(chǎn)生相移，使經(jīng)過(guò)耦合器的光的干涉狀態(tài)發(fā)生改變。光電探測(cè)器可以把光信號(hào)轉(zhuǎn)化成電信號(hào)進(jìn)行信號(hào)采集，對(duì)采集信號(hào)進(jìn)行處理即可得到擾動(dòng)信息。然而，基于后向散射光的光纖傳感系統(tǒng)對(duì)光源、探測(cè)器等硬件要求較高，光路較為復(fù)雜，成本較高。同時(shí)，受脈沖重復(fù)頻率的限制，其在長(zhǎng)距離傳感中可測(cè)量的頻率響應(yīng)范圍較低，這也極大地限制了其在實(shí)際中的應(yīng)用。

干涉型分布式光纖傳感系統(tǒng)可以彌補(bǔ)后向散射型分布式光纖傳感系統(tǒng)的缺點(diǎn)。同時(shí)，干涉型分布式光纖傳感系統(tǒng)具有頻率響應(yīng)寬、對(duì)光源要求較低、光路較為簡(jiǎn)單、成本低等優(yōu)勢(shì)。其中，基于雙M-Z結(jié)構(gòu)的分布式光纖傳感系統(tǒng)因其相位敏感性高、響應(yīng)速度快、定位精確，逐漸成為分布式光纖傳感系統(tǒng)的研究熱點(diǎn)，并已在周界安全防范、石油管道安全預(yù)警、輸電線網(wǎng)安全監(jiān)控及海底線纜監(jiān)控等領(lǐng)域有廣泛的應(yīng)用[7]。國(guó)內(nèi)外有關(guān)公司已對(duì)定位系統(tǒng)產(chǎn)品進(jìn)行設(shè)計(jì)及應(yīng)用，但由于定位精度低，很難進(jìn)行大面積推廣。隨著雙MZ應(yīng)用需求的逐漸增加，對(duì)其定位精度也有了越來(lái)越高的要求。

傳統(tǒng)的雙M-Z振動(dòng)定位方法是直接對(duì)信號(hào)做互相關(guān)運(yùn)算，通過(guò)計(jì)算時(shí)延得到振動(dòng)位置，陳偉民等使用傳統(tǒng)M-Z振動(dòng)定位方法在長(zhǎng)20 km光纜、10 MHz采樣率下進(jìn)行計(jì)算，定位誤差為149 m[8]。直接互相關(guān)計(jì)算會(huì)把大量的噪聲信息也計(jì)算其中，使得振動(dòng)定位精度下降。此外，在信號(hào)中難以確定振動(dòng)起始點(diǎn)也會(huì)對(duì)精度和處理時(shí)間帶來(lái)不利的影響[9-11]。為了提高雙M-Z振動(dòng)系統(tǒng)的定位精度，陳沁楠等采用過(guò)零率得到擾動(dòng)的起始點(diǎn)，并通過(guò)互相關(guān)算法處理，在2.25 km光纖實(shí)現(xiàn)了±20 m定位精度[7]。該方法的準(zhǔn)確度過(guò)于依賴于對(duì)過(guò)零窗口大小的選擇，窗口選擇較大會(huì)引入噪聲誤差，反之則會(huì)丟失重要信息。馬春宇等采用雙波長(zhǎng)M-Z系統(tǒng)結(jié)合連續(xù)小波變換技術(shù)通過(guò)提取采集信號(hào)包絡(luò)尋找瞬時(shí)頻率較高部分以減小定位誤差，在85 km的傳感光纖上平均誤差為24.3 m[12]。謝尚然等采用離散小波尋找振動(dòng)起始點(diǎn)進(jìn)行定位，在長(zhǎng)度為30.498 km、2 MHz采樣率下，定位誤差達(dá)50 m[13]。然而由于小波分析的結(jié)果主要取決于小波變換的基本功能，無(wú)法以相同的精度對(duì)時(shí)間或頻率相關(guān)信息進(jìn)行準(zhǔn)確分析。因此尋找有效的檢測(cè)方法不僅可以提高定位精度，而且可以減少處理時(shí)間,對(duì)于干涉行儀分布式光纖傳感具有重要意義。上述研究的定位精度還需要提高，限制了其應(yīng)用。

為了實(shí)現(xiàn)高精度的振動(dòng)定位，本文提出了一種基于希爾伯特-黃變換（Hilbert-Huang Transform，HHT）的雙M-Z型分布式光纖傳感振動(dòng)定位方法。通過(guò)HHT方法，可以準(zhǔn)確判斷瞬時(shí)高頻率信息，從而提高互相關(guān)計(jì)算的準(zhǔn)確度，減小雙M-Z傳感系統(tǒng)誤差。

2 實(shí)驗(yàn)系統(tǒng)及原理

雙M-Z型分布式光纖傳感實(shí)驗(yàn)系統(tǒng)裝置圖如圖1所示，由兩個(gè)M-Z干涉儀按照順時(shí)針和逆時(shí)針?lè)较蚪M合而成。光源發(fā)出的光經(jīng)隔離器到達(dá)耦合器C1分為兩束，其中一束光經(jīng)耦合器C2分為兩束后進(jìn)入傳感光纖的兩個(gè)干涉臂，在耦合器C3處干涉，并經(jīng)傳導(dǎo)光纖進(jìn)入耦合器C4被光電探測(cè)器PD2接收；另一束光經(jīng)耦合器C4進(jìn)入傳導(dǎo)光纖，在經(jīng)過(guò)耦合器C3進(jìn)入干涉儀的兩個(gè)干涉臂,兩束光在耦合器C2處發(fā)生干涉，被探測(cè)器PD1接收。兩個(gè)探測(cè)器把探測(cè)到的光信號(hào)轉(zhuǎn)為電信號(hào)后進(jìn)入采集處理系統(tǒng)進(jìn)行定位計(jì)算。

圖1 雙M-Z分布式光纖傳感系統(tǒng)原理圖Fig.1 Schematic diagram of dual M-Z distributed optical fiber sensor

當(dāng)有振動(dòng)作用在M點(diǎn)時(shí)，攜帶振動(dòng)的信號(hào)沿光纖順時(shí)針和逆時(shí)針?lè)较虻竭_(dá)兩個(gè)探測(cè)器存在時(shí)間延遲τ：

式中，tL?x為信號(hào)從振動(dòng)位置傳播到耦合器C3的時(shí)間，tL為振動(dòng)信號(hào)經(jīng)過(guò)傳導(dǎo)光纖到達(dá)探測(cè)器PD1的時(shí)間，tx為振動(dòng)信號(hào)從振動(dòng)位置傳到探測(cè)器PD2的時(shí)間。從而可以計(jì)算出振動(dòng)的位置：

式中，L為傳感光纖和傳導(dǎo)光纖的長(zhǎng)度，τ為順逆兩束光到達(dá)平衡探測(cè)器間的時(shí)間延遲，c為真空中的光速，n為光纖折射率，x為振動(dòng)位置M到耦合器C1的距離。

由于兩路探測(cè)器探測(cè)到的振動(dòng)信號(hào)是由同一振動(dòng)源產(chǎn)生，因此兩路信號(hào)具有高度的相關(guān)性，通過(guò)對(duì)兩路信號(hào)做互相關(guān)運(yùn)算可以得出時(shí)間延遲τ，進(jìn)而可以得到擾動(dòng)位置。

分布式光纖傳感技術(shù)中空間分辨率是衡量一個(gè)系統(tǒng)性能優(yōu)良的關(guān)鍵指標(biāo)。光干涉后被兩個(gè)光電探測(cè)器接收，光信號(hào)轉(zhuǎn)化為電信號(hào)，后將電信號(hào)經(jīng)過(guò)采集卡采樣之后送入計(jì)算機(jī)進(jìn)行處理?；诜菍?duì)稱雙 M-Z結(jié)構(gòu)的光纖擾動(dòng)傳感利用兩路干涉信號(hào)的時(shí)延計(jì)算擾動(dòng)位置。因此，采集卡的采樣頻率直接決定著傳感器的理論定位精度。假設(shè)采樣頻率為fs，此時(shí)相鄰采樣數(shù)據(jù)的時(shí)間間隔為?=1/fs， 任意時(shí)延均可以寫(xiě)成d=N?t的形式（N為正整數(shù)）。因此可以得到理想條件下傳感器的空間分辨率為：

由于可以認(rèn)為光纖折射率是不變的，因此定位空間分辨率只與采樣率有關(guān)，采樣率越高，理論空間分辨率就越高。本文使用的采集卡的采樣率為10 MS/s，因此，本文的理論空間分辨率為10 m。

3 希爾伯特-黃變換算法

希爾伯特-黃變換包含經(jīng)驗(yàn)?zāi)B(tài)分解（Empirical Mode Decomposition，EMD）和希爾伯特譜分析（Hilbert Spectrum Analysis，HSA）[14-18]。

3.1 經(jīng)驗(yàn)?zāi)B(tài)分解

EMD的基本思想是把一個(gè)頻率變化無(wú)規(guī)律的波轉(zhuǎn)化成多個(gè)單一頻率波+殘余波的形式。具體是將一個(gè)原始信號(hào)x(t)的極大值點(diǎn)用3次Hermite差值方式擬合出上包絡(luò)線，再用同樣的方式得出下包絡(luò)線；上下包絡(luò)線取均值，得到平均包絡(luò)線，記為pl；用原始信號(hào)x(t)減去平均包絡(luò)線pl，得到一個(gè)新的信號(hào)m(t)，即：

若新信號(hào)m(t)還存在負(fù)的局部極大值和正的局部極小值，則繼續(xù)按上述操作分解，由此原始信號(hào)x(t)分解為n個(gè)本征模函數(shù)（Intrinsic Mode Function，IMF）和余波r(t)，可表示為：

IMF滿足以下兩個(gè)條件：(1) 信號(hào)極值點(diǎn)的數(shù)量與零點(diǎn)數(shù)相等或者相差1；(2) 由極大值定義的上包絡(luò)和極小值定義的下包絡(luò)的局部均值是零。

作為燒結(jié)配料：鋼渣主要化學(xué)成分為CaO、FeO等，當(dāng)CaO含量較高時(shí)，鋼渣可作燒結(jié)礦助熔劑替代部分石灰熔劑，適量配入鋼渣可改善燒結(jié)礦質(zhì)量，使轉(zhuǎn)鼓指數(shù)和燒結(jié)率提高，并且由于鋼渣中FeO的氧化放熱，降低燒結(jié)礦燃耗。

3.2 希爾伯特譜分析

信號(hào)x(t)的希爾伯特變換y(t)表示為：

式中，p代表柯西主值，由此確定的解析信號(hào)為：

式中，a(t)為瞬時(shí)幅值， θ(t)為瞬時(shí)相位，由式(8)決定。

相位函數(shù)求導(dǎo)得到信號(hào)的瞬時(shí)頻率函數(shù)：

對(duì)所有IMF進(jìn)行Hilbert變換，得：

這里省略了余波r(t)，式中，H(f,t)是隨時(shí)間和頻率變化的信號(hào)幅度。定義Hilbert譜：

4 實(shí)驗(yàn)過(guò)程

搭建如圖1所示的實(shí)驗(yàn)系統(tǒng)，采用1550 nm窄線寬激光器(NLL)作為傳感光源，經(jīng)過(guò)隔離器(ISO)后進(jìn)入50∶50耦合器分束，分別沿順時(shí)針和逆時(shí)針?lè)较騻鞑ィ髯园l(fā)生干涉后被光電探測(cè)器(PD)采集，進(jìn)入采集處理系統(tǒng)進(jìn)行數(shù)據(jù)分析。信號(hào)采集使用NI PCI-5114采集卡，通過(guò)LabVIEW采集信號(hào)，信號(hào)處理使用Matlab進(jìn)行定位計(jì)算。采用的光纖為2 km雙芯鎧裝光纜，抗環(huán)境噪聲的能力較強(qiáng)，同時(shí)更適合實(shí)際應(yīng)用。

實(shí)驗(yàn)中，在距離C2耦合器1800 m處用腳踩光纖的方式施加振動(dòng)，采樣率選用10 MHz/s，采樣時(shí)間為0.1 s，振動(dòng)時(shí)間選取0.01 s。采集到的信號(hào)如圖2（彩圖見(jiàn)期刊電子版）所示，可見(jiàn)兩路信號(hào)具有很強(qiáng)的相關(guān)性。由于光纖的雙折射因素、探測(cè)器增益帶寬和光源噪聲造成了兩路信號(hào)大小的不一致。由于兩路信號(hào)反映的都是同一振動(dòng)源的振動(dòng)，因此兩路信號(hào)形狀上具有高度相關(guān)性，通過(guò)互相關(guān)計(jì)算可以得到兩路相關(guān)信號(hào)的時(shí)間延遲，這個(gè)結(jié)果主要考慮兩路信號(hào)的形狀相似度，信號(hào)幅值的大小對(duì)其影響較小。雙M-Z的定位結(jié)果與得到的時(shí)間延遲有關(guān)，因此還是可以通過(guò)式（2）對(duì)振動(dòng)位置最終進(jìn)行定位。

圖2 干涉信號(hào)時(shí)域圖Fig.2 Time domain diagram of interference signal

該信號(hào)處理方法流程圖如圖3所示，具體運(yùn)算如下：

圖3 定位計(jì)算方法流程圖Fig.3 Flow chart of positioning calculation method

（1）對(duì)振動(dòng)光纖信號(hào)按分段3次Hermite差值多項(xiàng)式做EMD，得到IMF，分解層數(shù)選取10層，如圖4所示，這里給出了5層分解。

圖4 干涉信號(hào)EMD分解圖Fig.4 EMD decomposition diagram of interference signal

（2）對(duì)分解得到的每層IMF做Hilbert變換，結(jié)果反映了原始信號(hào)的瞬時(shí)頻率。

（3）將所有IMF做Hilbert變換后的結(jié)果進(jìn)行疊加，得到Hilbert譜，清晰直觀地反映原始信號(hào)的頻率變化情況，如圖5所示。

（4）根據(jù)Hilbert譜峰值位置提取原始信號(hào)的信號(hào)段，取最高峰值處附近0.01 s的數(shù)據(jù)做互相關(guān)運(yùn)算及定位運(yùn)算得出振動(dòng)位置信息。從圖5可以看出，880000采樣點(diǎn)處頻率最高，選取800000～900000采樣點(diǎn)處的區(qū)域信號(hào)作為振動(dòng)的有效信號(hào)，對(duì)此做互相關(guān)運(yùn)算得出兩路信號(hào)的時(shí)間延遲,通過(guò)計(jì)算得出定位為1810 m，與實(shí)際振動(dòng)位置1800 m相比，誤差為10 m。

圖5 干涉信號(hào)的希爾伯特譜圖Fig.5 Hilbert spectrum of interference signal

為了說(shuō)明該方法的有效性，分別取不同頻率的區(qū)域信號(hào)做定位運(yùn)算，比較定位誤差的大小，如表1所示。通過(guò)對(duì)不同區(qū)域信號(hào)定位誤差的對(duì)比可知，本文所提出的定位方法比傳統(tǒng)的直接互相關(guān)方法定位誤差小，同時(shí)利用該方法可精確定位出最優(yōu)振動(dòng)區(qū)域信號(hào)，取高頻段的信號(hào)定位誤差比低頻段信號(hào)定位誤差小。

表1 不同區(qū)域信號(hào)定位誤差對(duì)比Tab.1 Comparison of signal positioning errors in different regions

5 結(jié)論

本文利用HHT來(lái)確定振動(dòng)的區(qū)域信號(hào)位置，從而提取出包含振動(dòng)的有效區(qū)域信號(hào)，振動(dòng)區(qū)域?yàn)轭l率最高的信號(hào)區(qū)域，對(duì)該區(qū)域信號(hào)進(jìn)行互相關(guān)運(yùn)算，得出兩路信號(hào)的時(shí)間延遲，再進(jìn)行定位計(jì)算，減少了運(yùn)算量，并且極大地提高了定位精度，在2 km的傳感長(zhǎng)度、10 MHz采樣率下，定位精度達(dá)10 m。