基于狀態/參數交替優化估計的齒輪時變嚙合剛度辨識方法研究

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

關鍵詞： 參數辨識； 齒輪傳動系統； 時變嚙合剛度； 交替優化

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

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

引 言

齒輪結構被廣泛應用于機械傳動系統。在齒輪動力學模型中，時變嚙合剛度是反映齒輪系統動態特性的重要參數，其數值變化可導致機械傳動系統產生振動和噪聲。對齒輪副嚙合剛度進行建模與識別，有助于齒輪傳動系統的優化設計和狀態監測［1?5］。

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

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

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

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

（1）

其中：

（2）

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

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

（3）

（4）

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

（5）

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

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

2 齒輪系統狀態及參數交替優化

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

2.1 時變嚙合剛度辨識方法

為辨識齒輪傳動系統時變嚙合剛度，以嚙合頻率為基頻構建傅里葉基函數庫，用來表征齒輪副的時變嚙合剛度與嚙合阻尼：

（6）

假設有M個采樣點，為辨識系統參數，根據方程（5），可采用最小二乘法估計回歸系數：

（7）

其中：

（8）

然而，最小二乘法容易產生過擬合，為了克服這一缺陷，引入正則項優化參數估計：

（9）

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

相應地，齒輪系統的參數估計可表示為：

（10）

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

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

（11）

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

式（11）所示優化問題的解可表示為：

（12）

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

（13）

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

（14）

式中 σ為高斯核函數的標準偏差。

然后，基于RKHS的降噪方法被進一步應用于估計系統狀態參數：

（15）

其中：

（16）

在式（16）對應的目標函數中，第一項J1確保了系統狀態的平滑性，第二項J2限制了傳動誤差估計值與其測量值間的偏離；，為權重系數。

然后，回歸系數的估計值更新為：

（17）

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

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

2.2 時變嚙合剛度辨識流程

齒輪系統狀態和模型參數具體迭代更新流程如圖2所示。

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

步驟2：交替優化系統狀態及參數。

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

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

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

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

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

3 仿真與實驗驗證

3.1 仿真算例

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

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

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

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

3.2 實驗驗證

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

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

（18）

解耦得到，其中為柔性軸扭轉剛度。

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

（19）

其中：

（20）

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

4 結 論

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

參考文獻：

［1］馬志賽，丁千，劉莉，等. 線性時變結構模態參數時域辨識方法的研究進展［J］. 機械工程學報，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］蔡宇，劉旭，程英豪. 線性變參數振動系統的全局辨識［J］. 機械工程學報，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］周國峰，劉博. 某空氣舵系統非線性建模及參數識別［J］. 振動工程學報，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，FU 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，FENG 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―），男，博士，副教授，博士生導師。電話：（021）34206831-828；E-mail：ccming@sjtu.edu.cn。