基于连续小波变换的多普勒信号瞬时速度提取方法

摘要

本发明公开了一种基于连续小波变换的多普勒信号瞬时速度提取方法，属于信号处理技术领域。本发明采用连续小波变换得到输出信号的瞬时频率，然后对各个时刻的瞬时频率进行多项式拟合以减小误差，最后根据拟合后各个时刻点的瞬时频率解算出各个时刻的瞬时速度，适用于没有使用延迟光路和参考光路的干涉仪。本发明可以更好地做到同时保留时域和频域信息，得到的信号平滑性更好，且不受边界效应限制，应用范围更广，精度更高，可以较准确地提取到各个时刻的速度值。

权利要求

1.一种基于连续小波变换的多普勒信号瞬时速度提取方法，其特征在于，该方法适用于没有使用延迟光路和参考光路的干涉仪，采样点不超过10000个，采样率在最大频率的5-

7倍，通过分析输出信号各个时刻的瞬时频率得到各个时刻的多普勒频移量，进而求出各个时刻的瞬时速度；所述方法的实现步骤如下：步骤1)设第k个采样点中的第i次迭代尺度值为ai(k)，ai(k)＝ai(t0+kTs)；其中，t0为初始时刻，Ts为采样周期；给定初始时刻的初始尺度值为a0(t0)；

步骤2)设当前迭代的尺度值为ai(k)，利用小波变换得到该尺度值对应的小波系数WTa(k)，设 为WTa(k)的相位，则得到第k个采样点的第i+1次迭代的尺度值ai+1(k)为：其中，Db表示离散微分算法，ω0为基小波的中心频率；

步骤3)根据设定的精度ε判断迭代的尺度值是否满足下式：

当满足该式时，停止对第k个采样点的尺度值的迭代，标记所得到的第k个采样点的收敛尺度值ai+1(k)为ac(k)，进入步骤4；否则，将ai+1(k)作为新的迭代值ai(k)，进行步骤2执行；

步骤4)将第k个采样点的收敛尺度值ac(k)作为第k+1个采样点的迭代初始值a0(k+1)；

步骤5)k自增1，然后重复步骤2～4，直至把所有的采样点都遍历到，获得各采样点的收敛尺度值；

步骤6)根据下式求取待估信号在各采样点的瞬时频率

其中， 为基小波g(t)在初始时刻的瞬时频率；

步骤7)对各个时刻的瞬时频率进行多项式拟合，t时刻的瞬时频率拟合为多项式表示如下：其中，p为多项式拟合中的最高次幂；采用最小二乘估计准则估计系数kj；确定系数kj后根据拟合的多项式重新估计各个时刻的瞬时频率f(t)；

步骤8)确定t时刻的瞬时速度 其中，λ为探测光的波长，θ为探测光照射到运动物体的入射角度。

2.根据权利要求1所述的一种基于连续小波变换的多普勒信号瞬时速度提取方法，其特征在于，步骤3中所述的精度ε设置范围为[0.05,0.5]。

3.根据权利要求1所述的一种基于连续小波变换的多普勒信号瞬时速度提取方法，其特征在于，所述的步骤3中，设置最大迭代次数为10，控制迭代次数。

4.根据权利要求1所述的一种基于连续小波变换的多普勒信号瞬时速度提取方法，其特征在于，所述的步骤7中，设置多项式拟合中的最高次幂p为3。

5.根据权利要求1所述的一种基于连续小波变换的多普勒信号瞬时速度提取方法，其特征在于，所述的多普勒信号瞬时速度提取方法用于马赫-泽德干涉仪。

说明书

基于连续小波变换的多普勒信号瞬时速度提取方法

技术领域

[0001] 本发明涉及一种基于连续小波变换的多普勒信号瞬时速度提取方法，尤其是对马赫-泽德干涉仪输出的多普勒信号进行瞬时速度提取，属于信号处理技术领域。

背景技术

[0002] 测速技术在日常生活中的应用非常广泛，如汽车行驶速度的测量、电机转速的测量。激光多普勒测速是一种常见的测速方法，其中双光束多普勒测速可达到较高精度，应用较为广泛。但是，这种方法难以进行离焦测量，在测量固体运动的速度时会比较困难，通常只能用来测量流体的运动以及固体的微小振动。激光干涉测速的出现很好地解决了测量固体速度的问题，其原理是将激光器发出的光分成两束，即信号光与参考光，信号光照射到运动物体上后发生漫反射，反射的光与参考光产生干涉现象，通过分析干涉仪输出的多普勒信号得到物体运动的速度。然而，如何准确地从多普勒信号中提取瞬时速度一直是激光干涉测速领域的重点和难点。

[0003] 目前，虽然用来测量速度的激光干涉仪种类很多，但是其信号处理方法主要有三种，第一种方法是将输出信号分成两路，然后通过示波器或数据采集卡测量两路信号的相位差来实现速度的测量，这种方法最大的弊端就是增加了一路输出信号，需要额外的通道采集数据，而激光干涉测速法的输出信号往往频率较高，因此增加一路通道带来的成本较大；第二种方法是通过测量干涉条纹的变化来实现速度的测量，这种方法虽然不需要再将输出信号分成两路，但是直接观测条纹的变化并不容易，因此难以达到较高的精度，虽然采用CCD相机可在一定程度上解决这一问题，但成本也会增加很多；最后一种方法就是直接通过信号处理算法提取信号的瞬时频率，然后根据瞬时频率提取瞬时速度，这种方法虽然不会带来成本的增加，但目前并没有很成熟的算法可以准确地从多普勒信号中提取出速度信号，虽然可以使用FFT(快速傅里叶变换)算法来处理多普勒信号，但是由于多普勒信号的频率随时间变化较快，采用傅里叶变换的方式难以同时保留信号的频率和时间信息，因此测量精度会受到限制，而且计算速度仍有待提高。

发明内容

[0004] 本发明针对上述现有存在的问题，提供了一种基于连续小波变换的多普勒信号瞬时速度提取方法，采用连续小波变换得到输出信号的瞬时频率，然后对各个时刻的瞬时频率进行多项式拟合以减小误差，最后根据拟合后各个时刻点的瞬时频率解算出各个时刻的瞬时速度。小波变换与傅里叶变换相比，可以自动调节自身的尺度来对信号作分解，因此可以更准确地获得信号在频率和时间上的信息，而且对噪声有较好的抑制作用，非常适合多普勒信号的处理。

[0005] 本发明的一种基于连续小波变换的多普勒信号瞬时速度提取方法，适用于没有使用延迟光路和参考光路的干涉仪，最适合测量较短时间内速度变化的过程，采样点通常不超过10000个。为了达到更好的效果，通常设置采样率在最大频率的5-7倍。

[0006] 多普勒信号瞬时速度提取方法的实现步骤如下：

[0007] 步骤1)设第k个采样点的第i次迭代尺度值为ai(k)，ai(k)＝ai(t0+kTs)；其中，t0为初始时刻，Ts为采样周期；首先给定初始时刻t0的初始尺度值a0(t0)。

[0008] 步骤2)设当前迭代的尺度值为ai(k)，利用小波变换得到尺度值ai(k)在b＝kTs时刻对应的小波系数WTa(k)＝WT(kTs,ai(k))，设 为WTa(k)的相位，则得到第k个采样点的第i+1次迭代的尺度值ai+1(k)为：

[0009]

[0010] 其中，Db表示离散微分算法，ω0为采用的基小波的中心频率。

[0011] 步骤3)根据指定的精度ε判断当前迭代的尺度值ai+1(k)是否满足下式：

[0012]

[0013] 若满足该式，停止对第k个采样点的尺度值的迭代，标记所得到的第k个采样点的收敛尺度值ai+1(k)为ac(k)，进入步骤4执行；否则，进行步骤2执行。迭代精度ε的设置范围为[0.05,0.5]，因为过高的迭代精度会使计算速度变慢。考虑到个别尺度值有可能不收敛，步骤3还可设置最大迭代次数为10，来控制迭代过程。

[0014] 步骤4)将第k个采样点的尺度值ac(k)作为第k+1个采样点的迭代初始值a0(k+1)；

[0015] 步骤5)k自增1，重复步骤2～4，直至把所有的采样点都遍历到，并获得各采样点的收敛尺度值。

[0016] 步骤6)根据下式求取待估信号在各采样点的瞬时频率

[0017]

[0018] 为待估信号在第k个采样点的相位， 为基小波在初始时刻的相位。

[0019] 步骤7)对各个时刻的瞬时频率进行多项式拟合，t时刻的瞬时频率拟合为多项式表示如下：

[0020]

[0021] 其中，p为多项式拟合中的最高次幂，通常根据运动的情况选择合适的p值，通常情况下取p＝3即可达到较好的效果；采用最小二乘估计准则估计系数kj；确定系数kj后根据拟合的多项式重新估计各个时刻的瞬时频率f(t)。

[0022] 步骤8)根据瞬时速度与瞬时频率的关系 将波长λ与探测光照射到运动物体的入射角度θ代入，可求得瞬时速度u(t)。

[0023] 本发明所具有的有益效果是：

[0024] (1)本发明可直接通过信号处理算法从多普勒信号中提取到瞬时速度，既不需要将两路信号分开处理，也不需要额外的条纹监视系统，因此可大大节约成本。

[0025] (2)与基于傅里叶变换的各种多普勒信号处理算法相比，本发明采用了连续小波变换处理多普勒信号，因此可以更准确地保留多普勒信号在频率和时间上的局部信息，计算更简单，处理速度更快。

[0026] (3)本发明采用了多项式拟合的方式来减小噪声信号的影响，减小了噪声信号干扰对测量结果的影响，可以达到更高的精度。

附图说明

[0027] 图1为本发明的多普勒信号瞬时速度提取方法的流程图；

[0028] 图2为基于马赫-泽德干涉测速的装置原理图；

[0029] 图3为马赫-泽德干涉仪输出的多普勒信号示意图；

[0030] 图4为采用本发明方法使用连续小波变换估计出的一个频率曲线示意图；

[0031] 图5为采用本发明方法拟合得到的一个频率曲线示意图；

[0032] 图6为采用本发明方法得到的速度曲线与理论速度曲线的对比示意图。

具体实施方式

[0033] 下面将结合附图和实例对本发明作进一步的详细说明。

[0034] 本发明是一种基于连续小波变换的多普勒信号瞬时速度提取方法，整体流程如图1所示。首先通过连续小波变换解算出各个时刻点的瞬时频率，然后采用多项式拟合的方法对各个时刻的瞬时频率进行修正，最后根据瞬时速度与瞬时频率之间的关系得到瞬时速度。

[0035] 首先，说明本发明的多普勒信号瞬时速度提取方法的实现原理，图2为一种马赫-泽德干涉测速装置，激光器(Laser)发出的光经分束镜M1分成探测光和参照光两束，探测光经反射镜M2和分束镜M4后照射到运动物体上发生漫反射，反射光经分束镜M4反射后与参照光发生叠加产生干涉现象，其中波片主要用来调整探测光的强度，不考虑探测不到的高频分量，该装置输出信号的光强I为：

[0036]

[0037] 其中，Is、Il分别为信号光和参照光的光强，ωs、ωl分别为信号光和参照光的角频率， 分别为信号光和参照光的相位。

[0038] 根据光电探测器与光强成正比这一特性可知，光电探测器检测到的多普勒信号是一个带直流分量的变频余弦信号，采用高通滤波器可以将第一项的直流分量滤掉，根据多普勒测速的基本原理可知，多普勒频移量

[0039] 因此，通过分析输出信号各个时刻的瞬时频率即可得到各个时刻的多普勒频移量大小，进而根据公式 即可求出各个时刻的瞬时速度，为了测量方便以及减小误差，通常可使探测光垂直照射到运动物体上，此时入射角度θ＝90°，λ为探测光的波长。

[0040] 图3为马赫-泽德干涉仪输出的多普勒信号，采样频率为250MHz，图3中横坐标为时间，单位为s，纵坐标为信号振幅，单位为V。通常多普勒信号可看作渐进信号，即相位变化的速度远大于其幅度变化的速度，可设多普勒信号为s(t)。同时基小波g(t)也为渐进信号。s(t)和g(t)如下：

[0041]

[0042]

[0043] 其中，As(t)、Ag(t)分别为待估信号s(t)和基小波g(t)的幅值， 分别为待估信号s(t)和基小波g(t)的相位。

[0044] 则进行连续小波变换后得到的小波系数WTs(b,a)为

[0045]

[0046] 其中，g(a,b)(t)为使用的小波函数(基小波)， Φ(b,a)(t)为g(a,b)(t)的相位函数， a为尺度参数，b为时间偏移量。

[0047] 当Φ(b,a)(t)趋向于稳定的时候，Φ'(b,a)(ts(b,a))＝0将得到满足，即[0048]

[0049] 其中，ts(b,a)被定义为相位稳定点，上式中将ts(b,a)简写为ts。如果g(t)满足在g(0)处取得极大值，同时ts(b,a)无限接近于b的时候或者等于b的时候，则WTs(b,a)将在b时刻附近得到一个局部最大值，此时，将ts(b,a)＝b代入 则可获得信号s(t)在b时刻的瞬时频率值 为：

[0050]

[0051] 其中，满足ts(b,a)＝b的点(a,b)被称之为小波脊。ar(b)是b时刻的尺度参数值。因此，要求目标信号的瞬时频率，必须要求解小波脊(a,b)，即ar(b)。首先必须在一段连续的时间范围内搜索使得WTs(b,a)取得极大值对应的b；其次，要根据信号的整体变化趋势，在众多小波中选择合适的基小波，满足小波的在g(0)处取得极大值且还要反复调节基小波的频率参数以获得最大的性能。

[0052] 本发明方法采用离散序列的迭代算法求解信号的瞬时频率。假定Ts＝1/fs为采样周期，其中fs表示采样频率。那么可以得到信号sk＝s(kTs)的小波变换WTa(k)＝WTs(kTs,a)。此处k表示第k个采样点，在第k个采样周期采集，下标a表示尺度参数，WTs(kTs,a)表示b＝kTs时得到的小波系数。再令Ψa(k)为WTa(k)的相位，Db表示对kTs＝b时间点的差分算子，则可得到差分方程：

[0053]

[0054] 其中，ω0为采用基小波的中心频率(角频率)。

[0055] 采用本发明的多普勒信号瞬时速度提取方法，从马赫-泽德干涉仪的输出信号中提取速度信号包括以下几个步骤：

[0056] 步骤1)首先给定初始时刻t0的初始尺度值a0(t0)。

[0057] 步骤2)已知ai(k)是第k个采样点的第i次迭代尺度值。首先，利用小波变换，得到该尺度值对应的小波系数。然后，根据式 得该时刻点的第i+1次迭代的尺度值如下：

[0058]

[0059] Db表示离散微分算法，即 可通过上面的几个公式计算得到然后得到 最后可以求得ai+1(k)。

[0060] 步骤3)取迭代精度ε为0.05，当满足指定的精度ε时，即满足下面公式时，则解得第k个采样点的收敛尺度值ac(k)，ac(k)＝ac(t0+kTs)＝ai+1(k)，同时停止该时刻点的迭代，转到步骤4执行。否则，将ai+1(k)最为新的迭代值ai(k)，转步骤2执行。

[0061]

[0062] 所述的精度ε设置范围为[0.05,0.5]。考虑到个别尺度值有可能不收敛，还可设置最大迭代次数，可设置为10，来控制迭代次数。

[0063] 步骤4)令第k个采样点收敛的尺度值ac(k)为第k+1个采样点的迭代初始值，即[0064] a0(k+1)＝ac(k)

[0065] a0(k+1)为第k+1个采样点的初始迭代值。

[0066] 步骤5)更新k值，使k自增1，重复步骤2-4，直至把所有的采样点都遍历到，并获得各采样点对应的收敛尺度值。

[0067] 步骤6)由 可得到各个时刻点的瞬时频率曲线。

[0068] 在具体计算时，第k个采样点的瞬时频率 ac(k)为第k个采样点的收敛尺度， 为待估信号在第k个采样点的相位， 为基小波在初始时刻的相位，为基小波g(t)在初始时刻的瞬时频率。

[0069] 本发明实施例得到的一个频率曲线，如图4所示，该图是指某个时刻点的瞬时频率曲线，该曲线是对待估信号进行变换后得出来的，横轴表示时间，纵轴表示瞬时频率，即t-f图。

[0070] 步骤7)根据 对各个时刻的瞬时频率进行多项式拟合，取p＝3，采用最小二乘估计准则确定系数kj，然后根据上式得到各个时刻瞬时频率的估计值f(t)。

[0071] 在本步骤中，将步骤6中得到的每个 进行多项式拟合，第k个采样点对应的时刻t为t0+kTs。对第k个采样点的 在确定系数kj后根据拟合的多项式重新估计该时刻t的瞬时频率f(t)。本发明实施例得到的一个频率拟合曲线，如图5所示。

[0072] 从图3中可以看出采集到的信号具有较多毛刺，采用连续小波变换提取到各个点的瞬时频率也会带有毛刺，如图4所示。从图4可以看出，本发明方法的抗噪能力较强，并没有因为噪声而产生较大的毛刺，和短时傅里叶变换法相比可以更好地保留时域和频域的信息，由于速度理论上不能发生突变，因此图4所示的结果转换为速度值后并不符合实际情况，只代表一个大致的趋势。本发明方法采用了一个多项式拟合，将带有毛刺的曲线拟合成较平滑的曲线，如图5所示，更符合实际运动的效果，因为速度不能突变。图5为经过多项式拟合去噪后得到的结果，经过多项式拟合去噪后，可以看到毛刺基本被去除，可以得到较平滑的速度-时间曲线。本发明方法的优点在于可以更好地做到同时保留时域和频域信息(反例是短时傅里叶变换)，得到的信号平滑性更好，且不受边界效应限制(反例是希尔伯特变换)，应用范围更广，精度更高。

[0073] 本发明通过多项式拟合，获得的瞬时频率的估计值更加符合实际，与实际误差也比较小。

[0074] 步骤8)根据公式 可得到目标在各个时刻点的速度。

[0075] 本发明方法适用于没有使用延迟光路和参考光路的干涉仪，这种干涉仪发生干涉的两束光只有一束带有多普勒频移，输出信号的频率和速度成正比，根据公式可获得各个时刻的瞬时速度，f(t)为各个时刻的多普勒频移量。对于使用了延迟光路和参考光路的干涉仪，发生干涉的两束光都带有多普勒频移，只是频移的时刻和大小都不相同，一束光是延迟臂-运动物体-参考臂，另一束光是参考臂-物体-延迟臂，其余光的相干长度较大不会发生干涉，因此 这一公式将不适用于这种装置，对于这种装置的分析方法是通过光程差与相位进行分析。

[0076] 本发明实施例得到的一个目标速度曲线与理论速度曲线对比，如图6所示，从中可以看出，本发明方法拟合的曲线和理论曲线相差非常小，可以较准确地提取到各个时刻的速度值。