1988-2017年黃山自然與文化遺產(chǎn)地森林?jǐn)_動數(shù)據(jù)集

張少宇，朱嵐巍，霍艷輝，楊旭，唐世林

1. 中國科學(xué)院遙感與數(shù)字地球研究所，數(shù)字地球重點(diǎn)實驗室，北京 100094

2. 東華理工大學(xué)，南昌 330013

3. 華北理工大學(xué)，河北唐山 063210

4. 清華同方威視技術(shù)股份有限公司，北京 100094

5. 中國科學(xué)院南海海洋研究所廣東省海洋遙感重點(diǎn)實驗室，廣州 510220

數(shù)據(jù)庫（集）基本信息簡介

數(shù)據(jù)庫（集）組成本數(shù)據(jù)集主要包括1 9 8 8-2 0 1 7年長時間序列3 0年的黃山自然遺產(chǎn)地的森林?jǐn)_動結(jié)果數(shù)據(jù)。這些數(shù)據(jù)保存為一個壓縮文件（1 9 8 8-2 0 1 7年黃山自然與文化遺產(chǎn)地森林?jǐn)_動數(shù)據(jù)集.z i p），其中包含4個文件夾（地理要素、源數(shù)據(jù)、符號系統(tǒng)、制圖）。

引 言

世界遺產(chǎn)主要包括自然遺產(chǎn)、文化遺產(chǎn)、混合遺產(chǎn)3種，是由聯(lián)合國教科文組織和世界遺產(chǎn)委員會發(fā)起并確認(rèn)的世界人民公認(rèn)的、目前無法替代的財富，是全人類最具有代表性和最具有普遍價值的文物古跡及自然景觀。其中自然遺產(chǎn)包含奇特的自然景觀、完整的生態(tài)系統(tǒng)、原始的生物保護(hù)圈，突出體現(xiàn)了大自然的鬼斧神工、動植物的良好生活家園、生物進(jìn)化的顯著特征等。因此世界遺產(chǎn)逐漸受到人類的廣泛重視，越來越多的國家和人民加入保護(hù)自然遺產(chǎn)、保護(hù)世界遺產(chǎn)的行列當(dāng)中。

黃山在1990年被評為世界自然與文化遺產(chǎn)，在2015年被評為世界地質(zhì)公園，在2018年被列入世界生物圈保護(hù)區(qū)，因此黃山在中國乃至世界上都占據(jù)非常高的地位。本文利用遙感空間技術(shù)，在自然遺產(chǎn)的研究課題下針對中國的自然與文化遺產(chǎn)地-黃山進(jìn)行長時間序列的監(jiān)測研究，并利用云計算平臺及森林?jǐn)_動提取算法得到1988-2017年30年的森林?jǐn)_動數(shù)據(jù)產(chǎn)品，并進(jìn)行數(shù)據(jù)共享。相關(guān)的研究成果作為自然遺產(chǎn)的研究成果的一部分，反映了自然遺產(chǎn)的歷史及現(xiàn)在的生存狀況，同時也可為相關(guān)學(xué)者研究其他自然遺產(chǎn)地提供相關(guān)參考。

1 數(shù)據(jù)采集和處理方法

1.1 數(shù)據(jù)來源

本研究使用的云計算在線處理平臺，具有高效、快速、方便等多個特點(diǎn)，解決了數(shù)據(jù)年份過多、研究區(qū)過大、數(shù)據(jù)下載繁瑣等眾多問題。在線獲取的原始數(shù)據(jù)集是美國地質(zhì)勘探局（http://glovis.usgs.gov/）提供的Landsat產(chǎn)品數(shù)據(jù)集。Landsat具有較高的分辨率（30 m）、較長的時間跨度、開放式使用等眾多特點(diǎn)[1]，對長時間序列監(jiān)測研究具有顯著優(yōu)勢。區(qū)域尺度的植被變化監(jiān)測研究，主要使用的是Landsat影像數(shù)據(jù)[2-4]。通過人工預(yù)覽選取30年中較好的數(shù)據(jù)產(chǎn)品（表1），并進(jìn)行相應(yīng)的預(yù)處理得到可在線使用的1988-2017年的基礎(chǔ)數(shù)據(jù)集。

表1 黃山自然遺產(chǎn)地遙感影像數(shù)據(jù)列表

序號 成像日期 衛(wèi)星 傳感器 軌道號9 1 9 9 6-0 4-2 2 L a n d s a t 5 T M 1 2 0-0 3 9 1 0 1 9 9 7-1 1-0 3 L a n d s a t 5 T M 1 2 0-0 3 9 1 1 1 9 9 8-0 4-2 8 L a n d s a t 5 T M 1 2 0-0 3 9 1 2 1 9 9 9-1 2-2 7 L a n d s a t 5 T M 1 2 0-0 3 9 1 3 2 0 0 0-1 0-1 0 L a n d s a t 5 T M 1 2 0-0 3 9 1 4 2 0 0 1-0 2-1 5 L a n d s a t 5 T M 1 2 0-0 3 9 1 5 2 0 0 2-0 9-3 0 L a n d s a t 5 T M 1 2 0-0 3 9 1 6 2 0 0 3-0 2-0 5 L a n d s a t 5 T M 1 2 0-0 3 9 1 7 2 0 0 4-0 4-2 0 L a n d s a t 7 E T M 1 2 0-0 3 9 1 8 2 0 0 5-0 3-0 6 L a n d s a t 7 E T M 1 2 0-0 3 9 1 9 2 0 0 6-0 5-2 0 L a n d s a t 5 T M 1 2 0-0 3 9 2 0 2 0 0 7-0 1-0 7 L a n d s a t 7 E T M 1 2 0-0 3 9 2 1 2 0 0 8-0 2-2 7 L a n d s a t 7 E T M 1 2 0-0 3 9 2 2 2 0 0 9-0 3-0 9 L a n d s a t 5 T M 1 2 0-0 3 9 2 3 2 0 1 0-0 3-2 8 L a n d s a t 5 T M 1 2 0-0 3 9 2 4 2 0 1 1-0 5-1 8 L a n d s a t 5 T M 1 2 0-0 3 9 2 5 2 0 1 2-0 4-2 6 L a n d s a t 7 E T M 1 2 0-0 3 9 2 6 2 0 1 3-0 4-0 7 L a n d s a t 8 O L I 1 2 0-0 3 9 2 7 2 0 1 4-1 0-2 5 L a n d s a t 7 E T M 1 2 0-0 3 9 2 8 2 0 1 5-0 2-0 6 L a n d s a t 8 O L I 1 2 0-0 3 9 2 9 2 0 1 6-0 3-2 8 L a n d s a t 8 O L I 1 2 0-0 3 9 3 0 2 0 1 7-1 0-0 9 L a n d s a t 8 O L I 1 2 0-0 3 9

1.2 數(shù)據(jù)處理方法

近幾十年來，許多遙感變化檢測方法被提出、應(yīng)用、對比分析和評價[2]。同時也出現(xiàn)了很多監(jiān)測森林?jǐn)_動的指數(shù)，如MODIS全球擾動指數(shù)（MGDI）、基于纓帽變化的擾動指數(shù)（DI）、歸一化濕度指數(shù)（NDMI）、綜合森林特征指數(shù)（IFZ）、歸一化植被指數(shù)（NDVI）、歸一化燃燒比指數(shù)（NBR）等[3]。但是每個指數(shù)具有不同的特點(diǎn)，針對不同的擾動類型具有不同的監(jiān)測能力，同時眾多指數(shù)對于不同的研究區(qū)也會出現(xiàn)一定的差異性。本文參考眾多文獻(xiàn)進(jìn)行相關(guān)總結(jié)分析，選取DI指數(shù)作為研究方法的基礎(chǔ)，并采取因地制宜的方式進(jìn)行改進(jìn)，對黃山自然遺產(chǎn)地進(jìn)行監(jiān)測。

根據(jù)森林?jǐn)_動監(jiān)測算法原理，基于森林?jǐn)_動DI算法加以改進(jìn)，DI指數(shù)在根據(jù)對影像波段的進(jìn)行纓帽變化后，提取其主要的分量進(jìn)行數(shù)學(xué)運(yùn)算，以獲取主要的森林?jǐn)_動信息，其改進(jìn)后DI算法的主要原理如下[5]：

式（1）中Bt、Gt、Wt分別為標(biāo)準(zhǔn)化后的亮度、綠度、和濕度。式（2）中B為纓帽變化后第一個分量波段代表亮度，Bμ表示第一分量波段平均值，Bσ為第一分量波段方差。式（3）、式（4）同理分別對應(yīng)纓帽變化后第二分量和第三分量，代表綠度和濕度以及其平均值和方差。式（5）中DIμ代表森林?jǐn)_動指數(shù)平均值，根據(jù)其亮度平均值、綠度平均值以及濕度平均值計算而來；式（6）中DIm表示由多年份森林?jǐn)_動指數(shù)平均值計算而來，n為時間年份。

擾動像元以及恢復(fù)像元判別：

θ為閾值，閾值調(diào)整首先參考DIm平均值的大小，一般往大于DIm的方向選擇。連續(xù)三年DI值大于θ，即定義為擾動像元；即DIn、DIn+1、DIn+2均大于θ，則DIn為第n年擾動值，此時DIn對應(yīng)的該像元為擾動像元。DIn、DIn+1、DIn+2中DIn為擾動像元，DIn+1、DIn+2均小于θ，則DIn對應(yīng)的該像元為恢復(fù)像元。本次實驗通過計算DIm值，以及參考持久森林樣本DI值，綜合考慮后本實驗初次θ閾值選取為4.0，以步長為0.5，分別進(jìn)行多次實驗；對檢測結(jié)果的范圍大小、變化程度進(jìn)行比較，經(jīng)過多次試驗后確定閾值。

本文的研究路線如圖1，利用Google Earth Engine云計算在線處理平臺，先進(jìn)行原始數(shù)據(jù)的大氣校正、GEE中Landsat.simple Cloud Score函數(shù)及image.updateMask函數(shù)進(jìn)行去云掩膜；然后進(jìn)行監(jiān)督分類、計算NDVI濾波；然后結(jié)合Landsat TOA數(shù)據(jù)產(chǎn)品進(jìn)行纓帽變換、擾動提取、結(jié)果導(dǎo)出等。以上工作大部分可依靠編程完成。在數(shù)據(jù)平臺導(dǎo)出數(shù)據(jù)后，然后通過ArcGIS軟件進(jìn)行處理，同時利用 2017年哨兵數(shù)據(jù)、Google Earth高分辨率影像、實地驗證進(jìn)行結(jié)果分析和結(jié)果驗證，最終生成1988-2017年共30年逐年森林?jǐn)_動產(chǎn)品數(shù)據(jù)集。

圖1 森林?jǐn)_動研究技術(shù)路線

2 數(shù)據(jù)樣本描述

2.1 數(shù)據(jù)組成

本數(shù)據(jù)集主要包括1988-2017年長時間序列30年的黃山自然遺產(chǎn)地的森林?jǐn)_動結(jié)果數(shù)據(jù)，這些數(shù)據(jù)保存為一個壓縮文件。壓縮文件中包含4個文件夾（地理要素、源數(shù)據(jù)、符號系統(tǒng)、制圖），總數(shù)據(jù)量為8.19 MB。地理要素文件夾中包含建筑、水體、邊界、索道、公路、步道對應(yīng)的矢量數(shù)據(jù)，文件為shp格式，可直接進(jìn)行查看編輯；此邊界為相關(guān)課題邊界研究人員參考多方資料確定的黃山自然遺產(chǎn)地的邊界范圍；其他要素為本研究工作人員搜集數(shù)據(jù)，結(jié)合地形圖，規(guī)劃圖矢量化而成。源數(shù)據(jù)文件夾中每個文件夾單獨(dú)存儲每一年的對應(yīng)的tif柵格數(shù)據(jù)文件，共有6個屬性值，分別為未擾動森林、擾動森林、恢復(fù)森林、建筑裸巖、耕地和水體，采用的UTM-WGS84坐標(biāo)系；符號系統(tǒng)文件夾中為ArcGIS軟件中屬性值中的符號系統(tǒng)，此為模板方便賦值使用，格式為lyr可直接使用或編輯；mxd格式文件為ArcGIS軟件格式文件，包含所有的源文件信息可直觀進(jìn)行查看；制圖文件夾中為1988-2017年30年的制圖成果，格式為jpg。

2.2 數(shù)據(jù)樣本

圖2為黃山自然遺產(chǎn)地1988-2017年中6期的森林?jǐn)_動監(jiān)測產(chǎn)品。通過GEE平臺中的監(jiān)督分類的方式，選取足夠多的樣本進(jìn)行分類，主要分為森林、水體、建筑、裸巖、耕地；然后通過NDVI濾波再次篩選森林像元。森林?jǐn)_動信息的提取是對森林像元做變化監(jiān)測，根據(jù)連續(xù)3年的處理像元值的比較來判斷是否為擾動像元或恢復(fù)像元。使用的方法是GEE平臺中Join函數(shù)方法，查看當(dāng)前年份前后年份的像元值是否大于θ，只要存在連續(xù)3年大于θ，均定義為擾動像元。例如，起始年份1988年只查看當(dāng)前年份及后一年數(shù)據(jù)即1989年像元值是否大于θ，以及終止年份只查看當(dāng)前年份及左年份即2016年像元值是否大于θ，以此定義是否為擾動像元。森林?jǐn)_動點(diǎn)位主要集中在旅游路線的主干道兩側(cè)、山頂?shù)馁e館附近、水系附近、較高海拔、以及森林與裸巖的交界處。主要的擾動類型為森林-裸巖、森林-建筑、森林-土地以及植被的自然變遷。例如，2006年黃山自然遺產(chǎn)地的保護(hù)面積為160.6平方公里，緩沖區(qū)面積為490平方公里，2006年的森林?jǐn)_動面積接近7.6萬平方米，恢復(fù)面積為6.8萬平方米，擾動面積占據(jù)總覆蓋面積的0.0473%，恢復(fù)面積占據(jù)總覆蓋面積的0.0423%。

圖2 黃山自然與文化遺產(chǎn)地6期森林?jǐn)_動圖

3 數(shù)據(jù)質(zhì)量控制和評估

3.1 數(shù)據(jù)結(jié)果分析

根據(jù)結(jié)果統(tǒng)計可得1988-2017年30年中每年擾動面積和恢復(fù)面積見圖3。根據(jù)30年的森林?jǐn)_動趨勢可看出，早期及2011年前后具有大面積的擾動，2008年前后出現(xiàn)大面積擾動，隨后未發(fā)生較好的恢復(fù)。最近10年雖一直頻繁出現(xiàn)不同程度的擾動，但整體占據(jù)比例很小，仍在自然變化的可接受范圍之內(nèi)。在每次擾動之后均有一定幅度的面積恢復(fù)，分析可知主要以森林自然變化為主。選取每年的具體擾動點(diǎn)位進(jìn)行比較來看，大多擾動點(diǎn)位出現(xiàn)在低等植物與高等植物、裸巖交界之處；由于數(shù)據(jù)質(zhì)量的限制，不得不采用每年份相對較好的影像，因此就滿足不了同一月份進(jìn)行比較；同時由于黃山季節(jié)變化也會出現(xiàn)很大的差異，數(shù)據(jù)選取問題給提取精度帶來影響，造成一定程度的擾動信息誤提取，導(dǎo)致提取精度一般；但鑒于此方法高效、快捷、可同時進(jìn)行長時間序列分析，因此具有一定的普遍適用性和可擴(kuò)展性，可優(yōu)選其使用數(shù)據(jù)并用于其他地區(qū)研究。

圖3 森林恢復(fù)和擾動面積列表

3.2 精度驗證方式

采用輔助數(shù)據(jù)，2017年哨兵數(shù)據(jù)、Landsat原始數(shù)據(jù)、高分影像Google Earth數(shù)據(jù)進(jìn)行目視解譯。采用目視解譯的方法，根據(jù)隨機(jī)分層抽樣原理[6]；由于Google高分影像限制，選取2008-2017年共10年每年的擾動點(diǎn)，每層隨機(jī)抽樣30個樣本；未擾動森林、非森林每層隨機(jī)抽取樣本100個，共500個像元點(diǎn)。根據(jù)精度評價混淆矩陣，計算制圖精度、用戶精度、總體精度、誤檢率、漏檢率、Kappa系數(shù)加以評價?？傮w精度為80.66%，Kappa系數(shù)為0.77，見表2。

表2 精度混淆矩陣及評價表

預(yù)測值真實值持續(xù)森林 非森林 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 總計 用戶精度 誤檢率%%漏檢率% 30.37 11.82 3.80 0.00 3.70 0.00 0.00 2.90 1.71 9.70 9.40 0.00 80.66 0.77

通過實地驗證的方式，共選取 4種典型樣本點(diǎn)（圖 4），驗證點(diǎn)位分別為 a（118°09'42.27″E、30°04'00.78″N），b（118°10'38.27"E、30°4'30.89"N），c（118°14'15.68"E、30°8'32.67"N），d（118°9'32.237"E、30°8'36.161"N）。圖4a為早期擾動，現(xiàn)已經(jīng)完全恢復(fù)，森林密集程度走勢較好，左側(cè)為該樣本點(diǎn)DI值的走勢圖，由圖可知當(dāng)?shù)厣指采w良好；圖4b為早期多次擾動，近幾年逐漸恢復(fù)，由實地情況可知為較新樹種覆蓋良好。圖4c為持續(xù)未擾動森林，一直處于較好的狀態(tài)，由圖可知此樣本相比圖4a、4b兩樣本明顯有更好的覆蓋效果。圖 4d為近幾年發(fā)生擾動但沒有完全恢復(fù)樣本，經(jīng)過現(xiàn)場查看可知此處多為植被生長變化造成的自然擾動，由于是裸巖與植被交界處，受到季節(jié)植被變化影響，容易出現(xiàn)監(jiān)測擾動值持續(xù)較高，由此可見此處的覆蓋度較差。

圖4 4種典型樣本點(diǎn)

4 數(shù)據(jù)使用方法和建議

數(shù)據(jù)使用的tif格式、shp格式的文件、lyr格式文件、mxd格式文件，便于研究工作者查看和進(jìn)行編輯。本文采用ArcGIS軟件進(jìn)行編輯繪制，同時數(shù)據(jù)集中也包含了本文所設(shè)置的圖層模板，便于更好的查看和制圖。本數(shù)據(jù)集具有很強(qiáng)的代表性，可以為眾多研究學(xué)者研究森林變化、植被變化、自然遺產(chǎn)地監(jiān)測提供一定程度的參考。同時本研究數(shù)據(jù)集可直接用于黃山遺產(chǎn)地保護(hù)、規(guī)劃、可持續(xù)發(fā)展等研究。

致 謝

非常有幸能夠參與到自然遺產(chǎn)的保護(hù)研究工作當(dāng)中，能夠有機(jī)會致力于遺產(chǎn)的保護(hù)和發(fā)展；同時感謝楊瑞霞副研究員在驗證工作中給予指導(dǎo)，感謝杜文杰、王普、陳亞亞等同學(xué)在研究過程中提供了很多幫助。