空域抽樣與相干因子融合的超聲陣列自適應(yīng)波束形成算法

宋壽鵬，李 棋

(江蘇大學(xué) 機(jī)械工程學(xué)院，江蘇 鎮(zhèn)江 212013)

自適應(yīng)波束形成技術(shù)是超聲成像的關(guān)鍵技術(shù)之一，是根據(jù)陣列中陣元接收的信號動態(tài)地計(jì)算權(quán)值并對信號進(jìn)行加權(quán)以控制波束進(jìn)行成像，直接決定著成像的質(zhì)量[1-3]。

自適應(yīng)算法中典型的最小方差（minimum variance，MV）算法可以在保持期望信號能量的同時(shí)，使波束形成器輸出的信號功率最小化[4-6]。為了進(jìn)一步提高M(jìn)V算法的成像性能，不少學(xué)者將波束形成算法與相干因子（coherence factor，CF）相融合，以達(dá)到提高成像分辨率的目的[7-9]。Wang等[10]提出將CF與最小方差無畸變（minimum variance distortionless response，MVDR）算法結(jié)合應(yīng)用于高幀率成像的算法，提高成像分辨率和對比度以降低旁瓣等級。在此基礎(chǔ)上，Asl等[11]指出該方法可以有效減小聚焦誤差，在聲速不均勻時(shí)表現(xiàn)出良好的魯棒性。Xu等[12]提出了MV算法與CF和相位相干因子（phase coherence factor，PCF）結(jié)合的方法，獲得了高對比度和高分辨率的醫(yī)學(xué)超聲圖像。雖然這些改進(jìn)和研究可以有效提升MV算法的成像分辨率，但是在計(jì)算復(fù)雜度上卻沒有改善。對于一個(gè)陣元數(shù)為N的線性陣列，自適應(yīng)MV算法的計(jì)算復(fù)雜度高達(dá)O(N3)，嚴(yán)重影響了超聲成像的時(shí)效性。

針對MV算法中由于協(xié)方差矩陣反演造成的計(jì)算復(fù)雜度高且計(jì)算效率低的問題，不少學(xué)者在這一領(lǐng)域開展了研究。Park等[13]提出一種結(jié)合FPGA的波束形成方法，通過減少計(jì)算量達(dá)到降低成像計(jì)算復(fù)雜度的目的，加快了成像速度。Chen等[14]提出一種基于坐標(biāo)旋轉(zhuǎn)數(shù)字計(jì)算機(jī)（coordinate rotation digital computer，CORDIC）處理器的低復(fù)雜度、高吞吐量的波束形成算法。這些研究雖然有效提高了超聲成像效率，但本質(zhì)上是將軟件算法與硬件結(jié)合，借助硬件提高成像速度。一些學(xué)者進(jìn)一步改進(jìn)MV算法中的協(xié)方差矩陣，通過降低矩陣反演的計(jì)算復(fù)雜度等級來提高成像速度。Kim等[15]利用主成分分析（principal component analysis，PCA）實(shí)現(xiàn)降維，通過預(yù)先計(jì)算的常規(guī)MV權(quán)值離線估計(jì)主成分并通過選定的主成分的線性組合逼近MV權(quán)值。Park等[16]提出運(yùn)用QR分解技術(shù)將空間協(xié)方差矩陣轉(zhuǎn)換為標(biāo)量矩陣，無需對矩陣反演即可得到自適應(yīng)加權(quán)值和波束形成器的輸出。這些方法雖然可以使協(xié)方差矩陣反演的計(jì)算復(fù)雜度降低，但是計(jì)算得到的加權(quán)值是近似值并非實(shí)際值。Nilsen等[17]研究了波束空間（beam space，BS）域中的MV波束形成方法，把數(shù)據(jù)從陣列空間轉(zhuǎn)換至波束空間，這一操作可將矩陣反演計(jì)算復(fù)雜度等級由3維降至2維，但存在方向矢量誤差時(shí)，該方法的性能可能會下降。Asl等[18]通過對采樣協(xié)方差矩陣作近似，將其為Toeplitz矩陣后再進(jìn)行矩陣反演。雖然上述方法中矩陣反演的計(jì)算復(fù)雜度等級降至2維，但由于參與運(yùn)算的數(shù)據(jù)量大，因此計(jì)算復(fù)雜度仍然很高。為此，一些學(xué)者開展了減少采樣陣元數(shù)目的研究以減少處理數(shù)據(jù)量。

Sakhaei[19]描述了一種抽樣MV波束形成器算法，首先結(jié)合接收波束的波束模式的分析對全部數(shù)據(jù)抽樣，然后利用全部數(shù)據(jù)計(jì)算出加權(quán)系數(shù)并對抽樣數(shù)據(jù)進(jìn)行加權(quán)，最后在醫(yī)學(xué)成像上驗(yàn)證了其可行性。在此基礎(chǔ)上，Shamsian等[20]提出級聯(lián)結(jié)構(gòu)的快速波束形成方法，將低通濾波器與最小方差波束形成器級聯(lián)。其中：第1階段，通過低通濾波器去除回波信號中的離軸噪聲，再對數(shù)據(jù)進(jìn)行抽樣；第2階段，將抽樣數(shù)據(jù)作為MV算法的輸入以抑制軸上干擾，并提出抽取子波束的MV算法（decimated sub-beam MV，DSMV）。該方法實(shí)現(xiàn)了數(shù)據(jù)量和計(jì)算復(fù)雜度上的降低，但對陣元數(shù)據(jù)進(jìn)行抽樣這一過程不可避免會丟失部分信息，進(jìn)而影響成像質(zhì)量。

為了在保證成像分辨率的同時(shí)，降低算法計(jì)算復(fù)雜度進(jìn)而提高計(jì)算效率，本文提出了一種空域抽樣與相干因子融合的超聲陣列自適應(yīng)波束形成算法（decimated minimum variance combined with coherence factor，DCFMV）。該算法先利用數(shù)據(jù)空域抽樣以減少數(shù)據(jù)量，然后基于MV原則計(jì)算自適應(yīng)加權(quán)值，并融合相干因子對期望信號進(jìn)行增強(qiáng)。為了進(jìn)一步簡化運(yùn)算，在對數(shù)據(jù)處理時(shí)，將數(shù)據(jù)協(xié)方差矩陣構(gòu)造為Toeplitz矩陣的結(jié)構(gòu)。

1 空域抽樣與相干因子融合的自適應(yīng)波束形成算法

1.1 陣列回波信號的空域抽樣

圖1 陣列波束對比Fig.1 Contrast of the array beam

對完備陣列數(shù)據(jù)進(jìn)行空域抽樣，那么其接收波束模式B(u)可以表示為：

若抽取因子D取值較小，表明陣列數(shù)據(jù)信息丟失少，成像效果好，但數(shù)據(jù)量減少不明顯，成像時(shí)效性也得不到明顯提升；若抽取因子D取值較大，則陣列數(shù)據(jù)丟失的信息較多，雖然成像時(shí)效性可以明顯提高，但犧牲了成像質(zhì)量。因此，為解決成像質(zhì)量與時(shí)效性的矛盾，在空域數(shù)據(jù)抽取中選擇合適的抽取因子D至關(guān)重要。

如果抽取的非完備數(shù)據(jù)通過柵瓣泄露的功率與完備數(shù)據(jù)的總功率相比可以忽略不計(jì)，那么，最大抽取因子Dmax近似滿足[20]：

將u1代入式（2）可以得到：

考慮到抽取因子D和最大抽取因子Dmax均為正數(shù)，則存在：

確定抽取因子D之后，對接收陣元按抽取因子進(jìn)行空域等間隔抽樣。以N=16、D=2為例，空域抽樣數(shù)據(jù)原理示意圖如圖2所示。由圖2可知：發(fā)射陣元采用全矩陣發(fā)射模式；各陣元發(fā)射時(shí)，接收陣元按D=2等間隔抽樣，得到空域抽樣數(shù)據(jù)。

圖2 空域抽樣數(shù)據(jù)原理示意圖Fig.2 Schematic diagram of spatial sampling data

1.2 最小方差波束形成算法

得到抽樣數(shù)據(jù)后，采用最小方差波束形成方法使陣列通道所接收到的噪聲信號能量和非期望方向上的干擾信號能量達(dá)到最小，以此獲得最優(yōu)權(quán)值。將空域抽樣后的信號作為波束形成器的輸入，那么，波束形成器的輸出y(k)可以表示為：

式中：k為不同采樣點(diǎn)的序列號；w為加權(quán)列向量；[·]H為共軛轉(zhuǎn)置；x(k)=[x1(k-Δ1),x2(k-Δ2),···,xi(k-Δi),···,xM(k-ΔM)]T， 為延時(shí)后的空域抽樣信號，其中，xi(k)為第i個(gè)接收陣元在k時(shí)刻的接收回波， Δi為第i個(gè)接收陣元的延時(shí)時(shí)間， [·]T為轉(zhuǎn)置，下標(biāo)M為空域抽樣后的接收陣元數(shù)。

在約束條件為wHa=1的情況下，通過最小化輸出功率即 min(wHRw)求解最優(yōu)權(quán)向量。利用拉格朗日乘子法求解最優(yōu)權(quán)向量wopt，得到：

在實(shí)際檢測中，由于陣列中各陣元接收到的回波信號來自同一區(qū)域內(nèi)的檢測對象的反射或散射，因此，信號具有高相關(guān)性使得計(jì)算所得協(xié)方差矩陣R奇異，從而導(dǎo)致加權(quán)值w難以求得。為了解決加權(quán)值w求解這一難題，需要對回波信號解相關(guān)。解相關(guān)就是先將協(xié)方差矩陣進(jìn)行劃分，分為P=M-L+1個(gè)重疊子陣，子陣長度L＜M/2；然后，計(jì)算求得每個(gè)子陣的協(xié)方差矩陣；最后將每個(gè)子陣的協(xié)方差矩陣相加取平均值。最終得到平滑處理后的協(xié)方差矩陣：

式中，xp(k) 為第p個(gè)子陣的接收數(shù)據(jù)。

為增強(qiáng)MV算法對建模誤差的魯棒性，對協(xié)方差矩陣R進(jìn)行修正，在協(xié)方差矩陣中引入一個(gè)加載因子ε，也就是用R+εI代替R。 其中：I為單位陣列，從而起到壓縮干擾信號的作用；ε的取值遵循[21]：

式中，t r{·}為 矩陣的跡，γ為一恒定常數(shù)。

回波信號經(jīng)解相關(guān)和協(xié)方差矩陣修正后，波束形成器的輸出y(k)可改寫為：

式中，第p個(gè) 子陣在k點(diǎn)的加權(quán)值wp(k)為：

1.3 相干因子計(jì)算

將空域抽樣后的數(shù)據(jù)作為最小方差波束形成算法的輸入，可以有效減少算法中的數(shù)據(jù)處理量，但是，相比于完備數(shù)據(jù)，此時(shí)的輸入信號中不可避免地丟失了部分有用信息。

為增強(qiáng)輸入信號中的有效信息，引入相干因子。相干因子是通過計(jì)算目標(biāo)點(diǎn)處相干能量與總能量之比，給相干能量較大的目標(biāo)點(diǎn)賦予更大的權(quán)值Cf(k)，從而提高超聲回波圖像中相干能量，通過提高有效信息的區(qū)分度改善超聲回波圖像質(zhì)量。相干因子的計(jì)算公式為[10]：

式中，Sxg為空域抽樣信號中的相干能量，Ssum為空域抽樣信號的總能量。Cf(k)取值在0～1之間，取值越大，表示回波信號中相干能量占比越高。

對于采樣點(diǎn)k，空域抽樣與相干因子融合的波束形成算法的輸出修正為：

2 協(xié)方差矩陣反演的計(jì)算過程與復(fù)雜度分析

理想情況下協(xié)方差矩陣R具有Toeplitz矩陣的結(jié)構(gòu)，但實(shí)際中通常會出現(xiàn)相干源和低信噪比的情況，因此不滿足該結(jié)構(gòu)。在對R反演時(shí)，為了降低計(jì)算復(fù)雜度等級將其構(gòu)造為Toeplitz矩陣，該方法在有效降低求取加權(quán)值的計(jì)算復(fù)雜度等級的同時(shí)，也可以使R接近真實(shí)的數(shù)據(jù)協(xié)方差矩陣。

協(xié)方差矩陣R為M×M維矩陣，將R下三角和上三角部分各條對角線上元素取平均得到Toeplitz矩陣RT。RT每一條對角線上的任意兩個(gè)值都是相等的，任意元素的計(jì)算公式如下：

式中，ri,(i+q)為矩陣R中 第i行 (i+q) 列的元素，r(i+q),i為矩陣R中 第 (i+q) 行i列的元素，r-q為矩陣RT上三角部分第q條次對角線上的任一元素的值，rq為 矩陣RT下三角部分第q條次對角線上的任一元素的值。

遞推求出RT的 逆矩陣R-T1的各列，具體流程如下：

1）根據(jù)Zohar算法[22-23]，求解線性方程組：

式中，eM=(0,0,···,1)T，f=(r1-r1-M,r2-r2-M,···,rM-1-r-1,0)T。求得方程組（15）的唯一解為x與y，這一步驟共需 4M2+4 次乘除運(yùn)算、4M2-5M+8次加減運(yùn)算。

2）根據(jù)Toeplitz矩陣定理，R-T1的列向量zj(j=1,2,···,M)滿足如下關(guān)系[24]：

遞推求出矩陣的逆的各列元素，需要M2-M次乘除運(yùn)算、M2-M次加減運(yùn)算。

綜合步驟1）和2），一共需要 5M2-M+4次乘除運(yùn)算、5M2-M+4 次加減運(yùn)算。因此，對于M×M維矩陣R，經(jīng)過Toeplitz結(jié)構(gòu)化后反演的計(jì)算復(fù)雜度為O(M2)，即O(N2/D2)。

將本文提出的空域抽樣與相干因子融合的最小方差波束形成方法（decimated minimum variance combined with coherence factor，DCFMV）與傳統(tǒng)MV算法及融合相干因子的最小方差波束形成算法[11]（minimum variance combined with coherence factor，CFMV）計(jì)算協(xié)方差矩陣上的計(jì)算復(fù)雜度進(jìn)行對比，結(jié)果如表1所示。

表1 協(xié)方差矩陣反演的計(jì)算復(fù)雜度對比Tab. 1 Comparison of computational complexity of covariance matrix inversion

由于抽取因子D≥1，因此從表1中可直觀看出DCFMV方法比MV和CFMV算法計(jì)算復(fù)雜度明顯降低。

3 仿真與討論

為了研究DCFMV算法的性能，利用Field Ⅱ仿真軟件對本文提出的空域抽樣與相干因子融合的最小方差波束形成方法（DCFMV）、傳統(tǒng)MV算法、融合相干因子的最小方差波束形成算法（CFMV）[11]這3種方法進(jìn)行成像對比。仿真試驗(yàn)中選用裂紋和橫通孔兩種缺陷，結(jié)構(gòu)和分布示意圖分別如圖3（a）、（b）所示。試塊上設(shè)置裂紋尺寸為 0.2×10×20mm3，橫通孔為4個(gè)直徑為2 mm的通孔缺陷，缺陷之間間隔為20 mm，后續(xù)分析中從上至下依次命名為孔1、孔2、孔3和孔4。

圖3 缺陷分布示意圖Fig.3 Schematic diagram of the defect distribution

試驗(yàn)中，超聲換能器為線性陣列，陣元數(shù)N=64，各陣元中心距為0.6 mm，陣元寬度為0.5 mm，中心頻率為5 MHz，采樣頻率為100 MHz，換能器采用垂直入射方式，超聲波在被測試件中的聲速為5 900 m/s。子陣列長度L=N/2 ，對角加載系數(shù) γ=1/10L，抽取因子D=4，采用動態(tài)孔徑的成像方法，成像深度與孔徑的比值為定值2。設(shè)置超聲圖像動態(tài)顯示范圍為60 dB。使用的計(jì)算機(jī)配置為Intel（R）Core（TM）i5-8250U，運(yùn)行內(nèi)存8 GB。

3.1 數(shù)據(jù)信息不對等時(shí)成像性能分析

在DCFMV、MV、CFMV3種算法下，裂紋和橫通孔缺陷成像結(jié)果分別如圖4和5所示，其中，MV和CFMV這兩種算法采用完備數(shù)據(jù)，DCFMV算法采用空域抽樣后的非完備數(shù)據(jù)。

圖4 數(shù)據(jù)不對等時(shí)裂紋缺陷成像效果對比Fig.4 Comparison of imaging effects of crack defects under unequal data

由圖4可知，CFMV算法成像時(shí)裂紋與實(shí)際缺陷尺寸誤差較大，成像結(jié)果的量化精度沒有MV和DCFMV算法準(zhǔn)確。

用陣列成像質(zhì)量性能指標(biāo)（array performance indicator，API）表征成像分辨率。API是點(diǎn)擴(kuò)散函數(shù)空間大小的無量綱度量，定義為點(diǎn)擴(kuò)散函數(shù)從其最大值向下大于-6 dB的區(qū)域A-6dB與波長的平方之比，計(jì)算式[25]如下：

由圖5可知，3種方法中，CFMV的成像效果最好，MV的成像效果最差。通過式（17）計(jì)算圖5中橫通孔缺陷圖像的API值，發(fā)現(xiàn)CFMV算法API值最小，DCFMV算法次之，MV算法效果最差。對比圖5中3種成像算法的主瓣寬度，發(fā)現(xiàn)CFMV算法主瓣最窄，DCFMV算法次之，MV算法主瓣最寬，但在數(shù)值上相差不明顯。

圖5 數(shù)據(jù)不對等時(shí)橫通孔缺陷成像效果對比Fig.5 Comparison of imaging effects of cross-drilled hole defects under unequal data

為了對比DCFMV、MV、CFMV這3種算法的計(jì)算復(fù)雜度，以成像耗時(shí)作為衡量指標(biāo)，分別對3種成像算法的成像時(shí)間進(jìn)行對比，15個(gè)樣本取平均，結(jié)果如表2所示。

表2 數(shù)據(jù)不對等時(shí)算法成像時(shí)間對比Tab. 2 Comparison of imaging time of algorithms under unequal data

由表2可知：通過參考15個(gè)樣本取平均，DCFMV算法成像時(shí)間大幅降低；對于裂紋缺陷，DCFMV算法的成像時(shí)間相較于MV算法和CFMV算法分別減小了85.12%和85.24%；對于橫通孔缺陷，DCFMV算法的成像時(shí)間相比于MV算法和CFMV算法分別縮短了86.11%和86.47%。

3.2 使用空域抽樣數(shù)據(jù)成像性能分析

利用空域抽樣數(shù)據(jù)分別對DCFMV、MV、CFMV這3種成像算法進(jìn)行了對比，其中抽取因子D=4，成像結(jié)果分別如圖6和7所示。

由圖6、7可知，在3種算法下，MV算法成像效果最差。為了更好地比較3種成像算法的成像質(zhì)量，以API作為成像質(zhì)量衡量指標(biāo)，在對裂紋缺陷成像時(shí)，MV算法、CFMV算法及DCFMV算法計(jì)算得到的API值分別為3.53、1.20、1.11，DCFMV算法的API值分別比MV算法和CFMV算法降低了68.77%和8.34%，表明在空域抽樣數(shù)據(jù)下DCFMV算法的成像質(zhì)量好于MV算法和CFMV算法；在對通孔缺陷成像時(shí)，DCFMV算法的API值都大幅度降低，4個(gè)孔類缺陷從上至下的API值與MV算法相比分別降低了70.91%、61.86%、66.62%、66.88%，與CFMV算法相比分別降低了17.77%、7.74%、30.07%、11.99%。綜上所述，DCFMV算法下的圖像質(zhì)量優(yōu)于其他兩種算法。

圖6 空域抽樣數(shù)據(jù)下裂紋缺陷成像效果對比Fig.6 Comparison of imaging effects of crack defects under spatial sampling data

圖7 空域抽樣數(shù)據(jù)下橫通孔缺陷成像效果對比Fig.7 Comparison of imaging effects of cross-drilled hole defects under spatial sampling data

以通孔類缺陷為例，分別計(jì)算采用DCFMV、MV、CFMV這3種算法的4個(gè)不同深度孔類缺陷主、旁瓣的分布，結(jié)果如圖8所示。

從圖8中可以看出，DCFMV算法的主瓣寬度最小。

圖8 4個(gè)通孔缺陷主、旁瓣分布Fig.8 Main and side lobes of four cross-drilled hole defects

為了比較DCFMV、MV、CFMV這3種算法在空域抽樣數(shù)據(jù)下的成像時(shí)間，分別計(jì)算了3種成像算法的運(yùn)行時(shí)間，15個(gè)樣本取平均，結(jié)果如表3所示。

表3 空域抽樣數(shù)據(jù)下算法成像時(shí)間對比Tab. 3 Comparison of imaging time of algorithms under spatial sampling data

由表3可以看出：CFMV算法耗時(shí)最長，MV算法次之，DCFMV算法耗時(shí)最少；對于裂紋缺陷，DCFMV算法相較于MV算法和CFMV算法成像時(shí)間分別減小了65.90%和66.10%；對于橫通孔缺陷，DCFMV算法相較于MV和CFMV算法成像時(shí)間分別縮短了67.49%和68.61%。

4 結(jié) 論

為了降低波束形成算法缺陷成像時(shí)的計(jì)算復(fù)雜度，提高成像運(yùn)行時(shí)效性，本文提出了一種超聲陣列波束形成方法。該方法根據(jù)陣列陣元數(shù)確定最大抽取因子，對陣元數(shù)據(jù)空域抽樣后，數(shù)據(jù)量得到大幅減少。將協(xié)方差矩陣R構(gòu)造為Toeplitz矩陣對其進(jìn)行修正，算法計(jì)算復(fù)雜度降低；融合相干因子提高抽樣數(shù)據(jù)中有效信息占比，進(jìn)而提高成像質(zhì)量。DCFMV算法的計(jì)算復(fù)雜度由O(N3) 降為O(N2/D2)。通過在試塊中裂紋和通孔缺陷的仿真成像，DCFMV算法的成像質(zhì)量較好，與MV算法和CFMV算法相比，DCFMV算法時(shí)間效率提高分別提高85%和86%以上。數(shù)據(jù)空域抽樣后，DCFMV算法下的成像API最小，成像質(zhì)量最好。該算法減小了成像算法對硬件性能的過度依賴，可為超聲波束形成及缺陷成像提供理論參考。下一步將結(jié)合等間隔抽樣的信息損失評估，相干因子融合的效果評估，在矛盾制約條件下挖掘出最佳抽樣相干融合策略。