四元數聯合傳遞率在結構狀態檢測中的應用

任同群,鄭培珍,徐征,王曉東

(1.大連理工大學精密與特種加工教育部重點實驗室,遼寧 大連116024;2.大連理工大學 微納米技術及系統遼寧省重點實驗室,遼寧大連116024)

0 引言

在實際工程中,多數情況下結構的激勵未知或不易測得,這導致無法準確識別結構的模態參數,進而限制了基于結構參數變化的損傷檢測方法的應用。相比之下,直接基于振動響應信號分析的狀態檢測方法不影響結構正常運行且費用較低,已成為目前結構健康檢測的研究熱點,也得到了廣泛的工程應用,其基本過程包括損傷特征提取與分類[1-4]。

在眾多所提取的損傷特征中,振動傳遞率反映了結構兩測點之間的全部力學特性,可表示為結構剛度、阻尼和質量等物理參數的函數。Fauziyah等人[5]利用加速度計陣列測量多泵間的振動傳遞率,實現了多泵的故障診斷。Dong等人[6]基于振動傳遞率對隨機激勵下的螺栓松動監測進行了研究,克服了常用的振動方法在發現螺栓局部松動方面的不足。楊斌[7]采用相空間重構方法進行振動傳遞率函數分析,重構振動傳遞率函數矩陣,并以其奇異值熵作為損傷指標,提高了檢測方法的抗干擾性。Luo等人[8]以加權傳遞率作為損傷指標評估結構狀態,提高了TF的信噪比。Chen等人[9]利用神經網絡,提取基于振動傳遞率的損傷診斷指標,并驗證了其正確性。

在上述研究中,振動傳遞率僅基于單通道信號計算獲得。在工程應用中,測試信號是實際振動在傳感器敏感方向的投影。因此,基于單通道信號的方法在處理空間振動問題時容易陷入困難。為解決此問題,通常采用多加速度計或三軸加速度計獲取空間振動信號,進行單通道信號處理及特征提取后再進行特征融合。但是,這種處理方法忽略了信號通道之間的關聯性,融合后會帶來失真甚至是虛假的結果。因此,有必要將多通道信號視為矢量信號,進行同時處理,以保證通道間的關聯性。四元數三個虛部的運算是等價的,因此非常適合表示三維矢量信號。例如,Rehman和Mandic[10]研究了三維信號的經驗模態分解,并證明三通道信息聯合處理時會獲得更好的精度。Tong[11]采用四元數K-LT對故障異步電機進行特征提取,但處理的數據僅限于時域四元數振動序列。Yi和Lv等人[12-13]描述了四通道信號的相關性,提出一種基于凸優化的四元數奇異譜分析方法。Ren等人[14-15]采用四元數,構建三軸信號的振動傳遞率,并驗證了K-LT法和DQWT法,取得了較好的結果。此外,Ell和Sangwine[16-19]針對四元數傅立葉變換和相關分析問題進行了合作研究,相關研究成果為四元數頻域分析奠定了理論基礎。

本文提出了一種基于QTJT的狀態檢測方法,先后采用K-LT和DQWT技術提取狀態特征指標向量,然后定義狀態特征指標向量之間的歐氏距離并將其作為狀態指標,針對軌道結構扣件松脫及內部溫度力的變化,進行了狀態檢測與識別。通過在無砟軌道試驗平臺上開展實際實驗,驗證該方法的有效性。

1 四元數聯合傳遞率相關理論

1.1 四元數傳遞率計算

對于單通道信號,振動傳遞率定義為結構兩測試點的時域序列xi(輸出)和xj(輸入)的頻譜比值。它描述了振動(振幅和相位)如何在結構的兩個測試點之間進行傳遞。假設一個n自由度線性結構受外力f(t)=[f1(t),f2(t),…,fn(t)]的激勵,如果f(t)具有形如白噪聲的均勻譜密度,則fk(t),k=1,2,…,n的傅立葉變換則為一個常數D。此時,由節點j到i的振動傳遞率計算公式為

式中:Ai(w),Aj(w)分別為節點i,j處振動加速度信號的頻譜;H(ω)為n×n的頻響函數矩陣,H(ω)=-(ω2M+iωC+K)-1;M為n×n的質量矩陣;C為n×n的阻尼矩陣;K為n×n的剛度矩陣。如果激振力只作用于一點k,式(1)變為

上述兩情形中,激振力僅提供了振動能量,并未參與傳遞率的計算。此時,振動傳遞率僅與激勵位置有關。根據式(1)和式(2),振動傳遞率與頻響函數有關,反映了結構的內在特性。因此,振動傳遞率可以用于結構狀態檢測。

采用三軸加速度計采集獲得結構三個正交方向的振動信號ax(t),ay(t),az(t),可以構建純虛四元數的時域序列aq(t)=[aq1,aq2,…,aqM],其中,aqr=iaxr+jayr+kazr,r∈[1,2,…,M],M為序列維度。此時,傳統的基于單通道信號頻譜之比的振動傳遞率不再適用。此時四元數傳遞率可定義為

式中:AQ(ω)為四元數譜;·為四元數點乘符號;?為共軛運算符號。

Ell提出了四元數傅立葉變換用于分析二維線性定常系統[16-17],其原始定義可表示為

式(4)定義了二維函數h(t,τ)從空間域(t,τ)向空間頻域(jω,kv)的轉換。四元數傅立葉變換最早用于彩色圖像的分析,具體形式[18]如下

式中:μ為隨機單位純虛四元數;f(m,n)為M×N彩色圖像矩陣(三顏色通道構建純虛四元數);(m,n)和(u,v)分別為空間域和空間頻域的坐標。對于純虛四元數振動時域序列aq(t),可視其為1×N的彩色圖像,因此,aq(t)的四元數譜計算公式為

由式(6)可知,AQ(ω)也是四元數序列,故由頻譜之比定義的四元數傳遞率也是在頻域內的四元數序列,和普通傳遞率一樣,也包含了全部的結構內在特性。

1.2 四元數傳遞率特性

參考實數信號的傅立葉變換,四元數時域序列亦可以由不同頻率分量的四元數信號的線性組合來表示。類似地,在四元數傅立葉變換中,可將變換后的模值作為縱坐標得到四元數幅頻譜。

不失一般性,構造如下三通道信號[x(t),y(t),z(t)],三個通道分別由不同頻率的正(余)弦信號組成

對x(t),y(t),z(t)三個單通道信號分別進行獨立的傅立葉變換,得到幅值譜如圖1(a)所示,然后將三個通道信號聯合成[x(t),y(t),z(t)],構建四元數序列并進行四元數傅立葉變換,可得到整體的頻域信息,其幅值譜如圖1(b)所示。對比可以發現,聯合信號的四元數頻譜包含了全部信號的頻率成分,是三通道信息的無損融合[20]。

圖1 四元數傅立葉變換Fig.1 Quaternion Fourier transform of three-channel joint signal

當振動方向偏離傳感器的敏感方向且僅使用單軸傳感器時,振動幅值會減小,但頻率不變。如果傳感器的安裝方向不能嚴格一致,傳統的基于單通道信號的振動傳遞率將偏離理論值。在實際工程應用中,很難保證檢測時傳感器的安裝方向嚴格一致。不僅如此,若振動方向是時變的,將直接導致傳遞率的時變,進而造成結構狀態的誤診斷。

圖2 (a)描述了任意向量v∈R3繞3D單位向量u旋轉2θ的情形。用上節中的定義將向量u和v擴展為純虛四元數,則旋轉之后得到的向量v′計算公式為

式中:euθ為單位四元數的指數形式,euθ=cosθ+usinθ。式(8)等價于三軸加速度計的測量坐標系原點不變,坐標系的方向繞向量u旋轉2θ,如圖2(b)所示。

圖2 四元數描述空間向量旋轉Fig.2 Quaternion describes vector rotation in space

向量v′的模值計算公式為

上述分析表明當采用三軸加速度計時,只要安裝位置固定,無論其安裝姿態如何變化,均可如實采集獲得空間振動信號。進而可推論,兩測點的三軸加速度計位置固定且不考慮噪聲干擾時,不論安裝姿態如何,式(3)的四元數傳遞率幅值均為其理論值。為證明上述推論,基于無砟軌道綜合試驗平臺,獲得兩測點振動信號的四元數時域序列,通過理論計算將測點A,B的信號分別繞u1,u2旋轉π/6和π/3,獲得三組四元數振動傳遞率,即初原始、僅旋轉A信號以及旋轉A/B信號的四元數振動傳遞率,如圖3所示。由圖3可知,當傳感器姿態發生變化時,四元數傳遞率的實部與虛部幅值均發生變化,但整體幅值不變。其中,u1,u2分別為

圖3 信號旋轉前后的四元數傳遞率幅值Fig.3 Amplitude of quaternion transmissibility before and after signal rotation

1.3 對偶樹四元數小波變換

對偶樹小波變換由離散小波變換(discrete wavelet transform,DWT)和復小波變換(complex wavelet transform,CWT)拓展而來,基于二維解析信號和希爾伯特變換(Hilbert transform,HT)構造而成。與實小波變換相比,對偶樹四元數小波變換可以獲得多角度、多尺度下的一個幅值和三個相位信息,同時具有良好的時移不變性[21]。

為了推廣到二維,Bülow[22]根據四元數傅里葉變換給出了四元數解析信號的定義,給定一個實數二維信號f(x,y),其四元數解析信號定義為

式中:fHi1(x,y)和fHi2(x,y)分別為沿x和y軸的希爾伯特變換;fHi(x,y)為完整的希爾伯特變換。

二維離散小波變換由一維離散小波變換沿兩個正交基的張量積構成,可表示為尺度函數φ(x)φ(y)、水平方向子帶函數φ(x)ψ(y)、豎直方向子帶函數ψ(x)φ(y)和對角線方向子帶函數ψ(x)ψ(y)的組合[23]。相似地,四元數小波同樣具有一個尺度函數和三個小波子帶函數,即將二維離散小波變換與三個實小波(由一維希爾伯特變換得到)組合構成一個四元數,從而得到二維解析小波。以對角子帶ψ(x)ψ(y)為例,其二維希爾伯特變換計算公式為

式中:(ψh,ψg)為特定濾波器h和g對應的小波函數,它遵循{ψh,ψg=HT(ψh)}的原則。式(12)中的每個分量均可由一維對偶樹復小波的組合計算。因此,對角線子帶四元數小波計算公式為

同理可得水平、垂直子帶四元數小波以及四元數尺度函數[24]分別為

對偶樹小波變換的計算細節可參考文獻[25]。

2 基于四元數聯合傳遞率的狀態檢測方法

2.1 卡洛變換(K-LT)法

K-LT是一種基于對象統計特征的最優正交變換。它可以消除統計相關性,但保留分類信息。將狀態已知條件下測量得到的四元數傳遞率模值向量作為訓練樣本,并構建為狀態矩陣SQ=[s1,s2,…,sN]∈RM×N,其中M為測量四元數傳遞率的長度,即參與運算的有效譜線數目,N為測量的傳遞率數目。狀態矩陣SQ的協方差矩陣CQ∈RM×N可以分解為

式中:W=[w1,w2,…,wN]∈RM×N為特征向量矩陣張成的特征子空間;Σ=diag(λ1,λ2,…,λN)∈RN×N為特征值矩陣,且λ1≥λ2≥…≥λN。對任意訓練四元數傳遞率模值向量sj,j=1,2,…,N,將其投影到特征子空間可獲得新的向量sj′=WT·sj。此時,s′j可作為狀態索引向量。將未知狀態下測得的四元數傳遞率模值向量以同樣方式投影到特征子空間,獲得投影向量s′t,然后計算s′j和s′t的歐氏距離,稱取得最小歐式距離的投影向量為匹配向量,即同一狀態。

2.2 對偶樹四元數小波變換(DQWT)法

如圖4所示,將同一狀態下4個測試四元數傳遞率模值向量疊加,構建狀態特征圖。圖4表明同一狀態下的不同樣本有很大的相似性而又有著局部的差異性。

圖4 四元數傳遞率灰度圖Fig.4 Grey-scale map of quaternion transmissibility

對狀態特征圖進行對偶樹小波變換,其系數可構成一個矩陣F,即

F的第一列包含圖像的低頻信息,第2~4列分別代表圖像在水平方向、垂直方向和對角線方向的高頻信息。

對偶樹小波變換的系數包含了一個幅值信息和三個相位信息。F的每一列可以構成式(17)的四元數,計算其每個四元數的幅值和相位,并重新組合得到幅值-相位矩陣

獲得每種已知狀態下的多組四元數傳遞率,構建狀態特征圖作為訓練樣本,對其進行對偶樹四元數小波變換得到低頻幅值相位向量AFi,同樣地,構建未知狀態下的狀態特征圖,計算其低頻幅值相位向量AFj,計算AFi和AFj之間的最小歐式距離并將其作為狀態識別指標,與上節方法類似。取低頻幅值相位信息,既綜合考慮了二維傳遞率圖像的幅值相位信息,又減小了單傳遞率的測量不確定性對檢測結果的影響,且降低了原始數據的大小。

3 狀態檢測實驗

實驗對象為總長20 m的CRTS-I型“無砟軌道綜合試驗平臺”,鋼軌兩端鎖緊,中間由配套的WJ-7型扣件鎖緊。實驗時在軌底隔跨兩組扣件安裝兩個三軸加速度傳感器PCB 356A33,采用力錘YD-5T在實驗區域左側某處的軌頂施加激勵,敲擊方向垂直于軌頂并保持激勵位置固定,以10 kHz的采樣頻率采集振動信號,采集設備為DHDAS-5902動態信號采集分析系統。實驗現場及傳感器布置如圖5所示。

圖5 實驗現場及傳感器布置Fig.5 Experiment site and sensor location

首先,利用鋼軌扣件的松脫與否模擬不同狀態,根據扣件松脫順序設置五種不同工況,分別為:扣件全部鎖緊,扣件1松脫,扣件1和2松脫,扣件1,2和3松脫以及4個扣件全部松脫。每種工況下,在左側激勵點用力錘激勵40次。根據前述四元數傳遞率函數的計算方法,計算得到每種工況下的40組傳遞率函數。圖6為工況1下的全部40次四元數傳遞率,取200~1700 Hz部分傳遞率函數,可見其重合性很好[26]。一方面,說明了傳遞率計算的正確性;另一方面,也驗證了四元數傳遞率函數僅由激勵位置和結構自身特性確定,而與激勵大小無關的結論。

圖6 工況1下40組原始四元數傳遞率函數幅值譜[26]Fig.6 Amplitude spectrums of 40 sets of original quaternion transmissibility under situation 1[26]

分別求取五種工況下的40組四元數傳遞率函數平均值,如圖7所示。在選取頻段內,傳遞率譜線大體重合,僅在部分頻率區間存在較明顯差別。這是因為相較實驗區域,試驗平臺在較長一段鋼軌上的扣件是鎖緊的,僅松脫局部相鄰扣件對整體結構的影響較小。

圖7 五種工況下平均四元數傳遞率函數幅值譜Fig.7 Average amplitude spectrum of quaternion transmissibility under five situations

采用2.1中基于K-LT的方法進行狀態識別,分別取五種工況的前30組傳遞率函數作為訓練樣本,后10組作為測試樣本,用來驗證該方法的有效性。檢測結果如表1所示,其中,測試QTJT(i)下的第一列數字表示與之取得最小歐氏距離的訓練樣本列序號,正確識別時第j測試組下的QTJT(i)對應序號應為30×(j-1)+1~30j;第二列表示狀態識別結果,‘√’表示識別正確,‘×’表示識別錯誤。由表1可知,該方法取得了98%的狀態識別正確率。

表1 基于K-LT方法扣件松脫狀態檢測結果[26]Tab.1 Detection results based on K-LT method

基于上述數據,采用2.2中DQWT方法進行狀態識別處理,分解層次取2,每種狀態下計算得到的40組QTJTs模值矩陣V∈R152×40,每4組QTJT作為一個樣本,每種狀態得到10組樣本,其中前6組作為訓練樣本,后4組作為測試樣本。識別結果如表2所示,狀態i(i=1,2,3,4,5)正確的匹配范圍應為6(i-1)+1~6i。由表2可知,該方法狀態識別正確率達到了100%。

表2 基于DQWT方法扣件松脫狀態檢測結果[26]Tab.2 Detection results based on DQWT method

采用LG-900型液壓鋼軌拉伸機拉伸鋼軌,模擬鋼軌內部縱向溫度力的變化,分別設置拉應力:0,10,20,30,40,50 MPa。首先,松脫全部扣件以釋放鋼軌內部縱向溫度力,穩定后拉伸機分別以上述設置應力動作進行鋼軌拉伸,而后再將全部扣件鎖緊,縮緊力矩為150 N·m。同樣地,每次拉伸之后按照前面的實驗方法,依次松脫相鄰兩組扣件,得到30組狀態,每種狀態采集20組振動信號,計算得到20組四元數傳遞率。在每種狀態下,以其前12組QTJT分別構建3個狀態特征圖作為訓練樣本,以后8組QTJT分別構建2個狀態特征圖作為測試樣本,表3為該實驗的狀態組織安排情況。基于DQWT方法進行狀態識別處理,每種狀態的匹配范圍應為3(i-1)+1~3i,所有樣本的識別結果如表4所示,該方法狀態識別正確率達到了100%。

表3 溫度力、扣件松脫綜合實驗結構狀態安排Tab.3 Structure state arrangement for comprehensive experiment of temperature force and fastener looseness

表4 的實驗結論表明,該方法可實現分辨力優于10 MPa的鋼軌溫度力變化識別。基于鐵路軌道相關標準,鋼軌的相關參數如下:截面積為7.745×10-3m2,線膨脹系數α=0.0118 mm/(m·℃),彈性模量E=2.1×105 MPa。則軌溫有1℃變化時所產生的內部縱向力計算公式為

表4 基于DQWT溫度力、扣件松脫綜合鋼軌狀態檢測結果[26]Tab.4 Detection results based on DQWT method for the comprehensive experiment

LG-900型液壓鋼軌拉伸機的油缸截面積為S=3.352×10-3m2,則其10 MPa的拉應力對應10 MPa×S=33520 N的拉力。進而,10 MPa的拉應力相對應的鋼軌溫度變化為33520 N/19192.11 N×1℃=1.75℃。因此,本方法可實現分辨力優于10 MPa的鋼軌溫度力變化識別,相當于鎖定狀態下1.75℃的軌溫變化。

4 結論

將傳統的基于單通道信號的振動傳遞率推廣到多通道信號的四元數振動傳遞率,并結合K-LT技術和DQWT,給出了基于四元數聯合傳遞率的的結構狀態檢測方法。通過在無砟軌道試驗平臺上的實際試驗,驗證了該方法的有效性,狀態識別結果正確率高于95%。上述方法提供了一種識別狀態變化是否存在的解決方案,未來可以通過增加測試點和依次識別相鄰四元數傳遞率,進一步發展為同時識別狀態變化的存在和位置。更重要的是,該方法保留了振動傳遞率的理論值,對傳感器的安裝方向沒有嚴格的要求,在實際應用中具有明顯優勢。特別是基于DQWT的方法,綜合考慮了多個被測傳遞率,降低了單次測量不確定性的影響。本文的研究為推動直接基于振動響應信號分析的狀態檢測方法發展起到了積極作用,具有重要技術應用價值。