基于針刺儀平穩卡爾曼濾波器的樹木年齡估計*

潘 虹 盧 軍 雷相東 郭旭展 姚建峰 唐守正

(1.中國林業科學研究院資源信息研究所 北京 100091；2.信陽師范學院數學與統計學院 信陽 464000)

樹木年齡是森林資源調查的重要因子，可為確定合理采伐量、制定林業規劃和更加科學的森林經營措施等提供重要依據，在考古學、氣候學、災害學等方面也有廣泛應用(Bollschweileretal.，2010；D′Arrigoetal.，2010；Fangetal.，2012)。活立木年齡微損測量是目前林業生產和研究中最重要、最基礎的共性難題，是解決很多林業問題的關鍵，如預測木材產量和森林碳儲量、確定樹木經濟價值、制定古樹保護措施以及古樹保護等級劃分等(Halletal.，2004；Thurigetal.，2005；Micheletal.，2009)，對森林管護、森林經營和古樹保護等均具有重要的理論和現實意義。

確定樹木年齡最準確的方法是識別樹木年輪數(Wongetal.，2001；Rozas，2003)，通常采用生長錐測定法、解析木圓盤測定法、顯微鏡獲取年輪法等，但這些方法取樣基于生長芯或樹干橫截面，會破壞樹木生長。無損的樹木年齡確定方法有查數輪生枝法(吳斡寧等，2013)、建立線性生長模型(Careyetal.，2000；Loewensteinetal.，2000；Trotsiuketal.，2012)或非線性數學模型(Oberhuberetal.，2000；Rohneretal.，2013；潘虹等，2020)，但該類方法不具有普遍性且準確性也較低(Loewensteinetal.，2000；Rohneretal.，2013)?；谀壳皹淠灸挲g確定方法的不足以及我國現行天然林保護的大環境，尋求一種微損測定樹木年齡的方法是亟待解決的問題。

德國RINNTECH公司生產的針刺儀可以估計樹木年齡、樹木生長率(Orozco-Aguilaretal.，2018；Ohetal.，2019)以及腐朽程度(岳小泉，2017)。針刺儀測定樹木年齡無需樣本處理，對活立木傷害小，是微損的；但由于研究環境、操作者、儀器等因素影響針刺儀測量的靈敏度(Ukrainetzetal.，2010)，觀測數據中存在干擾噪聲，利用針刺儀自帶DECOM軟件檢測年輪邊界時(Chantreetal.，1997)很多真實年輪邊界不能準確識別，使得針刺儀應用于樹木年輪微損測定具有一定的局限性。而對針刺儀獲取的抗鉆阻力值序列進行數據分析，快速直接地估計樹木年齡，可以推進樹木年齡微損測定的研究進程(潘虹等，2021)。

卡爾曼濾波器本質上是一組數學方程，為最小二乘法提供有效的遞歸解，支持過去、現在和未來狀態的估計?？柭鼮V波是統計估計理論的里程碑，是處理噪聲的有利武器(Kalman，1960)，可應用于白噪聲激勵的任何平穩或非平穩隨機向量的估計，所得估計在線性估計中精度最佳(黃小平等，2015)。利用卡爾曼濾波器的去噪功能，處理抗鉆阻力值序列的干擾噪聲，可為更準確地估計樹木年齡做好數據準備。

鑒于此，本研究以德國RINNTECH公司生產的針刺儀(Resistograph 4452P/S)鉆入活力木獲取的抗鉆阻力值序列為研究樣本，基于卡爾曼濾波器時間更新方程和觀測更新方程，在理論上推導出針刺儀平穩卡爾曼濾波器，將平穩卡爾曼濾波器應用于活力木抗鉆阻力值序列估計活立木年齡，并通過大量數據驗證算法的準確性，以期為活立木年齡估計提供方法和依據。

1 基本思想

針刺儀鉆入活立木獲得的序列號用進針深度表示，由于進針是勻速的，因此可將序列號看作進針時刻，這樣抗鉆阻力值序列就可以看作一組離散時間信號，每株活立木都對應一組離散時間信號。直觀上來看，信號波峰波谷振幅差別較大，有很多振幅很小的波峰波谷，并不是所有成對出現的波峰波谷均能代表年輪邊界，如圖1所示。這說明原始信號存在噪聲干擾，去除不代表年輪變化的波峰波谷，也就是去除離散時間信號的噪聲，成為判斷樹木年輪邊界的一個關鍵點。

圖1 活立木抗鉆阻力剖面與樹芯年輪對比Fig.1 Comparison of drilling resistance profile and core ring of standing tree

為了利用卡爾曼濾波器的去噪功能處理針刺儀抗鉆阻力值序列的干擾噪聲，本研究在卡爾曼濾波器原理的基礎上，將針刺儀看作一個系統，給出針刺儀卡爾曼濾波器，同時根據其特有性質進行理論推導，提出針刺儀簡單卡爾曼濾波器、針刺儀平穩卡爾曼濾波器。其基本思想是采用針刺儀信號與噪聲的狀態空間模型，利用針刺儀抗鉆阻力前一時刻的估計值和現時刻的觀測值更新對狀態變量的估計，求出抗鉆阻力現在時刻的最優估計值。

2 針刺儀卡爾曼濾波算法推導

2.1 針刺儀卡爾曼濾波器

將針刺儀看作一個系統，針刺儀鉆入活立木獲取的相對阻力剖面記錄為序列z={z1,z2,z3,···,zk,···,zn}，n為針刺儀測量點的序列號，也可以看作n時刻(潘虹等，2021)。這是一組離散時間信號，信號中第k時刻的抗鉆阻力觀測值為zk。設針刺儀系統的狀態方程和觀測方程分別為：

xk=Axk-1+wk-1；

(1)

zk=Hxk+vk。

(2)

式中：xk∈Rn為第k時刻的抗鉆阻力狀態向量；zk∈Rm為第k時刻的抗鉆阻力實測值；wk為隨機噪聲向量；vk為觀測噪聲向量。稱式(1)為針刺儀的狀態方程、式(2)為針刺儀的觀測方程，稱A為狀態轉移矩陣、H為觀測矩陣。

基本假設1：隨機信號wk和vk分別表示針刺儀過程激勵噪聲和觀測噪聲，假設wk和vk相互獨立，且wk-1～N(0,Q)，vk-1～N(0,R)。

基本假設2：設相對抗鉆阻力的初始狀態值為x0，且x0不相關于wk和vk，

E(x0)=μ0,E[(x0-μ0)(x0-μ0)T]=P0。

(3)

式中：E表示均值。

基本假設3：設針刺儀的狀態轉移矩陣和觀測矩陣均為1，即A=1、H=1。

(4)

(5)

定理1(針刺儀卡爾曼濾波器) 在式(1)和(2)以及基本假設1～3下，遞推針刺儀卡爾曼濾波器如下。

1)針刺儀時間更新方程——向前推算當前時刻抗鉆阻力值誤差協方差估計的值，構造下一時刻抗鉆阻力值的先驗估計：

(6)

(7)

2)針刺儀觀測更新方程——根據抗鉆阻力先驗估計值和新的抗鉆阻力實測值得到改進后的抗鉆阻力后驗估計值：

(8)

(9)

(10)

圖2 針刺儀卡爾曼濾波器算法流程Fig.2 Kalman filter flow chart of calculating k time estimated value from k-1 time estimated value

2.2 簡化針刺儀卡爾曼濾波器

(11)

(12)

(13)

(14)

(15)

定理2(簡化針刺儀卡爾曼濾波器) 在式(1)和(2)以及基本假設1～3下，簡化針刺儀卡爾曼濾波器遞推公式如下。

(16)

(17)

(18)

2.3 針刺儀平穩卡爾曼濾波器

在簡化針刺儀卡爾曼濾波中，抗鉆阻力后驗估計值和后驗協方差的計算僅依賴于后驗信息，去掉其前驗估計值和前驗協方差的自回歸計算，使得繁瑣的計算過程得到簡化，但仍需要設置4個初始參數：Q、R、x0、P0，每個參數與最終結果之間的關系難以掌握。因此考慮從后驗協方差存在上下確界這條性質入手，進一步研究簡化針刺儀卡爾曼濾波器，提出針刺儀平穩卡爾曼濾波器的概念，簡化參數輸入，且給出抗鉆阻力值序列最優估計值的準確表達式。

(19)

將式(19)代入式(12)，可得后驗估計值的另外一個表達式：

(20)

定理3(平穩卡爾曼濾波器) 在式(1)和(2)以及基本假設1～3下，假設初始條件為：

簡化遞推卡爾曼濾波器可以表示為：

(21)

(22)

(23)

定理4 平穩卡爾曼濾波器等價于離散差分方程：

x[k]-sx[k-1]=(1-s)z[k]。

(24)

其解為：

(25)

移項得：

得到平穩卡爾曼濾波器的離散差分方程如下：

x[k]-sx[k-1]=(1-s)z[k]。

利用卷積法，求以上常系數差分方程的解為式(25)。

推論 針刺儀測量阻力值{z1,z2,…,zn}經平穩卡爾曼濾波器處理后，得到相應的最優解：

(26)

2.4 算法步驟

第一步，賦初始值,Q=1,R=Rat，Rat為常數，

(27)

第二步，對于k=2,3,4,…,n，令

s=(P+Q)(P+Q+R)-1；

(28)

(29)

依次計算zk的后驗估計值：

(30)

以上算法步驟通過MATLAB編程實現，相應的偽代碼如下：

輸入：研究樣本z={z1,z2,···,zn}；參數Rat；

過程：

2)fori=2,3,…,ndo；

4)end；

5)賦初始值t=0；

6)forl=2,3,…,ndo；

7)賦初始值tmax=0,tmin=0；

11)if |tmax-tmin|>0；

12)t=t+1；

13)end；

14)end；

15)end。

輸出：樣本集峰點個數t。

3 案例分析

3.1 數據獲取

以山西省羊圈溝林場華北落葉松(Larixprincipis-rupprechtii)為研究對象。2017年10月，在林場華北落葉松人工林(初植密度3 300 株·hm-2)中按徑階大、中、小至少各選3株，優勢木各選2株，共52株樣木進行針刺試驗。同一樣區相同徑階樣木按坡位上、中、下至少各選1株用噴漆標記編號，并詳細記錄樣木相對位置，所選取樣木盡可能樹干通直、飽滿、無病蟲害、不斷梢且生長良好。在每株樣木胸徑1.3 m處用針刺儀從4個方向鉆入，獲取有效抗鉆阻力數據208組。2018年6月，在相同樣木上用針刺儀在0.2 m和0.5 m處從2個方向鉆入，獲取有效抗鉆阻力數據115組。2次試驗共獲取有效抗鉆阻力數據323組作為研究樣本。伐倒針刺樣木，在針刺位置5 cm內截取圓盤，共獲取104個有效圓盤作為參考樣本，圓盤進行拋光、打磨、掃描，用WinDENDRO年輪分析系統結合人工判讀方式獲取圓盤年輪數(潘虹等，2020)。

3.2 試驗結果

以針刺儀鉆入活立木獲取的抗鉆阻力數據作為輸入，經過針刺儀平穩卡爾曼濾波器后，輸出為相應的抗鉆阻力最優估計值序列，相比抗鉆阻力原始數據，可較好地去除噪聲，原來較小的波峰波谷進行了平滑處理，且保持了原始數據的波峰波谷走勢。以數據編號1702的阻力折線圖為例，將原始阻力數據折線圖與濾波后的折線圖進行對比，并將折線的部分圖進行放大，見圖3。

圖3 原始抗鉆阻力值與濾波后最優估計值的差值比較Fig.3 Comparison of the difference between the original drilling resistance value and the optimal estimated value after filtering

對104個圓盤進行試驗，起始圓盤直徑為8 cm，以6 cm為1個徑階，選擇相應參數Rat值的依據見表1。

表1 參數Rat選擇依據Tab.1 The basis of parameter Rat selection

根據每個圓盤直徑選擇相應的參數Rat值，每組抗鉆阻力值經過平穩卡爾曼濾波器后，成對波峰波谷數量減少，且與圓盤實際年輪數較接近。從323組抗鉆阻力值中隨機選取20組，試驗結果見表2。試驗過程中，每個圓盤針刺獲取的抗鉆阻力值為2組或4組，將對應的2組或4組值濾波后的波峰波谷數取平均值，作為圓盤的濾波算法估計年齡。詳細的試驗結果和實際年齡對比如圖4所示。

圖4 濾波算法估計年齡與DECOM自動判定年齡對比Fig.4 Comparisons of age estimation by filtering algorithms and DECOM judgment test

表2 隨機選取20組抗鉆阻力值序列濾波前后峰谷數對比Tab.2 Comparison of peak-valley numbers before and after filtering for 20 groups of random drilling resistance values

每個圓盤的試驗相對誤差分布見圖5，針刺儀自帶軟件DECOM自動判定年齡結果誤差較大，范圍在-25～2年之間，平均誤差為-12年；相對誤差大多集中在-20%～-60%之間，最小相對誤差為-7.69%，最大相對誤差為-84.78%，平均相對誤差為-40.49%。算法估計年齡誤差范圍在-6～5年之間，平均誤差為-0.25年；相對誤差分布大多集中在-10%～10%之間，最小相對誤差為0%，最大相對誤差為25.69%，平均相對誤差為0.75%。

圖5 濾波算法估計年齡與DECOM自動判定年齡相對誤差分布Fig.5 Relative error distribution of age estimation by filtering algorithms and DECOM judgment test

分別對實際年齡與軟件自動判定年齡(第1對)、實際年齡與算法估計年齡(第2對)進行成對數據t-檢驗。第1對數據檢驗后t為t1=20.254，給定顯著性水平δ=0.05，查表可得t1-δ/2(n-1)=t0.975(9)=2.262 2，由于|t1|>2.262 2，故拒絕原假設，即可認為DECOM軟件自動判定樹木年齡均值與實際年齡均值之間有顯著差異，此時檢驗的p為2.2e-6。第2對數據檢驗后t為t2=-0.468 16，由于|t2|<2.262 2，故不能拒絕原假設，即可認為濾波算法估計樹木年齡均值與實際年齡均值之間無顯著差異，此時檢驗的p為0.640 7。

4 討論

通過試驗可以看出，針刺儀平穩卡爾曼濾波器理論可以應用于華北落葉松，使得針刺儀成為活立木年齡微損測定的有效工具，且具有以下特點：

1)針刺儀平穩卡爾曼濾波算法是微損的。傳統樹木年齡的確定方法是查數樹木年輪數，需獲取樹木生長芯或截取樹木橫截面圓盤(Jeongetal.,2017；Seoetal.,2017)，而針刺儀平穩卡爾曼濾波算法數據是基于針刺儀獲取抗鉆阻力值，由于針刺后在樹木上留下的洞口直徑小于1.5 mm，對樹木損害很小，因此與傳統樹木年代學方法確定樹木年齡相比，針刺儀平穩卡爾曼濾波算法是微損或準無損的。相比儀器測定分析獲取樹木年齡的方法，如圖像處理法、X-射線法和同位素法等，針刺儀平穩卡爾曼濾波算法使用針刺儀獲取數據，由于針刺儀體積小、質量輕，屬于便攜式設備，克服了以往儀器測定法笨重、成本高、操作繁瑣等缺點，便于野外進行試驗。

2)針刺儀平穩卡爾曼濾波算法提高了針刺儀測定樹木年齡的精度。針刺儀平穩卡爾曼濾波算法通過MATLAB編程實現，將抗鉆阻力值作為針刺儀平穩卡爾曼濾波器的輸入，根據活立木胸徑選擇適當的參數Rat值，就可以得出樹木的估計年齡，且與實際年齡誤差較小。針刺儀自身測定樹木年齡的平均相對誤差為-40.49%，利用針刺儀平穩卡爾曼濾波算法后，平均相對誤差減少至0.75%，同時將算法估計的樹木年齡與實際年齡進行成對數據t-檢驗，二者之間無明顯差異，說明算法應用于針刺儀抗鉆阻力值序列可以估計華北落葉松年齡，且精度較好。這樣就克服了針刺儀自身攜帶軟件判定樹木年齡準確率低的缺點，使針刺儀可以成為活立木年齡微損測定的有效工具。

3)針刺儀平穩卡爾曼濾波器理論的創新性。在1個濾波周期內，卡爾曼濾波器存在2個明顯的信息更新過程：時間更新過程和觀測更新過程(Stepanov，2011)，后驗估計值由先驗估計值、濾波增益矩陣和先驗協方差矩陣確定。針刺儀平穩卡爾曼濾波算法將針刺儀看作一個系統，給出其狀態空間模型，以一般的針刺儀卡爾曼濾波理論為基礎，推導出簡化針刺儀卡爾曼濾波表達式，使得抗鉆阻力最優估計值只由前一時刻的后驗估計值確定，迭代過程中可省略先驗估計值和先驗協方差矩陣計算，大大簡化了迭代過程，減少了計算工作量，使抗鉆阻力值卡爾曼濾波過程更加直觀、簡潔，克服了一般卡爾曼濾波參數較多、迭代過程復雜的缺點。

本研究根據后驗協方差存在上下確界的性質，進一步求出后驗協方差矩陣的極限值，將簡化針刺儀卡爾曼濾波表達式中的后驗協方差矩陣用其極限值表示，從而給出針刺儀平穩卡爾曼濾波器的定義。對針刺儀平穩卡爾曼濾波器的性質分析發現，其等價于一個離散差分方程，這樣就將抗鉆阻力的后驗估計值即最優估計值的求解過程轉換成差分方程的求解過程。由于差分方程的解只與參數s有關，而s僅與Rat有關，Rat為Q和R的比值，因此所求抗鉆阻力的最優估計值只與Rat有關。通過大量試驗發現，Rat值的選取與活立木胸徑有關，進而根據徑階大小給出了選擇參數依據。

卡爾曼濾波是處理噪聲的經典算法之一，本研究基于針刺儀獲取的抗鉆阻力數據，利用卡爾曼濾波思想，設定極小化性能指標為線性最小方差估計值，將針刺儀看作一個系統，建立系統的狀態方程和觀測方程，寫出一般卡爾曼濾波遞推表達式，并對其進行簡化得到簡化針刺儀卡爾曼濾波器，進一步對簡化針刺儀卡爾曼濾波器的表達式進行分析，求出后驗協方差的極限值，給出針刺儀平穩卡爾曼濾波器的定義。對于每株活立木，使用針刺儀獲取1組抗鉆阻力值序列，將其作為針刺儀平穩卡爾曼濾波器的輸入，輸出為抗鉆阻力的最優估計值序列，最優估計值序列的波峰波谷個數確定樹木年齡，結果僅僅與參數Rat的取值有關。通過研究樣本數據試驗結果顯示，總是可以找到合適的參數Rat值，使得估計年齡與實際年齡比較接近。

5 結論

本研究在理論上突破了原有卡爾曼濾波的自回歸迭代過程，推導出針刺儀平穩卡爾曼濾波器，使得一般針刺儀卡爾曼濾波器的迭代回歸過程得到簡化，且初始參數設定變得更加簡單，只由1個參數確定。通過給出的針刺儀平穩卡爾曼濾波器，可將針刺儀獲取的抗鉆阻力數據作為輸入，得出樹木年齡估計值，改進了針刺儀自帶DECOM軟件識別樹木年輪邊界準確率低、過于依賴人工經驗的缺點，可以作為一種微損估計活立木年齡的有效途徑。