基于自適應(yīng)濾波算法的樹木年輪測(cè)量方法

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

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

0 引言

樹木年齡是森林資源調(diào)查中的基本調(diào)查因子之一，是建立樹木生長模型、評(píng)價(jià)立地質(zhì)量、制定森林經(jīng)營方案的基礎(chǔ)數(shù)據(jù)。隨著樹木年輪學(xué)的不斷發(fā)展，樹木年輪學(xué)已廣泛應(yīng)用于林學(xué)、森林經(jīng)理學(xué)、生態(tài)學(xué)、氣候?qū)W、環(huán)境學(xué)等領(lǐng)域[1-4]。

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

1 原理與方法

1.1 樹木年輪測(cè)量儀機(jī)械結(jié)構(gòu)和工作原理

年輪測(cè)量儀由2個(gè)電機(jī)控制，一個(gè)電機(jī)控制鉆針的旋轉(zhuǎn)，另一個(gè)電機(jī)控制鉆針的前進(jìn)。控制鉆針旋轉(zhuǎn)的電機(jī)稱為旋轉(zhuǎn)電機(jī)，旋轉(zhuǎn)電機(jī)一般使用直流電機(jī)?？刂沏@針前進(jìn)的電機(jī)稱為推進(jìn)電機(jī)，推進(jìn)電機(jī)一般使用步進(jìn)電機(jī)。鉆針通過鉆針夾與旋轉(zhuǎn)電機(jī)軸連接，旋轉(zhuǎn)電機(jī)固定在滑動(dòng)基座上，推進(jìn)電機(jī)軸與傳動(dòng)絲桿連接，滑動(dòng)基座與傳動(dòng)絲桿連接。當(dāng)推進(jìn)電機(jī)旋轉(zhuǎn)時(shí)，帶動(dòng)傳動(dòng)絲桿轉(zhuǎn)動(dòng)，從而帶動(dòng)滑動(dòng)基座在傳動(dòng)絲桿上移動(dòng)。年輪測(cè)量儀機(jī)械傳動(dòng)結(jié)構(gòu)圖如圖1所示，實(shí)物圖如圖2所示，試驗(yàn)操作圖如圖3所示。

圖1 年輪測(cè)量儀機(jī)械結(jié)構(gòu)圖

圖2 年輪測(cè)量儀實(shí)物圖

圖3 年輪測(cè)量儀操作圖

當(dāng)旋轉(zhuǎn)電機(jī)恒速旋轉(zhuǎn)時(shí)，鉆針阻力與旋轉(zhuǎn)電機(jī)電樞電流成正比[18]。因此，可利用旋轉(zhuǎn)電機(jī)電樞電流表示鉆針阻力。鉆針為德國Rinntech公司生產(chǎn)的樹木針刺儀鉆針[7]，鉆針形狀如圖4所示。鉆針針頭寬3 mm，鉆針針桿直徑1.5 mm，因此，鉆針阻力主要集中在鉆針針頭上。鉆針阻力與鉆針針頭切割處的樹木微密度正相關(guān)。當(dāng)鉆針鉆入晚材部分時(shí)，晚材密度較大，鉆針?biāo)茏枇^大，因而，旋轉(zhuǎn)電機(jī)電流較大；當(dāng)鉆針鉆入早材部分時(shí)，早材密度較小，鉆針?biāo)茏枇^小，因而，旋轉(zhuǎn)電機(jī)電流較小。因此，可以通過旋轉(zhuǎn)電機(jī)電流變化情況確定樹木年輪。當(dāng)鉆針鉆入樹木方向與樹木年輪線切線方向垂直時(shí)，鉆針切割晚材的面積最大，且切割面積平穩(wěn)變化，鉆針?biāo)茏枇ζ椒€(wěn)變化，鉆針切割樹木晚材部分時(shí)所受的阻力最大。

圖4 鉆針形狀示意圖

1.2 直流電機(jī)電流測(cè)量方法

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

圖5 λ對(duì)電流限幅預(yù)處理的影響

1.3 試驗(yàn)樣本與試驗(yàn)方法

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

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

2 數(shù)據(jù)處理

2.1 直流電機(jī)電流限幅預(yù)處理

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

2.2 低通濾波器對(duì)限幅后電流信號(hào)的濾波處理

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

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

(1)

其中

α=2πfl/fs

(2)

式中k——采樣序號(hào)

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

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

α——濾波系數(shù)

fl——濾波器截止頻率

fs——原始信號(hào)采樣頻率

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

2.3 自適應(yīng)低通濾波器對(duì)限幅后電流信號(hào)的濾波處理

普通低通濾波器調(diào)試好后，濾波系數(shù)不變，濾波器的截止頻率也不變。自適應(yīng)濾波器根據(jù)原始信號(hào)變化幅度和濾波后輸出信號(hào)變化幅度，自動(dòng)更新低通濾波器濾波系數(shù)。當(dāng)輸入值變化幅度增大時(shí)，減小濾波系數(shù)，以提高濾波器平穩(wěn)度；當(dāng)濾波后輸出值變化幅度增大時(shí)，增大濾波系數(shù)，以提高濾波器的靈敏度。自適應(yīng)低通濾波器濾波系數(shù)α的計(jì)算公式為

α=

(3)

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

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

b1——濾波器輸入值影響基數(shù)

b2——濾波器輸出值影響基數(shù)

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

3 試驗(yàn)與結(jié)果分析

在120個(gè)樣本數(shù)據(jù)中，按照樹木徑階大、中、小分別統(tǒng)計(jì)測(cè)量路徑處的去皮樹干平均直徑、測(cè)量路徑處的平均年輪數(shù)、平均年輪寬度以及普通低通濾波器和自適應(yīng)低通濾波器識(shí)別年輪的正確率，統(tǒng)計(jì)結(jié)果如表1所示。

表1 年輪分析結(jié)果統(tǒng)計(jì)

3.1 大徑階圓盤

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

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

3.2 中徑階圓盤

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

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

3.3 小徑階圓盤

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

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

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

4 結(jié)論

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

(2)提出的自適應(yīng)低通濾波算法可通過自動(dòng)更新濾波系數(shù)來調(diào)節(jié)低通濾波器的濾波特性。當(dāng)?shù)屯V波器的輸入值變化較大時(shí)，減小濾波系數(shù)，提高濾波器的平穩(wěn)度；當(dāng)?shù)屯V波器的輸出值變化較大時(shí)，增加濾波系數(shù)，提高濾波器的靈敏度。因此，自適應(yīng)低通濾波算法具有很強(qiáng)的自適應(yīng)能力。

(3)對(duì)年輪測(cè)量儀直流電機(jī)電流進(jìn)行濾波處理，使用自適應(yīng)低通濾波器的年輪平均識(shí)別正確率達(dá)到89.66%，普通低通濾波器的年輪平均識(shí)別正確率為75.84%。

(4)當(dāng)樹木年輪寬度變窄時(shí)，自適應(yīng)低通濾波器的年輪識(shí)別正確率降低了13.59個(gè)百分點(diǎn)，而普通低通濾波器的年輪識(shí)別正確率降低了24.07個(gè)百分點(diǎn)。因此，自適應(yīng)低通濾波器的年輪識(shí)別正確率比普通低通濾波器穩(wěn)定。