遙感植被指數(shù)和CASA模型估算山東省冬小麥單產(chǎn)

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

1. 青島大學(xué)自動(dòng)化學(xué)院，山東 青島 266071 2. 青島大學(xué)計(jì)算機(jī)科學(xué)技術(shù)學(xué)院遙感信息與數(shù)字地球研究中心，山東 青島 266071 3. 中國科學(xué)院大學(xué)地球與行星科學(xué)學(xué)院，北京 100049 4. 中國科學(xué)院空天信息創(chuàng)新研究院，北京 100094

引 言

冬小麥?zhǔn)鞘澜缛蠹Z食作物之一，準(zhǔn)確模擬冬小麥單產(chǎn)及其空間分布對(duì)保證國家糧食安全和挖掘區(qū)域可利用的農(nóng)業(yè)資源具有重要意義[1]。山東省是以農(nóng)業(yè)生產(chǎn)為主的省份之一，冬小麥?zhǔn)巧綎|省主要的糧食作物之一，研究山東省冬小麥單產(chǎn)的時(shí)空分布特征對(duì)明確區(qū)域農(nóng)業(yè)生產(chǎn)現(xiàn)狀及其發(fā)展變化十分重要。

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

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

1 實(shí)驗(yàn)部分

1.1 數(shù)據(jù)與預(yù)處理

1.1.1 遙感植被指數(shù)數(shù)據(jù)

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

1.1.2 土地利用數(shù)據(jù)

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

工作中還使用了中國科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心發(fā)布的2005年、2010年和2015年的土地利用數(shù)據(jù)，提取其中的旱地分布，應(yīng)用于研究區(qū)內(nèi)冬小麥種植面積提取。該數(shù)據(jù)空間分辨率為1 km。

1.1.3 氣象數(shù)據(jù)

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

將使用到的遙感植被指數(shù)數(shù)據(jù)、土地利用數(shù)據(jù)和氣象數(shù)據(jù)統(tǒng)一重采樣為1 km，以估算山東省冬小麥單產(chǎn)。

1.1.4 農(nóng)業(yè)氣象站點(diǎn)數(shù)據(jù)

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

1.1.5 統(tǒng)計(jì)年鑒數(shù)據(jù)

所使用的統(tǒng)計(jì)年鑒數(shù)據(jù)包括研究區(qū)內(nèi)各市2000年—2016年的冬小麥統(tǒng)計(jì)單產(chǎn)和冬小麥播種面積數(shù)據(jù)，數(shù)據(jù)來源于2001年—2017年山東省統(tǒng)計(jì)年鑒。

1.2 研究方法

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

1.2.1 冬小麥面積提取方法

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

(1)

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

1.2.2 CASA模型簡(jiǎn)介

首先使用光能利用率模型CASA模擬山東省植被NPP，模型的總體設(shè)計(jì)如圖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月份的實(shí)際光能利用率，單位為g·C·mJ-1。RSG(x,t)為太陽總輻射，單位為MJ·m-2·mon-1； FPAR(x,t)表示植被對(duì)入射光合有效輻射的吸收比例，無量綱； 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)的計(jì)算可見參考文獻(xiàn)[3]。

圖1 技術(shù)路線圖Fig.1 Flowchart

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

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

1.2.3 冬小麥單產(chǎn)估算、驗(yàn)證與分析

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

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

(5)

(6)

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

使用提取的冬小麥面積分布掩膜模擬結(jié)果，獲取Case1和Case2的冬小麥模擬單產(chǎn)，采用年鑒統(tǒng)計(jì)單產(chǎn)數(shù)據(jù)分別對(duì)兩種模擬結(jié)果進(jìn)行驗(yàn)證，對(duì)比驗(yàn)證結(jié)果，選擇精度較高的模擬結(jié)果作為研究區(qū)冬小麥單產(chǎn)，分析冬小麥單產(chǎn)的時(shí)空特征。

2 結(jié)果與討論

2.1 冬小麥面積提取結(jié)果

2.1.1 MODIS EVI光譜數(shù)據(jù)S-G濾波結(jié)果

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

圖3 MODIS EVI光譜數(shù)據(jù)濾波前后曲線，以中國科學(xué)院禹城綜合試驗(yàn)站為例

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

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

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

2.2 基于固定ξmax與變化ξmax的單產(chǎn)模擬結(jié)果驗(yàn)證

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

圖5 基于(a)固定ξmax和(b)變化ξmax的模擬單產(chǎn)與年鑒統(tǒng)計(jì)單產(chǎn)的對(duì)比

2.3 山東省冬小麥單產(chǎn)模擬結(jié)果

2.3.1 冬小麥單產(chǎn)時(shí)間變化特征

從圖6中可以看出，采用模型模擬得到的冬小麥單產(chǎn)(a)與冬小麥年鑒統(tǒng)計(jì)單產(chǎn)(b)在2000年至2016年間都表現(xiàn)出明顯的上升趨勢(shì)。不同的是，模擬結(jié)果表現(xiàn)出來的上升速率(149.79 kg·hm-2·a-1)高于年鑒統(tǒng)計(jì)數(shù)據(jù)表現(xiàn)出來的單產(chǎn)上升速率(93.12 kg·hm-2·a-1)。模型模擬結(jié)果中，冬小麥單產(chǎn)在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。冬小麥年鑒統(tǒng)計(jì)單產(chǎn)在2002年最低，為4 553.55 kg·hm-2，在2015年最高，為6 175 kg·hm-2。 2.3.2 冬小麥單產(chǎn)空間變化特征 2.00年—2016年山東省冬小麥單產(chǎn)模擬結(jié)果如圖7所示?？梢钥闯觯髂攴蓍g的單產(chǎn)水平存在差異，其中2000年、2002年、2004年與2006年冬小麥單產(chǎn)比較低，與圖6(a)的結(jié)果一致。山東省在2000年春季、2002年、2006年發(fā)生干旱，導(dǎo)致這幾年冬小麥單產(chǎn)偏低。從空間上看，冬小麥單產(chǎn)在整體上呈現(xiàn)出西部高于東部的特征。而在發(fā)生明顯干旱的年份(2000年、2002年和2006年)，冬小麥單產(chǎn)較高的區(qū)域主要集中在德州和聊城的部分地區(qū)； 在不干旱或無明顯干旱發(fā)生的其他年份，山東省冬小麥單產(chǎn)較高的區(qū)域主要分布在山東省東部的德州、聊城、菏澤、濱州、濟(jì)南、濟(jì)寧和棗莊的大部分地區(qū)。

圖6 (a)模型模擬與(b)年鑒統(tǒng)計(jì)的 冬小麥單產(chǎn)的時(shí)間變化 Yield_Sta： 年鑒統(tǒng)計(jì)單產(chǎn)； Yield_Est： 模型模擬單產(chǎn)

2.4 站點(diǎn)尺度驗(yàn)證冬小麥單產(chǎn)模擬結(jié)果

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

圖7 2000年—2016年冬小麥單產(chǎn)及多年平均單產(chǎn)空間分布 其中，2012年冬小麥種植面積使用2011年的提取結(jié)果， 2016年冬小麥種植面積使用2015年的提取結(jié)果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 對(duì)最大光能利用率的討論

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

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

3 結(jié) 論

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

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

(2)在市級(jí)統(tǒng)計(jì)尺度上，使用固定ξmax得到的模擬單產(chǎn)與年鑒統(tǒng)計(jì)單產(chǎn)之間的R2為0.23，使用變化ξmax得到的模擬單產(chǎn)與年鑒統(tǒng)計(jì)單產(chǎn)之間的R2提高到了0.32。使用變化最大光能利用率模擬冬小麥單產(chǎn)優(yōu)于固定最大光能利用率的模擬結(jié)果。

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