姚建峰,符利勇,宋新宇,王雪峰,趙燕東,鄭一力,業(yè)巧林,翟勝楠
(1. 信陽師范學(xué)院計算機與信息技術(shù)學(xué)院,信陽 464000; 2. 中國林業(yè)科學(xué)研究院資源信息研究所,北京 100091;3. 信陽師范學(xué)院數(shù)學(xué)與統(tǒng)計學(xué)院,信陽 464000; 4. 北京林業(yè)大學(xué)工學(xué)院,北京 100083;5. 南京林業(yè)大學(xué)信息科學(xué)技術(shù)學(xué)院,南京 210037)
樹木生長除受遺傳特性影響外,還受氣候、環(huán)境和森林經(jīng)營等綜合外界因子的影響[1]。生長輪是樹干橫截面上的同心圓環(huán),通常每年形成一輪,故又稱為年輪。通過建立樹木年輪寬度、年輪密度等年輪參數(shù)與氣候因子之間的數(shù)學(xué)模型,可重建過去的氣候環(huán)境變化[2-3]。年輪寬度反映了1年內(nèi)外界因子對樹木生長的綜合影響,僅通過年輪寬度指標(biāo)所獲取的樹木年輪信息有限。而年輪密度包括每個年輪的平均密度、晚材密度、早材密度、晚材最大密度、早材最小密度等指標(biāo),可反映出樹木生長季前期和后期外界綜合因子對樹木生長的影響[4];同時,年輪密度是表征木材材質(zhì)綜合性變異的重要指標(biāo),通過研究年輪密度動態(tài)變異規(guī)律,可為森林經(jīng)理、木材定向培育、木材加工與利用等提供理論依據(jù)[5-6]。
從20世紀40年代末開始,研究人員就開始嘗試各種不同的年輪密度測量方法,其中包括體積比重法、光透視法、硬度測定法、抗張強度測定法、β微粒法等,但這些方法對樹木年代學(xué)研究而言都存在著不同程度的限制條件,以至于不能被廣泛地應(yīng)用于年輪學(xué)研究中[7]。隨著X射線技術(shù)成功應(yīng)用于年輪密度的測量中,年輪密度的研究才逐漸興起[8]。目前,X射線法仍是年輪密度測量中使用最廣泛的一種測量方法[8]。但X射線法儀器設(shè)備昂貴,維護成本高,體積大、質(zhì)量大,操作過程煩瑣,對所測量樣本要求較高,不適合野外作業(yè),因此科研人員一直試圖尋找一種可攜帶、操作簡單的年輪密度測量方法來代替X射線法[9]。微鉆阻力儀是使用電機控制鉆針勻速鉆入樹木并實時記錄鉆針阻力的一種儀器[10-11],鉆針阻力與木材密度正相關(guān)。當(dāng)鉆針鉆入晚材部分時,鉆針阻力較大;當(dāng)鉆針鉆入早材部分時,鉆針阻力較小;當(dāng)鉆針沿徑向鉆入樹木時,鉆針阻力呈峰谷交替規(guī)律變化。因此,根據(jù)鉆針阻力曲線圖可獲取年輪寬度和年輪密度信息[12]。
目前,微鉆阻力儀主要由德國Rinntech和IML(Instrumenta Mechanik Labor)公司生產(chǎn)。由于鉆針的振動,這2個廠家生產(chǎn)的微鉆阻力儀鉆針阻力信號中都含有部分的噪聲信號,鉆針阻力曲線圖中的波峰不能與樹木年輪的晚材一一對應(yīng),年輪識別主要通過專家人工完成,年輪識別精度主要依賴專家經(jīng)驗。2017年,在唐守正院士的指導(dǎo)下,中國林業(yè)科學(xué)研究院資源信息研究所和北京林業(yè)大學(xué)工學(xué)院成立了微鉆阻力儀研制團隊,試制了3個微鉆阻力儀樣機[13-16],研究了鉆針阻力表達方法[17],提出了4種樹木年輪識別算法[18-21]。筆者使用該團隊最新研制的微鉆阻力儀和德國Resistograph 650-S微鉆阻力儀測量輻射松(Pinusradiata)木塊的早材密度和晚材密度,以期測試微鉆阻力法測量樹木生長輪早晚材密度的精度。
微鉆阻力儀由2個電機控制:一個是直流電機,控制鉆針旋轉(zhuǎn)速度;另一個是步進電機,控制鉆針進給速度。微鉆阻力儀的機械傳動結(jié)構(gòu)如圖1所示。鉆針與直流電機軸相連,鉆針旋轉(zhuǎn)速度與直流電機轉(zhuǎn)速相同,直流電機安裝在絲桿滑塊上。絲桿滑塊在傳動絲桿的帶動下沿導(dǎo)軌直線移動,步進電機軸與傳動絲桿軸相連,傳動絲桿的轉(zhuǎn)速與步進電機轉(zhuǎn)速相同。當(dāng)步進電機旋轉(zhuǎn)時,步進電機帶動傳動絲桿同步旋轉(zhuǎn),使絲桿滑塊在導(dǎo)軌上直線移動,從而帶動鉆針移動[16]。
1.步進電機;2.聯(lián)軸器;3.后絲桿支撐座;4.直流電機;5.鉆針夾;6.傳動絲桿;7.直線導(dǎo)軌;8.鉆針;9.鉆針套頭;10.前絲桿支撐座;11.絲桿滑塊。圖1 機械傳動結(jié)構(gòu)Fig. 1 Mechanical transmission structure
當(dāng)鉆針鉆入晚材部分時,晚材密度較大,鉆針阻力較大,直流電機轉(zhuǎn)速下降,直流電機控制器輸出電壓上升,使直流電機電流增大,輸出力矩增加,直流電機轉(zhuǎn)速上升;當(dāng)鉆針鉆入早材部分時,早材密度較小,鉆針阻力較小,直流電機轉(zhuǎn)速上升,直流電機控制器輸出電壓下降,使直流電機電流減小,輸出力矩下降,直流電機轉(zhuǎn)速降低。由于直流電機的輸出力矩不宜實時測量,且直流電機電流信號中存在大量噪聲信號,因此,使用直流電機的電壓、測定轉(zhuǎn)速和設(shè)定轉(zhuǎn)速來間接表示鉆針阻力[17]。鉆針阻力(F)計算公式為:
F=Un0/n
(1)
式中:U為直流電機電壓,V;n0為直流電機設(shè)定轉(zhuǎn)速,r/min;n為直流電機實際轉(zhuǎn)速,r/min。
試驗儀器主要有自制的微鉆阻力儀和德國Rinntech公司生產(chǎn)的Resistograph 650-S微鉆阻力儀。自制微鉆阻力儀的鉆針旋轉(zhuǎn)速度n0設(shè)置為3 500 r/min,鉆針進給速度設(shè)置為15 cm/min,鉆針阻力采樣間距為0.002 5 mm,最大測量長度為500 mm。Resistograph 650-S微鉆阻力儀鉆針進給速度是60 cm/min,鉆針阻力采樣間距為0.01 mm,最大測量長度為550 mm。該儀器使用直流電機的功率間接表示鉆針阻力[10],但是阻力單位不是功率單位W,而是該公司自定義的單位Resi。因此,仍不能確定該儀器鉆針阻力具體的表示方法。
2.2.1 樣品加工和處理方法
為了能用體積法測量年輪的早材密度和晚材密度,本試驗材料是從木材市場上購買的年輪寬度較寬的新西蘭輻射松(Pinusradiata)木材(20 cm×20 cm×200 cm)。新西蘭輻射松生長快,心材部分年輪寬度可達2 cm,年輪線清晰,木材強度適中,易于分割年輪的早材部分和晚材部分,是用體積法測量早、晚材密度的理想材料。在原木材兩端截取厚度為3.5 cm的板材,每個板材加工成5個寬度為2 cm的長方體木塊,共10個長方體木塊。長方體的長沿木材徑向方向,長方體的寬盡量與大部分年輪線平行,長方體的高與樹高方向相同。首先根據(jù)年輪線弧度在板材橫截面上標(biāo)記出每個木塊的位置,標(biāo)記方法如圖2所示;然后使用手鋸沿標(biāo)記線鋸成長方體木塊,木材長度為15~18 cm,并使用打磨機修整、打磨長方體木塊,使木塊高3 cm、寬2 cm;最后把長方體木塊置于烘箱中烘烤至絕干狀態(tài)。為了防止長方體木塊受熱變形,先把烤箱溫度設(shè)置為25 ℃,然后每間隔6 h左右增加20 ℃,最后把烘箱溫度設(shè)置為105 ℃,恒溫干燥48 h,使木塊干燥至絕干狀態(tài),關(guān)閉烘箱,使木塊自然冷卻至常溫狀態(tài)。
圖2 木塊加工方法Fig. 2 Method of making wood samples
2.2.2 鉆針阻力取樣方法
分別使用2個微鉆阻力儀沿木塊徑向鉆入,為了減小因鉆針鉆入年輪晚材方向的不同對鉆針阻力的影響,每個儀器取2次鉆針阻力,一次取樣方向是從邊材鉆入心材,另一次取樣方向是從心材鉆入邊材。為了防止鉆針路徑重合,每2個取樣點在水平和垂直方向上相距1 cm左右。記錄鉆針阻力序號、鉆針路徑長度和鉆針鉆入木塊方向。
2.2.3 年輪早材密度和晚材密度測量方法
6.dabusu qami arun-a bui ɡebel ,dalai tenɡkis ba na arun-a要問鹽是哪里來的,是從大海和湖泊里來的)
使用游標(biāo)卡尺測量長方體木塊長ai、寬b、高h,精確到0.00 1 cm,使用電子天平測量長方體木塊的絕干質(zhì)量mi,精確到0.01 g。然后使用工具刀從邊材方向開始依次切削每個年輪的早材或者晚材,每切削一次,使用砂紙打磨長方體木塊切削處,使切削面平整光滑,最后再次測量打磨后的長方體木塊長度ai+1和質(zhì)量mi+1,并根據(jù)式(2)計算被切除早材或晚材部分的木材絕干密度ρ。
(2)
式中:mi和mi+1分別為第i次切削木材前后的木塊質(zhì)量,g;ai和ai+1分別為第i次切削木材前后的木塊長度,cm;b為木塊寬度,cm;h為木塊高度,cm。
2.2.4 微鉆阻力儀早晚材平均阻力計算方法
根據(jù)每次切除的早材或晚材的木材寬度,計算該部分對應(yīng)的鉆針阻力平均值。
1)平滑去噪處理。由于自制微鉆阻力儀和德國Resistograph 650-S原始阻力信號中存在高頻噪聲信號,為了減小噪聲信號對年輪密度測量精度的影響,首先對鉆針原始數(shù)據(jù)進行平滑處理,平滑寬度均設(shè)置為2 mm。由于自制微鉆阻力儀鉆針阻力采樣間距為0.002 5 mm,因此平滑窗口寬度設(shè)置為800個數(shù)據(jù)采樣點;Resistograph 650-S鉆針阻力采樣間距是0.01 mm,因此平滑窗口寬度設(shè)置為200個數(shù)據(jù)采樣點。德國Resistograph 650-S平滑處理前和平滑處理后的阻力曲線圖與木塊年輪的對比分析見圖3。從圖3中可以看出,原始阻力數(shù)據(jù)經(jīng)過平滑處理后可較清晰地反映出樹木年輪信息,峰谷差值較大的波峰與年輪的晚材位置相對應(yīng)。
2)去鉆針針桿摩擦力。盡管微鉆阻力儀鉆針阻力主要集中在鉆針針頭上,但是鉆針針桿與鉆孔之間仍存在一定的摩擦力,且該摩擦力隨著鉆孔深度的增加而增加。鉆針鉆入木材之前其針頭處于空載狀態(tài),鉆針阻力主要是微鉆阻力儀機電系統(tǒng)空載時的摩擦力;鉆針鉆出木材之后鉆針針頭也處于空載狀態(tài),鉆針阻力主要包括針桿阻力和微鉆阻力儀機電系統(tǒng)空載時的摩擦力。為了減小針桿摩擦力對木材年輪密度測量精度的影響,采用以下步驟去針桿摩擦力:①計算鉆針鉆入木材之前5 mm微鉆阻力儀的鉆針平均阻力f0;②計算鉆針鉆出木材之后5 mm微鉆阻力儀的鉆針平均阻力f1;③按式(3)計算鉆針鉆入木材深度為d時去針桿摩擦力的鉆針阻力F1。
F1=F0-d(f1-f0)/l
(3)
式中:F0為原始鉆針阻力;d為鉆針鉆入木材的深度;l為鉆針路徑長度。
圖4 2個方向阻力曲線及其平均阻力曲線對比Fig. 4 Comparison between resistance curves in two directions and its average resistance curve
2.2.5 數(shù)據(jù)統(tǒng)計分析、建模與測試
首先對測試數(shù)據(jù)進行異常數(shù)據(jù)處理和統(tǒng)計分析,刪除測試數(shù)據(jù)中的異常數(shù)據(jù),并分析早晚材絕干密度,以及自制微鉆阻力儀和Resistograph 650-S所測早晚材阻力數(shù)據(jù)的最大值、最小值、平均值和標(biāo)準差,然后隨機抽樣4/5的試驗數(shù)據(jù)作為建模數(shù)據(jù)集,以每個早材或晚材的平均阻力為自變量,以每個早材或晚材的絕干密度為因變量,使用中國林業(yè)科學(xué)研究院資源信息研究所自主開發(fā)的ForStat統(tǒng)計軟件建立每個儀器早晚材鉆針平均阻力與木材絕干密度的線性模型。最后把剩余1/5的試驗數(shù)據(jù)作為測試數(shù)據(jù)集,使用式(4)分別計算2個微鉆阻力儀預(yù)測年輪密度的預(yù)測精度ξ。
(4)
由于原始測試數(shù)據(jù)中可能存在異常數(shù)據(jù),為了減小異常數(shù)據(jù)對數(shù)據(jù)分析和數(shù)據(jù)建模的影響,首先使用箱線圖法查找原始測試數(shù)據(jù)中的異常數(shù)據(jù),然后使用刪除法處理這些異常數(shù)據(jù)。10個長方體木塊被分割成133個早材木片和137個晚材木片(有4個木塊兩端均為晚材),因此,共測得133組早材數(shù)據(jù)和137組晚材數(shù)據(jù)。分別對早材木片和晚材木片的絕干密度、自制微鉆阻力儀鉆針阻力和Resistograph 650-S鉆針阻力進行歸一化處理,把箱線圖的延伸倍數(shù)設(shè)置為1.5,使用R語言分別繪制早材木片和晚材木片歸一化數(shù)據(jù)的箱線圖,如圖5所示。
圖5 箱線圖對比Fig. 5 Comparison by boxplot
從圖5中可以看出:①早材絕干密度的最大值、自制微鉆阻力儀早材木片阻力的最大值和次大值為異常數(shù)據(jù);②晚材絕干密度的最大值、自制微鉆阻力儀晚材木片阻力的最大值和次大值為異常數(shù)據(jù)。從原始測試數(shù)據(jù)中查找對應(yīng)的異常數(shù)據(jù),然后刪除包含異常數(shù)據(jù)的測試記錄。
對刪除異常數(shù)據(jù)后的測試數(shù)據(jù)進行統(tǒng)計分析,正常早材木片數(shù)據(jù)有130組,正常晚材數(shù)據(jù)有134組,分析各種數(shù)據(jù)的最大值、最小值、平均值和標(biāo)準差,統(tǒng)計分析結(jié)果如表1所示。
表1 試驗數(shù)據(jù)Table 1 Test data
從表1可以看出:①新西蘭輻射松所有晚材木片的平均密度是556.839 kg/m3,所有早材木片的平均密度是386.850 kg/m3,所有晚材平均密度是所有早材平均密度的1.439倍;②自制微鉆阻力儀所有晚材部分的平均阻力為4.301 V,所有早材部分的平均阻力是3.706 V,所有晚材部分平均阻力是所有早材部分的1.161倍;③Resistograph 650-S所有晚材部分平均阻力為228.351 Resi,所有早材部分平均阻力是175.577 Resi,所有晚材部分平均阻力是所有早材部分的1.301倍;④單片早材絕干密度的最大值是472.116 kg/m3,單片晚材絕干密度的最小值是377.235 kg/m3。因此,新西蘭輻射松晚材密度與早材密度差異較大,2個微鉆阻力儀所有晚材的平均阻力均高于所有早材的平均阻力,但由于有少量單片早材的密度高于少量單片晚材的密度,所以有少量單片早材的平均阻力高于少量單片晚材的平均阻力(單個年輪早材部分或晚材部分的平均阻力)。因此,僅使用鉆針阻力值不能確定木材的早晚材分界線。
隨機抽取正常試驗數(shù)據(jù)中104組早材木片數(shù)據(jù)和107組晚材木片數(shù)據(jù)作為建模數(shù)據(jù)集,分別建立自制微鉆阻力儀和德國Resistograph 650-S微鉆阻力儀所測早晚材鉆針平均阻力與木材絕干密度線性模型,2個線性模型的擬合曲線如圖6所示,模型參數(shù)如表2所示。
自制微鉆阻力儀所測早晚材平均阻力與木材絕干密度的線性相關(guān)系數(shù)為0.657,Resistograph 650-S所測早晚材平均阻力與木材絕干密度的線性相關(guān)系數(shù)為0.655。從擬合結(jié)果可以看出:早晚材鉆針阻力與木材絕干密度呈顯著正相關(guān);自制微鉆阻力儀所測早晚材平均阻力與木材絕干密度的線性相關(guān)系數(shù)高于Resistograph 650-S。
圖6 早晚材鉆針平均阻力與木材絕干密度線性模型擬合曲線Fig. 6 Linear fitting curves between average drill resistance and absolute dry wood density of earlywood and latewood
表2 早晚材平均阻力與絕干密度線性模型Table 2 Linear model between average resistance and absolute dry wood density of earlywood and latewood
把剩余的54個早晚材木片數(shù)據(jù)作為測試數(shù)據(jù)集,分別使用自制微鉆阻力儀和德國Resistograph 650-S微鉆阻力儀鉆針阻力與木材絕干密度之間的線性模型預(yù)測木片的絕干密度,2個儀器的測量精度如表3所示。
表3 測量精度Table 3 Measurement accuracies 單位:%
從表3可以看出,2個微鉆阻力儀測量早晚材木材絕干密度的平均測量精度高達86%。為了檢驗2個儀器預(yù)測早晚材絕干密度與實際測量絕干密度,以及2個儀器之間預(yù)測的木材絕干密度是否存在顯著性差異,做了2個儀器預(yù)測早晚材絕干密度與實際測量絕干密度間的t檢驗,檢驗結(jié)果如表4所示。
表4 t檢驗結(jié)果Table 4 Results of t test
從表4可以看出,2個儀器預(yù)測早晚材絕干密度與實際測量絕干密度、自制微鉆阻力儀和Resistograph 650-S預(yù)測的絕干密度在顯著性水平0.1上均無顯著性差異。
從測試結(jié)果可以看出,2個微鉆阻力儀測量的樹木年輪早晚材絕干密度與實際測量絕干密度在顯著性水平0.1上無顯著性差異,2個儀器測量的早晚材絕干密度平均精度高達86%,說明微鉆阻力法測量樹木年輪的早材絕干密度和晚材絕干密度是可行的。但有少量年輪密度測量試驗數(shù)據(jù)精度較低,這主要是由以下幾個原因造成的。①鉆針振動:微鉆阻力儀鉆針較細長,當(dāng)鉆針高速旋轉(zhuǎn)時,鉆針振動幅度較大,鉆針阻力信號中存在部分噪聲信號。②鉆針針桿摩擦力不穩(wěn)定:鉆針針桿變形、彎曲程度受鉆針針頭進給方向上阻力的影響。當(dāng)鉆針鉆入晚材部分時,木材密度增加,鉆針在進給方向上的阻力增加,鉆針針桿變形、彎曲程度增加,針桿阻力增加;當(dāng)鉆針鉆入早材部分時,木材密度減小,鉆針在進給方向上的阻力減小,鉆針針桿變形、彎曲程度減小,針桿阻力減小。③樹木年輪晚材的柱面不規(guī)則,當(dāng)鉆針鉆入晚材部分時,鉆針針頭與晚材的接觸面積有時增大,有時減小。
對比分析2個儀器的年輪密度測量結(jié)果可知,自制微鉆阻力儀測得的早晚材平均阻力與木材絕干密度的線性相關(guān)系數(shù)和早晚材密度的平均測量精度均高于Resistograph 650-S,說明自制微鉆阻力儀機械結(jié)構(gòu)設(shè)計合理,鉆針阻力表示方式可行。
微鉆阻力儀鉆針阻力與木材密度正相關(guān)[10,22-26],當(dāng)鉆針沿徑向鉆入樹木時,鉆針阻力呈峰谷規(guī)律變化,通過鉆針阻力可以反映出樹木年輪寬度和年輪密度等信息。文獻[27]中使用Resistograph微鉆阻力儀測量了30余株山楊的年輪寬度,使用COFECHA對獲取的年輪寬度進行交叉定年和控制檢驗,并使用ARSTAN分別建立山楊的標(biāo)準化年表、差值年表、自回歸年表,研究發(fā)現(xiàn)山楊年輪寬度指標(biāo)與該地區(qū)月平均氣溫顯著相關(guān)(P<0.05)。文獻[10]中分別使用Resistograph微鉆阻力儀和X射線密度儀測量了冷杉、落葉松、云杉、松樹、椴樹和楊樹等6種樹木的阻力數(shù)據(jù)和密度數(shù)據(jù),盡管這2種方式在測量路徑、木材含水率、取樣精度等方面可能有所相同,但研究發(fā)現(xiàn)鉆針阻力曲線圖與X射線微密度曲線圖的變化趨勢高度一致。由于鉆針阻力信號中含有部分噪聲信號,通過鉆針阻力曲線圖不易確定早晚材分界線,因此,微鉆阻力法在測量年輪密度方面尚未得到廣泛應(yīng)用。
筆者最初試圖尋找早材密度與晚材密度的分界值,并根據(jù)早晚材密度分界值來確定早晚材之間的鉆針阻力分界值,但研究發(fā)現(xiàn)同一樹木不同位置的木材密度差異較大,部分早材密度高于部分晚材密度,因此樹木早材密度與晚材密度之間可能尚無統(tǒng)一的分界值。通過對比分析鉆針阻力曲線圖和木塊年輪圖,發(fā)現(xiàn)鉆針阻力曲線中峰谷差值較大的波峰與樹木年輪的晚材位置相對應(yīng),說明通過鉆針阻力可測量年輪密度信息。最后通過人工方式切割早晚材并建立早晚材絕干密度與早晚材鉆針阻力平均值之間的線性模型,研究發(fā)現(xiàn)2個微鉆阻力儀的線性模型R2均在0.4以上,平均早晚材密度測量精度高達86%,微鉆阻力法測量結(jié)果與體積法測量結(jié)果在顯著性水平0.1上無顯著差異,研究結(jié)果表明微鉆阻力法測量年輪密度可行。但是從圖6中可以看出,鉆針阻力與木材絕干密度的散點圖比較分散,相關(guān)性仍不高,說明鉆針阻力還可能受其他因素的影響,如樹脂含量、外部推力、儀器晃動、鉆針進給速度、鉆針旋轉(zhuǎn)速度、鉆針鋒利程度等[28]。當(dāng)使用微鉆阻力儀測量活力木年輪密度時,還需要考慮木材含水率對鉆針阻力的影響。
下一步將研究影響鉆針阻力的因素、含水率與鉆針阻力之間的關(guān)系、早晚材分界線的確定方法等,以期直接通過微鉆阻力測量活力木年輪密度,為林業(yè)、森林經(jīng)理學(xué)、氣候?qū)W等學(xué)科研究提供一種精確、微損、快速、簡便的年輪測量儀器。
1)自制微鉆阻力儀測得的早晚材平均阻力與木材絕干密度的線性相關(guān)系數(shù)為0.657,Resistograph 650-S測得的早晚材平均阻力與木材絕干密度的線性相關(guān)系數(shù)為0.655。
2)2個微鉆阻力儀測量早晚材的絕干密度平均精度高達86%,微鉆阻力法測量結(jié)果與體積法測量結(jié)果在顯著性水平0.1上無顯著差異。
3)自制微鉆阻力儀測得的早晚材平均阻力與木材絕干密度的線性相關(guān)系數(shù)和早晚材密度的平均測量精度均高于Resistograph 650-S,說明自主研制的微鉆阻力儀機械結(jié)構(gòu)設(shè)計合理,鉆針阻力表示方式可行。