基于壓縮感知和空間插值的樹木內部缺陷應力波層析成像算法

李佳捷，杜曉晨

（浙江農林大學信息工程學院，杭州311300）

0 引言

樹木內部缺陷成像技術具有重要的研究價值，使用該技術，能精準、無損地對樹木內部健康情況進行檢測，從而有效保護各種樹木、林木、原木、木材和木質建筑[1]。在眾多檢測信號中，應力波信號對人體無害且成本低廉，特別適合用于樹木內部缺陷檢測與成像，已成為本研究領域的主流信號[2]。該成像技術的基本原理是：應力波經過木材缺陷區域的波速比健康區域慢[3]。根據該原理，將若干傳感器均勻布置于某樹木斷層，當敲擊其中一個傳感器產生應力波時，其余傳感器接收該信號并計算出傳播時間。由于通過不同材質時波的衰減程度不同，導致所接收的傳播速度存在差異，根據所接收傳播速度矩陣進行分析、計算和圖像重構，最終以二維圖像方式直觀地顯示樹木斷層內的缺陷情況。

顯然，傳感器的數量多少是影響應力波成像精度的重要因素[4]。在傳統的樹木內部缺陷應力波層析成像方法中，均采用固定數量（一般為12 個）傳感器組去獲得足夠進行圖像重建的信號。但是傳感器的布置過程相對耗時，因此，應減少傳感器數量用于信號采集，進行快速層析成像以提高該技術的適用性，同時保證較高的成像精度。本文提出一種基于壓縮感知與空間插值的應力波聯合層析成像算法，將壓縮感知理論與應力波成像相結合，發揮壓縮感知對應力波信號稀疏表示和求解的優勢，結合空間插值算法恢復非缺陷空間區域的優勢，實現樹木內部缺陷的快速應力波層析成像。

1 壓縮感知原理

壓縮感知理論是建立在信號稀疏表示和逼近理論基礎上的方法[5]，近年來在醫學成像、地震成像、雷達成像等眾多領域得到了廣泛應用。該原理可以表示為：

其中，x 代表原始信號，需滿足稀疏性（或在某個變換域內是稀疏的）；y 代表壓縮后的信號；Φ 代表測量矩陣。壓縮感知方法就是通過少量的觀測次數，在已知測量值y 的情況下，利用凸優化方法求解欠定方程組，高概率地恢復原始信號x。一般情況下，自然信號x 本身并不稀疏，所以需要在某種稀疏基上進行稀疏表示：

其中，Ψ 為稀疏基矩陣，a 是稀疏系數，只有K 個非0 值。通過公式（1）和（2），能得到：

其中，Φ'稱為傳感矩陣，求解a 的逼近值a'，則能得到原始信號的逼近值x'?；趬嚎s感知的信號采集過程就是一個線性投影過程，一個向量y（M×1）可以恢復具有K 個非0 值的a。這樣，通過少量次數的觀測，即可有效地恢復原始信號。

由以上壓縮感知原理可知，該方法的關鍵步驟為：原始信號的有效稀疏表示、測量矩陣的設計和信號恢復算法。其中，測量矩陣通?？梢允且粋€隨機矩陣，而欠定方程組的求解問題可轉換為L1 范數約束最優化問題進行求解[6]。再由公式（2）可知，壓縮感知原理中設計該步驟的目的就是確保原始信號具備稀疏性，從而保證能對欠定方程組進行有效求解。

2 基于壓縮感知和空間插值的應力波層析成像

其中，L 表示應力波射線的長度，T 表示應力波的傳播時間，V 表示應力波的傳播速度，而S 則是應力波波速值的倒數。顯然，公式（4）與壓縮感知中的公式（3）相似，當空白網格點的尺寸足夠小時，S 與a 一樣，具備稀疏性，因此，可以用上式表示一次應力波信號采集過程（敲擊某個傳感器一次后的信號采集過程）。進一步將上式以矩陣形式表示如下：

可以發現，壓縮感知原理與層析成像基本反演方程非常相似，它們之間的差異體現在壓縮感知對原始信號有稀疏表示過程，這可以理解為該方法的求解基礎是強制要求原始信號具備稀疏性。而在樹木內部缺陷應力波層析成像領域，用于圖像重建的空白網格圖本身就具有稀疏性，因此，基于壓縮感知原理，可定義如下波速空間分布求解公式：

其中，矩陣L 中的每一行對應一條應力波傳播射線，lij表示的是第i 條射線，穿過第j 號網格點的距離，N 表示區域中的網格總數，M 表示區域中的射線總數，顯然M 遠遠小于N，這與壓縮感知的原理相符，ti表示某條射線對應路徑上的應力波傳播時間。那么，si就是某網格點中的應力波傳播速度分布值（si在網格圖中按從左到右，從上到下的順序排列）。在確定了各個傳感器在樹木橫截面上的布置坐標后，可以得到矩陣L 中的每個值，而向量T 可以通過設備采集（觀測）得到，這也與壓縮感知的原理相符。

按照壓縮感知的原理，可以求解得到公式（4）和（5）中的向量S，依次敲擊傳感器進行觀測后，即可通過少量觀測數據和少量觀測次數，得到期望的應力波傳播速度平面分布。但是，公式（5）與（3）一樣，仍然是一個病態方程組，根據壓縮感知原理，求解得到的向量S中，僅有少于M 個非0 值。所以基于該方法并不能直接對所有空白網格點進行波速估計，也就不能對樹木內部缺陷進行完整的高質量圖像重建。

本文的研究目標是在傳感器數量減少、采樣信號量減少時，仍能對樹木內部缺陷進行高分辨率的應力波層析成像。雖然以上提出的基于壓縮感知的波速分布求解方法不能直接對缺陷進行高質量成像，但當傳感器數量減少時，壓縮感知方法能充分發揮其在稀疏采樣下的信號表達、求解能力。因此，可以將上述求解得到的波速分布作為應力波層析成像方法中的一個重要參考依據，通過壓縮感知方法幫助恢復缺陷區域的潛在空間結構。另外，在項目組的前期研究中，提出過一種基于加權橢圓空間插值的應力波層析成像算法（Ellipse-Based Spatial Interpolation，EBSI[7]），該算法能對樹木內部缺陷進行有效的層析成像，但成像結果往往過高估計缺陷的面積大小。通過結論也可以認為：基于EBSI 算法的層析成像結果中，正常木材區域的重建可信度較高，該方法更有助于恢復非缺陷區域的空間結構?；谝陨蟽煞矫婵紤]，本文提出了一種基于壓縮感知與空間插值的應力波層析成像算法，算法示意圖如圖1 所示。

圖1 提出的應力波層析成像算法示意圖

本算法的基本思想是充分利用壓縮感知成像與空間插值成像的各自優勢，分別提取兩種算法中可信度較高的重建區域內特征點，在統一的重建圖像空間中進行混合并最終進行聯合層析成像。提出的應力波層析成像算法的具體步驟如下：

步驟1：將所有采集到的應力波傳播速度矩陣進行歸一化，并完成信號修正[7]。

步驟2：使用EBSI 算法進行應力波層析成像，并進行圖像分割，將成像結果分為缺陷區域和非缺陷區域。

步驟3：將修正后的應力波信號按傳感器敲擊順序進行分解，每次敲擊按公式（5）進行應力波稀疏信號表示和求解。

步驟4：重復步驟3，直到每次敲擊都以壓縮感知方法進行處理，再將所有求解結果按成像區域網格圖進行合并。

步驟5：在壓縮感知成像結果中進行特征點提取。

步驟6：在EBSI 成像結果的缺陷區域和非缺陷區域中，分別進行特征點提取。

步驟7：將兩類特征點在統一的成像區域內混合，基于經典的IDW 算法進行空間插值。

步驟8：利用統一的顏色尺度可視化網格圖，重建樹木內部缺陷的斷層圖像。

3 實驗結果

實驗環節中，采用了項目組自主研發的應力波信號測量設備進行應力波信號采集[8]。該設備包括12 個應力波傳感器、數據線、敲擊錘和信號處理箱（如圖2所示），信號處理箱將采集到的應力波傳播時間矩陣傳輸到上位機作為算法的輸入部分。樹木樣本則是一個含有真實空洞缺陷的Camphor 樹種木樁樣本。

圖2 應力波層析成像設備

為了生成高分辨率的樹木內部缺陷重建圖像，每個成像網格點的尺寸設為1 個像素，重建圖像的分辨率則為200×200 個像素點。采用紅色表示最小波速值（缺陷區域）、黃色表示較小速度值、綠色表示最大速度值（健康區域），以顯示理想的傳播射線圖和木樁斷層重建圖像視覺效果。此外，為了檢驗提出方法在信號量不斷減少情況下的層析成像性能，分別使用12 個、10 個、8 個和6 個傳感器進行了應力波層析成像實驗，傳感器則按順時針方向均勻布置于樹木橫截面。實驗環節中，本文還重復了EBSI 方法，并與提出算法進行了對比，以檢驗壓縮感知原理的融入，能否讓提出算法在信號量減少的情況下進行高質量缺陷成像。對比成像實驗結果如圖3 所示。

圖3 傳感器數量遞減時的樹木內部缺陷應力波層析成像結果

以上對比層析成像實驗結果顯示，當傳感器數量相對充足時，兩種算法均能取得較為理想的層析成像結果。但當傳感器數量遞減時，信號量也逐漸減少（傳播射線圖即為可視化后的應力波原始波速矩陣），兩個算法的成像結果均開始變差，特別是EBSI 算法，當傳感器數量減至8 個或6 個時，缺陷區域的成像結果與真實情況差別很大，缺陷圖像甚至出現了扭曲、失真等嚴重錯誤情況。相比而言，本文提出的算法，當傳感器數量減至8 個或6 個時，仍能保證較好的成像效果。該實驗結果表明：在信號量不足的情況下，壓縮感知方法在稀疏表示、求解方面的優勢有助于應力波層析成像算法重建高質量的樹木內部缺陷圖像。

4 結語

本文提出了一種基于壓縮感知和空間插值的應力波層析成像算法，利用壓縮感知恢復缺陷區域的潛在空間結構，當信號量不足時，仍能實現樹木內部缺陷的高質量成像。對比實驗結果表明了本文方法的有效性。如何進一步提高成像算法的精度，特別是在缺陷情況不明顯的情況下，保證較好的成像效果，是項目組下一步的研究目標。