本發(fā)明涉及一種應用于數(shù)控機床載荷譜領域的技術方法，更確切地說，本發(fā)明涉及一種在對數(shù)控機床時域載荷擴展的基礎上，進行時域載荷pot外推的方法。

背景技術：

載荷譜是產(chǎn)品整機或零部件所承受的載荷時間歷程，由于能夠反映隨機載荷的變化規(guī)律，因此載荷譜被普遍應用到零部件結構的可靠性試驗、疲勞壽命預測、可靠性設計、加速試驗等方面。

數(shù)控機床在不同的切削工況和工作模式下具有不同的切削載荷，其載荷相比與其他機械產(chǎn)品要復雜多變，另外數(shù)控機床自身不具備載荷測試的功能，如果安裝相應的測試裝置，又會給現(xiàn)場切削加工帶來很多不便，所以只能在實驗室進行載荷測取試驗。然而在進行實際載荷測試過程中，受測試時間和成本等因素的影響，只能對產(chǎn)品全壽命周期內的一小部分時間進行載荷測試。為了能夠獲取全壽命周期內可能出現(xiàn)的載荷，必須對實測的載荷進行外推。

時域載荷外推是一種新興的載荷外推方法，它是針對實測的載荷時域信號，采用一維分布對載荷樣本數(shù)據(jù)進行擬合，利用建立的概率密度函數(shù)從而實現(xiàn)了載荷的外推。其中最常見的極值抽樣模型為超閾模型(peakoverthreshold，簡稱pot模型):在一定條件下確定一個閾值u，然后提取所有大于該閾值的極值點構成新的極值樣本，再進行極值分布的擬合。在進行載荷外推之前，需要根據(jù)外推指定的載荷作用時間或者作用頻次對時域載荷進行擴展，而傳統(tǒng)的時域擴展方法未考慮極值載荷間時間間隔的差異，時間間隔都采用采樣頻率f的倒數(shù)，其與實際情況存在很大的差異。同時時域外推中極值閾值的選取至關重要，通常將分析閾值的大小等效為極值樣本廣義帕累拖分布(generalizedparetodistribution，簡稱gpd)擬合的好壞，然而針對gpd沒有公認的最優(yōu)擬合優(yōu)度檢驗，所以對閾值的選擇帶來了一定的困難。因此，本發(fā)明針對傳統(tǒng)的時域載荷擴展方法未考慮極值載荷時間間隔的差異及選取gpd擬合的最優(yōu)檢驗量等問題，提出了一種基于時域載荷擴展的數(shù)控機床載荷外推方法。

技術實現(xiàn)要素：

本發(fā)明所要解決的技術問題是克服了目前傳統(tǒng)的時域載荷擴展方法未考慮極值載荷時間間隔的差異及選取gpd擬合的最優(yōu)檢驗量等問題，提出了一種基于時域載荷擴展的數(shù)控機床載荷外推方法。

為解決上述技術問題，本發(fā)明是采用如下技術方案實現(xiàn)的：

一種基于時域載荷擴展的數(shù)控機床載荷外推方法，包括以下步驟：

步驟一：基于mcmc的時域載荷擴展法；

步驟二：利用mef圖確定閾值范圍；

步驟三：基于灰關聯(lián)分析的最優(yōu)閾值的選??；

步驟四：獲取外推的時域載荷。

步驟一中所述基于mcmc的時域載荷擴展法包括以下步驟：

step1提取時域載荷的極值及升程；

step2劃分載荷等級水平；

step3統(tǒng)計載荷循環(huán)次數(shù)及升程；

step4求解mcmc的狀態(tài)轉移矩陣；

step5求解三個狀態(tài)轉移概率矩陣的累積概率矩陣pcum，pc'um，p”cum；

step6計算馬爾科夫鏈的平穩(wěn)分布π(xt)；

step7在(0,1)區(qū)間內隨機產(chǎn)生符合均勻分布的隨機數(shù)，進行載荷狀態(tài)的判斷，確定載荷水平等級；

step8生成載荷時間歷程的轉折點序列；

step9確定相鄰極值點之間的間隔時間。

步驟二中所述利用mef圖確定閾值范圍包括以下步驟：

step1繪制載荷mef圖；

step2確定閾值范圍[umin,umax]及劃分成m個閾值。

步驟三中所述基于灰關聯(lián)分析的最優(yōu)閾值的選取具體步驟如下：

step1對不同閾值u對應的超越量進行gpd擬合；

step2擬合優(yōu)度檢驗；

step3計算采樣分布和擬合分布序列的d1,d2值；

step4求增量序列；

step5計算時間段的關聯(lián)系數(shù)；

step6計算灰關聯(lián)度。

步驟四中所述獲取外推的時域載荷是指根據(jù)超越量的數(shù)量，模擬產(chǎn)生服從最優(yōu)gpd的新超越量，并用新產(chǎn)生的超越量來代替原有的超越量，從而獲取外推的載荷時域信號。

技術方案中所述提取時域載荷的極值及升程，具體步驟如下：

(1)對實測的載荷時間歷程采用變程閾值模型進行轉折點的取舍，剔除無效的載荷循環(huán)；

(2)提取數(shù)控機床實測載荷時間歷程的轉折點{xt}作為極值載荷時間序列；

(3)提取時域載荷信號中的轉折點中的極大值點m1和極小值點m1；

式中，m1，m2，…,mn代表局部極小值點，m1，m2，…,mn代表局部極大值點，n表示極值點的個數(shù)。

(4)從轉折點{xt}中提取相鄰極值點的差值及對應相鄰極值點的間隔時間，將第t個相鄰極值點的差值記為升程st，對應的相鄰極值點的間隔時間記為tt，則第t個升程和對應的間隔時間記為tst＝(tt,st)。

技術方案中所述劃分載荷等級水平是指：

將提取的載荷極大值和極小值通過固定的載荷等級水平劃分為nd組；

將升程st通過固定的載荷等級水平劃分為nd組；

間隔時間為離散值，是頻率的倒數(shù)的整數(shù)倍。

技術方案中所述統(tǒng)計載荷循環(huán)次數(shù)及升程是指：

統(tǒng)計轉折點{xt}中從時域載荷的極大值到極小值和極小值到極大值狀態(tài)轉移的循環(huán)次數(shù)，將循環(huán)矩陣分別記為f和f’；

統(tǒng)計求解出每個狀態(tài)的升程對應不同間隔時間的頻次，將升程st及對應的間隔時間tt的頻次矩陣記為f”。

技術方案中所述求解mcmc的狀態(tài)轉移矩陣包括以下步驟：

(1)求解載荷極大值到極小值和極小值到極大值的狀態(tài)轉移矩陣，將循環(huán)矩陣f(i,j)和f’(j,i)的各級頻數(shù)除以各級頻數(shù)的總和，得到狀態(tài)轉移概率矩陣p和p’；

(2)將升程及間隔時間的頻次矩陣f”(j,i)的各級頻數(shù)除以總的頻數(shù)得到狀態(tài)轉移概率矩陣p”；

技術方案中所述三個狀態(tài)轉移概率矩陣的累積概率矩陣pcum，p'cum，p”cum為：

技術方案中所述計算馬爾科夫鏈的平穩(wěn)分布π(xt)

技術方案中所述進行載荷狀態(tài)的判斷，確定載荷水平等級，具體步驟如下：

(1)將產(chǎn)生的隨機數(shù)跟累計概率矩陣pcum中的第i行元素分別進行對比，假設其位于第j-1和j之間的狀態(tài)，則可以確定極大值m1所處的載荷等級n1，

(2)在(0,1)區(qū)間產(chǎn)生隨機數(shù)，將隨機數(shù)與累計概率矩陣p’cum中的第j行元素分別進行對比，確定極小值m2所處的載荷等級n2；

(3)再生產(chǎn)隨機數(shù)，確定下一個極大值m2所處的載荷等級n3，以此類推，即可以獲得載荷的狀態(tài)等級序列n。

技術方案中所述生成載荷時間歷程的轉折點序列是指：

依據(jù)載荷狀態(tài)等級序列生成載荷時間歷程的轉折點序列，原始載荷序列xt，載荷劃分的等級數(shù)nd，載荷的狀態(tài)等級序列n，根據(jù)式1.9可以生成載荷時間歷程的轉折點序列；

技術方案中所述確定相鄰極值點之間的間隔時間，步驟如下：

(1)根據(jù)生成的載荷時間歷程的轉折點序列，計算出相鄰極值點的差值，即升程st；

(2)在(0,1)區(qū)間內隨機產(chǎn)生符合均勻分布的隨機數(shù)，根據(jù)升程st，將隨機數(shù)與p”cum的第i行元素進行比較，如果這個隨機數(shù)位于狀態(tài)j-1和j之間，則可以確定該升程st對應的間隔時間tt，以此類推，確定出所有相鄰極值點對應的間隔時間，完成時域載荷的擴展。

技術方案中所述繪制載荷mef圖，具體步驟如下：

(1)定義e(u)＝e(x–u|x>u)為載荷極值x的超越量均值函數(shù)mef，載荷極值樣本經(jīng)驗超越量均值函數(shù)為：

(2)繪制由超越量均值函數(shù)值e(u)與極值閾值u的散點圖構成的曲線，如果在大于某一閾值ud之后，mef的曲線接近于斜率為正的直線，則載荷極值樣本服從形狀參數(shù)ξ＞0的gpd分布；如果直線段的斜率為負，表示該載荷極值樣本為薄尾分布，服從尾部較短的gpd分布；如果直線段為水平線，則表明載荷極值樣本遵從指數(shù)分布。

技術方案中所述確定閾值范圍[umin,umax]及劃分成m個閾值，具體步驟如下：

(1)根據(jù)繪制的mef圖，確定mef曲線接近為直線部分的最左端對應的閾值ud；

(2)如果閾值為離散的值，則以閾值ud為中心，從閾值ud相鄰的閾值中選出m個閾值，構成閾值的范圍[umin,umax]，對應的m個閾值分別是[umin,u2,…,um-1,umax]；

(3)如果閾值為連續(xù)的值，則選擇ud的90％值作為閾值的最小值，選擇ud的110％值作為閾值的最大值，確定了閾值的范圍[umin,umax]＝[90％×ud,110％×ud]，并將此范圍內的閾值等間距的分為m個離散的閾值，即[90％×ud,u2,…,um-1,110％×ud]。

技術方案中所述對超越量進行gpd擬合是指：

從閾值范圍[umin,umax]分別選擇m個不同的閾值，確定超過閾值的載荷極值樣本，針對該閾值對應的超越量進行gpd擬合，選擇極大似然估計法進行gpd的參數(shù)估計。

技術方案中所述擬合優(yōu)度檢驗是指：

采用χ²檢驗或者ks檢驗對gpd擬合進行檢驗，對未通過檢驗的分布直接淘汰，通過檢驗的進行下一步計算。

技術方案中所述計算采樣分布和擬合分布序列的d1,d2值，具體步驟如下：

(1)假設n個極值樣本值t1,t2,…,tn，采樣分布序列fn＝{fn(t1),fn(t2),…,fn(tn)}，其中fn(ti)由中位秩公式得到，fn(ti)＝(i-0.3)/(n+0.4)；

(2)擬合分布序列fo＝{fo(t1),fo(t2),…,fo(tn)}；

(3)計算采樣分布和擬合分布序列的d1，d2值

技術方案中所述求增量序列：

δy1＝{δy1(tk)＝[fn(tk)-fn(tk-1)]/δtk,k＝2,3,…,n}(δtk＝tk-tk-1)

δy2＝{δy2(tk)＝[f0(tk)-f0(tk-1)]/δtk,k＝2,3,…,n}

技術方案中所述計算時間段的關聯(lián)系數(shù)是指：

為了辨別關聯(lián)系數(shù)變化趨勢的方向，對公式

引入了符號函數(shù)，則關聯(lián)系數(shù)計算如下：

技術方案中所述計算灰關聯(lián)度，是用下列公式計算：

灰關聯(lián)度r(x1,x2)越大，表明樣本分布與擬合分布的曲線態(tài)勢越接近，則選擇灰關聯(lián)度最大的擬合分布對應的閾值作為pot外推中最佳的閾值。

與現(xiàn)有技術相比本發(fā)明的有益效果是：

1.傳統(tǒng)的馬爾可夫鏈蒙特卡羅(markovchainmontecarlo，簡稱mcmc)的載荷時域擴展方法并未考慮相鄰極值間的間隔時間，一般認為相鄰極值間的間隔時間為采用頻率的倒數(shù)，但這明顯與實際的時域載荷情況不符，因為采集的載荷時域信號中有很多無效的幅值數(shù)據(jù)，通過信號降噪和數(shù)據(jù)壓縮之后保留了有效的極值數(shù)據(jù)，如果還認為相鄰極值間的時間間隔為采用頻率的倒數(shù)，則會改變相同時間內時域載荷中的循環(huán)數(shù)量。針對數(shù)控機床切削加工過程，由于每種工藝參數(shù)組中的主軸轉速和進給速度的不同，所以刀具切削工件一周的時間不同。因此，針對切削載荷研究載荷相鄰極值間的間隔時間，才能夠準確的對時域載荷進行擴展。

2.閾值選擇的是否恰當可以等效為分析超越量gpd擬合的好壞，即對可以通過擬合優(yōu)度檢驗的優(yōu)良來判斷選擇的閾值是否恰當。然而目前針對gpd分布卻沒有一種公認的最優(yōu)擬合優(yōu)度檢驗。本發(fā)明專利通過灰關聯(lián)分析法對比分析不同閾值對應的超越量的gpd擬合的灰關聯(lián)度，選擇灰關聯(lián)度最大的gpd擬合對應的閾值為最優(yōu)的閾值。

附圖說明

下面結合附圖對本發(fā)明作進一步的說明：

圖1為本發(fā)明所述的一種基于時域載荷擴展的數(shù)控機床載荷外推方法的流程圖；

圖2為本發(fā)明所述的一種基于時域載荷擴展的數(shù)控機床載荷外推方法中的提取時域載荷的極值及升程的示意圖；

圖3為本發(fā)明所述的一種基于時域載荷擴展的數(shù)控機床載荷外推方法中的實測載荷時間歷程與擴展1倍的載荷時域歷程的對比圖；

圖4為本發(fā)明所述的一種基于時域載荷擴展的數(shù)控機床載荷外推方法中的mef圖及確定的閾值范圍和最佳閾值的示意圖；

圖5為本發(fā)明所述的一種基于時域載荷擴展的數(shù)控機床載荷外推方法中的閾值及超越量的示意圖；

具體實施方式

下面結合附圖對本發(fā)明作詳細的描述：

參閱圖1，本發(fā)明所述的一種基于時域載荷擴展的數(shù)控機床載荷外推方法由基于mcmc的時域載荷擴展法、利用mef圖確定閾值的范圍、基于灰關聯(lián)分析的最優(yōu)閾值的選取以及獲取外推的時域載荷等步驟組成。

(1)基于mcmc的時域載荷擴展法

step1提取時域載荷的極值及升程

參閱圖2，提取數(shù)控機床實測載荷時間歷程的轉折點作為極值載荷時間序列，在提取載荷轉折點時，需要剔除無效的載荷循環(huán)。采用變程閾值模型進行取舍。剔除完無效的載荷循環(huán)之后，提取時域載荷信號中的轉折點{xt}，即極大值點和極小值點。

式中，m1，m2，…,mn代表局部極小值點，m1，m2，…,mn代表局部極大值點，n表示極值點的個數(shù)。

從轉折點{xt}中提取相鄰極值點的差值及對應相鄰極值點的間隔時間，將第t個相鄰極值點的差值記為升程st，對應的相鄰極值點的間隔時間記為tt，則第t個升程和對應的間隔時間可以記為tst＝(tt,st)；

step2劃分載荷等級水平

將提取的載荷極大值和極小值通過固定的載荷等級水平劃分為nd組，通常情況下劃分為32組或者64組，即極值點有32或64個狀態(tài)。

將升程st通過固定的載荷等級水平也劃分為nd組，而間隔時間tt不需要進行分組，因為間隔時間為離散值，是采用頻率的倒數(shù)(1/f)的整數(shù)倍。

step3統(tǒng)計載荷循環(huán)次數(shù)及升程

統(tǒng)計轉折點{xt}中從時域載荷的極大值到極小值和極小值到極大值狀態(tài)轉移的循環(huán)次數(shù)，將循環(huán)矩陣分別記為f和f’；

統(tǒng)計求解出每個狀態(tài)的升程對應不同間隔時間的頻次，將升程st及對應的間隔時間tt的頻次矩陣記為f”。

step4求解mcmc的狀態(tài)轉移矩陣

求解載荷極大值到極小值和極小值到極大值的狀態(tài)轉移矩陣。將循環(huán)矩陣f(i,j)和f’(j,i)的各級頻數(shù)除以各級頻數(shù)的總和就可以得到狀態(tài)轉移概率矩陣p和p’。

再將升程及間隔時間的頻次矩陣f”(j,i)的各級頻數(shù)除以總的頻數(shù)就是狀態(tài)轉移概率矩陣p”。

step5求解以上三個狀態(tài)轉移概率矩陣的累積概率矩陣pcum，p'cum，p”cum

step6計算馬爾科夫鏈的平穩(wěn)分布π(xt)

step7在(0,1)區(qū)間內隨機產(chǎn)生符合均勻分布的隨機數(shù)，進行載荷狀態(tài)的判斷，確定載荷水平等級。

例如，將產(chǎn)生的隨機數(shù)跟累計概率矩陣pcum中的第i行元素分別進行對比，假設其位于第j-1和j之間的狀態(tài)，則可以確定極大值m1所處的載荷等級n1，再次在(0,1)區(qū)間產(chǎn)生隨機數(shù)，再將隨機數(shù)與累計概率矩陣p’cum中的第j行元素分別進行對比，確定極小值m2所處的載荷等級n2；然后再生產(chǎn)隨機數(shù)，確定下一個極大值m2所處的載荷等級n3，以此類推，即可以獲得載荷的狀態(tài)等級序列n。

step8生成載荷時間歷程的轉折點序列

依據(jù)載荷狀態(tài)等級序列生成載荷時間歷程的轉折點序列，原始載荷序列xt，載荷劃分的等級數(shù)nd，載荷的狀態(tài)等級序列n，根據(jù)式1.9可以生成載荷時間歷程的轉折點序列。

step9確定相鄰極值點之間的間隔時間

根據(jù)生成的載荷時間歷程的轉折點序列，計算出相鄰極值點的差值，即升程st。在(0,1)區(qū)間內隨機產(chǎn)生符合均勻分布的隨機數(shù)，根據(jù)升程st，將隨機數(shù)與p”cum的第i行元素進行比較，如果這個隨機數(shù)位于狀態(tài)j-1和j之間，則可以確定該升程st對應的間隔時間tt，以此類推，確定出所有相鄰極值點對應的間隔時間，完成了時域載荷的擴展，參閱圖3，對實測的60秒切削力載荷時間歷程信號與擴展1倍的載荷時域歷程的對比圖，其中前60秒為實測的60秒切削力載荷時間歷程信號，后60秒為擴展1倍的載荷時域歷程信號。

(2)利用mef圖確定閾值的范圍

step1繪制載荷mef圖

定義e(u)＝e(x–u|x>u)為載荷極值x的超越量均值函數(shù)(meanexcessfunction，簡稱mef)，載荷極值樣本經(jīng)驗超越量均值函數(shù)為：

參閱圖4，繪制由超越量均值函數(shù)值e(u)與極值閾值u的散點圖構成的曲線，如果在大于某一閾值ud之后，mef的曲線接近于斜率為正的直線，則載荷極值樣本服從形狀參數(shù)ξ＞0的gpd分布；如果直線段的斜率為負，表示該載荷極值樣本為薄尾分布，服從尾部較短的gpd分布；如果直線段為水平線，則表明載荷極值樣本遵從指數(shù)分布。

step2確定閾值的范圍[umin,umax]及劃分成m個閾值

根據(jù)繪制的mef圖，確定mef曲線接近為直線部分的最左端對應的閾值ud。

如果閾值為離散的值，則以閾值ud為中心，從閾值ud相鄰的閾值中選出m個閾值，構成閾值的范圍[umin,umax]，對應的m個閾值分別是[umin,u2,…,um-1,umax]，一般情況下m的取范圍在8到12之間。

如果閾值為連續(xù)的值，則選擇ud的90％值作為閾值的最小值，選擇ud的110％值作為閾值的最大值，確定了閾值的范圍[umin,umax]＝[90％×ud,110％×ud]，并將此范圍內的閾值等間距的分為m個離散的閾值，即[90％×ud,u2,…,um-1,110％×ud]。

如表1所示，對一組實測的車削力載荷信號進行時域載荷外推，確定的極大值的閾值范圍為[396.4,412.5]，極小值的閾值范圍為[156.4,172.5]，并將閾值根據(jù)實際載荷極值分成m＝8個閾值。

(3)基于灰關聯(lián)分析的最優(yōu)閾值的選取

step1對不同閾值u對應的超越量進行廣義帕累托分布(generalizedparetodistribution，簡稱gpd)擬合

參閱圖5，從閾值范圍[umin,umax]分別選擇m個不同的閾值，確定超過閾值的載荷極值樣本，針對該閾值對應的超越量進行gpd擬合，選擇極大似然估計法進行gpd的參數(shù)估計。

step2擬合優(yōu)度檢驗

采用χ²檢驗或者ks檢驗對gpd擬合進行檢驗，對未通過檢驗的分布直接淘汰，通過檢驗的進行下一步計算。

step3計算采樣分布和擬合分布序列的d1，d2值

假設在[a,b]時間內，有n個極值樣本值t1,t2,…,tn，采樣分布序列fn＝{fn(t1),fn(t2),…,fn(tn)}，其中fn(ti)由中位秩公式得到，fn(ti)＝(i-0.3)/(n+0.4)；擬合分布序列fo＝{fo(t1),fo(t2),…,fo(tn)}，其中fo(ti)是由試驗數(shù)據(jù)求得的分布函數(shù)的對應值。

step4求增量序列

δy1＝{δy1(tk)＝[fn(tk)-fn(tk-1)]/δtk,k＝2,3,…,n}(δtk＝tk-tk-1)

δy2＝{δy2(tk)＝[f0(tk)-f0(tk-1)]/δtk,k＝2,3,…,n}

step5計算時間段的關聯(lián)系數(shù)

為了辨別關聯(lián)系數(shù)變化趨勢的方向，引入了符號函數(shù)，則關聯(lián)系數(shù)計算如下：

step6計算灰關聯(lián)度

灰關聯(lián)度r(x1,x2)越大，表明樣本分布與擬合分布的曲線態(tài)勢越接近，則選擇灰關聯(lián)度最大的擬合分布對應的閾值作為pot外推中最佳的閾值。

如表1所示，表1為本發(fā)明所述的一種基于時域載荷擴展的數(shù)控機床載荷外推方法中的一組實測車削力載荷pot外推的閾值及對應的灰關聯(lián)度。

表1

對一組實測的車削力載荷確定的8個極大值閾值和極小值閾值，分別計算每個閾值對應的灰關聯(lián)度，并選擇最大的灰關聯(lián)度對應的閾值為最佳閾值，最終確定極大值最優(yōu)閾值為405.6n，極小值最優(yōu)閾值為163.3n，極大值和極小值最優(yōu)閾值對應的超越量的gpd擬合參數(shù)估計值如表2所示。

表2為本發(fā)明所述的一種基于時域載荷擴展的數(shù)控機床載荷外推方法中的一組實測車削力載荷最優(yōu)閾值對應的超越量的gpd擬合參數(shù)的估計值。

表2

(4)獲取外推的時域載荷

根據(jù)超越量的數(shù)量，模擬產(chǎn)生服從最優(yōu)gpd的新超越量，并用新產(chǎn)生的超越量來代替原有的超越量，從而獲取外推的載荷時域信號。