一种基于残周期拟合的激光测冲击信号解调方法

梁志国

(北京长城计量测试技术研究所 计量与校准技术重点实验室，北京 100095)

1 引 言

激光测量冲击与激光测振[1～5]，均主要利用相对运动产生的激光多普勒效应，因而，调频(FM)信号解调一直是该方法的核心。其中的关键有两点：解调准确度和解调速度。在计量校准中，最重要的是解调准确度。

关于调频信号的解调，目前存在两类方法，硬件法和软件法。无论哪类方法用于激光测量，人们都不能直接在光频上进行解调操作，均需要通过变频方式将激光频率中的调制信息变换到较低的无线电频率上，再执行后续解调和处理。

经典的硬件解调法，主要通过模拟电路，鉴频、检波解调出频率调制信号波形。

软件法主要在低频的基带信号基础上，使用高速数据采集技术进行波形测量，获得数字化测量序列，然后再以软件算法计算出调制信号波形，完成解调。软件法中，比较常用的是过零点检测法，正交变换法、希尔伯特变换法、曲线拟合法、跟踪滤波器法等[6～12]。

其中，用于外差FM信号的过零检测法最为简单，缺点是时间分辨力较低，且解调结果的时域坐标为非线性的。因而仅仅适用于获得调制波形的某些参量，在获得低失真调制波形方面存在原理缺陷。

正交变换法是最重要的解调方法，它借助于两路相位相差90°的正弦信号，在与载波含有相同的频率调制信息时，通过对两路信号的反正切运算获得已调制信号的实时相位波形，经过相位展开获得连续相位曲线，该曲线的微分即是频率调制信号的解调结果。

关于两路相位相差90°的载波信号的制备与获取，有光学法、电路法、软件法、希尔伯特变换法等。由此产生了正交变换法频率解调的不同技术方案，分别适应于零差式激光干涉法测振信号解调和外差式激光干涉法测振信号解调。其缺点是需要制备两路瞬时幅度均衡的已调FM信号，并且需要后续滤波和微分环节，才能得到频率解调结果。

正弦波拟合法用于FM信号解调，可避免两路信号的制备问题，无需滤波与微分即可获得频率解调结果，对于外差式激光干涉信号的解调有较高的准确度和稳定性。其不足是，用曲线拟合实现解调，计算量巨大，耗时很长。另外，面对冲击信号的激光差动干涉测量序列时，其没有明确的载波频率，起始频率接近直流，以此跃升到较高频率，解调难度较大，其前几个震荡波形的形状往往并不规则，与一个完整的正弦周期波形相去甚远。使用近似一个周期的正弦波形进行解调运算，在起始点处不易获得良好结果，需要对原始测量波形进行截取后才能进行解调运算。

针对上述问题，本文提出一种基于残周期正弦拟合的FM信号解调算法，使用约1/5左右的正弦波形周期模型进行FM信号的解调运算[13]，以增加解调算法对起始端信号波形不规则状况的适应性，同时，由于残周期相对于整周期所用的波形点数减少为原来的1/5，可望缩短解调运算时间。

2 残周期拟合法FM信号解调原理

用四参数正弦曲线拟合法进行正弦载波调频信号的解调，基本思想为，调频信号是由一段段频率随调制信号变化的正弦波组成，每一小段独立运算的波形长度作为解调窗口长度。通过滑动正弦波模型，对每一段曲线的正弦波模型频率进行估计，获得解调的调制信号。

关于解调，越窄的解调窗口将能越真实地反映调制信号的瞬时状况，而越宽的解调窗口将因为平均滤波效应，降低解调灵敏度；但解调窗口的减窄，同时将使模型参数的拟合误差增大，影响解调准确度；另外，两个相邻估计值的最短时间间隔，将决定解调信号的时间分辨力，人们也希望它越短越好。

这里，使用1/5个周期左右的信号点数进行拟合；首先，对第1个信号点开始的约1/5个周期的信号的模型参数进行估计，将它作为该段数据中心点的模型参数；然后，以该组估计参数为初始值，对中心点后面一点为中心的约1/5个周期的信号的模型参数进行估计；依次类推，直至最后一个数据点，结束估计。

其后，对获得的模型参数序列进行波形分析，获得频率调制信号。

首先，介绍一种三参数正弦曲线拟合算法和由它改造获得的四参数正弦曲线拟合算法。

2.1 三参数正弦曲线拟合法[14]

设理想正弦信号为：

y(t)=A0cos(2 π ft)+B0sin(2 π ft)+C0

=Acos(2 π ft+θ0)+C0

(1)

数据记录序列为时刻t1,t2,…,tn的采集样本y1,y2,…,yn，采集速率v，采样间隔为Δt，ti=i×Δt =i/v，(i=1,…,n)，数字角频率ω=2 π f/v，则式(1)可表示成下列离散形式：

y(i)=A0cos(ω·i)+B0sin(ω·i)+C0

(2)

三参数正弦曲线拟合过程，即为ω已知，寻找A1、B1、C，使残差平方和ε最小：

(3)

则，参数A1、B1、C即为A0、B0、C0的最小二乘拟合值。为寻找出A1、B1、C，构造矩阵：

式(3)可用矩阵表示如下：

ε=ε(ω)=(y-D0x0)T(y-D0x0)

拟合函数如下：

其幅度和相位表达形式：

(4)

式中：

拟合残差为：

ri=yi-A1cos(ω·i)-B1sin(ω·i)-C

(5)

由于这是一种闭合算法，因而收敛是肯定的。

2.2 残周期四参数正弦曲线拟合的绝对收敛算法

三参数正弦曲线拟合是一种闭合的线性过程，不存在收敛问题，而四参数正弦曲线拟合则不然。多年以来，人们总希望能找到一种判据，当满足判据要求时，四参数正弦波曲线拟合过程绝对收敛。下述过程具有该特点。

假设待估计的正弦波数字角频率目标值为ω0，其采样序列所含信号周期个数为p；则有，Δωmax=ω0/p，在区间[ω0-Δωmax，ω0+Δωmax]内的任意角频率ω下，残差平方和ε(ω)的极值存在且唯一，这样，可将四参数正弦曲线拟合中，对幅度、频率、相位、直流分量4个参数的四维非线性搜索，转成对频率分量ω造成的ε(ω)的一维线性搜索，可保证在区间[ω0-Δωmax，ω0+Δωmax]内，用三参数拟合法实现的四参数正弦曲线拟合过程绝对收敛。

当所用的数据在一个波形周期以下时，若波形长度为T0，则波形周期T≥T0，波形频率f≤1/T0。则上述四参数拟合的收敛区间变为(0,4 π/(v·T0)]，这是一个非常宽的收敛区间，很容易工程实现。由此，获得残周期正弦波四参数拟合过程如下：

(1) 设定拟合迭代停止条件为he；

(2) 从正弦波测量序列yi(i=1,…,n)中，截取不足一个波形周期的残周期正弦波测量序列yi(i=1,…,m)，直接计算波形周期估计值T0=m/v；

使用黄金分割法确定中值频率：

否则，重复(4)～(6)的过程，进行下一次迭代。

3 激光测振法用于冲击加速度的测量过程

3.1 测量过程

如图1所示，典型的冲击加速度校准装置主要由Hopkinson冲击机、差动式光栅激光干涉仪、数字示波器，以及计算机等设备组成。

图1 冲击加速度校准装置框图Fig.1 Primary Shock Calibration Unit of Acceleration

1) 使用Hopkinson冲击机发射弹体，利用Hopkinson棒在校准端面产生波形良好的冲击加速度，直接加载到光栅和被校准的加速度计上。

2) 使用差动式激光干涉仪，利用光电探测器检测光栅运动产生的激光多普勒信号y(t)。

3) 使用数字示波器对激光多普勒信号y(t)及经信号适调仪调理后的被校加速度计输出信号进行数据采集。获得y(t)的采样信号yi(i=1,2,…,n)。

4) 使用上述调频信号解调方法，解算出多普勒频移信号Δf(i)=fi(i=1,2,…,n)。

使用图1所示差动式衍射激光干涉仪进行冲击测量时，其运动速度波形信号v(i)[7]：

(6)

式中：d为光栅的栅距；p、q分别为差动激光干涉仪所使用的光栅衍射级数。

5) 按式(6)计算冲击速度序列v1,v2,…,vn；它是一个典型的阶跃波形。

6) 用众数法计算冲击速度波形上升时间tr；[15]

7) 对上述速度波形微分得到冲击加速度a(i)：

(7)

(8)

3.2 判定峰值最优

为判断L取值是否最优，采取如下判定方式：

在第i点vi处，将平均值窗口内的速度值的邻域序列vi-L,vi-L+1,…,vi+L-1,vi+L,做端基直线拟合，得：

v(k)=(k-i+L)·g+v0

(9)

(10)

以局域平均窗口中的误差均值Δv(t)曲线作为微分的误差指针，当窗口变窄时，Δv(t)曲线在加速度峰值附近的量值呈下降趋势，且有明显的局部极小值，当该极小值下降不再明显，且局部极小值特征微弱时，即可判定窗口选择值L已经接近最佳。

4 实验验证

4.1 基本操作

在如图1所示的冲击加速度校准装置上进行校准实验。其中：数字示波器型号为TDS544A，带宽500 MHz，8 Bits A/D，双通道采样工作方式，垂直电压设置为0.1 V/div，通道采样速率500 MSps，采集数据点数n=231 001。差动激光干涉仪所用光栅栅距d=1/150 mm，式(6)所用衍射级数分别为p=1，q=-1。

用图1所示校准装置，对某型冲击加速度传感器进行校准实验，获得冲击过程激光多普勒干涉信号如图2所示。

图2 冲击过程中的激光多普勒干涉信号Fig.2 Laser Doppler Interference Signal of Impulse

将图2所示测量曲线按上述方法处理，获得如图3所示的冲击速度波形。该冲击速度波形的上升时间为[15]：tr=62.3 μs。

图3 冲击过程产生的速度波形曲线(tr=62.3μs)Fig.3 Velocity waveform of impulse

选取不同宽度L的局域平均窗口，按照第3.1节 7)所述计算加速度波形如图4～图10所示。它们的冲击峰值最大值和最小值相差18%。由此可见，局域平均窗口宽度的变化对加速度峰值的影响不可忽视。比较不同传感器特性时，使用相同宽度的平均值窗口执行运算操作才更有意义和价值。

图4 冲击波形曲线(L=tr/(4×Δt))Fig.4 waveform of impulse(L=tr/(4×Δt))

图5 冲击波形曲线(L=tr/(6×Δt))Fig.5 waveform of impulse(L=tr/(6×Δt))

图6 冲击波形曲线(L=tr/(8×Δt))Fig.6 waveform of impulse(L=tr/(8×Δt))

图7 冲击波形曲线(L=tr/(10×Δt))Fig.7 waveform of impulse(L=tr/(10×Δt))

图8 冲击波形曲线(L=tr/(12×Δt))Fig.8 waveform of impulse(L=tr/(12×Δt))

图9 冲击波形曲线(L=tr/(16×Δt))Fig.9 waveform of impulse(L=tr/(16×Δt))

图10 冲击波形曲线(L=tr/(20×Δt))Fig.10 waveform of impulse(L=tr/(20×Δt))

4.2 峰值判定讨论

表面上看，平均窗口宽的图4滤除了更多的细节特征，从而使得冲击曲线更加光滑，而平均窗口越窄的图则保留了更多的细节特征，图6～图10中甚至出现了子母峰现象。但是，由于物理机理掌握尚不够深入，很难判定图5和图10哪一个更加接近真实的冲击加速度波形。这一直是冲击计量校准中最难以解决的问题之一。

尽管如此，从图10(a)和图5(a)相比，可以看出，仅仅由于滤波参数的选择，峰值变化可以达到18%之多。

本文中，为了判定滤波窗口的优劣，特别引入窗口内局部速度偏差均值曲线Δv(t)，从图4(b)～图10(b)可见，随着滤波窗口的变窄，加速度峰值附近的局部小区间内的Δv(t)呈下降趋势，且存在局部极小值特征，若极小值底部比较尖锐，则表明平均值窗口宽度过宽，应减小窗口宽度，直至出现比较平缓的局部极小值特征，或者直到某一宽度以后，局部极小值便稳定在一个范围内波动，局部极小值特征不再明显，这时，可以认定，再减小滤波窗口已经无意义了，此时的窗口接近最佳。

令横坐标为上升时间tr与滤波窗口宽度2L之比r，则有图11所示的Δv(t)-r关系曲线图。

图11 Δv(t)-r关系曲线图Fig.11 Δv(t)-r Curve

由图可见，随着r值的增大，Δv(t)总体呈下降趋势，但下降速度到一定程度后趋于稳定而不再下降，本实验中，r=6以后，局部速度偏差Δv(t)变化趋稳，此时，加速度变化也趋于缓慢。可认为获得的加速度复现性较好，量值比较准确可靠。

通常，认为曲线急速下降至较为平缓的分界拐点，是滤波窗口的特殊取值点。本实验中，r=6为分界拐点，此时窗口宽度为上升时间的1/6。

上述实验所获冲击加速度波形如图12所示。其加速度测量序列，按ISO标准方法处理，获得峰值为1 048 301.679 8 m/s2，此时，Butterworth滤波器截止频率选取10/T，T为图12所示脉冲的脉宽。

图12 冲击过程产生的加速度波形曲线Fig.12 Acceleration waveform of impulse

比较平均值窗口取上升时间1/6的图8与图12所示曲线峰值[16,17]，两者相差为7%，是一个比较大的值，可见，此时，本文方法的分界拐点与ISO推荐的方法的差异主要体现在子母峰现象上。

而从图10(a)和图12相比，可以看出，本文所述方法与ISO标准所列方法可以有近11.7%的峰值差异。由于大多数滤波器对于冲击峰值的作用都是将峰值幅度降低，因此，有理由相信，本文方法获得的子母峰现象是真实的，且所获得峰值应该比图12所示用ISO标准获得的量值更接近真实峰值。

5 讨 论

FM信号解调属于基本工程问题，在无线电专业领域，往往有一个远高于调制频率的载波频率，它使得FM信号的频率变化看上去比较缓慢，且相对于载波而言，仅仅在一个比较小的范围内变化，因而解调参量的变化幅度较小，容易实现较高质量的解调。而在差动式激光测振过程中，载波频率已经被消掉，实际上获得的波形是一个低频端理论上接近于直流，而高频端的频率数十MHz量级的变频信号。从中将瞬时频率提取出来的过程，要求解调系统能适应从极低的频率一直到很高的频率。实际测量获得的激光多普勒冲击信号波形，如图2所示，在起始段往往并不规则，除了频率变化外，幅度波动、波形不平稳、测量噪声等都将对解调产生影响，甚至导致解调失败。多数情况下，需要人工进行测量波形的预处理，包括切除起始段和末尾段不规则部分波形，滤波、以及剔除非平稳趋势和进行幅度归一化等，然后再执行瞬时频率的解调，难度更大。

本文所述方法，以残周期正弦模型滑动方式进行的FM信号解调，本意即是用残周期正弦波模型适应所解调的FM波形的各种变化条件情况，试图不进行滤波、截取、归一化、剔趋势等操作，即可以在任何条件下实现自适应瞬时频率的解调估计。实际过程已经实现了该目标。在实验中遇到的问题主要有：1)滑动模型法解调FM波形属于自适应变参数波形拟合法，因而在较低频率点上拟合数据量非常大，每个周期有成千上万个测量点，占用了大量时间，而所获得的信息又最少；在高频点上，每个周期采样点数较少，可能仅有十几个点，导致以残周期法进行拟合可能会有运算点不足的问题。解决方式可以用拟合点数上下限判据进行，多于上线点时，以上线点代替，少于下限点时，以下限点代替。

另外，本文所用残周期拟合算法在0～2个波形周期情况下均能保证收敛[13]，故可以使用固定长度的数据段进行滑动拟合解调，并能保证给出正确结果。但需要预先估计一下所解调的FM信号波形最短的波形周期中每个周期所含采样点数，并据此设定解调数据段长度。

本文上述方法的弱点是计算量庞大，所花费的运算时间较长，以上述实验为例，231 001个测量点的FM波形序列解调过程花费346 min又26 s(3 GHz 主频，双核CPU，1.96 GB内存)。因而该方法仅适合用于计量校准和实验数据的精确分析，不适合实时测量。

由FM信号经解调后获得的信号波形仅是速度波形，若想在此基础上获得加速度波形，微分是必不可少的。若不进行数字滤波，以速度波形直接进行差分获得的加速度波形，受噪声、失真等影响极大，很容易得出错误结果。对此，ISO标准给出了巴特沃斯滤波器及阶次和参数的推荐值，但它们对微分带来的影响一直不清楚，本文从另外一个角度，以多点直线拟合方式，给出了微分后的加速度波形，并以变窗口方式展现了滤波窗口对微分后加速度信号波形的影响。这也是高分辨力FM解调带来的益处。伴随着滤波窗口的减小，冲击加速度波形出现了子母峰现象，这也是以前使用推荐的巴特沃斯滤波器未能发现的现象。以至于多宽的滤波窗口是最佳的问题得以被讨论和分析。可以看出，最佳滤波窗口的选取，直接与速度波形的噪声、失真等有关，它们的量级水平直接影响滤波的效果以及加速度波形的质量。

6 结 论

综上所述，本文提出了一种为增加解调适应性和自动化程度、降低运算量和运算时间，以残周期正弦波模型滑动拟合实现的FM信号解调方法，并讨论了冲击加速度计量校准中的微分过程及滤波问题，提出了一种基于阶跃上升时间为关联参数的局域平均方式，同时实现数据滤波和数值微分，具有过程简洁，不确定度明确的优点，并且具有局部寻优的条件判据。对于具有子母峰形状的超高冲击波形依然适用。理论与实践均证明了所述方法的有效性与可行性，可用于冲击加速度的计量校准。