范玉川¹, 陈晔², 赵春雨¹, 卢泽宸¹
1. 东北大学 机械工程与自动化学院, 辽宁 沈阳 110819;
2. 辽宁工业大学 机械工程与自动化学院, 辽宁 锦州 121001
收稿日期：2018-12-24
基金项目：国家自然科学基金资助项目(51775094)。
作者简介：范玉川(1988-)，男，河南新乡人，东北大学博士研究生;
赵春雨(1963-), 男, 辽宁黑山人，东北大学教授，博士生导师。

摘要：传统的L曲线法在使用的时候常常不容易获得准确的正则化参数, 基于此, 提出了一种基于Newmark-β的反算-对比-调整-逼近(inverse computation-contrast-adjustment-approach, ICAA)正则化参数选取方法.该算法相比传统的L曲线法使用起来更加直观、简便, 并且计算耗费的时间更短、效率更高.通过一个四自由度系统的仿真算例和一个悬臂梁的实验验证了本算法的有效性, 并把本算法的载荷识别结果与L曲线法的载荷识别结果进行了对比.结果表明:该算法相比L曲线法不仅在计算效率方面有显著优势, 而且利用前者的正则化参数进行载荷识别，计算精度更高.
关键词：正则化参数L曲线法载荷识别Newmark-β计算效率
Regularization Parameter Determination Method of ICAA Based on Newmark-β
FAN Yu-chuan¹, CHEN Ye², ZHAO Chun-yu¹, LU Ze-chen¹
1. School of Mechanical Engineering & Automation, Northeastern University, Shenyang 110819, China;
2. School of Mechanical Engineering and Automation, Liaoning University of Technology, Jinzhou 121001, China
Corresponding author: ZHAO Chun-yu, E-mail: chyzhao@mail.neu.edu.cn
Abstract: When the traditional L-curve method is applied, it is often not easy to obtain accurate regularization parameters. Therefore, a new regularization parameter determination method of ICAA(inverse computation-contrast-adjustment-approach)based on Newmark-β was proposed. Compared with the traditional L-curve method, the algorithm is more intuitive and simple to use, and the calculation time is shorter and more efficient. The validity of the algorithm was verified by a simulation example of a four-DOF system and an experiment of a cantilever beam, and the load identification results of this algorithm were compared with the load recognition results of the L-curve method. The results showed that compared with the L-curve method, the proposed algorithm not only has greater advantages in computational efficiency, but also has higher accuracy in load identification calculation using the regularization parameters obtained by the algorithm.
Key words: regularization parameterL-curve methodload recognition identificationNewmark-βcomputational efficiency
动载荷识别方法主要分为频域法^[1-2]和时域法^[3]两大类.由于时域法能够识别各种类型的载荷, 识别精度高, 并且其识别结果具有明确的物理意义, 因而, 时域法越来越受到专家学者们的青睐.
动载荷识别问题属于第二类反问题^[4-5], 在反问题的研究中常常会遇到“病态”问题^[6], 正则化方法是解决病态问题的有效手段^[7-8], 传统的正则化方法有Tikhonov正则化、截断奇异值分解、小波正则化方法等.
在正则化过程中, 正则化参数的选取至关重要, 如果λ值过大就会使识别误差增大, λ的值过小则起不到正则化的作用.目前, 最常用的正则化参数选取方法是L曲线法^[9], 而L曲线法常常在使用过程中很难找到“拐点”, 这就给它的应用造成了很大的不便, 此时运用L曲线法往往很难找到准确的正则化参数值.
为了解决正则化参数值难以获得的问题, 本文提出了一种基于Newmark-β的反算-对比-调整-逼近(ICAA)正则化参数选取方法, 并通过一个仿真算例和实验对其有效性进行验证.
1 动载荷识别模型由文献[10]可得:

(1)

式中:Y=[y(t₁)^T y(t₂)^T … y(t_nt)^T]^T为总体响应向量; P=[P(t₁)^T P(t₂)^T … P(t_nt)^T]^T为总体载荷向量; y(t_i)和P(t_i)分别表示第i个时刻的响应向量和载荷向量, i=1, 2, …, nt; y(t_i)的维度等于系统的自由度数m×3;P(t_i)的维度等于系统的自由度数m;
为核函数矩阵; L为映射矩阵.
考虑到式(1)可能存在不适定性, 可以通过Tikhonov正则化方法来确定载荷P, 在阻尼最小二乘意义下对目标函数进行优化:

(2)

式中, λ是正则化参数.工程中通常用L曲线法来确定正则化参数λ.
2 基于Newmark-β的ICAA正则化参数选取方法由于L曲线法在实际应用中存在使用繁琐, 有时难以找到最佳正则化参数值, 以及计算时间长的缺点.基于此, 本文提出了一种基于Newmark-β的ICAA正则化参数选取方法.
2.1 Newmark-β仿真算法Newmark-β法作为逐步积分法之一, 其原理是假设加速度在一个时间步内按照线性规律变化, 即有以下速度和位移关系：

(3)

(4)

式中:x_t, 和分别表示第t个时刻的位移向量、速度向量和加速度向量；Δt表示时间间隔；γ和β是可调参数, 取值不同, 就形成不同的逐步积分法.当γ=1/2, β=1/4得到平均加速度法; 当γ=1/2, β=0得到中心差分法; 当γ=1/2, β=1/6得到线性加速度法.
在t+Δt时刻, 有动力方程:

(5)

由式(4)可将用x_t+Δt, 和表示：

(6)

把式(6)代入式(3), 可将用x_t+Δt, 和表示：

(7)

把式(6)和式(7)代入式(5), 可得关于x_t+Δt的方程

(8)

式中:
由式(8)可解得x_t+Δt,再进一步由式(6)和式(7)求解得到和.
2.2 反算-对比-调整-逼近(ICAA)正则化参数选取方法以上是Newmark-β仿真算法的描述, 下面具体说明基于Newmark-β的反算-对比-调整-逼近(ICAA)的正则化参数选取方法的计算过程.
反算:假设任取一较大的正则化参数λ值作为初始值, 将其代入动载荷识别算法及正则化方法中进行载荷识别计算得到未知载荷, 然后代入Newmark-β仿真算法计算得到响应, 和.
对比:将计算得到的响应信息(假设为位移信息)与测量得到的响应信息(假设为位移信息)利用式(9)计算其相对积累误差, 这里可以取ε作为阈值, 当

(9)

所得到的正则化参数λ的值即为最佳正则化参数, 此时计算得到的载荷值即为较为准确的动载荷识别结果.
调整:如果在对比过程中未达到所设定的阈值, 则需要调节λ的值, 直到式(9)达到所设定的阈值为止.
逼近:逼近过程和调整过程是同时进行的, 逼近指的是将正则化参数λ的值往满足式(9)的方向调整.实际应用中, 可以先取一个较大的正则化参数初值(0＜λ＜1, 可取0.01), 每次递减一个数量级, 往较小的正则化参数值调整(例如:第一次取10^-2, 第二次取10^-3, 第三次取10^-4, 依次类推, 当相对积累误差接近阈值ε时, 可细化调整, 取5×10^-5, 4×10^-5, 3×10^-5, 依次类推, 直到满足式(9), 终止计算).
该方法可以给计算机一个搜索方向让计算机自动搜索最佳正则化参数λ, 同时由于本方法是采用位移响应信息对比的方法, 操作起来比较简便直观.该方法也可以采用人工对比调节的方法, 便于根据反算结果进行对比, 更灵活地对正则化参数值进行实时的调整, 这样有助于节省大量的计算时间, 因而更能适应工程实际的工作环境.
基于Newmark-β的ICAA正则化参数选取方法如图 1所示.
图 1(Fig. 1)

图 1 ICAA正则化参数选取方法流程图Fig.1 Flow chart for regularization parameter determination method of ICAA

3 仿真算例本节通过一个四自由度的算例详细描述基于Newmark-β的ICAA正则化参数选取方法的具体计算流程, 并把该算法的计算结果与L曲线法的计算结果进行对比.
四自由度系统如图 2所示, 各自由度的质量分别为m₁=4 258 kg, m₂=2 258 kg, m₃=2 437 kg, m₄=1 729 kg.刚度k=8.64×106 N/m, 阻尼形式为比例阻尼, 阻尼系数为α₁=0.1, α₂=3.615 7×10^-4.算例中, 采样频率设置为1 000 Hz, 采样时间设置为3 s.
图 2(Fig. 2)

图 2 四自由度结构图Fig.2 Four degrees of freedom structure diagram

算例中识别结果和真实载荷的误差通过式(10)和式(11)进行计算：

(10)

(11)

式中:e₁表示相对误差; e₂表示相对积累误差; F_id表示识别的载荷; F_real表示真实载荷; |·|表示绝对值; ‖·‖表示l₂-范数.
在该四自由度系统的第三个自由度上施加谐波载荷:

(12)

利用第二个节点上的位移响应数据对第三个节点上的载荷进行识别计算.由于实际生产中不可避免地会受到测量噪声的影响, 所以在本算例中, 根据式(13)在所得到的位移信号中添加3%的随机噪声信号.

(13)

式中:l_noise表示噪声水平; std(·)是MATLAB脚本函数, 表示向量的标准偏差; MATLAB脚本函数rand(n, 1)可以得到一个n×1的随机向量, 表示区间(0, 1)上的服从标准均匀分布中的伪随机值.
利用基于Newmark-β的ICAA正则化参数选取方法对正则化参数进行确定, 本算例中阈值ε设置为3%, 每次迭代计算得到的正则化参数值、反算得到的位移相对积累误差以及载荷的相对积累误差见表 1.
表 1(Table 1)

表 1 基于Newmark-β的ICAA正则化参数Table 1 ICAA regularization parameters based on Newmark-β 选择次数 正则化参数 位移相对积累误差/% 载荷相对积累误差/% 1 1e-10 99.70 99.96 2 1e-12 64.08 89.62 3 1e-14 10.67 16.42 4 1e-15 3.86 6.30 5 1e-16 8.66 19.35 6 2e-16 6.88 10.01 7 3e-16 5.51 7.03 8 4e-16 2.76 4.68 9 终止计算

表 1 基于Newmark-β的ICAA正则化参数 Table 1 ICAA regularization parameters based on Newmark-β

由表 1可知, 当正则化参数值为4e-16时, 位移的相对积累误差达到2.76%, 满足了阈值ε, 终止计算.此时, 载荷识别相对积累误差为4.68%.利用相同的数据, 用L曲线法进行相同的计算, 此时得到的最佳正则化参数为8e-16, 对应的载荷相对积累误差为6.80%, 相比前者载荷识别的结果, 相对积累误差大了不少.
为了更直观地对比两种算法的计算结果, 两种算法计算得到的载荷识别结果如图 3所示, 两种算法计算得到的载荷识别相对误差曲线如图 4所示.由图 3, 4可知, 利用两种正则化参数选取方法进行载荷识别, 识别得到的载荷与真实载荷吻合较好, 识别误差比较小, 但本文算法在识别准确性上略有优势.
图 3(Fig. 3)

图 3 利用第2个节点上的位移响应对第3个节点上的激励的识别结果Fig.3 Recognition results of excitation on the 3rd node by using the displacement response on the 2nd node

图 4(Fig. 4)

图 4 利用第2个节点上的位移响应对第3个节点上的激励的识别误差Fig.4 Recognition errors of the excitation on the 3rd node by using the displacement response on the 2nd node

为了进一步比较两算法的优劣, 在不同噪声水平下利用两种正则化参数选取方法进行载荷识别计算得到的识别误差和所耗费的时间见表 2.
表 2(Table 2)

表 2 不同噪声水平下两种算法的载荷识别结果Table 2 Load recognition results of two algorithms at different noise levels 噪声水平/% 正则化参数 相对积累误差/% 耗费时间/s ICAA L曲线法 ICAA L曲线法 ICAA L曲线法 0 7e-16 5e-17 0.79 0.88 52.4 566.8 10 8e-12 2e-13 2.66 2.78 44.8 571.5 20 3e-8 6e-8 4.21 4.87 60.6 544.8 30 4e-6 8e-6 5.10 6.80 37.6 577.3

表 2 不同噪声水平下两种算法的载荷识别结果 Table 2 Load recognition results of two algorithms at different noise levels

由表 2可知:在相对积累误差方面, 两种算法都有较好表现, 本文算法有略微的优势; 但在所耗费的时间方面, 本文算法有较大优势, 这也从数据上体现出了本文算法在正则化参数选择效率方面的高效性.在动载荷识别计算中, 利用本文算法能够快速直观地得到较为准确的正则化参数, 节省了大量的时间, 提高了动载荷识别计算的效率.
为了验证本算法在多激励多输入状态下的效果, 分别在m₂和m₄上施加载荷f₂和f₄, 且f₂和f₄分别为

(14)

(15)

利用第1和第3个节点上的位移进行载荷识别的计算, 计算前, 同样在响应数据中加入3%的噪声.
利用基于Newmark-β的ICAA正则化参数选取方法得到的正则化参数为8×10^-16, 计算得到的f₂的相对积累误差为16.5%, f₄的相对积累误差为11.2%, 所耗费的计算时间为86.5s;利用L曲线法得到的正则化参数为2×10^-16, 计算得到的f₂的相对积累误差为18.3%, f₄的相对积累误差为26.2%, 所耗费的计算时间为773.4 s.可知, 在多源载荷下基于Newmark-β的ICAA正则化参数选取方法得到的正则化参数在准确性和计算效率方面仍然优于L曲线法.
通过两种算法得到的f₂和f₄载荷识别结果如图 5和图 6所示, 可知两种算法正则化参数的准确性的直观对比.
图 5(Fig. 5)

图 5 载荷f2的识别结果对比Fig.5 Comparison of identification results of load f2

图 6(Fig. 6)

图 6 载荷f4的识别结果对比Fig.6 Comparison of identification results of load f4

4 实验验证本节通过一个悬臂梁实验台, 对悬臂梁施加谐波载荷, 采集位移响应信号, 分别对两种算法进行载荷识别的计算并对载荷识别结果进行对比.
实验所用的数据采集设备为National Instruments公司的NI 9234数据采集卡, 所用的数据采集软件是NI公司的著名软件LABWIEW的SignalExpress模块.
悬臂梁实验台如图 7所示, 悬臂梁的结构图如图 8所示, 悬臂梁被划分成了18个单元19个节点进行处理, 每个节点仅考虑径向和弯曲两个自由度, 共有38个自由度, 其中, 节点1上的两个自由度受到全约束的限制.悬臂梁的前四阶固有频率分别为15.9, 99.8, 279.5和547.7 Hz.
图 7(Fig. 7)

图 7 悬臂梁试验台Fig.7 Cantilever beam test bench

图 8(Fig. 8)

图 8 悬臂梁结构图Fig.8 Cantilever beam structure

该悬臂梁的弹性模量为200 GPa, 密度为7 840 kg/m³, 泊松比为0.33, 长、宽、高分别为0.64, 0.056和0.008 m.
在图 8的第10个节点处施加径向的谐波载荷, 利用位移传感器测量得到第15个节点上的位移响应, 利用文献[10]中的载荷识别方法进行计算, 计算中分别利用基于Newmark-β的ICAA正则化参数选取方法和L曲线法获得正则化参数, 前者得到的正则化参数为7e-11, 载荷识别得到的相对积累误差为2.34%, 后者得到的正则化参数为2e-11, 载荷识别计算得到的相对积累误差为4.57%.
载荷识别结果如图 9所示, 图中可以清晰地看到两种算法在实验中的载荷识别结果, 结果显示两种算法在悬臂梁实验中效果都较好, 相比之下, 本文算法得到的正则化参数的载荷识别误差更小, 精度更高.
图 9(Fig. 9)

图 9 悬臂梁载荷识别结果Fig.9 Load identification results of cantilever beam

5 结论1) 本文提出了一种基于Newmark-β仿真算法的ICAA正则化参数选取方法, 该算法相比传统的L曲线法使用起来更加直观、简便且计算效率更高.
2) 通过一个四自由度系统的仿真算例和一个悬臂梁实验, 验证了本文算法在正则化参数选取方面的有效性, 并把载荷识别结果与传统的L曲线法的载荷识别结果进行了对比.结果表明:本文算法相比L曲线法不仅在计算效率方面有显著的优势, 而且利用前者取得的正则化参数进行载荷识别计算得到的结果精度更高.
参考文献

[1]	Hu J, Zhang X. An optimization method of load identification in frequency domain[J].Noise & Vibration Control, 2009, 29(6): 34–36.
[2]	He Z C, Lin X Y, Li E. A novel method for load bounds identification for uncertain structures in frequency domain[J].International Journal of Computational Methods, 2017(3): 1850051.
[3]	Law S S, Chan T H T, Zhu Q X. Regularization in moving force identification[J].Journal of Engineering Mechanics ASCE, 2001, 127(2): 136–148.DOI:10.1061/(ASCE)0733-9399(2001)127:2(136)
[4]	Uhl T. The inverse identification problem and its technical application[J].Archive of Applied Mechanics(Ingenieur Archiv), 2007, 77(5): 325–337.
[5]	Khan A A, Rouhani B D. Iterative regularization for elliptic inverse problems[J].Computers and Mathematics with Applications, 2007, 54(6): 850–860.DOI:10.1016/j.camwa.2007.02.011
[6]	周盼, 张权, 率志君, 等. 动载荷识别时域方法的研究现状与发展趋势[J].噪声与振动控制, 2014, 34(1): 6–11. ( Zhou Pan, Zhang Quan, Shuai Zhi-jun, et al. Review of research and development status of dynamic load identification in time domain[J].Noise and Vibration Control, 2014, 34(1): 6–11.DOI:10.3969/j.issn.1006-1335.2014.01.002)
[7]	Thite A N, Thompson D J. The quantification of structure-borne transmission paths by inverse methods.Part 1:improved singular value rejection methods[J].Journal of Sound & Vibration, 2003, 264(2): 411–431.
[8]	刘继军. 不适定问题的正则化方法及应用[M]. 北京: 科学出版社, 2005. ( Liu Ji-jun. Regularization method and application of ill-posed problems[M]. Beijing: Science Press, 2005.)
[9]	Leary O, Prost D, Hansen P C. The use of the L-curve in the regularization of discrete ill-posed problems[J].SIAM Journal on Scientific Computing, 1993, 14(6): 1487–1503.DOI:10.1137/0914086
[10]	Liu K, Law S S, Zhu X Q, et al. Explicit form of an implycit method for inverse force identification[J].Journal of Sound and Vibration, 2014, 33: 730–744.