李方杰，任建強(qiáng)，吳尚蓉，張寧丹，趙紅偉

（1. 中國農(nóng)業(yè)科學(xué)院農(nóng)業(yè)資源與農(nóng)業(yè)區(qū)劃研究所，北京 100081；2. 農(nóng)業(yè)農(nóng)村部農(nóng)業(yè)遙感重點(diǎn)實(shí)驗(yàn)室，北京 100081）

0 引言

遙感技術(shù)獲得的作物種植空間分布信息和傳統(tǒng)的農(nóng)作物面積統(tǒng)計(jì)數(shù)據(jù)一直是資源環(huán)境、生態(tài)、氣候變化和糧食安全等領(lǐng)域重要的研究基礎(chǔ)[1-2]。其中，目前的遙感技術(shù)大多通過全覆蓋方式進(jìn)行農(nóng)作物遙感識別和作物面積估算，從而獲得作物分布信息及其面積數(shù)量；表格形式的作物種植面積統(tǒng)計(jì)數(shù)據(jù)主要由統(tǒng)計(jì)部門采用傳統(tǒng)抽樣外推技術(shù)或遙感和空間抽樣相結(jié)合的外推方式獲得[3-4]。由于上述兩種作物面積數(shù)量獲取方法不同，加之遙感受混合像元、大氣干擾和復(fù)雜自然條件等因素影響，導(dǎo)致通過單純遙感技術(shù)獲取的農(nóng)作物種植面積數(shù)量與農(nóng)作物面積統(tǒng)計(jì)數(shù)據(jù)之間往往存在不一致的問題，這在一定程度上影響了農(nóng)作物空間分布遙感信息在資源環(huán)境、生態(tài)、氣候變化和糧食安全等地學(xué)領(lǐng)域中基于行政單元的空間分析應(yīng)用與尺度轉(zhuǎn)換。因此，為滿足資源環(huán)境、生態(tài)、氣候變化和糧食安全等領(lǐng)域?qū)εc農(nóng)作物種植面積統(tǒng)計(jì)數(shù)據(jù)高度吻合的高精度作物空間分布信息的迫切需求，既保證遙感獲取的區(qū)域作物面積總量與作物面積統(tǒng)計(jì)數(shù)據(jù)保持一致，又能獲得較高的作物遙感識別精度，區(qū)域作物種植面積總量控制的作物空間分布遙感提取和制圖技術(shù)成為亟需研究的關(guān)鍵技術(shù)之一。

目前，常用的農(nóng)作物遙感識別方法主要包括基于單一影像的農(nóng)作物遙感識別、基于時(shí)序數(shù)據(jù)的農(nóng)作物遙感識別和基于多源數(shù)據(jù)融合的農(nóng)作物遙感識別[5-6]。其中，基于單一影像的農(nóng)作物遙感識別主要借助一景影像，利用農(nóng)作物關(guān)鍵物候期影像光譜特征差異來實(shí)現(xiàn)農(nóng)作物的識別和面積提取。常用的分類方法包括支持向量機(jī)、最大似然分類、光譜角制圖、面向?qū)ο蠓诸惡蚄-means分類等[7-8]。該類方法的特點(diǎn)是簡單、效率高、可操作性強(qiáng)，但這類方法容易受天氣的影響，難以獲得農(nóng)作物最佳識別期影像數(shù)據(jù)，或是受遙感影像分辨率的影響，導(dǎo)致“同物異譜”和“同譜異物”現(xiàn)象普遍存在，從而使農(nóng)作物識別和提取精度受到一定程度影響[9]。

基于時(shí)序數(shù)據(jù)的農(nóng)作物遙感識別主要通過分析多景影像數(shù)據(jù)時(shí)相變化規(guī)律實(shí)現(xiàn)農(nóng)作物空間分布提取，根據(jù)特征選擇的不同可分為基于單一特征和基于多特征的兩類農(nóng)作物識別方法。其中，單一特征方法常用遙感參數(shù)包括歸一化差值植被指數(shù)（Normalized Difference Vegetation Index, NDVI）和增強(qiáng)型植被指數(shù)（Enhanced Vegetation Index, EVI）等，主要利用特征量隨作物物候變化呈現(xiàn)的規(guī)律或特征量間的相似性，通過構(gòu)建閾值模型實(shí)現(xiàn)農(nóng)作物的識別與提取[10-12]。該類方法操作簡單、運(yùn)算速度快，在作物種植結(jié)構(gòu)相對簡單的地區(qū)可以得到較高的作物識別精度。但該類方法目前大多側(cè)重于作物生育期內(nèi)完整時(shí)序或部分關(guān)鍵生育期內(nèi)時(shí)序數(shù)據(jù)，對于不同時(shí)段植被指數(shù)時(shí)序相似性及其相似性組合對農(nóng)作物面積提取精度影響鮮有研究。近些年來，一些學(xué)者[13-17]在農(nóng)作物種植結(jié)構(gòu)較復(fù)雜地區(qū)也開展利用多種遙感特征（如植被指數(shù)特征與光譜特征、時(shí)相特征、紋理特征、高程等）來實(shí)現(xiàn)農(nóng)作物的提取，取得了一系列研究成果，對進(jìn)一步提高作物空間分布提取結(jié)果精度發(fā)揮了重要作用。

基于多源數(shù)據(jù)融合的農(nóng)作物遙感識別主要包括不同時(shí)空分辨率多源遙感影像融合、遙感信息和非遙感獲取信息融合等兩種方法[18-19]。其中，基于遙感信息和非遙感獲取信息（如地形、土壤、耕地信息等）融合方法在遙感作物分類中對提高農(nóng)作物種植面積提取精度發(fā)揮了重要作用[20-21]。同時(shí)，國內(nèi)外學(xué)者[22-23]也利用農(nóng)作物面積統(tǒng)計(jì)數(shù)據(jù)與遙感信息融合的方法進(jìn)行農(nóng)作物空間分布信息提取研究，取得了較高的作物識別精度結(jié)果。然而，大多已有研究[22-23]只將作物面積統(tǒng)計(jì)數(shù)據(jù)作為輔助信息進(jìn)行作物分布提取，導(dǎo)致最終獲取的區(qū)域作物面積估算結(jié)果與作物面積統(tǒng)計(jì)數(shù)據(jù)間存在一定的差異，這在一定程度上影響了提取結(jié)果的進(jìn)一步應(yīng)用。因此，在確保作物遙感識別精度較高的同時(shí)，如何獲得區(qū)域作物面積總量與作物面積統(tǒng)計(jì)數(shù)據(jù)高度一致的作物分布提取結(jié)果需要進(jìn)一步開展深入研究。

綜上，為滿足資源環(huán)境、生態(tài)、氣候變化和糧食安全等領(lǐng)域?qū)εc農(nóng)作物面積統(tǒng)計(jì)數(shù)據(jù)一致的高精度農(nóng)作物空間分布遙感信息的迫切需求，本研究以中國重要糧食生產(chǎn)基地黃淮海平原內(nèi)河北省衡水市武邑縣為研究區(qū)，以冬小麥為研究作物，將覆蓋冬小麥全生育期的時(shí)序Sentinel-2 NDVI作為研究數(shù)據(jù)，在提出區(qū)域作物種植面積總量控制下的NDVI時(shí)序相似性閾值優(yōu)化冬小麥分布提取和制圖方法和精度驗(yàn)證基礎(chǔ)上，開展不同生育階段NDVI時(shí)序相似性及其相似性組合冬小麥分布提取精度對比分析研究，以期為在完備遙感信息及不完備遙感信息支持下區(qū)域作物種植面積總量控制的冬小麥分布提取和制圖方法方案優(yōu)化提供依據(jù)，也可為遙感數(shù)據(jù)和作物面積統(tǒng)計(jì)數(shù)據(jù)融合的大范圍農(nóng)作物分布遙感制圖及統(tǒng)計(jì)數(shù)據(jù)空間化提供一定技術(shù)方法參考和思路借鑒。

1 材料與方法

1.1 研究區(qū)概況

研究區(qū)位于中國北方糧食主產(chǎn)區(qū)黃淮海平原內(nèi)河北省 衡 水 市 武 邑 縣（37°37'N～38°00'N，115°45'E～116°08'E），該縣位于衡水市東北部（圖1a）。全縣總面積830.10 km2，總耕地面積59 533 hm2。全縣屬暖溫帶半干旱季風(fēng)氣候，年平均氣溫約為12.60 ℃，年平均降水量約為518.50 mm，年平均日照時(shí)間2 575.60 h，全年無霜期195 d，主要農(nóng)作物輪作種植制度為冬小麥-夏玉米一年兩熟制。研究區(qū)內(nèi)冬小麥種植時(shí)間為10月中旬，分蘗時(shí)間在11月下旬至12月上旬之間，12月中旬進(jìn)入越冬期，次年3月上中旬為返青期，拔節(jié)期為3月下旬至4月上中旬，4月下旬至5月上旬為抽穗開花期，灌漿乳熟期在5月中下旬，成熟期為6月上中旬。

本研究使用的樣點(diǎn)數(shù)據(jù)主要包括地面調(diào)查數(shù)據(jù)和Google Earth高分辨率影像目視解譯所獲取的樣點(diǎn)數(shù)據(jù)。全部樣點(diǎn)數(shù)據(jù)分布如圖1b所示，其中地面樣點(diǎn)數(shù)據(jù)共125個(gè)（冬小麥樣點(diǎn)67個(gè)，非冬小麥樣點(diǎn)58個(gè)），Google Earth目視解譯樣點(diǎn)共1 391個(gè)（冬小麥樣點(diǎn)768個(gè)，非冬小麥樣點(diǎn)623個(gè)）。

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

1.2.1 遙感數(shù)據(jù)

結(jié)合武邑縣冬小麥物候信息，根據(jù)影像質(zhì)量較好且影像時(shí)間間隔盡量一致的原則，本研究選取覆蓋冬小麥整個(gè)生長期（2017年10月中旬至2018年6月中旬）的18景10 m空間分辨率的Sentinel-2高分辨率影像，該影像數(shù)據(jù)主要從美國地質(zhì)調(diào)查局（https://earthexplorer.usgs.gov/）獲取，如表1所示。

表1 Sentinel-2影像列表 Table 1 List of Sentinel-2 images

由于下載的Sentinel-2遙感影像是L1C級，因此需對原始影像進(jìn)行輻射定標(biāo)、大氣校正以獲得地表真實(shí)反射率數(shù)據(jù)。在此基礎(chǔ)上，通過式（1）計(jì)算得到NDVI數(shù)據(jù)。

式中NIR為近紅外波段的反射率，R為紅光波段的反射率，NDVI值域范圍為[-1，1]。對計(jì)算得到的NDVI遙感數(shù)據(jù)按照時(shí)間順序進(jìn)行波段合成得到NDVI時(shí)序數(shù)據(jù)，由于受云、氣溶膠等大氣影響會使NDVI時(shí)序數(shù)據(jù)出現(xiàn)噪聲，從而降低數(shù)據(jù)自身質(zhì)量和作物識別精度[24-25]。因此，為提高遙感數(shù)據(jù)質(zhì)量，本研究利用Savitzky-Golay濾波方法對NDVI時(shí)間序列進(jìn)行平滑處理，最終獲取較高質(zhì)量的NDVI時(shí)序數(shù)據(jù)。

1.2.2 樣點(diǎn)數(shù)據(jù)的選取

在全部樣點(diǎn)中隨機(jī)選取185個(gè)冬小麥樣點(diǎn)和130個(gè)非冬小麥樣點(diǎn)作為訓(xùn)練樣本建立NDVI時(shí)間序列曲線，剩余650個(gè)冬小麥樣點(diǎn)和551個(gè)非冬小麥樣點(diǎn)作為驗(yàn)證樣點(diǎn)對最終冬小麥空間分布提取結(jié)果進(jìn)行精度驗(yàn)證。

1.2.3 其他輔助數(shù)據(jù)

本研究輔助數(shù)據(jù)主要包括武邑縣縣級行政區(qū)劃矢量數(shù)據(jù)、2018年武邑縣冬小麥種植面積目視解譯數(shù)據(jù)以及冬小麥物候信息等。其中，冬小麥種植面積目視解譯結(jié)果主要通過對研究區(qū)內(nèi)2018年4月20日5 m空間分辨率RapidEye衛(wèi)星遙感影像進(jìn)行目視解譯獲得。通過地面樣點(diǎn)驗(yàn)證可知，該冬小麥目視解譯空間分布結(jié)果的總體精度為99.07%，Kappa系數(shù)為0.98。通過對該高精度目視解譯結(jié)果進(jìn)行統(tǒng)計(jì)，武邑縣2018年冬小麥種植面積總量數(shù)據(jù)為27 578.776 7 hm2。為降低作物面積統(tǒng)計(jì)數(shù)據(jù)誤差對本研究結(jié)果精度帶來的不確定性影響，增加研究的可操作性，本研究將武邑縣目視解譯冬小麥面積結(jié)果作為冬小麥種植面積統(tǒng)計(jì)數(shù)據(jù)的替代數(shù)據(jù)，并作為區(qū)域冬小麥種植面積總量控制參考。另外，目視解譯冬小麥面積結(jié)果也將作為冬小麥面積遙感提取總量精度驗(yàn)證數(shù)據(jù)，從而消除統(tǒng)計(jì)數(shù)據(jù)可能存在的誤差對本研究方法驗(yàn)證結(jié)果產(chǎn)生不確定性影響。

1.3 研究方法

1.3.1 技術(shù)路線

本研究以Sentinel-2 NDVI為研究數(shù)據(jù)，利用交叉光譜匹配算法和全局優(yōu)化算法，結(jié)合區(qū)域冬小麥種植面積解譯數(shù)據(jù)，提出區(qū)域作物種植面積總量控制下的NDVI時(shí)序相似性閾值優(yōu)化冬小麥分布提取和制圖方法，并開展不同生育階段NDVI時(shí)序相似性及其相似性組合對冬小麥分布提取精度對比分析研究，以期為完備及不完備遙感信息支持下區(qū)域作物種植面積總量控制的冬小麥分布提取和制圖方法方案優(yōu)化提供依據(jù)。

本研究的技術(shù)路線如圖2所示，具體研究步驟如下：首先，根據(jù)冬小麥生長發(fā)育過程中NDVI的變化情況，對冬小麥全生育期內(nèi)的NDVI1-18時(shí)序（W）進(jìn)行劃分，共劃分為9種不同生育階段NDVI時(shí)序（S1、S2、F1、F2、F3、V1、V2、V3、V4）。其次，利用交叉光譜匹配（Cross Correlogram Spectral Matching, CCSM）算法獲得冬小麥全生育期和不同生育階段的NDVI時(shí)序相似性指標(biāo)均方根誤差（RMSEW、R MSES1、R MSES2、R MSEF1、R MSEF2、RMSEF3、 RMSEV1、 RMSEV2、 RMSEV3、 RMSEV4），在此基礎(chǔ)上，對冬小麥不同生育階段NDVI時(shí)序相似性進(jìn)行組合，共得到16種相似性組合結(jié)果（RMS ES1S2、RMSEF1F2、……、 RMSEV1V2V3V4）。然后，根據(jù)不同生育階段NDVI時(shí)序相似性及其相似性組合構(gòu)建冬小麥空間分布提取模型。在武邑縣冬小麥種植面積總量數(shù)據(jù)控制下，通過引入復(fù)合型混合演化算法（Shuffled Complex Evolution-University of Arizona, SCE-UA）進(jìn)行冬小麥提取模型閾值優(yōu)選，當(dāng)利用遙感提取的武邑縣冬小麥種植面積與冬小麥種植面積總量控制數(shù)據(jù)接近時(shí)，確定冬小麥提取模型的最優(yōu)閾值。在此冬小麥提取模型最優(yōu)閾值驅(qū)動下，輸出武邑縣冬小麥空間分布提取和制圖結(jié)果。最后，利用驗(yàn)證樣點(diǎn)數(shù)據(jù)對提取結(jié)果進(jìn)行精度驗(yàn)證和比較。

1.3.2 NDVI時(shí)序相似性

本研究基于NDVI時(shí)序曲線相似性進(jìn)行武邑縣冬小麥空間分布提取，所用到的相似性度量方法是由van der 等[26]提出的交叉相關(guān)光譜匹配（Cross Correlogram Spectral Matching, CCSM）算法。該方法引入了光譜相對滑動的概念，通過不斷計(jì)算2個(gè)光譜在波段重疊范圍內(nèi)的交叉相關(guān)系數(shù)來繪制交叉相關(guān)曲線，用實(shí)際交叉相關(guān)曲線與參考交叉相關(guān)曲線的匹配度來刻畫光譜曲線的相似性。CCSM算法具有良好的抗噪聲性能，并且對光譜幅值的變異不敏感，是一種較好的相似性度量方法[27-28]。

本研究首先利用冬小麥訓(xùn)練數(shù)據(jù)樣點(diǎn)，在時(shí)序遙感影像中得到冬小麥NDVI時(shí)間序列數(shù)據(jù)集，通過對冬小麥NDVI時(shí)間序列數(shù)據(jù)集中各時(shí)相的NDVI值進(jìn)行平均計(jì)算，并繪制時(shí)間序列曲線，最終得到冬小麥NDVI時(shí)序標(biāo)準(zhǔn)曲線。在此基礎(chǔ)上，根據(jù)交叉相關(guān)光譜匹配原理，計(jì)算冬小麥NDVI時(shí)序標(biāo)準(zhǔn)曲線和目標(biāo)像元NDVI時(shí)序曲線在不同匹配位置上的交叉相關(guān)系數(shù)，通過交叉相關(guān)系數(shù)來構(gòu)建交叉相關(guān)曲線（圖3）。其中，交叉相關(guān)系數(shù)（Rm）定義如式（2）所示：

式中m為目標(biāo)像元NDVI時(shí)序曲線的移動位置，m=0表示不移動，m=-1表示目標(biāo)像元NDVI時(shí)序曲線相對于NDVI時(shí)序標(biāo)準(zhǔn)曲線向左移動一個(gè)位置，m=1表示目標(biāo)像元NDVI時(shí)序曲線相對于參考曲線向右移動一個(gè)位置；n為移動后2條曲線的重疊波段數(shù)；NDVIr表示冬小麥NDVI時(shí)序標(biāo)準(zhǔn)曲線r中各NDVI值；NDVIt表示目標(biāo)像元NDVI時(shí)序曲線t中NDVI值。

不同m處的交叉相關(guān)系數(shù)Rm可以點(diǎn)繪成交叉相關(guān)曲線，即實(shí)際交叉相關(guān)曲線。參考交叉相關(guān)曲線的獲取則是通過將冬小麥NDVI時(shí)序標(biāo)準(zhǔn)曲線與其自身進(jìn)行交叉相關(guān)計(jì)算，并將結(jié)果繪成交叉相關(guān)曲線，即得到參考交叉相關(guān)曲線。利用冬小麥全生育期NDVI時(shí)序計(jì)算得到的參考交叉相關(guān)曲線與影像中某像元實(shí)際交叉相關(guān)曲線的對比如圖3所示，目標(biāo)像元NDVI時(shí)序曲線移動位置m的取值范圍為[-16，16]，交叉相關(guān)系數(shù)Rm的取值范圍為[-1，1]。根據(jù)所得到的冬小麥參考交叉相關(guān)曲線與目標(biāo)像元實(shí)際交叉相關(guān)曲線，本研究利用相似性指標(biāo)均方根誤差（Root Mean Square Error, RMSE）來描述2條曲線之間的差異[29-30]，RMSE的值越小，表示2條曲線越相似。因此，當(dāng)RMSE在一定的范圍內(nèi)時(shí)，目標(biāo)像元所代表的地物類型可判別為冬小麥。

1.3.3 冬小麥空間分布提取模型建立

本研究利用相似性指標(biāo)均方根誤差分別計(jì)算得到冬小麥全生育期W時(shí)序的相似性（RMSEW）和不同生育階段S1、S2、F1、F2、F3、V1、V2、V3、V4時(shí)序的9種相似性（RMSES1、RMSES2、RMSEF1、RMSEF2、RMSEF3、RMSEV1、RMSEV2、RMSEV3、RMSEV4）。根據(jù)RMSEW和RMSES1、RMSES2、…、RMSEV4構(gòu)建基于NDVI時(shí)序曲線相似性的冬小麥空間分布提取模型如式（3）所示：

式中P(i,j)表示影像中第i行、第j列的值，其中1表示冬小麥，0表示非冬小麥；C表示NDVI時(shí)序W、S1、S2、…、V4， RMSEC(i,j)表示時(shí)序W、S1、S2、…、V4第i行、第j列的相似性RMSEW、RMSES1、RMSES2、…、RMSEV4；T為時(shí)序相似性冬小麥空間分布提取模型中的閾值，當(dāng)任意像元的相似性小于等于該閾值時(shí)，即可判斷為冬小麥像元，否則判斷為非冬小麥。

對相同NDVI時(shí)序劃分次數(shù)下的不同生育階段NDVI時(shí)序相似性進(jìn)行組合，共得到16種相似性組合結(jié)果。對于不同生育階段NDVI時(shí)序相似性組合RMSEQd…QuQk，需要滿足每個(gè)NDVI時(shí)序相似性RMSEdQ、…、RMSEuQ、RMSEkQ均分別小于等于閾值Ti（i=1，2，…，n），因此，基于不同生育階段NDVI時(shí)序相似性組合的冬小麥空間分布提取模型如式（4）所示：

式中RMSEQd(i,j)、…、RMSEQu(i,j)、RMSEQk(i,j)分別表示時(shí)序相似性組合RMSEQd…QuQk中每個(gè)不同NDVI時(shí)序Qd、…、Qu、Qk第i行、第j列的相似性RMSEQd、…、RMSEQu、RMSEQk；Ti（i=1，2，…，n）為時(shí)序相似性組合冬小麥空間分布提取模型中的閾值。本研究中，冬小麥不同生育階段 NDVI時(shí)序相似性組合RMSEQd…QuQk包 括 RMSES1S2、RMSEF1F2、…、RMSEV1V2V3V4等16種。

1.3.4 冬小麥空間分布提取模型閾值優(yōu)化

復(fù)合型混合演化算法（Shuffled Complex Evolution-University of Arizona, SCE-UA）是由Duan等[31-32]提出的一種有效解決非線性約束問題的全局優(yōu)化算法，該算法綜合了多種算法的優(yōu)點(diǎn)，且具有較高的收斂速度和處理效率，能夠有效的找到全局最優(yōu)解。為了實(shí)現(xiàn)區(qū)域作物種植面積總量控制下基于NDVI時(shí)序相似性閾值優(yōu)化的高精度冬小麥空間分布提取與制圖，本研究利用SCE-UA算法，以區(qū)域冬小麥種植面積總量數(shù)據(jù)作為外部優(yōu)化數(shù)據(jù)參考對象，當(dāng)提取的區(qū)域冬小麥遙感種植面積與區(qū)域冬小麥種植面積總量數(shù)據(jù)間差異最小時(shí)，便得到冬小麥提取模型的最優(yōu)閾值。主要優(yōu)化過程如下：

1）確定尋優(yōu)空間，即確定閾值的上下界。本研究將冬小麥訓(xùn)練樣點(diǎn)像元均方根誤差的最大值作為上界，0作為下界，以此來確定尋優(yōu)空間。

2）確定SCE-UA算法中的主要參數(shù)。SCE-UA算法中的參數(shù)較多，絕大多數(shù)的值均采用已有研究成果的默認(rèn)值[33-34]，如a=2b+1、q=b+1、α=1、β=（2b+1），其中a為每個(gè)復(fù)合型的頂點(diǎn)數(shù)，b為待優(yōu)化參數(shù)個(gè)數(shù)，q為子復(fù)合型的頂點(diǎn)數(shù)，α為每個(gè)子復(fù)合型進(jìn)化后產(chǎn)生的連續(xù)后代的個(gè)數(shù)，β為每個(gè)復(fù)合型的進(jìn)化次數(shù)。復(fù)合型的個(gè)數(shù)p是唯一需要確定的參數(shù)，p值越大，越適宜于高階非線性問題，本研究設(shè)置p=2。

3）建立目標(biāo)函數(shù)。為達(dá)到閾值優(yōu)化的目的，本研究利用遙感提取的區(qū)域冬小麥面積與區(qū)域冬小麥總量數(shù)據(jù)間差值構(gòu)建目標(biāo)函數(shù)，如式（5）所示：

式中y為遙感提取的區(qū)域冬小麥種植面積與區(qū)域冬小麥種植面積總量控制數(shù)據(jù)之間的差值，hm2；E為基于像元統(tǒng)計(jì)獲得的冬小麥遙感提取的種植面積總量，hm2；S為冬小麥種植面積總量控制數(shù)據(jù)，hm2，即目視解譯獲得的冬小麥種植面積；k為研究區(qū)域個(gè)數(shù)（本研究設(shè)置k=1）。

4）確定迭代停止規(guī)則。當(dāng)滿足以下3個(gè)條件之一時(shí)，認(rèn)為目標(biāo)函數(shù)y達(dá)到最優(yōu)，停止優(yōu)化，輸出作物提取模型最優(yōu)閾值下冬小麥分布遙感結(jié)果。該3個(gè)條件包括：①閾值連續(xù)迭代5次后無顯著改變；②目標(biāo)函數(shù)值已無法顯著改善；③循環(huán)次數(shù)達(dá)10 000次。

1.3.5 冬小麥空間分布提取與制圖精度評價(jià)

本研究對區(qū)域作物種植面積總量控制下基于NDVI時(shí)序相似性閾值優(yōu)化的冬小麥空間分布提取與制圖結(jié)果進(jìn)行精度評價(jià)，主要從作物遙感識別精度和作物面積遙感估算精度兩個(gè)方面進(jìn)行驗(yàn)證。其中，作物遙感識別精度主要是基于驗(yàn)證樣點(diǎn)數(shù)據(jù)對作物空間分布提取結(jié)果進(jìn)行精度評價(jià)，采用的指標(biāo)包括總體精度（%）、Kappa系數(shù)、生產(chǎn)者精度（%）、用戶精度（%）等4項(xiàng)指標(biāo)[35]；作物面積遙感估算精度評價(jià)，則選用總量精度作為指標(biāo)進(jìn)行評價(jià)，該作物面積遙感估算總量精度（TA，%）的計(jì)算如式（6）所示：

2 結(jié)果與分析

2.1 全生育期NDVI時(shí)序相似性冬小麥提取結(jié)果

本研究利用冬小麥全生育期10 m空間分辨率Sentinel-2 NDVI時(shí)序數(shù)據(jù)和冬小麥訓(xùn)練樣點(diǎn)數(shù)據(jù)，結(jié)合交叉光譜匹配算法，對冬小麥參考交叉相關(guān)曲線與影像中各像元實(shí)際交叉相關(guān)曲線進(jìn)行全生育期相似性指標(biāo)RMSEW計(jì)算，其結(jié)果如圖4所示。由圖4可知，相似性指標(biāo)均方根誤差越小像元隸屬于冬小麥的程度越大。在此基礎(chǔ)上，本研究以武邑縣冬小麥種植面積總量數(shù)據(jù)作為外部優(yōu)化數(shù)據(jù)參考對象，以冬小麥訓(xùn)練樣點(diǎn)的平均相似性為初始閾值，利用SCE-UA算法，得到冬小麥提取模型的最優(yōu)閾值為0.227 7。根據(jù)最優(yōu)閾值對2018年武邑縣冬小麥空間分布進(jìn)行提?。▓D5），由圖5可知，2018年武邑縣整個(gè)縣域內(nèi)冬小麥分布總體較為均勻，東部和南部較為集中，西北部冬小麥空間分布較為破碎。

通過統(tǒng)計(jì)，本研究遙感提取武邑縣冬小麥種植面積為27 578.780 0 hm2，縣域內(nèi)冬小麥目視解譯種植面積為27 578.776 7 hm2，兩者相差僅為0.003 3 hm2，冬小麥空間分布遙感結(jié)果總量精度達(dá)99.99%。根據(jù)驗(yàn)證樣點(diǎn)數(shù)據(jù)構(gòu)建混淆矩陣，得到全生育期NDVI時(shí)序相似性冬小麥提取結(jié)果的總體精度為98.08%，Kappa系數(shù)為0.96，冬小麥的生產(chǎn)者精度為96.92%，用戶精度為99.53%，漏分誤差為3.08%，錯(cuò)分誤差為0.47%。這說明在冬小麥全生育期完整NDVI時(shí)序數(shù)據(jù)支持下，本研究方法具有一定的有效性和可行性，在保證高水平總量精度外，還能獲得較高水平的冬小麥識別精度，可以實(shí)現(xiàn)完備NDVI時(shí)序數(shù)據(jù)支持下本研究方法的目標(biāo)。

2.2 不同生育階段NDVI時(shí)序相似性及其相似性組合的冬小麥提取結(jié)果

2.2.1 冬小麥不同生育階段NDVI時(shí)序相似性計(jì)算結(jié)果

對9種冬小麥不同生育階段NDVI時(shí)序相似性指標(biāo)均方根誤差進(jìn)行計(jì)算，結(jié)果如圖6所示。其中，根據(jù)NDVI時(shí)序V2計(jì)算出來的相似性指標(biāo)均方根誤差RMSEV2整體偏小，像元之間均方根誤差的值相差不大，這在一定程度上對冬小麥識別較為不利，影響冬小麥遙感提取精度的提高。其他8個(gè)不同生育階段時(shí)序相似性均能在一定程度上體現(xiàn)冬小麥與其他地物之間相似性的差異，可以較好地區(qū)分冬小麥與其他地物，從而實(shí)現(xiàn)冬小麥種植面積的高精度遙感提取。

2.2.2 冬小麥空間分布提取模型閾值優(yōu)化結(jié)果

本研究利用SCE-UA全局優(yōu)化算法對不同生育階段NDVI時(shí)序相似性及其相似性組合冬小麥空間分布提取模型中的閾值進(jìn)行優(yōu)化，最終得到25種不同冬小麥空間分布提取模型最優(yōu)閾值，結(jié)果如表2和表3所示。對冬小麥不同生育階段NDVI時(shí)序相似性來說，如NDVI時(shí)序相似性RMSES1，若影像某一像元的RMSES1≤0.216時(shí)，則將該像元劃分為冬小麥；對于冬小麥不同生育階段NDVI時(shí)序相似性組合來說，如NDVI時(shí)序相似性組合RMSES1S2，當(dāng)影像中某一像元的NDVI時(shí)序相似性RMSES1S2和RMSES2分別小于等于0.363和0.494的條件同時(shí)成立時(shí)，將該像元劃分為冬小麥，否則為非冬小麥；NDVI時(shí)序相似性組合RMSEF1F2F3，當(dāng)影像中某一像元的NDVI時(shí)序相似性RMSEF1、RMSEF2和RMSEF3分別小于等于0.308、0.424和0.430的條件同時(shí)成立時(shí)，將該像元劃分為冬小麥，否則為非冬小麥；NDVI時(shí)序相似性組合RMSEV1V2V3V4，當(dāng)影像中某一像元的NDVI時(shí)序相似性RMSEV1、RMSEV2、RMSEV3和RMSEV4分別小于等于0.656、0.554、0.348和0.450的條件同時(shí)成立時(shí)，將該像元劃分為冬小麥，否則為非冬小麥。

表2 冬小麥不同生育階段NDVI時(shí)序相似性提取模型最優(yōu)閾值 Table 2 Optimal threshold values of extraction model based on similarity of NDVI time series at different growth stages of winter wheat

表3 冬小麥不同生育階段NDVI時(shí)序相似性組合提取模型最優(yōu)閾值 Table 3 Optimal threshold values of extraction model based on similarity combinations of NDVI time series at different growth stages of winter wheat

2.2.3 冬小麥空間分布提取結(jié)果精度驗(yàn)證

本研究利用9種不同生育階段NDVI時(shí)序相似性和16種不同生育階段NDVI時(shí)序相似性組合的冬小麥空間分布提取模型最優(yōu)閾值完成研究區(qū)內(nèi)冬小麥空間分布提取，結(jié)果如圖7和圖8所示。

通過與冬小麥種植面積總量控制數(shù)據(jù)進(jìn)行比較可知，25種冬小麥空間分布提取結(jié)果總量精度均達(dá)99.99%以上。在此基礎(chǔ)上，本研究利用驗(yàn)證樣點(diǎn)數(shù)據(jù)對提取結(jié)果進(jìn)行精度驗(yàn)證，如表4所示。

表4 不同生育階段NDVI時(shí)序相似性及其相似性組合的冬小麥提取結(jié)果精度 Table 4 Accuracy of winter wheat extraction results based on similarity and similarity combinations of NDVI time series at different growth stages

13種NDVI時(shí)序相似性及其相似性組合（如RMSES1、RMSES1S2、RMSEF1F2、RMSEF1F3、RMSEF1F2F3、RMSEV1V3、RMSEV1V4、RMSEV2V3、RMSEV1V2V3、RMSEV1V2V4、RMSEV1V3V4、RMSEV2V3V4和RMSEV1V2V3V4）的冬小麥空間分布提取結(jié)果總體精度和Kappa系數(shù)分別在90%以上和0.80以上；10種NDVI時(shí)序相似性及其相似性組合（RMSES2、RMSEF1、RMSEF2、RMSEV1、RMSEV2、RMSEV3、RMSEF1F3、RMSEV1V2、RMSEV2V4、RMSEV3V4）的冬小麥空間分布提取結(jié)果總體精度在80%～90%之間，Kappa系數(shù)在0.60～0.80范圍內(nèi)；2種NDVI時(shí)序相似性（RMSEF3和RMSEV4）的冬小麥空間分布提取結(jié)果總體精度在75%～80%之間，Kappa系數(shù)在0.55～0.60范圍內(nèi)。

以上結(jié)果可知，本研究方法在全生育期完備NDVI時(shí)序數(shù)據(jù)支持下能夠獲得高精度冬小麥分布提取結(jié)果外，除拔節(jié)期-成熟期NDVI時(shí)序（F3）和抽穗期-成熟期NDVI時(shí)序（V4），其余不同生育階段NDVI時(shí)序相似性及其相似性組合的冬小麥分布提取結(jié)果也能使得總量精度和識別精度均達(dá)到較高水平，進(jìn)一步證明了在不完備時(shí)序遙感數(shù)據(jù)支持下，本研究方法也具有一定的有效性和可行性，能夠滿足區(qū)域作物種植面積總量控制下的冬小麥分布提取與制圖方法精度要求。

通過比較，全生育NDVI時(shí)序相似性冬小麥空間分布提取結(jié)果精度最高，不同生育階段NDVI時(shí)序進(jìn)行相似性冬小麥空間分布提取時(shí)，冬小麥提取結(jié)果的總體精度降低了7%～20%，Kappa系數(shù)降低了0.14～0.40。

此外，當(dāng)全生育期NDVI時(shí)序被分為多個(gè)不同生育階段NDVI時(shí)序時(shí)，利用冬小麥返青前（出苗期-越冬期）的NDVI時(shí)序（S1）、越冬前期之前（出苗期-越冬期）的NDVI時(shí)序（F1）、越冬前（出苗期-越冬期）NDVI時(shí)序（V1）、越冬前期至抽穗前（越冬期-拔節(jié)期）的NDVI時(shí)序（F2）、越冬之后（越冬期-成熟期）的NDVI時(shí)序（S2）、越冬期-抽穗期的NDVI時(shí)序（V3）、越冬期的NDVI時(shí)序（V2）、抽穗期-成熟期的NDVI時(shí)序（V4）和拔節(jié)期至成熟期的NDVI時(shí)序（F3）分別進(jìn)行冬小麥空間分布提取與制圖，其冬小麥遙感識別精度依次呈現(xiàn)下降趨勢。可見依靠冬小麥出苗期-分蘗期的NDVI時(shí)序和返青期-拔節(jié)期的NDVI時(shí)序數(shù)據(jù)，可以分別獲得高精度的冬小麥分布提取和制圖結(jié)果；利用抽穗期-成熟期的NDVI時(shí)序數(shù)據(jù)進(jìn)行冬小麥空間分布提取其結(jié)果精度相對較低。這為基于不同生育階段NDVI時(shí)序相似性及其相似性組合的冬小麥空間分布提取和制圖技術(shù)研究中關(guān)鍵生育階段的數(shù)據(jù)篩選提供了一定理論依據(jù)。另外，當(dāng)進(jìn)行相同劃分次數(shù)下不同生育階段NDVI時(shí)序相似性冬小麥提取結(jié)果與相似性組合冬小麥提取結(jié)果比較時(shí)可知，在時(shí)序相似性組合的情況下，綜合利用不同生育階段NDVI時(shí)序相似性可以在一定程度上提高冬小麥提取和制圖結(jié)果的精度。例如組合時(shí)序相似性RMSES1S2的提取結(jié)果精度均高于時(shí)序相似性RMSES1和RMSES2的提取結(jié)果精度。以上研究結(jié)果可為不完備NDVI時(shí)序支持下區(qū)域作物種植面積總量控制的高精度冬小麥分布提取和制圖以及技術(shù)方案優(yōu)化提供一定依據(jù)。

3 討 論

由于高空間分辨率遙感影像（如Sentinel-2、GF-1）能夠提供更加豐富和具體的地物結(jié)構(gòu)、紋理等信息，國內(nèi)外學(xué)者[36-38]對高空間分辨率遙感影像的農(nóng)作物空間分布提取與制圖進(jìn)行了大量研究，但已有研究未利用農(nóng)作物種植面積總量數(shù)據(jù)作為控制進(jìn)行農(nóng)作物空間分布提取模型閾值優(yōu)化。如劉學(xué)剛等[36]以16 m的GF-1遙感影像作為主要數(shù)據(jù)，利用決策樹分類法提取了青島市2017年冬小麥空間分布和種植面積，經(jīng)驗(yàn)證，冬小麥提取結(jié)果總體精度為94.30%，Kappa系數(shù)為0.86，總量精度為93.60%。孫姝娟等[39]以GF-1遙感影像作為主要數(shù)據(jù)，根據(jù)水稻與其他植被和非植被的EVI曲線特征差異構(gòu)建水稻種植面積提取決策樹模型，對四川省樂至縣的水稻種植面積進(jìn)行了提取，通過驗(yàn)證，樂至縣水稻提取結(jié)果的總體精度為96.52%，Kappa系數(shù)為0.93，總量精度為97.81%。上述研究僅將農(nóng)作物種植面積統(tǒng)計(jì)數(shù)據(jù)作為輔助信息進(jìn)行農(nóng)作物空間分布提取結(jié)果的驗(yàn)證與精度評價(jià)，作物面積總量統(tǒng)計(jì)數(shù)據(jù)并未參與作物分布提取與制圖過程的關(guān)鍵環(huán)節(jié)。與上述研究不同的是，本研究將武邑縣目視解譯冬小麥面積結(jié)果（冬小麥種植面積總量真值）作為冬小麥種植面積統(tǒng)計(jì)數(shù)據(jù)的替代數(shù)據(jù)，并將冬小麥種植面積目視解譯數(shù)據(jù)作為總量控制參考，當(dāng)利用遙感提取的武邑縣冬小麥種植面積與冬小麥種植面積總量控制數(shù)據(jù)間差異達(dá)到最小化時(shí)，得到冬小麥提取模型的最優(yōu)閾值，并利用該最優(yōu)閾值獲取武邑縣冬小麥空間分布提取和制圖結(jié)果。通過驗(yàn)證，本研究冬小麥空間分布提取結(jié)果總體精度為98.08%，Kappa系數(shù)為0.96，總量精度達(dá)99.99%。與上述研究相比本研究結(jié)果總體精度提高了1.56%～3.78%，Kappa系數(shù)提高了0.03～0.10，總量精度提高了2.18%～6.39%。以上結(jié)果說明本研究方法在冬小麥種植面積總量控制下進(jìn)行冬小麥空間分布提取與制圖時(shí)，不僅作物面積總量精度達(dá)到較高水平，保證了遙感提取的作物面積與總量控制數(shù)據(jù)間的統(tǒng)一，而且作物遙感識別精度也得到一定幅度的提高?？傮w看，該研究可為高精度冬小麥分布提取和制圖技術(shù)提供一定參考依據(jù)，也可為獲取與統(tǒng)計(jì)數(shù)據(jù)高度吻合的大范圍、長時(shí)間和高精度的區(qū)域農(nóng)作物分布遙感制圖信息及統(tǒng)計(jì)數(shù)據(jù)空間化提供一定技術(shù)方法參考和思路借鑒。

此外，在大多數(shù)基于閾值模型進(jìn)行農(nóng)作物提取的同類研究中[40-41]，通常需要根據(jù)已有的先驗(yàn)知識或簡單的光譜統(tǒng)計(jì)分析，不斷的設(shè)置和調(diào)整農(nóng)作物關(guān)鍵物候期的閾值，從而獲得較高的識別和提取精度。如張佳華等[40]根據(jù)冬小麥生育期內(nèi)EVI時(shí)序曲線變化特征，建立了冬小麥空間分布提取閾值模型，并對模型中的閾值進(jìn)行了賦值，提取了華北平原典型地區(qū)冬小麥種植空間分布。盡管近年來一些學(xué)者[12,42]對合理設(shè)置閾值也進(jìn)行了一系列研究，但閾值設(shè)置的自動化程度仍有待進(jìn)一步提高。如張健康等[12]利用TM/ETM+遙感影像和MODIS EVI遙感影像，采取基于生態(tài)分類法的監(jiān)督分類與決策樹分類相結(jié)合的方法，建立決策樹模型，對黑龍港地區(qū)的主要農(nóng)作物進(jìn)行遙感提取。其中，決策樹模型中各分支約束條件的閾值通常以某一類別的EVI最小值為閾值，然后再輔以其他約束條件進(jìn)行篩選。本研究利用SCE-UA全局優(yōu)化算法實(shí)現(xiàn)冬小麥空間分布提取模型閾值自動優(yōu)選，這一應(yīng)用顯著提高了閾值設(shè)置的準(zhǔn)確性和自動化程度，減少了以往研究人員進(jìn)行作物識別時(shí)不斷設(shè)置和調(diào)整閾值所花費(fèi)的時(shí)間，并且在一定程度上能夠提高作物面積空間分布提取與制圖的效率和精度。本研究利用不同生育階段NDVI時(shí)序數(shù)據(jù)相似性進(jìn)行冬小麥空間分布提取和精度對比，對不完備時(shí)序遙感信息下高精度作物空間分布提取中遙感數(shù)據(jù)的有效應(yīng)用具有重要指導(dǎo)意義，可為高精度冬小麥分布提取和制圖技術(shù)及其方案優(yōu)化提供一定參考依據(jù)。

本研究利用全生育期完整NDVI時(shí)序相似性和不同生育階段NDVI時(shí)序相似性及其相似性組合的冬小麥空間分布提取和制圖研究取得了一定進(jìn)展，今后可以在以下幾個(gè)方面開展進(jìn)一步深入研究。1）多種特征參數(shù)結(jié)合的農(nóng)作物遙感提取是目前常用的一種方法，今后可引入其他遙感特征參量和輔助信息（如光譜特征、紋理特征、植被指數(shù)和高程等），充分利用每個(gè)特征參量所表達(dá)的信息，從而進(jìn)一步提高作物空間分布制圖和面積提取精度；2）為了實(shí)現(xiàn)本研究方法的大范圍應(yīng)用，在后續(xù)研究中可以在更大范圍尺度（如市級、省級）和更復(fù)雜種植結(jié)構(gòu)區(qū)域開展基于其他相似性指標(biāo)（如曼哈頓距離、歐式距離、光譜角匹配等）的農(nóng)作物空間分布提取和制圖精度比較研究；3）后續(xù)研究可以進(jìn)一步開展不同誤差梯度模擬總量控制數(shù)據(jù)對本研究方法的精度影響評價(jià)，從而實(shí)現(xiàn)農(nóng)作物種植面積統(tǒng)計(jì)數(shù)據(jù)控制下的作物空間分布提取和制圖。

4 結(jié) 論

本研究以中國重要糧食生產(chǎn)基地黃淮海平原內(nèi)河北省衡水市武邑縣為典型研究區(qū)，以冬小麥為研究對象，利用Sentinel-2 NDVI時(shí)間序列數(shù)據(jù)，以區(qū)域冬小麥種植面積目視解譯數(shù)據(jù)作為總量控制數(shù)據(jù)，結(jié)合交叉光譜匹配算法和SCE-UA全局優(yōu)化算法開展基于NDVI時(shí)序相似性及其相似性組合對區(qū)域作物種植面積總量控制的冬小麥制圖精度影響和比較研究，主要得到以下結(jié)論：

1）通過利用冬小麥全生育期NDVI時(shí)序相似性對本研究方法進(jìn)行驗(yàn)證可知，本研究方法獲得的區(qū)域冬小麥空間分布制圖結(jié)果總量精度達(dá)99.99%以上，總體精度達(dá)98.08%，Kappa系數(shù)為0.96，證明本研究方法不僅可以保證遙感獲取的區(qū)域冬小麥面積總量與冬小麥面積總量控制數(shù)據(jù)間的高度一致性，而且還能夠獲得較高的冬小麥遙感識別精度。可見本研究方法具有一定有效性和可行性，可為完備遙感信息支持下區(qū)域作物種植面積總量控制的冬小麥空間分布提取和制圖技術(shù)提供一定借鑒。

2）根據(jù)冬小麥不同生育階段NDVI時(shí)序相似性及其相似性組合的冬小麥空間分布提取結(jié)果與精度驗(yàn)證可知，除了抽穗期-成熟期生育階段的NDVI時(shí)序相似性冬小麥空間分布提取結(jié)果作物識別精度相對較低外，其余不同生育階段（如出苗期-越冬期、返青期-拔節(jié)期）NDVI時(shí)序相似性及其相似性組合的冬小麥空間分布提取結(jié)果也能使得總量精度和識別精度均達(dá)到較高水平。其中，出苗期-分蘗期NDVI時(shí)序（V1）和返青期-拔節(jié)期的NDVI時(shí)序（V3）數(shù)據(jù)，可獲得高精度的冬小麥空間分布提取和制圖結(jié)果；在時(shí)序相似性組合的情況下，綜合利用不同生育階段NDVI時(shí)序相似性可以在一定程度上提高冬小麥空間分布提取和制圖結(jié)果的精度。以上研究這可為不完備NDVI時(shí)序支持下區(qū)域作物種植面積總量控制的冬小麥空間分布提取和制圖技術(shù)方案優(yōu)化中關(guān)鍵生育階段數(shù)據(jù)的篩選和高精度作物分布提取和制圖提供一定理論依據(jù)。