柴國奇， 王靜璞， 王光鎮(zhèn)， 韓柳， 王周龍

(魯東大學資源與環(huán)境工程學院，煙臺 264025)

0 引言

草地作為全球分布最廣的植被類型之一，在陸地生態(tài)系統(tǒng)碳循環(huán)中起著重要作用[1]。從功能角度看，草地植被可以分為光合植被(photosynthetic vegetation，PV)和非光合植被(non-photosynthetic vegetation，NPV)2類，其中NPV主要指枯草、枯葉、枯枝、枯干和凋落物等[2]。NPV作為草原生態(tài)系統(tǒng)的重要組成部分，是衡量地表植被覆蓋狀況的重要參數[3-4]。NPV不僅能夠減緩地表徑流和蒸發(fā)，保持水土和營養(yǎng)物質，增加土壤水分滲透和有機質含量，提高土壤濕度和溫度，改善土壤質量[5-7]，而且能影響生態(tài)系統(tǒng)的碳、水和能量的流動與循環(huán)[8]。另外，NPV在火災風險和頻率、風和水侵蝕等方面也起著關鍵作用[9]。因此，及時準確地掌握草原非光合植被覆蓋度(fractional cover of non-photosynthetic vegetation，fNPV)信息對草地資源的合理利用以及草原生態(tài)系統(tǒng)的監(jiān)測和保護具有重要的應用價值。

遙感技術的發(fā)展為大范圍快速、準確地獲取fNPV信息以及進行連續(xù)觀測提供了有效的技術手段。近幾十a來，國內外學者針對fNPV的遙感估算方法已開展了不少研究，并取得了一定的研究進展[2-5,10-15]。Daughtry等[5]針對NPV在2 100 nm附近波段的吸收特征構建的高光譜纖維素吸收指數(cellulose absorption index，CAI)表現較好，其與fNPV顯著相關(決定系數R2=0.89)。但由于高光譜數據難以獲取，無法滿足大規(guī)模估算fNPV的需求[16-17]，如EO-1 Hyperion傳感器已超期運行且幅寬僅有7.5 km，面臨著大范圍長期監(jiān)測的適用性難題[18-19]。因此，研究者們開始挖掘多光譜數據對fNPV的估算潛力，如利用TM數據構建歸一化差異指數(normalized difference index，NDI)[10]、歸一化差異耕作指數(normalized difference tillage index，NDTI)[20]以及歸一化差異衰老植被指數(normalized difference senescent vegetation index，NDSVI)[21]等來反演fNPV，然而這些指數容易受到土壤背景的影響，且不適用于草原地區(qū)PV，NPV和裸土(bare soil，BS)3種組分都存在情況下的fNPV估算[3,11]。因此Guerschman等[2]提出了基于MODIS數據B7和B6波段的比值指數(shortwave infrared ratio，SWIR32)，結合歸一化植被指數(normalized difference vegetation index，NDVI)構建三元線性混合模型比較準確地估算了澳大利亞草原地區(qū)的fNPV； Cao等[11]提出基于MODIS波段的干枯燃料指數(dead fuel index，DFI)估算PV，NPV和BS這3種組分混合情況下的fNPV，結果表明DFI指數與fNPV相關性較好(R2=0.95)。但Cao等[11]的研究中混合場景覆蓋度是利用端元平均光譜通過光譜混合模型模擬得到的，并不是真實的混合場景覆蓋度，因此DFI指數表征野外真實混合場景fNPV的有效性有待進一步驗證。

MODIS數據具備大范圍、高時效、低成本等優(yōu)勢，已成為草地資源監(jiān)測的重要數據源[1,22]。由于大多數基于MODIS數據構建的非光合植被指數(non-photosynthetic vegetation indices，NPVIs)相對較新，目前將這些指數應用于草原地區(qū)開展fNPV大范圍監(jiān)測的研究還較少。為此，本文選取錫林郭勒典型草原為研究區(qū)，以MCD43A4產品為數據源，比較基于MODIS數據構建的幾種NPVIs指數與fNPV之間的相關關系，建立NPVIs指數反演fNPV模型并進行精度驗證，以期能夠找到適用于典型草原地區(qū)大范圍應用的fNPV監(jiān)測方法。

1 研究區(qū)概況及數據源

1.1 研究區(qū)概況

研究區(qū)位于內蒙古自治區(qū)錫林郭勒盟(E113°30′～120°00′，N43°15′～46°41′)，包含錫林浩特市、東烏珠穆沁旗、西烏珠穆沁旗和阿巴嘎旗(圖1)，海拔在800～1 400 m之間。

圖1 研究區(qū)及實測樣地位置Fig.1 Study area and sampling location

錫林郭勒草原屬于中溫帶半干旱大陸性季風氣候區(qū)，四季分明，氣溫偏低，年平均氣溫0～3 ℃，無霜期90～120 d； 雨季短促，年平均降水量為295 mm，由東南向西北遞減，主要集中在夏季和秋季。研究區(qū)地貌以波狀高平原為主，由東南向西北方向傾斜，東、南部為大興安嶺和陰山山脈延伸余脈的交錯地段； 西、北部地形平坦，為高原草場； 土壤以風沙土為主，部分地區(qū)有栗鈣土、粽鈣土； 地帶性植被由草甸草原向典型草原、荒漠草原過渡，植被群落中主要有羊草、大針茅、克氏針茅和蒿類(冷蒿和變蒿)等。

1.2 野外數據獲取

2017年9月28日—10月3日在研究區(qū)內開展了野外地面觀測實驗。此時正處于植被黃枯期，采集的數據為實測樣地PV，NPV和BS覆蓋度數據。本次實驗共獲得49個樣地數據(圖1)，具體設置方法為： 首先在研究區(qū)相對平坦均質的地區(qū)選擇一個500 m×500 m大小的樣地，對樣地點編號并在樣地中心通過手持全球定位系統(tǒng)(global positioning system,GPS)記錄相應的地理位置信息以及樣地植被情況和地貌特征； 其后在每個樣地內布設3～5個0.5 m×0.5 m的小樣方，取所有小樣方植被覆蓋度的平均值作為該樣地最終的植被覆蓋度。與此同時，在樣地中心布設一個20 m×20 m的大樣方，用于與較高空間分辨率影像對應。每個0.5 m×0.5 m小樣方植被覆蓋度測算采用照相法[23]，利用數碼相機于樣方中心點正上方1.5 m處垂直拍攝2張照片，記錄每個小樣方對應的照片編號以及樣方內植被信息。

在測算小樣方植被覆蓋度前，首先通過ENVI5.3軟件切除照片中非樣方部分。照片剩余樣方部分用于進行植被覆蓋度的測算，根據照片的RGB值采集PV，NPV和BS這3類地物的訓練樣本(陰影劃分為BS類別)； 然后利用監(jiān)督分類中的最大似然分類法(maximum likelihood classification，MLC)對RGB照片進行分類，并結合目視解譯方法判定分類結果的精度，最后根據NPV的像元數量統(tǒng)計出各個樣方的fNPV，取2次統(tǒng)計結果的平均值作為該樣方最終的植被覆蓋度[13]。

1.3 遙感影像獲取及其預處理

1.3.1 MODIS數據

MODIS數據采用經過雙向反射分布函數(bidirectional reflectance distribution function，BRDF)校正處理為天頂方向觀測的地表反射率(nadir BRDF adjusted reflectance，NBAR)產品數據MCD43A4(V006)，NBAR使用16 d內的數據來計算BRDF效應并對反射率數據進行歸一化處理。MCD43A4(V006)提供了包括可見光、近紅外和短波紅外的B1—B7波段每日柵格化L3級數據產品，投影為正弦曲線投影，空間分辨率為500 m。通過LP DAAC網站(https: //lpdaac.usgs.gov)獲取到研究區(qū)2017年9月30日無云覆蓋的MCD43A4產品數據。

MCD43A4產品已經過幾何精糾正、輻射定標和大氣校正等處理，因此僅需對數據進行投影轉換、圖像拼接與裁剪處理即可。具體處理步驟如下： 首先，利用MRT (MODIS reprojection tools)工具將下載的數據進行投影轉換，將地圖投影轉換為WGS 84橢球的Albers Conical Equal Area投影，同時完成影像的拼接和重采樣； 然后，通過ENVI 5.3軟件利用研究區(qū)邊界矢量數據裁剪出研究區(qū)范圍； 最后，采用閾值法對影像中水體像元進行掩模處理[24]。

1.3.2 Landsat8 OLI數據

Landsat8 OLI數據從美國地質調查局網站下載(https: //glovis.usgs.gov)，獲取時間為2017年10月5日，研究區(qū)內基本無云。由于Landsat8 OLI數據已經過幾何糾正，因此只需對其進行輻射定標和大氣校正即可。影像預處理主要利用ENVI 5.3軟件，采用OLI傳感器的輻射定標參數對原始影像進行輻射定標，然后通過FLAASH模型對影像進行大氣校正，最終得到地表真實反射率。此外，為了去除水體對指數計算的影響，采用水體指數法對影像上的水體進行掩模處理。

2 研究方法

2.1 NPVIs

2.1.1 基于等效于TM的MODIS波段構建的指數

目前，常用的NPVIs大多數基于TM影像構建且多為NDI[25]，如NDI5，NDI7，NDTI和NDSVI等指數。而在可見光—近紅外—短波紅外波段(400～2 500 nm)范圍內，MODIS的B1，B2，B3，B4，B6，B7波段分別與Landsat5 TM的B3，B4，B1，B2，B5，B7波段對應。已有研究表明，基于不同TM波段構建的NDI5，NDI7，NDTI，NDSVI及土壤耕作指數(soil tillage index，STI)同樣適用于等效的MODIS波段[19]，計算公式分別為

NDI5=(B2-B6)/(B2+B6)，

(1)

NDI7=(B2-B7)/(B2+B7)，

(2)

NDTI=(B6-B7)/(B6+B7)，

(3)

NDSVI=(B6-B1)/(B6+B1)，

(4)

STI=B6/B7，

(5)

式中B1，B2，B6和B7分別為MODIS數據B1，B2，B6和B7波段的反射率。

2.1.2 基于MODIS波段構建的指數

DFI與SWIR32指數[2,11]計算公式分別為

(6)

(7)

其中，為了擴大DFI差異性，將DFI值擴大了100倍。

利用ENVI5.3軟件中的band math工具根據式(1)—(7)對各個NPVIs指數進行波段組合計算。在ArcGIS10.2軟件中利用地面實測樣地經緯度信息提取各樣地相對應的NPVIs指數值。

2.2 模型建立與驗證

2.2.1 相關分析

將提取的NPVIs指數的計算值與fNPV實測值進行相關性分析，檢驗二者之間關系的相關程度，通過差異顯著性選擇相關性較好的NPVIs指數用來建立遙感估算模型。本研究利用SPSS統(tǒng)計軟件進行相關性分析，再根據結果分析篩選最佳估算模型。

2.2.2 模型建立與驗證

對49個樣地的fNPV實測值從大到小排列后，按照隔二取一的模式選取建模樣本(33個)與驗證樣本(16個)； 然后采用線性回歸分析方法，建立以NPVIs指數為自變量和fNPV實測值為因變量的遙感估算模型。模型驗證選用R2、均方根誤差(root mean square error，RMSE)和相對誤差(relative error，RE)作為評價指標，R2越大、RMSE和RE越小，說明fNPV估算模型驗證結果越好。其計算公式分別為

(8)

(9)

2.2.3 最優(yōu)模型的穩(wěn)定性驗證

基于MODIS數據的fNPV估算模型穩(wěn)定性驗證的思路是利用線性回歸方法建立較高空間分辨率的影像數據(Landsat8 OLI)和地面實測fNPV (20 m×20 m大樣方)的關系，進而得到較高空間分辨率的fNPV分布結果； 再以此為基礎，通過ENVI5.3軟件將較高空間分辨率的數據按照像元簡單平均的重采樣方式尺度上推到MODIS數據500 m的空間分辨率[26-27]； 最后隨機布設驗證點，進行對比分析。另外，雖然Landsat8 OLI空間分辨率30 m略大于20 m大樣方，但由于樣方設置在相對均勻的區(qū)域，所以這種差別基本不會在fNPV估算中造成影響； 本文選用Landsat8 OLI數據估算fNPV適用性較好的NDTIOLI指數[14]，計算公式為

(10)

式中B6OLI和B7OLI分別為Landsat8 OLI數據B6和B7波段的反射率。

3 結果與分析

3.1 NPVIs指數與fNPV的相關性

本研究對選擇的7種NPVIs指數DFI，SWIR32，STI，NDTI，NDI7，NDI5，NDSVI及研究區(qū)內49個樣地的fNPV實測值進行相關性分析，結果如表1所示。NPVIs指數與研究區(qū)fNPV實測值之間的相關分析中除了NDI5和NDSVI相關性相對較差之外，其余5個指數與fNPV實測值均存在著較好的相關性，相關性系數R絕對值均大于0.64，達到極顯著相關水平(P<0.01)。其中DFI的相關性最好，R可達0.77； SWIR32次之，R為-0.72； NDI7相關性較低，R僅為0.64。DFI，SWIR32和NDTI與fNPV實測值的高相關性表明，利用NPVIs指數構建典型草原地區(qū)fNPV估算模型是基本可行的。

表1 NPVIs指數與fNPV實測值的相關性Tab.1 Correlation between NPVIs and fNPV

3.2 NPVIs指數構建fNPV反演模型

3.2.1 模型構建

采用線性回歸分析方法，將遙感影像提取的NPVIs指數(DFI，STI，NDTI，NDI7和SWIR32)與研究區(qū)33個建模樣本的fNPV實測值進行模型構建(表2)。由于NDI5和NDSVI與fNPV實測值的相關性相對較差，故在建立模型時去除。

表2 NPVIs指數與fNPV之間的反演模型Tab.2 Inversion models of non-photosyntheticvegetation index and fNPV

從表2可知，5種NPVIs指數與fNPV實測值所建立的線性回歸模型均能通過極顯著性相關水平(P<0.01)的F檢驗，R2在0.42～0.57之間。5種NPVIs中，DFI反演模型的擬合效果最好，R2為0.57； 其次是SWIR32指數； STI，NDTI以及NDI7指數R2均小于0.50，且都與DFI差距較大，其中NDI7的擬合效果較不理想，R2僅為0.42。通過各NPVIs指數與fNPV實測值間的相關系數及反演模型的擬合R2可以看出，相比于其他幾種指數，DFI指數反演模型的擬合效果更好，能夠有效模擬fNPV的變化規(guī)律。

3.2.2 模型精度評價

利用野外同期實測樣地的另外16個驗證樣本作為實測數據與各NPVIs指數反演模型得到的fNPV估算結果進行比較，評價不同NPVIs反演模型的估算精度，檢驗不同模型在典型草原地區(qū)的適用性(表3)。

從表3中可以看出，fNPV估算模型中利用DFI指數反演模型得到的估測值和實測值相關性最好，R2為0.63，RMSE與RE在5種NPVIs指數中也是最小的，RE為26.73%； SWIR32的估算效果稍次于DFI，RE為29.83%； NDTI介于SWIR32與STI之間，其RMSE與RE也均介于SWIR32與STI之間； NDI7指數的效果最差，其誤差也較大，RE為34.69%。

各NPVIs指數反演模型的預測值與實測值散點圖如圖2所示。從圖2中可以看出SWIR32，NDTI，STI及NDI7指數的反演模型高估了fNPV的低值部分，預測值位于Y=X線的上方較多，尤其是NDI7估算模型高估現象最為嚴重。從散點圖相關方程式可知，DFI指數斜率更接近于1，且R2較大，其估算結果比較接近實測值。通過比較分析驗證模型的R2和斜率，總體上看DFI指數反演模型(R2=0.63，RMSE=0.098 7，RE=26.73%)的估算效果優(yōu)于其他4種NPVIs指數，fNPV預測值和實測值比較接近于Y=X線附近，這說明該模型具有較好的適用性，在典型草原地區(qū)fNPV估算方面更具優(yōu)勢，能夠滿足估算的需要。

(a) DFI指數 (b) SWIR32指數 (c) NDTI指數

(d) STI指數 (e) NDI7指數

圖2 不同fNPV估算模型預測值與實測值散點圖

Fig.2ScatterplotsofdifferentfNPVestimatedmodelsbetweenpredictedvalueandmeasuredvalue

3.3 最優(yōu)NPVIs指數反演fNPV模型的穩(wěn)定性驗證

通過獲得的2017年10月5日Landsat8 OLI影像提取出NDTIOLI指數并與fNPV實測值(20 m×20 m大樣方)進行回歸分析，其擬合R2為0.58(P<0.01)，RMSE為0.092 9，表明該模型擬合程度較高，估算誤差較小，可保證Landsat8 OLI估算fNPV的準確性。

為進一步驗證最優(yōu)反演模型的穩(wěn)定性，采用NDTIOLI指數回歸模型計算得到Landsat8 OLI(30 m空間分辨率)的fNPV并將尺度上推到500 m空間分辨率； 隨機布設50個驗證點，將尺度一致的OLI-fNPV與DFI指數反演fNPV模型的估算結果MODIS-fNPV進行對比分析(圖3)。

圖3 研究區(qū)OLI-fNPV與MODIS-fNPV對比Fig.3 Comparison between OLI-fNPV and MODIS-fNPV

從散點圖可以看出散點均勻分布在Y=X線附近，誤差比較小，說明OLI-fNPV與MODIS-fNPV具有較好的一致性，進一步表明DFI指數反演fNPV模型的穩(wěn)定性較好。

3.4 fNPV的空間分布特征

通過上述分析，DFI指數反演fNPV模型在研究區(qū)估算精度最高且穩(wěn)定性較好。因此，利用其與fNPV間的線性回歸模型進行反演，即

fNPV=0.048 7DFI-0.188 3，

(11)

式中fNPV為fNPV預測值。

該模型R2為0.57。研究區(qū)2017年9月30日的fNPV空間分布如圖4所示。

圖4 研究區(qū)fNPV分布Fig.4 Distribution of fNPV in study area

從圖4中可以看出9月末fNPV總體的分布特征與本次野外調查基本一致，從西南到東北fNPV值逐漸增大，整體呈現東北高西南低的趨勢。分級統(tǒng)計典型草原地區(qū)的fNPV(表4)，主要介于(0.4，0.6]之間，像元個數大約占研究區(qū)的44.61%，這些像元主要分布在降水較多的東部； 而在降水相對較少的西部，大部分地區(qū)小于0.2。

表4 2017年9月30日各個fNPV等級所占百分比Tab.4 Area percent of each fNPV grade on Sep.30,2017

4 討論

4.1 不同NPVIs指數與fNPV間相關性差異

本研究探討了幾種由MODIS多光譜數據提取的NPVIs指數對典型草原地區(qū)fNPV的監(jiān)測效果。研究表明，常用的遙感估算模型中NPVIs指數與fNPV實測值的相關性均達到顯著相關水平，相關性由大到小為： DFI>SWIR32>NDTI>STI>NDI7>NDI5>NDSVI。NDI指數雖然與fNPV具有線性關系，但受BS背景的影響較大[10]； NDTI指數能夠增強NPV與BS的差異，具有較好的穩(wěn)定性[14,28]，但容易受到PV的影響，從而掩蓋了NPV的信息[29]； SWIR32與STI同樣受PV的干擾較大[19]； DFI指數的表現優(yōu)于NDI，NDTI，STI，SWIR32以及NDSVI等指數，原因在于DFI相較于其他幾種指數利用了更多的波段(B1，B2，B6和B7)參與運算，一定程度上減少了BS背景效應的影響，提高了估算精度； 并且DFI指數是基于PV，NPV和BS這3種組分都存在的條件下建立的，能夠較準確估算PV，NPV以及BS這3種組分混合情況下草原地區(qū)的fNPV[11]。

4.2 最優(yōu)反演fNPV模型選擇

通過線性回歸分析方法對錫林郭勒典型草原地區(qū)建立了一系列NPVIs指數反演fNPV模型，其中DFI指數反演模型精度較高(R2=0.63,RMSE=0.098 7)，RE僅為26.73%。這與Cao等[30]基于MODIS09GA建立的草原地區(qū)反演模型研究結果相似，但仍存在一些差異，需要進一步研究。NDI，NDTI，STI以及SWIR32指數用于本研究區(qū)估算結果并不理想，主要原因在于這幾種指數僅基于除PV以外的2種組分(NPV和BS)創(chuàng)建，或受BS背景影響大或受PV的干擾嚴重，限制了其在NPV，PV和BS這3種組分混合情況下的應用[31]。因此，選擇DFI指數用于草原地區(qū)最優(yōu)fNPV估算模型的構建。

DFI指數由Cao等[11]基于MODIS數據波段最先提出，用以估算PV，NPV和BS這3種組分混合情況下的fNPV，但是目前國內外利用該指數直接與MODIS影像構建經驗模型估算fNPV還鮮有研究。本次研究證明了基于MODIS數據DFI指數能夠有效估算典型草原地區(qū)的fNPV，并且精度較高。因此，借助于MODIS數據大范圍監(jiān)測、高時間分辨率的優(yōu)勢有必要進一步探究DFI用于fNPV時空動態(tài)變化監(jiān)測的能力。

4.3 fNPV估算的復雜性和誤差來源

利用遙感估算fNPV受到眾多外部因素的影響，如水分吸收特征幾乎能夠完全掩蓋NPV的光譜特征[19]。此外，不同BS類型、地表濕度和有機質含量等因素均影響NPV估算結果[25,32]。因此，遙感估算fNPV具有復雜性，需要對各種影響因素進一步分析，以更好地證明fNPV估算模型的穩(wěn)定性和魯棒性。

值得注意的是，由于野外實測數據與影像數據存在時間上不同步、樣方植被覆蓋度計算具有主觀性的問題，這些數據誤差也在一定程度上會影響估算精度。另外，本研究僅對一次實驗結果進行了精度檢驗，樣地數量相對較少且缺乏系統(tǒng)的地面真實性檢驗工作，再加上線性回歸模型對實測數據的依賴性較高，模型的普適性并未得到充分的驗證。因此，在后續(xù)的研究中將選取固定觀測點對長時間序列的植被覆蓋度情況進行連續(xù)觀測，獲取更多的地面實測樣地數據，以更好地評價NPVIs指數估算典型草原地區(qū)fNPV的精度。

5 結論

本文以錫林郭勒典型草原為研究區(qū)，MODIS 500 m空間分辨率產品數據為數據源，對幾種fNPV估算模型進行了對比分析和精度驗證，得到如下主要結論：

1)基于MODIS影像構建的NPVIs指數與fNPV的相關性較好，達到顯著相關水平(P<0.05)，相關性依次為： DFI>SWIR32>NDTI>STI>NDI7>NDI5>NDSVI。

2)野外樣地fNPV實測值驗證表明，相比于其他NPVIs指數反演fNPV模型，DFI指數估算效果較好，fNPV估算的R2為0.63，RMSE為0.098 7，RE為26.73%，能夠有效估算典型草原地區(qū)的fNPV。

3)基于DFI指數構建的fNPV反演模型為最佳模型，fNPV總體的分布規(guī)律與實際相吻合，可有效應用于典型草原地區(qū)黃枯期fNPV的快速監(jiān)測。