基于STARFM 的草地地上生物量遙感估測(cè)研究
——以甘肅省夏河縣桑科草原為例

張玉琢，楊志貴，于紅妍，張強(qiáng)，楊淑霞，趙婷，許畫(huà)畫(huà)，孟寶平*，呂燕燕

（1. 南通大學(xué)脆弱生態(tài)研究所，地理科學(xué)學(xué)院，江蘇 南通 226007；2. 祁連山國(guó)家公園青海服務(wù)保障中心，青海 西寧 810001；3. 甘肅省環(huán)境監(jiān)測(cè)中心站，甘肅 蘭州 730020）

草地地上生物量（above ground biomass，AGB）通常是指草地地上部分植被干重［1］，是陸地生態(tài)系統(tǒng)碳循環(huán)的重要組成部分［2］，其時(shí)空格局分布可以反映草地生態(tài)系統(tǒng)功能和狀態(tài)的變化特征［3］。在實(shí)際生產(chǎn)活動(dòng)中，草地AGB 的變化可以用來(lái)監(jiān)測(cè)草地放牧程度，土地利用變化及早期草地退化預(yù)警［4］。因此，構(gòu)建精確的地上生物量估測(cè)模型對(duì)草地管理、草畜平衡、草地生長(zhǎng)狀況評(píng)估和生態(tài)環(huán)境保護(hù)具有非常重要的科學(xué)意義［5?7］。

青藏高原地區(qū)地形復(fù)雜，氣候環(huán)境獨(dú)特，傳統(tǒng)地面實(shí)測(cè)方法難以獲取高寒草地AGB 真實(shí)空間分布狀況。遙感因其具有全局、經(jīng)濟(jì)、速度快、負(fù)載大量信息等特點(diǎn)［8］，被廣泛應(yīng)用于大尺度草地生物量估測(cè)研究［9］。草地地上生物量估測(cè)中常用到的中低分辨率遙感資料主要有AVHRR、MODIS、ASTER 和SPOT 衛(wèi)星上搭載的VEGETATION，中高分辨率遙感資料主要包括Landsat 系列衛(wèi)星上搭載的MSS、TM、ETM+及OLI，環(huán)境減災(zāi)系列HJ-2A/2B CCD，高分系列衛(wèi)星（GF1/GF6 WFV），資源衛(wèi)星和Sentinel 衛(wèi)星等。Williamson 等［10］結(jié)合Landsat MSS 衛(wèi)星遙感資料和實(shí)測(cè)草地地上生物量數(shù)據(jù)，反演了澳大利亞半干旱地區(qū)草地地上生物量，為該地區(qū)火災(zāi)預(yù)警和草地植被生長(zhǎng)監(jiān)測(cè)提供了科學(xué)依據(jù)。Verbesselt 等［11］利用SPOT VEGETATION 遙感資料估測(cè)了南非克魯格國(guó)家公園薩瓦拉草原的植被生物量和含水量情況，從而為南非熱帶草原生態(tài)系統(tǒng)的火災(zāi)風(fēng)險(xiǎn)評(píng)估提供參考依據(jù)。徐斌等［12］基于MODIS NDVI 數(shù)據(jù)和全國(guó)草地地上生物量調(diào)查數(shù)據(jù)，分別構(gòu)建了我國(guó)6 大牧區(qū)草地地上生物量估測(cè)模型，模型估測(cè)精度總體上接近80%。孟寶平等［13］用MODIS、Landsat 8 OLI 和HJ-1A/1B CCD 等遙感資料的NDVI 對(duì)高寒草甸生物量估測(cè)精度進(jìn)行對(duì)比分析，結(jié)果表明，基于Landsat 8 OLI NDVI 對(duì)數(shù)模型估算的估測(cè)誤差最小，RMSE 僅為31.33 kg·hm?2。

長(zhǎng)期以來(lái)，使用Landsat MSS、TM、ETM+和OLI 等較高分辨率的衛(wèi)星資料（30 m）在草地AGB 監(jiān)測(cè)中有較高的精度，但受過(guò)境時(shí)間（16 d）、覆蓋范圍（185 km×185 km）和云雨等天氣因素的嚴(yán)重制約，無(wú)法獲取到持續(xù)長(zhǎng)時(shí)間序列、大范圍高質(zhì)量衛(wèi)星影像。而MODIS 等遙感資料空間分辨率較低（250～1000 m），在草地AGB 的估測(cè)研究中尚存在較大誤差［14］，但其擁有較高的時(shí)間分辨率（每天）和相對(duì)廣闊的空間覆蓋面積（幅寬達(dá)2330 km），非常適合應(yīng)用在分布遼闊的草地及其動(dòng)態(tài)變化監(jiān)測(cè)上。近年來(lái)，隨著遙感數(shù)據(jù)融合算法的發(fā)展，多源遙感數(shù)據(jù)融合產(chǎn)品在植被遙感監(jiān)測(cè)中得到了廣泛應(yīng)用，尤其是時(shí)空適應(yīng)反射率融合算法（spatial and temporal adaptive reflection fusion model，STARFM）［15］。如趙艷麗等［16］基于Landsat TM 和MODIS 反射率產(chǎn)品，利用STARFM 算法獲取了高時(shí)序特征的TM 數(shù)據(jù)，對(duì)內(nèi)蒙古呼倫貝爾呼倫湖濕地不同地物類(lèi)別的分類(lèi)研究表明，該方法得到的TM 時(shí)序數(shù)據(jù)集對(duì)研究區(qū)不同地物類(lèi)別具有較高的預(yù)測(cè)能力，尤其適用于在濕地這類(lèi)反射特點(diǎn)實(shí)時(shí)變化相對(duì)不大的地區(qū)。田羅［17］以MODIS NDVI 和Landsat（TM/OLI）NDVI 為數(shù)據(jù)源，結(jié)合STARFM 算法研究了2000?2015年隔湟水流域春小麥（Triticum aestivum）、青稞（Hordeum vulgare）和油菜（Brassica napus）3 種作物空間分布信息，并對(duì)流域內(nèi)農(nóng)作物種植構(gòu)造的時(shí)間空間演變進(jìn)行了探究。Peters［18］基于Landsat 和MODIS 遙感資料，利用STARFM 算法合成2000?2015年30 m 空間分辨率和16 d 時(shí)間分辨率的NDVI 植被指數(shù)數(shù)據(jù)集，從而分析了入侵物種野葛（Pueraria montana）物候動(dòng)態(tài)變化情況。以上研究均表明，STARFM 將Landsat 高空間分辨率和MODIS 高時(shí)間分辨率的優(yōu)勢(shì)結(jié)合起來(lái)，生成了具有長(zhǎng)時(shí)間序列的高空間分辨率遙感反射率數(shù)據(jù)，解決了時(shí)空數(shù)據(jù)融合的難題。然而，STARFM 算法在高寒草地地上生物量遙感估測(cè)中的應(yīng)用尚鮮有報(bào)道。

此外，STARFM 算法采用了鄰近相似像元的信息作為參考，對(duì)高分辨率影像像元值進(jìn)行預(yù)測(cè)，即以預(yù)測(cè)像元為中心，根據(jù)像元值相似程度、空間距離以及時(shí)間距離作為選擇依據(jù)，在滑動(dòng)窗口內(nèi)尋找預(yù)測(cè)像元的相似點(diǎn)，并對(duì)所有相似像元附以一定的權(quán)值，最終獲得中間預(yù)測(cè)像元的值，在整個(gè)圖像上進(jìn)行卷積運(yùn)算，從而獲得整幅影像的預(yù)測(cè)值。在對(duì)多源遙感資料進(jìn)行融合運(yùn)算時(shí)，還存在數(shù)據(jù)量大，運(yùn)行速度慢等特點(diǎn)［19］。自2008年美國(guó)地質(zhì)勘探局（United States Geological Survey，USGS）開(kāi)放了Landsat 影像后，谷歌開(kāi)發(fā)了基于云運(yùn)算的集成開(kāi)發(fā)環(huán)境Google Earth Engine（GEE）以便于調(diào)用地球引擎API（application programming interface）。GEE 囊括了Landsat和MODIS 等免費(fèi)多源遙感資料，集成了上千臺(tái)服務(wù)器，所有的進(jìn)程和計(jì)算都即時(shí)運(yùn)行［20］，在并行運(yùn)算和大數(shù)據(jù)算法下使得數(shù)據(jù)分析速度較單一服務(wù)器快40～1000 倍［19］。因其優(yōu)越的運(yùn)算性能，GEE 已在全球植被遙感觀測(cè)中得到廣泛的應(yīng)用，Patel 等［21］基于GEE，使用不同的規(guī)范化光譜向量在全球范圍分類(lèi)城市大小。Lobell 等［22］通過(guò)GEE 應(yīng)用Landsat 數(shù)據(jù)預(yù)估了美國(guó)玉米（Zea mays）和大豆（Glycine max）的產(chǎn)量。Hansen 等［20］通過(guò)GEE 繪制了全球森林覆蓋地圖。

基于上述因素考慮，以青藏高原東部的甘南地區(qū)典型研究區(qū)高寒草甸AGB 為研究對(duì)象，基于時(shí)間序列Landsat 系列遙感資料和MODIS 植被指數(shù)產(chǎn)品MOD13Q1 數(shù)據(jù)，結(jié)合GEE 平臺(tái)，采用STARFM 融合算法，探索研究區(qū)高時(shí)空分辨率植被指數(shù)的快速構(gòu)建方法，并構(gòu)建高寒草甸草地AGB 遙感反演模型，分析研究區(qū)草地AGB時(shí)空動(dòng)態(tài)變化趨勢(shì)，從而為高寒草甸草地AGB 快速、精確監(jiān)測(cè)研究提供科學(xué)依據(jù)，同時(shí)也為研究區(qū)草畜平衡和可持續(xù)發(fā)展提供理論基礎(chǔ)。

1 材料與方法

1.1 研究區(qū)概況

研究區(qū)位于青藏高原東緣的甘南藏族自治州夏河縣桑科鄉(xiāng)央吉社區(qū)（圖1a，地理坐標(biāo)為102°23′?102°26′E，35°05′?35°07′N(xiāo)），東西長(zhǎng)2.77 km，南北長(zhǎng)3.86 km，總面積1.61 km2。研究區(qū)年平均溫度是2.10 ℃，海拔3050 m，年降水量為580 mm，隸屬大陸性高原溫帶季風(fēng)氣候。研究區(qū)包括5 塊樣區(qū)，每塊樣區(qū)因利用方式和處理措施各異，植被群落和草地生物量有較大差異（圖1b）。其中，1 號(hào)樣區(qū)為劃破草皮+補(bǔ)播的人工改良草地，2 號(hào)樣區(qū)為放牧利用草地，3 號(hào)樣區(qū)為未施加培育草地（即對(duì)照樣地），4 號(hào)樣區(qū)為鼠丘補(bǔ)播草地，5 號(hào)樣區(qū)為施肥草地。研究區(qū)內(nèi)原生態(tài)草地類(lèi)別為高寒草甸，垂穗披堿草（Elymus nutans）為其區(qū)域優(yōu)勢(shì)植被物種。放牧家畜以甘加羊、牦牛為主。

1.2 野外采樣數(shù)據(jù)

野外采樣點(diǎn)布設(shè)如圖1b 和表1 所示，在研究區(qū)內(nèi)選擇13 個(gè)與MOD13Q1 植被指數(shù)產(chǎn)品柵格所對(duì)應(yīng)的區(qū)域作為外業(yè)觀測(cè)樣地，范圍大小為250 m×250 m。每個(gè)觀測(cè)樣地中選取30 m×30 m 范圍作為樣方采樣點(diǎn)布設(shè)區(qū)（圖2）并且在每個(gè)設(shè)置區(qū)域4 個(gè)角點(diǎn)和中心點(diǎn)設(shè)置0.5 m×0.5 m 大小的樣方5 個(gè)。2013?2016年間，每年7月初至8月底開(kāi)展野外觀測(cè)工作，該時(shí)間段內(nèi)研究區(qū)草地處于生長(zhǎng)旺季。4年來(lái)開(kāi)展野外觀測(cè)4 次，觀測(cè)樣地48 個(gè)，樣方240 個(gè)（表1）。每個(gè)樣方中AGB 采集采用齊地面刈割的方法，收集草地群落全部地上部分，采集的樣品置于透氣的布袋，帶回實(shí)驗(yàn)室在64 ℃恒定溫度的烘箱中烘干至恒重，利用電子天平稱(chēng)重并記錄每個(gè)樣方草地生物量干重，以同一樣地中所有樣方的草地生物量干重的平均值表征樣地草地生物量情況，用于后期草地生物量預(yù)估模型的建立。除此之外，采樣記錄的內(nèi)容包括采樣點(diǎn)的經(jīng)度、緯度、地面高程、海拔高度、覆蓋面積、草地地上部分生物量等指標(biāo)。此外，在研究?jī)?nèi)隨機(jī)選取20 個(gè)地面控制點(diǎn)，用于后期衛(wèi)星遙感資料的幾何精校正。

表1 研究區(qū)遙感影像及外業(yè)調(diào)查時(shí)間Table 1 Date remote sensing image and field survey

圖1 研究區(qū)概況Fig.1 Overview of the study areaa 為研究區(qū)位置；b 為草地AGB 采樣點(diǎn)分布；c、d 和e 分別為研究區(qū)2013年8月8日Landsat 8 OLI 真彩色合成圖像以及Landsat 8 OLI 和MOD13Q1 NDVI 植被指數(shù)。a is the location of the study area；b is the distribution of AGB sampling sites in grassland；c，d and e are the real-color composite images of Landsat 8 OLI，NDVI of Landsat 8 OLI and MOD13Q1 in the study area on August 8，2013.

圖2 研究區(qū)采樣點(diǎn)分布Fig.2 Distribution of sample sites in the study area

1.3 遙感數(shù)據(jù)獲取及其預(yù)處理

本研究所涉及的遙感數(shù)據(jù)包括MODIS MOD13Q1 NDVI 植被指數(shù)產(chǎn)品和Landsat 8 OLI 衛(wèi)星遙感數(shù)據(jù)資料，兩種遙感數(shù)據(jù)資料收集和預(yù)處理均在谷歌地球引擎平臺(tái)（Google Earth Engine，GEE）上云處理完成，利用JavaScript 語(yǔ)言整理搜集2000?2016年所有Landsat 影像及與其過(guò)境時(shí)間對(duì)應(yīng)的MODIS NDVI 植被指數(shù)產(chǎn)品數(shù)據(jù)。在實(shí)現(xiàn)數(shù)據(jù)融合之前，對(duì)MODIS、Landsat 和數(shù)據(jù)質(zhì)量進(jìn)行檢測(cè)，以去除云及陰影對(duì)運(yùn)算結(jié)果的影響。研究區(qū)內(nèi)被云、雪或陰影區(qū)域作為低質(zhì)量的柵格被裁剪。本研究中，根據(jù)低質(zhì)量柵格比例，將低質(zhì)量柵格占比高于5%的MODIS 和Landsat 數(shù)據(jù)對(duì)不參與融合運(yùn)算。研究區(qū)數(shù)據(jù)統(tǒng)計(jì)結(jié)果如圖3 所示，共收集394 對(duì)數(shù)據(jù)對(duì)，其中Landsat 低質(zhì)量柵格數(shù)小于5%的影像144 景，大于5%的影像124 景，缺失影像126 景。最終將兩類(lèi)遙感資料的投影轉(zhuǎn)化為Albers，以供后期數(shù)據(jù)融合使用。

圖3 MODIS 和Landsat 影像數(shù)據(jù)Fig.3 MODIS and Landsat image data

1.4 STARFM 算法

遙感數(shù)據(jù)的融合主要采用時(shí)空適應(yīng)反射率融合算法（spatial and temporal adaptive reflectance fusion model，STARFM）［15］，該算法首先通過(guò)同一時(shí)間（t1）的MODIS 和Landsat 影像，獲取影像空間分布的差異，再結(jié)合另外一個(gè)時(shí)間（t2）的低分辨率影像MODIS 對(duì)Landsat 進(jìn)行預(yù)測(cè)，在預(yù)測(cè)過(guò)程中采用滑動(dòng)窗口來(lái)減少低分辨率遙感數(shù)據(jù)像元邊界的問(wèn)題，利用給空間距離、光譜距離和時(shí)間距離賦權(quán)重值的辦法，獲取滑動(dòng)窗口中心像元值。其計(jì)算過(guò)程可用如下公式表達(dá)：

式 中：(xi，yi)為 像 元 點(diǎn) 坐 標(biāo)，w表示活動(dòng)窗口的 柵 格 數(shù)，L(xw2，yw2，t2)表示t2時(shí)刻Landsat NDVI 的預(yù)測(cè)值，M(xi，yj，t2)表示t2時(shí)刻MODIS NDVI 的值，L(xi，yj，t1)表示t1時(shí)刻Landsat NDVI 的值，M(xi，yj，t1)表示t1時(shí)刻MODIS NDVI 的值。Wijk表示權(quán)重，即滑動(dòng)窗口w內(nèi)像元對(duì)預(yù)測(cè)值的貢獻(xiàn)，由光譜距離（Sijk）、時(shí)間距離（Tijk）、空間距離（Dijk）3 項(xiàng)來(lái)確定，其計(jì)算公式為：

本研究中使用R 語(yǔ)言下利用STARFM v.1.2.1 算法［15］生成2000?2016年空間分辨率為30 m，時(shí)間分辨率為16 d 的高分辨率NDVI 植被指數(shù)產(chǎn)品，所有數(shù)據(jù)的融合均在GEE 云端完成，供下一步草地地上生物量遙感估測(cè)模型的構(gòu)建和動(dòng)態(tài)反演使用。

1.5 生物量模型構(gòu)建及其精度檢驗(yàn)

將研究區(qū)實(shí)地采樣的生物量數(shù)據(jù)作為因變量，與生物量采樣時(shí)間所對(duì)應(yīng)的Landsat NDVI（NDVILandsat）、MODIS NDVI（NDVIMODIS）和融合后NDVI（NDVISTARFM）產(chǎn)品分別作為自變量，使用SPSS 統(tǒng)計(jì)分析軟件分別將基于NDVILandsat、NDVIMODIS和NDVISTARFM的線性、對(duì)數(shù)、指數(shù)和乘冪4 類(lèi)回歸模型進(jìn)行構(gòu)建。本研究使用留一法交叉驗(yàn)證（leave-one-out cross validation，LOOCV）的方法來(lái)評(píng)估單因素參數(shù)模型的表現(xiàn)能力，該方法可以有效地泛化模型。將所有的自變量及其對(duì)應(yīng)的因變量以等樣本數(shù)n的方式分為n份，進(jìn)行交叉驗(yàn)證，每次選取1/n的樣本作為測(cè)試數(shù)據(jù)集，用來(lái)驗(yàn)證模型的估測(cè)能力，剩下的樣本作為訓(xùn)練集用來(lái)構(gòu)建草地生物量估測(cè)模型。利用每次測(cè)試集來(lái)驗(yàn)證訓(xùn)練集模型的誤差，每次構(gòu)建模型后利用測(cè)試集數(shù)據(jù)計(jì)算模型預(yù)測(cè)值與觀測(cè)值之間R2和均方根誤差（root mean square error，RMSE），重復(fù)選取測(cè)試集和訓(xùn)練集n次，直到所有的樣本都在測(cè)試集和訓(xùn)練集中出現(xiàn)為止，模型的估測(cè)能力用n次運(yùn)算得到的R2和RMSE 的平均值表示，R2值越大，RMSE 越小，模型的精度越高。其中RMSE 的計(jì)算公式為：

式中：yi為測(cè)試集草地地上生物量觀測(cè)值，y′i為測(cè)試集草地地上生物量預(yù)測(cè)值，地上生物量值的單位為kg·hm?2。

此外，由于研究區(qū)面積較小，且各樣地內(nèi)草地生長(zhǎng)狀況比較單一，變化較小。因此，本研究以2013 與2016年全部樣點(diǎn)的生物量的平均值作為真實(shí)值計(jì)算出5 個(gè)研究區(qū)的總生物量。將處于不同遙感資料下的AGB 反演模型得出的結(jié)果作為反演值，分別計(jì)算絕對(duì)誤差和相對(duì)誤差，在不相同的3 種模型下分別評(píng)價(jià)整個(gè)研究區(qū)及5 個(gè)樣區(qū)的總生物量估測(cè)精度。絕對(duì)誤差和相對(duì)誤差的計(jì)算公式如下：

式中：Δ 代表生物量的絕對(duì)誤差，x代表模型計(jì)算的總生物量，l表示來(lái)源地面實(shí)測(cè)點(diǎn)計(jì)算的生物量，δ 表示相對(duì)誤差（%），生物量值的單位為kg·hm?2。

1.6年際動(dòng)態(tài)變化分析

基于1.3 和1.4 中NDVILandsat、NDVIMODIS和NDVISTARFM的研究區(qū)最優(yōu)草地生物量估測(cè)模型，計(jì)算2000?2016年每年生長(zhǎng)季最大AGB，并利用Slope 變化趨勢(shì)分析法分析了17年來(lái)草地最大生物量年際動(dòng)態(tài)變化趨勢(shì)［23］。Slope 計(jì)算公式如下：

式中：i表示年份，取值范圍為1～17，Biomassi表示第i年生長(zhǎng)季最大草地生物量，n表示累計(jì)觀察的年份為17。

2 結(jié)果與分析

2.1 研究區(qū)實(shí)測(cè)生物量與對(duì)應(yīng)植被指數(shù)統(tǒng)計(jì)分析

2013?2016年研究區(qū)實(shí)測(cè)草地AGB 及其所對(duì)應(yīng)3 種植被指數(shù)統(tǒng)計(jì)性分析結(jié)果如表2 所示，研究區(qū)內(nèi)NDVILandsat、NDVIMODIS和NDVISTARFM3 類(lèi)植被指數(shù)變化情況相似，標(biāo)準(zhǔn)偏差（standard deviation，STD）介于0.08～0.11 之間，變異系數(shù)（coefficient of variation，CV）介于0.13～0.17 之間，3 類(lèi)植被指數(shù)中，NDVILandsat的平均值最大，為0.77，其次依次為NDVIMODIS和NDVISTARFM。而研究區(qū)內(nèi)草地地上生物量隨著氣候和生長(zhǎng)期的變化產(chǎn)生劇烈的變化，4年來(lái)STD 高達(dá)972.60 kg·hm?2，CV 為0.42。生物量最高達(dá)3997.33 kg·hm?2，最小僅為745.52 kg·hm?2，平均值為2299.30 kg·hm?2。

表2 研究區(qū)植被指數(shù)和草地生物量統(tǒng)計(jì)性描述Table 2 Statistical description of vegetation index and grassland biomass in the study area

2.2 生物量估算模型及精度檢驗(yàn)

基于LOOCV 檢驗(yàn)的NDVILandsat、NDVIMODIS和NDVISTARFM3 類(lèi)植被指數(shù)所構(gòu)建的線性、對(duì)數(shù)、乘冪和指數(shù)模型的精度情況如表3 所示。在4 類(lèi)模型中，由乘冪模型所建立的草地生物量預(yù)估模型精度最高，該模型R2介于0.24～0.76，RMSE 介于634.83～937.79 kg·hm?2，其次精度由高到低依次為指數(shù)模型和對(duì)數(shù)模型，而線性模型精度最低，R2介于0.17～0.65，RMSE 介于586.72～917.21 kg·hm?2。3 類(lèi)植被指數(shù)中，基于NDVILandsat的4 類(lèi)模型精度均優(yōu)于NDVIMODIS和NDVISTARFM所對(duì)應(yīng)的模型，R2介于0.65～0.76，RMSE 介于576.76～665.32 kg·hm?2，其次為NDVISTARFM和NDVIMODIS。所有模型中基于NDVILandsat的乘冪模型精度最高，R2達(dá)0.76，RMSE 為634.83 kg·hm?2；其次為基于NDVISTARFM的乘冪模型，R2為0.58，RMSE 為795.62 kg·hm?2；而基于NDVIMODIS最佳草地生物量預(yù)估模型在3 類(lèi)植被指數(shù)中最低，R2僅為0.24，RMSE 高達(dá)937.79 kg·hm?2（表3）。

表3 研究區(qū)草地生物量回歸模型精度檢驗(yàn)Table 3 Accuracy validation of biomass regression models in study area

對(duì)3 種植被指數(shù)的最佳草地生物量估測(cè)模型參數(shù)估計(jì)進(jìn)行T 檢驗(yàn)和回歸顯著性F檢驗(yàn)，結(jié)果如表4 所示。3 種植被指數(shù)的最佳預(yù)估模型全通過(guò)了顯著性水平P＜0.001 的F檢驗(yàn)和T 檢驗(yàn)，3 種植被指數(shù)中最優(yōu)預(yù)估模型公式如表5 所示。

表4 最優(yōu)反演模型參數(shù)T 檢驗(yàn)和回歸顯著性F 檢驗(yàn)Table 4 T test and F test for optimal inversion model

表5 基 于NDVILandsat、NDVIMODIS 和NDVISTARFM 最 優(yōu)AGB 估 測(cè)模型Table 5 The optimal AGB estimation model based on NDVILandsat，NDVIMODIS and NDVISTARFM

2.3 研究區(qū)生物量反演模型精度驗(yàn)證

如表6 所示，對(duì)于整個(gè)研究區(qū)而言，基于NDVILandsat最優(yōu)估測(cè)模型估算的草地AGB 絕對(duì)誤差和相對(duì)誤差均最小，分別為3.45萬(wàn)kg 和7.75%，其次為基于NDVISTARFM的最優(yōu)估測(cè)模型反演結(jié)果，其絕對(duì)誤差和相對(duì)誤差分別為5.75萬(wàn)kg 和15.95%。而基于NDVIMODIS最優(yōu)估測(cè)模型的反演結(jié)果精度最低，其絕對(duì)誤差達(dá)9.89萬(wàn)kg，相對(duì)誤差達(dá)22.22%。就5 個(gè)樣區(qū)而言，1～4 號(hào)樣區(qū)內(nèi)基于NDVISTARFM最優(yōu)估測(cè)模型的反演誤差最小，其次依次為基于NDVILandsat和NDVIMODIS的最優(yōu)估測(cè)模型；而5 號(hào)樣區(qū)中基于NDVILandsat最優(yōu)估測(cè)模型反演誤差最小，其次為基于NDVISTARFM和NDVIMODIS的最優(yōu)估測(cè)模型反演誤差。

表6 基于NDVILandsat、NDVIMODIS和NDVISTARFM最優(yōu)估測(cè)模型的研究區(qū)草地產(chǎn)草量的精度評(píng)價(jià)Table 6 Accuracy evaluation of the inversion biomass based on NDVILandsat，NDVIMODIS and NDVISTARFM

2.4 研究區(qū)草地AGB年際動(dòng)態(tài)變化

基于NDVISTARFM和NDVIMODIS最優(yōu)估測(cè)模型反演2000?2016年研究區(qū)年最大草地AGB，并分析其時(shí)空動(dòng)態(tài)變化狀況，結(jié)果如圖4 所示。基于NDVISTARFM最佳估測(cè)模型的Slope 動(dòng)態(tài)變化分析結(jié)果顯示，17年來(lái)研究區(qū)大部分草地AGB年增加量大于30 kg·hm?2，面積約占整個(gè)研究區(qū)面積的91.90%，主要分布在人類(lèi)活動(dòng)干擾較少的區(qū)域；年增加量在10～30 kg·hm?2之間的區(qū)域占整個(gè)研究區(qū)面積的3.82%，主要位于離居民地較遠(yuǎn)的道路、河道和溝谷地區(qū)；區(qū)內(nèi)年變化量保持基本不變和減少的區(qū)域分別占整個(gè)研究區(qū)面積的1.96%和2.30%，位于研究區(qū)內(nèi)居民聚集區(qū)及其周邊地區(qū)。基于NDVIMODIS反演的草地AGB 在過(guò)去17年的變化趨勢(shì)空間分異誤差較大，且與NDVISTARFM最佳估測(cè)模型分析結(jié)果差異較大。基于NDVIMODIS最優(yōu)估測(cè)模型計(jì)算的Slope 結(jié)果表明，研究區(qū)內(nèi)35.60%的區(qū)域草地AGB 的年變化量處于10～30 kg·hm?2之間，年增量大于30 kg·hm?2的區(qū)域占3.53%。年變化量在?10～10 kg·hm?2之間的區(qū)域占46.71%，而年減少量在10 kg·hm?2以上的區(qū)域僅占14.14%。

圖4 基于NDVISTARFM和NDVIMODIS最優(yōu)估測(cè)模型反演的草地AGB年最大值slope 動(dòng)態(tài)變化Fig. 4 The slope dynamic variation of annual maximum grassland AGB based on the optimal estimation model of NDVISTARFM and NDVIMODISa 和b 分別為基于NDVIMODIS和NDVISTARFM 最優(yōu)估測(cè)模型反演的草地AGB年最大值slope 動(dòng)態(tài)變化；AGB年變化量單位為kg·hm?2·year?1。a and b are the slope dynamic variation of annual maximum grassland AGB based on the optimal estimation models of NDVIMODIS and NDVISTARFM，respectively；The unit of annual variation of aboveground biomass is kg·hm?2·year?1.

3 討論

3.1 MODIS NDVI 在高寒草地AGB 遙感估測(cè)中的局限性與展望

MODIS 數(shù)據(jù)因其具有較高的時(shí)間分辨率（每天）和空間覆蓋范圍（2330 km），在青藏高原高寒草地AGB 遙感估測(cè)中得到廣泛的應(yīng)用［12，24?27］。然而受空間分辨率的影響（250～1000 m），在草地AGB 的估測(cè)研究中尚存在較大誤差［9，14］。因此，探索基于MODIS 數(shù)據(jù)高精度的草地生物量遙感估測(cè)新途徑，對(duì)于提高大區(qū)域草地遙感反演的精度具有重要意義［9］。本研究中，基于時(shí)空融合算法獲取了NDVISTARFM數(shù)據(jù)集，從而提高了MODIS NDVI 植被指數(shù)的空間分辨率，減少了土壤背景值對(duì)植被指數(shù)的影響，使研究區(qū)草地AGB 的總體預(yù)測(cè)精度提高了28.35%。除此之外，基于MODIS NDVI 的草地AGB 遙感估測(cè)模型的改進(jìn)還可以從以下幾個(gè)方面開(kāi)展：1）通過(guò)增加樣地內(nèi)觀測(cè)樣方的數(shù)量和面積，來(lái)提高觀測(cè)樣地空間代表性，從而減少地面采樣與衛(wèi)星影像像元之間空間尺度的不匹配性［23，28］；2）根據(jù)衛(wèi)星影像成像時(shí)間，合理安排野外調(diào)查時(shí)間，盡可能減小地面調(diào)查與衛(wèi)星影像獲取之間的時(shí)間差異；3）結(jié)合新的遙感觀測(cè)技術(shù)（如無(wú)人機(jī)高光譜遙感監(jiān)測(cè)技術(shù)），加強(qiáng)草地植被群落光譜特征研究和窄波段遙感植被指數(shù)在草地AGB 遙感監(jiān)測(cè)中的應(yīng)用［29?30］；4）基于遙感植被指數(shù)、土壤、地形、氣候因子和生物物理指標(biāo)等因子，結(jié)合多源統(tǒng)計(jì)模型和機(jī)器學(xué)習(xí)算法構(gòu)建草地AGB 估算模型［14，31?32］。

3.2 STARFM 融合在本研究中的不足

常用的MODIS NDVI 數(shù)據(jù)時(shí)間分辨率較高，但空間分辨率低（分辨率低于250 m），生物量遙感估測(cè)模型精度低；Landsat NDVI 數(shù)據(jù)具有30 m 空間分辨率，生物量遙感估測(cè)模型精度高，但是重訪周期長(zhǎng)（16 d），且易受云雨等天氣因素影響，無(wú)法獲取連續(xù)時(shí)間序列的高質(zhì)量植被指數(shù)數(shù)據(jù)集。因此結(jié)合多源遙感資料和數(shù)據(jù)融合算法是提高草地生物量遙感估測(cè)的重要途徑。Gao 等［15］結(jié)合高空間分辨率（185 m×185 m）的Landsat 和高時(shí)間分辨率（每天）的MODIS 建立了STARFM 時(shí)空適應(yīng)性反射率融合模型，此方法將兩種來(lái)源的遙感數(shù)據(jù)進(jìn)行地球表面反射率的融合估測(cè)，已經(jīng)在監(jiān)測(cè)地物季節(jié)性變化［33］，提高蒸散、總初級(jí)生產(chǎn)力（gross primary productivity，GPP）分析精度［34］及提高分類(lèi)精度方面得到成功應(yīng)用［35］。從圖5 中可以看出MODIS 數(shù)據(jù)的空間分辨率較低（圖5a），而Landsat 數(shù)據(jù)又時(shí)常受到云和陰影影響（圖5b），在STARFM 融合下，提高了植被指數(shù)空間分辨率，并且一定程度上改善了云雨及陰影干擾的問(wèn)題（圖5c）。然而，就本研究而言，STARFM 算法在研究區(qū)草地AGB 估測(cè)中還存在一些不足和不確定性，首先本研究只探討了Landsat 和MODIS NDVI 融合數(shù)據(jù)產(chǎn)品，并未探討其他植被指數(shù)和波段對(duì)研究區(qū)草地AGB 的估測(cè)情況，Meng 等［9］的研究結(jié)果表明，研究區(qū)內(nèi)MODIS EVI 對(duì)草地AGB 的估測(cè)精度更高。此外，孟寶平等［13］研究表明基于Landsat 8 OLI band7/band5 比值植被指數(shù)對(duì)高寒草甸植被比較敏感，然而本研究并沒(méi)有討論此類(lèi)植被指數(shù)與MODIS 數(shù)據(jù)的融合產(chǎn)品對(duì)AGB 的估測(cè)情況，因此探索研究區(qū)內(nèi)其他植被指數(shù)或敏感波段融合數(shù)據(jù)對(duì)草地生物量預(yù)測(cè)精度的作用是下一步研究的重點(diǎn)方向。其次，利用STARFM 對(duì)長(zhǎng)時(shí)間序列數(shù)據(jù)進(jìn)行融合運(yùn)算時(shí)，數(shù)據(jù)量較大，運(yùn)算速度較慢，對(duì)數(shù)據(jù)運(yùn)算平臺(tái)要求較高，本研究基于前期研究經(jīng)驗(yàn)，選擇在GEE 平臺(tái)對(duì)數(shù)據(jù)進(jìn)行云運(yùn)算，大幅提高了數(shù)據(jù)處理進(jìn)程。

圖5 3 種NDVI 對(duì) 比Fig.5 Three types of NDVI comparison（a）為MODIS NDVI，（b）為L(zhǎng)andsat NDVI，（c）為融合后NDVI。（a）is MODIS NDVI，（b）is Landsat NDVI，（c）is NDVI after fusion.

3.3 樣區(qū)的絕對(duì)誤差和相對(duì)誤差的數(shù)據(jù)選擇

利用研究區(qū)各樣區(qū)內(nèi)草地生物量總差來(lái)檢驗(yàn)?zāi)Ｐ凸罍y(cè)精度時(shí)，應(yīng)選取草地生長(zhǎng)狀況均一，具有代表性的樣區(qū)（利用單產(chǎn)乘以樣區(qū)面積作為整個(gè)樣區(qū)草地AGB 的真實(shí)值）。在本研究中，各樣區(qū)2013?2016年草地地上生物量統(tǒng)計(jì)性分析結(jié)果如表7 所示。結(jié)果表明，2015年各樣區(qū)內(nèi)樣方的草地生物量STD 和CV 在4年內(nèi)最低，該年產(chǎn)草量總體偏低，可能與氣候影響有關(guān)，僅為2263.90 kg·hm?2，屬于異常年份。故該年產(chǎn)草量不能作為驗(yàn)證數(shù)據(jù)。此外，同時(shí)對(duì)比2013，2014 和2016年各樣地草地生物量狀況，2014年5 個(gè)樣區(qū)內(nèi)樣方生物量的STD 和CV 值多大于2013 和2016年，故該年研究區(qū)內(nèi)草地生長(zhǎng)均勻程度不如2013 和2016年。因此，本研究?jī)H選取2013 和2016年的實(shí)測(cè)數(shù)據(jù)去驗(yàn)證模型反演產(chǎn)草量精度。

表7 各年份的實(shí)測(cè)數(shù)據(jù)統(tǒng)計(jì)性分析Table 7 Statistical analysis of measured data in each year

4 結(jié)論

本研究通過(guò)結(jié)合Landsat 和MODIS NDVI 數(shù)據(jù)，利用STARFM 時(shí)空數(shù)據(jù)融合算法和GEE 云計(jì)算平臺(tái)，生成了2000?2016年甘南州夏河縣研究區(qū)時(shí)間分辨率為16 d，空間分辨率為30 m 的NDVI 數(shù)據(jù)集合。在此基礎(chǔ)上建立了研究區(qū)高寒草甸AGB 估測(cè)模型，探討了高寒草甸AGB 時(shí)空動(dòng)態(tài)變化狀況。主要得到以下結(jié)論：1）在基于NDVISTARFM的4 類(lèi)遙感估測(cè)模型中，乘冪模型的精度最高，其R2和RMSE 分別為0.58 和795.62 kg·hm?2，相較于MODIS NDVI 所構(gòu)建的最佳草地AGB 估測(cè)模型，R2提高了0.34，RMSE 減少了142.17 kg·hm?2；2）NDVISTARFM最優(yōu)估測(cè)模型對(duì)草地產(chǎn)草量的估測(cè)精度次于NDVILandsat最佳估測(cè)模型，而優(yōu)于NDVIMODIS最佳估測(cè)模型，其絕對(duì)誤差和預(yù)測(cè)精度分別為5.75萬(wàn)kg 和84.05%；3）2000?2016年研究區(qū)草地生物量超出90%的區(qū)域展現(xiàn)出增加的變化趨勢(shì)，相較于MODIS NDVI 最優(yōu)AGB 估測(cè)模型，基于NDVISTARFM最優(yōu)估測(cè)模型可以精確、詳細(xì)的反演研究區(qū)草地生物量時(shí)空變化特征。