周友鋒,謝秉樓,李明詩

(1. 南京林業(yè)大學(xué)林草學(xué)院、水土保持學(xué)院,南方現(xiàn)代林業(yè)協(xié)同創(chuàng)新中心,江蘇 南京 210037;2. 浙江省森林資源監(jiān)測(cè)中心,浙江 杭州 310020)

森林是陸地生態(tài)系統(tǒng)的重要組成部分,是陸地上最大的碳庫,儲(chǔ)存了大約 80% 的地表碳和 40% 的地下碳,在全球碳循環(huán)中扮演重要角色[1-2]。森林地上生物量(aboveground biomass,AGB)是森林生態(tài)系統(tǒng)固碳能力的直觀表達(dá),也是評(píng)估森林生態(tài)系統(tǒng)碳收支平衡的重要指標(biāo)[3]。作為森林生態(tài)系統(tǒng)碳匯潛力評(píng)估的重要要素,如何快速、準(zhǔn)確地獲取大尺度 AGB 信息對(duì)于宏觀掌握森林碳儲(chǔ)量及其分布,進(jìn)一步制定公平合理的碳排放政策具有重要意義[4-5]。

隨機(jī)森林算法(random forest,RF)作為一種優(yōu)秀的機(jī)器學(xué)習(xí)算法,近年來被廣泛應(yīng)用于基于遙感數(shù)據(jù)的 AGB 制圖研究中。RF 所構(gòu)建的模型屬于非參數(shù)模型,能夠應(yīng)對(duì) AGB 與遙感因子之間復(fù)雜的非線性關(guān)系[6-7]。而且它對(duì)訓(xùn)練樣本中存在的噪聲敏感度低,能較好應(yīng)對(duì)由于數(shù)據(jù)缺失所引起的精度降低問題,同時(shí)還能識(shí)別預(yù)測(cè)變量的重要性[8-9]。諸多的研究表明,RF 相較于其他機(jī)器學(xué)習(xí)算法和傳統(tǒng)統(tǒng)計(jì)回歸方法具有更高的預(yù)測(cè)精度[10-12]。而 RF 在預(yù)測(cè) AGB 時(shí)只考慮了 AGB 與遙感因子間的關(guān)系,而忽略 AGB 制圖時(shí)鄰近觀測(cè)數(shù)據(jù)的空間自相關(guān)性[13-14]。將隨機(jī)森林與克里金法進(jìn)行協(xié)同從而構(gòu)建隨機(jī)森林/克里金方法框架(random forest Kriging,RFK)將能夠有效應(yīng)對(duì)上述缺陷。此框架通過克里金插值法對(duì)RF預(yù)測(cè)的模型殘差值進(jìn)行建模,以分離殘差項(xiàng)中的結(jié)構(gòu)化成分(空間自相關(guān)的描述項(xiàng)),并將之疊加到隨機(jī)森林模型的預(yù)測(cè)結(jié)果上,從而達(dá)到改進(jìn)制圖精度的目的[15]。在過去近10年內(nèi),RFK模型已經(jīng)用于預(yù)測(cè)環(huán)境因子、土壤有機(jī)質(zhì)、樹木材積等研究中,其預(yù)測(cè)精度均優(yōu)于 RF 模型[13,16-17]。而在近兩年的 AGB 預(yù)測(cè)研究中也出現(xiàn)了這類方法,研究結(jié)果皆表明,RFK 模型在 AGB 預(yù)測(cè)中是可行的。Chen 等[18]和 Silveira等[19]分別在大興安嶺地區(qū)和大西洋沿岸熱帶山地森林區(qū)域,采用 RFK 模型對(duì) AGB 進(jìn)行預(yù)測(cè),所得到的預(yù)測(cè)精度及各項(xiàng)誤差指標(biāo)均優(yōu)于 RF 模型。

然而,上述研究所使用的克里金插值法均為普通克里金插值法(ordinary Kriging,OK),它是最常用于區(qū)域化變量的最優(yōu)無偏插值,但在環(huán)境復(fù)雜的地區(qū)預(yù)測(cè)能力有限[20]。協(xié)同克里金插值法(Co-Kriging,CK)是 OK 的一種擴(kuò)展方法,通過添加一個(gè)或多個(gè)協(xié)變量,能夠考慮多個(gè)變量之間的相關(guān)關(guān)系,因而可提高空間插值的精度效果[21-22]。目前已有的研究多為基于多源遙感數(shù)據(jù),利用隨機(jī)森林普通克里金(random forest ordinary Kriging,RFOK)模型對(duì) AGB 進(jìn)行預(yù)測(cè),很少有利用隨機(jī)森林協(xié)同克里金(random forest Co-Kriging,RFCK)模型進(jìn)行AGB 制圖研究的案例。因此,本研究采用國家連續(xù)清查野外樣地調(diào)查數(shù)據(jù),Landsat 5 TM 數(shù)據(jù)、ALOS-1 PALSAR-1 數(shù)據(jù)和 STRM DEM 數(shù)據(jù),利用 RFCK 模型來執(zhí)行 2012 年廣東省北部亞熱帶森林 AGB 的制圖任務(wù)。期望本研究提出的方法和制圖結(jié)果,可以為我國南方森林增匯和可持續(xù)森林管理實(shí)踐發(fā)展提供技術(shù)參考。

1 材料與方法

1.1 研究區(qū)概況

研究區(qū)位于廣東省北部(112.74°～115.09°E,23.59°～25.51°N),主要包括韶關(guān)、清遠(yuǎn)和河源等市(圖 1)。研究區(qū)的氣候?qū)儆趤啛釒駶櫦撅L(fēng)氣候,年平均降水量為 1 300～2 400 mm,年平均氣溫達(dá)18～21 ℃ ,水熱條件優(yōu)良,適合植被生長。海拔 10～1 710 m,地形起伏大,以山地和丘陵為主。研究區(qū)內(nèi)森林資源豐富,植被群落類型復(fù)雜多樣,涵蓋了亞熱帶地區(qū)典型的植被群落類型。研究區(qū)屬于典型的人工林區(qū)域,經(jīng)常性的森林采伐和更新事件廣泛存在[23],林分年齡時(shí)空異質(zhì)性較高,便于 AGB 估算模型的構(gòu)建與性能評(píng)價(jià)。樹種多為常綠速生樹種,優(yōu)勢(shì)樹種包括杉木(Cunninghamialanceolata)、毛竹(Phyllostachysedulis)、馬尾松(Pinusmassoniana)、桉樹(Eucalyptusrobusta)等。

底圖審圖號(hào):GS(2019)3333。圖1 研究區(qū)樣地點(diǎn)分布STRM 影像(上)及研究區(qū)縮略圖(下)Fig. 1 The field sampling sites STRM image (upper) and overview map (lower)

1.2 數(shù)據(jù)來源及預(yù)處理

1.2.1 遙感數(shù)據(jù)

本研究所使用的遙感數(shù)據(jù)主要包括2011 的 Landsat 5 TM 數(shù)據(jù)、2010 年的 ALOS-1 PALSAR-1 數(shù)據(jù)和 30 m 分辨率的SRTM DEM 數(shù)據(jù),以上數(shù)據(jù)的投影坐標(biāo)系均為WGS 1984 UTM Zone 49N。其中,Landsat 5 TM 和 SRTM DEM 數(shù)據(jù)來源于 USGS EROS 數(shù)據(jù)中心(https://glovis.usgs.gov/);ALOS-1 PALSAR-1 數(shù)據(jù)來源于日本宇宙航空研究開發(fā)機(jī)構(gòu)(JAXA)。

Landsat 5 TM自 2011 年 11 月起停止采集影像,于是下載了鄰近年份的 Landsat 5 TM 影像,并處于植被生長旺季內(nèi)(采集時(shí)間為 2011 年 8 月 20 日)且含云量最小。并對(duì)之進(jìn)行幾何校正、輻射定標(biāo)和大氣校正以及地形校正等預(yù)處理操作。

由于 Landsat 5 TM 僅能采集森林冠層的平面信息,在高郁閉度森林內(nèi)易發(fā)生光譜飽和現(xiàn)象[24]。因此本研究引入 ALOS-1 PALSAR-1 數(shù)據(jù),其長波信號(hào)對(duì)森林冠層具有一定的穿透力,能在一定程度上降低飽和效應(yīng)[25]。ALOS-1 PALSAR-1 數(shù)據(jù)自 2011 年 4 月起停止采集,故本研究采用鄰近年份 2010 年的影像數(shù)據(jù),數(shù)據(jù)產(chǎn)品等級(jí)為 Level 1.1,包括 HH 和 HV 雙極化信息,空間分辨率為 25 m。從影像中提取振幅數(shù)據(jù),利用公式(1)將其轉(zhuǎn)換為后向散射系數(shù)(σ0):

(1)

式中:ND為 HH 和 HV 后向散射體的振幅數(shù)值;FC為絕對(duì)校準(zhǔn)因子,為 -83 dB[26]。然后,對(duì)轉(zhuǎn)換為后向散射系數(shù)的 HH 和 HV 影像執(zhí)行 7×7 窗口的 Lee 濾波操作以減少散斑噪聲。最后將之重采樣為 30 m 分辨率,以匹配 Landsat 5 TM 數(shù)據(jù)。

1.2.2 樣地?cái)?shù)據(jù)

本研究所使用的樣地?cái)?shù)據(jù)為 2012 年的國家森林連續(xù)清查數(shù)據(jù),來源于廣東省森林資源監(jiān)測(cè)中心。預(yù)處理過程如下:首先利用廣東省森林資源監(jiān)測(cè)中心發(fā)布的該地區(qū)各樹種的異速生長方程[27](表 1),計(jì)算出連續(xù)清查數(shù)據(jù)的樣地 AGB 觀測(cè)值,單位為t/hm2;之后,對(duì)原始的 290 個(gè)樣地 AGB 觀測(cè)值進(jìn)行質(zhì)量控制。

表1 研究區(qū)內(nèi)優(yōu)勢(shì)樹種的異速生長方程

將遙感影像中被云及陰影所覆蓋的樣地點(diǎn)排除,對(duì)速生樹種(如桉樹)的樣地 AGB 通過生物量隨林齡增長模型進(jìn)行修正[28],再利用3倍標(biāo)準(zhǔn)差法剔除離群值,共剩下 245 個(gè)有效樣點(diǎn)。最后,對(duì)剩余的有效 AGB 觀測(cè)值依數(shù)值大小進(jìn)行分層抽樣,提取 80% 作為訓(xùn)練數(shù)據(jù),剩余 20% 作為驗(yàn)證數(shù)據(jù)。245 個(gè)有效樣點(diǎn)依據(jù)林分類型劃分的 AGB 分布狀況如表 2 所示。

表2 有效樣點(diǎn) AGB 分布狀況

1.3 特征變量提取

對(duì)于 Landsat 5 TM 數(shù)據(jù),通過對(duì)影像進(jìn)行多種光譜特征變換、紋理信息提取以及波段組合運(yùn)算,提取了共 80 個(gè)特征變量,包括植被指數(shù)、纓帽變換指數(shù)、主成分、紋理測(cè)度、原始單波段和波段組合。植被指數(shù)包含歸一化植被指數(shù)(NDVI),差值植被指數(shù)(DVI),比值植被指數(shù)(RVI),大氣阻力植被指數(shù)(ARVI),增強(qiáng)植被指數(shù)(EVI),土壤調(diào)節(jié)植被指數(shù)(SAVI);纓帽變換指數(shù)包含亮度(brightness)、綠度(greenness)與濕度(wetness);主成分為包含了 95% 以上原始影像光譜信息的第1主成分(PC1)、第2主成分(PC2)與第3主成分(PC3);紋理信息則通過基于灰度共生矩陣的紋理測(cè)度進(jìn)行提取,其在刻畫森林空間分布形態(tài)上是有效且重要的[29]。對(duì)第1主成分(PC1)采用 8 個(gè)紋理測(cè)度,包含均值(mean)、方差(variance)、均勻性(homogeneity)、對(duì)比度(contrast)、相異性(dissimilarity)、熵(entropy)、二階性(second moment)、相關(guān)性(correlation)。紋理測(cè)度提取時(shí)的滯后距離為1個(gè)像元,方向?yàn)橛蚁?移動(dòng)窗口大小為 3×3,5×5,7×7,9×9;原始波段包括 B1—B7(不包含熱紅外波段)6個(gè)波段,同時(shí)將對(duì)應(yīng)波段的反射率倒數(shù)作為備選變量,即 TM1_1,TM2_1,……,TM7_1;波段組合是對(duì)原始影像信息的組合,不同的波段組合會(huì)凸顯不同的影像特征從而豐富影像信息,且部分波段組合對(duì) AGB 相關(guān)性較高[30]。本研究提取 Landsat 5 TM 影像每個(gè)原始波段的地表反射率與其他 5 個(gè)波段反射率的比值,如 TM 75,即 B7/B5,共 30 個(gè)作為波段組合變量。

從 PALSAR-1 數(shù)據(jù)中提取的變量包括 HH 和 HV 的后向散射系數(shù)、HH/HV 和雷達(dá)森林退化指數(shù)(RFDI)[31]。此外,通過同樣的方式提取 HH 和 HV 極化數(shù)據(jù) 4 個(gè)移動(dòng)窗口大小(3×3,5×5,7×7,9×9)的8 個(gè)紋理測(cè)度。

另外,考慮到地形因子對(duì)植被生長有所影響,故本研究也從 STRM DEM 數(shù)據(jù)中提取坡度(slope)、坡向(aspect)、粗糙度(roughness)等地形因子作為建模備選變量。

1.4 隨機(jī)森林克里金建模

1.4.1 隨機(jī)森林建模

RF是一種基于決策樹的分類和回歸算法,通過多次 bootstrap 抽樣獲得多個(gè)隨機(jī)樣本,并通過這些樣本分別建立相對(duì)應(yīng)的決策樹,從而構(gòu)成隨機(jī)森林。該方法適用于解決分類和回歸問題,對(duì)于回歸問題,取所有決策樹預(yù)測(cè)結(jié)果的均值作為最終預(yù)測(cè)結(jié)果。本研究使用 R 語言中的 ‘randomForest’ 軟件包來實(shí)現(xiàn) RF 的建模過程。RF 有兩個(gè)重要參數(shù),分別為表示輸入變量數(shù)量的mtry和代表決策樹數(shù)量的ntree。mtry默認(rèn)為數(shù)據(jù)集中變量數(shù)的二次方根(分類模型)或1/3(回歸模型);ntree值是通過在模型誤差相對(duì)穩(wěn)定的情況下,經(jīng)過不斷測(cè)試能獲得多少個(gè)決策樹來確定的。本研究 RF建模所設(shè)定的參數(shù)值為:mtry= 3,ntree= 500。

本研究所提取的特征變量較多,需要對(duì)變量進(jìn)行篩選,選取與AGB相關(guān)性較高的變量進(jìn)行建模。使用randomForest包中的important命令進(jìn)行變量重要性分析,通過2個(gè)指標(biāo)均方誤差百分比增加量(%IncMSE)與節(jié)點(diǎn)純度增加量(IncNodePurity)來評(píng)估每個(gè)變量對(duì)建模性能的貢獻(xiàn)。%IncMSE與IncNodePurity值越大,表明對(duì)應(yīng)的預(yù)測(cè)變量的重要性越強(qiáng)[23]。為了保證各變量對(duì) RF 模型的綜合代表性,同時(shí)降低模型計(jì)算的復(fù)雜度,以 IncNodePurity 中前 25% 的變量為基礎(chǔ),若這些變量在 IncMSE% 中排名前 10,則選取作為建模變量,進(jìn)行 AGB 預(yù)測(cè)。

1.4.2 克里金插值法

克里金插值法是一種用于空間插值的地統(tǒng)計(jì)學(xué)方法,可用估計(jì)的預(yù)測(cè)誤差來評(píng)估預(yù)測(cè)的質(zhì)量,所輸入的數(shù)據(jù)集需要滿足正態(tài)分布假設(shè)。本研究采用的克里金插值法包括 OK 和CK 兩種方法。OK 以變異函數(shù)理論和結(jié)構(gòu)分析理論為基礎(chǔ),它通過基于區(qū)域化變量的變異函數(shù)生成最優(yōu)無偏估計(jì)[32]。計(jì)算如下:

(2)

式中:ROK(x0)為通過 OK 得到的殘差預(yù)測(cè)值,n是用于插值的樣點(diǎn)數(shù)量,Wi是點(diǎn)i的加權(quán)系數(shù),可根據(jù)最優(yōu)無偏估計(jì)原理與拉格朗日最小化原則確定[33],R(xi)是樣地點(diǎn)i的殘差值。

CK是 OK 的一種擴(kuò)展方法,通過添加一個(gè)或多個(gè)協(xié)變量,考慮了多個(gè)變量之間的相關(guān)關(guān)系,可提高結(jié)構(gòu)化成分的預(yù)測(cè)精度。由于研究區(qū)位于廣東北部山區(qū),AGB 空間分布受地形因素影響較大[19]。因此,本研究選擇高程作為協(xié)變量。CK 的插值公式如式(3)所示:

(3)

式中:RCK(x0)為通過 CK 得到的殘差預(yù)測(cè)值,R1(x1i)為樣地點(diǎn)i的殘差值,W1i為樣地點(diǎn)i殘差權(quán)重,R2(x2j)為樣地j的高程,W2j為樣地j的高程權(quán)重,N1為訓(xùn)練樣本個(gè)數(shù),N2為高程樣本點(diǎn)個(gè)數(shù),其中N1≥N2。

克里金法用變異函數(shù)測(cè)定空間自相關(guān)要素。變異函數(shù)描述的是區(qū)域化變量空間變化的特征和強(qiáng)度,可表示為隨著距離增加,兩樣點(diǎn)間半變異函數(shù)值或協(xié)方差函數(shù)值的變化情況。變異函數(shù)可供擬合模型較多,本研究采用地統(tǒng)計(jì)學(xué)軟件 GS+ 進(jìn)行變異函數(shù)擬合模擬,選擇最優(yōu)擬合模型。GS+ 所提供的擬合模型有指數(shù)函數(shù)(exponential)、球面函數(shù)(spherical)和高斯函數(shù)(Gaussian),所采用的評(píng)價(jià)指標(biāo)為決定系數(shù)(R2)與殘差平方和(RSS)。擬合模型的R2越大,RSS 越小,擬合性能越好。變異函數(shù)的3個(gè)模型參數(shù)是塊金(nugget)、變程(range)和基臺(tái)(sill)。塊金是距離為 0 時(shí)的變異函數(shù)值,表示測(cè)量或分析誤差;變程是變異函數(shù)值穩(wěn)定時(shí)的距離,即與空間自相關(guān)距離相對(duì)應(yīng);基臺(tái)是變異函數(shù)的最大值。塊金效應(yīng)是塊金值與基臺(tái)值的比值,能夠描述空間自相關(guān)性的強(qiáng)弱,塊金效應(yīng)越小,空間自相關(guān)性越強(qiáng)[34]。

1.4.3 隨機(jī)森林克里金模型

RFK模型的實(shí)現(xiàn)步驟分為2步:①通過 RF 建模,得到 AGB 預(yù)測(cè)值;②通過克里金插值法分離殘差中的結(jié)構(gòu)化成分,并將之疊加到隨機(jī)森林模型預(yù)測(cè)值上[式(4)]。RF 模型預(yù)測(cè)殘差值通過公式(5)計(jì)算得到。

R(xi)=BBF(xi)-B(xi);

(4)

BRFOK/RFCK(xi)=BRF(xi)-ROK/CK(Xi)。

(5)

式中:R(xi)是樣地點(diǎn)i的殘差值,B(xi)是樣地點(diǎn)i的AGB觀測(cè)值,BRF(xi)是基于 RF 模型的樣地點(diǎn)i的AGB預(yù)測(cè)值,BRFOK/RFCK(xi)是通過 RFOK 或 RFCK 模型得到的 AGB 預(yù)測(cè)值,ROK/CK(xi)是樣地點(diǎn)i通過 OK 或 CK 得到的殘差預(yù)測(cè)值。

1.5 模型精度評(píng)價(jià)

使用決定系數(shù)(R2)[式(6)]、平均絕對(duì)誤差[MAE,式中記為σ(MAE)][式(7)]、均方根誤差[RMSE,式中記為σ(RMSE)][式(8)]指標(biāo)量化模型的性能。采用式(9)計(jì)算 RFOK 模型和 RFCK 模型相對(duì)于 RF 模型的相對(duì)改進(jìn)指數(shù)(relative improvement)[RI,式中記為σ(RI)],評(píng)估模型的改進(jìn)效果。公式如下:

(6)

(7)

(8)

(9)

2 結(jié)果與分析

2.1 隨機(jī)森林模型構(gòu)建

構(gòu)建 RF 模型首先要求優(yōu)選特征變量,為RF模型中變量重要度排序結(jié)果見圖2,選擇了如下 10 個(gè)建模參數(shù):HV、HVcorrelation99、mean99、TM75、TM57、TM53、TM35、TM21、TM12、B2?？傮w而言,Landsat 5 TM 數(shù)據(jù)的原始波段和波段組合占了建模變量的很大比例。而紋理信息變量(例如:HVcorrelation99 和 mean99)重要度排序較高,對(duì) AGB 預(yù)測(cè)也有一定影響。

以上述 10 個(gè)建模參數(shù)構(gòu)建 RF 模型,同時(shí)通過訓(xùn)練集進(jìn)行模型訓(xùn)練?；谟?xùn)練集預(yù)測(cè) AGB 和實(shí)測(cè) AGB 間的R2為 0.96,MAE 為15.77 t/hm2,RMSE 為 19.98 t/hm2,表明模型的擬合效果優(yōu)良。

B1、B2、B3和B7表示Landsat 5 TM影像1號(hào)、2號(hào)、3號(hào)與7號(hào)波段的地表反射率;TM12、TM13、TM15、TM21、TM24、TM31、TM34、TM35、TM42、TM51、TM52、TM53、TM57、TM74、TM75表示Landsat 5 TM影像某一波段地表反射率與另一波段地表反射率的比值,如TM12,即B1/B2;TM1_1、TM2_1、TM3_1和TM7_1表示Landsat 5 TM影像1號(hào)、2號(hào)、3號(hào)與7號(hào)波段地表反射率的倒數(shù);PC1和PC3表示Landsat 5 TM影像的第1主成分與第3主成分;RVI和ARVI分別表示比值植被指數(shù)與大氣阻力植被指數(shù);Brightness表示Landsat 5 TM影像通過纓帽變換得到的亮度值;mean77、mean99和correlation99分別表示基于PC1采用7×7窗口得到的均值紋理特征,采用9×9窗口得到的均值與相關(guān)性紋理特征;HH和HV表示PALSAR-1 HH與HV極化信息的后向散射系數(shù);HHcorrelation55和HHcorrelation77分別表示HH采用5×5與7×7窗口得到的相關(guān)性紋理特征;HVmean55、HVmean77和HVmean99分別表示HV采用5×5、7×7與9×9窗口得到的均值紋理特征;HVcorrelation77和HVcorrelation99分別表示HV采用7×7與9×9窗口得到的相關(guān)性紋理特征。B1, B2, B3 and B7 represent the surface reflectance of Landsat 5 TM images at bands 1, 2, 3 and 7. TM12, TM13, TM15, TM21, TM24, TM31, TM34, TM35, TM42, TM51, TM52, TM53, TM57, TM74 and TM75 represent the ratio of surface reflectance in one band of Landsat 5 TM image to that in another band, such as TM12, that is B1/B2; TM1_1, TM2_1, TM3_1 and TM7_1 represent the reciprocal of surface reflectance in bands 1, 2, 3 and 7 of Landsat 5 TM images. PC1 and PC3 represent the first principal component and the third principal component of Landsat 5 TM images. RVI and ARVI represent ratio vegetation index and atmospheric resistance vegetation index, respectively. Brightness indicates the brightness value of Landsat 5 TM images obtained by hat transformation. mean77, mean99 and correlation99 indicate the mean texture features obtained in 7×7 window size based on PC1, and the mean value and correlation texture features obtained in 9×9 window size, respectively. HH and HV represent the backscattering coefficients of HH and HV polarization information of PALSAR-1. HHcorrelation55 and HHcorrelation77 indicate the correlation texture features obtained in the 5×5 and 7×7 window sizes of HH, respectively. HVmean55, HVmean77 and HVmean99 represent the mean texture features obtained by HV using 5×5, 7×7 and 9×9 window sizes, respectively. HVcorrelation77 and HVcorrelation99 indicate the correlation texture features obtained by HV using 7×7 and 9×9 window sizes, respectively. 圖2 RF 中變量重要度排序結(jié)果Fig. 2 The importance ranking of the variables for AGB mapping by using RF model

2.2 隨機(jī)森林預(yù)測(cè)殘差分析

對(duì)隨機(jī)森林預(yù)測(cè)殘差進(jìn)行統(tǒng)計(jì)分析,得到殘差的范圍為-59.73～39.64 t/hm2,均值為 0.53 t/hm2,標(biāo)準(zhǔn)差為20.04 t/hm2,絕對(duì)峰度值為3.50(接近 3),絕對(duì)偏度值為0.89(接近 1),說明殘差近似服從正態(tài)分布。因此,符合進(jìn)行克里金插值的前提假設(shè),可以對(duì) AGB 預(yù)測(cè)殘差進(jìn)行隨后的克里金插值分析。

OK和 CK 的變異函數(shù)模型模擬結(jié)果及相關(guān)參數(shù)見圖3及表3。依據(jù) GS+ 建模的結(jié)果,在 OK 和 CK 方法下,高斯函數(shù)模型擬合效果略優(yōu)于其他兩個(gè)模型,最終作為 OK 和 CK 最優(yōu)擬合變異函數(shù)模型?？傮w而言,CK 的變異函數(shù)模型擬合性能優(yōu)于 OK 的模型,前者的R2更大,RSS 更小(表 3)。圖4是依式(2)和式(3)建立的殘差空間插值結(jié)果。CK 較 OK 具有更大的殘差預(yù)測(cè)范圍,OK 的殘差插值范圍為-19.61～19.38 t/hm2,CK 的范圍則為-28.72～21.37 t/hm2。此外,本研究還通過布設(shè) 100 個(gè)隨機(jī)點(diǎn)提取得到不同海拔區(qū)間的平均插值結(jié)果。OK 和 CK 在低海拔區(qū)域(34～300 m)的平均插值結(jié)果分別為 1.24 和 2.09 t/hm2,在中海拔區(qū)域(301～600 m)為 0.13 和 0.30 t/hm2,而在高海拔區(qū)域(601～1 158 m)則為-0.45 和-2.01 t/hm2。可見,在3個(gè)海拔區(qū)域中,通過分析 AGB 高程空間位置關(guān)系,CK 相比 OK 均分離了更多殘差項(xiàng)中的結(jié)構(gòu)化成分,所得到的殘差插值效果更好。

表3 OK 和 CK 的變異函數(shù)擬合模型及其參數(shù)

圖3 基于 OK 和CK的殘差變異函數(shù)模擬Fig. 3 The variogram simulation of residuals derived from OK and CK models

圖4基于OK和CK的殘差插值結(jié)果和各模型生成的AGB專題圖Fig. 4 The inter polated spatial patterns of residuals derived from OK, CK models and AGB patterns derived from the models

2.3 RFOK和RFCK模型精度分析

依據(jù)式(5),從RFOK模型和RFCK模型構(gòu)建改進(jìn)后的AGB預(yù)測(cè)值(圖5)。如圖5所示,在模型的泛化能力方面,RFCK模型的AGB預(yù)測(cè)范圍為0.48～165.76 t/hm2,略高于RF和RFOK模型。3個(gè)模型的AGB預(yù)測(cè)值在高海拔的山地區(qū)域都較大,且呈現(xiàn)由西北往東南增多的趨勢(shì)。用20%的獨(dú)立樣本進(jìn)行驗(yàn)證,結(jié)果如下:RFCK模型的RI值為0.08,R2由0.46增加到0.57,MAE由27.28減少到25.12 t/hm2,RMSE由32.48減少到29.80 t/hm2。RFOK模型的RI值為0.03,R2由0.46增加到0.51,MAE由27.28減少到26.63 t/hm2,RMSE由32.48減少到31.58 t/hm2。兩種改進(jìn)模型的精度評(píng)價(jià)表現(xiàn)均優(yōu)于RF模型,且RFCK模型較RFOK模型更勝一籌。圖5展示了RF模型的驗(yàn)證數(shù)據(jù)集R2在0.46左右,擬合線(虛線)與1∶1線有較大差異,存在比較明顯的低值被高估,高值被低估的現(xiàn)象,但是這種現(xiàn)象隨著針對(duì)預(yù)測(cè)誤差的空間插值技術(shù)的引入而得到一定程度的減弱。

圖5 模型驗(yàn)證時(shí) AGB 觀測(cè)值與預(yù)測(cè)值的散點(diǎn)圖Fig. 5 Scatterplots of the observed AGB and the predicted AGB when validating the models

3 討 論

根據(jù)重要度排序結(jié)果,本研究選擇了 10 個(gè)變量進(jìn)行 RF 建模。其中,與兩個(gè)短波紅外波段(B5、B7)反射率的波段組合占據(jù)很大比例。Landsat 5 TM 的B5與B7波段的反射率與植被水分含量有關(guān),在夏季多雨的亞熱帶地區(qū)分辨植被較為有利,而通過設(shè)置與其他波段的比值更能夠使原始波段上不易區(qū)分的植被能夠識(shí)別[35]。HVcorrelation99和mean99與AGB有較高相關(guān)性,這與部分研究的特征優(yōu)選結(jié)果是一致的[36],表明在林分結(jié)構(gòu)復(fù)雜的亞熱帶森林區(qū)域引入紋理信息進(jìn)行AGB估測(cè)是可靠的。另外,由PALSAR-1數(shù)據(jù)發(fā)展的特征變量(后向散射系數(shù)HV等)在AGB估測(cè)上也有貢獻(xiàn),這主要由于PALSAR-1數(shù)據(jù)的長波信號(hào)能夠穿透冠層,獲取植被AGB的主體——樹干的信息。

本研究的結(jié)果表明,結(jié)合了殘差插值結(jié)果的RFOK與RFCK模型較RF模型具有更高的AGB制圖精度,低值高估和高值低估的現(xiàn)象在一定程度上得到改善,各項(xiàng)精度指標(biāo)均優(yōu)于RF模型,同時(shí)也優(yōu)于部分亞熱帶地區(qū)AGB估測(cè)研究結(jié)果[11,21]。但本研究不足在于RFOK模型與RFCK模型的RI分別為0.03和0.08,精度改進(jìn)效果不高。以往有關(guān)RFK模型在AGB預(yù)測(cè)應(yīng)用的研究,最高RI值均能夠達(dá)到0.10以上,本研究結(jié)果與之相比偏低[18-19],這與本研究AGB殘差空間自相關(guān)性較低有關(guān)。本研究OK和CK變異函數(shù)模型的塊金效應(yīng)值較高,均在0.9以上。而之前有關(guān)土壤特性的研究結(jié)果表明,當(dāng)塊金效應(yīng)值低于0.6[13],其精度提升效果十分顯著。這也表明相比土壤特性,AGB預(yù)測(cè)殘差受空間地理位置分布關(guān)系的影響并不大。而就數(shù)據(jù)源而言,一方面本研究采用的遙感數(shù)據(jù)包括Landsat 5 TM 和PALSAR-1數(shù)據(jù)。前者易存在光譜飽和現(xiàn)象,后者雖能減弱飽和現(xiàn)象,但飽和點(diǎn)通常在 150 t/hm2左右,在本研究研究區(qū)內(nèi)的高 AGB 區(qū)域作用相對(duì)有限[26]。另一方面可能與 RF 模型的建模變量結(jié)合了Landsat 5 TM 和 PALSAR-1 數(shù)據(jù),而非單一數(shù)據(jù)源有關(guān)[18]。此外,還可能與本研究遙感數(shù)據(jù)與樣地?cái)?shù)據(jù)時(shí)間不匹配有關(guān)。本研究所使用的 Landsat 5 TM 數(shù)據(jù)為 2011 年,PALSAR-1 數(shù)據(jù)為2010年,而樣地?cái)?shù)據(jù)則為 2012 年。研究區(qū)內(nèi)部分樣地屬于灌木林、幼齡林與跡地,且樹種多為速生樹種。雖然本研究對(duì)這些樣地進(jìn)行了部分修正,但這些樣地的遙感數(shù)據(jù)信息與實(shí)際樣地信息顯然仍有差距。這一方面部分解釋了 RF 模型的過擬合現(xiàn)象,另一方面可能對(duì)于殘差數(shù)據(jù)的空間自相關(guān)性有所影響。

RFCK模型通過將高程作為協(xié)變量,除考慮 AGB 距離與方位間的關(guān)系外,還將 AGB高程間的空間關(guān)系考慮在內(nèi),相比 RFOK 模型所得到的變異函數(shù)擬合效果更好,空間自相關(guān)性更強(qiáng),這與其他有關(guān) CK 的研究結(jié)果是一致的[20,22]。從結(jié)果上看,相較 RFOK 模型,RFCK 模型在本研究中高海拔區(qū)域(如西北、中部及東南區(qū)域)和低海拔地區(qū)(如東北、南部區(qū)域)的插值效果更好。相對(duì)應(yīng)地,RFCK 模型在刻畫 AGB 空間分布模式的效果時(shí)也要優(yōu)于 RFOK 模型,輪廓更加清晰。RFCK 模型在山地和平原的分布情況更加接近實(shí)際情況,特別是山地區(qū)域,這也說明飽和效應(yīng)得到進(jìn)一步削弱。因此在具有高 AGB 的亞熱帶山地區(qū)域,相比傳統(tǒng)的 RF 和隨后改進(jìn)的 RFOK 模型而言,RFCK 模型更加適合用來進(jìn)行亞熱帶山地區(qū)域森林 AGB 制圖,若結(jié)合長時(shí)間序列數(shù)據(jù)進(jìn)行分析,能夠更好掌握當(dāng)?shù)罔駱?、松樹等人工林與天然林的生長健康狀態(tài)與森林碳儲(chǔ)量動(dòng)態(tài)分布情況,以評(píng)價(jià)當(dāng)前森林經(jīng)營撫育措施是否合理,有助于進(jìn)一步制定針對(duì)性強(qiáng)的森林經(jīng)營方法和公平合理的碳排放政策,增強(qiáng)森林生態(tài)系統(tǒng)綜合效益。

雖然 RFCK 模型的 AGB 預(yù)測(cè)精度優(yōu)于 RFOK 模型,但是精度提升的效果并不明顯,這可能與本研究只采用了一個(gè)協(xié)變量有關(guān)。如果在選取協(xié)變量時(shí),僅僅選取一個(gè)與目標(biāo)變量相關(guān)性較高的因子,則不能全面反映協(xié)變量的信息,預(yù)測(cè)效果提升有限[20]。因此未來的研究可考慮添加多個(gè)協(xié)變量,以充分發(fā)揮 CK 的優(yōu)勢(shì)。同時(shí)可通過主成分分析法選取協(xié)變量,以應(yīng)對(duì) CK 在插值時(shí)協(xié)變量較多造成計(jì)算復(fù)雜度增加,而協(xié)變量較少引起插值精度降低的問題[22]。另外,RFCK 模型雖然能降低 RF 模型的過擬合現(xiàn)象,但是 RFCK 模型的過擬合依然存在。RF 模型出現(xiàn)過擬合,與訓(xùn)練數(shù)據(jù)中存在噪聲有關(guān),而本研究中的樣地?cái)?shù)據(jù)恰恰有此問題。因此,除了上述通過優(yōu)選協(xié)變量進(jìn)行改進(jìn),也可對(duì) RF 模型訓(xùn)練時(shí)的樣地?cái)?shù)據(jù)進(jìn)行優(yōu)化減少噪聲,抑或增大研究區(qū),擴(kuò)大數(shù)據(jù)量來削弱模型過擬合現(xiàn)象。

本研究以 Landsat 5 TM、PALSAR-1 遙感影像數(shù)據(jù)為主要數(shù)據(jù)源,并且結(jié)合國家森林連續(xù)清查數(shù)據(jù),比較了 RF、RFOK 與 RFCK 模型在 AGB 預(yù)測(cè)中的效果,同時(shí)比較了 RFOK 與 RFCK 模型相較于 RF 模型的預(yù)測(cè)精度改進(jìn)水平?？傮w而言,以高程作為協(xié)變量的 RFCK 模型對(duì) AGB 的預(yù)測(cè)效果更好,其次是 RFOK 和 RF 模型,且 RFCK 模型的精度改進(jìn)效果也優(yōu)于 RFOK 模型,能更好地描繪實(shí)際 AGB 的分布情況,特別是在山地區(qū)域。本研究所獲得的高精度亞熱帶大區(qū)域 AGB 專題圖,有助于了解當(dāng)?shù)?AGB 分布情況,為政府制定與完善相關(guān)政策提供依據(jù),并為早日實(shí)現(xiàn)碳達(dá)峰與碳中和目標(biāo)打下基礎(chǔ)。