基于時序MODIS和氣象數據的日喀則地區(qū)青稞估產

潘瑩，仙巍，陳春容，陳軍，徐薇

(1.成都信息工程大學 資源環(huán)境學院，成都 610225；2.日喀則市氣象局，西藏 日喀則 857000)

0 引言

青稞是藏區(qū)人民的主食，也是海拔4 500 m以上唯一可以正常成熟的糧食作物。青藏高原是全球青稞產量和種植面積最大的地區(qū)，而日喀則地區(qū)作為西藏的第二大城市，也是青稞的主產區(qū)。因此，快速準確地監(jiān)測青稞的產量，對藏區(qū)的產業(yè)發(fā)展與經濟發(fā)展具重要意義。

當前，對青稞估產的相關研究較少，而估產可分為作物的種植面積提取與估產模型建立兩部分。日喀則地區(qū)地形復雜，青稞種植區(qū)域分布零碎，傳統的調查估產手段難以得到全面準確的結果，但遙感技術可以解決其無法全覆蓋和成本過高的問題，具有較好的準確度和時空連續(xù)性的優(yōu)勢，在種植面積提取和產量預測中起著重要的作用[1]。利用遙感手段實現產量預測方法主要有統計預測法和生長模型預測法兩大類，其中，生長模型預測法是數值模擬作物的生長過程，需引入大量的影響作物生長的參數與同化后的遙感數據，故在大范圍的作物估產中難以實現[2-3]；而統計預測法只需選入少量的影響產量的相關參數進行估產，如氣象數據和遙感指數，對大面積的區(qū)域估產相對較易實現。遙感數據估產原理是利用作物的光譜反射特征提取其潛在的產量信息，并且光譜植被指數可以定量化，其中最為常見的歸一化植被指數(normalized difference vegetation index，NDVI)光譜組成與作物產量的物理關系在機理上已得到證實[4-6]。劉寶鋒等[7]運用遙感分析法，結合Landsat TM數據和歷年青稞產量數據，建立了甘孜縣青稞估產模型并取得了較好的結果。劉新杰等[8]基于時序衛(wèi)星的定量NDVI遙感數據實現了冬小麥的長勢監(jiān)測和精確估產。通過遙感技術進行農作物估產主要是建立經驗模型，而回歸模型是典型的代表，王長耀等[9]、任建強等[10]利用MODIS植被指數產品，采用逐步回歸法建立農作物估產模型，取得了較高的估產精度。

氣候變暖正影響著脆弱的高原生態(tài)環(huán)境，青藏高原作為中國乃至世界氣候變化的敏感區(qū)和啟動區(qū)，農業(yè)生產嚴重依賴于氣候資源[11-12]。日喀則地區(qū)的青稞產量對氣候變化存在明顯的響應，特別是對氣溫和降水量異常敏感[13]，而目前基于單純的遙感數據的產量預測模型較多，因此本文結合時間序列的MODIS遙感數據和氣象數據，擬建立一種適用于青藏高原日喀則地區(qū)的綜合的青稞遙感氣象估產模型。

1 數據與方法

1.1 研究區(qū)概況

本文以我國西藏日喀則地區(qū)為研究區(qū)，其地處青藏高原東南部，位于82°00′E～90°07′E、27°55′N～31°16′N，面積為18.2×104km2，地形復雜多樣，平均海拔在4 000 m以上。日喀則地區(qū)空氣稀薄，太陽輻射強，日照長，氣溫低，年較小，日較大，常年農作物播種面積為2.6×105ha。青稞是日喀則地區(qū)播種面積最大的農作物，主要集中在雅魯藏布江和年楚河流域，由拉孜—仁布寬谷和江孜—日喀則平原組成，還包括喜馬拉雅山脈北側、藏南高原上的朋曲河谷平原和其他零星的河谷平原。這些谷地坡度平緩，土層深厚，氣候宜人，水源較充足，有利于農作物的種植。

1.2 數據準備

本文所需數據有遙感數據和非遙感數據。其中遙感來源于美國國家航天局(National Aeronautics and Space Administration，NASA)的官方網站(https：//ladsweb.modaps.eosdis.nasa.gov)，獲取2015—2019年16 d合成的MOD13A1和8 d合成的MOD09A1的NDVI，分別用于青稞的種植區(qū)提取與估產。獲取行列號為h25v05與h25v06，其中MOD13A1共230景，MOD09A1共460景。使用MRT輸入命令行可對以上數據進行重采樣和鑲嵌的批處理。非遙感數據包括氣象數據和統計數據，分別來源于中國氣象數據網(http：//data.cma.cn/)和《日喀則統計年鑒》，包括日喀則地區(qū)2016—2019年各市(縣)的日平均氣溫、日降水量、青稞單產產量和該地區(qū)的播種總面積。

2 研究方法

2.1 青稞種植區(qū)提取

隨著農作物的物候期更替，近紅外波段因作物的葉面積系數和覆蓋率增大而增大，紅光反射率則逐漸減小，因此在某一時間斷面上，日喀則地區(qū)的青稞和冬小麥的狀態(tài)和結構會因物候期的不同而有明顯差別[14-15]。日喀則地區(qū)的冬小麥比青稞要早一個月進入生長最茂盛的時期即拔節(jié)期，利用這種不同作物間物候差異特征，通過MOD13A1的NDVI數據，經過反復的探究與實驗，得到式(1)所示的日喀則地區(qū)青稞種植區(qū)提取方法。該方法不僅可以提取青稞的種植區(qū)，還能有效剔除日喀則地區(qū)以冬小麥為主的其他農作物與綠色植物的影響。

(1)

式中：QK是青稞種植范圍，當QK為1時，為種植區(qū)域，當QK為0時，為非種植區(qū)；NDVI拔節(jié)為拔節(jié)期的NDVI；NDVI播種前為播種前的NDVI。

為提高青稞種植范圍的提取精度，將利用式(1)提取的結果與高清的天地圖進行人工目視解譯對比，剔除個別未在耕地上的錯誤值，最終得到日喀則市青稞2015—2019年的種植面積。

將青稞提取總播種面積與統計總播種面積相比，2019年的提取精度最高，為99.51%，2017年提取精度最低，為97.46%(表1、表2)。歷年來，位于年楚河流域的日喀則市、江孜縣和白朗縣青稞播種面積較多，約占日喀則地區(qū)青稞播種總面積的1/2；亞東縣的播種面積是所提取的縣中最少的，其中2018年和2019年由于青稞播種面積過小，難以提取。分析統計數據發(fā)現，2015年總播種面積最小，為1 293.1 km2，在2018年達最大，為1 395.8 km2，2019年總播種面積為1 391.90 km2。

表1 日喀則地區(qū)青稞種植面積提取結果 km2

表2 日喀則地區(qū)青稞種植面積提取精度

2.2 青稞產量預測方法

1)時間序列分析法。時間序列分析法是一種基于數理統計與隨機過程理論的方法，研究的是連續(xù)時間的動態(tài)數據所遵循的統計規(guī)律[16-19]。青稞生長過程曲線的變化可以反映其長勢的變化，而長勢會影響最終的產量[20]。利用時間序列植被指數曲線分析方法，提取遙感波段反射率，根據青稞生長的時相變化特點可確定遙感監(jiān)測的周期，從而選擇恰當的遙感數據進行估產(圖1)。植被指數時間序列濾波法是去除噪聲的一種常用方法，利用S-G(savizky-golay)濾波法可去除云和大氣對NDVI的影響，可得到去噪后的時間序列NDVI數據。

圖1 青稞時間序列植被指數曲線圖

2)產量預測模型。農作物的生產是一個復雜的綜合系統，不僅會受到農業(yè)生物、技術、經濟的影響，也會受到各氣象條件等環(huán)境因素的影響。分解模型法是一種同時考慮農作物產量的構成因素(農業(yè)技術、農業(yè)環(huán)境和隨機因素等)以及這些因素對產量的影響的方法。模型分解法可將農作物產量分解為趨勢產量、波動產量與隨機產量，如式(2)所示[21]。趨勢產量為農業(yè)技術進步導致的產量部分，波動產量為農作物的自然環(huán)境導致的產量部分，隨機產量是由突發(fā)自然災害和其他不確定因素導致的產量部分(隨機產量只能在模型外單獨處理，本研究不予考慮)。趨勢產量可通過區(qū)域產量序列擬合得出，表現為時間的正函數。波動產量為實際產量與趨勢產量之差，可通過波動產量與氣象參數進行相關性分析，以及反映作物生長情況的植被指數進行時間序列分析，并與產量進行相關系數分析，找出影響波動產量狀況的遙感氣象參數，建立波動產量預測模型，再結合趨勢產量，建立單產預測模型。

YE=XT+XF+E

(2)

式中：YE為預測單產；XT為趨勢產量；XF為基于遙感氣象因子的波動產量；E為隨機產量。

3 青稞產量預測模型的建立

3.1 植被指數時序特征分析

利用8 d的NDVI數據(MDO09A1)監(jiān)測青稞的生長情況，提取日喀則地區(qū)各縣(市)的青稞種植區(qū)所對應的NDVI，獲得日喀則地區(qū)2015—2019年原始NDVI時間序列數據集，利用S-G濾波法去除云、氣溶膠等對時間序列NDVI的影響，得到S-G去噪后的NDVI時間序列數據集并繪制對應的時間序列變化圖。以江孜縣為例(圖2)。從圖2可知，S-G去噪后的NDVI圖像的噪聲明顯減少，其時間序列的NDVI值變化更穩(wěn)定。分析去噪后的NDVI時間序列曲線可知，每年都具有一個明顯的大波峰，峰值出現時間與日喀則地區(qū)青稞的抽穗期相吻合。在整個時間序列上，日喀則地區(qū)的青稞抽穗期NDVI值相對穩(wěn)定，基本為0.6左右，說明青稞種植生長情況基本較穩(wěn)定。

圖2 日喀則地區(qū)江孜縣NDVI時間序列變化圖

3.2 植被指數的相關系數分析

NDVI反映的是作物生長的總體情況，與作物的產量有直接的關系，因此在進行相關性分析時，應該將NDVI與產量直接進行相關分析，然后將NDVI作為影響因素，引入波動產量的預測模型中，作為NDVI對產量的影響。圖3是研究區(qū)時序NDVI與單產的相關系數隨時間的變化曲線，當青稞種植區(qū)NDVI值越高，說明青稞長勢越好，所對應的單產也越高。通過時序的NDVI與單產的相關系數曲線，不僅可以看出何時的NDVI與單產的相關系數最高，也能了解日喀則地區(qū)青稞整體的生長狀況。

圖3 NDVI與單產的相關系數時序變化曲線

由圖3可知，NDVI與單產的相關系數隨時間變化總體呈現拋物線形狀，呈先上升后下降的趨勢。在時間序列曲線的起始部分，NDVI與單產的相關系數較低，這是因為在3月份青稞還沒有開始播種。而相關系數曲線大約在4月底五月初開始出現較大幅度的上升，即在137至161天之間，青稞播種之后開始出苗、分蘗、拔節(jié)、抽穗、孕穗，是青稞生長的關鍵時期，其NDVI值增大，相關系數也增大，對于單產多的地區(qū)，其青稞長勢自然也更好。而在169至233天之間(6月中旬至8月中旬)，日喀則市整體的NDVI值達到飽和并處于相較穩(wěn)定的階段，故其與單產的相關系數也比較穩(wěn)定。觀察該曲線，其最大值出現在209天(7月下旬)，即青稞的抽穗期，而日喀則地區(qū)的青稞也在7月下旬時生長達到最旺盛的階段，所以青稞抽穗期的NDVI是日喀則市青稞估產的最優(yōu)選擇。

3.3 氣象因子的確定

日喀則地區(qū)的青稞產量對氣溫和降水具有較好的響應特征。統計日喀則地區(qū)各市(縣)青稞生長階段的月積溫、月平均氣溫和累計降水量，并提取青稞生長重要物候期內的積溫，將以上氣象因子與波動產量進行相關性分析，計算結果如表3所示。選取與波動產量相關性最大的一個氣象指標可作為構建農業(yè)氣象估產模型的主要因子。由表3可知，日喀則地區(qū)的7月上旬積溫與青稞波動產量的相關性最高，因此7月上旬積溫可作為構建農業(yè)氣象估產模型的主要因子。進一步分析表3發(fā)現：波動產量與積溫和平均氣溫的相關系數均為負值，這是因為青稞屬于喜涼作物，溫度的升高會對青稞產量起到負效應[22]；青稞生長的關鍵期，累積降水與波動產量的相關系數也為負數，可能是因為降水過多時云量大、太陽輻射也減少，而當地過高的土壤含水量導致土壤蒸散增加會使近地表溫度降低，進而降低了青稞的光合作用能力降低，從而對青稞產量起到負效應[23]。

表3 日喀則地區(qū)各氣象因子與波動產量的相關關系

3.4 產量預測模型

選取最為敏感的遙感氣象參數作為構建青稞遙感氣象估產模型(式(2))的主要因子，并檢驗參數的穩(wěn)定性，從而建立一元線性或多元線性回歸方程。波動產量的回歸方程可以分為三種(表4)：氣象波動產量回歸方程、遙感波動產量回歸方程和遙感氣象波動產量回歸方程。其中，XF-m為只考慮氣象因子影響的波動產量；XF-rs為只考慮NDVI的遙感因子的遙感波動產量；XF為既考慮氣象因子又考慮NDVI遙感因子的遙感氣象波動產量；Xm為氣象因子中的7月下旬的積溫；Xrs為遙感因子的7月下旬的NDVI(MDO09A1的209天的NDVI)。

表4 三種波動產量回歸方程

基于以上三種計算波動產量的回歸方程，然后根據2019年的遙感氣象數據計算2019年日喀則地區(qū)各縣中青稞的遙感氣象波動產量，并基于式(2)結合趨勢產量得到青稞單產的預測結果(表5)。

表5 2019年日喀則地區(qū)各縣青稞單產預測結果 kg·hm-2

3.5 誤差分析

將預測的青稞單產與統計數據進行比較和相關性分析，驗證模型的顯著性和擬合精度。由表6可知：建立的青稞產量預測模型中，遙感氣象估產模型的預測結果最為穩(wěn)定，誤差最大的絕對值為16.19%，相對誤差控制在-14.34%～16.19%；氣象估產模型誤差最大的絕對值為19.66%，相對誤差范圍控制在-9.59%～19.66%；遙感估產模型誤差最大的絕對值為19.14%，相對誤差控制在15.41%～19.14%；平均估產誤差(即各市(縣)所有模型的相對誤差的絕對值之和的平均)為7.53%。由圖4可知，氣象、遙感、遙感氣象預測單產與統計單產的相關系數分別為0.943 6、0.946 1、0.938 2，這說明三種預測模型得到的預測單產與統計單產都具有良好的線性關系。

表6 2019年日喀則地區(qū)各縣青稞單產預測誤差 %

圖4 2019年日喀則地區(qū)青稞預測單產與統計單產對比

4 結束語

本文提出了一種利用MOD13A1快速準確地提取青稞的方法，2015—2019年青稞的提取精度均大于97.4%。日喀則地區(qū)青稞主要分布在雅魯藏布江和年楚河流域的河谷平原，西北降水少，海拔高，農田很少，提取結果與高子恒等[24]得到的日喀則地區(qū)的農田結果相似，具有一定的可靠性。

日喀則地區(qū)的遙感估產的敏感因子為7月下旬(抽穗期)的NDVI(MDO09A1的209天的NDVI)?；跁r間序列分析，對日喀則地區(qū)的青稞進行長勢監(jiān)測，近五年青稞生長曲線波峰的NDVI都在抽穗期，穩(wěn)定在0.6左右，說明青稞種植生長情況穩(wěn)定；青稞NDVI與單產的相關系數在生長期呈拋物線，7月下旬(抽穗期)達到最大值0.709，青稞抽穗期的NDVI是日喀則地區(qū)青稞估產的最優(yōu)選擇。

日喀則地區(qū)氣象估產的敏感因子為7月上旬的積溫。日喀則地區(qū)青稞產量對氣候變化存在明顯響應，且對氣溫的敏感性大于降水量。

本文構建的日喀則地區(qū)青稞單產預測模型的平均估產誤差為7.53%，平均相關系數為0.942 6，在一定程度上能滿足當地的青稞估產需求。利用模型分解法可建立青稞的氣象單產預測模型、遙感單產預測模型和遙感氣象單產預測模型。其中，遙感氣象單產預測模型最穩(wěn)定，最大的相對誤差控制16.19%以內。所有預測結果與統計單產都具有良好的線性關系，相關系數均大于0.938 0。因此，本文所構建的估產模型在日喀則地區(qū)青稞生產區(qū)具有一定的可操作性和適用性。

由于日喀則地區(qū)的范圍較廣，氣象數據受到區(qū)域范圍的影響，會對利用氣象因子預測的單產造成一定的誤差；而利用遙感數據提取各年份的青稞總播種面積精度均大于97.4%，但由于混合像元的原因，也會對利用遙感因子預測的單產造成誤差。另外，個別縣(市)的青稞單產數據年際差別較大，也是造成預測結果誤差偏大的原因，如定結縣、崗巴縣、拉孜縣和日喀則市。日喀則地區(qū)的遙感和氣象估產的敏感因子分別為7月下旬即青稞抽穗期的NDVI(MDO09A1的209天的NDVI)和7月上旬的積溫。氣象因子中積溫對波動產量的影響最大，不同的氣溫累積會直接影響物質的累積與蛋白質、淀粉的含量，從而影響產量[25]。日喀則地區(qū)處于高海拔的青藏高原西南部，在青藏高原強烈的選擇壓力下，海拔和系統發(fā)育綜合因子會對植物千粒重產生影響，進而會對青稞產量造成一定的影響。日喀則地區(qū)的青稞種植區(qū)的海拔落差大約有1 700 m，由于較大的海拔差異，對溫度影響顯著，主要表現隨著海拔的升高，溫度降低，作物生長季短，因此會造成一定的物候時間差異，影響產量預測的結果[26]。

本文彌補了當前針對高原地區(qū)大范圍青稞估產的空白，構建了適應于當前發(fā)展趨勢的綜合遙感氣象農作物估產模型，對高原地區(qū)的糧食生產與經濟發(fā)展有重要的指導意義。在模型的實際應用中，需要對不同海拔高度的縣進行實地調查，利用抽樣統計可對作物單產模型進行修正，從而提高模型精度。