基于狀態(tài)/參數(shù)交替優(yōu)化估計的齒輪時變嚙合剛度辨識方法研究

摘要： 齒輪副嚙合剛度是齒輪傳動動力學(xué)模型中的重要參數(shù)，嚙合剛度識別對于齒輪傳動系統(tǒng)振動與噪聲的分析以及狀態(tài)監(jiān)測具有重要意義。由于齒輪傳動系統(tǒng)的時變嚙合剛度很難用直接測量的方法得到，因此有必要發(fā)展基于數(shù)據(jù)驅(qū)動的齒輪時變嚙合剛度辨識方法。為了解決該難題，本文提出了一種系統(tǒng)狀態(tài)和模型參數(shù)交替優(yōu)化方法以辨識齒輪傳動系統(tǒng)時變嚙合剛度。該方法以齒輪嚙合頻率為基頻構(gòu)造傅里葉基函數(shù)庫表征齒輪時變嚙合剛度，并進一步提出了一種基于再生核希爾伯特空間的降噪方法來估計系統(tǒng)狀態(tài)和模型參數(shù)。系統(tǒng)狀態(tài)變量與時變嚙合剛度在動力學(xué)模型與數(shù)據(jù)雙重約束下進行交替迭代優(yōu)化，實現(xiàn)對齒輪傳動系統(tǒng)時變嚙合剛度的辨識。齒輪傳動動力學(xué)仿真和齒輪傳動系統(tǒng)實驗結(jié)果驗證了所提齒輪時變嚙合剛度辨識方法的有效性。

關(guān)鍵詞： 參數(shù)辨識； 齒輪傳動系統(tǒng)； 時變嚙合剛度； 交替優(yōu)化

中圖分類號： O327； TH113.1； TH132.41 文獻標(biāo)志碼： A 文章編號： 1004-4523（2024）10-1707-07

DOI：10.16385/j.cnki.issn.1004-4523.2024.10.008

引 言

齒輪結(jié)構(gòu)被廣泛應(yīng)用于機械傳動系統(tǒng)。在齒輪動力學(xué)模型中，時變嚙合剛度是反映齒輪系統(tǒng)動態(tài)特性的重要參數(shù)，其數(shù)值變化可導(dǎo)致機械傳動系統(tǒng)產(chǎn)生振動和噪聲。對齒輪副嚙合剛度進行建模與識別，有助于齒輪傳動系統(tǒng)的優(yōu)化設(shè)計和狀態(tài)監(jiān)測［1?5］。

確定齒輪副時變嚙合剛度的方法可分為兩種：數(shù)值計算［6?8］和實驗辨識［9?11］。齒輪嚙合剛度數(shù)值計算方法主要包括解析法、有限元法和混合法。解析法基于材料力學(xué)原理對齒輪嚙合剛度進行數(shù)學(xué)建模，如勢能法［12］和石川法［13］。有限元法則在計算中綜合考慮了齒輪幾何結(jié)構(gòu)、輪齒缺陷、材料性能等多種因素的影響。HU等［14］構(gòu)建了含有齒隙的有限元模型，并研究了輪齒修形對嚙合剛度的影響。與解析法相比，有限元法在復(fù)雜齒輪結(jié)構(gòu)仿真中具有更高的計算精度，但其計算時間隨著有限元模型的網(wǎng)格細化而大幅增加。為減少有限元法的計算量并保持其高計算精度，混合法（解析?有限元法）被用以替代有限元法計算齒輪的離散嚙合剛度。在混合法中，齒輪副結(jié)構(gòu)的整體變形由有限元法計算得到，而輪齒接觸區(qū)的局部變形由解析法確定。CHANG等［15］通過解析赫茲解推導(dǎo)齒輪副的局部接觸變形，有效避免了所構(gòu)建有限元模型對輪齒接觸區(qū)域的過度細化。CHEN等［16］在輪齒接觸點附近構(gòu)建了局部剛性區(qū)以避免模型的局部畸變，并利用所構(gòu)建的解析?有限元混合模型分析了不同裂紋擴展路徑對嚙合剛度的影響。實驗辨識一般通過編碼器或成對加速度傳感器預(yù)先測量齒輪系統(tǒng)傳動誤差來估計嚙合剛度，也可應(yīng)用模態(tài)分析、光彈性法、數(shù)字圖像相關(guān)法等技術(shù)辨識齒輪副嚙合剛度。例如RAGHUWANSHI等［11］在全嚙合周期內(nèi)的齒輪副多個嚙合點處進行模態(tài)分析實驗，得到離散的齒輪副嚙合剛度估計。

上述數(shù)值方法在計算嚙合剛度時需要精確的齒輪幾何結(jié)構(gòu)先驗信息，而模態(tài)分析法［11］等實驗辨識法則需要多次改變齒輪副嚙合角進行密集測試，這為齒輪時變嚙合剛度辨識帶來了困難。為了解決上述難題，本文提出一種系統(tǒng)狀態(tài)和模型參數(shù)交替優(yōu)化的辨識方法，該方法可在齒輪結(jié)構(gòu)先驗信息缺失的條件下通過系統(tǒng)輸入和輸出信號辨識齒輪時變嚙合剛度，仿真和實驗結(jié)果驗證了該方法的有效性。

1 齒輪時變嚙合剛度辨識問題

在齒輪傳動系統(tǒng)動力學(xué)模型中，齒輪副間的相互作用由齒輪嚙合界面描述。典型的齒輪傳動動力學(xué)模型如圖1所示，I為轉(zhuǎn)動慣量，r為基圓半徑，θ為旋轉(zhuǎn)角位移，T為扭矩，下標(biāo)1，2分別表示主、從動輪。通過引入等效質(zhì)量和傳動誤差的概念，該模型的動力學(xué)方程可表示為：

（1）

其中：

（2）

式中 m為齒輪副等效質(zhì)量；c和k分別為齒輪副的嚙合阻尼與嚙合剛度；為嚙合力；e為幾何傳動誤差（也稱空載傳動誤差）。

若直接對式（1）所示的齒輪傳動系統(tǒng)模型進行嚙合剛度辨識，需要預(yù)知幾何傳動誤差e。測量幾何傳動誤差要求齒輪系統(tǒng)以極低的轉(zhuǎn)速與扭矩（近似空載，扭矩大小僅使齒輪副在整個嚙合周期恰好保持接觸）運轉(zhuǎn)，從而避免系統(tǒng)動態(tài)特性和齒輪副彈性變形的影響。然而，齒輪副在低轉(zhuǎn)速下的運行平穩(wěn)性往往并不理想，不利于載荷的有效加載和幾何傳動誤差的測量。考慮到低速測量的實際缺陷，本文的輸入、輸出信號在中高轉(zhuǎn)速下獲得：齒輪系統(tǒng)分別在兩次不同的扭矩輸入下運行，則與兩嚙合力（fa ，fb）相對應(yīng)的動力學(xué)方程可表示為：

（3）

（4）

將式（3）與式（4）作差，可得：

（5）

式中 狀態(tài)變量表示兩個不同扭矩下的動態(tài)傳動誤差（Dynamic Transmission Error，DTE）測量值之差。系統(tǒng)輸入表示對應(yīng)嚙合力之差。DTE定義為。

值得注意的是，式（5）對應(yīng)的嚙合剛度辨識問題不再需要預(yù)知幾何傳動誤差e，有效避免了低速空載（或低速輕載）的嚴(yán)苛測量條件。

2 齒輪系統(tǒng)狀態(tài)及參數(shù)交替優(yōu)化

嚙合剛度是齒輪動力學(xué)模型中的時變參數(shù)，其數(shù)值隨著傳動過程中嚙合點的移動而變化。嚙合剛度曲線與輪齒的幾何形態(tài)及齒輪副的嚙合狀態(tài)密切相關(guān)，可以反映齒輪系統(tǒng)的動力學(xué)行為。為從系統(tǒng)輸入、輸出信號中辨識嚙合剛度，本節(jié)提出了一種狀態(tài)/參數(shù)交替優(yōu)化方法以辨識齒輪時變嚙合剛度。

2.1 時變嚙合剛度辨識方法

為辨識齒輪傳動系統(tǒng)時變嚙合剛度，以嚙合頻率為基頻構(gòu)建傅里葉基函數(shù)庫，用來表征齒輪副的時變嚙合剛度與嚙合阻尼：

（6）

假設(shè)有M個采樣點，為辨識系統(tǒng)參數(shù)，根據(jù)方程（5），可采用最小二乘法估計回歸系數(shù)：

（7）

其中：

（8）

然而，最小二乘法容易產(chǎn)生過擬合，為了克服這一缺陷，引入正則項優(yōu)化參數(shù)估計：

（9）

式中 ，為對時間的第i階導(dǎo)數(shù)。為得到嚙合剛度曲線的平滑估計，本文將i取為2。超參數(shù)可通過交叉驗證法［17］確定，并采用自適應(yīng)矩估計（Adaptive Moment Estimation，ADAM）方法［18］求解優(yōu)化問題。

相應(yīng)地，齒輪系統(tǒng)的參數(shù)估計可表示為：

（10）

嚙合剛度辨識需要傳動誤差數(shù)據(jù)及其第一、第二階時間導(dǎo)數(shù)（）。若僅傳動誤差可測，則可通過對序列的差分計算近似得到。然而，測量過程的噪聲信號將在差分運算中被大幅放大，從而降低嚙合剛度的辨識精度。

為了降低噪聲干擾，本文提出一種基于再生核希爾伯特空間（Reproducing Kernel Hilbert Space，RKHS）［19?20］的降噪方法，用于處理含有噪聲的質(zhì)量項及阻尼耦合項（）：

（11）

式中 為RKHS正則項；γ為權(quán)重系數(shù)；Z由再生核希爾伯特空間映射得到。

式（11）所示優(yōu)化問題的解可表示為：

（12）

式中 為Gram矩陣，向量為質(zhì)量及阻尼耦合項（）在再生核希爾伯特空間的映射：

（13）

式中 單位矩陣的維度為。得益于高斯核的通用性和簡潔性，本文選取高斯核來構(gòu)建矩陣G：

（14）

式中 σ為高斯核函數(shù)的標(biāo)準(zhǔn)偏差。

然后，基于RKHS的降噪方法被進一步應(yīng)用于估計系統(tǒng)狀態(tài)參數(shù)：

（15）

其中：

（16）

在式（16）對應(yīng)的目標(biāo)函數(shù)中，第一項J1確保了系統(tǒng)狀態(tài)的平滑性，第二項J2限制了傳動誤差估計值與其測量值間的偏離；，為權(quán)重系數(shù)。

然后，回歸系數(shù)的估計值更新為：

（17）

可按式（10）更新嚙合剛度的估計值。

將式（13）中的數(shù)據(jù)替換為優(yōu)化后的進行估計，再通過式（12）更新系統(tǒng)質(zhì)量與阻尼項的組合估計。進一步地，通過式（15）和（16）更新系統(tǒng)的狀態(tài)估計。最后，通過式（10）和（17）相應(yīng)地更新嚙合剛度估計。交替更新系統(tǒng)的狀態(tài)與參數(shù)估計直至收斂或達到更新次數(shù)上限。

2.2 時變嚙合剛度辨識流程

齒輪系統(tǒng)狀態(tài)和模型參數(shù)具體迭代更新流程如圖2所示。

步驟1：測量并記錄齒輪副在兩次不同輸入扭矩下的轉(zhuǎn)角qIkb9Fl8gma2tat0HB1qwg==及扭矩信號。通過式（9）初始化回歸系數(shù)，通過式（10）確定嚙合剛度的初始估計。

步驟2：交替優(yōu)化系統(tǒng)狀態(tài)及參數(shù)。

步驟2.1：通過式（12）和（13）估計質(zhì)量及阻尼耦合項Z；

步驟2.2：通過式（15）和（16）估計傳動誤差x；

步驟2.3：通過式（10）和（17）重新估計嚙合剛度k；

步驟2.4：用更新后的傳動誤差及嚙合剛度估計構(gòu)造替換式（13）中的。

步驟3：重復(fù)步驟2直至Z，k，x估計值的變化量小于設(shè)定值或達到迭代次數(shù)上限。

3 仿真與實驗驗證

3.1 仿真算例

為驗證所提齒輪時變嚙合剛度辨識方法的有效性，采用表1中的參數(shù)構(gòu)建齒輪系統(tǒng)動力學(xué)模型（式（1））進行動態(tài)仿真。齒輪副嚙合剛度的真值由MA等［21］所提解析法計算得到，嚙合阻尼依據(jù)瑞利阻尼模型按10%阻尼比取值。幾何傳動誤差參考DONG等［22］采用的建模方式，假設(shè)其由主、從動輪的跳動誤差（偏心）引起，誤差參數(shù)如表1所示。

依據(jù)所提嚙合剛度辨識方法，齒輪動力系統(tǒng)的輸入扭矩設(shè)為10 N·m和30 N·m，輸入軸轉(zhuǎn)速為600 r/min。齒輪系統(tǒng)輸出的傳動誤差信號由四階龍格?庫塔法計算得到。采樣頻率為23 kHz，采樣時間為10 s。為模擬實際測量環(huán)境，輸出的傳動誤差信號按信噪比（Signal?to?Noise Ratio，SNR） =30，25，20 dB混入高斯白噪聲。在回歸矩陣中，傅里葉級數(shù)的基頻（齒輪副嚙合頻率）為230 Hz，頻率最高倍數(shù)N設(shè)為40。同時，非線性動力學(xué)稀疏識別（Sparse Identification of Nonlinear Dynamics，SINDy）方法［23］被用于辨識同一齒輪傳動系統(tǒng)的嚙合剛度。

辨識得到的時變嚙合剛度曲線及其真值如圖3所示。表2為所提方法與SINDy方法在不同信噪比下辨識結(jié)果的均方根誤差（Root Mean Square Error，RMSE）指標(biāo)，其定義為剛度估計與其真值間的均方根誤差。

識別齒輪系統(tǒng)參數(shù)需要同時用到DTE信號及其一階、二階導(dǎo)數(shù)數(shù)據(jù)。當(dāng)僅DTE信號可測時，SINDy方法直接使用對噪聲敏感的傳動誤差信號高階導(dǎo)數(shù)進行剛度估計。與SINDy方法相比，本文方法通過狀態(tài)信號在再生核希爾伯特空間的映射變換有效削弱噪聲影響。此外，本文方法引入了齒輪動力系統(tǒng)的物理模型與測量數(shù)據(jù)雙重約束，在保證所估計剛度曲線和狀態(tài)曲線平滑性的同時，限制了狀態(tài)估計與其測量值間的偏差。本節(jié)嚙合剛度的辨識結(jié)果表明，在相同噪聲環(huán)境下，所提方法相較于SINDy方法具有更高的辨識精度。

3.2 實驗驗證

為進一步驗證所提嚙合剛度辨識方法的有效性，以圖4所示的齒輪箱為對象開展實驗。安裝在輸入軸和輸出軸端的編碼器每轉(zhuǎn)產(chǎn)生1024個脈沖信號。為了簡化實驗配置，嚙合力近似為，其中輸入扭矩可由安裝在輸入軸側(cè)的單個動態(tài)扭矩傳感器檢測得到。在輸入扭矩大小分別為7.6 N·m和1.2 N·m的條件下開展測試，輸入軸轉(zhuǎn)速為600 r/min。數(shù)據(jù)采集儀的采樣頻率為500 kHz，采樣時間為5 s。圖5為兩次不同輸入扭矩下動態(tài)傳動誤差信號的測量結(jié)果。對兩次不同嚙合力對應(yīng)的實測DTE信號按相關(guān)性最大化進行對齊，然后計算其差值，得到系統(tǒng)的輸出信號x。

齒輪箱參數(shù)與仿真中的齒輪副參數(shù)（見表1）一致，區(qū)別在于齒輪實驗臺的輸入、輸出軸是柔性件（因此會參與構(gòu)成齒輪副的嚙合）。相應(yīng)地，齒輪嚙合剛度kGM（t）需要從總剛度估計中按：

（18）

解耦得到，其中為柔性軸扭轉(zhuǎn)剛度。

此外，傳動誤差信號的直流分量無法通過編碼器測量得到。為了確定未知直流分量a，系數(shù)向量的初始化過程被改寫為：

（19）

其中：

（20）

傳動誤差信號x被修正為用以估計系統(tǒng)狀態(tài)及參數(shù)。圖6所示為采用所提方法和SINDy方法辨識得到的齒輪嚙合剛度曲線。作為參考，依照表1所示齒輪箱的先驗信息，利用MA等［20］所提解析法計算嚙合剛度的理論值。實驗結(jié)果顯示，所提方法辨識的嚙合剛度與解析法計算結(jié)果的一致性明顯優(yōu)于SINDy方法，表明所提齒輪時變嚙合剛度方法具有較高的辨識精度。

4 結(jié) 論

本文針對齒輪傳動系統(tǒng)時變嚙合剛度建模與辨識問題，提出了一種狀態(tài)/參數(shù)交替優(yōu)化估計方法。首先，為表征嚙合剛度時變特性，構(gòu)建以齒輪嚙合頻率為基頻的傅里葉基函數(shù)庫來擬合時變剛度，并提出了兼顧齒輪動力系統(tǒng)物理模型約束與剛度曲線平滑性要求的參數(shù)初始化方法。進一步地，通過引入物理模型與數(shù)據(jù)雙重約束，實現(xiàn)對系統(tǒng)狀態(tài)與模型參數(shù)的交替優(yōu)化。同時，狀態(tài)測量信號的噪聲在再生核希爾伯特空間的映射變換中被削減，有效提高了嚙合剛度辨識精度。齒輪傳動的數(shù)值仿真和動態(tài)實驗驗證了所提齒輪時變嚙合剛度辨識方法的有效性。需要注意的是，當(dāng)系統(tǒng)轉(zhuǎn)速較高時，實驗測試的采樣頻率不宜設(shè)置得過低，否則會影響對輸入、輸出軸轉(zhuǎn)角的測量精度，從而降低嚙合剛度的辨識精度。

參考文獻：

［1］馬志賽，丁千，劉莉，等. 線性時變結(jié)構(gòu)模態(tài)參數(shù)時域辨識方法的研究進展［J］. 機械工程學(xué)報，2018，54（23）： 137-159.

MA Zhisai，DING Qian，LIU Li，et al. Research progress on time-domain modal parameter estimation methods for linear time-varying structures［J］. Journal of Mechanical Engineering，2018，54（23）： 137-159.

［2］蔡宇，劉旭，程英豪. 線性變參數(shù)振動系統(tǒng)的全局辨識［J］. 機械工程學(xué)報，2022，58（23）： 114-122.

CAI Yu，LIU Xu，CHENG Yinghao. Global identification of linear parameter-varying vibration systems［J］. Journal of Mechanical Engineering，2022，58（23）： 114-122.

［3］周國峰，劉博. 某空氣舵系統(tǒng)非線性建模及參數(shù)識別［J］. 振動工程學(xué)報，2023，36（5）： 1216-1222.

ZHOU Guofeng，LIU Bo. Nonlinear modeling and parameter identification of an air rudder system［J］. Journal of Vibration Engineering，2023，36（5）： 1216-1222.

［4］SHI J F，GOU X F，ZHU L Y. Five-state engaging model and dynamics of gear-rotor-bearing system based on time-varying contact analysis considering gear temperature and lubrication［J］. Applied Mathematical Modelling，2022，112： 47-77.

［5］NATALI C，BATTARRA M，DALPIAZ G，et al. A critical review on FE-based methods for mesh stiffness estimation in spur gears［J］. Mechanism and Machine Theory，2021，161： 104319.

［6］CHEN Z G，ZHOU Z W，ZHAI W M，et al. Improved analytical calculation model of spur gear mesh excitations with tooth profile deviations［J］. Mechanism and Machine Theory，2020，149： 103838.

［7］COOLEY C G，LIU C G，DAI X，et al. Gear tooth mesh stiffness： a comparison of calculation approaches［J］. Mechanism and Machine Theory，2016，105： 540-553.

［8］YANG L T，ZENG Q，YANG H S，et al. Dynamic characteristic analysis of spur gear system considering tooth contact state caused by shaft misalignment［J］. Nonlinear Dynamics，2022，109（3）： 1591-1615.

［9］PANDYA Y，PAREY A. Experimental investigation of spur gear tooth mesh stiffness in the presence of crack using photoelasticity technique［J］. Engineering Failure Analysis，2013，34： 488-500.

［10］RAGHUWANSHI N K，PAREY A. Experimental measurement of gear mesh stiffness of cracked spur gear by strain gauge technique［J］. Measurement，2016，86： 266-275.

［11］RAGHUWANSHI N K，PAREY A. A new technique of gear mesh stiffness measurement using experimental modal analysis［J］. Journal of Vibration and Acoustics，2019，141（2）： 021018.

［12］DAI H，LONG X H，CHEN F，et al. An improved analytical model for gear mesh stiffness calculation［J］. Mechanism and Machine Theory，2021，159： 104262.

［13］SHI J L，MA X G，XU C L，et al. Meshing stiffness analysis of gear using the Ishikawa method［J］. Applied Mechanics and Materials，2013，401-403： 203-206.

［14］HU Z H，TANG J Y，ZHONG J，et al. Effects of tooth profile modification on dynamic responses of a high speed gear-rotor-bearing system［J］. Mechanical Systems and Signal Processing，2016，76-77： 294-318.

［15］CHANG L H，LIU G，WU L Y. A robust model for determining the mesh stiffness of cylindrical gears［J］. Mechanism and Machine Theory，2015，87： 93-114.

［16］CHEN K K，HUANGFU Y F，MA H，et al. Calculation of mesh stiffness of spur gears considering complex foundation types and crack propagation paths［J］. Mechanical Systems and Signal Processing，2019，130： 273-292.

［17］ARLOT S，CELISSE A. A survey of cross-validation procedures for model selection［J］. Statistics Surveys，2010，4： 40-79.

［18］KINGMA D P，BA J. ADAM： a method for stochastic optimization［J］. arXiv preprint arXiv：1412.6980，2014.

［19］LIN M M，CHENG C M，PENG Z K，et al. Nonlinear dynamical system identification using the sparse regression and separable least squares methods［J］. Journal of Sound and Vibration，2021，505： 116141.

［20］CHENG C M，ZHAO B X，F(xiàn)U C K，et al. A two-stage sparse algorithm for localization and characterization of local nonlinear structures［J］. Journal of Sound and Vibration，2022，526： 116823.

［21］MA H，ZENG J，F(xiàn)ENG R J，et al. An improved analytical method for mesh stiffness calculation of spur gears with tip relief［J］. Mechanism and Machine Theory，2016，98： 64-80.

［22］DONG X J，HUANGFU Y F，YU X L，et al. Identification of the error excitation in gear systems： a mediator algorithm between simulation and experiment［J］. Journal of Sound and Vibration，2024，568： 118060.

［23］BRUNTON S L，PROCTOR J L，KUTZ J N. Discovering governing equations from data by sparse identification of nonlinear dynamical systems［J］. Proceedings of the National Academy of Sciences，2016，113（15）： 3932-3937.

Time-varying mesh stiffness identification of gear systems based on an alternating state/parameter optimization estimation

LIN Miao?miao1，CHENG Chang-ming1，PENG Zhi?ke1，2

（1.State Key Laboratory of Mechanical System and Vibration，Shanghai Jiao Tong University，Shanghai 200240，China； 2.School of Mechanical Engineering，Ningxia University，Yinchuan 750021，China）

Abstract： The time-varying mesh stiffness is a core parameter of gear systems，and the mesh stiffness identification is of great significance for the dynamic analysis and condition monitoring of gear transmission systems. Since it is difficult to directly measure the mesh stiffness，it is necessary to develop a data-driven time-varying mesh stiffness identification method. To deal with this problem，an alternating state-parameter optimization method is proposed to identify the time-varying mesh stiffness of gear systems. The Fourier series with the fundamental frequencipayzdxpwuQ7Cwg3lQ1RjA==y of the mesh frequency is constructed to characterize the mesh stiffness. Furthermore，a Reproducing Kernel Hilbert Space （RKHS）-based de-noise method is further proposed to estimate the system state and parameter. The system state and stiffness parameter are alternately optimized with the joint constrains of dynamic model and data to realize the time-varying mesh stiffness identification of gear transmission systems. The simulation and experimental studies validate the effectiveness of the new mesh stiffness identification method for gear systems.

Key words： parameter identification；gear transmission system；time-varying mesh stiffness；alternating optimization

作者簡介： 林苗苗（1997―），女，博士研究生。E-mail：linmmiao@sjtu.edu.cn。

通訊作者： 程長明（1987―），男，博士，副教授，博士生導(dǎo)師。電話：（021）34206831-828；E-mail：ccming@sjtu.edu.cn。