本發(fā)明屬于農(nóng)業(yè)信息
技術(shù)領(lǐng)域：
：，尤其是一種基于r語言的小麥生育期特征參數(shù)估算方法。
背景技術(shù)：
：：自工業(yè)革命以來，全球氣候變化已經(jīng)成為一個不爭的事實。全球范圍的氣候變化不僅影響人類賴以生存的環(huán)境，也給農(nóng)業(yè)生產(chǎn)帶來了嚴(yán)峻的挑戰(zhàn)。隨著對全球氣候變化研究的不斷深入，有關(guān)氣候變化對作物生長發(fā)育及產(chǎn)量形成影響的研究也開始受到廣泛的關(guān)注。在未來氣候變化情景下，及時準(zhǔn)確的預(yù)測作物生長情況，對指導(dǎo)及調(diào)控農(nóng)業(yè)生產(chǎn)、提高管理水平、制定適宜的糧食生產(chǎn)政策等意義重大。作物生產(chǎn)系統(tǒng)是受氣候、土壤、品種以及栽培管理措施等因素影響的系統(tǒng)，從而使農(nóng)業(yè)生產(chǎn)者難以綜合考慮多因子互作、預(yù)測農(nóng)業(yè)生產(chǎn)趨勢并量化各因子對農(nóng)作物生產(chǎn)的影響。作物生長模擬模型為定量描述作物生長發(fā)育過程及其與環(huán)境和技術(shù)的動態(tài)關(guān)系提供了新的方法和手段。但作物模擬模型在應(yīng)用時存在著一定的不確定性。這些不確定性主要來源于三個方面，即輸入變量導(dǎo)致的不確定性，參數(shù)值導(dǎo)致的不確定性以及模型公式導(dǎo)致的不確定性。不確定性分析是基于模型的風(fēng)險分析和決策的重要組成部分，它能夠給風(fēng)險分析人員和決策支持者提供模型輸出結(jié)果的準(zhǔn)確性程度。參數(shù)與作物模型密切相關(guān)，模型中的參數(shù)值一般難以直接獲得，并且由于獲得參數(shù)值的方法不同以及輸入數(shù)據(jù)的不同，從而導(dǎo)致模擬結(jié)果存在較大的不確定性，增加了使用者基于模型決策的不確定性。然而，已有對作物模擬模型不確定性的研究大都將參數(shù)值、輸入變量及模型公式這三種不確定性來源混合在一起進(jìn)行分析，沒有研究參數(shù)不確定性對模擬結(jié)果不確定性產(chǎn)生的影響。技術(shù)實現(xiàn)要素：本發(fā)明所解決的技術(shù)問題在于提供一種基于r語言的小麥生育期特征參數(shù)估算方法，利用r語言對小麥生長模擬模型進(jìn)行重新編寫，結(jié)合非線性最小二乘法，判斷開花期和成熟期的實測值與模擬值差方的和的最小值，計算得到差方的和的最小值對應(yīng)的參數(shù)調(diào)試結(jié)果，再利用r語言計算得到小麥品種的參數(shù)范圍。實現(xiàn)本發(fā)明目的的技術(shù)解決方案為：一種基于r語言的小麥生育期特征參數(shù)估算方法，包括以下步驟：步驟1：采集n組小麥生長模擬模型wheatgrow的生育期子模塊所需的生育期數(shù)據(jù)和氣象數(shù)據(jù)；步驟2：利用拉丁超立方體為小麥生育期相關(guān)的小麥品種特征參數(shù)調(diào)試設(shè)置10組不同的初始值；步驟3：采用r語言對小麥生長模擬模型的生育期子模塊進(jìn)行重新編寫，并將輸出的儒歷日期連續(xù)化；步驟4：將田間觀測的某一品種nsites個種植站點的氣象數(shù)據(jù)通過r語言讀取到小麥生長模擬模型中，調(diào)用第i組小麥品種特征參數(shù)初始值，運行小麥生長模擬模型n次，得到模擬的n次小麥主要生育期結(jié)果，其中，i＝1,2，…，10，nsites是數(shù)據(jù)的站點數(shù)；步驟5：將田間實測的n次小麥主要生育期數(shù)據(jù)和模擬的n次小麥主要生育期結(jié)果相結(jié)合，計算生育期中開花期和成熟期實測值與模擬值差方的和步驟6：若達(dá)到最小值，計算得到參數(shù)調(diào)試結(jié)果并保存，利用r語言中vcov()函數(shù)得到參數(shù)間的協(xié)方差，利用cov2cor(vcov())函數(shù)得到參數(shù)間的相關(guān)系數(shù)，若i＜10則令i＝i+1并轉(zhuǎn)到步驟4，若i＝10則轉(zhuǎn)到步驟7；若未達(dá)到最小值，則在第i組參數(shù)的基礎(chǔ)上產(chǎn)生新的第i組參數(shù)，返回步驟4；步驟7：計算i站點j年的模型誤差和其中為播種-開花的天數(shù)誤差，為開花-成熟的天數(shù)誤差，再利用r語言中的kendall秩相關(guān)函數(shù)計算和的相關(guān)性大小p；步驟8：若p＜0.5，則和是相關(guān)獨立的，步驟6中得到的參數(shù)調(diào)試結(jié)果為最終結(jié)果；若p≥0.5，則對各站點的原始數(shù)據(jù)進(jìn)行轉(zhuǎn)化，使轉(zhuǎn)化后的模型誤差相關(guān)性p＜0.5，再轉(zhuǎn)到步驟5進(jìn)行參數(shù)調(diào)試并得到最終調(diào)試結(jié)果；步驟9：選取10套參數(shù)調(diào)試結(jié)果中最小的參數(shù)組作為最優(yōu)值，結(jié)合參數(shù)間協(xié)方差，利用r語言中的正太分布函數(shù)mvrnorm(1000,mu,sigma)得到1000套參數(shù)值并計算其概率分布，取80％區(qū)間為置信區(qū)間，該區(qū)間內(nèi)的參數(shù)值為該品種的參數(shù)范圍，其中，mu為最優(yōu)值，sigma為參數(shù)間協(xié)方差。進(jìn)一步的，本發(fā)明的基于r語言的小麥生育期特征參數(shù)估算方法，步驟2中的小麥生育期相關(guān)的小麥品種特征參數(shù)包括：生理春化天數(shù)、光周期敏感性、基本早熟性、溫度敏感性和灌漿期因子。。進(jìn)一步的，本發(fā)明的基于r語言的小麥生育期特征參數(shù)估算方法，步驟3中采用在開花期和成熟期用一個帶小數(shù)位的天數(shù)替代整數(shù)天數(shù)來實現(xiàn)小麥生長模擬模型輸出的儒歷日期連續(xù)化。進(jìn)一步的，本發(fā)明的基于r語言的小麥生育期特征參數(shù)估算方法，步驟5中開花期和成熟期實測值與模擬值差方的和的計算公式為：其中，ni是i站點的數(shù)據(jù)年份數(shù)，和分別是i站點j年的播種-開花天數(shù)和開花-成熟天數(shù)的實測值，和分別是i站點j年的播種-開花天數(shù)和開花-成熟天數(shù)的模擬值。進(jìn)一步的，本發(fā)明的基于r語言的小麥生育期特征參數(shù)估算方法，步驟6中采用非線性最小二乘法判斷是否達(dá)到最小值。進(jìn)一步的，本發(fā)明的基于r語言的小麥生育期特征參數(shù)估算方法，步驟7中的和的計算公式為：其中，μant和μmat分別是導(dǎo)致播種-開花和開花-成熟模擬誤差的固定值，和分別是由站點效應(yīng)導(dǎo)致播種-開花和開花-成熟模擬誤差的隨機(jī)值，和分別是播種-開花和開花-成熟的殘差。進(jìn)一步的，本發(fā)明的基于r語言的小麥生育期特征參數(shù)估算方法，步驟8中的對各站點的原始數(shù)據(jù)進(jìn)行轉(zhuǎn)化的步驟包括：步驟8-1：將和的值利用r語言中的lmer函數(shù)計算出誤差的方差-協(xié)方差矩陣∑；步驟8-2：采用r語言中的cholesky對∑進(jìn)行分解：∑＝rtr，得到rt，利用rt對原始數(shù)據(jù)y進(jìn)行轉(zhuǎn)化：y*＝(rt)-1y，得到轉(zhuǎn)化后的數(shù)據(jù)y*。本發(fā)明采用以上技術(shù)方案與現(xiàn)有技術(shù)相比，具有以下技術(shù)效果：1、本發(fā)明的方法通過基于r語言的非線性最小二乘法和小麥生長模型相結(jié)合，有效反演出小麥品種特征參數(shù)；2、本發(fā)明的方法在計算生育期時，不僅將各個生育期看作獨立分布的，而且考慮到各生育期之間的相互作用，從而確保參數(shù)調(diào)試結(jié)果的準(zhǔn)確性；3、本發(fā)明的方法采用的非線性最小二乘法提高了抽樣的效率，避免了參數(shù)陷入局部最優(yōu)；4、本發(fā)明的方法基于貝葉斯理論框架，以較大值為搜索起始點，不斷從兩端向中間尋找誤差最小的值，有效地考慮了參數(shù)的先驗分布，通過轉(zhuǎn)化原始數(shù)據(jù)解決了站點之間誤差相關(guān)性的問題，使得參數(shù)調(diào)試結(jié)果準(zhǔn)確高效。附圖說明圖1是本發(fā)明的方法流程圖；圖2是小麥生長模擬模型的框架結(jié)構(gòu)圖；圖3是京冬8號小麥的生育期特征參數(shù)概率分布圖；圖4是京冬8號小麥的生育期模擬值與實測值的1：1關(guān)系圖。具體實施方式下面詳細(xì)描述本發(fā)明的實施方式，所述實施方式的示例在附圖中示出，其中自始至終相同或類似的標(biāo)號表示相同或類似的元件或具有相同或類似功能的元件。下面通過參考附圖描述的實施方式是示例性的，僅用于解釋本發(fā)明，而不能解釋為對本發(fā)明的限制。利用本發(fā)明的基于r語言的小麥生育期特征參數(shù)估算方法對1980-2010年在12個站點種植63次的冬小麥品種京冬8號的生育期數(shù)據(jù)進(jìn)行有效地估計，并結(jié)合圖1的流程進(jìn)行詳細(xì)說明，如圖1所示具體包括以下步驟：步驟1：采集63組小麥生長模擬模型wheatgrow(如圖2所示)所需的生育期數(shù)據(jù)和氣象數(shù)據(jù)，包括12個站點。其中，采集氣象數(shù)據(jù)：搜集整理1980-2010年研究站點的逐日氣象數(shù)據(jù)，涉及日最高氣溫、日最低氣溫、日平均氣溫、日照時數(shù)和日降水量5個氣象要素，來自于中國氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)(http://cdc.cma.gov.cn/)；采集小麥生育期數(shù)據(jù)：包括1980-2010京冬8號生育期資料(包括播種期、開花期、成熟期)。數(shù)據(jù)均來自國家氣象中心資料室。冬小麥京冬8號基本種植信息見表1所示：表1京冬8號實測生育期步驟2：利用拉丁超立方體為小麥生育期相關(guān)的小麥品種特征參數(shù)調(diào)試設(shè)置10組不同的初始值，小麥品種特征參數(shù)具體包括：生理春化天數(shù)(pvt)、光周期敏感性(ps)、基本早熟性(ie)、溫度敏感性(ts)和灌漿期因子(fdf)，具體公式為：paratable[i,1]＜-runif(1,0.6+c1[i]/(10/0.9)-0.12,0.6+c1[i]/(10/0.9))paratable[i,2]＜-runif(1,c2[i]/(10/3)-0.3,c2[i]/(10/3))paratable[i,3]＜-runif(1,c3[i]/(10/0.01)-0.001,c3[i]/(10/0.01))paratable[i,4]＜-runif(1,c4[i]/(10/60)-6,c4[i]/(10/60))paratable[i,5]＜-runif(1,0.8+c5[i]/(10/0.7)-0.12,0.8+c5[i]/(10/0.7))paratable[i,1]為ie，paratable[i,2]為ts，paratable[i,3]為ps，paratable[i,4]為pvt，paratable[i,5]為fdf，i＝1,2,…,10。得到的初始值如表2所示：表2京冬8號10組初始參數(shù)值步驟3：采用r語言對小麥生長模擬模型的生育期子模塊進(jìn)行重新編寫，并采用在開花期和成熟期用一個帶小數(shù)位的天數(shù)替代整數(shù)天數(shù)將小麥生長模擬模型輸出的儒歷日期連續(xù)化。在本模型中，存在著幾個臨界溫度變化的轉(zhuǎn)折點，臨界溫度變化導(dǎo)致模型公式變化，從而導(dǎo)致模擬值輸出的不連續(xù)性，據(jù)此進(jìn)行連續(xù)化的點見表3所示：表3公式改變時的pdt值步驟4：將田間觀測的某一品種12個種植站點的氣象數(shù)據(jù)通過r語言讀取到小麥生長模擬模型中，調(diào)用第i組小麥品種特征參數(shù)初始值，運行小麥生長模擬模型63次，得到模擬的63次小麥主要生育期結(jié)果，其中，i＝1,2，…，10。步驟5：將田間實測的63次小麥主要生育期數(shù)據(jù)和模擬的63次小麥主要生育期結(jié)果相結(jié)合，如圖4所示，計算生育期中開花期和成熟期實測值與模擬值差方的和的計算公式為：其中，ni是i站點的數(shù)據(jù)年份數(shù)，和分別是i站點j年的播種-開花天數(shù)和開花-成熟天數(shù)的實測值，和分別是i站點j年的播種-開花天數(shù)和開花-成熟天數(shù)的模擬值。步驟6：采用非線性最小二乘法判斷是否達(dá)到最小值，若達(dá)到最小值，計算得到參數(shù)調(diào)試結(jié)果并保存，如表4所示：表4京冬8號10組參數(shù)估計結(jié)果然后利用r語言中vcov()函數(shù)得到參數(shù)間的協(xié)方差，如表5所示：表5京冬8號參數(shù)間協(xié)方差利用cov2cor(vcov())函數(shù)得到參數(shù)間的相關(guān)系數(shù)，如表6所示：表6京冬8號參數(shù)間相關(guān)系數(shù)若i＜10則令i＝i+1并轉(zhuǎn)到步驟4，若i＝10則轉(zhuǎn)到步驟7；若未達(dá)到最小值，則在第i組參數(shù)的基礎(chǔ)上產(chǎn)生新的第i組參數(shù)，返回步驟4；步驟7：計算i站點j年的模型誤差和其中，為播種-開花的天數(shù)誤差，為開花-成熟的天數(shù)誤差，再利用r語言中的kendall秩相關(guān)函數(shù)計算和的相關(guān)性大小p；和的計算公式為：其中，μant和μmat分別是導(dǎo)致播種-開花和開花-成熟模擬誤差的固定值，和分別是由站點效應(yīng)導(dǎo)致播種-開花和開花-成熟模擬誤差的隨機(jī)值，和分別是播種-開花和開花-成熟的殘差。r語言中kendall的秩相關(guān)函數(shù)表達(dá)為：和的計算結(jié)果如表7所示：表7京冬8號實測值與模擬值的差值步驟8：若p＜0.5，則和是相關(guān)獨立的，步驟6中得到的參數(shù)調(diào)試結(jié)果為最終結(jié)果，即表4所示；若p≥0.5，則對各站點的原始數(shù)據(jù)進(jìn)行轉(zhuǎn)化，使轉(zhuǎn)化后的模型誤差相關(guān)性p＜0.5，再轉(zhuǎn)到步驟5進(jìn)行參數(shù)調(diào)試并得到最終調(diào)試結(jié)果，具體的轉(zhuǎn)化步驟包括：步驟8-1：將和的值利用r語言中的lmer函數(shù)計算出誤差的方差-協(xié)方差矩陣∑；步驟8-2：采用r語言中的cholesky對∑進(jìn)行分解：∑＝rtr，得到rt，利用rt對原始數(shù)據(jù)y進(jìn)行轉(zhuǎn)化：y*＝(rt)-1y，得到轉(zhuǎn)化后的數(shù)據(jù)y*。京冬8號品種的相關(guān)性大小p＝0.26<0.5，因此不需要對該品種原始數(shù)據(jù)進(jìn)行轉(zhuǎn)化，以上通過非線性最小二乘法調(diào)試出的結(jié)果即為最終的結(jié)果，即表4所示。從中選取為最小時的參數(shù)值作為最優(yōu)值，如表8所示：表8京冬8號最優(yōu)值步驟9：選取10套參數(shù)調(diào)試結(jié)果中最小的參數(shù)組作為最優(yōu)值，結(jié)合參數(shù)間協(xié)方差，利用r語言中的正太分布函數(shù)mvrnorm(1000,mu,sigma)得到1000套參數(shù)值并計算其概率分布，取80％區(qū)間為置信區(qū)間，該區(qū)間內(nèi)的參數(shù)值為該品種的參數(shù)范圍，如圖3所示，其中，mu為最優(yōu)值，即表8所示，sigma為參數(shù)間協(xié)方差，即表5所示。以上所述僅是本發(fā)明的部分實施方式，應(yīng)當(dāng)指出，對于本
技術(shù)領(lǐng)域：
：的普通技術(shù)人員來說，在不脫離本發(fā)明原理的前提下，還可以做出若干改進(jìn)，這些改進(jìn)應(yīng)視為本發(fā)明的保護(hù)范圍。當(dāng)前第1頁12當(dāng)前第1頁12