二阶运动范文

二阶运动范文（精选7篇）

二阶运动 第1篇

关键词：活塞销偏置,活塞二阶运动,热变形,弹性变形

1 引言

活塞是往复活塞式内燃机的关键零部件之一, 其工作情况不仅对内燃机的工作可靠性、经济性有直接影响, 还对内燃机的振动噪声特性产生一定影响。活塞敲击噪声是由活塞和气缸套的接触引起, 即由活塞二阶运动引起。因此, 从控制内燃机噪声的角度来讲, 探究活塞结构对活塞二阶运动的影响特性很有必要。鉴于活塞结构参数较为复杂, 本文仅就活塞销偏置对活塞二阶运动的影响特性进行研究。

2 基本理论

2.1 活塞二阶运动理论

内燃机活塞运动研究始于1970年代至1980年代, 研究人员认识到活塞二阶运动的存在。1987年, K.P.Oh等建立了第一个包括弹性流体润滑的活塞动力学模型, 使得活塞二阶运动研究理论得到了显著发展。

1992年, Zhu D和Cheng H S等[1,2]在分析活塞二阶运动过程中, 提出了混合润滑模型。该研究使活塞二阶运动研究得到进一步发展, 并成为之后活塞二阶运动研究的主要依据。

依据Zhu D和Cheng H S提出的混合润滑模型, 奥地利AVL公司开发了针对活塞二阶运动研究的Glide模块, 并得到广泛应用。AVL Glide在计算过程中能够充分考虑活塞裙部轮廓、活塞销偏置、活塞热变形和弹性变形等因素。

2.2 计算机辅助技术

计算机辅助技术是采用计算机为工具, 将计算机用于产品的设计和制造等过程, 辅助人们在特定应用领域内完成任务的理论、方法和技术。包括CAD、CAE和CAM等。

CAD技术是近几十年来发展起来的一项新技术。目前已在制造业等行业获得广泛应用。当前, 著名的CAD软件主要有AutoCAD、SolidWorks、Pro/E和UG等。

CAE的理论基础主要是有限元方法。其基本思想是将连续物体离散为有限个单元, 将连续无限自由度问题转化为离散的有限自由度问题。CAE技术现已在内燃机活塞研究中获得广泛应用。当前, 著名的CAE软件包括ANSYS、NASTRAN和IDEAS等。

3 研究过程

3.1 研究对象

本文以某型柴油机活塞为研究对象。该柴油机活塞缸径为128mm。该柴油机在2000r/min转速下的展开示功图如图1所示。

3.2 活塞弹性变形分析

在Glide软件中, 活塞的弹性变形用活塞刚度矩阵表示。按照Glide软件中活塞刚度矩阵的编写方式, 取0° (主推力面一侧型线) 、30°和60°三条纵向型线上的若干点, 作为载荷施加点和变形值读取点, 并对各点进行编号, 见图2。

活塞刚度矩阵的具体求解过程为:首先, 在1点分别施加2000N和5000N的径向力, 并分别记录图2所标记各点的径向变形值;之后, 依次分别对其余各点施加相同大小的力, 同样记录各点变形值;最后按软件要求格式编写活塞刚度矩阵文件。

以下仅详述在1点施加2000N径向力的情况下, 各点径向变形值的求解过程。

本文所分析活塞的结构具有对称性, 但考虑到载荷的分布不具有对称性, 所以仍取全活塞模型进行有限元分析。利用Pro/E软件建立全活塞实体模型, 并将其导入ANSYS软件, 使用SOLID 95单元对其进行网格划分, 得结果如图3所示。在柱坐标系状态下, 限制活塞顶部凹坑的上下运动以及约束对称面的旋转, 之后在1点施加2000N的径向力。然后进行求解, 读取求解结果并绘制活塞径向变形云图, 如图4所示。其中, 蓝色部分表示活塞径向尺寸变小区域, 而红色部分表示活塞径向尺寸变大区域。由图可知, 蓝色区域的变形量明显大于红色区域的变形量, 因此可判断分析结果与实际较为相符。最后, 记录所需各点的变形值。

其余情况的求解过程与之相近, 可依据上述方法分别进行求解, 此处不再赘述。最终, 由求得结果编写活塞刚度矩阵。

3.3 活塞热变形分析

依据文献[3]所采用方法, 对活塞热变形进行分析。所得活塞温度场分布情况和活塞径向热变形情况分别如图5和图6所示。所需各点的径向变形量可从活塞热变形分析结果中读取。

3.4 活塞动力学建模

利用AVL Glide建立活塞动力学模型。分别取活塞销偏置量为-2mm、-1mm、-0.5mm、0mm、0.5mm和1mm (活塞销偏置量以偏向主推力面一侧为负值) 进行研究。针对不同的活塞销偏置量, 分别建立相应的事件 (case) 。通过对比分析各事件中活塞二阶运动结果, 即可分析活塞销偏置对活塞二阶运动的影响。

4 研究结果及分析

依据计算结果, 分别绘制不同活塞销偏置量下的活塞横向位移对比图如图7所示;不同活塞销偏置量下的活塞摆角对比图如图8所示。

分析上述两图, 可得以下结论:

(1) 活塞销偏置对活塞摆动的影响比活塞销偏置对活塞横向运动的影响更为显著, 所以主要考察活塞销偏置对活塞摆动的影响;

(2) 在压缩和做功行程, 活塞销偏置对活塞摆动影响较显著, 在压缩上止点前后影响最大, 且在压缩上止点处, 活塞销偏置方向影响活塞摆动方向, 活塞销偏置量与活塞摆角呈正相关关系;

(3) 活塞在压缩上止点换向时, 活塞销采用零偏置活塞以几近平移的方式运动到主推力面侧, 会产生拍击对噪声控制产生不利影响;活塞销采用正偏置, 活塞摆角为负, 较硬的活塞头部最先与气缸套接触, 对噪声控制亦会产生不利影响;活塞销采用负偏置, 活塞摆角为正, 较软的活塞裙底部最先与气缸套接触。因而, 相对零偏置和正偏置, 负偏置对活塞敲击噪声的控制最有利。

5 结语

综上分析可知, 仅从活塞敲击噪声控制的角度考虑应优先选择活塞销负偏置方式。对本文研究对象而言, 选择-0.5mm作为最优偏置量。

实际设计过程中, 活塞偏置量的选择不仅涉及换向时的噪声控制问题, 偏置量对活塞组的摩擦损耗等影响也应考虑在内。

参考文献

[1]ZHU D, CHENG H S, et al.A numerical analysis for piston skirts in mixed lubrication:partⅠ-basic modeling[J].Journal of Tribology, 1992, 114:553-562.

[2]ZHU D, HU Y Z, et al.A numerical analysis for piston skirts in mixed lubrication:partⅡ-Deformation Considerations[J].Journal of Tribology, 1993 (115) :125-133.

二阶运动 第2篇

内燃机噪声按产生机理分类主要包括机械噪声、燃烧噪声和空气动力学噪声,机械噪声中活塞敲击噪声是最重要的成分之一[1]。活塞在惯性力和缸内最高燃烧压力的作用下对气缸壁产生明显的敲击现象,缸内最高燃烧压力的周期性及往复惯性力方向不断变化使得活塞对缸套的冲击呈现主推力侧 TS (thrust side)和副推力侧ATS(anti thrust side)交替作用的特点;周期性的敲击振动传递到气缸体结构表面,从而产生明显的结构辐射噪声。近年来国内外学者对活塞敲击进行了相关的研究工作[2,3,4,5],取得了一定的成果:采用计算机仿真与试验相结合的方法是进行活塞与缸套间动力学二阶运动分析的有效方法。

针对某6缸四冲程柴油机的活塞与缸套系统,利用多体动力学软件Glide软件计算了活塞二阶运动及活塞动态敲击力,利用有限元Nastran软件计算了气缸体结构的振动响应;在试验样机上布置测量系统对气缸体结构两侧在活塞敲击力作用下的时、频振动响应进行测量,并对仿真与试验结果进行了分析和对比验证。

1 活塞二阶运动理论

1.1 活塞体控制方程

活塞体在缸套中运动需要满足的动态载荷相平衡的原则,建立活塞与缸套系统分析坐标系,X方向为活塞往复运动方向;Y方向垂直于气缸中心线,为活塞副推力侧指向主推力侧。活塞体的受力平衡关系如式(1)～(3)及图1所示。

X方向平动:

Y方向平动:

绕活塞销轴线转动:

式中,mc为活塞质量;Θ为活塞体绕活塞销轴线的转动惯量;Fgc为活塞质量力;Fgas为气缸内气体力;Fr为活塞环组与活塞的接触力;Fci为活塞在主推力和副推力侧接触力;Fpin为活塞销上的载荷;Flink为连杆小头处附加载荷;Mc为侧向力引起的活塞体绕销轴扭矩;Mgc为质量引起的活塞体绕销轴扭矩;Mr为活塞环的轴向和径向力引起的活塞体扭矩;Mpc为活塞销处的摩擦力矩。

1.2 活塞缸套间油膜流体动力学方程

活塞与缸套连接副润滑油膜的弹性流体动力学特性对活塞二阶运动及敲击力有重要影响[6]。活塞裙部与气缸之间的流体润滑特性与发动机滑动轴承的流体润滑特性相似:活塞裙部桶形型面在沿气缸壁面滑动时产生剪切效应;活塞体在上下止点发生换向运动时对润滑油膜产生挤压效应。因此,活塞裙润滑油膜压力分布可借用滑动轴承的雷诺方程描述,与滑动轴承油膜压力方程之间的主要区别在于,需要将滑动轴承油膜方程的轴颈旋转坐标系转化为活塞缸套间往复运动的坐标系。液力润滑偏微分方程为:

undefined

式中,x为活塞沿缸套运动方向;ϕ为活塞与缸套挤压周向角度;h为油膜厚度;p为油膜压力;μ为润滑油动黏度;t为时间。

1.3 迭代计算步骤

活塞二阶运动的计算在时间域内采用龙格库塔积分方法进行,时间步长为均匀的曲轴转角步长,燃烧上止点前0～90 °CA的预计算确定出活塞缸套系统动力学的初始条件;预计算开始时刻活塞运动的初始条件是:速度和加速度均为零;活塞缸套之间的弹性力,活塞销与连杆小头的接触计算采用隐式牛顿-拉斐逊法进行。活塞缸套系统动力学特性迭代求解步骤为:(1)活塞连杆部件运动学位置的确定;(2)准静定法求解作用在活塞及连杆上的载荷;(3)牛顿-拉佛逊法迭代求解活塞缸套接触载荷;(4)力及力矩平衡判断与迭代残差收敛判断;(5)迭代求解连杆大小头轴承载荷;(6)力及力矩平衡判断与迭代残差收敛判断;(7)下一时间载荷步开始。

2 活塞体动力学及二阶运动边界条件

2.1 活塞体的径向刚度

由前述的活塞缸套计算理论模型可知:活塞体在运动方程中以集总参数模型描述;工作条件下活塞在上下止点换向时与缸套壁面在油膜压力的作用下发生弹性挤压与接触。因此,需要考虑活塞体在横向挤压载荷下的弹性变形,即径向刚度的影响。

活塞体主要的承载部位在裙部承压面上。采用有限元法建立活塞体有限元模型,采用圆柱坐标系,其原点位于活塞顶当量直径的圆心处,将对称截面绕活塞轴向的旋转自由度进行约束,同时约束住顶面节点的沿活塞轴线的移动自由度;在活塞裙部有限元模型的各节点处加载分布力,进行静态变形计算。

本文采用有限元方法计算活塞径向刚度,从主推面至副推面,从头部到裙部依次布置若干个加载点,在主要的承载面内进行布置。在活塞1/4裙面的子午线方向和周向45°范围内共选取18个加载点Ni,i=1,2,,18,外载取3 600 N;各加载点依次加载,加载面积为8 mm8 mm的小正方形,得到各测点对应的径向变形Si,j,i,j=1,2,,18,表示第j测点加载时,第i测点的径向变形量。由此得到所有测点的刚度,并组成活塞体径向刚度矩阵Kij,如图2所示。

图3为有限元计算得到活塞体刚度结果。活塞体局部刚度从头部至裙部逐渐减小,活塞头部主要承受缸内热负荷,结构刚度较大;裙部主要起到导向作用,且内部需要保证连杆小头的运动空间,因此裙部结构通常为刚度较弱的桶型。对裙部三个截面的刚度进行对比:主、副推力侧垂直于活塞销轴线,此侧刚度最弱;从主、副推力侧转过-30°的截面刚度较大;而转过60°的截面刚度在三者中最大。这是由于该截面靠近活塞销座区域,销座处金属的堆积显著增大了结构刚度。

2.2 活塞体裙部型面

现代柴油机活塞的裙部截面通常设计为变椭圆度的桶形型面;垂直于活塞体中心线截面的椭圆长轴指向为主、副推力侧,椭圆短轴指向活塞销轴线方向;椭圆度的定义通常分为单椭圆度、双椭圆度和四椭圆度;为了改善裙部与缸套的配合以及满足数控车床加工的工艺性,我国柴油机行业一般采用双椭圆公式来表示截面的椭圆度。

undefined

式中,ΔR为活塞裙部半径缩减量;G为直径椭圆度;k为双椭圆度修正系数。

对活塞裙部冷态型面的描述只需主、副推力面的型线及其他各截面高度上的变椭圆度即可。本文分析的柴油机活塞主、副推力面的型线采用对称形式,其他角度下的型线由变椭圆度决定,活塞型面由接触区域内的各角度下的型线组成。计算活塞的主、副推力侧及分别转过-30°与60°侧的型线,如图4所示。

2.3 活塞体与气缸套温度场的热态型面

柴油机实际运行条件下活塞体和缸套均承受缸内的热负荷;在缸内热场的作用下活塞体与缸套的热膨胀导致型线的径向增量,从而构成热态型面,对活塞缸套的分析将产生重要的影响。利用硬度塞法测量了活塞体的温度场数据,采用试验热电偶法及热流体动力学计算得到缸套在主副推面侧的温度分布[7]。限于篇幅不再赘述。活塞体及缸套在主副推力侧的温度变形与冷态型面及缸套安装变形按式 (6) 叠加,构成了活塞体裙部和缸套的热态型面。

式中,ΔRt为热态型面的半径缩减量;ΔR为冷态半径缩减量;α为材料工作温度下导热系数;ΔT为工作温度与室温间温度差。

图5为活塞与缸套冷热态配缸间隙分布情况。由图5可见:活塞主、副推力侧型线在热态和冷态下与气缸套的配合间隙差异较大,后续计算中须采用热态型面;缸套温度场自上而下递减梯度明显,导致热态型线下活塞裙部在上止点处与缸套间隙较大,上止点处桶形裙部与缸套的热态配缸间隙为35～40 μm之间;下止点时配缸间隙减小,但未发生下止点卡死活塞的现象,下止点时热态配缸间隙保持在20～25 μm之间,可保证活塞裙部润滑油膜厚度的需要。

3 活塞二阶运动及敲击力计算结果

3.1 活塞径向位移与绕活塞销转动倾角

活塞在气缸中除了沿气缸中心线的往复运动之外,在上下止点以及中间行程处,在外部载荷的作用下发生接触面瞬间换向的二次运动现象;配缸间隙的存在,活塞二次运动包括两个不同形式的运动:沿缸孔径向的平动和绕活塞销轴线的转动;活塞径向平动以径向位移表示,活塞绕活塞销的转动以倾角表示。图6和图7分别显示了一个发动机循环内活塞体二阶运动的径向位移与倾角曲线。

3.2 活塞二阶运动姿态

活塞的二阶运动径向位移曲线(图6)显示,径向位移在燃烧上止点过后最大,压缩行程末期活塞裙部副推面压在缸套表面上行,在上止点附近最高燃烧压力加载在活塞顶部,活塞从副推面向主推面换向,且最高燃烧压力产生较大侧向力,因此活塞径向位移曲线在上止点附近除了配缸间隙较大外,活塞主推面发生的弹性变形也有重要贡献。

活塞的倾角曲线(图7)显示,在0、180、360、450、540 °CA附近区域发生活塞绕活塞销的转动,如图8所示在0、180、360和450 °CA时刻活塞倾角的姿态图。该活塞结构未使用活塞销偏置设计,因此,在0 °CA附近倾角较小;而其他换向点处均有明显的倾角发生,特别是360 °CA换气上止点附近的倾角主要由往复惯性载荷导致,同时上止点处的热态配缸间隙较大,引起倾角最大达0.12°;中间行程时发生的倾角由往复惯性力改变方向引起,但中间行程时活塞速度较大,这些时刻发生的二次运动通常对活塞缸套之间的流体动力润滑特性影响较小。

3.3 活塞二次运动的动态敲击力

在缸内最高燃烧压力和往复惯性载荷的激励,以及连杆摆角的影响下,活塞与缸套之间的侧向力通常可按照曲轴转角下的几何条件及参数计算出侧向力,如图9所示。没有考虑活塞缸套之间的流体润滑特性和配缸间隙等因素,计算出的侧向力只是准静态载荷(图9a)。实机工作条件下,除了准静态载荷外,活塞二阶运动导致的动态敲击力和准静态侧向力载荷共同构成活塞对缸套结构的侧向力,本文在多体动力学仿真环境考虑了活塞二阶运动的边界条件和间隙,从而得到活塞缸套间实际侧向力(图9b),时域曲线上明显叠加了冲击载荷信号,这些正是活塞对缸套的动态敲击载荷。将动态敲击力从总的侧向力中分离出来(图9c),其中燃烧上止点处的动态敲击力最大,这主要是此时活塞径向运动引起;其他换向时刻产生的动态敲击力较小,主要是由活塞换向时转动倾角和径向运动产生。

4 活塞敲击气缸体振动响应的试验与仿真

对活塞敲击引起的结构表面振动响应采用试验方法进行了测量。试验机型为6缸四冲程柴油机,分别在气缸体第6缸结构主敲击侧和副敲击侧布置振动加速度传感器,测量了发动机在最大扭矩工况(1 200 Nm、1 500 r/min)下的振动加速度时域信号;同时,在试验样机上利用缸压传感器测量了最高燃烧压力曲线以供仿真模型使用;布置上止点信号传感器可判断敲击时刻的出现。振动测量布置如图10所示。

利用多体动力学仿真得到缸套主、副推力面高度方向各9个截面的总侧向力(包含动态敲击力),将其作为外载施加在气缸体结构的有限元模型上,气缸体有限元模型截面加载如图11所示。缸套在主副推侧由上至下各取9个加载点,各截面上的敲击力按活塞与缸套间瞬态接触压力及接触面积积分得到。除了侧向力外,主轴承孔施加最大扭矩工况下的轴承载荷。气缸体的位移边界条件在上、下底面采用Z向对称约束,台架对气缸体的支撑处约束三向平动位移。有限元计算利用Nastran软件,采用直接响应法进行,计算得到气缸体有限元结构对应测量部位的振动响应,包括时间域和频率域的振动加速度结果。

图12为活塞缸套主推面测点的振动加速度时域结果。计算结果以第6缸上燃烧0 °CA时为起始点标记。通过对比试验与仿真曲线可知:活塞在燃烧上止点发生对缸套主推面的冲击强烈,较好地表现了发生强烈敲击的时刻和幅值;活塞经过燃烧上止点后还会在惯性力的作用下发生换向,以及其他各缸活塞在上止点处的敲击也会传递到测量点处,因此,在时域曲线上出现其他多处幅值较小的冲击;但对主推力面振动响应而言,测量点对应的气缸燃烧上止点处发生的冲击最为剧烈。

图13为气缸体副推面测点的加速度时域结果。副推面发生的振动响应在换气上止点处最大,此时活塞发生较大的倾角;同时裙部从主推面瞬间敲击副推面,产生较大冲击。副推面冲击幅值比主推面的明显要小,同时副推面测点布置位置更易受到其他各缸的影响。副推面时域振动加速度清楚地反映出发动机循环内各个缸内产生的敲击载荷。试验曲线和仿真曲线对应较好,可以反映出一个工作循环下缸体副推面处的振动响应。

对主推力面振动时域曲线进行FFT频谱转换后得到频率域振动响应结果曲线,如图14所示。频谱曲线试验结果与仿真结果在500～3 500 Hz吻合得较好,敲击力对结构振动响应在中高频率(1 000～3 000 Hz)产生了重要的贡献。由此可见:活塞敲击导致气缸体结构在中高频域产生辐射噪声,对整车通过噪声和车内噪声均有重要影响。

副推面频域振动响应试验与仿真曲线如图15所示。在500～3 500 Hz两者吻合较好。副推面较大的振动响应出现在1 500～2 500 Hz频段范围内,这与活塞侧向敲击力的主要频率成分及对应气缸体结构副推面处的局部模态特性密切相关。在主、副推面频域振动响应曲线中,低于500 Hz的频率范围不是活塞敲击的主要频率域;同时仿真模型有限元计算时忽略了其他结构的影响,如缸体、缸头、支架等连接条件,从而导致低于500 Hz时仿真结果与试验结果之间有一定误差。

5 结论

(1) 活塞动力学及二阶运动分析中,活塞体径向刚度、活塞缸套的热态型面及热态配缸间隙是仿真计算的重要输入条件。

(2) 活塞对缸套的总侧向力包含准静态分量和动态分量,动态敲击力包含重要的中高频冲击与噪声成分,因此,气缸体在敲击力的响应和结构辐射噪声主要受到动态敲击力的影响。

(3) 多体动力学仿真得到活塞在热态配缸间隙的影响下发生的二阶运动,可通过二阶运动曲线与动态判断活塞动态敲击力的产生因素及各曲轴转角下的活塞运动轨迹。

(4) 结合有限元时域瞬态响应和频率响应计算,可得到气缸体上部结构两侧振动响应;通过试验实测结构振动加速度响应,验证了仿真方法的有效性。

参考文献

[1]庞剑,湛刚.汽车与发动机振动噪声理论与应用[M].北京:科学出版社,2006.

[2]Rebecca M,Sudjianto A,Du X P,et al.Robust piston designand optimization using piston secondary motion analysis[C]//SAE 2003-01-0148,2003.

[3]Kazuhiro N,Yuji Y,Kyochi S.Approch to minimization ofpiston slap force for noise reduction by numberical simulation[C]//JSAE Review 9930720,1999.

[4]Nobutaka T,Takayuki K,Keita H,et al.Optimization of pro-file for reduction of piston slap excitation[C]//SAE 2004-01-4309,2004.

[5]Mansouri S H,Wong V W.Effects of piston design on pistonsecondary motion and skirt-liner friction[C]//SAE 2004-01-2911,2004.

[6]Goenka P K,Meernik P R.Lubrication analysis of piston skirts[C]//SAE 920490,1992.

二阶波浪力的仿真计算 第3篇

当今的水下机器人很多都是在深海中运动, 深海中是比较平静的, 受到海浪的影响并不大, 而浅水域中水下机器人的受力情况却大不一样。当水下机器人在浅水域航行时, 主要的受力是由环境引发的载荷, 如风、波浪和水流, 其中最重要的是波浪力。波浪引起的力可以分解为两部分:①具有波频率的一阶振荡力, 引起一阶振荡运动。一阶波浪干扰力是由高频小振幅具有振荡特性的波浪产生的, 它主要引发水下机器人的纵摇和垂荡运动, 对横摇的影响稍次之[1]。②二阶波浪漂移力 (简称二阶波浪力) , 是由不规则波的长周期漂移运动引起, 由于运动的平衡位置不在浮体原先的平衡位置上产生, 所以它与一阶波浪干扰力相比, 通常具有较小的数量级。但是, 这个力特别是在谐振的情况下, 会引起相当大的附加应力[2], 二阶波浪力的研究对于浅水水域水下机器人的姿态保持是有重要意义的。

计算波浪力的方法多种多样, 著名的有STF方法[3], STF方法是从三维物体在波浪中运动的一般场合出发, 对物体在流场中的流体动力和波浪扰动力进行推导, 在运动方程流体力学系数的表达式上反映了航速影响和物体对波浪流场的绕射影响。后来有STF理论和Frank源分布紧密结合的一种算法[4,5], 可以方便地在任意三维物体剖面上布置计算节点。其中二阶波浪力与波高平方成正比, 而波高随深度呈指数规律衰减。研究二阶波浪力的计算是离不开对二阶波浪的研究的, 本研究主要探讨二阶波浪力的仿真计算。

1 二阶波浪的仿真计算

水面波浪可以看做是一个平稳的随机过程, 它由多个 (理论上是无穷多个) 不同周期和不同随机初始相位的余弦波叠加而成。研究水平面上某一点的波高, 在仿真条件下, 考虑到实际应用, 可以忽略高次谐波, 则有:

式中 η (x, t) 固定点x处的波动水面相对于静止水面在t时刻的瞬时高度;ai第i个组成波的振幅;x, t位置和时间;εi第i个组成波的初始相位, 是在 (0～2π) 上均匀分布的随机变量;ωi, ki第i次组成波频率和波数, 分别为:ωi=2π/Ti;ki=2π/λi=ω${}_i^2$/g。

为了能够使平面上多个方向的波叠加形成波浪场, 引入了波浪频谱的概念, 在水下机器人运动控制研究中应用比较普遍的是ITTC推荐的波谱, 它是基于Pierson-Moscowitz谱 (简称PM谱) [6], 计算公式如下:

其中, A=0.78, B=3.11H${}_{1/3}^2$, H1/3为有义波高。

当水下机器人以一定的航速和航向在水上航行时, 波浪实际作用于船上的频率已不同于自然频率, 而是所谓的遭遇频率ωe。设水下机器人以速度U运动, 遭遇浪向角β (波浪传播方向与水下机器人航向间的夹角) 航行时, 遭遇频率ωe与自然频率ω有如下的关系:

根据能量等效原理, 自然频率下微元面积能量转换为遭遇频率下的微元面积能量并没有改变, 因此有Sζ (ω) dω=Sζ (ωe) dωe由式 (3) 可得遭遇频率波能谱为:

2 二阶波浪振幅的仿真

2.1二阶波浪的仿真过程

(1) 选择要仿真的波谱, 确定波浪的有关参数 (如有义波高) 和水下机器人的航行参数 (如航速, 遇浪向角) , 这里取H=0.8, U=5, β=90°, ω1=0.12, ω2=0.98, t=1∶1∶1 500;

(2) 对波谱进行离散化, 离散的方法很多, 这里采用等间隔法时;

(3) 将各个谐波按式 (1) 进行叠加即可得到仿真的海浪, 其中各谐波的初始相位角εi可用一个0～1之间均匀分布的伪随机数程序来产生。

对具有遭遇频率的海浪的仿真, 可由下式得到:

2.2二阶波浪仿真的结果

利用Matlab编程得到的仿真结果如图1、图2所示。

3 二阶波浪漂移力的计算

本研究引入流体力学中的势流理论来讨论波吸力的计算方法[7,8,9,10]。根据势流理论做如下假设:

(1) 波面具有稳定的各态历经性, 流场是均匀的、不可压缩的理想流体的无旋运动。

(2) 把水下机器人当作一个变截面圆柱体且每个微段圆截面视作无限圆柱的一段。

(3) 只考虑波长相对于水下机器人的截面半径为长波的波元, 尤其是波长接近艇长的长峰波元的影响, 不考虑短波的影响。

由势流理论, 通过流场的速度势可以得到流场的速度分布, 利用伯努利方程建立速度场和压力场的联系, 进一步通过积分运算就可求得二阶波浪力和力矩:

其中, Zs, o, Zz, Zθ, Ms, o, Mz, Mθ为波浪干扰力系数, 可分别表示为:

式中 R=R (x) 水下机器人的几何尺寸;β浪向角;h海平面到大地坐标系原点的距离;k波数ρ为海水密度。

设仿真条件:潜深11 m, 遭遇角为0°～180°, 利用式 (8) ～式 (13) , 可得二阶波浪力和力矩系数随浪向角的变化曲线如图3所示。

进步一验证了水下机器人在浅水水域某一深度下受不同浪向角的二阶波浪力的影响, 如图4所示。

4 结束语

本研究针对对于浅水水下机器人影响很大的二阶波浪力进行了分析计算, 先对二阶波浪用Matlab对其进行了仿真计算, 得出一系列的曲线, 得到二阶波浪的信息用于二阶波浪力的计算, 并得出一系列的仿真曲线。由以上的讨论可知, 深度越大, 二阶波浪力越小, 当深度一定而海况变化时, 二阶波浪力随波长的减小而增大, 并到达最大值, 然后二阶波浪力很快下降。另外, 二阶波浪力的大小与浪向角有关系。在以后水下机器人控制系统设计中有一定的参考意义。

摘要：水下机器人在浅水水域受到浪流的影响比较大, 其中波浪所产生的二阶波浪力对水下机器人的姿态影响很重要, 在二阶波浪力的扰动下, 将产生相应的低频非线性运动。而二阶波浪力的计算直接与它的振幅相关。首先研究了二阶波浪的仿真及运算, 然后利用得到的波浪信息进行二阶波浪漂移力的计算并对其进行了仿真, 给出了波浪力参数随浪向角变化的结果。该研究结果为水下机器人控制系统的设计奠定了基础。

关键词：浅水,水下机器人,二阶波浪,二阶波浪漂移力,仿真

参考文献

[1]HIRANO M, TAKASHINA J, TAKAISHI Y, et al.Shipturning trajectory in regular waves[D].Trans.West-Japan-Soeiety of Naval Architects, 1980:17-31.

[2] WU J S, SHEU J J. An exact solution for a simplified model of the heave and pitch motions of a ship hull due to a moving load and experimental results[J].Journal of Sound and Vibration, 1996, 192 (2) :495-520.

[3] SALCESEEN N, TUCKE O, FALTINSEN O. Ship motions and sea loads[M]. SNAME, 1970, 78.

[4]FRANK W.Oscillation of cylinders in or belowthe free sur-face of deep fluids[R].AD-668338.

[5] BEDEL J W, LEE C M. Normerecal calculation of cylinders oscillating in or below a free surface[R]. AD-722690.

[6]FLEISCHER D, PARK S K.Plane hydroelastic beam vibra-tions due to uniformly moving one axle vehicle[J].Journalof Sound and Vibration, 2004, 273 (3) :585-606.

[7] LEE C M, NEWMAN J N. First-and second-order wave effects on a submerged spheroid[J].J.S.R., 1991, 35 (3) :183-190.

[8]SALVESEN N.On Second Order Wave Theory for Sub-merged Two Dimensional Bodies[A].Paper presented at6thNaval Hydrodynamics Symposium[C].Washington, D.C., 1966.

[9] MUSKER A J. Prediction of wave force and moments on a near-surface submarine[J].I.S.P., 1984 (31) :2-12.

二阶切换线性系统的滑模控制 第4篇

切换动态系统是一类重要的混合动态系统，它是由几个连续时间子系统或离散时间子系统以及相应的切换规则构成。许多实际物理系统都可以采用切换系统模型描述，如化学过程、电力系统、汽车控制系统等。由于切换规则的作用，使得切换系统不同于一般的连续时间系统或离散时间系统,其动态特性变得非常复杂，并且具备一些特殊的性质。例如系统的状态轨迹可能会发生跳变，系统的稳定性可能会发生改变。因此，对切换规则的选择非常重要。而切换系统的稳定性的一个显著特点是，其子系统的稳定性不等于整个系统的稳定性，即可能存在这样的情形，切换系统的每个子系统是稳定的，但是在按照某一规则进行切换时，会导致整个系统不稳定;与此对应，也可能存在这样的情形，尽管每个子系统是不稳定的，但是可以通过某种切换规则使整个系统稳定。即使切换系统的每个子系统都是线性定常系统，其整体一般也不为线性系统，而属于非线性系统。切换系统研究已引起国内外学者的广泛关注，在过去三十年中，控制界对切换系统的建模、分析、综合与控制的研究兴趣不断升高。切换系统分析和切换控制的研究越来越受到人们的关注[1,2]。

滑模变结构控制理论(SMC)是变结构控制理论的主要理论体系。它已经形成了一整套综合系统的独立理论，包括滑动模态的设计方法，控制器的各种综合方法，系统的稳定性分析，系统的到达条件等。与此同时，工程实践已经对控制理论提出了更高的要求。控制对象的复杂程度日渐增加，控制系统的非线性日益严重。不确定性系统的控制问题已经被鲁棒控制、自适应控制等理论领域所研究，非线性系统的控制问题已经引起了许多学者的重视;变结构控制理论恰恰是解决非线性系统控制问题很好的一种方法。近年来滑模变结构控制理论研究备受关注，研究结果主要集中在：1)基于各种滑模面的研究;2)滑模到达条件;3)控制器的各种非线性结构;4)变结构自适应控制及其它交叉研究[3,4,5]。滑模控制策略通过控制量的切换使系统状态沿着滑模面运动，使得系统在受到外界扰动的时候具有不变性，因此滑模变结构控制可应用到处理各种非线性系统。滑模变结构控制系统的基本原理在于，当系统状态穿越状态空间的滑动超平面时，反馈控制的结构就发生变化，从而使系统性能达到某个期望指标。滑模变结构控制器的作用就是把系统的状态在有限的时间内驱动并维持在该子流形上。

本文基于滑模控制原理，在充分分析二阶切换线性系统的特性基础上，借助二阶滑模面推导出二阶切换线性系统采用滑模控制时的控制率求解公式。并将其应用于一个具体的仿真实例。实验结果表明，采用滑模控制方法处理此类二阶切换线性系统时具有较好的控制效果，同时显示出良好的鲁棒性和稳定性。

2 问题描述

考虑如下二阶切换线性系统：

其中X∈Rn;切换率σ(t)∈{1,2,,m}为时间t的右连续函数。Aσ(t)、Bσ(t)、Cσ(t)、Dσ(t)为二阶矩阵。

对于上述系统，采用常规控制方法控制这类过程对象时，当系统未发生切换之前，可以获得较好的控制效果，而一旦系统发生了切换，控制作用则可能变得很差甚至是失控。因此我们需要采用一种新方法去控制这类特殊的非线性系统，获得较好的稳定性和鲁棒性。

滑模变结构控制理论恰恰是解决非线性系统控制问题很好的一种方法。由于具有滑动模态的变结构控制在滑模上对外界干扰和系统摄动具有很强的鲁棒性与稳定性，而且可以通过滑动模态的设计获得满意的动态品质。因此，不失为一种很好地处理非线性系统控制问题的有效方法。

3 滑模控制器设计

一般来讲，滑模控制的设计分为两步。第一步即选择合适的滑模面方程，使系统在滑模上的运动符合预期要求;第二步即设计切换控制率，保证到达条件和滑动条件。

对于上述二阶切换线性系统，采用状态反馈滑模控制，定义滑模面为：

其中e=yd-y。由于系统为二阶线性系统，则：

由于Aσ(t)和Bσ(t)的存在，通常意义下的导数不可求，假设所求的Ueq也是随σ(t)的变化而变化，即为Ueqσ(t)，则上式变形为：

若Dσ(t)=0，可以求出：

若Dσ(t)≠0

求出的方程解为：

以上便是二阶切换线性系统滑模控制器控制率求解过程，从中可以看出：二阶切换线性系统的滑模控制率，随着切换率的变化而变化。当系统发生切换时，控制率也发生较大的切换。

4 实例研究

对于如下二阶切换线性系统：

对于上述二阶切换线性系统，分别依据经典控制理论设计PID控制器和滑模变结构理论设计滑模控制器，并逐一应用于上述二阶切换系统的控制。在Matlab环境下运行仿真程序，两种不同控制器的控制曲线如图1、2所示。

通过对比可知，系统在t=5s时，对象特性发生切换，传统PID控制器经过很长一段时间后，输出才稳定到设定值附近;而采用滑模控制器却能很快让系统输出跟踪到设定值。这说明在切换系统控制中，采用滑模控制理论进行分析、设计的滑模控制器可以获得满意的控制效果和较好的鲁棒性。

5 结束语

本文针对一类特殊的非线性系统切换线性系统，基于滑模控制理论，对二阶切换线性系统的滑模控制率进行推导，得出具体的求解公式。通过Matlab仿真，可以看出当系统特性发生切换时，系统输出最终能够稳定在设定值，控制效果较好，并具有良好的鲁棒性和抗干扰性。

由于在推导滑模控制率时，定义的滑模面随着系统阶数越高越复杂。本文只是推导出了二阶切换线性系统的滑模控制率求解。对于更高阶次切换线性系统的滑模控制率求解问题，还需进一步研究。同时，在设计滑模控制器时，如何避免或消除抖振，提高控制的品质，还需进一步探讨。

摘要：当前,切换系统的稳定控制问题作为控制领域的一个研究难点,引起了学术界和工业界的广泛关注。该文针对二阶切换线性系统,考虑切换时对象特性的变化,将其转化为一类特殊的非线性系统,并基于滑模控制理论,设计出适用于二阶切换线性系统的滑模控制器。最后通过具体的实例仿真验证、结果对比,表明滑模控制器比传统PID控制器在二阶切换线性系统中具有更好的控制效果和鲁棒性。

多模光纤线路的组合二阶失真分析 第5篇

近年来,在光通信网的短程应用中,如何提高多模光纤带宽距离积的应用研究很多[1,2,3,4]。为了提高带宽,研究精确的光纤信道传输模型[5]显得尤为重要。常用的方法是基于耦合功率流方程,采用这种方法可以描述数字脉冲的传输[6],但是考虑模拟信号或基带、射频信号传输时,则存在一定的局限。在这种情况下,需要考虑传输电场信号的方法,提出基于多模光纤的电场传播模型,模型中考虑不同的光源损耗,包括时域相干光源、空间相干光源、光源啁啾参数、模内、模间色散和模式耦合等诸多因素的影响[7]。本文基于电场传播模型方法得到了多模光纤线路的组合二阶失真,分析了在多模光纤上传输宽带射频信号的情况,具体分析估计了多模光纤线路中的二阶谐波失真和互调失真,并给出了分析结果,为光通信网络中多模光纤的应用提供了参考。

1 光接收功率分析

基于文献[7]中多模光纤线路的传输函数,线性频响表明了光纤传输宽带无线电信号的潜力。在实际工程应用中,对于在多模光纤线路中传输的模拟信号,分析其非线性的谐波、互调失真是必不可少的,其中组合二阶失真是主要的失真指标之一。采用副载波复用时,调制信号由多路射频信号组成。图1所示为多模光纤线路示意图,该线路由光源、光外部调制器、多模光纤和光检测器组成。假设调制信号由频率为f1和f2的两个射频信号组成,则已调波信号可以表示为

undefined

式中,P为平均传输光功率;m0为调制系数;α为光源啁啾参数。

图1中调制信号通过外部调制器调制之后经过多模光纤传输,光检测器接收到的射频信号中主要包括f1和f2分量,其中Ω1=2πf1,Ω2=2πf2,另外产生了谐波分量2f1和2f2,以及互调分量f2-f1和f2+f1,经过多模光纤线路传输后的接收功率可表示为

undefined

式中,m和n=1M代表模式数目;Cmn为光注入系数,考虑满注入情形则设为1/M;Gmn为不同模式的耦合系数;Umn与注入系数和耦合系数之和有关。式(3)中第1个指数项是由光源啁啾参数引起的载波抑制效应;第2个指数项取决于不同模式群组对信号传输引起的干扰;第3个指数项取决于群延时τn和τm。Amn(Ω1,Ω2)为二阶谐波失真和互调失真的共同影响,如式(4)所示,即

undefined

式中,Amn(Ω1)和Amn(Ω2)表示原信号的传输,展开为式(5);Amn(2Ω*)和Amn(Ω2∓Ω1)分别表示二阶谐波失真和互调失真分量。

undefined

式(4)中涉及到的二阶谐波失真Amn(2Ω*)和互调失真Amn(Ω2∓Ω1)分别由式(6)和式(7)表示,其中Ω*可以表示Ω1或Ω2。

undefined

2 组合二阶失真仿真分析

假设对于所有模式一阶色散参数都是相等的,即βundefined≈βundefined,∀m。由分析可知:一阶色散与调制系数m0成比例,那么二阶谐波失真与mundefined成比例。这就说明增大m0可以提高信噪比,但同时也增大了失真。当Ω1=Ω2时,Ω2+Ω1处的互调失真相当于2Ω2处的失真。

以下分析估计二阶失真的影响,采用在有线电视(CATV)和手机系统中常用的方法。对于CATV系统中每一个信道的组合二阶失真部分的功率与载波功率的比值(CSO/C)一般要求<-53 dBc。谐波失真的CSO1/C和互调失真的CSO2/C分别表示为

undefined

为了估计非线性部分对系统性能的影响,考虑一般频率的倍数,如f1=f和f2=4f1=4f。每一个频带接收功率的频率响应如图2所示,参数设置如下:5 km的多模光纤,分布反馈式激光器的发射波长为1 300 nm,光束的线宽均方根为10 MHz ,光源啁啾值为零,调制系数mo=0.01。

由图2可知,调制信号f1和f2的接收功率为0～-40 dB,而二次谐波分量和互调分量的接收功率为-60 ～ -100 dB。图中多模光纤线路射频信号传输的频率响应与不同自由光谱范围的微波光子横向滤波器的输出类似,没有载波抑制效应。

由式(8)和式(9)分析计算二阶谐波失真和互调失真,结果表明:光谱相当于在基带f1的线性频响,CSO/C的值<-60 dBc 。这就说明在这些频谱范围,谐波、互调失真的影响可以忽略,即当调制频率在线性频率响应的一个谐振频率带时,谐波、互调失真可以忽略。值得注意的是频率响应的3 dB带宽和谐振频率依赖于光纤传播的模式群组数目和模式延时的传播,所以通过改变线路的长度、纤芯渐变折射率剖面分布、发射波长和光束注入情况等可以调整频率响应的范围。另外,在光谱谐振频率之外组合二阶互调失真的阈值范围为-40 ～ -80 dBc,足以应用在工程实际中。通过计算可以得到调制系数m0=0.05时,CSO/C组合失真的取值范围是-30～-70 dBc。

3 结束语

光纤具有很宽的带宽和低的传输损耗,因而适合传输多频道射频信号。本文分析了在多模光纤上传输射频调制信号时,二阶谐波失真和互调失真的影响,分析表明当调制系数m0=0.01时,若CSO/C组合失真<-50 dBc,谐波失真和互调失真可以忽略。采用副载波复用技术,对于在中短距离多模光纤线路中提高传输容量具有广泛的应用前景。

参考文献

[1]王晋飞,胡贵军,曲仁慧,等.多模光纤模式群分集复用系统[J].中国激光,2010,(12):1966-1969.

[2]Kanprachar S,Jacobs I.Diversity of coding for subcar-rier multiplexing on multimode fibers[J].IEEE TransCommun,2003,(51):1546-1553.

[3]Gu X J,Mohammed W,Smith P W.Demonstration ofall-fiber WDM for multimode fiber local area networks[J].IEEE Photon Technol Lett,2006,(18):244-246.

[4]Tyler E J,Kourtessis P,Webster M.Toward terabit-per-second capacities over multimode fiber links usingSCM/WDM techniques[J].J Lightwave Technol,2003,(21):3237-3243.

[5]戴礼明,蔡祥宝.光纤中脉冲传输模型的研究及其软件设计[J].光通信研究,2010,(1):42-44.

[6]Yabre G.Comprehensive Theory of Dispersion in Gra-ded-Index Optical Fibers[J].J Lightwave Technol,2000,(18):166-177.

构造高维二阶可旋转设计点集 第6篇

首先假设实验中有k个因子 (因子已经被编码) , 用k维坐标表示k个因子的实验记录.在k维空间中, N次试验点表示为 (x1u, x2u, , xku) , u=1, 2N, xtu为在第u次实验中因子t的水平.当这些点集满足以下条件时形成一个k因子的二阶可旋转安排:

若点集形成一个二阶可旋转设计, 则还需满足条件

当取等号时, 所有的设计点会在k维的球面上.为方便起见, 当提出一个可旋转设计时, 我们通常令λ2=1.

2. 由k-1维二阶设计点转变到k维二阶可旋转设计点的方法

首先选择k-1维空间的点, 设共有N'个, 即 (x1u, x2u, , xk-1, u) , 1uN'并且满足 (1) , (2) , (3) 和 (4) , 则 (x1u, x2u, , xk-1, u, ±b) , u=1, 2, , N', 共2N'个点.再考虑点集 (0, 0, , 0, ±p) , (0, 0, , 0, ±q) , 加上n0个中心点可形成一个k维空间上的二阶可旋转安排, 共有点数为N=2N'+4+n0.则可得

其中u=1, 2, , N, 要使这些点集满足可旋转的条件, 即会有

此时令λ2=1, 利用 (1) , 可得

由于p2, q2为非负实数, 那么12,

3. 三维到四维空间点的选取

从本节开始我们通过具体的例子对上述方法进行阐述.

在假设没有加入中心的前提下, 我们由三维空间点集构造到四维空间可旋转的点集, 当四维空间点个数为32, 按照上述方法我们可以选取三维空间的点集 (±a, ±a, ±a) , (±c, 0, 0) , (0, ±c, 0) , (0, 0, ±c) , 然后转换为四维空间的点集 (±a, ±a, ±a, ±b) , (±c, 0, 0, ±b) , (0, ±c, 0, ±b) , (0, 0, ±c, ±b) , (0, 0, 0, ±p) , (0, 0, 0, ±q) .

那么=1.55∈[1, 2]满足可旋转条件.我们可以得到b, c, p, q与a的关系.

接下来我们尝试用三维空间中 (±a1, ±a1, ±a1) , (±a2, ±a2, ±a2) 和 (±c, 0, 0) , (0, ±c, 0) , (0, 0, ±c) 依照上述方法去构造四维空间点集.首先要满足

假设a1=a2=a时, 构成四维空间点集有两个 (±a, ±a, ±a, ±b) 和 (±c, 0, 0, ±b) , (0, ±c, 0, ±b) , (0, 0, ±c, ±b) , (0, 0, 0, ±p) , (0, 0, 0, ±q) , 在三维空间没有加入中心点的时候N'=22, =2.5[1, 2].因此这种选择点的方法是不可行的.当我们在三维空间加入两个中心点变换到四维空间时还是中心点, 有N'=24, =1.71∈[1, 2], 此时的点集可以构成可旋转安排.

我们还可以考虑三维空间的点集 (rcostα, rsintα, 0) , (rsintα, 0, rcostα) , (0, rcostα, rsintα) , (±a, ±a, ±a) , (±c, 0, 0) , (0, ±c, 0) , (0, 0, ±c) .其中α=2π/s, t=0, 1, s-1转换到四维空间点集 (rcostα, rsintα, 0, ±b) , (rsintα, 0, rcostα, ±b) , (0, rcostα, rsintα, ±b) , (±c, 0, 0, ±b) , (0, ±c, 0, ±b) , (0, 0, ±c, ±b) , (0, 0, 0, ±p) , (0, 0, 0, ±q) , 三维空间中点集满足

即有3sr4-128a4+16c4=0, 构成四维空间点集后可以得到

这类转换是否可行和s的取值有关, 最终目的要满足=[1, 2].

因此, 对于这种转换方法, 选取正确合适的k-1空间的点集是非常重要的.

4. 四维的二阶设计点集转变为五维二阶可旋转设计点集

考虑四维空间上的没有中心点的二阶设计的点集 (±a, ±a, ±a, ±a) , (±c, 0, 0, 0) , (0, ±c, 0, 0) , (0, 0, ±c, 0) , (0, 0, 0, ±c) 由二阶可旋转安排的要求得

3C=16a4+2c4=48a4c=2a

则C=16a4, A=16a2+2c2=24a2, N'=24

进而=1.78, 因此可以使用上述的方法.接下来我们假设五维空间中二阶可旋转设计的点为 (±a, ±a, ±a, ±a, ±b) , (±c, 0, 0, 0, ±b) , (0, ±c, 0, 0, ±b) , (0, 0, ±c, 0, ±b) , (0, 0, 0, ±c, ±b) , (0, 0, 0, 0, ±p) , (0, 0, 0, 0, ±q) , 根据小节2中的叙述我们得到p2=4.8165a2, q2=31835a2.即p=2.195a, q=1.784a, c=2a, b=0.817a.所以我们得到了一个五维空间上有52个点的二阶可旋转安排.这52个点为

a的值由条件λ2=1所决定.而这些点所在球面的半径分别为, p, q.即2.16a, 2.16a, 2.195a, 1.784a.因此这个安排不是球面的.然而

当n0=0时, N=52, 所以

此时k/ (k+2) =0.7143<λ4/λ22.

则这些点集可以构成我们的可旋转设计.

摘要：实验者在拟合一个响应曲面进行实验数据分析时, 需要对可操作变量的水平进行设定.Box和Hunter (1959) 提出了一种可旋转方式, 使得响应曲面容易被拟合并且可以提供等高线上的球面信息.本文在Draper (1960) 提出的由k-1维二阶设计点集转变到k维二阶可旋转设计点集方法的基础上, 从具体例子出发分析了最初选取点集的重要性, 并详细地描述了由四维点集构造五维点集的过程.

关键词：二阶可旋转设计,响应曲面,点集,转变

参考文献

[1]D.A.Gardiner, A.H.E.Grandage and R.J.Hader:Some third order rotatable designs Institute of Statistics MimeoSeries No.149, Raleigh, North Caroline, 1956.

[2]G.E.P.Box:Multi-factor designs of first order, Biometrika, Vol.39 (1952) , pp.49-57.

[3]G.E.P.Box and D.W.Behnken:Simplex-Sumdesigns:a class of second order Ratatable designs derivable fromthose of first order, Ann.Math.Stat., Vol.31 (1960) , pp.838-846.

[4]Norman R.Draper:Second order rotatable designs inthree dimensions, Ann.Math.Stat., Vol.31 (1960) , pp.23-33.

二阶运动 第7篇

近几十年来,人脸识别是计算机视觉和模式识别领域研究的一个热点问题。特征脸方法[1](PCA方法)是一个有效的特征提取和降低维数方法,被广泛地应用到人脸识别中。然而,由于真实的人脸是很复杂的,如光照变化较大的人脸图像,单一的特征脸集并不能有效地描述人脸图像。为解决此问题,Wang和Tan提出了二阶特征脸方法[2]。

2004年,二维主分量分析[3](即2DPCA方法)被提出,并被应用于人脸识别。随后,文献[4]的研究发现2DPCA仅提取了图像行方向的信息;事实上,也可以将2DPCA用于提取图像列方向的信息。从而文献[4]提出了(2D)2PCA方法。与2DPCA方法相比,该方法进一步降低了特征提取的时间,并且其识别率略高。

但是,二阶特征脸以及(2D)2PCA方法都有一定的局限性。二阶特征脸应用了两次PCA,而PCA基于高维的图像向量空间提取特征,往往需要较多的运行时间;而(2D)2PCA在光照变化较大时,受光照信息的干扰,识别精度不高,这一点可以从本文第4节的实验中看到。

基于二阶特征脸和(2D)2PCA两种方法,本文提出了二阶(2D)2PCA方法。该方法将(2D)2PCA技术分别应用到原始图像矩阵集和剩余图像矩阵集。这样,原始图像和剩余图像的特征矩阵都被提取为图像的代数特征并用到人脸识别中。

实验表明,相对于特征脸、二阶特征脸和(2D)2PCA方法,二阶(2D)2PCA得到了更高的识别精度并节省了计算时间。为了简便起见,我们将上述方法依次记为PCA、Sec-PCA、(2D)2PCA和Sec-(2D)2PCA。

1 Sec-PCA和(2D)2PCA方法回顾

1.1 Sec-PCA方法

文献[2]的实验表明,PCA方法中前几个主元特征主要体现了图像的光照信息。Sec-PCA方法在原始图像向量集中先执行一次PCA,提取特征脸集的前几个反映光照信息的主元后,在余像空间中再次执行PCA。

设X={x1,x2,,xM}是N维原始图像向量集,由特征脸方法可知,通过变换:

y=UT(xi-μ) i=1,2,,M

得p1维特征向量集Y={y1,y2,,yM}。设$B\{\widehat{x}{}_1,\widehat{x}{}_2,\cdots ,\widehat{x}{}_{{\rm M}}\}$对原始图像集X用了一次PCA方法后的重建图像向量集,令$x^{\prime }{}_i=x{}_i-\widehat{x}{}_i‚x^{\prime }{}_i$称为原始图像向量xi的剩余图像向量;作X′={x′1,x′2,,x′M},X′称为原始图像集X的剩余图像集。Sec-PCA方法是在剩余图像集X′上再次执行PCA。这里,一阶主元数p1事先给定,二阶主元数p2按准则:

${\displaystyle \sum _{i=1}^{p{}_2}\lambda }{}_i/{\displaystyle \sum _{i=1}^{{\rm N}}\lambda }{}_i\ge \theta$

选取,其中λ是样本图像集的协方差矩阵的特征值,θ是一给定的门限值。

我们知道,PCA方法需要将二维图像矩阵转化为一维图像向量后在向量空间中执行。通常这些一维图像向量的维数很高,从而PCA的计算量较大。而Sec-PCA在原始图像集和剩余图像集中分别执行一次PCA,所以所需的计算量大,计算时间较长。

1.2 (2D)2PCA方法

1.2.1 二维主分量分析-2DPCA

设I={A1,A2,,AM}是原始图像集,其中Ai(i=1,2,,M)是大小为mn的原始图像矩阵。设A是图像集I中的任一图像矩阵,2DPCA的思想是通过线性变换:

Yk=AXkk=1,,r

将图像矩阵A投影到向量Xk上,获得了一个m维的投影特征向量Yk。令X=[X1,,Xr],Y=[Y1,,Yr]。显然Y = AX,矩阵Y称为A的特征矩阵。最优投影轴Xk是图像协方差矩阵:

$G{}_{row}=\frac{1}{{\rm M}}{\displaystyle \sum _{j=1}^{{\rm M}}(}A{}_j-\overline{A}){}^{{\rm T}}(A{}_j-\overline{A})$

最大的r个特征值所对应的r个特征向量[3]。设θ是给定的门限值,最优轴个数r可按准则${\displaystyle \sum _{i=1}^r\lambda }{}_i/{\displaystyle \sum _{i=1}^n\lambda }{}_i\ge \theta$选取,其中λ是协方差矩阵Grow的特征值。

1.2.2 (2D)2PCA方法

文献[4]的研究表明,2DPCA是从反映图像行间信息的训练样本集中学习到最优投影矩阵X,通过矩阵变换Y=AX,得到md的特征矩阵Y。由此,文献[4]提出了alternative 2DPCA方法。该方法是从反映图像列间信息的训练样本集中学习到最优投影矩阵Z,通过矩阵变换:

B=ZTA

得到qn的特征矩阵B。其中Z=[Z1,,Zc]是图像协方差矩阵:

$G{}_{col}=\frac{1}{{\rm M}}{\displaystyle \sum _{j=1}^{{\rm M}}(}A{}_j-\overline{A})(A{}_j-\overline{A}){}^{{\rm T}}$

最大的前c个特征值所对应的c个特征向量构成。最优轴个数c按如下准则选取:

${\displaystyle \sum _{i=1}^c\lambda }{}_i/{\displaystyle \sum _{i=1}^m\lambda }{}_i\ge \theta$

综合2DPCA和alternative 2DPCA,文献[4]提出了(2D)2PCA方法。该方法作矩阵变换:

C=ZTAX

将mn的图像矩阵A同时投影到最优投影矩阵X和Z上,得到qd的特征矩阵C。识别时,将C作为图像A的代数特征送入分类器。与PCA和Sec-PCA相比,(2D)2PCA基于图像矩阵提取特征,大大降低了计算量。但是对于光照变化较大的人脸识别问题,(2D)2PCA的识别精度较低。这一点可以从第4节的实验中看到。

2 Sec-(2D)2PCA方法

2.1 光照信息

事实上,对于光照变化强烈的图像,(2D)2PCA方法的前几个最大的特征值所对应的主元向量主要体现了光照信息,而并不是人脸身份特征信息。这一点可以从下面的实验中看到。

实验中,我们给出了用(2D)2PCA行方向、列方向上前5个最大的特征值所对应的主元向量的重建图像;其次分别给出了(2D)2PCA的剩余图像,该剩余图像是原始图像与上述重建图像之差。实验结果如图1所示。

由图1可知,重建图像(b)主要体现了光照信息,而包含的的身份特征信息明显很少,而剩余图像(c)则包含了图像大部分的高频身份特征信息。

2.2 Sec-(2D)2PCA方法

本节我们提出了Sec-(2D)2PCA方法。该方法首先执行一次(2D)2PCA,提取前几个反映光照信息的主元向量;然后在剩余图像集中再次使用(2D)2PCA方法。这样的处理对于光照变化较大的人脸图像,可以避免光照信息的干扰。

设I={A1,A2,,AM}是原始图像集,其中Ai(i=1,2,,M)是原始图像。设A是图像集I中的任一图像矩阵。首先,在图像集I上用(2D)2PCA方法,学习到反映行特征信息的最优投影矩阵X1(大小为nr1)和反映列特征信息的最优投影矩阵Z1(大小为mc1)。对任一图像A,作变换 :

C=Z${}_1^{{\rm T}}$AX1

得大小为c1r1的矩阵C,称为图像矩阵A的一阶特征矩阵。这里r1、c1个主元特征主要体现了光照信息,一般是较小值。

易知A的重建图像为$\widehat{A}{}_i={\rm Z}{}_{1}C{}_{i}X{}_1^{{\rm T}}(i=1,2,\cdots ,{\rm M})$。令$A^{\prime }{}_i=A{}_i-\widehat{A}{}_i$,得到剩余图像集I′={A′1,A′2,,A′M}。Sec-(2D)2PCA方法就是在剩余图像集I′上再次使用(2D)2PCA。

由(2D)2PCA方法,在剩余图像集上,再次学习到反映行信息的最优投影矩阵X2(大小为nr2)和反映列信息的最优投影矩阵Z2(大小为mc2)。对任一余像A′,作变换:

C′=Z${}_2^{{\rm T}}$A′X2

得大小为c2r2的矩阵C′,称为图像矩阵A′的二阶特征矩阵。这里,r2和c2分别按准则:

${\displaystyle \sum _{i=1}^{r{}_2}\lambda }{}_i/{\displaystyle \sum _{i=1}^n\lambda }{}_i\ge \theta {\displaystyle \sum _{i=1}^{c{}_2}\lambda }{}_i/{\displaystyle \sum _{i=1}^m\lambda }{}_i\ge \theta$

选取,其中θ是一给定的门限值。图2给出了Sec-(2D)2PCA方法特征提取的流程图。

对于人脸识别问题,设训练样本集为I={A1,A2,,AM},由Sec-(2D)2PCA方法,得一阶特征矩阵集{C1,,CM}和二阶特征矩阵集{C′1,,C′M}。每一个样本被指定为一个特定的类ωk。对于训练样本图像Aj,设相应的特征矩阵为Cj和C′j(j=1,2,,M)。对于一个测试图像At,由Sec-(2D)2PCA方法,分别得到一阶和二阶特征矩阵Ct和C′t。按最小距离法分类,定义距离:

d(At,Aj)=a1 ‖Ct-Cj‖+a2‖C′t-C′j

其中权系数a1和a2分别反映了一阶和二阶特征矩阵的重要程度。而a1=1-a2,因此我们只需确定a2的值。对于直观上光照变化不大的图像库,a2∈[0,0.5];而对于直观上光照变化较大的图像库,a2∈[0.5,1]。如果:

$d(A{}_t,A{}_k)=\underset{j}{{\displaystyle \min }}d(A{}_t,A{}_j)$

且Ak∈ωk,则可以作出决策At∈ωk。

3 基于Sec-(2D)2PCA的图像重建

根据第2节,执行一次(2D)2PCA后,原始图像A的重建图像为$\widehat{A}{}_1={\rm Z}{}_{1}CX{}_1^{{\rm T}}$;执行第二次(2D)2PCA后,剩余图像A′的重建图像为$\widehat{A}{}_2={\rm Z}{}_2{C}^{\prime }X{}_2^{{\rm T}}$。所以,由Sec-(2D)2PCA方法,原始图像A的重建图像为:

$A\approx \tilde{A}=\widehat{A}{}_1+\widehat{A}{}_2={\rm Z}{}_{1}CX{}_1^{{\rm T}}+{\rm Z}{}_2{C}^{\prime }X{}_2^{{\rm T}}$

我们分别用PCA、Sec-PCA、(2D)2PCA和Sec-(2D)2PCA方法做了图像重建的实验。实验结果如图3所示。

实验中,PCA和(2D)2PCA均按θ=95%的准则提取特征。对于Sec-PCA,第一阶的主元数为5,第二阶的主元数按θ=95%的准则提取,值为59;对于Sec-(2D)2PCA,第一阶的主元数为5,第二阶的主元数按θ=95%的准则提取,行、列方向分别为20、24。这里,Sec-PCA和Sec-(2D)2PCA的一阶主元数均按使识别率最大的原则确定。从图3可见,Sec-(2D)2PCA的重建图像质量明显优于其他的重建图像。

4 实 验

我们选择常用于测试光照变化条件下人脸识别问题的extended Yale人脸库[5],比较了四种方法的识别性能:即PCA、Sec-PCA、(2D)2PCA和本文的Sec-(2D)2PCA。

Extended Yale人脸库含38人,由Yale人脸库B[6]和Yale人脸库C构成。Yale人脸库B包含了10个人的5760单光源人脸图像。每一个人的人脸图像均在576种视角条件下(9姿态 64光照条件)获取的。根据光源方向偏离摄像头的主轴方向的角度,这个库被分成了四个子库:子库1(12°),子库2(25°),子库3(50°)及子库4(77°)。Yale人脸库C包含了28个人的16128单光源人脸图像。该库实际上是Yale人脸库B的扩展,每个人的图像是在和Yale人脸库B相同的方式下获取的。这样,整个extended Yale人脸库被分成4个子库,每个人在子库1、2、3和4中分别含有7、12、12和14幅图像。由于本文主要处理光照变化条件下的人脸识别问题,我们仅仅选择了4个子库中的1710幅正面人脸图像,其中4个子库依次含266、456、456和532幅图像。作为示例,图4给出了一个人在4个子库中的部分正面样本图像。可见,从子库1到子库4,光照变化越来越强烈。

实验中,子库1被选作训练样本库,其余的3个子库分别被选作测试库。我们在Matlab 6.5环境下编程,机器配置为:HP xw4300 Workstation,Intel Pentium 4 CPU 3.2G(双CPU),内存2.0G。软件编程过程:

训练阶段:

1) 设定r1(r1从1到10取值,逐次增1,c1=r1);

2) 读取子库1的图像为样本图像,按照2.2节的方法,学习到行、列方向最优投影矩阵X1(nr1)和Z1(mc1),提取样本图像的一阶特征矩阵C,并分别保存为MAT文件;

3) 计算余像,按θ=95%的准则学习到行、列方向最优投影矩阵X2(nr2)和Z2(mc2),提取余像的二阶特征矩阵C′,并分别保存为MAT文件。

识别阶段

1) 从子库2(3,4)中读取测试图像,并加载训练阶段步骤2所得的最优投影矩阵X1、Z1,对测试图像在行、列方向上投影,得一阶特征矩阵Ct;

2) 计算测试图像的余像,并加载训练阶段步骤3所得的最优投影矩阵X2、Z2,对余像在行、列方向上投影,得二阶特征矩阵C′t;

3) 加载训练阶段步骤2和3所得的一阶和二阶特征矩阵C、C′,根据上述所得的Ct和C′t,按2.2节的距离对测试图像分类(a2∈[0.5,1],初值为0.5,步长为0.05,逐次搜索最低错误识别率)。

Sec-(2D)2PCA的识别性能与参数r1、a2的组合有关(因c1=r1,a1=1-a2)。在实验中,对每一个确定的r1,逐次搜索a2,从而选出最佳的r1、a2组合,使得错误识别率最低。

四种方法的错误识别率比较见表1,运行时间的比较见表2。表1、表2中记号的意义如下:PCA(p):p个主元,按θ=95%选取;Sec-PCA(p1,p2;a1,a2):p1、p2分别为一阶、二阶主元特征。当p1确定后,p2按θ=95%选取;a1和a2为一阶和二阶特征的权重系数。(2D)2PCA(rc):r、c为行、列方向上的主元特征向量数,均按θ=95%选取;Sec-(2D)2PCA(r1,c1; r2,c2; a1,a2):r1、c1为行、列方向一阶主元特征数,r2、c2为行、列方向二阶主元特征数;a1和a2为一阶和二阶特征的权重系数。

实验表明,当r1增至某一值时,对某一a2,Sec-(2D)2PCA达到最低错误识别率;而后当r1增加时,错误识别率又开始升高。图5在extended Yale人脸库上,对于r1、a2不同的组合,给出了Sec-(2D)2PCA方法在子库3上错误识别率的变化。从图5可见,当r1增至5,a2=0.95时,错误识别率最低。

5 结 论

本文提出了一种新的特征提取和人脸识别方法:二阶(2D)2PCA。在Extended Yale人脸库上的实验表明:

1) 在光照变化较大的条件下,Sec-(2D)2PCA的平均错误识别率为22.5252%。与PCA、Sec-PCA以及(2D)2PCA相比(其平均错误识别率分别为59.4%、30.68%、58.3233%),其识别精度明显提高。说明在光照变化较大的情形,Sec-(2D)2PCA是一种有效的人脸图像的描述和识别方法。

2) Sec-(2D)2PCA的平均计算时间为11.2243秒。相对于PCA、Sec-PCA相比(其平均计算时间分别为14.378秒、28.9110秒),Sec-(2D)2PCA只需要较少的计算时间。这是因为Sec-(2D)2PCA是直接基于图像矩阵的特征提取,而不需要像PCA和Sec-PCA那样必须将图像矩阵转化为向量来处理,所以节省了大量的特征提取时间。

摘要：研究光照变化条件下的人脸识别问题。基于最近提出的二阶特征脸方法和(2D)2PCA方法,提出了二阶(2D)2PCA方法。该方法将(2D)2PCA技术分别应用到原始图像矩阵集和剩余图像矩阵集。在extended Yale人脸库上的实验表明,在光照变化条件下,二阶(2D)2PCA方法是一种有效的人脸识别方法。该方法与传统的特征脸、二阶特征脸方法和(2D)2PCA相比,具有更高的识别精度;且比特征脸和二阶特征脸方法节省计算时间。

关键词：人脸识别,特征脸,二阶特征脸,(2D)2PCA,光照变化

参考文献

[1] Belhumeur P N,Hespanha J P,Kriegman D J.Eigenfaces vsfisherfaces: class-specific linear projection[J].IEEE TransPattern AnalMachIntell,1997,19(7) : 711-720.

[2] Wang L,TanA T K.New proposal for face feature description[R].ISO/IEC JTC1/SC29/WG11/m5750,Noordwijkerhout,2000.

[3] Yang J,Zhang D,Aledjandro F,et al.Two-dimensional PCA: a new approach to appearance-based face representation and recognition[J].IEEE TransPattAnalMachIntel,2004,26(1) : 131-137.

[4]Zhang D Q,Zhou Z H.2D(PC)2A:Two-directional two-dimensionalPCA for efficient face representation and recognition[J].Neurocomput-ing,2004(25):1173-1181.

[5]Lee Q-C,Ho J,Kriegman D.Acquiring linear subspaces for face recog-nition under variable lightingIEEE TransPattern AnalMachIntel[J].2005(275):684-698.