水氮數(shù)據(jù)聯(lián)合驅(qū)動(dòng)的冬小麥數(shù)據(jù)同化研究

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

（武漢大學(xué)水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢430072）

0 引 言

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

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

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

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

1 材料與方法

1.1 SWAP模型

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

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

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

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

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

模型觀測(cè)可以表示為：

因此有：

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

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

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

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

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

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

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

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

1.3.2 不確定性參數(shù)

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

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

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

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

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

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

2 結(jié)果與分析

2.1 更新LNA的重要性

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

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

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

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

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

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

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

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

2.3 觀測(cè)誤差的影響

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

2.4 觀測(cè)頻率的影響

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

3 結(jié) 論

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