吳 迪，李 冰，楊愛(ài)玲

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

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

吳 迪，李 冰，楊愛(ài)玲

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

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

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

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

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

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

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

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

圖1 研究區(qū)位圖

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

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

表1 擾動(dòng)監(jiān)測(cè)指數(shù)

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

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

2.1 研究方法

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

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

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

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

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

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

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

2.2 變量計(jì)算

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

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

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

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

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

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

3 樹(shù)高和生物量模型

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

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

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

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

4 結(jié)果與討論

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

圖3 模型預(yù)測(cè)值與實(shí)測(cè)值之間的相關(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é)束語(yǔ)

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

使用Landsat長(zhǎng)時(shí)間序列數(shù)據(jù)可以準(zhǔn)確的估算森林參數(shù)，但這其中存在著很大的不確定性，因?yàn)槭苌謹(jǐn)_動(dòng)歷史資料缺失的限制人們無(wú)法對(duì)擾動(dòng)變量進(jìn)行準(zhǔn)確的驗(yàn)證，目前常用的檢驗(yàn)方法都是基于高分辨率影像進(jìn)行目視解譯和人工判讀，由于面向的對(duì)象是大量的時(shí)間序列數(shù)據(jù)，使擾動(dòng)參數(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ū)動(dòng)力分析[J]. 測(cè)繪工程, 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遙感影像的鄱陽(yáng)湖水邊線提取方法研究[J].測(cè)繪與空間地理信息，2016,39(1)：144-146.

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

[9] 張飛，王建強(qiáng)，羅寒.基于時(shí)間序列的鄱陽(yáng)湖水位分析與預(yù)測(cè)[J].測(cè)繪與空間地理信息，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

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

10.19349/j.cnki.issn1006-7949.2017.06.001

2016-03-28

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

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

P237

A

1006-7949(2017)06-0001-05