楊澤至，舒清態(tài)，邱 霜，席 磊

(西南林業(yè)大學 林學院，云南 昆明 650224)

葉面積指數(shù)(leaf area index,LAI)作為影響植被光合作用、呼吸作用、蒸騰作用和碳氮水循環(huán)等生理過程的關鍵因素以及反映植被凈初級生產力和物種多樣性水平的重要因子，是利用大氣物理模型模擬植被與氣候相互作用關系的關鍵輸入?yún)?shù)[1-4].大尺度森林LAI 分布信息對于深入揭示地表植被的時間變化對氣候變化的響應機制，進而深化對植被與氣候相互作用關系的認識至關重要[1].因此，如何大面積科學、高效地獲取森林LAI 空間分布信息是地學領域長期以來的研究熱點.傳統(tǒng)的LAI地面直調查方法雖然能夠取得較高的測量精度，但存在耗時長、成本高、僅限于“點狀”尺度等難以克服的固有缺陷，難以獲取LAI 的空間格局信息，阻礙了LAI 的具體應用[5].當前，遙感技術憑借其快速獲取大范圍空間連續(xù)數(shù)據(jù)的技術優(yōu)勢已經成為快速提取森林LAI 空間分布信息的主流手段.包括光學遙感、機載激光雷達和合成孔徑雷達等在內的多種遙感技術在植被LAI 定量反演中得到了不同程度的應用，尤其以基于被動光學數(shù)據(jù)的LAI估測應用最為廣泛[5].然而，基于以上遙感數(shù)據(jù)的森林LAI 的估算方法均存在各自的局限性.光學遙感雖然能獲取大尺度的空間連續(xù)數(shù)據(jù)，被廣泛應用于多種生境條件下的LAI 估測，但植被指數(shù)等光譜變量常出現(xiàn)光飽和效應，進而容易導致LAI 的高值低估[6-7].機載激光雷達作為當前快速提取森林結構參數(shù)的先進遙感手段，盡管具有高精度獲取不同高度的植被結構信息的獨特優(yōu)勢，可實現(xiàn)LAI 的高精度反演，但其硬件昂貴、數(shù)據(jù)獲取成本高，且通常只能獲取范圍較小的帶狀數(shù)據(jù)，使其技術推廣受到制約[8-10].合成孔徑雷達雖然不受氣象條件限制，能夠全天時進行植被LAI 建模，但模型精度易受到傳感器特性、冠層結構特征及地表性質等諸多因素的影響，因而不同程度限制了該技術的廣泛應用[5].隨著遙感技術進一步發(fā)展，星載激光雷達為森林結構參數(shù)大尺度高效提取提供了新的技術途徑[10].然而，以星載激光雷達作為主要數(shù)據(jù)源進行森林LAI 估測的研究尚不多見，缺乏對星載激光雷達估算森林LAI 潛力和效果的充分探討.

新一代激光雷達衛(wèi)星ICESat-2（Ice，Cloud，and land Elevation Satellite-2）搭載的傳感器ATLAS（Advanced Terrain Laser Altimeter System）是全球首個光子計數(shù)激光雷達.作為新型星載激光雷達系統(tǒng)，與全波形激光系統(tǒng)GLAS 相比，ATLAS 光子數(shù)據(jù)具有更小的光斑直徑（17 m）、更高的脈沖頻率（10 kHz）和光子密度，在森林資源監(jiān)測中具有重要的應用潛力[11].目前，ICESat-2 數(shù)據(jù)已逐漸在極地冰蓋與海冰高程檢測、地表高程測量、湖泊水位動態(tài)監(jiān)測以及森林結構參數(shù)提取等多個領域發(fā)揮著關鍵作用[12].在森林結構參數(shù)提取方面，ICESat-2 數(shù)據(jù)的應用研究主要集中在區(qū)域尺度的林分樹高分布制圖和森林地上生物量反演兩個方向[12-15]，而應用于森林LAI 估測的研究非常有限，鮮見報道.因此，基于ICESat-2 數(shù)據(jù)估測森林LAI 的有效性和潛力仍有待深入探討.此外，與光學影像和機載LiDAR 點云數(shù)據(jù)不同，ICESat-2 光子數(shù)據(jù)缺乏空間連續(xù)性，無法直接形成連續(xù)的面狀或帶狀數(shù)據(jù)，而必須借助一定數(shù)據(jù)外推方法才能使目標反演參數(shù)由點狀分布拓展到面狀分布，進而形成其預測性空間分布.目前，實現(xiàn)遙感變量或反演結果外推的方法主要包括兩大類，即建立反演結果與衛(wèi)星光學影像特征因子的統(tǒng)計關系[16-17]和克里格插值[18].其中，克里格插值在遙感變量及其反演結果外推中的應用較少，作為對未抽樣點的線性無偏最優(yōu)估計[19]，其不僅可以有效克服光學影像所固有的易飽和缺陷，而且能夠實現(xiàn)ICESat-2 光子數(shù)據(jù)的充分利用.

綜上，利用ICESat-2/ATLAS 數(shù)據(jù)進行森林結構參數(shù)反演具有重要的應用潛力，其在大尺度森林LAI 估測中的有效性需待深入探討.對此，本文提出了基于ICESat-2 數(shù)據(jù)并結合克里格插值，實現(xiàn)森林LAI 連續(xù)估測的方法.其基本思路為：結合LAI 樣地的對應ATLAS 光斑特征因子與其47 個LAI 實測數(shù)據(jù)，分別利用逐步回歸和隨機森林回歸算法構建森林LAI 估測模型，隨后通過克里格插值實現(xiàn)光斑特征因子屬性的尺度外推，進而實現(xiàn)基于ICESat-2/ATLAS 數(shù)據(jù)的研究區(qū)森林LAI 的整體預測.基本研究過程如圖1 所示.

圖1 主要研究思路示意圖Fig.1 Flow chart of main research ideas

1 數(shù)據(jù)與方法

1.1 研究區(qū)概況本研究以地處滇西北的香格里拉市為試驗區(qū).香格里拉市位于云南省西北部、迪慶藏族自治州中東部，西鄰德欽縣和維西傈僳族自治縣，南接麗江市，東部與四川省稻城、木里二縣毗鄰，跨越99°8′2″E～100°21′15″E 和26°49′7″N～28°54′38″N，是青藏高原向云貴高原的過渡地帶和橫斷山區(qū)的腹地（圖2）.香格里拉市境內地形復雜多樣，以山地地形為主，地勢中間低，四周高，平均海拔3 000 m 以上，多種氣候帶鑲嵌交錯，土壤類型多樣.該地區(qū)森林覆蓋率高，地處亞熱帶常綠闊葉林植被帶與青藏高原高寒植被帶的交錯地帶，低緯高原的地理特征造就了其顯著的植被垂直分異特征，主要樹種為高山松Pinus densata、冷杉Abies fabri、云杉Piceaasperata、和落葉松Larix gmelinii等，并伴有楊樹Populus tremula等伴生樹種，尤其以高山松分布最為廣泛，是我國重要的生態(tài)單元之一，具有重要的研究意義.

圖2 研究區(qū)位置及LAI 樣地分布Fig.2 Location of the study area and distribution of LAI sample sites

1.2 數(shù)據(jù)獲取及處理研究數(shù)據(jù)主要包括ICESat-2/ATLAS 的2 種數(shù)據(jù)產品：全球地理定位產品ALTA03 數(shù)據(jù)陸地與植被冠層高度產品ATLA 08數(shù)據(jù)；試驗區(qū)ATLAS 光斑尺度下的森林LAI 樣地數(shù)據(jù).

1.2.1 ICESat-2/ATLAS 數(shù)據(jù)獲取 本研究采用ATLAS 采集的ATLA 03 和ATLA 08 數(shù)據(jù)作為研究區(qū)森林LAI 估測的遙感數(shù)據(jù)源，以探討其在區(qū)域尺度的森林LAI 估測中的潛力.目前，ATLAS 按照4 個級別（Level0～Level3）已生成了共計21 種產品數(shù)據(jù)集，依次編號為ATL00～ATL21[20-21].研究所采用的覆蓋整個試驗區(qū)的ATL03、ATL08 產品數(shù)據(jù)集均由美國國家冰雪中心官網（https://nsidc.org/data/ ATL03/versions/1，https://nsidc.org/data/ATL08/versions/1https://nsidc.org/data/icesat-2/datasets）下載得到，并與試驗區(qū)LAI 樣地實測數(shù)據(jù)保持時間同步.為覆蓋研究區(qū)整個林地區(qū)域，本試驗獲取了2020 年1 月至2021 年6 月之間香格里拉市境內的所有ATL03 和ATL08 產品數(shù)據(jù)集.其中，ATL03 數(shù)據(jù)共計118 條記錄，含354 條軌道、708條光子軌道波束；ATL08 數(shù)據(jù)共計118 條記錄，含354 條軌道、708 條光子軌道波束.

1.2.2 研究區(qū)ATLO8 光子數(shù)據(jù)去噪和分類處理

（1）光子點云去噪 ATLAS 光子數(shù)據(jù)在采集過程中易受太陽輻射、大氣散射、地形條件和傳感器特性等多種因素影響，包含大量非植被光子的背景噪聲.因此，需進行強背景噪聲的過濾，才能確保光子分類和地表信息提取的有效性[22-23].研究表明，基于改進OPTICS 的光子點云去噪算法具有較強的自適應性，且對聚類參數(shù)不敏感，尤其適用于存在強背景噪聲和地形條件復雜的山地，并在分離信號光子與噪聲光子、估算地面高程及樹高方面表現(xiàn)均優(yōu)于改進的DBSCAN 算法[12,24].本研究區(qū)地形復雜，山地廣布，嚴重影響信號光子信息的高精度提取.因此，本文采用基于改進OPTICS 的光子點云去噪算法對研究區(qū)的ATL08 數(shù)據(jù)進行噪聲濾除，以有效提取森林信號光子.該算法很好利用了噪聲光子點和信號光子點在空間分布規(guī)律上的顯著差異，具有較強機理性.其去噪過程如圖3所示.

圖3 基于改進OPTICS 的光子點云去噪過程[13]Fig.3 Photonic point cloud denoising process based on improved OPTICS[13]

（2）光子點云分類 經過噪聲濾除后的ATLAS 光子數(shù)據(jù)包含了植被冠層光子和地面光子等信號光子.因此，需對信號光子進行分類處理，以有效提取森林冠層光子.研究顯示，經改進的漸進三角不規(guī)則網絡加密[progressive triangular irregular network (TIN) densification，PTD]光子點云分類算法在去除低層物體上的非地面點和在復雜環(huán)境中保留地面特征點方面比經典PTD 方法表現(xiàn)得更好.與經典的PTD 算法相比，改進型PTD 方法具有更低的分類誤差和更高的光子分類精度，對地形條件復雜的區(qū)域具有較強的適用性[25-26].因此，本研究基于改進的PTD 光子分類算法對研究區(qū)去噪后的ATLAS 光子數(shù)據(jù)進行分類處理，進而有效提取植被冠層光子.該算法包括3 個核心步驟：參數(shù)指定、種子點選擇和初始TIN 構造以及TIN 的迭代加密.

通過光子去噪和分類處理后，獲得覆蓋香格里拉市的上千萬個ATL03 產品的有效光子，并基于ATL08 產品的數(shù)據(jù)組織結構，對所獲取的ATL03產品光子進行100 m 的抽稀采樣，進而提取得到研究區(qū)內有效光斑94 086 個.隨后根據(jù)研究區(qū)最新林地范圍Shp 文件提取得到覆蓋研究區(qū)內有效林地光斑74 230 個，非林地有效光斑19 856 個.最終的研究區(qū)ATLAS 光斑分布如圖4 所示.

圖4 研究區(qū)林地ATLAS 光斑分布Fig.4 ATLAS spot distribution in the study area

1.2.3 特征因子提取與優(yōu)選 本文采用Python語言，基于Jupyter 集成開發(fā)環(huán)境建立了ICESat-2/ATLAS 特征參數(shù)提取模塊程序，以提取模型特征因子，并采用皮爾遜相關分析法對ICESat-2/ATLAS-LAI 模型因子進行優(yōu)選.

1.2.4 樣地數(shù)據(jù)采集及處理

（1）LAI 的樣地設計 為了與ATLAS 光斑位置和大小完全相一致，本研究采用光斑腳印在地面上的投影作為調查樣地，即LAI 樣地是半徑為8.5 m 的樣圓，并以ATLAS光斑的圓心作為樣地的幾何中心.在調查過程中利用千尋星矩SR3 差分定位儀對光斑中心點坐標進行精確采集，以確保LAI 樣地與ATLAS 光斑在空間上的重合.本次樣地調查采用分層隨機抽樣法進行樣地布設，共調查了覆蓋研究區(qū)主要森林類型的LAI 樣地54 塊，整個樣本數(shù)據(jù)能夠有效反映研究區(qū)森林類型的分布情況.

（2）樣地實測LAI 的計算 采用半球攝影法實現(xiàn)樣地LAI 實測數(shù)據(jù)的獲取[10].首先確保光照條件穩(wěn)定，天色光較為均勻，且無陽光直射，隨后對每個樣地均勻選取13 個點位，并逐一在每個點位上將魚眼鏡頭呈豎直向上各取1 張林冠魚眼半球照片，形成13 張魚眼照片，即每個LAI 樣地均對應13 張林冠半球照片，以此作為計算樣地LAI 值的基礎.對照片質量合格的47 塊樣地的魚眼照片進行預處理，分別統(tǒng)一照片行列號數(shù)和格式為3 000×4 000 和JPG 格式.最后，利用基于MATLAB 程序開發(fā)的 CANEYE V6.1 的半球圖像分析軟件對47塊有效樣地的魚眼照片進行LAI 計算，進而得到47 塊樣地的LAI 實測值，其描述統(tǒng)計特征見表1.

表1 有效樣地實測LAI 描述統(tǒng)計Tab.1 Descriptive statistics of measured LAI in effective sample plots

1.3 研究方法

1.3.1 建模方法 研究分別利用成熟、穩(wěn)健的逐步回歸和隨機森林回歸[27-30]（random forest regression，RFR）模型構建森林LAI 估測模型.

1.3.2 變異函數(shù) 變異函數(shù)又稱半方差函數(shù)，其定量描述了地理數(shù)據(jù)的空間自相關性，可以表征區(qū)域化變量的空間分布格局，是地質統(tǒng)計學的核心工具和克里格插值的理論基礎[19].在地學變量滿足本征假設和二階平穩(wěn)的假設下，變異函數(shù)可定義為區(qū)域化變量在對應采樣間隔下樣本方差的數(shù)學期望，其函數(shù)表達式為：

式中:γ(h)為變異函數(shù)，h為成對樣點的采樣間距，n(h)是當樣點間距為h時的樣點總對數(shù)，z(xi)和z(xi+h)分別為區(qū)域化變量在空間位置xi和xi+h處的屬性值.

經變異函數(shù)計算得到的實驗變異函數(shù)需通過球狀模型、指數(shù)模型線性模型和高斯模型等方法進行擬合，建立其理論變異函數(shù)模型，從而為地學變量的結構分析及其空間局部估計即克里格插值提供理論基礎[19].變異函數(shù)理論模型的結構特征由函數(shù)類型、塊金方差（Nugget variance，記作C0）、基臺值（Sill，記作C0+C，C為偏基臺值）和變程（Range，記作a）等4 個參數(shù)來表征.其中，塊均方差指樣本間距h=0 時的變異函數(shù)值，來自測量誤差或局部微小變異；基臺值是當采樣間距增加到某一臨界值時變異函數(shù)趨于穩(wěn)定而所對應的常數(shù)；而變程則是變異函數(shù)達到基臺值時的樣本間隔距離，是地學變量空間自相關性質的作用距離.

1.3.3 克里格插值 克里格插值是建立在地學變量數(shù)據(jù)分布及其變異函數(shù)結構分析基礎上，對未知點屬性值的無偏最優(yōu)估計.其插值過程中既充分考慮未知點與所有已知樣本點之間的空間位置關系，又充分利用各樣本點數(shù)據(jù)的空間結構特征，對每個已知點賦予一定的插值權重，并采用加權平均的方法對未知點的值進行預測[19,31].本研究采用克里格插值法，以光斑采樣尺度（17 m×17 m）為插值尺度，對擬參與LAI 建模的ICESat-2 光子特征參數(shù)進行線性無偏最優(yōu)估計，進而形成對應的光子特征參數(shù)的柵格圖層.其插值表達式如下：

式中：Z0為擬參與建模的ICESat-2 光子特征參數(shù)在未知點的屬性值，λi為參與插值的樣地點對估測點屬性的貢獻權重，Xi為樣地點位置，Z(Xi)為ICESat-2 光子特征參數(shù)在位置Xi處的觀測值.在確定已知點的貢獻權重時不僅要參照估測點和已知點之間的距離，還需考慮到已知點屬性值及其空間分布格局.

1.3.4 模型評價方法 以可決系數(shù)R2和均方根誤差(root mean square error,RSME)作為模型精度評價指標，并采用留一交叉驗證法（leave-one-out method of cross-validation,LOOCV）對模型性能進行驗證評估，以有效避免建模數(shù)據(jù)集和測試數(shù)據(jù)集劃分過程中LAI 建模樣本選取的隨機性和不確定性[32]，提高小樣本建模結果評價的可靠性.相關統(tǒng)計指標的表達式如下所示：

式中:yi為因變量實測值，為模型的預測值，為因變量實測值的均值，n為樣本數(shù).

2 結果與分析

研究結果包括模型因子相關性分析、森林LAI 估測模型的建立、LAI 模型特征因子克里格插值和研究區(qū)森林LAI 空間分布預測等4 個部分.

2.1 ATLAS08 光斑特征參數(shù)與樣地實測LAI 相關性分析在提取47 個LAI 樣地對應的ICESat-2 光斑特征參數(shù)的基礎上，通過SPSS 25.0 版統(tǒng)計軟件，采用皮爾遜相關分析法（Pearson correlation analysis）對47 個樣地中的LAI 實測值與其對應的56 個光斑特征因子進行線性相關分析，以從中優(yōu)選出與樣地實測LAI 相關性較為顯著的光斑特征因子作為LAI 的建模因子，其優(yōu)選結果見表2.Pearson 相關性分析結果顯示,與樣地實測LAI 相關性最高光斑特征因子為Landsat percent，相關系數(shù)為0.346 (P＜0.05)；其次為冠層光子比例（Canopy photon rate），相關系數(shù)為0.327（P＜0.05）；相關性最低光斑特征因子為信噪比SNR，相關系數(shù)為0.01.共有4 個光斑特征因子與樣實測LAI 的相關性較為顯著，按其相關系數(shù)由大至小依次為排列為：Landsat percent ＞ Canopy photon rate ＞ Minimum canopy＞ Number of canopy photons.因此，本文將上述4個ICESat-2 光斑特征因子作為LAI 的模型因子.

表2 ATLAS 特征因子優(yōu)選結果Tab.2 ATLAS feature factor preference results

2.2 LAI 模型預測結果與分析根據(jù)2.1 節(jié)的光斑特征因子優(yōu)選結果，以Landsat percent、Canopy photon rate、Minimum canopy height 和Number of canopy photons 等4 個與實測LAI 相關性相對顯著ICESat-2/ATLAS 特征變量為建模因子，以樣地實測LAI 為模型因變量，分別采用逐步線性回歸和隨機森林回歸算法，建立基于ICESat-2 光子數(shù)據(jù)的森林LAI 遙感估測模型.其建模結果見表3～表7.

2.2.1 基于逐步回歸的LAI 估測模型 利用Minitab 統(tǒng)計分析軟件，以經優(yōu)選的4 個光斑特征因子為模型自變量，以47 個樣地實測LAI 真值為響應變量，建立森林LAI 估測模型，以檢測ICESat-2 光子數(shù)據(jù)與森林LAI 的線性關系水平.其建模結果表明（表3），逐步回歸分析一共產生了3 個LAI線性估測模型.最高LAI 擬合精度為0.389，表明入選建模的ICESat-2 光斑特征因子能夠解釋樣地LAI 的38.9%的方差變化.由其他因素引起LAI 方差變化占到了61.1%.由表4 并結合模型擬合優(yōu)度及顯著性可知，LAI 的最優(yōu)線性模型為：LAI=-0.074+0.006×Landsat percent+0.242×Canopy photon rate+0.11×Minimum canopy height，其LOOCV 評價見表5.

表4 LAI 逐步回歸模型參數(shù)Tab.4 Parameters of LAI stepwise regression model

表5 LAI 逐步回歸模型留一交叉評價Tab.5 LAI stepwise regression model leave-one-out crossevaluation

2.2.2 基于隨機森林回歸的LAI 估測模型 為探究森林LAI 與ICESat-2/ATLAS 光斑特征因子的非線性映射關系，研究借助RFR 算法構建了基于ICESat-2 光斑特征因子的森林LAI 反演模型.基于2.1 節(jié)的光斑特征因子優(yōu)選結果，結合47 份LAI-ICESat-2 樣本數(shù)據(jù)，利用Python 語言，借助Jupyter book 平臺，通過Scikit-learn 機器學習庫（https://scikit-learn.org.cn）調用RandomForestRegressor 模塊實現(xiàn)了LAI 的隨機森林建模.其模型訓練的關鍵超參數(shù)設置見表6.結果表明，留一交叉驗證下的LAI 隨機森林模型的R2和RMSE 分別為0.557 和0.175（圖5、圖6）.可見，RFR 模型在整體上較好地擬合了ICESat-2 光斑特征因子與LAI 的關系.

表6 LAI 隨機森林回歸模型參數(shù)取值Tab.6 Parameter values of LAI random forest regression model

圖5 RFR 預測LAI 的效果點線圖Fig.5 Point line diagram of the effect of RFR to predict LAI

圖6 基于RFR 模型的LAI 預測精度Fig.6 Accuracy of LAI prediction based on RFR model

表7 顯示，留一交叉驗證下，LAI 的逐步線性回歸模型和隨機森林回歸模型的可決系數(shù)R2分別為0.187 和0.557；RMSE 分別為0.231 和0.175；隨機森林回歸模型對LAI 的整體預測效果優(yōu)于逐步回歸模型.因此，本研究以隨機森林回歸模型作為估測整個研究區(qū)森林LAI 分布格局的建模方法.

表7 LAI 估測模型精度對比Tab.7 Comparison LAI estimation model accuracy

2.3 LAI 模型光斑特征因子克里格插值結果獲取研究區(qū)ICESat-2/ALTAS 光子特征參數(shù)的連續(xù)分布數(shù)據(jù)是實現(xiàn)研究區(qū)森林LAI 空間分布制圖的基礎條件.為此，本文基于光斑特征因子優(yōu)選結果，采用變異函數(shù)對經優(yōu)選并參與LAI 建模的4 個光斑特征因子（表8）進行空間結構分析，并通過克里格插值對4 個光斑特征因子的屬性值進行連續(xù)內插，繼而得到4 個光斑特征因子的柵格圖層.本節(jié)的變異函數(shù)分析通過GS+9.0 軟件實現(xiàn)，該軟件為用戶提供了地學變量的變異函數(shù)值計算及其理論模型的擬合工具.4 模型因子的變異函數(shù)結構分析及其克里格插值結果分別見表9 和圖7.

表8 待插值的4 個光斑特征因子Tab.8 The 4 spot feature factors to be interpolated

表9 LAI 模型光斑特征因子變異函數(shù)結構分析Tab.9 Structural analysis of the variance function of the LAI model light spot feature factor

圖7 研究區(qū)LAI 光斑特征因子屬性空間分布Fig.7 Spatial distribution of LAI spot feature factor attributes in the study area

變異函數(shù)結構分析表明，上述4 個ICESat-2 光子特征屬性具有高度的空間異質性,且表現(xiàn)出不同的變異特征.研究區(qū)4 個光斑特征因子中除了冠層光子數(shù)（Number of canopy photons）以外，其余3 個光斑特征因子的變異函數(shù)最優(yōu)擬合模型均為指數(shù)模型，其可決系數(shù)R2均在0.75 以上，且殘差平方RSS 和在4 個變異函數(shù)理論模型中均為最小值.對于Number of canopy photons，其球狀模型和高斯模型的決定系數(shù)R2均為0.68，其高斯模型的殘差略高于球狀模型，因而球狀模型是其變異特征的最佳表征模型.此外，塊金方差與基臺值之比[C0/(C0+C)]是描述地學變量總空間變異性基本構成的度量，直接反映了隨機因素對總空間異質性的影響強度，亦稱塊金效應.塊金效應可分3 個水平：[75%，100%]；[25%，75%）；[0，25%），依次分別代表地學變量呈強空間自相關、中等空間自相關和弱空間自相關.此效應越強，表明隨機因素對總空間異質性的影響越大，反之，表明總空間異質性中由變量空間自相關性引起的部分占比越高.塊金效應表明，4 個光斑特征因子在研究區(qū)內由其空間自相關性引起的異質性在總空間異質性的占比均達80%以上，而由隨機因素引起的空間變異性僅占到了總空間異質性的約6%～12%.4 個光斑特征因子在研究區(qū)內的空間自相關作用距離，即變程a處于5 700～7 500 m.可見，研究區(qū)4 個光斑特征因子的空間分布特征主要由ICESat-2 光子數(shù)據(jù)本身的結構性因素（如地、植被覆蓋和冠層結構特征等）控制，即空間結構特性是形成其異質性的決定性因素，而受隨機因素影響較小.

4 個光斑特征因子屬性的變異特征分析顯示，LAI 模型的4 個光斑特征因子屬性均具有顯著的空間自相關性，具備對其進行克里格插值的充分條件.因此，本文在變異結函數(shù)結果分析的基礎上，借助ArcGIS 10.7 中的地統(tǒng)計分析（Geostatistical Analyst）模塊分別對Landsat percent、Canopy photon rate、Minimum canopy height 和Number of canopy photons 等4 個光斑特征因子的屬性進行克里格插值，進而形成4 個特征因子的柵格屬性圖層，如圖7所示.

2.4 研究區(qū)森林LAI 空間分布預測基于2.2 節(jié)LAI 建模結果和2.3 節(jié)的ICESat-2 光斑特征因子屬性的克里格內插結果進行香格里拉市的森林LAI 的空間分布預測.以研究區(qū)林地范圍內的光斑特征因子連續(xù)分布數(shù)據(jù)為自變量，利用RFR 模型逐像元地對研究區(qū)森林LAI 進行預測，并以17 m輸出分辨率，以確保LAI 的制圖尺度與ICESat-2的光斑采樣尺度一致，進而實現(xiàn)研究區(qū)森林LAI的連續(xù)制圖（圖8）.

圖8 研究區(qū)森林LAI 預測分布圖Fig.8 Forest LAI prediction distribution map in the study area

圖8 所示的LAI 反演結果較好地反映了研究區(qū)內森林LAI 的空間分布，直觀、有效地展示了森林LAI 的空間分布差異，與ICESat-2 光斑特征因子一樣，具有高度的空間異質性.總體而言，研究區(qū)的LAI 整體分布不均衡，整體呈現(xiàn)“中間高、四周低”的分布格局.LAI 高值區(qū)集中分布在香格里拉市的中部偏北區(qū)域和南部區(qū)域，且多為森林腹地，而LAI 低值區(qū)主要分布在其中部、與城鎮(zhèn)毗鄰的林地區(qū)域，符合地理學第一定律（空間自相關性）和地理學第二定律（空間異質性）.

3 討論與結論

3.1 討論

（1）基于ICESat-2/ATLAS 光子數(shù)據(jù)估測森林LAI 具有重要的應用潛力.本研究基于新一代星載激光雷達ICESat-2 光子數(shù)據(jù)進行森林LAI 估測，以探討星載光子數(shù)據(jù)在森林LAI 建模中的潛力，是區(qū)域尺度森林LAI 反演方法的新嘗試.在橫向對比上，本研究得到的研究區(qū)森林LAI 空間分布格局與譚德宏等[33]和趙洪瑩等[34]等利用衛(wèi)星光學數(shù)據(jù)對同一區(qū)域森林LAI 的預測分布高度一致，尤其在LAI 低值區(qū)更為吻合.這從側面反映了ICESat-2 光子數(shù)據(jù)在森林LAI 估測中的可行性.此外，與Zhang 等[35]同樣利用ICESat-2 數(shù)據(jù)并結合Beer-Lambert 定律進行大尺度森林LAI 反演不同，本文以林地ICESat-2 光斑內的森林LAI 的實測樣本與ICESat-2 光子數(shù)據(jù)的統(tǒng)計關系進行建模，估測結果能夠進行更為科學的地面驗證，因而具有顯著的生態(tài)學意義.在縱向對比上，ATLAS 光子計數(shù)激光雷達的數(shù)據(jù)采集性能得到顯著提升，能獲得直徑更小且密度遠大于第一代星載激光雷達GLAS全波形數(shù)據(jù)的光斑數(shù)據(jù)，為大尺度森林結構參數(shù)的快速提取提供了更為有效的數(shù)據(jù)支撐[11].早期研究，如駱社周等[16]利用直徑為70 m 的GLAS 光斑數(shù)據(jù)，基于LPI 激光穿透指數(shù)建立LAI 估測模型，并結合TM 光學影像進行了大尺度森林LAI 反演，是基于星載激光雷達數(shù)據(jù)反演區(qū)域森林LAI 的一次重要嘗試.然而，該研究受制GLAS 全波形數(shù)據(jù)自身特性，存在地面采樣成本較高，反演精度地域性強等問題，其方法技術的可移植性不強.

（2）隨機森林回歸模型對LAI 的估測效果優(yōu)于逐步回歸模型.隨機森林和逐步回歸都是性能穩(wěn)健、并被廣泛應用的建模方法.本文分別采用逐步回歸和隨機森林模型進行研究區(qū)森林LAI 建模，并對比兩種模型對LAI 的預測效果.留一交叉驗證顯示，逐步回歸和隨機森林回歸模型對研究區(qū)森林LAI的預測精度R2分別為0.187 和0.557，RMSE 分別為2.34 和1.75，可見隨機森林回歸模型的預測效果優(yōu)于逐步回歸.然而，本文應用的兩種模型對研究區(qū)森林LAI 的估測精度均不高，其主要原因在于LAI 實測樣本容量有限，限制了隨機森林模型預測性能的充分發(fā)揮，不易深入發(fā)現(xiàn)ICESat-2 光子數(shù)據(jù)與實測LAI 的內在關系.然而，較大的實測樣本量，必然以大量財力物力的消耗為代價，降低了研究的可行性與遙感估測模型的實用性.因此，如何在有限條件下，充分發(fā)掘遙感數(shù)據(jù)與LAI 實測值的潛在關系、積極探討新興人工智能方法（如梯度提升GBTD[36]、極端隨機森林Extree[37]或深度學習算法[38]）在基于ICESat-2 數(shù)據(jù)估測森林LAI中的潛力將是今后的研究重點.

（3）克里格插值法能夠充分利用ICESat-2 數(shù)據(jù)對LAI 估測結果進行外推.在利用星載激光雷達進行森林LAI 反演的早期研究中，有關學者借助衛(wèi)星光學影像實現(xiàn)基于星載激光雷達的反演結果的尺度上推，以使反演結果呈連續(xù)分布.如駱社周等[16]利用TM 影像將基于GLAS 波形數(shù)據(jù)反演得到的森林LAI 外推至試驗區(qū)，從而形成了LAI 的預測分布圖.Narine 等[15]借助Landsat-8 影像將基于ICESat-2 數(shù)據(jù)的AGB 反演結果尺度上推至整個研究區(qū).然而，此類方法對估測模型的預測精度要求較高，且會造成估測結果外推過程中預測誤差的二次傳遞，并帶來光學數(shù)據(jù)的飽和問題，適用范圍有限.針對上述問題，本研究基于變異函數(shù)結構分析，對經優(yōu)選的4 個呈離散分布的ATLAS 光斑特征因子進行克里格插值，進而形成其各自屬性覆蓋研究區(qū)的連續(xù)數(shù)據(jù)，并通過RFR 實現(xiàn)了研究區(qū)森林LAI 的連續(xù)估測.本文方法不僅充分利用了ICESat-2 光子數(shù)據(jù)，而且有效避免了借助衛(wèi)星光學影像對LAI 估測結果進行尺度外推過程中新的誤差引入以及其帶來的光飽和效應.此外，研究僅采用了普通克里格法對4 個參與LAI 建模的光斑特征因子進行插值，而未對多種克里格法的插值效果進行比較，其它克里格插值法，如協(xié)同克里格、泛克里格法在光斑特征因子插值中的效果需進一步探討.

此外，ATLAS 的強、弱波束光子在森林LAI預測中的效果需分別探究；RFR 模型參數(shù)的優(yōu)化與否對LAI 模型精度是否有顯著提升也需進一步探究，是今后研究的重點之一.

3.2 結論針對傳統(tǒng)光學遙感易受光飽和影響、機載激光雷達與微波數(shù)據(jù)獲取成本高等問題，本文提出了一種基于ICESat-2/ATLAS 光子數(shù)據(jù)進行森林LAI 估測的新方法.該方法在研究區(qū)林地光子點云數(shù)據(jù)去噪、分類等預處理的基礎上，提取了各個林地光斑的特征因子，并結合47 份森林LAI 樣本數(shù)據(jù)，建立了基于樣地光斑尺度的森林LAI 隨機森林估測模型，其留一交叉驗證下的R2和RMSE分別為0.557 和1.75.由此表明，ICESat-2 光子數(shù)據(jù)在區(qū)域尺度的森林LAI 預測中具有重要潛力.為實現(xiàn)研究區(qū)森林LAI 的連續(xù)估測，本文基于變異函數(shù)結構分析，對經優(yōu)選的4 個參與LAI 建模的光斑特征因子進行克里格插值，形成其各自屬性的對應柵格圖層，并利用已訓練的RFR 模型實現(xiàn)了研究區(qū)森林LAI 的整體估測，繼而形成了17 m 空間分辨率的研究區(qū)森林LAI 空間分布圖.本研究估測得到的森林LAI 在空間梯度上與同類研究以及研究區(qū)森林覆蓋度基本一致，結論可信，是以較小成本實現(xiàn)大尺度森林LAI 快速反演的一個新思路.