基于遙感信息和WOFOST模型參數(shù)同化的冬小麥單產(chǎn)估算方法研究

陳艷玲, 顧曉鶴, 宮阿都, 胡圣武

(1.北京師范大學(xué)環(huán)境遙感與數(shù)字城市北京市重點實驗室，北京 100875；2.北京師范大學(xué)環(huán)境演變與自然災(zāi)害教育部重點實驗室，北京 100875；3.北京師范大學(xué)地理科學(xué)學(xué)部，北京 100875; 4.國家農(nóng)業(yè)信息化工程技術(shù)研究中心，北京100097; 5.河南理工大學(xué)測繪與國土信息工程學(xué)院，河南焦作 454000)

中國是農(nóng)業(yè)大國，冬小麥?zhǔn)俏覈饕募Z食作物之一，其播種面積占我國糧食總播種面積的1/5[1]。冬小麥長勢監(jiān)測的實時性和產(chǎn)量預(yù)報的準(zhǔn)確性是保障國家糧食安全、經(jīng)濟發(fā)展的基礎(chǔ)。因此，農(nóng)作物長勢監(jiān)測和估產(chǎn)研究已成為國內(nèi)外學(xué)者的研究熱點。迄今為止，衛(wèi)星遙感和作物生長模擬兩項新技術(shù)發(fā)展最為迅速，潛在應(yīng)用價值巨大，同時存在一些有待深入解決和研究的問題。

衛(wèi)星遙感能及時、準(zhǔn)確和客觀地反映農(nóng)作物在時間和空間上的分布信息，在大尺度作物長勢動態(tài)監(jiān)測和估產(chǎn)方面優(yōu)勢明顯。遙感估產(chǎn)研究已從試驗階段逐步進入到實際應(yīng)用階段。但是，地學(xué)遙感信息表征的只是地表或作物冠層的瞬時物理狀況，尚不能真正揭示作物生長發(fā)育和產(chǎn)量形成的內(nèi)在機理過程、個體生長發(fā)育狀況及其環(huán)境氣象條件對其的影響[2]。

近年來，作物生長模型模擬技術(shù)發(fā)展迅速，它綜合考慮外部條件、作物遺傳特性等因素的影響，通過給定一系列作物參數(shù)和環(huán)境參數(shù)，定量模擬作物的生長、發(fā)育和產(chǎn)量形成過程。迄今為止，以荷蘭和美國為代表所建立的作物生長模型多達百種。但是，模型中一些參數(shù)獲取困難，在外界因素的影響下，模擬作物生長過程的一些重要變量地域變異大，難以保證區(qū)域尺度上的模擬精度[3]。

綜上所述，遙感監(jiān)測作物宏觀狀況，表征環(huán)境因素對作物的影響，空間連續(xù)性強，而作物模型模擬作物生長過程，揭示原因和本質(zhì)，具有時間連續(xù)性的優(yōu)勢，兩者各有優(yōu)點和不足。近年來，將遙感與作物模型兩者優(yōu)勢相結(jié)合的研究已成為國內(nèi)外的研究熱點[4-9]，通過遙感數(shù)據(jù)的輔助不僅可以實時準(zhǔn)確地獲取區(qū)域尺度的模型初始輸入?yún)?shù)，同時有利于提高作物生長模型的模擬精度。目前，研究方法主要分為驅(qū)動法和同化法。驅(qū)動法是初期研究的主要方法，使用簡單，但是對反演精度要求高。同化法受遙感數(shù)據(jù)觀測時間、空間以及精度的影響較少，應(yīng)用越來越廣泛。同化方法又可分為兩種，一種是直接利用遙感數(shù)據(jù)作為模型的初始化參數(shù)，另一種是以植被指數(shù)或者反射率等作為同化變量。如Dente[6]通過調(diào)整小麥播種時間、萎蔫點和田間持水量三個參數(shù)，使利用ENVISAT ASAR和MERIS數(shù)據(jù)反演LAI和模擬LAI差異最小，實現(xiàn)雷達數(shù)據(jù)和作物模型的同化，并使同化后的CERES-Wheat模型預(yù)測小麥產(chǎn)量的能力顯著提高。Ma等[10]利用SCE-UA算法將MODIS-LAI數(shù)據(jù)產(chǎn)品同化進WOFOST模型中，通過MODIS-LAI值和模型模擬LAI值之間的差異，調(diào)整WOFOST模型中的關(guān)鍵參數(shù)，取得了較好的產(chǎn)量預(yù)測結(jié)果。Anderson等[11]分析了2003-2013年巴西主要作物不同生育時期的蒸發(fā)脅迫指數(shù)(ESI)與產(chǎn)量的關(guān)系，發(fā)現(xiàn)開花期和灌漿期的ESI與產(chǎn)量相關(guān)性最高。

WOFOST模型是由荷蘭瓦赫寧根農(nóng)業(yè)大學(xué)和世界糧食研究中心共同開發(fā)的具有較大影響力的作物生長模型之一，普適性強，是模擬特定氣候和土壤條件下一年生作物生長的動態(tài)解釋性模型，在農(nóng)作物生長評估、產(chǎn)量預(yù)測、土地定量評價等領(lǐng)域得到廣泛應(yīng)用[12]。本研究利用查找表優(yōu)化算法實現(xiàn)WOFOST模型和遙感反演LAI信息同化，校正模型敏感性參數(shù)，以期提高WOFOST模型區(qū)域尺度小麥單產(chǎn)模擬精度。

1 材料與方法

1.1 研究區(qū)概況

本試驗選取河北省藁城市作為研究區(qū)，該區(qū)地處河北省西南部，介于東經(jīng)114°38′45″～114°58′47″ 、北緯37°51′00″～38°18′44″之間。全境屬于暖溫帶半濕潤大陸性季風(fēng)氣候，四季分明，冬冷夏熱，平均氣溫12.5 ℃。年平均降水量494 mm，7-8月份降水量最多，約占全年的56.2%。年日照時數(shù)2 711.4 h，日照率高達61.2%，無霜期190 d。褐土類和潮土類是藁城市面積最大、最主要的土壤類型。藁城市地處華北平原，地勢平坦，氣候因素分布均勻，是冬小麥的主產(chǎn)區(qū)，因此選取該區(qū)作為研究區(qū)。

圖1 研究區(qū)概況及樣點分布

1.2 數(shù)據(jù)采集

1.2.1 氣象數(shù)據(jù)采集

本研究所需要的2009-2010年的氣象數(shù)據(jù)由中國氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)下載得到，主要包括每日的最高溫、最低溫、降水量、風(fēng)速、日照時數(shù)、水汽壓等氣象要素。WOFOST模型的輸入?yún)?shù)是總輻射量，因此利用FAO公式(1)將日照時數(shù)轉(zhuǎn)化為總輻射量(Rs)。

Rs=(as+bsn/N)Ra

(1)

式中，as和bs是由聯(lián)合國糧食及農(nóng)業(yè)組織FAO為不同地區(qū)提供的與大氣質(zhì)量有關(guān)的參數(shù)，由于本研究區(qū)位于溫帶地區(qū)，因此分別采用0.18和0.55；n為日照時長，可由石家莊氣象站獲取；N為可日照時數(shù)[公式(2)]；Ra為大氣上屆入射輻射[公式(3)]。

N=24ωs/π

(2)

Ra=37.6dr(ωssinφsinδ+cosφcosδsinωs)

(3)

式中，dr為日地距離系數(shù)；ωs為日面中心的時角；δ為地球赤道平面與太陽和地球中心的連線之間的夾角(太陽赤緯)；φ為測點維度，分別由公式(4)計算得到:

dr=1+0.033cos(0.017 2J)

δ=0.420 9sin(0.017 2J-1.39)

(4)

ωs=arcos(-tanφtanδ)

式中，J為日序。

最后將處理好的氣象要素按照WOFOST作物生長模型需要的格式建立對應(yīng)的數(shù)據(jù)庫文件。

1.2.2 農(nóng)學(xué)數(shù)據(jù)采集

2009-2010年在河北省藁城市進行冬小麥主要生育期的田間試驗，選取22個樣本點，樣本點的分布如圖1所示，田間獲取葉面積指數(shù)(LAI)和產(chǎn)量數(shù)據(jù)。

設(shè)定取樣面積為2 500 cm2(50 cm×50 cm)，采用比葉重法測定LAI，計算公式為：

LAI=(w1+w2+w3)/(w1×A×10 000)×S

式中，w1為標(biāo)葉干重(g)；w2為15株的余葉重；S表示標(biāo)葉面積(cm2)；A為取樣面積(cm2)；w3為剩余葉片干重。

在成熟期，收取每個樣本小區(qū)1 m2面積冬小麥，曬干后脫粒測產(chǎn)以獲取每個樣本點的產(chǎn)量數(shù)據(jù)。

1.2.3 遙感數(shù)據(jù)采集

2008年6月我國為環(huán)境監(jiān)測與災(zāi)害損失監(jiān)測發(fā)射了兩顆環(huán)境減災(zāi)小衛(wèi)星HJ-1A和HJ-1B。相比于高時間分辨率的MODIS數(shù)據(jù)和高空間分辨率的Landsat TM數(shù)據(jù)，CCD數(shù)據(jù)兼具高時間分辨率(2 d)和高空間分辨率(30 m)特點，可滿足農(nóng)情監(jiān)測多時相、大范圍的需求，有利于提高大尺度地表狀況動態(tài)監(jiān)測的準(zhǔn)確性。所以，本研究選取環(huán)境減災(zāi)小衛(wèi)星CCD數(shù)據(jù)作為研究所需遙感數(shù)據(jù)，分別選取時相為2010年5月20日和5月29日的兩景晴空過境影像，裁剪至研究區(qū)范圍，并進行輻射定標(biāo)、大氣校正、幾何配準(zhǔn)等處理。

1.3 研究方法

1.3.1 采用植被指數(shù)法反演葉面積指數(shù)(LAI)

在反映作物長勢的指標(biāo)中，LAI是最常用的參數(shù)之一，是衡量植被生長和發(fā)育狀態(tài)的重要因子。由于在可見光和近紅外波段的植被光譜反射率與LAI相關(guān)性高，因而利用LAI與作物植被指數(shù)間的關(guān)系反演LAI的方法已在區(qū)域和大尺度遙感反演中得到廣泛應(yīng)用。本研究基于寬波段反射率，選取了7個在反演冬小麥LAI效果較好的植被指數(shù)[13-18]推算LAI，再分別與實測數(shù)據(jù)進行回歸統(tǒng)計，篩選反演指示性、穩(wěn)定性和預(yù)測能力都較好的植被指數(shù)。

本研究從2010年河北省藁城市冬小麥開花期(5月20日)和灌漿期(5月29日)植被指數(shù)(表1)的22個地面實測點數(shù)據(jù)中，隨機選擇15個實測LAI與對應(yīng)點的植被指數(shù)構(gòu)建統(tǒng)計模型，用剩余的7個數(shù)據(jù)進行驗證。以決定系數(shù)(r2)和均方根誤差(RMSE)來評價模型的準(zhǔn)確性，篩選出與LAI擬合度最好的模型。

表1 冬小麥 LAI 反演各類植被指數(shù)公式Table 1 Formulas of the vegetation index for retrieving LAI of winter wheat

1.3.2 WOFOST模型參數(shù)敏感性分析

模型的“區(qū)域化”和“本地化”是遙感信息和作物模型耦合的首要問題，最關(guān)鍵的環(huán)節(jié)是模型參數(shù)的優(yōu)化和校正，模型參數(shù)的敏感性分析是解決該問題的有效和常用方法，其中擴展傅立葉振幅靈敏度檢驗法EFAST(extended fourier amplitude sensitivity test)具有較大的應(yīng)用潛力[19-20]。為解決WOFOST模型在我國的適用性問題，馬玉平等選取泰安、鄭州和固城三個站點，對該模型參數(shù)進行調(diào)整，使其更好地模擬華北地區(qū)冬小麥的生長過程[21]；郭建茂等對其也進行了適當(dāng)改進，以適應(yīng)山東禹城冬小麥生長模擬[12]。本研究基于前人的研究成果利用敏感性分析專業(yè)軟件Simlab 2.2-EFAST模塊對WOFOST模型中冬小麥26個作物參數(shù)進行全局敏感性分析。首先，采用均勻分布的方法對26個參數(shù)進行插值，得到5 000×26個樣本數(shù)據(jù)。將采樣的參數(shù)組輸入模型，獲得對應(yīng)的5 000個產(chǎn)量數(shù)據(jù)。最后，利用EFAST方法對各個參數(shù)進行敏感性分析，根據(jù)敏感性指數(shù)的大小篩選出對產(chǎn)量影響較大的敏感因子。

1.3.3 模型同化

查找表優(yōu)化算法(look-up table，LUT)是一種簡單易實現(xiàn)的算法，在模型反演核函數(shù)優(yōu)化等方面得到了廣泛應(yīng)用。因此，本研究選取該算法進行模型同化。葉面積指數(shù)LAI的準(zhǔn)確性直接影響作物長勢好壞和對產(chǎn)量高低的評價，因此選取LAI作為優(yōu)化比較對象，優(yōu)化目標(biāo)函數(shù)定義為：

表2 WOFOST模型參數(shù)及波動范圍Table 2 Range of input variables for WOFOST model

式中，LAIsi和LAIri分別為WOFOST模型模擬和遙感反演LAI值；N為外部同化LAI數(shù)據(jù)的個數(shù)。

本研究將目標(biāo)函數(shù)J作為一個判別條件，當(dāng)J值最小，即LAIsi和LAIri間差異最小時，確定待優(yōu)化參數(shù)的取值為“最佳”，完成同化。

本研究總的技術(shù)流程如圖2所示。

2 結(jié)果與分析

2.1 遙感反演LAI結(jié)果

將小麥開花期和灌漿期的8個植被指數(shù)與LAI進行線性擬合，結(jié)果(圖3)表明，兩個時期EVI的擬合效果均最好。兩個時期擬合方程的決定系數(shù)r2分別為0.913和0.798，RMSE分別為0.410和0.470；利用獨立數(shù)據(jù)進行檢驗，r2分別為0.858和0.861，RMSE分別為0.531和0.428(圖4)，說明利用EVI的反演效果較好。

2.2 WOFOST模型參數(shù)敏感性分析結(jié)果

WOFOST模型中冬小麥26個作物參數(shù)的一階敏感性指數(shù)和總敏感性指數(shù)表現(xiàn)基本一致(圖4)。在選取的26個作物參數(shù)中， TSUM1、SPAN、SLATB1、SLATB2、TSUM1和TSUM2等6個參數(shù)比較敏感，一階敏感性指數(shù)和總敏感性指數(shù)均超過0.10，其余參數(shù)敏感指數(shù)均不足0.10。因此，這6個參數(shù)是本研究中WOFOST模型作物參數(shù)“本地化”時需要優(yōu)先考慮調(diào)整和優(yōu)化，剩余參數(shù)可以認為是對產(chǎn)量形成敏感性小或者不敏感，由實測計算、查閱文獻獲得或者直接利用WOFOST模型默認值進行固定。

圖2 總體流程圖

圖3 開花期(a)和灌漿期(b)的EVI與LAI的建模結(jié)果

2.3 模型同化結(jié)果

將篩選出來的6個敏感性因子插值生成5 000×6個樣本數(shù)據(jù)組，分別代入模型生成5 000個相應(yīng)的LAI模擬值。根據(jù)查找表優(yōu)化算法，從WOFOST模型模擬的LAI值中篩選出與遙感反演的LAI值差異最小的數(shù)據(jù)組。對LAI的模型模擬值與遙感反演值進行擬合(圖6)，開花期和灌漿期的r2分別達到了0.912和0.839，說明模型模擬的LAI與遙感反演的LAI差異較小。再將同化后的模型模擬LAI、遙感反演LAI與實測LAI擬合，同化后的模型兩個時期的模擬精度均高于遙感反演精度(圖6和圖7)，說明加入遙感信息后WOFOST模型模擬的LAI更能反映冬小麥生長的真實狀況。

根據(jù)上述確定的敏感性參數(shù)的最佳取值，得到WOFOST模型的模擬產(chǎn)量數(shù)據(jù)。將同化前后的模擬產(chǎn)量、真實產(chǎn)量進行對比，模型同化前對產(chǎn)量的模擬精度較低(圖8a)，同化后模擬精度明顯提高，r2由原來的0.707提高至0.914，RMSE由426.74 kg·hm-2降低至253.67 kg·hm-2，表明模型同化效果較好，也說明利用查找表優(yōu)化算法進行遙感反演LAI與WOFOST作物生長模型同化對作物單產(chǎn)模擬的可行性。

圖4 開花期(a)和灌漿期(b)的反演結(jié)果

圖5 WOFOST模型參數(shù)一階和全局敏感性參數(shù)

圖6 開花期(a)和灌漿期(b)模型模擬LAI和遙感反演LAI擬合精度

圖7 開花期(a)和灌漿期(b)模型模擬LAI和實測LAI擬合精度

圖8 開花期(a)和灌漿期(b)遙感反演LAI和實測LAI擬合精度

圖9 同化前(a)和同化后(b)模型模擬產(chǎn)量和實測產(chǎn)量的相關(guān)性

3 討 論

遙感信息與作物模型同化是近年來進行區(qū)域作物產(chǎn)量估算的重要方法之一。但多數(shù)研究僅局限于田間尺度，在區(qū)域或國家尺度的研究較少。本研究以 LAI為變量，借助查找表同化算法將其與WOFOST作物模型進行耦合，通過調(diào)整敏感性參數(shù)，實現(xiàn)WOFOST產(chǎn)量預(yù)測。在所選取的7個植被指數(shù)中，EVI是預(yù)測LAI效果最好的指數(shù)。NDVI是最常用指示植被覆蓋度及植被生長狀況的植被指數(shù)之一。研究表明，大氣會降低紅光和近紅外反射信號的對比度，散射、上行路徑輻射中大氣的貢獻導(dǎo)致紅光信號的增強，而近紅外信號由于散射和水汽吸收等大氣衰減作用而減弱，最終導(dǎo)致NDVI信號的下降和地表植被量被低估。而EVI是一個優(yōu)化的植被指數(shù)，通過加入藍色波段以增強植被信號，矯正土壤背景和氣溶膠散射的影響，對冠層結(jié)構(gòu)變化冠層類型、植被相和冠層結(jié)構(gòu)更加敏感，提高對高生物量區(qū)域的敏感性。

在本研究中，所選的冬小麥26個作物參數(shù)對產(chǎn)量的全局敏感性指數(shù)和一階敏感性指數(shù)分析結(jié)果的趨勢基本一致。全局敏感性分析綜合考慮參數(shù)自身及參數(shù)間的相互作用對模擬結(jié)果的影響，得到的是各參數(shù)變化對模擬結(jié)果影響的總貢獻，相較于局部敏感性分析而言，篩選和提取的敏感參數(shù)更有利于模型的改進和優(yōu)化。在選取的26個作物參數(shù)中6個參數(shù)比較敏感，敏感性指數(shù)均超過0.10，這與馬玉平[21]、郭建茂等[12]的研究結(jié)果高度一致。其中，DVS=0.5時的比葉面積SLATB2是對冬小麥產(chǎn)量形成影響最大的參數(shù)，可解釋作物產(chǎn)量變化方差的49.0%；另外SLATB1的敏感性指數(shù)也達到了0.237。比葉面積SLATB是葉片單面面積與其干重之比，一般與受光照強弱程度有關(guān)，反映作物對不同環(huán)境的適應(yīng)特征。已有研究表明，SLATB與葉片的凈光合速率呈正相關(guān)，SLATB1和SLATB2的敏感性指數(shù)較大，可以變相表征出光合作用較強，有利于LAI的增長和干物質(zhì)的形成，這就解釋了比葉面積對產(chǎn)量形成最為敏感的原因。其次是出苗到開花期的有效積溫TSUM1的敏感指數(shù)，其值為0.271，而TSUM2的改變對產(chǎn)量形成的影響不顯著。這是因為WOFOST模型是利用實際累計溫度與TSUM1的比值表征冬小麥的生長發(fā)育進程，與LAI的動態(tài)變化和干物質(zhì)的分配等密切相關(guān)。另外，葉片衰老指數(shù)SPAN的敏感指數(shù)也達到了0.191，它是計算葉齡的重要參數(shù)，葉齡是表征作物生長發(fā)育進程的外部形態(tài)診斷指標(biāo)，是影響作物光合作用的內(nèi)部因素之一。20 ℃條件下的單葉光能初始利用效率EFFTB3和30 ℃條件下對最大CO2同化率的影響系數(shù)TMPF4也比較敏感，敏感性指數(shù)分別達到了0.126和0.149。這兩個是影響光合速率的重要參數(shù)，光合速率的大小決定生物量的形成，溫度影響葉片中酶的活性，進而影響光合速率的大小，因此這兩個參數(shù)對產(chǎn)量的形成較為敏感。

氣象條件、土壤狀況、田間管理等外部因素以及小麥品種遺傳參數(shù)的變化會引起相應(yīng)的產(chǎn)量變化。通過敏感性分析篩選出這些因素中對產(chǎn)量影響較大的參數(shù)。光照條件、溫度以及光合效率是影響生物量形成的主要因素。上述研究結(jié)果表明，6個敏感性因子中，比葉面積SLATB1和SLATB2表征光照強弱，TSUM1反映溫度的影響，EFFTB3和TMPF4代表光合速率大小。葉片衰老指數(shù)SPAN表征作物的生長發(fā)育進程，影響葉片光合能力。除敏感性參數(shù)的影響之外，WOFOST模型其他參數(shù)的取值大小對估產(chǎn)結(jié)果也有影響。本研究只考慮了作物參數(shù)的影響，未考慮土壤狀況的影響，計劃在下一步工作中完善。

本研究將遙感監(jiān)測信息同化到模型中，通過優(yōu)化算法調(diào)整敏感性參數(shù)，使其更接近遙感監(jiān)測的“客觀”狀態(tài)，WOFOST模型模擬的LAI更能反映冬小麥生長的“真實”狀況，使得模擬產(chǎn)量更接近“真實值”，模型模擬河北省藁城市的小麥產(chǎn)量與實測產(chǎn)量的r2高達0.914。LAI作為中間參數(shù)，同化后的擬合精度高于遙感監(jiān)測精度，說明同化目的達到。為同化方法應(yīng)用于區(qū)域尺度作物單產(chǎn)估測積累了理論基礎(chǔ)。冬小麥的主要生育時期包括返青期、拔節(jié)期、孕穗期、開花期和灌漿期等，本研究只考慮了開花期和灌漿期的影響，值得進一步的研究。

4 結(jié) 論

本研究通過對WOFOST模型參數(shù)進行適應(yīng)性調(diào)整，建立了適合河北省藁城市冬小麥生長模型，將遙感信息和WOFOST模型相結(jié)合，以LAI為狀態(tài)變量，利用遙感信息校準(zhǔn)模型中6個敏感性作物參數(shù)，模擬研究區(qū)冬小麥產(chǎn)量。研究表明，同化后的區(qū)域產(chǎn)量精度相較于未同化有顯著提高，說明基于查找表算法的遙感信息與作物模型的同化在區(qū)域產(chǎn)量預(yù)測方面具有廣闊的應(yīng)用前景。