基于自適應濾波算法的樹木年輪測量方法

姚建峰 趙燕東 盧 軍 鄭一力 高瑞東 唐守正

(1.中國林業科學研究院資源信息研究所，北京 100091；2.信陽師范學院計算機與信息技術學院，信陽 464000；3.北京林業大學工學院，北京 100083；4.山西省管涔山國有林管理局，忻州 034000)

0 引言

樹木年齡是森林資源調查中的基本調查因子之一，是建立樹木生長模型、評價立地質量、制定森林經營方案的基礎數據。隨著樹木年輪學的不斷發展，樹木年輪學已廣泛應用于林學、森林經理學、生態學、氣候學、環境學等領域[1-4]。

目前，測量樹木年輪最常用的方法是生長錐法。生長錐法測量精度高，但是取樣樹芯困難、速度慢，且對樹木生長有一定的負面影響[5]。本團隊設計的樹木年輪測量儀可微損、快速測量樹木年輪，其機械結構類似于德國Rinntech公司生產的Resistograph系列樹木針刺儀的結構[6-7]，通過直流電機帶動鉆針鉆入樹木，并用直流電機電樞電流表示鉆針阻力。由于直流電機電流信號中含有大量的高頻噪聲信號，使直流電機原始電流曲線圖不能直觀地反映樹木年輪信息，因此，需要對直流電機電流進行濾波處理。普通低通濾波器濾波參數固定[8]，當樹木年輪寬度發生變化時，使有效年輪信號頻率、噪聲信號頻率發生變化，因此，濾波后的電流信號中既可能保留了較多的噪聲信號，同時又可能濾去有效的樹木年輪信號。因此，普通低通濾波器不適合用于年輪測量儀直流電機電流信號的濾波處理。自適應低通濾波器可以通過自動更新濾波系數來調節濾波特性，以達到濾波要求，具有很強的自適應能力，且不需要精確的數學模型，已經廣泛應用于系統識別、預測、噪聲消除和反演建模等領域[9-15]。本文利用自適應算法的思想[16-17]，在普通低通濾波算法的基礎上，提出一種自適應低通濾波算法，根據原始信號變化幅度和濾波后輸出信號變化幅度自動更新低通濾波器濾波系數。

1 原理與方法

1.1 樹木年輪測量儀機械結構和工作原理

年輪測量儀由2個電機控制，一個電機控制鉆針的旋轉，另一個電機控制鉆針的前進。控制鉆針旋轉的電機稱為旋轉電機，旋轉電機一般使用直流電機。控制鉆針前進的電機稱為推進電機，推進電機一般使用步進電機。鉆針通過鉆針夾與旋轉電機軸連接，旋轉電機固定在滑動基座上，推進電機軸與傳動絲桿連接，滑動基座與傳動絲桿連接。當推進電機旋轉時，帶動傳動絲桿轉動，從而帶動滑動基座在傳動絲桿上移動。年輪測量儀機械傳動結構圖如圖1所示，實物圖如圖2所示，試驗操作圖如圖3所示。

圖1 年輪測量儀機械結構圖

圖2 年輪測量儀實物圖

圖3 年輪測量儀操作圖

當旋轉電機恒速旋轉時，鉆針阻力與旋轉電機電樞電流成正比[18]。因此，可利用旋轉電機電樞電流表示鉆針阻力。鉆針為德國Rinntech公司生產的樹木針刺儀鉆針[7]，鉆針形狀如圖4所示。鉆針針頭寬3 mm，鉆針針桿直徑1.5 mm，因此，鉆針阻力主要集中在鉆針針頭上。鉆針阻力與鉆針針頭切割處的樹木微密度正相關。當鉆針鉆入晚材部分時，晚材密度較大，鉆針所受阻力較大，因而，旋轉電機電流較大；當鉆針鉆入早材部分時，早材密度較小，鉆針所受阻力較小，因而，旋轉電機電流較小。因此，可以通過旋轉電機電流變化情況確定樹木年輪。當鉆針鉆入樹木方向與樹木年輪線切線方向垂直時，鉆針切割晚材的面積最大，且切割面積平穩變化，鉆針所受阻力平穩變化，鉆針切割樹木晚材部分時所受的阻力最大。

圖4 鉆針形狀示意圖

1.2 直流電機電流測量方法

年輪測量儀采用串聯電阻法測量電機電樞回路電流。在電機電樞回路中串聯一個10 mΩ高精密采樣電阻，根據歐姆定理，電樞回路電流與采樣電阻兩端的電壓成正比。采樣電阻兩端的電壓信號經放大器放大后，輸入到DSP(Digital signal processing)處理器的A/D(Analog/Digital)轉換接口中。DSP處理器的A/D轉換器是12位，轉換后的數據范圍是[0,4 095]。

圖5 λ對電流限幅預處理的影響

1.3 試驗樣本與試驗方法

試驗樣本是2018年6月在山西省忻州市五寨縣羊圈溝林場華北落葉松人工林采集的圓盤。該落葉松人工林造林時間是1988年，因此，采樣時樹齡為30 a。按樹木直徑大(胸徑大于等于16 cm)、中(胸徑在10～16 cm)、小(胸徑小于等于10 cm)各選10株，樣木隨機分布在落葉松人工林內，盡量涵蓋不同的坡向和海拔。在樹干高0.5 m和1.3 m處各截取厚度為5～7 cm圓盤，共截取60個圓盤。

用年輪測量儀對60個華北落葉松圓盤進行測試，每個圓盤分別測量2次。本次試驗電流采樣周期是1 ms,鉆針前進速度是10 cm/min，因此2個連續采樣點之間的距離是0.001 67 mm。由于鉆針路徑大部分沒有經過髓心，一般要偏離髓心2～4個年輪，所以測量時記錄測量路徑去皮樹干直徑和測量路徑年輪數。先調試電流限幅預處理系數，使電流經過限幅預處理后能濾去原始電流中的突變電流噪聲信號。然后調試普通低通濾波器濾波系數α，使濾波器對所有經過預處理后的電流信號濾波后，電流波形圖能反映出70%左右的圓盤年輪信息，且濾波后電流波形圖變化平滑，毛刺較少。最后調試自適應低通濾波器中各參數的系數，進一步優化濾波效果。

2 數據處理

2.1 直流電機電流限幅預處理

電機原始電流信號中存在時間極短、幅度變化很大的突變噪聲信號。為了消除突變噪聲信號，在對電流濾波處理之前，需要對電流進行限幅預處理，限制兩次相鄰采樣電流的最大變化幅度[19]。限幅預處理是把第k次電流采樣值x(k)與第k-1次采樣值x(k-1)相減，求出其變化量。如果變化量的絕對值大于λx(k-1)(λ為限幅系數)，則變化量取λx(k-1)。λ取值非常重要，如果λ取值過大，突變噪聲信號無法得到有效抑制，如果λ取值過小，可能濾除部分變化比較大的有效電流信號。圖5是用年輪測量儀測量一個直徑為14 cm圓盤的直流電機電流曲線圖，橫坐標表示鉆針前進路徑的長度L，縱坐標表示直流電機的電流I。從圖5可以看出，當λ取0.1時，限幅預處理后的電流仍含有大量的突變噪聲信號；當λ取0.001時，突變噪聲信號得到有效抑制，但是當電機真實電流信號急劇上升或者急劇下降時，限幅預處理后的電流信號變化緩慢。經大量測試，本試驗λ取0.01比較合適。

2.2 低通濾波器對限幅后電流信號的濾波處理

低通濾波算法是將硬件RC(Resistance capacitance)低通濾波器的微分方程用差分方程來表示[20-22]，其數學表達式為

y(k)=(1-α)y(k-1)+αx(k)

(1)

其中

α=2πfl/fs

(2)

式中k——采樣序號

x(k)——濾波器第k次濾波輸入值

y(k)——濾波器第k次濾波輸出值

α——濾波系數

fl——濾波器截止頻率

fs——原始信號采樣頻率

在大部分應用領域，噪聲信號的頻率和有效信號的頻率是動態變化的，所以很難確定濾波器的截止頻率。通常情況下，一般根據應用系統的需要調試參數α，使濾波效果滿足應用系統的要求。α的取值范圍是[0，1]，α越小，截止頻率越低，輸入值x(k)在輸出值y(k)中的比例越小，濾波器平穩度越高，靈敏度越低，濾波效果越好。但當α取值過小、濾波器靈敏度過低時，會產生濾波器的輸出滯后于輸入的變化趨勢。對落葉松圓盤的限幅預處理后的電流信號進行濾波分析，低通濾波器中α取0.01比較合適。

2.3 自適應低通濾波器對限幅后電流信號的濾波處理

普通低通濾波器調試好后，濾波系數不變，濾波器的截止頻率也不變。自適應濾波器根據原始信號變化幅度和濾波后輸出信號變化幅度，自動更新低通濾波器濾波系數。當輸入值變化幅度增大時，減小濾波系數，以提高濾波器平穩度；當濾波后輸出值變化幅度增大時，增大濾波系數，以提高濾波器的靈敏度。自適應低通濾波器濾波系數α的計算公式為

α=

(3)

式中λ1——濾波器輸入值變化量影響因子

λ2——濾波器輸出值變化量影響因子

b1——濾波器輸入值影響基數

b2——濾波器輸出值影響基數

當λ1越大時，α受濾波器輸入值變化量影響越大；當λ2越大時，α受濾波器輸出值變化量影響越大。b1、b2的數量級與輸入量變化量、輸出量變化量的數量級相近，且要求α0=b1/(b0+b1)，α0是沒有使用自適應濾波算法時調試好的濾波系數。經過大量測試，對各參數取不同的值，對華北落葉松圓盤進行濾波處理，對比各參數取不同值時濾波效果的差異，最后確定λ1=0.5、λ2=8、b1=99、b2=1。

3 試驗與結果分析

在120個樣本數據中，按照樹木徑階大、中、小分別統計測量路徑處的去皮樹干平均直徑、測量路徑處的平均年輪數、平均年輪寬度以及普通低通濾波器和自適應低通濾波器識別年輪的正確率，統計結果如表1所示。

表1 年輪分析結果統計

3.1 大徑階圓盤

對于大徑階圓盤，自適應低通濾波器的年輪識別正確率是95.61%，普通低通濾波器的年輪識別正確率是86.37%。盡管自適應低通濾波器的正確率比普通低通濾波器的正確率高9.24個百分點，但是普通低通濾波器的正確率也能達到86.37%，基本上能識別大部分年輪。主要因為對于大徑階圓盤，年輪比較寬，平均年輪寬度達3.8 mm，最窄年輪寬度達2 mm，由于年輪測量儀每毫米采樣600個電流數據，最窄的年輪也能采樣1 200個數據。由于本試驗普通低通濾波器的濾波系數α為0.01，所以，只需要100次采樣就能從電流的最小值跳變到最大值。所以，即使是普通低通濾波器，也能識別大部分的樹木年輪。對于大徑階圓盤，不能識別的年輪主要位于圓盤髓心附近。髓心部分年輪不容易識別的原因如下：① 髓心附近，樹木的晚材寬度比較窄，晚材密度比較小，導致鉆針鉆過晚材部分的電流值變化不大，波峰不明顯。② 一般情況下，鉆針路徑很難通過樹木髓心，越靠近髓心，鉆針路徑偏離年輪線切線的垂直方向越遠，導致鉆針鉆過晚材部分時阻力發生波動，電流值也發生波動。以一個直徑為172 mm的圓盤為例，分別分析普通低通濾波器濾波后的電流曲線圖和自適應低通濾波器濾波后的電流曲線圖與圓盤年輪的匹配情況，從而顯示自適應低通濾波器和普通低通濾波器對大徑階圓盤的濾波效果差異。圖6a為自適應低通濾波器的年輪分析圖，圖6b為普通低通濾波器的年輪分析圖。從圖6a中可以看出，在髓心附近，自適應低通濾波器濾波后的電流變化幅值較大，分析出來的年輪數可能比實際年輪數大；從圖6b中可以看出，在髓心附近，普通低通濾波器濾波后的波峰與波谷變化幅值較小，分析出來的年輪數可能比實際年輪數小。該圓盤右半部分沒有樹皮，當鉆針鉆出圓盤時，圓盤邊緣部分破裂，導致最后一個年輪不能識別。

圖6 大徑階圓盤年輪分析示例圖

3.2 中徑階圓盤

對于中徑階圓盤，自適應低通濾波器的年輪識別正確率是91.36%，普通低通濾波器的年輪識別正確率是78.84%。自適應低通濾波器的正確率比普通低通濾波器的正確率高12.52個百分點。自適應低通濾波器的識別正確率比較高，但普通低通濾波器的識別正確率偏低。盡管中徑階圓盤的年輪平均寬度達2.80 mm，但是由于有些年份氣候異常(如干旱)，有些年輪的寬度小于1 mm，導致普通低通濾波器難以識別。以一個直徑為121 mm的圓盤為例，分別分析普通低通濾波器濾波后的電流曲線圖和自適應低通濾波器濾波后的電流曲線圖與圓盤年輪的匹配情況，顯示自適應低通濾波器和普通低通濾波器對中徑階圓盤的濾波效果差異。圖7a為自適應低通濾波器濾波后的電流信號年輪分析圖，圖7b為普通低通濾波器濾波后的電流信號年輪分析圖。在圖7中，標有1、2的位置處，在自適應低通濾波器濾波后的電流波形圖中有較為明顯的波峰，而普通低通濾波器不明顯，從而導致普通低通濾波器年輪識別正確率偏低。

圖7 中徑階圓盤年輪分析示例圖

3.3 小徑階圓盤

圖8 小徑階圓盤年輪分析示例圖

對于小徑階圓盤，自適應低通濾波器的年輪識別正確率是82.02%，普通低通濾波器的年輪識別正確率是62.30%。盡管自適應低通濾波器的正確率比普通的低通濾波算法的正確率高19.72個百分點，但自適應低通濾波器的正確率仍偏低。對于小徑階樹木，當森林郁閉后，小徑階樹木是被壓木，生長緩慢，邊材部分的平均年輪寬度小于1 mm。由于鉆針針頭有0.3 mm長的針尖，當鉆針針頭沒有鉆出上一個年輪的晚材部分時，鉆針針尖已經進入了下一個年輪的晚材，再加上鉆針鉆入方向可能偏離年輪線切線方向的垂直方向、鉆針抖動等原因，使鉆針針尖阻力不能反映窄年輪的早、晚材之間的密度變化，導致2種濾波算法的年輪識別正確率都偏低。以一個直徑為70 mm的圓盤為例，分別分析普通低通濾波器濾波后的電流曲線圖和自適應低通濾波器濾波后的電流曲線圖與圓盤年輪的匹配情況，顯示自適應低通濾波器和普通低通濾波器對小徑階圓盤濾波效果的差異。圖8a為自適應低通濾波器年輪分析圖，圖8b為普通低通濾波器年輪分析圖。在圖8中，在標有1的位置處，自適應低通濾波器濾波后的電流波形圖中有小的波峰，而普通低通濾波器則幾乎沒有波峰，所以自適應低通濾波器的年輪識別正確率高于普通低通濾波器；在標有2的位置處，2種濾波算法都不能識別出有效的年輪信號，所以，對于小徑階圓盤，2種濾波算法的年輪識別正確率都偏低。

當平均年輪寬度從3.8 mm(大徑階圓盤)下降到1.9 mm(小徑階圓盤)時，自適應低通濾波器的年輪識別正確率從95.61%下降到82.02%，下降了13.59個百分點，而普通低通濾波器的年輪識別正確率從86.37%下降到62.30%，下降了24.07個百分點。

4 結論

(1)設計了一種便攜式樹木年輪微損測量儀，通過落葉松圓盤試驗驗證了該設備的可行性。

(2)提出的自適應低通濾波算法可通過自動更新濾波系數來調節低通濾波器的濾波特性。當低通濾波器的輸入值變化較大時，減小濾波系數，提高濾波器的平穩度；當低通濾波器的輸出值變化較大時，增加濾波系數，提高濾波器的靈敏度。因此，自適應低通濾波算法具有很強的自適應能力。

(3)對年輪測量儀直流電機電流進行濾波處理，使用自適應低通濾波器的年輪平均識別正確率達到89.66%，普通低通濾波器的年輪平均識別正確率為75.84%。

(4)當樹木年輪寬度變窄時，自適應低通濾波器的年輪識別正確率降低了13.59個百分點，而普通低通濾波器的年輪識別正確率降低了24.07個百分點。因此，自適應低通濾波器的年輪識別正確率比普通低通濾波器穩定。