基于灰度圖像的峰谷年輪識別方法研究

姚建峰，郭旭展，宋新宇，王靜賢，祁 柯，代琪琦

(信陽師范學(xué)院 a.計算機與信息技術(shù)學(xué)院； b.數(shù)學(xué)與統(tǒng)計學(xué)院， 河南 信陽 464000)

0 引言

在溫帶和亞熱帶地區(qū)，樹木在生長季早期形成的木材顏色較淺，通常把這部分木材稱為“早材”；在生長季晚期形成的木材顏色較深， 通常把這部分木材稱為“晚材”[1]。當(dāng)年的晚材與次年的早材之間界限分明，呈現(xiàn)紋輪，稱該紋輪為年輪線。相鄰年輪線之間的距離稱為年輪寬度。樹干基部的年輪個數(shù)記錄了樹木年齡，年輪寬度記錄了氣候、環(huán)境和森林經(jīng)營等綜合外界因子對樹木生長的影響[2-3]。通過研究年輪寬度與氣候因子之間的模型，可重建過去的氣候變化情況，甚至可預(yù)期未來的氣候變化規(guī)律[4]。同時，樹木年齡是森林資源調(diào)查的一項基本指標(biāo)，是林業(yè)經(jīng)營管理中重要的時間指標(biāo)，在確定林分齡級和齡組、調(diào)整林分齡級配置、評價森林碳匯潛力、制定森林經(jīng)營方案等方面具有重要意義[5-8]。

目前測量樹木年輪的方法主要有顯微鏡法、X射線法和圖像法[9]。

顯微鏡法是使用顯微鏡或者放大鏡，人工測量經(jīng)過打磨的圓盤或者生長錐木芯獲取樹木年輪信息。該方法是基于實物的測量年輪方法，對樣品的存儲環(huán)境要求較高、檢查或者復(fù)測樣品困難。

X射線法是使用X射線掃描經(jīng)過打磨的圓盤或者生長錐木芯的圖像，然后使用X射線圖像獲取樹木年輪信息[10]。該方法可獲取準(zhǔn)確的年輪寬度和年輪密度信息，但是該方法設(shè)備昂貴，使用成本和維護成本較高，操作復(fù)雜[11]。

圖像法是使用掃描儀獲取經(jīng)過打磨的圓盤或者生長錐木芯圖像，然后使用年輪圖像處理軟件測量樹木年輪信息[12]。年輪圖像處理軟件有商業(yè)化的軟件和科研人員自主開發(fā)的年輪識別軟件。目前最具有代表性的商業(yè)化年輪圖像處理軟件是加拿大Regent公司開發(fā)的WinDendro軟件[13]，該軟件可進(jìn)行年輪標(biāo)記、獲取年輪個數(shù)和年輪寬度等數(shù)據(jù)，但是，該軟件價格較高，國內(nèi)僅少數(shù)科研機構(gòu)可以承擔(dān)。

隨著計算機圖像技術(shù)[14-16]的發(fā)展，越來越多的科研人員使用圖像處理技術(shù)獲取樹木年輪信息[17-19]。圖像處理技術(shù)涉及圖像濾波算法、早晚材邊界分割算法和早晚材邊界增強算法等，不同圓盤圖像需要設(shè)置不同的參數(shù)才能達(dá)到較高的年輪識別精度。因此，基于圖像處理技術(shù)的年輪識別方法要求使用人員具有較好的數(shù)字圖像處理技術(shù)基礎(chǔ)，熟悉相關(guān)算法的原理，掌握相關(guān)參數(shù)調(diào)試方法。由于大部分基層林業(yè)工作人員不熟悉數(shù)字圖像處理技術(shù)，所以，該方法不適合廣大基層林業(yè)工作人員使用。潘虹等[20]提出了一種基于微鉆阻力的峰谷年輪識別算法，該算法只需設(shè)置1個參數(shù)(峰谷差值)的閾值,便可識別出鉆針阻力信號中的樹木年輪信號，使用簡單。本文使用峰谷年輪識別算法識別圓盤圖像中的年輪信息，并計算年輪識別的精度，為自動識別圓盤圖像年輪信息提供一種簡單、可行的方法。

1 灰度圖像峰谷年輪識別原理及實現(xiàn)過程

年輪圖像是由多個明暗相間的同心圓環(huán)構(gòu)成，圓環(huán)中心是樹木髓心。早材顏色較淺，圖像灰度值較大；晚材顏色較深，圖像灰度值較小。因此，圓盤灰度圖像中經(jīng)過樹木髓心的每行灰度數(shù)據(jù)的曲線呈峰谷交替規(guī)律變化。在基于微鉆阻力的峰谷年輪識別算法中，鉆針阻力波峰對應(yīng)年輪晚材，鉆針阻力波谷對應(yīng)年輪早材，與灰度曲線中的對應(yīng)關(guān)系相反。

為了便于使用峰谷年輪識別算法處理年輪灰度曲線，本文取灰度值的補碼進(jìn)行年輪識別處理。由于圓盤圖像中存在樹脂斑點、打磨擦痕等噪聲信號，如果只取1行灰度數(shù)據(jù)值進(jìn)行年輪識別處理，這些噪聲信號對年輪識別結(jié)果影響就較大，年輪識別精度不高。

為了減小噪聲信號對年輪識別的影響，本文選取經(jīng)過髓心、寬為3～5 mm的圓盤圖像，然后取每列灰度補碼的平均值作為該列在行方向上各像素點的灰度數(shù)據(jù)，從而把行、列兩個方向上的灰度數(shù)據(jù)轉(zhuǎn)化為行方向上的灰度數(shù)據(jù)。

為了減小行方向上噪聲信號對年輪識別的影響，首先在年輪識別前對行方向上的灰度數(shù)據(jù)進(jìn)行平滑處理，然后再對平滑后的灰度數(shù)據(jù)進(jìn)行年輪識別處理。

峰谷年輪識別算法是根據(jù)波形圖中波峰與波谷之間的差值來識別波形圖中的波峰是否為年輪信號。該算法首先設(shè)置波峰與波谷之間差值的閾值Det。當(dāng)波峰與相鄰波谷之間的差值大于或者等于設(shè)定的閾值時，則把該波峰記為樹木年輪信號；當(dāng)波峰與相鄰波谷之間的差值小于閾值時，則把該波峰記為噪聲信號。將年輪信號的個數(shù)作為選取圓盤圖像部分的年輪數(shù)，將年輪信號波峰極大值點作為年輪之間的分界點，將相鄰年輪信號波峰極大值點之間的距離作為年輪寬度。

峰谷年輪識別函數(shù)的處理過程如下：

(1)獲取圖像灰度數(shù)據(jù)：根據(jù)圖像名Name值，讀取圓盤圖像的數(shù)據(jù)，并去掉背景色的數(shù)據(jù)，然后將彩色圖像數(shù)據(jù)轉(zhuǎn)換成灰度圖像數(shù)據(jù)。

(2)計算每列像素點的灰度平均值，并對灰度平均值取補碼，得到行方向像素點的一維灰度補碼數(shù)據(jù)，使用平滑算法對一維灰度補碼數(shù)據(jù)進(jìn)行平滑處理，平滑窗口長度為Det，并記錄行方向上的像素點的個數(shù)m。

(3)根據(jù)圖像分辨率，計算行方向上各像素點的位置，并存儲在二維數(shù)組Value的第1列中；每個像素點的灰度補碼值存儲在數(shù)組Value的第2列中。

(4)比較數(shù)組Value的第2列中前2個數(shù)值的大小，確定像素數(shù)據(jù)變化趨勢trend是上升段還是下降段。若Value [1][2]< Value[2][2]，則數(shù)據(jù)變化趨勢trend處于上升階段，把trend賦值為1;若Value [1][2] ≤Value [2][2]，則阻力變化趨勢trend處于下降階段或者穩(wěn)定階段，把trend賦值為0。把數(shù)值較小點的數(shù)據(jù)值和位置作為局部極小值點的數(shù)據(jù)值vally_y和位置vally_x，把數(shù)值較大點的數(shù)據(jù)值和位置作為局部極大值點的數(shù)據(jù)值peak_y和位置peak_x；

(5)當(dāng)數(shù)據(jù)變化趨勢處于上升階段時，若新數(shù)據(jù)點的值大于局部極大值，則把新數(shù)據(jù)點作為局部極大值數(shù)據(jù)點；若新數(shù)據(jù)點的值比局部極大值點小，且它們之間的差值d小于設(shè)置的閾值Det，則認(rèn)為該數(shù)據(jù)點是一個噪聲信號；若新數(shù)據(jù)點的值比局部極大值點小，且它們之間的差值d大于或者等于設(shè)置的閾值Det，則把該數(shù)據(jù)點作為本年輪結(jié)束的標(biāo)志，數(shù)據(jù)變化趨勢從上升階段變?yōu)橄陆惦A段或者穩(wěn)定階段，把局部極大值點作為一個年輪信號進(jìn)行處理，有效波峰的個數(shù)n加1，把新數(shù)據(jù)點的值和位置分別作為局部極小值點的值vally_y和位置vally_x；

(6) 當(dāng)數(shù)據(jù)變化趨勢處于下降階段或者穩(wěn)定階段時，若新數(shù)據(jù)點的值比局部極小值點的阻力值小，則把新數(shù)據(jù)點作為局部極小值數(shù)據(jù)點；若新數(shù)據(jù)點的值比局部極小值大，且它們之間的差值d小于設(shè)置的閾值Det，則認(rèn)為該數(shù)據(jù)點是一個噪聲信號；若新數(shù)據(jù)點的值比局部極小值大，且它們之間的差值d大于或者等于設(shè)置的閾值Det，則把該數(shù)據(jù)點作為該年輪晚材部分的開始階段，數(shù)據(jù)變化趨勢從下降階段或者平穩(wěn)段變?yōu)樯仙A段，把局部極小值點作為一個年輪信號的谷值點進(jìn)行處理，把新數(shù)據(jù)點的值和位置分別作為局部極大點的值peak_y和位置peak_x。

灰度圖像峰谷法識別年輪的流程圖如圖1所示。

圖1 峰谷算法流程圖Fig. 1 Flow chart of peak-valley algorithm

使用Matlab語言編寫峰谷年輪識別函數(shù)實現(xiàn)峰谷年輪識別算法。識別函數(shù)的定義如下：

[Value, Value_Key, Wide, n]=Peak_Vally(Name, Resolution, Det)。

在形式參數(shù)中，參數(shù)Name是一個字符串，用來傳遞灰度圖像的文件路徑和文件名；參數(shù)Resolution用來傳遞圖像的分辨率，單位為dpi；參數(shù)Det用來傳遞波峰與波谷之間差值的閾值。函數(shù)返回值中，參數(shù)Value用于存儲轉(zhuǎn)化后的灰度圖補碼數(shù)據(jù)，是一個二維數(shù)組，第1列用于存儲像素點的位置，單位為mm，第2列用于存儲像素點像素值的補碼；參數(shù)Value_Key用于存儲峰谷識別后的關(guān)鍵點數(shù)據(jù)，即年輪信號波峰極大值點和波谷極小值點的位置和像素數(shù)據(jù)，是一個二維數(shù)據(jù)，第1列存儲像素點的位置，第2列像素值的補碼，第3列存儲波峰或者波谷的標(biāo)志，若該像素點是波峰，則該數(shù)組元素為1；若該像素點是波谷，則該數(shù)組元素為0；參數(shù)Wide是一維數(shù)組，用于存儲年輪寬度；參數(shù)n是樹木年輪信號的個數(shù)，即數(shù)組Value_Key中波峰的個數(shù)。

2 材料與方法

2.1 試驗材料及處理方法

試驗材料是2020年9月在信陽師范學(xué)院校內(nèi)的天然次生林中采樣的12棵近2個月內(nèi)枯死的馬尾松，在樹干高度為1.3 m處附近截取樹干通直、無明顯缺陷、厚約10 cm的圓盤，并記錄圓盤北向方向。使用打磨機打磨圓盤，使年輪線清晰可見。在圓盤光面東西、南北方向兩個方向上選取寬為3～5 mm木材，兩個方向的木材都要經(jīng)過圓盤髓心，并用記號筆做好標(biāo)記。使用Lintab6年輪分析儀測量每個圓盤東西、南北2個方向上的年輪線個數(shù)和年輪寬度。把掃描儀分辨率設(shè)置為400 dpi，使用掃描儀掃描圓盤獲取圓盤彩色圖像，使用畫圖軟件選取圓盤圖像上標(biāo)記區(qū)域的木材部分，并做好記錄。

2.2 閾值確定方法

閾值設(shè)置值要合理，若閾值設(shè)置值過大，則把波峰與波谷差值較小的年輪信號判別為噪聲信號；若閾值設(shè)置值過小，則把波峰與波谷差值較大的噪聲信號判別為年輪信號。本文根據(jù)3個圓盤圖像在不同閾值下識別結(jié)果與圓盤圖像的對比分析圖，選擇能正確識別大部分年輪的閾值作為所有圓盤圖像的閾值。

2.3 年輪識別誤差計算方法

對照圓盤圖像和像素數(shù)據(jù)年輪識別圖，若像素數(shù)據(jù)年輪識別圖中的年輪線與圓盤圖中的晚材位置相對應(yīng)，則把該年輪線標(biāo)記為正確的年輪線；若像素數(shù)據(jù)年輪識別圖中的年輪線與圓盤圖中的早材位置相對應(yīng)，則把該年輪線標(biāo)記為錯誤的年輪線；若圓盤圖的晚材位置沒有對應(yīng)的年輪線，則把該年輪標(biāo)記為未能識別的年輪。

保存圓盤圖和像素數(shù)據(jù)年輪識別圖的匹配圖，分別記錄每個像素數(shù)據(jù)的預(yù)估的年輪個數(shù)m1、誤判的年輪個數(shù)m2和漏判的年輪個數(shù)m3。每個取樣方向上以Lintab6測量的晚材個數(shù)為m為年輪數(shù)的真值，每個圖像數(shù)據(jù)的年輪誤判率e1、年輪漏判率e2、年輪識別錯誤率e和年輪個數(shù)測量精度ε1,計算公式分別如下:

e1=m2/m,

(1)

e2=m3/m，

(2)

e=e1+e2,

(3)

ε1=1-|m1-m|/m。

(4)

2.4 年輪寬度測量精度計算方法

對照圓盤圖和像素數(shù)據(jù)年輪識別圖的匹配圖，首先分析圓盤圖像測量的年輪寬度數(shù)據(jù)，把錯誤年輪線所對應(yīng)的年輪寬度進(jìn)行合并，計算相鄰的2個正確識別年輪線之間的寬度w1。然后分析Lintab6測量的年輪寬度，若2個正確識別年輪線之間沒有漏判的年輪，則把Lintab6測量該年輪的寬度作為該年輪寬度的真值w2；若2個正確識別年輪線之間存在漏測的年輪，則把2個正確識別年輪線之間所有年輪的寬度和作為這2個正確識別年輪線之間的年輪寬度真值w2，最后使用式(5)計算每相鄰2個正確識別年輪線測量年輪寬度的精度ε0:

ε0=1-|w1-w2|/w2。

(5)

把每個圓盤每2個正確識別年輪線測量年輪寬度的精度ε0的平均值，作為每個圓盤圖像測量年輪寬度的平均精度ε2。

3 結(jié)果與分析

根據(jù)3個圓盤圖像在不同閾值下的識別結(jié)果和圓盤圖像的對比分析圖可知，當(dāng)閾值設(shè)置為6時，能正確識別圓盤中大部分年輪信息。圖2是其中一個圓盤圖像的年輪識別結(jié)果與圓盤彩色圖像和灰度圖像的對比圖。從圖2中可以看出，盡管原始圖像中存在色斑、霉變、不清晰的年輪線等問題，但是峰谷年輪算法仍可正確識別圓盤中大部分年輪信息。

圖2 年輪識別結(jié)果與圓盤圖像對比Fig. 2 Comparision between tree-ring recognition results and disc image

通過峰谷年輪識別結(jié)果與圓盤圖像對比，在圓盤圖像晚材位置標(biāo)記漏判的年輪，在年輪識別圖中的年輪分界線上標(biāo)記誤判的年輪，然后計算年輪識別的漏判率、誤判率、錯誤率和年輪寬度測量精度，各圖像數(shù)據(jù)的年輪識別結(jié)果如表1所示。

表1 年輪識別結(jié)果

續(xù)表1

從表1中可以看出，峰谷年輪識別算法的平均漏判率為10.084%，平均誤判率為4.059%，平均錯誤率為14.143%，年輪個數(shù)的平均測量精度為90.096%，年輪寬度的平均測量精度為89.263%。從總體識別結(jié)果來看，峰谷年輪識別算法預(yù)測樹木年齡精度和年輪寬度精度高達(dá)90%左右。

部分圖像數(shù)據(jù)的誤判率或者漏判率較高，可能有以下原因：(1)設(shè)定的閾值對部分圓盤圖像不合適：本文對所有數(shù)據(jù)都使用相同的閾值進(jìn)行識別，對有些圓盤圖像閾值偏大，導(dǎo)致早材與晚材顏色差異較小的年輪信號被識別為噪聲信號，使漏判率偏高；對有些圓盤圖像閾值偏小，導(dǎo)致把部分峰谷差值較大的噪聲信號識別為年輪信號，使誤判率偏高。在實際使用過程中，可通過人工方式標(biāo)記不清晰的年輪線、對不同圓盤圖像設(shè)置不同的閾值等方式來進(jìn)一步提高年輪識別精度。

4 結(jié)語

首先使用列方向求均值法，將二維圓盤灰度數(shù)據(jù)轉(zhuǎn)化為一維灰度數(shù)據(jù)，然后對一維灰度數(shù)據(jù)使用峰谷年輪識別算法獲取年輪信息。經(jīng)測試，峰谷年輪識別算法預(yù)測樹木年齡精度和年輪寬度精度高達(dá)90%左右。該算法易于理解，使用簡單，可為林業(yè)工作人員提供一種簡單、可行的年輪識別方法。