基于哨兵-2A模擬反射率及其影像的冬小麥收獲指數估算

任建強 張寧丹 劉杏認 吳尚蓉

(1.中國農業科學院農業資源與農業區劃研究所， 北京 100081；2.中國農業科學院農業環境與可持續發展研究所， 北京 100081)

0 引言

收獲指數(Harvest index，HI)是評價作物單產水平和栽培成效的重要生物學參數，也是作物單產進一步提高的重要決定因素之一[1]。對糧食作物來說，一般的收獲指數是指成熟期作物籽粒產量占作物地上生物量的百分數，該指標本質反映了作物同化產物在籽粒和營養器官中的分配比例[2]。近年來，為了定量描述籽粒灌漿過程中作物收獲指數逐步形成至達到最大值的變化過程，部分學者以冬小麥為例提出了作物動態收獲指數(Dynamic harvest index，D-HI)的概念，這不僅豐富了收獲指數的概念內涵，也對作物生長過程和產量形成的定量模擬具有重要意義[3]。因此，快速準確地獲取成熟期作物收獲指數及其動態過程信息對于作物單產準確模擬[4-6]、作物生物量估算[7]、作物品種選育和表型信息獲取[8-9]、作物栽培技術優化與效果評價[10-11]具有重要科學意義，同時，對農業管理部門及時掌握農作物長勢和作物產量估算信息，有效開展農業生產管理也具有重要指導意義[12-14]。

一般的作物收獲指數研究主要基于田塊尺度的農學試驗層面，大多側重于作物收獲指數的數學模擬、收獲指數與相關農學參數關系、作物生長環境及其管理措施對收獲指數的影響評價等[15-16]。作物收獲指數獲取方法主要有直接法和間接法。直接法主要通過田間取樣計算獲得，該方法雖然準確，但費時費力且很難在大范圍內開展，也無法獲得收獲指數的連續空間分布信息。間接法主要包括2類：①通過模擬灌漿期作物收獲指數的逐步增加變化過程與時間之間函數關系實現作物收獲指數的準確估算[17-18]。②基于作物生長過程中的植被信息(如各種植被指數、生物量等)、作物生長環境影響因子(如溫度、光照、土壤含水率等)與HI建立函數關系實現HI的估算[19]。

近些年來，遙感技術憑借其覆蓋范圍大、快速和準確獲取地表作物參數信息的優勢，為準確獲取區域作物HI空間信息提供了可靠的技術手段[20-21]。其中，通過寬波段多光譜傳感器(如NOAA、MODIS、MERIS等)數據，國內外學者利用發芽-開花和開花-成熟2個階段能夠反映作物長勢狀況的時序植被遙感信息(如歸一化差值植被指數和葉面積指數等)開展了成熟期冬小麥HI遙感估算研究，對利用遙感信息獲取區域尺度作物HI具有重要借鑒意義[15,17,22-24]。然而，已有研究只進行了成熟期最終作物收獲指數估算，且使用開花后時序NDVI表征籽粒灌漿過程，對作物HI形成過程中動態收獲指數信息的考慮不足[25]。針對這一問題，張寧丹等[26]以河北省衡水市深州市冬小麥為研究對象，在構建花后累積地上生物量比例動態參數(Dynamic fG, D-fG)基礎上，提出了敏感波段中心構建歸一化差值光譜指數(Normalized difference spectral index，NDSI)估算D-fG的作物動態收獲指數估測方法，實現了田間冠層尺度D-HI的準確獲取。在敏感波段中心基礎上，通過波段擴展確定了冬小麥D-fG估算敏感波段最大寬度，實現了最大波寬下田間冠層尺度D-HI的遙感準確獲取。上述冠層尺度D-HI估算研究為基于高光譜和多光譜等遙感衛星開展區域尺度作物動態收獲指數獲取奠定了一定基礎。但是，該方法由田間冠層尺度擴展到區域尺度的應用效果還有待進一步評價，在已獲得的敏感波段中心和最大波寬基礎上如何開展基于遙感衛星的區域冬小麥收獲指數遙感估算還有待進一步深入研究，特別是如何充分利用反映作物生長狀況的紅光、紅邊、近紅外等波段信息進行大范圍D-HI遙感估算成為亟待研究的問題。

基于以上分析，考慮到哨兵-2A(Sentinel-2A)是具有較高重訪頻率的寬幅高分辨率多光譜成像衛星，是目前唯一在紅邊區域含有3個獨特紅邊波段的光學遙感數據，這對監測植被健康信息和植物生長狀況非常有效[27-29]，而且該數據已在大范圍作物監測中發揮了重要作用[30-33]。但是，基于Sentinel-2數據的作物HI估測研究尚未見報道。因此，為進一步驗證基于花后累積地上生物量比例的D-HI遙感估算方法在大范圍獲取作物HI空間信息的可行性和有效性，在前期田間冠層尺度作物HI估算研究確定的冬小麥D-fG和HI估算敏感波段中心和最大波段寬度基礎上，本研究利用近地高光譜和Sentinel-2A影像數據，開展基于Sentinel-2A遙感模擬數據和真實遙感數據的區域冬小麥HI估算研究，以期為大范圍作物HI空間信息衛星遙感準確獲取提供技術方法借鑒。

1 研究區與數據

1.1 研究區概況

圖1 研究區位置與采樣點空間分布Fig.1 Location of study area and spatial distribution of sampling points

研究區為中國北方糧食生產基地黃淮海平原內河北省衡水市(115.17°～116.57°E，37.05°～38.38°N)，典型試驗區為深州市(115.35°～115.80°E，37.70°～38.17°N)(圖1)。衡水市總面積8.836×107km2，地處河北沖積平原，地勢自西南向東北緩慢傾斜，海拔12～30 m。該研究區域屬于暖溫帶半干旱區季風氣候，年均降水量約480 mm，年均溫度約13.4℃，無霜期200 d左右，研究區為冬小麥-夏玉米一年兩熟輪作種植制度。其中，冬小麥種植時間為10月上中旬，返青期為次年3月上中旬，拔節期為3月下旬至4月上中旬，抽穗開花期為4月下旬至5月上旬，灌漿-乳熟期為5月中下旬，成熟期為6月上旬。

1.2 數據獲取與處理

1.2.1地面數據采集與處理

數據獲取內容主要包括冠層高光譜、作物地上干生物量、動態籽粒產量和GPS定位信息等。本研究分別在冬小麥開花期(2021年5月3日)、灌漿前期(2021年5月15日)、灌漿后期(2021年5月25日)、成熟期(2021年6月5日)共4個關鍵時期進行數據采集。深州市典型研究區共18個樣方，每個樣方布設5個采樣點，每次試驗獲取90個樣本數據。衡水市區域調查共布設50個調查點，每次調查獲取50個樣點數據，主要包括作物地上生物量、動態籽粒產量和GPS定位信息等。以冬小麥開花期為時間基準，在深州市共獲取了灌漿前期、灌漿后期和成熟期共3個時期270個地面樣本數據，根據比例3∶2 將其分為建模數據集(162個)和驗證數據集(108個)。衡水市10個縣50個調查點在灌漿前期、灌漿后期和成熟期3個時期共獲得150個地面樣本數據，主要作為后期衡水地區作物收獲指數的區域驗證數據。

(1)地上干生物量

在冬小麥關鍵生育期(包括開花期、灌漿前期、灌漿后期和成熟期)，首先根據各個采樣點測定的GPS位置信息，在每個采樣點分別割取1行長度20 cm的冬小麥地上部分作為樣本。在將冬小麥莖葉穗分離基礎上，將分離后冬小麥莖葉穗放入105℃干燥箱并進行30 min殺青處理。然后，在85℃條件下連續干燥48 h以上，直至樣本質量恒定再進行稱量。其次，在將冬小麥莖葉穗干質量進行相加的基礎上，根據種植密度和樣本干質量換算成單位面積冬小麥地上干生物量。最后，在對各樣點冬小麥穗脫粒處理基礎上，分別記錄各個采樣點的籽粒質量。

(2)冬小麥冠層高光譜

利用美國ASD(Analytical Spectral Devices Inc.)公司生產的Field Spec 4型光譜輻射儀(波長350～2 500 nm)采集冬小麥冠層高光譜。其中，波長范圍350～1 000 nm采樣間隔1.4 nm，波長范圍1 000～2 500 nm采樣間隔2 nm，重采樣后光譜間隔為1 nm。選擇天氣狀況良好、陽光照射充足條件下進行冠層高光譜測定，觀測時間為10:00—14:00。高光譜測量前用標準白板校正，測量時保證探頭垂直向下，光譜裝置探頭視角為25°，為了保證冬小麥樣本處于探測視場內，同時減少下墊面光譜反射對測定結果的影響，探頭距離作物冠層頂部高度約0.5 m。測量時，每個采樣點獲取10條光譜數據，取其均值作為該采樣點的光譜反射率。然后，對每個采樣點的光譜曲線利用9點加權移動平均法進行去噪平滑處理。圖2為深州市調查樣方采樣點經過光譜平均和光譜平滑處理后的不同生育期冬小麥冠層高光譜曲線。

圖2 深州市典型樣方冬小麥不同生育期冠層 高光譜曲線Fig.2 Canopy hyperspectral curves of winter wheat at different growth stages of quadrats in Shenzhou City

1.2.2遙感數據獲取與預處理

Sentinel-2A數據主要從歐空局哥白尼數據中心(https:∥scihub.copernicus.eu/)獲取，衛星數據產品級別為Level-1C(L1C)。Sentinel-2A遙感影像有13個光譜波段，包括可見光、紅邊、近紅外、水汽、卷云以及短波紅外波段等。其中，L1C級數據為已經過輻射定標和幾何校正的大氣頂反射率數據，因此，只需要進行大氣校正處理便可得到地表反射率數據。利用Sen2cor、SNAP(Sentinel application platform)和ENVI軟件對衛星影像分別進行大氣校正、格式轉換和鑲嵌裁剪等預處理工作。為了保證遙感數據各波段空間分辨率的一致性，在SNAP平臺中使用最鄰近內插法將遙感數據各波段重新采樣至10 m。結合冬小麥物候信息以及地面觀測試驗時間，根據選擇影像云量較低原則(小于10%)，獲得了衡水市2021年5月17日、2021年5月27日、2021年6月6日共3個時期Sentinel-2A影像數據。考慮到前期冠層高光譜的作物收獲指數估算敏感波段中心和最大波段寬度僅對應Sentinel-2A的藍光、綠光、紅光、紅邊-1、紅邊-2、紅邊-3、近紅外等8個波段，因此，僅使用了相應的B2波段至B8a波段信息，具體Sentinel-2A多光譜影像波段信息如表1所示[34-36]。

1.2.3其他輔助數據

(1)光譜響應函數

由于ASD冠層高光譜與Sentinel-2A衛星影像數據的波段寬度不同，為了將前期冠層尺度ASD高光譜確定的冬小麥D-fG估算敏感波段中心以及最大波段寬度應用到寬波段多光譜Sentinel-2A衛星影像數據，需要利用衛星的光譜響應函數將近地面高光譜數據轉換為衛星波段反射率，從而將近地面高光譜數據與Sentinel-2A衛星遙感各個波段反射率對應。Sentinel-2A光譜響應函數如圖3所示。

表1 本研究選取的Sentinel-2A衛星遙感主要 波段參數Tab.1 Main band information of Sentinel-2A selected in this study

圖3 Sentinel-2A光譜響應函數Fig.3 Curves of spectral response function of Sentinel-2A

(2)冬小麥空間分布信息

為了獲得研究區的冬小麥空間分布信息，采用被廣泛使用的監督分類中支持向量機(Support vector machine，SVM)方法進行研究區2021年冬小麥空間分布提取，利用的遙感影像為Sentinel-2A影像。支持向量機是一種建立在統計學基礎上的機器學習方法[37]，即在有限的分類樣本信息前提下，在模型的復雜性和學習能力之間尋求最佳結果，保證得到的極值解是全局的最優解[38]。

樣點數據主要包括訓練樣點數據和驗證樣點數據。其中，訓練樣點用于支持向量機分類進行作物分布提取，驗證樣點用于對冬小麥空間分布結果的精度驗證。根據冬小麥實際調查點和Google Earth上選取的冬小麥和非冬小麥樣點，在衡水研究區最終獲得冬小麥樣點和非冬小麥樣點共計4 845個。其中，選擇冬小麥樣點870個、非冬小麥樣點1 175個作為訓練樣本，其余1 300個冬小麥樣點和1 500個非冬小麥樣點作為驗證樣本對冬小麥空間分布結果進行精度驗證，具體地面樣點數據分布信息如圖4所示。最終，本研究基于2021年4月19日的Sentinel-2A影像和SVM分類方法實現冬小麥空間分布提取。通過驗證，冬小麥空間分布提取總體精度為91.77%，Kappa系數為0.835 6，精度達到較高水平，可以滿足本研究區域冬小麥收獲指數遙感估算所需作物空間分布的精度要求，衡水市冬小麥空間分布結果如圖5所示。

圖4 衡水市2021年冬小麥提取地面樣點數據分布Fig.4 Distribution of ground sampling points for winter wheat mapping in Hengshui City in 2021

圖5 衡水市2021年冬小麥空間分布結果Fig.5 Spatial distribution results of winter wheat in Hengshui City in 2021

2 研究方法

2.1 總體技術路線

在篩選出的敏感波段中心及其最大波段寬度基礎上，首先，確定對應Sentinel-2A遙感數據的波段范圍，并利用遙感數據的光譜響應函數和地面高光譜數據進行寬波段遙感反射率模擬；其次，根據模擬的寬波段反射率數據構建歸一化差值光譜指數(NDSI)進行D-fG和D-HI的估算和精度驗證。在此基礎上，篩選出冬小麥D-fG和D-HI估算最優波段組合。最后，在最優波段組合信息篩選結果下，基于Sentinel-2A真實衛星遙感數據進行區域冬小麥D-fG和D-HI的估算與驗證。主要技術路線如圖6所示。

圖6 技術路線Fig.6 Flowchart of the research

2.2 主要參數構建與計算

2.2.1D-fG參數構建

為了估算作物收獲指數，KEMANIAN等[39]提出了fG參數，該參數為作物開花期-成熟期累積地上生物量與成熟期地上生物量的比值，但一般fG參數只應用于成熟期fG計算，而未考慮開花期至成熟期之間fG參數的動態變化過程，這在一定程度上會降低利用fG參數進行作物收獲指數估算的精度。為了提高作物收獲指數估算模型的精度，增強模型的穩定性和可靠性，降低成熟期fG參數不穩定對收獲指數估算精度的不利影響，在已有研究基礎上，考慮了開花期-成熟期fG參數的動態變化過程，采用動態參數D-fG進行作物動態收獲指數估算，即作物開花期至采樣時期累積的地上生物量與對應采樣時期地上生物量間比值[26]。D-fG指標計算公式為

(1)

式中DfG——花后累積地上生物量比例動態參數(D-fG)指標

Wpost——冬小麥開花期至采樣日期累積地上生物量，kg/hm2

Wwhole——冬小麥播種至采樣日期累積地上生物量，kg/hm2

Wt——播種至采樣日期t累積地上干物質量，kg/hm2

Wa——開花期地上干物質量，kg/hm2

2.2.2動態收獲指數構建

一般糧食作物(如小麥、玉米、水稻等)開花至成熟期間，隨著籽粒不斷灌漿，作物籽粒產量占作物地上部干物質量百分比呈逐步增加狀態，直到成熟期收獲指數達到最大值，其中，用于定量描述作物灌漿過程中籽粒產量占作物地上部干物質量百分比逐步增加至最大值的收獲指數變化過程指標，稱為動態收獲指數[3]。因此，本文在已有研究基礎上，利用冬小麥籽粒灌漿過程中各個地面觀測時間的收獲指數動態變化信息獲得動態收獲指數信息，即在獲得各采樣點冬小麥莖葉穗干質量基礎上，分別對各采樣點小麥穗進行脫粒處理，并記錄各采樣點的籽粒質量。最后，計算冬小麥灌漿至成熟期期間各個地面觀測時間的冬小麥動態收獲指數(D-HI)。D-HI計算公式為

(2)

式中DHI——作物動態收獲指數(D-HI)

WZ,t——灌漿至成熟期期間采樣日期t的冬小麥籽粒干質量，kg/hm2

WA,t——采樣日期t的冬小麥地上干生物量，kg/hm2

2.2.3NDSI計算

光譜指數(Spectral index，SI)可以通過某些特定波段的組合來指示綠色植被內部的色素含量、水分變化和營養狀態等[40]。為更好地利用多光譜遙感各個波長所包含的信息，本文將冠層高光譜模擬的多光譜遙感反射率及其真實遙感反射率分別進行任意兩兩波段組合，從而構建歸一化差值光譜指數(Normalized difference spectral index，NDSI)，構建形式為

(3)

式中NDSI(λ1,λ2)——波長λ1、λ2計算的NDSI指數，值域為[-1,1]

Rλ1——波長λ1所對應光譜反射率

Rλ2——波長λ2所對應光譜反射率

考慮到作物光譜在波段1 350～1 415 nm和波段1 800～1 950 nm受大氣和水蒸氣影響較大[41]，因此，在可見光-近紅外波段350～1 000 nm內進行D-fG估算敏感波段篩選和動態作物收獲指數遙感估算，即式(3)中λ1、λ2表示在350～1 000 nm內的任意波長λ1和波長λ2。

2.3 地面高光譜模擬寬波段遙感數據方法

為了利用模擬遙感數據篩選作物D-fG和D-HI估算的最優波段，需利用灌漿前期、灌漿后期和成熟期的模擬Sentinel-2A多光譜數據開展冬小麥關鍵生育期動態收獲指數估算。在已有近地面冠層高光譜確定的冬小麥D-fG和D-HI估算敏感波段中心和最大波段寬度研究結果基礎上[26]，在對應最大波段寬度范圍內[26](表2)，借助Sentinel-2A光譜響應函數，將實測的近地面高光譜窄波段反射率轉換為Sentinel-2A多光譜波段反射率模擬數據，波段反射率轉換公式為

(4)

式中Rrs——模擬衛星波段反射率

λm——傳感器光譜探測起始波長

λn——傳感器光譜探測終止波長

S(λ)——傳感器在波長λ處的光譜響應函數值

R(λ)——冠層光譜在波長λ處光譜反射率

表2 D-fG和D-HI估算敏感波段中心最大波寬對應的 地面高光譜波段范圍Tab.2 Wavelength range of ground hyperspectral bands corresponding to center of sensitive band and its maximum bandwidth for estimation of D-fG and D-HI

2.4 歸一化光譜指數NDSI與D-fG間模型構建

基于模擬Sentinel-2A反射率數據和真實Sentinel-2A衛星反射率數據構建的歸一化光譜指數NDSI分別進行D-fG遙感估算，為冬小麥動態收獲指數獲取奠定基礎。NDSI與D-fG間的線性模型為

DfG=aNDSI(λ1,λ2)+b

(5)

式中a——一次項系數b——常數項

當利用模擬的寬波段反射率數據進行D-fG估算時，λ1、λ2分別為模擬多光譜遙感數據的波長，NDSI為模擬多光譜遙感波段λ1、λ2構建的NDSI光譜指數。當基于多光譜遙感衛星數據進行區域D-fG估算時，NDSI為真實寬波段遙感數據構建的光譜指數。

2.5 冬小麥動態收獲指數估算模型

在KEMANIAN等[39]提出的基于成熟期實測fG的作物收獲指數估算方法基礎上，提出了基于D-fG遙感信息的D-HI遙感估算方法。D-fG和動態收獲指數D-HI間統計關系模型為

DHI=HI0+kDfG

(6)

式中HI0——截距，即在作物開花期之后生物量不發生變化情況下的動態收獲指數

k——斜率常數

2.6 精度評價

精度評價指標包括決定系數(Coefficient of determination，R2)、均方根誤差(Root mean square error，RMSE)、歸一化均方根誤差(Normalized root mean square error，NRMSE)和平均相對誤差(Mean relative error，MRE)。其中，基于模擬Sentinel-2A的冬小麥收獲指數估算中，主要使用深州市灌漿前期、灌漿后期和成熟期共3個時期270個地面樣本數據，將其分為建模數據集(n=162)和驗證數據集(n=108)。基于真實Sentinel-2A衛星遙感波段的冬小麥收獲指數的估算中，共獲取2021年5月17日、2021年5月27日、2021年6月6日3個時期遙感影像，綜合使用深州市270個地面樣本數據和衡水區域調查的150個驗證數據，最終得到建模數據集(n=162)和驗證數據集(n=258)。

3 結果與分析

3.1 遙感數據間波段對應關系及NDSI計算

3.1.1Sentinel-2A模擬數據波段范圍確定

在綜合考慮敏感波段中心最大寬度對應的地面高光譜波段范圍和Sentinel-2A真實遙感數據波段范圍基礎上，進一步確定了每個敏感波段中心最大波段寬度所對應的模擬Sentinel-2A波段范圍。由于Sentinel-2A衛星遙感寬波段數據的波段寬度比地面高光譜波段數據寬，因此，構建NDSI的兩個敏感波段在同一個寬波段內的可能性比較大，篩選出的波段組合相對較少。考慮到篩選出的地面高光譜波段范圍與真實Sentinel-2A波段范圍并非完全一致，為了便于利用光譜響應函數模擬Sentinel-2A波段反射率，當地面高光譜波段范圍超出真實Sentinel-2A波段范圍時，應用Sentinel-2A波段范圍內的數據；當地面高光譜波段在Sentinel-2A波段范圍內時，應用地面高光譜波段內數據[42]。地面高光譜波段范圍與Sentinel-2A衛星遙感及其模擬數據之間波段對應關系如表3所示。其中，當D-fG敏感波段中心在(366 nm，489 nm)時，對應的地面高光譜波段范圍(351～381 nm)不在Sentinel-2A真實波段范圍內，因此，無法構建相應的NDSI；同理，當D-fG敏感波段中心在(443 nm，495 nm)時，地面高光譜波段(409～477 nm和461～529 nm)所對應的Sentinel-2A真實波段均處于藍光波段內，也無法構建相應的NDSI，故本研究對上述兩種情況進行了舍棄處理。

表3 地面高光譜波段范圍與Sentinel-2A衛星遙感及其模擬數據之間對應關系Tab.3 Correspondence between wavelength ranges of the ground hyperspectrum and bands of Sentinel-2A satellite remote sensing images and their simulated date

3.1.2NDSI計算

根據表3中高光譜敏感波段對應的模擬Sentinel-2A波段數據和真實Sentinel-2A波段數據，計算了冬小麥灌漿前期、灌漿后期、成熟期Sentinel-2A模擬數據及其真實遙感影像各個波段組合構建的NDSI。限于篇幅，本文僅展示地面高光譜波段λ1(672～680 nm)和λ2(855～875 nm)對應Sentinel-2A紅光波段和窄近紅外波段構建的NDSI的空間分布圖，冬小麥灌漿前期、灌漿后期和成熟期NDSI空間信息計算結果如圖7所示。通過分析可知，冬小麥灌漿前期NDSI范圍在0.40～0.95之間，研究區冬小麥NDSI平均值為0.88；冬小麥灌漿后期NDSI范圍在0.20～0.80之間，研究區冬小麥NDSI平均值為0.58；冬小麥成熟期NDSI范圍在0.10～0.50之間，研究區冬小麥NDSI平均值為0.35。

圖7 基于Sentinel-2A衛星紅光和窄近紅外波段構建的冬小麥NDSI空間分布圖Fig.7 NDSI spatial distribution map constructed by red band and NIR2 band of Sentinel-2A

3.2 動態收獲指數估算模型建立

根據深州市建模數據集(n=162)中的地上生物量數據和灌漿過程中籽粒產量動態數據，計算162個冬小麥樣本點的D-fG和動態收獲指數D-HI。在此基礎上，對實測D-fG和實測動態收獲指數D-HI間的線性關系進行擬合，得到D-fG參數和動態收獲指數D-HI間擬合方程

DHI=0.104 3+0.763 1DfG

(7)

其中，基于實測D-fG的動態收獲指數D-HI線性估算模型R2為0.940 9(圖8)，這為開展基于D-fG遙感參數信息的動態收獲指數遙感估算奠定了基礎。

圖8 基于D-fG參數的D-HI估算模型Fig.8 D-HI estimation model based on D-fG parameter

3.3 基于模擬Sentinel-2A波段數據的冬小麥收獲指數估算

3.3.1冬小麥D-fG參數估算和精度驗證

利用Sentinel-2A光譜響應函數模擬地面高光譜波段所對應的衛星遙感反射率數據，并以模擬反射率結果構建相應的NDSI；然后，構建模擬NDSI與實測D-fG之間的估算模型，并估算冬小麥D-fG。最后，得到基于模擬Sentinel-2A波段構建NDSI與D-fG間統計關系及D-fG估算精度，結果如表4所示。從表4可以看出，9個模擬Sentinel-2A波段構建的NDSI擬合D-fG在P<0.01水平上均達到極顯著水平，估算模型R2在0.663 6～0.836 2之間。通過驗證數據集(n=108)進行精度檢驗可知，其精度均達到了較高水平。其中，精度驗證R2在0.908 9～0.925 9之間，RMSE在0.038 6～0.048 2之間，NRMSE在11.18%～13.96%之間，MRE在10.07%～12.31%之間。可以看出，相較于紅、綠、藍波段的組合，紅光波段、紅邊波段與近紅外波段間組合的HI估算精度普遍較高，這主要由于紅邊波段是指示綠色植物生長狀況的敏感波段，能夠有效反映植被養分狀況、健康狀態和生理生化參數等信息。其中，λ1(672～680 nm)和λ2(855～875 nm)模擬Sentinel-2A反射率構建的NDSI進行冬小麥D-fG估算的精度最高，其RMSE、NRMSE和MRE分別為0.038 6、11.18%、10.07%，這主要是由于近紅外與紅光波段交界處快速變化的區域能夠對植被冠層結構和葉綠素含量等微小變化和植被生長狀況進行有效反映，同時，紅光和近紅外波段反射率有明顯的反差。因此，紅光和處于近紅外區域的窄近紅外波段組合構成的NDSI估算D-fG精度更高，效果更為理想。

3.3.2冬小麥動態收獲指數估算和精度驗證

在模擬波段構建的NDSI估算D-fG的條件下，根據式(6)分別計算冬小麥D-HI估算結果，并進行精度驗證。冬小麥動態收獲指數總體精度驗證結果如表5所示。從表5可知，估算結果驗證均達到了高精度水平，其精度驗證R2在0.889 9～0.909 7之間，RMSE在0.040 4～0.051 5之間，NRMSE在10.83%～13.81%之間，MRE在9.56%～12.38%之間。其中，紅光波段、紅邊波段與近紅外波段間組合進行D-HI估測精度均較高，且基于波段λ1(672～680 nm)和λ2(855～875 nm)模擬Sentinel-2A反射率數據估算D-fG參數的D-HI估測結果精度最高，RMSE、NRMSE和MRE分別為0.040 4、10.83%、9.56%。由于紅光和近紅外波段的反射率具有明顯的反差，故由上述紅光波段和窄近紅外波段構成的NDSI對不同植被光譜的變化更為敏感，能很好地反映作物長勢、生長狀況以及D-fG，加之D-fG和D-HI間具有較好的正相關關系，因此，基于波段λ1(672～680 nm)和λ2(855～875 nm)獲得的D-HI估算精度最高。上述基于模擬Sentinel-2A遙感數據波段組合進行D-HI估測的精度結果，對基于真實遙感數據的D-HI估測中遙感數據源選擇、波段組合優選具有重要指導意義。

表4 基于模擬Sentinel-2A波段構建NDSI的D-fG估算精度驗證Tab.4 Verification of estimation accuracy of D-fG based on NDSI constructed by simulated bands of Sentinel-2A

表5 基于模擬Sentinel-2A波段的D-HI估算模型總體 精度驗證Tab.5 Verification of overall accuracy of D-HI estimation model based on simulated bands of Sentinel-2A

3.4 基于Sentinel-2A衛星遙感數據的冬小麥收獲指數估算

3.4.1冬小麥D-fG空間信息估算和精度驗證

在基于模擬Sentinel-2A衛星遙感數據獲得的D-fG和D-HI估算最高精度波段優選結果基礎上，利用Sentinel-2A真實衛星遙感數據進行冬小麥收獲指數估算。在利用寬波段Sentinel-2A衛星遙感數據估算冬小麥D-fG過程中，首先根據地面采樣點的GPS定位信息對預處理后的遙感影像進行反射率提取，構建相應的NDSI與實測D-fG間的模型；然后，根據地面高光譜敏感波段所對應的Sentinel-2A波段，直接運用波段計算獲取相應NDSI的空間分布，在此基礎上實現區域冬小麥D-fG的空間提取并進行精度驗證。基于Sentinel-2A的紅光波段和窄近紅外波段反射率構建的NDSI與D-fG間擬合模型及D-fG估算精度如圖9、10所示。其中，紅光波段和窄近紅外波段構建NDSI進行冬小麥D-fG空間信息估算的R2為0.907 1，RMSE、NRMSE、MRE分別為0.044 3、13.13%、11.90%。本研究基于上述最高精度的Sentinel-2A影像寬波段構建的NDSI實現了衡水地區冬小麥D-fG空間信息分布獲取，其中，灌漿前期、灌漿后期、成熟期D-fG空間估算結果如圖11所示。通過分析可知，灌漿前期D-fG估算值在0.19～0.51之間，研究區D-fG平均值為0.23；灌漿后期D-fG估算值在0.27～0.62之間，研究區D-fG平均值為0.40；成熟期 D-fG估算值在0.44～0.68之間，研究區D-fG平均值為0.53。

圖9 基于Sentinel-2A衛星紅光和窄近紅外波段構建 NDSI的D-fG估算模型Fig.9 D-fG estimation model based on NDSI constructed by red band and NIR2 band of Sentinel-2A

圖10 基于Sentinel-2A紅光和窄近紅外波段衛星 遙感數據的冬小麥D-fG精度驗證結果Fig.10 Verification of D-fG of winter wheat based on red band and NIR2 band of Sentinel-2A

圖11 基于Sentinel-2A衛星遙感紅光和窄近紅外波段的冬小麥D-fG空間估算結果Fig.11 Spatial estimation results of winter wheat D-fG based on red band and NIR2 band of Sentinel-2A

3.4.2冬小麥收獲指數空間信息估算和總體精度驗證

在冬小麥D-fG空間信息獲取的基礎上，利用Sentinel-2A的紅光波段和窄近紅外波段構建的NDSI實現了衡水地區冬小麥D-HI空間信息估算，并進行總體精度驗證。不同時期D-HI空間信息估算結果如圖12所示，D-HI估算總體精度驗證結果如圖13所示。由圖12可知，冬小麥灌漿前期D-HI估算值在0.25～0.49之間，研究區D-HI平均值為0.28；冬小麥灌漿后期D-HI估算值在0.31～0.58之間，研究區D-HI平均值為0.41；冬小麥成熟期D-HI估算值在0.44～0.62之間，研究區D-HI平均值為0.51。通過驗證，基于Sentinel-2A衛星遙感紅光波段和窄近紅外波段估算D-fG參數的D-HI估算結果精度驗證R2為0.879 8，RMSE為0.050 2，NRMSE為13.81%，MRE為12.00%。上述基于紅光波段和處于近紅外區域的窄近紅外波段D-HI估算結果對寬波段多光譜衛星遙感數據選擇和傳感器波段設置具有重要指導意義。

圖12 基于Sentinel-2A衛星遙感紅光和窄近紅外波段的冬小麥D-HI空間信息估算結果Fig.12 Spatial estimation results of D-HI based on red band and NIR2 band of Sentinel-2A

圖13 基于Sentinel-2A衛星遙感紅光和窄近紅外 波段的冬小麥D-HI估算精度驗證結果Fig.13 Accuracy verification of D-HI estimation model based on red band and NIR2 band of Sentinel-2A

4 討論

與以往利用遙感信息獲取收獲指數精度指標RMSE一般在0.01～0.06之間相比[15,17,22,24]，本研究基于Sentinel-2A模擬遙感數據的作物收獲指數估算取得了較高的估算精度，且估算D-HI最高精度R2為0.908 3，RMSE為0.040 4，NRMSE為10.83%，而基于真實衛星遙感數據和相關最優波段的D-HI估算中，總體驗證精度R2為0.879 8，RMSE為0.050 2，NRMSE為13.81%。上述結果一方面說明本研究所提方法在田間冠層尺度擴展到區域尺度范圍內理論上是可行的，另一方面也說明在基于衛星遙感的區域收獲指數估算研究過程中，由于受大氣影響造成衛星遙感反射率數據出現一定誤差或波動，從而在一定程度上降低了作物收獲指數估算模型的應用精度。雖然基于真實衛星遙感數據估算D-HI的精度略有下降，但仍然能夠滿足區域農業遙感監測中收獲指數獲取的精度要求，對于大范圍區域農作物收獲指數空間信息獲取仍具有一定指導意義。同時，與以往僅實現成熟期收獲指數估算比較[15,17,22,24]，本研究所提基于花后累積地上生物量比例的區域冬小麥收獲指數遙感估算方法不僅可以實現成熟期收獲指數的準確估算，而且可以實現不同灌漿階段作物收獲指數動態變化空間信息獲取。從應用前景看，本研究所提收獲指數獲取技術方法不僅對提高精準農業農田作物管理能力具有重要意義，對開展作物生長過程定量模擬、作物品種選育和表型信息獲取、作物生長環境評價、作物災害損失定量評價、農業對氣候變化的響應等研究和應用也具有重要實用價值[43]。此外，本研究目前只針對冬小麥開展了收獲指數研究，但對于其他作物，特別是禾本科糧食作物(如水稻、玉米、高粱等)收獲指數遙感估算也具有重要的參考價值。

本研究在利用前期作物冠層高光譜確定的冬小麥D-fG敏感波段中心以及最大波段寬度結果基礎上，基于近地高光譜模擬Sentinel-2A數據及其真實衛星遙感數據開展冬小麥收獲指數遙感估算研究取得了較好的結果，但研究過程中忽略了一些外部因素的影響，且部分問題還需要進一步深入研究。如篩選出的敏感波段最大波寬所對應的波段范圍與遙感數據的波段不能完全對應，這在一定程度上影響了最大波段寬度的實際效果；其次，多光譜遙感數據構建的光譜指數與敏感波段構建的光譜指數在數量以及位置上難以一一對應，這使得部分光譜指數組合形式在作物收獲指數估算過程中無法進一步應用；再者，研究中Sentinel-2A遙感數據過境時間與地面實際采樣時間存在一定時間差，考慮到較短時間差內作物生長變化較小，因此，本研究尚未考慮遙感數據過境時間與地面采樣時間不完全一致對收獲指數估算精度的影響，但今后需進一步開展衛星遙感和地面觀測同步的作物收獲指數估算精度驗證研究。此外，本研究在冠層尺度冬小麥收獲指數敏感波段中心及其最大波寬優選結果基礎上，系統開展了Sentinel-2A模擬數據不同波段組合下的D-HI估測，對今后基于真實遙感數據的D-HI估測遙感數據選擇、波段組合效果評價、衛星傳感器波段應用與設置(如紅邊波段)等具有一定指導意義，但本文僅利用Sentinel-2A真實遙感波段數據最優波段組合獲取了高精度的區域收獲指數空間信息，然而基于真實衛星數據其他波段組合的區域作物收獲指數空間信息估測精度還有待進一步研究和評價。同時，除了Sentinel-2以外，在考慮作物物候信息基礎上，如何進一步利用具有紅邊的其他衛星(如RapidEye、國產GF-6等)進行更大范圍的作物收獲指數估算也將是今后的重要研究內容。

5 結論

(1)在作物收獲指數估算敏感波段中心和最大波段寬度基礎上，基于Sentinel-2A模擬遙感反射率數據，通過構建不同波段組合NDSI與實測D-fG間統計模型，實現了不同波段組合條件下冬小麥D-fG參數的遙感準確估算。在此基礎上，基于D-fG與D-HI間統計關系模型和D-fG參數遙感信息，實現了波段組合條件下灌漿期不同階段和成熟期冬小麥D-HI的準確估測和最優波段組合的篩選。其中，紅光波段、紅邊波段與近紅外波段間組合進行D-HI估測精度均較高，且基于模擬反射率波段λ1(672～680 nm)和λ2(855～875 nm)估算D-HI的精度最高，RMSE、NRMSE和MRE分別為0.040 4、10.83%、9.56%。上述模擬遙感數據不同波段組合的D-HI估測結果說明，Sentinel-2A衛星遙感數據在基于花后累積地上生物量比例的冬小麥收獲指數估算中具有一定的應用潛力。

(2)在基于模擬Sentinel-2A多光譜遙感數據進行作物動態收獲指數估算篩選出的最高精度波段組合信息基礎上，利用真實Sentinel-2A衛星遙感數據開展了花后累積地上生物量比例的作物D-HI遙感估測應用，實現了區域冬小麥D-fG參數遙感獲取和區域D-HI準確估測。其中，利用Sentinel-2A的紅光波段(R)和窄近紅外波段(NIR2)進行區域D-HI估測的RMSE、NRMSE和MRE分別為0.050 2、13.81%、12.00%。上述區域冬小麥動態收獲指數估算精度能夠滿足農業遙感監測中對收獲指數動態信息估測精度的要求，證明了基于Sentinel-2A衛星遙感數據和花后累積地上生物量比例的作物動態收獲指數估算方法在大范圍作物D-HI空間信息獲取中的可行性和有效性。