刚体动力学（精选4篇）
刚体动力学 篇1
充液系统是一种特殊的力学系统,需要6个常微分方程和4个偏微分方程共同描述.它们分别是描述质心运动的动量方程、描述刚体转动的Euler方程、描述黏性不可压液体的Navier-Stokes方程,以及液体不可压缩条件下的连续方程.采用Hamilton-Ostrogrdskiy变分原理[1]或Jourdain原理[2],再结合Lagrange乘子方法,可以推导出充理想液体的刚-液耦合系统的动力学方程.考虑流体的黏性摩擦力并用偏应力张量表示之后,这种推导方法可以推广到充黏性流体的情况[3].本文从黏性流体的力学性质出发,将流体微团上单位体积流体的表面力合力(流体的内力)作为主动力,无需引入Lagrange乘子,利用Jourdain原理可以直接推导出充黏性液体系统的动力学方程.该方程可以退化为充理想液体的情况.
1 充液系统的力学模型
设系统由刚体和液体组成,如图1所示.分别用V'和V表示刚体和液体占据的空间区域,Γw,和Γf分别为湿壁面、非湿壁面和自由液面,液体是连续介质,用Euler方法描述,流场中各质点相对刚体的流速u(r,t)和压强p(r,t)都是位置坐标和时间的函数.为了描述充液系统的运动,我们建立两个坐标系,即惯性坐标系OoXYZ和刚体固连坐标系Oxyz.O相对Oo的矢径为ro,系统内的任意质点ml对Oo的矢径为
其中ρl为ml对O点的矢径.质点ml的速度和加速度为式中ω为刚体的角速度,为刚体的角加速度,,分别为质点ml相对刚体的速度和加速度.对于刚体上的质点,ρl为常数,相对速度和相对加速度均为零.
另外,由牛顿型不可压缩黏性液体的本构关系,液体域内任意点处的应力张量P的分量形式如下
μ为液体动力黏性系数.液体微团上受到的表面力的合力,与该点应力张量及作用面的法向矢量n相关.根据广义Gauss定理,单.位体积流体的表面力的合力为[4]
其中D为液体微团体积.这个合力正是质点ml处单位体积微团所受的液体内力.
2 变分原理结合Lagrange乘子方法建模过程简介
文献[1～3]以变分原理建立系统的动力学方程时,都通过Lagrange乘子将液体的不可压缩条件作为约束引入,各自的形式如下:
Hamilton-Ostrogrdskiy积分变分原理
Jourdain微分变分原理
其中T是系统的动能,Fl为任意质点所受的主动力,u为液体相对刚体的速度,p是Lagrange乘子.该乘子对应的约束是不可压缩条件,这样p与流体的动压强相同,因此动压强可看作系统的约束力[1].显然,这样的分析对充理想液体的情况是合适的,因为理想流体的应力张量各向同性且液体微团的变形与运动无关,即应力张量为
这样微团表面的受力恰好是液体的动压强p在微团表面的积分.
利用Gauss定理计算式(6)和(7)中含Lagrange乘子的积分项,最终可以得到充液刚体的Newton-Euler方程及腔内理想液体的欧拉方程[1,2].
若充液刚体系统中液体是不可压缩的黏性液体,还必须考虑与速度梯度有关的黏性摩擦力[3]
其中.将X作为主动
力写进式(6)或(7)的Fl项中,采用与充理想流体刚体系统的相同的推导方法,就能获得充液刚体的Newton-Euler方程及腔内流体的Navier-Stokes方程.
3 充液系统的动力学方程推导
我们将液体微团看作质点,液体看作质点系.每个质点所受的主动力Fl由外力和内力两部分组成,即Fl=.将Jourdain变分原理写成
将式(2)两边取等时变分,代入式(10),第1项展开为
将式(11)右端最后1项变成刚体区域和液体区域分别求和
式中ρ为液体的密度.将系统的动量记为Q,系统对O点的动量矩记为H,式(11)化简为
由上面的分析知,液体质点除了受到质量力f,还有表面力合力Y,展开式(10)的第2,3项得到
式(14)右端的前两项为系统内力的虚功率,注意到,系统任意点处的内力(包含湿壁面都成对出现,且方向相等大小相反,因此有,.将式(14)右端最后1项变成刚体区域和液体区域分别求和
式(15)第1项是液体的体积力(外力)的虚功率,第2项是自由液面处的面积力(外力)的虚功率,其中po是自由表面处的环境压强.第3项是液体质点ml处微团所受合力Y(内力)的虚功率.将外力主矢、对O点的主矩分别记为F,M,整理式(13)～(15),有
由于,δω,δu均是独立变分,可以导出以下的动力学方程
和自由液面动力学边界条件
其中式(19)是液体的动量方程,若代入式(3)以及牛顿型不可压缩黏性液体的应力张量表达式(4),就得到NavierStokes方程
其中是液体的运动学黏性系数,上式左端第3,4,5项之和,可作为流体计算方程的源项,写到式子的右边.
注意到,如果液体是理想流体,将方程(19)中的应力张量用式(8)代替,就可以得到液体的欧拉方程,相应地,
式(20)为理想流体的自由面动力学边界条件.
4 结论
本文在充液刚体动力学方程推导过程中,用液体应力张量表示质点ml处单位体积微团的表面力合力Y,其中包括单位体积液体微团的流体压强合力以及黏性偏应力张量的合力.在用Jourdain变分原理建立刚-液耦合系统动力学方程的过程中,无需引入Lagrange乘子,而将合力Y统一作为主动力处理,出现在系统内力的虚功率中,物理意义清晰.得到的动力学方程可以退化到理想液体情况.
摘要：针对充黏性液体的刚体，给出了一种应用Jourdain变分原理推导系统动力学方程的方法．将流体微团上所受的表面力合力视为质点系的主动力，出现在系统内力的虚功率中，代入系统的动力学普遍方程，从而获得充液刚体系统的Newton-Euler方程及腔内液体的Navier-Stokes方程．该方法不仅可以导出充黏性液体系统的动力学方程，还可以退化为充理想液体的系统动力学方程．
关键词：充液系统,黏性液体,液体应力张量
参考文献
[1]莫依舍夫HH，鲁苗采夫BB．充液刚体动力学．韩子鹏译．北京：宇航出版社，1992．1～26
[2]王照林，刘延柱．充液系统动力学．北京：科学出版社，2002．1～24
[3]刘延柱．刚体动力学理论与应用．上海：上海交通大学出版社，2006．181～183
[4]张兆顺，崔桂香．流体力学．北京：清华大学出版社，2001．367
刚体动力学 篇2
在压杆稳定性问题的教学过程中, 以杆发生弯曲变形后的构型为切入点, 学生容易产生疑惑:为什么压缩杆件会存在稳定性问题?为什么特征值对应的压力就是临界压力?为什么试验结果与理论分析差距较大?本文在理论力学的基础上, 以弹簧刚体模型为例, 通过稳定性分析的能量法引入小变形、大变形、考虑缺陷等情况下结构稳定性的几个基本概念.
1 小变形理论
分析在集中外力作用下的一端受圆柱销和扭转弹簧并联约束的刚性杆, 理论模型如图1所示.刚性杆的长度为L, 扭转弹簧刚度系数为k, 杆水平时弹簧处于不受力状态, 在杆的自由端施加水平向左的集中力P, 分析系统的稳定性.
假设系统处于变形后的构型, 如图1所示.
则系统总能量为
式中, U为系统的变形势能, W为外力做的功.对于这里的问题
式中, θ为刚性杆转过的角度.将式 (2) , 式 (3) 代入式 (1) 可得
在稳定性分析的能量法中, 平衡对应于 稳定性分析则依赖于T对θ的高阶导数.于是可得系统平衡方程
对于小变形问题
当轴向压力P增加到一定程度时, 杆件直线平衡状态开始失去稳定产生变形, 这个力具有临界的性质, 因此称为临界力.将式 (6) 代入式 (5) 可得临界力
系统在临界力的作用下可处于变形后和变形前的平衡位置, 这取决于临界力的稳定性.在小变形状态, 对系统总能量函数T进一步求导数可判定系统的稳定性, 由于
所以如果 系统在临界力作用下稳定;如果 系统在临界力作用下失稳;如果 即P=Pcr, 系统在临界力作用下稳定性不确定.也就是分析在某种变形模式下, 结构内部变形势能和外力做功的大小, 如果外力功小于变形势能则结构稳定, 反之则不稳定, 如果两者相等, 需进一步判断.在小变形假设下, 临界力与系统状态变量无关, 所以难以对后屈曲形态进行判别, 无法得到明确的平衡路径分叉图.
2 大变形分析
为了对后屈曲形态进行判别, 放弃小变形假设, 直接通过求解式 (5) 得到系统平衡时的集中力
通过式 (9) 可直接制作平衡路径的分叉图, 进行大变形分析.平衡路径在Cr点 (参见图2) 出现了4个分支, Cr点称为平衡路径的分叉点, 分叉点对应的压力则为临界力.当构件承受的压力大于临界力后所表现出的稳定性能称为后屈曲稳定性能, 故此时平衡路径对应的稳定性分析称为非线性后屈曲平衡路径稳定性研究, 其研究通常要考虑非线性因素的影响.为了研究方便, 通常将工作压力除以临界力称为无量纲屈曲因子λ=P/Pcr[1,4], 其随θ的变化规律如图2所示.
式 (9) 除以式 (7) 可得
当θ趋向于0时λ趋向于1, 并且随着|θ|的增大, λ不断变大, 呈明显非线性规律.下面通过考虑几何非线性——大变形的影响来分析式 (7) 临界载荷的稳定性, 式 (4) 直接对θ求二阶导数
在后屈曲平衡路径上, 式 (9) 成立, 将其代入式 (11) 可得
分析式 (12) 发现对于所有不为0的θ, 其值永远大于0.所以θ=0时, 平衡路径是稳定的.当θ=0时, 系统平衡路径通过临界点, 这时的稳定性研究称为临界点稳定性分析, 并且总能量的二阶导数 所以不能通过二阶导数直接判定稳定性.在θ=0的情况下, 对总能量函数求更高阶导数, 寻找下一个不为0的项
通过Taylor级数, 取四阶项可近似表示系统在θ=0时的总能量
由于下一个非零项为偶次的, 并且大于0, 所以系统在θ=0时是稳定的, 也就是说压杆的临界力表达式 (7) 是稳定的, 即临界点是稳定的.事实上分析式 (14) 发现, θ在0位置的任何变化都会引起总能量T的增加, 说明变形势能对总能量的影响大于外力做的功, 即系统还有抵抗外力增加的能力, 所以系统在该位置是稳定的, 对应的该类稳定问题对应的平衡路径分叉称为稳定分叉.分析图2发现, 对于λ<1的情况, 由扰动导致系统偏离θ=0时, 当扰动消除后系统会自动恢复到θ=0的平衡位置, 所以这里的平衡路径是稳定的;对于λ>1的情况, θ=0时系统平衡需要的能量大于θ=0的位置, 所以θ≠0为稳定平衡路径.一般情况下, 系统不会停留在λ>1, θ=0的位置, 所以图上用虚线表示.
至此可以判定扭转弹簧作用下的刚性压杆屈曲前、临界点以及后屈曲过程中的平衡状态都是稳定的.
3 缺陷的影响
实际结构总会存在各种各样的缺陷, 为了引入缺陷的影响, 仍以扭转弹簧作用下的刚性压杆为例, 假设初始状态下, 杆与水平线之间存在夹角α, 如图3所示.系统总能量为
式 (15) 对θ求一阶导数, 并令其等于0可得平衡路径表达式
于是可以直接通过式 (16) 除以式 (7) 可得 无量纲屈曲因子λ随转角θ和缺陷α的变化规律如图4所示.
与图2相比, 考虑缺陷 (α=0) 的图4显然已经不存在明显的λ-θ路径分叉;但是随着α越来越小, (λ-θ) 路径趋向于不考虑缺陷系统的 (λ-θ) 路径;缺陷会明显影响临界载荷的大小.对于实际结构, 缺陷不可避免, 所以实验结果很难得到无缺陷结构的分叉点和临界力.
下面分析考虑缺陷时平衡路径的稳定性, 式 (15) 对θ求二阶导数可得
将式 (16) 代入式 (17) 可得
4 不稳定实例
考虑一端铰接, 另一端连接横向线弹簧的刚性杆, 该弹簧的另一端可在轴向自由移动, 横向分别固定在地面和刚性杆上, 如图5所示.
直接通过大变形理论分析, 系统总能量为
令式 (19) 对θ的一阶导数为0可得平衡载荷的表达式
式 (19) 确定的 (λ-θ) 曲线如图6所示.当θ=0时通过式 (20) 可得小变形状态下的临界力
观察图6发现, (λ-θ) 曲线在θ=0处, 对θ增加任意的小扰动都会使得屈曲稳定因子下降, 从而导致系统不稳定.事实上通过能量法容易判别, 随着P的增加, 在Pcr以后系统都是不稳定的, 该类稳定问题对应的平衡路径分叉称之为不稳定分叉.
式 (19) 对θ的二阶求导, 并令式 (20) 成立可得
当θ=0时, 式 (22) 小于0, 系统不稳定.下边分析θ=0的稳定性问题, 令式 (20) 成立, 并且θ=0, 寻求第一个使得 都不为0的i, 发现
所以θ=0的临界点的平衡状态是不稳定的, 这证实了式 (21) 确定的临界力是不稳定的.就是说在图6上, λ<1, θ=0平衡路径是稳定的;θ≠0平衡路径是不稳定的;特殊的分叉点λ=1, θ=0也是不稳定的.
如同第3小节一样考虑缺陷的影响, 此时系统总能量为
令式 (24) 对θ的一阶导数为0可得平衡载荷的表达式
令式 (25) 除以式 (21) 可得
令式 (25) 对θ的一阶导数为0可得各个路径上最大的压力Pmax=k L cos3θ对应的θ值
通过式 (26) 和式 (25) 确定的 (λ-θ) 曲线如图7所示, 虚线表示式 (27) 确定的最大压力.在式 (25) 成立的情况下, 式 (24) 对θ的二阶导数
容易证实, 当P
5 结论
本文针对简单的结构稳定性问题, 以刚性杆的静力学为基础, 应用能量法引出了结构稳定性分析的几个重要基本概念:临界力、平衡路径、平衡路径分叉点、后屈曲、缺陷、大变形、无量纲屈曲因子、极值点屈曲等, 由于避免了复杂的弹性力学公式, 使得问题概念清晰, 易于理解.
但稳定性问题由于研究对象和载荷的复杂性, 通常难以用本文的方法来描述, 实际上材料力学与结构力学中的研究对象是变形体, 为无穷多自由度系统, 其对应的稳定性问题为泛函的驻值问题, 而单自由度系统的稳定性问题则是函数的极值问题, 虽然通过单自由度系统能够引出稳定性问题的一些概念, 但稳定性完整理论分析要复杂得多.
参考文献
[1] Barber JR.Intermediate Mechanics of Materials.Oxford:Oxford University Press Inc, 2011
[2] Boresi AP, Schmidt RJ.Advances Mechanics of Materials.Laramie:John Wiley&Sons, 2003
[3] Simitses GJ, Hodges DH.Fundamentals of Structural Sta-bility.New York:Elsevier Inc, 2006
[4]范钦珊, 郭光林.工程力学2.北京:高等教育出版社, 2011
[5]龙驭球, 包世华.结构力学教程II.北京:高等教育出版社, 2001
[6]陈骥.钢结构稳定理论与设计.第4版.北京:科学出版社, 2008
刚体的平面运动和旋轮线 篇3
令, 上式变为
此式表示一个旋轮线.对式 (1) 求时间的导数, 得到M的速度投影
若刚体逆时针转动, 则式 (1) , (2) 变为
众所周知, 旋轮线可以认为是由沿一条直线做纯滚动的圆轮上 (或拓展部分) 的一点刻画而成, 但是式 (1) , (3) 表示的旋轮线具体形状显然随V, ω, l的不同而不同.因此, 以下讨论平面运动刚体上的质点在不同情况下的旋轮线轨迹是由怎样的圆轮作怎样的纯滚动形成的[1], 以期对质点的运动有更直观的理解.
1 旋轮线的形成过程
首先讨论式 (1) 所示的情形, 并且只研究ωyC>V, 即yC>R的情形.
第1种情况:运动的初始条件有lω>V, 以致l>R.此时, 旋轮线形成过程如图2所示:圆轮半径为R, 圆轮拓展部分上与轮心距离为l处有一点M, 起初在Oy轴上且在最高点;Ox轴上方有DE直线y=yC-R, 圆轮在DE直线上方沿着DE作顺时针方向匀角速纯滚动时, M点的轨迹方程即为式 (1) ;N是圆轮的瞬心, 由R=V/ω可知, N也是刚体的瞬心.这种情况下:曲线的斜率总是存在的, 曲线上没有“折点”, 是处处滑顺的曲线;当ωt= (2n-1) π (n=1, 2, ···) 时, vx<0, vy=0, 曲线斜率k=vy/vx=0, 表示质点在这些位置处作运动方向与Ox轴负向一致的水平运动, 运动方向与圆轮质心运动方向相反.
第2种情况:lω=V, 以致l=R.此时, 旋轮线形成过程如图3所示:圆轮半径为R, 圆轮边缘上有一点M, 起初在Oy轴上且在轮的最高点;Ox轴上方有DE直线y=yC-R, 圆轮在DE直线的上方沿着DE作顺时针方向的匀角速纯滚动时, M点的轨迹方程即为式 (1) ;N是圆轮的瞬心, 也是刚体的瞬心.这种情况下, vx≥0;但是, 当ωt= (2n-1) π (n=1, 2, ···) 时, vx=0, vy=0, 故曲线的斜率不存在且带电粒子的运动方向突变, 表示质点的运动轨迹在这些位置出现“拐点”.
第3种情况:lω
对于ωyC
再讨论式 (3) 所示情形, 并且也只研究ωyC>V, 即yC>R的情形.
第1种情况:lω>V, 以致l>R.此时, 旋轮线形成过程如图5所示:圆轮半径为R, 其拓展部分上与轮心距离为l (l>R) 处有一点M, 起初在Oy轴上且在最高点;Ox轴的上方有DE直线y=yC+R, 圆轮在DE下方沿DE作逆时针方向的匀角速纯滚动时, M点轨迹的方程即为式 (3) ;N是圆轮的瞬心, 也是刚体的瞬心.这种情况下:曲线斜率总是存在, 曲线无“拐点”, 处处光滑;但是, 当ωt=2nπ (n=0, 1, 2, ···) 时, 恒有vx<0, 曲线斜率k=0, 表示质点在这些位置作运动方向与Ox轴负向一致的水平运动, 运动方向与轮心的运动方向相反.
第2种情况:lω=V, 以致l=R.形成这种旋轮线的过程如图6所示:圆轮的半径为R, 边缘上与轮心距离为l (=R) 处有一点M, 起初在Oy轴上且在最高点;Ox轴的上方有DE直线y=yC+R, 圆轮在DE下方沿DE轴逆时针匀角速纯滚动时, M点的轨迹即为式 (3) , M点的轨迹是普通旋轮线;N是圆轮的瞬心, 也是刚体的瞬心.这种情况下恒有vx≥0;但是, 当ωt=2nπ (n=0, 1, 2, ···) 时, vx=0, vy=0, 曲线斜率不存在且运动方向突变, 质点轨迹在此处出现“折点”.
第3种情况:lω
2 实例
如图8所示, 光滑水平面上静置质量为M, 长度为2L的均质细直杆, 质量为m的质点以速度v0沿水平面垂直射入直杆一端.求:碰撞后质点的绝对运动方程.
研究整个系统.以质点的碰前速度方向为速度的正方向, 竖直向下为动量矩的正方向.碰撞过程中, 系统的水平方向动量守恒, 其质心C的水平速度不变, 有
并且, 质点与 C 之间和 C 与直杆质心 C2 (如图9所示) 之间的距离均不变, 分别为
碰撞前:质点绕C点的动量矩为L1=lm (v0-V) , 直杆绕C点的动量矩为L2=-a M (0-V) .碰撞之后:设系统绕通过质心C且垂直于细杆的轴的转动惯量为IC, 则系统绕C点的动量矩为L3=ICω.由动量守恒定律, 有
即
将-代入, 式 (5) , (6) 代入上式消去a, l, 得到
以光滑水平面上与直杆中心C2的初位置重合的固定点O为原点建立定系Oxy, 如图9, 图10所示, Ox轴平行于质点的初速度方向.由式 (3) , 质点的绝对运动方程为
其中
显然, 这表示一个旋轮线, 如图10所示:它是由半径为R的圆轮沿DE直线y=a-R做纯滚动时, 其内部或延拓部分上、距离轮心 C 为 l 的点所刻画的曲线; N是圆轮的瞬心, 也是系统的瞬心.
摘要：讨论了作平面运动的薄板上一个质点在不同条件下的旋轮线轨迹, 并讨论了这些不同形状的旋轮线是怎样由圆轮的纯滚动形成的.
关键词：平面运动,质心,旋轮线,曲率半径
参考文献
刚体动力学 篇4
本文来源于与故宫科技部合作研究青铜鼎碎片复原的项目。青铜鼎破碎成大小100多片, 人工匹配拼接消耗大量时间和精力, 于是我们利用计算机辅助拼接来虚拟修复青铜鼎。与此相关的工作有:Wolfson[1]提出先用多边形近似方法得到碎片2D边界曲线的特征串, 然后通过几何杂凑法寻找最长公共子串。这种方法不能给出曲线匹配的准确信息, 也不能完全排除重叠的情况。Ucoluk[2]讨论了3D物体的自动重建问题, 在局部匹配中采用了穷举方法, 时间复杂性非常高, 而且没有用实际数据进行测试。Kong[3]研究了2D和3D拼图问题, 采用动态规划法进行曲线匹配, 但是没有考虑碎片之间不完全匹配的情况 (例如物体在破碎过程中可能会丢失一些非常小的碎片) 。杨洛斌[4]等就形状匹配在文物复原方面的研究, 其主要思想为以其中任意一个碎片为模板和其他碎片进行一一匹配, 选择其中匹配曲线段最长的碎片作为匹配碎片并进行曲线拼接和曲面拼接, 将拼接后的碎片作为新的模板进行新一轮匹配, 该工作以最长匹配段作为判定。还有一系列的文章[5,6,7,8] , 使用弹性模型进行曲线匹配, 这种方法计算复杂度很大, 尤其在碎片比较多的情况下, 几乎就是不可解。本文使用的方法是采用分级分段的方法进行匹配计算, 使用分级可以加快计算速度, 利用分段不仅可以降低采样点的个数, 还可以利用相临的段进行验证, 提高匹配正确的概率。
1数据预处理
1.1扫描数据
利用三维激光扫描设备 (FastScan) 采集破碎文物的表面数据信息 (点云数据) , 对扫描碎片只扫描其中一面, 并且所有碎片都只扫描同一侧面的面。 (根据具体情况选取内面还是外面, 根据具体情况, 这里取外面, 因为青铜器外面比较光滑, 更容易拼接计算) 。通过扫描, 获得到的点云三维数据, 经三角网格剖分后生成其对应的三角网格曲面模型, 如图1所示。
1.2预处理碎片数据
三维网格曲面模型经过适当简化后, 即可进行后续的特征轮廓线提取、匹配等处理。在使用三维激光扫描仪进行数据采集的过程中, 对空间形体的最高分辨精度可以达到0.1mm, 因设备所引起的比如裂缝、噪声以及离散网格的非流形问题可以通过数据简化来解决 (该功能由设备的软件处理平台提供) 。但是扫描的模型中带有的断裂面的数据并不能消除掉, 我们需要把断裂面数据去掉, 只获得轮廓的信息。这里使用的消除断裂面并获得轮廓线的算法是参考文献[9]中的方法。该方法就是利用曲率, 提取外轮廓线和内轮廓线, 然后再根据物体的破碎特征得到真正的外轮廓线, 该方法与物体破碎程度有关, 所以这里加入了人工辅助修正。图2所示为经过预处理后提取的轮廓线。
2基于轮廓线的匹配技术
2.1轮廓线分段
如图3所示, 根据取得的轮廓线的各个点的曲率, 按照曲率变化, 把整体的轮廓线转化为关键点的集合, 这些关键点就描述了这个整体轮廓。在按照关键点中的拐点 (就是那些变化率大的关键点) 划分成不同的段。 (图3中分成了A、B、C三个段) 也就是以曲率变化大的那些点作为关键点对整体轮廓进行划分, 分成不同的段, 这样就得到所有原始点的关键点的集合, 各个段的集合以及段的关系。其中段和段的关系是使用段与段之间的关键点来描述的。如果使用曲线拟合重采样获得点的集合来描述这个曲线, 在重采样的过程中有可能丢失掉那些特征比较明显的点, 比如拐点。这里的分段分两级, 第一级为粗线条分, 就是按照曲率变化比较大的区域划分;第二级, 在已经划分好的段, 再次使用同样的方法对已划分的段再次进行划分, 新划分所使用的曲率变化比上一级曲率变化要小。
2.2轮廓线匹配技术
2.2.1 轮廓线的匹配标准
(1) 通过段匹配计算, 如果两个段之间存在的重叠区域很小并且段匹配的匹配长度最长, 那么匹配成功。
(2) 如果一个段匹配成功, 根据段与段之间的关系, 计算与匹配的段相互连接的段, 如果这个相邻接的段也匹配成功, 那么这个匹配是成功的。
各个碎片的轮廓线匹配拼接过程如下:
1) 选择一个碎片的一个段 (可以随机选择, 也可以指定) , 以它作为基本匹配段, 简称基段;
2) 在其他碎片的段中寻找可以与之匹配的段, 简称匹配段 (所使用的是匹配标准A) ;
3) 如果在其他碎片中找到一个与之匹配的段, 记录下他的匹配优先级, 匹配优先级按照匹配度排序;
4) 选择匹配优先级最高的进行拼接, 拼接完成后, 必然引起匹配段所在的轮廓线中的其他段的位置变化, 根据这些变化, 计算与基段相连接的其他的段与这个匹配段所在的碎片的其他段的匹配情况, 如果有新的匹配产生, 那么, 整体匹配成功。 (使用的是匹配标准B) 如果没有新的匹配产生, 可以认为是可能的匹配, 作为备选, 根据人工判断是否匹配成功;
5) 以上是第一级匹配, 如果匹配成功, 那么进入第二级匹配。第二级匹配与第一级匹配算法是一致的, 唯一的不同是不需要再次进行匹配验证。第二级匹配的目的是使匹配度更加精细, 拼接结果更理想。
2.2.2 轮廓线各分段匹配计算
如图4所示, 有2条要匹配的段, 假设上面的段为段A, 下面的段为段A′, 那么可以获得他们的原始特征点。图5为这2个段的原始特征点。
把各个段的原始特征点进行坐标无关化处理, 也就是各个取得采样点的坐标无关信息。坐标无关化使用的方法为参考文献[4]中的方法并作出改进。取起始点作为基点设坐标矩阵为Tb, 那么与它相邻接的特征点, 假设坐标矩阵为Tl, 那么它们之间的关系矩阵设为Tlb, 这个矩阵是一个与坐标旋转、平移无关的量。图6为。特征点和相应的坐标无关的方向矢量。
现在假设二个段的采样点集合为 { VA| VA[i], i=1, 2, … }和{ VB| VB[j], j=1, 2, … }, 那么他们的匹配过程如下:
(1) 以VA[i]为原点, 按照方向矢量建立坐标系, 并还原得到曲线;VB[j]的点为原点, 在VA[i]建立的坐标系, 还原得到曲线。
(2) 计算二个曲线之间的关系, 如果满足匹配规则就匹配成功, 不满足i++;或者j++;返回到步骤 (1) 。
在所有满足匹配规则的成功匹配中, 那个最长的匹配就是要求的结果, 相应的伪代码为:
match () 函数是计算以VA[i]、VB[j]为原点的两个曲线的匹配长度。
length () 函数为计算曲线含有的原始特征点的个数。
int m=length (VA) ;
int n=length (VB) ;
for (int i=0;i<m;i++)
for (int j=0;j<n;j++) {
if (Max<match (VA[i], VB[j]) )
Max= match (VA[i], VB[j]) ;}
match函数的原理:
以VA[0]点为原点, 那么VA[1]就是该点的下一个点的描述信息, 只要做VA[1]+ VA[0], 就是新坐标中的点, 同样, 可以得到以VB[0]为原点的曲线。
假设VA[i]为原点和VB[j]为匹配起点, 那么VA[i]与VA[i+1]就存在一个凸包, 把这个曲线包围起来, 如图7所示凸包的最大半径为R, 如果VB[j]与VB[j+1]在这个凸包中, 那么i=i+1, j=j+1, 继续计算, 测试新的点是否在凸包里, 如果VB[j]与VB[j+1]不在凸包中, 那么, 计算从起点到终点的曲线长度, 这个长度就是匹配的优先级。
3碎片拼接
如果分段VA和分段VB是匹配的, 那么就需要把VA和VB所在的曲面进行拼接。
VA和VB是坐标无关的, 所以, 首先把他们放到坐标系中, 建立坐标系, 那么VA和VB所在的碎片的拼接就是按照VA和VB中的特征关键点重合的规律, 通过平移和旋转来实现。
拼接过程如下:
首先, 建立局部坐标系。假设破碎曲面中的基段为VA, 匹配段为VB, 根据分段匹配计算的结果, VA, VB的匹配结果, 至少有二个点是重合或者是接近重合的, 假设这两个点为 (VA[i], VB[j]) , (VA[m], VB[n]) , 以起点VA[i]为原点, VA[i]和VA[m]之间的连线为x轴, 向x轴作垂线, 那么这个垂线的方向就是y轴的方向, 这样就得到了x和y轴, x轴, y轴的向量叉乘就是z轴。如图8所示。
因为所记录的各个特征点信息是与坐标无关的, 所以, 只要把特征无关向量矩阵代入所建立的坐标系, 就得到了拼接结果。剩下要做的就是曲线合并, 它的做法为, 找到重合的点, 把重合区域合并, 并把两个碎片拼接起来, 获得一个新的碎片。具体拼接算法如下:
假设有不同的两个曲线pi、pj, 对于pi、pj, 它的局部坐标系是不一致的。当两个曲面的轮廓曲线中有一段曲线段匹配, 那么这两个曲面将被拼合成为一个曲面。由于匹配曲线段近似认为是碎片的碎裂曲线, 属于共享边界, 所以会将两个碎片曲面都转换到局部坐标系中。我们根据每一个pi、pj来确定一个局部坐标系。计算可以得到两个局部坐标系的变换公式, 而在理论上匹配曲线段应该在一个局部坐标系中是一致的, 所以实际上这个坐标变换公式中的旋转矩阵、位移向量就是我们所要得到的。
为了实现两个曲面、曲线的拼合, 将世界坐标系中的每一个碎片曲面图形通过坐标旋转、平移变换使得两个曲面拼接在一块, 两个曲面中的轮廓曲线中的匹配曲线段重合。
假设两个曲面SA、SB, 那么它们的局部坐标系分别表示为:
IA: { (o′A, e′A1, e′A2, e′A3) }
IB: { (o′B, e′B1, e′B2, e′B3) }
则根据局部坐标系和世界坐标系的点坐标变换公式:
AA、AB分别IA、IB坐标向世界坐标系的过渡矩阵。从这一公式可以得到IA、IB中的点的坐标变换关系:
设定旋转矩阵R和平移向量T:
那么:
这里的rij (i, j=1, 2, 3) 是R中的元素, ti (i=1, 2, 3) 是T中第i个分量。
实际工作中, 我们在作曲面SA、SB的拼合时, 计算出旋转矩阵R和平移向量T, 然后将对曲面SB加以旋转、平移, 使得两个曲面的轮廓曲线匹配段合一, 从而实现了两个曲面的拼合, 曲线的拼接, 如图9所示。最后拼接结果如图10所示。
4结论
本文提出了利用轮廓线分段进行匹配的方法。在分段的过程中, 使用曲率比较大的点作为特征点, 并作为分段的依据, 把一个碎片的轮廓分解成多个段, 再利用不同碎片的轮廓段进行匹配计算, 然后按照匹配规则获得不同优先级的匹配方案, 进行匹配, 然后根据匹配结果进行拼接计算。利用这个方法, 解决了破碎青铜器拼接的问题。
参考文献
[1]Wolfson Hj.On curve matching[J].IEEE Transactions on Pattern A-nalysis and Machine Intelligence, 1990, 12 (5) :483-489.
[2]Ucoluk G, Toroslu I H.Automatic reconstruction of broken 3d surfaceobjects[J].Computers and Graphics, 1999, 23 (4) :573-582.
[3]Kong Weixi, Kimia B B.On solving 2d and 3d puzzles using curvematching[J].Proceedings of the IEEE Conference on Computer Visionand Pattern Recognition (CVPR) , 2001, 22 (5) , 583-590.
[4]杨洛斌.形状匹配技术在文物复原中的研究与应用[J].西北大学学报:自然科学版, 2004, 2 (2) .
[5]Leitao HC G, Stolfi J.A multi-scale method for the reassembly of frag-mented objects[J].Proc.BMVC’2000, September 2000.
[6]Leitao HC G, Stolfi J.A multiscale method for the reassembly of two-dimensional fragmented objects[J].IEEE Transactions on Pattern A-nalysis and Machine Intelligence, 2002, 24:1239-1251.
[7]Leutwyler K.Solving a digital jigsaw puzzle[J].Scientific American, June 2001.
[8]Reid S K.Marble fragments from tablets in the petra small temple[J].Technical Report Publication Pending, Department of Archaeology, Brown University, 2001-2002.