基于希尔伯特-黄变换的双马赫-曾德分布式光纤传感振动定位方法

杨文晨，秦增光，刘兆军，徐演平，李 钊，渠 帅，丛振华，王泽群

（1.山东大学 信息科学与工程学院，山东 青岛 266237；2.山东省激光技术与应用重点实验室，山东 青岛 266237；3.山东大学 激光与红外系统集成技术教育部重点实验室，山东 青岛 266237）

1 引言

近年来，光纤通信技术发展迅速并且得到广泛应用，光纤传感技术也受到了广泛的关注，在越来越多的领域得以应用，大型建筑、桥梁、管道的安全监测，飞机、高铁的运行安全监测，电缆的温度监测，周界安全的防护监测，海底地质灾害预测，甚至石油产业、电力行业、机械加工及医疗行业等领域都与其息息相关[1-6]。作为光纤传感领域最重要的一类，分布式光纤传感主要分为基于后向散射光的分布式光纤传感系统和基于干涉仪结构的分布式光纤传感系统。基于后向散射的分布式光纤传感系统又根据后向散射光的不同，可以分为拉曼光时域反射计（Raman Optical Time Domain Reflectometer，R-OTDR）、布里渊光时域反射计（Brillouin Optical Time Domain Reflectometer，B-OTDR）和相位敏感型光时域反射计（Phasesensitive Optical Time Domain Reflectometer， φ-OTDR）。基于干涉仪结构的分布式光纤传感系统根据其装置的不同可分为萨格纳克（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的距离。

由于两路探测器探测到的振动信号是由同一振动源产生，因此两路信号具有高度的相关性，通过对两路信号做互相关运算可以得出时间延迟τ，进而可以得到扰动位置。

分布式光纤传感技术中空间分辨率是衡量一个系统性能优良的关键指标。光干涉后被两个光电探测器接收，光信号转化为电信号，后将电信号经过采集卡采样之后送入计算机进行处理。基于非对称双 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。