應(yīng)用地基激光雷達(dá)對毛竹林有效葉面積指數(shù)的提取及敏感性分析1)

黃蘭鷹 鄭亞雄 張美曼 尹子旭 姜小雨 肖簫 官鳳英

(國家林業(yè)和草原局-北京市共建竹藤科學(xué)與技術(shù)重點(diǎn)實(shí)驗(yàn)室(國際竹藤中心)，北京，100102)

毛竹是我國南方重要的經(jīng)濟(jì)生態(tài)樹種，具有生態(tài)系統(tǒng)服務(wù)、生物多樣性保護(hù)和促進(jìn)社會經(jīng)濟(jì)發(fā)展等重要價(jià)值[1]。毛竹具有速生豐產(chǎn)、收獲周期短的特點(diǎn)[2]。毛竹與其他木材相比，竹材具有彈性強(qiáng)、韌性好等優(yōu)良力學(xué)性能，廣泛運(yùn)用于家具、膠板、竹碳和竹纖維等行業(yè)[3]。毛竹具有不同于其他森林類型的優(yōu)良固碳能力[4]，對毛竹的培育與經(jīng)營在緩解氣候變化中具有重要作用[5]。由于毛竹是典型的異齡結(jié)構(gòu)，擇伐仍然是毛竹林主要的采伐方式。但隨著我國經(jīng)濟(jì)快速的發(fā)展和人口老齡化的加劇，勞動力成本日益上升和勞動力資源缺乏，傳統(tǒng)竹林采伐方式已經(jīng)成為竹林經(jīng)營的瓶頸之一，也是限制竹產(chǎn)業(yè)發(fā)展的重要因素。為了降低竹材采伐成本，提高采伐效率，毛竹林帶狀采伐方式得到了極大的關(guān)注。目前對帶狀采伐后毛竹林恢復(fù)的研究主要集中在林下特征[6-8]，而對冠層的研究較少。冠層結(jié)構(gòu)控制著光照的質(zhì)量、數(shù)量、時(shí)間和分布，葉面積指數(shù)作為冠層結(jié)構(gòu)的參數(shù)之一，表征著植被生長狀況的關(guān)鍵物理量，控制植被的生物物理過程(如光合作用、呼吸作用、蒸騰作用、碳循環(huán)和降水?dāng)r截等[9-11])。同時(shí)，葉面積指數(shù)也是描述大氣和陸地生物地球化學(xué)交換過程的重要參數(shù)[12-14]，葉面積指數(shù)的大小不僅反映森林經(jīng)營管理的水平[15]，也可以為帶狀采伐后毛竹林恢復(fù)特征的評估提供數(shù)據(jù)支撐。

葉面積指數(shù)可以通過直接和間接兩種方法測量[16]。直接測量方法是通過破壞性采樣，對葉片進(jìn)行計(jì)數(shù)和測量，該方法可以得到高度準(zhǔn)確的結(jié)果，但缺點(diǎn)是具有破壞性、勞動強(qiáng)度大、耗時(shí)且適用范圍有限。針對葉面積指數(shù)直接測量方法的缺點(diǎn)，測量葉面積指數(shù)的間接方法應(yīng)運(yùn)而生，主要包括使用遙感影像(如MODIS、AVHRR等)和地面光學(xué)儀器(如LAI-2000植物冠層分析儀、半球攝影等)。然而，間接測量的方法也有其自身的局限性，遙感影像記錄葉面積指數(shù)較高的冠層存在光譜信號的飽和效應(yīng)，背景的光學(xué)特性在不同的空間和季節(jié)存在差異，以及對冠層下部結(jié)構(gòu)信息獲取的局限性，這些都是影響葉面積指數(shù)估算精度的不利因素[17]。光學(xué)儀器應(yīng)用貝爾定律，結(jié)合消光系數(shù)和孔隙率來估算植被的葉面積指數(shù)，其對曝光和光照射條件高度敏感，分辨率也是影響葉面積指數(shù)估測精度的重要因素[18]。

地基激光雷達(dá)(TLS)作為一種主動遙感技術(shù)，能夠直接、快速、精確地獲取研究對象的三維地理坐標(biāo)。因此，通過地基激光雷達(dá)提供的森林結(jié)構(gòu)的精確三維信息[19]，可以有效地評估樹冠特征、提高葉面積指數(shù)估算的準(zhǔn)確性，也避免了葉面積指數(shù)較高樹冠的飽和效應(yīng)。無論通過把三維TLS點(diǎn)云數(shù)據(jù)轉(zhuǎn)化為二維柵格數(shù)據(jù)估算葉面積指數(shù)[20-21]，或者利用點(diǎn)象限分析的方法計(jì)算葉面積指數(shù)(如形態(tài)變化[22]、局部表面擬合[23]、3D凸包算法[24]等)，這些方法無法重建樹木和進(jìn)行枝葉分離，也無法分析噪聲對葉面積指數(shù)估算的影響。因此，利用地基激光雷達(dá)信息計(jì)算葉面積指數(shù)，需要計(jì)算激光束與冠層的接觸率或間隙率，也需要對TLS點(diǎn)云進(jìn)行體元化[25]。體元化要求對原始點(diǎn)云數(shù)據(jù)進(jìn)行歸一化處理并為空間中的體元提供確定的體積，體元化的方法具有不需要對冠層組分的空間分布、大小或形狀進(jìn)行假設(shè)的特點(diǎn)。Hosoi et al.[25]利用TLS點(diǎn)云數(shù)據(jù)，根據(jù)體元化的方法在單木水平上計(jì)算了葉面積指數(shù)，認(rèn)為當(dāng)激光入射天頂角為57.5°時(shí)，葉傾角校正因子與葉傾角無關(guān)，可以近似等于常數(shù)；Zhu et al.[26]通過建立三維體元模型結(jié)合分層的方法計(jì)算葉面積指數(shù)，發(fā)現(xiàn)估算的葉面積指數(shù)比地面實(shí)測值偏低0.306，但是各分層空間變化曲線比較吻合；王洪蜀等[27]基于體元的冠層分析法對單棵玉蘭樹冠層的有效葉面積指數(shù)進(jìn)行反演，相對誤差為1.26%；Wu et al.[28]采用體元化的方法量化了澳洲芒果、堅(jiān)果和鱷梨的樹冠葉面積、葉面積密度和葉面積密度垂直分布的變化，其結(jié)果可以作為評價(jià)樹木健康的一般指標(biāo)，也可以幫助種植者改善灌溉和施肥決策。因此，基于體元估算葉面積指數(shù)是一種十分有效的方法。

然而，在計(jì)算葉面積指數(shù)時(shí)，需要仔細(xì)選擇體元大小，體元大小決定計(jì)算的效率和精度，影響葉面積指數(shù)估算的準(zhǔn)確性；噪音點(diǎn)也是影響葉面積指數(shù)估算準(zhǔn)確性的重要影響因子，目前很少有研究分析去噪對葉面積指數(shù)估算的影響；使用體元模型計(jì)算葉面積指數(shù)時(shí)，需要實(shí)際測量冠層葉片的傾斜角度，計(jì)算葉傾角的校正因子。一般假設(shè)冠層葉片隨機(jī)分布，將葉傾角的校正因子設(shè)為1.1，得到的葉面積指數(shù)是有效葉面積指數(shù)(LAIe)，而非真實(shí)葉面積指數(shù)。同時(shí)，目前使用體元方法估算有效葉面積指數(shù)多集中在喬木，應(yīng)用于毛竹方面的研究鮮見報(bào)道。因此，針對帶狀采伐毛竹林監(jiān)測和恢復(fù)特征評估的需要，以3種不同寬度(5、8、12 m)帶狀采伐毛竹林為研究對象，利用TLS多站掃描獲取毛竹林完整點(diǎn)云數(shù)據(jù)，采用基于體元的冠層分析方法對樣地毛竹林有效葉面積指數(shù)進(jìn)行提取，并進(jìn)一步分析有效葉面積指數(shù)的估測值對去噪和體元大小的敏感性。

1 研究區(qū)概況

研究區(qū)位于江蘇省宜興市宜興國有林場，該地區(qū)屬天目山余脈，南部為丘陵地區(qū)，北部為平原區(qū)。地處亞熱帶季風(fēng)氣候，降水豐沛，但降水量的時(shí)空分布不均勻，春夏降水比較集中，年降水量1 177 mm以上，年平均氣溫在17.5 ℃，年平均無霜期240 d。研究區(qū)毛竹林占80%，森林覆蓋率為97.5%。

2 研究方法

2.1 樣地設(shè)置

實(shí)驗(yàn)使用的地基激光雷達(dá)設(shè)備是Trimble TX8，該儀器具有360°×317°的視場角和每秒百萬點(diǎn)的數(shù)據(jù)獲取速度，在水平和垂直方向上均以0.036°為步長，Trimble TX8具體性能指標(biāo)見表1。樣地設(shè)置為5 m×20 m(P1)，8 m×20 m(P2)，12 m×20 m(P3)大小的3塊樣方，樣地地勢均較為平坦，且林下植被較少。數(shù)據(jù)采集時(shí)間為2019年12月，每個(gè)樣地掃描5站，采用“角落設(shè)置”的方法布站[28]，在幾乎無風(fēng)的條件下進(jìn)行掃描。在每塊樣地中布設(shè)了6個(gè)參考球，所有參考球不在同一個(gè)水平面上，保證掃描站每次至少能看到3個(gè)以上參考球，以便后期多站數(shù)據(jù)拼接處理。采用LAI-2000植被冠層分析儀在每塊樣地中心以及東、南、西、北4個(gè)方向共測量5次，葉面積指數(shù)平均值作為樣地的葉面積指數(shù)實(shí)測值[29]，以盡可能準(zhǔn)確表達(dá)樣地的葉面積指數(shù)，3塊樣地測得的葉面積指數(shù)值分別為2.78、2.94、2.54。

2.2 數(shù)據(jù)預(yù)處理

點(diǎn)云配準(zhǔn)：采用多站掃描的方式最大程度地減少竹枝和竹葉導(dǎo)致的遮擋效應(yīng)，盡可能獲取樣方完整的點(diǎn)云數(shù)據(jù)。利用Trimble TX8自帶的Trimble RealWorks 10.1軟件提供的自動識別和配準(zhǔn)功能，對樣地多站點(diǎn)云進(jìn)行配準(zhǔn)拼接，3塊樣地配準(zhǔn)的平均距離誤差均小于0.002 m，滿足林業(yè)應(yīng)用精度需求。

表1 Trimble TX8的性能指標(biāo)

點(diǎn)云去噪：點(diǎn)云數(shù)據(jù)中噪聲點(diǎn)的數(shù)量與掃描距離、掃描入射角、掃描速度和天氣狀況等多種因素有關(guān)。采用統(tǒng)計(jì)離群值分析的方法去除噪聲點(diǎn)[30]。該方法對每一個(gè)點(diǎn)搜索指定鄰域點(diǎn)個(gè)數(shù)的相鄰點(diǎn)，并計(jì)算該點(diǎn)到相鄰點(diǎn)的平均值距離、平均值距離的中值和標(biāo)準(zhǔn)差。如果這個(gè)點(diǎn)的平均值距離大于最大距離(最大距離=中值+標(biāo)準(zhǔn)差倍數(shù)×標(biāo)準(zhǔn)差)，則認(rèn)為是噪聲點(diǎn)，將被去掉。但該方法需要優(yōu)化相鄰點(diǎn)的數(shù)量(n)和標(biāo)準(zhǔn)差的倍數(shù)(K)兩個(gè)參數(shù)，通過對參數(shù)n∈[0,100]間隔為10、參數(shù)K∈[0,10]間隔為1進(jìn)行試驗(yàn)，當(dāng)參數(shù)n和K分別設(shè)置為50和5時(shí)效果最好。

2.3 葉面積指數(shù)提取

地形校正：為了去除地形對葉面積指數(shù)提取的影響，需要對樣方的激光雷達(dá)點(diǎn)云數(shù)據(jù)進(jìn)行歸一化處理，消除地形因子。首先采用改進(jìn)的漸進(jìn)加密三角網(wǎng)濾波算法從激光雷達(dá)點(diǎn)云數(shù)據(jù)中分離地面點(diǎn)[31]，然后利用反距離權(quán)重插值生成樣方的高程，最后對預(yù)處理的TLS數(shù)據(jù)進(jìn)行高程歸一化處理。

體元化：在樣地點(diǎn)云數(shù)據(jù)地形校正后，使用體元的冠層分析法計(jì)算有效葉面積指數(shù)。該方法是根據(jù)冠層與激光束的接觸頻率，并對點(diǎn)云數(shù)據(jù)進(jìn)行體元化[25]的基礎(chǔ)上計(jì)算有效葉面積指數(shù)。體元大小根據(jù)激光儀的測距范圍、光束的發(fā)散度和直徑來決定[32]。假設(shè)物體完全攔截了發(fā)射的激光束，那么所有的激光束都會提供有效的返回值，則可以通過計(jì)算兩個(gè)相鄰激光束之間的距離[33]，具體計(jì)算公式為：

dCD=(tanα×Di)-2R。

式中：R為半徑，Di為掃描儀的掃描范圍，BDv為光束發(fā)散度，BDm為激光儀出口處的光束直徑，dCD為兩個(gè)相鄰激光束之間的距離，α為角分辨率。

Trimble TX8在垂直和水平平面上以0.036°的固定步長發(fā)射激光束，每個(gè)激光束發(fā)散度為0.3 mrad，激光儀出口處的光束直徑為6 mm，最大的掃描范圍為120 m，則相鄰激光束之間的距離(最大未采樣距離)為33.40 mm。

體元大小確定后，將整個(gè)區(qū)域的點(diǎn)云劃分為l×w×h個(gè)體元，體元大小為Δl×Δw×Δh(即為dCD)，計(jì)算公式如下：

式中：int是取整符，表示將計(jì)算得到的小數(shù)取整為與其最接近的整數(shù)，(xmin,ymin,zmin)表示TLS點(diǎn)云數(shù)據(jù)(x,y,z)的最小值，(xmax,ymax,zmax)表示TLS點(diǎn)云數(shù)據(jù)(x,y,z)的最大值。

Δl×Δw×Δh是體元大小，l×w×h表示整個(gè)點(diǎn)云區(qū)域的體元數(shù)量。體元的體素值通過判斷體元內(nèi)是否有激光點(diǎn)來確定，如果體元內(nèi)激光點(diǎn)個(gè)數(shù)大于等于1，代表激光束被攔截，體素值賦為1，否則體素值賦為0。使用體元化的方法，可以將點(diǎn)云數(shù)據(jù)分組到每個(gè)體元，因此，體元化減少了點(diǎn)云數(shù)量，并提高了點(diǎn)云接觸頻率的計(jì)算效率，有利于信息的表達(dá)。

有效葉面積指數(shù)的計(jì)算：在點(diǎn)云數(shù)據(jù)體元化后，將其分為s層，每層厚度為Δs，將每個(gè)水平層厚度設(shè)置為1 m，即Δs=1，則s=h/Δs(h為樣地內(nèi)毛竹最高的竹高)。根據(jù)構(gòu)建的三維網(wǎng)格，統(tǒng)計(jì)每個(gè)水平層內(nèi)三維格網(wǎng)的總數(shù)和包含激光點(diǎn)的格網(wǎng)數(shù)量，計(jì)算每個(gè)水平層激光束的接觸頻率[25]，計(jì)算公式如下：

N(s)=nI(s)/nT(s)。

式中：N(s)為第s水平層激光點(diǎn)的接觸頻率，nI(s)為第s水平層包含激光點(diǎn)的格網(wǎng)數(shù)量(即體素值為1的體元數(shù)量)，nT(s)為第s水平層的三維格網(wǎng)總數(shù)。

根據(jù)第s水平層激光點(diǎn)的接觸頻率計(jì)算第s水平層的有效葉面積指數(shù)L(s)，計(jì)算公式如下：

L(s)=α(θ)×N(s)。

式中：α(θ)為葉片傾斜度校正因子，本文假設(shè)毛竹冠層葉片隨機(jī)分布，將葉傾角校因子設(shè)為1.1[25]。將各層的LAIe進(jìn)行累加，得到整個(gè)統(tǒng)計(jì)區(qū)域內(nèi)的有效葉面積指數(shù)，計(jì)算公式如下：

采用常用的絕對誤差和相對誤差對估測的有效葉面積指數(shù)進(jìn)行精度評價(jià)。

3 結(jié)果與分析

3.1 去噪對有效葉面積指數(shù)估算的影響

圖1顯示了單株毛竹去噪前后的TLS點(diǎn)云數(shù)據(jù)，大多數(shù)的噪聲點(diǎn)位于相鄰葉片之間的空間和冠層上部，通過統(tǒng)計(jì)離群值可以很好的去除這些噪聲。由表2可知，去噪前后有效葉面積指數(shù)有明顯的差異，去噪前，3塊樣地估測的有效葉面積指數(shù)分別是3.62、3.75、3.26，與地面實(shí)測值相比，3塊樣地有效葉面積指數(shù)估測的精度分別為69.78%、72.45%、71.65%。去噪前后樣地有效葉面積指數(shù)與實(shí)測葉面積指數(shù)的絕對誤差從0.84減少到了0.34，且3塊樣地有效葉面積指數(shù)估測的精度均達(dá)到了86%以上，其中P1樣地的精度提高了16.55%，P2和P3樣地的精度分別提高了15.59%、14.57%。

圖1 去噪前后點(diǎn)云的比較

表2 去噪對有效葉面積指數(shù)的估算

3.2 體元大小對有效葉面積指數(shù)估算的影響

體元大小會影響冠層結(jié)構(gòu)信息的細(xì)節(jié)以及接觸率的計(jì)算精度，也是獲取植被冠層信息的重要參數(shù)。為了探究體元大小對有效葉面積指數(shù)的影響，本文使用了6種體元大小(1、2、3、4、5、6 cm)構(gòu)建體元化模型，得到不同體元大小情況下有效葉面積指數(shù)的變化。

由圖2可知，當(dāng)體元大小從1 cm增加到6 cm時(shí)，3塊樣地的有效葉面積指數(shù)分別從1.18、0.98、0.72持續(xù)增加到5.08、5.80、6.00。且體元大小和有效葉面積指數(shù)之間呈高度對數(shù)相關(guān)，3塊樣地的決定系數(shù)分別為0.98、0.97以及0.94。這表明應(yīng)用體元冠層分析法估算的有效葉面積指數(shù)對體元大小十分敏感，體元大小的直接影響有效葉面積指數(shù)估算的精度。因此，選擇合適的體元大小構(gòu)建三維體元模型對準(zhǔn)確反演樣地有效葉面積指數(shù)具有重要意義。

與柵格圖像的像元類似，在構(gòu)建三維體元模型時(shí)，體元越小，體元數(shù)量越多，數(shù)據(jù)量越大，越能詳細(xì)表達(dá)冠層內(nèi)部的結(jié)構(gòu)信息，進(jìn)而能更加準(zhǔn)確模擬激光束穿過竹林冠層的內(nèi)部情況。但是當(dāng)體元大小小于激光采樣最小間距時(shí)，此時(shí)屬性為1的值將達(dá)到最大值，而屬性為0的值會隨著體元減小而繼續(xù)增加，導(dǎo)致計(jì)算結(jié)果偏小；體元越大，體元總數(shù)量越少，簡化了點(diǎn)云數(shù)據(jù)，體元越大包含更多的空隙，表達(dá)的冠層結(jié)構(gòu)更加模糊，導(dǎo)致結(jié)果高估。因此，應(yīng)用體元冠層分析法計(jì)算葉面積指數(shù)時(shí)，需要在體元大小與間隙數(shù)量之間進(jìn)行平衡，以提高葉面積指數(shù)計(jì)算精度。本文研究中，有效葉面積指數(shù)在一個(gè)數(shù)量級上變化，與地面實(shí)測葉面積指數(shù)相比，合適的體元大小為3 cm，這與利用激光儀的測距范圍、光束的發(fā)散度和直徑所推導(dǎo)出來的體元大小相差不大。

圖2 體元大小與有效葉面積指數(shù)的關(guān)系

4 討論

根據(jù)本文的實(shí)驗(yàn)結(jié)果可知，與實(shí)際葉面積指數(shù)測量值相比，應(yīng)用TLS數(shù)據(jù)估算的有效葉面積指數(shù)存在一定的誤差，這與前人的研究一致[27，34]。一是由于竹枝和密集的竹葉引起的遮擋效應(yīng)，激光束不能把所有的葉子記錄下來，遮擋是應(yīng)用TLS數(shù)據(jù)估算有效葉面積指數(shù)的主要限制因素[25]。二是3塊樣地假設(shè)的葉傾角校正因子均為1.1，沒有考慮毛竹真實(shí)葉傾角的差異；同時(shí)，點(diǎn)云的濾波方法、樣地冠層真實(shí)的邊界等因素也影響有效葉面積指數(shù)的估算精度。三是未直接對非光合組織(如莖干和竹枝)進(jìn)行分離，林分中非光合組織的點(diǎn)云參與了有效葉面積指數(shù)計(jì)算，影響了有效葉面積指數(shù)估算的準(zhǔn)確性，未來應(yīng)結(jié)合竹子的莖葉分離算法，進(jìn)一步探討竹稈和竹枝等非光合組織對有效葉面積指數(shù)估算的影響。

應(yīng)用點(diǎn)云的最大距離，采用統(tǒng)計(jì)離群值分析的方法去除噪聲點(diǎn)，并評估了噪聲點(diǎn)對有效葉面積指數(shù)估算精度的影響。點(diǎn)云去噪后，3塊樣地葉有效葉面積指數(shù)的估測精度分別提高了16.55%、15.59%、14.57%，表明噪聲點(diǎn)是影響有效葉面積指數(shù)估算精度的重要因素。地基激光雷達(dá)較高的發(fā)射頻率和較短的范圍掃描將激光脈沖入射到鄰近區(qū)域，導(dǎo)致在葉片之間檢測到更多的反射能量[35]。因此，需要在掃描精度、掃描距離和噪聲點(diǎn)數(shù)量之間達(dá)到一個(gè)平衡，以提高有效葉面積指數(shù)估算的準(zhǔn)確性。

體元大小是估算有效葉面積指數(shù)的重要影響因素。應(yīng)用體元的冠層分析方法是通過計(jì)算激光束與冠層的接觸頻率進(jìn)行，有效葉面積指數(shù)與體元大小的相關(guān)關(guān)系可以通過體元大小與間隙數(shù)量之間的平衡解釋。把三維格網(wǎng)中有激光點(diǎn)的體元賦值為1，體元越大，賦值為1的體元越多，計(jì)算的有效葉面積指數(shù)也大。但是賦值為1的體元中將包含更多潛在的空隙，造成結(jié)果存在高估，因此，在點(diǎn)云密度足夠高的條件下，體元大小應(yīng)該越小越好。實(shí)際上，由于脈沖之間的步長較小，屬于同一單個(gè)葉子的兩點(diǎn)之間仍然存在一定的間距，當(dāng)體元太小時(shí)，該間距很容易被視為空白體元，同時(shí)體元過小會增加數(shù)據(jù)量，這違背體元化“降維”的思想。體元大小通過影響間隙率進(jìn)而影響有效葉面積指數(shù)估算的準(zhǔn)確性，因此，計(jì)算有效葉面積指數(shù)時(shí)，需要選擇合適的體元大小。Wu et al.[15]研究發(fā)現(xiàn)，盡管體元大小對葉面積指數(shù)估算的精度有影響，但不同體元大小計(jì)算出的葉面積指數(shù)變化是一致的，表明如果關(guān)注相對葉面積指數(shù)的變化，則體元大小對葉面積指數(shù)估算的精度影響較小。

5 結(jié)論

利用TLS多站掃描對3種不同寬度(5、8、12 m)帶狀采伐毛竹林進(jìn)行了數(shù)據(jù)獲取，采用基于體元的冠層分析方法對樣地毛竹林有效葉面積指數(shù)進(jìn)行提取，并進(jìn)一步分析了有效葉面積指數(shù)的估測值對去噪和體元大小的敏感性。(1)體元化的TLS點(diǎn)云數(shù)據(jù)能夠?qū)崿F(xiàn)毛竹林有效葉面積指數(shù)的提取，3塊樣地有效葉面積指數(shù)的提取精度分別是86.33%、88.04%、86.22%；(2)去噪后，3塊樣地有效葉面積指數(shù)估算的精度分別提高了16.55%、15.59%、14.57%；(3)有效葉面積指數(shù)的估算精度與體元大小高度相關(guān)，通過與實(shí)測葉面積指數(shù)相比，適合有效葉面積指數(shù)估算精度的體元大小為0.03 m。本文的研究結(jié)果可為帶狀采伐毛竹林監(jiān)測和恢復(fù)特征評估提供數(shù)據(jù)支撐，對利用TLS數(shù)據(jù)提高有效葉面積指數(shù)的估算精度具有參考作用。