基于粒子群優化的超聲共振譜法測量材料彈性常數

沈 飛，樊 璠，王 蕊，張 強，牛海軍

(北京航空航天大學 生物與醫學工程學院，北京 100191)

超聲共振譜法(Resonant Ultrasound Spectroscopy, RUS)是20世紀90年代發展起來的一種材料力學特性測量方法，被物理學家認為是測量高Q值(品質因數)固體材料彈性系數的最準確方法[1-2]。其基本原理是運用超聲激勵樣本使其產生自由振動，再根據測得的多個固有共振頻率，利用反演方法計算材料的彈性系數[2-3]。方法的獨特之處在于可用于尺寸極小(小于1 mm)的樣本的無損檢測，一次試驗就能估計樣本的全部彈性系數和力學參數，且具有可重復性[1]。

RUS方法求取彈性常數的過程一般轉化為非線性最小二乘問題，利用梯度優化(如Levenberg-Marquardt，LM方法)對問題進行求解[1-4]。但是該方法對于高阻尼材料(Q值較低，如巖石材料、生物硬組織等)彈性常數的估計存在缺陷，主要原因是這類材料的譜線平緩、共振頻率求取困難。盡管一些研究者結合線性預測法提出了平緩共振譜線中共振頻率的提取技術，但是低Q值材料的試驗共振頻率丟失現象仍舊無法避免，這直接導致了試驗共振頻率與理論共振頻率的配對困難，也增加了優化過程的難度[4-5]。這對彈性常數初始值的設置提出了要求，偏差稍大就會導致算法無法收斂，彈性常數計算無效。

粒子群優化是一種基于數量的進化算法，具有優良的全局尋優能力，被廣泛應用于動態優化、信號處理等領域[6]。為了獲得彈性常數的點估計，筆者提出了基于粒子群優化算法的RUS求解方案，可快速實現寬泛彈性常數初值設置下的彈性常數估計；并利用低Q值骨模擬材料試驗進行了驗證。

1 檢測方法

RUS方法包括3個步驟：① 利用超聲掃頻試驗獲得材料的共振譜線；② 使用線性預測技術從共振譜線中提取共振頻率；③ 不斷調整彈性常數Cij，直到計算得到的共振頻率fcal與試驗提取的共振頻率fexp相一致(見圖1，圖中M和K分別為樣本的質量矩陣和剛度矩陣,α為共振模態，ω為角頻率)。

圖1 超聲共振譜方法流程圖

1.1 理論共振頻率計算與一致性準則

三維固體材料在零應力邊界條件下，其共振頻率由形狀、密度以及彈性常數決定。一般情況下，固體的形狀和密度易于測量，RUS技術通常采用Rayleigh-Ritz(瑞利-里茨)法獲取理論共振頻率f與彈性常數的關系。

(2πf)2Mα=Kα

(1)

已知彈性常數Cij，f以及相應的共振模態α可以通過求解式(1)實現。

但已知共振頻率f無法直接反算彈性常數Cij，因此基于逆問題理論對問題建模。定義彈性常數C的后驗概率密度函數，如式(2)所示。

(2)

式中：ρ為理論共振模態與試驗共振模態的配對關系，表示為離散向量的形式，如ρ=[1,3,6],說明試驗測量的前3個共振頻率分別與理論共振頻率的第1,3,6個對應;p為概率。

C與ρ的先驗相互獨立，似然函數p(f|C,ρ)表示試驗共振頻率與彈性常數和配對的關系，可以使用多元高斯分布來表示[7]。這里假設C的先驗服從多元對數高斯分布，進一步假設ρ的先驗如式(3)所示。

(3)

由于高階共振頻率一般較集中，故配對低階共振頻率比配對高階共振頻率更加可信，即ρ中元素應傾向于取小值。根據逆問題的貝葉斯理論，后驗概率密度最大值所對應的C便是材料彈性常數的最優點估計。

最大化式(2)的后驗概率密度，采用分步計算的方式，首先對后驗概率密度函數取負對數，估計理論共振頻率與試驗共振頻率的最優配對ρ，即ρ表示為C的函數，再優化彈性常數C。

考慮到ρ的變動對似然函數p(f|C,ρ)取值的影響要遠大于其對先驗函數p(ρ)的影響，故估計最優配對ρ近似于最小化-logp(f|C,ρ)。顯然這是一個線性和指派問題，可采用匈牙利算法進行求解[8]。彈性常數C的優化則可以借助全局優化器粒子群算法實現。

1.2 粒子群優化估計彈性常數

考慮到模型的高度非線性，以及負對數后驗概率密度函數具有眾多極值點，采用粒子群優化的方法尋找全局最小值。粒子群優化(PSO)是一種基于數量的隨機優化算法，函數上每一個可能的采樣點稱為一個粒子，基于有限粒子間的信息交流來尋找問題的最優解。為了加速粒子的探索能力，在原始PSO中引入了梯度下降和變異操作(見圖2)。

圖2 具有梯度和變異操作的粒子群優化算法流程圖

梯度下降：對后驗概率密度函數計算一階導數和二階導數，并將模擬牛頓法應用到迭代算法中。盡管這一方法有助于提高粒子的局部尋優能力，但這是以計算導數為代價的，因此會增加計算開銷。

每次迭代時，計算最優粒子所對應的理論共振頻率與試驗共振頻率的相對均方誤差，當均方誤差小于一定值時(試驗中設置為0.5%)停止迭代。

2 試驗過程

2.1 試驗系統以及試驗共振頻率的提取

為了驗證算法在實際材料測量上的有效性，選擇了一種皮質骨模擬材料(短纖維環氧樹脂，橫觀各向同性，SAWBONES?)，Q值約為25，將材料切割、打磨成標準長方體，骨樣品的尺寸(長×寬×高)為5.84 mm×6.74 mm×6.64 mm,質量為428.6 mg。樣品的對角被夾持在兩個超聲探頭中間，接收探頭端連接電荷放大器將接收到的振動信號放大后輸入網絡分析儀，并上傳到計算機,試驗系統框圖如圖3所示。計算機對共振譜線進行處理，提取試驗共振頻率。共振頻率提取細節參見文獻[9]。圖4展示了測量到的共振譜圖以及通過線性預測法獲得的試驗共振頻率。

圖3 試驗系統框圖

圖4 皮質骨模擬材料在70 kHz～300 kHz間的譜線以及提取的16個試驗共振頻率點

2.2 試驗結果

彈性常數的先驗來自于前人的皮質骨彈性常數測量結果。

將粒子數設置為50，參數搜索空間設置為先驗值的99.7%置信區間，最終優化得到骨模擬材料的彈性常數，并將計算出的楊氏模量與制造商進行拉壓試驗提供的楊氏模量進行對比，結果見表1。通過彈性常數計算楊氏模量的方法參見文獻[10]。結果顯示，利用提出的方法計算出的橫向楊氏模量Et和縱向楊氏模量Ea與制造商提供的數值具有很好的一致性。

表1 骨模擬材料的彈性常數以及楊氏模量 GPa

同時采用了Levenberg-Marquardt算法進行彈性常數的計算，將初值設置為C11=16 GPa，C33=22 GPa，C12=8 GPa，C13=8 GPa，C44=5 GPa，最終結果與采用粒子群算法計算出的結果相一致。以上兩個對比證明了利用粒子群算法可以有效估計出低Q值材料的彈性常數以及工程模量。

需要指出的是，使用LM算法時所設置的彈性常數初值與真實彈性常數非常接近，如果使用皮質骨的統計先驗均值作為初值，則無法收斂到該解。這進一步突顯出了粒子群算法在初值設置方面的優越性。此外，粒子群算法在計算速度方面表現良好，采用MATLAB軟件編程，算法整個優化過程耗時約60 s，遠小于基于馬爾科夫鏈蒙特卡洛方法所需要的時間(幾十分鐘)。

3 結語

RUS作為一種測量材料彈性常數高精度的試驗方法，由于模型本身的非凸性以及存在試驗共振頻率丟失的現象，彈性常數的計算結果嚴重依賴于初始值設置。針對此問題，提出了基于PSO的優化算法，可有效估計低Q值材料的彈性常數。該算法具有更高的魯棒性和全局優化能力，在初始值隨機設置時表現出了更優的效果，并且相比于馬爾科夫鏈蒙特卡洛方法，計算速度顯著提升。此外，算法具有普適性，不但可以估計貝葉斯最大后驗概率點，同樣適用于最小二乘模型的求解。在未來的工作中，算法有望在效率上進一步提升，如：引入遺傳算法的交叉運算，以記錄粒子探索過的局部區域；減少計算過程中的無效變異，以節省計算時間。