張錦水，趙光政，洪友堂，孫智虎，段雅鳴

基于像元物候曲線匹配的生長季內河北省冬小麥空間分布識別

張錦水1,3,4，趙光政2，洪友堂2，孫智虎2，段雅鳴1,3,4

（1. 北京師范大學遙感科學國家重點實驗室，北京 100875；2. 中國地質大學（北京）土地科學技術學院，北京 100083；3. 北京師范大學地理科學學部北京市陸表遙感數(shù)據(jù)產品工程技術研究中心，北京 100875；4. 北京師范大學地理科學學部遙感科學與工程研究院，北京 100875）

及時、準確的農作物空間分布信息是進行作物長勢監(jiān)測、災害評估與產量估計的基礎。傳統(tǒng)方法一般在作物收獲期前后進行作物的識別，時間上滯后，難以滿足農業(yè)生產的應用，時空泛化能力差，模型復用程度低。該研究以歷史知識為支撐，提出冬小麥像元匹配模型（Pixel-Matched Model，PMM）進行冬小麥空間分布提取，旨在生長季內實現(xiàn)冬小麥空間分布的快速提取。研究結果表明，PMM能充分利用作物物候特征變化，排除冬小麥種植物候空間異質性的影響，能夠在播種后2個月內實現(xiàn)冬小麥的準確提取，總體精度達到了95.49%，1分數(shù)為0.83，且不隨物候曲線的延伸而大幅提高精度。與傳統(tǒng)參考曲線模型（Reference Curve Model，RCM）相比，PMM在消除區(qū)域內冬小麥生長物候差異方面具有優(yōu)勢，可在年際間實現(xiàn)冬小麥的準確識別，具有較強的時間泛化能力，能夠實現(xiàn)冬小麥的自動化識別。

模型；物候；時間序列；時間泛化；像元相似度；先驗知識；冬小麥

0 引 言

及時獲取冬小麥空間分布不僅能夠支持農作物種植管理和產量預測，也有利于滿足農業(yè)市場、農業(yè)保險、期貨市場的需求，對于研究氣候變化和災害監(jiān)測等科學問題具有重要的應用價值與科學意義[1-3]。一般而言，遙感技術進行冬小麥識別包括單時相識別、多時相變化檢測和時間序列識別方法[4-6]。單時相識別方法是利用關鍵時期影像的光譜特征，采用機器學習的方法提取出冬小麥分布[7-8]。這種方法不可避免的會遇到“異物同譜”和“同物異譜”的問題，影響冬小麥的識別精度[9-10]。多時相變化檢測方法是利用作物不同時期在遙感影像上呈現(xiàn)出的光譜差異，進行冬小麥的提取。王建興等[11]使用兩期Landsat影像，通過NDVI閾值分割技術提取出冬小麥種植區(qū)域。Zhu等[12]基于冬小麥播種期和拔節(jié)期的光譜特征變化準確識別出冬小麥的分布。相對于多時相變化檢測方法，時間序列方法是通過冬小麥在生長期內物候特征的變化規(guī)律，構建單一特征量（如NDVI）表達作物的生長物候特征，設定閾值，提取出目標作物。楊小喚等[13]利用MODIS時序NDVI特征值提取出作物面積，總體精度可以達到95%；許青云等[14]基于多年時序NDVI數(shù)據(jù)和農作物物候信息，構建了年際間動態(tài)閾值實現(xiàn)了農作物種植模式和類型的識別?？傮w來看，由于冬小麥物候在大范圍內存在差異，單一模型或閾值設定難以滿足在大區(qū)域范圍內進行冬小麥的提取，空間泛化能力低，而且模型在時間推廣能力的探討也比較少。另外，分析大尺度范圍內作物識別方法發(fā)現(xiàn)[15-16]，主要通過搜集地面實況數(shù)據(jù)和參考數(shù)據(jù)集，不斷更新和豐富訓練樣本以提高分類模型的泛化能力，在全球范圍和國家范圍內得到很好的應用，但這種做法大大增加了人力成本，而且一般在作物成熟前的4～6周內提取出高精度的空間分布[5,17-19]，并且對于大尺度作物生長季內快速識別的研究，物候差異也易被忽略。

針對傳統(tǒng)作物分布提取方法時間上滯后，大尺度作物物候差異導致模型時空推廣困難的問題，本文提出像元匹配模型（Pixel-Matched Model，PMM）進行冬小麥的遙感識別。模型結合Landsat-8和MODIS遙感影像，對比分析歷史冬小麥分布和冬小麥標準NDVI曲線來定量刻畫冬小麥的相似度，進行生長季內大尺度冬小麥分布的快速提取，并驗證模型的跨年度應用能力。

1 材料與方法

1.1 研究區(qū)概況

研究區(qū)為河北省（36°05′N～42°40′N，113°27′E～119°50′E），地形復雜，平原、高原、山地、丘陵和湖泊交錯，面積約18.88×104km2，平原區(qū)面積為8.15×104km2。該地區(qū)屬溫帶大陸性季風氣候，四季分明，降水充沛（年降水量484.5 mm）；年日照時數(shù)2 303 h，無霜期較長（204 d），集中連片的平原宜耕區(qū)和適宜的氣候是進行農業(yè)生產得天獨厚的條件。冬小麥一般在10月上旬播種，大約7 d后出苗，11月進入分蘗期，11月下旬至次年2月下旬停長，3月上旬至5月下旬進入返青至乳熟期，于6月中旬進入成熟收獲期，其種植范圍大致分布在河北省中部、南部以及東部的平原區(qū)，唐山市和秦皇島市南部也有少量種植。圖1是Landsat-8 30 m 2016年12月—2017年4月無云RGB真彩色中值合成影像，平原區(qū)的綠色光譜范圍表明了冬小麥的分布情況。在耕地范圍內由北到南選擇4個直徑40 km的冬小麥密集種植樣本區(qū)用于構建測試集和驗證集，相鄰2個區(qū)域相距約160 km，保證選擇樣本的空間跨度（圖1）。

圖1 研究區(qū)Landsat遙感影像

1.2 研究方法

1.2.1 數(shù)據(jù)與處理

1）耕地數(shù)據(jù)

耕地數(shù)據(jù)提供了冬小麥種植的范圍，輔助尋找用于冬小麥識別模型閾值的樣本，并可用于掩膜冬小麥提取結果。本研究采用清華大學宮鵬教授團隊開發(fā)的30 m全球尺度的土地利用分布數(shù)據(jù)（http://data.ess.tsinghua.edu.cn/），編號為110W40N和110W50N，從中提取出耕地數(shù)據(jù)的專題信息，并聚合到250 m，與本研究使用的MODIS空間分辨率保持一致。

2）遙感數(shù)據(jù)

MODIS衛(wèi)星的回訪周期為1 d，這種高頻的回訪周期保證了對冬小麥物候特征的表達能力。在美國國家航空航天局（National Aeronautics and Space Administration，NASA）官網（https://www.nasa.gov/）下載2017、2018和2019年冬小麥生育期影像，行列號分別為h26v04、h26v05、h27v04和h27v05。MOD09GQ數(shù)據(jù)為陸地二級標準數(shù)據(jù)產品（L2G），利用IDL 8.5平臺編程完成預處理并計算NDVI值，形成從冬小麥播種期（10月）到收獲期（次年6月）的NDVI日時序數(shù)據(jù)，利用Savizky-Glolay濾波和變步長滑動均值濾波重構時序數(shù)據(jù)[20-21]，消除噪聲的影響并能夠有效的保留作物生長周期內的關鍵節(jié)點信息。

相對于MODIS，30 m分辨率Landsat-8 OLI在空間尺度具有優(yōu)勢，能夠精確的提取冬小麥，其識別結果在本研究中有2個用途：其一，作為先驗知識支撐作物的提??；其二，作為參照值用于評價MODIS識別結果[22]。從美國地質勘探局（https://earthexplorer.usgs.gov/）下載2016年12月—2017年4月期間云量低于15%的60余景Landsat-8 OLI影像，覆蓋了停長期和拔節(jié)期2個物候期。將兩期影像進行疊加，采用支撐向量機（Support Vector Machine，SVM）提取2017年的小麥分布，逐縣在耕地范圍內隨機抽選100個像元，在外業(yè)數(shù)據(jù)的支持下目視解譯進行精度評價，總體精度要達到90%以上[23]。依據(jù)上述方法進行2018和2019年的小麥分布提取。圖2為OLI數(shù)據(jù)提取的2018年的小麥分布，可作為參照值用于評價2018年MODIS識別結果，也可作為先驗知識支撐2019年冬小麥的提取。

圖2 OLI數(shù)據(jù)提取2018年研究區(qū)冬小麥分布

1.3 季內提取空間分布模型設計

冬小麥特有的“雙峰”物候特性[24]是提取冬小麥分布的基礎。本研究以歷史冬小麥的物候特征為參照，設計基于像元的冬小麥相似性指標，發(fā)展冬小麥像元匹配物候特征模型，實現(xiàn)冬小麥生長季早期提取冬小麥的分布。

1.3.1 像元相似度定義

時間序列圖像上的每個像元，在一段時間內的物候特征曲線能夠表達為特征向量=[1,2,3,… ,s][20]。像元類型是否歸屬為目標作物（=[1,2,3, … ,c]）可以通過相似度定量刻畫，代表物候曲線的特征點個數(shù)。本研究采用歐式距離(,)來度量待分類像元與目標作物之間的相似度SC（式（1）），二者之間歐式距離越小，則相似度越大，像元被歸為目標作物的可能性越大。

1.3.2 像元匹配模型

傳統(tǒng)大尺度的冬小麥識別方法對物候差異性的考慮較少，本研究提出像元匹配模型（Pixel-Matched Model，PMM），以歷史冬小麥空間分布為基礎，在耕地范圍內確定“先驗”冬小麥和非冬小麥像元，以像元之間的相似性為基礎，通過對比冬小麥對應像元年際間的變化識別冬小麥，以消除物候特征區(qū)域差異的影響。圖3展示了像元匹配模型的技術流程。對于當前年份的冬小麥像元識別分為2種情況：去年為冬小麥和去年為非冬小麥。

針對第1種情況，根據(jù)式（1）計算像元年際間的物候曲線向量=[s1,s2,s3, … ,s]和=[s1,s2,s3, … ,s]（=1, 2, 3, … ,Y；為像元號，為特征個數(shù)，Y為去年冬小麥的像元個數(shù)）間的相似度，設置閾值1，像元的小于該閾值，則為冬小麥。閾值確定是該方法的重點，本研究采用標準差迭代法達到最優(yōu)的檢驗精度來確定閾值1。標準差迭代法是指，以多個樣本的相似度均值作為基礎，在其上疊加不同倍數(shù)的標準差得到多個測試閾值，在試驗中獲得最優(yōu)測試精度的測試閾值將被確定為最優(yōu)閾值[3]。在耕地范圍內隨機選擇400個兩年均為冬小麥的樣本像元，計算均值和標準差。統(tǒng)計學上通常以2倍標準差作為置信區(qū)間，所以將t取值范圍設為[?2.0,2.0]，步長設為0.1，代入式（2）可產生系列的閾值T，其次，在耕地范圍內獨立選擇去年為冬小麥，當前為冬小麥、非冬小麥的像元各400個作為驗證集，計算像元對應前后兩年間的相似度S，并與系列閾值T一起代入式（3）確定出像元當前類型1，通過式（4）確定識別精度，得到最優(yōu)精度的閾值1。

式中“0”代表非冬小麥，“1”代表冬小麥。

式中0表示驗證集像元檢測成功率，1表示冬小麥像元被正確識別的個數(shù)，2表示非冬小麥像元被正確識別的個數(shù)，表示驗證集像元的總數(shù)，本研究為800。

針對第2種情況，計算當前像元時間序列曲線向量F=[f1,f2,f3, … ,f]（=1,2,3,…,N；N為上一年非冬小麥像元總數(shù)）與標準冬小麥曲線=[1,2,3,…s]的相似度SF，設置閾值2，確定冬小麥。標準冬小麥特征曲線是在去年選出400個冬小麥像元，計算時間序列特征平均值得到。閾值2確定與1的方法相同。

注：T1和T2分別對應當前年份冬小麥像元識別的第1種和第2種情況的閾值；N和Y分別對應判別過程中的“否”和“是”。

1.4 精度評價

混淆矩陣是精度評價的一種技術方法，表示為×的矩陣形式，計算總體精度（Overall Accuracy，OA）、用戶精度（User’s Accuracy，UA）和生產者精度（Producer's Accuracy，PA）等指標，從不同的方面定量表達圖像分類的精度[25-26]。具體公式如下：

式中x為分類結果中第類像元里正確分類的像元數(shù)；x為分類結果中第類像元的總數(shù)，x為真實情況中第類像元總數(shù)。

1分數(shù)（1-score）是統(tǒng)計學中用來衡量二分模型分類精度的一種指標，同時兼顧了分類模型的精確率（Precision，用表示，見式（8））和召回率（Recall，用表示，見式（9）），可視為2種精度指標的加權平均（見式（10）），有效反映模型對目標類的判別能力。1分數(shù)的取值范圍為[0,1]，值越大說明分類精度越高，模型越穩(wěn)健。