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

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

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

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

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

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

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

0 引 言

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

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

1 材料與方法

1.1 研究區概況

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

圖1 研究區Landsat遙感影像

1.2 研究方法

1.2.1 數據與處理

1）耕地數據

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

2）遙感數據

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

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

圖2 OLI數據提取2018年研究區冬小麥分布

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

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

1.3.1 像元相似度定義

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

1.3.2 像元匹配模型

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

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

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

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

針對第2種情況，計算當前像元時間序列曲線向量F=[f1,f2,f3, … ,f]（=1,2,3,…,N；N為上一年非冬小麥像元總數）與標準冬小麥曲線=[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為分類結果中第類像元里正確分類的像元數；x為分類結果中第類像元的總數，x為真實情況中第類像元總數。

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