基于改进希尔伯特振动分解的机械故障诊断方法研究
教务系统
JOURNAL OFVIBRATION AND SHOCK基于改进希尔伯特振动分解的机械故障诊断方法研究唐贵基、庞彬彬(华北电力大学机械工程系，河北保定071003 )。
摘要：针对多分量机器故障振动信号的特征提取问题，介绍一种基于希尔伯特振动分解(HVD )的时频分析方法法律。 该方法首先利用Hilbert变换获得原始振动信号的分析信号，然后通过分析信号的瞬时频率低通滤波器获得信号中幅度最大分量的瞬时频率同时通过同步检测得到对应的瞬时幅度和初始相位，最后通过迭代运算自适应地检测原信号的时频信息。 针对HVD方法的边界效应问题，提出了一种基于相关系数准则的波形匹配边界延拓法，并对其进行了改进。通过两组伪信号分析验证了HVD方法对多分量非平稳信号的分解能力，同时表明改进的HVD方法能较好地抑制边界效果。 给出了子系统间油膜涡动故障诊断的实例，验证了该方法的工程实用性。关键词：希尔伯特振动分解； 多成分信号； 时频分析； 波形匹配延拓； 机械故障诊断
中图分类号： TH133.3； TH17文献标识代码： adoi:10.13465/j.cn ki.JVs.2015.03.027amechanicalfaultdiagnosismethodbasedonimprovedhilbertvibrationdecompositionG Gui—ii. NGBin
(SchoolofMechanicalEngineering，noahchinaelectricpoweruniversity，Baoding071003，China ) )Abstract: A time一至equencyapproachbasedonhilbertvibrationdecomposition (hvd ) methodwasintroducedinordertoextractfaultfeaturesofmulti—componentmechanicalfaultvibrationsignalsaccurately.firstly，the analyticalsignalsoftheoriginalvibrationsignalswereobtainedthroughhilberttransformation.secondly，the instantaneousstationaryfrequencyofthelargestamplitudecomponentwasachievedusingalow-passfilteringofanalyticalsignals’instantaneousfrequencies，thecorespondingamplitudeandinitialphasewerealsoestimatedaccordingtothesynchronousdetecting，then the time一至equencyinformationofeachcomponentoftheoriginalsignalwasdetectedadaptivelywithiteration computation.aimingattheendeffectsofhvd，awavematchingextendingmethodbasedoncorelationcoeficientcriteriawasproposedtoimprovehvd.theanalysisoftwogroupsofsimulatedsignalshowedagoodcapacityofhvdindecomposingthenon—stationary multi—component signals，andtheresultsshowedthattheimprovedhvdsuppressesitsendeffects.finally，afaultdiagnosisinstanceofanoilwhirlofarotorsystemwasgiventovalidatethefeasibilityofthismethod。keywords:hilbertvibrationdecomposition； 多组件信号； time一至~equency analysis； wave matching扩展程序； 机械故障诊断
机械类零件发生故障时，振动信号通常为多点测量瞬态信号，采用适当的时频分析方法提取振动信号号码中各分量的时频分布信息是故障特征提取和故障类别判断的关键(j )。 机械故障诊断常用时频分析方法为Wigner-Ville分布、小波分析、希尔伯特黄变换(HilbertHuangTransforln，HHT )等。 其中Wigner—Ville分布理论上具有最高的时频聚集性，但处理较多对于分量信号，存在严重的交叉项干扰问题j。 小波分析没有足够的适应性，一旦变小波基函数选定 ，那 么所有
收稿 日期 ：2014一O1—02 修 改稿 收到 Et期 ：2014—02—03第一作者 唐贵基 男 ，教授 ，1962年生
数据都必须用 此小波函数分析 ，易造成信号 的能量泄漏 ，使定量 的时频分析变得 困难 J。HHT由 EMD和希尔伯特变换两部分组成 ，方法简单 ，应用较广。但 EMD理论上具有模 态混叠 、端点效应 、过包络 、欠 包络等缺陷 ，很大程度上影响了 HHT的分析精度 。希尔伯特振动分解是 Feldman等 在总结 上述 时频分析方法基础上提出的一种非平稳信号分析的新方法。该方法首先利用 Hilbea变换获得原始信号的解析信号 ，然后通过低通滤波对原信 号的非平稳频率成分进行分离 ，获得信号中幅值最大的分量 的瞬时频率 ，并经同步检测得到该分量 的瞬时 幅值 和初相位 ，最后通过迭代运 算 自适 应检 测 出原信 号各个 分量 的时频信振 动 与 冲 击 2015年第 34卷息。HVD和 HHT方法都 以希尔伯 特变换 为基础 ，但HVD方法避免了 HHT方法中复杂的 EMD过程 ，对 多分量非平稳信号进行分析 ，具有分解精度高 ，分解 速度快的优点。针对 HVD的边界效应 问题 ，提 出一种基 于相关系数 准则 的波形 匹配边 界延拓 法来对其 进行改进 。应用该方法对仿真信 号和转子系统油膜涡动实验测试信号进行分析 ，分析结果验证 了 HVD方法对非平稳信号良好 的时频分析效果 ，通过对 比分析表 明改进后 的 HVD方法消除了自身 的端点 问题 ，使得分析结果更加精确 。1 HVD 的基本原理1．1 Hilbert变换
任一非平稳 连续时间信号 (t)的 Hilbe~变换可表 示为 ]
()：研 ()]=7r Jr一* t — d丁 (1)(t)和 (t)构成共轭复数对 ，从而得到 (t)的解析信号z(t)= (t)+j(t)=A(t)e‘ (2)式 中：A(t)和 (t)分别对应 (t)的瞬时幅值和瞬时相位 ，表达式分别为A(t)=√ (t)+矿(t) (3))… ctan (4)
对应瞬时频率表达式为f(t)= 1 (5)
1．2 幅值最大分量的频率估计多分量非平稳信号 (t)可表示为¨叫()=∑n()cos(2盯 (t)dt+O1)=∑口。()cos( ()) (6)
式 中：／、r为 (t)中的分量个数 ，02(t) (t)、0、 (t)分别表示第 Z个分量 的瞬时幅值 、瞬时频率 、初相位和相位 。则经 Hilbe~变换 ，(t)的解析信号为z(f)=∑0f(￡)e (7)参考文献 [10]，(t)的瞬时幅值 A(t)可表示为4(t) = ( l(t)+A2(t)+… +AⅣ(t))丁 =Ⅳ ．
(∑A(￡))寺 (8)其中：
A(￡)=∑ (￡) (￡)cos( (￡)一 ()) (9)式 中：m=1，2，…Ⅳ同时 (t)的瞬时频率 t)可表示为Ⅳ
∑[()一 ()]A()
)= ()+坚 — _ 一 (10)
上式表明 f(t)可表 示 为两 部分 ，其 中第 1部 分(t)为幅值最大分量 的瞬 时频率 ，第 2部分为相对于()快速变化的振荡频率部分。因而实际应用 中可通过低通滤波滤除 f(t)的高频振荡频率 ，得到幅值最大分量 的瞬时频率 (t)。本文采用 FIR低通滤波器 ，为了实现不同瞬时频率成分的完全分离 ，截止频率 需小于相邻次分量瞬时频率差值绝对值的最小值 ，则 越小越能实现不同分量 的准确检测 。通常 设 为采样频率的0．01—0．05倍时 即可完全滤 除厂(t)的高频振荡频率 。实际分析 中在考虑滤波器稳定性 的前提 下 ，若需要设定更小的截止频率 ，可采用降采样的方法实现。1．3 同步检测求幅值和相位将以上估计的瞬时频率看作参考频率 (t)，构造
两正交信号 COS(2~rffr(t)dt)、sin(2"rffr(t)dt)。将(t)与 cos(2盯J (t)dt)相乘得下面表达式Ⅳ ．
。 (f)=∑ (f()COS(2丌 ()dt+r 1
icve智慧职教官网
Z1(t)=÷。()cos0 (13)1
之(t)=一÷0，(t)sin0， (14)从而同步检测得到瞬时频率为 (t)的分量的瞬时幅值 口，(t)和初相位 为 ：n(t)=2√((f))+(之(f)) (15)： 一 arcta 盟 (t6)
(t)
1．4 迭代运算
根据上述步骤提取 出 (t)中幅值最大的分量
1(t)=01(t)cos(2~r (t)dt+01) (17)将 (t)与 (t)的差作为初始信号 ，即Ⅳ一1= (t)一 1(t) (18)
重复 1．2和 1．3节 的步骤获得 (t)。通过迭代
第 3期 唐 贵基等 ：基于改进 的希 尔伯 特振动分解 的机械故 障诊 断方 法研 究 169运算依次获得其他分量 的时频信息 ，将式 (18)的归一化标准差作 为迭代 终止条件 ，当 or<0．001时迭代停止 。
2 基 于相 关 系数准则 的波形 匹配边界 延拓法HVD算法以 Hilbert变换为基础 ，并且在对各分量的频率 、幅值 、和相位估计 过程 中应用到 了低通滤波。由于 Hilbert变换和低通滤波均存在边界效应 ，使得 当前分量的检测产生误差 。并且 随着迭代运算的不断累积 ，误差会越来越 大 ，严重影 响 了 HVD算 法 的检测精度。 目前 ，解决边界效应 的一种有效方法 是对信号进行边界延拓 。为降低数据延拓 的偶然性误差 ，本文提 出一种基于相关系数准则的波形匹配边界延拓新方法 。具体步骤如下 ：(1) 假 设 样 本 数 量 为 11,的离 散 信 号 (t)为{(t)，(t)，(t)，… ，(t)}，共有 m个极大值点和P个极小值点 ，其中信号采样间隔At=t—t 。极大值
点 ：indmax={ (1)， (2)，… ，M(m)}，对 应时 间为t = {taxl，tm越，…，t }。极 小 值 点 ：indmin={Ⅳ (1)，Ⅳ (2)，…，N (P)}，对应 时 间为 ti={ti，tin2，… ，tminp}。(2)如 图 1所示 ，以左边界第一极值点 为极 大值点为例 ，选取左端点到第 2个极值点间的波形 {(t)，…
， (t… )，…，(t )}为研究对象 ，记为 W 。j型馨
o tjs
图 1 波形 匹配示 意图
Fig．1W aveform matchingschematicdiagram(3)依次以 M(2)，… ， (m)为中心定义与 W 同样 长度 的子 波 形 W(iI>2)，W ={(t i一(t 一t。))，…
， (t…i)，…， (t +(ti。一t 1))}。分另0求 W1与 W的相关系数值 作为两个波形的匹配系数 ，其 中
g ：= ———————■■—————————一 (‘1l9)， D 1式 中：E[·]表示数学期望 ， ．和 ， 分别 为波形 W 和W 的均值 ， ．和 分别为原信号 W 和 W 的标准差。(4)取 昭 绝对值最大的子波形 W为 W 的最佳匹配波形 ，将 W前 的信号数据延拓 到 (t)左侧 ，延拓点数根据需要选择。(5)按照同样 的方式对信号进行右边界延拓。(6)对延拓后的信号进行 HVD分析 ，根据原始信号的长度及在原始信号在延拓信号中对应 的位置截取分析结果 。
3 仿真信 号分析
仿真信号 1：考察一两分量调频信号
1 (t)= sin(2~rlOOt+COS(27r5t))+0．5sin(2,rr2OOt+COS(2-rr4Ot)) (20)式中：t∈[0，0．25]，信号采样频率fs=4096Hz。图 2为信号经 EMD分析得到的前两个 IMF分量 ，两个分 量分别 出现了一定程度 的模态混叠和幅值失真现象。图3(a)为 (t)直接用 HVD方法分析得到的两个分量 ，图 2 (t)的 EMD分析结果Fig．2EMDanalysisresultsof l(t)2
O
一 2
篓o．s0． O．5l罂 0
0 5 lO 15 2O 25t× lO ／s
(a) 。(f)基于HVD的检 测分量(b) X1(f)基 Zf-HVD的检测各分量时频分 布三维 图图 3 ()直接进行 HVD分析结果
Fig．3HVDanalysisof I(#)withoutimproving图 3(b)为对应 的时频分 布三维 图。从分析结果 中看出 HVD方法克服了 EMD方法的模态混叠和幅值失真问题 ，但存在 明显 的边界效应 问题 。并且 随着 分解过程的继续 ，第 2个分 量的端点效应 比第 1个分量更 为明显。图 4(a)、图4(b)分别为本文改进 的 HVD方法对 (t)分析得到的分量和时频分布三维 图，两个分量的边界效应得到抑制 ，使得时频分析结果更加准确 。仿真信号 2：考察一两分量调幅调频信号170 振 动 与 冲 击 2015年第 34卷2 (￡)= (2+0．5cos(10wt))sin(300'rrt+cos(20~t))+(1+0．4cos(20,rrt))sin(900~t+cos(40~t)) (21)分析时间 t∈[0，1]，信号采样频率 =2048Hz。
图 5、图 6分别 为直接用 HVD方法和改进 HVD方法对：(t)的分析结果 。可见 HVD方法有效提取了 ：(t)两个分量的时频信息 ，同时经过改进有效克服 了 HVD方2O
逼-2馨 lO1
全国职教大会
(a) ．(f)基于改进的HVD的检测 分量(b) x【f)基于改进 的HVD的各分量时频分布三维 图
图 4 (t)基于改进 的 HVD分析结果Fig．4AnMysisof l(t)basedon improved HVD4 油膜涡动故 障诊断 实例O． 5一些 20
法 自身的边界效应 。
仿真信号 以两组多分量调制信号为例验证 了 HVD方法对多分量非平稳信号 的时频 分析能力。齿 轮、轴承 、转子等旋转机械元件 的故 障振动信 号通 常为多分量非平稳信号 ，HVD方法为此类机械元件的故 障特征提取和故障类型判断提供一种新的途径 。0 0．2 0．4 0．6 0．8 1．0t／s(a) ，(f)基于HVD的检测分 量(b) ，(『)基于HVD的各分量时频分布三维图
图 5 (t)直接进行 HVD分析结果Fig．5HVDanalysisof 2(t)withoutimproving
利用 BentlyRK一4转子实验台模拟了 2800r／rain转速下转子系统 油膜涡动故 障。实验数据 由美 国 Io—tech公 司生产 的 ZonicBook／618E信号采集设 备采得 ，采样频率为 1280Hz，采样点数为 1024。油膜 涡动 是转 子轴 颈 在 轴 承运 动过 程 由于 受到 油膜力的作用 ，转子绕 自身轴心旋转 ，同时轴心又绕轴承中心连线 回转的一种运动形式 ，涡动频率约为转速 的0．47倍图 7 油膜 涡动故障信号时域波形
Fig．7 Timedomainwaveform ofoilwhirlt／s
(a) (f)基 于改进 的HVD的检测分量(b) ，(r)基于改进的HVD的各分量 时频 分布三维 图
图 6 (￡)基 于改进的 HVD分析结果Fig．6Analysisof 2(t)based on improved HVD图 7为转子油 膜涡动故 障信 号 的时域波形 。图8(a)为油膜涡动故障信号运用改进的 HVD方法分解得到的两个分量 ，第一个分量对应涡动信号 ，第二个分量对应基频信号 。图 8(b)为对应的时频分布三维 图，两个分量的瞬时频率分别为 23．25Hz、46．65Hz，涡动频率近似为基频的一半 。可见本文方法有效提取 了实验转 子系统油膜涡动 的故 障特征。图 9为直接 HVD方法的分析结果 ，虽然也提取到 了涡动信 号和基 频信号两个分量 ，体现 出油膜涡动的故障特征 ，但端点效应比较明显 ，极大影响了 HVD方法的分析精度。为突 出本文方法 的优点，同时对油膜涡动故障信号进行 了EMD分解 ，结果如图 10所示 ，产生 6个 IMF分量。其中 IMF1、IMF2分别对应基频信号和涡动信号。与改进的 HVD方法相比，两个包含故障特征信 息的 IMF分量均出现了一定程度 的模 态混叠 和幅值 失真现象 ，代表的物理过程被干扰 ，不能体现设 备的真实状态 。对 比分析表明改进的 HVD方法分解效率 高、分解速度快 。提取 的故障信息更加直观 、准确 。] ¨ ] _]～一嘲一～圈一蝴一一一 一 ～ ～0
第 3期 唐 贵基等 ：基于改进 的希尔伯特振动 分解 的机械故障诊断方法研究 171>一2
t／s
(a) 改进 的HVD方法的检测分 量(b)改进的HVD方法的检测分量的时频分布三维图
图 8 改进 的 HVD方法 的分析结果
Fig．8 Improved HVD analysis圣一罂88
o蝴 蝴40L — —
- 0● 2— — 0●4— — 0● 6 0● 8t／s
(a)直接HVD方法的检测分量(b)直接HVD方法的检 测分量的时频分布三维 图图 9 直 接进 行 HVD分 析结果8
Fig．9 HVD analysiswithoutimproving至塞 。堇．5Or_——————————————巫——————]_50E二二二二二 二 ==二E三
E三三三 5． 5
O 0．2 0．4 0．6 0．8t／s图 10 油膜涡动信号的 EMD分析结果
Fig．10 AnalysisresultsofEMD foroilwhirlsignal5 结 论本文介绍了一种多分量非平稳信号的时频分析新
方法——希尔伯特振动分解(HVD)。该方法 以 Hilbert变换为基础 ，可以 自适应地将原始振动 信号的各个分量分解出来 ，通过检测各个分量 的频率 、幅值和相位信息 ，提取振动信号 的故 障特征 ，具有分解精度 高、分解速度快 的优点。针对 HVD方法 自身的边界效应 问题 ，提 出一种基于相关系数准则 的波形匹配边界延拓法对HVD方法进行改进 。通过仿真信号分析验证 了 HVD方法对多分量信号 的分解 能力 和本 文方法对 HVD边界效应 的改进效果。最后给出转子系统油膜涡动故障诊断实例 ，取得很好 的分析效果 ，为机械故 障诊断提供一 种新 的途径 。参 考 文 献
[1]程军圣 ，杨宇 ，于德介．基 于广 义解 调时频分析的多分量信号分解方法[J]．振动工程学报，2006，19(6)：563—568．CHENG Jun—sheng， YANG Yu， YU De-jie． A multi—component signal decomposition method based on thegeneralized demodulation timefrequency analysis 『J]．JournalofVibrationEngineering，2006，19(6)：563—568．[2]来五星，轩建平，史铁林 ，等．Wigner—Ville时频分布研究及其在齿轮故 障诊 断 中的应 用 [J]．振 动 工程 学报 ，2003，16(2)：247—250．LAIW u—xing，XUAN Jian—ping，SHITie-lin，eta1．ResearchofW igner—Ville time frequency and application in detectinggear pinion fault[J]． JournalofVibration Engineering，2003，16(2)：247—250．
[3]张晓峰 ，李 功燕 ．应用小 波分析 提取故 障诊 断信号 的特定频段[J]．振动与冲击，2004，23(4)：47—50．
ZHANG Xiao—feng， LI Gong—yan． Review on structuralcontrolforoffshoreplatforms[J]．JournalofVibrationandShock，2004，23(4)：47—50．
[4]HuangN E，Shen Z， Long S R．The empiricalmodedecompositionandtheHilbertspectru m fornonlinearandnon—stationarytimeseriesanalysis[J]．Proc．R．Soc．Lond．A，1998，454：903—995．[5]郭琦 ，刘 卜瑜 ，史立波 ，等．基于二次 EEMD的 Wigner．Ville分布旋转机械故障信号分析及试验研究[J]．振动与冲
击 ，2012，31(13)：129—133．
GUOQi，LIU Bo—yu，SHILi—bo，eta1．Experimentalstudyand faultsignalsanalysisofmachinery based on dual EEMDandWigner．Villedistribution 『J]．Journal ofVibration andShock，2012，31(13)：129—133．[6]向玲 ，唐 贵基 ，胡爱军 ．旋转机械非平稳振动信号 的时频 分析 比较 [J]．振动与 冲击 ，2010，29(2)：42—45．
XIANG Ling，TANG Gui-ji，HU Ai-jun．Vibrationsignal Stime—frequency analysis and comparison for a rotatingmachinery [J]．Joumal ofVibration andShock，2010，29(2)：42—45．(下转第 182页)>里 覃l馨
182 振 动 与 冲 击 2015年第 34卷2001，21(2)：1105—1110．
[9]孔明．半速涡动 引起 的振动故 障诊 断及处理 [C]／／第三届安徽 自然 科 学 学术 年 会 安 徽 省 电机 工 程 学 会 论 文 集 ，2005．27—30．[10]朱瑜 ，张朋波 ，王雪 ．转 子系统 油膜 涡动及 油膜 振荡 故 障特征分析 [J]．汽轮机技术 ，2012，54(4)：306—308．
ZHU Yu， ZHANG Peng—bo， W ANG Xue． Fault featureanalysisofoilwhirland oilwhip forrotor system 『J]．TurbineTechnology，2012，54(4)：306—308．[11]唐贵基 ，向玲，朱永利．基于 HHT的旋转机械油膜涡动和油膜振荡故障特征分析 [J]．中国电机工程学报 ，2008，28(2)：77—81．TANGGui-ji，XIANGLing，ZHUYong-li．Faultanalysisofoilwhirland oilwhip based on Hilbert—Huang transform forrotorsystem 『J1． Proceedings ofthe Chinese Society forElectricalEngineering，2008，28(2)：77～81．[12]胡 爱军 ，朱瑜 ．基于阶 比分析 的转子 系统油 膜涡动故 障诊断研究[J]．汽轮机技术 ，2012，54(3)：223—225．HU Ai-jun．ZHUYu．Analysisofoilwhirldiagnosisbasedonorderanalysisforrotorsystem [J]．TurbineTechnology，2012，54(3)：223—225．
[13]孔庆鹏，宋开臣，陈鹰．发动机变速阶段振动信号时频分
析 阶比跟踪 研 究 [J]．振 动 工程 学报 ，2005，18(4)：448— 452．
LUO Song—rong，CHENG Jun—sheng，YANG Yu． Machinefault diagnosis method using ITD and variable predictivemodel—based classdiscrimination lJ1．Journal ofVibrationandShock．2013，32(13)：43—48．8 ] Feldman M． Time—varying vibration decomposition andanalysisbasedontheHilbe~transfoITnlJ1．JournalofSoundand Vibration．2006．295：518—530．[9]刘慧，刘 国海，沈跃．采用希尔伯特振动分解的非整数次谐
波检测新方法 [J]．高 电压技 术 ，2009，35(7)：1758—1762．LIU Hui．LIU Guo—ha1．SHEN Yue． Novelmethod for non2005，18(4)：448—452．[14]赵 晓平 ，侯荣涛．基 于 Viterbi算 法的 Gabor阶 比跟踪技 术[J]．机械工程学报 ，2009，45(11)：247—252．
ZHAO Xiao—ping， HOU Ran g—tao． Gabor order trackingbased on Viterbi algoithm [J]． Chinese Journal ofMechanicalEngineering，2009，45(11)：247—252．[15]金阳，郝志勇，郑旭，等．基于 Gabor重构的发动机振声信号阶次分 量 瞬 时 幅值 的 提 取 [J]．振 动 与 冲 击 ，2011，30(9)：15—20．JIN Yang，HAO Zhi—yong，ZHENG Xu，eta1．Extractionofinstantaneous amplitude of order components in enginevibroacoustic signals based on Gabor．reconstruction fJ1．JournalofVibrationandShock，2011，30(9)：15—20．[16] 申永军 ，张光明 ，杨绍普 ，等．基于 Gabor变换 的盲 信号分离方法[J]．振动与冲击 ，2010，29(10)：166—169．SHEN Yong-jun，ZHANG Guang—ming，YANG Shao—pu，eta1．New method ofblind source separation based on GabortransformationlJ】．Journal ofVibrationandShock，2010，29(10)：166—169．
[17]BastiaansM J，GeilenM C．OnthediscreteGabortransformandthe discrete Zak transform [J]． Signal Processing，1996，49(3)：151—166．
[18] SondergaardPL．EficientalgorithmsforthediscreteGabortransfoFinwith alongFIR window 『J1．JournalofFourierAnalysisand Application，2012，18：456～470．[19]张贤达．现代信号处理(第二版)[M]．北京：清华大学出版社 ，2002．[10][11]
integer harmonics measurement using Hilbert vibrationdecomposition[J]．HighVoltageEngineer，2009，35(7)：1758—1762．
Feldman M． Theoretical analysis and comparison of theHilbea transforln decomposition metl1ods『J]． MechanicalSystemsan d Signal Processing，2008，22：509—519．向玲 ，朱永利，唐贵基．HHT方法在转子振动故障诊断中的应用 [J]．中国电机工程学报 ，2007，27(35)：84—88．
XIANG Ling，ZHU Yong—li，TANG Gui-ji． Application ofHilbert—HuangTransform method in vibration faultsdiagn osisforrotorsystem [J]．ProceedingsoftheCSEE，2007，27(35)：84—88．