基于峰谷分析算法用針刺儀測定樹木年齡的可行性分析

潘 虹，盧 軍，郭旭展，唐守正，高瑞東，徐建軍

(1. 中國林業科學研究院資源信息研究所，北京 100091；2. 信陽師范學院數學與統計學院，河南 信陽 464000；3. 管涔山國有林管理局，山西 寧武 036700；4. 管涔山國有林管理局羊圈溝林場，山西 五寨 036200)

活立木年齡微損測量是解決很多林業問題的關鍵，對森林經營、生長模擬以及古樹保護等都具有重要的意義[1-4]。傳統的確定活立木年齡是通過識別樹木的年輪數量[5-6]，這種方法雖然比較準確，但是需要用生長錐取樹木生長芯，或者伐倒樹木截取圓盤，對樹木具有侵入性和破壞性；無損的確定活立木年齡的方法有查數輪生枝法和建立非線性數學模型[7-11]等，由于樹木有時一年形成二層輪枝，查數輪生枝法確定樹木年齡精度不高；通過建立數學模型來估算樹木的年齡方法不具有普遍性，并且準確性較低[12]，實際生產中并不能得到廣泛的應用。

樹木針刺儀是微損的測定樹木年齡的工具，將其應用于活立木，可以獲得相對阻力折線圖[13]，折線圖中峰谷的趨勢反映樹木年輪內早材和晚材邊界引起的變化[14]，可以估計樹木的年齡和樹木生長率[15-16]。樹木針刺儀自帶DECOM 軟件可以對抗鉆阻力折線圖進行分析，通過人工干預半自動檢測年輪邊界，估計樹木年齡[17]，由于針刺儀的鉆針軸直徑只有1.5 mm，對樹干的損傷很小，使用針刺儀確定樹木年齡的方法對活立木傷害小，是微損的設備，但是DECOM 軟件自動識別結果很不準確，很多真實年輪邊界不能正確識別，誤差很大，無法在林業生產中投入使用。

本研究以德國RINNTECH 公司生產的樹木針刺儀（Resistograph 4450P/S）鉆入活立木獲取的抗鉆阻力值序列折線圖為研究對象，給出峰谷分析算法，根據確定閾值Det，記錄關鍵點值和關鍵點序號，來分析利用峰谷分析算法估計活立木年齡的可行性，為針刺儀自帶DECOM 軟件改進提供理論支撐，推進微損測定樹木年齡方法的研究進程。

1 數據

從吉林汪清金溝嶺林場天然林中選取52 株活立木作為研究對象，試驗樹木選擇的都是樹干通直，沒有病蟲害，健康活立木。試驗樹木包括紅松(Pinus koraiensis Sieb. et Zucc.)20 株 、 冷 杉（Abies fabri（Mast.）Craib）32 株，紅松的年齡分布范圍是42 a 至88 a，冷杉的年齡分布范圍是35 a 至92 a。首先使用生長錐取出通過髓心的樹芯，共獲取52 個有效樹木生長芯作為參考樣本。然后將針刺儀垂直于樹干中心，分別在取樹芯鉆孔位置上下方5 cm 以內進行針刺，確保針頭盡可能鉆過髓心，抗鉆阻力值自動保存記錄，每個樹芯對應一組抗鉆阻力值，獲取紅松的抗鉆阻力數據為40 組，冷杉的抗鉆阻力數據64 組，共104 組有效抗鉆阻力數據作為研究樣本。對參考樣本進行拋光、打磨和掃描，借助WinDENDRO 年輪分析系統軟件確定樹木的年輪數。

2 方法

2.1 針刺儀測量原理

德國RINNTECH 公司開發的樹木針刺儀（Resistograph 4452P/S）是通過電腦控制電子傳感器的鉆刺針，測量樹木的鉆入阻抗，獲取抗鉆阻力折線圖，橫坐標軸刻度代表針頭鉆入的深度（單位： mm）， 縱坐標軸刻度為阻力值（ 單位：resi）。根據阻力和密度之間的線性關系，可以利用折線圖中波峰波谷的趨勢來確定活立木早材和晚材邊界[18]。使用針刺儀自身攜帶的DECOM 軟件可以自動分析折線圖，識別年輪邊界，估計樹木的年齡。

2.2 峰谷分析算法

由于樹木春材和秋材的密度不同，相應抗鉆阻力值總會有差異，在抗鉆阻力折線圖中可以很直觀的看到很多波峰和波谷出現，其中有很多躍度較小的峰谷，并不代表年輪變化。峰谷分析算法可以根據實測樹芯年輪數，確定躍度的閾值Det，去掉偽年輪，留下真年輪，并將有效的年輪峰值和谷值記錄下來，確定最終有效的峰谷數來估計樹木的年齡。

2.2.1 確定閾值和極值點 首先對獲取的抗鉆阻力剖面進行去首尾處理，去掉兩端無效數據，保留木質部的數據，處理后的針刺儀獲取樹木的相對阻力記錄為整數序列 x ={x1, x2,···, xn}， n為針刺儀測量點的序列號。峰谷分析算法是從抗鉆阻力序列 x起始點 k0開始，設定閾值Det，允許誤差 e = 0.2，確定峰值點和峰谷點序號，以及相應的極值。各項初始值如下：

極值點進入點：

K max = k0,K min = k0,

極大值點起始進入點和退出點分別記為：

K max sta = K max,K maxend = K max,

極小值點起始進入點和退出點分別記為：

K min sta = K min,K minend = K min,

極值點 X max = X(k0),X min = X(k0)，極值點起始點 K sta = 0。

分成3 種情況進行分析：從起始點 X (k0)以后未知升降段、升段、降段。

輸入為： x，Det，e， k0，升降段表示islor2，取值為0,1,2，其中islor2=1 表示升段，islor2=2 表示降段，islor2=0 表示未知升降；

輸出為：極值點起點序號 Ksta，極值點最后一點序 號 Kend ， 極 值 點值 Vmid； 升降 段 表示is12 = YZ 兩位數， Y 是本段代碼， Z是下段代碼，is12 =-1表示出錯；其中，

Y = 1表 示本段升， Y = 2表 示本段降， Y = 3表示本段是頭（不完整段）；

Z = 1表 示下段升， Z = 2表 示下段降， Z = 3表示下段是尾（不完整段）；

流程圖見圖1:

2.2.2 記錄關鍵點序號和值，統計峰谷點個數 根據2.2.1 節得到的極值點起點和退出序號，以及升降段表示，記錄峰點和谷點的數量，以及關鍵點值和關鍵點號。參數說明如下：

根據2.2.1 節輸出極值點起點序號 Ksta，極值點最后一點序號 Kend ，極值點值 Vmid；升降段表示 i s12 = YZ 兩位數，將參數 i s12 = YZ拆分成兩個參數 id1和 id2，

圖1 峰谷分析流程圖Fig. 1 Peak-valley analysis flow chart

則有：

id1 = 1表 示本段升， id1 = 2表 示本段降，id1 =3表示本段是頭（不完整段）；

id2 = 1表 示下段升， id2 = 2表 示下段降，id2 =3表示下段是尾（不完整段）；

用 id1 = 1表 示波峰， id1 = 2表 示波谷，id1 =3表示端點；

矩陣P 為 m行3 列矩陣，第一列記錄關鍵點號，第二列記錄關鍵點值，第三列記錄 id；

矩陣 N為1 行2 列的矩陣， N =(No1,No2)，其中 N o1統 計峰點的個數， N o2統計谷點個數。

2.3 步驟

輸入：針刺儀測量數據去首尾得到的列向量X ={x1, x2,···, xn}和閾值Det；

step1. 將原來針刺儀數據去首尾，得到的數據

step2. 分析初始段

此 時 取 參數islor2=0， k0= 1， e = 0.2， 根 據

2.2.1 得到相應的輸出Ksta,Kend,Vmid,is12。

若 id2 = 3， 記錄矩陣P 為 m行3 列矩陣，第一列記錄關鍵點號，第二列記錄關鍵點值，第三列記錄 id。X ={x1, x2,···, xn}， 加上一列序號 X1={1,2,···,n}, 分別以 X1作 為第一列， X作為第二列，構造矩陣的DX，則矩陣DX 表示帶序號的觀測值，可以表示為：

若i d2 = 1，分成兩種情況討論：

若i d1 = 3，此時有不完整前段，記錄

若 i d2 = 2，此時后段是完整段，分成兩種情況討論：

step3. 初始段之后，各段輪流升降

當 Kend ＜n時 ，取參數islor2=id2， k0= Kend，e = 0.2，根據2.1.1 節得到相應的輸出Ksta，Kend，Vmid，is12。記錄 P為分塊矩陣，將所有的關鍵點序號和關鍵點值依次循環寫入，即

輸出：矩陣 P ， 矩陣 N o =[No1,No2]。

峰谷分析算法的估計年齡 A可以由峰谷數的均值的二分之一取整數得到，即 A =(No1+ No2)/4取整數。

2.4 年輪寬度計算

根據2.3 節輸出的結果矩陣 P ， P 中的第一列是關鍵點序號，即峰點起點序號和結束點序號以及谷點起點序號和結束點序號，第二列是關鍵點值，第三列表示關鍵點序號是峰點或者谷點。由于針刺儀的取樣規律是每1 mm 取100 個點，說明抗鉆阻力兩個相鄰的序號之間的距離就為0.01 mm，假設矩陣 P每個波峰的第一個波峰點所在的行號為

提取 P 所有的 li,i = 1,2,···,No1行，重新組成新的矩陣，記為

將矩陣 Q中的第一列相鄰兩個元素做差分，得到相鄰兩個峰點之間序號之差，因針刺儀抗鉆阻力數據采取間隔是每1 mm 取100 個點，則相鄰序號差值的百分之一可以是對應年度樹木的估計年輪寬度，年輪寬度單位為mm。

3 結果

根據樹芯的實測年齡， 選擇適當的閾值Det 后，將峰谷分析算法應用于104 組活立木抗鉆阻力序列，能夠較好地去除樹木偽年輪。圖2 為研究樣本紅松3 的峰谷分析結果圖，當選擇閾值為54 時，得到的峰谷數均值為45，與實際年齡吻合，將峰谷關鍵點值和關鍵點序號連線標記可以直觀看出，峰谷值躍度低于54 的峰谷點都被忽略，只留下峰谷值躍度大于等于54 的峰谷點。可以看出，確定恰當的閾值后，峰谷分析算法估計的年齡與實際年齡非常接近。峰谷分析算法估計年齡平均絕對誤差是－2，范圍在－5 年至5 年之間，平均相對誤差為－2.69%，范圍在－6.73% 至6.73% 之間；每組數據的相對誤差分布見圖3。

對真實年齡與算法估計年齡進行成對數據t-檢驗，數據檢驗得到t 值為1.31，由于 |t2| ＜2.2622，故不能拒絕原假設，即可以為峰谷分析算法估計樹木的年齡均值與樹木的真實年齡均值無顯著差異，此時檢驗的 p值為0.192。

4 討論

4.1 峰谷分析方法的優點

本算法數據是基于樹木針刺儀獲取抗鉆阻力值，由于針刺儀的鉆針直徑小，對樹木損害小，與傳統的樹木年代學方法估計樹木年齡[19-20]相比較，本算法估計樹木年齡是屬于微損的，或者準無損的方法，適用于受保護的天然林管理和古樹年齡鑒定。由于針刺儀體積小，質量輕，屬于便攜式設備，克服了以往儀器測定樹木年齡法的笨重、操作繁瑣等缺點，便于野外試驗進行操作。

4.2 提高針刺儀測定樹木年齡精度的可行性

針刺儀自身攜帶有年輪分析軟件DECOM，能對獲取的抗鉆阻力數據自動劃分年輪邊界，但是通過與參考樣本獲取的實際年輪數對比發現，自動檢測的年輪邊界準確率較低，平均相對誤差為－92.27%。使用DECOM 軟件檢測年輪邊界時，需要靠人工經驗進行視覺判讀，手動添加年輪邊界，耗時，效率低。而本算法通過MATLAB 編程實現，將抗鉆阻力值作為峰谷分析算法的輸入，根據實測樹芯年齡來選擇適當的參數Det 值，可以得出樹木的估計年齡，并且與實際年齡誤差較小，平均相對誤差降低至－2.69%。樹木針刺儀通過峰谷分析算法測定樹木的年齡是可行的。

圖2 紅松3 峰谷分析結果Fig. 2 Peak-valley analysis results of Korean pine 3

圖3 峰谷分析算法估計年齡與實測年齡相對誤差分布Fig. 3 Relative error distribution between Peak-valley analysis algorithmand true age

4.3 閾值Det 選取

峰谷分析過程中，設置了很多過程參數來輔助描述分析細節，分析過程細致周密，并且使用起來非常方便。峰谷分析算法輸入端只有兩個變量：抗鉆阻力值序列X和閾值Det，隨著閾值Det的取值不同會得到不同的峰谷數。確定閾值Det后，當峰或者谷的躍度，或者高度，超過Det之后作為一個完整段。在本研究中，紅松的閾值取值范圍為13～67，冷杉的閾值取值范圍為6～55。由于Det的分布范圍較大，可以考慮用機器學習的方法，找出每棵樹的特征，來研究去掉樹木的假年輪對應的恰當閾值Det值。

4.4 年輪寬度的計算

本算法輸出結果中，將每棵樹所有的峰點和谷點的序號，以及對應的抗鉆阻力值記錄到矩陣中。記錄的過程采用分塊矩陣的形式，為計算樹木的年輪寬度提供了依據。雖然年輪寬度估計值和實測之間有差別，但是估算的樹木年輪寬度能基本上反映樹木的年生長趨勢，可以作為估算樹木連年生長量的一個參考。

5 結論

確定恰當的閾值后，通過給出的峰谷分析算法，可以將針刺儀獲取的抗鉆阻力數據作為輸入，得出樹木年齡的估計值，峰谷分析算法應用于針刺儀抗鉆阻力序列來估計樹木年齡是可行的，確定存在恰當的閾值使針刺儀估計樹木年齡精度很高，閾值的選擇依據是下一步研究的重點內容。