遙感植被指數和CASA模型估算山東省冬小麥單產

張 莎，白 雲，劉 琦，童德明，徐振田，趙 娜，王兆雪，王霄鵬，李詠沙1, ，張佳華

1. 青島大學自動化學院，山東 青島 266071 2. 青島大學計算機科學技術學院遙感信息與數字地球研究中心，山東 青島 266071 3. 中國科學院大學地球與行星科學學院，北京 100049 4. 中國科學院空天信息創新研究院，北京 100094

引 言

冬小麥是世界三大糧食作物之一，準確模擬冬小麥單產及其空間分布對保證國家糧食安全和挖掘區域可利用的農業資源具有重要意義[1]。山東省是以農業生產為主的省份之一，冬小麥是山東省主要的糧食作物之一，研究山東省冬小麥單產的時空分布特征對明確區域農業生產現狀及其發展變化十分重要。

光能利用率模型是估算植被生產力[2-3]和作物單產[4-5]的常用模型之一，如光能利用率(the Carnegie-Ames-Stanford approach，CASA)模型[3]等。該類模型具有一定的機理性，且所需輸入數據較少，在區域尺度進行植被生產力或作物單產模擬時易于使用。但是，目前的研究對模型中的關鍵參數之一，最大光能利用率ξmax，仍然存在較大的爭議和不確定性[2, 5-7]。在使用光能利用率模型模擬作物單產時，不同學者對作物最大光能利用率的大小也存在爭議。Lobell[4]根據實驗數據，測量了1993年—1994年和1999年—2000年小麥、玉米和大豆的實際光能利用率在2.1～2.6 g·MJ-1PAR之間。任建強等[5]在研究中取黃淮海平原內石家莊市、邢臺市和衡水市冬小麥3-5月份的實際光能利用率為固定常數1.25 g·C·mJ-1PAR。由于作物實際光能利用率受作物溫度、降水等影響[6]，由此推斷，研究文獻[4-5]中冬小麥的ξmax必高于1.25 g·C·mJ-1PAR。劉真真等[8]綜合分析大量國內外文獻，將河北省邯鄲市南部三縣冬小麥ξmax取值為2.8 g·C·mJ-1。有研究在兩個光能利用率模型中分別設置玉米ξmax為4.68和6.13 g·C·mJ-1。史曉亮等[9]根據實測數據估算了松嫩平原玉米ξmax為3.08 g·mJ-1。上述研究中，通常根據植被類型或作物種類將ξmax設置為固定常數，然而已有研究表明同一植被類型的ξmax也會隨時間變化而變化，比如農作物，由于品種更替和管理措施水平的提高，使得農作物的ξmax得到了一定水平的提升[9]。采用光能利用率模型模擬作物單產的模型僅模擬一年的作物單產，可以忽略ξmax隨時間的變化。對于一年或短時間段的模擬，這種忽略是可以接受的，但如果使用光能利用率模型模擬長時間序列的作物單產，就需要考慮ξmax隨時間的變化。然而目前這一方面的研究還比較缺乏。

本研究以山東省為研究區，采用固定最大光能利用率(固定ξmax)和隨時間變化的最大光能利用率(變化ξmax)分別驅動CASA模型模擬研究區內2000年—2016年植被凈初級生產力(net primary productivity，NPP)，結合收獲指數(harvest index, HI)模擬冬小麥單產； 將兩種情況下模擬得到的冬小麥單產在市級尺度上進行均值統計，與市級尺度統計單產做相關性分析，分析最大光能利用率對冬小麥單產估算的影響，并進一步分析冬小麥單產的時序變化特征。

1 實驗部分

1.1 數據與預處理

1.1.1 遙感植被指數數據

使用的遙感植被指數數據包括時序歸一化植被指數(normalized difference vegetation index，NDVI)與增強型植被指數(enhanced vegetation index，EVI)，均使用中等分辨率成像光譜儀(moderate resolution imaging spectro-radiometer，MODIS)植被指數產品數據。MODIS NDVI數據下載自地理空間數據云(http://www.gscloud.cn/)，是MODIS影像中國區合成數據，時間分辨率為月，空間分辨率為500 m，使用2000年—2016年間數據，用于輸入CASA模型。MODIS EVI數據下載自Earthdata(https://search.earthdata.nasa.gov/search)，時間分辨率為16 d，空間分辨率為1 km，行列號為h27v05，本文使用2000年—2015年數據，用于提取研究區內冬小麥空間分布。

1.1.2 土地利用數據

研究中使用了2000年—2016年的MODIS土地利用數據，MCD12Q1產品，其時間分辨率為年，空間分辨率為500 m。該數據包括5種分類體系，本文采用第一種分類體系。

工作中還使用了中國科學院資源環境科學數據中心發布的2005年、2010年和2015年的土地利用數據，提取其中的旱地分布，應用于研究區內冬小麥種植面積提取。該數據空間分辨率為1 km。

1.1.3 氣象數據

驅動CASA模型需要的氣象數據包括： 月平均氣溫、月降水量和月太陽總輻射。為避免站點數據插值過程帶來的不確定性，使用ERA-Interim氣象再分析資料的氣溫、降水和太陽總輻射數據。該數據時間分辨率為1 d，空間分辨率為0.125°。將日尺度氣象要素數據合成為月值，作為CASA模型的輸入數據。

將使用到的遙感植被指數數據、土地利用數據和氣象數據統一重采樣為1 km，以估算山東省冬小麥單產。

1.1.4 農業氣象站點數據

使用了包括研究區及周圍省份在內的農業氣象站冬小麥播種期和成熟期數據，構建了多元線性模型模擬冬小麥的播種期和成熟期：DoYS=0.031P+0.072T+261.5，DoYM=-0.009T+1.298R+119.7，其中DoYS和DoYM分別表示擬合的播種期DoY(day of year)和成熟期DoY，P，T和R分別表示生育期內的降水(mm)、氣溫(℃)與總輻射(0.001 W·m-2)。此外，還使用了山東省內德州、聊城、菏澤和曹縣共4個農業氣象站點的多年冬小麥單產記錄數據，用于單產模擬結果的補充驗證。這些數據均從中國氣象數據網(http://www.nmic.cn/)獲取。

1.1.5 統計年鑒數據

所使用的統計年鑒數據包括研究區內各市2000年—2016年的冬小麥統計單產和冬小麥播種面積數據，數據來源于2001年—2017年山東省統計年鑒。

1.2 研究方法

技術路線如圖1所示。首先使用差分法結合光譜突變法提取了研究區內冬小麥每年的種植面積分布； 其次使用CASA模型在兩種模擬情景下模擬每年的植被NPP； 最后，結合收獲指數和提取的冬小麥面積，獲取2000年—2016冬小麥單產的時空分布并進行分析。

1.2.1 冬小麥面積提取方法

在研究時段內，山東省逐年冬小麥面積分布的提取方法使用了先前研究工作提出的提取流程，以進一步驗證該提取方法的普適性。提取過程中，使用農業氣象站觀測播種期與成熟期和氣象再分析資料獲取的站點氣溫、降水和輻射建立多元線性回歸模型，以使用氣象再分析資料模擬冬小麥播種期和成熟期的空間分布； 采用Savitzky-Golay(S-G)濾波重構時序EVI序列，根據模擬的播種期和成熟期的空間分布，截取每個像元在播種期至成熟期之間的重構EVI序列，采用差分法計算截取后得到的EVI序列的峰值頻數； 獲取每個像元成熟期前后的EVI，采用光譜突變法計算其在冬小麥成熟后相較于成熟前的EVI突變情況(Slope)； 提取波峰個數為2和光譜下降斜率低于負0.02的像元，并與旱地分布取交集，獲取冬小麥的空間分布。其中，光譜突變法計算Slope的公式如式(1)

(1)

式(1)中，EVIm-16和EVIm+16分別為冬小麥成熟期前、后16 d的EVI值。

1.2.2 CASA模型簡介

首先使用光能利用率模型CASA模擬山東省植被NPP，模型的總體設計如圖2所示。CASA模型最早由Potter[2]提出，模型表述如式(2)—式(4)

NPP(x,t)=APAR(x,t)×ξ(x,t)

(2)

APAR(x,t)=RSG(x,t)×FPAR(x,t)×0.5

(3)

ξ(x,t)=Tξ1(x,t)×Tξ2(x,t)×Wξ(x,t)×ξmax

(4)

其中，APAR(x,t)為像元x在t月份吸收的光合有效輻射，單位為MJ·m-2·mon-1；ξ(x,t)為像元x在t月份的實際光能利用率，單位為g·C·mJ-1。RSG(x,t)為太陽總輻射，單位為MJ·m-2·mon-1； FPAR(x,t)表示植被對入射光合有效輻射的吸收比例，無量綱； 0.5表示植被所能利用的光合有效輻射占太陽總輻射的比例。Tξ1(x,t)和Tξ2(x,t)表示溫度脅迫因子；Wξ(x,t)為水分脅迫因子；ξmax為植被最大光能利用率。FPAR(x,t)，Tξ1(x,t)，Tξ2(x,t)和Wξ(x,t)的計算可見參考文獻[3]。

圖1 技術路線圖Fig.1 Flowchart

圖2 CASA模型框架圖Fig.2 Frame of CASA model

為探究最大光能利用率對作物單產模擬的影響，在驅動CASA模型時，設置了Case1和Case2兩種模擬情景，其中Case1設置ξmax為不隨時間變化的固定值2.8 g·C·mJ-1，Case2設置ξmax為隨時間變化的序列值。參照前人研究工作，在Case2中設定2014年冬小麥ξmax為2.8 g·C·mJ-1，按ξmax每年增加0.03 g·C·mJ-1·a-1的速率[9]推算2000年—2016年間其他年份的冬小麥最大光能利用率。

1.2.3 冬小麥單產估算、驗證與分析

光能利用率模型模擬作物單產的原理如式(5)和式(6)

Yield(x)=NPPsum(x)×Tc×HI

(5)

(6)

其中，Yield(x)為像元x的模擬單產(kg·hm-2)； NPPsum為冬小麥生長季內NPP的和，根據山東省冬小麥物候特征，取3月—5月NPP的和[5]；Tc為植物體內碳素轉換為干物質的轉換系數，取值為2.22[5]； HI為收獲指數，取其值為0.48，與已有研究[5]中的收獲指數取值接近； NPP(x,t)為像元x在t月份的NPP模擬值，單位為g·C·m-2·mon-1。

使用提取的冬小麥面積分布掩膜模擬結果，獲取Case1和Case2的冬小麥模擬單產，采用年鑒統計單產數據分別對兩種模擬結果進行驗證，對比驗證結果，選擇精度較高的模擬結果作為研究區冬小麥單產，分析冬小麥單產的時空特征。

2 結果與討論

2.1 冬小麥面積提取結果

2.1.1 MODIS EVI光譜數據S-G濾波結果

在提取冬小麥種植面積時，使用了光譜突變法和差分法，其中的差分法對小峰極其敏感，在使用差分法之前，需使用S-G對原始時序MODIS EVI數據進行濾波，以去除小峰波動對峰值頻數統計的影響。圖3展示了2014年—2015年冬小麥生長季內中國科學院禹城綜合試驗站所在像元的EVI經S-G濾波前后的曲線?？梢钥闯?，濾波前的光譜數據在冬小麥越冬期間存在一些小峰，濾波后的曲線平滑掉了這些小峰，更接近真實狀況。從濾波后光譜曲線中可以看出，冬小麥在一個生長季內有2個峰，分別出現在第305天和次年第113天，大約對應研究區內冬小麥的分蘗期和開花期，與實際情況相符。還可以看出，在第129—161天時，冬小麥

圖3 MODIS EVI光譜數據濾波前后曲線，以中國科學院禹城綜合試驗站為例

EVI明顯下降，是由于第161天前后大約為冬小麥成熟期，在冬小麥成熟收割后，EVI主要表現為裸地的光譜特征，與成熟前EVI相比發生明顯下降。生育期內峰值個數和成熟收割導致的光譜突降，這兩個光譜特征不受年份變化的影響，也不因區域不同而存在差異，具有較好的普適性，本文也正是利用以上這兩個光譜特征進行了冬小麥種植面積提取。 2.1.2 冬小麥面積提取結果及驗證

將所提取的2000年—2015年冬小麥按市進行匯總，使用2000年—2015年年鑒統計的市級冬小麥面積對其進行驗證，驗證結果如圖4(a)所示。由于數據質量問題，導致提取的2012年冬小麥面積出現較大誤差。剔除掉2012年17個市的提取結果與其他年份3個市的缺失值，共剩余252個數據對。在后續分析中，2012年使用了2011年的提取結果。從圖4(a)可以看出，本文提取的山東省2000年—2015年冬小麥種植面積與相應的年鑒統計面積之間吻合較好，兩者之間的決定系數(R2)達0.71，平均相對誤差(RMSE)為111.03 k·hm2，提取的冬小麥種植面積大部分在90%的置信區間內。為簡化文章的圖表內容，本文展示2013年僅一年的提取結果[圖4(b)]。

圖4 (a)2000年—2015年冬小麥提取結果驗證與(b)2013年冬小麥種植面積分布

2.2 基于固定ξmax與變化ξmax的單產模擬結果驗證

使用2000年—2016年市級統計單產分別驗證采用固定ξmax和變化ξmax模擬得到的冬小麥單產(圖5)?？梢钥闯?，在市級尺度上，固定ξmax模擬得到的單產與年鑒統計單產之間的決定系數為0.23，而變化ξmax模擬得到的單產與年鑒統計單產之間的決定系數提高到了0.32，高于已有研究結果[10-11]。采用變化ξmax的模擬結果優于固定ξmax的模擬結果，這說明在使用光能利用率模型模擬作物單產時，對于同一種作物采用固定ξmax會帶來較大誤差，應考慮ξmax在年際間的變化。以下均基于變化ξmax的模擬結果進行分析。

圖5 基于(a)固定ξmax和(b)變化ξmax的模擬單產與年鑒統計單產的對比

2.3 山東省冬小麥單產模擬結果

2.3.1 冬小麥單產時間變化特征

從圖6中可以看出，采用模型模擬得到的冬小麥單產(a)與冬小麥年鑒統計單產(b)在2000年至2016年間都表現出明顯的上升趨勢。不同的是，模擬結果表現出來的上升速率(149.79 kg·hm-2·a-1)高于年鑒統計數據表現出來的單產上升速率(93.12 kg·hm-2·a-1)。模型模擬結果中，冬小麥單產在2000年、2002年和2004年最低，分別為4 932.61，4 994.39和4 943.38 kg·hm-2，2016年最高，為8 168.58 kg·hm-2，其次為2014年，為7 631.78 kg·hm-2。冬小麥年鑒統計單產在2002年最低，為4 553.55 kg·hm-2，在2015年最高，為6 175 kg·hm-2。 2.3.2 冬小麥單產空間變化特征 2.00年—2016年山東省冬小麥單產模擬結果如圖7所示?？梢钥闯?，各年份間的單產水平存在差異，其中2000年、2002年、2004年與2006年冬小麥單產比較低，與圖6(a)的結果一致。山東省在2000年春季、2002年、2006年發生干旱，導致這幾年冬小麥單產偏低。從空間上看，冬小麥單產在整體上呈現出西部高于東部的特征。而在發生明顯干旱的年份(2000年、2002年和2006年)，冬小麥單產較高的區域主要集中在德州和聊城的部分地區； 在不干旱或無明顯干旱發生的其他年份，山東省冬小麥單產較高的區域主要分布在山東省東部的德州、聊城、菏澤、濱州、濟南、濟寧和棗莊的大部分地區。

圖6 (a)模型模擬與(b)年鑒統計的 冬小麥單產的時間變化 Yield_Sta： 年鑒統計單產； Yield_Est： 模型模擬單產

2.4 站點尺度驗證冬小麥單產模擬結果

在驗證兩種模擬情景的冬小麥單產模擬結果時，使用了所提取的冬小麥種植面積掩膜模擬結果之后，又在地級市區域尺度取平均值，是對冬小麥單產模擬結果的區域尺度驗證。為進一步驗證模擬結果，使用山東省德州、聊城、菏澤和曹縣共4個農業氣象站的冬小麥單產記錄對兩種模擬情景的模擬結果分別進行了驗證。由于農業氣象站多數座落在城鎮或城鎮與農田交界處，導致很多站點在所提取的1 km空間分辨率的冬小麥種植面積空間分布中處于非冬小麥種植區，因此僅使用了研究區內這4個處于冬小麥種植區的農業氣象站點數據。為了降低混合像元的影響，在提取這4個站點的單產遙感模擬結果時，使用了窗口半徑為7的圓形緩沖區，將緩沖區范圍內像元的均值作為相應站點的單產模擬值。刪除掉農業氣象站記錄缺失或遙感提取顯示無冬小麥種植的數據，最終得到德州站點2000年—2002年、2005年—2006年、2006年、2008年—2010年和2012年，聊城站點2001年、2005年—2010年和2012年，菏澤站點2001年、2005年、2009年和2011年，曹縣站點2001年—2010年和2013年，共32站—年的數據。這32站—年的農業氣象站單產記錄和兩種模擬情景(固定ξmax和變化ξmax)得到的單產之間的R2分別為0.17和0.29，均分別低于2.2節得到的兩種模擬情景下的R2(0.23和0.32)，這是由于混合像元造成的站點尺度與像元尺度不匹配，但依然呈現出變化ξmax模擬情景得到的R2高于固定ξmax模擬情景得到的R2，與2.2節的結果是一致的，進一步說明了采用變化ξmax能有效提高CASA模型模擬作物單產的精度。

圖7 2000年—2016年冬小麥單產及多年平均單產空間分布 其中，2012年冬小麥種植面積使用2011年的提取結果， 2016年冬小麥種植面積使用2015年的提取結果Fig.7 Distribution of estimated winter wheat yield in each year (a—q) and the average (r) The planted winter wheat areas in 2011 and 2015 are used to be as those in 2012 and 2016, respectively

2.5 對最大光能利用率的討論

從圖6(a)和圖7中可以看出，2016年冬小麥單產明顯偏高，這是由于假設ξmax隨時間的推移而增加，增加速率為0.03 g·C·mJ-1·a-1，由此計算得到的2016年冬小麥ξmax為2.86 g·C·mJ-1，該值可能高于當前條件下的冬小麥最大光能利用率。作物的最大光能利用率是作物的生理屬性之一，受基因的控制。隨著作物育種技術的提高，作物品種的更替，作物的ξmax也在增加。因此，作物ξmax會表現出隨時間變化而增加的特征。但是，作物的育種技術和最大光能利用率并非隨時間線性增加，因此使用線性關系描述冬小麥ξmax隨時間的變化存在不足，作物最大光能利用率的時間變化特征還需要進一步探討。

此外，本研究僅考慮了冬小麥ξmax的時間變化，對整個研究區來說，在某一年份使用的仍是一個常數值，缺乏對冬小麥ξmax空間差異性的考慮。由于最大光能利用率是作物的生理屬性，如果要考慮它在空間上的差異性，則需要獲取作物品種的空間分布，而目前這一數據較難獲取。在未來的工作中，如果能夠獲取不同區域(如地級市)冬小麥播種的主要品種，可以根據品種分區設定冬小麥的最大光能利用率，以在模擬冬小麥單產時考慮最大光能利用率的空間差異性。

3 結 論

采用MODIS EVI數據、ERA氣象再分析資料和農業氣象站觀測生育期數據，使用差分法結合光譜突變法，提取了山東省冬小麥種植面積； 使用MODIS NDVI數據和氣溫、降水與輻射數據，采用固定ξmax和變化ξmax分別驅動CASA模型，結合收獲指數與提取的冬小麥種植面積，模擬了山東省2000年—2016年冬小麥單產，主要結論如下：

(1)濾波后的MODIS EVI數據能夠較好地反映冬小麥生育期內的光譜特征。光譜突變法與差分法結合能有效提取2000年—2015年山東省冬小麥種植面積，具有較好的普適性。在研究時段內，與年鑒統計面積之間的R2達0.71，RMSE為111. 03 k·hm2。

(2)在市級統計尺度上，使用固定ξmax得到的模擬單產與年鑒統計單產之間的R2為0.23，使用變化ξmax得到的模擬單產與年鑒統計單產之間的R2提高到了0.32。使用變化最大光能利用率模擬冬小麥單產優于固定最大光能利用率的模擬結果。

(3)山東省冬小麥單產在2000年—2016年間呈明顯上升趨勢，基于統計數據的上升速率和基于模型模擬結果的上升速率分別為93.12和149.79 kg·hm-2·a-1。山東省冬小麥單產在空間上呈現西部高于東部的分布特征。