基于長時間序列遙感數據反演NPP的耕地質量評價

馬佳妮 張 超 呂雅慧 高璐璐 鄖文聚 朱德海

(1.中國農業大學信息與電氣工程學院， 北京 100083； 2.中國農業大學土地科學與技術學院， 北京 100083； 3.自然資源部農用地質量與監控重點實驗室， 北京 100035)

0 引言

耕地質量管護對耕地數量、質量、生態三位一體保護以及耕地占補平衡具有重要意義，而掌握耕地質量空間分布是耕地質量管護的重要前提。目前，耕地質量監測與評價技術方法仍以區劃—布點—調查—評價為主[1-2]，但在樣點布設科學性、數據與結果時效性、空間連續性等方面尚存在局限。遙感數據具有覆蓋面積大、獲取數據快、時空連續的特點，可有效彌補傳統方式的不足，為耕地質量監測與評價提供有效手段。

基于遙感的耕地質量監測與評價研究主要包括耕地質量指標遙感監測和耕地質量遙感反演。在耕地質量指標監測中，黃健熙等[3]利用多年蒸散發產品結合年度有效灌溉數據評價區域灌溉能力；呂雅慧等[4]采用0.5 m的GEO Eye-1 影像，構建歸一化差值植被指數(Normalized difference vegetation index, NDVI)與二維熵的決策樹，實現了耕地質量重要內容農田林網的自動識別；楊建鋒等[5]通過建立多光譜遙感影像波段與有機質含量、地形坡度、表層土壤質地的反演模型，獲得了耕地質量自然等。以上通過遙感數據信息提取與波段信息反演解決了部分指標獲取困難的問題。此外，還有通過遙感數據直接監測和評價耕地質量相關研究，林晨等[6]建立了基于MODIS數據的耕地質量自然等反演模型。以上耕地質量監測與評價多以一個時間斷點，忽略了時間斷點數據不穩定性與評價結果的時效性。

已有學者將多期遙感數據應用于土地監測與評價[7-11]。研究表明，凈初級生產力(NPP)較NDVI更能全面反映作物長勢，已廣泛應用于農田生產力評價[12-13]。以上研究采用時序數據彌補了時間斷點數據的不穩定性，但利用植被指數表征全域各類土地質量狀況，忽略了不同作物監測指標的差異，造成結果的不可比性。耕地質量各構成要素的特點和相互間的影響，決定了耕地質量的外在表現[14]。耕地質量包含自然質量、利用質量與經濟質量等[15]。而作物長勢是耕地質量利用的外在表現，產量與產能是反映耕地質量的重要指標。依據吉林省氣象站點物候統計，吉林省玉米與水稻生長期為5—9月。農田NPP在作物生長期變化明顯[16]，耕地質量影響因素如田間施肥、灌溉等人為因素和土壤理化性質在作物生長期變化明顯。多年作物生長期的NPP體現了多年影響耕地質量的各因素對作物的貢獻程度，間接評價了耕地質量。因此，本文選擇吉林省為研究區，采用2000—2010年5—9月500 m的MODIS MOD09A1遙感數據，結合氣象站點數據，構建省域耕地質量利用等評價方法，考慮不同作物植被指數差異，分別計算反映不同作物長勢的指標，經距平分析法消除量綱，利用反映多年作物長勢的NPP評價耕地質量。

1 研究區和數據

1.1 研究區概況

選擇吉林省為研究區，地理位置為東經121°38′～131°19′，北緯40°50′～46°19′，如圖1所示。吉林省屬于溫帶大陸性季風氣候，四季分明，雨熱同季。多年平均氣溫2～6℃，降水量400～900 mm。以中部大黑山為界分為東部山區與西部平原。全省耕地共553.5萬hm2，占吉林省面積28.98%。70%耕地分布在中、西部平原地區。全省耕地主要為旱地和水田。耕作制度為一年一熟制，種植結構穩定，主要農作物有玉米和水稻等，糧食單產地區間差異大，糧食平均單產在全國處于中等偏下水平。耕地質量利用等別分布在8～13等別。耕地后備資源相對充足，有31萬hm2耕地后備資源。

圖1 研究區位示意圖Fig.1 Geographical distribution map of study area

1.2 數據源

1.2.1遙感數據

采用2000—2010年玉米與水稻生長期5—9月的8 d合成500 m分辨率MODIS地表反射率產品MOD09A1(www.gscloud.cn)，共215景。經投影轉換、鑲嵌、重采樣預處理后，形成2000—2010年長時間序列遙感數據集。

1.2.2氣象數據

氣象數據來源于中國氣象科學數據共享網，包括吉林省內與周邊共51個國家氣象站點，包括2000—2010年間逐日最大溫度、最小溫度、總輻射量、日照時數等數據。

1.2.3光合有效輻射數據等

光合有效輻射數據(Photosynthetically active radiation, PAR)是利用晴空指數、太陽高度角與日照時數建立的估算模型而得[17]。其他數據包括2013年吉林省農用地等別成果中的國家利用等別數據以及吉林省行政邊界數據。

2 研究方法

NPP是植物光合作用所固定的有機物總量扣除植物自養呼吸后的剩余部分。單年NPP值易受氣候等因素影響，用于評價長期穩定的耕地質量誤差較大。計算2000—2010年的NPP均值評價耕地質量可減少氣候等因素引起的不確定性。基于影像計算每8 d的NPP值，經累加求得年NPP值，通過計算多年NPP均值，評價耕地質量。由于NPP受作物類型的影響[18-20]，需對研究區主要作物水稻(C3)與玉米(C4)進行分類，利用不同作物生長期對指數響應程度不同，構建基于時序指數的分類特征，并采用決策樹主要區分水稻與玉米，剩余少量大豆與水稻均為C3作物。不同作物類型(主要分C3和C4)分別利用植被光合作用模型(Vegetation photosynthesis model，VPM)計算NPP，并經累加計算2000—2010年NPP年累積量。后經距平分析法去除作物類型間差異，結果用于評價耕地質量利用等別。相應的技術路線如圖2所示。

圖2 技術路線圖Fig.2 Flow chart of method

2.1 基于時序遙感數據的水稻與玉米識別

由于水稻、玉米的光合作用過程、光能利用率、積溫等差異，導致兩種作物類型計算NPP不同，所以首先利用多時相進行農作物分類。通過構建增強植被指數(Enhanced vegetation index, EVI)、地表水指數(Land surface water index, LSWI)與NDVI時序數據集，選取水稻移栽期、移栽期后40 d及收獲后為關鍵時期，采用決策樹分類區分水稻和玉米。

2.1.1遙感指數計算

利用MODIS數據計算NDVI、EVI和LSWI，計算公式為

(1)

(2)

(3)

式中RB——藍波段波長(459～479 nm)

RR——紅波段波長(620～670 nm)

RNIR——近紅外波段波長(841～876 nm)

RSWIR——短波紅外波段波長(1 628～1 652 nm)

盡管8 d合成MODIS數據最大限度削弱低云、云陰影、氣溶膠等對圖像質量的影響，然而6、7月存在長時間云覆蓋情況，故利用S-G濾波消除云的影響[21]。

2.1.2基于時序指數的水稻玉米分類

依據水稻與玉米物候期時序變化差異(圖3)，采用關鍵生育期LSWI和EVI進行作物分類。選擇水稻移栽期、移栽期后40 d及收獲后作為分類關鍵期。移栽期水稻田因灌水存有2～15 cm積水[22]， LSWI時序曲線出現突增，EVI出現突降。采用EVI-LSWI≤0.05作為特征參量。移栽40 d后，水稻進入返青期，分蘗數逐漸增加，EVI迅速增大，該時段利用EVI>0.5EVImax作為分類特征。在收獲后期，利用NDVI<0.6，采用決策樹實現作物分類。

圖3 水稻與玉米時序圖Fig.3 Rice and corn timing chart

2.2 NPP計算

VPM是基于遙感數據進行NPP估算的光能利用率模型[23]。首先計算

GPP=εgFPARchlPAR

(4)

其中

FPARchl=aEVI

(5)

εg=ε0WscalarTscalarPscalar

(6)

(7)

式中εg——光能利用率，g/MJ，受溫度、葉齡與水分的影響

FPARchl——植被光合的部分吸收光合有效輻射比例

a——系數，取1

ε0——植被最大光能利用率，g/MJ，C4植物最大光能利用率高于C3植物[24]，玉米與水稻的ε0分別為3.65 g/MJ和1.15 g/MJ[25-26]

Tscalar——溫度系數

Wscalar——水分系數

Pscalar——葉齡系數

tlswimax——葉面完全伸展時期

Tmin、Tmax、Topt——作物的最低溫、最高溫、最適溫[27]，玉米、水稻的最低溫分別為0、-3℃，最高溫為45、42℃，最適溫為23、16℃

2.3 耕地質量利用等別評價

計算水稻與玉米年NPP累積量后，通過距平分析法消除作物類型間的量綱，即

(8)

式中NPP——11年NPP均值

D——像元距平百分率

最終根據D值參照農用地等別的分級標準得到耕地質量的評價結果。

3 實驗結果與分析

3.1 多時相遙感作物分類

由于水稻和玉米是吉林省主要作物，多年種植結構穩定，故采用2006年分類結果作為11年后續研究的基礎數據。分類樣本來源于實地采樣結合現有耕地圖斑，分類結果如圖4所示，總體分類精度為80.56%，Kappa系數為0.7。吉林省玉米面積為70.75萬hm2，水稻面積32.82萬hm2。玉米種植區集中在北部白城市、中部松原市。水稻種植區位于中東部通榆市、德惠市與中西部四平市，處于松花江與東遼河兩側。

圖4 2006年吉林省水稻玉米分布圖Fig.4 Rice and maize distribution map of Jilin Province in 2006

3.2 VPM模型計算NPP

利用VPM模型計算8 d NPP。后經累加得到年NPP值。并生成2000—2010年吉林省NPP空間分布圖，如圖5所示，全域NPP累積量在[0，507.249]之間。

3.3 耕地質量評價

利用距平分析法將歸一化后的D值作為表征耕地質量。并與2013年農用地分等中吉林省耕地質量利用等對比驗證。為與吉林省國家利用等別等級一致，按照吉林省國家利用等別各等別面積占比，將耕地質量反演結果分為9個等別，分別對應吉林省國家利用等的6～14等別，6為最高等別，按序等別降低。從圖6b可見，全域耕地質量呈階梯狀分布。北部白城市與松原市耕地質量處于11～14等別。中部長春市、四平市耕地質量處于9～10等別，少量分布8、12等別。南部地區在9～14等別都有分布，主要集中在10～13等別。北部與南部質量低，中部質量高。由圖6可得，反演結果與國家利用等別在整體空間分布較為一致，北部反演結果集中在9～12等別，說明南部整體反演等級偏高，主要是由于南部耕地零碎分散并存在于耕地與林地混合的像元導致反演結果偏高。對比發現，兩種結果在吉林省北部均處于9～12等別，但2013年吉林省農用地分等國家利用等別的結果位于北部的乾安縣整體處13～14等別，與周邊縣等級差異較大。圖6a中北部的評價結果空間連續，符合耕地質量實際分布。圖6b中縣域存在異常值的原因是現有耕地質量評價是以縣域為評價單元，各縣在指標獲取與計算的差異造成縣域結果整體偏高或偏低。本文方法可檢驗縣域評價結果異常。

圖5 2000—2010年吉林省累積NPP空間分布Fig.5 Cumulative NPP spatial distribution during 2000—2010 in Jilin Province

圖6 本文評價結果與2013年吉林省農用地分等國家利用等別對比Fig.6 Comparison of evaluation results with utilization of agricultural land in Jilin Province in 2013

采用反演結果等級與農用地分等成果等別之差，評價構建方法的可靠性。空間分布如圖7所示，級差大于2的區域主要集中在南部的敦化市、汪清縣、柳河縣、樺甸市，中部的公主嶺市及北部鎮賚縣。說明本文的評價結果高于農用地分等結果，主要因為耕地零散破碎，林地NPP高于耕地NPP，耕地與林地混合像元造成反演結果較高。級差小于-2的區域在中部的蛟河市與北部的洮南市。

如圖8所示，級差在-1～1級有73.72%耕地，級差-2～2級的耕地面積有17.01%，級差在3級及以上的耕地占9.27%，說明本文提出方法的評價結果與農用地分等中的分級結果具有較好的吻合度。

圖7 本文評價結果與2013年吉林省農用地分等國家利用等別級差分布Fig.7 Inversion results and spatial distribution of difference in utilization of agricultural land in Jilin Province in 2013

圖8 級差概率分布Fig.8 Level difference probability distribution

4 結論

(1)構建了利用長時間影像數據集耕地質量評價方法，通過監測多年反映作物長勢的NPP表征耕地質量。在耕地集中連片的平原地區評價結果與現有數據結果一致，同時表明，500 m的遙感數據可用于平原耕地集中連片區域的耕地質量監測與評價。

(2)利用多時序植被指數作物分類，選擇水稻

移栽期、移栽期后40 d及收獲后作為關鍵生育期，移栽期未被稻苗遮擋的田間積水，在影像上可作為區別水稻田與旱地的特征。移栽期后40 d，水稻處于急速生長的分蘗期。收獲后田間使用LSWI、EVI、NDVI指數特征參與分類，最終的分類精度達80.56%。

(3)采用2000—2010年MOD09A1影像數據，結合氣象站點數據，經VPM模型計算水稻和玉米NPP值，并得到多年NPP累積量。使用距平分析法消除作物類型差異，以此表征耕地質量，并對其結果的合理性進行驗證，結果表明，該方法適合評價吉林省中部與北部等平原地區，在耕地質量評價方面具有可行性。

(4)基于長時間序列遙感數據反演NPP的耕地質量評價，利用多年作物長勢間接反映耕地質量，避免單年評價的偶然性，提高了準確性。同時利用覆蓋全域的MODIS影像，有效解決了省域耕地質量不可比、評價由點及面與評價周期長等問題。在獲取耕地質量空間分布的同時，可以獲得耕地質量時空演變過程。