基于4DVAR和EnKF的遙感信息與作物模型冬小麥估產

劉正春 徐占軍 畢如田 王 超 賀 鵬 楊武德

(1.山西農業大學資源環境學院， 太谷 030801； 2.山西農業大學農學院， 太谷 030801)

0 引言

作物生長監測及產量預測對農業管理和國家糧食安全至關重要[1]。根據作物品種特性、氣象條件、土壤條件以及田間管理措施，作物生長模型能夠準確模擬作物在單點尺度的生長發育及作物產量[2]。然而，當作物模型應用到區域尺度時，由于地表環境的異質性，導致模型在參數獲取以及區域化方面存在困難[3]。遙感在大面積獲取地表信息方面具有優勢，能夠有效地解決作物模型區域參數獲取困難的問題[4-5]。因此，利用數據同化技術將遙感信息與作物模型相融合是當前提高區域作物生長模擬和產量估測精度的重要途徑[6]。而同化算法的性能直接影響遙感信息與作物模型同化系統的精度[7]。

目前，應用廣泛且較為成熟的數據同化算法主要有2類：以四維變分(Four-dimensional variation，4DVAR)為主的變分同化算法和以集合卡爾曼濾波(Ensemble Kalman filtering，EnKF)為主的順序同化算法[7-8]。4DVAR和EnKF算法同化遙感信息和作物模型均能有效提高區域作物估產精度[9-10]，但已有研究中對兩種同化算法性能的對比分析較少。

作物模型同化常用的遙感數據為Landsat和MODIS，其時間或空間分辨率較低。MODIS低空間分辨率(0.25～1 km)會出現混合像元現象，導致MODIS-LAI值較實測值偏低，對作物模型同化估產精度的提高有限[3,11]。Landsat影像低的時間分辨率(16 d)，在冬小麥生育期不能提供足夠頻率的觀測數據，對于模型模擬軌跡的調整有限[5,12]。歐洲航天局哥白尼計劃(GMES)中的地球觀測衛星Sentinel多源數據，因其高空間和時間分辨率而成為作物模型同化的理想數據[13-14]。目前，已有研究將Sentinel-2數據與作物模型同化結合進行估產，證明該數據有較好的應用前景[15]。本文基于Sentinel多源遙感數據(Sentinel-1雷達數據反演土壤含水率，Sentinel-2數據反演LAI)，利用4DVAR和EnKF兩種算法同化LAI和土壤含水率，對比分析兩種同化算法的同化性能，并利用同化精度高的LAI和土壤含水率構建估產模型，估測研究區的冬小麥產量。

1 材料與方法

1.1 研究區概況

研究區位于山西省晉南地區，包括聞喜縣、新絳縣和襄汾縣(110°59′33″～111°40′31″ E，35°9′38″～36°3′14″ N)(圖1)。氣候屬暖溫帶大陸性季風氣候，年均降水量在450～600 mm之間，年平均氣溫為12～14℃，全年無霜期160～190 d。該區域作物種植模式主要為冬小麥與夏玉米輪作，冬小麥種植面積108 061 hm2，占耕地面積的73.31%，是研究區的主要農作物[16]。研究區冬小麥一般于10月上旬播種，于6月上旬收獲。

1.2 數據來源

1.2.1田間實測數據

在研究區共設45個實測點，20個紅色樣點的數據用于作物模型同化估產的參數設定及優化，25個黃色樣點用于產量精度驗證(見圖1)。2017年9月冬小麥播種前對20個紅色樣點取土壤樣品進行理化性質分析。土壤剖面分為7層(0～10 cm、10～20 cm、20～50 cm、50～80 cm、80～120 cm、120～160 cm、160～200 cm)，測定了土壤剖面各層的理化性質，包括田間持水量、容重、有機碳含量、全氮質量比、有效鉀質量比、速效磷質量比、銨態氮質量比、pH值和土壤含水率(θ)，作為CERES-Wheat模型的土壤輸入參數。分別于2018年3月19日、2019年3月26日(拔節期)和2018年4月16日、2019年4月19日(抽穗-灌漿期)對全部實測點進行了冬小麥LAI和0～20 cm土壤含水率的測量。在2018年和2019年冬小麥成熟期對全部實測點以1 m×1 m的樣方采樣，將樣本小麥曬干、脫粒并稱量，用干燥后的籽粒計算單產(kg/hm2)。此外，在2017—2019年調查了田間管理數據，主要包括冬小麥品種、主要生育期、灌溉信息和施肥信息。本研究利用2018年的實測和田間管理數據進行CERES-Wheat模型的標定，使用2019年的數據進行研究區冬小麥作物模型同化估產研究。

1.2.2Sentinel多源數據

Sentinel-1和Sentinel-2同屬于Sentinel系列衛星[17-18]，分別為雷達和光學遙感衛星。在冬小麥關鍵生育期(返青期-成熟期)共獲取7景Sentinel-1影像(分別為2019年3月2日、3月14日、3月26日、4月7日、4月19日、5月1日和5月13日)，獲取云量小于10%的Sentinel-2光學影像共7景(分別為2019年3月17日、4月1日、4月16日、5月11日、5月21日、5月26日和6月10日)。

1.3 冬小麥種植區域提取

LI 等[19]利用冬小麥與其他植被不同的物候特征及不同時相的光譜差異采用決策樹法提取冬小麥的種植區域，本文冬小麥種植區域提取參考此算法，選擇2019年3月17日、4月16日和6月10日3景Sentinel-2影像，建立決策樹： 2019年3月17日NDVI大于0.25， 2019年4月16日NDVI大于2019年6月10日NDVI。提取了研究區2019年冬小麥種植區域，總面積為103 254 hm2，與山西省農業局統計數據(108 061 hm2)基本一致，冬小麥種植面積提取精度為96.0%。

1.4 Sentinel-2反演LAI

LAI和NDVI之間有良好的定量關系[20]。首先計算7景Sentinel-2遙感影像的NDVI，并根據田間實測樣點的地理坐標獲取樣點在2019年4月16日遙感影像像元的NDVI，并與2019年4月19日田間實測樣點的LAI建立回歸模型

LAI=8.804 9NDVI-0.986 6

(1)

回歸模型達到極顯著水平(P<0.01)，其決定系數R2為0.52，均方根誤差(Root mean square error, RMSE)為0.995 5 m2/m2。根據公式(1)反演7景Sentinel-2影像的LAI。

1.5 Sentinel-1反演土壤含水率

使用Sentinel-1雷達數據運用水云模型(Water-cloud model)反演冬小麥種植區土壤含水率。該模型假定植被層為一個均質散射體，農作物覆蓋地表的微波后向散射包含農作物直接反射的體散射和經過農作物雙次衰減后的地表散射兩部分[21]。

(2)

回歸模型達到極顯著水平(P<0.01)，其決定系數R2為0.42，RMSE為0.030 5 cm3/cm3。根據公式(2)反演7景Sentinel-1影像的土壤含水率。

1.6 CERES-Wheat模型

CERES-Wheat模型以天為時間步長模擬小麥生長發育、氮碳水平衡過程、產量形成等[2]。

CERES-Wheat模型運行需要4類輸入數據：氣象數據、土壤數據、田間管理數據和作物遺傳參數。其中，氣象數據的最小數據集包括研究樣點的逐日最高氣溫、最低氣溫、降水量和太陽輻射量，從中國氣象數據網(http:∥data.cma.cn/)下載。土壤數據主要包括土壤各層(0～10 cm、10～20 cm、20～50 cm、50～80 cm、80～120 cm、120～160 cm、160～200 cm)的理化參數，如粘粒質量分數、粉粒質量分數、容重、田間持水量、有機碳質量比、全氮質量比、有效鉀質量比、速效磷質量比、pH值等，通過田間實測獲得。田間管理數據為實地調查的冬小麥物候期、灌溉信息和施肥信息。遺傳參數控制著小麥的生長發育進程，直接關系到植株形態的發育與作物產量的高低，其標定后才能用于研究區冬小麥的產量估測[22-23]。遺傳參數的標定一般使用 “試錯法”[24]。

本研究使用“試錯法”，采用2017—2018年冬小麥生育期的氣象數據、土壤數據、田間管理數據和初始作物遺傳參數驅動CERES-Wheat模型，根據研究區20個紅色實測樣點的冬小麥物候期、LAI和單產調整模型的品種遺產參數，直到模型模擬的物候期、LAI、產量與實際調查值差距最小時，確定為本研究區的代表性品種遺傳參數。經過標定的CERES-Wheat模型，模擬收獲日期與實際收獲日期的差值小于4 d，模型模擬LAI與實測LAI的RMSE為1.124 3 m2/m2，模擬單產與實測單產的RMSE為622.23 kg/hm2，由此，標定的CERES-Wheat模型精度可靠，可用于研究區冬小麥生長發育及產量的模擬。

1.7 4DVAR算法

4DVAR算法通過構建代價函數描述模擬狀態變量和觀測值之間的差距，在時間同化窗口T內利用所有外部觀測數據迭代調整模型模擬狀態變量，直至模型的模擬軌跡逼近同化窗口內觀測數據的軌跡，達到優化模型模擬狀態變量的目的[25]。4DVAR的目標函數為

(3)

式中V0——同化窗口初始時刻的狀態變量

Vi——V0代入模型算子M運行到i時刻所得值

B——模擬值的誤差協方差矩陣

Oi——遙感反演值誤差協方差矩陣

n——同化窗口時長

1.8 EnKF算法

EnKF同化算法的核心是在機理模型的動力框架內，把外部觀測數據和模型模擬數據進行誤差加權，計算當前時刻狀態變量的最優估計值并代替模型模擬值，然后運行到下一時刻，直到所有外部觀測值被代入模型，完成模型模擬軌跡的優化[25-26]。EnKF濾波的目標函數為

(4)

(5)

Pk——預報矩陣的協方差矩陣

Rk——觀測矩陣的協方差矩陣

H——觀測算子

Dk——觀測矩陣

M——狀態變換方程

1.9 改進的層次分析法

改進的層次分析法通過構建比較矩陣、傳遞矩陣及擬優一致判斷矩陣確定各因素的權值[27]。本文根據農業經驗知識判定冬小麥主要生育期(返青期、拔節期、抽穗-灌漿期和乳熟期)的LAI或土壤含水率對產量影響的重要程度，建立比較矩陣

(6)

公式(6)是由冬小麥4個主要生育時期LAI對產量的影響程度構建的比較矩陣，0表示第i個生育期沒有第j個生育時期重要；1表示第i個生育期與第j個生育期同樣重要；2表示第i個生育期比第j個生育期重要。

(7)

式中，rmax=max{rj},rmin=min{rj}，k=rmax/rmin。根據C計算傳遞矩陣和擬優一致判斷矩陣，最后通過對擬優一致判斷矩陣的歸一化處理得到冬小麥各生育時期的權重。

2 結果與分析

2.1 LAI同化結果與分析

通過4DVAR和EnKF兩種算法同化CERES-Wheat模型模擬的LAI和Sentinel-2數據反演的LAI，分別得到4DVAR算法同化的LAI(4DVAR-LAI)和EnKF算法同化的LAI (EnKF-LAI)。以灌溉樣點襄汾縣南辛店村、新絳縣澤掌村、聞喜縣蘇村，以及旱作樣點襄汾縣東郭村、聞喜縣柏林村、聞喜縣戶頭村為例，將4DVAR-LAI、EnKF-LAI和CERES-Wheat模型模擬的LAI曲線進行對比(圖2)，4DVAR-LAI和EnKF-LAI保持了模擬LAI的趨勢變化特征，在返青期至抽穗-灌漿期LAI快速增加，至播種后200 d左右LAI達到峰值，之后冬小麥進入乳熟期，葉片逐漸變黃，LAI開始下降。總體上，4DVAR-LAI和EnKF-LAI均保持了模擬LAI趨勢變化特征，且更接近LAI實測值，更符合冬小麥生長的實際變化情況。但是，4DVAR算法對模擬LAI曲線的修正幅度較大，得到的4DVAR-LAI更接近遙感反演LAI，而通過EnKF算法得到的EnKF-LAI對模擬LAI的修正幅度較小；且EnKF-LAI達到最大值的時間與模擬LAI達到最大值的時間同步，而4DVAR-LAI達到最大值的時間受遙感反演LAI值的影響更大，與遙感反演LAI最大值的時間同步。

為檢驗4DVAR-LAI、EnKF-LAI和模擬LAI的精度，采用20個模型模擬樣點的田間實測LAI數據，分別建立4DVAR-LAI、EnKF-LAI、模擬LAI與實測LAI的回歸方程，并計算R2和RMSE。4DVAR-LAI、EnKF-LAI和模擬LAI的R2分別為0.706 8、0.644 6、0.618 6，RMSE分別為1.039 6、1.188 6、1.634 4 m2/m2。與模擬LAI相比，4DVAR-LAI和EnKF-LAI的R2均有所提高，RMSE均降低，說明兩種算法得到的同化LAI更接近實測LAI值，且4DVAR-LAI的精度高于EnKF-LAI(圖3)。

2.2 土壤含水率同化結果與分析

通過4DVAR和EnKF算法同化CERES-Wheat模型模擬的土壤含水率和Sentinel-1數據反演的土壤含水率，分別得到4DVAR算法同化的土壤含水率(4DVAR-θ)和EnKF算法同化的土壤含水率(EnKF-θ)。以灌溉樣點襄汾縣南辛店村、新絳縣澤掌村、聞喜縣蘇村以及旱作樣點襄汾縣東郭村、聞喜縣柏林村、聞喜縣戶頭村為例，對4DVAR-θ、EnKF-θ和CERES-Wheat模型模擬的土壤含水率變化特征進行對比(圖4)。模擬土壤含水率變化特征的每個峰值代表灌溉或降水引起的土壤含水率增加，而后隨著土壤蒸發及作物蒸騰下降，峰值的出現與灌溉日期或降水日期相對應，模擬土壤含水率變化趨勢與實際土壤含水率的變化趨勢相吻合。4DVAR-θ變化特征保持了模型模擬土壤含水率的趨勢變化特征，但更接近Sentinel-1反演土壤含水率，在灌溉樣點的灌溉日期至第1次降水期間(冬小麥種植后170～190 d)，模型模擬的土壤含水率偏高，4DVAR-θ的降低幅度明顯，更接近實測土壤含水率；在旱作樣點，返青期至拔節-抽穗期(冬小麥種植后150～190 d)，降水較少，模型模擬土壤含水率偏低，4DVAR-θ有所提高，更接近實測土壤含水率。灌溉樣點的EnKF-θ在保持模擬土壤含水率變化趨勢下有略微的調整，旱作樣點的EnKF-θ變化特征與模擬土壤含水率基本一致，受遙感反演土壤含水率的修正較小。

為檢驗4DVAR-θ、EnKF-θ和模擬土壤含水率的精度，采用20個模型模擬樣點的田間實測土壤含水率數據，分別建立4DVAR-θ、EnKF-θ、模擬土壤含水率與實測土壤含水率的線性回歸方程，并計算R2和RMSE。4DVAR-θ、EnKF-θ和模擬土壤含水率的R2分別為0.717 1、0.673 4、0.654 7，RMSE分別為0.032 1、0.041 2、0.065 6 cm3/cm3。4DVAR-θ、EnKF-θ的RMSE均低于模擬土壤含水率，說明兩種算法得到的同化土壤含水率更接近實測土壤含水率，且4DVAR-θ的精度高于EnKF-θ(圖5)。

綜上所述，4DVAR和EnKF算法同化的LAI和土壤含水率RMSE均低于CERES-Wheat模型模擬值，表明兩種同化方法都能較好地結合模型模擬值和遙感反演值的優勢，同化LAI和土壤含水率的精度高于模擬值。4DVAR算法同化的LAI和土壤含水率較EnKF更接近田間實測值，因此，本研究選擇4DVAR算法的LAI和土壤含水率同化值估測研究區冬小麥產量。

為利用4DVAR-LAI和4DVAR-θ進行冬小麥區域估產，需將單點尺度的LAI和土壤含水率同化結果擴展到區域尺度。綜合考慮冬小麥的物候期及遙感影像的成像時間等因素，選取2019年3月17日、4月1日、4月16日和5月11日Sentinel-2反演的LAI，2019年3月14日、4月7日、4月19日和5月13日Sentinel-1反演的土壤含水率，作為返青期、拔節期、抽穗-灌漿期和乳熟期的LAI和土壤含水率進行區域同化。將Sentinel多源數據反演的LAI和土壤含水率與20個紅色樣點相應日期的4DVAR-LAI、4DVAR-θ分別進行回歸分析，得到上述物候期冬小麥區域尺度的同化LAI和土壤含水率。

2.3 區域冬小麥產量估測

已有學者研究表明，同時同化LAI和土壤含水率的估產精度高于單獨同化LAI或土壤含水率。因此，本研究采用4DVAR算法同化后的LAI和土壤含水率雙變量進行冬小麥估產，首先，依據改進的層次分析法確定LAI和土壤含水率在冬小麥4個主要生育期(返青期、拔節期、抽穗-灌漿期和乳熟期)對冬小麥產量貢獻的權重分別為0.055 0、0.265 0、0.566 0、0.114 0和0.055 5、0.565 5、0.260 5、0.118 5，由此計算冬小麥主要生育期加權4DVAR-LAI、4DVAR-θ；之后與圖1中20個紅色實測點的產量數據進行回歸分析建立估產模型Y=1 169.74LAI+8 205.38θ-396.42，模型顯著性檢驗P<0.001，R2=0.886 4。采用圖1中的25個黃色實測點數據進行產量精度驗證，同化4DVAR-LAI、4DVAR-θ雙變量建立的估產模型精度RMSE為449.77 kg/hm2，平均相對誤差(Mean relative error，MRE)為7.85%，CERES-Wheat模型直接模擬產量的RMSE為641.55 kg/hm2，MRE為10.23%。4DVAR-LAI+4DVAR-θ建立的估產模型的RMSE和MRE較模型模擬產量的估產誤差小，表明同化遙感信息與作物模型能有效提高估產精度。

采用由4DVAR-LAI+4DVAR-θ建立的估產模型進行區域冬小麥估產，得到研究區2019年冬小麥產量如圖6a所示。結合研究區數字高程模型(Digital elevation model, DEM)圖(圖6b)及山西省耕地地力評價數據庫中灌溉能力數據(圖6c)可以發現，高程小于500 m且滿足灌溉條件的區域，冬小麥產量較高(大于5 500 kg/hm2)；在聞喜縣東、西兩側高程大于780 m且不滿足灌溉條件的區域，冬小麥產量較低(小于3 500 kg/hm2)；高程處于500～780 m且基本滿足灌溉條件的區域，冬小麥產量在3 500～5 500 kg/hm2之間。

3 討論

3.1 4DVAR和EnKF 算法在同化遙感信息和作物模型中性能對比分析

本研究通過4DVAR和EnKF算法同化Sentinel多源數據和CERES-Wheat模型，4DVAR算法同化的LAI和土壤含水率精度高于EnKF算法，且4DVAR算法同化的LAI的冬小麥物候期更接近冬小麥實際生長物候期，因此4DVAR算法在Sentinel多源數據和CERES-Wheat模型同化系統中表現更優，同化效果更好。這可能是由于所采用的高時空分辨率的遙感信息Sentinel數據反演的LAI和土壤含水率精度較高，4DVAR算法通過代價函數迭代調整模型模擬狀態變量，將模型的模擬軌跡逼近同化窗口內所有觀測數據的軌跡，達到優化模型模擬狀態變量的目的[28-29]，所以4DVAR算法同化的狀態變量與遙感觀測值更接近，當遙感觀測值精度較高時，其同化性能更優。解毅等[8]利用4DVAR和EnKF算法同化Landsat數據和CERES-Wheat模型，結果顯示EnKF同化的LAI精度高于4DVAR算法。這可能是由于Landsat較低的時空分辨率導致Landsat反演的LAI準確度較差，4DVAR同化結果逼近Landsat反演的LAI，以致4DVAR優化效果較差；而EnKF集成濾波器的狀態相關不確定性信息和集成構建能力能夠克服遙感觀測不確定性的影響[28-29]，所以Landsat-LAI與CERES-Wheat模型進行耦合時，采用EnKF算法得到的同化LAI較4DVAR有更高的精度。由此，本研究認為不同的同化算法應用于不同的遙感信息與作物模型同化系統中的表現不同，這與HUANG等[7]在研究作物模型同化方法中發現同化算法的表現與遙感影像的分辨率高度相關的結論相一致。在以后的研究中，將采用不同精度的遙感數據與作物模型進行同化估產研究，比較4DVAR和EnKF在不同分辨率遙感影像中的同化性能。

3.2 4DVAR算法同化Sentinel數據和CERES-Wheat模型的估產精度分析

本研究表明4DVAR算法同化Sentinel數據和CERES-Wheat模型的區域冬小麥產量估測精度較高(4DVAR同化后的LAI和土壤含水率建立的估產模型的RMSE和MRE較CERES-Wheat模型模擬的RMSE和MRE分別降低了191.78 kg/hm2、2.38%)。冬小麥的產量空間分布呈現出隨海拔高度的增加、灌溉條件的不足而遞減的趨勢，這與冬小麥實際產量空間分布相一致(圖6)。證明4DVAR-LAI+4DVAR-θ雙變量建立的估產模型估測的區域冬小麥產量空間分布合理。

本文采用高時空分辨率的遙感信息(Sentinel數據)，并采用4DVAR和EnKF兩種算法進行同化性能對比分析，為本研究區的Sentinel數據與CERES-Wheat模型同化系統選擇了精度更高的4DVAR算法，從遙感數據源與同化算法兩方面提高了遙感信息與作物模型同化系統的精度，得到高精度的同化LAI和土壤含水率，并用來構建估產模型，有效提高了研究區冬小麥估產精度。

4 結論

(1)4DVAR和EnKF兩種算法同化的LAI和土壤含水率θ的R2均高于模型模擬值，RMSE均低于模擬值。兩種同化算法都能較好地結合遙感觀測和模型模擬的優勢，提高了同化LAI和土壤含水率的精度。

(2)4DVAR-LAI和4DVAR-θ的RMSE分別比EnKF-LAI和EnKF-θ低0.149 0 m2/m2、0.009 1 cm3/cm3，4DVAR同化結果更符合實際LAI和土壤含水率實際變化情況，且根據遙感實際監測值，4DVAR算法同化的LAI更能精確識別冬小麥的物候期，與實際冬小麥生長發育的物候期更相符。這是由于4DVAR算法通過構建代價函數調整模型模擬狀態變量，將模型的模擬軌跡逼近觀測數據的軌跡，當采用的遙感信息精度高時，4DVAR算法的同化表現更優。

(3)由4DVAR同化后LAI和土壤含水率雙變量建立的估產模型精度(RMSE為449.77 kg/hm2，MRE為7.85%)高于CERES-Wheat模型模擬的單產精度(RMSE為641.55 kg/hm2，MRE為10.23%)。本研究從遙感數據源與同化算法兩方面提高了遙感信息與作物模型同化系統的精度，得到高精度的同化LAI和土壤含水率，并用來構建估產模型，有效提高了研究區冬小麥的估產精度。