楊 華，姜海燕 ，趙空暖，錢崢遠(yuǎn)

（南京農(nóng)業(yè)大學(xué)人工智能學(xué)院，南京 210095）

0 引言

作物生長(zhǎng)模擬模型（以下簡(jiǎn)稱“作物模型”）是作物精確栽培技術(shù)研究的重要工具，可以動(dòng)態(tài)模擬不同品種在生長(zhǎng)發(fā)育生理過程中與環(huán)境變量間的定量關(guān)系[1]。由于建模人員對(duì)作物生理生態(tài)過程認(rèn)識(shí)不足，無(wú)法考慮全部的限制因子，導(dǎo)致模型存在一定的不確定性[2]，會(huì)使使用者對(duì)模型預(yù)測(cè)結(jié)果不信任并影響應(yīng)用效果[3]。這一因素限制了作物模型的大規(guī)模推廣，因此在應(yīng)用時(shí)有必要進(jìn)行不確定性優(yōu)化。按照來(lái)源，作物模型的不確定性可分為結(jié)構(gòu)不確定性、參數(shù)取值不確定性、數(shù)據(jù)不確定性以及主觀的不確定性[4-5]。模型不確定性優(yōu)化方法亦是圍繞上述幾個(gè)方面展開，主要分為參數(shù)不確定性優(yōu)化、模型結(jié)構(gòu)不確定性優(yōu)化及綜合的不確定性優(yōu)化[6]。其中參數(shù)不確定性優(yōu)化是作物生長(zhǎng)模型應(yīng)用的核心環(huán)節(jié)。作物模型參數(shù)主要是指品種參數(shù)，雖然模型自身提供了參數(shù)的參考值，但部分參數(shù)隨生態(tài)點(diǎn)地理位置及作物品種等變化而變化[7]，因此，在應(yīng)用模型時(shí)先要對(duì)這部分參數(shù)進(jìn)行優(yōu)化并驗(yàn)證，提升其在某一地區(qū)的適用性。

校正作物模型品種參數(shù)的常用方法是基于統(tǒng)計(jì)理論的，主要包括廣義似然不確定性估計(jì)（generalized likelihood uncertainty estimation，GLUE）方法和馬爾科夫鏈蒙特卡羅（Markov chain Monte Carlo，MCMC）方法[8-9]。GLUE 方法存在計(jì)算時(shí)消耗資源過大，運(yùn)行周期較長(zhǎng)的問題，而且只考慮了參數(shù)的先驗(yàn)分布，未將樣本信息先驗(yàn)分布以似然函數(shù)（likelihood function，LF）形式化，這可能低估了模型的不確定性，導(dǎo)致參數(shù)校正結(jié)果及不確定度（uncertainty ratio，UR）分析時(shí)造成偏差[10]。而MCMC 方法充分利用了樣本信息先驗(yàn)分布，并將其用LF 表示，結(jié)合參數(shù)先驗(yàn)信息，以概率密度核函數(shù)的形式來(lái)表示模型參數(shù)分布[11]，并能夠量化參數(shù)不確定性，因此成為不確定條件下作物模型參數(shù)校正的主流方法。黃健熙等[12-13]利用MCMC 對(duì)WOFOST、ORYZA 系列等模型進(jìn)行參數(shù)標(biāo)定，均取得良好效果。WALLACH 等[14-17]采用MCMC 進(jìn)行模型不確定性量化和分析，為模型本地化應(yīng)用提供參考，TAN 等[18]初步比較了GLUE 和MCMC 方法對(duì)參數(shù)校正結(jié)果的影響，發(fā)現(xiàn)MCMC 方法一定程度上優(yōu)于GLUE。

利用MCMC 方法對(duì)模型參數(shù)校正的關(guān)鍵是LF 形式設(shè)計(jì)，LF 代表了模型殘差的分布特點(diǎn)，大部分作物模型研究中都是將LF 設(shè)為高斯正態(tài)似然函數(shù)（gaussian likelihood function，GLF）形式[19]，GLF 要求模型殘差的方差恒定，然而作物模型殘差具有異方差性，即模型在可控制變量條件下具有變化的方差，這種變化的方差是由觀測(cè)數(shù)據(jù)的不確定性和模型本身的復(fù)雜性造成的[20]。有研究表明，在農(nóng)業(yè)觀測(cè)數(shù)據(jù)中，由于測(cè)量手段、測(cè)量人員主觀性和環(huán)境因素影響，關(guān)鍵物候期和產(chǎn)量等測(cè)量值隨年份變化而波動(dòng)較大，會(huì)導(dǎo)致模型殘差平穩(wěn)性較差和離散程度較大，從而導(dǎo)致異方差性[21]。這會(huì)給參數(shù)校正的結(jié)果帶來(lái)偏差，影響作物模型的應(yīng)用。

本研究以RiceGrow 和Oryza2000 水稻物侯期模型為研究對(duì)象，以中國(guó)南方地區(qū)早、中、晚3 種不同熟性的水稻品種栽培試驗(yàn)數(shù)據(jù)為基礎(chǔ)，通過引入變異系數(shù)（coefficient of variation，CV）變換的高斯似然函（GLF with CV transformation，GLF-CV）和引入BC（Box-Cox）變換的高斯似然函數(shù)（GLF with BC transformation，GLFBC）對(duì)觀測(cè)數(shù)據(jù)和模型結(jié)構(gòu)造成的異方差進(jìn)行特征描述，并以參數(shù)后驗(yàn)分布UR 和模型預(yù)測(cè)UR 為評(píng)價(jià)標(biāo)準(zhǔn)進(jìn)行不同LF 下參數(shù)校正結(jié)果比較，以期為利用MCMC 方法進(jìn)行水稻生長(zhǎng)模型參數(shù)校正的LF 選擇提供參考，同時(shí)也為作物生長(zhǎng)模型本地化應(yīng)用提供指導(dǎo)。

1 材料與方法

1.1 站點(diǎn)數(shù)據(jù)

研究數(shù)據(jù)集選取中國(guó)廣東省肇慶市高要區(qū)2004—2009 年、江蘇省泰州市興化市2001—2004 年、安徽省六安市1991—2004 年3 個(gè)南方地區(qū)水稻種植生態(tài)點(diǎn)的早熟（雪花粘）、中熟（武育粳3 號(hào)）、晚熟（汕優(yōu)63 號(hào)）水稻品種各年份的田間試驗(yàn)資料，品種和栽培信息見表1。其中，雪花粘的播種期范圍為3 月7—13 日，成熟期范圍為7 月8—15 日左右；武育粳3 號(hào)播種期范圍為5 月4—13 日左右，成熟期范圍為10 月8—15 日左右；汕優(yōu)63 號(hào)播種期范圍為4 月18—25 日左右，成熟期范圍為9 月3—15 日左右。

表1 不同熟性水稻品種種植地點(diǎn)和年份Table 1 Planting place and year of rice varieties with different maturity

氣象數(shù)據(jù)來(lái)源：高要、興化、六安3 個(gè)地點(diǎn)各年份氣象數(shù)據(jù)均來(lái)源于國(guó)家氣象局氣象中心，基礎(chǔ)氣象數(shù)據(jù)包括日最高溫度（℃）、日最低溫度（℃），日照時(shí)數(shù)（h/d）和降雨量（mm）。

1.2 模型品種參數(shù)

研究采用2 個(gè)非線性程度不同的物侯期模擬模型，分別是由南京農(nóng)業(yè)大學(xué)國(guó)家信息農(nóng)業(yè)工程技術(shù)中心研制的RiceGrow 物候期模型[22]和由國(guó)際水稻研究所及荷蘭瓦格寧根大學(xué)聯(lián)合開發(fā)的Oryza2000 物候期模型[23-24]。各模型敏感性品種參數(shù)描述及范圍如表2 所示。

1.3 模型參數(shù)校正方法

1.3.1 貝葉斯統(tǒng)計(jì)推斷

對(duì)于水稻物侯期模型的驅(qū)動(dòng)數(shù)據(jù)、模型品種參數(shù)和模型輸出，公式描述為

其中yi代表模型第i個(gè)試驗(yàn)輸出關(guān)鍵物候期（如拔節(jié)期、成熟期等）儒歷天數(shù)（confucian calendar days，CCD）。為第i個(gè)驅(qū)動(dòng)數(shù)據(jù)，如氣象數(shù)據(jù)、土壤數(shù)據(jù)和田間管理數(shù)據(jù)等，θ為模型品種參數(shù)，εi為模型殘差，根據(jù)貝葉斯公式，模型參數(shù)的后驗(yàn)分布為

式中p(yi|θ,xi)為模型殘差εi先驗(yàn)分布。本研究中用實(shí)際關(guān)鍵物候期與模型模擬值CCD 之間差值的LF 表示，p(θ)為模型參數(shù)θ的先驗(yàn)分布。

1.3.2 MCMC 參數(shù)校正框架

在作物模型領(lǐng)域，MCMC 方法是基于貝葉斯推理的參數(shù)后驗(yàn)分布估計(jì)中最有效的方法。該方法進(jìn)行模型參數(shù)估算的基本思想是產(chǎn)生1 個(gè)馬爾科夫鏈，以目標(biāo)分布為平穩(wěn)分布，目標(biāo)分布一般為p(yi|θ,xi)，根據(jù)馬爾科夫理論，1 個(gè)馬爾科夫鏈從任意初值出發(fā)，都會(huì)收斂到平穩(wěn)分布。MCMC 方法的目的是在構(gòu)造1 個(gè)馬爾可夫鏈的基礎(chǔ)上生成參數(shù)集的樣本，通過參數(shù)集的平穩(wěn)分布即可求出每個(gè)參數(shù)的后驗(yàn)分布[25]?？蚣芤妶D1。首先設(shè)定品種參數(shù)先驗(yàn)分布、模型約束目標(biāo)變量殘差的LF 和MCMC采樣算法，從參數(shù)先驗(yàn)分布中選擇參數(shù)值進(jìn)行初始化，初始化值確定依據(jù)具體的MCMC 采樣算法而定。在田間管理數(shù)據(jù)、土壤數(shù)據(jù)和氣象數(shù)據(jù)等驅(qū)動(dòng)下，運(yùn)行生長(zhǎng)模型得到約束目標(biāo)變量值，將此值與LF 進(jìn)行符合程度比較，并依據(jù)符合程度構(gòu)造參數(shù)的馬爾科夫鏈，經(jīng)過蒙特卡羅采樣，獲取新的參數(shù)值，重新代入模型運(yùn)行，依次循環(huán)，當(dāng)達(dá)到馬爾科夫收斂條件或設(shè)置的指定馬爾科夫鏈個(gè)數(shù)時(shí)，停止計(jì)算，此時(shí)獲得參數(shù)的后驗(yàn)分布。

圖1 基于MCMC 方法的作物模型參數(shù)校正框架Fig.1 Parameter correction framework of crop model based on Markov chain Monte Carlo (MCMC) method

1.3.3 品種參數(shù)先驗(yàn)分布

模型參數(shù)的先驗(yàn)分布對(duì)用MCMC 方法進(jìn)行參數(shù)后驗(yàn)分布估計(jì)影響不大，一般設(shè)為均勻分布，參數(shù)先驗(yàn)均勻分布的范圍見表2。

1.3.4 模型殘差先驗(yàn)分布

模型殘差是模型模擬值和實(shí)測(cè)值之間的差值，從統(tǒng)計(jì)學(xué)的角度來(lái)看其服從一定的先驗(yàn)分布。模型殘差的先驗(yàn)分布可用概率密度函數(shù)來(lái)表示，概率密度函數(shù)是一種關(guān)于模型殘差特征的假設(shè)，而LF 是以概率密度函數(shù)中各參數(shù)為變量的函數(shù)，MCMC 方法通過不斷采樣模型參數(shù)值，使模型殘差的分布符合這個(gè)LF，因此LF 的形式對(duì)參數(shù)校正的結(jié)果影響較大，關(guān)于LF 形式的討論具體見1.4。

1.3.5 約束目標(biāo)變量

由于不同模型在物候期模擬上存在一定差異，且不同水稻品種的物候期試驗(yàn)數(shù)據(jù)在年份上不完全相同，因此，RiceGrow 物候期模型中，雪花粘品種選擇出苗、拔節(jié)期、抽穗期和成熟期作為目標(biāo)變量，武育粳3 號(hào)品種選擇拔節(jié)期、抽穗期和成熟期作為目標(biāo)變量，汕優(yōu)63 號(hào)品種選擇出苗、拔節(jié)期、抽穗期和成熟期作為目標(biāo)變量；在Oryza2000 物候期模型中，3 個(gè)品種均選擇孕穗期、抽穗期和成熟期作為目標(biāo)變量。

1.3.6 MCMC 采樣方法

研究采用仿射不變馬爾科夫鏈蒙特卡洛集成采樣算法（ensemble sampling for affine-invariant MCMC，EMCEE），該算法由GOODMAN 提出，并由FOREMAN等用python 工具實(shí)現(xiàn)[26]。它引入了仿射不變性采樣器，相對(duì)于傳統(tǒng)的MCMC 方法能產(chǎn)生更多的獨(dú)立樣本，自相關(guān)時(shí)間更短，且利用多個(gè)CPU 內(nèi)核，提高了計(jì)算并行性。目前已經(jīng)在Nature 和Science 等期刊上有關(guān)天體物理學(xué)文獻(xiàn)中被使用[27-28]，但是還未被應(yīng)用于作物模型中。源碼見https://github.com/dfm/emcee。

在EMCEE 算法中，有以下參數(shù)需要設(shè)置：

1）參數(shù)維度。依據(jù)模型參數(shù)而定，本研究中RiceGrow 與Oryza2000 物候期模型帶校正品種參數(shù)均為5 個(gè)，參數(shù)維度設(shè)置為5。

2）并行馬爾科夫鏈初始條數(shù)。EMCEE 算法采用多條并行馬爾科夫鏈進(jìn)行采樣，并行的馬爾科夫鏈初始條數(shù)一般是參數(shù)維度的4～6 倍[27]，本研究中取6 倍，即30。

3）參數(shù)初始值。1 條馬爾科夫鏈需要1 組參數(shù)值，品種參數(shù)初始值的維度等于初始馬爾科夫鏈條數(shù)，依據(jù)步驟2），本研究中取30。每1 組的參數(shù)初始值由參數(shù)先驗(yàn)分布區(qū)間上均勻采樣而來(lái)。

4）馬爾科夫鏈長(zhǎng)度。馬爾科夫鏈長(zhǎng)度為初始馬爾科夫鏈達(dá)到穩(wěn)定時(shí)的條數(shù)，在EMCEE 算法中作為收斂條件，一般依據(jù)經(jīng)驗(yàn)設(shè)定，本研究中設(shè)為1 000。

1.4 LF 形式

1.4.1 GLF

GLF 不考慮模型殘差的異方差性，其分布服從均值為0，方差為常數(shù)的正態(tài)高斯分布。并可用概率密度函數(shù)表示為

式中σ為標(biāo)準(zhǔn)差，各關(guān)鍵物候服從以觀測(cè)日期CCD 為期望的高斯分布。為了使值穩(wěn)定、函數(shù)形式簡(jiǎn)單、計(jì)算方便，所有概率密度的計(jì)算均通過取對(duì)數(shù)形式計(jì)算。模型殘差εi為

各關(guān)鍵物候期觀測(cè)值聯(lián)合GLF 為

式中n為關(guān)鍵物侯期個(gè)數(shù)，是第i個(gè)關(guān)鍵物侯期儒歷天數(shù)模型模擬值，Si是第i個(gè)關(guān)鍵物侯期儒歷天數(shù)模型模擬值，Oi是第i個(gè)關(guān)鍵物侯期儒歷天數(shù)觀測(cè)值。本研究中標(biāo)準(zhǔn)差σ取Oi的1%。

1.4.2 GLF-CV

在實(shí)際物侯期觀測(cè)中，即使觀測(cè)方法和觀測(cè)設(shè)備相同，不同水稻品種歷年關(guān)鍵物候期觀測(cè)值的方差也并不恒定，這是模型殘差異方差性重要來(lái)源之一。

本研究中，汕優(yōu)63 號(hào)品種成熟期的觀測(cè)數(shù)據(jù)出現(xiàn)較多異常值，且中位數(shù)分布在上、下四分位數(shù)之外，說(shuō)明該物候期觀測(cè)數(shù)據(jù)離散程度較大。受不同年份季節(jié)氣候的影響，同一水稻品種關(guān)鍵物侯期觀測(cè)值的方差和均值隨年份變化而變化，因此在使用GLF 時(shí)可能存在偏差。

平穩(wěn)性和離散程度是數(shù)據(jù)異方差性是否存在的直接反映，一般用變異系數(shù)CV 表示，雖然水稻物候期觀測(cè)數(shù)據(jù)的方差會(huì)隨年份變化，但方差和均值的比值在理想情況下應(yīng)該接近恒定的值[29]，本研究引入CV 來(lái)對(duì)模型殘差的異方差性進(jìn)行修正。CV 定義為方差和均值的比值，可得修正后的各關(guān)鍵物候期聯(lián)合GLF-CV：

CV 根據(jù)每個(gè)關(guān)鍵物候期的歷年觀測(cè)值確定，將觀測(cè)值變化的方差轉(zhuǎn)為恒定的CV 值與均值的乘積，一定程度上達(dá)到修正模型殘差異方差性的作用。

1.4.3 GLF-BC

水稻物候期模型的非線性、不連續(xù)、非凸的特點(diǎn)[6]也是模型殘差異方差性的重要來(lái)源。為了同時(shí)考慮模型結(jié)構(gòu)和觀測(cè)數(shù)據(jù)帶來(lái)的異方差性，本研究引入GLF-BC進(jìn)行改善。

GLF-BC 的思想是BC 變換，它是統(tǒng)計(jì)學(xué)中一種通過數(shù)學(xué)變換手段改善數(shù)據(jù)異方差性的有效方法，主要特點(diǎn)是引入一個(gè)參數(shù)，通過數(shù)據(jù)本身估計(jì)該參數(shù)進(jìn)而確定應(yīng)采取的數(shù)據(jù)變換形式，通過BC 變換可以將模型殘差轉(zhuǎn)化為獨(dú)立相同分布的特性[30]。BC 變換的一般形式為

式中λ是變換參數(shù)，y是模型模擬或觀察到的結(jié)果。經(jīng)過BC 變換后，模型的殘差表示為

將式（5）中的εi替換為式（9），可得各關(guān)鍵物候期觀測(cè)值聯(lián)合GLF-BC 為

當(dāng)λ=1 時(shí)，BC 變換無(wú)效，即模型殘差無(wú)變換，當(dāng)λ=0 時(shí)，BC 變換為對(duì)數(shù)變換，而通過觀測(cè)值確定時(shí)常常結(jié)果在接近邊界（通常靠近1）[31]，因此為了使BC 變換有效同時(shí)方便計(jì)算，本研究中使用BC 變換中的平方根轉(zhuǎn)換，即取λ=0.5。

1.5 試驗(yàn)環(huán)境與設(shè)計(jì)

1.5.1 試驗(yàn)環(huán)境

物侯期預(yù)測(cè)模型程序均為用Java 語(yǔ)言自主開發(fā)具有Restful 風(fēng)格的web 服務(wù)，主體算法程序使用python 語(yǔ)言開發(fā)，采樣算法使用EMCEE 工具包，試驗(yàn)運(yùn)行環(huán)境是Intel（R）Xeon（R）Platinum8163CPU@2.50 GHz，內(nèi)存16 GB，Windows10 64 位操作系統(tǒng)。

1.5.2 試驗(yàn)設(shè)計(jì)

設(shè)計(jì)了兩組試驗(yàn)分析比較了不同LF 對(duì)用MCMC 方法進(jìn)行模型參數(shù)校正的影響。數(shù)據(jù)源為3 個(gè)水稻品種歷年站點(diǎn)數(shù)據(jù)（表1），模型包括RiceGrow 和Oryza2000物候期模型。LF 分別為GLF、GLF-CV 和GLF-BC，結(jié)果分別對(duì)校正后的模型參數(shù)后驗(yàn)分布、模型參數(shù)UR 和模型預(yù)測(cè)UR 進(jìn)行比較。其中，GLF-CV 是在GLF 基礎(chǔ)上改善了σ，使觀測(cè)數(shù)據(jù)方差趨于穩(wěn)定，GLF-BC 在GLF 基礎(chǔ)上同時(shí)改善了σ和εi，從降低模型結(jié)構(gòu)非線性造成的異方差進(jìn)行修正。

1.5.3 試驗(yàn)評(píng)價(jià)指標(biāo)

評(píng)價(jià)方法包括參數(shù)校正后驗(yàn)分布、參數(shù)UR 和模型預(yù)測(cè)UR 評(píng)價(jià)。

參數(shù)后驗(yàn)分布用概率密度核函數(shù)圖表示，由多組參數(shù)向量通過曲線擬合而成。

參數(shù)屬于無(wú)量綱的值，參數(shù)UR 是一種參數(shù)校正結(jié)果可信賴程度的量化指標(biāo)，用均方根偏差（root mean square deviation，RMSD）和相對(duì)均方根偏差（relative root mean square deviation，RRMSD）表示[32]，

式中p為后驗(yàn)分布中參數(shù)集個(gè)數(shù)，θi表示第i個(gè)參數(shù)，為參數(shù)集均值。當(dāng)θ服從高斯分布時(shí)，RMSD 可以作為總標(biāo)準(zhǔn)差的無(wú)偏估計(jì)。為了比較不同參數(shù)集的離散程度和穩(wěn)定性，計(jì)算參數(shù)后驗(yàn)分布的RRMSD，計(jì)算式為

模型預(yù)測(cè)UR 用均方根誤差（root mean square error，RMSE）、標(biāo)準(zhǔn)均方根誤差（normalized root mean square error，NRMSE）來(lái)表示。

式中q為關(guān)鍵物侯期個(gè)數(shù)。

2 結(jié)果與分析

2.1 參數(shù)校正結(jié)果比較

2.1.1 RiceGrow 物候期模型參數(shù)后驗(yàn)分布與UR 比較

1）參數(shù)后驗(yàn)分布比較

GLF、GLF-BC 和GLF-CV 下RiceGrow 物侯期模型雪花粘、武育粳3 號(hào)、汕優(yōu)63 號(hào)品種參數(shù)后驗(yàn)分布的概率密度核函數(shù)圖見圖2。可以看出參數(shù)后驗(yàn)分布收斂區(qū)間相對(duì)于先驗(yàn)分布均有縮小，3 種LF 在RiceGrow 物候期模型參數(shù)校正均有一定的效果。

圖2 不同似然函數(shù)下RiceGrow 物候期模型參數(shù)后驗(yàn)分布概率密度核函數(shù)圖Fig.2 Posterior distribution probability density kernel function of RiceGrow phenological model parameters under different likelihood functions

品種參數(shù)TS、FDF、To、IE 與品種熟性無(wú)明顯關(guān)系，PS 代表了水稻的感光性，PS 值越高代表感光性越強(qiáng)，一般來(lái)說(shuō)早熟品種基本不感光或感光極弱，晚熟品種感光較強(qiáng)，而中熟品種感光性較復(fù)雜，雪花粘、武育粳3 號(hào)、汕優(yōu)63 號(hào)分別對(duì)應(yīng)早熟、中熟和晚熟品種。由圖2 可以看出，3 種LF 下雪花粘PS 概率密度函數(shù)峰值對(duì)應(yīng)的參數(shù)值均小于汕優(yōu)63 號(hào)，這與實(shí)際較為吻合。在汕優(yōu)63 號(hào)中，GLF-BC 下PS 概率密度函數(shù)的區(qū)間較GLF-CV 和GLF 更為收斂，表明GLF-BC 下得出的PS后驗(yàn)分布更為精確。

2）參數(shù)UR 比較

由表3 可以看出，GLF-BC 所得的各種參數(shù)RRMSD最小，其次是GLF-CV，GLF 表現(xiàn)最不理想。這是因?yàn)槎鳪LF-CV 只考慮觀測(cè)數(shù)據(jù)帶來(lái)的異方差，未考慮模型結(jié)構(gòu)帶來(lái)的異方差，GLF-BC 同時(shí)考慮了這兩種異方差的來(lái)源，因此GLF-BC 下的參數(shù)RRMSD 小于GLF-CV，GLF-CV 小于GLF，這表明了GLF-BC 和GLF-CV 在改善殘差異方差性上具有一定的作用，也從反面證明了RiceGrow 物候期模型的異方差性與模型結(jié)構(gòu)本身有關(guān)。

表3 不同似然函數(shù)下RiceGrow 物候期模型參數(shù)不確定度Table 3 Uncertainty ratio of RiceGrow phenological model parameters under different likelihood functions

2.1.2 Oryza2000 物候期模型參數(shù)后驗(yàn)分布與UR 比較

1）參數(shù)后驗(yàn)分布比較

3 種LF 下Oryza2000 物侯期模型參數(shù)后驗(yàn)概率密度核函數(shù)圖見圖3，品種參數(shù)DVRJ、DVRP、DVRR、PPSE與品種熟性無(wú)明顯關(guān)系，DVRI 是水稻感光性的倒數(shù)，DVRI 值越小代表感光性越強(qiáng)，圖中可以看出3 種LF 下雪花粘DVRI 概率密度函數(shù)峰值對(duì)應(yīng)的參數(shù)值均大于汕優(yōu)63 號(hào)，這與實(shí)際較為吻合。在汕優(yōu)63 號(hào)中，GLFBC 下DVRI 概率密度函數(shù)的區(qū)間較GLF-CV 和GLF 更為收斂，且峰值對(duì)應(yīng)的參數(shù)值更小，表明GLF-BC 下得出的DVRI 后驗(yàn)分布更為精確。

圖3 不同似然函數(shù)下Oryza2000 物候期模型參數(shù)后驗(yàn)分布概率密度核函數(shù)圖Fig.3 Posterior distribution probability density kernel function of Oryza2000 phenological model parameters under different likelihood functions

2）模型參數(shù)UR 比較

由表4 可知，GLF 所得的雪花粘品種參數(shù)RRMSD最小，GLF-BC 所得的汕優(yōu)63 號(hào)品種參數(shù)RRMSD 最小；GLF-CV 所得的武育粳3 號(hào)品種參數(shù)RRMSD 最小，3 種LF 所得的各品種參數(shù)RRMSD 均有差別。這可能是因?yàn)樯莾?yōu)63 號(hào)品種的觀測(cè)數(shù)據(jù)年份最長(zhǎng)（見表1），根據(jù)統(tǒng)計(jì)學(xué)理論觀點(diǎn)，觀測(cè)數(shù)據(jù)越多，其數(shù)據(jù)特征越趨于穩(wěn)定，GLF-CV 的假設(shè)即是觀測(cè)數(shù)據(jù)具有固定的變異系數(shù)，這與觀測(cè)數(shù)據(jù)具有穩(wěn)定性是一致的。同時(shí)由于在RiceGrow物候期模型利用二次曲線函數(shù)來(lái)描述每日光周期效應(yīng)，Beta 函數(shù)描述每日熱效應(yīng)，這兩種函數(shù)均是非線性函數(shù)，而Oryza2000 物候期模型結(jié)構(gòu)采用多個(gè)線性函數(shù)進(jìn)行級(jí)聯(lián)，其非線性較RiceGrow 弱。若其異方差大部分來(lái)源于模型結(jié)構(gòu)，GLF-BC 好于GLF-CV 是有可能的。若觀測(cè)數(shù)據(jù)的異方差性整體較小且Oryza2000 模型非線性較弱，GLF 好于GLF-CV、GLF-BC 也是有可能的。

2.2 模型預(yù)測(cè)UR 比較

將參數(shù)后驗(yàn)分布中的均值帶入模型運(yùn)行得到模擬結(jié)果與觀測(cè)值進(jìn)行比較分析。在RiceGrow 物侯期模型中，選取出苗期、拔節(jié)期、抽穗期、成熟期4 個(gè)關(guān)鍵物候期進(jìn)行比較；在Oryza2000 物侯期模型中，選取孕穗期、抽穗期、成熟期3 個(gè)關(guān)鍵物候期進(jìn)行比較。結(jié)果見表5。

表5 不同似然函數(shù)下2 種物候期模型預(yù)測(cè)UR 及參數(shù)UR 比較Table 5 Comparison of two phenophase models for predicting uncertainty ratio and parameters uncertainty ratio under different likelihood functions

GLF-BC 下雪花粘、武育粳3 號(hào)、汕優(yōu)63 號(hào)品種的RiceGrow 物候期模型預(yù)測(cè)RMSE 為3.34、3.49、2.66 d；整體上均小于GLF-CV 下的3.43、3.56、3.46 d 和GLF 下的4.54、3.70、2.73 d，這與前一節(jié)得出的參數(shù)UR 中GLFBC 小于GLF-CV 和GLF 對(duì)應(yīng)，說(shuō)明了在RiceGrow 物候期模型中，參數(shù)RRMSD 越小，模型預(yù)測(cè)RMSE 越小，整體上GLF-BC 表現(xiàn)最好，GLF-CV 其次，GLF 最后。

雪花粘、武育粳3 號(hào)、汕優(yōu)63 號(hào)品種的Oryza2000物候期模型預(yù)測(cè)RMSE 最小的分別是GLF、GLF-BC、GLF-CV，并沒有像RiceGrow 物候期模型中出現(xiàn)GLFBC 一直表現(xiàn)良好的現(xiàn)象，這是因?yàn)長(zhǎng)F 可以看作是一種描述模型預(yù)測(cè)UR 的方法，而模型預(yù)測(cè)UR 來(lái)源包括模型結(jié)構(gòu)UR、模型參數(shù)UR、觀測(cè)數(shù)據(jù)UR 等，不同的LF 體現(xiàn)的各類UR 來(lái)源權(quán)重不同，對(duì)模型的匹配程度也不盡相同，相對(duì)于RiceGrow，Oryza2000 物候期模型的結(jié)構(gòu)非線性程度和復(fù)雜性較低，對(duì)于參數(shù)校正而言，觀測(cè)數(shù)據(jù)UR 對(duì)模型預(yù)測(cè)UR 的影響較大，而不同品種和年份的觀測(cè)數(shù)據(jù)存在一定差異，因此在Oryza2000 物候期模型中，LF 對(duì)于不同品種參數(shù)校正的適應(yīng)性不同。

取每個(gè)LF 下的所有品種的參數(shù)RRMSD 和模型預(yù)測(cè)RMSE 的值，得到參數(shù)整體UR 和模型預(yù)測(cè)UR 的量化關(guān)系，從表5 可以看出，在RiceGrow 物候期模型中，整體上，參數(shù)RRMSD 越小，模型預(yù)測(cè)RMSE 越小，這可能是因?yàn)槠鋮?shù)UR 是模型預(yù)測(cè)UR 的主要來(lái)源。而在Oryza2000 物候期模型中，GLF 的模型預(yù)測(cè)RMSE 最小，但其參數(shù)RRMSD 卻最大，模型預(yù)測(cè)的UR 并未隨著參數(shù)UR 同方向變化，這可能是因?yàn)槠淠Ｐ蛥?shù)UR占據(jù)模型預(yù)測(cè)UR 的來(lái)源權(quán)重較小。

對(duì)比兩種模型結(jié)構(gòu)，在揭示水稻發(fā)育進(jìn)程的規(guī)律中，Oryza2000 物候期模型對(duì)于光周期效應(yīng)和熱效應(yīng)的描述均采用線性方程，相對(duì)于RiceGrow 更為簡(jiǎn)化，因此同樣的觀測(cè)數(shù)據(jù)，在Oryza2000 中對(duì)于參數(shù)校正結(jié)果的影響較大，因此其觀測(cè)數(shù)據(jù)UR 可能占據(jù)模型預(yù)測(cè)UR 的來(lái)源權(quán)重較大。而參數(shù)UR 對(duì)模型預(yù)測(cè)UR 的影響，只有在參數(shù)UR 占據(jù)主要來(lái)源權(quán)重時(shí)，才會(huì)呈同趨勢(shì)變化。

3 討論

在RiceGrow 物候期模型中，用MCMC 方法進(jìn)行參數(shù)校正時(shí)，GLF-BC 好于 GLF 和 GLF-CV，然而Oryza2000 物候期模型中不同LF 的效果更依賴于觀測(cè)數(shù)據(jù)，從模型結(jié)構(gòu)來(lái)看，這是由于Oryza2000 物候期模型的線性程度較RiceGrow 高。引入GLF-BC 和GLF-CV的目的是為了改善模型殘差的方差特征，使模型殘差符合正態(tài)高斯分布，這對(duì)非線性系統(tǒng)模型來(lái)說(shuō)是有一定效果的，在線性系統(tǒng)模型中，模型殘差一般均為隨機(jī)誤差，且符合正態(tài)高斯分布，一般不需要進(jìn)行變換。因此將GLF-BC 和GLF-CV 用于MCMC 方法進(jìn)行Oryza2000 物候期模型參數(shù)校正和不確定性分析效果可能不明顯，這需要更多的試驗(yàn)進(jìn)行驗(yàn)證。

在GLF-BC 中，本文假設(shè)的BC 變換系數(shù)值取值為1/2，這是參考一般的BC 均值變換，雖然在RiceGrow物候期模型中效果良好，但是若在MCMC 采樣過程中能夠結(jié)合觀測(cè)數(shù)據(jù)動(dòng)態(tài)調(diào)整，將進(jìn)一步提高LF 的適用性，比如在Oryza2000 物候期模型中使用自適應(yīng)的LF，這是下一步將要研究的問題。

模型結(jié)構(gòu)UR、參數(shù)UR 和觀測(cè)數(shù)據(jù)UR 是模型預(yù)測(cè)UR 的重要來(lái)源，本研究初步量化了參數(shù)UR 與預(yù)測(cè)UR 的關(guān)系，但并未對(duì)模型結(jié)構(gòu)UR、觀測(cè)數(shù)據(jù)UR 與模型預(yù)測(cè)UR 的關(guān)系進(jìn)行量化，這是下一步的工作。

4 結(jié)論

本研究通過引入引入變異系數(shù)（coefficient of variation，CV）變換的高斯似然函數(shù)（Gaussian likelihood function with CV transformation，GLF-CV）和 BC（Box-Cox）變換的高斯似然函數(shù)（GLF with BC transformation，GLF-BC）對(duì)水稻物候期模型殘差的異方差性進(jìn)行修正，以RiceGrow 和Oryza2000 物候期模型為研究對(duì)象，通過試驗(yàn)對(duì)比了3 種似然函數(shù)（likelihood function，LF）下用馬爾科夫鏈蒙特卡洛（Markov chain Monte Carlo，MCMC）算法進(jìn)行雪花粘、武育粳3 號(hào)、汕優(yōu)63 號(hào)水稻品種參數(shù)校正和不確定性量化分析結(jié)果。得出主要結(jié)論如下：

1）引入的GLF-BC 和GLF-CV 在用MCMC 方法進(jìn)行水稻物候期模型參數(shù)校正時(shí)均有效果，得出的雪花粘、武育粳3 號(hào)、汕優(yōu)63 號(hào)參數(shù)后驗(yàn)分布均值帶入RiceGrow和Oryza2000 進(jìn)行模型預(yù)測(cè)，均方根誤差范圍分別為2.66～4.54、2.30～4.41 d。

2）在RiceGrow 物候期模型中，3 個(gè)水稻品種參數(shù)相對(duì)均方根偏差RRMSD（relative root mean square deviation，RRMSD）和模型預(yù)測(cè)均方根誤差（root mean square error，RMSE）均是GLF-BC 最小，在GLFBC 下模型預(yù)測(cè) RMSE 比 GLF-CV 小 0.09、0.07、0.80 d，比GLF 小1.21、0.20、0.07 d，表明GLF-BC 對(duì) RiceGrow物候期模型具有良好的適應(yīng)性。在Oryza2000 物候期模型中，雪花粘、武育粳3 號(hào)、汕優(yōu)63 號(hào)3 個(gè)水稻品種的模型預(yù)測(cè)RMSE 最小的是GLF、GLF-BC 和GLF-CV，分別為2.30、4.17、3.50 d，3 種LF 各有優(yōu)勢(shì)。

3）初步量化了不同LF 下模型預(yù)測(cè)不確定度（uncertainty ratio，UR）和參數(shù)UR 之間的關(guān)系。在RiceGrow 物候期模型中，參數(shù)UR 是模型預(yù)測(cè)UR 的主要來(lái)源，在Oryza2000 物候期模型中，觀測(cè)數(shù)據(jù)UR 可能是模型預(yù)測(cè)UR 的主要來(lái)源。MCMC 通過對(duì)采樣得到的參數(shù)進(jìn)行統(tǒng)計(jì)分析，估計(jì)模型參數(shù)的后驗(yàn)分布，后驗(yàn)分布反映了參數(shù)UR。而LF 在參數(shù)校正中起到了關(guān)鍵的作用，參數(shù)后驗(yàn)分布的形狀和位置受LF 的影響，因此LF 的選擇和定義可能因具體問題而異，需要能夠與模型和觀測(cè)數(shù)據(jù)的特性相匹配。LF 的選擇與模型殘差異方差的主要來(lái)源有關(guān)，當(dāng)主要來(lái)源為觀測(cè)數(shù)據(jù)時(shí)，GLF-CV好于其他；當(dāng)主要來(lái)源為模型結(jié)構(gòu)本身時(shí)，GLF-BC 好于其他；當(dāng)模型殘差的異方差性較小時(shí)，可使用GLF。