基于壓縮感知的EMT探傷和圖像重建?

莊雨欣何 敏?孫海濱曾星星周 明馮肖維

(1.上海海事大學(xué)物流工程學(xué)院，上海 201306；2.中國船舶集團(tuán)有限公司第七○四研究所，上海 200031；3.上海寶信軟件股份有限公司，上海 201900)

近年來，隨著港口的高速發(fā)展，港口機(jī)械時(shí)常處于連續(xù)高強(qiáng)度的超負(fù)荷運(yùn)行狀態(tài)。這導(dǎo)致部分港機(jī)的金屬結(jié)構(gòu)因疲勞或者老化而出現(xiàn)不同類型的缺陷，其中裂紋是最常見的[1－4]。根據(jù)調(diào)查發(fā)現(xiàn)，這些裂紋的寬度一般在1 mm左右，而長度則在幾毫米到幾厘米不等。

基于電磁感應(yīng)原理的電磁層析成像(electromagnetic tomography，EMT)技術(shù)是二十年前發(fā)展起來的一種過程層析成像技術(shù)[5]。在港機(jī)金屬結(jié)構(gòu)裂紋的探測中，EMT具有非侵入、無輻射、成本低、可視化等[6－7]優(yōu)點(diǎn)。要實(shí)現(xiàn)EMT探傷和圖像重建關(guān)鍵在于EMT逆問題的求解，而求解EMT逆問題時(shí)由于受到硬件系統(tǒng)的限制以及外部檢測環(huán)境的影響，使得到的實(shí)際測量數(shù)據(jù)個(gè)數(shù)遠(yuǎn)少于要求解的電磁特征參數(shù)的個(gè)數(shù)，從而導(dǎo)致逆問題具有嚴(yán)重的欠定性。因此，亟需找到一種有效可行的方法來解決這個(gè)問題，從而實(shí)現(xiàn)高質(zhì)量的圖像重建。

文獻(xiàn)[8]中，作者通過利用數(shù)學(xué)方法對靈敏度矩陣作降維處理，并結(jié)合人群搜索算法降低靈敏度矩陣條件數(shù)，以此來改善逆問題的欠定性。在文獻(xiàn)[9]中，作者采用了統(tǒng)計(jì)概率的方法來處理逆問題的欠定性。這些方法的使用在一定程度上都有助于改善逆問題的欠定性，但也存在一些有待優(yōu)化的地方。比如，文獻(xiàn)[8]中作者用數(shù)學(xué)方法直接對靈敏度矩陣進(jìn)行降維處理，這很可能會(huì)導(dǎo)致有用信息的丟失，從而影響后續(xù)的圖像重建。

為有效解決這個(gè)問題，本文提出一種基于壓縮感知(compressive sensing，CS)原理的探傷和圖像重建方法。壓縮感知是近年來一種新發(fā)展起來的數(shù)據(jù)采樣技術(shù)，其過程包括選取恰當(dāng)?shù)南∈枳儞Q基對原始信號進(jìn)行稀疏表示，將稀疏變換后的原始信號投影到構(gòu)建的測量矩陣上進(jìn)而得到壓縮采樣信號。由上述過程可知，壓縮感知技術(shù)可以減少對測量信息采樣數(shù)量的需求，弱化逆問題求解所需的求解條件，從而有效改善EMT逆問題的欠定性。

1 EMT探傷系統(tǒng)的數(shù)學(xué)模型

通常，EMT探傷系統(tǒng)的數(shù)學(xué)模型可以用如式(1)所示的非線性方程來表示:

式中:v表示傳感器獲得的測量數(shù)據(jù)，x表示金屬結(jié)構(gòu)的電導(dǎo)率分布，F(xiàn)則是由系統(tǒng)決定的非線性函數(shù)。

式(1)所示的非線性方程往往無法直接求解，故對它進(jìn)行線性化處理，可得到如式(2)所示的線性方程[10]。

式中:x0表示空場電導(dǎo)率分布，表示靈敏度矩陣。

將式(2)進(jìn)一步簡化，則可得到EMT系統(tǒng)的數(shù)學(xué)模型，如式(3)所示。

式中:U為測量數(shù)據(jù)向量，G為電導(dǎo)率空間分布向量，為靈敏度矩陣，本文采用模型擾動(dòng)法獲取靈敏度矩陣，并用式(4)所示的計(jì)算公式進(jìn)行求解。

式中:uful(j)，uemp(j)，u(j)｜k分別表示滿場(裂紋覆蓋整個(gè)物場區(qū)域)，空場(無裂紋)以及第k個(gè)像素點(diǎn)為裂紋時(shí)，第j個(gè)檢測線圈內(nèi)的感應(yīng)電壓值。

2 基于壓縮感知的EMT探傷信號處理

2.1 壓縮感知原理

壓縮感知技術(shù)能夠在以遠(yuǎn)低于Nyquist采樣頻率采集特征信息的情況下，實(shí)現(xiàn)對原始信號的重構(gòu)。壓縮感知理論主要包括信號的稀疏表示、編碼測量和信號重構(gòu)三個(gè)方面[11]。

2.1.1 信號的稀疏表示

實(shí)現(xiàn)壓縮感知的先決條件是信號必須為稀疏的或在某種變換域下近似稀疏。而常見的自然信號一般都不是稀疏的，因此有必要對其進(jìn)行稀疏表示，以使信號達(dá)到稀疏或者近似稀疏。

對于一個(gè)非稀疏的N維離散信號u，其稀疏表示如式(5)所示。

式中:ψ為N×N維的稀疏變換基，x為離散信號u的稀疏表示形式。常見的稀疏變換基有[12]:離散傅里葉變換(DFT)基、離散余弦變換(DCT)基、離散正弦變換(DST)基等。

2.1.2 信號編碼測量

編碼測量的核心是測量矩陣的設(shè)計(jì)，它關(guān)系著信號的壓縮采樣。其具體過程為:原始信號經(jīng)過稀疏處理后，被投影到測量矩陣?M×N上，從而獲得原始信號u的壓縮采樣形式。該過程的數(shù)學(xué)表達(dá)如式(6)所示。

式中:λ＝?ψ，稱λM×N為傳感矩陣，y M×1為壓縮采樣信號，并且M?N。由此可知，信號的編碼測量處理弱化了式(6)的求解條件。

由于M?N，方程的個(gè)數(shù)遠(yuǎn)少于未知數(shù)的個(gè)數(shù)，故式(6)是一個(gè)欠定方程，求出的解也不唯一。而壓縮感知理論表明[13]，當(dāng)測量矩陣?與稀疏變換基ψ不相關(guān)，且它們的乘積傳感矩陣λ滿足一定的限制條件即約束等距特性(Restricted Isometry Property，RIP)時(shí)，式(6)存在唯一解。服從高斯隨機(jī)分布的矩陣與稀疏基的相關(guān)性很小且較高概率滿足RIP性質(zhì)[14]，因此可通過構(gòu)建一個(gè)服從高斯隨機(jī)分布的測量矩陣來實(shí)現(xiàn)壓縮采樣，然后結(jié)合相關(guān)的重構(gòu)算法來求解式(6)。

2.1.3 信號重構(gòu)

信號重構(gòu)是壓縮感知理論的核心，指的是由壓縮采樣信號y重構(gòu)原始信號u的過程。一般轉(zhuǎn)換成求解最小L0范數(shù)的優(yōu)化問題，如式(7)所示。

式中:xopt表示式(6)的最優(yōu)解，uopt表示重構(gòu)得到的原始信號最佳恢復(fù)。

式(7)的直接求解屬于NP-hard問題[15]。研究人員提出了一些求次最優(yōu)解的算法，主要有:最小L1范數(shù)法、匹配追蹤系列算法以及迭代閾值法等[16－17]。

2.2 基于壓縮感知的EMT探傷信號處理

基于壓縮感知的EMT探傷信號處理包括稀疏表示、編碼測量和重構(gòu)。

①待重構(gòu)信號的稀疏表示

在EMT系統(tǒng)的數(shù)學(xué)模型中，我們將所求的電導(dǎo)率G視為壓縮感知理論中待重構(gòu)的原始信號，對其進(jìn)行稀疏處理，如式(8)所示。

式中:x為電導(dǎo)率G的稀疏表示形式。

②信號的編碼測量

要實(shí)現(xiàn)EMT探傷中信號的編碼測量關(guān)鍵在于靈敏度矩陣能否滿足壓縮感知中測量矩陣的特性。根據(jù)前述分析可知，高斯隨機(jī)矩陣是一種符合條件的性能較好的測量矩陣。對于一個(gè)M×N的高斯隨機(jī)矩陣，其元素滿足期望為0，方差為1/M的高斯分布。而靈敏度矩陣的定義指的是檢測線圈感應(yīng)電壓變化與引起檢測線圈感應(yīng)電壓變化的成像空間中單個(gè)像素上電磁參數(shù)變化之比。根據(jù)此定義可知，靈敏度矩陣是被測區(qū)域所有各自相互獨(dú)立的像素點(diǎn)共同作用的結(jié)果，每個(gè)像素點(diǎn)只是貢獻(xiàn)了很小的一部分。文獻(xiàn)[18]指出，若一個(gè)量是由許多微小的獨(dú)立隨機(jī)因素影響的結(jié)果，則這個(gè)量近似服從高斯分布。由此可知，靈敏度矩陣也近似服從高斯分布。而通過計(jì)算課題組設(shè)計(jì)的如圖1所示的用于EMT探傷的傳感結(jié)構(gòu)的靈敏度矩陣的期望和方差也發(fā)現(xiàn)，其結(jié)果與高斯隨機(jī)矩陣的期望和方差非常接近，這進(jìn)一步表明靈敏度矩陣近似服從高斯分布。

圖1 七線圈傳感結(jié)構(gòu)圖

因此，本文將EMT探傷數(shù)學(xué)模型中的靈敏度矩陣S近似為壓縮感知理論中的高斯隨機(jī)測量矩陣，并結(jié)合式(8)，可以得到EMT探傷中信號的編碼測量過程如式(9)所示。

式中:測量數(shù)據(jù)向量U的維數(shù)為M(M?N)，傳感矩陣λ＝Sψ。

信號的編碼測量過程使得經(jīng)稀疏處理的電導(dǎo)率信號G，在測量矩陣S的作用下，得到了更低維數(shù)的測量向量U。這表明，壓縮感知可以弱化求解條件，從而改善EMT逆問題的欠定性。這也意味著EMT系統(tǒng)能夠以較少的測量數(shù)據(jù)實(shí)現(xiàn)檢測和圖像重建過程。

此外，為了使得靈敏度矩陣S與稀疏變換基ψ具有很小的相關(guān)性以及它們的乘積傳感矩陣λ較高概率滿足RIP性質(zhì)，以隨機(jī)順序去獲取測量數(shù)據(jù)可以使靈敏度矩陣具有較高的隨機(jī)性[19－20]。因此，在測量時(shí)先以隨機(jī)順序去采集各傳感器的檢測數(shù)據(jù)，然后將獲得的測量數(shù)據(jù)通過歸一化公式(4)計(jì)算直接得到隨機(jī)的靈敏度矩陣。這樣所得到的靈敏度矩陣就可近似為高斯隨機(jī)測量矩陣，從而實(shí)現(xiàn)壓縮感知中即信號的編碼測量過程。

③信號重構(gòu)

EMT探傷中的信號重構(gòu)其實(shí)質(zhì)就是對式(9)進(jìn)行求解。本文采用了兩種不同的方法來間接地求解最小L0范數(shù)問題。第一種是將L0范數(shù)松弛到最小L1范數(shù)的優(yōu)化求解，如式(10)所示。

式中:xopt表示式(9)的最優(yōu)解，Gopt表示重構(gòu)得到的最佳的電導(dǎo)率恢復(fù)。

另一種則是利用稀疏逼近的方法來間接地解決最小L0范數(shù)問題。利用匹配追蹤的思想，通過逐步地從傳感矩陣λ中選擇恰當(dāng)?shù)牧衼硐∈璞平鼫y量信號U，進(jìn)而逐步地確定索引集I?{1，…，N}，如式(11)所示。

3 基于壓縮感知的EMT探傷圖像重建算法

文獻(xiàn)[21－22]表明，最小L1范數(shù)法和匹配追蹤系列算法對于信號不是特別稀疏且測量數(shù)據(jù)包含一定的噪聲的情況比較有效，而迭代閾值法對初值敏感，很難直接用于信號的重構(gòu)，通常需與其他方法結(jié)合使用。因此，本文選取了兩種算法——基于最小L1范數(shù)的迭代重加權(quán)最小二乘法和基于匹配追蹤的壓縮采樣匹配追蹤算法進(jìn)行圖像重建。

3.1 迭代重加權(quán)最小二乘法

基于L1范數(shù)的迭代重加權(quán)最小二乘(iteratively re-weighted least squares，IRLS)法是一種稀疏算法。

該算法首先用加權(quán)L2范數(shù)替代式(10)中的目標(biāo)函數(shù)，如式(12)所示。

接著，使用Lagrange乘法得到每次迭代的x(τ)＝D N ST(SD N ST)－1U，其中D N為N×N的對角陣，其元素為并且ε(τ)＝ε(τ－1)/10。具體步驟如下:

①初始化:x0為U＝Sψx的最小二乘解，ε＝1，τ＝1；

②W(τ)＝{[x(τ－1)]2＋ε}p/2－1，v(τ)＝1./W(τ)，D(τ)＝diag[v(τ)，0]；

③x(τ)＝D N ST(SD N ST)－1U；

④若‖x(τ)－x(τ－1)‖2<ε(1/2)/100，ε＝ε/10；

⑤若ε>10－8，τ＝τ＋1。

其中“1./w(τ)”表示用1與矩陣W(τ)中各元素相除。

3.2 壓縮采樣匹配追蹤算法

壓縮采樣匹配追蹤(compressive sampling matching pursuit，CoSaMP)算法是經(jīng)典匹配追蹤(MP)算法的改進(jìn)。該算法在每次迭代過程中從測量矩陣?yán)镎页雠c信號相關(guān)性最大的2K列，同時(shí)剔除一些不相關(guān)的列以使得算法能夠快速的收斂。其具體步驟如下:

①初始化:初值x(0)＝0，初始?xì)埩縱(0)＝U，迭代次數(shù)n＝0，索引集I(0)為空集；

②計(jì)算殘量與測量矩陣的相關(guān)性:z＝(Sψ)Tv；

③從z中找出2K個(gè)最大值其相應(yīng)的位置集合Ω，然后更新索引集I(n)＝I(n－1)∪Ω；

上述兩種算法各有優(yōu)缺點(diǎn)，IRLS算法作為一種凸優(yōu)化算法，優(yōu)勢在于能夠在較少測量數(shù)據(jù)的情況下準(zhǔn)確重構(gòu)原信號，但耗時(shí)較長[23]；此外，IRLS算法涉及到多矩陣乘法的逆運(yùn)算，所以計(jì)算復(fù)雜度較高。而CoSaMP算法作為貪婪算法中的一種，其迭代求解過程與IRLS算法相比，計(jì)算效率大大提高，并且計(jì)算復(fù)雜度也較低，但它對測量數(shù)據(jù)的需求要高于IRLS算法。這兩種算法應(yīng)用于EMT探傷圖像重建的效果具體如何，將通過實(shí)驗(yàn)來驗(yàn)證。

4 基于壓縮感知的EMT探傷圖像重建仿真

本文在Ansoft Maxwell 16.0環(huán)境中建立EMT探傷的三維仿真模型。

4.1 稀疏基的選擇

原信號經(jīng)不同稀疏變換基處理后得到的稀疏信號有所不同，從而使得圖像重建的效果也不一樣。本文比較了三種不同的稀疏變換基(DFT基、DCT基和DST基)對原信號處理后，分別使用IRLS算法和CoSaMP算法進(jìn)行圖像重建的效果。

仿真中設(shè)置了30 mm和15 mm兩種不同長度的裂紋，在加入隨機(jī)噪聲的情況下進(jìn)行了仿真實(shí)驗(yàn)，其中裂紋的寬都為1 mm，深度為10 mm，仿真結(jié)果如圖2所示。

上述結(jié)果表明，信號經(jīng)不同稀疏基處理后圖像重建得到的結(jié)果完全不同，由圖2(a)可以得出，DFT基稀疏處理后得到的信號用于IRLS算法的重建效果更好。由圖2(b)可以看出，基于DCT基的CoSaMP算法能夠準(zhǔn)確且高質(zhì)量的重建目標(biāo)裂紋圖像。

圖2 不同稀疏變換基下IRLS算法和CoSaMP算法圖像重建結(jié)果比較

下面具體分析上述結(jié)果的原因，基于L1范數(shù)的IRLS算法作為一種稀疏算法，要求信號盡可能的稀疏，這樣重建效果就越好。而DFT基是這三種變換基中各元素?cái)?shù)值最接近于零的，用它來處理原信號，可以使原信號達(dá)到最稀疏，所以基于DFT基的IRLS算法的重建效果要更好。而關(guān)于CoSaMP算法，其利用的是匹配追蹤的思想，通過不斷的將求解范圍縮小，從而求得最優(yōu)解。若選擇非常接近于零的變換基對原信號作稀疏處理，則在上述算法流程的第二步計(jì)算相關(guān)性時(shí)可能會(huì)直接出現(xiàn)相關(guān)性為0的情況，導(dǎo)致后續(xù)求解過程失敗。但為了使處理后的信號稀疏以降低求解的復(fù)雜度，變換基元素的數(shù)值也不能太大，否則會(huì)增加數(shù)據(jù)處理的難度。而在上述三種變換基中，更符合條件的只有DCT基。因此，基于DCT基的CoSaMP算法重建效果較好。

4.2 壓縮感知算法與傳統(tǒng)算法圖像重建的對比仿真

為驗(yàn)證壓縮感知算法(即基于DFT基的IRLS算法和基于DCT基的CoSaMP算法)的性能，本文選取了EMT領(lǐng)域中常用的兩種傳統(tǒng)算法進(jìn)行比較。一種是Tikhonov正則化，另一種是Landweber迭代法。其中Tikhonov正則化屬于非迭代算法，耗時(shí)短，引入基于2范數(shù)的約束項(xiàng)來縮小解空間以求得最佳近似解[24]。而Landweber迭代法是一種迭代算法，圖像重建質(zhì)量較好，但耗時(shí)較長、存在半收斂特性[25]。這兩種算法各有優(yōu)缺點(diǎn)，在EMT領(lǐng)域非常具有代表性。實(shí)驗(yàn)中同樣也設(shè)置了30 mm和15 mm兩種不同長度的裂紋進(jìn)行仿真，裂紋的寬也都為1 mm，深度為10 mm。仿真結(jié)果分別如圖3和圖4所示。

由圖3和圖4可以看出，傳統(tǒng)Landweber迭代法和Tikhonov正則化法圖像重建的結(jié)果中出現(xiàn)了不少干擾且重建裂紋形狀存在失真。而壓縮感知算法即IRLS算法和CoSaMP算法圖像重建得到的目標(biāo)裂紋圖像質(zhì)量較好，明顯優(yōu)于Landweber迭代法和Tikhonov正則化法。其原因在于被測對象裂紋尺寸太小，使得引起的磁場變化非常微弱，加之測量數(shù)據(jù)本身就含有噪聲，導(dǎo)致測量信噪比大大降低，而上述兩種傳統(tǒng)算法都是直接求解，并沒有對信號進(jìn)行相關(guān)的去燥處理，所以造成重建效果較差。而在應(yīng)用壓縮感知算法的時(shí)候首先就選取了相應(yīng)的變換基對信號進(jìn)行稀疏處理，這在一定程度上達(dá)到了去燥的效果，因此重建效果要更好。

圖3 裂紋長度為30 mm時(shí)不同算法圖像重建結(jié)果對比

圖4 裂紋長度為15 mm時(shí)不同算法圖像重建結(jié)果對比

上述仿真實(shí)驗(yàn)表明，在相同數(shù)量測量數(shù)據(jù)的情況下，基于壓縮感知的EMT圖像重建效果遠(yuǎn)好于傳統(tǒng)的EMT圖像重建方法。實(shí)際上，由于各種條件的限制，獲得的測量數(shù)據(jù)是有限的，遠(yuǎn)少于要求解的物場參數(shù)，從而導(dǎo)致了EMT逆問題嚴(yán)重的欠定性。所以在傳感器結(jié)構(gòu)確定的情況下，傳統(tǒng)的EMT重建方法由于欠定性的影響，造成較差的圖像重建質(zhì)量。而壓縮感知能夠弱化求解條件，在測量數(shù)據(jù)一定的情況下，通過壓縮感知就相當(dāng)于達(dá)到了增加測量數(shù)據(jù)數(shù)量的效果，從而改善了欠定性，實(shí)現(xiàn)高質(zhì)量的圖像重建。

5 實(shí)驗(yàn)驗(yàn)證

本文作者開發(fā)的7線圈EMT探傷系統(tǒng)如圖5所示。

圖5 EMT探傷實(shí)驗(yàn)系統(tǒng)

在實(shí)際實(shí)驗(yàn)時(shí)，同樣選取了15 mm和30 mm兩種不同長度的裂紋進(jìn)行圖像重建，四種算法的重建結(jié)果如圖6和圖7所示。

圖6 裂紋長度為30 mm時(shí)不同算法圖像重建的實(shí)驗(yàn)對比

圖7 裂紋長度為15 mm時(shí)不同算法圖像重建的實(shí)驗(yàn)對比

從圖6和圖7的實(shí)驗(yàn)結(jié)果可以看出，傳統(tǒng)的Landweber迭代法和Tikhonov正則化法重建的圖像中存在較多干擾，失真也較嚴(yán)重。而壓縮感知的IRLS算法與CoSaMP算法所得到的重建圖像含干擾很小，圖像重建質(zhì)量較好。為定量比較這四種算法的圖像重建質(zhì)量，本文用式(13)所示的計(jì)算公式分別計(jì)算出它們各自的相對誤差Error，結(jié)果如圖8所示。

圖8 裂紋長度為30mm和15mm時(shí)不同算法圖像重建實(shí)驗(yàn)的相對誤差

式中:G為電導(dǎo)率實(shí)驗(yàn)分布，G?為電導(dǎo)率設(shè)定分布。

由圖8可得，壓縮感知的IRLS算法和CoSaMP算法的圖像重建的相對誤差明顯小于傳統(tǒng)算法。這進(jìn)一步說明，壓縮感知算法的圖像重建質(zhì)量要好于傳統(tǒng)算法。此外，通過圖8還可以發(fā)現(xiàn)，CoSaMP算法的相對誤差略小于IRLS算法，表明CoSaMP算法圖像重建的精度稍高于IRLS算法。

6 結(jié)論

本文研究了基于壓縮感知原理的EMT探傷和圖像重建。壓縮感知可以弱化逆問題的求解條件，有效改善EMT逆問題的欠定性，從而能夠在較少測量數(shù)據(jù)的情況下完成圖像重建。在實(shí)際EMT探傷中，由于硬件實(shí)現(xiàn)的限制，獲得的測量數(shù)據(jù)是有限的，遠(yuǎn)少于要求解的物場參數(shù)。所以在傳感結(jié)構(gòu)確定的情況下，壓縮感知弱化求解條件的這種特性，就相當(dāng)于達(dá)到了增加測量數(shù)據(jù)的效果。這樣就避免了通過直接增加傳感器數(shù)量來增加測量信息造成數(shù)據(jù)處理更復(fù)雜，累計(jì)測量誤差更大，并且對硬件系統(tǒng)性能的要求也更高。

在對EMT探傷信號進(jìn)行處理后，本文引入了兩種適用于EMT探傷圖像重建的壓縮感知算法，一種是基于最小L1范數(shù)的IRLS算法，另一種則是基于匹配追蹤思想的CoSaMP算法。

仿真和實(shí)驗(yàn)結(jié)果均表明，基于壓縮感知的IRLS算法和CoSaMP算法的圖像重建質(zhì)量均好于EMT圖像重建中常用的Landweber迭代法以及Tikhonov正則化法。其中，CoSaMP算法的圖像重建精度又比IRLS算法略勝一籌。在實(shí)際應(yīng)用時(shí)，如果獲得的測量數(shù)據(jù)有限，且對耗時(shí)要求不高，則可以采用IRLS算法進(jìn)行重建。如果能夠獲得較多的測量數(shù)據(jù)，則選擇CoSaMP算法，以得到精度更高的重建圖像。

綜上，基于壓縮感知的EMT探傷圖像重建可以有效提高EMT探傷中裂紋圖像重建的質(zhì)量。