袁 幸,朱永生,張優(yōu)云,洪 軍

(1.西安交通大學(xué)潤滑理論及軸承研究所,陜西 西安 710049;2.西安交通大學(xué)機械制造系統(tǒng)工程國家重點實驗室,陜西西安 710049)

引 言

飛機、電力機車、高檔機床等重大裝備的核心技術(shù)問題之一就是軸承,它對整個裝備制造業(yè)的發(fā)展水平有著舉足輕重的作用,其運行狀態(tài)往往直接影響整臺機器的精度、可靠性及壽命[1]。 Sawalhi N將最小熵盲反卷積和譜峭度應(yīng)用于滾動軸承故障檢測中[2]。Karthik K應(yīng)用互信息進行了故障特征選取[3]。譚繼勇針對滾動軸承沖擊信號的檢測問題,構(gòu)造了隨機共振自適應(yīng)檢測算法,提高了信噪比,避免了漏峰現(xiàn)象[4]。文獻[5]提出了基于經(jīng)驗?zāi)B(tài)分解和奇異值差分譜的軸承故障診斷方法。從已有研究結(jié)果來看,以先進信號分析為基礎(chǔ)的滾動軸承診斷技術(shù)研究較多,人們大都針對損傷信號處理與分類,研究相關(guān)的算法,卻很少關(guān)注服役軸承具體的損傷尺寸。對于滾動軸承這樣的基礎(chǔ)部件,應(yīng)用廣泛但規(guī)格種類繁多,工況差別大,不同損傷程度的樣本獲取困難。以機車為例,輪對軸承的故障樣本少,一個輪對檢修流水線在3個月時間檢測的1 200多個輪對軸承中僅有2例故障軸承[6]。再者,單純的基于信號處理的傳統(tǒng)方法較少考慮軸承系統(tǒng)的具體結(jié)構(gòu)特點,側(cè)重于對故障存在性的檢測,不易描述損傷演化過程,因而難以準確評價軸承可以持續(xù)運行的安全余量。

本文系統(tǒng)地提出根據(jù)正、反問題進行損傷量化診斷的原理及方法。建立多體接觸振動模型來描述滾動軸承空間運動特性,涵蓋軸承界面復(fù)雜相互作用及摩擦學(xué)系統(tǒng)行為,以接觸表面局部區(qū)域尺寸變化來刻畫損傷發(fā)展趨勢,實現(xiàn)復(fù)雜工況下軸承的損傷建模。通過精細劃分損傷區(qū)域由正問題求解標(biāo)準模式數(shù)據(jù)庫,將反問題計算轉(zhuǎn)化為多維優(yōu)化問題,運用灰色關(guān)聯(lián)度辨識軸承損傷位置和具體尺寸。

1 滾動軸承損傷接觸建模

1.1 接觸剛度系數(shù)

由 Hertzian理論,點接觸負荷-變形關(guān)系為[7]

式中 W為彈性趨近量,Eschmann給出了兩固體點接觸彈性變形表達式

式中 Δ,kb,_是由材料屬性決定的 Hertzian系數(shù);d分別為彈性模量、泊松比與接觸表面曲率和。點接觸的接觸力可表示為

因此軸承接觸剛度系數(shù)為

∑ d的計算方法由 Harris給出,詳見文獻 [7],一旦軸承設(shè)計尺寸和材料確定,K可由式(2)～(4)求出。

1.2 接觸關(guān)系的坐標(biāo)表達

圖1是滾動軸承坐標(biāo)示意圖,第i個鋼球-套圈接觸變形Wi是內(nèi)圈在X,Y方向位移(xs,ys),鋼球位置角θi和游隙c的函數(shù),如下式

設(shè)(xb,yb)為鋼球坐標(biāo),由于振動傳遞作用,考慮到鋼球自身的振動,局部接觸變形變?yōu)?/p>

θi表示第i個鋼球的位置角

式中 kc為軸承公轉(zhuǎn)速度,N為鋼球個數(shù)。設(shè)轉(zhuǎn)速為k,則,其中Db為鋼球直徑,Dp為軸承節(jié)圓直徑。

圖1 滾動軸承示意圖Fig.1 Schematic diagram of rolling bearing

設(shè)(xo,yo)為外圈處的運動坐標(biāo),接觸變形同樣可表示為

1.3 動態(tài)接觸建模

球軸承接觸阻尼隨平均負荷變化,其值很小,Kramer通過頻域測試方法給出了阻尼近似計算公式[8]

故xb方向包含阻尼的動態(tài)接觸力可寫成

1.4 拖動力模型

點線接觸中兩個滾動元件之間的相對滑動速度引起的潤滑劑與滾動元件之間的摩擦力被稱為拖動力。拖動特性是決定軸承元件動態(tài)特性的關(guān)鍵因素之一,是滾動軸承動力學(xué)分析中不可缺少的基本參數(shù)。引入摩擦系數(shù)_,yb和xb方向拖動力的函數(shù)關(guān)系為

摩擦系數(shù)_由彈流潤滑模型決定,設(shè)Λ為油膜厚度與表面粗糙度的比值,則[7]

式中 _bd受潤滑參數(shù)及溫度的影響,在本文使用的潤滑模型中,摩擦系數(shù)_bd是滑 /滾率s的函數(shù)。

通過以上分析,對深溝球軸承的平面運動,考慮到剛度、阻尼、拖動力和滑 /滾率的動態(tài)接觸力表達式為

1.5 運動方程式

滾動軸承的振動存在著振動源和振動傳遞路徑的交互,如圖2所示,系統(tǒng)響應(yīng)最終為傳感器接收。

圖2 滾動軸承系統(tǒng)振動模型Fig.2 Vibration model of rolling bearing system

軸承表面損傷形態(tài)決定了激振力的頻譜,界面動態(tài)接觸決定了振動系統(tǒng)的傳遞特性,最終的振動頻譜由上述二者共同決定。正問題分析中,外圈質(zhì)量和剛度將激勵系統(tǒng)較低的固有頻率?？紤]到軸承的高頻諧振頻帶,設(shè)置mp,cp,kp為附加自由度來模擬軸承-軸承座結(jié)合部的約束以激勵高頻響應(yīng)。將傳遞路徑的質(zhì)量、剛度、阻尼及幾何位置作為系統(tǒng)參數(shù),根據(jù)牛頓定律可建立內(nèi)圈-鋼球-外圈接觸振動微分方程組如下

式中ms,mb,mo分別為內(nèi)圈與軸的質(zhì)量,鋼球質(zhì)量,外圈質(zhì)量;ko,co為外圈剛度與阻尼;fx,fy為徑向載荷與不平衡力的總和。由于工程環(huán)境、安裝條件等因素的影響,使得軸承-軸承座結(jié)合部的材料特性、幾何尺寸和邊界條件具有微小的偏差。為了與實際工況等效,mp,cp,kp的具體數(shù)值可通過健康信號的頻譜加以修正。模型參數(shù)導(dǎo)致的誤差做了誤差置零處理[9],獲得數(shù)值計算與實驗結(jié)果一致的數(shù)值模型。

1.6 表面損傷建模

滾動軸承由疲勞誘發(fā)的表面損傷主要包括裂紋、凹點、剝落等。最早研究軸承損傷振動模型的是McFadden和 Smith,該模型通過軸承運轉(zhuǎn)過程中鋼球經(jīng)過故障位置產(chǎn)生的一系列沖擊與沖擊衰減函數(shù)的卷積來模擬振動[10]。Choudhury A提出外加沖擊力的形式模擬損傷所激勵的動態(tài)效應(yīng),建立了滾動軸承表面損傷的隨機振動數(shù)學(xué)模型[11]。Behzad M等通過基于模型的振動響應(yīng)與傳遞函數(shù)相結(jié)合的方法較好地模擬了軸承的振動行為[12]。

到目前為止,提出的大多數(shù)為滾動軸承損傷的數(shù)字信號仿真與模擬沖擊模型,未能有效地將損傷融入到滾動軸承自身的幾何學(xué)和動力學(xué)方程中。本文通過找出滾動軸承在運轉(zhuǎn)過程中鋼球通過損傷區(qū)的運動變化規(guī)律,在宏觀上建模并與接觸振動相結(jié)合為分析滾動軸承疲勞累積損傷提供一種可用的方法,損傷模型如圖 3所示。

損傷長度為L,深度為h,當(dāng)鋼球運行至損傷區(qū)域時,游隙值將隨之增加h,使得W=W+h,瞬間通過長度L后,游隙值恢復(fù)為W。當(dāng)鋼球處于損傷區(qū)域時,導(dǎo)致接觸面間的 Hertzian接觸力突然降低或變?yōu)榱?相當(dāng)于滾動軸承的非線性時變系統(tǒng)的激勵發(fā)生突變使得振動能量出現(xiàn)波動,損傷程度可用L×h或L/h來定量衡量。

圖3 表面損傷實物與幾何模型Fig.3 Physical and geometric model of surface damage

2 數(shù)值計算與分析

式(15)所代表的狀態(tài)方程組非線性強,難以解析,本文采用數(shù)值方法獲取系統(tǒng)響應(yīng)。在MATLAB中編制了相應(yīng)程序,計算方法為定步長 Runge-Kutta法。算例采用美國凱斯西儲大學(xué)(CWRU)軸承研究中心的深溝球軸承 (6205-2RS SKF),轉(zhuǎn)速1 797 r/min,采樣頻率 12 000 Hz,CWRU提供了軸承詳細設(shè)計參數(shù)[13]。軸承受到來自轉(zhuǎn)子不平衡激勵所產(chǎn)生的強迫振動,其振動頻率為轉(zhuǎn)子的旋轉(zhuǎn)頻率。軸承也將產(chǎn)生由于總剛度連續(xù)周期變化而形成的VC(Varying compliance)振動及其諧波組合[14,15]。這是因為,在徑向載荷的作用下,各鋼球的受力情況是不一樣的,且隨著鋼球上的某一點的運動位置不同受力情況亦不一樣,隨著鋼球相對于徑向載荷作用線的移動,軸承剛度以數(shù)倍于鋼球沿靜止套圈轉(zhuǎn)動的頻率呈周期性變化。這種振動為參數(shù)激勵的通過振動,VC振動頻率與外圈損傷頻率數(shù)值上相等屬于軸承自身本體振動成分,VC=kc×nb=106.7 Hz。圖4的功率譜清晰地顯示了VC振動及其高次諧波,初步表明文中計算具有一定的可靠性。

圖4 正常狀態(tài)波形及功率譜圖Fig.4 Waveform and spectrum of the normal bearing

圖5～7分別是損傷信號波形及功率譜圖,內(nèi)、外圈和鋼球損傷頻率依次為fi=161.3 Hz,fo=106.7 Hz,fb=140.4 Hz。與圖 4比較,損傷時波形具有明顯的周期沖擊特征,這是游隙瞬間變化引起的,波形重復(fù)出現(xiàn)的時間間隔為損傷出現(xiàn)的周期0.006 2,0.009 4,0.007 1 s。

由于接觸振動模型全面考慮了振動的傳遞路徑,誘發(fā)軸承各階固有頻率,使得高頻諧振頻率與損傷頻率調(diào)制,在高頻處譜圖出現(xiàn)以損傷頻率為主的邊頻帶,幅值遠離載波中心程逐漸減小的趨勢。由于內(nèi)圈損傷位置隨軸旋轉(zhuǎn)不斷改變,損傷引起的沖擊響應(yīng)幅值必然受到轉(zhuǎn)頻調(diào)制,會在損傷頻率各次諧波處出現(xiàn)邊頻譜線,造成內(nèi)圈損傷功率譜圖更加復(fù)雜,如圖5所示。試驗信號與仿真結(jié)果取得較好的一致性,也說明了本文建模方法的有效性。

圖5 內(nèi)圈損傷狀態(tài)波形及功率譜圖Fig.5 Waveform and spectrum of the bearing with inner race damage

圖6 外圈損傷狀態(tài)波形及功率譜圖Fig.6 Waveform and spectrum of the bearing with outer racedamage

圖7 鋼球損傷狀態(tài)波形及功率譜圖Fig.7 Waveform and spectrum of the bearing with rolling element damage

3 故障量化特征提取

3.1 故障演變過程描述

將L×h設(shè)為輕度、中度、嚴重 3個階段,分別取值為(0.1,0.1),(0.6,0.8),(1.2,1.5)mm。圖8中(a)～(c)為 3個階段的外圈損傷頻譜,可見譜圖結(jié)構(gòu)差異明顯。輕微損傷時損傷頻率被諧波分量所干擾,低頻調(diào)制邊頻能量較弱,隨著損傷程度加深低頻調(diào)制譜線逐漸清晰,并且譜圖能量分布變化顯著。圖8(d)為損傷發(fā)展時歸一化振動烈度變化曲線,損傷加劇沖擊越明顯,振動烈度值越大;沖擊越不明顯,甚至被VC及其諧波頻率所淹沒時,振動烈度越小。不難看出,通過剖分L×h的數(shù)值可以有效區(qū)分振動響應(yīng)。以L×h的大小量化損傷區(qū)域等效于局部不同柔度的植入,從而誘發(fā)振動響應(yīng)變化,形成模式的多樣性。

圖8 外圈不同損傷程度的頻譜及振動烈度值Fig.8 Spectrum and vibration intensity of the bearing with outer race damagefor different fault severities

3.2 峭度最優(yōu)Laplace小波特征提取

機械診斷屬于模式識別問題,為了增強特征向量的敏感性,達到精確提取軸承損傷時沖擊響應(yīng)信號的目的,本文研究了一種峭度最優(yōu) Laplace小波。R Lind等人構(gòu)造的 Laplace小波是一種復(fù)指數(shù)小波,解析表達式為[16]

式中 參數(shù)矢量γ=[a,kn]決定了小波的性能;成員變量a,kn為模態(tài)動力學(xué)參數(shù);A用來歸一化小波函數(shù)。峭度對沖擊信息非常敏感,又是區(qū)分非高斯分布的指標(biāo),高的峭度值意味著信號含有豐富的沖擊成分?；诓ㄐ纹ヅ涞乃枷?本文用峭度最大化準則定量描述參數(shù)矢量γ。在一定范圍內(nèi)變化a和kn,選擇使得峭度最大的a,kn值為參數(shù),經(jīng)計算獲得最優(yōu)值a=0.85,kn=16,這就是峭度最優(yōu) Laplace小波變換。本質(zhì)上,基于峭度最優(yōu)的 Laplace小波變換也就是控制小波濾波器的通帶帶寬,參數(shù)優(yōu)化相當(dāng)于調(diào)整品質(zhì)因子Q,最終結(jié)果使小波沖擊特征更加顯著,得到與信號最相似的基元函數(shù)。

時域峭度、脈沖等 5個指標(biāo)已被證明可較好區(qū)分多重狀態(tài)[17]。僅依靠時域信息對軸承損傷尺寸的精確識別還不夠理想,因為信息含量太少,頻譜結(jié)構(gòu)變化得不到體現(xiàn)。因此,結(jié)合頻域無量綱指標(biāo)組建標(biāo)準模式,頻域指標(biāo)定義為平均頻率、頻譜偏斜度、頻譜峭度、均方根比[18]。由于工況變動的隨機性和采集信息的噪聲影響,量化診斷具有一定的不確定性。模糊分析是解決識別過程中不確定性問題的有力工具。以 [0,1]之間的模糊隸屬度表示特征向量有助于提升特征敏感度,優(yōu)化特征空間布局,增強信息可靠性。根據(jù)經(jīng)驗可選擇升半嶺型函數(shù)表示隸屬度[17]。至此,構(gòu)造了量化診斷標(biāo)準模式向量

以外圈損傷尺寸 0.001,0.002,0.003 mm為例,根據(jù)模型計算特征向量,結(jié)果如圖9所示。可見,經(jīng)峭度最優(yōu) Laplace小波提取的 9個無量綱指標(biāo)對差值為 0.001 mm的損傷尺寸實現(xiàn)了良好的分類,也說明動力學(xué)響應(yīng)參數(shù)可以較好映射局部損傷尺寸,為量化診斷奠定了基礎(chǔ)。

圖9 外圈不同損傷尺寸的特征向量Fig.9 Feature vectors of the bearing with outer race damage for different fault size

4 局部損傷的定量檢測

4.1 定量識別的灰色關(guān)聯(lián)度分析

量化診斷正問題實際上是通過對含有不同損傷位置和尺寸的滾動軸承進行動力學(xué)求解,獲得損傷特征并轉(zhuǎn)化為模糊向量,由計算得到的離散值構(gòu)建損傷數(shù)據(jù)庫。本文將反問題轉(zhuǎn)化為軸承的實際振動輸出與接觸振動數(shù)值模型的數(shù)值輸出間差異最小化為優(yōu)化目標(biāo)的優(yōu)化問題。利用關(guān)聯(lián)尋優(yōu)求解出與輸入值最相似的樣本點,進而測出損傷位置及損傷長度L,深度h。正問題計算得到標(biāo)準模式數(shù)據(jù)庫

式中i為模式類別,k=9表示特征向量個數(shù)。令待檢模式為?j(k)={xj(1),xj(2),… ,xj(9)},其中j為模式類別,k的意義同上。記Si對Sj的關(guān)聯(lián)系數(shù)為Yij(k),Yij(k)由下式確定[19]

式中 Δij=|Si(k)-Sj(k)|稱為絕對差,d為分辨系數(shù),量化診斷時可取d=0.5,則Si與Sj的關(guān)聯(lián)度為

關(guān)聯(lián)度從物理意義講是對灰色過程的白化處理,從數(shù)學(xué)角度講是空間幾何曲線的比較。反問題辨識損傷尺寸可用如下數(shù)學(xué)關(guān)系式描述

正、反問題進行滾動軸承量化診斷的流程如圖10。

圖10 量化診斷流程Fig.10 Flow chart for the quantitative identification

4.2 案例驗證

案例同樣源自 CWRU軸承中心,采用電火花機床分別在測試軸承的內(nèi)、外圈和鋼球上加工深度為11 mils(0.279 4 mm)的盲孔,盲孔的直徑為 7 mils(0.177 8 mm),14 mils(0.355 6 mm),21 mils(0.533 4 mm),這樣就可以通過實驗數(shù)據(jù)完成量化診斷。測試軸承分別在 4種不同載荷(載荷分別為0～ 3馬力)以及相對應(yīng)的 1 797,1 772,1 750,1 730 r/min四種轉(zhuǎn)速下測試。振動信號由16通道DAT數(shù)據(jù)記錄儀采集,采樣頻率為12 k Hz。 CWRU未提供第4種工況下的外圈數(shù)據(jù),因此本文用內(nèi)、外圈和鋼球在前 3種工況下的樣本各 20個按照圖 10的流程進行量化診斷。

正問題計算區(qū)間定為L=0.1～1 mm,h=0.1～1 mm,搜索步長0.001 mm。表1給出了任意一種工況下h的反問題分析結(jié)果,表 2為 3種工況下L的分析結(jié)果。計算出損傷尺寸預(yù)測值的相對誤差,結(jié)果表明,損傷深度預(yù)測相對誤差不超過 4%,損傷長度的相對誤差不超過 7%。

表1 量化診斷結(jié)果(h/mm)Tab.1 Results for the quantitative identification(h/mm)

表 2 量化診斷結(jié)果(L/mm)Tab.2 Results f or thequantitative identif ication(L/mm)

5 結(jié) 論

(1)滾動軸承的振動響應(yīng)是軸承內(nèi)部幾何關(guān)系、運動關(guān)系在外部載荷作用下的綜合體現(xiàn)。結(jié)合軸承界面復(fù)合約束作用及摩擦學(xué)行為,建立了計及傳遞路徑的滾動軸承動態(tài)接觸振動模型,獲取了軸承不同損傷狀態(tài)下的本體振動信號。該模型可以敏銳地捕捉到軸承損傷時振動信號的沖擊信息且涵蓋了頻譜的中、高頻成分,理論結(jié)果與案例數(shù)據(jù)分析相吻合。

(2)局部損傷的出現(xiàn)引起內(nèi)部游隙的變化,進而影響到軸承的振動特性。不同損傷位置和尺寸形成了不同的損傷模式;反之,運用模式向量也可反求損傷位置和具體尺寸。

(3)峭度對信號中的沖擊成分非常敏感,同時峭度又是一種無量綱指標(biāo),它不受信號絕對水平的影響。根據(jù)特征波形匹配原理,研究峭度最優(yōu) Laplace小波復(fù)合特征提取方法,全面挖掘了振動信號時、頻域動態(tài)信息。

(4)將量化診斷中的反問題轉(zhuǎn)化為優(yōu)化辨識,以實測特征向量作為輸入,運用灰色關(guān)聯(lián)度辨識具體損傷尺寸。診斷實踐表明,識別誤差不超過7%,這為方法的進一步實際應(yīng)用提供了有力的依據(jù)。

[1] Robert B Randall,Jére Antoni.Rolling element bearing diagnostics—— A tutorial[J].Mechanical Systems and Signal Processing,2011,25(2):485—520.

[2] Sawalhi N,Randall R B,Endo H.The enhancement of fault detection and diagnosis in rolling element bearings using minimum entropy deconvolution combined with spectral kurtosis[J].Mechanical Systems and Signal Processing,2007,21(6):2 616—2 633.

[3] Karthik K,Nataraj C.Feature selection for fault detection in rolling element bearings using mutual information[J].Journal of Vibration and Acoustics,2011,133(6):061001—11.

[4] 譚繼勇,陳雪峰,何正嘉.沖擊信號的隨機共振自適應(yīng)檢測方法 [J].機械工程學(xué)報,2010,46(23):61— 67.Tan Jiyong,Chen Xuefeng,He Zhengjia.Impact signal detection method with adaptive stochastic resonance[J].Journal of Mechanical Engineering,2010,46(23):61—67.

[5] 張超,陳建軍,徐亞蘭.基于 EMD分解和奇異值差分譜理論的軸承故障診斷方法[J].振動工程學(xué)報,2011,24(5):539—545.Zhang Chao,Chen Jianjun,Xu Yalan.A bearing fault diagnosis method based on EMD and difference spectrum theory of singular value[J].Journal of Vibration Engineering,2011,24(5):539—545.

[6] 屈梁生.機械監(jiān)測診斷中的理論與方法 [M].西安:西安交通大學(xué)出版社,2009.Qu Liangsheng.Theory and Methodology for Mechanical Condition Monitoring and Fault Diagnosis[M].Xi'an:Xi'an Jiao Tong University Press,2009.

[7] Harris T A.Rolling Bearing Analysis[M].4th ed.New York:John Wiley and Sons,2001.

[8] Kramer E.Dynamic of Rotors and Foundations[M].Berlin:Springer,1993.

[9] 陳雪峰.小波有限元理論與裂紋故障診斷的研究[D].西安:西安交通大學(xué),2004.Chen Xuefeng.Research on the Theory of Wavelet Finite Element Method and Cracks Fault Diagnosis[D].Xi'an:Xi'an Jiao Tong University,2004.

[10]Mc Fadden P D,Smith J D.Model for the vibration produced by a single point defect in a rolling element bearing[J].Journal of Sound and Vibration,1984,96(1),69—82.

[11]Choudhury A, Tandon N.Vibration response of rolling element bearings in a rotor bearing system to a local defect under radial load[J].ASME,Journal of Tribology,2006,128(2):252—261.

[12]Behzad M,Bastami A R,Mba D.A new model for estimating vibrations generated in the defective rolling element bearings[J].ASME,Journal of Vibration and Acoustics,2011,133(4):041011—8.

[13]The Case Western Reserve University Bearing Data Center Website.Bearing data center seeded fault test data[EB/OL].http:http://www.eecs/cwru/edu/laboratory/bearing.

[14]TIWARI M,GUPTA K,PRAKASH O.Effect of radial internal clearance of a ball bearing on the dynamics of a balanced horizontal rotor[J].Journal of Sound and Vibration,2000,238(5):723—756.

[15]TIWARIM,GUPTA K,PRAKASH O.Dynamic response of an unbalanced rotor supported on ball bearings[J].Journal of Sound and Vibration,2000,238(5):757—779.

[16]何正嘉,訾艷陽,張西寧.現(xiàn)代信號處理及工程應(yīng)用[M].西安:西安交通大學(xué)出版社 ,2007.He Zhengjia,Zi Yanyang,Zhang Xining.Modern Signal Processing and Engineering Application[M].Xi'an:Xi'an Jiao Tong University Press,2007.

[17]徐敏,虞和濟,張瑞林.設(shè)備故障診斷手冊 [M].西安:西安交通大學(xué)出版社,1998.Xu Min,Yu Heji,Zhang Ruilin.Handbook of Machinery Fault Diagnosis[M].Xi'an:Xi'an Jiao Tong University Press,1998.

[18]Chen P,Masatoshi T,Toshio T.Fault diagnosis method for machinery in unsteady operating condition by instantaneous power spectrum and genetic programming[J].Mechanical Systems and Signal Processing,2005,19(1):175—194.

[19]鄧聚龍.灰色系統(tǒng)基本方法 [M].武漢:華中科技大學(xué)出版社,2005.Deng Julong.The Primary Methods of Grey System Theory[M].Wuhan: Huazhong University of Science& Technology Press,2005.