基于DIC实验的空间域信号去噪方法研究

朱海华，孟 帆，宋汉文

（同济大学 航空航天与力学学院，上海200092）

数字图像相关方法（Digital image correlation，DIC）作为一种视觉测量技术，通过相机系统采集物体在运动过程的图像，以结构表面预处理的散斑为信息载体，通过图像相关计算以及立体匹配，从图像坐标重构得到物理点的世界坐标。DIC相对于传统接触式测量方法而言，无需布置传感器，不会引入附加质量及刚度效应，且可实现全场同步测量。Phillip L 等[1]对比了基于DIC 方法及激光多普勒测量仪所辨识得到的模态参数，验证了DIC 方法在振动测试中的有效性及优越性。因DIC 具有非接触式、高空间分辨率测量的优势，国内外已有一定量的研究将其应用到动力学相关领域[2-5]。

自然环境激励技术（Natural excitation technique，NExT）于1992年被美国SADIA 国家实验室的James等[6]提出，指出白噪声激励下系统位移响应的相关函数与脉冲响应函数在表达上具有相似形式，皆为系列衰减谐波的叠加，包含了相同的模态频率、阻尼比及振型。因此可以利用相关函数代替脉冲响应函数代入到经典的时域辨识算法中，进行环境激励下的模态分析。目前，对于工况模态分析已有大量的相关理论及应用研究[7-8]。

DIC 在振动测试的应用中，通常相机系统与激励系统是独立的，二者之间非同时触发，激励信号与图像采样频率亦不相同，这会导致系统输入与响应之间并非一一对应。因此，可完全基于DIC 测量所得响应信号进行工况模态分析，得到结构模态参数。DIC 是对变形前后的图像划分子集进行相关计算，输出结果为每一帧图像所对应的结构表面变形场，即全场的面响应。而DIC 方法本身受相机本底噪声、光照、散斑等方面的影响，其计算结果所包含的噪声不可忽略。因此，需先对DIC 直接输出的时序图像所对应位移场结果进行空间滤波，以消除空间域噪声的影响。本文通过DIC获得薄板结构在白噪声激励下时序的面响应，并通过空间滤波消除二维信号中的噪声影响，重构位移场信息。再从时域信号出发，从响应的互相关函数出发进行工况模态分析，并将实验结果与有限元计算模态结果进行对比。

1 复模态下的自然激励技术[8]

对于一个N 自由度线性定常系统强迫振动，其运动的微分方程可表示为

其中：[M]、[C]及[K]分别表示系统的质量矩阵，阻尼矩阵以及刚度矩阵；{x(t)}为系统稳态响应；{f(t)}为系统所受强迫激励。

假设系统在k点受到白噪声激励力{fk(t)}，则在第i测点上的响应为[9]

其中：{ψ}r(r=1,2,3,…,N)为系统的第r阶模态振型 ；λr为系统的第r阶特征值；diag(a1…ar…aN a*1…a*N)为系统矩阵在2N空间中复振型加权正交所得矩阵；*表示共轭运算。

则第i,j测点间响应的相关函数Rkij(T)为

若系统在L个点上分别受到互不相关的多点白噪声激励，则第i，j测点间的相关函数可写为

式中：Θjr表示为

其中：αk是和激励强度有关的常数。

单位脉冲响应函数表达式为

通过观察式(4)和式(6)可以发现，相关函数的表达式与脉冲响应函数具有相似的表达形式，均为衰减谐波的叠加，且包含相同的模态频率、阻尼比及振型信息，不同之处在于模态参与因子。因此，可利用响应的互相关函数代替脉冲响应函数，代入到ERA算法中进行模态参数辨识。

2 空间信号去噪的基本原理

对于DIC 实验结果而言，相机系统所引入的噪声不可避免。传统测量方式的噪声是包含在单个测点的时域响应信号内，而DIC 的直接输出结果是每个时刻的位移场，噪声分布在空间域内，如图1中所示。因此，需对空间信号进行去噪预处理，再由时域信号出发进行工况模态分析。

图1 DIC结果实例

本节通过对位移场的空间傅里叶变换得到二维信号的波数域形式，在波数域中对其进行低通滤波，以实现对位移场的去噪重构。由FFT算法可知，图1中(b)的高频成分，不仅包含噪声，同时又存在因信号非周期边界条件的影响所引入的吉伯斯效应。本文提出改进的傅里叶连续延拓方法（Fourier extension）[10]，延拓二维信号并构造周期边界条件，消除了吉伯斯效应影响。再通过波数域的窗函数实现空间低通滤波，对位移场进行去噪重构。

2.1 傅里叶连续延拓

对于空间非周期信号f()x,y，其分布区间为J1=[0,ax]×[0,ay]上，现将该非周期信号进行空间域上的延拓使其为周期信号w()x,y，延拓后的区间为J2=[0,bx]×[0,by]，其中bx=2ax，by=2ay。将延拓后的周期信号用傅里叶指数进行拟合，其数学表达式为

其中：当Mx为偶数时，为奇数时，；当My为偶数时，My为奇数时，。而未知系数Amn可通过原区间上的已知信号得到最小二乘意义下的解：

其中：xp,yq分别表示x,y方向上的第p,q个坐标，f(xp,yq)表示区间J1=[0,ax]×[0,ay]上的原信号。将式(8)以矩阵形式表示，f(xp,yq)及Amn重构为一维向量：

则式(8)表示为求解方程[B]{α}≈{f}在区间J1=[0,ax]×[0,ay]的解，其中矩阵[B]中的元素为

通过对矩阵[]B进行SVD 分解，则系数矩阵重构后的一维向量的解为

其 中：[V]是MxMy×MxMy的酉矩阵，[S]†是MxMy×axay的对角矩阵，表示[S]矩阵的伪逆表示[U]矩阵的共轭转置矩阵。当系数矩阵求得后，利用二维傅里叶级数进行拟合并延拓原信号至区间J2=[0,bx]×[0,by]，重构后满足周期边界条件。

现通过MATLAB 给定一个理想二维信号f(x,y)=ex+y进行示例，其分布在区间[0,1]×[0,1]。使用上述傅里叶指数进行拟合并延拓，在x和y上的拟合阶数为Mx=My=25，将信号延拓至区间[0,2]×[0,2]。原信号与延拓后信号如图2所示。

2.2 波数域窗函数

在消除由于吉伯斯效应所引入的高频影响后，对位移场中的噪声进行低通滤波。将位移场通过空间傅里叶变换，转换到波数域，数学过程如式(13)所示：

式中：其中ℱ表示空间傅里叶变换算子，kx,ky分别表示在x和y方向上的波数，W(kx,ky)表示空间域信号w(x,y)的波数域表达。

本文在采用Li等[11]改进的波数域滤波方法进行低通滤波，再逆变换得到重构后的位移场，滤波窗函数表达式如下：

图2 延拓信号示意图

图3 空间低通滤波窗函数

2.3 空间域滤波信号的重构

上述内容阐述了本文对于空间信号的延拓方法与波数域滤波所涉及到的基本原理。通过对DIC所测位移场进行二维傅里叶延拓消除吉伯斯效应所引入的波数域高频影响，再通过波数域的低通滤波对位移场进行去噪重构，结果如图4所示。

3 基于数字图像相关的工况模态分析

3.1 实验装置及流程

实验对象选用一块纯铝薄板，长440 mm，宽300 mm，厚1 mm，密度为2 700 kg/m3，泊松比为0.3，杨氏模量为70 GPa。实验装置如图5所示。使用两块63 mm×40 mm×5 mm 标准角钢通过螺栓夹持铝板，再与隔振平台相连接，以此构成单边固支的边界条件。除去夹持段，铝板实际尺寸为400 mm×300 mm。本实验采用激振器作为激励方式，通过功率放大器与信号发生器相连接，激振器使用弹力绳悬挂于特制金属框架上。

实验采用双目高速相机系统进行图像采集，两台相机通过数据传输线与工控机相连接，工控机负责图像的存储。双目相机通过信号同步仪实现同时触发采集，保证两台相机所采集图像的同步性。考虑相机所拍实际视场以及所关心铝板的频率段，调整相机的分辨率以及拍摄帧率。本次实验采用图像分辨率为600 pixel×800 pixel，拍摄帧率为450 Frame/s。实验采用两台1 500 W 摄影灯作为可调光源。

图4 空间去噪示意图

图5 实验装置图

实验前使用Speckle generator 软件进行散斑生成，考虑散斑点大小在图像中所占大小为3～5个像素平方，故散斑直径大小统一为2.5 mm。散斑图纸生成后，通过纤薄纸张打印再均匀、紧密地粘贴于被测铝板表面。

在实验开始，首先对双目相机进行标定获得相机内外参数。记录下试件在静止状态下的图像，作为后期进行图像相关计算的参考图像。实验采用的输入信号为稳态白噪声，经由功率放大器传递给激振器。为减少激振器所引入的附加质量及刚度效应，激励点选择在靠近铝板夹持处的角点。双目相机在同步控制仪的外触发下，共采集10 000 帧作为参考图像，参考点响应时程如图6所示。

图6 参考点响应

图像相关计算中，所划分的子集大小为21 pixel×21 pixel，计算步长为7 pixel，对于图像中的ROI 共划分为3 975 个计算点，其中在x方向上53列，y方向上75 行。选取自由端角点处附近处测点作为参考点，全场测点与其计算互相关函数，见图7，代入ERA算法进行模态参数辨识，实验数据流程如图8所示。

图7 互相关函数

图8 实验数据处理流程示意图

4 实验结果

分别对DIC直接结果与经过空间去噪的时序位移场，提取测点时域响应信号。通过互相关函数代替脉冲响应函数，代入到ERA 算法中，使用模态置信准则（MAC）剃除伪模态之后辨识得到被测薄板的前5 阶模态参数，固有频率与阻尼比辨识结果见表1，振型辨识结果如图9所示。

对比基于含噪的时序位移场以及经过空间滤波的位移场所辨识得到的前5 阶固有频率、阻尼比及振型，频率辨识结果在第1 阶偏差较大，达到9.97%。阻尼比受DIC 图像噪声影响敏感，最高达到11.11%。经过空间去噪后所辨识得到的振型更加平滑，符合物理空间的连续性质。

5 有限元计算

5.1 计算模态分析

通过有限元软件ANSYS模态分析模块，进行模型的计算模态分析。有限元模型选用shell 181 单元，单元横截面厚度为1mm，杨氏模量为70 Gpa，泊松比为0.3，密度为2 700 kg/m3，单元大小5 mm×5 mm，模型共划分为4 800 个单元，如图10中(a)所示。

图9 前5阶振型

表1 前5阶模态频率与阻尼比

图10 有限元计算振型

单边固支状态。通过Block Lanczos 法得到结构前5阶固有频率见表2，振型如图10所示。

5.2 对比分析

通过对比有限元仿真与基于DIC 实验的结果，前5 阶固有频率吻合良好，误差小于7.77%。经过对时序位移场去噪重构后提取测点响应，再利用时域信号辨识得到的固有频率误差更小，对于低阶模态表现明显。通过模态置信准则（MAC）判断有限元计算所得振型与OMA 辨识振型的相关性，MAC值均在90%以上，振型结果高度吻合。观察图11，经过空间去噪预处理所辨识得到的振型MAC 值普遍高于未经空间滤波的，其振型表达更加连续平滑。

图11 MAC值比较

6 结语

本文基于数字图像相关方法测量薄板结构振动过程的时序位移场，提取全场的DIC测点时域响应，通过互相关函数代入ERA算法进行工况模态分析，辨识得到系统前5 阶模态频率、阻尼比及振型。通过对比ANSYS仿真与实验结果，有结论如下：

(1)改进了傅里叶连续延拓方法使其适用于二维信号，通过延拓位移场构成周期边界条件，消除空间傅里叶变换中的吉伯斯效应。

表2 前5阶模态频率与阻尼比

(2)引入波数域窗函数，联合傅里叶连续延拓方法，对DIC计算所得含噪位移场结果实现空间滤波；

(3）计算模态分析与基于DIC技术的OMA结果吻合良好，空间去噪方法提高了辨识振型的平滑性。