基于激光多普勒测振技术的动态分析方法研究

《Optics & Laser Technology》：Study of dynamic analysis method based on laser Doppler Vibrometry technology

　　向鹏|穆谷|李苏军|任启峰|赵如进
中国科学院空间光电精密测量科学与技术重点实验室，成都610209，中国

**摘要**
由于激光多普勒测量技术具有非接触性、响应速度快、灵敏度高和分辨率高等优异特性，其在现代高精度测量中受到了越来越多的关注和应用。本文研究了基于激光多普勒振动测量（LDVib）技术的动态数据分析方法，揭示了多普勒信号与目标运动参数之间的物理关系，并统一了平衡光电探测器输出信号的数学表示。一方面，针对采样率远高于振动频率的高采样率场景，提出了一种采样校正迭代DFT（SCI-DFT）方法。该方法通过“预测+补偿”方法确保了实时处理和高精度测量能力，有效解决了分层频率测量方案中多普勒速度序列采样不均匀的问题。另一方面，针对采样率接近振动频率的低采样率场景，将信号转换为贝塞尔函数表示的谐波展开形式，并提出了一种基于自相关的基本频率提取方法。该方法有效抑制了噪声干扰，实现了实时分析。实验结果表明，在这两种条件下，分析方法都能快速准确地获得高精度测量结果。本研究为激光多普勒测量技术的工程应用提供了理论基础和技术参考，对相关激光测量技术的发展做出了重要贡献。

**1. 引言**
对于运动物体或在动态条件下的物体，振动频率往往是最关键和直观的特征参数。传统的振动测量技术主要依赖于基于接触的方法，即将加速度计安装在目标表面，将电荷输出信号转换为加速度、速度或位移测量值。然而，诸如测量微振动、微小速度变化或无法接触的表面等操作限制，常常使得传感器安装变得困难。现代领域，包括机械制造、医疗诊断和航空航天，越来越需要高精度、高效且无损的部件检测。激光多普勒振动测量（LDVib）技术因其非接触性、简单的光学配置和成本效益而受到广泛关注——特别是在高风险或高压环境中进行精密测量时具有显著优势[1] [2] [3] [4] [5]。因此，开发具有高精度和实时能力的新非接触式振动测量技术已成为工程科学的重要研究方向[6] [7] [8]。

传统激光多普勒速度测量（LDV）技术的原理基于多普勒效应：当激光束照射到振动物体的表面时，由于目标的运动，反射激光的频率会发生改变（多普勒频移）。通过精确解调这种频率变化，可以直接获得物体的运动速度，然后通过进一步积分或微分可以得到位移和加速度。需要注意的是，LDV测量的是瞬时速度，而非平均速度。基于LDV，如果目标的运动具有周期性，可以通过分析一系列瞬时速度或相位序列的变化规律来获取振动频率，这就是LDVib技术[9] [10] [11] [12]。LDVib通常采用两种分析方法：
(1) 将整个测量信号波形记录下来并分割成短段，通过对每个段进行FFT处理得到一系列瞬时速度，然后通过分析这些速度序列提取振动信息。
(2) 应用相位生成载波（PGC）技术来解调中频（IF）载波信号的相位，通过对得到的相位序列进行二次分析来确定振动信息。

由于涉及二次解析操作，传统方法通常采用“先测量后分析”的数据处理方式。为了满足现代精密测量系统对实时处理和严格通道同步的要求，通常使用高速FPGA来进行计算任务并保证稳定输出[13]。然而，系统复杂性的增加不可避免地提高了成本并带来了限制。本文提出了一种简化的数据处理架构和实时动态数据分析方法，该方法采用“预测+补偿”策略以确保其实时性能和高精度测量能力。此外，还提出了一种基于自相关的分析方法，适用于采样率较低的场景，该方法可以绕过相位解包步骤，在低噪声条件下直接提取振动频率。这意味着仅通过主机计算机软件就能实现准确的振动分析，为测量目标的振动状态提供了一种低成本且简化的解决方案。

**2. 激光多普勒振动测量系统架构**
当被照物体的反射表面相对于观察者移动时，光的多普勒效应会导致干涉谱的频率偏移。通过将测量光束和参考光束同时投射到光电探测器的光敏表面上，会形成一个干涉图案[14] [15]。如图1所示，激光被光纤分光器分为参考光束和信号光束。信号光束通过光纤环行器传输，并通过光学天线投射到振动目标上，被调制后发生后向散射。同时，参考光束通过声光调制器（AOM）发生频率偏移。最终，参考光束和后向散射信号光束结合在一起，在光电探测器表面形成干涉图案[16]。

光电探测器对这两束光的干涉场做出响应，生成一个光电流信号，该信号不仅取决于入射光的强度，还取决于其频率和相位。通过检测干涉条纹强度的变化，可以确定目标的运动速度。为了减少光的回波反馈效应和环境噪声对干涉信号的影响，提高信噪比和光能利用率，实验测量系统中使用了光纤环形谐振器和平衡光电探测器。在光学混频器中与参考光束混合后，组合光束被送入平衡光电探测器。由DSP芯片和集成AD模块处理的IF电信号被放大，以提取目标的多普勒频移信息。这意味着通过两个高频激光的异频相干叠加，可以获得由多普勒频移调制的IF包络信号[17]。

入射参考光和信号光的电场分别为：
$$
E_1(t) = E_{10}\cos(2\pi f_1 t) + \phi_1 \\
E_2(t) = E_{20}\cos(2\pi f_2 t) + \phi_2
$$
当两者叠加时，得到的光电流为：
$$
I(t) = \alpha \cdot (E_1 + E_2) = I_1 + I_2 + \alpha \cdot E_{10}E_{20}\cos(2\pi f_1 t) + \phi_1 + \phi_2 + \alpha \cdot E_{10}E_{20}\cos(2\pi f_1 t) - \phi_2
$$
由于光电探测器的响应速率有限，它们只能在一定时间内检测高频光信号的平均值。在方程(3)中，$I_1 + I_2$是直流项；$f_1 + f_2$是数十亿赫兹的高频项，探测器无法响应，因此只能获得其平均值零。因此，探测器只能对信号光和参考光之间的差频$f_D = f_1 - f_2$产生有效响应。经过直流过滤后，测量端得到的信号电流仅包含频率偏移信息。通过分析干涉信号，可以获得目标的瞬时径向速度$V(t)$。根据多普勒频移原理，径向速度$f_D$与频率偏移$f$之间的关系为：
$$
f_D = f_{AOM} + 2V(t)\frac{\lambda}{c}
$$
其中，$f_{AOM}$表示由AOM引起的异频频移，用于确定被测目标的运动方向。在LDVib应用中，为了防止光电流信号产生零交叉反转效应，$f_{AOM}$应至少是目标振动频率的两倍。

平衡光电探测器由一对响应特性几乎相同的光电二极管组成，并与跨阻放大器（TIA）耦合。光电二极管将输入的光信号转换为电流信号。TIA阶段将两个光电二极管之间的差分输出电流转换为比例电压信号，然后该电压被放大并传送到后续信号处理单元。信号光和参考光在通过保偏光纤耦合器后发生相位反转，耦合器的输出为：
$$
i_1(t) = \alpha \cdot (E_1 + E_2)^2 = I_1 + I_2 + I_1I_2\cos(2\pi f_D t) + \Delta\phi + I_{noise} \\
i_2(t) = \alpha \cdot (E_1 - E_2)^2 = I_1 + I_2 - I_1I_2\cos(2\pi f_D t) + \Delta\phi + I_{noise}
$$
因此，平衡光电探测器的输出主要与参考光和信号光之间的差频项和相位差有关。参考光和信号光产生的直流信号可以相互抵消，随机噪声也可以被显著抑制（尽管不能完全消除）。对于周期性振动且速度连续变化的目标，探测器输出的瞬时电流信号为：
$$
i(t) = 2I_1I_2\cos(2\pi \int_0^{t} f_D dt + \Delta\phi + I_{noise} = 2I_1I_2\cos(2\pi \int_0^{t} f_{AOM} dt + 2\pi \int_0^{t} f_d dt + \Delta\phi + I_{noise} = 2I_1I_2\cos(2\pi \int_0^t 2f_V(t)\cos\lambda dt + 2\pi f_{AOM} t + \Delta\phi + I_{noise}
$$
其中，$\lambda$为激光波长，$st$为振动目标的瞬时振动位移，$\Delta\phi$包括信号光相对于参考光的回波传输延迟、延迟时间内的随机相位抖动以及照射目标表面后的相位反转等因素，导致探测器端出现未知的综合相位偏移。

**3. 激光多普勒信号采集硬件平台**
本文实验选用的数据采集硬件平台是Spectrum Instrumentation公司的M4i.4421-x8开发板。该板集成了高速数据采样、A/D转换、数据记录和算术处理等功能。板上的时钟晶体振荡器基本周期为4 ns，触发采样时间选择为2的幂次倍数（8、16、32、64等），因此该板的最小触发采样时间为4 ns，最大为1.049 ms（$2^20 ns$），对应的最小采样率为953 Hz，最大为250 MHz。如图2所示，数据采样率的选择范围相对较大。

实验系统配置如下：开发板通过PCIe接口安装在主机计算机上，通过其API函数由PC应用程序控制。如图3所示，测量系统内的物理连接主要依赖光纤。1550 nm激光源为平衡光电探测器生成参考光束和信号光束，开发板接收来自发生器的耦合光电流信号和频率偏移恢复信号。振动信号的核心高精度实时分析能力主要通过软件级算法实现。

在操作过程中，选择开发板的多段记录模式进行数据采集。该模式允许用户自定义触发数据段之间的时间间隔，从而支持均匀和非均匀触发。这种模式特别适用于低占空比测量应用，因为它优化了采集内存的使用。采样过程的详细示意图如图4所示。

系统的理论测量精度主要取决于光谱分辨率的大小，光谱分辨率由下式确定：
$$
\Delta f = \frac{Fs}{N}
$$
为了实现更高的测量精度，一种方法是降低数据采样率$Fs$，另一种方法是增加数据采样点数$N$。对于每个采集段内的采样点数$N_{seg}$，较大的$N_{seg}$不仅提高了处理精度，还会相应增加计算复杂性。选择$N_{seg}$主要考虑DFT或FFT的计算负担。即使使用FFT对每个段的数据进行主要频率分析，随着$N_{seg}$的增加，频率确定的计算负担也会呈非线性增长。过多的采样点$N_{seg}$不仅会降低计算速度，还会增加采样区间内目标运动状态的复杂性。因此，为了确保测量系统的实时能力，Nseg（采样点数）不应设置得过大，通常应低于4096点。根据所选的数据采样率Fs、采样点数Nseg以及被测物体的振动频率f0，我们在数据分析过程中主要考虑两种情况：(a) 采样率Fs远大于目标物体的振动频率f0（通常大于20·Nseg·f0），则认为在信号采集的短时间内目标物体的运动状态是恒定的；(b) 采样率Fs相对于目标物体的振动频率f0较低（通常小于12·Nseg·f0），则认为在信号采集的短时间内目标物体的运动状态是变化的。对于处于这两种情况中间状态的被测物体，可以通过调整采样率Fs将其转换成上述两种情况之一进行处理。

4. 高采样率下对低频振动的分层频率分析方案
当采样率Fs远大于被测物体的振动频率f0（即NsegFs<120·f0）时，可以假设在每个信号采集时段内，被测物体的运动状态保持恒定。这意味着公式(7)中的差频fD是一个固定值。随后可以通过FFT计算直接从原始采集的信号中提取多普勒频移（后续文本中称为“主频率测量”），并将其代入多普勒频移转换公式(4)中得到速度值。接着，基于一段时间内的速度变化序列波形提取速度变化频率（后续文本中称为“次频率测量”，以此来获取目标物体的振动信息[18][19]。

关于次频率测量操作，其输入数据采样率实际上是主频率测量的输出结果率。因此，当可用于处理的采样点数有限时，需要尽可能提高主频率测量的输出率。传统方法通常采用后处理方式，即先保存所有必需的原始测量数据，然后再分段进行FFT计算，从而最大化输出率，如图5所示。然而，这种方法的显著缺点是采样率通常远远超过数值处理的处理能力。对于动态变化的目标物体，这种“先测量后分析”的后处理方法难以提供实时结果。

为了确保数据的及时性和稳定性，需要同时进行原始测量数据的采集和处理。这要求在获取一个数据段后立即进行计算和处理，然后再获取下一个数据段，每个ΔTcollect之后必须跟着一个ΔTcalculate，如图6中的数据流所示。其中，每个数据段的采样时间ΔTcollect=NsegFs=1/Δf，是原始采样数据的主频率测量分辨率的倒数。相应的数据处理时间ΔTcalculate由CPU性能、数据长度Nseg甚至数据本身的浮点数决定。由于数据采集板上的时钟晶体振荡器及其精确的触发机制，主频率测量的处理相对简单。然而，次频率测量过程的输入是非均匀的，这意味着每个单独数据段进行FFT计算所需的ΔTcalculate可能不相等。传统方法是通过高速FPGA模块生成固定周期的中断信号来触发每个数据段的采样[20]，如图4所示的均匀触发方式。

然而，很难精确设置触发间隔Ctime以满足ΔTcalculate的处理时间要求。原则上，Ctime应尽可能小。但是，过小的Ctime可能会导致实际问题：如果需要调整Nseg或其它进程突然抢占计算资源，导致ΔTcalculate超出Ctime的范围，这将导致数据流同步混乱。因此，通常的做法是将触发间隔Ctime设置得严格大于所有可能Nseg值下的ΔTcalculate。尽管如此，这种方法仍会导致计算资源浪费和处理效率降低。

为了简化并优化数据处理工作流程，每个采样数据段都配置为在当前数据段的FFT计算完成后立即激活下一个数据段的采样触发。这确保了数据流的最小延迟。然而，在这种配置下，连续数据段之间的时间间隔可能不相等。基于这种非均匀触发机制，本文提出了一种采样校正迭代DFT（SCI-DFT）方法。这种方法使次频率测量过程能够适应非均匀采样，同时也能抵御其他进程在串行计算过程中的突然中断[21]。因此，它在软件层面实现了稳定的计算输出和精确的时间同步。

4.1. 非均匀触发数据的校正处理
由于次频率测量提取的是速度序列波形的变化，而DFT/FFT处理的前提是输入数据是均匀采样的，因此需要对数据间隔进行均质化校正。具体实现方法如下：根据处理一批次频率测量输入数据（例如500或1000个数据段）所花费的时间计算平均时间间隔ΔTaver。因此，次频率测量的数据采样率为Fs?=1/ΔTaver，这意味着对于振动目标来说，主采样率限制了可测量的振动幅度范围（速度量），次采样率限制了可测量的振动频率范围[22]。获得标准化的均质化数据间隔ΔTaver后，使用每个数据段的时间戳信息进行拉格朗日二次插值，以获得均匀间隔的次频率测量输入数据，如图7所示。其中，(x', y')数组表示非均匀采样的实际数据，(x, y)数组表示经过均匀间隔计算后的预测时间点。对于每个目标点x(i)，使用最接近x(i)的三个采样点进行拉格朗日二次插值，得到相应的插值结果y(i)。在整个采样数据段上重复这个过程，即完成了采样数据的均质化校正。

对于频率为f、幅度为A的正弦信号，以采样率Fs进行采样，使用拉格朗日二次插值重建的最大绝对误差|E|可以表示为：
E～OA·f/√(3Fs)

由于目标信号的频率通常满足f?Fs，可以得出结论，拉格朗日二次插值方法有效地解决了非均匀触发采样的问题，从而在串行计算过程中具有一定的抗干扰能力[23]。

4.2. 采样校正迭代DFT方法
另一方面，次频率测量的分辨率需要与要测量的目标振动频率的数量级相匹配，即：
Fs?N??f0/2
其中N*表示次频率测量过程中涉及的数据点数量，因此有：
Ttotal=ΔTaver·N?+ΔTcalculate > N?Fs??2T0
这意味着用于分析计算的全部速度波形数据所占用的总时间必须大于目标物体的一个完整振动周期。因此，如果再次使用FFT对完整样本进行后处理分析，将会引入几个周期的延迟。为了缩短计算输出与最新采样点之间的延迟，并减少每个ΔTcalculate的计算负担，本文设计了一种次频率测量分析的迭代数据处理工作流程，称为采样校正迭代DFT（SCI-DFT）方法。该方法建立了两个专用数据缓冲序列，分别用于实时采样映射和时间预测插值。通过将数据均匀性校正与滑动窗口迭代DFT（SWI-DFT）计算相结合，该方法使串行处理流程能够以最小的计算成本输出与最新采样点同步的信号特征。

DFT的计算公式为：
X[k]=∑n=0N-1x(n)·e^(-j2πkn/N)
其中WNk=e^(-j2πNk)=cos(2πk/N-i·sin(2πk/N)，是逆时针旋转因子。考虑到次频率测量过程中每次更新都会引入一个新的采样点，计算窗口不断向前滑动（如图6上半部分所示），因此无需每次更新都计算完整的DFT，使用SWI-DFT可以显著降低计算复杂度。SWI-DFT的原理是[24]：当新的采样点x(m+1)到达时，计算窗口变为x(m+1), x(m), …, x(m+2-N)，并丢弃旧的采样点x(m+1-N)，从而使新的DFT系数变为：
Xm+1[k]=∑n=0N-1x(m+2-N+n)·e^(-j2πkn/N)=∑n=0N-2x(m+2-N+n)·e^(-j2πkn/N)+x(m+1)·e^(-j2πk(N-1))/N
在这个过程中，Xm+1[k]是时间窗口m+1的DFT系数，Xm[k]是时间窗口m的DFT系数，x(m+1-N)是从计算窗口中移除的最旧采样点，x(m+1)是添加到计算窗口中的最新采样点，WN-k=e^({j*2π*k/N}是顺时针旋转因子，表示一个采样间隔的相位旋转，k=0,1,2,?,N-1是频率索引。SWI-DFT的巧妙之处在于利用了相邻窗口之间的高相似性（重叠N-1个采样点），通过利用前一个窗口的DFT结果、被移除的旧采样点(x_old)和新增的采样点(x_new)，以及简单的相位旋转因子，直接更新当前窗口的DFT结果。这种方法将每次窗口移动后的计算复杂度从O(N^2)降低到O(N)。

虽然考虑了次频率测量过程中非均匀触发的校正以及使用SWI-DFT来提高计算效率，但仍需考虑计算输出向量与最新输入数据之间的相位匹配问题。所提出的SCI-DFT方法通过引入两个额外的缓冲区——“实时数组”和“计算时间数组”——以及特定的数据处理和传输流程来解决这个问题。这确保了DFT分析阶段（在DFT段的首个点内）与最新的速度测量样本保持完全一个周期的同步，从而实现相位同步。在这种应用场景中，SCI-DFT方法适应了采样数据的不均匀性，对每个新动态输入的数据点进行偏移预测和时间间隔补偿，同时保持每个单独数据点计算的CPU资源消耗一致。具体工作流程如图8所示。

如图流程图所示，次频率测量的实现主要包括三个阶段：
(a) 预处理阶段：在实际计算之前，选择一个数据段（在点N0之前）并进行时间戳记录，用于统计分析。在样本索引N0处计算平均时间间隔ΔTaver。
(b) 数据均质化调整阶段：从样本索引N0+1开始，作为次频率测量的起点，建立两个缓冲数组：实时数组和计算时间数组。将每个带有“Last”标签的速度测量数据副本复制到实时数组中，并在计算时间数组中用时间间隔ΔTaver标记下一个预测的样本索引Qi。从实时数组中，选择最新的实际采样点Pi，其时间范围包含下一个预测的样本索引Qi，以及其前两个数据点Pi-1和Pi-2，然后使用拉格朗日二次插值法计算对应于样本索引Qi的速度值。
(c) 迭代计算阶段：将新生成的数据点Qi处理到计算时间数组中。如果数组长度不超过设定的DFT_Length，则执行单步DFT处理；否则执行SWI-DFT处理。同时，将标记为“Cur”的数据更新为“Last”标记，然后将“Cur”标记点前进到最新的速度测量数据。这种设计有两个关键优势。一方面，它将速度测量计算和振动频率测量计算分离开，防止偶尔的停滞或错位导致级联错误。另一方面，在理想条件下，“Last”标记下的数据是DFT计算中使用的最后一个点，而“Cur”标记下的数据与DFT计算中使用的第一个点之间恰好相差一个完整的DFT周期，这使得它们具有相同的相位状态。这确保了速度结果和振动频率结果输出的同步。通过这种方式，可以尽可能少地延迟输出与最新采样值同步的振动信号分析值。此外，在工程实施中应注意，当目标振动幅度太小时，信号较弱，容易受到环境噪声的影响，可能导致误差[25]；当目标振动幅度太大时，可能会超出测量范围，导致获得的信号出现“剪切”现象。此外，由于激光源不是严格单色的，其线宽会在相干叠加的光电流中引入额外的拍频分量。在我们的实验中，参考光噪声表现为1 MHz的信号，其频率值在光路差异改变时基本保持不变。因此，可以将其视为光源的固有特性，并需要在计算多普勒频移的过程中对其进行筛选和过滤。

5. 低采样率下的自相关分析方法用于可比阶数的振动频率
当设备的检测器采样率Fs与测试中的振动频率f0处于可比阶数（2f0< /> 2f0，即在直接从检测器获取的原始信号波形内可以观察到明显的频率变化，此时不能直接应用傅里叶变换方法来解方程（7）。相反，应根据方程（8）提取目标的振动频率。
传统的解调方法主要包括模拟和解调方法。模拟解调通常使用由压控振荡器（VCO）、环路滤波器和相位检测器（PD）组成的锁相环（PLL）电路来实现[26]。然而，由于模拟相位检测器的测量范围有限，以及模拟电子元件中的元件容差和热漂移等问题，模拟解调难以满足高频或快速响应应用的需求。因此，现代LDVib技术主要采用数字解调方法。由于存在异频移fAOM，从平衡检测器直接输出的光电流信号同时包含振动相位和线性相位分量。由于线性相位每360°周期性反转，它们容易干扰振动信号的测量。因此，引入一个参考信号与中频（IF）载波信号混合，从而消除由AOM移位器生成的射频（RF）载波分量。在传统的布拉格衍射模式下，AOM的驱动信号频率正好等于诱导的光频移。因此，通常使用供给AOM的驱动信号来生成正交基带信号。这些信号与光电流信号混合后，通过滤波提取IF振荡信号，然后使用相位生成载波（PGC）技术进行相位解调[27]。
PGC解调算法主要包括差分交叉相乘（DCM）、差分自相乘（DSM）和反正切方法。其中，反正切方法因其实现简单和能够抑制幅度噪声而成为主流的PGC方法。其原理如图9所示。

下载：下载高分辨率图像（64KB）
下载：下载全尺寸图像
图9. PGC正交解调原理示意图。

将检测器输出的中频光电流信号与正弦和余弦载波信号相乘得到混合信号，使用低通滤波器（LPF）滤除2fAOM+fd的高频分量，同时保留差频分量fd，从而得到两个正交信号：
(16) I(t) = i(t)·cos(2πfAOMt)·hLPF = 2I1I2cos(4πλ·st) + Δφ
(17) Q(t) = i(t)·sin(2πfAOMt)·hLPF = 2I1I2sin(4πλ·st) + Δφ
将这两个I&Q基带信号分开，并使用反正切运算和相位展开算法解调目标振动信息。即：
(18) φs(t) = 4πλ·st + Δφ = arctan(Q(t)/I(t) + mπ
其中m表示π的整数倍，通过相位展开算法获得适当的m值以确保连续的相位输出。然后通过对φs(t)序列进行傅里叶变换来获得测量目标的位移调制信息[28]。

然而，我们认为PGC解调算法中的相位展开状态转换很复杂，依赖于高速FPGA处理。在目标进行简谐振动的情况下，可以通过基于自相关的基本频率提取方法直接计算振动频率。这种方法消除了生成正交基带信号所需的额外频率移位器，并避免了相位展开的复杂性，从而为激光振动测量技术提供了简化的实现方式。

5.1 相位调制信号的白氏级数展开
假设目标是简谐振动，其运动方程为：
(19) s(t) = A·cos(ωt + φ0)
将其代入方程（16），我们得到：
(20) I(t) = 2I1I2cos(4πAλ·cosωt + φ0) + Δφ
也就是说，通过直接将图1中所示的AOM驱动信号与光电流信号混合，可以得到一个典型的相位调制信号。与完全从检测光路径中移除AOM以创建同频系统的做法不同，对光电流信号进行数字化处理后再进行混合和滤波可以保留目标运动的方向信息，从而避免了全波整流失真。

由于方程（20）中的时变项ωt可以覆盖φ0的值，但不能覆盖Δφ的范围，因此可以忽略φ0的影响。设置模拟参数：振幅A=1μm，角频率ω=2π*500，初始相位φ0=0，数据采样率Fs = 100 kHz，每个片段采样数Nseg = 4096，激光波长λ=1550 nm，光信号相位差Δφ=π/4。将简谐振动信号与图10所示的标准化多普勒信号进行比较。仿真结果显示，多普勒信号是一个特征性的相位调制信号，其包络显示出周期性密集-稀疏的变化。

下载：下载高分辨率图像（564KB）
下载：下载全尺寸图像
图10. 简谐振动及其对应的多普勒信号示意图。

设置B=4πAλ并展开方程（20），得到：
(21) I(t) = 2I1I2cos(B·cosωt) + Δφ = 2I1I2cos(B·cosωt)cos(Δφ) - sinB·cos(ωt)sin(Δφ)
因此，核心问题是从cosB·cosωt和sinB·cosωt的形式的时间域信号中计算ω的值。自然对数e的幂级数的泰勒级数展开为：
(22) e^x = ∑_{i=0}^∞ x^i!/i!
类似地，将贝塞尔函数的原函数展开为泰勒级数：
(23) φ_a,z = e^a·(z-z^-1) = e^a·z·e^(-a^2·z^-1) = ∑_{i=0}^∞ 1/i!·(a^2)^i·z^i * ∑_{j=0}^∞ (-1)^j/j!·(a^2)^j·z^-j = ∑_{i,j=0}^∞ (-1)^ji!j!·(a^2)^i+j·z^j
合并相同的幂级数得到：
(24) φ_a,z = ∑_{n=-∞}^∞ J_n(a)·z^n
当n≥0时：
(25) J_na = ∑_{k=0}^∞ (-1)^k/k!·(a^2)^k+n!
当n<0时：
(26) J_na = ∑_{k=-n}^∞ (-1)^k/k!·a^2^(2k+n)
因此，它满足J_na=(-1)^n·J_na，J_na是贝塞尔方程的第n阶解函数。
(27) x^2/dx^2 + y/dx + x^2/α^2y = 0
其积分表示为：
(28) J_na = 1/2π∫_{-π}^π e^(a·sinx-n·x)·dx = 1/π∫_0^π cosa·sinx-n·x·dx
由于贝塞尔原函数的展开适用于所有非零复数z，令z=e^ix，我们可以得到：
(29) e^a·sin(x) = cos(a·sin(x)) + i·sin(a·sin(x)) = ∑_{n=-∞}^∞ J_na·e^i·n·x
这是复数形式的傅里叶级数展开，展开式的右侧可以写为：
(30) ∑_{n=-∞}^∞ J_na·e^i·n·x = J_0a + ∑_{n=1}^∞ J_na·[e^i·n·x + (-1)^n·e^(-i·n·x] = J_0a + 2∑_{k=1}^∞ J_2^k(a)·cos(2kx) + 2i∑_{k=1}^∞ J_2^k-1·a·sin(2k-1x)
分别对应实部和虚部，即：
(31) cosa·sin(x) = J_0a + 2∑_{k=1}^∞ J_2^k·cos(2kx)
(32) sina·sin(x) = 2∑_{k=1}^∞ J_2^k-1·a·sin(2k-1x)
将x替换为x+π/2得到：
(33) cosa·cos(x) = J_0a + 2∑_{k=1}^∞ (-1)^k·J_2^k·cos(2kx)
(34) sina·sin(x) = 2∑_{k=1}^∞ (-1)^k-1·J_2^k-1·a·sin(2k-1x)
从方程（21）的描述和方程（33）及（34）的展开形式可以看出，It信号包含直流分量以及基于角频率ω的奇次调和偶次谐波分量[29]。然而，这些谐波分量之间的能量分布与光信号相位差Δφ、简谐振动的振幅A和激光波长λ有关。特别是在Δφ=0或Δφ=π附近，奇次谐波的贡献非常小，It信号的基本频率成为振动频率的两倍。由于It信号在其谐波分量中的能量分布存在不确定性，可以使用自相关方法来分析采样信号[30]，[31]。

5.2 自相关方法的原理
对于典型的相位调制信号x(t)=cos(ωt)+b，其自相关函数定义为：
(35) R_τ = lim_{T→∞} 1/2T∫_{-T/2}^T/2 x(t)·x(t+τ)·dt
= lim_{T→∞} 1/2T∫_{-T/2}^T/2 [cos(ωt) + b·cos(ω(t+τ)]·dt
使用积分和差分公式，上述方程可以重写为：
(36) R_τ = 1/2lim_{T→∞} 1/2T∫_{-T/2}^T/2 [cos(ωt) - cos(ω(t+τ)]·dt + 1/2lim_{T→∞} 1/2T∫_{-T/2}^T/2 [cos(ωt) + cos(ω(t+τ)]·dt
当b≠0时，上述方程的第二项是一个振荡频率远高于ω的复周期函数，该项的积分贡献趋于零。因此，方程（36）可以简化为：
(37) R_τ = 1/2lim_{T→∞} 1/2T∫_{-T/2}^T/2 [cos(2sin(ωτ^2)·sin(ωt+ωτ^2)]·dt
其中，设C(τ)=2sin(ωτ^2)是一个仅与延迟τ相关且与时间t无关的常数。将方程（31）中贝塞尔函数的展开形式代入方程（37）后，我们得到：
(38) R_τ = 1/2lim_{T→∞} 1/2T∫_{-T/2}^T/2 J_0C(τ) + 2∑_{k=1}^∞ J_2^kC(τ)·cos(2kωt+kωτ·dt = 1/2J_0C(τ) + ∑_{k=1}^∞ J_2^kC(τ)·lim_{T→∞} 1/2T∫_{-T/2}^T/2 [cos(2kωt+kωτ]·dt
上述方程第二项的极限为零，所以最终结果为：
(39) R_τ = 1/2J_0·2sin(ωτ^2)
可以看出，尽管J_0x不是周期函数，但其自变量C(τ)=2sin(ωτ^2)是一个周期为Tc=2πω/2=4πω的周期函数。然而，应注意J_0x是一个偶函数，因此自相关函数R_τ的实际周期是C(τ)周期的一半，即：
(40) TR = 1/2Tc = 2πω
这意味着自相关函数R_τ的周期等于原始信号中调制余弦函数的周期。由于J_0x在变量x=0时取最大值，R_τ会在τ=n·TR (n=0,1,2,…)处显示一系列峰值，第一个非零峰值出现在τ=TR=2πω时。

当b=0时，方程（36）的第二项不能忽略，类似地，可以解为：
(41) R_τ = 1/2J_0·2sin(ωτ^2) + 1/2J_0·2cos(ωτ^2)
结合方程（25）的值，我们可以得到自相关函数R_τ的实际周期为TR*=2π/ω，这是原始信号中调制余弦函数周期的一半。从上述方程（21）和方程（33）、（34）可以看出，如果Δφ→0，调制信号仅包含偶次谐波项，计算出的信号基本频率是其实际振动频率的两倍。从理论角度来看，自相关方法可以准确提取相位调制信号的振动频率ω。此外，如方程（36）和（38）的简化基础所示，这种方法对背景噪声和恒定相位偏移Δφ的变化具有很强的鲁棒性，除了在接近零的值附近。

5.3 自相关方法的实现步骤
自相关方法的具体实现步骤如下：
(1) 采样得到s[n]（n = 0, …, N-1）。记采样间隔为Ts。
(2) 计算R[m]=∑_{n=m}^N-1 s(n)s(n-m) 对于m = 0, 1, …, N-1。
(3) 归一化，计算R[m] = R[m] / R[0]。
(4) 设置一个阈值（例如，阈值=0.5）。从m=1搜索到N/2，找到第一个满足R[m] > 阈值且R[m]大于R[m-1]和R[m+1]的局部最大点，该m值即为采样延迟m0。
(5) 如果没有找到这样的点，则降低阈值或检查信号质量。
(6) 由于测量信号的周期性，其自相关函数将在m0, 2*m0, 3*m0等处显示峰值。如果自相关曲线中的峰值位置显示出这种周期性，则记录索引m0。
(7) 对出现在m0, 2*m0, 3*m0等的峰值进行排序，其中最大值的索引表示为M0，计算周期T = M0 * Ts和角频率ω = 2π / T，该角频率对应于多普勒信号的基本频率。
(8) 由于可能存在Δφ→0的情况，奇次谐波的权重可以忽略不计，因此信号主要由偶次谐波组成。因此，测量值M0可能对应于一个完整周期或半个周期的采样点数。为了解决这种不确定性，将之前获取的数据段延迟M0/4个点，导致相位移动至Δφ′=Δφ+90°或Δφ′=Δφ+45°。这确保了Δφ和Δφ′不能同时趋近于零。然后使用延迟的数据序列s*[n]重复步骤1)到7)，以获得新的估计值M0*。通过选择M0和M0*中的较大值来确定目标在简谐运动下的最终振动频率。根据图10对应的仿真数据，并使用自相关方法进行处理，结果显示在图11中。可以观察到提取的振动频率f与预设的初始值完全匹配。此外，即使在线性信号中加入约10%的随机噪声干扰，计算结果仍然准确。下载：下载高分辨率图像（239KB）下载：下载全尺寸图像

图11. 与多普勒信号对应的自相关函数和峰值检测。

6. 实验演示

激光多普勒振动测量系统的物理构造如图3所示，选择了波长为1550纳米的单频激光作为信号源。信号发生器向AOM提供40 MHz的驱动信号，该信号也作为混合参考信号输入到开发板上。压电陶瓷振动器作为振动目标，其可配置的最大振动频率为1 kHz。PC配置如下：Intel 3.2 GHz四核处理器、8 GB RAM、1 TB硬盘容量、300 W电源和Windows 7操作系统。

6.1. 高采样率振动测量实验

在振动测量实验中，使用数据采样率Fs = 31.25 MHz和每个数据段的采样点数Nseg = 1024进行验证，为压电陶瓷设置了不同的振动频率。由于数据采样率Fs远大于所测量的振动频率f0，与第4节描述的理论框架一致，因此采用了分层频率测量方案结合SCI-DFT方法进行信号处理。实际测量结果如表1所示。

表1. 不同频率下的振动测量结果（高采样率）。

序号 | 测量振动频率（Hz） | 实际测量频率结果（Hz） | 测量误差（Hz） | 相对误差
----- | ------------- | -------------- | -------------- | ---------------
1 | 50 | 50.10 | 0.00 | 0.00 |
2 | 80 | 80.24 | 0.04 | 0.04 |
3 | 100 | 100.16 | 0.16 | 0.00 |
4 | 150 | 150.00 | 0.00 | 0.00 |
5 | 200 | 200.14 | 0.04 | 0.00 |
6 | 250 | 250.09 | 0.05 | 0.05 |
7 | 300 | 300.14 | 0.04 | 0.00 |
8 | 350 | 350.09 | 0.05 | 0.05 |
9 | 400 | 400.14 | 0.04 | 0.00 |
10 | 450 | 450.09 | 0.05 | 0.05 |
11 | 500 | 500.15 | 0.01 | 0.01 |
12 | 550 | 550.10 | 0.01 | 0.01 |
13 | 600 | 600.16 | 0.04 | 0.04 |
14 | 650 | 650.09 | 0.05 | 0.05 |
15 | 700 | 700.18 | 0.02 | 0.02 |
16 | 750 | 750.42 | 0.02 | 0.02 |
17 | 800 | 800.42 | 0.02 | 0.02 |
18 | 850 | 850.94 | 0.05 | 0.05 |
19 | 900 | 900.99 | 0.06 | 0.06 |
20 | 950 | 950.20 | 0.02 | 0.02 |

测量结果表明，在高采样率条件下，测量目标振动频率的精度约为0.2 Hz。通过处理每个数据段的时间戳信息，计算出二次频率测量数据的平均时间间隔ΔTaver ≈ 0.472 ms（包括数据读写和存储所需的时间），相应的采样率（速度刷新率）为2118.6 Hz，如图12所示。可以看出，每个数据段的“采集+处理”时间明显长于一次频率测量所需的时间ΔTcollect = Nseg/Fs = 0.033 ms。此外，记录的时间戳差异表明实际处理时间在各个数据段之间有所不同，这证实了第4节中描述的“非均匀触发”采样现象，并验证了在解决过程中使用SCI-DFT方法的必要性。此外，二次频率测量期间超过2000 Hz的数据刷新率足以满足一般工程应用的动态测量要求。

下载：下载高分辨率图像（570KB）下载：下载全尺寸图像

图12. 高采样率振动测量实验数据记录（当f0 = 10 Hz时）。

6.2. 低采样率振动测量实验

第5章描述的情景旨在实现在硬件采样限制下准确检测高频振动的能力。由于实验成本的限制，实验中使用的压电陶瓷振动器只能产生低于1000 Hz的振动频率。因此，有意将M4i.4421-x8开发板的采样频率降低，以“相对提高”测量频率比，从而测试分析算法检测高频振动的能力。原始数据采样率设置为Fs = 30.517 kHz，每个数据段的采样点数Nseg = 2048，为压电陶瓷设置了不同的振动频率。应用自相关方法进行信号处理，测量结果如表2所示。

表2. 不同频率下的振动测量结果（低采样率）。

序号 | 测量振动频率（Hz） | 实际测量频率结果（Hz） | 测量误差（Hz） | 相对误差
----- | ------------- | -------------- | -------------- | ---------------
1 | 150 | 150.02 | 0.02 | 0.02 |
2 | 200 | 200.05 | 0.05 | 0.05 |
3 | 250 | 250.09 | 0.05 | 0.05 |
4 | 300 | 300.14 | 0.04 | 0.04 |
5 | 350 | 350.09 | 0.05 | 0.05 |
6 | 400 | 400.14 | 0.04 | 0.04 |
7 | 450 | 450.09 | 0.05 | 0.05 |
8 | 500 | 500.15 | 0.01 | 0.01 |
9 | 550 | 550.10 | 0.01 | 0.01 |
10 | 600 | 600.16 | 0.04 | 0.04 |
11 | 650 | 650.09 | 0.05 | 0.05 |
12 | 700 | 700.18 | 0.02 | 0.02 |
13 | 750 | 750.42 | 0.02 | 0.02 |
14 | 800 | 800.33 | 0.03 | 0.03 |
15 | 850 | 850.81 | 0.18 | 0.18 |
16 | 900 | 900.99 | 0.20 | 0.20 |
17 | 950 | 950.27 | 0.07 | 0.07 |

从测量结果可以看出，当Fs > 100 f0时，测量误差保持在1 Hz以内。然而，随着测量频率的增加，误差有明显的放大。这是因为，在固定采样率Fs的情况下，目标振动频率f0的增加减少了每个周期的采样点数，导致离散自相关函数峰值定位的量化偏差。理论上的最大相对误差为：(42) ER = Δf/f0max = 1/2·f0/Fs。我们可以看到表2中列出的测量误差与理论公式相符。通过对自相关峰值检测过程进行拉格朗日二次插值优化后，得到的测量结果如表3所示。可以观察到测量误差通常在0.3 Hz以内，特别是对于较高的振动频率f0，算法的准确性有显著提高。

表3. 不同频率下的振动测量结果（低采样率）。

序号 | 测量振动频率（Hz） | 实际测量频率结果（Hz） | 测量误差（Hz） | 相对误差
----- | ------------- | -------------- | -------------- | ---------------
1 | 150 | 150.02 | 0.02 | 0.02 |
2 | 200 | 200.05 | 0.05 | 0.05 |
3 | 250 | 250.09 | 0.05 | 0.05 |
4 | 300 | 300.14 | 0.04 | 0.04 |
5 | 350 | 350.09 | 0.05 | 0.05 |
6 | 400 | 400.14 | 0.04 | 0.04 |
7 | 450 | 450.09 | 0.05 | 0.05 |
8 | 500 | 500.15 | 0.01 | 0.01 |
9 | 550 | 550.10 | 0.01 | 0.01 |
10 | 600 | 600.16 | 0.04 | 0.04 |
11 | 650 | 650.09 | 0.05 | 0.05 |
12 | 700 | 700.18 | 0.02 | 0.02 |
13 | 750 | 750.42 | 0.02 | 0.02 |
14 | 800 | 800.33 | 0.03 | 0.03 |
15 | 850 | 850.81 | 0.18 | 0.18 |
16 | 900 | 900.99 | 0.20 | 0.20 |
17 | 950 | 950.98 | 0.09 | 0.09 |
18 | 1000 | 1000.99 | 0.18 | 0.18 |

在测量信号相对干净且阈值设置合理的情况下（即不考虑第5.3小节中描述的步骤5）的递归调用，并且禁用了UI绘图过程时，使用自相关方法进行一次计算的实际时间约为5.73 ms。考虑到显著降低的数据采样率（与第6.1小节中的参数相比）以及第4章描述的情景下不需要二次频率测量的迭代计算，这种方法也能够满足实时动态高精度测量的需求。

7. 结论

本文介绍了LDVib系统的工作原理，并分析了在不同条件下对振动目标干扰信号进行动态实时分析的方法。实验验证表明，所提出的分析方法能够实时输出精确的测量结果，无论是对于信号发生器生成的标准波形还是实际的标准振动目标。本研究的主要创新包括：

(a) 本文提出了一个统一的数学公式，突出了使用低成本硬件和有限采样率的实际测量场景。它解释了干扰信号曲线如何反映速度变化和位移变化，并阐明了它们在不同应用场景中的不同表现。具体来说，它们包括：高采样率测量（基于速度的），其中目标状态在测量段内可以忽略不计；以及低采样率测量（基于位移的），其中目标状态在测量段内变化显著。

(b) 在基于速度的测量模式下，我们提出了一种基于分层频率测量框架的SCI-DFT方法。这种方法可以动态补偿高速采样期间的非均匀触发，最小化滚动时间窗口内的频域变换计算负载，并确保满足动态测量要求的及时数据刷新率。

(c) 在基于位移的测量模式下，测量信号的数学表达式被转换为用贝塞尔函数描述的谐波展开形式，并提出了一种基于自相关方法的处理程序来解析信号的基频。贝塞尔级数展开清晰地表示了测量信号的潜在状态，而自相关方法对随机噪声具有一定的鲁棒性。这两种方法的结合使得能够实时动态分析变化的信号。

激光多普勒测量相关技术的应用需求在全球范围内不断增长，大多数现有研究集中在激光材料、光路设计和机械结构等领域的硬件改进和优化上。相比之下，本文从数字信号处理的角度分析了基本LDVib技术的基本原理，并提供了对平衡光电探测器输出的干扰信号的数据分析方法的全面介绍。与主流振动测量解决方案相比，所提出的方法不需要额外的高速FPGA组件，具有更简单的设备结构和更低的成本，并能快速准确地计算目标的振动频率。本文详细阐述了在两种测量状态（高采样率和低采样率）下处理输出信号波形的区别和统一性。它介绍了各种技术方法的理论基础和实现程序，如分层频率测量、SCI-DFT、贝塞尔级数谐波展开和自相关方法用于计算基频。这项工作可以作为LDVib技术工程应用的有价值参考，并促进国内光电精密测量技术的进一步研究和开发。

然而，本文提出的测量算法仍然存在局限性。具体来说，其检测原理主要依赖于信号本身的显著周期性特征。因此，它仅对典型的单频信号能够获得良好的测量结果，无法准确解析非周期或多频混合振动信号。提高这种广泛的适用性仍然是LDVib技术未来探索的方向。

CRediT作者贡献声明

彭翔：撰写 - 原始草稿、可视化、验证、软件、项目管理、方法论、研究、资金获取、数据管理、概念化。
顾穆：撰写 - 审稿与编辑、概念化。
李素军：形式分析。
任启峰：验证、监督。
赵茹珍：资源、资金获取。