魏智海，張志超，李崇貴*

(1.國(guó)家林業(yè)和草原局 西北調(diào)查規(guī)劃設(shè)計(jì)院，陜西 西安 710048；2.西安科技大學(xué) 測(cè)繪科學(xué)與技術(shù)學(xué)院，陜西 西安 710054)

郁閉度是森林資源調(diào)查中的重要林分因子，與此同時(shí)，它也是森林生長(zhǎng)狀況和結(jié)構(gòu)的重要指標(biāo)，對(duì)制定森林經(jīng)營(yíng)管理方案和森林總體規(guī)劃設(shè)計(jì)研究具有重要意義[1-3]。傳統(tǒng)森林郁閉度調(diào)查方法主要有目測(cè)法、樣線法、樣點(diǎn)法和樹冠投影法等，但隨著林業(yè)工程建設(shè)的穩(wěn)步推進(jìn)，傳統(tǒng)的郁閉度調(diào)查方法不僅耗時(shí)費(fèi)力，而且獲取的數(shù)據(jù)量非常有限，無(wú)法滿足大范圍內(nèi)郁閉度研究的需求[3]。

近年來(lái)，伴隨著衛(wèi)星技術(shù)的快速發(fā)展，利用遙感影像進(jìn)行森林郁閉度定量估測(cè)已經(jīng)成為一種新途徑。李崇貴等[4]基于Landsat TM影像采用最小二乘估計(jì)及嶺估計(jì)方法對(duì)以RS和GIS為基礎(chǔ)的森林郁閉度估測(cè)，研究發(fā)現(xiàn)當(dāng)影響郁閉度估測(cè)的變量間存在復(fù)共線性時(shí)，嶺估計(jì)具有較高精度；陳崇等[2]基于Landsat TM影像進(jìn)行亞熱帶森林郁閉度遙感估算，研究發(fā)現(xiàn)利用BP神經(jīng)網(wǎng)絡(luò)建立的森林郁閉度優(yōu)于逐步回歸模型；譚炳香[5]利用星載EO-1 Hyperion高光譜遙感數(shù)據(jù)進(jìn)行森林郁閉度定量估測(cè)，研究發(fā)現(xiàn)高光譜數(shù)據(jù)精細(xì)的光譜信息非常有利于提高植被參數(shù)的估測(cè)和反演精度；柏藝[6]基于3種不同數(shù)據(jù)源(Pleiades-1、SPOT-5和Landsat-8影像)，通過(guò)建立偏最小二乘法對(duì)川西常綠闊葉林的森林郁閉度進(jìn)行估測(cè)，研究發(fā)現(xiàn)對(duì)于川西南常綠闊葉林郁閉度的估測(cè)上，Pleiades-1影像更具有優(yōu)勢(shì)；吳飏等[1]采用逐步回歸方法，結(jié)合光譜、地形和紋理特征進(jìn)行郁閉度遙感估測(cè)，研究發(fā)現(xiàn)結(jié)合圖像紋理特征的方法比傳統(tǒng)的只基于光譜或地形特征的方法在估測(cè)精度上有很大提高；胡振華等[7]基于Hyperion數(shù)據(jù)采用偏最小二乘法對(duì)香格里拉森林郁閉度進(jìn)行遙感估測(cè)，研究發(fā)現(xiàn)基于樣地?cái)?shù)據(jù)郁閉度變化敏感性分析模型優(yōu)于基于Hyperion影像的分段主成分分析模型；李擎等[3]基于GF-2遙感影像采用逐步回歸分析對(duì)天山云杉林郁閉度定量估測(cè)，研究發(fā)現(xiàn)基于紋理特征、光譜信息、地形因子相結(jié)合所建立的模型效果最好。Pu等[8]利用小波變換的Hyperion影像，基于逐步回歸方法估測(cè)森林郁閉度，研究表明Hyperion影像估測(cè)精度可以達(dá)到近85%。然而，這些方法未對(duì)不同樣地?cái)?shù)據(jù)適用的模型進(jìn)行分析，并且未考慮不同的森林覆蓋類型對(duì)于郁閉度估測(cè)精度的影響。

因此，本研究以內(nèi)蒙古自治區(qū)東北部柴河林業(yè)局為研究區(qū)域，綜合利用該區(qū)域的GF-1遙感影像、二類調(diào)查小班數(shù)據(jù)、DEM數(shù)字高程模型以及林地一張圖數(shù)據(jù)，采用k-NN、穩(wěn)健估計(jì)以及偏最小二乘法建立以樣地為單位的森林郁閉度估測(cè)模型，輔之以二類調(diào)查小班數(shù)據(jù)和現(xiàn)地實(shí)測(cè)數(shù)據(jù)進(jìn)行模型精度檢驗(yàn)并挑選森林郁閉度最優(yōu)估測(cè)模型。

1 材料與數(shù)據(jù)

1.1 研究區(qū)概況

研究以位于內(nèi)蒙古自治區(qū)東北部柴河林業(yè)局為研究對(duì)象，范圍在120°36′-121°46′E， 47°05′-47°45′N(圖1)。柴河林業(yè)局屬寒溫帶大陸季風(fēng)氣候，年平均氣溫約為3 ℃，年平均降水量約為550 mm。柴河林業(yè)局地勢(shì)西高東低，地貌以中低山地貌為主，西部屬于中山地貌，東部屬于低山丘陵地貌[9]。該地區(qū)森林覆蓋率>75%，西部主要興安落葉松覆蓋為主，東部則以蒙古櫟、黑樺占絕對(duì)優(yōu)勢(shì)[10]。

圖1 研究區(qū)與樣地位置圖像Fig.1 Study area and locations oof sample sites

1.2 數(shù)據(jù)獲取與處理

1.2.1 遙感數(shù)據(jù)來(lái)源與預(yù)處理 遙感數(shù)據(jù)為通過(guò)中國(guó)資源衛(wèi)星應(yīng)用中心網(wǎng)站下載的一景GF-1遙感影像，傳感器及軌道號(hào)為：GF1_PMS1_E121.5_N47.4_20160920_L1A0001840656_PMS1。該衛(wèi)星上搭載了2臺(tái)2 m分辨率的全色/8 m分辨率的多光譜高分辨率相機(jī)[11]，并且它是國(guó)家高分辨率對(duì)地觀測(cè)系統(tǒng)中的首發(fā)星。因此，研究選用GF-1遙感影像的森林郁閉度定量估測(cè)對(duì)促進(jìn)國(guó)產(chǎn)遙感衛(wèi)星在森林資源調(diào)查和監(jiān)測(cè)中的應(yīng)用具有重要意義。在數(shù)據(jù)預(yù)處理方面，利用ENVI5.3SP1對(duì)影像進(jìn)行輻射定標(biāo)、大氣校正、正射校正、影像融合以及幾何精校正等預(yù)處理，最后通過(guò)疊加研究區(qū)域行政界線圖進(jìn)行裁剪。其中大氣校正使用FLAASH大氣校正，正射校正使用GF-1數(shù)據(jù)自帶的RPC文件和30M分辨率的DEM數(shù)據(jù)進(jìn)行，幾何精校正以17級(jí)谷歌衛(wèi)星圖(無(wú)偏移)為基準(zhǔn)。

地形特征數(shù)據(jù)為地理空間數(shù)據(jù)云(www.gscloud.cn)下載的DEM數(shù)據(jù)，其空間分辨率為30 m，該數(shù)據(jù)主要用于提取研究區(qū)的坡度(Slope)和坡向(Aspect)，亦可作為地形特征參與郁閉度定量估測(cè)。對(duì)研究區(qū)2017年林地“一張圖”數(shù)據(jù)進(jìn)行重分類，并通過(guò)ArcGIS生成有林地分布圖(圖2)。2019年9月在該影像范圍內(nèi)均勻選取具有代表性的47塊小班進(jìn)行外業(yè)現(xiàn)地實(shí)測(cè)，通過(guò)實(shí)測(cè)獲取到的小班郁閉度調(diào)查結(jié)果參與估測(cè)精度評(píng)價(jià)和實(shí)際應(yīng)用效果討論。

圖2 研究區(qū)林地分布Fig.2 Distribution of forest land in the study area

1.2.2 地面數(shù)據(jù)來(lái)源 用于研究的地面樣地?cái)?shù)據(jù)采用2016年全國(guó)森林資源二類調(diào)查小班數(shù)據(jù)。該數(shù)據(jù)通過(guò)在研究區(qū)內(nèi)按照4 km×4 km的間距布設(shè)解譯樣地，解譯樣地為0.06 hm2的長(zhǎng)方形(10 m×60 m)樣地，通過(guò)將解譯樣地所落小班的郁閉度擬定為樣地郁閉度，最終獲得有林地樣地54個(gè)。

2 研究方法

本研究依據(jù)樣地的基礎(chǔ)信息(位置、大小以及形狀)，通過(guò)疊加到處理好的GF-1影像提取得到樣地大小的影像，通過(guò)計(jì)算遙感信息、紋理信息和DEM提取的坡度、坡向等地形信息最終獲取到樣地大小的自變量，通過(guò)代入3種估測(cè)模型(穩(wěn)健估計(jì)、偏最小二乘法和k-NN法)計(jì)算每一個(gè)樣地的林分郁閉度估測(cè)值[12-14]，以此類推，最終獲得整個(gè)柴河林業(yè)局的森林郁閉度空間分布專題圖。

2.1 穩(wěn)健估計(jì)法

穩(wěn)健估計(jì)法不僅能夠有效消除或者減弱異常樣地對(duì)森林郁閉度估測(cè)模型的干擾能力，而且能夠有效剔除含有異常值的數(shù)據(jù)，從而使得回歸模型參數(shù)具有一定的“抗干擾性”和“穩(wěn)定性”[13,15]。人工有規(guī)則地挑選樣地會(huì)使得模型無(wú)法滿足普適性原則，因此，樣地挑選應(yīng)遵循隨機(jī)性原則。本研究采用隨機(jī)抽樣的方式進(jìn)行樣地抽取，并且采用m：n-m比例抽取訓(xùn)練數(shù)據(jù)集和驗(yàn)證數(shù)據(jù)集。假設(shè)樣地郁閉度為v，自變量個(gè)數(shù)為q，則

v=β1x1+β2x2+…+βqxq+β0

(1)

故m個(gè)樣地郁閉度V與q個(gè)自變量之間的關(guān)系可以表示為

V=Xβ+eV,(eV)=0,Cov(eV)=σ2I,V=(v1,v2,…,vm)′,eV=Xβ-V

(2)

穩(wěn)健估計(jì)需要滿足式(3)要求才能使用，即要求誤差的絕對(duì)最小值作為極值函數(shù)。

(3)

式中:ρ(ei)是極值函數(shù)，采用迭代函數(shù)法進(jìn)行求解，具體步驟及原理詳見文獻(xiàn)[15]。

2.2 偏最小二乘法

偏最小二乘法是一種新型的回歸方法。一方面，集成了主成分分析的優(yōu)點(diǎn)，能夠?qū)Χ鄠€(gè)自變量和因變量分別進(jìn)行分析；另一方面，采用了相關(guān)性最大化等原則進(jìn)行模型建立，可以有效解決變量多重相關(guān)性造成的不利影響，提高所模型模的可靠性[16-17]。具體步驟可參考文獻(xiàn)[15,18]。

2.3 k-最近鄰法

k-NN法是一種簡(jiǎn)單的機(jī)器學(xué)習(xí)算法，同時(shí)也是比較典型的非參數(shù)化方法。該方法是通過(guò)找出一個(gè)樣本的k個(gè)最近“鄰居”并將其屬性按照取平均的方式賦值給該樣本以此得到樣本屬性。同時(shí)，k-NN法與其他建模方法相比較，該建模方法不受自變量間不存在復(fù)共線性的限制[18-19]。公式如下

(6)

(7)

式中:Wpi,q為權(quán)重矩陣，本研究選用反距離函數(shù)計(jì)算得到；pi為訓(xùn)練樣地點(diǎn)，q為待估測(cè)點(diǎn)；Dpi,q為待估測(cè)點(diǎn)q和訓(xùn)練樣地點(diǎn)pi之間的距離，本研究選用歐氏距離進(jìn)行計(jì)算[18-20]；t為距離分解因子，本研究取t=2；upi為點(diǎn)pi的森林郁閉度；upi為待估測(cè)點(diǎn)q的森林郁閉度。

采用k-NN法建模時(shí)，選用交叉驗(yàn)證的方式消除研究區(qū)樣地?cái)?shù)量少的誤差，最終通過(guò)比較RMSE的方式優(yōu)選最佳k值[13,20-21]。RMSE(RMSE)計(jì)算公式如下：

(8)

3 結(jié)果與分析

3.1 自變量設(shè)置與變量?jī)?yōu)選

研究發(fā)現(xiàn)[1,7,3]，遙感、地形和紋理特征都會(huì)影響森林郁閉度估測(cè)，因此本研究綜合考慮這3個(gè)方面特征。地形特征包括坡度和坡向(陰坡和陽(yáng)坡)；遙感特征包括影像波段灰度信息和比值波段。影像波段灰度信息包括藍(lán)(B)、綠(G)、紅(R)、近紅外(Nir)4個(gè)波段的灰度信息；比值波段主要包括歸一化植被指數(shù)(NDVI)、增強(qiáng)型植被指數(shù)(EVI)、差值植被指數(shù)(DVI)、修正植被指數(shù)(mNDVI)、簡(jiǎn)單比值植被指數(shù)(SR)、GF42。比值波段構(gòu)建公式如表1所示。

表1 比值波段構(gòu)公式Table 1 Ratio band structure formula

本研究選用的紋理特征信息主要包括均值(Mean，ME)、方差(Variance，VAR)、協(xié)同性(Homogeneity，HOM)、對(duì)比度(Contrast，CON)、相異性(Dissimilarity，DIS)、熵(Entropy，ENT)、二階矩(Angular Second Moment ，ASM)和相關(guān)性(Correlation，COR)8種信息?，F(xiàn)有研究表明，不同窗口大小的紋理信息對(duì)郁閉度估測(cè)的影響不同，而且某些紋理信息之間存在顯著的相關(guān)性，不同紋理信息之間存在大量冗余信息[1,3]。本研究首先對(duì)GF-1影像進(jìn)行降維處理，提取其第1主成分；然后利用ENVI5.3SP1提取第1主成分圖像在不同窗口大小下的紋理圖像，窗口大小分別設(shè)置為3×3、5×5、…、41×41，方向設(shè)置為135°，步長(zhǎng)設(shè)置為1，灰度級(jí)壓縮至64級(jí)。分別以提取的3×3窗口至41×41窗口下的8個(gè)紋理信息為自變量，以樣地郁閉度為因變量，利用SPSS分別進(jìn)行多元逐步回歸分析，所得模型精度如圖3所示。

通過(guò)模型的R、R2以及RMSE來(lái)衡量不同窗口的紋理信息對(duì)研究區(qū)森林郁閉度定量估測(cè)的貢獻(xiàn)[22]。分析圖3可知，當(dāng)窗口大小為3×3時(shí)，R和R2均最大，RMSE最小，該窗口大小的紋理信息對(duì)森林郁閉度估測(cè)有較大的貢獻(xiàn)。

圖3 基于不同窗口紋理信息的逐步回歸分析模型精度Fig.3 Model accuracy of stepwise regression analysis based on texture information of different windows

同一窗口下不同紋理特征之間存在顯著的相關(guān)性(表2)。對(duì)3×3窗口下的各紋理特征信息以及各紋理特征信息與樣地郁閉度進(jìn)行相關(guān)性分析，不難發(fā)現(xiàn)，對(duì)比度與方差，相異性與協(xié)同性、對(duì)比度，對(duì)比度與熵，二階矩與熵之間的相關(guān)性均>0.9，存在顯著的相關(guān)性，綜合考慮各紋理特征信息之間以及各紋理特征信息與郁閉度之間的相關(guān)性，最終選取均值、方差、協(xié)同性、二階矩以及相關(guān)性共5個(gè)紋理特征作為自變量因子，在最大限度保留有效信息的基礎(chǔ)上減少冗余信息對(duì)建模的影響。

表2 相關(guān)性分析Table 2 Correlation analysis

綜上所述，本研究從遙感、地形以及紋理特征3個(gè)方面共得到15個(gè)原始特征，但過(guò)多的變量會(huì)造成數(shù)據(jù)冗余，使計(jì)算量大為增加，并且影響郁閉度估測(cè)模型的預(yù)報(bào)能力。本研究采用平均殘差平方和準(zhǔn)則和窮舉算法進(jìn)行建模變量?jī)?yōu)選[16,18,23]。平均殘差平方和表達(dá)式為

(9)

式中:q為變量的個(gè)數(shù)，Y為因變量矩陣，X為自變量矩陣。

3.2 模型構(gòu)建

為研究自變量因子設(shè)置不同對(duì)森林郁閉度估測(cè)的影響規(guī)律，本研究以遙感信息為基礎(chǔ)，通過(guò)對(duì)4種方案進(jìn)行建模變量?jī)?yōu)選，最終確定分3種方案設(shè)置可能影響郁閉度估測(cè)的自變量因子。方案1：遙感信息；方案2：遙感和地形信息；方案3：遙感、地形和紋理信息。對(duì)于每種方案，穩(wěn)健估計(jì)和偏最小二乘法均經(jīng)過(guò)變量?jī)?yōu)選獲取最終參與建模的自變量因子，k-NN法使用方案內(nèi)的所有自變量因子進(jìn)行建模。每種方案均進(jìn)行35次隨機(jī)抽樣建模，綜合分析所有建模結(jié)果從而減少隨機(jī)抽樣給建模帶來(lái)的偶然性影響。各方案建模效果如表3所示。

表3 3種方案建模綜合效果比較Table 3 Comparison of the comprehensive effects of three schemes of modeling %

分析表3，橫向比較可知，3種建模方案中k-NN法的模型精度標(biāo)準(zhǔn)差最小，故該方法建模效果是最好的；而相較于k-NN法而言，穩(wěn)健估計(jì)和偏最小二乘法的標(biāo)準(zhǔn)差均較大并且兩者建模效果相近，建模效果容易受到隨機(jī)抽樣的偶然性影響。縱向比較可知，在遙感信息的基礎(chǔ)上增加地形信息或紋理信息作為自變量因子能夠有效提升參數(shù)化方法的建模效果，并且標(biāo)準(zhǔn)差減小，建模穩(wěn)定性增加，但k-NN法的建模效果前后沒有發(fā)生明顯的變化。

為研究3種建模方法在森林郁閉度定量估測(cè)中的實(shí)際應(yīng)用效果，在方案3中本研究共對(duì)樣地進(jìn)行35次隨機(jī)抽樣，抽樣的原則按照模型殘差相對(duì)均方根誤差與預(yù)報(bào)偏差相對(duì)均方根誤差相近且兩者都在1倍標(biāo)準(zhǔn)差范圍內(nèi)，而后對(duì)模型結(jié)果進(jìn)行篩選，最終分別挑選出3種方法中的最具代表性的郁閉度估測(cè)模型。各模型精度如表4所示，參數(shù)化方法的最優(yōu)變量為：紅光波段(R)、方差(VAR)、二階矩(ASM)、EVI、GF42、mNDVI、NDVI、SR、坡度(S)。

表4 各模型精度Table 4 Accuracy of each model (%)

由圖4可見，對(duì)郁閉度估測(cè)模型中的k取不同值時(shí)發(fā)現(xiàn)，基于k-NN法的郁閉度估測(cè)模型的最佳k值為24。

圖4 郁閉度估測(cè)模型K取不同值時(shí)均方根誤差Fig.4 The root mean square error when the canopy closure estimation model K takes different values

采用穩(wěn)健估計(jì)法所建模型的估測(cè)方程如下

采用偏最小二乘法所建模型的估測(cè)方程如下

3.3 估測(cè)精度分析

利用GF-1PMS影像、坡度和坡向等地形信息以及紋理信息，獲取各樣地大小單元對(duì)應(yīng)自變量的值，建立以樣地大小為單位的郁閉度估測(cè)模型，最終獲得各單元對(duì)應(yīng)的郁閉度[13-14]。最后利用均值法對(duì)郁閉度估測(cè)值進(jìn)行統(tǒng)計(jì)，獲取小班郁閉度。分別以47個(gè)現(xiàn)地實(shí)測(cè)小班數(shù)據(jù)和隨機(jī)抽取的227個(gè)二類調(diào)查小班數(shù)據(jù)(約占研究區(qū)有林地小班的5%)為標(biāo)準(zhǔn)評(píng)價(jià)各郁閉度估測(cè)模型的實(shí)際應(yīng)用效果。最終3種郁閉度估測(cè)模型的實(shí)際應(yīng)用效果如表5所示。

表5 各郁閉度估測(cè)模型的實(shí)際應(yīng)用效果Table 5 The actual application effect of each canopy closure estimation model %

依據(jù)表5可知，非參數(shù)化估測(cè)方法的實(shí)際應(yīng)用效果優(yōu)于參數(shù)化估測(cè)方法，整體而言3種郁閉度估測(cè)方法的應(yīng)用效果均能滿足森林資源二類調(diào)查的實(shí)際應(yīng)用需求，但k-NN法和穩(wěn)健估計(jì)的實(shí)際應(yīng)用效果明顯優(yōu)于偏最小二乘法。

4 結(jié)論與討論

森林郁閉度定量估測(cè)是森林遙感估測(cè)研究中的一項(xiàng)主要內(nèi)容[24]，傳統(tǒng)的方法未對(duì)不同樣地?cái)?shù)據(jù)適用的模型進(jìn)行分析，未提取出針對(duì)不同樣地?cái)?shù)據(jù)最優(yōu)的模型，估測(cè)精度難以提高。本研究對(duì)傳統(tǒng)方法中存在的問(wèn)題，運(yùn)用GF-1衛(wèi)星遙感影像提取遙感信息以及紋理信息，結(jié)合樣地?cái)?shù)據(jù)和地形信息，建立以樣地大小為單位的森林郁閉度估測(cè)模型，最終反演獲得研究區(qū)森林郁閉度，并且模型估測(cè)精度有很大提高。

穩(wěn)健估計(jì)、偏最小二乘法和k-NN法3種郁閉度估測(cè)方法的應(yīng)用效果均能滿足實(shí)際需求，但k-NN法無(wú)論是以現(xiàn)地實(shí)測(cè)小班數(shù)據(jù)為標(biāo)準(zhǔn)還是以二類調(diào)查小班數(shù)據(jù)為標(biāo)準(zhǔn)，估測(cè)精度均達(dá)到85%以上，具有較好的估測(cè)效果，明顯優(yōu)于其他2種參數(shù)化方法。這說(shuō)明GF-1 PMS 多光譜影像在森林郁閉度定量估測(cè)方面具有一定的潛力。

2種參數(shù)化方法的模型精度相近，但實(shí)際郁閉度估測(cè)效果卻存在較大的差別，穩(wěn)健估計(jì)的實(shí)際估測(cè)效果明顯優(yōu)于偏最小二乘法。

關(guān)于估測(cè)穩(wěn)定性評(píng)估方面，k-NN法在森林郁閉度遙感估測(cè)中穩(wěn)定性較好，而穩(wěn)健估計(jì)法和偏最小二乘法估測(cè)模型不穩(wěn)定，估測(cè)效果波動(dòng)性較大。

在遙感信息的基礎(chǔ)上增加地形信息或紋理信息作為自變量因子能夠有效提升參數(shù)化方法的建模效果，增加建模穩(wěn)定性，但對(duì)于k-NN法并沒有明顯的影響。遙感信息完全可以滿足利用k-NN法進(jìn)行森林郁閉度估測(cè)的需求，因此在進(jìn)行郁閉度估測(cè)時(shí)，應(yīng)根據(jù)選用的方法選擇合適的建模變量因子組合。