张令春 姜海明 张俊喜 谢康

1) (合肥工业大学仪器科学与光电工程学院,合肥 230009)

2) (广东工业大学,省部共建精密电子制造技术与装备国家重点实验室,广州 510006)

3) (广东工业大学机电工程学院,广州 510006)

已有的椭圆拟合算法利用李萨如图形解决了非理想3×3 耦合器所引入的解调结果误差,但在研究中没有充分考虑电路噪声与相位噪声对李萨如图形的影响.针对噪声影响提出了基于局部异常因子优化的3×3 耦合器椭圆拟合解调方法,可以有效地消除李萨如图形中因噪声产生的离群点.仿真结果表明,该方法能够实现准确解调,利用蒙特卡洛分析得出综合误差率约为0.13%,验证了方法的稳定性.与传统最小二乘拟合方法相比,所提方法有效地提高了信号解调的稳定性和精确度.

1 背景

随着光纤通信技术的发展,光纤传感技术也成为国内外研究的热点[1-7].光纤传感器因其稳定性好、可靠性高、对电磁波的不敏感性等特点,广泛应用于高温、腐蚀性或危险性的环境中[8,9].信号解调是光纤传感器应用中的一个重要研究方向,而干涉法解调以其高分辨率的独特优势深受青睐,是最常用到的解调方法之一.其中,基于3×3 耦合器的解调算法对激光源的要求较低,稳定性高[10],在实践中得到广泛应用.常用的基于3×3 耦合器解调方法有反正切法、微分交叉相乘法(differential cross multiplication,DCM)[11]和椭圆拟合估计法等[12].微分交叉相乘法要求三路输出信号相差120°,但在实际应用中3×3 耦合器无法达到完全对称,各路输出功率有所差别,加上插入损耗等因素,需要对信号进行归一化处理[13],增加了算法难度.椭圆拟合估计法相比于DCM 不要求耦合器完全对称,增强了系统的适应性;无需对信号直流补偿,简化了算法.

Liu 等[14]利用最小二乘法对每两路输出信号构成的李萨如图形进行椭圆拟合反解出交直流项系数,进而求解相位变化量.王伟等[15]利用双参数椭圆拟合算法实现了法布里-珀罗腔长的解调.梅泽等[16]聚焦弱信号解调,改进后的算法可以准确解调弱振动信号,提升了探测分辨率.但以上研究只考虑了3×3 耦合器两路信号不完全相差120°这一情况,未考虑实际解调中电路噪声和相位噪声的干扰,噪声干扰会直接影响椭圆拟合算法中的核心—李萨如图形,从而影响解调精度.高晓文为降低信号畸变提出的改进方案中考虑了噪声干扰这一实际因素[17],但设定的干扰信号极其微弱,几乎相当于理想状态.综合考虑3×3 耦合器两路信号不完全相差120°以及实际中噪声干扰的问题,本文提出了一种基于局部异常因子算法(local outlier factor,LOF)的3×3 耦合器椭圆拟合解调方法,该方法相较于传统直接利用最小二乘拟合的方法,能够有效地降低电路噪声和相位噪声对椭圆拟合解调的影响,提高解调精度.

2 系统光路

光纤振动传感中基于3×3 耦合器的迈克尔逊干涉仪结构如图1 所示.它由环形器、3×3 耦合器、压电陶瓷(piezoelectric transducer,PZT)和法拉第反射镜组成.光信号通过环形器1—2 端口,首先经3×3 耦合器的两臂分为2 束光进入干涉仪,然后分别被法拉第反射镜反射之后,在 3×3 耦合器重新会合发生干涉,最终包含振动相位调制信息的光信号分别由3 个光电探测器接收.引入法拉第反射镜的目的是消除偏振态对干涉信号的影响.

图1 基于3×3 耦合器解调的迈克尔逊干涉仪结构Fig.1.The structure of Michelson interferometer based on 3×3 coupler.

3 椭圆拟合解调

3.1 基于椭圆拟合的3×3 光纤耦合器的解调原理

对于理想的3×3 耦合器,假设条纹可见性相等,则耦合器三路输出信号中的两路经光电转换后为[18]

其中θ=120°;I1,I2与耦合器性能有关.

实际工程应用中,理想的3×3耦合器并不存在,3×3 耦合器输出信号的条纹可见性不相等,信号之间的相位差θ也不稳定在 120°,因此耦合器的输出为:

检测信号系数A1,A2,B1,B2,C2的求解建立在椭圆拟合的基础上,椭圆合成遵循李萨如图的性质和规律.由(2)式可知,耦合器输出信号u1和u2对同一相位函数φ进行调制,其变化形式与简谐振动一致,变化频率相同,相位差θ′恒定,符合李萨如图的形成条件.此外,由于两路信号的相位差接近 120°,如果在互相垂直的两个坐标轴上同时画出u1,u2,其轨迹能形成一个椭圆,椭圆方程可表示为

式中,g=[a,b,c,d,e,f]T为椭圆的系数矢量.对于(4)式所示的二次曲线,当且仅当b2−4ac ＜0 时,曲线为椭圆.为方便计算,选取b2−4ac=0,采用一种适当的椭圆拟合方法解得椭圆方程系数,再与信号系数对应,计算出信号系数,传统主要采用最小二乘拟合.信号系数计算公式为

将(5)式解得的信号系数代入信号方程(3),可以获得 cosφ() 与 sinφ() 信号.采用 3×3 耦合器的微分交叉相乘还原算法,解得相位变化:

3.2 噪声对于椭圆拟合的影响

在实际解调中,噪声的干扰会影响最终解调结果.迈克尔逊干涉仪的噪声来源主要分为两个方面: 一是激光器的强度噪声和频率噪声;二是解调系统中通过解调干涉仪引入的相位噪声、环境噪声、信号采集电路中光电探测器引起的噪声等.这些噪声在所有基于3×3 光纤耦合器的解调结构中都存在,研究的泛用性较广.

系统解调过程中,输出信号经过光电转换变为电流信号,再通过电路转换成电压信号,这一过程不可避免地会引入电路噪声,电路噪声会形成随机噪声,随机噪声还与光电探测器性能和环境因素有关,而在解调过程中还存在相位偏移量即相位噪声,相位噪声与 3×3 耦合器的实际性能有关.这些噪声都会对解调产生不利影响.假设相位噪声nε(t)和加性电路噪声nc(t) 都用正态分布随机噪声模拟,两路信号可表示为

为了验证噪声对干涉信号的影响,在仿真中信号幅值的单位选为电压幅值,假设直流量A1,A2分别为和的值分别为9 V 和 6 V ;分别用幅值为0.5 V 和 1.5 V的随机噪声来模拟相位噪声nε(t) 和加性电路噪声nc(t)[19,20],在传感光纤上叠加扰动信号:

令仿真系统的采样频率为fs=20000Hz,采样数为N=10000,则两路仿真信号如图2 所示.

图2 两路仿真信号 (a) 无噪声干扰信号;(b) 有噪声干扰信号Fig.2.Two simulated signals: (a) Without noise interference;(b) with noise interference.

图2 显示了(7)式表示的两路信号,其中绿色曲线代表信号u1,蓝色曲线代表u2.通过图2(a)和图2(b)两图的对比可以看出,噪声对两路信号均有明显干扰.而两路信号形成的李萨如图是椭圆拟合的基础,两路信号受噪声干扰,由两路信号形成的李萨如图也势必受到影响.两路信号的李萨如图形如图3 所示.

从图3 可以看出,噪声产生的离群点在理想椭圆曲线周围不规则分布,椭圆内外轮廓的离群点使得椭圆边界变 “厚”,导致很难直接利用最小二乘等数值估计方法准确获取椭圆表达式,对比可以看出,最小二乘法拟合的椭圆曲线椭圆与理想椭圆曲线相比倾角有明显偏差.

图3 两路信号的李萨如图形及拟合曲线Fig.3.Lissajous figure and fitted curves of two signals.

4 利用局部异常因子算法(LOF)对椭圆拟合进行优化

4.1 LOF 算法原理

局部离群因子(local outlier factor)检测算法是一种无监督的异常检测方法,是基于密度的局部离群点检测算法[21].算法内核是比较目标点与相邻点的局部密度偏差,通过计算局部密度偏差的比值来判断目标点是否为离群点.由于算法在计算目标点的同时计算了相邻点,目标点的局部密度偏差比值包含了邻近点的信息,算法的精度较高.LOF算法的相关定义和计算方法如下.

4.1.1 样本点之间的距离

对于样本集合中的n 个监测样本,任意的两个样本数据点之间的距离可以采用欧氏(Eucildean)距离、汉明(Hamming)距离、马氏(Mahalanobis)距离、切比雪夫(Chebyshev)距离、闵可夫斯基(Minkowski)距离和球面距离等.其中,闵可夫斯基距离公式如(9)式所示:

当(9)式中q为2 时,即为欧式距离.

4.1.2 第k距离

如果数据集中至少有不包括p在内的k个数据点oi ∈C{x/=p},满足d(p,o′)≤d(p,o),并且至多有不包括p在内的k −1 个点o′ ∈C{x/=p},满 足d(p,o′)≤d(p,o),即数据点p的第k距离 dk(p) 就是指数据集中距离p第k远的点与p之间的距离.

4.1.3 第k距离领域

数据点p的第k距离领域Nk(p) 就是指在数据集中所有与p点之间的距离小于或等于 dk(p)的数据点的集合.所以Nk(p) 中的元素个数大于或等于k.

4.1.4 可达距离

如果点o是远离点p的,那么点o到点p的第k可达距离就是o,p间的真实距离,否则可达距离用o的第k距离表示(见图4).也就是说,o到距离它最近的k个点的可达距离可以认为是相等的,并且都等于dk(o).所以可达距离的公式定义是:

图4 k=5 时,对象 p,o1的可达 距离和对象 p,o2的可达距离的示意图Fig.4.The schematic diagram of the reach-distances of p,o1and p,o2 when k=5.

4.1.5 局部可达密度

点p的局部可达密度公式可以表示为

局部可达密度是指点p的第k邻域内的点到p的平均可达距离的倒数.局部可达密度越高则越有可能属于正常数据点,局部可达密度越低则越有可能是离群点.如果样本点与其近邻域的点不属于同一个聚簇,那么它的可达距离就是数值较大的d(p,o1),所以局部可达密度更小,但如果样本点与其近邻域的点属于同一个聚簇,那么可达距离就是值较小的dk(o),所以局部可达密度会更大.

4.1.6 局部异常因子

点p的局部异常因子公式可以表示为

点p的局部异常因子表示点p的领域点Nk(p)的局部可达密度与点p的局部可达密度之比的平均数.如果局部异常因子小于1,则表示点p的局部可达密度高于它的领域点的局部可达密度,那么点p属于密集点;如果局部异常因子等于1,则表示点p的局部可达密度与它的领域点的局部可达密度差不多,那么点p与它的邻近点可能属于同一聚簇;如果局部异常因子大于1,则表示点p的局部可达密度低于它的领域点的局部可达密度,那么点p属于离群点.

4.2 局部异常因子算法(LOF)对椭圆拟合进行优化

3.2 节讨论了噪声对于椭圆拟合的影响,指出由于噪声影响无法直接利用最小二乘等数值估计方法来准确获取椭圆表达式.由图3 可以看出,噪声会在标准拟合曲线内外侧形成无规则的离散点,利用LOF 算法可以判断这些离散点哪些为离群点,并且将这些离群点剔除.将离群点剔除后,可以使得李萨如图形趋向于无噪声干扰时所合成的椭圆.从而可以较准确地获取有效信号的李萨如图形椭圆表达式,再反解椭圆系数,从而消除电路噪声和相位偏移对信号解调的影响.LOF 算法处理李萨如图形中离群点可以看作图像处理类方法,所以首先要将两路信号生成的图像矩阵转化为灰度图像再转化为二值图,通过LOF 算法处理减少噪声的影响.最后利用最小二乘拟合求出椭圆系数并利用微分交叉相乘方法(DCM)完成扰动信号的相位还原.系统的总体设计流程图如图5 所示.

图5 系统流程图Fig.5.System flowchart.

5 仿真验证

利用数值仿真进行验证,按照3.2 小节部分设置信号参数,A1,A2分别为10V和12V,B1′,B2′的值分别为9 V 和 6 V ;分别用幅值为0.5 V 和1.5 V的随机噪声来模拟相位噪声nε(t) 和加性电路噪声nc(t),扰动信号设置为:φ(t)=5+sin(100πt) .两路信号形成的李萨如图像进行灰度化和二值化处理,再经LOF 算法处理,处理结果如图6(b)所示,处理前的李萨如图像如图6(a)所示.

图6 LOF 算法处理前后的李萨如图像 (a) 处理前的李萨如图像;(b)处理后的李萨如图像Fig.6.Lissajous figures with/without LOF processing: (a) Before processing;(b) after processing.

通过LOF 算法处理前后对比可以看出,经过LOF 算法处理后的李萨如图像离群点得到有效去除,图像左下部离群点较少,椭圆轮廓明显“收窄”,从而有利于后续椭圆拟合.经过LOF 算法处理并采用微分交叉相乘法得到的解调后的完整信号如图7(a)所示,其中红色曲线代表原始扰动信号即(8)式所代表的信号,绿色曲线代表利用最小二乘法得到的解调信号,蓝色曲线是利用LOF 算法得到的解调信号.为清晰起见,图7(b)给出了图7(a)中红色虚线框部分的放大图.

图7 解调信号 (a) 完整解调信号;(b) 解调信号局部Fig.7.The demodulated signals: (a) Full signals;(b) partially-amplified signals.

由图7 可见,经LOF 处理的椭圆拟合解调法与直接利用最小二乘算法的解调法相比,前者所得的解调信号更接近原始信号.故该方法可以较好地还原扰动信号,降低解调信号的误差.对比其他改进最小二乘解调法,可以看到文献[20]中系数误差最小为0.1%,最大误差为1.7%,平均误差为0.63%,而本文所提方法系数误差最小为0.001%,最大误差为0.4%,平均误差为0.13%,解调性能得到了极大提高.

6 蒙特卡洛分析

蒙特卡洛模拟法(Monte Carlo simulation)又称随机模拟法(random simulation)[22],随机抽样技术(random sampling)或统计实验方法(statistical testing).和传统数学方法相比,蒙特卡洛模拟法具有思想新颖、直观性强、简便易行的优点,便于解决难以用数学分析方法求解的动态系统复杂问题[23].当所要求解的问题是某种事件出现的概率,或者是某个随机变量的期望值时,可以通过某种“试验”的方法,得到这种事件出现的频率,或者这个随机变数的平均值,并用它们作为问题的解.这就是蒙特卡罗方法的基本思想.

由于信号中包含了相位噪声nε(t) 和加性电路噪声nc(t),所以每次仿真信号中的离群点分布都不尽相同,最后由椭圆拟合求出的椭圆系数也不同.通过蒙特卡洛分析将椭圆系数定量化评定可以得到该方法的规律,同时也验证了方法的稳定性.

假设经过n次仿真,A1的误差值为P1=[x1,x2,x3···xi···xn],误差概率密度函数为p1,同理A2,C1,B1,B2的误差值分别为P2,P3,P4,P5,误差概率密度函数分别为p2,p3,p4,p5.由于各系数的误差都会影响最后的解调,综合评定系数的误差为E=5/(P1+P2+P3+P4+P5),可以通过相关函数拟合得到综合评定系数误差概率密度函数e.本文所提方法在第5 节仿真参数下各系数的误差概率密度以及综合评定系数误差概率密度如图8所示.

从图8 可以看出,综合评定系数呈正态分布,误差范围在0.13%附近.试验样本数为20000,即通过20000 次试验验证,该方法的误差率约为0.13%.

图8 综合评定系数误差概率密度Fig.8.Error probability density of comprehensive evaluation coefficient.

7 结论

本文以光纤振动传感中基于3×3 耦合器的迈克尔逊干涉仪结构为例,分析了在实际应用中电路噪声和相位偏移对3×3 耦合器椭圆拟合解调的影响,提出了利用LOF 算法来处理噪声干扰的椭圆拟合问题,相比传统直接利用最小二乘拟合的方法,本文所提方法使信号解调的精度得到了较大提高.通过仿真解调结果对比,验证了所提方法的有效性.然后进一步利用蒙特卡洛方法,计算得到了该方法的误差率约为0.13%.证明了所提方法的稳定性.