本發(fā)明屬于農(nóng)業(yè)遙感領(lǐng)域，具體涉及一種遙感、作物模型與氣象預(yù)報融合的作物成熟期預(yù)測方法。

背景技術(shù)：

傳統(tǒng)的成熟期預(yù)測方法是以野外觀測為基礎(chǔ)的目視觀察法，即直接通過地面定點觀測作物的生長狀況來預(yù)判成熟情況，但由于區(qū)域局限性強，同時需要消耗大量的時間及人力物力，難以獲得尺度范圍作物成熟期的空間分布信息。

近年來，大量研究通過建立基于氣象數(shù)據(jù)的統(tǒng)計模型，研究溫度、光周期、降水等因子對作物的影響，進而預(yù)測作物的成熟期，該模型雖簡單易用，但驅(qū)動參數(shù)少，難以描述地塊性質(zhì)差異的影響，不適宜大范圍推廣。除此之外，由于氣候變化的影響，基于歷史氣象數(shù)據(jù)建立的統(tǒng)計分析模型不能準確地描述未來的氣候情況。

還有一些研究利用衛(wèi)星遙感數(shù)據(jù)監(jiān)測作物的一個或多個物候指標的動態(tài)變化來預(yù)測其成熟期，但單純的遙感方法對遙感數(shù)據(jù)的時間和空間分辨率都有嚴格要求，受到云和衛(wèi)星軌道的影響，獲取大區(qū)域作物關(guān)鍵生育期所需的中等分辨率遙感數(shù)據(jù)難以在數(shù)據(jù)質(zhì)量上滿足要求。

另一種方法是通過作物生長模型的模擬，從作物生長機理出發(fā)描述作物生長發(fā)育與產(chǎn)量形成的過程，并以作物產(chǎn)量或品質(zhì)(或兩者綜合)為目標構(gòu)建代價函數(shù)，反向求解優(yōu)化作物的收獲時間，實現(xiàn)作物成熟期的預(yù)測，但大區(qū)域范圍的作物生長模型標定與校準困難使其存在一定的局限性。

雖然氣象統(tǒng)計模型、作物生長模型與遙感監(jiān)測這三種方式都存在一定弊端，但是如果將三者優(yōu)勢互補、有機集成，利用遙感獲取區(qū)域的農(nóng)田參數(shù)來校準作物生長模型敏感參數(shù)集，并將集合預(yù)報數(shù)據(jù)作為作物生長模型未來時段的氣象輸入，驅(qū)動作物模型預(yù)測冬小麥成熟期，就能夠有效地克服當前基于統(tǒng)計模型的成熟期預(yù)測難以大范圍推廣應(yīng)用的局限性，適合于區(qū)域范圍冬小麥成熟期預(yù)測。

目前，遙感數(shù)據(jù)和作物生長模型數(shù)據(jù)同化以及引入集合預(yù)報數(shù)據(jù)來實現(xiàn)區(qū)域冬小麥成熟期預(yù)測的研究尚未見報道。

技術(shù)實現(xiàn)要素：

為解決現(xiàn)有技術(shù)中存在的如下問題：“如何將大范圍高時效性的遙感數(shù)據(jù)和準確模擬作物生長的機理模型進行同化，并結(jié)合氣象預(yù)報數(shù)據(jù)，進而能夠在大范圍內(nèi)準確地對冬小麥成熟期提前動態(tài)預(yù)測”，本發(fā)明提供一種遙感、作物模型與氣象預(yù)報融合的作物成熟期預(yù)測方法。

本發(fā)明提供一種遙感、作物模型與氣象預(yù)報融合的作物成熟期預(yù)測方法，具體步驟如下：

s1、根據(jù)氣象數(shù)據(jù)、種植模式和產(chǎn)量水平等對作物種植區(qū)域進行分區(qū)；

s2、基于農(nóng)業(yè)氣象站點觀測數(shù)據(jù)，標定各分區(qū)作物生長的wofost(worldfoodstudies：世界糧食研究)模型參數(shù)；

s3、收集反映作物關(guān)鍵生育期的空間分辨率500米、步長4天的modislai產(chǎn)品并合成時間序列曲線，運用濾波方法消除云的污染；

s4、在全局敏感性分析的基礎(chǔ)上，篩選出對成熟期敏感的參數(shù)，構(gòu)建基于歸一化的modislai和wofost模擬lai的代價函數(shù)，采用sce-ua算法最小化代價函數(shù)，獲得優(yōu)化參數(shù)數(shù)據(jù)集；

s5、用s4所得優(yōu)化參數(shù)數(shù)據(jù)集驅(qū)動wofost作物生長模型，并將tigge(交互式全球大集合預(yù)報系統(tǒng))集合預(yù)報數(shù)據(jù)作為wofost未來時段的氣象輸入，提前10天預(yù)測作物成熟期；

s6、逐格網(wǎng)單元運行s5，生成區(qū)域作物成熟期預(yù)測空間分布圖，指導(dǎo)農(nóng)作物的及時采收。

步驟s1中所述種植結(jié)構(gòu)是指通過modisndvi和lai影像分類獲得的種植結(jié)構(gòu)圖。

所述步驟s1中對冬小麥種植區(qū)進行分區(qū)的方法如下：

先分別建立區(qū)域多年環(huán)境指標均值與每年環(huán)境指標值數(shù)據(jù)庫，經(jīng)比較r方和半偏r方值，確定區(qū)域聚類數(shù)目；將區(qū)域劃分為多個固定形狀大小的地理網(wǎng)格，對這些地理網(wǎng)格每年的環(huán)境指標值均進行k-means屬性聚類，多年間出現(xiàn)次數(shù)最高的類作為地理網(wǎng)格最終歸屬類別，然后進行空間連續(xù)性調(diào)整，得到時空型區(qū)劃結(jié)果。

其中，地理網(wǎng)格屬于類別i的概率pi如下公式(1)所示：

所述步驟s3中所采用的濾波方法是savitzky-golay濾波，保留高的數(shù)值而減少異常值，濾波過程按公式(2)求得:

式(2)中y^*表示濾波后的lai值；y表示lai的原始值；m為窗口的半徑；c表示第i個lai值的濾波系數(shù)；n為卷積數(shù)目；窗口的寬度為2m+1。

步驟s4所述全局敏感性分析，采用擴展的傅里葉振幅敏感度檢驗法(efast)對模型作物參數(shù)和土壤與管理參數(shù)進行全局敏感性分析，輸出一階敏感性指數(shù)和總敏感性指數(shù)，按如下公式(3)求得:

st.i＝si+sij+sijm+…+s12…i…k

v(y)＝∑ivi+∑i≠jvij+∑i≠j≠mvijm+v12...k(3)

式中，v(y)為模型結(jié)果y的總方差；vi為參數(shù)xi的方差；vij～v12…k為各參數(shù)相互作用的方差。si、sij、sijm為參數(shù)xi的一階、二階、三階敏感性指數(shù)，st.i為參數(shù)xi的總敏感性指數(shù)。

所述步驟s4中采用sce_ua算法對代價函數(shù)進行最小化，并結(jié)合如權(quán)利要求4所述的歸一化處理，代價函數(shù)按如下公式(4)求得：

其中，j為代價函數(shù)值，laioi代表遙感l(wèi)ai值，laimi代表作物生長模型模擬lai值，laiomax和laismax分別代表遙感或模型模擬lai時間序列曲線上的最大值，laiomin和laismin分別代表其最小值，n代表遙感l(wèi)ai值的觀測次數(shù)。

所述步驟s5中處理tigge集合預(yù)報數(shù)據(jù)時采用encu模型，計算方法如下公式(5)所示：

其中，ecmwf、ncep、cma、ukmo分別為歐洲、美國、中國、英國相關(guān)氣象機構(gòu)的預(yù)報數(shù)據(jù)。

其中，所述作物優(yōu)選為冬小麥。

本發(fā)明還提供所述遙感、作物模型與氣象預(yù)報融合的作物成熟期預(yù)測方法在指導(dǎo)農(nóng)作物生產(chǎn)中的應(yīng)用。

本發(fā)明與現(xiàn)有技術(shù)相比，有益效果為：

本發(fā)明所述的一種基于遙感、作物生長模型與氣象集合預(yù)報融合的區(qū)域冬小麥成熟期預(yù)測方法，融合了遙感數(shù)據(jù)和作物模型的優(yōu)勢，把遙感觀測和模型模擬的時間序列l(wèi)ai變化趨勢提取出來，通過優(yōu)化算法實現(xiàn)遙感和作物模型的同化，既能在區(qū)域尺度上預(yù)測冬小麥的成熟期，又能提高作物模型的預(yù)測精度，而且氣象預(yù)報數(shù)據(jù)的應(yīng)用又提高了動態(tài)預(yù)測成熟期的精度，實用性更強。

附圖說明

圖1為本發(fā)明遙感、作物模型與氣象預(yù)報融合的作物成熟期預(yù)測方法流程示意圖；

圖2為歸一化后modislai和wofost模型模擬lai時間序列曲線對比圖；

圖3為河北、河南、山東三省冬小麥種植區(qū)成熟期提前10天預(yù)報的空間分布圖。

具體實施方式

實施例1

本發(fā)明遙感、作物模型與氣象預(yù)報融合的作物成熟期預(yù)測方法的流程示意圖參見附圖1，以下面以河北、河南、山東三省研究區(qū)為例進一步闡述本發(fā)明的技術(shù)方案。

s1根據(jù)氣象數(shù)據(jù)、種植模式和產(chǎn)量水平等對作物種植區(qū)域進行分區(qū)；

選擇河北、河南、山東三省作為研究區(qū)域，該區(qū)地處黃淮海平原，屬暖溫帶半濕潤季風氣候。獲取以下數(shù)據(jù)：根據(jù)研究區(qū)的外包絡(luò)范圍選取62個氣象站點，獲取氣象預(yù)報數(shù)據(jù)(最高氣溫、最低氣溫、總輻射、水汽壓、風速、降水模型共6個氣象要素的數(shù)據(jù))；從研究區(qū)內(nèi)農(nóng)業(yè)氣象站點獲取采集的土壤參數(shù)、作物參數(shù)和物候數(shù)據(jù)；從歐洲中期天氣預(yù)報中心獲取氣象預(yù)報數(shù)據(jù)(最高氣溫、最低氣溫、總輻射、水汽壓、風速、降水模型共6個氣象要素的數(shù)據(jù))；獲取經(jīng)緯度、高程等控制參數(shù)；獲取modisl4級標準數(shù)據(jù)產(chǎn)品mod15a3。

利用研究區(qū)2011至2015年同時相的tm影像(指上述獲取的mod15a3)解譯得到冬小麥種植區(qū)域，生成1公里的數(shù)據(jù)掩膜，然后利用modis-lai數(shù)據(jù)中包含的mod12土地利用信息，并設(shè)定閾值剔除不符合冬小麥物候規(guī)律的像素生成新的掩膜，將兩種掩膜進行與運算，生成最終的運行區(qū)域。對模型的產(chǎn)量水平和氣象數(shù)據(jù)(6個氣象要素的數(shù)據(jù))統(tǒng)一度量，將氣象數(shù)據(jù)插值生成1公里每像素的柵格數(shù)據(jù)。

分別建立研究區(qū)2011至2015年環(huán)境指標(生育期內(nèi)累積日照時數(shù)、累積活動積溫、累積降雨、高程)均值與每年環(huán)境指標值數(shù)據(jù)庫，比較r方和半偏r方值，確定區(qū)域聚類數(shù)目；將區(qū)域劃分為1km×1km的地理網(wǎng)格，對這些地理網(wǎng)格每年的環(huán)境指標值均進行k-means屬性聚類，多年間出現(xiàn)次數(shù)最高的類作為地理網(wǎng)格最終歸屬類別，然后進行空間連續(xù)性調(diào)整，得到時空型區(qū)劃結(jié)果。

其中，地理網(wǎng)格屬于類別i的概率pi如下公式(1)所示：

s2wofost模型參數(shù)標定

基于農(nóng)業(yè)氣象站點的觀測數(shù)據(jù)，根據(jù)作物參數(shù)的敏感性程度調(diào)整參數(shù)，具體方案如下：

(1)對于不敏感的作物參數(shù)，或者是敏感性較高但參數(shù)變化范圍較小的參數(shù)，通過查閱文獻或直接采用wofost模型的默認值來確定。

(2)對于敏感性較高且參數(shù)變化范圍較大的參數(shù)，通過查閱文獻或利用實測值計算得到，也可以根據(jù)文獻或先驗知識確定參數(shù)的范圍，然后通過wofost模型附帶的優(yōu)化程序(fseopt)確定。本次實施例具體參數(shù)確定如下：tsum1(出苗至開花有效積溫)、idem(出苗日期)、efftb(對應(yīng)生育期單片光能利用效率)、fl01(葉分配系數(shù))、fs01(莖分配系數(shù))、fotb(貯存器官分配系數(shù))。

s3收集反映作物關(guān)鍵生育期的空間分辨率500米、步長4天的modislai產(chǎn)品并合成時間序列曲線，運用savitzky-golay濾波方法消除云的污染。具體如下：

s3.1生成時間序列曲線

利用nasa提供的mrt投影轉(zhuǎn)換工具，將研究區(qū)對應(yīng)遙感數(shù)據(jù)進行鑲嵌、投影轉(zhuǎn)換和格式轉(zhuǎn)換，選取主要生育期1到177天(儒略日)空間分辨率500米、步長4天的modislai數(shù)據(jù)，共45景影像進行疊加，生成每個格網(wǎng)點對應(yīng)的時間序列曲線。

s3.2savitzky-golay濾波處理

為了消除云的污染，利用savitzky-golay濾波解決葉面積指數(shù)變化過程不連續(xù)的問題，其原理表達如下：

式中：y表示lai的原始值；y^*表示濾波后的lai值；c表示第i個lai值的濾波系數(shù)；m為窗口的半徑；n為卷積數(shù)目；窗口的寬度為2m+1，這個多項式的設(shè)計是為了保留高的數(shù)值而減少異常值。具體操作如下：

(1)對初始的lai時間序列進行savitzky-golay濾波，得到的平滑后的結(jié)果。并分別保存平滑前和平滑后的序列；

(2)對比上一步的兩個序列，用上述公式生成新的序列，作為初始序列；

其中o和n分別為初始和濾波后的lai值,t代表迭代次數(shù),并且i是lai時間序列索引。

(3)重復(fù)(1)(2)兩步，直到整個序列σni^t小于指定的閾值0.1。

s3.3模型曲線擬合，提取特征

在區(qū)域范圍運行作物模型。針對生成的最終運行區(qū)域中的所有像素運行作物模型，對結(jié)果中的lai時間序列利用logistic曲線擬合。

s4在全局敏感性分析的基礎(chǔ)上，篩選出對成熟期敏感的參數(shù)，構(gòu)建基于歸一化的modislai和wofost模擬lai的代價函數(shù)，采用sce-ua算法最小化代價函數(shù)，獲得優(yōu)化參數(shù)數(shù)據(jù)集

s4.1全局敏感性分析

根據(jù)wofost模型在研究區(qū)冬小麥的初步標定結(jié)果，借助于敏感性與不確定性專業(yè)分析軟件simlab的efast模塊實現(xiàn)wofost模型的作物參數(shù)全局敏感性分析：首先，根據(jù)設(shè)定各參數(shù)在取值范圍內(nèi)均勻分布，利用蒙特卡洛方法分別隨機采樣3000次和1000次；然后，將生成的參數(shù)組合輸入到wofost模型中獲得冬小麥模擬產(chǎn)量；最后，利用efast方法計算各個參數(shù)的一階敏感型指數(shù)和總敏感性指數(shù)，得到對于冬小麥模擬產(chǎn)量敏感度較高的wofost作物參數(shù)，利用地面樣點實測數(shù)據(jù)對模型敏感參數(shù)進行標定，完成模型的本地化調(diào)整。

輸出一階敏感性指數(shù)和總敏感性指數(shù)，按如下公式(3)求得:

st.i＝si+sij+sijm+…+s12…i…k

v(y)＝∑ivi+∑i≠jvij+∑i≠j≠mvijm+v12...k(3)

式中，v(y)為模型結(jié)果y的總方差；vi為參數(shù)xi的方差；vij～v12…k為各參數(shù)相互作用的方差。si、sij、sijm為參數(shù)xi的一階、二階、三階敏感性指數(shù)，st.i為參數(shù)xi的總敏感性指數(shù)。

s4.2構(gòu)建基于歸一化的modislai和wofost模擬lai的代價函數(shù)

歸一化后modislai和wofost模型模擬lai時間序列曲線對比圖見圖2。

建立優(yōu)化代價函數(shù)如下：

其中，j為代價函數(shù)值，laioi代表遙感l(wèi)ai值，laimi代表作物生長模型模擬lai值，laiomax和laismax分別代表遙感或模型模擬lai時間序列曲線上的最大值，laiomin和laismin分別代表其最小值，n代表遙感l(wèi)ai值的觀測次數(shù)。

對每個像素通過不斷重新初始化wofost模型敏感參數(shù)，使得模型輸出的lai數(shù)值不斷發(fā)生變化。然后在代價函數(shù)中對比modis提取和模型輸出的lai數(shù)值，sce_ua算法用來在初始參數(shù)空間中搜索全局最優(yōu)解，使得代價函數(shù)快速收斂。當以下三個收斂條件被滿足時，同化過程結(jié)束，按行政邊界匯總，輸出成熟期預(yù)測結(jié)果(見附圖3)。

s5、用s4所得優(yōu)化參數(shù)數(shù)據(jù)集驅(qū)動wofost作物生長模型，并將tigge(交互式全球大集合預(yù)報系統(tǒng))集合預(yù)報數(shù)據(jù)作為wofost未來時段的氣象輸入，提前10天預(yù)測作物成熟期；

所述步驟s5中處理tigge集合預(yù)報數(shù)據(jù)時采用encu模型，計算方法如下公式(5)所示：

其中，ecmwf、ncep、cma、ukmo分別為歐洲、美國、中國、英國相關(guān)氣象機構(gòu)的預(yù)報數(shù)據(jù)。

s6、逐格網(wǎng)單元運行s5，生成區(qū)域作物成熟期預(yù)測空間分布圖，指導(dǎo)作物收獲。

本發(fā)明實施例所述的遙感、作物生長模型與氣象集合預(yù)報融合的區(qū)域冬小麥成熟期預(yù)測方法，融合了遙感數(shù)據(jù)和作物模型的優(yōu)勢，把遙感觀測和模型模擬的時間序列l(wèi)ai變化趨勢提取出來，通過優(yōu)化算法實現(xiàn)遙感和作物模型的同化，既能在區(qū)域尺度上預(yù)測冬小麥的成熟期，又能提高作物模型的預(yù)測精度，而且氣象預(yù)報數(shù)據(jù)的應(yīng)用又提高了動態(tài)預(yù)測成熟期的精度，實用性更強。

雖然，上文中已經(jīng)用一般性說明及具體實施方案對本發(fā)明作了詳盡的描述，但在本發(fā)明基礎(chǔ)上，可以對之作一些修改或改進，這對本領(lǐng)域技術(shù)人員而言是顯而易見的。因此，在不偏離本發(fā)明精神的基礎(chǔ)上所做的這些修改或改進，均屬于本發(fā)明要求保護的范圍。