基于遙感和地面觀測(cè)數(shù)據(jù)的水稻生長(zhǎng)季長(zhǎng)度空間建模*

胡文君，葉立明※

(1.中國(guó)農(nóng)業(yè)科學(xué)院農(nóng)業(yè)資源與農(nóng)業(yè)區(qū)劃研究所，北京　100081；2.中國(guó)農(nóng)業(yè)科學(xué)院-比利時(shí)根特大學(xué)全球變化與糧食安全聯(lián)合實(shí)驗(yàn)室，北京　100081)

?

·技術(shù)方法·

基于遙感和地面觀測(cè)數(shù)據(jù)的水稻生長(zhǎng)季長(zhǎng)度空間建模*

胡文君1， 2，葉立明1， 2※

(1.中國(guó)農(nóng)業(yè)科學(xué)院農(nóng)業(yè)資源與農(nóng)業(yè)區(qū)劃研究所，北京100081；2.中國(guó)農(nóng)業(yè)科學(xué)院-比利時(shí)根特大學(xué)全球變化與糧食安全聯(lián)合實(shí)驗(yàn)室，北京100081)

農(nóng)作物生長(zhǎng)的季節(jié)周期對(duì)環(huán)境變化敏感，研究農(nóng)作物生長(zhǎng)季長(zhǎng)度的時(shí)空變化規(guī)律對(duì)于農(nóng)業(yè)應(yīng)對(duì)氣候變化、保障國(guó)家糧食安全有重要意義。在區(qū)域尺度上衛(wèi)星遙感是揭示作物生長(zhǎng)季長(zhǎng)度時(shí)空特征的有效手段，但遙感數(shù)據(jù)須經(jīng)地面觀測(cè)校正才能更準(zhǔn)確反映作物的真實(shí)生長(zhǎng)狀況。該研究以黑龍江省水稻生長(zhǎng)的開(kāi)始和結(jié)束日期為研究對(duì)象，通過(guò)在氣象參數(shù)、遙感數(shù)據(jù)和地面站記錄數(shù)據(jù)之間建立轉(zhuǎn)換函數(shù)的方法，將站點(diǎn)觀測(cè)的水稻生長(zhǎng)季數(shù)據(jù)在空間上外推到站點(diǎn)以外區(qū)域，實(shí)現(xiàn)地面觀測(cè)數(shù)據(jù)由點(diǎn)到面的擴(kuò)展。結(jié)果表明：(1)僅包含積溫因子的模型優(yōu)于同時(shí)考慮積溫和降水的模型； (2)利用模型校正后的黑龍江省水稻生長(zhǎng)季長(zhǎng)度，在空間上存在自西南向東北逐漸延長(zhǎng)的趨勢(shì)。該研究為區(qū)域尺度農(nóng)作物生長(zhǎng)季長(zhǎng)度觀測(cè)數(shù)據(jù)的融合提供一種新方法。

生長(zhǎng)季遙感空間建模積溫水稻

0　引言

作為全球變化領(lǐng)域的一個(gè)熱點(diǎn)，農(nóng)作物生長(zhǎng)季研究為認(rèn)識(shí)氣候變化生物響應(yīng)機(jī)理機(jī)制提供重要依據(jù)。在全球變化背景下[1]，農(nóng)作物通過(guò)對(duì)呼吸和光合作用的調(diào)節(jié)來(lái)適應(yīng)氣候變化，而這種調(diào)節(jié)也會(huì)影響陸地—大氣之間的碳循環(huán)，進(jìn)而對(duì)氣候產(chǎn)生反饋[2-5]。農(nóng)作物生長(zhǎng)季的3個(gè)重要指標(biāo)—開(kāi)始期(Start of the Season，SOS)、結(jié)束期(End of the Season，EOS)、生長(zhǎng)季長(zhǎng)度(Length of the Season，LOS)，是作物產(chǎn)量及糧食安全等模型的輸入?yún)?shù)[6-8]，因此生長(zhǎng)季起止日期的確定是作物產(chǎn)量定量分析的基礎(chǔ)[9-10]，對(duì)保障國(guó)家糧食安全具有重要意義[11-13]。

農(nóng)作物的生長(zhǎng)季起止日期可通過(guò)地面觀測(cè)和衛(wèi)星遙感2種方法獲取。傳統(tǒng)的地面觀測(cè)針對(duì)小范圍植株個(gè)體，采用定點(diǎn)定時(shí)記錄方式，具有準(zhǔn)確、精細(xì)和長(zhǎng)時(shí)間序列的特點(diǎn)。然而，除了歐洲的Hubbard Brook和美國(guó)的Harvard Forest等地區(qū)外[14]，全球大部分地區(qū)缺乏有效的地面觀測(cè)數(shù)據(jù)。稀疏的地面站點(diǎn)分布使得大尺度的空間分析難以實(shí)現(xiàn)。隨著遙感技術(shù)的出現(xiàn)與發(fā)展，其覆蓋范圍廣、時(shí)間連續(xù)和成本低廉的特點(diǎn)彌補(bǔ)了地面數(shù)據(jù)的不足，為空間分析的實(shí)現(xiàn)提供了可能[4， 15-17]。但遙感數(shù)據(jù)需要使用地面觀測(cè)數(shù)據(jù)進(jìn)行校正后，才能更準(zhǔn)確的揭示地表植被生長(zhǎng)規(guī)律[18]、補(bǔ)充地面觀測(cè)站點(diǎn)未覆蓋區(qū)域的植被生長(zhǎng)季信息。

目前，用于解決上述問(wèn)題的傳統(tǒng)研究思路可歸納為2種，即“自下而上”法[19-20]和“自上而下”法[21-23]。自下而上法是先確定地面站點(diǎn)的生長(zhǎng)季起止日期，再找出遙感影像上對(duì)應(yīng)區(qū)域的歸一化差值植被指數(shù)(Normalized Difference Vegetation Index，NDVI)，進(jìn)而以NDVI值外推站點(diǎn)外區(qū)域的植被生長(zhǎng)季起止日期； 與之相對(duì)，自上而下法是先確定遙感植被生長(zhǎng)季起止日期，再建立遙感數(shù)據(jù)與地面觀測(cè)數(shù)據(jù)之間的關(guān)聯(lián)，進(jìn)而對(duì)全區(qū)域遙感數(shù)據(jù)進(jìn)行校正。2種校正方法都是僅從一種數(shù)據(jù)源出發(fā)，通過(guò)分析遙感數(shù)據(jù)與地面數(shù)據(jù)之間的相似度來(lái)進(jìn)行偏差校正，未考慮二者之外的獨(dú)立數(shù)據(jù)源，導(dǎo)致校正結(jié)果精度不高、校正方法普適性受限等缺陷。針對(duì)該問(wèn)題，文章在自上而下法研究思路的基礎(chǔ)上，引入氣象因子這一外生變量，并在它和遙感—地面數(shù)據(jù)的差異值之間建立統(tǒng)計(jì)模型，然后將得到的模型應(yīng)用于黑龍江省水稻生長(zhǎng)區(qū)，反演站點(diǎn)外區(qū)域水稻生長(zhǎng)季起止日期，并分析其空間分布特征。該文為區(qū)域尺度單一類型農(nóng)作物生長(zhǎng)季信息提取提供一種新方法。

1　研究區(qū)域與數(shù)據(jù)

1.1研究區(qū)域概況

黑龍江省地處歐亞大陸東部、太平洋西岸、中國(guó)最東北部，介于北緯43°26′～53°33′、東經(jīng)121°11′～135°05′之間，盛行溫帶大陸性季風(fēng)氣候。西北—東南走向的小興安嶺位于該地區(qū)的北部，西北部為大興安嶺，西南部屬松嫩平原，東北部為三江平原。黑龍江省地勢(shì)平坦、土壤肥沃、光水資源優(yōu)越，耕地面積及水稻播種面積均居我國(guó)各省區(qū)首位。近年來(lái)黑龍江省水稻面積擴(kuò)增迅猛， 2004～2011年平均每年增加27萬(wàn)hm2[24]。

圖1　黑龍江省農(nóng)氣觀測(cè)站點(diǎn)和氣象站點(diǎn)分布

1.2數(shù)據(jù)

1.2.1衛(wèi)星遙感數(shù)據(jù)

該文采用2007年中分辨率成像光譜儀(Moderate-resolution Imaging Spectroradiometer，MODIS)全球陸表動(dòng)態(tài)變化產(chǎn)品(MODIS Global land Cover Dynamics Product，MLCD)中的MOD12Q2，該產(chǎn)品是美國(guó)宇航局(National Aeronautics and Space Administration，NASA)陸地觀測(cè)系統(tǒng)生產(chǎn)并發(fā)布的一個(gè)全球物候產(chǎn)品(https：//earthdata.nasa.gov/)。該產(chǎn)品的時(shí)間跨度為2001年至今，每年記錄兩個(gè)生長(zhǎng)周期、每個(gè)生長(zhǎng)周期包含7個(gè)數(shù)據(jù)集，用于監(jiān)測(cè)植被生長(zhǎng)過(guò)程的4個(gè)關(guān)鍵物候期：生長(zhǎng)起點(diǎn)、生長(zhǎng)峰值點(diǎn)、生長(zhǎng)峰值終點(diǎn)和生長(zhǎng)終點(diǎn)[14， 25- 26]。該文使用其中的生長(zhǎng)起點(diǎn)和生長(zhǎng)終點(diǎn)2個(gè)數(shù)據(jù)集。

1.2.2地面觀測(cè)數(shù)據(jù)

地面物候觀測(cè)數(shù)據(jù)采用中國(guó)氣象局資料中心提供的2007年《中國(guó)農(nóng)作物生長(zhǎng)發(fā)育狀況資料數(shù)據(jù)集》，該數(shù)據(jù)詳細(xì)記錄了水稻生長(zhǎng)各個(gè)時(shí)期的發(fā)育程度、生育時(shí)間等，全國(guó)共348站，其中黑龍江省21站(圖1)。

采用2套氣象數(shù)據(jù)。一是《中國(guó)地面氣候資料日值數(shù)據(jù)集》，由中國(guó)氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)提供(http：//www.escience.gov.cn/metdata/page/index.html)，包含發(fā)育期名稱、發(fā)育期日期、生長(zhǎng)狀況等24項(xiàng)指標(biāo)屬點(diǎn)數(shù)據(jù)(圖1)； 二是黑龍江區(qū)域氣溫和降水柵格數(shù)據(jù)，空間分辨率1 km，由國(guó)家氣候中心[27-28]提供。黑龍江省2007年水稻本底數(shù)據(jù)由黑龍江省農(nóng)業(yè)科學(xué)院遙感中心提供。該數(shù)據(jù)以SPOT星上專題制圖儀(Thematic Mapper，TM)所拍攝的影像為數(shù)據(jù)源，通過(guò)專家目視解譯和實(shí)地調(diào)查的方式得到30m分辨率水稻空間分布。這是當(dāng)前可靠性最高的水稻空間分布數(shù)據(jù)集。

2　研究方法

2.1數(shù)據(jù)預(yù)處理

2.1.1水稻生長(zhǎng)季起止日期確定

由于研究區(qū)域水稻均為單季稻，且遙感技術(shù)只能監(jiān)測(cè)到水稻的大田生育期，因此水稻生長(zhǎng)季增強(qiáng)型植被指數(shù)(Enhanced Vegetation Index，EVI)時(shí)序表現(xiàn)為一個(gè)從移栽期至成熟期的單峰曲線[29]。水稻在移栽期前，大田注水，遙感數(shù)據(jù)反映的主要是空曠水面，此時(shí)EVI處于最小值。秧苗移栽到秧苗恢復(fù)生長(zhǎng)，水稻進(jìn)入返青期，EVI此時(shí)開(kāi)始增加，對(duì)應(yīng)于MLCD數(shù)據(jù)中的SOS； 之后兩周，水稻植株分蘗、葉片增多、根系增長(zhǎng)，為生殖生長(zhǎng)積累必要的營(yíng)養(yǎng)物質(zhì)，這期間EVI快速增加。此后水稻經(jīng)歷拔節(jié)、孕穗期，EVI繼續(xù)增加。到抽穗期，營(yíng)養(yǎng)生長(zhǎng)如莖節(jié)伸長(zhǎng)、上位葉生長(zhǎng)和根系發(fā)生仍在進(jìn)行，但水稻營(yíng)養(yǎng)生長(zhǎng)已漸漸達(dá)到頂峰并開(kāi)始生殖生長(zhǎng)，此時(shí)EVI達(dá)到最大值。生殖生長(zhǎng)繼續(xù)進(jìn)行，水稻的養(yǎng)分逐漸轉(zhuǎn)入到籽粒中，EVI逐漸下降。成熟后，植株葉片變黃、谷粒變黃、米質(zhì)變硬，EVI迅速下降，至收獲期EVI達(dá)到最小值[29-34]，對(duì)應(yīng)于MLCD數(shù)據(jù)的EOS。

2.1.2緩沖區(qū)建立與異常值去除

圖2　緩沖區(qū)內(nèi)像元隨緩沖區(qū)大小變化的箱線

地面物候數(shù)據(jù)通常以野外觀測(cè)為基礎(chǔ)，通過(guò)定點(diǎn)目視判別，記錄單個(gè)植株或種群的生長(zhǎng)期節(jié)點(diǎn)[35-37]，而遙感數(shù)據(jù)往往覆蓋范圍廣[36]，以面數(shù)據(jù)的形式進(jìn)行存儲(chǔ)，兩組數(shù)據(jù)存在尺度差異。另外，由于農(nóng)氣站點(diǎn)的設(shè)置并不總是在緊鄰大田的位置，農(nóng)氣站位置并不一定落在水稻像元上。以上2個(gè)問(wèn)題可以通過(guò)建立緩沖區(qū)(buffer)的辦法加以克服。首先，為保證所有站點(diǎn)的緩沖區(qū)內(nèi)均包含水稻像元，需逐步擴(kuò)大緩沖區(qū)范圍； 其次，為避免緩沖區(qū)過(guò)大，無(wú)關(guān)像元干擾增多，使用箱線圖的中位線變化來(lái)確定緩沖區(qū)最大范圍，即利用緩沖區(qū)內(nèi)所有像元值做出的箱線圖，反映緩沖區(qū)內(nèi)水稻生長(zhǎng)季日期的總體信息，隨著緩沖區(qū)范圍的擴(kuò)大，包含的像元數(shù)目增多，箱線圖的中位線位置會(huì)發(fā)生變化，以中位線變化波動(dòng)小于2天來(lái)確定最大的緩沖區(qū)范圍(圖2)。最終，確定最大的緩沖區(qū)范圍為4km。異常值是指一批數(shù)據(jù)中個(gè)別觀測(cè)值明顯偏離其所屬整體的其余觀測(cè)值。異常值的出現(xiàn)可能是由數(shù)據(jù)質(zhì)量問(wèn)題[38]，數(shù)據(jù)錄入錯(cuò)誤引起[39]。若不加判斷和剔除，直接將其代入數(shù)據(jù)的計(jì)算分析中，會(huì)對(duì)結(jié)果產(chǎn)生不良影響[38， 40-41]，因此該研究從統(tǒng)計(jì)圖形角度對(duì)異常值進(jìn)行檢測(cè)和剔除，以提升數(shù)據(jù)質(zhì)量。首先，通過(guò)箱線圖的上(Q3)下(Q1)四分位數(shù)和四分位距(IQR)確定異常值，即小于Q1-1.IQR或大于Q3+1.5IQR的值[42-43]； 其次，基于R語(yǔ)言的箱線圖異常值剔除程序編譯； 最后，數(shù)據(jù)檢查，異常值全部剔除完畢則停止，否則，回到上一步繼續(xù)剔除異常值。

2.2MLCD算法

該文采用的NASA算法[26， 44]處理MLCD數(shù)據(jù)，具體如下。對(duì)已預(yù)處理的時(shí)序數(shù)據(jù)，使用5個(gè)時(shí)相EVI值組成的滑動(dòng)窗口判斷農(nóng)作物生長(zhǎng)的持續(xù)上升和持續(xù)下降階段[44]。由于受到云和氣溶膠等因素的影響，EVI曲線可能有小范圍上升或者下降，而這種波動(dòng)與農(nóng)作物的生長(zhǎng)周期無(wú)關(guān)，為去除無(wú)關(guān)因素的影響，對(duì)持續(xù)上升和下降階段做如下限制：(1)區(qū)間EVI振幅大于年EVI振幅的0.35； (2)區(qū)間極大值大于年內(nèi)極大值的0.7倍。滿足上述兩條件，則此區(qū)間可確定為農(nóng)作物持續(xù)的生長(zhǎng)上升階段或持續(xù)的下降衰老階段。

使用Logistic函數(shù)模擬農(nóng)作物生長(zhǎng)的持續(xù)上升和持續(xù)下降兩階段，Logistic函數(shù)表達(dá)式為：

(1)

式中，t為時(shí)間，y(t)是t時(shí)刻對(duì)應(yīng)的EVI值，a和b是擬合參數(shù)，d是無(wú)云狀況下最小EVI值即EVI初始背景值，c+d是EVI最大值。

Logistic函數(shù)擬合持續(xù)上升和下降區(qū)間后，再利用擬合曲線的曲率變化極值點(diǎn)確定4個(gè)重要的生長(zhǎng)季節(jié)點(diǎn)：(1)生長(zhǎng)起點(diǎn)SOS：EVI增長(zhǎng)階段擬合曲線的曲率第一個(gè)極大值； (2)成熟起點(diǎn)：EVI增長(zhǎng)階段擬合曲線曲率第二個(gè)極大值點(diǎn)； (3)衰落起點(diǎn)：EVI下降階段擬合曲線第一個(gè)極大值點(diǎn)； (4)生長(zhǎng)終點(diǎn)EOS：EVI下降階段擬合曲線第二個(gè)極大值點(diǎn)。

2.3統(tǒng)計(jì)建模

首先，利用站點(diǎn)位置將遙感數(shù)據(jù)與地面觀測(cè)數(shù)據(jù)一一對(duì)應(yīng)，計(jì)算每個(gè)站點(diǎn)遙感數(shù)據(jù)與地面觀測(cè)值之間的差值，并用T檢驗(yàn)檢測(cè)遙感數(shù)據(jù)是否有顯著的提前或者推遲于地面數(shù)據(jù)。差值公式為：

ΔS=SOSMOD-SOSOBS

(2)

ΔE=EOSMOD-EOSOBS

(3)

式中，SOSOBS和EOSOBS分別表示地面物候數(shù)據(jù)記錄的水稻生長(zhǎng)季開(kāi)始日期和結(jié)束日期； SOSMOD和EOSMOD表示對(duì)應(yīng)于地面站點(diǎn)位置的遙感數(shù)據(jù)提取得到的水稻生長(zhǎng)季開(kāi)始日期和結(jié)束日期，以緩沖區(qū)內(nèi)水稻像元均值表達(dá)； Δs為遙感數(shù)據(jù)與地面數(shù)據(jù)計(jì)算所得的生長(zhǎng)季開(kāi)始日期差異值，ΔE為生長(zhǎng)季結(jié)束日期差異值。

其次，選擇積溫和降水兩因子作為10個(gè)候選模型的輸入變量(表1)。由于水稻生長(zhǎng)期間的發(fā)育狀況與氣象條件密切相關(guān)，與其他氣象條件如光照時(shí)長(zhǎng)、風(fēng)速相比，積溫和降水是影響其物候現(xiàn)象發(fā)生的兩個(gè)最重要因子[9， 18， 31， 45-51]。該研究選用水稻生長(zhǎng)期內(nèi)5～9月的累計(jì)降水(P)和大于0℃積溫(T)作為自變量，地面數(shù)據(jù)與遙感數(shù)據(jù)之間的差值作為因變量，建立回歸方程。利用已建立的模型，輸入1km的柵格氣象數(shù)據(jù)和MLCD遙感數(shù)據(jù)，反推無(wú)農(nóng)氣站分布區(qū)域的水稻地面生長(zhǎng)季起止日期，實(shí)現(xiàn)空間上的外推。

2.4統(tǒng)計(jì)檢驗(yàn)

分別對(duì)10個(gè)候選模型進(jìn)行T檢驗(yàn)、F檢驗(yàn)和正態(tài)分布檢驗(yàn)。F檢驗(yàn)用來(lái)檢驗(yàn)自變量與因變量之間線性關(guān)系的顯著性，具體表征指標(biāo)是p值(P-value)和R2。p值越小于0.05，顯著性越強(qiáng)，線性回歸越顯著； R2表示擬合優(yōu)度，R2越大，方程擬合越好。T檢驗(yàn)用來(lái)檢驗(yàn)回歸系數(shù)的顯著性，系數(shù)的顯著性也是通過(guò)p值來(lái)表示，不同自變量對(duì)應(yīng)各自的p值，也就是p值越小(小于0.05就是顯著)則該自變量對(duì)因變量的影響越顯著。除此之外，判斷一個(gè)模型的優(yōu)良程度還需要考察模型的系統(tǒng)性偏差。理想狀況下，表現(xiàn)優(yōu)良的模型殘差是隨機(jī)分布在數(shù)字0兩側(cè)的，并且這些正負(fù)殘差值服從正態(tài)分布，因此需要對(duì)回歸方程殘差的正態(tài)性進(jìn)行檢驗(yàn)，即Shapiro-Wilk檢驗(yàn)。檢驗(yàn)所得的p值越高，則表明模型越優(yōu)； 而當(dāng)p值小于0.1時(shí)，則說(shuō)明殘差不服從正態(tài)分布，模型需要進(jìn)一步改進(jìn)。通過(guò)R2和p值對(duì)模型進(jìn)行比較選取最優(yōu)模型。

表1　模型參數(shù)及檢驗(yàn)結(jié)果

因變量模型自變量參數(shù)估計(jì)回歸系數(shù)顯著性回歸方程顯著性Shapiro-Wilk檢驗(yàn)p-valueR2p-valuep-valueΔEOSY=β0+β1X1X1β0-30.092640.795-0.10040.68370.9497β10.016270.684Y=β0+β1X2X2β00.784230.9520.075820.22380.633β10.019020.224Y=β1X1X1β10.0059434.08E-06***0.90554.09E-06***0.9521Y=β1X2X2β10.0199041.78E-06***0.92131.78E-06***0.5969Y=β0+β1X1+β2X2X1,X2β0120.791970.4390.035920.36510.219β1-0.045420.44β20.033150.188Y=β1X1+β2X2X1,X2β1-5.44E-050.9910.91152.51E-05***0.5889β22.01E-020.241Y=β1X1+β2X12X1β1-4.56E-030.9090.89475.04E-05***0.9486β23.60E-060.793Y=β1X2+β2X22Pβ12.22E-020.1680.91182.48E-05***0.673β2-2.52E-060.879Y=β1X1+β2X2+β3X12X1,X2β13.83E-020.4380.90790.0001544***0.2296β23.34E-020.186β3-1.45E-050.435Y=β1X1+β2X2+β3X12+β4X22X1,X2β1-3.39E-020.7330.90480.000713***0.8701β22.78E-010.357β3-2.40E-060.919β4-1.37E-040.413 注:其中X1和X2分別指水稻生長(zhǎng)期內(nèi)5～9月大于0℃積溫(T)和累計(jì)降水(P);“***”表示極顯著,“**”表示高度顯著,“*”表示顯著,“.”表示不太顯著,沒(méi)有記號(hào)為不顯著

圖3　MLCD和地面觀測(cè)的SOS/EOS散點(diǎn)

3　結(jié)果

3.1最優(yōu)模型選擇

對(duì)大田作物如水稻生長(zhǎng)季的開(kāi)始和結(jié)束日期的觀測(cè)，無(wú)論是人工觀測(cè)記錄還是星載傳感器遙測(cè)，反應(yīng)的都是同一個(gè)地物。天地觀測(cè)在區(qū)域尺度上通常差異，但也可能差異不顯著。因此該文用統(tǒng)計(jì)學(xué)的方法來(lái)檢驗(yàn)差異是否顯著。如圖3，EOS的點(diǎn)均落在1： 1比例線左上方，且對(duì)EOS進(jìn)行T檢驗(yàn)，p值為4.336e-06，說(shuō)明遙感提取的水稻生長(zhǎng)季結(jié)束日期明顯晚于地面農(nóng)氣站觀測(cè)的生長(zhǎng)季結(jié)束日期，天地觀測(cè)的EOS差異顯著，需要引入外生變量來(lái)校正； 而SOS的點(diǎn)均勻分布在1： 1比例線兩側(cè)，即相對(duì)于地面觀測(cè)值而言，遙感方法確定的生長(zhǎng)季開(kāi)始日期并沒(méi)有明顯的提前或者推遲，天地觀測(cè)差異不顯著，意味著將衛(wèi)星觀測(cè)結(jié)果校正到地面觀測(cè)值的需求不存在。因此該研究只對(duì)EOS建立模型，MLCD數(shù)據(jù)的SOS能夠基本反映地面生長(zhǎng)狀況。

對(duì)10個(gè)模型分別進(jìn)行T檢驗(yàn)、F檢驗(yàn)和回歸方程正態(tài)檢驗(yàn)所得的顯著性結(jié)果如表1。對(duì)比10個(gè)模型的檢驗(yàn)結(jié)果可看出，以5～9月積溫為自變量的方程Y=β1X表現(xiàn)最優(yōu)，其回歸系數(shù)和回歸方程的p值均遠(yuǎn)遠(yuǎn)小于0.05，顯著性極強(qiáng)，且R2達(dá)到0.9055，殘差p值為0.9521，各個(gè)指標(biāo)相對(duì)于其他方程都表現(xiàn)出最優(yōu)結(jié)果，故選擇此最優(yōu)回歸方程進(jìn)行水稻生長(zhǎng)季結(jié)束期的空間校正。方程如下：

Y=0.005943X

(4)

其中 X為5月到9月的積溫； Y為MLCD與地面數(shù)據(jù)的水稻EOS差值ΔE。

圖4　模型校正前后星—地觀測(cè)的LOP散點(diǎn)圖對(duì)比

3.2模型精度驗(yàn)證

該文校正SOS和EOS的目的在于得到最接近地面真值的水稻生長(zhǎng)季長(zhǎng)度，為作物生長(zhǎng)模型、糧食安全評(píng)價(jià)服務(wù)。通過(guò)以下步驟對(duì)模型的結(jié)果進(jìn)行驗(yàn)證。

首先，針對(duì)水稻生長(zhǎng)季長(zhǎng)度，該文使用10個(gè)農(nóng)氣站的記錄值對(duì)模型校正后的值進(jìn)行驗(yàn)證。留一交叉驗(yàn)證(Leave-one Out Cross Validation，LOOCV)[52]是該文的主要驗(yàn)證方法，其基本思路是每次從10個(gè)樣本集—10個(gè)農(nóng)氣站記錄數(shù)據(jù)集中留取一個(gè)樣本作為驗(yàn)證集，剩下的9個(gè)樣本作為模型的訓(xùn)練樣本，重復(fù)進(jìn)行10次操作，依次取遍所有的10個(gè)樣本集作為驗(yàn)證樣本，最后將10個(gè)數(shù)據(jù)結(jié)果作為泛化的誤差估計(jì)。

其次，利用農(nóng)氣站觀測(cè)的水稻生長(zhǎng)季長(zhǎng)度與模型校正前后的LOP值對(duì)比，并分別用平均相對(duì)誤差(Relative Rrror，RE)、均方根誤差(Root Mean Square Error，RMSE)和R2、相關(guān)系數(shù)(Correlation coefficient)4個(gè)指標(biāo)進(jìn)行誤差分析和相關(guān)性分析，對(duì)比結(jié)果見(jiàn)表2、圖4。2種顏色的點(diǎn)分別表示模型校正前的LOP和經(jīng)過(guò)最優(yōu)模型校正后的LOP相對(duì)于地面農(nóng)氣站觀測(cè)生長(zhǎng)季長(zhǎng)度的情況，如圖4，校正后的LOP值相對(duì)于校正前LOP值明顯靠近于1： 1比例線，且均勻分布在比例線兩側(cè)； 從表2中可以看出，模型校正后RE由16.9%下降至3.2%，RMSE由19.9下降至4.23，且R2由0.002上升至0.53，相關(guān)系數(shù)由0.05上升至0.73。

根據(jù)模型校正前后的誤差分析以及相關(guān)性比較可知，模型的校正后的值與農(nóng)氣站的地面觀測(cè)值的相關(guān)性更高，誤差更小、更穩(wěn)定，更能反映水稻的真實(shí)生長(zhǎng)季長(zhǎng)度。

表2　模型前后的誤差分析和相關(guān)性分析

RE(%)RMSER2Correlationcoefficient校正后3.24.20.530.73校正前16.919.90.0020.05

3.3水稻生長(zhǎng)季空間分布特征

圖5是MLCD產(chǎn)品直接計(jì)算得到的2007年黑龍江省生長(zhǎng)季開(kāi)始日期空間分布圖。從該區(qū)域總體的空間分布特征來(lái)看，水稻開(kāi)始期集中在第140～157d，即5月下旬，但沒(méi)有明顯的南北差異。

圖6為模型校正的2007年黑龍江省水稻生長(zhǎng)季結(jié)束期空間分布圖。從圖6中可以發(fā)現(xiàn)，黑龍江水稻結(jié)束期存在十分明顯的空間差異，清楚地呈現(xiàn)出一個(gè)從西南至東北逐漸推遲的空間特征，三江平原的水稻結(jié)束期要明顯晚于松嫩平原的水稻結(jié)束期。松嫩平原北部的呼蘭河流域和嫩江上游地區(qū)，水稻生長(zhǎng)季結(jié)束最早，發(fā)生在第239～255d，約9月上旬。三江平原東北部的完達(dá)山兩側(cè)，水稻生長(zhǎng)結(jié)束普遍偏晚，多結(jié)束于第269～274d，約9月下旬，但也有少部分地區(qū)，如佳木斯市的水稻區(qū)，在9月中旬結(jié)束生長(zhǎng)。

結(jié)束期日期減去開(kāi)始期日期所得天數(shù)即為水稻的生長(zhǎng)季時(shí)間長(zhǎng)度，由圖7表示。該圖表明2007年黑龍江水稻生長(zhǎng)季長(zhǎng)度集中在100～120d，且從西南向東北表現(xiàn)出逐漸延長(zhǎng)的趨勢(shì)，三江平原的水稻生長(zhǎng)季長(zhǎng)度明顯長(zhǎng)于松嫩平原的水稻生長(zhǎng)季長(zhǎng)度，三江平原長(zhǎng)度大致118～129d，松嫩平原長(zhǎng)度較短為96～106d。

圖5　利用MLCD數(shù)據(jù)提取的黑龍江省2007年　　　　　　圖6　經(jīng)模型校正的黑龍江省2007年水稻地面水稻生長(zhǎng)季開(kāi)始空間分布　　　　　　　　　　　　　　　　　生長(zhǎng)季結(jié)束空間分布

圖7　2007年黑龍江省水稻地面生長(zhǎng)季長(zhǎng)度空間分布

4　結(jié)論與討論

該文通過(guò)建立氣象參數(shù)、遙感數(shù)據(jù)和地面站點(diǎn)記錄數(shù)據(jù)之間空間模型的方式，外推得到無(wú)農(nóng)氣站分布區(qū)域的水稻地面生長(zhǎng)季起止日期，從而實(shí)現(xiàn)了水稻的地面生長(zhǎng)季數(shù)據(jù)從點(diǎn)到面的尺度轉(zhuǎn)換。主要結(jié)論如下：

(1)分別建立了積溫、降水與遙-地差異值之間的線性回歸模型及非線性回歸模型，經(jīng)過(guò)對(duì)這些模型的比較發(fā)現(xiàn)：①總體上看，積溫和降水都表現(xiàn)出顯著相關(guān)，但積溫顯著性更高； 僅考慮積溫單因子的模型校正效果優(yōu)于同時(shí)考慮積溫和降水的模型校正效果。對(duì)于水稻發(fā)育期而言，降水并不成為主要的發(fā)育期顯著因子，水稻的發(fā)育受溫度影響更顯著，這可能是跟地域差異、環(huán)境條件和管理方式有關(guān)，黑龍江省的水稻主要以灌溉方式為水稻給水，因此對(duì)自然降水的要求相對(duì)偏低。②線性回歸模型較非線性回歸模型擬合優(yōu)度更高。理論上復(fù)雜的線性模型可以通過(guò)增加如具有生理意義的參數(shù)，來(lái)提高模型的描述能力，但在不同的地區(qū)，不同年份，校正效果會(huì)有所不同，另外，隨著模型的復(fù)雜性的增加，模型的一些參數(shù)耐受性會(huì)變低[9]，且擬合會(huì)變得更加困難。

(2)黑龍江省水稻的生長(zhǎng)季結(jié)束期和生長(zhǎng)季長(zhǎng)度存在十分明顯的空間差異。東北部水稻生長(zhǎng)結(jié)束時(shí)間較晚，生長(zhǎng)季長(zhǎng)度較長(zhǎng)，西南部結(jié)束時(shí)間早，生長(zhǎng)季長(zhǎng)度短。造成黑龍江省水稻生長(zhǎng)季結(jié)束期差異性的原因主要與當(dāng)?shù)氐淖匀坏乩憝h(huán)境密切相關(guān)，特別是溫度、降水等氣候條件。緯度越高，溫度越低，農(nóng)作物為保證其正常的發(fā)育成熟，需要通過(guò)延長(zhǎng)生長(zhǎng)時(shí)間來(lái)滿足對(duì)溫度的需求，因此導(dǎo)致了黑龍江東北水稻區(qū)結(jié)束期較晚，生長(zhǎng)季長(zhǎng)度較長(zhǎng)。

該研究為區(qū)域尺度單一類型的農(nóng)作物生長(zhǎng)季信息提取提供了一種新方法，但在模型建立方面，還可對(duì)其進(jìn)行進(jìn)一步優(yōu)化，如加入其他自變量，以提高模型的校正精度。

[1]IPCC.Climate Change 2014 Synthesis Report Summary for Policymakers Chapter.2014

[2]Keenan TF.Phenology：Spring greening in a warming world.Nature， 2015, 526(7571): 48～49

[3]Wu C，Chen JM，Black TA，et al.Interannual variability of net ecosystem productivity in forests is explained by carbon flux phenology in autumn.Global Ecology and Biogeography， 2013, 22(8): 994～1006

[4]Hmimina G，Dufrêne E，Pontailler JY，et al.Evaluation of the potential of MODIS satellite data to predict vegetation phenology in different biomes：An investigation using ground-based NDVI measurements.Remote Sensing of Environment， 2013， 132： 145～158

[5]Morisette JT，Richardson AD，Knapp AK，et al.Tracking the rhythm of the seasons in the face of global change：phenological research in the 21st century.Frontiers in Ecology and the Environment， 2009, 7(5): 253～260

[6]FEWS FSEW.An index to characterize temporal patterns in NDVI at the end of the growing season.Ouagadougou，Burkina Faso, 1993

[7]Nijhuis WAS，Groten SME.TNO/ITC Project proposal to SRON/BCRS for use of NDVI assessment of LGP for climatic modeling of ozone production(LOTUS).1997

[8]李正國(guó)， 唐華俊，楊鵬，等.植被物候特征的遙感提取與農(nóng)業(yè)應(yīng)用綜述.中國(guó)農(nóng)業(yè)資源與區(qū)劃， 2012, 33(5): 20～28

[9]李榮平， 周廣勝，王笑影，等.不同物候模型對(duì)東北地區(qū)作物發(fā)育期模擬對(duì)比分析.氣象與環(huán)境學(xué)報(bào)， 2012, 28(3): 25～30

[10]Ye L，Verdoodt A，Moussadek R，et al.Assessment of China′s food producing capacities using a Web-based land evaluation engine and a grid-based GIS.Advances in Geoecology， 2008， 39： 699～718

[11]Ye L，Van Ranst E.Production scenarios and the effect of soil degradation on long-term food security in China.Global Environmental Change， 2009, 19(4): 464～481

[12]Ye L，Xiong W，Li Z，et al.Climate change impact on China food security in 2050.Agronomy for Sustainable Development， 2013, 33(2): 363～374

[13]Ye L，Tang H，Wu W，et al.Chinese Food Security and Climate Change：Agriculture Futures.Economics：The Open-Access，Open-Assessment E-Journal， 2014， 8(2014-1)： 1～39

[14]Ganguly S，F(xiàn)riedl MA，Tan B，et al.land surface phenology from MODIS：Characterization of the Collection 5 global land cover dynamics product.Remote Sensing of Environment， 2010, 114(8): 1805～1816

[15]Brown ME，de Beurs KM，Marshall M.Global phenological response to climate change in crop areas using satellite remote sensing of vegetation，humidity and temperature over 26 years.Remote Sensing of Environment， 2012， 126： 174～183

[16]Groten SME，Ocatre R.Monitoring the length of the growing season with NOAA.International Journal of Remote Sensing， 2002, 23(14): 2797～2815

[17]Reed BC，Brown JF，VanderZee D，et al.Measuring phenological variability from satellite imagery.Journal of Vegetation Science， 1994, 5(5): 703～714

[18]陳效逑， 胡冰，喻蓉.中國(guó)東部溫帶植被生長(zhǎng)季節(jié)的空間外推估計(jì).生態(tài)學(xué)報(bào)， 2007, 27(1): 65～74

[19]Chen X，Xu C，Tan Z.An analysis of relationships among plant community phenology and seasonal metrics of Normalized Difference Vegetation Index in the northern part of the monsoon region of China.Int J Biometeorol， 2001, 45(4): 170～177

[20]Chen X，Tan Z，Schwartz MD，et al.Determining the growing season of land vegetation on the basis of plant phenology and satellite data in Northern China.Int J Biometeorol， 2000, 44(2): 97～101

[21]You X，Meng J，Zhang M，et al.Remote Sensing Based Detection of Crop Phenology for Agricultural Zones in China Using a New Threshold Method.Remote Sensing， 2013, 5(7): 3190～3211

[22]Badeck F，Bondeau A，Bottcher K.Responses of spring phenology to climate change.New Phytologist， 2004, 162(2): 295～309

[23]White MA，Thornton PE，Running SW.A continental phenology model for monitoring vegetation responses to interannual climatic variability.Global Biogeochemical Cycles， 1997， 11(2)： 217

[24]矯江， 許顯斌，卞景陽(yáng)，等.氣候變暖對(duì)黑龍江省水稻生產(chǎn)影響及對(duì)策研究.自然災(zāi)害學(xué)報(bào)， 2008，(03)： 41～48

[25]馬明國(guó)， 宋怡，王旭峰，等.AVHRR、VEGETATION和MODIS時(shí)間系列遙感數(shù)據(jù)產(chǎn)品現(xiàn)狀與應(yīng)用研究進(jìn)展.遙感技術(shù)與應(yīng)用， 2012, 27(5): 663～670

[26]Zhang X，F(xiàn)riedl MA，Schaaf CB，et al.Monitoring vegetation phenology using MODIS.Remote Sensing of Environment， 2003, 84(3): 471～475

[27]Gao X，Shi Y，Zhang D，et al.Climate change in China in the 21st century as simulated by a high resolution regional climate model.Chinese Science Bulletin， 2012, 57(10): 1188～1195

[28]Gao X，Shi Y，Giorgi F.A high resolution simulation of climate change over China.Science China Earth Sciences， 2011, 54(3): 462～472

[29]于成龍， 劉丹，張志國(guó).基于FY-3的黑龍江省水稻關(guān)鍵發(fā)育期識(shí)別.中國(guó)農(nóng)學(xué)通報(bào)， 2014, 9(30): 55～60

[30]Sun H，Huang J，Peng D.Detecting major growth stages of paddy rice using MODIS data.Journal of Remote Sensing， 2009, 6(13): 1122～1129

[31]官春云. 現(xiàn)代作物栽培學(xué).北京：高等教育出版社， 2011

[32]Li Z，Tang H，Yang P，et al.Spatio-temporal responses of cropland phenophases to climate change in Northeast China.JOURNAl of GEOGRAPHICAL SCIENCES， 2012, 22(1): 29～45

[33]李正國(guó)， 唐華俊，楊鵬，等.東北三省耕地物候期對(duì)熱量資源變化的響應(yīng).地理學(xué)報(bào)， 2011, 66(7): 928～939

[34]Sakamoto T，Yokozawa M，Toritani H，et al.A crop phenology detection method using time-series MODIS data.Remote Sensing of Environment， 2005， 96(3-4)： 366～374

[35]武永峰， 李茂松，宋吉青.植物物候遙感監(jiān)測(cè)研究進(jìn)展.氣象與環(huán)境學(xué)報(bào)， 2008, 24(3): 51～58

[36]Chuanfu X，Jing L，Qinhuo L.Review of advances in vegetation phenology monitoring by remote sensing.Journal of Remote Sensing， 2013, 1(17): 1～8

[37]宛敏渭， 劉秀珍.中國(guó)物候觀測(cè)方法.北京：科學(xué)出版社， 1987

[38]程開(kāi)明. 統(tǒng)計(jì)數(shù)據(jù)預(yù)處理的理論與方法述評(píng).統(tǒng)計(jì)與信息論壇， 2007,(06)： 98～103

[39]王懷亮. 統(tǒng)計(jì)數(shù)據(jù)異常值的識(shí)別及R語(yǔ)言實(shí)現(xiàn).電子技術(shù)， 2012,(05)： 6～7

[40]薛毅， 陳立萍.統(tǒng)計(jì)建模與R軟件.北京：清華大學(xué)出版社， 2009

[41]缺失數(shù)據(jù)統(tǒng)計(jì)分析(中文版). 北京：中國(guó)統(tǒng)計(jì)出版社， 2004

[42]Ye L，Tang H，Zhu J，et al.Spatial patterns and effects of soil organic carbon on grain productivity assessment in China.Soil Use and Management， 2008, 24(1): 80～91

[43]Yao Y，Ye L，Tang H，et al.Cropland soil organic matter content change in Northeast China， 1985-2005.Open Geosciences， 2015, 7(1): 234～243

[44]Zhang X，F(xiàn)riedl MA，Schaaf CB.Global vegetation phenology from Moderate Resolution Imaging Spectroradiometer(MODIS)：Evaluation of global patterns and comparison with in situ measurements.Journal of Geophysical Research， 2006

[45]李正國(guó)， 楊鵬，唐華俊，等.氣候變化背景下東北三省主要作物典型物候期變化趨勢(shì)分析.中國(guó)農(nóng)業(yè)科學(xué)， 2011, 44(20): 4180～4189

[46]齊尚紅， 王冰潔，武作書(shū).農(nóng)業(yè)生產(chǎn)與溫度的關(guān)系.河南科技學(xué)院學(xué)報(bào)(自然科學(xué)版)， 2007,(04)： 20～23

[47]Chen X，Pan W.Relationships among phenological growing season，time-integrated normalized difference vegetation index and climate forcing in the temperate region of eastern China.International Journal of Climatology， 2002, 22(14): 1781～1792

[48]黃淑娥， 田俊，吳慧峻.江西省雙季水稻生長(zhǎng)季氣候適宜度評(píng)價(jià)分析.中國(guó)農(nóng)業(yè)氣象， 2012, 33(4): 527～533

[50]Lambers H，Stuart CIF，Pons TL.Plant Physiological Ecology.Springer， 2008

[51]徐雨晴， 陸佩玲，于強(qiáng).氣候變化對(duì)我國(guó)刺槐紫丁香始花期的影響.北京林業(yè)大學(xué)學(xué)報(bào)， 2004, 6(26): 94～97

[52]Geisser S.A predictive approach to the random effect model.Biometrika， 1974, 61(1): 101～107

SPATIAL MODELING OF THE LENGTH OF GROWING SEASON OF RICE BASED ON REMOTE SENSING AND GROUND OBSERVATION*

Hu Wenjun1，2，Ye Liming1，2※

(1.Institute of Agricultural Resources and Regional Planning，Chinese Academy of Agricultural Sciences，Beijing 100081；2.CAAS-UGent Joint Laboratory of Global Change and Food Security，Beijing 100081)

The seasonal dynamics of crop growth are sensitive to environmental change. The characterization of spatial-temporal patterns of the crop growing season is essential for climate change adaptation and food security improvements. Although satellite remote sensing is an effective means of detecting growing season changes at regional scales, satellite data has to be calibrated before being used. Here we spatially extrapolate the start and end of the growing season of rice in Heilongjiang province of Northeast China by establishing a transfer function between satellite-ground observed difference of growing season dates and climatic parameters. Our results show that: (1) Accumulative temperature is a better parameter, than rainfall, in transfer function establishment; (2) length of the rice growing season tends to increase from southwest to northeast in Heilongjiang province. This paper provides a new method of data fusion in growing season characterization at regional scales.

growing season; remote sensing; special model; accumulated temperature; rice

10.7621/cjarrp.1005-9121.20160802

2015-11-30

胡文君(1991—)，女，湖南衡陽(yáng)人，碩士研究生。研究方向：農(nóng)業(yè)遙感。 ※通訊作者：葉立明(1966—)，男，湖北十堰人，副研究員。研究方向：氣候變化、糧食安全。Email：yeliming@caas.cn

S511; S162.5; S127

A

1005-9121[2016]08-0012-09

*資助項(xiàng)目：中央級(jí)公益性科研院所基本科研業(yè)務(wù)費(fèi)專項(xiàng)“CAAS-UGent全球變化與糧食安全聯(lián)合實(shí)驗(yàn)室科技創(chuàng)新能力建設(shè)”(IARRP-2015-28)