基于頻譜分析用針刺儀測定樹木年齡的算法

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

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

活立木年齡微損測定是目前林業(yè)生產(chǎn)研究中最重要和最基礎的共性難題，對森林管護、森林經(jīng)營以及古樹保護等都具有重要的理論和實踐意義[1-4]。傳統(tǒng)的測定樹木年齡的方法具有破壞性[5-6]，并增加樹木致病感染的風險[7-10]；無損測定樹木年齡的方法如建立數(shù)學模型和查數(shù)輪生枝法等[11-15]，準確性較低。目前，在我國實行天然林保護的環(huán)境中，亟需一種微損的測定樹木年齡的方法。

針刺儀是估計樹木年齡的微損工具，將針刺儀應用于活立木，可以獲得相對阻力剖面[16]。針刺儀的基本工作原理是基于抗鉆阻力與木材密度之間的線性關系，根據(jù)剖面曲線波峰波谷的趨勢能反映年輪內(nèi)部早材和晚材的界限[17]，峰值和谷值分別代表早材和晚材，可以用相對阻力剖面圖來估計樹木年齡以及樹木生長率[18-20]，由于針葉樹和闊葉樹的樹干材質(zhì)不同，導致針刺儀鉆入阻力變化不同，針葉樹種的阻力變化比闊葉樹種的阻力變化更明顯，峰谷的區(qū)分度更好，能更準確地表示樹木生長的變化，所以本研究選取針葉樹中的華北落葉松（Larix principis-rupprechtiiMayr）作為研究對象。雖然使用針刺儀自帶DECOM 軟件對抗鉆阻力剖面圖進行分析，可以自動檢測年輪邊界，估計樹木年齡，但是誤差很大，不能應用于林業(yè)實際生產(chǎn)，基于抗鉆阻力值序列來估算樹木年齡的方法鮮見報道。

數(shù)字信號處理是將復雜的信號用簡單的基本的信號表示，從而將復雜的信號分析轉(zhuǎn)化成對基本信號的分析[21]。已經(jīng)實際應用于廣泛的學科領域，如語音處理、電話信道傳輸、圖像處理和傳輸?shù)萚22-23]。傅里葉變換是數(shù)字信號頻域分析的一種重要方法，可以將時域的信號用虛指數(shù)信號表示，反映了信號的時域與頻域之間的關系[24-25]，是信號分析中不可或缺的重要工具，在各個領域都有廣泛應用[26-27]。

針刺儀鉆入活立木獲取的抗鉆阻力值序列，可以看做離散周期信號，從而可以考慮利用數(shù)字信號處理對其進行分析，來尋找樹木年度變化的信息。本研究以德國RINNTECH 公司生產(chǎn)的樹木針刺儀（Resistograph 4452P/S）鉆入華北落葉松獲取的抗鉆阻力值序列為研究對象，利用數(shù)字信號處理中的傅里葉變換，對抗鉆阻力值序列進行離散譜分解，通過確定代表樹木年齡變化的諧波的周期數(shù)來估計樹木的年齡。

1 樣品采集

2017 年10 月，在山西省羊圈溝林場的華北落葉松人工林中，按徑階大、中、小至少各選3 株，優(yōu)勢木各選2 株，共52 株樣木進行針刺試驗。用針刺儀在每株樣木胸徑1.3 m 處4 個方向鉆入，獲取有效抗鉆阻力值序列208 組。2018 年6 月，在相同的樣木上，使用針刺儀在距離地面0.2 m 處和0.5 m 處分別從兩個不同方向鉆入，獲取有效抗鉆阻力數(shù)據(jù)115 組，兩次試驗共獲取有效抗鉆阻力數(shù)據(jù)323 組作為研究對象。

最后將樣木伐倒，在樹干上標明南北向，并分別在0.2 m、0.5 m 和1.3 m 針刺位置5 cm 內(nèi)截取圓盤，在圓盤非工作面上標明南北方向，并以分式形式注記，分子為樣地號和解析木號，分母為圓盤號和斷面高度，共獲取有效圓盤104 個。對圓盤進行拋光，打磨，掃描，用WinDENDRO 年輪分析系統(tǒng)結(jié)合人工判讀方式獲取圓盤年輪數(shù)，圓盤的年齡分布和徑階分布見表1 和表2。

表1 圓盤年齡分布匯總Table 1 Summery of disc age distribution

表2 圓盤徑階分布匯總Table 2 Summeryof disc radial distribution

2 利用頻譜分析測定樹木年齡的原理與方法

2.1 頻譜分析原理

假設每年生長近似相同，針刺儀鉆入活立木獲取的相對阻力剖面記錄為時間序列z={z1,z2,···,zk,···,zn}，n為針刺儀測量點的序列號，假設時間間隔為δ，則序列z是總時間長度為T=nδ的離散周期信號，以下將針刺儀抗鉆阻力序列稱為離散周期信號z。

離散時間信號z中 包含3 個部分的信息：趨勢、年度變化、噪音。利用頻譜分析估計樹木年齡的基本思路：將離散時間信號z去趨勢后，進行譜分析，確定代表年度變化的諧波信息。根據(jù)傅里葉分析可知，在時間T上，信號z可以分解成各種頻率三角函數(shù)的和，由于樹木生長1 年對應1 個峰，1 個周期內(nèi)有1 個峰，估計樹木的年齡取決于峰的數(shù)量，峰值和谷值分別代表早材和晚材。所以在這些三角函數(shù)中找到代表年度變化的周期數(shù)f，也就確定了周期P和年齡A，周期P在本研究中是指樹木生長1 年所對應的抗鉆阻力值序列長度，即由T=PA,T=f A，可得P=f。因次，估計樹木年齡關鍵是確定代表年度變化的周期數(shù)。

2.2 利用頻譜分析估計樹木年齡的方法

估計樹木年齡的基本思想是將抗鉆阻力序列利用平滑器去趨勢，然后再用傅里葉變換做頻譜分析，給出k次諧波的系數(shù)，求出k次諧波的振幅，樹木的年齡規(guī)律在振幅較大的諧波中得以體現(xiàn)。設定振幅比ε，尋找與最大振幅比值大于ε的振幅所對應的k次諧波的周期數(shù)k。

2.2.1 Savitzky-Golay 平滑器 Savitzky A 和Golay M 在1964 年提出了一種基于多項式擬合的方法來設計的形式簡單的濾波器[20]，對抗鉆阻力值所組成的離散周期信號z，把其中任一段數(shù)據(jù)測量點位置記為向量m=?M,···,0,···,M，測量點處對應的數(shù)據(jù)用一個p階多項式擬合：

共有2M+1點數(shù)據(jù)，多項式擬合階數(shù)為p，這種濾波器稱為Savitzky-Golay 平滑器，該平滑器允許窗長可以較大，且對于大部分數(shù)據(jù)的平滑都比較有效。使用該平滑濾波器對信號z濾波，可以求取趨勢項，用信號z減去趨勢項，就得到所期望的去趨勢后的信號，記為x={x0,x1,···,xk,···,xn?1}。

2.2.2 傅里葉級數(shù)與傅里葉變換 利用離散傅里葉級數(shù)，實現(xiàn)抗鉆阻力序列經(jīng)過去趨勢后所得到的新的序列x從時域到頻域的映射，利用離散的傅里葉變換，將長度為n的 時域序列x表示成n項虛指數(shù)信號的加權和，從而反映離散時間信號x中不同諧波分量的分布規(guī)律。

離散時間信號x的 離散傅里葉級數(shù)系數(shù)表示為：

離散時間信號x的逆離散傅里葉級數(shù)系數(shù)表示為：

2.2.3 離散時間信號的頻譜分析 針刺儀獲取的抗鉆阻力值序列z經(jīng)過去趨勢后得到離散周期信號x，記為：

其諧波分解（譜分解）表示為

系數(shù)的估計公式為：

當n=2q+1時，q=r，上面3 個估計系數(shù)的等式成立。

當n=2q時，r=q?1對于前q?1項系數(shù)，即k

用離散傅里葉級數(shù)系數(shù)表示離散譜分解系數(shù)，估計公式為：

用逆離散傅里葉級數(shù)系數(shù)表示離散譜分解系數(shù)，估計公式為：

記：α=(α1,α2,···,αq)，其中αi=Re(bi),i=1,2,···,q；β=(β1,β2,···,βq)，其中βi=Re(bi),i=1,2,···,q。

將第k諧波的振幅記為Mk，即

諧波周期數(shù)為k，說明波峰出現(xiàn)的次數(shù)為k。若代表年齡變換為第k次諧波，則樹木的估計年齡為諧波的次數(shù)k。將所有諧波的振幅及對應的周期數(shù)用矩陣表示：

矩陣P的每一行表示振幅和對應的周期數(shù)。

將向量M1和N1作為矩陣P1的第一列和第二列，即

將向量M2和N2作為矩陣P2的第一列和第二列，即

則矩陣P2是由矩陣P1的前j行組成，向量N2是與最大振幅比大于ε的所有的振幅按照從大到小排列組成。向量M2是每個振幅所對應的周期數(shù)。

將向量M2中所有分量取均值，記為δ，即

則δ表示樹木的年齡估計值。頻譜分析估計樹木年齡的流程圖（圖1）：

圖1 頻譜分析估計樹木年齡流程Fig.1 Flow chart of tree age estimation by spectrum analysis

2.2.4 步驟輸入：針刺儀鉆入活立木獲取的相對阻力值序列z={z1,z2,···,zk,···,zn}。

step1 給定窗口大小Wid，將針刺儀獲取的抗鉆阻力序列z為離散時間信號，經(jīng)過Savitzky-Golay平滑器，進行去趨勢后得到離散時間信號x。

step2 利用公式（3）求出逆傅里葉變換求出離散時間信號x的逆傅里葉級數(shù)系數(shù)。

step3 利用公式（11）求出離散時間信號x的諧波分解系數(shù)。

step4 求出所有諧波的振幅，并按照從大到小進行排列，與所對應的諧波的周期數(shù)作為行向量寫入矩陣P1。

step5 給定振幅比ε=0.85，將所有與最大振幅比值大于ε的振幅以及對應的周期數(shù)作為行寫入矩陣P2。

step6 求出矩陣P2的第一列向量的均值δ，作為樹木年齡的估計值。

輸出：樹木年齡估計值δ。

以上過程通過MATLAB 編程實現(xiàn)。

3 結(jié)果

對104 個圓盤所對應的抗鉆阻力值序列進行頻譜分析，首先根據(jù)圓盤直徑大小，選擇相應窗口值Wid，如表3 所示，對原始抗鉆阻力值序列進行去趨勢得到相應的離散周期信號，然后對其進行離散譜分析，計算得到所有諧波的振幅，將所有的振幅與對應的周期數(shù)作圖可以看出全頻譜分布情況，設定振幅比為0.85，求出與最大振幅比大于0.85的所有振幅都對應的周期數(shù)。以數(shù)據(jù)1 699 為例，數(shù)據(jù)1 699 是使用針刺儀在落葉松6 胸徑處獲取的數(shù)據(jù)，其胸徑為14.7 cm，設定相應窗口為501，進行去趨勢獲取相應的離散時間信號記為x1699，對x1699進行離散譜分析，可以得到所有諧波的振幅，諧波次數(shù)與對應振幅做直方圖可以看出全頻譜分布情況，見圖2。設定振幅比為0.85 后，計算出該數(shù)據(jù)第49 次諧波所對應的振幅滿足條件，說明第49 次諧波的周期變化代表樹木年齡變化，因第49 次諧波所對應的周期數(shù)為49，從而可以將周期數(shù)的一半作為該株落葉松胸徑處的估計年齡，為24.5 a，實測年齡為25 a，精度較好。

表3 參數(shù)Wid 選擇依據(jù)Table 3 The selection basis of parameter Wid

圖2 數(shù)據(jù)編號1 699 頻譜分析過程Fig.2 Spectrum analysis process of 1 699 data number

根據(jù)每個圓盤的直徑選擇相應的窗口大小，對原始抗鉆阻力值序列去趨勢后進行離散譜分析，可以得到每組抗鉆阻力值所對應的的周期數(shù)的平均數(shù)δ，試驗過程中，每個圓盤針刺獲取的抗鉆阻力值為2 組或者4 組，將對應的2 組或者4 組值的頻譜分析結(jié)果取均值，作為圓盤的頻譜分析算法估計年齡，所有圓盤頻譜分析算法估計年齡與實測年齡對比見圖3。

圖3 頻譜分析算法估計年齡與DECOM判定年齡對比Fig.3 Comparisons of age estimation by spectrum analysis algorithms and DECOM judgment test

每個圓盤的頻譜分析算法估計年齡與實測年齡的相對誤差分布見圖4，針刺儀自帶軟件DECOM自動判讀年齡結(jié)果誤差較大，誤差范圍是?25 a 至2 a 之間，平均誤差是?12 a，相對誤差大多集中在?20%至?60%之間，最小相對誤差為?7.69%，最大相對誤差達到?84.78%，平均相對誤差達到?49.98%。圓盤的實測年齡范圍是18～27 a，頻譜分析算法估計年齡誤差范圍是?5～6 a 之間，平均誤差是?0.25 a，平均絕對誤差是2 a；相對誤差分布大多集中在?10%至10%之間，最小相對誤差為0，最大相對誤差為25.69%，平均相對誤差為?0.35%。

圖4 頻譜分析算法估計年齡與DECOM 判定年齡絕對誤差分布Fig.4 Relative error distribution map of age estimation by spectrum analysis algorithms and DECOM judgment test

分別對真實年齡與軟件自動判別年齡（第1對）、真實年齡與算法估計年齡（第2 對）進行成對數(shù)據(jù)t 檢驗。第一對數(shù)據(jù)檢驗得到的t值為20.25，給定顯著性水平δ=0.05，查表可得t1?δ/2(n?1)=t0.975(9)=2.262 2，由于 |t|>2.262 2，故拒絕原假設，即可視為DECOM 判定樹木年齡均值與真實年齡均值之間有顯著差異，此時檢驗的p值為2.2 × 10?6。第二對數(shù)據(jù)檢驗得到t值為0.85，由于|t|<2.262 2，故不能拒絕原假設，說明頻譜分析算法估計樹木的年齡均值與樹木的真實年齡均值無顯著差異，此時檢驗的p值為0.394 9。

4 結(jié)論

通過給出的頻譜分析算法，將針刺儀獲取的華北落葉松抗鉆阻力值序列作為輸入，根據(jù)樹木的徑階給定窗口參數(shù)后，可得出活立木年齡的估計值，且精度較好，提高了針刺儀測定華北落葉松年齡的精度，改進了針刺儀自帶DECOM 軟件識別樹木年輪邊界準確率低，過于依賴人工經(jīng)驗的缺點，可以作為一種微損的估計華北落葉松年齡的有效方法，為活立木年齡的微損測定開辟了新的途徑。