基于參數(shù)降階模態(tài)的分層貝葉斯在線裂紋檢測

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

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

0 前言

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

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

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

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

1 問題描述

1.1 參數(shù)化降階模型

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

(1)

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

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

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

(2)

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

V(θ)TSpp(t)

(3)

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

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

(4)

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

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

1.2 狀態(tài)空間表示

為了將建模的系統(tǒng)動力學(xué)與實際振動測量值融合，需要將參數(shù)化的降階運動方程即公式(3)轉(zhuǎn)換為連續(xù)的時間-狀態(tài)-空間表示，其表示為

(5)

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

(6)

(7)

(8)

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

G(θ)=

(9)

(10)

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

2 遞階狀態(tài)輸入?yún)?shù)估計

2.1 狀態(tài)輸入估計

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

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

(11)

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

(12)

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

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

(16)

(17)

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

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

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

(18)

(19)

圖1 提出算法的結(jié)構(gòu)框圖

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

(20)

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

(21)

(22)

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

(23)

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

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

2.2 參數(shù)估計

(24)

當(dāng)處于當(dāng)前步驟k，可以假設(shè)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)，根據(jù)貝葉斯規(guī)則，在式(24)中進行替換后，得出以下遞推公式：

(25)

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

(26)

2.3 更新策略

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

(27)

(28)

(29)

算法1：用于參數(shù)和輸入狀態(tài)估計僅輸出的分層濾波器。

(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μ分別是更新協(xié)方差矩陣的排名第一和第nθ，p的學(xué)習(xí)率。算法1中記錄了實現(xiàn)分層濾波器的詳細步驟。

3 仿真案例

3.1 機身面板

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

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

圖2 面內(nèi)應(yīng)變(a)和面外加速度振動測量(b)

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

圖3 參考網(wǎng)格和裂紋配置

表1 調(diào)整機身面板應(yīng)用程序的參數(shù)值

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

圖4 估計的參數(shù)值

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

圖6 預(yù)測位移響應(yīng)

圖7 輸入估計曲線

圖8 圖形化描述

3.2 翼箱面板

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

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

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

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

表2 調(diào)整機翼箱面板應(yīng)用程序的參數(shù)

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

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

圖11 參數(shù)估計結(jié)果(a)和粒子權(quán)重(b)

圖12 權(quán)值曲線

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

圖13 預(yù)測響應(yīng)時間

圖14 輸入估計

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

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

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

4 結(jié)論

為了保證檢測的穩(wěn)定性和準確性，并且有效處理參數(shù)的關(guān)聯(lián)性與強非線性，提出了一種基于參數(shù)降階模態(tài)的分層貝葉斯在線裂紋檢測方法。最后通過航空航天應(yīng)用中的真實組件對提出的方法進行仿真，得如下結(jié)論：

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

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

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