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

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

0 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(4)比較數(shù)組Value的第2列中前2個(gè)數(shù)值的大小，確定像素?cái)?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ù)值較小點(diǎn)的數(shù)據(jù)值和位置作為局部極小值點(diǎn)的數(shù)據(jù)值vally_y和位置vally_x，把數(shù)值較大點(diǎn)的數(shù)據(jù)值和位置作為局部極大值點(diǎn)的數(shù)據(jù)值peak_y和位置peak_x；

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

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

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

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

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

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

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

2 材料與方法

2.1 試驗(yàn)材料及處理方法

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

2.2 閾值確定方法

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

2.3 年輪識(shí)別誤差計(jì)算方法

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

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

e1=m2/m,

(1)

e2=m3/m，

(2)

e=e1+e2,

(3)

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

(4)

2.4 年輪寬度測量精度計(jì)算方法

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

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

(5)

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

3 結(jié)果與分析

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

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

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

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

續(xù)表1

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

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

4 結(jié)語

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