鐘 莉 陳蕓芝 汪小欽

(福州大學(xué)空間數(shù)據(jù)挖掘與信息共享教育部重點(diǎn)實(shí)驗(yàn)室 衛(wèi)星空間信息技術(shù)綜合應(yīng)用國(guó)家地方聯(lián)合工程研究中心數(shù)字中國(guó)研究院(福建) 福州 350108)

森林是陸地生態(tài)系統(tǒng)的重要組成部分，是地球上最大的碳儲(chǔ)存庫(kù)，其動(dòng)態(tài)變化直接影響全球陸地森林生態(tài)系統(tǒng)的碳循環(huán)和碳蓄積(葛全勝等， 2008； 趙金龍等， 2013)。森林動(dòng)態(tài)變化包括森林突發(fā)性變化和植被時(shí)序趨勢(shì)變化。森林干擾是森林突發(fā)性變化，根據(jù)干擾的起因，通常可分為人為干擾和自然干擾，其形式主要表現(xiàn)為森林火災(zāi)、森林采伐以及因受雨雪和病蟲(chóng)災(zāi)害等所引起的森林改變(Gowardetal., 2013)。近年來(lái)，隨著經(jīng)濟(jì)發(fā)展，全球各地森林植被頻繁遭到嚴(yán)重破壞，人為干擾(包括森林采伐、城市化、改林為耕)也逐漸增加(李維長(zhǎng)， 2000)。

遙感技術(shù)的不斷進(jìn)步，衛(wèi)星遙感數(shù)據(jù)作為時(shí)間尺度變化監(jiān)測(cè)的重要手段越來(lái)越受到研究者重視(Coppinetal.， 1996)。森林干擾研究的數(shù)據(jù)源通常分為基于區(qū)域和全球尺度遙感數(shù)據(jù)2類(lèi)，其中，AVHRR、MODIS和SPOT VEG等數(shù)據(jù)空間分辨率較低，多用于大尺度的森林變化監(jiān)測(cè)(宋富強(qiáng)等， 2011； 王強(qiáng)等， 2012)； 相比而言，Landsat系列影像數(shù)據(jù)具有較高精度，多用于區(qū)域尺度的森林干擾監(jiān)測(cè)及其造成的森林碳通量變化研究(楊辰等， 2015)。美國(guó)地質(zhì)調(diào)查局于2008年免費(fèi)開(kāi)放了Landsat堆棧數(shù)據(jù)(Woodcocketal.， 2008)，該數(shù)據(jù)具有良好的穩(wěn)定性和持續(xù)性，基于Landsat時(shí)序數(shù)據(jù)的研究取得了長(zhǎng)足發(fā)展(Townshendetal.， 2012； Wulderetal.， 2015)。以往研究中，森林變化監(jiān)測(cè)通常是對(duì)影像分類(lèi)后進(jìn)行兩兩比較(Cohenetal.， 2010)，不僅費(fèi)時(shí)費(fèi)力導(dǎo)致效率低下，而且針對(duì)長(zhǎng)時(shí)序分析精度不高，不能滿足應(yīng)用需求(楊辰等， 2013)；影像差異法利用差值、比值、典型相關(guān)性變化等方法構(gòu)建差異影像，雖然消除了分類(lèi)后帶來(lái)的誤差積累，但是受輻射差異影響較大，除進(jìn)行嚴(yán)格的輻射歸一化外，只能通過(guò)選擇同一傳感器和同一季相數(shù)據(jù)以減少噪聲(Deroseetal.， 2011； Maseketal.， 2008)。干擾監(jiān)測(cè)中，運(yùn)用幾期影像很難區(qū)分突發(fā)干擾和背景噪聲，而時(shí)間序列分析中，可利用背景噪聲和突發(fā)變化的時(shí)間序列軌跡不同將其分開(kāi)(Hostertetal.， 2003)。國(guó)外學(xué)者利用時(shí)間序列影像聯(lián)合分析，通過(guò)提取時(shí)序特征，進(jìn)而識(shí)別森林干擾，大大提高了工作效率和監(jiān)測(cè)精度(Huangetal.， 2009； Jamalietal.， 2015； Kennedyetal.， 2010； Sulla-Menasheetal.， 2014； Verbesseltetal.， 2010)。

福建是我國(guó)南方重點(diǎn)林區(qū)之一，森林覆蓋率約66%，居全國(guó)首位。森林在全球碳循環(huán)和碳管理中占有重要位置(Faheyetal.， 2010)。本研究以福建省長(zhǎng)汀縣為研究區(qū)，基于Landsat時(shí)序數(shù)據(jù)，采用時(shí)間序列軌跡分析(LandTrendr)方法進(jìn)行森林干擾監(jiān)測(cè)，以期為陸地森林生態(tài)系統(tǒng)碳循環(huán)和碳蓄積以及氣候變化研究提供參考。

1 材料與方法

1.1 研究區(qū)概況

長(zhǎng)汀縣地處福建省西部，25°18′40″—26°02′05″N，116°00′45″—116°39′20″E，南與廣東省近鄰，西與江西省接壤，呈四周凸起、中間下凹的盆地特征，地形破碎，嶺谷相間。該縣處中亞熱帶季風(fēng)氣候區(qū)，平均氣溫18.3 ℃，無(wú)霜期年均260天，四季分明。花崗巖發(fā)育的紅壤抗蝕能力弱，原生植被多遭破壞，現(xiàn)有植被主要包括馬尾松(Pinusmassoniana)、灌叢和荒草坡在內(nèi)的次生植被和人工植被。研究區(qū)地理位置見(jiàn)圖1。

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

選用長(zhǎng)汀縣2000—2016年15期Landsat影像數(shù)據(jù)(表1)，條帶號(hào)121/42，數(shù)據(jù)來(lái)源于美國(guó)地質(zhì)調(diào)查局(USGS)。所選影像盡量保持同一季相，以減少植被物候差異對(duì)光譜識(shí)別的影響。長(zhǎng)汀縣夏季雨水豐富，為獲取高質(zhì)量影像，將時(shí)間集中在9、10月，個(gè)別年份由于云覆蓋等因素不滿足要求，將時(shí)間間隔延長(zhǎng)至2年。數(shù)字高程資料來(lái)源于地理空間數(shù)據(jù)云平臺(tái)提供的全球30 m分辨率GDEM數(shù)字高程產(chǎn)品。

對(duì)每年影像進(jìn)行輻射定標(biāo)和集成于ENVI的基于MODTRAN4+輻射傳輸模型的FLAASH大氣校正，得到長(zhǎng)汀縣2000—2016年時(shí)序地表真實(shí)反射率影像。為消除云和云陰影對(duì)結(jié)果的干擾，利用Fmask算法(Zhuetal.， 2012)逐一去除云像元。對(duì)掩膜后留下的空值采用最鄰近年份的無(wú)云像元時(shí)間內(nèi)插方法進(jìn)行填充，對(duì)于時(shí)間序列影像的第1年和最后1年，使用最鄰近年份的無(wú)云像元替代：

xi=xp+(i-p)(xn-xp)/(n-p)。

(1)

式中：xi為第i年的填充值；xp和xn分別為第p年和第n年的指數(shù)值。

圖1 研究區(qū)地理位置

表1 遙感影像數(shù)據(jù)介紹

1.3 光譜指數(shù)選取

選用歸一化燃燒率(normalized burn ratio，NBR)構(gòu)建時(shí)序數(shù)據(jù)進(jìn)行干擾監(jiān)測(cè)。國(guó)外學(xué)者(Kennedyetal.， 2010； Cohenetal.， 2010)基于大量數(shù)據(jù)研究表明，NBR對(duì)多種類(lèi)型干擾事件具有最大的敏感性。NBR由對(duì)植被葉綠素含量較敏感的近紅外(NIR)和對(duì)水汽含量較敏感的短波紅外(SWIR)組合構(gòu)成，健康植被因高NIR值和低SWIR值，使其具有較高的NBR值，而受干擾的像元由于呈現(xiàn)土壤特征，SWIR值增大，使NBR值降低。據(jù)此，可根據(jù)時(shí)序軌跡中NBR值的突變識(shí)別森林變化信息：

(2)

1.4 時(shí)間序列軌跡擬合算法

LandTrendr算法(Kennedyetal.， 2010)的核心是將復(fù)雜的時(shí)序軌跡簡(jiǎn)化為一組相連的線段，并對(duì)分割后的分段進(jìn)行線性擬合，旨在消除噪聲，突出重要的信息。算法主要參數(shù)設(shè)置如表2所示，關(guān)鍵流程如下：

1) 時(shí)序軌跡提取及噪聲消除 逐像元提取時(shí)序軌跡(圖2a)，并通過(guò)設(shè)置噪聲值進(jìn)一步去除預(yù)處理過(guò)程中殘留的細(xì)小噪聲(圖2b)。

2) 時(shí)序軌跡分割 對(duì)整個(gè)時(shí)間序列，首先將時(shí)序中的第1年和最后1年作為第1次分割的起始點(diǎn)，得到觀測(cè)年份與指數(shù)值的一階回歸方程【圖2c(1)】，并計(jì)算每個(gè)觀測(cè)年指數(shù)值與預(yù)測(cè)值的絕對(duì)差，選取絕對(duì)差最大的一年作為下次分段的分段點(diǎn)，據(jù)此將時(shí)間序列分為2部分【圖2c(2)】。然后對(duì)2部分分別計(jì)算回歸方程和每段的均方誤差(MSE)，選取MSE較大的一段參與下次分割，以此類(lèi)推，直到分段數(shù)達(dá)到設(shè)定的初次分段數(shù)【圖2c(3)～(5)】。為防止分段數(shù)過(guò)多造成過(guò)擬合現(xiàn)象，選用角度閾值進(jìn)一步判斷，使分段數(shù)達(dá)到最大分段數(shù)。通過(guò)計(jì)算對(duì)比原始時(shí)序軌跡中段與段的夾角，去除最“淺”角對(duì)應(yīng)的分段【圖2c(6)】。

3) 時(shí)序軌跡擬合 聯(lián)合基于點(diǎn)與點(diǎn)連線及基于回歸連線2種方法進(jìn)行時(shí)序軌跡擬合。對(duì)于第1個(gè)分段，分別計(jì)算2種方法的MSE，選取MSE較小的方法作為分段的連接線。從下一分段開(kāi)始，每一個(gè)連接線必須與前一個(gè)相連，同樣計(jì)算2種方法的MSE，選取MSE較小的方法連接，最后得到一條完整的擬合軌跡。

4) 模型簡(jiǎn)化 采用恢復(fù)率閾值(1/n)剔除不合理的分段點(diǎn)，對(duì)上一過(guò)程擬合的時(shí)序軌跡通過(guò)迭代方式簡(jiǎn)化模型并重新計(jì)算擬合。在研究區(qū)內(nèi)，絕大多數(shù)森林干擾為突發(fā)性事件，但植被恢復(fù)是個(gè)緩慢的過(guò)程,通過(guò)設(shè)置恢復(fù)率閾值(1/n)，將植被恢復(fù)時(shí)間小于n年的森林干擾視為噪聲。根據(jù)長(zhǎng)汀縣植被恢復(fù)狀況，將恢復(fù)率閾值設(shè)為0.5，恢復(fù)時(shí)間小于2年，即干擾發(fā)生1年后，植被恢復(fù)到干擾發(fā)生前狀態(tài)的區(qū)域視為噪聲。

5) 光譜濾波 通過(guò)回歸模型將NBR值轉(zhuǎn)化為植被覆蓋率進(jìn)而生成干擾量，以植被變化強(qiáng)度閾值(表2濾波過(guò)程)為濾波條件，剔除原始分割算法中小的光譜分段點(diǎn)，以減少光譜異常信號(hào)或植被物候差異造成的時(shí)序軌跡過(guò)度擬合。

基于以上過(guò)程，可識(shí)別時(shí)間序列中NBR下降、上升和保持不變3種特性。經(jīng)過(guò)噪聲去除及對(duì)由光譜異常信號(hào)或植被物候差異造成的偽變化進(jìn)行濾波后，NBR值下降即由森林干擾導(dǎo)致。

圖2 LandTrendr分割流程

表2 LandTrend參數(shù)設(shè)置

1.5 森林干擾制圖

制圖前需先排除水體、裸地等非森林地物。基于NBR指數(shù)的時(shí)序分割算法獲取頂點(diǎn)的濕度、亮度和綠度影像，計(jì)算每年每個(gè)像元的纓帽變換角(TCA)。TCA值越大，對(duì)應(yīng)的森林覆蓋率越高，選擇平均值最高的TCA值，其對(duì)應(yīng)的像元最像森林像元(Huangetal.， 2002)。通過(guò)TCA閾值提取森林像元，提取的像元在整個(gè)時(shí)序中被識(shí)別為持續(xù)森林。由于研究所用時(shí)間序列較短，可監(jiān)測(cè)到的二次干擾區(qū)域很少，故只提取最大干擾作為監(jiān)測(cè)結(jié)果。TCA計(jì)算公式為：

式中： TCG和TCB分為纓帽變換綠度和濕度分量。

2 結(jié)果與分析

2.1 森林干擾監(jiān)測(cè)結(jié)果

圖3為長(zhǎng)汀縣2000—2016年森林干擾分布，圖中非森林和持續(xù)森林已被腌膜，以白色背景表示，其余不同顏色代表不同年份的干擾分布。

2.2 森林干擾時(shí)空分布特征

圖4為長(zhǎng)汀縣森林干擾面積隨時(shí)間的變化特征。經(jīng)統(tǒng)計(jì)分析， 2000—2016年長(zhǎng)汀縣森林干擾總面積192.49 km2，平均每年受擾動(dòng)森林12.83 km2； 2001年干擾量最小，不足1 km2； 2004、2008和2009年森林受干擾較為嚴(yán)重，均在30 km2以上，約占當(dāng)年森林面積的1.3%，共占干擾總面積的50%，其中2004年干擾面積高達(dá)32.85 km2； 2003、2006、2007和2010—2011年干擾面積略大于10 km2，其余年份干擾面積小于或遠(yuǎn)小于10 km2。

圖4 長(zhǎng)汀縣森林干擾面積隨時(shí)間的變化特征

從長(zhǎng)汀縣氣象局獲取的數(shù)據(jù)得知， 2004年是長(zhǎng)汀縣近20年來(lái)森林火災(zāi)次數(shù)最多的一年， 2008和2009年也爆發(fā)了大面積森林火災(zāi)(江帆等， 2015)。2008年受“拉尼娜”特大雨雪冰凍災(zāi)害侵襲，南方森林資源損失慘重(曹坤芳等， 2010)。2002—2004年、2006—2008年干擾面積呈持續(xù)上升趨勢(shì)， 2001和2005年干擾降低，特別是從2009年起干擾面積持續(xù)減少直到2016年近乎最低。2000—2016年森林干擾面積在個(gè)別年份波動(dòng)較大，但總體上隨時(shí)間呈下降趨勢(shì)。

對(duì)每年影像進(jìn)行森林提取，獲得長(zhǎng)汀縣歷年森林面積及所占土地總面積比例。由圖5可知， 2004、2008和2009年森林面積下跌較大， 2009年起呈緩慢增長(zhǎng)趨勢(shì)，與擾動(dòng)情況相呼應(yīng)。對(duì)比圖4和圖5，2004、2008和2009年的森林干擾面積占當(dāng)年森林面積的1.3%左右，其余各年均小于當(dāng)年森林面積的0.6%。

圖5 長(zhǎng)汀縣歷年森林面積及所占總面積比例

由圖6可知，長(zhǎng)汀縣森林干擾面積隨持續(xù)時(shí)間呈極速下降趨勢(shì)，后保持低值平緩，說(shuō)明干擾持續(xù)時(shí)間保持較短的年數(shù)，主要為1～2年，少數(shù)為3年，其中發(fā)生在1年的干擾面積比例最大，達(dá)82%； 發(fā)生在2和3年的干擾面積分別占13%和4%。結(jié)合Google影像目視解譯，長(zhǎng)汀縣森林干擾主要是由森林火災(zāi)和人工砍伐造成的急劇干擾事件。

圖6 長(zhǎng)汀縣森林干擾面積隨持續(xù)時(shí)間變化分布

由圖3各年份干擾分布可知，長(zhǎng)汀縣森林干擾面積具有較強(qiáng)的空間分布特征。2004、2008和2009年森林干擾斑塊面積較大，主要集中在長(zhǎng)汀縣中東部非森林區(qū)域附近，受人為干擾影響。2008年除分布在非森林區(qū)域附近的干擾外，有大量因受雨雪災(zāi)害影響的細(xì)小斑塊分布在長(zhǎng)汀縣西部。其他年份亦如此，主要分布在長(zhǎng)汀縣非森林區(qū)域附近，四周分布零散細(xì)小的斑塊。圖7清晰展現(xiàn)了長(zhǎng)汀縣森林干擾面積在各海拔段的分布情況，隨海拔升高，干擾呈明顯下降趨勢(shì)，其中海拔380～460 m之間的干擾面積最大，超過(guò)60%的干擾發(fā)生在中低海拔地區(qū)，海拔大于940 m的干擾僅占0.45%。

圖7 長(zhǎng)汀縣森林干擾面積隨海拔變化分布

2.3 基于地面調(diào)查的精度評(píng)估

圖8 長(zhǎng)汀縣森林干擾實(shí)地調(diào)查

結(jié)合森林干擾監(jiān)測(cè)結(jié)果和地面調(diào)查資料對(duì)4個(gè)典型樣區(qū)(樣區(qū)分布見(jiàn)圖3)進(jìn)行驗(yàn)證，獲取干擾發(fā)生的時(shí)間和位置。圖8a-d從左到右分別為樣區(qū)干擾年TM影像、監(jiān)測(cè)結(jié)果、實(shí)地調(diào)查照片和NBR指數(shù)時(shí)序軌跡，實(shí)地調(diào)查照片均拍攝于2013年7月。當(dāng)干擾發(fā)生時(shí)，NBR指數(shù)大幅度降低，之后逐漸增加，植被慢慢恢復(fù)。經(jīng)目視解譯并對(duì)照提取結(jié)果可知，干擾斑塊可被完整提出，且邊界準(zhǔn)確清晰，細(xì)小干擾也能逐一識(shí)別。圖8a、b干擾分別發(fā)生在2010和2008年，到2013年植被有所恢復(fù)但不是特別理想，從實(shí)地調(diào)查照片可知植被主要以灌木為主； 從圖8c、d可以看出，植被恢復(fù)效果良好。實(shí)地調(diào)查結(jié)果與研究得到的干擾監(jiān)測(cè)結(jié)果一致，表明該方法適用于長(zhǎng)汀縣森林干擾監(jiān)測(cè)，且具有較好的監(jiān)測(cè)能力。

2.4 基于像元尺度的精度驗(yàn)證

采用分層隨機(jī)抽樣方法，以像元為單位提取森林干擾樣本進(jìn)行精度驗(yàn)證。選用每年作為一層(包括15期森林干擾區(qū)域，森林和非森林各為1層)，共17層。以最少像素層30為基準(zhǔn)，對(duì)森林干擾層進(jìn)行等比例隨機(jī)抽樣，15期干擾樣本數(shù)量共9 075個(gè)，森林和非森林像元數(shù)偏多，適當(dāng)減少，分別為3 200和1 000個(gè)，即共13 275個(gè)像元樣本參與驗(yàn)證。驗(yàn)證結(jié)果統(tǒng)計(jì)為混淆矩陣，見(jiàn)表3。

由表3可知，時(shí)間序列軌跡算法可以很好提取長(zhǎng)汀縣2000—2016年森林干擾，總體精度達(dá)96.26%，Kappa系數(shù)為0.92，具有較高的監(jiān)測(cè)精度。各年份用戶精度均在80%以上， 2014和2016年相對(duì)其他年份出現(xiàn)一定漏提。2016年處于時(shí)序軌跡的末端且前1年影像缺失，加之目視解譯發(fā)現(xiàn)該年干擾范圍和斑塊面積較小，分布零星，導(dǎo)致2016年出現(xiàn)45%左右的漏提。

3 討論

LandTrendr算法根據(jù)數(shù)據(jù)本身決定變化軌跡，利用干擾與背景噪聲不同，通過(guò)參數(shù)控制進(jìn)行分割擬合，消除時(shí)序中的噪聲和小的光譜變化點(diǎn)，以突出重要的信息，最終得以識(shí)別干擾。該方法不僅可提取多種類(lèi)型的干擾，還可監(jiān)測(cè)同一區(qū)域多次擾動(dòng)變化，同時(shí)獲取植被干擾和恢復(fù)相關(guān)信息。

由于云和云陰影的影響，不能保證每年都有滿足要求的高質(zhì)量影像。某些年份影像缺失，可能導(dǎo)致干擾監(jiān)測(cè)時(shí)間延誤1年，或某些區(qū)域出現(xiàn)植被恢復(fù)過(guò)快現(xiàn)象導(dǎo)致輕微漏提。后續(xù)工作可重構(gòu)高時(shí)間分辨率影像，降低因影像缺失造成的漏提、誤提現(xiàn)象。經(jīng)目視解譯可知，長(zhǎng)汀縣森林干擾類(lèi)型主要為森林火燒和人工采伐，本研究所用方法雖可提取多種類(lèi)型的干擾，但未能劃分出具體的干擾類(lèi)型，后續(xù)工作可結(jié)合其他技術(shù)方法分別提取出火燒跡地和采伐跡地。

4 結(jié)論

針對(duì)我國(guó)森林干擾頻繁、干擾類(lèi)型復(fù)雜多樣的特點(diǎn)，采用LandTrendr時(shí)序軌跡分析方法研究福建省長(zhǎng)汀縣森林干擾狀況，并進(jìn)行干擾時(shí)空特征分析。結(jié)論如下：

1) 2000—2016年，長(zhǎng)汀縣森林干擾面積在個(gè)別年份波動(dòng)較大，但總體上隨時(shí)間呈下降趨勢(shì)。

2) 長(zhǎng)汀縣森林干擾主要是由森林火災(zāi)和人工砍伐造成的急劇干擾事件，且主要發(fā)生在非森林區(qū)域附近的低海拔地區(qū)。

3) 基于目視解譯和實(shí)地調(diào)查結(jié)果與研究得到的干擾監(jiān)測(cè)結(jié)果一致，干擾斑塊可被完整提出，且邊界準(zhǔn)確清晰，細(xì)小干擾也能逐一識(shí)別； 基于像元尺度精度驗(yàn)證的總體精度達(dá)96.26%，Kappa系數(shù)為0.92，各年份用戶精度均在80%以上，除個(gè)別年份外，生產(chǎn)者精度均在75%以上，具有較高的監(jiān)測(cè)精度。