基于EWT-FastICA的斜拉橋監(jiān)測撓度溫度效應(yīng)分離*

譚冬梅，姚歡，吳浩，甘沁霖

（1.武漢理工大學(xué)土木工程與建筑學(xué)院 武漢，430070）

（2.華中師范大學(xué)城市與環(huán)境科學(xué)學(xué)院 武漢，430079）

引言

隨著橋梁結(jié)構(gòu)健康監(jiān)測系統(tǒng)的發(fā)展，橋梁的健康狀態(tài)及安全評(píng)估越來越方便高效，但大量的數(shù)據(jù)積累造成很多有用的信息被掩蓋。撓度作為橋梁安全評(píng)估最直觀的參數(shù)之一，其中包括車載作用以及溫度作用引起的撓度，混凝土收縮徐變等引起的長期撓度［1］。將這些不同的效應(yīng)準(zhǔn)確分離出來有利于對(duì)橋梁進(jìn)行損傷識(shí)別及健康監(jiān)測，目前有許多學(xué)者針對(duì)這方面做過一系列的研究。

劉綱等［2］從不同時(shí)間尺度展開，以日溫差效應(yīng)頻率為中心頻率，利用粒子群優(yōu)化算法自適應(yīng)的調(diào)整濾波器的帶寬從而達(dá)到分離日溫差效應(yīng)。孫雅瓊等［3-4］通過應(yīng)變監(jiān)測數(shù)據(jù)，運(yùn)用時(shí)變的平均值法剔除出溫度效應(yīng)。陳國良等［5］根據(jù)撓度響應(yīng)在時(shí)間尺度上不耦合的特點(diǎn)，運(yùn)用中心移動(dòng)平均法從時(shí)間序列分析來剔除年溫差效應(yīng)和日溫差效應(yīng)，并用差分整合移動(dòng)平均自回歸模型（autoregressive integrated moving average model，簡稱ARIMA）預(yù)測了撓度變形的長期趨勢。劉夏平等［6］將溫度輸入最小二乘支持向量機(jī)中得到撓度溫度效應(yīng)來定量分析溫度與溫度效應(yīng)之間的關(guān)系，根據(jù)關(guān)系式來從溫度變化得到撓度溫度效應(yīng)的變化。梁宗保等［7］建立了應(yīng)變與溫度之間的經(jīng)驗(yàn)回歸方程，將溫度效應(yīng)視作回歸方程得到的值與應(yīng)變測量值之間的差值。譚冬梅等［8］通過采用改進(jìn)的集合經(jīng)驗(yàn)?zāi)B(tài)分解（modified ensemble empirical mode decomposition，簡稱MEEMD）算法進(jìn)行多次分解最終對(duì)溫度效應(yīng)進(jìn)行了分離。

筆者提出一種基于EWT-FastICA撓度成分分離算法。首先，采用經(jīng)驗(yàn)小波變化根據(jù)極大值選取原則，在頻譜上將日溫差效應(yīng)分割并構(gòu)造Meyer小波濾波器快速準(zhǔn)確分離日溫差效應(yīng)；其次，針對(duì)年溫差和長期撓度頻率相近難分離或者分離效果差的問題，修改傅里葉頻譜分割邊界，將年溫差和長期撓度成分充分分解，形成一系列IMF分量；然后，采用PCA降維，提取主元；最后，運(yùn)用FastICA算法解混分離出年溫差和長期撓度。EWT對(duì)比傳統(tǒng)的經(jīng)驗(yàn)?zāi)?態(tài) 分 解（empirical mode decomposition，簡 稱EMD）和集合經(jīng)驗(yàn)?zāi)B(tài)分解（ensemble empirical mode decomposition，簡稱EEMD）算法，可以極大地提升運(yùn)算速度，受噪聲影響小，自適應(yīng)性好，且結(jié)合FastICA算法解決了一般算法較難分離年溫差和長期撓度的問題。

1 理論簡介

1.1 經(jīng)驗(yàn)小波變換

經(jīng)驗(yàn)小波變換［9］是一種新穎的時(shí)頻分析變換方法，它是對(duì)傳統(tǒng)的小波變換以及EMD［10］的改進(jìn)。EWT實(shí)質(zhì)是通過自適應(yīng)分割頻譜，提取出原始信號(hào) 的 調(diào) 幅（amplitude modulation，簡稱AM）-調(diào)頻（frequency modulation，簡稱FM）成分，從而成功地分解出各IMF［11-12］。

EWT首先把信號(hào)的傅里葉譜歸一化得到頻率ωn∈(0，π)，再在歸一化后的區(qū)間上自適應(yīng)劃分將之分割成N份連續(xù)的小區(qū)間Λn，每個(gè)小區(qū)間長度為ωn-ωn-1，可以表示為

其中：ω0=0，ωN=π。

以每個(gè)分割界限ωn為中心設(shè)置一個(gè)寬度為2τn的過渡帶，τn如式2所示，自適應(yīng)分割后在各區(qū)間Λn構(gòu)造窄帶濾波器。根據(jù)傳統(tǒng)小波變換的定義，經(jīng)驗(yàn)小波變換的細(xì)節(jié)系數(shù)和近似系數(shù)可以表示為

因此，得到的各階IMF則可以表示為

1.2 主成分分析

主成分分析法（PCA）常運(yùn)用于數(shù)據(jù)降維，在不損失主要信息的情況下精簡了數(shù)據(jù)，極大提升了計(jì)算速度，同時(shí)可提取信號(hào)中的主元，降低信號(hào)中無用成分的干擾以及消除噪聲等。通過運(yùn)用n維向量代替m維向量，其中n＜m，PCA的處理步驟可以分為：

1）數(shù)據(jù)中心化；

2）計(jì)算信號(hào)的協(xié)方差矩陣C=E(XXT)；

3）特征值分解協(xié)方差矩陣；

4）將分解后的特征值λi按降序排列，選取占比較大的前幾個(gè)特征值；

5）根據(jù)步驟4所求特征值得到其對(duì)應(yīng)的特征向量ω1，ω2，…，ωn；

6）將選取的特征向量構(gòu)成新的矩陣乘上原信號(hào)矩陣實(shí)現(xiàn)降維。

1.3 快速獨(dú)立分量分析

快速獨(dú)立分量分析（FastICA）［13］是利用快速尋優(yōu)迭代的方法把采集到的盲信號(hào)進(jìn)行批量分解，找出統(tǒng)計(jì)獨(dú)立的源信號(hào)。盲信號(hào)指的是混合矩陣以及原始信號(hào)均無法知道的情況下，僅僅利用采集到的信號(hào)來估計(jì)源信號(hào)。

設(shè)有m個(gè)源信號(hào)Sn×m=[s1(t)，s2(t)，…，sm(t)]存在某個(gè)混合矩陣An×n，使得

其中：X為觀測信號(hào)。

在任一時(shí)刻t，X均符合式（7）。因此，關(guān)鍵目標(biāo)則是找出分離矩陣W，其中W=A-1使得Y(t)=為S(t)的估計(jì)。

快速獨(dú)立分量分析可以分為如下步驟：

1）將觀測信號(hào)X去均值；

2）進(jìn)行白化預(yù)處理；

3）引入定點(diǎn)迭代方式。

定點(diǎn)迭代方式為

其中：E｛｝為特征值分解；v為白化處理后的數(shù)據(jù)；g（y）為非線性二次函數(shù)的3種形式；a1為常數(shù)，通常取為1。

通過近似牛頓法選擇系數(shù)α，對(duì)式（8）兩邊各自加上αw；根據(jù)Kuhn-Tucker條件，在約束條件下，E{G(wTv)}在E{vg(wTv)}+βw=0點(diǎn)取得最優(yōu)。令K=E{vg(wTv)}+βw，其中β為常數(shù)。

此時(shí)式（10）變成了可逆的對(duì)角矩陣，于是得到如下牛頓近似迭代

兩邊同乘上E{g′(wTv)+β}，最終得到

FastICA雖然能很好地解決頻率混疊的問題，但是也有著幅值不確定性的缺點(diǎn)。當(dāng)混合矩陣w某行乘上某個(gè)常數(shù)k，對(duì)應(yīng)位置源信號(hào)同時(shí)乘上常數(shù)1/k，結(jié)果將會(huì)不變。

1.4 EWT-FastICA分離流程

首先，利用EWT將橋梁撓度信號(hào)分解成為一系列IMF，其中由于日溫差效應(yīng)相比年溫差效應(yīng)以及長期撓度頻率有明顯的差異，便可根據(jù)頻率有效地剔除出其中的日溫差效應(yīng)；其次，把剔除日溫差效應(yīng)后剩余的IMF運(yùn)用PCA降維，根據(jù)特征值貢獻(xiàn)占比來選取其對(duì)應(yīng)的特征向量，將選取的特征向量與原始矩陣重構(gòu)得到降維后矩陣；最后，把PCA降維后的矩陣運(yùn)用FastICA盲源分離算法分解出長期撓度與年溫差效應(yīng)。分離流程如圖1所示。

圖1 分離流程圖Fig.1 Separation flow chart

2 橋梁仿真撓度信號(hào)分離

2.1 橋梁模型仿真信號(hào)

橋梁撓度作為多種因素共同作用的結(jié)果，其長期監(jiān)測的數(shù)據(jù)中包括車載和風(fēng)載等引起的動(dòng)撓度、溫度作用撓度以及混凝土收縮徐變等引起的長期撓度［14］，其中溫度作用引起的撓度主要包括日溫差效應(yīng)和年溫差效應(yīng)。運(yùn)用Midas civil軟件建立了如圖2所示的仿真模型，該模型為武漢某斜拉橋的簡要模型。

圖2 武漢某斜拉橋的模型圖Fig.2 Model drawing of a cable-stayed bridge in Wuhan

根據(jù)改變模型溫度來觀測主跨跨中的撓度變化，其中每升高溫度1℃，跨中下?lián)?.53 mm；每降低溫度1℃，跨中上撓1.53 mm。再對(duì)跨中截面升高1℃，跨中下?lián)?.34 mm；截面降低1℃，跨中上撓0.34 mm。根據(jù)武漢天氣，假設(shè)每天平均日溫差為12℃，梁橫截面日溫差可假設(shè)為6℃，武漢1年的年溫差可以假設(shè)為36℃。日溫度變化以及年溫度變化均可以用正弦函數(shù)來表示，由于溫度與溫度效應(yīng)之間簡化為簡單的線性關(guān)系［15］，所以溫度效應(yīng)亦可表示成正弦函數(shù)。

整體日溫差效應(yīng)f11=-9.18sin(πt/12)，截面日溫差效應(yīng)f12=-1.02sin(πt/12)，因此日溫差效應(yīng)f1(t)=f11(t)+f12(t)，年溫差效應(yīng)f2(t)=-27.54sin(πt/4 380)。長期撓度f3(t)根據(jù)《公路鋼筋混凝土及預(yù)應(yīng)力混凝土橋涵設(shè)計(jì)規(guī)范》（JTG 3362—2018），運(yùn)用指數(shù)型函數(shù)擬合得到，則總溫度效應(yīng)以及長期撓度可表示為

年溫差1個(gè)周期為8 760 h。為了方便顯示日溫差曲線圖，日溫差只取前4 kh的數(shù)據(jù)，年溫差曲線圖取2個(gè)周期17 520 h的數(shù)據(jù)。各個(gè)效應(yīng)時(shí)域曲線及頻域圖如圖3～5所示。

圖3 日溫差效應(yīng)時(shí)域曲線及頻譜圖Fig.3 Time domain curve and frequency spectrogram of daily temperature difference effect

圖4 年溫差效應(yīng)時(shí)域曲線及頻譜圖Fig.4 Time domain curve and frequency spectrogram of annual temperature difference effect

圖5 長期撓度時(shí)域曲線及頻譜圖Fig.5 Time domain curve and frequency spectrogram of long-term deflection

2.2 溫度效應(yīng)分離

各效應(yīng)合成的總撓度時(shí)程曲線如圖6所示。由于年溫差效應(yīng)的頻率為1/8 760 Hz及長期撓度頻率過低接近于0，因此一般的分離方法較難分離出年溫差效應(yīng)和長期撓度。

圖6 總撓度時(shí)域曲線Fig.6 Time domain curve of total deflection

采用EWT對(duì)合成撓度信號(hào)進(jìn)行分離，首先，選用Locmax分割策略，在頻譜上自動(dòng)判斷選出極大值；其次，在兩相鄰極大值中間分割，對(duì)每個(gè)極大值建立Meyer小波濾波器。當(dāng)取N值為3時(shí)，由于0和π總是第1個(gè)和最后1個(gè)邊界，此時(shí)在信號(hào)頻譜上根據(jù)極大值還會(huì)再次產(chǎn)生2個(gè)分割邊界。EWT頻譜切割圖和分割后接近0處的局部放大圖分別如圖7和圖8所示，表明日溫差頻率部分得到較好的切割，但年溫差與長期撓度頻率接近，切割效果差。圖9為分解后的IMF圖，表明對(duì)日溫差進(jìn)行了有效的分離（為顯示方便，僅顯示前4 kh），但長期撓度與年溫差存在混疊。

圖7 EWT頻譜切割圖Fig.7 EWT segmentation of spectrogram

圖8 EWT分割后接近0處的局部放大圖Fig.8 Enlarged local image near 0 after EWT segmentation

圖9 EWT分解后IMF圖Fig.9 IMF diagram after EWT decomposition

即使選用Locmaxmin分割策略，分割界限位于年溫差和長期撓度交匯處，分離得到的年溫差與長期撓度也嚴(yán)重混疊。因此，筆者提出在日溫差分割界限與年溫差長期撓度界限間設(shè)定多個(gè)分割界限，將年溫差和長期撓度混疊部分分離為多個(gè)IMF，取自定分割區(qū)間數(shù)為7，分解如圖10所示。

圖10中IMF7便是分離出來的日溫差效應(yīng)，由IMF1～I(xiàn)MF6可看出長期撓度與年溫差混疊。從原始信號(hào)f(t)中剔除掉IMF7，得到長期撓度與年溫差效應(yīng)的混合信號(hào)S1(t)，將前6階IMF以及S1(t)組成多通道的混合信號(hào)，再運(yùn)用PCA算法降維提取主元。PCA數(shù)據(jù)特征值處理結(jié)果如表1所示，其中占比為該特征值與所有特征值之和的比值。由表1可知，前2階特征值遠(yuǎn)大于其余特征值，占比較大，于是選取該2階特征值所對(duì)應(yīng)的特征向量來確定降維后的數(shù)據(jù)，再把降維后的數(shù)據(jù)運(yùn)用FastICA算法盲源分離，然后對(duì)分離后的結(jié)果與原始信號(hào)進(jìn)行幅值對(duì)比，把比值作為系數(shù)乘上分離后的結(jié)果，最終分離出日溫差、年溫差和長期撓度，各分離結(jié)果如圖11所示。

圖10 EWT分解圖Fig.10 EWT decomposition diagram

表1 PCA處理后特征值數(shù)據(jù)Tab.1 Eigenvalue data after PCA processing

圖11 最終分解圖Fig.11 Final decomposition diagram

通過相關(guān)系數(shù)（r）來評(píng)價(jià)分解效果，其中相關(guān)系數(shù)式為

其中：x為實(shí)際撓度信號(hào)；y為分離后的信號(hào)。

當(dāng)相關(guān)系數(shù)r越接近于1，則表明分離的效果越好。將分離前后的信號(hào)代入式（15）計(jì)算，結(jié)果如表2所示。

表2 各效應(yīng)分離前后的相關(guān)系數(shù)Tab.2 Correlation coefficient before and after effect separation

從表2可以看出，日溫差效應(yīng)以及年溫差效應(yīng)的相關(guān)系數(shù)都接近于1，均實(shí)現(xiàn)了很好的分離，但其中年溫差端點(diǎn)效應(yīng)較為明顯；長期撓度由于受年溫差效應(yīng)干擾，分離效果相對(duì)較差，但總體相關(guān)系數(shù)都較高，證明了該方法對(duì)于各效應(yīng)分離的有效性。

3 實(shí)測橋梁撓度信號(hào)各成分分離

3.1 實(shí)際工程背景

武漢某斜拉橋主橋?yàn)殡p塔雙索面鋼箱梁與預(yù)應(yīng)力混凝土箱梁組合型斜拉橋結(jié)構(gòu)，如圖12所示。該斜拉橋采用北斗實(shí)時(shí)在線監(jiān)測技術(shù)，在斜拉橋關(guān)鍵位置布置撓度測點(diǎn)，用來監(jiān)測橋梁運(yùn)營狀態(tài)的安全情況，達(dá)到實(shí)時(shí)監(jiān)控?fù)隙茸兓男Ч１倍吠ㄟ^全自動(dòng)化的高精度監(jiān)測技術(shù)對(duì)斜拉橋進(jìn)行了全年撓度連續(xù)的監(jiān)測，每秒記錄1次撓度值。當(dāng)撓度變化大于閾值時(shí)，將會(huì)自動(dòng)產(chǎn)生告警，能為橋梁后期的運(yùn)營服務(wù)提供很好的幫助。測點(diǎn)位置布置如圖13所示。

圖12 武漢某斜拉橋Fig.12 A cable-stayed bridge in Wuhan

圖13 武漢某斜拉橋主橋監(jiān)測測點(diǎn)布置圖(單位：m)Fig.13 Monitoring map of main bridge of a cable-stayed bridge in Wuhan(unit:m)

3.2 實(shí)際橋梁撓度分離

選取主跨跨中測點(diǎn)BD12的撓度數(shù)據(jù)，時(shí)間間隔為2017年9月 至2018年8月，每秒采集數(shù)據(jù)1次，每小時(shí)取樣1次，BD12撓度曲線如圖14所示。

圖14 BD12全年撓度數(shù)據(jù)Fig.14 BD12 annual deflection data

首先，對(duì)實(shí)測撓度數(shù)據(jù)用EWT分解；其次，將分解的低頻信號(hào)進(jìn)行PCA降維，根據(jù)特征值占比大小選取所對(duì)應(yīng)的特征向量重組信號(hào)；最終，將重組后的信號(hào)輸入到FastICA算法，分離出實(shí)測撓度中的各效應(yīng)作用，并且考慮到實(shí)際分離信號(hào)受到環(huán)境因素等的影響，將分離結(jié)果擬合，得到BD12全年撓度數(shù)據(jù)分離結(jié)果如圖15所示。

圖15 BD12全年撓度數(shù)據(jù)分離結(jié)果Fig.15 Separation results of BD12 annual deflection data

由于相鄰測點(diǎn)之間受溫度效應(yīng)的影響基本相同，為了驗(yàn)證實(shí)測撓度信號(hào)分離的效果，取該測點(diǎn)對(duì)稱處跨中BD35測點(diǎn)的數(shù)據(jù)進(jìn)行分離，將各自分離效果進(jìn)行對(duì)比，圖16和圖17分別為BD35的全年撓度時(shí)程和分離后各撓度效應(yīng)。

圖16 BD35全年撓度數(shù)據(jù)Fig.16 BD35 annual deflection data

圖17 BD35全年撓度數(shù)據(jù)分離結(jié)果Fig.17 Separation results of BD35 annual deflection data

BD12和BD35測點(diǎn)分離效果對(duì)比如圖18所示，從圖中可以看出，BD12與BD35分解的效果相似度較好，證明準(zhǔn)確率較高，分離效果好。將此方法與EEMD結(jié)合改進(jìn)PCA算法進(jìn)行對(duì)比，兩側(cè)點(diǎn)的分離結(jié)果相關(guān)系數(shù)如表3所示。

圖18 對(duì)稱測點(diǎn)分離結(jié)果對(duì)比圖Fig.18 Comparison chart of separation results of symmetrical measuring points

由表3可知，兩測點(diǎn)相關(guān)系數(shù)均高于0.9，說明基于EWT-FastICA與EEMD和改進(jìn)PCA方法都能較好地分離出實(shí)際橋梁撓度信號(hào)中的溫度效應(yīng)，但基于EWT-FastICA算法精度較高。同時(shí)由于EWT是直接通過在頻譜上自適應(yīng)切割并構(gòu)造濾波器，而EEMD則是通過多次添加高斯白噪聲和多次求平均，因此EWT相比EEMD計(jì)算速度具有較大提升，可顯著提高分離效率。為對(duì)比計(jì)算速度，選擇2017年9月至2018年9月實(shí)測數(shù)據(jù)，考慮溫度效應(yīng)周期較長特性，故每小時(shí)提取1個(gè)數(shù)據(jù)，共8 928個(gè)數(shù)據(jù)分別運(yùn)用EWT和EEMD算法進(jìn)行信號(hào)分離，運(yùn)行速度結(jié)果如表4所示。從表4可以看出，EWT運(yùn)行速度遠(yuǎn)高于EEMD，而橋梁監(jiān)測海量數(shù)據(jù)實(shí)時(shí)在線分析中，實(shí)測橋梁采集頻率通常達(dá)到幾十赫茲甚至幾百赫茲，將會(huì)產(chǎn)生龐大的數(shù)據(jù)量，故EWT算法體現(xiàn)出了更大的優(yōu)越性。

表3 兩測點(diǎn)分離后的相關(guān)系數(shù)Tab.3 Correlation coefficient before and after effect separation

表4 運(yùn)行速度對(duì)比表Tab.4 Running speed comparison table

4 結(jié)論

1）采用EWT算法通過在頻譜上自適應(yīng)分割并構(gòu)造小波濾波器將單通道的混合信號(hào)分解成多通道信號(hào)，解決了FastICA算法需要觀測信號(hào)數(shù)大于等于源信號(hào)數(shù)的先決條件。

2）EWT算法相對(duì)于EMD和EEMD算法做到了直接分離日溫差效應(yīng)，對(duì)于難以分離的長期撓度和年溫差，則通過PCA降維提取主元向量后，再基于FastICA算法較好地解決了頻率相近信號(hào)模態(tài)混疊的問題，成功分離出長期撓度與年溫差效應(yīng)。

3）從模擬信號(hào)和實(shí)測橋梁撓度信號(hào)的相關(guān)系數(shù)以及EWT運(yùn)行速度來看，EWT算法較EEMD算法極大增強(qiáng)了運(yùn)算速度，對(duì)于橋梁實(shí)時(shí)在線監(jiān)測海量數(shù)據(jù)分析，EWT算法體現(xiàn)出較大的優(yōu)越性。