基于同化葉面積指數和條件植被溫度指數的冬小麥單產估測

張樹譽，孫輝濤，王鵬新，景毅剛，李 俐

(1.陜西省氣象局， 陜西 西安 710014； 2.中國農業大學信息與電氣工程學院， 北京 100083)

傳統的作物估產主要采用統計調查、農業氣象預報等方法，成本偏高、效率較低，難以實現區域尺度高精度作物產量的估測[1]。目前，作物生長模型已成功應用于單點尺度的作物生長發育過程模擬和產量估測，但由于難以獲取能夠準確描述大范圍非均勻性地表、近地表環境空間的模型輸入參數，使其難以應用于區域尺度[2]；而衛星遙感具有及時、宏觀、信息量大等特點，并能夠定量地監測和描述作物在區域尺度的生長狀況和反映環境因子對作物生長的綜合影響[3]。數據同化方法能夠有效融合作物生長模型和衛星遙感觀測信息，對現代化農業的發展有著十分重要的意義，已經成功應用于區域產量估測[4-5]。

同化變量(參數)的選取對同化結果精度至關重要，基于多變量數據同化策略能夠綜合多因子對作物產量共同作用的影響，已經成為農業數據同化的研究趨勢[6-7]。葉面積指數(LAI)是評估作物籽粒潛在產量的重要指標[8]，土壤水分(SM)由于控制著作物水分的脅迫信息，與作物產量變化密不可分[9]。Ines等[10]采用改進的集合卡爾曼濾波(EnKF)算法對CERES-Maize模型模擬和遙感數據觀測的LAI和SM實施同化，結果表明，同時同化LAI和SM雙變量的模型估測結果精度更高。Huang等[11]采用四維變分(4DVAR)算法對SWAP模型模擬和遙感數據觀測的LAI和ET實施同化，結果表明，與單獨同化LAI或ET單變量相比，同時同化LAI和EI雙變量明顯提高了冬小麥單產的估測精度。

目前，獲取同化所需的觀測土壤水分的方式主要是微波遙感技術[10，12]，但由于其在高植被覆蓋度條件下對土壤淺層水分的反演精度較低[13]，因此會影響同化結果的精度。王鵬新等[14]在歸一化植被指數(NDVI)和地表溫度(LST)的散點圖呈三角形區域分布的基礎上，提出了條件植被溫度指數(VTCI)的干旱監測方法，并在陜西省關中平原冬小麥的干旱監測、預測及其影響評估等研究中得到了廣泛應用[15-18]。以往的研究表明，VTCI與土壤淺層水分具有顯著相關性[19]，適合高植被覆蓋度下的干旱監測[20]。基于此，本研究以LAI和VTCI為同化系統狀態變量，采用粒子濾波(PF)同化算法對CERES-Wheat模型模擬和遙感觀測的LAI和VTCI實施同化，運用熵的組合賦權方法獲取加權觀測和加權同化變量值，結合樣點實測單產構建基于觀測和同化變量的估產模型，并分別評價其單點和區域尺度的估測結果精度，以期更準確地進行作物長勢監測、干旱監測預測及其影響評估。

1 材料與方法

1.1 試驗區概況

關中平原位于陜西省中部，素有“八百里秦川”之稱，其行政區域包括西安、寶雞、咸陽、渭南、銅川5個省轄市和楊凌國家農業高新技術產業示范區。區域內地勢平坦，土層深厚，蓄水性好，但是區域內水資源相對不足，干旱是造成糧食減產的主要原因[21]。作物種植模式在灌溉區域主要為冬小麥和夏玉米輪作，而在旱作區域多為夏季休閑式的冬小麥連作。冬小麥播種時間一般為10月上、中旬，播種5～6 d后出苗，4月下旬進入抽穗期，乳熟期為5月中、下旬，5月下旬或6月上旬收獲[22]。研究選擇2008—2014年冬小麥生長季內，關中平原12個典型的小麥種植區作為試驗樣點(圖1)，其中眉縣常興鎮、扶風縣城北、三原縣魯橋鎮和臨渭區藺店鎮為灌溉樣點，其余為旱作(雨養)樣點。

圖1研究區域

Fig.1 Study area

1.2 數據來源及預處理

1.2.1 CERES-Wheat模型模擬LAI和VTCI 在CERES系列模型中，CERES-Wheat模型是面向小麥生長和發育過程的決策等級模型。CERES-Wheat模型的輸入數據主要包括氣象數據、土壤數據、田間管理數據和小麥品種遺傳特性參數。除小麥品種遺傳特性參數外，均由實地調研和實測獲取。冬小麥的遺傳特性參數控制著其生長發育進程，直接關系到植株物理形態的發育與作物產量的形成，因此模型在實際使用前先結合當地實際測量的LAI、生物量、單產、收獲日期等對這些參數進行本地化標定[22-23]。

同化所需模擬LAI和VTCI時間序列數據的獲取。運行標定后的模型得到樣點模擬LAI時間序列數據。研究所需的樣點模擬VTCI不能由模型直接模擬得到，但以往的研究結果表明，該地區多年旬尺度的VTCI與土壤淺層(0～20 cm)水分具有顯著相關性[19]，因此本文將觀測VTCI與模擬淺層水分數據進行線性回歸分析，建立經驗關系，獲取樣點模擬VTCI時間序列數據。

1.2.2 MODIS遙感數據反演LAI和VTCI 基于Aqua-MODIS日地表反射率產品(MYD09GA)和日地表溫度產品(MYD11A1)，利用MODIS重投影工具MRT(MODIS re-projection tool)對原始數據進行裁剪、投影轉換、重采樣、格式轉換等預處理操作，得到研究區日LST和日地表反射率數據。利用日地表反射率數據計算日NDVI，應用最大值合成技術分別生成旬NDVI和旬LST最大值合成產品，并依據VTCI計算方法獲取旬VTCI[14,24]。本研究LAI選取2008—2014年冬小麥主要生育期內所有MCD15A3產品，時間分辨率較高，每4 d合成一次，空間分辨率為1 km，更有利于農作物長勢和物候的監測。利用MRT將所有數據轉換為統一的Lambert等面積投影，對投影后的數據作裁剪處理，得到關中平原的遙感觀測LAI。以樣點所在像素為中心3像素×3像素內所有像素的LAI和VTCI的均值作為同化所需觀測LAI和VTCI時間序列數據。

1.3 研究方法

(1)

(2)

根據敏感性分析結果，本研究將粒子數設置為200，對于LAI和VTCI，基于田間實測和以往經驗，CERES-Wheat模型模擬誤差標準差分別設定為0.5和0.1，分別設定遙感觀測的13%和3%作為觀測誤差標準差。

1.3.2 熵的組合賦權方法 將冬小麥越冬后的生育期劃分為返青期(3月上旬—3月中旬)、拔節期(3月下旬—4月中旬)、抽穗～灌漿期(4月下旬—5月上旬)和乳熟期(5月中旬—5月下旬)4個生育時期[18]。應用研究區域2008—2014年冬小麥4個生育時期(n)的LAI和VTCI數據構建數據矩陣A(amn)N×4(N為研究選取的總年份數)，計算各個生育時期的熵值hn[26]：

(3)

計算第n個生育時期的差異性系數gn，并對其進行歸一化處理得到第n個生育時期的權重wn：

(4)

(5)

1.4 時間尺度的匹配與估產模型的構建

由于研究選用的觀測(未同化)和同化VTCI均是以旬為尺度，而觀測LAI是4 d合成產品，同化LAI以1 d為步長。為了解決兩變量時間尺度不一致的問題，在冬小麥的4個生育時期將觀測LAI采用S-G濾波方法[27]插值成以1 d為步長，然后將各生育時期S-G濾波后的觀測LAI和經過PF算法的同化LAI累加求和分別作為該生育時期的累積觀測LAI和累積同化LAI。取各生育時期內所包含的多旬觀測和同化VTCI的均值分別作為該生育時期的觀測和同化VTCI。運用熵的組合賦權方法獲取各生育時期的加權LAI和VTCI，將觀測和加權LAI和VTCI分別與樣點實測單產進行一元線性回歸分析構建冬小麥單產估測模型(研究選用2012年用于單產估測模型的精度驗證，故2012年各樣點的相關數據均未參與模型的構建)。

2 結果與分析

2.1 LAI和VTCI的同化結果與分析

通過對2008—2014年12個樣點的綜合對比分析，以2014年代表性較好的乾縣石牛鄉旱作樣點為例，進行LAI和VTCI的同化結果分析(圖2)。關中平原的冬小麥一般在10月上旬播種，之后進入越冬階段，次年的3月上旬進入返青期，返青期后LAI逐漸增大，至抽穗～灌漿期達到最大。經過PF算法的同化LAI變化趨勢(圖2a)，可以看出，同化LAI最大值出現時間在4月25日左右，在合理范圍內。在曲線上升階段，模擬LAI變化趨勢從3月2日至3月14日左右出現強烈的振蕩現象，由于綜合了遙感觀測LAI變化趨勢，同化LAI變得更為平滑；在曲線下降階段，遙感觀測LAI變化趨勢從5月17日至5月24日左右出現先緩慢下降再緩慢上升又迅速下降的過程，而同化LAI從5月12日至5月18日左右變化趨勢與遙感信息基本一致。由于CERES-Wheat模型是機理模型，在針對研究區進行參數和模型準確標定后，其LAI絕對值更接近于實際情況，根據遙感觀測特點，如果不受到嚴重的云污染影響，遙感LAI變化趨勢可直接反映作物LAI的實際變化情況。因此，經過PF算法的同化LAI變化趨勢與冬小麥實際生長變化情況更為相符。

研究已經證明VTCI是一種實時的干旱監測方法，VTCI的值越小，干旱程度越嚴重，土壤水分含量較小，作物受水分脅迫的程度就愈嚴重；反之，則干旱程度越輕，土壤水分含量較大，作物受水分脅迫程度較輕，其與降水量數據有顯著的相關性[19]。為驗證同化VTCI是否能較觀測VTCI和模擬VTCI更好地與降水量數據結合，將觀測VTCI、模擬VTCI和同化VTCI與旬累積降水量數據對比分析(圖2b)，可以看出，整體上同化VTCI均能較好地表達模型模擬和遙感觀測VTCI，具體表現為：乾縣石牛鄉3月上旬和5月下旬的降水量僅為6.4 mm和5.1 mm，但觀測VTCI分別達到0.76和0.80，屬于不旱范圍[28]，受到模擬VTCI的影響，同化VTCI分別調整為0.34和0.32，調整幅度較大，結果更能反映樣點實際水分脅迫程度。另外3月中旬、4月上旬、4月下旬等同化VTCI在原觀測基礎上均有不同程度的調整，結果能更好地與旬累積降水量數據結合，更加符合關中平原實際干旱情況。

圖2 2014年乾縣石牛鄉葉面積指數(LAI)和條件植被溫度指數(VTCI)變化趨勢

Fig.2 Comparison of seasonal changes of leaf area index (LAI) and vegetation temperature condition index (VTCI) in Shiniu of Qian County in 2014

2.2 估產模型的構建與分析

將2008—2014年(除2012年)每年主要生育時期的累積觀測LAI和VTCI以及累積同化LAI和VTCI分別運用熵的組合賦權方法計算獲取每年加權LAI和VTCI，并將加權變量值分別與樣點每年實測單產進行線性回歸分析，構建基于觀測和同化LAI和VTCI的冬小麥單產估測模型(表1)。結果表明，各估產模型均通過了1%的顯著性水平檢驗，無論是基于觀測還是同化變量，基于LAI和VTCI雙變量構建的估產模型的估測結果精度均高于單獨基于LAI或VTCI單變量構建的結果，且基于同化變量構建的估產模型的估測結果精度明顯優于基于觀測變量構建的結果。其中，基于觀測LAI單變量構建的估產模型的估測結果精度最低(R2=0.279)，可能的主要原因是MODIS遙感觀測LAI受云和混合像素的影響會導致整體偏低，與冬小麥實際情況偏差較大。通過對比分析，可以看出，基于同化LAI和VTCI雙變量構建的估產模型的估測結果精度最高(R2=0.547)，確定其為最優的單產估測模型。

應用基于同化LAI和VTCI雙變量構建的最優單產估測模型對2012年12個研究樣點進行單產估測，并將估測單產與實測單產進行一元線性回歸分析。結果表明，回歸模型的決定系數較高(R2=0.594)，達顯著水平(P<0.01)；估測單產與實測單產間間的均方根誤差為527.10 kg·hm-2，表明應用基于PF算法的同化變量構建的單產估測模型的單點尺度估產結果誤差較低，精度較高。

2.3 區域冬小麥單產的估測與分析

應用最優單產估測模型對2008—2014冬小麥生長年關中平原冬小麥區域單產進行估測(圖3)。結果表明，從年際變化上看，關中平原2008—2014年的單產呈現個別年份波動，但整體平穩增長的趨勢，其中2013年由于旱情嚴重而減產較為嚴重，這與關中平原實際單產及旱情的分布特點較為一致。

表1 基于觀測和同化LAI和VTCI的冬小麥單產估測模型

注：Y表示估測的冬小麥單產，單位為kg·hm-2。 Note:Yrepresents the estimated yield (kg·hm-2).

圖3 2008—2014年關中平原冬小麥單產的估測結果

Fig.3 Estimated yields of winter wheat in the Guanzhong Plain in the years from 2008 to 2014

通過對區域冬小麥估測單產的空間分布(圖3)進行分析，可以看出，關中平原冬小麥單產整體呈現中部單產高于東部和西部，西部單產高于東部的空間分布規律，這與關中平原實際情況相符。綜上所述，基于PF算法的同化LAI和VTCI雙變量構建的估產模型的區域單產估算結果在時空分布上均與關中平原實際情況一致，說明基于PF算法的同化變量構建的估產模型的精度高，適用于區域冬小麥單產的估測。

3 結論與討論

粒子濾波(PF)同化算法通過構建一組具有權重的隨機樣本，以樣本均值代替復雜的后驗概率積分運算，隨著粒子數目的不斷增加，這些粒子的概率密度函數逐漸逼近最優貝葉斯估計的效果，可用于非線性非高斯動態系統[29]。研究應用PF算法對CERES-Wheat模型模擬和遙感數據反演的LAI和VTCI實施同化。結果表明，相比于模型模擬和遙感反演的LAI，同化LAI具有良好的時間和空間連續性，在MODIS-LAI偏低的情況下同化LAI有所提升，且同化LAI變化趨勢更加符合關中平原冬小麥實際變化情況。同化VTCI能綜合表達模擬和觀測VTCI的變化趨勢，且與旬累積降水量的相關性更高，能更準確地反映冬小麥水分脅迫的程度。

運用熵的組合賦權方法分別構建基于觀測LAI和VTCI以及經過PF算法的同化LAI和VTCI的冬小麥單產估測模型。結果表明，基于同化LAI和VTCI構建的估產模型的單產估測精度明顯優于基于觀測LAI和VTCI構建估產模型的單產估測精度，且基于LAI和VTCI雙變量構建估產模型的估測結果優于基于LAI或VTCI單變量構建的結果。采用最優估產模型分別對單點和區域尺度的冬小麥單產進行估測，結果表明，2012年12個樣點的估測單產與實測單產間的相關性達顯著水平，兩者間的均方根誤差為527.10 kg·hm-2；區域冬小麥單產估測結果在時空分布上均與關中平原實際情況相符，說明無論是單點還是區域尺度，基于PF算法同化LAI和VTCI雙變量構建的估產模型的估測精度均較高。

[1] 黃健熙,李昕璐,劉帝佑,等.順序同化不同時空分辨率LAI的冬小麥估產對比研究[J].農業機械學報,2015,46(1):240-248.

[2] 邢雅娟,劉東升,王鵬新.遙感信息與作物生長模型的耦合應用研究進展[J].地球科學進展,2009,24(4):444-451.

[3] Wiegand C L, Richardson A J, Kanemasu E T. Leaf area index estimates for wheat from Landsat and their implications for evapotranspiration and crop modeling[J]. Agronomy Journal, 1979,71(2):336-342.

[4] 王鵬新,孫輝濤,王 蕾,等.基于4D-VAR和條件植被溫度指數的冬小麥單產估測[J].農業機械學報,2016,47(3):263-271.

[5] 靳華安,王錦地,柏延臣,等.基于作物生長模型和遙感數據同化的區域玉米產量估算[J].農業工程學報,2012,28(6):162-173.

[6] 黃健熙,馬鴻元,田麗燕,等.基于時間序列LAI和ET同化的冬小麥遙感估產方法比較[J].農業工程學報,2015,31(4):197-203.

[7] Nearing G S, Crow W T, Thorp K R, et al. Assimilating remote sensing observations of leaf area index and soil moisture for wheat yield estimates: An observing system simulation experiment[J]. Water Resources Research, 2012,48(5):W5525.

[8] Huang J, Sedano F, Huang Y, et al. Assimilating a synthetic kalman filter leaf area index series into the WOFOST model to improve regional winter wheat yield estimation[J]. Agricultural and Forest Meteorology, 2016,216:188-202.

[9] He B, Li X, Quan X, et al. Estimating the aboveground dry biomass of grass by assimilation of retrieved LAI into a crop growth model[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2015,8(2):550-561.

[10] Ines A V M, Das N N, Hansen J W, et al. Assimilation of remotely sensed soil moisture and vegetation with a crop simulation model for maize yield prediction[J]. Remote Sensing of Environment, 2013,138(6):149-164.

[11] Huang J, Ma H, Su W, et al. Jointly assimilating MODIS LAI and ET products into the SWAP model for winter wheat yield estimation[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2015,8(8):4060-4071.

[12] Zhang H, Qin S, Ma J, et al. Using residual resampling and sensitivity analysis to improve particle filter data assimilation accuracy[J]. IEEE Geoscience and Remote Sensing Letters, 2013,10(6):1404-1408.

[13] Bindlish R, Jackson T J, Gasiewski A J, et al. Soil moisture mapping and AMSR-E validation using the PSR in SMEX02[J]. Remote Sensing of Environment, 2006,103(2):127-139.

[14] 王鵬新,龔健雅,李小文.條件植被溫度指數及其在干旱監測中的應用[J].武漢大學學報:信息科學版,2001,26(5):412-418.

[15] 張樹譽,李 慧,王鵬新,等.條件植被溫差指數干旱監測方法的研究與應用[J].干旱區資源與環境,2014,28(4):118-122.

[16] 王 蕾,王鵬新,田 苗,等.效率系數和一致性指數及其在干旱預測精度評價中的應用[J].干旱地區農業研究,2016,34(1):229-235.

[17] 林 巧,王鵬新,張樹譽,等.不同時間尺度條件植被溫度指數干旱監測方法的適用性分析[J].干旱區研究,2016,33(1):186-192.

[18] 李 艷,王鵬新,劉峻明,等.基于條件植被溫度指數的冬小麥主要生育時期干旱監測效果評價Ⅲ—干旱對冬小麥產量的影響評估[J].干旱地區農業研究,2014,32(5):218-222.

[19] 王 維,劉翔舸,王鵬新,等.條件植被溫度指數的四維變分與集合卡爾曼同化方法[J].農業工程學報,2011,27(12):184-190.

[20] 王鵬新,孫輝濤,解 毅,等.基于LAI和VTCI及粒子濾波同化算法的冬小麥單產估測[J].農業機械學報,2016,47(4):248-256.

[21] 何慧娟,卓 靜,李紅梅,等.基于MOD16產品的陜西關中地區干旱時空分布特征[J].干旱地區農業研究,2016,34(1):236-241.

[22] 劉驍月,王鵬新,張樹譽,等.基于作物模型模擬年際生物量變化的冬小麥干旱監測研究[J].干旱地區農業研究,2013,31(1):212-218.

[23] 賀 鵬,王鵬新,解 毅,等.基于動態模擬的冬小麥水分脅迫敏感性研究[J].干旱地區農業研究,2016,34(1):213-219.

[24] 孫 威,王鵬新,韓麗娟,等.條件植被溫度指數干旱監測方法的完善[J].農業工程學報,2006,22(2):22-26.

[25] 姜志偉,陳仲新,任建強,等.粒子濾波同化方法在CERES-Wheat作物模型估產中的應用[J].農業工程學報,2012,28(14):138-146.

[26] 蘇 濤,王鵬新,劉翔舸,等.基于熵值組合預測和多時相遙感的春玉米估產[J].農業機械學報,2011,42(1):186-192.

[27] Savitzky A, Golay M J E. Smoothing and differentiation of data by simplified least squares procedures[J]. Analytical Chemistry, 1964,36(8):1627-1639.

[28] 張樹譽,孫 威,王鵬新.條件植被溫度指數干旱監測指標的等級劃分[J].干旱區研究,2010,27(4):600-606.

[29] Weerts A H, El Serafy G Y H. Particle filtering and ensemble kalman filtering for state updating with hydrological conceptual rainfall runoff models[J]. Water Resources Research, 2006,42(9):W9403.