一種利用形態相似性的水稻信息提取算法

高永康，王臘紅，陳家贏，李建民

(華中農業大學 資源與環境學院，武漢 430070)

0 引言

中國65%人口的主食為稻米，2013年中國稻米產量已占世界稻米產量的28%。及時、有效地獲取水稻分布信息，對優化糧食生產布局、指導糧食生產，保障國家糧食安全具有重要意義[1]。

針對水稻種植分布的調查，傳統方法如農業報表和實地抽樣法，工作量大，成本高，時間跨度大，不能及時有效地獲取水稻種植信息以指導糧食生產。而由于遙感技術快速發展，通過遙感手段來快速、方便地獲取地面目標區域大范圍覆蓋數據，實現地表覆被的識別已成為大范圍作物監測的重要手段[2-5]。其中，歸一化植被指數(normalized vegetation index，NDVI)作為一種植被監測領域的重要指標，在作物識別、農作物產量估測、植被覆蓋動態監測等方面[6-10]均有著廣泛的應用。

利用作物的NDVI時序數據來提取作物空間分布，最直接的方式即判斷未知像元點NDVI時序曲線與目標作物標準NDVI時序曲線間的相似性，選擇恰當的閾值來判斷該未知點是否屬于目標作物類別。時間序列由于自身具有形態多變性，容易出現振幅平移、伸縮、線性漂移及時間軸伸縮和彎曲等形變[11]，使得在保證效率及準確性的基礎上對時間序列的相似性度量尤為困難。特別地，對于通過遙感手段，基于地表覆被時序數據來識別作物類別來說，由于受到天氣影響較大及作物物候空間分布差異，更增添了時序曲線的突變性。

為了提高NDVI時序數據的質量，陸俊等[12]利用時空數據融合算法(spatiotemporal data fusion algorithms，STDFA)，獲取了兼具Landsat-8 OLI影像空間分辨率和MODIS 8 d地表反射率合成產品(MOD09Q1)時間分辨率的融合影像?；谠摂祿?，使用支持向量機(support vector machine，SVM)方法得到的水稻制圖精度達到94%，效果較好，但算法實現較為復雜，融合以小窗形式逐步進行，計算資源消耗較大。郭昱杉等[13]基于歐式距離方法對冬小麥、玉米和棉花3種作物進行空間信息提取，總體制圖精度達到86.90%，但由于歐氏距離對時間序列的突變和噪聲較為敏感[14]，適用性較差，如玉米的制圖精度僅有79%。一般地，對歐氏距離的使用應基于相應的轉化，以降低其對突變和噪聲的敏感性。如動態時間彎曲距離(dynamic time warping distance，DTW)，其通過基于累積距離矩陣的動態規劃方法，求解序列之間的動態時間彎曲距離[15-16]，已有研究表明動態時間彎曲距離方法在提取地表覆被信息方面具有較好的適用性。韓曉勇等[17]使用DTW方法對秦巴山區水田、旱地等植被信息進行提取，總體制圖精度為83.80%。管續棟等[18]基于DTW方法對泰國東北部進行水稻的提取，其制圖精度為83.38%。但DTW算法復雜度較高，在等長序列比較中，相對于SVM、歐式距離等方法并沒有明顯的優勢，使其的使用受到限制[17-19]。

為找尋一種既容易實現，又具有較高的地物提取精度的方法，本研究基于歐式距離，考慮曲線的形態分布特征，提出了一種根據對應元素特征距離大小，對特征間距離進行差異性放大的算法，以此對水稻信息進行提取，并驗證算法的有效性。

1 研究區域

仙桃市位于湖北省中部，地處江漢平原，在112°55′E～13°49′E、30°04′N～30°32′N之間。市境內地形大部為沖擊平原，地勢平坦；氣候屬于亞熱帶季風性氣候，常年氣候溫和，雨量充沛，日照充足，年平均日照時數在2 000 h以上，無霜期一般256 d，是傳統的農業區。仙桃市總面積2 538 km2，2017年全市耕地面積900 km2，占比35.5%，其中水田580 km2，旱地323 km2，糧食作物以小麥、水稻為主，經濟作物以棉花、油菜為主。

2 數據與方法

2.1 MODIS數據及其預處理

本研究所用遙感數據是從NASA官網(https：//reverb.echo.nasa.gov/reverb/)獲取的2005—2015年仙桃地區MOD13Q1數據。圖幅行列號為h27v05，一年計23幅圖像，11年共計253幅。該數據已做過輻射校正、幾何校正等基本處理。MOD13Q1數據集包含NDVI圖層、EVI圖層、紅光波段、近紅外波段反射率等數據，它們的時間分辨率均為16 d，空間分辨率為250 m。

作為一種高時間分辨率的植被監測數據，MODIS NDVI數據有助于實現植被季節間及年際變化比較，研究表明，NDVI對植被的生長變化較為敏感[20-22]。MODIS NDVI數據本身已經過處理，極大降低了噪聲的存在，但時序噪聲仍然不可避免，因此考慮到噪聲的影響，本文對MODIS NDVI數據進行經驗模態分解處理，去除小尺度噪聲，重新構造時間序列，其分解過程如圖1所示。圖1中，IMF1～IMF3分別為第1～3次模態分解分量,res表示模態分解最后所得的趨勢項。根據經驗模態分解原理，IMF(intrinsic mode function)分量越靠前，其體現的變化越細微，時間尺度越??；但由于尺度越小，觀測信號受到噪聲影響越大，因此IMF分量信噪比越低。據對樣本數據點的抽樣觀測，發現基本分解到2～3次即停止分解，而且第一個IMF分量多成鋸齒狀，信噪比過低，噪聲影響過大，難以發現其變化規律，而其余分量則趨勢性比較明顯(圖1)。因此，本文決定舍棄第一個IMF分量，根據剩余分量和趨勢項重建NDVI時序數據。

圖1 經驗模態分解示意圖

2.2 統計數據及樣本獲取

小麥、油菜、中稻等作物物候信息(表1)來自中國農業部種植業管理司分省農時數據庫(http：//202.127.42.157/moazzys/nongshi.aspx)。江漢平原耕作制度基本為一年兩季，主要以2種種植模式為主：其一為小麥/油菜-中稻/棉花模式(中稻種植面積遠大于棉花，以下簡稱單季稻模式)；其二為早稻-晚稻模式(以下簡稱雙季稻模式)。

由于植被生長周期的存在，植被對紅光波段的吸收利用同樣具有周期性，而近紅外波段比較穩定，光合作用越強，紅光吸收率越高，反射率越低，NDVI值也就越高。不同植被物候期有差異，因此光合作用最強的時期也就不同，表現在NDVI時序曲線上，即是波峰所處時期有所不同。

另一方面，對于非植被來說，近紅外波段與紅光波段吸收率季節性變化不大，波動較小，故非植被反射率在時間上較為穩定，且與光照強度成正相關，表現在NDVI時序曲線上，波動較植被明顯更為平緩，NDVI值隨季節更迭做小幅波動。這些特性使得基于NDVI時序曲線對地表覆被進行信息的提取成為可能。

由表1 可以看到，單季稻模式中，在第一季作物中，小麥的抽穗期出現在3月下旬至4月中旬，油菜的開花期處于3月上旬至4月中旬，此時生長狀態最為旺盛，NDVI值會達到最大；在第二季作物中，中稻在6月上旬至7月下旬之間處于抽穗期，棉花6月下旬至8月中旬處于開花期，在此階段，二者NDVI值會達到峰值。對于雙季稻模式，從表1中可以看到，早稻會在4月下旬開始移栽，6月中旬至6月下旬左右開始抽穗，晚稻則在9月中旬至9月下旬開始抽穗，NDVI在此時期出現波峰。

表1 一年中各類型作物生長物候期分布

注：顏色深淺代表該種作物覆蓋度高低；無色表示該時期未種植該作物。

為了對文中所提形態相似性測量算法進行有效性檢驗，本文共提取了5種地物樣本數據，其中根據上述分析，由于單季稻和雙季稻具有明顯的獨特雙波峰特征，因此利用該物候特征對單季稻和雙季稻進行樣本數據的初步提取，具體篩選條件如表2和表3所示。之后，對初步提取樣本進行相關性分析，選擇樣本彼此間相關性大于0.9的樣本作為最終樣本。

表2 單季稻NDVI時序曲線特征

表3 雙季稻NDVI時序曲線特征

同時，由于研究區面積較小，為了避免采樣出現極端情況，即樣本的分布過于偏離正常區間，表4中所有樣本的獲取均是在湖北省尺度上進行。雖然表1中所列作物物候每年將會有所不同，但由于本研究所使用NDVI數據時間分辨率為16 d，物候的年際變動幅度將不足以對樣本的提取產生較大的影響。

表4 用于算法檢驗的各地表覆被類型樣本數目 個

注：以上均為11 a總體樣本。

在天地圖的幫助下，目視選取城鎮用地、水域、森林3種地表覆被類型樣本參與算法的檢驗，其中每個樣本選取區域面積基本大于4 km2，均大于1 km2。因為存在研究區11 a的NDVI數據，故為了豐富樣本總量，以減少異常數據的產生，每個樣點位置均重復采集11次，各類型樣本數目如表4所示。為了方便算法的有效性測試，表中各類樣本均是在省級水平上采樣，并經過篩選得到的最佳樣本。

仙桃市水稻種植主要為中稻，雙季稻種植面積少且分散，使用MODIS NDVI數據提取雙季稻結果客觀上可能精度不夠理想。因此在檢驗算法對實際情況的適應性時，僅提取中稻，即單季稻空間分布信息，提取時根據表2所示條件，在仙桃市范圍內對單季稻樣本進行再次取樣，結果如表5所示。

表5 2005—2015年仙桃市單季稻樣本數目

2.3 形態相似性度量算法(morphological similarity measurement algorithm，MSMA)

本研究通過NDVI時間序列之間相似性的計算，并基于合理的閾值劃分，獲取研究區水稻種植空間分布信息。為了計算NDVI時序曲線之間的相似性，首先需要獲取水稻標準時間序列，其為對樣本數據集對應元素取平均得到，見式(1)。

Sst(i)=mean(Sset(i))，i∈(1，2，3，…，n)

(1)

式中：Sst(i)表示水稻NDVI標準時序數據的第i個元素；Sset(i)表示樣本NDVI時序數據集的第i個元素集合；n為時間序列長度；mean表示取平均函數。

考慮到空間位置導致的異常值的影響，在計算標準時序曲線之前，首先通過閾值法剔除異常值，上下限閾值的確定如式(2)～式(3)所示。

up_thd(i)=pri_75(i)+1.5×pri_75_25(i)

(2)

low_thd(i)=pri_25(i)-1.5×pri_75_25(i)

(3)

s.t.i∈(1，2，3，…，n)

式中：up_thd(i)、low_thd(i)分別表示樣本數據集中第i個元素集合閾值上限、閾值下限；pri_25(i)、pri_75(i)、pri_75_25(i)則分別為樣本數據集中第i個元素集合的25%分位數、75%分位數，及二者之差絕對值。修正后的標準時序曲線獲取函數如式(4)所示。

(4)

為確定未知像元時間序列在各個時間點相對于標準時間序列的變化幅度，計算出樣本數據集各個特征元素與標準時間序列的對應元素平均絕對值距離Dmean(i)作為參考，稱之為特征均距，如式(5)所示。

(5)

式中：abs表示對元素取絕對值。

對于給定未知像元NDVI時間序列，S1=(k1，k2，k3…，kn)，及標準像元NDVI時間序列S2=(l1，l2，l3…，ln)，考慮到曲線的形態特征分布對相似性的影響，本算法不對整個時間序列計算歐式距離，而是針對2個時間序列對應元素之間計算歐氏距離Deuclid(i)，見式(6)。

(6)

于是，可得未知像元與標準像元的各個特征元素距離與相應平均距離之差，稱之為特征偏移距離，偏移距離最小為-Dmean(i)，此時認為未知像元第i個特征與標準像元第i個特征完全相同，如式(7)所示。

Drel(i)=Deuclid(i)-Dmean(i)

(7)

最終形態相似性距離計算形式如(8)式所示，其中α為偏移距放大項。當α為1時，形態相似性距離等同于絕對值距離，因此為了保持算法的穩定性，α取值應大于等于1。

(8)

為了盡量減少外界干擾誤差，偏移距放大項α應根據數據本身生成?？紤]到偏移距本身大小即是特征相似性的重要度量，簡單地以線性方式對偏移距進行放大，單位偏移距放大效果相同，缺乏對偏移距本身大小的重要性的突出考慮，因此對于α項，應考慮非線性方式放大。本文取指數函數作為放大基本式，如式(9)所示。

(9)

式中：對各特征元素特征偏移距與特征均距做除法，結果即為偏移距相較于特征均距的倍數，即偏移倍數，取偏移倍數作為基礎指數。對基礎指數，即偏移倍數加1以確保α項取值范圍為大于等于1的實數。由于α項的確定僅依賴于數據本身，因此無人為及其他因素干擾。

綜上所述，本研究中形態相似性度量距離計算如式(10)所示。

(10)

其中，為了增強算法的適應性，防止過度地對偏移距進行放大，算法中添加β為偏移倍數調整因子。β越大，偏移放大效應越明顯；當β為0時，該算法等同于絕對值距離。由于噪聲的存在，適宜的β值應是使算法在有著恰當的噪聲放大效果的同時，可以最大化地物時序曲線類別的差異。

2.4 閾值選擇與精度控制

針對使用相似性度量來提取地物類別的方法，閾值選取對精度至關重要。關于閾值的選取，存在標準差法、分位數法、試湊法等。對于標準差法，一般選擇標準樣本曲線加減2倍或3倍樣本標準差(可以不為整數倍數)所形成的上下界曲線與標準樣本曲線的距離，作為判斷待分類像元是否歸屬目標類別的標準?？梢酝ㄟ^改變標準差倍數來改善結果精度，但如果需要達到滿意的精度，該方法所使用的上下界曲線范圍內包含的樣本數目并不確定，可以很小，此時上下界曲線并沒有代表性；可以很大，甚至所包含樣本數等于樣本總數，但此時上下界曲線同樣有過于遠離樣本的可能，相應的所計算出來的距離閾值不再可靠。分位數法是使用可以包含大部分樣本的距離作為閾值，相比標準差法，分位數法可以確切知道距離閾值可以囊括的樣本比例，所獲得的閾值較之更具有代表性，但它的缺點仍舊是所獲閾值一般不能囊括所有樣本，這點與標準差法較為相似。

試湊法常與其他閾值選取方法結合使用，如標準差法。分位數法也可以看作試湊法。單純的試湊法所獲閾值無法判斷其可靠性，盡管可以得到極高的用戶或制圖精度，但結果的準確性難以判定。

為了解決上述問題，基于本研究提出的MSMA方法，結合分位數法和試湊法，在以樣本距離值域上限作為閾值的基礎上，做出一些改進：選擇樣本距離的百分之百分位數作為距離閾值。同時由于本形態相似性距離測量方法存在β偏移倍數調整因子，不同的因子值會有不同的差異放大效果，調整因子大小，閾值同時會相應變化。該方法閾值獲取方式不變，始終為該形態距離度量方法計算的樣本距離的最大值，保證囊括了所有參與計算的合格樣本，具有較強的合理性，同時在調整β參數以獲得最佳MSMA模型的過程中，閾值始終處于動態調整之中，增強了閾值的適用性。

2.5 算法實現

為了測試算法的有效性，本研究使用了單季稻、雙季稻、森林、城鎮、水域共5種不同類型的地表覆被NDVI時序數據進行算法的檢驗。其中，各自選取樣本數量的50%作為各自類型的訓練數據集，使用剩余的樣本NDVI時序數據作為驗證數據集。

模型有效性驗證過程如下。

1)模型訓練。①對于某種地物，在訓練中，隨機選取該地物訓練集的20%，在每次迭代中，使用MSMA算法獲取一次標準樣本數據、樣本集與標準數據的平均絕對值距離(計算方法如公式(5)所示)及閾值。②將標準樣本數據、平均絕對值距離作為輸入，結合剩余的所有地物的各自80%訓練集的集合，再次使用MSMA算法，獲取相似性度量值Dshape。③使用①中所得閾值，對結果進行類別判定，精度計算。④從訓練結果中篩選出使得用戶精度最大化的β值。

2)模型驗證。①將模型訓練中④所得β值作為已知條件，隨機選取相應地物驗證集的20%。使用MSMA算法獲取樣本標準數據、驗證樣本集與標準數據的平均絕對值距離(計算方法如公式(5)所示)及閾值。②將標準樣本數據、平均絕對值距離及模型訓練中所得β值作為輸入，結合剩余所有地物類別各自80%驗證集的集合，再次使用MSMA算法，獲取相似性度量值Dshape。③使用①中所得閾值，對結果進行類別判定，并計算精度。④模型可行性評估。

模型訓練過程和模型驗證過程相似，但模型驗證過程不需要進行迭代計算。

3 結果與分析

3.1 算法有效性檢驗

訓練結果如表6所示，其中β為算法偏移倍數調整因子。由于并非全部訓練樣本參與訓練結果的統計分析，因此參考樣本數為80%的訓練集樣本。

由于本步中最優β值對應下的各類地物歸屬判斷均為獨立的算法實現，故表6、表7均為非混淆矩陣。同時算法的初步有效性檢驗是在5種已經確定類別的樣本點集間進行，以算法對5種地物的分辨能力來判斷算法有效性與否，因此此處認為不存在混合像元?；旌舷裨獑栴}在本文中僅存在于對仙桃市單季稻的提取中。單季稻是仙桃市最為主要的糧食作物，通過僅對單季稻信息提取，可以最大限度地減小混合像元的影響。

表6 形態相似性度量算法訓練結果

注：表中每行表示對訓練中所有地物類別對應的各自80%訓練集的集合進行該行對應地物的樣本類別歸屬判斷；括號中數字為該列對應類別樣本錯分為該行對應類別的樣本所占百分比。

表7 形態相似性度量算法驗證結果

注：表中每行表示對驗證中所有地物類別對應的各自80%驗證集的集合進行該行對應地物的樣本類別歸屬判斷；括號中數字為該列對應類別樣本錯分為該行對應類別的樣本所占百分比。

在算法實現中，可以通過調整β值來改變檢測精度。由于算法是采用指數映射形式來增大測試曲線與標準曲線的差異，因此β的值不能過大，否則偏移映射值將超出算法檢測范圍，導致精度的降低甚至算法無效。

本研究中β值設置為0～3之間，間隔值為0.05，理論上每個地表覆被類型測試集將進行60次迭代測試(但當精度提前達到設定精度閾值，將終止迭代)，選取使用戶精度最為接近100%的β值，作為算法最終偏移倍數調整因子。

表6中可以看到，算法的偏移倍數調整因子最小為0.05，最大為0.35，此時5類地表覆被類型算法判分結果的用戶精度基本均接近100%。其中，森林參考集中完全沒有被錯分為其他4類地物的樣本，其他4類地物同樣完全沒有樣本被錯分為森林，說明MSMA算法對森林與其他4類地物之間具有完全的分辨能力。在單季稻類別歸屬判斷時，雙季稻參考集中有11個樣本被錯分為單季稻，錯分比例為1.67%，其他3類完全沒有被錯分樣本，說明MSMA算法對單季稻與森林、城鎮、水域之間具有完全的分辨能力，單季稻和雙季稻之間具有較強的分辨能力。在雙季稻類別歸屬判斷時，單季稻參考集中錯分樣本最多，但錯分比例僅為0.92%，城鎮參考集中僅有2個樣本被錯分為雙季稻，錯分比例依然較低，為0.85%，水域和森林完全無錯分樣本，說明此時MSMA算法對雙季稻與森林、水域之間具有完全的分辨能力，雙季稻和城鎮、單季稻之間具有較強的分辨能力。在城鎮類別歸屬判斷時，僅有水域參考集中存在錯分樣本，錯分比例為0.85%，在水域類別歸屬判斷時，僅有城鎮參考集中存在錯分樣本，錯分比例為0.42%，說明MSMA算法對城鎮、水域與單季稻、雙季稻、森林之間具有完全的分辨能力，城鎮與水域之間具有較強的分辨能力。對于制圖精度而言，5種地表覆被類型測試結果相對用戶精度較低，但同樣達到了較高水平，其中單季稻的制圖精度最高，達到97.76%，雙季稻制圖精度最低，為90.24%。

對錯分情況分析發現，錯分事件在5種地表覆被中基本僅發生在單季稻和雙季稻，城鎮和水域之間，盡管存在樣本數量過少和偶然情況(如測試集的選取是隨機的，每次的選取結果均不盡相同，這很有可能導致測試結果的差異)的影響，但這也同樣說明了單季稻和雙季稻之間，城鎮和水域之間的混淆概率相較其他組合較大。

對每種地物樣本NDVI時序數據進行異常值剔除，之后進行取平均操作，獲取5種地表覆被類型標準NDVI曲線如圖2所示。分析發現，單季稻和雙季稻所在耕地均為一年兩季耕作模式，一年中會有2個時期為作物生長旺盛期，此階段NDVI值最大，因而其NDVI時序曲線存在2個波峰；而由于森林的生長不同于作物的種植，其生長期較長，受人類活動影響較小，因此在一年中僅存在一個階段為生長旺盛期，且該階段較長，因此其NDVI時序曲線僅存在一個波峰，與作物差異明顯，并且森林的NDVI值長期處于高值；相較于作物與森林，城鎮和水域在時間上的光譜響應最為遲鈍，其光譜響應的敏感程度較低，基本僅與光強的大小成正相關，表現為NDVI值整體較小，NDVI時序曲線趨勢與光照強度變化一致。因此，5種地表覆被中，城鎮和水域的NDVI時序曲線最為相似。這可能為單季稻和雙季稻、城鎮和水域易混淆的原因之一。

圖2 各地表覆被類型NDVI標準時序曲線

表7為MSMA算法驗證結果。可以看到，各類地物錯分事件分布基本與訓練情況一致，僅是數目有所差別。在β值固定的條件下，單季稻和雙季稻用戶精度有所上升，分別為99.85%和96.50%，森林用戶精度保持不變，城鎮和水域用戶精度有所下降，但基本穩定；對制圖精度而言，5種地物均有所下降，其中雙季稻最低，制圖精度為88.13%，但仍然較為理想，單季稻制圖精度最高，為96.84%?？傮w而言，MSMA模型訓練和模型驗證前后用戶精度和制圖精度保持穩定，說明了MSMA算法具有一定的穩健性，在地物分類領域具有操作的可行性。

為了進一步驗證MSMA算法可靠性，本文同時檢驗了動態時間彎曲距離方法(DTW)對5種地物的分辨能力，使用的數據與MSMA算法驗證樣本相同，并使用分位數法作為閾值獲取方式，在不同分位數水平下，各類樣本錯分情況如圖3所示，用戶精度及制圖精度結果如表8所示。從表8可以看出，隨著分位數的增長，制圖精度隨之增大，當分位數為95%時，制圖精度均在90%以上；而對于用戶精度，水域和單季稻表現較好，其他3種地物則基本在50%以下；5種地物中，唯有水域的用戶精度和制圖精度均保持較高水平，森林的用戶精度最低，一直保持在10%以下，雙季稻次之，保持在20%左右。

圖3 不同分位數水平下各類地物樣本錯分比例(DTW方法)

表8 不同分位數水平下DTW方法分類精度對比 %

由圖3分析知，在從樣本集中對森林進行識別時，存在大量的單季稻和雙季稻樣本被分為森林類別，在對單季稻進行識別時，則有大量的森林和雙季稻被分為單季稻。而對于水域，各類地物則有著相對較小的錯分比例，當分位數為90%以下時，水域錯分為各類地物的比例均低于10%，各類地物被錯分為水域的比例同樣較低?？梢姡悩颖镜拇罅垮e分是各類別用戶精度低下的直接原因。

綜上所述，在使用DTW方法對5種地物NDVI時序數據進行分類時，水域與其他4類有著最大的時序差異；而除水域之外的4種地物則具有較小的時序差異。這是因為，DTW方法原理即是通過對時間序列的伸縮來最小化未知時序數據與標準時序數據差異，其注重時序空間上的時序曲線的形態，而對位置要求不高。由圖2可知，在本次實驗中，5種地物的NDVI時序曲線均呈現單峰或雙峰特征，而由于DTW度量2條曲線的相似性時，并不是嚴格的特征一對一匹配，而是伴隨著錯位匹配，導致了除非地物時序曲線在取值上有著較大差異；否則在使用DTW方法度量時，單峰和單峰曲線、單峰和雙峰曲線間及雙峰曲線間的差異性很容易被DTW方法的錯位匹配極大程度地降低，最終導致結果的不可靠。

3.2 單季稻提取及精度評價

從表7可知，當β值固定為0.05時，單季稻用戶精度和制圖精度均達到較高水平，可以最大限度地將單季稻與其他地物區分開來，因此在單季稻的提取中，使用0.05作為β值的最優值。仙桃市各年單季稻樣本數量如表5所示。單季稻提取步驟如下。

1)隨機將各年樣本分為80%和20%的樣本集1和樣本集2。

2)基于樣本集1，第一次實現MSMA算法，以得到相應年份的標準單季稻NDVI曲線、平均絕對值距離Dmean(i)，i∈(1，2，3，…，n)，及閾值。

3)將標準單季稻時序曲線、平均絕對值距離與β值作為輸入，第二次實現MSMA算法，計算各年單季稻相似性度量值Dshape。

4)使用各年20%的樣本驗證結果的可靠性。

5)制圖。

表9為按照以上步驟，得到的精度驗證結果。從表9可以看出，11年中，單季稻提取制圖精度基本達到90%以上，其中低于90%的有2年，2009年制圖精度最低，為84.21%，2008年和2014年制圖精度最高，均為97.30%，11 a總體精度達到93.29%，說明基于MSMA的單季稻提取效果較好。

表9 2005—2015年仙桃市單季稻空間提取精度

表10為仙桃市11 a間單季稻面積提取結果?？梢钥吹?，與空間提取精度相比，面積提取精度明顯偏低，11 a中，有6個年份面積精度低于80%。2015年面積提取精度最高，接近86%，2008面積提取精度最低，為74.57%。面積精度偏低的原因主要源于混合像元的存在，有3種形式：①本文假設不同地表覆被具有不同的NDVI時序曲線，通過衡量不同覆被NDVI時序曲線的相似性來區分不同的覆被類型，但由于作物物候期重疊現象的存在，如棉花與中稻幾乎具有相同的物候特征(表1)，同一種NDVI時序曲線可能表示多種地表覆被類型，導致提取的面積數據一般大于統計數據；②本文所用遙感數據時間分辨率為16 d，這使得本來不具有相同物候特征的植被類型，在遙感影像上反而可能表現出相同的物候特征，從而導致了混合像元的產生；③本文所用遙感數據空間分辨率為250 m，空間異質性低，每個像元可能包含多種地表覆被類型，是混合像元產生的最主要原因。因此，雖然表10所示面積精度提取結果較低，但符合實際情況，也從面積上驗證了MSMA算法的有效性。

表10 2005—2015年仙桃市單季稻面積提取精度

注：由于仙桃市年鑒收集不全，統計面積中帶*號的數據為基于當年稻谷總產量、中稻單產、中稻占稻谷總產量一般比例推算得到。

使用MSMA算法對仙桃市單季稻提取結果分析發現，仙桃市單季稻種植基本分布全境，其中，仙桃境內中西部、中北部、東南部在11 a的提取結果中均有較大的空白區域出現。中西部空白區域為仙桃稻蝦養殖區域，由于水稻收割之后需要灌水以形成蝦類的生存環境，一年中除了中稻生長時期，基本被水覆蓋，因此該類地表覆被類型NDVI時序曲線既不同于單季稻，也不同于水域。本研究中對單季稻樣本的提取是基于中稻與其他作物(小麥、油菜等)連作的假設(這正是一般情況)，因此雖然稻蝦養殖模式仍然可看作單季稻模式，但在提取結果中卻未能把稻蝦養殖區域判分為單季稻區域。同時在2005—2010年，該地區稻蝦養殖區域逐漸擴大，之后逐漸縮小。

11 a間，在仙桃各地區中，西北部單季稻種植變化最為明顯，2005年、2009年、2014年均有較為顯著的變化。根據天地圖判斷，東南部為仙桃市灘湖所在，靠近長江，周圍魚塘較多，并且同樣存在稻蝦養殖模式，使得單季稻判分錯誤。至于中北部空白區域，則為仙桃市區所在，在11 a中，該區域有明顯的擴展，與城市化進程相符。整體而言，11 a間仙桃市單季稻種植空間分布基本穩定。

4 結束語

基于本研究所提出算法，對仙桃市單季稻區域進行提取，總體制圖精度為93.29%，總體面積精度為80.57%，達到較高水平，證明了MSMA算法的可靠性。MSMA算法相較于其他相似性度量算法有以下優點：

1)相較于動態時間彎曲距離算法，MSMA更易于實現，復雜度較低，分類精度更高。

2)相較于歐式距離的特征差異簡單累加，MSMA算法基于特征差異動態調整思想，增強了對不同地類的區分能力。

3)陸俊等[12]基于STDFA數據融合模型和SVM算法得到的水稻制圖精度達到94%;而MSMA算法以更低的實現復雜度，達到了相似的效果。

雖然MSMA具有諸多優點，但本研究仍然存在一些問題：

1)本研究中單季稻樣本的獲取是基于小麥/油菜-中稻模式的物候期特征獲取的，雖然這種模式在研究區為主要作物種植模式，但仍然存在其他模式具有相似的物候期特征(如小麥-棉花模式)，因此提取的單季稻樣本中包含一定的其他作物模式的樣本，這使得結果的可靠性下降。

2)本研究使用的數據為MOD13Q1數據集，其中NDVI數據空間分辨率為250 m，因此空間分辨率低下導致的混合像元的存在將是影響提取結果可靠性的另一重要因素，混合像元的分解將能進一步改善分類效果。

3)MOD13Q1數據集中NDVI數據時間分辨率為16 d，時間分辨率較低，不能保證可以獲取水稻淹水期NDVI數據;而水稻淹水期NDVI特征是區分水稻與其作物的重要物候特征，這使得單純使用作物NDVI時序曲線，依靠作物時序上的特征提取作物空間分布存在局限。

在算法有效性檢驗和后續單季稻提取中，MSMA均有著較好的表現，說明了該算法的有效性。一旦上述問題得到解決，相信使用該算法獲得的目標類別提取結果會具有更高的可靠性。