吳 迪，李 冰，楊愛玲

(國家測繪地理信息局黑龍江基礎(chǔ)地理信息中心，黑龍江 哈爾濱 150086)

基于Landsat長時(shí)間序列數(shù)據(jù)估算樹高和生物量

吳 迪，李 冰，楊愛玲

(國家測繪地理信息局黑龍江基礎(chǔ)地理信息中心，黑龍江 哈爾濱 150086)

以Landsat長時(shí)間序列數(shù)據(jù)為研究對象，旨在以光譜序列信息反演森林參數(shù)為視角，應(yīng)用Landtrendr算法從時(shí)間序列數(shù)據(jù)中提取森林?jǐn)_動變量，使用隨機(jī)森林計(jì)算方法建立擾動變量、反射率和GLAS激光點(diǎn)森林參數(shù)之間的關(guān)系模型，獲取樹高和生物量的空間分布信息。為多源遙感數(shù)據(jù)反演森林參數(shù)提供參考，研究證明基于Landsat長時(shí)間序列數(shù)據(jù)獲得的森林?jǐn)_動變量能夠增強(qiáng)反射率和森林參數(shù)之間的相關(guān)性，可提高預(yù)測精度。

Landsat長時(shí)間序列數(shù)據(jù)；樹高；生物量；森林?jǐn)_動變量；Landtrendr算法

森林在生長過程中會受到各種擾動(火災(zāi)、病蟲害、砍伐)和非擾動(溫度、濕度和降水)因素的影響。導(dǎo)致林木中的碳以不同的形式流動于碳庫之中。森林?jǐn)_動信息與碳儲量密切相關(guān)，對林木參數(shù)具有很好的解釋能力。NASA陸地衛(wèi)星(Landsat)計(jì)劃自1972年以來共發(fā)射8顆衛(wèi)星，對地球進(jìn)行了長達(dá)40多年的觀測，積累了大量的對地觀測數(shù)據(jù)，它的高分辨率長時(shí)間序列信息為重構(gòu)過去幾十年的森林?jǐn)_動和修復(fù)歷史提供了獨(dú)一無二的資源。Li 等[1]使用VCT算法從Landsat時(shí)間序列中提取擾動信息，應(yīng)用回歸樹模型對密西西比地區(qū)森林樹高進(jìn)行了準(zhǔn)確的估算(R2=0.91，RMSE=1.97 m)。Pflugmacher等[2]應(yīng)用雷達(dá)數(shù)據(jù)和時(shí)間序列擾動因子準(zhǔn)確獲取美國藍(lán)山保護(hù)區(qū)森林生物量。大量的實(shí)驗(yàn)證明使用長時(shí)間序列數(shù)據(jù)估算森林樹高和生物量具有很大的潛力。

按計(jì)算方式可以將時(shí)間序列分析方法歸納為兩類：“求差異”和“求趨勢”，其中“求差異”是檢測時(shí)間序列中突變信號，從背景噪聲中分離微小變化[3-4]?！扒筅厔荨笔菑脑肼曋蟹蛛x出持續(xù)狀態(tài)較長的趨勢信號[5-6]，前者主要用于識別突發(fā)性的擾動事件如火災(zāi)、風(fēng)災(zāi)、砍伐等，后者用于識別持續(xù)性的擾動事件如病蟲害，兩種技術(shù)只適用于特定的擾動類型，在應(yīng)用上存在局限性。Kennedy[7-10]等提出了landtrendr(Landsat-based detection of Trends in Disturbance and Recovery)時(shí)間序列分析方法可以用于識別不同類型的擾動事件并且能夠提供擾動后的森林修復(fù)信息。

本研究嘗試使用Landtrendr方法獲取Landsat 長時(shí)間序列數(shù)據(jù)的森林?jǐn)_動信息，利用這些信息對森林參數(shù)(樹高和地上生物量)進(jìn)行解釋，協(xié)同GLAS波形數(shù)據(jù)建立森林參數(shù)空間擴(kuò)展模型，獲取森林樹高和生物量的空間分布信息。

1 研究區(qū)域和數(shù)據(jù)準(zhǔn)備

研究區(qū)域位于小興安嶺中心，地理坐標(biāo)為(46°53′34″～48°03′40″N，127°28′05″～129°22′05″E)，森林類型以針闊葉混交林為主，見圖1。

圖1 研究區(qū)位圖

選用L3D(2005年10—11月)和L3G(2006年10—11月)2期版本號為33的GLA01、GLA05和GLA14產(chǎn)品計(jì)算森林參數(shù)(見圖1)，結(jié)合野外實(shí)測數(shù)據(jù)對ICESat/GLAS波形進(jìn)行平滑濾波、分解擬合、噪聲估計(jì)、確定有效信息范圍、坡度校正處理，獲得激光點(diǎn)樹高和生物量。

使用1984～2006年33期條帶號為117/27的Landsat TM/ETM+影像計(jì)算森林?jǐn)_動變量，所有影像均為經(jīng)過LEDAPS系統(tǒng)預(yù)處理后的反射率數(shù)據(jù)，數(shù)據(jù)的獲取時(shí)間集中在森林生長期(7～8月)平均每年對應(yīng)兩景影像，為擾動參數(shù)計(jì)算提供充足的數(shù)據(jù)。對影像進(jìn)行纓帽變換、波段計(jì)算獲取TCA和TCD兩個(gè)擾動監(jiān)測指數(shù)為識別擾動變量做準(zhǔn)備(見表1)。

表1 擾動監(jiān)測指數(shù)

注：TCB，TCW，TCG為纓帽變換三變量

2 研究方法和變量計(jì)算

2.1 研究方法

使用Landtrendr算法從Landsat長時(shí)間序列數(shù)據(jù)中提取森林?jǐn)_動變量的計(jì)算過程包括提取光譜時(shí)間序列、移除噪聲、識別潛在頂點(diǎn)、擬合軌跡、簡化模型和確定最佳模型6個(gè)部分。

1)提取光譜時(shí)間序列：使用移動窗口求平均值的方法從擾動監(jiān)測指數(shù)影像中提取光譜信息，獲取每個(gè)像元的時(shí)間序列軌跡。

2)移除噪聲，根據(jù)相似度閾值迭代刪除軌跡中具有噪聲特征的突起，每次迭代刪除最差的突起直到滿足閾值要求。

3)識別潛在頂點(diǎn)：使用時(shí)間序列分割方法確定頂點(diǎn)，將軌跡上起止年份對應(yīng)的光譜值作為初始頂點(diǎn)進(jìn)行最小二乘回歸計(jì)算，再根據(jù)擬合結(jié)果對每一年的光譜值預(yù)測，將真實(shí)值與預(yù)測值相差最大的那個(gè)年份設(shè)置為新頂點(diǎn)。新頂點(diǎn)將時(shí)間序列軌跡分割成兩部分，再分別對這兩部分執(zhí)行回歸計(jì)算獲取均方誤差，針對均方誤差值相對較大的部分重復(fù)上面的操作尋找新頂點(diǎn)。反復(fù)進(jìn)行迭代運(yùn)算直到獲取的頂點(diǎn)和分段數(shù)量滿足控制文件中的相應(yīng)設(shè)置為止。為了防止過度擬合還需要根據(jù)角度準(zhǔn)則對分段結(jié)果做收斂判斷，計(jì)算各分段之間的角度，迭代剔除角度最淺的頂點(diǎn)和線段直到滿足閾值設(shè)定。

4)軌跡擬合：確定頂點(diǎn)后使用回歸和點(diǎn)對點(diǎn)兩種連接方式進(jìn)行時(shí)間序列擬合，從第一個(gè)頂點(diǎn)開始對比兩種連接方式的擬合結(jié)果，使用均方誤差小的連接方式擬合線段，第二條線段的起點(diǎn)是第一條線段的終點(diǎn)，使用兩種方法分別計(jì)算，對回歸連接法和點(diǎn)到點(diǎn)連接法的計(jì)算結(jié)果進(jìn)行比較擇優(yōu)擬合。對所有頂點(diǎn)依次計(jì)算，最后得到一組前后連接的線段代替原始時(shí)間序列軌跡。計(jì)算擬合結(jié)果p值，如果滿足控制文件中的閾值要求則進(jìn)行下一步計(jì)算，否則需要對軌跡進(jìn)行重新擬合。

5)簡化模型，上述步驟計(jì)算得到最復(fù)雜的分割模型，對模型進(jìn)行簡化。根據(jù)控制文件中的修復(fù)率閾值和均方誤差準(zhǔn)則迭代計(jì)算，每次移除最弱的頂點(diǎn)對剩下的頂點(diǎn)重新進(jìn)行軌跡擬合直到分段數(shù)等于1時(shí)停止。

6)確定最佳模型。簡化模型的過程中每刪除一個(gè)頂點(diǎn)就會重新擬合一次軌跡產(chǎn)生一個(gè)F統(tǒng)計(jì)的p值，對這些p值進(jìn)行對比，選取最優(yōu)的模型得到最終的時(shí)間序列重構(gòu)結(jié)果。

2.2 變量計(jì)算

使用Landtrendr算法描繪林地像元光譜變化趨勢，通過Landtrendr時(shí)間序列分割計(jì)算對光譜響應(yīng)做分段和線性擬合可以獲得一系列首尾連接的線段，根據(jù)線段的變化趨勢可以捕獲擾動年份和森林修復(fù)信息。用下降趨勢線表示擾動，上升趨勢線表示擾動后的修復(fù)，水平線表示未受到干擾(見圖2)。為了對軌跡進(jìn)行定量化描述，研究從擬合曲線上提取了一系列擾動參數(shù)(見表2)，以圖2為例DBYEA和DBVAL變量分別對應(yīng)圖中點(diǎn)A的年份和光譜值；DDUR變量對應(yīng)A、B兩點(diǎn)之間的時(shí)間差；DMAG對應(yīng)a值；RBYEAR和RBVAL變量分別對應(yīng)點(diǎn)B的年份和光譜值；RMAG對應(yīng)b值；RDUR對應(yīng)A、C兩點(diǎn)之間的時(shí)間差。每個(gè)像元的時(shí)間序列可能會出現(xiàn)多個(gè)擾動修復(fù)事件，只研究有林地類像元計(jì)算最大的擾動信息。通過計(jì)算可以得到4種類型的時(shí)間序列擬合軌跡，分別為擾動后修復(fù)、擾動后未修復(fù)、只修復(fù)和沒有受到擾動影響。對所有類型的軌跡按以下規(guī)則進(jìn)行參數(shù)設(shè)置。

1)如果一個(gè)像元的時(shí)間序列只包含修復(fù)信息，那么將DBYEA設(shè)置為時(shí)間序列的起始年；DBVAL設(shè)置為起始年對應(yīng)的光譜值；DDUR和DMAG定義為0。

2)如果一個(gè)像元的時(shí)間序列只包含擾動信息，那么將RBYEAR設(shè)置為時(shí)間序列截止年份；RBVAL設(shè)置為截止年份對應(yīng)的光譜值；RMAG和RDUR定義為0。

3)如果一個(gè)像元的時(shí)間序列不包含擾動和修復(fù)信息，那么將所有參數(shù)設(shè)置為0。

圖2 TCA擬合軌跡(像元坐標(biāo)：128.15486E，47.41262N)

表2 擾動修復(fù)參數(shù)

3 樹高和生物量模型

為了評價(jià)擾動參數(shù)估算能力，研究使用兩組預(yù)測變量建立森林參數(shù)模型(見表3)。其中一組預(yù)測變量是與GLAS數(shù)據(jù)獲取時(shí)間(L3G：2006年10～11月)最接近的一景Landsat光譜數(shù)據(jù)。另一組預(yù)測變量是在單時(shí)相光譜信息的基礎(chǔ)上添加了擾動和修復(fù)信息。

使用隨機(jī)森林計(jì)算方法(Random Forest，RF)建立森林參數(shù)(樹高、地上生物量)與預(yù)測變量(見表3)之間的關(guān)系模型。通過遍歷分析找到每個(gè)RF模型的最優(yōu)模型參數(shù)組合，根據(jù)訓(xùn)練集均方誤差MSEOOB確定反演模型。從表4可知增加擾動信息可以提高森林樹高和地上生物量模型的擬合精度。

表3 模型預(yù)測變量

表4 模型擬合驗(yàn)證

4 結(jié)果與討論

使用隨機(jī)森林反演模型對試驗(yàn)區(qū)的森林參數(shù)進(jìn)行估算，應(yīng)用實(shí)測樣地?cái)?shù)據(jù)對反演結(jié)果進(jìn)行檢驗(yàn)。圖3從上到下依次顯示了樣地?cái)?shù)據(jù)驗(yàn)證樹高和生物量的結(jié)果，通過對兩組變量的預(yù)測結(jié)果進(jìn)行對比可知，使用單時(shí)相光譜信息可以估算森林參數(shù)，但是精度不高，其中地上生物量預(yù)測值與實(shí)測值之間的相關(guān)性較差(R2=0.305),增加擾動信息可以增強(qiáng)光譜反射率和森林參數(shù)之間的相關(guān)性，提高森林樹高和地上生物量模型的預(yù)測精度(樹高的預(yù)測精度從82.72%提高到88.27%，地上生物量的預(yù)測精度從37.63%提高到69.39%)。由圖4可知研究區(qū)域生物量范圍介于0～150 t/ha之間，集中分布在50～80 t/ha范圍；研究區(qū)樹高范圍介于0～27 m之間，大部分樹木高于12 m。

圖3 模型預(yù)測值與實(shí)測值之間的相關(guān)性

圖4 研究區(qū)域模型反演結(jié)果(地理坐標(biāo)：46°53′34″～48°03′40″N，127°28′05″～129°22′05″E，時(shí)間：2006年7月)

5 結(jié)束語

使用時(shí)間軌跡分析方法(Landtrendr)對33期Landsat TCA和TCD影像進(jìn)行時(shí)間序列重構(gòu)提取森林?jǐn)_動參數(shù)，再分別使用單時(shí)相的Landsat光譜信息和單時(shí)相光譜信息加擾動參數(shù)兩組變量對森林參數(shù)進(jìn)行估算，對比結(jié)果顯示增加擾動信息可以增強(qiáng)光譜反射率和森林參數(shù)之間的相關(guān)性，提高樹高和生物量的預(yù)測精度，使用含有擾動參數(shù)的Landsat光譜信息對GLAS波形數(shù)據(jù)進(jìn)行空間擴(kuò)展可以獲得準(zhǔn)確的森林參數(shù)信息。

使用Landsat長時(shí)間序列數(shù)據(jù)可以準(zhǔn)確的估算森林參數(shù)，但這其中存在著很大的不確定性，因?yàn)槭苌謹(jǐn)_動歷史資料缺失的限制人們無法對擾動變量進(jìn)行準(zhǔn)確的驗(yàn)證，目前常用的檢驗(yàn)方法都是基于高分辨率影像進(jìn)行目視解譯和人工判讀，由于面向的對象是大量的時(shí)間序列數(shù)據(jù)，使擾動參數(shù)的驗(yàn)證工作充滿了挑戰(zhàn)。

[1] LI A, HUANG C, SUN G, et al. Modeling the height of young forests regenerating from recent disturbances in Mississippi using Landsat and ICESat data[J] Remote Sensing of Environment, 2011, 115(8): 1837-1849.

[2] PFLUGMACHER D, COHEN W B, KENNEDY R E. Using Landsat-derived disturbance history (1972～2010) to predict current forest structure[J].Remote Sensing of Environment, 2012, 122: 146-165.

[3] 鄧超, 程鋼, 孫震輝. 基于 TM 影像的焦作市城市擴(kuò)張與驅(qū)動力分析[J]. 測繪工程, 2015,24(8):57-64.[4] HUANG C, GOWARD S N, MASEK J G, et al. An automated approach for reconstructing recent forest disturbance history using dense Landsat time series stacks[J]Remote Sensing of Environment, 2010, 114(1): 183-198.

[5] Hostert P, R?DER A, HILL J. Coupling spectral unmixing and trend analysis for monitoring of long-term vegetation dynamics in Mediterranean rangelands[J]Remote sensing of environment, 2003, 87(2): 183-197.

[6] VOGELMANN J E, KOST J R, TOLK B, et al. Monitoring landscape change for Landfire using multi-temporal satellite imagery and ancillary data[J] Selected Topics in Applied Earth Observations and Remote Sensing, 2011, 4(2): 252-264.

[7] 王海君，許捍衛(wèi)，金文韜.基于Landsat-8遙感影像的鄱陽湖水邊線提取方法研究[J].測繪與空間地理信息，2016,39(1)：144-146.

[8] 劉利.國際地理信息經(jīng)濟(jì)效益定量估算初步研究[J].測繪與空間地理信息，2016,39(5)：4-7.

[9] 張飛，王建強(qiáng)，羅寒.基于時(shí)間序列的鄱陽湖水位分析與預(yù)測[J].測繪與空間地理信息，2016,39(8)：35-37.

[10] KENNEDY R E, BRAATEN J, YANG Z Q, et al. Land Trendr Users Guide[C]The Laboratory for Applications of Remote Sensing in Ecology (LARSE),2011.

[責(zé)任編輯：李銘娜]

Estimation of tree height and biomass based on long time series data of landsat

WU Di, LI Bing, YANG Ailing

(Heilongjiang Geomatics Center,China Bureau of Surveying and Maping, Harbin 150086,China)

This paper, taking the long time series data of Landsat as the research object, aims at retrieving forest parameters from spectral sequence information. The variables of forest disturbance are extracted from the time series data using the Landtrendr algorithm. With the random forest method to establish the model of the relationship among the disturbance variables, the reflectivity and the GLAS laser point, the spatial distribution information of the tree height and biomass is obtained, which provides a reference for retrieving forest parameters from multi source remote sensing data. The research proves that the forest disturbance variables obtained from the Landsat long time series data can enhance the correlation between the reflectance and forest parameters to improve the prediction accuracy.

long time series data of Landsat；tree height；biomass；forest disturbance variables；Landtrendr algorithm

引用著錄：吳迪，李冰，楊愛玲.基于Landsat長時(shí)間序列數(shù)據(jù)估算樹高和生物量[J].測繪工程，2017，26(6)：1-5.

10.19349/j.cnki.issn1006-7949.2017.06.001

2016-03-28

地理空間信息工程國家測繪地理信息局重點(diǎn)實(shí)驗(yàn)室贊助項(xiàng)目

吳 迪(1987-),女，工程師,博士.

P237

A

1006-7949(2017)06-0001-05