基于導波能量特征PCA的結冰定量檢測方法研究

周世圓， 趙明華， 胡曉丹， 于全朋， 徐春廣

(北京理工大學 機械與車輛學院，北京 100081)

飛機在飛行過程中，當溫度較低時，飛機表面與云層發生碰撞，使得過冷水迅速發生相變，導致在飛機的許多部件上結冰[1]。飛機結冰會造成飛機質量增加，飛行阻力增大，影響飛機的氣動性能、操作穩定性和起降特性，嚴重降低飛機飛行的安全性，甚至導致飛行事故的發生[2]。因此，及時檢測和清除飛機表面上的結冰對于飛行安全非常重要。結冰檢測是除冰系統工作的基礎，可靠的結冰檢測對于后續的除冰是至關重要的。

針對飛機結冰檢測，研究人員提出了許多方法。Jarvinen[3]研究了通過測量阻抗、熱導率和溫度變化來檢測飛機上是否存在結冰并確定其厚度的方法。Zhuge等[4]提出了基于近紅外圖像處理技術的飛機機翼結冰檢測方法。Zou等[5]研究了使用傾斜端面光纖技術檢測結冰類型的方法。Ikiades等[6]研究了通過測量散射和反射光強度確定機翼結冰厚度和類型的方法。Liu等[7]應用超聲脈沖回波法測量頻率相關衰減系數，進而對結冰類型進行表征。Schlegl等[8]基于空氣、水和冰之間的電特性差異對飛機表面結冰進行檢測。然而，這些方法只能在一個點或一個小區域內探測結冰。對于飛機表面結冰，這些方法需要大量的分布式傳感器，實際上未在飛機結冰檢測中得到很好的應用。

近年來，超聲導波方法因其具有高靈敏度和大面積檢測能力而受到了許多研究人員的關注和研究。Vellekoop等[9]研究了基于Love波的結冰檢測方法，能夠測量1～1.5 mm厚度的結冰。Gao和Rose[10]將鋁層和冰層視為兩層結構，對SH波傳播進行建模，通過分析結冰的頻散曲線來檢測和分類結冰。Mlyniec等[11]研究表明鋁板中A0模態導波幅值和S0模態導波飛行時間可以用于結冰檢測。Zhao和Rose[12]通過結合Lamb波和概率重構算法，對結冰進行成像。趙偉偉等[13]基于橢圓定位算法研究了不同壓電片布局對導波結冰檢測成像質量的影響。

現有的導波檢測方法大多是對結冰進行定性檢測和分類，對結冰進行定量檢測的研究較少，因此很難在結冰初期探測到結冰，無法及時啟動除冰系統去除結冰，給飛機飛行帶來一定安全隱患。結冰定量檢測較難實現的主要原因是超聲導波具有多模態特性和頻散特性，結冰探測中導波信號波形較為復雜，給信號識別尤其是定量檢測帶來極大挑戰。飛機的機翼是常見的結冰部位，通常機翼蒙皮材料為鋁，為解決飛機結冰定量檢測難題，盡早發現飛機初期結冰，本文采用鋁板作為基底材料，采集結冰探測導波信號，提取導波信號多個特征，利用主成分分析(Principal Component Analysis,PCA)方法處理導波信號特征，獲得導波信號綜合得分，并建立綜合得分與結冰長度(結冰區域沿導波傳播方向尺寸)的關系，在結冰檢測時計算導波信號綜合得分，依據綜合得分與結冰長度的關系對鋁板上結冰長度進行定量檢測。通過仿真實驗和實測實驗驗證了該方法的可行性和可靠性。

1 結冰長度定量檢測PCA方法

PCA是一種通過降維技術把多個變量轉化為少數幾個主成分的統計方法，是基于降維思想產生的處理高維數據的方法[14]。PCA的主要目的是用較少的特征去解釋原來數據，將相關性很高的特征轉化成彼此相互獨立或不相關的特征[15]。主成分綜合得分可以作為表征數據的綜合性指標。

對于一組數據，共有n個樣本，每個樣本有p項評價指標，構成一個n×p階矩陣，

(1)

式中,xj=[x1j,x2j,…,xnj]T，j=1,2,…,p。PCA就是將p個觀測變量綜合成p個新的綜合變量，即主成分：

(2)

由于每個主成分的方差是遞減的，包含的信息量也是遞減的，所以實際分析時，一般不是選取p個主成分，而是根據各個主成分累積貢獻率的大小選取前m個主成分。根據學者Jolliffe的研究，一般累積貢獻率達到85%以上，可保證綜合變量包括原始變量的絕大部分信息[16]。要求累積貢獻率達到95%，以利于更好地保留原始變量的信息。

超聲導波檢測結冰時，由于導波信號復雜，采用多信號特征取代單信號特征進行檢測可提高檢測精度與可靠性。但由于信號各特征之間存在相關性，過多的特征會增加分析問題的難度與復雜性。因此考慮使用一個綜合特征對信號進行描述，從而得到結冰長度與該綜合特征的關系。PCA方法可以通過降維技術把多個變量轉化為少數幾個主成分，也容易根據主成分獲得綜合得分，從而獲得上述的“關系”。因此，采用PCA方法對導波信號特征進行降維，獲得信號特征主成分綜合得分，主成分綜合得分與結冰長度有對應關系，以此來對結冰長度進行定量檢測。

綜上所述，結冰長度定量檢測PCA的步驟流程如圖1所示。首先，提取導波信號特征，構建信號特征矩陣。由于各信號特征量級相差較大，采用Z標準化方法對數據進行標準化處理。Z標準化方法是最常用的標準化方法之一，將一組數據分別減去均值，除以標準差，標準化后的數據均值為0、方差為1[17]。然后計算標準化各特征之間的相關系數，構建相關系數矩陣，并通過雅可比方法求解矩陣特征值和特征向量。再根據累積貢獻率確定主成分個數，貢獻率為該主成分特征值占全部特征值總和的比例，累積貢獻率為貢獻率的累加和。確定主成分個數后，以特征值對應的特征向量為權值，對各標準化特征進行加權計算得到主成分，再以主成分貢獻率為權值，對各主成分進行加權計算得到綜合得分。最后，建立綜合得分與結冰長度的關系，依據此關系進行結冰檢測。

圖1 結冰長度定量檢測PCA步驟流程

2 結冰長度定量檢測仿真實驗

為驗證利用結冰探測導波信號主成分綜合得分測量結冰長度的可行性，對結冰鋁板結構進行有限元仿真，獲得多個結冰長度導波信號時域和時頻域特征，基于PCA方法對多信號特征進行分析，采用綜合得分對結冰長度進行評價。

2.1 仿真模型建立

在COMSOL中分別建立鋁板和結冰鋁板二維時域仿真模型，如圖2所示，選擇固體力學物理場進行分析。鋁板尺寸為160 mm×1 mm，材料選用COMSOL的內置材料Aluminum，楊氏模量為70 GPa，泊松比為0.33，密度為2700 kg/m3。結冰厚度為1 mm，長度分別為5 mm、10 mm、15 mm、20 mm、25 mm、30 mm、35 mm和40 mm，冰的類型為明冰[12]，楊氏模量為8.3 GPa，泊松比為0.35，密度為900 kg/m3。模型兩端面為低反射邊界，模型左端垂直于斜楔加載位移載荷，采用加漢寧窗的5個周期正弦信號激勵產生S0模態導波，中心頻率為2.5 MHz，在距離左端面100 mm處接收導波信號。網格剖分采用自由三角形網格，并且保證在1個波長范圍內至少包含8個網格。

圖2 鋁板和結冰鋁板仿真模型

2.2 仿真數據分析

在不同結冰長度下進行仿真，采集的導波信號如圖3所示。由圖3可知，信號幅值隨結冰長度的增加明顯減小，時域能量特征變化顯著，主要原因是導波從鋁板進入冰層，會發生反射、透射與折射，損失部分能量。因此，提取信號均方根值、標準差、絕對均值、方根幅值、峰峰值、波形因子、峰值因子、脈沖因子、裕度因子和峭度共10個與能量相關的時域特征。

結冰還會導致鋁板剛度、阻尼和結構固有頻率等發生變化，影響導波信號的頻率成分，使不同結冰長度下的各頻帶能量不同。為了分析結冰長度對導波信號各頻帶能量的影響，選擇db8小波基函數對導波信號進行3層小波包分解，得到各頻帶能量如圖4所示，從圖4中可以看出，頻帶2和頻帶4能量較高，其余頻帶能量接近于0，并且在結冰長度不同時，這2個頻帶能量變化較明顯。因此，提取這2個頻帶能量，作為時頻域特征。綜上，對導波信號共提取12個特征參數，如表1所示。

對提取的信號特征進行PCA，首先對數據進行標準化處理，用Zi(i=1,2,…,12)表示標準化后的12個特征參數，結果如表2所示。對標準化后的特征參數，計算兩兩之間的相關系數，構成相關系數矩陣，用雅可比方法求解相關系數矩陣的特征值和特征向量，并按λ1>λ2>…>λ12排列，計算主成分的貢獻率和累積貢獻率，如表3所示。主成分1和主成分2的貢獻率分別為89.515%和9.929%，累積貢獻率達到99.444%。根據累積貢獻率e≥95%能夠較好地解釋原始數據，確定2個主成分，并使用主成分1和主成分2計算綜合得分，對結冰長度進行評價。根據主成分特征值對應的特征向量得到主成分F1和F2表達式:

圖3 不同結冰長度仿真信號

F1=0.2984×Z1+0.2984×Z2+0.2664×Z3+0.2511×Z4+

0.3045×Z5+0.2868×Z6+0.2874×Z7+0.2877×Z8+

0.2880×Z9+0.2877×Z10+0.3030×Z11+0.2996×Z12

(3)

F2=0.1860×Z1+0.1860×Z2+0.4453×Z3+0.5168×Z4-

0.0513×Z5-0.2914×Z6-0.2721×Z7-0.3033×Z8-

0.3015×Z9-0.2914×Z10+0.1008×Z11+0.1549×Z12

(4)

根據主成分和貢獻率得到綜合得分F表達式:

F=0.89515×F1+0.09929×F2

(5)

表1 仿真信號特征參數

表2 標準化仿真信號特征參數

表3 主成分貢獻率和累積貢獻率(仿真)

經過計算得到綜合得分F隨結冰長度的變化關系如圖5所示。從圖5中可以看出，隨著結冰長度的增加，利用提取的結冰探測導波信號特征，基于PCA獲得的綜合得分逐漸減小。根據此關系曲線，通過計算未知長度結冰探測導波信號綜合得分，可以得到結冰長度。仿真實驗驗證了利用PCA方法處理結冰探測導波信號和定量檢測結冰長度的可行性。

圖5 綜合得分隨結冰長度變化關系

3 結冰長度定量檢測實驗驗證

為進一步驗證利用結冰探測導波信號主成分綜合得分測量結冰長度的可行性，搭建結冰探測實驗系統，開展不同長度結冰探測實驗，基于上述PCA方法對結冰探測導波信號進行分析，建立綜合得分與結冰長度的關系。進而對預制結冰區域開展結冰長度定量檢測驗證實驗，驗證PCA方法的有效性和可靠性。

3.1 實驗系統組成

實驗系統如圖6所示，包括導波檢測系統、結冰鋁板、換能器、斜楔和冰箱。被測鋁板尺寸為300 mm×300 mm×3 mm，采用斜入射換能器一發一收的形式進行結冰探測，激勵S0模態導波，激勵的中心頻率為2.5 MHz。在鋁板上設置7個結冰區域，如圖7所示，冰的厚度為1 mm、寬度為20 mm，長度分別為5 mm、10 mm、15 mm、20 mm、25 mm、30 mm、40 mm。

3.2 綜合得分與結冰長度關系建立

首先采集鋁板上無冰時的導波信號，然后采集鋁板上冰長分別為10 mm、20 mm、30 mm、40 mm時的導波信號。采集的5個導波信號如圖8所示，由圖8可知，信號的幅值隨結冰尺寸的增加明顯減小，時域能量特征變化顯著。根據上述信號特征提取方法，提取導波信號10個與信號能量相關的時域特征和2個與能量相關的時頻域特征，共12個特征參數，如表4所示。

圖6 結冰探測實驗系統

圖7 結冰區域示意圖

圖8 不同結冰長度實驗信號

用Zi(i=1,2,…,12)表示標準化后的12個特征參數，對數據進行標準化處理的結果如表5所示。計算標準化后的數據兩兩之間的相關系數，構成相關系數矩陣，并求解矩陣的特征值和特征向量，計算得到主成分1和主成分2的貢獻率分別為85.273%和11.192%，累積貢獻率達到96.465%，如表6所示。根據累積貢獻率e≥95%確定2個主成分，使用這2個主成分計算綜合得分，對結冰長度進行評價，主成分F1、F2及綜合得分F的表達式為

F1=0.3120×Z1+0.3102×Z2+0.2944×Z3+0.2734×Z4+

0.3122×Z5+0.3020×Z6+0.3123×Z7+0.3123×Z8+

0.3123×Z9-0.0485×Z10+0.2928×Z11+0.2750×Z12

(6)

表4 實驗信號特征參數

表5 標準化實驗信號特征參數

表6 主成分貢獻率和累積貢獻率(實驗)

F2=0.0485×Z1-0.0948×Z2+0.2573×Z3+0.3764×Z4+

0.0076×Z5-0.1944×Z6+0.0165×Z7-0.0121×Z8-

0.0211×Z9+0.8334×Z10-0.2180×Z11+0.0052×Z12

(7)

F=0.85273×F1+0.11192×F2

(8)

經過計算得到綜合得分F隨結冰長度的變化關系如圖9所示。從圖9中可以看出，隨著結冰長度的增加，基于PCA獲得的綜合得分逐漸減小，采用線性回歸進行擬合，表達式為

y=-0.1703x+3.4061

(9)

式中，y為綜合得分;x為結冰長度。擬合優度R2=0.9732，擬合效果較好，說明主成分綜合得分與結冰長度具有良好的線性關系，進一步驗證了利用導波信號主成分綜合得分測量結冰長度的可行性。在結冰檢測時，激勵并接收導波信號，提取所接收導波信號的能量特征，采用PCA法對導波信號能量特征進行處理并計算綜合得分，根據綜合得分與結冰長度的對應關系，可以獲得結冰長度。

圖9 綜合得分隨結冰長度變化關系

3.3 結冰長度定量檢測實驗驗證

為驗證基于導波能量特征PCA的結冰定量檢測方法的有效性，以預制的結冰長度分別為5 mm、15 mm和25 mm的結冰區域為例進行驗證，采集的結冰檢測導波信號如圖10所示。提取信號能量特征，并進行標準化處理，如表7所示。代入式(6)～式(8)計算得到其能量特征主成分綜合得分分別為2.5318、0.9762和-0.6532。根據綜合得分與結冰長度的關系式(9)，計算得到結冰長度分別為5.14 mm、14.27 mm和23.84 mm，與實際結冰長度誤差較小，證明了基于導波能量特征PCA的結冰長度定量檢測方法的有效性。

圖10 結冰長度定量檢測驗證實驗信號

表7 標準化驗證實驗信號特征參數

此外，采用傳統的信號能量分析法(信號幅值平方和)對上述鋁板上無冰和冰長分別為10 mm、20 mm、30 mm、40 mm時的導波信號進行了處理，得到導波信號能量隨結冰長度的變化關系如圖11所示。從圖11中可以看出，信號能量隨結冰長度的增大而減小，采用線性回歸進行擬合，表達式為,

y=-7420x+692669

(10)

式中，y為信號能量;x為結冰長度;擬合優度R2=0.9949。以上述預制的結冰長度分別為5 mm、15 mm和25 mm的結冰檢測導波信號進行驗證，計算得到信號能量分別為683583、542848和479732，代入式(10)計算得到結冰長度分別為1.22 mm、20.19 mm和28.70 mm，與實際結冰長度誤差較大。

圖11 信號能量隨結冰長度變化關系

兩種方法的檢測結果如圖12所示，由圖12可知，PCA法的檢測值更接近于實際值。檢測結果對比如表8所示，PCA法對預制5 mm、15 mm和25 mm冰長結冰的檢測誤差分別為2.8%、4.9%和4.6%，而傳統能量法的檢測誤差為75.6%、34.6%和14.8%。PCA法的檢測誤差較小，檢測精度明顯高于傳統的能量法。并且在結冰長度較小時，傳統能量法的檢測誤差非常大，基本無法進行定量檢測，而PCA法在結冰長度較小時仍有較高的檢測精度。因此，PCA法更適用于結冰長度的定量檢測，對飛機初期結冰有更好的檢測能力。

圖12 結冰長度定量檢測結果

表8 結冰長度定量檢測結果對比

4 結束語

提出了基于導波能量特征PCA的結冰定量檢測方法，通過仿真實驗對鋁板結冰進行了探測，提取了導波信號時域和時頻域中與能量相關的特征，利用PCA方法對導波信號特征進行降維，獲得了信號主成分綜合得分與結冰長度的關系，即隨著結冰長度的增加，信號特征的綜合得分逐漸減小，驗證了利用信號主成分綜合得分檢測結冰長度的可行性。并采用上述方法進行實測實驗，實驗結果表明主成分綜合得分與結冰長度具有良好的線性關系，并建立了綜合得分與結冰長度的變化關系式。對預制結冰區域的結冰長度定量檢測結果表明PCA法的檢測誤差較小，檢測精度高于傳統的能量法，尤其是在結冰長度較小時，PCA法的優勢更為突出。該方法對小尺寸結冰的探測能力較強，能夠用于飛機初期結冰的探測，從而在結冰初期啟動除冰系統，保證飛機飛行的安全性。

在結冰寬度和厚度一定的條件下實現了結冰長度的定量檢測，但未分析結冰寬度和厚度對導波信號綜合得分的影響，對結冰尺寸的檢測不完善。因此，下一步工作的重點是研究結冰寬度和厚度對導波信號綜合得分的影響，建立結冰尺寸與綜合得分的關系，實現結冰尺寸或體積的定量檢測。