(1.中國飛行試驗(yàn)研究院測試技術(shù)研究所，陜西 西安 710089；2.中國飛行試驗(yàn)研究院飛機(jī)飛行試驗(yàn)技術(shù)研究所，陜西 西安 710089)

0 引言

近幾十年來，隨著航空航天領(lǐng)域的快速發(fā)展，復(fù)雜工程結(jié)構(gòu)內(nèi)部的應(yīng)力變得越來越高，加之使用環(huán)境也越來越惡劣，結(jié)構(gòu)中的微小損傷有可能會導(dǎo)致結(jié)構(gòu)損傷失效事故層。因此對于這些結(jié)構(gòu)需要進(jìn)行實(shí)時(shí)在線監(jiān)控。

由于Lamb波具有傳播距離遠(yuǎn)、對結(jié)構(gòu)微小損傷敏感等特點(diǎn)，基于Lamb波的主動監(jiān)測技術(shù)被認(rèn)為是有效的結(jié)構(gòu)損傷檢測方法，并成為了國際研究熱點(diǎn)之一[1]。在基于Lamb波的結(jié)構(gòu)健康監(jiān)測中，通過對損傷散射信號的處理分析，能夠?qū)崿F(xiàn)對結(jié)構(gòu)損傷的監(jiān)測與評估[2]。Lamb波在結(jié)構(gòu)中傳播時(shí)具有無方向性、頻散及多模式等特點(diǎn)，這就為損傷特征的提取增加了很多困難。同時(shí)，目前結(jié)構(gòu)健康監(jiān)測技術(shù)發(fā)射的波形固定，檢測環(huán)境發(fā)生改變時(shí)，僅依靠對結(jié)果信息的分析難以達(dá)到理想的檢測效果。因此激勵(lì)信號波形選擇在提升結(jié)構(gòu)的損傷檢測效率方面有著重要的研究意義。

1 Lamb波激勵(lì)信號

基于Lamb波的主動結(jié)構(gòu)健康監(jiān)測方法中，通常采用寬帶激勵(lì)或窄帶激勵(lì)的方式激發(fā)出被測結(jié)構(gòu)中的Lamb波。本文的研究工作使用窄帶Lamb波信號，因此下文中將對窄帶激勵(lì)進(jìn)行闡述。

通常情況下，選擇調(diào)制正弦信號作為Lamb波窄帶激勵(lì)信號，其表達(dá)式[3]如下：

sin2πfct

(1)

式中，n為信號波峰個(gè)數(shù)；fc為信號中心頻率；A為信號的幅度調(diào)制；H為Heaviside階梯函數(shù)。

圖1為中心頻率為300 kHz的五周期正弦調(diào)制信號的時(shí)域波形及其傅里葉變換后的頻域波形。

圖1 五周期激勵(lì)信號及其傅里葉變換

根據(jù)Lamb波頻散特性，由于部分模態(tài)具有截止頻率，為了減少激發(fā)的模態(tài)數(shù)量，避免多模態(tài)的混疊，因此需要選擇合適的激勵(lì)信號中心頻率，以確保被測結(jié)構(gòu)中的Lamb波只存在S0與A0兩種模式。激勵(lì)信號的周期數(shù)一般在3.5～13.5之間較為合適[4]。

由于Lamb波的自身特點(diǎn)會造成其傳播過程較為復(fù)雜，因此通過數(shù)值計(jì)算的方法模擬Lamb波激勵(lì)信號傳播一段時(shí)間后的波形變得十分有意義，本節(jié)將對此展開介紹。

設(shè)激勵(lì)信號為F(t)，激勵(lì)信號傳播一段距離后的信號為u(x，t)，其中，x為傳播距離，t為傳播時(shí)間。激勵(lì)信號在t=0時(shí)刻作用于坐標(biāo)系上x=0的位置，二者之間的關(guān)系可由以下初始條件確定：

u(x，0)=0

(2)

u(0，t)=F(t)

(3)

由文獻(xiàn)[5]可推導(dǎo)出結(jié)構(gòu)內(nèi)任意一點(diǎn)的信號為：

(4)

式中，k(ω)為角頻率為ω時(shí)的波數(shù)；F(ω)是f(t)的傅里葉變換信號。

當(dāng)激勵(lì)信號遇到類似于邊界或損傷的反射源時(shí)，反射后的信號為g(t)，信號接收點(diǎn)接收到的第j個(gè)反射源的反射信號為gj(t)，設(shè)激勵(lì)點(diǎn)到反射源的距離為dj，反射源到接收點(diǎn)的距離為xj，反射系數(shù)為Aj(ω)，則gj(t)可表示為：

(5)

在MATLAB里進(jìn)行編程模擬激勵(lì)信號傳播一段時(shí)間后的波形，以中心頻率為300 kHz的五周期正弦調(diào)制信號為激勵(lì)信號，模擬激勵(lì)信號傳播一段距離(10.7 cm)，波形如圖2所示，現(xiàn)只討論其中S0直達(dá)波的波達(dá)時(shí)間。

從圖3可以看出，該算法所得數(shù)值模擬得到的直達(dá)波波形與在ANSYS中進(jìn)行仿真模擬得到的直達(dá)波波形，除個(gè)別峰值外，每個(gè)周期峰值的波達(dá)時(shí)間基本吻合，可認(rèn)為該算法在一定條件下，可準(zhǔn)確模擬激勵(lì)信號在頻散情況下的傳播。這為后續(xù)的波形選擇算法奠定了基礎(chǔ)。

圖2 模擬傳播一段距離后的波形

圖3 數(shù)值解與仿真模擬S0模式對比

2 波形選擇

無論是在雷達(dá)領(lǐng)域還是在結(jié)構(gòu)健康監(jiān)測系統(tǒng)，環(huán)境噪聲的干擾一直是重要的影響因素，為了排除噪聲干擾，提高對目標(biāo)的測量精度，近年來雷達(dá)領(lǐng)域興起了對認(rèn)知雷達(dá)的研究，而結(jié)構(gòu)健康監(jiān)測也引入眾多的濾波方法。同時(shí)，結(jié)構(gòu)健康監(jiān)測系統(tǒng)本身與雷達(dá)系統(tǒng)具有相似的特點(diǎn)，所以可以在結(jié)構(gòu)健康監(jiān)測領(lǐng)域中引入雷達(dá)系統(tǒng)中的先進(jìn)方法。

本節(jié)旨在給出一種提高靜態(tài)目標(biāo)檢測精度的方法，將新興的認(rèn)知雷達(dá)跟蹤動態(tài)目標(biāo)[6]轉(zhuǎn)化為結(jié)構(gòu)損傷檢測的靜態(tài)目標(biāo)。首先從信號的回波模型出發(fā)，由信號估計(jì)理論推導(dǎo)測量噪聲協(xié)方差的卡拉美羅下界(Cramer-Rao Low Bound，CRLB)，得出波形參數(shù)對其影響與意義，然后在Kalman濾波的基礎(chǔ)上增加波形選擇，選擇合適波形使濾波協(xié)方差最小，從而提高對目標(biāo)狀態(tài)的檢測精度。

2.1 信號參量估計(jì)

2.1.1 回波信號模型

Lamb波在板結(jié)構(gòu)中傳播形式較為復(fù)雜，為了降低計(jì)算量，本節(jié)提出激勵(lì)信號傳播一段距離后波形的簡化計(jì)算模型。激發(fā)出的Lamb波信號f(t)，所激發(fā)的場量用u(x，t)表示，其中t與x分別為激勵(lì)信號傳播時(shí)間與傳播距離，激勵(lì)器處可表示為：

u(x，t)=f(t)

(6)

板中任意一點(diǎn)處：

(7)

根據(jù)圖3數(shù)值解與仿真模擬S0模式對比，MATLAB仿真得出的傳播一段距離后的S0波形與ANSYS仿真情況下得出的S0波形相似度極高，即MATLAB仿真波形與及ANSYS仿真波形的頻散特性基本相同。這里我們使用MATLAB數(shù)值模擬傳播一段距離后的S0波形，并與激勵(lì)信號作對比。圖4為傳播0.27 cm的S0波形與激勵(lì)信號的對比。

圖4 MATLAB數(shù)值模擬的S0波形與激勵(lì)信號的對比

同時(shí)，使用ANSYS有限元軟件仿真模擬激勵(lì)信號傳播一段距離后的S0波形與激勵(lì)信號對比如圖5所示。

圖6顯示的是在損傷長度為4 cm裂縫的條件下，激勵(lì)信號與ANSYS仿真模擬的損傷散射信號進(jìn)

圖5 ANSYS仿真模擬S0波形與激勵(lì)信號的對比

行的對比。在對幅值進(jìn)行歸一化處理后，能夠看出損傷散射信號與激勵(lì)信號的波形較為相似。

圖6 ANSYS仿真模擬S0波形與激勵(lì)信號的對比

通過圖5、圖6對比可以得出，在小頻散的情況下，傳播一段距離后S0模式波形和損傷散射波與激勵(lì)信號的波形十分相似。這里提出一個(gè)簡化的計(jì)算模型，接下來在本文的研究中我們使用式(8)代替式(7)。

(8)

其中，θ1表示損傷或者邊界反射后幅值的衰減程度，θ2為反射信號的時(shí)延。

2.1.2 無偏估計(jì)量的CRLB

在信號處理中，我們通常要根據(jù)觀測信號來估計(jì)信號的某些特征參量，這就需要用到信號參量估計(jì)的相關(guān)知識。在結(jié)構(gòu)的損傷檢測中，本文將著重研究信號的振幅、時(shí)延的聯(lián)合估計(jì)，即多參量估計(jì)。本文采用的是最大似然估計(jì)，最大似然估計(jì)對于解決復(fù)雜估計(jì)問題具有較好的性能。

在本文中，對接收到的信號進(jìn)行觀測時(shí)，我們只研究信號的幅值及時(shí)延，故設(shè)觀測矢量x=[x1x2]T，由于噪聲的干擾造成的觀測結(jié)果不準(zhǔn)確，待估計(jì)參量設(shè)為θ=[θ1θ2]T，則以θ為參量的概率密度函數(shù)為f(x|θ)。在已觀測到的觀測值x=x1的條件下，似然函數(shù)f(x=x1|θ1)反映了θ1取各個(gè)值的可能性大小。估計(jì)問題其本質(zhì)是根據(jù)觀測值確定未知量，換言之，根據(jù)已知觀測值x=x1來估計(jì)未知量θ1的值。

(9)

若似然函數(shù)是可導(dǎo)函數(shù)，那么最大似然估計(jì)的必要條件為：

(10)

我們假定真實(shí)信號為s(t，θ)，s(t，θ)是關(guān)于t的連續(xù)函數(shù)，θ=[θ1，θ2，…，θM]T，包含M個(gè)未知參量，θi可以是信號的頻率、時(shí)延、幅度等。在一段觀測時(shí)間[0，T]內(nèi)，觀測到的混入噪聲的信號為：

x(t)=s(t，θ)+n(t)，0≤t≤T

(11)

式中，n(t)是功率譜密度為N0/2的高斯白噪聲。

在得到一個(gè)估計(jì)量后，為了衡量其性能是否達(dá)到最佳，且使參量估計(jì)為有效估計(jì)，我們引入克拉美羅下界[7]來揭示無偏估計(jì)量的估計(jì)誤差的最小值。

假定概率密度f(x|θ)滿足正則條件：

(12)

(13)

式中，

(14)

(15)

I-1(θ)即為無偏估計(jì)量的CRLB，其倒數(shù)I(θ)稱之為費(fèi)希爾信息。

如果有效估計(jì)量存在，即存在達(dá)到CRLB的估計(jì)量，那么這個(gè)估計(jì)量必定是最大似然估計(jì)，此時(shí)的參量估計(jì)具有很好的性能。這里CRLB會隨發(fā)射的波形參數(shù)的變化而改變，是下文中波形選擇的關(guān)鍵。

2.1.3 測量噪聲方差的CRLB

由信號參量估計(jì)理論，我們接下來推導(dǎo)基于簡化模型的測量噪聲方差的卡拉美羅下界。

將(8)式帶入(11)式，信噪比表示為SNR。

x(t)=

(16)

由(10)式：

(17)

(18)

由(14)式：

(19)

則無偏估計(jì)量的CRLB為：

(20)

計(jì)算可得測量噪聲協(xié)方差的CRLB：

(21)

根據(jù)測量噪聲協(xié)方差的CRLB可以得出以下結(jié)論：

1)信噪比越大，信號時(shí)域脈沖寬度(這里對應(yīng)于n/fc)越小，信號的時(shí)延估計(jì)方差就越小，精度也就越高。

2)對于信號的幅值，信噪比越大，信號的帶寬越小(時(shí)域脈沖寬度越大，即n/fc越大)，其估計(jì)方差就越小，精度也就越高。

2.2 基于Kalman濾波的波形選擇

Kalman濾波理論是由Wiener波理論的發(fā)展而發(fā)展起來的。隨著Kalman濾波理論的發(fā)展，它在隨機(jī)過程參數(shù)估計(jì)中的應(yīng)用越來越廣泛，其后來廣泛的應(yīng)用于解決各種最優(yōu)濾波問題[8]。在本節(jié)中，在Kalman濾波基礎(chǔ)上添加目標(biāo)狀態(tài)的一步預(yù)測是算法的關(guān)鍵內(nèi)容。

2.2.1Kalman濾波系統(tǒng)模型

首先引入一個(gè)離散控制過程的系統(tǒng)，其狀態(tài)方程可表示為：

X(k)=AX(k-1)+W(k)

(22)

離散時(shí)間系統(tǒng)的測量方程可表示為：

Z(k)=HX(k)+V(k)

(23)

其中，A為狀態(tài)轉(zhuǎn)移矩陣，由于本文的研究對象為檢測結(jié)構(gòu)的損傷狀況，且結(jié)構(gòu)的損傷狀態(tài)不隨時(shí)間而改變，故A為單位矩陣。H為測量矩陣，由于測量值反映了結(jié)構(gòu)的損傷狀況，故同樣為單位矩陣。W(k)和V(k)對應(yīng)的協(xié)方差分別為Q和R，分別表示預(yù)測和測量過程中的噪聲。

首先利用系統(tǒng)的過程模型預(yù)測下一狀態(tài)，基于系統(tǒng)的上一狀態(tài)預(yù)測現(xiàn)在狀態(tài)：

X(k|k-1)=AX(k-1|k-1)

(24)

X(k-1|k-1)是上一狀態(tài)最優(yōu)的結(jié)果。用P表示X(k|k-1)的預(yù)測協(xié)方差為：

P(k|k-1)=P(k-1|k-1)+Q

(25)

Kalman增益K為：

K(k)=P(k|k-1)/(P(k|k-1)+R)

(26)

當(dāng)前時(shí)刻目標(biāo)狀態(tài)的最優(yōu)化估算值X(k|k)：

X(k|k)=

X(k|k-1)+K(k)(Z(k)-X(k|k-1))

(27)

更新X(k|k)的協(xié)方差以為保證Kalman濾波器繼續(xù)運(yùn)行下去，則：

P(k|k)=(E-K(k))P(k|k-1)

(28)

而如果要在Kalman濾波的基礎(chǔ)上增加波形選擇部分，我們需要更進(jìn)一步計(jì)算出濾波協(xié)方差的一步預(yù)測。

預(yù)測協(xié)方差的一步預(yù)測為：

P′(k+1|k)=P(k|k)+Q

(29)

Kalman增益的一步預(yù)測為：

K′(k+1)=P′(k+1|k)/(P′(k+1|k)+R)

(30)

濾波協(xié)方差的一步預(yù)測為：

P′(k+1|k+1)=(E-K′(k+1))P′(k+1|k)

(31)

濾波協(xié)方差的一步預(yù)測滿足一定波形選擇準(zhǔn)則后，所選取的波形參數(shù)會對下一時(shí)刻的目標(biāo)狀態(tài)的最優(yōu)化估算值產(chǎn)生影響，進(jìn)而使得基于Kalman濾波的目標(biāo)狀態(tài)的最優(yōu)化估計(jì)循環(huán)計(jì)算下去。

2.2.2 波形選擇準(zhǔn)則

為使每一時(shí)刻狀態(tài)估計(jì)誤差的均方差最小，引入均方誤差最小準(zhǔn)則作為波形選擇的依據(jù)，這一準(zhǔn)則使?fàn)顟B(tài)空間每個(gè)維數(shù)上誤差的平方和最小。即：

(32)

Zk為k時(shí)刻的測量值集合，濾波協(xié)方差定義為：

P(k|k)=

(33)

等號兩邊同時(shí)取跡，則：

Tr{P(k|k)}

(34)

則：

(35)

由(34)式與(35)式，可得：

Tr{Pk|k(αk)}=

Tr{Pk|k-1-Pk|k-1(Pk|k-1+R(αk))-1Pk|k-1}

(36)

根據(jù)(36)式可以得出，由預(yù)測協(xié)方差和測量噪聲方差能夠求得k時(shí)刻的濾波協(xié)方差的跡。由此可知，在均方誤差最小準(zhǔn)則這一波形選擇的依據(jù)中，測量噪聲方差是與波形的參數(shù)相關(guān)的變量，也正是前文中所提到的是波形選擇的關(guān)鍵所在。

2.2.3 基于Kalman濾波的波形選擇流程圖

如圖7所示，圖中左邊的兩部分為常規(guī)Kalman濾波，右邊的部分為所增加的波形選擇模塊。

圖7 基于Kalman濾波的波形選擇流程圖

根據(jù)流程圖，可以看出在整個(gè)波形選擇過程中，包含有兩個(gè)循環(huán)計(jì)算過程，其一是狀態(tài)估計(jì)協(xié)方差的循環(huán)計(jì)算過程，這其中包括了協(xié)方差的一步預(yù)測與波形選擇；其二是目標(biāo)狀態(tài)的估計(jì)，包括了目標(biāo)狀態(tài)的預(yù)測方程與測量方程，其狀態(tài)估計(jì)的對象為前文中所提到的信號的時(shí)延與幅值，信號的時(shí)延與幅值精度是反映本算法是否可行的關(guān)鍵指標(biāo)。

3 基于MATLAB的模擬計(jì)算

本節(jié)通過在MATLAB中進(jìn)行模擬計(jì)算，驗(yàn)證基于Kalman濾波的波形選擇這一算法的有效性，并與不帶有波形選擇模塊的Kalman濾波結(jié)果對比，驗(yàn)證上文所述算法的計(jì)算精度。

在模擬計(jì)算環(huán)境時(shí)，我們將上一時(shí)刻最優(yōu)值作為該時(shí)刻的預(yù)測，并且假定測量目標(biāo)的真實(shí)狀態(tài)對應(yīng)的信號時(shí)延為2.4×10-5s，對應(yīng)的幅值變化為0.4。目標(biāo)初始狀態(tài)設(shè)為：時(shí)延3.5×10-5s、幅值變化為0.5。為了檢驗(yàn)所研究的波形選擇方法在不同強(qiáng)度的噪聲下的濾波效果，信噪比分別設(shè)為15 dB與25 dB，初始發(fā)射的波形為五周期漢寧窗調(diào)制函數(shù)，其波形參數(shù)為n=5以及fc=300 kHz。根據(jù)Lamb波的傳播特性以及頻散曲線，波形庫的建立周期數(shù)的取值區(qū)間為3～8，中心頻率取值區(qū)間為[250 kHz，350 kHz]，步長為10 kHz，整個(gè)波形庫包含有66種波形。

在經(jīng)過50次迭代計(jì)算后，將不含有波形選擇的Kalman濾波與帶有波形選擇模塊的Kalman濾波作對比，比較信號的時(shí)延誤差與幅值變化的誤差。如圖8與圖9所示，虛線、點(diǎn)線、實(shí)線分別表示帶有波形選擇的Kalman濾波、不帶有波形選擇的Kalman濾波、參量的真實(shí)值。其中圖10為信號的時(shí)延估計(jì)誤差的對比，圖11為信號的幅值變化估計(jì)誤差的對比。

圖8 時(shí)延估計(jì)對比圖(SNR=25dB)

圖9 幅值估計(jì)對比圖(SNR=25dB)

圖10 多次平均的時(shí)延估計(jì)對比圖(SNR=25dB)

圖11 多次平均的幅值估計(jì)對比圖(SNR=25dB)

為了檢驗(yàn)算法在不同強(qiáng)度噪聲下的適用性，當(dāng)信噪比為15 dB時(shí)，信號的時(shí)延估計(jì)與幅值估計(jì)顯示在圖12和圖13中。

圖12 時(shí)延估計(jì)對比圖(SNR=15dB)

圖13 幅值估計(jì)對比圖(SNR=15dB)

根據(jù)圖14和圖15，在經(jīng)過Kalman濾波后，信號的時(shí)延估計(jì)與幅值估計(jì)的精度都較初始值有了較明顯的提高。在經(jīng)過一定次數(shù)的迭代計(jì)算之后，在初始值誤差比較大的情況下，Kalman濾波仍能將誤差顯著的減小。并且，無論是信號的時(shí)延估計(jì)還是信號的幅值估計(jì)，虛線(帶有波形選擇的Kalman濾波)比點(diǎn)線(不帶有波形選擇的Kalman濾波)更接近實(shí)線(真實(shí)值)，這表明了帶有波形選擇的Kalman濾波的計(jì)算結(jié)果要好于不帶有波形選擇的Kalman濾波。

圖14 多次平均的時(shí)延估計(jì)對比圖(SNR=15dB)

圖15 多次平均的幅值估計(jì)對比圖(SNR=15dB)

根據(jù)圖8與圖12、圖9與圖13的對比，在噪聲強(qiáng)度改變時(shí)，經(jīng)過50次的迭代計(jì)算之后，帶有波形選擇的Kalman濾波精度均較初值有所提高，但是當(dāng)噪聲強(qiáng)度比較大時(shí)，波形選擇后結(jié)果的誤差相對就會大一些。信號的時(shí)延估計(jì)與幅值估計(jì)精度的提升驗(yàn)證了波形選擇算法在信噪比不同時(shí)均具有良好的效果。

4 結(jié)論

研究基于Kalman濾波的波形自適應(yīng)選擇，從建立簡化的回波信號模型出發(fā)，經(jīng)過信號估計(jì)理論中的最大似然估計(jì)，詳細(xì)推導(dǎo)了基于簡化回波模型的測量噪聲協(xié)方差的卡拉美羅下界，然后在Kalman濾波的基礎(chǔ)上增加波形選擇模塊，并以均方誤差最小準(zhǔn)則作為波形選擇的依據(jù)，最后通過模擬計(jì)算驗(yàn)證了經(jīng)過波形選擇后信號的時(shí)延與幅值精度有了較大的提高。波形選擇依據(jù)的是均方誤差最小準(zhǔn)則，由于是在原有激勵(lì)信號波形的基礎(chǔ)上進(jìn)行波形選擇，且在波形選擇的準(zhǔn)則下，所選擇的新波形對應(yīng)的Kalman濾波的協(xié)方差會比原有激勵(lì)信號波形對應(yīng)的Kalman濾波協(xié)方差更小，這是帶有波形選擇模塊的Kalman濾波精度更高的根本原因。證明了波形選擇方法是提高噪聲干擾下的結(jié)構(gòu)損傷定位精度的另外一種思路。