謝 杰 ，邢艷秋 ，尤號田 ，田 昕 ，安立華，姚松濤

（1. 東北林業(yè)大學 森林作業(yè)與環(huán)境研究中心，黑龍江 哈爾濱 150040；2. 中國林業(yè)科學研究院 資源信息研究所， 北京 100091）

森林通過光合作用、呼吸作用和蒸騰作用實現(xiàn)了物質(zhì)能量的交換[1]，而葉片是光合作用的主要載體，是生物生產(chǎn)力轉(zhuǎn)換的關鍵[2-3]。葉面積指數(shù)（Leaf Area Index，Lai）用來表征葉片的疏密程度[4]，是反映植被生長狀況的重要指標[5]，被定義為單位地表面積上單面植物光合作用面積的總和[6-8]。

激光雷達（Light Detection and Ranging，LiDAR）作為一種先進的主動遙感技術，既能彌補光學遙感在反演植被結構參數(shù)上的不足[9]，又能夠穿透植被冠層，獲取地面的三維結構信息及植被的垂直結構信息[10]，現(xiàn)已廣泛運用于林業(yè)研究[11]。近年來，國內(nèi)外學者利用LiDAR數(shù)據(jù)對森林Lai進行了大量的研究[12-13]。如：李文娟等[14]基于機載LiDAR點云數(shù)據(jù)，利用樣方內(nèi)地面點數(shù)量與樣方內(nèi)總點云數(shù)量的比值，求得激光穿透指數(shù)（Laser Penetration Index，Lpi）用于估測森林Lai，估測精度R2=0.730。Luo等[15]利用機載LiDAR強度數(shù)據(jù)，提取樣方Lpi用于估測森林Lai，結果表明基于點云強度Lpi的估測精度R2為0.825，RMSE為0.165。通過對先前類似研究進行分析發(fā)現(xiàn)，由LiDAR數(shù)據(jù)計算得到的Lpi能夠有效估測森林Lai，但在計算過程中常采用地面點云強度總和與樣方點云強度總和之比，未曾考慮多回波強度數(shù)據(jù)間的差異。因為激光脈沖在穿透森林冠層時，會造成部分能量損失，因而多次回波強度數(shù)據(jù)間的可比性相對很小。

本研究以內(nèi)蒙古依根農(nóng)場機載LiDAR點云數(shù)據(jù)為基礎研究數(shù)據(jù)，通過對點云數(shù)據(jù)進行去噪、分類、強度校正及高程歸一化等處理，從中依據(jù)回波類型的不同分別提取了點云數(shù)量和點云強度的激光穿透指數(shù)，之后利用提取的激光穿透指數(shù)借助多元線性回歸算法對森林Lai進行估測，以期實現(xiàn)森林Lai的高精度估測。

1 研究區(qū)概況與數(shù)據(jù)獲取

1.1 研究區(qū)概況

研究區(qū)位于內(nèi)蒙古呼倫貝爾市西北部、額爾古納市東南部的上庫力農(nóng)場（120°36′50.48″～120°52′56.53″E， 50°21′11.08″～ 50°24′32.00″N），東與牙克石市為鄰，東南與陳巴爾虎旗相連，西和西南與拉布大林農(nóng)牧場接壤，北與三河馬場交界，如圖1所示。地處較復雜的山岳地形，屬于寒溫帶大陸性季風氣候，海拔在600～700 m。山脈丘陵陰坡廣泛分布著以白樺Betula platyphylla為主的天然次生林，其中混生樹種主要有落葉松Larix gmelinii和樟子松Pinus sylvestris等，林下灌層主要由石棒繡線菊Spiraea media、筐柳Salix linearistipularis等組成。

圖1 研究區(qū)位置Fig. 1 Location of study area

1.2 機載LiDAR數(shù)據(jù)獲取

2012年8月26日由飛機搭載Leica ALS60機載激光雷達掃描系統(tǒng)，相對飛行高度1 300 m，平均速度為160 km/h，當天天氣狀況良好。激光脈沖波長為1 064 nm，最大掃描頻率200 kHz，掃描方式為線性掃描，平均點云密度在2～4個脈沖點/m2，光斑大小0.22 m，共獲取東西方向22個條帶的數(shù)據(jù)。數(shù)據(jù)采用標準LiDAR存儲格式las1.2，該數(shù)據(jù)格式存儲了點云的坐標值、高程值、回波強度、回波次數(shù)、自定義分類等信息。

1.3 野外數(shù)據(jù)獲取

為了避免野外數(shù)據(jù)采集時偶然因素造成的影響以保證模型的穩(wěn)定性，于2012年8月在研究區(qū)中部連續(xù)布設了80塊大小為10 m×10 m的方形樣地，樣地內(nèi)分布樹種均為白樺。在研究區(qū)開闊地帶架設基站，采用型號為Juno SB的手持式GPS對樣地的4個角點和中心點進行定位，之后利用基站數(shù)據(jù)對手持GPS坐標進行差分運算，使最后差分定位精度在厘米級。與此同時，利用Lai-2000冠層分析儀分別對每塊樣地的中心點及角點進行測量，取多次測量的平均值作為樣地Lai的實測值，為了避免林下灌木、雜草的影響，野外測量Lai時，Lai-2000冠層分析儀始終處于距離地面1.3 m的位置。野外Lai采集統(tǒng)計結果如表1所示。

表1 野外Lai數(shù)據(jù)基本統(tǒng)計量Table 1 Basic statistics of measured Lai data

2 研究方法

2.1 LiDAR數(shù)據(jù)去噪

激光雷達飛行時，收集的數(shù)據(jù)不僅包含所需要的有用數(shù)據(jù)，同時也會產(chǎn)生噪聲點，這部分噪聲點通常為孤立點，往往是高于地物點或者是低于地面點。本研究在TerraSolid軟件平臺上，利用孤立點算法，以某一點為球心，給定相應的閾值半徑，若此球內(nèi)無其他點，即認為此點為噪聲點，最后人工檢查修正。

2.2 LiDAR點云能量強度校正

由于大氣傳輸過程中的衰減、目標物與發(fā)射器之間的距離等原因會產(chǎn)生誤差，為減小誤差累積，在研究前對點云能量進行近似校正，利用Baltsavias等[16]研究激光雷達能量方程進行近似校正。對于本研究采用同一系統(tǒng)進行數(shù)據(jù)采集，飛機水平飛行，數(shù)據(jù)獲取時天氣晴朗，近似認為硬件參數(shù)不變，只考慮目標的面積、反射率以及距離對接收能量的影響，傳感器接收的能量與距離的四次方成反比[17-18]。為了排除距離對計算結果產(chǎn)生影響，按照簡化公式（1）進行強度校正。

式中，Rs是參考距離（m），I是回波強度值（w），I(Rs)是校正到參考距離處的強度值（w），R是發(fā)射器與目標之間的距離（m），Hf、Hp分別為飛行絕對航高和回波點的高程（m），α是瞬時掃描角（°）。

2.3 LiDAR數(shù)據(jù)分類

采用不規(guī)則三角網(wǎng)濾波算法進行點云分類，插值算法生成數(shù)字高程模型，利用所有點云高程值與數(shù)字高程模型做差，得到歸一化后的點云。然后對歸一化后的數(shù)據(jù)進行進一步的劃分，將高于1.3 m的點記為植被點，低于1.3 m的點記為地面點。分別統(tǒng)計樣方內(nèi)的所有點云、植被點、地面點。在研究區(qū)中，大部分為白樺林，白樺生長呈圓錐形，堆疊枝葉較少，通透性較好，使得LiDAR數(shù)據(jù)多為首次回波（包括單次回波）和末次回波點，中間回波數(shù)量極少，因此本研究將點云按照首次回波（包括單次回波）、末次回波進行統(tǒng)計，同時統(tǒng)計各類點云的能量值。

2.4 基于點云數(shù)量的參數(shù)提取

本研究根據(jù)點云的回波類型，對激光穿透指數(shù)進行了細分，得到基于點云數(shù)量的多回波次數(shù)的激光穿透指數(shù)nLPIfirst、nLPIlast、nLPIcan。其中nLPIfirst為首次回波數(shù)量激光穿透指數(shù)，nLPIlast為末次回波數(shù)量激光穿透指數(shù)，nLPIcan為冠層回波數(shù)量激光穿透指數(shù)，計算公式見式3～5所示：

式中，nFg為樣方內(nèi)首次回波（包括單次回波）地面點的數(shù)量，nFc為樣方內(nèi)首次回波點（包括單次回波）冠層點的數(shù)量，nLg為樣方內(nèi)末次回波地面點的數(shù)量。

2.5 基于點云能量的參數(shù)提取

冠層孔隙的大小與激光穿透指數(shù)有著直接的關系，經(jīng)過校正后的點云強度值，對于冠層和地面，其強度值與物體表面積成正比，見公式（6），點云強度的值能夠較好的反映出冠層孔隙的大小。

式中I表示強度值，A表示面積，ρ為物體反射率。

設地表反射率為ρg，植被層反射率為ρc，因此植被反射率比值系數(shù)來代替真實地表反射率，即由于地面單次回波強度和冠層單次回波強度值有一定的波動范圍，但仍然在一個定值附近變化，通過計算兩種單次回波強度值的最大概率分布時對應的強度值的比值，其計算公式如公式（7）：

本研究區(qū)與邢艷秋等[19]研究區(qū)相同，且源數(shù)據(jù)也相同，根據(jù)邢艷秋等[19]研究結果，確定研究區(qū)內(nèi)α值為3.0。

當激光點云出現(xiàn)首次回波和末次回波時，首次回波點和末次回波點的的能量強度值由反射面積決定，面積不同，回波的強度值不同?；邳c云數(shù)量的Lpi不能夠很好的將由于面積不同而引起的差異反映出來。而激光在穿透冠層時，會有一部分的能量損耗，多次回波中不同點云之間的能量可比性很小。因此本研究基于多回波類型提取點云能量的激光穿透指數(shù)iLPIfirst、iLPIlast、iLPIcan，其中iLPIfirst為首次回波數(shù)量激光穿透指數(shù)，iLPIlast為末次回波數(shù)量激光穿透指數(shù)，iLPIcan為冠層回波數(shù)量激光穿透指數(shù)，計算公式詳見式8～10所示：

式中，iFg為樣方內(nèi)首次回波（包括單次回波）地面點的強度值，iFc為樣方內(nèi)首次回波（包括單次回波）冠層點的強度值，iLg為樣方內(nèi)末次回波地面點的強度值，α為植被反射率比值系數(shù)。

iLPIfirst排除了末次回波對激光穿透指數(shù)的影響，LiDAR首次回波包含的信息較多，將直接反映出首次回波對Lai的影響。為了增加對冠層小孔隙的靈敏度，iLPIlast的地面回波點采用首次回波地面點和末次回波地面點之和。iLPIcan反映冠層的信息，利用末次回波地面點和首次回波植被點對冠層穿透指數(shù)進行計算。

2.6 LAI模型建立及精度評價

本研究擬對提取的6個變量建立模型。在野外采集的80個樣方中隨機選取50個用于建模，余下30個用于預測模型精度評價。假定激光穿透指數(shù)與Lai具有線性關系，模型公式如式（11）：

式中：k為斜率；e為誤差項；Lpi為激光穿透指數(shù)。

本研究在單變量回歸模型的基礎上，嘗試采用多變量回歸模型，多變量回歸模型往往存在變量共線性問題，本研究采用方差擴大因子（Variance In fl ation Factor，Vif）來評價自變量之間的共線性問題，見公式（12）。方差擴大因子（Vif）描述了自變量之間共線性的嚴重程度，值越大表明共線性越嚴重，一般以小于10為判斷依據(jù)，即當Vif＜10時，表明自變量之間不存在共線性問題，反之當Vif＞10則存在嚴重的共線性問題[12]。

式中：Vifj為自變量Xj的方差擴大因子，R2j為自變量Xj對其余p-1個自變量的復決定系數(shù)。

預測模型精度評價主要用決定系數(shù)R2，其值越大，說明建模精度越好。預測模型精度采用平均絕對偏差（Mean Absolute Deviation，Mad）來評價，其值越小，說明預測精度越高[20]。

3 結果與分析

3.1 LiDAR點云處理結果

激光點云預處理后，對點云進行分類及高程歸一化處理，得到樣方原始點云，根據(jù)回波次數(shù)得到樣方首次回波地面點、末次回波地面點、首次回波植被點、末次回波植被點，如圖2所示：

圖2 激光雷達點云數(shù)據(jù)Fig.2 Airborne LiDAR point cloud data

從圖1中可以看出，激光雷達能夠穿透森林冠層到達地面，能夠很好地反映出森林的垂直結構信息。地面點被分為首次回波地面點和末次回波地面點，首次回波地面點和末次回波地面點各占有一定比例，可以間接反映出森林冠層情況。首次回波植被點明顯多于末次回波植被點，這對于激光束來說末次回波點多落在地面上，只有少量較密集冠層能夠阻擋激光束末次回波點，這也與樹種有關，在研究區(qū)中，大部分為白樺林，白樺生長呈圓錐形，堆疊枝葉較少，通透性較好，使得末次回波植被點相對較少。

3.2 單變量模型估測結果

對6個變量分別建立模型，擬合結果如圖3所示。

圖3 單變量模型建模精度Fig.3 Univariate empirical model estimation results

從圖3中可以看出，實測Lai與各激光穿透指數(shù)之間大多為負相關關系，R2在0.004～0.836之間。對于變量nLPIcan、iLPIcan而言，決定系數(shù)R2低，兩者與Lai之間并無明顯的線性關系，而其余4個變量均表現(xiàn)出一定的線性關系。Lai隨著nLPIfirst、nLPIlast、iLPIfirst、iLPIlast的增加而減小，這種規(guī)律符合實際情況，Lai越小，激光穿透率越大。對剩下的30個樣方數(shù)據(jù)進行預測，nLPIfirst、nLPIlast、iLPIfirst、iLPIlast模型的預測精度見圖4。

通過表2可以看出，模型精度最高的為iLPIfirst模型，其R2達到了0.836，其模型預測精度也最高，Mad=0.091。iLPIlast、iLPIfirst、nLPIlast反演精度次之，預測精度也隨之降低。綜合單變量模型結果，單變量iLPIfirst模型最好，即：Lai=7.911-7.327iLPIfirst。

圖4 單變量參數(shù)模型驗證樣方預測值與實測值Fig.4 Univariate parameter model veri fi cation of sample prediction and measured value

表2 單變量模型反演精度評價Table 2 Accuracy evaluation of single variable empirical model inversion

3.3 多變量模型估測結果

為了進一步的挖掘變量之間的信息，本研究嘗試利用多變量模型來反映與Lai的關系，多變量回歸模型能夠包含更多有解釋力的變量，避免遺漏變量偏差問題。為了解決解釋變量共線性的問題，將共線性變量排除，采取逐步回歸法得到多變量回歸模型，剔除嚴重共線性的變量組合，最終得到兩組符合條件的結果，如表3所示。

表3 多元回歸分析變量多重共線性檢驗Table 3 Multivariate linear regression analysis of variables

從表2中可以看出，剔除線性相關較大變量后，得到兩組多元回歸變量，iLPIfirst，nLPIcan和iLPIfirst，nLPIcan，iLPIcan，對兩組變量進行擬合，得到的擬合結果如表3所示，預測精度如圖5所示。

對比表2和表4可以看出，多變量回歸模型的精度都高于單變量回歸模型，R2均大于0.85，Mad都低于0.09。相比于單變量回歸模型，多變量回歸模型能夠更好的反映森林葉面積指數(shù)，得到了較好的精度。對比單變量iLPIfirst模型和兩個多變量模型。單變量iLPIfirst模型利用iLPIfirst參數(shù)很好的估計了森林葉面積指數(shù)，而多變量iLPIfirst、nLPIcan模型，相比于單變量iLPIfirst模型，加入了nLPIcan參數(shù)，反演精度有提高。單看nLPIcan單變量與Lai的關系，他們沒有顯著相關性，說明變量nLPIcan包含了森林Lai的相關信息，但是與Lai不是簡單的線性關系。變量iLPIcan的加入，使得多變量iLPIfirst、nLPIcan、iLPIcan模型相比于多變量iLPIfirst、nLPIcan模型精度提高了，iLPIcan包含許多森林Lai信息，能夠很好地和iLPIfirst共同反演森林Lai。

表4 多元回歸分析變量模型及精度評價Table 4 Multivariate regression analysis model and precision evaluation

4 結論與討論

圖5 多變量參數(shù)模型驗證樣方預測值與實測值Fig. 5 Multivariate parameter model veri fi cation of sample prediction and measured value

目前，利用遙感技術估測森林葉面積指數(shù)已經(jīng)得到廣泛運用，小光斑激光雷達作為一種主動遙感技術，為森林葉面積指數(shù)的估測提供新的方法。本研究通過對機載激光雷達點云數(shù)據(jù)進行去噪、強度校正、高程歸一化、多回波類型分類等處理，充分利用激光雷達點云數(shù)據(jù)中多回波類型之間所含信息的差異，提取基于多回波類型的激光穿透指數(shù)，應用這些參數(shù)對森林Lai進行反演，得到以下結果：

1）單變量模型中基于點云能量首次回波的iLPIfirst反演的精度最高，其次是基于點云能量末次回波的iLPIlast，二者均高于基于點云數(shù)量提取的Lpi的反演結果。結果表明：相比于點云數(shù)量，LiDAR點云能量信息包含了更多有關森林葉面積指數(shù)的信息，能夠更好地反演森林葉面積指數(shù)。此結論與駱社周等[21]研究結果類似，駱社周等[21]利用校正能量提取Lpi用于反演森林Lai，結果估測精度（R2=0.825）高于點云數(shù)量提取Lpi反演森林Lai的估測精度（R2=0.803）。

2）反演模型中多變量模型反演精度要高于單變量模型反演精度。多變量模型中，利用iLPIfirst、nLPIcan、iLPIcan3個變量的反演模型精度最高，R2為0.883。邢艷秋等[19]利用經(jīng)過強度校正后的LiDAR點云，通過對樣方內(nèi)點云進行分束，得到樣方內(nèi)所有單束Lpi后求其均值LPImean，用于反演森林Lai，結果R2為0.80，低于本研究所用多變量模型的估測精度但高于部分單精度模型的精度。分析其主要原因，本研究進一步細分了首次回波、末次回波、冠層回波的Lpi，得到了多變量回歸模型，彌補了單變量解釋能力不足的缺陷。多變量模型估測結果從另一個方面說明單變量參數(shù)存在一定的局限性，若要得到高精確度的森林葉面積指數(shù)，可以聯(lián)合應用多個變量參數(shù)以消除單變量解釋能力的不足，進而實現(xiàn)變量間的優(yōu)勢互補。

本研究利用機載LiDAR離散點云數(shù)據(jù)多回波間的能量信息，估測了白樺林Lai，獲得了較好的結果。但由于本研究所用研究區(qū)樹種結構較單一，所得模型對其它樹種的適用性結果未知，因此在未來的研究中應逐漸加入其他樹種以驗證模型對其他樹種的適用性。

[1]Carlson TN, Ripley DA. On the relation between NDVI,fractional vegetation cover, and leaf area index [J].Remote Sensing of Environment. 1998,62(3):241-252.

[2]Chen JM, Rich PM, Gower ST,et al.Leaf area index of boreal forests: Theory, techniques, and measurements[J].Journal of Geophysical Research Atmospheres, 1997,102(D24):29429-29443.

[3]凌成星, 鞠洪波, 張懷清,等. 基于植被指數(shù)比較的濕地區(qū)域LAI遙感估算研究[J].中南林業(yè)科技大學學報, 2016, 36(5):11-18.

[4]孫 華, 羅朝沁, 林 輝,等. 基于k-NN算法的葉面積指數(shù)遙感反演[J].中南林業(yè)科技大學學報,2016,36(12):11-17.

[5]Bréda NJJ. Ground-based measurements of leaf area index: a review of methods, instruments and current controversies[J].Journal of Experimental Botany. 2003,54(392):2403-2417.

[6]Jonckheere I, Fleck S, Nackaerts K,et al. Review of methods for in situ leaf area index determination. I. Theories, sensors and hemispherical photography[J].Agricultural & Forest Meteorology, 2004,121(1):19-35.

[7]王希群,馬履一,賈忠奎,等.葉面積指數(shù)的研究和應用進展[J].生態(tài)學雜志,2005,24(5):537-541.

[8]Watson DJ. Comparative Physiological Studies on the Growth of Field Crops: I. Variation in Net Assimilation Rate and Leaf Area between Species and Varieties, and within and between Years[J].Annals of Botany, 1947,11(1):41-76.

[9]Wulder M. Optical remote-sensing techniques for the assessment of forest inventory and biophysical parameters[J].Progress in Physical Geography, 1998,22(4):449-476.

[10]Myneni RB, Maggion S, Iaquinta J,et al. Optical remote sensing of vegetation: Modeling, caveats, and algorithms [J].Remote Sensing of Environment, 1995,51(51):169-188.

[11]Lim K, Treitz P, Wulder M,et al. LIDAR remote sensing of forest structure[J].Progress in Physical Geography, 2003,27(1):88-106.[12]Peduzzi A, Wynne RH, Fox TR,et al. Estimating leaf area index in intensively managed pine plantations using airborne laser scanner data[J].Forest Ecology & Management, 2012,270(4):54-65.

[13]Solberg S. Mapping gap fraction, LAI and defoliation using various ALS penetration variables[J].International Journal of Remote Sensing, 2010,31(5):1227-1244.

[14]李文娟,趙傳燕,別 強,等.基于機載激光雷達數(shù)據(jù)的森林結構參數(shù)反演[J].遙感技術與應用, 2015,30(5):917-924.

[15]Luo SZ, Wang C, Zhang GB,et al.Forest Leaf Area Index(LAI) Estimation Using Airborne Discrete-Return Lidar Data[J].Chinese Journal of Geophysics, 2013,56(3):233-242.

[16]Baltsavias EP. Airborne laser scanning: basic relations and formulas[J].Isprs Journal of Photogrammetry & Remote Sensing,1999, 54(2-3):199-214.

[17]賴旭東.機載激光雷達基礎原理與應用[M].北京: 電子工業(yè)出版社,2010.

[18]Lerma JL. A new methodology to estimate the discrete-return point density on airborne lidar surveys[J].International Journal of Remote Sensing, 2014,35(4):1496-1510.

[19]邢艷秋,霍 達,尤號田,等.基于機載LiDAR單束激光穿透指數(shù)的白樺林LAI估測[J].應用生態(tài)學報, 2016, 27(11):3469-3478

[20]Kvalseth TO. Cautionary Note about R2[J].The American Statisti cian,1985,39(4):279-285.

[21]駱社周,王 成,張貴賓,等.機載激光雷達森林葉面積指數(shù)反演研究[J].地球物理學報, 2013,56(5):1467-1475.