水氮數據聯合驅動的冬小麥數據同化研究

李金敏，韓景曄，余丹陽，史良勝

（武漢大學水資源與水電工程科學國家重點實驗室，武漢430072）

0 引 言

中國是世界上最大的小麥生產國和消費國，小麥也是中國第3大農作物，僅次于水稻和玉米[1,2]。因此對冬小麥進行準確的生長模擬和產量預測是十分重要的。作物模型從系統科學的觀點出發，運用數學物理方法和計算機技術，對作物生長過程中的光合、呼吸、蒸騰等重要生理生態過程進行定量描述和預測[3]，是作物生長發育及其對環境影響評估的重要工具[4]。

在過去的半個世紀中，研究者們開發出很多作物模型來模擬冬小麥的生長過程，比如AquaCrop[5,6]，DSSAT[7,8]，WOFOST[9]，SWAP[10]和STICS[11]等。然而，作物模型的精度通常會受到模型結構、參數和輸入等因素引起的不確定性的影響，導致模擬結果和實際情況有偏差[11]。

數據同化是一種集成觀測和模型這2種基本科學研究手段的重要方法，通過融合多源觀測和帶有不確定性的模型，從而改善對作物生長過程的模擬[13,14]。目前，數據同化方法在作物生長模擬和產量估計研究中受到廣泛關注[15-19]。Wang 等結合遙感數據反演的LAI和葉片氮累積量（LNA），利用粒子群優化算法對作物生長模型RiceGrow 的播種日期、播種密度和施氮量等參數進行初始化參數校準，提高了水稻產量的預測精度[20]。Huang 等采用集合均方根濾波（EnSRF）方法，將遙感反演的LAI 和LNA作為觀測加入到WheatGrow 模型數據同化系統中，成功地估計了冬小麥LAI和LNA，提高了產量的估計精度[21]。然而，這些作物模型數據同化系統大多關注作物葉面積指數（LAI），對土壤和作物水氮狀態的數據同化研究相對較少。實際上，水和氮是作物生長的2大關鍵因素，對于土壤和作物水氮狀態的準確估計是十分重要的[22-26]。因此，本研究中嘗試利用在數據同化系統中同時更新土壤含水量（SM）和作物葉片氮累積量（LNA），探究水氮數據聯合驅動更新對作物生長模擬的影響。

本研究中的冬小麥生長模型為SWAP-WOFOST 模型，該模型采用WOFOST 的作物生長模塊，兼顧土壤-作物水氮平衡，已經成功地應用于作物的生長模擬和產量估計研究中[17，18，27]。基 于SWAP-WOFOST 模 型 構 建 虛 擬 觀 測 模 擬（OSSE）系統，采用集合卡爾曼濾波（EnKF）方法進行順序性數據同化，圍繞SWAP-WOFOST-EnKF 系統開展系統研究：土壤含水量和葉片氮累積量的聯合更新對模擬結果的影響、觀測數據的價值分析和不同數據同化更新策略對模擬結果的影響。

1 材料與方法

1.1 SWAP模型

SWAP（Soil-Water-Atmosphere-Plant）模型是荷蘭瓦赫寧恩大學開發的一個垂向一維農業水文模型，它可以用來模擬非飽和區域內土壤水、溶質和熱量的運移與作物生長之間的關系。SWAP 采用Richards 方程模擬土壤水運動，通過對流、彌散和擴散方程來描述溶質運移轉化等過程。SWAP模型中的作物生長部分采用WOFOST 模塊，通過量化光合作用、呼吸消耗等過程，并考慮水肥脅迫等環境影響，形成的干物質在不同的器官間進行分配，而這一切都受到作物生育期（DVS）的影響，DVS在種子萌發時為0，在開花時取1，最終成熟時DVS=2。通過耦合土壤氮模塊，計算出土壤中氨態氮、硝態氮的動態平衡。土壤中氮素通過作物根系吸收進入作物體內，并根據作物氮需求進行氮素分配，通過計算氮素脅迫因子NNI來影響作物光合、葉生長、氮素分配等過程。

1.2 數據同化方法

集合卡爾曼濾波（EnKF）是一種順序性數據同化算法，可以通過不斷獲取觀測數據對模型狀態和參數進行持續更新。由于其具有較低的計算成本，并且具有一定程度非線性問題適應能力，在農業水文領域有著廣泛的應用[16,28-30]。

進行同化的狀態向量通常由模型參數、狀態和觀測共同組成：

式中：yt為t時刻的狀態向量；為模型參數，是模型狀態變量，如土壤含水量和作物葉面積指數等；是觀測值，上標T表示矩陣轉置。

模型觀測可以表示為：

因此有：

式中：H為觀測算子，將狀態向量矩陣轉換為觀測值相對應的模型預報值。

采用EnKF對模型狀態向量進行更新：

式中：j為樣本序號為數據同化更新后的狀態向量為模型預測的狀態向量；d t,j為擾動后的觀測值；Kt是卡爾曼增益，記為如下公式。

式中：為模型預測狀態向量的協方差矩陣；為狀態向量與觀測向量的互協方差；為觀測向量的自協方差；為觀測誤差。

1.3 虛擬觀測模擬系統

本研究中，采用虛擬觀測模擬（Observing Simulation System Experiment， OSSE）系統完成相關研究內容。OSSE 系統共包括N+1次樣本模擬，其中的1個樣本模型結果作為參照值，余下的N個樣本由帶有偏差的先驗參數抽樣后模擬得來，采用參照值來評價模擬效果。在本研究中，樣本量N=50，與de Wit等和Curnel等的研究保持一致[31,32]。

1.3.1 SWAP-WOFOST 模型輸入數據

氣象數據來源于武漢大學灌溉排水試驗場的氣象站，包括氣溫（最高和最低溫度）、太陽輻射、降雨、空氣水汽壓和風速，數據頻率為1 d，模擬生長時間從2017年的11月18日至2018年的5月10日，共174 d。模型上邊界為大氣邊界，下邊界采用隔水邊界。施肥情況為在作物播種當天施加氮素100 kg/hm2。

1.3.2 不確定性參數

在進行數據同化模擬時，需要通過擾動不確定性參數來生成樣本，因此敏感參數的選擇是非常重要的一步。本研究中采用Morris篩選法[33]，主要分析模型輸出的產量（Yield）、LAI和葉片氮累積量（LNA）對50個參數的敏感性。Morris方法首先計算參數的基效應：

式中：EEi為第i個模型參數的基效應數值；x1，x2，…，xm為m個模型參數的取值；Δ為參數增量；f為模型狀態輸出。

通過計算基效應EEi絕對值的均值來得到敏感性指標μ*，μ*越大，表示參數越敏感。本研究采用的樣本數量為1 000，模型參數為50 個，共需計算51 000 次運算，計算的敏感性結果見圖1。由圖1 可知，模型的50 個參數并不都是敏感參數，只有部分參數在模型運算過程中起到重要的影響。其中土壤水參數ALFA、NPAR、OSAT，土壤氮參數TCSF_N 和作物參數EFF、CVL、FS_090 對于模型輸出的產量、LAI和LNA都有影響，尤其是TCSF_N，敏感性指標數值很高；有些參數只對產量敏感，比如CVO 和FO_110。結合HU 等的研究，本文中選取了ALFA、NPAR、TSUMEA、EFF、CVL、CVO、TCSF_N等關鍵參數[16]。其中新加入的TCSF_N 為蒸騰濃度流系數，與作物根系吸收土壤氮的速率有關，可以反映作物吸收氮的能力，從而影響作物生長。假設所有不確定性參數服從對數正態分布，初始均值和方差為：ln（ALFA）～（-3.0，0.10），ln（NPAR）～（0.53，0.05），ln（TSUMEA）～（-6.91，0.03），ln（EFF）～（-0.60，0.08），ln（CVL）～（-0.29，0.03），ln（CVO）～（-0.22，0.03），ln（TCSF_N）～（-0.36，0.10）。

1.3.3 算例設計

本研究基于SWAP-WOFOST 模型構建冬小麥作物生長的數據同化系統，利用多種觀測聯合估計作物水氮狀態。選擇合適的狀態變量是進行數據同化模擬的關鍵。傳統研究大都將土壤含水量、LAI作為狀態變量，本研究在此基礎上加入了葉片氮累積量（LNA）這一狀態變量，進行水氮狀態的聯合更新。其中觀測時間為作物生長的第10到第160 d，觀測頻率為10 d一次，全生育期一共16次觀測。土壤水觀測為0～10 cm的平均土壤含水量（SM），觀測誤差為±0.04 m3/m3，LAI的觀測誤差為±0.4 m2/m2，LNA的觀測誤差為±4 kg/hm2。

算例設計見表1。其中，C1 為利用傳統的土壤含水量和LAI觀測進行數據同化，狀態變量為負壓（H）、土壤含水量（SM）、葉干重（LV）、葉齡和葉片氮累積量（LNA）。在C1 的基礎上，C2 的觀測中加入葉片氮累積量（LNA），試圖探索加入LNA觀測對數據同化系統的影響。C2 與C3 的對比是為了研究2種同化策略對土壤水分剖面和作物生長模擬的價值，其中SSPE 指的是同時更新作物模型的狀態和參數，而USO 只更新模型狀態。在C4中，觀測誤差由固定值改變為觀測值的10%，通過和C2 結果對比，研究觀測誤差對數據同化系統模擬的影響。在C5 中，通過將觀測頻率從10 d 增加到20 d，從而探究不同觀測頻率對作物生長模擬的影響。

表1 數據同化算例設計Tab.1 Case design of data assimilation

2 結果與分析

2.1 更新LNA的重要性

圖2展示了算例C1和C2模擬的LAI、LNA、作物產量以及5、20 和40 cm 處土壤含水量隨時間的變化，其中橫坐標為作物播種后天數，“Reference”表示OSSE 系統模擬的真值，“Open-loop”表示不進行數據同化的確定性模型結果。由圖2（a）可以發現，通過在算例C1 中加入LAI觀測可以實現較好的LAI模擬，但是仍有一定改善空間；算例C2 中加入LNA觀測后，作物LAI的模擬更加接近真實值，表明加入作物葉片氮信息有助于LAI的更新，但是總體來看，C1和C2對于LAI的模擬結果相差不大。由圖2（b）中可知，算例C1 中由于僅加入LAI的觀測，不能很好地模擬出作物葉片氮累積量的信息，尤其是在作物生長的中間階段（第100～130 d）。相比于C1，C2可以更好地模擬出作物生長后期的LNA變化，從而可以更加精準地模擬作物的產量。然而，無論是C1 還是C2，對于作物早期的狀態都沒有能夠很好地重現，這主要是由于數據同化過程中給定的觀測誤差較大，尤其在早期作物LAI和LNA都比較小的情況下，給出較大的觀測誤差會影響狀態向量的更新，因此對于作物早期狀態的估計都存在一定的改善空間。由圖2（c）可以發現，C1 中對于作物產量的估計效果很差，而加入LNA觀測的C2算例能夠較好地重現產量的積累過程。

圖2（d）～圖2（f）展示了算例C1 和C2 模擬的5、20 和40 cm 處的土壤含水量隨時間的變化，可以發現C1 和C2 都能較好地重現出土壤含水量的剖面分布。整體來看，C1 對于土壤含水量的模擬效果優于C2，可能是由于加入LNA觀測后，引入了新的觀測誤差，導致C2 中對于土壤含水量的模擬精度略低于C1。

對比算例C1 和C2 對于4 種狀態變量（SM、LAI、LNA和作物產量）的更新結果可以發現，僅加入LAI這一個作物狀態的觀測可以改善模型對LAI的模擬，但是對于作物產量的模擬效果不佳，主要是由于沒有很好地更新模型的氮素信息；在模型能夠準確模擬水分狀態的情況下，通過加入LNA的觀測，可以有效提升作物對氮素狀態的模擬，從而提升作物LAI和產量的估計精度。

2.2 SSPE和USO數據同化策略對比研究

圖3 展示了不同數據同化策略下算例C2 和C3 模擬的LAI、LNA、作物產量以及5、20 和40 cm 處土壤含水量隨時間的變化。相比于算例C1，可以發現C2 和C3 都能很好地重現作物LAI和LNA隨時間的變化，主要是因為C2 和C3 中都加入了LNA的觀測，能夠很好捕捉作物氮素信息。然而，相比于同時更新模型狀態和參數的C2，僅更新狀態的C3在作物產量的估計上效果不佳。

圖3（d）～圖3（f）展示了算例C2 和C3 模擬的5、20 和40 cm 處的土壤含水量隨時間的變化，可以發現不更新參數的C3 算例不能很好地模擬5 cm 處和20 cm 處的土壤含水量，對于土壤水分狀態的模擬存在偏差。

對比算例C2 和C3 對于4 種狀態變量（土壤含水量、LAI、LNA和作物產量）的更新結果可以發現，只更新模型狀態的C3 算例盡管能夠很好地模擬作物狀態，但是由于缺乏對土壤水分狀態的正確估計，導致其產量估計依舊出現很大的誤差；而同時更新狀態和參數的SSPE 方法可以很好重現土壤—作物水分和氮素的變化，因此能夠更好地實現作物的產量估計。圖4展示了在SSPE方法中土壤水參數NPAR、作物參數EFF和氮素參數TCSF_N 的變化過程，可以發現3 類參數在數據同化模擬過程中都會往參照值靠近，不確定性逐漸減小，即SSPE方法可以同時減小模型參數和過程的不確定性，從而獲得更加準確的模擬效果。

對比C1、C2和C3這3種算例可以發現，C1很好地重現了土壤含水量的變化，但是不能準確模擬作物氮素狀態；而C3很好地重現了作物葉片氮素累積的變化，但是不能準確模擬土壤的水分狀態，因此這2個算例對于最終作物產量的估計都存在一定偏差。C2 由于同時更新了土壤—作物的水氮狀態，有效提升了產量估計的精度。因此，準確的土壤—作物水氮狀態的模擬對于作物生長的模擬和產量估計至關重要。

2.3 觀測誤差的影響

圖5 展示了在不同觀測誤差下，算例C2 和C4 模擬的LAI、LNA和作物產量隨時間的變化。C2 采用固定的觀測誤差值，SM為±0.04 m3/m3，LAI為±0.4 m2/m2，LNA為±4 kg/hm2，而C4采用相對的觀測誤差，即任一時刻的觀測誤差為對應觀測值的10%。由圖5可知，相比于C2，C4可以獲得更加精準的LAI和LNA模擬效果，尤其在作物生長的早期，由于C2 采用固定的觀測誤差，而作物生長早期LAI和LNA值都較小，觀測誤差就相對較大，因此同化效果不佳；C4 中采用相對觀測誤差，在整個模擬過程中都能較好地重現作物生長的狀態。盡管C4能夠完成更好的作物LAI和LNA模擬，C2 的模擬效果也不算差，而且C2估計的最終產量比C2更接近參照值。考慮到觀測誤差通常取決于觀測手段和觀測儀器，采用固定觀測誤差可能更符合實際應用情景。

2.4 觀測頻率的影響

圖6 展示了在不同觀測頻率下，算例C2 和C5 模擬的LAI、LNA和作物產量隨時間的變化。其中C2 采用10 d 一次的觀測頻率，而C5 為20 d 一次。由圖6 可知，減小觀測頻率的C5 在LAI和LNA模擬過程中波動比C2 要大一些，但是在2 種算例中，最終產量的估計卻是一致的，都非常接近參照值。這表明，適當地減小觀測頻率并不會降低模型的預測能力，數據同化系統在不同觀測頻率下具有一定的穩健性。對于冬小麥而言，間隔20 d 一次的觀測幾乎可以覆蓋所有的生育階段，因此從權衡模型表現和觀測成本的角度來看，每個生育期1次觀測即可滿足本研究中的數據同化模擬要求。

3 結 論

傳統的數據同化系統只有土壤水和LAI的觀測，無法有效估計作物氮素狀態，通過加入葉片氮累積量的觀測，可以很好地重現作物葉片氮素變化過程。相比于只更新狀態的USO方法，同時更新模型狀態和參數的SSPE 方法在更新土壤含水量方面具有更加重要的價值。水分和氮素狀態的錯誤估計都會引起作物產量模擬存在偏差，因此水氮數據的聯合驅動對于作物生長模擬和產量估計是十分必要的。此外，本研究中的數據同化系統在不同觀測誤差和觀測頻率下都具有良好的穩健性。