基于高斯回歸預(yù)測的超聲成像高分辨率重建*

劉皓，譚超，董峰

(天津大學(xué)電氣自動化與信息工程學(xué)院 天津市過程檢測與控制重點實驗室 天津 300072)

多相流廣泛存在于自然界和人類生產(chǎn)過程，如能源、動力、冶金、化工、宇航、醫(yī)藥等現(xiàn)代工程領(lǐng)域[1]。其中，油水兩相流是在石油、化工等生產(chǎn)過程中廣泛存在的一種多相流形式，對其流動狀態(tài)的掌握和流動參數(shù)的獲取，是保障生產(chǎn)過程安全穩(wěn)定運行的前提[2-3]。過程層析成像技術(shù)作為一種可視化檢測技術(shù)，可實現(xiàn)多相流流動狀態(tài)和流動參數(shù)的在線獲取。根據(jù)激勵源能量的種類不同，過程層析成像技術(shù)主要有光電法[4]、超聲法[5]及射線法[6]等。與光電法和射線法相比，超聲層析成像技術(shù)由于具有非侵入、無輻射、安裝方便等優(yōu)點，在油水兩相流檢測中有很好的應(yīng)用前景。

超聲層析成像技術(shù)的圖像重建是利用超聲波在被測場域內(nèi)基本上是沿直線傳播并線性衰減的特性，在接收信號和被測物體間形成基于路徑的線性關(guān)系。通過比較激勵信號和接收信號的幅值差異，獲得收發(fā)探頭間投影路徑上聲衰減的平均估計，并采用圖像重建算法，通過對場域內(nèi)各位置衰減系數(shù)的估計，實現(xiàn)圖像重建。

傳統(tǒng)的圖像重建方法包括線性反投影(linear back projection，LBP)[7]、濾波反投影(filter back projection，F(xiàn)BP)[8]等算法，這一類算法成像實時性高，但重建精度較差，無法實現(xiàn)對油水兩相流等弱聲阻抗比介質(zhì)的有效重建。另一類圖像重建方法為基于迭代和約束項的算法，如總變差正則化(total variation，TV)[9]、代數(shù)重建(algebraic reconstruction technique，ART)[10]、同步迭代重建(simultaneous iteration reconstruction technique，SIRT)[11]等算法，該類算法將被測物場的均勻分布狀態(tài)作為成像的初始估計值，在每次迭代中加入約束項，實現(xiàn)場域內(nèi)介質(zhì)分布的最優(yōu)搜索，該類算法對噪聲有一定的魯棒性，并能穩(wěn)定收斂至最優(yōu)解。

對于廣泛存在的弱聲阻抗比兩相介質(zhì)，由于其界面透射因數(shù)遠(yuǎn)大于反射因數(shù)，通常采用透射模態(tài)的成像方法。在透射模態(tài)基礎(chǔ)上，對低聲阻抗比內(nèi)含物進(jìn)行較高精度超聲重建主要依賴兩方面因素：首先是獲得包含準(zhǔn)確衰減信息的高信噪比邊界測試數(shù)據(jù)，其次是設(shè)計高精度的圖像重建算法，在滿足實時性要求下達(dá)到盡可能高的精度和分辨率。事實上，超聲成像圖像重建的質(zhì)量一直嚴(yán)重依賴于探頭安裝數(shù)量，過多的探頭數(shù)量對超聲信號收發(fā)裝置的設(shè)計提出了較高的要求，不適用于工業(yè)多相流的測量；另一方面，換能器數(shù)量的增加大幅提高了圖像重建算法中系數(shù)矩陣的維度，對系數(shù)矩陣的求逆過程占用大量的計算資源和過長的計算時間，使其無法動態(tài)反映管道內(nèi)油水兩相流的實時變化情況。

針對采用超聲層析成像實現(xiàn)低對比度聲學(xué)特性的油、水兩相介質(zhì)分布重建問題，提出一種基于高斯回歸預(yù)測的高分辨率相介質(zhì)分布圖像重建算法。在獲取采用超聲連續(xù)波激勵方式的邊界測量數(shù)據(jù)的基礎(chǔ)上，提取邊界測量聲壓幅值中包含的超聲衰減信息，通過采用能夠同時兼顧重建速度和重建精度的基于截斷奇異值分解的同步代數(shù)迭代重建算法，并構(gòu)建高斯回歸模型，將低分辨率重建結(jié)果轉(zhuǎn)化為高分辨率圖像。仿真模擬測試結(jié)果表明，所提出的算法計算速度快、重建誤差小、圖像分辨率高，為采用超聲層析成像技術(shù)實現(xiàn)油水兩相流動態(tài)可視化檢測提供了有效的解決方案。

1 超聲層析成像原理與數(shù)學(xué)模型

1.1 超聲層析成像原理

超聲層析成像的基本原理是對收發(fā)信號間幅值、相位、時間的差異進(jìn)行分析，確定內(nèi)含物的密度、聲速、聲衰減系數(shù)。超聲層析成像系統(tǒng)一般包含：超聲換能器陣列、超聲信號激勵采集裝置、主控計算機及圖像重建算法等部分，組成結(jié)構(gòu)如圖 1 所示。

圖1 超聲層析成像基本原理Fig.1 Basic schematic of ultrasound tomography

圖1中，超聲換能器陣列由多個超聲換能器組成，等間距排布在被測場域的外圍并由統(tǒng)一的多路復(fù)用器控制激勵和采集時序。當(dāng)選擇某一特定探頭進(jìn)行激勵時，采用一定幅值的正弦電信號加載在壓電晶片兩端，產(chǎn)生柱面波形式的超聲波振動，接收探頭記錄相應(yīng)激勵下的聲壓幅值響應(yīng)。按照“一發(fā)全收”的方式進(jìn)行循環(huán)激勵與采集，對每個超聲探頭分別進(jìn)行激勵以獲取不同投影角度的超聲幅值響應(yīng)信息，總計獲得N×(N-1)組時變超聲信號(N為探頭數(shù)量)，通過采集裝置將解調(diào)數(shù)據(jù)發(fā)送給主控計算機，重建被測場域衰減系數(shù)分布。

1.2 超聲層析成像數(shù)學(xué)模型

超聲在場域內(nèi)含物中的衰減過程較為復(fù)雜，耦合了散射衰減、擴(kuò)散衰減和吸收衰減等多種機制。作為各種衰減機制的綜合表征，衰減系數(shù)是指單位距離內(nèi)收發(fā)聲壓信號的幅值比例差異。一般意義上，對收發(fā)探頭間特定路徑的平均聲衰減的估計可以表示為

(1)

式中：As是單介質(zhì)場(全水)時的接收信號，Ar是內(nèi)含物(油泡)存在時的接收信號，α0表示單介質(zhì)場的聲衰減系數(shù)。

超聲層析成像的重建過程為：被測場域共包含n個像素，測試數(shù)據(jù)共構(gòu)成m條路徑，隨著第j個像素在第i條路徑上位置的確定，像素上的聲衰減系數(shù)用式(1)進(jìn)行求解。將其簡化為矩陣表示的形式：

R·a=τ，

(2)

R·a=τ+τnoise，

(3)

式中：τnoise表示加性噪聲。超聲成像圖像重建算法需根據(jù)預(yù)先計算的系數(shù)矩陣R和從傳感器陣列中得到的邊界測量數(shù)據(jù)τ計算衰減系數(shù)分布a。在式(2)和式(3)中，系數(shù)矩陣R的計算過程一般稱為“正問題求解”過程。本文采用面積占比法[5]對系數(shù)矩陣進(jìn)行求解，系數(shù)矩陣中的元素rij為第i條投影路徑和第j個像素之間的重疊面積占該像素面積的比例。

2 圖像重建算法

2.1 同步代數(shù)迭代重建方法

在使用傳統(tǒng)的反投影方法對式(2)進(jìn)行逆問題求解時，往往將矩陣R-1等效成RT進(jìn)行重建計算，這一方法嚴(yán)重影響重建的精度。1984年，Anderson和Kak[12]提出同步代數(shù)重建(simultaneous algebraic reconstruction technique,SART)的方法，有效融合了ART和SIRT的優(yōu)點，在保持算法魯棒性的基礎(chǔ)上可大幅提升成像質(zhì)量和算法收斂速度，其迭代過程表示為

(4)

式中：rij代表系數(shù)矩陣R中的元素；aj表示第j個像素的聲衰減系數(shù)；τi表示第i條路徑的信號衰減；λ表示迭代過程中的松弛因子；P為高通濾波器模板，表示為

(5)

式中：δ(x-xi,y-yi)代表二維狄拉克函數(shù)，w為頻率系數(shù)，(x,y)表示濾波器模板中像素坐標(biāo)，(xi,yi)表示濾波器模板目標(biāo)像素的坐標(biāo)。將SART的迭代過程轉(zhuǎn)化成矩陣形式，表示為

a(k+1)=a(k)+λSp(SrR)T(τ-Ra(k))，

(6)

式中：

(7)

2.2 基于截斷奇異值分解的SART算法

SART算法在迭代過程中，計算速度受每步運算中的矩陣維度影響較大。為提升算法速度以滿足成像實時性要求，需要對SART算法中系數(shù)矩陣的維度進(jìn)行壓縮。根據(jù)系數(shù)矩陣的構(gòu)建方式及其稀疏性特點，選取截斷奇異值分解(truncated singular value decomposition，TSVD)方法對系數(shù)矩陣R進(jìn)行降維處理，進(jìn)而提升重建算法速度以滿足實時性要求。

TSVD作為一種常見的特征提取方法，廣泛用于主成分分析及大型復(fù)雜系統(tǒng)識別，其數(shù)學(xué)基礎(chǔ)為奇異值分解(singular value decomposition,SVD)。對于給定的系數(shù)矩陣R，其奇異值分解形式表示為

R=UΣVT，

(8)

Vq=span{v1,v2,…,vq}.

(9)

式中：q表示截斷系數(shù)，{v1,v2,…,vq}為矩陣R最大的q個奇異值對應(yīng)的矩陣V中的列。令Vq=V,得到TSVD的基本形式為

R=UΣVT≈UqΣqVqT.

(10)

在對系數(shù)矩陣R進(jìn)行截斷的同時，需要給出SART每步迭代計算的TSVD形式。首先在子空間中，對a進(jìn)行降維處理：

(11)

將式(10)、式(11)代入式(4)并將V替換成Vq，可以得到

(12)

(13)

式中：Gq是正交變換后的系數(shù)矩陣，Uq和Σq分別是變換后的矩陣U和矩陣S。

對于涉及濾波的SART算法而言，濾波器模板P同樣也需要進(jìn)行降維處理，對應(yīng)低維的迭代算法表示為

(14)

2.3 高斯回歸預(yù)測圖像重建方法

為提升圖像重建精度及分辨率，在采用基于TSVD的SART算法實現(xiàn)圖像初步重建的基礎(chǔ)上，進(jìn)一步采用高斯回歸模型(Gaussian process regression，GPR)對低分辨率重建結(jié)果進(jìn)行處理，實現(xiàn)高分辨率重建圖像。

GPR是一種實現(xiàn)各種統(tǒng)計學(xué)習(xí)問題的常用工具，它通過對訓(xùn)練數(shù)據(jù)樣本，獲得輸入輸出間的聯(lián)合概率分布和預(yù)測數(shù)據(jù)的先驗概率分布，并計算預(yù)測數(shù)據(jù)輸出的后驗概率分布。與卷積神經(jīng)網(wǎng)絡(luò)(convolutional neural network，CNN)相比，GPR的訓(xùn)練集與測試集的輸入數(shù)據(jù)不必保證數(shù)據(jù)維度一致。與線性回歸為代表的擬合方法相比，GPR可以有效減少局部測試集輸入數(shù)據(jù)誤差所引入的波動，對訓(xùn)練樣本中的噪聲和擾動有較好的魯棒性[13]。因此，本文中采用GPR將迭代求解的低分辨率重建結(jié)果映射到高分辨率重建圖像，在保持成像主體對象高分辨率重建的同時，將重建結(jié)果中的偽影和圖像噪聲在映射過程中有效去除，提升圖像純凈度。

給定訓(xùn)練集D{(xi,yi),i=1,…n}={X,y}，其中X表示低分辨率重建結(jié)果各個節(jié)點坐標(biāo)，y表示各個節(jié)點重建的衰減系數(shù)值。假定xi至yi的映射f(xi)符合正態(tài)分布，則x至y的映射可以表示成為高斯過程(Gaussian process)，其概率分布表示為

(15)

式中：m(x)表示均值，k(x,x′)表示協(xié)方差矩陣。設(shè)數(shù)據(jù)集中要預(yù)測的數(shù)據(jù)輸入為X*，輸出為f*，則根據(jù)貝葉斯公式

(16)

GPR的核心是根據(jù)訓(xùn)練數(shù)據(jù)的概率分布和訓(xùn)練及預(yù)測數(shù)據(jù)的聯(lián)合概率分布，計算預(yù)測數(shù)據(jù)輸出的后驗概率分布。訓(xùn)練數(shù)據(jù)的概率分布p(f)由式(15)給出，訓(xùn)練數(shù)據(jù)與預(yù)測數(shù)據(jù)的聯(lián)合概率分布為

(17)

式中：K(X*,X)為半正定的協(xié)方差矩陣。根據(jù)式(15)、式(17)，預(yù)測數(shù)據(jù)輸出的聯(lián)合后驗概率分布為

(18)

對于半正定協(xié)方差矩陣的選取，考慮到訓(xùn)練數(shù)據(jù)于測試數(shù)據(jù)的輸入均為坐標(biāo)值，具有數(shù)值上的周期性和范數(shù)上的單調(diào)性。因此選用平方指數(shù)協(xié)方差函數(shù)作為高斯回歸模型的核函數(shù)，其計算方式為

(19)

高斯回歸預(yù)測高分辨率成像結(jié)果的具體實現(xiàn)步驟為：

1)使用前述降維SART方法獲得低分辨率成像結(jié)果(X,y)，其中X表示低分辨率成像結(jié)果各像素坐標(biāo)，y表示低分辨率成像結(jié)果像素值。

2)基于選定的核函數(shù)，根據(jù)式(15)計算高斯回歸模型參數(shù)。

3)基于高分辨率圖像各像素點坐標(biāo)X*，根據(jù)式(18)計算高分辨率像素的概率密度分布f*，并求得各像素取值y*。

采用高斯回歸預(yù)測的超聲層析成像圖像重建算法(SART-GPR)流程圖如圖 2所示。算法中，首先對邊界測量數(shù)據(jù)進(jìn)行處理并根據(jù)面積占比法進(jìn)行系數(shù)矩陣R的計算；然后使用TSVD降維并基于SART算法進(jìn)行迭代計算；當(dāng)殘差小于設(shè)定閾值時，通過GPR進(jìn)行預(yù)測并給出最終的高分辨率成像結(jié)果。

3 仿真結(jié)果和討論

3.1 仿真設(shè)置

為驗證算法在成像精度、計算速度方面的優(yōu)勢，采用仿真實驗的方法驗證對不同情況下的內(nèi)含物結(jié)構(gòu)進(jìn)行成像測試。仿真實驗的數(shù)值計算通過COMSOL多物理場仿真軟件實現(xiàn)，以使其具有強大的偏微分波動方程求解能力和最優(yōu)解迭代快速收斂的優(yōu)點。

圖2 高斯回歸預(yù)測重建算法流程Fig.2 Flow chart of the proposed reconstruction method

邊界測試數(shù)據(jù)求解過程中，場域邊界的聲壓分布由狄利克雷方程控制，物理場的控制方程為基于動量守恒和質(zhì)量守恒的等熵二維波動方程為

(20)

式中：ρ0為背景介質(zhì)的密度，p為質(zhì)點的聲壓分布，c代表質(zhì)點處的聲速。

使用有限元法對超聲傳播過程進(jìn)行數(shù)值計算的過程中，場域內(nèi)網(wǎng)格剖分的尺寸為波長的1/6，保證在計算精度的同時盡可能減小計算資源的占用。仿真計算的相關(guān)參數(shù)設(shè)置詳見表 1。

表1 仿真參數(shù)設(shè)置 Table 1 Simulation parameter setting

在采用仿真測試數(shù)據(jù)的重建中，采用在超聲層析成像中廣泛應(yīng)用的圖像重建算法與本文所提出的SART-GPR方法進(jìn)行比較(其他算法使用線性插值將分辨率提升到與SART-GPR相同的水平)。為定量比較分析重建效果，定義相對誤差和相關(guān)系數(shù)兩個重建指標(biāo)：

(21)

(22)

3.2 成像結(jié)果

圖 3 給出TV、FBP、未經(jīng)優(yōu)化的SART與本文提出的SART-GPR算法的單內(nèi)含物和雙內(nèi)含物分布成像結(jié)果，所有成像結(jié)果經(jīng)歸一化處理以便比較。從圖中可以看出，傳統(tǒng)的TV、FBP及未經(jīng)優(yōu)化的SART算法成像結(jié)果偽影均較嚴(yán)重，圖像純凈度低；SART-GPR算法的成像結(jié)果幾乎無偽影，且圖像純凈度相比其他算法有顯著的提升。

圖3 不同算法重建結(jié)果Fig.3 Reconstruction results using different methods

進(jìn)一步根據(jù)式(21)、式(22)計算成像結(jié)果的相關(guān)系數(shù)和相對誤差，結(jié)果如圖 4 所示。圖 4 中，SART-GPR算法的平均相對誤差減小至其他3種算法的1/2以下，達(dá)到25.3%；平均相關(guān)系數(shù)與其他3種算法相比也有所提升，達(dá)到88.6%。

圖4 圖像重建定量指標(biāo)對比Fig.4 Quantitative comparison of image reconstruction

SART-GPR算法的成像結(jié)果受投影數(shù)量(即探頭數(shù)目)影響，隨著探頭數(shù)量的增加，成像結(jié)果的邊緣梯度(保邊性)及內(nèi)含物尺寸的重建精度均有明顯提高。針對圖 3 中的模型1至模型4，SART-GPR算法在采用16、32、48和64探頭時，成像的平均相對誤差分別為0.427、0.253、0.204和 0.181；平均相關(guān)系數(shù)分別為0.635、0.886、0.907和 0.914。

3.3 成像速度

針對用于油水兩相流流動過程成像的實時性要求，為提升算法的計算速度，使用TSVD對系數(shù)矩陣R及SART算法過程進(jìn)行降維處理，其中截斷系數(shù)q的選取對重建結(jié)果有著較大影響。較大的截斷系數(shù)可以保留系數(shù)矩陣中的主要奇異值，重建精度相對較高，但矩陣維度較大，重建過程較為耗時；較小的截斷系數(shù)可以減小矩陣維度、加快計算速度，但其重建精度將有明顯下降。

圖 5 為不同截斷參數(shù)q下的圖像重建結(jié)果的相對誤差和計算時間曲線。為分析不同探頭數(shù)目對成像結(jié)果的影響，分別采用16、32、48和64共4組超聲換能器數(shù)目進(jìn)行圖像重建。圖中，所計算的定量評估指標(biāo)為圖 3 中的模型1至模型4計算結(jié)果的平均值。

圖5中，隨著探頭數(shù)目增加，重建精度會明顯提高，但相應(yīng)計算時間也會增加。因此在實際測試時需要根據(jù)測試目標(biāo)的要求和測試實時性要求合理選擇探頭數(shù)量。此外，在確定探頭數(shù)目下，隨著截斷系數(shù)p的增大，重建誤差明顯減小，計算時間顯著增大；當(dāng)截斷系數(shù)較小時，重建誤差較高，而計算時間較短。為平衡計算時間與重建誤差之間的矛盾關(guān)系，采用L-曲線法[14]對截斷參數(shù)q進(jìn)行優(yōu)化，L曲線分布及優(yōu)化結(jié)果如圖 6 所示。經(jīng)優(yōu)化，截斷參數(shù)q選取為0.4，相應(yīng)的不同算法計算時間對比如表 2 所示。

圖5 不同截斷參數(shù)下誤差及時間Fig.5 Errors and speeds at different truncation parameter values

圖6 L曲線分析結(jié)果Fig.6 Results of L-curve analysis

表2 不同算法成像時間對比Table 2 Computing time comparison among different methods s

表2中，F(xiàn)BP算法成像速度最快，TV次之，傳統(tǒng)SART算法成像速度較慢，不能做到實時在線成像。SART-GPR算法的成像速度相比傳統(tǒng)SART算法提升約9倍，其成像速率達(dá)到35幅/s，可以滿足實時成像需求。

4 總結(jié)

針對采用陣列式超聲傳感器實現(xiàn)以油水兩相流為代表的弱聲阻抗比兩相介質(zhì)分布重建問題，提出一種基于高斯回歸預(yù)測的高分辨率相介質(zhì)分布圖像重建方法。在獲取采用超聲連續(xù)波激勵的方式的邊界測量數(shù)據(jù)基礎(chǔ)上，提取邊界測量聲壓幅值的衰減信息。在逆問題成像算法方面，通過采用能夠同時兼顧重建速度和重建精度的基于截斷奇異值分解的同步代數(shù)迭代重建算法，并通過構(gòu)建高斯回歸模型，將低分辨率重建結(jié)果轉(zhuǎn)化為高分辨率圖像。

通過對仿真模型測試數(shù)據(jù)的分析，所提出的算法成像精度及圖像純凈度相比傳統(tǒng)方法明顯提高，算法平均相對誤差減小至25.3%，平均相關(guān)系數(shù)提升至88.6%。經(jīng)TSVD降低維度后，算法成像速率達(dá)到35幅/s，滿足實時成像需求。

未來工作需要在構(gòu)建超聲信號激勵與采集系統(tǒng)，對所提出的算法進(jìn)行實驗驗證；并進(jìn)一步探索超聲層析成像技術(shù)實現(xiàn)油水兩相流介質(zhì)分布高效重建的方法。