貝葉斯優(yōu)化TQWT 參數(shù)在軸承故障診斷中的應用

張樂，彭先龍，朱華雙

（西安科技大學 機械工程學院，西安 710054）

軸承是旋轉(zhuǎn)機械的核心部件，研究早期故障信號特征，識別故障類型，采取合理的應對措施，對維護設備安全穩(wěn)定運行十分重要。

可調(diào)品質(zhì)因子小波變換（Tunable Q-factor wavelet transform，TQWT）是Selesnick 提出的一種參數(shù)可控可調(diào)節(jié)的過完備小波變換[1]，在機械故障診斷中得到廣泛應用。其優(yōu)點是通過調(diào)節(jié)品質(zhì)因子Q、冗余度r和分解層數(shù)J這3 個關鍵參數(shù)，可以靈活地構造小波基函數(shù)，匹配求解旋轉(zhuǎn)機械振動信號中存在的周期性沖擊特征。為了合理地確定（Q、r、J）參數(shù)，孔運等[2]利用故障特征評價指標，以提升子帶信號故障特征顯著性為目標，遍歷搜索參數(shù)組成空間，通過檢測和整合不同參數(shù)組合條件下故障特征指標變化趨勢，確定合理參數(shù)，構造自適應可調(diào)品質(zhì)因子小波。由于三元參數(shù)存在龐大的候選組合，為提高網(wǎng)格搜索執(zhí)行的效率，賀王鵬等[3]以（Q,r）為小波基函數(shù)自適應性調(diào)控主導參數(shù)，構建網(wǎng)格搜索空間。任學平等[4]在分析已有研究結果的基礎上，通過調(diào)節(jié)（Q,J）參數(shù)組合獲得信號最佳分解效果。合理優(yōu)化搜索空間提高了可用參數(shù)的搜索效率。

應用優(yōu)化算法在故障特征指標導引下實現(xiàn)參數(shù)自動尋優(yōu)是目前研究的方向，已在變分模態(tài)分解關鍵參數(shù)自適應優(yōu)化研究方面，取得了大量研究成果。例如：唐貴基等[5]采用粒子群算法以包絡熵為優(yōu)化目標，獲得分解個數(shù)k和懲罰參數(shù)α最優(yōu)解。劉嘉敏等[6]提出應用遺傳算法對目標函數(shù)在解空間進行全局并行隨機搜索，求取二元參數(shù)（k,α）最佳組合。王振亞等[7]以最小包絡均熵為優(yōu)化目標，用樽海鞘群算法實現(xiàn)變分模態(tài)分解的參數(shù)尋優(yōu)。通過網(wǎng)格搜索和優(yōu)化算法實現(xiàn)參數(shù)尋優(yōu)，需要依靠評估函數(shù)分析某樣本參數(shù)下，故障信號TQWT 分解獲得子帶的故障特征顯著性。因評估函數(shù)無法用簡單函數(shù)表示，所以樣本參數(shù)評估代價大，在分析TQWT 三元參數(shù)時，遍歷樣本會耗費大量的計算資源。因此，尋找高效的優(yōu)化算法是十分必要的。

貝葉斯優(yōu)化是一種基于概率代理模型的序貫優(yōu)化，可以有效地利用歷史觀測信息主動選擇評估點，減少尋優(yōu)迭代次數(shù)[8-10]，適用于求解存在觀測噪音、多峰、非凸、黑箱及評估代價高等問題的優(yōu)化目標，廣泛用于機器學習和深度學習超參數(shù)優(yōu)化[11-12]。本文針對網(wǎng)格搜索和算法優(yōu)化實現(xiàn)TQWT 調(diào)參過程中目標函數(shù)評估代價高的問題，提出基于貝葉斯優(yōu)化TQWT 參數(shù)的故障診斷算法，以平均香農(nóng)熵和峭度聯(lián)合建立故障特征評價指標，在故障特征指標導引下提取最顯著故障特征子帶信號，經(jīng)TQWT 逆變換后由重構信號包絡譜判別軸承故障類型。

1 可調(diào)品質(zhì)因子小波變換（TQWT）

TQWT 采用單位離散傅里葉變換求得輸入信號的頻域表示，經(jīng)由一系列高、低通濾波器，不斷對低尺度濾波器的輸出進行雙通道迭代分解，獲得多尺度小波系數(shù)序列，其低、高通濾波器頻響函數(shù)H0（ω）和H1（ω）定義為：

式中：α為低通尺度變換系數(shù)；β為高通尺度變換系數(shù)；θ（ω）為函數(shù)過渡區(qū)間設計函數(shù)。

具有2 階消失矩的德比契斯（Daubechies）頻率響應函數(shù)可定義為

可以證明θ（ω） 在過渡頻帶上滿足

則低、高通濾波器頻響函數(shù)H0（ω）和H1（ω）在各自通帶和過渡帶上滿足如下關系

且為了保證完美重構，α與β應滿足

尺度內(nèi)濾波器參數(shù)（α,β）可由TQWT 品質(zhì)因子Q和冗余度r確定[2]，即

因此，TQWT 多尺度分解與重構由三元參數(shù)（Q，r，J）決定，其中J為分解層數(shù)。為了節(jié)約算力，同時避免過渡分解導致的頻帶信息冗余，分解層數(shù)J存在理論上限Jmax，計算公式為

式中：N為待分析信號長度；“”為向下取整運算符。

2 貝葉斯優(yōu)化TQWT 參數(shù)

TQWT 算法匹配故障特征構造可調(diào)品質(zhì)因子小波，可以看作一個優(yōu)化問題，即

式中：x為一個三維決策向量，由TQWT 算法預設參數(shù)（Q,r,J）定義x的決策空間Ω；xbest是最優(yōu)參數(shù)設計，f（·）為優(yōu)化目標函數(shù)。如優(yōu)化目標為最大化問題，可對目標函數(shù)取負號轉(zhuǎn)化求解。

2.1 概率代理模型

貝葉斯優(yōu)化通常采用高斯過程建立非參數(shù)模型。以目標函數(shù)f（x）評估TQWT 在某參數(shù)x下分解信號所得子帶故障特征的顯著性，會得到確定的計算值。為滿足高斯過程條件，在目標函數(shù)f（x）計算中添加觀察噪聲得到觀察值

式中ε為觀察誤差，ε～N（0,σ2）。

目標函數(shù)f（x）作為觀察對象，其對應的觀察值y為服從正態(tài)分布的隨機變量，則“觀察”某參數(shù)條件下 TQWT 子帶故障特征成為具有一定成功概率的獨立事件，滿足高斯過程，可通過概率模型估計目標函數(shù)f（x），即

式中：GP為高斯分布；m（x） 為均值函數(shù)，m（x）=E[f（x）]；k（x,x'）為協(xié)方差函數(shù),k（x,x'）=E{[f（x）-m（x）][f（x'）-m（x'）]}。

由高斯過程性質(zhì)推導可得概率代理模型預測公式為

2.2 基于預期改進的采集函數(shù)

當決策空間進行了n次采樣，目標函數(shù)最小值的改進量定義為

新觀察位置x應滿足目標函數(shù)f（x）改進量I（x）取最大值，即

在選擇下一觀察位置之前，新位置的目標函數(shù)值未知，可由歷史觀察數(shù)集y, 根據(jù)f（x）后驗分布估計預期的提升量EIn（x）（Expected improvement，EI），則式（15）更新為

式中En[·]為函數(shù)f（x）在已給定的n次采樣下所得后驗分布的期望。

由式（12）可知f（x）后驗分布服從均值為μ(x)和方差為σ2(x)的正態(tài)分布，求得期望[13]為

式中：Δn(x)為目標函數(shù)在點x的歷史最優(yōu)值與期望值之間的差值，Δn(x)=-μ(x) ； Φ(·)為標準正態(tài)分布的累積函數(shù)；φ(·)為標準正態(tài)分布的概率密度函數(shù)。使用1 階或2 階優(yōu)化算法可以求解式（17）中滿足式（16）預期改進量EIn（x）最大化的解，即新的“潛在”觀察評估點。

2.3 優(yōu)化目標函數(shù)

貝葉斯優(yōu)化TQWT 參數(shù)需要有效的故障特征導引。峭度對信號振動變化十分敏感，能反映信號中蘊含的沖擊特征，其值越大信號的沖擊特征就越顯著，廣泛用于顯著故障特征子帶抽取[14-15]，定義為

香農(nóng)熵可以用于故障特征評價，振動信號周期性越明顯，熵值越低；噪聲干擾越多，熵值越大[16]，但是TQWT 分解所得小波系數(shù)序列的長度不同，經(jīng)典香農(nóng)熵無法評估此類小波系數(shù)序列，可采用平均香農(nóng)熵，即:

在參考已有綜合特征指標建立方法的基礎上，本文以平均香農(nóng)熵與峭度之比定義熵-峭指標，評價信號的故障特征顯著性，其定義為

式中：Ej和Kurtj分別為尺度j下子帶信號對應的平均香農(nóng)熵和峭度。由熵-峭指標EK定義可知，EK值越小，沖擊特征越顯著，據(jù)此建立的綜合目標函數(shù)為

式中：x為TQWT 三元參數(shù)（Q,r,J）定義的決策向量；EK（x）j為參數(shù)x作用下尺度j子帶信號熵-峭指標計算值。目標函數(shù)f（x）以多尺度子帶熵-峭指標最小值評估當前參數(shù)組合x的性能。

貝葉斯優(yōu)化TQWT 參數(shù)的流程如圖1 所示。

圖1 貝葉斯優(yōu)化過程Fig.1 Bayesian optimization process

具體步驟如下：

步驟1 根據(jù)TQWT 算法三元參數(shù)（Q,r,J）取值范圍初始化參數(shù)集X={x1,x2, ···,xn}，由式（10）和式（22）得到X對應的故障特征觀測數(shù)據(jù)集y={y1,y2, ···,yn}。

步驟2 根據(jù)式（12）更新目標函數(shù)f（x）高斯代理模型。

步驟3 根據(jù)式（16）求解滿足采集函數(shù)EIn（x）最大化的下一“潛在”觀察目標位置xnew。

步驟4 由式（10）和式（22）計算新觀察點xnew觀察值ynew，加入歷史觀測數(shù)據(jù)集。若新點xnew目標函數(shù)值f（xnew）取得最小值，則更新當前最優(yōu)位置xbest=xnew。

步驟5 重復執(zhí)行步驟2 ～ 步驟4，循環(huán)迭代指定次數(shù)，輸出最優(yōu)位置xbest。

與粒子群優(yōu)化、遺傳算法等進化方法相比，貝葉斯優(yōu)化算法可以顯著提高優(yōu)化求解效率，降低計算資源消耗。進化方法在迭代過程中需要更新種群內(nèi)部全體或部分粒子（個體）。如文獻[5]中粒子群優(yōu)化通過迭代更新每個粒子的速度和位置搜索參數(shù)組合空間。文獻[16] 中遺傳算法采用單點交叉和變異，每代中有1/4 的個體被淘汰。粒子（個體）每次更新后都要重新計算對應的適應度值，而評估更新粒子（個體）的適應度值會耗費一定的計算資源，當種群規(guī)模擴大時會影響算法的分析性能。

貝葉斯優(yōu)化算法更新采樣點的策略與進化方法不同，在迭代過程中根據(jù)概率代理模型，將采集函數(shù)預期改進最優(yōu)的位置作為下一觀察點，提高了趨向目標函數(shù)極值更新參數(shù)的準確性，且每次迭代只需評估一個新位置的目標函數(shù)值，減少了算法耗費的計算資源。通過初始觀察數(shù)據(jù)集建立的概率統(tǒng)計模型隨迭代會不斷更新，提高預測“潛在”觀察位置的準確性， 使優(yōu)化目標函數(shù)快速收斂穩(wěn)定，可在更新少量樣本的條件下求得可用的優(yōu)化參數(shù)解。

3 仿真信號分析

為驗證貝葉斯優(yōu)化算法利用本文熵-峭指標優(yōu)化TQWT 參數(shù)提取軸承故障特征的有效性，模擬構建一組含噪周期沖擊信號。仿真信號由周期沖擊信號、諧振干擾和背景噪聲疊加構成，可以表示為

式中：I（t）為沖擊振動信號；h（t）為系統(tǒng)諧振干擾；n（t）為工作背景噪音。

仿真信號及各組成成分如圖2 所示，其采樣頻率為12 kHz，采樣點數(shù)12 000 點。

圖2 仿真信號時域波形及組成成分Fig.2 Time domain waveform and components of simulated signals

各組成部分如下：

1）周期性沖擊振動信號I（t），由單點沖擊信號延時疊加構成，即

式中： ζ為阻尼系數(shù)，ζ=0.02 ； ω為系統(tǒng)共振頻率，ω=2πfn，fn=2 000 Hz；A為振動的幅值，A=1.2；T為沖擊周期，T=0.05 s。

2）諧波干擾信號h（t）頻率為50 Hz，則

3）背景噪音n（t）選擇服從均值μ為0，方差σ為0.6 的高斯白噪聲。合成仿真信號噪比為-13 dB。從圖2 中可以看出，合成仿真信號受到諧振和噪聲的嚴重污染，已無法直接觀察振動信號的變化規(guī)律。原信號所含周期性沖擊特征在圖3 包絡譜中無法有效識別。應用貝葉斯優(yōu)化算法優(yōu)化TQWT參數(shù)，設置品質(zhì)因子和冗余度的初始值QL和rL均為2，上限QU和冗余度rU分別為8 和6。分解層數(shù)的初始值JL為5，對應上限JU根據(jù)最大分解層數(shù)公式（8）確定，由QL和rL計算得到JU為17。后續(xù)實驗采用相同設置，如需要增大分解層數(shù)上限，需要匹配增大品質(zhì)因子Q和冗余度r的下限值。

圖3 仿真信號Hilbert 包絡譜Fig.3 Hilbert envelope spectrum of simulated signal

根據(jù)圖1 所示的貝葉斯優(yōu)化算法執(zhí)行100 次迭代計算，目標函數(shù)值變化過程如圖4 所示，在40 次迭代運算后趨于穩(wěn)定，得到TQWT 參數(shù)組合為（Q= 7,r= 2,J= 13）。利用文獻[2]TQWT 濾波算法分解仿真信號，以本文熵-峭指標評估TQWT 分解獲得子帶信號故障特征的顯著性，結果如圖5 中熵-峭指標計算值變化曲線所示，確定的最顯著故障特征子帶編號j=8（最小值）。

圖4 目標函數(shù)變化過程Fig.4 Objective function variation processes

圖5 最顯著特征子帶選擇結果Fig.5 The most significant feature subband selection results

選擇最顯著故障特征子帶進行單支重構，得到的故障特征信號波形如圖6 所示，圖7 為其包絡譜。比較圖6 中故障特征信號與原始沖擊信號，可以明顯觀察到周期為0.05 s 的沖擊特征。圖7 包絡譜中可以清晰辨認沖擊故障特征頻率1 倍頻（實際19.7 Hz）及高倍頻成分。圖7 中：“f”為曲線對應的理論故障特征頻率，“1f”為1 倍頻，“2f”為2 倍頻，其余及后續(xù)圖示采用相同標記。因此，通過貝葉斯優(yōu)化算法合理確定TQWT 參數(shù)，在熵-峭指標導引下，選擇最顯著特征子帶信號重構后進行包絡譜分析，可以成功識別噪聲混疊狀態(tài)下的周期性沖擊特征。

圖6 提取故障特征時域波形與原沖擊特征對比Fig.6 Comparison of extracted fault characteristic time domain waveform with original impact characteristics

圖7 提取故障特征Hilbert 包絡譜（f=20 Hz）Fig.7 Extracting Hilbert envelope spectrum of fault characteristics (f=20 Hz)

4 實測故障信號分析

4.1 軸承故障信號分析

為驗證本方法分析實測機械故障振動信號的有效性，采用美國凱斯西儲大學（CWRU）滾子軸承故障檢測數(shù)據(jù)進行驗證。軸承型號6 205-2RS JEM SKF 深溝球軸承，電火花加工模擬單點故障直徑為0.177 8 mm，深度0.279 4 mm，電機轉(zhuǎn)速1 791 r/min（轉(zhuǎn)頻fr=30 Hz），采樣頻率fs=12 000 Hz，數(shù)據(jù)文件編號IR007_0（內(nèi)圈故障），其故障信號時域波形如圖8 所示。

圖8 軸承故障信號時域波形Fig.8 Time domain waveform of bearing fault signal軸承內(nèi)圈和外圈故障特征頻率[17]計算公式為：

式中：Z為滾珠個數(shù)；d為滾珠直徑；D為軸承滾道節(jié)徑；α為軸承接觸角，深溝球軸承為0°；fin和fout分別為軸承內(nèi)圈和外圈故障特征頻率；fr為軸承內(nèi)外圈相對轉(zhuǎn)頻。查閱軸承結構參數(shù)，根據(jù)式（26）可以得到對應內(nèi)圈故障特征頻率fin=162.2 Hz。

從圖8 原始振動信號時域波形可以看出信號波形變化復雜，含有大量的噪聲干擾。采用本文貝葉斯優(yōu)化算法，經(jīng)過圖9 所示的100 次迭代計算求得TQWT 參數(shù)組合為（Q=6,r=4,J=15）。

圖9 目標函數(shù)變化過程Fig.9 Objective function variation process of inner ring fault signals

利用熵-峭指標評估TQWT 分解獲得子帶信號故障特征的顯著性，結果如圖10 所示。確定最顯著故障特征子帶編號j=4，由TQWT 逆變換進行單支重構獲得降噪特征信號及其包絡譜，如圖11 所示。

圖10 最顯著特征子帶選擇結果Fig.10 The most significant feature subband selection results of inner ring fault signals

圖11 故障特征提取結果Fig.11 Fault feature extraction results of inner ring fault signals

將圖8 原始故障信號和圖11a）提取故障特征信號時域波形相比較，可知算法降噪效果良好，提取特征信號的周期性沖擊特征顯著。分析圖11b）特征信號包絡譜可以準確獲得故障特征頻率1 倍頻（實際161.9 Hz）及其2 倍頻和3 倍頻，還可以觀察轉(zhuǎn)頻fr（實際30.0 Hz）。試驗結果與理論分析結果一致，揭示了原故障信號中包含的周期性沖擊特征，據(jù)此可準確識別軸承故障類型。

4.2 全壽命信號早期故障分析

為驗證貝葉斯優(yōu)化TQWT 濾波算法與熵-峭指標提取早期軸承故障所含微弱沖擊特征的能力，對NSF I/UCR智能維護系統(tǒng)中心（IMS）滾動軸承實驗數(shù)據(jù)進行了分析研究。IMS 軸承數(shù)據(jù)共記錄了3 組滾動軸承全壽命振動加速度信號，選擇第2 組1 號軸承實驗數(shù)據(jù)（外圈故障），位置如圖12 所示。軸承型號ZA-2115，采樣頻率fs= 20 000 Hz，采樣間隔10 min，持續(xù)時間164 h，共采集數(shù)據(jù)文件984 個，每個文件長度為20 480 個采樣點，電機轉(zhuǎn)速2 000 r/min（轉(zhuǎn)頻fr= 33 Hz）。查閱軸承結構參數(shù)，根據(jù)式（27）得到此軸承外圈故障特征頻率fout= 235.3 Hz。

振動信號均方根值的變化趨勢可以反映軸承的工作狀態(tài)。圖13 為1 號軸承在全壽命工作過程中均方根值的變化。從圖中可以看出，軸承在5 000 min之前運行穩(wěn)定，到5 200 min 后均方根值開始逐漸增大，約至7 000 min 時發(fā)生突變，標志軸承狀態(tài)發(fā)生改變，7 000 ～ 9 000 min 時段均方根值波動上升，運行9 000 min 后均方根值開始大幅增加，在9 800 min時取得最大值，軸承壽命達到極限。試驗選擇分析3 600 min 時刻振動數(shù)據(jù)，此時故障仍處于潛伏發(fā)展階段，故障信號包絡譜分析無法直接有效提取振動數(shù)據(jù)中可能存在的故障特征[5]。

圖13 1 號軸承振動信號均方根值變化趨勢Fig.13 Root mean square (RMS) value variation trend of vibration signals of bearing 1

為了準確提取早期故障中存在微弱振動信號，采用本文貝葉斯優(yōu)化算法，執(zhí)行500 次迭代運算，確定的TQWT 參數(shù)組合為（Q=3,r=6,J=11），如圖14所示。

圖14 目標函數(shù)變化過程Fig.14 Objective function variation process of outer ring fault signals

利用熵-峭指標評估TQWT 分解子帶信號故障特征的顯著性，結果如圖15 所示，確定最顯著故障特征子帶編號j=1，通過單支重構獲得的故障特征信號及其包絡譜如圖16 所示。由圖16a） 可以發(fā)現(xiàn)提取故障特征信號時域波形的沖擊特征變化具有較強的規(guī)律性，在圖16b）包絡譜中可以準確獲得其故障特征頻率主頻峰值點（實際231.0 Hz）。理論分析與實驗分析結果一致，但是倍頻特性與其他頻率成分混雜，無法清晰辨識，所得的主頻峰值點可作為早期故障預判的支撐數(shù)據(jù)。

圖15 最顯著特征子帶選擇結果Fig.15 The most significant feature subband selection results of outer ring fault signals

圖16 故障特征提取結果（f = 235 Hz）Fig.16 Fault feature extraction results of outer ring (f = 235 Hz)

應用本文貝葉斯優(yōu)化算法，分別以平均香農(nóng)熵和峭度建立優(yōu)化目標函數(shù)，分析相同時刻（3 600 min）軸承振動數(shù)據(jù)，執(zhí)行500 次迭代求取合理設置參數(shù)，在獨立特征指標導引下分析早期故障特征，如圖17a）和圖18a）所示。

圖17 以平均香農(nóng)熵為目標函數(shù)提取故障特征Fig.17 Utilizing average Shannon entropy as objective function for fault feature extraction

圖18 以峭度為目標函數(shù)提取故障特征Fig.18 Extracting fault features with kurtosis as objective function

以平均香農(nóng)熵作為優(yōu)化目標得到的TQWT 參數(shù)組合為（Q=7,r=6,J=12），其對應最顯著特征子帶編號j= 11。當以峭度作為優(yōu)化目標時，因子帶信號故障特征對應峭度最大值，在算法求解時需要取負號建立目標函數(shù)，得到的TQWT 參數(shù)組合為（Q=8,r=6,J=16），其對應最顯著特征子帶編號j= 13。分別選擇兩種目標函數(shù)對應的最顯著特征子帶進行單支重構，獲得對應故障特征信號的包絡譜如圖17b）和圖18b）所示，從包絡譜圖中無法定位軸承故障特征頻率及其他相關頻率信息。因此，基于平均香農(nóng)熵和峭度分別建立目標函數(shù)，無法有效提取軸承早期故障特征，熵-峭指標對早期軸承故障特征的適應性優(yōu)于獨立特征評價指標。

5 結論

貝葉斯優(yōu)化算法近些年在多領域研究中得到廣泛應用，本文將貝葉斯優(yōu)化算法與可調(diào)品質(zhì)因子小波變換（TQWT）相結合，提出一種基于貝葉斯優(yōu)化TQWT 參數(shù)的故障診斷算法，得到以下結論：

1）貝葉斯優(yōu)化算法通過概率代理模型和采集函數(shù)提高了樣本迭代更新的效率，在更新少量樣本的條件下目標函數(shù)值快速收斂穩(wěn)定，降低了算法執(zhí)行過程中樣本總體的評估代價。

2）本文以平均香農(nóng)熵與峭度之比定義熵-峭指標，評估信號故障特征顯著性，在小波變換域上可準確識別最顯著故障特征子帶。實測信號分析表明熵-峭指標對早期故障特征的適應性優(yōu)于獨立評價指標。

3）綜合仿真實驗與實測軸承振動信號分析結果表明，利用貝葉斯優(yōu)化TQWT 參數(shù)，采用熵-峭指標引導最顯著故障特征子帶選擇，可以提取軸承早期故障中存在的微弱故障特征，為故障診斷和早期預警提供有價值的支撐數(shù)據(jù)。