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

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

（1.山東大學 信息科學與工程學院，山東 青島 266237；2.山東省激光技術與應用重點實驗室，山東 青島 266237；3.山東大學 激光與紅外系統集成技術教育部重點實驗室，山東 青島 266237）

1 引言

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

基于后向散射型的分布式光纖傳感技術原理是當光纖中的光受到外界擾動時會產生相移，使經過耦合器的光的干涉狀態發生改變。光電探測器可以把光信號轉化成電信號進行信號采集，對采集信號進行處理即可得到擾動信息。然而，基于后向散射光的光纖傳感系統對光源、探測器等硬件要求較高，光路較為復雜，成本較高。同時，受脈沖重復頻率的限制，其在長距離傳感中可測量的頻率響應范圍較低，這也極大地限制了其在實際中的應用。

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

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

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

2 實驗系統及原理

雙M-Z型分布式光纖傳感實驗系統裝置圖如圖1所示，由兩個M-Z干涉儀按照順時針和逆時針方向組合而成。光源發出的光經隔離器到達耦合器C1分為兩束，其中一束光經耦合器C2分為兩束后進入傳感光纖的兩個干涉臂，在耦合器C3處干涉，并經傳導光纖進入耦合器C4被光電探測器PD2接收；另一束光經耦合器C4進入傳導光纖，在經過耦合器C3進入干涉儀的兩個干涉臂,兩束光在耦合器C2處發生干涉，被探測器PD1接收。兩個探測器把探測到的光信號轉為電信號后進入采集處理系統進行定位計算。

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

當有振動作用在M點時，攜帶振動的信號沿光纖順時針和逆時針方向到達兩個探測器存在時間延遲τ：

式中，tL?x為信號從振動位置傳播到耦合器C3的時間，tL為振動信號經過傳導光纖到達探測器PD1的時間，tx為振動信號從振動位置傳到探測器PD2的時間。從而可以計算出振動的位置：

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

由于兩路探測器探測到的振動信號是由同一振動源產生，因此兩路信號具有高度的相關性，通過對兩路信號做互相關運算可以得出時間延遲τ，進而可以得到擾動位置。

分布式光纖傳感技術中空間分辨率是衡量一個系統性能優良的關鍵指標。光干涉后被兩個光電探測器接收，光信號轉化為電信號，后將電信號經過采集卡采樣之后送入計算機進行處理?；诜菍ΨQ雙 M-Z結構的光纖擾動傳感利用兩路干涉信號的時延計算擾動位置。因此，采集卡的采樣頻率直接決定著傳感器的理論定位精度。假設采樣頻率為fs，此時相鄰采樣數據的時間間隔為?=1/fs， 任意時延均可以寫成d=N?t的形式（N為正整數）。因此可以得到理想條件下傳感器的空間分辨率為：

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

3 希爾伯特-黃變換算法

希爾伯特-黃變換包含經驗模態分解（Empirical Mode Decomposition，EMD）和希爾伯特譜分析（Hilbert Spectrum Analysis，HSA）[14-18]。

3.1 經驗模態分解

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

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

IMF滿足以下兩個條件：(1) 信號極值點的數量與零點數相等或者相差1；(2) 由極大值定義的上包絡和極小值定義的下包絡的局部均值是零。

作為燒結配料：鋼渣主要化學成分為CaO、FeO等，當CaO含量較高時，鋼渣可作燒結礦助熔劑替代部分石灰熔劑，適量配入鋼渣可改善燒結礦質量，使轉鼓指數和燒結率提高，并且由于鋼渣中FeO的氧化放熱，降低燒結礦燃耗。

3.2 希爾伯特譜分析

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

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

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

相位函數求導得到信號的瞬時頻率函數：

對所有IMF進行Hilbert變換，得：

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

4 實驗過程

搭建如圖1所示的實驗系統，采用1550 nm窄線寬激光器(NLL)作為傳感光源，經過隔離器(ISO)后進入50∶50耦合器分束，分別沿順時針和逆時針方向傳播，各自發生干涉后被光電探測器(PD)采集，進入采集處理系統進行數據分析。信號采集使用NI PCI-5114采集卡，通過LabVIEW采集信號，信號處理使用Matlab進行定位計算。采用的光纖為2 km雙芯鎧裝光纜，抗環境噪聲的能力較強，同時更適合實際應用。

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

圖2 干涉信號時域圖Fig.2 Time domain diagram of interference signal

該信號處理方法流程圖如圖3所示，具體運算如下：

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

（1）對振動光纖信號按分段3次Hermite差值多項式做EMD，得到IMF，分解層數選取10層，如圖4所示，這里給出了5層分解。

圖4 干涉信號EMD分解圖Fig.4 EMD decomposition diagram of interference signal

（2）對分解得到的每層IMF做Hilbert變換，結果反映了原始信號的瞬時頻率。

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

（4）根據Hilbert譜峰值位置提取原始信號的信號段，取最高峰值處附近0.01 s的數據做互相關運算及定位運算得出振動位置信息。從圖5可以看出，880000采樣點處頻率最高，選取800000～900000采樣點處的區域信號作為振動的有效信號，對此做互相關運算得出兩路信號的時間延遲,通過計算得出定位為1810 m，與實際振動位置1800 m相比，誤差為10 m。

圖5 干涉信號的希爾伯特譜圖Fig.5 Hilbert spectrum of interference signal

為了說明該方法的有效性，分別取不同頻率的區域信號做定位運算，比較定位誤差的大小，如表1所示。通過對不同區域信號定位誤差的對比可知，本文所提出的定位方法比傳統的直接互相關方法定位誤差小，同時利用該方法可精確定位出最優振動區域信號，取高頻段的信號定位誤差比低頻段信號定位誤差小。

表1 不同區域信號定位誤差對比Tab.1 Comparison of signal positioning errors in different regions

5 結論

本文利用HHT來確定振動的區域信號位置，從而提取出包含振動的有效區域信號，振動區域為頻率最高的信號區域，對該區域信號進行互相關運算，得出兩路信號的時間延遲，再進行定位計算，減少了運算量，并且極大地提高了定位精度，在2 km的傳感長度、10 MHz采樣率下，定位精度達10 m。