基于改進的層析成像圖像重建ART算法的木材應力波檢測研究

摘要

首先介紹了成像的物理基礎，并對重建算法的各種方法進行了簡要闡述并對ART算法及其改進算法進行了主要說明。在之前學者對成像原理的研究下，闡述了應力波層析成像技術應用于木材缺陷檢測當中的原理及測試方法，并通過先驗的介質模型及實際模型測試，分別用ART算法和改進后的重建算法得到木材缺陷圖像，通過比較兩種方法重建的圖像清晰度，驗證改進后算法的合理性與可行性。

關鍵詞 層析成像;圖像重建算法;ART;改進的ART;木材檢測;應力波

中圖分類號 S126 "文獻標識碼 A "文章編號 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

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

與其他無損檢測方法相比，木材應力波無損檢測有很多優勢，如對人體傷害小，容易操作，設備體積小且攜帶方便等。在木材應力波無損檢測技術跟隨時代，向數字化、系統化、直觀化的方向發展中，應力波與CT技術的結合占據了主導地位。為此，筆者開展了基于改進的層析成像圖像重建ART算法的木材應力波檢測研究。

1 "層析成像簡介及數學基礎

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

CT技術通過重建木材內部圖像來對其進行檢測從而判斷或確知有否缺陷，從而提高木材利用率。以射線理論為基礎的波速走時成像是應力波層析成像的特點，是近幾年來發展起來的，一個木材檢測方面比較新的成像方法。

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

對于任意的函數f，記其對應的Radon變換為Rf，其中Rf是f沿一直線l的線積分，將這個積分值亦稱投影值。函數的Radon變換定義為：

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

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

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

（2）

式中，P（l，θ）=Rf（l，θ），Pl（l，θ）表示P（l，θ）關于l的偏導數。

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

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

2.1 分析法

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

2.2 代數重建法

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

2.2.1

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

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

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

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

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

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

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

2.2.2

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

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

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

該重建算法用于圖像重建是存在兩個缺點：①圖像在投影計算時會產生重建場的嚴重噪聲;②需要很大迭代次數，重建效率不高。

2.3 聯合迭代算法

聯合迭代算法（Simultaneous Iterative Reconstruction Techniques，簡稱SIRT）與ART相似，不同的地方在于ART是每次的修正都只考慮一條射線，而SIRT方法中，一個方格的平均修正值的確定是將方格內通過的所有的射線修正值都考慮在內進行計算的。SIRT的這種迭代算法的誤差修正計算方法的改進，可以有效地抑制這些缺點，其迭代公式如下[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個像素第n次迭代后的慢度值;Mj為通過第j個像素的射線條數;ΔS（n）j為第j個像素在第n次迭代中求得的慢度修正量;lij為第i根射線通過第j個像素的射線長度;lik為第i根射線在第k個像素中的長度;Ti為第i條射線的觀測走時;m為離散的像素總數;S（n）k為第k像素在n次迭代后的慢度值。通過以上公式，根據初始的慢度模型及觀測走時即可求出所有像素慢度值。

2.4 阻尼最小二乘法

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

AX=B

min‖AX=B‖2 （8）

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

2.5 改進的ART算法

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

2.5.1

確定缺陷單元網格。

vi=∑jaij/ti （9）

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

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

（11）

式中，vi為第i條射線波速;v為波速均值;sv為波速標準差;n為射線總數。

假設同一測區射線波速服從正態分布，則可得到臨界波速公式：

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

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

2.5.2

確定迭代初值。把沒有穿過缺陷的射線長度除以走時值即作為正常射線平均波速，把最小的正常射線波速賦值給預判的缺陷單元。

2.5.3

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

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

（13）

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

3 "ART算法及其改進算法的實現與比較

3.1 試驗模型的構造及數據采集

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

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

3.2 ART和SIRT成像算法實現及結果

對得到的慢度向量，使用軟件SURFER 8.0進行直接畫圖。實物圖像見圖1，ART及改進的ART算法的圖像如圖2、3所示。

圖1 實物圖像

圖2 ART算法圖像

圖3 改進的ART算法圖像

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

參考文獻

[1]

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

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

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

[4] 王學勝.超聲層析技術中射線追蹤方法的研究與應用[D].北京：中國地質大學（北京），2005.

[5] 沈寬.工業CT窄角扇束的迭代圖像重建算法研究[D].重慶：重慶大學，2002.

[6] 方璽.混凝土結構的超聲探傷信號處理及二維CT軟件設計[D].武漢：武漢理工大學，2006.