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

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

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

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

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

1 檢測方法

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

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

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

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

(2πf)2Mα=Kα

(1)

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

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

(2)

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

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

(3)

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

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

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

1.2 粒子群優(yōu)化估計彈性常數(shù)

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

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

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

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

2 試驗過程

2.1 試驗系統(tǒng)以及試驗共振頻率的提取

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

圖3 試驗系統(tǒng)框圖

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

2.2 試驗結果

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

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

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

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

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

3 結語

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