基于線性雙源遙感蒸散模型的作物估產
——以河套灌區為例

杜鑫， 蘇濤*

(1. 安徽理工大學空間信息與測繪工程學院， 淮南 232001; 2. 安徽理工大學礦山采動災害空天地協同監測與預警安徽普通高校重點實驗室， 淮南 232001; 3. 安徽理工大學礦區環境與災害協同監測煤炭行業工程研究中心， 淮南 232001)

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

水分生產函數的重要參數之一是作物在各個生育期內的蒸散量，目前主要通過地表水分平衡或者地表能量平衡來估算作物在生育期蒸散量[17-18]?，F通過一定時序的NDVI數據反演河套灌區主要3種作物的種植結構，首次利用改進線性雙源遙感蒸散模型和水分生產函數進行主要作物玉米估產，分析統計作物不同生育期的蒸散量與最大生產潛力，從而有效節約農業用水并提高經濟效益，為研究區科學灌溉提供依據。

1 材料與方法

1.1 研究區概況

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

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

1.2 數據來源

為研究河套灌區的主要農作物種植結構及利用改進線性雙源遙感蒸散模型和常用水分生產函數模型進行作物估產，用到的數據如下。

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

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

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

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

1.3 研究方法

1.3.1 作物識別

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

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

(1)

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

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

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

1.3.2 構建線性雙源遙感蒸散模型

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

(2)

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

篩選出研究區臨河、烏拉特中旗兩個氣象站點數據，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)

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

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

ET=fvETv+(1-fv)ETs

(6)

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

(7)

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

(8)

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

(9)

結合NDVI對蒸散量和fv的影響，利用NDVI繼續簡化模型參數公式為

(10)

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

1.3.3 水分生產函數

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

(11)

Stewart模型的公式為

(12)

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

2 結果與分析

2.1 主要作物空間分布特征

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

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

圖3 主要作物種植結構Fig.3 Major crop planting structures

表2 3個地區種植玉米面積統計Table 2 Statistics on the area of corn cultivated in the three regions

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

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

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

圖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

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

2.3 玉米產量估算及精度驗證

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

表3 兩種模型作物估產結果Table 3 Results of two model crop estimates

3 結論

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

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

利用遙感數據反演出的作物種植結構結合改進線性雙源遙感蒸散模型進行研究區的作物估產，此種方法受作物反演精度的影響較大。研究區作物套種現象較為普遍，在農用地范圍內混合像元較多，對監測作物生育期NDVI變化具有一定影響。在定性分析的基礎上，結合影像分辨率進行數理統計及作物自動識別分類精度相對較低，后續應在更高分辨率的影像上展開進一步的研究。