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

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

0 引 言

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

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

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

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

1 材料與方法

1.1 SWAP模型

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

1.2 數(shù)據(jù)同化方法

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

進行同化的狀態(tài)向量通常由模型參數(shù)、狀態(tài)和觀測共同組成：

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

模型觀測可以表示為：

因此有：

式中：H為觀測算子，將狀態(tài)向量矩陣轉(zhuǎn)換為觀測值相對應(yīng)的模型預(yù)報值。

采用EnKF對模型狀態(tài)向量進行更新：

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

式中：為模型預(yù)測狀態(tài)向量的協(xié)方差矩陣；為狀態(tài)向量與觀測向量的互協(xié)方差；為觀測向量的自協(xié)方差；為觀測誤差。

1.3 虛擬觀測模擬系統(tǒng)

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

1.3.1 SWAP-WOFOST 模型輸入數(shù)據(jù)

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

1.3.2 不確定性參數(shù)

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

式中：EEi為第i個模型參數(shù)的基效應(yīng)數(shù)值；x1，x2，…，xm為m個模型參數(shù)的取值；Δ為參數(shù)增量；f為模型狀態(tài)輸出。

通過計算基效應(yīng)EEi絕對值的均值來得到敏感性指標μ*，μ*越大，表示參數(shù)越敏感。本研究采用的樣本數(shù)量為1 000，模型參數(shù)為50 個，共需計算51 000 次運算，計算的敏感性結(jié)果見圖1。由圖1 可知，模型的50 個參數(shù)并不都是敏感參數(shù)，只有部分參數(shù)在模型運算過程中起到重要的影響。其中土壤水參數(shù)ALFA、NPAR、OSAT，土壤氮參數(shù)TCSF_N 和作物參數(shù)EFF、CVL、FS_090 對于模型輸出的產(chǎn)量、LAI和LNA都有影響，尤其是TCSF_N，敏感性指標數(shù)值很高；有些參數(shù)只對產(chǎn)量敏感，比如CVO 和FO_110。結(jié)合HU 等的研究，本文中選取了ALFA、NPAR、TSUMEA、EFF、CVL、CVO、TCSF_N等關(guān)鍵參數(shù)[16]。其中新加入的TCSF_N 為蒸騰濃度流系數(shù)，與作物根系吸收土壤氮的速率有關(guān)，可以反映作物吸收氮的能力，從而影響作物生長。假設(shè)所有不確定性參數(shù)服從對數(shù)正態(tài)分布，初始均值和方差為：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 算例設(shè)計

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

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

表1 數(shù)據(jù)同化算例設(shè)計Tab.1 Case design of data assimilation

2 結(jié)果與分析

2.1 更新LNA的重要性

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

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

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

2.2 SSPE和USO數(shù)據(jù)同化策略對比研究

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

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

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

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

2.3 觀測誤差的影響

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

2.4 觀測頻率的影響

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

3 結(jié) 論

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