基于線性雙源遙感蒸散模型的作物估產(chǎn)
——以河套灌區(qū)為例

杜鑫， 蘇濤*

(1. 安徽理工大學(xué)空間信息與測繪工程學(xué)院， 淮南 232001; 2. 安徽理工大學(xué)礦山采動災(zāi)害空天地協(xié)同監(jiān)測與預(yù)警安徽普通高校重點實驗室， 淮南 232001; 3. 安徽理工大學(xué)礦區(qū)環(huán)境與災(zāi)害協(xié)同監(jiān)測煤炭行業(yè)工程研究中心， 淮南 232001)

河套灌區(qū)是中國七大農(nóng)業(yè)主產(chǎn)區(qū)之一，具有多種種植結(jié)構(gòu)。灌區(qū)內(nèi)降雨稀少且蒸發(fā)量大，黃河引渠灌溉是作物生長用水主要的來源之一[1]，因此估算灌區(qū)作物蒸散量和產(chǎn)量對節(jié)約用水、優(yōu)化灌溉制度具有重要意義。近年來隨著遙感技術(shù)的發(fā)展，利用遙感手段監(jiān)測不同種類的作物種植結(jié)構(gòu)并觀察大范圍作物的生長分布情況，為提高作物產(chǎn)量，增加土地利用率和提高生產(chǎn)效益提供數(shù)據(jù)支撐[2]。許多學(xué)者利用深度學(xué)習(xí)等技術(shù)建立了多種土地覆蓋分類的算法，目前常見的作物分類方法主要有遙感影像分類、單波段閾值法和多波段譜間關(guān)系法等[3]，其中利用多時相時間序列影像數(shù)據(jù)并結(jié)合歸一化植被指數(shù)(normalized difference vegetation index， NDVI)監(jiān)測作物的豐度和長勢可以有效區(qū)分作物類別，實現(xiàn)作物生育期的持續(xù)監(jiān)測[4]。大量的研究證實NDVI是植被生育期最敏感的特征參數(shù)之一，它與植被覆蓋率、光合作用指數(shù)、葉面積指數(shù)、生物量等植被特征指數(shù)密切相關(guān)[5]。NDVI 的值域為-1～1，一般情況下根據(jù)經(jīng)驗閾值，當(dāng)NDVI<0.15時，可將該區(qū)域看作沒有綠色植被或者少有綠色植被覆蓋[4]。而在作物生育期對不同時間序列NDVI閾值分類可以較為簡便的提取作物種植結(jié)構(gòu)[6-7]。灌區(qū)內(nèi)主要農(nóng)作物有玉米、葵花和小麥，其生育期長勢監(jiān)測對作物估產(chǎn)具有重要意義。目前使用較多的是利用作物水分生產(chǎn)函數(shù)進行作物監(jiān)測，因其能夠反映作物耗水與產(chǎn)量兩者之間的相對關(guān)系而被廣泛地應(yīng)用在農(nóng)田的試驗上[8]。王仰仁等[9]利用水分變化對小麥的影響規(guī)律首次提出了水分敏感指數(shù)。司昌亮等[10]使用玉米膜下滴灌非充分灌溉試驗數(shù)據(jù)探究Jensen、Minhas、Blank、Stewart、Singh 5 種水分生產(chǎn)函數(shù)模型，得出了各模型水分敏感指數(shù)。李中愷等[11]通過篩選田間數(shù)據(jù)，修正了不同地區(qū)獲得的小麥水分生產(chǎn)函數(shù)，結(jié)果表明發(fā)展機制模型是作物水分生產(chǎn)函數(shù)研究最重要的部分。鄔佳賓等[12]、張淑杰等[13]結(jié)合蒸散量進行了不同地區(qū)作物生長需水量等研究，得到了不同種類作物的作物系數(shù)?；粜堑萚14]、李生勇等[15]、蔣磊等[16]分別利用不同的水分生產(chǎn)函數(shù)模型估算了河套灌區(qū)不同作物的產(chǎn)量，結(jié)果證明利用水分生產(chǎn)函數(shù)進行作物估產(chǎn)在干旱和半干旱區(qū)具有較好的結(jié)果。

水分生產(chǎn)函數(shù)的重要參數(shù)之一是作物在各個生育期內(nèi)的蒸散量，目前主要通過地表水分平衡或者地表能量平衡來估算作物在生育期蒸散量[17-18]。現(xiàn)通過一定時序的NDVI數(shù)據(jù)反演河套灌區(qū)主要3種作物的種植結(jié)構(gòu)，首次利用改進線性雙源遙感蒸散模型和水分生產(chǎn)函數(shù)進行主要作物玉米估產(chǎn)，分析統(tǒng)計作物不同生育期的蒸散量與最大生產(chǎn)潛力，從而有效節(jié)約農(nóng)業(yè)用水并提高經(jīng)濟效益，為研究區(qū)科學(xué)灌溉提供依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

河套灌區(qū)位于中國內(nèi)蒙古自治區(qū)中西部。該地區(qū)為溫帶大陸性氣候，日照充足，作物年可一熟。灌區(qū)劃分為五原縣、臨河區(qū)、杭錦后旗和磴口縣4個區(qū)縣。因磴口縣所種植與本研究相關(guān)的作物較少，因此主要研究區(qū)域為杭錦后旗、臨河區(qū)和五原縣(圖1)。灌區(qū)全年平均氣溫為7～9 ℃，氣溫由南向北、由西向東遞減，1月份平均氣溫達到全年最低，達到-15～-6 ℃，7月份平均氣溫達到全年最高，達到23～26 ℃。

圖1 河套灌區(qū)地理位置圖Fig.1 Geolocation of the river cover irrigation area

1.2 數(shù)據(jù)來源

為研究河套灌區(qū)的主要農(nóng)作物種植結(jié)構(gòu)及利用改進線性雙源遙感蒸散模型和常用水分生產(chǎn)函數(shù)模型進行作物估產(chǎn)，用到的數(shù)據(jù)如下。

(1)Landsat8 OLI數(shù)據(jù)。在USGS網(wǎng)站(https：//earthexplorer.usgs.gov/)下載2016年5月30日、6月15日、7月1日、8月2日和9月19日行列號分別為129/31和129/32的分辨率為30 m的十景遙感影像。

(2)MODIS數(shù)據(jù)。包括2016年5—9月(玉米生長期)8 d合成的蒸散產(chǎn)品MOD16A2。

(3)氣象數(shù)據(jù)。溫室數(shù)據(jù)共享平臺(data.sheshiyuanyi.com)的臨河、烏拉特中旗兩個站點的平均氣溫、最大最低氣溫、日照時數(shù)、水氣壓等數(shù)據(jù)。

(4)作物生育期水分生產(chǎn)函數(shù)參數(shù)及產(chǎn)量信息來自于巴彥淖爾市農(nóng)牧局(https：//nmj.bynr.gov.cn/)和統(tǒng)計年鑒數(shù)據(jù)。

1.3 研究方法

1.3.1 作物識別

NDVI具有對植被變化敏感性高、適用于大范圍的植被監(jiān)測等優(yōu)點，在作物不同生育期具有較大的差異。因此可充分根據(jù)不同作物在一定時序的生長期內(nèi)光譜變化、紋理變化各不相同的特征，建立對應(yīng)的訓(xùn)練樣本庫來區(qū)分不同作物，再根據(jù)同種作物在不同時相的遙感影像上具有不同的光譜響應(yīng)及紋理差異，并結(jié)合一定時序的NDVI對作物進行分類。NDVI可用近紅外波段(near infrared，NIR)和紅波段(red，R)表示為

NDVI=(ρNIR-ρR)/(ρNIR+ρR)

(1)

式(1)中：ρNIR為近紅外波段的反射率；ρR為紅波段的反射率。

采用最大似然分類法將土地覆蓋類型分為農(nóng)用地、水體、建筑用地和裸地4類，其中裸地主要包括裸荒地、鹽荒地等沒有地物覆蓋的土地。由于研究區(qū)小麥、葵花和玉米3種主要作物具有不同的生育期(表1)，選擇一景影像區(qū)分作物植被較為困難，因此使用了5景Landsat8 OLI影像對作物覆蓋進行分類。在經(jīng)過輻射校正和大氣校正后的影像，植被覆蓋率高于30%的地表光譜反射率才能呈現(xiàn)出植被的光譜特征，NDVI的臨界值約為0.3[6]，結(jié)合研究區(qū)經(jīng)驗閾值，本文將像元有無作物植被覆蓋的NDVI臨界值選取為0.35。由于在多時相的遙感影像上具有不同作物的物候信息，利用不同作物的物候差異，結(jié)合經(jīng)過去除異常值后的時間序列NDVI數(shù)據(jù)進行農(nóng)作物的精細分類。

對研究區(qū)進行作物識別分類，重點是確定不同種類作物在不同生育期的NDVI閾值。而在不同的月份，不同的作物也處在不同的生長期，例如在五月份玉米進行播種沒有長苗，也要將其劃定到玉米這一類別中。研究區(qū)小麥4月上旬播種，5月下旬和6月上旬處于拔節(jié)期，NDVI值高于其他兩種作物，7月上處于乳熟期，NDVI值較高，而8月下旬已經(jīng)被收割，因此NDVI值低于其他兩種作物?？úシN時間相對較晚為5月下旬或6月上旬，NDVI值低于其他兩種作物；開花期在8月上旬；9月下旬處于成熟期，NDVI值最大。玉米5月下旬和6月上旬處于出苗期，NDVI值高于葵花低于小麥，7月上處于拔節(jié)期，8月下旬處于乳熟期，NDVI值較高；9月下旬處于成熟期。因此可根據(jù)NDVI經(jīng)驗閾值結(jié)合不同作物生長期對作物進行分類(圖2)。

1.3.2 構(gòu)建線性雙源遙感蒸散模型

地表凈輻射(Rn)代表地表凈收入或凈支出的輻射量，是遙感蒸散模型的重要參數(shù)之一。其中太陽輻射值Ra是Rn的重要組成部分，計算公式[19-20]為

(2)

式(2)中：Ra為太陽輻射值；Gsc為太陽常數(shù)；dr為日地距離系數(shù)；ωs為太陽時角；φ為緯度；δ為太陽赤緯。

篩選出研究區(qū)臨河、烏拉特中旗兩個氣象站點數(shù)據(jù)，Rn計算公式為

(3)

表1 3種主要作物生育期Table 1 Fertility periods for three major crops

注：5.30NDVI、6.15NDVI、7.1NDVI、8.2NDVI分別代表5月30日、6月15日、7月1日和8月2日的作物NDVI值圖2 3種作物NDVI閾值分類Fig.2 Three crop NDVI threshold classifications

(4)

Rs0=(0.75+2×10-5Z)Ra

(5)

式中：α為地表反射率;σ為斯蒂芬-波爾茲曼常數(shù)，取值4.903×10-9MJ/(k4·m2);Tmin，k4為月平均最低溫度，K；Tmax，k4為月平均最高溫度，K;ea為實際水汽壓，kPa；n為實際日照時數(shù)，h;N為理想狀態(tài)下日照時數(shù)，h;as為陰天地表短波輻射與天文輻射的比例系數(shù)；bs為晴天地表短波輻射與天文輻射的比例系數(shù)；Z為高程，m;Rs0為晴天作物表層短波輻射，MJ/m2;Rs為作物表層短波輻射，MJ/m2。

大量研究表明NDVI與植被覆蓋度(fv)具有很強的相關(guān)性，且fv是地表總蒸散量(ET)的主要影響因素之一[21]。使用NDVI表示植被覆蓋率fv，將總蒸散量分成植被蒸散量(ETv)和土壤蒸散量(ETs)兩部分，結(jié)合Nishida等[22]提出的線性雙源遙感蒸散模型，并參考柳錦寶[23]提出的改進的線性雙源遙感蒸散模型的建模及簡化參數(shù)過程，具體簡化的公式為

ET=fvETv+(1-fv)ETs

(6)

假設(shè)植被蒸騰最主要的影響因素是地表凈輻射Rn，其他的影響因素包括平均溫度T、空氣溫差的倒數(shù)1/(Tmax-Tmin)，引入回歸系數(shù)b1和b2簡化植被蒸騰模型，即

(7)

式(7)中：Tmax-Tmin為空氣晝夜溫差，可直接從站點數(shù)據(jù)獲取?？諝鉁夭畹箶?shù)1/(Tmax-Tmin)同時可表示土壤含水量，再引入回歸系數(shù)b3，優(yōu)化土壤蒸發(fā)部分的表達式為

(8)

因為Rn是總蒸散量最主要控制因素之一，假定b0Rn為整個表達式的更正項。結(jié)合以上方程將總蒸散量ET表示為

(9)

結(jié)合NDVI對蒸散量和fv的影響，利用NDVI繼續(xù)簡化模型參數(shù)公式為

(10)

式(10)中：λ0、λ1、λ2、λ3均為回歸系數(shù)。將氣象站點數(shù)據(jù)以及經(jīng)過MRT軟件投影轉(zhuǎn)換和重采樣的MOD16A2蒸散產(chǎn)品數(shù)據(jù)代入模型，利用狼群算法[24]求取回歸系數(shù)。

1.3.3 水分生產(chǎn)函數(shù)

幾種水分生產(chǎn)函數(shù)主要思想是利用實際產(chǎn)量(Xa)與理想條件下產(chǎn)量 (Xb) 的比值與實際蒸散量(ETa)和理想條件下蒸散量(ETb)的比值兩者之間的相互關(guān)系建立函數(shù)方程。在前人的研究基礎(chǔ)上，選擇Jensen模型和Stewart模型對研究區(qū)的作物進行估產(chǎn)。其中Jensen模型公式為

(11)

Stewart模型的公式為

(12)

式中：N為作物生育期階段總數(shù)；λi為第i階段作物水分敏感指數(shù)，一般來說水分敏感指數(shù)不能為負值，負值代表著缺水對作物具有增產(chǎn)的效果[8]。

2 結(jié)果與分析

2.1 主要作物空間分布特征

基于時序NDVI的決策樹分類法主要依靠的是玉米、小麥及葵花3種主要作物的時序NDVI平均曲線特征的差異，同時也要分析它們在同一時期由于長勢不同作物本身的NDVI值變化范圍[6]。將Landsat8 OLI影像不同月份不同種類作物的NDVI值與其所處的生育階段相對應(yīng)，在遙感影像上對作物種類進行監(jiān)測和識別，得到不同作物的空間分布特征(圖3)。對分類后結(jié)果進行斑塊去除、分類疊加、柵格轉(zhuǎn)換、地物掩膜等處理。最后將分類結(jié)果進行各種作物面積統(tǒng)計及精度評定，得到這3種主要作物的面積占比，經(jīng)過精度評定得到2016年影像分類的Kappa系數(shù)為0.953 3，表明分類結(jié)果較好。

其中，研究區(qū)主要作物中玉米的種植量最多，葵花其次，小麥最少。數(shù)理統(tǒng)計分析得臨河區(qū)的玉米種植面積最大約為47 798 hm2，杭錦后旗其次約為42 178 hm2，五原縣最少約為27 425 hm2。3個地區(qū)遙感反演玉米種植面積與實際統(tǒng)計數(shù)據(jù)對比如表2所示。研究區(qū)主要作物中小麥的種植分布相對較為分散，大面積連片的種植較少，種植面積旗縣分布不均勻，總體呈西多東少的趨勢。葵花的種植也相對分散，主要分布臨河區(qū)和五原縣，總體上呈現(xiàn)南多北少。玉米的種植較為廣泛，各個旗縣均有大面積分布。

圖3 主要作物種植結(jié)構(gòu)Fig.3 Major crop planting structures

表2 3個地區(qū)種植玉米面積統(tǒng)計Table 2 Statistics on the area of corn cultivated in the three regions

2.2 日尺度下玉米4個生長期蒸散量

玉米是研究區(qū)種植量最多的作物，一般是在4月下旬或5月上旬播種，9月下旬收割。按照不同月份的影像數(shù)據(jù)，將玉米的生育期劃分為初期、中期、旺盛期和末期4個階段，由改進線性雙源遙感蒸散模型得玉米在4個生長期的蒸散量(圖4)，通過分析得知在玉米生長4個時期1 d的平均蒸散量分別為2.00、3.09、4.20和2.79 mm，再結(jié)合研究區(qū)玉米4個時期的平均生育期天數(shù)和總生育期天數(shù)，得到玉米在4個生育期蒸散量分別為100、154.5、168和83.7 mm，研究區(qū)2016年玉米蒸散總量約為506.20 mm。

分別提取4個時期同一經(jīng)緯度下的預(yù)測蒸散值和經(jīng)過預(yù)處理的MODIS 實際蒸散產(chǎn)品MOD16A2蒸散值各50個代表點來驗證模型精度，將兩者進行擬合得到4個時期內(nèi)R2分別為0.660 1、0.761 7、0.770 1和0.721 9，證明兩者在生育期內(nèi)均具有較強相關(guān)性(圖5)。且4個生育期內(nèi)蒸散預(yù)測值與實際值擬合結(jié)果均通過了0.05顯著性水平的驗證。時間尺度上，線性雙源遙感蒸散模型擬合最好的是作物生長中期和旺盛期，其次是末期，雖然初期擬合效果最差，但R2仍高于0.6，因此在整個生育期內(nèi)反演結(jié)果均具有一定的可信度。空間尺度上，杭錦后旗大部分地區(qū)和臨河部分地區(qū)在作物生育期蒸散量較高，而五原縣相對較低，呈現(xiàn)從東到西遞減的趨勢；在作物生長旺盛期時擬合程度最好，未形成相關(guān)性系數(shù)值偏低區(qū)域；初期和末期時臨河縣中部地區(qū)和五原縣北部部分地區(qū)相關(guān)性系數(shù)擬合程度相對較低。

圖4 玉米4個時期蒸散量Fig.4 The amount of maize vapor over four periods

圖5 玉米4個時期蒸散量精度驗證Fig.5 Verification of the accuracy of the evapotranspiration of maize over four periods

由實驗結(jié)果結(jié)合作物生育期可知，在旺盛期即7月中下旬到8月中下旬作物蒸散量最高，占整個生育期的35%，初期最小占整個生育期的17%。一定程度上說明了各個生育期的蒸散量對作物生長狀態(tài)起著重要作用。

2.3 玉米產(chǎn)量估算及精度驗證

將玉米的生育期具體劃分為苗期到拔節(jié)、拔節(jié)到抽雄、抽雄到灌漿、灌漿到成熟4個階段，以臨河區(qū)為主要研究對象，提取玉米在4個不同生長期的蒸散量，查詢統(tǒng)計數(shù)據(jù)知臨河區(qū)歷年玉米平均生長期天數(shù)為156 d，平均潛在產(chǎn)量為12 500.0 kg/hm2，且2016年臨河區(qū)玉米實際產(chǎn)量為10 315.5 kg/hm2。結(jié)合水分估產(chǎn)函數(shù)Jensen模型和Stewart模型在研究區(qū)經(jīng)驗參數(shù)對臨河區(qū)玉米進行估產(chǎn)分析，得到的結(jié)果如表3所示。其中Jensen模型相對誤差為9.35%，Stewart模型相對誤差為 6.82%，Stewart模型模擬效果更好。

表3 兩種模型作物估產(chǎn)結(jié)果Table 3 Results of two model crop estimates

3 結(jié)論

基于Landsat8影像時序NDVI圖像的決策樹分類法對研究區(qū)的主要農(nóng)作物和土地利用分類的總體精度較高，結(jié)果說明這種方法對河套灌區(qū)地物識別和農(nóng)作物空間分布提取和反演比較可靠，該研究方法和結(jié)果為監(jiān)測干旱區(qū)和半干旱區(qū)農(nóng)作物種植結(jié)構(gòu)提供了一定的參考。目前研究區(qū)小麥種植面積較少，主要問題是效益較低，在主要農(nóng)作物中畝收益排最后，研究區(qū)農(nóng)民種植的積極性不高。因此可增加企業(yè)的帶動能力，增加市場占有率和產(chǎn)品附加值，改善研究區(qū)灌溉排水體系，提高收益和調(diào)動種植的積極性。

基于改進線性雙源遙感蒸散模型反演的研究區(qū)蒸散量的值與MODIS實際蒸散產(chǎn)品值相比較，R2在作物生育期內(nèi)均超過0.6，證明兩者具有較強的相關(guān)性，表明此模型適用于干旱地區(qū)大面積尺度的地表蒸散監(jiān)測。并且通過各個生育期內(nèi)蒸散量結(jié)合玉米各個生育期的水分敏感指數(shù)，初步估算出研究區(qū)玉米在該年份的產(chǎn)量，相對誤差最小為6.82%，對在缺乏地面實驗數(shù)據(jù)的情況下作物估產(chǎn)提出了一種新的思路。

利用遙感數(shù)據(jù)反演出的作物種植結(jié)構(gòu)結(jié)合改進線性雙源遙感蒸散模型進行研究區(qū)的作物估產(chǎn)，此種方法受作物反演精度的影響較大。研究區(qū)作物套種現(xiàn)象較為普遍，在農(nóng)用地范圍內(nèi)混合像元較多，對監(jiān)測作物生育期NDVI變化具有一定影響。在定性分析的基礎(chǔ)上，結(jié)合影像分辨率進行數(shù)理統(tǒng)計及作物自動識別分類精度相對較低，后續(xù)應(yīng)在更高分辨率的影像上展開進一步的研究。