多源遙感數據協同的景洪橡膠林葉面積指數光飽和特征研究

羅洪斌 舒清態 王 強 王冬玲 字 李 謝福明

( 西南林業大學林學院，云南 昆明 650233)

葉面積指數（LAI）定義為單位地面面積上所有葉片單面面積的總和[1]，是表征植物冠層結構生理的一個重要參數，目前對其測定的方法大致可分為直接測量法和間接測量法兩大類[2]。直接測量借助一定的測量工具對獲取的樣本葉片直接測量其面積。雖然結果真實可靠，但會對植被造成一定的破壞性，測量效率低，測量范圍小。間接測量通常以測定植被冠層內外輻射的參數為基礎，然后進行植被冠層結構特征的推斷，借助一定的數學模型推算葉面積指數；另一種是直接使用葉面積指數測量儀得到葉面積指數。以上測量方法僅能獲取地面有限點的LAI值，不能滿足大面積LAI監測。遙感技術的發展使LAI無損、大面積的監測成為了可能[3]，但仍需通過獲取的地面的LAI樣本點數據，并與遙感影像中提取的光譜植被指數之間建立回歸統計模型，從而預測植被LAI的空間分布。

目前較多使用冠層分析儀來獲取地面葉面積指數，但只能以點進行樣地尺度的葉面積指數獲取，且測量值比真實值要偏小。機載激光雷達（LiDAR）能夠同時獲得較大范圍的植被冠層和垂直結構信息，目前已廣泛應用于郁閉度、樹高、葉面積指數等森林參數的提取[4]，但是機載激光雷達數據獲取成本高，且覆蓋區域有限。因此，LiDAR數據與光學遙感數據的多數據源協同就起到了不可代替的作用。然而，研究葉面積指數與光學遙感信號響應時發現，當葉面積指數增大到一定程度時，遙感信號響應值不再增大，即出現飽和現象。這是定量遙感估測中碰到的世界性難題，也是影響到葉面積指數估測精度的重要原因。武紅敢等[5]運用TM影像提取的植被指數監測馬尾松（Pinus massoniana）林分葉面積指數的動態變化時發現，當葉面積指數值大于6時，葉面積指數變化的靈敏度降低。Zhao等[6]研究發現，當葉面積指數超過3時，冠層反射率達到飽和。Baret等[7]根據模型推斷植被葉面積指數值大于3時，光譜信息不能精確反演葉面積指數。對比相關研究來看，針對遙感影像原始波段反射率與LAI的光飽和特征的研究較少，Landsat 8 OLI影像數據作為當前LAI遙感估測的重要信息源，對其原始波段光飽和特征研究具有重要意義。

本研究以云南省景洪市機載LiDAR飛行范圍為研究區，以陸地衛星影像Landsat 8 OLI為主要信息源，在利用機載LiDAR點云數據提取地面LAI的基礎上，借助地統計學中的半方差函數，對葉面積指數各波段光譜飽和特性進行了分析，并對光譜“飽和點”進行求算。研究可為LAI遙感估測中的光譜飽和現象研究提供參考。

1 研究區概況

研究區為云南省景洪市，位于東經100°41′35″～101°00′16″，北緯 21°59′26″～22°04′28″（圖1）。該區地處云貴高原、橫斷山系南端、瀾滄江大斷裂帶兩側，地勢北高南低，溝壑縱橫山地眾多，最高海拔2 196.8 m，最低海拔485 m。屬熱帶濕潤型季風氣候，干濕季分明，常年高溫酷暑，日溫差較大，年溫差小，基本無霜。平均氣溫在18.6～21.9 ℃，降水充沛，年平均降水量達1 200～1 700 mm，年平均日照時間為1 800～2 300 h。在熱帶氣候和地形的作用下，形成以赤紅壤、磚紅壤為主的土壤類型，且土層深厚，自然肥力高。森林植被隨海拔和地形起伏變化呈明顯的垂直分布，植被主要以熱帶季節雨林、熱帶山地雨林、南亞熱帶季風常綠闊葉林為主。橡膠（Hevea brasiliensis）林在景洪山地的各海拔梯度均有分布，其生長分布的海拔上線是1 200 m左右，是景洪地區極具代表性的喬木植被。

圖1 研究區位置Fig. 1 Location of the study area

2 材料與方法

2.1 數據來源與處理

2.1.1 機載激光雷達數據

機載激光雷達數據來源于中國林科院的LiCHy（LiDAR，CCD and Hyperspectral）系統[8]。數據于 2014年 4月由 RIEGL LMS-Q680i傳感器獲取，平均飛行高度為851 m，掃描角度為±30°，最大脈沖發射頻率400 kHz，點云密度為3.9個/m2。離散點云最終以LAS1.2格式存儲，飛行區域長32.1 km，寬3.7 km，6個航帶總面積118 km2。

數據的前期處理由數據提供商完成，后期處理運用北京數字綠土公司開發的LiDAR 360軟件對點云數據進行去噪、植被點和地面點的分離（0～2 m設為低矮植被點，2～6 m設為中等植被點，6～30 m設為高等植被點，30～150 m為建筑物點）、歸一化等處理。

2.1.2 Landsat 8 OLI影像數據

研究區所需的2景Landsat 8 OLI數據來源于地理空間數據云網站（http：//www.gscloud.cn/），拍攝時間為2014年4月，具體參數見表1。為滿足研究需要，本研究使用ENVI 5.3軟件對其進行了輻射定標和大氣校正、地形校正等處理。

表1 Landsat 8 OLI影像基本信息Table 1 Basic information of Landsat 8 OLI image

2.1.3 輔助數據

景洪市森林資源二類調查矢量小班圖層數據，利用ArcGIS對橡膠林分布區進行裁剪。

2.2 有效葉面積指數的提取

有效葉面積指數在LiDAR 360軟件下進行提取，利用地面點生成的DEM對點云進行歸一化處理，即將植被點的高程轉變為該點相對于地面的高度值。為了與Landsat 8 OLI柵格大小匹配，有效葉面積指數以30 m×30 m分辨率輸出。由于沒有考慮聚集情況，也沒有考慮光合作用部分和非光合作用部分的分離，LiDAR 360提取的葉面積指數指的是有效植被葉面積指數。葉面積基于間隙率模型進行計算，結合用戶設定的葉傾角分布計算得到有效葉面積指數（已有研究表明葉傾角的橢圓分布適用于大部分森林的實際情況，值為0.5[9]），其計算方法如下：

式中：a是平均掃描角，ai是第i個點的掃描角度；GF是間隙率，k是消光系數，消光系數與樹冠的葉傾角分布緊密相關；ng是提取的激光點云高度值低于高度閾值（2 m）的地面點數，n是總點數。

2.3 訓練樣本的選取

根據已有的研究區橡膠樹分布矢量圖裁剪得到條帶區橡膠林的有效葉面積指數。在ArcGIS軟件中創建漁網產生隨機樣點，并提取樣點所對應的有效葉面積指數。為保證提取樣本的代表性，本研究結合葉面積指數大小與光譜反射率散點圖對變異值較大的點進行剔除，最終剩下1 649個點作為先驗樣本數據（圖2），樣本最大值、最小值、平均值、標準差分別為9.36、0.01、2.50、1.83。

圖1 LiDAR條帶樣地點分布Fig. 2 Location distribution of LiDAR strips

2.4 橡膠林葉面積指數光飽和值的求算

在研究森林葉面積指數與遙感信號響應時發現，當葉面積指數增大到一定程度時遙感信號（包括光學和微波）的響應值不再增大，達到“飽和點”，即出現飽和現象。因此，本研究以葉面積指數為自變量，以特征波段反射率為因變量，借助地統計學中半方差函數理論求解不同波段的飽和值[10]。求解模型有球狀模型、高斯模型、指數模型和線性模型等，當前使用較普遍的是球狀模型，模型如下：

式中：Y（LAI）為光譜反射率值；C0為葉面積指數在LAI=0時的光譜反射率值；C是供高，即光譜變化率；C0+C為葉面積指數飽和時光譜反射率最大或最小值；BS即為各個波段的光譜飽和值。

若令b0=C0，b1=3C/2BS，b2=-C/2BS3則球狀模型變為一般的多元線性回歸模型：y=b0+b1×LAI+b2×LAI3，運用最小二乘法即可對模型參數進行求解，從而獲得不同波段所對應的葉面積指數飽和值。

3 結果與分析

3.1 相關性分析

使用SPSS 19對葉面積指數和不同波段反射率值進行皮爾遜相關性雙側檢驗，通過分析相關系數和顯著性水平，篩選出與葉面積指數相關性較好的變量，進行“飽和點”的求算，結果見表2。由表2可知，葉面積指數與Landsat 8 OLI原始波段反射率值均呈極顯著相關（P＜0.01），相關系數大小介于0.28～0.45之間，其中近紅外波段呈正相關，其余波段呈負相關。

表2 葉面積指數與Landsat 8 OLI影像各波段的相關性分析Table 2 Correlation analysis between leaf area index and each band of Landsat 8 OLI image

3.2 葉面積指數光譜響應分析

Landsat 8 OLI的7個波段反射率與葉面積指數的相關性分析結果表明，各個波段反射率與葉面積指數均呈極顯著相關（P＜0.01），相關系數均大于0.2。因此，本研究對7個波段光譜響應特征進行分析，Landsat 8 OLI的7個波段反射率與葉面積指數變化見圖3。

圖3 Landsat 8 OLI原始波段光譜反射率信息Fig. 3 Landsat 8 OLI original spectral reflectance information for each band

由圖3可知，B5的反射率變化隨著葉面積指數的增大而增大，進而達到飽和，其余波段的反射率均隨著葉面積指數的增加而減小，最后達到飽和。結合植被光譜特征分析可知，在可見光范圍內（B1、B2、B3、B4）葉綠素是影響植物光譜的重要因素，葉綠素大量吸收攝入的能量，因此隨著葉面積指數的增加光譜反射率逐漸減小；在近紅外波段（B5）隨著葉面積指數的增多，大量的葉片在近紅外波段產生較高的反射率，能量透過表層的葉片后，被下一層葉片反射，從而增強了表層葉片的反射能量；而在短波紅外（B6、B7）范圍內，植被的光譜響應表現為水汽的吸收，葉面積指數增大，葉片水分吸收越大。從光譜反射率變化散點圖看出，總體雖然呈現出飽和的趨勢，但是在LAI較小范圍有部分點并不呈現此變化規律。

3.3 識別葉面積指數光譜飽和點

使用SPSS 19軟件分別對半方差函數轉化成的二元一次方程進行參數求解，并對回歸系數進行t檢驗，擬合的回歸方程系數見表3。

表3 回歸系數的顯著性Table 3 Significance of regression coefficients

由表3可知，除了近紅外波段（B5）的回歸系數b2顯著性小于0.05以外，其余波段所有自變量回歸系數的顯著性均小于0.01，表明在7個波段的回歸方程在擬合過程中所有自變量對因變量均有顯著影響，根據回歸方程系數對不同波段的葉面積指數光飽和值進行求解，并根據回歸方程擬合得到橡膠林葉面積指數在不同波段反射率下的光譜飽和曲線（圖4），從而得到不同波段光譜“飽和點”，即葉面積指數飽和值見表4。

由圖4、表4可知，Landsat 8 OLI的7個單波段對應的橡膠林葉面積指數的光譜反射率大小在0～0.3 μm之間，葉面積指數飽和值大小介于5.08～7.51之間，最大值出現在近紅外波段（B5），飽和值為7.51；最小值出現在海岸波段（B1），飽和值為5.08；除近紅外波段（B5）飽和值較大之外，其余波段的飽和值均介于5～6之間；另外，在可見光范圍內（B1～B4）葉面積指數飽和值隨著波長的增加而增加，進而達到飽和；在近紅外（B5）和短波紅外波段（B6、B7）葉面積指數飽和值隨著波長的增加而減小，最后達到飽和。

圖4 光譜反射率變化曲線Fig. 4 Spectral reflectance change curve

表4 Landsat 8 OLI單波段對應的葉面積指數飽和值Table 4 Leaf area index saturation values for the single band of Landsat 8 OLI

4 結論與討論

本研究以Landsat 8 OLI數據為信息源，在對LiDAR激光云預處理的基礎上利用LiDAR 360軟件提取點云條帶區橡膠林有效葉面積指數，并提取1 649個點作為地面樣本數據，借助地統計學半方差函數理論，對Landsat 8 OLI各原始波段葉面積指數光譜飽和特征進行分析，結果表明：

1）Landsat 8 OLI近紅外波段的反射率隨著葉面積指數的增大而增大進而飽和，其余波段反射率隨著葉面積指數的增大而減小進而飽和，在可見光范圍內葉面積指數飽和值隨著波長的增加而增加；在近紅外和短波紅外波段葉面積指數飽和值隨著波長的增加而減小，雖然呈現出一定的規律性，但是差異并不是很明顯。

2）Landsat 8 OLI的7個原始波段對應的葉面積指數飽和值大小介于5.08～7.51之間，最大值出現在近紅外波段（B5），飽和值為7.51；最小值為海岸波段（B1），飽和值為5.08，除近紅外波段飽和值較大之外，其余波段的飽和值均介于5～6之間。

目前，地面葉面積指數的獲取主要通過直接測量和間接測量。直接測量為破壞性取樣法、異速生長方程法和凋落物收集法，此方法測量精確，但費時費力，但也存在一定誤差；間接測量主要借助于光學儀器，方便快捷且不具破壞性，但存在明顯的低估現象。本研究以LiDAR點云提取的LAI作為先驗數據樣本，此方法具有一定的創新性和可行性；從樣本數量上來說，本次提取了1 649個點作為樣本點，樣本數量有了一定的保證，研究結果雖然呈現出一定的規律性，但是差異并不是很明顯。本次研究對象為橡膠林，有關Landsat 8 OLI橡膠林LAI光譜飽和特征分析的研究較少，但是對比武紅敢等[5]、Zhao等[6]、Baret等[7]的研究結果來看，本研究的結果具有一定的合理性。同時，本研究也存在不足之處：1）結合數據源來看，地面LAI是從LiDAR點云中提取，此方法雖然已經成熟并且被廣泛應用，但是缺乏人工地面實測LAI數據對LiDAR提取的LAI以及LAI光譜響應特征的變化進行驗證；2）樣本點的選擇過程存在不可控因素，Landsat 8 OLI和LiDAR數據源在像元空間匹配上并不能做到絕對的匹配，只能以周圍像元平均值代替，從而造成光譜反射率不能完全反應地面LAI特征；3）LiDAR點云計算葉面積指數時，參與計算的激光點均是高于2 m的植被點，若在橡膠林郁閉度較低的情況下，灌草的滲入導致反射率增加或者減小，此時就造成橡膠林反射率偏大或者偏小，這也是影響橡膠林LAI飽和特征的重要因素；4）在飽和值求算過程中本研究只使用半方差函數理論，方法較為單一。

多源數據的協同相較于傳統的的森林資源監測節省了大量的人力、物力、財力，有著廣闊的應用前景，在今后的研究中應該考慮添加地面實測數據分析LAI的光譜飽和特征，從而對LiDAR點云中提取LAI光譜飽和特征進行驗證，另外，使用微分或者其他的飽和點求算方法進行對比。