动态载荷辨识问题的正则化求解新方法*

王　锋， 武　龙， 张小庆， 贺　伟

(1.中国空气动力研究与发展中心超高速空气动力研究所高超声速冲压发动机技术重点实验室， 四川 绵阳 261000；2.中国空气动力研究与发展中心吸气式高超声速技术研究中心， 四川 绵阳 261000)



动态载荷辨识问题的正则化求解新方法*

王锋1， 武龙2， 张小庆2， 贺伟2

(1.中国空气动力研究与发展中心超高速空气动力研究所高超声速冲压发动机技术重点实验室， 四川 绵阳 261000；2.中国空气动力研究与发展中心吸气式高超声速技术研究中心， 四川 绵阳 261000)

Tikhonov正则化方法是求解载荷辨识问题的有效方法，正则化参数的确定是影响求解准确性的关键问题。提出了一种新方法——虚拟载荷法，以虚拟载荷向量的范数最小作为正则化参数确定准则，将载荷辨识问题转化为单参数优化问题，用一维搜索法求得最优解。通过将法方程的系数矩阵预先进行三对角化，建立了高效的求解算法。通过两个数值算例对新方法进行了检验，结果表明新的参数确定准则是有效的，基于该准则得到的载荷辨识结果具有满意的精度。

载荷辨识； 反问题； Tikhonov正则化； 正则化参数

引　言

结构的动态载荷辨识是根据测量的结构响应(位移、速度、加速度或应变等)信息，结合已知的结构动力学特性，反求结构所受激励的过程，属于反问题的一种，有重要的科学研究价值，也有着广泛的工程应用。载荷辨识方法大致可分为频域法和时域法，前者研究较早，方法较成熟，主要通过对关心频段上的系统频响函数矩阵求逆来实现载荷识别[1-2]，比较适用于稳态或平稳随机载荷的识别[3-4]。时域方法直接从结构动力学方程(微分或积分形式)出发，建立整个待辨识时段上载荷输入与结构输出的关系方程，即辨识方程，通过求逆运算，一次得到整个时段上的载荷。然而，辨识方程通常是病态的，解对输出测量噪声十分敏感，几乎不可能用常规的矩阵求(广义)逆的方法直接求解，必须进行适当的正则化处理，基于原来不适定的问题构建一个参数化的适定问题，通过适当选择控制参数来使适定问题的解落在原问题解的一个可接受的邻域内。针对科学和工程中遇到的各种反问题，人们已提出了多种正则化方法[5-7]，如Tikhonov正则化方法、奇异值截断方法、迭代正则化方法等等，其中Tikhonov方法是比较常用的。对每种正则化方法，其最优正则化参数的选取准则都是影响辨识结果的关键因素，是非常重要的研究内容。

本文采用Tikhonov正则化方法求解载荷辨识问题，提出了一种新的准则来选取正则化参数，新准则基于力学概念提出，物理意义明确，而且计算简单。

1　载荷辨识的数学模型

考虑含测量传感器的线性结构动力学系统，其一般形式的动力学方程为

(1)

y=Rx

(2)

此即载荷辨识问题的数学模型。为了保证解的唯一性并有利于提高辨识精度，通常要求结构的测点数不少于待辨识的载荷数目。方程(2)通常是病态的，即使y中仅含有微小的误差，直接求解y=Rx得到的解也会严重偏离真值，从而失去意义。而实际中，测量误差是不可避免的。

2　Tikhonov正则化方法

Tikhonov正则化方法将问题(2)变为参数化的变分问题，用如下最小化问题的解作为其近似解[10-12]

(3)

式中α(>0)是正则化参数。等号右边的第一项使近似解xα尽量满足等式(2)，但又不能使其严格成立，因为y中噪声会使解x严重偏离真解，呈现高频振荡的形态，第二项的作用就是为了抑制x中的振荡，使解尽量光滑，参数α用来在两者之间进行折中、平衡。对式(3)中的目标函数关于x求驻值，可得

(4)

式中I为单位矩阵。由于正则化参数α的取值决定着解的好坏，因此如何选取α是个重要的研究课题。很多研究者从不同的观点出发，已经建立了一些选取准则和算法[7,13]，常用的有L曲线方法、交叉检验方法、偏差准则等，各种方法都有一定的优缺点和适用条件，Bauer和Lukas对此作了比较全面的对比分析[14]。各种方法都是要尽量挖掘原问题可供利用的先验知识或设立某些假设，从而导出确定正则化参数的某种准则。数学上，人们希望正则化解在真解的某个邻域内，但在缺少真解信息的情况下，其邻域也就无从谈起。现有准则通常主要是发掘利用输出的有关信息。针对载荷辨识问题，下面给出一种新的准则，该准则利用解的信息，具有明确的物理意义，而且计算简单。

(5)

该准则还可以按下式的意义来理解

(6)

即将虚拟载荷x2与真实载荷x1的能量之比和x1本身同时稳定在尽可能小的水平。

不过‖x‖2随α是单调减小的[13]，当α→∞时，x(α)→0，从而‖x2‖2→0。于是，可以猜想，若α从0开始增加，‖x2‖2随α的变化趋势应该是先减小(从大的振荡趋向于真值0向量)，后增大(背离真值0向量)，再减小(因α的过正则化)，所以最佳正则化参数应该位于使‖x2‖2从递减转为递增的那个点。

3　求解算法

在求解问题(5)的过程中，需要反复解方程(4)，有必要对其进行预处理，以降低计算量。RTR为实对称矩阵，存在正交矩阵H将其变为对称的三对角矩阵，具体算法可采用Householder变换[15]

HT(RTR)H=Q

(7)

令

(8)

带入式(4)，并用HT左乘式(4)，再利用式(7)并注意H为正交矩阵，可得

(9)

上式中系数矩阵Q+αI为对称三对角矩阵，可以用追赶法快速求解，比直接求解式(4)的计算量大幅下降。

问题(5)是个单参数优化问题，采用一维搜索算法来求解[16]。首先寻找目标函数的单谷区间，然后在此区间上用黄金分割法搜索最优的α。根据前面的讨论，因为当α从0开始增大时，‖x2‖2先递减再递增，所以，在搜索单谷区间时应从一个充分小的α开始，以避免直接进入后面的单调递减区间，导致搜索失败。此外，因为在最优的正则化参数附近，正则解对参数的变化十分敏感，也就是说‖x2‖2随α的变化曲线在最优解附近可能十分陡峭且狭窄，为了易于捕捉到包含最优解的单谷区间，应该在对数尺度下进行搜索，即令α=10β，将β作为搜索变量，这样，‖x2‖2随β的变化曲线在最优解附近会变得“平坦”和“宽阔”。

因为R=[R1R2]，将方程(4)分解

(10)

而原问题对应的方程是

(11)

当x2不为零时，由式(10)解出的x1中包含了x2产生的附加干扰，为了得到更好的解，在确定了最优的α之后，由式(11)来求出x1作为最终解。

综上所述，问题的求解过程如下：

(1)将RTR变换为三对角矩阵，得到H和Q，计算式(9)的右端项；

(2)令α=10β，以β为自变量，搜索‖x2‖2随β变化曲线的单谷区间，β的初始值为一个适当的负数，在搜索过程中需针对不同的β求解方程(9)，再根据式(8)得到目标函数值‖x2‖2；

(3)确定单谷区间[β1，β2]后，用黄金分割法在该区间上搜索最优的β，得到对应的正则化参数α；

在今年年初的时候，大部分人都无法准确地预料得到今年的行情会如此之“惨”。如果是从春节作为行情划分的调整时点算起，那么截止到2018年末，市场整整调整了11个月，这在A股历史上都是较为“罕见”的。所以，2018年，带给投资者的几乎都是损失，只是每个人亏多亏少而已。所有的公募股票型主动管理基金，都是以亏损告终。那只在今年6月30日才由原来保本型货币基金转型而来的股票型基金，因为一直空仓到年末，反而拿了公募股票型基金的冠军。这足以说明A股市场2018年的惨淡。

(4)由方程(11)解出原问题的解x1。

4　算例分析

4.1平面桁架载荷辨识

考虑一平面桁架结构，如图1所示，由79根长0.5 m、截面直径15 mm的钢质杆构成，在两端简支。假设在第31节点沿y方向作用一激励载荷，如下

f(t)=100(1-cos2πω0t)sin6πω0t

(12)

其中，ω0=10 Hz。先积分求解系统的位移响应，然后对测点的位移计算结果叠加一定强度的高斯白噪声,模拟实际测量结果，噪声的均方差按下式确定

(13)

式中yj表示第j个测量输出向量，l是其长度，ny是选取的测点个数，Δe是比例系数。

图1　平面桁架结构Fig.1　Planar truss structure

在此，取第5和16节点的y向位移输出来对f(t)进行辨识，而虚拟载荷假设沿y方向作用于第6个节点上，噪声比例系数Δe取0.01。输出采样步长取0.001 s，时长0.5 s。载荷辨识结果与真实值的对比如图2所示，二者吻合良好。图3是计算得到的虚拟载荷历程，其幅值相对于真实载荷较小，而且看上去类似零均值的噪声，这是所期望的，说明真实载荷的信息没有“泄露”到虚拟载荷上而影响到其辨识精度。

图2　载荷辨识结果与真实值的对比Fig.2　Comparison between true load and identified load

图3　虚拟载荷辨识结果Fig.3　Identified fictitious load

作为对比，不加虚拟载荷，而基于L曲线方法确定正则化参数，用Hanson开发的工具箱[17]来求解，得到最佳参数为6.8924×10-16，图4中标出了其位置，比本文方法确定的参数值约小两个量级，在该参数下得到的结果与真实载荷的对比如图5所示，对噪声的滤除效果稍差，辨识误差比本文方法要大。

图4　目标函数随正则化参数的变化Fig.4　Curve of objective function vs. regularization parameter

图5　基于L曲线方法的载荷辨识结果Fig.5　Load identification results based on L-Curve method

可以想象，取不同的虚拟载荷作用点，得到的辨识结果应该会有所不同。极端的情况是，虚拟载荷施加在真实载荷的作用点上并且方向相同时，两者的作用效果相同，应该无法辨识出真实载荷。因此，虚拟载荷的施加点应避免与真实载荷作用点距离太近。为考察虚拟载荷作用点的影响，定义辨识误差如下

(14)

图6　虚拟载荷作用于不同节点时的辨识误差Fig.6　Relative error history vs. action point of fictitious load

可见，当选择距离实际载荷作用点较近的节点时，辨识误差较大，如第30，32节点处约6%，第10，11节点处约9%，而选择实际载荷作用点第31节点时，误差为100%，辨识失败。而其余节点对应的辨识误差都差别不大，均在3.5%左右，这说明方法对虚拟载荷作用点的选择并不敏感。即使将作用点选在较差的第10节点，辨识结果也是可以接受的，其与真实载荷的对比如图7所示，甚至优于基于L曲线方法的结果。

图7　虚拟载荷作用点选在节点10时的辨识结果Fig.7　Identification results when the 10th node is selected as action point of the fictitious load

4.2四边固支板上的载荷辨识

考虑一矩形钢板，长0.8 m，宽0.6 m，厚度2.5 mm，四边固支。用有限元方法建模，单元划分如图8所示。假设板上同时作用2个载荷，选取4个节点的挠度测量值来对其进行辨识，使用1个虚拟载荷。图8中，●表示真实载荷作用点，○表示虚拟载荷作用点，⊗为位移输出测点。

图8　矩形板网格及载荷作用点和测点分布Fig.8　Meshes of the plate, action points and measuring points

2个激励载荷分别取间断的三角形脉冲和连续的正弦信号，作用时长0.4 s。4个测点的位移响应如图9所示，采样步长取0.001 s。可见其主要由正弦激励主导，两次脉冲激励仅引起相对较小的扰动。

图9　4个测点的位移输出Fig.9　History of outputs from the four measurement points

图10　载荷辨识结果Fig.10　Load identification result

虚拟测量噪声强度的定义仍为式(13)，当取Δe等于0.01时，2个载荷的辨识结果分别如图10所示，与真实载荷吻合很好。此时，目标函数‖x2‖2随α的变化曲线如图11所示，变化规律依然符合预期。而采用L曲线准则确定的最佳正则化参数为8.06×10-15，标于图11中，比本文方法确定的参数约小两个数量级。基于该参数得到的辨识结果如图12所示，与图10相比，精度要差一些。

图11　目标函数随正则化参数的变化Fig.11　Curve of objective function vs. regularization parameter

图12　基于L曲线方法的载荷辨识结果Fig.12　Load identification results based on L-Curve method

现考察噪声强度对辨识结果的影响。取不同大小的噪声，辨识结果的误差和对应的最优正则化参数列于表1，随着噪声的增大，辨识误差随之增大，最优正则化参数也在增大，在对数标度下均与噪声大小近似成线性关系，如图13所示。

表1　噪声水平对辨识结果的影响

图13　辨识误差和最优正则化参数随噪声水平的变化Fig.13　Variations of relative error and the optimal regularization parameter vs. level of noise

5　结　论

本文采用Tikhonov正则化方法求解线性结构的载荷辨识问题，通过增加一个虚拟载荷，使增广问题的解部分已知，从而给出确定最优正则化参数的准则，该准则物理意义明确，计算简单。在此基础上给出了优化求解算法，算例表明，本文方法能够稳健地辨识结构所受动态载荷，具有一定的应用价值。

不过，目前的工作仍属于初步探索，尚需开展进一步的理论研究，为新方法的合理性提供支撑，至少包括两方的深入工作：(1)从数学上证明虚拟载荷范数随正则化参数变化曲线极小点的存在性和唯一性，或者给出极小点存在的条件，即什么情况下本文提出的准则是适用的；(2)建立虚拟载荷作用点优劣的判定函数，用于选择最佳的虚拟载荷作用点。

本文的工作也为其他反问题求解提供了新的思路，即，设法将问题进行增广，注入已知信息，然后在求解过程中充分利用已知信息来约束增广问题的解，使其逼近真解。

[1]胡杰,张希农.一种频域载荷识别的优化方法[J].噪声与振动控制, 2009,(6):34—36.

HU Jie, ZHANG Xi-nong. An optimization method of load identification in frequency domain[J]. Noise and Vibration Control, 2009,(6):34—36.

[2]智浩,文祥荣,缪龙秀,等.动态载荷的频域识别方法[J].北方交通大学学报,2000,24(4):5—10.

ZHI Hao, WEN Xiang-rong, MIAO Long-xiu, et al. Dynamic loading identification in frequency domain[J]. Jounal of North Jiaotong University, 2000,24(4):5—10.

[3]周林,郑四发,王彬星,等.动态载荷识别位置优化的传递函数相干法[J].振动工程学报,2011,24(1):14—19.

ZHOU Lin, ZHENG Si-fa, WANG Bin-xing, et al. Coherence analysis method for dynamic force identif ication[J]. Journal of Vibration Engineering, 2011,24(1):14—19.

[4]朱涛,肖守讷,阳光武.载荷识别研究进展及其运用于铁道轮-轨载荷研究概述[J].铁道学报,2011,33(10):29—36.

ZHU Tao, XIAO Shou-ne, YANG Guang-wu. State-of-the-art development of load identification and it′s application in study on wheel-rail forces[J]. Journal of the China Railway Society, 2011,33(10):29—36.

[5]Sun X, Liu J, Han X, et al. A new improved regularization method for dynamic load identification[J]. Inverse Problems in Science and Engineering, 2014,22(7):1062—1076.

[6]Wang L, Han X, Liu J, et al. An improved iteration regularization method and application to reconstruction of dynamic loads on a plate[J]. Journal of Computational and Applied Mathematics, 2011,235:4083—4094.

[7]Vogel C R. Computational Methods for Inverse Problems[M]. Philadelphia: SIAM, 2002.

[8]郭杏林,毛玉明,赵岩,等.基于Markov参数精细积分法的载荷识别研究[J].振动与冲击,2009,28(3):27—31.

GUO Xing-lin, MAO Yu-ming, ZHAO Yan, et al. Load identification based on precise time-step integration for Markov parameters[J]. Journal of Vibration and Shock, 2009,28(3):27—31.

[9]韩旭,刘杰,李伟杰,等.时域内多源动态载荷的一种计算反求技术[J].力学学报,2009,41(4):595—602.

HAN Xu, LIU Jie, LI Wei-jie, et al. A computational inverse technique for reconstruction of multisource loads in time domain[J]. Chinese Journal of Theoretical and Applied Mechanics, 2009,41(4):595—602.

[10]Lampe J, Reichel L, Voss H. Large-scale Tikhonov regularization via reduction by orthogonal projection[J]. Linear Algebra and its Applications, 2012,436(8):2845—2865.

[11]Reichel L, Rodriguez G. Old and new parameter choice rules for discrete ill-posed problems[J]. Numerical Algorithms, 2012:1—23.

[12]梅立泉,崔维庚.面载荷识别的TSVD正则化方法[J].应用力学学报,2010,27(1):140—146.

MEI Li-quan, CUI Wei-geng. TSVD regularization method for area load reconstruction[J]. Chinese Journal of Applied Mechanics, 2010,27(1):140—146.

[13]Hansen P C, O Leary D. The use of L-curve in the regularization of discrete ill-posed problems[J]. SIAM Journal on Scientific Computing,1993,14:1487—1503.

[14]Bauer F, Lukas M A. Comparing parameter choice methods for regularization of ill-posed problems[J]. Mathematics and Computers in Simulation, 2011,81(9):1795—1841.

[15]张贤达.矩阵分析与应用[M].北京:清华大学出版社,2004:248—250.

ZAHNG Xian-da. Matrix Analysis and Applications[M]. Beijing: Tsinghua University Press, 2004:248—250.

[16]袁亚湘,孙文瑜.最优化理论与方法[M].北京:科学出版社,1997:56—75.

YUAN Ya-xiang, SUN Wen-yu. Optimization Theory and Method[M]. Beijing: Science Press, 1997:56—75.

[17]Hansen P C. Regularization tools version 4.0 for Matlab 7.3[J]. Numerical Algorithms, 2007,46:189—194.

A technique for solving dynamical force identification problems by Tikhonov regularization method

WANGFeng1,WULong2,ZHANGXiao-qing2,HEWei2

(1.Science and Technology on Scramjet Laboratory, Hypervelocity Aerodynamics Institute of CARDC, Mianyang 621000, China；2.Air-breathing Hypersonic Technology Research Center of CARDC, Mianyang 621000, China)

Tikhonov regularization method is effective in solving load identification problems, and proper selection of the regularization parameter is a key point that determines the solution accuracy. In this work, a new method which is named fictitious load method is proposed. Taking the minimization of the norm of the fictitious load is as the criterion of regularization parameter selection, the load identification problem isthus converted into a one-variable optimization problem which can be solved by one dimension search technique. By transforming the coefficient matrix of the normal equation into tridiagonal matrix, an efficient algorithm is developed. Two numerical examples are presented to validate the new method, and the results show that the new regularization parameter selection criterion is effective to revover the loads on structures with good accuracy.

load identification; inverse problem; Tikhonov regularization; regularization parameter

2014-06-24；

2015-11-16

国家自然科学基金资助项目(11372339)；高超声速冲压发动机技术重点实验室基金资助项目(STSKFKT2012001)

O327

A

1004-4523(2016)01-0031-07

10.16385/j.cnki.issn.1004-4523.2016.01.005

王锋(1976—),男，博士，副研究员。电话： 15908208358； E-mail： fwang_cardc@163.com