摘要

首先介紹了成像的物理基礎(chǔ)，并對(duì)重建算法的各種方法進(jìn)行了簡(jiǎn)要闡述并對(duì)ART算法及其改進(jìn)算法進(jìn)行了主要說(shuō)明。在之前學(xué)者對(duì)成像原理的研究下，闡述了應(yīng)力波層析成像技術(shù)應(yīng)用于木材缺陷檢測(cè)當(dāng)中的原理及測(cè)試方法，并通過(guò)先驗(yàn)的介質(zhì)模型及實(shí)際模型測(cè)試，分別用ART算法和改進(jìn)后的重建算法得到木材缺陷圖像，通過(guò)比較兩種方法重建的圖像清晰度，驗(yàn)證改進(jìn)后算法的合理性與可行性。

關(guān)鍵詞 層析成像;圖像重建算法;ART;改進(jìn)的ART;木材檢測(cè);應(yīng)力波

中圖分類號(hào) S126 "文獻(xiàn)標(biāo)識(shí)碼 A "文章編號(hào) 0517-6611（2014）32-11585-03

Wood Stresswave Detection Based on Improved Tomography Image Reconstruction Algorithim of ART

LIU Jiaxin， CAI Yanli， YUE Qiang

（College of Machinery Electricity of Northeast Forestry University， Harbin， Heilongjiang 150040）

Abstract This paper firstly introduces the physical basis of imaging， "has a briefly introduction of various methods of reconstruction algorithms "and the ART algorithm and its improved algorithm are illustrated. On this basis， combining the above study of imaging principle， expounds the principle and test method of stresswave tomography technique which is applied to the wood defects detection， and through the medium model and actual test model， "the image definition comparison of two methods of reconstruction verify the rationality and feasibility of this algorithm respectively by the reconstruction algorithm and improved ART algorithm to get the wood defect images.

Key words "Tomography; Image reconstruction algorithm; ART; Improved ART; Wood detection; Stresswave

隨著現(xiàn)代的人類科技進(jìn)步、文化程度和生活水平的提高，人們對(duì)木材需求量如房子建造、一次性筷子使用等日益增多。但是由于對(duì)森林資源的過(guò)度砍伐，沒(méi)有規(guī)律地使用，導(dǎo)致全球森林日漸減少。在阻礙國(guó)家生產(chǎn)力水平的因素中，木材的供需沖突成為第一大因素，尤其是對(duì)于森林資源均衡不一的國(guó)家中，我國(guó)的森林資源是相對(duì)貧乏的。因此，森工企業(yè)要想滿足自身需求，就必須要想到解決提高木材的生產(chǎn)水平及木材質(zhì)量的辦法。

與其他無(wú)損檢測(cè)方法相比，木材應(yīng)力波無(wú)損檢測(cè)有很多優(yōu)勢(shì)，如對(duì)人體傷害小，容易操作，設(shè)備體積小且攜帶方便等。在木材應(yīng)力波無(wú)損檢測(cè)技術(shù)跟隨時(shí)代，向數(shù)字化、系統(tǒng)化、直觀化的方向發(fā)展中，應(yīng)力波與CT技術(shù)的結(jié)合占據(jù)了主導(dǎo)地位。為此，筆者開(kāi)展了基于改進(jìn)的層析成像圖像重建ART算法的木材應(yīng)力波檢測(cè)研究。

1 "層析成像簡(jiǎn)介及數(shù)學(xué)基礎(chǔ)

層析成像也稱為圖像重建，是利用某種射線源從研究“對(duì)象”外部用檢測(cè)設(shè)備所獲得的投影數(shù)據(jù)，依照物理和數(shù)學(xué)公式，利用computer來(lái)反演“對(duì)象”內(nèi)部未知的某種物理量的分布，生成圖像，重建“對(duì)象”的內(nèi)部特征[1]。

CT技術(shù)通過(guò)重建木材內(nèi)部圖像來(lái)對(duì)其進(jìn)行檢測(cè)從而判斷或確知有否缺陷，從而提高木材利用率。以射線理論為基礎(chǔ)的波速走時(shí)成像是應(yīng)力波層析成像的特點(diǎn)，是近幾年來(lái)發(fā)展起來(lái)的，一個(gè)木材檢測(cè)方面比較新的成像方法。

所有的層析成像的數(shù)學(xué)問(wèn)題是基本一致的，但是不同的實(shí)際問(wèn)題大都帶有所屬學(xué)科領(lǐng)域的特點(diǎn)。反演過(guò)程即是一個(gè)利用數(shù)學(xué)的一些方法來(lái)求逆的過(guò)程。反演問(wèn)題中存在兩個(gè)非常顯著的數(shù)學(xué)問(wèn)題：一是不適定問(wèn)題;二是非線性問(wèn)題。面對(duì)上述兩個(gè)關(guān)鍵問(wèn)題，歷代數(shù)學(xué)家和工程師們積極考慮各種解決辦法，截至目前，已經(jīng)有了各種各樣的CT，其算法和反演線性模型被提出且在工程中被應(yīng)用[2]。1917年Radon發(fā)表的關(guān)于積分變換和逆變換的數(shù)學(xué)論文后來(lái)被人們發(fā)現(xiàn)，并作為了Hounsfield和Cormack層析成像開(kāi)創(chuàng)性工作的理論基礎(chǔ)。Radon變換及其逆變換成為了CT層析成像技術(shù)得以實(shí)現(xiàn)的數(shù)學(xué)基礎(chǔ)[3]。Radon變換從某種意義上來(lái)說(shuō)是一種泛函算子，當(dāng)這種泛函算子作用在一個(gè)函數(shù)上的時(shí)候，產(chǎn)生另外的一個(gè)實(shí)數(shù)。根據(jù)投影值來(lái)重建對(duì)象內(nèi)部某物理量的分布圖像即是層析成像的反演過(guò)程，這也是Radon逆變換理論化公式的實(shí)現(xiàn)過(guò)程。

對(duì)于任意的函數(shù)f，記其對(duì)應(yīng)的Radon變換為Rf，其中Rf是f沿一直線l的線積分，將這個(gè)積分值亦稱投影值。函數(shù)的Radon變換定義為：

[Rf]（l，θ）=∫Lf（r，）ds（|l|≤E，0≤θ≤π） " " " " " （1）

其逆變換算子R-1，它滿足R-1Rf=f。對(duì)于逆變換的解，有如下公式：

[R-1P]（r，）=12π2∫π0 ∫E-E1rcos（θ-）-1Pl（l，θ）dldθ

（2）

式中，P（l，θ）=Rf（l，θ），Pl（l，θ）表示P（l，θ）關(guān)于l的偏導(dǎo)數(shù)。

2 "圖像重建算法的基本算法

層析成像從本質(zhì)上是根據(jù)對(duì)函數(shù)的某種積分資料，多數(shù)情況下是沿著一系列射線路徑的積分資料來(lái)反推函數(shù)分布的一種算法或者數(shù)據(jù)處理方法。圖像重建算法主要是投影重建方法，大部分投影重建算法是根據(jù)投影數(shù)據(jù)進(jìn)行重建的，而這些算法都是以Radon的理論為基礎(chǔ)發(fā)展起來(lái)的，目前的重建算法中，依據(jù)其特點(diǎn)大致可以分為兩類[4]。第一類是變換重建的方法，第二類是級(jí)數(shù)展開(kāi)重建法。以下介紹的重建算法分別為分析法、代數(shù)重建算法、聯(lián)合迭代算法、阻尼最小二乘法及改進(jìn)的ART算法。

2.1 分析法

在該方法類別中，Radon變換方法、Fourier方法等是主要的方法。該方法要求投影的數(shù)據(jù)完整全面，足夠精確且射線路徑最好為直線的前提下，才能夠比較準(zhǔn)確地重建對(duì)象內(nèi)部圖像，醫(yī)學(xué)CT常常采用此類變換法進(jìn)行圖像重建。但是該方法抗噪能力差，通常應(yīng)力波層析成像不采用該方法[1]。

2.2 代數(shù)重建法

代數(shù)重建法（Algebraic Reconstruction Techniques，簡(jiǎn)稱ART）是一種迭代的級(jí)數(shù)展開(kāi)法，是最早提出并為人們所熟悉的一種迭代型算法[5]。此算法是在一定的生成函數(shù)下，選擇一組列向量使得該序列逐漸逼近滿足一定最優(yōu)化準(zhǔn)則和約束條件的圖像向量估計(jì)值。

2.2.1

迭代法的基本公式。在這里采用網(wǎng)格劃分的方法，將待重建區(qū)域分為M×N個(gè)格子即像素，每個(gè)格子選擇一個(gè)基圖像bj，所有的J=M×N個(gè)基圖像構(gòu)成一列基向量{b1，b2，…，bJ}，它的線性組合可逼近重建函數(shù)f[1]，即用

f*=∑Jj=1（xj×bj） " " " " " " " " " " " （3）

逼近f，其中f*表示在第j個(gè)像素內(nèi)的平均值。{xj}（j=1，2，…，J）稱為圖像向量并用X表示。引入投影觀測(cè)值yi，并用它代替f*的Radon變換且令rij=ki×bj，則得到：

yi=∑Jj=1rij×xj " " "（4）

令R表示由元素rij所定義的I×J維矩陣即投影矩陣，Y表示由元素yi定義的一維向量，成為測(cè)量向量。同時(shí)e為I維列向量，它的第i個(gè)分量ei是式（4）兩端之差，式（4）可改寫(xiě)為：

Y=R×X+e " " " " " " " " " " " "（5）

式（5）即為迭代重建法基本公式。待重建圖像的估計(jì)值f在有一定的優(yōu)化準(zhǔn)則和約束條件下（以取代誤差向量求解線性方程組）計(jì)算出圖像向量的解X*，再通過(guò)公式（4）、（3）得到。

2.2.2

迭代算法在木材檢測(cè)圖像重建方面的計(jì)算公式。該算法是根據(jù)木材檢測(cè)各方面問(wèn)題考慮在內(nèi)的情況下提煉出來(lái)的算法，其基本的思想是將一個(gè)初始值先賦值給需要重建的場(chǎng)的反演物理量，然后開(kāi)始計(jì)算投影值殘差，再將其一個(gè)個(gè)沿射線方向均勻地反投影，進(jìn)入迭代計(jì)算循環(huán)來(lái)校正重建圖像，直到滿足精度要求才停止。通過(guò)求解方程組（5），即可得到符合該方向的ART的迭代公式為：

S（n+1）=S（n）+Ti-[li，S（n）][li，li]lTi " " " " " " " （6）

式中，S（n）是求解慢度向量;li代表距離矩陣中，第n次的迭代值，第i條射線的分量;[]代表向量的內(nèi)積;T i是第i條射線的實(shí)際時(shí)間值。ART算法由于大量的迭代計(jì)算，較高的儲(chǔ)存空間等需要其能快速地收斂。

該重建算法用于圖像重建是存在兩個(gè)缺點(diǎn)：①圖像在投影計(jì)算時(shí)會(huì)產(chǎn)生重建場(chǎng)的嚴(yán)重噪聲;②需要很大迭代次數(shù)，重建效率不高。

2.3 聯(lián)合迭代算法

聯(lián)合迭代算法（Simultaneous Iterative Reconstruction Techniques，簡(jiǎn)稱SIRT）與ART相似，不同的地方在于ART是每次的修正都只考慮一條射線，而SIRT方法中，一個(gè)方格的平均修正值的確定是將方格內(nèi)通過(guò)的所有的射線修正值都考慮在內(nèi)進(jìn)行計(jì)算的。SIRT的這種迭代算法的誤差修正計(jì)算方法的改進(jìn)，可以有效地抑制這些缺點(diǎn)，其迭代公式如下[6]：

S（n+1）j=S（n）j+ΔS（n）j

ΔS（n）j=1M∑Mi=1[lijTi-∑mk=1likS（n）k）/∑mk=1l2ik] （7）

式中，S（n）j為第j個(gè)像素第n次迭代后的慢度值;Mj為通過(guò)第j個(gè)像素的射線條數(shù);ΔS（n）j為第j個(gè)像素在第n次迭代中求得的慢度修正量;lij為第i根射線通過(guò)第j個(gè)像素的射線長(zhǎng)度;lik為第i根射線在第k個(gè)像素中的長(zhǎng)度;Ti為第i條射線的觀測(cè)走時(shí);m為離散的像素總數(shù);S（n）k為第k像素在n次迭代后的慢度值。通過(guò)以上公式，根據(jù)初始的慢度模型及觀測(cè)走時(shí)即可求出所有像素慢度值。

2.4 阻尼最小二乘法

阻尼最小二乘法是一種基于最小二乘法的投影方法，其利用Lanczos方法進(jìn)行方程求解[1]。其主要是利用迭代的方法來(lái)求解如下公式的一種算法：

AX=B

min‖AX=B‖2 （8）

式中，A為大型稀疏矩陣。

2.5 改進(jìn)的ART算法

傳統(tǒng)的迭代重建算法中，各個(gè)成像單元慢度迭代初值的選取，只依據(jù)先驗(yàn)知識(shí)給出一個(gè)統(tǒng)一的值作為初值。而迭代的過(guò)程中，需要將每條射線走時(shí)殘差均勻分給各個(gè)網(wǎng)格單元。但在實(shí)際過(guò)程中，射線的走時(shí)非常重要，而影響其精確度的是由缺陷單元引起的，若能提前預(yù)知缺陷單元網(wǎng)格的可能，在迭代初值選擇中區(qū)別對(duì)待，便可以提高算法精確度。具體算法如下。

2.5.1

確定缺陷單元網(wǎng)格。

vi=∑jaij/ti （9）

v=∑ni=1vi/n （10）

Sv=∑i（vi-v）2/（n-1）

（11）

式中，vi為第i條射線波速;v為波速均值;sv為波速標(biāo)準(zhǔn)差;n為射線總數(shù)。

假設(shè)同一測(cè)區(qū)射線波速服從正態(tài)分布，則可得到臨界波速公式：

V0=v-λgs*v/n （12）

若某條射線波速v小于臨界波速V0，則認(rèn)為射線穿過(guò)缺陷單元，并假定它穿過(guò)的單元格為缺陷格，但其實(shí)也有正常單元格。先求出一個(gè)激發(fā)點(diǎn)到其他所有接收點(diǎn)的射線的臨界波速，并找出穿過(guò)缺陷單元的射線，確定一個(gè)缺陷單元集合，其中該單元集合也有正常單元格。再根據(jù)式（12）計(jì)算其他所有接收點(diǎn)作為激發(fā)點(diǎn)時(shí)射線的臨界波速，同時(shí)判斷出正常無(wú)缺陷單元集合，求兩個(gè)集合交集，即可求得最終的缺陷單元格集合。

2.5.2

確定迭代初值。把沒(méi)有穿過(guò)缺陷的射線長(zhǎng)度除以走時(shí)值即作為正常射線平均波速，把最小的正常射線波速賦值給預(yù)判的缺陷單元。

2.5.3

迭代計(jì)算松弛因子的選取。公式為：

fq∧，i+1j=fq∧，ij+μgaij∑mj=1a2ij（ti-tq∧i）

（13）

式中，fq∧，ij為q次迭代后第i條射線對(duì)第j個(gè)單元格的波慢估值;當(dāng)j為缺陷時(shí)，取0≤μ=μ1≤1，當(dāng)j為正常單元時(shí)，取0lt;μ=μ2lt;1。完成一次迭代后判斷是否達(dá)精度要求，若滿足，停止，否則按式（13）繼續(xù)迭代。

3 "ART算法及其改進(jìn)算法的實(shí)現(xiàn)與比較

3.1 試驗(yàn)?zāi)Ｐ偷臉?gòu)造及數(shù)據(jù)采集

該試驗(yàn)采用Fakopp Microsecond Time 應(yīng)力波木材檢測(cè)儀，截取一段圓柱形的木材，并在測(cè)試木材上均勻插入8個(gè)測(cè)試傳感器，在這里測(cè)量的誤差中，忽略傳感器針腳釘入深度的影響。連接好設(shè)備后，開(kāi)始測(cè)量，并記錄下數(shù)據(jù)。

數(shù)據(jù)處理，先是使用C程序計(jì)算出距離矩陣，然后利用ART及改進(jìn)的ART算法Matlab程序計(jì)算各像元慢度值。

3.2 ART和SIRT成像算法實(shí)現(xiàn)及結(jié)果

對(duì)得到的慢度向量，使用軟件SURFER 8.0進(jìn)行直接畫(huà)圖。實(shí)物圖像見(jiàn)圖1，ART及改進(jìn)的ART算法的圖像如圖2、3所示。

圖1 實(shí)物圖像

圖2 ART算法圖像

圖3 改進(jìn)的ART算法圖像

比較圖1、2的重建效果可以看出，木材中存在的缺陷區(qū)域都可以由兩種算法檢測(cè)出，但是重建圖像缺陷區(qū)域的形狀和大小存在誤差。改進(jìn)的ART算法有效地抑制了ART算法帶來(lái)的模糊性，因此，其重建圖像精度也比ART高。

參考文獻(xiàn)

[1]

閆再興.基于應(yīng)力波原木內(nèi)部缺陷二維圖像重建的初步研究[D].哈爾濱：東北林業(yè)大學(xué)，2007.

[2] 章毓晉.圖像處理和分析[M].北京：清華大學(xué)出版社，1999.

[3] HERMAN G T.Image reconstrution from projections：The fundamentals of computerized tomography[M].New York：Academic，1980.

[4] 王學(xué)勝.超聲層析技術(shù)中射線追蹤方法的研究與應(yīng)用[D].北京：中國(guó)地質(zhì)大學(xué)（北京），2005.

[5] 沈?qū)?工業(yè)CT窄角扇束的迭代圖像重建算法研究[D].重慶：重慶大學(xué)，2002.

[6] 方璽.混凝土結(jié)構(gòu)的超聲探傷信號(hào)處理及二維CT軟件設(shè)計(jì)[D].武漢：武漢理工大學(xué)，2006.