一種基于時序遙感影像混合像元分解的耕地種植強度估算方法
——以湖北省為例

劉軼青， 陶建斌， 陳 曦， 陳 濂， 衛詩琪

(地理過程分析與模擬湖北省重點實驗室/華中師范大學城市與環境科學學院， 武漢 430079)

農業是人類社會生存和發展的基石，耕地是進行農業生產的基本資源和條件.隨著城市化迅速發展，大量耕地被轉化成建設用地[1-2].同時，隨著人口持續增長、飲食結構改變和非食物用途(如生物能源)農產品消耗的增加，人類對農產品的需求不斷.而在我國，市場經濟的發展對農民種糧收益的沖擊、農戶個體追求利益最大化而增加經濟作物生產，使得部分地區農作物種植頻率呈現下降趨勢(表現為耕地撂荒和冬閑田).如何緩解這種人地矛盾，提高現有耕地的集約化利用程度，是國際社會普遍關注的焦點問題[6].

現有耕地集約化利用的指標定義主要有種植頻率和復種指數[7-12].種植頻率是對種植模式的硬劃分(如單季作物和雙季作物)，缺乏詳細的空間信息[13-15].復種指數的計算多基于統計數據，忽略了行政區內作物分布狀況的空間異質性[16]，且數據存在時效性和準確性問題[17].隨著遙感技術的發展，利用時間序列數據，進行大尺度的耕地集約化利用估算成為可能[18-20].但在我國由于復雜地形和復雜種植結構導致的混合像元的影響，估算結果存在不確定性.總之，現有耕地集約化利用指標及其估算方法不能滿足精細作物制圖的要求，如何獲得高時空分辨率的數據集亟待解決.

本文基于MODIS(Moderate Resolution Imaging Spectroradiometer)時間序列數據，在時間維度構建一個特有的特征空間，基于線性光譜分解的方法進行混合像元分解來估算種植強度.本文方法為復雜地形、復雜種植結構下高時空分辨率的種植強度估算提供了一種新方法，對于耕地的集約化研究有著重要現實意義.

1 研究區域及數據

湖北省地處長江中游(108°21′42″E～116°07′50″E、29°01′53″N～33°6′47″N)，是我國主要的糧食生產基地，素有“魚米之鄉”的美譽.全省70%的耕地主要集中分布在江漢平原、鄂東沿江平原及鄂中丘陵地區；兼有水田和旱地，主要作物類型有水稻、小麥、玉米、油菜等.

圖1 研究區域土地覆蓋狀況Fig.1 The main land-cover types in the research area

本文采用的主要數據為：2015年MOD13Q1(MODIS/Terra Vegetation Indices 16-Day L3 Global 250m SIN Grid)EVI(Enhanced Vegetation Index)數據，來源于美國國家航空航天局MODIS WEB網站(https://modis.gsfc.nasa.gov/).基于湖北省矢量行政邊界，研究區所涉及圖幅編號為h27v05、h27v06、h28v05、h28v06.研究數據可滿足大尺度、長時間連續監測地表植被的需要.

在數據預處理上，采用MRT工具(MODIS Reprojection Tool)進行影像的拼接、重采樣、重投影處理.同時，考慮到數據受到云、雪、陰影等諸多因素的影響，原始EVI曲線季節變化不明顯，利用數學方法，使用TIMESAT工具重建原始數據.TIMESAT作為平滑長時間序列數據的專業濾波處理工具，集成了高斯濾波以及S-G濾波等濾波方式.本文采用的是S-G濾波算法，除達到去云、消除離異值的目的以外，還使得EVI數據在時間維度上的時序變化更加穩定.

除此之外，高分辨率遙感數據采用了Landsat 8多光譜遙感數據，空間分辨率為30 m，重返周期為16 d，數據來源于美國地質調查局(https://glovis.usgs.gov).結合農作物的物候規律，選擇時相在1月中下旬至2月上旬的遙感影像，可以較好區分二者，獲取較為純凈的雙季農作物端元.在該時間段內，雙季農作物中的第一季(主要是冬小麥和冬油菜)長勢較好，而單季農作物尚未播種.在標準假彩色合成方式下，雙季農作物在影像中呈現紅色，單季農作物或者撂荒地呈現灰黃色.本文中，Landsat 8遙感數據主要用于端元選取的參考和樣本數據的制作.相比于MODIS數據，能在較大尺度上滿足精細觀測的需要，提供更為豐富的地表信息.

2 研究方法

2.1 樣本數據制作

本文利用Landsat 8遙感影像制作與MODIS空間分辨率一致的耕地種植頻率的樣本數據(圖2).

利用ISODATA方法對樣本區Landsat 8遙感影像進行分類，最小分類數定為10，最大分類數定為20，迭代20次.通過分類后處理，將分類結果合并成單季農作物、雙季農作物和非耕地這3個類別.將樣本數據與估算得到的種植強度數據在像元尺度上進行比較，以用于精度檢驗與評價.

圖2 種植頻率樣本數據制作示意圖Fig.2 The schematic diagram of preparing cropping frequency samples

2.2 特征空間構建

根據農作物物候規律，雙季農作物的EVI曲線在一年中的第65 d和第225 d左右會各自出現一個波峰，在第145 d左右出現一個波谷.水體及建設用地在全年內EVI曲線總體數值較低，沒有較大波動.天然植被在第145 d之后長勢最好，EVI曲線出現一個明顯的峰值.對時相EVI數據分析后發現，不同土地覆蓋類型具有不同EVI時間譜曲線(圖3)，農作物具有鮮明的區別于其他地物類型的季節性物候特征.這些明顯的物候特征，是區分不同地物類型的基礎.

本文選取主成分變換后的前兩個主分量構建特征空間，特征空間呈現為三角形(圖4).利用ENVI軟件的N維可視化工具，建立從影像空間到特征空間的投影，從而定位、識別特征空間的端元.通過查看其時間譜曲線后，發現這3類點群分別為：雙季農作物、天然植被、水體.綠色一端為天然植被，從頂點到中心，土地覆蓋類型向灌叢和草地過渡；紅色一端為雙季農作物，趨近中心，向單季農作物過渡；藍色一端為水體，從頂點至中心，土地覆蓋類型過渡為建設用地.

圖3 主要土地覆蓋類型時間譜曲線Fig.3 Time series profiles of the main land-cover types

圖4 利用PCA 1和PCA 2構建二維特征空間Fig.4 The two-dimensional feature space based on PCA 1 and PCA 2

2.3 端元選擇

端元(End Member)代表某種純凈的地物，所選端元的質量直接影響混合像元分解的精度[21-22].根據所構建的特征空間，端元的選擇應覆蓋雙季農作物、天然植被、水體這3類主要土地覆蓋類型.考慮到不同地域同一地物類型的差異性，所選端元在研究區內盡量均勻分布且端元總數達到像元總數的3%[23].

本文利用不同農作物的物候規律，以關鍵物候期增強型植被指數的假彩色合成影像為底圖，在時間序列 MODIS 影像上選擇端元，雙季農作物在影像中顯示為品紅色.為保證端元的純凈性，本文還利用 Landsat8遙感影像以及Google Earth影像上大塊的連續的純凈地物作為參考，輔助進行端元選取.

2.4 豐度估算

本文利用全限制性線性光譜分解的方法進行豐度估算.不同于在光譜維度上進行混合像元分解，本文方法是在時間維度上進行的.通過對一年23期的MODIS時間序列遙感影像進行混合像元分解，得到雙季農作物、天然植被、水體的豐度影像.線性光譜分解模型表達式(式1)及約束條件(式2、式3)如下：

(1)

(2)

(3)

式中，Rib為一年中的第b期影像上的第i像元的像素值；fki為對應于第i像元的第k個基本組分所占分量值；Ckb為第k個基本組分在一年中的第b期影像的像素值；n為像元i所包含的基本組分數目；εib為誤差項.

通過解混可得到3類主要地物的豐度數據，考慮到豐度的物理意義，本文利用MATLAB軟件將豐度影像的值域控制在0～1.雙季農作物的豐度可代表耕地種植強度信息，0值代表撂荒地，1代表較為純凈的雙季農作物像元，根據豐度影像的像素值高低可直觀判斷耕地種植強度大小.

圖5 3種主要地物類型的端元選擇(影像RGB合成方案：DOY 065， 145， 065)Fig.5 Endmember collections of the main land-cover types (RGB composite: DOY 065， 145， 065)

3 結果與分析

通過豐度估算，本文得到湖北省2015年耕地種植強度結果(圖6).對關鍵物候期的MODIS假彩色合成影像、湖北省種植頻率數據、種植強度估算結果進行目視比較，三者具有較好的空間一致性(圖7).相比于MODIS原始數據，基于混合像元分解所得的種植強度數據，降低了混合像元的影響，空間分辨率得到提高.相比于湖北省種植頻率數據，種植強度數據避免將耕地種植模式硬性劃分為單季、雙季或非耕地，所蘊含的信息更真實、豐富.

圖6 湖北省2015年耕地種植強度Fig.6 Cropping intensity in Hubei Province in 2015

(a)和(b)為MODIS假彩色合成數據(影像RGB合成方案：DOY 065， 145， 065)，(c)和(d)為湖北省種植頻率數據，(e)和(f)為種植強度估算結果(a) and (b) are MODIS false color composite( DOY 065，145，065)， (c) and (d) are cropping frequency in Hubei， (e) and (f) are the estimated cropping intensity圖7 2015年種植強度估算結果目視比較Fig.7 Visual Comparison of cropping intensity in 2015

除目視比較之外，本文也采用定量方法進行精度檢驗.首先，將本文得到的種植強度估算結果與2015年中國農田熟制遙感監測數據集進行縣級尺度的比較[24]，R2達到0.90.

圖8 湖北省2015年縣級尺度的精度檢驗Fig.8 Accuracy validation at county level in Hubei Province in 2015

其次，選擇和樣本數據同一區域的影像進行像元尺度上的驗證.利用相關性分析，驗證二者間的線性關系，R2達到0.84(圖9).

圖9 湖北省2015年像元尺度的精度檢驗Fig.9 Accuracy validation at pixel level in Hubei Province in 2015

4 討論與結論

時間序列遙感影像使得耕地資源集約化利用的精確表達成為可能.本文引入“種植強度”這一概念，種植強度是亞像元尺度的農作物分布和種植頻率狀況的綜合表征.相較于“種植頻率”、“復種指數”等傳統耕地集約化利用指標，種植強度避免將種植模式硬性劃分為單季或雙季農作物，得到更高時空分辨率的子像元水平的豐度信息，更準確地表征了農作物的種植“強度”.

利用“種植強度”這一概念，在時間維度構建特征空間，并使用線性光譜分解的方法估算耕地種植強度，在亞像元尺度對耕地集約利用程度進行表達.經檢驗，估算結果精度較高，證明了本文方法提取種植強度的可靠性.

本文主要探討的是一種基于時間序列遙感數據進行混合像元分解的種植強度估算方法.不同于多數已有的在光譜維度上進行的混合像元分解，本文的混合像元分解方法是在時間維度上進行.不同的土地覆蓋類型具有其特有的EVI時間譜曲線，因此，農作物具有鮮明的不同于天然植被、水體的物候特征.本文構造的特征空間，能夠較好地區分農作物與其它土地覆蓋類型，以及不同種植強度的農作物.基于此提取的耕地種植強度能夠在亞像元級別實現對耕地集約利用程度的表達，優于傳統的種植頻率和復種指數，同時，可滿足高時間分辨率和高空間分辨率的耕地集約利用監測要求.

農作物與天然植被的季節性物候差異是本文進行混合像元分解的理論基礎.但考慮到湖北省境內農作物物候的地域性差異帶來的“同物異譜”的現象[25]，以及所收集到的端元不夠純凈的情況，精度在一定程度上受到影響.其次，精度驗證時所利用的Landsat遙感數據，受到獲取時相的限制，也無法完全反映驗證區域的真實種植情況.