变形算法

变形算法（精选7篇）

变形算法 篇1

0 引 言

随着网络技术的迅速发展, 逐渐扩张的逆向工程技术使不法用户可以通过静态分析与动态跟踪来分析软件的核心算法和重要数据, 严重侵害了开发者的知识产权[1]。因此, 软件保护越来越为重要, 而代码变形是软件保护的一个重要手段[2]。压缩变形是代码变形的重要技术之一, 经普通代码变形技术作用后的程序体积会由于插入了大量的指令代码而远超过原始文件[3], 空间资源的浪费会对软件和系统造成严重影响, 所以对变形文件进行压缩变得尤为重要。

由于可执行文件, 即PE ( Portable Executable) 文件的特定格式[4], 以及代码变形的一些特殊要求, 如文件的可执行性[5]与保密性, 一般的数据压缩方法在代码变形中并不能完全适用。目前的代码变形压缩技术通常只对PE文件中的某一部分采用现有的开源数据压缩算法。例如游程编码 ( RLE) 、哈夫曼编码、Zlib、a PLib、LZMA等算法在移植应用过程中并没有对压缩算法进行修改, 虽然拥有较高的压缩效果, 但保密性低且易被破解, 达不到软件保护的目的[6,7,8,9]。本文分析了LZMA算法在代码变形中的应用方式, 针对其中的不足, 提出进一步的增强方法, 使其在不改变文件可执行性的前提下, 更安全、更高效地完成PE文件的压缩。

1 代码变形的压缩机制

代码变形技术是面向PE文件的一种代码混淆技术, 其中的文件压缩机制与一般的数据压缩机制不同, 在代码变形领域, 经过压缩的PE文件仍为可执行文件。若按照传统的数据压缩策略, 根据Windows操作系统及其装载、执行机制的复杂性, 不加区分地对PE文件中的所有数据进行压缩, 极有可能导致压缩后的文件无法正常执行[10]。这表明不能将PE文件作为一个整体进行压缩, 而需分块处理, 最后添加解压缩代码并修改文件头, 生成一个新的PE文件。根据PE文件结构[4]可知,

PE = { DOSheader, PEheader, Section table, Section, Other} , 其中Section = { . data, . text, . idata, . reloc, . tls, . edata, . bss, . rda-ta, . rsrc} 对PE文件各部分的压缩策略进行分析如下。

( 1) 可压缩数据

. data, . text, . idata, . reloc均为可压缩的Section。

( 2) 不可压缩数据

{ DOS header, PE header, Section table, Other} 均在系统加载文件时使用, 不可压缩。. tls块和线程局部存储有关。系统加载PE文件时若. tls块处于压缩状态, 系统无法正确初始化线程局部变量。. edata块存放引出数据与函数的列表, 以供其他模块引用。操作系统可能在不装载PE文件的情况下使用引出信息, 因此. tls块和. edata块均不可压缩。

. rdata块可用于多处, 包括存放TLS目录、调试目录、OLE程序的GUID、DEF文件的描述信息等, 由于系统加载PE文件时需用到TLS目录与调试信息, 因此对. rdata块也采取不压缩策略。

. bss块最初用于在进程空间中保留一块存放未初始化的静态与全局变量的地址, 但现在使用较少, 一般情况下, PE文件中并无. bss块的原始数据, . bss块表中的Pointer ToRaw Data字段总是0。对. bss块亦采取不压缩策略。

( 3) 特殊数据

在Section中最特殊的是. rsrc块, 它用于存放程序资源。通常资源信息只有在PE文件运行时才会被用到, 但有些特殊类型的资源, 在PE文件没有运行的情况下, 仍被系统读取和使用。此类资源包括Version information ( 版本信息) 、Icon ( 图标) 、Group Icon ( 组图标) 等。因此将资源块分为两部分, 对可压缩的部分进行压缩, 对不可压缩的部分则保持原状。

令. rsrc块中可压缩部分集合为Crsrc, PE文件中可压缩数据集合为C, 不可压缩数据集合为NC, 则:

2 基于增强算法的代码变形模型

根据代码变形的压缩机制分析及其领域的特殊要求, 可对LZMA算法做进一步的增强处理, 本文提出的一种基于增强LZ-MA算法的代码变形系统模型如图1所示。

非压缩变形使用的变形技术包括插入跳转指令、插入无效指令、等效指令替换、指令乱序、反调试等[3]。

压缩变形模块采用的是本文提出的增强LZMA算法实现。包括加密、并行编码以及加壳。并行编码采用了分块并行策略, 增强的字典搜索策略, 以及增强的输出元组策略。

3 代码变形的增强 LZMA 算法

3. 1 加 密

传统LZMA算法开源、保密性较低, 无法达到软件保护的目的, 所以可在压缩过程中加入加密处理, 以提高压缩算法的安全性。

1) 密钥

密钥是包含4个32位二进制串的数组Key[4]。令BTime为当前系统时间的32位二进制串, BName为当前文件名的二进制串 ( 超过32位低位截断, 少于32位低位补零) , 可根据如式 ( 1) 所示规则得到一个32位BKey。

将BKey串按8位划分得 到4个子串, 分别是BKey1、BKey2、BKey3、BKey4, 则密钥Key[4]各成员如式 ( 2) :

2) 加密策略

代码变形中的加密与压缩类似, 均不可影响PE文件的正常执行, 故加密也需分块执行, 仅加密可加密的部分。令输入数据流为In[], 长度为l, 为向下取整) , 那么输出流Out[]可用式 ( 3) 计算, 将Out[]循环右移l mod32位得到密文。

加密后仅将BTime存入文件中, 在加密算法未知的情况下解密难度较大, 安全性较高。同时该加密算法均为位运算, 运算十分快速, 对压缩时间影响较低。

3. 2 分块并行

按传统的LZMA串行策略, 对PE文件中可压缩数据块 ( 总数为N) 进行编码的时间Ttotal为:

式 ( 4) 中TI ( I = 1, 2, …, N) 表示处理第I块可压缩数据的时间。为了提高算法的压缩时效, 本文按照PE文件可压缩数据集合C分块并行。此时的并行处理时间T'total为:

式 ( 5) 中TMAX为处理PE文件中可压缩数据块的最大时间, T'I ( I= 1, 2, …, N) 表示在并行策略中处理第I块可压缩数据的时间, 包括数据编码时间与并行线程间通信时间。可知T'total< Ttotal, 故采用分块并行策略可以提高LZMA算法压缩PE文件的效率。

3. 3 增强的字典搜索策略

LZMA算法的搜索缓存很大, 目前定位搜索符号的普遍做法是利用两个临近字节的散列值, 但有时会造成散列数组空间的浪费, 搜索时间也不理想。本文设计一种策略, 在遍历待编码区域的同时, 即可完成对字典的搜索查找以及字典的扩充, 不会造成空间的浪费, 压缩时效也较高。

1) 字典的存储结构

由于PE文件是二进制程序文件, 故字典的存储可采用基数树结构, 该基数树的结构特点为:

( 1) 每个节点包含一个三元组 < 索引, 权值, 二进制串 > , 即 < Index, Weight, String > 。Index项与插入字典的顺序有关, 将当前最大Index记作IMAX, 为避免与输出字符混淆, Index最小为2。初始化时树中有3个节点, 根节点Root为 < NULL, - 1, NULL > , 左子节点 < 2, 0, '0' > , 右子节点 < 3, 1, '1' > 。

( 2) 树中每个节点的String都是其左右子树String的前缀。令当前节点的String为SFATHER, 其左子节点的String为SLEFT, 其右子节点的String为SRHIGHT, SOTHER与S'OTHER表示剩余子串, 则:

( 3) 权值计算:

等比数列求和公式 ( 首项a1, 项数为m, 公比为q, k为整数) 。

令二进制串位数为M, 最低位默认为第零位, 根据式 ( 4) 计算出权值w为:

令当前遍历到第i位的子串权值为wi, 根据式 ( 8) , 可推导出i + 1位子串权值wi + 1为:

令wf为父节点的权值, 其左右节点的权值分别为wl和wr, 该树中节点权值必须符合wf< wl< wr。

2) 字典搜索算法

根据字典结构, 本文提出一种基于该结构的字典搜索算法, 算法描述为:

3. 4 增强的输出元组

传统LZMA算法有两种输出, 分别是“长度”, “位置”, 以及“下一个符号”, 这两种输出均经过区间编码[6]。但如果完全匹配某一项时, “长度”便跟字典树中该项二进制串长度相同, 故可简化完全匹配时的编码, 只需输出经区间编码的“位置”即可。这样可以在一定程度上减少编码输出, 提高压缩比。区间编码的编码模式如下。

对于任意的一个可用基数b的w次方表示的整数, 设有一个整数区间 ( i| L≤i≤H) , 令fs为符号S的频度, 符号值小于S的所有符号的频度总和为符号S的累积频度, 记作Fs。则:

可根据符号频度、累积频度以及总频度, 计算出该符号在区间[L, H]上的映射区间[L', H'] ( R'表示映射区间[L', H']的范围, N为所有符号的总频度) 。

编码后输出一个数值V ( 属于[L', H']) , 若V可分解为V =V'bn, 则V'为数据的区间编码。

3. 5 加 壳

为保证压缩变形后PE文件的可执行性, 数据压缩后必须进行加壳, 即添加解压解密代码。在PE文件加载时壳优先获得控制权, 先对文件解压解密后才能正常执行。本文采取加壳策略的模式是将各压缩块统一存储于生成目标PE文件的. data块内, 而将解压解密代码以及构造原PE文件代码块运行环境的代码, 程序入口点均放在. text块中。

4 实验结果与分析

本实验仅测试压缩变形模块, 暂不做非压缩变形处理。实验分为3组, 分别观察增强LZMA算法的保密性, 对文件的可执行性以及压缩效果。实验环境为Intel Core i5 - 2450M CPU2. 50 GHz, 4 GB RAM, Windows XP操作系统。

4. 1 保密性

随机选取本地的一个PE文件, 文件名为360sd_x64_4. 0. 0.4033A. exe, 即360杀毒软件的某一版安装包。暂不进行压缩编码, 仅观察加密前后的文件内容变化, 截取部分代码段如图2、图3所示。

当前系统时 间为2013 /9 /3 15∶39∶11, 用十进制20130903153911转化成二进 制数据并 进行低位 截断, 得到BTime。将文件名360sd_x64_4. 0. 0. 4033A转化成二进制串, 从而得到BName, 根据式 ( 1) 、式 ( 2) , 计算出密钥key[4], 该代码段长度lmod32 = 0, 可由式 ( 3) 得到密文, 恰如图3所示。

从实验结果可以看出, 经过加密之后的代码与原始代码差别较大, 在未经过解密的情况下, 文件不具备原始文件的功能, 甚至不能正常执行。该组实验表明, 增强LZMA算法提升了传统开源LZMA算法的保密性。

4. 2 对文件可执行性的影响

随机选取本地的10个PE文件, 观察压缩变形前与压缩变形后的PE文件的可执行性, 实验结果如表1所示。从实验结果来看, 压缩变形前后的PE文件均能正常执行, 表明本文提出的增强LZMA算法不影响文件的可执行性。

4. 3 压缩效果

分别采用RLE、LZW、LZSS、bzip2 ( - 9) 、gzip ( - 9) 、传统LZMA算法以及增强LZMA算法对上述10个PE文件的可压缩数据部分进行处理, 比较各压缩方法的压缩比与压缩时间。压缩比R计算公式为:

其中Osize为压缩前数据大小, Csize为压缩后数据大小。压缩比、压缩时间分别如表2、表3所示。

从上述两表中可以看出, 在压缩比方面, 增强LZMA算法与RLE、LZW、LZSS、bzip2 ( - 9) 、gzip ( - 9) 及传统LZMA算法相比, 压缩比最小, 说明增强的输出元组与搜索策略可改善传统LZMA算法的压缩比; 在压缩时间方面, 增强LZMA算法与LZSS及传统LZMA算法相比, 压缩时间更短, 说明采用分块并行与增强的搜索策略可降低LZMA算法的压缩时间, 但压缩时间并没有如初期预计的大幅度缩短, 说明增强LZMA算法在分析提取可压缩数据、加密及并行通信方面还是产生了一部分时延, 与gzip与bzip2仍存在一定差距。RLE与LZW的压缩时间虽然较高, 但是压缩比较低, RLE更存在压缩后PE文件体积变大的现象。

5 结 语

本文提出了一种能够应用于代码变形领域的增强LZMA算法, 在传统算法的基础上, 针对加密、分块并行、字典搜索策略、输出元组以及加壳等五个方面做出了进一步的改进。该算法避免了传统数据压缩算法对PE文件的可执行性的影响, 同时兼顾了代码变形所要求的保密性, 在压缩比与压缩时间上也有了明显的提高。

摘要：为了使LZMA压缩算法满足代码变形所要求的文件可执行性与保密性, 并进一步提升压缩效果, 提出了在加密、分块并行、字典搜索策略、输出元组及加壳等五个方面的增强方法, 设计了一个基于增强LZMA算法的代码变形系统模型。根据增强LZMA算法的保密性、压缩效果以及对文件可执行性影响的实验结果表明, 该算法有很好的应用效果。

关键词：软件保护,代码变形,数据压缩,LZMA算法

参考文献

[1]王朝坤, 付军宁, 王建民, 等.软件防篡改技术综述[J].计算机研究与发展, 2011, 48 (6) :923-933.

[2]赵玉洁, 汤战勇, 王妮, 等.代码混淆算法有效性评估[J].软件学报, 2012, 23 (3) :700-711.

[3]吴丹飞, 王春刚, 郝兴伟.恶意代码的变形技术研究[J].计算机应用与软件, 2012, 29 (3) :74-77.

[4]李露, 刘秋菊, 徐汀荣.PE文件中脱壳技术的研究[J].计算机应用与软件, 2010, 27 (9) :279-282.

[5]Li L, Liu Q J, Xu T R.Research and Implementation of Compression Shell Unpacking Technology for PE File[C].Information Technology and Applications, 2009.IFITA'09.International Forum on.IEEE, 2009, 1:438-442.

[6]Tu Z J, Zhang S Y.A Novel Implementation of JPEG 2000 lossless coding based on LZMA[C].Computer and Information Technology, 2006.CIT'06.The Sixth IEEE International Conference, 2006:140-144.

[7]Mohammed M H, Dutta A, Bose T, et al.Deliminate-a fast and efficient method for loss-less compression of genomic sequences[J].Bioinformatics, 2012, 28 (19) :2527-2529.

[8]张中华.PE程序加壳的研究与实现[D].北京:北方工业大学, 2009.

[9]刘坚, 李胜乐, 王子影.基于LZMA的数据库压缩存储应用研究[J].大地测量与地球动力学, 2009, 29 (6) :144-147.

[10]庞立会, 胡华平.PE文件加密保护技术研究[J].计算机研究与发展, 2006, 43 (S2) :249-254.

曲线驱动的网格变形算法与应用 篇2

关键词：鞋楦,定制CAD,拉普拉斯网格变形,特征线

0 引言

鞋楦是制作鞋子的模具,其设计的好坏直接决定鞋子的合体性与穿着的舒适性。传统的鞋楦设计需要由设计师对物理样楦进行反复修改,设计周期长,对操作人员个人能力的依赖性强。随着三维激光扫描仪的推广普及和逆向工程技术的进步,快速获得标准样楦的数字模型成为可能,这样我们就可以基于数字样楦进行鞋楦的再设计。为此必须提出相应的楦型设计方法并开发配套的鞋楦定制CAD系统[1,2]。

文献[1]基于距离映射构建楦型编辑系统,并提出了一种鞋楦局部修改方法,但鞋楦变形不够自然。文献[2]从脚型数据出发,在抽取脚型参数的基础上,通过经验公式获得对应的鞋楦参数,从而构建个性化的鞋楦模型,但该方法难于通过界面交互直接对楦型进行编辑修改。

本文基于笔者自主研制开发的三维楦型/脚型激光扫描系统,以曲线驱动三角网格曲面发生变形为主体框架,构建楦型快速定制CAD系统。系统通过编辑特征线、截面轮廓曲线和侧影轮廓线等实现鞋楦再设计,不需点选许多约束变形点,操作方便,实时直观;直接通过交互编辑B-Spline曲线来修改轮廓外形,不需要计算源曲线与目标曲线间的对应点;可以利用成熟的B-Spline曲线编辑技术,在曲线形状调整满意后,再实施变形物体的自由变形,减少计算量,提高设计效率和显示速度。相对于自由变形(free form deformation,FFD)技术,笔者采用的基于微分坐标的变形方法,不但可以保留细节,还可以交互指定任意变形区域,控制楦型设计的影响区域,使设计更加灵活。

1 离散拉普拉斯网格变形

针对三维网格模型的变形与形状修改,Olga等[3]提出了基于离散拉普拉斯坐标的形状编辑算法,该算法将三维网格顶点的拉普拉斯坐标定义为顶点与其一阶邻域顶点集的质心之间的偏差矢量,在编辑过程中努力保持多边形的拉普拉斯坐标,从而保持原多边形的视觉特征。该方法具有许多优点,如保留细节变形、线性求解,可以很容易地添加各种线性约束条件,变形时尽可能保持原多边形的视觉特征,使得变形满足用户的各种需要等。

1.1离散拉普拉斯算子

给定一个具有n个顶点的三角网格模型M=(V,E,F),V为顶点集,E为边集,F为三角面集合。对于网格中任意一点,其相对坐标(δ坐标)可以线性表示为其绝对坐标与其直接邻域内顶点质心坐标的差[4,5],即

式中,wij为vj点的权值;δi为点vi的拉普拉斯坐标(δ坐标);L为网格M的拉普拉斯算子。

为了便于表示和计算,网格模型的绝对坐标和相对坐标之间的转换采用矩阵形式表示。设网格模型M的邻接矩阵为A,I为单位矩阵,D为对角矩阵,其元素Dii =di,di为vi 1-ring邻域内的顶点个数。网格M的拉普拉斯算子矩阵L可以表示为

L=I-D-1A (2)

网格M上各个顶点的拉普拉斯坐标矩阵Δ为

Δ=LV

1.2基于δ坐标的网格变形

我们的目标在于给定一初始网格模型,将其中一些顶点移动到目标位置,在尽可能保持其各顶点相对坐标和初始网格模型视觉特征的基础上,求解其变形后各顶点的绝对坐标。因为网格的拉普拉斯算子矩阵L的秩小于n,是奇异矩阵,也就是说δ坐标具有平移变换的不变性,矩阵L是不可逆的,不能仅由L唯一确定网格模型各顶点的绝对坐标。要唯一确定网格模型各顶点的绝对坐标,必须求解一个满秩的线性方程系统,为此必须增加相应的约束条件。

设Q为初始网格模型,变形后的网格为Q′。Q中的部分顶点vj∈Cb作为约束条件,其变形后的绝对坐标已知,v′j=uj,j∈C,C={1,2,…,m}。这些点我们称之为“锚点”,这时的线性方程组扩展为

其中,V′为n×3阶矩阵,表示网格模型变形后各个顶点的绝对坐标;G为m×n阶矩阵,其中每行只有一个非零元素,对应于每个约束点,其值为

$g{}_{kj}=\{\begin{array}{l}\omega j=i{}_k\in C\\ 0\text{其}\text{他}\end{array}1\le k\le m,1\le j\le n$

式中,ω为约束点的权值。

F为m×3阶矩阵,每列代表一个坐标分量,这样以此即可求出3个坐标分量Fk:

Fk=ω uj

约束条件的添加使系统变为超静定,但扩展后的拉普拉斯矩阵$\tilde{\mathit{L}}$是满秩的,因而可获得最小二乘意义上的唯一解:

${\tilde{V}}^{\prime }=\underset{\begin{array}{l}\\ V^{\prime }\end{array}}{{\displaystyle \text{a}\text{r}\text{g}\text{m}\text{i}\text{n}}}(\parallel \mathit{L}{\mathit{V}}^{\prime }-\mathit{\Delta }\parallel {}^2+{\displaystyle \sum _{j\in C}\omega }{}^2\parallel {\mathit{v}}^{\prime }{}_j-\mathit{u}{}_j\parallel {}^2)$ (4)

${\mathit{V}}^{\prime }=(\tilde{\mathit{L}}{}^{\text{Τ}}\tilde{\mathit{L}}){}^{-1}\tilde{\mathit{L}}{}^{\text{Τ}}\mathit{B}$

矩阵$\tilde{\mathit{L}}$T$\tilde{\mathit{L}}$是稀疏且正定的,可以通过Cholesky分解快速求解该线性系统,并改善系统求解的稳定性。

2 曲线的获取

对楦型的修改主要是通过对鞋楦特征线(如鞋楦底板线、统口线等)和一些截面轮廓线(如跖趾围线、后跟弧线等)的编辑实现的,下面主要针对三角网格模型,说明其上的截面轮廓曲线和特征线的生成方法。

2.1截面轮廓线

首先,我们用点、面、半边和边的数据结构来表示三角网格模型,并建立其拓扑邻接关系。在进行三角网格模型拓扑重构后,即可利用模型中各个三角网格之间的拓扑信息计算平面与三角网格模型的交线。对于一个平面Γ,首先计算第一个与该平面相交的三角片T,得到交点坐标,然后根据局部邻接信息找到相邻的三角片,并求出交点,依次追踪,直到回到T,并得到一条有向封闭的轮廓环;重复上述过程,直到所有轮廓环计算完毕,并最终得到该层完整的截面轮廓[6]。

设已知平面Γ的方程为

A(x-x0)+B(y-y0)+C(z-z0)=0

其中,平面的法矢量为(A,B,C),且经过点(x0,y0,z0)。建立三维点数组Intersection来存储交点信息,其中,交点的ID记为交点所在边的ID,则求平面Γ与三角网格模型交点的具体步骤如下。

第一步,在完成模型拓扑重建工作后,遍历模型的实体边表,按序取出边表中的边,逐一判断,求出第一条与Γ相交的边,具体过程为:

(1)取第r条边,获取当前边的两端点Vs和Ve。

(2)分别将Vs和Ve的坐标值代入平面函数F(x0,y0,z0)=A(Vx-x0)+B(Vy-y0)+C(Vz-z0),可得对应的函数值分别为value_s和value_e,Vx、Vy和Vz为点的坐标分量。

(3)对value_s和value_e进行判断:若value_s=0且value_e≠0,说明平面与边的交点即是端点Vs,将Vs加入数组Intersection;若value_e=0且value_s≠0,说明平面与边的交点即是端点Ve,将Ve加入数组Intersection;若value_e=0且value_s=0,说明整条边都位于平面上,此时认为边的中点为交点,将其加入数组Intersection;若value_e与value_s异号,说明该边与平面有交点,求交,并将求得的交点加入数组Intersection;若value_e与value_s同号,说明该边与平面没有交点,则从边表读入下一条边,进行求交判断。

(4)当边表遍历结束还未找到与平面相交的边时,这说明平面与模型不相交。

第二步,根据截线的封闭性,可记录第一条与平面相交的边ID为FirstEdge,将其作为求交循环结束的标志;取该边的两条半边中的任意一条作为求交循环的初始半边He_init。具体过程如下:

(1)由初始半边开始,将其作为当前半边。

(2)取当前半边的上一条半边He_p与下一条半边He_n进行判断:若He_p所在的边与平面有交点,根据平面最多与当前三角面片的两条边相交的原则,无需再判断He_n所在边是否与平面相交,只需直接求出He_p所在的边与平面的交点,并将交点加入Intersection即可;若He_p所在的边与平面没有交点,则取He_n所在边与平面求交,将求得交点加入Intersection。

(3)判断交点所在边是否为初始边FirstEdge,若不是,则令求得交点的半边的伙伴半边He_a为当前半边,转第(2)步继续求交;若是,则求交结束。

以上求交方法的优点在于:①利用拓扑关系,使截面得到的交点集合是有序的,无须再进行排序,可直接获得首尾相连的有向封闭截线;②在三角面片与平面求交时,对某个三角面片只需计算一个边的交点,由面的邻接关系,可继承邻接面片的一个交点。

在得到截面轮廓上的有序交点集合后,采用三次非均匀B-Spline曲线对其拟合,便可利用成熟的B-Spline曲线编辑技术对截面曲线进行修改了。

2.2楦型特征线

这里的楦型特征线,主要是指楦底板线和鞋楦统口线,它们的形状和尺寸决定楦底和楦口的形状和尺寸,因而需要首先把这两条特征线提取出来,通过对它们的交互修改实现楦型的再设计。特征线的提取主要分为两步,即特征点的提取与特征线生成。

2.2.1 特征点提取

特征线的生成以提取出模型的特征顶点为前提,这里我们采用主曲率极值法来判断特征顶点[6,7]。首先计算点的法矢,对三角网格曲面上一顶点P,由其邻接三角形法矢的加权叠加计算该顶点的法矢;然后由P点及其邻接顶点在一局部坐标系内拟合出一个二次曲面S(u,v)=(u,v,h(u,v)),其中h(u,v)=au2+buv+cv2,如图1所示。曲面在P点的法曲率k中的极小值k1和极大值k2称为主曲率。k1、k2对应的法截线在P点的切线分别为m1、m2,这两条切线的方向称为主方向,两切线总是互相垂直的。

在计算出各顶点的主曲率之后,可通过对主曲率的极值进行判断得出特征点集,具体方法如图2所示。设m1及其反向延长线与三角形的交点分别为A、B,当Pi点主曲率k1的绝对值大于A、B两点在切线m1的法曲率k的绝对值时,则Pi点就是在m1方向上的曲率极值点,标示为特征点。A、B点的法曲率可由Vj、Vj+1在曲线m1的k1线性组合求得。该方法同样适用于k2,可判断Pi在曲线m2是否为特征点。经过上述判断得出的特征点还需经过进一步过滤,得出最后可应用于生成特征线的特征点集。

图3所示为对鞋楦三角网格模型进行特征点提取的示意,网格顶点数初始为18 064,最后提取的特征点数为3824。

(a)鞋楦三维模型 (b)提取的特征点

2.2.2 特征线生成

在得出鞋楦上的特征点集后,我们根据需要将其分离为统口线点集和楦底板线点集两部分。由于特征点集的分布类似于带噪声的测量点云,我们首先采用移动最小二乘法(moving least-squares,MLS)对特征点集进行细化和光顺;然后对细化后的特征点集进行排序和稀疏处理,使之有序化,便于采用B-Spline构造特征曲线。

设特征点集R经过MLS细化后所得点集为R′={P′i=(x′i,y′i)|i=1,2,…,N},从中随机选出一点P′j,根据需要,指定一个整数K,将其作为点P′j的K-邻近顶点个数,该值的大小决定了特征点集稀疏的程度。如果点集R′足够小,点P′j及其邻域内的点近似在一条直线上,搜寻P′j的K-邻域中距离P′j最远的点作为下一个搜寻点。点集排序与稀疏的算法流程如下:

(1)从点集中随机选出一点P′j作为初始点,计算其K-邻域Al,计算Al内各点到P′j的距离,存入数组dj,并按照从大到小的顺序对dj排序。把P′j存入排序后的点列数组NewP。

(2)设Al中距离P′j最远的点为P′j+1,向量F=P′j+1-P′j表示点P′j到P′j+1的方向(P′j+1、P′j分别为点P′j+1和P′j的坐标)。把P′j+1存入排序后的点列数组NewP。

(3)计算点P′j+1的K-邻域NP,同样计算NP内各点到P′j+1的距离,并按照从大到小的顺序排序。

(4)从已排序的NP中循环找出一点P′l(坐标为P′l),计算向量FN=P′l-P′j+1与向量F的内积e=FN·F,判断e值的大小,如果e>0,则将P′l存入排序后的点列数组NewP,以P′l作为新的P′j+1转入(3)进行迭代;如果e≤0,则继续在NP中循环。如果点P′l位于初始点P′j的K-邻域A中,说明点集构成一个封闭曲线,迭代终止。

对点云进行有序化稀疏后,即可利用比较成熟的有序点曲线重构技术对其进行曲线重构,本文采用3次非均匀B样条作为重构曲线类型。首先采用积累弦长参数化方法对有序化的点列进行参数化,从而确定曲线的节点矢量,由节点矢量反算三次B样条插值曲线的控制顶点,求解控制顶点时采用抛物线边界条件。在计算出曲线的控制顶点后,即可用德布尔算法对其进行正算,得出曲线上的点,其具体计算细节请参照文献[8]。图4所示为由楦底板特征点生成的特征曲线。

3 曲线驱动的网格变形

3.1曲线与网格模型间对应点的计算

要使曲线的变形能够驱动鞋楦模型做相应的变形修改,必须使模型与曲线相关联,使模型上的一些点与曲线上的采样点建立一一的映射关系,这样在曲线修改变形的时候就能带动网格模型做相应的变形。对于截面轮廓曲线,在构建曲线过程中,已自然建立了模型与曲线的关联,而对于特征线,由于对特征点进行了稀疏,曲线不完全经过特征点集,必须重新计算曲线上采样点与模型上顶点之间的映射关系,使特征线与楦型相关联。

首先以模型的平均边长为参考值,在特征线上按照等参变化取一系列的采样点c(ti),在模型上找其对应的最近点Vi,如果c(ti)与Vi的距离小于某一指定阈值,认为Vi为模型上的曲线关联点之一,记录Vi点的ID和其在曲线上的对应参数值ti。

这里有两个问题需要注意,一是对于曲线上的两个甚至三个相邻采样点,它们在模型上的最近点有可能是同一点,这时只记录距离最小的那一组对应点;二是对于曲线上的采样点c(ti)而言,其在模型上的最近点为Vi,但反过来,Vi在曲线上的投影点并不一定是c(ti),须进一步计算其在曲线上的f oot point点,修正其在曲线上的对应位置坐标。

如图5所示,设曲线上采样点c(ti)在模型上的最近点为Vi,通过牛顿迭代计算点Vi在曲线上的对应参数值t*。要计算t*,可通过最小化点Vi到曲线的距离获得,即使下式最小

g(t)=‖c(t)-Vi‖2 (5)

其中,c(t)为曲线上任一点的坐标;Vi为点Vi的坐标,采用牛顿迭代计算式(5)使g(t)最小以获得t*。设t(i)为当前迭代值,下一次迭代值为[9]

迭代的初始值取Vi在曲线上的初始对应点c(ti)的参数值ti。

建立了曲线与模型间的关联之后,这些与曲线关联的模型顶点将作为位移约束加入到求解网格变形的线性系统中。

3.2变形区域指定

在对楦型进行编辑时,根据设计意图,必须在鞋楦模型上指定一个变形区域,也称为兴趣区域(region of interest,ROI)。一方面这是鞋楦编辑的需要,即需要设计者按照设计意图制定所编辑区域的影响范围,以保证鞋楦的其他部分不被修改;另一方面,由于网格模型编辑的计算量与顶点数有关,指定区域后,可以减小计算量。

变形区域的指定按照鞋楦编辑的功能模块指定,重点在于如何构建其邻接矩阵和处理其边界。鞋楦编辑中,楦底编辑、后跟弧线编辑、顶弧编辑以及底弧编辑等的变形区域相对比较规范,其指定方法在后文叙述。要对楦型进行局部编辑,需要交互指定一个局部变形区域,指定过程如下:

(1)建立每个三角面包含的三个半边结构,记录每个半边所对应的三角形的ID[10],如图6所示。

(2)设网格曲面的法矢方向指向模型外部,通过人机交互,用鼠标按逆时针方向点选模型上的一系列顶点,其中最后一点需和第一个点重合。

(3)由相邻的两点生成dijkstra近似测地线,得到一个由网格三角形的边组成的封闭环,构成所选区域的边界,标记边界点为已访问。

(4)从边界点中选择两相邻顶点Vi和Vi+1,则ViVi+1构成一条半边,选择该半边所在三角形的第三个顶点P作为选择区域搜寻的种子点,如图7所示。

(5)从种子点P开始,依托模型的拓扑邻接信息,按照广度优先算法搜索边界内的所有顶点及其对应的三角形,直到边界内所有点都标记为已访问。

在获得选择区域内的顶点后,就可以根据其拓扑邻接信息建立该选择区域的邻接矩阵,并进一步建立其拉普拉斯算子矩阵L。其中需要特别处理的是所选区域的边界点,因为一方面在构建邻接矩阵时,边界上的点需要特殊处理,只计那些处于选择区域内的邻接点;另一方面,需要对这些边界点作出标记,以便在做楦型修改时将其作为边界约束加入到求解网格变形的线性系统中。

3.3曲线约束的网格编辑

我们的目标在于通过交互修改曲线来驱动鞋楦模型做相应的变形,达到按设计意图进行楦型定制的目的。在获得指定区域的拉普拉斯算子矩阵L后,我们将把边界约束和曲线修改引起的位移约束加入求解系统,构建式(3)。

所选区域边界上的顶点将作为“锚点”约束加入求解系统,即这些边界顶点在变形过程中保持位置不变,以保证变形区与非变形区域的平滑过渡。可以通过调整边界点的权值来调整过渡区的变化。

对于模型中与曲线直接关联的顶点,将作为位移约束加入求解系统。设模型中顶点v与曲线上的点c(ti)相关联,曲线修改后,点c(ti)变为q(坐标为q),则认为v点修改后的位置为

v′=v+(q-c(ti))

在边界约束和位移约束条件下求解式(4),即可求的修改后的鞋楦模型的顶点位置。

4 鞋楦CAD

基于前述方法,我们采用VC++6.0和OpenGL技术在微机平台上实现了基于样楦的鞋楦定制CAD系统。系统主要分为文件管理模块、楦型重构模块、楦型编辑设计模块、数据库管理模块和视图管理模块等。其中楦型编辑设计模块是本系统的主体模块,分为楦底板编辑、后跟编辑、底弧编辑、顶弧编辑、楦头编辑和局部区域编辑以及楦型放样等子模块。其他功能模块不是本文的主要内容,这里不做详细介绍。

各个模块的功能和操作面板基本相似,以楦型底板编辑模块为例,首先针对楦底编辑的要求,选定一变形区域作为设计变更的影响区域,如图8所示。

而对于特征线的生成,可以采用灵活的方法,可以按照前述方法由系统自动生成,也可以通过手工交互输入,甚至可以通过外部导入。曲线的修改方法如下:主要通过交互拖动曲线上的型值点不断调整曲线的形状,在曲线形状修改满意后,再预览整个鞋楦的变形情况。这样把曲线修改和网格变形分步处理,可以减少计算量,提高设计效率和实时显示速度。如果对预览模型的编辑结果不满意,还可以继续修改曲线,直到获得满意的设计效果,最后确认设计结果。图9所示为通过调整楦底板线进行楦型再设计的结果。

相较于自由变形技术,本文方法可以指定变形区域,不会发生变形干涉现象,在变形的同时可保持细节特征。

5 结论

(1)采样曲线驱动三维网格变形方法实现楦型的编辑修改,不需点选许多约束点,易操作,交互性好,符合设计师设计习惯,容易表示设计意图。

(2)可以任意指定设计影响区域,设计灵活;采样先修改曲线,满意后再驱动三维模型作相应的更新,计算效率高,实时显示。

(3)采用离散拉普拉斯网格变形方法,容易增加各种线性约束条件,使得变形满足用户的各种需要。

参考文献

[1]Leng J,Du R.A CAD Approach for Designing Cus-tomized Shoe Last[J].Computer-aided Design&Applications,2006,3(1/4):377-384.

[2]徐从富,刘勇,蒋云良,等.个性化鞋楦CAD系统的设计与实现[J].计算机辅助设计与图形学学报,2004,16(10):1437-1441.

[3]Olga S,Yaron L,Daniel C O,et al.Laplacian Surface Editing[C]//Proc.of Eurographics Symposium on Geometry Processing.Nice,France,2004:179-188.

[4]Marc A.Differential Coordinates for Local Mesh Morphing and Deformation[J].The Visual Comput-er,2003,19(2):105-114.

[5]Nealen A,Sorkine O,Alexa M,et al.A Sketch-based Interface for Detail-preserving Mesh Editing[J].ACM Transactions on Graphics,2005,24(3):1142-1147.

[6]上官宁.基于三维网格模型的鞋楦CAD关键技术研究[D].泉州:华侨大学,2008.

[7]刘胜兰,周儒荣,张丽艳.三角网格模型的特征线提取[J].计算机辅助设计与图形学学报,2003,15(4):445-448.

[8]施法中.计算机辅助几何设与非均匀有理B样条[M].北京:高等教育出版社,2001.

[9]Hu S M,Wallner J.A Second Order Algorithm for Orthogonal Projection onto Curves and Surfaces[J].Computer Aided Geometric Design,2005,22(3):251-260.

变形算法 篇3

1 块匹配搜索算法

数字散斑相关搜索法中多采用基于区域匹配的灰度匹配法,所以视频图像处理中的基于块匹配的运动估计快速搜索法的有关思想可以运用到该搜索算法中[4]。视频图像处理中是跟踪最小块误差(Minimum Block Deviation,MBD)点,在DICM中改为跟踪相关系数最大(Maximum Correlation Coefficient,MCC)点。本文选用标准化协方差相关函数,因为其相关系数矩阵呈单峰分布,且相关峰更尖锐,定位精度高。具体的形式如式(1)所示。

式(1)中,f(xi,yi)为参考图像中图像子区中点(xi,yj)处的灰度值;为变形后图像中图像子区中点()处的灰度值;为参考图像中图像子区的平均灰度;为变形后图像中图像子区的平均灰度。

典型的块匹配快速搜索算法有全搜索法(Full Search Method,FS)、三步搜索算法(Three Steps Search Method,TSS)、四步搜索法(Four Steps Search Method,FSS)和菱形搜索法(Diamond Search Method,DS)[5],这四种算法都是通过改变搜索方式,从而达到减小搜索区域提高搜索速度的目的。

2 仿真计算

2.1 仿真散斑的生成

为了对比不同算法的搜索效率,本文根据Peng Zhou提出的算法[6]以及上述四种块匹配搜索算法的基本原理利用MATLAB软件自行编制了一套针对木材变形测量的分析程序,该程序能够生成精确控制位移的数学仿真散斑图,并且可以准确地控制图像大小、散斑数量及散斑亮度;同时可以选取不同的块匹配搜索算法对生成的散斑图进行数字散斑相关分析,得到每像素点的位移。本文采用的仿真散斑大小为256×256 pixels,散斑数量为800,散斑大小为4 pixel。图1 (a)、(b)分别为MATLAB软件生成的散斑图及该散斑图在x正方向平移1 pixel后的散斑图。

2.2 仿真结果及分析

为了对比分析单一方向存在位移以及x方向和y方向均存在位移时的不同情况,首先生成如图1所示的四组散斑图,每组散斑图的x方向和y方向位移如表1,第一组和第二组表示单一方向存在位移的情况,第三组和第四组表示x方向和y方向均存在位移的情况。每组分别有10对平移前后的模拟散斑图,以考察算法的稳定性。

分别使用全搜索法、三步搜索法、四步搜索法和菱形搜索法计算每像素点的整像素位移,可以得到不同位移大小下四种搜索算法的搜索时间。表2给出了四种搜索算法搜索每组仿真散斑图对所用的平均搜索时间,由表2可以看出由于全搜索法的搜索区域最大,其搜索时间远远大于其他三种算法,如果处理较大的图像,全搜索法计算效率低的缺点将会更为明显。

3 试验验证

为了验证四种块匹配算法在木材弯曲变形测量中的有效性,本文首先进行了刚体平移试验,验证不同算法的精度,然后进行了木材三点弯曲验证试验。

3.1 试验材料及设备

试验材料为白桦,含水率约为10%,尺寸为200mm×20 mm×20 mm,设备为电子万能试验机(型号为RGM—2050),精度可达0.001 mm。本文利用黑色与白色的自喷漆,在试件表面交替喷涂形成随机分布的人工散斑场。试验过程中,试件表面由LED光源照明,被测试件的图像由放置在其正前方约0.1 m的图像采集系统采集。图像采集系统包括一个分辨率为1 280 pixel×1 024 pixel的CMOS数字摄像机(型号MV—130UM)和变焦镜头(AFT—0850ZMP)。调焦清晰成像后,图像中每像素所对应的空间尺寸约为25μm。

3.2 试验结果与分析

利用自行编制的MATLAB程序进行数字散斑相关处理,选取有效区域为研究对象,有效区域位于图像中间位置的128 pixel×128 pixel区域(如图2),得到搜索时间与搜索精度如表3(表中计算位移为通过块匹配算法计算得到的位移,试验位移为由试验机测得的实际位移),其中刚体平移试验中位移为128×128个像素点的平均位移,三点弯曲试验中位移为试件弯曲的挠度值。

对于刚体平移,其前后图像如图3所示。由表3可以看出,全搜索算法的搜索时间为3.179 9 s,远远大于其他三种算法;三步搜索算法虽然搜索时间最短,但是68.6%的相对误差说明其匹配错误率较高;四步搜索法和菱形搜索法位移相对误差仅为1.49%,表明在误差允许范围内,四步搜索算法和菱形搜索算法可以用于木材弯曲变形的测量。

木材三点弯曲前后的图像如图4所示,由表3可以看出,全搜索法搜索时间为2.991 3 s,远远大于其他三种算法,这点与仿真试验一致,表明全搜索法搜索效率低;三步搜索算法的搜索时间为0.141 9 s,在四种算法中计算速度最快,但是三步搜索算法的位移相对误差(69.7%)远大于全搜索法(1.99%)、四步搜索算法(1.99%)和菱形搜索算法(1.99%);四步搜索算法的搜索时间(0.286 3 s)略大于菱形搜索算法的搜索时间(0.261 9 s),但当处理较大的图像时,这种搜索时间上的差距会更加明显。因此综合考虑搜索时间与搜索精度,对于木材弯曲变形测量,菱形搜索算法最优,四步搜索算法次之。

4 结论

本文考察了四种常用的块匹配搜索算法在利用数字散斑相关方法进行木材弯曲变形测量的有效性。综合考虑仿真计算与木材三点弯曲验证试验,结果表明,对于木材弯曲变形测量,全搜索算法虽然搜索精度很高,但搜索时间远远大于三步搜索算法、四步搜索算法和菱形搜索算法,其搜索效率低;三步搜索算法虽然计算效率最高,但搜索精度较低;四步搜索算法和菱形搜索算法搜索精度较高,菱形搜索算法的搜索时间略少于四步搜索算法。因此在进行木材弯曲变形测量时,菱形搜索算法最优,四步搜索算法次之。

参考文献

[1]江泽慧,费本华,张东升,等.数字散斑相关方法在木材科学中的应用及展望.中国工程科学,2003;5(11):1-5

[2]孙艳玲,赵东,高继河.数字散斑相关方法在木材断裂力学上的应用分析.北京林业大学学报,2009;31(1):206-210

[3] Chen Yongsheng,Hung Yiping,Fuh Chioushann.Fast block matching algorithm based on the winner-update strategy.IEEE Transactions on Image Processing,2001;10(8 ):1212-1222

[4]伍卫平.图像相关技术的亚像素位移算法与实验研究.武汉:华中科技大学,2009

[5]杨高波,杜青松.MATLAB图像/视频处理应用及实例.北京:电子工业出版社,2010

[6] Zhou P,Goodson K E.Subpixel displacement and deformation gradient measurement using digital image/speckle correlation.Opt Eng, 2001;40(8):1613-1620

[7]王怀文,亢一澜,谢和平.数字散斑相关方法与应用研究进展.力学进展,2005;35(2):195-203

变形算法 篇4

1 系统组成

1.1 硬件选择

图1为本系统结构框图, 其组成由采集模块、数据传输和数据处理单元组成。采集模块主要是安装在井壁的压力和温度传感器, 数据传输单元包括有单总线、RS485总线和RS232总线, 数据处理单元包括有上位机、监测仪表R T U;监测仪表R T U包括有主控器MCU以及与主控器MCU连接的网络接口控制、液晶电路、报警电路, 主控器MCU内置有RTU管理控制程序。图中1为多路温度传感器;2为多路压力传感器;3为上位机;4为485转232接口;5为电源;6为监测仪表RTU。

1.2 工作原理

温度传感器和压力传感器采集的数据通过单总线发送到监测仪表RTU的主控器MCU, 并将所得数据显示在液晶显示器上, 同时转换成RS485数据送到RS485总线上;上位机上设有RS232接口电路和相应软件、485总线上传输的数据通过232/485转换器转换成PC机能识别的电平, 将数据经RS232总线送到上位机的RS232接口, 随即送入上位机中的萤火虫算法模型软件中, 经过演算可以获得较准确的未来发展趋势预警, 并将数据实时显示在上位机上。便于工作人员查看。

压力传感器包括检测部分和输出芯片DS2438, 由于输出为数字量, 故可以将所得压力信号直接挂到单总线上。对传感器所得信号在井下采用单总线传输, 地面采用485总线。

2 萤火虫算法的理论基础

2.1 模型的说明

井壁在坍塌前都会有一定的变形, 尽管井壁变形有多种表现类型, 但井壁坍塌的直接原因都是由于井壁所受压力过大, 此外温度对井壁的承受能力也有一定的影响, 所以本系统需要检测井壁各处的压力和温度即可实现对井壁变形程度和状态的掌握, 对这些数据进行在上位机上进行萤火虫算法演算, 得出井壁是不是处于危险状态。

2.2 井壁预警的萤火虫算法实现

采用萤火虫算法对矿井井壁的安全进行预测, 通过上位机软件计算后, 密集度最高的区域对应着比较危险的部位, 需要加大维护力度, 防患于未然。萤火虫算法的基本步骤是, 首先将萤火虫群体随机散布在解空间, 处于不同位置的不同萤火虫发出的荧光亮度也不同, 其次, 根据荧光亮度的不同, 荧光强度低的萤火虫会飞向比自己亮的萤火虫, 飞行的路程由吸引度的大小来确定。在实际应用中为了能够在更大范围内搜索, 找到全局的最优解防止局部最优解干扰, 导致算法不收敛, 在计算萤火虫更新后的位置时需要在位置更新过程中增加扰动项。最后, 通过多次迭代后, 算法收敛于亮度最高的萤火虫, 也就是说该萤火虫会把其他萤火虫都吸引到自己的位置处, 便完成了对最优解的寻找, 达到精确预测的目地。

3 系统仿真测试

在计算机中, 对建立的模型, 对采集的部分数据进行萤火虫算法防真, 运行得到如图2的结果, 所用萤火虫收敛于几个点, 也即在些地方需要加强防范, 很有可能会在这些地方出现事故。

4 结语

本文完成了一种基于萤火虫算法的井壁变形监测系统, 该系统具有如下特点:完全实现自动化。从数据采集到数据传输到数据处理, 整个过程完全由计算机控制, 不需人的参与。 (2) 实时监测。位于井下的传感器可以实现连续采集数据, 所得数据被传输至地面中心计算机, 整个系统可以实现工作。 (3) 报警功能。工作人员在地面上位机可以设定报警阈值, 软件的输出达到或者超过该设定值是系统便会发出报警信号。 (4) 速度快。整个系统从数据采集到最终结果显示在终端设备和PC上, 仅需几分钟。

摘要：本文针对煤矿井壁变形问题设计了一种实时监控系统, 来预防井壁事故的发生。结合现代检测、通讯技术和萤火虫算法, 采用温度传感器, 多路压力传感器测量井壁的变形信息, 并将各检测单元所得数据经由MCU演算, 同时把数据传送到地面中心计算机, 地面工作人员即可实时地掌握井壁的形变情况。该系统可实时测量和监控, 其检测精度高、鲁棒性好, 且可根据需要接入Ethernet或者LAN, 实现远程监控和数据共享。

关键词：井壁变形,萤火虫算法,煤矿预警,实时监测

参考文献

[1]刘延强, 吕英民, 蔡强康.井壁变形和摩擦阻力对底部钻具组合的影响[J].石油大学学报 (自然科学版) , 1991, 01:32-41.

[2]梁化强.矿井井壁变形监测及治理研究[J].矿业安全与环保, 2008, 05:45-46+55.

[3]张黎明, 王在泉, 任禀洁.井壁安全远程自动监测及井壁变形的灰色马尔柯夫预测[J].煤矿安全, 2005, 12:50-52.

[4]刘鹏飞.井壁变形全自动监测系统的工程应用[J].煤炭科学技术, 2005, 02:7-9+6.

变形算法 篇5

关键词：变形监测,小波去噪,GPS,多路径

0 引 言

目前, 用于测量结构物的振动频率和振幅的仪器主要是加速度仪。加速度仪的局限在于需要进行双重积分, 且其无法获得绝对的位移;并且在大多数情况下不能实现实时或准实时的处理, 因而无法完全获得结构物的动态特征。而GPS则弥补了这一缺点。目前的GPS采样率已经最高达到了100 Hz, 而且基于单历元的解算算法也日趋完善, 并且GPS具有可以全天候全天时工作, 实时获得三维的动态坐标的优点。因此GPS用于结构物的变形监测也越来越受到国内外学者的重视。Kijewski[1], Nakamura[2], Ashkenzai[3]等学者研究如何将GPS技术应用到来 测量台风对高层建筑和桥梁的影响。德国的Geo++公司[4]还开发出基于GPS的大桥振动变形监测软件。西南交通大学的熊永良教授[5]和香港理工大学的丁晓利教授在这一方面也作了很多的研究。

GPS用于结构物的变形监测时, 基线一般都比较短, 与卫星, 接收机, 大气相关的误差均可以在差分定位中消除, 影响其定位的主要因素是多路径效应。在正常情况下, 卫星信号沿光程最小的路径从卫星直接到达接收机天线。如果天线附近有反射物体, 到达接收机天线的信号为直接信号和反射信号的叠加信号。反射信号对观测值的影响产生一个附加的时延量, 这种现象称为多路径 (如图1所示) 。在载波相对定位过程中, 多路径效应的影响可以达到5 mm, 严重时还会对模糊度的固定产生影响。

在静态测量过程中, 一般是通过选择远离反射体的地点架设天线来减弱多路径效应的影响。然而, 进行结构物的变形监测过程中, 往往无法避开反射源, 因此多路径效应在结构物的变形监测过程中更加明显, 如在桥梁的变形监测过程中, 附近的高楼, 桥梁的拉索都很容易导致多路径效应。由于卫星位置的周期重复性 (运行周期约为11小时58分02秒) , 所以其多路径效应也存在周期的相关性, 我们把这一部分多路径称为系统性的多路径, 也称为系统性的背景噪声。Wubbena提出利用多路径的周期性的相关性, 采用多日观测来削弱多路径的方法。该方法能够消除大部分的系统性的多路径误差。本文的提出来的方法是首先利用多路径的周期性高相关性, 消除大部分的系统误差, 然后再利用小波去噪的方法消除剩余的误差的影响, 从而提高结构振动监测的精度。实验证明该方法是非常可行的。

1 小波去噪

对于待分析的信号f, 其连续小波变换为:

$W{}_f(a,b)=|a|{}^{-\frac{1}{2}}{\displaystyle {\int }_{R}f}(t)\overline{\psi (\frac{t-b}{a})}dt(1)$

经常遇到的信号为离散信号 (GPS测量序列也是离散信号) , 此时采用离散小波变换的方法, 常用的算法是 Mallat金字塔算法, 设f (1) = (x${}_1^{(0)}$, x${}_2^{(0)}$, …, x${}_{{\rm N}}^{(1)}$) 为离散信号, 利用如下变换可将其分解为近似信号f (1) = (x${}_1^{(1)}$, x${}_2^{(1)}$, …, x${}_{{\rm N}}^{(1)}$) 和细节信号d (1) = (d${}_1^{(1)}$, d${}_2^{(1)}$, …, d${}_{{\rm N}}^{(1)}$)

$\begin{array}{l}f{}_k^{(1)}={\displaystyle \sum _{n}h}{}_{0(n-2k)}f{}_n^{(0)}(2)\\ d{}_k^{(1)}={\displaystyle \sum _{n}h}{}_{1(n-2k)}x{}_n^{(0)}(3)\end{array}$

其中

$\begin{array}{l}h{}_{0k}=<\text{ϕ}{}_{10}(t),\text{ϕ}{}_{0k}(t)>=\frac{1}{\sqrt{2}}{\displaystyle \int \text{ϕ}}(\frac{t}{2})\text{ϕ}{}^{*}(t-k)dt(4)\\ h{}_{1k}=<\psi {}_{10}(t),\psi {}_{0k}(t)>\frac{1}{\sqrt{2}}{\displaystyle \int \psi }(\frac{t}{2})\text{ϕ}{}^{*}(t-k)dt(5)\\ x{}_k^{(2)}={\displaystyle \sum _n<}\text{ϕ}{}_{1n}(t),\text{ϕ}{}_{2k}(t)>x{}_n^{(1)}(6)\\ d{}_k^{(2)}={\displaystyle \sum }^n<\text{ϕ}{}_{1n}(t),\psi {}_{2k}(t)>x{}_n^{(1)}(7)\end{array}$

其中ψ (·) 和ϕ (·) 分别为小波函数和尺度函数。

从上式可知基于Mallat算法的小波分解变换, 实际上是将信号通过低通和高通两组滤波器, 近似分量主要反映信号的低频部分, 细节分量主要反映信号的高频部分。小波重构算法是小波分解算法的逆过程, 依据先验信息, 将与噪声相应的高频细节信号置零, 然后进行重构, 得到的信号是去噪后的信号。

在本次实验中选用的是db4小波, 去噪采用的是软阈值法该方法的原理是:对一有限长度信号, 某一尺度时小波变换系数yj=Ajf (x) 可以表示为如下的形式:

$y{}_j=x{}_j+\delta z{}_j(8)$

其中xj为原是信号的小波变换系数, δ>0为常数, 是噪声级, zj为一个标准的白噪声, 服从正态分布N (0, σ2) 。若要求是混有噪声的信号yj=Ajf (x) 中获得xj这一原始信号, 首先计算含有噪声心好的正交小波变换, 选择合适的小波分解层数j, 然后将信号进行小波分解至j层, 得到小波分解系数。然后进行阈值处理。

$\begin{array}{l}\widehat{x}{}_j=sgn(y{}_j)(|y{}_j|-\delta )\\ =\{\begin{array}{cc}sng(y{}_j)(|y{}_j|-\delta )& |y{}_j|＞\delta \\ 0& |y{}_j|＜\delta \end{array}(9)\end{array}$

设共有N个离散点, 通常取$\delta =\sigma \sqrt{2\log {\rm N}}$, 上式实际含义为当|yj|>δ时, yj被认为是有用的信号, 予以保留, 否则认为yj主要是由噪声引起的, 应当消除。然后通过小波逆变换, 得到恢复的原始信号估计$\widehat{x}$j, 即去除噪声的信号。

2 结合小波去噪的振动监测方法

该方法的基本思想是将GPS天线固定在振动台上, 首先在振动台不振动时用GPS进行一天连续不间断的测量, 获得一天的连续观测值, 然后进行单历元解算, 获得静态时的坐标序列。由于GPS处于静止状态, 此时获得的坐标序列与该点实际坐标的差值即为系统性的多路径。然后在第二天让振动台振动, 此时用GPS测量, 并且用单历元结算软件进行求解, 获得的坐标序列与前一天对应于同一卫星分布的时刻的系统性多路径求差, 即后一天的坐标序列中t时刻坐标值与前一天t-3 min 56 s的坐标值求差, 得到消除系统性误差后的坐标序列, 然后再进行坐标域的小波去噪处理, 最终获得干净的坐标序列即为振动台的实际振动序列。具体的过程如图2所示。

3 振动试验及其分析

3.1 试验设计

模拟结构物振动的振动台能够模拟频率从0.025 Hz到1.8 Hz, 振幅从2 mm到50 mm 的振动。它既可以模拟x方向或y方向的单独的振动, 也可以做水平方向内的圆周运动。GPS天线可以固定在振动台上。同时该振动台内与一台电脑相连, 电脑可以控制振动台模拟各种振动, 并且以50 Hz的采样频率记录振动台的x, y方向的振动位移。

试验采用徕卡两台GX1230型双频GPS, 采样率为20Hz, 使用的天线为AT504型扼流圈天线, 用以减弱多路径效应的影响。基线长度为86 m。这样与大气、卫星和接收机相关的误差就可以在差分钟基本消除。在试验过程中, 第一天使振动台保持静止, 第二天的观测中GPS天线置于振动台上, 由计算机控制模拟各种频率和振幅的振动。试验的频率为0.025, 0.1, 0.5, 1 Hz。各种频率的均进行了振幅为2, 5, 10, 20, 40 mm的正弦运动。试验数据处理采用的是GPSSHM单历元解算软件。

3.2 多路径分析

当GPS天线, 周围的反射源均保持静止不动时, 它的多路径效应存在周期性的相关性。Radovanovic et al[6]通过实验发现相关度可以达到80%。为了确定在振动台振动过程中的GPS多路径效应是否和静止时的依然高度相关, 进行一组试验, 将静态时的多路径效应与振动时的进行比较, 结果如图3所示。上图为静止时的多路径效应图, 下图为振动时的多路径效应图, 它们对应于卫星运动周期的同一位置。动态时的多路径是将GPS测量的结果与电脑记录的振动台的振动相减求得的。从图4可以看出, 振动时测量出来的多路径效应与静态时的多路径是非常相似的, 其相关度为84.1%。

在消除系统性多路径效应后, 我们把测量值与实际振动值相减, 得到残余的噪声, 如图5所示。经假设检验证实, 其分布服从N (0, σ2) 的正态分布。

3.3 精度评定标准

用如下的公式定义GPS的测量误差:

$\begin{array}{l}Error\min (\%)=\frac{x{}_{g\text{m}\text{i}\text{n}}-x{}_{t\text{m}\text{i}\text{n}}}{x{}_{t\text{m}\text{i}\text{n}}}100\%(10)\\ Error\max (\%)=\frac{x{}_{g\text{m}\text{a}\text{x}}-x{}_{t\text{m}\text{a}\text{x}}}{x{}_{t\text{m}\text{a}\text{x}}}100\%(11)\\ Errorstd(\%)=\frac{x{}_{gstd}-x{}_{tstd}}{x{}_{tstd}}100\%(12)\end{array}$

式中:xgmin表示GPS测量出来的波谷的平均值;xtmin表示振动台振动实际振动的波谷的平均值, 它是由与振动台相连的电脑记录的;Error min表示测量波谷的相对误差;xgmax表示GPS测量出来的波峰的平均值;xtmax表示振动实际的波峰的平均值;error max表示测量波峰的相对误差;xgstd表示GPS测量得到坐标序列的标准差 (std) , xtstd表示振动实际的坐标序列的标准差, Errorstd表示测量的标准差的相对误差。

3.4 实验结果及其分析

3.4.1 频率为0.025 Hz的振动实验

实验中振动台的振动频率为0.025 Hz, 振幅从分别为2, 5, 10, 20, 40 mm。图5中虚线部分为消除了系统背景噪音的数据 (没有进行小波去噪) , 实线部分为振动台实际振动曲线。图7中虚线部分为GPS测量数据经过滤波后的最终的结果, 实线部分为振动台实际振动曲线, 从图中可以看出经过滤波后的数据更加接近真实振动, 表1为精度分析。

3.4.2 频率为0.1 Hz的振动实验

实验中振动台的振动频率为0.025 Hz, 振幅从分别为2, 5, 10, 20, 40 mm。表2为精度分析。

3.4.3 频率为0.5 Hz的振动实验

实验中振动台的振动频率为0.5 Hz, 振幅从分别为2, 5, 10, 20, 40 mm。表3为精度分析。

3.4.4 频率为1 Hz的振动实验

实验中振动台的振动频率为1 Hz, 振幅从分别为2, 5, 10, 20。表4为精度分析。

3.4.5 同一振幅不同频率的数据的比较分析

表5为振幅为20 mm时, 频率分别为0.025 Hz, 0.1 Hz, 0.5 Hz, 1 Hz的比较分析从表中可以看出随着振幅的增大, 频率的减少, 其测量误差也逐渐的减少。

4 结 论

从以上的实验和分析可以看出, GPS的背景噪音存在周期性的相关性, 在本次的实验中, 其相关性高达84.1%。基于系统性多路径的周期性重复这一特性, 我们设计实验中利用多周期的观测消除系统性多路径。通过对残余的背景噪音的分析发现, 背景噪音服从正态分布, 我们采用小波去噪的方法可以很好的消除这些误差。从实验分析中可知, GPS测量的精度与振动频率, 振幅紧密相关。精度随着振幅的增大, 频率的减少而增加。从上面的分析中也可以看出, 当振幅大于2mm, 频率小于1Hz时, GPS测量的精度是令人满意的。

参考文献

[1]Kijewski, T., Full Scale Measurements and System Identifica-tion:A Time-Frequency Perspective[M].PhD Thesis, Universi-ty of Notre Dame, United States, 2003.

[2]Nakamura, S., 2000.GPS measurement of wind-induced sus-pension bridge girder displacements[J], J.Struct.Eng.126, 1413-1419.

[3]Ashkenazi, V., and Roberts, G.W.Experimental monitoring ofthe Humber Bridge Using GPS[J].Proc.of the Institution ofCivil Engineers, 120, 177-182.

[4]Gerhard Wubbena.Permanent object monitoring with GPS with1 millimeter accuracy[J].ION GPS-2001.

[5]熊永良.基于小波变换的单历元解算算法及其在结构振动检测中的应用研究[J].测绘学报.2005, 34 (3) :202-207.

[6]Radovanovic R, W F Teskey, R Fox, et al.Short Baseline De-formation Detection[J].Technical Report, Department of Geo-matics Engineering, The University of Calgary.

[7]Ogaja, C.On-line GPS integrity monitoring and deformation a-nalysis for structural monitoring applications.Tech.Meeting ofthe Satellite Division of the U.S.Inst.of Navigation, Salt LakeCity.

变形算法 篇6

制造企业的实际生产中,许多产品需要大批量生产,如汽车零部件、轴承等,其主要生产组织方式为流水生产。在流水生产中如何确定产品的生产顺序和节拍至关重要,此类问题可归结为流水车间排序问题。不合理的生产顺序会导致设备得不到有效利用,降低企业的生产效率。

流水车间排序问题是典型的NP-hard问题,至今没有可以精确求得最优解的多项式时间算法。随着计算机技术的发展,各种智能技术在作业排序问题的研究中得到了应用[1~4]。遗传算法具有不依赖问题的具体领域、搜索结果全局近优等特点,一直是研究的热点;但它是一种搜索算法,在寻优过程中,防止早熟、局部收敛非常重要。遗传操作中的交叉操作在遗传算法中起着关键作用,是产生新个体的主要方法[5]。为了提高其求解流水车间排序问题的寻优性能,许多学者在交叉算子方面做了研究,如在交叉操作中混合使用LOX和单点交叉两种交叉算子来改善遗传算法求解Flow-Shop问题的质量[6];利用四种新的交叉算子:SJOX、SBOX、SJ2OX和SB2OX,解决交叉操作中子代适应度值小于父代的问题[7]。本文针对Flow-Shop问题,通过分析典型单点交叉算子的优缺点,根据其进化缓慢的缺陷,提出一种改进型单点交叉算子,以提高遗传算法的寻优性能。

1 Flow-Shop问题的数学描述

流水车间排序问题也常称为Flow-Shop问题,一般可以描述为:n个工件在m台机器上加工,每个工件需要经过m道工序,每道工序要求不同的机器。n个工件在m台机器上的加工顺序相同,工件在机器m上的加工时间是给定的。问题的目标是确定n个工件在每台机器上的最优加工顺序,使最大流程时间达到最小。

令Pij为工件i在机器j上的加工时间,Cij为工件i在机器j上的完工时间,不失一般性,假设各工件按机器1至m的顺序进行加工,令为所有工件的一个排序,以上定义和假设用数学模型表示为[8]:

式(1)为作业σn在机器m上的完工时间,表示以最大完工时间最小化为指标,即Makespan指标;式(2)表示作业σ1在机器1上的完工时间;式(3)表示作业σi在机器1上的完工时间;式(4)表示作业σ1的第j道工序必须在第j-1道工序完成后才能开始;式(5)表示作业σi-1在机器j上加工完成后和作业σi在机器j-1上加工完成后,才能在机器j上加工作业σi。

2 Flow-Shop问题的求解

2.1 遗传算法求解Flow-Shop的基本流程

求解流水车间排序的遗传算法基本操作流程如图1所示。

第一步,随机生成N个个体组成初始种群。其中,每个个体表示染色体的基因编码,对于流水车间排序问题,编码方式可采用基于工件的自然数编码。例如,对于六个工件的流水车间排序问题,假设工件的加工顺序是j2,j1,j5,j4,j3,j6可将其编码为:[2,1,5,4,3,6]。

第二步,计算个体的适应度。根据遗传算法目标函数的方向与适应度值变化方向相同的原则,本次采用的适应度函数为式(1)目标函数的倒数。

第三步,选择操作。本文采用最常用的轮盘赌选择算子,用正比于个体适应度值的概率来选择个体,这种选择思想与生物的自然选择相似。

第四步,交叉操作。交叉操作是产生新个体的主要方法,在遗传算法中起着关键作用。本文将在下一节重点对其进行展开分析。

第五步,变异操作。常用的变异方法有交换变异、倒位变异、插入变异等。本文采用交换变异,即交换染色体中两个基因座的基因值。例如一条染色体为[2,1,5,4,3,6],随机选取两位置2和4,则交换变异的结果是[2,4,5,1,3,6]。

2.2 交叉算子

单点交叉操作与生物的交配重组类似,是基本遗传算法中最常用的交叉算子之一。基于工件编码的单点交叉操作为:1)随机产生一交叉点;2)把两父代交叉点前的基因复制给子代1和子代2;3)从父代2(或1)中找出子代1(或2)中缺少的基因,并将其按顺序复制给两子代交叉点后的位置。如图2所示。

从单点交叉的操作可以看出,单点交叉使子代继承了父代交叉点前的基因,交叉操作只改变了交叉点后的部分基因,而交叉点前的基因一直未发生改变。这样,交叉点前的基因进化只能依靠变异操作进行,而变异的概率是很低的,显然,进化对初始种群依赖较大,进化速度缓慢。

为了提高寻优质量,本文提出了一种改进型单点交叉算子。取上述父代1、父代2为例,其操作为:首先随机产生一交叉点;然后产生一个随机数(rand),确定交叉段的保留部分,具体如下。

1)当rand小于等于0.5时,把两父代交叉点后的基因复制给两子代,然后从父代2(或1)中找出子代1(或2)中缺失的基因,并将其按顺序复制给两子代交叉点前的位置。此时得到的子代如图3所示。

2)当rand大于0.5时,把两父代交叉点前的基因复制给子代,然后从父代2(或1)中找出子代1(或2)中缺失的基因,并将其按顺序复制给两子代交叉点后的位置。此时得到的子代同图2。

从以上操作可知,改进型单点交叉使父代交叉点前后的基因以等概率的方式得以保留或改变。与典型的单点交叉相比,改进型单点交叉有机会改变染色体中的所有基因,增加了染色体基因的改变范围,提高了基因重组的程度,这样就扩大了解空间的搜索范围,有利于提高进化速度和优解的质量。

3 实验分析

本实验用标准测试案例Car05[9](10×6,即10个工件,6台机器)作为测试数据,其标准最优解为7720。针对Flow-Shop问题,对上述两种交叉算子进行比较分析。在计算实验时,使用相同的初始种群,种群个数60,交叉率0.8,变异率0.05,进化代数300;分别独立运行50次。

图4、图5是这两种交叉算子具有一定代表性的进化历程图。

从图4、图5中的进化历程图可知,两种交叉算子的平均适应度曲线都有爬升趋势,表示种群整体都有往收敛方向搜索的趋势。

图5中,改进型单点交叉进化历程前期平均适应度曲线上升速度比单点交叉快;从平均适应度曲线可以看出,改进型单点交叉算子种群整体进化历程明显优于单点交叉算子。

得到的最好染色体是[5 4 2 1 3 8 6 10 9 7],最优解为7720,图6是相应的的甘特图,图中的数字表示工件。

表1中的数据为50次计算的统计结果。

由表中的数据可知,改进型单点交叉算子的全局收敛能力比单点交叉算子强,全局收敛次数是单点交叉算子的近2倍;且首次出现最优解的平均代数小,平均在130余代就能找到最优解。实验结果显示,改进型单点交叉算子性能良好。

本实验对大规模问题也进行了测试,种群整体进化趋势与上述案例类似,改进型单点交叉算子无论是求解质量还是寻优速度都优于单点交叉算子。限于篇幅,不再具体展开。

4 结论

本文提出一种改进型单点交叉算子,通过改变染色体基因的变化范围,扩大解空间的搜索范围,提高了进化速度。标准测试案例表明,改进型单点交叉算子在求解流水车间排序问题时的性能优于典型的单点交叉算子。

参考文献

[1]赵良辉,邓飞其.用于作业车间调度的模拟退火算法[J].制造业自动化,2006,28(3):10-13.

[2]王娟,王建.自适应蚁群算法在流水车间调度的应用[J].计算机应用与软件,2008,25(11):117-119.

[3]伊华伟,张秋余.求解置换Flow_shop调度问题的改进遗传算法[J].计算机工程与应用,2007,43(22):41-43.

[4]Bassem Jarboui,Mansour Eddaly,Patrick Siarry.A hybridgenetic algorithm for solving no-wait flowshop[J].Int JAdv Manuf Technol,2010,11.

[5]周明.遗传算法原理及应用[M].北京:国防工业出版社,1996.

[6]涂雪平,施灿涛,李铁克.求解置换流水车间调度问题的改进遗传算法[J].计算机工程与应用,2009,45(36):50-53.

[7]Rubén Ruiz,Concepción Maroto,Javier Alcaraz.Two newrobust genetic algorithms for the flow-shop schedulingproblem[J].OMEGA,2006,34:461-476.

[8]武维,管晓宏,卫军胡.Flow shop问题的嵌套分区优化调方法[J].控制理论与应用,2009,26(3):233-237.

变形算法 篇7

关键词：精确算法,谐波齿轮,柔轮,共轭,装配变形

0 引言

谐波齿轮传动是由美国学者Musser[1]于1955年提出、在薄壳弹性变形理论的基础上发展而来的传动技术。谐波齿轮通常由柔轮、刚轮和波发生器组成,装配前,柔轮的原始剖面呈圆形,柔轮和刚轮的周节相同,但柔轮的齿数比刚轮齿数略少;波发生器的形状要能保证柔轮齿和刚轮齿在柔轮最大变形区域全齿啮合,在柔轮最小变形区域完全脱开。波发生器装入柔轮后,其长轴两端的齿与刚轮齿完全啮合,短轴处的齿则完全脱开;长轴与短轴之间不同区段内的齿处于啮入或啮出的状态。在波发生器为主动、柔轮从动、刚轮固定的传动状态中,当波发生器转动时,波发生器迫使柔轮变形不断变换,使柔轮齿处于啮合-啮出-脱开-啮入-啮合的循环状态,从而使柔轮产生与波发生器相反方向的旋转。波发生器转过一周,柔轮转过两者的齿数差,从而获得减速传动。

Musser[1]用来确定谐波齿轮传动齿形的计算模型只考虑了波发生器作用下柔轮中性层变形曲线的径向位移,没有考虑柔轮中性层变形曲线的曲率变化引起的柔轮轮齿对称线的转动,并基于满足定传动比和使轮齿实现面接触以提高承载能力,采用直线齿形。但这种齿形不能保证柔轮和刚轮齿形啮合时相互共轭,无法得到最佳的啮合性能[2]。柔轮在波发生器作用下产生的弹性变形是谐波齿轮获得变速传动的前提,同时对谐波传动的运动和啮合起决定作用,有必要对谐波传动中柔轮的变形进行深入研究,建立更精确的齿廓共轭方法,以提高啮合性能和啮合精度。

柔轮的变形形状是以柔轮中性层变形曲线来表征的,称之为原始曲线[2]。在装配状态下,柔轮齿的位置除了径向和切向位移外,其对称线相对于径矢转过一个角度。谐波齿轮传动共轭齿廓的求解方法中[2,3,4,5,6,7],为了便于计算,对柔轮中性层变形曲线的切向位移引起的转角和齿廓对称线相对于径矢的转角采用近似算法进行求解。本文对现有柔轮径向变形的求解方法进行分析,提出了一种基于柔轮弹性变形的共轭精确算法:在已知径向位移的基础上,根据中性层变形曲线不伸长条件,计算输出端的转角;在轮齿啮合端,准确计算切向位移引起的转角和齿廓对称线相对于径矢的转角,利用迭代法计算变形后齿廓共轭位置。对常见五种波发生器[2,8]作用下轮齿的共轭区间和共轭齿廓进行了比较计算。

1 定义坐标系

假定刚轮固定(即刚轮齿形曲线$\tilde{G}$固定),固定坐标系OXY与刚轮固连,Y轴与刚轮齿槽对称线重合,原点O位于刚轮回转中心。动坐标系S1{O1x1y1}和S2{Ox2y2}分别与啮合端柔轮齿和波发生器固连,y1轴与柔轮齿的对称线重合,原点O1位于中性层变形曲线$\tilde{C}$上;y2轴与波发生器的长轴重合,原点O位于波发生器的回转中心。

谐波齿轮传动中,柔轮、刚轮和波发生器之间的坐标关系及其柔轮中性层变形曲线的各位移示意图见图1,图中所示位置为波发生器相对于Y轴逆时针旋转至φ2,柔轮非变形端(传动输出端)相对于Y轴顺时针旋转θE。根据谐波齿轮传动摩擦模型的运动学理论,φ2与θE之间满足

$\theta {}_{\text{E}}=\frac{z{}_2-z{}_1}{z{}_1}\phi {}_2$

式中,z1、z2分别为柔轮和刚轮的齿数。

柔轮的中性层变形曲线由波发生器的形状决定,变形曲线的自变量φ为波发生器长轴相对于柔轮未变形端转过的角度,以图1所示逆时针方向为正,从图1中可以看出,φ=φ2+θE。波发生器迫使柔轮啮合端中性层曲线产生径向位移w(φ)、切向位移v(φ)和相对于径矢的转角μ(φ),这些位移引起柔轮齿啮合曲线$\tilde{R}$所在的坐标系S1移动到图1所示位置,O1点相对于Y轴转过的角度为γ(φ),动坐标系S1的y1轴相对于Y轴转过的角度为ϕ。图1中的角度位移满足如下关系:

$\gamma =\theta {}_{\text{E}}+\frac{v}{\rho }=\phi {}_1-\phi {}_2$

$\phi =\phi {}_2+\theta {}_{\text{E}}=\frac{z{}_2}{z{}_1}\phi {}_2$

ϕ=μ+γ

2 共轭齿廓算法

在坐标系S1内,柔轮齿啮合曲线$\tilde{R}$以参数形式给出,即

$\begin{array}{l}x{}_1=x{}_1(u)\\ y{}_1=y{}_1(u)\end{array}\}(1)$

式中,u为齿廓上坐标点的位置参数,u=tanαk-tanα0;αk为任意点的压力角;α0为基准齿形角。

中性层变形曲线$\tilde{C}$极坐标方程为

$\begin{array}{l}\rho =r{}_{\text{m}}+w(\phi )\\ v(\phi )=-{\displaystyle \int }{}_0^{\phi }w(\phi )\text{d}\phi \\ \text{ϕ}=\mu (\phi )+\gamma (\phi {}_1)\end{array}\}(2)$

式中,ρ为中性层变形曲线$\tilde{C}$之极径(由波发生器相对位置φ确定);w(φ)为变形后柔轮中性层变形曲线上点的径向位移;v(φ)为变形后柔轮中性层变形曲线上点的切向位移;rm为变形前柔轮中性层圆曲线的半径。

2.1 变形后柔轮中性层变形曲线转角位移的近似计算

柔轮变形时,除轮齿产生径向和切向位移外,其对称轴线相对于中性层变形曲线径矢转过角度μ,文献[2]中,考虑到μ值较小(一般μ<2°),以及w(φ)较之rm要小得多,故

$\begin{array}{l}\mu =\mu (\phi )=\arctan (\dot{\rho }/\rho )\approx \dot{\rho }/\rho \approx \\ \frac{\dot{w}(\phi )}{r{}_{\text{m}}}=\frac{1}{r{}_{\text{m}}}\frac{\text{d}w(\phi )}{\text{d}\phi }(3)\end{array}$

第一个近似认为ρ的导数很小,用$\dot{\rho }/\rho$代替$\arctan (\dot{\rho }/\rho )$;第二个近似认为径矢ρ与未变形柔轮中面半径rm非常接近,用rm代替ρ。

根据中性层变形曲线不伸长条件,未变形端圆弧的弧长等于变形端变形后的弧长,即

$\begin{array}{l}r{}_{\text{m}}\phi ={\displaystyle \int }{}_0^{\phi {}_1}\sqrt{\rho {}^2+\dot{\rho }{}^2}\text{d}\phi {}_1\approx {\displaystyle \int }{}_0^{\phi {}_1}\rho \text{d}\phi {}_1={\displaystyle \int }{}_0^{\phi {}_1}r{}_{\text{m}}\text{d}\phi {}_1+\\ {\displaystyle \int }{}_0^{\phi {}_1}w(\phi )\text{d}\phi {}_1\approx r{}_{\text{m}}\phi {}_1+{\displaystyle \int }{}_0^{\phi }w(\phi )\text{d}\phi {}_1=r{}_{\text{m}}\phi {}_1-v(\phi )(4)\end{array}$

第一个近似认为,相比ρ本身,ρ的导数值很小,将其忽略;第二个近似认为φ和φ1差别很小,将其看作相等。由此可将φ1表示成φ的函数:

φ1≈φ+v(φ)/rm (5)

由此得出γ和切向位移v(φ)的近似关系:

$\gamma \approx \phi +\frac{v(\phi )}{r{}_{\text{m}}}-\phi {}_2=\frac{z{}_2-z{}_1}{z{}_2}\phi +\frac{v(\phi )}{r{}_{\text{m}}}(6)$

采用了近似算法,可大大简化数学计算。对一般精度的谐波齿轮传动,近似算法引起的误差不大,能满足工程需要。但是,对于高精度谐波传动,这种近似算法的偏差需要准确计算,才能为高精度、零回差的谐波齿轮传动设计提供理论指导。

2.2 变形后柔轮中性层变形曲线转角位移的精确算法

柔轮轮齿对称线相对于径矢转过角度μ,由微分几何可知:

$\mu (\phi )=-\arctan (\dot{\rho }/\rho )=-\arctan \frac{\dot{w}(\phi )}{r{}_{\text{m}}+w(\phi )}(7)$

根据中性层变形曲线不伸长条件,φ由下面积分关系确定:

$\phi ={\displaystyle \int }{}_0^{\phi {}_1}\sqrt{(1+\frac{w(\phi )}{r{}_{\text{m}}}){}^2+(\frac{\dot{w}(\phi )}{r{}_{\text{m}}}){}^2}\text{d}\phi =F(\phi {}_1)(8)$

则

$\gamma (\phi {}_1)=\phi {}_1-\phi {}_2=\phi {}_1-\frac{z{}_1}{z{}_2}F(\phi {}_1)(9)$

比较式(3)和式(7),式(7)直接解微分方程,精确求解柔轮轮齿对称线相对于径矢的转角。而在式(8)中,如果不采用式(5)的近似计算,计算复杂程度增加,求解φ1关于φ的表达式需求解积分方程,无法给出解析表达式。故采用变量代换的方法求解,因为柔轮装配变形使得两齿廓的实际共轭位置是在φ1处,共轭算法变为以φ1为自变量,求解φ1位置上共轭齿廓,所有参数均表示为φ1的函数。

φ关于φ1的导数为

$\frac{\text{d}\phi }{\text{d}\phi {}_1}=\sqrt{(1+\frac{w(\phi {}_1)}{r{}_{\text{m}}}){}^2+(\frac{\dot{w}(\phi {}_1)}{r{}_{\text{m}}}){}^2}(10)$

按照包络理论,与柔轮齿形曲线$\tilde{R}$共轭的刚轮齿形曲线$\tilde{G}$的表达式为

$\left[\begin{array}{c}x{}_2(u,\phi )\\ y{}_2(u,\phi )\\ 1\end{array}\right]=\left[\begin{array}{ccc}\cos \text{ϕ}& \sin \text{ϕ}& \rho \sin \gamma \\ -\sin \text{ϕ}& \cos \text{ϕ}& \rho \cos \gamma \\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}x{}_1(u)\\ y{}_1(u)\\ 1\end{array}\right](11)$

它满足方程

$\begin{array}{l}\frac{\partial x{}_2(u,\phi )}{\partial u}\frac{\partial y{}_2(u,\phi )}{\partial \phi }\frac{\text{d}\phi }{\text{d}\phi {}_1}-\\ \frac{\partial x{}_2(u,\phi )}{\partial \phi }\frac{\text{d}\phi }{\text{d}\phi {}_1}\frac{\partial y{}_2(u,\phi )}{\partial u}=0(12)\end{array}$

$\frac{\partial x{}_2}{\partial u}=\frac{\partial x{}_1}{\partial u}\cos \text{ϕ}+\frac{\partial y{}_1}{\partial u}\sin \text{ϕ}$

$\frac{\partial y{}_2}{\partial u}=-\frac{\partial x{}_1}{\partial u}\sin \text{ϕ}+\frac{\partial y{}_1}{\partial u}\cos \text{ϕ}$

$\frac{\partial x{}_2}{\partial \phi {}_1}=(-x{}_1\sin \text{ϕ}+y{}_1\cos \text{ϕ})\frac{\text{d}\text{ϕ}}{\text{d}\phi {}_1}+\frac{\text{d}\rho }{\text{d}\phi {}_1}\sin \gamma +\rho \cos \gamma \frac{\text{d}\gamma }{\text{d}\phi {}_1}$

$\frac{\text{d}\rho }{\text{d}\phi {}_1}=\dot{w}(\phi )\frac{\text{d}\phi }{\text{d}\phi {}_1}$

$\frac{\text{d}\gamma }{\text{d}\phi {}_1}=1-\frac{z{}_1}{z{}_2}\frac{\text{d}\phi }{\text{d}\phi {}_1}$

以上公式包含w(φ)及其一阶、二阶导数。柔轮中性层变形曲线的径向位移w(φ)随波发生器的形式不同而不同。

3 柔轮的中性层变形曲线的径向位移

文献[8]给出了四力作用下等效圆环变形后的径向位移w(φ)、切向位移v(φ)的计算公式。按照圆环理论,考虑柔轮筒对齿圈的加强作用,建立具有等效刚度的齿圈模型。通过求解二维齿圈弯曲模型,获得齿圈弯曲变形后的变形曲线。该解答用分段函数表示,在对称位置满足相应的对称变形条件,在滚轮的集中力作用处满足变形连续条件。利用等效圆环理论模型可以得到装配状态下的圆环的变形和齿圈应力计算公式,为啮合分析和应力强度计算提供依据。圆环理论简化了柔轮筒体的空间变形特征。

文献[2]给出了四力作用下柔轮中性层变形曲线的级数解,通过建立等厚度壳体模型的应力-应变关系,用能量法给出近似解。等厚度壳体的简化模型将柔轮简化为等厚度的薄圆柱壳,按照半无矩理论计算柔轮的变形和应力。壳体理论的变形和应力结果为准确了解柔轮筒体的空间变形特征和应力分布提供了极大的帮助。特别是,齿圈部位的壳体变形特征是齿圈变形后啮合计算的基本参数。等厚度壳体模型没有考虑齿圈局部增厚对变形的影响。

等效圆环理论模型和等厚度壳体的简化模型作为计算柔轮变形的等效模型,都可以对柔轮在装配状态下的变形和应力进行计算,在等效模型上,两者的相互影响按照等效刚度的办法,给齿圈增加一个等效长度,或者给壳体增加当量厚度来体现。等厚度壳体级数解有截断误差,运算速度比较慢,但应用级数解编程较简单。由于未考虑沿齿长方向的歪斜变形,所讨论的谐波齿轮啮合作为平面啮合问题研究,柔轮的变形特征为平面变形。本文中,柔轮的中性层变形曲线的径向位移采用分段函数表示的等效圆环理论解,利用分段函数可使得运算速度大大提高。

式(12)的柔轮中性层变形曲线的径向位移取:

对0≤φ≤β的区间

$\begin{array}{l}w{}_1=\frac{w{}_0}{A-4/\text{π}}(A\text{c}\text{o}\text{s}\phi +\phi \sin \beta \sin \phi -\frac{4}{\text{π}})\\ \dot{w}{}_1=\frac{w{}_0}{A-4/\text{π}}(-A\text{s}\text{i}\text{n}\phi +\phi \sin \beta \cos \phi )\\ \ddot{{\displaystyle w}}{}_1=\frac{w{}_0}{A-4/\text{π}}(-A\text{c}\text{o}\text{s}\phi -\phi \sin \beta \sin \phi )\end{array}\}(13)$

对β<φ≤π/2的区间

$\begin{array}{l}w=\frac{w{}_0}{A-4/\text{π}}[B\text{s}\text{i}\text{n}\phi +(\frac{\text{π}}{2}-\phi )\cos \beta \cos \phi -\frac{4}{\text{π}}]\\ \dot{w}{}_2=\frac{w{}_0}{A-4/\text{π}}[B\text{c}\text{o}\text{s}\phi -\cos \beta \cos \phi -(\frac{\text{π}}{2}-\phi )\cos \beta \sin \phi ]\\ \ddot{{\displaystyle w}}{}_2=\frac{w{}_0}{A-4/\text{π}}[-B\text{s}\text{i}\text{n}\phi +\cos \beta \sin \phi +\cos \beta \sin \phi -\\ (\frac{\text{π}}{2}-\phi )\cos \beta \cos \phi ]\end{array}\}(14)$

$A=\sin \beta +(\frac{\text{π}}{2}-\beta )\cos \beta B=\cos \beta +\beta \sin \beta$

式中,β为集中力作用方向与变形长轴之间的夹角;w0为柔轮的中性层变形曲线的最大径向位移。

通过求解满足式(12)的φ1,以φ1为变量,由式(8)确定出φ的值及相应的共轭齿廓。

4 计算实例

算例一

谐波齿轮参数如下:z1=200,z2=202,模数m=0.5mm,α0=20°,柔轮变位系数x1=3.0,β=30°。选择渐开线齿廓

$\{\begin{array}{l}x{}_1=r{}_1[-\sin (u-\theta )+u\text{c}\text{o}\text{s}\alpha {}_0\cos (u-\theta +\alpha {}_0)]\\ y{}_1=r{}_1[\cos (u-\theta )+u\text{c}\text{o}\text{s}\alpha {}_0\sin (u-\theta +\alpha {}_0)]-r{}_{\text{m}}\end{array}$

θ=s/(2r1)=e/(2r1)

式中,r1为柔轮分度圆半径;s(或e)为变位后齿轮分度圆齿厚(或槽宽),s=e=(π/2±Δ)m;Δ为分度圆齿厚改变系数。

用本文算法和近似算法计算w0分别取0.9m、1.0m、1.1m和1.2m时的共轭区间。

近似算法和本文算法给出四滚轮波发生器作用下的结果如表1所示。

从表1看出,w0取不同值,本文算法计算出的共轭区间都偏移了一个角度。w0取0.9m,最小偏移角度0.362 02°,共轭区间大小4.782 93°。定义偏移角度与共轭区间大小比值的百分比为偏差,则偏差7.569%;w0取1.0m,最大偏移角度0.491 286°,共轭区间大小4.654 922°,偏差10.554%。

算例二

与算例一参数相同的谐波齿轮,计算在常见的五种波发生器[2,8]作用下的共轭区间,两种算法的结果如图2所示。从图2可看出,不同波发生器作用下,共轭区间几乎都平移了一个角度。

1.四滚轮(本文算法) 2.四滚轮(近似算法) 3.余弦凸轮(本文算法) 4.双圆盘(本文算法) 5.标准椭圆(本文算法) 6.余弦凸轮(近似算法) 7.双圆盘(近似算法) 8.标准椭圆(近似算法) 9.双滚轮(本文算法) 10.双滚轮(近似算法)

四滚轮波发生器作用下,两种方法结果的差别随φ角的增大而减小,最大偏移角度0.491 66°,共轭区间大小4.782 93°,偏差10.279%,最小偏移角度0.235 10°,偏差4.915%。双滚轮波发生器作用下,两种方法结果的差别随φ角的增大略有减小,最大偏移角度0.344 87°,共轭区间大小4.361 00°,偏差7.908%,最小偏移角度0.304 07°,偏差6.972%。

双圆盘、余弦凸轮、标准椭圆等三种波发生器作用下,两种算法的结果差别随φ角的增大而增大;双圆盘波发生器作用下,最小偏移角度0.420 72°,共轭区间大小4.674 20°,偏差9.001%,最大偏移角度0.803 69°,偏差17.194%。余弦凸轮波发生器作用下,最小偏移角度0.538 53°,共轭区间大小4.571 81°,偏差11.779%,最大偏移角度0.778 64°,偏差17.031%。标准椭圆波发生器作用下,最小偏移角度0.525 72°,共轭区间大小4.565 09°,偏差11.516%,最大偏移角度0.774 09°,偏差16.957%。

总体上看,两种算法的绝对偏差介于0.21°～0.81°之间;由于共轭区间较小(介于4°～5°之间),相对偏差比较大,介于4%～18%之间。

图3所示是五种波发生器作用下的共轭齿廓的计算结果。从图3看出,两种算法结果差别很小。

算例三

谐波齿轮参数如下:z1=140,z2=142,m=0.2mm,α0=20°,w0=1.0m,x1=2.13,刚轮变位系数x2=1.925,柔轮壁厚δ=0.3mm。计算四滚轮波发生器作用下谐波齿轮的共轭区间和共轭齿廓。

图4所示是共轭区间计算结果。可以看出,四滚轮波谐波齿轮的两种结果的差别随φ角的增大而减小,最大偏移角度0.634 33°,共轭区间大小0.856 27°,偏差74.081%;最小偏移角度0.593 97°,偏差69.367%。

该谐波齿轮主要用于传递运动,模数很小,共轭区间范围不到1°,所以偏差很大。图5所示是共轭齿廓计算结果。可以看出,两种算法结果差别很小。

5 结论

谐波齿轮传动的正确啮合只可能发生在共轭区间内。渐开线谐波齿轮的共轭区间在长轴区域比较小的区间内分布;对于传递动力的谐波齿轮,因受载变形,可使共轭区间外的轮齿也参与啮合,但应为尖点啮合。由本文实例看出,近似算法得出的共轭区间比精确算法都有一个提前量。对于小模数的传递运动的谐波齿轮,近似算法引起的偏差很大,不能忽略。从计算实例得出:

(1)不同w0值引起的偏差不同。当w0等于1倍模数时,偏差最大;随齿数和模数的减小,共轭区间减小,偏差增大;齿廓对称线相对于径矢的转角的偏差较小;切向位移引起的转角的偏差较大。

(2)不同波发生器作用下的中性层变形曲线引起的偏差不同。四滚轮、双滚轮波发生器作用下的共轭区间,两种方法的差别随着φ角的增大而减小;双圆盘、余弦凸轮、标准椭圆波发生器作用下的共轭区间,两种方法的差别随着φ角的增大而增大;总体来看,近似算法对共轭区间影响较大,对形成的共轭齿廓影响很小。

本文算法也适用于对波发生器轮廓线的设计计算。采用柔轮中性层曲线的等距曲线作为波发生器轮廓线的形状,在波发生器轮廓计算中可采用本文算法精确计算切向位移引起的转角。

参考文献

[1]Musser C W.Strain Wave Gearing:US,2906143[P].1959-09-29.

[2]沈允文,叶庆泰.谐波齿轮传动的理论和设计[M].北京:机械工业出版社,1985.

[3]张丽华.基于有限元分析的谐波齿轮传动变形协调研究[D].大连:大连理工大学,2004

[4]刘书海.基于运动几何学的谐波齿轮传动齿形的研究[D].大连:大连理工大学,2003.

[5]辛洪兵.研究谐波齿轮传动啮合原理的一种新方法[J].中国机械工程,2002,13(3):181-184.

[6]辛洪兵,何惠阳.谐波齿轮传动共轭齿廓的研究[J].长春光学精密机械学院学报,1996,19(2):22-26.

[7]范元勋,王华坤,宋德锋.谐波传动共轭齿廓的运动学仿真研究[J].南京航空航天大学学报,2002,34(5):447-450.