姚 軍 王 歡 蘇 泉

（北京航空航天大學(xué) 可靠性與系統(tǒng)工程學(xué)院，北京 100191）

為盡快預(yù)測(cè)高可靠、長(zhǎng)壽命產(chǎn)品的壽命，通常采用加速試驗(yàn).它通過(guò)提高應(yīng)力的方式加速產(chǎn)品的失效，從而根據(jù)高應(yīng)力下產(chǎn)品的壽命來(lái)外推正常應(yīng)力下產(chǎn)品的壽命.為了保證這種統(tǒng)計(jì)推斷的準(zhǔn)確性，產(chǎn)品在正常應(yīng)力和加速應(yīng)力下必須具有相同的失效機(jī)理，這是進(jìn)行加速壽命試驗(yàn)設(shè)計(jì)的重要前提［1］.

目前，關(guān)于加速試驗(yàn)失效機(jī)理一致性檢驗(yàn)方法的研究，主要分為三類［2］.第一類是基于加速模型參數(shù)不變的一致性檢驗(yàn)方法.其原理是在激活能不變的前提下，可利用產(chǎn)品壽命與溫度應(yīng)力之間的對(duì)數(shù)線性關(guān)系來(lái)檢驗(yàn)阿倫尼斯、Coffin-Man-son、逆冪率、艾林模型的失效機(jī)理一致性.郭春生等人對(duì)二極管3DG130進(jìn)行150～310℃的序進(jìn)應(yīng)力加速試驗(yàn)時(shí)發(fā)現(xiàn)在失效機(jī)理不發(fā)生改變時(shí)，電子器件失效敏感參數(shù)的退化速率與施加應(yīng)力的負(fù)倒數(shù)遵從指數(shù)關(guān)系，從而確定了失效機(jī)理一致的應(yīng)力范圍［3］.該方法缺點(diǎn)在于：第一，失效機(jī)理不變的情況下，激活能是否隨溫度變化，學(xué)術(shù)界尚未有一個(gè)統(tǒng)一的認(rèn)識(shí)；第二，該方法只適用于某幾類包含激活能的加速模型，適用范圍較局限.

第二類是基于試驗(yàn)觀察的失效機(jī)理一致性檢驗(yàn)方法，主要包括直接觀察法、化學(xué)分析法、顯微鏡觀察法3種.例如，包亦望等人在研究鈦化物陶瓷的高溫蠕變行為與失效機(jī)理時(shí)，對(duì)鈦化物樣品的常溫和高溫800，1 000，1 100，1 200℃斷口進(jìn)行了掃描電鏡微觀分析比較，來(lái)判別鈦化物陶瓷的失效機(jī)理一致性［4］.與加速試驗(yàn)失效機(jī)理一致性檢驗(yàn)方法相比，試驗(yàn)觀察的方法適用范圍更廣泛，對(duì)失效機(jī)理變化的判別更為直接，這種基于試驗(yàn)觀察判別方法也有缺點(diǎn)，主要體現(xiàn)在以下幾個(gè)方面.①有些產(chǎn)品的失效機(jī)理變化不容易被觀察到，化學(xué)分析和顯微觀察的方法有可能由于技術(shù)水平不足或放大倍數(shù)不夠等原因而沒(méi)有發(fā)現(xiàn)失效機(jī)理發(fā)生變化的真正應(yīng)力水平.②該方法在使用上有一定的局限性，化學(xué)分析的方法只適用于加速腐蝕試驗(yàn)；顯微觀察的方法不適用于復(fù)雜產(chǎn)品，因?yàn)閺?fù)雜產(chǎn)品失效機(jī)理的變化也很復(fù)雜，如果將產(chǎn)品每個(gè)部分都進(jìn)行顯微分析，工作量會(huì)非常大，成本也會(huì)非常高，甚至是不可行的，且難以找出由產(chǎn)品各部分的聯(lián)系引起的失效機(jī)理變化；直接觀察的方法雖然可以用于復(fù)雜產(chǎn)品，但某些產(chǎn)品的失效模式和失效機(jī)理不一定一一對(duì)應(yīng)，直接觀察法觀察到的失效模式變化并不一定代表失效機(jī)理發(fā)生了變化，沒(méi)有觀察到失效模式的變化也不能說(shuō)明失效機(jī)理就沒(méi)有變化，這使得直接觀察法缺乏說(shuō)服力.

第三類是基于統(tǒng)計(jì)方法的一致性檢驗(yàn)方法，如針對(duì)威布爾分布的產(chǎn)品形狀參數(shù)一致性檢驗(yàn)等.這種方法的原理是若產(chǎn)品的失效機(jī)理不變，則加速系數(shù)是與可靠度值無(wú)關(guān)的常數(shù)，由加速系數(shù)與可靠度值無(wú)關(guān)可以推出產(chǎn)品壽命分布的變異系數(shù)不變.這樣，只要對(duì)變異系數(shù)進(jìn)行假設(shè)檢驗(yàn)，若變異系數(shù)發(fā)生變化，則可推出失效機(jī)理發(fā)生了變化.文獻(xiàn)［5］給出了常見(jiàn)19種壽命分布的失效機(jī)理不變條件，其中關(guān)于Weibull分布、對(duì)數(shù)正態(tài)分布、Gamma分布的失效機(jī)理不變條件與文獻(xiàn)［6］根據(jù)試驗(yàn)結(jié)果與工程經(jīng)驗(yàn)得出的結(jié)果一致.與前兩種方法相比，該方法適用范圍更廣泛.這種方法也有缺點(diǎn)，主要表現(xiàn)為：①統(tǒng)計(jì)方法的假設(shè)物理意義不明確，沒(méi)有從物理的角度給出變異系數(shù)與機(jī)理的關(guān)系；②目前統(tǒng)計(jì)方法只能用于加速試驗(yàn)數(shù)據(jù)，只能用于事后檢驗(yàn)，不能在加速試驗(yàn)之前對(duì)其進(jìn)行理論指導(dǎo).

材料是任何產(chǎn)品的物質(zhì)基礎(chǔ)，產(chǎn)品的失效在廣義上一般都可歸結(jié)為材料失效.失效機(jī)理的變化是材料本身屬性的變化，在某種確定的外因（應(yīng)力、溫度、介質(zhì)等）條件下，只呈現(xiàn)出某種確定的材料屬性.因此，若外因到達(dá)失效機(jī)理變化點(diǎn)，材料內(nèi)部會(huì)發(fā)生某種物理化學(xué)或結(jié)構(gòu)上的變化（如合金的同素異構(gòu)過(guò)程、樹(shù)脂的玻璃化過(guò)程等），致使材料屬性發(fā)生變化，引起產(chǎn)品的失效機(jī)理變化，表征到宏觀結(jié)構(gòu)就是產(chǎn)品被測(cè)性能參數(shù)的突變［7－8］.灰色理論是一種處理少數(shù)據(jù)不確定性問(wèn)題的理論，GM（1，1）建模方法是灰色理論的一個(gè)重要分支，可解決少數(shù)據(jù)壽命預(yù)測(cè)問(wèn)題［9－10］.等維新息模型是一種在選擇建模數(shù)據(jù)時(shí)以最新的實(shí)際數(shù)據(jù)作為參考點(diǎn)，去掉老數(shù)據(jù)，保存數(shù)列等維，依次建模的一種模型.本文采用GM（1，1）建模方法與等維新息模型結(jié)合的一種預(yù)測(cè)方法，對(duì)預(yù)試驗(yàn)數(shù)據(jù)進(jìn)行建模，得出各應(yīng)力點(diǎn)性能預(yù)測(cè)值后與實(shí)際值進(jìn)行殘差計(jì)算，觀察殘差點(diǎn)的突變情況，以此來(lái)找尋失效機(jī)理的變化點(diǎn).

1 理論分析

1.1 基于GM（1，1）方法的灰色預(yù)測(cè)模型

設(shè)預(yù)試驗(yàn)中得到的性能退化序列為

其中，ξ（0）（1）為最小應(yīng)力對(duì)應(yīng)的性能數(shù)據(jù)；ξ（0）（n）為最大應(yīng)力對(duì)應(yīng)的性能數(shù)據(jù).

為將不同極性的性能退化序列轉(zhuǎn)化為同極性（最大值極性）的序列，將ξ（0）按式（1）進(jìn)行一階累加，灰生成新序列［10］

將性能退化序列累加后，弱化了其隨機(jī)性，由于建立的模型是1階1變量序列，可采用GM（1，1）模型建模方法，則可認(rèn)為ξ（1）滿足一階線性微分方程［4］.

式中，a為發(fā)展灰數(shù)，反映了性能退化序列ξ（0）和累加性能退化序列ξ（1）的發(fā)展趨勢(shì)；u稱為內(nèi)生控制灰數(shù)，反映了數(shù)據(jù)間的變化關(guān)系.

為求解a與u，令＝［a，u］T為待估向量，由于性能退化序列是離散的，將式（2）中離散化，則有

其中，Z（1）（k）為式（2）的背景值；μ為權(quán)重系數(shù)，μ∈［0，1］.

根據(jù)灰色預(yù)測(cè)方法，一般μ取值為0.5，則

將式（2）離散化后，可得到累加生成性能退化序列的離散化公式

利用最小二乘法對(duì)式（6）進(jìn)行求解，可得出待估參數(shù)為

根據(jù)得出的待估參數(shù)，代入式（2），可求解出累加生成性能退化序列ξ（1）的預(yù)測(cè)值為

對(duì)式（8）累減還原可得性能退化序列的預(yù)測(cè)模型為

1.2 基于等維新息模型的預(yù)測(cè)方法

雖然，式（9）可以作為長(zhǎng)期預(yù)測(cè)模型，但真正具有實(shí)際意義、精度較高的僅僅是ξ（0）（n）以后的一兩個(gè)數(shù)據(jù)，其他更遠(yuǎn)的數(shù)據(jù)不是預(yù)測(cè)數(shù)據(jù)而是規(guī)劃性的數(shù)據(jù).可以認(rèn)為越往未來(lái)發(fā)展，模型計(jì)算的預(yù)測(cè)數(shù)據(jù)，其預(yù)測(cè)意義就越小.隨著系統(tǒng)的發(fā)展，老數(shù)據(jù)的信息意義將逐步降低，在不斷補(bǔ)充新信息的同時(shí)，及時(shí)地去掉老信息，建模序列更能反映系統(tǒng)在目前的特征.

基于GM（1，1）方法的建模數(shù)據(jù)允許做不同的取舍，但必須等距、相鄰、不得有跳躍.在選擇建模數(shù)據(jù)時(shí)以最新的實(shí)際數(shù)據(jù)作為參考點(diǎn)，去掉最老的數(shù)據(jù)，保持?jǐn)?shù)列等維，依次建模，這樣建立的模型就是等維新息模型.

在預(yù)測(cè)前，必須確定模型的維度，為使模型達(dá)到最高的精度，可選擇維度分別為4～11維（至少為4）.

1.3 基于灰色理論的失效機(jī)理一致性檢驗(yàn)方法

1）選擇適當(dāng)?shù)木S度，對(duì)強(qiáng)化試驗(yàn)數(shù)據(jù)建立等維新息模型；

2）根據(jù)數(shù)據(jù)的特點(diǎn)選擇適當(dāng)?shù)幕疑Ｐ蛯?duì)等維新息模型中的每一段數(shù)據(jù)進(jìn)行預(yù)測(cè)，得到一個(gè)預(yù)測(cè)值；

3）將預(yù)測(cè)值與對(duì)應(yīng)的真實(shí)值進(jìn)行比較，計(jì)算殘差；

4）觀察殘差是否發(fā)生顯著的變化，若有，找出第一個(gè)顯著變化點(diǎn)，則認(rèn)為該點(diǎn)對(duì)應(yīng)的地方機(jī)理發(fā)生了變化，若沒(méi)有顯著變化，則認(rèn)為產(chǎn)品的機(jī)理在強(qiáng)化試驗(yàn)中沒(méi)有發(fā)生變化.

2 分析實(shí)例

2.1 試驗(yàn)樣品

本次預(yù)試驗(yàn)圍繞某型號(hào)衛(wèi)星雷達(dá)天線驅(qū)動(dòng)機(jī)構(gòu)上用光電編碼器展開(kāi)，影響編碼器可靠性的環(huán)節(jié)主要是發(fā)光二極管、光敏三極管在空間熱環(huán)境及輻射環(huán)境下性能漂移.本次試驗(yàn)樣品為配對(duì)使用的紅外固體發(fā)光二極管（型號(hào)：OP224TXV）及光敏三極管 （型號(hào)：OP604TXV），圖 1 為OP224TXV及OP604TXV的外觀圖，圖2為兩器件配對(duì)使用連接示意圖.考察的性能參數(shù)為接收元件光敏三極管信號(hào)幅值.

圖1 OP224TXV（左）及OP604TXV（右）外觀圖

2 OP224TXV（左）及OP604TXV（右）配對(duì)使用連接圖

2.2 試驗(yàn)步驟

試驗(yàn)開(kāi)始前，常溫下測(cè)量樣品的性能參數(shù)，作為試驗(yàn)中比較的基準(zhǔn).試驗(yàn)中每個(gè)溫度臺(tái)階駐留時(shí)間為1h，前10min（t1）為保溫階段，后50min（t2）初始和末端各測(cè)量一次性能參數(shù).初始環(huán)境溫度為95℃，初始步進(jìn)量級(jí)為10℃，每次試驗(yàn)將試驗(yàn)中測(cè)試數(shù)據(jù)與初始測(cè)試數(shù)據(jù)進(jìn)行比較，觀察其性能退化趨勢(shì)，試驗(yàn)剖面如圖3所示.試驗(yàn)設(shè)備如圖4所示，測(cè)試設(shè)備如圖5所示.

圖3 工作電應(yīng)力下溫度極限試驗(yàn)剖面

圖4 試驗(yàn)設(shè)備

圖5 測(cè)試設(shè)備

2.3 試驗(yàn)數(shù)據(jù)分析及處理

預(yù)試驗(yàn)數(shù)據(jù)如表1所示.

表1 光電編碼器強(qiáng)化試驗(yàn)數(shù)據(jù)

試驗(yàn)數(shù)據(jù)共10個(gè)，選擇等維新息模型的維度為4，對(duì)各數(shù)據(jù)段進(jìn)行建模，計(jì)算出發(fā)展系數(shù)a，見(jiàn)表2.

表2 各數(shù)據(jù)段灰色預(yù)測(cè)的發(fā)展系數(shù)

用各數(shù)據(jù)段建立的預(yù)測(cè)模型推測(cè)后一個(gè)性能參數(shù)值，例如，用1～4數(shù)據(jù)建立的預(yù)測(cè)模型來(lái)推測(cè)第5個(gè)性能參數(shù)值，并與實(shí)際觀測(cè)值進(jìn)行對(duì)比，計(jì)算殘差，如表3.用表3數(shù)據(jù)，作預(yù)測(cè)值及觀測(cè)值走勢(shì)圖，如圖6，為檢驗(yàn)殘差突變點(diǎn)，作殘差圖，見(jiàn)圖7.從圖6可以看出，預(yù)測(cè)值及觀測(cè)值的整體趨勢(shì)呈下降趨勢(shì)，在165℃之前，預(yù)測(cè)值及觀測(cè)值下降幅度很小，在各應(yīng)力觀測(cè)點(diǎn)的預(yù)測(cè)值及觀測(cè)值大小交替；從圖7可以看出，135，145，155和165℃的殘差值在0的上下浮動(dòng)，在175℃附近，殘差出現(xiàn)突變點(diǎn)，可以判定此處性能參數(shù)發(fā)生了突變，失效機(jī)理發(fā)生了改變.

表3 信號(hào)幅值的預(yù)測(cè)值、真實(shí)值及殘差

圖6 各溫度應(yīng)力下觀測(cè)值及預(yù)測(cè)值走勢(shì)

圖7 各溫度應(yīng)力下殘差

2.4 基于統(tǒng)計(jì)方法的失效機(jī)理一致性檢驗(yàn)

為驗(yàn)證試驗(yàn)件在175℃附近發(fā)生失效機(jī)理變化，采用加速壽命試驗(yàn)數(shù)據(jù)進(jìn)行基于統(tǒng)計(jì)方法的失效機(jī)理一致性檢驗(yàn)分析，加速壽命試驗(yàn)數(shù)據(jù)如表4，產(chǎn)品失效判據(jù)為光敏三極管的信號(hào)幅值降至初值的60%以下.

表4 光電編碼器加速壽命試驗(yàn)數(shù)據(jù)

假定光電編碼器在各溫度應(yīng)力臺(tái)階下其失效時(shí)間服從威布爾分布，威布爾的變異系數(shù)為形狀參數(shù)m，即從形狀參數(shù)的變化可以表征失效機(jī)理的變化.計(jì)算各溫度應(yīng)力下的威布爾形狀參數(shù)m，可以得出：160℃時(shí)，m＝29.3；175℃ 時(shí)，m＝16.017 2；190℃時(shí)，m＝13.719 2；220℃時(shí)，m＝13.16.很明顯，在175，190及220℃時(shí)，m值基本一致，其失效機(jī)理一致，而160℃時(shí)，m值明顯大于其他應(yīng)力水平，其失效機(jī)理與其他應(yīng)力情況下的不一致，也就是說(shuō)，在175℃附近，失效機(jī)理發(fā)生了改變.與基于灰色理論失效機(jī)理一致性檢驗(yàn)方法做出的判斷一致.

2.5 失效機(jī)理分析

影響信號(hào)幅值變化最根本的原因在于發(fā)光二極管的光功率的衰退.對(duì)于發(fā)光二極管來(lái)說(shuō)，造成光功率的衰退可能有兩方面的原因：LED結(jié)溫的產(chǎn)生和LED熱阻的生成.

LED結(jié)溫是指P-N結(jié)區(qū)的溫度，其形成是由于LED空穴、電子運(yùn)動(dòng)，一部分能量產(chǎn)生有效的光電效應(yīng)，發(fā)出光子；另一部分是以發(fā)熱的形式消耗掉了，從而導(dǎo)致P-N結(jié)區(qū)芯片發(fā)熱.對(duì)于一個(gè)封裝好的LED發(fā)光管來(lái)說(shuō)，產(chǎn)生結(jié)溫最重要的因素在于LED封裝的散熱能力.LED熱阻是LED散熱通道上的兩個(gè)節(jié)點(diǎn)之間熱功率流的阻值，其封裝的材質(zhì)及形式對(duì)LED熱阻的產(chǎn)生有很大的影響.

試驗(yàn)用LED發(fā)光管其封裝環(huán)氧膠在高溫下會(huì)發(fā)生變性、發(fā)黃，出光效率下降衰減.且LED用的封裝環(huán)氧膠存在著一個(gè)重要特性，即當(dāng)環(huán)氧膠溫度超過(guò)一個(gè)特定溫度時(shí)，封裝環(huán)氧的特性將從一種剛性的類玻璃狀態(tài)轉(zhuǎn)變成一種柔軟的似橡膠態(tài)狀的物質(zhì).此時(shí)材料的膨脹系數(shù)急劇增加，形成一個(gè)明顯的拐點(diǎn)，這個(gè)拐點(diǎn)所對(duì)應(yīng)的溫變即為環(huán)氧樹(shù)脂的玻璃狀轉(zhuǎn)化溫度，當(dāng)器件在此溫度附近或高于此溫度變化時(shí)，將發(fā)生明顯的膨脹或收縮，致使二極管發(fā)光功率出現(xiàn)明顯衰退.

對(duì)試驗(yàn)用二極管環(huán)氧膠為雙酚A型環(huán)氧樹(shù)脂膠，圖8為雙酚A型環(huán)氧樹(shù)脂膠在140℃（左上）、160℃（右上）、170℃（左下）、180℃（右下）溫度下 斷 裂 面 （SEM，Scanning Electron Microscope）照片.從圖8可以看出，在140℃及160℃，膠粘劑成剛性狀態(tài)，而170℃時(shí)，膠粘劑成絮狀物的橡膠形態(tài)，180℃則更加明顯.由此可知，170℃左右為封裝用環(huán)氧膠材料屬性變化溫度點(diǎn)，該材料屬性的變化致使產(chǎn)品失效機(jī)理發(fā)生改變，所以從宏觀上能觀察到在175℃處，接受元件光敏三極管信號(hào)幅值殘差出現(xiàn)突變點(diǎn).

圖8 膠粘劑斷裂面SEM圖

3 結(jié)束語(yǔ)

本文提出了一種基于灰色理論的失效機(jī)理一致性判別方法，它克服了傳統(tǒng)方法的缺陷，與基于統(tǒng)計(jì)的方法相比，此方法僅需要預(yù)試驗(yàn)數(shù)據(jù)就可對(duì)失效機(jī)理一致性進(jìn)行檢驗(yàn)，可指導(dǎo)加速壽命試驗(yàn).與基于試驗(yàn)觀察的方法對(duì)比，此方法僅需要少量數(shù)據(jù)就可做出準(zhǔn)確的檢驗(yàn)，簡(jiǎn)單明了.與基于加速模型參數(shù)不變的方法相比，本方法從本質(zhì)上來(lái)說(shuō)是一種統(tǒng)計(jì)方法，與失效機(jī)理、激活能等物理化學(xué)理論無(wú)關(guān)，應(yīng)用范圍更為廣泛.

（References）

［1］茆詩(shī)松，王玲玲.可靠性統(tǒng)計(jì)［M］.上海：華東師范大學(xué)出版社，1984

Mao Shisong，Wang Lingling.Reliability statistics［M］.Shang hai：East China Normal University Press，1984（in Chinese）

［2］Pan Xiaoxi，Huang Xiaokai，Chen Yunxia，et al.Connotation of failure mechanism consistency and identification method for accelerated testing［C］／／Prognostics and System Health Management Conference.［S.l.］：IEEE，2011：1－7

［3］郭春生，謝雪松，馬衛(wèi)東，等.加速試驗(yàn)中失效機(jī)理一致性的判別方法［J］.半導(dǎo)體學(xué)報(bào)，2006，27（3）：560－563

Guo Chunsheng，Xie Xuesong，Ma Weidong，et al.Failuremechanism identification method in accelerated testing［J］.Journal of Semiconductors，2006，27（3）：560－563（in Chinese）

［4］包亦望，蘇盛彪，王毅敏，等.鈦化物陶瓷的高溫蠕變行為與失效機(jī)理［J］.硅酸鹽學(xué)報(bào)，2002，30（3）：300－304

Bao Yiwang，Su Shengbiao，Wang Yimin，et al.High-temperature creep and failure mechanism of ceramics containing titaniferous compound［J］.Journal of the Chinese Ceramic Society，2002，30（3）：300－304（in Chinese）

［5］周源泉，翁朝曦，葉喜濤.論加速系數(shù)與失效機(jī)理的不變的條件（I）-壽命型隨機(jī)變的情況［J］.系統(tǒng)工程與電子技術(shù)，1996（1）：55－67

Zhou Yuanquan，Weng Chaoxi，Ye Xitao.Study on accelerated factor and conditionfor constant failure mechanism（Ⅰ）—the case for lifetime is a random variable［J］.Journal of Systems Engineering and Electronics，1996（1）：55－67（in Chinese）

［6］Nelson W.Analysis of accelerated life test data-part I：the arrhenius model and graphical methods［J］.Transactions on E-lectrical Insulation，1971，6（4）：165－181

［7］鹽見(jiàn)弘.失效物理基礎(chǔ)［M］.楊家鏗譯.北京：科學(xué)出版社，1982

Shiomi Hiroshi.Basis of failure physics［M］.Translated by Yang Jiakeng.Beijing：Science Press，1982（in Chinese）

［8］Hu J M，Barker D，Dasgupta A，et al.Role of failure-mechanism identification in accelerated testing［C］／／Proceedings Annual Reliability and Maintainability Symposium.Las Vegas.NV：IEEE，1992：181－188［9］Atwa Y M，El-Saadany E F.Annual wind speed estimation utilizing constrained grey predictor［J］.Transaction on Energy Conversion IEEE，2009，24（2）：548－550

［10］鄧聚龍.灰色系統(tǒng)基本方法［M］.武漢：華中科技大學(xué)出版社，2005

Deng Julong.Basic methods of gray system［M］.Wuhan：Huazhong University of Science and Technology Press，2005（in Chinese）

［11］Chang B R，Tsai H F.Forecast approach using neural network adaptation to support vector regression grey model and generalized auto-regressive conditional heteroscedasticity［J］.Expert Systems with Applications，2008，34（2）：925－934

［12］Huang K Y，Jane C J.A hybrid model for stock market forecasting and portfolio selection based on ARX grey system and RS theories［J］.Expert Systems with Applications，2009，36（3）：5387－5392

［13］Truong D Q，Ahn K K.Force control for hydraulic load simulator using self-tuning grey predictor-fuzzy PID ［J］.Mechatronics，2009，19（2）：233－246

［14］Kayacan E，Ulutas B，Kaynak O.Grey system theory-based models in time series prediction［J］.Expert Systems with Applications，2010，37（2）：1784－1789