張 偉，李 瑋，陶冠宏，李愛農(nóng)，覃志豪，雷光斌，陳藝曦

·農(nóng)業(yè)信息與電氣技術(shù)·

利用STARFM模型提高復(fù)雜地表下復(fù)種指數(shù)遙感提取精度

張 偉1,2，李 瑋1,2，陶冠宏1,2，李愛農(nóng)3，覃志豪4，雷光斌3，陳藝曦1,2

（1. 中電科大數(shù)據(jù)研究院有限公司，貴陽 550000；2. 提升政府治理能力大數(shù)據(jù)應(yīng)用技術(shù)國家工程實(shí)驗(yàn)室，貴陽 550000；3. 中國科學(xué)院水利部成都山地災(zāi)害與環(huán)境研究所，成都 610041；4. 中國農(nóng)業(yè)科學(xué)院農(nóng)業(yè)資源與農(nóng)業(yè)區(qū)劃研究所，北京 100081）

復(fù)種指數(shù)是表征耕地利用程度的重要參數(shù)。然而，傳統(tǒng)方法存在對影像獲取條件要求較高，或在地表復(fù)雜區(qū)域提取精度較低等問題。高時(shí)空分辨率數(shù)據(jù)融合算法（如Spatial and Temporal Adaptive Reflectance Fusion Model, STARFM）能有效地結(jié)合不同數(shù)據(jù)的優(yōu)勢，有望被應(yīng)用于提高復(fù)雜地表區(qū)域復(fù)種指數(shù)的提取精度。該研究以Landsat TM（Thematic Mapper）及MODIS（Moderate-resolution Imaging Spectroradiometer）為數(shù)據(jù)源，基于STARFM模型，構(gòu)建了川東丘陵某區(qū)域內(nèi)2010－2011年的Landsat-like 時(shí)序NDVI（Normalized Difference Vegetation Index）數(shù)據(jù)集，進(jìn)而提取了該區(qū)域2010年冬季作物種植區(qū)及鹽亭縣2011年耕地復(fù)種指數(shù)的空間分布情況。利用目視解譯樣點(diǎn)（1509個(gè)）驗(yàn)證及多尺度（30～4 000 m）驗(yàn)證方法，對不同方法提取的2010年冬季作物種植區(qū)進(jìn)行了對比分析。結(jié)果表明：1）在30m空間尺度上，基于Landsat影像分類法的總體驗(yàn)證精度為89.73%，高于基于Landsat-like時(shí)序NDVI峰值法的54.94%；2）在250～4 000 m空間尺度上，基于Landsat-like時(shí)序NDVI峰值法的總體驗(yàn)證精度比基于MODIS時(shí)序NDVI峰值法高3%～7%。利用統(tǒng)計(jì)年鑒及調(diào)查樣點(diǎn)（73個(gè)）數(shù)據(jù)，對基于新方法提取的鹽亭縣2011年耕地復(fù)種指數(shù)結(jié)果進(jìn)行了驗(yàn)證，在縣域尺度上其與統(tǒng)計(jì)數(shù)據(jù)非常接近；其與調(diào)查樣點(diǎn)的總體驗(yàn)證精度達(dá)到73.97%。綜上，基于數(shù)據(jù)融合算法提高數(shù)據(jù)源空間分辨率的方式，不僅能夠提高復(fù)雜地表復(fù)種指數(shù)結(jié)果的空間精細(xì)程度和提取精度，而且在實(shí)際應(yīng)用中也有很好的實(shí)用性。

遙感；NDVI；時(shí)空數(shù)據(jù)融合；復(fù)種指數(shù)；STARFM；地表復(fù)雜區(qū)域

0 引 言

復(fù)種指數(shù)指某耕作單元在單位時(shí)間內(nèi)的耕作次數(shù)[1]。弄清某區(qū)域的耕地復(fù)種信息，對于該區(qū)域耕地潛力挖掘、種植制度優(yōu)化、糧食產(chǎn)量估算及作物固碳研究等具有重要的意義。耕地的復(fù)種信息通常是通過統(tǒng)計(jì)匯總的方法得到的，然而這種方法通常需要耗費(fèi)大量人力和時(shí)間。同時(shí)，獲取的結(jié)果是以各級行政區(qū)為單位的統(tǒng)計(jì)數(shù)據(jù)，不僅展現(xiàn)形式不直觀，而且空間精細(xì)程度也較差。基于遙感技術(shù)獲取復(fù)種信息的方法因具備客觀、高時(shí)效和低成本等優(yōu)點(diǎn)，正逐漸成為該領(lǐng)域的普遍研究手段。

基于不同的數(shù)據(jù)源及監(jiān)測原理，不同學(xué)者已提出了多種不同的復(fù)種指數(shù)遙感監(jiān)測方法?；诟呖臻g分辨率的遙感影像（如Landsat），可以通過觀測生長季內(nèi)種植作物像元與休耕像元之間的光譜差異，然后進(jìn)行分類，得到耕地的復(fù)種情況[2]。這類方法原理簡單，且運(yùn)算速度快。已有研究表明，在數(shù)據(jù)質(zhì)量較高的情況下，對于小農(nóng)耕作區(qū)域（小于2 hm2），基于Landsat閾值分類法較常用的基于MODIS時(shí)序NDVI峰值法能夠得到更高的復(fù)種指數(shù)提取精度[3]。但是，該方法對于高空間分辨率遙感影像的獲取條件要求較高，需同時(shí)滿足3個(gè)條件：在生長季內(nèi)，且作物的綠葉在光譜上能夠明顯被識別；每個(gè)生長季至少一期；無云霧等的影響。上述條件通常難以同時(shí)滿足，尤其對于復(fù)種指數(shù)較高的區(qū)域，從而在一定程度上阻礙了此類方法的應(yīng)用[3]。

基于高時(shí)間分辨率的植被指數(shù)數(shù)據(jù)提取復(fù)種指數(shù)信息是更普遍的方法，且已被廣泛應(yīng)用于大區(qū)域尺度的復(fù)種指數(shù)研究中[4-7]。早期有基于重構(gòu)時(shí)序植被指數(shù)的峰值點(diǎn)探測法及基于重構(gòu)時(shí)序植被指數(shù)的曲線匹配度檢驗(yàn)法。前者利用各種不同的算法計(jì)算出重構(gòu)后時(shí)序植被指數(shù)曲線的峰值點(diǎn)個(gè)數(shù)，從而反推該耕地像元的復(fù)種信息[8-10]；后者則是通過基于不同的算法測算與復(fù)種信息已知的時(shí)序植被指數(shù)曲線（標(biāo)準(zhǔn)曲線）之間的匹配程度，并認(rèn)為某耕地像元的復(fù)種信息與匹配度最高的標(biāo)準(zhǔn)曲線復(fù)種信息相同，以獲取該耕地像元的復(fù)種信息[11]。除這2種較常用的方法之外，近年來還發(fā)展了一些新的方法。如基于重構(gòu)時(shí)序植被指數(shù)的時(shí)域混合模型分解法[12]，回歸模型法[3]及連續(xù)小波變換法[13]等。盡管這些基于高時(shí)間分辨率植被指數(shù)數(shù)據(jù)的算法都已被證實(shí)在空間異質(zhì)性較低的區(qū)域具有較高的驗(yàn)證精度。但是，由于常用高時(shí)間分辨率數(shù)據(jù)源的空間分辨率較低（小于250 m），這些方法目前仍難以有效地應(yīng)用于復(fù)雜地表區(qū)域[14]（如川東丘陵區(qū)），因?yàn)樵谶@些區(qū)域耕地的空間尺度較小，存在較多混合像元的問題[15]。

因此，在復(fù)雜地表區(qū)域應(yīng)用上述2類復(fù)種指數(shù)遙感提取方法均受限于不同數(shù)據(jù)源的時(shí)空特征。基于高空間分辨率數(shù)據(jù)的方法（如Landsat閾值分類法）受限于數(shù)據(jù)源的時(shí)間分辨率較低；而基于高時(shí)間分辨率數(shù)據(jù)的方法（如基于MODIS時(shí)序NDVI峰值法）受限于數(shù)據(jù)源的空間分辨率較低。基于高時(shí)空分辨率數(shù)據(jù)融合算法，構(gòu)建具有“雙高”（高時(shí)間分辨率及高空間分辨率）特征的遙感數(shù)據(jù)，是提高地表復(fù)雜區(qū)域復(fù)種指數(shù)提取精度的一種思路[16]。近年來，已有多種高時(shí)空分辨率數(shù)據(jù)融合算法被提出[17-19]。其中，最常用的是STARFM算法。該算法已被多項(xiàng)研究結(jié)果證實(shí)不僅能夠獲取高精度的時(shí)間序列Landsat-like影像[20]，還能夠通過提高數(shù)據(jù)源的時(shí)間或空間分辨率，在具體應(yīng)用中提高某些遙感產(chǎn)品的提取精度[21]。已有研究基于STARFM算法獲取的高時(shí)空分辨率數(shù)據(jù)集，獲取了更高精度的復(fù)種指數(shù)結(jié)果[22-23]。但這些研究區(qū)域的地表復(fù)雜度較低，難以體現(xiàn)其在高空間分辨率方面的優(yōu)勢。本文以典型的復(fù)雜地表區(qū)域——川東丘陵區(qū)為例，擬通過該算法構(gòu)建高時(shí)空分辨率NDVI數(shù)據(jù)集，以提高復(fù)種指數(shù)的提取精度。

1 研究區(qū)概況

本文選取四川省東部面積大小約8 100 km2（90 km×90 km）的區(qū)域作為研究區(qū)（103°40′28″～106°11′49″E，30°48′36″～32°42′6″N）。該區(qū)域?yàn)榈湫偷那鹆陞^(qū)，海拔大約在300～700 m之間。該區(qū)域處于四川盆地亞熱帶濕潤季風(fēng)氣候區(qū)。降水分干濕兩季。濕季（約5－10月）降水量占全年降水量的80%以上。氣溫四季分明：夏季溫度最高，平均溫度達(dá)25 ℃左右；春秋兩季溫度適中，約16～17 ℃；冬季溫度較低，僅約6 ℃。適宜的氣候條件能夠滿足冬（約為11月至第二年5月）、夏（約為5月至10月）兩季作物的生長。盡管耕地分布廣泛（圖1c），但受地形起伏影響，這些耕地地塊面積大都小于250 m，有的甚至小于30 m，這給利用現(xiàn)有高時(shí)間分辨率植被指數(shù)數(shù)據(jù)（如MODIS NDVI）提取該區(qū)域復(fù)種指數(shù)的相關(guān)方法帶來了困難。

圖1 研究區(qū)地理位置與土地覆蓋類型

2 數(shù)據(jù)及預(yù)處理

2.1 數(shù)據(jù)源

1）Landsat數(shù)據(jù)

選取了該區(qū)域2010年12月8日的Landsat TM影像。該數(shù)據(jù)幾乎沒有受到云雪等的影響，質(zhì)量較好（圖1b）。水體、林地、居民區(qū)、耕地等不同的土地覆蓋類型的光譜特征也十分明顯。在耕地區(qū)域，通過目視解譯，能很容易地區(qū)分種植區(qū)域（圖1b中呈現(xiàn)出粉紅色特征區(qū)域）和非種植區(qū)域。該數(shù)據(jù)有2個(gè)作用：一方面在利用STARFM算法構(gòu)建高時(shí)空分辨率NDVI數(shù)據(jù)集時(shí)，該數(shù)據(jù)可作為STARFM算法的輸入數(shù)據(jù)（基期Landsat影像）；另一方面通過對其進(jìn)行分類，可得到該區(qū)域2010年冬季作物的種植分布情況，用于提取結(jié)果的驗(yàn)證。

2）MODIS數(shù)據(jù)

相較于NOAA AVHRR（Advanced Very High Resolution Radiometer）及SPOT VGT（Vegetation）等時(shí)間序列植被指數(shù)數(shù)據(jù)，由于具有較高的空間分辨率，MODIS時(shí)序植被指數(shù)數(shù)據(jù)(即NDVI)成為了近年來該研究領(lǐng)域最為常用的數(shù)據(jù)源[9,12]。在本文的研究中，該數(shù)據(jù)作為高時(shí)間分辨率NDVI數(shù)據(jù)源。為盡量減小云雪等的負(fù)面影響，本文選擇的數(shù)據(jù)集為16 d最優(yōu)觀測值合成數(shù)據(jù)（即MOD13Q1），在時(shí)間上包含了2010年冬季以及2011年全年（即2010年第20期至2011年第23期）。該數(shù)據(jù)具有2個(gè)作用：其一作為基于STARFM算法構(gòu)建高時(shí)空分辨率NDVI數(shù)據(jù)集的輸入數(shù)據(jù)；其二基于該數(shù)據(jù)利用傳統(tǒng)方法提取復(fù)種指數(shù)，并將其與基于數(shù)據(jù)融合算法形成的高時(shí)空分辨率數(shù)據(jù)集的提取結(jié)果進(jìn)行對比分析。

3）數(shù)據(jù)預(yù)處理

本文使用MODIS重投影工具M(jìn)ODIS Reprojection Tool（MRT）對MOD13Q1數(shù)據(jù)做了投影轉(zhuǎn)換、重采樣和裁剪等預(yù)處理，得到UTM投影，像元大小為250 m的數(shù)據(jù)。為進(jìn)一步降低云雪等的負(fù)面影響，本文應(yīng)用SG 濾波算法[24-25]對預(yù)處理后的時(shí)序MOD13Q1 NDVI數(shù)據(jù)進(jìn)行了平滑處理。另外，STARFM算法要求輸入的Landsat與MODIS數(shù)據(jù)具有相同的投影和像元大小。因此，本文還將平滑后的250m時(shí)序MODIS NDVI數(shù)據(jù)重采樣為30m像元大小。對于Landsat TM數(shù)據(jù)，本文采用了Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS)[26]數(shù)據(jù)預(yù)處理程序，對其進(jìn)行輻射定標(biāo)和大氣校正等預(yù)處理，最終得到地表反射率和NDVI數(shù)據(jù)。

2.2 驗(yàn)證數(shù)據(jù)

首先，利用該區(qū)域2010年土地覆被分類數(shù)據(jù)（圖1c）生成一個(gè)耕地的掩膜文件，將2010年12月8日Landsat TM影像上的非耕地部分去除；然后，在去除非耕地后的TM影像上，隨機(jī)選取1509個(gè)樣點(diǎn)（圖2，含未種植樣本點(diǎn)627個(gè)；有作物種植樣本點(diǎn)882個(gè)）；最后，對每個(gè)樣點(diǎn)逐個(gè)進(jìn)行目視判定，得到用于驗(yàn)證的樣本數(shù)據(jù)。利用該樣本數(shù)據(jù)，可以在30 m空間尺度上對基于Landsat TM數(shù)據(jù)分類和基于STARFM算法形成的Landsat-like 時(shí)序NDVI數(shù)據(jù)集得到的作物種植區(qū)結(jié)果進(jìn)行精度驗(yàn)證。

圖2 隨機(jī)樣本點(diǎn)空間分布

由于Landsat-like 時(shí)序NDVI數(shù)據(jù)（分辨率30 m）與MODIS時(shí)序NDVI數(shù)據(jù)（分辨率250 m）的空間分辨率不同，本文采用空間聚合的方法，在多個(gè)空間尺度下（250～4 000 m）對比2種方法的提取精度。通過對2010年12月8日Landsat影像進(jìn)行分類，得到2010年冬季作物種植區(qū)空間分布圖，并將此數(shù)據(jù)作為多尺度驗(yàn)證的參考數(shù)據(jù)。

2013年8月上旬，在鹽亭縣境內(nèi)的搜集了17塊樣地的耕地復(fù)種信息。調(diào)查樣地的空間分布如圖3所示。盡管由于交通不便，在空間上樣地未能在整個(gè)研究區(qū)域呈均勻分布，但是從地形上看，從起伏較小到起伏較大的區(qū)域都有樣地分布。因此，調(diào)查的樣地具有一定的空間代表性。另外，由于每塊樣地有3～6個(gè)不等的樣點(diǎn)，各樣點(diǎn)之間相距均在30 m以上。因此，總計(jì)可用于分析的樣點(diǎn)數(shù)量大于樣地?cái)?shù)量，為73個(gè)。利用該數(shù)據(jù)對基于新方法提取的鹽亭縣2011年復(fù)種指數(shù)結(jié)果進(jìn)行逐樣點(diǎn)對比，從而驗(yàn)證新方法的實(shí)用性。

圖3 四川省鹽亭縣地理位置示意圖及野外調(diào)查樣地空間分布

2.3 輔助數(shù)據(jù)

其他的輔助數(shù)據(jù)包括該區(qū)域2010年土地覆被分類數(shù)據(jù)（圖1c）。該數(shù)據(jù)來源于中國科學(xué)院水利部成都山地災(zāi)害與環(huán)境研究所數(shù)字山地與遙感應(yīng)用中心。數(shù)據(jù)的生產(chǎn)基于中國環(huán)境減災(zāi)衛(wèi)星（HJ）影像并輔助Landsat TM影像作為數(shù)據(jù)源，采用了面向?qū)ο笈c決策樹分類器相結(jié)合的方式[27]，并經(jīng)過多輪人工目視修正，共提取了38類土地覆被類型，經(jīng)獨(dú)立樣點(diǎn)驗(yàn)證，該數(shù)據(jù)一級類精度超過95%，二級類精度超過85%，滿足本研究對土地覆被數(shù)據(jù)精度的要求。

3 研究方法

該區(qū)域?yàn)橐荒陜墒旆N植區(qū)。冬季作物生長季通常為11月至第二年5月；夏季作物生長季通常為6－10月。在夏季，該區(qū)域云霧等較多，導(dǎo)致難以獲取高質(zhì)量Landsat影像。此外，夏季作物的光譜與林地及灌木等難以區(qū)分。因此，本文在驗(yàn)證精度和實(shí)用性時(shí)采用了不同的策略。精度驗(yàn)證時(shí)僅使用了2010年冬季的數(shù)據(jù)對不同方法提取的作物種植區(qū)結(jié)果進(jìn)行驗(yàn)證。同時(shí)，為驗(yàn)證方法的實(shí)用性，利用本文方法提取了綿陽市鹽亭縣2011年全年的復(fù)種指數(shù)空間分布情況，并利用統(tǒng)計(jì)年鑒及2013年在部分區(qū)域的實(shí)地調(diào)查數(shù)據(jù)，對提取結(jié)果進(jìn)行了驗(yàn)證。

本文使用3種不同的方法獲取研究區(qū)域內(nèi)2010年冬季的作物種植區(qū)域。如圖4所示，第1種方法為Landsat影像分類法：該方法的結(jié)果作為多尺度驗(yàn)證的參考數(shù)據(jù)。具體實(shí)現(xiàn)過程如下：首先，對去除非耕地后的TM影像進(jìn)行多分辨率分割[28]；然后，在分割后的影像上人工選取種植作物耕地與閑置耕地的分割單元作為樣本；最后，利用選取的樣本對去除非耕地后的TM影像進(jìn)行分類，得到分類圖。利用選取的樣本數(shù)據(jù)進(jìn)行驗(yàn)證，該分類圖的總體精度為：89.73%（表1）。第2種為基于Landsat-like時(shí)序NDVI峰值法（圖4）：該方法與傳統(tǒng)方法的不同之處在于，其需要利用STARFM算法融合基期Landsat TM NDVI（如2010年12月8日）與濾波后的MODIS NDVI數(shù)據(jù)（即2010年第20期至2011年第9期），構(gòu)建時(shí)序Landsat-like NDVI數(shù)據(jù)集；然后基于該數(shù)據(jù)集設(shè)定合理的閾值（閾值1），以獲取作物種植區(qū)。第3種方法為常用的基于MODIS時(shí)序NDVI峰值法（圖4）：首先，基于濾波后的時(shí)序MODIS NDVI數(shù)據(jù)，利用二次差分算法（見3.5部分）獲取各像元可能的波峰；然后，設(shè)定合理的閾值（閾值2）剔除其中的“偽波峰”；最后，統(tǒng)計(jì)“真波峰”的個(gè)數(shù)即可得到相應(yīng)的種植區(qū)域。

圖4 本研究中使用的3種作物種植區(qū)提取方法

3.1 STARFM算法

假定經(jīng)幾何精校正、坐標(biāo)轉(zhuǎn)換及重采樣等預(yù)處理后，某t時(shí)刻MODIS地表反射率影像與該時(shí)刻TM地表反射率影像具有同樣大小的像元尺度及坐標(biāo)系統(tǒng)。由于MODIS地表反射率與TM地表反射率具有很好的一致性[29]。若某時(shí)刻（0）MODIS像元地表反射率(0)已知，TM地表反射率(0)未知；而另一時(shí)刻（t）MODIS地表反射率(t)與TM地表反射率(t)均已知。那么，0時(shí)刻未知的TM地表反射率可由下式求得[17]：

式中表示搜索窗口的大?。ū疚脑O(shè)置為3 000 m×大3 000 m）；（x,y）表示某像元的位置；W為某鄰近光譜相似像元的權(quán)重值，因考慮混合像元、地表覆被或物候改變及BRDF效應(yīng)等的影響而引入鄰近光譜相似像元信息而設(shè)置。該權(quán)重取決于以下3個(gè)因素：1）由其與對應(yīng)MODIS像元反射率之間的差異所表示的光譜距離(S)；2）其與中心像元之間的空間距離(d)；3）由與之對應(yīng)的兩時(shí)刻（0及t）MODIS像元反射率差異所表示的時(shí)間距離（T）。而對鄰近相似像元的篩選條件是：其必須比中心像元能提供更多的光譜信息。另外，考慮Landsat與MODIS地表反射率數(shù)據(jù)在預(yù)處理過程中存在一定的不確定性，最終的篩選條件還加入了不確定性參數(shù)（本文設(shè)置為0.025）。有關(guān)各參數(shù)的具體計(jì)算方式及算法原理的更詳細(xì)描述請參考原始文獻(xiàn)[17]。

3.2 時(shí)序NDVI構(gòu)建方案

STARFM算法為單波段的運(yùn)算模式[17]。因此，存在2種不同的方案可用于構(gòu)建高時(shí)空分辨率NDVI數(shù)據(jù)集。方案1：首先基于STARFM算法逐波段地融合基期影像的地表反射率，生成高時(shí)空分辨率的紅波段與近紅外波段地表反射率；然后通過波段之間的代數(shù)運(yùn)算間接得到預(yù)測期的NDVI數(shù)據(jù)。方案2：首先基于NDVI與紅及近紅外波段之間的轉(zhuǎn)換關(guān)系，分別計(jì)算出基期影像的NDVI；然后應(yīng)用STARFM算法直接融合NDVI得到預(yù)測期NDVI數(shù)據(jù)。對于高時(shí)空分辨率Landsat-like NDVI影像的預(yù)測精度而言，方案2優(yōu)于方案1[30-31]；并且在復(fù)種指數(shù)提取的應(yīng)用中，2種方案的結(jié)果高度相似[31]。而由于僅運(yùn)算NDVI一個(gè)波段，方案2所需的存儲空間及運(yùn)算時(shí)間相對方案1更少。因此，在復(fù)種指數(shù)遙感反演的應(yīng)用中，選擇方案2構(gòu)建高時(shí)空分辨率NDVI數(shù)據(jù)集是更合理的選擇[31]。

3.3 樣本點(diǎn)驗(yàn)證

基于隨機(jī)選取的1509個(gè)樣本點(diǎn)，逐個(gè)地對基于Landsat影像分類法和基于Landsat-like時(shí)序NDVI峰值法提取的結(jié)果進(jìn)行對比驗(yàn)證，并將驗(yàn)證結(jié)果形成混淆矩陣?；诨煜仃?，計(jì)算生產(chǎn)者精度、用戶精度和總體精度，以及Kappa系數(shù)等值，定量地對2種方法提取的冬季作物種植區(qū)結(jié)果進(jìn)行精度驗(yàn)證。

3.4 多尺度驗(yàn)證

鑒于Landsat、Landsat-like與MODIS數(shù)據(jù)的空間分辨率不同，本文采用了多尺度驗(yàn)證的方法對作物種植區(qū)提取結(jié)果進(jìn)行驗(yàn)證。首先，將3種方法得到的提取結(jié)果分別聚合至250～4 000 m的不同尺度，然后，以Landsat影像分類法提取的結(jié)果為驗(yàn)證參考數(shù)據(jù)，在不同的空間尺度上，對比基于Landsat-like時(shí)序NDVI峰值法及基于MODIS時(shí)序NDVI峰值法的提取結(jié)果。另外，為更直觀地表達(dá)不同方法提取作物種植區(qū)的對比效果，本文引入了重合率的概念，即在不同尺度上，基于不同方法的提取結(jié)果與驗(yàn)證數(shù)據(jù)結(jié)果相同的像元個(gè)數(shù)占像元總數(shù)的百分比。例如：首先將不同方法的提取結(jié)果和驗(yàn)證數(shù)據(jù)均聚合到250 m尺度；再基于“贏者通吃”的策略（即聚合后大于等于0.5為1；小于0.5則為0）將聚合后的結(jié)果轉(zhuǎn)換為用“0”和“1”分別表示休耕與種植的“二值”影像；然后，逐像元地判定不同方法轉(zhuǎn)換后的結(jié)果與驗(yàn)證數(shù)據(jù)之間的重合率，進(jìn)而對比兩者的精度。

3.5 復(fù)種指數(shù)提取

基于時(shí)間序列NDVI數(shù)據(jù)集，本文選擇較為常用的“二次差分”算法提取可能的波峰[8]：首先，將一年內(nèi)重構(gòu)后時(shí)間序列植被指數(shù)的個(gè)值按時(shí)間順序形成數(shù)列；然后，依次用后面的植被指數(shù)值減去前面與其相鄰的植被指數(shù)值，形成-1個(gè)新值；之后，對這-1個(gè)新值按以下規(guī)則重新賦值：如果是正數(shù)則為1，反之為?1；最后，對新賦值的-1個(gè)值按上面的方法再進(jìn)行一次差分運(yùn)算，得到-2個(gè)由?2、0、2組成的數(shù)列。其中，元素為?2且前后元素皆為零的點(diǎn)就是可能的波峰。

波峰閾值及作物物候等信息被認(rèn)為能有效地用于“偽峰值點(diǎn)”的剔除[32-33]。本文參考了近10 a該區(qū)域鄰近物候觀測站點(diǎn)的數(shù)據(jù)確定了冬季和夏季作物的物候規(guī)律，判定了2個(gè)峰值點(diǎn)可能出現(xiàn)的時(shí)期大致分別為每年的2－3月和7－8月（即2011年第4～7期和14～17期）；采用了目視觀察及統(tǒng)計(jì)判定的方法確定了時(shí)序NDVI的波峰閾值：首先參考Landsat影像分類結(jié)果圖，選取純耕地像元；再統(tǒng)計(jì)這些耕地像元時(shí)序NDVI的波動情況，以確定兩種時(shí)序NDVI數(shù)據(jù)集（即時(shí)序Landsat-like NDVI數(shù)據(jù)與MODIS NDVI數(shù)據(jù)）的波峰閾值。

4 結(jié)果與分析

4.1 波峰閾值

2種數(shù)據(jù)集中，耕地像元時(shí)序NDVI波動的值域不同。就2010年冬季生長季的耕地而言，Landsat-like時(shí)序NDVI數(shù)據(jù)各期均值為0.43～0.56；而MODIS時(shí)序NDVI數(shù)據(jù)各期均值為0.54～0.65。總體而言，后者平均高于前者約0.11。這一方面可能與MODIS及TM傳感器的紅及近紅外波段的波譜設(shè)置有關(guān)，MODIS傳感器的波譜設(shè)置較窄，其對植被信息則更為敏感，NDVI值也相對較高；另一方面可能由于MODIS像元較大，更容易受到鄰近常綠森林或灌木等植被信息的影響，從而導(dǎo)致NDVI值偏高。兩序列NDVI數(shù)據(jù)均在2011年3月（2011年第6期）達(dá)到峰值。但峰值相差近0.10。因此，在應(yīng)用二次差分算法提取兩種數(shù)據(jù)的可能波峰后，需利用不同的波峰閾值對其進(jìn)行篩選。設(shè)置該閾值，既要保證不會將真正的峰值“濾除”，還需考慮盡量消除雜草的影響。本文中選取的樣本數(shù)據(jù)可以看出，2種數(shù)據(jù)在峰值點(diǎn)處的均值分別為0.65±0.06和0.55±0.07（圖5）。該均值減去一倍標(biāo)準(zhǔn)差后的值分別為0.59及0.48。綜合考慮，本文將兩個(gè)峰值點(diǎn)的閾值分別設(shè)定為0.55與0.45。

注：2010年冬季包含MOD13Q1數(shù)據(jù)時(shí)間軸2010年第21期－2011年第9期。

4.2 目視評價(jià)

不同方法提取的研究區(qū)2010年冬季耕地種植分布情況如圖6。圖6a是基于該區(qū)域2010年12月8日的Landsat TM影像進(jìn)行分類得到的結(jié)果?；跁r(shí)序NDVI數(shù)據(jù)的2種方法得到的結(jié)果與基于Landsat影像分類法的有較大差異（圖6a與圖6b；圖6a與圖6c）。2種峰值法的結(jié)果顯示：在該區(qū)域的左部和下部，大部分耕地均未耕種，而在右部，大部分區(qū)域有作物種植（圖6b和圖6c）；但分類法的結(jié)果顯示：閑置的耕地與種植作物的耕地分布比較隨機(jī)和均勻（圖6a）。后者的結(jié)果明顯更加符合實(shí)際情況。出現(xiàn)這種異常分布的原因可能是用于提取復(fù)種指數(shù)時(shí)的數(shù)據(jù)源受到了“污染”。由于原始數(shù)據(jù)經(jīng)過了濾波處理，因此，“污染”來自兩方面：其一是云霧及其陰影；其二是SG濾波算法。2種峰值法的結(jié)果在空間上表現(xiàn)出較高的整體一致性（圖6b與圖6c），則進(jìn)一步印證了這種可能性。進(jìn)一步的檢視中，本文發(fā)現(xiàn)在某些區(qū)域，基于Landsat-like時(shí)序NDVI峰值法的提取結(jié)果不僅精細(xì)度較高，而且分布情況合理，在空間上也能較好地與參考數(shù)據(jù)吻合（圖6a與圖6b）。

4.3 樣本點(diǎn)驗(yàn)證

利用1509個(gè)樣本點(diǎn)，對基于Landsat影像分類法和Landsat-like時(shí)序NDVI峰值法提取的作物種植區(qū)結(jié)果進(jìn)行了精度驗(yàn)證（表1）。總體而言，兩種方法提取的結(jié)果與目視評價(jià)的情況基本一致?；贚andsat影像分類法取得了很高的驗(yàn)證精度（89.73%），一方面因?yàn)橛懈哔|(zhì)量的Landsat數(shù)據(jù)；另一方面使用了高精度土地利用覆被圖對非耕地區(qū)域進(jìn)行了掩膜處理。高精度的提取結(jié)果也為多尺度驗(yàn)證提供了高質(zhì)量的參考數(shù)據(jù)。而Landsat-like時(shí)序NDVI峰值法的驗(yàn)證精度較低（54.94%）。主要原因是用于提取的時(shí)序NDVI數(shù)據(jù)受到了“污染”，未能表現(xiàn)出正常的波峰，從而導(dǎo)致利用峰值法不能正確提取出作物種植信息。這從目視評價(jià)的結(jié)果中也能得到印證（圖6）。另外，由于數(shù)據(jù)質(zhì)量較差，較多實(shí)際有種植像元被錯(cuò)誤地判定為休耕，也直接導(dǎo)致了對于休耕像元提取精度要明顯低于種植像元。

注：空白區(qū)域?yàn)榉歉亍?/p>

4.4 多尺度驗(yàn)證

表2統(tǒng)計(jì)了不同空間尺度下2種方法的提取結(jié)果與驗(yàn)證參考數(shù)據(jù)之間的重合率。與目視評價(jià)及樣本點(diǎn)驗(yàn)證結(jié)果相似，在30m空間尺度上，基于Landsat-like時(shí)序NDVI峰值法的提取結(jié)果與基于Landsat影像分類法的提取結(jié)果之間的重合率較低（0.57）。類似地，有較多的種植像元被判定為休耕（0.28），是致使重合率較低的主要原因（表1和表2）。這也進(jìn)一步印證了前文的分析：基于Landsat-like時(shí)序NDVI數(shù)據(jù)提取種植區(qū)時(shí)，當(dāng)數(shù)據(jù)受到了“污染”，無法表現(xiàn)出正常的波峰，將導(dǎo)致不能正確地提取出種植信息。當(dāng)然，如果將2種方法提取的結(jié)果聚合到較低的空間尺度，能夠在一定程度上提高提取結(jié)果的精度。當(dāng)空間尺度下降到4 000 m時(shí)，兩者的總體重合率分別提高了0.17和0.13。

表1 基于Landsat影像分類法和Landsat-like時(shí)序NDVI峰值法的2010年冬季作物種植區(qū)提取精度對比

注：數(shù)字上標(biāo)a代表某種類別像元的個(gè)數(shù)。例如，505表示基于Landsat影像分類法提取為休耕像元，驗(yàn)證樣點(diǎn)也為休耕像元的個(gè)數(shù)。122表示基于Landsat影像分類法提取為休耕像元，但是驗(yàn)證樣點(diǎn)為種植像元的個(gè)數(shù)。

Note: The values labeled with “a” are the number of pixels judged as one certain category. For example, 505 is the number of the pixels judged as “fallow” based on the Landsat data, which also verified as “fallow” by the validation samples. 122 is the number of the pixels judged as “fallow” based on the Landsat data, yet verified as “planting” by the validation samples.

表2 基于Landsat-like和MODIS時(shí)序NDVI峰值法提取的2010年冬季作物種植區(qū)結(jié)果在不同空間尺度上的重合率對比

注：a列表示對比驗(yàn)證參考數(shù)據(jù)，基于某方法提取的不同類別像元的重合率。如“種植-種植”和“種植-休耕”表示在參考影像上的種植區(qū)域內(nèi)，基于某方法判定結(jié)果則分別為“種植”像元、“休耕”像元。余同?？傮w重合率表示基于某方法判定正確的所有像元(含種植與休耕)所占的比例。如在30 m空間尺度上，基于Landsat-like時(shí)序NDVI峰值法的總體重合率為0.57。其值等于“種植-種植”的重合率0.30，加上“休耕-休耕”的重合率0.27。

Note: The column “category” indicates that the coincidence rate of different categories of pixels based on one method, compared with the reference map. e.g., “Planting- Planting” and “Planting- Fallow” mean one pixel determined as planting pixels, and one pixel yet determined as fallow pixels based on one method, in the planting area of the reference map, respectively. Same as below. The overall coincidence rate is the proportion of all pixels judged correctly, planting and fallow pixels included. e.g., the value of 0.57, which is the overall coincidence rate by the peak-counting method based on the Landsat-like time-series NDVI dataset at 30 m spatial scales, is exactly equal to the sum of the coincidence rates of “Planting- Planting” and “Fallow- Fallow” (i.e., 0.57=0.30+0.27).

對比2種峰值法提取結(jié)果，基于Landsat-like時(shí)序NDVI峰值法的總體重合率在各個(gè)空間尺度上均高于基于MODIS時(shí)序NDVI峰值法（表2）。這與2種數(shù)據(jù)集之間的空間尺度不確定性有關(guān)[23]。由于四川丘陵地區(qū)的地表異常復(fù)雜，MODIS數(shù)據(jù)的空間分辨率較低，在該區(qū)域內(nèi)存在嚴(yán)重的混合像元問題，難以獲取高精度的種植信息。通過數(shù)據(jù)融合算法，將MODIS時(shí)序NDVI數(shù)據(jù)的空間分辨率提高，能夠獲取更精細(xì)的作物種植分布信息（圖6），從而提高了提取精度。另外，前者對于種植像元的判定精度要低于后者；而前者對于休耕像元的判定精度卻明顯高于后者。MODIS數(shù)據(jù)較低的空間分辨率降低了像元內(nèi)的空間異質(zhì)性，也使得時(shí)間序列NDVI的波峰形態(tài)更趨于“平均”：對于作物種植的像元，其時(shí)序NDVI波峰將趨于“平坦”。這就導(dǎo)致基于其判定為“有作物種植”的結(jié)果更為可靠。而對于無作物種植的像元，其時(shí)序NDVI波峰將被“抬升”。這就導(dǎo)致基于其判定為“無作物種植”的結(jié)果較為不可靠。

4.5 復(fù)種指數(shù)結(jié)果

基于Landsat-like時(shí)序NDVI峰值法提取了鹽亭縣2011年全年復(fù)種指數(shù)的空間分布情況（圖7）?？梢钥闯?，鹽亭縣2011年耕地仍以一年兩熟為主，且集中分布于北部和東部。在西部和南部區(qū)域，有大量耕地為一年一熟。統(tǒng)計(jì)年鑒結(jié)果顯示，鹽亭縣2011年耕地的復(fù)種指數(shù)為1.69。而本文的方法得到的該縣復(fù)種指數(shù)結(jié)果為1.67。因此，在縣域尺度上，基于Landsat-like時(shí)序NDVI峰值法提取的復(fù)種指數(shù)結(jié)果具有很高的精度。

利用2013年8月中旬的野外調(diào)查樣點(diǎn)，對基于Landsat-like時(shí)序NDVI峰值法的提取結(jié)果進(jìn)行了逐點(diǎn)驗(yàn)證（表3）。提取結(jié)果與實(shí)地調(diào)查樣本的總體驗(yàn)證精度達(dá)到了73.97%。與本文新方法提取的復(fù)種指數(shù)結(jié)果一致，實(shí)地調(diào)查結(jié)果中多數(shù)樣點(diǎn)也為一年兩熟制。表中數(shù)據(jù)顯示，有較多的一年兩熟樣點(diǎn)被判定為一年一熟。這是因?yàn)樵搮^(qū)域受云霧干擾較為嚴(yán)重，部分耕地像元的時(shí)序NDVI數(shù)據(jù)的波峰未被“探測”出來。

圖7 四川省鹽亭縣2011年復(fù)種指數(shù)空間分布

表3 基于Landsat-like時(shí)序NDVI峰值法提取的鹽亭縣2011年復(fù)種指數(shù)精度驗(yàn)證

注：上標(biāo)a的值代表對某種類別像元進(jìn)行驗(yàn)證后為某種類別像元的個(gè)數(shù)。

Note: The values labeled with “a” are the number of pixels judged as one certain category, and verified as the other category by the validation samples.

5 討 論

5.1 數(shù)據(jù)融合算法

不同學(xué)者基于不同的理論基礎(chǔ)提出了多種高時(shí)空分辨率數(shù)據(jù)融合算法。除少數(shù)基于半物理模型[34]、數(shù)據(jù)同化[35]、貝葉斯理論[36]、及神經(jīng)網(wǎng)絡(luò)[37]提出的幾種融合算法，大部分算法可歸為3類：第1類為基于混合像元分解算法[18,38]；第2類為基于學(xué)習(xí)的算法[19,39]；第3類即為STARFM系列算法[17,40-41]。相比較而言，前者能夠更好地保留低空間分辨率影像上的多光譜特性；而STARFM系列算法能夠使預(yù)測的影像較大限度地貼近真實(shí)影像的空間分辨率[42]。對于復(fù)雜地表區(qū)域，在需要構(gòu)建高時(shí)空分辨率NDVI數(shù)據(jù)集時(shí)（如復(fù)種指數(shù)或作物物候研究），由于涉及到的波段數(shù)目較少，STARFM系列算法具有更廣泛的應(yīng)用前景。

盡管多位學(xué)者對原始算法進(jìn)行了改進(jìn)[40,43]，改進(jìn)后的算法能夠獲取更接近真實(shí)的遙感影像，但在實(shí)際應(yīng)用中，通常不僅僅關(guān)注算法對影像的預(yù)測能力，更加關(guān)注構(gòu)建的高時(shí)空分辨率數(shù)據(jù)集能否有助于提高某種參數(shù)的反演精度。因此，在對STARFM系列算法進(jìn)行選擇時(shí)，應(yīng)該綜合考慮算法的預(yù)測精度，算法的效率，實(shí)現(xiàn)的難易程度，對基期高空間分辨率影像數(shù)量的要求，以及數(shù)據(jù)集中數(shù)據(jù)的一致性等問題。例如：若需要對森林砍伐情況進(jìn)行監(jiān)測，使用STARFM原始算法構(gòu)建數(shù)據(jù)集將會存在較大誤差，因?yàn)樵撍惴▽τ谕蝗话l(fā)生的細(xì)小變化難以捕捉[17,41]；而在分類制圖的研究中，使用單基期Landsat影像構(gòu)建的NDVI數(shù)據(jù)集，或許能有更好的分類精度，因?yàn)楸M管其預(yù)測精度較低[21]，但相較于基于多基期Landsat影像構(gòu)建的NDVI數(shù)據(jù)集而言，該數(shù)據(jù)集具有更好的一致性[21]。在本文的研究中，僅選擇了單基期Landsat TM影像構(gòu)建時(shí)序Landsat-like NDVI數(shù)據(jù)。一方面由于高質(zhì)量Landsat影像較少；另一方面則是為保證數(shù)據(jù)具有更好的一致性，因?yàn)榉逯捣ㄊ峭ㄟ^計(jì)算波峰的個(gè)數(shù)以獲取像元的復(fù)種情況，數(shù)據(jù)不一致則很可能會導(dǎo)致出現(xiàn)更多的“偽波峰”。

5.2 時(shí)間序列數(shù)據(jù)濾波算法

復(fù)種指數(shù)遙感監(jiān)測方法較多。除基于高空間分辨率遙感影像的閾值分類法外，大部分算法是基于高時(shí)間分辨率遙感數(shù)據(jù)而提出。為減弱云霧等的干擾，通常需要應(yīng)用一些濾波算法對高時(shí)間分辨率遙感數(shù)據(jù)集進(jìn)行重構(gòu)。盡管濾波算法較多，但是沒有哪一種算法在任何情況下都能達(dá)到最優(yōu)的重構(gòu)效果[44]。根據(jù)具體的應(yīng)用目的對濾波算法進(jìn)行篩選是較合理的策略。本文應(yīng)用的是SG濾波算法[24]。該算法是一種局部擬合算法，已被證明在數(shù)據(jù)污染較重的時(shí)候，可能會因?yàn)椤斑^度擬合”而導(dǎo)致濾波效果較差[45]。在本文的研究中，基于重構(gòu)時(shí)序NDVI數(shù)據(jù)兩種方法的提取結(jié)果在空間上的分布基本一致，并且都很不均勻（圖6）。原因可能是SG濾波算法不適合受云霧污染較頻繁的四川丘陵區(qū)。在實(shí)際應(yīng)用中，若使用半局部擬合（如logistic函數(shù)擬合法[46]）或全局?jǐn)M合算法（如Harmonic Analysis of NDVI Time-Series, HANTS[47]），濾波效果應(yīng)該會更好。

5.3 復(fù)種指數(shù)提取算法閾值

基于重構(gòu)后的高時(shí)空分辨率植被指數(shù)數(shù)據(jù)，有多種不同的復(fù)種指數(shù)提取算法，如二次差分法[8]；交叉擬合度檢驗(yàn)法[11]；峰值法[9]；基于LMM（linear mixture model）分類法[12]；連續(xù)小波變換法[13]；自適應(yīng)波峰閾值法[48]；迭代滑動窗口法[4]。但是，由于對這些提取算法的對比研究較少，因此，對于這些方法的選擇沒有一定的標(biāo)準(zhǔn)。本文使用的是峰值法中的二次差分算法。該方法原理簡單，且容易實(shí)現(xiàn)。其難點(diǎn)在于對剔除“偽波峰”閾值的設(shè)定。在以往的研究中，通常設(shè)定該值為0.4[49]。但是，鑒于該地區(qū)具有耕地與常綠林地相間分布的特征（圖1c），在250 m空間尺度上，無論是作物像元或是雜草像元都可能會受到常綠林地的影響，導(dǎo)致NDVI波峰值偏高。若將該值簡單地設(shè)定為0.4是不合理的（本文通過統(tǒng)計(jì)分析的方法得出兩種數(shù)據(jù)集波峰閾值分別設(shè)定為0.55和0.45）。這是因?yàn)長andsat-like NDVI影像并不是真實(shí)的觀測值，其在基于STARFM算法合成的過程中，可能會因融合了低空間分辨率NDVI像元的信息而偏高。因此，由于合成的時(shí)序Landsat-like NDVI數(shù)據(jù)集具有與原始時(shí)序低空間分辨率植被指數(shù)數(shù)據(jù)不同的空間和光譜特征，在應(yīng)用不同的復(fù)種指數(shù)提取算法時(shí)，應(yīng)特別注意不能將在原提取算法中使用的參數(shù)或構(gòu)建的模型等簡單地應(yīng)用于該數(shù)據(jù)集，以免造成較大的誤差。

5.4 數(shù)據(jù)融合算法應(yīng)用成效

在地形簡單區(qū)域，基于MODIS時(shí)序NDVI數(shù)據(jù)提取復(fù)種指數(shù)，通常能夠獲得較高的精度（如88.5%[49]；90.1%[12]；97%[48]）。但是，對于地形復(fù)雜區(qū)域或者作物種植結(jié)構(gòu)復(fù)雜區(qū)域而言，由于在MODIS數(shù)據(jù)空間尺度下存在較多混合像元，難以獲取較準(zhǔn)確的復(fù)種指數(shù)信息[23,50]（如表1）。本文基于STARFM數(shù)據(jù)融合算法，通過提高M(jìn)ODIS時(shí)序NDVI數(shù)據(jù)的空間分辨率，將復(fù)種指數(shù)的提取精度相對地提高了約3%～7%（表2）。這可能代表了通過數(shù)據(jù)融合提高數(shù)據(jù)源空間分辨率的方式以提高復(fù)種指數(shù)提取精度的數(shù)量級。因?yàn)樵撝蹬c在其他研究中基于相同思路提高某參數(shù)遙感提取精度的上升幅度相當(dāng)（如2%[21]；5%[20]）。另外，盡管可能受地表異常復(fù)雜及云霧干擾等因素影響，本文基于STARFM數(shù)據(jù)融合算法得到的復(fù)種指數(shù)提取精度較低（54.94%，表1）。但是，在其他相關(guān)研究中，取得了較高的提取精度（如82.7%～89.0%[22]；88.89%23]）。因此，在諸如四川丘陵區(qū)等容易受到云霧干擾且地形復(fù)雜的區(qū)域，數(shù)據(jù)融合算法在提升復(fù)種指數(shù)提取精度方面有廣闊的應(yīng)用前景和提升空間。

5.5 不同區(qū)域?qū)嶋H應(yīng)用策略

基于遙感技術(shù)，已經(jīng)發(fā)展了多種復(fù)種指數(shù)提取方法。在實(shí)際應(yīng)用中，針對不同氣候和地理特點(diǎn)的區(qū)域，應(yīng)當(dāng)采用不同的策略。對于地表簡單區(qū)域，基于MODIS等低空間分辨率時(shí)序NDVI數(shù)據(jù)，利用傳統(tǒng)方法提取復(fù)種指數(shù)就能夠取得較高的精度（如90.1%[12]；97%[48]）。對于地形復(fù)雜，但受云霧干擾較低的區(qū)域，可以結(jié)合同尺度多源衛(wèi)星數(shù)據(jù)（如Landsat and Sentinel-2）[51-52]，或者利用多年時(shí)間跨度的數(shù)據(jù)[53]，以獲取更多可用的高空間分辨率影像。因?yàn)榛跁r(shí)間序列高空間分辨率影像提取復(fù)種指數(shù)，通常能夠達(dá)到很好的效果（如78.88%[23]；93%[51]；81.2%[22]）。對于地形復(fù)雜，且受云霧干擾較嚴(yán)重的區(qū)域，如果能夠在每個(gè)生長季的合適時(shí)間內(nèi)獲取1期高質(zhì)量的高空間分辨率影像，基于該影像進(jìn)行分類，也能獲得較高的提取精度[23]。但是，這通常是難以滿足的。有效地利用數(shù)據(jù)融合算法得到高時(shí)空分辨率數(shù)據(jù)集，不僅可能獲取更精確的復(fù)種指數(shù)信息[22-23]，而且對高空間分辨率影像獲取條件的要求也較容易達(dá)到。此外，在允許的情況下，空間聚合也是提高復(fù)種指數(shù)提取精度的一種實(shí)用方法[3]。在本文的研究中，將30 m尺度的結(jié)果聚合到4 000 m尺度時(shí)，相對精度提高了0.17。

6 結(jié) 論

受財(cái)政和技術(shù)條件的限制，目前單一傳感器獲取的遙感影像僅具有高空間分辨率或者高時(shí)間分辨率的特征。在諸如四川丘陵區(qū)等地表復(fù)雜且常受云霧等干擾的區(qū)域，僅使用具備某一種特征的遙感數(shù)據(jù)獲取耕地復(fù)種信息難以取得良好的效果。如何有效地融合不同特征遙感數(shù)據(jù)的信息，提高復(fù)種指數(shù)的提取精度，是值得研究的一個(gè)課題。本文以川東丘陵區(qū)部分區(qū)域?yàn)檠芯繉ο?，基于STARFM數(shù)據(jù)融合算法構(gòu)建了該區(qū)域2010－2011年高時(shí)間分辨率Landsat-like NDVI數(shù)據(jù)集，并基于其提取了該區(qū)域2010年冬季作物種植區(qū)及鹽亭縣2011年耕地復(fù)種指數(shù)空間分布信息。在30m空間尺度上，利用1509個(gè)目視解譯的隨機(jī)樣點(diǎn)進(jìn)行驗(yàn)證，基于Landsat-like時(shí)序NDVI峰值法的總體精度為54.94%。將結(jié)果聚合到250～4 000 m的不同空間尺度后，與傳統(tǒng)的基于MODIS時(shí)序NDVI峰值法提取結(jié)果進(jìn)行了多尺度對比驗(yàn)證，基于Landsat-like時(shí)序NDVI峰值法的提取精度提高了3%～7%。將基于Landsat-like時(shí)序NDVI峰值法的種植區(qū)提取結(jié)果從30 m聚合到4 000 m空間尺度，總體驗(yàn)證精度提高了0.17。同時(shí)，利用基于Landsat-like時(shí)序NDVI峰值法提取的鹽亭縣2011年復(fù)種指數(shù)結(jié)果，對比了統(tǒng)計(jì)年鑒及實(shí)地調(diào)查結(jié)果（總體精度為73.97%），從而驗(yàn)證了基于Landsat-like時(shí)序NDVI峰值法的實(shí)用性。本文的研究結(jié)果表明：基于數(shù)據(jù)融合算法提高數(shù)據(jù)源的空間分辨率，不僅能夠提高復(fù)雜地表區(qū)域復(fù)種指數(shù)的提取精度，而且具有較好的實(shí)用性。這將有助于解決在類似地表復(fù)雜且受云霧干擾嚴(yán)重的區(qū)域，應(yīng)用傳統(tǒng)復(fù)種指數(shù)提取方法精度較差的問題。

隨著遙感技術(shù)的進(jìn)一步發(fā)展，一些商業(yè)遙感衛(wèi)星星座（如“吉林一號”等）能夠直接獲取高時(shí)空分辨率數(shù)據(jù)集。限于數(shù)據(jù)的預(yù)處理技術(shù)還不成熟，以及數(shù)據(jù)的獲取費(fèi)用較高等，這些數(shù)據(jù)的廣泛使用還需要一段時(shí)間。但是，這些數(shù)據(jù)的應(yīng)用將極大地推動地表復(fù)雜區(qū)域復(fù)種指數(shù)遙感提取研究。未來基于數(shù)據(jù)融合算法獲取復(fù)雜地表區(qū)域復(fù)種指數(shù)方法的研究價(jià)值更多地在于：1）利用歷史數(shù)據(jù)，構(gòu)建長時(shí)間序列高時(shí)空分辨率數(shù)據(jù)集，與直接獲取的高時(shí)空分辨率數(shù)據(jù)集銜接，從而研究某區(qū)域復(fù)種指數(shù)在較長歷史時(shí)間內(nèi)的連續(xù)變化情況；2）利用高時(shí)空分辨率數(shù)據(jù)集在土地覆蓋分類方面的優(yōu)勢，形成高精度土地利用分類圖獲取耕地信息，并基于其獲取高精度的耕地復(fù)種及作物物候等農(nóng)情信息。

[1] Siebert S, Portmann F T, D?ll P. Global patterns of cropland use intensity [J]. Remote Sensing, 2010, 2(7): 1625-1643.

[2] Heller E, Rhemtulla J M, Lele S, et al. Mapping crop types, irrigated areas, and cropping intensities in heterogeneous landscapes of Southern India using multi-temporal medium-resolution imagery: implications for assessing water use in agriculture[J]. Photogrammetric Engineering & Remote Sensing, 2012, 78(8): 815-827.

[3] Jain M, Mondal P, DeFries R S, et al. Mapping cropping intensity of smallholder farms: A comparison of methods using multiple sensors[J]. Remote Sensing of Environment, 2013, 134: 210-223.

[4] Li L, Friedl M A, Xin Q C, et al. Mapping crop cycles in China using MODIS-EVI time series[J]. Remote Sensing, 2014, 6(3): 2473-2493.

[5] Estel S, Kuemmerle T, Levers C, et al. Mapping crop landuse intensity across Europe using MODIS NDVI time series[J]. Environmental Research Letters, 2016, 11(2): 024015.

[6] Qiu B, Lu D, Tang Z, et al. Mapping cropping intensity trends in China during 1982-2013[J]. Applied Geography, 2017, 79: 212-222.

[7] 王學(xué)，李秀彬，談明洪，等. 華北平原 2001-2011年冬小麥播種面積變化遙感監(jiān)測[J].農(nóng)業(yè)工程學(xué)報(bào)，2015，31(8)：190-199.

Wang Xue, Li Xiubin, Tan Minghong, et al. Remote sensing monitoring of changes in winter wheat area in North China Plain from 2001 to 2011[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2015, 31(8): 190-199. (in Chinese with English abstract)

[8] 范錦龍，吳炳方.復(fù)種指數(shù)遙感監(jiān)測方法[J].遙感學(xué)報(bào)，2004，8(6)：628-636.

Fan Jinlong, Wu Bingfang. A methodology for retrieving cropping index from NDVI profile[J]. Journal of Remote Sensing, 2004, 8(6): 628-636. (in Chinese with English abstract)

[9] Sakamoto T, Van N N, Ohno H, et al. Spatio-temporal distribution of rice phenology and cropping systems in the Mekong Delta with special reference to the seasonal water flow of the Mekong and Bassac rivers[J]. Remote Sensing of Environment, 2006, 100(1): 1-16.

[10] 張偉. 川東丘陵區(qū)耕地復(fù)種指數(shù)遙感反演研究：以綿陽市鹽亭縣為例[D]. 北京：中國農(nóng)業(yè)科學(xué)院，2014.

Zhang Wei. Study on the Remotely Sensed Method for Monitoring Multiple Cropping Index in Hilly Region in Eastern Sichuan: A Case Study for Yanting County in Mianyang[D]. Beijing: Chinese Academy of Agricultural Sciences, 2014. (in Chinese with English abstract)

[11] 辜智慧. 中國農(nóng)作物復(fù)種指數(shù)的遙感估算方法研究：基于SPOT/VGT多時(shí)相NDVI遙感數(shù)據(jù)[D].北京：北京師范大學(xué)，2003.

Gu Zhihui. A Study of Calculating Multiple Cropping Index of Crop in China Using SPOT/VGT Muti-temporal NDVI Data[D]. Beijing: Beijing Normal University, 2003. (in Chinese with English abstract)

[12] Chen C F, Son N T, Chang L Y. Monitoring of rice cropping intensity in the upper Mekong Delta, Vietnam using time-series MODIS data[J]. Advances in Space Research, 2012, 49(2): 292-301.

[13] Qiu B, Zhong M, Tang Z, et al. A new methodology to map double-cropping croplands based on continuous wavelet transform[J]. International Journal of Applied Earth Observation and Geoinformation, 2014, 26(1): 97-104.

[14] Yu Q, Hu Q, Van vliet J, et al. GlobeLand30 shows little cropland area loss but greater fragmentation in China[J]. International Journal of Applied Earth Observation and Geoinformation, 2018, 66: 37-45.

[15] Liu W, Dong J, Xiang K, et al. A sub-pixel method for estimating planting fraction of paddy rice in Northeast China[J]. Remote Sensing Environment, 2018, 205: 305-314.

[16] 張偉，李愛農(nóng)，雷光斌. 復(fù)種指數(shù)遙感監(jiān)測研究進(jìn)展[J]. 遙感技術(shù)與應(yīng)用，2015，30(2)：199-208.

Zhang Wei, Li Ainong, Lei Guangbin. Recent progresses in monitoring multiple cropping index by using remote sensing technologies[J]. Remote Sensing Technology and Application, 2015, 30(2): 199-208. (in Chinese with English abstract)

[17] Gao F, Masek J, Schwaller M, et al. On the blending of the Landsat and MODIS surface reflectance: Predicting daily Landsat surface reflectance[J]. IEEE Transactions on Geoscience and Remote Sensing, 2006, 44(8): 2207-2218.

[18] Zhang W, Li A N, Bian J H, et al. An enhanced spatial and temporal data fusion model for fusing Landsat and MODIS surface reflectance to generate high temporal Landsat-like data[J]. Remote Sensing, 2013, 5(10): 5346-5368.

[19] Huang B, Song H. Spatiotemporal re?ectance fusion via sparse representation[J]. IEEE Transactions on Geoscience and Remote Sensing, 2012, 50(10): 3707-3716.

[20] Walker J J, deBeurs K M, Wynne R H. Dryland vegetation phenology across an elevation gradient in Arizona, USA, investigated with fused MODIS and Landsat data [J]. Remote Sensing of Environment, 2014, 144: 85-97.

[21] Watts J D, Powell S L, Lawrence R L, et al. Improved classification of conservation tillage adoption using high temporal and synthetic satellite imagery[J]. Remote Sensing of Environment, 2011, 115(1): 66-75.

[22] Li L, Zhao Y, Fu Y, et al. High resolution mapping of cropping cycles by fusion of Landsat and MODIS data[J]. Remote Sensing, 2017, 9(12): 1232.

[23] Xiang M T, Yu Q Y, Wu W B. From multiple cropping index to multiple cropping frequency: Observing cropland use intensity at a ?ner scale[J]. Ecological Indicators, 2019, 101: 892-903.

[24] Chen J, J?nsson P, Tamura M, et al. A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golay filter[J]. Remote Sensing of Environment, 2004, 91(3): 332-344.

[25] 邊金虎，李愛農(nóng)，宋孟強(qiáng)，等. MODIS植被指數(shù)時(shí)間序列 Savitzky-Golay 濾波算法重構(gòu)[J]. 遙感學(xué)報(bào)，2010，14(4)：725-741.

Bian Jinhu, Li Ainong, Song Mengqiang, et al. Reconstruction of NDVI time-series datasets of MODIS based on Savitzky-Golay filter[J]. Journal of Remote Sensing, 2010, 14(4): 725-741. (in Chinese with English abstract)

[26] Masek J G, Vermote E F, Saleous N E, et al. A Landsat surface reflectance dataset for North America, 1990-2000[J]. IEEE Geoscience and Remote Sensing Letters, 2006, 3(1): 68-72.

[27] 張正健，李愛農(nóng)，雷光斌，等. 基于多尺度分割和決策樹算法的山區(qū)遙感影像變化檢測方法：以四川攀西地區(qū)為例[J]. 生態(tài)學(xué)報(bào)，2014，24(34)：7222-7232.

Zhang Zhengjian, Li Ainong, Lei Gugangbin, et al. Change detection of remote sensing images based on multiscale segmentation and decision tree algorithm over mountainous area: a case study in Panxi region, Sichuan Province[J]. Acta Ecologica Sinica, 2014, 34(24): 7222-7232. (in Chinese with English abstract)

[28] Baatz M, Schape A. Multiresolution segmentation-An optimization approach for high quality multi-scale image segmentation[M]. In Angewandte Geographische Informations-Verarbeitung XII, Wichmann Verlag: Karlsruhe, Germany, 2000: 12-23.

[29] Feng M, Huang C, Channan S, et al. Quality assessment of Landsat surface reflectance products using MODIS data[J]. Computers and Geosciences, 2012, 38(1): 9-22.

[30] Tian F, Wang Y, Fensholt R, et al. Mapping and evaluation of NDVI trends from synthetic time series obtained by blending Landsat and MODIS data around a coalfield on the loess plateau[J]. Remote Sensing, 2013, 5(10): 4255-4279.

[31] Li A N, Zhang W, Lei G B, et al. Comparative analysis on two schemes for synthesizing the high temporal Landsat-like NDVI dataset based on the STARFM algorithm[J]. ISPRS International Journal of Geo-Information, 2015, 4(3): 1423-1441.

[32] 王鵬新，荀蘭，李俐，等. 基于時(shí)間序列葉面積指數(shù)傅里葉變換的農(nóng)作物種植區(qū)域提取[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2017，33(21)：207-215.

Wang Pengxin, Xun Lan, Li Li, et al. Extraction of planting areas of main crops based on fourier transformed characteristics of time series leaf area index products[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(21): 207-215. (in Chinese with English abstract)

[33] 徐昔保，楊桂山.太湖流域1995－2010年耕地復(fù)種指數(shù)時(shí)空變化遙感分析[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2013，29(3)：148-155.

Xu Xibao, Yang Guishan. Spatial and temporal changes of multiple cropping index in 1995-2010 in Taihu Lake basin, China[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2013, 29(3): 148-155. (in Chinese with English abstract)

[34] Roy D P, Ju J, Lewis P, et al. Multi-temporal MODIS- Landsat data fusion for relative radiometric normalization, gap filling, and prediction of Landsat data[J]. Remote Sensing of Environment, 2008, 112(6): 3112-3130.

[35] Cai W, Song J, Wang J, Xiao Z. High spatial-and temporal-resolution NDVI produced by the assimilation of MODIS and HJ-1 data[J]. Canadian Journal of Remote Sensing, 2011, 37(6): 612-627.

[36] Xue J, Leung Y, Fung T. A bayesian data fusion approach to spatio-temporal fusion of remotely sensed images[J]. Remote Sensing, 2017, 9(12): 1310.

[37] Li Y F, Li J, He L, et al. A sensor-bias driven spatio-temporal fusion model based on convolutional neural networks[J]. Science China Information Sciences, 2020, 63(4): 140302.

[38] Lu M, Chen J, Tang H, et al. Land cover change detection by integrating object-based data blending model of Landsat and MODIS[J]. Remote Sensing of Environment, 2016, 184: 374-386.

[39] Zhao Y, Huang B, Song H. A robust adaptive spatial and temporal image fusion model for complex land surface change[J] Remote Sensing of Environment, 2018, 208: 42-62.

[40] Zhu X, Chen J, Gao F, et al. An enhanced spatial and temporal adaptive reflectance fusion model for complex heterogeneous regions[J]. Remote Sensing of Environment, 2010, 114(11): 2610-2623.

[41] Hilker T, Wulder M A, Coops N C, et al. A new data fusion model for high spatial-and temporal-resolution mapping of forest disturbance based on Landsat and MODIS[J]. Remote Sensing of Environment, 2009, 113(8): 1613-1627.

[42] Gevaert C M, Garca-Haro F J. A comparison of STARFM and an unmixing-based algorithm for Landsat and MODIS data fusion[J]. Remote Sensing of Environment, 2015, 156: 34-44.

[43] Xia H, Chen Y, Li Y, et al. Combining kernel-driven and fusion-based methods to generate daily high-spatial- resolution land surface temperatures[J]. Remote Sensing of Environment, 2019, 224: 259-274.

[44] Atkinson P M, Jeganathan C, Dash J, et al. Inter-comparison of four models for smoothing satellite sensor time-series data to estimate vegetation phenology[J]. Remote Sensing of Environment, 2012, 123: 400-417.

[45] J?nsson P, Eklundh L. TIMESAT-a program for analyzing time-series of satellite sensor data[J]. Computers and Geosciences, 2004, 30(8): 833-845.

[46] Zhang X Y, Friedl M A, Schaaf C B, et al. Monitoring vegetation phenology using MODIS[J]. Remote Sensing of Environment, 2003, 84(3): 471-475.

[47] Verhoef W. In Fourier Analysis of temporal NDVI in the Southern African and American Continents[M]. Wageningen (The Netherlands): DLO Winand Staring Centre, 1996.

[48] Fan C, Zheng B, Myint S W, et al. Characterizing changes in cropping patterns using sequential Landsat imagery: An adaptive threshold approach and application to Phoenix, Arizona[J]. International Journal of Remote Sensing, 2014, 35(19/20): 7263-7278.

[49] Galford G L, Mustard J F, Melillo J, et al. Wavelet analysis of MODIS time series to detect expansion and intensification of row-crop agriculture in Brazil[J]. Remote Sensing of environment, 2008, 112(2): 576-587.

[50] Biradar C M, Xiao X. Quantifying the area and spatial distribution of double- and triple-cropping croplands in India with multi-temporal MODIS imagery in 2005[J]. International Journal of Remote Sensing, 2011, 32(2): 367-386.

[51] Liu L, Xiao X M, Qin Y W. Mapping cropping intensity in China using time series Landsat and Sentinel2 images and Google Earth Engine[J]. Remote Sensing of Environment, 2020, 239: 111624.

[52] Liu C, Zhang Q, Tao S Q et al. A new framework to map fine resolution cropping intensity across the globe: Algorithm, validation, and implication[J]. Remote Sensing of Environment, 2020, 251: 112095.

[53] Kontgis C, Schneider A, Ozdogan M. Mapping rice paddy extent and intensi?cation in the Vietnamese Mekong River Delta with dense time stacks of Landsat data[J]. Remote Sensing of Environment, 2015, 169: 255-269.

Improvement of extraction accuracy for cropping intensity in complex surface regions using STARFM

Zhang Wei1,2, Li Wei1,2, Tao Guanhong1,2, Li Ainong3, Qin Zhihao4, Lei Guangbin3, Chen Yixi1,2

(1. China Electronics Technology Group Corporation CETC Big Data Research Institute Co., Ltd, Guiyang 550000, China; 2. National Engineering Laboratory for Big Data Application on Improving Government Governance Capabilities, Guiyang 550000, China;3. Institute of Mountain Hazards and Environment, Chinese Academy of Sciences, Chengdu 610041, China; 4. Institute of Agro-Resources and Regional Planning, Chinese Academy of Agricultural Sciences, Beijing 100081, China)

Cropping intensity (CI) is essentially related to the utilization condition of arable land. A classification method based on high spatial resolution data (TM), and a peak-counting method based on high temporal resolution NDVI data (MODIS NDVI) are mainly used to extract CI in remotely sensed field. However, it is difficult to acquire the high spatial resolution data with high quality characteristics in the most classification method, as three requirements cannot be satisfied concurrently. The first condition is the time requirement, where the high spatial resolution data can be observed in the growing season. The second one is the frequency requirements, where at least one period of data can be obtained in the growing season. The third one is the quality requirement, where the data cannot be polluted by the cloud, mist and other environmental surroundings. Furthermore, the performance was relatively low, when the peak-counting method was applied in some complex surface regions, since the high temporal resolution data has the low spatial resolution (>250 m). A feasible strategy can be made to obtain a high extraction accuracy of CI, in order to improve the spatial resolution of high temporal resolution data in the use of spatiotemporal data fusion (Spatial and Temporal Adaptive Reflectance Fusion Model, STARFM). Taking the hilly region in the eastern Sichuan as a research area, the 2010-2011 time series Landsat-like NDVI dataset with high temporal and spatial resolution was collected to extract the spatial distribution of cropping area for the winter-season in2010 and the CI for summer- and winter-seasons in 2011 in Yanting County, using the STARFM algorithm, Landsat TM and MODIS data. Precision validation between different data in various methods was done at a series of spatial scales (30-4 000 m), referring to the visual interpretation of 1509 random samples, and the spatial distribution of cropping area for winter-season in2010 using Landsat data. The results showed that: 1) Compared to the peak-counting method using Landsat-like time-series NDVI dataset, the overall accuracy of verification was higher for the classification method using Landsat data (89.73% vs. 54.94%), at the scale of 30 m. 2) Compared to the peak-counting method using MODIS time-series NDVI dataset, the overall accuracy of verification was higher 3%-7% for the peak-counting method using Landsat-like time-series NDVI dataset, at the scales from 250 m to 4 000 m. 3) The overall accuracy of verification was higher 0.17, when the data from the new method was spatially aggregated from 30 m to 4 000 m. In the CI result of 2011 in Yanting County, the data in the Mianyang Statistical Yearbook of 2011, and 73 samples from the field survey data in August 2013 were used to validate the practicability of new method. There were very close CI values from the statistical data (1.69) and new method (1.67), indicating that the overall accuracy of verification was 73.97% in the field survey. A high accuracy of extraction and fine spatial details of CI can be acquired, when using the fusion of spatiotemporal data in some complex surface regions. The findings can contribute to the evaluation of food production and carbon sequestration potential in complex surface regions.

remote sensing; NDVI; spatiotemporal data fusion; cropping intensity; STARFM; complex surface regions

張偉，李瑋，陶冠宏，等. 利用STARFM模型提高復(fù)雜地表下復(fù)種指數(shù)遙感提取精度[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2020，36(21)：175-185.doi：10.11975/j.issn.1002-6819.2020.21.021 http://www.tcsae.org

Zhang Wei, Li Wei, Tao Guanhong, et al. Improvement of extraction accuracy for cropping intensity in complex surface regions using STARFM[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(21): 175-185. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2020.21.021 http://www.tcsae.org

2020-07-03

2020-10-04

國家自然科學(xué)基金青年基金項(xiàng)目（41701433）

張偉，遙感算法工程師，主要從事時(shí)空數(shù)據(jù)融合及農(nóng)情參數(shù)反演算法研究。Email：zhang106wei@sina.com.

陶冠宏，博士，高級工程師。主要從事信號處理、深度學(xué)習(xí)研究。Email：taoguanhong1986@163.com.

10.11975/j.issn.1002-6819.2020.21.021

P237

A

1002-6819(2020)-21-0175-11