多源遙感數(shù)據(jù)協(xié)同的景洪橡膠林葉面積指數(shù)光飽和特征研究

羅洪斌 舒清態(tài) 王 強(qiáng) 王冬玲 字 李 謝福明

( 西南林業(yè)大學(xué)林學(xué)院，云南 昆明 650233)

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

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

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

1 研究區(qū)概況

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

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

2 材料與方法

2.1 數(shù)據(jù)來源與處理

2.1.1 機(jī)載激光雷達(dá)數(shù)據(jù)

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

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

2.1.2 Landsat 8 OLI影像數(shù)據(jù)

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

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

2.1.3 輔助數(shù)據(jù)

景洪市森林資源二類調(diào)查矢量小班圖層數(shù)據(jù)，利用ArcGIS對橡膠林分布區(qū)進(jìn)行裁剪。

2.2 有效葉面積指數(shù)的提取

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

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

2.3 訓(xùn)練樣本的選取

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

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

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

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

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

若令b0=C0，b1=3C/2BS，b2=-C/2BS3則球狀模型變?yōu)橐话愕亩嘣€性回歸模型：y=b0+b1×LAI+b2×LAI3，運用最小二乘法即可對模型參數(shù)進(jìn)行求解，從而獲得不同波段所對應(yīng)的葉面積指數(shù)飽和值。

3 結(jié)果與分析

3.1 相關(guān)性分析

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

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

3.2 葉面積指數(shù)光譜響應(yīng)分析

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

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

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

3.3 識別葉面積指數(shù)光譜飽和點

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

表3 回歸系數(shù)的顯著性Table 3 Significance of regression coefficients

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

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

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

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

4 結(jié)論與討論

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

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

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

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

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