光纤振动传感相位解调的局部异常因子优化椭圆拟合方法

1.本发明涉及光纤通信技术领域，尤其涉及一种光纤振动传感相位解调的局部异常因子优化椭圆拟合方法。


背景技术：

2.随着光纤通信技术的发展，光纤传感技术也成为国内外研究的热点
[1-7]
。光纤传感器因其稳定性好、可靠性高、对电磁波的不敏感性等特点，因而广泛应用于高温、腐蚀性或危险性的环境中
[8-9]
。信号解调是光纤传感器应用中的一个重要研究方向，而干涉法解调以其高分辨率的独特优势深受青睐，是最常用到的解调方法之一。其中，基于3
×
3耦合器的解调算法对激光源的要求较低，稳定性高
[10]
，在实践得到广泛应用。常用的于3
×
3耦合器解调方法有反正切法、微分交叉相乘法(differential cross multiplication,dcm)
[11]
和椭圆拟合估计法等
[12]
.微分交叉相乘法要求三路输出信号相差120
°
，但在实际应用中3
×
3耦合器无法达到完全对称，各路输出功率有所差别，加上插入损耗等因素，需要对信号进行归一化处理
[13]
，增加了算法难度。椭圆拟合估计法相比于不要求耦合器完全对称，增强了系统的适应性；无需对信号直流补偿，简化了算法。
[0003]
liu等利用最小二乘法对每两路输出信号构成的李萨如图形进行椭圆拟合反解出交直流项系数，进而求解相位变化量
[14]
。王伟利用双参数椭圆拟合算法实现了法布里-珀罗腔长的解调
[15]
。梅泽等聚焦弱信号解调，改进后的算法可以准确解调弱振动信号，提升了探测分辨率
[16]
。但以上研究只考虑了3
×
3耦合器两路信号不完全相差120
°
这一情况，未考虑实际解调中电路噪声和相位噪声的干扰，噪声干扰会直接影响椭圆拟合算法中的核心——李萨如图形，从而影响解调精度。高晓文为降低信号畸变提出的改进方案中考虑了噪声干扰这一实际因素
[17]
，但设定的干扰信号极其微弱，几乎相当于理想状态。综合考虑3
×
3耦合器两路信号不完全相差120
°
以及实际中噪声干扰的问题。


技术实现要素：

[0004]
本发明目的就是为了弥补已有技术的缺陷，提供一种光纤振动传感相位解调的局部异常因子优化椭圆拟合方法，相较于传统直接利用最小二乘拟合的方法能有效降低电路噪声和相位噪声对椭圆拟合解调的影响，提高解调精度。
[0005]
本发明是通过以下技术方案实现的：
[0006]
一种光纤振动传感相位解调的局部异常因子优化椭圆拟合方法，具体包括以下步骤：
[0007]
(1)采用椭圆拟合的方法对光纤耦合器进行解调，解得椭圆方程系数，再与信号系数对应，计算出信号系数；
[0008]
(2)通过仿真系统验证噪声对于椭圆拟合的影响；
[0009]
(3)利用局部异常因子算法对椭圆拟合进行优化；
[0010]
(4)通过数值仿真进行验证利用局部异常因子算法对椭圆拟合进行优化的有效
光纤耦合器的微分交叉相乘还原算法，解得相位变化：
[0025][0026]
步骤(2)所述的通过仿真系统验证噪声对于椭圆拟合的影响，具体如下：
[0027]
假设相位噪声n
ε
(t)和加性电路噪声nc(t)都用正态分布随机噪声模拟，两路信号表示为：
[0028][0029]
为了验证噪声对干涉信号的影响，在仿真中信号幅值的单位选为电压幅值，假设直流量a1、a2分别为和10v、12v，b1′
、b2′
的值分别为9v和6v；分别用幅值为0.5v和1.5v的随机噪声来模拟相位噪声n
ε
(t)和加性电路噪声nc(t)，在传感光纤上叠加扰动信号：
[0030][0031]
令仿真系统的采样频率为fs＝20000hz，采样数为n＝10000，绘出两路仿真信号图；通过两路仿真信号图两图的对比看出，噪声对两路信号均有明显干扰，而两路信号形成的李萨如图是椭圆拟合的基础，两路信号受噪声干扰，由两路信号形成的李萨如图也势必受到影响，绘出两路信号的李萨如图，从李萨如图看出，最小二乘法拟合的椭圆曲线椭圆与理想椭圆曲线相比倾角有明显偏差。
[0032]
步骤(3)所述的利用局部异常因子算法对椭圆拟合进行优化，具体如下：
[0033]
对于样本集合中的n个监测样本，任意的两个样本数据点之间的距离采用闵可夫斯基距离公式计算，如公式(9)所示：
[0034][0035]
当公式(9)中q为2时，即为欧式距离；
[0036]
如果数据集中至少有不包括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之间的距离；
[0037]
数据点p的第k距离领域nk(p)指在数据集中所有与p点之间的距离小于或等于dk(p)的数据点的集合，所以nk(p)中的元素个数大于或等于k；
[0038]
如果点o是远离点p的，那么点o到点p的第k可达距离就是o、p间的真实距离，否则可达距离用0的第k距离表示，即o到距离它最近的k个点的可达距离认为是相等的，并且都等于dk(o)，所以可达距离的公式定义是：
[0039]
reach_dk(p，o)＝max{d(p，o)，dk(o)}
ꢀꢀꢀ
(10)
[0040]
点p的局部可达密度公式表示为：
[0041][0042]
局部可达密度是指点p的第k邻域内的点到p的平均可达距离的倒数，局部可达密度越高则越有可能属于正常数据点，局部可达密度越低则越有可能是离群点；如果样本点与其近邻域的点不属于同一个聚簇，那么它的可达距离就是数值较大的d(p，o1)，所以局部
可达密度更小，但如果样本点与其近邻域的点属于同一个聚簇，那么可达距离就是值较小的dk(o)，所以局部可达密度会更大；
[0043]
点p的局部异常因子公式表示为：
[0044][0045]
点p的局部异常因子表示点p的领域点nk(p)的局部可达密度与点p的局部可达密度之比的平均数，如果局部异常因子小于1，则表示点p的局部可达密度高于它的领域点的局部可达密度，那么点p属于密集点；如果局部异常因子等于1，则表示点p的局部可达密度与它的领域点的局部可达密度差不多，那么点p与它的邻近点可能属于同一聚簇；如果局部异常因子大于1，则表示点p的局部可达密度低于它的领域点的局部可达密度，那么点p属于离群点；
[0046]
由于噪声影响无法直接利用最小二乘等数值估计方法来准确获取椭圆表达式，噪声会在标准拟合曲线内外侧形成无规则的离散点，利用局部异常因子算法判断这些离散点哪些为离群点，并且将这些离群点剔除；将离群点剔除后，使得李萨如图形趋向于无噪声干扰时所合成的椭圆，从而较准确地获取有效信号的李萨如图形椭圆表达式，再反解椭圆系数，从而消除电路噪声和相位偏移对信号解调的影响；将局部异常因子算法处理李萨如图形中离群点看作图像处理类方法，首先将两路信号生成的图像矩阵转化为灰度图像再转化为二值图，通过局部异常因子算法处理减少噪声的影响，最后利用最小二乘拟合求出椭圆系数并利用微分交叉相乘方法完成扰动信号的相位还原。
[0047]
步骤(4)所述的通过数值仿真进行验证利用局部异常因子算法对椭圆拟合进行优化的有效性，具体如下：
[0048]
利用数值仿真进行验证，按照步骤(2)所述的设置信号参数，a1、a2分别为10v、12v，b1′
、b2′
的值分别为9v和6v；分别用幅值为0.5v和1.5v的随机噪声来模拟相位噪声n
ε
(t)和加性电路噪声nc(t)，扰动信号设置为：(t)，扰动信号设置为：两路信号形成的李萨如图像进行灰度化和二值化处理，再经局部异常因子算法处理，通过局部异常因子算法处理前后对比看出，经过局部异常因子算法处理后的李萨如图像离散点得到有效去除，图像左下部离散点较少，椭圆轮廓明显收窄，从而有利于后续椭圆拟合；经局部异常因子算法处理的椭圆拟合解调法与直接利用最小二乘算法的解调法相比，前者所得的解调信号更接近原始信号。
[0049]
步骤(5)所述的利用蒙特卡洛方法计算利用局部异常因子算法对椭圆拟合进行优化后的误差率，并验证其稳定性，具体如下：
[0050]
假设经过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；绘制各系数的误差概率密度以及综合评定系数误差概率密度图，得出综合评定系数呈正态分布，误差率为0.13％。
[0051]
本发明的优点是：本发明以光纤振动传感中基于3
×
3耦合器的迈克尔逊干涉仪结构为例，分析了在实际应用中电路噪声和相位偏移对3
×
3耦合器椭圆拟合解调的影响，提
出了利用lof算法来处理噪声干扰的椭圆拟合问题，相比传统直接利用最小二乘拟合的方法，本发明所提方法使信号解调的精度得到了较大提高；通过仿真解调结果对比，验证了所提方法的有效性。然后进一步利用蒙特卡洛方法，计算得到了该方法的误差率约为0.13％.证明了所提方法的稳定性。
附图说明
[0052]
图1为基于3
×
3耦合器解调的迈克尔逊干涉仪结构。
[0053]
图2为两路仿真信号图，(2a)为无噪声干扰信号图；(2b)为有噪声干扰信号图。
[0054]
图3为两路信号的李萨如图形及拟合曲线图。
[0055]
图4为k＝5时，对象p，o1的可达距离和对象p，o2的可达距离的示意图。
[0056]
图5为系统流程图。
[0057]
图6为lof算法处理前后的李萨如图像，(6a)为处理前的李萨如图像；(6b)为处理后的李萨如图像。
[0058]
图7为解调信号图，(7a)为完整解调信号图，(7b)为解调信号局部图。
[0059]
图8为综合评定系数误差概率密度图。
具体实施方式
[0060]
1、系统光路
[0061]
光纤振动传感中基于3
×
3耦合器的迈克尔逊干涉仪结构如图1所示。它由环形器、3
×
3耦合器、压电陶瓷(piezoelectric transducer，pzt)和法拉第反射镜组成。光信号通过环形器1-2端口，首先经3
×
3耦合器的两臂分为2束光进入干涉仪，然后分别被法拉第反射镜反射之后，在3
×
3耦合器重新会合发生干涉，最终包含振动相位调制信息的光信号分别由3个光电探测器接收。引入法拉第反射镜的目的是消除偏振态对干涉信号的影响。
[0062]
2、椭圆拟合解调
[0063]
2.1基于椭圆拟合的3
×
3光纤耦合器的解调原理
[0064]
对于理想的3
×
3耦合器，假设条纹可见性相等，则耦合器三路输出信号中的两路经光电转换后为
[18]
：
[0065][0066]
其中，θ＝120
°
，i1，i2与耦合器性能有关。
[0067]
实际工程应用中，理想的3
×
3耦合器并不存在，3
×
3耦合器输出信号的条纹可见性不相等，信号之间的相位差θ也不稳定在120
°
，因此耦合器的输出为：
[0068][0069]
其中：a1≠a2，b1′
≠b2′
，θ
′
＝120
°‑
δ，δ是由耦合器非对称性造成的微小相位偏差。由于光纤传感器受到外界环境如温度改变、大地蠕变等的影响，a1，a2，b1′
，b2′
，δ会随时间发生缓慢的变化。但由于在本系统中，每次解调计算都是针对1s所采集的数据，而温度和大地条件的变化则要缓慢得多，可以认为是同一次解调运算，a1，a2，b1′
，b2′
，δ均是常数。式(2)可展开成：
[0070][0071]
检测信号系数a1，a2，b1，b2，c2的求解建立在椭圆拟合的基础上，椭圆合成遵循李萨如图的性质和规律。由式(2)可知，耦合器输出信号u1和u2对同一相位函数进行调制，其变化形式与简谐振动一致，变化频率相同，相位差θ
′
恒定，符合李萨如图的形成条件。此外，由于两路信号的相位差接近120
°
，如果在互相垂直的两个坐标轴上同时画出u1，u2，其轨迹能形成一个椭圆，椭圆方程可表示为：
[0072][0073]
式中，g＝[a，b，c，d，e，f]
t
为椭圆的系数矢量.对于式(4)所示的二次曲线，当且仅当b
2-4ac＜0时，曲线为椭圆.为方便计算，选取b
2-4ac＝0，采用一种适当的椭圆拟合方法解得椭圆方程系数，再与信号系数对应，计算出信号系数，传统主要采用最小二乘拟合.信号系数计算公式为：
[0074][0075]
将式(5)解得的信号系数代入信号方程(3)，可以获得与信号.采用3
×
3耦合器的微分交叉相乘还原算法，解得相位变化：
[0076][0077]
3.2噪声对于椭圆拟合的影响
[0078]
在实际解调过程中，输出信号经过光电转换为电流信号，再通过电路转换成电压信号，这一过程不可避免的会引入电路噪声，而在解调过程中还存在相位偏移量。这些噪声都会对解调产生不利影响。假设相位噪声n
ε
(t)和加性电路噪声nc(t)都用正态分布随机噪声模拟，两路信号可表示为：
[0079][0080]
为了验证噪声对干涉信号的影响，在仿真中信号幅值的单位选为电压幅值，假设直流量a1、a2分别为和10v、12v，b1′
、b2′
的值分别为9v和6v；分别用幅值为0.5v和1.5v的随机噪声来模拟相位噪声n
ε
(t)和加性电路噪声nc(t)，在传感光纤上叠加扰动信号：
[0081][0082]
令仿真系统的采样频率为fs＝20000hz，采样数为n＝10000，则两路仿真信号如图2所示。
[0083]
通过图(2a)、(2b)两图的对比可以看出，噪声对两路信号均有明显干扰。而两路信
号形成的李萨如图是椭圆拟合的基础，两路信号受噪声干扰，由两路信号形成的李萨如图也势必受到影响。两路信号的李萨如图形如图3所示。
[0084]
从图3可以看出，噪声产生的离散点在理想椭圆曲线周围不规则分布，椭圆内外轮廓的离散点使得椭圆边界变“厚”，导致很难直接利用最小二乘等数值估计方法准确获取椭圆表达式，对比可以看出，最小二乘法拟合的椭圆曲线椭圆与理想椭圆曲线相比倾角有明显偏差。
[0085]
3、利用局部异常因子算法(lof)对椭圆拟合进行优化
[0086]
3.1lof算法原理
[0087]
局部离群因子(local outlier factor)检测算法是一种无监督的异常检测方法，是基于密度的局部离群点检测算法
[19]
。算法内核是比较目标点与相邻点的局部密度偏差，通过计算局部密度偏差的比值来判断目标点是否为离群点。由于算法在计算目标点时同时计算了相邻点，目标点的局部密度偏差比值包含了邻近点的信息，算法的精度较高。lof算法的相关定义和计算方法如下。
[0088]
3.1.1样本点之间的距离
[0089]
对于样本集合中的n个监测样本，任意的两个样本数据点之间的距离可以采用欧氏(eucildean)距离、汉明(hamming)距离、马氏(mahalanobis)距离、切比雪夫(chebyshev)距离、闵可夫斯基(minkowski)距离和球面距离等。其中，闵可夫斯基距离公式如公式(9)所示：
[0090][0091]
当公式(9)中q为2时，即为欧式距离。
[0092]
3.1.2第k距离
[0093]
如果数据集中至少有不包括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之间的距离。
[0094]
3.1.3第k距离领域
[0095]
数据点p的第k距离领域nk(p)就是指在数据集中所有与p点之间的距离小于或等于dk(p)的数据点的集合。所以nk(p)中的元素个数大于或等于k。
[0096]
3.1.4可达距离
[0097]
如果点o是远离点p的，那么点o到点p的第k可达距离就是o、p间的真实距离，否则可达距离用o的第k距离表示。也就是说，o到距离它最近的k个点的可达距离可以认为是相等的，并且都等于dk(o)。所以可达距离的公式定义是：
[0098]
reach_dk(p，o)＝max{d(p，o)，dk(o)}
ꢀꢀ
(10)
[0099]
3.1.5局部可达密度
[0100]
点p的局部可达密度公式可以表示为：
[0101][0102]
局部可达密度是指点p的第k邻域内的点到p的平均可达距离的倒数。局部可达密
度越高则越有可能属于正常数据点，局部可达密度越低则越有可能是离群点。如果样本点与其近邻域的点不属于同一个聚簇，那么它的可达距离就是数值较大的d(p，o1)，所以局部可达密度更小，但如果样本点与其近邻域的点属于同一个聚簇，那么可达距离就是值较小的dk(o)，所以局部可达密度会更大。
[0103]
3.1.6局部异常因子
[0104]
点p的局部异常因子公式可以表示为：
[0105][0106]
点p的局部异常因子表示点p的领域点nk(p)的局部可达密度与点p的局部可达密度之比的平均数。如果局部异常因子小于1，则表示点p的局部可达密度高于它的领域点的局部可达密度，那么点p属于密集点；如果局部异常因子等于1，则表示点p的局部可达密度与它的领域点的局部可达密度差不多，那么点p与它的邻近点可能属于同一聚簇；如果局部异常因子大于1，则表示点p的局部可达密度低于它的领域点的局部可达密度，那么点p属于离群点。
[0107]
3.2局部异常因子算法(lof)对椭圆拟合进行优化
[0108]
上面讨论了噪声对于椭圆拟合的影响，指出由于噪声影响无法直接利用最小二乘等数值估计方法来准确获取椭圆表达式。由图3可以看出噪声会在标准拟合曲线内外侧形成无规则的离散点，利用lof算法可以判断这些离散点哪些为离群点，并且将这些离群点剔除。将离群点剔除后，可以使得李萨如图形趋向于无噪声干扰时所合成的椭圆。从而可以较准确地获取有效信号的李萨如图形椭圆表达式，再反解椭圆系数，从而消除电路噪声和相位偏移对信号解调的影响。lof算法处理李萨如图形中离群点可以看作图像处理类方法，所以首先要将两路信号生成的图像矩阵转化为灰度图像再转化为二值图，通过lof算法处理减少噪声的影响。最后利用最小二乘拟合求出椭圆系数并利用微分交叉相乘方法(dcm)完成扰动信号的相位还原。系统的总体设计流程图如图5所示.
[0109]
4、仿真验证
[0110]
利用数值仿真进行验证，按照上面部分设置信号参数，a1、a2分别为和10v、12v，b1′
、b2′
的值分别为9v和6v；分别用幅值为0.5v和1.5v的随机噪声来模拟相位噪声n
ε
(t)和加性电路噪声nc(t)，扰动信号设置为：(t)，扰动信号设置为：两路信号形成的李萨如图像进行灰度化和二值化处理，再经lof算法处理，处理结果如图(6b)所示，处理前的李萨如图像如图(6a)所示。
[0111]
通过lof算法处理前后对比可以看出，经过lof算法处理后的李萨如图像离散点得到有效去除，图像左下部离散点较少，椭圆轮廓明显“收窄”，从而有利于后续椭圆拟合。经过lof算法处理并采用微分交叉相乘法得到的解调后的完整信号如图(7a)所示。为清晰起见，图(7b)给出了图(7a)中红色虚线框部分的放大图。
[0112]
由图7可见，经lof处理的椭圆拟合解调法与直接利用最小二乘算法的解调法相比，前者所得的解调信号更接近原始信号。故该方法可以较好地还原扰动信号，降低解调信号的误差。对比其他改进最小二乘解调法，可以看到文献
[20]
中系数误差最小为0.1％，最大误差为1.7％，平均误差为0.63％，而本文所提方法系数误差最小为0.001％，最大误差为
0.4％，平均误差为0.13％，解调性能得到了极大提高。
[0113]
5蒙特卡洛分析
[0114]
蒙特卡洛模拟法(monte carlo simulation)又称随机模拟法(random simulation)
[21]
，随机抽样技术(random sampling)或统计实验方法(statistical testing)。和传统数学方法相比，蒙特卡洛模拟法具有思想新颖、直观性强、简便易行的优点，解决难以用数学分析方法求解的动态系统复杂问题
[22]
。当所要求解的问题是某种事件出现的概率，或者是某个随机变量的期望值时，可以通过某种“试验”的方法，得到这种事件出现的频率，或者这个随机变数的平均值，并用它们作为问题的解。这就是蒙特卡罗方法的基本思想。
[0115]
由于信号中包含了相位噪声n
ε
(t)和加性电路噪声nc(t)，所以每次仿真信号中的离散点分布都不尽相同，所以最后由椭圆拟合求出的椭圆系数也不同。通过蒙特卡洛分析将椭圆系数定量化评定可以得到该方法的规律，同时也验证了方法的稳定性。
[0116]
假设经过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所示。
[0117]
从图8可以看出综合评定系数呈正态分布，误差范围在0.13％附近。试验样本数为20000，即通过20000次试验验证，该方法的误差率约为0.13％。
[0118]
本发明以光纤振动传感中基于3
×
3耦合器的迈克尔逊干涉仪结构为例，分析了在实际应用中电路噪声和相位偏移对3
×
3耦合器椭圆拟合解调的影响，提出了利用lof算法来处理噪声干扰的椭圆拟合问题，相比传统直接利用最小二乘拟合的方法，本发明所提方法使信号解调的精度得到了较大提高。通过仿真解调结果对比，验证了所提方法的有效性。然后进一步利用蒙特卡洛方法，计算得到了该方法的误差率约为0.13％.证明了所提方法的稳定性。