高 伟， 于开平， 盖晓男

(哈尔滨工业大学 航天学院，哈尔滨 150001)

L∞范数拟合正则化方法在飞行器动态载荷识别中的应用

高 伟， 于开平， 盖晓男

(哈尔滨工业大学 航天学院，哈尔滨 150001)

基于变限积分理论，构造了响应函数最小二乘意义下的加权变限积分。通过适当次数的积分滑动平均，有效过滤测量噪声中的高频噪声。针对测量响应中残留的低频噪声，使用L∞范数拟合正则化方法识别载荷，提出了一种选取L∞范数拟合正则化方法最优正则化参数的单调性检验方法。数值仿真及试验验证说明单调性检验方法可以有效确定L∞范数拟合正则化方法的最优正则化参数，得到比传统L2范数正则化方法拟合性质更好精度更高的识别载荷；针对遥测数据中噪声特点，使用模拟遥测数据利用L∞范数拟合正则化方法对冲击载荷进行了有效识别。

载荷识别；L∞范数正则化方法；滑动平均；遥测数据

载荷识别问题和系统辨识问题是反问题理论中的两类问题[1-2]。很多实际工程问题中，直接测量工程结构所受动态载荷非常困难甚至是不可能的，比如火箭点火和分离时飞行器结构所受冲击载荷、导弹在发射和飞行时某些关键部位的动态载荷信息等。飞行器关键部位所受动态载荷信息在飞行器模态分析、内力识别、强度校核等领域扮演关键角色。所以载荷识别技术越来越受到人们的重视，载荷识别技术是在结构响应及系统模型已知的情况下识别外部动态载荷[3-5]。载荷识别问题往往是不适定的，伴随正则化方法[6]的提出大量专家投入到载荷识别技术的研究中来。

目前，载荷识别技术主要有频域和时域两类[7]。载荷识别频域方法发展较早，基本思想是在频域内建立输入、系统、输出三者之间的系统模型。频域方法载荷识别系统模型的建立及求解相对容易，对于周期形式的载荷可以得到较好的识别结果。但是频域方法也存在自身的局限性，比如在共振频率附近频响函数矩阵呈现病态。另外，频域方法还要求结构响应数据具有一定的长度以便得到高精度的频谱，所以频域方法对于瞬态脉冲载荷的识别具有较大的局限性。与频域方法相比较，载荷识别时域方法对数据长度并无要求，所以时域方法更适用于瞬态脉冲载荷识别问题。Doyle[8]提出的小波反卷积方法识别的冲击载荷比频域方法识别结果有更高的精度。载荷识别频域方法得到的结果为关心频带内的载荷谱。而时域方法得到的结果是更加直观的时域内的载荷时间历程，所以载荷识别时域方法更容易研究载荷在时域内的表现形式。载荷识别时域方法中，正则化方法是目前使用最普遍的载荷识别方法。Sanchez等[9]总结了一些目前得到充分研究及成功用于载荷识别的L2范数正则化方法，指出系统矩阵的病态性及测量噪声是引起载荷识别误差的主要因素。L2范数是高斯白噪声的最佳最小二乘近似，所以传统L2范数正则化方法针对高斯白噪声可以得到理想的载荷识别结果[10]。但是在实际工程应用中存在高斯白噪声水平较高的情况，积分滑动平均可以有效抑制高水平高斯白噪声中的高频噪声部分，残留的低水平噪声部分不再是白噪声，而是较小幅值的低频形式的噪声。另外，在飞行器关键部位动态载荷识别所使用的遥测数据采集过程中，由于环境噪声及数据采集系统噪声的影响，遥测数据中的噪声不具有高斯白噪声特性，而是近似服从均匀分布的随机噪声。前面提到的较小幅值的低频形式噪声及遥测数据中近似服从均匀分布的随机噪声，不是L2范数正则化方法的最佳适用噪声形式。L∞范数是均匀噪声的最佳最小二乘近似，Clason[11]提出的L∞范数拟合正则化方法针对均匀噪声可以得到理想的载荷识别结果，所以本文引入L∞范数拟合正则化方法进行飞行器动态载荷识别并与传统L2范数正则化方法载荷识别结果相比较。

本文主要研究L∞范数拟合正则化方法在实际工程载荷识别问题中的应用。下文组织如下。第1部分建立用于载荷识别的离散系统方程积分滑动平均模型。第2部分介绍了L∞范数拟合正则化方法；第3部分基于最优化问题的目标函数，提出了用于确定L∞范数拟合正则化方法最优正则化参数的单调性检验方法。第4部分通过数值仿真算例及试验验证说明了L∞范数拟合正则化方法及单调性检验方法的合理性及有效性，并且给出了与L2范数正则化方法的比较结果。第5部分中给出使用新方法基于模拟遥测数据对加载在飞行器结构上关键部位冲击载荷的无量纲化识别结果。最后在第6部分中给出结论。

1 离散系统方程模型

载荷识别离散系统方程模型的建立方法之一是基于Duhamel积分的离散化。

(1)

(2)

式中，n=1,2,…为积分次数。适当进行多次积分以后可以得到方程为

(3)

(4)

再令

(5)

则由(5)式中矩阵可逆有

(6)

由式(6)可得权函数为

hk(t)=u[t+(m-k)Δt]

(7)

式中：4≤k≤Q-2；(k-1)Δt

(8)

同理有

(9)

再利用最小二乘拟合形函数方法[13]，基于前面整个时间域的离散化方式可得离散系统方程为yn=Gp，再由5.2次平滑系数令矩阵M为

(10)

由式(9)及式(10)有

yn=Myn-1=…=Mny*

(11)

式中，y*=[y*(t1),y*(t2),…,y*(tQ)]T。取定合适的积分次数n=n0，则有离散系统方程模型

Mn0y*=Gp

(12)

2 L∞范数拟合正则化方法

针对式(12)中离散系统方程，L∞范数正则化方法考察最优化问题

(13)

(14)

式中，c∈R为一实数。约束条件中L∞范数是不可导的，所以利用Moreau-Yosida近似来解决这个问题。对于γ>0，式(14)的Moreau-Yosida近似为

(15)

式中，max,min两个函数均为逐点运算。式(15)中Moreau-Yosida近似所考察最优化问题的解强收敛到式(14)中所考察的最优化问题的解。式 (15)中最优化问题目标函数为严格凸及弱下半连续。由式(15)中最优化问题的目标函数可导，可以使用半光滑牛顿法求解式(15)中最优化问题[15-16]。

3 最优正则化参数选取方法

由式(13)中最优化问题目标函数，我们定义函数为

(16)

(2) 当α*≈α时，式(13)中拟合项与正则项处于较好的折中状态，所以在参数α*附近非常小的范围内，函数H(α)的值几乎不变。

(4) 所以我们选定最优正则化参数值αop满足H(αop)=minH(α)。

单调性检验方法算法实现

步骤1 由系统矩阵G的SVD (Singular Value Decomposition)[17]分解得系统矩阵最大奇异值σ1，给出初始参考向量α=[σ1/50, 2σ1/50, …,σ1]。

步骤2 计算H(α)函数值构成的向量H=[H(α1),H(α2),…,H(α50)]。

步骤3 如果H中的最小分量恰为第一个分量，更新参数向量为α1=α/2。

步骤4 重复步骤2、步骤3直到向量H中的第一个分量不再是最小分量(假设此时H中的最小分量为第j个分量，其中2

4 数值仿真及试验验证

采用悬臂梁有限元模型用来模拟将在试验验证中使用的图7所示真实钢制悬臂梁结构。将此悬臂梁划分为10个单元，10个节点从左到右依次编号1-10。假设悬臂梁结构所对应系统为线性系统，且在载荷加载到结构上之前初始条件为0，如图1所示。

图1 悬臂梁有限元模型

另外用于载荷识别的加速度响应数据假设受到环境测量噪声的污染，所以建立给数值计算的加速度响应中添加噪声的公式模型为

y*=ytrue+lnonoisestd(ytrue)

(17)

式中：y*为用于载荷识别的加速度响应；ytrue为数值计算得到的响应；lno为噪声水平；noise为均值为0方差为1的白噪声序列；std(ytrue)为真实响应ytrue的标准差。识别载荷与真实载荷之间的相对误差及相关系数计算公式为

(18)

(19)

图2 积分滑动平均后残留低频噪声

在数值仿真算例中取基函数向量为q(t)=[1,t,t2]T构造式(12)中离散系统方程，使用L∞范数拟合正则化方法和L2范数正则化方法进行载荷识别，分别使用单调性检验和GCV (Generalized Cross-Validation)方法来确定最优正则化参数，最后将两种正则化方法所得载荷识别结果进行比较分析。

4.1 正弦载荷识别

在第8节点处施加正弦载荷p1(t)=40sin(40πt)，载荷作用时间为0.3 s。采样频率为1 000 Hz，使用第4节点处的加速度响应分别在5%、10%及15%噪声水平下进行载荷识别，积分滑动平均次数为30。15%噪声水平下，函数H(α)图像，如图3所示，图3中星号处表示最优正则化参数所在位置。两种正则化方法所得识别载荷，如图4所示。不同噪声水平下正弦载荷识别相关结果见表1。由图4可知，在正弦载荷识别过程中，单调性检验方法及GCV方法均可以有效确定两种正则化方法的最优正则化参数，相应识别载荷较好的反映了真实载荷的时间历程。另外，不难看出在正弦载荷的峰值附近L∞范数拟合正则化方法所得识别载荷与真实载荷有更好的拟合性质。

图3 15%噪声水平下函数H(α)

图4 15%噪声水平下识别正弦载荷

表1 正弦载荷识别结果

Tab.1 Results of the identified sinusoidal load

噪声水平L∞范数L2范数RErr/%CCoeRErr/%CCoe5%8.050.999413.900.998610%8.870.998413.970.998015%10.330.996814.440.9969

4.2 冲击载荷识别

半正弦波可以用来表示典型的冲击载荷[18]，所以在第8节点处施加冲击载荷

(20)

整个时间历程为0.03 s。采样频率为5 000 Hz，使用第6节点处长度为0.03 s的加速度响应，分别在5%，10%及15%噪声水平下识别冲击载荷，积分滑动平均次数为30。冲击载荷的主要特征为最大幅值[19]，所以文中冲击载荷对应的RErr指的是识别载荷最大幅值与真实载荷最大幅值之间的RErr。

5%噪声水平下两种正则化方法所得识别载荷，如图5所示。不同噪声水平下对冲击载荷识别相关结果见表2。由图5可知，L∞范数拟合正则化方法所得识别载荷较好的反映了冲击载荷的时间历程，与真实载荷有较好的拟合性质，识别载荷最大幅值与真实冲击载荷最大幅值非常接近。L2范数拟合正则化方法对应的识别载荷最大幅值与真实幅值差别相对较大。相关识别结果见表2。经计算真实冲击载荷最大幅值为39.92。

图5 5%噪声水平下识别冲击载荷

表2 冲击载荷识别结果

Tab.2 Results of the identified impact load

噪声水平5%10%15%L∞范数L2范数幅值39.2539.1639.04RErr/%1.671.912.20幅值38.1738.1338.09RErr/%4.394.494.58

综合数值仿真中使用两种不同正则化方法针对不同形式载荷的识别结果可得，单调性检验方法能够有效确定L∞范数拟合正则化方法的最优正则化参数，得到理想精度的稳定数值解。利用积分滑动平均有效抑制白噪声中的高频噪声部分，对于残留的低频噪声部分L∞范数拟合正则化方法与L2范数拟合正则化方法相比较有更好的适应性，L∞范数拟合正则化方法能够得到精度更高光滑性质更好的识别载荷。同时也说明与L2范数正则化方法相比较，L∞范数拟合正则化方法对于结构响应中的测量噪声有更好的稳定性。

图6 不同模型误差下载荷识别结果

由式(12)中离散系统方程模型中系统矩阵G是通过有限元方法计算的形函数响应所构造，形函数响应中可能会包含结构有限元模型建立过程中所带来的模型误差。因此，考察两种正则化方法对于模型误差的稳定性是非常重要的。初始条件保持不变，模型误差添加公式为Gerror=G(1+Δ)，其中Δ表示模型误差水平[20]。结果如图6所示，如果以冲击载荷识别误差在10%以内为满足我们需要的标准，则当模型误差在7%以内时L2范数正则化方法可以得到理想的载荷识别结果，而L∞范数拟合正则化方法对于10%以内的模型误差均可以得到理想的冲击载荷识别结果。图6中结果说明L∞范数拟合正则化方法对于模型误差有比L2范数正则化方法更好的稳定性。

4.3 试验验证

悬臂梁结构实验图片如图7所示，基本假设及初始条件与数值模拟中相同。悬臂梁左端固定右端自由。规格0.9 m×0.05 m×0.009 m，密度7.8×103kg/m3，弹性模量E=200 Gpa。将此悬臂梁划分为10个单元，如图1所示，从左到右依次为节点1～10。在第8节点处施加频率为20 Hz的正弦形式载荷并同时测量真实载荷数据，载荷作用时间为0.5 s，在第3～第6节点处放置加速度传感器测量结构响应，采样频率为1 024 Hz。

图7 悬臂梁结构

使用第6节点处加速度响应识别载荷，如图8所示。L∞范数正则化方法识别载荷RErr及CCoe分别为

图8 试验验证识别载荷

17.46%及0.986 2。L2范数正则化方法识别载荷RErr及CCoe分别为17.51%及0.986 0。识别载荷初始时刻附近的波动原因为载荷的突然加载产生的瞬态脉冲响应造成的干扰。实验室环境下，环境带来的测量噪声及数据采集系统噪声都处于非常低的水平，所以使用两种正则化方法进行载荷识别所得结果的误差及精度几乎没有差别。

5 模拟遥测数据冲击载荷识别

飞行器结构在点火、分离瞬间关键部位所受冲击载荷可能导致飞行器结构严重损伤，所以飞行器结构冲击载荷识别问题是一类典型的非常重要的载荷识别问题。建立飞行器的三维集中质量梁模型，如图9所示，集中质量三维梁模型共34个质量站点，建立34个0D单元属性，mass_1到mass_34，分别赋给相对应的站点建立34个point单元。按照给定的舱段划分，再建立12个1D单元属性，liang_1～liang_7及liang_9～liang_13分别赋给相对应的舱段建立33个bar单元。构造形如图13中实线所代表的模拟冲击载荷，用此模拟冲击载荷来模拟飞行器级间分离时刻结构y方向所受典型冲击载荷。将模拟冲击载荷作用于飞行器级间分离位置所对应的质量站点，使用飞行器模拟遥测加速度响应来识别模拟冲击载荷。

图9 集中质量梁模型

遥测数据采集过程中，由于数据采集系统精度的限制导致采集到的响应数据只能分布在一些特定的数值上，如图10所示。所以遥测数据中所包含的噪声主要来源就是数据采集系统带来的系统噪声。由于受到环境噪声及采集系统噪声的综合影响，导致遥测数据中所包含的噪声并不是高斯白噪声，而是一种水平较高的近似服从均匀分布的随机噪声，如图11中概率密度分布所示。实际情况中飞行器载荷识别问题所使用的响应数据是遥测数据。既往的飞行器冲击载荷识别问题研究中所考虑的噪声基本为高斯白噪声，鲜有考虑噪声为均匀分布或者近似服从均匀分布情况下的飞行器冲击载荷识别问题研究。由Clason提出的L∞范数拟合正则化方法适用与均匀分布的噪声情形，所以我们使用L∞范数拟合正则化方法及模拟遥测数据进行载荷识别。模拟遥测数据采样频率为6 400 Hz，类似数值仿真算例中方法计算形函数响应。使用第17质量站点处y方向的形函数响应及模拟遥测数据构造形如式(12)中离散系统方程。

图10 模拟遥测数据

图11 系统噪声概率密度分布

这里由于系统噪声并非高斯白噪声，而是水平较高的近似服从均匀分布的随机噪声。对于此种噪声积分滑动平均并不能有效的抑制噪声，所以取积分滑动平均次数n0=0，也即不对模拟遥测数据进行积分滑动平均。分别用两种正则化方法进行载荷识别。单调性检验函数H(α)图像，如图12所示，图12中星号代表最优正则化参数所在位置。两种正则化方法识别模拟冲击载荷，如图13所示，L∞范数拟合正则化方法识别冲击载荷在两处最高峰值处的相对误差由左至右分别为9.9%和10.4%。L2范数正则化方法识别结果不理想，说明对于此种近似均匀分布噪声L2范数正则化方法并不适用。L∞范数拟合正则化方法可以有效识别动态载荷，所以与L2范数拟合正则化方法相比较L∞范数拟合正则化方法对于噪声形式有更大的适用范围。

图12 冲击载荷识别单调性检验函数H(α)

图13 识别模拟冲击载荷

6 结 论

在飞行器动态载荷识别问题中，针对较高水平高斯白噪声情形下载荷识别问题，L∞范数拟合正则化方法结合积分滑动平均可以有效识别动态载荷。针对遥测数据中主要由数据采集系统带来的近似服从均匀分布的较高水平系统噪声，L2范数拟合正则化方法具有自身的局限性。然而，L∞范数拟合正则化方法可以有效识别冲击载荷并得到拟合性质较好精度较高的识别结果，所以与L2范数拟合正则化方法相比较L∞范数拟合正则化方法对于结构响应中的测量噪声有更好的适用范围及稳定性。与L2范数拟合正则化方法相比较，L∞范数拟合正则化方法同样对于系统模型误差有更好的稳定性。因此，综合数值模拟及试验验证中动态载荷识别结果，与L2范数拟合正则化方法相比较L∞范数拟合正则化方法具有更大的适用范围及更好的稳定性。

[1] 郭荣, 房怀庆, 裘剡, 等. 基于Tikhonov正则化及奇异值分解的载荷识别方法[J]. 振动与冲击, 2014, 33(6): 53-58.

GUO Rong, FANG Huaiqing, QIU Shan, et al. Novel load identification method based on the combination of Tikhonov regularization and singular value decomposition[J]. Journal of Vibration and Shock, 2014, 33(6): 53-58.

[2] 杨武, 刘莉, 周思达, 等. 前后向时间序列模型联合估计的时变结构模态参数辨识[J]. 振动与冲击, 2015, 34(3):129-135.

YANG Wu, LIU Li, ZHOU Sida, et al. Modal parameter identification of time-varying structures using a forward-backward time series model based on joint estimation[J]. Journal of Vibration and Shock, 2015, 34(3):129-135.

[3] 马超, 华红星. 基于改进正则化方法的状态空间载荷识别技术[J]. 振动与冲击, 2015, 34(11): 146-149.

MA Chao, HUA Hongxing. State space load identification technique based on an improved regularized method[J]. Journal of Vibration and Shock, 2015, 34(11): 146-149.

[4] GHAJARI M, SHARIF-KHODAEI Z, ALIABADIL M H, et al. Identification of impact force for smart composite stiffened panels[J]. Smart Materials and Structures, 2013, 22(8): 085014-085026.

[5] 马超, 华红星. 一种基于新的正则化技术的冲击载荷识别方法[J]. 振动与冲击, 2015, 34(12): 164-168.

MA Chao, HUA Hongxing. Impact force identification based on improved regularization technique[J]. Journal of Vibration and Shock, 2015, 34(12): 164-168.

[6] INNOUE H, HARRIGAN J J, REID S R. Review of inverse analysis for indirect measurement of impact force[J]. Applied Mechanics Reviews, 2001, 54(6): 503-524.

[7] THITE A N, THOMPSON D J. The quantification of structure-borne transmission paths by inverse methods, Part1: improved singular value rejection methods[J]. Journal of Sound and Vibration, 2003, 264(2): 411-431.

[8] DOYLE J F. A wavelet deconvolution method for impact force identification[J]. Experimental Mechanics, 1997, 37(4): 403-408.

[9] SANCHEZ J, BENAROYA H. Review of force reconstruction techniques[J]. Journal of Sound and Vibration, 2014, 333(14): 2999-3018.

[10] MAO Y M, GUO X L, ZHAO Y. A state space force identification method based on Markov parameters precise computation and regularization technique[J]. Journal of Sound and Vibration, 2010, 329(15): 3008-3019.

[11] CLASON C.L∞fitting for inverse problems with uniform noise[J]. Inverse Problem, 2012, 28(10): 104007.

[12] GORRY P A. General least-squares smoothing and differentiation by the convolution (savitzky-golay) method[J]. Analytical Chemistry, 1990, 62(6): 570-573.

[13] LIU J, SUN X S, HAN X, et al. A novel computational inverse technique for load identification using the shape function method of moving least square fitting[J]. Computers and Structures, 2014, 144: 127-137.

[14] CLASON C, CHRISTIAN, KUNISCH K, et al. Minimal invasion: an optimalL∞state constraint problem[J]. ESAIM: Mathematical Modelling and Numerical Analysis, 2010, 45(3): 505-522.

[15] HINTERMÜLLER, LIER M, ITO K, et al. The primal-dual active set strategy as a semi-smooth Newton method[J]. Journal on Optimization, 2003, 13(3): 865-888.

[16] ULBRICH M. Semi-smooth Newton methods for operator equations in function space[J]. Journal on Optimization, 2003, 13(3): 805-842.

[17] GOLUB G H, HEATH M, WAHBA G. Generalized cross-validation as a method for choosing a good ridge parameter[J]. Technometrics, 1979, 21(2): 215-223.

[18] CHOI K, CHANG F K. Identification of impact force and location using distributed sensors[J]. AIAA Journal, 1996, 34(1): 136-142.

[19] YAN G, ZHOU L. Impact load identification of composite structure using genetic algorithms[J]. Journal of Sound and Vibration, 2009, 319(3/4/5): 869-884.

[20] ZHOU K, DOYLE J C, GLOVER K. Robust and optimal control[D]. New Jersey: Prentice Hall, 1996.

Application ofL∞norm fitting regularization method in dynamic load identification of space-crafts

GAO Wei, YU Kaiping, GAI Xiaonan

(Department of Astronautic Science and Mechanics，Harbin Institute of Technology，Harbin 150001, China)

Based on the variable limit integral theory, the weighted variable limit integral of response function was constructed with the least squares estimation technique. By performing the integral moving average method several times, the high-frequency noise in the measured noise was filtered effectively. Aiming at the residual low-frequency noise in the measured response, theL∞norm fitting regularization method was used for load identification. In order to realize this regularization method, a monotonicity validation method was proposed to select the optimal regularization parameters. The numerical simulation and the test verification showed that the proposed monotonicity validation method can determine the optimal regularization parameters of theL∞norm fitting regularization method effectively; in addition, the identified load is more accurate and more reliable than that using the traditionalL2norm regularization method; aiming at the characteristics of noise in telemetry data, the impact load is identified effectively by using theL∞norm fitting regularization method and simulated telemetry data.

load identification;L∞norm regularization method; moving average; telemetry data

国家自然科学基金资助(11372084)

2015-11-13 修改稿收到日期：2016-02-25

高伟 男，博士生，1981年生

于开平 男，教授，1968年生

O327

A

10.13465/j.cnki.jvs.2017.09.016