蒋翌超，顾劭傑，张刚，许林广，葛强，吴许强，俞本立

（1 安徽工程大学 电气工程学院，芜湖 241000）

（2 安徽工程大学 数理与金融学院，芜湖 241000）

（3 安徽大学 信息材料与智能感知安徽省实验室，合肥 230601）

0 引言

激光相干测振技术具有测量精度高、非接触测量和动态范围大等优点［1］，在引力波探测、发动机诊断、语音监听和结构安全监测等领域得到了应用［2-4］。该技术将被测物体的振动信号转化成激光干涉仪的相位变化，利用光电探测器将干涉光的光强变化转化为电信号，然后进行相位解调处理，进而提取出相位变化，获得振动信息。相位解调方案是激光相干测振技术的核心技术之一，主要包括零差解调和外差解调两种方案类型。在零差解调方案类型中，美国海军实验室DANDRIDGE A 等［5］首次提出的相位生成载波（Phase Generated Carrier，PGC）解调方案具有灵敏度高、动态范围大、线性度高和硬件开销小等优点，已成为最常用的解调技术之一［6］。但是在PGC 解调方案中，光强扰动、载波相位延迟、相位调制深度（C值）偏差和器件的非理想性能会引起非线性误差，影响解调方案的稳定性和准确性。因此，在激光相干测振应用中迫切需要一种可以抑制非线性失真的PGC 解调算法。

PGC 解调方案主要采用微分交叉相乘法（Differential and Cross Multiplying，DCM）和反正切法（Arctangent，Arctan）对正交信号中的相位变化进行提取。PGC-DCM 解调方案对C值偏差和载波相位延迟不敏感，但易受光强扰动的影响。相比之下，PGC-Arctan 解调方案可以克服光强扰动的影响，但易受C值偏差和载波相位延迟的影响。当C值偏离2.63 rad 或存在载波相位延迟时，解调出的信号会存在严重的谐波失真，甚至会出现相位解调失败情况。为了消除非线性失真，许多学者对PGC 解调算法进行了改进研究。CHEN Benyong 等［7］提出了一种融合了PGC-Arctan 和PGC-DCM 两种方案的主动线性PGC 解调算法，可以有效消除C值波动和相位载波延迟引起的非线性失真，但是系统较为复杂。常天英团队［8-9］在PGC-DCM 解调算法的基础上引入了微分自相乘和除法器，消除了光强扰动的影响。NIKITENKO A N 等［10］对PGC-Arctan解调算法进行了改进，引入了载波相位延迟补偿模块，可以消除载波相位延迟引起的非线性失真。孙抗等［11］在PGC-Arctan 解调算法中加入了微分自相除、相除和开方运算，消除了C值偏差引起的非线性失真。但是，以上算法均无法同时消除所有因素引起的非线性失真。近年来，由于椭圆拟合算法可以同时消除各种因素引起的非线性失真，被越来越多地应用于相位解调方案以改善系统的性能［12-17］。POZAR T 等［12］提出将椭圆特殊拟合算法应用于零差正交激光干涉仪以抑制系统的非线性失真，该算法计算效率高且输出的参数不会给出损坏的、非椭圆的二次曲线。严利平等［13］利用基于卡尔曼滤波的椭圆拟合算法对PGC 解调中正交分量的幅值和偏置进行最优估计和修正，减小了正弦相位调制激光干涉仪的非线性误差。畅楠琪等［14］提出了一种基于扩展卡尔曼滤波参数估计的PGC 解调方案，总谐波失真降低了19.25 dB，信噪谐波比提高了17.16 dB，系统性能显著提升。QU Zhiyu 等［15］使用基于最小二乘法的椭圆拟合算法和比例、积分和微分闭环控制模块对PGC 解调进行了改进，不仅保证系统工作在最佳相位调制深度，还进一步抑制了系统的非线性失真，输出信号的信纳比可以达到61.57 dB。HOU Changbo 等［16］利用高斯-牛顿迭代实现的椭圆拟合改进了PGC 解调算法，该算法具有良好的抗噪声和抗伴生调制能力，输出信号的信纳比波动小于3 dB。程坤等［17］将时分复用与椭圆拟合参考补偿技术应用于PGC 解调方案，实现了高分辨率的静态应变传感。然而，在被测信号幅值较小的情况下，以上改进相位解算法中的椭圆拟合均无法正常工作，甚至会导致解调失败，这限制了算法的小信号解调能力和鲁棒性。

本文提出一种面向激光测振应用的改进PGC 解调算法，该算法使用迭代重加权优化技术提高椭圆特殊拟合的精度，减小离群数据的影响。而椭圆特殊拟合是在直接最小二乘法的基础上引入椭圆约束矩阵以避免拟合结果退化为双曲线的算法，鲁棒性好，计算效率高。此外，通过低频三角波驱动压电换能器引入低频相位调制，确保在小信号情况下椭圆拟合结果的准确性，克服传统椭圆拟合算法的缺点。经过混频和低通滤波后的两路含有非线性误差的正交信号，经迭代重加权椭圆特殊拟合算法处理后被校正为一对完全正交的信号，可以显著抑制解调信号中的非线性失真。

1 原理与仿真

1.1 基本原理

面向激光测振应用的改进PGC 解调算法如图1所示，算法在高频载波调制信号上加一个低频三角波调制信号，经数模转换器（Digital to Analog Converter，DAC）输出后作用在压电换能器（Piezoelectric Transducer，PZT）2 上而促动反射镜M2，进而额外引入一个低频、大幅值的相位变化。与此同时，外部激励通过驱动PZT1 引起反射镜M1 的振动，进而模拟待测的振动信号。迈克尔逊干涉仪输出的干涉光被光电探测器（Photodetector，PD）转化为电压信号并被模数转换器（Analog to Digital Converter，ADC）采集，采集到的信号被送入改进PGC 解调算法，经过乘法器（Multiplier，MUL）、低通滤波器（Low Pass Filter，LPF）、迭代重加权椭圆特殊拟合、DCM 和高通滤波器（High Pass Filter，HPF）处理后，最终输出与激励信号线性相关的期望信号φ1（t）。

图1 面向激光测振应用的改进PGC 解调算法的示意图Fig.1 Schematic of the improved PGC demodulation algorithm for laser vibrometers

激光干涉仪输出的干涉光被PD 探测后，输出的信号可以表示为［18］

式中，A和B是与输入激光功率相关的常数，此外，B还与激光干涉仪的混合效率相关；C是载波调制深度，ω0是载波频率，θ是载波相位延迟；φ（t）=φ1（t）+φm（t）+ψ（t），包括与振动信号线性相关的期望信号φ1（t），低频调制引起的相位变化φm（t）和环境扰动ψ（t）。

ADC 采集到的干涉信号被分成两路，分别与幅值为G的基频信号和幅值为H的倍频信号混频，然后经LPF 低通滤波后，得到一对正交信号，即［19］

式中，J1（C）和J2（C）分别是一阶贝塞尔函数和二阶贝塞尔函数。

由于光强扰动、相位载波延迟、电子元器件的噪声以及C值偏差等非线性因素的存在，式（2）可以改写为［20］

式中，h和k是直流偏置，a和b是交流幅值，δ是I1（t）和I2（t）之间的相位差。

经过混频和低通滤波而获得的两路正交信号的李萨如图通常是一个椭圆弧，可以拟合为一个椭圆。以x、y代替I1（t）和I2（t），椭圆方程的一般式可表示为［21］

式中，A1、B1、C1、D1、E1和F1是椭圆拟合参数。对ADC 采集的N组I1（t）和I2（t）的数据点进行椭圆拟合，拟合得到的目标函数可以用向量内积形式表示［22］

式中，矩阵ξ=[x2，xy，y2，2x，2y，1]T，向量θ=[A1，B1，C1，D1，E1，F1]T。

式中，D为6×6 的椭圆约束矩阵，令初始权重Wα=1，θ0=［0，0，0，0，0，0］T，求解式（6）可以得到广义特征向量θ，即为椭圆特殊拟合算法的解。

迭代重加权椭圆特殊拟合算法将椭圆特殊拟合算法的解θ作为下次迭代的输入，比较θ和θ0，若两次解的差值大于设置的阈值，则将θ值赋给θ0，更新权重，然后求解更新后的式（6）。特征向量θ0和更新之后的权重Wα可以表示为［20］

式中，V0［ξα］为ξα的协方差矩阵，含有数据点的噪声信息。数据点的噪声越大，点乘的结果也越大，相应的倒数则越小。在迭代过程中，高噪声点的权重逐渐减小，降低了噪声对椭圆拟合结果的影响。若两次解的差值小于设置的阈值，则停止迭代，返回迭代重加权椭圆特殊拟合的解θ。迭代重加权椭圆特殊拟合算法是在椭圆特殊拟合算法的基础上结合了迭代重加权优化技术，其初始解是椭圆特殊拟合算法的解，然后进行迭代优化。

将式（4）中正交信号的参数h、k、a、b和δ用迭代重加权求解得出的返回特征向量θ的参数表示为［19］

接着，用式（8）计算得到的参数h、k、a、b和δ对正交信号I1（t）和I2（t）进行校正，得到一对新的正交信号，即［22］

最后，使用微分交叉相乘法（DCM）对校正后的正交信号进行处理并经过高通滤波器滤除低频调制信号和环境干扰，即可得到期望信号φ1（t）。

1.2 算法仿真

为了验证提出的迭代重加权椭圆特殊拟合算法的性能，使用MATLAB 软件进行仿真。令I1（t）=3cosφ（t），I2（t）=2cos［φ（t）-π/2］，然后在数据中加入高斯噪声以模拟式（3）中的正交信号。其中，φ（t）是一个频率为1 kHz 的正弦波信号，当φ（t）幅值分别为π/4 和π/2 时，正交信号对应的李萨如图分别为1/4 弧长和1/2 弧长。椭圆特殊拟合算法使用式（5）和（6）对数据进行拟合（每个数据点的权重Wα=1），迭代重加权椭圆特殊拟合算法使用式（5）、（6）和（7）对数据进行拟合。两种算法求得的解θ是椭圆方程的参数，代入式（5）即可得到拟合椭圆，拟合结果如图2所示。在添加σ=0.02 的高斯噪声情况下，两种算法对1/4 弧长和1/2 弧长的数据拟合结果如图2（a）和（b）所示；在添加σ=0.05 的高斯噪声情况下，两种算法对1/4 弧长和1/2 弧长的数据拟合结果如图2（c）和（d）所示。根据两种算法拟合出的椭圆和真实椭圆的偏差来判断更优的算法。由图2可知：信号φ（t）的幅值越大，李萨如图的弧长越大，两种算法拟合得到的椭圆与真实椭圆相差越小，但迭代重加权椭圆特殊拟合算法的拟合结果优于椭圆特殊拟合算法的拟合结果。高斯噪声越大，迭代重加权椭圆特殊拟合算法的拟合结果越优。此外，迭代重加权椭圆特殊拟合算法迭代3 至4 次即可得到收敛结果。

图2 两种算法对模拟数据的拟合结果Fig.2 Fitting results of the two algorithms for the simulated data

2 实验装置与结果

2.1 实验装置

实验装置如图3所示，主要由氦氖激光器（He-Ne laser）、迈克尔逊干涉仪、PD、数据采集卡（Data Acquisition，DAQ）、信号发生器（Signal generator）、功率放大电路（Power amplifier circuit）和写入解调算法的计算机（Personal Computer，PC）组成。迈克尔逊干涉仪由光阑（Aperture）、分光棱镜（Beam splitter）、聚焦透镜（Len）和两个顶部安装有反射镜的PZT 组成，其中光阑用于控制光束直径和消除杂散光对干涉光的影响，分光棱镜用于将入射激光分成光强相等的反射和透射激光，聚焦透镜用于将干涉光汇聚于PD 的探测窗口。氦氖激光器出射的激光经过光阑后被分光棱镜分成两束，分别到达安装在PZT1 和PZT2 顶端的反射镜后被反射回来。改进PGC 解调算法利用计算机控制DAQ（NI USB 6363）中的DAC 输出一个包含1 Hz 三角波和20 kHz 载波的复合调制信号，信号发生器（Rigol DG4202）输出一个正弦激励信号，功率放大电路用于放大调制信号和激励信号的电压幅值，从而驱动PZT 以产生振动。调制信号经功率放大电路放大了100 倍用于驱动PZT2，以产生低频相位调制和高频相位载波。正弦激励信号经功率放大电路放大了10 倍用于驱动PZT1，以模拟激光测振仪的待测信号。调制后的两束反射激光再次汇聚在分光棱镜处并发生干涉，最后干涉光经过光阑和聚焦透镜到达PD。PD 探测到的干涉光强被DAQ 中的ADC 采集并送入计算机，经基于LabVIEW 软件编写的解调程序处理而得到解调信号。

图3 实验装置Fig.3 Experimental setup

2.2 实验结果

2.2.1 低频调制幅值和算法影响

在不对PZT1 施加激励信号的情况下，设置DAQ 输出的载波电压幅值为0.16 V，根据文献［5］可知，通过椭圆拟合输出的参数可以计算出此时的相位调制深度（C值）为1.9 rad。调节1 Hz 三角波信号的幅值，不同低频调制（Low Frequency Modulation，LFM）幅值下混频滤波后的正交信号的李萨如图如图4所示。当LFM 的幅值分别为0.035 V、0.085 V、0.013 5 V 和0.185 V 时，正交信号合成的李萨如图约为四分之一个椭圆、二分之一个椭圆、四分之三个椭圆和整个椭圆。从图4 可以看出，改变低频三角波的调制电压幅值可以有效引入不同大小的相位变化幅值，使正交信号合成的李萨如图的椭圆弧长发生改变，进而确保椭圆拟合结果的准确性。

图4 不同低频调制幅值下的李萨如图Fig.4 Lissajous figures under different amplitudes of the LFM

使用信号发生器对PZT1 施加一个频率为1 kHz、幅值为100 mV 的正弦激励信号，与此同时，将LFM 的电压幅值设置为0.185 V。ADC 采集到的干涉信号如图5所示，图5（a）是干涉信号的时域波形图，干涉信号的最大电压为2.41 V，最小电压为0.12 V，根据（Vmax-Vmin）/（Vmax+Vmin）可以算出迈克尔逊干涉仪的条纹衬比度为0.91，具有较好的干涉效率；图5（b）是干涉信号的频谱图，从图中可以清晰地看出待测信号被转换到了载波（20 kHz）及倍频（40 kHz）的两侧边带上，与低频噪声实现了分离。

图5 PD 探测的干涉光对应的电压信号的波形和频谱Fig.5 Waveform and spectrum of the voltage signal corresponding to interference light detected by PD

接着，干涉信号分别与基频信号和倍频信号进行混频，然后经LPF 低通滤波后得到一对正交信号I1（t）和I2（t），它们的时域波形如图6所示。从图6 可以明显看出，这对信号的直流偏置h和k不为零且交流幅值a和b也不相等，在信号中存在非线性误差。

图6 低通滤波后两路正交信号的时域波形Fig.6 Time domain waveform of two quadrature signals after low pass filtering

在不同LFM 幅值的情况下，分别对正交信号I1（t）和I2（t）进行迭代重加权椭圆特殊拟合，其原始李萨如图和拟合后的李萨如图如图7所示。图7（a）显示此时LFM 的幅值为0 V，拟合前正交信号的李萨如图是一小段椭圆弧，说明待测信号的幅值较小；拟合后的李萨如图与标准圆有较明显的偏差，这是在小信号情况下椭圆拟合结果不精确而导致的结果；虽然有偏差，但是迭代重加权椭圆特殊拟合总会输出一组最优的椭圆拟合参数，不会出现解调失败的情况，算法的健壮性强。图7（b）～（d）显示迭代重加权椭圆特殊拟合前后正交信号的李萨如图，由于不同幅值LFM 的存在，迭代重加权椭圆特殊拟合校正后圆弧的弧度大小不一，但它们都与幅值为1 V 的标准圆重叠，拟合结果具有较高的精确度。

图7 不同低频调制幅值下椭圆拟合前后的李萨如图Fig.7 Lissajous figures before and after ellipse fitting under different LFM

迭代重加权椭圆特殊拟合算法校正后的正交信号被送入DCM 和HPF，输出解调信号的频谱如图8所示。从图8 可以看出在无LFM 的情况下，解调信号存在明显非线性失真，其信纳比（Signal to Noise and Distortion Ratio，SINAD）和总谐波失真（Total Harmonics Distortion，THD）分别为27.49 dB 和3.88%；加入LFM 后，当LFM 幅值分别为0.085 V、0.135 V 和0.185 V 时，改进PGC 解调算法输出解调信号的SINAD 和THD 分别为41.66 dB 和0.600%，41.60 dB 和0.593%，和42.45 dB 和0.535%。从实验结果可知，LFM 的加入可以有效改善解调算法的性能，但当正交信号原始李萨如图的图形大于半个椭圆后，SINAD 和THD 相差不大。因此，LFM 的电压幅值可以设置在0.085～0.185 V 范围内，在后续实验中，LFM 的幅值被设置为0.135 V。

图8 在100 mV 激励幅值、不同低频调制幅值下解调信号的频谱Fig.8 Spectra of demodulated signals under the stimulating amplitude of 100 mV and different amplitudes of LFM

保持PZT1 的激励电压为100 mV，基于椭圆特殊拟合和迭代重加权椭圆特殊拟合的两种改进PGC 解调算法的输出信号的频谱如图9所示。由于LFM 的幅值为0.135 V，此时正交信号的李萨如图大于半个椭圆，两种椭圆拟合算法的拟合结果均具有较高的精度。但迭代重加权椭圆特殊拟合算法是在椭圆特殊拟合算法的基础上使用迭代重加权优化技术减小了高噪声点的权重，减小了噪声对椭圆拟合结果的影响，具有更高的精度，进而可以更好地抑制非线性失真。从图9 可以看出：基于迭代重加权椭圆特殊拟合的PGC 解调算法输出信号的SINAD 和THD 均优于基于椭圆特殊拟合的PGC 解调算法，分别提升了1.99 dB 和0.27%，与仿真结果相符。

图9 两种算法的解调信号频谱对比Fig.9 Spectrum comparison of demodulated signals of two algorithms

2.2.2 解调信号稳定性测试

载波调制幅值直接和C值线性相关，PGC-Arctan 解调算法受C值偏差影响显著，当C值偏离最佳值2.63 rad 时，正交信号的交流幅值会不相等，其李萨如图会由正圆变为椭圆，进而导致解调结果中出现严重的谐波失真［8］。从式（2）可知，相位载波延迟也影响着正交信号的交流幅值，从而影响解调结果。利用迭代重加权椭圆特殊拟合实时对正交信号的幅值进行归一化，消除C值偏差和相位载波延迟等对解调结果的影响。此外，还采取了DCM 方法从校正后的正交信号中提取期望信号，使得算法对正交信号交流幅值的波动更加不敏感。将信号发生器的激励信号幅值从100 mV 增加到400 mV。接着，调整载波调制信号的幅值，以0.2 rad 为间隔逐渐增大C值，使其由0.8 rad 增大至3.4 rad。不同C值下正交信号的李萨如图和拟合校正后正交信号的李萨如图如图10所示。图10（a）是正交信号I1（t）和I2（t）的原始李萨如图，此时C值分别为1.4 rad、2.0 rad、2.6 rad 和3.2 rad，随着C值增加，李萨如图由椭圆变为正圆，然后再变为椭圆；从图中也可以看到当C值为3.2 rad 时，椭圆向左侧倾斜，表明两路信号的相位差不是严格的π/2。图10（b）是经过迭代重加权椭圆特殊拟合算法校正后的李萨如图，由图可以看出不同C值下的李萨如图均被校正为一个完整的正圆，这说明正交信号的直流偏置被消除，交流幅值被归一化，相位差被校正为π/2。实验结果表明：迭代重加权椭圆特殊拟合算法可以有效消除正交信号中由于C值偏差和相位延迟引起的非线性误差。

图10 不同C 值下正交信号的李萨如图Fig.10 Lissajous figures of the quadrature signals with different C values

不同C值下输出的解调信号的SINAD 和THD 如图11所示。改进PGC 解调算法在不同C值下SINAD 的平均值为42.99 dB，对应的标准差为0.55 dB；THD 的平均值为0.44%，对应的标准差为0.03%。实验结果表明：提出的改进PGC 解调方案具有优秀的稳定性，可以有效消除C值偏移和载波相位延迟的影响，在0.8～3.4 rad 的C值范围内，SINAD 和THD 波动很小。

图11 不同C 值下解调信号的SINAD 和THDFig.11 SINAD and THD of demodulated signal under different C values

2.2.3 动态范围、频率响应和振动信号实测

图12 系统的响应线性度、底噪及THD=1.00%情况下解调信号的波形和频谱Fig.12 Response linearity and noise floor of the system and waveform and spectrum diagrams of the demodulated signal when THD=1%

在载波频率为20 kHz 的PGC 解调系统中，理论最大可探测信号的频率为10 kHz。为了防止信号混叠，将信号的最大可探测频率设计为8 kHz，因此LPF的截至频率设置为8 kHz。考虑到系统中HPF 的截止频率为20 Hz，系统可探测信号的频率范围为20～8 000 Hz。设置载波信号幅值、LFM 幅值和PZT1 上激励信号幅值分别为0.16 V、0.135 V 和100 mV，调节PZT1 上激励信号的频率并记录下解调信号的幅值，实验结果如图13所示。由图可知：在1 000～6 000 Hz 的频率范围内系统的输出幅值波动较小；在6 000～8 000 Hz 的频率范围内，由于受到LPF 过渡带的影响，系统的输出幅值有较小的下降。此外，对两种椭圆拟合算法的运行时间进行了800 次连续测试，它们运行时间的均值分别为0.010 s和0.031 s，与仿真结果中的迭代次数基本相符。迭代重加权椭圆特殊拟合算法虽然延长了计算时间，但是0.021 s的延时对信号的实时解调影响不大，仍能满足动态信号的解调需求。

图13 系统输出幅值的频率响应Fig.13 Frequency response of the system′s output amplitude

关闭驱动PZT1 的激励电压，先后对M1 施加振动信号1 和振动信号2，系统解调出信号的波形和频谱如图14所示。对比图14（a）和（b）中解调信号的波形，可以看出两个振动信号存在明显差异，振动信号1 的幅值和周期小于振动信号2 的幅值和周期。对比图14（c）和（d）中解调信号的频谱，可以发现振动信号1 包含更多的高频成分而振动信号2 存在较多的低频成分。实验结果表明，提出的改进PGC 解调算法结合低频调制和迭代重加权椭圆特殊拟合算法消除了正交信号中的非线性误差，提高了系统响应的稳定性、线性度和动态范围，具备解调不同幅值和频率的振动信号的能力，可以较好地恢复出振动信号。

图14 系统解调出的两种振动信号的波形和频谱Fig.14 Waveform and spectrum diagrams of the two demodulated vibration signals of the system

3 结论

本文提出了一种面向激光测振应用的改进PGC 解调算法，利用低频调制技术和迭代重加权椭圆特殊拟合算法改善PGC 解调算法的性能。低频调制的加入确保了算法在小信号场景下正常工作，迭代重加权椭圆特殊拟合算法减小离群数据对椭圆拟合结果的影响，提高了椭圆特殊拟合的准确度。实验结果表明：低频调制电压大于0.085 V 可以有效保证拟合结果的准确性；改进的PGC 解调算法在不同C值下（0.8～3.4 rad）测得的平均SINAD 和THD 分别为42.99 dB 和0.44%，对应的标准差分别为0.55 dB 和0.03%，输出解调信号具有高稳定性；激光测振实验系统的响应线性度优于99.99%且动态范围达到了103.90 dB @ 500 Hz &THD=1%。提出的改进PGC 解调算法精度高、鲁棒性强且计算效率高，在激光测振领域具有良好的应用前景。