基于k-NN非參數(shù)模型的高山松生物量遙感估測研究

謝福明，舒清態(tài)，字 李，吳 榮，吳秋菊，汪 紅，劉 延，吉一濤

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

在應(yīng)對全球氣候變化背景下，森林碳匯的相關(guān)研究成為科學(xué)界關(guān)注的熱點。具體研究內(nèi)容涉及森林固碳能力與計量、監(jiān)測方法與技術(shù)研究等方面。其中，森林碳匯計量和監(jiān)測是當(dāng)前研究的核心內(nèi)容之一。我國森林生物量/碳儲量研究剛開始不久，存在基礎(chǔ)資料不夠系統(tǒng)完整和估算方法不統(tǒng)一等問題，估算的結(jié)果差異較大，可比性較差，估算精度還需進一步提高；由于遙感數(shù)據(jù)具有很大的優(yōu)勢，在大尺度的森林生物量估測中，遙感估測法成為大面積預(yù)估森林生物量最主要的方法[1]。傳統(tǒng)的光學(xué)遙感技術(shù)通過處理分析衛(wèi)星影像的像元亮度值，結(jié)合地面調(diào)查，間接估測較大區(qū)域的連續(xù)森林生物量分布。目前，基于以線性、非線性回歸為主的參數(shù)模型和以神經(jīng)網(wǎng)絡(luò)(NN)、支持向量機(SVM)、最近鄰法為代表的非參數(shù)化模型的估測方法被普遍使用。針對大多數(shù)基于遙感數(shù)據(jù)或其指數(shù)與實測生物量數(shù)據(jù)之間的相關(guān)性建立經(jīng)驗?zāi)Ｐ瓦M行區(qū)域尺度的空間反演中，非參數(shù)化方法比較有優(yōu)勢[1]。其中k-最近鄰法(k-NN)以其可用于單變量和多變量預(yù)測、不需要假設(shè)自變量與因變量的分布、容易把外部信息運用在感興趣的地理區(qū)域和適用于各種數(shù)據(jù)集等優(yōu)勢使它特別受歡迎[2-3]。1951年，F(xiàn)ix和Hodges首次將最近鄰法作為一種非參數(shù)識別技術(shù)用于對分布未知的種群進行分類[4]。1990年，Tomppo[5]在芬蘭國家森林資源調(diào)查中提及了k-最近鄰法的應(yīng)用，為該技術(shù)在自然資源領(lǐng)域的應(yīng)用研究開創(chuàng)了里程碑。此后，諸多國家的許多研究者使用k-最近鄰回歸模型對生物量、蓄積量、碳儲量、樹高、胸高斷面積等森林參數(shù)進行估測，并且在采用遺傳算法(GA)加權(quán)優(yōu)化特征變量間的相似度評估后，該模型的精度得到了提升[6-10]。結(jié)合Landsat影像、MODIS衛(wèi)星數(shù)據(jù)以及機載激光掃描數(shù)據(jù)的k-NN模型在大面積尺度的森林變量估測中都相繼取得了較好的結(jié)果[10-15]。國內(nèi)的部分研究學(xué)者也將k-最近鄰法運用于區(qū)域森林生理參數(shù)估計，如陳爾學(xué)等[16]運用Landsat數(shù)據(jù)和k-NN法對小面積統(tǒng)計單元森林蓄積量估測，其結(jié)果表明采用k-NN法對縣市級統(tǒng)計單元森林參數(shù)的估測效果明顯優(yōu)于只利用固定樣地數(shù)據(jù)的傳統(tǒng)參數(shù)估測方法。2011年，郭穎[17]運用k-最近鄰法(k-NN)等非參數(shù)回歸模型對甘肅省西水林場的森林地上生物量進行估測，并用隨機森林算法(RF)進行特征選擇后估測精度得以提升，優(yōu)化后的算法在處理錯誤樣本時具有良好的容錯能力。然而，針對不同的森林類型或樹種、不同規(guī)模尺度的遙感數(shù)據(jù)與樣地數(shù)據(jù)以及不同的模型等，其估測效果都有所不同，而且沒有可比性。因此，本文結(jié)合Landsat數(shù)據(jù)與非參數(shù)k-NN模型實現(xiàn)了香格里拉市的高山松地上生物量估測反演，獲取高山松生物量空間分布格局，為低緯度高海拔地區(qū)的喬木森林資源保護及開發(fā)提供了技術(shù)案例。

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

1.1 研究區(qū)概況

圖1 研究區(qū)地理位置與樣地點分布圖Fig.1 Location of study area and the distribution of field data

研究區(qū)香格里拉市(99°20′～100°19′E、26°52′～28°52′N)隸屬迪慶藏族自治州，位于云南省西北部、迪慶州東北部(圖1)。境內(nèi)地勢高聳，熱量不足，氣溫偏低，多年平均氣溫5.5 ℃,歷年平均降雨量618.4 mm，屬山地寒溫帶季風(fēng)氣候。由于地處云南亞熱帶常綠闊葉林植被區(qū)向青藏高原高寒植被區(qū)過渡地帶，植被分布南北差異明顯，山體垂直分布完整而典型。森林面積大，覆蓋率高，天然林占絕對優(yōu)勢，主要分布有10種植被類型，常見的樹種有云杉(Piceaasperata)、冷杉(Abiesfabri)、高山松(Pinusdensata)、云南松(Pinusyunnanensis)、高山櫟(Quercussemicarpifolia)等。其中，高山松是研究區(qū)境內(nèi)喜光、耐寒、耐貧瘠的主要優(yōu)勢樹種，適應(yīng)性廣，更新能力強，一般分布在海拔3 300～3 800 m，占全市喬木林面積的22.7%。

1.2 數(shù)據(jù)采集

1.2.1 遙感數(shù)據(jù) Landsat-8衛(wèi)星于2013年2月11日發(fā)射，是美國陸地探測衛(wèi)星系列的后續(xù)衛(wèi)星。Landsat-8衛(wèi)星搭載了陸地成像儀(OLI)和熱紅外傳感器(TIRS)。OLI被動感應(yīng)地表反射的太陽輻射和散發(fā)的熱輻射，有9個波段的感應(yīng)器，覆蓋了從紅外到可見光的不同波長范圍。與Landsat-7衛(wèi)星的ETM+傳感器相比，OLI增加了一個藍(lán)色波段(0.433～0.453 μm)和一個短紅外波段(band 9：1.360～1.390 μm)，藍(lán)色波段主要用于海岸帶觀測，短紅外波段具有水汽強吸收特征；可用于云檢測的TIRS是有史以來最先進，性能最好的熱紅外傳感器。本文采用OLI傳感器3景影像數(shù)據(jù)完全覆蓋香格里拉市區(qū)域(見圖1(b))，衛(wèi)星數(shù)據(jù)參數(shù)詳見表1。

表1 香格里拉研究區(qū)衛(wèi)星數(shù)據(jù)參數(shù)表

1.2.2 樣地數(shù)據(jù) 本研究中樣地數(shù)據(jù)采集時間為2014年10月24日—11月14日，在香格里拉市境內(nèi)的高山松分布范圍內(nèi)，布設(shè)了50塊不同齡級的30 m×30 m樣地(圖1)，使模型的建立更符合現(xiàn)實意義。主要記錄樣地內(nèi)每株立木的徑階、樹高、胸徑以及樣地坐標(biāo)和海拔等因子。樣地中心GPS坐標(biāo)是使用差分GPS(DGPS)測得樣地4個角點坐標(biāo)后計算得到。由于DGPS精度相對手持GPS精度較高，其相對誤差均小于1 m，遠(yuǎn)小于遙感因子的30 m分辨率，所以本次樣地采樣在很大程度上減少了一般手持GPS誤差較大而帶來的樣地位置和真實位置之間的差異。

2 數(shù)據(jù)處理與研究方法

2.1 數(shù)據(jù)處理

2.1.1 遙感數(shù)據(jù)預(yù)處理 遙感影像前期預(yù)處理在圖像像元上最大化地還原地物的原始反射率等諸多信息，對于構(gòu)建的反演模型至關(guān)重要，包括輻射定標(biāo)、大氣校正、幾何精校正、鑲嵌與裁剪等工作。利用ENVI 5.3實現(xiàn)對影像數(shù)據(jù)的預(yù)處理，輻射定標(biāo)工具(Radiometric Calibration)能自動從元數(shù)據(jù)文件中讀取參數(shù)，將圖像的數(shù)字化量值(DN)轉(zhuǎn)化為輻射亮度值或反射率等物理參數(shù)；FLASSH大氣校正模塊是基于MODTRANS輻射傳輸模型的較先進的校正算法，能夠排除氣溶膠、地形等因素造成的反射信息誤差；幾何精校正則基于SPOT高分辨率影像，利用地面控制點GCP和多項式幾何校正數(shù)學(xué)模型，實現(xiàn)對非系統(tǒng)因素產(chǎn)生的誤差的校正；最后，利用矢量邊界完整的裁剪出研究區(qū)域的影像，見圖1。

2.1.2 樣地數(shù)據(jù)處理 利用高山松立木材積公式計算每株立木的材積，求和得樣地高山松每木檢尺蓄積量，樣地實測數(shù)據(jù)包含：平均樹高、平均胸徑、樣地蓄積和中心點地理坐標(biāo)；其中，蓄積量和生物量轉(zhuǎn)換模型則引用于前人的研究結(jié)果,見式(1),實測生物量描述見表2，該模型基于統(tǒng)計學(xué)原理與方法，實現(xiàn)了林分蓄積量與生物量之間良好的線性轉(zhuǎn)換關(guān)系。最后，將樣地數(shù)據(jù)依據(jù)地理坐標(biāo)疊加至遙感數(shù)據(jù)圖層，提取樣地點處對應(yīng)的像元值作為建模因子。

B=0.516 8×V+33.237 8[18]

(1)

式(1)中：B為高山松的林分生物量，單位：t/hm2；V為林分蓄積量，單位：m3/hm2。

表2 高山松實測樣地地上生物量描述統(tǒng)計分析

2.1.3 建模因子提取 對預(yù)處理后的影像數(shù)據(jù)通過波段運算工具提取單波段、多波段組合、植被指數(shù)和地形特征等共計18個因子，作為建模因子備選參數(shù)：

(1)原始波段因子。一般研究植被的波長范圍為400～2 500 nm，故本次采用的原始波段為可見光波段(400～700 nm)：Band 1-4、近紅外波段(700～1 300 nm)：Band 5和短紅外波段(1 300～2 500 nm)：Band6-7。

(2)波段指數(shù)因子。波段指數(shù)的信息來源于原始波段，卻以2個或多個波長范圍內(nèi)的地物反射率的組合運算增強了地物某一特性或者細(xì)節(jié)。本次研究涉及的波段指數(shù)因子及其計算公式如下：

歸一化植被指數(shù)：NDVI=(ρNIR-ρR)/(ρNIR+ρR)；

(2)

比值植被指數(shù)：SR=(ρNIR/ρR)；

(3)

差值植被指數(shù)：DVI=ρNIR-ρR；

(4)

(5)

(6)

(7)

(8)

其中，ρNIR、ρR、ρBLUE分別為近紅外波段、紅外波段、藍(lán)光波段的反射率；Bandi為Landsat OLI 傳感器的第i波段；L為土壤調(diào)節(jié)系數(shù)，因香格里拉市植被覆蓋面積為89%，本文中L取0.25[19]。

(3)地形因子。由DEM派生的坡度(slope)、坡向(aspect)和海拔(elevation)因子組成。

2.2 k-最近鄰法(k-NN)

(9)

(10)

dij2=(Xi-Xj)′M(Xi-Xj)

(11)

2.2.2 精度評價 運用樣本置留法，按3∶2的比例隨機選取30個樣本作為參考集，20個樣本作為目標(biāo)集。并采用決定系數(shù)(R2)、均方根誤差(RMSE)對模型精度進行評價，二者計算公式為：

(12)

(13)

式中，n為樣本數(shù)，xi為目標(biāo)集響應(yīng)變量的模型估測值，yi為目標(biāo)集響應(yīng)變量的實測值。

3 結(jié)果與分析

本研究中的k-NN模型算法在MATLAB軟件平臺編程實現(xiàn)，主要算法流程為數(shù)據(jù)的讀取與歸一化、計算特征變量間的空間相似度、依據(jù)k值選取最近鄰個體、確定近鄰加權(quán)大小和計算變量的預(yù)測值等。同時，考慮到模型反演數(shù)據(jù)量大，采用分塊處理的方法將研究區(qū)進行分塊后逐一估測，實現(xiàn)對整個研究區(qū)像元尺度上的高山松生物量估測反演。

3.1 最優(yōu)參數(shù)的選擇

3.1.1k值的優(yōu)化選擇k值決定了計算響應(yīng)變量預(yù)測值時需要考慮的最近鄰個數(shù)。本次研究中考慮了0～20范圍內(nèi)的k值，在不同的模型結(jié)構(gòu)配置下采用了精度評價指標(biāo)分析選取最佳值：當(dāng)k<5時，估測模型的精度隨k值呈現(xiàn)增加趨勢，在k=1時模型的精度最低；當(dāng)k≥10時，模型的精度趨向平穩(wěn)，受k值的影響較?。辉趉值為4或5時，模型精度達(dá)到最佳，均方根誤差(RMSE)達(dá)到最低值。圖2是不同參數(shù)組合下的精度評價指標(biāo)隨k值波動的變化曲線，直觀地反映了k-NN模型參數(shù)對其預(yù)測結(jié)果的影響情況。

圖2 不同參數(shù)組合下的k-NN模型精度變化曲線Fig.2 Accuracy curve of the k-NN model with different parameters

表3 不同參數(shù)組合下的k-NN模型精度對比(部分)

圖3 高山松實測生物量與預(yù)測生物量的擬合線圖(模型參數(shù)：歐式距離、k=4、t=2)Fig.3 Observations versus predictions:Pinus densata aboveground biomass(Model metrics:Euclidean distance,k=4,t=2)

綜合上述參數(shù)的擇優(yōu)分析，確定本次研究最終采用的k-NN模型參數(shù)結(jié)構(gòu)為：采用歐式距離度量特征變量間的相似度，k、t值分別取4和2；圖3是基于參數(shù)優(yōu)化的k-NN模型的生物量預(yù)測值與實測值的線性擬合結(jié)果。

3.2 生物量估測反演

目前，結(jié)合遙感數(shù)據(jù)和實測樣地數(shù)據(jù)，利用回歸模型對生物量等森林生理參數(shù)和屬性估測的研究中一般采用基于像元尺度或基于樣地尺度的2種估測方法。本文基于像元尺度的模型估算研究區(qū)高山松生物量并用預(yù)測的生物量值進行反演。圖4a是依據(jù)二類調(diào)查數(shù)據(jù)將高山松分布區(qū)域隨機分成27塊的編號圖，利用已構(gòu)建的模型依次對每一塊區(qū)域的生物量進行估測，其對應(yīng)的像元數(shù)與生物量估測結(jié)果見表4。統(tǒng)計得到整個研究區(qū)的高山松生物量估計值為22 496 140.79(約0.22億t)。估測的高山松生物量在空間分布上與研究區(qū)實際情況相一致，與岳彩榮[20]基于SVM的高山松生物量估測結(jié)果(約為0.17億t)和王金亮等[21]利用遙感信息模型估測的結(jié)果(約為0.20億t)等有可對比性。

(a)分塊編號；(b)生物量空間分布圖(a)Block number;(b) The spatial distribution of biomass圖4 像元尺度下基于k-NN模型的香格里拉市高山松地上生物量反演結(jié)果Fig.4 The spatial distribution of Pinus densata aboveground biomass in Shangri-la at the pixel scale based on k-NN model

編號Number像元個數(shù)Number of pixel生物量/tAboveground biomass編號Number像元個數(shù)Number of pixel生物量/tAboveground biomass編號Number像元個數(shù)Number of pixel生物量/tAboveground biomass197 0611 113 706.251069 531807 940.921957 616671 981.94227 263304 072.2011169 4791 958 356.972049 767569 445.943102 3601 173 724.961287 8811 016 597.732165 713780 240.74448 066550 318.6313103 0101 203 704.102266 890779 716.65567 590778 537.251475 249866 274.592350 723589 941.746111 4831 269 431.551524 887285 493.502429 258339 609.41748 328558 960.921653 199615 392.442566 588774 198.578117 7881 363 803.001750 877594 906.542629 323341 042.579165 1111 900 327.331870 132792 035.852743 692496 379.45

像元個數(shù)總計為1 939 779；生物量總計為22 496 140.79(t)

研究區(qū)的高山松生物量預(yù)測值最低為62 t/hm2，最高為157 t/hm2，主要集中在125～135 t/hm2，其次是115～125 t/hm2。數(shù)值介于實際生物量的測量值(49～195 t/hm2)，正如Chirici[3]所說，k-NN法的一個缺陷就是估測的值只能局限于實測值范圍內(nèi)進行估測，不會對超出實測值范圍的值進行預(yù)測，但也能客觀的反映出研究區(qū)高山松生物量的分布格局(圖4b)。

4 討論與結(jié)論

4.1 討論

本研究基于非參數(shù)模型建立了高山松地上生物量遙感估測模型，并實現(xiàn)了研究區(qū)香格里拉市的高山松生物量反演：以Landsat8/OLI為信息源，結(jié)合地面50個實測樣地數(shù)據(jù)，建立了研究區(qū)高山松地上生物量k-最近鄰法(k-NN)遙感估測模型。對k-NN模型參數(shù)的優(yōu)化結(jié)果表明：(1)特征變量間的相似度評價方式中，歐式距離最穩(wěn)定，效果最好，其次是馬氏距離和絕對值距離；(2)最佳的近鄰個數(shù)k值為4，其對模型精度的影響呈現(xiàn)規(guī)律：當(dāng)k<5時，模型的精度隨k值呈現(xiàn)增加趨勢；當(dāng)k≥10時，模型的精度趨向平穩(wěn)，受k值的影響較小；在k值為4或5時，模型精度最高，k=1時模型的精度最低；(3)模型精度受距離分解因子t的影響較小，相比之下，t值取2時模型精度較好。參數(shù)優(yōu)化后的k-NN模型擬合精度為，R2等于0.71，RMSE為18.21 t/hm2，估測得到研究區(qū)像元尺度下的高山松生物量為0.22億t。

4.2 結(jié)論

(1)在建模因子選取多少的問題上，應(yīng)采取優(yōu)化算法實現(xiàn)對特征變量的加權(quán)或排除不相關(guān)變量，并通過與不同的傳感器以及非遙感變量(如降水或氣候圖)等數(shù)據(jù)的集成使用，來提高對森林結(jié)構(gòu)和功能等森林屬性的預(yù)測精度。

(2)k-NN模型的結(jié)構(gòu)配置優(yōu)化可以通過進行多種診斷測試，觀察參數(shù)和特征變量的變化來確定更準(zhǔn)確的結(jié)果，比如包括特征變量的篩選、選取復(fù)雜的距離度量方式，以及k和t的優(yōu)化。目前，沒有統(tǒng)一標(biāo)準(zhǔn)的優(yōu)化步驟順序使模型產(chǎn)生更高的精度，一個合理的順序應(yīng)為：首先選擇通過消除冗余特征變量來減少特征變量的維度災(zāi)難影響；然后根據(jù)特征變量集的性質(zhì)，選擇更復(fù)雜的距離度量(例如加權(quán)歐氏距離或典型對應(yīng)分析指標(biāo))進行測試；最后可以利用交叉驗證的留一法對k值進行優(yōu)化選擇。

(3)雖然遙感數(shù)據(jù)能夠較好的估算森林生物量，但普遍存在生物量飽和問題，即隨著生物量的增大，光譜反射率和后射散射系數(shù)對生物量的敏感性較小，不能有效的反應(yīng)生物量的變化情況，降低了高森林生物量的估測精度，故研究光譜“飽和點”的識別并在其基礎(chǔ)上的生物量遙感估算模型的構(gòu)建就變得更有意義。