基于像元物候曲線匹配的生長(zhǎng)季內(nèi)河北省冬小麥空間分布識(shí)別

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

基于像元物候曲線匹配的生長(zhǎng)季內(nèi)河北省冬小麥空間分布識(shí)別

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

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

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

模型；物候；時(shí)間序列；時(shí)間泛化；像元相似度；先驗(yàn)知識(shí)；冬小麥

0 引 言

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

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

1 材料與方法

1.1 研究區(qū)概況

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

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

1.2 研究方法

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

1）耕地?cái)?shù)據(jù)

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

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

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

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

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

1.3 季內(nèi)提取空間分布模型設(shè)計(jì)

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

1.3.1 像元相似度定義

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

1.3.2 像元匹配模型

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

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

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

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

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

注：T1和T2分別對(duì)應(yīng)當(dāng)前年份冬小麥像元識(shí)別的第1種和第2種情況的閾值；N和Y分別對(duì)應(yīng)判別過(guò)程中的“否”和“是”。

1.4 精度評(píng)價(jià)

混淆矩陣是精度評(píng)價(jià)的一種技術(shù)方法，表示為×的矩陣形式，計(jì)算總體精度（Overall Accuracy，OA）、用戶精度（User’s Accuracy，UA）和生產(chǎn)者精度（Producer's Accuracy，PA）等指標(biāo)，從不同的方面定量表達(dá)圖像分類的精度[25-26]。具體公式如下：

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

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