基于快速稀疏正則化的旋轉梁損傷識別*

楊達豪，呂中榮，汪利

中山大學航空航天學院，廣東廣州 510006

在工程中，風機葉片［1］、直升機葉片［2］等是相關結構的重要部件，并且可以簡化成旋轉梁模型。在長期的工作環境下，這些重要部件容易產生裂紋、發生損傷、甚至斷裂，嚴重地影響使用安全和經濟效益。為此，有必要對旋轉梁實行結構健康監測，維護整體結構正常運作。已有的研究表明，旋轉梁會發生大幅運動和微彈性變形［3］；在科氏力和離心力的作用下，容易產生剛度效應［4-5］。為此，引入非笛卡爾坐標系變量延伸變形s（stretch deformation）與笛卡爾坐標系變量翼弦變形v（chordwise deformation）、揮舞變形w（flapwise deformation），以建立可以反映旋轉梁更為真實運動的線性偏微分方程，并且使用有限單元法進一步得到離散方程［6］。

結構損傷識別是典型的反問題，在噪聲影響或者數據不足的情況下，該反問題往往是一個病態的待求解目標方程。為了克服該弊端，需要引入正則化約束，如Tikhonov 正則化［7-8］，稀疏正則化［9-11］。Tikhonov 正則化是L2 范數并具有連續可導性質。目標方程在引入Tikhonov 正則化之后，可以通過SVD 分解，獲得解析解，但正則化參數仍然由L-curve 曲線確定［7］。實際上，結構損傷數目是具有稀疏性的，即損傷數目應該是少量的，這與稀疏正則化所呈現的性質一致。為此，本文使用稀疏正則化求解目標函數，以呈現更好的魯棒性。稀疏正則化項即L1 正則化項，由于不具有連續可導性質，使得目標方程求解顯得較為困難。目前已有較為成熟的算法［10］，如LASSO 法（least absolute selection and shrinkage operator）、內點法（interior point method）、主動集法（active set method）。但這些方法都需要確定一個合適的稀疏正則化參數。已有的稀疏正則化參數確定方法有［11］：L-curve曲線法，差異原則法（method of discrepancy principle）以及廣義交叉驗證法（GCV，generalized cross-validation method）等。這些方法應用較為廣泛，但容易受到噪聲的影響且不可避免地有著計算代價高的缺點。為克服此類缺點，本文提出一種受摩擦模型啟發的快速稀疏正則化方法。該方法指出，引入稀疏正則化相當于引入摩擦效應，其中稀疏正則化參數對應于靜摩擦力。基于摩擦模型的物理意義，可直接選定某個摩擦力為正則化參數，從而快速求解目標函數，識別損傷程度和位置。數值算例表明，該方法具有效率高、魯棒性強的優勢。

1 旋轉梁動力模型

如圖1所示，旋轉梁的長度為L，質量密度為ρ，楊氏模量為E，截面面積為A，關于y軸和z軸的極慣性矩為Iy，Iz，且有Iy=Iz。給定旋轉速度Ω，旋轉梁上一點從P運動到P'，其位移向量設為dr。在笛卡爾坐標系下，該點的運動速度（證明見附錄）為［12］：

圖1 旋轉梁模型Fig. 1 Configuration of a rotating cantilever beam

其中dw、dsv分別是揮舞運動變形向量，翼弦變形向量；Sw、Ssv是由旋轉運動所產生的剛度矩陣，Gsv是陀螺矩陣（gyroscopic matrix）；Msv、Mw是質量矩陣；Kw、Ksv是剛度矩陣。顯然，方程（4）、（5）分別表征旋轉梁的揮舞運動方程、翼弦運動方程。考慮勻速旋轉的情況，即Ω? = 0，通過求解方程（4）、（5）的特征值問題［13］，可以獲得旋轉梁的模態數據，包括固有頻率λ（特征值）和模態φ（特征向量）。

2 目標函數

結構損傷通常表現為系統剛度的局部減小，即相關單元剛度發生折減。設損傷系數為pi，pi?(-1，0]，其中i表示的是第i個單元。pi= -1時表示第i個單元斷裂，即完全損傷；pi= 0時表示第i個單元完好。一般，損傷位置是稀疏的，即系統的損傷參數p=[p1，p2，p3，...，pn]具有少量非零元素，即具有稀疏性。

系統剛度可以寫成關于損傷參數的函數，即

取最小值，其中W是一個權重矩陣。

方程（9）是一個非線性目標函數，可以通過迭代法進行求解，即：

①給定初始值p0=[0，0，0，...，0]。

②基于第（k- 1）步迭代結果，求解第k步的迭代值，即pk=pk-1+ Δp，其中Δp是在第k步求解得到的增量，k=[1，2，3，...，n]。

③當滿足收斂條件的時候，結束迭代求解。

為了求得增量Δp或參數pk，在第k- 1步時，對R(pk)在pk-1處進行泰勒級數展開，即線性化處理

其中γ是正則化參數。通過求解方程（13）的最小值，即可求解得到增量Δp或參數pk.

3 快速稀疏正則化

為求解最終目標方程（13），首先要確定稀疏正則化參數γ。為此，可從一個由質量塊、庫倫摩擦單元、彈簧單元組成的多自由度摩擦模型出發，如圖2 所示?；谠撃Ｐ?，可推導出正則化參數γ的表達式，即快速稀疏正則化。

圖2 快速稀疏正則化的物理摩擦模型Fig.2 Illustration of fast sparse regularization as a friction model

設整個物理摩擦模型系統所作的功為fγ(α)，即其中α、A、b、γ分別是質量塊的位移向量、系統的剛度矩陣、外荷載向量、庫倫摩擦單元的靜摩擦力。方程（14）第一項表示彈簧存儲的能量，第二項為外力做的功。方程（15）的第二項為摩擦消耗的能量，也就是說稀疏正則項的物理意義為能量耗散或外力做功［14］。

通過以上“摩擦力”的選擇方程，即方程（20），可以有效控制解的稀疏性。方程（21）中的La在初始迭代時具有較大值，即正則化參數γ具有較小值，即“更多的滑塊可以滑動”，這可以避免損傷程度和位置信息的丟失。在迭代一定次數之后，La具有最小值，此時γ取得最大值，即“只有個別的滑塊可以滑動”，從而保證了最終結果的稀疏性。

4 數值算例

考慮一旋轉梁，其長度L= 0.5 m，質量密度ρ= 2 700 kg/m3，楊氏模量E= 71GPa，截面面積A=10-4m2，關于y軸和z軸的極慣性矩Iy=Iz= 8.33× e-12m4。轉臺半徑a= 0.1m。將旋轉梁離散成11 個單元，設第3個單元發生損傷，損傷稀疏為pi= -0.3?？紤]轉速Ω= 30 rad/s?？紤]特征值噪聲為0.25%，特征向量噪聲為5%。對于方程（21） 中La的選擇，令d= 1，Lmin= 1，Lmax= 11。在無噪聲情況下，N= 10；有噪聲情況下，N= 100。基于快速稀疏正則化識別方法，只取前四階模態數據即可取得較好的識別結果，基于揮舞模態數據和翼弦模態數據的損傷識別結果分別如圖3-4 所示。為了進一步突出本文所提方法的有效性、快速性，選取L2正則化方法［7］進行識別結果對比。為了取得較好的識別結果，需要使用前六階的模態數據?；趽]舞模態數據和翼弦模態數據的損傷識別結果分別如圖5-6所示。

圖3 基于揮舞模態數據的識別結果Fig. 3 The identification results based on the flapwise modal data

圖4 基于翼弦模態數據的識別結果Fig. 4 The identification results based on the chordwise modal data

圖5 基于揮舞模態數據的識別結果（使用L2正則化）Fig. 5 The identification results based on the flapwise modal data(with L2 method)

圖3- 4 表明，在沒有噪聲的情況下，基于兩種模態數據，快速稀疏正則化方法在第7 次迭代之后，識別結果已經接近于準確值。在有噪聲的情況下，基于揮舞模態數據的識別結果要優于基于翼弦模態數據的識別結果，這是因為翼弦運動方程含有延伸變形，同時含有廣義阻尼項2ΩGsv，使得翼弦運動方程更加復雜，容易受到轉速的影響，求解相應的模態靈敏度也更加困難。在受到噪聲的影響情況下，基于快速稀疏正則化的迭代收斂次數也僅20次。相比L-curve曲線方法等選擇不同的正則化參數，快速稀疏正則化方法可以基于摩擦模型，直接選取“摩擦力”為正則化參數，在效率上具有顯著優勢［11］。此外，如圖5-6 所示，L2 正則化在無噪聲情況下，可以取得較為準確的識別結果，但是所需要的模態數據和迭代求解次數顯著增加，這對數據采集的傳感器、計算機性能提出了更高的要求。在噪聲影響下，L2 正則化方法錯誤地識別若干單元具有損傷，并且p3的識別精度顯著下降。由于翼弦模態響應更加復雜，L2 正則化識別方法在識別過程中也更加容易受到噪聲數據的影響，如圖6所示。

圖6 基于翼弦模態數據的識別結果（使用L2正則化）Fig. 6 The identification results based on the chordwise modal data(with L2 method)

為了確定快速稀疏正則化方法的識別精度，各工況的識別結果列于表1。基于快速稀疏正則化的識別結果，最大的相對誤差發生在基于翼弦模態數據（噪聲）的第3個單元，僅為6.57%，表明所提方法具有較高的損傷識別精度和魯棒性?；贚2 正則化的識別結果，該方法容易錯誤地識別單元損傷，最大誤差發生在基于翼弦模態數據（噪聲）（L2）的第六個單元，錯誤識別值為-0.053 6；最大相對誤差發生在基于揮舞模態數據（噪聲）（L2）的第3 個單元，為16.4%。綜上，基于快速稀疏正則化方法不僅可以利用損傷位置和程度的稀疏性特性，提升識別結果的精度；還能準確地選擇合適的正則化參數，提升收斂速度。

表1 基于快速稀疏正則化的旋轉梁損傷識別結果Table 1 Identification result of rotating beam with the fast sparse regularization approach

5 小 結

本文提出一種快速稀疏正則化方法識別旋轉梁的損傷。該方法從摩擦模型出發，指出正則化參數可以通過“摩擦力”大小直接確定，并能控制解的稀疏性。數值算例表明，通過對比L2 正則化方法和L-curve 方法，所提方法不僅計算效率高且可以準確快速地識別出旋轉梁的損傷程度和位置。同時，所提的快速稀疏正則化方法表明，基于揮舞模態數據的識別結果要優于基于翼弦模態數據的識別結果，這是因為翼弦運動方程更為復雜，容易受到旋轉速度的影響，且模態的靈敏度分析更為困難。

附 錄

公式（1）是建立在旋轉坐標系下。設慣性坐標系為oij，由于軸保持不變，平面發生旋轉可表示成圖7。圖1中，P點運動到P'，P'的坐標rp為

圖7 旋轉梁平面示意圖Fig.7 Configuration of a rotating flexible cantilever beam in-plane deformation

其中α是轉過的夾角，既有α?=Ω，對公式（24）關于時間t求導，可得

將公式（25）帶入公式（23），可得公式（1）。