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

張 偉，李 瑋，陶冠宏，李愛農，覃志豪，雷光斌，陳藝曦

·農業信息與電氣技術·

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

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

（1. 中電科大數據研究院有限公司，貴陽 550000；2. 提升政府治理能力大數據應用技術國家工程實驗室，貴陽 550000；3. 中國科學院水利部成都山地災害與環境研究所，成都 610041；4. 中國農業科學院農業資源與農業區劃研究所，北京 100081）

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

遙感；NDVI；時空數據融合；復種指數；STARFM；地表復雜區域

0 引 言

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

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

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

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

1 研究區概況

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

圖1 研究區地理位置與土地覆蓋類型

2 數據及預處理

2.1 數據源

1）Landsat數據

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

2）MODIS數據

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

3）數據預處理

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

2.2 驗證數據

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

圖2 隨機樣本點空間分布

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

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

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

2.3 輔助數據

其他的輔助數據包括該區域2010年土地覆被分類數據（圖1c）。該數據來源于中國科學院水利部成都山地災害與環境研究所數字山地與遙感應用中心。數據的生產基于中國環境減災衛星（HJ）影像并輔助Landsat TM影像作為數據源，采用了面向對象與決策樹分類器相結合的方式[27]，并經過多輪人工目視修正，共提取了38類土地覆被類型，經獨立樣點驗證，該數據一級類精度超過95%，二級類精度超過85%，滿足本研究對土地覆被數據精度的要求。

3 研究方法

該區域為一年兩熟種植區。冬季作物生長季通常為11月至第二年5月；夏季作物生長季通常為6－10月。在夏季，該區域云霧等較多，導致難以獲取高質量Landsat影像。此外，夏季作物的光譜與林地及灌木等難以區分。因此，本文在驗證精度和實用性時采用了不同的策略。精度驗證時僅使用了2010年冬季的數據對不同方法提取的作物種植區結果進行驗證。同時，為驗證方法的實用性，利用本文方法提取了綿陽市鹽亭縣2011年全年的復種指數空間分布情況，并利用統計年鑒及2013年在部分區域的實地調查數據，對提取結果進行了驗證。

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

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

3.1 STARFM算法

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

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

3.2 時序NDVI構建方案

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

3.3 樣本點驗證

基于隨機選取的1509個樣本點，逐個地對基于Landsat影像分類法和基于Landsat-like時序NDVI峰值法提取的結果進行對比驗證，并將驗證結果形成混淆矩陣。基于混淆矩陣，計算生產者精度、用戶精度和總體精度，以及Kappa系數等值，定量地對2種方法提取的冬季作物種植區結果進行精度驗證。

3.4 多尺度驗證

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

3.5 復種指數提取

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

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

4 結果與分析

4.1 波峰閾值

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

注：2010年冬季包含MOD13Q1數據時間軸2010年第21期－2011年第9期。

4.2 目視評價

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

4.3 樣本點驗證

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

注：空白區域為非耕地。

4.4 多尺度驗證

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

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

注：數字上標a代表某種類別像元的個數。例如，505表示基于Landsat影像分類法提取為休耕像元，驗證樣點也為休耕像元的個數。122表示基于Landsat影像分類法提取為休耕像元，但是驗證樣點為種植像元的個數。

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時序NDVI峰值法提取的2010年冬季作物種植區結果在不同空間尺度上的重合率對比

注：a列表示對比驗證參考數據，基于某方法提取的不同類別像元的重合率。如“種植-種植”和“種植-休耕”表示在參考影像上的種植區域內，基于某方法判定結果則分別為“種植”像元、“休耕”像元。余同。總體重合率表示基于某方法判定正確的所有像元(含種植與休耕)所占的比例。如在30 m空間尺度上，基于Landsat-like時序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種峰值法提取結果，基于Landsat-like時序NDVI峰值法的總體重合率在各個空間尺度上均高于基于MODIS時序NDVI峰值法（表2）。這與2種數據集之間的空間尺度不確定性有關[23]。由于四川丘陵地區的地表異常復雜，MODIS數據的空間分辨率較低，在該區域內存在嚴重的混合像元問題，難以獲取高精度的種植信息。通過數據融合算法，將MODIS時序NDVI數據的空間分辨率提高，能夠獲取更精細的作物種植分布信息（圖6），從而提高了提取精度。另外，前者對于種植像元的判定精度要低于后者；而前者對于休耕像元的判定精度卻明顯高于后者。MODIS數據較低的空間分辨率降低了像元內的空間異質性，也使得時間序列NDVI的波峰形態更趨于“平均”：對于作物種植的像元，其時序NDVI波峰將趨于“平坦”。這就導致基于其判定為“有作物種植”的結果更為可靠。而對于無作物種植的像元，其時序NDVI波峰將被“抬升”。這就導致基于其判定為“無作物種植”的結果較為不可靠。

4.5 復種指數結果

基于Landsat-like時序NDVI峰值法提取了鹽亭縣2011年全年復種指數的空間分布情況（圖7）。可以看出，鹽亭縣2011年耕地仍以一年兩熟為主，且集中分布于北部和東部。在西部和南部區域，有大量耕地為一年一熟。統計年鑒結果顯示，鹽亭縣2011年耕地的復種指數為1.69。而本文的方法得到的該縣復種指數結果為1.67。因此，在縣域尺度上，基于Landsat-like時序NDVI峰值法提取的復種指數結果具有很高的精度。

利用2013年8月中旬的野外調查樣點，對基于Landsat-like時序NDVI峰值法的提取結果進行了逐點驗證（表3）。提取結果與實地調查樣本的總體驗證精度達到了73.97%。與本文新方法提取的復種指數結果一致，實地調查結果中多數樣點也為一年兩熟制。表中數據顯示，有較多的一年兩熟樣點被判定為一年一熟。這是因為該區域受云霧干擾較為嚴重，部分耕地像元的時序NDVI數據的波峰未被“探測”出來。

圖7 四川省鹽亭縣2011年復種指數空間分布

表3 基于Landsat-like時序NDVI峰值法提取的鹽亭縣2011年復種指數精度驗證

注：上標a的值代表對某種類別像元進行驗證后為某種類別像元的個數。

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 數據融合算法

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

盡管多位學者對原始算法進行了改進[40,43]，改進后的算法能夠獲取更接近真實的遙感影像，但在實際應用中，通常不僅僅關注算法對影像的預測能力，更加關注構建的高時空分辨率數據集能否有助于提高某種參數的反演精度。因此，在對STARFM系列算法進行選擇時，應該綜合考慮算法的預測精度，算法的效率，實現的難易程度，對基期高空間分辨率影像數量的要求，以及數據集中數據的一致性等問題。例如：若需要對森林砍伐情況進行監測，使用STARFM原始算法構建數據集將會存在較大誤差，因為該算法對于突然發生的細小變化難以捕捉[17,41]；而在分類制圖的研究中，使用單基期Landsat影像構建的NDVI數據集，或許能有更好的分類精度，因為盡管其預測精度較低[21]，但相較于基于多基期Landsat影像構建的NDVI數據集而言，該數據集具有更好的一致性[21]。在本文的研究中，僅選擇了單基期Landsat TM影像構建時序Landsat-like NDVI數據。一方面由于高質量Landsat影像較少；另一方面則是為保證數據具有更好的一致性，因為峰值法是通過計算波峰的個數以獲取像元的復種情況，數據不一致則很可能會導致出現更多的“偽波峰”。

5.2 時間序列數據濾波算法

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

5.3 復種指數提取算法閾值

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

5.4 數據融合算法應用成效

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

5.5 不同區域實際應用策略

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

6 結 論

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

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

[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] 王學，李秀彬，談明洪，等. 華北平原 2001-2011年冬小麥播種面積變化遙感監測[J].農業工程學報，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] 范錦龍，吳炳方.復種指數遙感監測方法[J].遙感學報，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] 張偉. 川東丘陵區耕地復種指數遙感反演研究：以綿陽市鹽亭縣為例[D]. 北京：中國農業科學院，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] 辜智慧. 中國農作物復種指數的遙感估算方法研究：基于SPOT/VGT多時相NDVI遙感數據[D].北京：北京師范大學，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] 張偉，李愛農，雷光斌. 復種指數遙感監測研究進展[J]. 遙感技術與應用，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] 邊金虎，李愛農，宋孟強，等. MODIS植被指數時間序列 Savitzky-Golay 濾波算法重構[J]. 遙感學報，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] 張正健，李愛農，雷光斌，等. 基于多尺度分割和決策樹算法的山區遙感影像變化檢測方法：以四川攀西地區為例[J]. 生態學報，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] 王鵬新，荀蘭，李俐，等. 基于時間序列葉面積指數傅里葉變換的農作物種植區域提取[J]. 農業工程學報，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年耕地復種指數時空變化遙感分析[J]. 農業工程學報，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模型提高復雜地表下復種指數遙感提取精度[J]. 農業工程學報，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

國家自然科學基金青年基金項目（41701433）

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

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

10.11975/j.issn.1002-6819.2020.21.021

P237

A

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