基于參數降階模態的分層貝葉斯在線裂紋檢測

魏保立 ，馮雅珊 ，羅坤，付偉

(1.鄭州鐵路職業技術學院機電工程學院，河南鄭州 451460；2.河南理工大學能源科學與工程學院，河南焦作 454003)

0 前言

在基于狀態的維修領域，面向模型的方法已成為解決退化和損傷識別問題的一種常用方法，不僅在檢測和定位層面產生了重要作用，而且在更精細的量化階段也能夠發揮其優勢[1]。與強烈依賴密集傳感信息的純數據驅動方案不同，系統模型的可用性存在一定程度的冗余，允許提取監測量和模型預測量之間的相互依賴關系[2]。如何設計一種精確的裂紋或損傷檢測方法成為健康管理領域的關鍵。

目前，研究者已經在傳統的裂紋或損傷檢測的基礎上提出了許多具有一定功能的替代方法。在各種方法中，遞歸貝葉斯估計器是適應實時性能要求的可靠手段，并且已經作為成熟的模型應用在航空航天工程、電信、金融和控制工程等領域[3-4]。在遞歸貝葉斯框架內，損傷識別本質上被歸類為一個狀態參數估計問題，這是一個非線性問題，即使對于實際動力學處于線性狀態的系統也是如此。文獻[5]引入貝葉斯分類方法，優化選擇裂紋形核的信息狀態變量預測因子，建立了一個簡單的標量裂紋形核指示器。文獻[6]針對當前疲勞裂紋擴展預測研究中較少考慮不確定因素而導致預測結果偏差大的問題，提出基于動態貝葉斯網絡的疲勞裂紋擴展預測方法，以變幅載荷作用下的疲勞裂紋擴展為研究對象。利用統一疲勞壽命預測模型構建疲勞裂紋擴展的物理狀態方程。雖然上述方法取得了一定效果，但是與這些濾波器相關的數值問題之一是樣本貧化。

為此，研究者提出了一系列基于狀態增強的參數估計方法。文獻[7]提出了一種基于自然梯度法的變分貝葉斯卡爾曼濾波的變分下界最大化技術，得到了狀態估計的變分超參數和相應的誤差協方差的估計。文獻[8]提出了一種新的綜合預測模型貝葉斯優化隨機森林-卡爾曼濾波，其中卡爾曼濾波和隨機森林算法分別用于預測裂紋的趨勢和周期。上述方法中的一個關鍵挑戰是過程和測量噪聲項的調整，因為在這些問題中，推理過程是基于初始參數值和互協方差項實現的，從本質上調節了觀測狀態的信息轉化為未觀測參數的修正過程。當參數本身表現出很強的相互關聯時，問題變得更具挑戰性。此外，過程噪聲協方差交叉項的靜態值可能不足以解決狀態參數問題，其中系統的動態特性取決于參數估計，因此需要進行調整，以避免不穩定性。

為解決上述方法中存在的缺點，本文作者提出一種基于參數降階模態的分層貝葉斯在線裂紋檢測方法。在實時性能的約束下，系統動力學由一個參數化的降階模型表示，該模型與所尋求的特征有關，隨后與分層貝葉斯方法融合，解決參數空間采樣后生成的多個模型配置解的輸入狀態估計問題。最后通過仿真驗證所提方法的有效性。

1 問題描述

1.1 參數化降階模型

文中的研究目的是僅輸出振動測量值，在線檢測結構或機械系統在運行條件下的裂紋，并進行輸入狀態估計?；旧峡赏ㄟ^假設實際振動響應測量的可用性來實現，這些測量是從高保真(High Fidelity，HF)有限元模型和所考慮的系統中綜合獲得的。為了適應不同的裂紋結構，后者根據裂紋特征即位置、長度和方向進一步參數化。因此，連續時間系統動力學由以下參數相關的二階微分方程表示：

(1)

其中：u∈Rn是位移向量，n表示有限元模型空間離散化產生的自由度數；M(θ)、C(θ)、K(θ)∈Rn×n，分別是質量矩陣、阻尼矩陣和剛度矩陣；θ∈Rm是參數向量，此處表示裂紋的幾何特征。等式(1)的右側由力矢量組成p(t)∈Rnp和相應的選擇矩陣Sp∈Rn×np，它指定p(t)選擇相應的自由度。應注意的是，上述方程通常是有限元離散化的結果，其應足夠密集，以詳細地捕捉系統的響應。因此，在在線設置中使用這些方程的積分可能會消耗大量資源。

為了使系統動力學預測的計算效率更高，從而允許在線解決輸入、狀態和參數估計問題，將式(1)的解投影到一個子空間上，該子空間的維數為r，比全訂單模型小得多(r?n)，由于系統動力學依賴于參數向量θ，減少的基礎也是θ，這意味著投影步長可以表示為

u(t)=V(θ)q(t)

(2)

其中：V(θ)∈Rn×r包含約化空間的基向量；q(t)∈Rr是對應廣義坐標的向量。將式(2)代入式(1)并結合uT得到降階運動方程：

V(θ)TSpp(t)

(3)

通過本征正交分解(POD)可以獲得全局約化基，從而求解多個參數樣本的運動方程θj(j=1,2,…,ns)，ns表示樣本的數量。然后，通過收集所有參數樣本位移向量的迭代次數構建Snapshot矩陣，如式(4)所示：

S=[U(θ1)U(θ2) …U(θns)]

(4)

其中：U(θj)∈Rn×nt指定從等式(1)的解中獲得的Snapshot矩陣θ=θj，nt表示解決方案步驟的數量，即Snapshot的數量。最后，可以通過最初提取的Snapshot矩陣的奇異值分解(Singular Value Decomposition，SVD)來執行縮減步驟，然后僅保留第r個顯著性來截斷相應的基向量，這些基向量用于形成全局約簡基V∈Rn×r。對于具有大量時間步長的動態問題，可以增加一個額外的步驟來減少Snapshot矩陣的存儲需求，該步驟涉及在Snapshot級別通過SVD進行縮減。

由于參數相關問題的全局POD方法通常會導致一個超大的基礎V，因此，為了消除計算增益，文中采用了一種基于聚類的POD方案，其中包含局部基和網格變形，作為表示不同裂紋形態的一種手段。根據前述的過程，對參數空間的某些區域提取約簡基V，這些區域通過k-均值聚類算法進行提取，并使用網格變形對未經訓練的裂紋配置的解進行插值。

1.2 狀態空間表示

為了將建模的系統動力學與實際振動測量值融合，需要將參數化的降階運動方程即公式(3)轉換為連續的時間-狀態-空間表示，其表示為

(5)

y(t)=G(θ)x(t)+J(θ)p(t)

(6)

(7)

(8)

輸出矩陣和饋通矩陣的結構G(θ)和J(θ)分別由觀測到的響應量確定。對于結構系統，這些量通常是應變和加速度信號，因此得出以下表達式：

G(θ)=

(9)

(10)

其中：Sa∈Rna×n是一個布爾矩陣，用于選擇加速度測量和；Sd，ε∈Rnd×n是與應變測量相關的選擇矩陣。在處理有限元模型時，應變是在單元級別計算的，因此Sd，ε選擇相應單元自由度的位移，然后使用矩陣將它轉換為應變De(θ)∈Rnε×nd。后者具有塊對角形式，每個塊條目包含元素的變形矩陣。

2 遞階狀態輸入參數估計

2.1 狀態輸入估計

式(5)(6)在稀疏和噪聲輸出的基礎上，僅對測量值進行離散化。為了適應離散的自然響應測量，將式(5)、(6)轉移到離散時間域，如下所示：

xk+1=Ad(θ)xk+Bd(θ)pk+vk

(11)

yk=Gd(θ)xk+Jd(θ)pk+wk

(12)

同時進一步補充了零均值高斯噪聲項vk∈R2r和wk∈Rny，其協方差矩陣分別為Qv∈R2r×2r和Qw∈Rny×ny。矩陣Ad(θ)、Bd(θ)、Gd(θ)和Jd(θ)是離散時間對應的A(θ)、B(θ)、G(θ)和J(θ)，其零階保持表達式如下所示：Ad(θ)=eA(θ)dt、Bd(θ)=(A(θ)-I)·A-1(θ)B(θ)、Gd(θ)=G(θ)和Jd(θ)=J(θ)。系統參數θ在估計過程中發生變化時，可以使用泰勒級數展開，以犧牲精度為代價，更有效地逼近狀態方程矩陣，從而得到Ad(θ)≈I+A(θ)dt+0.5A2(θ)dt2+O(dt3)及Bd(θ)≈B(θ)dt。

(13)

式(13)使用了過程模型的馬爾可夫特性：p(xk|xk-1,Yk-1)=p(xk|Yk-1)。在計算先驗知識后，可以使用貝葉斯定理進行更新：

(14)

其遞推公式如下：

(15)

用由狀態空間表示的似然函數p(yk|xk)的測量方程定義。所尋求的過濾概率密度函數的估計為p(xk|Yk)，在工程應用中，狀態的統計矩通常是有意義的，可以計算如下：

(16)

(17)

E[·]表示期望算子和上述表達式的積分極限(-∞，+∞)。在線性系統方程組和加性高斯噪聲的特殊情況下，由著名的卡爾曼濾波器獲得式(11)(12)的解。雖然該系統由線性表達式(9)(10)描述，并且假設加性高斯噪聲，即額外估計未知參數的問題θ和輸入pk在狀態xk下變為非線性，因此需要替代解決方案策略。序列蒙特卡羅(Sequential Monte Carlo，SMC)采樣器構成了這一大類方法，也將在文中使用。

盡管狀態增強是解決狀態參數估計問題的一種直接技術，但該方法的實現強烈依賴于初始狀態變量和追求的參數之間的耦合，以及初始狀態和參數值的均值和協方差。由于參數本身不可觀測，因此其估計取決于狀態變量(可能部分可觀測)。這意味著，參數動態由狀態參數協方差的耦合項控制，它本質上調節了狀態更新轉化為未觀測參數的修正。因此，狀態和參數估計的精度與系統的實際誤差統計密切相關，根據定義可知，這些誤差統計是未知的，只能得到近似值。

為了避免使用未知參數增加狀態向量，文中使用了一組濾波器，使它適應所尋求的系統參數。在這種情況下，未知的裂紋幾何特征或通常要估計的系統參數為一組隨機變量，狀態和參數估計問題可分解為

(18)

(19)

圖1 提出算法的結構框圖

式(19)中所描述的分解本質上是Rao-Blackwellization操作，然而，由于某些原因，文中方案與標準Rao-Blackwellization粒子過濾器(RBPF)有所不同。也就是說，所提方法基于輸入狀態、參數更新和校正步驟的不同時間尺度。此外，參數動力學不受虛構模型的控制，在此類問題中通常假設隨機游動。相反，通過進化策略(Evolutionary Strategies，ES)探索參數空間，即協方差矩陣自適應進化策略(Covariance Matrix Adaptation Evolutionary Strategies，CMA-ES)中使用的策略[9]。因此，收斂參數估計的關鍵要素由CMA-ES的步長保證，該步長強制參數核收縮。

(20)

聯合輸入狀態估計可以使用AKF遞歸獲得，AKF本質上是一個標準的卡爾曼濾波器(KF)，在增廣狀態空間模型上操作。

(21)

(22)

隨后通過如下遺忘因子更新測量噪聲的實際協方差矩陣：

(23)

其中：δQw表示用戶定義的遺忘因子，用于調節預期更新的級別。δQw會減弱括號中的術語，得到非自適應方案，而如果δQw的值較小，則更新是突然進行的。文中所采用的協方差更新方案優于更為成熟的噪聲識別方法，例如自協方差最小二乘(Autocovarianoe Least Squares，ALS)方法，由于可以通過遞歸實現該方法，因此能夠在線檢測。

通過試錯法調整遺忘因子δQw，直到濾波器產生次優或接近最優的結果。文中采用了同樣的啟發式方法，這證明了該方法可以提供足夠準確的結果，這可以歸因于這樣一個事實，即整個方法的關鍵因素不是每個輸入狀態估計器本身的最優性，而是濾波器組的公平權重，這也由次優自適應方案保證。然而，對于全自動校準δQw，更新策略可以與優化方案相結合，同時可以在LOO設置中進一步實施，以確保對創新序列的無偏估計，從而確保自適應的最佳性。

2.2 參數估計

(24)

當處于當前步驟k，可以假設yk-1滿足p(yk-1|θj)=p(yk-1)=1，同樣適用之前的所有觀測yk-1，那么p(Yk-1|θj)=p(Yk-1)=1。因此p(θj|yk-1)=p(θj|Yk-1)=p(θj)，根據貝葉斯規則，在式(24)中進行替換后，得出以下遞推公式：

(25)

參數向量的離散狀態θ應達到參數空間要求的分辨率，在離線應用中，可以通過參數空間密集離散化來解決，原因為計算時間不是一個關鍵因素。然而，在需要實時性能且只允許有限數量粒子的在線應用中，情況并非如此。

(26)

2.3 更新策略

參數的更新不是由隨機游走模型決定的，通常是在增廣貝葉斯方案中假設的。相反，參數通過CMA-ES進行更新，超過簡并閾值時調用CMA-ES。這意味著參數更新步驟是在與執行輸入和狀態校正的時間尺度不同的時間尺度下執行的。從這個層面上說，每個重采樣事件都會創建新的參數樣本，這些樣本來自多元正態分布，因此：

(27)

(28)

(29)

算法1：用于參數和輸入狀態估計僅輸出的分層濾波器。

(5)fork=1，2，…，Tdo

(6)forj=1，2，…，nθdo

(12)end for

(16)forj=1，2，…，nθdo

(20)end for

(21)end for

(30)

(31)

(32)

其中：

(33)

c1和cμ分別是更新協方差矩陣的排名第一和第nθ，p的學習率。算法1中記錄了實現分層濾波器的詳細步驟。

3 仿真案例

3.1 機身面板

通過2個仿真示例驗證所提出的方法，利用稀疏傳感信息識別裂紋特征。為了避免反演問題，每個示例的合成振動測量值的生成都基于高保真模型(HFM)，這2種應用中的高保真模型與用于反向問題的模型不同。即，HFM使用高度精細的網格構建，該網格獨立于ROM使用的網格，而其材料特性的最小空間變化遵循高斯分布，其平均值等于彈性模量的標稱值，其標準偏差等于平均值的1%。最后，向模擬的振動響應信號中隨機添加高斯白噪聲，每個信號的水平等于相應RMS的3%，以便創建一個盡可能接近實際的模擬裝置。所有應用程序都在MATLAB中實現，硬件設施為3.80 GHz Intel Xeon E3-1275四核處理器的工作站，內存為32 GB。

下面的仿真研究基于裂紋面之間沒有接觸的假設。在沒有拉伸載荷的情況下，這種假設可能不現實，拉伸載荷會打開裂紋，對于厚度為零的裂紋也不現實。然而，當前工作的主要目的不是準確地表示裂紋附近的物理現象，而是捕捉裂紋導致的局部剛度降低對系統整體響應的影響。此外，在這2種應用中，pROMs的訓練空間決定參數搜索空間，同時選擇參數樣本的數量，以便在不犧牲精度的情況下最小化計算成本。然而，最小參數樣本的數量與追求參數的數量有關。因此，參數的增加會相應地擴大參數樣本量。

圖2 面內應變(a)和面外加速度振動測量(b)

將系統動力學投影到一個由25個基向量構成的參數化降階空間中。模型的參數化是根據裂紋特性進行的，即裂紋位置坐標xc∈[-0.15,0.15]、yc∈[-0.1,0.1]和長度lc∈[0.03，0.12]，其定義在x-y系統的投影如圖2所示，參考裂紋與圖3的相應網格。文中認為壓力是面板表面上的均勻壓力，其大小通過估計來確定。但這并不是文中方法的局限性，當使用降階模型時，它可以通過高斯過程模型[13]或廣義力項[14]輕松擴展，以適應非均勻壓力的估計。

圖3 參考網格和裂紋配置

表1 調整機身面板應用程序的參數值

使用文中方案估計的參數值如圖4(a)所示，其中算法在模擬的前幾秒鐘內收斂到實際參數值?；谒阉骺臻g的隨機探索，算法運行5次得到的相應參數估計。圖4(a)(b)分別為參數樣本在2個不同時刻的演變情況，即t=0 s和t=7 s?？梢钥闯觯鹤畛跛阉鞣植荚诳臻g的很大一部分，并隨著時間的推移而收斂到由黃色三角形標記的實際參數向量區域。t=7 s后，對參數空間重新采樣，但是搜索分布已經收斂到非常接近目標值的位置，新粒子對估計的參數的改進不大。

圖4 估計的參數值

使用pROM獲得的狀態估計值與用于生成合成振動數據的模型狀態不可直接比較，因為后者是一個全階模型，其狀態與pROM的狀態不一致。因此，根據測量和未測量位置的響應估計，可以證明預測狀態的準確性。圖5展示了點1、3處的預測位移響應，其位置如圖2所示，從正向模擬中獲得相應的噪聲信號。因此，圖2中定義的未測量點A、B、C和D處的預測位移響應如圖6所示?？梢钥闯觯罕M管信號的平均值存在一定偏差，但在4個位置都能很好地捕捉到實際的系統動力學。上述曲線圖中，刪除模擬時間的前5 s，因為響應主要由瞬態動力學控制，其量級比穩態響應大一個數量級，并且當與穩態響應一起繪制時，后者不可見。

圖6 預測位移響應

圖7 輸入估計曲線

圖8 圖形化描述

3.2 翼箱面板

第二個示例是翼箱面板，如圖9所示。面板由一塊3 mm厚的板組成，該板由2個C形槽肋支撐，肋位于2個邊緣x方向，并由2個角截面縱梁進一步加固，由虛線表示其位置。實際面板在每個肋的兩端受到彈性支承的約束，彈性支承基本上代表了互連元件的剛度，但為了簡單起見，文中認為面板完全固定在肋的兩端。這種固定的設置使系統比實際情況稍硬，其動態特性對局部變化不太敏感，其識別問題比文獻[15]中研究的問題更具挑戰性。假定面板的材料特性為線彈性，且E=73 GPa、ν=0.3、密度ρ=2 700 kg/m3，在第一和第四振型中進一步引入了2%的比例阻尼。

圖9 機翼箱面板幾何圖形和裂紋參數定義的示意

圖9還顯示了面板的尺寸(單位：m)和傳感器位置，以及第一桁條上的裂縫，該裂縫垂直于構件的縱軸。用于識別問題的pROM由15個基向量組成，從y方向開始，裂紋相對位置xc進行參數化，范圍在距縱梁中點0.2 ～2 m，或笛卡爾坐標系下0.15～0.55 m之間，長度為lc∈[0.002 5，0.022 5] m。圖10描述了用于ROM構造的面板參考離散化以及參考裂紋。在操作條件下測試面板的振動響應，假設動力學由均勻壓力驅動，該壓力的大小為高斯白噪聲過程，并沿振動方向施加在板上z軸，與之前的應用類似，文中根據量級識別該輸入p(t)。

圖10 參考離散化(a)以及參考裂紋(b)

表2 調整機翼箱面板應用程序的參數

在計算性能方面，整個算法平均需要472 s，模擬周期為30 s。算法中最耗時的部分是新系統矩陣的生成，大約需要90%的執行時間，即419 s，而剩下的53 s用于依次執行所有粒子的預測、濾波和噪聲適應步驟。后者幾乎滿足實時處理，可以通過并行處理，從而實現數字孿生。在這種情況下，可以持續檢測物理系統，并根據測得的響應更新響應狀態。另外，降低階數建模部分的實現在計算時間方面沒有進行優化，這也可以在進一步并行化時顯著提高速度。

從算法的5次不同運行中獲得的參數估計如圖11所示，根據模擬的不同時刻的估計歷史時間和參數樣本。從圖11(a)可以觀察到：在前5 s內，參數估計非常接近實際值。橙色軌道的參數估計如圖11所示，其相應顏色的粒子權重在整個仿真時間內的更新如圖12所示。在前2.5 s內，粒子快速探索搜索空間，到達目標附近。在密集重采樣事件中，當簡并閾值達到時，整個權重集中到單個粒子，從而產生高度接近1的峰值。然后，在目標附近繼續搜索，直到t=8.5 s附近，粒子被略微拉開，經過一段密集的重采樣后，在t=11.5 s附近返回目標點。此后，算法保持接近目標，只進行有限的調整，對參數估計沒有任何顯著的改進。需要強調的是，參數值在整個模擬過程中不斷更新，因此圖12中所示的權值曲線不是恒定的參數值。

圖11 參數估計結果(a)和粒子權重(b)

圖12 權值曲線

運行系統的預測響應時間如圖13所示。將觀測到應變和加速度的點1處的預測位移和點2處的加速度與圖13(a)中的測量結果進行比較。這些信號再次出現在不包含瞬態動力學的時間窗口中，從2.5 s開始。值得注意的是，在測量點處較好地捕獲振動信號，這在很大程度上是由于噪聲調諧，下面將進一步討論。另外，在圖13(b)中觀察到略有不同的圖片，其中以位移和加速度表示未測量點A和B的響應。點A的動力學在位移方面得到了充分的追蹤，在響應的峰值處出現了一定的差異。這可以部分歸因于預測裂紋參數中的誤差以及系統模型的降階，該降階是基于實時性能選擇的，以犧牲準確性為代價。在實際監控場景中，可以相應地調整系統順序，以平衡計算需求，提高精度。相反，點B處預測加速度的誤差主要與估計輸入的誤差有關，如圖14所示。雖然位移和速度估計僅取決于未觀測狀態估計的質量，但加速度預測強烈依賴于輸入估計。

圖13 預測響應時間

圖14 輸入估計

最后，為了證明測量噪聲項的協方差矩陣自適應的效果，并進一步對遺忘因子δQw的選擇值進行推理，對δQw的不同值進行實驗。通常，噪聲協方差矩陣的調整旨在實現濾波器的最佳性能，可以通過使用不同的最佳性能測試及噪聲的白度進行驗證。在此例中，使用預測殘差的頻域表示在點1處的定性檢查最優性e1(t)，圖15(a)中繪制了遺忘因子的3個不同值。圖15(b)中顯示了相應的時域殘差。為了濾除密集重采樣事件引入的偏差，重新初始化濾波器并導致殘差突然增加，圖15所示的結果是基于12.5 s后的樣本，由充分收斂的模型生成這些樣本，如圖11(a)所示。

圖15 遺忘因子測量值和預測值之間的差異

遺忘因子通?？梢赃M行更貪婪地探索，文獻[15]中提出了遺忘因子，但是，這不在文中的研究范圍。文中旨在證明協方差矩陣自適應性對算法性能的影響，并為實際應用中遺忘因子的調整提供進一步的指導。圖15中顯示的輸入、狀態和參數結果是令δQw=4，如圖15(a)所示，與在試驗中獲得的光譜相比，產生了白色光譜δQw=3和5。另外，δQw意味著協方差矩陣的適應非常緩慢，測量噪聲仍然接近初始值。另一方面，接近1的值會根據計算的殘差突然改變協方差矩陣。在此例中，由于預測和測量之間存在較大差異，小于2.5的δQw在某些過濾器中產生不穩定性，即遠離目標。相反，當遺忘因子值較大時，這些問題并沒有表現出來，但是，其預測精度較低，反映在高度相關的殘差上，與δQw=5相似。

4 結論

為了保證檢測的穩定性和準確性，并且有效處理參數的關聯性與強非線性，提出了一種基于參數降階模態的分層貝葉斯在線裂紋檢測方法。最后通過航空航天應用中的真實組件對提出的方法進行仿真，得如下結論：

(1)該算法提供了穩定和準確的參數整定估計，并且能夠有效處理參數的關聯性與強非線性，實現高精度裂紋的檢測。

(2)閾值不會影響估計的準確性，但它是關系到算法計算性能和收斂速度的決定性參數。當使用較低的閾值時，該算法需要更多步驟才能達到簡并，從而在時間上產生更稀疏的重采樣事件。

(3)雖然位移和速度估計僅取決于未觀測狀態估計的質量，但加速度預測強烈依賴于輸入估計。另外，當遺忘因子較小時，在某些過濾器中產生不穩定性；當遺忘因子值較大時，預測精度較低。