基于改進果蠅優(yōu)化算法正交匹配追蹤的超聲信號降噪方法

董 明，李 敬，索永錄，唐恩賢，馬宏偉，陳 淵，張廣明，萬 翔

（1. 西安科技大學(xué)機械學(xué)院，陜西西安 710054；2. 陜西陜煤黃陵礦業(yè)集團有限責任公司技術(shù)中心，陜西延安 727307；3. 西安科技大學(xué)能源學(xué)院，陜西西安 710054；4. 陜西省礦山機電裝備智能監(jiān)測重點實驗室，陜西西安 710054；5. 英國利物浦約翰摩爾斯大學(xué)通用工程研究所，英國利物浦L3 3AF）

1 引言

超聲檢測具有操作方便、穿透能力強、靈敏度高等優(yōu)點，被廣泛應(yīng)用于材料內(nèi)部缺陷檢測和材料特性表征. 探頭接收到的回波微弱且波形復(fù)雜，易受到噪聲的干擾，嚴重時回波完全淹沒在噪聲中，影響評價的準確性［1］. 超聲波在材料晶界上發(fā)生折射、反射和波形轉(zhuǎn)換，引起雜亂回波，晶界散射與波長和晶粒尺寸有關(guān)，這類噪聲被稱為結(jié)構(gòu)噪聲；超聲波還受環(huán)境中電磁干擾以及測量系統(tǒng)產(chǎn)生的電子噪聲的影響［2，3］，這類噪聲通常是一種不相干的高斯白噪聲. 超聲信號降噪技術(shù)是準確提取缺陷信息的關(guān)鍵，降噪效果直接影響檢測結(jié)果的表征. 降噪通常是將信號在特定的基函數(shù)上進行分解，例如小波變換和傅里葉變換等，基函數(shù)的特征決定了信號的特征［4，5］. 超聲信號具有非平穩(wěn)性和非線性的特點，在單一的、有限的基函數(shù)上進行分解，得到的分量依賴基函數(shù)，信號自身的特征被忽略，容易引起重構(gòu)信號失真. 稀疏分解（sparse decomposition）將信號在過完備原子庫中分解，提高了降噪效果. 近幾年提出的稀疏分解算法有基追蹤、稀疏貝葉斯學(xué)習(xí)、Focuss、匹配追蹤（Matching Pursuits，MP）等［6～9］，其中MP 是一種“貪婪”的算法，在整個原子字典中遍歷式搜索最佳匹配原子，其計算效率較差，如果殘差與最佳原子非正交，則結(jié)果不是最優(yōu)的. 為此提出了正交匹配追蹤（Orthogonal Matching Pursuit，OMP），對每次匹配到的最佳原子進行正交化處理，使得殘差與選擇過的原子都是正交的，與MP 算法相比，OMP 算法每次選擇的都是最佳原子，但是OMP算法引入了正交化，每一次迭代的計算量較大，計算時間較長［10］. 朱會杰等人［11］用子空間追蹤法選擇最佳原子，并提出一種自適應(yīng)迭代停止標準，結(jié)果表明該方法在信噪比和均方誤差方面都優(yōu)于常規(guī)去噪算法. 魏東等人［12］提出了一種K-SVD（K-Singular Value Decomposition）和OMP 相結(jié)合的超聲信號降噪算法，利用K-SVD 算法訓(xùn)練反映信號特征的過完備字典，提高了重構(gòu)精度，但是訓(xùn)練字典增加了計算量. 采用群搜索算法可以將離散的有限的字典原子變?yōu)樵谶B續(xù)空間中無限的原子，在連續(xù)空間進行分解，提高了回波重構(gòu)精度.

果蠅優(yōu)化算法（Fruit Fly Optimization Algorithm，F(xiàn)OA）是一種基于果蠅覓食行為的仿生算法，相較于人工魚群算法［13］、人工蜂群算法［14］、粒子群算法［15］等群搜索算法，F(xiàn)OA 算法收斂速度快、全局搜索能力強、易于學(xué)習(xí)，但是其性能依賴搜索半徑，不易收斂到最優(yōu)值，全局搜索能力和局部搜索能力存在矛盾［16］. 石建平等人［17］提出2種策略隨機選擇的優(yōu)化算法，避免了算法過早收斂. 鳳麗州等人［18］引入細菌趨化理論，根據(jù)果蠅分布特點使用雙重驅(qū)動更新果蠅位置，避免無效搜索.王友衛(wèi)等人［19］將果蠅種群劃分為搜索果蠅和跟隨果蠅，并使用分區(qū)采樣的策略，提高了收斂的穩(wěn)定性.Wu等人［20］提出了搜索方向的自適應(yīng)選擇機制，果蠅以較大的概率飛向最佳搜索方向，以小概率朝遠離最佳搜索方向飛行.Zhang等人［21］提出了一種新的多尺度協(xié)同變異FOA算法，克服了局部最優(yōu)的局限性.Fu等人［22］提出了一種結(jié)合隨機模擬的多目標離散FOA，從解的表示、啟發(fā)式解碼規(guī)則、嗅覺搜索、視覺搜索和遺傳搜索等方面進行改進以增強搜索能力. 田旭等人［23］針對多維優(yōu)化問題，每一次搜索僅擾動其中一維，加快搜索速度.

本文針對正交匹配追蹤算法分解速度慢、重構(gòu)精度低等問題，提出了基于改進果蠅優(yōu)化算法和正交匹配追蹤（Improved Fruit Fly Optimization Algorithm based Orthogonal Matching Pursuit，IFOA-OMP）的超聲回波降噪方法，利用一種自適應(yīng)步長計算方法優(yōu)化全局遍歷性，采用高維廣義CAT 映射使算法能夠跳出局部最優(yōu)值，對仿真和實驗超聲信號進行處理，驗證了本文方法的降噪效果.

2 正交匹配追蹤算法

2.1 Gabor字典

Gabor 函數(shù)與超聲信號的特征最為接近，處理超聲信號效果較好，Gabor 字典是由Gabor 函數(shù)經(jīng)過離散化產(chǎn)生.Gabor函數(shù)為

其中，g(t)=e-πt2是高斯窗函數(shù)；s代表尺度影響參數(shù)，確定函數(shù)的能量分布情況；u為平移參數(shù)，確定回波在橫坐標上的位置；v為頻率參數(shù)，確定函數(shù)的中心頻率；w為相位參數(shù)，確定函數(shù)的初始相位.

分別將Gabor函數(shù)中的s，u，v和w這4個參數(shù)離散化，構(gòu)建過完備Gabor字典庫{gγ}γ∈Γ，即

其中，α=2，Δu=1/2，Δv=π，Δw=π/6，0 ＜j≤log 2N，0 ＜p＜N·2-j+1，0 ≤k＜2j+1，0 ≤i≤12；N為信號的采樣點數(shù).

按照上述準則對4 個參數(shù)進行離散，生成多個Gabor原子，構(gòu)成過完備原子庫.

2.2 正交匹配追蹤

正交匹配追蹤（OMP）在匹配追蹤（MP）的基礎(chǔ)上，引入了正交化的步驟，通過迭代計算從Gabor字典中找出與原始信號匹配度最高的原子，并對每一次匹配到的最佳原子進行正交化處理，避免重復(fù)選擇同一原子.OMP的具體流程如下.

（1）輸入初始信號f，初始化殘差R0=f，將殘差與過完備原子庫{gγ}γ∈Γ中的原子gγ求內(nèi)積，保留第一次迭代找到的內(nèi)積最大的原子gγ1，即

令gu1=gγ1，并對gu1進行歸一化處理，gb1=，第1次迭代的殘差R1=R0-＜R0，gb1＞gb1.

（2）進行迭代計算，第k次（k≥2）迭代得到的最佳原子gγk為

（3）對最佳原子gγk進行施密特正交化處理，即

對guk進行歸一化處理，gbk=.

（4）將殘差Rk-1分解為在最佳原子上的投影＜Rk-1，gbk＞gbk，更新殘差Rk，即

（5）判斷是否滿足迭代停止條件，若不滿足，則重復(fù)步驟（2）～（5），若滿足條件，則保存最佳原子，并重構(gòu)回波.

OMP 算法的收斂速度較快，但是由于引入了正交化，提高了算法的復(fù)雜度，增加了運行時間.

3 改進的果蠅優(yōu)化算法

3.1 果蠅優(yōu)化算法

果蠅優(yōu)化算法（FOA）源于果蠅的覓食行為，是一種全局優(yōu)化方法. 果蠅依靠敏銳的嗅覺發(fā)現(xiàn)遠處的食物，當靠近食物時，依靠視覺找到食物和同伴聚集的位置.FOA算法的基本流程如下.

（1）確定果蠅種群的迭代次數(shù)Maxgen、種群大小Sizepop 和搜索范圍LR，設(shè)果蠅種群的起始位置(X0，Y0)為

其中，rand(·)是產(chǎn)生［0，1］隨機數(shù)的函數(shù).

（2）果蠅開始尋找食物，給定隨機的搜尋方向和距離，得到新的目標位置，即

其中，randvalue1和randvalue2是果蠅搜索下一目標的隨機步長.

（3）計算每個果蠅個體到原點的距離Disti和味道濃度判別值Si為

（4）將味道濃度判別值Si代入濃度判別函數(shù)Function(·)，得到每個果蠅所在位置的味道濃度，即

（5）保存果蠅群體里味道濃度最高的個體此時所在的位置和味道濃度的大小，即

（6）保存最優(yōu)味道濃度smellbest=bestsmell，群體中的其他果蠅朝著味道濃度最高的位置飛去，更新下一次搜索的位置，即

（7）重復(fù)步驟（2）～（5），若搜尋到的味道bestsmell濃度大于之前迭代所得到的味道濃度smellbest，并且迭代次數(shù)小于最大迭代次數(shù)，則執(zhí)行步驟（6），反之，則結(jié)束算法.

上述FOA 算法是針對單參數(shù)的優(yōu)化問題，多參數(shù)時需要多個果蠅種群，每個果蠅種群搜索一個參數(shù)，多種群同時搜索解決多參數(shù)優(yōu)化問題. 對于超聲信號的降噪，需要優(yōu)化Gabor 函數(shù)中的s，u，v，w這4 個參數(shù)，因此需要4個果蠅種群. 在FOA 算法中，提高最佳原子的匹配度需要提高算法的局部搜索能力，提高全局遍歷性又要求算法能及時跳出局部最優(yōu)，這兩者是互相矛盾的，為了搜尋到最優(yōu)的原子，本文從以下兩方面對FOA算法進行了改進.

3.2 自適應(yīng)步長

果蠅搜索過程中，每迭代一次，果蠅就會按照隨機的方向和步長飛向下一個目標，即

其中，step1=step2，并且step1和step2的值是預(yù)先給定且不變的. 第i個果蠅的位置為Pi=(Xi，Yi)，此時該果蠅距離群體的距離. 通過改變X0，Y0，randvalue1和randvalue2的值可以調(diào)節(jié)果蠅和種群之間的距離，存在以下2種極端情況：

（1）若X0/randvalue1=0 和Y0/randvalue2=0，此時第i個果蠅距離群體無窮遠，代表果蠅可以在無窮遠處搜索；

（2）若|X0/randvalue1|→∞和|Y0/randvalue2|→∞，此時第i個果蠅距離群體無窮近，代表果蠅可以在無窮近處搜索.

因此果蠅的搜索半徑范圍是(0，∞). 搜索步長的選擇直接影響函數(shù)參數(shù)估計的準確性，在FOA 算法迭代尋優(yōu)的過程中，前期需要較大的步長使得果蠅能夠迅速收斂到最優(yōu)值附近，后期需要較小的步長來準確搜索到最優(yōu)值，適當調(diào)整搜索步長，可以保證果蠅在全局和局部都具有良好的搜索能力. 在FOA 算法的應(yīng)用中通常采用隨機步長，可能在迭代前期步長太小，導(dǎo)致收斂速度慢，在后期步長太長導(dǎo)致不能收斂到最優(yōu)值，不利于參數(shù)優(yōu)化. 為了解決這一問題，本文提出了自適應(yīng)步長，即

其中，W為初始步長；T為步長收斂系數(shù)，且T=1/λδ（λ為正交匹配追蹤算法的迭代次數(shù)，δ為步長收斂因子，δ取0.1～0.25）；γ為非均勻變異因子，γ取0.6；sign(·)是符號函數(shù).

取初始步長W=20，收斂系數(shù)δ=0.2，改進后的搜索步長如圖1 所示，可見，步長隨著種群迭代而不斷減小，滿足果蠅遍歷搜索的要求.

圖1 自適應(yīng)步長的演化過程

為了對比自適應(yīng)步長和隨機步長的搜索性能，分別對加入信噪比為-1 dB 高斯白噪聲的超聲信號進行處理，OMP 每一次搜索都會找到一個最佳原子，F(xiàn)OA 可以優(yōu)化OMP 的搜索原子過程，以均方根誤差（Root Mean Square Error，RMSE）表示處理的效果.RMSE反映了重構(gòu)信號和原始信號之間的平均誤差，輸出信噪比越大，RMSE 越小則說明降噪的效果越好. 如圖2所示，采用隨機步長經(jīng)過3 次搜索，重構(gòu)信號的RMSE 達到最小，而用自適應(yīng)步長需要經(jīng)過4 次搜索才能達到最小RMSE，收斂速度略慢，但相比于隨機步長，重構(gòu)信號的RMSE 更小，有效改善了FOA 算法搜索的精度. 2 種方法都存在過匹配的問題，當找到最優(yōu)值以后，繼續(xù)迭代反而會造成重構(gòu)信號失真，可以通過設(shè)置迭代停止條件解決此問題.

圖2 不同原子搜索步長性能對比

3.3 高維廣義CAT映射

FOA 算法與大多數(shù)群搜索算法一樣，都存在陷入局部最優(yōu)的問題. 為了提高算法跳出局部最優(yōu)的能力，引入跳脫參數(shù)c，設(shè)置跳脫參數(shù)閾值為f，若最佳味道濃度更新，則將c置為0. 跳脫參數(shù)記錄已選擇的最佳味道濃度未被替換的次數(shù)，當c>f時，表示算法陷入了局部最優(yōu). 為了跳出局部最優(yōu)，計算最后一次迭代得到的果蠅群體味道濃度的均方差，對小于均方差的果蠅進行混沌映射，再次進入尋優(yōu)的過程. 混沌映射產(chǎn)生的序列存在著不可預(yù)測、不可重復(fù)的特點，但是混沌映射有著遍歷性，能夠改善搜索算法的性能，提高搜索的精度.Logistic映射廣泛應(yīng)用于各種優(yōu)化算法中，但是該映射存在著遍歷性不足的缺點.Hariyanto 等人［24］提出了CAT 映射，該映射產(chǎn)生的序列是均勻分布的，相比于Logistic映射有著更好的遍歷性.

2種映射以相同的初值x0=0.3迭代4 000次，得到2種映射序列的遍歷分布圖和分布直方圖，如圖3 和圖4所示.Logistic映射集中分布在解空間的兩端，CAT映射在解空間內(nèi)的分布更均勻，Logistic 映射取值在［0，0.1］范圍內(nèi)的概率是0.204，取值在［0.9，1］范圍內(nèi)的概率是0.213 5，而CAT 映射取值范圍在［0，0.1］范圍內(nèi)的概率是0.093，取值在［0.9，1］范圍內(nèi)的概率是0.105 5.

圖3 2種映射序列的遍歷分布圖

圖4 2種映射序列的分布直方圖

傳統(tǒng)的CAT映射是二維的，不適用于有4個參數(shù)的Gabor 函數(shù)，高維廣義CAT 映射相比二維CAT 映射的Lyapunov 指數(shù)更大，遍歷性更好［25］，能進行多參數(shù)尋優(yōu). Gabor 函數(shù)的4 個參數(shù)中，對超聲回波波形影響較大的參數(shù)是尺度參數(shù)s和頻率參數(shù)v［26］，因此，只需對s和v進行映射即可，其動力學(xué)方程為

3.4 時間復(fù)雜度分析

在IFOA 算法中，已知果蠅種群大小為Sizepop，種群迭代次數(shù)為Maxgen，搜索空間的維度為m，設(shè)產(chǎn)生每一維果蠅所需的時間為t1，則果蠅位置初始化的時間復(fù)雜度為

設(shè)計算自適應(yīng)步長的時間為t2，計算果蠅味道濃度所需時間為f(m)，比較與更新最佳味道濃度的時間為t3，對果蠅進行高維廣義CAT 映射需要的時間為t4，更新果蠅起始位置需要的時間為t5，則IFOA-OMP 算法在迭代過程的時間復(fù)雜度為

于是IFOA算法的時間復(fù)雜度為

FOA 算法無需計算自適應(yīng)步長和映射，因此其時間復(fù)雜度為

由式（22）和式（23）可知，IFOA 算法和FOA 算法的時間復(fù)雜度相等，本文算法提高了FOA 算法的搜索精度，并沒有增加時間復(fù)雜度.

通過以上改進，得到改進的果蠅優(yōu)化算法流程圖如圖5所示.

圖5 改進的果蠅優(yōu)化算法流程圖

4 實驗驗證和結(jié)果分析

4.1 仿真超聲信號

為驗證本文方法對超聲回波信號的處理能力，首先對仿真的超聲信號進行處理. 超聲信號可用高斯調(diào)制的正弦波表示［27］，表達式為

其中，β是幅度系數(shù)，α是帶寬因子，fc是中心頻率，τ是回波到達時間，φ是相位.

材料的非均勻性會導(dǎo)致超聲波產(chǎn)生散射衰減，使接收信號的中心頻率低于探頭的中心頻率，因此本文模擬2 個不同頻率的超聲回波. 第1 個回波的fc1=5.8 MHz，β1=0.8，τ1=0.5 μs，第2 個回波的fc2=6 MHz，β2=1，τ2=1.2 μs，且α1=α2=50(MHz)2，φ1=φ2=8.89.仿真信號如圖6（a）所示，其頻譜如圖6（b）所示. 在超聲信號中加入-5 dB 的高斯白噪聲，時域波形和頻譜如圖6（c）和圖6（d）所示，可見，超聲信號完全淹沒在噪聲中.

圖6 仿真超聲信號時域波形和頻譜圖

采用IFOA-OMP 算法、果蠅優(yōu)化正交匹配追蹤（FOA-OMP）算法和OMP 算法對加噪仿真超聲信號進行處理，設(shè)置Maxgen 為120，Sizepop 為130，搜索步長初始值為20. 降噪結(jié)果如圖7 所示，經(jīng)過IFOA-OMP 算法處理后，重構(gòu)信號的波形沒有失真，但是2 個信號的頻譜幅值都有所下降，中心頻率偏移，這是由于原子與信號的殘差匹配存在誤差. 但相比于FOA-OMP 算法和OMP 算法，IFOA-OMP 算法降噪效果較好，F(xiàn)OA-OMP 算法和OMP算法僅能重構(gòu)信號的部分波形.

圖7 不同降噪方法結(jié)果對比

為了更有效地驗證本文方法降噪的效果，對不同信噪比的仿真超聲信號進行處理，重構(gòu)后回波信號的信噪比和均方根誤差如圖8 所示，可見IFOA-OMP 算法和FOA-OMP 算法在連續(xù)空間內(nèi)搜索最優(yōu)值，相較于OMP 算法，其降噪效果更好，證實了在連續(xù)空間搜索的精度是高于離散空間的，用參數(shù)尋優(yōu)的方法代替OMP算法中匹配追蹤的過程是可行的. 隨著噪聲的加大，3種方法的輸出信噪比均有所下降，均方根誤差逐漸增大，但IFOA-OMP算法在較低的信噪比情況下仍能輸出較大信噪比的信號，均方根誤差也更低.

圖8 不同信噪比超聲信號降噪效果對比

4.2 實際缺陷超聲信號

用SIUI CTS-4020數(shù)字超聲探傷儀采集鍛件試塊的超聲回波信號，試塊厚度40 mm，內(nèi)有直徑5 mm 的平底孔來模擬裂紋缺陷，平底孔距試塊表面30 mm. 探頭型號為5N10Z，中心頻率5 MHz，直徑10 mm，超聲波在試塊中的聲速為5 920 m/s，超聲儀采樣頻率為100 MHz.實驗信號如圖9（a）所示. 第一個回波是平底孔回波，回波幅值低；第二個回波是底面回波，底面反射面積大，回波幅值高. 經(jīng)本文算法處理后的結(jié)果如圖10 所示，從圖10（a）可見重構(gòu)信號波形更平滑. 對比圖9（b）和圖10（b），實驗信號中包含有雜亂的高頻噪聲，主要是儀器產(chǎn)生的電子噪聲，重構(gòu)信號中的高頻成份被抑制，低頻部分與原始信號基本吻合.

圖9 實驗超聲信號

圖10 IFOA-OMP處理結(jié)果

圖9（a）的實驗回波噪聲較小，為獲得低信噪比數(shù)據(jù)，減小超聲儀的發(fā)射功率，增大接收增益，采集到超聲回波信號如圖11 所示. IFOA-OMP 處理結(jié)果如圖12所示，可見該方法能有效去除信號中的噪聲成份，提取出超聲信號.

圖11 低功率超聲信號

圖12 IFOA-OMP處理低功率超聲信號結(jié)果

為模擬實際工業(yè)現(xiàn)場的強噪聲環(huán)境，對圖11（a）的超聲信號添加-5 dB 的高斯白噪聲，時域波形如圖13（a）所示，頻譜如圖13（b）所示，圖13（c）和圖13（d）分別是IFOA-OMP方法處理后的時域和頻譜圖.

圖13 加噪實驗超聲信號處理結(jié)果對比

不同算法去噪效果如表1 所示. OMP 算法處理分離回波信號的均方根誤差較大，輸出信噪比較低，分解精度不高. FOA-OMP 算法的搜索步長是隨機的，搜索精度不如IFOA-OMP 算法.IFOA-OMP 算法的均方根誤差最小，輸出信號的信噪比最大，去噪的效果最好，并且重構(gòu)的超聲信號中缺陷回波和底面回波與實驗超聲信號基本一致，信號頻譜特征得到了保留，而OMP算法和FOA-OMP算法處理后的信號已經(jīng)失真.

表1 不同算法去噪效果

5 結(jié)論

本文采用參數(shù)尋優(yōu)的思想，提出了基于改進果蠅優(yōu)化算法正交匹配追蹤（IFOA-OMP）的超聲信號降噪方法，在連續(xù)空間中尋找最佳原子，顯著提高了算法的降噪能力和計算效率.

IFOA-OMP 將正交匹配追蹤中的“貪婪”搜索轉(zhuǎn)換為Gabor函數(shù)的參數(shù)優(yōu)化問題，利用改進的果蠅優(yōu)化算法估計Gabor 函數(shù)的最優(yōu)值. 本文提出了自適應(yīng)步長的思想，前期采用較大步長快速收斂到最優(yōu)值附近，后期采用較小的步長以準確逼近最優(yōu)值，并可根據(jù)當前果蠅種群迭代次數(shù)調(diào)整步長的大小，保證算法在全局范圍內(nèi)的搜索精度. 同時為使算法能夠跳出局部最優(yōu)，本文引入高維廣義CAT 映射的方法改變Gabor 函數(shù)的尺度參數(shù)s和頻率參數(shù)v，以更新果蠅個體的位置，重新搜索到最優(yōu)值.

對仿真和實驗超聲信號的處理結(jié)果表明，IFOAOMP 能夠抑制強噪聲的干擾. 對添加-5 dB 的高斯白噪聲實驗信號的處理結(jié)果表明，采用IFOA-OMP重構(gòu)信號的信噪比為5.612 1 dB，均方根誤差為0.138 4，降噪效果明顯優(yōu)于FOA-OMP和OMP算法.